地月系平动点轨道动力学建模与控制策略研究

于登云 周文艳 高珊

(1中国航天科技集团有限公司,北京 100048)(2北京空间飞行器总体设计部,北京 100094)

平动点是第三体在受2个大天体的万有引力作用时在空间中的引力平衡点。1765年,欧拉发现在一个旋转二体重力场中存在3个共线的天平动点;1772年,拉格朗日指出在一个旋转二体重力场中存在另外2个天平动点。这5个点统称为拉格朗日点(L1～L5),也称平动点。其中,L1、L2和L3是共线平动点,L4和L5为三角平动点。平动点特殊的动力学特性和在三体问题中相对固定的几何位置,使其在中继通信、在轨服务、天文观测、星际导航、星际转移等深空探测任务中具备良好的工程应用价值。虽然此前开展了一些研究探索[1-5],但鉴于技术复杂性等因素,未有利用地月系平动点开展地月中继服务的先例。

嫦娥四号是世界首次着陆月球背面进行月球探测的任务,由于月球背面始终背向地球,探测器无法直接与地面通信,需要中继卫星作为地面站与月球背面探测器之间通信的桥梁。为保障中继任务顺利开展,其轨道设计不仅要考虑日-地-月-星相对几何关系,还要考虑各种工程约束条件,如测控要求、运载火箭发射要求、推进剂约束、光照条件等;同时,轨道振幅参数、相位等设计变量与速度增量、轨道光照等约束条件具有极为复杂的耦合关系,在工程设计中要对轨道特性进行分析来解耦这些约束。因此,必须要建立考虑三体引力的高精度动力学模型进行轨道计算、分析,寻找其中的规律和轨道特性,确定满足任务需求的中继卫星最佳轨道方案和最省推进剂的轨道控制策略。

1 轨道动力学建模与求解

1.1 非线性动力学通用模型

在惯性空间中,n个物体在无其他外力的作用下,公共质心的加速度为零,各质量体在相互作用下绕质心作圆周运动,示意见图1。设n个物体的质量分别为m i(i＝1,…,n),它们相对于惯性坐标系OXYZ的矢径分别为r i(i＝1,…,n),惯性坐标系的原点为它们的公共质心O。

图1 n体问题示意Fig.1 Illustration of multi-body problem

第i个物体的单位质量加速度为

式中:G为引力常数;第i个物体到第j个物体的矢径r ij＝r j－r i。

一般情况下,要研究卫星相对于某个主物体的运动,坐标系原点在主物体质心。把第n个物体当作卫星,第1个物体当作主物体,则卫星相对于主物体惯性坐标系中的加速度为

这就是n体运动方程。研究卫星在平动点的运动是一个四体问题,要包括日、地、月的中心引力,卫星在地心惯性坐标系下的运动方程如下。

式中:r和v为卫星的位置矢量和速度矢量;rm和rs为地心指向月球和太阳的矢量;rmd和rsd为月心和日心指向卫星的矢量;μe,μm,μs为地球、月球和太阳的引力常数;J2为地球带谐系数;矢量r z可表示为[0 0z]T;地球平均赤道半径Re＝6378.140 km。

1.2 微分修正方法

卫星在不同运行阶段的计算是选取不同的初始和末端状态,利用状态转移矩阵进行迭代求解,从而找到精确的结果[6-7]。初始状态参数矩阵p和末端状态参数矩阵q之间的关系可以表示为q＝f(p),两状态之间的误差关系可线性化表示为Δq＝利用这个关系,在合适的初值基础上,通过若干次迭代就可以求出十分精确的数值解。Δq就是误差传递矩阵,它的数值解可以通过数值积分得到,见式(4)。

式中:r0和v0为卫星在初始时刻的位置和速度;3个矩阵S,R,T的乘积结果为误差传递矩阵Ts,其中S和T分别为末端和初始状态转移矩阵,R为状态转移矩阵,见式(5),其微分形式见式(6)。

式中:g为。

首先,给出轨道积分的初值,将式(3)和式(6)一起积分,得到任意时刻卫星的位置、速度和状态转移矩阵R,再利用轨道积分的初值求出矩阵T,轨道积分到末端状态求出矩阵S,这3个矩阵相乘就得到了误差传递矩阵Ts;然后,把积分得到的末端状态减去要求达到的值,差值就是误差传递矩阵Δq,则初始状态的修正量Δp＝Ts－1Δq。用这个修正量来修正初始状态,然后再重复上述积分,这样反复迭代到末端状态达到一定的精度,最终求解出满足末端状态的轨道。

1.3 解析计算模型

考虑3个物体的动力学系统(即三体系统),其中2个为主引力天体,第3个为卫星,其质量远小于前2个天体,可认为前2个天体的运行不受卫星影响。进一步假定2个天体围绕该系统的公共质心作匀速圆周运动,则称为圆型限制性三体问题。以地月系为例描述这一问题,建立地月质心会合坐标系OBXBYBZB,其坐标原点位于三体系统的公共质心,XB轴为大天体指向小天体方向,ZB轴沿2个天体运行的角动量方向,YB轴满足右手坐标系定义。在此基础上定义平动点坐标系OLXLYLZL,将质心会合坐标系平移至相应的平动点,其余指向不变。

图2以地月系L2为例给出地月质心会合坐标系

OBXBYBZB、平动点坐标系OLXLYLZL和地心惯性坐标系OXYZ的关系。平动点坐标系(地月L2会合坐标系)的原点在L2,L2与地月系公共质心的距离为1.155682 110 911D(D为地月距离,由星历表DE405得到),L2与月心的距离为0.167 832 682 341D;地月L2会合坐标系的XL轴为瞬时地月矢量方向,ZL轴为白道面法向。

图2 3个坐标系的关系Fig.2 Relationship among three coordinate systems

圆型限制性三体问题动力学模型普遍表达为

式中:xB,yB,zB为卫星在地月质心会合坐标系中的位置分量;拟势函数U＝(1－μ)/r1B＋μ/r2B＋(xB2＋yB2)/2,其中μ为三体系统的2个天体中较小天体的质量占比,r1B和r2B分别为地月质心会合坐标系中2个天体的矢径大小;U xB,U yB,U zB表示拟势函数U对xB,yB,zB求偏导。

平动点轨道类型有很多,最常见和工程上已经实现的是Halo轨道和Lissajous轨道。Halo轨道分为南向轨道和北向轨道,以地月系为例,判断Halo轨道的最简单准则是,当卫星经过地月视线时,其速度近似平行于白道面并垂直于地月视线。Lissajous轨道的通常定义是,平面内振动频率和垂直平面的振动频率不同。平动点轨道还有Lyapunov轨道、垂直轨道、轴向轨道、蝶形轨道等类型,关于各种类型平动点轨道的定义,目前没有统一[8]。本文结合“鹊桥”卫星的任务特点,主要围绕Halo轨道开展研究,对于Lissajous轨道模型不再赘述,具体可参见文献[9]。下面主要针对Halo轨道的解析解进行介绍。

三体动力学方程无法直接得到解析解,但可以利用勒让德(Legendre)多项式将平动点附近的运动方程线性化近似,再通过求解线性化方程组给出周期轨道的近似解析解。将式(7)所采用的地月质心会合坐标系平移至相应的平动点(如L2),指向保持不变,并且进行无量纲化和归一化处理,可以得到第三体在共线平动点坐标系下的运动方程形式为

考虑运动方程中的三阶项影响,利用Legendre级数得到保留三阶小量的方程,见式(9)。采用Lindestedt-Poincaré方法,可进一步推导出Halo轨道的三阶近似解析解,见式(10)。

式中:A XL和A ZL分别为XL轴和ZL轴方向的振幅;τ1为相位角;a,b,k,d均为推导因子;zL取“＋”时为南向Halo轨道,取“－”时为北向Halo轨道。

需要说明的是,式(10)中利用Legendre级数对三体动力学方程进行展开,在平动点附近时,因高阶项均为小量o(x n),可以在近似计算中舍去。然而,当第三体位置与平动点的距离较大时,高阶项中因包含该距离及各方向分量的高次项,不再是小量,若采用舍去高阶项的近似处理,则由高次项后产生的差异会使解析近似计算结果与实际情况偏差较大。

2 中继卫星运行轨道的确定

2.1 飞行方案简介

“鹊桥”卫星飞行阶段分为发射段、地月转移段、月球到L2转移段、捕获段和使命段。

(1)发射段:从运载火箭发射到星箭分离的飞行阶段。

(2)地月转移段:从星箭分离到近月点的飞行阶段,1天2个发射窗口的转移轨道飞行时间分别按照112 h和136 h设计。

(3)月球到L2转移段:卫星到达近月点时进行近月制动,制动的速度增量不足以使卫星捕获进入月球轨道,而是进入稳定流形,飞行3～4天到达L2附近,即到达经过地月连线并垂直于白道面的平面。

(4)捕获段:卫星到达L2后的第1圈轨道。捕获段通过轨道捕获控制形成满足任务要求的使命轨道。

(5)使命轨道段:捕获段结束后,卫星进入使命轨道并开展中继工作,预期工作寿命为3年。

在轨道设计过程中,需要首先根据任务需求和约束确定使命轨道,然后以使命轨道为目标,针对从地球出发的转移轨道的任务参数和轨道控制策略进行设计。下文给出具体设计过程和结果。

2.2 使命轨道的确定

中继轨道可选择方案有月球冻结轨道(环月轨道)、地月L2轨道,以及介于两者之间的近直线晕轨道(NRHO),见图3。环月中继轨道可采用大椭圆冻结轨道,这样几乎不需要进行轨道维持,但是无法保证月球背面着陆点的全时段覆盖,存在通信盲区。地月L2中继轨道通过合理选择着陆区及轨道构型,可实现月球背面着陆点的全时段通信覆盖。NRHO也是近年来备受关注的一种轨道类型,其主要优点是通过设置轨道的共振周期避免月影,但也无法实现月球背面着陆点的全时段通信覆盖。3种中继轨道通信覆盖对比如表1所示。可以看出:以NRHO共振周期与月球会合周期比例9∶2为例,NRHO对嫦娥四号着陆区也有较好的覆盖性,但对近月点幅角要求约90°,直接地月转移轨道无法满足幅角要求,需要采用低能转移(弱稳定边界转移)的方式,转移时间3～4个月;另外,该类型轨道对运载火箭出发能量要求较高。通过对比,中继轨道选择地月L2轨道。

图3 中继轨道类型Fig.3 Types of relay orbit

表1 中继轨道通信覆盖对比Table 1 Comparison of communicationscoverage of relay orbits

对于平动点轨道类型的选择,Lissajous轨道在会合坐标系下呈Lissajous曲线的形式,在长时间飞行的情况下难以避免月掩的出现,对关键事件测控弧段的安排带来不利。因此,为保障长期运行对地通信的需求,使命轨道类型选择Halo轨道。由于嫦娥四号任务着陆区域位于南半球,为了使中继卫星能够有更多时间处于南半球上空,以利于中继任务的空间几何关系,因此使命轨道选择南向Halo轨道。

利用解析方法,可以计算出不同形状Halo轨道的位置和速度,见图4。一些典型的振幅对应的远月端和近月端的位置和速度见表2[10]。

图4 Halo轨道的位置和速度Fig.4 Halo orbit position and velocity

表2 不同振幅下的Halo轨道的位置和速度Table 2 Halo orbit position and velocity of different amplitudes

对2018年5－6月各窗口使命轨道本影时长进行统计,使命轨道的运行时长为3年,结果见表3。由表3可知:对于2018年5－6月的发射窗口,A ZL为13 000 km的轨道最长本影时长为4.7 h。

“鹊桥”卫星相对于嫦娥四号着陆器的距离和方位角的极坐标及仰角变化见图5。可以看出:对于主选着陆区中心点,使命轨道段运行3年内,“鹊桥”卫星相对着陆器的仰角均大于20°。

综合考虑光照、测控、速度增量等约束,L2环绕轨道选择A ZL为13 000 km的南向轨道。Halo轨道的基本参数为:A XL约为12 500 km,A YL约为37 000 km,A ZL约为13 000 km,轨道的平均周期约为14天。其轨道形状如图6所示。考虑轨道演化和误差的影响,使命轨道运行3年的振幅变化范围为:A XL不超过14 000 km,A YL不超过38 000 km,A ZL不超过17 000 km。该使命轨道与着陆器通信覆盖率达100%,且对于首发窗口在3年飞行过程中只有2次经历地球的本影。

表3 Halo轨道本影时长Table 3 Eclipse time of Halo orbit

图5 “鹊桥”卫星方位角和仰角Fig.5 Azimuth and elevation angle of Queqiao satellite

图6 Halo轨道形状Fig.6 Illustration of Halo orbit

2.3 转移轨道的确定

实现地月L2轨道飞行的转移轨道,有直接转移、低能转移和月球引力辅助转移3种方案可选,见图7。

直接转移轨道与地月转移轨道类似,只是其远地点更高,因此出发能量比地月转移轨道稍大。直接转移轨道在到达L2时要进行捕获控制,所需速度增量约1.0 km/s。低能转移轨道即行星际高速公路,通过从地球发射进入抛物线轨道,使卫星运行到距离地球1.5×106km左右后再返回地月系,进入地月L2轨道,其所需速度增量小,但是近地点出发速度较大,且转移时间一般为3～4个月。月球引力辅助转移轨道是将卫星从地球发射进入地月转移轨道,卫星飞行至近月点时通过近月制动完成减速变轨,所需速度增量约为0.2 km/s。卫星并不进入环月轨道,而是在飞越月球后进入地月L2附近的稳定流形,在到达L2附近时可不经变轨或只需很小的速度增量就进入环绕L2的轨道。3种转移轨道参数比较见表4。

图7 3种转移轨道Fig.7 Three transfer trajectories

表4 转移轨道对比Table 4 Comparison of transfer trajectories

月球引力辅助转移轨道继承性好,转移时间短,速度增量需求远小于直接转移轨道,而且考虑到发射能量对发射质量的影响,“鹊桥”卫星的剩余质量大。

2.3.1 地月转移段

“鹊桥”卫星的地月转移段设计充分继承嫦娥三号探测器的设计方法和特点,由长征-4C运载火箭送入近地点高度200 km、远地点高度380000 km、倾角28.5°降轨发射的地月转移轨道,经过5～6天的飞行,降轨到近月点高度100 km、轨道倾角15°。在设计过程中,根据任务需求对到达倾角进行优化,飞行时间根据近月制动时的测控条件在112 h和136 h附近进行调整。

2.3.2 月球到L2转移段

初始状态参数矩阵p和末端状态参数矩阵q在月球到L2转移段可表示为[6]

式中:Δv为近月点制动速度增量,卫星在地月转移段运行到近月点时,位置和速度是确定的,此时在近月点沿卫星相对于月球的速度反方向进行减速制动,使卫星进入地月L2稳定流形;v XL为卫星运行到地月L2会合坐标系OLXLYLZL的XLOLZL平面时XL轴方向速度。

Δv的初值为,其中,μm为月球引力常数,rp为卫星从月球出发时的月心距,月球到L2转移段的剩余速度平方C3取－0.175 km2/s2。

迭代计算时采用逐次迭代收敛的方法,具体过程如下。

(1)以近月制动速度增量为控制变量,目标为第1次经过地月L2会合坐标系XLOLZL平面时,相对于会合坐标系的速度v XL＝0,迭代收敛;

(2)继续以近月制动速度增量为控制变量,目标为第2次经过地月L2会合坐标系XLOLZL平面时,相对于会合坐标系的速度v XL＝0,迭代收敛;

(3)继续对轨道进行积分,目标为第3次经过地月L2会合坐标系XLOLZL平面时,相对于会合坐标系的速度v XL＝0,迭代收敛。

经过3次迭代收敛,求出近月制动速度增量。卫星按该增量减速制动后,从月球飞至地月L2,进入L2的稳定流形,近月点高度100 km,轨道倾角15°,地月转移时间112 h。不同发射窗口下近月制动速度增量见图8。可以看出,近月制动速度增量1个月变化1个周期,与月球公转周期相同。

图8 不同发射窗口的近月制动速度增量ΔvFig.8 Perilune decelerationΔv of different launch windows

3 中继卫星轨道控制策略

“鹊桥”卫星整个飞行过程需要进行的轨道控制可分为捕获控制、轨道修正和轨道维持3类。捕获控制的目的是使卫星捕获进入预定的使命轨道;轨道修正的目的是对飞行过程中的各类误差(包括运载火箭入轨误差、测控误差、轨道控制误差等)进行修正,使在误差情况下飞行轨道仍能满足任务要求;轨道维持则是由于地月系L2是不稳定点,在使命轨道飞行阶段需要定期进行轨道维持以保证卫星在Halo轨道上的长期稳定运行。

3.1 捕获控制

根据轨道特性和地月L2的Halo轨道构型需要,“鹊桥”卫星通过1次近月制动和2次Halo轨道捕获控制,变轨进入特定构型的Halo轨道(南向、A XL为12500 km、A YL为38 000 km、A ZL为13 000 km),并通过合理安排2次捕获脉冲的位置及Halo轨道进入点的相位,使3次脉冲的速度增量之和最小。

“鹊桥”卫星近月制动后进入月球到L2转移段,为确保在到达L2附近后能够顺利进入环绕L2的轨道并开展捕获控制,近月制动及后续转移段的控制机动均应以保障到达L2附近后第1圈轨道的稳定为目标。近月制动采用4个20 N推力器加8个5 N推力器组合执行,推力方向沿近月点速度反方向,采用惯性定向,1次完成,通过控制推力作用时间实现预定轨道控制目标。

“鹊桥”卫星到达地月L2附近后的第1圈为捕获段,进行2次捕获控制,用于调整轨道位置和速度,形成满足中继任务的使命轨道;同时,由于捕获控制速度增量较大,为保证轨道稳定性,在各次捕获控制后分别安排1次轨道修正,用于修正捕获控制残差。

第1次捕获控制为“3对2”的控制策略,即通过调整变轨点3个方向的速度分量瞄准第3次穿过地月L2会合坐标系XLOLZL平面时XL轴方向和ZL轴方向位置。该策略解不唯一,通过对控制变量进行优化获得推进剂最优解。

第2次捕获控制为“3对1”控制策略,需要首先以Halo轨道过XLOLZL平面时地月L2会合坐标系下各方向速度分量要求(对于A ZL为13000 km的南向Halo轨道,过远月端XLOLZL平面时v XL为0 m/s,v ZL为0 m/s,v YL为－165.40 m/s)为初始目标,计算满足各方向速度所需要的变轨速度增量,并将计算结果作为初值,通过迭代调整变轨点各方向速度,使控制后第3次穿过地月L2会合坐标系XLOLZL平面时沿XL轴方向的速度为零。

3.2 轨道修正

地月转移段中途修正安排3次,综合考虑各发射窗口测控弧段差异和后续可能的微调,第1次安排在星箭分离后的14 h±1 h,第2次安排在星箭分离后的38 h±1 h,第3次安排在到达近月点前24 h±1 h实施。地月转移中途修正为“3对3”的微分修正策略,由变轨速度增量3个方向分量瞄准近月点高度100 km、轨道倾角15°和真近点角0°的轨道控制目标。

在月球到L2转移段,根据定轨结果和测控弧段安排2次中途修正,分别在近月制动后24 h和48 h,用于修正近月制动残差。中途修正沿速度方向施加速度增量,轨道控制目标为第3次穿过地月L2会合坐标系XLOLZL平面时沿XL轴方向的速度为零。

“鹊桥”卫星在每次捕获后约2天进行1次轨道修正,具体修正时刻根据实际测控条件决定,并可根据每次捕获后的误差情况决定是否取消。轨道修正的目标与前一次捕获控制的目标一致。

3.3 轨道维持

“鹊桥”卫星在L2工作期间,计划每周进行1次轨道维持,3年共150次。轨道维持安排在“鹊桥”卫星过地月L2会合坐标系XLOLZL平面前后,具体时刻可根据国内地面站测控条件确定。轨道维持策略是通过调整控制时刻各方向速度增量,实现第3次穿过地月L2会合坐标系XLOLZL平面时沿XL轴方向速度为零的目标。

4 应用效果

“鹊桥”卫星于2018年5月21日由长征-4C运载火箭直接送入地月转移轨道,经过1次地月转移中途修正、1次近月制动、1次月球到L2转移中途修正、2次Halo轨道捕获控制,于6月14日进入距离月球约65 000 km的地月L2的Halo使命轨道,成为世界上首颗运行在地月L2的Halo轨道卫星,飞控实施的飞行轨迹见图9[11],与设计结果一致。

图9 “鹊桥”卫星飞行轨迹Fig.9 Trajectory of Queqiao satellite

进入使命轨道后至2019年4月初,“鹊桥”卫星围绕L2飞行20圈,共实施32次轨道维持。轨道维持情况见表5。针对已有的维持实施情况分析维持频率和速度增量消耗,“鹊桥”卫星轨道维持频率平均约9天1次,1年对应的速度增量约为18 m/s,满足前期的速度增量预算。卫星整个飞行过程中的参数和轨道控制效果符合预期,速度增量满足前期预算设计结果,充分验证了卫星轨道设计和轨道控制策略设计的正确性和合理性。

“鹊桥”卫星轨道设计开创了地月平动点Halo轨道设计的先河,也是国际上首次工程实现。与其他月球探测轨道设计不同,它涉及到发射轨道和地月转移轨道的拼接、地月转移轨道和月球到L2转移轨道的拼接、月球到L2转移轨道和L2使命轨道的拼接,设计过程复杂。另外,L2轨道的不稳定,在轨道计算时,初值、计算精度都会影响计算的收敛性,仅使用状态转移矩阵计算很难收敛。“鹊桥”卫星的轨道设计结果满足工程任务需求,为中继任务的成功实施提供了保证,并经过了飞行试验的验证,设计模型和设计方法在后续平动点相关任务中具有良好的应用前景。

表5 轨道维持统计Table 5 Summary of orbit stationkeeping