原始 Markdown文档、Visio流程图、XMind思维导图见：https://github.com/LiZhengXiao99/Navigation-Learning
文章目录
- 一、卫星天线相位中心改正
- 1、原理
- 2、文件读取
- 3、setpcv()：设置天线参数
- 4、satantoff()：卫星 PCO 改正
- 5、satantpcv()：卫星 PCV 改正
- 二、接收机天线相位中心改正
- 1、原理
- 2、antmodel()：接收机 PCO、PCV 改正
- 三、天线相位缠绕改正
- 1、model_phw()：计算天线相位缠绕改正
- 2、sat_yaw()：卫星姿态
- 3、sunmoonpos()：计算 ECEF 下日月坐标
- 4、sunmoonpos_eci()：计算 ECI 下日月坐标
- 4、yaw_angle()、yaw_nominal()：计算卫星名义姿态航偏角
- 四、潮汐改正
- 1、原理
- 2、tidedisp()：潮汐改正入口函数
- 3、获取 ERP 参数
- 1. ERP 参数介绍
- 2. readerp()：读取 ERP 参数
- 3. geterp()：插值获取当前时间的 ERP 参数
- 4、tide_solid()：计算固体潮
- 5、tide_oload()：计算海潮负荷
- 6、tide_pole()：计算极潮
- 五、地球自转效应改正
- 六、引力延迟改正
- 1、原理
- 2、gravitationalDelayCorrection()：引力延迟改正
一、卫星天线相位中心改正
1、原理
GNSS 的距离测量值为接收机天线至卫星天线的几何距离，而一般精密产品给出的卫星的坐标以卫星质量中心为参考，我们把相位中心和质量中心之间的差异称为卫星天线相位误差。实际测量中，天线相位引起的误差随时间变化，在对其处理时，我们一般将其分为常量和变量两个部分。常量部分称为卫星天线相位中心偏差（Phase Center Offset，PCO），表示卫星质量中心和卫星平均相位中心的差异。变量部分称为相位中心变化（Phase Center Variation，PCV），表示天线瞬时相位中心和平均相位中心之间的差异。从 2006 年 11 月 5 日，IGS 开始使用绝对相位中心改正模型 IGS_05，该模型给出了与天底角相关的卫星 PCV 和不同型号接收机的 PCO。tp://sopac-ftp.ucsd.edu/archive/garner/ gamit/tables/ 可以下载最新的 IGS 天线改正文件，其包括最新 GNSS 卫星和接收机的天线 PCO 和 PCV 改正信息。一般通过 igs14.atx 文件链接到最新的天线文件。
2、文件读取
听说 RTKLIB 的 readantex() 函数有bug，接收机端同时出现 GPS、GLO 的PCO、PCV时，会用 GLO 系统的值覆盖 GPS，不知道 GAMP改了没有。
3、setpcv()：设置天线参数
查找各颗卫星的天线参数，并放在 nav->pcvs 中；查找接收机的天线参数，并放在 popt->pcvr 中，需要注意的是，只有当 popt->anttype[i] 是 *(auto) 的时候，通过 RINEX 头识别到的接收机天线类型，以及天线的 H/E/N 偏差才会起作用。
4、satantoff()：卫星 PCO 改正
分别计算 L1 和 L2 的卫星端 PCO。如果 PPP 中使用的是无电离层组合，因此会计算无电离层组合后的 PCO 修正值。之后在 peph2pos() 函数中，会在由精密星历计算的卫星位置上进行 PCO 修正。
e
z
s
=
−
r
s
∣
r
s
∣
e
s
=
r
sun
−
r
s
∣
r
sun
−
r
s
∣
e
y
s
=
e
z
s
×
e
s
∣
e
z
s
×
e
s
∣
e
x
s
=
e
y
s
×
e
z
s
E
s
=
(
e
x
s
,
e
y
s
,
e
z
s
)
\begin{array}{l}\boldsymbol{e}_{z}^{s}=-\frac{\boldsymbol{r}^{s}}{\left|\boldsymbol{r}^{s}\right|} \\ \boldsymbol{e}_{s}=\frac{\boldsymbol{r}_{\text {sun }}-\boldsymbol{r}^{s}}{\left|\boldsymbol{r}_{\text {sun }}-\boldsymbol{r}^{s}\right|} \\ \boldsymbol{e}_{y}^{s}=\frac{\boldsymbol{e}_{z}^{s} \times \boldsymbol{e}_{s}}{\left|\boldsymbol{e}_{z}^{s} \times \boldsymbol{e}_{s}\right|} \\ \boldsymbol{e}_{x}^{s}=\boldsymbol{e}_{y}^{s} \times \boldsymbol{e}_{z}^{s} \\ \boldsymbol{E}_{s}=\left(\boldsymbol{e}_{x}^{s}, \boldsymbol{e}_{y}^{s}, \boldsymbol{e}_{z}^{s}\right)\end{array}
ezs=−∣rs∣rses=∣rsun −rs∣rsun −rseys=∣ezs×es∣ezs×esexs=eys×ezsEs=(exs,eys,ezs)
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 */
for (i=0;i<3;i++) r[i]=-rs[i];
if (!normv3(r,ez)) return; // (E.8.5)
for (i=0;i<3;i++) r[i]=rsun[i]-rs[i];
if (!normv3(r,es)) return; // (E.8.6)
cross3(ez,es,r);
if (!normv3(r,ey)) return; // (E.8.7)
cross3(ey,ez,ex); // (E.8.8)
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;
// 把 L1 频率转到 L2
gamma=SQR(lam[k])/SQR(lam[j]);
C1=gamma/(gamma-1.0);
C2=-1.0 /(gamma-1.0);
if (sys==SYS_GPS) {
j=0;
k=1;
}
else if (sys==SYS_GLO) {
j=0+NFREQ;
k=1+NFREQ;
}
else if (sys==SYS_CMP) {
j=0+2*NFREQ;
k=1+2*NFREQ;
}
else if (sys==SYS_GAL) {
j=0+3*NFREQ;
k=1+3*NFREQ;
}
else if (sys==SYS_QZS) {
j=0+4*NFREQ;
k=1+4*NFREQ;
}
/* 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;
}
}
5、satantpcv()：卫星 PCV 改正
θ = arccos e r s T r s ∣ r s ∣ \theta=\arccos \frac{\boldsymbol{e}_{r}^{s T} \boldsymbol{r}^{s}}{\left|\boldsymbol{r}^{s}\right|} θ=arccos∣rs∣ersTrs
static void satantpcv(int sat, 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); // (E.8.10)
antmodel_s(sat,pcv,nadir,dant);
}
二、接收机天线相位中心改正
1、原理
和卫星相似，接收机端也存在由天线相位中心引起的误差 PCO 和 PCV。GNSS观测量是相对于接收机天线的平均相位中心而言的，而接收机天线对中是相对于几何中也而言的，这两种中心一般不重合，两者之差就称为平均相位中心偏差（PCO），其大小可达毫米级或厘米级。且接收机天线的相位中也会随卫星信号输入的方向和强度的变化而变化，此时观测时刻的瞬时相位中也与平均相位中心的偏差称为平均相位中心变化（PCV），它与卫星的高度角和方位角有关。因此接收机天线相位偏差由接收机天线PCO和PCV两部分组成。
对于常见型号的接收机，IGS 给出了 GPS 和 GLONASS 的改正信息，可从天线改正文件中获取。由于缺少 BDS 的接收机天线相位中心改正值，在 PPP 数据处理中，一般将 GPS 的接收机 PCO 和 PCV 信息用于 BDS。
NGS 提供的 ANTEX 格式天线模型，包含了卫星天线模型以及部分接收机天线模型。使用天线模型的目的包括：
-
修正天线参考点和天线相位中心的之间的偏差。
-
修正和仰角有关的误差。
-
修正 L1 和 L2 之间的相位中心偏差（这个误差可能对整周模糊度固定造成影响）。
-
RTKLIB 支持 NGS PCV 以及 ANTEX 格式的天线模型，其中包括了 PCO 和 PCV 修正参数。通过手册 E.8 章节可知，接收机天线修正如下：
-
PCO修正通常是当地坐标系 ENU 参数，因此需要利用转换矩阵转到 ECEF 坐标系:
d r , p c o , i = E r T d r , p c o , i , e n u \boldsymbol{d}_{r, p c o, i}=\boldsymbol{E}_{r}{ }^{T} \boldsymbol{d}_{r, p c o, i, e n u} dr,pco,i=ErTdr,pco,i,enu -
PCV修正则通过对高度角进行插值得到：
d r , p c o , i = E r T d r , p c o , i , e n u d r , p c v , i ( E l ) = ( E l − E l i ) d r , p c v , i ( E l i ) + ( E l i + 1 − E l ) d r , p c v , i ( E l i + 1 ) E l i + 1 − E l i \boldsymbol{d}_{r, p c o, i}=\boldsymbol{E}_{r}{ }^{T} \boldsymbol{d}_{r, p c o, i, e n u}d_{r, p c v, i}(E l)=\frac{\left(E l-E l_{i}\right) d_{r, p c v, i}\left(E l_{i}\right)+\left(E l_{i+1}-E l\right) d_{r, p c v, i}\left(E l_{i+1}\right)}{E l_{i+1}-E l_{i}} dr,pco,i=ErTdr,pco,i,enudr,pcv,i(El)=Eli+1−Eli(El−Eli)dr,pcv,i(Eli)+(Eli+1−El)dr,pcv,i(Eli+1)
-
-
2、antmodel()：接收机 PCO、PCV 改正
- 计算 LOS 视向量在 ENU 中的单位矢量
e - 频段不同，天线的相位中心偏移(PCO)不同，先计算出每个频段天线在东、北、天三个方向总的偏移,即相位中心偏移
pcv->off与del[j]之和。 - 计算相位中心偏移(PCO)在观测单位矢量
e上的投影dot(off,e,3) - 计算天线相位中心变化量(PCV)：不同的高度角，相位中心变化不同，因此根据高度角对
pcv->var[i]进行插值计算。 - PCO 和 PCV 两部分求和为
dant[]
del为相对天线参考点偏移值
azel为方位角和俯仰角，
pcv->off为 phase center offset（PCO）
pcv->var为 phase center variations （PCV）
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;
// 计算视线向量在 ENU 中的单位矢量 e
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;
}
// 相位中心偏移(PCO)，pcv->off[i][j] 中的值来自于天线 PCV 文件
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);
}
// 相位中心偏移(PCO)在观测单位矢量 e 上的投影 dot(off,e,3)
// 计算天线相位中心变化量(PCV)：不同的高度角，相位中心变化不同，因此根据高度角对 pcv->var[i] 进行插值计算。
// dant[] 为上面两部分相加
dant[i]=-dot(off,e,3)+(opt?interpvar0(0,90.0-azel[1]*R2D,pcv->var[ii],0):0.0);
}
return ;
}
// 频段不同，天线的相位中心偏移(PCO)不同。
// 先计算出每个频段天线在东、北、天三个方向总的偏移，即相位中心偏移 pcv->off[i][j] 与 del[j] 之和
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;
}
// 相位中心偏移(PCO)，pcv->off[i][j] 中的值来自于天线 PCV 文件
for (j=0;j<3;j++) off[j]=pcv->off[ii][j]+del[j];
if (pcv->dazi!=0.0)
dant[i]=-dot(off,e,3)+interpvar1(azel[0]*R2D,90-azel[1]*R2D,pcv,ii);
else
dant[i]=-dot(off,e,3)+interpvar0(0,90.0-azel[1]*R2D,pcv->var[ii],0);
}
}
三、天线相位缠绕改正
GNSS 载波相位是右旋圆极化电磁波，当接收机天线和卫星天线发生相对旋转时，载波相位观测值会因此产生误差，即相位缠绕（phase wind-up）。相位缠绕误差最大可达载波相位的一个波长，需要进行改正。
这部分算法模型摘自：GPS从入门到放弃（二十三） — 相位缠绕
如下图所示是卫星、地球与太阳的位置关系：
在卫星天线上建立卫星天线坐标系，以卫星的天线相位中心为原点；
Z
Z
Z 轴沿卫星天线方向指向地心；
X
X
X 轴在地球、太阳和卫星组成的平面内，指向太阳；Y轴与Z轴和X轴构成右手系。于是可以求得三轴方向的单位矢量
e
x
s
,
e
y
s
,
e
z
s
e_{x s}, e_{y s}, e_{z s}
exs,eys,ezs 分别为：
e
z
s
=
−
r
s
a
t
∣
r
sat
∣
e
y
s
=
e
z
s
×
e
s
u
n
e
x
s
=
e
z
s
×
e
y
s
\begin{array}{c} e_{z s}=\frac{-\boldsymbol{r}_{s a t}}{\left|\boldsymbol{r}_{\text {sat }}\right|} \\ e_{y s}=e_{z s} \times \boldsymbol{e}_{s u n} \\ \boldsymbol{e}_{x s}=\boldsymbol{e}_{z s} \times \boldsymbol{e}_{y s} \end{array}
ezs=∣rsat ∣−rsateys=ezs×esunexs=ezs×eys
其中
r
s
a
t
\boldsymbol{r}_{\mathrm{sat}}
rsat 和
r
s
u
n
\boldsymbol{r}_{\mathrm{sun}}
rsun 分别是 ECEF 坐标系中卫星和太阳的位置矢量，而
e
s
u
n
\boldsymbol{e}_{\mathrm{sun}}
esun 是卫星至太阳方向的单位矢量：
e
sun
=
r
sun
−
r
sat
∣
r
sun
−
r
sat
∣
e_{\text {sun }}=\frac{\boldsymbol{r}_{\text {sun }}-\boldsymbol{r}_{\text {sat }}}{\left|\boldsymbol{r}_{\text {sun }}-\boldsymbol{r}_{\text {sat }}\right|}
esun =∣rsun −rsat ∣rsun −rsat
计算相位缠绕时，在卫星和接收机处各定义一个有效偶极
D
s
\boldsymbol{D}_{s}
Ds 和
D
r
\boldsymbol{D}_{\boldsymbol{r}}
Dr ，且分别对应于星固坐标系和接收机所在位置的站心坐标系，有：
D
s
=
e
x
s
−
e
k
(
e
k
⋅
e
x
s
)
−
e
k
×
e
y
s
D
r
=
e
x
r
−
e
k
(
e
k
⋅
e
x
r
)
+
e
k
×
e
y
r
ϕ
f
=
sign
(
e
k
⋅
(
D
s
×
D
r
)
)
arccos
(
D
s
⋅
D
r
∣
D
s
∣
∣
D
r
∣
)
\begin{array}{c} \boldsymbol{D}_{s}=e_{x s}-e_{k}\left(e_{k} \cdot e_{x s}\right)-e_{k} \times e_{y s} \\ \boldsymbol{D}_{r}=e_{x r}-e_{k}\left(e_{k} \cdot e_{x r}\right)+e_{k} \times e_{y r} \\ \phi_{\mathrm{f}}=\operatorname{sign}\left(e_{k} \cdot\left(\boldsymbol{D}_{s} \times \boldsymbol{D}_{r}\right)\right) \arccos \left(\frac{\boldsymbol{D}_{s} \cdot \boldsymbol{D}_{r}}{\left|\boldsymbol{D}_{s}\right|\left|\boldsymbol{D}_{r}\right|}\right) \end{array}
Ds=exs−ek(ek⋅exs)−ek×eysDr=exr−ek(ek⋅exr)+ek×eyrϕf=sign(ek⋅(Ds×Dr))arccos(∣Ds∣∣Dr∣Ds⋅Dr)
其中
e
k
e_{k}
ek 为卫星至接收机方向的单位向量，
e
x
r
、
e
y
r
e_{x r} 、 e_{y r}
exr、eyr 为接收机所在位置的站心坐标系的坐标轴方向的单位矢量，而
ϕ
f
\phi_{\mathrm{f}}
ϕf 为相位缠绕值的小数部分。
e
k
=
r
r
−
r
s
a
t
∣
r
r
−
r
s
a
t
∣
e_{k}=\frac{\boldsymbol{r}_{r}-\boldsymbol{r}_{s a t}}{\left|\boldsymbol{r}_{\boldsymbol{r}}-\boldsymbol{r}_{s a t}\right|}
ek=∣rr−rsat∣rr−rsat
总结可以得出完整公式如下：
E
r
=
(
e
r
,
x
T
,
e
r
,
y
T
,
e
r
,
z
T
)
T
E
s
=
(
e
x
s
T
,
e
y
s
T
,
e
z
s
T
)
T
D
s
=
e
x
s
−
e
u
s
(
e
u
s
⋅
e
x
s
)
−
e
u
s
×
e
y
s
D
r
=
e
r
,
x
−
e
r
s
(
e
r
s
⋅
e
r
,
x
)
+
e
r
s
×
e
r
,
y
ϕ
p
w
=
sign
(
e
r
s
⋅
(
D
s
×
D
r
)
)
arccos
D
s
⋅
D
r
∥
D
s
∥
∥
D
r
∥
/
2
π
+
N
\begin{array}{l}\boldsymbol{E}_{r}=\left(\boldsymbol{e}_{r, x}{ }^{T}, \boldsymbol{e}_{r, y}{ }^{T}, \boldsymbol{e}_{r, z}{ }^{T}\right)^{T} \\ \boldsymbol{E}^{s}=\left(\boldsymbol{e}_{x}^{s T}, \boldsymbol{e}_{y}^{s T}, \boldsymbol{e}_{z}^{s T}\right)^{T} \\ \boldsymbol{D}^{s}=\boldsymbol{e}_{x}^{s}-\boldsymbol{e}_{u}^{s}\left(\boldsymbol{e}_{u}^{s} \cdot \boldsymbol{e}_{x}^{s}\right)-\boldsymbol{e}_{u}^{s} \times \boldsymbol{e}_{y}^{s} \\ \boldsymbol{D}_{r}=\boldsymbol{e}_{r, x}-\boldsymbol{e}_{r}^{s}\left(\boldsymbol{e}_{r}^{s} \cdot \boldsymbol{e}_{r, x}\right)+\boldsymbol{e}_{r}^{s} \times \boldsymbol{e}_{r, y} \\ \phi_{p w}=\operatorname{sign}\left(\boldsymbol{e}_{r}^{s} \cdot\left(\boldsymbol{D}^{s} \times \boldsymbol{D}_{r}\right)\right) \arccos \frac{\boldsymbol{D}^{s} \cdot \boldsymbol{D}_{r}}{\left\|\boldsymbol{D}^{s}\right\|\left\|\boldsymbol{D}_{r}\right\|} / 2 \pi+N\end{array}
Er=(er,xT,er,yT,er,zT)TEs=(exsT,eysT,ezsT)TDs=exs−eus(eus⋅exs)−eus×eysDr=er,x−ers(ers⋅er,x)+ers×er,yϕpw=sign(ers⋅(Ds×Dr))arccos∥Ds∥∥Dr∥Ds⋅Dr/2π+N
式中，
e
r
s
\mathbf{e}_{r}^{s}
ers 表示卫星指向接收机的单位向量；
x
,
y
\mathbf{x}, \mathbf{y}
x,y 和
x
′
,
y
′
\mathbf{x}^{\prime}, \mathbf{y}^{\prime}
x′,y′ 分别表示接收机和卫星的两个有效偶极矢量；sign 表示符号函数。
1、model_phw()：计算天线相位缠绕改正
static int model_phw(gtime_t time, int sat, const char *type, int opt,
const double *rs, const double *rr, double *phw)
{
double exs[3],eys[3],ek[3],exr[3],eyr[3],eks[3],ekr[3],E[9];
double dr[3],ds[3],drs[3],r[3],pos[3],cosp,ph;
int i;
if (opt<=0) return 1; /* no phase windup */
// 首先调用 sat-yaw 函数，根据卫星的姿态模型计算出卫星本体坐标系 X,Y 方向的单位矢量exs、eys，即上面公式里的SX、SY
/* satellite yaw attitude model */
if (!sat_yaw(time,sat,type,opt,rs,exs,eys)) return 0;
// 计算卫星至接收机的单位矢量
/* unit vector satellite to receiver */
for (i=0;i<3;i++) r[i]=rr[i]-rs[i];
if (!normv3(r,ek)) return 0;
// 计算接收机天线在当地坐标系的北向、西向单位矢量
/* unit vectors of receiver antenna */
ecef2pos(rr,pos);
xyz2enu(pos,E);
exr[0]= E[1]; exr[1]= E[4]; exr[2]= E[7]; /* x = north */
eyr[0]=-E[0]; eyr[1]=-E[3]; eyr[2]=-E[6]; /* y = west */
// 根据公式以及前一次的相位缠绕误差计算当前时刻相位缠绕误差
/* phase windup effect */
cross3(ek,eys,eks);
cross3(ek,eyr,ekr);
for (i=0;i<3;i++) {
ds[i]=exs[i]-ek[i]*dot(ek,exs,3)-eks[i];
dr[i]=exr[i]-ek[i]*dot(ek,exr,3)+ekr[i];
}
cosp=dot(ds,dr,3)/norm(ds,3)/norm(dr,3);
if (cosp<-1.0) cosp=-1.0;
else if (cosp> 1.0) cosp= 1.0;
ph=acos(cosp)/2.0/PI;
cross3(ds,dr,drs);
if (dot(ek,drs,3)<0.0) ph=-ph;
*phw=ph+floor(*phw-ph+0.5); /* in cycle */
return 1;
}
2、sat_yaw()：卫星姿态
由于不同类型卫星制造商的星固坐标系定义不同，为保持一致性，IGS 定义星固系如下
- Z 轴平行于卫星天线信号发射方向并指向地心。
- Y 轴平行于太阳帆板并垂直于太阳、地球和卫星构成的平面。
- X 轴垂直于 Y 轴和 Z 轴并构成右手坐标系并指向太阳入射方向。
GNSS 卫星名义姿态在星固系下 3 轴单位向量
e
x
、
e
y
、
e
z
\boldsymbol{e}_{x} 、 \boldsymbol{e}_{y} 、 \boldsymbol{e}_{z}
ex、ey、ez 可由下式确定：
e
x
=
e
y
×
e
z
e
y
=
e
⊗
×
r
∣
e
⊗
×
r
∣
e
z
=
−
r
∣
r
∣
}
\left.\begin{array}{l} e_{x}=e_{y} \times e_{z} \\ e_{y}=\frac{e_{\otimes} \times r}{\left|e_{\otimes} \times r\right|} \\ e_{z}=-\frac{r}{|r|} \end{array}\right\}
ex=ey×ezey=∣e⊗×r∣e⊗×rez=−∣r∣r⎭
⎬
⎫
式中，
e
⊗
e_{\otimes}
e⊗ 为卫星至太阳方向的单位向量；
r
r
r 为地心指向卫星方向的单位向量；
∣
∗
∣
\mid*\mid
∣∗∣ 表示向量取模运算符。GNSS卫星偏航角
φ
\varphi
φ 定义为沿轨道切线方向与星固系
X
X
X 轴之间的夹角：
φ
=
arccos
(
e
T
⋅
e
x
)
\varphi=\arccos \left(e_{T} \cdot e_{x}\right)
φ=arccos(eT⋅ex)
式中，
e
T
、
e
X
\boldsymbol{e}_{T} 、 \boldsymbol{e}_{X}
eT、eX 分别沿轨道切线方向、卫星星固系 X 轴单位向量；
arccos
(
⋅
)
\arccos (\cdot)
arccos(⋅) 为反余弦函数。根据太阳高度角、轨道角与下式的几何关系，名义姿态偏航角可以表示为：
φ
=
arctan
2
(
−
tan
β
,
sin
μ
)
\varphi=\arctan 2(-\tan \beta, \sin \mu)
φ=arctan2(−tanβ,sinμ)
式中，
β
\beta
β 为太阳高度角；
μ
\mu
μ 为轨道角 (以远日点为起点)。对 GNSS 卫星而言，由于卫星的信号发射方向始终指向地心，因此不存在俯仰角与横滚角，卫星姿态仅用偏航姿态角
φ
\varphi
φ 确定 ，如图所示，将卫星在轨切线方向
e
T
\boldsymbol{e}_{T}
eT 绕星固系的
Z
Z
Z 轴旋转
φ
\varphi
φ 角度，即可确定星固系
X
X
X 轴的指向，因此，卫星在姿态异常时期，偏航姿态模型的建立主要是确定偏航姿态角
φ
\varphi
φ 的变化。
代码中：
- 调用
sunmoonpos()计算日月 ECEF 坐标rs、rm。 - 计算太阳高度角
beta、轨道角mu。 - 调用
yaw_angle()计算名义姿态偏航角yaw。 - 计算卫星名义姿态在星固系下
exs、eys。
static int sat_yaw(gtime_t time, int sat, const char *type, int opt,
const double *rs, double *exs, double *eys)
{
double rsun[3],ri[6],es[3],esun[3],n[3],p[3],en[3],ep[3],ex[3],E,beta,mu;
double yaw,cosy,siny,erpv[5]={0};
int i;
// 调用 sunmoonpos() 计算日月 ECEF 坐标 rs、rm
sunmoonpos(gpst2utc(time),erpv,rsun,NULL,NULL);
// 计算太阳高度角 beta、轨道角 mu
/* beta and orbit angle */
matcpy(ri,rs,6,1);
ri[3]-=OMGE*ri[1];
ri[4]+=OMGE*ri[0];
cross3(ri,ri+3,n);
cross3(rsun,n,p);
if (!normv3(rs,es)||!normv3(rsun,esun)||!normv3(n,en)||
!normv3(p,ep)) return 0;
beta=PI/2.0-acos(dot(esun,en,3));
E=acos(dot(es,ep,3));
mu=PI/2.0+(dot(es,esun,3)<=0?-E:E);
if (mu<-PI/2.0) mu+=2.0*PI;
else if (mu>=PI/2.0) mu-=2.0*PI;
// 调用 yaw_angle() 计算名义姿态偏航角 yaw
/* yaw-angle of satellite */
if (!yaw_angle(sat,type,opt,beta,mu,&yaw)) return 0;
// 计算卫星名义姿态在星固系下 exs、eys
/* satellite fixed x,y-vector */
cross3(en,es,ex);
cosy=cos(yaw);
siny=sin(yaw);
for (i=0;i<3;i++) {
exs[i]=-siny*en[i]+cosy*ex[i];
eys[i]=-cosy*en[i]-siny*ex[i];
}
return 1;
}
3、sunmoonpos()：计算 ECEF 下日月坐标
- 根据 ERP 参数计算 UT1 时间
- 调用 sunmoonpos_eci() 计算日月 ECI 坐标
- 调用 eci2ecef() 计算 ECI 到 ECEF 的转换矩阵 U
- 用转换矩阵 U 将日月坐标转到 ECEF
extern void sunmoonpos(gtime_t tutc, const double *erpv, double *rsun,
double *rmoon, double *gmst)
{
gtime_t tut;
double rs[3],rm[3],U[9],gmst_;
// 根据 ERP 参数计算 UT1 时间
tut=timeadd(tutc,erpv[2]); /* utc -> ut1 */
// 调用 sunmoonpos_eci() 计算日月 ECI 坐标
/* sun and moon position in eci */
sunmoonpos_eci(tut,rsun?rs:NULL,rmoon?rm:NULL);
// 调用 eci2ecef() 计算 ECI 到 ECEF 的转换矩阵 U
/* eci to ecef transformation matrix */
eci2ecef(tutc,erpv,U,&gmst_);
// 用转换矩阵 U 将日月坐标转到 ECEF
/* sun and moon postion in ecef */
if (rsun ) matmul("NN",3,1,3,1.0,U,rs,0.0,rsun );
if (rmoon) matmul("NN",3,1,3,1.0,U,rm,0.0,rmoon);
if (gmst ) *gmst=gmst_;
}
4、sunmoonpos_eci()：计算 ECI 下日月坐标
太阳的位置是通过历书时、太阳的视运动、恒星际位置以及太阳辐射量等信息来计算的；而月亮的位置则是通过观察月亮的视运动等信息来计算的。
- 确定观测时间：通过输入的时间变量（如
tut）确定观测时间。 - 计算天文参数：通过函数
ast_args(t, f)等计算与观测时间相关的天文参数。 - 确定赤经和赤纬：根据观测时间和相关天文参数计算出太阳或月亮的赤经和赤纬。
- 计算ECI坐标：使用计算出的赤经和赤纬，结合ECI坐标系的定义，计算出太阳或月亮在ECI坐标系中的位置。
static void sunmoonpos_eci(gtime_t tut, double *rsun, double *rmoon)
{
const double ep2000[]={2000,1,1,12,0,0};
double t,f[5],eps,Ms,ls,rs,lm,pm,rm,sine,cose,sinp,cosp,sinl,cosl;
// 从2000年1月1日12时到输入时间当前
t=timediff(tut,epoch2time(ep2000))/86400.0/36525.0;
// 天文参数计算
/* astronomical arguments */
ast_args(t,f);
/* obliquity of the ecliptic */
eps=23.439291-0.0130042*t; // 黄赤交角
sine=sin(eps*D2R); cose=cos(eps*D2R);
/* sun position in eci */
if (rsun) {
Ms=357.5277233+35999.05034*t;
ls=280.460+36000.770*t+1.914666471*sin(Ms*D2R)+0.019994643*sin(2.0*Ms*D2R);
rs=AU*(1.000140612-0.016708617*cos(Ms*D2R)-0.000139589*cos(2.0*Ms*D2R));
sinl=sin(ls*D2R); cosl=cos(ls*D2R);
rsun[0]=rs*cosl;
rsun[1]=rs*cose*sinl;
rsun[2]=rs*sine*sinl;
}
/* moon position in eci */
if (rmoon) {
lm=218.32+481267.883*t+6.29*sin(f[0])-1.27*sin(f[0]-2.0*f[3])+
0.66*sin(2.0*f[3])+0.21*sin(2.0*f[0])-0.19*sin(f[1])-0.11*sin(2.0*f[2]);
pm=5.13*sin(f[2])+0.28*sin(f[0]+f[2])-0.28*sin(f[2]-f[0])-
0.17*sin(f[2]-2.0*f[3]);
rm=RE_WGS84/sin((0.9508+0.0518*cos(f[0])+0.0095*cos(f[0]-2.0*f[3])+
0.0078*cos(2.0*f[3])+0.0028*cos(2.0*f[0]))*D2R);
sinl=sin(lm*D2R); cosl=cos(lm*D2R);
sinp=sin(pm*D2R); cosp=cos(pm*D2R);
rmoon[0]=rm*cosp*cosl;
rmoon[1]=rm*(cose*cosp*sinl-sine*sinp);
rmoon[2]=rm*(sine*cosp*sinl+cose*sinp);
}
}
4、yaw_angle()、yaw_nominal()：计算卫星名义姿态航偏角
yaw_angle() 就是直接调用 yaw_nominal() 用下面公式计算：
φ
=
arctan
2
(
−
tan
β
,
sin
μ
)
+
π
\varphi=\arctan 2(-\tan \beta, \sin \mu)+\pi
φ=arctan2(−tanβ,sinμ)+π
extern int yaw_angle(int sat, const char *type, int opt, double beta, double mu,
double *yaw)
{
*yaw=yaw_nominal(beta,mu);
return 1;
}
static double yaw_nominal(double beta, double mu)
{
if (fabs(beta)<1E-12&&fabs(mu)<1E-12) return PI;
return atan2(-tan(beta),sin(mu))+PI;
}
四、潮汐改正
潮汐改正这部分，公式的具体参数与 RTKLIB 略有不同，但原理和计算方法一样。
1、原理
由于地球不是理想刚体，形状会因其他天体的引力发生形变，所以即使将接收机固定在地球表面，接收机和地心的相位位置也会发生变化，接收机在地固系中的坐标随之改变，由此生的测量误差称为地球潮汐效应，影响如下：
GNSS 数据处理一般将潮汐分为三个部分，地球固体潮、海洋负荷潮和极潮。
- 地球固体潮：是指地球固体在其他天体的引力作用下发生周期变化的现象，对接收机水平和高程方向的影响分别可达数厘米和数分米。
- 海洋负荷潮：是指在太阳和月球引力作用下地球海洋发生的周期性变化，对接收机水平和高程方向的影响都在厘米级。
- 极潮：是指地球自转轴发生的瞬时变化引起的测量误差，对接收机影响在厘米级。
关于以上潮汐改正，目前主要通过国际地球自转协议（International Earth Rotation Service，IERS）的 IERS Convention 技术文档中的模型改正。
2、tidedisp()：潮汐改正入口函数
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;
// 如果有 erp 参数，调用 geterp() 获取
if (erp) {
geterp(erp,utc2gpst(tutc),erpv);
}
// UTC 加上 erpv[2] 得到 UT 时间 tut
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 */
// 调用 sunmoonpos() 计算日月 ECEF 坐标 rs、rm
/* sun and moon position in ecef */
sunmoonpos(tutc,erpv,rs,rm,&gmst);
// 调用 tide_solid() 计算固体潮改正 drt，改正到 dr
tide_solid(rs,rm,pos,E,gmst,opt,drt);
for (i=0;i<3;i++) dr[i]+=drt[i];
}
if ((opt&2)&&odisp) { /* ocean tide loading */
// 调用 tide_oload() 计算海洋潮改正 drt，改正到 dr
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() 极潮改正到 drt，改正到 dr
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];
}
}
3、获取 ERP 参数
1. ERP 参数介绍
内容摘自博客：GPS从入门到放弃（二十一）地球自转参数
地球自转参数（ERP: Earth rotation parameters）主要包括地球极点的位移和速率、UT1-UTC的时间差、以及由天文观测确定的一天的时间长度与 86400 秒之间的差值 LOD。
地球自转参数可以从ftp服务站 ftp://cddis.nasa.gov/gnss/products/ 下载。IGS 提供的 ERP 数据与精密星历数据放在一个目录中。此路径下的数据是以 GPS 周数（GPS Week）为目录名整理放置的。比如想找 2020 年元旦的 ERP 数据，经过计算知道那一天是GPS周第 2086 周，所以进入2086 目录下去下载相应数据。
**类似于 IGS 精密星历，ERP 参数文件也分为三种：**最终 ERP 参数（IGS Final，标识为 IGS）、快速ERP参数（IGS Rapid，标识为 IGR）、以及超快速ERP参数（IGS Ultra-Rapid，标识为 IGU）。其中超快速ERP参数又分为观测的部分和预测的部分。他们的延时、精度等指标如下表所示。在实际工作中，我们可以根据项目对时间及精度的要求，选取不同类型的文件来使用。
2. readerp()：读取 ERP 参数
ERP格式非常简单，几乎不用解释，一看就明白。这里附上2020年1月1日的快速ERP参数文件igr20863.erp如下：
文件中除了说明，就是一个表格，每一行表示一天的ERP数据。当然，上面这个文件中表格内容只有一行。第一列是时间，这个时间是用简化的儒略日来表示的，儒略日 58849.50 就是 2020 年 1 月 1 日。
extern int readerp(const char *file, erp_t *erp)
{
FILE *fp;
erpd_t *erp_data;
double v[14]={0};
char buff[256];
// 以读的方式打开文件
if (!(fp=fopen(file,"r"))) {
printf("erp file open error: file=%s\n",file);
return 0;
}
// 循环读取每一行
while (fgets(buff,sizeof(buff),fp)) {
// sscanf 格式化取值，文件头不符合格式，自然就被跳过
if (sscanf(buff,"%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf",
v,v+1,v+2,v+3,v+4,v+5,v+6,v+7,v+8,v+9,v+10,v+11,v+12,v+13)<5) {
continue;
}
// 如果 erp 参数量超出容量，就扩大一倍容量
if (erp->n>=erp->nmax) {
erp->nmax=erp->nmax<=0?128:erp->nmax*2;
erp_data=(erpd_t *)realloc(erp->data,sizeof(erpd_t)*erp->nmax);
if (!erp_data) { // 空间开辟失败，就退出不继续读取
free(erp->data); erp->data=NULL; erp->n=erp->nmax=0;
fclose(fp);
return 0;
}
erp->data=erp_data;
}
erp->data[erp->n].mjd=v[0];
erp->data[erp->n].xp=v[1]*1E-6*AS2R;
erp->data[erp->n].yp=v[2]*1E-6*AS2R;
erp->data[erp->n].ut1_utc=v[3]*1E-7;
erp->data[erp->n].lod=v[4]*1E-7;
erp->data[erp->n].xpr=v[12]*1E-6*AS2R;
erp->data[erp->n++].ypr=v[13]*1E-6*AS2R;
}
fclose(fp);
return 1;
}
3. geterp()：插值获取当前时间的 ERP 参数
- 计算当前 mjd。
- 若当前时间早于 ERP 参数中最早的时间，采用最早的时间来计算。
- 若当前时间晚于 ERP 参数中最晚的时间，采用最晚的时间来计算。
- 若当前时间在 ERP 参数中最早与最晚的时间之间，则先找到最接近的两个时间，然后用插值。
extern int geterp(const erp_t *erp, gtime_t time, double *erpv)
{
const double ep[]={2000,1,1,12,0,0};
double mjd,day,a;
int i=0,j,k;
// 如果没有 ERP 参数直接返回
if (erp->n<=0) return 0;
// 计算当前 mjd
mjd=51544.5+(timediff(gpst2utc(time),epoch2time(ep)))/86400.0;
// 若当前时间早于 ERP 参数中最早的时间，采用最早的时间来计算
if (mjd<=erp->data[0].mjd) {
day=mjd-erp->data[0].mjd;
erpv[0]=erp->data[0].xp +erp->data[0].xpr*day;
erpv[1]=erp->data[0].yp +erp->data[0].ypr*day;
erpv[2]=erp->data[0].ut1_utc-erp->data[0].lod*day;
erpv[3]=erp->data[0].lod;
return 1;
}
// 若当前时间晚于 ERP 参数中最晚的时间，采用最晚的时间来计算
if (mjd>=erp->data[erp->n-1].mjd) {
day=mjd-erp->data[erp->n-1].mjd;
erpv[0]=erp->data[erp->n-1].xp +erp->data[erp->n-1].xpr*day;
erpv[1]=erp->data[erp->n-1].yp +erp->data[erp->n-1].ypr*day;
erpv[2]=erp->data[erp->n-1].ut1_utc-erp->data[erp->n-1].lod*day;
erpv[3]=erp->data[erp->n-1].lod;
return 1;
}
// 若当前时间在 ERP 参数中最早与最晚的时间之间，则先找到最接近的两个时间，然后用插值。
for (j=0,k=erp->n-1;j<=k;) {
i=(j+k)/2;
if (mjd<erp->data[i].mjd) k=i-1;
else if (mjd>erp->data[i+1].mjd) j=i+1;
else break;
}
if (erp->data[i].mjd==mjd-erp->data[i+1].mjd) {
a=0.5;
}
else {
a=(mjd-erp->data[i+1].mjd)/(erp->data[i].mjd-mjd-erp->data[i+1].mjd);
if (i+1>=erp->n||i<0) {
printf("i+1>=erp->n || i<0 %d\n",i);
getchar();
}
}
erpv[0]=(1.0-a)*erp->data[i].xp +a*erp->data[i+1].xp;
erpv[1]=(1.0-a)*erp->data[i].yp +a*erp->data[i+1].yp;
erpv[2]=(1.0-a)*erp->data[i].ut1_utc+a*erp->data[i+1].ut1_utc;
erpv[3]=(1.0-a)*erp->data[i].lod +a*erp->data[i+1].lod;
return 1;
}
4、tide_solid()：计算固体潮
固体潮为在太阳、月球等星体的引潮力的作用下，固体地球产生的同期性形变现象；虽然太阳的质量比月球大，但是由于太阳距离地球比月球距离地球远，因此太阳引起的固体潮相对月球较小。其他距离地球更远的天体，其固体潮效应也可忽略不计。随着作用点的位置的变化和太阳、月球位置的变化，固体潮的大小和方向也在不断变化，垂直方向最高可达 30～40 cm。在距离较短的相对定位中，固体潮的影响可以通过差分的方式消除，但精密单点定位不能通过差分的方式消除固体潮的影响，因此需要利用模型对其进行改正。固体潮二阶潮汐项对测站坐标影响的计算公式为：
Δ
r
=
∑
j
=
2
G
M
G
M
r
4
R
j
3
{
[
3
l
2
(
R
^
,
⋅
r
^
)
]
R
^
,
+
[
3
(
h
2
2
−
l
2
)
(
R
^
,
⋅
r
^
)
2
−
h
2
2
]
r
^
}
\Delta \boldsymbol{r}=\sum_{j=2} \frac{G M}{G M} \frac{\boldsymbol{r}^{4}}{\boldsymbol{R}_{j}^{3}}\left\{\left[3 l_{2}(\hat{\boldsymbol{R}}, \cdot \hat{\boldsymbol{r}})\right] \hat{\boldsymbol{R}}_{,}+\left[3\left(\frac{h_{2}}{2}-l_{2}\right)(\hat{\boldsymbol{R}}, \cdot \hat{\boldsymbol{r}})^{2}-\frac{h_{2}}{2}\right] \hat{\boldsymbol{r}}\right\}
Δr=j=2∑GMGMRj3r4{[3l2(R^,⋅r^)]R^,+[3(2h2−l2)(R^,⋅r^)2−2h2]r^}
式中：
-
G M G M GM：地球的引力常数
-
G M 2 G M_{2} GM2：月球的引力常数
-
G M 3 G M_{3} GM3：太阳的引力常数
-
r r r：地球在地心坐标系中的视线向量的模
-
r ^ \hat{r} r^ ： r r r 的単位向量
-
R 2 \boldsymbol{R}_{2} R2：月球在地心坐标系中的位置向量的模， R ^ 2 \hat{R}_{2} R^2 为 R 2 R_{2} R2 的单位向量
-
R 3 \boldsymbol{R}_{3} R3：月球在地心坐标系中的位置向量的模， R ^ 3 \hat{R}_{3} R^3 为 R 3 R_{3} R3 的单位向量
-
l 2 ， h 2 l_{2}，h_{2} l2，h2：二阶 Love 数和 Shida 数：
h 2 = 0.6087 − 0.0006 P 2 ( sin φ ) l 2 = 0.0847 − 0.0002 P 2 ( sin φ ) } \left.\begin{array}{l} h_{2}=0.6087-0.0006 P_{2}(\sin \varphi) \\ l_{2}=0.0847-0.0002 P_{2}(\sin \varphi) \end{array}\right\} h2=0.6087−0.0006P2(sinφ)l2=0.0847−0.0002P2(sinφ)}
式中， φ \varphi φ 为测站的地心纬度； P 2 ( sin φ ) = 3 sin 2 φ − 1 2 P_{2}(\sin \varphi)=\frac{3 \sin ^{2} \varphi-1}{2} P2(sinφ)=23sin2φ−1 。固体潮二阶项为固体潮改正的主要部分,最大可达分米级。
固体潮三阶潮汐项对测站位置影响的公式为：
Δ
r
‾
=
∑
j
=
2
3
G
M
G
M
r
5
R
j
4
{
[
h
3
r
^
(
5
2
(
R
^
j
⋅
r
^
)
3
−
3
2
(
R
^
,
r
^
)
]
R
^
j
+
l
3
(
15
2
(
R
^
j
⋅
r
^
2
−
3
2
)
[
R
^
j
−
(
R
^
j
⋅
r
^
)
r
^
]
}
\begin{aligned} \Delta \overline{\boldsymbol{r}}=\sum_{j=2}^{3} \frac{\mathrm{GM}}{\mathrm{GM}} \frac{\boldsymbol{r}^{5}}{\boldsymbol{R}_{j}^{4}}\left\{\left[h_{3} \hat{\boldsymbol{r}}\left(\frac{5}{2}\left(\hat{\boldsymbol{R}}_{j} \cdot \hat{\boldsymbol{r}}\right)^{3}-\frac{3}{2}(\hat{\boldsymbol{R}}, \hat{\boldsymbol{r}})\right] \hat{\boldsymbol{R}}_{j}+\right.\right. \\ l_{3}\left(\frac{15}{2}\left(\hat{\boldsymbol{R}}_{j} \cdot \hat{\boldsymbol{r}}^{2}-\frac{3}{2}\right)\left[\hat{\boldsymbol{R}}_{j}-\left(\hat{\boldsymbol{R}}_{j} \cdot \hat{\boldsymbol{r}}\right) \hat{\boldsymbol{r}}\right]\right\} \end{aligned}
Δr=j=2∑3GMGMRj4r5{[h3r^(25(R^j⋅r^)3−23(R^,r^)]R^j+l3(215(R^j⋅r^2−23)[R^j−(R^j⋅r^)r^]}
式中，
h
3
=
0.292
;
l
3
=
0.015
h_{3}=0.292 ; l_{3}=0.015
h3=0.292;l3=0.015。
一般由上式计算出的由太阳引起的固体湖改正较小，可以忽略不计，月球引起的固体潮改正为毫米级。
此外，固体潮的影响除了使测站产生周期性的位移外，其引湖位的零频项还会引起测站的永久性位移。设测站在径向和北向的永久性位移分别为
U
r
U_{\mathrm{r}}
Ur 和
U
n
U_{\mathrm{n}}
Un，则有：
U
r
=
[
−
0.1206
+
0.001
P
2
(
sin
φ
)
]
P
2
(
sin
φ
)
U
n
=
[
−
0.0252
−
0.001
P
2
(
sin
φ
)
]
sin
(
2
φ
)
}
\left.\begin{array}{l} U_{\mathrm{r}}=\left[-0.1206+0.001 P_{2}(\sin \varphi)\right] P_{2}(\sin \varphi) \\ U_{\mathrm{n}}=\left[-0.0252-0.001 P_{2}(\sin \varphi)\right] \sin (2 \varphi) \end{array}\right\}
Ur=[−0.1206+0.001P2(sinφ)]P2(sinφ)Un=[−0.0252−0.001P2(sinφ)]sin(2φ)}
代码中改正项比上面介绍的好像多一点，先调用两次 tide_pl() 计算日月潮汐二阶三阶改正项，然后计算频域改正项，再计算永久形变项。
static void tide_pl(const double *eu, const double *rp, double GMp,
const double *pos, double *dr)
{
const double H3=0.292,L3=0.015; // 三阶 Love、Shida 数
double r,ep[3],latp,lonp,p,K2,K3,a,H2,L2,dp,du,cosp,sinl,cosl;
int i;
if ((r=norm(rp,3))<=0.0) return;
for (i=0;i<3;i++) ep[i]=rp[i]/r;
K2=GMp/GME*SQR(RE_WGS84)*SQR(RE_WGS84)/(r*r*r);
K3=K2*RE_WGS84/r;
latp=asin(ep[2]); lonp=atan2(ep[1],ep[0]);
cosp=cos(latp); sinl=sin(pos[0]); cosl=cos(pos[0]);
// 二阶项
/* step1 in phase (degree 2) */
p=(3.0*sinl*sinl-1.0)/2.0;
H2=0.6078-0.0006*p; // 二阶 Love 数
L2=0.0847+0.0002*p; // 二阶 Shida 数
a=dot(ep,eu,3);
dp=K2*3.0*L2*a;
du=K2*(H2*(1.5*a*a-0.5)-3.0*L2*a*a);
// 三阶项
/* step1 in phase (degree 3) */
dp+=K3*L3*(7.5*a*a-1.5);
du+=K3*(H3*(2.5*a*a*a-1.5*a)-L3*(7.5*a*a-1.5)*a);
/* step1 out-of-phase (only radial) */
du+=3.0/4.0*0.0025*K2*sin(2.0*latp)*sin(2.0*pos[0])*sin(pos[1]-lonp);
du+=3.0/4.0*0.0022*K2*cosp*cosp*cosl*cosl*sin(2.0*(pos[1]-lonp));
dr[0]=dp*ep[0]+du*eu[0];
dr[1]=dp*ep[1]+du*eu[1];
dr[2]=dp*ep[2]+du*eu[2];
}
static void tide_solid(const double *rsun, const double *rmoon,
const double *pos, const double *E, double gmst, int opt,
double *dr)
{
double dr1[3],dr2[3],eu[3],du,dn,sinl,sin2l;
// 时域
/* step1: time domain */
eu[0]=E[2]; eu[1]=E[5]; eu[2]=E[8];
tide_pl(eu,rsun, GMS,pos,dr1);
tide_pl(eu,rmoon,GMM,pos,dr2);
// 频域，只有 K1 半径
/* step2: frequency domain, only K1 radial */
sin2l=sin(2.0*pos[0]);
du=-0.012*sin2l*sin(gmst+pos[1]);
dr[0]=dr1[0]+dr2[0]+du*E[2];
dr[1]=dr1[1]+dr2[1]+du*E[5];
dr[2]=dr1[2]+dr2[2]+du*E[8];
// 消除永久形变
/* eliminate permanent deformation */
if (opt&8) {
sinl=sin(pos[0]);
du=0.1196*(1.5*sinl*sinl-0.5);
dn=0.0247*sin2l;
dr[0]+=du*E[2]+dn*E[1];
dr[1]+=du*E[5]+dn*E[4];
dr[2]+=du*E[8]+dn*E[7];
}
}
5、tide_oload()：计算海潮负荷
海洋在日月引力等作用下产生潮汐变化，使得实际海平面相对于平均海平面发生周期性涨落,称为海洋潮。固体地球对海水质量潮汐分布产生的弹性响应称为海洋负荷效应。海洋潮的影响随测站与海洋距离的增加而不断减弱，近海地区可达厘米级，离海洋较远的地区约为毫米级。对于高精度精密单点定位而言，沿海地区应考虑海洋潮的影响，距离海洋 1000 k m 1000 \mathrm{~km} 1000 km 以上的测站，海洋潮影响可不予考虑。
由 IERS 2010 可知，测站的海洋潮瞬时位移
Δ
c
k
(
k
=
1
\Delta c_{k}(k=1
Δck(k=1 表示为东方向、
k
=
2
k=2
k=2 为北方向、
k
=
3
k=3
k=3 为垂直方向) 可表示为 11 个潮波 (4 个半日潮波
M
2
,
S
2
,
N
2
,
K
2
,
4
M_{2}, S_{2}, N_{2}, K_{2}, 4
M2,S2,N2,K2,4 个周日潮波
K
1
K_{1}
K1,
O
1
,
P
1
,
Q
1
O_{1}, P_{1}, Q_{1}
O1,P1,Q1 和 3 个长周期潮波
M
f
,
M
m
,
S
s
a
)
\left.M_{f}, M_{m}, S_{s a}\right)
Mf,Mm,Ssa) 海洋潮共同影响的矢量叠加：
Δ
c
k
=
∑
j
=
1
11
f
j
A
k
,
j
cos
(
ω
j
t
+
χ
j
+
μ
j
−
Φ
k
,
j
)
\Delta c_{k}=\sum_{j=1}^{11} f_{j} A_{k, j} \cos \left(\omega_{j} t+\chi_{j}+\mu_{j}-\Phi_{k, j}\right)
Δck=j=1∑11fjAk,jcos(ωjt+χj+μj−Φk,j)
式中，
A
k
,
j
A_{k, j}
Ak,j 和
Φ
k
,
j
\Phi_{k, j}
Φk,j 分别表示潮波
j
j
j 在
k
k
k 方向的振幅和 Greenwich 相位；
ω
j
\omega_{j}
ωj 和
χ
j
\chi_{j}
χj 分别为潮波
j
j
j 的角频率和天文幅角；
f
j
f_{j}
fj 和
μ
j
\mu_{j}
μj 为顾及月亮轨道升交点调节作用的参数 (周期约为
18.6
a
18.6 \mathrm{a}
18.6a )，分别称为节点因数和天文相角。在
1
∼
3
m
m
1 \sim 3 \mathrm{~mm}
1∼3 mm 的海洋负荷位移改正精度下，可以认为
f
j
=
1
f_{j}=1
fj=1 和
μ
j
=
0
\mu_{j}=0
μj=0，则式可以简化为：
Δ
c
k
=
∑
j
=
1
11
A
k
,
j
cos
(
ω
j
t
+
χ
j
−
Φ
k
,
j
)
\Delta c_{k}=\sum_{j=1}^{11} A_{k, j} \cos \left(\omega_{j} t+\chi_{j}-\Phi_{k, j}\right)
Δck=j=1∑11Ak,jcos(ωjt+χj−Φk,j)
不同研究机构根据卫星测高、船测等数据建立了多个海洋潮汐模型，如 CSR、OSU、FES、GOT 等都有多个版本，不同模型之间的差异可以达到 3mm 左右。
GAMP 中先定义了 11 个潮波参数 args，然后计算当前时间的天文参数，再计算由 11 个潮波引起的位移。
static void tide_oload(gtime_t tut, const double *odisp, double *denu)
{
// 11 个潮波定义
const double args[][5]={
{1.40519E-4, 2.0,-2.0, 0.0, 0.00}, /* M2 */
{1.45444E-4, 0.0, 0.0, 0.0, 0.00}, /* S2 */
{1.37880E-4, 2.0,-3.0, 1.0, 0.00}, /* N2 */
{1.45842E-4, 2.0, 0.0, 0.0, 0.00}, /* K2 */
{0.72921E-4, 1.0, 0.0, 0.0, 0.25}, /* K1 */
{0.67598E-4, 1.0,-2.0, 0.0,-0.25}, /* O1 */
{0.72523E-4,-1.0, 0.0, 0.0,-0.25}, /* P1 */
{0.64959E-4, 1.0,-3.0, 1.0,-0.25}, /* Q1 */
{0.53234E-5, 0.0, 2.0, 0.0, 0.00}, /* Mf */
{0.26392E-5, 0.0, 1.0,-1.0, 0.00}, /* Mm */
{0.03982E-5, 2.0, 0.0, 0.0, 0.00} /* Ssa */
};
const double ep1975[]={1975,1,1,0,0,0};
double ep[6],fday,days,t,t2,t3,a[5],ang,dp[3]={0};
int i,j;
// 角度参数
/* angular argument: see subroutine arg.f for reference [1] */
time2epoch(tut,ep);
fday=ep[3]*3600.0+ep[4]*60.0+ep[5];
ep[3]=ep[4]=ep[5]=0.0;
days=timediff(epoch2time(ep),epoch2time(ep1975))/86400.0+1.0;
t=(27392.500528+1.000000035*days)/36525.0;
t2=t*t; t3=t2*t;
a[0]=fday;
a[1]=(279.69668+36000.768930485*t+3.03E-4*t2)*D2R; /* H0 */
a[2]=(270.434358+481267.88314137*t-0.001133*t2+1.9E-6*t3)*D2R; /* S0 */
a[3]=(334.329653+4069.0340329577*t-0.010325*t2-1.2E-5*t3)*D2R; /* P0 */
a[4]=2.0*PI;
// 计算由 11 个潮波引起的位移
/* displacements by 11 constituents */
for (i=0;i<11;i++) {
ang=0.0;
for (j=0;j<5;j++) ang+=a[j]*args[i][j];
for (j=0;j<3;j++) dp[j]+=odisp[j+i*6]*cos(ang-odisp[j+3+i*6]*D2R);
}
denu[0]=-dp[1];
denu[1]=-dp[2];
denu[2]= dp[0];
}
6、tide_pole()：计算极潮
极移 (Polar Motion) 是地球瞬时自转轴在地球本体内的运动，表现为以极点为中心，在
±
4
′
′
\pm 4^{\prime \prime}
±4′′ (相当于
24
m
×
24
m
24 \mathrm{~m} \times 24 \mathrm{~m}
24 m×24 m ) 范围内，按与地球自转相同的方向形成一条时伸时缩的螺旋形曲线，即由于地壳对极移的弹性响应造成的地球自转轴变化。
假设测站的坐标为
(
φ
,
λ
,
h
)
(\varphi, \lambda, h)
(φ,λ,h)，极潮在径向与平面对该测站的影响分别为：
S
r
=
−
32
sin
(
2
θ
)
(
m
1
cos
λ
+
m
2
sin
λ
)
S
θ
=
−
9
cos
(
2
θ
)
(
m
1
cos
λ
+
m
2
sin
λ
)
S
λ
=
9
cos
θ
(
m
1
sin
λ
−
m
2
cos
λ
)
}
\left.\begin{array}{l} S_{r}=-32 \sin (2 \theta)\left(m_{1} \cos \lambda+m_{2} \sin \lambda\right) \\ S_{\theta}=-9 \cos (2 \theta)\left(m_{1} \cos \lambda+m_{2} \sin \lambda\right) \\ S_{\lambda}=9 \cos \theta\left(m_{1} \sin \lambda-m_{2} \cos \lambda\right) \end{array}\right\}
Sr=−32sin(2θ)(m1cosλ+m2sinλ)Sθ=−9cos(2θ)(m1cosλ+m2sinλ)Sλ=9cosθ(m1sinλ−m2cosλ)⎭
⎬
⎫
式中，
S
S
S，为径向变形 (向上为正)，
S
θ
S_{\theta}
Sθ 为平面变形 (南北方向，向南为正) 、
S
λ
S_{\lambda}
Sλ 为平面变形 (东西方向，向东为正)，单位均为毫米。其中
θ
=
π
/
2
−
φ
\theta=\pi / 2-\varphi
θ=π/2−φ,称为余纬。假设计算极潮改正时刻对应的 ERP 极移参数为
(
x
p
,
y
p
)
\left(x_{p}, y_{p}\right)
(xp,yp)，则有如下辅助公式：
m
1
=
x
p
−
x
ˉ
p
,
m
2
=
−
(
y
p
−
y
ˉ
p
)
x
ˉ
p
(
t
)
=
x
ˉ
p
(
t
0
)
+
(
t
−
t
0
)
x
ˉ
˙
p
(
t
0
)
y
ˉ
p
(
t
)
=
y
ˉ
p
(
t
0
)
+
(
t
−
t
0
)
y
ˉ
˙
p
(
t
0
)
}
\left.\begin{array}{c} m_{1}=x_{p}-\bar{x}_{p}, m_{2}=-\left(y_{p}-\bar{y}_{p}\right) \\ \bar{x}_{p}(t)=\bar{x}_{p}\left(t_{0}\right)+\left(t-t_{0}\right) \dot{\bar{x}}_{p}\left(t_{0}\right) \\ \bar{y}_{p}(t)=\bar{y}_{p}\left(t_{0}\right)+\left(t-t_{0}\right) \dot{\bar{y}}_{p}\left(t_{0}\right) \end{array}\right\}
m1=xp−xˉp,m2=−(yp−yˉp)xˉp(t)=xˉp(t0)+(t−t0)xˉ˙p(t0)yˉp(t)=yˉp(t0)+(t−t0)yˉ˙p(t0)⎭
⎬
⎫
式中：
x
ˉ
p
(
t
0
)
=
0.054
,
x
ˉ
˙
p
(
t
0
)
=
0.00083
y
ˉ
p
(
t
0
)
=
0.357
,
y
ˉ
˙
p
(
t
0
)
=
0.00395
}
\left.\begin{array}{l} \bar{x}_{p}\left(t_{0}\right)=0.054, \dot{\bar{x}}_{p}\left(t_{0}\right)=0.00083 \\ \bar{y}_{p}\left(t_{0}\right)=0.357, \dot{\bar{y}}_{p}\left(t_{0}\right)=0.00395 \end{array}\right\}
xˉp(t0)=0.054,xˉ˙p(t0)=0.00083yˉp(t0)=0.357,yˉ˙p(t0)=0.00395}
式中，
x
ˉ
p
,
y
ˉ
p
\bar{x}_{p}, \bar{y}_{p}
xˉp,yˉp 单位为弧秒，对应的速度
x
ˉ
p
,
y
ˉ
p
\bar{x}_{p}, \bar{y}_{p}
xˉp,yˉp 单位为弧秒/年，
m
1
,
m
2
m_{1}, m_{2}
m1,m2 单位也为弧秒，
t
0
=
t_{0}=
t0= 2000.0 ，对应的时间为 2000 年 1 月 1 日 12 时 0 分 0 秒。
考虑到 m 1 、 m 2 m_{1} 、 m_{2} m1、m2 最大可能至 0.8 弧秒，因此极潮改正在径向最大可到 25 m m 25 \mathrm{~mm} 25 mm，而在平面上变形最大约为 7 m m 7 \mathrm{~mm} 7 mm 。因此，在高精度精密单点定位中应考虑此项极潮的影响。
实际应用中，如需将以站心地平坐标系中表示的极潮位移改正转化到空间直角坐标系中，可以采用如下公式进行：
[
d
X
d
Y
d
Z
]
T
=
R
T
[
S
θ
S
λ
S
r
]
T
\left[\begin{array}{lll} \mathrm{d} X & \mathrm{~d} Y & \mathrm{~d} Z \end{array}\right]^{\mathrm{T}}=\boldsymbol{R}^{\mathrm{T}}\left[\begin{array}{lll} S_{\theta} & S_{\lambda} & S_{r} \end{array}\right]^{\mathrm{T}}
[dX dY dZ]T=RT[SθSλSr]T
式中：
R
=
[
cos
θ
cos
λ
cos
θ
sin
λ
−
sin
θ
−
sin
λ
cos
λ
0
sin
θ
cos
λ
sin
θ
sin
λ
cos
θ
]
\boldsymbol{R}=\left[\begin{array}{ccc} \cos \theta \cos \lambda & \cos \theta \sin \lambda & -\sin \theta \\ -\sin \lambda & \cos \lambda & 0 \\ \sin \theta \cos \lambda & \sin \theta \sin \lambda & \cos \theta \end{array}\right]
R=
cosθcosλ−sinλsinθcosλcosθsinλcosλsinθsinλ−sinθ0cosθ
GAMP 中先调用 iers_mean_pole() 计算平均极点坐标 xp_bar、yp_bar，后面用的时候会再乘 1E-3，平均极点坐标计算与上面公式略有不同：
- 使用三次多项式来拟合 2000 年到 2010 年的平均极坐标
- 线性函数模型，拟合 2010 年以后的平均极坐标
然后用 ERP 参数和平均极点坐标计算 m1、m2，再套潮汐改正公式。
static void iers_mean_pole(gtime_t tut, double *xp_bar, double *yp_bar)
{
const double ep2000[]={2000,1,1,0,0,0};
double y,y2,y3;
y=timediff(tut,epoch2time(ep2000))/86400.0/365.25;
// 使用三次多项式来拟合 2000 年到 2010 年的平均极坐标
if (y<3653.0/365.25) { /* until 2010.0 */
y2=y*y; y3=y2*y;
*xp_bar= 55.974+1.8243*y+0.18413*y2+0.007024*y3; /* (mas) */
*yp_bar=346.346+1.7896*y-0.10729*y2-0.000908*y3;
}
// 线性函数模型，拟合 2010 年以后的平均极坐标
else { /* after 2010.0 */
*xp_bar= 23.513+7.6141*y; /* (mas) */
*yp_bar=358.891-0.6287*y;
}
}
static void tide_pole(gtime_t tut, const double *pos, const double *erpv,
double *denu)
{
double xp_bar,yp_bar,m1,m2,cosl,sinl;
// 计算平均极点坐标 xp_bar、yp_bar，后面用的时候会再乘 1E-3
/* iers mean pole (mas) */
iers_mean_pole(tut,&xp_bar,&yp_bar);
// 用 ERP 参数和平均极点坐标计算 m1、m2
/* ref [7] eq.7.24 */
m1= erpv[0]/AS2R-xp_bar*1E-3; /* (as) */
m2=-erpv[1]/AS2R+yp_bar*1E-3;
/* sin(2*theta) = sin(2*phi), cos(2*theta)=-cos(2*phi) */
cosl=cos(pos[1]);
sinl=sin(pos[1]);
denu[0]= 9E-3*sin(pos[0]) *(m1*sinl-m2*cosl); /* de= Slambda (m) */
denu[1]= -9E-3*cos(2.0*pos[0])*(m1*cosl+m2*sinl); /* dn=-Stheta (m) */
denu[2]=-33E-3*sin(2.0*pos[0])*(m1*cosl+m2*sinl); /* du= Sr (m) */
}
五、地球自转效应改正
GNSS 的 MEO 卫星轨道高度大约 20000km。导航信号由卫星发出到接收机接受需要数十微秒，对于 GEO 卫星信号传播时间更长。在导航信号传播过程中，由于地球自转的影响，地固坐标系已随地球旋转了一定角度，由此给观测值造成的误差可达数十米，其影响不可忽略。因此需要对观测值进行距离改正，改正数的计算方法如下：
Δ
D
=
ω
c
[
y
S
(
x
R
−
x
S
)
−
x
S
(
y
R
−
y
S
)
]
\Delta D=\frac{\omega}{c}\left[y^{S}\left(x_{R}-x^{S}\right)-x^{S}\left(y_{R}-y^{S}\right)\right]
ΔD=cω[yS(xR−xS)−xS(yR−yS)]
式中，
x
S
,
y
S
,
x
R
,
y
R
x^{S}, y^{S}, x_{R}, y_{R}
xS,yS,xR,yR 分别表示信号发射时刻卫星和接收机在地固坐标系下的坐标，
ω
\omega
ω 表示地球自转角速度，单位为弧度每秒。
地球自转效应改正在 geodist() 函数中计算站心几何距离的时候就已经改正了。GAMP 又写了 sagnac() 函数专门计算地球自转效应改正量，并不用于定位解算，只是赋值到 PPP_Info 里，用于生成 RCVEX 文件，代码与公式完全对应：
extern double sagnac(const double *rs, const double *rr)
{
return OMGE*(rs[0]*rr[1]-rs[1]*rr[0])/CLIGHT;
}
六、引力延迟改正
1、原理
广义相对论效应还会产生引力延迟。这里不加推导，直接给出引力延迟的计算公式如下：
Δ
D
g
=
2
μ
c
2
ln
r
+
R
+
ρ
r
+
R
−
ρ
\Delta \mathrm{D}_{\mathrm{g}}=\frac{2 \mu}{\mathrm{c}^{2}} \ln \frac{\mathrm{r}+\mathrm{R}+\rho}{\mathrm{r}+\mathrm{R}-\rho}
ΔDg=c22μlnr+R−ρr+R+ρ
式中：
μ
\mu
μ 为万有引力常数
G
\mathrm{G}
G 与地球总质量
M
\mathrm{M}
M 之乘积；
c
\mathrm{c}
c 为真空中的光速；
r
\mathrm{r}
r 为卫星至地心的距离；
R
R
R 为测站至地心的距离；
ρ
\rho
ρ 为测站至卫星的距离。当卫星接近地平面时引力延迟取得最大值，约为
19
m
m
19 \mathrm{~mm}
19 mm 。当卫星在测站天顶方向时引力延迟取得最小值，约为
13
m
m
13 \mathrm{~mm}
13 mm 。引力延迟引起的相对测距误差不足
1
0
−
9
10^{-9}
10−9，一般的单点定位中无需顾及，但在精密单点定位 (PPP) 中应予以考虑。在相对定位中，当站间距离为
1000
k
m
1000 \mathrm{~km}
1000 km 时，两站的引力延迟之差最大可达
1.3
m
m
1.3 \mathrm{~mm}
1.3 mm；当站间距离为
3000
k
m
3000 \mathrm{~km}
3000 km 时，两站的引力延迟之差最大可达
3.6
m
m
3.6 \mathrm{~mm}
3.6 mm 。只有在高精度相对定位中才需顾及此项改正。
GAMP 中不同系统的引力常数定义如下，在小数点后 13 位有区别：
#define MU_GPS 3.9860050E14 /* gravitational constant ref [1] */
#define MU_GLO 3.9860044E14 /* gravitational constant ref [2] */
#define MU_GAL 3.986004418E14 /* earth gravitational constant ref [7] */
#define MU_CMP 3.986004418E14 /* earth gravitational constant ref [9] */
2、gravitationalDelayCorrection()：引力延迟改正
gravitationalDelayCorrection() 代码中，先计算接收机到地心的距离 receiverModule、卫星到地心的距离 satelliteModule、卫星到接收机距离 distance，然后套公式计算。
因为函数传入的是 ECEF 坐标，地心坐标就是原点(0,0,0)，所以到地心距离直接就是 ECEF 坐标的模。
extern double gravitationalDelayCorrection (const int sys, const double *receiverPosition,
const double *satellitePosition)
{
double receiverModule; // 接收机到地心的距离
double satelliteModule; // 卫星到地心的距离
double distance; // 接收机到卫星的矩阵
double MU=MU_GPS;
// 先计算接收机到地心的距离 receiverModule、卫星到地心的距离 satelliteModule、卫星到接收机距离 distance
// 传入的是 ECEF 坐标，地心坐标是原点(0,0,0)，所以到地心距离直接就是 ECEF 坐标的模
receiverModule=sqrt(receiverPosition[0]*receiverPosition[0]+receiverPosition[1]*receiverPosition[1]+
receiverPosition[2]*receiverPosition[2]);
satelliteModule=sqrt(satellitePosition[0]*satellitePosition[0]+satellitePosition[1]*satellitePosition[1]+
satellitePosition[2]*satellitePosition[2]);
distance=sqrt((satellitePosition[0]-receiverPosition[0])*(satellitePosition[0]-receiverPosition[0])+
(satellitePosition[1]-receiverPosition[1])*(satellitePosition[1]-receiverPosition[1])+
(satellitePosition[2]-receiverPosition[2])*(satellitePosition[2]-receiverPosition[2]));
// 不同系统引力常数在小数点后 13 位有区别
switch (sys) {
case SYS_GPS:
MU=MU_GPS;
break;
case SYS_GLO:
MU=MU_GLO;
break;
case SYS_GAL:
MU=MU_GAL;
break;
case SYS_CMP:
MU=MU_CMP;
break;
default:
MU=MU_GPS;
break;
}
// 套公式计算
return 2.0*MU/(CLIGHT*CLIGHT)*log((satelliteModule+receiverModule+distance)/(satelliteModule+receiverModule-distance));
}
MU_GAL 3.986004418E14 /* earth gravitational constant ref [7] /
#define MU_CMP 3.986004418E14 / earth gravitational constant ref [9] */
### 2、gravitationalDelayCorrection()：引力延迟改正
`gravitationalDelayCorrection()` 代码中，先计算接收机到地心的距离 `receiverModule`、卫星到地心的距离 `satelliteModule`、卫星到接收机距离 `distance`，然后套公式计算。
> 因为函数传入的是 ECEF 坐标，地心坐标就是原点(0,0,0)，所以到地心距离直接就是 ECEF 坐标的模。
```c
extern double gravitationalDelayCorrection (const int sys, const double *receiverPosition,
const double *satellitePosition)
{
double receiverModule; // 接收机到地心的距离
double satelliteModule; // 卫星到地心的距离
double distance; // 接收机到卫星的矩阵
double MU=MU_GPS;
// 先计算接收机到地心的距离 receiverModule、卫星到地心的距离 satelliteModule、卫星到接收机距离 distance
// 传入的是 ECEF 坐标，地心坐标是原点(0,0,0)，所以到地心距离直接就是 ECEF 坐标的模
receiverModule=sqrt(receiverPosition[0]*receiverPosition[0]+receiverPosition[1]*receiverPosition[1]+
receiverPosition[2]*receiverPosition[2]);
satelliteModule=sqrt(satellitePosition[0]*satellitePosition[0]+satellitePosition[1]*satellitePosition[1]+
satellitePosition[2]*satellitePosition[2]);
distance=sqrt((satellitePosition[0]-receiverPosition[0])*(satellitePosition[0]-receiverPosition[0])+
(satellitePosition[1]-receiverPosition[1])*(satellitePosition[1]-receiverPosition[1])+
(satellitePosition[2]-receiverPosition[2])*(satellitePosition[2]-receiverPosition[2]));
// 不同系统引力常数在小数点后 13 位有区别
switch (sys) {
case SYS_GPS:
MU=MU_GPS;
break;
case SYS_GLO:
MU=MU_GLO;
break;
case SYS_GAL:
MU=MU_GAL;
break;
case SYS_CMP:
MU=MU_CMP;
break;
default:
MU=MU_GPS;
break;
}
// 套公式计算
return 2.0*MU/(CLIGHT*CLIGHT)*log((satelliteModule+receiverModule+distance)/(satelliteModule+receiverModule-distance));
}
