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

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

0 引 言

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

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

1 主要算法基本原理

1.1 布谷鳥優(yōu)化算法

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

(1)布谷鳥一次產(chǎn)一個(gè)卵，并隨機(jī)選擇鳥窩位置孵化。

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

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

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

(1)

1.2 粒子濾波基本原理

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

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

(2)

(3)

(4)

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

(5)

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

2.1 改進(jìn)粒子的目標(biāo)函數(shù)

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

(6)

2.2 粒子權(quán)值更新

將布谷鳥算法融入PF的核心思想是，將粒子濾波中的每個(gè)粒子都進(jìn)行迭代尋優(yōu)，從而，粒子向后驗(yàn)概率密度值高的區(qū)域移動(dòng)，進(jìn)而提高估計(jì)結(jié)果的準(zhǔn)確性。但是，各粒子在狀態(tài)空間的位置也發(fā)生了改變，各粒子所表示的分布密度函數(shù)p(xk|y1:k-1)被改變，那么粒子濾波基于貝葉斯的理論基礎(chǔ)就丟失。所以，在粒子更新時(shí)，對(duì)粒子權(quán)值進(jìn)行更新。參照文獻(xiàn)[3]，利用式(6)進(jìn)行粒子位置更新時(shí)更新粒子權(quán)值。公式如下

(7)

2.3 CS-PF算法和步驟

基于2.1和2.2節(jié)的改進(jìn)，提出的CS-PF算法流程圖，如圖1所示。

圖1 CS-PF算法流程

CS-PF算法的詳細(xì)過程如下。

步驟1 初始化N個(gè)粒子

步驟2 為了讓粒子盡量向真實(shí)區(qū)域移動(dòng)，又要避免最終收斂，這里利用小規(guī)模多種群思想，將k時(shí)刻N(yùn)個(gè)粒子平均劃分為m個(gè)子集，每個(gè)子集粒子為

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

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

步驟3 模擬布谷鳥鳥窩尋優(yōu)

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

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

(8)

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

(9)

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

(10)

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

FOR-2end

FOR-1end

步驟4 歸一化粒子權(quán)值

(11)

步驟5 輸出狀態(tài)

(12)

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

3 實(shí)驗(yàn)仿真

通過實(shí)驗(yàn)仿真的方法，確定重要參數(shù)pa的取值范圍，并驗(yàn)證算法的有效性。

3.1 實(shí)驗(yàn)內(nèi)容

利用以上實(shí)驗(yàn)環(huán)境和第2節(jié)提出的CS-PF進(jìn)行以下3項(xiàng)實(shí)驗(yàn)，以充分說明算法的有效性。

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

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

(13)

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

以上3項(xiàng)仿真實(shí)驗(yàn)基于以下過程狀態(tài)方程和觀測(cè)方程

(14)

(15)

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

3.2 實(shí)驗(yàn)結(jié)論

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

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

圖2 pα=0.05濾波效果

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

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

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

圖4 pα= 0.10濾波效果

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

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

圖6 pα=0.15濾波效果

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

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

為了進(jìn)一步充分說明pα對(duì)算法性能和粒子多樣的影響，以及選取pα合適的取值范圍，采用同樣的方法重復(fù)以上實(shí)驗(yàn)直到pα= 0.60，經(jīng)統(tǒng)計(jì)得出以下結(jié)論，見表1。

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

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

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

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

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

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

程序運(yùn)行環(huán)境：處理器，Intel(R) Core(TM) i5-4200U CPU @ 1.60 GHz 2.30 GHz；內(nèi)存，4.00 G；操作系統(tǒng)，64位Windows 7。實(shí)驗(yàn)結(jié)果見表3。

表3 不同算法RMSE與運(yùn)算量對(duì)比

在表3中，本文算法CS-PF的RMSE=2.1735，濾波精度低于文獻(xiàn)[15]的imp-WOPF算法，卻優(yōu)于標(biāo)準(zhǔn)PF，以及FA-PF和ADE-PF，但運(yùn)行時(shí)間高于標(biāo)準(zhǔn)PF，而少于其它幾種濾波算法。主要是因?yàn)镃S-PF中融入了基于布谷鳥算法的粒子尋優(yōu)過程，但是，該過程的運(yùn)算量卻低于螢火蟲算發(fā)、自適應(yīng)差分進(jìn)化的粒子濾波和權(quán)值改進(jìn)算法。

綜合上述實(shí)驗(yàn)結(jié)果，得出以下結(jié)論。

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

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

3.3 實(shí)際應(yīng)用

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

4 結(jié)束語

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