彈道目標(biāo)的變結(jié)構(gòu)多模型數(shù)據(jù)濾波方法*

黃晶晶 陳世友

（武漢數(shù)字工程研究所 武漢 430074）

1 引言

隨著彈道導(dǎo)彈的發(fā)展，彈道導(dǎo)彈防御技術(shù)［1］的研究也愈加受到重視，彈道目標(biāo)數(shù)據(jù)濾波技術(shù)是彈道導(dǎo)彈防御的核心技術(shù)之一，對(duì)導(dǎo)彈落點(diǎn)預(yù)測(cè)的準(zhǔn)確度和彈道導(dǎo)彈攔截的成功率至關(guān)重要［2］。根據(jù)彈道導(dǎo)彈運(yùn)動(dòng)過(guò)程中的受力情況，可將其彈道分為主動(dòng)段、自由段和再入段三個(gè)飛行階段。由于各階段運(yùn)動(dòng)特征差別較大，單模型彈道目標(biāo)數(shù)據(jù)濾波方法不能實(shí)現(xiàn)對(duì)全彈道的連續(xù)估計(jì)。使用交互多模型算法可以解決上述問(wèn)題，但也存在一些缺點(diǎn)，交互多模型算法的估計(jì)性能與所選模型集密切相關(guān)［3～4］，而模型集的選擇通常處于兩難的境地，為提高估計(jì)精度需要盡可能多地增加模型的數(shù)量，而過(guò)多的模型會(huì)增加模型之間的競(jìng)爭(zhēng)并且增加計(jì)算量，導(dǎo)致估計(jì)精度降低。另外，由于彈道目標(biāo)各階段的運(yùn)動(dòng)特征明顯不同，每個(gè)階段均使用所有模型來(lái)濾波會(huì)增加不必要的計(jì)算量。

為克服上述缺點(diǎn)，本文基于變結(jié)構(gòu)多模型算法［5］，提出了一個(gè)適用于彈道目標(biāo)全階段的濾波方法。該方法采用多模型的結(jié)構(gòu)，相較于單模型跟蹤算法，可用于跟蹤任意階段的彈道目標(biāo)；相較于交互多模型算法，可以減少不必要的計(jì)算，同時(shí)提高精度。本文基于三個(gè)不同階段分別建立運(yùn)動(dòng)模型，每個(gè)階段對(duì)應(yīng)一個(gè)變結(jié)構(gòu)多模型的濾波器。其中，主動(dòng)段采用Singer模型［6］和重力轉(zhuǎn)彎模型［7］兩種模型的交互輸出作為主動(dòng)段濾波器的輸入，自由段采用橢球體地球重力模型［8］，再入段采用基于橢球體地球重力模型和空氣阻力模型的運(yùn)動(dòng)模型［9～10］。各階段運(yùn)動(dòng)模型都建立在混合坐標(biāo)系［10］下，在雷達(dá)站坐標(biāo)系下建立彈道目標(biāo)各運(yùn)動(dòng)階段的狀態(tài)方程，在雷達(dá)站球坐標(biāo)系下建立各運(yùn)動(dòng)階段的量測(cè)方程。本文基本濾波算法采用擴(kuò)展卡爾曼濾波（EKF）算法［11］。

2 彈道目標(biāo)運(yùn)動(dòng)模型

本文根據(jù)各階段受力情況應(yīng)用了不同的運(yùn)動(dòng)模型，由于各階段運(yùn)動(dòng)模型的狀態(tài)變量的維度不同，無(wú)法直接用于變結(jié)構(gòu)多模型算法中模型的交互，所以將各個(gè)狀態(tài)變量分為包括目標(biāo)位置和速度的六維主狀態(tài)變量和包括其他狀態(tài)變量的次狀態(tài)變量，主狀態(tài)變量用于多模型之間的交互，次狀態(tài)變量?jī)H在濾波器內(nèi)部循環(huán)。本文中所有運(yùn)動(dòng)模型都建立在混合坐標(biāo)系下，在雷達(dá)站坐標(biāo)系下建立彈道目標(biāo)狀態(tài)方程，在雷達(dá)站球坐標(biāo)系下建立量測(cè)方程。

2.1 主動(dòng)段運(yùn)動(dòng)模型

主動(dòng)段主要受到發(fā)動(dòng)機(jī)推力、地球重力和空氣動(dòng)力的作用。使用Singer模型和重力轉(zhuǎn)彎模型分別描述該階段的運(yùn)動(dòng)情況，運(yùn)動(dòng)模型分別如下。

2.1.1 Singer模型

設(shè)彈道目標(biāo)Singer模型［12］的狀態(tài)變量如式（1），各變量依次表示目標(biāo)在各個(gè)方向的位置、速度和加速度。

Singer模型在雷達(dá)站坐標(biāo)系下的狀態(tài)方程為

將式（1）的狀態(tài)向量分離為主狀態(tài)向量X1和次狀態(tài)向量X2：

該運(yùn)動(dòng)模型中其他矩陣向量參考文獻(xiàn)［12］。

2.1.2 重力轉(zhuǎn)彎模型

設(shè)彈道目標(biāo)的狀態(tài)變量為

彈道目標(biāo)在雷達(dá)站坐標(biāo)系下的狀態(tài)方程如式（8）所示，式中各變量含義為，其中T表示發(fā)動(dòng)機(jī)推力，D表示空氣動(dòng)力，m表示彈道目標(biāo)質(zhì)量；β=δm/m，其中δm>0，表示導(dǎo)彈質(zhì)量減小的速率；μ=3.986005×1014m3/s2，為地球重力常數(shù)；Reh=rem+h，rem=6371010m，為地球半徑，h為雷達(dá)站的大地高程。

速度v和位移r0表達(dá)式如下：

對(duì)狀態(tài)方程進(jìn)行泰勒級(jí)數(shù)展開(kāi)并忽略高階項(xiàng)有：

狀態(tài)方程的雅可比矩陣求解公式為

將狀態(tài)向量分離為主狀態(tài)向量X(1)和次狀態(tài)向量X(2)，表達(dá)式如下所示：

將狀態(tài)方程、泰勒級(jí)數(shù)展開(kāi)的函數(shù)和雅可比矩陣也相應(yīng)分離，分別在濾波器內(nèi)迭代計(jì)算。

2.2 自由段運(yùn)動(dòng)模型

自由段是指彈頭從頭體分離到再次進(jìn)入大氣層的一段彈道，這一階段目標(biāo)近似于在真空中運(yùn)動(dòng)，可認(rèn)為僅受到地球重力的作用，為使運(yùn)動(dòng)模型更接近真實(shí)情況，地球重力采用橢球地球模型。

狀態(tài)方程的雅可比矩陣求解公式與式（7）相同。

2.3 再入段運(yùn)動(dòng)模型

再入段是指彈頭從再次進(jìn)入大氣層到落地為止的一段彈道，這一階段彈頭除受到地球重力之外還受到巨大的氣動(dòng)阻力。因此，將氣動(dòng)阻力α增廣到狀態(tài)向量中得到表達(dá)式如下：

雷達(dá)站坐標(biāo)系下再入段彈道目標(biāo)的狀態(tài)方程為

式中，α表示阻力參數(shù)，ρ(h)=ρ0e-kh為空氣密度函數(shù)，其中ρ0=1.22kgm-3。

對(duì)上式進(jìn)行泰勒級(jí)數(shù)展開(kāi)并忽略高階項(xiàng)有：

狀態(tài)方程的雅可比矩陣求解公式與式（7）相同。

將狀態(tài)向量分離為主狀態(tài)向量X(1)和次狀態(tài)向量X(2)，表達(dá)式如下所示：

對(duì)狀態(tài)方程、泰勒級(jí)數(shù)展開(kāi)的函數(shù)和雅可比矩陣的處理與主動(dòng)段中的重力轉(zhuǎn)彎模型相同。

3 雷達(dá)量測(cè)模型

狀態(tài)向量在雷達(dá)站坐標(biāo)系下描述，量測(cè)值在雷達(dá)站球坐標(biāo)系獲取，其中R為距離，A為方位角，R為俯仰角。則用狀態(tài)變量表達(dá)的量測(cè)方程如下：

對(duì)上式進(jìn)行泰勒展開(kāi)并忽略高階項(xiàng)有：

式（11）中H(Xk|k)為量測(cè)方程的雅可比矩陣，表達(dá)式如下：

4 彈道目標(biāo)的變結(jié)構(gòu)多模型算法

模型集切換算法是變結(jié)構(gòu)多模型的一種模型集自適應(yīng)算法，該算法的優(yōu)勢(shì)是模型集之間的切換通過(guò)兩階段進(jìn)行，在激活候選模型集后不立刻終止當(dāng)前模型集，而是運(yùn)行二者的并集直至當(dāng)前模型集達(dá)到設(shè)定的終止條件。這種切換方式可以有效地避免虛假切換。本文在模型集切換算法的基礎(chǔ)上，設(shè)計(jì)了一個(gè)適用于彈道目標(biāo)全階段的數(shù)據(jù)濾波方法，相較于單模型跟蹤算法，該算法可用于跟蹤任意階段的彈道目標(biāo)。該方法的一個(gè)周期由以下七個(gè)步驟組成。

1）設(shè)計(jì)模型集的覆蓋。彈道目標(biāo)主動(dòng)段是整個(gè)飛行階段受力最復(fù)雜的階段，為了更準(zhǔn)確地對(duì)主動(dòng)段數(shù)據(jù)濾波，本文采用Singer模型和重力轉(zhuǎn)彎模型的交互多模型算法描述該階段運(yùn)動(dòng)情況，兩個(gè)模型分別用m1、m2表示；彈道目標(biāo)自由段在三個(gè)運(yùn)動(dòng)階段中受力情況最簡(jiǎn)單，運(yùn)動(dòng)軌跡可預(yù)測(cè)性更強(qiáng)，該階段采用橢球體地球重力模型，用m3表示；彈道目標(biāo)再入段主要受到地球重力和空氣阻力的作用，采用基于橢球地球模型和空氣動(dòng)力模型的運(yùn)動(dòng)模型描述，用m4表示。

2）設(shè)計(jì)候選模型集。模型集的設(shè)計(jì)與整體模型集合以及它們的拓?fù)渚o密相關(guān)。由各個(gè)模型代表的系統(tǒng)模式很接近時(shí)，這些模型可以聚成一個(gè)模型集。彈道目標(biāo)飛行所經(jīng)歷的運(yùn)動(dòng)階段依次為主動(dòng)段、自由段和再入段。設(shè)計(jì)各個(gè)階段模型集的集合為C={M1，M2，M3}。其中，M1={m1，m2}，為主動(dòng)段模型集；M2={m3}，為自由段模型集；M3={m4}，為再入段模型集。

3）判斷彈道目標(biāo)當(dāng)前飛行階段。當(dāng)雷達(dá)探測(cè)到彈道目標(biāo)的數(shù)據(jù)時(shí)，首先要做的事情是判斷彈道目標(biāo)當(dāng)前處于哪個(gè)階段。初始模型集的激活只能依靠從量測(cè)序列中獲得的后驗(yàn)信息，通過(guò)后驗(yàn)信息總結(jié)出模型概率和似然，從而判斷當(dāng)前彈道目標(biāo)處于哪個(gè)運(yùn)動(dòng)階段。初始模型集的激活需要同時(shí)滿足以下兩個(gè)條件：

其中，模型集概率是k時(shí)刻基于計(jì)算的模型集中所有模型的概率之和；模型集似然為k時(shí)刻模型集中模型的所有（邊緣）似然的概率加權(quán)和。

4）時(shí)間計(jì)數(shù)k增加1，運(yùn)行變結(jié)構(gòu)交互多模型算法遞歸周期［13］。

5）判斷是否激活候選模型集。如果是，則進(jìn)入下一步；如果不是，則輸出從變結(jié)構(gòu)交互多模型算法周期獲得的目標(biāo)狀態(tài)估計(jì)值k|k、狀態(tài)誤差協(xié)方差Pk|k和模型概率，返回4）。

設(shè)Mk為當(dāng)前運(yùn)行的模型集，當(dāng)候選模型集Mk+1同時(shí)滿足以下兩個(gè)條件時(shí)，激活Mk+1。

t0為該激活邏輯中的門(mén)限值。

6）如果有候選模型集被激活，下一時(shí)刻的模型集變?yōu)楫?dāng)前時(shí)刻運(yùn)行的模型集M0和新激活的候選模型集Ma的并集Mk，并且運(yùn)行變結(jié)構(gòu)交互多模型算法的一個(gè)周期。

（1）計(jì)算并集Mk的模式概率、狀態(tài)矩陣和誤差協(xié)方差矩陣。

（2）判斷是否終止模型集M0和Ma，對(duì)模型集Ml=Mo，Ma計(jì)算下列模型概率和模型似然。

7）選擇模型轉(zhuǎn)移概率矩陣。綜合考慮彈道目標(biāo)的飛行特性，定義模型轉(zhuǎn)換矩陣如下：

模型飛行階段切換示意圖如圖1。

圖1 彈道目標(biāo)飛行階段模型轉(zhuǎn)換圖

5 仿真結(jié)果

以某兩級(jí)彈道導(dǎo)彈為仿真對(duì)象，假設(shè)發(fā)射點(diǎn)的大地坐標(biāo)為(130°E，30°N )，發(fā)射方位角為 40°。仿真后彈道目標(biāo)的落點(diǎn)大地坐標(biāo)為(95°W，29°55′N(xiāo) )。假設(shè)雷達(dá)站位于地球表面，大地坐標(biāo)為(95°W，30°N )。設(shè)定量測(cè)向量的距離量測(cè)方差為1m，方位角和俯仰角的量測(cè)方差均為10-3rad，跟蹤數(shù)據(jù)率為1Hz。將本文提出的基于變結(jié)構(gòu)多模型數(shù)據(jù)濾波算法與單模型EKF算法進(jìn)行對(duì)比分析，采用位置和速度的均方根誤差（Root-Mean-Square Error，RMSE）作為標(biāo)準(zhǔn)對(duì)比分析濾波效果。

5.1 主動(dòng)段仿真結(jié)果對(duì)比

從圖2和圖3主動(dòng)段的彈道目標(biāo)數(shù)據(jù)濾波算法的估計(jì)精度對(duì)比可以看到，基于變結(jié)構(gòu)多模型算法的數(shù)據(jù)濾波方法的位置和速度的RMSE都優(yōu)于基于EKF算法的單模型數(shù)據(jù)濾波方法，在50s左右變結(jié)構(gòu)多模型算法逐漸收斂，而另外兩種單模型算法的誤差值仍然是波動(dòng)上升的趨勢(shì)。從表1可以得出，變結(jié)構(gòu)多模型算法與基于Singer的EKF算法的位置誤差和速度誤差的比值分別為0.803和0.787，與基于重力轉(zhuǎn)彎模型的EKF算法的位置誤差和速度誤差比值的分別為0.835和0.673。

圖2 主動(dòng)段的位置均方根誤差對(duì)比

圖3 主動(dòng)段的速度均方根誤差對(duì)比

表1 主動(dòng)段的均方根誤差均值對(duì)比

5.2 自由段仿真結(jié)果對(duì)比

從圖4和圖5自由段的彈道目標(biāo)數(shù)據(jù)濾波算法的估計(jì)精度對(duì)比可以看到，基于變結(jié)構(gòu)多模型算法的數(shù)據(jù)濾波方法的位置和速度的均方根誤差都優(yōu)于基于EKF算法的單模型數(shù)據(jù)濾波方法，在50s左右變結(jié)構(gòu)多模型算法逐漸收斂，而EKF單模型算法在200s時(shí)才逐漸收斂。從表2可以得出，變結(jié)構(gòu)多模型算法與EKF算法的位置誤差和速度誤差比值的分別為0.601和0.654。

圖4 自由段的位置均方根誤差對(duì)比

圖5 自由段的速度均方根誤差對(duì)比

表2 由段的均方根誤差均值對(duì)比

5.3 再入段仿真結(jié)果對(duì)比

從圖6和圖7再入段的彈道目標(biāo)數(shù)據(jù)濾波算法的估計(jì)精度對(duì)比可以看到，基于變結(jié)構(gòu)多模型算法的數(shù)據(jù)濾波方法的位置和速度的均方根誤差都優(yōu)于基于EKF算法的單模型數(shù)據(jù)濾波方法，在20s左右變結(jié)構(gòu)多模型算法逐漸收斂，而另外兩種單模型算法的誤差值仍然是波動(dòng)上升的趨勢(shì)。從表3可以得出，變結(jié)構(gòu)多模型算法與EKF算法的位置誤差和速度誤差比值的分別為0.536和0.740。

表3 再入段的均方根誤差均值對(duì)比

圖6 再入段的位置均方根誤差對(duì)比

圖7 再入段的速度均方根誤差對(duì)比

6 結(jié)語(yǔ)

本文提出了一種基于變結(jié)構(gòu)多模型算法的彈道目標(biāo)數(shù)據(jù)濾波方法，可以實(shí)現(xiàn)對(duì)彈道目標(biāo)全彈道的數(shù)據(jù)濾波。通過(guò)仿真實(shí)驗(yàn)得出，與基于單模型的EKF算法相比，本文所提方法有以下優(yōu)勢(shì)：各運(yùn)動(dòng)階段的位置和速度均方根誤差均值都更小，即該算法估計(jì)精度更高；當(dāng)濾波數(shù)據(jù)出現(xiàn)波動(dòng)時(shí)，本文方法在各運(yùn)動(dòng)階段的收斂速度都更快，即該算法更穩(wěn)定；在主動(dòng)段和自由段、自由段和再入段的交替處基于變結(jié)構(gòu)多模型算法的彈道目標(biāo)數(shù)據(jù)濾波方法的均方根誤差值抖動(dòng)范圍更小，抖動(dòng)時(shí)間更短，這表明該算法在彈道目標(biāo)出現(xiàn)機(jī)動(dòng)時(shí)能夠更快適應(yīng)。