趙凱麗，高火濤，曹婷

(武漢大學(xué) 電子信息學(xué)院，湖北 武漢 430072)

0 引言

臨近空間一般泛指距離海平面高度20～100 km的高空區(qū)域，臨近空間高超聲速飛行器因其機動特性強、軌跡變化快、強非線性等特點[1-2]，使這類飛行器很難被雷達跟蹤，有很好的隱身性，極難防御，如近幾年出現(xiàn)的X-51A[3]高超聲速飛行器， 對空天安全造成了巨大威脅。因此，臨近空間高超聲速目標(biāo)跟蹤相關(guān)算法的研究對空天安全的保障、促進空天防御技術(shù)發(fā)展都具有重大意義。

臨近空間目標(biāo)跟蹤對于普通機動目標(biāo)跟蹤中使用的單一模型跟蹤效果較差，難以覆蓋其可能采取的機動形式，而文獻[4-5]則證明多種目標(biāo)模型相互之間交互耦合使用的交互式多模型(interacting multiple model, IMM)算法對臨近空間目標(biāo)跟蹤有更好的效果。經(jīng)典IMM算法中使用的是Kalman濾波來對目標(biāo)狀態(tài)做濾波預(yù)測，但是在臨近空間中，環(huán)境惡劣，飛行軌跡受到各種擾動和空氣阻力的影響，這使得狀態(tài)方程和量測方程具有很高的非線性，而Kalman濾波算法僅對線性方程有較好的濾波效果，非線性方程濾波誤差會增加。

針對臨近空間的這一特點，國內(nèi)外已有學(xué)者相繼提出了一些基于IMM算法的改進算法，如秦雷使用擴展卡爾曼濾波(extend Kalman filter,EKF)[6]、Julier S使用無跡卡爾曼濾波(unscented Kalman filter,UKF)[7]、曹亞杰使用粒子濾波[8](particle filter,PF)分別來代替經(jīng)典IMM算法中的Kalman濾波來對目標(biāo)狀態(tài)進行估計。雖然這些改進算法都在一定程度上對目標(biāo)跟蹤精度有所提高，但是EKF算法僅在非線性較弱時效果較好，當(dāng)系統(tǒng)為強非線性環(huán)境時，狀態(tài)的估計精度和穩(wěn)定性都會有所影響，使目標(biāo)跟蹤誤差較大[9]。UKF算法雖然濾波精度高于EKF算法，但是它必須設(shè)置比例因子等相關(guān)參數(shù)的值，而這些參數(shù)值的選取目前并沒有確切的理論依據(jù)，并且這些參數(shù)的選擇對濾波的效果影響極大。PF算法的濾波精度要高于EKF或者UKF，在非線性系統(tǒng)中屬于濾波效果較好的算法，但是標(biāo)準(zhǔn)的PF算法在濾波過程中并沒有融入當(dāng)前時刻的量測信息，進而導(dǎo)致濾波精度下降甚至出現(xiàn)濾波發(fā)散現(xiàn)象，而引入了容積卡爾曼濾波的容積粒子濾波算法( cubature particle filter,CPF)，因加入了當(dāng)前時刻的觀測信息，使其重要性密度函數(shù)與后驗概率函數(shù)更加接近，從而使算法的濾波性能進一步提升[10-11]，目標(biāo)跟蹤效果更好。

本文將CPF算法與IMM算法相結(jié)合，形成交互式多模型容積粒子濾波算法(interacting multiple model cubature particle filter,IMMCPF)。IMMCPF算法使用IMM算法對目標(biāo)進行跟蹤，并使用CPF算法做IMM算法中的濾波預(yù)測，得到目標(biāo)狀態(tài)的估計值，即使用容積數(shù)值積分原則來計算非線性隨機變量的均值和協(xié)方差，這種方法既避免了求導(dǎo)運算，同時也加入了當(dāng)前時刻的觀測信息，使其重要性密度函數(shù)與其后驗概率函數(shù)更加接近，從而使算法的濾波性能得到提升。從上面的分析中可看出，IMMCPF算法能很好的實現(xiàn)對臨近空間高超聲速目標(biāo)的高精度跟蹤，仿真實驗也驗證了本文算法的可行性。

1 臨近空間高超聲速目標(biāo)飛行軌跡研究

根據(jù)美國公開的X-51A飛行器的部分?jǐn)?shù)據(jù)[12]分析可知，高超聲速飛行器的飛行過程大致如下：飛行器和助推器由轟炸機搭載，攜帶至約16 km的高空，并脫離飛行器，與此同時助推級開始點火，把飛行器推行到19 km的高空，飛行Ma約為4.5；然后助推級與飛行器脫離，飛行器在無動力情況下滑行，此時，超燃沖壓發(fā)動機開始點火，飛行Ma達到5，進入距離地面20 km的臨近空間，關(guān)閉發(fā)動機，進入無動力滑翔階段；最后發(fā)動機重新提供動力，將飛行器推送到新的高度。重復(fù)上述過程，一直到達預(yù)定的目的地上空，進入俯沖攻擊階段，攻擊目標(biāo)，完成作戰(zhàn)任務(wù)。

根據(jù)臨近空間高超聲速飛行器的飛行軌跡特點，本文近似擬合出了飛行器的整個運動過程大致如下圖1所示。圖2和圖3分別為飛行器的速度和加速度曲線。

圖1 飛行器軌跡Fig.1 Aircraft trajectory

圖2 飛行器速度變化曲線Fig.2 Velocity curve of the aircraft

圖3 飛行器加速度變化曲線Fig.3 Acceleration curve of the aircraft

從圖1的軌跡曲線看，飛行器飛行在上萬米的高空中，其軌跡變化復(fù)雜，且有很強的非線性，難以跟蹤；從圖2的速度變化曲線看，飛行器的平均飛行Ma能達到5～6，飛行速度極快；而從圖3中能看出，臨近空間高超聲速目標(biāo)的機動特性非常強。所以傳統(tǒng)的目標(biāo)跟蹤方法在臨近空間中跟蹤誤差較大，需要研究新的算法來對臨近空間高超聲速目標(biāo)進行跟蹤。

2 交互式多模型算法

根據(jù)上述仿真得到的臨近空間高超聲速目標(biāo)的飛行軌跡曲線，本文選擇使用對復(fù)雜的機動目標(biāo)跟蹤有較好效果的IMM算法來跟蹤臨近空間高超聲速目標(biāo)[13]。它采用多個單一模型來共同描述目標(biāo)的運動狀態(tài)，然后通過對不同模型分配相應(yīng)的權(quán)值來反映目標(biāo)運動狀態(tài)，最后對各模型加權(quán)融合即為所估計結(jié)果，因此該算法能很好的跟蹤臨近空間高超聲速目標(biāo)。

交互式多模型算法[14-15]是1988年Blom和Bar-Shalom在廣義偽貝葉斯算法的基礎(chǔ)上進一步提出的，該算法中有一個包括目標(biāo)整個運動狀態(tài)的模型集M，其中的任意一個模型都有一個對應(yīng)的濾波器與之相匹配，每個模型之間又可相互轉(zhuǎn)換，模型之間兩兩對應(yīng)的相互轉(zhuǎn)換與馬爾可夫過程一致，每個模型對應(yīng)的濾波器都并行運行，算法最終的結(jié)果是各個模型跟蹤結(jié)果的加權(quán)平均值[16-17]。

IMM算法的一個循環(huán)大致可分為4個部分：輸入交互、濾波預(yù)測、模型概率更新、數(shù)據(jù)融合[18]。

假定在IMM算法中一共有m個模型集，各模型集對應(yīng)的離散狀態(tài)方程可表示為

X(k+1)=FjX(k)+ωj(k)，j=1,2,…,m，

(1)

式中：Fj為模型j的狀態(tài)轉(zhuǎn)移矩陣；ωj(k)是均值為0，協(xié)方差為Qj的過程噪聲。

模型j的量測方程為

Z(k)=HjX(k)+vj(k)，

(2)

式中：Hj為模型j的觀測矩陣；vj(k)為觀測噪聲，其協(xié)方差值為R。

系統(tǒng)的馬爾可夫轉(zhuǎn)移矩陣為

(3)

式中：pij為從模型i到模型j的轉(zhuǎn)移概率。

(1) 輸入交互

μij(k-1k-1),

(4)

(5)

式中：

(6)

(7)

(2) 濾波預(yù)測

(8)

(9)

(10)

Sj(k)=HjPj(kk-1)HT+R，

(11)

(12)

(13)

(14)

(3) 模型概率更新

通常模型概率更新都會選擇使用最大似然函數(shù)法，根據(jù)該方法得到的似然函數(shù)值來計算各個模型與當(dāng)前時刻目標(biāo)運動模式的匹配值。

k時刻時，與第j個模型相匹配的最大似然函數(shù)可由式(15)計算得出：

(15)

則模型j的后驗概率為

(16)

式中：c為歸一化常數(shù)，

(17)

(4) 數(shù)據(jù)融合

在計算出各個模型的輸出概率和相應(yīng)的輸出結(jié)果后即可求得k時刻的交互式最終輸出結(jié)果：

(18)

(19)

以上即為IMM算法的一次完整的循環(huán)過程，當(dāng)進行下一時刻的交互循環(huán)時，就可把上面求得的輸出結(jié)果當(dāng)作下一時刻交互輸入的輸入值代入運算，計算新一輪的交互輸出。

3 IMMCPF算法

由于臨近空間高超聲速目標(biāo)的飛行軌跡具有很強的非線性特性，而交互式多模型算法的濾波預(yù)測中使用的是卡爾曼濾波，卡爾曼濾波對強非線性目標(biāo)的跟蹤效果不理想，有時甚至可能會出現(xiàn)濾波發(fā)散的現(xiàn)象。針對這個問題，本文選擇使用容積粒子濾波算法來代替Kalman濾波完成IMM算法的濾波預(yù)測過程，即提出一種新的IMMCPF算法來跟蹤臨近空間高超聲速目標(biāo)。

IMMCPF算法將IMM算法和CPF算法相結(jié)合，它使用一組權(quán)值大小一樣的容積點集來代替貝葉斯估計中復(fù)雜的積分運算，降低了算法復(fù)雜度，并且該算法在產(chǎn)生PF的重要性密度函數(shù)時進一步加入了當(dāng)前時刻最新觀測值，從而使其生成的重要性密度函數(shù)與系統(tǒng)真實的后驗概率密度函數(shù)更加接近。因此，IMMCPF算法在非線性系統(tǒng)中有更好的濾波性能，對臨近空間高超聲速目標(biāo)的跟蹤也有很好的效果。

(20)

(21)

時間更新：

(22)

(23)

(24)

(25)

量測更新：

(26)

(27)

(28)

(29)

(30)

(31)

(32)

(33)

更新粒子權(quán)值并歸一化：

(34)

(35)

(36)

(37)

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

本文擬根據(jù)X-51的試飛軌跡來建立臨近空間目標(biāo)仿真模型，并利用本文算法對臨近空間高超聲速目標(biāo)進行跟蹤，驗證所提算法性能。

根據(jù)臨近空間高超聲速目標(biāo)的飛行軌跡特征，本文選擇的交互式多模型算法中包含3個經(jīng)典的模型，2個Singer模型和1個Current模型。

(1) Singer模型：對一般的機動目標(biāo)(目標(biāo)機動不大時)有較好的跟蹤效果。

模型的狀態(tài)方程：

X(k+1)=FX(k)+ω(k)

，

(38)

式中：F為狀態(tài)轉(zhuǎn)移矩陣，

(39)

ω(k)為過程噪聲，協(xié)方差為

(40)

式中：

q11= 1/2α5(1-e-2αT+2αT+2α3T3/3-

2α2T2-4αTe-αT),

q12= 1/2α4(e-2αT+1-2e-αT+

2αTe-αT-2αT+α2T2),

q13=1/2α3(1-e-2αT-2αTe-αT),

q22=1/2α3(4e-αT-3-e-2αT+2αT),

q23=1/2α2(e-2αT+1-2e-αT),

q33=1/2α(1-e-2αT),

(41)

(2) Current模型：Current模型能自適應(yīng)地調(diào)節(jié)加速度范圍，對加速度變化較大的目標(biāo)有更好的跟蹤效果。

模型的狀態(tài)方程：

(42)

G(k)=

(43)

式中：T與α的定義與Singer模型相同。

本文仿真使用的3個模型的初始概率為(1/3,1/3,1/3)，3個模型之間的馬爾可夫轉(zhuǎn)移矩陣為

(44)

(45)

目標(biāo)位置的均方根誤差定義為

(46)

利用IMMEKF算法、IMMPF算法、IMMCPF算法仿真得到的目標(biāo)跟蹤軌跡與目標(biāo)真實運動軌跡對比圖如圖4所示。圖5是粒子數(shù)為100時各濾波算法在各方向上經(jīng)歷100次MC仿真后的目標(biāo)位置均方根誤差圖，圖6是相對應(yīng)的目標(biāo)速度均方根誤差圖。

圖4 各算法濾波結(jié)果對比圖Fig.4 Filter results comparison of eaeh algorithm

圖5 位置RMSEFig.5 Position RMSE

從圖5,6中可以看出，無論x方向位置RMSE、y方向位置RMSE，還是x方向速度RMSE、y方向速度RMSE，IMMEKF算法的估計精度都是最低的，IMMPF算法稍微高些，IMMCPF算法跟蹤得到的估計精度最高。這是因為臨近空間高超聲速目標(biāo)的飛行軌跡具有很強的非線性，而IMMEKF算法僅在非線性較弱時效果較好；IMMPF算法對強非線性有較好的效果，但是容易產(chǎn)生粒子退化問題導(dǎo)致濾波精度下降，而IMMCPF算法的重要性密度函數(shù)中因加入了當(dāng)前時刻最新的觀測信息來產(chǎn)生預(yù)測粒子集，這樣其生成的重要性密度函數(shù)與其最優(yōu)的重要性密度函數(shù)更加接近，所以IMMCPF算法的跟蹤結(jié)果精度最高。

表1為3種算法經(jīng)歷100次MC仿真后各狀態(tài)量的RMSE對比，從表1的數(shù)據(jù)中可以清楚的看出本文所使用算法對目標(biāo)跟蹤精度的提升。相比于IMMEKF算法，IMMCPF算法位置RMSE提升了13.12%，速度RMSE提升了28.21%；相比于IMMPF算法，IMMCPF算法位置RMSE提升了8.32%，速度RMSE提升了8.51%。

圖6 速度RMSEFig.6 Velocity RMSE

RMSE算法IMMEKFIMMPFIMMCPF位置RMSE/mx方向103.251 197.839 389.697 2y方向119.158 398.352 389.633 6速度RMSE/(m·s-1)x方向99.140 377.788 471.167 4y方向122.245 480.849 576.496 2

由此可以看出，本文提出的IMMCPF算法與IMMEKF算法和IMMPF算法相比具有更高的精度和更好的魯棒性，對臨近空間高超聲速目標(biāo)的跟蹤也有更好的效果。

5 結(jié)束語

本文根據(jù)臨近空間高超聲速目標(biāo)的飛行軌跡特性，選擇使用IMM算法來跟蹤這類目標(biāo)，但是經(jīng)典IMM算法中使用的卡爾曼濾波在非線性系統(tǒng)中濾波精度不高，針對這個問題，本文選擇使用CPF算法來代替IMM算法中的Kalman濾波算法，從而形成了一個新的IMMCPF算法。該算法使用CPF算法對目標(biāo)進行濾波預(yù)測，它融入了當(dāng)前時刻的觀測信息，使其重要性密度函數(shù)更加逼近后驗概率函數(shù)，從而使其濾波精度提高，對臨近空間高超聲速目標(biāo)的跟蹤有更好的效果。仿真結(jié)果也表明，IMMCPF算法各狀態(tài)量的估計精度均優(yōu)于IMMEKF算法、IMMPF算法。由此可看出，IMMCPF算法在臨近空間高超聲速目標(biāo)跟蹤中有更好的跟蹤效果，說明該算法有一定的實用價值。