胡玉梅 潘 泉 胡振濤 郭 振,5

非線性濾波器設(shè)計(jì)與實(shí)現(xiàn)因其普適性及重要性一直是國內(nèi)外學(xué)者研究的熱點(diǎn)問題,近年來非線性狀態(tài)估計(jì)理論已成功應(yīng)用于陸、海、空、天中運(yùn)動(dòng)目標(biāo)的預(yù)警與防御,智能交通的精確導(dǎo)航與制導(dǎo),無人機(jī)定位與遙感監(jiān)測、工業(yè)過程監(jiān)控與故障診斷等眾多領(lǐng)域[1-3].考慮到運(yùn)動(dòng)目標(biāo)跟蹤系統(tǒng)機(jī)動(dòng)、隱身等復(fù)雜多變的人為對抗特征以及非視距、干擾、遮擋等環(huán)境因素?zé)o法避免,其系統(tǒng)建模、估計(jì)與辨識(shí)過程中越來越無法回避非線性、非高斯以及參數(shù)未知等復(fù)雜系統(tǒng)特征的影響[4].在對復(fù)雜環(huán)境下運(yùn)動(dòng)目標(biāo)系統(tǒng)噪聲先驗(yàn)信息進(jìn)行建模時(shí),建模誤差存在于狀態(tài)演化模型中,并通常假設(shè)其滿足一定的參數(shù)化統(tǒng)計(jì)特性[5].然而,在實(shí)際工程中,由于先驗(yàn)建模信息的不足導(dǎo)致難以對此類參數(shù)進(jìn)行準(zhǔn)確賦值.例如,在現(xiàn)代目標(biāo)跟蹤系統(tǒng)中出現(xiàn)的欺騙、干擾、雜波、未知分布特性的量測噪聲和系統(tǒng)噪聲等情形,尤其在非合作運(yùn)動(dòng)目標(biāo)強(qiáng)機(jī)動(dòng)場景中,因難以對其運(yùn)動(dòng)過程進(jìn)行精細(xì)化建模,常造成目標(biāo)跟蹤航跡間斷甚至無法正常起始的現(xiàn)象.

針對系統(tǒng)噪聲未知問題,其解決思路通常選取逆伽馬分布和逆威沙特(Inverse-Wishart,IW)分布等具有統(tǒng)一參數(shù)表達(dá)形式的指數(shù)族分布作為共軛先驗(yàn)分布.S?rkk? 等[6]在變分貝葉斯 (Variational Bayes,VB)推斷框架下利用指數(shù)族分布構(gòu)建共軛先驗(yàn)分布近似未知量測噪聲的后驗(yàn)概率密度函數(shù)(Probability density function,PDF),進(jìn)而結(jié)合貝葉斯濾波機(jī)理實(shí)現(xiàn)時(shí)變噪聲方差和目標(biāo)狀態(tài)的聯(lián)合估計(jì),其濾波效果達(dá)到與交互式多模型方法相接近的估計(jì)精度,并且憑借 VB 便于結(jié)合平均場近似解耦理論將高維變量求解轉(zhuǎn)化為多個(gè)低維變量計(jì)算的特點(diǎn),使其在解決多未知擾動(dòng)問題方面更具有優(yōu)勢.隨后,文獻(xiàn)[7] 采用IW 分布近似多變量噪聲后驗(yàn)PDF 的思想,給出變分推斷框架下噪聲方差估計(jì)的一般實(shí)現(xiàn)形式.在此基礎(chǔ)上,變分貝葉斯方法在未知量測/系統(tǒng)噪聲估計(jì)方面得到了進(jìn)一步的發(fā)展[8-9],并成功推廣應(yīng)用于多目標(biāo)跟蹤環(huán)境[10-11].在文獻(xiàn)[12-13]中濾波器的構(gòu)建均建立在量測噪聲和系統(tǒng)噪聲統(tǒng)計(jì)特性未知的假設(shè)條件下,其中,文獻(xiàn)[12] 在狀態(tài)和量測擴(kuò)維基礎(chǔ)上通過采用批處理方式實(shí)現(xiàn)對未知噪聲矩陣的估計(jì),然而這種擴(kuò)維方式不可避免造成計(jì)算復(fù)雜度的急劇增加;文獻(xiàn)[13] 給出兩類噪聲統(tǒng)計(jì)特性均未知條件下噪聲方差的在線估計(jì)方法,并且為避免系統(tǒng)噪聲和狀態(tài)相互獨(dú)立假設(shè)的不合理性,采用獨(dú)立于當(dāng)前狀態(tài)的前一時(shí)刻狀態(tài)預(yù)測誤差協(xié)方差代替系統(tǒng)噪聲的統(tǒng)計(jì)特性,進(jìn)而利用平均場理論近似解耦計(jì)算邊緣后驗(yàn) PDF 的變分分布,以迭代遞推方式求解其估計(jì)變量期望的解析解,同時(shí),將狀態(tài)一步預(yù)測誤差協(xié)方差的先驗(yàn)分布建模為參數(shù)化IW 分布,以實(shí)現(xiàn)狀態(tài)后驗(yàn)概率分布的更新.

考慮到外界環(huán)境干擾以及非線性傳播等因素的影響,量測數(shù)據(jù)概率分布往往呈現(xiàn)出重尾和非對稱等非高斯特征.例如,在基于電磁信號(hào)的距離估計(jì)中,障礙物遮擋造成的非視距量測誤差往往遠(yuǎn)大于其他來源服從對稱分布的誤差,導(dǎo)致量測分布呈現(xiàn)非對稱非高斯現(xiàn)象.其次,受傳感器精度和靈敏度的限制,當(dāng)運(yùn)動(dòng)目標(biāo)發(fā)生強(qiáng)機(jī)動(dòng)時(shí),雷達(dá)極化反射不穩(wěn)定性將導(dǎo)致的量測隨機(jī)缺失現(xiàn)象,也將造成量測產(chǎn)生非高斯特征.針對量測噪聲非高斯問題的處理,文獻(xiàn)[14]采用萊斯分布構(gòu)建非共軛指數(shù)族變換點(diǎn)檢測模型,并將其成功應(yīng)用于雷達(dá)目標(biāo)跟蹤系統(tǒng)中.文獻(xiàn)[15]則提出一種通過構(gòu)建狀態(tài)演化模型和量測模型的條件矩實(shí)現(xiàn)非線性非高斯系統(tǒng)的狀態(tài)估計(jì)方法.為了進(jìn)一步提升濾波器對不同分布形狀噪聲的魯棒性,文獻(xiàn)[16] 將 Skew-t分布作為非高斯量測似然的近似分布形式,在此基礎(chǔ)上設(shè)計(jì)了一種魯棒變分貝葉斯估計(jì)器.文獻(xiàn)[17] 考慮在α散度最小化準(zhǔn)則下,利用變分分布實(shí)現(xiàn)對后驗(yàn) PDF 的近似.α散度對異常數(shù)據(jù)具有較好的抑制作用,但是也打破了對數(shù)邊緣概率與變分置信下界 (Evidence lower bound,ELBO)、KL 散度(Kullback-Leibler divergence)之和的等式約束關(guān)系,其實(shí)質(zhì)是一種偽 VB 推斷方法.考慮到學(xué)生t分布和 Skew-t分布對因量測異常導(dǎo)致的量測重尾現(xiàn)象和非對稱性具有較好魯棒性[18-19],文獻(xiàn)[20] 針對系統(tǒng)噪聲和量測噪聲均具有重尾特性的線性系統(tǒng)采用學(xué)生t分布對狀態(tài)預(yù)測概率密度函數(shù)和量測似然函數(shù)進(jìn)行建模,進(jìn)而提出了一種基于高斯-學(xué)生t混合分布的線性濾波器,與傳統(tǒng)高斯假設(shè)條件下濾波器估計(jì)效果相比進(jìn)一步提升了系統(tǒng)的狀態(tài)估計(jì)的精度和魯棒性.在系統(tǒng)噪聲未知且時(shí)變和量測噪聲重尾條件下,文獻(xiàn)[21] 構(gòu)建參數(shù)化IW 分布作為狀態(tài)一步預(yù)測誤差協(xié)方差的共軛先驗(yàn)分布,同時(shí),選取參數(shù)化學(xué)生t分布刻畫具有厚尾特點(diǎn)的量測似然函數(shù).然而,上述噪聲自適應(yīng)方法未考慮非線性濾波器自身的優(yōu)化問題和如何衡量后驗(yàn)分布近似程度的問題.

考慮到在非線性濾波器設(shè)計(jì)過程中,參數(shù)化變分分布對系統(tǒng)狀態(tài)后驗(yàn)近似的程度是提高估計(jì)精度的關(guān)鍵因素之一,基于采樣的隨機(jī)優(yōu)化和基于擬牛頓、高斯牛頓和梯度上升等確定性優(yōu)化方法的非線性濾波器相繼提出.在隨機(jī)性優(yōu)化方面,MCMC(Markov chain Monte Carlo)[22-23]和序貫蒙特卡羅(Sequential Monte Carlo,SMC)利用隨機(jī)樣本來逼近后驗(yàn)概率,并以樣本分布近似未知變量后驗(yàn)分布.隨機(jī)優(yōu)化方法的優(yōu)點(diǎn)是在大量樣本的條件下能夠保證較高精度的估計(jì),但要承擔(dān)較大的計(jì)算負(fù)擔(dān)[24].在確定性優(yōu)化方面,文獻(xiàn)[25] 以最小二乘準(zhǔn)則為目標(biāo)函數(shù),采用牛頓法推導(dǎo)給出一種迭代卡爾曼濾波器.文獻(xiàn)[26] 綜合梯度法和牛頓法提出一種無噪聲條件下的滾動(dòng)時(shí)域估計(jì)器,并證明其漸進(jìn)收斂性和穩(wěn)定性.文獻(xiàn)[27] 針對非線性狀態(tài)空間的狀態(tài)估計(jì)及其估計(jì)誤差協(xié)方差的迭代更新問題,利用正則化的非線性最小二乘思想提出一種隨機(jī)增量近端梯度算法.VB 方法通過求解參數(shù)化的目標(biāo)函數(shù),利用參數(shù)優(yōu)化結(jié)果逼近后驗(yàn)分布,是一種確定性近似方法.因此,VB 具有確定性優(yōu)化方法特有的計(jì)算量較小的優(yōu)勢,同時(shí)由于 VB 便于結(jié)合平均場理論近似解耦聯(lián)合概率分布,進(jìn)一步減小了計(jì)算代價(jià).

非線性狀態(tài)估計(jì)優(yōu)化的實(shí)質(zhì)是對多維狀態(tài)后驗(yàn)PDF 的近似逼近,并且其近似程度不能簡單地采用歐氏距離進(jìn)行度量.因此,選取合理度量準(zhǔn)則將有利于提高后驗(yàn) PDF 的近似程度,KL 散度作為分布之間“差異”的度量在后驗(yàn)分布近似性衡量中具有天然優(yōu)勢.文獻(xiàn)[17]給出狀態(tài)估計(jì)的變分迭代優(yōu)化實(shí)現(xiàn)形式,獲得α散度和 KL 散度下對后驗(yàn) PDF 更緊密的近似.同時(shí),從信息幾何角度出發(fā),概率分布是統(tǒng)計(jì)流形上點(diǎn),在一定條件下兩概率分布之間的 KL 散度與作為統(tǒng)計(jì)流形度規(guī)的 Fisher 信息滿足一定的數(shù)學(xué)關(guān)系.基于信息幾何理論,Amari[28]利用自然梯度優(yōu)化方式實(shí)現(xiàn)統(tǒng)計(jì)流形空間中目標(biāo)函數(shù)的最速梯度下降(或上升).文獻(xiàn)[29] 結(jié)合自然梯度策略和卡爾曼濾波框架設(shè)計(jì)一種非線性狀態(tài)估計(jì)方法,在文獻(xiàn)[30] 中,該方法進(jìn)一步推廣應(yīng)用于傳感器網(wǎng)絡(luò)目標(biāo)跟蹤系統(tǒng).文獻(xiàn)[31] 中作者證明了針對非線性狀態(tài)估計(jì)優(yōu)化的自然梯度方法在克拉美羅下界意義下是漸進(jìn)最優(yōu)的.考慮自然梯度優(yōu)化的優(yōu)勢,以變分置信下界最大化條件下狀態(tài)估計(jì)及其估計(jì)誤差協(xié)方差的自然梯度為切入點(diǎn),從信息幾何角度實(shí)現(xiàn)對狀態(tài)后驗(yàn)概率密度函數(shù)的“緊密”逼近,進(jìn)而提高狀態(tài)估計(jì)精度.

在過程噪聲先驗(yàn)不準(zhǔn)確和由于量測隨機(jī)丟失導(dǎo)致的量測噪聲分布非高斯的情況下,針對非線性系統(tǒng)狀態(tài)估計(jì)精度和魯棒性提升問題,本文提出一種基于自然梯度的噪聲自適應(yīng)變分貝葉斯濾波算法.本文結(jié)構(gòu)如下:第 1 節(jié)結(jié)合IW 分布和學(xué)生t分布分別實(shí)現(xiàn)對狀態(tài)預(yù)測誤差協(xié)方差和量測似然的參數(shù)化建模;第 2 節(jié)結(jié)合平均場近似和坐標(biāo)上升方法給出了變分貝葉斯框架下未知變量變分分布的迭代更新方式;第 3 節(jié)推導(dǎo)給出 ELBO 關(guān)于系統(tǒng)狀態(tài)估計(jì)及其誤差協(xié)方差的自然梯度,進(jìn)而構(gòu)建并設(shè)計(jì)一種基于自然梯度的噪聲自適應(yīng)非線性變分貝葉斯濾波器;第 4 節(jié)給出仿真驗(yàn)證與分析;第 5 節(jié)總結(jié)全文,并展望了后續(xù)研究方向.

1 預(yù)備知識(shí)

假設(shè)動(dòng)態(tài)系統(tǒng)隱變量Ξk的量測為zk,根據(jù)變分貝葉斯理論可知,在變分貝葉斯框架下以變分分布q(Ξk|ψk)作為橋梁可將難以積分的問題轉(zhuǎn)化為ELBO 的優(yōu)化問題.

其中,L(ψk)表示具有單調(diào)遞增特性的變分 ELBO,DKL[q(Ξk|ψk)‖p(Ξk|zk)]表示變分分布q(Ξk|ψk)與后驗(yàn)分布p(Ξk|zk)之間的 KL 散度,其表達(dá)式分別為

其中,q(Ξk|ψk)表示以ψk為參數(shù)的變分分布.

未知變量的估計(jì)實(shí)際是對狀態(tài)后驗(yàn)分布p(Ξk|zk)的近似逼近,當(dāng)非線性狀態(tài)后驗(yàn)分布難以獲取時(shí),變分貝葉斯方法能夠通過構(gòu)建簡單的變分分布實(shí)現(xiàn)對p(Ξk|zk)的近似逼近.變分分布選取需考慮先驗(yàn)分布和后驗(yàn)分布具有相同的函數(shù)形式,一般情況下具有共軛特性的指數(shù)族分布滿足這一條件.本文考慮系統(tǒng)噪聲未知情況,根據(jù)標(biāo)準(zhǔn)卡爾曼濾波實(shí)現(xiàn)結(jié)構(gòu)可知,系統(tǒng)噪聲的統(tǒng)計(jì)特性僅影響狀態(tài)預(yù)測誤差協(xié)方差,因此可直接對狀態(tài)預(yù)測誤差協(xié)方差進(jìn)行先驗(yàn)建模.并且當(dāng)系統(tǒng)噪聲假設(shè)服從均值已知但方差未知的多變量高斯分布假設(shè)時(shí),可選取IW 分布作為其分布方差矩陣的共軛先驗(yàn)分布;此外,針對量測缺失導(dǎo)致量測噪聲出現(xiàn)重尾問題,考慮選取能夠有效表征這一現(xiàn)象的學(xué)生t分布對量測似然函數(shù)進(jìn)行建模.

1.1 IW 分布共軛先驗(yàn)?zāi)Ｐ?/h3>

在貝葉斯濾波框架下,對于系統(tǒng)噪聲高斯分布參數(shù)未知問題的處理,可轉(zhuǎn)化為狀態(tài)預(yù)測誤差協(xié)方差的估計(jì)問題.這里選取IW 分布作為共軛先驗(yàn)分布,以保證后驗(yàn)分布與先驗(yàn)分布具有相同的函數(shù)形式.一個(gè)服從IW 分布的n×n維隨機(jī)對稱正定矩陣A的概率密度函數(shù)可以表示為

其中,t和T分別表示IW 分布的自由度參數(shù)和n×n維逆尺度矩陣,tr(·)表示矩陣的跡,Γn(·)表示n維的伽馬函數(shù).進(jìn)而狀態(tài)預(yù)測誤差協(xié)方差Pk|k-1的先驗(yàn)分布表示為

對于逆威沙特分布,當(dāng)tk≥n+1 時(shí),其變量期望即狀態(tài)預(yù)測誤差協(xié)方差的均值E[Pk|k-1]表示為

令先驗(yàn)分布參數(shù)tk表示為

式 (6)進(jìn)一步轉(zhuǎn)化為

其中,τ表示調(diào)諧參數(shù)[7].根據(jù)式 (5)～(8)計(jì)算得到狀態(tài)預(yù)測協(xié)方差的先驗(yàn)分布及其分布參數(shù).

1.2 學(xué)生t分布量測似然模型

當(dāng)量測噪聲分布由于量測值隨機(jī)異?；蛉笔С尸F(xiàn)重尾特征時(shí),傳統(tǒng)高斯噪聲分布模型難以對噪聲分布特性進(jìn)行有效刻畫.考慮學(xué)生t分布能夠更好地體現(xiàn)這一特征,并且對異常值具有較好的魯棒性.當(dāng)量測噪聲υk滿足均值為0、方差為Rk的分布參數(shù)時(shí),建立參數(shù)化學(xué)生t分布模型,即

其中,v表示學(xué)生t分布的自由度參數(shù),噪聲方差Rk是該分布的尺度矩陣.在此條件下量測似然函數(shù)可表示為

通常學(xué)生t分布的概率密度函數(shù)難以求得封閉解,為了解決這個(gè)問題,引入服從伽馬分布的輔助隨機(jī)變量λk,從而進(jìn)一步將學(xué)生t分布轉(zhuǎn)化為服從參數(shù)化高斯分布的伽馬分布的積分,式 (10)中量測似然函數(shù)可轉(zhuǎn)化為

其中,G(λk;αk,βk)表示輔助變量λk服從伽馬分布,αk和βk分別為形狀參數(shù)和逆尺度參數(shù)[20].綜合上述特點(diǎn),在傳感器量測存在隨機(jī)異常值情況下量測似然函數(shù)表示為如下結(jié)構(gòu)化形式:

并且G(λk;αk,βk)及其均值E[λk]分別表示為

最后,通過對輔助變量λk及其超參數(shù)αk和βk更新確定量測似然的表達(dá)式.

1.3 平均場近似

根據(jù)以上分析建模,式 (2)中的系統(tǒng)隱變量Ξk定義為Ξk=(xk,Pk|k-1,λk).考慮到聯(lián)合后驗(yàn)分布形式的復(fù)雜性以及參數(shù)相互耦合的問題,無法直接求得其解析解.在假設(shè)各未知變量相互獨(dú)立的前提下,通過引入平均場近似策略實(shí)現(xiàn)聯(lián)合變分分布的近似解耦,即

其中,q(xk),q(Pk|k-1)和q(λk)分別表示狀態(tài)xk,狀態(tài)一步預(yù)測協(xié)方差Pk|k-1和輔助變量λk的變分分布,通過平均場理論對待估計(jì)變量聯(lián)合分布進(jìn)行分解,有助于在變分貝葉斯迭代框架結(jié)合坐標(biāo)上升和各類優(yōu)化方法實(shí)現(xiàn)多變量的聯(lián)合優(yōu)化.結(jié)合上述第1.1 節(jié)和第1.2 節(jié)中的共軛先驗(yàn)?zāi)Ｐ秃土繙y似然模型,相應(yīng)的變分分布分別表示為

其中,αk和βk為輔助隨機(jī)變量λk的超參數(shù),同理,tk|k-1和Tk|k-1為隨機(jī)未知變量Pk|k-1的超參數(shù);xk|k和Pk|k為隱變量xk的超參數(shù).在變分貝葉斯優(yōu)化過程中,可結(jié)合平均場理論將聯(lián)合后驗(yàn)的變分分布近似解耦為多個(gè)單變量變分分布,并利用坐標(biāo)上升方法分別對每個(gè)變量進(jìn)行迭代更新.

2 未知變量的變分分布更新

根據(jù)上述構(gòu)建的共軛先驗(yàn)?zāi)Ｐ秃土繙y似然模型,變分置信下界中隱變量Ξk包括狀態(tài)xk,狀態(tài)一步預(yù)測協(xié)方差Pk|k-1和輔助變量λk,Ξk可定義為Ξk=(xk,Pk|k-1,λk).采用變分分布近似上述多未知變量的聯(lián)合后驗(yàn)分布,即

進(jìn)而,相應(yīng)變分置信下界轉(zhuǎn)化為

在變分置信下界最大化約束下,結(jié)合坐標(biāo)上升方法,變量xk,Pk|k-1和λk概率分布的對數(shù)形式的更新分別表示為

結(jié)合建模信息和各變量變分分布的特點(diǎn),聯(lián)合后驗(yàn)分布及其對數(shù)形式分別表示為

其中,cΞk表示迭代過程中與變量相關(guān)的實(shí)常數(shù),Dz和n分別表示向量zk和矩陣Pk|k-1的維數(shù),(xk-表示xk|k-1)的縮寫形式(下文中類似情況不再一一說明).根據(jù)式(26)可知,由于高斯分布、學(xué)生t分布和伽馬分布均屬于指數(shù)分布族而具有簡單的對數(shù)形式,因此便于采用迭代更新的方式計(jì)算其解析解.以下給出λk、Pk|k-1以及xk的迭代估計(jì)的具體實(shí)現(xiàn)原理和步驟.

2.1 λk迭代更新

根據(jù)式 (22)和式 (26),在第i+1 次變分迭代更新過程中,λk的變分分布對數(shù)形式可表示為

其中,cλk表示xk和Pk|k-1積分后相關(guān)的實(shí)常數(shù).根據(jù)貝葉斯無偏估計(jì)理論,,其中和分別表示第i次估計(jì)值及其估計(jì)偏差.同時(shí),結(jié)合局部線性化技術(shù),式 (27)中期望部分可進(jìn)一步改寫為

計(jì)算式 (14)中伽馬分布的對數(shù)表達(dá)式,并結(jié)合式 (27)和式 (28),不難看出,第i+1 次更新后參數(shù)的變分分布中超參數(shù)的更新過程可表示為

根據(jù)式 (15)中伽馬分布的均值表達(dá)式,同時(shí)結(jié)合式 (29)和式 (30),參數(shù)λk的第i+1 次變分迭代更新可表示為

由式 (31)可以清晰地看出,輔助隨機(jī)變量λk的更新與負(fù)相關(guān),并且在迭代過程中的變化與上一次迭代的狀態(tài)估計(jì)值及量測函數(shù)在此估計(jì)值處的局部線性化矩陣有關(guān).根據(jù)表達(dá)式,第i次迭代的估計(jì)值(可理解為第i+1 次迭代估計(jì)的先驗(yàn))越接近真實(shí)狀態(tài),值越小,根據(jù)式 (31)可知?jiǎng)t越大;相反,第i次的估計(jì)狀態(tài)距真實(shí)狀態(tài)較遠(yuǎn),值越大,則越小.其物理意義是,當(dāng)發(fā)生量測隨機(jī)缺失時(shí)導(dǎo)致濾波更新的新息增大,同時(shí)值將隨之增大,從而造成式 (11)中實(shí)際量測噪聲方差增大,這種情況將反饋于濾波過程中狀態(tài)估計(jì)的自適應(yīng)更新.為了直觀地解釋這種現(xiàn)象,結(jié)合式 (11)和式 (13),采用一維高斯分布和伽馬分布說明其相關(guān)關(guān)系.圖1 給出伽馬分布參數(shù)對量測似然的影響示意圖.在圖1 中,假設(shè)高斯分布的均值和方差分別為25 和2/E[λ],并且E[λ]=α/β,伽馬分布參數(shù)α和β的值分別在圖1(a)和圖1(b)中注明.當(dāng)量測未發(fā)生丟失時(shí),β的值較小,高斯分布比較集中(如圖1(a)所示);當(dāng)量測發(fā)生隨機(jī)丟失,β的值增大,高斯分布比較分散(如圖1(b)所示).

圖1 伽馬分布參數(shù)對量測似然函數(shù)的影響示意圖Fig.1 The diagram of the influence of Gamma distribution parameters on likelihood

2.2 Pk|k-1迭代更新

假設(shè)第i次迭代狀態(tài)估計(jì)值已知,根據(jù)式 (23)和式 (26),在第i+1 次變分貝葉斯迭代更新過程中,Pk|k-1的變分分布的對數(shù)形式表示為

因此,根據(jù)式 (6)中Pk|k-1均值表達(dá)式,同時(shí)結(jié)合式 (34)和式 (35),狀態(tài)預(yù)測誤差協(xié)方差的第i+1 次變分迭代更新表示為

根據(jù)卡爾曼濾波實(shí)現(xiàn)原理可知,當(dāng)過程噪聲方差建模不準(zhǔn)確時(shí),狀態(tài)預(yù)測誤差協(xié)方差較大,并且導(dǎo)致量測預(yù)測不精確而產(chǎn)生較大的量測新息zkhk(xk|k-1).同時(shí),由式 (33)可知,的迭代更新與量測新息正相關(guān),進(jìn)而從物理意義上說明狀態(tài)預(yù)測協(xié)方差建模的合理性.

2.3 xk|k迭代更新

根據(jù)式 (24)和式 (26),在第i+1 次迭代更新過程中,系統(tǒng)狀態(tài)的變分分布對數(shù)形式表示為

其中,cxk表示與Pk|k-1和λk積分后相關(guān)性的實(shí)常數(shù).在第i+1 次迭代中,依據(jù)式 (31)和式 (36)已經(jīng)獲取狀態(tài)預(yù)測誤差協(xié)方差Pk|k-1和輔助變量λk的估計(jì)值

由貝葉斯濾波理論可知,狀態(tài)后驗(yàn)分布的近似變分分布應(yīng)具有以下形式

結(jié)合高斯分布表達(dá)式,式 (40)可進(jìn)一步改寫為

綜合式 (37)～(42),目標(biāo)狀態(tài)后驗(yàn)分布的近似變分分布分別是以為期望和協(xié)方差的高斯分布,即

針對非線性系統(tǒng)狀態(tài)估計(jì)問題,現(xiàn)有貝葉斯濾波方法在實(shí)現(xiàn)式 (43)中狀態(tài)估計(jì)及其誤差協(xié)方差的迭代更新過程中,往往難以獲得與真實(shí)后驗(yàn)“緊密”近似的變分分布.雖然在多個(gè)變量坐標(biāo)迭代優(yōu)化過程中能夠彌補(bǔ)部分由于線性化和參數(shù)不精確帶來的誤差,但估計(jì)精度提升效果受限,有時(shí)難以滿足實(shí)際工程需求.因此,考慮在獲得λi+1和估計(jì)值的基礎(chǔ)上,在置信下界最大化條件下結(jié)合自然梯度方法和 Fisher 信息矩陣,推導(dǎo)給出迭代更新解析表達(dá)式的具體實(shí)現(xiàn).

3 基于自然梯度的噪聲自適應(yīng)變分貝葉斯濾波器

3.1 變分置信下界的自然梯度

第i次變分迭代后,在分別獲取參數(shù)λk和Pk|k-1的估計(jì)值的基礎(chǔ)上,式 (2)可轉(zhuǎn)化為僅包含未知變量xk的形式

定義此時(shí)變分參數(shù)ψk=(xk|k,Pk|k),通過最大化式(44)中的置信下界即可獲得狀態(tài)估計(jì)及其誤差協(xié)方差.自然梯度通過信息幾何方法尋找統(tǒng)計(jì)流形空間的目標(biāo)函數(shù)最速下降/上升方向,給出了一種實(shí)現(xiàn)式 (44)置信下界求解的可行性方案[25].下面以置信下界L(ψk)最大化為目標(biāo)函數(shù),同時(shí)結(jié)合Fisher 信息矩陣推導(dǎo)關(guān)于變分參數(shù)ψk的自然梯度.式 (44)的置信下界可以進(jìn)一步表示為

假設(shè)第i次迭代是第i+1 次迭代的先驗(yàn),即式(45中的先驗(yàn)信息p(xk)可以近似為,則變分參數(shù)ψk的最優(yōu)估計(jì)值表示為

通過計(jì)算式 (46)等號(hào)右邊關(guān)于ψk線性化函數(shù),可獲取置信下界的自然梯度以及相應(yīng)的最優(yōu).

定理 1.假設(shè)系統(tǒng)狀態(tài)演化滿足高斯分布特性,并且對數(shù)似然期望Eq[logp(zk|xk)]在的鄰域內(nèi)一階可導(dǎo),同時(shí),在的鄰域內(nèi)二階可導(dǎo),則第i+1 次迭代過程中變分置信下界的自然梯度為

證明.根據(jù)變分參數(shù)的最優(yōu)表達(dá)式 (46)可知,通過對其求關(guān)于ψk的偏導(dǎo)數(shù)即可獲得置信下界最大化條件下的極值點(diǎn).由于對數(shù)似然期望一階可導(dǎo),令,根據(jù)泰勒展開上式對數(shù)似然期望線性化后可以近似表示為

其中,Fψk為具有半正定性 Fisher 信息矩陣,即

此時(shí),式 (46)進(jìn)一步轉(zhuǎn)化為

計(jì)算式(51)關(guān)于ψk的偏導(dǎo),并令其取值為零,可獲得在第i+1 次變分迭代中ELBO 關(guān)于ψk的自然梯度

因此,在置信下界最大化條件下變分參數(shù)ψk的更新表達(dá)式為

需要說明的是,VB 迭代更新等價(jià)于在統(tǒng)計(jì)流形空間沿以 Fisher 信息矩陣作為權(quán)重的自然梯度方向移動(dòng),從而獲得后驗(yàn)分布的最佳近似,并且具有以下性質(zhì):

2)Δψk →0 時(shí),式 (46)中的第2 項(xiàng)KL 散度表達(dá)式轉(zhuǎn)化為式 (49),因此從局部意義上迭代過程中的變分分布與后驗(yàn)分布的近似誤差是正定二次型.

3)自然梯度是Fisher 有效的[25],因此上述優(yōu)化過程是 Fisher 有效的估計(jì)器.特別地,上述ψk的迭代估計(jì)是無偏的.

4)在自然梯度的推導(dǎo)過程中,在無限小的鄰域內(nèi) (即Δψk →0 )給出 Fisher 信息矩陣與 KL 散度之間的近似關(guān)系,并且的計(jì)算依賴于參數(shù)化的變分分布q(xk|ψk),因此,隨著變分迭代的進(jìn)行,Fihser 信息不斷變化.

3.2 基于自然梯度的濾波更新

系統(tǒng)狀態(tài)的變分分布包含兩個(gè)超參數(shù),即當(dāng)前時(shí)刻狀態(tài)估計(jì)xk|k及其誤差協(xié)方差Pk|k,為實(shí)現(xiàn)其迭代更新需分別計(jì)算相應(yīng)的自然梯度.定理 2 及其證明給出了變分置信下界分別關(guān)于xk|k和Pk|k的自然梯度及其推導(dǎo)過程.

定理 2.假設(shè)第i次迭代更新的狀態(tài)估計(jì)及其估計(jì)誤差協(xié)方差分別為,在滿足定理 1 的條件下,對數(shù)似然期望Eq[logp(zk|xk)]分別在xk=的鄰域內(nèi)一階可導(dǎo);同時(shí),KL散度分別在和的鄰域內(nèi)二階可導(dǎo);并且在已經(jīng)分別獲取第i次迭代狀態(tài)預(yù)測誤差協(xié)方差以及輔助變量估計(jì)值的前提下,變分置信下界關(guān)于xk和Pk|k的自然梯度分別表示為

證明.首先推導(dǎo)在定理 2 中假設(shè)條件下 KL 散度關(guān)于xk|k和Pk|k的二階導(dǎo)數(shù).若高斯分布q1～N(ξ1;μ1,C1)和q2～N(ξ2;μ2,C2)具有相同維數(shù)d,則兩分布之間的KL 散度表示為

因此,迭代過程中變分分布之間的 KL 散度表示為

結(jié)合式 (50),關(guān)于xk|k的 Fisher 信息矩陣表示為

關(guān)于Pk|k的 Fisher 信息矩陣的推導(dǎo)如下:

其中,哈達(dá)瑪積?表示取矩陣的對角線元素構(gòu)成單位陣.忽略誤差協(xié)方差矩陣非對角線元素的影響,式(58)可以近似簡化為

由矩陣求導(dǎo)知識(shí)可知,矩陣X逆函數(shù)X→X-1的導(dǎo)數(shù)表示為,其中,m和p分別表示矩陣X的不同行,n和l分別表示其不同列.因此,在迭代優(yōu)化過程中關(guān)于估計(jì)誤差協(xié)方差的 Fisher 信息矩陣表示為

其中,符號(hào)?表示矩陣克羅內(nèi)克積運(yùn)算.

其次,推導(dǎo)Eq[logp(zk|xk)]關(guān)于xk|k和Pk|k的一階偏導(dǎo)數(shù).根據(jù) Bonnet 定理,Eq[logp(zk|xk)]關(guān)于xk|k的梯度可以表示為

結(jié)合服從高斯分布的似然函數(shù)表達(dá)和局部線性化方法,式(61)可轉(zhuǎn)化為

同時(shí),根據(jù) Price 定理可知

由于Eq[zk-hk(xk)]=0,局部線性化處理后,式(63)可表示為

綜合上述過程,在第i+1 次迭代中變分置信下界關(guān)于xk|k和Pk|k的自然梯度分別表示為

結(jié)合式 (53)中變分參數(shù)迭代更新以及式 (66)和式(67)中自然梯度的具體表達(dá)式,第i+1 次迭代時(shí)狀態(tài)估計(jì)值及其誤差協(xié)方差表示為

綜合上述式(68)和式 (69)的構(gòu)建過程,本文給出了一種基于自然梯度的噪聲自適應(yīng)變分貝葉斯濾波算法,為便于理解算法整體實(shí)現(xiàn)流程,下面以偽代碼的形式給出算法具體實(shí)現(xiàn)步驟,詳見算法1.

算法1.基于自然梯度的噪聲自適應(yīng)變分貝葉斯濾波算法

由算法1 中狀態(tài)估計(jì)以及估計(jì)誤差協(xié)方差的更新過程可知,算法1 具有以下特點(diǎn)和優(yōu)勢.

3)自然梯度的方向是目標(biāo)函數(shù)在統(tǒng)計(jì)流形空間中 ELBO 的最速上升方向,即 ELBO 沿基于Fisher 信息矩陣的自然梯度向使其增大的方向移動(dòng),以從信息幾何角度獲得狀態(tài)后驗(yàn)分布的最優(yōu)近似.

為實(shí)現(xiàn)迭代優(yōu)化過程的自適應(yīng),采用以下相對估計(jì)誤差er作為判斷迭代是否終止的條件.

其中,?表示一個(gè)較小的正實(shí)數(shù),其取值大小可根據(jù)實(shí)際應(yīng)用場景的精度需求合理設(shè)定.

4 實(shí)驗(yàn)驗(yàn)證

4.1 仿真算例

為驗(yàn)證本文提出算法的可行性和有效性,采用二維笛卡爾坐標(biāo)系下三雷達(dá)量測的目標(biāo)跟蹤系統(tǒng)進(jìn)I2表示2 維單位陣.初始估計(jì)值和其行500 次 Monte Carlo 仿真實(shí)驗(yàn),目標(biāo)做轉(zhuǎn)彎率未知的勻速轉(zhuǎn)彎運(yùn)動(dòng),并假設(shè)過程噪聲先驗(yàn)信息不準(zhǔn)確且緩慢變化[13],同時(shí)三個(gè)雷達(dá)傳感器均存在量測隨機(jī)缺失現(xiàn)象.本文所提算法記為 VBAKF-NG,其三個(gè)傳感器分布式估計(jì)誤差協(xié)方差融合方法記為DVBAKF-NG.仿真實(shí)驗(yàn)中選取對比算法分別為:結(jié)合本文構(gòu)建的逆威沙特分布的共軛先驗(yàn)分布模型和服從學(xué)生t分布的量測似然模型,分別利用擴(kuò)展卡爾曼濾波器 (Extended Kalman filter,EKF)和無跡卡爾曼濾波 (Unscented Kalman filter,UKF)濾波器實(shí)現(xiàn)狀態(tài)的預(yù)測與更新,其中變分?jǐn)U展卡爾曼自適應(yīng)濾波算法記為 EKF-VB,其對應(yīng)的三個(gè)傳感器分布式估計(jì)誤差協(xié)方差融合方法記為DEKF-VB,變分無跡卡爾曼自適應(yīng)濾波方法記為 UKF-VB,采用傳統(tǒng)自適應(yīng)結(jié)構(gòu)的無跡卡爾曼濾波方法記為 UKF-FN,基于虛擬量測采樣的集合卡爾曼濾波器記為 EnKF,其對應(yīng)的三個(gè)傳感器分布式估計(jì)誤差協(xié)方差融合方法分別記為 DUKF-VB,DUKF-FN 和 DEnKF;同時(shí),將采用真實(shí)量測噪聲和過程噪聲的 EKF 和 UKF 實(shí)現(xiàn)方式分別記為DEKF-TN 和 DUKF-TN.

4.1.1 場景設(shè)置

在k時(shí)刻目標(biāo)狀態(tài)表示為,其中,[xk yk]T和分別表示k時(shí)刻目標(biāo)在水平方向和豎直方向的位置與速度,=-0.052 rad/s表示轉(zhuǎn)彎率.相應(yīng)地,線性化后的狀態(tài)轉(zhuǎn)移矩陣Fk表示為

其中,T=1 表示采樣周期.過程噪聲ωk隨時(shí)間緩慢變化且先驗(yàn)知識(shí)不準(zhǔn)確,假設(shè)系統(tǒng)噪聲協(xié)方差的變化滿足方程

其中,q1=0.3 m2/s2和q2=1×10-4m2/s2分別為噪聲方差的參數(shù),ts表示仿真的總時(shí)長,誤差協(xié)方差分別為x0|0=[2 100 m,0 m/s,1 800 m,30 m/s,-π/50 rad/s]和P0|0=diag{150,10,150,10,0.0001}.直角坐標(biāo)系內(nèi)三個(gè)雷達(dá)的位置分別為(0 m,0 m)、(500 m,0 m)和(0 m,500 m),其量測方程表示為

仿真實(shí)驗(yàn)中,迭代次數(shù)Niter=5,傳感器掃描周期T=1 s,其余仿真實(shí)驗(yàn)參數(shù)取值如表1 所示.

表1 仿真參數(shù)Table 1 Simulation parameters

4.1.2 結(jié)果與分析

圖2 和圖3 分別給出基于傳感器 1 量測的狀態(tài)預(yù)測誤差協(xié)方差對角線上位置分量和速度分量的估計(jì)值,從圖2 和圖3 可以看出,基于本文建模方法的3 種濾波算法均能對預(yù)測誤差協(xié)方差進(jìn)行有效估計(jì),避免了因過程噪聲先驗(yàn)建模不準(zhǔn)對狀態(tài)估計(jì)精度造成的不利影響.同時(shí),VBAKF-NG 在更新階段采用自然梯度優(yōu)化獲得更高精度的估計(jì)和相應(yīng)的誤差協(xié)方差,并作為下一時(shí)刻的先驗(yàn)信息,因而其預(yù)測誤差協(xié)方差最小,并且更有利于在時(shí)序上形成一個(gè)“當(dāng)前時(shí)刻高精度估計(jì)?下一時(shí)刻高精度預(yù)測”的良性循環(huán).

圖2 傳感器1 位置狀態(tài)預(yù)測誤差協(xié)方差的估計(jì)值Fig.2 The expectation of the position state prediction error covariance from Sensor 1

圖3 傳感器1 速度狀態(tài)預(yù)測誤差協(xié)方差的估計(jì)值Fig.3 The expectation of the velocity state prediction error covariance from Sensor 1

圖4～7 分別給出傳感器 2 和傳感器 3 的狀態(tài)預(yù)測誤差協(xié)方差對角線上位置量和速度量的估計(jì)值,可以看出,基于傳感器 2 和傳感器 3 量測的3種算法對狀態(tài)預(yù)測誤差協(xié)方差估計(jì)的曲線與傳感器1 相似,即VBAKF-NG 預(yù)測誤差協(xié)方差位置分量和速度分量的估計(jì)值曲線始終處于最優(yōu)的位置且更加平滑,體現(xiàn)出相對于 EKF-VB 和 UKF-VB 能夠獲取更高的估計(jì)精度和魯棒性.

圖4 傳感器2 位置狀態(tài)預(yù)測誤差協(xié)方差的估計(jì)值Fig.4 The expectation of the position state prediction error covariance from Sensor 2

圖5 傳感器2 速度狀態(tài)預(yù)測誤差協(xié)方差的估計(jì)值Fig.5 The expectation of the velocity state prediction error covariance from Sensor 2

圖6 傳感器3 位置狀態(tài)預(yù)測誤差協(xié)方差的估計(jì)值Fig.6 The expectation of the position state prediction error covariance from Sensor 3

圖7 傳感器3 速度狀態(tài)預(yù)測誤差協(xié)方差的估計(jì)值Fig.7 The expectation of the velocity state prediction error covariance from Sensor 3

圖8給出傳感器 1 在徑向距和方位角上的量測噪聲方差.圖8 中黑色跳變表示該傳感器此時(shí)出現(xiàn)量測丟失情況,從圖8 中可以看出,UKF-VB 和 VBAKF-NG 能夠準(zhǔn)確地識(shí)別出量測噪聲方差的跳變,以便及時(shí)調(diào)整分布參數(shù)更新.同時(shí),從圖8 中可以清晰地看出,VBAKF-NG 對噪聲方差跳變的識(shí)別精度更高.圖9 和圖10 分別給出傳感器 2 和傳感器 3 在徑向距和方位角上的量測噪聲方差.結(jié)果同樣表明與 UKF-VB 和 EKF-VB 相比,VBAKFNG 對噪聲方差跳變具有更高的識(shí)別精度.實(shí)際上,圖2～7 中 VBAKF-NG 之所以在濾波精度和魯棒性上優(yōu)于 EKF-VB 和UKF-VB,其原因正是由于VBAKF-NG 可以有效辨識(shí)量測隨機(jī)缺失現(xiàn)象產(chǎn)生,從而使得算法具有較好的容錯(cuò)特性.

圖8 傳感器1 量測噪聲方差的估計(jì)值Fig.8 The expectation of the measurement noise variance from Sensor 1

圖9 傳感器2 量測噪聲方差的估計(jì)值Fig.9 The expectation of the measurement noise variance from Sensor 2

圖10 傳感器3 量測噪聲方差的估計(jì)值Fig.10 The expectation of the measurement noise variance from Sensor 3

圖11和圖12 分別給出上述6 種方法狀態(tài)估計(jì)的徑向距均方根誤差 (Root mean squared error,RMSE)和徑向速度均方根誤差的對比,從圖11 和圖12 中可以看出,采用IW 分布和學(xué)生t分布分別建模預(yù)測誤差協(xié)方差和量測噪聲的 DEKF-VB,DUKF-VB 和 DVBAKF-NG 三種方法均優(yōu)于傳統(tǒng)自適應(yīng)濾波算法 DUKF-FN 和集合卡爾曼濾波算法 DEnKF.同時(shí),由于 DVBAKF-NG 在更新過程中采用自然梯度的方法更新目標(biāo)狀態(tài),與 DEKFVB,DUKF-VB 相比,取得了相對更優(yōu)的估計(jì)精度.需要說明的是,DUKF-TN 在濾波過程中采用真實(shí)的過程噪聲,因此具有最高的精度.

圖11 徑向距估計(jì)RMSE 的對比Fig.11 The RMSE comparison of radical range estimation

圖12 徑向速度估計(jì)RMSE 的對比Fig.12 The RMSE comparison of radical velocity estimation

4.2 毫米波雷達(dá)數(shù)據(jù)驗(yàn)證

實(shí)驗(yàn)場景為采用 ARS408-21 毫米波雷達(dá)對地面沿直線運(yùn)動(dòng)的行人目標(biāo)進(jìn)行跟蹤.在量測過程中考慮量測噪聲非高斯特點(diǎn),行人分別在 17 時(shí)刻和27 時(shí)刻被遮擋,導(dǎo)致出現(xiàn)兩次量測丟失.同時(shí)考慮系統(tǒng)噪聲方差的未知時(shí)變情況,行人的速度大小不斷變化.在不考慮雜波和數(shù)據(jù)關(guān)聯(lián)的情況下,假設(shè)一個(gè)時(shí)刻返回一個(gè)量測,量測數(shù)據(jù)和濾波結(jié)果如圖13 所示.

圖13 毫米波雷達(dá)地面行人跟蹤實(shí)驗(yàn)Fig.13 Tracking a pedestrian by using a millimeter wave radar

需要說明的是,由于真實(shí)系統(tǒng)噪聲未知,在此實(shí)驗(yàn)中采用真實(shí)系統(tǒng)噪聲方差的 UKF-TN 算法未參與對比.從圖13 中可以看出,由于系統(tǒng)噪聲未知,在濾波初始階段,EnKF 和 UKF-FN 算法的估計(jì)結(jié)果與量測偏離較大,采用IW 分布建模的 EKFVB、UKF-VB 和 VBAKF-NG 與量測偏離較小.當(dāng)發(fā)生量測缺失時(shí),對比算法均存在不同程度的發(fā)散現(xiàn)象,VBAKF-NG 能夠較好地抑制濾波發(fā)散.因此,毫米波雷達(dá)地面行人跟蹤實(shí)驗(yàn)驗(yàn)證了本文所提 VBAKF-NG 算法的優(yōu)越性.

5 結(jié)束語

針對過程噪聲時(shí)變且先驗(yàn)信息建模不準(zhǔn)確和量測噪聲非高斯問題,首先,引入?yún)?shù)化的IW 分布作為狀態(tài)預(yù)測誤差協(xié)方差矩陣的共軛分布建立先驗(yàn)信息模型,同時(shí)構(gòu)建參數(shù)化學(xué)生t分布模型以近似由于量測隨機(jī)丟失造成的非高斯量測似然函數(shù);繼而,結(jié)合平均場理論近似解耦狀態(tài)和多參數(shù)的聯(lián)合概率密度函數(shù)以獲取相應(yīng)的邊緣變分分布,并采用坐標(biāo)上升的迭代方法更新狀態(tài)和各參數(shù)的邊緣變分分布;在此基礎(chǔ)上,結(jié)合 Fisher 信息矩陣推導(dǎo)給出置信下界關(guān)于狀態(tài)估計(jì)及其誤差協(xié)方差的自然梯度,進(jìn)而給出一種基于自然梯度的噪聲自適應(yīng)變分貝葉斯濾波算法,實(shí)現(xiàn)過程噪聲時(shí)變和由于量測隨機(jī)缺失導(dǎo)致的量測噪聲非高斯條件下非線性系統(tǒng)的高精度估計(jì).本文算法處理的對象局限在單目標(biāo)跟蹤系統(tǒng),如何將其拓展到機(jī)動(dòng)多目標(biāo)跟蹤系統(tǒng),在變分貝葉斯框架下,結(jié)合自然梯度策略有效地解決模型自適應(yīng)選取以及多目標(biāo)跟蹤中數(shù)據(jù)關(guān)聯(lián)問題是下一步將開展的工作.