王國慶 楊春雨 馬 磊 代 偉

以卡爾曼濾波 (Kalman filter,KF) 為代表的狀態(tài)估計(jì)算法是信息融合、導(dǎo)航定位、目標(biāo)跟蹤、智能電網(wǎng)等領(lǐng)域的關(guān)鍵技術(shù)[1-3].經(jīng)典的KF 針對線性高斯系統(tǒng)設(shè)計(jì),能夠獲得多種統(tǒng)計(jì)意義下的最優(yōu)估計(jì)值.隨后KF 擴(kuò)展到非線性系統(tǒng)的狀態(tài)估計(jì),代表性算法包括擴(kuò)展卡爾曼濾波 (Extended KF,EKF)、無跡卡爾曼濾波 (Unscented KF,UKF) 和容積卡爾曼濾波 (Cubature KF,CKF).EKF 直接采用一階泰勒展開將非線性函數(shù)線性化,當(dāng)非線性程度較強(qiáng)時(shí)估計(jì)精度和穩(wěn)定性難以保證.UKF 和CKF 均是采用確定性采樣點(diǎn)來計(jì)算非線性積分,從形式上講CKF 是UKF 算法的一種特例,兩者較EKF 有更好的估計(jì)精度[4].以上算法均是針對高斯噪聲系統(tǒng)設(shè)計(jì),然而在很多具體應(yīng)用中傳感器噪聲會(huì)偏離高斯分布.在室內(nèi)定位、水聲導(dǎo)航、機(jī)動(dòng)目標(biāo)跟蹤、電力系統(tǒng)等應(yīng)用中[5-8],經(jīng)常會(huì)出現(xiàn)受到大量野值、干擾等引起傳感器具有厚尾特性非高斯噪聲的情況.直接使用這些針對高斯噪聲設(shè)計(jì)的濾波方法會(huì)導(dǎo)致估計(jì)精度下降甚至算法發(fā)散[8].

針對非高斯噪聲系統(tǒng)狀態(tài)估計(jì)的經(jīng)典方法包括利用蒙特卡洛思想的粒子濾波和基于有限高斯分布逼近非高斯噪聲的高斯求和濾波算法[6],由于計(jì)算復(fù)雜度高,限制了兩者使用范圍.近年來,針對非高斯噪聲系統(tǒng)狀態(tài)估計(jì)算法的研究主要包括基于距離[9]和基于分布[10]兩種思路[11].基于距離的方法包括基于極大相關(guān)熵的濾波算法[12]、基于最小誤差熵的濾波算法[5,13]、基于M 估計(jì)的濾波算法[9]、基于統(tǒng)計(jì)相似度的魯棒濾波方法[14-15]等.以基于極大相關(guān)熵的濾波算法為例,此類方法利用相關(guān)熵概念替代常用的二范數(shù)形式來定義代價(jià)函數(shù),據(jù)此求解得到針對線性系統(tǒng)[12]、非線性系統(tǒng)[16-18]和傳感器網(wǎng)絡(luò)系統(tǒng)[19]的各類濾波算法.此類算法不針對特定類型的非高斯分布,但是其抑制非高斯噪聲的能力與核寬度有關(guān),且缺乏嚴(yán)格的自適應(yīng)更新依據(jù),其估計(jì)精度提升有限.與之類似,基于最小誤差熵和基于M估計(jì)的濾波算法均有對估計(jì)性能影響顯著的參數(shù)需要調(diào)整,基于統(tǒng)計(jì)相似度的濾波方法需要考慮相似度函數(shù)選擇的問題.第二種思路是采用典型的分布來建模常見的厚尾非高斯噪聲,常用的有學(xué)生t 分布[20-22]、多變量Laplace 分布[11,23]、橢球分布[24-25]等,一般是利用此類分布的高斯尺度混合表達(dá)形式借助變分貝葉斯等方法得到后驗(yàn)估計(jì).

根據(jù)中心極限定理,常見傳感器的噪聲可以建模為高斯分布,且在工程實(shí)踐中一般能滿足需求.但是在特定的場景下傳感器不可避免面臨各種異常擾動(dòng)、使用環(huán)境變化等引起噪聲分布發(fā)生改變.以基于傳感器網(wǎng)絡(luò)的目標(biāo)跟蹤為例,在開闊的理想場景下,傳感器噪聲建模為高斯分布是合適的;當(dāng)目標(biāo)進(jìn)入室內(nèi)、森林、地下等復(fù)雜場景,在傳感器和目標(biāo)之間存在干擾、遮擋、折射、多路徑等情況下傳感器測量誤差會(huì)增大,存在大量異常噪聲,此時(shí)可以建模為厚尾非高斯分布;當(dāng)目標(biāo)再次進(jìn)入空曠的環(huán)境中,則可以繼續(xù)利用高斯噪聲假設(shè).能夠適應(yīng)這種非平穩(wěn)厚尾非高斯噪聲的魯棒狀態(tài)估計(jì)算法更具有一般的實(shí)用價(jià)值.現(xiàn)有的魯棒狀態(tài)估計(jì)算法理論上能夠適應(yīng)這種變化,但是很難在兩種情況下均取得較好的估計(jì)效果.針對線性系統(tǒng),文獻(xiàn)[26]和文獻(xiàn)[27]考慮將這種噪聲建模為高斯-學(xué)生t 混合分布,文獻(xiàn)[28]進(jìn)一步考慮了噪聲存在偏態(tài)的情形,但是其混合函數(shù)的精確后驗(yàn)分布沒有給出.針對非線性系統(tǒng),文獻(xiàn)[29]中采用了高斯分布和學(xué)生t 分布來混合建模非平穩(wěn)厚尾量測噪聲,利用UKF 計(jì)算非線性積分項(xiàng),但是算法沒有考慮量測噪聲矩陣的變化,表征高斯分布和厚尾噪聲的切換參數(shù)也缺乏自適應(yīng)更新策略.

針對現(xiàn)有方法存在的問題,需要選擇適當(dāng)?shù)姆植急ＷC其既能建模不同厚尾程度的非高斯噪聲又能精確得到其后驗(yàn)參數(shù)概率密度函數(shù).此外還需要解決如何自適應(yīng)更新噪聲參數(shù)以及切換參數(shù)的問題.考慮到多種常見厚尾分布為廣義雙曲分布特例且其尺度混合函數(shù)為共軛分布的特性[30-32],本文基于廣義雙曲分布構(gòu)建高斯-廣義雙曲混合分布來研究帶非平穩(wěn)厚尾量測噪聲的非線性系統(tǒng)狀態(tài)估計(jì)問題,通過對超參數(shù)選擇合適的先驗(yàn)分布并借助變分貝葉斯學(xué)習(xí)實(shí)現(xiàn)參數(shù)的自適應(yīng)更新.本文主要貢獻(xiàn)包括:1) 通過引入伯努利變量構(gòu)建高斯-廣義雙曲混合分布建模非平穩(wěn)厚尾噪聲,利用廣義雙曲分布的高斯分層形式建立系統(tǒng)的概率形式模型;2) 選取合適的先驗(yàn)分布,利用變分貝葉斯方法實(shí)現(xiàn)系統(tǒng)狀態(tài)、噪聲參數(shù)以及切換參數(shù)后驗(yàn)分布的聯(lián)合估計(jì),建立了此類問題的估計(jì)框架,現(xiàn)有的幾種魯棒濾波算法[22-23,29]均是其特例;3) 基于傳感器網(wǎng)絡(luò)的機(jī)器人目標(biāo)跟蹤仿真實(shí)驗(yàn)表明,本文所提算法與現(xiàn)有方法相比有更好的估計(jì)精度和數(shù)值穩(wěn)定性.

1 問題描述

考慮如下一般的非線性加性噪聲系統(tǒng)

其中,xi ∈Rm表示系統(tǒng)在i時(shí)刻的狀態(tài)向量,過程噪聲wi ∈Rm是均值為零、方差為Qi的白噪聲,zi ∈Rn表示傳感器在i時(shí)刻的量測數(shù)據(jù),狀態(tài)轉(zhuǎn)移函數(shù)fi(·) 和量測函數(shù)hi(·) 為維數(shù)適當(dāng)?shù)囊阎蔷€性函數(shù),量測噪聲vi是均值為零、方差為Ri的白噪聲,且與過程噪聲wi不相關(guān).當(dāng)過程噪聲和量測噪聲均服從高斯分布時(shí),(1)和(2)可以表示為條件概率形式[29]

其中,p(·)表示概率密度函數(shù).本文考慮量測噪聲為非平穩(wěn)厚尾噪聲時(shí)的狀態(tài)估計(jì)問題,即量測噪聲在厚尾非高斯噪聲和高斯噪聲之間切換.在利用高斯-廣義雙曲混合分布建模此類噪聲的基礎(chǔ)上,進(jìn)而利用變分貝葉斯的方法求解狀態(tài)和參數(shù)后驗(yàn)估計(jì),建立此類問題的卡爾曼濾波算法框架.

2 高斯-廣義雙曲混合分布

本節(jié)首先介紹雙曲分布,用來建模一般的厚尾非高斯噪聲;進(jìn)而構(gòu)建高斯-廣義雙曲混合分布建模這類在高斯噪聲和厚尾噪聲之間切換的非平穩(wěn)厚尾量測噪聲.

2.1 廣義雙曲分布

如果一個(gè)變量y∈Rny服從廣義雙曲分布,則其概率密度函數(shù)可以表示為如下形式[30-31]

其中,μ是與變量y同維度的向量,P∈Rny×ny為正定矩陣,混合函數(shù)p(τ) 服從廣義逆高斯分布GIG(τ;δ,ω,η),其參數(shù)取值以及詳細(xì)特性在附錄A 中給出,逆高斯分布、伽瑪分布和逆伽瑪分布均是其特例.

此外,p(y|τ) 是均值為 (μ+τγ)、方差為τP的高斯分布,即

由上式可以建立廣義雙曲分布與高斯分布的聯(lián)系.

廣義雙曲分布可以借助 (μ+τγ) 建模偏態(tài)特性,本文為簡單起見,僅考慮μ=0 且γ=0 的零均值無偏態(tài)的情景,此時(shí)廣義雙曲分布記作GH(y;δ,ω,η,P).本文提出的算法框架也很容易擴(kuò)展到噪聲存在偏態(tài)的情形.

p(y) 分布特性取決于其混合函數(shù)p(τ). 當(dāng)p(τ)取不同的參數(shù)時(shí),廣義雙曲分布將退化成多種常見分布.除高斯分布外,還包括正態(tài)逆高斯分布、雙曲分布、K 分布 (又稱廣義Lapalce 分布)、廣義雙曲學(xué)生t 分布 (常見的學(xué)生t 分布是其特例) 等多種具有厚尾特性的分布 (詳見表1)[30-32].當(dāng)參數(shù)δ,ω和η變化時(shí),廣義雙曲分布的概率密度函數(shù)形狀可以參閱文獻(xiàn)[31].可以發(fā)現(xiàn)其形狀對于參數(shù)δ最為敏感,選取適當(dāng)?shù)膮?shù)該分布能表征各種不同厚尾程度的非高斯噪聲.基于廣義雙曲分布建模所得到的結(jié)果適用于一系列的厚尾分布噪聲,結(jié)果具有普適性.

表1 廣義雙曲分布的幾種特殊分布Table 1 Several special cases of generalized hyperbolic distribution

2.2 高斯-廣義雙曲混合分布

為建模非平穩(wěn)厚尾量測噪聲,引入輔助變量si,量測噪聲則可以表示為如下高斯-廣義雙曲混合分布形式

其中,si ∈{0, 1}. 當(dāng)si=1 時(shí),表示量測噪聲服從典型的高斯分布;當(dāng)si=0 時(shí),表示量測噪聲服從廣義雙曲分布,用來建模大量的尖峰脈沖干擾等誘導(dǎo)的厚尾類型非高斯噪聲.采用此建模形式,可以很方便地在高斯分布與厚尾非高斯噪聲之間實(shí)現(xiàn)切換.

由于si ∈{0, 1}是離散的二值變量,因此可以利用參數(shù)為λi的伯努利分布來建模[26-27],其概率質(zhì)量函數(shù)可以表示為

伯努利分布參數(shù)λi的共軛先驗(yàn)為貝塔分布[23]

其中,0＜κ0＜1 和 1-κ0為貝塔分布的參數(shù).

根據(jù)廣義雙曲分布的特性,其混合函數(shù)τi服從廣義逆高斯分布

進(jìn)一步考慮量測噪聲矩陣Ri存在不確定性,對于這種正定矩陣常用逆威沙特分布來建模[3]

其中,v0和V0分別為自由度參數(shù)和逆尺度矩陣.

綜上所述,(3) 和 (7) 分別是本文所研究系統(tǒng)概率形式的狀態(tài)方程和量測方程.為便于理解,圖1中給出了相應(yīng)變量關(guān)系的圖模型.

圖1 本文系統(tǒng)的圖模型Fig.1 Graph model of the system used in this paper

3 本文提出的算法

為得到遞推的狀態(tài)估計(jì)算法,首先假設(shè)上一時(shí)刻的狀態(tài)后驗(yàn)為高斯分布p(xi-1|Zi-1)=N(xi-1;其中為i-1 時(shí)刻狀態(tài)估計(jì)值,Pi-1|i-1為誤差協(xié)方差矩陣,Zi-1={z1,z2,···,zi-1}表示量測數(shù)據(jù)的集合.算法初始化時(shí)假設(shè)已知狀態(tài)均值和誤差協(xié)方差矩陣P0|0.接下來將利用模型信息和i時(shí)刻的量測信息zi設(shè)法求得當(dāng)前i時(shí)刻后驗(yàn)估計(jì)與常規(guī)卡爾曼濾波算法一致,算法也包括預(yù)測更新和量測更新兩步.

3.1 預(yù)測更新

由于上一時(shí)刻狀態(tài)后驗(yàn)為高斯分布,考慮到系統(tǒng)噪聲也服從高斯分布,則算法的時(shí)間更新與常規(guī)的非線性估計(jì)方法一致.狀態(tài)一步預(yù)測的概率密度函數(shù)表示為[23,29]

其中,狀態(tài)一步預(yù)測可以通過下式得到

其誤差協(xié)方差矩陣Pi|i-1為

3.2 量測更新

根據(jù)貝葉斯公式,系統(tǒng)狀態(tài)以及未知參數(shù)后驗(yàn)估計(jì)可以根據(jù)下式求得

其中,Θi={xi,Ri,τi,si,λi}.由于量測噪聲不再服從高斯分布且未知變量之間存在耦合,上式解析解不存在,不能直接進(jìn)行求解.

針對系統(tǒng)狀態(tài)與參數(shù)聯(lián)合估計(jì)問題,常用的方法有EM 算法和變分貝葉斯方法,其中變分貝葉斯方法能適應(yīng)隱變量維數(shù)更高的狀態(tài)估計(jì)問題,且能夠得到更好的估計(jì)結(jié)果[1],也是EM 算法和極大后驗(yàn)估計(jì)的推廣[31].因此本文借助變分貝葉斯方法進(jìn)行求解,將聯(lián)合后驗(yàn)分布近似為每個(gè)變量后驗(yàn)分布的乘積[1,3]

其中,q(·) 表示變量的后驗(yàn)概率密度函數(shù)或者概率質(zhì)量函數(shù).

每個(gè)變量后驗(yàn)分布通過最小化其與真實(shí)后驗(yàn)之間的Kullback-Leibler 散度 (Kullback-Leibler divergence,KLD) 來求解

根據(jù)變分貝葉斯框架,上式的最優(yōu)解為[3,23,31]

其中,?表示集合 Θi中的任意一個(gè)元素,Θ-i?表示集合中除?以外其他元素的集合,E (·) 表示期望,c?表示與?有關(guān)的常數(shù),l n(·) 表示對數(shù)運(yùn)算.

根據(jù)系統(tǒng)概率模型,結(jié)合貝葉斯公式和條件概率公式,聯(lián)合概率密度函數(shù)p(Θi,Zi) 可以分解為以下乘積形式

代入相關(guān)分布的概率密度函數(shù)表達(dá)式,可以得到(19)的對數(shù)形式

其中,t r(·) 表示矩陣的跡,|·| 表示行列式.

借助 (18) 和 (20),可以利用變分貝葉斯迭代分別求解各個(gè)未知變量的后驗(yàn)分布.接下來將給出各個(gè)變量具體的后驗(yàn)分布形式以及參數(shù)計(jì)算方法.

3.2.1 狀態(tài)變量 xi 的后驗(yàn)分布

根據(jù) (18) 和 (20),令?=xi,合并與變量xi有關(guān)的參數(shù),其余待求解變量用其上一次迭代的期望替代,可以得到狀態(tài)xi后驗(yàn)分布的對數(shù)表達(dá)式

其中,上標(biāo) (l) 表示第l次固定點(diǎn)迭代,cxi表示與xi相關(guān)的常數(shù) (下同).

整理可得

其中,量測預(yù)測值與常規(guī)非線性濾波方法一致,即

3.2.2 變量 R i 的后驗(yàn)分布

令?=Ri,利用 (18) 和 (20) 可以得到

根據(jù)后驗(yàn)分布的對數(shù)形式,可以看出其依然為逆威沙特分布.通過參數(shù)匹配可以得到變量Ri的后驗(yàn)分布為

利用威沙特分布的均值特性,可以得到如下期望[3]

3.2.3 變量 τi 的后驗(yàn)分布

令?=τi,代入 (18) 和 (20),可以得到

由于廣義逆高斯分布為共軛分布[30-31],通過匹配參數(shù),可以發(fā)現(xiàn)其后驗(yàn)分布依然為廣義逆高斯分布,即

相應(yīng)的參數(shù)分別為

3.2.4 變量 si 的后驗(yàn)分布

令?=si,根據(jù) (18) 和 (20) 有

根據(jù) l nq(l+1)(si) 的形式,可以發(fā)現(xiàn)si的后驗(yàn)分布依然服從伯努利分布,其后驗(yàn)概率質(zhì)量函數(shù)表示為

利用伯努利分布的特性,si的期望為

3.2.5 變量 λi 的后驗(yàn)分布

令?=λi,根據(jù) (18) 和 (20) 可得

匹配后驗(yàn)分布,可以發(fā)現(xiàn)λi依然為貝塔分布,即

其中,ψ(·) 表示雙伽瑪函數(shù).以上期望用于 (47) 和(48) 的計(jì)算.

至此,已經(jīng)得到所有未知參數(shù)的后驗(yàn)分布以及參數(shù)計(jì)算方法,完整算法將在第 3.3 節(jié)給出.

3.3 算法流程

本文提出的基于高斯-廣義雙曲混合分布的非線性卡爾曼濾波算法在i時(shí)刻完整的步驟如下:

步驟 1.輸入?yún)?shù):v0,δ0,ω0,η0,κ0;

步驟 2.預(yù)測更新: 根據(jù) (13) 和 (14) 計(jì)算和Pi|i-1;

步驟 3.量測更新:

1) 迭代初始化:

算法利用雙曲分布建模厚尾非高斯噪聲,通過引入二值變量si來表征系統(tǒng)在高斯噪聲和厚尾噪聲之間的切換.由于雙曲分布的混合概率密度函數(shù)τi為共軛的廣義逆高斯分布,τi的后驗(yàn)估計(jì)是封閉的,利用精確的后驗(yàn)密度得到的估計(jì)結(jié)果更為準(zhǔn)確.此外算法考慮了Ri,si和λi的自適應(yīng)更新,算法的魯棒性更好.算法中存在多處非線性積分項(xiàng),可以利用無跡變換、容積變換等基于有限采樣點(diǎn)的方式來計(jì)算[3],或者借助一階泰勒展開的思路計(jì)算非線性積分的均值和方差信息[16].換而言之,本文算法的具體實(shí)現(xiàn)可以借助EKF、UKF 等常規(guī)非線性濾波算法中處理非線性的方法來處理非線性.本文算法顯然也適用于線性系統(tǒng).

另一方面,由于利用廣義雙曲分布建模厚尾噪聲,現(xiàn)有常用的幾種厚尾分布均是其特例,此外本文算法還考慮了參數(shù)的自適應(yīng)估計(jì),因此現(xiàn)有的幾種常見的魯棒濾波和自適應(yīng)濾波算法均是本文算法的特例.具體而言,在si=1 條件下,取η0=0,δ0＜0,ω0＞0,且δ0=-ω0,則本文提出的算法退化為文獻(xiàn)[22]中基于學(xué)生t 分布建模非高斯量測噪聲的魯棒濾波算法;當(dāng)ω0=0,η0=1,δ0=1 時(shí),則雙曲分布退化為多變量Laplcae 分布,本文提出的算法退化為文獻(xiàn)[23]中基于多變量Laplace 分布的RSEML 濾波算法.當(dāng)si=0,本文算法退化為文獻(xiàn)[33]中高斯噪聲下自適應(yīng)估計(jì)量測噪聲的VBAKF-R算法.文獻(xiàn)[29]中基于切換思想的SMKF 算法是本文算法中si=0.5,η0=0,δ0＜0,ω0＞0,Ri為固定值的一種特例.

所提算法的計(jì)算量主要受固定點(diǎn)迭代次數(shù)N影響,后面的仿真實(shí)驗(yàn)可以看到僅需要幾次迭代即可收斂,因此算法的計(jì)算量相對經(jīng)典CKF 和UKF算法增加較少.本文算法中參數(shù)v0和V0可以參考文 獻(xiàn)[23] 設(shè) 為v0=ρ和其中常數(shù)ρ ＞為量測噪聲主導(dǎo)方差矩陣.后文仿真結(jié)果表明,由于算法能夠自適應(yīng)調(diào)整噪聲參數(shù)的后驗(yàn)估計(jì),在一定范圍內(nèi)本文所提算法對于初始參數(shù)具有較好的魯棒性,具體而言,算法對于δ0,ω0和η0的初始值具有一定的魯棒性,對于參數(shù)κ0和E(0)(si)不敏感.實(shí)際使用時(shí),可以按照表1 中典型分布選取δ0,ω0和η0的值,其余參數(shù)可以參考后文的仿真效果選取.

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

4.1 仿真設(shè)置

為驗(yàn)證本文所提算法的有效性和優(yōu)越性,下面以一個(gè)二維平面內(nèi)目標(biāo)跟蹤的例子來比較本文算法與同類算法的估計(jì)性能,并考察算法參數(shù)對估計(jì)效果的影響.

考慮利用3 個(gè)測距傳感器 (如超聲波傳感器、UWB 模塊等) 跟蹤一個(gè)機(jī)器人的場景 (如圖2 所示)[19,34].

圖2 基于傳感器網(wǎng)絡(luò)的機(jī)器人跟蹤示意圖Fig.2 The illustration of tracking a robot with the sensor network

二維平面內(nèi)跟蹤問題的系統(tǒng)方程表示為[19,23]

其中,取q=10.

量測方程為目標(biāo)到3 個(gè)傳感器的測距信息,即[19,23,34]

其 中,(x(1)=0 m,y(1)=0 m),(x(2)=1 000 m,y(2)=0 m),(x(3)=0 m,y(3)=1 000 m) 分別為傳感器的坐標(biāo).

為便于比較不同算法的估計(jì)效果,定義常用的位置和速度的均方誤差 (Root mean square error,RMSE) 來度量估計(jì)誤差.以位置RMSE 為例,其定義為[19,22-23]

其中,上標(biāo)s表示第s次蒙特卡洛仿真實(shí)驗(yàn),為i時(shí)刻真實(shí)的位置,表示i時(shí)刻濾波算法的估計(jì)結(jié)果,M=500 是蒙特卡洛仿真的總實(shí)驗(yàn)次數(shù),R MSEvel的定義也類似.

為公平起見,各類算法中的非線性積分項(xiàng)均采用容積規(guī)則 (一種特殊的無跡變換,故部分對比算法后綴用UKF 結(jié)尾) 來計(jì)算.需要變分貝葉斯迭代的算法默認(rèn)迭代次數(shù)N=10,仿真時(shí)長設(shè)為 4 00 s.與同類算法[19,23,26-28]類似,量測噪聲按照如下方式生成

其中,w1,i～U(10, 100) 表示服從 [ 10, 100] 之間連續(xù)的均勻分布;w2,i=1+0.5 cos(πi) 用來建模噪聲方差矩陣時(shí)變; w.p.表示出現(xiàn)概率;R0=diag{100 m2, 50 m2, 200 m2}.后續(xù)仿真中量測噪聲主導(dǎo)方差矩陣設(shè)為即在 [ 0 s, 300 s] 內(nèi)模擬典型的厚尾非高斯噪聲,量測噪聲70%概率為 N (0,w2,iR0),30%概率為N(0,w1,iR0). 在 (300 s, 400 s] 內(nèi)量測噪聲服從參數(shù)時(shí)變的高斯分布 N (0,w2,iR0).為便于理解,以一維情況為例,圖3 中給出了噪聲的幅值序列以及各段的概率密度函數(shù),其中 [ 0 s, 300 s] 為厚尾分布,(300 s, 400 s]為高斯分布,[ 0 s, 400 s] 為非平穩(wěn)厚尾分布.可以看到,盡管非平穩(wěn)噪聲的幅值發(fā)生較為明顯的變化,但是概率密度函數(shù)變化不明顯,本文算法考慮到此種噪聲分布的切換,能夠得到更好的估計(jì)效果.

圖3 仿真中產(chǎn)生一維噪聲的幅值以及概率密度函數(shù)Fig.3 The amplitude and probability density functions of the one-dimensional noise used in the simulation

4.2 算法對比

接下來比較本文提出算法與同類算法的估計(jì)精度.實(shí)驗(yàn)中對比以下10 種算法: 現(xiàn)有的6 種算法,包括利用量測噪聲主導(dǎo)方差矩陣的經(jīng)典UKF 算法、文獻(xiàn)[22]中的STUKF 算法、文獻(xiàn)[18]中的MCUKF算法、文獻(xiàn)[23]中的RSE-ML 算法、文獻(xiàn)[29]中的SEUKF 算法、利用真實(shí)量測噪聲方差矩陣的最優(yōu)UKF 算法 (OUKF);本文提出算法的4 種典型特例,即算法中廣義雙曲分布分別為正態(tài)逆高斯分布(δ0=-0.5,ω0=2,η0=2,記作G-GHUKF1)、雙曲分布 (δ0=1,ω0=2,η0=2,記作G-GHUKF2)、K 分布 (δ0=2,ω0=0,η0=2,記作G-GHUKF3)和廣義雙曲學(xué)生t 分布 (δ0=-2,ω0=2,η0=0,記作G-GHUKF4).STUKF 算法中先驗(yàn)參數(shù)設(shè)置為exp(-5); RSE-ML 算法中ρ2=5;MCUKF 算法中核寬度σ=8; SEUKF 算法中選擇a0=2,b0=2,λ0=0.5(為比較更為公平,逆伽瑪分布的先驗(yàn)信息與本文算法選取一致);本文提出的幾種算法中設(shè)置為比較不同算法的計(jì)算復(fù)雜度,本節(jié)還給出了不同算法總的運(yùn)行時(shí)間.仿真采用MATLAB R2020b 軟件進(jìn)行,所用的臺(tái)式機(jī)電腦CPU 型號為英特爾酷睿i9-10900,主頻為2.80 GHz,內(nèi)存為16 GB.

仿真結(jié)果在圖4 和圖5 中給出.利用噪聲主導(dǎo)方差信息的UKF 在量測噪聲為非高斯噪聲時(shí)估計(jì)誤差比高斯噪聲時(shí)顯著增大,因此在實(shí)際使用中當(dāng)非高斯噪聲嚴(yán)重時(shí)必須予以考慮.現(xiàn)有MCUKF 在非高斯噪聲時(shí)能夠顯著改善UKF 的估計(jì)精度,在噪聲恢復(fù)為參數(shù)時(shí)變的高斯噪聲時(shí)其估計(jì)效果依然較好,但是由于該算法沒有充分利用噪聲統(tǒng)計(jì)特性,且其核寬度難以自適應(yīng)更新,其估計(jì)精度還有提升空間.STUKF 算法和RSE-ML 算法分別采用了學(xué)生t 分布和多變量Laplace 分布來建模厚尾非高斯噪聲,兩者在非高斯噪聲情況下估計(jì)精度較好,但是在噪聲切換到正常的高斯噪聲時(shí),算法估計(jì)精度相對較差.SEUKF 算法建模中考慮了噪聲切換,但是由于其不能自適應(yīng)更新量測噪聲方差矩陣以及噪聲切換信息,在本文的場景下容易出現(xiàn)數(shù)值計(jì)算錯(cuò)誤,圖中僅給出了前51 s 的誤差曲線 (后續(xù)部分由于頻繁出現(xiàn)數(shù)值計(jì)算不穩(wěn)定,軟件自動(dòng)剔除).本文所提算法的4 種特例在非高斯噪聲情況和高斯噪聲下均能夠得到較好的估計(jì)精度,特別是在由非高斯噪聲切換到參數(shù)時(shí)變的高斯噪聲后,算法能夠自動(dòng)適應(yīng)這種分布變化,精度顯著高于同類的STUKF和RSE-ML 算法,體現(xiàn)所提算法的優(yōu)越性.此外本文所提的算法當(dāng)雙曲分布參數(shù)選為4 類不同的具體分布時(shí),估計(jì)精度相差不大.

圖4 本文所提算法與同類方法的RMSEposFig.4 The R MSEpos of the proposed algorithms and related ones

圖5 本文所提算法與同類方法的RMSEvelFig.5 The R MSEvel of the proposed algorithms and related ones

根據(jù)表2 中不同算法的計(jì)算時(shí)間來看,本文所提算法計(jì)算耗時(shí)與同類基于變分貝葉斯的方法處于同一量級.由于每個(gè)采樣時(shí)刻進(jìn)行了10 次固定點(diǎn)迭代,因此本文算法計(jì)算耗時(shí)大約是經(jīng)典UKF 的10 倍,后續(xù)仿真可知利用更少的迭代次數(shù)依然能夠得到較好的估計(jì)效果,相應(yīng)的計(jì)算耗時(shí)也會(huì)隨之減小.另外G-GHUKF4 中由于參數(shù)選擇的原因,后驗(yàn)分布不需要計(jì)算第二類修正貝塞爾函數(shù),因此計(jì)算量相比其他幾種情況稍低.

表2 不同算法總運(yùn)行時(shí)間Table 2 The total simulation time of different algorithms

4.3 算法參數(shù)選擇

接下來將詳細(xì)探討本文所提算法參數(shù)對于估計(jì)效果的影響.

首先,考察變分貝葉斯迭代次數(shù)對于估計(jì)精度的影響.固定其他參數(shù) (δ0=1,ω0=2,η0=2),迭代次數(shù)N分別取3、6、9、12、15,估計(jì)結(jié)果在圖6和圖7 中給出.可以看到當(dāng)N≥6 時(shí)估計(jì)誤差重合,即算法僅需要幾次就可以收斂,估計(jì)精度不再提高,因此變分貝葉斯迭代帶來的計(jì)算負(fù)擔(dān)較小.

圖6 本文算法在迭代次數(shù) N 不同時(shí)的RMSEposFig.6 The R MSEpos of the proposed algorithms with different iteration numberN

圖7 本文算法在迭代次數(shù) N 不同時(shí)的RMSEvelFig.7 The R MSEvel of the proposed algorithms with different iteration numberN

其次,固定δ0=-0.5,令ω0=η0分別取1、2、3、4、5,仿真結(jié)果在圖8 和圖9 中給出.圖中可以看到,當(dāng)兩者參數(shù)過小 (均為1) 時(shí)估計(jì)精度略微降低,兩者參數(shù)設(shè)置更大時(shí)估計(jì)結(jié)果基本不變.

圖8 本文算法在 ω0 和 η0 不同時(shí)的RMSEposFig.8 The R MSEpos of the proposed algorithms with different ω0 andη0

圖9 本文算法在 ω0 和 η0 不同時(shí)的RMSEvelFig.9 The R MSEvel of the proposed algorithms with different ω0 andη0

隨后,固定ω0=2 和η0=0,δ0分別取-2、-3、-4、-5、-6,估計(jì)誤差在圖10 和圖11 中給出,可以看到當(dāng)δ0過大時(shí)估計(jì)精度會(huì)略微下降.

圖10 本文算法在 δ0 不同時(shí)的RMSEposFig.10 The R MSEpos of the proposed algorithms with different δ0

圖11 本文算法在 δ0 不同時(shí)的RMSEvelFig.11 The R MSEvel of the proposed algorithms with different δ0

最后,還考察了 E(0)(si) 和κ0對算法估計(jì)結(jié)果的影響,改變兩者取值對于估計(jì)基本沒有影響,因此不再給出其仿真結(jié)果.

5 結(jié)論

針對帶非平穩(wěn)厚尾量測噪聲非線性系統(tǒng)的狀態(tài)估計(jì)問題,本文引入高斯-廣義雙曲混合分布進(jìn)行噪聲建模,利用變分貝葉斯方法求解其狀態(tài)以及各種噪聲參數(shù)的后驗(yàn)分布,得到了基于高斯-廣義雙曲混合分布的非線性卡爾曼濾波算法框架.分析表明,現(xiàn)有的幾種魯棒 (自適應(yīng)) 濾波算法均是其特例,且本文提出的算法框架相比現(xiàn)有方法具有更好的估計(jì)精度和數(shù)值穩(wěn)定性.未來的研究方向包括考慮量測噪聲存在偏態(tài)、系統(tǒng)噪聲也存在非平穩(wěn)厚尾噪聲、分布式的多傳感器信息融合等復(fù)雜場景.

附錄A.廣義逆高斯分布及其特性

廣義逆高斯分布GIG(τ;δ,ω,η) 的概率密度函數(shù)為[30-31]

可以看到這是一種三參數(shù)連續(xù)型分布,其參數(shù)選取需要滿足[30]:當(dāng)δ<0 時(shí),ω >0,η≥0; 當(dāng)δ=0 時(shí),ω >0,η>0; 當(dāng)δ >0時(shí),ω≥0,η>0.

常見的逆高斯分布、伽瑪分布和逆伽瑪分布分別是廣義逆高斯分布取δ=-0.5,ω=0 和η=0 時(shí)的特例.對于一般的廣義逆高斯分布p(τ) 而言,當(dāng)ω>0 且η>0 時(shí),其r階原點(diǎn)矩可以表示為[30-31]

當(dāng)ω=0 或η=0 時(shí),上式不再適用,此時(shí)p(τ) 分別退化為伽瑪分布和逆伽瑪分布[30-31].

當(dāng)η=0,δ<0,ω>0 時(shí),廣義逆高斯分布退化為逆伽瑪分布,即

其中,伽瑪分布的兩個(gè)參數(shù)分別為形狀參數(shù)和尺度參數(shù).當(dāng)τ服從伽瑪分布時(shí),則τ-1服從相同形狀參數(shù)和尺度參數(shù)的逆伽瑪分布.由此可以得到本文中用到的τ和τ-1的期望分別為(δ/=-1)

當(dāng)ω=0,η>0,δ>0 時(shí),廣義逆高斯分布退化為伽瑪分布,即

則可以得到(δ/=1)