在ppp中，rtklib如此对pco和pcv改正，简要mark
开始之前先要把概念搞清楚，根据antex14，如图，
卫星PCO
卫星的pco改正在根据精密星历计算卫星位置时进行，即位于rtklib中的peph2pos函数，通过satantoff求得改正量，计算流程简要说一下，前两步参考rtklib手册p173.
1、计算太阳的ecef坐标；
2、计算星固系到ecef的转换矩阵；
3、计算消电离层组合系数；
4、得到消电离层组合的pco。
程序贴一下：
/* satellite antenna phase center offset ---------------------------------------
* compute satellite antenna phase center offset in ecef
* args : gtime_t time I time (gpst)
* double *rs I satellite position and velocity (ecef)
* {x,y,z,vx,vy,vz} (m|m/s)
* int sat I satellite number
* nav_t *nav I navigation data
* double *dant I satellite antenna phase center offset (ecef)
* {dx,dy,dz} (m) (iono-free LC value)
* return : none
*-----------------------------------------------------------------------------*/
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 [ex, ey, ez] 星固系到ecef旋转矩阵 */
for (i=0;i<3;i++) r[i]=-rs[i];
if (!normv3(r,ez)) return; /* 星固系z轴单位矢量 */
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; /* 星固系y轴单位矢量 */
cross3(ey,ez,ex); /* 星固系x轴单位矢量 */
/* 卫星所属系统为伽利略或sbas */
if (NFREQ>=3&&(satsys(sat,NULL)&(SYS_GAL|SYS_SBS))) k=2;
if (NFREQ<2||lam[j]==0.0||lam[k]==0.0) return;
/* 消电离层组合系数C1 C2 */
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;/* 消电离层组合得到的ecef下的pco */
}
}
注意，求出卫星的pco后加到卫星坐标，注意跟前面图中的公式比对。
卫星PCV
卫星的pcv改正在ppp中计算载波和伪距的残余函数ppp_res中进行，改正函数为satantpcv，该函数的作用时计算卫星的天底角"Nadir Angle"，
/* satellite antenna phase center variation ----------------------------------*/
/* rs 卫星坐标
rr 接收机坐标
pcv 读取的atx文件中有效的改正值
dant 改正量
*/
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];
}
/* eu为ru单位矢量，指向接收机；ez为rz单位矢量，指向地心 */
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);
}
然后调用antmodel_s，计算各个频率的天线相位中心变化改正值，
/* satellite antenna model ------------------------------------------------------
* compute satellite antenna phase center parameters
* args : pcv_t *pcv I antenna phase center parameters
* double nadir I nadir angle for satellite (rad)
* double *dant O range offsets for each frequency (m)
* return : none
*-----------------------------------------------------------------------------*/
extern void antmodel_s(const pcv_t *pcv, double nadir, double *dant)
{
int i;
trace(4,"antmodel_s: nadir=%6.1f\n",nadir*R2D);
/* 计算各个频率的天线相位中心变化改正值 */
for (i=0;i<NFREQ;i++) {
dant[i]=interpvar(nadir*R2D*5.0,pcv->var[i]); /* 插值 */
}
trace(5,"antmodel_s: dant=%6.3f %6.3f\n",dant[0],dant[1]);
}
这里天底角（关于天底角见rtklib手册p173, F.8-2）由弧度转化为角度，因为atx文件中都是按照角度给定变化值的。然后调用计算pcv：
/* interpolate antenna phase center variation --------------------------------*/
static double interpvar(double ang, const double *var)
{
double a=ang/5.0; /* ang=0-90 */
int i=(int)a;
if (i<0) return var[0]; else if (i>=18) return var[18];
return var[i]*(1.0-a+i)+var[i+1]*(a-i);
}
关于返回值pcv，是这样计算的，对于如下给出的值：
如果天底角为10.2度，那么pcv为
0.8∗v10+0.2∗v11=0.8∗0.70+0.2∗0.000.8*v_{10}+0.2*v_{11}=0.8*0.70+0.2*0.000.8∗v10+0.2∗v11=0.8∗0.70+0.2∗0.00
另外注意计算出来后这样改正到伪距载波：
L[i]=obs->L[i]*lam[i]-dants[i]-dantr[i]-phw*lam[i];
P[i]=obs->P[i] -dants[i]-dantr[i];
还有，rtklib只顾及卫星信号的天顶距而不考虑信号方位角变化时天线相位中心的变化。也就是说，NOAZI下面这些都没有用！
接收机PCO和PCV
接收机的这两项改正都是放在antmodel函数改正的。贴：
/* receiver antenna model ------------------------------------------------------
* compute antenna offset by antenna phase center parameters
* args : pcv_t *pcv I antenna phase center parameters
* double *azel I azimuth/elevation for receiver {az,el} (rad)
* int opt I option (0:only offset,1:offset+pcv)
* double *dant O range offsets for each frequency (m)
* return : none
* notes : current version does not support azimuth dependent terms
*-----------------------------------------------------------------------------*/
/* del 用户在处理之前设定的天线偏差 */
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++) {
/* 天线文件中的pco + 用户设定的天线偏差 */
for (j=0;j<3;j++) off[j]=pcv->off[i][j]+del[j];
/* pco同pcv归化到距离来改正 */
dant[i]=-dot(off,e,3)+(opt?interpvar(90.0-azel[1]*R2D,pcv->var[i]):0.0);
}
trace(5,"antmodel: dant=%6.3f %6.3f\n",dant[0],dant[1]);
}
这里可能会不明白怎么把相位中心偏差归算到距离。atx中给出的是东北天下的接收机天线相位中心偏差，那么用东北天下接收机天线到卫星的单位矢量，自然就可以把偏差归算到接收机天线到卫星的连线上。至于为什么为负值，可以这样理解，根据前面图中的公式，我们将其简化为：
rcenter=rARP+rpco
r_{center} = r_{ARP} + r_{pco}
rcenter=rARP+rpco
假设卫星位置无偏差，为rsatr_{sat}rsat，那么：
rsat−rcenter=rsat−rARP−rpco
r_{sat} - r_{center} = r_{sat} - r_{ARP} - r_{pco}
rsat−rcenter=rsat−rARP−rpco
化为距离，显然方程左侧前两项就是观测的伪距（假设没有其他误差），则：
ρtrue=ρmeas−ρpco
\rho_{true} = \rho_{meas} - \rho_{pco}
ρtrue=ρmeas−ρpco
pcv改正与卫星一样，不再赘述。
补充
现在感觉写的挺乱的，再补充一下atx文件相关内容吧，如图
1、DAZI的值表示方位角从0度每隔多少度到360度，比如截图中的是5度，看NOAZI下面是不是0，5，10，…；
2、ZEN1 / ZEN2 / DZEN表示NOAZI后面的值从0度(ZEN1)每隔5度(DZEN)变化到90度(ZEN2)，共19个值，没错吧！
1729
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
