周跳是影响GPS监测数据可靠性的因素之一，高精度的监测成果必须保证载波相位观测值中无周跳存在^[1-2]。周跳探测方法主要有相位多次求差法^[3]、电离层残差法^[4-6]、非差法^[7-8]、双频码相组合法^[9]、拟合法^[10-11]和TurboEdit方法^[12-14]等。上述方法在实际应用中均存在一定局限性，如相位多次求差法对小周跳探测不敏感; 双频码相组合法与电离层残差法不能确定周跳发生的频率位置及多值问题，并且对特殊组合周跳失效; 拟合法对数据个数及拟合阶数要求较高，同时也会产生拟合误差和舍入误差; TurboEdit方法中MW组合观测值易受伪距观测值噪声的影响，部分周跳容易被噪声湮没。针对上述问题，本文组建相位方程ΔΦ=ΔN₁－(λ₂/λ₁)ΔN₂，考虑到周跳在ΔΦ、ΔN₁、ΔN₂序列上表现为异常值的特点，利用Grubbs准则探测ΔN₁、ΔN₂异常值，通过约束ΔN₁、ΔN₂的范围来探测不同组合周跳。

1 双频周跳探测模型

在双频条件下，设L₁和L₂上周跳大小分别为ΔN₁、ΔN₂，根据相位和伪距观测值可得:

$ \left\{\begin{array}{l} \Delta \mathit{\Phi}=\Delta N_{1}-\frac{\lambda_{2}}{\lambda_{1}} \Delta N_{2}+\left[\Delta_{\text {ion }}(t+1)-\Delta_{\text {ion }}(t)\right] \\ \Delta N_{1}=\Delta \mathit{\Phi}_{1}-\frac{\Delta P_{1}}{\lambda_{1}}+\Delta \varepsilon_{1} \\ \Delta N_{2}=\Delta \mathit{\Phi}_{2}-\frac{\Delta P_{2}}{\lambda_{2}}+\Delta \varepsilon_{2} \end{array}\right. $

(1)

式中，ΔΦ为电离层残差在历元间的变化值; λ为波长; Δ_ion(t)为某时刻的电离层残差; Δε为总误差项，主要受电离层和多路径影响。

若采样间隔较小，电离层在历元间的变化Δ_ion(t)很小，则式(1)可简化为:

$ \left\{\begin{array}{l} \Delta \mathit{\Phi}=\Delta N_{1}-\frac{77}{60} \Delta N_{2} \\ \Delta N_{1}=\Delta \mathit{\Phi}_{1}-\frac{\Delta P_{1}}{\lambda_{1}} \\ \Delta N_{2}=\Delta \mathit{\Phi}_{2}-\frac{\Delta P_{2}}{\lambda_{2}} \end{array}\right. $

(2)

从式(2)可以看出，针对2个未知数ΔN₁、ΔN₂，1个方程无法解算出周跳值，因此可通过计算附有限制条件ΔΦ=ΔN₁－(λ₂/λ₁)ΔN₂来确定ΔN₁、ΔN₂的取值区间。

2 双频周跳特征分析

假设ΔΦ服从正态 $ N\left( {0,\sigma _{\Delta \mathit{\Phi} }^2} \right)$分布，在观测过程中相位观测值的中误差均为0.01周，则σ_ΔΦ=0.023周，按3倍中误差理论计算，则ΔΦ的限差为|3σ_ΔΦ|，即0.069周。因此在无周跳时，式(2)中ΔΦ=ΔN₁－(λ₂/λ₁)ΔN₂应为[-0.07, 0.07]。如超出该取值区间，可认为存在周跳，因此可利用该值进行周跳探测，结果记为一般周跳。若ΔN₁－(λ₂/λ₁)ΔN₂≈0，即ΔN₁=9 k、ΔN₂=7 k或ΔN₁=77 k、ΔN₂=60 k(k∈z)时，采用[-0.07, 0.07]约束ΔΦ=ΔN₁－(λ₂/λ₁)ΔN₂无法检测出周跳，记为(9, 7)型和(77, 60)型特殊周跳。因此周跳值均可看作ΔΦ、ΔN₁、ΔN₂序列的异常值，探测出的异常值则为周跳值。本文基于周跳特征提出基于Grubbs准则和区间搜索的周跳探测方法。

3 Grubbs检验法

设存在观测序列Y=(Y₁, Y₂, …, Y_n)，构建Grubbs检验统计量:

$ G = \frac{{\mathop {\max }\limits_{i = 1, 2, \cdots , n} \left| {{Y_i} - \bar Y} \right|}}{S} $

(3)

式中，$ {\bar Y}$和S分别为序列Y的均值和标准差。

Grubbs检验统计量为样本标准偏差单位与样本均值的最大绝对偏差。在显著性水平α下，若

$ G>\frac{n-1}{\sqrt{n}} \sqrt{\frac{t_{\alpha /(2 n), n-2}^{2}}{n-2+t_{\alpha /(2 n), n-2}^{2}}} $

(4)

表明存在异常值，即存在周跳。式中，t_{α/(2n), n－2}表示具有(n-2)自由度和显著性水平α/(2n)的T分布临界值。

4 区间搜索法周跳探测流程

对于式(2)，由于周跳ΔN₁、ΔN₂的整数特性，若能限制ΔN₁、ΔN₂范围，则可利用$\Delta \mathit{\Phi} = \Delta {N_1} - \frac{{77}}{{60}}\Delta {N_2} $解算出ΔN₁、ΔN₂，即周跳值。由于特殊周跳的成对出现特性，ΔN₁、ΔN₂序列存在异常值，而与之对应历元位置的ΔΦ未出现异常值。各类型周跳都需探测出序列中的异常值，而异常值探测的最核心问题为ΔN₁、ΔN₂的整数区间。如果可确定搜索区间，在区间内满足式(2)中$\Delta \mathit{\Phi} = \Delta {N_1} - \frac{{77}}{{60}}\Delta {N_2} $的整数解就为周跳值。区间搜索法周跳探测流程为:

1) 组建观测序列ΔΦ、ΔN₁、ΔN₂。

2) 判断序列ΔΦ是否超出范围[-0.07, 0.07]。若是，则ΔΦ存在普通周跳，记录ΔΦ对应的位置ti、ΔΦ(ti)，并执行3);若否，则ΔΦ无普通周跳，执行4)。

3) 剔除ti位置对应的ΔN₁(ti)、ΔN₂(ti)，将余下序列组成新序列ΔYN₁、ΔYN₂，利用Grubbs准则分别探测新序列ΔYN₁、ΔYN₂是否存在奇异值; 若是，分别探测ΔYN₁、ΔYN₂奇异值的位置及t1i、ΔYN₁(t1i)、t2i、ΔYN₂(t2i)大小，并计算剔除ΔYN₁、ΔYN₂中奇异值后的序列中误差m₁、m₂，然后执行5);若否，则无特殊周跳。

4) 利用Grubbs准则分别探测序列ΔN₁、ΔN₂是否存在奇异值; 若是，分别探测ΔN₁、ΔN₂奇异值的位置及t1i、ΔN₁(t1i)、t2i、ΔN₂(t2i)大小，并计算剔除ΔN₁、ΔN₂奇异值后的序列中误差m₁、m₂，然后执行5);若否，则无特殊周跳。

5) 判断t1i、t2i是否存在相同值t12i; 若是，则存在特殊周跳，计算剔除ti、t1i、t2i位置ΔΦ后的中误差m，分别执行6)~8)以探测普通周跳、(9, 7)型特殊周跳和(77, 60)型特殊周跳; 若否，则无特殊周跳。

6) 计算普通周跳:

$ \left\{\begin{array}{l} \left|\Delta N_{1}-\frac{77}{60} \Delta N_{2}-\Delta {\Phi}(t i)\right|<k m, k=3 \\ \Delta N_{1} \in \mathrm{INT} \\ {\left[\Delta N_{1}(t i)-3 m_{1} \quad \Delta N_{1}(t i)+3 m_{1}\right]} \\ \Delta N_{2} \in \mathrm{INT} \\ {\left[\Delta N_{2}(t i)-3 m_{2} \quad \Delta N_{2}(t i)+3 m_{2}\right]} \end{array}\right. $

(5)

7) 计算(9, 7)型特殊周跳:

$ \left\{\begin{array}{l} \left|\Delta N_{1}-\frac{77}{60} \Delta N_{2}\right|<0.07 \\ \Delta N_{1} \in \mathrm{INT} \\ {\left[\Delta N_{1}(t 12 i)-3 m_{1} \quad \Delta N_{1}(t 12 i)+3 m_{1}\right]} \\ \Delta N_{2} \in \mathrm{INT} \\ {\left[\Delta N_{2}(t 12 i)-3 m_{2} \quad \Delta N_{2}(t 12 i)+3 m_{2}\right]} \end{array}\right. $

(6)

8) 计算(77, 60)型特殊周跳:

$ \left\{\begin{array}{l} |\Delta N_{1}-\frac{77}{60} \Delta N_{2}|=0 \\ \Delta N_{1} \in \mathrm{INT} \\ [\Delta N_{1}(t 12 i)-3 m_{1} & \Delta N_{1}(t 12 i)+3 m_{1} ] \\ \Delta N_{2} \in \mathrm{INT} \\ [\Delta N_{2}(t 12 i)-3 m_{2} & \Delta N_{2}(t 12 i)+3 m_{2}] \end{array}\right. $

(7)

5 实验验证分析

本文对提出的基于Grubbs准则和区间搜索的周跳探测方法进行验证，实验数据采样率为1 s。在原采样率1 s数据基础上，提取采样率为5 s、15 s的双频观测数据进行分析。实验时人为加入连续或非连续的不同组合周跳(ΔN₁、ΔN₂)，图 1~6为加入周跳前后ΔΦ、ΔN₁、ΔN₂序列对比，探测结果与直接取整法^[15]对比见表 1~3(单位周)。

从图 1、图 3和图 5可以看出，加入周跳前后序列波动发生明显变化，周跳大小对序列ΔΦ、ΔN₁、ΔN₂均有影响; 从图 2、图 4和图 6可以看出，本文方法能准确剔除普通周跳和特殊型周跳。

从表 1~3可以看出，本文方法探测的周跳与人为加入的周跳完全一致; 而取整法能探测出位置，但准确性较低。随着采样间隔的增加，ΔΦ、ΔN₁、ΔN₂剔除周跳后中误差m、m₁、m₂均增大，但本文方法仍能准确探测周跳，取整法在双频上会产生约1~2周误差，为原噪声湮没所致。

6 结语

本文基于Grubbs准则和区间搜索进行双频GPS周跳探测，得出以下结论:

1) 以ΔΦ=ΔN₁－(λ₂/λ₁)ΔN₂为主方程，同时利用周跳的整数特性，通过解附有限制条件方程可确定唯一周跳值，避免多值问题。

2) 通过Grubbs准则判断ΔN₁、ΔN₂异常值，剔除异常值后计算序列中误差，利用3倍中误差理论约束双频周跳的搜索范围，可避免噪声湮没问题。

3) 本文方法能准确探测出周跳发生的历元、各频率上周跳大小，同时能准确探测普通型周跳和(9, 7)型、(77, 60)型特殊周跳，以及连续周跳和非连续周跳。

4) 实验数据分析表明，对于采样率为15 s以内的数据，本文方法探测准确率高，计算成功率为100%。