王龍生， 顧 浩， 余云智

(江蘇自動化研究所，江蘇 連云港 222006)

0 引言

近年來，粒子濾波(Particle Filter，PF)［1－2］因適用于任何能用狀態(tài)空間模型表示的非高斯背景的非線性隨機(jī)系統(tǒng)，精度可以逼近最優(yōu)估計(jì)，是一種很有效的非線性濾波技術(shù)而備受關(guān)注，被廣泛應(yīng)用于計(jì)算機(jī)視覺、信號處理、定位導(dǎo)航及目標(biāo)跟蹤等諸多工業(yè)領(lǐng)域［3］。

粒子濾波的核心思想是通過從后驗(yàn)概率中抽取的隨機(jī)狀態(tài)粒子來表達(dá)其分布，是一種順序重要性采樣方法。該算法存在的最主要的問題是需要用大量的樣本數(shù)量才能很好地近似系統(tǒng)的后驗(yàn)概率密度。另外，重采樣階段會造成樣本有效性和多樣性的損失，導(dǎo)致出現(xiàn)樣本貧化現(xiàn)象［4－5］。如何保持粒子的有效性和多樣性，克服樣本貧化，一直是研究的重點(diǎn)。

為了解決多樣性損失問題，文獻(xiàn)［6］針對重采樣算法的改進(jìn)提出了遺傳重采樣粒子濾波算法，通過引入遺傳交叉變異操作，在一定程度上解決了樣本貧化問題，但是遺傳算法實(shí)時性不好，在工程應(yīng)用上還有待提高。微分進(jìn)化(Differential Evolution，DE) 算法［7］，是一種采用浮點(diǎn)矢量編碼的連續(xù)空間中進(jìn)行隨機(jī)搜索的優(yōu)化算法，已被證明在許多優(yōu)化問題中都表現(xiàn)出優(yōu)于自適應(yīng)模擬退火算法、粒子群算法、遺傳算法的性能［8］。本文通過引入微分進(jìn)化算法進(jìn)行重采樣，利用其交叉變異操作創(chuàng)造更多的粒子，增加樣本多樣性，抑制退化現(xiàn)象，并將微分進(jìn)化重采樣和殘差重采樣進(jìn)化組合優(yōu)化，解決算法實(shí)時性的問題，蒙特卡羅仿真表明此方法的有效可行性。

1 粒子濾波

粒子濾波是一種基于蒙特卡羅方法和貝葉斯估計(jì)的統(tǒng)計(jì)濾波方法，它首先依據(jù)系統(tǒng)狀態(tài)向量的經(jīng)驗(yàn)條件分布在狀態(tài)空間產(chǎn)生一組隨機(jī)樣本(粒子)，然后根據(jù)量測不斷調(diào)整粒子的權(quán)重和位置，通過調(diào)整后的粒子信息修正最初的經(jīng)驗(yàn)條件分布。

為了得到正確的狀態(tài)估計(jì)，通常希望粒子權(quán)值的方差盡可能趨近于零，然而，序貫蒙特卡羅模擬方法一般都存在權(quán)值退化問題。在實(shí)際計(jì)算中，隨著濾波迭代次數(shù)的增加，粒子權(quán)值的方差隨著時間增大，狀態(tài)空間中的有效粒子數(shù)較少，即當(dāng)濾波一段時間后，粒子的權(quán)值會出現(xiàn)兩極分化，少數(shù)粒子具有很大的權(quán)值，對結(jié)果起主導(dǎo)作用，而其余大多數(shù)的粒子權(quán)值很小，對結(jié)果幾乎沒有影響，從而使得大量計(jì)算浪費(fèi)在求解幾乎不起任何作用的粒子的更新上，結(jié)果粒子集無法表達(dá)實(shí)際的后驗(yàn)概率分布，這就是粒子濾波算法的退化問題。一般通過引入重采樣算法基本上消除了粒子濾波算法的退化現(xiàn)象，但是重采樣算法帶來的另一個問題是粒子樣本枯竭。以上諸多因素嚴(yán)重貧化了樣本支撐集，極大地影響了估計(jì)結(jié)果，不利于狀態(tài)估計(jì)與跟蹤，當(dāng)估計(jì)那些在很長時間維持不變的量時，退化現(xiàn)象更為嚴(yán)重。因此，一個好的重采樣算法應(yīng)該在減少權(quán)值較小的粒子數(shù)目和增加粒子樣本多樣性之間進(jìn)行有效折衷。

2 組合進(jìn)化重采樣算法

2.1 微分進(jìn)化重采樣

使用一般重采樣算法，粒子多樣性損失是不可避免的，本文通過引入微分進(jìn)化算法進(jìn)行重采樣，經(jīng)過交叉變異操作在一定程度上可以克服粒子濾波的樣本貧化問題。微分進(jìn)化算法具有高效性、收斂性、魯棒性等優(yōu)點(diǎn)［9］，相比于遺傳算法，它保留了基于種群的全局搜索策略，采用實(shí)數(shù)編碼、基于差分的簡單變異操作和一對一的競爭生存策略，降低了遺傳操作的復(fù)雜性。微分進(jìn)化重采樣算法中交叉、變異操作定義如下［10］。

1)交叉算子。

在DE中，執(zhí)行交叉運(yùn)算可以增加種群多樣性。常見的有二項(xiàng)交叉和冪交叉。由下式計(jì)算交叉向量

式中:CR∈［0，1］;jrand∈{1，2，…，N}，j=1，2，…，N。CR為交叉概率，CR越大，DE收斂越快;CR值越小，DE魯棒性越好，但增加了執(zhí)行時間。

2)變異算子。

目標(biāo)向量表示為 Xi=(x0，i，…，xN－1，i)T;變異向量表示為 Vi=(v0，i，…，vN－1，i)T，i=0，1，…，NP－ 1;N 為目標(biāo)向量的維數(shù);NP為種群數(shù)。常用的變異操作為(下標(biāo)G表示第G代)

式中:r1，r2，r3∈{1，…，NP}，為互不相等的整數(shù)，且它們與i不相等;F∈(0，1)，用于控制變異比例。

微分進(jìn)化粒子濾波流程如圖1所示。

圖1 微分進(jìn)化粒子濾波流程圖Fig.1 Flow chart of the different evolution particle filter

微分進(jìn)化重采樣算法步驟如下。

①初始化采樣。

從初始分布p(x0)中采樣N個粒子樣本集xi0，i=1，2，…，N。

②微分進(jìn)化操作。

將N個粒子作為粒子初始種群，利用式(1)和式(2)進(jìn)行交叉變異處理，得到N個新樣本，將原始樣本集與新樣本合并為大小為2N的候選集。

③狀態(tài)估計(jì)。

將權(quán)值大于w0的樣本組成一個集合，取其均值作為預(yù)測值。w0為常數(shù)。

④粒子重采樣。

計(jì)算候選集中各樣本的權(quán)值wi，并降序排序，選擇前一半作為下一代種群。

⑤粒子更新。

進(jìn)入下一時刻，粒子進(jìn)行系統(tǒng)狀態(tài)轉(zhuǎn)移，然后轉(zhuǎn)至步驟②。

2.2 殘差重采樣

針對標(biāo)準(zhǔn)粒子濾波算法存在的權(quán)值退化現(xiàn)象，文獻(xiàn)［11］提出了自舉粒子濾波算法，該算法在遞推過程中引入了重采樣的思想以克服退化問題。重采樣的基本思想是:通過對用粒子和相應(yīng)權(quán)值表示的后驗(yàn)概率密度函數(shù)重采樣產(chǎn)生新的粒子集，即在滿足條件下，將粒子集合更新為。在常用的重采樣算法中，殘差重采樣法具有效率高、易于實(shí)現(xiàn)的特點(diǎn)。設(shè)其中為取整操作。殘差重采樣采用新的權(quán)值選擇余下的個粒子。殘差重采樣的主要過程有:

3)對于每個μi，尋找歸一化權(quán)值累計(jì)量大于或等于 μi的最小標(biāo)號 m，即。當(dāng) μi落在區(qū)間［λm－1，λm］時，被復(fù)制一次。這樣，每個粒子經(jīng)重采樣后的個數(shù)為步驟3)中被選擇的若干粒子數(shù)目與Ni之和。

2.3 組合算法

樣本貧化是粒子濾波算法的一個主要問題，微分進(jìn)化算法重采樣由于引入了微分進(jìn)化算法的交叉變異操作，獲得了更多的可選擇粒子，使得粒子的多樣性得到保障，從而抑制了粒子退化現(xiàn)象。但同時，由于交叉變異操作帶來了一定的實(shí)時性問題，因此加大了濾波時間。常用的重采樣算法具有簡潔、實(shí)時性高的特點(diǎn)，其中以殘差重采樣算法效率最高，本文通過將微分進(jìn)化重采樣算法和殘差重采樣算法進(jìn)行組合重采樣，從而既保持了粒子的多樣性又提高了算法實(shí)時性。

對于樣本貧化現(xiàn)象，引入有效粒子采樣尺度

Neff變小，就意味著粒子存在退化現(xiàn)象。如果Neff=1，即一個粒子的權(quán)值為1，其他粒子的權(quán)值為0，說明粒子退化嚴(yán)重;如果Neff=1/N，則表示每個粒子的權(quán)值為1/N，粒子未退化。因此，為了克服樣本貧化現(xiàn)象，當(dāng)Neff低于某一門限值時需要進(jìn)行多樣性恢復(fù)。

組合算法流程如圖2所示。

圖2 組合算法流程圖Fig.2 Flow chart of the combined algorithm

組合算法步驟如下。

1)初始化采樣。

從初始分布p(x0)中采樣N個粒子樣本集xi0，i=1，2，…，N。

2)重要性權(quán)值計(jì)算。

在k時刻，更新粒子權(quán)值

并且歸一化

3)有效粒子采樣尺度計(jì)算。

4)重采樣算法選擇。

設(shè)Neff門限值為 N0，當(dāng) Neff低于 N0時，轉(zhuǎn)至步驟5);否則轉(zhuǎn)至步驟6)。

5)微分進(jìn)化重采樣。

使用交叉變異操作，增加粒子樣本數(shù)，進(jìn)行狀態(tài)估計(jì)和粒子更新，轉(zhuǎn)至步驟8)。

6)殘差重采樣。

使用殘差重采樣算法，進(jìn)行粒子重采樣，并行重新分配粒子權(quán)值，所有粒子權(quán)值為1/N。

7)狀態(tài)估計(jì)。

8)粒子更新。

進(jìn)入下一時刻，粒子進(jìn)行系統(tǒng)狀態(tài)轉(zhuǎn)移，然后轉(zhuǎn)至步驟2)。

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

本文選用粒子濾波器性能評測時常用的一個基準(zhǔn)(benchmark)［7］進(jìn)行仿真，即單變量的非靜態(tài)經(jīng)濟(jì)增長模型，它是計(jì)量經(jīng)濟(jì)學(xué)中常見的問題。該系統(tǒng)具有高度非線性，似然函數(shù)呈雙峰，使得用傳統(tǒng)方法來解決非常困難。

該問題的過程模型和量測模型為

式中，w(t)和v(t)為零均值高斯噪聲，方差分別為10與1，初始狀態(tài)取x(0)=10。

粒子權(quán)值計(jì)算式為

式中:y為觀測值;y'為觀測的估計(jì)。

對標(biāo)準(zhǔn)的基于殘差重采樣的粒子濾波算法和微分進(jìn)化粒子濾波算法以及本文的組合算法進(jìn)行實(shí)驗(yàn)，取粒子數(shù)目為1000個，采樣時間序列為100，交叉概率和變異概率分別為0.95和0.85，種群數(shù) NP為5。程序運(yùn)行于Windows XP系統(tǒng)，電腦配置為酷睿雙核Q8400處理器，2 G內(nèi)存，運(yùn)行環(huán)境為Matlab 2008。

表1為500次蒙特卡羅仿真結(jié)果，其中DE－PF為微分進(jìn)化粒子濾波算法，RDPF為本文新算法，PF為殘差重采樣粒子濾波算法。由表1可知，本文算法在精度上要明顯優(yōu)于PF算法，相對于DE－PF精度也有所提高，而且在實(shí)時性上，要明顯優(yōu)于DE－PF算法。新算法使用300個粒子可以超過PF算法1000個粒子的精度，并且實(shí)時性更好。

表1 運(yùn)行500次蒙特卡羅仿真的濾波效果對比Table 1 The tracking results of 500 times of Monte Carlo simulation

圖3和圖4分別為3種濾波算法的濾波效果圖和誤差對比圖。

圖3 濾波效果對比圖Fig.3 The tracking results of different filter

圖4 濾波誤差圖Fig.4 The tracking error of filter estimation

由圖可以看出，本文的新算法不僅繼承了微分進(jìn)化重采樣算法在精度方面的優(yōu)點(diǎn)，而且也獲得了殘差重采樣算法的實(shí)時性優(yōu)點(diǎn)，在一定程上產(chǎn)生了更多粒子，抑制了粒子退化現(xiàn)象。新算法具有精度高，實(shí)時性好的雙重優(yōu)點(diǎn)。

4 總結(jié)

本文的RDPF算法將微分進(jìn)化重采樣和殘差重采樣有機(jī)地結(jié)合起來，由于微分進(jìn)化交叉變異操作創(chuàng)造了新的粒子，在一定程度上增加了粒子的多樣性，抑制了粒子退化現(xiàn)象，因此濾波性能有所提高。同時，RDPF算法由于組合了殘差重采樣算法，算法的實(shí)時性得到了一定的提高。仿真結(jié)果表明，本文的算法具有精度高和實(shí)時性好的優(yōu)點(diǎn)，應(yīng)用前景樂觀。

［1］ DOUCET A，DE FREITAS J，GORDON N.Sequential Monte Carlo methods in practice［M］.New York:Springer，2001.

［2］ DOUCET A，GODSILL S，ANDRIEU C.On sequential Monte Carlo sampling methods for Bayesian filtering［J］.Statist Computer，2000(10):197－208.

［3］ 姚智穎，劉冬，劉光斌.基于粒子濾波的INS/TAN組合導(dǎo)航［J］.電光與控制，2006，13(3):54－55.

［4］ 夏克寒，許化龍，張樸睿.粒子濾波的關(guān)鍵技術(shù)及應(yīng)用［J］.電光與控制，2005，12(6):1－4.

［5］ 楊德貴，張長城，鄭濤.粒子濾波算法的性能分析與改進(jìn)［J］.電光與控制，2008，15(5):72－75.

［6］ 葉龍，王京玲，張勤.遺傳重采樣粒子濾波器［J］.自動化學(xué)報(bào)，2007，33(8):885－887.

［7］ STORN R，PRICE K.Differential evolution—a simple and efficient heuristic for global optimization over continuous spaces ［J］.JournalofGlobalOptimization，Kluwer Academic Publishers，1997，11:341－359.

［8］ 馮琦，周德云.基于微分進(jìn)化算法的時間最優(yōu)路徑規(guī)劃［J］.計(jì)算機(jī)工程與應(yīng)用，2005，41(12):74－75.

［9］ STORN R.Sytem design by constraint adaptation and differential evolution［J］.IEEE Transactions on Evolutionary Computation，1999，3(1):22－34.

［10］ 蘇海軍，楊煜普，王宇嘉.微分進(jìn)化算法的研究綜述［J］.系統(tǒng)工程與電子技術(shù)，2008，30(9):1793－1797.

［11］ GORDON N，SALMOND D J，SMITH A F M.Novel approach to nonlinear/non－Gaussian Bayesian state estimation［J］.IEEE Proceedings F Radar and Signal Processing，1993，140(2):107－113.