基于自然梯度的噪聲自適應變分貝葉斯濾波算法

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

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

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

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

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

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

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

1 預備知識

假設動態系統隱變量Ξk的量測為zk,根據變分貝葉斯理論可知,在變分貝葉斯框架下以變分分布q(Ξk|ψk)作為橋梁可將難以積分的問題轉化為ELBO 的優化問題.

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

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

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

1.1 IW 分布共軛先驗模型

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

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

對于逆威沙特分布,當tk≥n+1 時,其變量期望即狀態預測誤差協方差的均值E[Pk|k-1]表示為

令先驗分布參數tk表示為

式 (6)進一步轉化為

其中,τ表示調諧參數[7].根據式 (5)～(8)計算得到狀態預測協方差的先驗分布及其分布參數.

1.2 學生t分布量測似然模型

當量測噪聲分布由于量測值隨機異常或缺失呈現重尾特征時,傳統高斯噪聲分布模型難以對噪聲分布特性進行有效刻畫.考慮學生t分布能夠更好地體現這一特征,并且對異常值具有較好的魯棒性.當量測噪聲υk滿足均值為0、方差為Rk的分布參數時,建立參數化學生t分布模型,即

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

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

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

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

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

1.3 平均場近似

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

其中,q(xk),q(Pk|k-1)和q(λk)分別表示狀態xk,狀態一步預測協方差Pk|k-1和輔助變量λk的變分分布,通過平均場理論對待估計變量聯合分布進行分解,有助于在變分貝葉斯迭代框架結合坐標上升和各類優化方法實現多變量的聯合優化.結合上述第1.1 節和第1.2 節中的共軛先驗模型和量測似然模型,相應的變分分布分別表示為

其中,αk和βk為輔助隨機變量λk的超參數,同理,tk|k-1和Tk|k-1為隨機未知變量Pk|k-1的超參數;xk|k和Pk|k為隱變量xk的超參數.在變分貝葉斯優化過程中,可結合平均場理論將聯合后驗的變分分布近似解耦為多個單變量變分分布,并利用坐標上升方法分別對每個變量進行迭代更新.

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

根據上述構建的共軛先驗模型和量測似然模型,變分置信下界中隱變量Ξk包括狀態xk,狀態一步預測協方差Pk|k-1和輔助變量λk,Ξk可定義為Ξk=(xk,Pk|k-1,λk).采用變分分布近似上述多未知變量的聯合后驗分布,即

進而,相應變分置信下界轉化為

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

結合建模信息和各變量變分分布的特點,聯合后驗分布及其對數形式分別表示為

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

2.1 λk迭代更新

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

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

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

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

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

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

2.2 Pk|k-1迭代更新

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

因此,根據式 (6)中Pk|k-1均值表達式,同時結合式 (34)和式 (35),狀態預測誤差協方差的第i+1 次變分迭代更新表示為

根據卡爾曼濾波實現原理可知,當過程噪聲方差建模不準確時,狀態預測誤差協方差較大,并且導致量測預測不精確而產生較大的量測新息zkhk(xk|k-1).同時,由式 (33)可知,的迭代更新與量測新息正相關,進而從物理意義上說明狀態預測協方差建模的合理性.

2.3 xk|k迭代更新

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

其中,cxk表示與Pk|k-1和λk積分后相關性的實常數.在第i+1 次迭代中,依據式 (31)和式 (36)已經獲取狀態預測誤差協方差Pk|k-1和輔助變量λk的估計值

由貝葉斯濾波理論可知,狀態后驗分布的近似變分分布應具有以下形式

結合高斯分布表達式,式 (40)可進一步改寫為

綜合式 (37)～(42),目標狀態后驗分布的近似變分分布分別是以為期望和協方差的高斯分布,即

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

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

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

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

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

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

通過計算式 (46)等號右邊關于ψk線性化函數,可獲取置信下界的自然梯度以及相應的最優.

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

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

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

此時,式 (46)進一步轉化為

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

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

需要說明的是,VB 迭代更新等價于在統計流形空間沿以 Fisher 信息矩陣作為權重的自然梯度方向移動,從而獲得后驗分布的最佳近似,并且具有以下性質:

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

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

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

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

系統狀態的變分分布包含兩個超參數,即當前時刻狀態估計xk|k及其誤差協方差Pk|k,為實現其迭代更新需分別計算相應的自然梯度.定理 2 及其證明給出了變分置信下界分別關于xk|k和Pk|k的自然梯度及其推導過程.

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

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

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

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

關于Pk|k的 Fisher 信息矩陣的推導如下:

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

由矩陣求導知識可知,矩陣X逆函數X→X-1的導數表示為,其中,m和p分別表示矩陣X的不同行,n和l分別表示其不同列.因此,在迭代優化過程中關于估計誤差協方差的 Fisher 信息矩陣表示為

其中,符號?表示矩陣克羅內克積運算.

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

結合服從高斯分布的似然函數表達和局部線性化方法,式(61)可轉化為

同時,根據 Price 定理可知

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

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

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

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

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

由算法1 中狀態估計以及估計誤差協方差的更新過程可知,算法1 具有以下特點和優勢.

3)自然梯度的方向是目標函數在統計流形空間中 ELBO 的最速上升方向,即 ELBO 沿基于Fisher 信息矩陣的自然梯度向使其增大的方向移動,以從信息幾何角度獲得狀態后驗分布的最優近似.

為實現迭代優化過程的自適應,采用以下相對估計誤差er作為判斷迭代是否終止的條件.

其中,?表示一個較小的正實數,其取值大小可根據實際應用場景的精度需求合理設定.

4 實驗驗證

4.1 仿真算例

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

4.1.1 場景設置

在k時刻目標狀態表示為,其中,[xk yk]T和分別表示k時刻目標在水平方向和豎直方向的位置與速度,=-0.052 rad/s表示轉彎率.相應地,線性化后的狀態轉移矩陣Fk表示為

其中,T=1 表示采樣周期.過程噪聲ωk隨時間緩慢變化且先驗知識不準確,假設系統噪聲協方差的變化滿足方程

其中,q1=0.3 m2/s2和q2=1×10-4m2/s2分別為噪聲方差的參數,ts表示仿真的總時長,誤差協方差分別為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}.直角坐標系內三個雷達的位置分別為(0 m,0 m)、(500 m,0 m)和(0 m,500 m),其量測方程表示為

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

表1 仿真參數Table 1 Simulation parameters

4.1.2 結果與分析

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

圖2 傳感器1 位置狀態預測誤差協方差的估計值Fig.2 The expectation of the position state prediction error covariance from Sensor 1

圖3 傳感器1 速度狀態預測誤差協方差的估計值Fig.3 The expectation of the velocity state prediction error covariance from Sensor 1

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

圖4 傳感器2 位置狀態預測誤差協方差的估計值Fig.4 The expectation of the position state prediction error covariance from Sensor 2

圖5 傳感器2 速度狀態預測誤差協方差的估計值Fig.5 The expectation of the velocity state prediction error covariance from Sensor 2

圖6 傳感器3 位置狀態預測誤差協方差的估計值Fig.6 The expectation of the position state prediction error covariance from Sensor 3

圖7 傳感器3 速度狀態預測誤差協方差的估計值Fig.7 The expectation of the velocity state prediction error covariance from Sensor 3

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

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

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

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

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

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

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

4.2 毫米波雷達數據驗證

實驗場景為采用 ARS408-21 毫米波雷達對地面沿直線運動的行人目標進行跟蹤.在量測過程中考慮量測噪聲非高斯特點,行人分別在 17 時刻和27 時刻被遮擋,導致出現兩次量測丟失.同時考慮系統噪聲方差的未知時變情況,行人的速度大小不斷變化.在不考慮雜波和數據關聯的情況下,假設一個時刻返回一個量測,量測數據和濾波結果如圖13 所示.

圖13 毫米波雷達地面行人跟蹤實驗Fig.13 Tracking a pedestrian by using a millimeter wave radar

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

5 結束語

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