基于布谷鳥算法優化的粒子濾波

白曉波，邵景峰，和 征，田建剛

(1.西安工程大學 管理學院，陜西 西安 710048；2.陸軍邊海防學院 工程基礎系，陜西 西安 710108)

0 引 言

很多學者對粒子濾波(particle filter，PF)[1,2]進行了深入研究。文獻[3]利用螢火蟲算法迭代尋優的思想優化粒子分布，但螢火蟲算法本身存在著對于初始解分布的依賴性；文獻[4]將粒子濾波算法用于嵌入式系統進行目標跟蹤；K.Madhan Kumar等[5]提出混合多核偏最小二乘(PLS)粒子濾波方法，基于多核偏最小二乘法，結合粒子濾波算法提高遙感圖像中道路提取的準確度。其它如文獻[6-8]提出了相關的改進算法，都集中于粒子采樣和粒子權值計算的改進；文獻[9]主要解決粒子數的自適應選擇問題，在保證適當粒子數的情況下，保證濾波精度，提高粒子濾波的實時性；陳世明等[10]提出一種基于引力場的粒子濾波算法以避免傳統粒子濾波算法中粒子貧化與退化現象。

對粒子濾波的研究在一定程度上解決了粒子權值退化和粒子貧化問題，但增加了算法運算量，影響算法實時性。融入新的智能優化算法，在解決粒子權值退化和粒子貧化時，有效降低運算量有著重要的理論意義和應用價值。因此，提出了基于布谷鳥算法改進的粒子濾波(CS-PF)，并實驗仿真得出重要參數的取值范圍。

1 主要算法基本原理

1.1 布谷鳥優化算法

劍橋大學學者YANG Xinshe和DEB Suash通過對布谷鳥尋窩產卵的行為的模擬，提出了一種全局搜索的Cuckoo Search(CS)算法[11]。其基本思想基于3個理想假設。

(1)布谷鳥一次產一個卵，并隨機選擇鳥窩位置孵化。

(2)在隨機選擇的一組鳥窩中，最好的鳥窩會被保留到下一代。

(3)可利用的鳥窩數量n為常量，宿主發現布谷鳥蛋的概率為pα∈[0,1]。

基于以上3個理想條件，布谷鳥尋找孵化鳥窩的路徑和位置更新公式如下

(1)

1.2 粒子濾波基本原理

粒子濾波基本原理是利用貝葉斯濾波理論解決狀態估計問題，計算粒子權重，使用蒙特卡羅序列[12]方法確定狀態的后驗概率，其過程如下。

已知狀態的初始概率密度函數為p(x0|z0)=p(x0)，狀態更新方程為

(2)

(3)

(4)

然而，在實際應用時存在抽取有效樣本困難的問題，為此，引入了重要性采樣方法[12]，主要計算公式如下

(5)

2 布谷鳥算法優化的粒子濾波(CS-PF)

2.1 改進粒子的目標函數

通過增加擾動項來進一步保證粒子的多樣性。從式(1)分析，第i個鳥窩在第t+1代的鳥窩位置與第t代位置加上Levy隨機搜索路徑，由于Levy飛行是一個非高斯的隨機過程[13]。理論上，如果將其應用于粒子進化過程，短步長搜索有利于粒子向高似然區移動，偶爾大步長搜索，避免粒子陷入局部最優，且能保證粒子向下迭代時的多樣性。但是，仍然存在局部最優的可能性。因此，增加擾動項對式(1)進行優化，如下所示

(6)

2.2 粒子權值更新

將布谷鳥算法融入PF的核心思想是，將粒子濾波中的每個粒子都進行迭代尋優，從而，粒子向后驗概率密度值高的區域移動，進而提高估計結果的準確性。但是，各粒子在狀態空間的位置也發生了改變，各粒子所表示的分布密度函數p(xk|y1:k-1)被改變，那么粒子濾波基于貝葉斯的理論基礎就丟失。所以，在粒子更新時，對粒子權值進行更新。參照文獻[3]，利用式(6)進行粒子位置更新時更新粒子權值。公式如下

(7)

2.3 CS-PF算法和步驟

基于2.1和2.2節的改進，提出的CS-PF算法流程圖，如圖1所示。

圖1 CS-PF算法流程

CS-PF算法的詳細過程如下。

步驟1 初始化N個粒子

步驟2 為了讓粒子盡量向真實區域移動，又要避免最終收斂，這里利用小規模多種群思想，將k時刻N個粒子平均劃分為m個子集，每個子集粒子為

FOR-1:m=1,2,…,log2N

FOR-2:j=1,2,…,N/log2N

步驟3 模擬布谷鳥鳥窩尋優

(1)生成隨機數r∈[0,1]，若r>pα，執行步驟(2)，否則，執行步驟(3)。

(2)根據式(6)，對第m子集第k時刻第j個粒子的位置更新公式重新表示，如下所示

(8)

參數?(通常?=1)和運算符⊕與式(1)相同，L(λ)為Levy隨機搜索路徑，且L～u=s-λ(1<λ≤3)，s為Levy飛行隨機步長，為了便于計算，通常使用下列公式計算Levy隨機數

(9)

u,v為服從正態分布的隨機數，u～N(0,σ2u),v～N(0,σ2v)， 0<β<2，通常β=1.5，σu、σv的取值如下

(10)

Γ為伽瑪函數。返回步驟(1)。

FOR-2end

FOR-1end

步驟4 歸一化粒子權值

(11)

步驟5 輸出狀態

(12)

由于布谷鳥算法具有全局最優特點，為了讓粒子群整體向真實值附近移動，且避免最終收斂，我們通過設置隨機數r是否小于閾值pα作為迭代的終止條件。因為通過粒子位置迭代更新的目的，并不是讓所有粒子都集中于高似然區，只是為了讓粒子集在高似然區合理分布，否則就會降低粒子的多樣性。但是，pα的取值必須合適，如果迭代次數過多，增加算法時間復雜度，CS-PF的實時性能會受到很大影響；如果迭代次數過少，粒子集發散，不能合理集中于高似然區，就不能提高濾波精度。因此，在下一節的主要內容之一就是通過實驗仿真的方法確定閾值pα的取值范圍。

3 實驗仿真

通過實驗仿真的方法，確定重要參數pa的取值范圍，并驗證算法的有效性。

3.1 實驗內容

利用以上實驗環境和第2節提出的CS-PF進行以下3項實驗，以充分說明算法的有效性。

(1)確定參數pα。在CS-PF算法中，是否更新粒子位置和算法運算量，主要取決于是否r>pα，r為隨機數，r∈[0,1]，pα∈[0,1]，從理論上分析pα越小，通過式(8)粒子位置更新算法更新的粒子越多，較多粒子分布于高似然區，從而濾波精度較高，但也會較大的增加算法運算量；若pα越大，被更新位子的粒子越少，也就只有較小的粒子分布于高似然區，粒子保持了多樣性，但濾波精度降低，算法運算量也小，實時性較高。因此，需通過實驗仿真的方法，在濾波精度和算法實時性之間找到pα合適的取值范圍。

(2)檢驗在高斯，非高斯噪聲下CS-PF的性能。分別在量測噪聲符合高斯分布、均勻分布、伽馬分布、瑞利分布的噪聲序列對CS-PF重復25次實驗，對本文算法運算性能進行分析，并計算結果的均方根誤差均值，計算公式如下

(13)

(3)CS-PF與其它算法對比分析。本文算法與解決非線性問題的經典算法，如文獻[3]的基于螢火蟲算法智能優化的粒子濾波(FA-PF)、文獻[14]的基于自適應差分進化的粒子濾波(ADE-PF)、文獻[15]的imp-WOPF和標準PF，在相同的過程方差和量測方差下進行濾波效果對比分析。

以上3項仿真實驗基于以下過程狀態方程和觀測方程

(14)

(15)

式(14)為狀態方程，式(15)為觀測方程，式中：Q為制造過程噪聲方差，R為量測噪聲方差，rand為高斯分布的隨機數。由于此系統為強非線性系統，且似然函數為雙峰狀[3]，一般濾波方法很難對其進行有效處理。參考現有文獻，式(8)中，參數?=1，步長因子α=0.1。式(9)中β=1.5。

3.2 實驗結論

(1)參數pα取值范圍。初始狀態x=1，過程噪聲Q和量測噪聲R的方差均為5，模擬時序k=1,2,…,600，粒子數N=100。為了取得合適的pα，取pα=0.05+0.05i,(i=0,1,2,…,11)，并記錄下pα等于不同值時，被更新的粒子平均數和粒子的運動趨勢，結果如下。

(a)pα=0.05。即粒子在向下迭代時，理論上，r為均勻分布的隨機數，pα越小，被更新的粒子越多，算法運算量越大，通過多次迭代，多數粒子將會向高似然區域移動，故不能保證粒子的多樣性。算法濾波結果如圖2所示，k=232時的粒子分布情況如圖3所示。

圖2 pα=0.05濾波效果

圖3 pα= 0.05,k=232粒子分布情況

在圖2中，平均每個粒子位置更新19.02次。RMSE=2.12。圖3中，濾波進入后續階段，PF算法的粒子大多數集中于少數狀態值，從而失去了粒子多樣性，不利于濾波的狀態估計。而CS-PF在高似然區和低似然區都有一定數量的粒子分布，可以在一定程度上保證粒子的多樣性。但是，在當前狀態下，算法運算量較高。

(b)pα=0.10。平均每個粒子位置更新9.23次，RMSE=2.29，如圖4所示，在k=232時，粒子分布情況，如圖5所示。

圖4 pα= 0.10濾波效果

圖5 pα=0.10,k=232粒子分布情況

(c)pα=0.15。平均每個粒子位置更新6.04次，RMSE=1.67，如圖6所示。在k=232時，粒子分布情況，如圖7所示。

圖6 pα=0.15濾波效果

圖7 pα=0.15,k=232粒子分布情況

從以上3步實驗可以初步進行以下推測。在pα逐漸增大時，CS-PF濾波精度略有變化，RMSE逐漸增大，但變化較小，而分布于高似然區的粒子數在逐漸減少，低似然區的粒子數在增加，從實驗環節說明pα的增大有利于促進粒子的多樣性，而其減小對算法濾波精度的降低有限，卻較大增加了算法運算量。因此，參數pα的取值是影響算法精度和算法運算量的重要因素，同時，也極大地影響到粒子的多樣性。

為了進一步充分說明pα對算法性能和粒子多樣的影響，以及選取pα合適的取值范圍，采用同樣的方法重復以上實驗直到pα= 0.60，經統計得出以下結論，見表1。

合適的pα∈[0.20,0.35]。從以上實驗得出以下結論，算法運算量、濾波精度都和pα有一定的相關性。若pα<0.2，平均21.59個粒子在高似然區分布，平均每個粒子更新達20.66次，但同時也喪失了粒子多樣性，且算法運算量大增，若粒子數較大，就極大影響了算法實時性；反之，若pα>0.5，絕大多數粒子位置都沒有更新，平均更新次數在0.65左右，也就是每個粒子更新不到一次，分布于高似然區的粒子數就相對少很多，不能有效提高算法濾波精度，但算法運算量較低。因此，結合表1中高似然區平均粒子數的均值為10，綜合考慮，既要保證足夠多的粒子分布于高似然區，以保證濾波精度，又能以較低的運算量保證算法實時性，所以最為合適的pα∈[0.20,0.35]之間。

表1 pα取值、粒子位置更新情況與RMSE

(2)利用MatLab R2012a的gamrnd函數、rand函數、raylrnd函數和randn函數，生成伽馬分布、均勻分布、瑞利分布和高斯分布的量測噪聲序列。取pα=0.25，設初始值為x=1，噪聲序列服從均值為0，過程噪聲方差Q=5，量測噪聲方差R=5，粒子數為N=100，濾波時序k=600，分析噪聲在高斯分布、均勻分布、伽馬分布、瑞利分布下，對算法的性能影響。針對不同噪聲分布分別做25次獨立實驗，計算不同噪聲分布下的RMSE均值。結果見表2。

從表2可以得出以下結論，過程噪聲與量測噪聲不變，高斯分布、均勻分布、伽馬分布和瑞力分布的均方根誤差分別為2.1722、2.1963、2.2491和2.2687，說明CS-PF針對強非線性系統，在固定的過程方差和量測方差、高斯分布和非高斯分布的量測噪聲下濾波效果幾乎一致。

表2 CS-PF在不同分布下的RMSE均值

(3)取pα=0.25，設初始值為x=1，噪聲序列服從均值為0，過程噪聲方差Q=5，量測噪聲方差R=5，粒子數為N=100，濾波時序k=600，噪聲分布為高斯分布。對文獻[3]的基于螢火蟲算法智能優化的粒子濾波(PA-PF)，文獻[14]的ADE-PF和文獻[15]的imp-WOPF對比研究。針對每種算法進行了25次濾波，并計算每種算法濾波的均方根誤差均值，計算方法和式(13)相同，并統計分析每種濾波各自的運行時間均值。

程序運行環境：處理器，Intel(R) Core(TM) i5-4200U CPU @ 1.60 GHz 2.30 GHz；內存，4.00 G；操作系統，64位Windows 7。實驗結果見表3。

表3 不同算法RMSE與運算量對比

在表3中，本文算法CS-PF的RMSE=2.1735，濾波精度低于文獻[15]的imp-WOPF算法，卻優于標準PF，以及FA-PF和ADE-PF，但運行時間高于標準PF，而少于其它幾種濾波算法。主要是因為CS-PF中融入了基于布谷鳥算法的粒子尋優過程，但是，該過程的運算量卻低于螢火蟲算發、自適應差分進化的粒子濾波和權值改進算法。

綜合上述實驗結果，得出以下結論。

(1)將布谷鳥算法的路徑尋優思想融入粒子濾波代替粒子重采樣，能夠解決粒子權值退化問題，并保持粒子多樣性。

(2)相比標準粒子濾波，在運算效率降低12%的情況下，通過合理的粒子更新閾值pα=0.25，能夠提高濾波精度18%。

3.3 實際應用

為了進一步驗證算法的可靠性，在紡織實際生產環境中對CS-PF算法進行驗證。計算機軟硬件環境為一臺Dell服務器(處理器為2個Intel(R) Xeon(R) CPU E5-2407 2.40 GHz，Windows Server 2008(64位)操作系統，32 GB內存，8 TB硬盤容量，通信帶寬：1 Gb/s)，對織布車間300臺織機的原始數據進行實時在線監測。通過實際應用表明：在針對非平穩、非線性的織造過程數據時，具有很好的適用性，結果誤差小于3%，誤差范圍小于生產環境5%的要求。結果表明：本文算法能夠有效地處理織機聲非線性織造過程的非平穩數據，有效地保證織造過程中采集坯布的質量和產量數據的準確性，從而，很好地滿足生產管理的需要，以及生產管理、統計分析數據的準確性也大大提高。

4 結束語

基于布谷鳥路徑尋優思想融入粒子濾波，以解決粒子濾波中的權值退化和粒子貧化問題。實驗結果表明：CS-PF在利用布谷鳥尋優思想代替粒子重采樣后，能夠保持粒子多樣性和濾波精度。與其它改進的粒子濾波算法相比，CS-PF具有以下有點：①粒子位置的迭代更新采用Levy飛行模式，通過閾值控制迭代次數，能夠保持粒子多樣性，克服陷入局部最優。使得部分粒子分布于高似然區，部分分布于低似然區。②CS-PF調節參數，易于編程實現，容易應用于實際工程問題。③在運算效率降低較少的情況下，運算精度有了較大提高。未來主要的研究應著眼于，結合現有研究基礎，在實際工程問題中自適應的控制粒子位置更新閾值，以及粒子數的自適應選擇上，以進一步降低運行時間，提高算法的實時性。