[学习笔记]周跳探测

GNSS观测中的周跳探测

1. 周跳定义:
  • 整周跳变是载波相位观测值特有的问题
  • 周跳:引起整周计数暂时中断,发生整周跳变。
  • 主要原因:卫星信号被某障碍物阻挡而无法到达接收机,或者由于外界干扰或接收机所处的动态条件恶劣而引起卫星信号的暂时失锁等。
2. 周跳探测与修复方法
2.1 基于观测值随时间变化规律的方法:
  • 高次差法、多项式拟合法
  • 载波相位观测值随时间变化主要受站星距离影响,而后者取决于卫星与接收机的运动状态。相比之下,接收机运动状态较为稳定,故此法适用于静态数据处理,需考虑卫星钟差、接收机钟差、对流层折射及电离层折射随时间的变化。
2.2 基于不同观测值组合的方法:
  • 单频/双频码相组合法、电离层残差法、多普勒积分法
  • 利用不同观测值之间的关系来探测周跳,组合后的观测值一般与站星几何距离无关,该法不受卫星钟差、接收机钟差及接收机运动状态影响
2.3 基于观测值估值残差的方法
  • 此类方法根据参数估计后或得到的观测值的估计残差来确定周跳
3. 常用的周跳探测方法
3.1 双频码相组合法
  • 指利用双频载波相位和测码伪距组合观测值来探测/修复周跳

  • 常用的为Melbourne-Wubbena组合(M-W组合)

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

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

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

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

φ 2 = ρ λ 2 + A c f 2 2 − N 2 (4) \varphi_2=\frac{\rho}{\lambda_2}+\frac{A}{c{f_2}^2}-N_2\tag{4} φ2=λ2ρ+cf22AN2(4)

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

(1)-(2)得:
P 1 − P 2 = A f 1 2 − f 2 2 f 1 2 f 2 2 (5) P_1-P_2=A\frac{{f_1}^2-{f_2}^2}{{f_1}^2{f_2}^2}\tag{5} P1P2=Af12f22f12f22(5)

A = f 1 2 f 2 2 f 1 2 − f 2 2 ( P 1 − P 2 ) (6) A=\frac{{f_1}^2{f_2}^2}{{f_1}^2-{f_2}^2}(P_1-P_2)\tag{6} A=f12f22f12f22(P1P2)(6)

(3)-(4)得:
φ 1 − φ 2 − ρ ( 1 λ 1 − 1 λ 2 ) − A c ( 1 f 1 − 1 f 2 ) = N 2 − N 1 (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_2-N_1\tag{7} φ1φ2ρ(λ11λ21)cA(f11f21)=N2N1(7)
根据无电离组合观测值有:
ρ = f 1 2 f 1 2 − f 2 2 P 1 − f 2 2 f 1 2 − f 2 2 P 2 (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} ρ=f12f22f12P1f12f22f22P2(8)
把(6)(8)带入(7)中得到:
φ 1 − φ 2 − f 1 − f 2 f 1 + f 2 ∗ ( P 1 λ 1 + P 2 λ 2 ) = N 2 − N 1 (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φ2f1+f2f1f2(λ1P1+λ2P2)=N2N1(9)
其中, φ 1 − φ 2 \varphi_1-\varphi_2 φ1φ2即为宽巷观测值 φ Δ \varphi_\Delta φΔ,将 φ Δ = c f Δ = c f 1 − f 2 \varphi_\Delta=\frac{c}{f_\Delta}=\frac{c}{f_1-f_2} φΔ=fΔc=f1f2c带入(9)中可得到:
φ Δ λ Δ − f 1 P 1 + f 2 P 2 f 1 + f 2 + 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组合也可以写成:
M W = f 1 − f 2 f 1 + f 2 ( P 1 λ 1 + P 2 λ 2 ) − ( φ 1 − φ 2 ) = N 1 − N 2 MW=\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+f2f1f2(λ1P1+λ2P2)(φ1φ2)=N1N2
M-W组合消除了电离层、对流层、钟差和几何距离的影响,等式右边仅剩下宽巷模糊度,未发生周跳时,M-W检验量应在一常数附近波动,当有周跳发生时,检验量会发生突变。因此,逐历元解算各卫星的MW组合观测值,并通过历元间差分获得周跳探测检验量 Δ M W \Delta MW ΔMW,若 Δ M W \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观测值 M W = c ∗ ( λ 1 − λ 2 ) f 1 − f 2 − f 1 P 1 + f 2 P 2 f 1 + f 2 MW=\frac{c*(\lambda_1-\lambda_2)}{f_1-f_2}-\frac{f_1P_1+f_2P_2}{f_1+f_2} MW=f1f2c(λ1λ2)f1+f2f1P1+f2P2,当 Δ M W \Delta MW ΔMW大于阈值THRES_MW_JUMP=10.0时,在slip中标记周跳

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

G F = λ 1 φ 1 − λ 2 φ 2 = − λ 1 N 1 + λ 2 N 2 + A f 1 2 − A f 2 2 = − λ 1 N 1 + λ 2 N 2 − ( P 1 − P 2 ) 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+f12Af22A=λ1N1+λ2N2(P1P2)

L I = φ 1 λ 1 − φ 2 λ 2 L_I=\varphi_1\lambda_1-\varphi_2\lambda_2 LI=φ1λ1φ2λ2,即为电离层残差组合,有:
L I + ( P 1 − P 2 ) − N 1 λ 1 + N 2 λ 2 = 0 L_I+(P_1-P_2)-N_1\lambda_1+N_2\lambda_2=0 LI+(P1P2)N1λ1+N2λ2=0
上式消除了电离层延迟、卫星至接收机几何距离、卫星钟差和接收机钟差的影响,也可用于确定 N 1 N_1 N1 N 2 N_2 N2

逐历元计算各卫星的GF组合观测值,并通过历元间求差或的周跳探测检验量 Δ G F \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观测值: G F = λ 1 φ 1 − λ 2 φ 2 = ( φ 1 f 1 − φ 2 f 2 ) / c GF=\lambda_1\varphi_1-\lambda_2\varphi_2=(\frac{\varphi_1}{f_1}-\frac{\varphi_2}{f_2})/c GF=λ1φ1λ2φ2=(f1φ1f2φ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组合的方法探测周跳,同样顾及采样率和卫星高度角,给定周跳探测的阈值,其公式为:
    请添加图片描述
    函数对应部分为:
    请添加图片描述

### 关于探测的概念与实现 #### 的定义及其影响 是指在连续观测过程中,由于某些外部因素(如卫星信号被遮挡、外界干扰或环境恶劣引起接收机失锁),导致整计数发生错误的现象[^4]。这种现象会使相位观测值比正常值差几个整数期的跃,尽管不足一的部分仍能继续工作,但整体上会显著降低GPS定位的精度[^1]。 #### 探测的重要性 的存在直接影响到GPS定位解算中的整模糊度求解过程,从而进一步削弱最终定位结果的可靠性。因此,在实际应用中,尤其是涉及高精度定位的任务时,必须对进行有效的探测和修复。 #### 常见的探测方法 以下是几种常用的探测技术: ##### 1. **双频组合法** 通过利用两个不同频率的载波相位观测值构建无电离层组合或其他线性组合来增强噪声抑制能力,并检测可能存在的异常变化。这种方法基于假设同一颗卫星在同一时刻不会同时出现多个独立的随机误差源,故任何突兀的变化都可视为潜在的事件[^2]。 ##### 2. **残差分析法** 此方法依赖于建立一个初步估计模型并计算相应的残差序列。如果某个特定时间段内的残差突然增大,则可以推测在此期间可能发生了一次或多起。具体操作包括但不限于最小二乘平差后的单位权标准偏差统计检验等手段[^3]。 ##### 3. **时间序列分析** 通过对一段时间范围内的载波相位数据做滑动窗口平均或者指数加权移动平均处理后观察是否存在明显的不连贯点作为判断依据之一。此外还可以借助傅里叶变换之类工具寻找高频成分所对应的间断位置来进行辅助确认。 ##### 4. **几何无关组合 (Geometry-Free Combination)** 这是一种专门针对电离层延迟效应设计的方法,它能够有效消除大部分系统性的偏移项,从而使剩余部分更加敏感地反映出真实的物理状态改变情况——即是否有新的整数目加入到了当前累积总和之中去。 ```python def detect_cycle_slip(phase_data, threshold=0.5): """ Detect cycle slips using simple difference method. Args: phase_data (list): List of float values representing carrier phase measurements over time. threshold (float): Threshold value to determine a potential slip. Returns: list: Indices where possible cycle slips are detected. """ differences = [abs(b-a) for a, b in zip(phase_data[:-1], phase_data[1:])] return [i+1 for i, diff in enumerate(differences) if diff > threshold] # Example usage with dummy data dummy_phase_measurements = [1.0, 1.1, 1.2, 2.3, 2.4, 2.5] cycle_slipts_indices = detect_cycle_slip(dummy_phase_measurements) print(f"Cycle slips detected at indices {cycle_slipts_indices}") ``` 上述代码片段展示了一个简单的差异阈值算法用于识别可能发生的索引位置。 --- ###
评论 6
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值