Skip to content

8.1 pntpos():算法入口

图10-1 pntpos() 函数调用关系

8.1.1 pntpos():函数主体

1. 参数列表

c
/* args */
obsd_t   *obs     I   observation data
int       n       I   number of observation data
nav_t    *nav     I   navigation data
prcopt_t *opt     I   processing options
sol_t    *sol     IO  solution
double   *azel    IO  azimuth/elevation angle (rad) (NULL: no output)
ssat_t   *ssat    IO  satellite status              (NULL: no output)
char     *msg     O   error message for error exit
/* return */
int       status  -   1:ok,0:error

2. 执行流程

  • sol->time赋值第一个观测值(顺序存储的第一颗卫星)的时间
  • 如果处理选项不是 SPP(针对 PPP),电离层校正选 Klobuchar 模型 ,对流层校正采用 Saastmoinen 模型
  • 调用satposs()计算卫星位置、速度和钟差(ECEF 坐标系)
    • rs [(0:2)+i*6]:obs[i] 卫星位置 {x, y, z} (m)
    • rs [(3:5)+i*6]:obs[i] 卫星速度 {vx, vy, vz} (m/s)
    • dts[(0:1)+i*2]:obs[i] 卫星钟差与钟漂 {bias, drift} (s, s/s)
    • var[i]:卫星的误差方差 (m^2),代表了卫星位置和钟差所包含的误差水平,最终会合入伪距的方差中
    • svh[i]:卫星健康标志 (-1:correction not available)
  • 调用estpos()使用伪距进行位置估计,采用加权最小二乘,其中会调用 valsol 进行卡方检验和 GDOP 检验
  • estpos()valsol检验失败,即位置估计失败,会调用raim_fde接收机自主完好性监测重新估计,前提是卫星数 > 6、对应参数解算设置opt->posopt[4]=1
  • 调用estvel()使用多普勒进行速度估计
  • 保存方位角和俯仰角信息,赋值卫星状态结构体ssat

3. 注意事项

  • 关于钟漂信息:这里只计算了接收机的钟差,而没有计算接收机的频漂,原因在于 estvel 函数中虽然计算得到了接收机频漂,但并没有将其输出到 sol_t:dtr中。
  • C语言内存使用问题:C语言中用 malloc 申请的内存需要自己调用 free 来予以回收,源码中的 mat 、imat 、zeros 等函数都只是申请了内存,并没有进行内存的回收,在使用这些函数时,用户必须自己调用 free 来回收内存!源码中将使用这些函数的代码放置在同一行,在调用函数结尾处也统一进行内存回收,位置较为明显,不致于轻易忘记。

4. 问题思考

  • RTKLIB观测数据中的时间:源码中将 obs[0].time作为星历选择时间传递给 satposs函数,这样对于每一颗观测卫星,都要使用第一颗观测卫星的数据接收时间作为选择星历的时间标准。是否应该每颗卫星都使用自己的观测时间?或者应该使用每颗卫星自己的信号发射时间?还是说这点差别对选择合适的星历其实没有关系?该问题的分析参考附录B.1
  • raim_fde 对卫星数目的要求:这里规定能够执行 raim_fde 函数的前提是数目大于等于 6,感觉不是只要大于等于 5就可以了吗?该问题的分析参考附录B.2

5. 源码注释

点击查看代码
c
extern int pntpos(const obsd_t *obs, int n, const nav_t *nav,
                  const prcopt_t *opt, sol_t *sol, double *azel, ssat_t *ssat,
                  char *msg)
{
    prcopt_t opt_=*opt;
    double *rs,*dts,*var,*azel_,*resp;
    int i,stat,vsat[MAXOBS]={0},svh[MAXOBS];
    
    trace(3,"pntpos  : tobs=%s n=%d\n",time_str(obs[0].time,3),n);
    sol->stat=SOLQ_NONE;
    
    if (n<=0) {
        strcpy(msg,"no observation data");
        return 0;
    }

    sol->time=obs[0].time; // sol->time赋值第一个观测值的时间
    msg[0]='\0';
    
    rs=mat(6,n); dts=mat(2,n); var=mat(1,n); azel_=zeros(2,n); resp=mat(1,n);
    
    if (opt_.mode!=PMODE_SINGLE) { /* for precise positioning */
        opt_.ionoopt=IONOOPT_BRDC; // 电离层校正选 Klobuchar 广播星历模型
        opt_.tropopt=TROPOPT_SAAS; // 对流层校正采用 Saastmoinen 模型
    }
    
    // 计算卫星位置、速度和钟差
    /* satellite positons, velocities and clocks */
    satposs(sol->time,obs,n,nav,opt_.sateph,rs,dts,var,svh);
    
    // 用伪距位置估计,加权最小二乘,其中会调用 valsol 进行卡方检验和 GDOP 检验
    /* estimate receiver position with pseudorange */
    stat=estpos(obs,n,rs,dts,var,svh,nav,&opt_,sol,azel_,vsat,resp,msg);
    
    /* RAIM FDE */
    if (!stat&&n>=6&&opt->posopt[4]) {  
        // estpos 中 valsol 检验失败,即位置估计失败,则会调用 RAIM 接收机自主完好性监测重新估计,
        // 前提是卫星数 > 6,且对应参数解算设置 opt->posopt[4]=1,即pos1-posopt5=on(RAIM FDE)
        stat=raim_fde(obs,n,rs,dts,var,svh,nav,&opt_,sol,azel_,vsat,resp,msg);
    }
    // 调用 estvel 进行多普勒速度估计
    /* estimate receiver velocity with Doppler */
    if (stat) {
        estvel(obs,n,rs,dts,nav,&opt_,sol,azel_,vsat);
    }
    if (azel) {
        for (i=0;i<n*2;i++) azel[i]=azel_[i];
    }
    if (ssat) { // 保存卫星状态信息至ssat结构体
        for (i=0;i<MAXSAT;i++) {
            ssat[i].vs=0;
            ssat[i].azel[0]=ssat[i].azel[1]=0.0;
            ssat[i].resp[0]=ssat[i].resc[0]=0.0;
            ssat[i].snr[0]=0;
        }
        for (i=0;i<n;i++) {
            ssat[obs[i].sat-1].azel[0]=azel_[  i*2];
            ssat[obs[i].sat-1].azel[1]=azel_[1+i*2];
            ssat[obs[i].sat-1].snr[0]=obs[i].SNR[0];
            if (!vsat[i]) continue;
            ssat[obs[i].sat-1].vs=1;
            ssat[obs[i].sat-1].resp[0]=resp[i];
        }
    }
    free(rs); free(dts); free(var); free(azel_); free(resp);
    return stat;
}