趙凱麗，高火濤，曹婷

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

0 引言

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

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

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

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

1 臨近空間高超聲速目標飛行軌跡研究

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

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

2 交互式多模型算法

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

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

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

假定在IMM算法中一共有m個模型集，各模型集對應的離散狀態(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) 濾波預測

(8)

(9)

(10)

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

(11)

(12)

(13)

(14)

(3) 模型概率更新

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

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

(15)

則模型j的后驗概率為

(16)

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

(17)

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

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

(18)

(19)

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

3 IMMCPF算法

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

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

(20)

(21)

時間更新：

(22)

(23)

(24)

(25)

量測更新：

(26)

(27)

(28)

(29)

(30)

(31)

(32)

(33)

更新粒子權值并歸一化：

(34)

(35)

(36)

(37)

4 仿真結果與分析

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

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

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

模型的狀態(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模型能自適應地調(diào)節(jié)加速度范圍，對加速度變化較大的目標有更好的跟蹤效果。

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

(42)

G(k)=

(43)

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

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

(44)

(45)

目標位置的均方根誤差定義為

(46)

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

圖4 各算法濾波結果對比圖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算法跟蹤得到的估計精度最高。這是因為臨近空間高超聲速目標的飛行軌跡具有很強的非線性，而IMMEKF算法僅在非線性較弱時效果較好；IMMPF算法對強非線性有較好的效果，但是容易產(chǎn)生粒子退化問題導致濾波精度下降，而IMMCPF算法的重要性密度函數(shù)中因加入了當前時刻最新的觀測信息來產(chǎn)生預測粒子集，這樣其生成的重要性密度函數(shù)與其最優(yōu)的重要性密度函數(shù)更加接近，所以IMMCPF算法的跟蹤結果精度最高。

表1為3種算法經(jīng)歷100次MC仿真后各狀態(tài)量的RMSE對比，從表1的數(shù)據(jù)中可以清楚的看出本文所使用算法對目標跟蹤精度的提升。相比于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算法相比具有更高的精度和更好的魯棒性，對臨近空間高超聲速目標的跟蹤也有更好的效果。

5 結束語

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