劃分交互式箱粒子概率假設密度濾波法

王 海，楊小軍

(長安大學 信息工程學院，陜西 西安 710064)

0 引 言

多目標跟蹤是民用自動駕駛、軍用雷達探測等方面的關鍵技術，其利用傳感器接收到的測量信息，通過預測更新依次確定多個目標的數量和狀態的過程。傳統的多目標跟蹤基于“小目標”假設[1]來處理復雜問題，常用方法有：聯合概率數據關聯算法(joint probabilistic data association，JPDA)[2]、多假設跟蹤算法(multiple hypotheses tracking，MHT)[3]、概率多重假設跟蹤(probabilistic multiple hypothesis tracking，PMHT)[4]、隨機有限集(random finite set，RFS)[5]等。RFS將目標和量測建立為隨機集集合來處理問題，文獻[6]提出概率假設密度(probability hypothesis density，PHD)濾波器，降低了多目標貝葉斯濾波器的組合計算困難問題。進一步提出了PHD濾波器的混合高斯(Gaussian mixture，GM)[7]實現，序貫蒙特卡羅(sequence Monte Carlo，SMC)方法[8-9]等。

在PHD濾波中，采用粒子對后驗狀態進行擬合，可以解決非線性、非高斯條件下的多目標跟蹤問題，但是為了達到良好的擬合效果，往往需要大量的粒子，這增加了計算量和時間。文獻[10-11]中，A. Gning等利用區間包含粒子濾波，提出箱粒子濾波算法，箱粒子可以有效減少粒子數，提高運算效率。文獻[12]提出了箱粒子與PHD結合的(Box-PHD)方法，驗證了箱粒子在計算效率上的優勢。進一步，文獻[13]運用箱粒子處理機動目標跟蹤問題，文獻[14]利用箱粒子劃分技術提升多目標狀態的估計性能。文獻[15]將交互式多模[16-17]的思想與箱粒子結合提出交互多模型箱粒子濾波(interacting multiple model box particle filter，IMM-BPF)來處理多機動目標。對于單傳感器的多機動目標跟蹤問題，該文在上述研究的基礎上，提出了劃分交互式概率假設密度濾波PIMM-Box-PHD算法，通過仿真實驗與IMM-Box-PHD算法進行比較，驗證了該濾波算法在跟蹤橢圓形多機動目標的良好性能。

1 系統模型

根據RFS理論，在多目標系統中，將所有測量時刻，所有目標的狀態和量測信息分別表示為兩個集合。由于狀態和量測信息的變化，可以用隨機有限子集來表示多目標狀態集Xk和量測集Zk。

Xk={xk,1,xk,2,…,xk,Nk}∈F(χ)

(1)

Nk表示k時刻目標數目，集合元素xk,i(i=1～Nk)表示單個目標i在單目標狀態空間F(χ)中的狀態矢量。

Zk={zk,1,zk,2,…,zk,Mk}∈F(ψ)

(2)

Mk表示k時刻量測數目，集合元素zk,j(j=1～Mk)表示量測j在單目標量測空間F(ψ)中的量測矢量，當然在量測信息中包含雜波量測。該文采用如下系統：

xk=fk-1(x)+uk-1

(3)

zk=hk(x)+vk

(4)

式(3)、式(4)分別表示k時刻的狀態方程和量測方程(線性或非線性)，過程噪聲uk-1和量測噪聲vk均為加性噪聲。在經典的貝葉斯框架下，將濾波過程分為預測和更新兩個部分，預測是根據馬爾可夫轉移矩陣，通過預測對象狀態來評估后驗密度，而更新是利用當前時刻的量測信息修正上一時刻的預測結果。

為了減少多目標貝葉斯遞推公式的計算復雜度，Mahler提出了一種次最優的解決方案─PHD濾波器。PHD濾波器通過傳播后驗強度函數完成目標數的估計，然后再根據預測目標數加權融合估計出每個目標的狀態，將雜波和不確定關聯中估計多個目標的問題被放到貝葉斯濾波框架中。PHD過濾器由兩個步驟組成：預測和更新，根據文獻[12]：

k時刻多目標后驗PHD預測：

其中，b(xk)表示新生目標強度，b(xk|xk-1)表示與前一狀態有關的衍生概率，s(xk-1)表示前一時刻到下一時刻目標存活概率，f(xk|xk-1)表示目標轉移概率。

k時刻多目標后驗PHD更新：

pk|k(xk)=[1-pD(xk)]pk|k-1(xk)+

(6)

其中，pD表示檢測概率，g(z|xk)表示測量似然，λ表示雜波均值，c(z)表示雜波分布先驗概率。上式分為兩部分，第一部分為未檢測到目標對更新PHD的影響，后一部分為檢測到的目標對其更新的貢獻。

2 箱粒子劃分

箱粒子濾波(Box-PF)在文獻[10]中被首次提出，是粒子濾波與區間分析相結合的一種濾波方法。在二維空間中利用區間形式畫出一定大小范圍的“箱”來表示點粒子，用箱粒子代替傳統的點粒子來進行濾波的預測和更新，稱為箱粒子濾波。相較于傳統的粒子濾波需要大量粒子來擬合后驗概率密度以得到良好的跟蹤效果，使用箱粒子可以有效提高計算效率，降低運算成本。由于區間測度的不確定性，在區間量測噪聲大雜波率高等復雜場景中，傳感器獲得的量測信息包含很多虛假量測，根據前一時刻量測信息產生的新生箱粒子偏大，用偏大且包含過多虛假量測信息的箱粒子去修正預測信息，會導致估計目標狀態有偏的情況[18-19]。

針對該問題，為了提高估計精度，提出箱粒子劃分的思想：該方法將預測后得到的持續箱粒子和新生箱粒子，在速度維和位置維進行箱的劃分，劃分之后的箱粒子變小，然后利用劃分后的小箱粒子集進行量測更新，通過預測步的劃分，去除冗余信息，找到包含更多量測信息的小箱粒子，就能在更新步更好地利用區間量測信息，達到提高目標估計精度的效果。

3 劃分交互式概率假設密度濾波算法

本節將箱劃分技術和IMM算法與BOX-PHD結合，提出處理多機動目標的PIMM-BOX-PHD算法。

使用IMM算法，首先需要確定目標的運動模型。理論上，建立的運動模型需要覆蓋目標的所有運動模式，對于各模型，產生一組箱粒子，這組箱粒子經過輸入交互、箱粒子濾波后進行重采樣，最后進行輸出交互。如此不斷循環遞推傳播更新這些箱粒子以完成對目標狀態的估計。

首先，設馬爾可夫模型轉移概率矩陣為：π={pij,i,j=1,2,…,Nr}，各模型初始概率為：U={μi,k,i=1,2,…,Nr}，Nr為模型數。該算法具體流程如下：

3.1 初始化箱粒子產生

3.2 交互式預測

計算k-1時刻模型i到j的混合概率:

(7)

得到模型j的交互箱粒子輸入：

(8)

目標狀態預測如下：

(9)

根據式(5)中的存活概率進行權值更新：

(10)

3.3 劃分箱粒子

在交互式預測之后，對預測的箱粒子進行劃分，可表示為：

3.4 更 新

(11)

其中：

(12)

λk|k-1([zm])=λc([zm])+

(13)

其中，[hcp]表示量測函數的包含函數，pD表示檢測概率，[gk|k(·)]表示箱粒子量測似然，λk表示雜波均值數，c(·)表示雜波分布先驗概率密度。因為量測是區間的形式，所以似然函數就不能采用傳統的目標集中式融合的方法，此處的似然函數通過量測箱粒子和預測箱粒子的重疊面積比，再采用收縮算法得到。

3.5 約束箱粒子

3.6 重采樣

3.7 模型概率更新

模型似然：

(14)

更新后的模型概率為：

(15)

3.8 目標狀態估計

(16)

最后將更新后的概率模型和提取的目標狀態混合輸出。

4 仿真實驗

4.1 仿真儀器

仿真場景中使用軟、硬件設施如下：

操作系統：Win10 64位操作系統；運行軟件：MATLAB R2012a；計算機品牌：Lenovo；處理器：Intel(R) Core(TM) i5-4258U CPU @ 2.40 GHz 2.10 GHz；RAM：8.00 GB。

4.2 仿真場景設置

假設一個雷達扇面觀測場景，傳感器雷達位于坐標[0,0]處，多個橢圓機動目標在[-200 m,200 m]×[0 m,200 m]的二維觀測區域內運動，運動模型由一個勻速直線運動模型、一個協同右轉運動模型和一個協同左轉運動模型構成。

勻速直線(CV)運動模型：

協同轉彎(CT)運動模型：

噪聲標準差為σv=0.02 m/s，三個模型之間的模型轉換概率矩陣為:

考慮到傳感器存在系統誤差，會造成一定的觀測偏差，把區間量測定義為：

在監控區域內檢測時長為100 s，目標存活概率ps=0.99，目標檢測概率pD=0.98，均勻雜波的泊松平均率λ=5，每個時刻產生新生箱粒子1個，存活箱粒子數8個。箱粒子擴展范圍為10 m。

4.3 實驗結果分析

假設雷達扇面中有三個機動目標，目標的初始狀態及出生死亡時刻如表1所示。

表1 目標初始狀態和時間

在雷達扇面檢測區域內，目標1在5 s～20 s做勻速直線運動，21 s～65 s協同右轉運動，66 s～85 s協同左轉運動；目標2在15 s～35 s做勻速直線運動，36 s～65 s協同左轉運動，66 s～95 s協同右轉運動；目標3在20 s～35 s做勻速直線運動，36 s～66 s協同左轉運動，67 s～86 s協同右轉運動。

圖1所示為三個橢圓形機動目標的運動狀態軌跡及其在單次蒙特卡洛實驗下的目標區間量測圖。

圖1 區間量測分布

該文采用最優子模型分配距離[22](optimal subpattern assignment，OSPA)對算法的估計精度進行評估。c為目標數目估計懲罰參數，c值越大，在OSPA距離中目標數目估計誤差所占比重越大，p為階數，也稱目標位置懲罰參數，p值越大，在OSPA距離中對目標位置估計誤差越敏感，位置估計誤差所占比重越大。本實驗中p=1，c=30。圖2所示為PIMM-Box-PHD算法單次蒙特卡洛實驗效果跟蹤圖，拆分為x和y兩個坐標，可以看出該方法能有效地對多機動目標的狀態和運動軌跡進行估計。

圖2 目標狀態估計跟蹤效果

表2 三種算法單次運行時間比較

從表2可看出，BOX-PHD方法平均運行時間最短，在加入IMM模型匹配之后運行時間變為原來的6倍，同時可以看出在加入箱粒子劃分技術之后運行時間比未劃分前增加1/6，由于PIMM-Box-PHD算法采用劃分技術使得箱粒子個數增多，計算效率略有下降。

進行100次蒙特卡洛實驗，圖3和圖4為IMM-Box-PHD與PIMM-Box-PHD兩種算法對多目標狀態估計的OSPA距離對比圖和OSPA位置分量對比圖。

圖3 OSPA距離估計

圖4 OSPA位置估計

圖3所示為綜合考慮目標個數估計和目標位置估計的OSPA距離，從圖中可以直觀看出PIMM-Box-PHD算法在所有采樣時刻，其OPSA距離都明顯小于IMM-Box-PHD，驗證了PIMM-Box-PHD在目標估計方面的優越性能。進一步，拆分出對多目標的位置估計分量，圖4所示為OSPA位置分量，可以看出在前80 s采樣時刻，PIMM-Box-PHD算法的OSPA距離明顯小于IMM-Box-PHD算法，前者對位置的估計更準確更穩定。仿真結果驗證了將箱粒子劃分為更小的箱粒子集，有助于目標的精確定位。

圖5 多目標的勢估計

圖5所示為多機動目標的勢估計，可以進一步看出，兩種算法的勢估計在采樣時間內與實際目標數進行比較，有新目標出生和原目標死亡時，對目標的勢估計存在一定的延時和誤差，不過很快修正，能夠準確檢測到目標個數的變化，兩種方法都能對目標的個數進行穩定的估計，在目標數目增加或減少時能夠迅速進行調整，記錄目標變化情況。PIMM-Box-PHD算法略有優勢，但效果不明顯。然后使用均值將兩種算法的OSPA距離進行量化比較，如表3所示。

表3 100次MC平均OSPA比較

通過表3可以看出在距離和位置估計方面，PIMM-Box-PHD算法的平均OSPA明顯小于IMM-Box-PHD，前者性能更優，說明劃分后的箱粒子能更好地利用量測信息，更準確地估計目標狀態和位置。在多目標個數估計方面，PIMM-Box-PHD算法的性能并不明顯，但其平均OSPA略小，說明箱劃分技術對目標個數的估計精度略有提升。

5 結束語

將箱劃分技術與IMM算法加入到BOX-PHD濾波器中，利用交互式實現多個模型對目標狀態進行聯合描述，實現對橢圓多機動目標跟蹤。通過仿真實驗對比分析，證明所提PIMM-Box-PHD算法在略微犧牲計算效率的情況下，能有效提高目標的跟蹤性能，對目標的定位更加精確，改進目標跟蹤過程中狀態中箱粒子冗余和估計有偏的情況，對最終多機動目標的狀態估計起到了很好的修正作用。