基于EKF及彈道方程的彈道目標(biāo)跟蹤濾波器設(shè)計

孫照強，王志貴，孟 飛，李陸雨，于 中，陳 燕

（北京無線電測量研究所，北京 100854）

0 引 言

彈道導(dǎo)彈是現(xiàn)代高技術(shù)戰(zhàn)爭中的重要作戰(zhàn)武器，為了實現(xiàn)有效的攻擊，彈道導(dǎo)彈常采用以下突防方式：飽和式、多波次攻擊，多彈頭突防，釋放輕、重誘餌等。為了實現(xiàn)對彈道導(dǎo)彈的有效防御，地面跟蹤制導(dǎo)雷達必須具備跟蹤多批次目標(biāo)的能力以及單批次團目標(biāo)跟蹤能力。另外，為了對來襲目標(biāo)進行識別并對制導(dǎo)攔截彈進行攔截，又要求雷達具備較高的跟蹤精度。然而，由于雷達資源有限，當(dāng)雷達跟蹤目標(biāo)數(shù)目較多時，需要降低跟蹤數(shù)據(jù)率以保證對當(dāng)前照射區(qū)域內(nèi)所有目標(biāo)的穩(wěn)定跟蹤，這就出現(xiàn)了跟蹤數(shù)據(jù)率和跟蹤精度之間的矛盾。例如，地面制導(dǎo)雷達當(dāng)前時刻對10批目標(biāo)進行跟蹤，而每批目標(biāo)團內(nèi)又有10個目標(biāo)，此時需要對100個目標(biāo)同時進行濾波處理，若采用常規(guī)的濾波算法，計算量將急劇上升，不利于工程實現(xiàn)。因此，迫切需要研究一種在低數(shù)據(jù)率下具有高精度跟蹤能力、且計算量適中的彈道目標(biāo)濾波方法。

彈道導(dǎo)彈從發(fā)射點到落點的整個軌跡通常分為兩個階段：主動段和被動段［1］。在導(dǎo)彈整個飛行過程中，其會受重力、推力、氣動阻力、地球自轉(zhuǎn)偏向力的公共作用，其運動特征表現(xiàn)為高度非線性，對其跟蹤的核心問題是設(shè)計精確的非線性濾波器，對其運動狀態(tài)進行估計。高精度彈道導(dǎo)彈跟蹤需要滿足兩個條件：精確的目標(biāo)運動模型以及性能優(yōu)良的濾波方法［2］。

在彈道導(dǎo)彈運動模型方面，動力學(xué)建模方法是目前公認的較為準(zhǔn)確的一種描述目標(biāo)飛行過程中受力情況的方法［3-5］。在對被動段彈道導(dǎo)彈跟蹤中，使用較多的是球形模型［6-7］，然而當(dāng)雷達采樣率較低并且彈道導(dǎo)彈目標(biāo)射程較遠時，常用的球形模型會無法滿足高精度跟蹤需求，此時需要采用基于橢球J2修正的重力加速度模型［8-10］。

在彈道導(dǎo)彈濾波方法方面，一般采用一些對于非線性系統(tǒng)具有良好適應(yīng)性的濾波算法［11-13］，其中具有代表性的算法有擴展卡爾曼濾波（extended Kalman filter，EKF）［14-15］、無跡卡爾曼濾波（unscented Kalman filter，UKF）［16-18］和粒子濾波（particle filter，PF）［19-21］。EKF算法采用泰勒級數(shù)展開方法將非線性濾波過程轉(zhuǎn)化為一個近似線性問題，工程實現(xiàn)簡單，應(yīng)用范圍廣。UKF算法采用Sigma點進行非線性逼近，可以取得三階泰勒展開精度［22］，但是其濾波穩(wěn)定性易受中心采樣點的權(quán)值影響，產(chǎn)生較大的波動［23］。PF方法一般用于非高斯系統(tǒng)，但是其計算量較大，無法進行工程應(yīng)用。文獻［24］對彈道目標(biāo)跟蹤的4種濾波器的性能進行了比較，然而，目標(biāo)運動模型采用了分段常加速模型，對彈道導(dǎo)彈運動過程描述不夠精確，導(dǎo)致濾波精度受限。文獻［25］針對彈道系數(shù)己知的再入彈道目標(biāo)進行了跟蹤，對不同濾波器性能進行了仿真比較，但是其僅對目標(biāo)運動過程中受到的重力與空氣阻力進行了建模，未考慮地球自轉(zhuǎn)帶來的影響，目標(biāo)運動模型不夠精確。文獻［26］提出了一種迭代線性化彈道導(dǎo)彈被動段跟蹤方法，在混合坐標(biāo)系下完成了彈道導(dǎo)彈的跟蹤濾波。文獻［27］利用多模型對彈道導(dǎo)彈在不同飛行階段的運動狀態(tài)進行描述，并采用粒子濾波方法對其進行跟蹤濾波，但是其具有較大的運算量。文獻［28］針對彈道導(dǎo)彈目標(biāo)提出一種基于高斯粒子濾波的狀態(tài)依賴交互多模型方法，對彈道導(dǎo)彈3個飛行階段進行跟蹤，取得了較好的效果。

綜上所述，考慮到彈道導(dǎo)彈跟蹤時數(shù)據(jù)率和精度之間的矛盾，以及目前所用模型不夠精確的問題，本文針對被動段飛行的彈道導(dǎo)彈跟蹤展開研究。首先建立了彈道導(dǎo)彈被動段的精確跟蹤模型，隨后詳細推導(dǎo)了基于彈道運動方程的EKF濾波過程，使雷達在較低數(shù)據(jù)率下仍然保持較高的跟蹤精度，同時具有較小的計算量，最后通過仿真實驗對所提的運動模型及濾波方法進行了驗證，結(jié)果表明本文所提算法計算量小，且跟蹤精度滿足當(dāng)前制導(dǎo)雷達跟蹤需求。

1 導(dǎo)彈被動段質(zhì)心運動方程

本文主要研究導(dǎo)彈在被動段的跟蹤過程，為此本節(jié)首先給出彈道目標(biāo)在被動段的運動模型。

考慮到雷達測量值為雷達大地直角坐標(biāo)系下的目標(biāo)徑向距離、方位角、俯仰角，因此將導(dǎo)彈運動方程建立在雷達大地直角坐標(biāo)系下，這樣可以減少坐標(biāo)變換過程，提高濾波效率。為了提高低數(shù)據(jù)率濾波精度，本文將考慮導(dǎo)彈在被動段所受到的諸多作用力，以盡可能建立與實際運動情況相符的運動模型。

如圖1所示，以地心OE為坐標(biāo)原點建立地心慣性坐標(biāo)系。雷達位于O點，P點為彈道上的一點，其在雷達大地直角坐標(biāo)系中的坐標(biāo)為（x，y，z），將點P的地心矢徑r表示為r＝R0＋ρ，其模r＝r，其中R0為地球半徑，ρ為雷達距P點的徑向距離。

圖1 彈道上點P和雷達點的地心矢徑關(guān)系圖Fig.1 Relation between the trajectory point P and geocentric arrow diameter of the radar

設(shè)導(dǎo)彈相對雷達大地直角坐標(biāo)系的速度V＝（V x，V y，V z）T，其模V＝V。在地面雷達看來，該階段導(dǎo)彈受到地心引力、離心慣性力、科氏慣性力、空氣阻力的共同影響，其在雷達大地直角坐標(biāo)系中的導(dǎo)彈質(zhì)心動力學(xué)方程［29-30］為

和gωe分別是引力加速度g在地心矢徑方向和地球自轉(zhuǎn)角速度方向的分量；?為點P的地心緯度值；J2為攝動常數(shù)；ωe為地球自轉(zhuǎn)角速度；μ為引力常數(shù)；ρ為目標(biāo)位置的大氣密度值，按照美國標(biāo)準(zhǔn)大氣76模型［31］計算；β為目標(biāo)的質(zhì)阻比。

2 基于彈道運動方程的EKF濾波算法模型

雷達量測值為雷達大地球坐標(biāo)系下的R、Az、E。為了估計目標(biāo)質(zhì)阻比，將質(zhì)阻比表示為β＝β0eγ，β0為設(shè)置的初始質(zhì)阻比大小，γ為一待估分量。將導(dǎo)彈的運動方程建立在雷達大地直角坐標(biāo)系下，狀態(tài)變量為

導(dǎo)彈在雷達大地直角坐標(biāo)系中的狀態(tài)微分方程為

為了應(yīng)用EKF的濾波算法，需要求出狀態(tài)函數(shù)和量測函數(shù)的Jacobi矩陣。

狀態(tài)函數(shù)Jacobi矩陣為

當(dāng)?shù)玫綘顟B(tài)函數(shù)和量測函數(shù)的Jacobi矩陣后，結(jié)合EKF［32］過程即可得到最終的基于彈道運動方程的EKF算法，詳細濾波過程如下：

（1）狀態(tài)預(yù)測

式中：Q（k）為系統(tǒng)噪聲方差陣；Φ（k）為狀態(tài)轉(zhuǎn)移矩陣，在本文中可表示為

（2）狀態(tài)估計

式中：增益K（k＋1）可表示為

式中：R（k＋1）為量測噪聲方差陣。

本文研究的基于彈道運動方程的EKF算法，采用了基于橢球地球模型的J2修正重力加速度模型，同時考慮了該被動段導(dǎo)彈受地心引力、離心慣性力、科氏慣性力、空氣阻力的共同影響，建立了精確的導(dǎo)彈質(zhì)心動力學(xué)方程。將該精確模型與非線性EKF濾波結(jié)合從而實現(xiàn)了彈道目標(biāo)的跟蹤精度。

3 BEKF濾波算法的仿真和驗證

3.1 彈道目標(biāo)運動場景

仿真用的彈道目標(biāo)其射程為1 000 km，地面雷達觀測的徑向距離、方位角、俯仰角以及全速度變化曲線如圖2～圖5所示，目標(biāo)再入速度最大值接近3 km／s。

圖2 目標(biāo)距離與時間的變化情況Fig.2 Variation of target range with time

圖3 目標(biāo)方位角與時間的變化情況Fig.3 Variation of target azimuth with time

圖4 目標(biāo)仰角與時間的變化情況Fig.4 Variation of target pitch angle with time

圖5 目標(biāo)全速度與時間的變化情況Fig.5 Variation of target velocity with time

3.2 實驗驗證

采用上述目標(biāo)場景，設(shè)定目標(biāo)質(zhì)阻比初值為6 000 kg／m2。設(shè)雷達距離、方位和俯仰量測噪聲服從相互獨立的零均值正態(tài)分布，其方差根據(jù)雷達威力、目標(biāo)散射界面等計算得到。為了驗證本文提出算法的有效性和優(yōu)點，在雷達跟蹤數(shù)據(jù)率為5 Hz情況下分別對傳統(tǒng)的基于常加速模型的EKF算法、基于彈道運動方程的UKF算法以及本文所提的基于彈道運動方程的EKF算法進行比較，其中基于彈道運動方程的UFK方法采樣點數(shù)為13，尺度參數(shù)設(shè)為0.1。以上3種方法分別進行30次蒙特卡羅仿真，結(jié)果如圖6～圖8所示。

圖6 距離的均方根誤差曲線Fig.6 Range root mean square error

圖7 方位角的均方根誤差曲線Fig.7 Azimuth root mean square error

圖8 俯仰角的均方根誤差曲線Fig.8 Pitch angle root mean square error

從實驗結(jié)果中可以看出：

（1）本文提出的基于彈道運動方程的EKF算法，濾波精度和收斂速度明顯優(yōu)于傳統(tǒng)的基于常加速模型的EKF算法，與基于彈道運動方程的UKF算法相當(dāng)；

（2）本文所提的算法在濾波精度上與基于彈道運動方程的UKF算法相當(dāng)，考慮到計算時間和計算穩(wěn)定性的雙重優(yōu)勢，基于彈道運動方程的EKF更具有工程應(yīng)用價值。

4 結(jié) 論

針對地面跟蹤雷達在低數(shù)據(jù)率下對多目標(biāo)的高精度需求，本文提出了基于彈道運動方程的EKF算法，基于精確的彈道導(dǎo)彈被動段質(zhì)心運動方程，推導(dǎo)了EKF濾波過程。通過與傳統(tǒng)的基于常加速模型的EKF算法和基于彈道運動方程的UKF算法比較，驗證了基于彈道運動方程的EKF具有低數(shù)據(jù)率下濾波精度高、計算量小等優(yōu)點，解決了地面跟蹤雷達實際中遇到的問題，具有工程應(yīng)用價值。