RTKLIB专题学习(七)—精密单点定位实现初识（三）
上两篇我们介绍了RTKLIB中精密单点定位的大致流程，今天我们对照RTKLIB学习手册，来学习相应改正公式和误差源
1.在PPP模式中
RTKLIB也使用类似于附录E.7中描述的方案的EKF估计过程。与基线处理的不同，PPP采用ZD(非差)测量方程，就像单点定位模型，而不是双差模型。
2.利用E3.2和E3.4，可以推导出非差观测值的无电离层组合模型：
组合后的非差观测模型如下：
图中
为L1和L2载波相位组成的无电离层组合的改正项，包括天线相位中心偏差和变化、固体潮引起的位置偏移、相位缠绕。
2.对于接收机天线相位中心模型
分为PCV和PCO改正；PCO是天线参考点与天线相位中心的相对位置；PCV定义为天线取决于高度角和方位角的多余相位延迟；
示意图如下：
计算公式如下：
其中公式中的PCV和PCO值可以来自ANTEX format的数据，IGS定期发布天线相位中心改正文件，提供卫星在星固系中的天线相位中心偏差值。
相应源码如下：
/* receiver antenna model ------------------------------------------------------
* compute antenna offset by antenna phase center parameters
* args : pcv_t *pcv I antenna phase center parameters
* double *del I antenna delta {e,n,u} (m)
* 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
*-----------------------------------------------------------------------------*/
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(5,"antmodel: dant=%6.3f %6.3f\n",dant[0],dant[1]);
}
通过option可以设置选择修改天线相位中心偏差或变化改正；pco是指天线质心或几何中心与相位中心的偏差；pco改正需要将其由站心地平坐标系转换到地心地固坐标系下
3.卫星天线相位中心模型
类似于接收机的PCV和PCO的定义，示意图如下：
由于igs提供的PCO为星固坐标系下的 值，因此需要转换为地心地固坐标系下的数值，示意图及计算公式如下：
源码部分如下：
/* 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);
}
其中调用的函数为：
/* 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]);
}
4.固体潮引起的测站位移
固体潮现象使得测站的实际坐标随时间周期性变化，测站垂直方向位移最大可达数十厘米
源代码中的潮汐引起的位移改正如下：
/* tidal displacement ----------------------------------------------------------
* displacements by earth tides
* args : gtime_t tutc I time in utc
* double *rr I site position (ecef) (m)
* int opt I options (or of the followings)
* 1: solid earth tide
* 2: ocean tide loading
* 4: pole tide
* 8: elimate permanent deformation
* double *erp I earth rotation parameters (NULL: not used)
* double *odisp I ocean loading parameters (NULL: not used)
* odisp[0+i*6]: consituent i amplitude radial(m)
* odisp[1+i*6]: consituent i amplitude west (m)
* odisp[2+i*6]: consituent i amplitude south (m)
* odisp[3+i*6]: consituent i phase radial (deg)
* odisp[4+i*6]: consituent i phase west (deg)
* odisp[5+i*6]: consituent i phase south (deg)
* (i=0:M2,1:S2,2:N2,3:K2,4:K1,5:O1,6:P1,7:Q1,
* 8:Mf,9:Mm,10:Ssa)
* double *dr O displacement by earth tides (ecef) (m)
* return : none
* notes : see ref [1], [2] chap 7
* see ref [4] 5.2.1, 5.2.2, 5.2.3
* ver.2.4.0 does not use ocean loading and pole tide corrections
*-----------------------------------------------------------------------------*/
extern void tidedisp(gtime_t tutc, const double *rr, int opt, const erp_t *erp,
const double *odisp, double *dr)
{
gtime_t tut;
double pos[2],E[9],drt[3],denu[3],rs[3],rm[3],gmst,erpv[5]={0};
int i;
#ifdef IERS_MODEL
double ep[6],fhr;
int year,mon,day;
#endif
trace(3,"tidedisp: tutc=%s\n",time_str(tutc,0));
if (erp) {
geterp(erp,utc2gpst(tutc),erpv);
}
tut=timeadd(tutc,erpv[2]);
dr[0]=dr[1]=dr[2]=0.0;
if (norm(rr,3)<=0.0) return;
pos[0]=asin(rr[2]/norm(rr,3));
pos[1]=atan2(rr[1],rr[0]);
xyz2enu(pos,E);
if (opt&1) { /* solid earth tides */
/* sun and moon position in ecef */
sunmoonpos(tutc,erpv,rs,rm,&gmst);
#ifdef IERS_MODEL
time2epoch(tutc,ep);
year=(int)ep[0];
mon =(int)ep[1];
day =(int)ep[2];
fhr =ep[3]+ep[4]/60.0+ep[5]/3600.0;
/* call DEHANTTIDEINEL */
dehanttideinel_((double *)rr,&year,&mon,&day,&fhr,rs,rm,drt);
#else
tide_solid(rs,rm,pos,E,gmst,opt,drt);
#endif
for (i=0;i<3;i++) dr[i]+=drt[i];
}
if ((opt&2)&&odisp) { /* ocean tide loading */
tide_oload(tut,odisp,denu);
matmul("TN",3,1,3,1.0,E,denu,0.0,drt);
for (i=0;i<3;i++) dr[i]+=drt[i];
}
if ((opt&4)&&erp) { /* pole tide */
tide_pole(tut,pos,erpv,denu);
matmul("TN",3,1,3,1.0,E,denu,0.0,drt);
for (i=0;i<3;i++) dr[i]+=drt[i];
}
trace(5,"tidedisp: dr=%.3f %.3f %.3f\n",dr[0],dr[1],dr[2]);
}
其中包括，固体潮、大洋负荷以及极潮改正
5.相位缠绕改正
GNSS卫星信号采用右极化方式，这种极化方式的信号使得观测到的载波相位与卫星、接收机天线的朝向有关，其对精密单点定位的影响可达半个波长
6.最后，由下面函数进行码和相位测量值的校正
可以看到，对于伪距测量值，只进行天线相位中心改正；载波相位还进行相位缠绕改正
/* antenna corrected measurements --------------------------------------------*/
static void corr_meas(const obsd_t *obs, const nav_t *nav, const double *azel,
const prcopt_t *opt, const double *dantr,
const double *dants, double phw, double *L, double *P,
double *Lc, double *Pc)
{
double freq[NFREQ]={0},C1,C2;
int i,sys=satsys(obs->sat,NULL);
for (i=0;i<NFREQ;i++) {
L[i]=P[i]=0.0;
freq[i]=sat2freq(obs->sat,obs->code[i],nav);
if (freq[i]==0.0||obs->L[i]==0.0||obs->P[i]==0.0) continue;
if (testsnr(0,0,azel[1],obs->SNR[i]*SNR_UNIT,&opt->snrmask)) continue;
/* 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];
/* P1-C1,P2-C2 dcb correction (C1->P1,C2->P2) */
if (sys==SYS_GPS||sys==SYS_GLO) {
if (obs->code[i]==CODE_L1C) P[i]+=nav->cbias[obs->sat-1][1];
if (obs->code[i]==CODE_L2C) P[i]+=nav->cbias[obs->sat-1][2];
}
}
/* iono-free LC */
*Lc=*Pc=0.0;
if (freq[0]==0.0||freq[1]==0.0) return;
C1= SQR(freq[0])/(SQR(freq[0])-SQR(freq[1]));
C2=-SQR(freq[1])/(SQR(freq[0])-SQR(freq[1]));
if (L[0]!=0.0&&L[1]!=0.0) *Lc=C1*L[0]+C2*L[1];
if (P[0]!=0.0&&P[1]!=0.0) *Pc=C1*P[0]+C2*P[1];
}
总结：
本篇概略介绍了一些RTKLIB中用到的对伪距和载波相位观测值做出的改正，可以了解到一些方法和原理
1089
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
