星歷模型地月系統平動點擬周期軌道設計研究

劉 剛，黃 靜，郭思巖

(1.上海市空間智能控制技術重點實驗室,上海 201109； 2.上海航天控制技術研究所,上海 201109)

星歷模型地月系統平動點擬周期軌道設計研究

劉 剛1,2，黃 靜1,2，郭思巖1,2

(1.上海市空間智能控制技術重點實驗室,上海 201109； 2.上海航天控制技術研究所,上海 201109)

為改進使用圓形限制性三體模型設計軌道時缺乏攝動分析的不足，提高軌道設計的精度，對星歷模型地月系統平動點擬周期軌道設計方法進行了研究。在不設置假設條件的前提下，考慮月球真實軌道及地球、月球和太陽的影響，在地心J2000慣性坐標系中建立平動點附近航天器高精度的星歷模型。以圓形限制性三體中的周期軌道作為迭代初值，用星歷表數據對軌道進行拼接獲得所需的擬周期軌道；用多步打靶法替代單步微分修正進行迭代，對軌道上各節點進行校正以獲得所求的擬周期軌道，給出了軌道設計步驟。仿真結果表明：所提方法可有效獲得地月系統平動點附近擬周期軌道，提供滿足真實動力學環境的軌道，有效節約軌道保持所需的燃料。該方法有較大的工程應用價值。

地月系統； 平動點； 限制性三體問題； 太陽光壓； 引力攝動； 擬周期軌道； 多步打靶法； 星歷表

0 引言

平動點是圓形限制性三體問題(CRTBP)中的5個平衡位置，在平動點附近存在大量穩定、不穩定、周期和擬周期運動。平動點附近無殘余大氣、空間碎片、原子氧等干擾因素的影響，地球紅外射線、重力梯度和地球磁場對航天器的影響也非常微弱，非常適于開展各種科學探測任務。因此，對平動點附近軌道設計問題的研究有廣泛的應用前景和價值。目前，對平動點附近動力學和軌道設計進行了大量研究。文獻[1]用Linstedt-Poincaré(LP)攝動法推導了地月系統L2點附近的周期及擬周期軌道三階及四階近似解析解。文獻[2]提出了用LP法構建一般Halo軌道的三階近似解析解的方法，該結果目前仍被用于求解Halo軌道數值計算的初值。隨著現代計算機技術的發展，數值計算方法被廣泛用于周期軌道的求解，極大地促進了平動點軌道相關研究的進展。文獻[3]在FARQUHAR給出的地月系統共線平動點近似解析解的基礎上，用數值方法得到了地月系統共線平動點附近的多族Halo軌道。文獻[4]同樣在FARQUHAR研究的基礎上，論述了CRTBP三維Halo軌道及相應軌道簇的數值求解方法，并在原單步微分校正的基礎上，針對Lissajous軌道不具備周期性的特點，采用多步打靶法實現了擬周期軌道的數值求解。上述研究成果大多是基于CRBTP模型得到的，但在實際的地月空間中，CRTBP模型不能完全反映真實的動力學環境，用該模型設計的軌道與實際存在較大差別，會明顯增加軌控所需的燃料。因此，對CRTBP模型的改進又進行了大量研究。文獻[1]在地月系統共線平動點軌道設計中考慮了太陽引力的作用，用月球軌道平均要素及引力展開成無窮級數后的有限項逼近。文獻[5]在CRTBP基礎上提出限制性四體模型，并以CRTBP的Halo軌道為初值，用數值法得到了四體問題的周期軌道。此后，Howell課題組在研究CRTBP問題時采用噴氣推進實驗室(JPL)的星歷表(Ephemeris)數據，建立了高精度的動力學模型，并將其用于多個任務設計，但未給出更詳盡具體的軌道設計細節[6-9]。目前，國內在該領域的研究和應用較少，僅少數高校開展過相關工作，研究在真實力環境中設計平動點軌道有一定的理論意義和較大的工程價值。

到目前為止，對傳統CRTBP的理論研究已較成熟，而針對平動點附近更精確動力學模型的研究及軌道設計仍存在挑戰。尤其是地月系統中，平動點附近動力學環境非常復雜，軌道穩定性較差，CRTBP的結果與實際軌道差別較大。如考慮的假設條件過多，設計的軌道精度不高，將嚴重影響航天器在軌飛行的時間和保持所需的燃料。本文針對上述問題，以地月系統共線平動點為背景，在不設置任何假設條件的條件下直接在慣性系中建立動力學模型，并采用軌道拼接和多步打靶法，調用星歷表數據設計更符合真實的地月空間動力學環境的擬周期軌道。

1 坐標系統

建立高精度星歷模型和數值仿真需使用坐標系4個：描述CRTBP的旋轉坐標系、J2000地球赤道面慣性坐標系、地心瞬時旋轉坐標系和平動點瞬時旋轉坐標系[10]。如圖1所示：CRTBP坐標系未單獨繪出，若假設月球公轉軌道為正圓形且轉速恒定，則該坐標系與平動點瞬時旋轉坐標系平行。具體定義如下。

a)CRTBP的旋轉坐標系O-xyz：假設月球繞地球公轉軌道為正圓，且軌道角速度恒定，以地球和月球公共質心為原點O，Ox軸從地心指向月心，Oy軸與月球繞地球公轉的角速度方向重合，Oz軸與Ox、Oz軸構成右手垂直坐標系。

b)J2000地球赤道面慣性坐標系OI-xIyIzI(J2000 GECS)：以地心為原點OI，J2000歷元時刻地球平赤道面為參考平面，OIxI軸從原點OI指向J2000歷元時刻的春分點，OIzI軸與參考平面法線正向一致，OIyI軸與OIxI、OIzI軸構成右手垂直坐標系。

c)地心瞬時坐標系OE-xEyEzE(GRCS)：以地心為原點OE，OExE軸從地心指向月心，OEzE軸與月球繞地球公轉的瞬時角速度方向重合，OEyE軸與OExE、OEzE軸構成右手垂直坐標系。

d)平動點(L2)瞬時旋轉坐標系OL-xLyLzL(LCRCS)：以所選取的平動點，如L2點為原點OL，OLxL、OLyL、OLzL軸分別與OE-xEyEzE系的OExE、OEyE、OEzE軸平行，僅原點移動到至瞬時平動點的位置。

2 星歷模型

目前在限制性三體問題的研究中，用星歷表數據可選擇的建模方式主要有在J2000慣性系中或在地心旋轉系中建立動力學模型兩種。因J2000系中的模型相對簡單，故本文用前者進行建模。J2000系中航天器m、地球、太陽，以及月球間的位置矢量如圖2所示。

`本文只考慮地球、月球、太陽引力和太陽光壓力的影響。因太陽系中各天體均圍繞太陽系質心(SSB)運動，故J2000坐標系中太陽系質心可視為慣性點，航天器相對地球(原點)的位置矢量可表示為

rp=rpssb-rEssb

(1)

式中：rpssb，rEssb分別為航天器和地心相對太陽系質心的位置矢量。本文所有天體位置信息均通過SPICE在MATLAB中的工具箱MICE從星歷表bsp文件DE421中提取。根據牛頓定律，航天器在J2000坐標系中的星歷模型方程為

(2)

rpS=rp+rES

(3)

rpM=rp+rEM

(4)

若考慮其他天體的影響，則在式(2)中添加相應的引力項即可。

3 CRTBP周期軌道設計

在設計高精度擬周期軌道時需以CRTBP的軌道作為迭代初值，本文先給出CRTBP周期軌道的設計方法。旋轉系中單位歸一化后平動點附近的航天器動力學方程為

(5)

(6)

(7)

式中：x，y，z為航天器的位置坐標；μ為月球與地月系統總質量的比值；

將該動力學方程寫成狀態空間的形式，有

(8)

(9)

式中：Df(x(t))為f(x)在x(t)處的Jacobi矩陣。設Φ(t,t0)=?x(t)/?x(t0)為式(8)描述的系統從時刻t0至時刻t的狀態轉移陣，則x(t)=Φ(t,t0)x(t0)，即狀態轉移矩陣反映兩個時刻狀態變量間的線性映射關系。Φ(t,t0)對時間求導可得

A(t)Φ(t,t0)

(10)

根據Φ(t0,t0)=I，給定初始條件x(t0)，用數值積分即可求得任意時刻的狀態轉移矩陣。

目前常用的CRTBP周期軌道數值設計方法為微分校正方法，其原理就是前時刻初始狀態到目標時刻進行數值積分得到末端時刻狀態與理想狀態的偏差，通過狀態轉移矩陣對初始狀態進行修正，使末端狀態的偏差不斷減小直至滿足約束條件。數值積分的初值可通過LP法等攝動法求得[2]。

設目標狀態為xdes(tf)，對x(t0)進行數值積分至tf,得到時刻tf的狀態x(tf)。x(tf)求變分可得

(11)

對式(8)、(10)進行1個周期的數值積分得?x(tf)/?x(t0)=Φ(tf,t0)，δx(tf)=xdes(tf)-x(tf)，再求解式(11)即可得初值的修正量δx(t0)和周期的修正量δtf。式(11)有未知量7個和方程6個，為減少計算量，可用軌道關于某一平面的對稱性質，如Halo軌道在旋轉坐標系內均關于平面XOZ對稱，即軌道穿越XOZ平面時，有

y=0

(12)

(13)

(14)

4 星歷模型的擬周期軌道

以上述CRTBP的周期軌道作為星歷模型中擬周期軌道設計的迭代初值，通過調用星歷表數據對軌道進行拼接可得所需的擬周期軌道。具體軌道設計步驟如下。

步驟1：選取星歷模型需設計的擬周期軌道的起始歷元時刻t0和終端歷元時刻tf，歷元時刻由世界時(UTC)給出并轉化成星歷時刻(ET)。起始時刻選定后，以該時刻地月間的距離為無量綱化的單位距離，在CRTBP模型中用微分校正方法生成相應的周期軌道。

步驟2：在歷元時刻t0，tf間選取N個拼接點對應時刻ti，i=0，…，N-1。本文用時間等間距方法選取，將離散時刻ti轉換至CRTBP模型中，求得對應的周期軌道上的狀態，再根據這些離散狀態對應的歷元時刻和相應的天體位置信息將其轉換至平動點瞬時旋轉坐標系中，構成星歷模型中的軌道節點狀態初值。

步驟3：將前一步得到的平動點瞬時旋轉坐標系中的節點狀態作為迭代初值，用固定時間多步打靶法進行位置和速度校正，得到N-1段軌道弧段。當這些弧段在拼接點處的位置誤差小于1×10-5km、速度誤差小于1×10-8km/s時，可認為實現了位置連續和速度連續。

5 CRTBP系至J2000系的轉換

(15)

(16)

(17)

歷元時刻t時地月旋轉坐標系相對J2000系的角速度為

(18)

(19)

(20)

式中：TEM為月球繞地球的公轉周期。根據式(19)、(20)可得從地月旋轉系至J2000系的轉換矩陣為

(21)

6 多步打靶法

在上述擬周期軌道設計中用多步打靶法替代單步微分校正進行迭代，這是因為與CRTBP模型不同，星歷模型中平動點附近存在的軌道為擬周期軌道，不具備x(t+T)=x(t)的性質，故單步微分修正的方法已不適用，需用多步打靶法對軌道上的各節點進行校正以獲得所求的軌道。

設軌道的使用壽命為T=tf-t0，將T等分為N個時間段，各時間段的節點分別為ti，i=0，1，2，…，N，其中tN=tf，時間節點處的狀態為x(ti)。設Δt=ti+1-ti，以x(ti)為初值積分Δt時間后的狀態為φ(x(ti))。當所有x(ti)均位于同一軌道上時，有

φ(x(ti))=x(ti+1)

(22)

由式(22)可得

φ(x(ti)+δx(ti))≈φ(x(ti))+Φ(ti+1,ti)δx(ti)=

x(ti+1)+δx(ti+1)

(23)

式(23)在N個時間段上均滿足。聯立N個方程組可得

(24)

式中：I為單位陣。令DF為-F在ΔX處的Jacobi矩陣，因矩陣DF行數小于列數，故該方程組有無數組解，其偽逆解為

δX=-(DF)T(DF·(DF)T)-1F

(25)

當N較大時，直接求解式(25)計算量會非常大，本文用文獻[11]的方法簡化求解過程。

設M=DF·(DF)T，引入變量Z=[Z1…ZN-1]T=M-1F，則式(25)變為

MZ=F

(26)

δX=-(DF)TZ

(27)

對M進行Cholesky分解，可得

(28)

式中：

D1=I+Φ(t1,t0)(Φ(t1,t0))T

Li=-Φ(ti,ti-1)(Di-1)-1

Di=I+Φ(ti,ti-1)(Φ(ti,ti-1))T-LiDi-1(Li)T

引入中間變量Y=[Y1…YN]T，則式(26)的解Z為

Y1=x(t1)-φ(x(t0))

Yk=x(tk)-φ(x(tk-1))-Lkδx(tk-1),

k=2,3,…,N-1

ZN-1=(DN-1)-1δx(tN-1)

Zk=(Dk)-1δx(tk)-(Lk+1)TZk+1,

k=2,3,…,N-1

再由式(27)即可求得各節點處的狀態修正量δx(ti)。

7 數值仿真

用本文的共線平動點附近擬周期軌道的構造和分析方法，研究星歷模型中地月系共線平動點附近的擬周期軌道。相關參數見表1。

用本文方法在地月系統L2點附近拼接的某條擬周期Halo軌道如圖4～6所示。軌道運行時間為2014年7月1日0時至2015年7月1日0時，即1年的時間，所用拼接點數為200，形成弧段199個。

表1 地月系平動點基本參數

其中作為迭代初始軌道的Halo軌道的Z向幅值約10 000 km。為更清晰地給出擬周期軌道的形狀、瞬時平動點的位置及航天器運行在該軌道上時與地球和月球的位置關系，并與CRTBP的軌道進行對比，分別在平動點瞬時旋轉系、地心瞬時旋轉系和J2000系中給出了設計的軌道。由圖4～6可知：在地心瞬時坐標系中，由于月球偏心率等因素的影響，地月間的距離隨時間變化，導致瞬時平動點沿X軸也呈現周期性變化，周期軌道的形狀與CRTBP中的軌道差別較大，隨平動點沿X軸作周期性擺動，但在平動點瞬時旋轉系中，周期軌道的形狀與CRTBP的軌道較為相近。

8 結束語

本文針對地月系統共線平動點附近動力學環境復雜，軌道穩定性差，對設計精度要求高的特點，考慮月球真實軌道以及地球、月球和太陽的影響，在慣性系中建立了地月系統平動點附近航天器的高精度星歷模型。針對在CRTBP設計的軌道精度不高的問題，以CRTBP的周期和擬周期軌道為初值，用軌道拼接方法和多步打靶法設計了星歷模型中的擬周期軌道，并用數值仿真對算法進行了驗證。仿真結果表明：本文設計的算法可有效獲得所需的軌道。

本文所研究內容可為平動點附近的航天器提供滿足真實動力學環境的軌道，有效節省軌道保持所需的燃料。

[1] FARQUHAR R W, KAMEL A A. Quasi-periodic orbits about the translunarlibration point[J]. Celestial Mechanics, 1973, 7(4): 458-473.

[2] RICHARDSON D L. Analytic construction of periodic orbits about the collinear points[J]. Celestial mechanics, 1980, 22(3): 241-253.

[3] BREAKWELL J V, BROWN J V. The Halo family of 3-dimensional periodic orbits in the Earth-moon restricted 3-body problem[J]. Celestial Mechanics, 1979, 20(4): 389-404.

[4] HOWELL K C. Three-dimensional periodic ‘Halo’ orbits[J]. Celestial Mechanics, 1984, 32(1): 53-71.

[5] HOWELL K C, Spencer D B. Periodic orbits in the restricted four-body problem[J]. Acta Astronautica, 1986, 13(8): 473-479.

[6] PAVLAK T, HOWELL K. Strategy for optimal, long-term stationkeeping of libration point orbits in the Earth-moon system[C]// AIAA/AAS Astrodynamics Specialist Conference. Minneapolis: AIAA/AAS, 2012: 13-16.

[7] OZIMEK M T. Low-thrust trajectory design and optimization of lunar south pole coverage missions[D]. Purdue: Purdue University, 2010.

[8] WAWRZYNIAK G G. The dynamics and control of solar-sail spacecraft in displaced lunar orbits[D]. Purdue: Purdue University, 2011.

[9] PAVLAK T A. Trajectory design and orbit maintenance strategies in multi-body dynamical regimes[D]. Purdue: Purdue University, 2013.

[10] 錢霙婧. 地月空間擬周期軌道上航天器自主導航與軌道保持研究[D]. 哈爾濱: 哈爾濱工業大學, 2013.

Quasi-PeriodicOrbitDesignaroundEarth-MoonLibrationPointsUsingEphemerisModel

LIU Gang1, 2， HUANG Jing1, 2， GUO Si-yan1, 2

(1. Shanghai Key Laboratory of Aerospace Intelligent Control Technology, Shanghai 201109, China;2. Shanghai Institute of Spacecraft Control Technology, Shanghai 201109, China)

In order to improve the insufficient of lacking perturbation analysis during the orbit design process by using the circular restricted three-body model and increase the accuracy of the orbit design, the method of quasi-periodic orbit design around the Earth-moon libration points using ephemeris model was studied in this paper. Without any assumptions, the high-accuracy ephemeris model for spacecraft around the Earth-moon libration points was established in the J2000 inertial coordinate system with the consideration of real moon orbit and effect of the Earth, moon and sun. The periodic orbit in circular restricted three-body problem (CRTBP) was treated as the iterative initial value. The quasi-periodic orbit was obtained by orbit conjoint using ephemeris parameters. The single differentiation modification was taken placed by multiple-shooting method for iteration. The points in the orbit were corrected to obtain the quasi-periodic orbit. The design steps were given. The simulation results showed that the quasi-periodic orbit around libration points in the Earth-moon system would be obtained, which could provide the orbit that satisfied for real dynamic environment and save fuel for orbit keeping. The method has some application valuable in enginerring.

Earth-moon system; libration point; restricted three-body problem; solar radiation; gravity perturbation; quasi-periodic orbit; multiple shooting method; ephemeris

1006-1630(2017)05-0016-07

2016-11-19；

2016-12-15

劉 剛(1985—)，男，博士，主要研究方向為航天器軌道與姿態控制、航天器軌道設計、非線性控制、最優控制等。

V412.41

A

10.19328/j.cnki.1006-1630.2017.05.003