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

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

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

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

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

在彈道預(yù)報方面,研究的關(guān)鍵主要有兩方面:一是獲取準(zhǔn)確的預(yù)報初值;二是建立精確的動力學(xué)模型或運動學(xué)模型。文獻(xiàn)[15,16]采用彈道參數(shù)辨識和狀態(tài)估計的思想獲取預(yù)報初值,通過建立修正動力學(xué)模型,實現(xiàn)對彈道的預(yù)報。文獻(xiàn)[17]通過構(gòu)建深度神經(jīng)網(wǎng)絡(luò),基于傳統(tǒng)模型彈道數(shù)據(jù)進(jìn)行訓(xùn)練,實現(xiàn)了導(dǎo)彈彈道的規(guī)劃和快速預(yù)報,但這種方法對強(qiáng)機(jī)動目標(biāo)效果較差。

在已有文獻(xiàn)研究的基礎(chǔ)上,本文采用“濾波估計—彈道預(yù)報”的思想解決紅外探測條件下的攔截彈制導(dǎo)律辨識與彈道預(yù)報問題。首先,根據(jù)矩陣?yán)碚搶ο到y(tǒng)可觀性進(jìn)行證明分析。在此基礎(chǔ)上,將制導(dǎo)律系數(shù)增廣為濾波狀態(tài)分量,構(gòu)造濾波器對制導(dǎo)律系數(shù)與運動狀態(tài)進(jìn)行辨識估計,并給出一種濾波穩(wěn)定判斷方法,將穩(wěn)定收斂時刻的制導(dǎo)律系數(shù)及運動狀態(tài)作為彈道預(yù)報初值。最后,構(gòu)建攔截彈彈道預(yù)報方程,采用無跡變換的方法進(jìn)行彈道預(yù)報與誤差傳播分析。

1 制導(dǎo)律辨識的可觀性分析

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

攔截彈的狀態(tài)矢量為:

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

XMk=XIk+ΔRk

(1)

假設(shè)加速度符合指數(shù)模型,此時飛行器與攔截彈的狀態(tài)矢量分別為:

(2)

(3)

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

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

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

(4)

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

(5)

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

(6)

式中:

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

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

BX=C

(7)

式中:

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

①當(dāng)視線方位角β為固定值時,矩陣B經(jīng)過初等變換,有:

B→(B0Om×2)

(8)

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

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

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

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

2 制導(dǎo)律系數(shù)辨識模型

2.1 系統(tǒng)狀態(tài)方程

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

(9)

(10)

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

(11)

式中:ζ為發(fā)射方位角,φ為發(fā)射點的緯度,λ為發(fā)射點的經(jīng)度。

(12)

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

(13)

對式(13)進(jìn)行離散化處理,可以得到制導(dǎo)段系統(tǒng)狀態(tài)模型。

Xk+1=ΦkXk+Ukak

(14)

式中:

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

2.2 系統(tǒng)觀測方程

當(dāng)飛行器利用紅外探測器對攔截彈進(jìn)行探測跟蹤時,測量數(shù)據(jù)為觀測方位角ζ、觀測俯仰角ξ,與攔截彈狀態(tài)矢量X間的關(guān)系為:

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

(15)

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

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

(16)

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

(17)

式(16)～式(17)構(gòu)成了觀測方程的解析表達(dá)式。

3 制導(dǎo)律系數(shù)辨識濾波器設(shè)計

3.1 過程噪聲方差設(shè)計

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

(18)

式中:τ為機(jī)動時間常數(shù)。

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

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

(19)

式中:φk為用于辨識制導(dǎo)律系數(shù)的離散化驅(qū)動白噪聲,噪聲方差為:

(20)

(21)

對于攔截彈的制導(dǎo)律系數(shù),其取值一般為有限區(qū)間Ni∈[Nmin,Nmax],假設(shè)制導(dǎo)律系數(shù)在區(qū)間內(nèi)滿足均勻分布,則其方差為:

(22)

在制導(dǎo)段攔截彈一直處于機(jī)動狀態(tài),機(jī)動頻率取飛行時間的倒數(shù)。

3.2 濾波算法流程

①使用式(9)計算ak/k,并進(jìn)行狀態(tài)預(yù)測。

(23)

②使用3.1節(jié)方法計算過程噪聲方差Qk,用于補(bǔ)償預(yù)測協(xié)方差矩陣。

(24)

③采用修正比例采樣規(guī)則構(gòu)造Sigma采樣點χk+1/k(l)和相應(yīng)權(quán)重Wl。

(25)

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

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

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

(26)

⑤計算觀測值的預(yù)估值和協(xié)方差矩陣。

(27)

⑥計算新息矢量和新息協(xié)方差矩陣。

(28)

⑦狀態(tài)估計與協(xié)方差更新。

(29)

4 彈道預(yù)報與誤差傳播分析

4.1 預(yù)報初值確定

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

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

(30)

4.2 彈道預(yù)報方程

(31)

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

(32)

根據(jù)濾波穩(wěn)定收斂的結(jié)果,采用四階龍格庫塔積分求解微分方程得到彈道預(yù)報方程。

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

(33)

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

4.3 誤差傳播分析

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

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

(34)

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

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

(35)

(36)

具體流程如圖2所示。

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

5 仿真實驗與結(jié)果分析

5.1 仿真場景

根據(jù)飛行器動力學(xué)方程和數(shù)值積分算法,優(yōu)化一條飛行器滑翔彈道。攔截彈發(fā)射點地理坐標(biāo)為19.6°E、0.4°N,高程為5 m,發(fā)射方位角為-117.65°,在發(fā)射后26.2 s進(jìn)入制導(dǎo)段,采用比例導(dǎo)引方式,制導(dǎo)律系數(shù)為4。以飛行器飛行時間為基準(zhǔn),攔截彈在575.4 s發(fā)射,在601.6 s進(jìn)入制導(dǎo)段,構(gòu)建三維場景如圖3所示。

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

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

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

可以看出,視線角有較大程度的變化,根據(jù)第2節(jié)分析可知,系統(tǒng)是可觀的。

5.2 制導(dǎo)律辨識

①參數(shù)設(shè)置

假設(shè)以攔截彈助推段濾波估計結(jié)束后的數(shù)據(jù)作為攔截彈初始狀態(tài)矢量,制導(dǎo)律系數(shù)取值范圍N∈[2,6],根據(jù)經(jīng)驗取值為3,則濾波初始估計值為:

初始協(xié)方差矩陣可設(shè)為:

②結(jié)果分析

根據(jù)上面的參數(shù)設(shè)置,本文算法的制導(dǎo)律系數(shù)的辨識結(jié)果如圖5所示,對比方法的制導(dǎo)律辨識結(jié)果如圖6所示。

圖5 攔截彈制導(dǎo)律系數(shù)估計值Fig.5 Estimation of guidance law coefficient for interceptor

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

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

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

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

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

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

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

5.3 彈道預(yù)報

①參數(shù)設(shè)置

根據(jù)濾波估計得到的攔截彈制導(dǎo)律系數(shù)及狀態(tài)估計值,采用第4節(jié)中方法進(jìn)行彈道預(yù)報。濾波數(shù)據(jù)處理過程中,w=20,閾值κ=5。經(jīng)計算,在607.8 s時攔截彈運動狀態(tài)濾波結(jié)果整體穩(wěn)定收斂,取該時刻數(shù)據(jù)作為預(yù)報初值,此時的狀態(tài)估計值為:

對應(yīng)的估計協(xié)方差矩陣為:

對應(yīng)的3σ誤差橢球如圖10所示。

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

②結(jié)果分析

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

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

彈道預(yù)報位置的標(biāo)準(zhǔn)差與實際預(yù)報誤差的變化曲線如圖12所示,其中虛線表示彈道預(yù)報的標(biāo)準(zhǔn)差,實線表示彈道預(yù)報的實際誤差。

圖12 彈道預(yù)報標(biāo)準(zhǔn)差與實際誤差曲線Fig.12 Trajectory prediction standard deviation and actual error

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

通過分析可知,引入無跡變換的彈道預(yù)報算法能夠?qū)崿F(xiàn)攔截彈彈道的準(zhǔn)確預(yù)報,并且能夠給出彈道預(yù)報的誤差范圍,但彈道預(yù)報的精度主要取決于對攔截彈制導(dǎo)律系數(shù)與運動狀態(tài)的估計精度。

6 結(jié)論

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