基于不变流形的地月L2点Halo轨道转移轨道设计

彭 坤，李明涛，王 平，田 林，果琳丽，杨 雷

（1.中国空间技术研究院载人航天总体部，北京100094；2.中国科学院空间科学与应用研究中心，北京100190）

基于不变流形的地月L2点Halo轨道转移轨道设计

彭 坤1，李明涛2，王 平1，田 林1，果琳丽1，杨 雷1

（1.中国空间技术研究院载人航天总体部，北京100094；2.中国科学院空间科学与应用研究中心，北京100190）

针对Halo轨道转移轨道设计不易收敛的问题，结合月球探测背景，分析了地月L2点Halo轨道及其不变流形可靠近月球的运动轨迹特性，给出基于二分法的先粗选后精选的零消耗转移轨道搜索方法并研究得出Halo轨道全相位点入轨的零消耗转移轨道不超过2条，提出一种最小x轴约束的近月点终止条件和自适应退步搜索的改进微分修正算法，对Halo轨道全相位点入轨的转移轨道设计进行求解。仿真结果表明，该改进微分修正算法收敛速度快，能有效避免奇异，且适应性强，能搜索出全相位点入轨的所有转移轨道。

月球探测；平动点；Halo轨道；不变流形；转移轨道

1 引言

近年来，随着月球探测热潮的兴起，地月系统L2点（简称EML2点）附近的Halo轨道逐渐受到航天界的重视。早在阿波罗任务时期，Farquhar就提出在EML2点附近Halo轨道放置一颗卫星，为阿波罗月球背面探测任务提供测控通信支持［1］，由于资金原因最终计划取消，但该方案为月球背面探测任务的测控通信方案提供了一种思路；也有学者提出可以将登月航天器停泊在EML2点附近Halo轨道，航天员在航天器上遥控着陆月面的机器人进行探测［2］。对于后者，往往需要两个航天器在Halo轨道进行交会，因此需研究月球到Halo轨道全相位位置的转移轨道特性。

在平动点Halo轨道转移轨道研究方面，How⁃ell等利用不变流形和微分修正算法求解了最佳Halo轨道入轨点的日地系统L1点Halo轨道与地球之间的转移轨道［3］。Zazzera等［4］人基于不变流形和庞加莱截面设计了无约束转移轨道。胡少春等［5］将序优化理论与微分修正相结合求解日地系统L1点Halo轨道到地球的转移轨道，增加了微分修正的收敛范围。李明涛［6⁃8］系统研究了日地系统和地月系统平动点Halo轨道转移轨道设计，推导了各轨道约束的微分修正公式，求解了多约束条件下的精确转移轨道。乔栋等［9］以嫦娥二号扩展任务为背景设计了环月轨道到日地L2点的转移轨道。张景瑞和曾豪［10⁃11］在李明涛的基础上，设计了求解日地系统Halo轨道的多约束转移轨道分层搜索算法以及不同月球借力条件的地月系统Halo轨道转移轨道设计。

虽然上述文献对平动点Halo轨道转移轨道搜索算法进行了多方面研究，但对于Halo轨道到次天体的转移轨道搜索算法并未进行详细讨论，且大多研究为Halo轨道上单个入轨点的转移轨道，并非针对Halo轨道全相位点进行的转移轨道设计。本文在圆型限制性三体模型下，以月球到地月系统L2点Halo轨道的转移轨道设计为例，提出一种基于不变流形的终止条件快速识别和自适应退步搜索的平动点转移轨道改进微分修正算法，减少转移轨道搜索迭代次数和轨道仿真时间，同时提高转移轨道收敛性，并分析了月球到Halo轨道不同相位的转移轨道特性，为未来实施基于平动点的月球探测计划提供参考。

2 EML2点Halo轨道及其不变流形设计分析

2.1 圆型限制性三体模型

对Halo轨道转移轨道的轨道特性初步分析时，无需采用复杂的高精度星历模型［2］，本文采用圆型限制性三体模型来设计Halo轨道及其转移轨道。令地球m1和月球m2为主天体，假设m1和m2绕其地月系统质心做角速度为ω的匀速圆周运动，且航天器质量m3≪m2＜m1。以地月系统质心为原点，地月连线由地球指向月球方向为x轴，z轴沿角动量方向，建立会合坐标系O－xyz。利用拉格朗日方程可推导出会合坐标系下航天器的动力学方程，并进行归一化处理可得式（1）［12］：

式中变量定义如式（2）：

2.2 Halo轨道设计与分析

将动力学方程从会合坐标系下转化到平动点L2坐标系下，并除以距离尺度 γ2，即 ρ＝（ξ，η，ζ）T＝1/γ2·（x－x2，y，z）T，则可将其平动点附近运动构造为式（3）形式：

其中Ax和Az分别表示平面内、外的振幅，λ和υ分别表示平面内、外的频率，φ和ψ为相位角。当λ和υ相等时则形成Halo轨道。采用Richardson三阶近似值作为Halo轨道初值［13］。

式中Φij为状态转移矩阵Φ分量。根据以上方法可求解出EML2点附近的Halo轨道在归一化单位下的飞行轨迹，如图1和图2。图1为Az＝8000 km的Halo轨道三视图，图2为Az＝2000～32000 km的南北向Halo轨道簇。由图2可知，Halo轨道按照其z轴偏向，可分为南、北向两簇。

图1 EML2点单个Halo轨道轨迹Fig.1 Single Halo orbit near EML2 point

图2 EML2点南北向Halo轨道簇Fig.2 North and south Halo orbit cluster near EML2 point

2.3 Halo轨道不变流形设计与分析

不变流形是与Halo轨道光滑连接的一簇空间轨道，分为稳定流形和不稳定流形。其中稳定流形逐渐趋近Halo轨道，不稳定流形逐渐远离Halo轨道［12］。因此，可以将Halo轨道的不变流形作为Halo轨道转移轨道的初值，在其基础上进行调整可得到Halo轨道转移轨道。

不变流形计算过程主要有三步［12］：1）计算Halo轨道某点状态及其积分一个周期后状态所对应的状态转移矩阵Φ（0，T）（也称单值矩阵）；2）根据不变流形的类型和方向计算该点Φ（0，T）的特征向量、扰动量，从而求出该点不变流形初始状态；3）利用动力学方程进行积分，得到不变流形的轨迹。

图3给出Az＝8000 km的Halo轨道的不稳定流形和稳定流形在归一化单位下的轨迹，扰动值ε＝50 km，积分时间ΔT＝18 d。由图3可知，EML2点Halo轨道的左向流形逐渐靠近月球，可以作为转移轨道到达环月轨道；而其右向流形则延伸到远离月球的方向，无法作为转移轨道到达环月轨道。

图3 EML2点Halo轨道的不变流形Fig.3 Invariant manifolds of Halo orbit near EML2 point

3 基于不变流形的零消耗转移轨道分析

由2.3节指出的不变流形轨迹特性可知，利用左向稳定流形可以设计从环月轨道到L2点Halo轨道的转移轨道。若不变流形的最小月心距等于环月轨道的月心距，则该不变流形可直接作为转移轨道，称为零消耗转移轨道（入轨点处无需变轨）［8］。以下以Az＝8000 km的Halo轨道为例，设环月轨道高度为HL＝100 km，即近月点月心距为RLP＝1837.400 km，分析其零消耗转移轨道特点。

1）将Halo轨道按时间等分为360段，初始点为1，与重点361重合，如图4所示。

图4 EML2点Halo轨道分段Fig.4 Different segments of Halo orbit near EML2 point

2）生成Halo轨道上每个节点的左向稳定流形，积分时间取为ΔT＝18 d，如图5所示。计算每条不变流形距离月球的最小月心距，看是否有流形最小月心距等于环月轨道月心距。由图5可得，整圈Halo轨道有两处（约为86节点和186节点）不变流形的最小月心距与环月轨道月心距相等，可以作为零消耗转移轨道。

3）按照文献［8］的思路，可采用二分法求解精确的Halo轨道不变流形节点。本文采用先粗选，再精选的流程搜索精确的Halo轨道不变流形节点。首先，选取Halo轨道上一定间隔的采样点，计算其不变流形最小月心距与环月轨道月心距的差值，记录月心距差值突变的采样点数。本文间隔10个节点采样，既保证采样精度，同时减少计算量。其次，以差值突变处的前后两个采样点作为边界，设置月心距误差精度，采用二分法精确搜索满足月心距约束条件的节点。

4）根据以上算法，可以精确求解出2条满足月心距约束的零消耗转移轨道，其结果如图6～7和表1所示。由搜索结果可知，若利用天然不变流形设计零消耗的转移轨道，则满足近月点月心距约束的转移轨道个数较少，无法满足Halo轨道全相位（1～361节点）的轨道转移需求。

图5 EML2点Halo轨道左向稳定流形及其最小月心距Fig.5 Left stable invariant manifolds of Halo orbit near EML2 point and its minimum selenocentric dis⁃tance

表1 零消耗转移轨道搜索结果Table 1 Search results of zero consumption transfer trajectory

4 基于不变流形的两脉冲转移轨道改进微分修正设计

4.1 问题描述

为满足Halo轨道上全相位轨道转移的要求，需要在天然不变流形的基础上，在Halo轨道入轨点加入脉冲变轨，使稳定流形的近月点月心距满足约束，此时形成的流形称为受摄流形［12］。同时在近月点施加脉冲，使航天器从环月轨道进入受摄流形。通过这两次脉冲变轨，可得到满足Halo轨道全相位点入轨的转移轨道［8］。由于Halo轨道处于三体混沌区内，其动力学特性具有强非线性，转移轨道搜索容易发散，本文考虑以其不变流形为初值，进行微小速度增量调整，得到入轨速度增量较小的两脉冲转移轨道。

图6 零消耗转移轨道精确入轨点Fig.6 Precise injection point of zero consumption transfer trajectory

图7 零消耗转移轨道轨迹Fig.7 Zero consumption transfer trajectory

4.2 微分修正算法

本文重点研究转移轨道特性及其设计方法，故主要考虑近月点月心距约束a1＝和航迹角约束 a4＝sinγ－sinγdes，其中RfT＝Rf－RT为相对月球中心的位置矢量，Rf为航天器位置矢量，RT为月球位置矢量，Rdes为月心距目标值，γ为航天器航迹角，γdes＝0°为航迹角目标值。为减少修正约束量，采用航迹角为终止条件，则入轨点速度增量的微分修正算法如式（5）［7］：

通过求导链式法则可得式（6）：

令状态转移矩阵为Φ，则可得终端位置速度与初始位置速度的转换关系，如式（7）所示：

由式（7）可知终端位置速度与初始速度的转换关系分别对应Φ的一部分，如式（8）：

4.3 快速终止条件设计

在微分修正的轨道积分过程中，需要判断a4＝sinγ－sinγdes＝0。若直接判断a4＝0，则有时会收敛到远月点；同时对于Halo轨道上n＝1或181的节点，其起点位置即满足a4＝0，会导致奇异。为此，有些学者对航迹角终止条件进行处理，以增加算法收敛性。文献［11］中设置固定积分时间ΔT，找出所有a4＝0的点，并比较其月心距，以最小月心距作为其积分终止点，该方法大大改善了算法收敛性。但该算法每次需要积分固定时间ΔT，且需要比较月心距，使每次轨道积分时间较长；同时该算法也会将远月点作为积分终止点，造成修正矩阵不准。

本文考虑近月点的特征，加入航迹角导数，快速确定近月点，节省轨道积分时间；同时加入最小x轴约束，避免n＝1或181节点处产生奇异以及收敛到月心距较大的局部近月点。其终止条件判断公式如式（9）：

4.4 自适应退步搜索

由于Halo轨道周围相空间的强非线性，即使以不变流形为初值，在微分修正过程中也会出现终止条件不满足a4＝0，从而给出错误的修正矩阵信息，导致转移轨道发散或收敛到入轨速度增量较大的转移轨道上。为此，本文提出一种自适应退步搜索。当采用微分修正算法给出的修正值时，若轨道积分停止到固定时间ΔT且月心距较大，则视为迭代错误返回上一步速度增量值，并将修正值减半后再进行积分，如此不断循环直至迭代跳出错误区域。其具体算法如式（10）所示：

式中，RfT＝‖RfT‖，ξ（2）为L2点距离月球的无量纲化长度。

4.5 设计流程

总结以上策略，可得到两脉冲转移轨道设计流程，如图8所示。

图8 两脉冲转移轨道改进微分修正算法Fig.8 Improved differential correction algorithm for designing two⁃impulse transfer trajectory

5 两脉冲转移轨道仿真验证

5.1 仿真参数

考虑对月球的覆盖性要求，Halo轨道设为地月L2点的北向Halo轨道，其z轴幅值Az＝8000 km［2］。设转移轨道的近月点月面高度100 km，即目标月心距为Rdes＝1837.400 km，目标航迹角为γdes＝0°。由2.3节分析可知，不变流形积分18 d左右就可以到达月球附近，故设转移轨道最长飞行时间ΔT＝20 d。

5.2 快速终止条件验证仿真

取Halo轨道上节点n＝1为入轨点，仿真比较本文快速终止条件（终止条件1）与γ＝0°的最小月心距终止条件（终止条件2）的计算速度和结果，如图9所示。本文快速终止条件一次轨道积分所用仿真时间为t1＝0.79 s，γ＝0°的最小月心距终止条件的仿真时间为t2＝0.98 s，前者每次轨道积分比后者快0.19 s，说明本文终止条件轨道仿真速度更快。

图9 两种终止条件的转移轨道飞行轨迹Fig.9 The transfer trajectories with two terminal conditions

5.3 自适应退步搜索验证仿真

取Halo轨道上节点n＝136为入轨点，搜索其对应的两脉冲转移轨道并验证自适应退步搜索策略，搜索结果如图10所示。图10中深蓝色轨迹为正常微分修正迭代轨迹，浅绿色轨迹为退步搜索轨迹，品红色粗实线轨迹为最终搜索得到的转移轨道。由图10可明显看出，微分修正迭代到第2步时就出现了迭代错误的情况，此时积分终止时间为ΔT＝20 d，但其航迹角γ≠0°，该状态对应的修正矩阵信息是错误的（图中最上面一条轨迹）。通过自适应退步搜索，经过5步退步搜索后（图中浅绿色轨迹）找到满足γ＝0°的状态，再进行微分修正迭代，最后搜索到满足约束条件的转移轨道（图中品红色轨迹）。

图10 自适应退步搜索策略的迭代轨迹Fig.10 Iteration trajectories of adaptive backstep⁃ping search strategy

5.4 全相位入轨点转移轨道设计仿真

采用改进的微分修正算法，每隔3个节点采样对全相位Halo轨道入轨点设计转移轨道，其仿真结果如图11～13所示。图11绘制了搜索到的全相位Halo轨道入轨点对应的转移轨道。由图11可知，所有的转移轨道都收敛到两条轨道附近，这2条轨道对应零消耗转移轨道。图12和图13分别为Halo轨道上不同节点作为入轨点的转移轨道的转移时间和入轨点速度增量。由图12和图13可知，对于Az＝8000 km的北向Halo轨道，本算法搜索到的以不变流形为基础的转移轨道转移时间范围为9～19 d；入轨点速度增量范围为0～8 m/s，属于小速度增量。

6 结论

1）利用天然不变流形可设计环月轨道到Ha⁃lo轨道的零消耗转移轨道，但其入轨点较少，一般最多有2个；

2）本文提出的利用最小x轴约束的近月点终止条件和自适应退步搜索的改进微分修正算法可快速准确地搜索全相位Halo轨道入轨点对应的两脉冲转移轨道；

图11 全相位Halo轨道入轨点对应的转移轨道Fig.11 Transfer trajectories for all phase injection points of Halo orbit

图12 Halo轨道上不同入轨点转移轨道的转移时间Fig.12 Transfer time of transfer trajectories for all phase injection points of Halo orbit

图13 Halo轨道上不同入轨点转移轨道的入轨速度增量Fig.13 Velocity increment of transfer trajectories for all phase injection points of Halo orbit

3）Halo轨道上不同相位点对应的基于不变流形的两脉冲转移轨道基本收敛到2条零消耗转移轨道附近，入轨点变轨速度增量较小。

（

）

［1］Farquhar R W.Lunar communications with libration⁃point satellites［J］.Journal of Spacecraft and Rockets，1967，4（10）：1383⁃1384.

［2］Hopkins J B，Pratt W，Buxton C，et al.Proposed orbits and trajectories for human missons to the Earth⁃Moon L2 region［C］//64th International Astronautical Congress，Beijing，China，2013.

［3］Howell K C，Barden B T，Lo M W.Application of dynamical systems theory to trajectory design for a libration point mission［J］.Journal of Astronautical Sciences，1997，45（2）：161⁃178.

［4］Zazzera F B，Topputo F，Massari M.Assessment of mission design including utlizatlon of libration points and weak stabili⁃ty boundaries［R］.ESTEC contract No.18147/04/NL/MV，2004

［5］胡少春，孙承启，刘一武.基于序优化理论的晕轨道转移轨道设计［J］.宇航学报，2010，31（3）：662⁃668.Hu S C，Sun C Q，Liu Y W.Transfer trajectory design for Halo orbit based on ordinal optimization theory［J］.Journal of Astronautics，2010，31（3）：662⁃668.（in Chinese）

［6］Li M T，Zheng J H.The optimization of transfer Irajectory for small amplitude Halo orbits［J］.Measurement and Control，2008，41（3）：81⁃84.

［7］李明涛，郑建华，于锡峥，等.约束条件下的Halo轨道转移轨道设计［J］.宇航学报，2009，30（2）：437⁃441.Li M T，Zheng J H，Yu X Z，et al.Transfer trajectory design for Halo orbit with multiple constraints［J］.Journal of Astro⁃nautics，2009，30（2）：437⁃441.（in Chinese）

［8］李明涛.共线平动点任务节能轨道设计与优化［D］.北京：中国科学院空间科学与应用研究中心，2010.Li M T.Low energy trajectory design and optimization for col⁃linear libration points missions［D］.Beijing：Center for Space Science and Applied Research，Chinese Academy of Sci⁃ences，2010.（in Chinese）

［9］Qiao D，Cui P Y，Wang Y M，et al.Design and analysis of an extended mission of CE⁃2：from lunar orbit to Sun⁃Earth L2 region［J］.Adv Space Res，2014，54（10）：2087⁃2093.

［10］张景瑞，曾豪，李明涛.日地Halo轨道的多约束转移轨道分层微分修正设计［J］.宇航学报，2015，36（10）：1114⁃1124.Zhang J R，Zeng H，Li M T.Hierarchical differential correc⁃tion based transfer trajectory design for Halo orbit with multi⁃ple constraints［J］.Journal of Astronautics，2015，36（10）：1114⁃1124.（in Chinese）

［11］张景瑞，曾豪，李明涛.不同月球借力约束下的地月Halo轨道转移轨道设计［J］.宇航学报，2016，37（2）：159⁃168.Zhang J R，Zeng H，Li M T.A design method for Earth⁃Moon Halo orbit transfer trajectory under different constraints to Moon gravity⁃assisted maneuvers［J］.Journal of Astronau⁃tics，2016，37（2）：159⁃168.（in Chinese）

［12］李言俊，张科，吕梅柏，等.利用拉格朗日点的深空探测技术［M］.西安：西北工业大学出版社，2015：21⁃129.Li Y J，Zhang K，Lv M B，et al.Deep Space Exploration Technology by using Lagrange Point［M］.Xi'an：Northwest⁃ern Polytechnical University Press，2015：21⁃129.（in Chi⁃nese）

［13］刘林，侯锡云.深空探测器轨道力学［M］.北京：电子工业出版社，2012：86⁃157.Liu L，Hou X Y.Orbit Mechanics of Deep Space Probe［M］.Beijing：Publishing House of Electronics Industry，2012：86⁃157.（in Chinese）

Transfer Trajectory Design for EML2 Halo Orbit Based on Invariant Manifolds

PENG Kun1，LI Mingtao2，WANG Ping1，TIAN Lin1，GUO Linli1，YANG Lei1
（1.Institute of Manned Space System Engineering，China Academy of Space Technology，Beijing 100094，China；2.Technology and Engineering Center for Space Utilization，Chinese Academy of Sciences，Beijing 100190，China）

To cope with the problem that the transfer trajectory design for Halo orbit is difficult to convergence，and considering the background of lunar exploration，the trajectory characteristic of Halo orbit around Earth⁃Moon L2 point and its invariant manifolds that can reach the moon were ana⁃lyzed in this paper.Then a hierarchical bisection method to design zero consumption transfer trajec⁃tory for the Halo orbit was given and the number of zero consumption transfer trajectory for all phase injection points was no more than two.Finally，an improved differential correction method including the perilune terminal condition with minimum x⁃axis constraint and adaptive backstepping search strategy was proposed to design transfer trajectories for all phase injection points of the Halo orbit.The simulation results indicate that the improved differential correction method has fast convergence speed and can effectively avoid the singularity，moreover it has strong adaptability and can search out transfer trajectories for all phase injection points of the Halo orbit.

lunar exploration；libration point；Halo orbit；invariant manifolds；transfer trajectory

V412.4

A

1674⁃5825（2016）06⁃0673⁃07

2016⁃06⁃15；

2016⁃10⁃14

载人航天预先研究项目（010701）

彭坤（1984－），男，博士，高级工程师，研究方向为航天器总体与轨道设计。E⁃mail：bhkpeng＠126.com