parfor 模式在粒子概率假設密度濾波中的應用

陳金廣，馬全海

(西安工程大學計算機科學學院，西安 710048)

針對多目標跟蹤，Mahler 在貝葉斯框架的基礎上提出了概率假設密度濾波(probability hypothesis density filter，PHD)，通過遞推多目標后驗的一階統計矩，從一階統計矩中提取目標的狀態［1］。目前，PHD 濾波主要有2 種實現方式:一種是基于粒子概率假設密度濾波(Particle PHD，P-PHD)(也叫序貫蒙特卡羅(Sequential Monte Carlo，SMC)的實現方式)［2］;另外一種方式是高斯混合概率假設密度濾波器(Gaussian Mixture PHD，GM-PHD)［3］。這2 種方式雖然都是由Vo 提出的，但是適用的情況卻不同。前者用諸多的粒子來近似PHD 濾波的一階統計矩，能適用于非線性非高斯系統;而后者僅適用于線性高斯系統。

基于粒子的PHD 濾波雖然有其自身的優點，但是隨著粒子數目的增多會使得計算量相應增大，從而導致時間復雜度也隨之增加。文獻［4］中提出一種非序貫重采樣算法，能夠實現P-PHD 濾波算法中權值更新步驟和重釆樣步驟的并行操作，在同樣的濾波效果下，降低了時間復雜度。文獻［5］中根據相對熵來動態計算每個目標所分配的粒子數目，而不是固定的給每個目標分配固定的粒子數，從而提高了濾波效率，減少了計算量。文獻［6］中則采用擬蒙特卡羅方法并行處理多目標跟蹤的數據關聯，使得運算效率得到提高。文獻［7］中間接地使用并行算法，將目標分為已有目標和新生目標兩類，利用極大似然門限來分配兩類目標對應的最優觀測數據，然后分別進行目標狀態更新，從而提高了濾波性能。文獻［8］中則是通過自適應門限濾掉雜波，然后僅利用落入門限內的量測進行更新，大大降低了時間復雜度。文獻［9］中闡述了粒子濾波具有可并行化的特點，并用OpenMP 實現了粒子濾波目標跟蹤的并行化處理，提高了粒子濾波的計算效率。文獻［10］中提出了獨立性分布式粒子濾波并行計算。并驗證了在處理時間上比粒子濾波更有優勢。文獻［11］中將粒子濾波和GPU 結合起來，充分發揮GPU(Graphics Processing Unit)的并行計算功能，在濾波效果相當的情況下，顯著減少了處理時延。文獻［12］中把并行算法引入到重采樣中，降低了重采樣的處理時延。

上述這些算法表明了并行計算在粒子濾波中具有可行性。文獻［4］中只是在理論上闡述，沒有實現粒子濾波的并行處理。文獻［11］中用OpenMP 語言實現了基于粒子濾波的圖像并行處理跟蹤，減少了處理時延。本研究把Matlab 中的parfor 并行計算結構引入到基于粒子的概率假設密度濾波的更新步驟中，實現了并行實驗。并行算法和傳統的PPHD 濾波算法相比在濾波精度上基本相同，隨著并行核數的增加，更新步驟的時間復雜度明顯減少。

1 背景知識

1.1 多目標隨機集模型

在基于隨機有限集(RFS)多目標濾波中，分為多目標狀態集和多目標量測集。k 時刻的狀態集和量測集可表示為:Xk={xk，1，…，xk，N}，Zk={ zk，1，…，zk，M}。其中M，N 對應表示量測目標個數和目標個數。也可表示為

式中:Kk表示所有雜波的量測集;表示所有真實目標所產生的量測集。

1.2 基于粒子的概率假設密度濾波(P-PHDF)算法

對應的，PHD 的粒子實現也對應為預測步驟和更新步驟。預測分為對新生目標的預測以及對存活目標的預測兩部分，需要預測的是對應粒子的狀態和權重。對于新生目標，其目標狀態和權值的預測方程為

對于存活目標，目標狀態和權重的預測由方程式(4)得到

式中:Lk-1、Jk分別為存活目標和新生目標對應粒子的數目;對應表示k -1 時刻和k 時刻第i 個粒子的狀態;為k 時刻的第i 粒子對應的權重;p(·|Zk)和q(·|Zk)分別對應于新生目標和存活目標的重要性函數為新生目標的概率假設密度函數。其中:ek|k-1(·|·)表示k-1 時刻的目標在k 時刻仍然存活的概率，fk|k-1(·|·)表示單個目標的狀態轉移概率密度，βk|k-1(·|·)表示k -1 時刻目標衍生的概率假設密度。

在基于粒子的PHD 濾波的更新步驟中，只需更新粒子的權重，粒子的狀態無需更新。權重更新方程為:

2 parfor 在P-PHD 中的應用

并行計算是一種許多指令同時進行的計算模式。它的思想是在同時進行的前提下，將計算的過程分解成小部分，之后以并發方式來加以解決。并行計算將一個復雜的任務分配給多個處理器去處理，而不是由一個處理器去順序執行，從而能夠減少處理的時延，大大提高運行的效率。Matlab自帶的工具箱可以在多核計算機上解決計算問題，其并行處理結構為任務并行(parfor 循環)和數據并行。①使用parfor循環執行任務并行時，Matlab 會將獨立的循環迭代自動分配給多個Matlab worker。工具箱中的parfor 循環提供了一種在多個Matlab worker 間分配任務的方式。使用該循環，可以將獨立的循環迭代自動分配給多個Matlab worker。parfor(并行for 循環)結構管理著Matlab 客戶端會話與worker 之間的數據和代碼傳輸。它會自動檢測是否有worker，如果沒有，則會還原為串行方式。②執行數據并行時，對于需要大型數據集處理的Matlab 算法，Parallel Computing Toolbox 提供了分布式數組、并行函數，以及使用spmd (單程序多數據)關鍵詞注釋代碼區段的功能，可用于在數個worker 上并行執行。這些并行結構可處理worker 間通信，并協調后臺的并行計算。

Matlab 的并行計算實質還是主從結構的分布式計算。當初始化Matlab 并行計算環境時，最初的Matlab 進程自動成為主節點，同時初始化多個(具體個數手動設定)Matlab 計算子節點。Parfor 的作用就是讓這些子節點同時運行Parfor語句段中的代碼。Parfor 運行之初，主節點會將Parfor 循環程序之外變量傳遞給計算子節點。子節點運算過程時互不干擾，運算完畢，則應該有相應代碼將各子節點得到的結果組合到同一個數組變量中，并返回到Matlab 主節點。在Matlab 并行工具箱中，已經提到了parfor 并行循環對于蒙特卡羅效果較好，本研究將parfor 并行循環引入到基于粒子的概率假設密度濾波中，利用其并行計算來減少運行時延。

在P-PHD 濾波框架中，其主要步驟分為初始化、預測、更新、估計目標數、重采樣和聚類等。在時間復雜度方面，重采樣和聚類較大，更新步驟次之，預測和其他步驟的時延基本可以忽略。傳統的P-PHD 濾波中，除更新步驟外，其他步驟暫不具備并行計算的條件，或者采用并行計算效果不明顯。理由如下:預測步驟中，分為存活目標的預測和新生目標的預測，對于存活目標可以引入parfor 并行計算，但是因為其占用的時延很小，故加速效果不明顯。重采樣的目的是為了去除權重小的粒子，復制權重大的粒子，在傳統的PPHD 的框架下不能用parfor 并行循環計算。聚類步驟是為了提取目標狀態，其方法本身限制了parfor 循環的使用。

由式(5)～式(7)得出，在P-PHD 濾波的串行更新步驟中，需要計算所有量測對于每個粒子權重的影響，由于每個粒子都是相互獨立的，這就滿足了并行計算的條件。在預測步驟中，已經得到了Lk=Lk-1+Jk個粒子的狀態和權重，即接下來是要對每個粒子權重的所有量測進行并行更新。對于Lk=Lk-1+Jk個粒子，使用任務并行計算parfor 循環并行更新粒子權重，步驟如下:①開啟并行parfor循環執行步驟②～⑨;②對于n 個量測值Zm×n，循環③～⑦來計算其對粒子權重的影響; ③根據公式θm= arctan計算粒子對應3 個傳感器的集中式融合量測值; ④根據③得出的融合量測值，計算似然g =N(θ1，Z(1，m)，σw(1))N(θ2，Z(2，m)，σw(2))N(θ3，Z(3，m)，σw(3));⑤由雜波平均數計算雜波強度kz=r/(2π)3;⑥ 由②中的似然 g 計算權值 w1= pdg/(kz+⑦累加所有量測對粒子權值的影響w2=w2+w1;⑧由公式得出更新后粒子權值;⑨最后得出

3 仿真結果與分析

假設在一個二維的場景中，監視區域范圍為［-40，60］×［-40，60］。采樣周期選為T =1，共仿真40 步。目標運動滿足線性模型，其運動方程為

式中，xk表示目標運動狀態，其表達式為xk=［x，vx，y，vy］T，其中(x，y)表示目標的位置信息，(vx，vy)表示目標的速度信息;過程噪聲并且雜波平均數為λ=3，服從泊松分布，并假設其均勻地分布在整個觀測空間。本實驗中設置3 個傳感器來跟蹤目標，其位置坐標分別為(-40，-40)，(10，60)，(60，-40)。假設傳感器只能獲得目標的方位角θk，并且對方位角進行集中式融合處理，觀測方程為

其中:x，y 分別對應目標的真實運動狀態的x 軸和y 軸坐標，si表示第i 個傳感器。量測噪聲diag(［0.000 1 0.000 1 0.000 1］)，vk，wk相互獨立。雜波服從均值為r 的泊松分布。設置目標檢測概率pD=1，無虛警，目標的存活概率pS=0.99。假設沒有目標衍生。表1 為4個目標的初始狀態和存活時刻。

表1 目標運動狀態

針對粒子濾波的更新步驟引入并行計算功能，并行更新每個粒子的權重，粒子數從500 ～5 000 分別開啟1 ～4 核做并行計算，其中1 核即為串行執行方式。各種情況的運行時間統計如表2 所示。綜合表2 可以看出，并行模式下，核數越多，實時性越好。這是因為在更新步驟中，對于每個粒子，都需要獨立計算其對于所有量測值更新之后的權重。并行計算可以并發的計算其權重，而不是順序執行，進而提高了實時性。需要指出的是，并行模式所用的核數與時間性能的提高并不是呈倍數關系，這是因為并行模式中，并行計算的同步和結果的匯總也需要消耗時間。

表2 并行更新運行時間

仿真環境為操作系統為win 7 旗艦版; 處理器為Intel(R)Core i5 -2400 CPU@3.10GHz; 內存為2G，32 位操作系統。Matlab 版本為R2013a。本算法采用Wasserstein 距離作為性能指標。以2 000 個粒子，開啟雙核并行，對并行算法和串行粒子概率假設密度濾波進行對比，仿真結果如圖1 ～圖4 所示。

圖1 給出了目標做線性運動時的真實運動軌跡，其中位于坐標軸上的3 個三角形表示傳感器的位置，標明了其視場范圍;圖2 則給出了真實坐標與所跟蹤的X 和Y 坐標，在第10、11、28 時刻略差于串行P-PHD 濾波，第7 時刻略優于串行P-PHD 濾波，其他時刻2 種算法所跟蹤的目標位置基本相同。圖3 給出了整個觀測時間內SMC-PHD 濾波算法以及并行算法對目標數的估計和實際目標數的對比，第7 時刻串行的P-PHD 濾波誤估計了目標數。從圖4 可以看出，在第7、29 時刻并行算法的估計誤差比傳統的P-PHD 濾波算法大，在第6 時刻較串行P-PHD 稍小，其他時刻基本相同。實驗結果可以表明并行算法和串行算法運算結果稍有差異，這是因為實驗過程中每次仿真所產生的數據具有隨機波動性，因而盡管并行算法滿足運算結果的封閉性，但是不同次的仿真結果是會出現差異的。

圖1 目標運動軌跡

圖2 估計目標X 和Y 坐標

圖3 目標數

圖4 Wasserstein 距離

4 結束語

針對傳統P-PHD 更新步驟中采用串行更新粒子的權重所存在的不足，通過引入并行計算更新粒子的權重。大量的實驗表明，在不失濾波性能的情況下，隨著核數的增加，更新步驟的運算時間得以減少，運算效率較串行計算有明顯提高。但在P-PHD 濾波中引入并行計算，要求數據必須是相互獨立的，否則不能實現并行計算。

［1］Mahler R. Multitarget Bayes filtering via first-order multitarget moments［J］.IEEE Transactions on Aerospace and Electronic Systems，2003，39(4):1152-1178.

［2］Vo B N，Singh S，Doucet A.Sequential Monte Carlo implementation of the PHD filter for multi-target tracking［C］//In Proceedings of the International Conference on Information Fusion.Cairns，Australia，2003:792-799.

［3］Vo B N，Ma W K.The Gaussian mixture probability hypothesis density filter［J］.IEEE Transactions on Signal Processing，2006，54(11):4091-4104.

［4］鄭云美.實時概率假設密度粒子濾波的算法研究［D］.浙江:浙江大學，2013.

［5］李威，韓崇昭，閆小喜.基于相對熵的概率假設密度濾波器序貫蒙特卡羅實現方式［J］. 控制與決策，2014，29(6):997-1002.

［6］張俊根，姬紅兵，蔡紹曉.并行高斯粒子濾波被動多目標跟蹤新算法［J］.西安電子科技大學學報:自然科學版，2011，38(1):117-122.

［7］章濤，來燃，吳仁彪.觀測最優分配的GM-PHD 多目標跟蹤算法［J］.信號處理，2014，30(12):1419-1426.

［8］章濤，吳仁彪. 自適應門限GM-PHD 多目標跟蹤算法［J］.數據采集與處理，2014，29(4):549-554.

［9］王愛俠，李晶皎，王青，等.基于多核的并行粒子濾波運動目標跟蹤［J］.計算機科學，2012，38(8):296-299.

［10］王美.粒子濾波并行處理的算法研究［D］.西安:西安電子科技大學，2014.

［11］孫偉平，向杰，陳加忠，等.基于GPU 的粒子濾波并行算法［J］.華中科技大學學報:自然科學版，2011，39(5):63-66.

［12］Gong P，Basciftci Y O，Ozguner F. A Parallel Resampling Algorithm for Particle Filtering on Shared-Memory Architectures［C］//Parallel and Distributed Processing Symposium Workshops ＆ PhD Forum (IPDPSW). Shanghai，China，2012:1477-1483.