星歷模型地月系統(tǒng)平動點擬周期軌道設(shè)計研究

劉 剛，黃 靜，郭思巖

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

星歷模型地月系統(tǒng)平動點擬周期軌道設(shè)計研究

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

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

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

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

0 引言

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

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

1 坐標(biāo)系統(tǒng)

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

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

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

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

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

2 星歷模型

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

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

rp=rpssb-rEssb

(1)

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

(2)

rpS=rp+rES

(3)

rpM=rp+rEM

(4)

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

3 CRTBP周期軌道設(shè)計

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

(5)

(6)

(7)

式中：x，y，z為航天器的位置坐標(biāo)；μ為月球與地月系統(tǒng)總質(zhì)量的比值；

將該動力學(xué)方程寫成狀態(tài)空間的形式，有

(8)

(9)

式中：Df(x(t))為f(x)在x(t)處的Jacobi矩陣。設(shè)Φ(t,t0)=?x(t)/?x(t0)為式(8)描述的系統(tǒng)從時刻t0至?xí)r刻t的狀態(tài)轉(zhuǎn)移陣，則x(t)=Φ(t,t0)x(t0)，即狀態(tài)轉(zhuǎn)移矩陣反映兩個時刻狀態(tài)變量間的線性映射關(guān)系。Φ(t,t0)對時間求導(dǎo)可得

A(t)Φ(t,t0)

(10)

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

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

設(shè)目標(biāo)狀態(tài)為xdes(tf)，對x(t0)進行數(shù)值積分至tf,得到時刻tf的狀態(tài)x(tf)。x(tf)求變分可得

(11)

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

y=0

(12)

(13)

(14)

4 星歷模型的擬周期軌道

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

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

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

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

5 CRTBP系至J2000系的轉(zhuǎn)換

(15)

(16)

(17)

歷元時刻t時地月旋轉(zhuǎn)坐標(biāo)系相對J2000系的角速度為

(18)

(19)

(20)

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

(21)

6 多步打靶法

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

設(shè)軌道的使用壽命為T=tf-t0，將T等分為N個時間段，各時間段的節(jié)點分別為ti，i=0，1，2，…，N，其中tN=tf，時間節(jié)點處的狀態(tài)為x(ti)。設(shè)Δt=ti+1-ti，以x(ti)為初值積分Δt時間后的狀態(tài)為φ(x(ti))。當(dāng)所有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個時間段上均滿足。聯(lián)立N個方程組可得

(24)

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

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

(25)

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

設(shè)M=DF·(DF)T，引入變量Z=[Z1…ZN-1]T=M-1F，則式(25)變?yōu)?/p>

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)即可求得各節(jié)點處的狀態(tài)修正量δx(ti)。

7 數(shù)值仿真

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

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

表1 地月系平動點基本參數(shù)

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

8 結(jié)束語

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

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

[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ǎo)航與軌道保持研究[D]. 哈爾濱: 哈爾濱工業(yè)大學(xué), 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—)，男，博士，主要研究方向為航天器軌道與姿態(tài)控制、航天器軌道設(shè)計、非線性控制、最優(yōu)控制等。

V412.41

A

10.19328/j.cnki.1006-1630.2017.05.003