“嫦娥2号”(中国第2颗月球探测器)于2010年10月1日发射，2010年10月9日进入高度为100 km的环月轨道。2011年6月初，“嫦娥2号”在完成了月球全景测绘任务后，择机飞离月球并于2011年8月抵达日地系统L2点附近的拟周期停泊轨道。当时，“嫦娥2号”仍剩余少量燃料可供其后续飞行。Gao等(2012)针对“嫦娥2号”的后续飞行提出了一系列参考方案^[1]，比如：撞击月球或重新捕获到月球轨道、返回地球轨道或大气再入、飞往地月系统L1/L2点Halo轨道或日地系统L1点Halo轨道、飞越近地小行星等。最终，“嫦娥2号”于2012年12月13日近距离飞越了Toutatis小行星(一颗对地球有潜在危险的近地小行星)。在文献[2]中，Gao分析并设计了“嫦娥2号”飞越近地小行星的轨道方案，设计结果表明：“嫦娥2号”可从日地系统L2拟周期停泊轨道直接转移到一条日心拟开普勒轨道，并将在这条轨道上与Toutatis小行星相遇，总共所需消耗燃料的等效速度脉冲约为70 m/s(2012年3月25日出发)或者90 m/s&(2012年7月4日出发)，转移期间不经过地月系统空间。“嫦娥2号”的实际飞行轨道类似于文献[2]所展示的飞行方案。但如果“嫦娥2号”先飞入地月系统，设法利用月球的引力辅助实现轨道机动，然后再飞往Toutatis小行星，探测器所需的总速度脉冲是否可以小于70 m/s是一个值得探讨的问题。

在文献[2]的设计方案中，“嫦娥2号”无需进入地月系统直接从拟周期停泊轨道变轨飞往Toutatis小行星，在这种情况下，月球引力仅被认为是摄动力。利用拟周期停泊轨道的不稳定流形是轨道设计与小行星目标筛选的一种有效方法。拟周期停泊轨道不稳定流形偏离黄道面的幅值大约为4 × 10⁵ km (也可称作为z轴偏&量)，因此，若小行星轨道z轴偏量大于4 × 10⁵ km，则探测器需要进行额外的轨道机动来改变z轴偏量以实现小行星飞越。以飞往Toutatis小行星为例，如图 1所示，探测器在施加了一个小速度脉冲$\Delta {{v}}_{\text{d}}^*$后飞向小行星，中途的速度脉冲$\Delta {{{v}}_{\text{z}}}$的z轴分量大约为60 m/s(针对70 m/s的设计结果)，可见大部分的燃料消耗用于匹配探测器与小行星位置的z轴分量上。为了减少z轴方向上的速度脉冲，本文研究通过借力月球降低燃料消耗的飞行方案。

图 1(Fig. 1)

图 1 速度脉冲改变Δvz轴偏量的示意图[2](飞行轨道表示在以地球为原点的日地旋转坐标系) Fig. 1 Illustration of an orbital maneuver Δvz to change the z-axis bias [2](The trajectory is plotted in the Sun-Earth rotating coordinate frame with the Earth at the origin)

对于穿越月球公转轨道的航天器而言，月球引力辅助很自然地成为提升飞行性能的潜在方案，尤其体现在节省燃料方面。在以往的文献中，已有很多利用月球引力辅助的轨道设计案例，其中大量的工作考虑地月系统内的飞行，比如：建立地月系统共振轨道^[3]、载人探月的自由返回轨道^[4]以及通过近距离飞越月球转移到地月系统L2点的飞行轨道^[5]。利用月球引力辅助飞出地月系统的轨道设计案例包括地球逃逸、俘获行星际人造卫星、基于日地月系统四体问题的轨道机动策略^[6]以及从地球轨道转移到日地系统拉格朗日点。NASA提出的若干拉格朗日点空间任务就利用了月球的引力辅助，其中包括ISEE-3^[7]、Artimis、WIND、MAP和 GAIA 等任务。

显然，“嫦娥2号”拓展飞行需要研究在多体系统下(日-地-月系统)的低能耗轨道设计问题。目前，利用Halo轨道及其不变流形被认为是解决这类问题的有效途径，Halo轨道不变流形代表了一大簇轨道，它们渐近地接近或者远离Halo轨道。在圆形限制性三体轨道模型中，Halo轨道的不变流形形成了一条管道的外表面，不同的流形相互拼接即可应用于轨道设计。Farquhar(1985)^[7]、Belbruno(2004)^[8]、Lo(1998)^[9]、Koon(2000)^[10]、Gomez(2001)^[11]都开展了代表性研究工作。利用Halo轨道不变流形实现轨道转移是将人造卫星从地球停泊轨道转移到地月或日地系统的Halo轨道的经典方法^{[12, 13, 14, 15, 16]}。Howell和Kakoi(2006)还分析了地月系统与日地系统之间的转移^[17]，他们提出了一类适用于稳定流形与不稳定流形之间的拼接方法。“嫦娥2号”首先需要离开日地系统L2拟周期停泊轨道，因此可以考虑合理地利用该轨道的不稳定流形，同时，如何利用月球引力辅助去设计“嫦娥2号”飞越近地小行星的全局最优解也是一个引人入胜的关注点。

本文提出一种方法用于设计低能耗飞越Toutatis小行星的飞行轨道，该轨道始于日地系统L2拟周期停泊轨道并利用一次月球借力，还将与文献[2]的直接转移方案进行对比。本文的轨道设计没有考虑实际的工程约束。通过一次月球借力的小行星飞越轨道可以分为两段：1)飞向月球的改进不稳定不变流形段；2)从月球附近到Toutatis小行星的转移段。我们需要考虑如何将该两段轨道拼接起来。值得注意的是，简单地采用不变流形拼接方法并不能保证得到飞行时间、燃料质量等最优的性能指标，而是需要更加复杂的数值方法实现最优化拼接，同时要考虑多体轨道模型，并通过适当的脉冲机动来实现。我们将侧重描述一种设计逻辑来寻找全局最优的转移轨道，本文所谓的全局最优解可以看成基于本设计逻辑基础上的最优解。

考虑太阳、地球与月球的引力，“嫦娥2号”在地球赤道惯性参考系的运动方程为

其中：探测器的位置和速度矢量分别为r_I = ${[{x_{\text{I}}} \ \; {y_{\text{I}}} \ \; {z_{\text{I}}}]^{\text{T}}}$和v_I = ${[{\dot x_{\text{I}}} \ \; {\dot y_{\text{I}}} \ \; {\dot z_{\text{I}}}]^{\text{T}}}$；地球、太阳和月球的引力常数分别为μ_e、μ_s和μ_m。在地球赤道惯性参考系中，太阳和月球的位置矢量分别用r_s和r_m表示，r_s和r_m可通过JPL DE405星历^[18]计算得到。论文后续所用到的地心惯性坐标系、日地旋转坐标系以及太阳、月球和探测器的位置和速度矢量如图 2所示。相对于利用三体模型及其相互拼接的设计策略，直接在四体轨道模型中设计出的飞行轨道更接近实际飞行。

图 2(Fig. 2)

图 2 地心惯性坐标系和日地旋转坐标系(日地系统质心作为坐标原点)示意图 Fig. 2 Illustration of the geocentric equatorial inertial reference frame and the Sun-Earth rotatingcoordinate framewith the barycenter of the Sun-Earth system as the origin

同时，在日地旋转坐标系中(太阳为较大天体，地球为较小天体，月球引力可忽略)，探测器的运动可以采用圆形限制性三体轨道模型^[19]近似表达

式(2)～式(4)中，探测器的位置和速度矢量分别表示为r = $ {[x\; y\; z]^{\text{T}}}$和v = ${[\dot x\;\dot y\;\dot z]^{\text{T}}}$，μ为地球和太阳的质量比常数(μ = 3.040 36×10^–6)。探测器与太阳和地球的距离分别表示为r₁和r₂。距离、速度和时间均为归一化参数，其单位分别为D(日地平均距离D = 1.495 978 706 6×&10⁸km)、Dn(日地连线在惯性空间的旋转角速度，n = 1.990 990×10^-7 rad/s)和1/n。日地旋转坐标系的xyz轴如图 2所示。式(2)中，U对位置矢量的偏导数为

设定r_I = [x_I y_I z_I]$ {^{\text{T}}}$和v_I = ${[{\dot x_{\text{I}}}\;{\dot y_{\text{I}}}\;{\dot z_{\text{I}}}]^{\text{T}}}\!$的单位分别为km和km/s，r = [x y z]$ {^{\text{T}}}$和v = ${[\dot x\;\dot y\;\dot z]^{\text{T}}}\!$为归一化参数，则在日地旋转坐标系中探测器位置矢量为

式(9)～式(10)中，r_se和v_se分别为地球相对于太阳的位置和速度矢量。日地旋转坐标系中探测器速度矢量v可以简化表示为

相反，r_I = [x_I y_I z_I]$ {^{\text{T}}}$和v_I = ${[{\dot x_{\text{I}}}\;{\dot y_{\text{I}}}\;{\dot z_{\text{I}}}]^{\text{T}}}\!$可根据式(9)～式(11)由r = [x y z]$ {^{\text{T}}}$和v = ${[\dot x\;\dot y\;\dot z]^{\text{T}}}\!$计算得到

式(11)～式(13)成立的关键因素是矩阵Q为正交矩阵。由于文中设定地球相对于太阳的公转轨道为圆轨道，因此矩阵Q必定为正交矩阵。

初始历元时刻t_ref以及相应的“嫦娥2号”探测器位置矢量和速度矢量可参考文献[2]。在t_ref时刻，只需施加一个小速度脉冲即可将探测器维持在拟周期停泊轨道上，且维持时间不小于400天，该速度脉冲通常小于 2 m/s，因此其燃料消耗不计入小行星飞越任务。在t_ref时刻，探测器轨道状态(考虑轨道维持之后)如表 1所示。“嫦娥2号”的拟周期停泊轨道可在四体轨道模型下数值积分得到，在日地旋转坐标系中的飞行轨道如图 3 所示。从图 3可以看出，拟周期停泊轨道至少在两圈之内呈现周期性(轨道周期大约为180天)。值得说明的是，拟周期停泊轨道并非严格的周期轨道。本文设定“嫦娥2号”离开拟周期停泊轨道的时刻从t_ref算起不超过200天。

表 1(Table 1)

表 1 参考历元时刻及其相应的“嫦娥2号”探测器位置与速度矢量 Table 1 Reference epoch and corresponding position and velocity of the spacecraft 历元时刻 位置/km 速度/ (km·s-1) tref = JD 2 456 012.000 766 018 (Mar-25-2012) xI = -1 604 674.532 721 yI = 479 521.542 806 4 zI = 267 254.457 717 4 ${\dot x_{\text{I}}}$ = -18.154 958 19 ${\dot y_{\text{I}}}$ = -31.239 422 35 ${\dot z_{\text{I}}}$ = -142.194 596 64 注：历元时刻为动力学时，JD表示儒略日。

表 1 参考历元时刻及其相应的“嫦娥2号”探测器位置与速度矢量 Table 1 Reference epoch and corresponding position and velocity of the spacecraft

图 3(Fig. 3)

图 3 日地系统L2拟周期停泊轨道，表示在以地球为原点的日地旋转坐标系 Fig. 3 The Change’E-2’s nominal Sun-Earth L2 quasi-periodic parking orbit in the Sun–Earth rotating coordinate frame with the Earth at the origin

经过初步筛选后的小行星目标可以参考文献[2]。最终，Toutatis、2005 NZ6、2002 AW以及2010 CL19选定为最可能在2012下半年与2013年实现飞越的潜在目标。这些潜在目标都是阿波罗型(Apollo-type)小行星，其轨道预报的不确定性均比较小^[20]。飞越小行星的位置设定为地心距2 × 10⁷  km之内且与日地黄道面偏差在±10⁶ km之内。探测器轨道直接转移飞越Toutatis、2005 NZ6、2002 AW，与2010 CL19的速度脉冲分别为90 m/s (Jul-4-2012出发)、55 m/s、367 m/s与37 m/s。本文主要论述飞越Toutatis的转移轨道设计方法，该方法同样可以应用于飞越其他目标。

设定x = [r v]$ {^{\text{T}}}$为圆形限制性三体轨道模型中探测器的轨道状态，其一阶导数$\dot {x} = {f}({x})$参见公式(2)～公式(8)。状态转移矩阵Φ(t,t₀)(t₀为“嫦娥2号”拟周期停泊轨道的初始时刻)为一个6×6的矩阵。Φ(t,t₀)可以反映轨道状态偏差随时间的变化特性，表达式为

状态转移矩阵Φ(t,t₀)的一阶微分方程表达式为

Φ(t,t₀)的初值为Φ(t₀,t₀) = I_6×6，参考轨道状态为x(t) = x^*(t)，雅可比矩阵A(t)的表达为(U的下标表示对相应变量求二阶偏导)

参考轨道状态x(t) = x^*(t)先在四体轨道模型中计算，然后根据式(9)～式(11)转化到日地旋转坐标系。因此，根据公式(15)数值积分x^*(t)可以得到Φ(t,t₀)。设定探测器在拟周期停泊轨道上的始末时刻分别为t₀和t_f，则可以得到任意时刻t_i(t_i = t₀ + a(t_f - t₀)，a∈(0 1)为预先设定的变量)的状态转移矩阵Φ(t_i + T_h,t_i)(T_h为积分时间间隔)。

然后，可以计算Φ(t_i + T_h,t_i)的特征值和特征向量。显然，状态转移的时间间隔T_h不必是一个周期或者一个常值，可取如下经验值

取Φ(t_i + T_h,t_i)特征值中实部绝对值最大者，得到相应的特征向量为[β_r β_v]。对于真实的飞行任务，在t_i时刻施加一个小速度脉冲可使得探测器稳定地离开日地系统L2拟周期停泊轨道。若沿β_v方向施加一个小速度脉冲Δv_d(施加前速度矢量为v^*(t_i))，即可得到探测器离开拟周期停泊轨道的初始速度

由探测器初始位置r^*(t_i)和速度v_un(t_i)(r^*(t_i)和v_un(t_i)应当转化到地心惯性坐标系)在四体轨道模型(参见式1)下数值积分可以得到从t_i至t_un(t_un > t_i)的飞行轨道，即为改进不稳定流形，此定义不要求停泊轨道为严格周期轨道。

总之，改进不稳定流形的设计参数为t_i,Δv_d,s(+1 or -1)和t_un(其中$\Delta {v_{\text{d}}} = \left\| {\Delta {{{v}}_{\text{d}}}} \right\|$)。尽管离开拟周期停泊轨道的速度脉冲的方向不是优化解，但是它不失为确定速度脉冲方向的一种有效方法。值得说明的是，初始状态r^*(t_i)和v_un(t_i)是根据圆形限制性三体轨道模型中关于周期轨道及其稳定性的概念计算得到。

二体轨道模型中的月球引力辅助可简化为位置矢量不变的无燃耗速度脉冲，设定飞越时刻以及飞越前后探测器相对于月球的速度矢量分别表示为t_LGA，${v}_{\infty}^{-}$与${v}_{\infty}^{+}$，${v}_{\infty}^{-}$与${v}_{\infty}^{+}$需要满足如下约束条件

其中：δ为${v}_{\infty}^{-}$和${v}_{\infty}^{+}$的夹角，可由下式计算得到

其中：μ_moon(=4 902.801 km³/s²) 为月球引力常数；R_LGA为近月点的月心距。借力月球可获得的最大速度脉冲由R_LGA的最小值确定

设定$R_{{\text{LGA}}}^{\min }$ = 1 838 km(飞越高度为100 km)，则由月球引力提供的最大速度脉冲多达1.513 km/s，一般而言可以显著地改变飞行轨道。

相对速度矢量${v}_{\infty}^{-}$和${v}_{\infty}^{+}$可由B-平面模型来表达，如图 4所示。设定坐标系$[{\hat {S}}\;{\hat {T}}\;{\hat {R}}]$，其定义如下：$\hat {S}$平行于${v}_{\infty}^{-} $,$\hat {T}=\hat {S} \times \hat {N}$,$\hat {R} = \hat {S}\times \hat {T}$,$\hat {N}={{[0\ 0\ 1]}^{\rm T}}$与月球质心参考坐标系的z轴平行。月球质心参考坐标系可定义如下：z轴由地球指向月球，z轴垂直于月球公转轨道面，y轴满足右手定则，指向月球轨道速度。

图 4(Fig. 4)

图 4 月球引力辅助的几何示意图 Fig. 4 Geometry of lunar gravity assist

因此，${v}_{\infty}^{+}$的表达式为

其中：旋转角度δ可由公式(20)计算得到；θ为任意旋转角度。

然而，上述模型是在二体轨道模型中的近似，为了建立四体轨道模型中月球引力辅助模型，本文建立了一个辅助坐标系$[{\hat {i}} \; \; {\hat {j}} \;\; {\hat {k}}]$，如图 4所示。

首先，${\hat {k}}$为飞越轨道面法方向单位矢量，其表达式为

然后，${\hat {j}}$为${\hat {S}}$轴绕${\hat {k}}$轴顺时针旋转(π - δ)/2得到，${\hat {j}}$满足${\hat {j}} = {\hat {k}} \times {\hat {i}}$。因此，近月点处探测器位置和速度矢量可以在$[{\hat {i}} \; \; {\hat {j}} \;\; {\hat {k}}]$坐标系中表示为

其中：r_p为引力辅助飞越轨道半径，r_p不小于的$R_{{\text{LGA}}}^{\min }$，v_p可由二体轨道模型下轨道能量近似计算得到

综上，月球引力辅助模型从二体近似模型过渡到四体轨道模型。在$[{\hat {i}} \; \; {\hat {j}} \;\; {\hat {k}}]$坐标系中的轨道状态最终转化到地心赤道惯性坐标系。根据近月点处(即月心系下双曲线轨道的近月点)探测器的位置和速度矢量，根据式(1)向前和向后数值积分可得到飞越月球前后的飞行轨道。

总之，月球引力辅助前后决定转移轨道的参数为t_LGA(基于该历元时刻可以得到月球位置矢量r_moon和速度矢量v_moon)、${v}_{\infty}^{-}$(或${v}_{\infty}^{+}$)、r_p和θ。需要注意：${v}_{\infty}^{-}$的概念来源于二体轨道模型中双曲线轨道，可用于计算${v}_{\infty}^{+}$，然后再计算出四体轨道模型中的r_p与v_p。

在这一步骤中，采用改进不稳定流形(沿改进不稳定流形的飞行轨道)检测探测器是否飞越临近月球的区域，这是确定是否能够利用月球引力辅助的必要条件。实际上，“嫦娥2号”从拟周期停泊轨道出发后的转移轨道分为两种方案：1)直接逃离地月系统，飞向小行星；2)先飞向地球，然后逃离地月系统，飞向小行星。前者可参考文献[2]，本文仅考虑后者。

首先，对改进不稳定流形的设计参数进行网格搜索，设计参数网格的设定情况为：t_i = [0:1:200]d，Δv_d = [1:1:10]m/s([x₁:dx:x₂]表示一组离散数值(x₁+ ndx)，n = 0,1,2… 且 (x₁+ndx)≤x₂)。改进不稳定流形的飞行时间(t_un-t_i)设定为300 d。这样，总共需要数值积分计算 2 010条转移轨道。改进不稳定流形上的每条轨道都有可能多次飞临月球，近月点均不难求出。

能够飞临月球的转移轨道在三维空间中可分解成平面内和平面外轨道(平面指的是黄道面)。设定两个约束条件来鉴别探测器是否飞临月球：探测器与月球的平面内最近距离不大于6万 km；平面外z轴方向最近距离不大于20万 km。“嫦娥2号”从拟周期停泊轨道出发且满足以上两个约束条件的转移轨道均被记录。记录的所有转移轨道在飞临月球时刻的探测器和月球的位置，如图 5所示。这些飞行轨道均有可能利用月球引力变轨，不仅可以应用于近地小行星探测任务，也可以应用于其他的飞行任务，如图 6所示。

图 5(Fig. 5)

图 5 月球与探测器在近月点时刻的位置(表示在以地球为原点的日地旋转坐标系) Fig. 5 The Moon’s and spacecraft’s positions at lunar close approaches (expressed in the Sun-Earth rotating coordinate frame with the Earth at the origin)

图 6(Fig. 6)

图 6 两种接力月球后的转移轨道方案(表示在以地球为原点的日地旋转坐标系) Fig. 6 Two types of trajectories with lunar gravity assist in the Sun-Earth rotating coordinate frame with the Earth at the origin

施加合适的中途速度脉冲可以消除近月点时刻平面内和平面外的位置偏差。基于以往轨道设计经验，消除z轴方向位置偏差通常比消除平面内位置偏差需要较多的速度脉冲。因此，在探测器飞离拟周期停泊轨道到飞越月球时间段内需要施加一个经验脉冲Δv₁。不过，Δv₁的取值并无解析解，一种保守估值方法可参考文献[2]：$\left\| {\Delta {{{v}}_1}} \right\| \approx ({v_{{\text{SE2}}}}/{r_{{\text{SE2}}}})\max (\delta z)$，其中r_SE2为太阳到探测器的距离，v_SE2为探测器日心系速度，max(δz)为最大z轴方向偏差。比如说，在z轴方向施加33 m/s的速度脉冲可以引起最大2 × 10⁵ km的位置偏差。因此，设定近月点处z轴方向的位置偏差小于20万 km 的剪枝条件。相反，若消除z轴方向20万 km的位置偏差则需要33 m/s的速度脉冲。此外，消除平面内位置偏差需用的速度脉冲可以暂时忽略。

探测器每次飞临月球时，可以求得四体轨道模型下的轨道状态。同时，探测器相对于月球的位置和速度矢量很容易计算得到。探测器相对于月球的飞行轨道在月球引力范围内可以近似为月心双曲线轨道。月心双曲线轨道在近月点处的速度矢量可通过将近月点位置矢量(${{r}_{\text{p}}} = {r_{\text{p}}}\hat {i}$，如图 4所示)绕月球飞越平面法向量旋转(π-δ)/2(${v}_{\infty}^{-}$由${{v}_{\text{p}}} = {r_{\text{p}}}\hat {i}$和${{v}_{\text{p}}} = {v_{\text{p}}}\hat {j}$求得)得到。相对速度矢量${v}_{\infty}^{-}$用于求解月球引力辅助后的飞行轨道簇。

给定${v}_{\infty}^{-}$值，采用2.3节给出的月球引力辅助模型，r_p和θ是确定探测器后续飞行轨道的设计参数。在每个近月点处，月球引力辅助设计参数网格设置如下：r_p = [1 838:400:50 000] km，θ = [0:0.1:2π] rad。根据3.3节所建立模型，可以得到一簇近月点轨道状态(r_p和v_p)，r_p和v_p与图 5所示状态并不完全一致。基于近月点轨道状态，对探测器运动方程进行数值积分，积分时间段从近月点时刻到飞越小行星最晚时刻。基于文献[2]，飞越小行星最晚时间为JD = 2 456 291，此时Toutatis已飞离近地点，距离地球约为2 × 10⁷ km。月球引力辅助对探测器转移轨道的影响有两种代表性的转移方式，如图 6所示。第一种是探测器重新飞回行星际空间(即逃离地月系统)；第二种是探测器被地月系统俘获。第二种轨道转移方式说明探测器可以转移到地球轨道或者用较少的燃耗实现大气再入。

对于小行星飞越任务，我们只关注逃离地月系统且能够到达小行星附近的转移轨道。对于这些转移轨道，我们主要关注探测器与小行星的相对距离。如果在某时刻相对距离小于图 1所示距离，则月球引力辅助是一种可行的节能方案。我们最终保留了相对距离小于20万 km的转移轨道，并记录了飞越小行星时刻。这样，最终我们发现了两个不同的飞越月球时刻，其对应于不同的月球引力辅助时间窗口如图 7所示。图 7也展示了“嫦娥2号”从拟周期停泊轨道到月球、从月球飞往小行星的轨道段的飞行轨道。所展示的该两段未经优化的转移轨道在近月点处并不连续。

图 7(Fig. 7)

图 7 两组粗略的转移轨道(表示在以地球为原点的日地旋转坐标系) Fig. 7 Two distinct transfer manners (coarse trajectory solutions) in the Sun-Earth rotating coordinate frame with the Earth at the origin

显然，为了消除飞越小行星的位置偏差(包括平面内位置偏差和z轴方向位置偏差)，在探测器转移过程中需要施加一个额外的速度脉冲。在月球引力辅助后探测器转移轨道施加速度脉冲位置及其初值的设计方法参考文献[2]或者2.1节。同样，这种猜测方法也是基于探测器与小行星距离不大于20万 km的剪枝条件。最终，我们得到了粗略的未经优化的探测器转移轨道，包括探测器从拟周期停泊轨道到月球的转移轨道和从月球到小行星的转移轨道，如图 7所示。

按照第3部分的描述，整条飞行轨道可以分为两个主要部分：沿着改进不稳定流形的飞行轨道与借力月球后的飞行轨道。除了在t_d时刻用于离开日地系统L2拟周期停泊轨道的速度脉冲Δv_d，还需考虑另外两个速度脉冲Δv₁与Δv₂，Δv₁在t₁时刻施加，用于拼接改进不稳定流形与借力月球后的飞行轨道，Δv₂在时刻t₂施加，用于消除飞越小行星时刻的位置误差。假设飞越月球时经过近月点的时刻为t_c，小行星飞越时刻为t_flvby，t₁与t₂时刻的约束为

以经过近月点的轨道状态(由r_c,v_c表示)作为边界条件,借力月球前后的飞行轨道可通过对方程(1)的数值积分得到，从t_c到t₂逆向积分得到r_LGA(t₁),v_LGA(t₁)，从t_c到t₂正向积分得到r_LGA(t₂),v_LGA(t₂)。如果改进不稳定流形上一条轨道在t₁时刻的轨道状态r_UMD(t₂),v_UMD(t₂)与r_LGA(t₂),v_LGA(t₂)不相同，则轨道状态在t₁时刻不连续。类似地，如果在t_flvby时刻行星与探测器位置r_ast(t_flvby)与r(t_flvby)不相同，则会存在探测器与小行星之间的位置误差。为了得到可行的飞行轨道，状态不连续性与位置误差应该最终被消除。上述描述可由图 8 展示，揭示了从日地系统L2借力月球飞往Toutatis的基本飞行机制。

图 8(Fig. 8)

图 8 从日地系统L2拟周期轨道借力月球飞越Toutatis的示意图 Fig. 8 An illustration of the transfer from the Sun-Earth L2 quasi-periodic parking orbit to Toutatis via lunar gravity assist

轨道优化问题最终转换为参数优化问题(见表 2，简称为NLP问题)，并利用非线性规划方法求解。在NLP问题中，设计变量包括：t_d,Δv_d,t₁,Δv₁,r_c,v_c,t_c,t₂,Δv₂,t_flvby。r_c,v_c可以表达为${v}_{\infty}^{-}$(or ${v}_{\infty}^{+}$),r_p,与θ的函数，因此可由${v}_{\infty}^{-}$,r_p,θ替代。为了避免撞击月球，需要考虑最小的飞越轨道半径。最终，NLP问题总共包含19个参数，9个等式约束。

表 2(Table 2)

表 2 参数优化问题 Table 2 The parameteroptimization problem 参数 值 优化变量 td, Δvd, t1, Δv1, ${v}_{\infty}^{-}$ , rp, θ, tc, t2, Δv2, tflyby 最小化 $\left\| {\Delta {{{v}}_{\text{d}}}} \right\| + \left\| {\Delta {{{v}}_{\text{1}}}} \right\| + \left\| {\Delta {{{v}}_{\text{2}}}} \right\|$ 等式约束 r(tflyby) = rast(tflyby) rLGA(t1) = rUMD(t1) vLGA(t1) = rUMD(t1) 初始状态 表1与图3给出的日地系统L2拟周期停泊轨道

表 2 参数优化问题 Table 2 The parameteroptimization problem

实际上，图 7所示飞行轨道的粗略解为除了t₁与t₂的所有设计变量提供了一套初值。设计变量t₁与t₂并无先验知识，但由于被限制在固定的时间段内，我们可以为其随机选择一系列初值，从而可以得到一系列局部优化解，再从其中选择一个最好的结果。该结果是在四体轨道模型中优化得到的，如图 9所示，主要的设计参数如表 3所示。在该最优解中，总速度脉冲为58.47 m/s(优化结果表明$\left\| {\Delta {{{v}}_{\text{2}}}} \right\| = 0$)，比直接飞越小行星所需的70 m/s更小。最初认为借力月球可以进一步节省燃料的想法也得以验证。对于图 7中所示的方 案2，我们目前没有得到比70 m/s更优的设计结果。

表 3(Table 3)

表 3 最佳局部最优解的主要设计结果 Table 3 The main solution parameters of the “best” trajectory 设计参数 借力月球 直接转移 [2] 飞越时间 Aug. 29, 2012 N/A 飞越轨道半径/km 6 238.03 N/A 飞越旋转角θ/rad 4.750 1 N/A 月球飞越相对速度/km·s-1 1.185 9 N/A 最小地球距离/km 54 665.721 N/A Toutatis飞越时间 Dec. 13, 2012 Dec. 13, 2012 Toutatis飞越速度/km·s-1 10.83 10.73 总速度脉冲/m·s-1 58.47 70

表 3 最佳局部最优解的主要设计结果 Table 3 The main solution parameters of the “best” trajectory

图 9(Fig. 9)

图 9 从日地系统L2拟周期轨道借力月球飞越Toutatis的最佳局部最优解(太阳-地球旋转坐标系，地球为坐标原点) Fig. 9 The “best” locally optimized transfer trajectory from the Sun-Earth L2 quasi-periodic parking orbit to Toutatis via a single lunar flyby (in the Sun-Earth rotating coordinate frame with the Earth at the origin)

图 10给出了借力月球与直接飞越两组设计结果的飞行轨道图，不难发现借力月球的作用。可以看出，探测器借力月球可以在z轴方向抵达-10⁶ km，探测器无需再施加图 1中所示的Δv_z。对此，优化结果$\left\| {\Delta {{{v}}_{\text{2}}}} \right\| = 0$也是一种解释。当然，飞往月球本身也需要合适的轨道机动，在本飞行方案中，所需的轨道机动Δv₁代价不大，从而最终导致较小的总速度脉冲。

图 10(Fig. 10)

图 10 考虑与不考虑借力月球的Toutatis小行星飞越轨道的对比(图 10(a)即为图 1) Fig. 10 A comparison of the designed Toutatis flyby trajectories with and without lunar flyby (or lunar gravity assist, Fig. 10(a) also shown in Fig. 1)

首先，我们设计了从日地系统L2拟周期停泊轨道出发借力月球飞越Toutatis的低能耗飞行轨道。与不借力月球的飞越方案相比，借力月球有望进一步减小探测器燃料消耗。轨道设计方法也揭示了该飞行方案的基本飞行机制。

同时，我们展示了如何寻找全局最优解的飞行轨道设计逻辑，虽然并未从理论上证明解的全局最优性。改进不稳定流形与借力月球模型以及速度脉冲的经验分配，可用于设计与优化四体模型飞行轨道的有效方法。当然，速度脉冲的经验分配以及速度脉冲的个数，还可以进一步研究。

本研究中，轨道设计直接在四体模型中进行，但是利用了圆形限制性三体与二体模型的先验知识。这也启示我们可以无需事先设计三体模型飞行轨道，然后再利用四体轨道模型进行修正。这个设计策略可以应用于更多的多天体引力场飞行轨道设计问题。

最后，作者向“嫦娥2号”探测任务及其科学与工程团队以及所取得的技术成就^[21]致敬！除本研究之外，已经陆续发表了一系列研究成果，涵盖了轨道设计与测控(文献[22]～文献[25]仅为部分代表性成果，引发了对飞行方案的广泛探讨)、月球探测以及Toutatis小行星探测等诸多领域。按照飞越Toutatis小行星之后的飞行规律，“嫦娥2号”还有望在十几年后再次飞临地月系统^[27]，如果它能够返回地月系统并回到环月轨道，它将超越ISEE-3^[7]，成就一个太空飞行的传奇。