RTKlib相对定位源码解析 - 周跳检测函数 detslp_mw，detslp_gf

GNSS RTKlib 相对定位

在此输入正文
首先是ppp.c中的两个周跳检测函数（detslp_mw和detslp_gf），由于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相似，只不过这里使用的是双差后的载波相位。另外对detslp_ll和detslp_dop提出两点疑问：
1. detslp_dop函数由于clock-jump的问题，所以暂时没有使用。Clock-jump是指多普勒测量值中的时钟漂移吗？还是历元的跳变？
2. rtkpost.c中detslp_ll比较重要的特点是考虑了半周跳变化的情况，如果半周跳标志（LLI=2）前后两个历元发生变化，则认为有周跳。如果一直维持半周跳情况，并没有将slip的第0位置1，在之后的处理中仍然会使用该数据，这样似乎不太合理？
3. rtkpost.c中detslp_ll，如果是后向滤波，除了2中提到的半周跳处理，感觉原始数据中LLI似乎没有起到作用？

我所基于的代码版本是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} (\phi_1 - \phi_2)-\frac{\lambda_2 P_1+\lambda_1P_2}{\lambda_2 + \lambda_1}$，其中$\phi_1$、$\phi_2$分别为L1，L2的载波相位测量值（单位:周），$P_1$、$P_2$为伪距测量值(单位:m);
  2. 将当前历元的MW组合值与上一次的组合值进行比较，如果差值的绝对值大于了阈值THRES_MW_JUMP（默认值10.0），则认为有周跳。

• 注意事项：

  1. $MW$组合本质上是计算$\lambda_wN_w$，其中$\lambda_w=\frac{c}{f1-f2}$为宽巷组合的波长，$N_w = N_1-N_2$是宽巷组合的整周模糊度。可想而知，如果没有周跳发生，那么上一历元和当前历元的$\lambda_wN_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\phi_1-\lambda_2\phi_2$，其中$\phi_1$、$\phi_2$分别为L1，L2的载波相位测量值（单位:周），$\lambda_1$、$\lambda_2$为L1、L2的波长;
  2. 将当前历元的GF组合值与上一次的组合值进行比较，如果差值的绝对值大于了阈值（默认值0.05），则认为有周跳。

• 注意事项：

  1. GF组合本质上是计算$GF = \lambda_1\phi_1-\lambda_2\phi_2 =\lambda_1N_1-\lambda_2N_2$，如果没有周跳发生，那么上一历元和当前历元的GF值应该是一致的。
  2. GF组合根据其表达式可以推导其噪声：$\sigma_{GF}=\sqrt{2}\sigma_L$。该周跳检测的阈值也可以根据噪声水平进行设定，GAMP中同样也根据采样间隔和高度角对阈值进行了调整，感兴趣可以参考其论文和代码。