基于共面圆轨道假设的木卫引力辅助捕获轨道设计

孟雅哲，胡海霞，郭建新，韩 冬，孟 斌

（1. 北京控制工程研究所，北京 100190；2. 空间智能控制技术重点实验室，北京 100190）

0 引言

引力辅助（gravity assist, GA)，又称引力甩摆或借力飞行，是指探测器在绕一个大天体飞行时，以相对小的时间尺度经过另一小天体并改变相对大天体轨迹的情况。当大天体为太阳时，小天体通常指行星，如借用金星和地球引力辅助实现木星探测［1］；当大天体为行星时，小天体通常指行星的卫星，如采用月球引力辅助改变地月空间中探测器的轨道［2］及采用伽利略木卫引力辅助改变木星系统中探测器的轨道［3-4］。不考虑经过小天体时所受气动力，引力辅助可用模型包括多体动力学模型［5］、B平面要素模型［6］和引力辅助简化模型［7-8］。其中，最简单的引力辅助简化模型可将复杂的多体问题转化为圆锥曲线拼接问题，从而提升包含引力辅助的轨道的设计效率。

引力辅助轨道的设计包括两大类：引力辅助序列的设计和确定序列下的引力辅助轨道优化。引力辅助序列设计方法主要包括3 类：其一是用物理意义分析轨道的存在性，如Tisserand 图，其应用轨道能量匹配方法确定引力辅助轨道的存在性［7］；其二是采用组合优化方法，如深度广度搜索［9］、进化方法等；其三是构造神经网络模型等智能方法进行计算［7］。确定序列下的引力辅助轨道优化方法包括间接法和直接法。间接法是将引力辅助作为内点约束推导的必要条件，将含有引力辅助的轨迹优化问题转化为微分方程边值问题进行求解［10］；直接法则是将包含引力辅助的轨迹求解问题转化为数值问题，采用基于梯度的方法［11］、粒子群方法等数值优化方法进行求解［12］。

木星捕获轨道为木星探测器轨道的重要组成部分，可通过应用木卫引力辅助来有效节省捕获燃料消耗。文献［13］中假定可以采用脉冲和小推力补充相位差，并在引言中进行了木卫引力辅助捕获的研究情况概述。文献［14］中分析了木卫一、木卫二、木卫三及木卫四各自单次引力辅助捕获所需的剩余速度。Lynam首先给出了从地球出发，引力辅助金星到达木星，并通过连续2个木卫、3个木卫和4个木卫引力辅助捕获卫星轨道的方法［15］，而后对多次连续伽利略木卫引力辅助捕获航天器进行了一系列的研究，包括给定序列的3次引力辅助轨道研究［16-17］；采用拉普拉斯共振和共面圆轨道假设快速搜索各种木卫引力辅助序列下3次和4次引力辅助对应的窗口出现的频率［13］；采用“木卫四→木卫一→近木点→木卫三”引力辅助的捕获轨道［9］。探测器从木星赤道平面进入木卫系统时，应用伽利略木卫引力辅助进行捕获的一般性规律尚无论文进行分析。

为得到木星赤道附近引力辅助捕获轨道的一般性结论，本文假定伽利略木卫运行在共面圆轨道上且航天器轨道和木卫轨道共面，给出一种从给定双曲线轨道出发，按照设定引力辅助次序，通过木卫引力辅助将航天器捕获到绕木轨道的方法。本文首先简要给出了所采用的引力辅助简化模型；然后说明了共面圆轨道假设下的引力辅助时刻可以从木卫相位差反算得到，航天器轨道可通过几何旋转被放置在合适的相位上，并分析了确定共面连续引力辅助轨道的参数；而后给出了求解引力辅助轨道的广度优先搜索方法；最后给出了算例和搜索结果。

1 引力辅助简化模型

本节简要给出所用引力辅助简化模型，即航天器绕天体O 运行时，其轨道被同样绕O 运动的天体G 的引力改变的情况，如图1 所示。航天器沿轨道一运行时，在某个位置近距离飞越天体G，飞越过程中，航天器受G的引力影响，转为沿轨道二运行。假设引力辅助过程在一个时刻点t完成，用t-、r(t-)和v(t-)分别表示引力辅助前的时刻以及该时刻航天器的位置矢量和速度矢量，用t+、r(t+)和v(t+)分别表示引力辅助后的时刻以及该时刻航天器的位置矢量和速度矢量，则在引力辅助时刻前后，航天器的位置与天体G的位置rG(t)重合，见式（1），但航天器的速度从v(t-)变为v(t+)。v(t-)和v(t+)的关系可由航天器和天体G 的相对运动情况确定。引力辅助前后，航天器相对天体G的速度分别为双曲线入射渐近线方向的速度和出射渐近线方向的速度，二者大小相等但方向不同。

图1 引力辅助示意图Fig.1 Diagram of gravity assist

根据圆锥曲线的性质，给定双曲线的入射速度以及双曲线轨道的近地点矢径rp，即可确定双曲线的出射速度。

式中：v-—— 引力辅助前航天器相对于天体O 的速度；vG—— 引力辅助时刻天体G 相对于天体O 的速度；v∞——t时刻航天器与天体G 的相对速度大

小；δ——入射速度和出射速度的夹角；e—— 双曲线轨道的偏心率；rp—— 矢径长度，rp=，其通常受rp≥rG+hmin的限制，其中hmin=50 km，hmin对应的转角被记为δmax；μ—— 天体G 引力常数；φ——自由设定的方向角；i、j和k—— 由v-∞和vG求得的直角坐标系3个轴的方向矢量，见式（6）。

式（3）确定的航天器相对天体G 的出射速度与引力辅助时刻天体G 速度的矢量和，即为引力辅助后航天器相对天体O的速度：

通常来说，我们希望通过设计v-和t，使得v+满足任务要求。任务输入中，可能包含v+的取值，也可能只包含对v+的约束，如要求半长轴变短。通常来说，v-和t的取值应使得航天器和天体在不外加脉冲的情况下尽量可以相遇，即使得满足式（1）的时刻t存在。

2 基于共面圆轨道假设的分析

连续多次木卫引力辅助捕获需要找到一个时间窗口，窗口内各木卫相位差很特殊，可以使得航天器在半圈内依次经过多个木卫。共面圆轨道假设可以带来的便利是：

1） 时间窗口可以用木卫相位差反算；

2） 一条借力轨道可以通过轨道旋转被放置在不同的木卫相位处；

3） 一条多次引力辅助轨道可以由第一次引力辅助的入射速度、入射速度角和每次引力辅助的速度偏转角决定。

2.1 根据木卫相位差反算时间窗口

本节介绍利用木卫相位差反算出发时刻t0的方法。基于共面圆轨道假设，真经度L被用于表征木卫的相位。设给定时刻t*处4个伽利略木卫的相位分别为Lj，*，（其中，j为木卫序号，且j=1，2，3，4），伽利略木卫的圆轨道角速度为nj，t0时刻伽利略木卫的相位为Lj，则有

ΔLj=Lj+1-Lj， ΔLj，*=Lj+1，*-Lj，*，Nj=-N′j，ΔLk=Lk+1-Lk，ΔLk，*=Lk+1，*-，Nk=-，其中j=1，2，3，k=1，2，3，则有

若连续2 次木卫引力辅助的相位已知，可以设定Nj，用圆轨道的性质直接计算得到t0。若连续3 次或4次木卫引力辅助，则需在一定范围内遍历整数Nj，应用式（9）求得相应的Nk，若其为整数，则用Nj和圆轨道的性质可直接计算得到t0；若其不为整数，说明相应Nj取值下不存在连续引力辅助轨道。

2.2 轨道的几何旋转

设给定轨道的引力辅助木卫的相位为0，即Ω+ω+f=0，其中Ω为升交点赤经，ω为近地点幅角，f为真近点角。若实际上引力辅助木卫相位为L，可通过旋转引力辅助捕获轨道，即将rG、vG、r、v-、v+均乘以Givens旋转矩阵得到真实的运行轨道。实际的L取值由时间窗口决定。Givens旋转矩阵如下：

2.3 确定共面多次引力辅助轨道所需的变量

确定一条多次引力辅助捕获轨道所需要确定的参数包括引力辅助木卫序列、首次引力辅助的入射速度大小v-及其方向φ′以及各次引力辅助的转角δi，其中i=1，2，…，N，N为引力辅助次数。

第一次引力辅助前航天器的速度可以由其大小v-和角度φ′∈[0，2π]确定，即

设定木卫的相位L，即可得借力时刻的木卫速度，从而可以得到借力前航天器与木卫的相对速度v-∞，而借力后的相对速度可以用式（12）表示。

式（12）是应用φ=270°简化式（3）得到的。其中sinδ＞0，可由式（5）保证。共面假设使得φ=90°或φ=270°。若引力辅助后轨道周期变小，则φ=270°。这是由于i为的方向，j为vG的与i垂直的分量方向，而向j的负向偏转可以保证结合式（13）即可保证式（13）中为vG与的夹角的最大值为δmax。

确定了引力辅助后航天器的轨道。本次借力和下一次借力的木卫相位差、下一次借力的入射速度，均可依据开普勒轨道的性质得到。

3 广度优先搜索

考虑木卫的入射速度和大小往往由日心二体轨道搜索确定，本文假定v-已知，对每个φ′值，确定各次引力辅助的δ，以得到引力辅助轨道。考虑前一次引力辅助的δ取值直接影响下一次引力辅助的δ的取值范围以及相应的引力辅助参数，本文采用广度优先树搜索计算、查找满足要求的引力辅助轨道。搜索中，每个节点表示一次引力辅助，搜索内容包括引力辅助时刻、该时刻航天器位置及引力辅助前后航天器的速度。搜索树的第n层节点，为第n次引力辅助的可行情况。第一层节点用遍历的方式决定：在φ′的可行区间中，按照设定步长选取一系列φ′值，对每个φ′在[0，δmax]内遍历所有δ值，若引力辅助后轨道和下一次引力辅助木卫的轨道可以相交，见式（14），则将该次引力辅助存为树节点，其中a和e分别为引力辅助后轨道的半长轴和偏心率，rnext为下一次引力辅助木卫的轨道半径。

生成第一层节点后，启动广度优先搜索方法遍历树中的节点，遍历过程也是树中节点的增删过程。访问某个节点时，需要判定搜索该节点是否可以生成下一层节点：若可，则在树中保留此节点；否则，从树中删除当前节点。若当前节点代表一条满足要求的捕获轨道，则保存之。广度优先搜索流程可简述如下：

1） 访问某节点，利用该节点的位置和输出速度计算引力辅助后航天器运行的开普勒轨道根数，包括半长轴a、偏心率e和当前真近点角L1，*。

2） 计算该开普勒轨道到达第二次引力辅助木卫半径r2时的真近点角L1，#。若L1，*沿逆时针旋转小于π的角度可以到达L1，#，即L1，*和L1，#满足式（15）中3 个方程之一，则认为可以连续引力辅助；若无法找到满足上述要求的L1，#，则认为无法连续引力辅助，转到步骤7）处理当前访问节点。

式（15）中，第三式为从θ1，#逆时针越过2π到达θ1，*的情况，其中p=a(1-e2)，为航天器轨道的半通径。

3） 用开普勒轨道知识求得L1，*到L1，#所需的时间Δt12。

4） 计算两次引力辅助天体在首次引力辅助时刻的相位差ΔL12，见式（16）；n2表示第二个引力辅助天体的角速度，可以用计算得到。若为第三次和第四次引力辅助，则给出所有引力辅助木卫在前一次引力辅助时刻的相位，并给出相应相位差。

5） 在给定时间区间内，用2.1 节的内容查找该相位差对应的时刻。若存在这样的时刻，说明所访问节点的子节点存在。

6） 查找相应子节点引力辅助后的开普勒轨道，若该轨道和下一次引力辅助木卫的轨道有交点，则为有效节点，保存到树中；若不相交，则将当前节点的子节点转入步骤7）进行处理。

7） 查看引力辅助后的轨道周期是否满足要求，若是，则保存引力辅助轨道，并删除该节点；否则，继续遍历下一节点。

4 仿真过程和结果

4.1 算例说明

本文采用第十届全国空间轨道竞赛中给定的航天器双曲线轨道，即在出发时刻t（0对应简化儒略日数MJD 为[62502，63867]），木心距r=675RJ（其中RJ=71 492 km，为木星半径），速度大小v=4 km s，航天器质量为2 500 kg。假定伽利略木卫的偏心率和倾角均为0，伽利略木卫的其余轨道参数如表1 所示；同时假设共面圆轨道是可以接受的，因为伽利略木卫实际的最大偏心率为0.009 384，最大倾角为0.465 °。

表1 伽利略木卫基本参数Tab.1 Basic parameters of Galilean moons

4.2 仿真结果

考虑航天器初始位置大于木卫四轨道半径，且轨道行进过程中木心距由大变小，航天器经过木卫轨道的顺序为“木卫四（Callisto）→木卫三（Ganymede）→木卫二（Europa）→木卫一（Io）”，用卫星首字母表示为“C-G-E-I”，这是本文设定的第一个引力辅助序列；考虑木卫二的引力系数较小，本文还将计算“C-G-I-E”的序列形式。

由于开普勒轨道能量不变，因此，可以得到半径rG处的航天器速度：

为保证得到多次引力辅助轨道，可以在首次引力辅助前施加脉冲以减小航天器速度，故首次引力辅助前航天器的速度v-不大于12.054 km s。遍历搜索v-∈[10.000 0，12.054 0]和φ′∈[0，2π]工况下满足式（14）的φ′取值，如图2 所示，其中搜索步长分别为0.005 和0.0001。而后，对每种输入速度进行广度优先搜索，其中v-从12.054 km s 开始以0.005 km s 为步长逐步减小，共搜索410个速度值；φ′在图2所示上、下限范围内进行遍历，步长为0.01。每个引力辅助的偏转角在[0，δmax]中平均取31个数据进行遍历。计算2次引力辅助轨道的时间时，设定Nj=0；计算3 次引力辅助轨道的时间时需应用2.1节中的式（9），计算时需遍历Nj（即第1和第3个引力辅助木卫的圈数差），Nj=0的一个时刻是已搜索的2 次引力辅助时刻；计算4 次引力辅助时，所遍历的Nj被设定为是第1 和第2 个引力辅助木卫的圈数差。3次和4次引力辅助中Nj的搜索范围均被设定为0～700，当所得Nk为“整数±0.01”时，认为存在引力辅助解。两种引力辅助序列下，不同初始速度v-搜索得到的轨道周期结果如图3 和图4 所示。图中结果基本符合如下规律，即横坐标v-越小，较小周期的轨道解越多。部分区域轨道有较大波动的原因一是搜索的步长较大，会漏掉一部分解；二是轨道的1～3次引力辅助的结果意味着该序列下，航天器无法在一圈内与所设定的下一个引力辅助木卫相遇，但并不代表没有相应的轨道解。

图2 φ′取值范围随v-的变化Fig.2 Variation of the value range of φ′ with respect to v-

图3 “C-G-I-E”引力辅助所得不同周期的轨道解的个数分布Fig.3 Distribution of the trajectories with different periods for the "C-G-I-E" gravity assists(GA)

图4 “C-G-E-I” 引力辅助所得不同周期的轨道解的个数分布Fig.4 Distribution of the trajectories with different periods for the "C-G-E-I" GA

由图3 和图4 可知，调换了最后两个木卫的借力顺序后，周期为50～200 天的捕获轨道对应的最大v-有较大不同，但4 次引力辅助对应所需的最大v-基本相同。

两种引力辅助序列下，引力辅助所得轨道的最小周期如图5和图6所示。图中，v-较大时捕获轨道的周期出现了远大于相邻值的点，这是因为v-越大，捕获的轨道越少，更容易因为搜索步长较大而遗漏了一些轨道。两种引力辅助序列下的轨道如图7和图8所示，可以看出，多次引力辅助可持续减小航天器轨道能量，从而将双曲线轨道捕获为椭圆轨道。

图5 “C-G-E-I”引力辅助后的最小轨道周期Fig.5 The minimum periods after "C-G-E-I" GAs

图6 “C-G-I-E”引力辅助后的最小轨道周期Fig.6 The minimum periods after "C-G-I-E" GAs

图7 “C-G-E-I”引力辅助（GA）结果示意，v-=11.70 km/s，4 次GA 后轨道周期为28.66 天Fig.7 Trajectories with "C-G-E-I"GA, v-=11.70 km/s,the obtained period is 28.66 days

图8 “C-G-I-E”引力辅助（GA）结果示意，v-=11.76 km/s，4 次GA 后轨道周期为31.72 天Fig.8 Trajectories with "C-G-I-E"GA, v-=11.76 km/s,the obtained period is 31.72 days

5 结束语

本文假设伽利略木卫运行在共面圆轨道上，采用广度优先搜索求解该平面上的连续多次引力辅助的轨道。共面圆轨道假设和轨道共面假设降低了广度优先搜索的复杂度。仿真结果说明，可以利用伽利略木卫引力辅助得到周期小于50天的绕木轨道。文中方法所得结果可以作为初值，求解实际伽利略木卫轨道下的小倾角绕木引力辅助轨道，作为木星探测轨道的一部分。为增加搜索效率，可以进一步分析文中搜索式的解的存在性。