关于RTKLIB中用多普勒检测周跳方法的改进和思考

熟悉RTKLIB中周跳检测方法的同学应该知道，RTKLIB中利用多普勒积分检测周跳的方法，由于存在接收机时间跳变的问题，因此该方法一直没有进行使用。前段时间看到rtklibexplore博主对多普勒积分方法进行了改进（具体参见博客Google Smartphone Decimeter Challenge），因此写博客记录分析下他的主要改进思想以及我的一些思考。

RTKLIB中的周跳检测方法

在之前撰写的关于周跳检测的博客 – RTKlib源码解析：ppp和rtkpost中的周跳检测函数中，介绍了RTKLIB在PPP和RTK中使用的周跳检测函数。

RTKLIB中的周跳检测主要包括了以下几种检测方法：

• 通过RINEX文件中标注的LLI进行检测。我觉得LLI如果标注得很好，理论上其实其他检测方法不用都行，因为LLI一般是通过锁定时间（lock time）来确定的，不过这取决于各个厂家数据转换软件的优劣程序。
• 通过MW方法来检测周跳。MW组合中包括了L1，L2的载波相位和伪距，从MW组合的推导中可以看出该组合消除了电离层延迟误差、卫星、接收机钟差、卫星至接收机的几何距离。该组合虽有较高的检测精度，但是如果L1，L2同时发生相同大小的周跳，则不能成功检测。MW组合的噪声与伪距噪声相关：${\sigma }_{MW}^2\approx \frac{1}{2}{\sigma }_p^2$
• 通过GF组合检测周跳。GF组合中包括了L1，L2的载波相位，组合噪声只与载波相位噪声相关：${\sigma }_{GF}=\sqrt{2}{\sigma }_L$

我感觉目前周跳检测中存在的短板有：

• 从目前在用的方法可以看出，MW组合和GF组合需要使用双频数据，如果只有单频数据，那么这两种方法都无法使用。
• 此外，MW组合的噪声和伪距噪声相关，我们都知道伪距噪声比载波的噪声要大得多，因此MW的阈值设定得比较松，比如RTKLIB现在默认是10.0m，但是GF的阈值可以设定得比较紧，比如现在设定的是0.05m。这就导致MW方法对一些比较小的周跳，可能检测效果没有那么好。对于多数接收机，如果发生失锁，重新锁定卫星，载波相位虽然发生了周跳，但是一般都会用伪距去初始化载波相位，所以即使周跳了，其实周跳数也不会太大，MW方法的检测效果可能不会太理想。不过GF方法，我觉得还是非常灵敏的。

多普勒积分检测周跳中存在的问题

理论上，在没有周跳的情况下，历元间多普勒的积分值就等于载波相位的变化值。因此可以对多普勒进行积分计算，并计算前后两个历元的载波相位变化量，将两者进行比较，如果差值超过阈值，则认为有周跳发生。多普勒积分检测周跳的优点在于，只需要单频数据就可以进行周跳检测。

但是部分接收机等钟差累积到一定程度后会进行调整，例如septentrio的接收机在钟差达到0.5ms时，会调整1ms，使得钟差成为-0.5ms ，此时伪距和载波相位会跳变0.001*velcity_of_light = 299792.458m，如下图所示【2】。显然，这会造成我们周跳检测出错，将时钟跳变误检测为周跳。

从参考【3】来看，大部分接收机调整时钟时都是用的整数ms。

对多普勒积分检测周跳方法的改进和思考

1. rtklibexplore的改进

作者rtklibexplore暂时还没有把这部分代码正式提交到github，不过可以通过这个链接RTKLIB demo5 b34e下载得到源码。

作者的主要改进思想为：

• 时钟跳变对当前历元所有的（多普勒积分 - 载波相位变化值）影响相同，让它们都增加了一个共同的偏差。
• 对于每颗卫星、每个频点，计算（多普勒积分 - 载波相位变化值），将不异常的（多普勒积分 - 载波相位变化值）加起来计算平均值，将这个平均值看作时钟跳变带来的偏差。目前的代码中对“不异常”的判定条件是（多普勒积分 - 载波相位变化值）没有超过3倍设定的阈值。然后从每个（多普勒积分 - 载波相位变化值）中减去这个平均值，再和阈值进行比较。

显然，这个改进方法也并不完美，如果时钟跳变很大，那么所有（多普勒积分 - 载波相位变化值）看起来都是异常值，这个周跳检测方法也就不起作用了。

在RTKLIB demo5 b34e中rtkpos.c中的detslp_dop函数如下所示，也就是上述改进方法的实现。

/* detect cycle slip by doppler and phase difference -------------------------*/
static void detslp_dop(rtk_t *rtk, const obsd_t *obs, const int *ix, int ns,
                       int rcv, const nav_t *nav)
{
    int i,ii,f,sat,ndop=0;
    double dph,dpt,lam,thres,mdop=0;
    double dopdif[MAXSAT][NFREQ], tt[MAXSAT][NFREQ];

    trace(4,"detslp_dop: rcv=%d\n", rcv);
    if (rtk->opt.thresdop<=0) return;  /* skip test if doppler thresh <= 0 */

    /* calculate doppler differences for all sats and freqs */
    for (i=0;i<ns;i++) {
        ii = ix[i];
        sat=obs[ii].sat;

        for (f=0;f<rtk->opt.nf;f++) {
            dopdif[i][f]=0;tt[i][f]=0.00;
            if (obs[ii].L[f]==0.0||obs[ii].D[f]==0.0||rtk->ssat[sat-1].ph[rcv-1][f]==0.0) continue;
            if (fabs(tt[i][f]=timediff(obs[ii].time,rtk->ssat[sat-1].pt[rcv-1][f]))<DTTOL) continue;
            if ((lam=sat2freq(obs[ii].sat,obs[ii].code[f],nav))<=0.0) continue;

            /* calc phase difference and doppler x time (cycle) */
            dph=obs[ii].L[f]-rtk->ssat[sat-1].ph[rcv-1][f];
            dpt=-obs[ii].D[f]*tt[i][f];
            dopdif[i][f]=dph-dpt;

            /* if not outlier, use this to calculate mean */
            if (fabs(dopdif[i][f])<3*rtk->opt.thresdop*fabs(tt[i][f])) {
                mdop+=dopdif[i][f]/fabs(tt[i][f]);
                ndop++;
            }
        }
    }
    /* calc mean doppler diff, most likely due to clock error */
    mdop=ndop>0?mdop/ndop:0;

    /* set slip if doppler difference with mean removed exceeds threshold */
    for (i=0;i<ns;i++) {
        sat=obs[ix[i]].sat;

        for (f=0;f<rtk->opt.nf;f++) {
            if (dopdif[i][f]==0.00) continue;
            thres=rtk->opt.thresdop*fabs(tt[i][f]);
            errmsg(rtk,"slip detected doppler (sat=%2d rcv=%d dL%d=%.3f thres=%.3f off=%.3f)\n",
                   sat,rcv,f+1,dopdif[i][f]-mdop,thres,mdop);
            if (fabs(dopdif[i][f]-mdop)>thres) {
                rtk->ssat[sat-1].slip[f]|=1;
            }
        }
    }
}

2. 其他博主改进

另外，我也看到博主PhD.Knight，在CSDN发表的博客【干货】提升RTK模糊度固定率的建议之周跳探测中，对RTKLIB的doppler检测周跳方法进行了改进。

在原本的detslp_dop函数上增加以下几行代码，我感觉这种方法针对整ms时钟跳变应该是非常有效的。

int cj;
cj=floor(fabs(dph*lam/CLIGHT)*1000+0.5); // clock jump, integer ms
if(cj!=0) continue;

参考文献

[1] Google Smartphone Decimeter Challenge
[2] asterx-m2a_firmware_v4.8.0_reference_guide
[3] gnss receiver clocks
[4]【干货】提升RTK模糊度固定率的建议之周跳探测