方安然,李 旦,張建秋

(1.復(fù)旦大學(xué)電磁波信息科學(xué)教育部重點(diǎn)實(shí)驗(yàn)室,上海 200433；2.復(fù)旦大學(xué)電子工程系,上海 200433)

0 引 言

在雷達(dá)、聲吶、通信、語(yǔ)音識(shí)別等應(yīng)用中,觀測(cè)中隨機(jī)出現(xiàn)異常值是一種十分常見(jiàn)的現(xiàn)象[1-2]。例如,無(wú)線通信中,電路通斷暫態(tài)過(guò)程產(chǎn)生的脈沖干擾;在雷達(dá)或聲吶系統(tǒng)中,人為或自然因素產(chǎn)生的沖擊性干擾,而引起雷達(dá)或聲納觀測(cè)的隨機(jī)異常波動(dòng)等[3]。文獻(xiàn)[4]表明:這種觀測(cè)噪聲是非高斯的,其概率密度分布有一個(gè)很厚重的“尾部”,因此通常稱(chēng)其為長(zhǎng)尾分布/噪聲?？柭鼮V波(Kalman filter,KF)[5],由于其成立的條件是觀測(cè)噪聲為高斯白噪聲,因此當(dāng)觀測(cè)存在隨機(jī)異常值時(shí),其性能會(huì)顯著下降,甚至失效[2,6]。

針對(duì)具有隨機(jī)觀測(cè)異常值的噪聲環(huán)境,人們研究了不同的模型方法和濾波算法。就模型方法而言,文獻(xiàn)[7]發(fā)現(xiàn):高斯分布與均勻分布的疊加混合,可以反映含異常值觀測(cè)噪聲的“長(zhǎng)尾”特性。然而,這種建模方法只適合于觀測(cè)數(shù)據(jù)有界的系統(tǒng)。針對(duì)觀測(cè)數(shù)據(jù)無(wú)界的系統(tǒng),文獻(xiàn)[3,8]則分別報(bào)道了利用高斯與學(xué)生t分布疊加混合,或與拉普拉斯分布疊加混合的模型方法。可是,它們的數(shù)學(xué)描述都比較復(fù)雜,不便于參數(shù)估計(jì)。為了解決上述問(wèn)題,文獻(xiàn)[2]報(bào)道了利用學(xué)生t分布來(lái)對(duì)含異常值的噪聲進(jìn)行建模的方法。但由于其完全依賴(lài)于學(xué)生t分布,因此有可能無(wú)法適應(yīng)觀測(cè)噪聲可能存在的各種不同分布。

針對(duì)非高斯系統(tǒng)的濾波,主要算法有高斯和濾波(Gaussian sum filtering,GSF)[8]和粒子濾波(particle filtering,PF)[7]等。盡管這兩種算法能處理非高斯噪聲環(huán)境中的線性系統(tǒng)濾波問(wèn)題,但都存在運(yùn)算復(fù)雜度偏高的問(wèn)題[7]。為了解決這一問(wèn)題,文獻(xiàn)[9]基于新息正交性,提出了一種修正KF (modified Kalman filtering,MKF)算法。盡管這個(gè)濾波算法對(duì)噪聲先驗(yàn)信息是否準(zhǔn)確并不敏感,但其收斂過(guò)程較慢。針對(duì)這一問(wèn)題,文獻(xiàn)[10-11]提出了一類(lèi)基于相關(guān)熵的濾波方法,包括相關(guān)熵濾波(correntropy filtering,CF)[10]、修正CF(modified CF,MCF)[11]和最大相關(guān)熵KF(maximum correntropy criterion KF,MCC-KF)[11]等方法,加快了算法的收斂速度。不過(guò),其收斂性高度依賴(lài)濾波參數(shù)選擇是否合適。如何選擇濾波參數(shù),目前還沒(méi)有一般性的辦法。為此,文獻(xiàn)[12]通過(guò)引入假設(shè)檢驗(yàn),提出了一種自適應(yīng)魯棒卡爾曼濾波(adaptive robust KF,ARKF)方法?？捎捎诩僭O(shè)檢驗(yàn)需要噪聲協(xié)方差的準(zhǔn)確先驗(yàn),因此應(yīng)用受限。為了進(jìn)一步提高濾波算法的魯棒性,文獻(xiàn)[13]在最小化最大魯棒估計(jì)的準(zhǔn)則下,提出了一種中心化融合的魯棒KF (robust centralized fusion KF,RCFKF)方法。該方法僅適用于未知參數(shù)值有界的情況,且估計(jì)誤差較大。

針對(duì)存在觀測(cè)異常值的線性系統(tǒng),本文提出了一種新的濾波算法,稱(chēng)之為對(duì)異常值魯棒KF (outlier-robust KF,ORKF)算法。分析表明:當(dāng)最大后驗(yàn)(maximum a posterior,MAP)準(zhǔn)則中加權(quán)觀測(cè)誤差的l2范數(shù),用Huber損失函數(shù)代替后,就構(gòu)造了一個(gè)新的最優(yōu)化準(zhǔn)則函數(shù)。由于Huber損失函數(shù)可同時(shí)描述l1和l2范數(shù)[14],因此借助于這個(gè)新的最優(yōu)化準(zhǔn)則函數(shù),本文推導(dǎo)的卡爾曼濾波器就具有了l1范數(shù)對(duì)異常值的魯棒性。而當(dāng)含異常觀測(cè)值的未知統(tǒng)計(jì)分布,利用具有未知參數(shù)的高斯混合模型(Gaussian mixture model,GMM)描述時(shí),那么就可根據(jù)變分貝葉斯(variational Bayes,VB)的逼近思想,利用貪婪期望最大化(expectation-maximization,EM)算法[15]來(lái)對(duì)GMM模型中的未知參數(shù)進(jìn)行估計(jì)。如此,本文借助于這一思想,進(jìn)一步提出了對(duì)異常值和未知觀測(cè)噪聲同時(shí)魯棒的KF算法,并稱(chēng)其為帶VB的ORKF (ORKF-VB)算法。在仿真和實(shí)測(cè)實(shí)驗(yàn)驗(yàn)證分析結(jié)果的同時(shí),也展示了提出算法的濾波性能在異常值和觀測(cè)噪聲統(tǒng)計(jì)分布未知時(shí),均優(yōu)于現(xiàn)有文獻(xiàn)報(bào)道魯棒卡爾曼濾波類(lèi)的算法。

1 異常值的統(tǒng)計(jì)模型

異常值通常用具有長(zhǎng)尾分布特性的非高斯噪聲來(lái)進(jìn)行描述[4],這是因?yàn)楹惓Ｖ档碾S機(jī)噪聲的分位數(shù)-分位數(shù)圖(quantile-quantile plot,QQ-Plot)在原點(diǎn)附近近似線性,即在均值附近,其分布近似于高斯分布;而在尾部區(qū)域,QQ-Plot偏離了線性,那么長(zhǎng)的尾部區(qū)域就是非高斯分布[16]。位于尾部區(qū)域的數(shù)據(jù)也就稱(chēng)為異常值,它可以描述為一個(gè)具有很大方差的零均值高斯分布[16]。這樣,包含異常值的觀測(cè)噪聲就可由兩個(gè)零均值高斯分布構(gòu)成的GMM來(lái)表示,其中一個(gè)是背景高斯噪聲,另一個(gè)就是異常值的統(tǒng)計(jì)分布[16]。在一個(gè)動(dòng)態(tài)系統(tǒng)中,若vk表示k時(shí)刻的隨機(jī)觀測(cè)噪聲,而p(vk)表示隨機(jī)觀測(cè)噪聲的概率密度分布,那么就有[16]:

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

(1)

式中,權(quán)值ε是異常值出現(xiàn)的概率,表示了長(zhǎng)尾分布的非高斯性,ε越大,該分布的尾部就越“厚”,偏離高斯分布的程度也就越高,反之亦然;N(vk;0,Σ1)表示背景高斯噪聲,Σ1是背景高斯分布的協(xié)方差;N(vk;0,Σ2)則描述了異常值的高斯分布,Σ2是異常值高斯分布的協(xié)方差。

2 KF的回顧

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

xk=Fk-1xk-1+wk

(2)

yk=Hkxk+vk

(3)

當(dāng)式(3)中觀測(cè)不存在異常值時(shí),即觀測(cè)噪聲為高斯白噪聲時(shí),KF的預(yù)測(cè)[7]為

(4)

(5)

從最大后驗(yàn)的角度,KF算法的狀態(tài)更新就是最大化了如下函數(shù)[15]:

(6)

(7)

如此,KF的更新步驟就如下[8]:

(8)

(9)

(10)

3 對(duì)異常值魯棒的KF

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

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

(11)

式中,評(píng)價(jià)函數(shù)ρ(·)的定義為

(12)

(13)

式中,A1/2是對(duì)稱(chēng)正定矩陣A的Cholesky分解;AT/2是A1/2的轉(zhuǎn)置,滿(mǎn)足A=A1/2AT/2;A-1/2是A1/2的逆矩陣;A-T/2是AT/2的逆矩陣。

當(dāng)將式(7)中第一項(xiàng)用式(11)的Huber損失函數(shù)替代時(shí),本文給出了一個(gè)新的損失/評(píng)價(jià)函數(shù)如下:

(14)

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

(15)

(16)

定義矩陣:

(17)

那么據(jù)式(13),就有

(18)

將式(16)～式(18)代入式(15),并利用文獻(xiàn)[14,17]中的矩陣恒等式,就有

(19)

這樣再將式(13)代入式(19),就可得

(20)

(21)

再整理得

(22)

利用矩陣和求逆公式[18]:

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

(23)

那么就有

(24)

將式(24)代入式(22),得

(25)

整理得

(26)

再次利用式(23)的矩陣求逆公式,就有

(27)

將式(27)代入式(26),得

(28)

(29)

(30)

(31)

(32)

(33)

式(29)和式(33)就是引入Huber損失函數(shù)后,狀態(tài)及其協(xié)方差的迭代更新公式,也就是本文提出的ORKF算法,其計(jì)算步驟如下。

算法 1 ORKF算法

步驟 1預(yù)測(cè)

(34)

(35)

步驟 2計(jì)算尺度函數(shù)

(36)

式中,

(37)

(38)

步驟 3更新

(39)

(40)

(41)

(42)

4 參數(shù)估計(jì)

4.1 利用3σ法則確定Huber函數(shù)調(diào)諧參數(shù)μ

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

判斷觀測(cè)是否屬于異常值,一方面取決于觀測(cè)的真實(shí)值與預(yù)測(cè)值之間的歸一化殘差向量ek,另一方面取決于調(diào)諧參數(shù)μ的取值。據(jù)式(13)知:歸一化殘差向量ek與觀測(cè)方差有關(guān)。也就是說(shuō),觀測(cè)yk是否屬于異常值同時(shí)也取決于觀測(cè)方差Rk的取值,而Rk則是背景高斯噪聲的協(xié)方差。

在Huber函數(shù)中,調(diào)諧參數(shù)μ是判斷觀測(cè)是否屬于異常值的閾值,殘差的絕對(duì)值超過(guò)此閾值的判定為異常值,低于此閾值的判定為非異常值。據(jù)3σ原則:在正態(tài)分布中,數(shù)值分布在(m-3σ,m+3σ)中的概率約為99.74%,其中m是分布的均值,σ是分布的標(biāo)準(zhǔn)差。因此,對(duì)于正態(tài)分布,若將調(diào)諧參數(shù)μ取為3倍標(biāo)準(zhǔn)差,就能很好地判斷ek是否來(lái)自異常值。由于ek是經(jīng)過(guò)歸一化的殘差,因此就可取μ=3。

4.2 觀測(cè)噪聲分布未知時(shí)的參數(shù)估計(jì)

本文前面提出的ORKF算法,是針對(duì)隨機(jī)觀測(cè)噪聲分布已知的線性系統(tǒng)。若隨機(jī)觀測(cè)噪聲的分布未知,那么就需要引入VB的思想,同時(shí)對(duì)系統(tǒng)的狀態(tài)和噪聲參數(shù)做出估計(jì),這就是本文提出的ORKF-VB算法。

據(jù)式(1)知隨機(jī)觀測(cè)噪聲的GMM模型為

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

(43)

而未知統(tǒng)計(jì)分布的隨機(jī)觀測(cè)噪聲,就指式(43)中異常值出現(xiàn)的概率ε,以及背景高斯噪聲和異常值的協(xié)方差Σ1和Σ2都未知。

為了估計(jì)式(12)和式(43)中的未知參數(shù)μ、ε、Σ1和Σ2,本文引入了一種貪婪EM算法[15]。這種算法的優(yōu)點(diǎn)主要是在提供了一個(gè)簡(jiǎn)便易行參數(shù)估計(jì)方法的同時(shí),不容易陷入局部極小值,且估計(jì)誤差比較小。

據(jù)文獻(xiàn)[17-19]的在線EM算法,對(duì)觀測(cè)序列進(jìn)行分塊處理,以提高參數(shù)估計(jì)的穩(wěn)健性。即對(duì)觀測(cè)序列塊k=jL-L+1(j=1,2,…,L,其中L是序列塊的長(zhǎng)度),進(jìn)行狀態(tài)估計(jì),并執(zhí)行貪婪EM算法就可得到當(dāng)前待求的未知參數(shù)。對(duì)下一個(gè)觀測(cè)序列塊,以上一個(gè)觀測(cè)序列塊中估計(jì)的未知參數(shù)值為初始值再次開(kāi)始執(zhí)行EM算法,這樣就避免了更新時(shí)可能發(fā)生較大的異變,也就保證了參數(shù)估計(jì)的穩(wěn)健性。

在系統(tǒng)模型式(2)和式(3)中,隨機(jī)觀測(cè)噪聲的分布未知,也就意味著其參數(shù)Rk未知,而Rk則是隨機(jī)觀測(cè)噪聲中排除異常值后的協(xié)方差。由于含異常值的隨機(jī)觀測(cè)噪聲是一個(gè)長(zhǎng)尾分布,其在均值附近主要表現(xiàn)為近似于高斯分布的背景高斯噪聲,而在遠(yuǎn)離均值的部分則表現(xiàn)為異常值,因此Rk就是隨機(jī)觀測(cè)噪聲分布的中心區(qū)域背景高斯噪聲的協(xié)方差。背景高斯噪聲的方差小于異常值的方差,而方差實(shí)際上就是協(xié)方差的跡。因此,在參數(shù)估計(jì)中,Rk就應(yīng)更新為跡較小的協(xié)方差。如此,本文ORKF-VB算法的流程就如下。

算法2 ORKF-VB算法

當(dāng)j=1,2,…,執(zhí)行以下操作:

步驟 1對(duì)時(shí)刻k=jL-L+1:jL,執(zhí)行算法1所示的ORKF算法,得到狀態(tài)xjL-L+1:xjL的分布,其中L是時(shí)間塊的長(zhǎng)度;

步驟 3更新Rk

(44)

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

本節(jié)將進(jìn)行仿真實(shí)驗(yàn),以驗(yàn)證本文提出的ORKF和ORKF-VB算法的有效性。針對(duì)一種基準(zhǔn)導(dǎo)航系統(tǒng)的模型[11],本文將比較提出的兩種算法與MKF[9]、MCC-KF[11]、ARKF[12]和RCFKF[13]算法的性能。其中,MKF和ARKF算法都是針對(duì)異常值問(wèn)題的經(jīng)典算法;MCC-KF是相關(guān)熵類(lèi)濾波方法中最新的算法;RCFKF則是針對(duì)噪聲方差未知的通常算法。

基準(zhǔn)導(dǎo)航系統(tǒng)的動(dòng)態(tài)模型可描述[11]為

(45)

(46)

式中,狀態(tài)變量為

(47)

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

(48)

式中,ε表示異常值出現(xiàn)的概率；Σ1=diag([2,2])是背景高斯噪聲的協(xié)方差;Σ2=diag([200,200])是異常值的協(xié)方差。

本文將采用目標(biāo)位置坐標(biāo)的均方根誤差(root mean square error,RMSE)[11]來(lái)作為算法性能的評(píng)價(jià)指標(biāo),以衡量濾波結(jié)果的準(zhǔn)確性。

5.1 不含異常值的仿真

首先來(lái)看不含異常值的情況,即異常值出現(xiàn)的概率ε=0。此時(shí),對(duì)每一種算法,給定觀測(cè)噪聲的協(xié)方差為R=diag([2,2])。仿真結(jié)果如圖1所示。

圖1 不含異常值的仿真結(jié)果Fig.1 Simulation result without outlier

在圖1中可以發(fā)現(xiàn),RCFKF算法的RMSE最大,其余6種算法的RMSE都能收斂到2.5 m以下,其中ARKF的RMSE稍大,而MKF算法收斂較慢,而此時(shí)本文算法ORKF和ORKF-VB的性能與傳統(tǒng)KF的性能相當(dāng)。

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

接下來(lái)仿真觀測(cè)噪聲含有異常值的情況,異常值出現(xiàn)的概率ε=0.5。

根據(jù)文獻(xiàn)[9]和文獻(xiàn)[11],MKF和MCC-KF兩種算法都是把異常值和背景高斯噪聲視為同一個(gè)統(tǒng)計(jì)分布,因此要求的噪聲協(xié)方差是異常值與背景高斯噪聲混合分布的協(xié)方差,據(jù)式(1)知,這種非高斯噪聲的二階統(tǒng)計(jì)量Σv為

Σv=(1-ε)Σ1+εΣ2

(49)

將ε=0.5,Σ1=diag([2,2],Σ2=diag([200,200])代入式(49),得Σv=diag([101,101])。也就是說(shuō),在這個(gè)仿真實(shí)驗(yàn)中,只有將MKF和MCC-KF算法的噪聲協(xié)方差給定為R=Σv=diag([101,101])時(shí),才能達(dá)到最佳的性能。因此,該方差也就是MKF和MCC-KF算法的準(zhǔn)確方差。

根據(jù)文獻(xiàn)[12],ARKF算法要求的噪聲協(xié)方差是除去異常值后的背景高斯噪聲的協(xié)方差,即R=diag([2,2])。

這樣,對(duì)每一種算法給定正確的噪聲協(xié)方差即為:對(duì)于MKF和MCC-KF算法是R=diag([101,101]),而對(duì)于ARKF和本文算法ORKF與ORKF-VB則是R=diag([2,2])。

圖2中給出了這些算法在已知準(zhǔn)確噪聲協(xié)方差時(shí)的RMSE,可以發(fā)現(xiàn):RCFKF仍是性能最差的一種算法,MKF、MCC-KF和本文算法ORKF與ORKF-VB的RMSE性能則基本相當(dāng)。不過(guò)MKF算法收斂較慢,而ARKF的RMSE算法的性能比上述四者略差。

圖2 已知準(zhǔn)確噪聲協(xié)方差時(shí)的仿真結(jié)果Fig.2 Simulation result with correct noise covariance

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

在本節(jié)的仿真實(shí)驗(yàn)中,仍取異常值出現(xiàn)的概率為ε=0.5。然而,這里先驗(yàn)的噪聲協(xié)方差都存在誤差。若給定3組噪聲協(xié)方差,它們分別是背景高斯噪聲協(xié)方差的2倍、4倍、10倍,則得到如圖3所示的結(jié)果。

圖3 噪聲協(xié)方差存在先驗(yàn)誤差的仿真1Fig.3 Simulation 1 of a priori error in noise covariance

在圖3中,可以發(fā)現(xiàn)本文提出的兩種算法仍然保持了與圖1相近的性能,這表明它們具有很強(qiáng)的魯棒性。不難發(fā)現(xiàn),MKF仍是收斂最慢的算法,RCFKF仍是性能最差的算法。當(dāng)給定噪聲方差的誤差增大到10倍時(shí),ARKF的性能顯著變差,表明該算法的魯棒性差。

當(dāng)給定的噪聲協(xié)方差分別是背景高斯噪聲與異常值的混合噪聲協(xié)方差的2倍、4倍、10倍時(shí),全部算法的性能如圖4所示。由于此時(shí)ARKF算法發(fā)散,因此就沒(méi)有在圖4中繪制出其RMSE曲線。

圖4 噪聲協(xié)方差存在先驗(yàn)誤差的仿真2Fig.4 Simulation 2 of a priori error in noise covariance

在圖4(a)中可以發(fā)現(xiàn),MKF、MCC-KF和本文提出的ORKF-VB算法,收斂后性能幾乎相當(dāng),但MKF收斂較慢。本文提出的ORKF算法性能稍遜于以上3種算法,這是因?yàn)閷?duì)于ORKF來(lái)說(shuō),正確的噪聲協(xié)方差是背景高斯噪聲的協(xié)方差;而對(duì)于MKF和MCC-KF來(lái)說(shuō),正確的噪聲協(xié)方差是背景高斯噪聲與異常值混合噪聲的協(xié)方差。因此，對(duì)本文提出的ORKF算法而言,始終在利用存在誤差的協(xié)方差進(jìn)行估計(jì),且該誤差比給出倍數(shù)還要大。如圖4(a)中給每一種算法的噪聲協(xié)方差都是混合噪聲協(xié)方差的2倍,即R=diag([202,202]),它與背景高斯噪聲的協(xié)方差的差距達(dá)101倍,而與混合噪聲的協(xié)方差的差距僅有2倍,更大的先驗(yàn)誤差造成了算法性能更大地下降。而本文提出的ORKF-VB算法盡管最初也使用了錯(cuò)誤的協(xié)方差,但其在濾波的同時(shí)可對(duì)噪聲協(xié)方差進(jìn)行估計(jì),因此當(dāng)估計(jì)的噪聲協(xié)方差接近正確后,就能表現(xiàn)出更優(yōu)的性能。在圖4(a)中可以看到,ORKF-VB的RMSE在200 s之前與ORKF極為相近,在200 s之后則明顯下降。這是因?yàn)镺RKF-VB中進(jìn)行一次高斯噪聲協(xié)方差估計(jì)的觀測(cè)序列塊的長(zhǎng)度是L=200,由于觀測(cè)的時(shí)間間隔為Δt=1 s,因此每200 s ORKF-VB算法就對(duì)高斯噪聲協(xié)方差進(jìn)行一次估計(jì),而200 s則是第一次估計(jì)結(jié)果改變發(fā)生的時(shí)刻。因此隨著估計(jì)參數(shù)的不斷更新,ORKF-VB算法就獲得了越來(lái)越準(zhǔn)確的噪聲協(xié)方差,后續(xù)的濾波性能也就更好,這一點(diǎn)同樣可在圖4(b)和圖4(c)中發(fā)現(xiàn)。

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

為了比較這些算法的計(jì)算復(fù)雜度,對(duì)第5.1節(jié)中的仿真實(shí)驗(yàn),記錄了每一算法運(yùn)行所需的時(shí)間,如表1所示。仿真實(shí)驗(yàn)的平臺(tái)為64位win10操作系統(tǒng),內(nèi)存8 GB,Intel處理器,內(nèi)核i7-4790,CPU 3.6 G,IDLE為Python 3.8。

表1 算法運(yùn)行時(shí)間比較Table 1 Algorithm run time comparison ms

在表1中,本文算法ORKF的平均運(yùn)行時(shí)間高于KF、RCFKF和MKF,而低于ARKF和MCC-KF。本文算法ORKF-VB的平均運(yùn)行時(shí)間遠(yuǎn)遠(yuǎn)高于其他所有算法,這是因?yàn)槎嗔藚f(xié)方差估計(jì)的步驟,這表明該算法的高魯棒性能是以更高的計(jì)算復(fù)雜度為代價(jià)的。

6 實(shí)際驗(yàn)證

本節(jié)將以鋰電池荷電狀態(tài)監(jiān)測(cè)問(wèn)題為例,進(jìn)行實(shí)際驗(yàn)證,以證明本文算法的有效性。針對(duì)該問(wèn)題,本文利用文獻(xiàn)[21]提出的一種描述鋰電池具有時(shí)變內(nèi)阻和開(kāi)路電壓滯回的系統(tǒng)模型,對(duì)鋰電池荷電狀態(tài)(state of charge,SOC)進(jìn)行估計(jì)。與第5節(jié)相同,本節(jié)將比較提出的兩種算法與KF[5]、MKF[9]、MCC-KF[11]、ARKF[12]和RCFKF[13]算法的性能。

文獻(xiàn)[21]給出的電池模型如圖5所示。在圖5中,電池的容量Cn、電池內(nèi)阻R0、R1和C1、R2和C2構(gòu)成的2個(gè)RC網(wǎng)絡(luò)都是電池的內(nèi)部參數(shù);1V·Z表示電池電量,其中Z表示鋰電池的SOC;I表示流過(guò)電池的電流,充電時(shí)為正,放電時(shí)為負(fù);VOC(Z)表示電池的開(kāi)路電壓。

圖5 鋰電池等效電路模型Fig.5 Equivalent circuit model of lithium-ion battery

文獻(xiàn)[21]描述鋰電池狀態(tài)的系統(tǒng)模型為

(50)

Vk=VOCA,k+VH,k+R0,kIk+V1,k+V2,k+vk

(51)

式中,Zk是k時(shí)刻電池的SOC;V1,k和V2,k反映了k時(shí)刻電池所同時(shí)具有的短和長(zhǎng)時(shí)間常數(shù)的動(dòng)態(tài)特性;R0,k表示k時(shí)刻電池的內(nèi)阻;Ik表示k時(shí)刻流過(guò)電池的電流;狀態(tài)變量xk=[Zk,V1,k,V2,k,R0,k]T;Vk表示電池兩端的電壓;VOC(Z)=VOCA,k+VH,k是k時(shí)刻電池的開(kāi)路電壓,由Zk及電池的充放電狀態(tài)決定,文獻(xiàn)[21]給出了其與Zk之間的函數(shù)關(guān)系;wk是狀態(tài)轉(zhuǎn)移噪聲,vk是隨機(jī)觀測(cè)噪聲。在電子電路中,存在熱噪聲、散彈噪聲、粉紅噪聲等多種隨機(jī)噪聲[22],而電路通斷產(chǎn)生的脈沖及外部干擾等又易引起觀測(cè)的異常波動(dòng),因此隨機(jī)觀測(cè)噪聲中極有可能含有異常值且其分布未知。

實(shí)驗(yàn)中使用了一種鋰離子聚合物電芯LGABF1L18650電池,額定容量3 350 mA·h,額定電壓3.7 V。電壓Vk和電流Ik的數(shù)據(jù)由BQ40Z50芯片實(shí)際采集,采樣周期為1 s。所有實(shí)驗(yàn)均在25℃下進(jìn)行。實(shí)驗(yàn)中電池模型參數(shù)R1=0.001 Ω,C1=618 F;R2=0.025 7 Ω,C2=707.7 F。模型參數(shù)的確定方法及電容Cn的更新方法由文獻(xiàn)[21]給出。電池SOC的真實(shí)值是通過(guò)改進(jìn)的安時(shí)積分法計(jì)算獲得的[23]。實(shí)驗(yàn)中,MKF、MCC-KF、ARKF的參數(shù)設(shè)定分別與文獻(xiàn)[9]、文獻(xiàn)[11]、文獻(xiàn)[12]中狀態(tài)估計(jì)誤差最小的方法相同;對(duì)于所有的算法,狀態(tài)轉(zhuǎn)移噪聲和觀測(cè)噪聲的相關(guān)參數(shù)均由文獻(xiàn)[21]給定。

實(shí)驗(yàn)結(jié)果如圖6和表3所示。圖6(a)顯示了采用SOC估計(jì)的絕對(duì)誤差。圖6(b)顯示了不同SOC區(qū)間內(nèi)的平均相對(duì)誤差,其中SOC區(qū)間的劃分方法由文獻(xiàn)[24]給出,如表2所示。表3則給出了不同方法估計(jì)SOC的均方根誤差RMSE。

圖6 鋰電池SOC估計(jì)的實(shí)驗(yàn)結(jié)果Fig.6 Experimental result of SOC estimation of lithium-ion battery

表2 SOC區(qū)間劃分Table 2 Partition of SOC blocks

表3 鋰電池SOC估計(jì)的RMSETable 3 RMSE of SOC estimation of lithium-ion battery

圖6(a)顯示了SOC估計(jì)的絕對(duì)誤差,圖6(b)顯示了不同SOC區(qū)間內(nèi)的平均相對(duì)誤差。從圖6(b)中可以看到,ARKF性能受SOC區(qū)間的影響最明顯,當(dāng)SOC在0.9～1范圍內(nèi)時(shí)性能最好。然而,隨著SOC逐漸減小,ARKF的相對(duì)估計(jì)誤差就急劇增大。KF、MKF、MCC-KF及本文算法ORKF的性能相差無(wú)幾,其在各個(gè)SOC區(qū)間內(nèi)的相對(duì)誤差都小于RCFKF,同時(shí)又都大于本文算法ORKF-VB的相對(duì)誤差。在表3中,本文算法ORKF-VB具有最小的RMSE。綜上,在鋰電池SOC估計(jì)問(wèn)題中,本文算法ORKF-VB是性能最佳的估計(jì)方法。

7 結(jié) 論

針對(duì)存在異常觀測(cè)值和/或未知觀測(cè)噪聲環(huán)境中的線性動(dòng)態(tài)濾波問(wèn)題,本文利用Huber損失函數(shù)代替推導(dǎo)卡爾曼濾波器MAP準(zhǔn)則中觀測(cè)誤差的l2范數(shù),構(gòu)造了一種新的準(zhǔn)則函數(shù),并由此推導(dǎo)出了一種ORKF算法。由于Huber函數(shù)兼顧了l1范數(shù)的魯棒性,由此本文推導(dǎo)出的卡爾曼濾波器對(duì)異常值具有魯棒性。

當(dāng)觀測(cè)噪聲分布未知時(shí),本文將噪聲建模為具有未知參數(shù)的GMM,并依據(jù)VB思想,引入一種貪婪EM算法,進(jìn)一步提出了ORKF-VB算法。

仿真驗(yàn)證和實(shí)際驗(yàn)證證明了本文分析結(jié)果的有效性。同時(shí)也表明:在含異常值和統(tǒng)計(jì)分布未知的觀測(cè)噪聲環(huán)境中,本文提出算法性能優(yōu)于現(xiàn)有文獻(xiàn)報(bào)道的魯棒類(lèi)KF算法。