帶負載失衡調節的分布式粒子濾波算法研究＊

查巍巍 張 勇 應曉慧 周 源 李知杰

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

1 引言

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

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

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

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

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

2 粒子濾波算法

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

遞推過程可分為預測和更新兩步

預測：

更新：

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

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

一般選擇狀態轉移概率密度函數p（xk｜xk－1）作為重要密度函數。

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

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

2）權值計算。

并且歸一化

則可得k時刻的狀態估計為

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

4）時刻k＝k＋1，轉到第1）步。

3 分布式處理

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

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

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

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

Step1：中心節點根據子節點的權重和進行第一次一次比例分配重采樣。

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

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

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

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

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

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

根據Step1得

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

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

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

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

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

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

一般狀態空間模型分為狀態轉移模型和量測模型：

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

PSO更新過程為

經過M次迭代后，粒子的最優值符合某閾值ε時，說明粒子群已經分布在真實狀態附近，則停止優化。

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

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

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

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

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

1）子節點采樣。

2）子節點權值計算。

聚合數據如下：

4）中心節點進行狀態估計。

5）中心節點進行重采樣。判斷N＿resample，若N＿resample（j）／N＜α，則置位PSO優化標志，將數據｛u′，Ck，F｝發送到子節點，其中α為判斷負載是否失衡的閾值，F為PSO優化的標志位。

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

6 仿真結果與分析

本文在Matlab平臺下進行仿真，假設目標在二維平面內做勻速直線運動。

系統狀態轉移方程為

觀測方程為

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

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

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

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

表1 不同粒子數、子節點數時的誤差

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

圖4 1個節點的誤差曲線

圖5 5個節點的誤差曲線

圖6 10個節點的誤差曲線

對于GDPF誤差在粒子數增加到2000時可以降到很小，但是該算法需要傳輸大量的粒子，增加通信損耗；對于LDPF算法，由于局部重采樣導致誤差比同粒子數時GDPF的要大，但是不需要傳送任何粒子信息，實現方便；本文所采用的算法也不需要傳送任何粒子信息，并且由于引入了PSO算法，優化了用于估計的粒子，因此誤差降低，但是PSO運行的次數越多會帶來額外的計算開銷。可以發現在粒子數為100時，容易出現負載失衡，也就是調用PSO算法的次數會很多，誤差減小的同時帶來很大的PSO計算開銷，當粒子數上升為2000時，調用PSO算法的次數為0，也就是沒有出現負載失衡的情況，此時誤差值和節點數為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.