圆型限制性三体问题相对运动解析研究

周 敬，胡 军，张 斌

(1. 北京控制工程研究所， 北京 100190; 2. 空间智能控制技术国防科技重点实验室， 北京 100190)

0 引 言

近年来，随着航天技术的快速发展，深空探测受到世界各航天大国与组织越来越多的重视。按照我国航天事业发展的指南，未来航天的三个重点方向之一就是深空探测[1]。因此，与深空探测密切相关的各项基础和应用研究便成为了我国航天事业中的一个热点研究方向。

不同于经典二体问题下的近地空间航天活动，深空探测所涉及的运动模型更多的是三体问题。由于动力学本质上的区别，三体问题相对二体问题更加复杂。在三体问题的研究中，最简单、应用最广泛的模型是圆型限制性三体问题(Circular restricted three-body problem，CRTBP)。在CRTBP中存在5个动平衡点，也称为平动点(Libration points)或拉格朗日点(Lagrangian points)。这5个平动点相对主天体的位置是固定的，其中L1,L2和L3与两主天体共线(L1位于主天体之间的连线上，L2和L3位于延长线上)，称为共线平动点；L4和L5均与两个主天体构成等边三角形，称为三角平动点。平动点(尤其是共线平动点)及其附近的周期/拟周期轨道所具有的独特位置优势使其在深空探测任务设计中占据着重要的地位，已经成为宇宙观测、天文研究的理想场所和星际高速公路(Interplanetary superhighway，IPS)的门户站。

为了实现重返月球，美国航天局和其他太空研究机构已经考虑在地-月L1点建立空间站，以支持载人登月、载人火星探测甚至宇宙航行。与近地轨道空间站一样，平动点空间站的建设和运行同样离不开交会对接(Rendezvous and docking，RVD)技术的支持，而RVD技术的基础则是相对运动。此外，与近地轨道航天任务发展趋势一样，深空探测任务发展趋势也必将由单颗航天器独立工作向多颗航天器协同工作转变，如TPF、Darwin、MAXIM计划等，而多颗航天器协同工作所采用的编队飞行技术更需要相对运动研究作为基石。因此，对三体问题平动点附近相对运动的研究兼具理论研究和工程应用价值。

目前，涉及三体问题相对运动的研究工作主要有：侯锡云等[2]系统地阐述了限制性三体问题中共线平动点的动力学特征，给出了这类平动点附近的周期轨道和拟周期轨道的计算方法；徐明[3]详细介绍了平动点轨道的发展历史，并深入剖析了平动点附近的相空间结构；文献[4]通过建立改进型庞加莱截面图，同时结合状态转移矩阵和打靶法，建立了一种计算三体问题周期轨道的新方法；Richardson等[5]将推导Halo轨道三阶解析解的思想应用到编队飞行，得到了非线性相对运动方程的三阶展开式；Gómez等[6]将其关于Lissajous轨道的Fourier算法等研究成果推广到相对运动动力学，得到了平动点相对运动的高阶展开式；文献[7-9]基于Richardson三阶解析解研究了Halo轨道上编队飞行航天器的相对运动构型特性；Howell等[10]运用瞄准法获得了期望的日-地平动点附近航天器的非自然编队构型；文献[11]基于单值矩阵和Fourier级数，将编队设计问题从微分方程转化为代数方程，分别设计了沿Halo轨道的长期和短期自然松散编队构型；Luo等[12]提出了一种新型的、基于遍历Poincaré映射的数值搜索方法，获得了CRTBP拟Halo轨道附近的自然编队飞行构型；Luquett等[13]以全非线性形式给出了限制性三体问题下的自然相对运动动力学模型，然后将其线性化并转化为状态空间的形式，获得了相对运动的近似解析解；文献[14]应用Lindstedt-Poincaré方法和多项式展开的方法，构造了限制性三体问题平动点附近周期相对运动构型的解析解；Cont等[15]通过线性化的方法对限制性三体问题平动点附近的相对运动和交会问题进行解析分析，并与数值仿真结果进行了对比验证；Ferrari等[16]通过分析初始构型和关键参数的方法，研究了圆型和椭圆型限制性三体问题下共线平动点周期轨道附近无控编队飞行的演化过程；黄海滨等[17]运用Legendre伪谱法和协同进化粒子群(CPSO)方法研究了深空环境下卫星编队飞行队形重构实时重规划问题；Arnot等[18]在CRTBP模型下，应用最小位置反馈方法构造了新的稳定的、调频的相对运动轨道；Case[19]运用线性化相对运动动力学和非线性微分修正方法实现了CRTBP下平动点轨道的交会；Prioroc等[20]利用极值搜索方法对时延反馈机构和PD控制器的增益进行优化，实现了小推力下的日-地CRTBP平动点附近的编队飞行；Jung等[21]运用开关-Hamiltonian保结构控制方法实现了地-月L2点Halo轨道上的编队飞行控制；文献[22]设计了基于输出调节理论的位置保持控制器，实现了CRTBP平动点附近航天器的位置保持和编队飞行控制；Prioroc等[23]运用非线性控制的方法研究了Halo轨道编队飞行的相对运动控制问题；张楷田等[24]基于CRTBP模型推导了航天器编队飞行的非线性动力学方程，然后应用线性自抗扰控制技术实现了编队飞行控制；姜春生等[25]结合扩张观测器和自抗扰控制方法，实现了对日-地平动点附近的航天器编队飞行的控制；Yin等[26]运用自抗扰控制方法解决了考虑不确定性和摄动情况下的日-地CRTBP下L2平动点附近的多航天器编队飞行控制问题；Wang等[27]结合反馈控制和分布式自适应同步方法，实现了速度不可测和质量不确定情况下平动点轨道附近的多航天器编队飞行；王峰等[28]建立日-地L2平动点编队飞行相对运动QLPV动力学模型，并提出一种改进的多项式特征结构配置方法实现日-地L2平动点编队飞行高精度相对位置保持；文献[29]使用多重打靶修正方法改进了中心航天器的拟Halo参考轨道，并使用线性化的变分方程来进行位置保持控制，最终运用Hamiltonian保结构控制方法实现了L1点拟周期轨道附近编队飞行的稳定控制。

从上述研究中可以发现，有关三体问题下的相对运动研究多以编队飞行为研究背景，通过将轨道绝对运动的近似解析解直接作差，得到相对运动的近似解析解；或者将相对运动动力学方程线性化，然后运用线性系统理论对其求解，获得相对运动的近似解析解；亦或者将相对运动问题具化为实现某一具体任务而引申出的控制问题加以解决。这三种方式在多数情况下可以满足任务需求，但不便于进一步获得三体问题下相对运动的本质规律。基于此，本文借鉴Floquet理论的研究思想，在对单值矩阵进行分析的基础上，从更本质的角度出发，对三体问题中的共线平动点附近周期/拟周期轨道下的相对运动进行研究。

1 圆型限制性三体问题

限制性三体问题研究的是一个小天体或航天器在两个大天体(主天体)引力场中的运动情况，其中小天体或航天器的质量相比于两主天体可忽略不计，不会对两主天体之间的相互运动产生影响。当主天体围绕公共质心作圆轨道运动时，即为圆型限制性三体问题(CRTBP)。为方便起见，本文以地-月系统CRTBP中的L1点作为研究对象开展相对运动研究。本文研究方法同样适用于其他三体系统共线平动点附近的相对运动。

1.1 动力学模型

为方便描述航天器运动和简化计算，在CRTBP中通常需要对相关物理量进行无量纲化处理，并以时间作为独立变量。相应的质量[M]、长度[L]和时间[T]的归一化单位取为

(1)

式中:m1，m2为两主天体质量；L12为两主天体之间的距离；G为万有引力常数。

为了更直观、更清晰地描述航天器在CRTBP中的运动，通常选择在会合坐标系(也称为质心旋转坐标系)下进行研究。如图1所示，会合坐标系的原点位于两主天体的公共质心，x轴由较大主天体指向较小主天体，z轴与主天体系统角动量方向平行，y轴满足右手坐标系定则。

图1 地-月系统CRTBP下的会合坐标系和L1会合坐标系

(2)

式中：Ω表示伪势能函数。

(3)

式中：X1表示在会合坐标系下L1点距离主天体公共质心的距离，γ1为L1点与最近主天体之间的距离。

因此，L1合坐标系下航天器运动的动力学方程为

(4)

1.2 近似解析解

在会合坐标系下，假设[δX,δY,δZ]T表示航天器相对L1平动点的位置矢量。根据文献[25]，共线平动点附近航天器运动的近似解析解为

(5)

式中：Ai(i=1,…,6)是由初始条件确定的积分常数，并在本文中作为Lissajous轨道和Lyapunov轨道的轨道参数，参数k1,k2以及λ1,λ3,λ5的详细计算公式可以参考文献[30]。

由于Halo轨道是共线平动点附近一类特殊的三维周期轨道，且常作为深空探测航天器的工作轨道，具有更加重要的研究和应用价值。因此，Richardson[31]基于Lindstedt-Poincaré方法获得了Halo轨道的三阶近似解析解，如式(6)所示，相关参数的含义及计算过程可以参考文献[31]。

(6)

此外，关于Halo轨道三阶近似解析解的改进可以参考文献[32]。

1.3 微分修正

由于CRTBP具有强非线性、不稳定性和“混沌”性，第1.2节获得的近似解析解在真实动力学模型下积分一段时间后会很快发散。为了获得精确、收敛的周期轨道，可以采用微分修正方法。有关微分修正方法的详细介绍可以参考文献[1]。微分修正方法主要利用周期/拟周期轨道的对称性来修正近似解析解提供的轨道初值，使修正后的初值进行数值积分后生成的周期轨道可以自然维持多圈而不发散。下面分别介绍微分修正方法在Halo轨道、Lyapunov轨道和Lissajous轨道上的具体应用过程。

1.3.1Halo轨道

由于Halo轨道关于XZ平面对称，根据参考文献[1]和近似解析解(6)，微分修正的初始状态x0设置为

(7)

(8)

式中：Φ(i,j)表示矩阵Φ中第i行，第j列元素。式(8)即为Halo轨道初值的微分修正公式。详细修正过程可以参考文献[1]。

1.3.2Lyapunov轨道

由于Lyapunov轨道是一种XY平面轨道且关于x轴对称，根据参考文献[1]和近似解析解(5)，微分修正的初始状态x0设置为

(9)

(10)

式(10)即为Lyapunov轨道初值的微分修正公式。详细修正过程可以参考文献[1]。

1.3.3Lissajous轨道

Lissajous轨道是一种拟周期轨道，而非严格意义上的周期轨道，因此其微分修正过程较Halo轨道和Lyapunov轨道更为复杂。Howell等[33]提出使用两层微分修正方法来求解CRTBP下的Lissajous轨道，考虑到两层微分修正方法的复杂性，本文提出一种较简单的修正方法。鉴于Lissajous轨道在XY平面内的运动与Z轴的运动是解耦的，且两部分运动的周期也不相同，因此可以在XY平面和Z轴上分别、同时进行微分修正，最终获得精确的Lissajous轨道初值。

首先考虑XY平面内的运动，Lissajous轨道在XY平面内的运动轨道实际上为Lyapunov轨道，因此在XY平面上的微分修正可以沿用Lyapunov轨道的修正方法。

根据近似解析解(5)，微分修正的初始状态x0设置为

(11)

对应的第一个期望状态设置为

(12)

经过时间t1后的真实状态x1为

(13)

则真实状态x1和期望状态xd1之差δx1为

(14)

(15)

然后，考虑Z轴上的运动，对应的第二个期望状态设置为

(16)

经过时间t2后的真实状态x2为

(17)

则真实状态x2和期望状态xd2之差δx2为

(18)

对式(18)进行反解，得到修正量δz0和δt2的计算公式

(19)

式(15)和式(19)即为Lissajous轨道初值的微分修正公式。上述两部分微分修正同时进行，便可以获得精确的Lissajous轨道初值。

2 相对运动

相对运动通常涉及两航天器，即目标航天器和追踪航天器。假设目标航天器的运行轨道为地-月L1点附近的一条周期/拟周期轨道，追踪航天器运行在目标航天器附近。故本文相对运动特指，以某个Halo/Lyapunov/Lissajous轨道为参考轨道，研究航天器在此参考轨道附近相对此参考轨道的运动。

图2 相对运动示意图

本文研究相对运动的目的是获得与近地轨道CW方程类似的CRTBP下相对运动的近似解析解。下面以Halo轨道作为目标航天器的运行轨道为例，介绍相对运动的研究方法。

在CRTBP中，单值矩阵(Monodromy Matrix)M是指一个轨道周期T对应的周期轨道的状态转移矩阵(State transition matrix,STM)Φ(t0+t,x0)，即

M=Φ(t0+T,x0)

(20)

虽然单值矩阵没有精确的解析表达式，但可以通过数值积分的方式获得精确的数值解，详细的计算方法可以参考文献[25]。

单值矩阵共存在6个特征值，表示为(λ1,λ2,λ3,λ4,λ5,λ6)，可以划分为三对，即

(21)

第一对特征值(λ1,λ2)是互为倒数的正实数，与参考轨道的双曲特性有关。λ1>1是绝对值最大的特征值，其对应的特征向量代表着参考轨道的发散方向，因此λ1及其对应的特征向量在一定程度上代表着参考轨道的不稳定特性。与之相反，由于λ2=1/λ1<1，故λ2及其对应的特征向量在一定程度上表示参考轨道的稳定特性。

第二对特征值(λ3,λ4)的数值均为1，在一定程度上代表着参考轨道的既不发散也不收敛特性，即参考轨道的中性特性。

第三对特征值(λ5,λ6)为模值为1的共轭复数，在一定程度上代表着参考轨道的中心流形，即参考轨道的周期特性。

综上所述，单值矩阵的六个特征值和对应的特征向量提供了一种具有几何意义的表示，即不同的特征值和对应的特征向量可以代表相对运动中具有不同特性的分量模态，如不稳定模态、稳定模态、中性模态以及周期模态，这也是Floquet理论[34]的核心思想。因此，通过对单值矩阵的分析，可以建立由一组独立基向量所形成的一个局部坐标系统，并以此局部坐标系统来描述相对运动。借鉴Floquet理论的思想，这组基向量可以由单值矩阵的特征向量形成，而这组基向量也就是本文提出的相对运动模态。

假设相对运动模态矩阵E(t)由6个相对运动模态组成

E(t)=[e1(t),e2(t),e3(t),e4(t),e5(t),e6(t)]

(22)

式中：ei(t),i=1,2,3,…,6即为6个相对运动模态，分别与单值矩阵的6个特征值和特征向量相对应。

通过对单值矩阵的分析，可以令Halo轨道上当前状态对应的单值矩阵的特征向量组成相对运动模态矩阵E(t)的初始值E0，即

E0=[v1,v2,v3,v4,real(v5),imag(v6)]

(23)

式中：vi(i=1,…,6)即为单值矩阵的6个特征向量，可以通过求解单值矩阵获得，real(v5)表示特征向量的实部，imag(v6)表示特征向量的虚部。不难看出rank(E0)=6，即这组基向量是线性无关的，可以用来描述6维空间中的任意向量，即可以用来描述任意初始状态的相对运动。

根据Halo轨道的性质[30]以及状态转移矩阵性质，相对运动模态矩阵满足如下关系

E(t)=Φ(t)E0

(24)

在获得相对运动模态变化规律E(t)的基础上，便可以将相对运动看做上述6个相对运动模态的线性组合。假设追踪航天器和目标航天器的运动状态分别记为x′和x，则两者之间的差值即为航天器间的相对运动δx，

(25)

将相对运动δx表示成相对运动模态的线性组合形式，即

(26)

式中：C(t)为各个模态在相对运动中所占的权值矩阵，不难发现其为一个对角矩阵，即

C(t)=diag(c1(t),c2(t),c3(t),c4(t),

c5(t),c6(t))

(27)

而相对运动的初始状态δx0可以表示为

δx0=C0E0

(28)

根据文献[35]，由于相对运动的尺度相对参考轨道尺度来说为一小量，因此可以将追踪航天器相对于目标航天器的相对运动δx看做目标航天器状态x的一种摄动，则根据状态转移矩阵的性质，有

δx(t)=Φ(t)δx0

(29)

将式(29)与式(24)和式(28)联立，有

(30)

从式(30)可以得到

C(t)=C0

(31)

即权值矩阵C(t)为一常值矩阵。

这样，便获得了Halo轨道为参考轨道的相对运动近似解析解

δx(t)=C0Φ(t)E0

(32)

由于任意时刻的状态转移矩阵Φ(t)可以事先获得，即作为已知量，依据此近似解析解，便可以获得任意时刻的相对运动状态。

由于单值矩阵M是非线性系统周期解的一种共有的特性，而Lyapunov轨道同Halo轨道一样均为非线性CRTBP动力学方程的周期解；此外，虽然Lissajous轨道仅为拟周期轨道，但实际上其XY平面内的运动周期与Z轴上的运动周期相差并不大，在一定精度条件下也可以视为周期轨道。因此，上述Halo轨道下的相对运动分析方法同样适用于Lyapunov轨道和Lissajous轨道下的相对运动。

根据上述分析，本文提出的CRTBP共线平动点附近相对运动近似解析解的适用条件如下：

1)目标航天器运行的参考轨道为Halo轨道、Lyapunov轨道时的相对运动，以及精度要求不高的情况下，也适用于Lissajous轨道为参考轨道的相对运动。

2)相对运动相比参考轨道的尺度较小，即相对运动可以作为参考轨道的摄动来处理。

对比参考文献[11]，本文提出的相对运动研究方法不仅适用于Halo轨道，也适用于Lyapunov轨道和Lissajous轨道，是一种在Floquet理论基础上改进的、通用的、基于单值矩阵的CRTBP共线平动点附近周期/拟周期轨道下相对运动的解析研究方法。

3 仿真校验

为校验本文所提出的基于单值矩阵的CRTBP共线平动点附近周期/拟周期轨道下的相对运动解析研究方法的正确性，本节分别以地-月系统L1点附近的Halo轨道、Lissajous轨道和Lyapunov轨道作为目标航天器的参考轨道，研究追踪航天器与目标航天器之间的相对运动。

3.1 Halo轨道

假设目标航天器和追踪航天器均运行在各自的Halo轨道上，相应的Z轴振幅分别记Az1和Az2，相对运动时间记为tTOF，xT0为目标航天器的初始状态，xC0为追踪航天器的初始状态，以上各参数在L1会合坐标系下的无量纲单位(Dimensionless unit, DU)取值如表1所示。

Halo轨道相对运动的6个模态的初值以及一个轨道周期内的变化趋势分别如表2和图3所示。从仿真结果可以看出，模态1呈发散趋势，证明了其为不稳定模态；模态2呈收敛趋势，证明了其为稳定模态；模态3和模态4呈既不发散也不收敛趋势，证明了其为中性模态；模态5和模态6呈周期变化趋势，证明了其为周期模态。

表1 Halo轨道相对运动相关参数

表2 Halo轨道相对运动6个模态初始值

图3 Halo轨道相对运动6个模态位置空间的变化规律

图4表示真实Halo轨道相对运动与反演Halo轨道相对运动的对比。真实相对运动为追踪航天器的轨道数据与目标航天器的轨道数据之差；反演相对运动是由相对运动近似解析解式(32)生成。从图4可以看出，反演相对运动与真实相对运动相一致，虽然由于线性化的原因会不可避免地出现一定的误差，但依然可以说明本方法的有效性。

图4 Halo轨道相对运动仿真结果

3.2 Lissajous轨道

假设目标航天器和追踪航天器均运行在各自的Lissajous轨道上，相应的轨道参数[A1,A2,A3,A4,A5,A6]分别记A1和A2，相对运动时间记为tTOF，xT0为目标航天器的初始状态，xC0为追踪航天器的初始状态，以上各参数在L1会合坐标系下的无量纲单位取值如表3所示。

表3 Lissajous轨道相对运动相关参数

Lissajous轨道相对运动的6个模态的初值以及一个周期内的变化趋势分别如表4和图6所示。从仿真结果可以看出，6个相对运动模态同样均与其自身的性质相符合。

图5表示真实Lissajous轨道相对运动与反演Lissajous轨道相对运动的对比。从图5可以看出，反演相对运动与真实相对运动相一致，误差很小，说明了本方法的有效性。

图5 Lissajous轨道相对运动仿真结果

表4 Lissajous轨道相对运动6个模态初始值

3.3 Lyapunov轨道

假设目标航天器和追踪航天器均运行在各自的Lyapunov轨道上，相应的轨道参数[A1,A2,A3,A4,A5,A6]分别记A1和A2，相对运动时间记为tTOF，xT0为目标航天器的初始状态，xC0为追踪航天器的初始状态，以上各参数在L1会合坐标系下的无量纲单位取值如表5所示。

Lyapunov轨道相对运动的6个模态的初值以及一个周期内的变化趋势分别如表6和图7所示。从仿真结果可以看出，6个相对运动模态同样均与其自身的性质相符合。

图8表示真实Lyapunov轨道相对运动与反演Lyapunov轨道相对运动的对比。从图8可以看出，反演相对运动与真实相对运动相一致，误差很小，说明了本方法的有效性。

表5 Lyapunov轨道相对运动相关参数

表6 Lyapunov轨道相对运动6个模态初始值

图7 Lyapunov轨道相对运动6个模态位置子空间的变化规律

图8 Lyapunov轨道相对运动仿真结果

4 结 论

本文针对圆型限制性三体问题共线平动点附近周期/拟周期轨道下的相对运动问题，给出了一种通用的、基于单值矩阵分析的相对运动近似解析解。该方法通过对三体问题下周期/拟周期轨道单值矩阵的分析，依据单值矩阵的特征值和对应的特征向量建立了6个相对运动模态，具体分为不稳定模态、稳定模态、中性模态以及周期模态，然后将航天器间的相对运动表示为上述模态的线性组合，最终获得了相对运动的近似解析解。最后，以地-月系统L1点为研究对象，对参考轨道分别为Halo轨道、Lissajous轨道和Lyapunov轨道时的相对运动模态和相对运动进行了仿真分析，说明了方法的有效性。