異常值和未知觀測噪聲魯棒的非線性濾波器

方安然，李旦,3,*，張建秋

1.復(fù)旦大學(xué) 智慧網(wǎng)絡(luò)系統(tǒng)研究中心，上海 200433 2.復(fù)旦大學(xué) 電子工程系，上海 200433 3.上海市空間智能控制技術(shù)重點實驗室，上海 201101

在雷達、聲吶、通信和語音識別等應(yīng)用中，觀測中隨機出現(xiàn)異常值是一種十分常見的現(xiàn)象[1-2]。例如，無線通信中，電路通斷暫態(tài)過程產(chǎn)生的脈沖干擾；在雷達系統(tǒng)或聲吶系統(tǒng)中，人為或自然因素產(chǎn)生的沖擊性干擾，而引起雷達或聲納觀測的隨機異常波動[3]等。文獻[4]表明：這種觀測的隨機異常波動是一種非高斯噪聲，其概率密度分布有一個比高斯分布更厚重的“尾部”，因此通常稱其為長尾分布。傳統(tǒng)非線性系統(tǒng)的濾波方法一般聚焦于非線性過程的近似描述，如：擴展卡爾曼濾波(Extended Kalman Filter, EKF)[5]，就是利用一階泰勒展開對非線性方程近似線性化，然而，由于它的線性化誤差較大，因此就造成濾波性能不佳；迭代擴展卡爾曼濾波(Iterated Extended Kalman Filter, IEKF)[6]對EKF的改進，就是通過迭代的逐步逼近，以減小線性化誤差，卻也因此造成了運算復(fù)雜度的上升；利用Sigma點方法逼近非線性過程，在避免迭代計算的同時降低了狀態(tài)估計誤差，這一類算法包括利用無跡變換(Unscented Transform, UT)的無跡卡爾曼濾波(Unscented Kalman Filter, UKF)[7]、利用球面容積積分的容積卡爾曼濾波(Cubature Kalman Filter, CKF)[8]和利用高斯-埃爾米特積分的正交卡爾曼濾波(Quadrature Kalman Filter, QKF)[9]等。但是，以上非線性濾波方法，都是假設(shè)觀測噪聲為高斯白噪聲。因此在含異常值的噪聲環(huán)境中，即存在隨機觀測異常值或非高斯的噪聲環(huán)境中，這些算法的性能都會顯著下降，甚至失效[5]。

針對非高斯非線性系統(tǒng)的濾波，傳統(tǒng)算法主要有以高斯混合模型(Gaussian Mixture Model, GMM)來逼近非高斯分布的高斯和濾波(Gaussian Sum Filter, GSF)[10]，以及利用蒙特卡洛模擬近似描述非高斯分布的粒子濾波(Particle Filter, PF)[11]等兩種。盡管這兩種算法都能處理含異常值或非高斯噪聲環(huán)境中的動態(tài)濾波問題，但它們都需要有非高斯噪聲的準確先驗知識，以及都存在運算復(fù)雜度偏高等問題[5]。為了解決這些問題，文獻[2]提出了對異常值(離群點)魯棒的卡爾曼濾波(Outlier-Robust Kalman Filter, ORKF)。它通過對卡爾曼濾波算法進行修正，以滿足非高斯非線性系統(tǒng)濾波的要求。不過由于其對先驗噪聲協(xié)方差準確與否十分敏感，因此魯棒性欠佳且應(yīng)用受限。為了進一步提高濾波算法的魯棒性，文獻[12]以Kullback-Leibler散度(Kullback-Leibler Divergence, KLD)為目標(biāo)函數(shù)，對非線性觀測方程進行線性化，提出了一種后驗線性化濾波(Posterior Linearization Filtering, PLF)算法。盡管該算法對噪聲先驗信息的準確與否并不敏感，但其估計性能欠佳。針對這一問題，文獻[13]提出了最大相關(guān)熵的無跡卡爾曼濾波器(Maximum Correntropy Unscented Kalman Filter, MCUKF)，改善了狀態(tài)估計的性能，不過其收斂較慢，且高度依賴濾波參數(shù)的選擇是否合適。不幸的是：如何選擇濾波參數(shù)，目前還沒有令人信服的解決辦法[14]。為此，以Huber函數(shù)[15]作為損失函數(shù)，文獻[14]又提出了魯棒微分無跡卡爾曼濾波器(Robust Derivative Unscented Kalman Filter，RDUKF)，這是因為基于Huber函數(shù)的濾波參數(shù)選擇，目前已有較成熟的方法[14]。然而，RDUKF要求動態(tài)系統(tǒng)的狀態(tài)轉(zhuǎn)移方程為線性，因此不能應(yīng)用于非線性系統(tǒng)的濾波。文獻[16-18]則報道了基于Huber函數(shù)的另一類濾波算法，稱為基于Huber的無跡濾波器(Huber-based Unscented Filter，HUF)，盡管它們適用于非線性系統(tǒng)的濾波，但必須有觀測噪聲二階統(tǒng)計量的先驗。為了克服這一缺點，文獻[19]報道了一種R-自適應(yīng)無跡卡爾曼濾波器(R-Adaptive Kalman Filter, RUKF)，不幸的是，該方法一旦受到異常值的干擾，則極易發(fā)散[20]。

針對含異常觀測值非線性系統(tǒng)的濾波，本文提出了一種新的濾波算法，稱之為對異常值魯棒的后驗線性化濾波器(Outlier-Robust Posterior Linearization Filter, ORPLF)。分析表明：以Huber損失函數(shù)代替最大后驗(Maximum a Posterior, MAP)準則中加權(quán)觀測誤差的l2范數(shù)，就得到了一個新的最優(yōu)化準則函數(shù)。由于Huber損失函數(shù)同時具有l(wèi)1和l2范數(shù)的性質(zhì)，因此借助于這個新的最優(yōu)化準則函數(shù)，本文推導(dǎo)的濾波器就兼具了l1范數(shù)對異常值的魯棒性，以及l(fā)2范數(shù)對函數(shù)擬合的精確性。當(dāng)觀測噪聲的分布未知時，則可通過箱線圖法[21]估計其分布的方差，這樣就進一步提出了對異常值和未知觀測噪聲魯棒的后驗線性化濾波器(outliers and unknown Observation Noise-Robust Posterior Linearization Filter, ONRPLF)。仿真實驗驗證了分析結(jié)果的有效性，它們也表明：本文算法的濾波性能(包括誤差和魯棒性)在含異常值的(未知)觀測噪聲中優(yōu)于現(xiàn)有文獻報道的非線性濾波類算法。

本文內(nèi)容組織如下：第1節(jié)回顧了PLF算法，并給出了利用Huber損失函數(shù)的MAP準則，進而利用該準則推導(dǎo)了適用于非線性系統(tǒng)的ORPLF算法；第2節(jié)給出了Huber函數(shù)中調(diào)諧參數(shù)μ的確定方法，也給出了利用箱線圖法估計噪聲分布方差的方法；第3節(jié)給出了目標(biāo)跟蹤的仿真實驗；第4節(jié)總結(jié)了全文。

1 異常值魯棒的非線性濾波器

1.1 后驗線性化濾波器

假設(shè)一離散非線性系統(tǒng)的狀態(tài)空間模型為[20]

xk=fk-1(xk-1)+wk-1

(1)

yk=hk(xk)+vk

(2)

(3)

式中：Ak∈Rd×n為線性化的觀測矩陣；bk∈Rd為線性化的偏差；ek∈Rd和ek～N(ek|0,Ωk)表示線性化時的誤差；Ωk為線性化誤差ek的協(xié)方差；Ak、bk、Ωk為在時刻k對hk(xk)線性化的參數(shù)。

文獻[12]給出了一種通過最小化KLD來進行非線性觀測方程線性化的方法，其KLD定義為

(4)

由文獻[22]知，KLD反映的是兩個分布的匹配程度。匹配程度越高，KLD就越小。若兩個分布完全匹配，則它們之間的KLD就為0。據(jù)此，數(shù)學(xué)上最小化KLD就可表示為[12]

p(xk|yk,Ak,bk,Ωk))

(5)

當(dāng)狀態(tài)xk的分布p(xk)和非線性方程hk(·)都已知時，利用hk(·)關(guān)于p(xk)的統(tǒng)計線性回歸(Statistical Linear Regression, SLR)，可得到式(5)的解為[12]

(6)

(7)

(8)

(9)

(10)

(11)

式(9)～式(11)中的積分可以通過不同的數(shù)值方法進行逼近[5]，例如：無跡變換[7]、球面容積積分[8]、高斯-埃爾米特積分[9]等。這樣，近似線性化后的新觀測方程就為

yk=Akxk+bk+rk

(12)

式中：rk=ek+vk為一個含異常值的觀測噪聲。

由于ek與vk是相互獨立的，而vk去除異常值之后的協(xié)方差為Rk且ek～N(ek|0,Ωk)，因此rk去除異常值之后的協(xié)方差就為(Ωk+Rk)。這樣，將含異常值的非線性觀測方程，變換成了含異常值的線性觀測方程。若采用廣義高斯濾波[5]的方法進行狀態(tài)預(yù)測，那么據(jù)線性化之后的觀測方程(12)，執(zhí)行卡爾曼濾波的更新步驟[5]，就可得到文獻[12]的PLF算法。

1.2 異常值魯棒的后驗線性化濾波器

從最大后驗(Maximum a Posterior, MAP)的角度，在高斯噪聲環(huán)境中PLF算法的狀態(tài)更新就是最大化如下函數(shù)[12]

(13)

對式(13)取負對數(shù)，則MAP就等價于最小化如下?lián)p失函數(shù)：

(14)

從式(14)中，可以發(fā)現(xiàn)它采用的是l2范數(shù)。由文獻[15]知，l2范數(shù)的預(yù)測值距離真實值越遠，則其懲罰力度越大，這就意味著它對異常值比較敏感。也就是說，l2范數(shù)容易受到隨機異常值的干擾，甚至有可能會導(dǎo)致算法失效。因此，l2范數(shù)不適用于含有隨機異常觀測值系統(tǒng)的濾波。由文獻[15]知，l1范數(shù)相較于l2范數(shù)對異常值具有更高的魯棒性?？墒?，l1范數(shù)可能存在不可導(dǎo)的奇異點，這為最小化l1范數(shù)的計算帶來了困難。為此，本文期望通過引入Huber損失函數(shù)，在降低異常值對濾波干擾的同時，又可保證等效l1范數(shù)處處可導(dǎo)。

Huber損失函數(shù)的定義為[15]

(15)

(16)

υk=(Ωk+Rk)-1/2[yk-(Akxk+bk)]

(17)

式中：A1/2為對稱正定矩陣A的Cholesky分解；AT/2為A1/2的轉(zhuǎn)置，滿足A=A1/2AT/2,A-1/2是A1/2的逆矩陣，A-T/2是AT/2的逆矩陣。

利用式(15)的Huber損失函數(shù)，將式(14)改寫為

(18)

對式(18)中的xk求導(dǎo)并令該導(dǎo)數(shù)等于零，有

(19)

(20)

定義矩陣：

(21)

那么由式(17)，有

(22)

將式(20)～式(22)代入式(19)，并利用文獻[15,18]中的矩陣恒等式，有

(23)

再將式(17)代入式(23)，可得

Γk(Ωk+Rk)-1/2[yk-(Akxk+bk)]=0

(24)

(25)

整理得

(26)

利用如下矩陣求逆公式[23]：

(A-1+BC-1D)-1=A-AB(DAB+C)-1DA

(27)

有

(28)

將式(28)代入式(26)，得

(29)

整理得

(30)

再用一次矩陣和求逆公式，有

(31)

將式(31)代入式(30)，有

(32)

(33)

其后驗分布協(xié)方差為

(34)

將式(33)代入式(34)，整理得

(35)

(36)

又因為：

(37)

(38)

將式(37)和式(38)代入式(36)，有

(39)

(40)

式(33)和式(40)就是引入Huber損失函數(shù)之后，本文重新推導(dǎo)的PLF的更新步驟。也就是說，得到了對異常值魯棒的后驗線性化濾波器(Outlier-Robust Posterior Linearization Filter, ORPLF)。

下面給出ORPLF算法的步驟：

算法1ORPLF算法

步驟1預(yù)測步驟：與傳統(tǒng)的廣義高斯濾波[23]算法完全相同：

(41)

(42)

式(41)和式(42)中的積分可采用無跡變換(Unscented Transform, UT)[7]、球面容積積分[8]、高斯-埃爾米特積分[9]等進行近似計算。

步驟2更新步驟：

1)觀測方程線性化：依式(6)～式(12)，將觀測方程hk(xk)線性化。

2)計算尺度函數(shù)：依式(17)，式(20)和式(21)，計算對角矩陣Γk。

3)狀態(tài)更新：

(43)

(44)

2 Huber函數(shù)的相關(guān)參數(shù)估計

2.1 Huber函數(shù)調(diào)諧參數(shù)μ的確定方法

式(15)和式(16)表明：Huber損失函數(shù)是一個分段函數(shù)，調(diào)諧參數(shù)μ為閾值，用于判斷觀測是否屬于異常值。若觀測不是異常值，那么Huber函數(shù)就是l2范數(shù)，最小化式(18)就等價于MAP。若觀測是異常值，那么Huber函數(shù)就是l1范數(shù)。在算法中的直觀表現(xiàn)就是依據(jù)真實值與預(yù)測值之間的歸一化殘差υk，動態(tài)地調(diào)整系統(tǒng)模型中的觀測協(xié)方差：歸一化殘差越大，相應(yīng)的觀測協(xié)方差就越大，反之亦然。

判斷觀測是否屬于異常值，一方面取決于觀測的真實值與預(yù)測值之間的歸一化殘差向量υk，另一方面取決于調(diào)諧參數(shù)μ的取值。據(jù)式(17)知：歸一化殘差向量υk與線性化后修正的觀測方差(Ωk+Rk)有關(guān)。由于線性化誤差Ωk通常很小[12]，因此主要取決于Rk。也就是說，觀測yk是否屬于異常值同時也取決于觀測方差Rk的取值。而在Huber函數(shù)中，調(diào)諧參數(shù)μ是判斷觀測是否屬于異常值的閾值，殘差的絕對值超過此閾值時，就判定為異常值，低于此閾值則判定為非異常值。

若隨機觀測噪聲服從高斯分布，那么，據(jù)3σ原則：在高斯分布中，數(shù)值分布在(m-3σ,m+3σ)中的概率約為99.74%，其中m是分布均值，σ是分布方差。因此，對于高斯分布，若將調(diào)諧參數(shù)μ取為3倍方差，就能很好地判斷υk是否屬于異常值。由于υk是經(jīng)過歸一化的殘差，方差為1，因此可取μ=3。

若隨機觀測噪聲屬于非高斯分布，據(jù)文獻[24]，理論上，Huber函數(shù)的調(diào)諧參數(shù)μ的最優(yōu)值μ*與分布中的異常值概率ε*(0<ε*<1)之間存在如下函數(shù)關(guān)系：

(45)

式中：當(dāng)ε*→0時，μ*→∞，此時Huber函數(shù)就完全是l2范數(shù)；當(dāng)ε*→1時，μ*→0，此時Huber就完全是l1范數(shù)。在Huber函數(shù)中，調(diào)諧參數(shù)μ的取值影響異常值的判定，直觀上來說，調(diào)諧參數(shù)μ越小，就會有越多的樣本被判定為異常值，反之亦然。因此，Huber函數(shù)的調(diào)諧參數(shù)μ應(yīng)當(dāng)與噪聲中異常值的占比相匹配。若調(diào)諧參數(shù)μ取值過大，就會導(dǎo)致一些異常值被判定為非異常值。這樣再據(jù)l2范數(shù)進行優(yōu)化時，會導(dǎo)致較大的狀態(tài)估計誤差，造成算法魯棒性不佳；若調(diào)諧參數(shù)μ取值過小，則會導(dǎo)致一些非異常值被判定為異常值，若對其進行l(wèi)1范數(shù)優(yōu)化，就失去了l2范數(shù)的低誤差擬合性，而造成算法性能下降。

在式(45)中，μ*不可顯式求解。為了近似求解μ*，定義函數(shù)g(μ*)：

(46)

對式(46)求導(dǎo)，得

(47)

因此，可知函數(shù)g(μ*)是單調(diào)減函數(shù)。

在式(46)中，注意到，當(dāng)μ*=3時，ε*≈2.547×10-4?？紤]到函數(shù)的單調(diào)性，若ε*≤2.547×10-4，可近似地認為，理論上調(diào)諧參數(shù)的最優(yōu)值μ*≥3。因此，在實際應(yīng)用中，在近似認為隨機觀測噪聲服從高斯分布的場合，可將調(diào)諧參數(shù)確定為μ=3。而當(dāng)ε*∈(2.547×10-4,1)時，由于μ*→0時g(μ*)→∞，且μ*=3時g(μ*)≈1.003-1/(1-ε*)<1.003-1/(1-2.547×10-4)≈0，因此可以采用二分法近似求解調(diào)諧參數(shù)的最優(yōu)值μ*，并將此最優(yōu)值確定為Huber函數(shù)中的調(diào)諧參數(shù)μ。

在實際應(yīng)用中，若隨機觀測噪聲中的異常值概率未知，那么可以根據(jù)文獻[14]推薦的調(diào)諧參數(shù)典型值，確定Huber函數(shù)的調(diào)諧參數(shù)μ=1.345。

2.2 箱線圖法確定噪聲分布參數(shù)

1.2節(jié)的ORPLF算法，針對的是觀測噪聲分布已知的非線性系統(tǒng)。若觀測噪聲分布未知，需要在此基礎(chǔ)上再加上參數(shù)估計算法。采用箱線圖法[21]檢測異常值，并估計分布參數(shù)，就進一步得到了對異常值和未知觀測噪聲魯棒的后驗線性化濾波器(outlier and unknown Observation Noise-Robust Posterior Linearization Filter, ONRPLF)。

箱線圖是一種顯示一組數(shù)據(jù)分布情況的方法。在統(tǒng)計學(xué)中，將所有樣本的數(shù)值從小到大排列，并分成四等分，其中處于3個分割點位置的數(shù)值就稱為四分位數(shù)，按照數(shù)值從小到大的順序，依次是下四分位數(shù)Q1、中位數(shù)Q2和上四分位數(shù)Q3，上四分位數(shù)與下四分位數(shù)的差值就稱為四分位數(shù)差(Inter Quartile Range, IQR)。可利用箱線圖來檢測異常值[21]，即：大于上四分位數(shù)1.5倍四分位數(shù)差的值，以及小于下四分位數(shù)1.5倍四分位數(shù)差的值，都可劃為異常值。由于四分位數(shù)受異常值影響較小，因此箱線圖法具有較高的魯棒性[21]。

利用箱線圖法檢測出異常值之后，將數(shù)據(jù)中的異常值除去，得到了經(jīng)清洗后的數(shù)據(jù)，再對它們計算其統(tǒng)計特征，可得到濾波需要的統(tǒng)計分布參數(shù)。箱線圖確定噪聲分布參數(shù)的算法，可由以下計算步驟組成：

算法2箱線圖法確定統(tǒng)計分布的參數(shù)

若一組數(shù)據(jù){Di|i=1,2,…,ND}中存在異常值，執(zhí)行以下步驟確定其去除異常值后統(tǒng)計分布的參數(shù)：

步驟1計算上四分位數(shù)Q3和下四分位數(shù)Q1。

步驟2計算四分位數(shù)差Riq=Q3-Q1。

步驟3檢測異常值：若Q1-1.5Riq≤Di≤Q3+ 1.5Riq，則Di不是異常值；否則，Di是異常值。

步驟4計算分布參數(shù)：根據(jù)異常值占總樣本數(shù)的比例確定異常值概率；將樣本去除異常值之后，計算其分布參數(shù)，包括均值和協(xié)方差等。

由于箱線圖法確定分布參數(shù)需要一組具有一定樣本量的噪聲數(shù)據(jù)，因此在確定噪聲分布參數(shù)時，可對觀測序列進行分塊處理，這樣就能得到一組樣本。即對觀測序列塊k=jL-L+1:jL(j=1,2,…,L, 其中L是序列塊的長度)，進行狀態(tài)估計，然后采用箱線圖法就可得到當(dāng)前待求的未知參數(shù)，再進行下一個觀測序列塊的狀態(tài)估計。

本文ONRPLF算法的步驟如下：

算法3ONRPLF算法

當(dāng)j= 0, 1, …，執(zhí)行以下操作：

步驟1對時刻k=jL+ 1：jL+L，執(zhí)行算法1所示的ORPLF算法，得到狀態(tài)xjL+1:xjL+L的分布。

觀測序列塊的長度L，既會影響箱線圖法推斷方差的準確性，也會影響ONRPLF的收斂速度。根據(jù)文獻[25]，對于一組采樣，若要求方差推斷的誤差更低，應(yīng)適當(dāng)增加樣本量。因此，若L過小，則方差推斷的準確性就會降低，進而影響到濾波性能，從而使ONRPLF的誤差增大。然而，由于ONRPLF必須等全部采集完下一個觀測序列塊數(shù)據(jù)后才能批處理進行估計，因此，若L過大，將使ONRPLF的收斂速度下降，且影響算法的實時性。因此，在實際應(yīng)用中，應(yīng)根據(jù)實際情況和需求合理選擇觀測序列塊的長度L。

根據(jù)文獻[25]，方差推斷所需的樣本容量L-與事先給定的置信水平1-α，及估計方差偏離置信區(qū)間中點的相對誤差ξ之間存在如下關(guān)系：

(48)

式中：Z1-α/2表示標(biāo)準正態(tài)分布分位數(shù)值。

若取置信水平1-α=95%，相對誤差ξ=0.2，則由式(48)可得L-≈42.4。考慮到箱線圖法在方差推斷之前已經(jīng)剔除了異常值，因此在確定樣本量時應(yīng)當(dāng)留有一定的余量，即觀測序列塊的長度L應(yīng)當(dāng)大于L-。又因為過大的L將造成收斂速度和實時性的下降，因此在本文的仿真實驗中，將觀測序列塊的長度取為L=50。

此外，由于Huber函數(shù)的調(diào)諧參數(shù)μ是依據(jù)異常值出現(xiàn)的概率εk確定的，而εk的精確度又受限于參數(shù)估計的樣本數(shù)，即觀測序列塊的長度L，因此求μ近似解的精度也不宜定得過高?？紤]到L=50，因此εk的精度是±0.02，在本文中，可將二分法求近似解的精度定為±0.02。

3 仿真實驗

本節(jié)將進行仿真實驗，以驗證本文提出的ORPLF和ONRPLF算法的有效性和魯棒性。針對一種機動目標(biāo)追蹤(Maneuvering Target Tracking, MTT)[12]問題，提出兩種算法與PF、ORKF、PLF、MCUKF和HUF算法的性能進行比較。這是因為PF是非線性濾波問題的經(jīng)典算法[5]；文獻[20]的研究表明：ORKF算法是目前若干種對異常值魯棒濾波算法中性能最好的算法；文獻[13]表明：MCUKF算法是目前對長尾噪聲最魯棒的算法；而HUF則是Huber類型濾波中的經(jīng)典算法[17]。

該MTT模型的狀態(tài)轉(zhuǎn)移方程為

xk+1=F(Δk)xk+wk

(49)

式中：

F(Δk)=

(50)

式中：τ=0.2 s為采樣間隔；wk∈Rn為隨機狀態(tài)轉(zhuǎn)移噪聲，服從零均值協(xié)方差為Qk的高斯分布：

(51)

式中：q1和q2為運動模型的參數(shù)，q1=0.5 m2/s3，q2=10-6rad2/s3。

MTT模型的觀測方程為

(52)

式中：arctan2(·,·)為四象限反正切函數(shù)；h=50 m是目標(biāo)的高度；vk為一種含異常值的隨機觀測噪聲，服從長尾分布，可用如下多維高斯混合分布來描述[12]：

p(vk)=(1-ε)N(vk|0,Σ1)+εN(vk|0,Σ2)

(53)

式中：ε為一個參數(shù)，表示異常值出現(xiàn)的概率；Σ1=diag([1 000,(3π/180)2,(30π/180)2, 1])是背景高斯噪聲的協(xié)方差；Σ2=diag([1 000,(30π/180)2,(30π/180)2, 100])是異常值的協(xié)方差。

在仿真實驗中，本文算法ORPLF和ONRPLF中的積分均采用UT變換進行逼近，UT變換中的參數(shù)與文獻[7]中的相同。ONRPLF算法中L=50。更新Huber函數(shù)的調(diào)諧參數(shù)μ時二分法求近似解的精確度為±0.02。PF中采樣粒子數(shù)為5 000。ORKF、PLF、MCUKF和HUF中的參數(shù)分別采用文獻[20]、文獻[12]、文獻[13]和文獻[17]中性能最佳的參數(shù)。

本文將采用目標(biāo)位置坐標(biāo)的均方根誤差(Root Mean Square Error, RMSE)作為評價指標(biāo)，來量化各種濾波器的性能。為獲取RMSE，每組實驗進行100次蒙特卡洛仿真。

3.1 不含異常值的仿真

首先來看不含異常值的情況，即異常值出現(xiàn)的概率ε= 0。此時，對每一種算法，給定觀測噪聲的協(xié)方差為R0=Σ1=diag([1 000,(3π/180)2,(30π/180)2, 1])。由于UKF是一種廣泛應(yīng)用的非線性濾波器[7]，因此在仿真實驗中又加上了這種濾波算法。由于觀測噪聲服從高斯分布，因此在本文算法ORPLF中，更新調(diào)諧參數(shù)μ時二分法求近似解的精確度為±0.02。Huber函數(shù)的調(diào)諧參數(shù)μ=3。ONRPLF算法中L=50，Huber函數(shù)的調(diào)諧參數(shù)賦初值μ=3。仿真結(jié)果如圖1和表1所示。

在表1中，本文兩種算法ORPLF與ONRPLF是平均RMSE最小的算法。在圖1中，小圖展示了RMSE范圍在7～15 m之間的圖形細節(jié)?？梢园l(fā)現(xiàn)：本文兩種算法ORPLF和ONRPLF，以及算法MCUKF的RMSE最終都收斂到了與UKF相近的性能，其中本文算法收斂最快，MCUKF收斂最慢。此外，在圖1中，ORKF和PLF算法收斂較快，但最終并沒有收斂到與UKF相近的RMSE。PF和HUF的RMSE較大，其中HUF收斂更慢。在不含異常值的條件下，除去本文的兩種算法，這些適用于非高斯環(huán)境的算法性能都不如UKF，這是因為它們都為了保證對異常值的魯棒，因而犧牲了一部分高斯噪聲下的性能。

圖1 不含異常值的仿真

表1 不含異常值的平均RMSE

3.2 已知準確噪聲協(xié)方差的仿真

接下來仿真觀測含有異常值的情況，異常值出現(xiàn)的概率ε=0.4。對于算法RKF、PLF、MCUKF、HUF，及本文算法ORPLF和ONRPLF，仍給定觀測噪聲的協(xié)方差為R0=Σ1=diag([1 000,(3π/180)2,(30π/180)2, 1])。對于算法PF1，給定真實的觀測噪聲分布；對于算法PF2，給定觀測噪聲分布是協(xié)方差為R0=Σ1=diag([1 000,(3π/180)2,(30π/180)2, 1])的高斯分布。由于觀測噪聲服從非高斯分布，因此在ORPLF算法中，Huber函數(shù)的調(diào)諧參數(shù)為μ=1.345。ONRPLF算法中L=50，Huber函數(shù)的調(diào)諧參數(shù)的初始值同樣取μ=1.345，更新調(diào)諧參數(shù)μ時二分法求近似解的精度為±0.02。圖2給出了這些算法在已知準確噪聲協(xié)方差時的RMSE，表2給出了它們的RMSE平均值。

圖2 已知準確噪聲協(xié)方差時的仿真

表2 已知準確噪聲協(xié)方差的算法平均RMSE

從表2中可以發(fā)現(xiàn)：PF2的平均RMSE最大，這是因為粒子濾波利用了錯誤的噪聲先驗信息；獲得了準確先驗的粒子濾波算法PF1的平均RMSE最?。欢疚乃惴∣RPLF和ONRPLF的性能最接近PF1；本文算法ONRPLF的性能略差于ORPLF，這是因為箱線圖法估計的噪聲協(xié)方差參數(shù)不可避免有一定誤差，因此也就制約了濾波器的性能。

在圖2中，小圖展示了RMSE范圍在10～25 m 之間的圖形細節(jié)。從圖2中可以發(fā)現(xiàn)：本文算法ORPLF和ONRPLF收斂較之PF1較慢；MCUKF最終也達到了接近ONRPLF的性能，但是它收斂較之本文算法更慢。在10～40 s的時間區(qū)間內(nèi)，ONRPLF與ORPLF的性能出現(xiàn)了差距，在40 s之后這一差距逐漸消失。這是因為ONRPLF中進行一次參數(shù)估計的觀測序列塊的長度是L=50，由于采樣率為5 Hz，因此每10 s本文ONRPLF算法就對Huber函數(shù)相關(guān)參數(shù)進行一次估計，而10 s則是第一次估計結(jié)果改變發(fā)生的時刻。隨著估計參數(shù)的不斷更新，本文ONRPLF算法就獲得了越來越準確的噪聲協(xié)方差和異常值概率，因此后續(xù)的濾波性能也就更好了。

3.3 噪聲協(xié)方差存在誤差的仿真

在本節(jié)的仿真實驗中，仍取異常值出現(xiàn)的概率為ε=0.4。然而，此時先驗的噪聲協(xié)方差都存在誤差。由于觀測噪聲屬于非高斯分布，因此在ORPLF算法中，Huber函數(shù)的調(diào)諧參數(shù)取值為μ=1.345。在ONRPLF算法中，觀測序列塊的長度L=50，Huber函數(shù)的調(diào)諧參數(shù)的初始值為μ=1.345，更新μ時二分法求近似解的精確度為±0.02。

若給定3組觀測噪聲協(xié)方差，它們分別是真實的觀測噪聲協(xié)方差的λ(λ=2,4,10)倍，則得到如圖3和表3所示的結(jié)果。

圖3 噪聲協(xié)方差存在先驗誤差的仿真1

在表3中，可以發(fā)現(xiàn)：當(dāng)先驗的噪聲協(xié)方差存在誤差時，本文兩種算法ORPLF和ONRPLF具有最小的平均RMSE，表現(xiàn)出了最佳的性能。此外，給定噪聲協(xié)方差與真實值的誤差越大，ON-RPLF算法相較于ORPLF算法的優(yōu)勢就越明顯，這是因為該算法中參數(shù)估計步驟起到了關(guān)鍵的作用。而當(dāng)給定噪聲協(xié)方差僅為真實值2倍時，ORPLF的平均RMSE甚至優(yōu)于ONRPLF，這是由于ONRPLF中以箱線圖法估計噪聲協(xié)方差的誤差，導(dǎo)致了其濾波性能的下降更顯著于ORPLF使用的有誤差的噪聲協(xié)方差。在表3中，PLF、ORKF、MCUKF和HUF算法受觀測噪聲協(xié)方差誤差的影響，都存在不同程度的性能下降；PF算法則極易受先驗信息誤差的影響。

表3 噪聲協(xié)方差存在先驗誤差的算法平均RMSE(仿真1)

在圖3(a)和圖3(b)中，小圖展示了RMSE范圍在10～25 m之間的圖形細節(jié)。從圖3中可以清晰地看到，隨著給定噪聲協(xié)方差與真實值的差距增大，算法收斂之后，ONRPLF相較于ORPLF的優(yōu)勢越來越明顯，而其他算法的性能皆不如本文的兩種算法，這與表3中得到的結(jié)論一致。從時間上來看，在10 s之前，ORPLF與ONRPLF之間的差距不明顯，在10 s之后開始出現(xiàn)差距。在圖3(a)中，ONRPLF的性能在10～20 s區(qū)間內(nèi)弱于ORPLF，在20 s之后逐漸趕上，在30 s 之后又出現(xiàn)了差距；在圖3(b)中，ONRPLF的性能在10～20 s區(qū)間內(nèi)略弱于ORPLF，在20 s 之后性能超越了ORPLF；在圖3(c)中，ONRPLF在10 s之后性能就超越了ORPLF，并且差距越來越大：總之，在10 s之后，隨著時間的推移，觀測數(shù)據(jù)量不斷增加，此時，相較于ORPLF，ONRPLF的性能越來越好。正如3.2節(jié)中的分析，這是因為在本仿真實驗中，每10 s本文ONRPLF算法就對背景高斯噪聲的協(xié)方差及異常值概率進行一次估計，并據(jù)此更新Huber函數(shù)相關(guān)參數(shù)，而10 s正是進行第一次參數(shù)估計的時刻。隨著估計參數(shù)的不斷更新，本文ONRPLF算法就獲得了越來越準確的噪聲協(xié)方差和異常值概率，因此后續(xù)的濾波性能也就更好了。

再給定3組觀測噪聲協(xié)方差，它們分別是真實的觀測噪聲協(xié)方差的0.5、0.2、0.1倍，則得到如圖4和表4所示的結(jié)果。

圖4 噪聲協(xié)方差存在先驗誤差的仿真2

在表4中，本文兩種算法的性能仍是最佳的。與表3不同的是，ONRPLF算法的性能始終稍弱于ORPLF算法，這一方面是因為ONRPLF算法中采用箱線圖法進行參數(shù)估計不可避免存在誤差；另一方面是觀測噪聲協(xié)方差給定為真實值的0.5、0.2、0.1倍時，實際上與真實值的差異并不大，而ORPLF算法的魯棒性確實極佳，在這樣的先驗信息差異下仍能保持較好的性能。

表4 噪聲協(xié)方差存在先驗誤差的算法平均RMSE(仿真2)

在圖4(a)和圖4(b)中，小圖展示了RMSE范圍在10～25 m之間的圖形細節(jié)。在圖4中，相較于ORPLF、ONRPLF在10～20 s區(qū)間內(nèi)都出現(xiàn)了不同程度的性能下降，在20 s之后兩種算法的差距逐漸減小。正如3.2和3.3節(jié)中的分析，這是由于在本仿真實驗中ONRPLF算法更新Huber函數(shù)相關(guān)參數(shù)的時間間隔是10 s，而10 s 正是進行第一次更新的時刻。隨著估計參數(shù)的不斷更新，Huber函數(shù)相關(guān)參數(shù)獲得了更優(yōu)的估計，因此ONRPLF的性能將逐漸提高。在圖4中，PLF和HUF在先驗信息存在差異時都表現(xiàn)出了不同程度的性能下降，但其魯棒性較之ORKF算法更好；PF算法受先驗信息的誤差影響較大；MCUKF算法則極易受觀測噪聲協(xié)方差誤差的影響，當(dāng)先驗信息誤差較大時性能下降十分明顯。

3.4 計算復(fù)雜度分析

為了比較這些算法的計算復(fù)雜度，在3.1節(jié)的仿真實驗中，記錄了每一種算法運行所需的平均運行時間，其結(jié)果如表5所示。仿真實驗的平臺為64位Win10操作系統(tǒng)，內(nèi)存8 GB，Intel處理器，內(nèi)核i7-4790，CPU3.6 G，IDLE為Python3.8。

在表5中，PF的平均運行時間遠遠超過其他算法，這是該算法為了達到低誤差付出的運算復(fù)雜度代價。此外，與其他的濾波算法不同，PF適用于任何分布的噪聲，而非像其他算法一樣僅適用于高斯噪聲與異常值混合這種特定的分布。本文兩種算法ORPLF和ONRPLF的平均運行時間與MCUKF相近，而遠低于ORKF、PLF和HUF。在本文兩種算法中，ONRPLF的平均運行時間又略高于ORPLF，這是因為相較于ORPLF算法，ONRPLF算法多了參數(shù)估計的步驟。

表5 算法運行時間比較

4 結(jié) 論

針對含異常觀測值的非線性系統(tǒng)濾波問題，本文提出了兩種魯棒的濾波算法ORPLF和ONRPLF。

1)分析表明：由于Huber函數(shù)兼顧了l1范數(shù)的魯棒性，因此用Huber損失函數(shù)代替推導(dǎo)PLF的MAP準則中觀測誤差的l2范數(shù)，構(gòu)造出了一種新的準則函數(shù)，并由此推導(dǎo)出的ORPLF算法對異常值具有魯棒性。

2)ORPLF僅適用于隨機噪聲的分布參數(shù)已知的情況，濾波器的參數(shù)由先驗給定；當(dāng)隨機噪聲的分布參數(shù)未知時，可利用箱線圖法檢測異常值，并進行Huber函數(shù)的相關(guān)參數(shù)估計，動態(tài)地調(diào)整濾波器的參數(shù)，這樣就得到了ONRPLF算法。

3)仿真驗證了分析結(jié)果的有效性，同時也表明：在觀測噪聲統(tǒng)計分布未知且含異常觀測值的環(huán)境中，本文算法性能優(yōu)于現(xiàn)有文獻報道的非線性濾波類算法。