RTKLIB源码学习之单点定位pntpos.c

本博客是转载学习之用，感谢：http://t.csdnimg.cn/rC8iC 

pntpos流程图

/* single-point positioning ----------------------------------------------------
* compute receiver position, velocity, clock bias by single-point positioning
* with pseudorange and doppler observables
* 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 : status(1:ok,0:error)
* notes  : assuming sbas-gps, galileo-gps, qzss-gps, compass-gps time offset and
*          receiver bias are negligible (only involving glonass-gps time offset
*          and receiver bias)
*-----------------------------------------------------------------------------*/
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)

pntpos函数利用多普勒频移测量和伪距进行单点定位，给出接收机位置，速度和钟差

主程序

if (n<=0) {strcpy(msg,"no observation data"); return 0;}//判断卫星个数是否大于零

   if (opt_.mode!=PMODE_SINGLE) { /* for precise positioning */
#if 0
        opt_.sateph =EPHOPT_BRDC;//可自己进行调试
#endif
        opt_.ionoopt=IONOOPT_BRDC;
        opt_.tropopt=TROPOPT_SAAS;
    }

处理选项opt中的模式不是单点模式时，电离层校正采用broadcast model，对流层校正采用SSaastamoinen model

satposs函数

/* satellite positions and clocks ----------------------------------------------
* compute satellite positions, velocities and clocks
* args   : gtime_t teph     I   time to select ephemeris (gpst)
*          obsd_t *obs      I   observation data
*          int    n         I   number of observation data
*          nav_t  *nav      I   navigation data
*          int    ephopt    I   ephemeris option (EPHOPT_???)
*          double *rs       O   satellite positions and velocities (ecef)
*          double *dts      O   satellite clocks
*          double *var      O   sat position and clock error variances (m^2)
*          int    *svh      O   sat health flag (-1:correction not available)
* return : none
* notes  : rs [(0:2)+i*6]= obs[i] sat position {x,y,z} (m)
*          rs [(3:5)+i*6]= obs[i] sat velocity {vx,vy,vz} (m/s)
*          dts[(0:1)+i*2]= obs[i] sat clock {bias,drift} (s|s/s)
*          var[i]        = obs[i] sat position and clock error variance (m^2)
*          svh[i]        = obs[i] sat health flag
*          if no navigation data, set 0 to rs[], dts[], var[] and svh[]
*          satellite position and clock are values at signal transmission time
*          satellite position is referenced to antenna phase center
*          satellite clock does not include code bias correction (tgd or bgd)
*          any pseudorange and broadcast ephemeris are always needed to get
*          signal transmission time
*-----------------------------------------------------------------------------*/
extern void satposs(gtime_t teph, const obsd_t *obs, int n, const nav_t *nav,
                    int ephopt, double *rs, double *dts, double *var, int *svh)

 satposs函数位于ephemeris.c，按所观测到的卫星顺序计算每颗卫星的位置和速度（rs）、钟差和频漂（dts）

for (i=0;i<n&&i<MAXOBS;i++) {
        //按观测数据的顺序，将当前观测卫星的rs、dts、var和svh数组的元素设置为零
        for (j=0;j<6;j++) rs [j+i*6]=0.0;
        for (j=0;j<2;j++) dts[j+i*2]=0.0;
        var[i]=0.0; svh[i]=0;
        
        /* 通过判断某一频率下信号的伪距是否为 0，来得到此时所用的频率个数。
注意，频率个数不能大于 NFREQ（默认为 3）*/
        for (j=0,pr=0.0;j<NFREQ;j++) if ((pr=obs[i].P[j])!=0.0) break;
        if (j>=NFREQ) {
            trace(2,"no pseudorange %s sat=%2d\n",time_str(obs[i].time,3),obs[i].sat);
            continue;//处理下一个观测数据
        }

        /* 用数据接收时间减去伪距信号传播时间，得到卫星信号的发射时间*/
        time[i]=timeadd(obs[i].time,-pr/CLIGHT);
        
        /* 调用 ephclk函数，由广播星历计算出当前观测卫星的钟差。
注意，此时的钟差是没有考虑相对论效应和TGD的*/
        if (!ephclk(time[i],teph,obs[i].sat,nav,&dt)) {
            trace(2,"no broadcast clock %s sat=%2d\n",time_str(time[i],3),obs[i].sat);
            continue;
        }

        /*用得到的信号发射时间减去卫星钟差，得到 GPS时间下的卫星信号发射时间 */
        time[i]=timeadd(time[i],-dt);
        
        /* 调用satpos函数，计算信号发射时刻卫星的 P(ecef,m)、V(ecef,m/s)、C((s|s/s))。注意，这里计算出的钟差是考虑了相对论效应，没考虑TGD */
        if (!satpos(time[i],teph,obs[i].sat,ephopt,nav,rs+i*6,dts+i*2,var+i,
                    svh+i)) {
            trace(2,"no ephemeris %s sat=%2d\n",time_str(time[i],3),obs[i].sat);
            continue;
        }

        /* 如果上一步计算出的钟偏为0，就再次调用ephclk函数，将其计算出的卫星钟偏作为最终的结果 */
        if (dts[i*2]==0.0) {
            if (!ephclk(time[i],teph,obs[i].sat,nav,dts+i*2)) continue;
            dts[1+i*2]=0.0;
            *var=SQR(STD_BRDCCLK);
        }
    }

注：TGD（time group delay）：卫星内部信号从产生到发射的时延之差 ，采用双频接收机时不必要采用。

          clock {bias,drift} (s|s/s) 频漂drift单位？？？

ephclk 函数利用广播星历计算卫星钟差

static int ephclk(gtime_t time, gtime_t teph, int sat, const nav_t *nav,
                  double *dts)
{
    eph_t  *eph;
    geph_t *geph;
    seph_t *seph;
    int sys;
    
    trace(4,"ephclk  : time=%s sat=%2d\n",time_str(time,3),sat);
    
    sys=satsys(sat,NULL);//satsys函数根据卫星编号确定该卫星所属的导航系统和该卫星的PRN编号
    
    if (sys==SYS_GPS||sys==SYS_GAL||sys==SYS_QZS||sys==SYS_CMP) {
        if (!(eph=seleph(teph,sat,-1,nav))) return 0;
        //调用seleph函数来选择toe值与星历选择时间标准teph最近的那个星历
        *dts=eph2clk(time,eph);//eph2clk函数通过广播星历和信号发射时间计算卫星钟差
    }
    else if (sys==SYS_GLO) {
        if (!(geph=selgeph(teph,sat,-1,nav))) return 0;
        *dts=geph2clk(time,geph);
    }
    else if (sys==SYS_SBS) {
        if (!(seph=selseph(teph,sat,nav))) return 0;
        *dts=seph2clk(time,seph);
    }
    else return 0;
    
    return 1;
}

satpos函数计算卫星位置

/* satellite position and clock ------------------------------------------------
* compute satellite position, velocity and clock
* args   : gtime_t time     I   time (gpst)
*