平動點雙脈沖轉移軌道的快速計算方法

潘 迅，楊 瑞，泮斌峰，唐 碩

(1. 西北工業大學航天學院，西安710072；2. 航天飛行動力學技術重點實驗室，西安710072)

?

平動點雙脈沖轉移軌道的快速計算方法

潘 迅1,2，楊 瑞1,2，泮斌峰1,2，唐 碩1

(1. 西北工業大學航天學院，西安710072；2. 航天飛行動力學技術重點實驗室，西安710072)

針對限制性三體問題中的平動點雙脈沖轉移，提出一種高效的計算方法。通過利用基于二維插值的數值流形近似方法對流形進行近似計算，同時利用二體模型下的圓錐曲線近似流形拼接段，根據經典軌道要素推導得到完成拼接所需的速度增量，避免在優化過程中對流形的重復積分計算，以及在三體模型下對拼接段的迭代計算，從而顯著提高計算效率。然后推導得到三體問題下的主矢量理論，可將其用于對優化所得的雙脈沖轉移軌道進行燃料最優性的驗證。最后，以航天器從近地圓軌道到地月系L1點的halo軌道的雙脈沖轉移為例進行數值仿真，驗證數值流形近似算法和二體模型近似脈沖的有效性，并表明該方法在優化過程中具有高效性。

圓限制性三體問題；平動點；不變流形近似；主矢量理論；軌道優化；地月轉移

0 引 言

隨著嫦娥計劃的進行，我國的月球探測正在有條不紊的推進過程中。相比于傳統的二體模型，圓限制性三體問題模型能更精確地描述地月空間，且存在平動點、周期軌道等眾多動力學特性。月球附近存在兩個平動點L1點和L2點，前者位于地月之間，可作為地月轉移的中轉站，后者位于地月連線的延長線上，可對地球和月球背面的連續不間斷通信提供支持[1]，也可作為深空探測的基站，對這兩個平動點的合理利用有利于后續深空探測的開展。文獻[2-4]對平動點附近的空間結構，動力學特性及其在深空探測中應用進行了詳細的介紹。 2010年ARTEMIS任務的P1和P2航天器分別進入了地月L2點和地月L1點的擬周期軌道，成為世界上第一次地月平動點任務[5]。我國的嫦娥五號試驗器也于2014年對地月系L2點進行了探測。

脈沖推進作為目前的主要推進方式，短期內在實際任務的開展和應用中的地位不可動搖。對于航天器從地球到平動點的脈沖轉移軌道的優化設計，目前有較為豐富的研究成果。文獻[6]利用三體Lambert問題的迭代算法，對地球到平動點周期軌道的雙脈沖轉移進行了相關研究；文獻[7]結合網格搜索算法和遺傳算法對地月系平動點轉移軌道的雙脈沖轉移進行優化；文獻[8]利用同構映射將平面圓限制性三體問題的維度從4維降為3維，然后運用幾何方法對地球到平動點的轉移軌道進行優化；文獻[9]基于擾動流形和軌道角動量對地月系平面限制性三體問題的地球停泊軌道和地月L1點之間的多脈沖轉移軌道進行設計；文獻[10]利用限制性三體問題的動力學系統理論研究了從地球到地月系L3點halo軌道的雙脈沖轉移；文獻[11]對地月系下雙脈沖轉移進行了詳細研究，并對不同流形拼接點和不同halo軌道的轉移進行了對比分析；此外，文獻[12]針對在日地系和地月系耦合雙三體模型下的脈沖轉移進行了研究。

在三體模型下對平動點轉移軌道進行設計時，為節省能量，均利用了穩定不變流形能無動力趨向于周期軌道的特性。不變流形可通過對流形初始點進行數值積分得到。在優化過程中，需要對從近地停泊軌道進入流形的位置進行優化，即選取合適的流形拼接點。在優化過程中需要多次積分得到不同的流形狀態點作為拼接點，計算量較大，增加了優化難度。文獻[13]提出了數值流形近似計算方法，先利用二維三次卷積插值得到流形近似值，然后根據能量對其進行修正，得到更為精確的流形狀態點，避免了重復積分過程；文獻[14]也對該流形近似計算方法進行了研究，并將其應用于基于形狀法的低能轉移軌道。

在得到脈沖轉移軌道后，需要對該軌道的最優性進行驗證。Lawden[15]在1963年提出二體動力學的主矢量理論，可用于確定是否為最優飛行軌道。文獻[16]對三體問題下的月地返回軌道進行優化，并利用主矢量理論驗證了其優化結果的最優性。

本文對地月系圓限制性三體問題下的平動點雙脈沖轉移軌道進行研究。對于近地停泊軌道與不變流形的拼接段，以二體模型下的圓錐曲線進行近似，得到優化結果后再進行修正，減小對拼接段脈沖的優化難度及優化時間，同時結合數值流形近似計算方法，縮短了優化時間。與文獻中的研究相比，采用本文的方法能在極短時間內得到優化結果，顯著提高優化效率，并經主矢量理論驗證，優化結果為最優脈沖轉移軌道。

1 限制性三體問題

對于地月系圓限制性三體問題，航天器在轉移過程中同時受到兩個主天體地球和月球的影響，其在旋轉坐標系中的動力學模型為

(1)

式中：μ表示主天體的質量比，是三體系統的唯一參數，在地月系中取值為0.01215；r1、r2表示航天器與地球和月球之間的距離，分別為

(2)

平動點、周期軌道和不變流形是限制性三體問題的主要特征。不變流形分為趨近于平動點周期軌道的穩定流形和遠離周期軌道的不穩定流形，利用不變流形可進行主天體與平動點周期軌道之間的轉移軌道設計。在計算過程中，通過施加微小擾動，使航天器偏離原周期軌道，得到流形積分初始點，對其進行積分則可得到相應不變流形，其初始點選取方法為

(3)

式中：x0為周期軌道上的點，Y為對應x0的單位特征向量，上角標s和u分別表示穩定流形不穩定流形，Xs(x0)和Xu(x0)分別為穩定流形和不穩定流形的積分初始點，ε表示小擾動量，在本文中取值為2×10-6。對流形積分初始點進行積分，即可得到相應不變流形。L1點halo軌道的穩定流形如圖1所示。

圖1 L1點halo軌道的穩定流形Fig.1 The stable invariant manifolds of L1 halo orbit

2 不變流形數值近似計算

由于不變流形不能到達地球附近區域，在對轉移軌道進行設計時，航天器從近地停泊軌道出發后需要選擇合適的位置進入流形，即對流形拼接點進行優化。傳統的計算方法為通過積分得到流形點，由于在優化過程中，對選取不同拼接點時的性能進行分析對比，需要大量計算流形狀態點，而基于數值積分的計算方法需要對整條流形進行積分才能得到某個狀態點，計算量較大，計算效率較為低下。Topputo針對該問題提出了數值近似快速計算方法，通過二維插值，對流形狀態點進行插值近似，在優化過程中不需要對流形重復積分，從而在滿足計算精度的前提下，大大提高了計算效率[13]。

在積分得到流形狀態點的過程中，只與參數τ和t有關，其中τ為對周期軌道初始點進行積分，并得到式中x0的積分時間，t為對穩定流形初始點Xs(x0)的積分時間，因此通過二維插值則可得到流形狀態點。首先通過積分得到流形插值的原始數據，對積分流形的網格劃分為

(4)

式中：T1為平動點周期軌道的周期，T2為流形的最大積分時間，對于穩定流形有T2<0，對不穩定流形則有T2>0，N和M為離散點數量。根據τi和tj則可積分得到相應的穩定流形點xs,ij(τ,t)。

得到插值的原始數據xs,ij后，對任意給定的(τ,t)進行計算的插值公式為

(5)

式中：ci+n,j+m為插值樣點，h1=T1/(N-1)和h2=T2/(M-1)為離散步長，u為插值內核。

對于插值內核，設u(0)=1，u(n)=0，對其他值的計算式為

(6)

插值樣點c為六維向量，分別對應流形狀態點的三維位置矢量和三維速度矢量，計算公式為

(7)

(8)

(9)

3 二體模型脈沖近似

在得到流形拼接點后，需要對近地停泊軌道和流形進行拼接。在限制性三體問題中，三體Lambert轉移需要進行多次迭代計算，計算過程較為復雜。由于月球引力遠小于地球引力，因此在計算拼接段時，可先利用地球-航天器二體模型進行近似計算，然后在三體模型下進行修正，從而降低計算難度和減小計算量。

在二體模型下，根據經典軌道六要素，對于拼接段圓錐曲線，有

(10)

(11)

rp=a(1-e)

(12)

(13)

然后可計算得到航天器在拼接點所需施加的速度脈沖為

ΔvMI=(1-k)vt

(14)

到達近地停泊軌道后需要施加的速度脈沖為

(15)

可得到二體模型近似下所需的速度增量為

Δvapp=ΔvE+ΔvMI

(16)

Φ(tp+Δt,xpr+Δxr)=Φ(tp,xpr)+

(17)

(18)

兩個修正量對應兩個約束方程，因此可求解得到Δt和p，從而計算得到三體模型下的拼接段以及完成拼接所需的速度增量ΔvEm和ΔvMIm。

4 主矢量理論

在1963年提出的主矢量理論中，Lawden給出的最優脈沖的必要條件為[15]：

4)對于所有內點脈沖(初始時刻脈沖和終端時候脈沖除外)成立。

本文對限制性三體問題下的主矢量進行推導。在三體模型中，推力加速度近似為脈沖表示

(19)

(20)

式中：ti為施加脈沖時刻。

根據最優控制理論，構造哈密爾頓方程為

(21)

(22)

(24)

得到脈沖優化結果后，根據所施加的脈沖計算得到λv(t0)和λv(tf)，再以λr(t0)為迭代變量，根據邊界條件λv(tf)求解出λr(t0)，然后積分得到主矢量λv在時間區間[t0,tf]上的值，判斷是否滿足所有Lawden必要條件，從而可判斷優化得到雙脈沖轉移軌道是否為最優脈沖轉移軌道。

5 數值仿真

考慮航天器從近地軌道出發，通過施加第一次脈沖進入拼接段，在流形拼接點處施加第二次脈沖進行L1點halo軌道的穩定流形，然后沿流形無動力的趨近于halo軌道。地球平均半徑為6378.137km，近地停泊軌道為高度300km的圓軌道，L1點halo軌道的幅值Az為5000km。本算例中的計算在CPU為i7-4700MQ，主頻為2.40GHz，內存為8G的計算機上進行，計算平臺為MATLAB2015a。

考慮到實際任務中的飛行時間限制，流形積分時間上限為50天，同時為保證數值近似流形的精度，在halo軌道上均勻選取200個點，在每條流形上選取2001個點，作為流形插值數據庫，即式中N=200，M=2001。為驗證數值近似流形的精度，將其與積分得到的流形點進行誤差對比分析。定義流形誤差為

(25)

誤差結果如表1所示，誤差的對數lg(ξ)的分布情況如圖2所示，可知平均精度為1.6529×10-6，相當于是0.6353km，在地月轉移軌道設計過程中，該精度已屬于較高精度，因此認為數值近似流形的精度足夠滿足轉移軌道初步設計的要求。

表1 數值近似流形的誤差情況Table 1 The error of the approximation invariant manifolds

選取τ=0時的流形，在該流形積分時間為[-50, -10]天的范圍內，根據時間均勻選擇160個點，對二體模型近似脈沖和三體模型修正后脈沖進行對比，如圖3所示，其中ΔvEm和ΔvMIm為三體模型修正后所需的脈沖，對比分析可知，二體模型下的近似脈沖能很好地估計完成拼接段轉移所需的速度增量。

圖2 數值近似流形的誤差分析lg(ξ)Fig.2 The error pattern lg(ξ) of approximation invariant manifolds over the evaluation grid

對于拼接點的確定，利用MATLAB自帶函數fmincon進行優化，優化變量為[τ,t]。在優化過程中，只利用二體模型脈沖近似，在得到優化結果后再在三體模型下進行修正，修正精度要求為1km。由于fmincon是局部優化函數，選取不同的初始值會得到不同的優化結果。多次選取不同的初始值，可得到4組優化結果，與圖3中的單條流形的4個極小值點相對應。優化結果數據如表2所示，其中優化時間為只考慮fmincon函數優化的時間，4組結果所需優化時間均小于0.1s。這4組結果對應的轉移軌道分別如圖4～7所示，其中左圖為旋轉坐標系下的轉移軌道，右圖為地心慣性系下的轉移軌道。從圖4～7的(a)圖可以看出，優化得到的4個拼接點均為流形的遠地點，這是由于在遠地點受到的地球引力較小，克服引力所消耗的能量較小，脈沖的利用率較大，因此完成轉移所需的總速度增加較小；從右圖可以看出，拼接點在流形上的位置相差約一個周期，結合表2中航天器沿流形的運動時間可知，流形繞地球運動一周所需時間約為10天。

圖3 二體模型近似脈沖與三體模型修正后脈沖對比Fig.3 Comparison of the two-body approximate impulse and the impulse modified in CRTBP

優化結果ΔvMImΔvLEOmΔvtot轉移段時間/d流形段時間/d優化時間/sA0．61633．06793．68423．858417．60340．0599B0．56103．06333．62443．447028．00000．0614C0．56793．06333．63123．416638．00450．0495D0．59073．06363．65443．433248．03220．0835Larsen等[7]0．663．063．723．018．8———Pontani等[8]0．5613．0163．57720．183．21062．8?

圖4 雙脈沖轉移優化結果AFig.4 The results of two-impulse transfers: A

Larsen等[7]和Pontani等[8]均對航天器從地球到L1點周期軌道的雙脈沖轉移進行了研究，但僅考慮xy平面內的運動，目標軌道為L1點Lyapunov軌道。前者的近地停泊軌道高度為300km，后者的近地停泊軌道高度為500km，優化結果在表2中描述。文獻[7]所需的速度增量小于本文優化結果主要是因為其近地停泊軌道的高度較高，克服地球引力所需能量較小；但對于相同高度的近地停泊軌道，本文的B組結果不論從轉移時間還是所需速度增量，都要優于文獻[8]中的結果。文獻[11]對地月系下的雙脈沖轉移進行了詳細研究，當目標halo軌道的橫坐標為319052km時，對拼接點進行優化時，得到多個極小值點，該拼接點均為流形遠地點，所需最小速度增量約為3.62km/s，與本文的優化結果類似。更為重要的是，利用文獻[7]中的網格法進行優化時需耗費較長時間，文本作者復現文獻[8]中的結果需用時1062.8s，利用[11]中計算單條轉移軌道時也需要對施加脈沖進行打靶求解，而利用數值流形近似和二體模型脈沖近似方法，積分得到流形插值數據需耗時26.7s，然后利用fmincon函數求得上文中一個極小值解所需優化時間不到0.1s，因此本文的快速計算方法顯著提高了計算效率。

圖5 雙脈沖轉移優化結果BFig.5 The results of two-impulse transfers: B

圖6 雙脈沖轉移優化結果CFig.6 The results of two-impulse transfers: C

圖7 雙脈沖轉移優化結果DFig.7 The results of two-impulse transfers: D

圖8 B組優化結果的主矢量曲線Fig.8 The primer vector curve of result B

6 結 論

本文對限制性三體問題的平動點雙脈沖轉移軌道的快速計算方法進行研究。利用數值流形近似方法對流形進行插值并修正，從而避免了在優化過程中重復積分計算流形。對于確定的流形拼接點，以二體模型下的圓錐曲線作為拼接段的近似，根據經典軌道要素推導得到完成拼接所需的速度增量，再在三體模型下進行修正。最后利用推導得到的三體模型下的主矢量理論對優化結果是否滿足燃料最優進行驗證。通過數值仿真表明，結合這兩種方法能在較短時間內得到優化結果，計算效率得到顯著提升；并選取其中一條轉移軌道，利用主矢量理論驗證該脈沖轉移軌道的燃料最優性，顯示利用該方法在一定程度上能得到燃料最優轉移軌道，從而表明該方法的合理性和高效性。

[1] Farquhar R, Kamel A. Quasi-periodic orbits about the translunar libration point [J]. Celestial Mechanics and Dynamical Astronomy, 1973, 7(4): 458-473.

[2] 徐明. 平動點軌道的動力學與控制研究綜述[J]. 宇航學報, 2009, 30(4):1299-1313. [Xu Ming. Overview of orbital dynamics and control for libration point Orbits [J]. Journal of Astronautics, 2009, 30(4):1299 -1313.]

[3] 徐明, 徐世杰. 地-月系平動點及halo軌道的應用研究[J]. 宇航學報, 2006, 27(4):695-699. [Xu Ming, Xu Shi-jie. The application of libration points and halo orbits in the Earth-Moon system to space mission design [J]. Journal of Astronautics, 2006, 27(4):695 -699.]

[4] 侯錫云, 劉林. 共線平動點的動力學特征及其在深空探測中的應用[J]. 宇航學報, 2008, 29(3):736-747. [Hou Xi-yun, Liu Lin. The dynamics and applications of the collinear libration points in deep space exploration [J]. Journal of Astronautics, 2008, 29(3):736 -747.]

[5] Cosgrove D, Frey S, Bester M, et al. Navigating THEMIS to the ARTEMIS low-energy lunar transfer trajectory[C]. SpaceOps 2010 Conference, Alabama, USA. April 25-30, 2010.

[6] 連一君. 基于三體平動點的低能轉移軌道設計研究[D]. 長沙：國防科學技術大學, 2008. [Lian Yi-jun. Research on low-cost transfer trajectories design based on three-body libration points [D]. Changsha: National University of Defense Technology, 2008. ]

[7] Larsen A, Anthony W, Critz T, et al. Optimal transfers with guidance to the Earth-Moon L1and L3libration points using invariant manifolds: a preliminary study [C]. AIAA/AAS Astrodynamics Specialist Conference, Minnesota, USA. August 13-16, 2012.

[8] Pontani M, Teofilatto P. Low-energy Earth-Moon transfers involving manifolds through isomorphic mapping [J]. Acta Astronautica, 2013, 91(10):96-106.

[9] 張漢清, 李言俊. 地月三體問題下L1-地球低能轉移軌道設計[J]. 哈爾濱工業大學學報, 2011, 43(5):84-88. [Zhang Han-qing, Li Yan-jun. Design of L1-Earth low energy transfer trajectory in Earth-Moon three-body problem [J]. Journal of Harbin Institute of Technology, 2011, 43(5):84-88.]

[10] Davis K, Born G, Butcher E. Transfers to Earth-Moon L3 halo orbits [J]. Acta Astronautica, 2013, 88(3):116 -128.

[11] Parker J S, Born G H. Direct lunar halo orbit transfers [J]. The Journal of the Astronautical Sciences, 2008, 56(4):441-476.

[12] Parker J S, Anderson R L, Peterson A. Surveying ballistic transfers to low lunar orbit [J]. Journal of Guidance Control & Dynamics, 2013, 36(5):1501 -1511.

[13] Topputo F. Fast numerical approximation of invariant manifolds in the circular restricted three-body problem [J]. Communications in Nonlinear Science & Numerical Simulation, 2016, 32:89-98.

[14] 張仁勇. 深空探測小推力低能轉移軌道設計方法研究[D]. 西安：西北工業大學, 2015. [Zhang Ren-yong. Research on low-thrust low-energy transfer trajectory design for deep space exploration [D]. Xi’an: Northwestern Polytechnical University, 2015. ]

[15] Lawden D F. Optimal trajectories for space navigation [J]. Mathematical Gazette, 1964, 48(366):478.

[16] Shen H, Casalino L. Indirect optimization of three-dimensional multiple-impulse Moon-to-Earth transfers [J]. The Journal of the Astronautical Sciences, 2014, 61(3):1-20.

通信地址：陜西省西安市碑林區友誼西路127號西北工業大學251信箱(710072)

電話：15191433469

E-mail：xpan2012@gmail.com

An Efficient Calculation Method for Two-Impulse Transfers to Libration Point

PAN Xun1,2, YANG Rui1,2, PAN Bin-feng1,2,TANG Shuo1

(1. School of Astronautics, Northwestern Polytechnical University, Xi’an 710072, China;2. National Key Laboratory of Aerospace Flight Dynamics, Xi’an 710072, China)

An efficient calculation method is proposed for the optimization of the two-impulse transfers to the librations in circular restricted three-body problem. With the invariant manifolds approximation based on a two-dimensional interpolation, the computation of the great number of manifold insertions is much more efficient. While the conic curve in a two-body problem used as an approximation of the transfer segment connecting the low Earth orbit and the stable manifold, the maneuvers can be deduced with the classical orbital elements, thus the iterative computation of the maneuvers is avoided in the optimization process. Then the primer vector of CRTBP is deduced to verify the optimality of the two-impulse transfer. At last, a numerical simulation of the two-impulse transfer from a low Earth orbit to aL1point halo orbit is presented, verifying the validity of the two-body impulse approximation and the manifold approximation, and demonstrating the efficiency of the method.

Circular restricted three-body problem; Libration point; Invariant manifolds approximation; Primer vector theory; Trajectory optimization; Earth-Moon transfers

2017-01-09;

2017-04-21

國家自然科學基金(11672234)

V448.2

A

1000-1328(2017)06-0574-09

10.3873/j.issn.1000-1328.2017.06.003

潘 迅(1990-)，男，博士生，主要研究方向為飛行動力學與控制，深空探測軌道優化設計。