基于多項式約束的三角平動點平面周期軌道設計方法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