閃爍噪聲條件下基于交互多?？蚣艿睦走_(dá)目標(biāo)跟蹤算法*

占榮輝,李祖檢,滕書華

(1. 國防科技大學(xué) 電子科學(xué)學(xué)院, 湖南 長沙 410073; 2. 湖南第一師范學(xué)院 電子信息學(xué)院, 湖南 長沙 410205)

通常情況下,雷達(dá)傳感器的測量誤差主要受系統(tǒng)噪聲和處理方法的影響,相應(yīng)的噪聲可稱為“常規(guī)噪聲”。但當(dāng)雷達(dá)對復(fù)雜運(yùn)動目標(biāo)進(jìn)行探測時,由于目標(biāo)與傳感器之間的幾何關(guān)系將動態(tài)發(fā)生變化,目標(biāo)上不同散射中心的相對相位也會隨機(jī)發(fā)生變化,從而造成回波相位波前的畸變。波前在觀測方向上的傾斜和隨機(jī)擺動現(xiàn)象被稱為閃爍效應(yīng),由此導(dǎo)致的測量噪聲常被稱為“閃爍噪聲”。理論上,凡尺度能與雷達(dá)波長相比擬,具有兩個或兩個以上散射中心的任何體目標(biāo)都可能產(chǎn)生閃爍效應(yīng),且目標(biāo)與傳感器越近,因閃爍效應(yīng)產(chǎn)生的量測噪聲也就越大。這對諸如導(dǎo)彈制導(dǎo)等應(yīng)用場合極為不利,將嚴(yán)重影響制導(dǎo)控制精度,甚至?xí)饘?dǎo)彈脫靶。

與常規(guī)噪聲的高斯分布特性不同,閃爍噪聲具有明顯的拖尾特性,是典型的非高斯噪聲。目前,關(guān)于閃爍噪聲的分布主要有高斯-拉普拉斯混合模型[1]、混合高斯模型[2-4]、t分布模型[5-7]、高斯-t混合分布模型[8-9]等。文獻(xiàn)[10]采用QQ-plot方法對閃爍噪聲進(jìn)行擬合,得到了模型參數(shù)的估計(jì)結(jié)果。值得一提的是,雷達(dá)的測量結(jié)果是在極坐標(biāo)(或球坐標(biāo))下得到的,而觀測方程具有非線性特點(diǎn),因此閃爍噪聲條件下的目標(biāo)跟蹤問題是典型的非線性、非高斯問題。

眾所周知,卡爾曼濾波器(Kalman filter, KF)是線性、高斯假設(shè)條件下的最優(yōu)狀態(tài)估計(jì)器,在目標(biāo)跟蹤領(lǐng)域得到廣泛的應(yīng)用。在非線性、非高斯條件下,最優(yōu)化濾波假設(shè)不再成立,國內(nèi)外學(xué)者對此展開了廣泛的研究。在現(xiàn)有針對閃爍噪聲的跟蹤算法中,根據(jù)原理框架的不同,可粗略分為以下三類:

第一類是以非線性評分函數(shù)(score function)為基礎(chǔ)的跟蹤算法。在文獻(xiàn)[11]中,作者最早通過引入評分函數(shù)對KF進(jìn)行修正。作為一種早期提出的方法,該算法涉及煩瑣的混合分布函數(shù)卷積操作,在系統(tǒng)狀態(tài)維數(shù)較高時異常復(fù)雜,因此在實(shí)際中應(yīng)用較少。為此,文獻(xiàn)[1]在評分函數(shù)的基礎(chǔ)上,通過將量測噪聲分布進(jìn)行正交展開來避免卷積操作,然而該算法需要對量測向量進(jìn)行解耦合,同時也存在復(fù)雜的矩生成函數(shù)(moment generating function, MGF)計(jì)算問題,因此實(shí)用性也較差。

第二類是以高斯混合(Gaussian mixture, GM)或高斯和(Gaussian sum, GS)濾波為基礎(chǔ)的跟蹤算法[12-13]。這類算法將目標(biāo)狀態(tài)、過程噪聲、測量噪聲分別用多個高斯分量來近似,并采用多個并行的高斯近似濾波器完成狀態(tài)估計(jì)。高斯近似濾波器的具體實(shí)現(xiàn)形式有擴(kuò)展卡爾曼濾波器(extended KF, EKF)、不敏卡爾曼濾波器(unscented KF, UKF)和容積卡爾曼濾波器(cubature KF, CKF)等[14-15]。該類算法的優(yōu)點(diǎn)是實(shí)現(xiàn)結(jié)構(gòu)比較簡單,不足之處是:高斯分量數(shù)會隨遞推步數(shù)(時間周期數(shù))呈指數(shù)增長,因此在每步遞推完成后,需要進(jìn)行高斯分量剪枝和合并處理;對高斯分量數(shù)的確定缺乏嚴(yán)格的理論指導(dǎo),只能通過人工經(jīng)驗(yàn)進(jìn)行選擇優(yōu)化。

第三類是以粒子濾波(particle filter, PF)為基礎(chǔ)的跟蹤算法[16-18]。這類算法主要利用了粒子濾波器對非線性、非高斯系統(tǒng)具有強(qiáng)大的近似能力這一特點(diǎn),采用概率分布近似(而不是函數(shù)近似)手段進(jìn)行目標(biāo)狀態(tài)估計(jì)。粒子的采樣通常有兩種方式:一是直接通過先驗(yàn)分布(狀態(tài)轉(zhuǎn)移函數(shù))進(jìn)行采樣,這種采樣方式比較簡單,但粒子質(zhì)量不高;二是從融合了量測信息的重要采樣函數(shù)或建議分布[19-21]中進(jìn)行采樣,該采樣方式雖能在一定程度上改善粒子質(zhì)量,但每一個粒子都需一個類卡爾曼濾波器(卡爾曼濾波器及各種衍生形式)來實(shí)現(xiàn)完整的狀態(tài)遞推,導(dǎo)致算法復(fù)雜度急劇上升。除此之外,對于高維目標(biāo)狀態(tài)(相對一維標(biāo)量狀態(tài)),若要達(dá)到理想的密度近似效果和良好的跟蹤性能,粒子規(guī)模必須足夠大,相應(yīng)地,算法的復(fù)雜度也會非常高。

事實(shí)上,雖然閃爍噪聲被建模為各種混合分布的形式,但閃爍噪聲的發(fā)生本質(zhì)上是概率事件。若將出現(xiàn)“常規(guī)噪聲”和“閃爍噪聲”視為兩個不同的事件,對應(yīng)兩種不同的觀測誤差模型,同時將兩種噪聲出現(xiàn)的概率視為模型概率,并用一階馬爾可夫模型對相鄰時刻噪聲分布的變化特性進(jìn)行建模,則可利用交互多模(interacting multiple model, IMM)框架[22-23]進(jìn)行處理,解決傳統(tǒng)高斯混合濾波算法中對高斯分量選擇缺乏理論指導(dǎo)、存在模型適配性差的問題。

受此思路啟發(fā),本文以閃爍噪聲的高斯混合分布建模為基礎(chǔ),提出了一種交互多模框架下的高性能跟蹤濾波算法。該算法用一個高斯分量對應(yīng)一種噪聲分布模型,同時通過一階馬爾可夫模型建立轉(zhuǎn)移概率矩陣,并用免導(dǎo)數(shù)運(yùn)算、對非線性系統(tǒng)具有強(qiáng)適應(yīng)能力的容積卡爾曼濾波器與不同高斯分量相對應(yīng)的模型匹配濾波,最后對各匹配濾波結(jié)果進(jìn)行綜合得到最終的跟蹤結(jié)果。仿真結(jié)果證明了所提方法的有效性,并通過與傳統(tǒng)高斯混合濾波算法、粒子濾波算法等進(jìn)行比較說明了本文方法的性能優(yōu)勢和執(zhí)行效率優(yōu)勢。

1 問題描述

考慮典型的目標(biāo)跟蹤系統(tǒng),目標(biāo)離散時間狀態(tài)演化方程可描述為

xk=f(xk-1)+wk-1

(1)

其中:xk∈nx為nx維狀態(tài)矢量,通常包括目標(biāo)位置、速度等分量;f(·)為狀態(tài)轉(zhuǎn)移方程;wk-1為高斯過程噪聲,其均值為零,協(xié)方差矩陣為Qk-1,簡記為wk-1～N(wk-1;0,Qk-1)。

相應(yīng)地,目標(biāo)觀測方程可描述為

zk=h(xk)+vk

(2)

式中:zk∈nz為nz維觀測矢量,通常包含目標(biāo)距離、角度等分量;h(·)為觀測方程;vk為觀測噪聲。

通常情況下,vk被視為均值為零、協(xié)方差矩陣為Rk的高斯觀測噪聲。但在雷達(dá)受到干擾、復(fù)雜目標(biāo)姿態(tài)發(fā)生急變等條件下,測量誤差會出現(xiàn)跳變,稱為“閃爍”,測量噪聲的高斯分布特性不再滿足。此時,常用的處理方法是將閃爍噪聲建模為混合高斯分布,即

p(vk)=(1-ε)p1(vk)+εp2(vk)

(3)

2 混合高斯濾波算法

1)預(yù)測:假定k-1時刻狀態(tài)的后驗(yàn)PDF為p(xk-1|z1:k-1),則其一步前向預(yù)測的PDF可根據(jù)Chapman-Kolmogorov方程計(jì)算為

(4)

式中,p(xk|xk-1)稱為狀態(tài)的轉(zhuǎn)移密度函數(shù),用于描述一階馬爾可夫過程。

2)更新:在獲得p(xk|z1:k-1)的基礎(chǔ)上,綜合新近的觀測zk,可通過式(5)所示的貝葉斯規(guī)則計(jì)算k時刻狀態(tài)的后驗(yàn)PDF:

(5)

理論上,任何一個復(fù)雜的分布函數(shù)都可以用多個高斯分布求和的形式進(jìn)行近似,混合高斯濾波算法正是基于這一原理發(fā)展起來的。

(6)

又設(shè)k時刻的過程噪聲為包含J個分量的混合高斯分布,即

(7)

根據(jù)式(1),k時刻狀態(tài)轉(zhuǎn)移方程的先驗(yàn)分布為

(8)

由此可得k時刻更新(預(yù)測)狀態(tài)分布為

p(xk|z1:k-1)

(9)

同理,假定k時刻的觀測噪聲可建模為如下的混合高斯分布形式

(10)

由式(2)可知,k時刻觀測方程的似然分布為

(11)

由此可得k時刻測量更新(修正)狀態(tài)分布為

p(xk|z1:k)

(12)

由此可見,經(jīng)過一步完整的遞推過程,高斯分量已由k-1時刻的I個變?yōu)閗時刻的I·J·L個,隨著遞推的進(jìn)行,高斯分量數(shù)將爆炸式增長。為此需要在每一步遞推完成后,進(jìn)行剪枝處理,主要是舍棄權(quán)值較小的高斯分量,同時對相似性較強(qiáng)(空間距離小)的高斯分量進(jìn)行合并,具體方法可參考文獻(xiàn)[25]。經(jīng)剪枝后,重新恢復(fù)成I個高斯分量,即

(13)

3 交互多模框架下的CKF算法

3.1 閃爍噪聲等效建模

zk=h(xk)+vk(rk)

(14)

式中,rk∈。

由模型條件可知,Pr(rk=M1)=1-ε,Pr(rk=M2)=ε,其中Pr(·)表示概率函數(shù)。

假定用πij表示k-1時刻由模式Mi向k時刻模式Mj轉(zhuǎn)移的條件概率,則有

πij?Pr{rk=Mj|rk-1=Mi}i,j=1,2

(15)

相鄰時刻的模式變化關(guān)系可用一階馬爾可夫鏈來描述,其轉(zhuǎn)化關(guān)系如圖1所示。

圖1 常規(guī)噪聲和閃爍噪聲的模式轉(zhuǎn)移關(guān)系Fig.1 Mode transition relation of normal noise and glint noise

根據(jù)圖1中所示的模式轉(zhuǎn)化關(guān)系,可得轉(zhuǎn)移概率矩陣為

(16)

需要說明的是,式(16)與傳統(tǒng)為機(jī)動目標(biāo)跟蹤所設(shè)計(jì)的轉(zhuǎn)移概率矩陣是不同的,因?yàn)楹笳咴谥鲗蔷€上的元素體現(xiàn)的是下一個時刻保持與當(dāng)前時刻一致運(yùn)動模式的概率,因此其取值往往較大[26]。

對于每一個Mi,Mj∈,初始模型概率表示為

(17)

(18)

3.2 交互多模遞推框架

考慮由式(1)和式(14)組成的包含多種不同模型的混合系統(tǒng)估計(jì)問題,根據(jù)全概公式,系統(tǒng)狀態(tài)的后驗(yàn)PDF為

為表示方便,式中用rk(j)代表rk=Mj。

以各模型為條件的狀態(tài)后驗(yàn)PDF為

p(xk|rk(j),z1:k-1,zk)

(20)

將全概公式再次用于式(20)中的第二項(xiàng),則

p(xk|rk(j),z1:k-1)

(21)

在高斯近似假設(shè)條件下,式(21)可進(jìn)一步表示為

p(xk|rk(j),z1:k-1)

(22)

式中,μk-1[i|j]?Pr(rk-1(i)|rk(j),z1:k-1)稱為混合概率。

進(jìn)一步假定式(22)所示的混合形式為高斯混合分布,而后采用單個高斯分量對其進(jìn)行矩匹配近似,由此可得

p(xk|rk(j),z1:k-1)

(23)

利用式(23)實(shí)現(xiàn)的多模濾波算法稱為IMM算法,圖2給出了該算法的實(shí)現(xiàn)結(jié)構(gòu),包含2個并行的交互濾波器,混合過程在每個濾波器的輸入階段完成,并以z1:k-1為條件?；旌细怕实挠?jì)算將在下文的具體算法流程中詳細(xì)說明。

圖2 IMM算法框架Fig.2 Framework of IMM algorithm

3.3 基于CKF匹配濾波的算法實(shí)現(xiàn)

算法1 交互多模濾波Alg.1 Interacting multiple model filtering

算法2 基于CKF的模型匹配濾波Alg.2 Model-matched filtering based on CKF

需要指出的是,算法1中Step 6只用于結(jié)果輸出,并不參與遞推計(jì)算。

4 仿真實(shí)驗(yàn)

4.1 仿真條件

(33)

式中:

導(dǎo)彈的運(yùn)動方程為

(34)

雷達(dá)傳感器位于彈載觀測坐標(biāo)系的原點(diǎn),其觀測量包括目標(biāo)的距離和方位角,觀測方程的具體形式為

(35)

式中:距離和方位觀測噪聲vk的分布形式同式(3),閃爍噪聲出現(xiàn)的概率為ε。

4.2 仿真結(jié)果

仿真中假定目標(biāo)和導(dǎo)彈的初始狀態(tài)分別為

圖3 導(dǎo)彈攔截目標(biāo)實(shí)驗(yàn)場景Fig.3 Experimental scenario for target intercept

分別采用CKF、高斯混合容積卡爾曼濾波器GM-CKF、標(biāo)準(zhǔn)粒子濾波器(standard PF, SPF)以及本文提出的交互多模容積卡爾曼濾波器IMM-CKF對目標(biāo)進(jìn)行跟蹤,目標(biāo)初始狀態(tài)估計(jì)誤差的協(xié)方差矩陣為

P0=diag[(200 m)2(100 m/s)2(200 m)2(100 m/s)2]

單次實(shí)驗(yàn)跟蹤結(jié)果如圖4所示,由此可見,盡管四種算法都能實(shí)現(xiàn)對目標(biāo)的跟蹤,但其估計(jì)效果各有差異,這一點(diǎn)在圖4(b)可以更清楚地看出來;同時不難發(fā)現(xiàn),相較其他幾種算法,CKF的誤差最大,尤其是在閃爍噪聲出現(xiàn)的時刻更加明顯。

(a) 原圖(a) Original view

(b) 局部放大圖(b) Partially zoomed view圖4 單次實(shí)驗(yàn)跟蹤結(jié)果示例Fig.4 Target tracking for a single run

為了得到定量化的分析評估結(jié)果,對500次Monte Carlo實(shí)驗(yàn)結(jié)果進(jìn)行統(tǒng)計(jì),得到目標(biāo)位置估計(jì)的均方根誤差(root mean square error,RMSE)如圖5所示。

(a) x軸(a) Axis-x

(b) y軸(b) Axis-y圖5 目標(biāo)位置估計(jì)均方根誤差Fig.5 RMSE for target position estimation

由圖5中的結(jié)果可以清楚地看出,在存在較大初始誤差的情況下,四種跟蹤算法經(jīng)數(shù)個周期的遞推濾波以后,目標(biāo)位置估計(jì)誤差快速下降,并逐漸趨于收斂狀態(tài)。相比之下,盡管SPF算法在前幾個觀測周期以最快的速度使估計(jì)誤差降至一定的范圍,但在6 s之后誤差下降緩慢(見y軸)或已趨于不變(見x軸),最終的跟蹤精度相對偏低;CKF、GM-CKF和IMM-CKF三種算法的跟蹤誤差雖呈現(xiàn)一致的收斂趨勢,但其誤差大小各不相同,且漸次減小。整體上,IMM-CKF算法表現(xiàn)出最佳的跟蹤效果,這種優(yōu)勢在目標(biāo)跟蹤持續(xù)6 s之后更加明顯。

為考察不同閃爍噪聲對跟蹤性能的影響,在前述仿真條件不變的情況下,ε分別取0.05,0.10,0.15,0.20,0.25,0.40,重新開展實(shí)驗(yàn),并對6 s以后各觀測周期的均方根誤差取均值,得到平均均方根誤差(averaged RMSE,ARMSE),如表1所示。表中各個不同ε取值條件下第一行數(shù)據(jù)對應(yīng)x軸方向ARMSE,第二行數(shù)據(jù)對應(yīng)y軸方向ARMSE。

表1 不同閃爍噪聲條件下的位置估計(jì)ARMSETab.1 ARMSE for target position estimation with different glint probability 單位:m

由表1可知,隨著閃爍概率的增大,四種算法的ARMSE也隨之增大;相比之下,IMM-CKF算法的性能最優(yōu),且閃爍噪聲出現(xiàn)的概率越大,IMM-CKF相對于其他算法的性能增量越明顯。

由前文交互多模算法推導(dǎo)過程可知, IMM-CKF除了能跟蹤目標(biāo)狀態(tài),還可根據(jù)模型概率對閃爍噪聲出現(xiàn)與否進(jìn)行在線估計(jì)。圖6給出了ε=0.10時,單次仿真中在不同時刻出現(xiàn)真實(shí)閃爍噪聲的示例,以及由IMM-CKF得到的閃爍噪聲出現(xiàn)概率的估計(jì)結(jié)果。由圖6中的結(jié)果不難看出,對于大部分真實(shí)出現(xiàn)過閃爍噪聲的時間點(diǎn),IMM-CKF也相應(yīng)呈現(xiàn)較高(或接近1)的閃爍噪聲模型概率,說明估計(jì)結(jié)果與實(shí)際情況吻合良好。

圖6 閃爍噪聲出現(xiàn)時刻估計(jì)結(jié)果Fig.6 Estimation result for the existence of the glint noise

為定量評估IMM-CKF對閃爍噪聲出現(xiàn)時刻的估計(jì)性能,表2中給出了不同ε條件下對閃爍噪聲出現(xiàn)時刻估計(jì)的準(zhǔn)確率,其結(jié)果通過500次Monte Carlo實(shí)驗(yàn)得到。

表2 不同仿真條件下對閃爍噪聲出現(xiàn)時刻估計(jì)的準(zhǔn)確率Tab.2 Estimation accuracy for the existence of glint noise under different simulation conditions

由表2可以看出,總體上,所提IMM-CKF能以較高的準(zhǔn)確率對閃爍噪聲出現(xiàn)的時刻進(jìn)行估計(jì),這對目標(biāo)姿態(tài)突變檢測和干擾判斷應(yīng)用等具有重要的參考價值;同時可以看出,隨著閃爍概率的增大,對閃爍噪聲出現(xiàn)時刻估計(jì)的準(zhǔn)確率也有所下降,這與實(shí)際情況也是相符的。究其原因,主要是因?yàn)闉V波算法是通過遞推更新的形式實(shí)現(xiàn)的,即當(dāng)前時刻的模式概率估計(jì)結(jié)果同時受當(dāng)前時刻量測和前一時刻估計(jì)結(jié)果的影響。當(dāng)閃爍概率增大時,模型跳轉(zhuǎn)變得更為頻繁,導(dǎo)致短時間內(nèi)(或瞬時)得到穩(wěn)定的模型概率估計(jì)結(jié)果難度增大,從而造成閃爍噪聲出現(xiàn)時刻估計(jì)精度降低的現(xiàn)象。

表3 四種算法單次仿真運(yùn)行時間Tab.3 Running time of the four algorithms for a single run

由表3中的結(jié)果可以看出,由于CKF算法僅采用了單個高斯分量,因此耗時最少。GM-CKF和IMM-CKF的運(yùn)行時間分別是CKF的20.61倍和2.19倍,這與GM-CKF采用了20個高斯分量、IMM-CKF采用了2個高斯分量進(jìn)行近似的實(shí)際情況也是完全吻合的。相比之下,SPF耗時最長,且是在粒子并行采樣條件下得到的結(jié)果。由此可見,IMM-CKF算法復(fù)雜度適中,實(shí)時性強(qiáng),非常便于工程實(shí)現(xiàn)與應(yīng)用。

5 結(jié)論

閃爍噪聲是影響目標(biāo)跟蹤精度的一種重要因素,在制導(dǎo)雷達(dá)應(yīng)用中,如何克服閃爍噪聲的影響直接關(guān)系到導(dǎo)彈對目標(biāo)的命中精度和作戰(zhàn)效能。本文以閃爍噪聲的混合高斯建模為基礎(chǔ),將閃爍噪聲發(fā)生概率建成一階馬爾可夫模型,并在混合系統(tǒng)理論框架下導(dǎo)出了高斯近似解的表達(dá)式,同時通過容積點(diǎn)求解非線性積分方程,從而得到交互多模容積卡爾曼濾波器IMM-CKF。結(jié)合典型制導(dǎo)跟蹤示例,對所提算法的有效性進(jìn)行了驗(yàn)證。仿真結(jié)果表明,本文算法復(fù)雜度僅為傳統(tǒng)CKF的2倍,明顯低于GM-CKF和SPF算法,且在不同閃爍概率條件下均取得了一致最優(yōu)的跟蹤性能;除此之外,所提算法還能對閃爍噪聲出現(xiàn)時刻進(jìn)行有效的估計(jì),具有很好的實(shí)用價值。