一种基于BDS-3四频数据的周跳探测与修复方法

吕 震 王振杰 聂志喜 徐晓飞 张远帆

1 中国石油大学(华东)海洋与空间信息学院，青岛市长江西路66号, 266580

载波相位观测是GNSS高精度定位及应用的关键[1]。由于接收机自身原因或GNSS信号受物理遮蔽等因素的影响，载波相位观测可能发生周跳，因此当载波相位观测量参与GNSS高精度定位及相关应用时，需要进行有效的周跳探测和修复[2]。目前常用的GNSS双频观测值周跳探测方法有TurboEdit方法、卡尔曼滤波法、多普勒积分法等，其中TurboEdit方法具有探测精度高、程序容易实现等优点，应用最为广泛[3]。随着三频GNSS的发展，相关学者针对三频周跳探测与修复的特性及方法展开了大量研究[4-6]。

如今中国北斗三号全球卫星导航系统BDS-3和欧盟Galileo系统已经可以播发4个频率的信号[7-8]，四频周跳处理方法相较于GNSS双频和三频具有更广阔的应用前景[9]，而目前针对GNSS四频周跳探测与修复特性及方法的研究较少。限于篇幅，本文仅利用BDS-3提供的四频信号研究周跳探测与修复的新方法，选用B1C、B1I、B2a、B3I四个频率信号，根据电离层延迟最小和组合观测值噪声最小原则对组合系数进行优选，联合四频无几何相位组合和四频无几何无电离层组合两种方法进行周跳探测，然后基于最小二乘法估计周跳浮点解，最后利用MGEX提供的BDS-3四频数据验证本文算法的有效性和可靠性。

1 四频组合观测值系数选择

1.1 四频无几何相位组合观测值特性分析

GNSS原始伪距和载波相位观测值的观测方程可以写为[9]：

Pi=ρ+c(dtr-dts)+T+ηiI1+εPi

(1)

λiφi=ρ+c(dtr-dts)+T-ηiI1+λiNi+λiεφi

(2)

对式(2)进行组合[10]，可以得到四频无几何GF相位组合：

φGF=αλ1φ1+βλ2φ2+γλ3φ3+δλ4φ4=

-ηGFI1+NGF+εGF

(3)

为满足无几何条件[10]，提高周跳探测的灵敏度，令

(4)

对相邻历元的GF组合观测值进行历间差分可得：

ΔNGF=NGF(k)-NGF(k-1)=

ΔφGF-ηGFΔI1+ΔεGF

(5)

式中，ΔφGF为相邻历元差分后的GF相位组合，ΔNGF为GF相位组合周跳探测量，k和k-1为当前历元和前一历元。

由式(5)可知，ΔNGF会受到ηGFΔI1和ΔεGF的影响，因此在选择组合系数时应尽量选择ηGFΔI1和ΔεGF较小的组合。在高采样率条件下，ΔI1的值极小，当ηGF也极小时，ηGFΔI1可忽略不计。根据误差传播定律，ΔNGF的标准差可表示为式(6)，此处σφ取0.01周：

σΔNGF=

(6)

以4倍周跳探测量标准差(对应正态分布假设下99.9%的置信水平)作为周跳探测阈值，GF相位组合探测出周跳的条件为：

|ΔNGF|≥4σΔNGF

(7)

上述公式中的4个频点分别对应BDS-3的B1C(1 575.420 mHz)、B1I(1 561.098 mHz)、B2a(1 176.450 mHz)和B3I(1 268.520 mHz)。在[-10,10]范围内BDS-3较优的GF相位组合系数及10周以内不敏感周跳个数如表1所示，由表可见，对于任一GF相位组合均存在不敏感周跳，但不同组合的不敏感周跳是不同的。为此，从表1中选取3个GF相位组合同时进行周跳探测，以减少不敏感周跳的组合数并提高周跳探测的灵敏度。以式(4)作为组合系数选择条件，同时保证联合使用3个GF相位组合时不存在10周以内的不敏感周跳。通过筛选，BDS-3组合[1,-1,0,0]、[0,1,0,-1]和[0,0,-1,1]满足条件，因此本文选择上述3个GF相位组合作为设定条件下的最优组合。由于四频GF相位组合中的任意4个都是线性相关的，因此无论使用多少个GF相位组合，总存在不敏感周跳。为解决这一问题，本文选择四频无几何无电离层GIF组合联合四频GF相位组合进行周跳探测。

1.2 四频无几何无电离层组合系数选择

根据式(1)和式(2)可知，四频GIF组合可以表示为：

φGIF=iφ1+jφ2+kφ3+lφ4-

iN1+jN2+kN3+lN4+εGIF

(8)

式中，i、j、k、l为GIF组合载波相位组合观测值系数，a、b、c、d为GIF组合伪距组合观测值系数，λGIF=c/(if1+jf2+kf3+lf4)为组合波长，iN1+jN2+kN3+lN4=NGIF为组合模糊度，组合观测噪声为εGIF=iεφ1+jεφ2+kεφ3+lεφ4-(aεp1+bεp2+cεp3+dεp4)/λGIF。

表1 BDS-3四频无几何相位组合及不敏感周跳个数

为消除几何误差项和电离层误差一阶项，尽可能降低伪距噪声的影响[4]，需满足：

(9)

当发生周跳时，通过历元间差分可构造周跳探测方程：

ΔNGIF=NGIF(k)-NGIF(k-1)=ΔφGIF+ΔεGIF

(10)

式中，ΔφGIF为相邻历元差分后的GIF组合，ΔNGIF为GIF组合周跳探测量，k和k-1为当前历元和前一历元。

根据误差传播定律，ΔNGIF的标准差可表示为式(11)，此处σφ取0.01周，σP取0.1 m：

σΔNGIF=

(11)

同样以4倍周跳探测量标准差作为周跳探测阈值，GIF组合探测出周跳的条件为：

|ΔNGIF|≥4σΔNGIF

(12)

在GIF组合系数筛选过程中，首先要确定载波相位组合系数[i,j,k,l]，然后在满足式(9)前两式的条件下在[-1,1]范围内对伪距组合系数a、b、c、d进行搜索，最后将满足a2+b2+c2+d2=min的伪距组合系数a、b、c、d及对应的载波相位组合系数[i,j,k,l]作为该设定条件下的最优GIF组合。通过筛选，BDS-3四频较优的GIF组合如表2所示，由表可知，观测噪声σΔNGIF越小，周跳探测灵敏度越高，因此为了减小σΔNGIF，应选择λGIF更长的超宽巷或宽巷组合。以4σΔNGIF作为探测阈值时，周跳探测组合均可实现对1周小周跳的探测，为了表达简便，以载波相位组合系数[i,j,k,l]表示GIF组合。当BDS-3组合[1,-1,0,0]具有较长的波长时，观测噪声达到最小，且周跳探测阈值均不超过0.1周，因此选择组合[1,-1,0,0]作为BDS-3四频GIF组合的周跳探测组合。针对3个GF相位组合中存在不敏感周跳的问题，采用一个GIF组合构成4个线性无关的组合，可以保证无公共不敏感周跳组合，从而实现对各类周跳的探测。

表2 BDS-3四频无几何无电离层组合

2 四频数据的周跳探测与修复

本文采用3个四频GF相位组合和1个四频GIF组合构造4个线性无关的周跳探测组合，以4倍周跳探测量标准差作为周跳探测阈值。当满足周跳探测量大于探测阈值时，表明相应历元发生了周跳。通过对GF相位组合系数和GIF组合系数进行优选，得到3个GF相位组合和1个GIF组合，其中BDS-3组合系数对应[1,-1,0,0]、[0,1,0,-1]、[0,0,-1,1]和[1,-1,0,0]，周跳探测阈值分别为0.015 3、0.017 2、0.019 7和0.081 1。

当载波相位观测值中的周跳被成功探测后，需要对周跳进行修复。在每个历元下，通过联立3个GF相位组合观测值和1个GIF组合观测值得到组合观测值L，这里L表示为[LGF1,LGF2,LGF3,LGIF]T，其中LGF1、LGF2、LGF3为3个不同的GF相位组合观测值，LGIF为GIF组合观测值；然后在相邻历元间对L进行求差，可得历元间差分后的组合观测值ΔL。当ΔL探测出周跳后，根据式(13)即可得到单个频点上的待估周跳值ΔX：

AΔX+ε=ΔL

(13)

式中，ε为观测噪声，ΔX=[ΔN1,ΔN2,ΔN3,ΔN4]T为4个不同频点上的待估周跳浮点解，A为其系数矩阵，此处表示为：

(14)

图1 BDS-3四频周跳探测与修复算法流程Fig.1 Procedure of quad-frequency cycle-slip detection and repair algorithm for BDS-3

3 实验分析

选取MGEX提供的WUH2测站(30.53°N，114.36°E)2020-10-26的观测数据，采样间隔为1 s，卫星截止高度角为10°。其中，BDS-3系统选择C39和C40两颗卫星，每颗卫星选择连续500个没有周跳和粗差的历元作为本次实验的观测时间序列。

由图2可见，3个GF载波相位组合探测量的波动均保持在[-0.02,0.02]范围内，而GIF组合周跳探测量在[-0.1,0.1]范围内。由此可知，相较于GIF组合，GF相位组合的周跳探测能力更高，这是因为GF相位组合比GIF组合受到的噪声影响更小。

图2 BDS-3无周跳组合探测值Fig.2 Combination detection values without cycle-slip for BDS-3

由于原始测站非差观测值中无周跳，因此本文在不同历元时刻加入一般周跳和特殊周跳(一般周跳指4个组合均可以探测出的1～10周的周跳，特殊周跳指单个频率上发生1～2周的小周跳或者在不同频率发生相同跳变的小周跳)，以验证本文所选线性组合的周跳探测能力及修复效果。

3.1 一般周跳的探测与修复

分别在BDS-3 C39和C40卫星每隔100个历元加入不同的10周以内一般小周跳(3,1,2,1)、(6,7,4,3)、(1,4,3,2)、(5,7,4,7)，探测结果如图3所示。由图可见，4个组合均可同时检测到一般周跳，BDS-3具体周跳探测与修复结果列于表3。由表可知，所有一般周跳的ratio值均大于阈值，同时利用LAMBDA算法得到的周跳计算值与周跳模拟值均保持一致，从而验证了采用LAMBDA算法对一般周跳进行修复的正确性。

图3 BDS-3一般周跳组合探测值Fig.3 Normal cycle-slip combination detection values for BDS-3

表3 BDS-3一般周跳的周跳探测与修复结果

3.2 特殊周跳的探测与修复

为进一步验证本文方法的周跳探测与修复效果，分别在BDS-3 C39和C40卫星不同历元时刻加入特殊周跳，探测结果如图4所示。

图4 BDS-3特殊周跳组合探测值Fig.4 Particular cycle-slip combination detection values for BDS-3

首先考虑当4个频率信号中某一单独频点发生周跳的情况，分别在第100、200、300、400历元中加入周跳(1,0,0,0)、(0,1,0,0)、(0,0,1,0)、(0,0,0,1)。由图4(a)和4(b)可见，GF相位组合[1,-1,0,0]和GIF组合均可探测出周跳(1,0,0,0)；对于周跳(0,1,0,0)，仅有GF相位组合[0,0,1,-1]对其不敏感；采用GF相位组合[0,0,-1,1]可探测出周跳(0,0,1,0)和(0,0,0,1)。考虑其他特殊的不敏感周跳，在第100和200历元加入4个频点发生相同周跳的周跳组合(1,1,1,1)和(5,5,5,5)，在第300历元加入频点2和频点3发生相同周跳的周跳组合(0,2,2,0)，在第400历元加入频点1、频点3、频点4同时发生相同周跳的周跳组合(3,3,0,3)。由图4(c)和4(d)可见，对于周跳组合(1,1,1,1)，C39和C40卫星的GF相位组合[0,1,0,-1]的探测量均超过探测阈值；对于周跳(5,5,5,5)、(0,2,2,0)和(3,3,0,3)，C39和C40卫星的GF相位组合[0,1,0,-1]和[0,0,-1,1]均可将其探测出。

通过分析可知，特殊的不敏感周跳可能对部分探测组合不敏感，但4个组合联合后均能很好地探测出周跳，极大地减少了探测盲区。BDS-3四个组合特殊周跳的周跳探测及修复具体信息列于表4，由表可知，所有特殊周跳的ratio值均大于阈值，且计算值均与模拟值相同，表明采用LAMBDA方法同样可以对特殊周跳进行准确修复。

表4 BDS-3特殊周跳的周跳探测与修复结果

4 结 语

本文对BDS-3四频数据的周跳探测与修复方法进行研究，首先基于无几何条件下的观测噪声最小准则优选出3个GF相位组合，再基于无几何无电离层条件下的观测噪声最小准则优选出1个GIF组合，并联合4个组合观测值进行周跳探测。在探测出周跳后，将4个周跳探测组合进行方程组联立，采用最小二乘法确定周跳浮点解，基于LAMBDA降相关搜索法获得周跳整数解，并采用ratio检验进一步对周跳整数解进行验证。最后采用MGEX提供的BDS-3四频观测数据进行实验验证，结果表明，本文提出的方法可以实现对单站非差观测值中各类周跳的准确探测及快速高效修复。值得说明的是，本文提出的方法同样可以用于Galileo四频数据的周跳探测与修复，限于篇幅，本文不进行具体介绍。