RTKLIB中的卫星天线与接收机天线修正

天线相位中心即天线接收信号的电气中心，其空间位置在出厂时往往不在天线的几何中心上。天线所辐射出的电磁波在离开天线一定的距离后，其等相位面会近似为一个球面，该球面的球心即为该天线的等效相位中心，即天线相位中心（Antenna Phase Center ）。

NGS提供的ANTEX格式天线模型，包含了卫星天线模型以及部分接收机天线模型。使用天线模型的目的包括：

• 修正天线参考点和天线相位中心的之间的偏差；
• 修正和仰角有关的误差；
• 修正L1和L2之间的相位中心偏差（这个误差可能对整周模糊度固定造成影响）

接收机天线相位中心模型

一般选取接收机天线底部与天线中轴的交点作为参考点（称天线参考点，ARP）：

• ARP与实际相位中心的几何偏差值称为天线相位中心偏差（PCO）
• 由不同高度角、方位角测得的距离产生系统性的测量偏差为天线相位中心变化（PCV）

目前RTKLIB支持NGS PCV以及ANTEX格式的天线模型，其中包括了PCO和PCV修正参数。通过手册E.8章节可知，接收机天线修正如下：

• PCO通常是当地坐标系NED参数，因此需要利用转换矩阵转到ECEF坐标系
  
• PCV修正则通过对高度角进行插值得到
  

卫星天线相位中心模型

卫星天线偏移（Satellite Antenna Offsets）是指卫星的天线相位中心与卫星质心的偏移。在PPP中我们需要用到精密星历。IGS 提供的精密星历(卫星位置)是相对卫星质心的，而载波相位、伪距测量值是相对天线相位中心的，因此我们需要计算PCV和PCO值，将所有值统一转换成相对天线相位中心的。

具体计算公式可以参考手册的E.8章节。

RTKLIB中相关源码

一、读取天线文件、存储相关参数

1. openses函数：位于postpos.c中，调用readpcv函数读取卫星天线参数和接收机天线参数

    /* read satellite antenna parameters */
    if (*fopt->satantp&&!(readpcv(fopt->satantp,pcvs))) {
        showmsg("error : no sat ant pcv in %s",fopt->satantp);
        trace(1,"sat antenna pcv read error: %s\n",fopt->satantp);
        return 0;
    }
    /* read receiver antenna parameters */
    if (*fopt->rcvantp&&!(readpcv(fopt->rcvantp,pcvr))) {
        showmsg("error : no rec ant pcv in %s",fopt->rcvantp);
        trace(1,"rec antenna pcv read error: %s\n",fopt->rcvantp);
        return 0;
    }

2. readpcv函数代码如下：

/* read antenna parameters ------------------------------------------------------*/
extern int readpcv(const char *file, pcvs_t *pcvs)
{
    pcv_t *pcv;
    char *ext;
    int i,stat;
    
    trace(3,"readpcv: file=%s\n",file);
    
    if (!(ext=strrchr(file,'.'))) ext="";
    
    if (!strcmp(ext,".atx")||!strcmp(ext,".ATX")) {
        stat=readantex(file,pcvs);
    }
    else {
        stat=readngspcv(file,pcvs);
    }
    for (i=0;i<pcvs->n;i++) {
        pcv=pcvs->pcv+i;
        trace(4,"sat=%2d type=%20s code=%s off=%8.4f %8.4f %8.4f  %8.4f %8.4f %8.4f\n",
              pcv->sat,pcv->type,pcv->code,pcv->off[0][0],pcv->off[0][1],
              pcv->off[0][2],pcv->off[1][0],pcv->off[1][1],pcv->off[1][2]);
    }
    return stat;
}

3. setpcv函数（设置天线参数）：位于postpos.c中，查找各颗卫星的天线参数，并放在nav->pcvs中；查找接收机的天线参数，并放在popt->pcvr中，需要注意的是，只有当popt->anttype[i]是*(auto)的时候，通过RINEX头识别到的接收机天线类型，以及天线的H/E/N偏差才会起作用。

/* set antenna parameters ----------------------------------------------------*/
static void setpcv(gtime_t time, prcopt_t *popt, nav_t *nav, const pcvs_t *pcvs,
                   const pcvs_t *pcvr, const sta_t *sta)
{
    pcv_t *pcv,pcv0={0};
    double pos[3],del[3];
    int i,j,mode=PMODE_DGPS<=popt->mode&&popt->mode<=PMODE_FIXED;
    char id[64];
    
    /* set satellite antenna parameters */
    for (i=0;i<MAXSAT;i++) {
        nav->pcvs[i]=pcv0;
        if (!(satsys(i+1,NULL)&popt->navsys)) continue;
        if (!(pcv=searchpcv(i+1,"",time,pcvs))) {
            satno2id(i+1,id);
            trace(3,"no satellite antenna pcv: %s\n",id);
            continue;
        }
        nav->pcvs[i]=*pcv;
    }
    for (i=0;i<(mode?2:1);i++) {
        popt->pcvr[i]=pcv0;
        if (!strcmp(popt->anttype[i],"*")) { /* set by station parameters */
            strcpy(popt->anttype[i],sta[i].antdes);
            if (sta[i].deltype==1) { /* xyz */
                if (norm(sta[i].pos,3)>0.0) {
                    ecef2pos(sta[i].pos,pos);
                    ecef2enu(pos,sta[i].del,del);
                    for (j=0;j<3;j++) popt->antdel[i][j]=del[j];
                }
            }
            else { /* enu */
                for (j=0;j<3;j++) popt->antdel[i][j]=stas[i].del[j];
            }
        }
        if (!(pcv=searchpcv(0,popt->anttype[i],time,pcvr))) {
            trace(2,"no receiver antenna pcv: %s\n",popt->anttype[i]);
            *popt->anttype[i]='\0';
            continue;
        }
        strcpy(popt->anttype[i],pcv->type);
        popt->pcvr[i]=*pcv;
    }
}

二、卫星端PCO修正

satantoff函数：位于preceph.c中，分别计算L1和L2的卫星端PCO。由于PPP中使用的是无电离层组合，因此会计算无电离层组合后的PCO修正值。之后在peph2pos函数中，会在由精密星历计算的卫星位置上进行PCO修正。

extern void satantoff(gtime_t time, const double *rs, int sat, const nav_t *nav,
                      double *dant)
{
    const double *lam=nav->lam[sat-1];
    const pcv_t *pcv=nav->pcvs+sat-1;
    double ex[3],ey[3],ez[3],es[3],r[3],rsun[3],gmst,erpv[5]={0};
    double gamma,C1,C2,dant1,dant2;
    int i,j=0,k=1;
    
    trace(4,"satantoff: time=%s sat=%2d\n",time_str(time,3),sat);
    
    /* sun position in ecef */
    sunmoonpos(gpst2utc(time),erpv,rsun,NULL,&gmst);
    
    /* unit vectors of satellite fixed coordinates */
    for (i=0;i<3;i++) r[i]=-rs[i];
    if (!normv3(r,ez)) return;
    for (i=0;i<3;i++) r[i]=rsun[i]-rs[i];
    if (!normv3(r,es)) return;
    cross3(ez,es,r);
    if (!normv3(r,ey)) return;
    cross3(ey,ez,ex);
    
    if ((NFREQ>=3)&&(satsys(sat,NULL)&SYS_SBS)) k=2;  /* always use L5 for SBAS */
    if ((NFREQ>=3)&&(satsys(sat,NULL)&SYS_GAL)&&(lam[1]==0)) k=2; /* use E5a if no E5b for GAL */
    if (NFREQ<2||lam[j]==0.0||lam[k]==0.0) return;
    
    gamma=SQR(lam[k])/SQR(lam[j]);
    C1=gamma/(gamma-1.0);
    C2=-1.0 /(gamma-1.0);
    
    /* iono-free LC */
    for (i=0;i<3;i++) {
        dant1=pcv->off[j][0]*ex[i]+pcv->off[j][1]*ey[i]+pcv->off[j][2]*ez[i];
        dant2=pcv->off[k][0]*ex[i]+pcv->off[k][1]*ey[i]+pcv->off[k][2]*ez[i];
        dant[i]=C1*dant1+C2*dant2;
    }
}

三、卫星端PCV修正

satantpcv函数： 计算卫星天线端PCV修正值，具体计算方法可以参考手册E.8：

/* satellite antenna phase center variation ----------------------------------*/
static void satantpcv(const double *rs, const double *rr, const pcv_t *pcv,
                      double *dant)
{
    double ru[3],rz[3],eu[3],ez[3],nadir,cosa;
    int i;
    
    for (i=0;i<3;i++) {
        ru[i]=rr[i]-rs[i];
        rz[i]=-rs[i];
    }
    if (!normv3(ru,eu)||!normv3(rz,ez)) return;
    
    cosa=dot(eu,ez,3);
    cosa=cosa<-1.0?-1.0:(cosa>1.0?1.0:cosa);
    nadir=acos(cosa);
    
    antmodel_s(pcv,nadir,dant);
}

四、接收机端PCO和PCV修正

antmodel函数计算接收机端修正值，需要注意的是，代码中实际计算方法和手册不太一致。代码先计算line of sight在NED中的单位矢量e，然后计算PCO在line of sight上的投影长度，并没有将PCO从NED转换到ECEF这一过程。

extern void antmodel(const pcv_t *pcv, const double *del, const double *azel,
                     int opt, double *dant)
{
    double e[3],off[3],cosel=cos(azel[1]);
    int i,j;
    
    trace(4,"antmodel: azel=%6.1f %4.1f opt=%d\n",azel[0]*R2D,azel[1]*R2D,opt);
    
    e[0]=sin(azel[0])*cosel;
    e[1]=cos(azel[0])*cosel;
    e[2]=sin(azel[1]);
    
    for (i=0;i<NFREQ;i++) {
        for (j=0;j<3;j++) off[j]=pcv->off[i][j]+del[j];
        
        dant[i]=-dot(off,e,3)+(opt?interpvar(90.0-azel[1]*R2D,pcv->var[i]):0.0);
    }
    trace(4,"antmodel: dant=%6.3f %6.3f\n",dant[0],dant[1]);
}

五、PPP中的修正流程

ppp_res函数：位于ppp.c中，通过satantpcv函数计算卫星天线端PCV修正，通过antmodel函数计算接收机天线端修正值，最后在corr_meas函数中进行修正。

        /* satellite and receiver antenna model */
        if (opt->posopt[0]) satantpcv(rs+i*6,rr,nav->pcvs+sat-1,dants);
        antmodel(opt->pcvr,opt->antdel[0],azel+i*2,opt->posopt[1],dantr);
        
        /* phase windup model */
        if (!model_phw(rtk->sol.time,sat,nav->pcvs[sat-1].type,
                       opt->posopt[2]?2:0,rs+i*6,rr,&rtk->ssat[sat-1].phw)) {
            continue;
        }
        /* corrected phase and code measurements */
        corr_meas(obs+i,nav,azel+i*2,&rtk->opt,dantr,dants,
                  rtk->ssat[sat-1].phw,L,P,&Lc,&Pc);

在corr_meas函数中进行修正，如下所示，在伪距、载波相位上减去了dants和dantr。

        /* antenna phase center and phase windup correction */
        L[i]=obs->L[i]*CLIGHT/freq[i]-dants[i]-dantr[i]-phw*CLIGHT/freq[i];
        P[i]=obs->P[i]       -dants[i]-dantr[i];

六、rtkpost中的修正流程

zdres函数:位于rtkpos.c中，调用antmodel函数计算接收机端修正值，在zdres_sat函数中进行修正。需要注意的是，相对定位只需要计算接收机端的天线相位中心修正值。这是由于，相对定位进行单差时，已经将卫星端的天线误差消除了。但是PPP则需要同时处理接收机和卫星端误差。

        /* calc receiver antenna phase center correction */
        antmodel(opt->pcvr+index,opt->antdel[index],azel+i*2,opt->posopt[1],
                 dant);
        
        /* calc undifferenced phase/code residual for satellite */
        zdres_sat(base,r,obs+i,nav,azel+i*2,dant,opt,y+i*nf*2,freq+i*nf);

参考文献

[1] GPS从入门到放弃（二十） — 天线偏移
[2] 天线高改正和天线相位中心改正
[3] 101-RTKLIB中的天线相位中心改正