基于多项式约束的三角平动点平面周期轨道设计方法1)

钱霙婧 翟冠峤 张 伟

(北京工业大学机电学院，北京100124)

基于多项式约束的三角平动点平面周期轨道设计方法1)

钱霙婧 翟冠峤 张 伟2)

(北京工业大学机电学院，北京100124)

平动点是圆型限制性三体问题中的五个平衡解.其中，三角平动点在平面问题中具有“中心×中心”的动力学特性，其附近存在着大量的周期轨道，研究这些周期轨道的构建方法在深空探测中具有理论及工程意义.本文从振动角度分析周期轨道，通过多项式展开法构建出主坐标下周期轨道两个运动方向之间的渐近关系.从新的角度分析了系统的动力学特性和平面周期运动两个方向内在关联以及物理规律.这种多项式形式的关系式，可以作为约束条件用于数值微分修正算法中，通过迭代的方式寻找周期轨道.数值仿真算例验证了方法的正确性及精确性.文章从振动的角度对周期轨道进行分析，改进了微分修正算法.提出的方法可以被拓展至圆型/椭圆型限制性三体问题的三维周期轨道构建中.

平面圆型限制性三体问题，平动点，多项式展开法

引言

平动点是圆型限制性三体问题中的 5个平衡解，包括3个共线平动点L1,L2和L3以及2个三角平动点L4和L5.其中，三角平动点具有“中心×中心”的动力学特性，其附近存在着大量的周期轨道，可以被用于构建空间中转站，编队导航等[1-4].研究这些轨道对深空探测具有理论价值及工程意义.

目前，已有学者对于平动点附近的周期轨道构建问题进行了解析以及数值研究.在解析求解方面，Richardson[5]应用Lindstedt-Poincar´e(L-P)法给出了圆型限制性三体问题共线平动点附近Halo周期轨道的三阶解析解.Erdi[6]和Zagouras[7]基于小参数展开法分别推导了三角平动点附近周期轨道的三阶和四阶解析解.这为动平衡点附近轨道的分析和研究奠定了基础.近期，Lei等[8]和Zhao[9]通过L-P法构建了三角平动点附近周期轨道的任意阶解析解.在发展解析解的同时，应用数值方法设计平动点附近周期轨道的思想也在飞速发展.Goodrich[10]，Bray等[11]，Zagouras等[12]在解析解的基础上，使用微分修正的方法分别求得了太阳系中不同系统下平动点附近的轨道数值解.Howell[13]则详细推导论述了圆型限制性三体问题中拟周期轨道的数值求解方法，即多步打靶法(multiple shooting).这种方法适用性强，可以根据不同的终端约束需求做出改动[14-15].随后，Andreu等[16-17]应用多步打靶法，通过初值积分的方式来得到多步打靶的拼接点，最终完成轨道设计.许多其他学者也研究了通过利用数值方法来设计轨道[18-20].

综上所述，在解析分析的研究中，传统摄动方法着重于修正线性条件下的振幅与频率，使其更加接近非线性条件下的真实运动，但是却很少关注运动中各维度之间的联系.数值分析的研究则主要以解析解为轨道的理想猜测初值，通过微分修正以及轨道本身的对称性来得到周期轨道，或者是以解析解轨道上的点为拼接点(patch points)为约束条件采用多步打靶法来求解连续轨道.然而，拼接点的求解本身就是迫切需要解决的难点问题[21-22].

在振动理论最新的进展中，由Shaw等[23-25]基于模态分析的思想提出了的一种基于多项式展开理论的求解方法，为周期轨道求解问题提供了新的思想.这种多项式展开方法定义了一种不变的相空间关系[26]，从而得到两自由度之间的多项式关系.能为数值求解真实力学模型下的周期轨道提供满足物理规律的约束条件[27].受启发于这种思想，本文首先采用多项式展开求解运动方程，得到平面圆型限制性三体问题三角平动点附近周期轨道两个维度之间的运动关系.利用微分修正的思想，采用两个维度之间的运动关系作为约束条件，进行多次迭代得到满足误差精度的周期轨道.

文中所提出的采用多项式展开方法得到的运动关系可以清晰地反映平面圆型限制性三体问题模型中三角平动点附近周期轨道两个自由度之间的关系，为分析其轨道动力学特性提供理论依据，并且可以为数值迭代求解周期轨道提供约束关系，为设计真实力学模型下的飞行器轨道提供借鉴.

1 基本动力学模型

1.1 平面圆型限制性三体问题

在平面圆型限制性三体问题中，一个质量相对无限小的第三体在两个围绕其公共质心做圆周运动的主天体的引力作用下做运动.

假设质量较大的主天体P1质量为m1，质量较小的主天体P2质量为m2.两个主天体绕其共同的质心C做匀速圆周运动.选取质心会合坐标系进行问题的研究，记为C-XY.其原点为C，XY平面为两个主天体相对运动平面，X轴由主天体P1指向主天体P2.如图1所示.

图1 坐标系示意图Fig.1 Schematic for coordinate systems

为了计算方便，通常取两个主天体质量之和、两个主天体间的距离分别为质量与长度的单位，即定义µ=m2/(m1+m2)为质量参数.则主天体P1质量为1-µ，坐标为(-µ,0).主天体P2质量为µ，坐标为(1-µ,0).小天体在此会合坐标系中的运动方程为

Ω为系统中的拟势能函数，通常表示为[28]

其中，R1与R2分别代表小天体到主天体P1与P2的距离

1.2 三角平动点附近的运动方程展开

为了描述三角平动点附近的运动，将坐标系移动到三角平动点.本文以L4点为例，将坐标系原点移动到L4，新的坐标轴x与y和原坐标系X与Y轴平行，如图1所示.在此坐标系下将原运动方程(1)按Legendre展开，可以表示为[29]

式中，Pn为n阶的Legendre多项式，

由于式(4)的左边存在线性耦合项，不利于后续计算，同时为了更清晰地体现三角平动点附近运动的几何特征，选择将原L4-xy坐标系旋转θ角，得到一个新坐标系L4-ξη[28]，如图1所示.L4-ξη为这一系统的主坐标(principle coordinate system).

引入新变量(ξ,η)代替(x,y)，其关系可以表示为在新坐标下，式(4)可以表示为

由此，消除了方程(4)左端的线性耦合项，线性椭圆运动轨道的长轴和短轴就分别位于这一新坐标系的ξ轴和η轴上.

2 周期运动的多项式展开分析

本节使用多项式展开法研究两个方向运动之间的关系，得到它们对系统动力学特性的影响，为分析三角平动点附近周期运动的运动形式以及动力学特性提供参照.这种方法的核心思想在于选取一个方向为基方向，基方向上的运动状态(位置和速度)为一组空间的二维基状态，将另一方向的运动描述为与基方向状态相关的多项式形式[23-25].通过求解多项式系数的方式寻找两个方向运动之间的关系.

将式(7)Legendre展开后表示成如下形式

这里的ξ与η分别代表两个运动方向的位移.为了得到比较精确的结果，将式(8)保留前3次项，可将原运动方程改写为

其中，A2i,B2i,A3i和B3i的具体表达式见附录.

选取ξ方向的位置和速度为周期运动时的基状态，即令

根据文献[23-25]，η方向的运动可以被描述为与ξ方向相关的如下多项式形式

其中,ai与bi是待定系数.通过对这些系数的求解，可以得到两个方向上位移与速度的关系.

将式(15)与式(16)分别代入式(9)的两个方程中得到

将式(15)对时间求导，可得

将式(16)对时间求导，可得

对比式(16)与式(20)、式(18)与式(22)的一次项系数可得到

对比二次项系数可得到

对比三次项系数可得到

2.1 线性项关系

求解式(23)，可以得到两组线性化运动方程的关系，得到

上述方程表明，对于平面问题，飞行器在三角平动点附近的线性运动是椭圆形的.其一个方向上的位移与另一个方向上的速度成比例，即存在π/2的相位差，当其中一个方向上的速度达到最大值时，另一个方向则刚好经过其零点.即如果取初始条件满足式(28)时，会使得原运动方程的其中一个运动分量为零，而只有另外一个运动分量，此时xy平面内的运动对应为周期运动.而式(26)与式(27)所表示的两个模态就对应两类周期轨道：长周期轨道和短周期轨道.

将式(28)代入到式(9)中，可以得到运动方程的线性化频率

这里的负号对应式(26)，正号对应式(27).根据平动点附件周期运动线性稳定性分析[30]，在质量参数µ小于µc=0.038520896(Routh极限值)时，上述两个频率分别代表了三角平动点附近运动的短周期和长周期运动.

2.2 二次非线性项关系

继续求解式(24)，可以得到运动方程二次非线性系数的解析表达式.令

2.3 三次非线性项关系

解得运动方程三次非线性系数方程(25)的解析形式，可得

其中，S1～S7的表达式见附录.

由此，用多项式展开的方法得到了平面限制性三体问题中三角平动点附近运动的两个自由度之间关系表达式中所有系数.三角点附近周期运动的两个方向ξ和η满足如下方程

式(33)和式(34)是描述平面周期运动两个方向内在关联物理规律的方程，在后文中将用于数值求解中微分修正的约束条件.

3 微分修正方法

本节主要介绍轨道微分修正方法，通过前文求得的运动方程两自由度之间的关系修正初值，为后文高精度轨道的求解提供理论基础.

微分修正法在本质上是一种迭代的打靶法[31]，是对广义约束和自由变量理论的应用.对于这样的问题，存在着控制变量和约束变量，而这两种变量则是通过状态转移矩阵联系起来的.微分修正方法通过状态转移矩阵描述约束变量相对控制变量微小改变的敏感性，并通过调整控制变量使约束变量达到期望值.

求解式(7)的线性部分，容易得到其在L4-ξη坐标系下线性解的形式

其中，ω如式(29)所示，α0和β分别代表幅值与相角.ξ和η代表两个方向上的位移，vξ和vη代表两个方向上的速度.

确定幅值和相角后，由式(35)得到L4-ξη坐标系下微分修正的迭代初值，记为固定时间进行数值积分得到终端状态值运用前文求得的两个自由度之间的关系式(33)和式(34)以确定终端目标的状态值为即

为了保证轨道设计过程中轨迹满足设计目的，Xf应该等于的值.即需要满足如下4个约束条件

对于初始猜测的特解X，可以将F(X)在X0附近进行泰勒展开，并保留其一阶项之后表示为

其中，∂F(Xj)/∂X为状态转移矩阵，δX0则为初始猜测X0需要进行修正的量.对于绝大多数非线性问题而言，需要通过多次迭代的方式得到最终的修正量，因而，联立式(37)和式(38)，并用j表示迭代的次数，可以得出

直到Xj+1=X∗或者||F(Xj+1)||≤ε而终止，其中ε为需要的收敛精度.

显然约束变量等于控制变量，迭代方程有唯一解，可以得到修正量δXj的表达式

这样，通过多次迭代修正初值量X使其满足精度要求.

4 数值仿真与对比

本节采用地--月--飞行器圆型限制性三体问题中三角平动点附近的周期轨道为例进行仿真计算.其质量参数为µ=0.012150568.取线性解式(35)中，幅值α0=0.05，相角β=π.以此为初值，在平面限制性三体问题全模型，即式(1)条件下，进行微分修正.

4.1 短周期轨道

首先对短周期轨道的情况进行算例仿真，即取ω=0.9545009306377使用线性化模型得到的状态初值进行数值积分，固定积分时间为5个单位化时间，得到积分终值.迭代8次后积分初值满足精度，迭代完毕.图2所示为C-XY坐标系下，8次修正过程积分初值点的变化情况(仅选取第1次、第3次、第6次、第8次修正过程示意).曲线代表线性解析解轨道，而点划线表示数值修正初值所得轨道，星号代表积分初值点.

分别对修正前以及修正后的初值进行50个单位时间积分.图3(a)所示为线性解积分后得到的轨道图，图3(b)所示为微分修正后初值所积分得到的轨道图.通过对比可以发现，由线性解积分得到的轨道不能保持周期特性，但由于三角平动点的“中心”特性，其轨迹不会远离三角平动点，而是呈现拟周期形式的运动.而由修正后的初值进行积分所得到的轨道则可以保持周期运动.数值仿真验证了文中方法的有效性.

图2 短周期轨道微分修正结果图Fig.2 Dif f erential correction for the short-period periodic motion

图3 微分修正前后短周期轨道图Fig.3 Comparison of the trajectories numerically integrates with the linearized initial condition and corrected initial condition for the short-period motion

4.2 长周期轨道

与短周期轨道的验证方式相同，通过式(29)得到长周期轨道频率ω=0.2982079365337，使用线性解(35)得到的状态初值.对于长周期轨道，仿真中取固定积分时间为15个单位化时间，得到积分终值.迭代8次后积分初值满足精度.图4所示为C-XY坐标系下，8次修正过程积分初值点的变化情况(仅选取第1次、第3次、第6次、第8次修正过程示意).曲线代表线性解析解轨道，而点划线表示数值修正初值所得轨道，星号代表积分初值点.

图4 长周期轨道微分修正结果图Fig.4 Dif f erential correction for the long-period periodic motion

分别对修正前以及修正后的初值进行50个单位时间积分.图5(a)所示为线性解积分后得到的轨道图，图5(b)所示为微分修正后初值所积分得到的轨道图.与短周期轨道类似，通过对比可以发现，由线性解积分得到的轨道不能保持周期特性，但由于三角平动点的“中心”特性，其轨迹不会远离三角平动点，而是呈现拟周期形式的运动.而由修正后的初值进行积分所得到的轨道则可以保持周期运动.数值仿真验证了文中方法的有效性.

图5 微分修正前后长周期轨道图Fig.5 Comparison of the trajectories numerically integrates with the linearized initial condition and corrected initial condition for the long-period motion

图5 微分修正前后长周期轨道图(续)Fig.5 Comparison of the trajectories numerically integrates with the linearized initial condition and corrected initial condition for the long-period motion(continued)

值得注意的是，文中方法的主要目标是寻找平面运动两个方向上的非线性关系，是一类几何约束.这一类约束可以维持周期轨道的几何特征，但却无法约束轨道的振幅.长/短周期轨道的仿真中，修正后的轨道振幅都相比较修正前放大.这一特性可以通过添加振幅约束得以改善.

5 结论

本文提出了基于多项式约束构建平面圆形三角平动点附件周期轨道的方法.通过多项式展开方法得到的运动关系可以清晰地反映平面圆型限制性三体问题模型中三角平动点附近周期轨道两个自由度之间的运动关系.为分析其轨道动力学特性提供理论依据，并揭示了周期运动的物理规律.文章创新性的将运动规律作为微分修正的约束条件，通过多次迭代构建平动点附近周期轨道.

本文所使用的研究方法同样适用于共线平动点附近的轨道构建问题，并可以被拓展至圆型/椭圆型限制性三体问题的三维周期轨道构建中.

1 Brouwer D,Clemennce GM.Methods of Celestial Mechanics.2nd edn.Academic Press,1985

2 Beuler G.Methods of Celestial Mechanics.Berlin,Heideberg：Springer-Verlag,2005

3 刘林.人造地球卫星轨道力学.北京：高等教育出版社,1992(Liu Lin.Orbit Dynamics of the Artificia Earth Satellite.Beijing：The Higher Education Press,1992(in Chinese))

4 刘林.航天器轨道理论.北京：国防工业出版社,2000(Liu Lin. The spacecraft orbit theory.Beijing：National Defence Industry Press,2000(in Chinese))

5 Richardson DL.Analytic construction of periodic-orbits about the collinear points.Celestial Mechanics,1980,22(3)：241-253

6 Erdi B.3-dimensional motion of trojan asteroids.Celestial Mechanics,1978,18(2)：141-161

7 Zagouras CG.3-dimensional periodic-orbits about the triangular equilibrium points of the restricted problem of 3 bodies.Celestial Mechanics,1985,37(1)：27-46

8 Lei HL,Xu B.High-order solutions around triangular libration points in the elliptic restricted three-body problem and applications to low energy transfers.Communications in Nonlinear Science and Numerical Simulation,2014,19(9)：3374-3398

9 Zhao L.Quasi-periodic solutions of the spatial lunar three-body problem.Celestial Mechanics&Dynamical Astronomy,2014, 119(1)：91-118

10 GoodrichE.Numericaldeterminationofshort-periodtrojanorbitsin the restricted three-body problem.The Astronomical Journal,1966, 71(2)：88-93

11 Bray T,Goudas L.Doubly symmetric orbits about the collinear Lagrangian points.The Astronomical Journal,1967,72(2)：202-213

12 Zagouras C,Kazantzis P.Three-dimensional periodic oscillations generating from plane periodic ones around the collinear Lagrangian points.Astrophysics and Space Science,1979,61(2)：389-409

13 Howell KC.Three-dimensional periodic“halo”orbits.Celestial Mechanics,1984,32(1)：53-71

14 Hughes S,Cooley D,Guzman JA.direct method for fuel optimal maneuversofdistributedspacecraftinmultiplefligh regimes//Space Flight Mechanics Meeting,Copper Mountain,Colorado,2005

15 Marchand B,Howell KC.A spherical formations near the libration points in the sun-earth/moon ephemeris system//14th AAS/AIAA Space Flight Mechanics Conference,Maui,Hawaii,2004

16 Andreu MA,Simo C.Translunar halo orbits in the quasi-bicircular problem//NATO ASI,1997：309-314

17 Andreu MA.Dynamic in the center manifold around L2in quasibicircular problem.Celestial Mechanics and Dynamical Astronomy, 2002,84(2)：105-133

18 Parker J,Born GH.Direct lunar halo orbit transfers//In 17th AAS/AIAA Space Flight Mechanics Meeting,2007,January 28-February 1(AAS)：07-229

19 Perozzi E,Salvo A.Novel spaceways for reaching the moon：an assessment for exploration.Celestial Mechanics&Dynamical Astronomy,2008,102(1-3)：207-218

20 Parker J,Martin L.Shoot the Moon 3D//AAS/AIAA Space Flight Mechanics Meeting,2005,August 7-11(Paper AAS-05-383)

21 Qian YJ,Liu Y,Zhang W,et al.Station keeping strategy for quasiperiodic orbit around Earth-Moon L2 point//Proceedings of the Institution of Mechanical Engineers,Part G：Journal of Aerospace Engineering,2016

22 侯锡云.平动点的动力学特征及其应用.[博士论文].南京：南京大学,2008(Hou Xiyun.Dynamic characteristics and applications of the libration point.[PhD Thesis].Nanjing：Nanjing University, 2008(in Chinese))

23 Shaw SW,Pierre C.Normal-modes for nonlinear vibratory-systems.Journal of Sound and Vibration,1993,164(1)：85-124

24 ShawSW,PierreC.Normal-modesofvibrationfornonlinearcontinuous systems.Journal of Sound and Vibration,1994,169(3)：319-347

25 Shaw SW.An invariant manifold approach to nonlinear normalmodesofoscillation.JournalofNonlinearScience,1994,4(5)：419-448

26 Arquier R.Two methods for the computation of nonlinear modes of vibrating systems at large amplitudes.Computers&Structures, 2006,84(24-25)：1565-1576

27 Vakakis AF,Manevitch L,Mikhlin Y,et al.Normal Modes and Localization in Nonlinear Systems.New York：John Wiley,1996

28 Szebehely V.Theory of Orbits.New York and London：Academic Press,1967

29 Jorba A,Masdemont J.Dynamics in the centre manifold of the collinear points of the restricted three body problem.Physica D, 1999,132：189-213

30 刘林,侯锡云.深空探测器轨道力学.北京：电子工业出版社,2012 (Liu Lin,Hou Xiyun.Orbital Mechanics of the Deep Space Probe. Beijing：Publishing House of Electronics Industry,2012(in Chinese))

31 Keller HB.Numerical solution of two point boundary value problems.Society for Industrial and Applied Mathematics,1976

附录

PLANAR PERIODIC ORBIT CONSTRUCTION AROUND THE TRIANGULAR LIBRATION POINTS BASED ON POLYNOMIAL CONSTRAINTS1)

Qian Yingjing Zhai Guanqiao Zhang Wei2)
(College of Mechanical Engineering，Beijing University of Technology，Beijing100124，China)

Libration points are the fve equilibrium solutions in the circular restricted three-body problem(CRTBP). The linearized motions around triangular libration points are typical center×center type.Studies about probes moving around orbits in the vicinity of the libration points have theoretical significance From the vibrational point of view,the polynomial series are used to derive approximately the relations in dif f erent directions during periodic motions,which provides a new point of view to exploring the dynamics and analyzing the overall characteristics of the whole system with general rules.The nonlinear relations in polynomial form between the directions of the planar motions can be treated as constraints to obtain the solutions by numerical integration.Numerical simulations verify the efficiency of the proposed method.The methodology of deriving topological relations has the potential to be extended to circular/elliptical R3BP in three dimensional cases.

planar circular restricted three-body problem，libration point，polynomial expansion method

V412.4

A doi：10.6052/0459-1879-16-215

2016-08-01收稿，2016-10-11录用,2016-10-18网络版发表.

1)国家自然科学基金资助项目(11402007).

2)张伟，教授，主要研究方向：非线性动力学.E-mail：sandyzhang0@163.com

钱霙婧,翟冠峤,张伟.基于多项式约束的三角平动点平面周期轨道设计方法.力学学报,2017,49(1)：154-164

Qian Yingjing,Zhai Guanqiao,Zhang Wei.Planar periodic orbit construction around the triangular libration points based on polynomial constraints.Chinese Journal of Theoretical and Applied Mechanics,2017,49(1)：154-164