基于阻力加速度倒數剖面的再入軌跡規劃與制導方法

黃漢斌，梁祿揚，楊業

北京航天自動控制研究所,北京 100854

以航天飛機[1-2]為代表的升力式飛行器具有巨大優勢，高升阻比的氣動布局，使其可以充分利用氣動力實現大航程和高精度飛行任務。再入段[3]通常是約束最復雜、耗時最長、飛行距離最遠的飛行段，再入飛行軌跡規劃與制導算法是實現再入精確飛行的核心技術[3-4]。因此，再入軌跡規劃與制導算法受到廣泛關注。

升力式飛行器再入軌跡規劃與制導算法[4]，分為標準軌跡跟蹤制導和數值預測校正制導兩種，其中數值預測校正制導算法需要迭代數值積分，由于計算效率不高，不能在線應用，而傳統的標準軌跡制導算法基于離線規劃的軌跡，不夠靈活，因此需要研究在線軌跡規劃與制導算法。由于再入過程中阻力加速度可以準確測量，同時可以克服模型等多種誤差帶來的影響，因此采用阻力加速度剖面作為標準軌跡[3]具有一定的優勢。在生成阻力加速度剖面時，需要滿足多種約束，因此需要建立阻力加速度約束走廊，同時由于再入終端時間不確定，采用能量作為阻力加速度自變量[5]，描述阻力加速度走廊，將高度和速度信息包含在能量中。為了獲得滿足待飛航程要求的阻力加速度-能量剖面，文獻[5-11]采用多段線性函數將阻力加速度剖面設定為參數化的形式，進而利用優化的方法對參數進行求解，以得到標準軌跡；但線性函數轉折點不光滑，不利于跟蹤制導，文獻[12] 則重新定義能量高度，采用分段三次函數描述剖面，利用簡單實時積分的方法迭代獲得光滑標準軌跡。阻力加速度標準軌跡必須滿足阻力加速度走廊約束，文獻[13-14]直接將阻力加速度剖面設定為走廊上下界的函數來進行優化，以使得阻力加速度剖面始終在走廊內。為了加快優化速度，簡化優化過程，Mease等[15]采用降階的動力學模型，在進行縱向阻力加速度剖面設計時，采用三段線性阻力加速度剖面進行優化，該方法在一定程度上加快了優化速度。

獲得標準軌跡之后，Dukeman[16]以能量作為自變量，建立小擾動線性化模型，通過反饋航程、高度和航跡傾角的偏差設計LQR控制器對阻力加速度進行跟蹤；Hanson[17]、Leavitt[18]等利用阻力加速度的動態誤差特性，根據標準軌跡和實際阻力加速度求出制導律，實現縱向制導。橫向制導，通過改變傾側角符號[16-18]將航向角誤差限制在走廊內。

標準軌跡生成時，優化的過程通常利用數值積分計算待飛航程，且初值隨機選取，因而計算時間較長，Yang等[19]將阻力加速度-能量剖面分成多段，分段采用三次樣條函數擬合，然后利用卡爾丹分解，推導出待飛航程的解析計算公式，大大提高了計算效率，但是卡爾丹分解推導過程復雜，本文對該解析航程預測方法進行了改進，消除卡爾丹分解過程，同時，基于該方法，提出了一種滿足終端和過程約束的阻力加速度軌跡在線生成與更新的新方法，并將其與跟蹤制導相結合，完成了再入軌跡在線規劃和跟蹤制導的閉環仿真。

本文首先利用過程約束、平衡滑翔約束建立阻力加速度倒數-能量走廊，將飛行器再入終端約束加入走廊中，選取阻力加速度倒數上下界的多個點，用三次樣條函數對其進行擬合。然后，利用簡化的航程計算公式，得到待飛航程的上下界，當待飛航程在約束走廊對應的航程上下界以內時，滿足待飛航程和約束的阻力加速度倒數剖面一定在走廊上下界之間，因此，阻力加速度剖面可以利用待飛航程以及待飛航程的上下界的近似線性關系生成阻力加速度剖面，并利用二分法對剖面進一步優化，從而得到阻力加速度剖面，之后根據待飛航程偏差，利用當前剖面來更新阻力加速度剖面。最后，利用阻力加速度的動態誤差特性建模，對其進行跟蹤。橫側向制導根據航向角偏差，改變傾側角符號，來消除橫向誤差。

1 運動模型

1.1 有量綱運動模型

假設地球為勻質圓球，忽略地球自轉，微分方程的6個變量為高度H，經度λ，緯度φ，速度Vd，航跡傾角γ，航向角ψ，建立式(1)中以時間t為自變量的運動模型。其中，re為地球半徑，g為重力加速度，控制量為傾側角σ和攻角α，升力加速度和阻力加速度采用式(2)和式(3)計算：

(1)

(2)

(3)

式中：m、Sref分別為飛行器質量和參考面積；ρ為大氣密度；CL、CD分別為升力系數和阻力系數。

1.2 阻力加速度倒數約束

(4)

(5)

(6)

(7)

由于飛行器初始速度、高度，以及終端速度、高度給定，則初始能量、阻力加速度，以及終端能量、阻力加速度即給定，阻力加速度走廊需要滿足初始和終端約束，因此需要對阻力加速度走廊進行處理。

利用同樣的方法對倒數第2個能量點進行處理，得到圖1中由N個節點(紅色點)組成的阻力加速度能量的上下界。

為了簡化計算待飛航程，將N個節點構成的無量綱阻力加速度走廊轉換成圖2中N個無量綱阻力加速度倒數節點三次樣條插值后的走廊。得到式(8)中三次樣條函數描述的無量綱阻力加速度倒數約束走廊：

(8)

圖1 無量綱阻力加速度走廊Fig.1 Dimensionless drag acceleration corridor

圖2 無量綱阻力加速度倒數走廊Fig.2 Inverse dimensionless drag acceleration corridor

2 再入制導算法

基于阻力加速度倒數剖面的標準軌跡跟蹤制導算法分為2個步驟：① 生成阻力加速度標準軌跡；② 計算跟蹤制導律。其中在線生成并更新滿足各種約束的阻力加速度剖面是本文算法的核心。本文算法結構框圖如圖3所示。

圖3 本文算法結構框圖Fig.3 Structure diagram of proposed algorithm

2.1 在線軌跡生成與更新

2.1.1 待飛航程預測

標準軌跡制導算法的核心之一在于軌跡的生成，軌跡生成方法主要依據一般是待飛航程的預測，待飛航程定義為當前點地心矢徑到目標點地心矢徑在地球表面所構成的大圓弧的長度，計算公式為

(9)

在平衡滑翔條件下航跡傾角γ較小，cosγ≈1，同時由于航向角偏差Δψ較小，cos Δψ≈1，采用能量作為待飛航程微分公式的自變量，則可以推導出待飛航程計算公式為

(10)

則將式(8)代入式(10)，可以得到待飛航程上下界的計算公式為

Smax=

(11)

Smin=

(12)

2.1.2 阻力加速度剖面在線解析生成算法

利用阻力加速度倒數剖面，通過式(10)計算待飛航程，則滿足待飛航程要求的阻力加速度倒數剖面與阻力加速度倒數上下界近似成線性關系，則可以利用這種關系，首先得到一條相近的阻力加速度剖面，然后利用二分法得到精確的阻力加速度剖面。

待飛航程與其上下界對應的比例系數為

(13)

因此，利用式(14)得到近似滿足待飛航程的阻力加速度倒數剖面上的N個節點。

然后式(15)利用三次樣條函數插值，得到倒數剖面。

i=0,1,…,N-1

(14)

i=1,2,…,N-1

(15)

對與上下邊界有比例系數η的節點進行插值后，節點之間的點不一定成比例系數，利用式(11)或者式(12)估計的結果不一定等于待飛航程，因此，需要進一步迭代求取剖面，這里采用二分法進行迭代，迭代過程如圖4所示。

經過迭代，得到滿足待飛航程要求的阻力加速度倒數剖面，該剖面與阻力加速度倒數上下界仍然成線性關系，該比例系數仍然采用η描述。得到圖5中的阻力加速度倒數剖面。

同時，阻力加速度倒數剖面可以無損轉化為阻力加速度剖面，式(16)給出轉化方式，由此便得到了滿足航程要求的阻力加速度剖面。

圖4 倒數剖面迭代過程Fig.4 Reciprocal profile iteration process

圖5 無量綱阻力加速度倒數剖面Fig.5 Inverse dimensionless drag acceleration profile

i=1,2,…,N-1

(16)

2.1.3 剖面更新算法

阻力加速度剖面依據當前點到目標點所構成的縱向平面進行規劃，由于控制量傾側角以及模型誤差和氣動誤差的存在，飛行器存在橫向機動來消耗能量使得飛行器跟蹤阻力加速度剖面，所以飛行器最終不能夠到達終點(圖6)，因此需要在一定的周期下對阻力加速度剖面進行更新。

由式(10)～式(12)可知，待飛航程在阻力加速度倒數走廊內近似成線性分布，因此可以根據待飛航程的偏差更新η，從而得到新的阻力加速度剖面。由當前飛行器位置到目標之間的待飛航程S(e)，同時，根據阻力加速度-能量剖面以及待飛航程計算公式，可以得到當前能量點處，沿規劃的阻力加速度剖面的剩余待飛航程Sego(e)，由于線性關系的存在，令

(17)

則滿足剩余待飛航程的阻力加速度倒數剖面與上下界之間的比例系數為

η(j+1)=η(j)+Δη

(18)

由式(14)、式(15)和式(18)，可得到滿足待飛航程和約束的新的阻力加速度剖面。

圖6 阻力加速度剖面Fig.6 Drag acceleration profile

2.2 阻力加速度跟蹤制導

2.2.1 縱向制導

縱向制導對規劃的阻力加速度剖面進行跟蹤，再入過程中通常采用與速度相關的固定的攻角剖面，僅通過傾側角進行跟蹤制導，利用阻力加速度動態誤差特性推導制導律[14]。動態誤差特性公式為

(19)

式中：ξ為阻尼系數；ωn為自然頻率；k1為積分系數，實際阻力加速度一階導D′和D″二階導采用文獻[18]得到的推導公式，代入式(19)，得到縱向制導律的表達式為

|σ|=

(20)

(21)

2.2.2 橫向制導

側向制導旨在消除側向誤差，不去跟蹤某一標準軌跡，因此，僅需要將航向角偏差的限定在以速度為函數的范圍內，即當航向角偏差超出邊界，改變傾側角符號來實現橫向制導。

航向角偏差的計算公式為

Δψ=ψ-az

(22)

式中：ψ為當前航向角；az為目標點關于當前點的方位角，計算公式如式(23)所示:

(23)

式中：λf為終端時刻經度；φf為終端時刻緯度。

航向角偏差的上下界±Δψup，采用插值函數給出，如圖7和表1所示，當航向角偏差超出邊界，傾側角反向，使得航向角偏差減小。應當注意的是，由于傾側角反向時，有速度限制，并不能立即使得航向角偏差減小，因此在設計航向角偏差邊界時，要留有一定的余量。

圖7 航向角偏差走廊Fig.7 Azimuth deviation corridor

3 仿真結果

3.1 仿真模型及環境

對縱向制導律即傾側角的大小，設定上界，|σ|≤90°，傾側角反向速度小于15 (°)/s，攻角剖面采用圖8中與速度相關的固定剖面。

仿真環境：Intel Core i7-4790、Windows XP系統、MATLAB R2013b。

圖8 攻角剖面Fig.8 Angle of attack profile

表1 航向角誤差上下界插值表Table 1 Azimuth deviation corridor interpolation table

Vd/(km·s-1)Δψup/(°)87775144162140.64

3.2 標準任務仿真

利用表2中給出的標準任務下飛行器再入初始和終端要求，根據表中條件可以算出初始和終端能量與阻力加速度。

仿真過程中采用的阻力加速度剖面更新周期為1 s，再入段終點為終端能量管理段的起點，即當距離目標小于終端能量管理段航程STAEM時，進入終端能量管理段。選取STAEM=100 km，式(24)為待飛航程計算公式：

S=re·arccos[sinφsinφf+

cosφcosφfcos(λf-λ)]-STAEM

(24)

表2 標準任務初始和終端要求Table 2 Standard task initial and end requirements

仿真結果顯示，首次軌跡規劃周期小于0.023 s，每更新一次軌跡所需時間小于0.005 s，(因仿真環境不同，存在一定差異)，而以文獻[19]中未改進的待飛航程預測方法為基礎來生成軌跡的時間多次仿真均大于0.09 s，由于本文方法省去了卡爾丹分解過程，減少了待飛航程計算時的矩陣求逆和解方程的次數，軌跡生成時間得到明顯減少，滿足工程在線應用要求。

圖9給出了阻力加速度走廊以及阻力加速度指令和實際阻力加速度的曲線，起始段有一個小的波動，是由于阻力加速度指令變大，而高空大氣密度較小，初始航跡傾角不理想，因此跟蹤阻力加速度指令，需要調整過程；在能量大約為0.92時，實際阻力加速度有一個小的波動，這是由于在此時速度Vd=3 km/s，攻角指令在此處非線性，對于升力系數和阻力系數產生影響而造成的，調整由非線性造成的影響需要一定的時間。

圖10為制導給出的傾側角指令，初始傾側角直接達到最大值，為使得實際阻力加速度盡快調節到與跟蹤指令一致，傾側角直接達到了最大限幅值。傾側角反向，是由于航向角達到了調整的邊界，在橫向制導律的作用下，傾側角改變符號。

圖9 無量綱阻力加速度Fig.9 Dimensionless drag acceleration

圖10 傾側角Fig.10 Angle of bank

圖11給出了路徑約束值，在走廊內的阻力加速度滿足路徑約束的限制。圖12給出了高度-航程的曲線，終端高度Hf=26 975 m，滿足要求。

圖11 熱流、過載和動壓變化曲線Fig.11 Variation curves of heating rate, normal load and dynamic pressure

圖13分別給出了速度、彈道傾角(航跡傾角)、航向角的曲線，終端速度Vf=899.54 m/s滿足要求，彈道傾角小于0°，并且絕對值維持在較小水平，所以高度沒有發生跳躍；為了消除橫向誤差，航向角會在初始航向角附近變化，來減小航向角偏差，仿真終端航向角誤差為0.478°。

圖14給出了經緯度結果，飛行結束點準確的落在再入段與終端能量管理段交班的圓上，誤差為1.076 km，算法精度較高。

圖12 高度-航程關系Fig.12 Altitude vs range

圖13 速度、航道傾角和航向角變化曲線Fig.13 Variation curves of speed, flight path angle and azimuth

圖14 經緯度Fig.14 Longitude and latitude

由于傾側角作為制導律輸出變量，會產生較大的橫向偏差，因此軌跡更新周期要小，同時要考慮到計算效率問題，更新周期不能過小，表3給出了不同周期下的終端位置誤差，可以看到采用1 s作為更新周期，可以兼顧算法精度與計算效率。

表3 不同更新周期下終端位置誤差Table 3 Terminal position error in different update cycles

3.3 隨機誤差適應性仿真

為了驗證本文方法對于不確定性誤差的適應性，采用蒙特卡羅打靶方法進行多次仿真，選取的打靶仿真偏差項如表4所示，仿真過程中10個參數的誤差隨機選取，為了保證阻力加速度剖面調整不超出邊界，選取的待飛航程應該距離阻力加速度上下界所對應的待飛航程有一定的界限，以留出足夠的余量來更新阻力加速度剖面來實現任務。

再入段結束點的位置、速度、高度、航跡傾角、以及航向角偏差對于末端制導具有重要影響，表5給出蒙特卡羅打靶下終端的仿真結果，與終端要求進行對比，可以看出，本文方法能夠應用在有一定偏差的情況下，有較強的適應性。

表4 蒙特卡羅打靶偏差項Table 4 Deviation term of Monte Carlo simulation

表5 蒙特卡羅打靶仿真結果Table 5 Results of Monte Carlo simulation

3.4 飛行任務適應性仿真

本文提出的方法靈活易用，能夠適應不同任務的需求，分別采用初始待飛航程為2 000、3 000、 4 000 km的任務進行仿真，即相對于表2中的標準任務的要求，只改變終端經緯度，使之對應相應的待飛航程，其他仿真條件及算法參數不變。

圖15給出了3種任務對應的飛行過程中阻力加速度剖面，可以看到算法能夠適應3種任務的要求，使飛行器達到相應的阻力加速度；圖16給出3種任務下，飛行器的經緯度曲線，3種任務下算法均能達到較高的精度，使飛行器完成再入任務。

仿真結果表明，本文提出的方法靈活度高，相對于前人研究，適應不同的飛行任務需求，可以充分利用飛行器在走廊內的飛行能力，完成相應的飛行任務。

圖15 不同任務下的無量綱阻力加速度剖面Fig.15 Dimensionless drag acceleration profile corresponding to different tasks

圖16 不同任務下的經緯度Fig.16 Longitude and latitude corresponding to different tasks

4 結 論

本文提出了一種基于阻力加速度倒數走廊的再入制導與跟蹤方法。基于阻力加速度倒數走廊求出滿足航程要求和約束條件的解析阻力加速度倒數剖面，然后將其轉化成阻力加速度剖面，分別設計縱向和橫向制導律進行跟蹤制導。仿真結果表明：

1) 由于待飛航程在阻力加速度倒數走廊內近似成線性關系，因此，本文提出的算法效率較高，可以在線應用。

2) 由于生成和更新的阻力加速度剖面與阻力加速度走廊約束相關，因此，可以充分利用飛行器的飛行能力。

3) 本文的算法精度較高，有較強的隨機誤差適應性和飛行任務適應性，能夠實現工程在線應用。