CJA
设一个三体系统中，第三体的质量远小于其他2个天体，因而不考虑第三体对于其他2个大天体的引力作用，其他2个大天体在相互引力作用下作Kepler运动，这样的三体系统构成限制性三体问题。若2个大天体作圆运动，则称为圆型限制性三体问题，若作椭圆运动，则称为椭圆型限制性三体问题。圆限制性三体问题中可求得5个平动点，包含2个三角平动点和3个共线平动点，其中，平动点附近存在着大量周期轨道或拟周期轨道，地月三角平动点可以被应用于导航星座、空间中转站等深空探测任务中，具有重要的应用价值。
近似解析法和数值法是目前求解平动点附近周期轨道和拟周期轨道的常用方法。Suryawanshi［1］构造了共线平动点附近Halo轨道的近似解析解，Richardson［2］利用连续逼近法进一步给出了其三阶解析解。翟冠峤等［3-4］和钱霙婧等［5］研究了三角平动点平面和垂直周期轨道各个运动方向之间的近似解析关系。Liang等［6］利用李级数法研究了三角平动点的动力学结构，并给出了长周期轨道的近似解析解。Meng和Chen［7］、Connor Howell［8］以及Pan和Hou［9］基于圆型限制性三体模型和平动点周期轨道的对称性，利用微分修正法求解了Halo轨道等共线平动点周期轨道。为进一步提高算法收敛域，Gómez等［10］提出了多步打靶法。雷汉伦［11-12］等基于多步打靶法求解得到了真实力模型下的共线平动点拟周期轨道，但在求解过程中发现平动点附近轨道较为敏感，长时间打靶难以收敛。除上述方法，针对近月空间，Canales等［13］利用插值与多项式拟合的方法给出了李亚普诺夫轨和远距离逆行轨道等周期轨道计算方法，比传统迭代算法效率提高50%。Peterson和Scheeres［14］通过构造局部角轨道根数，简化了平动点附近的动力学，对于周期轨道和拟周期轨道的求解问题提供了一种对数值方法的实用补充。Beom和Kathleen［15］对比了圆型限制性三体模型、椭圆型限制性三体模型和双圆型限制性四体模型下地月L2点Halo轨道在摄动项上与高精度星历模型的接近程度。有学者［16-18］考虑太阳光压等摄动力，构造了平动点周期轨道近似解析解与数值解。Peng等［19］通过参数化方法逼近了在太阳光压摄动下地月三角平动点附近的稳定流形和不稳定流形。El Motelp和Radwan［20］研究了摄动力对三角平动点周期轨道的周期、方向和偏心率的影响。
部分学者研究了限制性三体问题平动点的稳定性。Vincent等［21］分析了系统参数对共线平动点与三角平动点的稳定性的影响，并基于数值法给出了2001SN263小行星三体系统共线平动点周围的周期轨道。刘林等［22-23］通过研究发现如果要选取相应的拟周期轨道作为探测器定点在三角平动点附近的目标轨道，需要尽量消除其中的长周期运动分量且必须根据状态偏离的规律考虑在轨运行过程中进行轨道控制。Li和Hou［24］、Limebeer和Sabatta［25］针对地月系统中三角平动点的不稳定运动，分别提出了控制策略。为了降低控制频率，部分学者着手研究限制性三体问题的有界运动。Sosnyts’Kyi［26］、Liu等［27］给出了三角平动点空间稳定区边界形成的机制。Rodnikov［28］建立了在平动点附近进行有界运动的初始条件。He等［29］提出了三角平动点附近可长时间保持且不发散的自然有界相对轨道的数值搜索方法。目前，也有学者着手研究三角平动点周期轨道的应用问题。Meng和Chen［7］提出了地月空间星座设计的要素和步骤。Xin等［30］和Xu等［31］利用平动点附近周期轨道将组合导航应用于地月平动点附近的星座，实现了地月空间导航服务。张磊［32］得到了满足地月空间导航全覆盖要求的三角平动点周期轨道，刘斌等［33］考虑太阳光压摄动，引入动力学替代轨道，利用三角平动点周期轨道星间测距数据自主定轨。Qi等［34］利用三角平动点附近的轴向轨道设计了天基望远镜的观测策略，用以监测和搜索近地天体。
综上所述，平动点周期轨道的近似解析解与数值解计算已经较为成熟，然而，上述学者大都在圆型限制性三体模型下计算和分析三角平动点周期轨道，在高精度模型下三角平动点周期轨道的实现和特点分析还较少。高精度模型下短周期轨道计算的经典方法是并行打靶法，后来学者又提出了连续环绕拼接和多段连续拼接等许多类型的计算方法。然而，并行打靶法在计算多圈短周期轨道时计算量很大，收敛性难以保证，对求解器要求很高。有的学者尝试采用序列二次规划算法等优化算法来求解，但这样计算量更大。除此之外，三角平动点虽是稳定平衡点，但在摄动条件下轨迹仍易于发散，部分学者研究了轨道维持控制策略，但在高精度模型下，如何使得短周期轨道保持长时间有界运动的相关研究较少。本文将针对三角平动点的短周期轨道的实际应用轨道设计进行分析。
1 动力学模型
1.1 圆型限制性三体动力学模型
限制性三体问题仅研究第三体在两大天体引力作用下的运动，故至多是3自由度系统。
如式（1） 所示，为了分析问题和计算上的方便，采用归一化的单位。质量单位为天体质量之和m 1+m 2，长度单位为2个天体的距离a，时间单位为小天体绕大天体飞行的角速度的倒数1/n。
设 ，则在归一化单位下，万有引力常数G=1。大天体的质量为 ，小天体的质量为 。
如图1所示，质心惯性系原点位于两个天体质心处，其X轴指向春分点，Z轴垂直于天体公转轨道平面，方向与两天体总角动量方向一致，Y轴由右手法则确定，研究三体问题时，通常用到质心旋转坐标系，其X c轴始终由大天体指向小天体，Z c轴与质心惯性系重合，Y c轴垂直于X c轴指向小天体飞行方向。
在圆型限制性三体问题下，质心旋转坐标系中三体运动的微分方程为
式中：
1.2 高精度动力学模型
式中： 为航天器位置矢量，g和gi 分别表示地球和太阳等大天体i的引力常数， i 和 i 分别表示大天体位置矢量和航天器相对大天体的位置矢量， 表示光压摄动、地球非球形摄动、月球非球形摄动等摄动力。其中太阳引力摄动、光压摄动和地球J2项影响较大，其他因素相对较小。因此在轨道设计和分析时，高精度模型仅考虑日地月中心引力、光压和地球J2摄动。则式（3） 进一步可以写为
式中： 和 分别为航天器所受太阳光压力和地球J2摄动力，Pr 是距离太阳中心1个天文单位处的太阳辐射压力，数值约为4.6×10-6 N/m2，Cr 为航天器反射系数，A为航天器光照面积， s 为航天器相对于太阳的位置矢量，J 2为地球扁率摄动常数数值约为0.001 082 63，Re 为地球平均半径数值约为6 378 137 m，ψ为地心纬度。
2 圆型限制性三体问题下短周期轨道计算
当幅值很小时，短周期轨道近似为绕着三角平动点的椭圆，随着幅值增大，轨道逐渐扭曲成肾形。如图2所示，计算三角平动点轨道时建立三角平动点旋转坐标系（以L4为例）。
坐标系的原点为大天体M 1，坐标系X t轴的正方向由M 1指向L4，Y t轴与X t轴垂直指向坐标旋转方向。质心旋转坐标系到三角平动点旋转坐标系的转换公式为
式中： 为质心旋转坐标系下的位置速度状态， 为三角平动点旋转坐标系下的位置速度状态。对于L4中心的旋转坐标系θ取60°，对于L5中心的旋转坐标系θ取-60°。
在计算三角平动点短周期轨道时，在以L4为中心的旋转坐标系下，对于小幅值轨道仍可用穿越X t Z t面的x t速度为0初步计算，但三角平动点周期轨道并不是关于X t Z t面对称，精确计算需要积分一周，终点的x t方向位置和y t方向速度要与初始点相等。因此计算三角平动点短周期轨道时初值设置在X t轴上，给定x t坐标（y t坐标为0），待解变量为：x t和y t方向速度和轨道周期，即 。目标为一个周期后的x t方向位置、y t方向位置和y t方向速度与初始的偏差为0，即 等于0。采用微分改正法求解，初值用数值搜索得到，周期初值为2π。三角平动点为稳定平衡点，幅值为0则速度为0，对于幅值很小的情况下，初始速度很小，xt 方向的速度初值固定设置为0，y t方向速度初值用0.001作为步长搜索很快能够得到收敛解。在得到一条小幅值的短周期轨道后，采用延拓方法得到离散的所有短周期轨道参数，后续计算时根据xt 方向幅值插值计算即可。
三角平动点短周期轨道在白道面内，在以L4或L5中心的旋转坐标系下，当确定了X t轴坐标时，短周期轨道即被确定，因此三角平动点短周期轨道通常以X t轴坐标为延拓参数，如图2所示。
延拓计算是以易于求得的轨道出发，通过一点点改变轨道的某个参数（改变量即延拓步长），以最近点的解为初值微分改正计算精确解，直到延拓变量到达目标边界或者在延拓步长小于某个值的情况下微分改正仍无法收敛。为提高延拓的效率，根据微分改正动态调整延拓步长，并且可以通过曲线拟合插值计算下一步的初值。延拓流程如图3所示，具体如下：
1） 用近似解析解或者数值搜索结果作为初值。
2） 设置初始延拓步长为0.001，设置周期变化门限为0.5，设置符号系数为1。
3） X t轴方向位置增加一个步长乘以符号系数。
4） 以上一个轨道的解为初值微分改正计算当前的轨道。
5） 如果求解不收敛或周期变化量大于预设值，则延拓步长缩短一半且符号系数设置为-1并回到3），否则符号系数设置为1并继续往下。
6） 计算X t轴方向位置的变化量。
7） 如果连续n max次求解正常则延拓步长加倍；
8） 回到3）继续延拓。
如图4所示，是以X t轴坐标为延拓参数所得地月L4点平面短周期轨道族，延拓过程中沿X t轴方向延拓能力较为有限，观察延拓轨道发现：轨道逐渐绕大天体向L5延伸。因此尝试用坐标系的旋转角度为延拓参数，即角度延拓法。延拓方法是建立新的旋转系，坐标原点选择在距离大天体M 1长度为1的点
式中： 为新旋转坐标系下的位置速度状态，α为X a轴与2个天体连线的夹角，是延拓参数，初值为60°，延拓方向为逆时针。
由于L5点与L4点关于质心旋转系的XcZc 平面对称，因此只需要延拓L4点的短周期轨道即可。在计算L5点的短周期轨道时，首先根据幅值计算得到M 1中心L4点旋转坐标系的轨道参数，然后将X t方向速度反号即得到M 1中心L5点旋转坐标系的轨道参数。
3 高精度模型下短周期轨道计算
在高精度模型下，由于轨道摄动和月球轨道偏心率的影响，圆型限制性三体问题下的轨道基本不再是周期轨道，需要采用不同的算法计算连续多圈的轨迹以逼近周期轨道，以此来尽量保持原轨道形状特点。
如图6所示，是基于圆型限制性三体轨道初值并用并行打靶法求解5圈轨道，然后分别外推1年、5年和20年的地月系L4点短周期轨道的轨迹，黑色圆点是轨迹初始点，红色圆点是轨迹终点。
可见，如果没有修正所有圈的位置，在摄动因素影响下三角平动点短周期轨道会从较小的幅值（<1×104 km）开始逐渐发散到很大的幅值上（>2.5×105 km），这将给工程应用带来很多问题，如星座构型的破坏和不稳定性，并行打靶法不能约束轨迹的范围。因此本文提出修正节点速度使得轨迹包络面积最小的方法，以限制短周期轨道的运行范围。
设计变量为［x， y， vx， vy ］或［vx， vy ］，如图7所示，优化目标为外包络矩形的面积S或对角线L，因此在高精度模型下短周期轨道求解问题可以转化为非线性规划问题，其中优化目标如式（9） 或式（10） 所示。
具体计算步骤如下：
1） 用圆型限制性三体问题的短周期轨道作为初值。
2） 在轨道初值上增加修正量。
3） 积分计算整个目标时长内的轨迹，按一定间隔采样（2 h~1 d），计算采样点的Xt 和Yt 轴坐标，并记录其最大值和最小值x max、x min、y max和y min。
4） 计算矩形包络面积S或对角线L。
这种方法优化变量只有2个或4个，通过大量数值计算结果分析，速度的求解范围设为±30 m/s，位置修正范围设为5 000~10 000 km，可以得到基于一个节点的最紧凑短周期轨道，其中式（7） 不改变节点的位置，主要用于后续节点的计算或工程中的轨道控制。
用优化面积包络法可以计算10年（约120圈）的三角平动点短周期轨道，如果外推时间再增加，由于积分截断误差会导致后续轨迹发散，因此用优化面积包络法计算时，需要确保轨道能够在三角平动点附近，否则会造成求解空间扩大，轨道在地月空间里穿梭，难以得到最优解。然而，直接用圆型限制性三体模型下的短周期轨道作为初值不能保证轨迹在三角平动点附近环绕12圈以上，既不能确保优化包络面积法能够进行下去，因此为确保轨迹就能够在较长时间内保持在三角平动点周围，将并行打靶法计算4圈轨道的起始点作为初值。此外为了计算更长的轨迹，采用多段拼接的方法，其示意图如图8所示，其中每条线段代表环绕三角平动点若干圈的轨迹，黑色圆点代表节点表示上一次优化后的轨迹距离上一次起始点为4~6圈的点。
第1个节点采用k个优化变量进行计算，后续节点只优化部分优化变量。每次都优化计算N圈，然后将第M圈（一般为N/4或N/6）的位置速度作为一个节点，再次优化计算N圈，如此循环直到计算的圈数达到指定的圈数。虽然中间有很多优化计算且一次外推的时间很长，但是与并行打靶法相比，没有内层的两点边值问题求解，而且优化求解的变量很少，需计算的次数很少，所以计算量比并行打靶少，重要的是收敛性可以得到保证。如此计算多个节点即可拼接成一条10年以上的短周期轨道，节点的速度增量通常在0.2 m/s左右，即每年各个节点的速度修正量总和在0.5 m/s左右，这基本可以作为设计轨道使用。该方法理论上可无限计算短周期轨道，而且可以实时记录已完成节点的轨迹。
如果要求节点的速度差更小，则可以在此基础上进一步采用并行打靶法修正。由于节点的速度偏差已经很小，所以进一步的修正收敛速度很快，可以修正到每年0.001 m/s以下，而且节点的长度较长，所以节点数比直接并行打靶法少很多，相应的计算量也少很多。通过数值仿真可得，直接采用并行打靶法，用牛顿迭代求解15圈以上轨迹就难以收敛，用序列二次规划求解120圈轨迹则需要1 h以上，而采用优化面积包络法，采用经典的坐标轮换法等优化算法就能够在常规的单核电脑中，4 min内计算出一条120圈轨迹的短周期轨道。
如图9所示，是用并行打靶法+优化面积包络法+多段拼接法计算的30年的L4点和L5点短周期轨道，黑色圆点是轨迹初始点，红色圆点是轨迹终点。可见采用该方法可以稳定计算长时间保持稳定的三角平动点短周期轨道。此外通过数值仿真可知，该方法还可用于三角平动点中小幅值的垂直轨道的计算。
4 地月系短周期轨道特性分析
由图12分析可得，起始点x t坐标为5×104 km的轨道采用优化面积包络法逐渐趋于小幅值短周期轨道，x t方向在5×104 km内，y t方向在1×105 km内，最小值则接近零，即在幅值范围内广泛分布；此外轨迹的主要部分的外包络基本与圆型限制性三体问题轨道一致。起始点x t坐标为1×105 km的轨道，最后收敛结果是外包络基本与圆型限制性三体问题轨道一致，同样在包络内部广泛分布。起始点x t坐标为1.5×105 km的轨道，其外包络仍与圆型限制性三体问题轨道一致，但内部出现空腔，类似含有x t幅值从4×104~1.5×105 km的短周期轨道集合。此外随着外推时间增加，幅值为1.5×105 km的轨道会逐渐发散并飞离L4点。
从x t幅值为1.5×105 km的轨道的空腔得到启发，是否能够在限定外包络的同时使得内部空腔最大化，从而得到周期性较好的轨迹。为了更好地验证计算的正确性，采用遍历的方法来计算，并且预报时长2年以上。为避免轨迹发散，设定距离三角平动点的最大距离（大于限定距离的轨道终止计算并抛弃），计算每条轨道的最小距离，然后取其中最小距离最大的。采用该方法对幅值1×105 km的轨道进行计算，结果如图13所示。
进一步对小幅值短周期轨道对比发现，虽然轨迹散布范围有3.5×104 km，但不同初始历元的轨道最终稳定在1 000 km左右的范围内。如图14所示，是分别以2024年1月1日和2027年1月1日为起点计算的200圈地月系L4点短周期轨道的相对距离和位置曲线。
相对于3.5×104 km的轨道幅值，1 000 km的相对位置差是很小的量，说明最后稳定阶段各短周期轨道将在一个较小的范围内运动，这对于编队来说较为有利，除此之外，这个距离也能保证编队之间的航天器不发生碰撞。
5 结 论
首先建立了地月空间三体问题圆型限制性三体动力学模型与高精度动力学模型，并基于圆限制性三体动力学模型在旋转坐标系下给出了三角平动点短周期轨道延拓的具体步骤，以L4点为例，沿X t方向延拓能力较为有限，但在延拓过程中发现延拓轨道逐渐绕大天体向L5点延伸，为提高延拓效率，提出了用坐标系旋转角度为延拓参数的角度延拓法，并基于上述方法得到了更大范围的地月三角平动点短周期轨道族。进一步考虑日地月中心引力、光压和地球J2摄动，提出了通过修正节点速度使得轨迹包络面积最小的方法，并给出了并行打靶法法、多段拼接法与优化包络面积法相结合的混合优化方法，高效得到了30年的三角平动点短周期轨道。
根据上面的数值仿真可以得出，地月系三角平动点短周期轨道为散布在给定幅值范围的轨迹，优化包络面积法是计算三角平动点短周期轨道的非常有效和可靠的方法。由于地月系三角平动点大幅值短周期轨道散布范围大，不利于星座构型的保持，因此三角平动点轨道的幅值应该尽量地小。