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

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;
        }
    }
}

• 处理过程：

1. 先调用mwmeas函数计算Melbourne-Wubbena组合值$MW=\frac{{\lambda }_1{\lambda }_2}{{\lambda }_2-{\lambda }_1}({\varphi }_1-{\varphi }_2)-\frac{{\lambda }_2{P}_1+{\lambda }_1{P}_2}{{\lambda }_2+{\lambda }_1}$，其中${\varphi }_1$、${\varphi }_2$分别为L1，L2的载波相位测量值（单位:周），$P_1$、$P_2$为伪距测量值(单位:m);
2. 将当前历元的MW组合值与上一次的组合值进行比较，如果差值的绝对值大于了阈值THRES_MW_JUMP（默认值10.0），则认为有周跳。

• 注意事项：

1. $MW$组合本质上是计算${\lambda }_w{N}_w$，其中${\lambda }_w=\frac{c}{f1-f2}$为宽巷组合的波长，$N_w=N_1-N_2$是宽巷组合的整周模糊度。可想而知，如果没有周跳发生，那么上一历元和当前历元的${\lambda }_w{N}_w$应该是一致的。
2. $MW$组合的推导可以参考《GPS测量与数据处理》书中4.3.3节“不同类型观测值线性组合”一节。从MW组合的推导中可以看出该组合消除了电离层延迟误差、卫星、接收机钟差、卫星至接收机的几何距离。该组合虽有较高的检测精度，但是如果L1，L2同时发生相同大小的周跳，则不能成功检测。
3. $MW$组合的噪声水平可根据其表达式进行推导，结果如下：
   $${\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}$$
   由于载波相位的噪声水平远低于伪距，因此：
   $${\sigma }_{MW}\approx \sqrt{[(\frac{f_1}{f_1+f_2}{)}^2+(\frac{f_2}{f_1+f_2}{)}^2]{\sigma }_p^2}$$
   上式进一步近似，则：${\sigma }_{MW}^2\approx \frac{1}{2}{\sigma }_p^2$
4. $MW$组合周跳检测的阈值设定，可以根据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;
        }
    }
}

• 处理过程：

1. 先调用gfmeas函数计算GF组合值$GF={\lambda }_1{\varphi }_1-{\lambda }_2{\varphi }_2$，其中${\varphi }_1$、${\varphi }_2$分别为L1，L2的载波相位测量值（单位:周），${\lambda }_1$、${\lambda }_2$为L1、L2的波长;
2. 将当前历元的GF组合值与上一次的组合值进行比较，如果差值的绝对值大于了阈值（默认值0.05），则认为有周跳。

• 注意事项：

1. GF组合本质上是计算$GF={\lambda }_1{\varphi }_1-{\lambda }_2{\varphi }_2={\lambda }_1{N}_1-{\lambda }_2{N}_2$，如果没有周跳发生，那么上一历元和当前历元的GF值应该是一致的。
2. GF组合根据其表达式可以推导其噪声：${\sigma }_{GF}=\sqrt{2}{\sigma }_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 

• 处理过程：

1. 由于前向处理和后向处理，在利用LLI进行周跳判断上有所不同。因此在处理过程中，仅以前向为例；
2. LLI=getbitu(&rtk->ssat[sat-1].slip[f],0,2); 首先将上一历元该卫星的周跳标志取出存在LLI变量中，该变量之后会用来判断上一历元和当前历元，半周跳标志（LLI&2，即LLI的bit 1位）是否发生变化；
3. slip=obs[i].LLI[f];将当前历元原始观测量中的LLI赋值给slip变量；
4. 如果上一历元和当前历元的半周跳标志不同，则认为有周跳，将slip的第0位置1；
5. 存放当前历元的LLI和周跳。

• 注意事项：

1. 后向处理和前向处理之间的区别：仔细观察代码，会发现后向处理中，判断周跳利用的是前一历元的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。

• 我的疑惑：

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

• 处理过程：

1. 因为这个函数在RTKlib2.4.3中已经不再使用，所以有些变量的定义不太清楚具体的含义，但是大致处理方法和流程还是比较清晰的。
2. 首先计算周跳检测阈值，thres=MAXACC*tt *tt/2.0/lam+rtk->opt.err[4]*fabs(tt)*4.0;
3. 计算前后两个历元的载波相位差值：dph=obs[i].L[f]-rtk->ph[rcv-1][sat-1][f];计算多普勒的积分值，dpt=-obs[i].D[f]*tt;；
4. 将上述两个计算值相减，如果差值的绝对值超过阈值，则认为有周跳。

• 注意事项：

1. RTKlib注释中提到，detslp_dop函数由于clock-jump的问题，所以暂时没有使用。我猜测clock-jump可能是历元间时间的跳变，因为理论上 应该对多普勒进行连续积分，如果历元时间跳变，直接用dpt=-obs[i].D[f]*tt计算积分值，误差是很大的。
2. 不考虑clock-jump的问题，利用前一历元和当前历元，这两个历元的平均多普勒进行梯形积分，可能会比仅使用当前历元的多普勒值进行矩形积分的效果更好。即多普勒积分值进行如下计算：
   $$d_D(k)=-[D(k)+D(k-1)]\ast 0.5\ast dt$$