目 录
21星历计算的时间和坐标系统
21.1 有关的时间系统与坐标系统
21.1.1 时间系统及其换算
41.1.2 坐标系统及其换算
71.2 计算单位和有关常数
122 轨道动力学计算的基本数学模型
122.1 二体问
题
快递公司问题件快递公司问题件货款处理关于圆的周长面积重点题型关于解方程组的题及答案关于南海问题
122.2 地球非球形引力摄动
152.3 日、月摄动
162.4 太阳直接辐射压摄动
192.5 地球固体潮摄动
192.6 大气阻力摄动
202.7 Y轴偏差加速度摄动
202.8 巡航姿态控制动力摄动
212.9 其它摄动影响
21附录：日月位置计算
243 轨道计算方法
243.1 Runge_Kutta积分法
253.2 Adams_Cowell积分
273.3 轨道计算
283.4 星历的快速插值
324 轨道根数与位置矢量、速度矢量的关系
324.1 由位置矢量和速度矢量计算轨道根数
334.2 由轨道根数计算位置矢量和速度矢量
1星历计算的时间和坐标系统
1.1 有关的时间系统与坐标系统
轨道计算过程重要涉及到不同的时间系统和坐标系统，下面将空间战场环境系统中所涉及到的时间系统和坐标系统进行定义，并说明各系统之间的相互关系。一般情况下，仿真系统采用的是TDT时间系统和J2000地心惯性坐标系。
1.1.1 时间系统及其换算
在轨道计算中，时间是独立变量。但是，在计算不同的物理量时，却使用不同的时间系统。例如：在计算恒星时用世界时UT1；定位解算时采用GPS时GPST；岁差和章动量的计算采用TDB时等。所以必须清楚各时间系统的定义和各时间系统之间的转换，下面给出各种时间系统的定义及它们之间的转换
公式
小学单位换算公式大全免费下载公式下载行测公式大全下载excel公式下载逻辑回归公式下载
。
格林尼治恒星时
格林尼治恒星时为春分点对格林尼治平天文子午面的时角。由于岁差、章动原因，它由格林尼治真恒星时（GAST）和平恒星时（GMST）之分。
两者的关系是：
其中：
为赤经章动
EMBED Equation.3 为自
起算至观测
时刻的儒略世纪数，即
世界时
是以平北极（国际习惯用原点）为统一
标准
excel标准偏差excel标准偏差函数exl标准差函数国标检验抽样标准表免费下载红头文件格式标准下载
的观测世界时，是反映地球实际自转的时间，恒星时计算与此有关。
国际原子时
时以铯原子
基态两能级间跃迁辐射的9192631770周所经历的时间作为1秒长的均匀时间，起点在1958年1月1日
。
国际协调时
是经跳秒修改后的国际原子时，它与世界时
的差
，观测纪录时间是以此为准的。
质心动力学时
（Barycentric Dynamical Time）
为相对于太阳质心的运动方程给出的历表、引数等所用的时间尺度，岁差及章动量的计算是以此为依据的。
地球动力学时
（Terrestrial Dynamical Time）
为视地心历表所用的时间尺度，它具有均匀连续的特性，卫星运动方程就是以此为独立的时间变量。
时间
是由系统定义和应用的一种时间尺度，于1980年1月6日
GPST与UTC相等，在此以后由系统主控站密切跟踪UTC以保持高度统一。但GPST不作跳秒修正，因此它与UTC具有整秒的差异（1997年1月至6月相差为
）。在计算GPS卫星轨道的初值时将涉及到GPST，GPS精密星历的参考时间为GPST。
以上各时间尺度的相互关系如下：
其中：
可从地球自转参数文件中获得；
；
。
上式中的T为自
年起算至观测TDB时刻的儒略世纪数，即：
不同时间系统间的关系如下图所示：
图1 几种时间系统之间的关系
1.1.2 坐标系统及其换算
卫星轨道计算和实际定位解算分别是在J2000.0惯性坐标系与WGS-84地固系中进行的，此外，卫星加速度计算中还涉及到星固坐标系。下面给出与本课题有关的主要坐标系的定义及相互转换关系。
地心惯性系
原点：地球质心
Z轴：向北指向
年平赤道面（基面）的极点
X轴：指向
平春分点
Y轴：符合右手系法则
位置矢量：
星固坐标系
原点：卫星质心
Z轴：指向卫星的天线方向，即指向地心
X轴：在轴与太阳构成的平面内，完成右手系法则
Y轴：沿卫星太阳能翼的支轴
位置矢量：
星固坐标系坐标轴
在惯性系中的方向余弦分别为：（
分别为太阳和卫星在
地心惯性系中的位置矢量）
WGS-84坐标系
WGS-84为1984年世界大地坐标系（WGS为World Geodetic System的简称），WGS-84的坐标定义及其采用的椭球参数为：
原点：地球质心
Z轴：指向BIH1984.0定义的协议地球极（CTP）方向
X轴：指向BIH1984.0的零子午面和CTP赤道的交点
Y轴：与X、Z轴成右手系
地球椭球长半径：ae=6378137 m
地球引力常数(含大气层)： GM=3986005(108 m3/s2
正常化二阶带球谐系数：
=-484.16685(10-6
地球自转角速度：(=7292115(10-11 rad/s2
地球椭球扁率：f = 1/298.257223563
地固坐标系与惯性坐标系的转换模型
a. 惯性坐标系
地固坐标系 模型
b. 地固坐标系
惯性坐标系 模型
式中：[A]为极移矩阵；
[B]为自转矩阵；
[C]为章动矩阵；
[D]为岁差矩阵。
上述各矩阵的意义及具体定义如下：
极移：由于地球不是刚体及其它一些地球物理因素的影响，地球自转轴相对于地球的位置随时间而变化从而引起观察者的天顶在天球上的位置发生变化，称为极移，矩阵为[A]：
式中：
为地极坐标，可从地球自转参数文件中给出的极移值插值得到。
自转：即地球公转的同时也在绕自转轴旋转。矩阵[B]：
式中：
为格林尼治恒星时
EMBED Equation.3
章动：是指外力作用下，地球自转轴在空间运动的短周期摆动部分，即同一瞬间真天极相对平天极的运动，月球对地球引力的变化是形成章动现象的主要外力作用，其次是太阳。矩阵[C]：
式中：
——交角章动
——黄经章动
其中：
都为常数，可自章动系数表1中查出。
而月亮平近点角：
太阳平近点角：
月亮平升交角：
日月平角矩：
月亮轨道对黄道平均升交点的黄经：
其中：
T0为自J2000.0起算至t的儒略世纪数
岁差：地球在太阳、月球和行星的引力作用下，地球的自转轴在空间不断发生变化，其长期运动称为岁差，矩阵[D]：
式中：
T0的意义同上。
旋转矩阵的求法分别为
1.2 计算单位和有关常数
轨道计算采用人卫单位系统，具体定义为：
长度单位：地球赤道半径
质量单位：地球的总质量
时间单位：
在以上人卫单位系统中，引力常数G=1。为完整起见，给出以下常数：
地球赤道半径：6378137m
地球扁率：
地球总质量：
地球自转角速度：
地心引力常数：
日心引力常数：
天文单位长度：
月球地球质量比：
光速：
引力常数：
太阳常数：
地球引力位系数采用WGS-84 EGM的规格化值。参见表2。
表1 黄经和倾角章动序列表
引 数
黄经章动
倾角章动
i
L
L’
F
D
0".0001
0".0001T
0".0001
0".0001T
1
0
0
0
0
1
-171996
-174.2
92025
8.9
2
0
0
0
0
2
2062
0.2
-895
0.5
3
-2
0
2
0
1
46
0
-24
0
4
2
0
-2
0
0
11
0
0
0
5
-2
0
2
0
2
-3
0
1
0
6
1
-1
0
-1
0
-3
0
0
0
7
0
-2
2
-2
1
-2
0
1
0
8
2
0
-2
0
1
1
0
0
0
9
0
0
2
-2
2
-13187
-1.6
5736
-3.1
10
0
1
0
0
0
1426
-3.4
54
-0.1
11
0
1
2
-2
2
-517
1.2
224
-0.6
12
0
-1
2
-2
2
217
-0.5
-95
0.3
13
0
0
2
-2
1
129
0.1
-70
0
14
2
0
0
-2
0
48
0
1
0
15
0
0
2
-2
0
-22
0
0
0
16
0
2
0
0
0
17
-0.1
0
0
17
0
1
0
0
1
-15
0
9
0
18
0
2
2
-2
2
-16
0.1
7
0
19
0
-1
0
0
1
-12
0
6
0
20
-2
0
0
2
1
-6
0
3
0
21
0
-1
2
-2
1
-5
0
3
0
22
2
0
0
-2
1
4
0
-2
0
23
0
1
2
-2
1
4
0
-2
0
24
1
0
0
-1
0
-4
0
0
0
25
2
1
0
-2
0
1
0
0
0
26
0
0
-2
2
1
1
0
0
0
27
0
1
-2
2
0
-1
0
0
0
28
0
1
0
0
2
1
0
0
0
29
-1
0
0
1
1
1
0
0
0
30
0
1
2
-2
0
-1
0
0
0
31
0
0
2
0
2
-2274
-0.2
977
-0.5
32
1
0
0
0
0
712
0.1
-7
0
33
0
0
2
0
1
-386
-0.4
200
0
34
1
0
2
0
2
-301
0
129
-0.1
35
1
0
0
-2
0
-158
0
-1
0
36
-1
0
2
0
2
123
0
-53
0
37
0
0
0
2
0
63
0
-2
0
38
1
0
0
0
1
63
0.1
-33
0
39
-1
0
0
0
1
-58
-0.1
32
0
40
-1
0
2
2
2
-59
0
26
0
41
1
0
2
0
1
-51
0
27
0
42
0
0
2
2
2
-38
0
16
0
43
2
0
0
0
0
29
0
-1
0
44
1
0
2
-2
2
29
0
-12
0
45
2
0
2
0
2
-31
0
13
0
46
0
0
2
0
0
26
0
-1
0
47
-1
0
2
0
1
21
0
-10
0
48
-1
0
0
2
1
16
0
-8
0
49
1
0
0
-2
1
-13
0
7
0
50
-1
0
2
2
1
-10
0
5
0
51
1
1
0
-2
0
-7
0
0
0
52
0
1
2
0
2
7
0
-3
0
53
0
-1
2
0
2
-7
0
3
0
54
1
0
2
2
2
-8
0
3
0
55
1
0
0
2
0
6
0
0
0
56
2
0
2
-2
2
6
0
-3
0
57
0
0
0
2
1
-6
0
3
0
58
0
0
2
2
1
-7
0
3
0
59
1
0
2
-2
1
6
0
-3
0
60
0
0
0
-2
1
-5
0
3
0
61
1
-1
0
0
0
5
0
0
0
62
2
0
2
0
1
-5
0
3
0
63
0
1
0
-2
0
-4
0
0
0
64
1
0
-2
0
0
4
0
0
0
65
0
0
0
1
0
-4
0
0
0
66
1
1
0
0
0
-3
0
0
0
67
1
0
2
0
0
3
0
0
0
68
1
-1
2
0
2
-3
0
1
0
69
-1
-1
2
2
2
-3
0
1
0
70
-2
0
0
0
1
-2
0
1
0
71
3
0
2
0
2
-3
0
1
0
72
0
-1
2
2
2
-3
0
1
0
73
1
1
2
0
2
2
0
-1
0
74
-1
0
2
-2
1
-2
0
1
0
75
2
0
0
0
1
2
0
-1
0
76
1
0
0
0
2
-2
0
1
0
77
3
0
0
0
0
2
0
0
0
78
0
0
2
1
2
2
0
-1
0
79
-1
0
0
0
2
1
0
-1
0
80
1
0
0
-4
0
-1
0
0
0
81
-2
0
2
2
2
1
0
-1
0
82
-1
0
2
4
2
-2
0
1
0
83
2
0
0
-4
0
-1
0
0
0
84
1
1
2
-2
2
1
0
-1
0
85
1
0
2
2
1
-1
0
1
0
86
-2
0
2
4
2
-1
0
1
0
87
-1
0
4
0
2
1
0
0
0
88
1
-1
0
-2
0
1
0
0
0
89
2
0
2
-2
1
1
0
-1
0
90
2
0
2
2
2
-1
0
0
0
91
1
0
0
2
1
-1
0
0
0
92
0
0
4
-2
2
1
0
0
0
93
3
0
2
-2
2
1
0
0
0
94
1
0
2
-2
0
-1
0
0
0
95
0
1
2
0
1
1
0
0
0
96
-1
-1
0
2
1
1
0
0
0
97
0
0
-2
0
1
-1
0
0
0
98
0
0
2
-1
2
-1
0
0
0
99
0
1
0
2
0
-1
0
0
0
100
1
0
-2
-2
0
-1
0
0
0
101
0
-1
2
0
1
-1
0
0
0
102
1
1
0
-2
1
-1
0
0
0
103
1
0
-2
2
0
-1
0
0
0
104
2
0
0
2
0
1
0
0
0
105
0
0
2
4
2
-1
0
0
0
106
0
1
0
1
0
1
0
0
0
表2 计算地球引力位加速度的引力位系数（GEM-T3）
2
0
-0.48416510e-3
0.0
2
1
-0.17e-9
1.19e-9
2
2
0.24390658e-5
-0.14000946e-5
3
0
0.95720109e-6
0.0
3
1
0.20277142e-5
0.24921712e-6
3
2
0.90447073e-6
-0.61944767e-6
3
3
0.72034249e-6
0.14138845e-5
4
0
0.53952118e-6
0.0
4
1
-0.53615108e-6
-0.47343598e-6
4
2
0.35021814e-6
0.66301523e-6
4
3
0.99093372e-6
-0.20092742e-6
4
4
-0.18877065e-6
0.30942370e-6
5
0
0.68343345e-7
0.0
5
1
-0.58280231e-7
-0.96083937e-7
5
2
0.65271099e-6
-0.32386369e-6
5
3
-0.45233006e-6
-0.21529578e-6
5
4
-0.29558409e-6
0.49690346e-7
5
5
0.17376349e-6
-0.66890704e-6
6
0
-0.14951352e-6
0.0
6
1
-0.76894161e-7
0.26998384e-7
6
2
0.48734482e-7
-0.37401311e-6
6
3
0.57203194e-7
0.93727543e-8
6
4
-0.86826474e-7
-0.47130637e-6
6
5
-0.26733038e-6
-0.53678021e-6
6
6
0.96846267e-7
-0.23713482e-6
7
0
0.91300885e-7
0.0
7
1
0.27486873e-6
0.97465920e-7
7
2
0.32779498e-6
0.93246712e-7
7
3
0.25122012e-6
-0.21529269e-6
7
4
-0.27556101e-6
-0.12376718e-6
7
5
0.13261698e-8
0.18620005e-7
7
6
-0.35883140e-6
0.15173875e-6
7
7
0.97027653e-9
0.24083642e-7
8
0
0.48883181e-7
0.0
8
1
0.23628239e-7
0.58847236e-7
8
2
0.77598534e-7
0.66008696e-7
8
3
-0.17785247e-7
-0.86346983e-7
8
4
-0.24633984e-6
0.70179584e-7
8
5
-0.25041147e-7
0.89426831e-7
8
6
-0.64923712e-7
0.30912257e-6
8
7
0.67462214e-7
0.75094831e-7
8
8
-0.12419836e-6
0.12017218e-6
2 轨道动力学计算的基本数学模型
卫星在轨道上运行要受到各种力因素的影响，产生的摄动是多方面的。国内外一些学者对卫星轨道的受摄问题作了详细的研究与分析，尤以澳大利亚的C.Rizos、A.Stolz和美国的H.F.Fliegel等人为代表。统筹考虑精度的需要和时间耗费，通过大量试算，本课题考虑了卫星所受的以下作用力来进行轨道计算（注：时间系统采用TDT时间系统、坐标系统采用J2000.0惯性坐标系），主要包括：地球质心引力F0、除质心外的地球引力FE、太阳和月球引力FN、太阳辐射压力FA、大气阻力（低卫星轨卫星）、Y轴偏差FY、地球潮汐附加力FT。
(2.1)
其中地球质心引力F0是最主要的，其次是地球的非质心引力FE，称为地球非球形摄动力。如果将地球质心引力视为1，地球非球形摄动力可达10-3量级，而其它摄动力则大多在10-6以下。
2.1 二体问题
在惯性系中，卫星运动方程被描述为
(2.2)
其中：
和
分别表示时刻卫星在惯性系中的位置、速度和加速度矢量；
和
为分别为引力常数和地球总质量。
（3.2）式的第一项为地球质心引力项，称为二体运动，是力模型的主项；第二项为摄动力引起的总摄动项，是
的
函数
excel方差函数excelsd函数已知函数 2 f x m x mx m 2 1 4 2拉格朗日函数pdf函数公式下载
矢量。
2.2 地球非球形引力摄动
在地固坐标系中，地球引力位函数作为拉普拉斯方程的解，其非球形部分为：
(2.3)
式中：
(2.4)
其中：λ和φ表示单位质点在地固坐标系中的地心经、纬度；
表示地球平均半径；
为规格化的缔合勒让德多项式；
和
为规格化的地球引力位系数；
n和m为多项式的阶和次，N为取的最高阶数。
非球形引力摄动的求解，按如下步骤进行：
a. 首先求解
和
。用下列公式逐阶次推算得到，即：
(2.5)
式中：
(2.6)
递推方法如下所示：
递推初值为：
(2.7)
递推方法采用先求对角线，再按行递推(m不变，n增加)进行。
如果采用非规格化的系数U和V进行计算，递推公式为：
(2.8)
递推初值为：
递推方法同上。此方法完成递推后，还需将U和V规格化，其公式为：
(2.9)
b. 计算
和
的偏导数
和
i. 首先计算
(2.10)
其中乘数因子为
(2.11)
ii. 计算
和
(2.12)
而：
c. 求非球形引力摄动引力位
EMBED Equation.3 (2.13)
则摄动加速度为：
(2.14)
式中：矩阵[W]为地固坐标系转换至惯性坐标系的旋转矩阵。
2.3 日、月摄动
考虑卫星的N体影响，只需顾及太阳和月球的引力作用已满足精度要求。N体摄动模型的建立是基于牛顿第二运动定律和万有引力定律。由图2所示，分析卫星的受力可得卫星
图2 三体的几何关系
的日、月摄动加速度矢量为：
(2.15)
其中：
，
为摄动体至卫星的中心距，即矢量
的模；
分别为摄动体
和卫星在惯性系中的位置矢量，
的计算参见本节附录。
这里S和L分别代表太阳和月球。
太阳和月球对卫星的摄动影响主要呈长周期变化，且与卫星轨道对太阳和月球的定向有关。
2.4 太阳直接辐射压摄动
照射在卫星上的太阳光，一部分被其吸收转化为热能，另一部分被反射向太空。因此，卫星会受到照射时的辐射力和反射时的反射力的作用，这里统称为直接辐射压摄动。直接辐射压摄动与光压强度、卫星表面的反射系数和光照面积有关。
由光压和牛顿第二运动定律建立的直接辐射压对卫星产生的运动加速度矢量为：
(2.16)
其中：
；
对应卫星表面反射系数项，是为补偿模型不足而引入的拟合参数；
；
为天文单位长度（地球轨道长半径）；
为人卫时间单位；
为卫星的日心距，即矢量
的模，
分别为卫星和太阳在惯性系中的位置矢量；
为蚀因子，具体定义为：
(2.17)
其中：
分别为太阳的被蚀视面积和视面积。要确定蚀因子，需计算下列诸量：
(2.18)
(2.19)
(2.20)
(2.21)
式中：
为月球和地球的视面积；
分别代表太阳和地球的扁率，
；
为考虑地球大气衰减以及地球扁率效应的有效地球赤道半径；
=1738000m为月球半径；
为考虑太阳扁率的太阳有效半径；
为地心纬度；
是地球—卫星—太阳张角；
是月球—卫星—太阳张角。
地影和月影判断过程如图3所示。
当卫星在地球半影中时：
其中：
当卫星在地球伪本影中时：
当卫星在月球半影中时：
其中：
当卫星在月球伪本影中时：
则蚀因子为：
图3 地影和月影判断流程
太阳直接辐射摄动对卫星轨道的影响是十分复杂的，它与卫星表面的反射特性、卫星轨道相对太阳的定向以及太阳活动等有关。卫星是由各种不同折射性质的原材料构成的不规则形体，在其运行过程中，太阳对它照射的面积也在不断地改变（它的太阳能翼始终是面向太阳的）。此外，由于太阳活动的变化，所谓太阳常数
也并非常数。因此，给出卫星受太阳直接辐射压摄动的精确模型是很困难的。所以采用较简单的平面模型计算太阳辐射加速度影响。
2.5 地球固体潮摄动
地球并非刚体，它受日月引力作用会产生弹性形变，称为潮汐现象。这种形变使得地球内部物质发生小的变化，随之导致引力位函数产生小的形变位差——潮汐位。卫星运动的地球固体潮摄动就是潮汐位效应的结果。
已知日（或月）对地面点的引力位球谐展开式为：
式中：
为日（或月）的质量，
分别为引力体至地心距和与地面点的地心角，
为n阶勒让德多项式。
从上式中排除不产生形变位差的0或1阶项，且只取n=2阶项，可得日月潮汐形变对卫星产生的摄动位为：
其中：
为二阶Love数（取
）。
顾及关系式
，最后得到卫星的固体潮摄动加速度矢量为：
(2.22)
式中：
分别为
的单位矢量。
2.6 大气阻力摄动
高层大气对卫星的运行将产生阻力，这种阻力对低轨道卫星是主要摄动力之一，但大气密度变化机制非常复杂，不但随高度变化，也与太阳活动、时间、季节、纬度和地磁活动的变化有关。本文采用了静止球型大气密度模型（HP模型）。若只考虑大气分子对卫星表面的法向作用力而忽略其切向作用力，大气阻力使卫星产生的摄动加速度为:
EMBED Equation.3
EMBED Equation.3 (32.23)
其中：
为卫星星体部分的大气阻力摄动加速度；
为卫星太阳帆板的大气阻力摄动加速度；
为大气阻力系数；
是大气密度，由HP模型计算得到；
是卫星相对大气的速度矢量
；
是卫星参考面积与卫星质量之比；
为太阳帆板的大气阻力系数；
为太阳帆板的面积；
为太阳帆板的法向与卫星相对于大气的速度方向的夹角；
为
向量的单位向量。
2.7 Y轴偏差加速度摄动
Y轴偏差加速度主要是GPS卫星的结构失调和卫星体的热辐射而产生的一个附加加速度在星固坐标系Y轴方向上的分量。在
设计
领导形象设计圆作业设计ao工艺污水处理厂设计附属工程施工组织设计清扫机器人结构设计
上，为使卫星的太阳翼以最大面积朝向太阳，两翼的支轴应保持在一条直线上，并要求垂直于卫星和太阳方向的连线，用于控制两翼俯仰的太阳传感器也应完全垂直于卫星翼支轴，卫星的偏航高度控制应保持正确，但事实上并不完全这样；另一方面，由卫星本身产生的超高温要从Y轴方向的百叶孔排出，这样使处在不稳定状态中的卫星体也会产生结构失调现象。由此导致了Y轴偏差加速度的存在。H.F.Fliegel等人把结构失调的影响表示为[14]：
(2.24)
其中：
为考虑模型剩差而引入的待估参数；
为星固坐标系的Y轴在J2000.0惯性坐标系中的方向余弦，由下式决定：
2.8 巡航姿态控制动力摄动
有些卫星在巡航过程中需要保持三轴稳定姿态，需通过姿态控制实现，有的卫星其姿态控制的动力来源于高压气瓶的喷气。这样，在姿态控制的同时也影响了卫星质心的运动。该摄动加速度矢量可以表示为:
(2.25)
其中：
；
为卫星在RTN坐标系中姿控动力摄动加速度的常数分量；
为卫星在RTN坐标系中姿控动力摄动加速度的时间变化率；
和
为周期项的系数；
对全局参数，为由历元时刻起算的相对时间；对弧段相关参数，为观测时刻t在相应弧段内的相对时间。
2.9 其它摄动影响
在轨道运行的卫星除受到上述摄动作用外，还受其它一些摄动的影响，如反照辐射摄动、地球自转形变摄动、广义相对论摄动、海潮摄动、大气潮摄动等。这些摄动影响对卫星轨道摄动非常小，但其计算却要耗费大量的时间片。考虑到2.1节所述的摄动已能满足课题的精度需求和时间消耗的限制，因此，本课题中忽略了其它摄动的影响。
附录：日月位置计算
a. 太阳位置矢量
计算
其中：[D]T为计算历元时刻的平春分点到J2000.0平春分点的转换矩阵，即岁差矩阵。
计算当时平春分点的几何平黄经为：
平近点角为：
平黄赤交角为：
以上式中
b. 月球位置矢量
计算
其中：
系数
、
、
和
的值以及幅角
的计算列入表5中，表内D为日月平角距，
，其计算式为：
表5 计算月球位置的系数表
j
1
0.022230
0.010020
0.000820
0.000000
2D-M'
2
0.011610
0.008250
-0.002300
-0.000190
2D
3
0.003240
0.000120
0.002830
0.000000
M+180°
4
0.000000
0.000000
-0.002130
-0.000190
2F+180°
5
0.001030
0.000090
0.000450
0.000000
2M'-2D+180°
6
0.001000
0.000420
0.000000
0.000000
2D-M'-M
7
0.000930
0.000900
-0.000080
0.000000
2D+M'
8
0.000800
0.000560
-0.000120
0.000000
2D-M
9
0.000720
0.000340
0.000000
0.000000
M'-M
10
0.000610
0.000290
0.000000
0.000000
D+PI
11
0.000530
0.000280
0.000000
0.000000
M'+M+180°
12
0.000270
0.000000
-0.027240
-0.002410
2F-2D+180°
13
0.000000
0.000000
-0.000080
0.000000
2F+M'+180°
14
0.000410
0.000200
0.000700
0.000000
2F-M'+180°
15
0.000180
0.000180
0.000000
0.000000
4D-M'
16
0.000150
0.000120
0.000000
0.000000
4D-2M'
17
0.000140
0.000070
0.000000
0.000000
2D+M-M'+180°
18
0.000120
0.000090
0.000000
0.000000
2D+M+180°
19
0.000090
0.000000
0.000000
0.000000
M'-D
20
0.000090
0.000000
0.000000
0.000000
D+M
21
0.000070
0.000070
0.000000
0.000000
2D-M+M'
22
0.000070
0.000090
0.000000
0.000000
2D+2M'
23
0.000070
0.000080
0.000000
0.000000
4D
24
0.000000
0.000000
0.000190
0.000180
2F-2M'
25
0.000000
0.000000
0.001100
0.000100
2F-2D+M
26
0.000000
0.000000
0.000470
0.000000
2D-2F+M
27
0.000000
0.000000
0.000230
0.000000
2F-2D-M'
28
0.000000
0.000000
0.000023
0.000000
2D-2F-M'
3 轨道计算方法
卫星运动方程的解有分析法、数值法和半分析半数值方法等。分析法是将力模型展开取有限项，给出任意时刻解的表达式，通常称之为一般摄动法。其优点是便于定性地给出轨道的特征，如长期项、长周期项影响等。但对摄动力模型处理限制较大，使精度受到影响。数值法即常微分方程的初值解问题，在轨道计算中称之为特别摄动法。其优点是能完整地顾及所选择的力模型、处理简单、计算精度高，是高精度卫星轨道计算最实际有效的方法。
经研究，Runge_Kutta型和Cowell型数值积分方法都可用于产生高精度的卫星轨道，尤其是后者通过预报——校准公式进行迭代计算，收敛快（一般只需迭代2~3次），累积误差影响要比前者小的多。所以这一方法是轨道计算最常用的方法之一。不足的是它为多步型积分器，必须用其它方法为它准备起步值。
本课题提供的积分模式有Runge_Kutta 4、Runge_Kutta 8和高精度的Adams_Cowell方法，高精度计算时采用Runge_Kutta单步积分起步，用Adams_Cowell多步积分求解卫星运动方程，确定卫星在J2000.0惯性坐标系中的位置和速度矢量。
图4 Runge_Kutta 积分框图
3.1 Runge_Kutta积分法
Runge_Kutta法是一种单步积分方法，公式形式简单，应用方便，将它用于解算卫星运动方程和变分方程也具有很好的稳定性。但该方法随着积分式阶数的增高，计算右函数的次数也会随之增多，无疑会导致累积误差增加，同时也会降低计算效率。因此，在大规模的轨道计算中，通常只将其作为多步型积分器的起步方法。流程图如图4所示。
设有初值问题：
(3.1)
则求解该问题的Runge_Kutta 积分公式为：
(3.2)
式中：
，为积分步长；
k为积分式所取的阶数；
,
均为已知的系数，具体参见表6所示。
由运动方程和变分方程形成的矩阵常微分方程的初值问题为：
上式中：左端
；
右函数
；
初值
3.2 Adams_Cowell积分
Adams_Cowell积分适用显含1阶导的二阶微分方程，初始值需要位置项和速度项。积分的预报公式为：
(3.3)
其中：
为Adams_Cowell积分公式的预报系数，按下列公式计算，结果参见表4.2。
校正公式为：
(3.4)
式中：
为Adams_Cowell积分公式的校正系数，按下列公式计算，结果参见表4.2。
(3.3)、(3.4)式中的
和
按以下公式递推计算：
和
分别为
的一次和分和二次和分，按下列公式递推计算：
和
的初值定义为：
用Adams_Cowell积分解卫星运动方程积分流程图如图5所示：
图5 Adams_Cowell积分流程图
3.3 轨道计算
预测卫星轨道实际上就是由一组正确的初始值数值积分卫星运动方程，外推出未来的卫星轨道。显然，在卫星运动力模型较完备的情况下，初始值的精确与否将直接影响外推轨道的精度。且在推算过程中，由于受初始值误差传播影响，推算区间越长，积分累积的误差就越大。轨道计算步骤为：
图6 卫星轨道计算流程图
3.4 星历的快速插值
通过分布式计算提供的星历数据间隔太大，而在实际定位计算中，有时需要间隔为1秒甚至步长更小的精确星历坐标。因此，高精度、快速的精密星历插值就成为一项重要工作。以精度和计算耗费小为出发点，本课题在研究的过程中，分析和比较了以下几种多项式插值方法：拉格朗日多项式插值、牛顿多项式插值、三角函数插值、样条函数插值、切比雪夫多项式插值等。
拉格朗日多项式插值的多项式构造和插值时间耗费非常大。其时间耗费可由下式表示：
其中：
表示已知点数；
A表示加法运算、M和D分别表示乘、除运算。
构造插值多项式还需额外的
运算。因此，总运算耗费为：
比拉格朗日多项式插值更有效的方法是利用数据均差的牛顿多项式插值[1][2]。假设有
个控制点
，定义在
处的零阶均差为
，在
处的1阶均差定义为
，k阶均差定义为：
则牛顿插值多项式可表示为：
(4.12)
上述表达式在每级中都由2个加法和1个乘法运算组成，其时间耗费为：
。与经典的拉格朗日多项式插值相比，牛顿多项式插值经济得多。给定
个不同的点
和相应的均差系数
，对任意时刻
插值多项式的结果由下列迭代产生（Horner算法）。
上述算法的精度和拉格朗日的精度相同，但时间耗费得到了明显的改善，而且各阶均差可以预先运算好储存起来，需要时直接调用即可。
插值点等间隔条件下的插值是牛顿多项式插值的一种特例，可以导出更简洁的形式，而不降低精度。IGS提供的精密星历数据是等间隔的，步长为15分钟。在实际应用中，我们采用牛顿前向差分形式。
牛顿前向均差的定义为：
同理可推出各高阶前向均差为：
(4.13)
假设插值多项式的阶数为
，则插值多项式为：
(4.14)
令
，则上式可表示为：
，
用(4.13)式构造插值多项式，再用(4.14)式用作插值，该方法在时间耗费上比Horner算法有一定的改善。
与经典的拉格朗日插值方法相比，三次样条多项式插值速度快，但精度低，三角函数多项式插值必须用DFT求解系数，时间耗费大，且这两种方法都将星历数据近似为周期函数，插值精度低，不能满足高精度插值的需要。切比雪夫多项式插值方法在端点处抑制了误差的扩大，但计算量较大。牛顿多项式插值方法的精度高，速度快，充分利用了星历数据的特点，且克服了拉格朗日方法的不足。在精密星历插值上，牛顿多项式插值方法比拉格朗日多项式插值方法更为有效。由拉格朗日插值多项式和牛顿插值多项式的余项公式可推出截断误差在端点处将迅速增长。为克服这种不足、提高插值的精度，在实现时提出并采用了“滑动式牛顿多项式插值”方法[19]，即采用滑动方式进行插值：每次取13个点，生成12阶多项式，始终让被插值点在已知点的第6个点到第7个点之间。第1个点到第13个点为第一个插值区间，仅用来插值第6个点到第7个点之间的时间段；第2个点到第14个点为第二个插值区间，仅用来插值第7个点到第8个点之间的时间段；依次类推下去。这种滑动式牛顿多项式插值精度非常高，可达到cm级精度，且易于实现，时间耗费小，确保了精密轨道解算时的精度和实时要求。
表6 8阶Runge_Kutta积分公式系数表
1
2
3
4
5
6
7
8
9
10
1
0
2
4
3
1
3
4
1
0
3
5
1
0
0
3
6
13
0
-27
42
8
7
389
0
-54
966
-824
243
8
1
-231
0
81
-1164
656
-122
800
9
-127
0
18
-678
456
-9
576
4
10
1
1481
0
-81
7104
-3376
72
-5040
-60
720
41
0
0
27
272
27
216
0
216
41
表7 12阶Adams_Cowell积分公式系数表
0
4.259634
0.877861
0.264351
0.054794
1
-22.691026
-4.800775
0.823067
0.165534
2
80.019378
16.861107
-2.071621
-0.526813
3
-196.294318
-41.187179
4.414893
1.189914
4
349.412940
73.086989
-7.283104
-2.009198
5
-462.480193
-96.516821
9.192754
2.566623
6
460.087276
95.852539
-8.853279
-2.489666
7
-343.735149
-71.516859
6.460362
1.825385
8
190.394308
39.571065
-3.514964
-0.996493
9
-75.976728
-15.777189
1.383094
0.393084
10
20.680772
4.291462
-0.372243
-0.105997
11
-3.441245
-0.713663
0.061367
0.017501
12
0.264351
0.054794
-0.004677
-0.001336
4 轨道根数与位置矢量、速度矢量的关系
4.1 由位置矢量和速度矢量计算轨道根数
已知某时刻
的位置矢量
和速度矢量
，计算轨道根数的步骤如下：
由
（5.1）
其中
，分别计算
和
（5.2）
而
（5.3）
计算
（5.4）
计算
（5.5）
类似可得
(5.6)
所以：
(5.7)
计算
，由
得
(5.8)
4.2 由轨道根数计算位置矢量和速度矢量
若已知任何时刻
时天体的
和
，在轨道坐标系
中有
（5.9）
根据惯性坐标系与个坐标系的转换关系可得位置矢量的表达式
将旋转矩阵展开可得
(5.10)
其中
和
分别表示
轴和
的单位矢量，即
(5.11)
(5.12)
速度矢量表达式为：
(5.13)
� EMBED Equation.3 ���(� EMBED Equation.3 ���(� EMBED Equation.3 ���(� EMBED Equation.3 ���(� EMBED Equation.3 ���(� EMBED Equation.3 ���
(
0 � EMBED Equation.3 ���(� EMBED Equation.3 ���(� EMBED Equation.3 ���(� EMBED Equation.3 ���(� EMBED Equation.3 ���
(
0 � EMBED Equation.3 ���(� EMBED Equation.3 ���(� EMBED Equation.3 ���(� EMBED Equation.3 ���
(
0 � EMBED Equation.3 ���(� EMBED Equation.3 ���(� EMBED Equation.3 ���
(
0 � EMBED Equation.3 ���(� EMBED Equation.3 ���
( (
� EMBED Equation.3 ���
17
_982483215.unknown
_982558933.unknown
_983446077.unknown
_1112723604.unknown
_1115617741.unknown
_1143373967.unknown
_1143374474.unknown
_1143374840.unknown
_1143375363.unknown
_1143375519.unknown
_1143375634.unknown
_1143374907.unknown
_1143375188.unknown
_1143374863.unknown
_1143374742.unknown
_1143374760.unknown
_1143374599.unknown
_1143374224.unknown
_1143374440.unknown
_1143374447.unknown
_1143374420.unknown
_1143374081.unknown
_1143374167.unknown
_1143374046.unknown
_1143373177.unknown
_1143373607.unknown
_1143373752.unknown
_1143373774.unknown
_1143373694.unknown
_1143373372.unkno
浙公网安备 33021202002483