2. 北京航天飞行控制中心,北京 100094;
3. 探月与航天工程中心,北京 100037
2. Beijing Aerospace Control Center,Beijing 100094,China;
3. Lunar Exploration and Space Program Center,Beijing 100037,China
“嫦娥4号”探测器将实现人类首次月球背面软着陆并开展巡视探测。由于月球的自转周期和绕地球的公转周期同步,月球始终只有一面朝向地球,在地球上无法看到月球背面,在月球背面着陆的探测器也无法与地球直接通讯,因此要实现月球背面软着陆探测必须首先解决通讯问题。
1765年欧拉(Euler)发现了在一个旋转二体重力场中存在3个共线的天平动点,1772年拉格朗日(Lagrange)指出在一个旋转二体重力场中存在另外两个天平动点,后人将这5个点统称为拉格朗日点,也称平动点[1-2]。平动点是第三体在受两个大天体的万有引力作用时,在空间中的引力平衡点。运行于平动点的飞行器可以保持该双星系统的公转角速度而几乎不用消耗推进剂。由于平动点特殊的动力学特性和在三体问题中相对固定的位置,使其在停泊中转、中继通信、星际转移等未来深空探测任务中具备较好的工程应用价值。地月L2平动点位于地球至月球连线的延长线上,与地球、月球的位置相对固定。中继星在地月L2点附近可进行长期环绕飞行,对月球背面进行连续观察。因此,发射一颗环绕地月L2点的中继星,即可建立地球与月球背面探测器的中继通信链路,从而有效支持月球背面着陆探测任务。
“嫦娥4号”中继星将是我国首颗正式任务采用平动点轨道的中继卫星,为合理地设计飞行任务轨道,必须对从地球到地月L2点的转移轨道和L2点附近的使命轨道进行比较全面的研究。本文对从地球出发,利用月球引力辅助变轨,形成地月L2点的飞行轨道进行研究,分析了发射窗口、地月转移时间、近月点高度、近月点倾角和轨道振幅等因素对转移轨道和使命轨道特性的影响,寻求满足月球探测L2点中继任务需求的飞行轨道。
1 研究方法 1.1 飞行方式从地球到地月L2点的转移轨道通常有3种方式:直接转移、低能转移和月球引力辅助变轨转移[3-5],飞行轨迹示意图分别对应图 1的3条曲线。
直接转移轨道与地月转移轨道类似,只是其远地点更高,因此出发能量比地月转移轨道稍大,直接转移轨道在到达L2点时要进行捕获控制,所需速度增量约1 km/s;低能转移即行星际高速公路,通过从地球发射抛物线轨道,使卫星运行到日地平动点附近后返回地月系统,进入地月L2点轨道,低能转移轨道所需消耗的速度增量小但是近地点出发速度较大,且转移时间一般长达3~4个月[4,6-7]。月球引力辅助变轨转移方式是从地球发射直接地月转移轨道,飞行至近月点时通过月球借力完成减速变轨,所需速度增量约200 m/s,卫星不进入环月轨道,而是在飞越月球后进入地月L2点附近的稳定流形,在到达L2点附近时可不经变轨或只需很小的速度增量就进入环绕L2点的轨道[3-5]。3种转移方式的参数比较见表1。
本文主要对第3种飞行方式进行研究,这种方式继承性好、转移时间短、速度增量需求远小于直接转移,而且考虑到发射能量对发射质量的影响,探测器的剩余质量大。这种方式的转移轨道包括两段,一段是地月转移,另一段是从月球到地月L2点的转移。
在地月L2点飞行的轨道(平动点轨道)通常有Lissajou轨道和halo轨道。Lissajous轨道在地月旋转平面内和平面垂直方向均有振动,且平面和垂直方向的振动频率不同,而halo轨道在平面和垂直方向的振动频率相等,目前工程应用中已实现的主要是这两种轨道。平动点轨道还有其他几种类型,如垂直轨道、轴向轨道、蝶形轨道等,本文不进行讨论。
1.2 运动方程惯性空间中,n个物体在无其它外力的作用下,公共质心加速度为零,各物体在相互作用下绕质心作圆周运动。设n个物体的质量分别为
第i个物体的单位质量的加速度为
${\ddot r_i} = G\sum\limits_{\begin{array}{*{20}{c}} {j = 1}\\ {j \ne i} \end{array}}^n {\frac{{{m_j}}}{{r_{ij}^3}}{r_{ij}}} $ | (1) |
其中rij是第个i物体到第j个物体的矢径。
将第n个物体作为飞行器,第1个物体作为主物体,坐标系原点在主物体质心,则飞行器相对于主物体惯性坐标系中的加速度为
$\begin{aligned}{{{{\ddot r}}}_{1n}} & = {{{{\ddot r}}}_n} - {{{{\ddot r}}}_1} \\ &\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! = - G\frac{{\left( {{m_1} + {m_n}} \right)}}{{r_{1n}^3}}{{{r}}_{1n}} + \sum\limits_{j = 2}^{n - 1} {{m_j}\left( {\frac{{{{{r}}_{{{nj}}}}}}{{r_{nj}^3}} - \frac{{{{{r}}_{{{1j}}}}}}{{r_{1j}^3}}} \right)} \end{aligned}$ | (2) |
式(2)为n体运动方程。对于飞行器运动是四体问题时,考虑日、月的中心引力,根据式(2)和图 2,可以写出如下飞行器运动方程
$\left\{ \begin{array}{l}\!\!\! {{\dot r}} = {{v}} \\\!\!\! \begin{aligned}{{\dot v}} = & - \frac{{{\mu _{\rm e}}}}{{{r^3}}}{{r}} - {\mu _{\rm m}}\left( {\frac{{{{{r}}_{\rm{md}}}}}{{r_{\rm{md}}^3}} + \frac{{{{{r}}_{\rm m}}}}{{r_{\rm m}^3}}} \right) - {\mu _s}\left( {\frac{{{{{r}}_{\rm{sd}}}}}{{r_{\rm{sd}}^3}} + \frac{{{{{r}}_{\rm s}}}}{{r_{\rm s}^3}}} \right)\\ & - \left( {\frac{{3{\mu _{\rm e}}{J_2}R_{\rm e}^2}}{{2{r^5}}} - \frac{{15{\mu _{\rm e}}{J_2}R_{\rm e}^2{z^2}}}{{2{r^7}}}} \right){{r}} - \frac{{3{\mu _{\rm e}}{J_2}R_{\rm e}^2}}{{{r^5}}}\left( {\begin{array}{*{20}{c}}\!\!\! 0 \!\!\! \\\!\!\! 0 \!\!\! \\\!\!\! z \!\!\!\end{array}} \right)\end{aligned}\end{array} \right.$ | (3) |
其中:μe、μm和μs是地球、月球和太阳的引力常数;J2为地球带谐系数;Re是地球平均赤道半径。
1.3 坐标系在轨道的积分和计算中用到了两个坐标系:地心J2000坐标系和地月L2会合坐标系。地心J2000坐标系的原点在地心,三轴符合IAU2000的定义;地月L2点会合坐标系,原点在地月L2点,
限制性问题方程首先由欧拉提出,行星绕太阳运动近似为圆、小行星相对于这些大天体的质量很小。将一般的三体问题简化为限制性问题的基本假设是两个物体是主体(用P1和P2表示),相对于第三体(用P3表示)的质量非常大,第三体相对于主体的质量很小,不影响主体的运动。如果忽略P3的引力,主体P1和P2的运动简化成二体问题的解。P1和P2的运动是相对于它们质心的开普勒运动。二体运动进一步简化为主体相对质心的圆轨道,P3的运动即变成圆限制性三体问题的解。
如图 3所示定义坐标系xyz,以P1、P2二体的质心为坐标原点,x由大天体P1指向P2,y位于P1、P2的运动平面内并垂直于x,z垂直于运动平面。记P3在该坐标系下的坐标为(x,y,z)。对单位进行归一化,P3二阶微分运动方程的标量形式如下[8]
$\ddot x - 2\dot y - x = - \frac{{\left( {1 - \mu } \right)\left( {x + \mu } \right)}}{{r_1^3}} - \frac{{\mu \left( {x - \left( {1 - \mu } \right)} \right)}}{{r_2^3}}$ | (4) |
$\ddot y + 2\dot x - y = - \frac{{\left( {1 - \mu } \right)y}}{{r_1^3}} - \frac{{\mu y}}{{r_2^3}}$ | (5) |
$\ddot z = - \frac{{\left( {1 - \mu } \right)z}}{{r_1^3}} - \frac{{\mu z}}{{r_2^3}}$ | (6) |
定义拟势函数U
$U = \frac{{\left( {1 - \mu } \right)}}{{{r_1}}} + \frac{\mu }{{{r_2}}} + \frac{1}{2}\left( {{x^2} + {y^2}} \right)$ | (7) |
$\left\{ {\begin{array}{*{20}{c}}{\ddot x - 2\dot y = {U_x}}\\{\ddot y + 2\dot x = {U_y}}\\{\quad\quad\ddot z = {U_z}}\end{array}} \right.$ | (8) |
当拟势函数的偏导数(8)等于零时,即
共线的3个天平动点由方程
$x - \frac{{\left( {1 - \mu } \right)\left( {x + \mu } \right)}}{{{{\left| {x + \mu } \right|}^3}}} - \frac{{\mu \left( {x - \left( {1 - \mu } \right)} \right)}}{{{{\left| {x - \left( {1 - \mu } \right)} \right|}^3}}} = 0$ | (25) |
的解确定,对于地月系统x1= 0.837,x2= 1.156,x3= –1.005。L4和L5位于r1=r2= 1的点
${{\rm L}_4} = \left[ {\begin{array}{*{20}{c}}{\displaystyle\frac{1}{2} - \mu }\\{\displaystyle\frac{{\sqrt 3 }}{2}}\\0\end{array}} \right],\;{{\rm L}_5} = \left[ {\begin{array}{*{20}{c}}{\displaystyle\frac{1}{2} - \mu }\\{ - \displaystyle\frac{{\sqrt 3 }}{2}}\\0\end{array}} \right]$ |
对平动点轨道及转移轨道的求解除需要考虑太阳、地球、月球外,还要考虑太阳光压等的影响,没有解析解,只能采用数值积分对轨道进行精确计算。对于地月转移轨道采用文献[11]中提出的对轨道进行数值积分,然后用误差传递矩阵修正初值的方法找到地球出发的初始状态;对于月球到L2转移轨道采用二分法调整近月点的初值,找到满足地月L2点状态要求的近月点初始速度;而平动点轨道的求解是采用解析方法求得需要的轨道位置和速度,然后对轨控位置进行优化达到目标轨道,在平动点的长期运行及轨道维持同样采用二分法计算轨控的速度增量[12-13]。
2 轨道设计 2.1 使命轨道本节针对使命轨道设计过程中的任务几何、轨道振幅、阴影条件及轨道维持等问题进行探讨。
2.1.1 任务几何由于中继星承担着着陆器与地球之间的通讯任务,因此使命轨道设计中应充分考虑到月掩、着陆器可见性等任务几何条件对中继任务的影响。
1)月掩
掩星是指一个天体在另一个天体和观测者之间通过而产生的遮蔽现象,由月球遮蔽引起掩星现象称为月掩。当中继星在L2点运行时,可能会面临由于月掩而导致的无法实现中继任务的问题,如图 4所示。
在轨道设计时,应通过轨道类型和振幅的选择尽量避免月掩情况的出现。对于Lissajous轨道,其轨道在会合坐标系下呈Lissajous曲线的形式,虽然可以通过选取初始相位和轨道振幅使航天器尽量避开月掩带,但在长时间飞行的情况下仍然难以避免月掩的出现。因此考虑到保障长期运行对地通讯的需求,中继星轨道类型采用halo轨道。
halo轨道在会合坐标系下的各向振幅分别记为Ax、Ay、Az,三者单调相关,利用解析方法,可以计算出不同形状halo在近月端、远月端位置以及各向振幅大小,见图 5。
不同振幅的轨道,其任务几何条件和光照条件均有一定差异,因此halo轨道振幅的选择需考虑中继任务需求、光照条件以及轨道维持消耗等因素综合确定。为轨道描述方便,后文中统一选用Az作为轨道特征变量。
针对不同振幅下的halo轨道进行了MEP角(月心–地心–中继星夹角)、MFP(月心–测站–中继星夹角)进行了分析,其中测站以喀什站为例进行了计算。对于MEP角,Az振幅为1.2万 km的轨道不小于1.2°,Az振幅为1.3万 km的轨道不小于1.3°,Az为1.5万 km的轨道不小于1.5°。对于MFP角,Az振幅为1.2万 km、1.3万 km和Az= 1.5万 km的轨道分别不小于1.3°、1.4°和1.6°。按照测站与月球的距离关系,当MFP角大于0.3°时,月球不会遮挡测站与中继星的通信。
根据分析可知,振幅越大的halo轨道,MEP角和MFP越大,出现月掩的概率越小,对于Az= 1.2万~1.5万 km的轨道,均不会出现月掩问题。
2)着陆器可见性
针对不同振幅轨道的中继星相对着陆器的高度角进行了分析,对于Az= 1.2万 km的轨道,高度角范围为21~64°;对于Az= 1.3万 km的轨道,高度角范围为20~65°;对于Az= 1.5万 km的轨道,高度角范围为18~66°。按照着陆器仰角5°考虑,统计了不同振幅轨道下中继星运行3年内的几何可见性及占总时间的百分比,结果如表 2所示。可以看出,使命轨道振幅越大,中继星相对着陆器的最小高度角越低,而着陆器对中继星的几何可见性也越差。
因此,在进行使命轨道设计时,应当结合任务需求,综合考虑月掩和着陆器可见性因素来选择任务轨道振幅。
2.1.2 阴影分析地月L2点附近的阴影主要是由于受到地球或月球遮挡而造成的,效果如图 6所示。由于平动点空间位置的特殊性,运行在其附近的航天器光照条件与环月轨道差异很大,其阴影时长往往较长。
从图中可知:即使日–地–月不在同一直线上(即不是月食时),也有可能对平动点轨道产生阴影影响。根据太阳、地球、月球的形状位置参数,计算在地月L2点YZ平面(经过L2点且垂直地月连线的平面)的地月阴影投影半径如表 3所示。
相较万千米级别的中继轨道包络振幅,阴影半径已低一个数量级,阴影轨迹在中继轨道包络投影范围内呈窄带状分布。选取南北向不同振幅(Az= 5 000 km、1.15万 km、1.5万 km)的各3条轨道,一共6条halo轨道,分析不同初始相位下轨道在3.5年内的地月影情况,地月阴影的计算均选取本影。统计结果如表 4所示。
通过对阴影情况的统计分析可以得出以下结论:
1)地影平均时长2 h左右,且地影时长不随振幅、相位、南北向的改变而发生较明显的变化,3.5年长期飞行过程中无法完全规避。
2)halo轨道振幅越大,恶劣月影时长越短。
3)恶劣月影时间长,有可能威胁航天器安全。
4)月影时长随轨道初始相位改变较大,在一定的初始相位范围内,可以实现3.5年长期飞行过程中无月影。
5)3种振幅在月影方面均表现出南向优于北向的情况,这与所选分析阶段的天体空间分布有关。理论上,南北向的长期阴影分布情况应相同。
6)分析尚未考虑地月全影连续遮挡,但从halo轨道和日地月的空间分布可以初步判断,地月全影连续遮挡的情况比较罕见。
2.1.3 轨道维持由于地月L2点为不稳定点,且中继星在飞行过程中受到多种摄动和误差因素的影响,因此需要定期进行轨道控制才能保证中继星在L2点附近的长期稳定运行。halo轨道的维持选择在过旋转坐标系的XOZ面附近进行,根据测控系统对变轨测控弧段的需求,可在过XOZ面前后进行适当调整。维持控制方式考虑以下两种[14-15]。
1)拟halo轨道方式,维持目标是轨道在halo轨道附近,与halo轨道面的幅值偏差保持在一定范围内;
2)halo轨道方式,其维持目标是严格的halo轨道。
以各向速度误差0.01 m/s为例,计算不同控制方式和频率对应的速度增量如表 5所示。
通过计算可以得出如下结论:
1)同样误差条件下,halo轨道控制方式比拟halo方式所需的速度增量更多。
2)同样误差条件下,半圈一次的控制频率比一圈一次所需的速度增量大幅下降。
考虑到推进剂约束,在保证halo轨道形状的同时尽量减小轨控所需的速度增量,中继星轨道维持选择拟halo轨道的控制方式,通过调整控制时刻速度三方向分量,使控后第三次过XOZ面时X方向的速度为零,控制频率为半圈(约7 d)一次。
2.2 转移轨道对于地球至L2点的转移轨道设计,需要考虑各类工程要求和约束来确定轨道参数,其中速度增量需求是参数选择中一项重要的约束条件。中继星采用月球引力辅助变轨的转移方式,由运载火箭发射直接进入地月转移轨道,到达近月点时进行近月制动,减速的速度增量不足以捕获成月球卫星,中继星进入稳定流形,飞行约3~4 d到达L2附近,在到达L2后的第一圈轨道上中继星将执行2~3次的捕获机动,进入halo轨道。在进入使命轨道前的飞行过程中,中继星总速度增量为近月制动和捕获控制速度增量之和。
针对近月制动和捕获的速度增量的影响因素开展了研究。近月制动速度增量受发射日期、转移时间、近月点高度和捕获轨道倾角等参数的影响,捕获控制的速度增量受近月点状态和目标轨道振幅等参数影响。本节通过采用单变量分析方法,对比不同参数对轨道特性和速度增量的影响,分析讨论了满足探测任务需求条件下的转移轨道参数选择问题。
2.2.1 近月制动1)发射窗口。计算2个月内不同发射窗口下地月转移轨道和地月L2点转移轨道近月制动所需的速度增量,结果见图 7(a)。图中下方实线表示捕获成环月圆轨道所需的速度增量(802~852 m/s)(右纵坐标),上方实线表示形成地月L2平动点轨道所需要的近月制动速度增量(160~216 m/s),两者变化趋势一致。进入地月L2点轨道所需的近月制动速度增量最小和最大相差约60 m/s。
由于发射窗口的选择受到运载近地点幅角的约束,用虚线给出了出发近地点幅角的变化曲线。可以看出,2018年发射地月L2点探测器时,地月转移轨道近地点幅角最小,即月球位于地球最南纬附近时到达月球的转移轨道,所需的速度增量最小。在确定发射窗口时,需要综合考虑运载和大系统约束条件,适当选择发射窗口来满足任务要求的同时使得到达月球的速度最小。
2)近月点高度。不同近月点高度与近月制动速度的关系曲线如图 7(b)所示,近月点高度与近月制动速度增量基本成单调递增关系,近月点越高,近月制动速度增量越大。当近月点高度为100 km,所需近月制动速度增量200 m/s;近月点高度1 000 km时,近月点制动速度增量比100 km增加约42 m/s。因此,为节省探测器的速度增量,近月点高度的选择不应过高。
3)近月点倾角。近月点倾角影响航天器到达月球时的速度,考虑到探测器到达月球时有升轨和降轨两种飞行方向,图 7(c)给出了不同到达方式下的近月点倾角与近月制动速度的关系。可以看出,近月点轨道倾角越大,近月制动速度增量越小,近月点90°倾角到达月球比较小倾角到达月球,近月制动速度增量减少约15~20 m/s。近月点轨道倾角大于90°后,近月制动一次很难形成地月L2的平动点轨道,所以地月转移轨道的近月点是顺行轨道。
4)地月转移时间。不同转移时间对应的近月点制动速度增量见图 7(d),由图可以看出,当转移时间在4~5 d范围内时,近月制动速度增量达到最小;相比地月转移时间约5 d的轨道,约6 d的转移轨道近月制动速度增量增加不到20 m/s,约3 d的转移轨道近月制动速度增加约50 m/s。
为保证中继星准确进入满足任务要求的使命轨道,在到达L2后的第一圈轨道上安排了2~3次变轨,飞行过程示意如图 8所示。探测器飞越月球后,位置“1”“2”和“3”分别表示探测器飞越月球后的3次捕获控制位置。其中第一次捕获用于修正之前的飞行误差,第二次捕获控制通过在位置“2”进行机动,使得位置“3”到达目标位置,随后在位置“3”执行进入halo轨道的捕获控制。
针对不同的近月点参数(倾角、高度)和目标轨道振幅参数对捕获控制速度增量的影响进行了分析,结果如表 6所示。
通过计算可得以下结论。
1)相同倾角、高度到达月球,使命轨道振幅越大,所需要的捕获速度增量越大;
2)相同倾角到达月球,形成相同振幅的使命轨道,近月点高度越高,所需要的速度增量越大,升轨到达比降轨到达所需的速度增量大;
3)对于相同振幅的使命轨道、相同的近月点高度,倾角对速度增量的影响并非单调,对于Az= 1万 km的轨道,最小速度增量出现在最小倾角附近,而对于Az= 1.5万 km的轨道,最小速度增量出现在约降轨10°倾角附近,因此对于不同振幅的轨道,可通过优化倾角来减小速度增量的需求。
3 结 论本文主要针对地月L2点的中继星任务,分析了发射窗口、地月转移时间、近月点高度、近月点倾角、轨道振幅等因素对地月L2转移轨道和L2点使命轨道特性的影响,通过分析得到以下结论:
1)对于转移轨道,采用月球引力辅助变轨转移方式可兼顾转移时间和速度增量,地月转移时间约为4~5 d,到达近月点时为顺行轨道;近月制动速度增量随发射窗口呈周期性变化,且近月点高度越低、倾角越大,近月点制动的速度增量越小。
2)对于使命轨道,中继星采用halo轨道,维持方式采用拟halo方式,控制频率为半圈一次;halo振幅越大对避免月掩和阴影越有利,但中继星对着陆器可见性越差,捕获所需速度增量也越大。
因此在轨道设计时应综合考虑任务需求选择合适的轨道参数。本文通过各项参数的计算分析,明确了对轨道特性的影响,可为地月L2点轨道参数的选择、设计和优化提供参考。
[1] | Barrow-Green J. Poincaré and the three-body problem [M]. Providence, Rhode Island: American Mathematical Society-London Mathematical Society, 1997. |
[2] | Farquhar R W. The Flight of ISEE-3/ICE: origins, mission history, and a legacy [C]//AIAA/AAS Astrodynamics Specialist Conference and Exhibit. Boston, Massachusetts: AIAA, 1998. |
[3] | Farqubar R W. The Utilization of halo orbits in advanced lunar operations [R]. Greenbelt, Maryland: Goddard Space Flight Center, 1971. |
[4] | Miller J K. Targeting an optimal lunar transfer trajectory using ballistic capture [R]. Pasadena, California: Jet Propulsion Laboratory, California Institute of Technology, 1994. |
[5] | Gordon D P. Transfers to Earth-Moon L2 Halo orbits using lunar proximity and invariant manifolds [D]. Indiana: Purdue University, 2008.http://www.tdx.cat/bitstream/handle/10803/2121/EMA_THESIS.pdf.txt?sequence=2 |
[6] | Robin Biesbroek. Study on lunar trajectory from GTO by using weak stability boundary transfers and swing-by’s [R].[S.l.]: ESTEC Working Paper, 2014. |
[7] | Belbruno E A, Carrico J P. Calculation of weak stability boundary ballistic lunar transfer trajectories [C]//AIAA/AAS Astrodynamics Specialist Conference. Denver, Colorado, USA:AIAA, 2000. |
[8] | McInnes A. An Introduction to Libration Point Orbits [EB/OL]. (2009)[2017]. http://www.ibrarian.net/navon/paper/An_Introduction_To_ Libration_Point_Orbits.pdf?paperid=4651888. |
[9] | Martin W L. Libration point trajectory design[J]. Numerical Algorithms, 1997, 14: 153-164.doi: 10.1023/A:1019108929089 |
[10] | Szebehely V. Theory of orbits: the restricted problem of three bodies [M]. New York: Academic Press, 1967. |
[11] | 杨维廉. 发射极月卫星的转移轨道研究[J].航天器工程, 1997, 6(3): 19-33. |
[12] | 周文艳, 黄昊, 刘德成, 等. 嫦娥二号卫星日地L2点扩展任务轨道设计[J].中国科学: 技术科学, 2013, 43(6): 609-613. |
[13] | Dunham D W, Jen S J, Roberts C E, et al. Transfer trajectory design for the SOHO libration-point mission[C]//43rd International Astronautical Congress. Washington, USA: [s. n.], 1992. |
[14] | Folta D C, Pavlak T A, Howell K C, et al. Station keeping of lissajous trajectories in the earth-moon system with applications to ARTEMIS [C]//20TH AAS/AIAA Space Flight Mechanics Meeting. San Diego, CA: AAS/AIAA, 2010. |
[15] | Howell K C, Pernicka H J. Stationkeeping method for libration point trajectories[J].Journal of Guidance, Control, and Dynamics, 2013, 43(6): 609-613. |