龔宇蓮，何英姿，李毛毛，李克行

(1. 北京控制工程研究所，北京 100190；2. 空間智能控制技術重點實驗室，北京 100190)

0 引 言

航天器離軌制動，屬于變軌問題的研究范疇，它特指在軌航天器通過施加推力進行制動，進而返回或進入行星大氣層的軌道改變過程[1]。離軌制動需要選擇合適的軌道位置，施加恰當?shù)耐屏κ噶浚购教炱髯罱K滿足再入/進入點約束條件。

早期的離軌制動研究主要以沖量假設和二體軌道假設為前提[2]。文獻[2]研究了針對再入角需求選擇不同的制動速度及制動角策略。文獻[3]提出了一種求解離軌制動脈沖解析算法，并給出了確定最優(yōu)制動點位置的一種迭代算法。近年來，離軌研究開始考慮固定推力、變推力情況下的變軌問題，文獻[4]以制動過程攻角和側(cè)滑角的變化率為控制量，設計了燃料最優(yōu)的離軌制動方法，文獻[5]分別分析了固定推力和可調(diào)推力情況下的離軌優(yōu)化問題，文獻[4-5]都以偽譜法為基礎求解優(yōu)化問題。文獻[6]設計了一種組合制導方法，采用了燃料最省方案計算初始制動參數(shù)進行開環(huán)制動，而后采用相對運動控制的方式，跟蹤標稱軌跡實現(xiàn)高精度的末制導。文獻[7]以大升阻比飛行器為研究對象，將再入接口條件取為地心距、速度和速度傾角，采用有限推力離軌制動方式，研究了單次“推-滑”和多次“推-滑”的離軌制動問題。文獻[8-9]根據(jù)再入軌跡分析和落點經(jīng)緯度，估計再入軌道最優(yōu)升交點地理經(jīng)度。

需要確定的制動參數(shù)，通常包含制動點火位置、制動速度大小及方向，這些參數(shù)的確定一般是在地面規(guī)劃完成。為提高飛行器再入點精度，文獻[10]通過在地面大量仿真訓練BP神經(jīng)網(wǎng)絡建立起控制變量和返回特征參數(shù)的映射關系，改進了在軌快速預報算法。隨著飛控計算機性能的發(fā)展，以及對飛行器快速機動及實時性要求的提高，在軌自主制動方法的研究也引起了學者的廣泛興趣[11-14]。此外制動過程中的制導策略對制動精度的影響同樣十分顯著，阿波羅飛船采用的PEG方法，在行星探測器著陸問題中得到廣泛研究[15-16]。針對再入飛行器，文獻[17]中采用狀態(tài)空間攝動法將J2項攝動的影響引入制導律來保證再入角和再入點的精度。

以往的制動策略規(guī)劃研究中，往往忽略地球引力攝動對軌道的影響。一般基于二體假設和沖量假設，規(guī)劃制動脈沖、制動角及制動起始點,使得估算的制動參數(shù)難以滿足實際工程應用的高精度需求。為了保證再入點精度，往往需要引入末端修正等輔助手段。但實際工程應用中，制動發(fā)動機一般對于連續(xù)開關的時間間隔有約束，且軌道再入飛行器的軌控發(fā)動機一般只能提供單方向的推力，要滿足高精度修正還需通過姿態(tài)機動來確保推力矢量方向。

基于以上問題和已有方法的不足，本文考慮地球非球形引力以及有限推力的實際情況，提出一種具備在軌自主決策、高精度落點狀態(tài)控制的離軌制動方案。本方法的核心思路如下：

(1)首先將制動問題轉(zhuǎn)化為軌道控制問題，將再入點緯度、速度以及傾角等需求轉(zhuǎn)化為對離軌過渡軌道的半長軸、偏心率以及近地點幅角需求，從而便于利用軌道機動方程實現(xiàn)初值求解。

(2)在飛行器運行于離軌待命圓軌時利用機載計算機空閑時間，根據(jù)自主測定軌結果，通過數(shù)值遞推軌跡，迭代求解制動參數(shù)，使得軌控完成后目標軌道根數(shù)誤差滿足精度需求。

(3)制動開始后，根據(jù)瞬平根轉(zhuǎn)換的近似解析公式，預估再入點飛行狀態(tài)。當估計的再入點滿足約束條件時制動結束。

1 離軌制動問題

完整的軌道根數(shù)包含了對再入點軌道面內(nèi)與軌道面外的所有約束，軌道面外的約束往往需要通過回歸軌道的設計[19]或軌道相位的調(diào)整來使得再入點與離軌待命圓軌共面。軌道面外的調(diào)整問題需要根據(jù)燃料與時間的約束進行規(guī)劃，離軌階段一般僅考慮軌道面內(nèi)的制動問題。因此，對制動點的約束可以簡化為對再入點緯度、再入點速度及速度傾角的約束。再入點地心緯度可以表達為如下：

φE=arcsin[sin(iE)sin(uE)]

(1)

式中：uE為再入點緯度幅角，iE為再入點瞬時軌道傾角，φE為再入點地心緯度。

可見，在忽略軌道面外運動(即iE為固定值)的情況下，再入點地心緯度與再入點緯度幅角是一一對應的關系，對再入點緯度的需求可轉(zhuǎn)化為對再入點緯度幅角的需求。

再入大氣飛行階段一般以確定的海拔高度作為轉(zhuǎn)段點，地心距rE為確定值，可得再入點緯度幅角表達式：

(2)

再入點速度與速度傾角一般定義為地速及地速的傾角，而空間軌道是反映飛行器慣性系下的特性。由于地速與慣性速度存在確定的轉(zhuǎn)換關系，因此，本文將速度及傾角的需求轉(zhuǎn)化為慣性系下速度及傾角需求。軌道系下速度與航跡傾角分別為：

(3)

綜上可以看出，忽略軌道面外運動，若飛行器所在的軌道滿足再入點對aE，eE，ωE的需求，即可滿足再入點的約束條件。由此離軌制動問題即可等效為軌道面內(nèi)半長軸、偏心率及近地點幅角調(diào)整的變軌問題。

根據(jù)橢圓軌道運動理論，再入點與制動點滿足能量守恒和動量矩守恒定理，可求得制動點的目標速度和目標傾角，給出脈沖假設下離軌制動的脈沖解,但并非任意情況下均有解。文獻[7]中給出了一個確定地心距上界的公式，并提出當?shù)匦木啻笥谂R界地心距時，需要兩次以上的沖量才能滿足再入接口條件。從軌道控制的角度，可以給出幾何意義直觀的描述方法，即離軌待命軌道與離軌過渡軌道必須存在交點，才存在滿足所有再入點約束的單脈沖解。當離軌待命軌道為理想圓軌時，可以簡化為待命圓軌的地心距不能大于離軌過渡軌道的遠地點，即：

(4)

但以上條件均是基于理想二體假設和脈沖假設而得到。實際工程中制動是一個有限推力持續(xù)作用的過程，不易求得嚴格的解析解。此外由于地球非球形引力使得軌道的瞬時根數(shù)在一個軌道周期內(nèi)周期性波動，對于近地圓軌道衛(wèi)星，半長軸波動可達到10 km以上，文獻[13]中明確提出制動問題不能忽略J2項的影響。因此，本文在考慮地球非球形引力以及有限推力的前提下通過在線數(shù)值迭代求解制動參數(shù)。

2 制動參數(shù)在線確定方法

固定推力制動需要確定三個制動參數(shù):制動角θ、制動速度增量ΔVset以及制動起始點的緯度幅角uopen。其中θ為推力矢量與飛行器速度矢量負向的夾角?？梢钥闯?，制動問題是一個多參數(shù)尋優(yōu)問題。在地面規(guī)劃中，多參數(shù)尋優(yōu)的算法眾多，本文不再一一贅述。但對于在線求解快速決策而言，解的收斂性、算法的可靠性以及算法復雜度的問題，嚴重制約了傳統(tǒng)多參數(shù)尋優(yōu)算法的應用。因此，本文將制動參數(shù)求解問題分解為多次單參數(shù)搜索問題，并且從工程實用的角度，在再入點速度不嚴格限制的情況下，將算法簡化為單參數(shù)搜索問題。

2.1 基于平根數(shù)的制動參數(shù)初值求解

大量文獻研究中以再入點速度、傾角以及位置為約束條件，基于能量守恒和動量矩守恒，求解了制動增量的解析解。本文基于軌道根數(shù)進行研究，因此首先給出一個基于平均軌道根數(shù)的制動初值求解算法。上文已經(jīng)說明，對于制動問題，可轉(zhuǎn)化為軌道面內(nèi)半長軸、偏心率及近地點幅角調(diào)整的變軌問題。令待命軌道的平均根數(shù)為：a0,e0,ω0；過渡軌道的平均根數(shù)為：af,ef,ωf。本文對平均軌道根數(shù)的求解方式采用文獻[9]的方法，在本文2.3節(jié)中給出了短周期項的一階近似公式。

單脈沖有解的條件是兩條軌道存在交點。交點處，兩條共面軌道的緯度幅角u和地心距r一致，據(jù)此列寫方程：

(5)

展開式(5),并令：

efx=efcos(ωf)，efy=efsin(ωf)，e0x=e0cos(ω0)，e0y=e0sin(ω0)。

式(5)整理為：

k1cos(uopen)+k2sin(uopen)=k3

(6)

式中:

根據(jù)式(6)可求得制動點緯度幅角。在式(6)有解的條件下，可以解出平均軌道情況下兩條軌道理論的交點在軌道面內(nèi)所在的位置uopen以及該位置處對應的地心距r。式(6)一般存在兩個解，即過渡軌道與待命軌道存在兩個交點?？紤]到本公式只是提供一個可行的初值解，可選擇離再入點近的交點作為理想制動點初值。

變軌前徑向速度vor及法向速度vof分別為：

(7)

脈沖假設下，變軌后的徑向速度vfr及法向速度vff分別為：

(8)

則制動脈沖大小為：

(9)

制動角大小為：

(10)

綜上即可得到脈沖假設下制動角θ、制動點緯度幅角uopen、制動增量大小ΔVset。在制動弧段較小或精度要求低的軌道控制任務中，可直接利用脈沖假設的制動參數(shù)，將uopen設置在制動弧段的中心點。本節(jié)的計算僅提供制動參數(shù)的初值。

2.2 數(shù)值遞推方程

確定飛行器的位置、速度初值和制動參數(shù)初值后，根據(jù)如下的微分方程進行制動軌道的數(shù)值遞推：

(11)

為簡化參數(shù)搜索的復雜度，本文在數(shù)值遞推軌道過程中根據(jù)軌道根數(shù)實時預報再入點狀態(tài)，通過關機時機確保每條彈道均滿足再入角需求，從而避免對制動參數(shù)ΔVset的搜索。在Kepler軌道假設下，橢圓軌道未來任何時刻的位置、速度都可以通過簡單的幾何轉(zhuǎn)換獲得。但J2項的影響使得理想橢球近似下的軌道預報精度受到影響。文獻[12]中給出了考慮J2項與理想Kepler軌道在再入點航跡傾角上的差異，可以看出不考慮J2項將嚴重影響再入點傾角精度。以下給出在考慮J2項影響下的再入點預報算法。

2.3 基于瞬平根轉(zhuǎn)換的再入點狀態(tài)預報

考慮J2項攝動的情況下，難以得到精確解析解，但在軌道控制理論中，軌道根數(shù)被分解為平均軌道根數(shù)和攝動項。攝動項分解為長周期項和短周期項，軌道根數(shù)預報公式可以按式(12)給出：

(12)

根據(jù)瞬平根轉(zhuǎn)換的基本公式，近似解析的t時刻的瞬時軌道根數(shù)預報公式為：

(13)

式中：σ(t)為t時刻的瞬時軌道根數(shù)；σS(t)為t時刻的短周期攝動項。

本文主要參考文獻[9]給出的一種無奇異點的平根數(shù)迭代法。采用的根數(shù)為σ=[a,i,Ω,ξ,η,λ]T。其中：a為軌道半長軸、i為軌道傾角、Ω為升交點赤經(jīng)、ξ=ecosω，η=-esinω為軌道偏心率e在軌道面內(nèi)的兩個分量，λ=ω+M為近地點幅角ω和平近點角M之和。

表達式中采用了時刻t作為自變量，但對于短周期項，計算中使用的是緯度幅角，與絕對時間并無直接關系。在任一時刻，已知飛行器瞬時軌道根數(shù)時，可以通過有限次迭代計算出其短周期項，從而估計其平均軌道根數(shù)。文獻[9]中給出了考慮偏心率高階小量的短周期項計算公式，復雜度較高。對于近地軌道返回式飛行器，其離軌過渡軌道的偏心率一般較小，因此從工程實際的角度，本文給出僅考慮偏心率一階小量的短周期項計算公式,見式(14)～式(19)。

(14)

(15)

(16)

(17)

(18)

(19)

式中：

(20)

軌道根數(shù)的長周期攝動，是隨時間積分的量，平根數(shù)的一階長周期項為[9]：

(21)

(22)

長周期項需要知道待飛時間，數(shù)值遞推過程按下式估算待飛時間，從而補償軌道根數(shù)的長周期攝動項。

(23)

根據(jù)以上平根數(shù)及長周期項，可以計算出再入點的平均軌道根數(shù)：

(24)

將式(24)得到的平均軌道根數(shù)，代入式(14)～式(19)可得預期再入點位置所對應的軌道根數(shù)短周期項σS(t)=[aS(t),iS(t),ΩS(t),ξS(t),ηS(t),λS(t)]T;進而得到預期再入點位置所對應的軌道瞬時根數(shù)σ(t)=[a(t),i(t),Ω(t),ξ(t),η(t),λ(t)]T。

根據(jù)無奇異點的瞬時根數(shù)，解算出預報的偏心率和近地點幅角:

(25)

2.4 參數(shù)迭代方案

數(shù)值遞推過程中，制動開始后，首先判斷當前軌道的近地點高度rp是否低于理想的再入點地心距rE。如果rp≤rE，則開始按照上一節(jié)給出的公式預報再入點的瞬時軌道根數(shù)。本文以再入角滿足約束為停止制動的條件。根據(jù)軌道瞬根預報值，按下式估算再入點的真近點角：

(26)

圖1 uopen迭代流程圖

將真近點角f(t)=ω(t)+u(t)，代入式(3)得到再入點速度傾角γ(t)。當再入角滿足γ(t)≤γE時，結束制動，并記錄下累計制動速度增量ΔVset。

制動結束后，按慣性滑行繼續(xù)數(shù)值遞推動力學方程，至地心距滿足|r|≤rE時，結束本次遞推，并記錄下遞推結束時的緯度幅角u(t)及速度大小。

遞推過程中以再入角約束確定制動增量，因此再入角滿足約束，而再入點緯度幅角則與目標值存在偏差。因此一次遞推結束后，根據(jù)再入點緯度幅角u(t)與理想再入點緯度幅角uE的偏差，按下式修正制動起始緯度幅角:

uopen←uopen-ku(u(t)-uE)

(27)

式中:ku為修正系數(shù)，取值范圍為0.1～0.9。

由此可以看出，一次軌跡遞推過程中，ΔVset是根據(jù)再入角滿足條件確定，無需迭代求解。需要搜索的參數(shù)只有制動起始點緯度幅角uopen。綜上，一次制動開始點的緯度幅角搜索流程如圖1所示。

參考航天飛機返回任務，給出一個具體算例。假定飛行器位于345 km×345 km的圓軌道，離軌質(zhì)量約95 t，配備兩臺6000 N推力器，設置目標再入點為海拔120 km,慣性速度7882.8 m/s，慣性傾角-1.1°緯度幅角-14.3°。參數(shù)迭代結果如圖2～圖4。

圖2 制動開機點緯度幅角迭代曲線

圖3 制動速度增量迭代曲線

圖4 再入點緯度誤差迭代曲線

仿真算例可以看出，5次迭代后再入點緯度幅角誤差收斂至0.01度范圍內(nèi)，對應的縱向位置誤差在1 km以內(nèi)。由于迭代過程是以軌道根數(shù)預估的慣性傾角作為制動關機點，數(shù)值遞推得到的再入點慣性傾角誤差僅來源于預報誤差。

圖5 再入點慣性速度誤差迭代曲線

圖6 制動角與再入速度誤差的關系

圖7 制動角與制動速度增量的關系

2.5 落點速度與制動角的影響

制動問題概括為通過調(diào)整ΔVset，uopen以及θ使得再入點位置、速度以及航跡傾角滿足uE，VE，rE的約束條件。上一節(jié)的參數(shù)迭代方案，通過滿足uE,rE約束條件得到了ΔVset,uopen，但并沒有保證再入速度VE滿足約束。圖5即為上一節(jié)迭代過程中再入點慣性速度誤差變化曲線，仿真算例顯示再入點速度誤差未收斂至0附近。

由于上一節(jié)制動角θ采用2.1節(jié)中脈沖假設下解析計算的初值，在迭代過程沒有調(diào)整θ，本節(jié)分析調(diào)整制動角θ對速度約束的影響。

設置不同的制動角，重復采用2.4節(jié)的參數(shù)迭代方案，每次迭代收斂后記錄下速度誤差與制動角，以及此制動角對應的制動速度增量。圖6為選擇不同制動角重復采用參數(shù)迭代后對應的再入點速度誤差，圖7為選擇不同制動角所需的制動速度增量。

從圖7可以看出，對于2.4節(jié)的算例，當θ=±27°時，經(jīng)過迭代求解制動速度和制動起點，可以滿足再入點速度誤差到0附近。但此時對應的制動速度增量約為ΔVset≈97 m/s。對比圖6，在θ=0°時制動速度增量ΔVset≈88 m/s，比θ=±27°可省約10%的燃料，而再入速度誤差僅增加1.5 m/s。

可見為了保證所有再入點狀態(tài)約束，可能需要付出極大的能源代價。實際工程應用中，對于再入飛行器，特別是大升力體再入飛行器，對再入點初始速度偏差具備較強的適應能力。航天飛機的離軌制動任務中，再入點約束條件是使再入點徑向與切向速度滿足一定的線性關系[20]。對于再入點速度的限制，一般可以在任務規(guī)劃階段通過設計合理的待命軌道來保證再入點速度在一定范圍之內(nèi)。因此，本文不考慮對制動角進行搜索，θ按2.1節(jié)給出的計算公式求解。

3 制動關機策略

理想的制動軌跡應該是緯度幅角及航跡傾角同時得到滿足。但實際飛行中，由于推力誤差、比沖誤差，飛行器質(zhì)量誤差，姿態(tài)控制誤差等因素，會使得真實的軌跡與預測軌跡存在一定的差異。實際軌跡很難保證再入點落點、速度、傾角同時滿足需求。因此，需要考慮以何種指標決定制動關機時機。

對于軌道再入飛行器而言，最關鍵的再入點狀態(tài)量為航跡傾角和落點位置。本文2.3節(jié)的公式給出了基于預報再入點瞬時根數(shù)求解再入點速度傾角的算法。制動過程，同樣可以根據(jù)軌道預報的瞬時根數(shù)a(t),e(t),ω(t)，代入式(2)得到再入點對應的緯度幅角u(t)?？梢姡苿娱_機過程隨時可以根據(jù)當前的飛行狀態(tài)，計算出此時結束制動所對應的再入點緯度幅角及再入點速度傾角。

若不考慮制動過程中調(diào)整制動角，則制動一旦開始，控制系統(tǒng)只有惟一的控制量——制動關機時間。單一的控制量不能保證再入點緯度幅角和速度傾角同時滿足約束。如果以再入點緯度幅角滿足約束為條件結束制動，則再入點傾角會存在一定范圍的散布；反之，如果以再入點傾角滿足約束為條件結束制動，則再入點緯度幅角會存在一定范圍的散布。飛行任務中，以何種指標決定制動關機時間，主要取決于該指標對偏差的敏感度。下一章中，針對2.4節(jié)提出的參考航天飛機再入返回軌道的算例，通過數(shù)值仿真對比兩種不同的制動關機策略對應的再入點散布情況。

4 仿真分析

針對2.4節(jié)提出的算例，考慮如下的狀態(tài)散布：

1)初始軌道半長軸：±3 km隨機誤差；

2)初始軌道平均偏心率：0.004隨機誤差；

3)制動過程姿態(tài)：±1°隨機擾動；

4)飛行器質(zhì)量：1%隨機誤差(-950～950 kg)；

5)制動推力：3%隨機誤差(-360～360 N)。

制動參數(shù)在線確定，閉環(huán)制動仿真過程分別以預報的再入點緯度幅角滿足約束關機(策略1)和以預報的再入點航跡傾角滿足約束關機(策略2)。分別進行200次蒙特卡洛打靶仿真，仿真結果如圖8～圖17。其中圖8～圖10為20條按策略1仿真的制動過程速度、高度和航跡傾角變化曲線，圖11～圖13為20條按策略2仿真的制動過程速度、高度和航跡傾角變化曲線。圖14～圖17為兩種策略200組蒙特卡洛仿真的再入點緯度幅角和慣性傾角散布。不同的關機策略產(chǎn)生的誤差散布范圍如表1所示：

表1 不同關機策略的誤差散布Table 1 Errors dispersion with different shutdown strategies

圖8 策略1慣性速度變化曲線

圖9 策略1高度變化曲線

圖10 策略1慣性傾角變化曲線

圖11 策略2慣性速度變化曲線

圖12 策略2高度變化曲線

圖13 策略2慣性傾角變化曲線

圖14 策略1再入點緯度幅角偏差散布

圖15 策略2再入點緯度幅角偏差散布

圖16 策略1再入點慣性傾角偏差散布

圖17 策略2再入點慣性傾角偏差散布

從仿真結果可以看出，策略1以預報緯度幅角滿足約束條件為制動關機條件時，再入點縱向位置散布可以達到極高的精度(小于3.5 km)，誤差源包含預報誤差及主引擎關機后滑行過程中姿控推力器產(chǎn)生的軌跡擾動。策略2以預報再入點傾角滿足約束為制動關機條件，再入點傾角誤差散布極小(≤0.0015°)，但緯度幅角的散布則可達到120 km。

可見在制動角不控的工況下，通過實時預報再入點狀態(tài)確定制動關機時機，至少可以保證再入點緯度幅角或再入點傾角中的一個量具有較高的精度。工程應用中，可根據(jù)實際任務對再入角和緯度幅角誤差的容忍程度來確定制動關機的時機。

5 結 論

本文研究了一種在考慮地球非球形引力J2項影響下，基于瞬平根轉(zhuǎn)換關系對再入點狀態(tài)進行預報的自主離軌制動方法。算法首先根據(jù)再入點緯度和傾角的約束，通過在線數(shù)值迭代搜索制動參數(shù)，確定制動角、制動開機時間和制動速度增量大小。制動過程中通過近似解析的預報算法，計算再入點傾角和緯度幅角實時確定制動關機時刻。

在考慮離軌待命圓軌半長軸、偏心率存在初始散布，飛行器質(zhì)量、推力、姿態(tài)均存在偏差的情況下，按照本文設計的離軌制動方案，以保證再入點航跡傾角約束和保證再入點緯度幅角約束的兩種制動策略進行蒙特卡洛打靶仿真，仿真結果顯示出再入點飛行狀態(tài)精度較高?？梢詾榇髿鈱觾?nèi)再入飛行創(chuàng)造良好的初始條件。

離軌制動工程實踐中還存在碰撞、敏感區(qū)域過境風險等與飛行任務頂層規(guī)劃密切相關的問題，難以完全依靠在軌自主決策。實際應用中可通過實時下傳自主規(guī)劃的制動參數(shù)，通過地面設備對軌跡進行復核復算。當飛行軌跡確實存在風險時，可通過注入新的目標再入點狀態(tài)或直接注入制動角調(diào)整軌跡，從而增強本方法在工程實踐中的靈活性。