随着我国北斗卫星导航系统(BeiDou navigation satellite system，BDS)的建设与完善，卫星数据的易用性和精度都得到了大幅改善^[1]。数据预处理是BDS高精度定位与导航的重要前提和保障，其中周跳的处理是必不可少的一环。BDS三频信号相较于双频信号可为周跳处理提供众多电离层延迟小、噪声低、波长长的探测组合^[2-3]，具有更广阔的应用前景。

针对三频周跳探测与修复的特性，文献[4-5]研究了无电离层伪距相位周跳探测与修复方法，但组合探测敏感性易受到伪距噪声的影响，探测效果较差；文献[6]对几个无几何相位组合进行最优组合，但该优选组合不能构成3个线性无关探测量，为此使用了最小二乘模糊度降相关平差(least-square ambiguity decorrelation adjustment，LAMBDA)法进行搜索处理，虽然取得了较好的效果，但过程较为复杂；文献[7]对比分析了无几何相位组合法和伪距相位组合法的算法复杂度、探测效果以及受电离层影响程度等，综合表明伪距相位组合法更具优势。而当采用混合组合法联合两个无几何相位组合以及一个伪距相位组合处理周跳时，往往会受到法方程组病态的影响，使得众多学者不得不避开这些性质优良的探测组合^[8-9]，针对这一问题，目前的相关研究较少，为了更好利用组合观测值性质进行周跳数据处理，有必要对此展开进一步研究。

本文首先构建了两个无几何相位组合和一个无几何伪距相位组合，分析了其对各类周跳的探测能力，然后针对该探测组合存在的方程组病态的问题，从有偏估计角度出发，提出了一种基于拓展岭估计的周跳实时修复方法，最后利用北斗三频实测数据进行了验证与分析，得出相应结论。

1 周跳探测量构造 1.1 无几何相位组合

根据文献[10-11]，忽略多路径影响，可将相位的观测方程简写为：

$ {L}_{i}={\lambda }_{i}{\phi }_{i}=\rho -\frac{I}{{f}_{i}^{2}}+{\lambda }_{i}{N}_{i}+{\lambda }_{i}\varepsilon _{{\phi }_{i}} $

(1)

式中，λ_i、$ {\phi }_{i}\left(t\right)$分别为第i个频点的波长和以周为单位的相位观测值；$ \rho $为与频率无关项；f_i为频率；I为电离层延迟值，$ I=40.3\cdot \mathrm{S}\mathrm{T}\mathrm{E}\mathrm{C}$，其中$ \mathrm{S}\mathrm{T}\mathrm{E}\mathrm{C} $(slant total electron content)表示斜向电离层电子含量；${N}_{i} $为整周模糊度；$\varepsilon _{{\varphi }_{i}} $表示观测误差。于是，可将三频相位组合表示为：

$ \begin{array}{c}{L}_{\alpha \beta \gamma }=\alpha {L}_{1}+\beta {L}_{2}+\gamma {L}_{3}=\alpha {\lambda }_{1}{\phi }_{1}+\beta {\lambda }_{2}{\phi }_{2}+\gamma {\lambda }_{3}{\phi }_{3}=\\ (\alpha +\beta +\gamma )\rho -{\eta }_{\alpha \beta \gamma }{I}_{1}+{N}_{\alpha \beta \gamma }+\varepsilon _{{L}_{\alpha \beta \gamma }}\end{array} $

(2)

式中，$ \alpha $、$ \beta $、$ \gamma $为组合系数；电离层延迟系数$ {\eta }_{\alpha \beta \gamma }=\alpha +\beta \cdot \frac{{f}_{1}^{2}}{{f}_{2}^{2}}+\gamma \cdot \frac{{f}_{1}^{2}}{{f}_{3}^{2}}$；$ {I}_{1}$为第一频点电离层延迟值；组合模糊度$ {N}_{\alpha \beta \gamma }=\alpha {\lambda }_{1}{N}_{1}+\beta {\lambda }_{2}{N}_{2}+\gamma {\lambda }_{3}{N}_{3}$；组合噪声$ \varepsilon _{{L}_{\alpha \beta \gamma }}=\alpha {\lambda }_{1}\varepsilon _{{\phi }_{1}}+\beta {\lambda }_{2}\varepsilon _{{\phi }_{2}}+\gamma {\lambda }_{3}\varepsilon _{{\phi }_{3}}$。为得到无几何组合观测值，令：

$ \alpha +\beta +\gamma =0 $

(3)

再对历元间作差可得：

$ \mathrm{\Delta }{L}_{\alpha \beta \gamma }=-{\eta }_{\alpha \beta \gamma }\mathrm{\Delta }{I}_{1}+\mathrm{\Delta }{N}_{\alpha \beta \gamma }+\mathrm{\Delta }\varepsilon _{{L}_{\alpha \beta \gamma }} $

(4)

于是，得到组合周跳实数估值：

$ \mathrm{\Delta }{\hat N}\alpha \beta \gamma =\mathrm{\Delta }{L}_{\alpha \beta \gamma }+{\eta }_{\alpha \beta \gamma }\mathrm{\Delta }{I}_{1}-\mathrm{\Delta }\varepsilon _{{L}_{\alpha \beta \gamma }} $

(5)

由式(5)可知，无几何相位组合周跳探测会受到历元间电离层残差项$ {\eta }_{\alpha \beta \gamma }\mathrm{\Delta }{I}_{1}$和历元间组合噪声$ \mathrm{\Delta }\varepsilon _{{L}_{\alpha \beta \gamma }} $的影响。相较于组合噪声，在电离层延迟放大系数$ {\eta }_{\alpha \beta \gamma } $接近于0的情况下，可忽略电离层残差对组合探测量灵敏度的影响，从而得到周跳估值的标准差：

$ {\sigma }_{\mathrm{\Delta }{\hat N}_{\alpha \beta \gamma }}=\sqrt{2\left({\left(\alpha {\lambda }_{1}\right)}^{2}+{\left(\beta {\lambda }_{2}\right)}^{2}+{\left(\gamma {\lambda }_{3}\right)}^{2}\right){\sigma }_\varepsilon } $

(6)

式中，$ {\sigma }_\varepsilon $为载波相位观测噪声。由式(6)可知，要想提高周跳探测灵敏度，还需满足：

$ {\left(\alpha {\lambda }_{1}\right)}^{2}+{\left(\beta {\lambda }_{2}\right)}^{2}+{\left(\gamma {\lambda }_{3}\right)}^{2}=\mathrm{m}\mathrm{i}\mathrm{n} $

(7)

判断周跳发生的阈值为：

$ \left|\mathrm{\Delta }{L}_{\alpha \beta \gamma }\right|>m{\sigma }_{{}_{\mathrm{\Delta }{\hat N}_{\alpha \beta \gamma }}} $

(8)

式中，m取3、4分别对应99.7%、99.9%的置信水平探测周跳。本文$ {\sigma }_\varepsilon $取0.01周, 探测阈值取$ 4{\sigma }_{\mathrm{\Delta }{\hat N}_{\alpha \beta \gamma }}$。为得到较小的周跳估值标准差，依据式(3)、式(7)及$ {\eta }_{\alpha \beta \gamma }$趋近于0的标准，表 1列出了一些性质较好的系数组合。

1.2 无几何伪距相位组合

在历元$ t $时刻，$ i$($ i$=1，2，3)频点的伪距与相位原始观测方程为：

$ {P}_{i}\left(t\right)=\rho \left(t\right)+\frac{I}{{f}_{i}^{2}}+\varepsilon _{{P}_{i}} $

(9)

$ {\phi }_{i}\left(t\right)=\frac{\rho \left(t\right)}{{\lambda }_{i}}-\frac{I}{{f}_{i}^{2}}\cdot \frac{1}{{\lambda }_{i}}+{N}_{i}+\varepsilon _{{\phi }_{i}} $

(10)

式中，$ {P}_{i}\left(t\right)$、$ {\phi }_{i}\left(t\right) $分别是$ t $时刻的伪距观测量(单位：m)和相位观测量(单位：周)；$ \varepsilon _{{P}_{i}}$、$ \varepsilon _{{\phi }_{i}} $分别为伪距及载波观测噪声。于是，三频伪距组合可表示为：

$ \begin{array}{c}{P}_{abc}=a{P}_{1}+b{P}_{2}+c{P}_{3}=(a+b+c)\rho +\\ {\eta }_{abc}{I}_{1}+\varepsilon _{{P}_{abc}}\end{array} $

(11)

式中，$ {I}_{1}=\frac{I}{{f}_{1}^{2}} $为载波$ {f}_{1} $路径上的电离层影响量；$ {\eta }_{abc}=a+b\frac{{f}_{1}^{2}}{{f}_{2}^{2}}+c\frac{{f}_{1}^{2}}{{f}_{3}^{2}} $为$ {I}_{1}$的系数；$ \varepsilon _{{P}_{abc}}=a\varepsilon _{{P}_{1}}+b\varepsilon _{{P}_{2}}+c\varepsilon _{{P}_{3}}$为伪距组合观测值的噪声。为消去几何项影响，令：

$ a+b+c=1 $

(12)

则三频相位组合可表示为：

$ {\phi }_{\alpha \beta \gamma }=\frac{1}{{\lambda }_{\alpha \beta \gamma }}\rho -{\eta }_{\alpha \beta \gamma }{I}_{1}+{N}_{\alpha \beta \gamma }+\varepsilon _{{\phi }_{\alpha \beta \gamma }} $

(13)

式中，组合波长为：

电离层延迟系数为：

由式(11)、式(13)可进一步组合，得到无几何伪距相位组合表达式：

$ {\phi }_{\alpha \beta \gamma }-\frac{{P}_{abc}}{{\lambda }_{\alpha \beta \gamma }}=-({\eta }_{\alpha \beta \gamma }+\frac{{\eta }_{abc}}{{\lambda }_{\alpha \beta \gamma }}){I}_{1}+{N}_{\alpha \beta \gamma }+\varepsilon _{{\phi }_{\alpha \beta \gamma }}-\frac{\varepsilon _{{P}_{abc}}}{{\lambda }_{\alpha \beta \gamma }} $

(14)

经分析，$ {\phi }_{\alpha \beta \gamma }-\frac{{P}_{abc}}{{\lambda }_{\alpha \beta \gamma }}$为弱电离层组合，对式(14)进行历元间作差可得：

$ \begin{array}{c}\mathrm{\Delta }{\phi }_{\alpha \beta \gamma }-\frac{\mathrm{\Delta }{P}_{abc}}{{\lambda }_{\alpha \beta \gamma }}=-({\eta }_{\alpha \beta \gamma }+\frac{{\eta }_{abc}}{{\lambda }_{\alpha \beta \gamma }})\mathrm{\Delta }{I}_{1}+\mathrm{\Delta }{N}_{\alpha \beta \gamma }+\\ \mathrm{\Delta }\varepsilon _{{\phi }_{\alpha \beta \gamma }}-\frac{\mathrm{\Delta }\varepsilon _{{P}_{abc}}}{{\lambda }_{\alpha \beta \gamma }}\end{array} $

(15)

于是，可得到周跳实数解估值：

$ \mathrm{\Delta }{\hat N}\alpha \beta \gamma =(\mathrm{\Delta }{\phi }_{\alpha \beta \gamma }-\frac{\mathrm{\Delta }{P}_{abc}}{{\lambda }_{\alpha \beta \gamma }})+({\eta }_{\alpha \beta \gamma }+\frac{{\eta }_{abc}}{{\lambda }_{\alpha \beta \gamma }})\mathrm{\Delta }{I}_{1} $

(16)

在经过历元间作差之后，电离层一阶项影响基本被消除，此时影响探测灵敏度的因素为二阶电离层延迟及观测噪声影响。式(16)的标准差为：

$ \begin{array}{c}{\sigma }_{\mathrm{\Delta }{\hat N}_{\alpha \beta \gamma }}=\sqrt{2\left(\left({\alpha }^{2}+{\beta }^{2}+{\gamma }^{2}\right){\sigma }_{\phi }^{2}+\left({a}^{2}+{b}^{2}+{c}^{2}\right)\frac{{\sigma }_{P}^{2}}{{{\lambda }_{\alpha \beta \gamma }}^{2}}\right)}+\\ ({\eta }_{i\alpha \beta \gamma }+\frac{{\eta }_{abc}}{{\lambda }_{i\alpha \beta \gamma }}){\sigma }_{\mathrm{\Delta }{I}_{1}}\end{array} $

(17)

判断周跳发生的阈值为：

$ \left|\mathrm{\Delta }{\hat N}\alpha \beta \gamma \right|=m{\sigma }_{\mathrm{\Delta }{\hat N}_{\alpha \beta \gamma }} $

(18)

式中，$ m$取3、4分别对应99.7%和99.9%的置信水平。本文$ m$取4，$ {\sigma }_{\phi }=0.01 $周，$ {\sigma }_{P}=0.3\mathrm{ }\mathrm{m} $。分析式(17)可知，由于$ {\sigma }_{P}$对$ {\sigma }_{\mathrm{\Delta }{\hat N}_{\alpha \beta \gamma }} $影响远大于$ {\sigma }_{\phi } $，因此，要想得到较小的$ {\sigma }_{\mathrm{\Delta }{\hat N}_{\alpha \beta \gamma }} $值，需使得组合波长$ {\lambda }_{\alpha \beta \gamma } $较长，并且满足：

$ {a}^{2}+{b}^{2}+{c}^{2}=\mathrm{m}\mathrm{i}\mathrm{n} $

(19)

结合式(12)可知，$ a=b=c=\frac{1}{3} $；组合后的电离层延迟放大系数为$ \theta =({\eta }_{\alpha \beta \gamma }+\frac{{\eta }_{abc}}{{\lambda }_{\alpha \beta \gamma }}) $；若要忽略电离层二阶电离层延迟项影响，则需电离层延迟放大系数$ \theta $趋向于0。满足上述3个条件才能得到较小的周跳估值标准差。表 2给出了一些性质较优的组合(表中$ {\sigma }_{\mathrm{\Delta }{\hat N}_{\alpha \beta \gamma }} $均已忽略二阶电离层延迟影响)。

由表 1、表 2可以看出，北斗三频电离层延迟放大系数与周跳估值标准差是一个较难调和的矛盾，往往$ \theta $较小的组合对应的周跳估值标准差却较大，这与文献[3]的研究保持一致。虽然无几何相位组合与无几何伪距相位组合中[1 0 -1]、[0 1 -1]等周跳估值标准差与$ \theta $值均较小，但是这类组合只能探测两个频点上的周跳，其不敏感周跳组合较多，因此与其他周跳探测组合联合探测时效果不佳。考虑到组合系数应满足$ \alpha +\beta +\gamma =0 $^[12]，所以即便选择3个性质较优的伪距相位组合仍无法探测出该历元3个频点发生相等周跳的情况，而无几何相位组合至多构成两个线性无关的观测量，若单独使用3个无几何相位组合探测周跳，则需要借助于LAMBDA搜索法^[6]，过程较复杂。综合考虑，选择联合两个无几何相位组合[1 1 -2]、[1 -2 1]以及一个无几何伪距相位组合[1 3 -4]组成3个线性无关的探测量。

2 基于拓展岭估计的周跳修复

将3个线性无关探测方程联立，可以得到：

$ \mathit{\boldsymbol{A}}\mathrm{\Delta }\mathit{\boldsymbol{N}}=\mathit{\boldsymbol{L}} $

(20)

式中，

$ \mathit{\boldsymbol{A}}=\left[\begin{array}{ccc}{\lambda }_{1}& {\lambda }_{2}& -2{\lambda }_{3}\\ {\lambda }_{1}& -2{\lambda }_{2}& {\lambda }_{3}\\ 1& 3& -4\end{array}\right] $

(21)

$ \mathrm{\Delta }\mathit{\boldsymbol{N}}={\left[\begin{array}{ccc}\mathrm{\Delta }{N}_{1}& \mathrm{\Delta }{N}_{2}& \mathrm{\Delta }{N}_{3}\end{array}\right]}^{\mathrm{T}} $

(22)

$ \mathit{\boldsymbol{L}}={\left[\begin{array}{ccc}{l}_{1}& {l}_{2}& {l}_{3}\end{array}\right]}^{\mathrm{T}} $

(23)

式中, $ {l}_{1} $、$ {l}_{2}$、$ {l}_{3}$分别是无相位组合[1 1 -2]、[1 -2 1]以及无几何伪距相位组合[1 3 -4]的探测量，条件数达到866.35, 可见$ \mathit{\boldsymbol{L}}$很小的扰动就将引起$ \mathrm{\Delta }\mathit{\boldsymbol{N}}$解发生很大的变化，病态性较严重；$ \mathrm{\Delta }{N}_{i} $$ \left(i=\mathrm{1, 2}, 3\right)$分别是3个频段上发生的周跳值。直接求解$ \mathrm{\Delta }{N}_{i}$时，首先将$ {l}_{3} $取整得$ \stackrel{-}{{l}_{3}} $，然后得到新的观测量$ L\prime ={\left[\begin{array}{ccc}{l}_{1}& {l}_{2}& \stackrel{-}{{l}_{3}}\end{array}\right]}^{\mathrm{T}}$, 接着解算法方程：

$ {\mathit{\boldsymbol{A}}}^{\mathrm{T}}\mathit{\boldsymbol{A}}\mathrm{\Delta }\mathit{\boldsymbol{N}}={\mathit{\boldsymbol{A}}}^{\mathrm{T}}L\prime $

(24)

可以得到：

$ \mathrm{\Delta }{\hat N}={\left({\mathit{\boldsymbol{A}}}^{\mathrm{T}}\mathit{\boldsymbol{A}}\right)}^{-1}{\mathit{\boldsymbol{A}}}^{\mathrm{T}}L\prime $

(25)

最后对式(25)解得的$ \mathrm{\Delta }{N}_{i}^{̂}$直接求整，即可得到整数周跳解。但由于此处$ \mathit{\boldsymbol{A}} $矩阵中各列向量之间存在近似线性关系，使得法方程系数阵$ {\mathit{\boldsymbol{A}}}^{\mathrm{T}}\mathit{\boldsymbol{A}} $条件数过大，导致法方程存在病态问题，其求解得到的周跳值极不稳定。为了解决这一问题，本文从有偏估计角度出发，提出一种基于拓展的岭估计的周跳修复方法。

岭估计^[13-14]实质上是一种改良的最小二乘算法，通过放弃最小二乘的无偏性换取均方误差更小、更可靠以及更符合实际的估计结果。传统的岭估计方法是在法方程系数阵对角线上加一个很小的正数，得到岭估计解：

$ \mathrm{\Delta }{\hat N}\prime ={\left({\mathit{\boldsymbol{A}}}^{\mathrm{T}}\mathit{\boldsymbol{A}}+k\mathit{\boldsymbol{I}}\right)}^{-1}{\mathit{\boldsymbol{A}}}^{\mathrm{T}}\mathit{\boldsymbol{L}} $

(26)

式中，$ k$值为一很小正数；$ \mathit{\boldsymbol{I}}$为3阶单位阵。但在本算法实际应用中，最优解对应$ k$值往往是一个很小的负数，故本文将岭参数$ k $值拓展到实数域，所进行的岭参数估计称为拓展的岭估计。

要想得到最优的岭估计量，关键在于$ k$值的选取。现阶段$ k$值选取方法有很多，包括岭迹法、方差扩大因子法以及特征根法等，但这些方法或不适用于编程实现，或不适用于实数域$ k$值解算。因此，本文借鉴残差平方和$ {\sigma }^{2}={\mathit{\boldsymbol{V}}}^{\mathrm{T}}\mathit{\boldsymbol{P}}\mathit{\boldsymbol{V}} $($ \mathit{\boldsymbol{V}}$为残差向量)，提出一种以$ {v}_{1} $、$ {v}_{2}$最小为原则、在区域$ R $中搜索$ k $值的方法。其中$ {v}_{1} $反映修复的整数周跳与观测量的相符程度, $ {v}_{2}$反映修复的实数周跳与取整得到的整数周跳之间的相符程度。计算如下：

$ {v}_{1}={\mathit{\boldsymbol{V}}}_{1}^{\mathrm{T}}{\mathit{\boldsymbol{P}}}_{1}{\mathit{\boldsymbol{V}}}_{1} $

(27)

$ {v}_{2}={\mathit{\boldsymbol{V}}}_{2}^{\mathrm{T}}{\mathit{\boldsymbol{P}}}_{2}{\mathit{\boldsymbol{V}}}_{2} $

(28)

式中, $ {\mathit{\boldsymbol{V}}}_{1}=\mathit{\boldsymbol{A}}\mathrm{\Delta }\stackrel{-}{N}-\mathit{\boldsymbol{L}}$，$ {\mathit{\boldsymbol{V}}}_{2}=\mathrm{\Delta }{\hat N}\prime -\mathrm{\Delta }\stackrel{-}{N}$, 其中$ \mathrm{\Delta }\stackrel{-}{N} $为修复的整数周跳，$ \mathrm{\Delta }{\hat N}\prime $为实数周跳值；为简化研究，将$ {\mathit{\boldsymbol{P}}}_{1} $、$ {\mathit{\boldsymbol{P}}}_{2} $设为单位阵。在正确修复周跳的情况下，$ {v}_{1}$值应取得最小，同时，由岭估计得到的最优实数周跳与取整得到的整数周跳之间残差的$ {v}_{2} $值应最小，否则不能保证该岭估计实数解取整的可靠性。

对于$ k$值搜索区域$ R $的选取，算法采用统计-验证-优化的方式。首先在BDS卫星三频观测数据中加入单历元模拟周跳，以$ {v}_{1}$和$ {v}_{2}$最小为原则确定各历元周跳解对应最优$ k $值。然后将统计的$ k$值集合$ {R}_{0}$适当扩大成区域$ {{R}_{0}}^{\prime } $，应用$ {{R}_{0}}^{\prime } $为下一观测数据$ k $值搜索范围，验证新的$ k$值是否属于$ {{R}_{0}}^{\prime } $。若不属于，则优化原有搜索区域$ {{R}_{0}}^{\prime } $，得到$ {R}_{0}″$，如此循环，直到$ {R}_{0}^{\left(n\right)}$达到稳定状态，取最终稳定的$ {R}_{0}^{\left(n\right)} $为$ k$值搜索区域$ R $。对2015—2018年间的20余个北斗观测测站，包含35个北斗卫星的三频观测数据，使用统计-验证-优化的方式确定区域$ R $之后，综合考虑搜索效率，得到本算法适用的$ k$值搜索区域$ R $和搜索步长$ {l}_{0} $。

设置周跳修复准确性检验指标$ {v}_{1} $和$ {v}_{2} $，当周跳修复完成后，若某历元$ {v}_{1} $或$ {v}_{2}$出现异常值，则认为该历元周跳修复失败，需要扩大岭参数$ k$的搜索区域$ R $，进一步搜索直到该异常值消失，该指标可有效确保周跳修复的可靠性。本文算法具体实现流程如图 1所示(算法中异常值历元指的是该历元$ {v}_{1} $或$ {v}_{2}$值相较于其他历元大得多或该历元自身明显异常，形成显著孤值)。

3 周跳实验与结果分析

实验采用实测BDS三频数据，在GNSSer软件(http://www.gnsser.com)上实现了本文提出的算法，实验前已采用文献[15]的方法对所选卫星历元进行了钟跳探测与修复。

3.1 模拟周跳实验

实验选取2017-02-14 CUT0站使用TRIMBLE NETR9接收机观测的北斗二号系统(BDS-2)C01卫星30 s采样间隔的静态观测数据，共2 800个历元，所选弧段卫星高度角均大于10°，该卫星数据未参与岭参数$ k$值搜索区域$ R $的确定。单历元模拟周跳实验在3个频段的第500、1 000、1 500、2 000历元处分别加入[4 0 1]、[1 2 1]、[0 1 2]、[1 10 15]组合周跳，其中[4 0 1]是伪距相位组合的不敏感周跳组合，而[0 1 2]是[1 -2 1]无几何相位组合探测量较小的周跳组合，[1 2 1]为一般小周跳，[1 10 15]为一般大周跳，单历元模拟周跳探测以及修复情况分别如表 3、图 2所示。

由图 2可知，本文所选的无几何相位组合的探测量在未发生周跳时基本稳定在0.02周以内，具有很小的观测噪声水平，伪距相位组合的探测量也稳定在0.2周左右，三者组合使用具有很高的探测灵敏度。同时，3个周跳探测组合虽然各自存在不敏感周跳组合，但并不相互重叠，理论上任意周跳在3个观测量的结合使用中均可被有效探测。由图 2(a)可以看出，[1 3 -4]伪距相位组合未能探测出自身的敏感周跳组合[4 0 1]，但该周跳组合均能够被其他两个无几何相位组合有效识别；同样图 2(c)中，[1 -2 1]无几何相位组合未能有效探测其自身不敏感周跳组合[0 1 2]，但该周跳组合被其余两个组合观测量成功探测，这与前文理论分析保持一致。图 2(b)、图 2(d)则分别说明了该探测组合对于一般小周跳、一般大周跳的探测能力。在周跳修复方面，由表 3可以看出各个历元添加的模拟周跳均被有效修复，算法设置的修复正确性判断指标均处于正常水平，其中，$ {v}_{1} $值反映出了修复的整数跳值与探测量之间具有非常良好的相符程度，同时，处于$ {10}^{-4}$量级的$ {v}_{2} $表明通过算法得到的实数周跳与对其取整得到的整数周跳之间的舍入误差非常小，取整过程具有极高的可靠性。

为了检验算法处理连续周跳的能力，在第500、501、502历元处分别添加[4 0 1]、[1 2 1]、[0 1 2]小周跳组合；在第1 000、1 001、1 002历元处分别添加[4 0 1]、[1 2 1]、[0 1 2]小周跳组合；在第1 500、1 501、1 502历元处分别添加[1 0 1]、[1 2 3]、[0 2 4]小周跳组合；在第2 000、2 001、2 002历元处分别添加[1 3 4]、[2 0 5]、[2 2 3]小周跳组合。模拟连续周跳探测以及修复情况分别如图 3、表 4所示。

图 3反映了本文算法探测连续周跳的能力。实验在第500、501、502历元处以及第1 000、1 001、1 002历元处分别加入了3个相同的连续周跳组合，由图 3(a)、图 3(b)及表 4可知，加入的周跳组合在3个组合观测量的结合使用中均被有效探测和修复，反映了算法对随机历元连续周跳的探测以及修复能力。图 3(c)、图 3(d)则反映了算法对随机连续小周跳的探测能力，结合表 4可知，所有添加的连续小周跳组合均被有效探测和修复，且相应修复正确性判断指标均处于正常范围。

3.2 随机周跳实验

为检验算法的周跳处理能力，在所选C01卫星3个频点(B1、B2、B3)每间隔5历元，对其之后历元随机逐个加入0~9周的小周跳，共576个周跳发生历元。由于算法使用了伪距观测值，而BDS-2和北斗三号系统(BDS-3)伪距观测值噪声水平不一致，为进一步测试算法的适用性，添加BDS-3卫星随机周跳实验。实验选取由国际全球连续监测评估系统(international GNSS monitoring & assessment system，IGMAS)分析中心提供的2018-10-01 ALGR站使用GNSS_ggr接收机观测的BDS-3 C32卫星静态数据，弧段共360个历元，卫星高度角均大于10°。在C32卫星的3个频点(B1C、B2a、B3I)每间隔5历元，对其之后历元随机逐个加入0~9周的小周跳，共72个历元发生周跳。对两组随机实验实时探测出来的周跳立即修复，处理结果如图 4~图 8所示。可以看出，两组实验均成功探测并修复了所有历元随机周跳。

图 4(a)、图 5(a)分别为C01卫星和C32卫星的周跳探测量，两组实验探测量的量级相近，探测组合中的无几何相位组合的探测量在未发生周跳时基本均在0.02周以内，伪距相位组合的探测量也稳定在0.2周左右; 图 4(b)、图 5(b)分别为两组实验添加的随机周跳值。图 6(a)、图 7(a)中, 两组$ {v}_{1}$值量级相近，均反映了修复的整数周跳与实际探测量之间的非常良好的相符程度；由图 6(b)、图 7(b)可知，对拓展岭估计方法得到的实数周跳解取整，造成的舍入误差较小，取整具有高度的可靠性，但比较来看，C01卫星较C32卫星的$ {v}_{2}$值量级小得多，一定程度表明了BDS-2卫星C01数据质量更优。综合$ {v}_{1} $与$ {v}_{2}$指标，可有效确保整个周跳修复的可靠性。由图 8可知，绝大部分$k $值均匀分布在0值周围，具有较小的波动幅度，两者具有较好的一致性。

多站多星随机周跳实验选取2018-01-12 CUT0站、2016-05-06 FTNA站、2018-04-16 XMIS站、2018-04-16 SAMO站、2018-10-01 ALGR站的北斗静态30 s采样率的观测数据，这些选中的测站均未参与岭参数$ k$值的搜索区域$ R $的确定，防止出现自洽。对选中卫星采用如上的随机实验策略，得到如表 5所示的实验结果。

实验结果表明，算法对于不同测站不同北斗卫星均同样具有良好的周跳处理效果，能准确探测出实验所有模拟随机周跳，并且除了SAMO测站C01号卫星的一个未被成功修复的历元之外，其余历元周跳均得到有效修复。进一步分析该未被修复历元，发现该历元$ {v}_{1} $值异常，形成显著孤值，扩大搜索区域之后，仅对该历元进行新一轮$ k$值搜索解算得到正确整数周跳。由于小周跳正确探测与修复的实现难度远高于大周跳的探测与修复, 本文不再进行更大周跳的探测与修复实验。

为了验证该周跳修复方法对其他同存在病态问题的优选周跳组合探测量的适用性，另选取[1 -1 0]、[0 -1 1]无几何相位组合与一个无几何伪距相位组合[1 4 -5]组成3个线性无关的周跳探测量, 其条件数达到7 066.47，法方程系数阵也存在严重的病态问题，同样采用统计-验证-优化的方式得到相应岭参数$ k$值搜索区域$ W$以及搜索步长$ {l}_{1} $, 再进行多站多星随机周跳实验。结果表明，该周跳修复算法对于此周跳探测组合同样具有良好的适用效果，使其能够有效探测和修复所有模拟随机周跳，但由于$ k $值搜索区域$ W$要大于本文讨论的搜索区域$ R$，导致搜索时间相比原组合要长，效率稍低。

4 结语

本文算法联合优选的两个无几何相位组合以及一个伪距相位组合构建了3个线性无关的周跳探测量，模拟实验表明, 该组合能够准确探测所有模拟单个以及连续周跳，无不敏感周跳。首先针对周跳修复法方程系数阵病态的问题，提出一种基于拓展的岭估计的周跳修复方法，该方法能够准确修复所有模拟周跳, 并具有较高的处理效率，通过设置$ {v}_{1}$以及$ {v}_{2}$两个周跳修复准确性检验指标，可有效确保周跳修复的可靠性。实验表明，该周跳修复算法原理可以应用到其他法方程存在病态问题的优选三频周跳探测组合中，但各类周跳探测组合之间$ k $值的联系以及$ k $值搜索区域$ R $的确定优化还需作进一步研究；周跳修复准确性指标中权阵为非单位阵情况也需作进一步考虑。本文算法实用性较强，可应用在BDS-2和BDS-3卫星三频数据实时预处理中。