天线相位中心改正
文章目录
1. 相关概念：
- 天线参考点：Antenna Reference Point, ARP
- 天线瞬时相位中心： Instantaneous Phase Center, IPC
- 相位中心偏差：Phase Center Offset, PCO
- 天线相位中心变化：Phase Center Variation, PCV
- 天线相位中心分为两部分：
- 天线的平均相位中心（天线瞬时相位中心的平均值）与天线参考点ARP之间的偏差，成为天线相位中心偏差（PCO）
- GNSS信号实际测量的时天线瞬时相位中心，与平均相位中心的差值，被称为天线相位中心变化（PCV）
2. 天线相位中心改正
2.1 卫星天线相位中心改正
-
GPS测量测定的是从卫星发射天线相位中心至接收机相位中心间的距离，而IGS给出的星历是卫星质心的坐标，这之间的偏差被称为卫星天线相位中心偏差（PCO）
-
改正方法：
X ⃗ p c = X ⃗ m c + [ x ^ , y ^ , z ^ ] ∗ r ⃗ \vec{X}_{pc}=\vec{X}_{mc}+\begin{bmatrix}\hat{x},\hat{y},\hat{z}\end{bmatrix}*\vec{r} Xpc=Xmc+[x^,y^,z^]∗r
其中，X p c X_{pc} Xpc：卫星天线相位中心的位置矢量
X m c X_{mc} Xmc：卫星质心的位置向量
r ⃗ = [ Δ x , Δ y , Δ z ] \vec{r}=\begin{bmatrix}\Delta{x},\Delta{y},\Delta{z}\end{bmatrix} r=[Δx,Δy,Δz]：卫星天线相位中心的偏差矢量
[ x ^ , y ^ , z ^ ] \begin{bmatrix}\hat{x},\hat{y},\hat{z}\end{bmatrix} [x^,y^,z^]：卫星的姿态矩阵， z ^ \hat{z} z^是单位向量，从卫星质心指向地心， x ^ \hat{x} x^为单位向量，在地-星-日平面内垂直于 z ^ \hat{z} z^并指向太阳， y ^ \hat{y} y^与 x ^ \hat{x} x^、 z ^ \hat{z} z^组成右手系
-
当卫星高度角、方位角不断变动时，其天线相位中心也发生改变，造成的误差即为天线相位中心变化量（PCV）
-
PCO改正方式：卫星质心的位置=卫星天线相位中心的位置-PCO
-
PCV改正方式：卫星至接收机间的几何距离=观测距离-PCV+其他改正
2.2 接收机天线相位中心改正
-
接收机的天线相位中心与天线的参考点ARP不一致，另外仪器在对中、量天线高时参考的是天线参考点，故接收机天线相位中心改正也分为PCO和PCV
-
接收机天线PCO改正：天线参考点位置=天线相位中心的位置-天线相位中心偏差PCO
求得天线参考点后，再通过天线对中数据（是否有偏心）及仪器高等数据求得标石中心位置
-
接收机相位中心变化PCV通过两种形式给出
- 只顾及卫星信号的天顶距而不考虑信号方位角变化（PCV NOAZI）
- 同时顾及信号的天顶距以及方位角时的天象相位中心变化（PCV AZEL）
-
天线相位中心变化PCV通常是用来改正距离观测值的，公式如下：
几何距离=观测的距离-PCV+其他改正项
3. IGS提供的天线相位改正文件
3.1 绝对天线相位中心模型
- 从2006年11月，GPS周1400周，IGS开始使用绝对天线相位中心模型，该模型不仅考虑了天线相位中心偏差，还考虑了方位角、高度角等引起的天线相位变化。在这以前一直采用的是相对天线相位改正模型。
- 绝对相位中心模型中的天线相位中心偏差和相位中心变化通常可采用下列两种方法来进行测定：
- 第一种方法是在微波暗室中用微波信号发生器所产生的模拟GPS信号来对接收机天线进行检测；
- 第二种方法是在室外利用真正的GPS信号，通过自动机器人将接收机天线倾斜、旋转，从而来测定接收机天线的相位中心偏差和相位中心变化。
3.2 atx文件介绍
上图分别为一条卫星天线和接收机记录，其中：
-
每条记录开始于START OF ANTENNA ，结束于END OF ANTENNA（接收机天线记录部分由于篇幅原因未展示）
-
TYPE / SERIAL NO
- 对接收机而言分别为：Antenna type、Serial number(blank: all representatives of the specified antenna type）、blank 、blank
- 对卫星而言分别为：Antenna type、Satellite code “sNN”、Satellite code “sNNN”（其中s表示satellite system flag， NNN表示SVN number ）
-
METH / BY / # / DATE
- 分别表示：Calibration method、Name of agency、Number of individual antennas calibrated、Date
- 其中标定方法Calibration method有：‘CHAMBER’, ‘FIELD’, ‘ROBOT’,COPIED’ from other antenna, ‘CONVERTED’ from igs_01.pcv or blank
-
DAZI表示方位角步长，0.0表示与方位角无关
图2中为5，表示给出以5°为间隔0-360°随方位角变化的天线相位改正，对应数据在NORTH / EAST / UP行下，第一行NOAZI表示不顾及方位角时PCV，实际上就是方位角0-360°时的PCV平均值，
-
ZEN1 / ZEN2 / DZEN
- 对于接收机和天线均为：天顶角起始、终止以及增量
- 图1 中为天顶角从0-17°，每1°一个数据，一共有18个数据，对应下边在NORTH / EAST / UP行下一行
- 图2 中为天顶角从0-90°，每5°一个数据，一共有19列数据，方位角从0-360°，每5°一个数据，是73行，加上NOAZI一共是74行；可采用双线性内插法获取任意高度角的PCV
-
OF FREQUENCIES表示该条记录中的频率数
-
VALID FROM VALID UNTIL 为记录的有效时间
-
START OF FREQUENCY 和END OF FREQUENC标识对应频率的记录
-
NORTH / EAST / UP（PCO）
-
对接收机：平均相位中心与接收机天线参考点NEU方向的偏差，单位为毫米
为测站地平坐标系的三个分量，可转化为XYZ。
-
对卫星：平均相位中心与卫星天线质心在XYZ方向的偏差，单位为毫米
-
4. GAMP中相关函数：
4.1 读取atx文件：（RTKLIB/GAMP）
static int readantex(const char *file, pcvs_t *pcvs)
{
FILE *fp;
static const pcv_t pcv0={0};
pcv_t pcv;
double neu[3];
int i,f,freq=0,state=0,freqs[]={1,2,5,0};
char buff[256];
trace(3,"readantex: file=%s\n",file);
if (!(fp=fopen(file,"r"))) {
trace(2,"antex pcv file open error: %s\n",file);
return 0;
}
while (fgets(buff,sizeof(buff),fp)) {
if (strlen(buff)<60||strstr(buff+60,"COMMENT")) continue;
if (strstr(buff+60,"START OF ANTENNA")) {
pcv=pcv0;
state=1;
}
...
else if (strstr(buff,"NOAZI")) {
if (freq<1||NFREQ<freq) continue;
if ((i=decodef(buff+8,19,pcv.var[freq-1]))<=0) continue;
for (;i<19;i++) pcv.var[freq-1][i]=pcv.var[freq-1][i-1];
}
}
fclose(fp);
return 1;
}
- 读取的天线文件内容存储在pcvs_t结构体中，该结构体为：
typedef struct { /* antenna parameters type */
int n,nmax; /* number of data/allocated */
pcv_t *pcv; /* antenna parameters data */
} pcvs_t;
其中，pcv_t结构体对应为：
typedef struct { /* antenna parameter type */
int sat; /* satellite number (0:receiver) */
char type[MAXANT]; /* antenna type */
char code[MAXANT]; /* serial number or satellite code */
gtime_t ts,te; /* valid time start and end */
double off[NFREQ][ 3]; /* phase center offset e/n/u or x/y/z (m) */
double var[NFREQ][19]; /* phase center variation (m) */
/* el=90,85,...,0 or nadir=0,1,2,3,... (deg) */
} pcv_t;
-
RTKLIB中没有顾及PCV随方位角变化的情况
-
GAMP中相应函数与RTKLIB中基本一致，下面仅对比不一致之处。在readantex()函数中，GAMP添加了考虑PCV随方位角变化的部分，该部分对应代码为：
/* read antex file ----------------------------------------------------------*/
static int readantex(const char *file, pcvs_t *pcvs)
{
...
//该函数开始声明部分需要修改char buff[1024]否则报错
...
//接收机PCV考虑随方位角的变化
else if (strstr(buff,"NOAZI")) {
//if (SYS_CMP==PPP_Glo.sFlag[pcv.sat-1].sys) continue;
//if (freq<1||NFREQ<freq) continue;
if (freq<1) continue;
dd=(pcv.zen2-pcv.zen1)/pcv.dzen+1;
if (dd!=myRound(dd)||dd<=1) {
printf("*** WARNING: zen in atx file error (d!=round(d)||d<1)!\n");
continue;
}
if (pcv.dazi==0.0) {
//pcv.var[freq-1]=new double[int(dd)];
i=decodef(buff+8,(int)dd,pcv.var[freq-1]);//解码天线参数，调用decodef函数，形参中的v赋给了pcv.var
if (i<=0) {
printf("*** ERROR: error in reading atx (i<=0)!\n");
continue;
}
else if (i!=(int)dd) {
printf("*** ERROR: error in reading atx (i!=(int)dd)!\n");
continue;
}
}
else {
id=(int)((360-0)/pcv.dazi)+1;
//pcv.var[freq-1]=new double[int(dd)*id];
for (i=0;i<id;i++) {
fgets(buff,sizeof(buff),fp);
j=decodef(buff+8,(int)dd,&pcv.var[freq-1][i*(int)dd]);
if (j<=0) {
printf("*** ERROR: error in reading atx (j<=0)!\n");
continue;
}
else if (j!=(int)dd) {
printf("*** ERROR: error in reading atx (j!=(int)dd)!\n");
continue;
}
}
}
}
- readantex()函数运行完之后，读取的数据保存在pcvs_t结构体中，此处与RTKLIB不同的，是GAMP中pcv_t结构体加了随方位角变化的记录（增加了dazi、zen1,zen2,dzen，对应var大小变成了[80*20]）：
typedef struct { /* antenna parameter type */
int sat; /* satellite number (0:receiver) */
char type[MAXANT]; /* antenna type */
char code[MAXANT]; /* serial number or satellite code */
gtime_t ts,te; /* valid time start and end */
double off[NSYS_USED*NFREQ][3]; /* phase center offset e/n/u or x/y/z (m) */
double var[NSYS_USED*NFREQ][80*20]; /* phase center variation (m) */
/* el=90,85,...,0 or nadir=0,1,2,3,... (deg) */
double dazi; //Increment of the azimuth：0 to 360 with increment 'DAZI'(in degrees).
double zen1,zen2,dzen; //Receiver antenna:Definition of the grid in zenith angle.
//Satellite antenna:Definition of the grid in nadir angle.
} pcv_t;
4.2 天线相位改正
在pppos()函数中，通过ppp_res()函数计算中对接收机和卫星天线相位中心进行了改正：
/* 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);
4.2.1 接收机天线相位改正
extern void antmodel(int sat, 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,ii=0,sys;
e[0]=sin(azel[0])*cosel;
e[1]=cos(azel[0])*cosel;
e[2]=sin(azel[1]);
sys=PPP_Glo.sFlag[sat-1].sys;
if(strlen(pcv->type)==0) {
for (i=0;i<NFREQ;i++) {
if (sys==SYS_GPS||sys==SYS_CMP||sys==SYS_GAL||sys==SYS_QZS) {
ii = i;
if (i==2) ii=1;
}
else if (sys==SYS_GLO) {
ii = i+NFREQ;
if (i==2) ii=1+NFREQ;
}
for (j=0;j<3;j++) off[j]=pcv->off[ii][j]+del[j];
if (norm(pcv->off[ii],3)>0.0) {
sprintf(PPP_Glo.chMsg,"norm(pcv->off[ii],3)>0.0\n");
outDebug(OUTWIN,OUTFIL,0);
}
dant[i]=-dot(off,e,3)+(opt?interpvar0(0,90.0-azel[1]*R2D,pcv->var[ii],0):0.0);
}
return ;
}
/* 考虑方位角因素 */
for (i=0;i<NFREQ;i++) {
if (sys==SYS_GPS||sys==SYS_CMP||sys==SYS_GAL||sys==SYS_QZS) {
ii = i;
if (i==2) ii=1;
}
else if (sys==SYS_GLO) {
ii = i+NFREQ;
if (i==2) ii=1+NFREQ;
}
for (j=0;j<3;j++) off[j]=pcv->off[ii][j]+del[j];//将从文件头中读取的天线高加到接收机天线PCO上
if (pcv->dazi != 0.0) {
dant[i] = -dot(off,e,3)+interpvar1(azel[0]*R2D, 90-azel[1]*R2D, pcv, ii);//顾及方位角的PCV
}
else
dant[i]=-dot(off,e,3)+interpvar0(0,90.0-azel[1]*R2D,pcv->var[ii],0);
此处接收机的PCO和PCV均做了改正，-dot(off,e,3)为PCO，后面通过插值计算PCV，当dazi≠0时顾及方位角，调用interpvar1进行插值，当dazi=0时调用interpvar0进行插值。
4.2.2 卫星天线相位改正
- 卫星PCV改正：
/* 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]);
}
-
卫星PCO改正：（在利用精密星历计算卫星坐标时）
计算步骤大致为：
- 计算太阳位置ecef坐标
- 星固系到ecef坐标转换的旋转矩阵
- 计算无电离层组合系数
- 计算无电离层组合下的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,sys;
/* 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);
sys=satsys(sat,NULL);
//if (NFREQ>=3&&(sys&(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;
}
}
计算出的dant回到peph2pos()函数中，加到卫星位置上：
for (i=0;i<3;i++) {
rs[i ]=rss[i]+dant[i];
rs[i+3]=(rst[i]-rss[i])/tt;
}
参考资料：
- 101-RTKLIB中的天线相位中心改正
- 我的rtklib不完全学习笔记之四 .atx文件
- 《GPS测量与数据处理》（第三版）李征航黄劲松
扫一扫
2万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
