ECEF坐標(biāo)系彈道導(dǎo)彈跟蹤研究

許登榮，程水英，包守亮

(電子工程學(xué)院 504教研室, 合肥 230037)

?

·數(shù)據(jù)處理·

ECEF坐標(biāo)系彈道導(dǎo)彈跟蹤研究

許登榮，程水英，包守亮

(電子工程學(xué)院 504教研室, 合肥 230037)

針對單模型跟蹤算法無法實(shí)現(xiàn)對彈道導(dǎo)彈連續(xù)精確跟蹤的問題，提出了一種地心地固坐標(biāo)系下彈道導(dǎo)彈全階段連續(xù)跟蹤的交互式多模型(IMM)算法。根據(jù)導(dǎo)彈不同飛行階段的受力情況建立IMM算法的兩個濾波模型集：助推段采用恒軸向力模型，中段、再入段采用被動段跟蹤模型，兩個模型都利用濾波精度高、數(shù)值穩(wěn)定性較好的求容積卡爾曼濾波算法進(jìn)行濾波。針對該IMM算法在混合估計(jì)過程中引入偏差的問題，采用無偏混合的方法，以兩種典型的跟蹤場景進(jìn)行仿真校驗(yàn)，實(shí)驗(yàn)結(jié)果驗(yàn)證了該文算法的優(yōu)越性。

機(jī)動目標(biāo)跟蹤；彈道導(dǎo)彈；交互式多模型算法；求容積卡爾曼濾波

0 引 言

彈道導(dǎo)彈自問世以來，因其具有射程遠(yuǎn)、速度快、突防能力強(qiáng)、打擊精度高、殺傷威力大、效費(fèi)比高等優(yōu)點(diǎn)受到全世界的關(guān)注，針對它的防御研究也從未停止過。彈道導(dǎo)彈跟蹤是彈道導(dǎo)彈防御的核心問題之一，跟蹤的精度決定了落點(diǎn)預(yù)測的精度以及導(dǎo)彈攔截的成功率。彈道導(dǎo)彈的飛行階段通?？梢苑譃橹鲃佣魏捅粍佣?，其中，主動段又叫作助推段，被動段又可分為中段和再入段。大多數(shù)有關(guān)導(dǎo)彈跟蹤的文獻(xiàn)通常只對導(dǎo)彈某一飛行階段進(jìn)行研究[1]，但采用這種普通方法不能實(shí)現(xiàn)導(dǎo)彈的連續(xù)精確跟蹤。近些年來，一些文獻(xiàn)提出利用交互式多模型算法對彈道導(dǎo)彈進(jìn)行連續(xù)跟蹤[2-4]，減小了導(dǎo)彈不同飛行階段過渡處的跟蹤誤差，提高了導(dǎo)彈的跟蹤精度。但多模型算法性能的好壞很大程度上取決于所選的模型集[5]，針對大多數(shù)導(dǎo)彈跟蹤文獻(xiàn)中建立的模型較為粗糙的問題，本文根據(jù)導(dǎo)彈不同飛行階段的受力情況建立交互式多模型(IMM)算法較為精確的濾波模型集，實(shí)現(xiàn)對彈道導(dǎo)彈全階段的連續(xù)高精度跟蹤。

此外，彈道導(dǎo)彈跟蹤存在一個坐標(biāo)系選擇的問題。合理而恰當(dāng)?shù)剡x擇參考坐標(biāo)系會使描述導(dǎo)彈運(yùn)動規(guī)律的數(shù)學(xué)模型大為簡化，否則會導(dǎo)致問題復(fù)雜化，甚至陷入無法處理的困境[6]。目標(biāo)運(yùn)動的描述通常在笛卡爾坐標(biāo)系中，但量測值卻在傳感器坐標(biāo)系中描述。因而，對于導(dǎo)彈跟蹤來說，坐標(biāo)系的選擇存在四種可能性[7]，即混合坐標(biāo)系、笛卡爾坐標(biāo)系、傳感器坐標(biāo)系，以及其他坐標(biāo)系。文獻(xiàn)[7-8]對不同坐標(biāo)系中跟蹤的優(yōu)缺點(diǎn)進(jìn)行了分析。目前大部分的導(dǎo)彈跟蹤算法都選擇在混合坐標(biāo)系中進(jìn)行處理[8]，即在東-北-天(ENU)坐標(biāo)系中建立目標(biāo)的狀態(tài)方程，量測值在雷達(dá)球坐標(biāo)系中進(jìn)行描述。但是在ENU坐標(biāo)系中描述導(dǎo)彈運(yùn)動的數(shù)學(xué)模型較為復(fù)雜，而在地心地固(ECEF)坐標(biāo)系中不僅對導(dǎo)彈目標(biāo)進(jìn)行受力分析非常的方便，而且ECEF坐標(biāo)系中目標(biāo)的狀態(tài)方程也十分簡潔。此外，全局的戰(zhàn)場態(tài)勢演示等在ECEF坐標(biāo)系中更為形象和直觀。鑒于此，本文選擇在ECEF坐標(biāo)系中對導(dǎo)彈目標(biāo)進(jìn)行跟蹤。

本文采用IMM算法實(shí)現(xiàn)無先驗(yàn)信息條件下對彈道導(dǎo)彈全階段的連續(xù)跟蹤，在對導(dǎo)彈各階段受力分析的基礎(chǔ)上，在ECEF坐標(biāo)系中建立了IMM算法的兩個濾波器模型。主推段目標(biāo)的軸向力近似保持不變[1，4]，采用恒軸向力(CAF)[1-2]模型；中段和再入段采用將彈道系數(shù)擴(kuò)展為目標(biāo)狀態(tài)分量進(jìn)行聯(lián)合估計(jì)的被動段跟蹤模型，記為PM模型。反導(dǎo)雷達(dá)的量測在雷達(dá)站球坐標(biāo)系獲得，因而需要將量測轉(zhuǎn)換到ECEF坐標(biāo)系中。這個過程分為兩步：先將雷達(dá)球坐標(biāo)系量測轉(zhuǎn)換到ENU直角坐標(biāo)系中，再轉(zhuǎn)換到ECEF坐標(biāo)系中。但若球坐標(biāo)系轉(zhuǎn)換到直角坐標(biāo)系是有偏的，會造成濾波器的性能下降。針對此問題，本文先采用修正的無偏量測轉(zhuǎn)換(MUCM)方法[9-10]將雷達(dá)量測轉(zhuǎn)換為ENU坐標(biāo)系偽線性量測，再將ENU坐標(biāo)系量測和相應(yīng)的協(xié)方差轉(zhuǎn)換到ECEF中，并給出了轉(zhuǎn)換過程的推導(dǎo)，從而建立了對導(dǎo)彈全階段連續(xù)跟蹤所需的狀態(tài)方程和量測方程。由于狀態(tài)方程非線性，本文在濾波實(shí)現(xiàn)上采用數(shù)值精度較高、穩(wěn)定性較好的求容積卡爾曼濾波(CKF)算法[11]，同時比較了幾種主流的非線性濾波算法在彈道目標(biāo)跟蹤的性能。另外，由于該IMM算法的兩個濾波模型狀態(tài)維數(shù)不同，而要現(xiàn)實(shí)混合交互，需要將各模型的狀態(tài)向量轉(zhuǎn)化一致[12]，但在混合的過程中引入了偏差，針對這個問題，采用無偏混合[13]的方法。仿真結(jié)果表明：該算法能實(shí)現(xiàn)對彈道導(dǎo)彈連續(xù)精確跟蹤，且不同飛行階段過渡處的跟蹤誤差較小。

1 彈道導(dǎo)彈跟蹤模型

1.1 運(yùn)動模型

由于彈道導(dǎo)彈受力復(fù)雜且不同階段的受力不一樣，很難用統(tǒng)一的數(shù)學(xué)模型進(jìn)行描述，本文在受力分析的基礎(chǔ)上分別建立了彈道導(dǎo)彈不同階段的跟蹤濾波模型。

1.1.1 主動段跟蹤模型

ECEF坐標(biāo)系中，彈道導(dǎo)彈在主動段受到的作用力主要有推進(jìn)力、空氣阻力、地心引力和外在力(包括離心力和柯氏力)。由于推進(jìn)力和空氣阻力的一些參數(shù)與目標(biāo)具體型號有關(guān)，在無先驗(yàn)的情況下難以建立彈道導(dǎo)彈準(zhǔn)確的運(yùn)動模型。但推進(jìn)力和空氣阻力與目標(biāo)的速度幾乎在同一條直線上，這樣除重力以外的其他力可以合成一個軸向力，且基本保持不變。設(shè)ECEF坐標(biāo)系中目標(biāo)的狀態(tài)向量為xCAF=[xyzxyzab]T，位置矢量和速度矢量分別為p=[xyz]T，v=[vxvyvz]T。地球模型為橢球模型(WGS-84)，則助推段可以采用如下CAF模型對目標(biāo)進(jìn)行跟蹤

(1)

1.1.2 被動段跟蹤模型

被動段包括中段和再入段，和主動段相比，被動段不受推進(jìn)力的作用，即受到的作用力主要是地心引力、空氣阻力和外在力。其中，中段目標(biāo)受到的空氣阻力通?？梢院雎裕瑑H考慮重力和外在力的影響。而目標(biāo)重返大氣進(jìn)入再入段時，受到的作用力主要是地心引力和空氣動力，可以忽略外在力的影響，但由于空氣阻力與彈道系數(shù)有關(guān)且未知，需要對彈道系數(shù)進(jìn)行估計(jì)。因此，大多數(shù)有關(guān)文獻(xiàn)中段采用六狀態(tài)的跟蹤模型，再入段采用將彈道系數(shù)擴(kuò)展為目標(biāo)狀態(tài)的一個分量進(jìn)行實(shí)時估計(jì)的七狀態(tài)模型對目標(biāo)進(jìn)行跟蹤。為了實(shí)現(xiàn)對彈道目標(biāo)更加精確的跟蹤，本文中段和再入段考慮所有作用力的影響，同時為了減小IMM算法模型的數(shù)量，降低計(jì)算量，中段和再入段都采用如下七狀態(tài)的被動段模型對目標(biāo)進(jìn)行跟蹤。即設(shè)目標(biāo)的狀態(tài)變量為xPM=[xyzxyzα]T，則ECEF坐標(biāo)系中被動段彈道導(dǎo)彈的運(yùn)動方程為[1]

(2)

式中：ρ(h)=ρ0e-k0h為空氣密度，其中，h為目標(biāo)離地面的高度，ρ0=1.22 kg/m3，k0=1.414 1×10-4m-1；α=1/β為空氣阻力參數(shù)，β為彈道系數(shù)；wα(t)為零均值的高斯白噪聲。

到此，建立了彈道導(dǎo)彈各階段的運(yùn)動方程，可以看到，上述方程都是用微分方程進(jìn)行描述的，且都是非線性的。而為了利用相關(guān)的非線性濾波技術(shù)進(jìn)行遞推估計(jì)，通常需要將微分方程離散化，狀態(tài)方程離散化的方法有歐拉逼近法[14]、一階泰勒級數(shù)展開法[13]和四階龍格庫塔積分法[1,4]等。本文利用的是一階泰勒級數(shù)展開法。

1.2 量測方程

反導(dǎo)雷達(dá)的量測在球坐標(biāo)中獲得，量測值的表達(dá)式為

(3)

式中：xENU,k、yENU,k、zENU,k為k時刻目標(biāo)在雷達(dá)ENU坐標(biāo)系下的位置坐標(biāo)；zSPH, k=[Rm,Am,Em]T為雷達(dá)的觀測值，其中，Rm為目標(biāo)到雷達(dá)的距離，Am為方位角，Em為俯仰角；vSPH,k為零均值的高斯白噪聲序列，其協(xié)方差矩陣為RSPH,k。要實(shí)現(xiàn)遞推濾波估計(jì)，必須建立量測值與狀態(tài)變量之間的關(guān)系。本文先采用MUCM方法將球坐標(biāo)量測轉(zhuǎn)換為ENU坐標(biāo)系線性量測，然后將ENU坐標(biāo)系量測和相應(yīng)的協(xié)方差轉(zhuǎn)換到ECEF中，從而建立了濾波所需要的量測方程。設(shè)經(jīng)過MUCM方法轉(zhuǎn)換后的量測向量及協(xié)方差矩陣分別為zENU，k=[xMUCM,kyMUCM,kzMUCM,k]T，RENU,k=E((zENU,k-E(zENU,k))(zENU,k-E(zENU,k))T)，其轉(zhuǎn)換詳細(xì)過程可參考文獻(xiàn)[8]。

假設(shè)雷達(dá)站所在大地經(jīng)度為L，大地緯度為B，高程為H，地球?yàn)闄E球模型，則雷達(dá)站中心位置Or在ECEF坐標(biāo)系下的坐標(biāo)為

(4)

(5)

(6)

設(shè)由ENU坐標(biāo)系轉(zhuǎn)換到ECEF坐標(biāo)系的量測及協(xié)方差矩陣分別為zECEF,k，RECEF,k，則

(7)

RECEF,k= E((zECEF,k-E(zECEF,k))(zECEF,k-E(zECEF,k))T)=

(8)

式中：ρ=[xre,yre,zre]T。為了方便起見，省略下標(biāo)ECEF，即量測和協(xié)方差分別為zk，Rk。如無特殊申明，下文的量測及量測協(xié)方差矩陣均指轉(zhuǎn)換到ECEF坐標(biāo)系的量測及協(xié)方差。由此可建立如下線性的量測方程

zk=Hkxk+vk

(9)

式中：xk為k時刻的狀態(tài)向量；Hk為量測矩陣。二者都由跟蹤模型決定。vk為轉(zhuǎn)換后的量測噪聲，其協(xié)方差矩陣為Rk。采用橫軸向力模型時，狀態(tài)向量為xk=xCAF,k=[xyzxyzab]T，此時Hk的表達(dá)式為

(10)

至此，建立了導(dǎo)彈各飛行階段濾波所需的狀態(tài)方程和量測方程。由式(1)、式(2)以及式(9)可知，通過將量測轉(zhuǎn)換到ECEF坐標(biāo)系后，各跟蹤濾波器模型的狀態(tài)方程是非線性的，量測方程線性。針對非線性濾波問題，現(xiàn)有的方法有擴(kuò)展卡爾曼濾波(EKF)、無味卡爾曼濾波[15](UKF)、粒子濾波[16](PF)、CKF以及它們的改進(jìn)算法等。PF計(jì)算量較大，離工程實(shí)現(xiàn)較遠(yuǎn)；EKF通過在當(dāng)前狀態(tài)對非線性系統(tǒng)進(jìn)行一階泰勒級數(shù)展開來近似非線性函數(shù)本身，使非線性問題變?yōu)榫植烤€性問題，計(jì)算量小，但系統(tǒng)非線性較強(qiáng)時，濾波精度較差；當(dāng)系統(tǒng)非線性很強(qiáng)時，UKF精度好于EKF，且不用求Jacobian矩陣，但是對于高維非線性系統(tǒng)，對參數(shù)的選擇具有一定的依賴性；CKF與UKF類似，但CKF的參數(shù)(求容積點(diǎn)和權(quán)值)由狀態(tài)的維數(shù)唯一確定，實(shí)現(xiàn)更加方便，且計(jì)算量較小。綜合工程實(shí)現(xiàn)和濾波精度，本文采用CKF進(jìn)行濾波。

2 基于無偏混合估計(jì)的IMM算法

2.1 無偏混合估計(jì)

可以看到，IMM算法所采用的兩個跟蹤濾波器狀態(tài)維數(shù)不相同，且狀態(tài)變量的元素不一致，而IMM算法要實(shí)現(xiàn)狀態(tài)和協(xié)方差的混合交互，所有狀態(tài)向量必須轉(zhuǎn)換成一致的狀態(tài)變量，混合交互后再轉(zhuǎn)換為原來維數(shù)的狀態(tài)[12]。由于目標(biāo)的軸向力加速度和彈道系數(shù)是彈道目標(biāo)類型識別的一個重要參數(shù)，因此全系統(tǒng)的狀態(tài)變量為

(11)

假設(shè)助推段CAF模型、被動段(PM)模型(分別為模型1、2)k時刻的狀態(tài)估計(jì)值分別為

(12)

(13)

(14)

設(shè)混合概率為μi|jPr{rk=i|rk+1=j，Zk}，i,j∈{1,2,3}，Zk為直到k時刻的累積量測，則CAF模型的混合估計(jì)為

(15)

同理，PM模型的混合估計(jì)為

(16)

(17)

同時，根據(jù)修改后的狀態(tài)向量對協(xié)方差矩陣進(jìn)行相應(yīng)的調(diào)整，做上述修改后，各模型的混合估計(jì)為

(18)

這樣，混合后再轉(zhuǎn)換為原來維數(shù)的狀態(tài)作為濾波器的輸入初值進(jìn)行濾波就不會帶來偏差，IMM算法估計(jì)也具有更好的一致性。

2.2 算法步驟

對IMM算法做如上修改后，其他與傳統(tǒng)IMM算法相同。即可按照模型條件重初始化、模型條件濾波、模型概率更新、估計(jì)融合等步驟進(jìn)行遞推估計(jì)[5]，實(shí)現(xiàn)對導(dǎo)彈的穩(wěn)定跟蹤。圖1為本文提出的ECEF坐標(biāo)系下導(dǎo)彈連續(xù)跟蹤的IMM算法原理圖。

圖1 IMM算法原理圖

3 仿真結(jié)果與分析

實(shí)驗(yàn)采用MATLAB軟件進(jìn)行仿真計(jì)算，設(shè)計(jì)了兩個典型的跟蹤場景，即助推段和中段交接處的跟蹤以及中段和再入段交接處的跟蹤。彈道仿真采用四階龍格庫塔積分法外推生成標(biāo)準(zhǔn)彈道，雷達(dá)量測采用在彈道數(shù)據(jù)基礎(chǔ)上加量測噪聲。雷達(dá)的測距精度為σR=50 m，測角精度σA=σE=1 mrad，且量測噪聲服從高斯分布，整個跟蹤過程中雷達(dá)的測距和測角精度保持不變。

3.1 仿真場景1

以某近短程單級助推彈道導(dǎo)彈為跟蹤對象，導(dǎo)彈發(fā)射點(diǎn)為東經(jīng)0°，北緯30°，導(dǎo)彈發(fā)射方位角為80°，導(dǎo)彈的關(guān)機(jī)點(diǎn)時刻為85 s，目標(biāo)的彈道系數(shù)為5 000 kg/m2，并假設(shè)保持不變。雷達(dá)布站為東經(jīng)3.6°，北緯31°。圖2～圖4分別給出ECEF坐標(biāo)系下導(dǎo)彈的彈道曲線，以及導(dǎo)彈的加速度大小和速度大小變化曲線。雷達(dá)在導(dǎo)彈發(fā)射35 s后發(fā)現(xiàn)目標(biāo)，跟蹤雷達(dá)數(shù)據(jù)率為1 Hz，采樣點(diǎn)數(shù)為265。

圖2 ECEF坐標(biāo)系下導(dǎo)彈的彈道曲線

圖3 目標(biāo)的加速度變化曲線

圖4 目標(biāo)的速度變化曲線

由于助推段的機(jī)動性較強(qiáng)，為了驗(yàn)證本文算法的優(yōu)越性，在同等條件下與跟蹤機(jī)動目標(biāo)性能較好的Singer模型算法、CS模型算法進(jìn)行了比較，同時與文獻(xiàn)[3]提出的彈道導(dǎo)彈連續(xù)跟蹤方法進(jìn)行跟蹤效果對比分析。根據(jù)導(dǎo)彈的運(yùn)動特點(diǎn)，IMM算法的概率轉(zhuǎn)移矩陣采用如下設(shè)計(jì)

(19)

其他參數(shù)設(shè)置為：模型先驗(yàn)概率u0=[1/2, 1/2]T；CAF模型參數(shù)設(shè)置：軸向力加速度過程噪聲qa=10 m2·s-4，相對質(zhì)量燃燒率過程噪聲為qb=10-8，運(yùn)動加速度過程噪聲為qv1=1m2·s-4。PM模型的加速度過程噪聲為qv2=0.01 m2·s-4，空氣阻力參數(shù)噪聲為10-10。Singer模型算法參數(shù)設(shè)置：自相關(guān)時間常數(shù)為1/20，最大加速度為amax=100 m/s2，最大加速度概率為0.25，最小加速度概率為0.45。CS模型的自相關(guān)時間常數(shù)為1/20，最大加速度為amax=100 m·s-1，a-max=-50 m·s-2。

采用上述四種方法分別對仿真場景中的彈道導(dǎo)彈進(jìn)行跟蹤，并采用三點(diǎn)初始化方法[17]，其中，CAF模型中的軸向力加速度a的初值可以由初始狀態(tài)加速度的模值確定，相對質(zhì)量損耗率b的初始值設(shè)為0。PM模型以及文獻(xiàn)[3]方法的空氣阻力參數(shù)初值設(shè)為1/6 000。以均方根誤差為性能評價標(biāo)準(zhǔn)，Monte Carlo仿真實(shí)驗(yàn)100次。

圖5a)、圖5b)顯示的是本文方法和文獻(xiàn)[3]方法IMM算法單次實(shí)驗(yàn)時模型概率隨時間變化圖?？梢钥吹?，導(dǎo)彈在助推段飛行時，CAF模型起主導(dǎo)地位，模型概率接近1；而導(dǎo)彈關(guān)機(jī)后過渡到中段飛行，CAF模型概率迅速減小到接近0，而PM模型概率迅速增大到1附近，且一直保持不變。文獻(xiàn)[3]方法的模型概率變化與本文大致相同，但相對來說本文方法的模型概率更加合理，模型切換也更加迅速。圖中分別給出了兩種算法模型切換完畢的時間，可以計(jì)算出本文方法模型切換的過渡時間為6 s左右，而文獻(xiàn)[3]方法則需要15 s左右的時間。此外關(guān)機(jī)點(diǎn)檢測和估計(jì)具有重要意義，而通過本文IMM算法的模型概率的切換也可以用來檢測關(guān)機(jī)點(diǎn)時刻，即模型占優(yōu)切換前助推段CAF模型概率減小而中段PM模型概率同時增大的時刻可判為關(guān)機(jī)點(diǎn)時刻。

圖5 模型概率變化

接下來分析各算法的跟蹤誤差。圖6a)、圖6b)分別顯示了這四種算法的位置和速度均方根誤差曲線，表1給出了各算法觀測時間內(nèi)的平均估計(jì)誤差。表1表明：本文方法的平均誤差最小。從圖6中可以得出如下結(jié)論：

(1)助推段由于目標(biāo)距雷達(dá)較遠(yuǎn)，且機(jī)動性較強(qiáng)，四種算法的跟蹤精度都不是太高，但相對來說，本文方法的跟蹤精度遠(yuǎn)好于其他三種方法，模型與助推段匹配較好。

(2)導(dǎo)彈關(guān)機(jī)進(jìn)入中段飛行，本文方法和文獻(xiàn)[3]方法跟蹤誤差迅速減小，最后收斂到一個較小的值，且兩種方法在中段的跟蹤精度基本相當(dāng)。但由于本文方法采用精確的橢球模型對重力加速度建模，同時采用CKF作為跟蹤濾波器，跟蹤精度稍高于文獻(xiàn)[3]采用圓球模型對重力加速度建模，UKF作為濾波器的方法。單模型算法中，Singer模型算法和CS模型算法在中段的跟蹤誤差仍比較大，跟蹤誤差是多模型方法的幾倍。這是由于中段導(dǎo)彈目標(biāo)的機(jī)動性較弱，而仿真中Singer模型算法和CS模型算法的參數(shù)設(shè)置主要是針對助推段進(jìn)行設(shè)計(jì)的，保證在連續(xù)跟蹤中濾波器不會發(fā)散，而這導(dǎo)致過程噪聲協(xié)方差矩陣值較大，導(dǎo)致跟蹤精度低。這也是單模型算法的弊端所在，很難設(shè)置統(tǒng)一的參數(shù)，使其能適應(yīng)多種機(jī)動場景。另外CS模型算法在導(dǎo)彈關(guān)機(jī)處誤差增大明顯，容易導(dǎo)致濾波器發(fā)散甚至造成目標(biāo)丟失。

(3)對于助推段和中段的連續(xù)精確跟蹤來說，采用交互式多模型方法具有單模型方法不可比擬的優(yōu)勢，但綜合來說，本文方法具有更好的性能，尤其在助推段跟蹤精度遠(yuǎn)好于文獻(xiàn)[3]方法。

算法類型觀測時間內(nèi)估計(jì)平均誤差位置/m速度/(m·s-1)本文方法98．7518．36文獻(xiàn)[3]方法148．3736．22Singer模型295．84109．77CS模型224．2372．00

此外，本文算法還能對軸向力加速度大小進(jìn)行估計(jì)，圖7給出了單次實(shí)驗(yàn)的各時刻軸向力加速度大小的估計(jì)，可以看出，本文方法基本能準(zhǔn)確估計(jì)出軸向力加速度的大小。

圖7 軸向力加速度大小估計(jì)

3.2 仿真場景2

仍以仿真場景1中的目標(biāo)為跟蹤對象，雷達(dá)布站在導(dǎo)彈落點(diǎn)附近，跟蹤開始的時刻為彈道發(fā)射后的第600 s，跟蹤數(shù)據(jù)率為2 Hz，采樣點(diǎn)數(shù)為272。

由于再入段目標(biāo)狀態(tài)方程的非線性較強(qiáng)，此仿真場景下重點(diǎn)比較幾種非線性濾波算法的性能。記本文采用修正的無偏量測轉(zhuǎn)換，CKF作為濾波器的彈道導(dǎo)彈連續(xù)跟蹤的交互多模型算法為U-IMM-CKF算法，而采用EKF、UKF作為濾波器的算法記為U-IMM-EKF和U-IMM-UKF算法。同時為了分析采用無偏量測轉(zhuǎn)換和未采用無偏量測對濾波的影響，與不采用無偏量測轉(zhuǎn)換的被動段CKF(記為C-CKF)算法進(jìn)行了比較(由于不采用無偏量測轉(zhuǎn)換時，IMM算法性能很不穩(wěn)定，故直接與被動段模型進(jìn)行對比)。此外，與文獻(xiàn)[3]方法的跟蹤效果也進(jìn)行了對比分析。

采用上述五種方法對上述目標(biāo)進(jìn)行跟蹤，算法參數(shù)設(shè)置不變，濾波器初始化方法同仿真1，Monte Carlo仿真100次，仿真結(jié)果如圖8所示。

圖8的實(shí)驗(yàn)結(jié)果表明：

(1)當(dāng)空氣阻力較小時，系統(tǒng)非線性較弱，三種非線性濾波算法濾波性能基本相當(dāng)。而目標(biāo)重返大氣層進(jìn)入再入段，空氣阻力迅速增大，系統(tǒng)的非線性增強(qiáng)，EKF的濾波性能下降，濾波精度要差于CKF和UKF算法。CKF與UKF的濾波精度基本相當(dāng)，但UKF對參數(shù)選擇具有一定的依賴性，否則容易造成數(shù)值不穩(wěn)定問題。

(2)目標(biāo)距雷達(dá)較遠(yuǎn)時，不采用無偏量測轉(zhuǎn)換將產(chǎn)生較大的轉(zhuǎn)換誤差，造成跟蹤精度差。隨著目標(biāo)離雷達(dá)距離減小，轉(zhuǎn)換誤差減小，對濾波性能的影響減輕。

(3)文獻(xiàn)[3]采用彈道系數(shù)已知的六維再入段模型對目標(biāo)進(jìn)行跟蹤，可以看到，無目標(biāo)彈道系數(shù)的先驗(yàn)信息的情況下，此算法在再入段的跟蹤誤差急劇增大。而本文采用七維的再入段跟蹤模型，采用對目標(biāo)狀態(tài)和彈道系數(shù)進(jìn)行聯(lián)合估計(jì)的方法，即使無先驗(yàn)信息，也能實(shí)現(xiàn)對再入目標(biāo)的精確跟蹤。

圖8 均方根誤差

4 結(jié)束語

針對導(dǎo)彈連續(xù)跟蹤問題，提出了一種ECEF坐標(biāo)系中彈道導(dǎo)彈全階段連續(xù)跟蹤的交互式多模型算法。仿真結(jié)果表明，該算法相對其他算法具有良好的連續(xù)跟蹤性能，其采用的CKF濾波器具有較高的跟蹤精度和數(shù)值穩(wěn)定性，且通過IMM算法模型概率的變化能對導(dǎo)彈的關(guān)機(jī)點(diǎn)進(jìn)行檢測。此外，本文提出算法還能實(shí)現(xiàn)對目標(biāo)軸向力加速度和彈道系數(shù)的估計(jì)。當(dāng)然，本文算法仍存在不足，表現(xiàn)在估計(jì)誤差在關(guān)機(jī)點(diǎn)處有一定增大，這是下一步需要改進(jìn)的地方。

[1] LI X R, JILKOV V P. Survey of maneuvering target tracking, part II: motion models of ballistic and space targets[J]. IEEE Transactions on Aerospace and Electronic Systems, 2010, 46(1): 96-119.

[2] BENAVOLI A, CHISCI L, FARINA A. Tracking of a ballistic missile with a-priori information[J]. IEEE Transactions on Aerospace and Electronic Systems, 2007, 43(3): 1000-1016.

[3] 鈕俊清, 單 奇, 任清安, 等. 跟蹤彈道導(dǎo)彈全階段的可變多模型方法[J]. 雷達(dá)科學(xué)與技術(shù), 2011, 09(3): 224-231. NIU Junqing, SHAN Qi, REN Qingan, et al. A new algorithm with variable-structure multiple model for tracking ballistic missile at entire flight stages[J]. Radar Science and Technology, 2011, 09(3): 224-231.

[4] FARRELL W. Interacting multiple model filter for tactical ballistic missile tracking[J]. IEEE Transactions on Aerospace and Electronic Systems, 2008, 44(2): 418-426.

[5] LI X R, JILKOV V P. Survey of maneuvering target tracking, part V: Multiple-model methods[J]. IEEE Transactions on Aerospace and Electronic Systems, 2005, 41(4):1255-1321.

[6] 鮮 勇, 鄭曉龍. 彈道導(dǎo)彈攻防仿真系統(tǒng)建模[M]. 北京：國防工業(yè)出版社，2013. XIAN Yong, ZHENG Xiaolong. Modeling of ballistic missile attack defense simulation system[M]. Beijing: National Defense Industry Press, 2013.

[7] LI X R, JILKOV V P. A survey of maneuvering target tracking, part III: Measurement models[J]. Proceedings of SPIE Conference on Signal & Data Processing of Small Targets, 2001, 4473(11): 423-446.

[8] 趙艷麗, 李 宏, 高向東, 等. 不同坐標(biāo)系下中段彈道目標(biāo)跟蹤算法研究[J]. 現(xiàn)代雷達(dá), 2011, 33(5): 54-59. ZHAO Yanli，LI Hong，GAO Xiangdong, et al. A study on ballistic target tracking in midcourse of various coordinate systems[J]. Modern Radar, 2011, 33(5): 54-59.

[9] DUAN Z S, HAN C Z, LI R X. Comments on“Unbiased converted measurements for tracking”[J]. IEEE Transactions on Aerospace and Electronic Systems, 2004, 40(4):1374-1377.

[10] 李 為, 李一平, 封錫盛. 基于卡爾曼濾波預(yù)測的無偏量測轉(zhuǎn)換方法[J]. 控制與決策, 2015, 30(2): 229-234. LI Wei, LI Yiping, FENG Xisheng. Tracking with prediction-conditioned unbiased converted measurements[J]. Control and Decision, 2015, 30(2):229-234.

[11] ARASARATNAM I, HAYKIN S. Cubature Kalman filters[J]. IEEE Transactions on Automatic Control, 2009, 54(6): 1254-1269.

[12] BLACKMAN S S, POPOLI R. Design and analysis of modern tracking systems[M]. Boston: Artech House, 1999.

[13] GLASS J D, BLAIR W D, BAR-SHALOM Y. IMM estimators with unbiased mixing for tracking targets performing coordinated turns[C]// 2013 IEEE Aerospace Conference. [S.l.]: IEEE Press, 2013: 1-10.

[14] CARDILLO G P, MRSTIK A V, PLAMBECK T. A track filter for reentry objects with uncertain drag[J]. IEEE Transactions on Aerospace and Electronic Systems, 1999, 35(2): 394-409.

[15] 占榮輝, 郁春來, 萬建偉. 簡化UKF算法在單站無源目標(biāo)跟蹤中的應(yīng)用[J]. 現(xiàn)代雷達(dá), 2007, 29(3): 42-46. ZHAN Ronghui, YU Chunlai, WAN Jianwei. Aplication of smiplified UKF algorithm to single observer passive target tracking[J]. Modern Radar, 2007, 29(3): 42-46.

[16] 程水英, 張劍云. 粒子濾波評述[J]. 宇航學(xué)報, 2008, 29(4): 1099-1111. CHENG Shuiying, ZHANG Jianyun. Review on particle filters[J]. Journal of Astronautics, 2008, 29(4): 1099-1111.

[17] BAR-SHALOM Y, LI X R, KIRUBARAJAN T. Estimation with applications to tracking and navigation: theory algorithms and software[M]. New York: John Wiley & Sons, 2001.

許登榮 男，1991年生，碩士研究生。研究方向?yàn)槔走_(dá)數(shù)據(jù)處理。

程水英 男，1974年生，副教授。研究方向?yàn)闄C(jī)動目標(biāo)跟蹤、軟件無線電、信號與信息處理。

包守亮 男，1993年生，碩士研究生。研究方向?yàn)槔走_(dá)數(shù)據(jù)處理。

A Study on Ballistic Missile Tracking in ECEF Coordinate System

XU Dengrong，CHENG Shuiying，BAO Shouliang

(504 Laboratory of Electronic and Engineering Institute, Hefei 230037, China)

Aiming at the problem that single model tracking algorithm can not track the ballistic missile accurately, an interacting multiple model (IMM) algorithm is proposed in Earth-centered Earth-fixed coordinate for continuously tracking the missile for all the flight stage. According to the force condition in different stages of trajectory, the corresponding tracking model is established. The axial force of target almost keep constant in boost phase, the constant axis force model is adopted; coast and reentry phase use the passive phase tracking model to track the missile. Both the models use the cubature Kalman filter algorithm to filter. In view of the problem that a bias is caused in the mixing stage of this IMM algorithm, an unbiased mixing method is adopted. The simulation results verify the superiority of the algorithm.

maneuvering target tracking; ballistic missile; interacting multiple model algorithm; cubature Kalman filter

10.16592/ j.cnki.1004-7859.2016.11.012

國家自然科學(xué)基金資助項(xiàng)目(61201379)；安徽省自然科學(xué)基金資助項(xiàng)目(1608085MF123)

許登榮 Email：xudengrong@126.com

2016-08-15

2016-10-17

TN957.51

A

1004-7859(2016)11-0053-08