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

占榮輝劉盛啟 歐建平 張 軍

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

1 引言

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

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

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

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

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

2 系統模型

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

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

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

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

3 基于SMC-PHD的TBD算法

由此可得預測粒子的權重為

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

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

式中

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

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

得到k時刻的更新權值后，相應的后驗密度可表示為

在此基礎上通過粒子重采樣和聚類處理，即可實現遞推濾波和目標狀態提取。

4 SMC-PHD TBD算法改進

4.1 自適應粒子生成

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

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

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

并對其進行歸一化：

最后對粒子進行重采樣得

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

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

4.2 SMC高效實現

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

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

最后通過對粒子進行集合歸并處理實現子集分割

SMC高效實現方法的關鍵步驟可用表1中的偽代碼表示。

表1 粒子權值更新的高效實現方法

4.3 目標狀態提取

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

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

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

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

5 仿真實驗

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

表2 基于SMC-PHD的TBD仿真參數

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

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

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

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

圖1 基于峰值單元的新生粒子產生方法

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

圖3 目標個數估計結果

圖4 目標狀態提取結果

圖5 OSPA距離誤差比較

圖6 k=21時的目標狀態提取結果

圖7 平均OSPA距離比較

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

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

6 結束語

本文針對紅外傳感器目標探測應用背景，提出了一種基于SMC-PHD的時變多目標TBD改進算法。該算法通過引入自適應粒子產生機制改善了粒子的聚集程度及對新生目標狀態分布的近似精度，并利用粒子集分割技術提高了算法的執行效率，最后通過粒子動態聚類方法減小了目標狀態的提取誤差。需要指出的是，文中所提的方法僅利用了圖像的幀內信息，如何進一步挖掘圖像的幀間信息來提高多目標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目標狀態模型的多幀關聯動態規劃檢測前跟蹤算法[J]. 電子與信息學報,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]. 雷達學報, 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] 林再平, 周一宇, 安瑋. 改進的概率假設密度濾波多目標檢測前跟蹤算法[J]. 紅外與毫米波學報, 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.