黃漢斌，梁祿揚(yáng)，楊業(yè)

北京航天自動(dòng)控制研究所,北京 100854

以航天飛機(jī)[1-2]為代表的升力式飛行器具有巨大優(yōu)勢，高升阻比的氣動(dòng)布局，使其可以充分利用氣動(dòng)力實(shí)現(xiàn)大航程和高精度飛行任務(wù)。再入段[3]通常是約束最復(fù)雜、耗時(shí)最長、飛行距離最遠(yuǎn)的飛行段，再入飛行軌跡規(guī)劃與制導(dǎo)算法是實(shí)現(xiàn)再入精確飛行的核心技術(shù)[3-4]。因此，再入軌跡規(guī)劃與制導(dǎo)算法受到廣泛關(guān)注。

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

獲得標(biāo)準(zhǔn)軌跡之后，Dukeman[16]以能量作為自變量，建立小擾動(dòng)線性化模型，通過反饋航程、高度和航跡傾角的偏差設(shè)計(jì)LQR控制器對(duì)阻力加速度進(jìn)行跟蹤；Hanson[17]、Leavitt[18]等利用阻力加速度的動(dòng)態(tài)誤差特性，根據(jù)標(biāo)準(zhǔn)軌跡和實(shí)際阻力加速度求出制導(dǎo)律，實(shí)現(xiàn)縱向制導(dǎo)。橫向制導(dǎo)，通過改變傾側(cè)角符號(hào)[16-18]將航向角誤差限制在走廊內(nèi)。

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

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

1 運(yùn)動(dòng)模型

1.1 有量綱運(yùn)動(dòng)模型

假設(shè)地球?yàn)閯蛸|(zhì)圓球，忽略地球自轉(zhuǎn)，微分方程的6個(gè)變量為高度H，經(jīng)度λ，緯度φ，速度Vd，航跡傾角γ，航向角ψ，建立式(1)中以時(shí)間t為自變量的運(yùn)動(dòng)模型。其中，re為地球半徑，g為重力加速度，控制量為傾側(cè)角σ和攻角α，升力加速度和阻力加速度采用式(2)和式(3)計(jì)算：

(1)

(2)

(3)

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

1.2 阻力加速度倒數(shù)約束

(4)

(5)

(6)

(7)

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

利用同樣的方法對(duì)倒數(shù)第2個(gè)能量點(diǎn)進(jìn)行處理，得到圖1中由N個(gè)節(jié)點(diǎn)(紅色點(diǎn))組成的阻力加速度能量的上下界。

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

(8)

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

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

2 再入制導(dǎo)算法

基于阻力加速度倒數(shù)剖面的標(biāo)準(zhǔn)軌跡跟蹤制導(dǎo)算法分為2個(gè)步驟：① 生成阻力加速度標(biāo)準(zhǔn)軌跡；② 計(jì)算跟蹤制導(dǎo)律。其中在線生成并更新滿足各種約束的阻力加速度剖面是本文算法的核心。本文算法結(jié)構(gòu)框圖如圖3所示。

圖3 本文算法結(jié)構(gòu)框圖Fig.3 Structure diagram of proposed algorithm

2.1 在線軌跡生成與更新

2.1.1 待飛航程預(yù)測

標(biāo)準(zhǔn)軌跡制導(dǎo)算法的核心之一在于軌跡的生成，軌跡生成方法主要依據(jù)一般是待飛航程的預(yù)測，待飛航程定義為當(dāng)前點(diǎn)地心矢徑到目標(biāo)點(diǎn)地心矢徑在地球表面所構(gòu)成的大圓弧的長度，計(jì)算公式為

(9)

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

(10)

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

Smax=

(11)

Smin=

(12)

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

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

待飛航程與其上下界對(duì)應(yīng)的比例系數(shù)為

(13)

因此，利用式(14)得到近似滿足待飛航程的阻力加速度倒數(shù)剖面上的N個(gè)節(jié)點(diǎn)。

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

i=0,1,…,N-1

(14)

i=1,2,…,N-1

(15)

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

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

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

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

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

i=1,2,…,N-1

(16)

2.1.3 剖面更新算法

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

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

(17)

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

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

(18)

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

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

2.2 阻力加速度跟蹤制導(dǎo)

2.2.1 縱向制導(dǎo)

縱向制導(dǎo)對(duì)規(guī)劃的阻力加速度剖面進(jìn)行跟蹤，再入過程中通常采用與速度相關(guān)的固定的攻角剖面，僅通過傾側(cè)角進(jìn)行跟蹤制導(dǎo)，利用阻力加速度動(dòng)態(tài)誤差特性推導(dǎo)制導(dǎo)律[14]。動(dòng)態(tài)誤差特性公式為

(19)

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

|σ|=

(20)

(21)

2.2.2 橫向制導(dǎo)

側(cè)向制導(dǎo)旨在消除側(cè)向誤差，不去跟蹤某一標(biāo)準(zhǔn)軌跡，因此，僅需要將航向角偏差的限定在以速度為函數(shù)的范圍內(nèi)，即當(dāng)航向角偏差超出邊界，改變傾側(cè)角符號(hào)來實(shí)現(xiàn)橫向制導(dǎo)。

航向角偏差的計(jì)算公式為

Δψ=ψ-az

(22)

式中：ψ為當(dāng)前航向角；az為目標(biāo)點(diǎn)關(guān)于當(dāng)前點(diǎn)的方位角，計(jì)算公式如式(23)所示:

(23)

式中：λf為終端時(shí)刻經(jīng)度；φf為終端時(shí)刻緯度。

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

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

3 仿真結(jié)果

3.1 仿真模型及環(huán)境

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

仿真環(huán)境：Intel Core i7-4790、Windows XP系統(tǒng)、MATLAB R2013b。

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

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

Vd/(km·s-1)Δψu(yù)p/(°)87775144162140.64

3.2 標(biāo)準(zhǔn)任務(wù)仿真

利用表2中給出的標(biāo)準(zhǔn)任務(wù)下飛行器再入初始和終端要求，根據(jù)表中條件可以算出初始和終端能量與阻力加速度。

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

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

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

(24)

表2 標(biāo)準(zhǔn)任務(wù)初始和終端要求Table 2 Standard task initial and end requirements

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

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

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

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

圖10 傾側(cè)角Fig.10 Angle of bank

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

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

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

圖14給出了經(jīng)緯度結(jié)果，飛行結(jié)束點(diǎn)準(zhǔn)確的落在再入段與終端能量管理段交班的圓上，誤差為1.076 km，算法精度較高。

圖12 高度-航程關(guān)系Fig.12 Altitude vs range

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

圖14 經(jīng)緯度Fig.14 Longitude and latitude

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

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

3.3 隨機(jī)誤差適應(yīng)性仿真

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

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

表4 蒙特卡羅打靶偏差項(xiàng)Table 4 Deviation term of Monte Carlo simulation

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

3.4 飛行任務(wù)適應(yīng)性仿真

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

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

仿真結(jié)果表明，本文提出的方法靈活度高，相對(duì)于前人研究，適應(yīng)不同的飛行任務(wù)需求，可以充分利用飛行器在走廊內(nèi)的飛行能力，完成相應(yīng)的飛行任務(wù)。

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

圖16 不同任務(wù)下的經(jīng)緯度Fig.16 Longitude and latitude corresponding to different tasks

4 結(jié) 論

本文提出了一種基于阻力加速度倒數(shù)走廊的再入制導(dǎo)與跟蹤方法?；谧枇铀俣鹊箶?shù)走廊求出滿足航程要求和約束條件的解析阻力加速度倒數(shù)剖面，然后將其轉(zhuǎn)化成阻力加速度剖面，分別設(shè)計(jì)縱向和橫向制導(dǎo)律進(jìn)行跟蹤制導(dǎo)。仿真結(jié)果表明：

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

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

3) 本文的算法精度較高，有較強(qiáng)的隨機(jī)誤差適應(yīng)性和飛行任務(wù)適應(yīng)性，能夠?qū)崿F(xiàn)工程在線應(yīng)用。