[学习笔记]周跳探测

GNSS观测中的周跳探测

文章目录

GNSS观测中的周跳探测
- 1. 周跳定义：
- 2. 周跳探测与修复方法
- - 2.1 基于观测值随时间变化规律的方法：
  - 2.2 基于不同观测值组合的方法：
  - 2.3 基于观测值估值残差的方法
- 3. 常用的周跳探测方法
- - 3.1 双频码相组合法
  - 3.2 双频电离层残差法
  - 3.3 GF组合和MW组合各自的探测盲区：
- 4. GAMP中探测周跳

1. 周跳定义：

整周跳变是载波相位观测值特有的问题
周跳：引起整周计数暂时中断，发生整周跳变。
主要原因：卫星信号被某障碍物阻挡而无法到达接收机，或者由于外界干扰或接收机所处的动态条件恶劣而引起卫星信号的暂时失锁等。

2. 周跳探测与修复方法

2.1 基于观测值随时间变化规律的方法：

高次差法、多项式拟合法
载波相位观测值随时间变化主要受站星距离影响，而后者取决于卫星与接收机的运动状态。相比之下，接收机运动状态较为稳定，故此法适用于静态数据处理，需考虑卫星钟差、接收机钟差、对流层折射及电离层折射随时间的变化。

2.2 基于不同观测值组合的方法：

单频/双频码相组合法、电离层残差法、多普勒积分法
利用不同观测值之间的关系来探测周跳，组合后的观测值一般与站星几何距离无关，该法不受卫星钟差、接收机钟差及接收机运动状态影响

2.3 基于观测值估值残差的方法

此类方法根据参数估计后或得到的观测值的估计残差来确定周跳

3. 常用的周跳探测方法

3.1 双频码相组合法

指利用双频载波相位和测码伪距组合观测值来探测/修复周跳
常用的为Melbourne-Wubbena组合（M-W组合）

伪距测量和载波相位观测方程可简化为：

P1=ρ−Af12(1)P_1=\rho-\frac{A}{{f_1}^2}\tag{1} P1=ρ−f12A(1)

P2=ρ−Af22(2)P_2=\rho-\frac{A}{{f_2}^2}\tag{2} P2=ρ−f22A(2)

φ1=ρλ1+Acf12−N1(3)\varphi_1=\frac{\rho}{\lambda_1}+\frac{A}{c{f_1}^2}-N_1\tag{3} φ1=λ1ρ+cf12A−N1(3)

φ1=ρλ1+Acf12−N1(4)\varphi_1=\frac{\rho}{\lambda_1}+\frac{A}{c{f_1}^2}-N_1\tag{4} φ1=λ1ρ+cf12A−N1(4)

其中，ρ\rhoρ为卫星至接收机的距离与所有和频率无关的偏差改正项之和。

(1)-(2)得：
P1−P2=Af12−f22f12f22(5)P_1-P_2=A\frac{{f_1}^2-{f_2}^2}{{f_1}^2{f_2}^2}\tag{5} P1−P2=Af12f22f12−f22(5)

A=f12f22f12−f22(P1−P2)(6)A=\frac{{f_1}^2{f_2}^2}{{f_1}^2-{f_2}^2}(P_1-P_2)\tag{6} A=f12−f22f12f22(P1−P2)(6)

(3)-(4)得：
φ1−φ2−ρ(1λ1−1λ2)−Ac(1f1−1f2)=N1−N2(7)\varphi_1-\varphi_2-\rho(\frac{1}{\lambda_1}-\frac{1}{\lambda_2})-\frac{A}{c}(\frac{1}{f_1}-\frac{1}{f_2})=N_1-N_2\tag{7} φ1−φ2−ρ(λ11−λ21)−cA(f11−f21)=N1−N2(7)
根据无电离组合观测值有：
ρ=f12f12−f22P1−f22f12−f22P2(8)\rho=\frac{{f_1}^2}{{f_1}^2-{f_2}^2}P_1-\frac{{f_2}^2}{{f_1}^2-{f_2}^2}P_2\tag{8} ρ=f12−f22f12P1−f12−f22f22P2(8)
把(6)(8)带入(7)中得到：
φ1−φ2−f1−f2f1+f2∗(P1λ1+P2λ2)=N2−N1(9)\varphi_1-\varphi_2-\frac{f_1-f_2}{f_1+f_2}*(\frac{P_1}{\lambda_1}+\frac{P_2}{\lambda_2})=N_2-N_1\tag{9} φ1−φ2−f1+f2f1−f2∗(λ1P1+λ2P2)=N2−N1(9)
其中，φ1−φ2\varphi_1-\varphi_2φ1−φ2即为宽巷观测值φΔ\varphi_\DeltaφΔ，将φΔ=cfΔ=cf1−f2\varphi_\Delta=\frac{c}{f_\Delta}=\frac{c}{f_1-f_2}φΔ=fΔc=f1−f2c带入(9)中可得到：
φΔλΔ−f1P1+f2P2f1+f2+NΔλΔ=0\varphi_\Delta\lambda_\Delta-\frac{f_1P_1+f_2P_2}{f_1+f_2}+N_\Delta\lambda_\Delta=0 φΔλΔ−f1+f2f1P1+f2P2+NΔλΔ=0
另外，M-W组合也可以写成：
MW=f1−f2f1+f2(P1λ1+P2λ2)−(φ1−φ2)=N1−N2MW=\frac{f_1-f_2}{f_1+f_2}(\frac{P_1}{\lambda_1}+\frac{P_2}{\lambda_2})-(\varphi_1-\varphi_2)=N_1-N_2 MW=f1+f2f1−f2(λ1P1+λ2P2)−(φ1−φ2)=N1−N2
M-W组合消除了电离层、对流层、钟差和几何距离的影响，等式右边仅剩下宽巷模糊度，未发生周跳时，M-W检验量应在一常数附近波动，当有周跳发生时，检验量会发生突变。因此，逐历元解算各卫星的MW组合观测值，并通过历元间差分获得周跳探测检验量ΔMW\Delta MWΔMW，若ΔMW\Delta MWΔMW大于设定的阈值则判定其发生周跳，阈值的设置主要取决于伪距观测值的噪声水平。

MW组合法的特点：
- 不受采样间隔影响，其周跳探测能力主要取决于伪距观测的噪声水平。当卫星高度角较低时，由于伪距观测噪声较大，M-W组合对周跳的灵敏度较低，难以准确探测出2周以内的小周跳；
- MW组合不受接收机和卫星的几何距离、电离层折射以及卫星和接收机钟差影响，因而适用于动态、非差观测值周跳探测；
- 与载波相位相比，伪距观测水平较高，但是由于采用的宽巷观测值波长较长，因而可以探测出小周跳；
- 无法独立的区分发生周跳的频率，且当两个频率上发生的周跳数值相近或相等时，该检验量失效；
RTKLIB中对应函数：

其中，mwmeas为计算MW组合观测值函数：

其中，return (obs->L[0]-obs->L[1])*CLIGHT/(freq1-freq2)-(freq1*obs->P[0]+freq2*obs->P[1])/(freq1+freq2);

对应计算MW观测值MW=c∗(λ1−λ2)f1−f2−f1P1+f2P2f1+f2MW=\frac{c*(\lambda_1-\lambda_2)}{f_1-f_2}-\frac{f_1P_1+f_2P_2}{f_1+f_2}MW=f1−f2c∗(λ1−λ2)−f1+f2f1P1+f2P2，当ΔMW\Delta MWΔMW大于阈值THRES_MW_JUMP=10.0时，在slip中标记周跳

3.2 双频电离层残差法

又称为无几何距离（Geometery Free, GF）组合法，指利用双频载波相位观测值的电离层残差来探测/修复周跳（基本思路是观察不同历元间的电离层残差的变化）
无几何距离组合：

GF=λ1φ1−λ2φ2=−λ1N1+λ2N2+Af12−Af22=−λ1N1+λ2N2−(P1−P2)GF=\lambda_1\varphi_1-\lambda_2\varphi_2\\ =-\lambda_1N_1+\lambda_2N_2+\frac{A}{{f_1}^2}-\frac{A}{{f_2}^2}\\ =-\lambda_1N_1+\lambda_2N_2-(P_1-P_2) GF=λ1φ1−λ2φ2=−λ1N1+λ2N2+f12A−f22A=−λ1N1+λ2N2−(P1−P2)

令LI=φ1λ1−φ2λ2L_I=\varphi_1\lambda_1-\varphi_2\lambda_2LI=φ1λ1−φ2λ2，即为电离层残差组合，有：
LI+(P1−P2)−N1λ1+N2λ2=0L_I+(P_1-P_2)-N_1\lambda_1+N_2\lambda_2=0 LI+(P1−P2)−N1λ1+N2λ2=0
上式消除了电离层延迟、卫星至接收机几何距离、卫星钟差和接收机钟差的影响，也可用于确定N1N_1N1和N2N_2N2。

逐历元计算各卫星的GF组合观测值，并通过历元间求差或的周跳探测检验量ΔGF\Delta GFΔGF，一般而言相邻历元间电离层残差非常小，任何异常变化都可以表明在一个或两个频率的相位观测中发生了周跳，阈值的设置主要取决于历元间电离层延迟的残差。（对于30s以内的观测值，阈值通常可以取0.05~0.1m）

GF组合特点：
- GF组合由于消除了几何距离的影响，同样适用于静态和动态观测数据，且明显受到观测值采样间隔和卫星高度角影响。未发生周跳时，相邻历元电离层残差变化非常小，仅为几毫米每秒。当卫星高度角较低时，受载波相位观测值噪声影响，对应电离层残差变化相对较大，但其数值仍小于0.05m；
- 该组合不受卫星和接收机钟差影响；
- 仅由载波相位观测值构成检验量，精度较高可检测小周跳；
- 该方法是基于电离层平静期的假设（相邻历元间电离层折射延迟较小），当电离层活跃或者电离层环境差异较大时（如低轨卫星数据），该方法将难以准确探测周跳。
RTKLIB中对应函数：

其中，gfmea函数为计算GF组合， return (obs->L[0]/freq1-obs->L[1]/freq2)*CLIGHT;

对应GF观测值：GF=λ1φ1−λ2φ2=(φ1f1−φ2f2)/cGF=\lambda_1\varphi_1-\lambda_2\varphi_2=(\frac{\varphi_1}{f_1}-\frac{\varphi_2}{f_2})/cGF=λ1φ1−λ2φ2=(f1φ1−f2φ2)/c

3.3 GF组合和MW组合各自的探测盲区：

摘自《GNSS精密单点定位理论方法及其应用》P104

4. GAMP中探测周跳

在gampPos.c->execses()函数中，根据全局变量中的采样间隔PPP_Glo.sample，通过calCsThres()函数按照下式分别计算GF和MW的阈值：

其中，calCsThres()函数按照采样间隔分别设置对应检测量阈值：

对应函数为：

在pppos()函数中，通过detecs()->detslp_mw()+detslp_gf()函数进行周跳探测；
GAMP中，detslp_mw()函数实现M-W组合进行周跳探测，考虑卫星高度角+之前根据采样间隔计算的阈值，两种因素来计算计算dmw的阈值

阈值计算公式为：

在detslp()函数中相应部分为：

阈值取MIN(thres*fact,6.0)，如果wl1-wl0大于阈值则认为发生周跳。
detslp_gf()函数实现利用GF组合的方法探测周跳，同样顾及采样率和卫星高度角，给定周跳探测的阈值，其公式为：

函数对应部分为：