基于Kalman濾波的空天飛行器再入制導算法

2021-11-30 14:35:26尤志鵬楊勇劉剛曹曉瑞鄭宏濤

航空學報 2021年11期

關鍵詞：指令

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

中國運載火箭技術研究院,北京 100076

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

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

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

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

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

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

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

1 運動方程及約束

1.1 運動方程

假設地球是均質圓球,三維質點再入運動無量綱方程為

(1)

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

(2)

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

(3)

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

1.2 飛行約束

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

(4)

(5)

n=|Lcosα+Dsinα|≤nmax

(6)

(7)

(8)

(9)

(10)

(11)

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

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

h(Vf)=hf

(12)

θ(Vf)=θf

(13)

φ(Vf)=φf

(14)

式中:hf、θf和φf分別表示再入交班點高度、末端經度和末端緯度。

2 制導算法設計

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

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

本文在速度-高度剖面內設計再入制導律,將飛行高度擬合為速度的多項式形式。在第i個制導周期,已知狀態主要包括當前點的速度Vi、高度hi、航跡傾角γi,交班點的速度Vf、高度hf；在高度-速度剖面內,有

(15)

h=a4V4+a3V3+a2V2+a1V+a0

(16)

式中:a0～a4為系數。

在高度-速度剖面內,待飛航程隨速度變化表達式為

(17)

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

(18)

式中:R為待飛航程,即當前點至目標點對應的航程。

根據式(16),有

(19)

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

假設該點對應的高度為hm_i,記

S(hm_i)=-R+

(20)

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

2.2 縱向制導指令

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

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

建模過程如下：

狀態方程：

hm_i+1=hm_i+εi

(21)

觀測方程：

Zi=S(hm_i)+νi

(22)

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

濾波過程如下：

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

(23)

(24)

式中:δ為一小攝動量,取為10-5。

(25)

則

(26)

而

(27)

進而

(28)

式中:ηj,ωj分別是Gauss-Legendre節點和節點權值系數,參考文獻[24],通常六階Gauss-Legendre求積公式即可得到相當高的精度。

當hm_i+1確定后,即可通過式(29)求解該制導周期對應的速度-高度剖面擬合系數

(29)

為獲得傾側角指令,可將飛行高度對速度求二階導數,得

(30)

其中

(31)

(32)

(33)

這里下標i表示第i個制導周期,其中

(34)

當前飛行狀態接近終端狀態時,打靶仿真發現,式(29)可能出現奇異,導致擬合不準確。為解決該問題,當前飛行速度與終端速度之差小于100 m/s時,停止預測校正制導,不再更新飛行剖面,改為跟蹤最后一次產生的飛行剖面,即

(35)

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

(36)

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

2.3 飛行時間控制

(37)

剩余飛行時間為

(38)

(39)

T(ha_i,hb_i)=-tgo_i-

(40)

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

(41)

觀測方程

(42)

(43)

同樣,當前飛行狀態接近終端狀態時,式(43)可能出現奇異。打靶仿真表明,該式出現奇異時對應的當前速度相對于2.2節無飛行時間約束工況有所增大。當前飛行速度與終端速度之差小于470 m/s時,停止預測校正制導,不再更新飛行剖面,改為跟蹤最后一次產生的飛行剖面。

2.4 橫向制導

橫向制導采用經典的基于航跡偏角偏差走廊的傾側角反轉制導律,表達如下:

(44)

式中：Δψup和Δψdown表示航跡偏角走廊上邊界和下邊界；Δψ表示航跡偏角偏差,計算如下:

Δψ=ψ-ψLOS

(45)

ψLOS即當前位置至目標點的理想視線角,計算如下:

(46)

2.5 不確定性參數在線辨識

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

將偏差導數視為白噪聲形成狀態方程,ζ1,…,ζ4是互不相關高斯白噪聲,觀測量包括軸法向過載、動壓。狀態方程如下：

(47)

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

觀測方程為

(48)

式中：α,ρ,V,γ為當前標準參數;nx、ny、q為觀測得到的軸向過載、法向過載和動壓;L*,D*是受到擾動后的升力和阻力,分別計算如下:

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

(49)

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

(50)

通過擴展卡爾曼濾波(EKF),可以實現對不確定性參數的精確辨識。流程如圖1所示。

圖1 再入制導流程

3 仿真分析

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

(51)

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

3.1 無再入時間約束的標稱狀態仿真

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

表1 初始條件

表2 參數偏差

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

圖2 再入軌跡

圖3表示卡爾曼濾波估計的速度中點所對應的無量綱飛行高度及傾側角指令,各組算例濾波在少數幾個周期的振蕩后,即進入穩定變化狀態,具有收斂速度快的特點。

圖3 濾波估計結果及制導指令

圖4是不確定性參數在線辨識結果,可見算法對不確定性參數辨識收斂較快,能夠滿足對變化的不確定性參數的辨識需求。

圖4 不確定性參數辨識結果

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

表3 不同算例下制導指令生成所需時長

3.2 再入時間約束下的標稱狀態仿真

算例9～算例16主要對時間約束下飛行狀態進行仿真研究。其中,算例9～算例12初始條件與3.1節算例2相同,分別對再入飛行時間630 s、660 s、690 s、720 s進行仿真分析,主要驗證飛行時間控制及側向制導性能。算例13～算例16經緯度及航跡傾角均與3.1節算例2相同,但再入速度、再入高度不同,再入飛行時間均控制為690 s,主要驗證飛行時間約束下縱向制導性能。算例9～算例16初始速度、高度及預定飛行時間如表4所示。

表4 初始速度、高度及預定飛行時間

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

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