第七章 航天器轨道摄动
在二体问题的研究中,假设天体的引力场是与距离的平方成反比的中心引力场,航天器在其
作用下沿圆锥曲线运动,真实的情况与此是有差别的。天体的形状往往是不规则的球形,质量的
分布也是不均匀的,这将导致真实的引力场并非理想的中心引力场;大气阻力、其它天体的引力、
太阳辐射压力等也不能够完全忽略,这些因素都会导致航天器的实际运动偏离圆锥曲线。这种运
动偏离二体轨道的现象称为轨道摄动,相应的轨道称为摄动轨道或受摄开普勒轨道。
轨道摄动理论是经典天体力学的主要研究内容之一。在 17~19 世纪,为满足航海定位对精密
星表的需求,学者们对大行星和月球的轨道摄动开展了深入研究,主要工作是建立行星运动方程
近似解的分析理论。牛顿在《原理》中揭示了太阳摄动是导致月球运动周期差和近地点进动的主
要原因。欧拉首创了轨道要素变分法,开创了摄动理论的分析方法,分析了木星、土星和月球的
轨道摄动问题。拉格朗日是大行星运动理论的创始人,提出了著名的拉格朗日行星运动方程。克
莱洛、拉普拉斯、达朗贝尔等人也先后研究过行星、月球和彗星的摄动问题。在月球运动理论方
面,汉森、德洛内、希尔等人的工作尤为突出。汉森用一个大小和形状不变并在空间转动的椭圆
作为中间轨道,分析了月球的轨道摄动;德洛内利用分析力学中的正则变换方法分离出了月球运
动的周期部分,建立了纯分析的月球运动理论,奠定了天体力学变换理论的基础;希尔则以直角
坐标为基本变量,提出一套在旋转坐标系中的月球摄动理论。至 20 世纪初,关于轨道摄动的方法
已不下百种,这些方法同时也促进了分析力学和微分方程理论的发展。
研究各种摄动方法的主要目的是求解摄动运动方程,根据求解原理的不同,这些方法大致可
以分为两类:特殊摄动法和一般摄动法。一般摄动法通常先将摄动加速度表示成小参数的幂级数,
再逐项解析积分求解方程。由于是解析解,故得到的结果包含多种情况的解,能获得有关摄动轨
道的大量信息,这是它的优点,不足主要在于解析推导的过程过于复杂和繁琐,求高阶解时尤其
如此。一般摄动法的另一个优点是能从轨道数据中清楚地揭示出摄动源。比如,1846 年,根据勒
威耶和亚当斯的计算,由天王星的轨道摄动中发现了海王星的存在;1959 年,通过分析“先锋”
1 号卫星的轨道数据证实地球的形状是梨形的。特殊摄动法是通过数值积分求解摄动运动方程,
其优点是可用于计算任意轨道和任意摄动力,不足是积分结果只对应某个特定问题或某组特定的
初始条件。此外,为获得某个时刻的结果,必须计算出所需时刻以前的所有中间时刻的卫星坐标
和速度分量,效率较低,并且存在误差累积问题。特殊摄动法早期在分析小行星、彗星等小天体
的运动中发挥了重要作用,因为这些情形下一般摄动法不再适用。近年来,随着电子计算机的迅
速发展,特殊摄动法在轨道摄动分析、精密星历计算等任务中的应用越来越广泛。
虽然一般情况下,摄动力与中心引力相比非常小,但它们对航天器轨道的影响仍不可忽略。
例如,若忽略地球非球形引力的影响,就无法对人造地球卫星的运动作长期预报,也就无法设计
出满足要求的轨道;如果在星际航行中不考虑其它引力体的摄动,许多飞行任务就可能偏离目标
而无法完成。因此,人类进入航天时代以后,在工作中从天体力学继承了大量的经典摄动分析方
法;另一方面,航天技术的发展也推动了摄动理论的进步。例如,人造地球卫星与地球的距离较
近,地球非球形摄动的影响增强,同时卫星轨道观测的精度大幅提高,这都推动了地球非球形摄
动理论的快速发展;大气阻力、太阳光压等传统研究中很少涉及的摄动因素也纳入学者们的研究
范围。
1
本文档仅用于航天器轨道力学慕课课程参考教材本文档仅用于航天器轨道力学慕课课程参考教材
一. 摄动力分析
§7.1 特殊摄动法
近地航天器绕地球飞行过程中,受到的摄动力主要包括地球非球形引力、大气阻力、日月引
力、太阳辐射压力、小推力作用等。
(1)地球非球形摄动
地球非球形引力通常是影响近地航天器运动最主要的摄动力。由于地球形状的不规则和质量
分布的不均匀,造成存在与航天器的地心矢径相垂直的引力分量,它的作用会使航天器的实际运
动偏离二体轨道。显然,随着航天器地心距离的增大,地球非球形摄动的影响会减弱,因此当航
天器的地心距比地球半径大得多时,就可以用中心引力场来近似地球引力场。
(2)大气阻力摄动
航天器在地球高层大气中高速飞行时,会受到气动力的作用。气动力主要是阻力,作用方向
与航天器相对于大气的速度方向相反,但同时也会产生微弱的升力和侧力。当航天器轨道高度低
于 200km 时,大气阻力是最主要的摄动力。随着轨道高度的增加,气动力的影响急剧减弱。当高
度在 1000km 以上时,大气阻力可以忽略不计。大气阻力会导致轨道高度的衰减,是决定低轨航
天器轨道寿命的重要因素。
(3)日月引力摄动
航天器在地球附近运动时,日月引力是一种典型的第三体摄动力,它是由月球和太阳对航天
器与地球的引力加速度差造成的(见第十三章)。日月引力摄动主要取决于航天器的轨道高度、轨
道形状、轨道面位置和拱线相对于月地、日地连线的位置。航天器轨道高度越低,日月引力摄动
量越小;随着轨道高度的增加,影响也逐渐增大。当高度大于地球同步卫星轨道时,日月引力成
为主要的摄动因素。
(4)太阳光压摄动
根据量子力学理论,光是光子的流动,光子具有动量。当光子碰撞到航天器表面时,会有一
部分被反射回来,被反射的光子把一部分动量传递给航天器,造成辐射压强。因此,当航天器受
到太阳的直接照射时,会产生太阳辐射压强。由于太阳辐射的能量主要集中在可见光波段,因此
太阳辐射压亦称为太阳光压。太阳光压会影响航天器的运动,当轨道高度大于 800km 时,太阳光
压摄动将超过大气阻力摄动。特别对那些面积质量比很大的卫星,光压会对轨道产生实质性的影
响,在摄动分析中不可忽视。
(5)地球形变摄动
地球并不是一个刚体,在日、月引力作用和地球自转不均匀性的影响下,会产生弹性形变。
地球形变会使地球引力场产生附加变化,从而影响航天器受力,称为地球形变摄动。具体而言,
日月引力会导致潮汐现象,称为潮汐摄动,包括固体潮、海洋潮和大气潮;地球自转效应会导致
地球自转形变摄动。
(6)地球反照辐射摄动
地球受到太阳的直接辐射后,会产生次级辐射,包括两部分:一是地球反射光的影响,二是
地球吸收太阳光后转化成热辐射,也就是红外辐射的影响。这两类性质不同的次级辐射能量对航
天器运动同样构成一种摄动作用,其作用机理与光压摄动类似,称为地球反照辐射摄动。由于地
球反照辐射摄动比较小,除特殊需要外,一般可不予考虑。
2
本文档仅用于航天器轨道力学慕课课程参考教材本文档仅用于航天器轨道力学慕课课程参考教材
(7)广义相对论效应摄动
考虑广义相对论效应后,地心惯性系只是一个局部惯性参考系。在此参考系中,航天器的运
动方程不同于牛顿引力场中的形式,其差别相当于航天器受到一个附加摄动,也称为后牛顿效应。
(8)小推力摄动
航天器在轨运行过程中,需要实施姿态和轨道控制,因此会在某些时间段内开启发动机。这
些发动机的推力一般比较小,但也会对轨道运动产生微小影响,称为小推力摄动。
上述各摄动因素中,前四项的量级较大,一般不可忽略;其余项的量级相对较小,可根据精
度要求选择是否考虑。表 7-1 给出了近地航天器所受摄动力与当地中心引力加速度之比的数量级
估计。估计中,阻力系数 Cd = 2.2,面质比 A/m = 0.04m2/kg,太阳反射系数
算采用 US1976 标准大气模型。
,大气密度计
表 7-1 近地航天器所受摄动力的量级
作用力
J2 项
Jn 项(n > 2)
田谐项
大气阻力
太阳光压
日月引力
地球形变摄动
广义相对论效应
量级(g)
轨道高度(km)
10-3
10-5
10-6~10-9
10-6~10-9
10-6
10-10
10-8
10-7
10-7
10-5
10-8
10-9
300
35787
300
300
300
1000
300
35787
300
35787
300
—
在地心惯性坐标系中,考虑摄动力的影响时,航天器质心运动方程可写成如下形式
(7-1-1)
上式是以直角坐标表示的摄动运动方程,其中
是各种摄动力产生的总摄动加速度。根
据摄动力性质的不同,上式可改写成
(7-1-2)
式中 R 称为摄动函数,它包括所有可以用势函数表示的保守力,比如地球非球形引力、日月引力
等; 为耗散力产生的摄动加速度,比如大气阻力、太阳光压力等。根据机械能守恒可知,
不会引起轨道能量或半长轴 a 的长期变化, 则相反。
一般情况下,很难求得摄动运动方程(7-1-1)有限形式的解析解,只能求得数值解或近似级数
解。如前所述,两类求解方法分别称为特殊摄动法和一般摄动法。特殊摄动法通常按被积方程的
形式分类,目前常用的有科威尔法、恩克法和参数变分法,其中参数变分法还是一般摄动法的基
础。
二. 科威尔法
3
0.53prrra1nppkkaa3dRrrradaRda本文档仅用于航天器轨道力学慕课课程参考教材本文档仅用于航天器轨道力学慕课课程参考教材
科威尔法是所有摄动分析方法中最简单和最直接的方法,它是由英国天文学家科威尔(Philip
H. Cowell, 1870-1949)在 20 世纪初提出的,他用此方法计算出了木星的第八颗卫星的轨道,还精
确预测了哈雷彗星在 1910 年的回归。科威尔方法很简单,只要列出所研究对象的运动方程(7-1-1),
在方程的右端包含所有的摄动加速度,然后对运动方程逐步积分即可。随着电子计算机性能的不
断提升,科威尔方法的应用也越来越普遍。
为便于数值积分,将方程(7-1-1)写成一阶微分方程组的形式
其中 r 和 v 是航天器的位置和速度矢量。为进行数值积分,需要将上式写成坐标分量的形式
(7-1-3)
只要列出摄动加速度的解析表达式,就可以用数值积分方法求解上式,得到航天器在任意时
刻的位置和速度。
科威尔方法的优点是公式简单,运算方便,可以同时处理任意多个摄动。但它也有一些缺点,
主要表现在当航天器靠近引力天体运动时,为保证精度,积分步长必须取得很小,这就会大大增
加计算时间和舍入误差的累积。
三. 恩克法
虽然恩克法比科威尔法要复杂得多,但早在 1857 年德国天文学家恩克(Johann F. Encke,
1791-1865)就提出了这种方法,比科威尔法要早半个世纪。恩克将这种方法成功应用于计算短周
期彗星和小行星的轨道。
恩克法与科威尔法的主要区别是:科威尔法将所有的加速度放在一起积分,而恩克方法则对
主要加速度(中心引力加速度)与摄动加速度分别积分。这就意味着有一条基准轨道,没有任何
摄动加速度时,航天器将沿基准轨道运动,这条基准轨道称为密切轨道或吻切轨道(Osculating
Orbit),相应的轨道要素称为密切轨道要素。易知,在理想的中心引力场内,密切轨道就是二体轨
道。
密切轨道具有这样的性质:在任意瞬间,密切轨道与实际轨道在当前点处相切,且航天器在
密切轨道上的速度与实际轨道上的速度相同,如图 7-1 所示。恩克方法就是以密切轨道为基准轨
道,当摄动力与中心引力相比很小时,在较短的时间内,航天器在实际轨道上的位置与在密切轨
道上对应的位置只差一个很小的量。如图 7-2 所示,令 和
分别为 t 瞬时实际轨道与密切
轨道上航天器的位置矢量, 为两者之差。恩克法不直接计算 ,而是计算受摄动后的位置差 。
密切轨道上的位置 可以根据二体轨道公式解析计算,因此当用数值积分方法求得 后,实际的
位置矢量就可以用下式求出
(7-1-4)
若经过一段时间,位置矢量差 变得较大后,则通过校正的方法选择一个新的时刻和新的密切轨
4
3pddtddtrrvvra333,,,xxpxyypyzzpzdvdxvxadtdtrdvdyvyadtdtrdvdzvzadtdtrtrtρrrrρrrρrr本文档仅用于航天器轨道力学慕课课程参考教材本文档仅用于航天器轨道力学慕课课程参考教材
道,使积分继续进行。
图 7-1 密切轨道
图 7-2 恩克方法
下面讨论用恩克方法计算轨道的步骤。由于密切轨道的摄动加速度
,故由式(7-1-1)可得
其轨道方程为
(7-1-5)
上式即二体轨道运动方程。根据密切轨道的定义,在起始时刻 t0 有
设某时刻实际轨道与密切轨道的位置偏差为
可得
将式(7-1-1)、(7-1-5)和(7-1-7)代入上式,可得
(7-1-6)
(7-1-7)
(7-1-8)
(7-1-9)
上式即 的微分方程,从 t = t0 时刻开始对上式积分,即可获得
置和速度。但注意到
是一项几乎相等的两个量之差,为保持计算精度,需要额外增加计
,从而得到航天器的扰动位
算机的位数。为解决此问题,引进一小量 ε
由此可得
于是式(7-1-9)变成
5
(7-1-10)
(7-1-11)
起起起起起起起起起起起起vt0实际轨道密切轨道ρrrt校正点..p0a30ρρ00000tttttrρvrρrrρrrρ33333 1pprrrarrrarrrtr331/r2221r332312r本文档仅用于航天器轨道力学慕课课程参考教材本文档仅用于航天器轨道力学慕课课程参考教材
将上式中方括号内的量用二项式级数展开
(7-1-12)
(7-1-13)
在电子计算机出现之前,上式的计算是一项很繁重的工作。为此发展了一种逼近方法,首先
定义一个新的函数
则式(7-1-12)可改写成如下形式
(7-1-14)
(7-1-15)
由式(7-1-14)求得 q 后,就可以对上式进行数值积分,当然必须先求得 ε 的值。由于在直角坐标系
中有
故有
其中 ρx,ρy,ρz 是 ρ 的分量。式(7-1-16)也可以写成
将式(7-1-17)代入式(7-1-10),整理可得
(7-1-16)
(7-1-17)
(7-1-18)
考虑到实际轨道与密切轨道的位置偏差是很小的,故在上式中可略去 , , ,得到
(7-1-19)
虽然采用式(7-1-19)计算 ε 可能略快些,但在条件允许的情况下,还是推荐采用方程(7-1-13)。q 随
ε 变化的函数表可以事先编制出来,在轨道计算中通过查表获得 q 的值。
一般来说, 的变化要比 r 的慢得多,因此恩克法可以采用较大的积分步长,减少了积分次
数。虽然在每步积分中恩克法花费的计算时间要长一些,但总的来说,它的计算效率要高于科威
尔法。计算表明,计算行星际轨道时大约快 10 倍,计算地球卫星轨道时大约快 3 倍,因为后者
ap 的量级要大一些。在用恩克法计算轨道时还应注意,当
而应重新选择新的密切轨道作为基准轨道。因此,在实际计算中,当
变大后,不能再采用大的积分步长,
大于或等于某个小的常
数(一般可取 0.01)时,就应该进行校正,需要校正是恩克法的一个缺点。
前面给出了两种不同形式的积分方程,对方程积分还要选择相应的数值积分方法。目前已有
多种常微分方程的数值解法在轨道力学中得到了成功应用,这些方法总体上可以分为两大类,单
步法和多步法。单步法仅需一个自变量上的函数值就可以求解,比如 Runge- Kutta 法、Bulirsch-Stoer
外推法等;多步法则需要已知多个自变量上的函数值才能求解,比如 Adams-Bashforth 法、Cowell
法等,这分别是亚当斯发现海王星和科威尔计算哈雷彗星轨道时采用的方法。每种数值积分方法
6
323112prarr32323535711232!3!321112q3pqrarr2222rxyz2222xyzrxyz22222xyzrxxyyzz21222xyzxyzxyz2x2y2z21xyzxyzr/r/r本文档仅用于航天器轨道力学慕课课程参考教材本文档仅用于航天器轨道力学慕课课程参考教材
都有自身的优点与不足,很难单纯地评价哪种方法是最优的,各种方法的具体公式可参见文献
[4]-[6]。在文献[6]中,对各种积分方程和积分算法的效率与精度作了详细比较。
是 500km 高度的太阳同步圆轨道,倾角
利用特殊摄动法,图 7-3 给出了不同摄动力导致的轨道预报误差[7]。仿真中,航天器的轨道
。阻力系数 Cd = 2.2,面质比 A/m = 0.04m2/kg,
,大气模型采用 Jacchia-71 模型,地球引力场模型为 EGM-96 模型。
太阳辐射表面反射系数
计算预报误差时,以计及 12×12 阶地球引力场模型的轨道作为基准轨道,也即除二体轨道、考虑
J2 项以及 70×70 阶引力场模型的三种轨道外,其它摄动轨道的计算都考虑了 12×12 阶的地球引
力场模型。例如,太阳辐射压力(SRP)造成的轨道预报误差是考虑太阳辐射压力和 12×12 阶的
地球引力场模型后的轨道预报结果与基准轨道的偏差。由图可见,地球非球形引力和大气阻力是
对低地球轨道影响最大的摄动因素。
图 7-3 摄动力引起的轨道预报误差
§7.2 参数变分法
参数变分法又称轨道要素变分法,或称常数变易法。1748 年,欧拉首先研究了参数变分法,
1808 年,拉格朗日发表了完善的研究结果。参数变分法在天体力学摄动理论中占据核心地位,至
今仍是分析和计算航天器轨道摄动最重要的方法之一。
参数变分法也是建立在密切轨道的概念基础上。它与恩克法的不同之处在于,恩克法只是使
用某一时刻的密切轨道作为参考基准,在校正之前基准轨道是不变的,而参数变分法的基准轨道
是连续变化的。参数变分法的基本思想是,航天器的实际轨道与密切轨道在密切点处相切,且具
有相同的速度,因此只要建立起密切轨道要素随时间变化的规律,得到任一时刻的密切轨道要素,
航天器的实际运动也就随之确定。Equation Section (Next)
实际上,作为三维空间中的质点运动轨迹,二体轨道可以由任意六个相互独立的参数来描述,
轨道要素只不过是许多可能的参数中几何意义较明确的一组。参数变分法可以选用任意一组参数,
揭示它们在扰动加速度作用下随时间的变化规律,找出参数变化率的解析表达式,并进行积分以
求出未来某个时刻的参数值。本节将分别讨论用开普勒轨道要素和正则参数描述的变分方程。
一. 轨道要素变分方程
7
97.6i0.5012243648607210-410-2100102104106108Time /hDifference /mSRP70x70Thrid-BodyDrag J71J2Solid TideOcean TideTwoBody本文档仅用于航天器轨道力学慕课课程参考教材本文档仅用于航天器轨道力学慕课课程参考教材
建立轨道要素变分方程的方法有两种,一种称为摄动函数法,一种称为力分解法。前者建立
的方程称为拉格朗日型摄动运动方程,后者称为高斯型摄动运动方程。
1. 拉格朗日型摄动运动方程
(1)常数变易法
摄动函数法假定所有的摄动力都能以摄动函数表示,此时航天器的摄动运动方程(7-1-2)具有
如下形式
(7-2-1)
假如不存在摄动力,即 R = 0 时,则上式的解可以写成
(7-2-2)
式中 代表六个轨道要素。对开普勒轨道 为常数,故有
(7-2-3)
根据常微分方程理论中常数变易法的原理,只需要将式(7-2-2)中的轨道要素作为变量,也即看作
时间 t 的函数,其解就能满足方程(7-2-1)。故对式(7-2-2)求微分可得
(7-2-4)
根据密切轨道的特性,在任意瞬间实际轨道与密切轨道对应的位置和速度相等,故上式应同时满
足式(7-2-3),可得到
(7-2-5)
将式(7-2-4)再次对时间微分,并注意到式(7-2-5)成立,可得
(7-2-6)
实际轨道与密切轨道在对应点处的加速度并不相同,对密切轨道有
将上式及式(7-2-1)代入式(7-2-6),可得
(7-2-7)
(7-2-8)
方程(7-2-5)与(7-2-8)组成了关于六个密切轨道要素的一阶微分方程组。若将密切轨道要素对时间
的导数写成显式的形式使用会更加方便,为此将式(7-2-5)和(7-2-8)分别点乘
和
,
然后相减,可以得到
将上式等号左边括号内的项记作
(7-2-9)
8
223dRdtrrr,,1,2,,6itirriiddttrr61iiidddttdtrrr610iiiddtr2226221iiidddtttdtrrr22223ddttrrrr261iiidRtdtr2/jtr/jr2261,1,2,,6iijijijdRjttdtrrrrr本文档仅用于航天器轨道力学慕课课程参考教材本文档仅用于航天器轨道力学慕课课程参考教材