基于JMS-SMC-PHD濾波的檢測前跟蹤算法

薛秋條,寧巧嬌,吳孫勇,2,蔡如華,伍雯雯

(1.桂林電子科技大學 數(shù)學與計算科學學院，廣西 桂林 541004;2.廣西密碼學與信息安全重點實驗室，廣西 桂林 541004)

0 引言

多機動目標跟蹤是目標通過在不同的模型之間的轉(zhuǎn)換，利用模型轉(zhuǎn)移概率矩陣來實現(xiàn)的一種方法，因為基于單一模型的跟蹤是不能準確預測目標的運動軌跡的[1]。跳躍馬爾可夫系統(tǒng)(jump Markov system，JMS)模型是解決機動目標跟蹤的有效方法[2-4]。JMS 通過使用多個描述目標運行軌跡的模型并行實現(xiàn)機動目標跟蹤。對于多目標的跟蹤，Mahler 利用隨機有限集理論(random finite sets, RFS)，提出了有限集統(tǒng)計學理論(finite set statistics, FISST)[5]，隨后提出了概率假設(shè)密度(probability hypothesis density，PHD)濾波[6]。PHD濾波是用一階統(tǒng)計矩近似后驗概率密度函數(shù)，該算法將復雜的多目標狀態(tài)空間運算轉(zhuǎn)化為單目標狀態(tài)空間運算，有效避免了復雜數(shù)據(jù)聯(lián)合問題。國內(nèi)外學者將PHD、cardinalized probability hypothesis density(CPHD)和Bernoulli 濾波使用點量測與JMS 模型結(jié)合，實現(xiàn)對多機動目標的跟蹤[7-10]。

檢測前跟蹤(track before detect，TBD)是解決低信噪比條件下弱小目標檢測與跟蹤的有效方法[11]。TBD 不用事先設(shè)定檢測門限對目標有無進行檢測，而是直接利用傳感器量測數(shù)據(jù)進行目標跟蹤，可以最大限度地保留目標的信息。近年來，國內(nèi)外學者將TBD技術(shù)與隨機有限集濾波進行結(jié)合，有效地實現(xiàn)對多個弱小目標的檢測與跟蹤問題[12-14]。對于機動弱小目標的檢測與跟蹤問題，國內(nèi)外學者采用多模的隨機有限集濾波與TBD 技術(shù)進行結(jié)合，也能有效地實現(xiàn)[15-16]。本文針對低信噪比條件下機動弱小目標的檢測與跟蹤問題，提出基于JMS-SMC-PHD 濾波的檢測前跟蹤算法。該算法以SMC-PHD 濾波為基礎(chǔ)，在目標跟蹤的預測階段利用跳躍馬爾可夫系統(tǒng)模型實現(xiàn)對目標狀態(tài)的轉(zhuǎn)移，在權(quán)值更新階段直接使用紅外傳感器量測數(shù)據(jù)實現(xiàn)對所有粒子的權(quán)值進行更新。最后通過仿真實驗證明，所提算法可以有效地實現(xiàn)多個弱小機動目標的檢測與跟蹤。

1 系統(tǒng)模型

1.1 目標運動模型和量測模型

針對多目標的跟蹤問題，假設(shè)在k(k=1,2,…)時刻，在觀測區(qū)域內(nèi)存在Nk個目標，并用一個隨機有限集表示為:

式中:E(χ)為目標狀態(tài)的所有有限子集。目標的狀態(tài)集包含目標的存活、新生和衍生部分，k時刻目標s(k=1,2,…,Nk)的運動方程可表示為:

紅外傳感器以一定的幀頻對監(jiān)視區(qū)域進行監(jiān)測，得到一幅關(guān)于監(jiān)測區(qū)域的二維圖像，圖像包括nXm的感應像元每一個分辨單元對應一個矩形區(qū)域ΔxXΔx，且每一分辨單元(i,j)的中心為(iΔx,jΔx)，i=1,…,n，j=1,…,m。則分辨單元(i,j)的觀測強度zk(i,j)可以表示為:

式中:ωk(i,j)為零均值高斯白噪聲，即并且各個單元之間是相互獨立的。C表示目標擴散影響區(qū)域其中為離目標狀態(tài)xks最近的整數(shù)單元。hk(i,j)(xks)為目標xks對分辨單元(i,j)的貢獻強度，一般采用點擴散函數(shù)的形式:

式中:∑為傳感器引入的模糊參數(shù)，一般是已知的。點擴散函數(shù)產(chǎn)生的強度大小以目標xks影響最大單元向周圍單元減小，所以目標的影響區(qū)域是有限的，在最靠近目標的單元取得最大值。從而k時刻的量測可以表示為:以及從1～k時刻的所有量測數(shù)據(jù)可以表示為:Z1:k={zl:l=1,…,k}。

1.2 跳躍馬爾可夫系統(tǒng)模型

跳躍馬爾可夫系統(tǒng)通過一組參數(shù)化的狀態(tài)空間模型去描述，各個狀態(tài)空間模型之間轉(zhuǎn)移是服從馬爾可夫鏈[5]。JMS 提供了一種不再是通過單一模型的方法去模擬機動目標，是解決機動目標跟蹤的有效方法。

k-1 時刻多目標狀態(tài)隨機有限集Xk-1中的目標增廣狀態(tài)(xk-1s,rk-1s)∈E(χ)=RnXN(s=1,…,Nk-1)要么以概率為Ps,k(xk-1s,rk-1s)存活到下一個時刻，要么以概率為1-Ps,k(xk-1s,rk-1s)死亡，其中xk-1s∈Rn和rk-1s∈N分別表示k-1 時刻運動狀態(tài)和模型標簽，N為所有模型標簽組成的離散集。模型轉(zhuǎn)移概率是服從離散馬爾可夫鏈，且存活目標狀態(tài)從k-1 時刻轉(zhuǎn)移到k時刻的轉(zhuǎn)移密度為:所以對于增廣狀態(tài)向量(xks,rks)的轉(zhuǎn)移概率可表示為:

紅外傳感器量測的似然函數(shù)可表示為:

式中:H1表示存在Nk個目標，H0表示不存在目標。pN(zk(i,j))表示分辨單元(i,j)處背景噪聲的概率密度函數(shù)，則表示目標加噪聲的概率密度函數(shù)。對于高斯背景噪聲條件下，上式可以表示為:

由于目標影響區(qū)域是有限僅能對其所在單元和鄰近區(qū)域有顯著影響，因此似然函數(shù)可進一步描述為:

2 基于JMS-SMC-PHD 濾波的TBD 算法

根據(jù)貝葉斯理論，濾波的實現(xiàn)可以分為預測階段和更新階段，預測是實現(xiàn)目標狀態(tài)的轉(zhuǎn)移，更新則是在預測的基礎(chǔ)上利用當前量測來做出修正。根據(jù)目標的狀態(tài)轉(zhuǎn)移服從一階馬爾科夫過程，所以預測和更新過程可以分別描述為:

式中:fk|k-1(Xk|Xk-1)是狀態(tài)轉(zhuǎn)移概率密度;gk(Zk|Xk)為似然函數(shù)。最后再利用后驗概率密度去估計目標的狀態(tài)。

2.1 概率假設(shè)密度濾波

PHD 濾波是貝葉斯濾波的一種智能近似，通過傳遞后驗概率密度的一階統(tǒng)計矩來代替后驗概率密度函數(shù)，可以有效避免數(shù)據(jù)聯(lián)合問題。

k時刻的預測密度和后驗密度分別用vk|k-1和vk表示，則PHD 濾波的預測可描述為:

式中:fk|k-1(?|x′)為已知目標狀態(tài)x′轉(zhuǎn)移到k時刻的轉(zhuǎn)移概率密度;pS,k|k-1(x′)為目標x′存活的概率;βk|k-1(?|x′)為目標x′在k時刻衍生RFS 的強度;γ(?)為新生RFS 的強度。

PHD的更新利用k時刻多目標的量測信息和預測的概率密度，所以后驗密度可描述為:

式中:pD,k(x)為k時刻的檢測概率;Zk為k時刻多目標的量測;κk(?)為雜波RFS 的密度。

2.2 JMS-SMC-PHD 濾波實現(xiàn)TBD

PHD 的序貫蒙特卡羅實現(xiàn)也稱PHD 的粒子濾波實現(xiàn)，是用一組帶權(quán)值的粒子去近似估計目標的后驗概率密度。跳躍馬爾可夫系統(tǒng)下的 PHDvk-1(xk-1,rk-1|Z1:k-1)用一組帶權(quán)值ωk-1(s)的粒子表示為:

式中:δ(?)為狄拉克函數(shù);rk-1為模型標簽;Lk-1是存活粒子數(shù)。然后通過預測和更新階段實現(xiàn)濾波:

1)預測:包含了新生粒子的產(chǎn)生和持續(xù)存活粒子的傳遞，新生粒子是由k-1 時刻的量測集zk-1={zk(i,j),i=1,…,n,j=1,…m}通過降序排列后，強度較大的前M個分辨單元所產(chǎn)生的Jk個新生的粒子。由于在強度較大的位置是真實目標所在位置的概率比較大，所以這樣產(chǎn)生的粒子更加接近真實目標。假設(shè)qk( ?|xk-1,Z1:k)和pk(?|Z1:k)分別為存活目標和新生目標的重要密度函數(shù)。則預測PHD 可近似為:

式中:

式中:

其中，衍生密度函數(shù)為:

馬爾可夫轉(zhuǎn)移概率密度為:

2)對上述得到的粒子根據(jù)量測進行粒子權(quán)值的更新，對于在量測zk中，僅有那些鄰近的分辨單元受到較大影響，因此似然比可計算為:

假定由狀態(tài)xk引起的有效量測為:從而對于狀態(tài)的權(quán)值更新為:

得到了更新權(quán)值后，后驗密度可表示為:

5)提取目標的狀態(tài):根據(jù)估計目標數(shù)zk利用k_means 方法對重采樣后的粒子進行聚類劃分，從而提取目標的狀態(tài)。

3 仿真分析

本章給出多機動目標TBD 仿真實驗，仿真中目標在3 個模型中相互切換，分別為直線勻速，左轉(zhuǎn)彎運動和右轉(zhuǎn)彎運動。目標的運動模型為:

vk為零均值高斯白噪聲，其協(xié)方差為:

式中:q1=2，q2=1 分別表示目標運動和強度的過程噪聲大小。模型r=1 時，為勻速直線運動模型，模型旋轉(zhuǎn)角度ω=0 rad/s;模型r=2 時，模型旋轉(zhuǎn)角度ω=10 rad/s;模型r=3 時，模型旋轉(zhuǎn)角度ω=-10 rad/s。馬爾可夫轉(zhuǎn)移概率矩陣和新生RFS 模型概率密度分別為:

傳感器以間隔T=1 s 接收30 幀圖像，其他參數(shù)值為:Δx=Δy=1，n=m=20，σ=1，∑=0.7。目標的存活概率為Ps=0.98。采用平面中M=20 個最大強度單元作為新生目標的潛在位置，速度均服從U[-0.4,0.4]的均勻分布，強度值服從U[I-1,I+1]的均勻分布，不考慮衍生情況即bk=[0;0;0]。采用最優(yōu)子模式分配度量(optimal subpattern assignment,OSPA)[17]來檢驗目標的跟蹤精度，OSPA 距是描述集合間的誤差的一種有效方法，其最大的優(yōu)勢是同時考慮了兩個集合間的勢誤差(檢測精度)和位置誤差(跟蹤精度)。OSPA 距的參數(shù)為:階數(shù)p=1，截斷誤差c=10。紅外傳感器強度量測的信噪比為:

式中:σ為量測圖中噪聲強度協(xié)方差。

將PHD 濾波用于多個弱小目標的檢測前跟蹤中有一個前提條件，即圖像中的每一個單元不能同時被多個目標影響，也就是說多個目標之間不能彼此靠近。假設(shè)共有 3 個目標，初始狀態(tài)為x1=[5,0.3,1,0.4,I]T，從1～25 s 時刻存活。x2=[19,-0.4,9,0.5,I]T，從3～30 s 時刻存活。x3=[5,0.3,18,-0.4,I]T，從6～27 s 時刻存活，其中I=15。圖1為真實目標的運行軌跡圖。圖2為目標位置狀態(tài)的估計圖。由圖2可知，估計的目標狀態(tài)都能近似的分布在真實目標位置附近。

使用Matlab R2013a 運行100 次蒙特卡羅實驗。圖3為目標狀態(tài)估計的平均OSPA 誤差。本文方法(JMS-PHD-TBD)與目標模型為單一模型的概率假設(shè)密度濾波的檢測前跟蹤方法(SING-PHD-TBD)作比較。由圖可知，由于在初始時刻，目標的初始狀態(tài)選取為空集的情況，所以算法在開始還不能估計出目標的狀態(tài)，而在目標出現(xiàn)的時候，估計都出現(xiàn)了一定的延時。但JMS-PHD-TBD 的估計誤差總體上要優(yōu)于SING-PHD-TBD。而單模型估計機動目標也能估計出目標的狀態(tài)是由于新生粒子時采用的是高強度位置產(chǎn)生粒子機制，這樣能保證真實目標在所采樣的粒子中的概率較高。并且由于弱小目標速度比較小，模型的轉(zhuǎn)換影響會稍微較小，但這也是不容忽視的。

圖1 目標真實軌跡Fig.1 Real trajectory of the target

圖2 目標位置狀態(tài)估計Fig.2 Target location status estimation

圖3 平均OSPA 誤差估計Fig.3 Mean OSPA error estimation

圖4為100 次蒙特卡羅實驗得到的平均目標數(shù)目的估計圖。由圖可知，在新目標出現(xiàn)時，目標數(shù)目的估計出現(xiàn)了延時，這也是導致圖3中估計的OSPA 誤差出現(xiàn)較大峰值的原因。JMS-PHD-TBD 方法對與目標的數(shù)目估計更加準確。出現(xiàn)目標跟蹤丟失的情況比較少，更加能準確地估計目標的狀態(tài)。

圖4 平均目標數(shù)目估計Fig.4 Average target cardinality estimation

4 結(jié)論

本文對跳躍馬爾可夫系統(tǒng)模型下的概率假設(shè)密度濾波的檢測前跟蹤方法進行深入的分析和研究。在序貫蒙特卡羅概率假設(shè)密度濾波的基礎(chǔ)上，融合跳躍馬爾可夫系統(tǒng)，用于解決多機動弱小目標的檢測與跟蹤。仿真實驗表明，通過跳躍馬爾可夫系統(tǒng)對模型進行狀態(tài)轉(zhuǎn)移并直接對單幀量測數(shù)據(jù)更新權(quán)值的處理，可以有效地實現(xiàn)多機動弱小目標的檢測與跟蹤。