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

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

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

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

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

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

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

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

1.1 后驗(yàn)線性化濾波器

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

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

(1)

yk=hk(xk)+vk

(2)

(3)

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

文獻(xiàn)[12]給出了一種通過(guò)最小化KLD來(lái)進(jìn)行非線性觀測(cè)方程線性化的方法，其KLD定義為

(4)

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

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

(5)

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

(6)

(7)

(8)

(9)

(10)

(11)

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

yk=Akxk+bk+rk

(12)

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

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

1.2 異常值魯棒的后驗(yàn)線性化濾波器

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

(13)

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

(14)

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

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

(15)

(16)

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

(17)

式中：A1/2為對(duì)稱正定矩陣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)

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

(19)

(20)

定義矩陣：

(21)

那么由式(17)，有

(22)

將式(20)～式(22)代入式(19)，并利用文獻(xiàn)[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)

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

(34)

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

(35)

(36)

又因?yàn)椋?/p>

(37)

(38)

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

(39)

(40)

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

下面給出ORPLF算法的步驟：

算法1ORPLF算法

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

(41)

(42)

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

步驟2更新步驟：

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

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

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

(43)

(44)

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

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

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

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

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

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

(45)

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

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

(46)

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

(47)

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

在式(46)中，注意到，當(dāng)μ*=3時(shí)，ε*≈2.547×10-4。考慮到函數(shù)的單調(diào)性，若ε*≤2.547×10-4，可近似地認(rèn)為，理論上調(diào)諧參數(shù)的最優(yōu)值μ*≥3。因此，在實(shí)際應(yīng)用中，在近似認(rèn)為隨機(jī)觀測(cè)噪聲服從高斯分布的場(chǎng)合，可將調(diào)諧參數(shù)確定為μ=3。而當(dāng)ε*∈(2.547×10-4,1)時(shí)，由于μ*→0時(shí)g(μ*)→∞，且μ*=3時(shí)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ù)μ。

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

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

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

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

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

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

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

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

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

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

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

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

本文ONRPLF算法的步驟如下：

算法3ONRPLF算法

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

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

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

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

(48)

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

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

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

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

本節(jié)將進(jìn)行仿真實(shí)驗(yàn)，以驗(yàn)證本文提出的ORPLF和ONRPLF算法的有效性和魯棒性。針對(duì)一種機(jī)動(dòng)目標(biāo)追蹤(Maneuvering Target Tracking, MTT)[12]問(wèn)題，提出兩種算法與PF、ORKF、PLF、MCUKF和HUF算法的性能進(jìn)行比較。這是因?yàn)镻F是非線性濾波問(wèn)題的經(jīng)典算法[5]；文獻(xiàn)[20]的研究表明：ORKF算法是目前若干種對(duì)異常值魯棒濾波算法中性能最好的算法；文獻(xiàn)[13]表明：MCUKF算法是目前對(duì)長(zhǎng)尾噪聲最魯棒的算法；而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為隨機(jī)狀態(tài)轉(zhuǎn)移噪聲，服從零均值協(xié)方差為Qk的高斯分布：

(51)

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

MTT模型的觀測(cè)方程為

(52)

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

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

(53)

式中：ε為一個(gè)參數(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é)方差。

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

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

3.1 不含異常值的仿真

首先來(lái)看不含異常值的情況，即異常值出現(xiàn)的概率ε= 0。此時(shí)，對(duì)每一種算法，給定觀測(cè)噪聲的協(xié)方差為R0=Σ1=diag([1 000,(3π/180)2,(30π/180)2, 1])。由于UKF是一種廣泛應(yīng)用的非線性濾波器[7]，因此在仿真實(shí)驗(yàn)中又加上了這種濾波算法。由于觀測(cè)噪聲服從高斯分布，因此在本文算法ORPLF中，更新調(diào)諧參數(shù)μ時(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之間的圖形細(xì)節(jié)。可以發(fā)現(xiàn)：本文兩種算法ORPLF和ONRPLF，以及算法MCUKF的RMSE最終都收斂到了與UKF相近的性能，其中本文算法收斂最快，MCUKF收斂最慢。此外，在圖1中，ORKF和PLF算法收斂較快，但最終并沒(méi)有收斂到與UKF相近的RMSE。PF和HUF的RMSE較大，其中HUF收斂更慢。在不含異常值的條件下，除去本文的兩種算法，這些適用于非高斯環(huán)境的算法性能都不如UKF，這是因?yàn)樗鼈兌紴榱吮ＷC對(duì)異常值的魯棒，因而犧牲了一部分高斯噪聲下的性能。

圖1 不含異常值的仿真

表1 不含異常值的平均RMSE

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

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

圖2 已知準(zhǔn)確噪聲協(xié)方差時(shí)的仿真

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

從表2中可以發(fā)現(xiàn)：PF2的平均RMSE最大，這是因?yàn)榱Ｗ訛V波利用了錯(cuò)誤的噪聲先驗(yàn)信息；獲得了準(zhǔn)確先驗(yàn)的粒子濾波算法PF1的平均RMSE最小；而本文算法ORPLF和ONRPLF的性能最接近PF1；本文算法ONRPLF的性能略差于ORPLF，這是因?yàn)橄渚€圖法估計(jì)的噪聲協(xié)方差參數(shù)不可避免有一定誤差，因此也就制約了濾波器的性能。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

表5 算法運(yùn)行時(shí)間比較

4 結(jié) 論

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

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

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

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