基于序貫蒙特卡羅概率假設(shè)密度濾波的多目標檢測前跟蹤改進算法

占榮輝劉盛啟 歐建平 張 軍

(國防科學(xué)技術(shù)大學(xué)電子科學(xué)與工程學(xué)院 長沙 410073)

1 引言

在傳統(tǒng)的目標跟蹤中，通常先進行門限處理，從原始的觀測中提取目標的量測信息，這種處理方法常被稱為檢測后跟蹤(Detect-Before-Track,DBT)。DBT處理方法的潛在不足是，在進行門限處理過程中代表目標的有用信息可能被部分或全部舍棄，導(dǎo)致量測信息不完整(甚至有可能全部丟失)，這一問題在信噪比較低的條件下更為突出。為了提高低信噪比條件下的目標檢測與跟蹤能力，通常不事先設(shè)定門限，而直接對原始的觀測數(shù)據(jù)進行處理，同時完成目標的檢測與跟蹤，這種處理方法又稱為檢測前跟蹤(Track-Before-Detect, TBD)。

隨著技術(shù)的發(fā)展，已有多種TBD算法被相繼提出，這些算法按實現(xiàn)方式可粗略分為兩大類，即批處理方法和遞推型方法。作為批處理方法的典型代表如最大似然估計法[1]、動態(tài)規(guī)劃法[2,3]等需要利用多個連續(xù)幀的數(shù)據(jù)，其算法復(fù)雜度通常較大。遞推型方法(如貝葉斯推演)可有效降低運算開銷，且不需要處理和存儲多幀歷史數(shù)據(jù)，這些優(yōu)勢使其具有更強的適應(yīng)能力。貝葉斯推演的一種有效實現(xiàn)途徑是序貫蒙特卡羅(Sequential Monte Carlo, SMC)，也稱粒子濾波方法[4,5]，這種方法在單目標TBD背景下已有廣泛深入的研究。

文獻[6]給出了一種利用粒子濾波實現(xiàn)多目標TBD的典型應(yīng)用實例，文中通過設(shè)計轉(zhuǎn)移矩陣對目標個數(shù)的變化進行建模，并通過實驗證實該算法對具有兩個目標的特定應(yīng)用條件是非常有效的。其不足之處是很難推廣到更一般的多目標TBD應(yīng)用條件，且不能適用于最大目標數(shù)目事先未知的場合。

在目標數(shù)目未知且具有時變特性的應(yīng)用條件下，有限集統(tǒng)計量為多目標濾波提供了完整的貝葉斯推演手段，且這種推演是建立在隨機有限集(Random Finite Set, RFS)理論[7,8]基礎(chǔ)之上的。作為多目標貝葉斯濾波的一種矩近似求解手段，概率假設(shè)密度(Probability Hypothesis Density, PHD)濾波器[911]-自提出以來就受到持續(xù)的關(guān)注，且其SMC實現(xiàn)方式因具有理論上可證明的收斂特性[12]而被廣泛應(yīng)用。

在文獻[13]中，首次將PHD濾波思想引入到紅外圖像多目標TBD應(yīng)用中，通過建立目標運動模型和觀測模型，利用SMC方法實現(xiàn)對未知數(shù)目的多目標進行檢測與跟蹤，算法中完整融入了跟蹤的思想。不過，盡管文中的算法對非標準點目標模型具有適應(yīng)性，對目標個數(shù)的估計卻是有偏的。針對這一問題，文獻[14]對粒子的權(quán)值更新算式進行了修正，改善了估計效果。但是，作為一種新的理論方法，現(xiàn)有的基于SMC-PHD算法在TBD應(yīng)用條件下還存在以下3個突出問題：(1)用于近似目標狀態(tài)的有效粒子數(shù)少，效率低；(2)算法復(fù)雜度高，時間消耗大；(3)目標狀態(tài)提取不準確，穩(wěn)定性差。為此，本文在現(xiàn)有研究的基礎(chǔ)上提出了一種基于SMC-PHD的多目標TBD改進算法，該算法通過自適應(yīng)粒子產(chǎn)生機制使粒子在潛在的目標區(qū)域聚集，進而通過粒子子集分割手段加速運算，最后經(jīng)動態(tài)聚類方法實現(xiàn)多目標狀態(tài)的準確提取。仿真結(jié)果驗證了所提方法的優(yōu)越性。

2 系統(tǒng)模型

在多目標條件下，目標t的運動狀態(tài)可表示為

目標的強度可用傳感器的點擴散函數(shù)建模為

在高斯噪聲假設(shè)下，式(4)可進一步記為

由于目標的響應(yīng)僅能對其所處位置及其鄰近網(wǎng)格單元有較顯著的影響，因此似然函數(shù)可進一步近似為

3 基于SMC-PHD的TBD算法

由此可得預(yù)測粒子的權(quán)重為

進一步通過 PHD的更新方程可完成粒子權(quán)值的更新。注意到在量測中，僅那些緊鄰的分辨單元受到較大的影響，因此似然比可計算為

在文獻[13]最初給出的基于SMC-PHD的TBD算法中，提供了一種權(quán)值更新的經(jīng)驗算式。

式中

由式(12)～式(14)的計算式可以看出，文獻[13]在計算更新權(quán)時僅考慮了那些與待處理粒子具有相同支撐域的粒子，這對于傳統(tǒng)的點目標模型來說是合適；但對于具有擴散效應(yīng)的觀測模型，這種處理方法會導(dǎo)致目標數(shù)目的過估計。為了詳細說明這一點，假定由狀態(tài)引起的有效量測可表示為，顯然在TBD處理中，應(yīng)將整個而不是某一特定的分辨單元作為目標引起的量測。因此，對于某一特定的粒子而言，其似然比需用更準確的式子表示為

相應(yīng)地，式(13)中的ψ值應(yīng)修正為

得到k時刻的更新權(quán)值后，相應(yīng)的后驗密度可表示為

在此基礎(chǔ)上通過粒子重采樣和聚類處理，即可實現(xiàn)遞推濾波和目標狀態(tài)提取。

4 SMC-PHD TBD算法改進

4.1 自適應(yīng)粒子生成

在基于SMC的TBD算法中，一個重要的問題是如何有效地產(chǎn)生粒子來近似目標的真實狀態(tài)分布。通用的做法是在整個監(jiān)視區(qū)域均勻地散布粒子。這種方法的缺點是大量的粒子將被浪費在無效的近似上，因為這些粒子只占有微小的權(quán)重。另一種是數(shù)據(jù)驅(qū)動型方法，即采用量測數(shù)據(jù)來產(chǎn)生建議分布。在文獻[13]中，將粒子均勻分布在像平面的感興趣區(qū)域，感興趣區(qū)域定義為觀測強度大于某一預(yù)設(shè)門限Th的分辨單元，即。相比前面提及的通用做法，這種方法在某種程度上提高了粒子使用效率，但是量測信息仍沒有被充分利用。

事實上，為了更好地近似目標的狀態(tài)分布，人們總希望能自適應(yīng)地產(chǎn)生粒子，一種可行的方案是根據(jù)目標的存在概率按比例分配粒子。在這種思路啟發(fā)下，本文引入了一種新的代表新生目標的粒子產(chǎn)生策略，即首先通過強度量測信息對目標所在的分辨單元進行初始定位，而后通過粒子權(quán)值計算來評估目標的存在概率，最后通過重采樣實現(xiàn)粒子的自適應(yīng)分配。

具體地，假設(shè)將一部分觀測強度較大的分辨單元(代表了潛在目標可能的位置)記為，則在這些分辨單元附近均勻分配粒子，且位置、速度和強度分量分別為

并對其進行歸一化：

最后對粒子進行重采樣得

式中可采用不同的重采樣方法，如系統(tǒng)重采樣、殘差重采樣等。

經(jīng)重采樣后，代表新生目標的粒子將自動地向真實目標所處的位置聚集，從而能更好地近似多目標狀態(tài)的后驗密度。為了防止粒子過于集中可能導(dǎo)致的濾波性能下降問題，可在粒子重采樣之后添加隨機噪聲擾動來保持粒子的擴散特性。

4.2 SMC高效實現(xiàn)

通常，為了得到多目標狀態(tài)后驗密度的良好近似，所需的粒子數(shù)目較大，相應(yīng)地算法的運算開銷也會非常大。這是因為在式(16)中，計算每一個粒子的更新權(quán)時，都需要對粒子集里的所有粒子進行一次似然比計算。不過，由PHD濾波的收斂性不難推知，經(jīng)PHD預(yù)測和更新后，粒子將向真實目標所在的位置聚集。這就意味著，由粒子引起的有效量測(包含潛在目標的量測)數(shù)目不會太大，這是因為代表同一目標狀態(tài)的粒子將近似占有相同的圖像分辨單元。在這種假設(shè)前提下，通過對粒子全集進行子集分割處理，可達到TBD算法快速實現(xiàn)的目的。下面將給出所提方法的具體實現(xiàn)過程。

再對式(22)的位置分量進行離散化處理得到其在像平面中的分辨單元，即

最后通過對粒子進行集合歸并處理實現(xiàn)子集分割

SMC高效實現(xiàn)方法的關(guān)鍵步驟可用表1中的偽代碼表示。

表1 粒子權(quán)值更新的高效實現(xiàn)方法

4.3 目標狀態(tài)提取

在基于SMC-PHD的TBD算法中，需要從后驗密度中提取目標狀態(tài)，這主要通過粒子聚類來實現(xiàn)。常用的聚類方法為k-means算法[15]，即先通過更新后的粒子權(quán)值得到目標個數(shù)的估計為整數(shù))，再根據(jù)類間距離最小準則將全體粒子集聚成類，取每一類的聚心作為目標狀態(tài)的估計。這種方法的不足主要體現(xiàn)在兩方面：一是需要利用聚類數(shù)目先驗信息，目標個數(shù)估計不準確將直接導(dǎo)致狀態(tài)提取錯誤；二是需初始化聚心，聚心選擇不恰當(dāng)同樣將導(dǎo)致聚類錯誤。為此本文提出了一種動態(tài)聚類算法，該算法可不事先指定聚類數(shù)目，利用粒子的聚集特性及其在像平面中的分布規(guī)律實現(xiàn)自動聚類。

具體地，對于k時刻經(jīng)重采樣后的粒子，分別提取出其x軸、y軸坐標分量和，并按升序進行排列得到對應(yīng)的序列和，而后對其進行一階差分處理得到,。

顯然，對于同一個聚類，由于其擴散方差有限，因此將可容許的類內(nèi)間隔定義為和。若同時滿足，則將所有的粒子聚成一類；否則以和中超過門限,的元素個數(shù)作為各分量的類別數(shù)進行聚類。在此基礎(chǔ)上判斷原始粒子的坐標分量屬于哪個聚類，得到類屬矩陣。該矩陣中行向量不同的元素個數(shù)即為2維位置坐標總的聚類數(shù)，而行向量相同的元素序號(類屬矩陣的行號)所對應(yīng)的粒子即可組成一個獨立的聚類。

通常情況下，經(jīng)上述處理后即可完成粒子的準確聚類，實際應(yīng)用中為了進一步提高目標狀態(tài)提取的可靠性，可對每一個粒子聚類進行再次判斷，將粒子個數(shù)小于某一門限的聚類剔除。這是因為，經(jīng)重采樣之后，所有的粒子具有相等的權(quán)重，若某一聚類中的粒子個數(shù)過少，說明其權(quán)值之和較小，不能代表一個真實的目標，因此需將其視為虛假聚類進行剔除。

5 仿真實驗

本節(jié)將給出多目標TBD的應(yīng)用實例，仿真中假定目標運動滿足近勻速模型，而目標強度則采用隨機游走模型，其狀態(tài)方程可描述為

表2 基于SMC-PHD的TBD仿真參數(shù)

自適應(yīng)粒子采樣策略的有效性可通過圖1中的結(jié)果進行說明，圖中給出的是條件下的仿真示例。在圖 1(a)中，粒子被均勻散布在潛在目標所處的 20個網(wǎng)格分辨單元附近(擴散方差取為，粒子占據(jù)了大部分監(jiān)視區(qū)域，說明該粒子分配方法的效率是較低的。這是因為大量的粒子被分配到了并未出現(xiàn)真實目標的分辨單元，而代表真實目標狀態(tài)的粒子(有效粒子)數(shù)目極少，對目標狀態(tài)的近似效果不佳。相比之下，采用文中的自適應(yīng)粒子產(chǎn)生機制以后，絕大部分粒子聚集在真實目標所處的位置附近，如圖1(b)所示，該結(jié)果是先通過對粒子權(quán)值進行評估，進而采用重采樣技術(shù)實現(xiàn)粒子再分配得到的。由此可見，自適應(yīng)采樣方法存在兩大明顯的優(yōu)勢：(1)幾乎所有的粒子都用于近似真實的目標狀態(tài)，近似精度高；(2)粒子具有良好的聚集性，為后續(xù)基于SMC-PHD的TBD算法高效實現(xiàn)創(chuàng)造了有利條件。

為了進一步說明峰值像素點選擇方法的合理性(即不會使?jié)撛诘哪繕诉z漏)，仍采用k=1時的場景設(shè)置，動態(tài)調(diào)整峰值點的個數(shù)，考察10000次Monte Carlo仿真中所選擇的峰值像素點成功包含兩真實目標的次數(shù)與總的仿真次數(shù)的比率(這里視為目標的成功檢測率)，所得結(jié)果如圖2所示。由此可見，在本文給出的仿真條件下，當(dāng)峰值點個數(shù)取為 20時，成功檢測到目標的概率為 1，說明通過該方法來產(chǎn)生粒子不會導(dǎo)致目標漏檢。

圖3～圖5給出了目標初始強度I=30條件下的單次實驗結(jié)果，其中圖3為目標個數(shù)估計，圖4為不同方法提取的目標狀態(tài)，圖5為相應(yīng)條件下的最優(yōu)子模式分配(Optimal Sub-Pattern Assignment,OSPA[16])距離誤差。由圖4中的結(jié)果可以清楚地看出，即便是在目標個數(shù)估計完全正確的情況下，通過傳統(tǒng)k均值(k-means)聚類方法提取到的狀態(tài)仍存在部分錯誤，從而導(dǎo)致圖5中某些時刻的OSPA距離誤差顯著增大。究其原因，是因為k-means算法存在初始化環(huán)節(jié)，初始聚心選擇不當(dāng)將直接導(dǎo)致聚類錯誤。為詳細說明這一點，圖6進一步給出k=21時的單步仿真結(jié)果，其中虛線箭頭指向k-means算法得到的聚類結(jié)果，實線箭頭指向本文算法得到的聚類結(jié)果。由此可以看到，由于在隨機初始化過程中k-means算法將2個聚心錯誤地選擇在同一粒子云中，最終得到錯誤的聚類結(jié)果。而文中所提的聚類方法因利用了目標位置的有序差分信息，可有效克服這一問題，聚類結(jié)果完全正確(每一粒子云團得到一個與其完全對應(yīng)的、代表真實目標的聚類)。

圖7進一步給出了I=20條件下100次實驗得到的平均OSPA距離誤差。由圖可見，在新目標出現(xiàn)時刻，由于目標檢測可能存在一定的時延(在目標信號強度較弱的情況下這一問題將更加突出)，導(dǎo)致估計誤差出現(xiàn)短暫上升(這是時變多目標跟蹤中的一種普遍現(xiàn)象)。不過，文中所提算法的性能一致優(yōu)于傳統(tǒng)方法，這種優(yōu)勢的獲得主要源于兩方面原因：(1)當(dāng)兩種算法都能得到目標個數(shù)的準確估計時，所提算法因具有更好的目標狀態(tài)近似能力而呈現(xiàn)出更小的狀態(tài)估計誤差；(2)與傳統(tǒng)的k-means聚類方法相比，文中所提算法降低了錯誤提取目標狀態(tài)的可能性，從而降低了OSPA距離誤差。

圖1 基于峰值單元的新生粒子產(chǎn)生方法

圖2 不同條件下的目標檢測概率

圖3 目標個數(shù)估計結(jié)果

圖4 目標狀態(tài)提取結(jié)果

圖5 OSPA距離誤差比較

圖6 k=21時的目標狀態(tài)提取結(jié)果

圖7 平均OSPA距離比較

兩種不同實現(xiàn)方式下的算法運算開銷(CPU時間)如表3所示。表3中的傳統(tǒng)實現(xiàn)方式，指的是根據(jù)式(16)對所有的粒子逐一地計算概率生成泛函ψ，并通過式(12)進行粒子權(quán)值更新的計算方法。而文中所提的快速實現(xiàn)方式指的是先通過自適應(yīng)機制產(chǎn)生粒子，再利用 4.2節(jié)所述的子集分割策略進行ψ值計算和粒子權(quán)值更新的方法。仿真是在通用PC機中完成的，系統(tǒng)配置為2.0 GHz雙核CPU, 2.8 GB RAM，仿真軟件為Matlab7。由表中的結(jié)果可以看出，本文所提快速算法的執(zhí)行效率(單次仿真的平均耗時)比傳統(tǒng)實現(xiàn)方式要高得多，且在仿真參數(shù)設(shè)置條件下可以做到實時處理。運算效率的大幅提高主要是通過自適應(yīng)粒子產(chǎn)生機制使粒子聚集在某些特定的單元，在此基礎(chǔ)上根據(jù)圖像分辨單元將全體粒子集劃分成有限子集來處理而實現(xiàn)的。

表3 不同實現(xiàn)方式的時間消耗(s)

6 結(jié)束語

本文針對紅外傳感器目標探測應(yīng)用背景，提出了一種基于SMC-PHD的時變多目標TBD改進算法。該算法通過引入自適應(yīng)粒子產(chǎn)生機制改善了粒子的聚集程度及對新生目標狀態(tài)分布的近似精度，并利用粒子集分割技術(shù)提高了算法的執(zhí)行效率，最后通過粒子動態(tài)聚類方法減小了目標狀態(tài)的提取誤差。需要指出的是，文中所提的方法僅利用了圖像的幀內(nèi)信息，如何進一步挖掘圖像的幀間信息來提高多目標TBD性能是個值得探索的問題，也是下一步的研究方向。

[1] Tonissen S M and Bar-Shalom Y. Maximum likelihood trackbefore-detect with fluctuating target amplitude[J]. IEEE Transactions on Aerospace and Electronic Systems, 1998,34(3): 796-809.

[2] 鄭岱望, 王首勇, 楊軍, 等. 一種基于二階Markov目標狀態(tài)模型的多幀關(guān)聯(lián)動態(tài)規(guī)劃檢測前跟蹤算法[J]. 電子與信息學(xué)報,2012, 34(4): 885-890.

[3] Grossi E, Lops M, and Venturino L. A novel dynamic programming algorithm for track-before-detect in radar systems[J]. IEEE Transactions on Signal Processing, 2013,61(10): 2608-2619.

[4] Arulampalam S, Maskell S, Gordan N, et al.. A tutorial on particle filter for on-line nonlinear/non-Gaussian Bayesian tracking[J]. IEEE Transactions on Signal Processing, 2002,50(2): 174-188.

[5] Davey S J, Rutten M G, and Cheung B. Using phase to improve track-before-detect[J]. IEEE Transactions on Aerospace and Electronic Systems, 2012, 48(1): 832-849.

[6] Boers Y and Driessen J N. Multitarget particle filter track before detect application[J]. IEE Radar, Sonar and Navigation, 2004, 151(6): 351-357.

[7] 熊波, 甘露. MM-CBMeMBer濾波器跟蹤多機動目標[J]. 雷達學(xué)報, 2012, 1(3): 238-245.

[8] Mahler R. Multi-target Bayes filtering via first-order multi-target moments[J]. IEEE Transactions on Aerospace and Electronic Systems, 2003, 39(4): 1152-1178.

[9] Lin L K, Xu H, Sheng W D, et al.. Multi-target state-estimation technique for the particle probability hypothesis density filter[J]. Science China (Information Sciences), 2012, 55(10): 2318-2328.

[10] Nannuru S, Coates M, and Mahler R. Computationallytractable approximate PHD and CPHD filters for superpositional sensors[J]. IEEE Journal of Selected Topic in Signal Processing, 2013, 7(3): 410-420.

[11] Habtemariam B K, Tharmarasa R, and Kirubarajan T. PHD filter based track-before-detect for MIMO radars[J]. Signal Processing, 2012, 92(1): 667-678.

[12] Clark D and Bell J. Convergence results for the particle PHD filter[J]. IEEE Transactions on Signal Processing, 2006,54(7): 2652-2661.

[13] Punithakumar K, Kirubarajan T, and Sinha A. A sequential Monte Carlo probability hypothesis density algorithm for multitarget track-before-detect[J]. SPIE, 2005, 5913: 1-8.

[14] 林再平, 周一宇, 安瑋. 改進的概率假設(shè)密度濾波多目標檢測前跟蹤算法[J]. 紅外與毫米波學(xué)報, 2012, 31(5): 475-480.

[15] Clark D and Bell J. Multi-target state estimation and track continuity for the particle PHD filter[J]. IEEE Transactions on Aerospace and Electronic Systems, 2007, 43(4):1441-1453.

[16] Schuhmacher D, Vo B-T, and Vo B-N. A consistent metric for performance evaluation of multiobject filters[J]. IEEE Transactions on Signal Processing, 2008, 56(8): 3447-3457.