GPS卫星位置的计算

1. GPS卫星位置的计算

1.1 用广播星历计算卫星位置

1.1.1. 计算卫星运动的平均角速度 $n$

首先根据广播星历中给出的参数 $\sqrt{A}$ 计算参考时刻 $t_{\mathrm{oe}}$ 的平均角速度 $n_0$ :
$$n_0=\frac{\sqrt{GM}}{(\sqrt{A}{)}^3}$$
式中, $GM$ 为万有引力常数 $G$ 与地球总质量 $M$ 之乘积, 其值为 $GM=3.986005\times 10^{14}{\mathrm{m}}^3/{\mathrm{s}}^2$ 。 然后根据广播星历中给定的摄动参数 $\Delta n$ 计算观测时刻卫星的平均角速度 $n$ :
$$n=n_0+\Delta n$$

1.1.2. 计算观测睶间卫星的平近点角 $M$

$$M=M_0+n\left(t-t_{oe}\right)$$
式中, $M_0$ 为参考时刻 $t_{oe}$ 时的平近点角, 由广播星历给出。
为什么要用参考时刻 $t_{oe}$ 来替代卫星过近地点时刻 $t_0$ 来计算呢? 原因很简单, 因为广播 星历每 $2\mathrm{h}$ 更新一次, 将参考时刻设在中央时刻时, 外推间隔小于等于 $1\mathrm{h}$ 。而卫星的运行周 期为 $12\mathrm{h}$ 左右, 采用卫星过近地点时刻 $t_0$ 来计算时, 外推间隔最大有可能达 $6\mathrm{h}$ 。用 $t_{oe}$ 来取 代卫星过近地点时刻 $t_0$ 后, 外推间隔将大大减小, 用较简单的模型也能获得精度较高的结 果。

1.1.3. 计算偏近点角

用弧度表示的开普勒方程为:
$$E=M+e\sin E$$
用角度表示的开普勒方程为:
$$E^{\circ }=M^{\circ }+{\rho }^{\circ }\cdot e\sin E^{\circ }$$
解上述方程可用迭代法或微分改正法。

1.1.4. 计算真近点角 $f$

$$\{\begin{array}{l}\cos f=\frac{\cos E-e}{1-e\cos E}\\ \sin f=\frac{\sqrt{1-e^2}\sin E}{1-e\cos E}\end{array}$$
式中, $e$ 为卫星轨道的偏心率, 由广播星历给出。
所以,
$$f=\arctan \frac{\sqrt{1-e^2}\sin E}{\cos E-e}$$

1.1.5. 计算升交角距 $u^{\prime }$

$$u^{\prime }=\omega +f$$
式中, $\omega$ 为近地点角距, 由广播星历给出。

1.1.6. 计算摄动改正项 ${\delta }_{\mu }、{\delta }_r、{\delta }_i$

广播星历中给出了下列 6 个摄动参数: $C_{uc},C_{\mathrm{we}},C_{rc},C_r,C_{ic},C_{is}$, 据此可求出由于 $J_2$ 项 而引起的升交角距 $u$ 的改正项 ${\delta }_u$ 、卫星矢径 $r$ 的改正项 ${\delta }_r$ 和卫星轨道倾角 $i$ 的摄动改正项 ${\delta }_i$ 。 计算公式如下:
$$\{\begin{array}{l}{\delta }_u=C_{uk}\cos 2u^{\prime }+C_w\sin 2u^{\prime }\\ {\delta }_r=C_{rc}\cos 2u^{\prime }+C_{rr}\sin 2u^{\prime }\\ {\delta }_i=C_{ic}\cos 2u^{\prime }+C_{iu}\sin 2u^{\prime }\end{array}$$

1.1.7. 对 $u^{\prime }、r^{\prime }、i_0$ 进行掇动改正

$$\{\begin{array}{l}u=u^{\prime }+{\delta }_u\\ r=r^{\prime }+{\delta }_r=a(1-e\cos E)+{\delta }_r\\ i=i_0+{\delta }_i+\frac{\mathrm{d}i}{\mathrm{d}t}\left(t-t_{\partial t}\right)\end{array}$$
式中, $a$ 为卫星轨道的长半径, $a=(\sqrt{A}{)}^2,\sqrt{A}$ 由广播星历给出; $i_0$ 为 $t_{\mathrm{oe}}$ 时刻的轨道倾角, 由 广播星历中的开普勒六参数给出; $\frac{\mathrm{d}i}{\mathrm{d}t}$ 为 $i$ 的变化率, 由广播星历中的摄动九参数给出。

1.1.8. 计算卫星在轨道面坐标系中的位置

在轨道平面直角坐标系中 (坐标原点位于地心, $X$ 轴指向升交点) 卫星的平面直角坐 标为:
$$\{\begin{array}{l}x=r\cos u\\ y=r\sin u\end{array}$$

1.1.9. 计算观测瞬间升交点的经度 $L$

若参考时刻 $t_{\mathrm{oe}}$ 时升交点的赤经为 ${\Omega }_{t_{\alpha e}}$, 升交点对时间的变化率为 $\dot{\Omega }$, 那么观测瞬间 $t$ 的 升交点赤经 $\Omega$ 应为:
$$\mathbf{\Omega }={\Omega }_{t_{ee}}+\dot{\Omega }\left(t-t_{ot}\right)$$
$\dot{\Omega }$ 可从广播星历的摄动参数中给出。
设本周开始时刻 (星期日 0 时) 格林尼治恒星时为 ${\mathrm{GAST}}_{\text{veek }}$, 则观测瞬间的格林尼治恒 星时为:
$$\text{ GAST }={\mathrm{GAST}}_{\text{week }}+{\omega }_{e}t$$
式中, ${\omega }_e$ 为地球自转角速度, 其值为 ${\omega }_e=7.292115\times 10^{-5}\mathrm{rad}/{\mathrm{s}}_{\text{。 }}$

$t$ 为本周内的时间 $(\mathrm{s})$ 。显然, 上述 @埓中把地球自转看成是完全匀速的, 末顾及地球 自转的不均匀性。这样就可求得观测䀰间昲交点的经度值为:
$$\begin{array}{c}L=\Omega -\mathrm{GAST}={\Omega }_{t_{oe}}-{\mathrm{GAST}}_{\text{week }}+\dot{\Omega }\left(t-t_{oe}\right)-{\omega }_{e}t\\ {\Omega }_0={\Omega }_{t_{oe}}-{\mathrm{GAST}}_{\text{week }}\end{array}$$
则有:
$$L={\Omega }_0+\dot{\Omega }\left(t-t_{oe}\right)-{\omega }_{e}t={\Omega }_0+\left(\dot{\Omega }-{\omega }_e\right)t-\dot{\Omega }\cdot t_{oe}$$
注意: 广播星历中给出的并不是参考时刻 $t_{oe}$ 的升交点赤经 ${\Omega }_{t_{ec}}$, 而是该值与本周起始时 刻的格林尼治恒星时 ${\mathrm{GAST}}_{\mathrm{weck}}$ 之差。

1.1.10. 计算卫星在睬时地球坐标系中的位置

已知升交点的大地经度 $L$ 以及轨道平面的倾角 $i$ 后, 就可通过两次旋转方便地求得卫 星在地固坐标系中的位置:
$$\left(\begin{array}{l}X\\ Y\\ Z\end{array}\right)={\mathbf{R}}_Z(-L){\mathbf{R}}_X(-i)\left(\begin{array}{l}x\\ y\\ 0\end{array}\right)=\left(\begin{array}{c}x\cos L-y\cos i\sin L\\ x\sin L+y\cos i\cos L\\ y\sin i\end{array}\right)$$

1.1.11. 计算卫星在协议地球坐标系中的位畳

观测瞬间卫星在协议地球坐标系中的位置为:
$$\begin{array}{rl}\left(\begin{array}{l}x\\ y\\ z\end{array}\right)& ={\mathbf{R}}_Y\left(-x_p\right){\mathbf{R}}_x\left(-y_p\right)\cdot \left(\begin{array}{l}X\\ Y\\ Z\end{array}\right)\\ & =\left(\begin{array}{ccc}1& 0& x_p\\ 0& 1& -y_p\\ -x_p& y_p& 1\end{array}\right)\cdot \left(\begin{array}{l}X\\ Y\\ Z\end{array}\right)\end{array}$$
在 GPS 定位中, 常常需要多次计算卫星的位置和速度, 如用上述方法计算需占用较多的内存和计算时间。为此, 常将卫星星历用一个时间多项式来表示, 在内存中仅保存该多项式的系数, 供计算时调用。在各种多项式中, 切比雪夫多项式的拟合效果最佳, 即使在该时 间段的两端近似性也很好。用 $n$ 阶切比雪夫多项式来逼近时间段 $\left[t_0,t_0+\Delta t\right]$ 中的卫星星历时, 先将变量 $t\in \left[t_0,t_0+\Delta t\right]$ 变换为变量 $\tau \in [-1,1]:$
$$\tau =\frac{2}{\Delta t}\left(t-t_0\right)-1$$
于是卫星坐标可表示为:
$$X(t)=\sum _{i=0}^n{C}_{x_i}{T}_i(t)$$
式中, $n$ 为多项式的阶数; $C_{x_i}$ 为切比雪夫多项式的系数。根据已知的卫星坐标, 用最小二乘 法拟合出多项式系数 $C_{x_i}$, 就可用式(3-27) 计算出该时段中任一时刻的卫星位置。切比雪 夫多项式 $T_i$ 的递推公式如下:
$$T_0(\tau )=1$$

$$\begin{array}{rl}& T_1(\tau )=\tau \\ & T_n(\tau )=2\tau T_{n-1}(\tau )-T_{n-2}(\tau ),|\tau |⩽1,n⩾2\end{array}$$

2 用精密星历计算卫星位置

精密星历是按一定的时间间隔 (通常为 $15\min$ ) 来给出卫星在空间的三维坐标、三维运 动速度及卫星钟改正数等信息。著名的 IGS 综合精密星历需 $1\sim 2$ 周后才能获得。由 NGS 提出的格式被广泛采用,其中 ASCII 格式的 SP1 和二进制格式的 ECF1 格式不仅给出了卫 星的三维位置信息 $(\mathrm{km})$, 也给出了卫星的三维运动速度信息 $(\mathrm{km}/\mathrm{s})$ 。而 SP2(ASCII 格式) 和 ECF2 (二进制格式) 则仅给出了卫星三维位置信息。速度信息需通过位置信息用数值微 分的方法来求出。采用这种格式时存储量可减少一半左右。在 SP3(ASCII 格式) 和 ECF3 (二进制格式) 中增加了卫星钟的改正数信息。
观测瞬间的卫星位置及运动速度可采用内揷法求得。其中拉格朗日 (Lagrange) 多项式 内揷法被广泛采用, 因为这种内揷法速度快且易于编程。拉格朗日揷值公式十分简单: 已知 函数 $y=f(x)$ 的 $n+1$ 个节点 $x_0,x_1,x_2,\cdots ,x_n$ 及其对应的函数值 $y_0,y_1,y_2,\cdots ,y_n$, 对揷值区 间内任一点 $x$, 可用下面的拉格朗日揷值多项式来计算函数值:
$$f(x)=\sum _{k=0}^n\prod _{i=0,i\ne k}^n\left(\frac{x-x_i}{x_k-x_i}\right)y_k$$
Remondi 的研究表明, 对 GPS 卫星而言, 如果要精确至 $10^{-8}$, 用 $30\min$ 的历元间隔和 9 阶内插已足够保证精度。