利用SSR修正的RTKLIB PPP测试

欢迎关注个人公众号：导航员学习札记

一、前言

第一次看到PPP-RTK这个名词，还挺疑惑和好奇的。这几天工作之余，整理了一下PPP-RTK、PPP、RTK的区别，把RTKLIB中PPP-RTK（即SSR修正）相关的代码也都看了一下，感觉还是比较简单的。另外也用实际数据测试了一下RTKLIB PPP使用SSR修正的效果，如果你感兴趣的话，就请继续看下去吧！

如果有小伙伴们知道组合导航、测绘方向的宝藏博主或者公众号，希望可以留言推给我，让我共同学习！

二、PPP-RTK概念

PPP-RTK实质是利用CORS网，生成卫星轨道、卫星钟差、卫星伪距偏差、卫星相位偏差、电离层延迟、对流程延迟等改正数，通过卫星向地面播发，从而加快PPP收敛速度的定位技术。

武汉大学张小红教授2021年5月的报告《PPP-RTK技术面临的机遇和挑战》中，下面两张图比较清晰的展现了PPP、RTK、PPP-RTK的不同之处，以及优劣势。

如果在PPP中能得到所有的修正，其实和RTK也就差不多了，因为RTK本质上也是和基站双差后消除各种误差，使得模糊度具有整周特性。下图比较好的展示了在RTK、PPP-RTK、PPP中分别消除了哪些误差。

三、SSR修正概念

SSR（State Space Representation,状态空间域）修正其实就是一种修正方式，将观测值误差分解为卫星轨道、卫星钟差、卫星伪距偏差、卫星相位偏差、电离层延迟、对流程延迟等“状态量误差”；与之相对是RTK中的OSR（Observation Space Representation，观测值域），基站直接提供融合了各种修正的观测值。

根据参考[3]，RTCM标准委员会把SSR修正数据开发分成三个阶段：

• 第一阶段：包括卫星轨道、卫星时钟、伪距偏差修正数，以及URA
• 第二阶段：包括天顶TEC电离层修正，卫星相位偏差修正
• 第三阶段：包括电离层倾斜修正，对流程修正

SSR修正数据的修正方法在参考[3]或者其他文献中很容易就可以找到，这里不再赘述。第一阶段的RCTM SSR修正信息如下：

四、RTKLIB中SSR修正相关代码

RTKLIB目前主要支持第一阶段的SSR修正，即卫星轨道、卫星时钟、伪距偏差修正（也在试验性支持第二阶段中的相位偏差修正）。与SSR修正相关的源码主要包括一下几个部分：

1. readpreceph函数：位于postpos.c中，设置RTCM文件，以及初始化rtcm结构体，核心代码如下：

    /* set rtcm file and initialize rtcm struct */
    rtcm_file[0]=rtcm_path[0]='\0'; fp_rtcm=NULL;
    
    for (i=0;i<n;i++) {
        if ((ext=strrchr(infile[i],'.'))&&
            (!strcmp(ext,".rtcm3")||!strcmp(ext,".RTCM3"))) {
            strcpy(rtcm_file,infile[i]);
            init_rtcm(&rtcm);
            break;
        }
    }

2. update_rtcm_ssr函数：位于postpos.c中，每读取一个历元数据，检查是否需要进行SSR数据的更新，核心代码如下：

/* read rtcm file until current time */
    while (timediff(rtcm.time,time)<1E-3) {
        if (input_rtcm3f(&rtcm,fp_rtcm)<-1) break;
        
        /* update ssr corrections */
        for (i=0;i<MAXSAT;i++) {
            if (!rtcm.ssr[i].update||
                rtcm.ssr[i].iod[0]!=rtcm.ssr[i].iod[1]||
                timediff(time,rtcm.ssr[i].t0[0])<-1E-3) continue;
            navs.ssr[i]=rtcm.ssr[i];
            rtcm.ssr[i].update=0;
        }
    }

3. corr_phase_bias_ssr函数：位于postpos.c中，进行相位偏差修正，核心代码如下：

/* carrier-phase bias correction by ssr --------------------------------------*/
static void corr_phase_bias_ssr(obsd_t *obs, int n, const nav_t *nav)
{
    double freq;
    uint8_t code;
    int i,j;
    
    for (i=0;i<n;i++) for (j=0;j<NFREQ;j++) {
        code=obs[i].code[j];
        
        if ((freq=sat2freq(obs[i].sat,code,nav))==0.0) continue;
        
        /* correct phase bias (cyc) */
        obs[i].L[j]-=nav->ssr[obs[i].sat-1].pbias[code-1]*freq/CLIGHT;
    }
}

4. corr_meas函数：位于ppp.c中，进行伪距偏差修正，核心代码如下：

        if (opt->sateph==EPHOPT_SSRAPC||opt->sateph==EPHOPT_SSRCOM) {
            /* select SSR code correction based on code */
            if (sys==SYS_GPS)
                ix=(i==0?CODE_L1W-1:CODE_L2W-1);
            else if (sys==SYS_GLO)
                ix=(i==0?CODE_L1P-1:CODE_L2P-1);
            else if (sys==SYS_GAL)
                ix=(i==0?CODE_L1X-1:CODE_L7X-1);
            /* apply SSR correction */
            P[i]+=(nav->ssr[obs->sat-1].cbias[obs->code[i]-1]-nav->ssr[obs->sat-1].cbias[ix]);
        }

5. satpos_ssr函数：位于ephemeris.c中，利用星历文件和SSR修正数据进行卫星位置计算，核心代码如下：

    /* ssr orbit and clock correction (ref [4]) */
    if (fabs(t1)>MAXAGESSR||fabs(t2)>MAXAGESSR) {
        trace(2,"age of ssr error: %s sat=%2d t=%.0f %.0f\n",time_str(time,0),
              sat,t1,t2);
        *svh=-1;
        return 0;
    }
    if (ssr->udi[0]>=1.0) t1-=ssr->udi[0]/2.0;
    if (ssr->udi[1]>=1.0) t2-=ssr->udi[1]/2.0;
    
    for (i=0;i<3;i++) deph[i]=ssr->deph[i]+ssr->ddeph[i]*t1;
    dclk=ssr->dclk[0]+ssr->dclk[1]*t2+ssr->dclk[2]*t2*t2;
    
    /* ssr highrate clock correction (ref [4]) */
    if (ssr->iod[0]==ssr->iod[2]&&ssr->t0[2].time&&fabs(t3)<MAXAGESSR_HRCLK) {
        dclk+=ssr->hrclk;
    }
    if (norm(deph,3)>MAXECORSSR||fabs(dclk)>MAXCCORSSR) {
        trace(3,"invalid ssr correction: %s deph=%.1f dclk=%.1f\n",
              time_str(time,0),norm(deph,3),dclk);
        *svh=-1;
        return 0;
    }
    /* satellite postion and clock by broadcast ephemeris */
    if (!ephpos(time,teph,sat,nav,ssr->iode,rs,dts,var,svh)) return 0;
    
    /* satellite clock for gps, galileo and qzss */
    sys=satsys(sat,NULL);
    if (sys==SYS_GPS||sys==SYS_GAL||sys==SYS_QZS||sys==SYS_CMP) {
        if (!(eph=seleph(teph,sat,ssr->iode,nav))) return 0;
        
        /* satellite clock by clock parameters */
        tk=timediff(time,eph->toc);
        dts[0]=eph->f0+eph->f1*tk+eph->f2*tk*tk;
        dts[1]=eph->f1+2.0*eph->f2*tk;
        
        /* relativity correction */
        dts[0]-=2.0*dot(rs,rs+3,3)/CLIGHT/CLIGHT;
    }
    /* radial-along-cross directions in ecef */
    if (!normv3(rs+3,ea)) return 0;
    cross3(rs,rs+3,rc);
    if (!normv3(rc,ec)) {
        *svh=-1;
        return 0;
    }
    cross3(ea,ec,er);
    
    /* satellite antenna offset correction */
    if (opt) {
        satantoff(time,rs,sat,nav,dant);
    }
    for (i=0;i<3;i++) {
        rs[i]+=-(er[i]*deph[0]+ea[i]*deph[1]+ec[i]*deph[2])+dant[i];
    }
    /* t_corr = t_sv - (dts(brdc) + dclk(ssr) / CLIGHT) (ref [10] eq.3.12-7) */
    dts[0]+=dclk/CLIGHT;
    
    /* variance by ssr ura */
    *var=var_urassr(ssr->ura);

五、SSR修正的RTKLIB PPP测试

我使用rtklibexplore博客中A first look at the u-blox ZED-F9T dual frequency receiver的数据和配置进行测试。SSR数据来自CLK93，里面包含了第一和第二阶段的修正数，但是没有第三阶段的数据。CLK93数据中包含：

使用广播星历+SSR修正后，PPP大约不到20分钟收敛。由于RTKLIB现阶段基本只支持第一阶段修正（卫星轨道、星钟、伪距偏差），精度应该和直接用精密星历、精密星钟、DCB文件效果差不多。如果能够支持第二、第三阶段，相信精度会有更好的提升，期待RTKLIB的未来更新！

参考文献

[1] 从RTK、PPP到PPP-RTK
[2] PPP-RTK market and technology report
[3]RTCM State Space Representation Messages, Status and Plans
[4] Using SSR corrections with RTKLIB for PPP solutions
[5] A first look at the u-blox ZED-F9T dual frequency receiver