C语言计算GPS卫星位置
1 概述
在用GPS信号进行导航定位以及制订观测计划时，都必须已知GPS卫星在空间的瞬间位置。卫星位置的计算是根据卫星电文所提供的轨道参数按一定的公式计算的。本节专门讲解观测瞬间GPS卫星在地固坐标系中坐标的计算方法。
2 卫星位置的计算
1. 计算卫星运行的平均角速度n
根据开普勒第三定律，卫星运行的平均角速度n0可以用下式计算：
n0?GM/a3??(a)3?4?11?
143
2
式中μ为WGS-84坐标系中的地球引力常数，且μ=3.986005×10m/s。平均角速度n0加上卫星电文给出的摄动改正数Δn，便得到卫星运行的平均角速度n
n=n0+Δn (4-12)
2. 计算归化时间tk
首先对观测时刻t′作卫星钟差改正 t=t′-Δt
?t?a0?a1(t'?toc)?a2(t'?toc)2
然后对观测时刻t归化到GPS时系
tk=t-toc (4-13)
式中tk称作相对于参考时刻toe的归化时间（读者注意：toc≠toe）。 3. 观测时刻卫星平近点角Mk的计算
Mk=M0+ntk (4-14)
式中M0是卫星电文给出的参考时刻toe的平近点角。 4. 计算偏近点角Ek
Ek=Mk+esinEk(Ek,Mk以弧度计) （4-15）
上述方程可用迭代法进行解算，即先令Ek=Mk,代入上式，求出Ek再代入上式计算，因为GPS卫星轨道的偏心率e很小，因此收敛快，只需迭代计算两次便可求得偏近点角Ek。 5. 真近点角Vk的计算
由于:
cosVk?(cosEk?e)(1?ecosEk)?4?16?
sinVk?arctg?1?e2?sinEk(cosEk?e)?????因此:
Vk?arctg?(1?e2)?sinEk(cosEk?e)??????4?17?
?4?18?
6．升交距角Φk的计算
?k?Vk??ω为卫星电文给出的近地点角距。
7. 摄动改正项δu,δr,δi的计算
(4?19)
?u?Cuc?cos(2?k)?Cus?sin(2?k)???rr?Crc?cos(2?k)?Crs?sin(2?k)??ri?Cic?cos(2?k)?Cis?sin(2?k)??(4?20)
δu,δr,δi分别为升交距角u的摄动量，卫星矢径r的摄动量和轨道倾角i
的摄动量。
8. 计算经过摄动改正的升交距角uk、卫星矢径rk和轨道倾角ik
uk??k??u??rk?a(1?ecosEk)??r??ik?i0??i?Itk?9. 计算卫星在轨道平面坐标系的坐标
(4?21)
卫星在轨道平面直角坐标系（X轴指向升交点）中的坐标为
xk?rkcosuk??yk?rksinuk?10. 观测时刻升交点经度Ωk的计算
(4?22)
升交点经度Ωk等于观测时刻升交点赤经Ω（春分点和升交点之间的角距）与格林泥治视恒星时GAST（春分点和格林尼治起始子午线之间的角距）之差，
Ωk=Ω-GAST （4-23） 又因为:
????oe??tk (4-24)
其中Ωoe为参与时刻toe的升交点的赤经；
?是升交点赤经的变化率，卫星电文每小时更新一次Ω和toe。
?此外，卫星电文中提供了一周的开始时刻tw的格林尼治视恒星时GASTw。由于地球自转作用，GAST不断增加，所以:
GAST=GASTw+ωet （4-25）
式中ωe=7.29211567×10-5rad/s为地球自转的速率；t为观测时刻。 由式（4-24）和（4-25），得:
?k??oe?GASTw??et由（4-13）式，得:
(4-26)??
?k??0?(???e)tk??etoe其中
(4-27)
?0??oe?GASTw，?o、?、toe的值可从卫星电文中获取。
11. 计算卫星在地心固定坐标系中的直角坐标
把卫星在轨道平面直角坐标系中的坐标进行旋转变换，可得出卫星在地心固定坐标系中的三维坐标：
?Xk??xkcos?k?ykcosiksin?k??Y???xsin??ycosicos??kkik??k??k???Zk????yksinik?12. 卫星在协议地球坐标系中的坐标计算
考虑极移的影响，卫星在协议地球坐标系中的坐标为
(4-28)
? 1 0 Xp??Xk??X??Y??? 0 1 Y??Y??p??k?????Z???Z??CTS? Xp Yp 1???k?
(4?29)
利用C语言程序实现
#include
#include
#include
#include
#define u 3.986004418e+14
#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 倾角的摄动改正项, };
void wxzbjx(struct canshu *pt) {
pt->nian = pt->nian + 2000;
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; toe 参考历元*/
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;
toc = pt ->a0 + pt ->a1 * (t - pt ->toe) + pt ->a2 * (t - pt->toe) * (t - pt ->toe);
tk = tk - pt ->toe; mk = pt ->mo + n * tk; ek = mk;
for(temp=0;temp<10;temp++) { }
vk = 2 * atan(sqrt((1+ pt ->e) / (1 - pt ->e))* (tan(ek)) / 2 ); fik = vk + pt ->w;
uk = fik + pt ->cuc * cos(2* fik) + pt ->cus * sin(2*fik);
ek=mk+ pt ->e * sin(ek);//利用迭代法求偏近点角ek
rk = pt ->ga * pt ->ga * (1 - pt ->e * ek) + pt ->crc * cos(2* fik) + pt ->crs * sin(2*fik);
printf(\年%d月%d号%d时%.2f秒 %d号卫星的坐标:\>yue ,pt->ri , pt->shi ,pt->miao, pt->prn);
ik = pt ->io + pt ->cic * cos(2 * fik) + pt ->cis * sin(2* fik) + pt ->ii * tk;
xk = rk * cos(uk); yk = rk * sin(uk); zk = 0;
lk = pt ->pio + (pt ->pid - WE) * tk - WE * pt->toe;
XK = xk * cos(lk) - yk * cos(ik) * sin(lk); YK = xk * sin(lk) + yk * cos(ik) * cos(lk); ZK = yk * sin(ik);
printf(\
最新c语言计算gps卫星位置演示教学
