一種基于參數(shù)在線辨識的最優(yōu)再入制導(dǎo)算法*

尤志鵬 楊 勇 鄭宏濤 劉 剛

中國運載火箭技術(shù)研究院，北京 100076

0 引言

近年來以X-37B、IXV等為代表的升力式可重復(fù)使用飛行器得到越來越廣泛的關(guān)注。這種飛行器既能夠?qū)崿F(xiàn)在軌長時間運行，又能夠?qū)崿F(xiàn)水平著陸，并且具備良好的可重復(fù)使用能力，有力帶動了先進(jìn)制導(dǎo)控制等技術(shù)的發(fā)展，成為近年來的研究熱點[1-2]。

良好的再入制導(dǎo)律設(shè)計對提升飛行性能具有顯著作用。基于標(biāo)準(zhǔn)軌跡的再入制導(dǎo)是當(dāng)前應(yīng)用最廣泛的升力式再入制導(dǎo)方法，眾多學(xué)者對基于標(biāo)準(zhǔn)軌跡的再入制導(dǎo)算法進(jìn)行了深入研究。文獻(xiàn)[3]通過設(shè)計阻力加速度-能量參考曲線和航跡偏航角-能量參考軌跡，同時規(guī)劃飛行器縱向運動和橫向運動，提出衍化加速度再入制導(dǎo)律?？稍诰€調(diào)整標(biāo)準(zhǔn)阻力加速度剖面，具有較強(qiáng)的自適應(yīng)能力，可以應(yīng)用于大橫程再入飛行的情況。文獻(xiàn)[4]在速度-高度剖面內(nèi)設(shè)計參考軌跡，可在線預(yù)測待飛航程及飛行時間，可實現(xiàn)時間協(xié)同再入制導(dǎo)。文獻(xiàn)[5]設(shè)計了一種基于阻力加速度倒數(shù)-能量的飛行剖面，利用三次樣條描述阻力加速度倒數(shù)剖面并周期性更新，具有良好的工程應(yīng)用潛力。文獻(xiàn)[6]設(shè)計了一種通過跟蹤阻力加速度-能量剖面，并在末端引入航跡傾角控制實現(xiàn)航跡傾角終端要求，能夠得到較高的制導(dǎo)精度。

飛行器制導(dǎo)控制系統(tǒng)正在走向智能化，“會學(xué)習(xí)”成為其顯著特征[7]。而飛行器感知能力是“會學(xué)習(xí)”的重要能力支撐，實現(xiàn)再入?yún)?shù)不確定性在線辨識，對提升再入制導(dǎo)效果具有較大影響。不確定性參數(shù)在線辨識技術(shù)經(jīng)過多年的發(fā)展，已經(jīng)形成較為完備的方法體系。文獻(xiàn)[8]研究了基于卡爾曼濾波的氣動參數(shù)辨識問題，對比了擴(kuò)展卡爾曼濾波算法(EKF)、無跡卡爾曼濾波(UKF)和簡化UKF三種濾波估計算法，結(jié)果顯示三種算法均可實現(xiàn)較高估計精度。文獻(xiàn)[9]利用擴(kuò)展卡爾曼濾波對飛行過程中三維風(fēng)場進(jìn)行辨識，實現(xiàn)了很高的辨識精度。文獻(xiàn)[10]利用飛試數(shù)據(jù)對發(fā)射段風(fēng)場進(jìn)行辨識，辨識得到的風(fēng)場與氣象觀測得到的風(fēng)場具有良好的相似度。文獻(xiàn)[11]在研究火星再入的過程中，利用嵌入式大氣數(shù)據(jù)系統(tǒng)(FADS)與慣導(dǎo)系統(tǒng)耦合，成功辨識出火星再入過程中大氣密度、動壓、風(fēng)速等參數(shù)。隨著器載傳感器及計算機(jī)技術(shù)發(fā)展，參數(shù)辨識輔助制導(dǎo)控制系統(tǒng)設(shè)計得到了較多關(guān)注。

但是，文獻(xiàn)[3-6]剖面規(guī)劃均更著重于剖面可行而非最優(yōu)。文獻(xiàn)[8]研究了氣動參數(shù)辨識問題，文獻(xiàn)[9-10]主要研究飛行環(huán)境參數(shù)的辨識問題，但對于大氣和氣動參數(shù)均存在偏差的參數(shù)辨識問題，文獻(xiàn)[8-10]均未開展研究，文獻(xiàn)[11]雖然能夠?qū)崿F(xiàn)對火星再入過程中大氣及氣動參數(shù)的聯(lián)合辨識，但算法設(shè)計較為復(fù)雜，精度尚待提升。借鑒參數(shù)優(yōu)化技術(shù)的設(shè)計思想，提升剖面規(guī)劃最優(yōu)性，并通過參數(shù)辨識，在線獲取不確定性參數(shù)偏差，可改善再入制導(dǎo)效果，但當(dāng)前，對基于氣動及大氣參數(shù)聯(lián)合在線辨識并通過跟蹤最優(yōu)飛行剖面實現(xiàn)再入制導(dǎo)的研究仍較少，相關(guān)算法尚待進(jìn)一步開發(fā)。

本文提出一種基于參數(shù)在線辨識的最優(yōu)再入剖面規(guī)劃與制導(dǎo)算法。首先在速度-高度剖面內(nèi)設(shè)計分段解析飛行剖面，通過參數(shù)優(yōu)化設(shè)計在線迭代獲取飛行剖面分段點及設(shè)計參數(shù)，獲得滿足總吸熱量最小的參考飛行剖面；其次，利用反饋線性化的方法求解制導(dǎo)指令并通過航跡傾角偏差進(jìn)行修正，實現(xiàn)飛行剖面跟蹤，并通過觀測的過載、動壓、俯仰角等參數(shù)，在線辨識飛行過程中狀態(tài)及氣動參數(shù)和大氣環(huán)境參數(shù)，利用辨識結(jié)果修正模型及制導(dǎo)指令中的參數(shù)偏差，提高制導(dǎo)精度；最后，通過待飛航程與待飛航程預(yù)測值偏差決定是否進(jìn)行剖面重新規(guī)劃，并通過傾側(cè)角符號反轉(zhuǎn)實現(xiàn)橫向制導(dǎo)。與前人研究相比，本文所規(guī)劃的飛行剖面不但能夠滿足航程要求，且總吸熱量最小，對飛行過程中參數(shù)偏差具有良好的辨識能力，可為升力式再入飛行器飛行軌跡設(shè)計及指標(biāo)論證提供參考。

1 最優(yōu)飛行剖面規(guī)劃

1.1 再入運動方程

假設(shè)地球是均質(zhì)圓球，三維質(zhì)點再入運動無量綱方程為：

(1)

L=ρ(VcV)2SrefCL/(2mg0)

(2)

D=ρ(VcV)2SrefCD/(2mg0)

(3)

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

1.2 再入約束

再入過程約束主要包括動壓約束、熱流約束、過載約束，以及平衡滑翔約束。如下所示：

(4)

(5)

n=|Lcosα+Dsinα|≤nmax

(6)

(7)

末端約束主要包含末端高度、末端速度、末端經(jīng)緯度約束，表示為h(Vf)=hf，θ(Vf)=θf，φ(Vf)=φf，這里hf，θf和φf分別表示末端飛行高度、末端經(jīng)度和末端緯度。

1.3 剖面擬合

在高度-速度剖面內(nèi)，可以將高度擬合為速度的多項式或多項式的分段。分段擬合能夠得到性質(zhì)更好的飛行剖面，通過合理的段數(shù)劃分，使擬合得到的飛行剖面變化平穩(wěn)，設(shè)計靈活性更大。本文在高度-速度剖面內(nèi)將高度擬合為兩段三次多項式形式。即

(8)

式中：V0和Vf表示剖面規(guī)劃的初始點及末端點，href(V)表示不同速度下參考高度待確定系數(shù)，一共有8個，且有一個分段點Vm待確定。高度對速度的導(dǎo)數(shù)如式(9)所示

(9)

由式(9)可見，在初始速度、高度、航跡傾角給定的條件下，初始點高度對速度導(dǎo)數(shù)即可確定。同樣，對于末端點，高度和速度滿足終端約束時，若末端航跡傾角確定，則末端高度對速度的導(dǎo)數(shù)也可通過上式確定。同時，為保證兩段參考軌跡過渡平穩(wěn)，需保證兩段剖面連接點處高度相等，且高度隨速度變化率相等。至此，約束條件確定了6項。而待求解參數(shù)包含8個擬合系數(shù)及一個切換點，需要進(jìn)一步考慮其他約束，從而完全確定參考剖面。

1.4 剖面參數(shù)確定

通過1.3節(jié)方法，可以獲得一組飛行剖面，待確定參數(shù)有3個。強(qiáng)約束僅剩下航程約束，為求解出所有待定參數(shù)，可選擇優(yōu)化某種指標(biāo)函數(shù)，使之達(dá)到最優(yōu)化。本文對總吸熱量進(jìn)行優(yōu)化，總吸熱量表達(dá)式為

(10)

由于本文飛行剖面設(shè)計過程中采用飛行速度作為自變量，因此將總吸熱量轉(zhuǎn)換成以速度為自變量的積分形式，即

(11)

使總吸熱量最小，即使式(11)達(dá)到最小。在設(shè)計剖面下，有

(12)

式中：γref表示標(biāo)準(zhǔn)飛行剖面對應(yīng)的航跡傾角。

再入過程中，待飛航程隨速度變化的表達(dá)式為

(13)

因而航程約束可以表達(dá)為

(14)

(15)

式中：hmin(v1),hmax(v1)分別表示速度v1下速度-高度再入走廊下邊界和上邊界。hmin(v2),hmax(v2)分別表示速度v2下速度-高度再入走廊下邊界和上邊界。

求解上述優(yōu)化問題，可得兩個設(shè)計點對應(yīng)的高度及兩個剖面的分界點速度。待求解問題中，式(11)和(15)均存在積分項，可以利用數(shù)值求積分公式轉(zhuǎn)化成解析表達(dá)式形式[12]。記V0

δV=(Vf-V0)/N,Vi=Vi-1+δV

(16)

記

(17)

(18)

則式(15)可表達(dá)為

(19)

記

(20)

(21)

(22)

這里Vmin和Vmax是Vm的搜索邊界。

可采用非線性數(shù)值規(guī)劃算法求解上述問題，本文通過序列二次規(guī)劃算法對該問題進(jìn)行求解。

2 最優(yōu)飛行剖面跟蹤制導(dǎo)

2.1 傾側(cè)角指令

利用反饋線性化進(jìn)行傾側(cè)角指令求解，首先將高度對速度二次微分，可以得到

(23)

式中：

(24)

(25)

(26)

式中：j表示第j個制導(dǎo)周期。標(biāo)準(zhǔn)軌跡跟蹤可通過設(shè)計二階環(huán)節(jié)跟蹤器實現(xiàn)，即

(27)

式中：ζ和ω分別是阻尼比和自然頻率。

通過式(27)形成的制導(dǎo)指令，可實現(xiàn)對標(biāo)準(zhǔn)軌跡的跟蹤，但仿真發(fā)現(xiàn)，其收斂較慢，導(dǎo)致剖面更新次數(shù)增多，制導(dǎo)誤差較大。分析式(12)可知，在標(biāo)準(zhǔn)H-V飛行剖面下，在任意速度下，航跡傾角與高度對速度變化率均是一一對應(yīng)的。在式(27)形成的制導(dǎo)律的基礎(chǔ)上，進(jìn)一步調(diào)整傾側(cè)角，使航跡傾角迅速收斂于標(biāo)準(zhǔn)飛行剖面對應(yīng)的航跡傾角γref，即迅速實現(xiàn)對dhref/dV的跟蹤，可改善制導(dǎo)的動態(tài)過程，避免跟蹤過程中的振蕩。使實際航跡傾角收斂于標(biāo)準(zhǔn)飛行剖面對應(yīng)的航跡傾角，即實現(xiàn)

(28)

即

(29)

式(29)中，γ即在制導(dǎo)律(27)作用下產(chǎn)生的航跡傾角。進(jìn)一步調(diào)整后的傾側(cè)角記為σcmd。由式(29)可得

(30)

式中：Lref和Dref表示參考剖面對應(yīng)的無量綱升力和阻力，K為大于0的正數(shù)，M表達(dá)為

(31)

式(30)得到的制導(dǎo)律可提高收斂速度，當(dāng)偏差較小時，M影響較小，當(dāng)實際飛行剖面收斂于標(biāo)準(zhǔn)飛行剖面時，式(30)產(chǎn)生的制導(dǎo)指令與式(27)相同。

2.2 不確定性參數(shù)在線辨識

從式(30)可見，需要反饋的參數(shù)中，升力、阻力、航跡傾角受氣動參數(shù)偏差擾動及大氣密度、風(fēng)速擾動較大。為進(jìn)一步提升制導(dǎo)精度，可對偏差量進(jìn)行在線辨識。影響較大的待辨識參數(shù)主要包括升力系數(shù)偏差ΔCL、阻力系數(shù)偏差ΔCD、大氣密度偏差Δρ和攻角偏差Δα。將偏差導(dǎo)數(shù)視為白噪聲增廣進(jìn)入縱向狀態(tài)方程即形成連續(xù)的狀態(tài)方程，其中，ε1,…,ε6是互不相關(guān)高斯白噪聲，L*，D*是受到擾動后的升力和阻力，分別計算如下

(32)

(33)

(34)

將(32)離散即得到濾波估計的狀態(tài)方程。觀測量包括高度、法向過載、軸向過載以及俯仰角，F(xiàn)ADS等技術(shù)的發(fā)展使得動壓可測。因此可得到觀測方程，如式(35)所示。

(35)

系統(tǒng)噪聲協(xié)方差矩陣記為Qk，測量噪聲矩陣為Rk。至此建立了辨識擾動不確定性參數(shù)的狀態(tài)方程和觀測方程，它們均是非線性的，利用EKF算法進(jìn)行參數(shù)辨識，由于偏差預(yù)測模型準(zhǔn)確度較差，為得到更準(zhǔn)確的辨識結(jié)果，本文采用漸消記憶的擴(kuò)展卡爾曼濾波算法實現(xiàn)不確定參數(shù)的辨識[13]。至此，可得到每個制導(dǎo)周期內(nèi)不確定性參數(shù)辨識的辨識結(jié)果，利用辨識得到的結(jié)果對運動模型、傾側(cè)角指令及裝訂的攻角指令進(jìn)行修正，提高制導(dǎo)精度。

2.3 橫向制導(dǎo)指令

橫向制導(dǎo)通過改變傾側(cè)角符號實現(xiàn)，表達(dá)如下

(36)

式中：航向角偏差計算如下

Δψ=ψ-ψLOS

(37)

ψLOS即當(dāng)前位置至目標(biāo)點的理想視線角，計算如下

(38)

至此，已經(jīng)完成基于不確定參數(shù)在線辨識的最優(yōu)再入剖面規(guī)劃與制導(dǎo)算法設(shè)計?；具^程如圖1所示。圖1中，R表示待飛航程，Rpred表示待飛預(yù)測航程。

圖1 算法流程

3 仿真校驗

基于X-33模型進(jìn)行仿真驗證，該飛行器質(zhì)量37363 kg，參考面積149.4 m2，仿真再入過程中，攻角(°)按馬赫數(shù)裝訂，如下

(39)

首先對4種不同航程下的算例進(jìn)行仿真驗證。它們初始再入位置不同，但是交班點經(jīng)緯度均為(112°，42°)，交班點高度為24.87 km。升力系數(shù)、阻力系數(shù)、大氣密度、攻角偏差如表2所示。

表1 初始條件

表2 參數(shù)偏差

情況1～4仿真結(jié)果如圖2～4所示，圖2表示經(jīng)緯度變化曲線，可見均能精確到達(dá)預(yù)定目標(biāo)。圖3展示了高度隨速度變化曲線，可直觀看出，若達(dá)到較遠(yuǎn)航程，相同速度下，需要飛行高度更高。情況1由于預(yù)定航程較短，曲線已非常接近約束下邊界，而情況4曲線由于航程較長，已經(jīng)非常接近于平衡滑翔約束。圖4是考慮了偏差及航跡傾角反饋后的制導(dǎo)指令。

圖2 地面軌跡

圖3 速度-高度曲線

圖4 制導(dǎo)指令

圖5 攻角偏差及航跡傾角辨識誤差

圖6 升阻力系數(shù)偏差及大氣密度偏差辨識誤差

圖5～6表示情況1～4對表2所示偏差的辨識結(jié)果，可見對航跡傾角、攻角、升阻力系數(shù)偏差辨識結(jié)果均能快速收斂，但對大氣密度偏差辨識誤差逐漸增大，主要原因是大氣密度偏差隨飛行高度降低呈指數(shù)式增加，濾波結(jié)果呈現(xiàn)出滯后效應(yīng)。進(jìn)一步分析可發(fā)現(xiàn)該辨識誤差占大氣密度的比例呈逐漸收斂至0的趨勢，因此該辨識誤差對制導(dǎo)精度的影響逐漸降低。

為全面檢驗算法性能，引入初始參數(shù)偏差和狀態(tài)參數(shù)偏差進(jìn)行蒙特卡洛打靶仿真。蒙特卡洛打靶加入的隨機(jī)偏差項如表3所示。

表3 蒙特打靶偏差項

加入不確定性辨識環(huán)節(jié)并根據(jù)式(30)產(chǎn)生制導(dǎo)指令，落點分布如圖7所示，可見交班點在目標(biāo)點附近，且偏差散步較小，展現(xiàn)出良好的適應(yīng)性。若不對不確定性參數(shù)進(jìn)行在線辨識，且制導(dǎo)指令根據(jù)式(27)產(chǎn)生，則偏差散布較大，且偏差較大，如圖8所示。

圖7 加入辨識及附加反饋后的打靶仿真交班點

圖8 未加入辨識和附加反饋的打靶仿真交班點

4 結(jié)論

對于最優(yōu)飛行剖面規(guī)劃及參數(shù)不確定條件下的制導(dǎo)問題，本文利用分段多項式擬合速度-高度飛行剖面，并利用序列二次規(guī)劃求解得到滿足總吸熱量最優(yōu)的擬合參數(shù)。通過反饋線性化方法并附加航跡傾角偏差反饋形成的制導(dǎo)指令能夠保證對所設(shè)計飛行剖面的跟蹤，且飛行過程中的參數(shù)偏差可以通過漸消記憶卡爾曼濾波進(jìn)行準(zhǔn)確辨識，提高了制導(dǎo)精確性。仿真表明本文算法適應(yīng)性強(qiáng)，不確定性參數(shù)在線辨識并修正制導(dǎo)指令可有效提升制導(dǎo)精度。