胡青龍

西昌學院工程技術學院，四川西昌 615000

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

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

1 多目標交會轉移軌道問題描述

1.1 轉移軌道設計的目標函數(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ā)射天體繞太陽的公轉速度，Ve∞為逃逸發(fā)射天體引力場所需要的速度，VLp為探測器在停泊軌道上運行的速度。

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

(4)

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

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

(5)

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

1.2 轉移軌道設計的約束條件

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

(1)能量匹配約束

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

V∞1=Vs--Vps

(6)

V∞2=Vs+-Vps

(7)

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

∈R)

(8)

(2)最小飛越高度約束

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

e=1/sin(δ/2)

(9)

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

(10)

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

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

(11)

2 多目標交會轉移軌道的設計方法

2.1 交會目標序列選擇

采用基于Tisserand原理的能量P-rp圖法選擇借力飛行的路徑和可能交會目標的序列。P-rp圖法是采用與能量等高線圖有關的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ù)的確定與選擇區(qū)域

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

2.3 轉移軌道的設計與優(yōu)化算法

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

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

(13)

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

(14)

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

(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之間的關系為：

(19)

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

(20)

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

(21)

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

(22)

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

(23)

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

(24)

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

(25)

綜上，可將復雜的軌道優(yōu)化問題轉化為一個多維無約束的參數(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 算例與分析

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

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

圖1 本文設計出的ROSETTA任務的轉移軌道

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

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

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

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

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

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

4 結束語

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

參 考 文 獻

[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.