尤志鵬,楊勇,劉剛,曹曉瑞,鄭宏濤

中國(guó)運(yùn)載火箭技術(shù)研究院,北京 100076

近年來以X-37B、IXV、X-33等為代表的升力式可重復(fù)使用空天飛行器得到越來越廣泛的關(guān)注,促進(jìn)了低成本進(jìn)出太空技術(shù)的發(fā)展,為先進(jìn)制導(dǎo)控制等相關(guān)技術(shù)的進(jìn)步和試驗(yàn)驗(yàn)證提供了可靠平臺(tái),成為近年來研究的熱點(diǎn)[1-2]。為應(yīng)對(duì)再入過程中熱流、過載等嚴(yán)苛的約束,將飛行器精確導(dǎo)引至目標(biāo)點(diǎn),以高精度、高自主性、高魯棒性為特征的新一代再入制導(dǎo)算法得到持續(xù)而深入的研究。但由于飛行過程中具有高動(dòng)態(tài)、參數(shù)不確定性大、飛行工況變化劇烈的特點(diǎn)[3],偏差在線修正能力強(qiáng)、自適應(yīng)性高的再入制導(dǎo)算法仍待進(jìn)一步開發(fā)。

再入制導(dǎo)從算法上主要包含標(biāo)準(zhǔn)軌跡制導(dǎo)和預(yù)測(cè)校正制導(dǎo)[4]。標(biāo)準(zhǔn)軌跡制導(dǎo)算法在航天飛機(jī)再入制導(dǎo)律設(shè)計(jì)中取得了巨大的成功,這種制導(dǎo)算法需要事先或在線生成參考軌跡剖面,通過偏差反饋形成制導(dǎo)指令,計(jì)算量小,容易滿足再入走廊要求,但這種算法對(duì)初值、參數(shù)偏差、控制飽和等因素比較敏感,制導(dǎo)偏差較大。預(yù)測(cè)校正制導(dǎo)算法在制導(dǎo)精度方面具有優(yōu)勢(shì),自適應(yīng)能力強(qiáng),魯棒性好,但是由于預(yù)測(cè)過程通常需要進(jìn)行多次彈道積分,計(jì)算資源需求較多。降低預(yù)測(cè)過程所需資源,提高預(yù)測(cè)快速性,對(duì)預(yù)測(cè)校正制導(dǎo)的應(yīng)用將產(chǎn)生諸多有利影響。

早期預(yù)測(cè)校正算法大多設(shè)計(jì)參數(shù)較多,迭代過程較為復(fù)雜。Powell[5]為火星探測(cè)器設(shè)計(jì)了預(yù)測(cè)校正再入制導(dǎo)律,以四階次龍格庫塔法積分三自由度質(zhì)點(diǎn)動(dòng)力學(xué)方程,尋求滿足約束條件的傾側(cè)角幅值指令和傾側(cè)方向指令。Draper實(shí)驗(yàn)室為K-1軌道飛行器設(shè)計(jì)了預(yù)測(cè)校正再入制導(dǎo)律[6],用于解算傾側(cè)角指令和反轉(zhuǎn)時(shí)刻,為降低迭代復(fù)雜度,傾側(cè)角指令僅進(jìn)行一次反轉(zhuǎn)。Lu[7]對(duì)數(shù)值預(yù)測(cè)-校正算法進(jìn)行改進(jìn),使數(shù)值預(yù)測(cè)校正制導(dǎo)每次僅需要迭代獲取傾側(cè)角剖面,使之滿足航程要求,并通過橫程門限判斷是否需要反號(hào),降低了預(yù)測(cè)復(fù)雜度,提高了預(yù)測(cè)效率。Xue和Lu[8]進(jìn)一步考慮舵面偏轉(zhuǎn)速度和加速度限制,并改進(jìn)了橫程控制,使文獻(xiàn)[7]所提出的算法更加完善,最終形成一種適應(yīng)于不同升阻比飛行器且具有較強(qiáng)軌跡約束滿足能力的通用預(yù)測(cè)校正制導(dǎo)算法[9-11]。但是,為了得到滿足航程約束的傾側(cè)角剖面,在每個(gè)制導(dǎo)周期,上述制導(dǎo)算法均需要在縱向運(yùn)動(dòng)平面內(nèi)執(zhí)行多次彈道積分,算法效率有待進(jìn)一步提高。

為降低傳統(tǒng)數(shù)值預(yù)測(cè)校正制導(dǎo)算法計(jì)算復(fù)雜度,提升算法效率,胡軍等[12-13]提出基于一階特征模型的全系數(shù)自適應(yīng)預(yù)測(cè)校正制導(dǎo)律。外環(huán)利用全系數(shù)自適應(yīng)方法估計(jì)一階特征模型系數(shù)并求取制導(dǎo)修正量,實(shí)現(xiàn)對(duì)規(guī)劃彈道的持續(xù)修正,內(nèi)環(huán)則采用短周期的彈道跟蹤制導(dǎo)。這種內(nèi)外環(huán)結(jié)合的制導(dǎo)算法對(duì)預(yù)測(cè)更新頻率要求有所降低,避免了多次彈道積分。陳萬春等[14-16]對(duì)再入彈道解析解進(jìn)行深入研究,以旋轉(zhuǎn)地球?yàn)閰⒖枷?基于三維再入彈道解析解規(guī)劃參考軌跡,并提出彈道阻尼控制技術(shù)抑制彈道振蕩,航跡規(guī)劃快速、靈活,飛行軌跡平穩(wěn),終端約束及過程約束均可精確滿足。這2類制導(dǎo)算法均具有良好的收斂性和工程應(yīng)用價(jià)值,但前者主要應(yīng)用于飛船這類低升阻比飛行器,而后者主要應(yīng)用于高升阻比飛行器,對(duì)X-37B、IXV這類中等升阻比飛行器,2類算法均研究較少。

不同于上述可直接產(chǎn)生制導(dǎo)指令的算法,通過參數(shù)優(yōu)化及數(shù)值迭代方式獲取飛行剖面,并據(jù)此產(chǎn)生制導(dǎo)指令的制導(dǎo)算法也得到關(guān)注。張慶振等[17]將控制剖面參數(shù)化,通過將傾側(cè)角走廊上下界進(jìn)行插值得到規(guī)劃的傾側(cè)角剖面,將規(guī)劃航程擬合為插值系數(shù)的一階線性函數(shù),采用校正信息在線收集與預(yù)測(cè)方法,對(duì)擬合模型進(jìn)行辨識(shí),從而達(dá)到快速預(yù)測(cè)的目的,但是這種算法得到的飛行軌跡不平滑,再入最大過載偏大。沈作軍等[18-19]利用四維多項(xiàng)式擬合速度-高度變化曲線,通過確定擬合系數(shù)得到可行的再入縱向運(yùn)動(dòng)飛行剖面,該算法設(shè)計(jì)簡(jiǎn)單,收斂效果好,但為得到滿足約束的參考剖面,仍需要進(jìn)行多次迭代。王肖等[20]在速度-高度剖面內(nèi)設(shè)計(jì)參考軌跡,可在線預(yù)測(cè)待飛航程及飛行時(shí)間,可實(shí)現(xiàn)時(shí)間協(xié)同再入制導(dǎo),但飛行剖面較為復(fù)雜,制導(dǎo)誤差偏大。

相對(duì)于傳統(tǒng)數(shù)值預(yù)測(cè)校正算法,上述參數(shù)化制導(dǎo)律設(shè)計(jì)算法可有效降低生成制導(dǎo)指令所需要的計(jì)算資源,但制導(dǎo)指令的產(chǎn)生又依賴于當(dāng)前飛行狀態(tài)、所處飛行環(huán)境及規(guī)劃的飛行剖面,因此通過參數(shù)在線辨識(shí)獲取當(dāng)前氣動(dòng)力系數(shù)及大氣環(huán)境偏差,對(duì)提升制導(dǎo)精度和魯棒性具有顯著作用。Schoenenberger等[21]在研究火星再入的過程中,利用嵌入式大氣數(shù)據(jù)系統(tǒng)與慣導(dǎo)系統(tǒng)耦合,成功辨識(shí)出火星再入過程中大氣密度、動(dòng)壓、風(fēng)速、攻角、側(cè)滑角等參數(shù)。邵干等[22]利用啟發(fā)式優(yōu)化算法,成功辨識(shí)了無人機(jī)氣動(dòng)參數(shù),精度較高。陳爾康等[23]通過優(yōu)化傳感器布置方案,提出正交三角分解更新到達(dá)代價(jià)的滾動(dòng)時(shí)域估計(jì)(MHE-QR)算法實(shí)現(xiàn)對(duì)彈性飛行器狀態(tài)/參數(shù)聯(lián)合估計(jì),精度和實(shí)時(shí)性均較高。但是,當(dāng)前對(duì)升力式可重復(fù)使用飛行器再入過程的氣動(dòng)及大氣環(huán)境辨識(shí)研究仍較少。

為提升預(yù)測(cè)校正制導(dǎo)響應(yīng)速度,實(shí)現(xiàn)對(duì)再入飛行時(shí)間的控制,本文將在速度-高度剖面內(nèi),在每個(gè)制導(dǎo)周期,將飛行高度擬合為飛行速度的四次多項(xiàng)式,在再入走廊約束下,考慮當(dāng)前點(diǎn)狀態(tài)及終端狀態(tài)約束,通過卡爾曼濾波得到滿足航程要求的飛行剖面,并基于當(dāng)前狀態(tài)及迭代得到的飛行剖面生成該制導(dǎo)周期的傾側(cè)角幅值指令,橫向制導(dǎo)采用根據(jù)航向誤差走廊確定制導(dǎo)指令符號(hào)的方式實(shí)現(xiàn)。為提高制導(dǎo)精度,削弱氣動(dòng)系數(shù)、大氣環(huán)境系數(shù)偏差等因素造成的影響,再入過程中,對(duì)升阻力系數(shù)偏差、大氣密度偏差、攻角偏差等進(jìn)行在線辨識(shí),并利用辨識(shí)結(jié)果修正制導(dǎo)指令,從而增強(qiáng)制導(dǎo)算法對(duì)參數(shù)偏差的自適應(yīng)能力。該算法制導(dǎo)精度高,計(jì)算速度快,再入飛行時(shí)間可控,實(shí)現(xiàn)簡(jiǎn)單,在再入制導(dǎo)律設(shè)計(jì)、協(xié)同再入等方面具有較強(qiáng)的工程應(yīng)用潛力。

1 運(yùn)動(dòng)方程及約束

1.1 運(yùn)動(dòng)方程

假設(shè)地球是均質(zhì)圓球,三維質(zhì)點(diǎn)再入運(yùn)動(dòng)無量綱方程為

(1)

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

(2)

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

(3)

其中:ρ為大氣密度；Sref和m分別為參考面積和飛行器質(zhì)量；CL和CD分別為升力系數(shù)和阻力系數(shù),與攻角α有關(guān),而攻角通常按速度進(jìn)行裝訂。

1.2 飛行約束

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

(4)

(5)

n=|Lcosα+Dsinα|≤nmax

(6)

(7)

(8)

(9)

(10)

(11)

根據(jù)式(7),可利用牛頓迭代求解出滿足平衡滑翔條件的軟約束上邊界h_max(V),構(gòu)成再入走廊。

末端約束主要包含末端高度、末端速度、末端經(jīng)緯度約束,在速度-高度剖面下,速度是自變量,可將其取為V0～Vf,其他3個(gè)末端約束表達(dá)如下:

h(Vf)=hf

(12)

θ(Vf)=θf

(13)

φ(Vf)=φf

(14)

式中:hf、θf和φf分別表示再入交班點(diǎn)高度、末端經(jīng)度和末端緯度。

2 制導(dǎo)算法設(shè)計(jì)

與傳統(tǒng)預(yù)測(cè)再入制導(dǎo)算法類似,本文將攻角裝訂為馬赫數(shù)的函數(shù),通過傾側(cè)角幅值指令控制縱向航程,通過傾側(cè)角符號(hào)反轉(zhuǎn)控制橫向偏差。再入過程中,無量綱速度具有良好的單調(diào)性,可替代無量綱時(shí)間作為自變量。與傳統(tǒng)預(yù)測(cè)制導(dǎo)不同,本文首先將飛行高度隨速度變化擬合為四次多項(xiàng)式,在每個(gè)制導(dǎo)周期,通過卡爾曼濾波迭代獲得滿足航程約束的飛行剖面擬合結(jié)果,并基于擬合結(jié)果生成傾側(cè)角指令。在每次制導(dǎo)周期,僅進(jìn)行一次一維函數(shù)積分運(yùn)算,因此預(yù)測(cè)速度大大加快。采用基于航向角偏差走廊確定傾側(cè)角符號(hào)的方式實(shí)現(xiàn)橫向制導(dǎo)。

2.1 縱向高度-速度飛行剖面擬合

本文在速度-高度剖面內(nèi)設(shè)計(jì)再入制導(dǎo)律,將飛行高度擬合為速度的多項(xiàng)式形式。在第i個(gè)制導(dǎo)周期,已知狀態(tài)主要包括當(dāng)前點(diǎn)的速度Vi、高度hi、航跡傾角γi,交班點(diǎn)的速度Vf、高度hf；在高度-速度剖面內(nèi),有

(15)

h=a4V4+a3V3+a2V2+a1V+a0

(16)

式中:a0～a4為系數(shù)。

在高度-速度剖面內(nèi),待飛航程隨速度變化表達(dá)式為

(17)

在高度-速度剖面下飛行航程與待飛航程相等,即

(18)

式中:R為待飛航程,即當(dāng)前點(diǎn)至目標(biāo)點(diǎn)對(duì)應(yīng)的航程。

根據(jù)式(16),有

(19)

將式(19)代入式(18),可知速度-高度擬合表達(dá)式?jīng)Q定是否在該剖面下的飛行航程與待飛航程相等。因此,在第i個(gè)制導(dǎo)周期,考慮引入當(dāng)前點(diǎn)與交班點(diǎn)的速度中點(diǎn)Vm_i,且Vm_i=(Vi+Vf)/2。

假設(shè)該點(diǎn)對(duì)應(yīng)的高度為hm_i,記

S(hm_i)=-R+

(20)

則滿足航程約束即選擇合適的hm_i,得到的速度-高度剖面使得S(hm_i)=0。

2.2 縱向制導(dǎo)指令

為確定制導(dǎo)指令,需要在每個(gè)制導(dǎo)周期獲得合適的hm_i,從而確定飛行剖面并在此基礎(chǔ)上得到傾側(cè)角幅值。傳統(tǒng)預(yù)測(cè)校正制導(dǎo)算法通常通過迭代的方式獲取該參數(shù),通常每個(gè)制導(dǎo)周期需要迭代2～10次。但是由式(20)可見,多次迭代計(jì)算需要多次積分,不利于實(shí)時(shí)生成制導(dǎo)指令。

至此,預(yù)測(cè)校正問題實(shí)際上轉(zhuǎn)化為一個(gè)單參數(shù)單等式約束規(guī)劃問題,其中hm_i和S(hm_i)存在特定的隱藏函數(shù)對(duì)應(yīng)關(guān)系。因此,本文通過卡爾曼濾波,對(duì)每個(gè)制導(dǎo)周期內(nèi)滿足式(20)的hm_i進(jìn)行辨識(shí)。辨識(shí)得到的hm_i可用于得到該制導(dǎo)周期對(duì)應(yīng)的飛行剖面,從而產(chǎn)生該制導(dǎo)周期所需要的制導(dǎo)指令。通過卡爾曼濾波,能夠破除周期內(nèi)需要多次迭代、有效軌跡預(yù)測(cè)信息存在噪聲的缺陷,模型信息可以跨周期地傳遞與共享,每個(gè)制導(dǎo)周期僅需要積分式(20),以獲取觀測(cè)預(yù)測(cè)值及其觀測(cè)矩陣。相比于通過數(shù)值迭代方式獲取速度-高度剖面,避免了多次迭代帶來的積分次數(shù)增加,因而實(shí)時(shí)性得到增強(qiáng)。

建模過程如下：

狀態(tài)方程：

hm_i+1=hm_i+εi

(21)

觀測(cè)方程：

Zi=S(hm_i)+νi

(22)

式中：εi及νi(i=1,2,3,…)為互不相關(guān)零均值高斯白噪聲,協(xié)方差分別記為Qi,Wi。

濾波過程如下：

為獲得hm_i,使得S(hm_i)=0,因此可以認(rèn)為觀測(cè)Zi+1≡0,同時(shí),為保證剖面不超過再入走廊約束,需要確保在當(dāng)前速度中點(diǎn)下,hm_i滿足再入走廊約束,且距離再入走廊邊界具有一定的距離,該無量綱距離取為κ=5×10-4,可得到修正后的辨識(shí)結(jié)果如下

(23)

(24)

式中:δ為一小攝動(dòng)量,取為10-5。

(25)

則

(26)

而

(27)

進(jìn)而

(28)

式中:ηj,ωj分別是Gauss-Legendre節(jié)點(diǎn)和節(jié)點(diǎn)權(quán)值系數(shù),參考文獻(xiàn)[24],通常六階Gauss-Legendre求積公式即可得到相當(dāng)高的精度。

當(dāng)hm_i+1確定后,即可通過式(29)求解該制導(dǎo)周期對(duì)應(yīng)的速度-高度剖面擬合系數(shù)

(29)

為獲得傾側(cè)角指令,可將飛行高度對(duì)速度求二階導(dǎo)數(shù),得

(30)

其中

(31)

(32)

(33)

這里下標(biāo)i表示第i個(gè)制導(dǎo)周期,其中

(34)

當(dāng)前飛行狀態(tài)接近終端狀態(tài)時(shí),打靶仿真發(fā)現(xiàn),式(29)可能出現(xiàn)奇異,導(dǎo)致擬合不準(zhǔn)確。為解決該問題,當(dāng)前飛行速度與終端速度之差小于100 m/s時(shí),停止預(yù)測(cè)校正制導(dǎo),不再更新飛行剖面,改為跟蹤最后一次產(chǎn)生的飛行剖面,即

(35)

Lcosσ′=Lcosσ-K(γ-γref)

(36)

式中：K為反饋系數(shù)；γref為參考剖面對(duì)應(yīng)的航跡傾角,可通過將最后一次擬合的速度-高度表達(dá)式代入式(19)計(jì)算。

2.3 飛行時(shí)間控制

(37)

剩余飛行時(shí)間為

(38)

(39)

T(ha_i,hb_i)=-tgo_i-

(40)

此即飛行時(shí)間誤差。而航程誤差記為S(ha_i,hb_i),可通過式(20)計(jì)算。至此可建立飛行時(shí)間約束下剖面參數(shù)估計(jì)狀態(tài)方程為

(41)

觀測(cè)方程

(42)

(43)

同樣,當(dāng)前飛行狀態(tài)接近終端狀態(tài)時(shí),式(43)可能出現(xiàn)奇異。打靶仿真表明,該式出現(xiàn)奇異時(shí)對(duì)應(yīng)的當(dāng)前速度相對(duì)于2.2節(jié)無飛行時(shí)間約束工況有所增大。當(dāng)前飛行速度與終端速度之差小于470 m/s時(shí),停止預(yù)測(cè)校正制導(dǎo),不再更新飛行剖面,改為跟蹤最后一次產(chǎn)生的飛行剖面。

2.4 橫向制導(dǎo)

橫向制導(dǎo)采用經(jīng)典的基于航跡偏角偏差走廊的傾側(cè)角反轉(zhuǎn)制導(dǎo)律,表達(dá)如下:

(44)

式中：Δψu(yù)p和Δψdown表示航跡偏角走廊上邊界和下邊界；Δψ表示航跡偏角偏差,計(jì)算如下:

Δψ=ψ-ψLOS

(45)

ψLOS即當(dāng)前位置至目標(biāo)點(diǎn)的理想視線角,計(jì)算如下:

(46)

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

通過式(33)可知,制導(dǎo)指令主要受當(dāng)前狀態(tài)與末端狀態(tài)及升力系數(shù)、阻力系數(shù)等不確定性參數(shù)影響,對(duì)不確定性參數(shù)進(jìn)行在線辨識(shí),可提高制導(dǎo)精度。對(duì)制導(dǎo)指令及飛行狀態(tài)影響較大的待辨識(shí)參數(shù)主要包括升力系數(shù)偏差ΔCL、阻力系數(shù)偏差ΔCD、大氣密度偏差Δρ和攻角偏差Δα,這里假設(shè)攻角偏差主要由水平方向的風(fēng)速產(chǎn)生。

將偏差導(dǎo)數(shù)視為白噪聲形成狀態(tài)方程,ζ1,…,ζ4是互不相關(guān)高斯白噪聲,觀測(cè)量包括軸法向過載、動(dòng)壓。狀態(tài)方程如下：

(47)

式中:ΔCL_i、ΔCD_i、Δρi、Δαi為待辨識(shí)偏差量,ΔCL_i-1、ΔCD_i-1、Δρi-1、Δαi-1為前一步辨識(shí)結(jié)果。

觀測(cè)方程為

(48)

式中：α,ρ,V,γ為當(dāng)前標(biāo)準(zhǔn)參數(shù);nx、ny、q為觀測(cè)得到的軸向過載、法向過載和動(dòng)壓;L*,D*是受到擾動(dòng)后的升力和阻力,分別計(jì)算如下:

L*=(ρ+Δρ)(VcV)2Sref(CL+ΔCL)/(2mg0)

(49)

D*=(ρ+Δρ)(VcV)2Sref(CD+ΔCD)/(2mg0)

(50)

通過擴(kuò)展卡爾曼濾波(EKF),可以實(shí)現(xiàn)對(duì)不確定性參數(shù)的精確辨識(shí)。流程如圖1所示。

圖1 再入制導(dǎo)流程

3 仿真分析

本文基于X-33模型進(jìn)行仿真驗(yàn)證,該飛行器質(zhì)量37 363 kg,參考面積149.4 m2,所能承受的最大熱流為600 kW/m2,最大動(dòng)壓為160 kPa,最大過載為4g。再入過程中,制導(dǎo)周期為無量綱速度每變化0.001,即速度每變化7.91 m/s產(chǎn)生一次新的制導(dǎo)指令。攻角剖面如下:

(51)

本文首先進(jìn)行標(biāo)準(zhǔn)飛行狀態(tài)下的仿真,對(duì)4種不同航程下的算例進(jìn)行仿真驗(yàn)證；其次,針對(duì)同一算例,在不同飛行時(shí)間下開展仿真分析；最后,對(duì)初始狀態(tài)及升力系數(shù)、阻力系數(shù)、大氣密度、攻角等參數(shù)進(jìn)行拉偏,并線辨識(shí)不確定性參數(shù)偏差,通過蒙特卡洛仿真驗(yàn)證算法有效性。

3.1 無再入時(shí)間約束的標(biāo)稱狀態(tài)仿真

首先對(duì)8種不同航程下的算例進(jìn)行仿真驗(yàn)證,算例1～算例4再入高度和速度相同,主要驗(yàn)證標(biāo)稱狀態(tài)下不同航程下算法適應(yīng)能力及側(cè)向制導(dǎo)性能；而算例5～算例8分別對(duì)應(yīng)不同初始速度和高度條件下算法的適應(yīng)性。它們初始再入位置不同,但是交班點(diǎn)經(jīng)緯度均為(112°,42°),交班點(diǎn)高度為25 km,初始狀態(tài)如表1所示。再入過程中,對(duì)升力系數(shù)、阻力系數(shù)、大氣密度、攻角設(shè)置表2所示預(yù)置偏差,通過對(duì)偏差參數(shù)的辨識(shí),體現(xiàn)算法的有效性。

表1 初始條件

表2 參數(shù)偏差

圖2表示算例1～算例8所對(duì)應(yīng)的再入軌跡,包括高度-速度變化曲線、速度-待飛航程變化曲線、軸向和法向過載曲線及地面軌跡。通過高度-速度變化曲線可見這種算法對(duì)不同再入航程下的縱向和橫向制導(dǎo)均具有較強(qiáng)的適應(yīng)性,所得到的速度-高度曲線變化平滑,無明顯的躍起,滿足再入走廊要求。通過法向過載和軸向過載變化曲線可見,法向過載(虛線所示)小于4。由地面軌跡可見,飛行器均能準(zhǔn)確到達(dá)目標(biāo)點(diǎn),算法制導(dǎo)精度較高。

圖2 再入軌跡

圖3表示卡爾曼濾波估計(jì)的速度中點(diǎn)所對(duì)應(yīng)的無量綱飛行高度及傾側(cè)角指令,各組算例濾波在少數(shù)幾個(gè)周期的振蕩后,即進(jìn)入穩(wěn)定變化狀態(tài),具有收斂速度快的特點(diǎn)。

圖3 濾波估計(jì)結(jié)果及制導(dǎo)指令

圖4是不確定性參數(shù)在線辨識(shí)結(jié)果,可見算法對(duì)不確定性參數(shù)辨識(shí)收斂較快,能夠滿足對(duì)變化的不確定性參數(shù)的辨識(shí)需求。

圖4 不確定性參數(shù)辨識(shí)結(jié)果

表3表示本文算法與復(fù)現(xiàn)的文獻(xiàn)[8]所述數(shù)值預(yù)測(cè)校正算法(NPC)每個(gè)制導(dǎo)周期內(nèi)形成制導(dǎo)指令所需要的最大時(shí)長(zhǎng)和最小時(shí)長(zhǎng)對(duì)比。仿真工具為MATLAB(2013b),計(jì)算機(jī)操作系統(tǒng)為Windows 7,內(nèi)存4 G,主頻2.4 GHz,數(shù)值預(yù)測(cè)校正算法制導(dǎo)周期為1 s,預(yù)測(cè)步長(zhǎng)為5 s,每制導(dǎo)周期最大迭代次數(shù)為5次。通過對(duì)比可見本文算法計(jì)算效率具有一定的優(yōu)勢(shì),主要得益于僅需要在計(jì)算式(24)時(shí)需要2次單變量積分,避免了傳統(tǒng)預(yù)測(cè)校正制導(dǎo)算法傾側(cè)角指令生成過程中多次迭代積分縱向軌跡帶來的計(jì)算量增多問題。

表3 不同算例下制導(dǎo)指令生成所需時(shí)長(zhǎng)

3.2 再入時(shí)間約束下的標(biāo)稱狀態(tài)仿真

算例9～算例16主要對(duì)時(shí)間約束下飛行狀態(tài)進(jìn)行仿真研究。其中,算例9～算例12初始條件與3.1節(jié)算例2相同,分別對(duì)再入飛行時(shí)間630 s、660 s、690 s、720 s進(jìn)行仿真分析,主要驗(yàn)證飛行時(shí)間控制及側(cè)向制導(dǎo)性能。算例13～算例16經(jīng)緯度及航跡傾角均與3.1節(jié)算例2相同,但再入速度、再入高度不同,再入飛行時(shí)間均控制為690 s,主要驗(yàn)證飛行時(shí)間約束下縱向制導(dǎo)性能。算例9～算例16初始速度、高度及預(yù)定飛行時(shí)間如表4所示。

表4 初始速度、高度及預(yù)定飛行時(shí)間

圖5表示算例9～算例16再入飛行軌跡。通過算例9～算例12可見,在不同飛行時(shí)間下,均能夠滿足再入走廊約束,飛行時(shí)間控制精準(zhǔn),制導(dǎo)精度高。對(duì)比過載變化曲線及速度-高度曲線,可以發(fā)現(xiàn)相同再入條件下,預(yù)定飛行時(shí)間越長(zhǎng),飛行初始段過載越大而飛行末段過載越小。通過算例13～算例16可見,對(duì)于不同初始速度和高度的再入,可實(shí)現(xiàn)將飛行時(shí)間控制為相同,且具有良好的制導(dǎo)精度,過載等約束均可滿足要求。

圖6是不同飛行時(shí)間對(duì)應(yīng)的卡爾曼濾波估計(jì)結(jié)果和制導(dǎo)指令,可見考慮飛行時(shí)間約束后,濾波算法估計(jì)結(jié)果仍可迅速收斂。末端引入軌跡跟蹤制導(dǎo),可避免當(dāng)前速度接近終端速度時(shí)式(43)中矩陣求逆可能產(chǎn)生的奇異現(xiàn)象,確保傾側(cè)角指令變化平緩。

圖6 不同飛行時(shí)間下濾波估計(jì)結(jié)果及制導(dǎo)指令

表5表示8種不同再入條件下,每個(gè)制導(dǎo)周期內(nèi)制導(dǎo)指令形成所需要的最大及最小時(shí)長(zhǎng),仿真環(huán)境及平臺(tái)同3.1節(jié),每個(gè)制導(dǎo)周期0.069 s的最大制導(dǎo)指令生成時(shí)間仍能保證良好的實(shí)時(shí)性。

表5 再入時(shí)間約束下制導(dǎo)指令生成時(shí)長(zhǎng)

3.3 無再入時(shí)間約束的參數(shù)擾動(dòng)狀態(tài)仿真

首先基于算例2對(duì)應(yīng)的工況,考察濾波估計(jì)速度中點(diǎn)無量綱高度時(shí),初值選擇對(duì)估計(jì)結(jié)果的影響,將卡爾曼濾波估計(jì)速度中點(diǎn)對(duì)應(yīng)的無量綱高度初值從0.006 5變化至0.011 3(對(duì)應(yīng)于高度濾波初值從42 km變化至72 km)。

此時(shí)辨識(shí)效果如圖7所示,可見經(jīng)過初期少數(shù)幾個(gè)制導(dǎo)周期(<5)的調(diào)整后,濾波器將能夠平穩(wěn)地估計(jì)出速度中點(diǎn)對(duì)應(yīng)的無量綱高度,對(duì)初值選擇不敏感。

圖7 速度中點(diǎn)對(duì)應(yīng)無量綱高度濾波估計(jì)結(jié)果

考慮參數(shù)初始狀態(tài)偏差及不確定性參數(shù)辨識(shí)初值偏差,利用蒙特卡洛仿真驗(yàn)證算法魯棒性,仿真次數(shù)設(shè)定為300次。初值偏差均取為隨機(jī)偏差(表6),升力系數(shù)、阻力系數(shù)、大氣密度、攻角偏差取為隨機(jī)偏差和固定偏差的組合(表7)。

表6 初值參數(shù)偏差

表7 過程參數(shù)偏差

圖8表示蒙特卡洛仿真得到的飛行軌跡,可見算法對(duì)參數(shù)不確定性分布具有良好的適應(yīng)性,高度隨速度變化仍然較為平滑,且滿足再入走廊約束,制導(dǎo)精度高。

圖8 飛行軌跡蒙特卡洛仿真曲線

圖9表示速度中點(diǎn)對(duì)應(yīng)無量綱高度的估計(jì)結(jié)果及制導(dǎo)指令。濾波結(jié)果對(duì)于參數(shù)擾動(dòng)具有良好的適應(yīng)性,濾波值變化平穩(wěn),無跳變。

圖9 速度中點(diǎn)對(duì)應(yīng)無量綱高度及制導(dǎo)指令蒙特卡洛仿真

圖10表示交班點(diǎn)經(jīng)緯度分布,均分布在理想交班點(diǎn)附近。

圖10 蒙特卡洛仿真交班點(diǎn)分布

參數(shù)辨識(shí)結(jié)果如圖11所示,偏差比定義為

圖11 蒙特卡洛仿真參數(shù)偏差辨識(shí)結(jié)果

(52)

3.4 考慮再入時(shí)間約束的參數(shù)擾動(dòng)狀態(tài)仿真

基于算例11對(duì)應(yīng)的工況進(jìn)行蒙特卡洛仿真。首先考察2個(gè)待估計(jì)高度初值選擇對(duì)估計(jì)結(jié)果的影響。圖12表示2個(gè)高度初值分別在0.007 1～

圖12 初值擾動(dòng)下的特征高度估計(jì)結(jié)果

0.011 8、0.004 7～0.009 4范圍內(nèi)隨機(jī)取值,仿真20組的濾波結(jié)果。經(jīng)過5個(gè)制導(dǎo)周期,初值擾動(dòng)可收斂。

參數(shù)偏差設(shè)置如表6和表7所示,圖13表示飛行軌跡,可見算法仍具有良好的適應(yīng)性,高度及速度隨飛行時(shí)間變化較為平滑,時(shí)間控制精準(zhǔn),落點(diǎn)精確,過載相對(duì)于無飛行時(shí)間約束情況有所增大,但仍滿足再入走廊約束。

圖13 再入時(shí)間約束下飛行軌跡變化蒙特卡洛仿真

圖14表示參數(shù)擾動(dòng)條件下2個(gè)待估計(jì)高度變化曲線及制導(dǎo)指令變化曲線。圖15表示交班點(diǎn)分布及時(shí)間誤差分布?？梢娫诒疚目紤]時(shí)間約束再入制導(dǎo)算法下,制導(dǎo)精度較高且飛行時(shí)間與預(yù)定飛行時(shí)間之差可控制在2.5 s以內(nèi)。通過圖16可見,不確定性參數(shù)辨識(shí)結(jié)果與3.3節(jié)非常類似,具有良好的辨識(shí)精度。

圖14 無量綱高度濾波估計(jì)和制導(dǎo)指令蒙特卡洛仿真

圖15 交班點(diǎn)及時(shí)間偏差分布蒙特卡洛仿真

圖16 不確定性參數(shù)辨識(shí)蒙特卡洛仿真

4 結(jié) 論

本文研究了利用卡爾曼濾波快速生成滿足航程要求的飛行剖面,并基于當(dāng)前狀態(tài)和濾波結(jié)果,

通過解析方式產(chǎn)生傾側(cè)角指令的預(yù)測(cè)校正制導(dǎo)問題。仿真結(jié)果表明：

1) 通過將高度擬合為速度的四次多項(xiàng)式,利用濾波估計(jì)指定點(diǎn)對(duì)應(yīng)的高度,即可快速生成滿足要求的飛行剖面,有效避免了傳統(tǒng)預(yù)測(cè)校正制導(dǎo)中需要多次迭代積分縱向彈道帶來的計(jì)算較為復(fù)雜及再入飛行時(shí)間難以控制的問題。

2) 通過增加待估計(jì)參數(shù),實(shí)現(xiàn)將飛行器導(dǎo)引至目標(biāo)點(diǎn)的同時(shí),可實(shí)現(xiàn)對(duì)再入飛行時(shí)間的精確控制。

3) 針對(duì)末端擬合精度較差的問題,通過設(shè)置更新終止條件并跟蹤最后一次規(guī)劃產(chǎn)生的標(biāo)準(zhǔn)飛行剖面,同時(shí)引入航跡傾角反饋,可有效提升了算法在飛行末端的收斂性。

4) 在飛行過程中,通過對(duì)升力系數(shù)、阻力系數(shù)、攻角、大氣密度偏差進(jìn)行在線辨識(shí),有效補(bǔ)償偏差,提高制導(dǎo)精度。