CN106682274A - 考虑振幅约束的一种Halo轨道在轨保持方法 - Google Patents

Abstract

本发明是一种考虑振幅约束的Halo轨道在轨保持方法，属于航空航天技术领域。本发明通过在两个主天体和探测器构成的限制性三体模型下建立动力学方程，确定两个主天体和探测器构成的三体系统平衡点位置。确定大小天体和探测器构成的三体系统下平衡点附近的Halo轨道。根据扰动变量设计Halo轨道保持的微分修正算法。根据微分修正算法设计Halo轨道在轨保持策略。在真实星历环境下，考虑测量误差和执行误差，按照所述的轨道保持策略采用微分修正算法得到每次轨道修正的速度增量，按照速度增量进行轨道修正控制，可以实现满足振幅约束的Halo轨道在轨保持，同时尽可能降低轨道保持所需的燃料消耗。

Description

考虑振幅约束的一种Halo轨道在轨保持方法

技术领域

本发明是一种考虑振幅约束的Halo轨道在轨保持方法，适用于真实星历环境并考虑测量与执行误差条件下的平衡点Halo轨道在轨保持控制，属于航空航天技术领域。

背景技术

平衡点是探测器和两个相互环绕的主天体构成的三体系统中具有特殊性质的点，位于平衡点的物体相对主天体可以保持相对静止。平衡点附近存在满足一定速度位置约束的围绕平衡点运动的周期轨道，位于周期轨道上的物体同样相对大、小质量天体保持基本不变的位置关系。同时相比将探测器捕获至环绕轨道，捕获至平衡点周期轨道所需的速度增量相对较小。因此，以Halo轨道为代表的平衡点周期轨道非常适合行星探测和中继通讯任务，是未来探测器的理想工作地点。而在真实星历环境下探测器会受到其他天体的引力扰动，因此探测器在Halo轨道每隔一段时间需要进行在轨保持修正。

在已发展的关于Halo轨道在轨保持方法中在先技术[1](参见Howell K C,Pernicka H J.Station-keeping method for libration point trajectories[J].Journal of Guidance Control&Dynamics,1990,16(1):713-723.)给出了一种平衡点轨道在轨保持方法，该方法借鉴了行星际多目标探测转移轨道中途修正的计算方法，先根据任务需求将轨道分为多个转移段，再通过微分修正算法将转移轨道修正到标称节点上，从而得到每次变轨所需的最优速度增量。该方法虽然可以实现平衡点周期轨道的在轨保持，但每次轨道修正计算都需要标称轨道节点，因此该方法对于Halo轨道长期在轨保持存在较大的不便。

在先技术[2](参见Ulybyshev Y.Long-Term Station Keeping of SpaceStation in Lunar Halo Orbits[J].Journal of Guidance Control&Dynamics,2015,38(6).)给出了另一种Halo轨道长期在轨保持方法，该方法不需要在修正时得到标称轨道节点，只需保证在旋转坐标系下轨道每次穿过X-Z平面时X方向的速度为零。该方法得到的速度增量较小，计算简单，但随着在轨保持时间增长，Halo轨道的振幅会发生较大的改变。

发明内容

本发明的目的是为了提供一种考虑振幅约束的一种Halo轨道在轨保持方法。该方法基于真实星历模型，同时考虑探测器的测量误差和执行误差，在实现满足振幅约束的Halo轨道在轨保持前提下尽可能降低轨道保持所需的燃料消耗。

本发明的目的是通过下述技术方案实现的。

考虑振幅约束的一种Halo轨道在轨保持方法，通过在两个主天体和探测器构成的限制性三体模型下建立动力学方程，确定两个主天体和探测器构成的三体系统平衡点位置。确定大小天体和探测器构成的三体系统下平衡点附近的Halo轨道。根据扰动变量设计Halo轨道保持的微分修正算法。根据微分修正算法设计Halo轨道在轨保持策略。在真实星历环境下，考虑测量误差和执行误差，按照所述的轨道保持策略采用微分修正算法得到每次轨道修正的速度增量，按照速度增量进行轨道修正控制，可以实现满足振幅约束的Halo轨道在轨保持，同时尽可能降低轨道保持所需的燃料消耗。

考虑振幅约束的一种Halo轨道在轨保持方法，包括如下步骤：

步骤一：在两个主天体和探测器构成的限制性三体模型下建立动力学方程，确定两个主天体和探测器构成的三体系统平衡点位置。

限制性三体模型动力学系统中考虑质量可以忽略的探测器S₃在两个主天体S₁和S₂的引力共同作用下的运动，两个主天体S₁、S₂和探测器S₃三者质量关系为m₁＞m₂＞＞m₃。选择两个主天体系统的质心作为原点建立旋转坐标系，x轴方向由原点指向小天体，z轴为系统旋转的角速度方向，y轴与x、z轴构成右手坐标系。则探测器的运动描述为式(1)：

其中μ＝m₂/(m₁+m₂)表示系统的质量系数，m₁为较大主天体S₁的质量，m₂为较小主天体S₂的质量， 分别为运动探测器到两主天体的距离。

为获得旋转坐标系平衡点，令式(1)中速度和加速度等于0，共可获得五个平衡点(拉格朗日点)；所述五个平衡点均在x-y平面上，包括三个共线的平衡点和两个三角平衡点。在质心旋转系下三个共线平衡点的位置(x_Li,y_Li)分别为：

L1平衡点：L2平衡点：L3平衡点：

两个三角平衡点的位置分别为：

L4平衡点：L5平衡点：

步骤二：确定大主天体、小主天体和探测器构成的三体系统下平衡点附近的Halo轨道。

为获得更简洁的表达式，引入变量(ξ,η,ζ)，根据步骤一所得的平衡点的位置对旋转坐标系坐标进行替换，变量定义如下：

ξ＝x-x_Li；η＝y-y_Li；ζ＝z

平衡点附近的线性化运动方程可描述为式(2)：

式(2)中，ρ²＝x²+y²+z²；c₂(μ)、c_n(μ)为仅与质量有关的常数，可表示为式(3)：

式(3)中，γ为平衡点与主天体S₂的距离，p_n为n阶Legendre多项式。对方程(2)进行分析，平衡点附近Halo轨道线性项可表示为式(4)：

式(4)中，ω₀分别为周期轨道运动的频率，κ为常数，φ为Halo轨道相位；α、β分别为轨道平面和垂直振幅；但为了满足Halo轨道的周期性，要求平衡点附近Halo轨道的非线性项可表示为式(5)：

式(5)中，θ＝ω₀t+φ，i+j＝N称为解的阶数。利用式(4)和(5)作为平衡点附近Halo轨道的初值，通过多重打靶方法可以得到Halo轨道的精确数值解，即三体系统下平衡点附近的Halo轨道。

步骤三：根据扰动变量设计Halo轨道保持的微分修正算法；根据微分修正算法得到Halo轨道精确数值解。

引入扰动变量，相对于参考轨道的变分为(δx,δy,δz)，则六维状态矢量为：

则变分方程的状态空间形式表示为方程(6)：

与时间有关的矩阵A(t)是由4个3×3的矩阵构成，即式(7)：

式(7)中B(t)是一个时变函数，因此A(t)也是一个时变函数，I₃是3×3的单位矩阵，C为常值矩阵，满足如下表达式：

方程(6)的解的形式为式(8)：

这里Φ(t,t₀)是状态转移矩阵。状态转移矩阵是从初始状态到t时刻状态的线性映射，通过状态转移矩阵可估算出初始状态变分对轨道递推的影响。状态转移矩阵必须满足矩阵微分方程(9)：

方程(9)的初始条件为Φ(t₀,t₀)＝I₆，这里I₆为6×6的单位阵。在t时刻的单位状态转移矩阵，通过初始条件的数值积分得到。因为状态转移矩阵是6×6，所以方程(9)的递推可以通过36个一阶微分方程的积分得到。B(t)的元素取决于参考轨道，因此求解A(t)，方程(9)需要同时积分非线性方程，从而得到状态转移矩阵Φ(t,t₀)。

根据微分修正算法求解Halo轨道精确数值解一般需要通过多次迭代过程来完成，归结如下：

根据当前末端状态与理想状态的偏差，通过状态转移矩阵对初始状态进行修正，使得最终的末端状态满足要求的约束条件。选取六维状态矢量寻找初始状态使周期轨道在t_f时刻满足约束条件末端状态偏差与初始状态偏差的关系为式(10)：

矩阵的偏导数等于在t_f时刻的状态转移矩阵，那么方程(8)可以简化为方程(11)：

末端状态的偏差量为方程(12)：

为目标终端状态，为实际终端状态，将方程(12)代入方程(11)，可以解出一个初始状态修正量的估计值

由于Halo轨道关于x-z平面对称，周期轨道必定与x-z平面垂直，因此在交叉处有满足式(13)：

利用垂直平面相交的特点进一步对问题进行简化，为了满足Halo轨道在半个周期后满足垂直交叉的条件则需选择y_f＝0作为积分终止的条件，在利用方程(11)迭代过程中需要调整初始状态使得末端的状态和接近于0，从而得到三维周期Halo轨道的精确数值解。

步骤四：根据步骤三中的微分修正算法设计Halo轨道的在轨保持策略。

由步骤三可知理想Halo轨道保持需要满足垂直交叉条件，即满足但在真实星历环境下，严格满足约束的Halo轨道是不存在的，只需保证Halo轨道的振幅尺寸满足规定的范围即可。因此对于Halo轨道的在轨保持修正，可以不严格满足约束条件针对修正约束条件的严格性提出两种不同的修正策略：宽松修正方法和严格修正方法。

宽松修正方法是指：轨道积分穿过x-z平面时x方向的速度为零，即严格修正方法是指：轨道积分穿过x-z平面时x方向和z方向的速度均为零，也即

当轨道穿过x-z平面时z方向的振幅超过任务约束的任务上下边界时执行一次严格修正方法，否则执行宽松修正方法。同时在每两个轨道周期内执行的4次修正中至少执行1次严格修正策略。

步骤五：在真实星历环境下，考虑测量误差和执行误差，采用步骤四中的在轨保持策略进行Halo轨道在轨保持。

轨道保持修正采用步骤三中的微分修正算法。在求解修正脉冲时，积分初值为考虑测量误差的实际状态；而在轨道修正时，需要在实际状态进行脉冲修正时考虑修正的执行误差。在考虑测量误差和执行误差的条件下，按照轨道保持策略采用微分修正算法得到每次轨道修正的速度增量，按照速度增量进行轨道修正控制可以实现满足振幅约束的Halo轨道在轨保持，同时尽可能降低轨道保持所需的燃料消耗。

有益效果

1、本发明公开的考虑振幅约束的一种Halo轨道在轨保持方法，该方法不需要在修正时得到标称轨道节点，计算简单，适合探测器在Halo轨道的长期在轨保持。

2、本发明公开的考虑振幅约束的一种Halo轨道在轨保持方法，该方法根据两种不同修正策略的各自优势，设计了一种组合修正策略，可以实现满足振幅约束的Halo轨道在轨保持，同时尽可能降低轨道保持所需的燃料消耗。

3、本发明公开的考虑振幅约束的一种Halo轨道在轨保持方法，该方法基于真实星历模型，考虑探测器的测量误差和执行误差，Halo轨道在轨保持更加接近工程应用环境。

附图说明：

图1本发明方案流程示意图。

图2限制性三体模型旋转坐标系示意图。

图3限制性三体系统平衡点分布示意图。

图4本发明实例中在轨保持三维轨道示意图。

图5本发明实例中在轨保持轨道x-y平面投影图。

图6本发明实例中在轨保持轨道x-z平面投影图。

图7本发明实例中在轨保持轨道y-z平面投影图。

图8实例中先技术在轨保持三维轨道示意图。

具体实施方式

为了更好地说明本发明的目的和优点，下面通过对一个探测器在Halo轨道的在轨保持进行实例分析，对本发明做出详细解释。

实施例1：

本发明公开的考虑振幅约束的一种Halo轨道在轨保持方法，这里选择由地球-月球-探测器构成的限制性三体系统为例，平衡点选择L2点，方案流程示意图如图1所示，Halo轨道的在轨保持设计方法包括如下步骤：

步骤一：在地球-月球-探测器构成的限制性三体模型下建立动力学方程，确定地球-月球-探测器构成的三体系统平衡点位置。

选择地球和月球分别作为主天体S₁和S₂，在质心旋转坐标系下建立动力学方程。该坐标系选择地月系统的质心为原点，x轴方向由原点指向月球，z轴为月球公转方向，y轴与x、z轴构成右手坐标系，坐标系示意图如图2所示。则探测器的运动可以描述为式(1)：

其中μ＝m₂/(m₁+m₂)表示系统的质量系数，m₁为地球质量，m₂为月球质量，分别为运动探测器到地球和月球的距离。

令式(1)中速度和加速度等于0，经一系列推导和求解，可得五个平衡点(拉格朗日点)，这5个平衡点均在x-y平面上。这里选择L2平衡点作为Halo轨道部署位置，L2平衡点在x-y平面上位置坐标为：

平衡点的位置分别示意图如图3所示。

步骤二：确定大主天体、小主天体和探测器构成的三体系统下平衡点附近的Halo轨道。

为获得更简洁的表达式，引入变量(ξ,η,ζ)，根据步骤一所得的平衡点的位置对旋转坐标系坐标进行替换，变量定义如下：

ξ＝x-x_Li；η＝y-y_Li；ζ＝z

平衡点附近的线性化运动方程可描述为式(2)：

式(2)中，ρ²＝x²+y²+z²；c₂(μ)、c_n(μ)为仅与质量有关的常数，可表示为式(3)：

式(3)中，γ为平衡点与主天体S₂的距离，p_n为n阶Legendre多项式。对方程(2)进行分析，平衡点附近Halo轨道线性项可表示为式(4)：

式(4)中，ω₀分别为周期轨道运动的频率，κ为常数，φ为Halo轨道相位；α、β分别为轨道平面和垂直振幅；这里选择Halo轨道的垂直振幅为12000km，但为了满足Halo轨道的周期性，要求平衡点附近Halo轨道的非线性项可表示为式(5)：

式(5)中，θ＝ω₀t+φ，i+j＝N称为解的阶数。利用式(4)和(5)作为平衡点附近Halo轨道的初值，通过多重打靶方法可以得到Halo轨道的精确数值解，即三体系统下平衡点附近的Halo轨道。

步骤三：根据扰动变量设计Halo轨道保持的微分修正算法；根据微分修正算法得到Halo轨道精确数值解。

引入扰动变量，相对于参考轨道的变分为(δx,δy,δz)，则六维状态矢量为：

则变分方程的状态空间形式表示为方程(6)：

与时间有关的矩阵A(t)是由4个3×3的矩阵构成，即式(7)：

式(7)中B(t)是一个时变函数，因此A(t)也是一个时变函数，I₃是3×3的单位矩阵，C为常值矩阵，满足如下表达式：

方程(6)的解的形式为式(8)：

这里Φ(t,t₀)是状态转移矩阵。状态转移矩阵是从初始状态到t时刻状态的线性映射，通过状态转移矩阵可估算出初始状态变分对轨道递推的影响。状态转移矩阵必须满足矩阵微分方程(9)：

方程(9)的初始条件为Φ(t₀,t₀)＝I₆，这里I₆为6×6的单位阵。在t时刻的单位状态转移矩阵，通过初始条件的数值积分得到。因为状态转移矩阵是6×6，所以方程(9)的递推可以通过36个一阶微分方程的积分得到。B(t)的元素取决于参考轨道，因此求解A(t)，方程(9)需要同时积分非线性方程，从而得到状态转移矩阵Φ(t,t₀)。

根据微分修正算法求解Halo轨道精确数值解一般需要通过多次迭代过程来完成，归结如下：

根据当前末端状态与理想状态的偏差，通过状态转移矩阵对初始状态进行修正，使得最终的末端状态满足要求的约束条件。选取六维状态矢量寻找初始状态使周期轨道在t_f时刻满足约束条件末端状态偏差与初始状态偏差的关系为式(10)：

矩阵的偏导数等于在t_f时刻的状态转移矩阵，那么方程(8)可以简化为方程(11)：

末端状态的偏差量为方程(12)：

为目标终端状态，为实际终端状态，将方程(12)代入方程(11)，可以解出一个初始状态修正量的估计值从而得到三维周期Halo轨道的精确数值解。

步骤四：根据步骤三中的微分修正算法设计Halo轨道的在轨保持策略。

由于Halo轨道关于x-z平面对称，周期轨道必定与x-z平面垂直，因此在交叉处有满足式(13)：

利用垂直平面相交的特点进一步对问题进行简化，为了满足Halo轨道在半个周期后满足垂直交叉的条件则需选择y_f＝0作为积分终止的条件，在利用方程(11)迭代过程中需要调整初始状态使得末端的状态和接近于0，

由步骤三可知理想Halo轨道保持需要满足垂直交叉条件，即满足但在真实星历环境下，严格满足约束的Halo轨道是不存在的，只需保证Halo轨道的振幅尺寸满足规定的范围即可。因此对于Halo轨道的在轨保持修正，可以不严格满足约束条件针对修正约束条件的严格性提出两种不同的修正策略：宽松修正方法和严格修正方法。

宽松修正方法是指：轨道积分穿过x-z平面时x方向的速度为零，即严格修正方法是指：轨道积分穿过x-z平面时x方向和z方向的速度均为零，也即

当轨道穿过x-z平面时z方向的振幅超过任务约束的任务上下边界时执行一次严格修正方法，否则执行宽松修正方法。同时在每两个轨道周期内执行的4次修正中至少执行1次严格修正策略。

步骤五：在真实星历环境下，考虑测量误差和执行误差，采用步骤四中的在轨保持策略进行Halo轨道在轨保持。

轨道保持修正采用步骤三中的微分修正算法。在求解修正脉冲时，其积分初值为考虑测量误差的实际状态；而在轨道修正时，需要在实际状态进行脉冲修正时考虑修正的执行误差。这里星历模型选择DE405星历，测量误差位置精度1.6km(3σ)，速度精度1cm/s(3σ)，执行误差为推力方向误差为1°(3σ)，推力大小误差0.02m/s(3σ)。

在考虑测量误差和执行误差的条件下，按照轨道保持策略采用微分修正算法得到每次轨道修正的速度增量，按照速度增量进行轨道修正控制可以实现满足振幅约束的Halo轨道在轨保持，同时尽可能降低轨道保持所需的燃料消耗。选择在轨保持时间范围为2018年6月1日～2021年6月1日。在轨保持3年的三维轨道示意图如图4所示，图5～图7分别为三维轨道示意图在x-y，x-z和y-z平面的投影。

若采用先技术[2]中所述Halo轨道在轨保持控制方法，则在轨保持3年的三维轨道示意图如图8所示。本发明方法和先技术[2]中方法均以垂直振幅(z向振幅)12000km的Halo轨道为例，振幅上下边界约束范围为11500～12500km，分别采用相同在轨时间范围和测量执行误差，打靶100次得到在轨保持3年的结果如表1所示。对比图4和图8和表1中数据可以看出，本发明方法相比先技术[2]中方法所需速度增量略多，但可以有效控制Halo轨道振幅范围基本满足约束，即长期维持Halo轨道振幅不发散偏移。

表1本发明方法和先技术[2]方法结果对比

以上所述的具体描述，对发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例，用于解释本发明，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (2)

1.考虑振幅约束的一种Halo轨道在轨保持方法，其特征在于：通过在两个主天体和探测器构成的限制性三体模型下建立动力学方程，确定两个主天体和探测器构成的三体系统平衡点位置；确定大小天体和探测器构成的三体系统下平衡点附近的Halo轨道；根据扰动变量设计Halo轨道保持的微分修正算法；根据微分修正算法设计Halo轨道在轨保持方法。

2.考虑振幅约束的一种Halo轨道在轨保持方法，其特征在于：包括如下步骤：

步骤一：在两个主天体和探测器构成的限制性三体模型下建立动力学方程，确定两个主天体和探测器构成的三体系统平衡点位置；

限制性三体模型动力学系统中考虑质量可以忽略的探测器S₃在两个主天体S₁和S₂的引力共同作用下的运动，两个主天体S₁、S₂和探测器S₃三者质量关系为m₁＞m₂＞＞m₃；选择两个主天体系统的质心作为原点建立旋转坐标系，x轴方向由原点指向小天体，z轴为系统旋转的角速度方向，y轴与x、z轴构成右手坐标系；则探测器的运动描述为式(1)：

$\left\{\begin{array}{c}\stackrel{\·\·}{x}-2\stackrel{\·}{y}-x=-\frac{1-\mathrm{\μ}}{r_1^3}(x+\mathrm{\μ})-\frac{\mathrm{\μ}}{r_2^3}(x-1+\mathrm{\μ})\\ \stackrel{\·\·}{y}+2\stackrel{\·}{x}-y=-\frac{1-\mathrm{\μ}}{r_1^3}y-\frac{\mathrm{\μ}}{r_2^3}y\\ \stackrel{\·\·}{z}=-\frac{1-\mathrm{\μ}}{r_1^3}z-\frac{\mathrm{\μ}}{r_2^3}z\end{array}\right.---\left(1\right)$

其中μ＝m₂/(m₁+m₂)表示系统的质量系数，m₁为较大主天体S₁的质量，m₂为较小主天体S₂的质量，分别为运动探测器到两主天体的距离；

为获得旋转坐标系平衡点，令式(1)中速度和加速度等于0，共可获得五个平衡点(拉格朗日点)；所述五个平衡点均在x-y平面上，包括三个共线的平衡点和两个三角平衡点；在质心旋转系下三个共线平衡点的位置(x_Li,y_Li)分别为：

L1平衡点：L2平衡点：L3平衡点：

两个三角平衡点的位置分别为：

L4平衡点：L5平衡点：

步骤二：确定大主天体、小主天体和探测器构成的三体系统下平衡点附近的Halo轨道；

为获得更简洁的表达式，引入变量(ξ,η,ζ)，根据步骤一所得的平衡点的位置对旋转坐标系坐标进行替换，变量定义如下：

ξ＝x-x_Li；η＝y-y_Li；ζ＝z

平衡点附近的线性化运动方程可描述为式(2)：

$\left\{\begin{array}{c}\stackrel{\·\·}{\mathrm{\ξ}}-2\stackrel{\·}{\mathrm{\η}}-(1+2c_2)\mathrm{\ξ}=\frac{\∂}{\∂\mathrm{\ξ}}\underset{n\≥3}{\Σ}{c}_n\left(\mathrm{\μ}\right){\mathrm{\ρ}}^n{P}_n\left(\frac{\mathrm{\ξ}}{\mathrm{\ρ}}\right)\\ \stackrel{\·\·}{\mathrm{\η}}+2\stackrel{\·}{\mathrm{\ξ}}+(c_2-1)\mathrm{\η}=\frac{\∂}{\∂\mathrm{\η}}\underset{n\≥3}{\Σ}{c}_n\left(\mathrm{\μ}\right){\mathrm{\ρ}}^n{P}_n\left(\frac{\mathrm{\ξ}}{\mathrm{\ρ}}\right)\\ \stackrel{\·\·}{\mathrm{\ζ}}+c_2\mathrm{\ζ}=\frac{\∂}{\∂\mathrm{\ζ}}\underset{n\≥3}{\Σ}{c}_n\left(\mathrm{\μ}\right){\mathrm{\ρ}}^n{P}_n\left(\frac{\mathrm{\ξ}}{\mathrm{\ρ}}\right)\end{array}\right.---\left(2\right)$

式(2)中，ρ²＝x²+y²+z²；c₂(μ)、c_n(μ)为仅与质量有关的常数，可表示为式(3)：

$\left\{\begin{array}{c}{c}_2\left(\mathrm{\μ}\right)=\frac{1}{{\mathrm{\γ}}^3}\[\mathrm{\μ}+(1-\mathrm{\μ})\frac{{\mathrm{\γ}}^3}{{(1+\mathrm{\γ})}^3}\]\\ c_n\left(\mathrm{\μ}\right)=\frac{1}{{\mathrm{\γ}}^3}\[{(-1)}^n\mathrm{\μ}+{(-1)}^n(1-\mathrm{\μ}){\left(\frac{\mathrm{\γ}}{1+\mathrm{\γ}}\right)}^{n+1}\],n\≥3\end{array}\right.---\left(3\right)$

式(3)中，γ为平衡点与主天体S₂的距离，p_n为n阶Legendre多项式；对方程(2)进行分析，平衡点附近Halo轨道线性项可表示为式(4)：

$\left\{\begin{array}{c}\mathrm{\ξ}\left(t\right)=\mathrm{\α}cos({\mathrm{\ω}}_{0}t+\mathrm{\φ})\\ \mathrm{\η}\left(t\right)=\mathrm{\κ}\mathrm{\α}sin({\mathrm{\ω}}_{0}t+\mathrm{\φ})\\ \mathrm{\ζ}\left(t\right)=\mathrm{\β}\cos ({\mathrm{\ω}}_{0}t+\mathrm{\φ})\end{array}\right.---\left(4\right)$

式(4)中，ω₀分别为周期轨道运动的频率，κ为常数，φ为Halo轨道相位；α、β分别为轨道平面和垂直振幅；但为了满足Halo轨道的周期性，要求平衡点附近Halo轨道的非线性项可表示为式(5)：

$\left\{\begin{array}{c}\mathrm{\ξ}\left(t\right)=\underset{i,j}{\overset{\mathrm{\∞}}{\Σ}}\left(\underset{\left|k\right|\≤i+j}{\Σ}{\mathrm{\ξ}}_{ijk}cos\right(k\mathrm{\θ}\left)\right){\mathrm{\α}}^i{\mathrm{\β}}^j\\ \mathrm{\η}\left(t\right)=\underset{i,j}{\overset{\mathrm{\∞}}{\Σ}}\left(\underset{\left|k\right|\≤i+j}{\Σ}{\mathrm{\η}}_{ijk}sin\right(k\mathrm{\θ}\left)\right){\mathrm{\α}}^i{\mathrm{\β}}^j\\ \mathrm{\ζ}\left(t\right)=\underset{i,j}{\overset{\mathrm{\∞}}{\Σ}}\left(\underset{\left|k\right|\≤i+j}{\Σ}{\mathrm{\ζ}}_{ijk}cos\right(k\mathrm{\θ}\left)\right){\mathrm{\α}}^i{\mathrm{\β}}^j\end{array}\right.---\left(5\right)$

式(5)中，θ＝ω₀t+φ，i+j＝N称为解的阶数；利用式(4)和(5)作为平衡点附近Halo轨道的初值，通过多重打靶方法可以得到Halo轨道的精确数值解，即三体系统下平衡点附近的Halo轨道；

步骤三：根据扰动变量设计Halo轨道保持的微分修正算法；根据微分修正算法得到Halo轨道精确数值解；

引入扰动变量，相对于参考轨道的变分为(δx,δy,δz)，则六维状态矢量为：

$\mathrm{\δ}\stackrel{\‾}{x}\≡{\left[\begin{array}{cccccc}\mathrm{\δ}x& \mathrm{\δ}y& \mathrm{\δ}z& \mathrm{\δ}\stackrel{\·}{x}& \mathrm{\δ}\stackrel{\·}{y}& \mathrm{\δ}\stackrel{\·}{z}\end{array}\right]}^T$

则变分方程的状态空间形式表示为方程(6)：

$\mathrm{\δ}\stackrel{\·}{\stackrel{\‾}{x}}\left(t\right)=A\left(t\right)\mathrm{\δ}\stackrel{\‾}{x}\left(t\right)---\left(6\right)$

与时间有关的矩阵A(t)是由4个3×3的矩阵构成，即式(7)：

$A\left(t\right)=\left[\begin{array}{cc}0& I_3\\ B\left(t\right)& C\end{array}\right]---\left(7\right)$

式(7)中B(t)是一个时变函数，因此A(t)也是一个时变函数，I₃是3×3的单位矩阵，C为常值矩阵，满足如下表达式：

$C=\left[\begin{array}{ccc}0& 2& 0\\ -2& 0& 0\\ 0& 0& 0\end{array}\right]$

方程(6)的解的形式为式(8)：

$\mathrm{\δ}\stackrel{\‾}{x}\left(t\right)=\mathrm{\Φ}(t,t_0)\mathrm{\δ}\stackrel{\‾}{x}\left(t_0\right)---\left(8\right)$

这里Φ(t,t₀)是状态转移矩阵；状态转移矩阵是从初始状态到t时刻状态的线性映射，通过状态转移矩阵可估算出初始状态变分对轨道递推的影响；状态转移矩阵必须满足矩阵微分方程(9)：

$\stackrel{\·}{\mathrm{\Φ}}(t,t_0)=A\left(t\right)\mathrm{\Φ}(t,t_0)---\left(9\right)$

方程(9)的初始条件为Φ(t₀,t₀)＝I₆，这里I₆为6×6的单位阵；在t时刻的单位状态转移矩阵，通过初始条件的数值积分得到；因为状态转移矩阵是6×6，所以方程(9)的递推可以通过36个一阶微分方程的积分得到；B(t)的元素取决于参考轨道，因此求解A(t)，方程(9)需要同时积分非线性方程，从而得到状态转移矩阵Φ(t,t₀)；

根据微分修正算法求解Halo轨道精确数值解一般需要通过多次迭代过程来完成，归结如下：

根据当前末端状态与理想状态的偏差，通过状态转移矩阵对初始状态进行修正，使得最终的末端状态满足要求的约束条件；选取六维状态矢量寻找初始状态使周期轨道在t_f时刻满足约束条件末端状态偏差与初始状态偏差的关系为式(10)：

$\mathrm{\δ}\stackrel{\‾}{x}\left(t_f\right)=\frac{\∂\stackrel{\‾}{x}\left(t_f\right)}{\∂\stackrel{\‾}{x}\left(t_0\right)}\mathrm{\δ}\stackrel{\‾}{x}\left(t_0\right)+\stackrel{\·}{\stackrel{\‾}{x}}\left(t_f\right)\mathrm{\δ}(t_f-t_0)---\left(10\right)$

矩阵的偏导数等于在t_f时刻的状态转移矩阵，那么方程(8)可以简化为方程(11)：

$\mathrm{\δ}x\left(t_f\right)=\mathrm{\Φ}(t_f,t_0)\mathrm{\δ}\stackrel{\‾}{x}\left(t_0\right)+\stackrel{\·}{\stackrel{\‾}{x}}\left(t_f\right)\mathrm{\δ}(t_f-t_0)---\left(11\right)$

末端状态的偏差量为方程(12)：

$\mathrm{\δ}\stackrel{\‾}{x}\left(t_f\right)=\stackrel{\‾}{x}{\left(t_f\right)}_{des}-\stackrel{\‾}{x}\left(t_f\right)---\left(12\right)$

为目标终端状态，为实际终端状态，将方程(12)代入方程(11)，可以解出一个初始状态修正量的估计值

由于Halo轨道关于x-z平面对称，周期轨道必定与x-z平面垂直，因此在交叉处有满足式(13)：

$y=\stackrel{\·}{x}=\stackrel{\·}{z}=0---\left(13\right)$

利用垂直平面相交的特点进一步对问题进行简化，为了满足Halo轨道在半个周期后满足垂直交叉的条件则需选择y_f＝0作为积分终止的条件，在利用方程(11)迭代过程中需要调整初始状态使得末端的状态和接近于0，从而得到三维周期Halo轨道的精确数值解；

步骤四：根据步骤三中的微分修正算法设计Halo轨道的在轨保持策略；

由步骤三可知理想Halo轨道保持需要满足垂直交叉条件，即满足但在真实星历环境下，严格满足约束的Halo轨道是不存在的，只需保证Halo轨道的振幅尺寸满足规定的范围即可；因此对于Halo轨道的在轨保持修正，可以不严格满足约束条件针对修正约束条件的严格性提出两种不同的修正策略：宽松修正方法和严格修正方法；

宽松修正方法是指：轨道积分穿过x-z平面时x方向的速度为零，即严格修正方法是指：轨道积分穿过x-z平面时x方向和z方向的速度均为零，也即

当轨道穿过x-z平面时z方向的振幅超过任务约束的任务上下边界时执行一次严格修正方法，否则执行宽松修正方法；同时在每两个轨道周期内执行的4次修正中至少执行1次严格修正策略；

步骤五：在真实星历环境下，考虑测量误差和执行误差，采用步骤四中的在轨保持策略进行Halo轨道在轨保持；

轨道保持修正采用步骤三中的微分修正算法；在求解修正脉冲时，积分初值为考虑测量误差的实际状态；而在轨道修正时，需要在实际状态进行脉冲修正时考虑修正的执行误差；在考虑测量误差和执行误差的条件下，按照轨道保持策略采用微分修正算法得到每次轨道修正的速度增量，按照速度增量进行轨道修正控制可以实现满足振幅约束的Halo轨道在轨保持，同时尽可能降低轨道保持所需的燃料消耗。