改進的煙花算法優化粒子濾波研究*

白曉波，邵景峰，田建剛

1.西安工程大學 管理學院，西安 710048

2.陸軍邊海防學院 工程基礎系，西安 710108

1 引言

粒子濾波（particle filter）[1]基于蒙特卡洛思想，對非線性非高斯系統進行狀態估計。但是粒子濾波存在權值退化問題，針對該問題，主要采用重采樣的方法解決[2]，但是重采樣時只保留權值較大粒子，而舍棄小權值粒子，這又導致粒子貧化問題。這就必定影響到系統狀態估計的精度。針對粒子濾波中粒子權值退化和粒子貧化問題，很多學者進行了大量研究，如文獻[3-5]提出了相應的重采樣或改進算法。陳世明等人[6]提出一種基于引力場的粒子濾波算法以避免傳統粒子濾波算法中粒子貧化與退化現象。王向前等人[7]基于粒子濾波和高斯-鄰域搜索濾波的兩個實際特征，提出鄰域迭代重采樣粒子濾波。陳波等人[8]提出了改進的高斯粒子濾波，從而狀態變量服從線性變化，觀測方程為非線性變化的系統模型有顯著效果。這些改進方法都在很大程度上提高了粒子濾波的精度，但是，上述方法仍然基于重采樣方法，未能從根本上解決粒子貧化問題。王法勝等人[9]在粒子濾波的綜述中指出，樣本匱乏或粒子貧化問題依然是該研究領域的熱點。

為了徹底解決權值退化和粒子貧化問題，有學者在粒子濾波中融入了群體智能優化算法。同時，文獻[10]指出基于群體智能算法優化粒子濾波，是粒子濾波新的發展方向。核心思想是，利用群體尋優思想替代重采樣。如文獻[11]根據蟻群算法的基本思想，將樣本粒子隨時間傳遞的過程等效為螞蟻行進路線，避免了重采樣方法引起的樣本匱乏問題。文獻[12]對基于傳統優化算法以及機器學習方法的樣本匱乏解決技術進行了綜述。而比較有代表性的是田夢楚等人[13-14]和朱文超等人[15]利用螢火蟲算法對粒子濾波進行智能優化，以及陳志敏等人[16]利用蝙蝠算法優化粒子濾波。利用螢火蟲和蝙蝠算法的迭代尋優，代替標準粒子濾波的重采樣，使粒子在迭代過程中合理地分布于高似然區和低似然區，保證了粒子的多樣性，且提高了粒子濾波的精度。

然而，螢火蟲算法和蝙蝠算法有其局限性，蝙蝠算法的尋優能力主要依靠蝙蝠個體之間的相互作用和影響，但個體本身缺乏變異機制，一旦受到某個局部極值約束后自身很難擺脫；而且在進化過程中，種群中的超級蝙蝠可能會吸引其他個體迅速向其周圍聚集，使得種群多樣性大幅下降，這必然會在一定程度上影響到粒子多樣性。螢火蟲算法雖然考慮到了個體間的相互影響，但是個體向最優值趨近的過程中，可能會出現螢火蟲移動距離大于個體與最優值間距的情況，導致個體更新自己位置時跳過了最優值，進而最優值發現率降低，影響算法的收斂速度。

綜上，鑒于螢火蟲算法和蝙蝠算法的局限性，結合煙花算法（fireworks algorithm，FWA）[17-18]的優勢，并對其進行改進，以更好地對粒子濾波進行優化有著重要的理論和工程意義。

2 粒子分布多樣性度量

在文獻[13-14]和文獻[16]中，均體現了群體智能啟發式算法對粒子多樣性的促進作用。但是，各文獻并未對算法中粒子多樣性的程度進行度量。因此，這里提出在不同時刻粒子群多樣性的定義。

定義1設粒子濾波過程中，在k時有粒子集合為I={x1k,x2k,…,xNk}，粒子集合I的中心粒子為xck，各粒子xik到中心粒子xck的距離，則粒子多樣性指標表示如下。

以上粒子濾波中粒子多樣性度量的定義是后續研究工作的基礎。

3 改進煙花算法優化粒子濾波

基于螢火蟲算法和蝙蝠算法的局限性，考慮粒子群個體之間的相互作用與影響、個體本身的變異機制，以及算法在迭代過程中出現的粒子個體跳出最優值，出現震蕩，降低濾波精度的問題，先對煙花算法進行改進，再對粒子濾波進行優化。

3.1 標準煙花算法

標準煙花算法[18]由譚營教授等人于2010年提出，算法過程如下。

步驟1隨機初始化N個煙花種群。

步驟2計算煙花Xi的爆炸半徑Ri，產生的火花數Si，計算公式如下。

步驟3產生爆炸火花。計算煙花Xi生成的第j（1≤j≤Si）個火花的第k維坐標公式如下。

步驟4變異算子，以增加種群多樣性。

步驟5映射規則。修正火花坐標，公式如下。

步驟6選擇進入下一次迭代的煙花種群。

3.2 標準FWA不能直接優化PF的原因

標準的煙花算法不能直接用于優化粒子濾波，主要體現在如下兩方面。

（1）火花變異操作。標準FWA的火花變異，采用均值和方差為1的高斯變異，如式（5）所示，雖然能夠在一定程度上促進爆炸火花的隨機性，若直接應用于FWA，促進粒子多樣性的程度較低，對PF的濾波精度提高也極其有限。而PF針對的非線性、非高斯問題，需要粒子具有較好的多樣性，即在狀態值的高、低似然區均分布合理。那么，單一的高斯變異就不能最大程度地滿足實際系統。

（2）火花選擇策略。在文獻[18]中標準FWA的煙花選擇策略需要計算當前個體到候選集合其他個體之間的距離，若候選集合的待選火花數為N，有K個候選集合，那么，選擇策略的時間復雜度就為O(N2)。這樣的時間復雜度若直接用于優化PF，必然在一定程度上降低PF實時性。

3.3 改進煙花算法

對煙花算法進行改進的具體過程如下。

步驟1隨機初始化N個煙花種群。

步驟2利用式（2）、式（3）對每個煙花分別計算Ri和Si。其中為常數，控制爆炸半徑，fmin=min(f(Xi))為當前火花集合中適應度最小值；m為常數；fmax=max(f(Xi))為適應度最大值，f(Xi)為煙花Xi的適應度值，ξ表示計算機中的極小正常數。

為了防止Si過多或過少，對Si進行修正。具體方法如下：

rou(?)為取整函數，a、b為給定常數。

步驟3利用式（4）產生Xi的第j個、第k維坐標，同標準FWA的步驟3。

步驟4混合高斯變異火花。為了增加種群多樣性，通過變異算子產生新的火花，標準FWA中采用高斯函數生成高斯隨機數，但是在PF中，由于系統狀態非線性、噪聲數據非高斯，以致利用高斯函數產生的變異算子不能很好地適應粒子濾波。因此，將標準FWA的高斯函數進一步改進為混合高斯變異算子，提高生成更優良煙花的概率。

設煙花Xi按照預設參數產生M個火花，k=1,2,…,z表示火花xj的第k維坐標。則混合高斯分布下的坐標如下。

式中，f(T)為概率密度函數，c為階段分布系數，N為設定的高斯分布個數，βi、μi和σi分別為混合分布中第i個分布的比例系數、均值和標準差。T（Tmin≤T≤Tmax）為設定煙花隨機爆炸半徑，其分布函數表示如下。

那么，根據分布函數性質得F(Tmax)=1，c(F(Tmax)-F(Tmin))=1，得截斷參數c=1/(F(Tmax)-F(Tmin))，轉換為標準正態分布函數，如式（11）所示。

Φ(?)為標準正態分布函數。

步驟5利用式（6）修正火花坐標。同標準FWA步驟5的映射規則。

步驟6選擇策略。標準煙花算法采用精英選擇策略將適應度最好的煙花留作下一次迭代，具體計算方法詳見文獻[18]，其較復雜的計算過程不適用于粒子濾波，否則必然降低算法的實時性。因此，對選擇方法進行簡化。具體過程如下。

（1）第g代爆炸煙花為，煙花Xi,g的第g+1代火花集合為表示第k維坐標。則第g+1代的候選火花集合為C={xl,c|1≤l≤Si}。基于文獻[19]中的CR公式，改進自適應交叉因子CR的計算方法，對每個執行如下操作。

式中，jrand為[1,2,…,Si]之間的隨機整數，T為迭代次數。CRi為交叉概率，是一條開口向上的拋物線。在種群向下變異的初始階段，CR值較小，且增速較慢，有利于煙花算法的全局搜索能力和多樣性；隨著迭代的進行，CR不斷增大，有利于煙花算法快速收斂，并提高局部搜索能力。進而，有利于PF的粒子在系統狀態的高、低似然區合理分布。

（2）對于候選火花集合C={xl,c|1≤l≤Si}，采用如下方式，用ui,g+1取代Xi,g向下迭代。

pα為隨機選取的火花向下傳遞的概率，rndF(?)表示從集合C中隨機選取一個xl,c。基于改進的選擇策略能夠在一定程度上將大權值粒子直接用于下一次迭代，也考慮了父代的遺傳信息傳遞給子代。且選擇策略的時間復雜度從標準FWA的O(N2)降為O(N)，從而不影響優化后粒子濾波的實時性。

3.4 改進煙花算法的收斂性分析

從改進的煙花算法（improved fireworks algorithm，imp-FWA）整個過程分析，主要分為3個步驟：一是爆炸產生火花；二是煙花混合高斯變異；三是選擇下一代種群。因此，為了分析imp-FWA的收斂性，按這3個步驟依次進行。首先，引入以下定義及性質。

定義2[20]對于方陣A∈Rn×n：

（1）若1≤ ?i,j≤n,aij≥ 0，則稱A為非負的(A≥ 0)。

（2）若1≤ ?i,j≤n,aij>0，則稱A為全正的(A>0)。

（3）若存在一個整數m≥1，使Am>0，則稱A為本原的。

（4）若A非負，且在相同的行、列初等變換后得到如下形式：

且C、T是方陣，則稱A是可約的。

（5）若A非負且，則稱A為隨機的。

（6）若A為隨機陣，且同一列中元素全相同，則稱A是穩定的。

（7）若A為隨機陣，且每列中至少一個正數，則稱A是列容的。

引理1[21]如果矩陣Y是本原的隨機陣，那么Yk收斂到一個全正穩定的隨機陣：

(y1,y2,…,yn)T唯一滿足：

為了研究imp-FWA的馬氏鏈的表示，設所有煙花的坐標空間為S，煙花位置為X，X∈S，馬氏鏈狀態空間為G，則狀態空間維數d=|G|=|S|N。

設imp-FWA每次迭代對應的隨機變量為：R(t)={Xi(t)|1≤i≤n,1≤t＜∞}，表示煙花的位置集合，每次迭代都基于上一步的位置信息，但與原狀態沒有關系。因此，imp-FWA是馬爾可夫過程。為了簡單起見，根據imp-FWA主要的3個步驟，對imp-FWA的收斂性進行分步證明。

（1）爆炸產生的火花。將爆炸火花的爆炸算子的概率矩陣表示為F=(fij)L×L，fij表示爆炸后狀態由i變為j。根據全概率公式得，F為隨機陣。

（2）煙花混合高斯變異。imp-FWA利用混合高斯變異，使得在每一次迭代時，每個煙花都具有一定概率向當前最優位置移動。因此，引入參數ρ＞0，表示個體向當前最優位置移動的概率。則一個群體的狀態由i到j的概率就定義為M=(mi,j)L×L。那么，如果兩個群體狀態i和j有相等的煙花數Ζij，則mij=，N為煙花種群數。從而，M為全正矩陣。

（3）選擇下一代種群。通過選擇算子，煙花狀態發生改變，設sij表示狀態i變為狀態j，則選擇運算對應的概率矩陣為一隨機矩陣，表示為S=(sij)L×L，根據全概率公式得。

定理1通過混合高斯變異改進的imp-FWA能夠以概率1全局收斂到最優解。

證明imp-FWA每一步進化過程的馬氏鏈轉移概率矩陣表示為P=FMS。

設U=FM，V=US。F是隨機矩陣，那么F各行中至少存在一個大于零的元素，因此：

從而U>0。而S亦為隨機陣，同理可得：

由定義2可知P為本原矩陣，再利用引理1可知最優爆炸火花的位置不在馬氏鏈極限分布中的概率為0，且包括最優位置的火花的極限分布之和為1。因此，定理1得證。

3.5 改進的FWA優化PF

基于3.3節改進的煙花算法和標準PF（詳見文獻[9]），提出 FWA-PF（fireworks algorithm particle filter）。核心思想是利用imp-FWA對標準粒子濾波中的每一個粒子進行迭代爆炸尋優，取代PF的重采樣，使粒子向目標狀態后驗概率值高的方向移動，且保持良好的多樣性，以提高狀態估計的準確性。但是，經過煙花算法的處理后，各粒子在狀態空間的分布被改變，其分布密度函數也不能用p(xk|y1,k-1)表示，從而喪失了貝葉斯濾波理論基礎。因此，還需改進粒子權值的表示方法。根據文獻[22]的思想，改進后的粒子權值表示如下：

基于改進煙花算法，優化后的粒子濾波算法過程如圖1所示。

FWA-PF算法的詳細實施過程如下。

步驟1初始化FWA參數，交叉因子CRmin、CRmax，煙花爆炸半徑控制參數，火花數控制參數m，迭代次數初值T=0，ξ為計算機中能夠表示的最小正值，爆炸火花數修正參數a、b，火花被選概率pα。

步驟2初始化粒子集合，并初始化粒子權值ωi，以權值代替式（2）、式（3）中的適應度值。

步驟3對每個粒子Xi進行煙花爆炸變異。

步驟4利用式（15）修正粒子權值。

步驟5式（16）權值歸一化，計算方法如下：

步驟6系統狀態輸出：

3.6 FWA-PF時間復雜度

從FWA-PF的詳細過程來分析，影響算法效率的因素主要有兩個：一是粒子數N；二是每個粒子爆炸所產生的火花數Si。但是Si的取值具有隨機性，由式（3）計算、式（7）修正得出，若式（7）中的控制常數b無限接近1，則Si?m，那么FWA-PF在最壞情況下的時間復雜度為O(3m×N)。而在實際應用時，式（3）中的fmax和f(Xi)為歸一化的粒子權值，則：

式（18）的值向0靠近，那么，若a取值較大，Si＜am，Si的取值大多為rou(am)，FWA-PF的時間復雜度為O(3rou(am)×N)，0＜a＜b＜1。從而在大多數情況下，FWA-PF具有較低的時間復雜度。

3.7 FWA-PF收斂性分析

為了分析FWA-PF的收斂性，首先對FWA-PF過程以及后續分析過程中用到的部分符號進行說明。

（1）狀態過程。FWA-PF依然遵循標準PF的基本特性，即狀態過程X為一階Markov過程，將初始分布和狀態轉移概率分布分別表述為p(dx0)和K(dxt|xt-1)，則動態系統的表述如下。

假設p(dyt|xt)概率密度存在，表示為p(yt|xt)[23]，則：

且p(yt|x0:t)=p(yt|xt)。

（3）FWA-PF新粒子的迭代生成，其實質是利用爆炸算子、混合高斯變異和選擇算子共同作用的結果，也就是每次迭代的粒子按照如下分布進行的分步生成策略：Gh(dxt|xt-1,y1:t)為煙花爆炸混合高斯分布函數，假設存在其概率密度，則用Gh(xt|xt-1,y1:t)表示。

以下是分析過程中用到的符號及說明。‖F‖表示有界函數的最大值，表示大于等于ε的最小整數。Xk:t稱作擴展狀態，表示k時到t時的狀態軌跡，Xk:t={Xk,Xk+1,…,Xt}，小寫字母的xk:t={xk,xk+1,…,xt}或zk:t={zk,zk+1,…,zt}表示其實現。Yk:t為k時到t時的量測軌跡，其實現為yk:t={yk,yk+1,…,yt}。π0:t|m(dx0:t)表示量測為y1:m時，x0:t的后驗概率分布p(dx0:t|y1:m)，π0:t|m(dx0:t)=p(X0:t∈dx0:t|Y1:m=y1:m)。(υ,τ)=∫τυ表示函數υ和τ內積。

表示狀態過程X由0到k時的轉移概率分布，又因為X為一階Markov過程[24]，則：

δ為Dirac-Delta函數。

3.7.1 定義與假設

定義3R(t+1)nx上的函數h(x0:t)，若 (π0:t|t(dx0:t),|h(x0:t)|p)＜∞,p≥1成立，則h(x0:t)屬于Lt,p空間。由測度論[25]定義Lt,p空間的模如下：

假設1若量測軌跡y1:t已知，且t>0時(π0:t|t-1(dx0:t),p(yt|x0:t))>0成立。

假設2假設在FWA-PF中的任意參數μj,t（1≤j≤n,n為參數個數），在t>0時滿足條件。(π0:t|t-1(dx0:t),p(yt|x0:t))≥μj,t>0。實際應用時，雖然難以提前找到符合條件的μi,t，但是若假設1成立，那么總能找到{μi,t|1≤i≤n,t>0}滿足假設2。

假設3t>0時，的模：

引理2假設1和假設3同時滿足，若h(x0:t)∈Lt.p，那么：

引理3若函數h(x0:t)有界，那么h(x0:t)∈Lt,p。

3.7.2 FWA-PF收斂性證明

基于3.7.1節的3個假設，得出如下結論。

引理4若假設1、假設2和假設3都滿足，那么對于?h(x0)∈Lt,4均有與N無關的一組數和，使以下兩式成立：

fw(?)表示獲得N個初始煙花樣本的分布函數，fw(dx0)=π*0N:0|0(dx0:0)。

引理5當假設1、假設2和假設3都滿足，若?H(x0:t-1)∈Lt-1,4均有與N無關的一組數b*t-1|t-1(H(x0:t-1))和η*t-1|t-1(H(x0:t-1))，使以下兩式成立：

那么，對于?h(x0:t)∈Lt,4就具有一組與N無關的數qt|t-1(h(x0:t))和ηt|t-1(h(x0:t))，使以下兩式成立：

引理6假設1、假設2和假設3滿足時，若?H(x0:t)∈Lt,4，均存在數qt|t-1(H(x0:t))和ηt|t-1(H(x0:t))（與N無關），使式（29）、式（30）成立，那么對于?h(x0:t)∈Lt,4均存在與N無關的一組數，使以下兩式成立：

引理7若假設1、假設2和假設3都滿足，且對于?H(x0:t)∈Lt,4均存在與N無關的一組數和，使式（31）和式（32）成立，那么就有數σt（σt與N無關），使得不成立的概率為，并且當。

引理8假設1、假設2和假設3都滿足，且對于均存在與N無關的一組數和，使式（31）和式（32）成立，且當時，對于均存在一組與N無關的數，使得以下兩式成立：

引理9假設1、假設2和假設3都滿足，若，均存在一組與N無關的數和使式（33）和式（34）成立，那么對于均存在與N無關的一組數和，使以下兩式成立：

引理10假設1、假設2和假設3都滿足，若均存在與N無關的一組數和，使式（35）、式（36）成立，則對于Lt,4，均存在一組與N無關的數，使下式成立。

假設4對于FWA-PF，假設其迭代次數有限。

基于引理4至引理10，可得：若假設1、假設2、假設3和假設4都滿足，則對于成立。表示在經驗分布FWA-PF的估計值；表示其最優估計值。下，h(x0:t)的均值作為

3.8 FWA-PF運行機制

FWA-PF的基本思想是利用FWA的煙花爆炸產生的火花，進行混合高斯變異，再通過選擇策略選擇新的火花用作下一次迭代，以取代標準PF的粒子重采樣過程，以實現粒子的多樣性。其運行機制如圖2所示。

在文獻[16]中描述了標準PF在經過多次重采樣后的粒子分布情況，即很多粒子都集中于高似然區，這就喪失了粒子的多樣性。因此，在圖2中，重點表述了FWA-PF中粒子的初始分布、優化后的運動趨勢和優化后的分布情況。以下是對圖2的進一步說明。

（1）圖2（a）是FWA-PF的粒子初始分布，粒子在高低似然區均有分布。

（2）圖2（b）是FWA-PF的各粒子經爆炸、混合高斯變異，以及選擇以后，除兩個方形藍色的粒子在一定概率下直接用作下一代，位置不變外，其他粒子均向各自的方向移動。爆炸的隨機性保證了粒子的多樣性，又由于使用了式（6）的火花坐標修正方法，因此各粒子的分布并不會超過火花爆炸邊界。

（3）圖2（c）是FWA-PF各粒子在圖2（b）所示的趨勢運動后的位置，由于其他粒子的運動，原來紅色的最優粒子離真實值較遠，黑色實心粒子成為最優粒子，而方形藍色粒子位置保持不變。黑色菱形粒子為上一次迭代時最優粒子的新位置。

4 實驗仿真

基于改進煙花算法imp-FWA，在第3章設計了FWA-PF算法，重點從以下四方面說明算法的可行性。

（1）改進煙花算法imp-FWA的BenchMark測試。主要通過進化算法的6個BenchMark標準函數（Sphere、Quadric、Ackley、Rosenbrock、Griewangk、Rastrigin）測試imp-FWA的收斂性和性能。

（2）FWA-PF粒子多樣性評價。FWA-PF基于改進的煙花算法，其性能受初始參數的影響，煙花爆炸半徑控制參數，火花數控制參數m，火花數修正參數a、b。若變異的火花數較少，可保留較多父代的信息，但降低了種群的多樣性，以致尋優效果不好；若變異的火花數較多，可增加種群多樣性，但最優煙花對變異個體產生的貢獻少，導致算法的收斂速度緩慢，以致實時性較低。因此，根據定義1，先對兩個重要參數對粒子多樣性的影響進行分析。

Fig.2 Particle mechanism optimization of FWA-PF圖2 FWA-PF粒子優化機制

M=50為實驗仿真次數；k為濾波次數；j為狀態估計值；xj為量測值。

（4）FWA-PF與ADE-PF[3]、FA-PF[14]和BA-PF[16]做性能和效率的對比分析。

實驗環境、狀態方程和量測方程。

（1）實驗仿真環境。CPU：Intel?CoreTMi5-4200U@160 GHz 2.30 GHz；內存：4 GB，Windows 7 64位操作系統，500 GB硬盤，Matlab R2012a。

（2）狀態方程和量測方程。

Q=10為系統過程噪聲方差，R=10為量測噪聲方差，randn為(0,1)之間的隨機數，粒子數N=200，修正參數a=0.2，b=0.8，即煙花爆炸產生的火花數在[20,80]之間，式（9）中的參數N=3。

4.1 imp-FWA的BenchMark測試

對于進化搜索算法，通常利用BenchMark標準函數測試其性能。這里主要利用6個常用BenchMark函數測試imp-FWA的性能，并與標準FWA進行對比。用于測試的6個BenchMark函數如表1所示。

Table 1 6 BenchMark test functions表1 6個BenchMark測試函數

由圖3的6個BenchMark測試函數的對比得出以下結論：imp-FWA的收斂速度較FWA慢，但是由圖3（b）、圖3（c）和圖3（f）可知，在Ackley、Quadric和Rastrigin函數上，即使到搜索后期，imp-FWA要優于FWA，這得益于改進的選擇策略。而其他圖上的結果在搜索后期FWA和imp-FWA具有相似的收斂效果。尤其在圖3（d）中，由于Rosenbrock函數的信息較少，其搜索方向不明，故imp-FWA和FWA都具有相對穩定的階段。綜上，imp-FWA可以取得較好的收斂性，且由于選擇策略的簡化，使其整個搜索效率提高，能夠用來優化粒子濾波。

4.2 簡化的火花選擇策略對PF精度的影響

Fig.3 Comparison of convergence of imp-FWAand FWA圖3 imp-FWA與FWA收斂性對比

標準FWA的選擇策略是將適應度最好的火花傳遞給下一代，其目的在于使得迭代過程中火花種群能夠快速地向最優值靠近，而改進煙花算法重點目的在于解決PF中粒子權值退化和粒子在迭代過程中喪失多樣性的問題。因此，在改進火花選擇策略時，重點是為了保證PF粒子的多樣性，火花的選擇具有一定概率下的隨機性。同時，在一定概率下將父代的優良火花可以直接傳遞給下一代。如式（12）～式（14）所示。與標準FWA相比，imp-FWA只是在一定概率下將適應度最好的火花傳遞給下一代，導致了一部分適應度好的火花的丟失，在尋優過程的前期收斂性低于標準FWA，但全局收斂性略好于FWA，符合4.1節BenchMark測試結果。若優化PF，粒子多樣性的增加能夠提高濾波精度。另外，FWA-PF通過混合高斯變異算子來增加火花的多樣性。

以下為實驗仿真的方式，對比imp-FWA優化PF（記為PF-1）和標準FWA優化PF（記為PF-2）的濾波效果。如圖4所示。

在圖4中不難發現PF-1在離散點上更接近真實值。imp-FWA和標準FWA優化PF，在不同粒子數和相同方差下（過程方差和量測方差Q=R=20）進行50次獨立測試的濾波性能運行時間對比。如表2所示。

Fig.4 Filter effect comparison of PF-1 and PF-2圖4 PF-1與PF-2濾波效果對比

從表2的結果分析，在相同的粒子數、過程方差和量測方差下，imp-FWA優化后的PF的精度略高于標準FWA優化的PF，且運行效率也高。綜合4.1節和4.2節的實驗結果，得出以下結論。火花選擇策略的簡化，雖然在一定程度上降低了FWA的收斂速度，但是應用于粒子濾波的優化時，卻能保證粒子的多樣性以及父代優良粒子在一定概率下直接傳給下一代。相比將標準FWA直接用來優化PF，imp-FWA優化PF更能提高PF的性能和運行效率。

Table 2 Performance comparison of optimized PF by imp-FWAand by FWA表2 imp-FWA和FWA優化PF性能對比

4.3 FWA-PF粒子多樣性指標

Fig.5 Particle distribution(k=100)圖5 粒子分布(k=100)

在圖5中，k時的真值為15.858 9。從粒子的分布來看，PF的粒子大多集中于高似然區，而FWA-PF的分布在高低似然區，從而FWA-PF比PF的粒子多樣性更好。由式（1）得100個粒子多樣性指標ρ=25.145 3。FWA-PF的粒子多樣性指標ρ=25.400 6。從而也說明了粒子多樣性對濾波精度的促進作用。理論上，煙花爆炸半徑的增大，增加了爆炸產生火花的搜尋范圍，爆炸產生的火花數的增加，有利于增加粒子的多樣性。同時，鑒于實驗的隨機性，對主要兩個參數和m在不同取值時，重復進行50次實驗，計算粒子多樣性指標ρ的均值，分析煙花半徑和爆炸火花數對粒子多樣性指標的影響。具體實驗結果如圖6所示。

Fig.6 Relation of parameter,mand particle diversity index圖6 參數、m與粒子多樣性指標關系

4.4 主要參數取值對FWA-PF性能影響

從理論上分析，粒子多樣性指標越大，粒子的分布就越為合理，且能提高濾波精度，在4.1節中，FWA-PF的估計值為15.553 0，標準PF的估計值為15.470 8，真值為15.858 9。那么，FWA-PF的濾波值更接近真值，濾波精度更高。同樣基于實驗存在的隨機性，利用式（20）計算各參數下FWA-PF的RMSE均值。

Fig.7 Filtering comparison of PF and FWA-PF(m=110,=5)圖7 PF與FWA-PF濾波對比(=5,m=110)

根據圖7的對比結果分析，FWA-PF比標準PF的估計值更接近真值，但是其效果并不明顯。FWA-PF的均方根誤差RMSE=1.826 3，標準PF的均方根誤差RMSE=2.177 7。從而說明煙花爆炸半徑參數不變，增加火花爆炸數，能夠提高粒子濾波的濾波精度，這一結論與圖4所示的實驗結論相吻合，即火花爆炸數的增大能夠增加粒子多樣性指標，而粒子多樣性指標大，則濾波精度也高。

Fig.8 Filtering comparison of FWA-PF and PF(m=120,=5)圖8 FWA-PF與PF濾波對比(=5,m=120)

圖8中，通過對火花數的增加，FWA-PF的估計值比m=110時更接近真值，原因與圖7的實驗相同。標準PF的RMSE為2.245 8，FWA-PF的RMSE為1.629 8。

由于實驗仿真過程中存在隨機性，針對每個火花數控制參數m進行50次獨立實驗，并求得RMSE均值，分析火花數控制參數m對FWA-PF濾波精度的影響，如表3所示。

Table 3 Effect ofmon precision of FWA-PF filtering(=5)表3 m對FWA-PF濾波精度影響(=5)

Table 3 Effect ofmon precision of FWA-PF filtering(=5)表3 m對FWA-PF濾波精度影響(=5)

序號1 2 3 4 5火花爆炸數m100 110 120 130 140 RMSE均值1.928 8 1.877 9 1.750 4 1.725 6 1.706 1

在表3中，在煙花爆炸半徑不變，增加火花爆炸數能夠提高FWA-PF的濾波精度，這和圖6所示的結果一致，即火花爆炸數的增加能夠提高粒子多樣性，而粒子多樣性指標的增加又促進了濾波精度的提高。

其他參數不變，分析煙花爆炸半徑R的增加對FWA-PF精度的影響。

（3）參數m=100，=10。實驗結果如圖9所示。

Fig.9 Filtering comparison of FWA-PF and PF(m=100,=10)圖9 FWA-PF與PF濾波效果對比圖(m=100,=10)

（4）參數m=100，=15。實驗結果如圖10所示。

Fig.10 Filtering comparison of FWA-PF and PF(m=100,=15)圖10 FWA-PF與PF濾波效果對比圖(m=100,=15)

對于FWA-PF，在圖9和圖10的仿真效果上，并未能很好地體現爆炸半徑參數對濾波精度的影響，但是從均方根誤差RMSE來分析，RMSE(R=10)=1.867 0，RMSE(R=15)=1.852 9。即增大煙花爆炸半徑時，FWA-PF的濾波精度有所提高，均高于標準粒子濾波的濾波精度，在圖9和圖10中，RMSE(PF)=1.946 0。鑒于實驗仿真時結果的隨機性，針對不同的爆炸半徑各進行50次仿真實驗，得出實驗結果如表4所示。

在表4中，隨著火花爆炸半徑的增加，RMSE逐漸減小，說明濾波精度增加，即火花爆炸半徑的增加能夠提高FWA-PF濾波精度。

Table 4 Effect ofon precision of FWA-PF filtering(m=100)表4 對FWA-PF濾波精度影響(m=100)

Table 4 Effect ofon precision of FWA-PF filtering(m=100)表4 對FWA-PF濾波精度影響(m=100)

序號1 2 3 4 5火花爆炸半徑 10 15 20 25 30 RMSE均值1.920 3 1.896 8 1.775 2 1.742 8 1.726 7

4.5 FWA-PF與其他算法的對比

為了進一步說明FWA-PF的有效性，利用相同的真實數據，將 FWA-PF（m=120,=20）與 ADE-PF（adaptive differential evolution particle filtering）[3]、FAPF（firefly algorithm optimized particle filter）[14]和 BAPF（bat algorithm optimized particle filter）[16]進行以下對比。

（1）粒子數N=100,K=100,Q=10,R=10。以上4種濾波算法濾波效果對比如圖11所示。

Fig.11 4 algorithms of filtering effect contrast diagram圖11 4種算法濾波效果對比圖

在圖11實驗結果中，4種濾波算法的RMSE均值分別為：RMSE(FWA-PF)=1.650 7，RMSE(BA-PF)=2.180 1，RMSE(FA-PF)=2.190 1，RMSE(ADE-PF)=2.219 0。從RMSE均值可得出以下結論：FWA-PF的濾波精度優于其他3種濾波算法。

（2）不同粒子數算法性能與效率對比。不同粒子數時各進行50次實驗，對比其4種濾波算法的性能和效率。K=200，實驗結果如表5所示。

在表5的實驗結果中，FWA-PF在相同參數條件下，濾波精度均高于其他3種算法。在算法運行效率上，略低于BA-PF，與FA-PF基本相當，但優于ADEPF。通過以上實驗，充分說明了改進的FWA優化粒子濾波的可行性。

4.6 FWA-PF的實施方法

根據4.3節和4.4節的實驗結論，FWA-PF中的粒子多樣性和濾波精度均受火花爆炸數和爆炸半徑的影響，即粒子多樣性和濾波精度與煙花爆炸半徑和火花爆炸數都成正相關關系，但是，火花爆炸數的增加必然引起算法時間復雜度增加，并降低FWA-PF的實時性。那么，在解決實際問題時，取得合適的參數就顯得尤為重要。而對于FWA-PF，對其性能的要求，主要體現在兩方面，一是濾波精度，二是運行效率。因此，下面結合這兩個因素給出FWA-PF中爆炸半徑R和火花爆炸數S的設置過程如下。

步驟1在特定軟硬件環境下（操作系統、處理器性能、內存和硬盤容量），以及基本的粒子數、狀態過程方差和量測過程方差，設定合理的濾波精度閾值p（常以RMSE均值表示）和運行時間閾值e。

步驟2設置初始火花爆炸半徑R和爆炸數S。

步驟3運行程序，經過多次采樣（一般不少于50次），計算均方根誤差均值RMSEavg和運行時間均值tavg。在具體調整參數時又分為以下4種情況。

Table 5 Results of experimental simulation表5 實驗仿真結果

（1）若RMSEavg≤p，且tavg≤e，說明系統達到了濾波精度和實時性要求，不需要再進行調整。進入步驟4，結束調整。

（2）若RMSEavg>p，但tavg≤e，說明實時性滿足要求，而精度不足。那么就優先考慮以步長l增加煙花爆炸半徑，其次以步長s增加爆炸粒子數，兩者同時增加，重復步驟3。

（3）若RMSEavg≤p，但tavg>e，說明濾波精度達到要求，而運行效率不足。對于這種情況，優先考慮以步長s減少爆炸粒子數，重復步驟3。

（4）若RMSEavg>p，但tavg>e，說明濾波精度和運行效率均未達到要求，那么就優先考慮以步長l增加煙花爆炸半徑，并以步長s減少火花爆炸數，重復步驟3。

步驟4結束調整，設定參數火花爆炸半徑R和火花爆炸數S。

通過以上4個步驟，實現FWA-PF兩個重要參數的設定。

5 結束語

針對粒子濾波中的粒子退化和喪失多樣性問題，將煙花算法（FWA）的火花變異算法改為混合變異算法，以此提高FWA在粒子濾波中對不同分布的適應性，從而保證了粒子的多樣性。通過實驗仿真的方法，從濾波的精度和程序運行效率兩方面，驗證了利用改進的FWA優化標準粒子濾波的可行性。然而，FWA算法的性能受制于兩個主要參數：煙花爆炸半徑和爆炸火花數。這兩個參數影響到粒子的多樣性和濾波性能。因此，針對實際應用時的狀態方程和量測方程，如何自適應地選擇最合適的FWA參數是下一步主要的研究內容。