GPS卫星轨道计算及其MATLAB仿真

GPS卫星轨道计算及其MATLAB仿真
黎奇，白征东，李帅，陈波波
（清华大学地球空间信息研究所，北京 100084）
一、程序设计思路
1. 读取RINEX文件（注意：文件路径）
2. 计算测量日周积秒（测量日的格里历→GPST）
3. 按卫星轨道计算步骤计算WGS-84坐标系坐标（内插）
4. 按需要将WGS-84坐标系下坐标转换为所需坐标系坐标
5. 画图输出
二、n 文件说明及读取程序
参考时刻oe t 的RINEX 格式的 “
”广播星历文件具体如下：（加粗部分为本次轨道
化Ω，
率i ，弧度/秒
4-22）
标svacc ，米
）
收到的卫星信号解，秒
）
文件名：RinexNreader.m 输 入：文件地址，卫星编号
三、计算测量日的周积秒
文件名：GCtoGPS.m （其中调用函数：GCtoJD.m）输入：指定公历的年、月、日、时、分、秒
文件名：GCtoJD.m
输 入：指定公历的 年、月、日
四、GPS 卫星轨道计算步骤及计算程序
1. 计算卫星运动的平均角速度n
平均角速度()
03
n =
经摄动参数n ∆改正后的平均角速度0n n n =+∆
3#61-79），n ∆（2#42-60）；14
3
2
3.98600510/GM m s =⨯ 2. 计算归化时间k t
说明：①广播星历是oe t 时刻的，对应的轨道参数也是oe t 时刻的，而观测时间在t 时刻，显然
oe t t ＜。

所以，要想获得t 时刻的轨道参数，需要知道t 与oe t 之间的差值即k t 。

以此，按照oe t 时刻
轨道参数，外推t 时刻轨道参数。

②k t 的起算时间是星期六/星期日子夜0点，当302400k t s ＞时，604800k t s -；当302400k t s ＜-时，+604800k t s 。

（604800s=1周） =k oe t t t -，且604800302400604800
302400
k k k k k k t t t t t t =-⎧⎨
=+
⎩＞ ＜-
已知：oe t （1#4-21）3. 计算观测时刻的平近点角k M
0k k M M nt =+
已知：0M （2#61-79），n （见1），k t （见2） 4. 计算观测时刻的偏近点角k E
sin k k k E M e E =+
已知：k M （见3），e （3#23-41）
方法：迭代解算，设初值0k k E M =，迭代2次基本收敛。

5. 计算观测时刻的真近点角k f
cos cos 1cos arctan cos sin 1cos k k k k k k
k k k E e f e E E f E e E f e E -⎫
=
⎪-⎪
⇒=⎬-⎪=
⎪-⎭ 已知：e （3#23-41），k E （见4） 6. 计算升交角距k Φ
k k f ωΦ=+
已知：k f 见（5），ω（5#42-60）
7. 计算因地球非球形引力和日月引力等因素而引起的升交角距u 的改正项u δ、卫星矢径r 的改正项r δ、轨道倾角i 的摄动改正项i δ
()()()()()()cos 2sin 2cos 2sin 2cos 2sin 2u uc k us k r rc k rs k i
ic k is k C C C C C C δδδ=Φ+Φ⎧⎪
=Φ+Φ⎨⎪=Φ+Φ⎩ 已知：uc C （3#4-22）、us C （3#42-60）、rc C （2#23-41）、rs C （5#23-41）、ic C （4#23-41）、is C （4#61-79），k Φ（见6）
8. 计算经改正后的升交角距、卫星矢径、轨道倾角
()01cos k k u k k r k i k u r a e E i i i t
δδδ=Φ+⎧⎪
=-+⎨⎪=++⋅⎩
已知：k Φ（见6），u δ、r δ、i δ（见7）
3#61-79），e （3#23-41），k E （见4），0i （5#4-22），i （6#4-22），k t （见2）
9. 计算卫星在轨道坐标系的位置(),,k k k x y z
cos sin 0
k k k
k k k k x r u y r u z =⎧⎪
=⎨⎪=⎩
已知：k u 、k r （见8）
10. 计算观测时刻t 的升交点赤径k Ω
()0k e k e oe t t ωωΩ=Ω+Ω--
已知：0Ω（4#42-60），Ω（5#61-79），地球自转角速度5
=7.292115146710e ω-⨯，k t （见2），
oe t （1#4-21）
11. 计算卫星在WGS-84坐标系的位置
cos cos sin sin cos cos sin k k k k k k k k
k k k k k k k X x y i Y x y i Z y i Ω-Ω⎡⎤⎡⎤⎢⎥⎢⎥=Ω+Ω⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦
已知：(),,k k k x y z （见9），k Ω（见10），k i （见8） 12. 计算卫星在协议地球坐标系中的空间直角坐标
10011P k p k P p k X X X Y Y Y Z X Y Z ⎡⎤⎡⎤⎡⎤
⎢⎥⎢⎥⎢⎥=-⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢
⎥-⎣⎦⎣⎦⎣⎦ 13. 注意：
① 以上计算，其时间均为GPST ；
② 卫星位置计算只考虑了地球自转的影响，没有顾及章动、岁差和极移的影响，因而是不严密的。

但对广播星历而言，其精度没有什么影响，而且计算简便。

文件名：orbitDetermine.m
输入：卫星编号data 、时间t
五、坐标转换程序设计
文件名：XYZtoBLH(X,Y,Z).m
输入：卫星WGS-84坐标系坐标x y z
六、运行绘图
运行命令：tracePlot。