同步双小行星系统共振轨道设计

杨雅迪，陈 奇，李翔宇，乔 栋(1. 北京理工大学宇航学院，北京 100081；2. 深空探测自主导航与控制工信部重点实验室，北京 100081)

0 引 言

自20世纪90年代以来，小行星探测成为深空探测的热点。在小行星探测任务中，具有独特动力学特性的双小行星系统引起了科学家的关注[1]。通过对双小行星系统展开实地探测与采样，可获取小行星内部构造、物质组成和动力学演化等方面的重要信息，能得到对单颗小行星探测时无法获得的独特科学回报。

目前各国提出了多项针对双小行星系统的探测计划。美国NASA与欧洲ESA联合提出了小行星撞击和偏转评估(AIDA)任务[2]，计划对双小行星65803 Didymos进行撞击，评估验证近地小行星防御的可行性；美国JPL同时提出双小行星定位探测(BASiX)任务[3]，以双小行星65803 Didymos为目标，研究微重力环境中瓦砾堆小行星的起源及演化。在双小行星系统的探测任务中，复杂探测任务的开展首先需要获得系统动力学环境相关的先验信息，而系统完整的引力场模型必须在探测器靠近探测目标，并展开长时间观测与绕飞后才能获得[4-5]。同时通过环绕探测，也将为着陆点选择和着陆轨迹设计提供参考[6-7]。因此选择合适的环绕探测轨道对探测任务的成功开展具有举足轻重的作用。为了满足小行星附近尤其是双小行星系统中的探测任务需求，需要多样性丰富、燃料消耗少且可保持长期稳定的探测轨道。共振轨道可以较好地满足这几大特性，因此研究小行星附近的共振轨道对未来探测任务具有较大的参考价值。

在对共振轨道的研究中，2013和2014年，Vaquero和Howell[8-9]对地月系统共振轨道的动力学结构进行了深入细致的分析，借助庞加莱截面设计了二维及三维共振轨道，并提出利用共振轨道实现低能量转移的轨道设计方法。2014年，Antoniadou和Voyatzis[10]研究了三体系统下地外行星系统平面及空间中共振比为4∶3,3∶2,5∶2,3∶1及4∶1的共振轨道，研究发现稳定共振轨道普遍具有较大的偏心率。同年，于洋[11]研究了单颗小行星赤道面附近的1∶1共振轨道，从能量角度揭示了赤道面附近的1∶1共振的动力学本质，确定了形成1∶1共振抛射轨道的临界条件。2016年，Anderson等[12]基于圆型限制性三体模型，利用网格搜索方法对木卫系统中不稳定共振轨道进行了搜索，并对不同共振比轨道进行了特性分析。研究发现在稳定共振轨道的基础上，可利用连续法计算得到不稳定的共振轨道。此后，Feng等[13]对接触式双小行星附近动力学环境进行了分析，发现了三类不同共振比的共振轨道。

此前对共振轨道的研究，主要集中在地月和木卫系统[14]等质量比较大的系统。双小行星系统具有质量比较小，形状不规则和相对距离接近等特点，导致其附近动力学特性复杂。而针对双小行星系统内的共振轨道尚未得到系统性研究。因此本文将以双小行星系统探测任务为背景，给出一种具有普适性的双小行星系统共振轨道设计方法，对不同共振比的平面及三维共振轨道进行设计，并分析共振轨道随轨道周期及能量的演化规律。本文的研究可以为未来双小行星系统探测轨道设计提供参考。

1 双小行星系统模型及共振轨道定义

1.1 双椭球体模型

本文采用双椭球体模型表征同步双小行星系统，建立质心旋转坐标系O-xyz如图1所示。其中A1为双体小行星系统中质量较大的小行星(主星)，A2为质量较小的小行星(从星)，B为航天器。x轴沿A1,A2之间的连线由A1指向A2，z轴平行于轨道角动量矢量方向，y轴与x轴、z轴构成右手坐标系。

图1 双小行星系统椭球体-椭球体模型Fig.1 The double ellipsoid model of binary asteroid

假设A1与A2同步自旋，且惯性主轴均位于x轴上。它们的半主轴分别为a1,b1,c1及a2,b2,c2，且满足0≤c1≤b1≤a1，0≤c2≤b2≤a2。记A1和A2的质量分别为M1,M2,m(m≪M1,M2)，则系统的质量比为：

(1)

(2)

Rb=[x,y,z]T为航天器相对于系统质心的位置矢量。

(3)

则探测器在质心旋转系下的运动方程可表示为：

(4)

Us=(1-μ)U1(ρ-R1)+μU2(ρ-R2)

(5)

其中，U1,U2分别代表双小行星系统中两颗均质三轴椭球体的外部势能，

(6)

(7)

式中:λ1,λ2分别为φ(ρ-R1,λ1)与φ(ρ-R2,λ2)的根，v为积分变量。

(8)

(9)

(10)

(11)

式(4)表示成分量形式为：

(12)

式中：Ea,Eb,Ec为代表双椭球质量分布的椭圆积分，具体求解过程可参考文献[15]。

1.2 共振轨道初值选取

首先忽略A2的质量，基于二体模型求解共振轨道，作为双小行星系统内共振轨道的初值，假设探测器B绕A1飞行p圈所需的时间与A2绕A1飞行q圈所需的时间相同，p和q均为正整数，B与A2的轨道周期比为：

(13)

其中，np,nq分别为B与A2的平均角速度，

(14)

Tp,Tq分别为B与A2的周期。μ1=GM1为A1的引力常数。根据此定义，在惯性系中，探测器初始状态可由式(15)～(18)确定：

(15)

p1=a(1-e2)

(16)

(17)

(18)

其中，p1为半正焦弦，r0为初始距离大小，v0为初始速度。θ为真近角，在近心点处，θ=0°，在远心点处，θ=180°。探测器B的周期Tp可以由共振比p:q,A2周期Tq和式(13)计算得到：

(19)

共振轨道的确定还需要知道偏心率e的大小，选择近心点或远心点作为探测器的初始状态，可更方便地得到平面共振轨道。因此本文将通过选取近心点位置rp(或远心点位置ra)和已求得的半长轴a来求解偏心率e。

若初始状态完全确定，即可解析求得共振轨道在惯性坐标系中的解。假设双小行星中主从星的半径分别为r1=12.38 km和r2=3.60 km，相对距离R=33 km，密度ρ=2.91 g/cm3，同步旋转周期为2.917 h。初始时刻为A1,A2和B三者共线的时刻，可求得仅考虑主星质量下的3∶4共振轨道，选择双星相对距离作为单位长度，在旋转坐标系下共振轨道如图2所示。

图2 二体系统中3∶4共振轨道Fig.2 3∶4 resonant orbit in binary asteroid system based on two-body model

2 双小行星系统中共振轨道族设计方法

将第1节二体模型中求得的共振轨道作为初值，利用改进并行打靶法修正可得双小行星系统下的共振轨道。修正分两步进行，第一次修正得到忽略小行星形状，仅考虑从星质量的共振轨道，并将它的初始状态作为第二次修正的初值，从而得到双椭球模型下具有周期性的共振轨道。然后，将该共振轨道作为初值，利用改进的连续法延拓得到一组共振比相同的共振轨道族。最后，利用动力学分岔理论，沿共振轨道的特征方向延拓生成新的空间共振轨道族。

2.1 改进并行打靶法

由于双小行星系统内动力学环境的复杂性，共振轨道对初值敏感，需要对二体模型下的轨道初值进行修正。本文采用并行打靶法进行轨道设计，该方法相对传统的打靶法而言，有更好的稳定性和更大的收敛域；同时相对于常用的差分方法，要求解的非线性方程组的规模要小得多、计算速度快。

图3 改进的并行打靶法示意图Fig.3 Improved parallel shooting method

(20)

引入松弛变量β，则约束方程可表示为

(21)

终端点的约束向量为

(22)

为了得到周期性的轨道，打靶过程中第一个节点一定在超平面上，因此有y1=yh=0。

将积分时间记作T，则改进的多次打靶法中的自由变量X为

(23)

X由6n+2个元素构成。约束函数F(X)为

F(X)由6n+1个分量组成。则Jacobi矩阵DF(X)为(6n+1)×(6n+2)维矩阵，具体表达式为：

DF(X)=

(24)

其中Φ(t2,t1)表示从时刻t1至t2的状态转移矩阵，

A=diag(-1, -1, -1, -1, 0, 1)

(25)

(26)

(27)

利用改进的并行打靶法，将初始轨道分成若干段，首先考虑双小行星系统中的子星的质量，将系统的质量比由0变为真实值进行第一次修正，图4给出不同质量比下双小行星系统中的3∶4共振轨道的演化情况。

图4 μ变化下共振轨道演化Fig.4 Evolution of resonance orbits under varying μ

系统质量比μ从0逐渐增大为真实值0.0245，图4(a)给出了不同系统质量比下的共振轨道演化情况，发现随着质量比μ的增大，共振轨道所对应的Jacobi常数增大，同时共振轨道与主星间距离逐渐增大，与子星的距离减小。μ=0.0245时球体-球体模型中的最终修正轨道如图4(b)所示。

进一步考虑小行星形状变化，可得双椭球模型下的共振轨道。假设双小行星的形状参数为

a1=7.25 km,b1=5.90 km,c1=5.55 km

a2=1.90 km,b2=1.75 km,c2=1.75 km

图5(a)给出形状变化下双小行星系统中的3∶4共振轨道的演化情况，随着主从星形状连续变化，3∶4共振轨道所对应的Jacobi常数值不断增大。轨道同时远离主星和子星，但偏移距离小于质量比变化造成的改变。利用改进并行打靶法和两步修正策略，可以实现双椭球模型下精确共振轨道的设计，最终的共振轨道如图5(b)所示。

图5 形状变化下共振轨道演化Fig.5 Evolution of resonance orbits under varying shape

2.2 共振轨道的延拓

采用连续法可以从已知共振轨道中延拓出相同共振比的轨道，降低了轨道设计难度，实现共振轨道族的快速生成。本文采用伪弧长连续法对轨道族进行求解[16-17]。伪弧长连续法在与轨道族相切的方向上取连续变量，保证所求解处在同一平面内，避免解的空间不确定性。

(28)

其中，Δk为选取的与轨道族相切方向的步长，伪弧长延拓法中的约束矩阵表达式如下：

(29)

(30)

(31)

伪弧长连续法具有鲁棒性强的优点，可保证每一次迭代所求解的唯一性。根据第2.1节确定的3∶4共振轨道，通过伪弧长连续法可求得共振轨道族如图6所示。

图6 双小行星3∶4共振轨道族Fig.6 3∶4 resonance orbit family

2.3 共振轨道的稳定性与分岔

轨道的稳定性是判断一条轨道是否适用于探测任务的重要依据之一，稳定共振轨道的轨道维持所需的燃料较省。此外，共振轨道族的稳定性改变对应于轨道特征结构的变化，可用于寻找新的共振轨道族，从而为任务轨道的设计提供更多思路。

利用单值矩阵可对轨道稳定性进行分析，定义单值矩阵M为一个轨道周期T处的状态转移矩阵值，即

(32)

F(λ)=λ4-αλ3+βλ2-αλ+1

(33)

则剩余4个特征值可令F(λ)=0求得，当|λi|≤1(i=3,…,6)时，初始状态加上其对应的特征向量方向上的微小变化，周期轨道稳定；当|λi|>1(i=3,…,6)时，周期轨道不稳定。引入稳定性指数ν，定义稳定性指数ν为互为倒数的一对特征根的平均值，即

(34)

稳定性指数ν可用来分析给定解的稳定性。对于给定的单值矩阵来说，除了均为1的一对特征根，由另外两对特征根可得到两个不同的稳定性指数νh,νv。其中，“水平”稳定性指数ν=νh反映了轨道平面中的稳定性，“垂直”稳定性指数ν=νv反映了轨道平面外的稳定性。同理，如果稳定性指数ν大于1，则轨道稳定，反之则轨道不稳定。通过引入稳定性指数，可以较为全面地考虑单值矩阵特征根对轨道稳定性带来的影响，特别地，考虑引入最大稳定性指数νmax=max{|νh|,|νv|}，则若νmax≤1，轨道稳定，反之则轨道不稳定。

对于一组轨道族的特征根来说，若其结构发生显著改变，轨道的稳定性突然发生改变，表明轨道可能出现分岔[18]，利用分岔理论可以从分岔轨道延拓得到新的共振轨道族。在共振轨道中常见的分岔类型包括正切分岔及倍增分岔。

图7为正切分岔与倍增分岔的示意图，当一对原本位于单位圆上(即|λ|=1)的共轭复数特征根在复平面上的(-1,0)处(即ν=-1处)交汇，之后沿着实数轴Re(λ)变化时，出现正切分岔。一般在正切分岔中，轨道稳定阶数将发生改变。

图7 不同类型的分岔形式Fig.7 Different types of bifurcation

当一对原本位于单位圆上的共轭复数特征根在复平面上的点(1,0)处(即ν=1处)交汇，之后沿着实数轴Re(λ)变化时，出现倍增分岔。在倍增分岔前后，轨道稳定阶数不变，新的轨道族的周期为分岔点处初始轨道周期的两倍。

利用分岔理论，可在二维平面共振轨道的基础上对双小行星系统中的三维空间中的共振轨道进行设计。选取出现分岔的共振轨道，计算其状态转移矩阵将这条分岔轨道初始状态的z分量加一个小的摄动量。固定摄动分量的值，利用改进的多次打靶法可找到一条三维共振轨道，利用连续法即可得到三维的共振轨道。

3 双小行星系统平面与三维共振轨道

根据第2节提出的共振轨道计算与延拓方法，本文选取的双体小行星系统1999 KW4[19]，求解其附近的平面和空间共振轨道族，并讨论共振轨道特性随轨道周期及能量的演化规律。双体小行星系统1999 KW4的三轴椭球体形状模型及相关参数如表1所示，系统质量比为0.0541。

3.1 平面共振轨道

图8给出椭球体-椭球体模型中的平面共振轨道，选取的共振比均小于等于1，因为大于1的共振比的轨道大部分与主星相交，即发生碰撞，且周期小、轨道速度快，对工程探测而言意义不大。

表1 1999 KW4三轴椭球体模型Table 1 Tri-axial ellipsoid model of 1999 KW4

图8 双体小行星不同共振比的平面轨道族Fig.8 Resonant orbit families with different resonant ratios

从图8所示的共振轨道图可以看出，环的个数总是和共振比p∶q中p的值一致，即共振比为1∶1，1∶2，1∶3，1∶4的共振轨道中均只有一个环，而共振比为2∶3的共振轨道有2个环，共振比为3∶4的共振轨道有3个环。

由于同一族轨道的各参数连续变化，所以选取各图中轨道延拓的内外边界轨道(非轨道族的边界)进行分析，表2列出了它们的非零初始状态、周期和雅克比常数C的值。

表2 平面共振轨道的初始状态、轨道周期和雅克比常数值Table 2 Initial state, period and Jacobi constant of the planar resonant orbits

3.2 三维共振轨道

由图9可知，三维共振轨道族与对应的平面共振轨道族的形状保持一定的相似性；共振比为1∶1的分岔轨道的初始位置与主星较近，共振比为1∶2，1∶3，1∶4的分岔轨道的初始位置则较靠近从星，共振比为2∶3的分岔轨道有两族，其中一族的初始位置靠近主星，另一族的初始位置靠近从星。

表3中列出了三维空间共振轨道的非零初始状态、周期和Jacobi常数C的值。从表3可以看出，随着三维共振轨道初始状态z0分量绝对值增大，轨道周期增大、Jacobi常数减小。从轨道周期的数据可以看出，探测器周期与双小行星系统的公转周期之比均不为整数比，随着轨道族的延拓(往z0分量的绝对值增大的方向)更接近于整数比。

4 结 论

本文主要研究了双体小行星中共振轨道设计方法。利用改进并行打靶法和分岔理论，提出双小行星系统平面共振轨道修正方法和三维共振轨道生成方法，并利用连续法实现共振轨道族的快速延拓。以双小行星系统1999 KW4为例，设计了共振比为1∶1，1∶2，1∶3，1∶4，2∶3的平面和空间共振轨道族，并分析了共振轨道的特性及轨道周期和轨道能量的变化规律。研究发现双小行星系统中共振轨道“轨道环”的个数和共振比p∶q中p的值总是保持一致。由于双小行星系统复杂的动力学环境与特性，探测器的轨道周期与较小天体的公转周期比不为整数，但随着轨道能量的增大，比值将趋近于整数。本文给出的共振轨道设计方法和轨道运动规律研究，可为未来双小行星系统探测的任务轨道设计提供一定的借鉴。

表3 共振轨道的初始状态、轨道周期和雅克比常数值Table 3 Initial state, period and Jacobi constant of three dimensional resonant orbits