盧敏 ， 陳菘 ， 張敏 ,2

（1. 江西理工大學(xué)理學(xué)院，江西 贛州 341000；2. 嘉興學(xué)院數(shù)理與信息工程學(xué)院，浙江 嘉興 314001）

0 引 言

二戰(zhàn)時期，Von Neumann 提出了蒙特卡羅方法(Monte Carlo Method)，該方法的主要思想是從概率分布中隨機抽樣并用這些抽樣值來近似表示原分布[1]。 貝葉斯估計結(jié)合先驗概率函數(shù)與觀測似然函數(shù)，利用貝葉斯公式來求解估計問題[2]。粒子濾波算法使用遞歸貝葉斯估計和蒙特卡羅方法實現(xiàn)了對估計問題的求解，將系統(tǒng)狀態(tài)的處理轉(zhuǎn)化成序列的估計問題[3]。 由于粒子濾波算法可以處理任意狀態(tài)空間描述的模型[4]，因此獲得了廣泛的應(yīng)用。除用于估計問題外，粒子濾波算法在視覺跟蹤、航空導(dǎo)航、通信與信號處理、手勢識別、故障檢查等[5-8]領(lǐng)域也發(fā)揮著重要作用。

粒子濾波算法需要從后驗分布中抽樣，并用這些抽樣值求解后驗分布的期望、 方差等數(shù)字特征。但通常情況是不知道分布函數(shù)的具體解析式，無法直接采樣，為此引入了提議分布。 提議分布形式已知且與分布函數(shù)近似，對提議分布的采樣間接實現(xiàn)了對后驗分布的采樣。最優(yōu)提議分布融入了當(dāng)前觀測數(shù)據(jù)，但實際采樣中無法實現(xiàn)，因此用狀態(tài)轉(zhuǎn)移函數(shù)作為次優(yōu)提議分布。次優(yōu)提議分布導(dǎo)致了粒子退化問題（Particles Degeneracy）[9]，當(dāng)似然函數(shù)位于先驗概率分布尾部時，退化問題更為嚴重。 粒子退化現(xiàn)象產(chǎn)生了大量低權(quán)重粒子，低權(quán)重粒子無助于后驗分布的估計，浪費了大量的計算時間。 退化現(xiàn)象的存在使得少量粒子無法保證估計的精度，需要采集更多粒子實現(xiàn)對狀態(tài)的估計。 Douceta 通過數(shù)學(xué)方法證明了若干輪迭代后粒子退化現(xiàn)象發(fā)生的必然性[10]，在自舉濾波器[11]（Bootstrap Filter）提出以后，退化問題得以解決。 自舉濾波器根據(jù)權(quán)值大小復(fù)制粒子，用大權(quán)值粒子替換小權(quán)值粒子。 這種機制導(dǎo)致粒子種類銳減，出現(xiàn)了粒子多樣性缺失的新問題。

解決粒子多樣性缺失的方法有改進重采樣算法或改進提議分布[12]。 蔡登禹等用遺傳算法替代了粒子濾波的重采樣過程，結(jié)合變異、交叉等步驟，提高了濾波的精度[13]。 但是這種方法參數(shù)較多，參數(shù)的設(shè)置有很強的隨機性。張威虎等將粒子群優(yōu)化算法引入到粒子濾波過程中，通過蝴蝶個體香味濃度的高低和吸引半徑引導(dǎo)種群向最優(yōu)個體附近移動[14]，該算法運行時間較短，跟蹤誤差較小，但存在個體被多次吸引和移動的問題。本文將螢火蟲算法引入到粒子濾波中，把粒子比作螢火蟲個體，粒子的權(quán)值比作螢火蟲發(fā)光的亮度。 為避免文獻[14]粒子在最優(yōu)個體值附近震蕩的問題，本文引入遞減函數(shù)[15]，該函數(shù)使吸引度隨著迭代次數(shù)的增加而不斷減小。由于螢火蟲算法中每個螢火蟲都要與其他螢火蟲個體進行對比，亮度低于其他個體時均需要進行移動尋優(yōu)操作，當(dāng)種群數(shù)目過多時影響算法的運行效率與收斂速度。 為減少算法的復(fù)雜度，本文利用最優(yōu)鄰居引導(dǎo)螢火蟲個體的移動。設(shè)置半徑參數(shù)將螢火蟲種群進行劃分，個體向所屬種群中亮度最高的螢火蟲移動完成位置更新，不同種群的個體再向最優(yōu)種群的最亮個體移動并分布在最亮個體周圍。此操作加快了迭代尋優(yōu)的速度并保證了螢火蟲種群的多樣性，種群的多樣性提高了濾波的準確度。 最后實驗驗證本文改進的算法在跟蹤誤差、跟蹤精度方面優(yōu)于粒子濾波算法。

1 粒子濾波模型

1.1 模型分析

非線性系統(tǒng)的狀態(tài)模型可以用式（1）和式（2）表示：

其中，xk與zk分別表示k 時刻的狀態(tài)量與觀測量；ω 與v 分別表示過程噪聲和測量噪聲；f 與h 分別表示狀態(tài)轉(zhuǎn)移函數(shù)與測量函數(shù)。 圖1 是式（1）和式（2）的狀態(tài)空間模型。

圖1 狀態(tài)空間模型

由圖1 可以看出： ①狀態(tài) xk與狀態(tài)xk-1相關(guān)，與 xk-2，xk-3，…等狀態(tài)不相關(guān)。 ②zk與狀態(tài) xk有關(guān)，與 zk-1，zk-2，…不相關(guān)，與 xk-1，xk-2，…也不相關(guān)。

在 k-1 時刻，p（x0∶k-1｜z1∶k-1）已知，根據(jù)式（1），并利用 Chapman-Kolmogorov 方程，由 p（xk-1｜z1∶k-1）得到 p（xk｜z1∶k-1）：

當(dāng)最新時刻觀測值到來時，利用最新時刻的觀測值 zk修正 p（xk｜z1∶k-1）。

其中，p（zk｜z1∶k-1）為歸一化常數(shù)。

引入提議分布后，權(quán)重更新公式為：

通過從提議分布中采樣獲得樣本， 并根據(jù)式（5）計算樣本權(quán)重即實現(xiàn)了粒子濾波算法。圖2 為粒子濾波過程的示意圖。

圖2 粒子濾波過程

1.2 重采樣算法

重采樣基于逆變換的思想，粒子歸一化權(quán)值與總權(quán)值的比例關(guān)系決定了粒子被復(fù)制的概率。目前常采用的重采樣算法有多項式重采樣、殘差重采樣和系統(tǒng)重采樣[16]。 用式（6）衡量有效粒子數(shù)退化的程度：

其中,vap（w）是所有粒子權(quán)重的方差。

式（6）計算復(fù)雜，通常采用另一種形式來近似表示：

其中，wi是粒子的歸一化權(quán)值。 完全退化情況下權(quán)重值完全集中在某個粒子上，此時Neff值為1。最佳情況是所有粒子權(quán)重值相等 （重采樣后所有粒子權(quán)重值相等，為 1/N），此時 Neff=N。 使用式（7）來衡量粒子退化的程度，當(dāng)有效粒子數(shù)小于設(shè)定的閾值時需要進行重采樣。 重采樣過程如下：首先產(chǎn)生一個[0，1]之間均勻分布的隨機數(shù)ui，然后計算粒子的累計權(quán)值概率之和時，則粒子xm被選中復(fù)制一次，重復(fù)N 次獲得所需的樣本數(shù)。

2 群智能優(yōu)化理論

2.1 螢火蟲算法

群智能算法 （Swarm Intelligence Algorithm）屬于自然啟發(fā)式算法，可以更好地解決大規(guī)模復(fù)雜非線性問題。 該算法主要模擬了昆蟲、鳥群等動物的聚集行為，群體中的成員通過協(xié)同操作不斷優(yōu)化搜尋的方向從而達到求最優(yōu)解的目的。目前常見的群智能算法有蟻群算法、蛙跳算法和螢火蟲算法等。

螢火蟲算法有2 個版本，firefly 算法和glowworm算法[17]。 效果類似、功能相近，本文采用 firefly 算法。 在firefly 算法中，螢火蟲依靠自身所發(fā)光的亮度與對其他螢火蟲個體的吸引度大小完成位置更新。 亮度與吸引度均與距離有關(guān)， 隨著距離的增加亮度與吸引度均減小。 位置更新的原則是亮度低的螢火蟲會被亮度高的螢火蟲吸引，亮度公式定義如下：

其中，L0表示個體自身的亮度；γ∈[0.1，2]為光強系數(shù)，γ；rij表示距離，公式如下：

螢火蟲之間的吸引度βij定義為：

其中，β0為螢火蟲的最大吸引度，通常取[0.8，1]。

位置更新公式定義如下：

其中，rand 是介于 0 到 1 間的隨機數(shù)；α 是步長因子；β 是式（10）中計算得到的相對吸引度。

2.2 改進算法

由式（8）和式（10）可知，螢火蟲的亮度和吸引度大小與距離負相關(guān)。迭代前期螢火蟲個體間距離遠，吸引度小，導(dǎo)致位置更新步驟中移動的距離短，需要進行多次迭代才能移動到最優(yōu)個體附近。迭代后期螢火蟲個體間距離近，相互間的吸引度增大導(dǎo)致每次移動的距離大，容易出現(xiàn)在最優(yōu)值附近反復(fù)震蕩的問題。 另一方面，粒子濾波過程中次優(yōu)提議分布 p（xk｜xk-1）忽略了最新時刻的觀測值數(shù)據(jù)，導(dǎo)致了粒子退化現(xiàn)象的產(chǎn)生。因此本文主要從以下三個方面對螢火蟲算法進行改進：

1） 在螢火蟲相對亮度的計算公式中引入最新的觀測值數(shù)據(jù)。 修改后的亮度計算公式如下：

其中，R 表示測量噪聲方差；zpre表示由式（1）得到的預(yù)測值；z 表示由式（2）得到的觀測值。 L 值越大說明該螢火蟲個體亮度高，位置優(yōu)，對應(yīng)的粒子權(quán)值大，越接近真實的位置。

2）在螢火蟲吸引公式中引入遞減函數(shù)δ，該函數(shù)控制吸引度的更新。 迭代初期δ 值應(yīng)比較大，加強算法的全局搜索能力，有利于低亮度的螢火蟲個體迅速收斂到高亮度個體附近。迭代后期δ 值應(yīng)比較小，提高算法的局部搜索能力，避免低亮度螢火蟲個體在高亮度螢火蟲個體附近反復(fù)震蕩，無法收斂。 基于上述思路，δ 函數(shù)設(shè)置如下所示：

其中，k 為迭代次數(shù)， 一般設(shè)置為30 到100 之間，改進后的吸引度公式如式（14）所示：

3）當(dāng)螢火蟲種群規(guī)模N 很大時，個體在迭代尋優(yōu)的過程中需要與其他螢火蟲個體兩兩比較然后進行移動。這個比較過程增加了運算的復(fù)雜度，為此本文使用最優(yōu)鄰居引導(dǎo)螢火蟲個體的移動，有效減少了移動的次數(shù)。 用圖3 中的圓形表示每個螢火蟲個體，空心圓形代表比中心螢火蟲亮度低的螢火蟲個體，實心圓形代表比中心螢火蟲亮度高的螢火蟲個體，黑色實心圓形是亮度最高的螢火蟲個體。

圖3 最優(yōu)鄰居吸引模型

首先需要確定中心螢火蟲的半徑，計算過程如下：

式（15）通過計算中心螢火蟲與其余螢火蟲距離的平均值得到螢火蟲i 的半徑，然后計算所有螢火蟲個體 j（j≠i）與螢火蟲 i 的距離 dij。 如果 dij≤rand×diave，那么 j 確定為中心螢火蟲 i 的鄰居，不滿足式 （15）， 則該螢火蟲個體不是中心螢火蟲的鄰居。 以下分兩種情況討論：

情況1：螢火蟲i 不存在鄰居時，即中心螢火蟲為遠離群體的獨立個體，此時i 直接向全局最優(yōu)處的螢火蟲個體進行移動。 公式如下：

式（16）中： β 為吸引度大小，可由式（14）計算得到；Bk為全局最優(yōu)螢火蟲個體的位置。

情況2：螢火蟲i 存在鄰居時，由式（12）分別確定螢火蟲i 和其鄰居螢火蟲的亮度，并找出亮度最優(yōu)的個體。 若螢火蟲i 比所有鄰居螢火蟲亮度都高，說明該螢火蟲個體權(quán)重大，此輪過程不進行移動。 若鄰居螢火蟲存在最亮個體，則螢火蟲i 依據(jù)式（17）向最優(yōu)鄰居進行移動：

式（17）中：xk為鄰居中最亮螢火蟲個體的位置。

綜上，本文改進后的粒子濾波算法流程如下：

3 算法仿真

3.1 非線性方程模型及參數(shù)設(shè)置

實驗環(huán)境如下：硬件配置為Intel i7 處理器，8 GB內(nèi)存，軟件配置為Matlab R2016b。

本實驗相關(guān)參數(shù)設(shè)置如下：粒子數(shù)目為60，跟蹤時長T 為50，步長因子α 為常數(shù)，本文取0.3。最大吸引度β0為0.9，迭代次數(shù)k 為50。 初始種群個體滿足均值為0.1，方差為2 的高斯分布，量測噪聲與過程噪聲服從均值為0， 方差為1 的高斯分布。實驗中所采用的非線性狀態(tài)模型與量測模型來自文獻[18]，如式（18）和式（19）所示：

式（18）和式（19）是一個非線性非高斯方程，后驗分布具有多個峰值。 x（k）與 z（k）分別表示 k 時刻的狀態(tài)值與量測值，wk與vk是均值為0，方差分別為R=1 與Q=1 的高斯白噪聲。

3.2 實驗結(jié)果分析

圖4、圖5 分別是兩種算法實驗誤差的對比與跟蹤效果。

圖4 誤差分布

圖5 狀態(tài)估計

分析圖4 可得：在粒子數(shù)目與跟蹤時長均相同的條件下，相比于粒子濾波算法，本文改進算法誤差分布的峰值得到了有效降低，大誤差現(xiàn)象得到了一定程度的改善。這主要是改進算法在亮度公式的計算中利用了最新時刻觀測值，相比于粒子濾波算法的次優(yōu)提議分布誤差得到了減小。從圖5 可以看出，本文改進算法跟蹤效果有了顯著提高，大部分時刻都圍繞在真實值附近，而粒子濾波算法估計值與真實值之間出現(xiàn)了較大的偏差。這主要是本文改進算法不對低權(quán)重粒子進行舍棄，低權(quán)重粒子不斷向高權(quán)重粒子移動從而提高了跟蹤的準確性。

為了定量評價本文改進算法的性能，采用均方根誤差（Root Mean Square Error, RMSE）作為標準計算真實值與預(yù)測值的偏差。均方根誤差RSME 公式如下所示：

取過程噪聲方差R 分別為1 和0.1，粒子數(shù)目分別為50 和100， 其余參數(shù)詳見3.1 節(jié)參數(shù)設(shè)置。為了排除實驗的隨機性，減少偶然誤差，表1 的數(shù)據(jù)是在進行100 次實驗的基礎(chǔ)上取平均值得到的。

表1 均方根誤差數(shù)據(jù)

分析表1，當(dāng)噪聲方差相同時，隨著采樣粒子數(shù)目的增加，粒子濾波算法與本文改進算法的均方根誤差均相應(yīng)減小。 以R=0.1 為例，粒子濾波算法的誤差隨著采樣數(shù)目的增加減小了6.7%， 本文改進算法的誤差減小了9.2%。 增加采樣粒子數(shù)目可以減小狀態(tài)估計的誤差，但與之伴隨的是計算時間的增加，影響了算法的時效性。 在噪聲與粒子數(shù)目均相同的條件下，本文改進算法誤差小于粒子濾波算法誤差，驗證了改進算法的性能。分析原因，粒子濾波重采樣過程對低權(quán)重粒子直接舍棄，高權(quán)重粒子的大量復(fù)制降低了有效粒子數(shù)目。而本文改進算法對低權(quán)重粒子進行移動尋優(yōu)操作，使其分布在真實后驗概率分布值附近，這個過程中沒有舍棄低權(quán)重粒子，降低了估計狀態(tài)與真實狀態(tài)之間的誤差。

3.3 二維目標跟蹤模型

對于平面內(nèi)做勻速運動的載體， 可以使用位置、速度信息來描述其運動狀態(tài)。本文使用的狀態(tài)變量為X（k）=[xp（k），xv（k），yp（k），yv（k）]T。 其中，xp（k）、 xv（k）為載體水平方向的位置與速度，yp（k）、yv（k）為載體豎直方向的位置與速度。載體的狀態(tài)運動方程如式（21）所示：

其中，A 為狀態(tài)轉(zhuǎn)移矩陣，W（k）為系統(tǒng)模型噪聲。

其中，ax和 ay是服從的高斯白噪聲。

觀測站可以使用激光、無線電、紅外等儀器對運動載體的狀態(tài)X（k）進行探測，探測得到的距離用 Z（k）來表示，探測噪聲用 V（k）來表示。 系統(tǒng)的觀測方程如式（22）所示：

3.4 仿真結(jié)果分析

圖6、圖7 是兩種算法的跟蹤效果圖與誤差對比，圖8 是時間消耗對比，圖9 是方差為1 的高斯噪聲。

圖6 跟蹤效果對比

圖7 誤差對比

圖8 時間對比

圖9 方差為1 時的高斯噪聲

圖6、圖7 是對二維平面內(nèi)的運動載體跟蹤了100 次所得出的。 對比來看，本文改進算法的跟蹤效果優(yōu)于粒子濾波算法。 隨著跟蹤時長的增加，粒子濾波算法誤差不斷增加， 而本文算法跟蹤誤差小于粒子濾波算法， 在20 到30 時刻之間誤差改善效果更為明顯。 但另一方面，本文算法的運算時間明顯大于粒子濾波算法。 如圖8 所示，本文算法的時間消耗是粒子濾波的2～3 倍。 為了減小跟蹤的誤差， 這種以增加運算時間來換取跟蹤準確性的操作是必然的。 重采樣算法只是簡單地選中并復(fù)制粒子，造成了粒子種類的減少。 本文算法不對粒子進行選擇與復(fù)制， 通過粒子的移動來保證粒子種類的多樣性， 因此提高了跟蹤準確性并減小了跟蹤誤差。

4 結(jié) 論

本文將螢火蟲算法引入到粒子濾波過程中，改進了螢火蟲算法中相對亮度與吸引度的計算方法，在相對亮度的計算過程中利用了當(dāng)前時刻的觀測值數(shù)據(jù)。 在吸引度的計算中引入遞減函數(shù)，粒子在向最優(yōu)值靠近的過程中避免了在最優(yōu)值附近震蕩。 利用最優(yōu)個體指引螢火蟲的移動，降低了算法的復(fù)雜度， 最后通過仿真實驗驗證改進算法的性能。 結(jié)果表明：

1） 粒子濾波算法提議分布的選取沒有融合當(dāng)前時刻的觀測值，當(dāng)似然函數(shù)遠離提議分布的峰值時，無法保證狀態(tài)估計的準確性。 本文改進算法融合了當(dāng)前時刻的觀測值，提高了估計的精度。

2） 重采樣過程用大權(quán)重粒子替代小權(quán)重粒子，解決了粒子退化問題，但是造成了樣本的匱乏。本文改進算法中小權(quán)重粒子通過一定的規(guī)則向大權(quán)重粒子移動， 從而分布在大權(quán)重粒子的周圍，保證了樣本的多樣性，從而減小了估計的誤差。

3） 通過對非線性過程的仿真得出本文改進算法可以保證估計的精度， 具有一定的實用價值，下一步將研究改進算法在智能跟蹤、導(dǎo)航等場景的具體應(yīng)用。