查巍巍 張 勇 應(yīng)曉慧 周 源 李知杰

（武漢船舶通信研究所 武漢 430064）

1 引言

粒子濾波是一種基于蒙特卡羅方法和遞推貝葉斯估計的統(tǒng)計方法［1］，適用于任何用到非線性、非高斯遞推貝葉斯估計的系統(tǒng)，精度可以逼近最優(yōu)估計。然而由于其計算量龐大，若采用集中式結(jié)構(gòu)，粒子信息計算和通信所產(chǎn)生的巨大能耗會導(dǎo)致中心節(jié)點的癱瘓，因此有必要探索分布式的處理方法。

現(xiàn)有的分布式粒子濾波（Distributed Particle Fifter，DPF）主要包括兩個方面：一種基于信息融合技術(shù)［2］即每個節(jié)點獨立地運行完整的粒子濾波形成一次狀態(tài)估計，然后交由融合中心進行融合；另一種基于并行計算，即將粒子濾波算法分解成可并行計算的進程，多個節(jié)點協(xié)作共同完成一次狀態(tài)估計。由于粒子濾波的計算量巨大，單個節(jié)點上獨立運行完整的粒子濾波算法將產(chǎn)生大量能耗，因此在硬件實現(xiàn)上多偏向于并行計算。圖1給出了分布式粒子濾波并行計算模型原理圖。

Bashi等設(shè)計了三種分布式粒子濾波算法［3］：全局DPF（Global DPF）、局部 DPF（Local DPF）和壓縮 DPF（Compressed DPF）。三種方法都是基于并行計算，不同主要在于重采樣，全局DPF和壓縮DPF采用的是全局重采樣，局部DPF采用的是局部重采樣，然而若采用全局重采樣則需要上傳大量粒子的信息，增加了通信損耗，若采用局部重采樣，會出現(xiàn)部分子節(jié)點由于粒子退化對全局狀態(tài)估計作用甚小，導(dǎo)致最后全局估計誤差偏大。

圖1 分布式粒子濾波并行計算原理圖

為了克服全局重采樣帶來的大通信量弊端，有人提出了基于GMM模型的DPF［4］，主要改進在于將傳輸?shù)牧Ｗ佑肎MM分布近似，然而每次狀態(tài)估計都需要計算GMM模型參數(shù)，會產(chǎn)生大量的計算能耗。因此，本文給出了一種精度能達到全局重采樣，而且不需要傳輸具體的粒子信息的分布式粒子濾波并行計算方法。其中重采樣算法來源于比例分配重采樣［5］，只是將該算法在粒子集上并行實現(xiàn)，并采用粒子群優(yōu)化算法解決部分子節(jié)點由于粒子退化帶來的計算量失衡的問題。仿真結(jié)果表明，本文的算法在保持負載均衡的同時減少了通信損耗，并且可以取得較好的估計精度。

2 粒子濾波算法

粒子濾波器的本質(zhì)是遞歸貝葉斯濾波，根據(jù)貝葉斯估計理論，一階馬 爾科夫序列｛xk｝k＝1，2，…的估 計問題是通過觀測序列｛zk｝k＝1，2，…和狀態(tài)轉(zhuǎn)移概率密度函數(shù)p（xk｜xk－1）得到后驗概率密度函數(shù)p（xk｜zk），即：

遞推過程可分為預(yù)測和更新兩步

預(yù)測：

更新：

根據(jù)蒙特卡羅原理，后驗概率分布可以用一系列帶權(quán)重的隨機采樣點逼近，即：

由于從后驗概率中很難直接采樣，所以引入容易采樣的已知分布的重要密度函數(shù)：q（x0：k｜z1：k）。粒子從重要密度函數(shù)中采樣得到，則權(quán)值為

一般選擇狀態(tài)轉(zhuǎn)移概率密度函數(shù)p（xk｜xk－1）作為重要密度函數(shù)。

粒子濾波算法由四步組成，描述如下：

1）采樣。從p（xk｜xk－1）中采樣N 個粒子

2）權(quán)值計算。

并且歸一化

則可得k時刻的狀態(tài)估計為

3）重采樣。得到新的粒子集合

4）時刻k＝k＋1，轉(zhuǎn)到第1）步。

3 分布式處理

分布式粒子濾波算法就是將式（6）～（8）在粒子集上并行計算，這三個式子可以很好的在粒子集上進行拆分計算，但是難點在于重采樣的并行處理，不同的處理思想直接影響到分布式粒子濾波的執(zhí)行效率和估計精度。

常見的重采樣算法有多項式重采樣算法［1］、殘差重采樣算法［6～7，9］、分層重采樣算法［6～8］、系統(tǒng)重采樣算法［6～9］。然而這些重采樣算法都建立在粒子權(quán)重的累積分布函數(shù)，需要對全部的權(quán)重信息集中排序，若將這些算法并行計算，則不可避免進行粒子權(quán)重信息的交換，將帶來額外的通信能耗。

Bolic提出比例分配重采樣算法，該算法按粒子歸一化的權(quán)重進行按比例采樣，建立在權(quán)重和之上，不需要交換每個粒子的權(quán)重信息，因此可以很好地進行并行計算。

接下來可以進行并行重采樣計算，偽代碼如下：

Step1：中心節(jié)點根據(jù)子節(jié)點的權(quán)重和進行第一次一次比例分配重采樣。

Step2：中心節(jié)點將 ｛u′ ，C｝變量傳給子節(jié)點，子節(jié)點對本地粒子進行第二次比例分配重采樣。

其中，中心節(jié)點產(chǎn)生的N＿resample（j）表示第j個節(jié)點上應(yīng)該重采樣N＿resample（j）個粒子，子節(jié)點上n＿resam－表示第j個子節(jié)點上第i個粒子應(yīng)該被重采樣n＿resa－次。而變量u′的計算剛好保證第j個子節(jié)點上通過step2得到的粒子數(shù)之和與中心節(jié)點計算的相同，即：

也就是說局部獨立重采樣的結(jié)果和全局重采樣的結(jié)果是一樣的，圖2為并行比例分配重采樣的計算原理圖。

圖2 并行比例分配重采樣原理圖

其中節(jié)點1和節(jié)點2的重采樣并行進行，U和U′是并行重采樣用到的隨機數(shù)，u為集中式重采樣用到的隨機數(shù)，n（i）為第i個粒子重采樣后重復(fù)的個數(shù)，若

則并行重采樣和集中式重采樣能得到相同的結(jié)果，即相同的n值。

根據(jù)Step1得

由集中式比例分配重采樣算法得

即U′（1）＝u（6），以此類推，式（11）得證。也就是說比例分配重采樣并行計算的關(guān)鍵在于每個子節(jié)點上第一個隨機數(shù)的產(chǎn)生，即變量u′的計算。

4 分布式處理帶來的負載失衡問題

通過第3節(jié)可以很好地將粒子濾波算法并行計算，但是在若干次估計之后，一些子節(jié)點由于粒子退化使得權(quán)重值和很小，中心節(jié)點比例分配重采樣之后使得這些節(jié)點只需要產(chǎn)生少量的粒子數(shù)，即N＿resample很小，極端情況下，N＿resample只有一個非零，也就是說粒子全部在一個節(jié)點上計算，分布式粒子濾波退化成集中式粒子濾波。

上述問題可以理解為并行計算系統(tǒng)中的負載失衡問題，即有些節(jié)點超載，有些節(jié)點輕載，必須采用一種有效的算法均衡各個負載之間的計算任務(wù)。本文采用粒子群優(yōu)化算法（PSO）對出現(xiàn)輕載的節(jié)點上的粒子進行優(yōu)化，使其總體權(quán)重和變大，那么在下一次狀態(tài)估計中，N＿resample變大，在一定程度上減輕負載失衡的問題。

粒子群優(yōu)化算法（PSO）采用兩個極值來更新粒子的速度和位置，使得整個粒子群向最優(yōu)位置靠近，一個為個體極值Pi，即粒子本身從初始態(tài)到當(dāng)前態(tài)經(jīng)歷的最優(yōu)位置，另一個為種群極值Pg，即整個粒子群經(jīng)歷過的最優(yōu)位置。利用粒子群優(yōu)化算法對采樣過程進行優(yōu)化，使得粒子向高似然概率區(qū)域運動。本文采用一種改進的PSO算法［10］，該方法基于一個高斯分布來不斷更新粒子的速度，其收斂性好于經(jīng)典的PSO算法。

一般狀態(tài)空間模型分為狀態(tài)轉(zhuǎn)移模型和量測模型：

其中wk和Vk分別為過程噪聲和量測噪聲，并且是相互獨立，協(xié)方差分別為Qk和Rk的零均值加性噪聲序列，因此真實值應(yīng)該在觀測值周圍服從均值為0，協(xié)方差Rk的噪聲分布，定義適應(yīng)度函數(shù)為

PSO更新過程為

經(jīng)過M次迭代后，粒子的最優(yōu)值符合某閾值ε時，說明粒子群已經(jīng)分布在真實狀態(tài)附近，則停止優(yōu)化。

圖3是一個用PSO算法解決負載失衡問題的示意，其中總粒子數(shù)為100，子節(jié)點數(shù)為5。其中大括號內(nèi)的值表示每個節(jié)點重采樣之后的粒子數(shù)。

圖3 利用PSO算法解決負載失衡的例子

引入PSO算法的目的是為了優(yōu)化節(jié)點上的粒子，在一定程度上解決負載失衡的問題，但同時帶來了額外的計算量，并且PSO優(yōu)化打破了比例分配重采樣的全局并行實現(xiàn)，而是用局部重采樣代替，這樣將會帶來估計誤差，在實際中需要權(quán)衡PSO優(yōu)化的次數(shù)。

5 完整的分布式粒子濾波算法

用于粒子濾波的總粒子數(shù)目為N，用于計算的子節(jié)點數(shù)目為S，下標(biāo)k表示時間點，i表示第i個粒子，j表示第j個子節(jié)點。則分布式粒子濾波算法描述如下：

1）子節(jié)點采樣。

2）子節(jié)點權(quán)值計算。

聚合數(shù)據(jù)如下：

4）中心節(jié)點進行狀態(tài)估計。

5）中心節(jié)點進行重采樣。判斷N＿resample，若N＿resample（j）／N＜α，則置位PSO優(yōu)化標(biāo)志，將數(shù)據(jù)｛u′，Ck，F(xiàn)｝發(fā)送到子節(jié)點，其中α為判斷負載是否失衡的閾值，F(xiàn)為PSO優(yōu)化的標(biāo)志位。

6）子節(jié)點進行重采樣。若F相應(yīng)位被置位，那么相應(yīng)子節(jié)點先調(diào)用PSO算法，然后采用局部DPF中的局部重采樣；若F沒有被置位，則子節(jié)點進行比例分配重采樣。然后回到第1）步。

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

本文在Matlab平臺下進行仿真，假設(shè)目標(biāo)在二維平面內(nèi)做勻速直線運動。

系統(tǒng)狀態(tài)轉(zhuǎn)移方程為

觀測方程為

其中，F(xiàn)＝［1T00 0100 001T 0001］；G＝［T2／20 T0 0T2／2 0T］；Wk為過程噪聲，其協(xié)方差矩陣Q＝diag（0.12，0.12）；Vk為觀測噪聲，這里采用閃爍噪聲，其密度函數(shù)可表示為

協(xié)方差P＝（1－eta）＊P1＋eta＊P2，其中

取σr1＝σr2＝50；σb1＝1°，σb2＝5°；eta＝0.1

初始狀態(tài)選為x0＝［50000 300 50000 －100］T，仿真點為50，采用不同的粒子數(shù)和用于計算的子節(jié)點數(shù)，分別進行100次蒙特卡羅實驗，本文的算法和其它兩種分布式粒子濾波算法的誤差比較如表1所示。

表1 不同粒子數(shù)、子節(jié)點數(shù)時的誤差

表1中的結(jié)果可以表示成圖4～圖6的形式。分別對應(yīng)節(jié)點數(shù)為1、5和10時三種算法的誤差曲線。

圖4 1個節(jié)點的誤差曲線

圖5 5個節(jié)點的誤差曲線

圖6 10個節(jié)點的誤差曲線

對于GDPF誤差在粒子數(shù)增加到2000時可以降到很小，但是該算法需要傳輸大量的粒子，增加通信損耗；對于LDPF算法，由于局部重采樣導(dǎo)致誤差比同粒子數(shù)時GDPF的要大，但是不需要傳送任何粒子信息，實現(xiàn)方便；本文所采用的算法也不需要傳送任何粒子信息，并且由于引入了PSO算法，優(yōu)化了用于估計的粒子，因此誤差降低，但是PSO運行的次數(shù)越多會帶來額外的計算開銷。可以發(fā)現(xiàn)在粒子數(shù)為100時，容易出現(xiàn)負載失衡，也就是調(diào)用PSO算法的次數(shù)會很多，誤差減小的同時帶來很大的PSO計算開銷，當(dāng)粒子數(shù)上升為2000時，調(diào)用PSO算法的次數(shù)為0，也就是沒有出現(xiàn)負載失衡的情況，此時誤差值和節(jié)點數(shù)為1時相近，即比例分配重采樣的并行計算可以達到與集中式計算同樣的精度，證明了并行處理過程是可行的。

［1］Gordon N，Salmond D.Novel appr oach to non－linear and non－Gaussian Bayesian state estimation［J］.Procof I nstitute E lectr ic E ngineering，1993，140（2）：107－113.

［2］Coates M.Distributed particle fiflters for sensor networks.PIPSN，2004：99－107.

［3］Bashi A S，Jilkov V P，Li X R，et al.Distributed implementations of particlefilters［C］／／Proceeding of the Sixth International Conference of Information Fusion，2003，2：1164－1171.

［4］Sheng X H，Hu Y H.Parameswaran Ramanathan.Distributed particle filter with GMM approximation for multiple targets［C］／／4th International Symposium on.Information Pocessing in Sensor Networks，2005：181－188.

［5］Bolic M，Djuric P M，Hong S J.Resamping alogorithms and architectures for distributed partical filters［J］.IEEE Trans.on Signal Processin，2005，53（7）：2442－2450.

［6］R V Merwe，A Doucet，Nando De Freitas，et al.The Unsented Particle Filter［R］／／Technical Report，CUED／F－INPENG／TR 380.UK：Engineering Department，Cambridge University，2000.

［7］J D Hol.Resampling in Particle Filters［R］／／Intership report，LiTH－ISY－EX－ET－0283－2004.Sweden：Link ping University，2004.

［8］Miodrag Bolic'，Peter M Djuric'，Sangjin Hong.New Resampling Algorithms for Particle Filters［J］.ICASSP2003（S1845－5921），2003，2：589－592.

［9］KROHLING R A.Gaussian swarm：a novel particle swarm optimization algorithm［C］／／Proceedings of the IEEE Conference on Cybernetics and Intelligent Systems（CIS）.Singapore：IEEE Press，2004：372－376.

［10］Peng zhang，Ming Li，Yan Wu.An improved particle filter algorithm based on Markov Random Field modeling in stationary wavelet domain for SAR image despeckling［J］.Pattern Recognition Letters，2012，33（10）：1316－1328.