利用载波相位观测值代替伪距与载波相位组合的方法构建GF组合，并对其进行相邻历元间求差可解得GF历元差^[6-9]，公式如下：

式中，$ {L}_{\mathrm{G}\mathrm{F}}(i+1) $、$ {L}_{\mathrm{G}\mathrm{F}}\left(i\right) $分别为$ i+1 $和$ i $历元的GF组合观测值；$ \mathrm{\Delta }A $为历元间电离层延迟变化；$ {f}_{1} $、$ {f}_{2} $表示对应频段的频率；$ {\lambda }_{1}\mathrm{\Delta }{N}_{1}-{\lambda }_{2}\mathrm{\Delta }{N}_{2} $为周跳影响项，在周跳未发生时，其理论值为0；$ {e}_{12}\left(i\right) $为相位噪声之间的差值。

历元间电离层延迟$ \mathrm{\Delta }A $主要取决于沿信号传播路径的总电子含量（total electron contents，TEC）与数据采样率的大小，当电离层延迟变化平稳时，可视其为低频信号并随时间平稳变化^[10]，给定严格的检验阈值，可以有效地分辨一周的小周跳。但是，当出现电离层闪烁时，会使得$ \mathrm{\Delta }A $产生剧烈的波动，从而使$ {e}_{12}\left(i\right) $变大，这两种误差的影响无法通过一阶或高阶差分的方法消除，使得式（1）失去对小周跳的敏感性，并导致高误判与漏判。基于此，本文采用滑动多项式拟合法提取历元间电离层延迟与残余多路径效应误差的趋势项，以抑制两种误差的影响。

对子区间$ [{t}_{s}, {t}_{i-1}] $进行多项式拟合，令$ {Q}_{k}\left(\beta \right)={\beta }_{0}+\sum\limits_{j=1}^{m}{\beta }_{j}{t}_{k}^{j} $，根据最小二乘原理，参数$ {\mathit{\boldsymbol{\beta }}} $的估计准则定义如下：

式中，$ {\mathit{\boldsymbol{P}}} $为历元差观测值的权阵；$ \left\{\mathrm{\Delta }{L}_{\mathrm{G}\mathrm{F}}\right\} $为时间序列，由$ \left\{{L}_{\mathrm{G}\mathrm{F}}\right\} $经过一次差获得，相比于原序列，其变化较平稳。根据试验分析，当给定时间间隔10~15 min，多项式系数$ j $取2或3时，即可满足去趋势的要求。

利用式（2）求解得到多项式参数$ \widehat{{\mathit{\boldsymbol{\beta }}}} $，则可得$ {t}_{i} $历元处的预测$ {Q}_{i}\left({\mathit{\boldsymbol{\beta }}}\right)={\widehat{\beta }}_{0}+\sum\limits_{j=1}^{m}{\widehat{\beta }}_{j}{t}_{i}^{j} $，则周跳检测量可表示为：

若$ D\left(i\right)>n\sigma \left(i\right) $，则认为$ i $历元出现了周跳或粗差，其中$ n $为整数，$ \sigma \left(i\right) $为单位权中误差。

阈值系数$ n $的选择是周跳能否被成功探测的关键，选取IGS跟踪站shao 2018年年积日第275天中G15卫星的一段干净数据进行周跳试验，在低卫星高度角与电离层延迟活跃历元加入较小的难以探测的周跳，根据文献[11]中的数据探测法原理，取$ n=4 $进行周跳探测，结果如图 1所示。

图  1  $ D $检测量周跳探测结果(n=4)

Figure 1.  Cycle Slip Detection Results of D Detector(n=4)

从图 1中可以看出，无周跳时的$ D $检测量是一个极小的值，有利于小周跳的探测。但$ n $取4时，$ D $检测量会出现周跳的误探现象，特别是在100、500和800等受卫星高度角和电离层延迟变化影响较大的历元附近，周跳误探现象严重。

在电离层延迟变化平稳的假设下，式（3）中的残差$ e $具有方差齐性，$ 4\sigma $可以用于判断是否发生了周跳。而当电离层活跃或出现电离层闪烁时，残差$ e $的分布会出现尖峰长尾的特征，此时$ 4\sigma $将会引起高误判。但是，周跳对$ D $检测量的影响远大于电离层闪烁，因此，可采取两个途径解决这一问题：（1）设定固定阈值，如0.05；（2）增大$ n $。

考虑到固定阈值的算法适应性较差，本文选择增大$ n $的策略。为了避免$ n $增大后会引入不敏感周跳，同时要保证$ n\sigma $能够成功消除电离层闪烁的影响，所以$ n $的选取需满足下式：

式中，$ {D}_{I} $与$ {D}_{c} $分别为电离层闪烁与最小周跳组合对$ D $探测量的影响大小。选择文献[12]中磁暴发生期间的数据，根据上述算法计算周跳检测量$ D $与$ n $的取值，其中受磁暴影响较大的G13、G15、G25、G26号卫星的计算结果如表 1所示。

根据表 1，选取$ 8\sigma $作为周跳探测阈值，图 2为$ n=8 $时G25卫星的周跳探测结果。从图 2中可以看出，$ n $取8时，可以准确探测出（1，1）这种小周跳组合，且不存在周跳的误探现象。同时可以预见，当数据质量较好时，选择$ n=8 $的探测效果将会更好。

图  2  $ D $检测量周跳探测结果(n=8)

Figure 2.  Cycle Slip Detection Results of $ D $Detector(n=8)

为了弥补两个频段发生特殊周跳，即出现$ {\lambda }_{1}\mathrm{\Delta }{N}_{1}-{\lambda }_{2}\mathrm{\Delta }{N}_{2}=0 $，GF组合失效的问题，本文将MW组合作为GF组合的辅助探测，只用于探测GF组合的探测盲点。

当数据中含有异常时，可将异常归入函数模型，把含有异常的数据看作与其他正常数据相比具有不同期望、相同方差的数据点，可构建均值漂移模型为^[13-14]：

式中，$ \gamma $表示第$ i $点处的异常大小，若$ \gamma $显著不等于0，则说明第$ i $点处存在异常，因此基于均值漂移模型的假设检验问题为：

对于式（5）中的均值漂移模型，若$ \Delta \sim N(0, {\sigma }^{2}) $，则假设检验问题（6）的Score统计量可表示为^[15-16]：

式中，$ \sigma $、$ {v}_{i} $分别为数据拟合中误差、拟合残差；$ {h}_{ii} $为帽子矩阵$ {\mathit{\boldsymbol{H}}}={\mathit{\boldsymbol{B}}}({{\mathit{\boldsymbol{B}}}}^{\mathrm{T}}{{\mathit{\boldsymbol{B}}})}^{-1}{{\mathit{\boldsymbol{B}}}}^{\mathrm{T}} $的第$ i $个对角线元素；$ \alpha $为显著性水平，若$ {S}_{i}>{\chi }_{\alpha }^{2}\left(1\right) $，则认为此点为异常点。

根据式（1），时间序列$ \left\{\mathrm{\Delta }{L}_{\mathrm{G}\mathrm{F}}\right\} $由$ \left\{{L}_{\mathrm{G}\mathrm{F}}\right\} $经过一次差获得，使得基于GF历元差模型的各种算法无法分辨数据中的周跳与粗差。通过对式（1）的分析，本文将$ D $检测量与MW组合探测出的历元标记为可疑历元，并利用探测出的连续可疑历元个数（此处的连续是指可疑历元是否处于相邻历元），将周跳与粗差分为以下3种情况：①无连续可疑历元；②连续2个可疑历元；③连续多个可疑历元（3个及以上）。对于情况①，可直接判定此历元为周跳历元；对于情况③，表明此段数据中出现了多维连续粗差或周跳，此时无需对粗差或周跳进行区分，而是直接对此段数据进行剔除、分段；当情况②发生时，表明数据中出现了单个不连续粗差或连续2个周跳，此时若不对粗差与周跳加以鉴别，将会大大降低数据利用率，并导致周跳的误判。

由式（7）可以看出，Score统计量的大小取决于拟合残差与拟合中误差的大小，其对异常点具有很高的灵敏度，因此可以利用Score检验区别情况②中的周跳与粗差。

考虑到实时周跳探测过程中只能使用当前历元与历元前的数据，本文对Score检验进行如下改进：（1）首先对$ i $历元前$ n $个GF组合值进行多项式拟合，求得多项式系数的参数解$ \widehat{{\mathit{\boldsymbol{X}}}} $和拟合中误差$ \sigma $；（2）根据$ \widehat{{\mathit{\boldsymbol{X}}}} $，求解$ i $与$ i+1 $历元的残差$ {v}_{i} $、$ {v}_{i+1} $与$ {h}_{ii} $、$ {h}_{i+1, i+1} $；（3）将（2）中的残差$ v $与杠杆值$ h $代入式（7），求解$ i $与$ i+1 $历元的Score值，若$ {S}_{i}>{\chi }_{\alpha }^{2} $，则此历元未通过检验。

从上述改进方法可以看出，利用系数解$ \widehat{{\mathit{\boldsymbol{X}}}} $求解出的$ i $与$ i+1 $历元的残差$ {\mathit{\boldsymbol{v}}} $实际上是预测残差，但$ i+1 $历元的预测值与实际值可能出现较大的差异。为消除预测带来的影响，第一次检验完成后，需再进行一次预测值偏离检验，预测值偏离检验设置如下：

式（8）用于消除预测值偏离的影响，式（9）则是在固定阈值失效的情况下，保证改进的Score检验仍能区分周跳与粗差。

综上所述，基于改进的Score检验与预测值偏离检验的周跳与粗差分离算法如下：（1）对于$ D $检测量与MW组合探测出的连续两个可疑历元$ i $与$ i+1 $，选定包括连续两个历元在内的前$ n+2 $个GF历元数据进行改进的Score检验；（2）若Score检验未通过个数为1，则$ i $历元为粗差历元；（3）若Score检验未通过个数为2，则根据式（8）与式（9）进一步进行预测值偏离检验，若检验后未通过个数为1，则$ i $历元为粗差历元，若为2，则$ i $与$ i+1 $历元为连续周跳历元。

上述算法只利用了待探测历元前的数据，适用于实时周跳与粗差的分离。

通过联合$ D $检测量与MW组合探测周跳，并利用基于Score检验的周跳与粗差分离模型剔除粗差后，将剩下历元标记为周跳发生历元，即可进行周跳修复。联合式（3）与MW组合宽巷周跳值对两个频段上的周跳值进行求解，分别为：

式（10）可简写为$ {\mathit{\boldsymbol{L}}}={\mathit{\boldsymbol{A}}}{\mathit{\boldsymbol{N}}} $，理论上直接对式（10）的解$ {\mathit{\boldsymbol{N}}} $进行取整运算即可获得周跳估值$ \widehat{{\mathit{\boldsymbol{N}}}} $，但当卫星高度角较低或电离层延迟变化活跃时，MW组合的精度会受到极大的影响，使得周跳引起的$ \mathrm{\Delta }{N}_{\mathrm{M}\mathrm{W}} $与理论值相差过大，在利用$ \widehat{{\mathit{\boldsymbol{N}}}}=({{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}{{\mathit{\boldsymbol{A}}})}^{-1}{{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}{\mathit{\boldsymbol{L}}} $解算周跳时，$ {\mathit{\boldsymbol{N}}} $的中误差较大，通过取整运算，会获得错误的周跳值。

借鉴多频周跳修复中的搜索思想^[17]，考虑到影响双频周跳修复的主要因素为$ \mathrm{\Delta }{{\mathit{\boldsymbol{N}}}}_{\mathrm{M}\mathrm{W}} $的误差过大，其误差约为两周周跳^[7]，无法通过有效手段予以消除。因此，本文以$ {\mathit{\boldsymbol{N}}} $为搜索中心、一周周跳为搜索步长，以两周周跳确定一个搜索范围$ \Omega $，并鉴于本文周跳探测量的灵敏性，不成功的周跳修复会使得$ {\mathit{\boldsymbol{L}}} $的1范数急剧增大，正确周跳修复造成的探测量变化值的1范数应该最小，周跳修复值$ \widehat{{\mathit{\boldsymbol{N}}}} $的选取如下：

采用IGS网站下载的上海站2018年年积日第275天的数据对本文算法进行验证，选取G15卫星，共计870个历元，采样间隔为30 s。为验证本文周跳探测与修复算法的效果，在低卫星高度角与电离层延迟变化活跃的历元加入较小的难以探测的周跳与粗差，具体模拟情况见表 2。

本文采用两种方案对表 2中的周跳与粗差进行探测，方案1：周跳探测的历元差阈值模型法^[9]；方案2：本文改进的历元差模型法。

考虑到本文将精度较低的MW组合作为GF组合的辅助探测，因此上述两种方案中的MW组合都采用滑动窗口滤波模型^[9]。GF组合周跳探测结果见图 3。

图  3  GF组合周跳探测结果

Figure 3.  Cycle Slip Detection Results of GF Combination

在图 3（a）中，采用历元差阈值模型法对周跳进行探测，对于30 s采样率的数据，其最小阈值设置为0.15 m，只探测出了300历元处较大的周跳。而在图 3（b）中，通过在历元差模型中加入滑动多项式拟合法提取历元间电离层延迟与多路径效应误差的趋势项，使得$ D $检测量对小周跳具有较高的探测灵敏度，并通过8倍中误差的阈值选择，可以准确探测出50、100、200、300、400、401、500、700和800历元处发生的周跳与粗差。

由图 3可知，对于在200、500与700历元处模拟的小粗差，GF历元差组合无法区别周跳与粗差，使得粗差的出现会造成周跳的误探现象。此外，由于GF历元差是对GF组合相邻历元进行求差获得的，在200、500与700历元出现的粗差会使得201、501与701历元也被当作周跳历元探测出来，因此，粗差的出现会极大地增加周跳的误探。

由式（3）可知，当两个频段上出现的周跳值近似满足$ {\lambda }_{1}\mathrm{\Delta }{N}_{1}-{\lambda }_{2}\mathrm{\Delta }{N}_{2}=0 $时，GF组合无法探测出600历元处的特殊周跳，为此，将MW组合作为GF组合的辅助探测，用于探测特殊周跳组合，探测结果如图 4所示，$ \mathrm{\Delta } $MW为MW组合周跳探测量。

图  4  MW组合周跳探测结果

Figure 4.  Cycle Slip Detection Results of MW Combination

由图 4可知，在历元前段与后段，$ \mathrm{\Delta }\mathrm{M}\mathrm{W} $波动比较剧烈，这是因为此时卫星高度角较低，受多路径效应误差与伪距噪声的影响较大，MW组合探测精度不高。对于600历元处的特殊周跳，由于其所处的观测时段较好，MW组合能够准确将其探测出来，若小于两周的非等周周跳发生在低卫星高度角时段，则MW组合是无法将其探测出来的。因本文只是将MW组合用于探测（9，7）等较大的特殊周跳，所以MW组合满足辅助探测的需要。

综上所述，通过联合改进的GF历元差组合与MW组合可以准确探测出表 2中模拟的所有小周跳、特殊周跳、等周周跳、大周跳与连续周跳，探测性能比较稳定和高效。

上述改进的历元差模型无法区别周跳与粗差，但探测出了连续可疑历元200、201、400、401、500、501、700和701，符合本文周跳与粗差分离算法中的情况②，因此采用基于Score检验的周跳与粗差分离算法进行周跳与粗差的分离。周跳与粗差的分离结果见表 3。

从表 3可知，若取显著性水平$ \alpha =0.01 $，临界值$ {{\chi }_{1}}^{2}\left(0.01\right)=6.634\mathrm{ }9 $，则200、201、400、401、500、700、701未通过检验，拒绝假设问题（6）中的$ {H}_{0} $假设，说明这些历元的$ \gamma $显著不等于0，因而为异常历元。此处的异常是由周跳与粗差和预测值偏离引起的，因此需再进行一次预测值偏离检验。

以200历元处加入的粗差与400、401历元加入的连续周跳为例，虽然粗差只出现在200历元，但历元差模型却会将201历元也当做周跳探测出来。通过改进的Score检验，200历元的Score值远大于201历元，并利用预测值偏离检验消除201历元的影响，检验未通过个数为1。同理，400、401历元发生的周跳使得其Score值都过大，未通过个数为2，符合本文周跳与粗差分离的设计情况。检验最后结果判断出200、500、700历元为粗差历元，400、401历元为连续周跳历元。因此，本文基于Score检验的周跳与粗差分离模型可以有效解决历元差模型无法区别周跳与粗差的问题，避免了周跳的误判，提高了周跳探测的可靠性。

在利用改进的历元差模型与MW组合进行周跳探测，并通过Score检验模型剔除粗差后，将剩下历元标记为周跳历元，并利用式（10）进行周跳值解算，以式（11）作为周跳修复成功的判断准则，周跳修复结果见表 4。

在周跳修复过程中，以800历元处的周跳组合（1，1）为例，该组合对$ \mathrm{\Delta }\mathrm{M}\mathrm{W} $的理论影响值为0，但在表 4中可以看出，$ \mathrm{\Delta }\mathrm{M}\mathrm{W} $达到了0.876 7周，与理论值相差过大，受其影响，800历元周跳值取整结果为（3，2），与真值相差两周，周跳修复错误。

采用本文周跳修复方法，以（3，2）为搜索中心，以两周周跳确定搜索范围，并通过式（11）求解1范数，搜索出的最小1范数为0.879 2，此时对应的周跳组合为（1，1），可以确定原始观测值中在800历元处发生的周跳为（1，1）。同理，对于50、100、300、401历元，其周跳值取整结果与实际值不符，通过搜索使$ \left|\right|{\mathit{\boldsymbol{L}}}|{|}_{1} $达到最小的周跳组合，成功纠正了周跳修复错误。

在周跳探测方面，改进的GF历元差法能够削弱电离层延迟变化和多路径趋势项的影响，提高了周跳探测的精度。同时，通过合适的阈值选择，即使在高电离层闪烁的情况下，也能实现小周跳、连续周跳的准确探测。

在周跳与粗差分离方面，本文将拟合残差替换为预测残差改进Score检验，使其适用于实时情况下的异常检验。基于此构建的周跳与粗差分离模型可准确将粗差从周跳中成功分离出来，避免了周跳的误判，提高了实时TurboEdit周跳探测方法的可靠性。

在周跳修复方面，本文算法以正确周跳修复造成的探测量变化值$ {\mathit{\boldsymbol{L}}} $的1范数最小作为周跳修复成功的判断准则。与直接求解周跳解并取整相比，该算法不仅能够验证周跳解的正确性，还可对错误的周跳解进行实时纠正，得到准确的周跳值，对历元后续观测值进行修复。