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

黃晶晶 陳世友

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

1 引言

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

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

2 彈道目標運動模型

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

2.1 主動段運動模型

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

2.1.1 Singer模型

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

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

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

該運動模型中其他矩陣向量參考文獻［12］。

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

設彈道目標的狀態(tài)變量為

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

速度v和位移r0表達式如下：

對狀態(tài)方程進行泰勒級數(shù)展開并忽略高階項有：

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

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

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

2.2 自由段運動模型

自由段是指彈頭從頭體分離到再次進入大氣層的一段彈道，這一階段目標近似于在真空中運動，可認為僅受到地球重力的作用，為使運動模型更接近真實情況，地球重力采用橢球地球模型。

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

2.3 再入段運動模型

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

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

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

對上式進行泰勒級數(shù)展開并忽略高階項有：

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

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

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

3 雷達量測模型

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

對上式進行泰勒展開并忽略高階項有：

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

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

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

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

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

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

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

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

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

設Mk為當前運行的模型集，當候選模型集Mk+1同時滿足以下兩個條件時，激活Mk+1。

t0為該激活邏輯中的門限值。

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

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

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

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

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

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

5 仿真結(jié)果

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

5.1 主動段仿真結(jié)果對比

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

圖2 主動段的位置均方根誤差對比

圖3 主動段的速度均方根誤差對比

表1 主動段的均方根誤差均值對比

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

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

圖4 自由段的位置均方根誤差對比

圖5 自由段的速度均方根誤差對比

表2 由段的均方根誤差均值對比

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

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

表3 再入段的均方根誤差均值對比

圖6 再入段的位置均方根誤差對比

圖7 再入段的速度均方根誤差對比

6 結(jié)語

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