紅外探測下的攔截彈制導律辨識與彈道預報方法

張雪松,王 鋒,李東澤,朱雅萌,儲思思

(1.解放軍信息工程大學 數據與目標工程學院,河南 鄭州 450000;2.中國人民解放軍32035部隊,陜西 西安 710000;3.軍事科學院 國防科技創新研究院,北京 100071;4.中國運載火箭技術研究院,北京 100071)

利用飛行器搭載的紅外傳感器對攔截彈進行探測,不僅豐富了探測手段、拓展了探測范圍,同時也保證了平臺的隱蔽性、增強了生存能力,對于機動突防和反攔截有著十分重要的意義[1-3]。而飛行器進行機動規避實施突防的關鍵是對攔截彈制導律的辨識、運動狀態的估計以及后續彈道的預報。

由于紅外探測器測量時僅能獲得角度信息,對攔截彈定位屬于不完備觀測,而系統的可觀性是后續制導律辨識以及運動狀態估計的前提。目前,對系統可觀性問題的研究較多,常用的方法有微分幾何方法、數值分析方法等[4,5]。文獻[6,7]采用微分幾何方法對制導律辨識的非線性系統進行了可觀性分析,給出了系統弱可觀的充要條件。文獻[8]構建了被動傳感器對再入階段彈道目標的運動參數估計模型,并基于Fisher信息矩陣對系統的可觀性進行了數值驗證。

根據零化視線角速率的思想,各類制導律都可以等價為不同制導律系數的比例導引[9,10],因此對制導律的辨識可以轉化為對比例導引模型中制導律系數的辨識估計。文獻[11-13]基于多模型自適應卡爾曼濾波算法,將制導律辨識轉化為對模型的辨識,能夠準確地辨識出制導律類型和制導參數,但當模型集內無匹配的制導律時,則辨識效果較差。文獻[14]以飛行器和攔截彈的相對運動數據為樣本輸入,攔截彈制導律為標簽,基于門循環單元(GRU)實現了制導律的準確辨識,但該方法未考慮雙方的動力學信息,辨識的準確性有待進一步研究。

在彈道預報方面,研究的關鍵主要有兩方面:一是獲取準確的預報初值;二是建立精確的動力學模型或運動學模型。文獻[15,16]采用彈道參數辨識和狀態估計的思想獲取預報初值,通過建立修正動力學模型,實現對彈道的預報。文獻[17]通過構建深度神經網絡,基于傳統模型彈道數據進行訓練,實現了導彈彈道的規劃和快速預報,但這種方法對強機動目標效果較差。

在已有文獻研究的基礎上,本文采用“濾波估計—彈道預報”的思想解決紅外探測條件下的攔截彈制導律辨識與彈道預報問題。首先,根據矩陣理論對系統可觀性進行證明分析。在此基礎上,將制導律系數增廣為濾波狀態分量,構造濾波器對制導律系數與運動狀態進行辨識估計,并給出一種濾波穩定判斷方法,將穩定收斂時刻的制導律系數及運動狀態作為彈道預報初值。最后,構建攔截彈彈道預報方程,采用無跡變換的方法進行彈道預報與誤差傳播分析。

1 制導律辨識的可觀性分析

為了便于理論分析,這里不考慮噪聲的影響。設初始時刻飛行器的狀態矢量為:

攔截彈的狀態矢量為:

在k時刻,相對距離矢量為ΔRk,則飛行器與攔截彈的位置矢量關系為:

XMk=XIk+ΔRk

(1)

假設加速度符合指數模型,此時飛行器與攔截彈的狀態矢量分別為:

(2)

(3)

式中:λi,ωi(i=x,y,z)為機動時間常數。

飛行器對攔截彈的視線方位角為βk,視線高低角為εk,則沿視線坐標系y軸的單位矢量為:

uk=(-cosβksinεkcosεksinβksinεk)T

(4)

根據矢量的正交性,由式(1)、式(4)可得:

(5)

將式(2)～式(4)代入,經過整理可以得到:

(6)

式中:

可觀性即在相應的觀測時段能夠唯一確定系統的初值。對于式(6),等式兩邊左乘AT,通過觀察分析可知,ATA為6×6對稱矩陣,且前3行與后3行對應成比例,即Rank(ATA)=3,可知ATA不滿秩。因此,對于k時刻的一次觀測,式(6)沒有唯一解,系統不可觀。

當飛行器對攔截彈連續觀測,即k=1,2…m,式(6)可以寫成:

BX=C

(7)

式中:

對于觀測矩陣B,考慮以下幾種情況:

①當視線方位角β為固定值時,矩陣B經過初等變換,有:

B→(B0Om×2)

(8)

式中:B0為m×4矩陣,Om×2為m×2零矩陣。

根據矩陣秩的定義可知,矩陣B的非零子式的最高階小于6,即Rank(B)<6。而式(7)的狀態變量維數為6,因此沒有唯一解,系統不可觀。

同理,對于視線高低角為固定值時,系統同樣不可觀。

②當視線方位角、視線高低角有變化時,對于矩陣B進行分析可知,矩陣B為列滿秩,即列向量線性無關,也就是說B的零空間僅包含零向量。對于BTB來說,構造零空間BTBv=0,則有vTBTBv=0,即(Bv)TBv=0,所以Bv=0。由于v是B的零空間向量,所以BTB也是列滿秩,而BTB為6×6對稱矩陣,行秩等于列秩,即Rank(BTB)=6,滿足唯一解條件,即系統可觀。

2 制導律系數辨識模型

2.1 系統狀態方程

在攔截彈的制導段,主要受到地球引力、氣動力、離心力、科氏力和制導力的作用?？偧铀俣萢由以下幾項組成,包括引力加速度ag、制導力加速度as,氣動力加速度aA、離心力加速度acf和科氏力加速度aco,即a=ag+as+aA+acf+aco。

(9)

(10)

式中:φ、ψ、γ分別為俯仰角、偏航角和滾轉角,P1、P2、P3分別為繞x、y、z旋轉一定角度的方向余弦矩陣。

(11)

式中:ζ為發射方位角,φ為發射點的緯度,λ為發射點的經度。

(12)

式中:ν、σ、θ分別為傾側角、速度偏角、速度傾角。

(13)

對式(13)進行離散化處理,可以得到制導段系統狀態模型。

Xk+1=ΦkXk+Ukak

(14)

式中:

其中I為單位矩陣,O為零矩陣,T為采樣周期。

2.2 系統觀測方程

當飛行器利用紅外探測器對攔截彈進行探測跟蹤時,測量數據為觀測方位角ζ、觀測俯仰角ξ,與攔截彈狀態矢量X間的關系為:

Z(ζ,ξ)=H(X)+V

(15)

式中:V為等效測量噪聲,其協方差矩陣為Cov=E(VVT),H表示攔截彈狀態矢量X與觀測角度的變換關系。

在地固坐標系中,設飛行器的位置矢量為(xMyMzM)T,攔截彈的位置矢量為(xIyIzI)T,則攔截彈在飛行器視場坐標系中的位置矢量為:

(16)

由此可以計算紅外探測器的測量值為:

(17)

式(16)～式(17)構成了觀測方程的解析表達式。

3 制導律系數辨識濾波器設計

3.1 過程噪聲方差設計

濾波器的性能主要取決于過程噪聲方差的設計,對于攔截彈的制導律系數來說,其作為制導控制指令參數通常具有非隨機性和相關性,因此可以采用自相關函數為指數衰減形式的一階Markov過程進行表示[18]。

(18)

式中:τ為機動時間常數。

將式(18)進行離散化,可以得到:

(Ni)k+1=(Ni)kexp(-T/τ)+φk

(19)

式中:φk為用于辨識制導律系數的離散化驅動白噪聲,噪聲方差為:

(20)

(21)

對于攔截彈的制導律系數,其取值一般為有限區間Ni∈[Nmin,Nmax],假設制導律系數在區間內滿足均勻分布,則其方差為:

(22)

在制導段攔截彈一直處于機動狀態,機動頻率取飛行時間的倒數。

3.2 濾波算法流程

①使用式(9)計算ak/k,并進行狀態預測。

(23)

②使用3.1節方法計算過程噪聲方差Qk,用于補償預測協方差矩陣。

(24)

③采用修正比例采樣規則構造Sigma采樣點χk+1/k(l)和相應權重Wl。

(25)

式中:n為狀態向量維度,κ表示Sigma采樣點的散布程度,一般取n+κ=3。

④基于觀測方程的非線性傳播:

Zk+1/k(l)=H[χk+1/k(l)]l=0,1,…,2n

(26)

⑤計算觀測值的預估值和協方差矩陣。

(27)

⑥計算新息矢量和新息協方差矩陣。

(28)

⑦狀態估計與協方差更新。

(29)

4 彈道預報與誤差傳播分析

4.1 預報初值確定

對攔截彈的彈道預報需要濾波結果達到穩定收斂后才能進行。對于濾波得到的數據X1,X2,X3…,從第j個數據Xj開始,依次對w個數據相鄰求差的值取模

|Xj-Xj-1|=δ1
|Xj-1-Xj-2|=δ2
?
|Xj-w+2-Xj-w+1|=δw-1

(30)

4.2 彈道預報方程

(31)

攔截彈在制導段總的加速度可由式(9)計算。同時,考慮攔截彈具有過載飽和,可以表示為:

(32)

根據濾波穩定收斂的結果,采用四階龍格庫塔積分求解微分方程得到彈道預報方程。

y=RK4[f(X),Sp,y0,h]

(33)

式中:Sp為積分區間,y0為積分初值向量,h為積分步長。

4.3 誤差傳播分析

由于預報初值是濾波穩定收斂時的估計值及對應的協方差,在利用有誤差的估計值進行彈道預報時,會得到的一個逐漸擴散的“誤差管道”,如圖1所示。因此,需要對攔截彈彈道報結果進行誤差傳播分析。常用的分析方法有協方差描述函數法和無跡變換法,而無跡變換不需要對非線性方程線性化,精度更高[19,20],因此本文采用基于無跡變換的方法進行誤差傳播分析。

圖1 飛行器對攔截彈探測示意圖Fig.1 Diagram of aircraft detecting interceptor

(34)

②基于式(33)的攔截彈彈道預報方程,得到sigma采樣點的輸出值。

γl=RK4(χl)l=0,1,…,2n

(35)

(36)

具體流程如圖2所示。

圖2 基于無跡變換的彈道預報流程Fig.2 Trajectory prediction process based on Unscented Transformation

5 仿真實驗與結果分析

5.1 仿真場景

根據飛行器動力學方程和數值積分算法,優化一條飛行器滑翔彈道。攔截彈發射點地理坐標為19.6°E、0.4°N,高程為5 m,發射方位角為-117.65°,在發射后26.2 s進入制導段,采用比例導引方式,制導律系數為4。以飛行器飛行時間為基準,攔截彈在575.4 s發射,在601.6 s進入制導段,構建三維場景如圖3所示。

圖3 三維彈道曲線Fig.3 Three-dimensional trajectory

飛行器搭載紅外探測器對攔截彈進行探測跟蹤,采樣周期為T=0.1 s,視線方位角、視線高低角如圖4所示。

圖4 視線角變化曲線Fig.4 Variation curve of line-of-sight

可以看出,視線角有較大程度的變化,根據第2節分析可知,系統是可觀的。

5.2 制導律辨識

①參數設置

假設以攔截彈助推段濾波估計結束后的數據作為攔截彈初始狀態矢量,制導律系數取值范圍N∈[2,6],根據經驗取值為3,則濾波初始估計值為:

初始協方差矩陣可設為:

②結果分析

根據上面的參數設置,本文算法的制導律系數的辨識結果如圖5所示,對比方法的制導律辨識結果如圖6所示。

圖5 攔截彈制導律系數估計值Fig.5 Estimation of guidance law coefficient for interceptor

圖6 對比方法的制導律辨識結果(有匹配制導律)Fig.6 The guidance law identification result of the comparison method(with matched guidance law)

從圖5辨識結果可以看出,約在607 s制導律系數開始穩定收斂于真值附近,辨識結果較好,而圖6中對比方法約在610 s左右辨識出制導律模型。本文算法能夠直接對制導律系數實時估計,不需要構建制導律系數模型集,且辨識速度較快。而當模型集中無正確的制導律系數模型時,即N∈{2,2.5,3,3.5,4.5,5,5.5,6},辨識結果如圖7所示,其辨識結果是錯誤的,可見該算法的辨識準確性受限于模型的數量。而本文算法解決了對比算法因模型集有限而導致的辨識估計不準確的問題,具有更強的實時性和適用性。

圖7 對比方法的制導律辨識結果(無匹配制導律)Fig.7 The guidance law identification result of the comparison method(without matched guidance law)

對于攔截彈運動狀態估計,采用均方根誤差RMSE表示跟蹤過程中的誤差變化。以x軸方向狀態估計均方根誤差為例,其它方向基本相同,進行100次蒙特卡羅仿真,仿真結果如圖8、圖9所示。

圖8 x方向位置均方根誤差Fig.8 RMSE of position for x-direction

圖9 x方向速度均方根誤差Fig.9 RMSE of velocity for x-direction

實線描述本文算法的估計結果,虛線描述對比方法的估計結果。從圖中可以看出,本文算法濾波穩定時的位置均方根誤差約為95 m,而對比算法的位置均方根誤差呈緩慢發散趨勢;本文算法濾波穩定時的速度均方根誤差約為3 m/s,而對比算法的速度均方根誤差在610 s前與本文算法基本一致,而后逐漸發散。由此可見,本文算法濾波穩定性更好,對攔截彈的濾波更為穩定,且在整體上濾波估計的精度更高。

5.3 彈道預報

①參數設置

根據濾波估計得到的攔截彈制導律系數及狀態估計值,采用第4節中方法進行彈道預報。濾波數據處理過程中,w=20,閾值κ=5。經計算,在607.8 s時攔截彈運動狀態濾波結果整體穩定收斂,取該時刻數據作為預報初值,此時的狀態估計值為:

對應的估計協方差矩陣為:

對應的3σ誤差橢球如圖10所示。

圖10 預報初始時刻誤差橢球Fig.10 Error ellipsoid of initial time

②結果分析

圖11 彈道預報誤差橢球變化曲線Fig.11 Error ellipsoid variation of trajectory prediction

由圖可知,攔截彈預報初始時刻的誤差橢球隨著預報時間的增加而逐漸擴散,并形成了“誤差管道”。而且在y方向上的預報誤差增速較大,在x方向、z方向上的增速較小。由于是3σ橢球,真實的彈道將以97.3%的概率從誤差橢球構成的“誤差管道”穿過。

彈道預報位置的標準差與實際預報誤差的變化曲線如圖12所示,其中虛線表示彈道預報的標準差,實線表示彈道預報的實際誤差。

圖12 彈道預報標準差與實際誤差曲線Fig.12 Trajectory prediction standard deviation and actual error

由圖可知:各方向彈道預報的實際誤差基本落在標準差范圍內,y方向彈道預報的實際誤差在起始階段低于誤差標準差下界,這是因為在“濾波估計—彈道預報”過程中,濾波穩定收斂時,在y方向攔截彈位置估計誤差大于其他兩個方向。在整個預報過程中,隨著預報時刻的推移,各方向實際誤差的變化趨勢與標準差相符,位置誤差均不斷累積,呈發散趨勢,但整體均在標準差范圍內。在預報終止時刻,x方向實際誤差為85 m,預報標準差為240.7 m;y方向實際誤差為292.6 m,預報標準差為537.5 m;z方向實際誤差為25 m,預報標準差為80 m,總體來說彈道預報效果較好。

通過分析可知,引入無跡變換的彈道預報算法能夠實現攔截彈彈道的準確預報,并且能夠給出彈道預報的誤差范圍,但彈道預報的精度主要取決于對攔截彈制導律系數與運動狀態的估計精度。

6 結論

本文針對紅外探測下的攔截彈制導率辨識與彈道預報問題,主要做了以下工作:一是分析了在僅有角度測量信息的條件下,系統的可觀性取決于飛行器與攔截彈的觀測角變化情況;二是將制導律系數增廣為濾波狀態分量,提出一種過程噪聲自適應的無跡卡爾曼濾波算法,實現了對攔截彈制導律系數與運動狀態的準確估計;三是采用無跡變換的方法對攔截彈的彈道進行預報和誤差傳播分析。仿真結果表明,本文提出的算法是有效的,對制導律系數和運動狀態的辨識與估計在收斂速度和精度上能夠滿足后續計算要求,同時對彈道預報的結果能夠為飛行器機動突防和規避攔截提供理論支撐。