RTKlib源码解析：ppp和rtkpost中的周跳检测函数

本文解析了RTKlib ppp.c中两个周跳检测函数detslp_mw和detslp_gf，以及rtkpos.c中的detslp_ll和detslp_dop函数。

由于PPP中的detslp_ll直接根据LLI进行周跳判断，比较简单，根据Rinex中对LLI的定义即可明白，因此不进行解析。

rtkpos.c中的周跳检测函数包括detslp_ll、detslp_gf_L1L2、detslp_gf_L1L5、detslp_dop，其中detslp_gf_L1L2、detslp_gf_L1L5与PPP中的detslp_gf相似，只不过这里使用的是双差后的载波相位。由于rtkpost中的detslp_ll比PPP中稍微复杂一些，因此对detslp_ll和detslp_dop进行了解析。

我所基于的代码版本是RTKlib 2.4.3的一个拓展版本RTKexplore Demo5，这个版本主要针对低成本的GNSS定位。该版本整体算法并未做较大更改，只是针对低成本接收机进行了完善。

文章目录

detslp_mw
detslp_gf
detslp_ll
detslp_dop

detslp_mw

static void detslp_mw(rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav)

所在文件：ppp.c
功能说明：通过MW组合检测周跳
参数说明：

* args   : rtk_t  *rtk       IO  gps solution structure
*          obsd_t *obs       I   satellite observations
*          int     n         I   number of the user(rover) observations at current epoch
*          nav_t  *nav       I   satellite navigation data
/* detect slip by Melbourne-Wubbena linear combination jump ------------------*/
static void detslp_mw(rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav)
{double w0,w1;int i,j;trace(4,"detslp_mw: n=%d\n",n);for (i=0;i<n&&i<MAXOBS;i++) {if ((w1=mwmeas(obs+i,nav))==0.0) continue;w0=rtk->ssat[obs[i].sat-1].mw;rtk->ssat[obs[i].sat-1].mw=w1;trace(4,"detslip_mw: sat=%2d mw0=%8.3f mw1=%8.3f\n",obs[i].sat,w0,w1);if (w0!=0.0&&fabs(w1-w0)>THRES_MW_JUMP) {trace(3,"detslip_mw: slip detected sat=%2d mw=%8.3f->%8.3f\n",obs[i].sat,w0,w1);for (j=0;j<rtk->opt.nf;j++) rtk->ssat[obs[i].sat-1].slip[j]|=1;}}
}

处理过程：

先调用mwmeas函数计算Melbourne-Wubbena组合值MW=λ1λ2λ2−λ1(ϕ1−ϕ2)−λ2P1+λ1P2λ2+λ1MW = \frac{\lambda_1\lambda_2}{\lambda_2 - \lambda_1} (\phi_1 - \phi_2)-\frac{\lambda_2 P_1+\lambda_1P_2}{\lambda_2 + \lambda_1}MW=λ2−λ1λ1λ2(ϕ1−ϕ2)−λ2+λ1λ2P1+λ1P2，其中ϕ1\phi_1ϕ1、ϕ2\phi_2ϕ2分别为L1，L2的载波相位测量值（单位:周），P1P_1P1、P2P_2P2为伪距测量值(单位:m);

将当前历元的MW组合值与上一次的组合值进行比较，如果差值的绝对值大于了阈值THRES_MW_JUMP（默认值10.0），则认为有周跳。

注意事项：

MWMWMW组合本质上是计算λwNw\lambda_wN_wλwNw，其中λw=cf1−f2\lambda_w=\frac{c}{f1-f2}λw=f1−f2c为宽巷组合的波长，Nw=N1−N2N_w = N_1-N_2Nw=N1−N2是宽巷组合的整周模糊度。可想而知，如果没有周跳发生，那么上一历元和当前历元的λwNw\lambda_wN_wλwNw应该是一致的。

MWMWMW组合的推导可以参考《GPS测量与数据处理》书中4.3.3节“不同类型观测值线性组合”一节。从MW组合的推导中可以看出该组合消除了电离层延迟误差、卫星、接收机钟差、卫星至接收机的几何距离。该组合虽有较高的检测精度，但是如果L1，L2同时发生相同大小的周跳，则不能成功检测。

MWMWMW组合的噪声水平可根据其表达式进行推导，结果如下：
σMW=[(f1f1−f2)2+(f2f1−f2)2]σL2+[(f1f1+f2)2+(f2f1+f2)2]σp2\sigma_{MW} = \sqrt{[(\frac{f_1}{f_1-f_2})^2+(\frac{f_2}{f_1-f_2})^2]\sigma_L^2+[(\frac{f_1}{f_1+f_2})^2+(\frac{f_2}{f_1+f_2})^2]\sigma_p^2}σMW=[(f1−f2f1)2+(f1−f2f2)2]σL2+[(f1+f2f1)2+(f1+f2f2)2]σp2
由于载波相位的噪声水平远低于伪距，因此：
σMW≈[(f1f1+f2)2+(f2f1+f2)2]σp2\sigma_{MW}\approx\sqrt{[(\frac{f_1}{f_1+f_2})^2+(\frac{f_2}{f_1+f_2})^2]\sigma_p^2} σMW≈[(f1+f2f1)2+(f1+f2f2)2]σp2
上式进一步近似，则：σMW2≈12σp2\sigma_{MW}^2\approx\frac{1}{2}\sigma_p^2σMW2≈21σp2

MWMWMW组合周跳检测的阈值设定，可以根据3中的噪声水平，结合接收机的噪声大小，进行设定；也可采用滑动窗口计算平均值和标准差的方法。RTKlib中使用固定值10，另外开源代码GAMP中也有根据采样间隔、高度角进行阈值设定的方法，感兴趣可以进一步参考其论文和开源代码。

detslp_gf

static void detslp_gf(rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav)

所在文件：ppp.c
功能说明：通过GF（几何无关）组合检测周跳
参数说明：

* args   : rtk_t  *rtk       IO  gps solution structure
*          obsd_t *obs       I   satellite observations
*          int     n         I   number of the user(rover) observations at current epoch
*          nav_t  *nav       I   satellite navigation data
/* detect cycle slip by geometry free phase jump -----------------------------*/
static void detslp_gf(rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav)
{`在这里插入代码片`double g0,g1;int i,j;trace(4,"detslp_gf: n=%d\n",n);for (i=0;i<n&&i<MAXOBS;i++) {if ((g1=gfmeas(obs+i,nav))==0.0) continue;g0=rtk->ssat[obs[i].sat-1].gf;rtk->ssat[obs[i].sat-1].gf=g1;trace(4,"detslip_gf: sat=%2d gf0=%8.3f gf1=%8.3f\n",obs[i].sat,g0,g1);if (g0!=0.0&&fabs(g1-g0)>rtk->opt.thresslip) {trace(3,"detslip_gf: slip detected sat=%2d gf=%8.3f->%8.3f\n",obs[i].sat,g0,g1);for (j=0;j<rtk->opt.nf;j++) rtk->ssat[obs[i].sat-1].slip[j]|=1;}}
}

处理过程：

先调用gfmeas函数计算GF组合值GF=λ1ϕ1−λ2ϕ2GF = \lambda_1\phi_1-\lambda_2\phi_2GF=λ1ϕ1−λ2ϕ2，其中ϕ1\phi_1ϕ1、ϕ2\phi_2ϕ2分别为L1，L2的载波相位测量值（单位:周），λ1\lambda_1λ1、λ2\lambda_2λ2为L1、L2的波长;

将当前历元的GF组合值与上一次的组合值进行比较，如果差值的绝对值大于了阈值（默认值0.05），则认为有周跳。

注意事项：

GF组合本质上是计算GF=λ1ϕ1−λ2ϕ2=λ1N1−λ2N2GF = \lambda_1\phi_1-\lambda_2\phi_2 =\lambda_1N_1-\lambda_2N_2GF=λ1ϕ1−λ2ϕ2=λ1N1−λ2N2，如果没有周跳发生，那么上一历元和当前历元的GF值应该是一致的。

GF组合根据其表达式可以推导其噪声：σGF=2σL\sigma_{GF}=\sqrt{2}\sigma_LσGF=2σL。该周跳检测的阈值也可以根据噪声水平进行设定，GAMP中同样也根据采样间隔和高度角对阈值进行了调整，感兴趣可以参考其论文和代码。

detslp_ll

static void detslp_ll(rtk_t *rtk, const obsd_t *obs, int i, int rcv)

所在文件：rtkpost.c
功能说明：通过LLI标志进行周跳检测
参数说明：

* args   : rtk_t  *rtk       IO  gps solution structure
*          obsd_t *obs       I   satellite observations
*          int     i         I   index of obs
*          int    rcv        I   1: rover receiver; 2: base receiver

处理过程：

由于前向处理和后向处理，在利用LLI进行周跳判断上有所不同。因此在处理过程中，仅以前向为例；

LLI=getbitu(&rtk->ssat[sat-1].slip[f],0,2); 首先将上一历元该卫星的周跳标志取出存在LLI变量中，该变量之后会用来判断上一历元和当前历元，半周跳标志（LLI&2，即LLI的bit 1位）是否发生变化；

slip=obs[i].LLI[f];将当前历元原始观测量中的LLI赋值给slip变量；

如果上一历元和当前历元的半周跳标志不同，则认为有周跳，将slip的第0位置1；

存放当前历元的LLI和周跳。

注意事项：

后向处理和前向处理之间的区别：仔细观察代码，会发现后向处理中，判断周跳利用的是前一历元的LLI，而不是当前历元的LLI，而前向处理，使用的是当前历元的LLI。原因：假设一组连续8个历元的整周模糊度为：1-1-1-1-50-50-50-50，在历元5时，LLI=1。那么在前向处理时，将第5个历元标记为周跳是没问题的，因为在历元5，整周模糊度由1跳变为50。但是在后向处理时，虽然原始数据中历元5的LLI为1，但是实际周跳应该在历元4，因为后向处理时，历元4的整周模糊度由50跳变为1。

我的疑惑：

rtkpost.c中detslp_ll比较重要的特点是考虑了半周跳变化的情况，如果半周跳标志（LLI=2）前后两个历元发生变化，则认为有周跳。如果一直维持半周跳情况，并没有将slip的第0位置1，在之后的处理中仍然会使用该数据，这样似乎不太合理？

detslp_dop

static void detslp_dop(rtk_t *rtk, const obsd_t *obs, int i, int rcv, const nav_t *nav)

所在文件：rtkpost.c
功能说明：通过多普勒值进行周跳检测
参数说明：

* args   : rtk_t  *rtk       IO  gps solution structure
*          obsd_t *obs       I   satellite observations
*          int     i         I   index of obs
*          int    rcv        I   1: rover receiver; 2: base receiver
*          nav_t  *nav       I   satellite navigation data

处理过程：

因为这个函数在RTKlib2.4.3中已经不再使用，所以有些变量的定义不太清楚具体的含义，但是大致处理方法和流程还是比较清晰的。

首先计算周跳检测阈值，thres=MAXACC*tt *tt/2.0/lam+rtk->opt.err[4]*fabs(tt)*4.0;

计算前后两个历元的载波相位差值：dph=obs[i].L[f]-rtk->ph[rcv-1][sat-1][f];计算多普勒的积分值，dpt=-obs[i].D[f]*tt;；

将上述两个计算值相减，如果差值的绝对值超过阈值，则认为有周跳。

注意事项：

RTKlib注释中提到，detslp_dop函数由于clock-jump的问题，所以暂时没有使用。我猜测clock-jump可能是历元间时间的跳变，因为理论上应该对多普勒进行连续积分，如果历元时间跳变，直接用dpt=-obs[i].D[f]*tt计算积分值，误差是很大的。

不考虑clock-jump的问题，利用前一历元和当前历元，这两个历元的平均多普勒进行梯形积分，可能会比仅使用当前历元的多普勒值进行矩形积分的效果更好。即多普勒积分值进行如下计算：
dD(k)=−[D(k)+D(k−1)]∗0.5∗dtd_D(k) = - [D(k) + D(k-1)] *0.5*dtdD(k)=−[D(k)+D(k−1)]∗0.5∗dt

RTKlib源码解析：ppp和rtkpost中的周跳检测函数相关推荐

RTKLIB源码解析（一）、单点定位(pntpos.c)
目录 pntpos satposs estpos raim_fde estvel ephclk satpos satsys seleph eph2clk ephpos eph2pos rescod ...
rtklib源码 rtk差分解算，rtkpos和replos函数流程梳理
rtklib源码 rtk差分解算,rtkpos和replos函数流程梳理 rtkpos函数梳理总体流程 replos函数梳理 replos总体流程 1. 通过satposs函数计算卫星的位置.速度等 ...
12、常见Conditional注解源码解析-ConditionalOnBean（写作中...）
ConditionalOnBean 源码解析 ConditionalOnBean 注解看看 ConditionalOnBean 注解 @Conditional(OnBeanCondition.cla ...
【backtrader源码解析7】backtrader中mathsupport中计算平均值、方差和标准差的函数的分析(含金量挺低的)
前面的几篇文章尝试通过优化backtrader几个时间处理函数来提高效率,使用cython能提高单个函数效率但是在整体中效率降低导致失败之后,我又尝试了使用numba来改进那几个时间处理函数,也失败了 ...
python3 socketserver源码解析_解读python中SocketServer源码
再看继承真正的大餐来之前,还是来点儿开胃菜!回顾一下关于类的继承的知识: 我们先看上面的代码,这是一个简单的类继承,我们可以看到父类Base和子类Son,它们中各有一个Testfunc方法,当我们实 ...
RTKLIB源码解析（二）、 RTK定位(rtkpos.c)
本博客是转载,感谢: rtklib代码详解--rtkpos.c - 博客园-哆啦A梦 - 博客园主函数:rtkpos 1. 设置基站位置 2. 统计基站和流动站的卫星数量 3. 单点定位解算 4. ...
SpringSession的源码解析（从Cookie中读取Sessionid,根据sessionid查询信息全流程分析）
前言上一篇我们介绍了SpringSession中Session的保存过程,今天我们接着来看看Session的读取过程.相对保存过程,读取过程相对比较简单. 本文想从源码的角度,详细介绍一下Sessi ...
RTKlib相对定位源码解析： udstate函数
最近阅读RTKlib开源代码,非常感谢"塔奇克敲代码"博主的博客(RTKLIB源码解析--单点定位),他将单点定位部分整理成函数小卡片,为我理解RTKlib提供了很大的帮助.我参照 ...
RTKlib相对定位源码解析：zdres函数
最近阅读RTKlib开源代码,非常感谢"塔奇克敲代码"博主的博客(RTKLIB源码解析--单点定位),他将单点定位部分整理成函数小卡片,为我理解RTKlib提供了很大的帮助.我参照 ...
RTKlib相对定位源码解析： ddres函数
最近阅读RTKlib开源代码,非常感谢"塔奇克敲代码"博主的博客(RTKLIB源码解析--单点定位),他将单点定位部分整理成函数小卡片,为我理解RTKlib提供了很大的帮助.我参照 ...