閃爍噪聲條件下基于交互多?？蚣艿睦走_目標跟蹤算法*

占榮輝,李祖檢,滕書華

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

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

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

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

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

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

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

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

受此思路啟發,本文以閃爍噪聲的高斯混合分布建模為基礎,提出了一種交互多?？蚣芟碌母咝阅芨櫈V波算法。該算法用一個高斯分量對應一種噪聲分布模型,同時通過一階馬爾可夫模型建立轉移概率矩陣,并用免導數運算、對非線性系統具有強適應能力的容積卡爾曼濾波器與不同高斯分量相對應的模型匹配濾波,最后對各匹配濾波結果進行綜合得到最終的跟蹤結果。仿真結果證明了所提方法的有效性,并通過與傳統高斯混合濾波算法、粒子濾波算法等進行比較說明了本文方法的性能優勢和執行效率優勢。

1 問題描述

考慮典型的目標跟蹤系統,目標離散時間狀態演化方程可描述為

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

(1)

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

相應地,目標觀測方程可描述為

zk=h(xk)+vk

(2)

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

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

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

(3)

2 混合高斯濾波算法

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

(4)

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

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

(5)

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

(6)

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

(7)

根據式(1),k時刻狀態轉移方程的先驗分布為

(8)

由此可得k時刻更新(預測)狀態分布為

p(xk|z1:k-1)

(9)

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

(10)

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

(11)

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

p(xk|z1:k)

(12)

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

(13)

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

3.1 閃爍噪聲等效建模

zk=h(xk)+vk(rk)

(14)

式中,rk∈。

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

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

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

(15)

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

圖1 常規噪聲和閃爍噪聲的模式轉移關系Fig.1 Mode transition relation of normal noise and glint noise

根據圖1中所示的模式轉化關系,可得轉移概率矩陣為

(16)

需要說明的是,式(16)與傳統為機動目標跟蹤所設計的轉移概率矩陣是不同的,因為后者在主對角線上的元素體現的是下一個時刻保持與當前時刻一致運動模式的概率,因此其取值往往較大[26]。

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

(17)

(18)

3.2 交互多模遞推框架

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

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

以各模型為條件的狀態后驗PDF為

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

(20)

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

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

(21)

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

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

(22)

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

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

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

(23)

利用式(23)實現的多模濾波算法稱為IMM算法,圖2給出了該算法的實現結構,包含2個并行的交互濾波器,混合過程在每個濾波器的輸入階段完成,并以z1:k-1為條件。混合概率的計算將在下文的具體算法流程中詳細說明。

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

3.3 基于CKF匹配濾波的算法實現

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

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

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

4 仿真實驗

4.1 仿真條件

(33)

式中:

導彈的運動方程為

(34)

雷達傳感器位于彈載觀測坐標系的原點,其觀測量包括目標的距離和方位角,觀測方程的具體形式為

(35)

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

4.2 仿真結果

仿真中假定目標和導彈的初始狀態分別為

圖3 導彈攔截目標實驗場景Fig.3 Experimental scenario for target intercept

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

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

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

(a) 原圖(a) Original view

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

為了得到定量化的分析評估結果,對500次Monte Carlo實驗結果進行統計,得到目標位置估計的均方根誤差(root mean square error,RMSE)如圖5所示。

(a) x軸(a) Axis-x

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

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

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

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

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

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

圖6 閃爍噪聲出現時刻估計結果Fig.6 Estimation result for the existence of the glint noise

為定量評估IMM-CKF對閃爍噪聲出現時刻的估計性能,表2中給出了不同ε條件下對閃爍噪聲出現時刻估計的準確率,其結果通過500次Monte Carlo實驗得到。

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

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

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

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

5 結論

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