parfor 模式在粒子概率假設(shè)密度濾波中的應(yīng)用

陳金廣，馬全海

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

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

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

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

1 背景知識

1.1 多目標(biāo)隨機集模型

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

式中:Kk表示所有雜波的量測集;表示所有真實目標(biāo)所產(chǎn)生的量測集。

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

對應(yīng)的，PHD 的粒子實現(xiàn)也對應(yīng)為預(yù)測步驟和更新步驟。預(yù)測分為對新生目標(biāo)的預(yù)測以及對存活目標(biāo)的預(yù)測兩部分，需要預(yù)測的是對應(yīng)粒子的狀態(tài)和權(quán)重。對于新生目標(biāo)，其目標(biāo)狀態(tài)和權(quán)值的預(yù)測方程為

對于存活目標(biāo)，目標(biāo)狀態(tài)和權(quán)重的預(yù)測由方程式(4)得到

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

在基于粒子的PHD 濾波的更新步驟中，只需更新粒子的權(quán)重，粒子的狀態(tài)無需更新。權(quán)重更新方程為:

2 parfor 在P-PHD 中的應(yīng)用

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

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

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

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

3 仿真結(jié)果與分析

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

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

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

表1 目標(biāo)運動狀態(tài)

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

表2 并行更新運行時間

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

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

圖1 目標(biāo)運動軌跡

圖2 估計目標(biāo)X 和Y 坐標(biāo)

圖3 目標(biāo)數(shù)

圖4 Wasserstein 距離

4 結(jié)束語

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

［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］鄭云美.實時概率假設(shè)密度粒子濾波的算法研究［D］.浙江:浙江大學(xué)，2013.

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

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

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

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

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

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

［11］孫偉平，向杰，陳加忠，等.基于GPU 的粒子濾波并行算法［J］.華中科技大學(xué)學(xué)報:自然科學(xué)版，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.