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

《卫星轨道动力学数值计算.doc》由会员分享，可在线阅读，更多相关《卫星轨道动力学数值计算.doc（33页珍藏版）》请在三一办公上搜索。

1、目 录1星历计算的时间和坐标系统21.1 有关的时间系统与坐标系统21.1.1 时间系统及其换算21.1.2 坐标系统及其换算41.2 计算单位和有关常数72 轨道动力学计算的基本数学模型122.1 二体问题122.2 地球非球形引力摄动122.3 日、月摄动152.4 太阳直接辐射压摄动162.5 地球固体潮摄动192.6 大气阻力摄动192.7 Y轴偏差加速度摄动202.8 巡航姿态控制动力摄动202.9 其它摄动影响21附录：日月位置计算213 轨道计算方法243.1 Runge_Kutta积分法243.2 Adams_Cowell积分253.3 轨道计算273.4 星历的快速插值284

2、 轨道根数与位置矢量、速度矢量的关系324.1 由位置矢量和速度矢量计算轨道根数324.2 由轨道根数计算位置矢量和速度矢量331星历计算的时间和坐标系统1.1 有关的时间系统与坐标系统轨道计算过程重要涉及到不同的时间系统和坐标系统，下面将空间战场环境系统中所涉及到的时间系统和坐标系统进行定义，并说明各系统之间的相互关系。一般情况下，仿真系统采用的是TDT时间系统和J2000地心惯性坐标系。1.1.1 时间系统及其换算在轨道计算中，时间是独立变量。但是，在计算不同的物理量时，却使用不同的时间系统。例如：在计算恒星时用世界时UT1；定位解算时采用GPS时GPST；岁差和章动量的计算采用TDB时等

3、。所以必须清楚各时间系统的定义和各时间系统之间的转换，下面给出各种时间系统的定义及它们之间的转换公式。格林尼治恒星时格林尼治恒星时为春分点对格林尼治平天文子午面的时角。由于岁差、章动原因，它由格林尼治真恒星时（GAST）和平恒星时（GMST）之分。两者的关系是：其中：为赤经章动为自起算至观测时刻的儒略世纪数，即世界时是以平北极（国际习惯用原点）为统一标准的观测世界时，是反映地球实际自转的时间，恒星时计算与此有关。国际原子时时以铯原子基态两能级间跃迁辐射的9192631770周所经历的时间作为1秒长的均匀时间，起点在1958年1月1日。国际协调时是经跳秒修改后的国际原子时，它与世界时的差，观测纪

4、录时间是以此为准的。质心动力学时（Barycentric Dynamical Time）为相对于太阳质心的运动方程给出的历表、引数等所用的时间尺度，岁差及章动量的计算是以此为依据的。地球动力学时（Terrestrial Dynamical Time）为视地心历表所用的时间尺度，它具有均匀连续的特性，卫星运动方程就是以此为独立的时间变量。时间是由系统定义和应用的一种时间尺度，于1980年1月6日GPST与UTC相等，在此以后由系统主控站密切跟踪UTC以保持高度统一。但GPST不作跳秒修正，因此它与UTC具有整秒的差异（1997年1月至6月相差为）。在计算GPS卫星轨道的初值时将涉及到GPST，G

5、PS精密星历的参考时间为GPST。以上各时间尺度的相互关系如下：其中：可从地球自转参数文件中获得； ； 。 上式中的T为自年起算至观测TDB时刻的儒略世纪数，即：不同时间系统间的关系如下图所示：图1 几种时间系统之间的关系1.1.2 坐标系统及其换算卫星轨道计算和实际定位解算分别是在J2000.0惯性坐标系与WGS-84地固系中进行的，此外，卫星加速度计算中还涉及到星固坐标系。下面给出与本课题有关的主要坐标系的定义及相互转换关系。地心惯性系原点：地球质心Z轴：向北指向年平赤道面（基面）的极点X轴：指向平春分点Y轴：符合右手系法则位置矢量：星固坐标系原点：卫星质心Z轴：指向卫星的天线方向，即指向

6、地心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=3986005108 m3/s2正常化二

7、阶带球谐系数：=-484.1668510-6地球自转角速度：w=729211510-11 rad/s2地球椭球扁率：f = 1/298.257223563地固坐标系与惯性坐标系的转换模型a. 惯性坐标系地固坐标系 模型b. 地固坐标系惯性坐标系 模型式中：A为极移矩阵； B为自转矩阵； C为章动矩阵； D为岁差矩阵。上述各矩阵的意义及具体定义如下：极移：由于地球不是刚体及其它一些地球物理因素的影响，地球自转轴相对于地球的位置随时间而变化从而引起观察者的天顶在天球上的位置发生变化，称为极移，矩阵为A：式中：为地极坐标，可从地球自转参数文件中给出的极移值插值得到。自转：即地球公转的同时也在绕自转轴

8、旋转。矩阵B：式中：为格林尼治恒星时 章动：是指外力作用下，地球自转轴在空间运动的短周期摆动部分，即同一瞬间真天极相对平天极的运动，月球对地球引力的变化是形成章动现象的主要外力作用，其次是太阳。矩阵C：式中：交角章动黄经章动其中：都为常数，可自章动系数表1中查出。而月亮平近点角：太阳平近点角：月亮平升交角：日月平角矩：月亮轨道对黄道平均升交点的黄经：其中： T0为自J2000.0起算至t的儒略世纪数 岁差：地球在太阳、月球和行星的引力作用下，地球的自转轴在空间不断发生变化，其长期运动称为岁差，矩阵D：式中：T0的意义同上。旋转矩阵的求法分别为1.2 计算单位和有关常数轨道计算采用人卫单位系统，

9、具体定义为：长度单位：地球赤道半径质量单位：地球的总质量时间单位：在以上人卫单位系统中，引力常数G=1。为完整起见，给出以下常数：地球赤道半径：6378137m地球扁率：地球总质量：地球自转角速度：地心引力常数：日心引力常数：天文单位长度：月球地球质量比：光速：引力常数：太阳常数：地球引力位系数采用WGS-84 EGM的规格化值。参见表2。表1 黄经和倾角章动序列表引 数黄经章动倾角章动iLLFD0.00010.0001T0.00010.0001T100001-171996-174.2920258.920000220620.2-8950.53-20201460-240420-200110005

10、-20202-301061-10-10-300070-22-21-2010820-20110009002-22-13187-1.65736-3.110010001426-3.454-0.111012-22-5171.2224-0.6120-12-22217-0.5-950.313002-211290.1-70014200-204801015002-20-22000160200017-0.1001701001-1509018022-22-160.170190-1001-1206020-20021-6030210-12-21-503022200-2140-2023012-2140-2024100-

11、10-400025210-2010002600-22110002701-220-10002801002100029-10011100030012-20-10003100202-2274-0.2977-0.532100007120.1-703300201-386-0.420003410202-3010129-0.135100-20-1580-1036-102021230-5303700020630-203810001630.1-33039-10001-58-0.132040-10222-5902604110201-5102704200222-3801604320000290-1044102-22

12、290-1204520202-3101304600200260-1047-10201210-10048-10021160-8049100-21-1307050-10221-1005051110-20-7000520120270-30530-1202-70305410222-80305510020600056202-2260-305700021-60305800221-703059102-2160-3060000-21-5030611-100050006220201-503063010-20-40006410-20040006500010-40006611000-3000671020030006

13、81-1202-301069-1-1222-301070-20001-20107130202-3010720-1222-3010731120220-1074-102-21-2010752000120-107610002-201077300002000780021220-1079-1000210-1080100-40-100081-2022210-1082-10242-201083200-40-100084112-2210-108510221-101086-20242-101087-104021000881-10-20100089202-2110-109020222-10009110021-10

14、0092004-22100093302-22100094102-20-10009501201100096-1-102110009700-201-100098002-12-10009901020-100010010-2-20-10001010-1201-1000102110-21-100010310-220-100010420020100010500242-1000106010101000表2 计算地球引力位加速度的引力位系数（GEM-T3）20-0.48416510e-30.021-0.17e-91.19e-9220.24390658e-5-0.14000946e-5300.95720109e

15、-60.0310.20277142e-50.24921712e-6320.90447073e-6-0.61944767e-6330.72034249e-60.14138845e-5400.53952118e-60.041-0.53615108e-6-0.47343598e-6420.35021814e-60.66301523e-6430.99093372e-6-0.20092742e-644-0.18877065e-60.30942370e-6500.68343345e-70.051-0.58280231e-7-0.96083937e-7520.65271099e-6-0.32386369e-

16、653-0.45233006e-6-0.21529578e-654-0.29558409e-60.49690346e-7550.17376349e-6-0.66890704e-660-0.14951352e-60.061-0.76894161e-70.26998384e-7620.48734482e-7-0.37401311e-6630.57203194e-70.93727543e-864-0.86826474e-7-0.47130637e-665-0.26733038e-6-0.53678021e-6660.96846267e-7-0.23713482e-6700.91300885e-70.

17、0710.27486873e-60.97465920e-7720.32779498e-60.93246712e-7730.25122012e-6-0.21529269e-674-0.27556101e-6-0.12376718e-6750.13261698e-80.18620005e-776-0.35883140e-60.15173875e-6770.97027653e-90.24083642e-7800.48883181e-70.0810.23628239e-70.58847236e-7820.77598534e-70.66008696e-783-0.17785247e-7-0.863469

18、83e-784-0.24633984e-60.70179584e-785-0.25041147e-70.89426831e-786-0.64923712e-70.30912257e-6870.67462214e-70.75094831e-788-0.12419836e-60.12017218e-62 轨道动力学计算的基本数学模型卫星在轨道上运行要受到各种力因素的影响，产生的摄动是多方面的。国内外一些学者对卫星轨道的受摄问题作了详细的研究与分析，尤以澳大利亚的C.Rizos、A.Stolz和美国的H.F.Fliegel等人为代表。统筹考虑精度的需要和时间耗费，通过大量试算，本课题考虑了卫星所受的

19、以下作用力来进行轨道计算（注：时间系统采用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）

20、式的第一项为地球质心引力项，称为二体运动，是力模型的主项；第二项为摄动力引起的总摄动项，是的函数矢量。2.2 地球非球形引力摄动在地固坐标系中，地球引力位函数作为拉普拉斯方程的解，其非球形部分为： (2.3)式中： (2.4)其中：和表示单位质点在地固坐标系中的地心经、纬度； 表示地球平均半径； 为规格化的缔合勒让德多项式； 和为规格化的地球引力位系数； n和m为多项式的阶和次，N为取的最高阶数。非球形引力摄动的求解，按如下步骤进行：a. 首先求解和。用下列公式逐阶次推算得到，即： (2.5)式中： (2.6)递推方法如下所示：递推初值为： (2.7)递推方法采用先求对角线，再按行递推(m不变

21、，n增加)进行。L0 L0 L0 L0 LO M如果采用非规格化的系数U和V进行计算，递推公式为： (2.8)递推初值为：递推方法同上。此方法完成递推后，还需将U和V规格化，其公式为： (2.9)b. 计算和的偏导数和 i. 首先计算 (2.10)其中乘数因子为 (2.11) ii. 计算和 (2.12)而：c. 求非球形引力摄动引力位 (2.13)则摄动加速度为： (2.14)式中：矩阵W为地固坐标系转换至惯性坐标系的旋转矩阵。2.3 日、月摄动考虑卫星的N体影响，只需顾及太阳和月球的引力作用已满足精度要求。N体摄动模型的建立是基于牛顿第二运动定律和万有引力定律。由图2所示，分析卫星的受力可

22、得卫星图2 三体的几何关系的日、月摄动加速度矢量为： (2.15)其中：，为摄动体至卫星的中心距，即矢量的模； 分别为摄动体和卫星在惯性系中的位置矢量，的计算参见本节附录。 这里S和L分别代表太阳和月球。太阳和月球对卫星的摄动影响主要呈长周期变化，且与卫星轨道对太阳和月球的定向有关。2.4 太阳直接辐射压摄动照射在卫星上的太阳光，一部分被其吸收转化为热能，另一部分被反射向太空。因此，卫星会受到照射时的辐射力和反射时的反射力的作用，这里统称为直接辐射压摄动。直接辐射压摄动与光压强度、卫星表面的反射系数和光照面积有关。由光压和牛顿第二运动定律建立的直接辐射压对卫星产生的运动加速度矢量为： (2.1

23、6)其中：； 对应卫星表面反射系数项，是为补偿模型不足而引入的拟合参数； ； 为天文单位长度（地球轨道长半径）； 为人卫时间单位；为卫星的日心距，即矢量的模，分别为卫星和太阳在惯性系中的位置矢量；为蚀因子，具体定义为： (2.17)其中：分别为太阳的被蚀视面积和视面积。要确定蚀因子，需计算下列诸量： (2.18) (2.19) (2.20) (2.21)式中：为月球和地球的视面积； 分别代表太阳和地球的扁率，； 为考虑地球大气衰减以及地球扁率效应的有效地球赤道半径； =1738000m为月球半径； 为考虑太阳扁率的太阳有效半径； 为地心纬度； 是地球卫星太阳张角； 是月球卫星太阳张角。地影和月

24、影判断过程如图3所示。当卫星在地球半影中时：其中：当卫星在地球伪本影中时：当卫星在月球半影中时：其中：当卫星在月球伪本影中时：则蚀因子为：图3 地影和月影判断流程太阳直接辐射摄动对卫星轨道的影响是十分复杂的，它与卫星表面的反射特性、卫星轨道相对太阳的定向以及太阳活动等有关。卫星是由各种不同折射性质的原材料构成的不规则形体，在其运行过程中，太阳对它照射的面积也在不断地改变（它的太阳能翼始终是面向太阳的）。此外，由于太阳活动的变化，所谓太阳常数也并非常数。因此，给出卫星受太阳直接辐射压摄动的精确模型是很困难的。所以采用较简单的平面模型计算太阳辐射加速度影响。2.5 地球固体潮摄动地球并非刚体，它受

25、日月引力作用会产生弹性形变，称为潮汐现象。这种形变使得地球内部物质发生小的变化，随之导致引力位函数产生小的形变位差潮汐位。卫星运动的地球固体潮摄动就是潮汐位效应的结果。已知日（或月）对地面点的引力位球谐展开式为：式中：为日（或月）的质量，分别为引力体至地心距和与地面点的地心角， 为n阶勒让德多项式。从上式中排除不产生形变位差的0或1阶项，且只取n=2阶项，可得日月潮汐形变对卫星产生的摄动位为：其中：为二阶Love数（取）。顾及关系式，最后得到卫星的固体潮摄动加速度矢量为： (2.22)式中：分别为的单位矢量。2.6 大气阻力摄动高层大气对卫星的运行将产生阻力，这种阻力对低轨道卫星是主要摄动力之

26、一，但大气密度变化机制非常复杂，不但随高度变化，也与太阳活动、时间、季节、纬度和地磁活动的变化有关。本文采用了静止球型大气密度模型（HP模型）。若只考虑大气分子对卫星表面的法向作用力而忽略其切向作用力，大气阻力使卫星产生的摄动加速度为: (32.23)其中：为卫星星体部分的大气阻力摄动加速度； 为卫星太阳帆板的大气阻力摄动加速度；为大气阻力系数；是大气密度，由HP模型计算得到；是卫星相对大气的速度矢量；是卫星参考面积与卫星质量之比； 为太阳帆板的大气阻力系数；为太阳帆板的面积；为太阳帆板的法向与卫星相对于大气的速度方向的夹角；为 向量的单位向量。2.7 Y轴偏差加速度摄动Y轴偏差加速度主要是G

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

28、方向余弦，由下式决定：2.8 巡航姿态控制动力摄动有些卫星在巡航过程中需要保持三轴稳定姿态，需通过姿态控制实现，有的卫星其姿态控制的动力来源于高压气瓶的喷气。这样，在姿态控制的同时也影响了卫星质心的运动。该摄动加速度矢量可以表示为: (2.25)其中：； 为卫星在RTN坐标系中姿控动力摄动加速度的常数分量； 为卫星在RTN坐标系中姿控动力摄动加速度的时间变化率； 和为周期项的系数；对全局参数，为由历元时刻起算的相对时间；对弧段相关参数，为观测时刻t在相应弧段内的相对时间。2.9 其它摄动影响在轨道运行的卫星除受到上述摄动作用外，还受其它一些摄动的影响，如反照辐射摄动、地球自转形变摄动、广义相对

29、论摄动、海潮摄动、大气潮摄动等。这些摄动影响对卫星轨道摄动非常小，但其计算却要耗费大量的时间片。考虑到2.1节所述的摄动已能满足课题的精度需求和时间消耗的限制，因此，本课题中忽略了其它摄动的影响。附录：日月位置计算a. 太阳位置矢量计算其中：DT为计算历元时刻的平春分点到J2000.0平春分点的转换矩阵，即岁差矩阵。计算当时平春分点的几何平黄经为：平近点角为：平黄赤交角为：以上式中b. 月球位置矢量计算其中： 系数、和的值以及幅角的计算列入表5中，表内D为日月平角距，其计算式为： 表5 计算月球位置的系数表j10.0222300.0100200.0008200.0000002D-M20.011

30、6100.008250-0.002300-0.0001902D30.0032400.0001200.0028300.000000M+18040.0000000.000000-0.002130-0.0001902F+18050.0010300.0000900.0004500.0000002M-2D+18060.0010000.0004200.0000000.0000002D-M-M70.0009300.000900-0.0000800.0000002D+M80.0008000.000560-0.0001200.0000002D-M90.0007200.0003400.0000000.000000

31、M-M100.0006100.0002900.0000000.000000D+PI110.0005300.0002800.0000000.000000M+M+180120.0002700.000000-0.027240-0.0024102F-2D+180130.0000000.000000-0.0000800.0000002F+M+180140.0004100.0002000.0007000.0000002F-M+180150.0001800.0001800.0000000.0000004D-M160.0001500.0001200.0000000.0000004D-2M170.0001400

32、.0000700.0000000.0000002D+M-M+180180.0001200.0000900.0000000.0000002D+M+180190.0000900.0000000.0000000.000000M-D200.0000900.0000000.0000000.000000D+M210.0000700.0000700.0000000.0000002D-M+M220.0000700.0000900.0000000.0000002D+2M230.0000700.0000800.0000000.0000004D240.0000000.0000000.0001900.0001802F

33、-2M250.0000000.0000000.0011000.0001002F-2D+M260.0000000.0000000.0004700.0000002D-2F+M270.0000000.0000000.0002300.0000002F-2D-M280.0000000.0000000.0000230.0000002D-2F-M3 轨道计算方法卫星运动方程的解有分析法、数值法和半分析半数值方法等。分析法是将力模型展开取有限项，给出任意时刻解的表达式，通常称之为一般摄动法。其优点是便于定性地给出轨道的特征，如长期项、长周期项影响等。但对摄动力模型处理限制较大，使精度受到影响。数值法即常微分方

34、程的初值解问题，在轨道计算中称之为特别摄动法。其优点是能完整地顾及所选择的力模型、处理简单、计算精度高，是高精度卫星轨道计算最实际有效的方法。经研究，Runge_Kutta型和Cowell型数值积分方法都可用于产生高精度的卫星轨道，尤其是后者通过预报校准公式进行迭代计算，收敛快（一般只需迭代23次），累积误差影响要比前者小的多。所以这一方法是轨道计算最常用的方法之一。不足的是它为多步型积分器，必须用其它方法为它准备起步值。本课题提供的积分模式有Runge_Kutta 4、Runge_Kutta 8和高精度的Adams_Cowell方法，高精度计算时采用Runge_Kutta单步积分起步，用Ad

35、ams_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。校正公式为：