卫星轨道动力学数值计算.

目 录

1星历计算的时间和坐标系

统 . ....................................................................................................... 2

1.1 有关的时间系统与坐标系

统 . ........................................................................................... 2

1.1.1 时间系统及其换

算 . ................................................................................................ 2

1.1.2 坐标系统及其换

算 . ................................................................................................ 4

1.2 计算单位和有关常

数 . ....................................................................................................... 7

2 轨道动力学计算的基本数学模

型 . ............................................................................................ 12

2.1 二体问

题 . ................................................................................................................. 12

2.2 地球非球形引力摄

动 . ............................................................................................. 12

2.3 日、月摄

动 . ............................................................................................................. 15

2.4 太阳直接辐射压摄

动 . ............................................................................................. 16

2.5 地球固体潮摄

动 . ..................................................................................................... 19

2.6 大气阻力摄

动 . ......................................................................................................... 19

2.7 Y轴偏差加速度摄

动 ............................................................................................... 20

2.8 巡航姿态控制动力摄

动 . ......................................................................................... 20

2.9 其它摄动影

响 . ......................................................................................................... 21

附录：日月位置计

算 . .................................................................................................... 21

3 轨道计算方

法 ............................................................................................................................ 24

3.1 Runge_Kutta积分

法 ........................................................................................................ 24

3.2 Adams_Cowell积

分 ......................................................................................................... 25

3.3 轨道计

算.......................................................................................................................... 27

3.4 星历的快速插

值 . ............................................................................................................. 28

4 轨道根数与位置矢量、速度矢量的关

系 . ................................................................................ 32

4.1 由位置矢量和速度矢量计算轨道根数 . ......................................................................... 32

4.2 由轨道根数计算位置矢量和速度矢量 . ......................................................................... 33

1星历计算的时间和坐标系统 1.1 有关的时间系统与坐标系统

轨道计算过程重要涉及到不同的时间系统和坐标系统，下面将空间战场环境系统中所涉及到的时间系统和坐标系统进行定义，并说明各系统之间的相互关系。一般情况下，仿真系统采用的是TDT 时间系统和J2000地心惯性坐标系。

1.1.1 时间系统及其换算

在轨道计算中，时间是变量。但是，在计算不同的物理量时，却使用不同的时间系统。例如：在计算恒星时用世界时UT1；定位解算时采用GPS 时GPST ；岁差和章动量的计算采用TDB 时等。所以必须清楚各时间系统的定义和各时间系统之间的转换，下面给出各种时间系统的定义及它们之间的转换公式。

格林尼治恒星时

格林尼治恒星时为春分点对格林尼治平天文子午面的时角。由于岁差、章动原因，它由格林尼治真恒星时（GAST ）和平恒星时（GMST ）之分。

两者的关系是： εψcos ∆+=GMST GAST 其中：εψcos ∆为赤经章动

3521062. 0093104. 0 876600812866. 80184(841. 67310u s u s u h s s T T T GMST -⨯-+++=u T 为自 0. 24515(0. 2000

JD J 起算至观测1UT 时刻的儒略世纪数，即 0 . 365250. 24515 1(-=UT JD T u 世界时1UT

1UT 是以平北极（国际习惯用原点）为统一标准的观测世界时，是反映地球实际自转的时间，恒星时计算与此有关。

国际原子时TAI

TAI 时以铯原子133s C 基态两能级间跃迁辐射的9192631770周所经历的时间作为1秒长的均匀时间，起点在1958年1月1日UT h

0。

国际协调时UTC

UTC 是经跳秒修改后的国际原子时，它与世界时1UT 的差9. 0s ，观测纪录时间是以此为准的。

质心动力学时TDB （Barycentric Dynamical Time）

TDB 为相对于太阳质心的运动方程给出的历表、引数等所用的时间尺度，岁差及章动量的计算是以此为依据的。

地球动力学时TDT （Terrestrial Dynamical Time）

TDT 为视地心历表所用的时间尺度，它具有均匀连续的特性，卫星运动方程就是以此为的时间变量。

GPS 时间GPST

GPST 是由系统定义和应用的一种时间尺度，于1980年1月6日h 0GPST 与UTC 相等，

在此以后由系统主控站密切跟踪UTC 以保持高度统一。但GPST 不作跳秒修正，因此它与UTC 具有整秒的差异（1997年1月至6月相差为s

11）。在计算GPS 卫星轨道的初值时将涉及到GPST ，GPS 精密星历的参考时间为GPST 。

以上各时间尺度的相互关系如下： GPST UTC GPST TD TDT TDB TAI TDT AT

UTC TAI UT UTC UT s ∆+=∆+=+=∆+=∆+=184. 321 1

其中：1UT ∆可从地球自转参数文件中获得； GPST AT s ∆+=∆19；

(3019501. 628240040768. 6, sin 0167. 0sin(001658. 0rad T v v v TD s +=+=∆。 上式中的T 为自0. 2000J 年起算至观测TDB 时刻的儒略世纪数，即： 0

. 365250. 24515 (-=TDB JD T 不同时间系统间的关系如下图所示：

图1 几种时间系统之间的关系 1.1.2 坐标系统及其换算

卫星轨道计算和实际定位解算分别是在J2000.0惯性坐标系与WGS-84地固系中进行的，此外，卫星加速度计算中还涉及到星固坐标系。下面给出与本课题有关的主要坐标系的定义及相互转换关系。

0. 2000J 地心惯性系 原点：地球质心

Z 轴：向北指向0. 2000J 年平赤道面（基面）的极点 X 轴：指向0. 2000J 平春分点 Y 轴：符合右手系法则 位置矢量：r 星固坐标系 原点：卫星质心

Z 轴：指向卫星的天线方向，即指向地心

X 轴：在轴与太阳构成的平面内，完成右手系法则

Y 轴：沿卫星太阳能翼的支轴 位置矢量：a r

星固坐标系坐标轴 ˆ, ˆ, ˆ(a a a Z Y X 在惯性系中的方向余弦分别为：（r r s , 分别为太阳和卫星在0. 2000J 地心惯性系中的位置矢量） r r r r Z r r Y Y r Y r X s s a s s a a a a -=∆⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧-=⨯∆⨯∆=⨯⨯=, ˆˆˆ WGS-84坐标系

WGS-84为1984年世界大地坐标系（WGS 为World Geodetic System的简称），WGS-84的坐标定义及其采用的椭球参数为：

原点：地球质心

Z 轴：指向BIH1984.0定义的协议地球极（CTP ）方向 X 轴：指向BIH1984.0的零子午面和CTP 赤道的交点 Y 轴：与X 、Z 轴成右手系 地球椭球长半径：a e =6378137 m

地球引力常数(含大气层 ： GM=3986005⨯108 m3/s2 正常化二阶带球谐系数：0. 2=-484.16685⨯10-6 地球自转角速度：ω=7292115⨯10-11 rad/s 2 地球椭球扁率：f = 1/298.257223563 地固坐标系与惯性坐标系的转换模型

a. 惯性坐标系−−−→−转换到地固坐标系 模型 0 . 20000. 2000][]][][][[J J D Z Y X W Z Y X D C B A Z y X

⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡ b. 地固坐标系−−−→−转换到

惯性坐标系 模型 D

T D T T T T J Z Y X W Z Y X A B C D Z y X

⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡][][][][][0. 2000 式中：[A]为极移矩阵；

[B]为自转矩阵； [C]为章动矩阵； [D]为岁差矩阵。

上述各矩阵的意义及具体定义如下：

极移：由于地球不是刚体及其它一些地球物理因素的影响，地球自转轴相对于地球的位置随时间而变化从而引起观察者的天顶在天球上的位置发生变化，称为极移，矩阵为[A]：

( (][p X p Y y R x R A --=

式中：p p y x , 为地极坐标，可从地球自转参数文件中给出的极移值插值得到。 自转：即地球公转的同时也在绕自转轴旋转。矩阵[B]：

(][G θZ R B =

式中：G θ为格林尼治恒星时

cos(1062. 0093104. 0 876600812866. 80184(841. 67310352εεϕθ∆+∆+⨯-+++=-M u s u s u

h s s G T T T 0. 365250. 24515 1(-=UT JD T u 章动：是指外力作用下，地球自转轴在空间运动的短周期摆动部分，即同一瞬间真天极相对平天极的运动，月球对地球引力的变化是形成章动现象的主要外力作用，其次是太阳。矩阵

[C]：

( ( (][M X Z M X R R R C εϕεε∆-∆--= 式中：

30200001813'. ' 000059'. ' 08150'. ' 448'. ' 21' 2623T T T M +--= ε ∑∑===∆1061 5

1 cos(j k k jk j A n d ε——交角章动 ∑∑===∆106 1 5 1

sin(j k k jk j A n c ϕ——黄经章动

其中：jk j j n d c , , 都为常数，可自章动系数表1中查出。 而月亮平近点角： 3020010'. ' 0310'. ' 31 633'. ' 7159221325(733'. ' 485866T T T A r ++++= 太阳平近点角：302002012'. ' 0577'. ' 0 224'. ' 129258199(804'. ' 1287099T T T A r --++=

月亮平升交角：30200'

' 3011'. ' 0257'. ' 13 137'. 2952631342(877. 335778T T T A r +-++= 日月平角矩：30200' ' 4019'. ' 01'. ' 6 328'. ' 11056011236(307. 1072261T T T A r +-++=

月亮轨道对黄道平均升交点的黄经： 30200'

' 5008'. ' 0455'. ' 7 539'. ' 48205(280. 450160T T T A r +++-= 其中： 3601=r

T0为自J2000.0起算至t 的儒略世纪数 0 . 365255 . 514 (0-= TDB MJD T

岁差：地球在太阳、月球和行星的引力作用下，地球的自转轴在空间不断发生变化，其长期运动称为岁差，矩阵[D]：

( ( (][p Z p Y p Z R R Z R D ξθ--= 式中：

30200017998'. ' 030188'. ' 02181'. ' 2306T T T p ++=ξ 30200041833'. ' 042665'. ' 03109'. ' 2004T T T p --=θ

30200018203'. ' 009468'. ' 12181'. ' 2306T T T Z p ++= T 0的意义同上。 旋转矩阵的求法分别为 ⎥⎥⎥⎦

⎤

⎢⎢⎢⎣⎡-=θθθθ θcos sin 0sin cos 0001 (X R ⎥⎥ ⎥⎦⎤

⎢⎢⎢⎣⎡-=θθθθθcos 0sin 010sin 0cos (Y R ⎥⎥ ⎥⎦ ⎤

⎢⎢⎢⎣⎡-=1000cos sin 0sin cos (θθθθθZ R 1.2 计算单位和有关常数

轨道计算采用人卫单位系统，具体定义为：单位：地球的总质量e M

时间单位：2 2 1

310790870681112412. 8⨯=⎪⎪⎭⎫ ⎝ ⎛=e e

长度单位：地球赤道半径e a 质量

as GM a

T 在以上人卫单位系统中，引力常数G=1。为完整起见，给出以下常数： 地球赤道半径：6378137m 地球扁率：257. 298=f 地球总质量：kg M e 2410974. 5⨯=

地球自转角速度：s rad e /10292115. 75-⨯=ω 地心引力常数：2314/10986005. 3s m GM e ⨯= 日心引力常数：2320/1032712438. 1s m GM s ⨯= 天文单位长度：m a sun 11104959787. 1⨯= 月球地球质量比：01230002. 0=e L GM 光速：s m c /299792458=

引力常数： /(10 672. 62311

s kg m G ⋅⨯=- 太阳常数： /(105605. 426s m kg P s ⋅⨯=- 地球引力位系数采用WGS-84 EGM的规格化值。参见表2。 表1 黄经和倾角章动序列表

表2 计算地球引力位加速度的引力位系数（GEM-T3）

11 12

2 轨道动力学计算的基本数学模型

卫星在轨道上运行要受到各种力因素的影响，产生的摄动是多方面的。国内外一些学者对卫星轨道的受摄问题作了详细的研究与分析，尤以澳大利亚的C.Rizos 、A.Stolz 和美国的H.F.Fliegel 等人为代表。统筹考虑精度的需要和时间耗费，通过大量试算，本课题考虑了卫星所受的以下作用力来进行轨道计算（注：时间系统采用TDT 时间系统、坐标系统采用J2000.0惯性坐标系），主要包括：地球质心引力F 0、除质心外的地球引力F E 、太阳和月球引力F N 、太阳辐射压力F A 、大气阻力（低卫星轨卫星）、Y 轴偏差F Y 、地球潮汐附加力F T 。

Y T A N E F F F F F F F +++++=0 (2.1 其中地球质心引力F 0是最主要的，其次是地球的非质心引力F E ，称为地球非球形摄动力。如果将地球质心引力视为1，地球非球形摄动力可达10-3

量级，而其它摄动力则大多在10-6 以下。 2.1 二体问题

在惯性系中，卫星运动方程被描述为 , , ( , , (3 0r r t r r r

GM r r t r r r e +-=+= (2.2 其中：r r , 和r

分别表示时刻卫星在惯性系中的位置、速度和加速度矢量； G 和e M 为分别为引力常数和地球总质量。

（3.2）式的第一项为地球质心引力项，称为二体运动，是力模型的主项；第二项为摄

动力引起的总摄动项，是r r t

, , 的函数矢量。 2.2 地球非球形引力摄动

在地固坐标系中，地球引力位函数作为拉普拉斯方程的解，其非球形部分为： ∑∑==+=N n n m m n nm m n nm U 20 ( (2.3 式中： 13 1 1

sin (sincos (sin++= = n m n n e e m n n m n n e e m n

r m a GM r m a GM λφλφ (2.4

其中：λ和φ表示单位质点在地固坐标系中的地心经、纬度； e a 表示地球平均半径；

(sinφm n 为规格化的缔合勒让德多项式； nm 和nm 为规格化的地球引力位系数； n和m 为多项式的阶和次，N 为取的最高阶数。 非球形引力摄动的求解，按如下步骤进行：

a. 首先求解m n 和m n 。用下列公式逐阶次推算得到，即： ⎪⎪⎩⎪⎪⎨⎧⎥⎥⎦⎤⎢⎢⎣⎡⎥⎥⎦⎤ ⎢⎢⎣⎡=⎥⎥⎦⎤⎢⎢⎣⎡⎥⎥⎦

⎤⎢⎢⎣⎡-⎥⎥⎦⎤⎢⎢⎣⎡+=⎥⎥⎦⎤⎢⎢⎣⎡++++++---+++}sin cos cos cos {}sin 12{11111

1111111λφλφφn n n n n n n n n n e n n n n m n m n m

n m n m n m n e m n m n K r a J n I r a (2.5 式中： (, (22y x z

arctg x y arctg +==φλ ⎪ ⎩⎪ ⎨⎧>++==⋅ --+=

+-+++= ++-+0, 2 23 20 , 12 (( 1(1(3 211 11n n n n K r

a n m n m n J m n m n n I n n e m n m n (2.6

递推方法如下所示： 递推初值为： 0, 0, sin , 0100000100==== r a

r GM e e φ (2.7 递推方法采用先求对角线，再按行递推(m不变，n 增加 进行。 14

如果采用非规格化的系数U 和V 进行计算，递推公式为： ⎪⎪⎩⎪⎪⎨⎧⎥⎥⎦⎤⎢⎢⎣⎡⎥⎥⎦⎤⎢⎢⎣⎡+=⎥⎥⎦⎤

⎢

⎢⎣⎡⎥⎥⎦⎤⎢⎢⎣ ⎡+-⎥⎥⎦⎤⎢⎢⎣⎡++-= ⎥⎥⎦⎤⎢⎢⎣⎡++++--++}

sin cos cos cos { 12(} (sin 12{( 1(111 111

11λφλφφn n n

n n n n n e n n n n m n m

n e m n m n e m n m n U V V U r a n V U V U m n r a V U n m n r a V U (2.8 递推初值为：

0, 0, sin , 010000010 0==== V V r a

U r GM U e e φ 递推方法同上。此方法完成递推后，还需将U 和V 规格化，其公式为：

m n m n

n n m n m n n U n ! (! ( 12(212+-+=+= (2.9

b. 计算m n 和m n 的偏导数m n ∇和m n ∇ i. 首先计算 32/( 12(1(1 32/( 12(2(12

32/( 12(2(1(2+++++-=++++++=+++-+-=n n m n m n n n m n m n k n n m n m n k m n b m n a

m n (2.10 其中乘数因子为

15 ⎪ ⎩⎪

⎨⎧==>=1, 20, 2/21

, 1m m m K a ⎩⎨

⎧=>=0, 2/20, 1m m K b (2.11 ii. 计算m n ∇和m n ∇ ⎪⎪⎪⎪ ⎩

⎪⎪⎪⎪⎨⎧⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-+-=∇⎥⎥ ⎥⎦

⎤⎢⎢⎢⎣⎡----=∇+++-+++-++++-+++-+m n m n m n m n m n m n m n m n m n m n e m n

m

n m n m n m n m n m n m n m n m n m n e m n a a 11 1111 111111111

11111 (2.12 而：⎪⎩⎪⎨⎧-=-=+--m n m m n m n

m m n 1 1( 1( c. 求非球形引力摄动引力位

∑∑==∇+∇=∇N n n m m n nm m n nm U 20 ( (2.13

则摄动加速度为： U W r T E ∇=][ (2.14

式中：矩阵[W]为地固坐标系转换至惯性坐标系的旋转矩阵。 2.3 日、月摄动

考虑卫星的N 体影响，只需顾及太阳和月球的引力作用已满足精度要求。N 体摄动模型的建立是基于牛顿第二运动定律和万有引力定律。由图2所示，分析卫星的受力可得卫星

图2 三体的几何关系 的日、月摄动加速度矢量为： 16 ∑=+ ∆ ∆-=L

S j j j j j j N r r GM r , 33 ( (2.15 其中：j j r r

-=∆，j ∆为摄动体至卫星的中心距，即矢量j ∆ 的模； r r j , 分别为摄动体j M 和卫星在惯性系中的位置矢量，r r j , 的计算参见本节附录。

这里S 和L 分别代表太阳和月球。

太阳和月球对卫星的摄动影响主要呈长周期变化，且与卫星轨道对太阳和月球的定向有关。

2.4 太阳直接辐射压摄动

照射在卫星上的太阳光，一部分被其吸收转化为热能，另一部分被反射向太空。因此，卫星会受到照射时的辐射力和反射时的反射力的作用，这里统称为直接辐射压摄动。直接辐射压摄动与光压强度、卫星表面的反射系数和光照面积有关。

由光压和牛顿第二运动定律建立的直接辐射压对卫星产生的运动加速度矢量为：

( 2sin 025. 0cos( 2sin 015. 0sin(232 2ααααγ+-+-⎪⎪⎭

⎫ ⎝⎛∆=z x s e as sun A G G a T a k r (2.16 其中：ααα6sin 100. 44cos 103. 12cos 104. 210013. 110

997----⨯+⨯+⨯-⨯=k ；

z x G G , 对应卫星表面反射系数项，是为补偿模型不足而引入的拟合参数； ⎪⎪⎭

⎫ ⎝⎛∆∆-=s s r r , (arccos α； m a sun 11104959787. 1⨯=为天文单位长度（地球轨道长半径）；

210790870681112412. 8⨯=as T 为人卫时间单位； s ∆为卫星的日心距，即矢量r r s s

-=∆的模，s r r , 分别为卫星和太阳在惯性系中的位 置矢量；

γ为蚀因子，具体定义为： ⎪⎪⎩

⎪⎪⎨⎧-=本影、半影之内卫星在地球或月球的伪影之内卫星在地球或月球的本卫星在地影和月影之外s s A A '

101γ (2.17

其中：s s A A , ' 分别为太阳的被蚀视面积和视面积。要确定蚀因子，需计算下列诸量：

⎪⎩⎪⎨⎧===2 2

2e e

m m s s A A A παπαπα (2.18 17 2

22222/ (cos

/(r y x y x z arctg +=+=φφ ⎪⎪⎪

⎩⎪⎪⎪⎨⎧=∆=∆=--- (sin (sin (sin ' 1' 1' 1r a a a e e m m m s s s ααα (2.19 ⎪⎩⎪⎨⎧-+=-=φ2 ' ' cos 20000(e e e s s s f a a f a a (2.20 ⎪⎪⎩ ⎪⎪⎨

⎧∆⋅∆∆⋅∆=∆⋅∆⋅=--

(cos (cos 11s m s m ms s s es r r

θθ (2.21 式中：e m A A , 为月球和地球的视面积；

e s f f , 分别代表太阳和地球的扁率，310352813178. 3-⨯=e f ； '

e a 为考虑地球大气衰减以及地球扁率效应的有效地球赤道半径； m a =1738000m为月球半径；

m a s 8' 1096. 6⨯=为考虑太阳扁率的太阳有效半径； φ为地心纬度；

es θ是地球—卫星—太阳张角； ms θ是月球—卫星—太阳张角。 地影和月影判断过程如图3所示。 当卫星在地球半影中时：

221212' (cos ( cos e s es e e

es e s e s e A θθαθθααθα---+=-- 其中：es e s es e θααθθ22 22 -+=

当卫星在地球伪本影中时：2'

e e A πα=

当卫星在月球半影中时： 2 21212' (cos ( cos m s ms m

m ms m s m s m A θθαθθααθα---+=-- 其中：ms m

s ms m θααθθ22 22-+=

当卫星在月球伪本影中时：2 ' m m A πα= 18

则蚀因子为：2 ' ' ' , max(11s

m e s s A A A A παγ-=-=

图3 地影和月影判断流程

太阳直接辐射摄动对卫星轨道的影响是十分复杂的，它与卫星表面的反射特性、卫星轨道相对太阳的定向以及太阳活动等有关。卫星是由各种不同折射性质的原材料构成的不规则形体，在其运行过程中，太阳对它照射的面积也在不断地改变（它的太阳能翼始终是面向太

阳的）。此外，由于太阳活动的变化，所谓太阳常数s P 也并非常数。因此，给出卫星受太阳直接辐射压摄动的精确模型是很困难的。所以采用较简单的平面模型计算太阳辐射加速度影响。

19

2.5 地球固体潮摄动

地球并非刚体，它受日月引力作用会产生弹性形变，称为潮汐现象。这种形变使得地球内部物质发生小的变化，随之导致引力位函数产生小的形变位差——潮汐位。卫星运动的地球固体潮摄动就是潮汐位效应的结果。

已知日（或月）对地面点的引力位球谐展开式为： ∑+= (cos1j n n j n

e j j P r a GM v ψ

式中：j M 为日（或月）的质量，j j r ψ, 分别为引力体至地心距和与地面点的地心角， (cos

j n P ψ为n 阶勒让德多项式。 从上式中排除不产生形变位差的0或1阶项，且只取n=2阶项，可得日月潮汐形变对卫星产生的摄动位为：

∑==M S j j e

j j et P k r a r GM U , 2 235

3 (cosψ 其中：2k 为二阶Love 数（取3. 02=k ）。 顾及关系式]1 ˆˆ(3[2 1 (cos22-=

r r P T j j ψ，最后得到卫星的固体潮摄动加速度矢量为： } (6] (153{[2, 24

5

32j T j L S j T j e j j T et T r r r r r r r a r GM k r U r +-=⎪⎭ ⎫ ⎝⎛∂∂=∑= (2.22 式中： , (ˆ, ˆL S j r r j =分别为j r r

, 的单位矢量。 2.6 大气阻力摄动

高层大气对卫星的运行将产生阻力，这种阻力对低轨道卫星是主要摄动力之一，但大气

密度变化机制非常复杂，不但随高度变化，也与太阳活动、时间、季节、纬度和地磁活动的变化有关。本文采用了静止球型大气密度模型（HP 模型）。若只考虑大气分子对卫星表面的法向作用力而忽略其切向作用力，大气阻力使卫星产生的摄动加速度为:

D P D B D r r r +=R R D V V m A C ⎪⎭⎫ ⎝⎛-=ρ21s R P DP V m A C ∆⎪⎭ ⎫ ⎝⎛- 2cos 21ϕρ (32.23 其中：DB r

为卫星星体部分的大气阻力摄动加速度； DP r

为卫星太阳帆板的大气阻力摄动加速度； D C 为大气阻力系数； ρ是大气密度，由HP 模型计算得到； R V 是卫星相对大气的速度矢量r r V R ⨯-=ω；

20 m A

是卫星参考面积与卫星质量之比； DP C 为太阳帆板的大气阻力系数； P A 为太阳帆板的面积；

ϕ为太阳帆板的法向与卫星相对于大气的速度方向的夹角； s ∆ 为 s ∆ 向量的单位向量。 2.7 Y轴偏差加速度摄动

Y 轴偏差加速度主要是GPS 卫星的结构失调和卫星体的热辐射而产生的一个附加加速度在星固坐标系Y 轴方向上的分量。在设计上，为使卫星的太阳翼以最大面积朝向太阳，两翼的支轴应保持在一条直线上，并要求垂直于卫星和太阳方向的连线，用于控制两翼俯仰的太阳传感器也应完全垂直于卫星翼支轴，卫星的偏航高度控制应保持正确，但事实上并不完全这样；另一方面，由卫星本身产生的超高温要从Y 轴方向的百叶孔排出，这样使处在不稳定状态中的卫星体也会产生结构失调现象。由此导致了Y 轴偏差加速度的存在。H.F.Fliegel 等人把结构失调的影响表示为

[14] ：

a y Y y G r = (2.24

其中：y G 为考虑模型剩差而引入的待估参数； a y

为星固坐标系的Y 轴在J2000.0惯性坐标系中的方向余弦，由下式决定： r r s s -=∆ r r y s s a ⨯∆⨯∆=

2.8 巡航姿态控制动力摄动

有些卫星在巡航过程中需要保持三轴稳定姿态，需通过姿态控制实现，有的卫星其姿态控制的动力来源于高压气瓶的喷气。这样，在姿态控制的同时也影响了卫星质心的运动。该摄动加速度矢量可以表示为:

[](

u S u C T C C W r o T A sin cos 1 +++= (2.25 其中：f u +=ω； o C

为卫星在RTN 坐标系中姿控动力摄动加速度的常数分量； 1C

为卫星在RTN 坐标系中姿控动力摄动加速度的时间变化率； C 和S

为周期项的系数；

对全局参数，为由历元时刻起算的相对时间；对弧段相关参数，为观测时刻t 在相应弧段内的相对时间。

2.9 其它摄动影响

在轨道运行的卫星除受到上述摄动作用外，还受其它一些摄动的影响，如反照辐射摄动、地球自转形变摄动、广义相对论摄动、海潮摄动、大气潮摄动等。这些摄动影响对卫星轨道摄动非常小，但其计算却要耗费大量的时间片。考虑到2.1节所述的摄动已能满足课题的精度需求和时间消耗的，因此，本课题中忽略了其它摄动的影响。

附录：日月位置计算 a. 太阳位置矢量s 计算 ⎥⎥ ⎥⎦ ⎤

⎢⎢⎢⎣⎡-=0sin cos (][0s s s s x T s r r R D r θθε 其中：[D]T

为计算历元时刻的平春分点到J2000.0平春分点的转换矩阵，即岁差矩阵。 2

22210126. 010418. 001675104. 02sin 4 5 sin 2] 2cos 2

1

cos 211[B B s sun s T T e M e M e L M e M e e a r --⨯-⨯-=++=--+=θ 计算当时平春分点的几何平黄经为： 2' ' '

' ' ' ' 0. 113. 12960276804. 4841279B B T T L ++= 平近点角为：

2' ' ' ' ' ' ' . 010. 12959657904. 3328358B B T T M -+= 平黄赤交角为：

2' ' ' ' ' ' ' 00059. 0845. 4626. 082723B B T T --= ε 以上式中 . 365252415020 (-= TDB JD T B

b. 月球位置矢量m 计算 ⎥⎥⎥⎦ ⎤

⎢⎢⎢⎣⎡Ω-Ω--Ω--=0 ~~sin( ~

~cos( ~( ~( (][0λλεm m x z x T m r r J R R R D r 其中：∑=-Ω+⨯+Ω=Ωn

j j j a M 1 '

3sin sin 1047. 0~ ∑=+=n j j j a J J J 1 cos ~ ∑=+ +=n j j j a K C L 1 1' sin ~ λ

]cos 1[26. 601112∑=++=n j j j a P C r e m a r r ⋅= 14. 5 =J

208\". 431\". 1732537999\". 02' 26270' B B T T L -+=

2' ' ' ' ' ' ' 48. 723. 696291179. 5910259B B T T +-=Ω 2' ' ' ' ' ' ' ' 0912. 3379. 171791585659. 1606296B B T T M ++=

' 3' ' 13sin 1017. 02sin 00373. 0sin 10976. 0M M M C -⨯++= ' 3' ' 23cos 1018. 02cos 00297. 0cos 05. 0M M M C -⨯++=

系数j K 、j P 、j Ω和j J 的值以及幅角j a 的计算列入表5中，表内D 为日月平角距，

Ω-=' L F ，其计算式为：

2' ' ' ' ' ' ' 17. 518. 16029161195. 1444350B B T T D -+= 2' ' ' ' ' ' ' 56. 11. 173952729020. 031511B B T T F -+=

表5 计算月球位置的系数表

3 轨道计算方法

卫星运动方程的解有分析法、数值法和半分析半数值方法等。分析法是将力模型展开取有限项，给出任意时刻解的表达式，通常称之为一般摄动法。其优点是便于定性地给出轨道的特征，如长期项、长周期项影响等。但对摄动力模型处理较大，使精度受到影响。数值法即常微分方程的初值解问题，在轨道计算中称之为特别摄动法。其优点是能完整地顾及所选择的力模型、处理简单、计算精度高，是高精度卫星轨道计算最实际有效的方法。

经研究，Runge_Kutta型和Cowell 型数值积分方法都可用于产生高精度的卫星轨道，尤其是后者通过预报——校准公式进行迭代计算，收敛快（一般只需迭代2~3次），累积误差影响要比前者小的多。所以这一方法是轨道计算最常用的方法之一。不足的是它为多步型积分器，必须用其它方法为它准备起步值。

本课题提供的积分模式有Runge_Kutta 4、Runge_Kutta 8和高精度的Adams_Cowell方法，高精度计算时采用Runge_Kutta单步积分起步，用Adams_Cowell多步积分求解卫星运动方程，确定卫星在J2000.0惯性坐标系中的位置和速度矢量。

3.1 Runge_Kutta积分法

Runge_Kutta法是一种单步积分方法，公式形式简单，应用方便，将它用于解算卫星运动方程和变分方程也具有很好的稳定性。但该方法随着积分式阶数的增高，计算右函数的次数也会随之增多，无疑会导致累积误差增加，同时也会降低计算效率。因此，在大规模的轨道计算中，通常只将其作为多步型积分器的起步方法。流程图如图4所示。

设有初值问题： ⎩⎨ ⎧==n n y t y y t f t y

( , ( ( (3.1 则求解该问题的Runge_Kutta 积分公式为： ∑=++=k

i i i n n f C h y y 0 1 (3.2

图4 Runge_Kutta 积分框图 ⎪⎪⎪⎩ ⎪

⎪⎪⎨⎧++=++==∑= , ( , ( , (00101110k

i i ki k n k n k n n n n f b h a y h a t f f hf b a y h a t f f y t f f 式中：n n t t h -=+1，为积分步长； k为积分式所取的阶数；

i C , ij i b a , 均为已知的系数，具体参见表6所示。

由运动方程和变分方程形成的矩阵常微分方程的初值问题为： , ( (y t F t y

= 上式中：左端⎪⎪⎭⎫ ⎝⎛=ψψ r r t y (； 右函数⎪⎪⎭ ⎫ ⎝⎛= , , ( , , ( ( (r r t r r t r t t r F ψψ ； 初值⎪⎪⎭⎫ ⎝⎛=00 000 (ψψ

r r t y 3.2 Adams_Cowell积分

Adams_Cowell积分适用显含1阶导的二阶微分方程，初始值需要位置项和速度项。积分的预报公式为：

⎪⎪⎩ ⎪ ⎪⎨⎧

+=+=∑∑=-+=-+][][02101k

j j n j sn n k j j n j sn n y II h y y I h y

βα (3.3 其中：j j βα, 为Adams_Cowell积分公式的预报系数，按下列公式计算，结果参见表4.2。

⎪⎪⎩

⎪⎪⎨⎧''⎪⎪⎭⎫ ⎝⎛-='⎪⎪⎭ ⎫ ⎝⎛-=∑∑=+=+k j m m j j k j m m j

j j m j m 1( 1( 1( 1(21γβγα 校正公式为： ⎪⎪⎩ ⎪

⎪⎨⎧+=+=∑∑=+-+=+-+][][01' 2101' 1 k

j j n j sn n k

j j n j sn n y II h y y I h y βα (3.4 式中：' ' , j j βα为Adams_Cowell积分公式的校正系数，按下列公式计算，结果参见表4.2。

⎪⎪⎪⎪ ⎩ ⎪⎪

⎪⎪⎨⎧''⎪⎪⎭

⎫ ⎝⎛-='⎪⎪⎭⎫ ⎝⎛-=⎪⎩⎪⎨ ⎧>=+=∑∑=+=+k j m m j j k j m m j

j j j j j m j m j j 1( 1( 1( 1(0, 0, 12' 1\" \" \"

' γβγαααα (3.3、(3.4式中的 1(γ'和 1(γ''按以下公式递推计算： ⎪⎪⎪⎪ ⎪⎪

⎪ ⎩ ⎪ ⎪⎪ ⎪⎪⎪⎪

⎨⎧''='''='''='''+--='='∑∑∑∑===--=i j j i i j j i i j j i j i i j j i j i 000100 0( 1(

0( 1( 0( 0( 0( 0(11 0(1 0(γγγγγγγγγγ sn I 和sn II 分别为n y

的一次和分和二次和分，按下列公式递推计算：⎧+=+=--sn

sn sn n sn sn I II II y I I 11 sn I 和sn II 的初值定义为：⎪⎪⎩ ⎪

⎪⎨⎧-=-=∑∑=--=--k j j n j n sn k

⎩⎨

j j n j n sn y h y II y h y I 0' 210

' 1 βα 用Adams_Cowell积分解卫星运动方程积分流程图如图5所示： 3.3 轨道计算

预测卫星轨道实际上就是由一组正确的初始值数值积分卫星运动方程，外推出未来的卫星轨道。显然，在卫星运动力模型较完备的情况下，初始值的精确与否将直接影响外推轨道的精度。且在推算过程中，由于受初始值误差传播影响，推算区间越长，积分累积的误差就越大。轨道计算步骤为：

图5 Adams_Cowell积分流程图

图6 卫星轨道计算流程图 3.4 星历的快速插值

通过分布式计算提供的星历数据间隔太大，而在实际定位计算中，有时需要间隔为1秒甚至步长更小的精确星历坐标。因此，高精度、快速的精密星历插值就成为一项重要工作。以精度和计算耗费小为出发点，本课题在研究的过程中，分析和比较了以下几种多项式插值方法：拉格朗日多项式插值、牛顿多项式插值、三角函数插值、样条函数插值、切比雪夫多项式插值等。

拉格朗日多项式插值的多项式构造和插值时间耗费非常大。其时间耗费可由下式表示：

D n M n A n 1( 2( 22(-+-+- 其中：n 表示已知点数；

A表示加法运算、M 和D 分别表示乘、除运算。

构造插值多项式还需额外的nM A n +- 1(运算。因此，总运算耗费为： D n M n A n 1( 22( 33(-+-+-

比拉格朗日多项式插值更有效的方法是利用数据均差的牛顿多项式插值[1][2] 。假设有

1+n 个控制点n t t t , , , 10 ，定义在i t 处的零阶均差为 (][i i t f t f =，在j i t t , 处的1阶均差定义为[]

[][]j i j

i j i t t t f t f t t f a --= =, 0，k 阶均差定义为： [][][]0, , , , , , , , , 02111010>--= =-k t t t t t f t t t f t t t f a k k k k k

则牛顿插值多项式可表示为： 1

11221100}}} ({({({( (-----+-++-+-+=n n n n n n a t t a t t a t t a t t a t p (4.12 上述表达式在每级中都由2个加法和1个乘法运算组成，其时间耗费为： 2(M A n +。与经典的拉格朗日多项式插值相比，牛顿多项式插值经济得多。给定1+n 个不同的点

n t t t , , , 10 和相应的均差系数n a a a , , , 10 ，对任意时刻[]n t t t , 0∈插值多项式的结果由下列迭代产生（Horner 算法）。

end

b t t a b by step to n k for a b k k k k n n 1 (1 01+-+=--==

上述算法的精度和拉格朗日的精度相同，但时间耗费得到了明显的改善，而且各阶均差可以预先运算好储存起来，需要时直接调用即可。

插值点等间隔条件下的插值是牛顿多项式插值的一种特例，可以导出更简洁的形式，而不降低精度。IGS 提供的精密星历数据是等间隔的，步长为15分钟。在实际应用中，我们采用牛顿前向差分形式。

牛顿前向均差的定义为：

0, , ( ( ( ( (01>+=-+=-=∆+k kh t t t f h t f t f t f t f k i i i i i 同理可推出各高阶前向均差为： ( 1( ( ( (( (01 11 1 j i k

j j k i k i k i k i k t f j k t f t f t f t f +=--+--∑⎪⎪⎭ ⎫

⎝⎛-=∆ -∆ =∆ ∆=∆ (4.13

假设插值多项式的阶数为n ，则插值多项式为： (! ( (( ( ( ( (0110000t f h n t t t t t t t f h t t t f t p n n

n n ∆---++∆-+

=- (4.14 令h t t (0-=τ，则上式可表示为： ( (00t f j p j n j n ∆⎪⎪⎭ ⎫ ⎝⎛=∑=ττ，

用(4.13式构造插值多项式，再用(4.14式用作插值，该方法在时间耗费上比Horner 算法有一定的改善。

与经典的拉格朗日插值方法相比，三次样条多项式插值速度快，但精度低，三角函数多项式插值必须用DFT 求解系数，时间耗费大，且这两种方法都将星历数据近似为周期函数，插值精度低，不能满足高精度插值的需要。切比雪夫多项式插值方法在端点处抑制了误差的扩大，但计算量较大。牛顿多项式插值方法的精度高，速度快，充分利用了星历数据的特点，且克服了拉格朗日方法的不足。在精密星历插值上，牛顿多项式插值方法比拉格朗日多项式插值方法更为有效。由拉格朗日插值多项式和牛顿插值多项式的余项公式可推出截断误差在

端点处将迅速增长。为克服这种不足、提高插值的精度，在实现时提出并采用了“滑动式牛顿多项式插值”方法

[19]

，即采用滑动方式进行插值：每次取13个点，生成12阶多项式，始 终让被插值点在已知点的第6个点到第7个点之间。第1个点到第13个点为第一个插值区间，仅用来插值第6个点到第7个点之间的时间段；第2个点到第14个点为第二个插值区间，仅用来插值第7个点到第8个点之间的时间段；依次类推下去。这种滑动式牛顿多项式插值精度非常高，可达到cm 级精度，且易于实现，时间耗费小，确保了精密轨道解算时的精度和实时要求。

表6 8阶Runge_Kutta积分公式系数表

表7 12阶Adams_Cowell积分公式系数表

2 3 4 5 6 7 8 9 10 11 12 80.019378 -196.294318 349.412940 -462.480193 460.087276 -343.735149 190.394308 -75.976728 20.680772 -3.441245 0.2351 16.861107 -41.187179 73.0869 -96.516821 95.852539 -71.516859 39.571065 -15.7771 4.291462 -0.713663 0.0794 -2.071621 4.4143 -7.283104 9.1927 -8.853279 6.460362 -3.5149 1.383094 -0.372243 0.061367 -0.004677 -0.526813 1.1914 -2.009198 2.566623 -2.4666 1.825385 -0.9993 0.393084 -0.105997 0.017501 -0.001336 31

4 轨道根数与位置矢量、速度矢量的关系 轨道根数与位置矢量、 4.1 由位置矢量和速度矢量计算轨道根数 & 已知某时刻 t 的位置矢量 r 和速度矢量 r ，计算轨道根数的步骤如下： 由 r r 其中 r 2 2 2 h = h = h X + hY + hZ ，分别计算 i 和

（5.1） （5.2） 而

（5.3） 计算 ω tan ω = 计算 a eZ (eY sin

Ω + e X cos Ω sin i （5.4） a= 类似可得 h2 µ (1 − e 2 （5.5） tan u = 所以： Z (Y sin Ω + X cos Ω sin i f = u −ω (5.6 (5.7 计算 τ ，由 tan 得 E 1− e f = tan 2 1+ e 2 32

n(t − τ = E − e sin E (5.8 4.2 由轨道根数计算位置矢量和速度矢量 若已知任何时刻 t 时天体的 r 和 f ，在轨道坐标系 O − x ′′y ′′z ′′ 中有

据惯性坐标系与个坐标系的转换关系可得位置矢量的表达式

（5.9） 根 将旋转矩

阵展开可得 r r 其中 P 和 Q 分别表示 x ′′ 轴和 y ′′ 的单位矢量，即 r r r r = r cos f ⋅ P + r sin f ⋅

速度矢量表达式为： 2 2 r r r & = − a n sin E ⋅ P

+ a n 1 − e 2 cos E ⋅ Q r r r (5.13 33