胡青龍

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

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

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

1 多目標(biāo)交會(huì)轉(zhuǎn)移軌道問題描述

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

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

飛行器在整個(gè)星際飛行中的總的速度增量ΔV包括探測(cè)器離開停泊軌道進(jìn)行星際航行需要的速度增量ΔVL、探測(cè)器與到達(dá)星體交會(huì)時(shí)所需的制動(dòng)速度增量ΔVα、探測(cè)器飛越第i顆目標(biāo)星時(shí)需要增加的深空機(jī)動(dòng)補(bǔ)償ΔVmi等3個(gè)部分，則飛行器在整個(gè)星際飛行中的總的速度增量ΔV可以表示為：

(1)

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

(2)

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

(3)

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

探測(cè)器與到達(dá)星體交會(huì)時(shí)所需的制動(dòng)速度增量ΔVα的表達(dá)式為：

(4)

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

探測(cè)器飛越第i顆目標(biāo)星時(shí)需要增加的深空機(jī)動(dòng)補(bǔ)償ΔVmi為：

(5)

上式中，Vmi+和Vmi-分別為深空機(jī)動(dòng)前后的日心速度。

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

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

(1)能量匹配約束

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

V∞1=Vs--Vps

(6)

V∞2=Vs+-Vps

(7)

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

∈R)

(8)

(2)最小飛越高度約束

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

e=1/sin(δ/2)

(9)

則此時(shí)雙曲線軌道的近心點(diǎn)半徑rp為：

(10)

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

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

(11)

2 多目標(biāo)交會(huì)轉(zhuǎn)移軌道的設(shè)計(jì)方法

2.1 交會(huì)目標(biāo)序列選擇

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

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

(12)

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

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

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

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

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

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

(13)

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

(14)

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

(15)

(16)

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

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

(17)

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

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

(18)

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

(19)

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

(20)

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

(21)

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

(22)

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

(23)

同理，總的速度增量ΔV相對(duì)于第i次深空機(jī)動(dòng)位置Rmi的偏導(dǎo)數(shù)為：

(24)

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

(25)

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

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

X(k+1)=X(k)+h▽?duì)?X(k))

(26)

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

(27)

3 算例與分析

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

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

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

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

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

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

表2 本文設(shè)計(jì)的ROSETTA任務(wù)深空機(jī)動(dòng)參數(shù)

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

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

4 結(jié)束語

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

參 考 文 獻(xiàn)

[1] 徐明.平動(dòng)點(diǎn)軌道的動(dòng)力學(xué)與控制研究綜述[J].宇航學(xué)報(bào),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.