胡青龍

西昌學院工程技術(shù)學院，四川西昌 615000

目前，世界主要軍事大國均把星際探測任務作為本國航天發(fā)展的重點。在星際探測任務中，任務的實現(xiàn)通常需要大量的能量，利用常規(guī)推進裝置采用直接轉(zhuǎn)移軌道的方法難以完成探測任務。因此，既要節(jié)省能量又要便于技術(shù)實現(xiàn)的星際多目標轉(zhuǎn)移軌道方法成為各國深空探測研究的熱點[1-3]。當前倍受關(guān)注的深空探測軌道轉(zhuǎn)移技術(shù)有借力飛行技術(shù)[4]、小推力轉(zhuǎn)移技術(shù)[5]和太陽帆轉(zhuǎn)移技術(shù)[6]。小推力轉(zhuǎn)移技術(shù)雖然已在美國的DS-1任務、日本的MUSES-C任務以及歐空局的SMART-I任務中得到了試驗和應用，但是對長達數(shù)年乃至數(shù)十年的遠景星際探測任務而言，高效、超長壽命的小推力發(fā)動機的研發(fā)仍是困擾任務設(shè)計人員的一大難題。而具有光明前景的太陽帆轉(zhuǎn)移技術(shù)，目前在技術(shù)實施上還相對困難。與小推力轉(zhuǎn)移技術(shù)和太陽帆轉(zhuǎn)移技術(shù)相比，借力飛行技術(shù)不但可以有效降低探測任務所需的發(fā)射能量和總的速度增量，而且具有很好的工程實施性，已成為深空星際探測任務的基本手段。

本文針對星際探測任務中多目標交會轉(zhuǎn)移軌道設(shè)計問題，利用借力飛行技術(shù)，結(jié)合Tisserand原理[7]、Pork-Chop圖法[8]，提出一種多目標交會的混合設(shè)計方法。采用P-rp曲線[9]確定借力目標的序列，利用Pork-Chop圖確定設(shè)計參數(shù)的可行域和時序，避免了傳統(tǒng)方法對借力交會目標序列和初始軌道段設(shè)計不收斂的問題，通過“軟匹配”策略尋找最優(yōu)設(shè)計參數(shù)。最后，以國際上經(jīng)典的小天體探測任務(ROSETTA任務)為例，對其轉(zhuǎn)移軌道的初始方案進行設(shè)計和分析，驗證該方法的正確性和有效性。

1 多目標交會轉(zhuǎn)移軌道問題描述

1.1 轉(zhuǎn)移軌道設(shè)計的目標函數(shù)

對于一個多目標交會的星際探測任務而言，轉(zhuǎn)移軌道設(shè)計的實質(zhì)是搜索一條既能滿足各種約束條件，又能使得任務所需的燃料消耗量總和最小的轉(zhuǎn)移軌道。在轉(zhuǎn)移軌道計算中，燃料消耗通常用速度增量表示。假設(shè)探測器從出發(fā)星體發(fā)射，飛越n個目標星后，與到達星體交會；在飛越第i個目標星時，若能量不能完全匹配，則增加深空機動補償，那么多天體交會發(fā)射機會搜索的目標函數(shù)應為飛行器所消耗的能量最小，即飛行器在整個星際飛行中的總的速度增量最小。

飛行器在整個星際飛行中的總的速度增量ΔV包括探測器離開停泊軌道進行星際航行需要的速度增量ΔVL、探測器與到達星體交會時所需的制動速度增量ΔVα、探測器飛越第i顆目標星時需要增加的深空機動補償ΔVmi等3個部分，則飛行器在整個星際飛行中的總的速度增量ΔV可以表示為：

(1)

多天體交會發(fā)射機會搜索的目標函數(shù)則為：

(2)

這里，X=[VL,Vα,VA,Vl∞,Vlp,Va∞,Vap,Vm1+,Vm1-,…,Vmi+,Vmi-,…]T∈(R2k+8)，其中探測器離開停泊軌道進行星際航行需要的速度增量ΔVL的表達式如下：

(3)

上式中，VL為探測器日心引力場的出發(fā)速度，VEL為發(fā)射天體繞太陽的公轉(zhuǎn)速度，Ve∞為逃逸發(fā)射天體引力場所需要的速度，VLp為探測器在停泊軌道上運行的速度。

探測器與到達星體交會時所需的制動速度增量ΔVα的表達式為：

(4)

上式中，Vα為探測器日心引力場的到達速度，VA為到達星體繞太陽公轉(zhuǎn)的速度，Va∞為逃逸到達星體引力場所需的速度，Vap為探測器在繞飛軌道上運行的速度。

探測器飛越第i顆目標星時需要增加的深空機動補償ΔVmi為：

(5)

上式中，Vmi+和Vmi-分別為深空機動前后的日心速度。

1.2 轉(zhuǎn)移軌道設(shè)計的約束條件

在星際探測多目標交會轉(zhuǎn)移軌道設(shè)計中，除了固有的星歷約束外，主要考慮飛越目標星體時的能量匹配約束和最小飛越高度約束。

(1)能量匹配約束

由于探測器在行星重力場內(nèi)飛行的時間比在日心參考系內(nèi)飛行的時間小很多，所以假設(shè)探測器飛入和飛出行星重力場是在同一時刻，行星的位置和速度在飛越過程中不改變。探測器在飛越行星時距離行星質(zhì)心的距離遠遠小于行星到太陽的距離，可以認為探測器在飛越時刻日心參考系的位置矢量Rs與行星在飛越時刻的位置矢量Rps相等。所以，探測器飛入和飛出行星重力場的雙曲線逃逸速度V∞1和V∞2可表示為：

V∞1=Vs--Vps

(6)

V∞2=Vs+-Vps

(7)

其中，Vs-和Vs+分別為飛越前后探測器相對于日心的速度。

∈R)

(8)

(2)最小飛越高度約束

為了避免探測器在飛越過程中飛入行星大氣(氣動-借力除外)或與行星相撞，所以對飛越的高度也提出要求。設(shè)曲線逃逸速度矢量轉(zhuǎn)過的角度為δ，根據(jù)雙曲線軌道的基本關(guān)系可以計算得到飛越時雙曲線軌道的偏心率e：

e=1/sin(δ/2)

(9)

則此時雙曲線軌道的近心點半徑rp為：

(10)

用Hmin表示最小的飛越高度，則應滿足的最小飛越高度約束為：

χ(x)=rp-Hmin≥0 (χ∈R)

(11)

2 多目標交會轉(zhuǎn)移軌道的設(shè)計方法

2.1 交會目標序列選擇

采用基于Tisserand原理的能量P-rp圖法選擇借力飛行的路徑和可能交會目標的序列。P-rp圖法是采用與能量等高線圖有關(guān)的Tisserand原理描述借力飛行軌道的方法。Tisserand原理中的常數(shù)T定義為：

T=rP/a+2[a(1-e2)/rP]1/2cosi

(12)

其中，rP為行星的平均軌道半徑。由Tisserand原理、圓錐曲線拼接方法[10]和借力飛行的機理可以得到P-rp圖。采用P-rp曲線圖方法可以確定借力飛行的序列和基本的能量要求。

2.2 設(shè)計參數(shù)的確定與選擇區(qū)域

星際探測多目標交會轉(zhuǎn)移軌道設(shè)計參數(shù)的選擇區(qū)域，可通過Pork-chop圖法或等高線圖獲得，這里不再贅述。根據(jù)Lambert定理，如果始末端的位置矢量R1和R2、飛行時間和飛行方向已知，則飛行軌道就可以確定[11]。若已知發(fā)射和飛越時間，通過行星星歷的計算可以得到發(fā)射和飛越時刻行星的位置，求解Lambert問題就可以得到初始和末端的速度矢量V1，V2，從而確定該段軌道，由此可見時間是確定軌道的關(guān)鍵參數(shù)。對于整個多天體交會轉(zhuǎn)移軌道而言，若已知發(fā)射時間、飛越時間和到達時間，通過行星星歷的計算和Lambert問題的求解就可以確定各軌道段。因此，這里選擇時間作為設(shè)計尋優(yōu)的變量。

2.3 轉(zhuǎn)移軌道的設(shè)計與優(yōu)化算法

多目標交會轉(zhuǎn)移軌道設(shè)計問題可歸結(jié)為一種非線性多約束多變量搜索尋優(yōu)問題。對于該問題的初始參數(shù)可以通過群體搜索策略和種群中個體之間的信息交換確定，而精確的設(shè)計與搜索可通過梯度下降法得到。對于一個非線性多維參數(shù)優(yōu)化問題，如何求解性能指標對自由變量的偏導數(shù)是一個困難。本文基于變分和主矢量原理，對性能指標相對于自由變量偏導數(shù)的解析形式進行推導，從而使轉(zhuǎn)移軌道優(yōu)化問題得到簡化。

假設(shè)X(t)是相對于飛行器某段軌道的狀態(tài)軌線，定義M=X(t0)和N=X(t1)，假設(shè)狀態(tài)軌線有微小的改變δX(t)，那么泛函數(shù)M,N也會隨之改變。δX0，δX1分別為M,N的變分。狀態(tài)轉(zhuǎn)移矩陣可以描述初始狀態(tài)微小變化和終止狀態(tài)微小變化之間的關(guān)系，即M,N的變分之間的關(guān)系。對于多天體交會的飛行軌道而言，整個軌道可以分為若干個軌道段，每一段軌道的狀態(tài)確定后，該段軌道的狀態(tài)轉(zhuǎn)移矩陣就可以確定。為了便于討論，這里定義4個3×3的矩陣Φrr，Φrv，Φvr和Φvv，即

(13)

在軌道弧段末端點的位置狀態(tài)變量δX2、速度狀態(tài)變量δV2與初始點的位置狀態(tài)變量δX1、速度狀態(tài)變量δV1之間的關(guān)系為：

(14)

將(13)式整理，可以推導出δV1，δV2與δX1，δX2之間的關(guān)系，即

(15)

(16)

如果端點是深空機動點或端點受到行星星歷的約束，則上面的方程(16)必須求解帶有端點位置固定的約束，即：

dX(t)=δX(t)+V(t)δt

(17)

由于端點位置固定則有dX(t)=0，故方程(17)可簡化為：

δX(t)=-V(t)δt

(18)

聯(lián)合求解方程(16)和(18)，可以得到軌道段初始速度V1，末端速度V2與初始點時間t1和末端點時間t2之間的關(guān)系為：

(19)

如果端點是借力天體，則V1和V2為雙曲線超速；如果端點是深空機動，則V1和V2為飛行器的日心速度。由方程(19)式可以推導如下：

(20)

總的速度增量ΔV相對于發(fā)射時間tL的偏導數(shù)為：

(21)

假設(shè)轉(zhuǎn)移軌道可以分為n段，則總的速度增量ΔV相對于到達時間ta的偏導數(shù)為：

(22)

假設(shè)第i次深空機動前的軌道段為k，深空機動后的軌道段為k+1，則總的速度增量ΔV相對于第i次深空機動時間tmi的偏導數(shù)為：

(23)

同理，總的速度增量ΔV相對于第i次深空機動位置Rmi的偏導數(shù)為：

(24)

由于在優(yōu)化過程中，借力飛行的時間不作為自由變量，而是作為滿足借力飛行匹配條件的變量，假設(shè)第i次借力飛行，借力飛行前的軌道為第q段，借力飛行后的軌道為q+1段，所以這里有：

(25)

綜上，可將復雜的軌道優(yōu)化問題轉(zhuǎn)化為一個多維無約束的參數(shù)優(yōu)化問題，這里選用梯度下降法，選取目標函數(shù)的負梯度方向作為每步迭代的搜索方向，逐步逼近函數(shù)的極小值點。

利用梯度下降的迭代公式，對于第k+1次迭代，則有

X(k+1)=X(k)+h▽φ(X(k))

(26)

這里從X(k)出發(fā)，沿梯度方向，取步長參數(shù)h，下一步可到達點X(k+1)，其中φ為自由變量X的非線性函數(shù)，則φ的梯度定義為：

(27)

3 算例與分析

下面以國際上經(jīng)典的小天體探測任務(ROSETTA任務)為例，結(jié)合其約束條件(參見文獻[12])，對其轉(zhuǎn)移軌道的初始方案進行設(shè)計和分析，以驗證本文設(shè)計方法的有效性。

Churyumov-Gerasimenko彗星的近日點為1.24AU，遠日點為5.68AU，軌道傾角為7°。利用本文的多目標交會轉(zhuǎn)移軌道設(shè)計方法，設(shè)計出ROSETTA任務的轉(zhuǎn)移軌道包括3次地球借力飛行和1次火星借力飛行，飛越了2顆小天體，進行了5次較大的深空機動，最終與Churyumov-Gerasimenko彗星實現(xiàn)交會。任務飛行軌跡如圖1所示。

圖1 本文設(shè)計出的ROSETTA任務的轉(zhuǎn)移軌道

本文設(shè)計出任務發(fā)射和借力飛行時的雙曲線超速參數(shù)如表1所示。從表1可以看出，本文設(shè)計的ROSETTA任務結(jié)果與歐空局公布的結(jié)果相比，發(fā)射和借力飛行時的雙曲線超速偏差分別為0.003km/s，0.002km/s，0.011km/s，0.007km/s和0.004km/s。

表1 本文設(shè)計的ROSETTA任務雙曲線超速參數(shù)

設(shè)計出ROSETTA任務的深空機動參數(shù)如表2所示。從表2可以看出，本文設(shè)計的ROSETTA任務結(jié)果與歐空局公布的結(jié)果相比，深空機動和交會時所需速度增量的偏差分別為0.035km/s，0.042km/s，0.005km/s，0.0066km/s，0.043km/s和0.002km/s。

表2 本文設(shè)計的ROSETTA任務深空機動參數(shù)

設(shè)計出任務借力飛越高度結(jié)果如表3所示。在表3中，借力飛行高度均為探測器距離借力天體質(zhì)心的高度，實際距離借力天體表面的高度應減去借力天體的半徑。本文設(shè)計的ROSETTA任務結(jié)果與歐空局公布的結(jié)果相比，借力飛行高度的偏差分別為3.66%，1.73%，6.63%和6.09%。

表3 本文設(shè)計的ROSETTA任務借力飛行飛越高度參數(shù)

4 結(jié)束語

針對深空探測中的轉(zhuǎn)移軌道問題，提出一種星際多目標交會轉(zhuǎn)移軌道設(shè)計方法。該方法基于Tisserand原理，采用能量曲線確定交會目標的序列，由Pork-Chop圖確定設(shè)計參數(shù)的可行域和時序，避免了傳統(tǒng)方法對交會目標序列和初始軌道段的假設(shè)，通過優(yōu)化算法尋找最優(yōu)的設(shè)計參數(shù)。最后，將本文的設(shè)計方法用于解決歐空局的ROSETTA任務深空轉(zhuǎn)移軌道的設(shè)計，設(shè)計結(jié)果與歐空局公布的結(jié)果一致，從而驗證了該設(shè)計方法的可行性和正確性。

參 考 文 獻

[1] 徐明.平動點軌道的動力學與控制研究綜述[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.)

[2] Koom W S,Lo M W,Marsden J E,Ross S D. Dynamical Systems, the Three-body Problem and Space Mission Design[M].Springer-Verlag New York Inc.,2007.

[3] Bookless J,McInnes C.Control of Lagrange Point Orbits Using Solar Sail Propulsion [J].Acta Astronautica,2008,6(2):159-176.

[4] XU Ming,XU Shi-jie.Structure-preserving Stabilization for Hamiltonian System and Its Applications in Solar Sail[J].Journal of Guidance,Control and Dynamics,2009,32(3):997-1004.

[5] Mehradad Saif, Weitian Chen, Qing Wu.High Order Sliding Mode Observers and Differentiators Application to Fault Diagnosis Problem[J]. Lecture Notes on Control and Information Science，2008,37(5):321-344.

[6] Qing Wu,Mehrdad Saif. Robust Fault Detection and Diagnosis for a Multiple Satellite Formation Flying System Using Second Order Sliding Mode and Wavelet Networks[C]//Proceeding of the 2007 American Control Conference, New York,IEEE,2007:426-431.

[7] Waters T J,McInnes C R. Periodic Orbits Above the Ecliptic Plane in the Solar Sail Restricted 3-body Problem [J].Journal of Guidance,Control and Dynamics,2007,30(3):687-693.

[8] LIU Tao, HE Zhao-wei, ZHAO Yu-shan. Continuous-thrust Orbit Maneuver Optimization Using Modified Robust Algorithm [J].Chinese Journal of Astronautics,2008,29(4):1216-1221.

[9] Mayer H, Carter T. Orbits and Relative Motion in the Gravitational Field of an Oblate Body [J]. Journal of Guidance, Control and Dynamics,2008,31(3):522-532.

[10] Matthew C,VanDyke,Christopher D Hall. Decentralized Coordinated Attitude Control within a Formation of Spacecraft [J].Journal of Guidance, Control and Dynamics,2006,29(5):1101-1109.

[11] Bevilacqua R, Romano M.Rendezvous Maneuvers of Multiple Spacecraft Using Differential Drag under J2Perturbation [J].Journal of Guidance, Control and Dynamics,2008,31(6):1595-1607.

[12] Biele J,Ulamec S.Capabilities of Philae,the Rosetta Lander [J]. Space Science Reviews, 2008,138 (1):275-289.