长期拟周期近直线晕轨道高效设计方法

秦理民，杨洪伟，李 爽

（南京航空航天大学航天学院，南京 211106）

0 引言

NASA 的阿尔忒弥斯（Artemis）项目［1］计划将人类重新带回月球，使得目前国际上对于近月空间的探索成为了新的重点。同时，一个被称为月球网关［2］的新型空间站将被放置在月球轨道上，近月空间站建设可以大大提高人类对深空探测的长期探索能力，包含太阳系其他行星和小行星探测任务。未来我国也将进入对地月空间大规模开发与探索的新阶段以及进行新的深空探测任务［3-4］。

月球空间站的计划轨道之一是近直线晕轨道（Near rectilinear Halo orbit，NRHO）［2］，它是地月共线平动点halo 轨道的一个子族［5］。在圆型限制性三体问题（Circular restricted three-body problem，CRTBP）中NRHO 被确认是在月球附近有良好稳定性的周期轨道族，且这一轨道族相比于halo 轨道的其余部分有更低的近月点半径。NASA 目前选定的月球网关轨道为L2 点南族9∶2 会合共振NRHO［2］，文献［6］中分析了NRHO 作为月球空间站的适用性。近年来，NRHO 与halo 轨道在地月空间的转移任务也得到了广泛的研究［7-9］。

不同动力学模型下会合共振NRHO 的研究受到广泛关注，Boudad 等［10］分析了会合共振NRHO 在双圆限制性四体问题下对应周期轨道的动力学特性。Villegas-pinto 等［11］在考虑月球偏心率和太阳引力的椭圆-圆限制性四体问题下分析共振拟周期NRHO。此外，Zimovan-spreen 等［12］对NRHO 和其周围的高倍周期分岔轨道族的动力学特性进行了分析。这些研究都表明，会合共振NRHO 在多种动力学模型下有效保留了避免月球和地球遮挡以及近线性稳定的动力学特性。

但是，在高精度的星历模型下，NRHO 将演变为拟周期轨道。文献［13］中给出了圆型限制性三体问题下大幅值拟周期轨道的数值计算方法，但未考虑多摄动影响下的高精度模型。利用数值方法在高精度星历模型下构造拟周期轨道，包含二级微分修正方法［14-15］和多点打靶法［16-17］。但这些研究多关注于平动点周围的halo轨道和Lissajous轨道且往往只计算了少圈次的拟周期轨道［14-17］。目前对于在高精度模型下的NRHO 维持和导航已有相关研究分析［18-20］，其结果验证了在真实星历模型下维持拟NRHO 运动的可行性。虽然其中也有提及长期拟周期NRHO 的生成过程，但其中计算初始多圈次拟周期NRHO 包含了数值敏感区域的离散点，增加了多点打靶的求解难度。Williams 等［18］提出了滚动时域方法计算长期拟周期NRHO，但其中只提到从一个良好收敛的初始解开始可以保证滚动时域方法计算长期轨道的有效性，并未明确指出不同NRHO 初始生成解对于长期拟周期NRHO 计算的影响，同时并未给出滚动时域算法中具体的修正策略。且目前多关注于L2点南族NRHO［18-20］，缺少对于L2点北族和L1点NRHO的适用性分析。

为此，针对现有的长期拟周期NRHO 数值计算方法存在的不足，本文提出了新的离散点选取策略，提高了多点打靶法计算多圈次拟周期NRHO 的收敛性能并能有效降低迭代打靶所需的计算时间。进一步利用NRHO 的对称性给出了初始拟周期NRHO 圈次延续方法和滚动时域长期拟周期NRHO生成方法，同时提出依赖拟周期轨道近月点半径判断的修正策略，以保证滚动时域方法的可延续性和更少的迭代修正次数。在仿真算例中，本文将对地月L2 点NRHO 和L1 点NRHO 进行分析，评估拟周期NRHO 的近月点半径振荡对长期生成的影响，以及长期拟周期NRHO生成方法的有效性。

1 动力学建模与轨道选取

1.1 圆型限制性三体问题

圆型限制性三体问题通常在会合坐标系下描述航天器的运动，航天器的状态量为，其中rsyn和vsyn分别为航天器在会合坐标系下的位置和速度，其动力学方程为［21］：

其中，U代表系统的伪势能：

式中：μ为系统的质量参数；r1和r2分别为航天器相对于地球和月球的距离：

1.2 高精度星历模型

高精度星历模型很好地描述了航天器在多天体引力下的运动，此模型建立在以行星为中心的J2000 惯性系下，天体的位置和速度通过读取NASA的导航和辅助信息设施（NAIF）提供的SPICE 星历数据库［22］获得。

以地球为中心天体，航天器的状态量为X=其中r和v分别为航天器在惯性系下的位置和速度，航天器在多个天体引力作用下的动力学方程为：

式中：G为万有引力常数；me为地球质量；mi和ri分别为月球、太阳以及其他星体的质量和位置矢量。

定义瞬时归一化长度单位l*、时间单位t*和速度单位v*为

则航天器在地球J2000 惯性系下的状态量rin和vin转换为归一化的瞬时地月会合系下的位置和速度：

式（7）中rsyn，e=[-μ，0，0]T为瞬时地月会合坐标系下地球的位置坐标。通过上述转换的逆过程可将航天器在地月会合坐标系下的坐标转移至地球J2000惯性系下。

1.3 轨道选取

NRHO 是月球附近具有良好稳定性的周期轨道族，它被认为是月球门户和月球空间站建设的重要潜在应用轨道。地球-月球-太阳的会合周期为29.5 天，文献［10-12］论证了会合共振NRHO 能有效避免月球和地球遮挡。本文所选取的NRHO 参数如表1所示，所对应的轨道分别为L2点9∶2共振、L2 点4∶1 共振和L1 点11∶3 共振。表中NRHO 的状态量以会合坐标系下的无量纲单位表示。

表1 所选NRHO的轨道参数Table 1 Orbital characteristics of the selected NRHOs

2 长期拟周期NRHO生成方法

虽然NRHO 在圆型限制性三体问题下表现出良好的稳定性，但在真实情况下，受到月球偏心率和其他天体引力的影响，拟周期NRHO 在一定圈次后即发散。但生成长期的拟周期NRHO 标称轨迹是轨道保持的关键，故本小节给出了三级修正的长期拟周期NRHO设计方法。

1）第一级修正：从圆型限制性三体问题转到高精度星历模型。

以圆型限制性三体问题下生成的NRHO 作为初始猜测，利用多点打靶法［16］构造拟周期NRHO。此过程的自由变量为各个离散点的状态量Xq(q=1，…，M)：

为保证轨道的连续性，约束表示为

约束变量F1对于自由变量C1的雅可比矩阵为

式中：Φ为6 × 6的状态转移矩阵；I表示单位矩阵。

由于自由变量C1与约束变量F1的维度不相等，故使用最小范数解对自由变量进行迭代更新：

式中：K=0，1，2，…为迭代次数。

本文考虑到月球附近离散点对参数的敏感性，提出了离散点删除策略，对初始离散点选取进行改进，提高多点打靶过程对于多圈次拟周期NRHO 的收敛性。首先考虑将每圈次的NRHO 等时间离散为10 段，然后通过x轴坐标逐步删除靠近近月点位置的离散点，离散点选取方案示意图如图1 所示。

图1 离散点选取示意图Fig.1 Schematic diagram of discrete point selection

2）第二级修正：初始拟周期NRHO圈次延续。

为延续第一级修正获得的拟NRHO 圈次，第二级修正利用圆型限制性三体问题下NRHO 的对称特性，即NRHO 每次通过近月点位置时=0（此速度为在地月会合坐标系下的速度），以第一级修正获得的NRHO 作为初始解，通过修正X0使N圈次后的近月点速度满足对称性条件。此过程的自由变量为初始的三轴速度量：

约束为N圈次后近月点的x轴速度，表示为：

N圈次后的近月点在会合坐标系下的速度vsyn，p，N对于自由变量C2的偏导数为：

式中：Φi，j为状态转移矩阵Φ中的元素；Rt，和分别为末端近月点时刻所对应的旋转矩阵、角速度矩阵和瞬时归一化速度单位；Q为3 × 3的矩阵。

约束变量F2对于自由变量C2的雅可比矩阵为：

式中：下标i，j为Q中对应元素的位置。同样使用最小范数解对自由变量C2进行更新：

通过该修正过程可以获得N圈次下拟稳定的轨迹，将该轨迹延续至N+1 圈次，进一步修正初始速度即可获得N+1 圈次的拟周期NRHO 轨迹，重复上述过程以延续初始拟周期NRHO圈次。

3）第三级修正：滚动时域方法生成长期拟周期NRHO。

当获得理想圈次的拟周期NRHO 轨迹后，采用滚动时域方法计算长期拟周期NRHO。滚动时域方法开始于初始近月点Xp，0，其向前积分N圈次到新的近月点Xp，N，通过修正初始近月点的速度使=0。然后将修正后得到的初始近月点向前积分到下一个近月点Xp，1，重复上述过程直到获得理想圈次的轨迹。

由此，第三级修正的自由变量为k圈次时的近月点的三轴速度量：

约束为k+N圈次后的x轴速度，表示为

约束变量F3对于自由变量C3的雅可比矩阵同式（16）。

然而，近月点速度的修正只能保证后续N圈次是拟稳定的，同时会导致N圈次后近月点位置和历元时刻的改变，只包含一次修正的滚动时域方法对部分NRHO 会在一定圈次后发散。对近月点速度的多次修正可以解决轨迹发散的问题，但多次修正意味着计算量的成倍增加。本文引入当前修正N+1圈次后近月点半径Rp，k+N+1的约束来对后续发散程度较大的近月点Xp，k进行多次修正来保证算法的收敛性，同时避免不必要修正所带来的计算量的增加。设置允许的最大近月点半径Rp，max和最小近月点半径Rp，min，根据其发散特性给出了如下不同的修正策略：

1）Rp，k+N+1＞Rp，max。

此情况下N+1 圈次后的轨迹已存在较大发散，首先通过多次修正Xp，k使N圈次后近月点满足式（19）来减小Rp，k+N+1，同时设置最大修正迭代次数nit-max1。若多次修正后Rp，k+N+1仍不满足约束，则多次修正Xp，k使N+1 圈次后近月点满足式（19），同时设置最大修正迭代次数nit-max2。

2）Rp，k+N+1＜Rp，min。

此情况表明N+1 圈次内的轨迹没有发散，但近月点的降低会导致N+2 圈次后的轨迹发散而使滚动时域方法失败。此时通过多次修正Xp，k使N+1 圈次后近月点满足式（19）来抬升Rp，k+N+1，同时设置最大修正迭代次数nit-max3。

多数圈次在增加一次修正后可以满足近月点半径要求。但由于在真实星历模型下，拟周期NRHO 近月点半径发生振荡的动力学特性，对于初始近月点状态的多次修正可能也无法使多圈次后的近月点半径满足约束，故需设置最大修正迭代次数以结束修正，滚动时域方法流程图如图2 所示。综上所述，长期拟周期NRHO 生成方法图解如图3所示。

图2 滚动时域方法中修正策略流程图Fig.2 Flow chat of correction strategy in the receding horizon method

图3 长期拟周期NRHO生成方法图解Fig.3 Diagram of long-term quasi-NRHO generation method

3 仿真试验和结果分析

仿真以1.3 节所选取的NRHO 为例计算拟周期NRHO。起始时间为“2025.01.01.000 UTD”，考虑日地月木四天体引力作用下的高精度星历模型。

3.1 第一级修正算例分析

所选目标轨道在第一级修正过程中离散点选取与收敛圈次的比较如表2 所示。此结果表明，减少单圈次离散点的策略能有效增加多点打靶法的收敛圈次。同时该策略能有效减少迭代次数和计算时间，对于计算三圈拟周期NRHO（1）收敛轨迹，使用MATLAB 中的ode 45 积分器进行轨道积分，绝对误差和相对误差分别设置为1×10-14和1×10-12。采用单圈次5个离散点的策略求解时间和迭代次数分别为13.191 s 和6 次，而采用单圈次10 个离散点则需要48.898 s和19次迭代。

表2 离散点选取与收敛圈次的比较Table 2 Comparison of discrete point selection and converged revolutions

以NRHO（1）为例的第一级修正收敛解如图4所示，此结果绘制在归一化的地月会合坐标系下，虚线轨迹为NRHO，圆点对应不同收敛解下所选取的离散点，实线轨迹为拟周期NRHO。

图4 NRHO（1）第一级修正收敛解Fig.4 Converged solutions for the first-level correction of NRHO（1）

3.2 第二级修正算例分析

虽然第一级修正给出的离散点选择策略能有效增加多圈次拟周期NRHO的收敛性能，但由于删除近月点周围离散点导致NRHO部分信息的缺失，有可能导致拟周期NRHO近月点漂移。第二级修正可以有效解决这一问题，选择圈次较少但近月点偏离较少的收敛解进行延续可以很好地保证近月点半径的变化范围。对于NRHO（1），选择10 圈收敛解和选择15 圈收敛解生成维持21 圈轨迹的近月点半径变化如图5（a）所示；对于NRHO（2），选择7圈收敛解和选择9圈收敛解生成维持18圈轨迹的近月点半径变化如图5（b）所示。此结果表明，第二级修正策略维持了近月点半径更小幅度的变化。与此同时，近月点半径更少的偏移也保证了第三级修正过程不易发散。

图5 不同初始解下近月点半径变化Fig.5 The change of perilune radius with different initial solutions

从上述第一级修正中选择合适圈次的拟周期NRHO收敛解作为初始值，利用第二级修正方法延续拟周期NRHO 圈次，最终生成维持22 圈140 天的拟周期NRHO（1）轨迹、24圈170天的拟周期NRHO（2）轨迹和14 圈110 天的拟周期NRHO（3）轨迹如图6所示，实线轨迹为拟周期NRHO，圆点为月球。上述结果表明，第二级修正能有效弥补第一级修正收敛圈次的不足，生成更多圈次的拟周期NRHO。

图6 NRHOs第二级修正收敛解Fig.6 Converged solutions for the second-level correction of NRHOs

3.3 第三级修正算例分析

生成了足够圈次的初始拟周期NRHO，利用第三级修正滚动时域方法进行长期拟周期NRHO 计算。为验证本文中所提出的修正策略的有效性，同时利用单次修正、多次修正和本文修正策略对滚动时域方法可生成圈次进行了比较。其中单次修正指只对Xp，k进行一次修正使N圈次后近月点满足式（19）；多次修正是指对Xp，k进行一次或二次修正使N圈次后近月点满足式（19），再对Xp，k进行一次修正使N+1圈次后近月点满足式（19）（当有一种情况不发散时即停止）。分别对3.2 节中初始10 圈延续得到的拟周期NRHO（1）、初始7圈延续得到的拟周期NRHO（2）和初始9 圈延续得到的拟周期NRHO（2）进行了长期计算，其结果如表3所示。

表3 不同修正策略下滚动时域方法性能对比Table 3 Performance comparison of the receding-horizon method with different correction strategies.

由上述结果可以看出，对于L2 点9：2 会合共振的NRHO（1）来说，单次修正即可保证轨迹不会发散。而对于4∶1 共振的NRHO（2），即使选取了近月点漂移较少的第二级修正收敛解，单次修正依然会在少圈次后发散。这样的结果同时表明了9：2会合共振轨道作为网关轨道更好的适用性。对于近月点半径漂移较少的拟周期NRHO（2），多次修正可以保证轨迹不发散，但至少需要在每圈次近月点实施3 次修正计算，这意味着计算500 圈的长期拟周期轨道需要1 500 次迭代修正过程，而本文提出策略则只需要1 040 次，可大幅度减少迭代次数，提高计算效率。在利用滚动时域方法计算近月点半径漂移较大的拟周期NRHO（2）的长期结果时，多次修正策略仍会失败，这是因为当Rp，k+N+1＜Rp，min时，对N圈后近月点的修正会使N+1圈次后的近月点半径Rp，k+N+1进一步降低从而导致后续轨迹的发散，而根据近月点半径判断的本文策略仍可以延续至所需圈次。

除此之外，本文提出的依赖近月点半径的修正策略因为约束了近月点半径，确保了多数圈次内拟周期NRHO 近月点的振荡范围，使拟周期NRHO 的轨道特性更接近圆型限制性三体问题中的NRHO。

最终生成维持500圈约10年的拟周期NRHO 如图7 所示，实线轨迹为拟周期NRHO，圆点为月球。对于NRHO（1）选择N=12维持500圈所需的速度增量为58.55 mm/s，近月点处最大脉冲为0.358 9 mm/s。对于NRHO（2）选择N=10维持500圈所需的速度增量为29.65 mm/s，近月点处最大脉冲为0.291 2 mm/s。此结果比以往研究中L2点南族NRHO维持500圈次所需速度增量更小［18］，说明了本文方法的有效性，同时可能说明L2 点北族拟周期NRHO 具有比L2 点南族拟周期NRHO 更接近圆型限制性三体问题下NRHO的轨道特性。

图7 长期拟周期NRHOsFig.7 Long-term quasi-periodic NRHOs

4 结论

本文针对地月空间中高精度星历模型下长期拟周期NRHO 设计问题，给出了一种新型三级修正策略的长期拟周期NRHO 设计方法，对L1 和L2 点的NRHO 进行了仿真计算。计算结果表明，提出的离散点策略不仅降低了多点打靶法的收敛难度，同时提升了计算效率。第二级修正策略克服了初始圈次不足的困难，第三级修正通过引入近月点半径判断的修正策略保证了滚动时域方法不发散且减少了修正迭代次数。最终生成了维持500圈约10年的拟周期NRHO 轨迹，速度消耗在几十毫米量级。以上设计方法和设计结果可以为长期拟周期NRHO设计提供重要参考。