大气制动是一种气动力辅助的航天器变轨方式^[1]。将转移轨道的近拱点置于环绕天体的大气内，利用大气阻力降低航天器的轨道速度，通过多次穿越大气，使轨道远拱点高度逐步降低至目标高度，此时在远拱点执行脉冲变轨，即可实现环绕目标天体的任务轨道。这个过程所消耗的推进剂质量远小于直接变轨消耗的推进剂，可以大量节省深空探测器宝贵的燃料，从而延长寿命周期。到目前为止，已经有3次火星探测任务成功使用了大气制动技术。2016年欧洲太空局发射的ExoMars任务的火星微量气体轨道器(TGO)目前正在实施大气制动。应用结果表明大气制动技术对于深空探测任务的价值极大^[2]。在未来也是深空探测技术的重点组成部分。

近年来，随着航天技术的发展，大气制动在轨道建模、优化和控制等问题的研究越来越深入。大量的研究集中在大气层内飞行轨迹的优化和制导方面，如在满足热流的条件下使得探测器消耗的燃料最少^[3-4]，为防止热流过大对近拱点高度进行控制^[5]，在气动参数变化的情况下能跟踪设计轨迹等^[6-7]。还有一部分研究涉及到轨道建模，如近拱点选择对大气制动过程的影响^[8]，大气制动期间的气动模拟方法^[9-10]，时间消耗与探测器、轨道参数的关系^[11]，以及旋转和静止大气环境下的轨道特性^[12]等。但都没有涉及到气动力辅助变轨期间三体引力的影响。

众所周知，太阳系中最有科学价值的天体是木卫二、土卫六和土卫二，而土卫六是太阳系中唯一拥有浓密大气层的卫星，同地球一样富含氮气，成为当前国际深空探测的热点之一^[13]。"卡西尼-惠更斯"计划的顺利进行更推动了对土卫六的探测活动。土卫六具有稠密的大气，在未来探测任务中可以在其附近实施大气制动以节约变轨的燃料消耗。

在过去火星等天体的大气制动轨道研究中，三体引力摄动比较微弱，或摄动引起的偏心率/近拱点震荡周期与大气制动的时间相比很长，三体摄动基本被忽略。然而，土卫六附近的探测器轨道受到土星周期为8 d(土卫六轨道周期的一半)的显著三体引力摄动^[14]，引起偏心率/近拱点的大幅震荡，会对探测器所处环境的大气密度产生显著影响，可能会加速大气制动的过程，也可能导致大气制动失败或探测器坠毁。因此，有必要对土星三体引力摄动对土卫六附近大气制动的影响进行详细研究。

本文采用基于Milankovitch参数的半解析轨道方程来研究探测器在环绕土卫六大气制动过程中土星三体引力摄动的影响。首先，考虑扁率摄动、三体引力摄动以及大气阻力，基于无奇异的Milankovitch参数建立了半解析轨道模型。然后，通过数值计算分析土星初始方位角选择对土卫六大气制动轨道运动的影响，为地外天体大气制动变轨方案提供参考。

1 轨道动力学模型及分析 1.1 轨道动力学建模

对于大气制动的轨道来说，大气阻力与三体引力摄动均为长期影响，适用于平均轨道要素建模。将长周期和短周期摄动隔离开，就得到摄动影响和轨道变化的直接关系。而Milankovitch轨道要素具有简洁的向量形式和对物理关系的显式表达，没有奇异，有利于摄动影响的讨论^[15]。

大气制动过程中，探测器围绕卫星运行，将反复穿越卫星的大气层，同时卫星绕行星运动。卫星、探测器和行星的几何关系示意图如图 1所示，p为卫星自转轴方向矢量，H为探测器轨道角动量矢量，e为探测器轨道偏心率矢量，d_S为行星相对于卫星的位置矢量，θ为行星相对于轨道拱线的方位角，r和v分别为探测器的位置和速度向量。

将角动量H和偏心率矢量e写成位置r和速度v运算的形式

式中：μ为卫星的引力常数；r为探测器距卫星的距离。

对H和e建立高斯方程，并且将其平均化，有^[15]

式中：M为平近点角；a_p=a_J2+a_S+a_atm为摄动力，包括扁率摄动力a_J2、行星三体引力摄动力a_S和大气阻力a_atm，定义分别为^[15-16]

式中：J₂为卫星扁率；R_E为卫星平均半径；和分别为位置和卫星自转轴单位向量；μ_S为行星引力常数；ρ为大气密度，此处采用指数模型；δ=SC_D/m为探测器的弹道系数，S、C_D、m分别为探测器的有效面积、气动系数和质量；v_atm为大气速度矢量；ρ_p0为r_p0处的大气密度；H_atm为大气标高。

将式(5)~式(8)代入式(3)和式(4)经过运算化简后可得到与文献[15]类似的结果：

式中：和分别为d_S和e的单位向量，，e为偏心率；，a为轨道半长轴；z=ae/H_atm；；为轨道运动的平均角速度；。

1.2 摄动力影响分析

考虑式(9)、式(10)等号右边第1项和，均与e垂直，因此扁率只会影响轨道平面和半长轴进动，而不影响轨道大小变化。

同样地，考虑式(9)、式(10)等号右边第3项和，大气阻力只会引起半长轴和偏心率的减小，而不影响轨道平面方向变化。

而行星三体引力摄动的分析就比较复杂。由可以看到偏心率变化的大小和符号由土星相对于椭圆轨道拱线的方位角θ=〈d_S, e〉(见图 1)决定。以土星为例，该依赖关系可近似由正弦曲线描述，如图 2所示。

土星的运动导致其相对于轨道拱线的方位角不断增大，从而导致椭圆轨道的偏心率发生震荡。由图 2可知，该震荡的周期为土卫六轨道周期的一半，且震荡的初始相位由土星的初始方位角决定：若土星初始方位角为0°，则初始偏心率为震荡的最高点；而若土星初始方位角为90°，则初始偏心率为震荡的最低点。

因此土星三体引力摄动会引起轨道偏心率的震荡，从而使轨道的近拱点抬高或者降低。对于有大气的天体来说，这影响着轨道的寿命。而对于大气制动过程来说，近拱点的抬高或降低意味着探测器穿越大气层时利用大气制动效果的减弱或者加强。

2 数值仿真及分析

按照轨道模型式(9)、式(10)，本文给出一组不同初始方位角情况下的土卫六探测器大气制动轨道数值仿真，并分析土星三体引力摄动对大气制动过程的影响。大气制动可以节省的燃料消耗主要与轨道的初始和结束状态相对目标轨道的关系有关，不在本文的讨论范围，因此在仿真中并未考虑。仿真参数如表 1所示，其中大气密度参数根据文献[17]数据拟合得到。

表 1中轨道根数与Milankovitch参数的转换关系可参见文献[15]。仿真时长为200 d，仿真步长为10 s，并设置偏心率和近拱点高度阈值分别为0.008和100 km，如果偏心率或轨道近拱点高度达到此值，中断仿真。

下面暂不考虑大气阻力，给出不同土星初始方位角时的轨道运动如图 3所示。

从图 3可以发现，不同的土星初始方位角将导致明显不同的偏心率运动范围和近拱点高度运动范围。当土星初始方位角从30°增加到50°时，近拱点高度震荡中心明显变小。这意味着如果在有大气的情况下，初始方位角为50°的轨道将在穿越大气过程中经历更高的大气密度，减速效果更强烈。

考虑大气阻力后，不同土星初始方位角情况下的轨道运动如图 4所示。

从图 4可以发现，在三体引力摄动和大气阻力的双重作用下，不同土星初始方位角的轨道偏心率运动和近拱点高度运动范围和趋势明显不同。图 4(a)表明在土星初始方位角为50°时，轨道偏心率很快达到了结束阈值，而30°情况下时大气阻力作用并不明显。从图 4(b)可以看出土星初始方位角为50°时近拱点高度震荡谷值最低，气动作用最明显，同时近拱点高度震荡的峰值也降低，这是由于气动阻力的作用，轨道半长轴也在减小，相应地近拱点也会降低。图 4(c)显示随着时间的变化，大气制动作用越来越快。考虑大气制动最关注的指标热流率，此处v为探测器速度大小，结合表 2可以看到，土星初始方位角为50°的情况下，虽然大气制动经历的时间最短，但是热流率远高于其他情况，这将指导大气制动过程初始方位的选择。同时可以看到，远拱点和近拱点高度的波动状态，与之前在二体运动假设下的大气制动过程中单调下降的趋势不同，将影响大气制动退出条件以及燃料消耗。

在大气制动过程中，如果不考虑土星三体引力摄动，相同的轨道初始条件，期间不加入远拱点机动降低近拱点高度，大气制动过程将超过数十年。

3 结论

通过轨道动力学建模和数值仿真讨论了土星三体引力对土卫六大气制动过程的影响，仿真结果表明：

1) 土星的三体引力摄动对土卫六大气制动过程的影响较大，会极大地影响大气制动过程。在土卫六大气制动设计中不能忽略土星三体引力摄动。

2) 三体引力摄动对大气制动的影响主要取决于土星的初始方位角。初始方位角的不同导致椭圆轨道的偏心率震荡范围的不同，从而引起近拱点高度震荡范围的变化，导致不同的大气制动效果。

3) 本文采用了Milankovitch参数和平均轨道要素的方法来计算大气制动过程，拓宽了大气制动的仿真手段，为未来进一步开展地外天体气动力辅助变轨技术的研究奠定了良好的基础。