直接再入大气的月地转移轨道设计

本立言，严玲玲，谢祥华，张锐，*，王国际

（1.中国科学院微小卫星创新研究院，上海201203； 2.上海微小卫星工程中心，上海201203）

自中国探月工程开展以来，在地月转移轨道设计、中途修正等方面取得了丰富的研究成果，相关技术已成功应用于工程实践［1-3］。根据探月工程总体规划，中国将继续开展月球无人采样返回任务，轨道设计及轨道控制面临着新的挑战和难点。月地转移并不是地月转移简单的逆过程，相关约束条件更加复杂。根据月球采样返回任务的需求，从月球停泊轨道出发直接再入大气的月地转移轨道是月球采样返回任务的首选方案。

月地转移轨道设计不同于地月转移轨道。首先通过调整发射窗口，可容易获取共面的地月转移轨道，但由于月球停泊轨道面的固定，月地转移的初始状态不一定满足共面转移的条件，则需要进行轨道面的调整，并且探测器再入大气返回地球时，受回收场的约束，需要匹配落点的地理位置，避免轨道面的调整，保证任务可靠实施。

目前，国内外针对月地转移轨道设计展开了深入的研究。Ocampo和Saudemont［4］研究了利用单脉冲和三脉冲机动方式从环月轨道出发进入月地转移轨道的问题，但未考虑再入点约束。高玉东等［5］基于双二体模型假设，利用二分法搜索从月球停泊轨道到地球停泊轨道的共面月地转移轨道，但未涉及直接再入大气的约束。张磊等［6-7］以近月点速度、月球停泊轨道升交点赤经及近月点幅角为设计参数，以再入点地心距、再入角和再入倾角为终端参数，通过两级修正设计满足约束的月地转移轨道，在此基础上通过分别调整再入倾角和飞行时间完成对落点地理位置匹配，求解过程复杂。郑爱武等［8-10］以探测器出月球影响球的时刻，位置和速度为中间变量，将月地转移轨道分为地心段和月心段，通过调整影响球边界处的位置和速度获得一条连续的月地转移轨道，在此基础上给出了精确月地转移轨道设计方法和月地返回窗口搜索策略，但初值选择依赖于设计者的经验且未涉及月球停泊轨道面调整。汪中生等［11］分析了多种月地转移方案，并给出了相应的基本设计方法，但其求解效率尚有改进和优化的空间。

由于三体问题不存在解析解，月地转移轨道设计必须采用数值方法。一般的思路是：根据边值约束，采用简化模型设计初值，之后考虑精确摄动模型，采用微分修正等数值方法求解精确转移轨道。在目前的研究成果中，伪状态理论［12］是一种计算三体问题的近似解析形式，由于考虑地月球中心引力共同作用，预报精度相对较高，广泛应用于三体轨道的初值设计。

针对以上研究现状，本文基于伪状态理论，以月球停泊轨道为出发点，考虑再入点的边界约束条件，通过简单迭代设计单脉冲月地转移轨道的初值，在此基础上利用微分修正方法求解精确转移轨道，并对月地返回窗口进行搜索和分析。

1 动力学模型

在地球惯性系下，建立探测器动力学模型，模型中考虑了地球非球形项、月球及太阳的影响，其中，地球非球形引力取8×8阶，月球和太阳星历均采用DE405。

式中：r和v分别为探测器的位置和速度矢量；μe、μm和μh分别为地球、月球和太阳的引力常数；rh和rm分别为太阳和月球的位置；rhd和rmd分别为探测器到太阳和月球的位置；R为地球引力摄动位函数。

2 设计约束

由于探测器运行在月球停泊轨道上，机动时刻探测器的位置无法改变，只能将速度和转移时间作为自由变量，故初始轨道参数有速度和转移时间。则初始设计参数选取为对于直接再入大气的月地转移轨道，通常对再入时刻的地心距、再入角、再入倾角有严格要求，为匹配落点的地理位置，再入轨道面与回收场必须共面，可采用回收场位置与再入轨道角动量间夹角来描述，则目标终端参数选取为

式中：Re为再入时刻的地心距；γe为再入时刻的再入角；ie为再入时刻的再入倾角；φe为再入时刻回收场位置与再入轨道角动量的夹角。

3 伪状态理论

伪状态理论是W ilson［12］提出的一种考虑地月中心引力作用下轨道近似计算方法。

考虑地球、月球及探测器三体系统，认为探测器对地球和月球运动没有影响。给定探测器tI时刻相对地心的位置速度ΣI=（RI，VI），利用伪状态理论可以计算tk时刻相对地心的位置速度Σk=（Rk，Vk）和相对月心的位置速度σk=（rk，vk）。定义以tk时刻月球位置为中心的伪状态转换球（Pseudostate Transformation Sphere，PTS），半径为RPTS。假设RI在PTS外部，Rk在PTS内部，如图1所示。

图1 伪状态理论示意图Fig.1 Schemetic of pseudostate theory

具体计算步骤如下：

3）月心二体段。根据tc时刻相对月心的位置速度σxc，按月心二体轨道求解tk时刻相对月心的位置速度σk，即可得到相对地心的位置速度Σk。

4 初值设计

本节给出了月地转移轨道的初值设计方法。若机动时刻探测器的位置在近月点前，选择近月点对应的伪状态位置作为迭代变量，转移轨道分为两段处理，转移轨道在近月点前按月球二体轨道处理。若机动时刻探测器的位置在近月点后，选择机动时刻探测器对应的伪状态位置作为迭代变量。

定义机动时刻t0探测器在月球惯性系下位置为r0，近月点前转移时间Δt1，近月点后转移时间Δt2，伪状态位置为Rs。

步骤1计算t0时刻月球位置速度Rm和Vm，令Rs=Rm，Δt1=0。

步骤2定义回收场的经度为λl，纬度为φl，位置为Rl，he为地心段角动量，ke为地轴的单位方向。根据余弦定理，可得

式中：α的符号对应了升降轨方式，β的符号对应了轨道面的方位，决定了再入航程大小（见图2）。

图2 月地转移轨道的几何约束Fig.2 Geometric constraints of Moon-to-Earth transfer orbit

定义t0时刻Rs在地球固连系下的经度为Gs，再入时刻Rl在地球固连系下的经度为Gl=Gs+α+β，记G=Gl-λl，则地心段的飞行时间为

式中：ωe为地球自转角速度；N为指定的天数。

则Rs为返回轨道地球段的起始点，探测器到达再入点时地心距为Re，再入角为γe，采用Luo等［13］提出的飞行路径角约束的Lambert算法可以求解该类边值问题，这里不再赘述。

再入点速度为

令Rs处的飞行角为γs，速度大小为Vs，Vs由式（9）计算：

步骤3令探测器在PTS边界处相对月球的速度为vs，位置为rs，其大小为PTS半径，在月球无穷远处的速度大小为

轨道半长轴为

定义Δf为初始位置r0和v∞间的夹角，由几何关系可知：

式中：f∞为月球无穷远处真近点角；f0为机动时刻的真近点角。

令em为偏心率，存在以下关系：

为完成月地转移，机动时刻探测器速度应为

若f0≥0，根据开普勒飞行时间公式计算r0到rs的转移时间ts，否则转向步骤5。

步骤5此时机动时刻探测器位置在近月点前，根据开普勒飞行时间公式分别计算r0到rp的转移时间Δt1和rp到rs的转移时间ts，同时更新Rm为近月点时刻的月球位置。

步骤6更新Rs。根据伪状态理论得到

5 精确解求解

初始设计参数P和目标终端参数Q存在确定的函数关系为

则Q和P的误差之间的关系可线性近似为

为提高求解算法的速度，将轨道积分停止条件设置为到达给定再入角γe，则初始设计参数P中剔除了转移时间Δt，目标终端参数Q中解除了再入角γe的约束，则月地转移轨道设计降价为一个三维的两点边值问题。

令X=［r，v］T，则式（23）可以重写为

当积分停止条件设置为到达给定再入角，则存在以下关系：

式（27）给出了消除时间变量后偏差传递矩阵形式：

式中：除矩阵∂Xf/∂Xi外均有确定的形式，需要注意的是，对目标终端参数cosφe偏导数计算时，需要考虑由于再入时间变化引起的返回场惯性系下位置变化。由于三体轨道动力学特性导致矩阵∂Xf/∂Xi没有解析表达式，本文采用文献［2-3，14-15］提出的方法求解。

得到误差传递矩阵后，可利用各种迭代算法搜索精确解，本文采用微分修正方法求解精确解。

6 算例与分析

本节首先给出了一个月地转移轨道设计算例。相关算法均采用C++编程，程序运行环境：CPU为Intel Core 2.53G。

假设探测器在2022年1月1日进行单脉冲机动，此时探测器在月球惯性系下位置为［1 937.4，0，0］km，转移时间为3～4 d，再入点地心高度为120 km，再入倾角为45°，再入角为-6°，返回场的经度为110°，纬度为40°。

根据第4节的初值设计方法，存在4条月地转移轨道满足给定约束条件，下面以其中一条为例，表1给出了伪状态位置误差的3个分量迭代过程。在初值设计时，迭代6次，耗时1 ms，在精确值求解时，迭代4次，耗时3 s，即可满足设定精度要求，表2给出了再入点状态偏差的迭代过程。比较初值［-337.339，2 218.060，914.856］m/s与精确解［-346.252，2 219.854，914.926］m/s，两者模值的差为2.925m/s，两者矢量差的模值为9.092m/s，为初始速度大小的0.37%，在任务分析阶段，可用初值代替精确解，为任务分析提供参考。

表1 初值设计时伪状态位置误差迭代过程Table 1 Iteration process of pseudostate position error during initial solution design

表2 精确解求解时再入点状态偏差迭代过程Table 2 Iteration p rocess of reentry point state error during exact solution search

表3给出满足给定约束条件的4条月地转移轨道的关键特征参数，图3和图4给出相应的3D和2D示意图。分析结果可知，4条月地转移轨道通过调整转移时间满足了再入轨道与返回场的共面要求，再入点速度大小基本相同，再入航程因轨道方位不同而出现差别。

下面计算2022年1月1日到2022年1月6日从月球停泊轨道出发的月地返回窗口，其中，月球停泊轨道在2022年1月1日的具体参数：高度为200 km，偏心率为0，倾角为22.4°，升交点赤经为0°，纬度幅角为0°，目标轨道参数与上述算例相同。其中2022年1月1日满足共面转移条件，这是月地转移能量最省方式之一，以此为基准分析速度增量变化趋势。由于月球停泊轨道周期较小，角速度较快，速度方向变化较大，在每个轨道周期内变轨速度增量变化剧烈，如图5所示，在每个轨道周期内，存在一纬度幅角使变轨速度增量最小。

表3 月地转移轨道特征参数Table 3 Characteristic param eters of M oon-to-Earth transfer orbit

图3 月地转移轨道3D示意图Fig.3 Three-dimensional illustration of Moon-to-Earth transfer orbit

图4 月地转移轨道2D示意图Fig.4 Two-dimensional illustration of Moon-to-Earth transfer orbit

在窗口搜索时选取每个轨道周期内最小变轨速度增量，分析其变化趋势。由图6可知，结果有明显的周期性。在每一天内变轨速度增量逐渐增大，变轨速度增量差别最大为40m/s。随着探测器出发时刻延迟，月球停泊轨道与月地转移轨道面夹角增大，如图7所示。若实现月地转移，需要调整轨道平面，变轨速度增量逐渐增大，再入航程同时逐渐增大，若月球出发时间延迟6 d，变轨速度增量增加130m/s。

图5 一个轨道周期内变轨速度增量变化趋势Fig.5 Track change velocity increment versus time in one orbit period

图6 变轨速度增量变化趋势Fig.6 Track change velocity increment versus time

图7 轨道平面改变量变化趋势Fig.7 Change of orbit plane versus time

通过对月地返回窗口搜索可知，单脉冲转移方式通过调整轨道面，增加了月地返回机会，保证任务顺利进行。但随着对月球探测的继续深入，月球探测任务更加复杂，对探测器中止任务能力的需求相应增加，显然单脉冲转移方式在轨道平面调整量较大时，需要消耗更多能量，无法满足任务的要求，因此对多脉冲月地转移轨道的研究是后续工作的重点。

7 结 论

本文给出了从月球停泊轨道出发直接再入大气的月地转移轨道设计方法。

1）初值设计算法精度高，耗时少，在任务分析阶段可作为代替精确解，为任务分析提供参考依据。精确解求解算法求解效率高，可以作为设计月地转移轨道及搜索月地返回窗口的一个有力工具。

2）给出了满足给定约束的4条月地转移轨道，并对比了关键特征参数，为月地转移轨道的选择提供依据。

3）通过月地返回窗口分析可知，由于匹配返回场，每天内变轨速度增量逐渐增大。当月球停泊轨道与月地转移轨道面夹角逐渐增大时，使用单脉冲机动，则轨道面调整量增大，导致燃料消耗增大，限制了月地转移机会。