-
近年来,精密单点定位(precise point positioning,PPP)逐渐成为卫星导航定位中的研究热点,在PPP数据预处理中,周跳探测与修复是实现高精度定位和参数快速收敛的关键问题[1]。
目前,非差观测值的周跳探测方法主要有高次差法[2]、多项式拟合法[2]、卡尔曼滤波法[3]、TurboEdit法[4]等,其中文献[4]提出的TurboEdit法具有单站探测、不依赖于卫星和接收机的运动状态、不受钟差以及对流层延迟影响等优点,特别适用于非差数据的周跳探测[5-7]。但TurboEdit法容易受到伪距噪声、高电离层延迟和低高度角的影响,只能探测较大的周跳,为此,不少学者针对TurboEdit法的缺陷,提出了一些改进的方法。文献[5]通过估算总电子含量变化率(total electron contents rate,TECR)来削弱电离层延迟的影响,并联合Melbourne-Wübbena(MW)组合探测周跳;文献[6]提出了联合前后向滑动窗口的MW组合和二次历元间差分的电离层残差组合(second-order time difference phase ionospheric residual,STPIR)进行周跳探测的方法;文献[7]为处理数据中频繁出现的小周跳,提出了稳健的多项式拟合算法以获得自适应周跳检测阈值,提高了对小周跳的探测能力。然而,这些TurboEdit改进算法在周跳发生前后,需要几分钟连续并且干净的相位数以满足周跳探测的需要,因此,它们并不适用于实时PPP[8]。文献[9]在载波相位无几何距离(geometry-free,GF)组合中加入历元差模型,解决了TurboEdit法无法用于实时周跳探测的问题。但历元差模型存在以下问题:(1)周跳探测采用经验阈值,算法适应性较差;(2)历元差模型无法分辨周跳与粗差,容易将粗差误探为周跳。另外,由于信号噪声的影响,周跳难以准确修复。
为了解决上述TurboEdit法进行实时周跳探测与修复的问题,本文首先在GF历元差模型中加入滑动窗口拟合算法,以削弱电离层延迟与多路径效应对周跳探测的影响;然后引入基于均值漂移模型的Score统计量,并构建基于Score检验的周跳与粗差分离模型,避免了粗差存在时的周跳误探现象;最后联合改进的GF历元差模型与MW组合进行周跳修复,在修复过程中使用探测量最小1范数准则来选取正确的周跳值。
-
利用载波相位观测值代替伪距与载波相位组合的方法构建GF组合,并对其进行相邻历元间求差可解得GF历元差[6-9],公式如下:
$$ \begin{array}{c}\mathrm{\Delta }{L}_{\mathrm{G}\mathrm{F}}\left(i\right)={L}_{\mathrm{G}\mathrm{F}}(i+1)-{L}_{\mathrm{G}\mathrm{F}}\left(i\right)=\\ \frac{\mathrm{\Delta }A({f}_{2}^{2}-{f}_{1}^{2})}{{f}_{1}^{2}{f}_{2}^{2}}+{\lambda }_{1}\mathrm{\Delta }{N}_{1}-{\lambda }_{2}\mathrm{\Delta }{N}_{2}+{e}_{12}\left(i\right)\end{array} $$ (1) 式中,$ {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 }}} $的估计准则定义如下:
$$ \mathrm{m}\mathrm{i}\mathrm{n}\left\{\right(\mathrm{\Delta }{L}_{\mathrm{G}\mathrm{F}}-Q{\left({\mathit{\boldsymbol{\beta }}}\right))}^{\mathrm{T}}{\mathit{\boldsymbol{P}}}(\mathrm{\Delta }{L}_{\mathrm{G}\mathrm{F}}-Q({\mathit{\boldsymbol{\beta }}}\left)\right)\} $$ (2) 式中,$ {\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} $,则周跳检测量可表示为:
$$ \begin{array}{c}D\left(i\right)=\left|\mathrm{\Delta }{L}_{\mathrm{G}\mathrm{F}}\right(i)-Q(i\left)\right|={\lambda }_{1}\mathrm{\Delta }{N}_{1}\left(i\right)-\\ \mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }{\lambda }_{2}\mathrm{\Delta }{N}_{2}\left(i\right)+{e}_{D}\end{array} $$ (3) 若$ 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时,$ D $检测量会出现周跳的误探现象,特别是在100、500和800等受卫星高度角和电离层延迟变化影响较大的历元附近,周跳误探现象严重。
在电离层延迟变化平稳的假设下,式(3)中的残差$ e $具有方差齐性,$ 4\sigma $可以用于判断是否发生了周跳。而当电离层活跃或出现电离层闪烁时,残差$ e $的分布会出现尖峰长尾的特征,此时$ 4\sigma $将会引起高误判。但是,周跳对$ D $检测量的影响远大于电离层闪烁,因此,可采取两个途径解决这一问题:(1)设定固定阈值,如0.05;(2)增大$ n $。
考虑到固定阈值的算法适应性较差,本文选择增大$ n $的策略。为了避免$ n $增大后会引入不敏感周跳,同时要保证$ n\sigma $能够成功消除电离层闪烁的影响,所以$ n $的选取需满足下式:
$$ \mathrm{m}\mathrm{a}\mathrm{x}{D}_{I}<n\sigma <\mathrm{m}\mathrm{i}\mathrm{n}{D}_{c} $$ (4) 式中,$ {D}_{I} $与$ {D}_{c} $分别为电离层闪烁与最小周跳组合对$ D $探测量的影响大小。选择文献[12]中磁暴发生期间的数据,根据上述算法计算周跳检测量$ D $与$ n $的取值,其中受磁暴影响较大的G13、G15、G25、G26号卫星的计算结果如表 1所示。
表 1 阈值系数$ {\mathit{\boldsymbol{n}}} $选取结果
Table 1. Selection Results of Threshold Coefficient n
根据表 1,选取$ 8\sigma $作为周跳探测阈值,图 2为$ n=8 $时G25卫星的周跳探测结果。从图 2中可以看出,$ n $取8时,可以准确探测出(1,1)这种小周跳组合,且不存在周跳的误探现象。同时可以预见,当数据质量较好时,选择$ n=8 $的探测效果将会更好。
为了弥补两个频段发生特殊周跳,即出现$ {\lambda }_{1}\mathrm{\Delta }{N}_{1}-{\lambda }_{2}\mathrm{\Delta }{N}_{2}=0 $,GF组合失效的问题,本文将MW组合作为GF组合的辅助探测,只用于探测GF组合的探测盲点。
-
当数据中含有异常时,可将异常归入函数模型,把含有异常的数据看作与其他正常数据相比具有不同期望、相同方差的数据点,可构建均值漂移模型为[13-14]:
$$ \left\{\begin{array}{l}{L}_{j}={b}_{j}^{\mathrm{T}}{X}_{j}+{\Delta }_{j}, j=\mathrm{1, 2}\cdots n, \mathrm{且}j\ne i\\ {L}_{i}={b}_{i}^{\mathrm{T}}{X}_{i}+\gamma +{\Delta }_{i}, j=i\end{array}\right. $$ (5) 式中,$ \gamma $表示第$ i $点处的异常大小,若$ \gamma $显著不等于0,则说明第$ i $点处存在异常,因此基于均值漂移模型的假设检验问题为:
$$ \left\{\begin{array}{c}{H}_{0}:\gamma =0\\ {H}_{1}:\gamma \ne 0\end{array}\right. $$ (6) 对于式(5)中的均值漂移模型,若$ \Delta \sim N(0, {\sigma }^{2}) $,则假设检验问题(6)的Score统计量可表示为[15-16]:
$$ {S}_{i}=\frac{{v}_{{}_{i}}^{2}}{{\sigma }^{2}(1-{h}_{ii})}\sim {\chi }_{\alpha }^{2}\left(1\right) $$ (7) 式中,$ \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 $历元的预测值与实际值可能出现较大的差异。为消除预测带来的影响,第一次检验完成后,需再进行一次预测值偏离检验,预测值偏离检验设置如下:
$$ {S}_{i+1}>500 $$ (8) $$ \frac{{S}_{i}}{{S}_{i+1}}<3 $$ (9) 式(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组合宽巷周跳值对两个频段上的周跳值进行求解,分别为:
$$ \left[\begin{array}{c}D\\ \mathrm{\Delta }{N}_{\mathrm{M}\mathrm{W}}\end{array}\right]=\left[\begin{array}{cc}{\lambda }_{1}& -{\lambda }_{2}\\ \frac{c}{{f}_{1}-{f}_{2}}& -\frac{c}{{f}_{1}-{f}_{2}}\end{array}\right]\left[\begin{array}{c}\mathrm{\Delta }{N}_{1}\\ \mathrm{\Delta }{N}_{2}\end{array}\right] $$ (10) 式(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}}}} $的选取如下:
$$ \underset{\widehat{{\mathit{\boldsymbol{N}}}}\in \Omega }{\widehat{{\mathit{\boldsymbol{N}}}}=\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}\left(\right|\left|{\mathit{\boldsymbol{L}}}\right|{|}_{1})} $$ (11) -
采用IGS网站下载的上海站2018年年积日第275天的数据对本文算法进行验证,选取G15卫星,共计870个历元,采样间隔为30 s。为验证本文周跳探测与修复算法的效果,在低卫星高度角与电离层延迟变化活跃的历元加入较小的难以探测的周跳与粗差,具体模拟情况见表 2。
表 2 周跳与粗差模拟情况/周
Table 2. Situation of Cycle Slip and Gross Errors Simulation/cycle
历元 周跳 粗差 50 (1, 1) 100 (5, 4) 200 (1, 1) 300 (1, 0) 400 (-1, -1) 401 (-1, -1) 500 (1, 1) 600 (-9, -7) 700 (0, 0.5) 800 (1, 1) -
本文采用两种方案对表 2中的周跳与粗差进行探测,方案1:周跳探测的历元差阈值模型法[9];方案2:本文改进的历元差模型法。
考虑到本文将精度较低的MW组合作为GF组合的辅助探测,因此上述两种方案中的MW组合都采用滑动窗口滤波模型[9]。GF组合周跳探测结果见图 3。
在图 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可知,在历元前段与后段,$ \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 周跳与粗差分离结果
Table 3. Separation Results of Cycle Slip andGross Errors
历元 Score值 未通过
个数是否进行预测值偏离检验 未通过
个数检验
结果200 16 424.315 5 2 是 1 粗差 201 170.242 9 400 5 035.961 8 2 是 2 连续
周跳401 69 444.983 8 500 275.990 6 1 否 1 粗差 501 2.387 8 700 4 210.110 3 2 是 1 粗差 701 435.842 1 从表 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。
表 4 周跳修复结果
Table 4. Results of Cycle Slip Repair
历元 周跳组合 $ \mathrm{\Delta }\mathrm{M}\mathrm{W} $ 周跳浮点解 取整结果 最小1范数 估计周跳 50 (1, 1) 0.270 2 (2.209 3, 1.939 1) (2, 2) 0.271 0 (1, 1) 100 (5, 4) 1.585 4 (7.776 5, 6.191 1) (8, 6) 0.897 8 (5, 4) 300 (1, 0) 0.821 2 (0.342 1, -0.479 1) (0, 0) 0.180 4 (1, 0) 400 (-1, -1) 0.031 0 (-0.837 3, -0.868 3) (-1, -1) 0.032 2 (-1, -1) 401 (-1, -1) 0.186 9 (-0.086 4, -0.273 3) (0, 0) 0.188 6 (-1, -1) 600 (-9, -7) -1.933 0 (-8.754 9, -6.821 9) (-9, -7) 0.067 0 (-9, -7) 800 (1, 1) 0.876 7 (2.925 6, 2.0489) (3, 2) 0.879 2 (1, 1) 在周跳修复过程中,以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范数最小作为周跳修复成功的判断准则。与直接求解周跳解并取整相比,该算法不仅能够验证周跳解的正确性,还可对错误的周跳解进行实时纠正,得到准确的周跳值,对历元后续观测值进行修复。
An Improved Algorithm for Real-Time Cycle Slip Detection and Repair Based on TurboEdit Epoch Difference Model
-
摘要: 为了解决电离层扰动和低高度角下非差数据的实时周跳探测与修复问题,提出了一种基于TurboEdit历元差模型的优化算法。首先在无几何距离组合历元差中加入滑动多项式拟合法提取电离层延迟与多路径效应误差的趋势项,通过选择合适的阈值,增强对小周跳的探测能力;然后针对历元差模型无法区别周跳与粗差的问题,引入统计诊断中的Score检验量并构建了基于Score检验的周跳与粗差分离模型,避免了粗差引起的周跳误探现象;最后在成功探测出周跳后,使用探测量最小1范数准则来选取正确的周跳修复值。采用国际GNSS服务站GPS非差数据对所提算法进行验证,实验结果表明,该优化算法能够有效抑制低高度角时的多路径效应和电离层延迟误差,实时剔除粗差并准确修复周跳,对小周跳具有较高的探测灵敏度与较低的误判率。
-
关键词:
- 实时周跳探测与修复 /
- TurboEdit历元差模型 /
- 滑动多项式拟合法 /
- Score检验模型 /
- 最小1范数
Abstract:Objectives Carrier phase is the key to achieve high precision positioning, and its quality determines the reliability of positioning results. Due to the development and needs of society, the data acquisition environment is becoming more and more complex, and the real-time requirements are also constantly improving, which greatly increases the probability of cycle slip and increases the difficulty of cycle slip detection. However, some cycle slip detection algorithms commonly used at present have poor detection effect in such environments as active ionosphere, low elevation angle, occlusion, etc., which lead to false detection and missed detection. Methods In order to solve these problems, firstly, an improved algorithm based on TurboEdit epoch difference model is proposed. The trend terms of ionospheric delay and multipath effect errors are extracted by adding sliding polynomial fitting method to the epoch difference of GF (Geometry-Free), and by choosing appropriate threshold, the detection ability of small cycle slips is enhanced. Then, to solve the problem that epoch difference model can not distinguish cycle slip and gross errors, Score test in statistical diagnosis is introduced and constructed. A cycle slip and gross errors separation model based on Score test is built to avoid the phenomenon of cycle slip detection caused by gross errors. Finally, the improved GF epoch difference model is combined with Melbourne Wübbena(MW) combination to repair cycle slip, and the minimum one norm criterion of detection is used to select the correct repair value of cycle slip. The algorithm is verified by using the GPS un-differenced data of IGS (International GNSS Service) stations. Results The cycle slips obtained by the improved GF epoch difference still have high detection sensitivity for small cycle slips. The simulated small cycle slip, continuous cycle slip and special cycle slip are detected, such as (1, 1), (-1, -1), (5, 4), etc. Two successive abnormal epochs detected by GF epoch difference optimization algorithm have successfully separated the simulated (1, 1), (0, 0.5) small gross errors after a Score test and deviation test.When GF and MW probes are used to calculate cycle slip directly, there are errors in cycle slip calculation due to the excessive noise of MW probe, and the calculation error reaches 3 weeks. By using the improved minimum one norm of GF epoch differential survey and MW survey, the correct cycle slip solution is found successfully. Conclusions In the aspect of cycle slip detection, the improved GF epoch difference method can weaken the influence of ionospheric delay variation and multipath effect, and improve the accuracy of cycle slip detection. At the same time, through the appropriate threshold selection, even in the case of high ionospheric scintillation, small cycle slip and continuous cycle slip can be detected accurately. In the separation of cycle slip and gross errors, the fitting residual is replaced by the prediction residual to improve the Score test, which makes it suitable for anomaly test in real-time situation. Based on this, the separation model of cycle slip and gross errors can separate the gross errors from cycle slip accurately, avoid the misjudgment of cycle slip, and improve the reliability of real-time TurboEdit cycle slip detection method. In the aspect of cycle slip repair, the criterion of cycle slip repair success is to minimize the one norm of the probe change caused by correct cycle slip repair. Compared with solving cycle slip solution directly and rounding, this proposed algorithm can not only verify the correctness of cycle slip solution, but also correct the wrong cycle slip solution in real-time, get the accurate cycle slip value, and repair the subsequent observation value of epoch. -
表 1 阈值系数$ {\mathit{\boldsymbol{n}}} $选取结果
Table 1. Selection Results of Threshold Coefficient n
表 2 周跳与粗差模拟情况/周
Table 2. Situation of Cycle Slip and Gross Errors Simulation/cycle
历元 周跳 粗差 50 (1, 1) 100 (5, 4) 200 (1, 1) 300 (1, 0) 400 (-1, -1) 401 (-1, -1) 500 (1, 1) 600 (-9, -7) 700 (0, 0.5) 800 (1, 1) 表 3 周跳与粗差分离结果
Table 3. Separation Results of Cycle Slip andGross Errors
历元 Score值 未通过
个数是否进行预测值偏离检验 未通过
个数检验
结果200 16 424.315 5 2 是 1 粗差 201 170.242 9 400 5 035.961 8 2 是 2 连续
周跳401 69 444.983 8 500 275.990 6 1 否 1 粗差 501 2.387 8 700 4 210.110 3 2 是 1 粗差 701 435.842 1 表 4 周跳修复结果
Table 4. Results of Cycle Slip Repair
历元 周跳组合 $ \mathrm{\Delta }\mathrm{M}\mathrm{W} $ 周跳浮点解 取整结果 最小1范数 估计周跳 50 (1, 1) 0.270 2 (2.209 3, 1.939 1) (2, 2) 0.271 0 (1, 1) 100 (5, 4) 1.585 4 (7.776 5, 6.191 1) (8, 6) 0.897 8 (5, 4) 300 (1, 0) 0.821 2 (0.342 1, -0.479 1) (0, 0) 0.180 4 (1, 0) 400 (-1, -1) 0.031 0 (-0.837 3, -0.868 3) (-1, -1) 0.032 2 (-1, -1) 401 (-1, -1) 0.186 9 (-0.086 4, -0.273 3) (0, 0) 0.188 6 (-1, -1) 600 (-9, -7) -1.933 0 (-8.754 9, -6.821 9) (-9, -7) 0.067 0 (-9, -7) 800 (1, 1) 0.876 7 (2.925 6, 2.0489) (3, 2) 0.879 2 (1, 1) -
[1] 叶世榕. GPS非差相位精密单点定位理论与实现[D]. 武汉: 武汉大学, 2002 Ye Shirong.Theory and Its Realization of GPS Precise Point Position Undifferenced Phase Observation[D].Wuhan: Wuhan University, 2002 [2] Lichtenegger H, Hofmann-Wellenhof B. GPS-Data Preprocessing for Cycle-Slip Detection[M]// Global Positioning System: An Overview.New York: Springer, 1990 [3] Kee C, Walter T, Enge P, et al. Quality Control Algorithms on WAAS Wide-Area Reference Stations[J]. Navigation, 1997, 44(1): 53-62 doi: 10.1002/j.2161-4296.1997.tb01939.x [4] Blewitt G.An Automatic Editing Algorithm for GPS Data[J]. Geophysical Research Letters, 1990, 17(3): 199-202 doi: 10.1029/GL017i003p00199 [5] Liu Z. A New Automated Cycle Slip Detection and Repair Method for a Single Dual-Frequency GPS Receiver[J]. Journal of Geodesy, 2011, 85(3): 171-183 doi: 10.1007/s00190-010-0426-y [6] Cai C, Liu Z, Xia P, et al. Cycle Slip Detection and Repair for Undifferenced GPS Observations Under High Ionospheric Activity[J]. GPS Solutions, 2013, 17(2): 247-260 doi: 10.1007/s10291-012-0275-7 [7] Ju B, Gu D, Chang X, et al. Enhanced Cycle Slip Detection Method for Dual-Frequency BeiDou GEO Carrier Phase Observations[J]. GPS Solutions, 2017, 21(3): 1 227-1 238 doi: 10.1007/s10291-017-0607-8 [8] Zhang X H, Li X X.Instantaneous Re-Initialization in Real-Time Kinematic PPP with Cycle Slip Fixing[J]. GPS Solutions, 2012, 16(3): 315-327 doi: 10.1007/s10291-011-0233-9 [9] 张小红, 曾琪, 何俊, 等. 构建阈值模型改善TurboEdit实时周跳探测[J]. 武汉大学学报·信息科学版, 2017, 42(3): 285-292 doi: 10.13203/j.whugis20150045 Zhang Xiaohong, Zeng Qi, He Jun, et al. Improving TurboEdit Real-Time Cycle Slip Detection by the Construction of Threshold Model[J]. Geomatics and Information Science of Wuhan University, 2017, 42(3): 285-292 doi: 10.13203/j.whugis20150045 [10] Banville S, Langley R B. Mitigating the Impact of Ionospheric Cycle Slips in GNSS Observations[J]. Journal of Geodesy, 2013, 87(2): 179-193 doi: 10.1007/s00190-012-0604-1 [11] Baarda W. A Testing Procedure for Use in Geodetic Networks[M]. Delft: Netherlands Geodetic Commission, 1968 [12] 黄令勇, 翟国君, 欧阳永忠, 等. 三频GNSS电离层周跳处理[J]. 测绘学报, 2015, 44(7): 717-725 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201507003.htm Huang Lingyong, Zhai Guojun, Ouyang Yongzhong, et al. Ionospheric Cycle Slip Processing in Triple-Frequency GNSS[J]. Acta Geodaetica et Cartographica Sinica, 2015, 44(7): 717-725 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201507003.htm [13] 潘雄, 吕玉婷, 汪耀, 等. 基于半参数平差模型的粗差定位与定值研究[J]. 武汉大学学报·信息科学版, 2016, 41(11): 1 421-1 426 doi: 10.13203/j.whugis20140493 Pan Xiong, Lü Yuting, Wang Yao, et al. Research of the Location and Valuation of Gross Error Based on Semi-parametric Adjustment Model[J]. Geomatics and Information Science of Wuhan University, 2016, 41(11): 1 421-1 426 doi: 10.13203/j.whugis20140493 [14] Boos D. On Generalized Score Tests[J]. American Statistician, 1992, 46(4): 327-333 http://search.ebscohost.com/login.aspx?direct=true&db=aph&AN=9605216176&site=ehost-live [15] Weisberg C S. Diagnostics for Heteroscedasticity in Regression[J]. Biometrika, 1983, 70(1): 1-10 doi: 10.1093/biomet/70.1.1 [16] 韦博成, 林金官, 解锋昌. 统计诊断[M]. 北京: 高等教育出版社, 2009 Wei Bocheng, Lin Jinguan, Xie Fengchang. Statistical Diagnosis[M]. Beijing: Higher Education Press, 2009 [17] Lacy M C D. Real-Time Cycle Slip Detection in Triple-Frequency GNSS[J]. GPS Solutions, 2012, 16(3): 353-362 doi: 10.1007/s10291-011-0237-5 -