耿光有，王 玨，張志國，王建明，田繼超

(1.北京航空航天大學宇航學院，北京 100191；2.北京宇航系統(tǒng)工程研究所，北京 100076；3.中國運載火箭技術研究院，北京 100076)

0 引 言

運載火箭發(fā)射火星探測器等深空任務時，通常需要先進入低高度地球停泊軌道，待火箭滑行至預定位置后再次點火，從而將探測器送入奔赴火星的轉移軌道，新一代低溫長征火箭從海南發(fā)射火星探測器時，除了受到航落區(qū)安全及地面測控等因素的限制外，還有低溫入軌級兩次起動間最長允許滑行時間的嚴格限制[1-2]，使入軌點位置(或近地點幅角)的調節(jié)范圍較為有限，為滿足2～3周的有效發(fā)射日期窗口[3]，通常解決方法是加大射向(對應出發(fā)軌道傾角)或增加出發(fā)點速度能量，但加大射向會受到航落區(qū)限制，增加出發(fā)點速度能量一般也將增加探測器抵達火星時需要的制動速度增量，因此工程實現效果并不十分理想；而成功解決上述問題獲得更多優(yōu)質發(fā)射機會，是促成火星等深空探測工程順利實施的重要衡量指標；如果能在探測器從地球至火星耗時數月的轉移軌道途中，通過施加一合適的速度增量以解決上述問題，工程意義將是重大的；本文研究化學推進的火星探測器深空機動指可以等效為速度脈沖增量的中途機動。Fimple[4]曾提出在深空軌道與黃道面的升(或降)交點處，增加一次軌道面機動可以減小總速度增量，耿長福[5]、戴光明[6]介紹并推廣了這一思想，但考慮多天體攝動等影響因素后，最優(yōu)深空機動點位置可能已偏離軌道升(或降)交點很遠，而由于機動點具體位置和時刻一般相互獨立，因此很難用簡單的一維搜索予以解決；注意到深空機動的位置完全可能大幅度變化，因此在獲得一個良好的猜測初值前，綜合機動點位置與時刻四維搜索需要的計算量一般難以被接受。Lawden[7]提出主矢量法用于轉移軌道機動分析，Lion等[8]據此分析了非最優(yōu)軌道，文獻[9-10]推廣了應用，Conway[11]指出當不滿足Lawden必要條件時，通過增加深空機動是有用的，Iorfida等[12]采用極坐標開展了軌道面內、外中途修正(或深空機動)的進一步研究；Glandorf[13]將主矢量法拓展到圓錐曲線[14]軌道拼接，Navagh[15]、Olympio等[16]采用主矢量法分析了多天體借力飛行時深空機動的優(yōu)化問題，Olympio等給出了利用深空機動減小行星借力飛行中總速度增量的算例，并指出地球到火星借力飛行中一般僅需要一次深空機動，只有與黃道面近乎垂直的轉移軌道需要兩次深空機動，基本不需要兩次以上深空機動；喬棟等[17]采用主矢量法分析了深空多脈沖機動下的發(fā)射機會搜索，沈紅新[18]對脈沖推力最優(yōu)軌跡的Hamilton邊值統(tǒng)一性問題進行了分析，潘迅等[19]將主矢量法應用于月地平動點雙脈沖轉移軌道的尋優(yōu)快速確認中；但上述文獻均沒有行星引力影響球內飛行段的詳細分析研究，而對于長征運載火箭發(fā)射火星探測任務來講，深空機動能否降低運載火箭出發(fā)與探測器抵達火星的總速度增量，能否通過深空機動解決受運載火箭最長允許滑行時間與航落區(qū)限制下的發(fā)射日期窗口拓展等工程設計問題，亟需深入開展研究。

綜上，結合工程實際，本文采用主矢量法結合序列二次規(guī)劃[20]尋優(yōu)算法，完成了包括運載火箭地面發(fā)射起飛至探測器抵達近火點目標軌道，通過火星探測器轉移軌道深空中途增加一次速度脈沖機動(以下稱深空機動)下的系統(tǒng)優(yōu)化研究。

1 分析模型及算法

當確定了運載火箭的發(fā)射日期與探測器抵達火星日期，一般即可以開展無深空機動下，運載火箭發(fā)射火星探測轉移軌道的優(yōu)化分析計算。

1.1 火星探測發(fā)射軌道設計基本力學模型

采用火星探測器近火點目標軌道參數，由運載火箭完成抵達近火點的直接轉移軌道詳細設計分析可參見本人新近另外著述等，限于篇幅，此處僅給出扼要介紹。

在起飛時刻地心慣性系下，運載火箭穿越大氣飛行段的質心動力學矢量方程[21-22]為：

(1)

(2)

從探測器分離開始，整個地火轉移軌道都在日心J2000坐標系下完成積分，動力學方程如下：

(3)

采用雙向微分修正算法進行計算，其思路見文獻[23]，此算法與文獻[25]的直接打靶法有相通之處，優(yōu)點是便于增加深空速度脈沖后的深入優(yōu)化。

(4)

在近地出發(fā)點，存在3個約束：

(5)

在火星進入段，雙曲線軌道的三個約束：

(6)

正反雙向積分的相遇歷元時刻為Tm時：

(7)

式中:上標“+”表示正向積分軌道，上標“-”表示反向積分軌道，分析表明，連接點可以較大自由度選擇，初值不妨取Tm=(T1+T2)/2，即中間歷元時刻。

在方程(5)～(7)中給出了12項約束，采用牛頓迭代求解雅可比矩陣就可以完成運載火箭火星探測發(fā)射軌道分析設計。

1.2 主矢量算法猜測深空機動最優(yōu)初值

1.2.1主矢量理論及特性

Bryson和Ho[26]給出了動力學方程(8)下的哈密頓方程(9)：

(8)

式中：X為狀態(tài)量，rt為t時刻位置矢量，Vt為t時刻速度矢量，g(rt)為引力加速度，Γ為脈沖推力加速度，u(t)為推力的單位方向矢量。

(9)

式中：λr，λv為協(xié)態(tài)變量。

對于脈沖推力情況，增加一次深空機動后，成本函數為沿轉移軌道脈沖速度增量累計和[10-11]：

(10)

根據極小值原理，方程(9)的伴隨方程為：

(11)

λ(t)?-λv(t)

(12)

(13)

綜合式(11)～(13)，有：

(14)

對于最優(yōu)轉移軌道，主矢量滿足以下四個必要條件[7,10]：

(1)主矢量與它的微分處處連續(xù)；

(2)轉移期間主矢量模不超過1；

(3)如果存在中間推力脈沖時刻，主矢量方向與推力一致，其模為1；

(4)所有中間推力脈沖時刻的主矢量導數為0。

通過沿轉移軌道積分獲取λ(t)曲線，可以快速判斷轉移軌道是否最優(yōu)，按照條件(1)～(2)判定是否需要增加中途機動完成發(fā)射轉移軌道進一步優(yōu)化，對符合條件的采用主矢量法可快速完成最優(yōu)機動初值猜測，因此主矢量法已逐漸演化成一種重要的軌道中途機動優(yōu)化方法。

根據主矢量必要條件可以給出主矢量邊界條件[11]，定義主矢量初值及終態(tài)分別見式(15)、(16)：

(15)

(16)

式中:下標0表示起始狀態(tài)，f表示終端狀態(tài)，下同。

由方程(14)可得：

(17)

由式(17)可知，

(18)

1.2.2分析參考軌道主矢量曲線及轉移矩陣

依據主矢量特性，采用式(19)～(21)沿參考軌道得狀態(tài)轉移矩陣φ(t,t0)。

φ(t0,t0)=I

(19)

(20)

(21)

根據主矢量定義，結合式(15)～(18)，采用式(22)得λ(t)曲線。

(22)

1.2.3增加一次脈沖機動后最優(yōu)初值猜測

對于兩脈沖轉移軌道，成本函數為：

(23)

采用Jezewski[10]的思路，當增加一次中途脈沖機動后，成本函數變?yōu)椋?/p>

(24)

式中:右上角“+”表示脈沖機動后，“-”表示脈沖機動前，下同。

dJ=J1-J0

(25)

結合伴隨方程式(26)：

(26)

當起止端位置固定，在一階攝動下，可得式(27)：

(27)

將式(25)采用冪級數展開至二階項，再結合式(26)，寫成c的二次方程，并取?(dJ)/?c=0，得攝動速度c：

c=

(28)

(29)

(30)

(31)

此時，脈沖機動后最優(yōu)軌道所需位置rdsm與原指定參考軌道位置rnom和攝動量?rm見式(32)～(33)：

rdsm=rnom(tm)+?rm

(32)

(33)

對于工程中感興趣的地火直接轉移軌道，一般可以通過式(28)～(33)獲得深空機動最優(yōu)猜測初值。

1.2.4多中心天體引力下狀態(tài)轉移矩陣

由于工程中需要考慮地球或火星引力影響球內飛行段，經推導分析，對地球至火星轉移軌道，采用主矢量算法需要補充以下內容。

引力梯度矩陣式(21)需要修正為式(34)，式中ri為運載火箭或探測器相對地球(i=E)或火星(i=M)的位置矢量：

(34)

數值計算情況下，穿越影響球邊界時的6×6階狀態(tài)轉移矩陣：

(35)

式(35)對引力影響球邊界具體尺度不是特別敏感，以地球為例，設置距地心1000000 km處，注意到此式引起的主矢量一階導數不連續(xù)是由于引力中心坐標轉換引起，因此與Lawden主矢量特性必要條件并不抵觸；但新情況下需要深入具體分析。

當飛出地球影響球時，設t=t1，此時：

(36)

當飛入火星影響球時，設t=t2，此時：

(37)

其中,

(38)

式(36)～(38)中：P*,Q*分別指探測器在地球(或火星)引力邊界處相對日心的位置及速度；式(38)中“+”表示離開影響球，“-”表示進入影響球；式中右上角標“+”表示過影響球邊界后，“-”表示之前。

這樣，從地球停泊軌道出發(fā)到進入環(huán)火星目標軌道，涉及到的總過渡轉移矩陣為：

ψ(tf,t0)=φ(tf,t2)w(t2)φ(t2,t1)w(t1)φ(t1,t0)

(39)

與式(39)相匹配時，式(24)代表的累計速度指相對于地球停泊軌道出發(fā)速度增量、抵達火星目標軌道速度增量和中途深空機動脈沖速度增量?；诋斍拔墨I中多以忽略引力影響球內飛行段，即以引力影響球邊界處速度增量為關注點，為便于對比，需將式(24)變?yōu)槭?40)：

(40)

式中:下標E0指相對于地球中心，下標Mf指相對于火星中心。

1.2.5起始、抵達時刻與深空機動時刻的調整優(yōu)化

對于滿足主矢量必要條件，但總速度增量仍明顯高于所在年份最低速度增量[27]的轉移軌道，說明增加深空機動難以進一步降低成本函數，此時需要進一步優(yōu)化起始與抵達時刻，分析如下。

對于包含一次深空脈沖機動下，t時刻狀態(tài)量的Lambert求解問題：

(41)

在數值尋優(yōu)狀態(tài)下，成本攝動函數表示為：

(42)

式中:下標m表示脈沖機動狀態(tài)，下同。

對于從初始軌道轉移至目標軌道的一般問題，根據式(9)哈密頓函數，結合主矢量定義及伴隨方程式(26)，在滑行段，有：

(43)

沿著Jezewski[10]的思路，對于地球到火星轉移軌道，速度脈沖作用下，得成本攝動函數式：

(44)

1.2.6采用主矢量猜測最優(yōu)初值算法步驟

采用主矢量算法獲得深空機動最優(yōu)猜測初值的算法步驟見圖1，圖中關于“需調整出發(fā)抵達日期”主要指式(44)描述的轉移軌道情況。

圖1 采用主矢量算法獲得深空機動最優(yōu)猜測Fig.1 Obtaining initial guess of DSM via primer vector

1.3 序列二次規(guī)劃數值尋優(yōu)算法

利用深空機動最優(yōu)猜測初值，進一步采用序列二次規(guī)劃(SQP)尋優(yōu)算法[6,20]完成數值優(yōu)化；其非線性規(guī)劃的描述為：

(45)

基于關注的工程目標關系，式中t1，t2分別表示地球與火星引力影響球邊界跨越時刻。

約束條件為：

(46)

式中:ΔVmB表示允許的深空機動最大速度值。當限定運載火箭在停泊軌道的最長滑行時間時，需增加：

Th≤ThB

(47)

式中:ThB表示為限定運載火箭在停泊軌道的最長滑行時間。

由于機動點速度與機動時刻相互獨立，故取優(yōu)化變量為：

(48)

因為不同機動時刻所需要的ΔVm變化幅度較大，所以tm變化后，軌道基準點需要根據上一時刻地火轉移軌道獲得：

(49)

式中:上標“′”表示軌道新時刻基準點。

2 仿真分析與算例

本文分析算例中，初始軌道以長征火箭從海南發(fā)射場起飛，設定一初始射向A0(對應相應軌道傾角)，停泊軌道暫設為高度200 km圓，探測器分離時真近點角約25°；探測器抵達近火點的目標軌道高度為500 km圓，近火軌道傾角93°，到達近火點成為圓軌道時優(yōu)化結束。

2.1 地球或火星逃逸速度與近地點或近火點速度

根據近地點或近火點與引力影響球半徑處速度關系：

(50)

式中：下標p表示近地點或近火點，∞表示引力影響球半徑處。

相對近行星點時刻圓軌道的速度增量為：

(51)

由式(51)可知，近地點(或近火點)速度增量與地球(或火星)引力邊界逃逸速度不是簡單的線性關系，故采用不同優(yōu)化目標J2或J1得到的結果將存在一定差異。

2.2 地球到火星探測軌道深空機動分析算例

以2020年7月23日出發(fā)，約340天轉移，2021年6月28日中午抵達近火點目標軌道的計算為例，主矢量曲線圖如下。

圖2 參考軌道的λ～t曲線Fig.2 History of primer vector on reference trajectory

從圖2可以看出，算例中從地球出發(fā)到抵達近火點軌道轉移過程中的主矢量峰值小于1，但結合主矢量定義及穿越引力影響球對主矢量參數曲線計算的影響，再從地球引力影響球邊界到火星引力影響球邊界的過程來看，很可能存在通過深空機動進一步降低成本函數的可行性，具體需要采用數值尋優(yōu)算法進一步分析確認。

圖曲線Fig.3 Primer history of on reference trajectory

從圖3結合式(44)來看，調整出發(fā)與抵達日期均可以使轉移軌道發(fā)射總能量進一步降低，尤其是調整出發(fā)日期可以更明顯降低轉移軌道發(fā)射總能量,篇幅所限,此分析從略;下面重點關注深空機動后的效果影響,在主矢量法最優(yōu)猜測初值基礎上，依據不同優(yōu)化目標及約束，采用SQP算法數值尋優(yōu)結果，具體優(yōu)化結果見表1。

表1中，J1=ΔVEp+ΔVm+ΔVMp，指相對地球及火星停泊軌道速度增量與深空機動速度三者累計和；J2=ΔVE0+ΔVm+ΔVM0，指地球及火星引力影響球邊界速度增量與深空機動速度三者累計和；表頭中“J1目標Th≤400 s”指采用J1目標進行優(yōu)化且限定運載火箭停泊軌道滑行時間限定為不超過400 s；表中計算約束除采用射向A0=107°，對應出發(fā)軌道傾角25.5°外，其它如近地點高度等均同前述。

從表1可以看出，在設定的近地點與近火點軌道約束，以及限制運載火箭停泊軌道最長滑行時間條件下：

1)最優(yōu)猜測初值給出的優(yōu)化時機在出發(fā)后第101.9天，具體設計軌道約束下，采用SQP算法數值尋優(yōu)結果為141.8天到149.8天，即采用主矢量法獲得了良好初值。

表1 采用SQP算法數值尋優(yōu)結果Table 1 The numerical optimal results using SQP

2)采用J2為優(yōu)化目標，可將初始軌道的J2=7562.8 m/s降低到J2=6758.0 m/s，降低了804.8 m/s，雖然J1僅降低了8.6 m/s，但總的優(yōu)化效果巨大。

3)采用J1為優(yōu)化目標，可將初始軌道的J1=6377.1 m/s降低到J1=6347.7 m/s，降低了29.4 m/s，J2也降低了514.4 m/s，故優(yōu)化效果明顯。

4)采用J1為優(yōu)化目標且同時限定火箭停泊軌道最長滑行時間不超過400 s，則可以優(yōu)化到J1=6367.5 m/s，相比初始軌道可以降低9.6 m/s，相比不約束狀態(tài)多消耗了19.8 m/s，即通過深空機動完全可以調整運載火箭發(fā)射探測器入軌點的位置，而且消耗的總發(fā)射能量在可接受范圍內，故通過深空機動有效拓展發(fā)射日期窗口的優(yōu)化效果明顯。

注意到上述優(yōu)化求解是在嚴格限定近地點與近火點高度、傾角與真近點角約束下獲得的，而大多數行星際借力飛行等任務，只需要近星點最低高度約束，因此一般更易獲得進一步優(yōu)化的結果。

此外，研究分析發(fā)現：對地球到火星轉移軌道，只要主矢量曲線形狀能如圖2所示，即從地球引力影響球邊界t1到火星影響球邊界t2間，主矢量曲線中有高于t1或t2處的峰值，則如果以J2為優(yōu)化目標，一般均能找到更優(yōu)目標；而同等情況若以J1為優(yōu)化目標，除非總速度增量已接近(調整出發(fā)與抵達日期后得到的)最優(yōu)解[27]外，則一般也總能通過深空機動找到更優(yōu)解(如算例所示)。通過上述分析表明：深空機動對于降低運載火箭出發(fā)與探測器抵達火星的總速度增量，解決運載火箭受最長滑行時間與航落區(qū)限制下的發(fā)射日期窗口拓展等問題，意義顯著，雖然考慮引力影響球內飛行過程的優(yōu)化求解是非常復雜的。由于進一步的具體分析是非常復雜的工程問題，限于篇幅等因素，這里不再贅述。

綜上分析，可以得出結論：

1)利用主矢量法判斷并獲取猜測初值、再進一步由序列二次規(guī)劃算法完成精確數值尋優(yōu)，可以很好地實現多中心天體引力下的軌道優(yōu)化。

2)采用J1、J2或J1且限定運載火箭停泊軌道最長滑行時間等不同優(yōu)化目標，所得到的最終優(yōu)化點存在差異，故具體優(yōu)化目標需要視工程情況而確定。

3)采用深空機動可以降低總速度增量，但要想大幅度降低近地點出發(fā)與抵達近火點的速度增量是困難的，原因與衛(wèi)星軌道近地點加速時更易提高飛行軌道能量本質原理一致。

4)采用深空機動，通過深空中途小幅速度機動，可以調整對運載火箭入軌點位置的需求，即深空機動可以擴展發(fā)射日期窗口機會。

3 結 論

文中給出的多中心天體引力下，利用主矢量法獲取深空機動最優(yōu)初值、再由序列二次規(guī)劃算法完成精確數值尋優(yōu)的方法，很好地解決了當前長征火箭發(fā)射火星探測器工程任務中軌道優(yōu)化時遇到的復雜約束條件下，有效拓展發(fā)射日期窗口和降低總發(fā)射能量需求的難題，研究表明：

1)深空機動可以用來實現最小總速度增量的飛行軌道優(yōu)化，為了獲得期望的最小變軌優(yōu)化速度，可以根據具體約束，綜合總速度增量J1與J2優(yōu)化目標，以更好解決工程中的火星探測發(fā)射軌道優(yōu)化、軌道中途修正和行星際借力飛行優(yōu)化等問題。

2)以總速度增量J1為優(yōu)化目標，有效地解決了運載火箭最長滑行時間限制下的火星探測轉移軌道優(yōu)化設計的問題，通過深空機動，確保了長征火箭2～3周的火星探測發(fā)射窗口機會，實現了預期工程目標。

工程應用表明此方法穩(wěn)定、可靠好用，除了可以用于火星探測的發(fā)射軌道設計外，本方法還可用于運載火箭的其它深空探測任務的發(fā)射軌道設計。