多禁飛區(qū)高超聲速滑翔飛行器再入機動制導(dǎo)

高 楊,蔡光斌,2,張勝修,徐 慧

(1.火箭軍工程大學(xué) 導(dǎo)彈工程學(xué)院, 西安 710025; 2.西北工業(yè)大學(xué) 航天學(xué)院, 西安 710072)

近年來，高超聲速滑翔飛行器的再入機動制導(dǎo)研究逐漸成為世界各大國航天工業(yè)發(fā)展的一個熱點[1-5]。高超聲速滑翔飛行器再入飛行必須滿足復(fù)雜的約束條件，其中包括傳統(tǒng)的熱流率、動壓、過載約束和終端約束，也包括自然、軍事等因素造成的復(fù)雜禁飛區(qū)約束[4]。飛行器的禁飛區(qū)是一種路徑約束，在飛行器最優(yōu)軌跡的解算中，如果禁飛區(qū)的數(shù)量越多或者模型比較復(fù)雜，那么飛行器在軌跡優(yōu)化的路徑約束就越多，優(yōu)化軌跡的解算難度就越大。現(xiàn)今研究禁飛區(qū)突防的成果較少，且主要集中在離線求取最優(yōu)解，其中，偽譜法是數(shù)值求解方法的代表，采用不同的配點和多項式插值方法達到全局優(yōu)化的目的，這些方法無法應(yīng)用到在線的機動制導(dǎo)中，尤其是在禁飛區(qū)信息較復(fù)雜或者信息缺失的情況下[6]。因而，預(yù)測反饋的制導(dǎo)機制在飛行中作用明顯[7]，一些基于誤差校正的神經(jīng)網(wǎng)絡(luò)與系統(tǒng)辨識方法被廣泛應(yīng)用[8-11]，然而這些研究中對于未知形狀的禁飛區(qū)約束條件研究并不成熟，當(dāng)預(yù)測數(shù)據(jù)因為禁飛區(qū)等因素急劇變化時，會導(dǎo)致反饋誤差變大而導(dǎo)致制導(dǎo)錯誤。部分學(xué)者在研究昆蟲生物學(xué)觸角反饋機制的基礎(chǔ)上，使用較少的“觸角”探測前方未知區(qū)域并進行反饋，從而制定制導(dǎo)方案[12-14]。文獻[12]在移動機器人路徑規(guī)劃研究中，首先提出了基于觸角的方法，其中觸角代表了機器人可能的運動方向。文獻[14]使飛行器連續(xù)產(chǎn)生左右雙觸角對前方進行探測，通過跟蹤標(biāo)準(zhǔn)軌跡并使用側(cè)滑角正負(fù)反轉(zhuǎn)策略規(guī)避不同形狀的禁飛區(qū)。

鑒于現(xiàn)有研究結(jié)果中存在問題，本文改進了軌跡預(yù)測時觸角的設(shè)計與側(cè)滑角的瞬變策略，引入了航向角誤差限制條件和延時計數(shù)開關(guān)，計算了規(guī)避策略的優(yōu)先級，軌跡能夠順利規(guī)避禁飛區(qū)并到達終點。通過魯棒性仿真驗證了該策略預(yù)測過程計算時間足夠小，滿足機動制導(dǎo)策略的要求。

1 再入模型構(gòu)建

高超聲速滑翔飛行器的動力學(xué)模型以及再入飛行過程中的約束條件模型是高超聲速滑翔飛行器機動制導(dǎo)過程的基礎(chǔ)，下面對這些模型與條件分別進行定義。

1.1 高超聲速滑翔飛行器再入運動學(xué)建模

針對高超聲速滑翔飛行器在再入飛行階段的機動制導(dǎo)過程，如圖1所示，在地心赤道旋轉(zhuǎn)坐標(biāo)系內(nèi)，忽略地球自轉(zhuǎn)角的影響，建立高超聲速滑翔飛行器的動力學(xué)方程[20]：

其中，r是高超聲速滑翔飛行器地心距，V是飛行器的地球相對速度，ψ與γ是飛行器的航向角與航跡角，θ和φ分別是飛行器所處的經(jīng)度緯度，m是飛行器的質(zhì)量，g是重力加速度，σ是飛行器的側(cè)滑角。D=ρV2SrefCD/2與L=ρV2SrefCL/2是飛行器在飛行過程中氣動阻力與升力，其中ρ是空氣密度，Sref是飛行器的參考橫截面積，CL與CD分別是與飛行器攻角有關(guān)的空氣動力學(xué)參數(shù)。不同的飛行器由于設(shè)計的不同，導(dǎo)致升力阻力的空氣動力學(xué)參數(shù)不同以及升阻比的不同。

圖1 地心赤道旋轉(zhuǎn)坐標(biāo)系下飛行器動力學(xué)方程建模

1.2 常規(guī)約束條件建模

在臨近空間的長距離、無動力飛行過程中，高超聲速滑翔飛行器會受到多種約束的影響，這些因素都會影響制導(dǎo)的結(jié)果。本文橫向制導(dǎo)策略在預(yù)測的過程中，每條預(yù)測“觸角”上預(yù)測時的每一個積分步長端點為一個飛行器約束檢測點，預(yù)測計算時每到一個檢測點都要檢測是否遵循常規(guī)約束條件：

一是起始狀態(tài)與終止?fàn)顟B(tài)的端點約束，它是對飛行器起始與終止位置、速度和姿態(tài)的約束，也是飛行器不同狀態(tài)轉(zhuǎn)換的交接班條件，起始條件通常為給定的任務(wù)起始狀態(tài)[19]，令X=[r,θ,φ,V,γ,ψ]為飛行器狀態(tài)矩陣，εX為較小的常值矩陣，終端約束模型建立如下：

|X(tf)-Xf|≤εX

二是飛行器的常規(guī)路徑約束，分別是飛行器的熱流率約束、過載約束和動壓約束[10]，建立模型分別為

Q=KQρ0.5V3.15≤Qmax

三是航向角約束。在高超聲速滑翔飛行器因禁飛區(qū)等限制原因機動飛行的過程中，為了保證飛行器轉(zhuǎn)彎規(guī)避禁飛區(qū)后仍能到達目的地，需要對飛行過程中的航向角進行約束。這種約束使得飛行器能夠有足夠的距離通過對側(cè)滑角和攻角的控制，把飛行器制導(dǎo)到終點位置。首先將航向角約束定義為在軌跡檢測點上的飛行器實時航向與目標(biāo)地點在經(jīng)緯度平面的夾角。圖2中先求出要研究的飛行器的最大轉(zhuǎn)彎軌跡，在該軌跡上求出各檢測點Pi(i=1,2,…,m)的速度方向與終點連線在經(jīng)緯度平面上的夾角，為速度航向角誤差限制Ψi(i=1,2,…,m)，即

Ψi=ψi-arctan((θi-θf)/(φi-ff))

式中，[θi,φi,ψi]為檢測點Pi位置的經(jīng)緯度坐標(biāo)和航向角，[θf,φf]是終點的經(jīng)緯度坐標(biāo)。

圖3中是兩種飛行器Shuttle與Common Aero Vehicle-HPMARV(CAV_H)根據(jù)最大轉(zhuǎn)彎軌跡計算的航向角約束以及在文獻[4，14]中按照經(jīng)驗設(shè)定的對CAV-H飛行器的航向角約束。

圖2 最大轉(zhuǎn)彎軌跡與檢測點上的航向角

圖3 速度-航向角約束變化

1.3 多禁飛區(qū)約束建模

在飛行器實際飛行中禁飛區(qū)的形狀是多樣的，用解析式來表示很難，而用多個圓柱形禁飛區(qū)組合表示是比較方便的做法。假設(shè)這些禁飛區(qū)都是如圖4所示的圓柱形的區(qū)域，在俯視視角下如圖5所示，高超聲速滑翔飛行器面對的多個小型禁飛區(qū)的復(fù)雜約束情形以及各禁飛區(qū)與飛行器的距離。

圖4 圓柱體禁飛區(qū)示意圖

圖5 高超聲速滑翔飛行器多禁飛區(qū)約束情形

將圓柱形禁飛區(qū)約束簡化為在經(jīng)緯度上的二維約束，針對n個圓柱形禁飛區(qū)中的第i(i=1,2,…,n)個禁飛區(qū)，令禁飛區(qū)的中心經(jīng)緯度坐標(biāo)為(θi,φi)，半徑為Ri，其在經(jīng)緯度上的禁飛區(qū)邊界表示為

飛行器與各禁飛區(qū)中心之間的距離定義為{d1,d2,…,dn}，當(dāng)di≤Ri(i=1,2,…,n)時，飛行器被判定進入第i個禁飛區(qū)，飛行任務(wù)失敗。在未知禁飛區(qū)具體形狀的情況下，飛行器只能探測到禁飛區(qū)的邊界。

2 基于“觸角”預(yù)測的側(cè)滑角瞬變機動制導(dǎo)策略

2.1 基于線性二次調(diào)節(jié)器的縱向制導(dǎo)方法

在機動制導(dǎo)過程中，制導(dǎo)算法分為縱向和橫向兩個制導(dǎo)方法，縱向平面跟蹤方法通過線性二次調(diào)節(jié)器(Linear Quadratic Regulator,LQR)實現(xiàn)[17-19]，即

其中，Δα和Δσ為攻角和側(cè)滑角的調(diào)整量，[Δr,ΔV,Δγ]為實際軌跡與標(biāo)準(zhǔn)軌跡的地心距、速度和航跡角差值，因為是縱向平面的軌跡跟蹤，所以狀態(tài)量上只選擇了這3個關(guān)系密切的變量。

橫向制導(dǎo)方法是一種側(cè)滑角符號轉(zhuǎn)換方法，它與縱向制導(dǎo)策略同時進行，橫向制導(dǎo)策略的周期一般大于縱向制導(dǎo)策略周期，使側(cè)滑角在飛行器縱向制導(dǎo)的跟蹤過程中保持不變，即

ΔTLat=kΔT

其中，ΔTLat與ΔT分別是橫向和縱向制導(dǎo)策略的計算周期，k取整數(shù)，k值越大，橫向機動策略計算的次數(shù)就越少，總機動策略計算的平均時間就越小，k值越小，對禁飛區(qū)的探測效果越好。本文中兩個計算周期都取1 s。

2.2 基于動力學(xué)方程的三“觸角”探測模型

在對高超聲速滑翔飛行器機動制導(dǎo)的過程中，只有最開始解算出的標(biāo)準(zhǔn)軌跡是已知的，而前方飛行通道的實時禁飛狀況由于不能提前探知，需要飛行器采取預(yù)測校正的制導(dǎo)方法進行探測。本文采取的基于“觸手”的方法需要對前方發(fā)出3個觸手來預(yù)測軌跡的飛行狀態(tài)。這種預(yù)測軌跡是由計算機根據(jù)飛行器的動力學(xué)方程積分產(chǎn)生的，其中積分的步長越小，積分產(chǎn)生的預(yù)測路徑越符合實際軌跡。而高超聲速滑翔飛行器動力學(xué)方程較為復(fù)雜，預(yù)測軌跡的積分速度較慢，因而產(chǎn)生的預(yù)測軌跡越少越好。在文獻[14]中，每一個預(yù)測周期產(chǎn)生兩條“觸角”軌跡線，極大地縮減了計算量，但是缺少了直線前進的“觸角”預(yù)測軌跡線，導(dǎo)致在可以無機動地直線飛行時，采取了左右擺動的復(fù)雜制導(dǎo)策略，側(cè)滑角正負(fù)瞬變對控制系統(tǒng)造成的壓力也較大。

本文采取的是前向“觸角”與左右兩側(cè)“觸角”結(jié)合的三“觸角”預(yù)測方式，如圖6，當(dāng)飛行器在前向觸角遇到禁飛區(qū)時，即飛行段OA段，采取三觸角與雙觸角的反饋結(jié)果是基本相同的，飛行器只能在兩側(cè)的觸角中選擇其中的較優(yōu)路線進行機動。當(dāng)飛行器在前向觸角沒有遇到禁飛區(qū)而停止時，即飛行段AB段，飛行器需要判斷前向觸角和另外兩側(cè)觸角的反饋信息，文獻[14]中提到的方法，在類似A點的飛行狀態(tài)時，由于左右觸角給出的反饋信息的實時變化，飛行器會進行左右搖擺的機動，直到因為航向角約束的逐漸減小而最后使軌跡逼近終點；而采取類似于B點的三觸角探測策略，飛行器的前向觸角在AB段沒有禁飛區(qū)約束，因而觸角的終點更貼近終點，飛行器在經(jīng)過短暫的左右機動調(diào)整后可以沿直線飛行到終點。左右兩側(cè)觸角預(yù)測軌跡的側(cè)滑角符號一正一負(fù)，數(shù)值上使用飛行器的最大側(cè)滑角度數(shù)。如果使用度數(shù)較小的側(cè)滑角進行兩側(cè)軌跡的預(yù)測，大多時候會因為機動角度過小而無法使飛行器飛躍禁飛區(qū)；如果在原有三條“觸角”的基礎(chǔ)上，再加兩條小角度的預(yù)測線，會增加計算的時間，為橫向制導(dǎo)策略的選擇增加困難，同時因為高超聲速滑翔飛行器高靈敏的特性而增加控制的難度。

圖6 三“觸角”預(yù)測示意圖

飛行器每一個預(yù)測周期的三個“觸角”的起點是相同的，它們都以飛行器當(dāng)前時刻仿真積分一個時間步長的狀態(tài)為預(yù)測起點，當(dāng)預(yù)測結(jié)束并且獲得相應(yīng)的控制指令后，飛行器也剛剛飛至下一時刻，由于積分時間步長固定且較小，預(yù)測起點與實際狀態(tài)誤差較小。三個“觸角”的終止條件如下：

C1：積分時間超過了理論上飛行的最大任務(wù)時間，本文采取的最大任務(wù)時間是3 000 s；

C2：“觸角”的末端觸碰到了禁飛區(qū)的邊緣或者因為積分步長的原因進入了禁飛區(qū)；

C3：預(yù)測軌跡線上的航向角超出了該速度所對應(yīng)的最大航向角約束；

C4：預(yù)測軌跡順利到達終點位置。

其中，前3個條件屬于禁忌類條件，飛行器應(yīng)避免飛行到這三類“觸角”的終點，第4個終止條件屬于允許類條件，飛行器應(yīng)優(yōu)先選擇這樣終止的“觸手”制導(dǎo)飛行。其中的禁飛區(qū)終止條件是禁忌條件中禁忌程度最高的，而時間終止條件的成立幾率較小，由此得到終止條件所對應(yīng)的“觸角”選擇優(yōu)先級Ka，如表1所示。

表1 終止條件優(yōu)先級

3條軌跡預(yù)測只能以這4種情況結(jié)束，并將終止時的狀態(tài)和終止原因反饋給橫向制導(dǎo)方法，由橫向制導(dǎo)方法選擇三“觸角”的其中一條進行機動制導(dǎo)。

2.3 基于側(cè)滑角瞬變的橫向制導(dǎo)方法

文獻[4]和文獻[14]等提到的側(cè)滑角逆變策略是一種提供側(cè)滑角符號的方法，側(cè)滑角的絕對值是基本不變的，其中文獻[14]在這一基礎(chǔ)上根據(jù)LQR方法其側(cè)滑角的數(shù)值進行了微調(diào)，但都屬于側(cè)滑角逆變策略，即在機動時將飛行器的側(cè)滑角變號，改變飛行器的空氣動力方向。本文將側(cè)滑角歸零這一選擇加入到側(cè)滑角的瞬變策略中，使得飛行器能夠在不需要機動的時候選擇前向飛行，或者在側(cè)滑角進行正負(fù)轉(zhuǎn)換時，作為穩(wěn)定控制系統(tǒng)的過渡狀態(tài)，因此飛行器的側(cè)滑角度數(shù)有3種可能，即

σ=sgn(σ)·σmax,sgn(σ)∈{-1,0,1}

其中，σmax是側(cè)滑角最大值，左右觸角分別對應(yīng)符號函數(shù)sgn(σ)的-1與1。符號函數(shù)sgn(σ)的賦值由觸角反饋的總優(yōu)先級K決定。總優(yōu)先級K由終止條件優(yōu)先級Ka和觸角終止點與終點的距離優(yōu)先級Kb共同決定。因此，每一個觸角的總優(yōu)先級如下式所示：

Ki=-Ka·Kb,i=Left,Mid,Right

其中，Kb根據(jù)觸角終止位置的經(jīng)緯度坐標(biāo)與飛行路徑終點經(jīng)緯度坐標(biāo)的距離排序，距離越小，該觸角的終止位置越靠近終點，該觸角越應(yīng)該被選擇，即根據(jù)這一距離的遠(yuǎn)中近3個等級，Kb分別賦值為{1.1,1.0,0.9}。式中Ki的優(yōu)先級越大，代表該觸角更應(yīng)該被選擇作為下一個時刻的制導(dǎo)指令，即：

為了避免飛行器側(cè)滑角瞬變角度過大或者因為優(yōu)先級的振蕩變化造成的側(cè)滑角振蕩變化而控制系統(tǒng)不能完成這樣的變化，加入一種延時鎖定開關(guān)。如圖 7所示，由側(cè)滑角符號函數(shù)計算出來的結(jié)果進入延時鎖定開關(guān)后，如果開關(guān)關(guān)閉，延時計數(shù)開關(guān)開始計數(shù)，控制系統(tǒng)中的側(cè)滑角不變；如果開關(guān)打開，則重新判斷。

圖7 側(cè)滑角計算流程框圖

當(dāng)前時刻為j時刻，預(yù)測周期為ΔTLat，則控制系統(tǒng)獲得的側(cè)滑角σ′由下式計算得到：

其中，σj、σj-ΔTLat與σj+ΔTLat代表j、j-ΔTLat與j+ΔTLat三個時刻的側(cè)滑角符號函數(shù)輸入值，延時開關(guān)的延時時間最小值根據(jù)飛行器自身特性來決定，本文中的延時時間統(tǒng)一取10 s。在延時過程中，統(tǒng)計式(17)中出現(xiàn)的側(cè)滑角同類振蕩變化次數(shù)，振蕩次數(shù)超過一定限制，則打開延時開關(guān)，并將σj+ΔTLat輸入控制系統(tǒng)。

3 仿真算例與分析

本節(jié)采用CAV-H和Shuttle的高超聲速滑翔飛行器模型在多禁飛區(qū)條件下驗證基于三觸手的側(cè)滑角瞬變機動制導(dǎo)策略。飛行器Shuttle和CAV-H的模型數(shù)據(jù)分別來源于文獻[19]和文獻[20]，表2是兩個飛行器標(biāo)準(zhǔn)軌跡的飛行數(shù)據(jù)[14,19]：

根據(jù)表2初始數(shù)據(jù)獲得兩飛行器的初始標(biāo)準(zhǔn)軌跡高度曲線，如圖8所示。本文兩飛行器的標(biāo)準(zhǔn)軌跡中側(cè)滑角設(shè)為0，攻角的變化幅度較小，導(dǎo)致飛行時間較長。

3.1 未知多禁飛區(qū)情況下的機動制導(dǎo)仿真

為檢測機動制導(dǎo)算法的有效性，設(shè)置了11個半徑為1°的小型圓柱禁飛區(qū)，它們的中心經(jīng)緯度坐標(biāo)如表3所示，組成了復(fù)雜的多禁飛區(qū)模型。

表2 飛行器飛行數(shù)據(jù)

圖8 標(biāo)準(zhǔn)軌跡高度曲線

序號經(jīng)度緯度序號經(jīng)度緯度130-3765-1.523038651.5350-3980345001080055031180-36400

根據(jù)表3的禁飛區(qū)信息，按照圖9的算法流程采取兩種方法求解并比較。一種是利用GPOPS II工具包求解，該求解方法是在離線狀態(tài)下進行的，因此使用已知禁飛區(qū)具體信息求解，同時側(cè)滑角在最小值與最大值構(gòu)成的區(qū)間內(nèi)自由取值。另一方面采取本文的機動制導(dǎo)策略進行在線求解，求解過程只能通過觸角的探測，反饋給制導(dǎo)策略禁飛區(qū)部分輪廓信息，側(cè)滑角的取值只有最大值、最小值和零。如圖10和圖11所示，基于“觸角”的方法與基于GPOPS II的方法都可以使兩種飛行器順利規(guī)避禁飛區(qū)到達終點，同時，如圖12和圖13所示，由于基于“觸角”的方法的側(cè)滑角的選取區(qū)間相對簡單，路徑圖中的軌跡相對平滑。

圖9 仿真算法流程框圖

圖10 CAV-H飛行器規(guī)避禁飛區(qū)路徑

圖11 Shuttle飛行器規(guī)避禁飛區(qū)路徑

圖12 CAV-H控制變量變化

圖14和圖15顯示，基于“觸角”的機動制導(dǎo)策略在規(guī)避禁飛區(qū)的同時，由于速度和高度嚴(yán)格跟蹤標(biāo)準(zhǔn)軌跡，雖然在算法中觸角終止條件沒有判別是否符合3個經(jīng)典路徑約束條件，機動制導(dǎo)的結(jié)果依舊滿足約束條件。

圖14 CAV-H熱流率、動壓、過載變化

圖15 Shuttle熱流率、動壓、過載變化

由圖16中可以看出：飛行過程中的單次預(yù)測周期計算時間是呈下降趨勢的，其中的大小變化主要是由預(yù)測時飛行器位置決定的，飛行器距離禁飛區(qū)、航向角約束或終點越近，計算的時間就越短。

圖16 飛行過程中單次預(yù)測周期計算時間

3.2 模型的魯棒性與計算時間

為了驗證本文基于三“觸角”的側(cè)滑角瞬變機動制導(dǎo)模型，對兩個飛行器進行魯棒性仿真實驗。采用蒙特卡洛方法，對初始狀態(tài)中的緯度、航跡角、航向角以及飛行器的質(zhì)量、面積實施3°、3°、10°和5%、5%的離散，經(jīng)仿真獲得了每個飛行器的1 000次結(jié)果，從圖17和圖18中可以看出雖然有一些仿真路徑的機動幅度較大，但是所有仿真軌跡均能順利規(guī)避本文所設(shè)計的復(fù)雜禁飛區(qū)到達各自的標(biāo)準(zhǔn)軌跡終點。

圖17 CAV-H蒙特卡洛仿真結(jié)果

圖18 Shuttle蒙特卡洛仿真結(jié)果

本文進行實驗時使用的計算機CPU主頻為2.6 GHz，軟件為Matlab。在仿真實驗中，統(tǒng)計每一個仿真實驗的單次預(yù)測周期計算時間的最大值、最小值、平均值，得到圖19和圖20，可以看出，本文提出的機動制導(dǎo)策略在預(yù)測時最大預(yù)測時間不超過0.25 s，而且平均預(yù)測時間集中在0.01 s，完全符合機動制導(dǎo)的時間要求。同時，圖21與圖22展示了仿真中計算得到的熱流率、動壓和過載最大值均滿足表2的約束條件。

圖20 Shuttle單預(yù)測周期計算時間統(tǒng)計圖

圖21 CAV-H熱流率、動壓、過載最大值點集

圖22 Shuttle熱流率、動壓、過載最大值點集

4 結(jié)論

本文通過對基于“觸角”的預(yù)測校正制導(dǎo)方法的改進，綜合航向角約束和側(cè)滑角延時計數(shù)開關(guān)設(shè)計，提出了基于三“觸角”的側(cè)滑角瞬變機動制導(dǎo)策略，解決了高超聲速滑翔飛行器在多禁飛區(qū)信息未知情況下再入機動制導(dǎo)的問題。通過仿真驗證了該方法能夠有效解算出多禁飛區(qū)情況下高超聲速滑翔飛行器機動制導(dǎo)的在線飛行方案，通過蒙特卡洛仿真實驗驗證，突出體現(xiàn)了該方法計算時間上能夠滿足機動制導(dǎo)的要求。