C语言计算GPS卫星位置

C

语言计算G P S 卫星位置

1 概述 在用GPS 信号进行导航定位以及制订观测计划时，都必须已知GPS 卫星在空间的瞬间位置。卫星位置的计算是根据卫星电文所提供的轨道参数按一定的公式计算的。本节专门讲解观测瞬间GPS 卫星在地固坐标系中坐标的计算方法。

2 卫星位置的计算

1. 计算卫星运行的平均角速度n

根据开普勒第三定律，卫星运行的平均角速度n0可以用下式计算：

式中μ为WGS-84坐标系中的地球引力常数，且μ=3.986005×1014m 3/s 2。平均角速度n 0加上卫星电文给出的摄动改正数Δn ，便得到卫星运行的平均角速度n

n=n 0+Δn (4-12)

2. 计算归化时间t k

首先对观测时刻t ′作卫星钟差改正

t=t ′-Δt

然后对观测时刻t 归化到GPS 时系

t k =t-t oc (4-13)

式中t k 称作相对于参考时刻t oe 的归化时间（读者注意：toc ≠toe ）。

3. 观测时刻卫星平近点角M k 的计算

M k =M 0+n tk (4-14)

式中M 0是卫星电文给出的参考时刻toe 的平近点角。

4. 计算偏近点角E k

E k =M k +esinE k (E k ,M k 以弧度计) （4-15）

上述方程可用迭代法进行解算，即先令E k =M k ,代入上式，求出E k 再代入上式计

算，因为GPS 卫星轨道的偏心率e 很小，因此收敛快，只需迭代计算两次便可求得偏近点角E k 。

5. 真近点角V k 的计算

由于:

因此:

6．升交距角Φk 的计算

ω为卫星电文给出的近地点角距。

7. 摄动改正项δu,δr,δi 的计算

δu,δr,δi 分别为升交距角u 的摄动量，卫星矢径r 的摄动量和轨道

倾角i 的摄动量。

8. 计算经过摄动改正的升交距角u k 、卫星矢径r k 和轨道倾角i k

9. 计算卫星在轨道平面坐标系的坐标

卫星在轨道平面直角坐标系（X 轴指向升交点）中的坐标为

10. 观测时刻升交点经度Ωk 的计算

升交点经度Ωk 等于观测时刻升交点赤经Ω（春分点和升交点之间的角距）与

格林泥治视恒星时GAST （春分点和格林尼治起始子午线之间的角距）之差，

Ωk =Ω-GAST （4-23）

又因为:

tk oe Ω+Ω=Ω (4-24)

其中Ωoe 为参与时刻t oe 的升交点的赤经；

Ω

是升交点赤经的变化率，卫星电文每小时更新一次Ω和t oe 。 此外，卫星电文中提供了一周的开始时刻t w 的格林尼治视恒星时GAST w 。由于

地球自转作用，GAST 不断增加，所以:

GAST=GAST w +ωe t （4-25）

式中ωe ×10-5rad/s 为地球自转的速率；t 为观测时刻。

由式（4-24）和（4-25），得:

由（4-13）式，得:

其中0oe w GAST Ω=Ω-，o Ω、Ω、oe t 的值可从卫星电文中获取。

11. 计算卫星在地心固定坐标系中的直角坐标

把卫星在轨道平面直角坐标系中的坐标进行旋转变换，可得出卫星在地心固定坐标系中的三维坐标：

12. 卫星在协议地球坐标系中的坐标计算

考虑极移的影响，卫星在协议地球坐标系中的坐标为

利用C语言程序实现

#include

#include

#include

#include

#define WE 7.292115e-6

struct canshu

{

int prn, nian, yue, ri, shi, fen;//卫星PRN号，年，月，日，时，分double miao;//秒

long double adoe, a0, a1, a2, mo, dn, e, ga, pio, io, w, pid, ii, cuc, cus, cue, crs, crc, cis, cic, toe, aodc, wn;

/*参数说明： ADOE值，a0 卫星钟偏差, a1 卫星钟漂移, a2 卫星钟频率漂移, M0 平近点角, Δn平运动差, e偏心率, a1/2半长轴的平方根, Ω0 轨道平面升交点经度,

i0 倾角, ω近地点角距, *Ω升交点速率, IDot 倾角速率, Cuc Cus 升交角距的摄动改正项, Crc Crs 地心距的摄动改正项, Cic Cis 倾角的摄动改正项,

toe 参考历元*/

};

void wxzbjx(struct canshu *pt)

{

long double a, n0, n, t, tk, toc, mk, ek, vk, fik, uk, rk, ik;

long double xk, yk ,zk, lk;

long double XK, YK, ZK;

int temp;

pt->nian = pt->nian + 2000;

t = (long double)(((pt ->nian)- 1980) * 365 * 24 * 3600 + (pt -

>yue - 1) * 30 * 24 * 3600 + pt ->ri

* 24 * 3600 + pt->shi * 3600 + pt ->miao);

a = pt ->ga * pt ->ga;

n0 = sqrt(WE/(a*a*a));//平均角速度n0

n = n0 + pt ->dn;

tk = t - pt ->toe;