貝葉斯序貫重要性積分濾波器

張宏偉，張小虎，曹 勇

（1. 中山大學(xué)航空航天學(xué)院，廣東廣州 510725；2. 中國科學(xué)院空間精密測量技術(shù)重點實驗室，陜西西安 710119；3. 北京東方計量測試研究所，北京 100029）

1 引言

濾波是指從含有噪聲的傳感測量中估計目標(biāo)的未知狀態(tài)或參數(shù)，非線性濾波是解決該類不確定信號處理問題的核心技術(shù)，已被廣泛應(yīng)用于工程、統(tǒng)計和數(shù)學(xué)領(lǐng)域［1～8］. 貝葉斯推理為濾波算法提供了一種數(shù)學(xué)機制，然而在實際應(yīng)用中，由于非線性動態(tài)系統(tǒng)函數(shù)的復(fù)雜性，難以對方程的積分運算進行閉式求解，通常需利用數(shù)值方法近似解算，主要包括數(shù)值積分和蒙特卡羅兩大類.第一類數(shù)值積分方法基于確定性法則. 高斯匹配法［1］近似目標(biāo)分布簡單直接，然而數(shù)值精度較為一般. 擴展卡爾曼濾波（Extended Kalman Filter，EKF）算法［2］基于泰勒級數(shù)展開，在工程上最易實現(xiàn)且應(yīng)用最為廣泛，然而在復(fù)雜傳感測量環(huán)境中會引入截斷誤差和基點誤差. 無跡卡爾曼濾波（Unscented Kalman Filter，UKF）算法［3］選擇固定數(shù)量的sigma 點進行無跡變換以近似目標(biāo)原始分布，在滿足協(xié)方差矩陣為正定的前提條件下，可達到3階數(shù)值精度. 當(dāng)狀態(tài)空間的積分維度適中且滿足假定概率模型時，Gauss-Hermite積分算法能夠最大程度地減小近似誤差［4］. 通常，該類改進卡爾曼濾波（Kalman Filter，KF）算法的應(yīng)用局限于高斯噪聲擾動，因此制約了確定性數(shù)值積分方法在強非線性場景中的應(yīng)用［5］. 第二類蒙特卡羅方法基于隨機采樣的思想. 該方法能夠有效處理非線性非高斯濾波問題，無需對非線性函數(shù)進行線性化處理. 馬爾可夫鏈蒙特卡羅（Monte Carlo Markov Chain，MCMC）算法［6］通過模擬與目標(biāo)分布一致的馬爾可夫鏈生成有效的蒙特卡羅樣本，然而在樣本較少的情況下難以維持狀態(tài)為穩(wěn)態(tài)分布，增大數(shù)據(jù)采集量將導(dǎo)致計算復(fù)雜度和時間成本急劇增大. 與此不同，重要性采樣（Importance Sampling，IS）算法的優(yōu)勢在于可通過調(diào)整樣本權(quán)重以校正目標(biāo)重要性函數(shù)和目標(biāo)真實分布之間的偏差，因此，構(gòu)建能夠覆蓋有效似然區(qū)域的重要性函數(shù)至為關(guān)鍵［7］.Doucet等人［8］提出序貫蒙特卡羅采樣也稱為粒子濾波（Particle Filter，PF）方法，在貝葉斯濾波框架下推導(dǎo)了離散時間動態(tài)模型的重要性采樣法則. 理論上，該類蒙特卡羅方法的采樣誤差具有維數(shù)不變性，數(shù)學(xué)上的大數(shù)中心極限定理保證其計算收斂性. 然而由維數(shù)災(zāi)難造成的樣本急劇衰減和權(quán)值退化仍是亟待解決的技術(shù)難題.

從貝葉斯概率論的數(shù)學(xué)模型分析，復(fù)雜非線性高維狀態(tài)空間中的參數(shù)估計屬于模糊測度，難以精確解析［9，10］. 其根本原因不在于外在的統(tǒng)計變化量，而是源于復(fù)雜非線性動態(tài)系統(tǒng)的內(nèi)在模糊性，主要包括目標(biāo)模型屬性歧義和非線性不穩(wěn)定傳輸［11］. 該模糊性制約了許多單個非線性濾波算法的應(yīng)用，許多研究者提出了結(jié)合上述兩類數(shù)值方法的改進算法. 如通過EKF，UKF 及Gauss-Hermite積分等非線性濾波算法生成重要性函數(shù)，即擴展卡爾曼粒子濾波（Extended Kalman Particle Filter，EPF）［12］，無跡卡爾曼粒子濾波（Unscented Kalman Particle Filter，UPF）［13］，高斯-施密特粒子濾波（Gauss-Hermite Particle Filter，GHPF）［14］.Garcia-Fernandez等人［15］證明，當(dāng)前時刻觀測信息較為準(zhǔn)確時，改進第一、二階矩?zé)o法直接調(diào)制重要性函數(shù)和目標(biāo)分布之間的不匹配問題.

另外，多模型估計器廣泛應(yīng)用于非線性動態(tài)系統(tǒng)的狀態(tài)空間的細化［11］. 其中，交互多模型（Interacting Multiple Model，IMM）算法被成功應(yīng)用于目標(biāo)機動跟蹤，該算法通過融合有限個單濾波器的濾波輸出估計機動目標(biāo)狀態(tài). 事實上，無法通過有限個目標(biāo)狀態(tài)模型表征全域的目標(biāo)機動行為，同時不必要的模型競爭會加大計算復(fù)雜度［16］. Hostettler 等人［17］根據(jù)Rao-Blackwell理論對狀態(tài)空間降維，分別采用線性高斯逼近和粒子濾波方法估計系統(tǒng)的狀態(tài)模型參數(shù)和模型狀態(tài)演化. 王琪等人［18］采用Rao-Blackwell 粒子濾波方法，有效緩解了相控陣導(dǎo)引頭指向誤差斜率給導(dǎo)彈制導(dǎo)系統(tǒng)帶來的寄生回路振蕩問題. 該類模型驅(qū)動的濾波方法依賴具體問題，對于復(fù)雜非線性動態(tài)系統(tǒng)的狀態(tài)估計，難以平衡模型切換的準(zhǔn)確性和計算復(fù)雜度.

從香農(nóng)信息論的角度分析，充分利用系統(tǒng)的約束信息有助于實現(xiàn)狀態(tài)空間降維，提高模型預(yù)測的可靠性和狀態(tài)濾波的準(zhǔn)確性［10，19，20］. 近年來，約束非線性濾波算法發(fā)展迅速，其優(yōu)勢在于可根據(jù)應(yīng)用場景選取輔助條件. 現(xiàn)有公開文獻已提出截斷觀測驅(qū)動的約束UKF［19］、融入狀態(tài)約束的PF［20］、加載深度空間語義特征的PF［21］及基于注意力學(xué)習(xí)的正則化相關(guān)濾波［22］等非線性動態(tài)狀態(tài)估計和目標(biāo)跟蹤算法. 該類算法依賴隱馬爾可夫鏈的動力學(xué)模型［23］，在強非線性環(huán)境中統(tǒng)計迭代易陷入局部最優(yōu)，缺乏較強的泛化能力.

本質(zhì)上，觀測量和狀態(tài)變量之間存在因果映射［11，24］. 然而實際上，由于過程噪聲、觀測噪聲和目標(biāo)機動運動等不確定因素，在非線性濾波中存在觀測似然、重要性函數(shù)和真實目標(biāo)后驗分布失配的現(xiàn)象，本文作者在專利中分析了產(chǎn)生該現(xiàn)象的本質(zhì)原因［25］. 通常，確定的等式或不等式等硬約束條件會破壞常規(guī)的高斯分布假設(shè)，因此難以對嚴(yán)格的數(shù)值邊界條件進行閉式解析. 軟約束技術(shù)能夠有效表征非線性動態(tài)系統(tǒng)中不確定變量之間的模糊性［10，19，20］. 文獻［26］將統(tǒng)計熵的概念推廣到不確定信息處理的度量，理論上可通過相關(guān)信息熵測度多模態(tài)似然對狀態(tài)估計的影響. 基于此，本文推導(dǎo)了貝葉斯序貫重要性積分濾波（Bayesian Sequential Importance Quadrature Filter，SIQF）算法. 構(gòu)建系統(tǒng)軟空時約束下的有界似然，以消減由非線性動態(tài)系統(tǒng)的不可預(yù)測性造成的似然分布和目標(biāo)分布之間的偏差. 通過并行Gauss-Hermite 積分將修正先驗和反饋補償融入狀態(tài)演化過程，引入最大相關(guān)熵調(diào)制重要性函數(shù)和真實目標(biāo)分布的一致性，從而增強重要性采樣的多樣性和預(yù)測更新的容錯性.

2 問題描述

2.1 系統(tǒng)模型

本文利用空時約束信息構(gòu)建多模態(tài)似然以優(yōu)化系統(tǒng)狀態(tài)估計，記為模型集M={M1，…，MK} ，K≥2. 該集合包含K個相互獨立的截斷先驗和后驗反饋等元素，將其融入貝葉斯遞推過程. 由此，離散非線性動態(tài)系統(tǒng)模型可表示為

其中，k∈N 是離散時刻；Xk∈RnX和Zk∈RnZ分別是狀態(tài)和觀測矢量，RnX和RnZ分別是nX維和nZ維空間模型；f和h分別是非線性狀態(tài)方程和觀測方程；υk和ek分別是隨機過程噪聲和觀測噪聲，其概率密度函數(shù)分別記為pυ(υk)和pe(ek)，假設(shè)為相互獨立的0 均值白噪聲，協(xié)方差分別為Συ和Σe；M i k表示第i個模型.

狀態(tài)變量的初始先驗記為p0(Xk)，轉(zhuǎn)移概率為隱馬爾可夫序列，似然函數(shù)關(guān)聯(lián)系統(tǒng)狀態(tài)Xk和觀測Zk之間的所有信息. 根據(jù)貝葉斯全概率公式［1，8］，由觀測量推斷非線性動態(tài)系統(tǒng)中未知狀態(tài)的后驗概率密度為

其中，分母中的歸一化常數(shù)p(Z1：k)的積分運算難以解算，因此多區(qū)域似然的后驗分布近似為多高斯分布；和分別是第i個模型下狀態(tài)分布的均值和方差，分別定義為

在實際應(yīng)用中，目標(biāo)狀態(tài)通常在直角坐標(biāo)系推演，而觀測數(shù)據(jù)在球面坐標(biāo)系或者像坐標(biāo)系傳輸. 由于類KF 算法受限于高斯假設(shè)，因此無法通過提高矩的精度提升整體濾波性能［2，5，11］. 由此可知，在非線性較強的濾波環(huán)境中，式（4）和式（5）難以解算. 從算法機理角度出發(fā)，重要性采樣方法能夠有效處理非線性非高斯問題. 然而，對于似然函數(shù)、重要性函數(shù)偏離目標(biāo)分布的情況，該算法會發(fā)生樣本急劇衰減和權(quán)值退化，造成濾波發(fā)散. 針對該問題，本文構(gòu)建有界似然并推導(dǎo)能夠覆蓋多似然區(qū)域的多模重要性積分算法.

2.2 提出算法框架

非線性濾波近似的主要困難在于更新步驟產(chǎn)生的誤差［5，11，16］，較為常見的例子是若當(dāng)前測量包含更多有效信息，而過多依賴先驗近似目標(biāo)后驗分布會造成較大偏差. 在實際的物理動態(tài)系統(tǒng)模型中，傳感觀測系統(tǒng)提供的測量信息是有限的，因此觀測似然應(yīng)為有界分布. 基于此，通過截斷觀測噪聲概率密度pe(·)構(gòu)建可行域Λ(Zk)，即

其中，h-1( ·)是式（2）中觀測方程h的反函數(shù).

對于存在多域似然的濾波問題，目標(biāo)狀態(tài)后驗應(yīng)為多模分布. 因此將其近似為加權(quán)混合高斯分布，即

其中，和分別是第i個似然高斯分布的權(quán)重、均值和協(xié)方差.

綜上所述，本文構(gòu)建了系統(tǒng)空時軟約束模型，將有界似然和后驗反饋等先驗信息融入狀態(tài)演化過程，以提高整體濾波性能. 考慮動態(tài)系統(tǒng)的隨機干擾和約束函數(shù)的模糊性，為了避免陷入局部最優(yōu)或鞍點值，采用全局牛頓內(nèi)點法遍歷可行域迭代求解帕累托解，以近似當(dāng)前有效觀測驅(qū)動的修正先驗. 為了處理非線性非高斯估計問題，通過并行Gauss-Hermite 遞推狀態(tài)以生成混合重要性分布，綜合改善重要性采樣的多樣性和預(yù)測協(xié)方差的容錯性. 圖1為算法流程.

圖1 算法流程框架

3 貝葉斯序貫重要性積分濾波算法的實現(xiàn)

3.1 近似有界似然先驗

根據(jù)式（6）定義的可行域Λ(Zk)可構(gòu)建當(dāng)前觀測驅(qū)動的有界似然. 同時，在貝葉斯全概率公式中引入輔助變量集a，目標(biāo)狀態(tài)的后驗概率密度函數(shù)修正為

其中，plik(·)和p1(·)分別是似然函數(shù)和修正先驗的概率密度 函數(shù)；?為 常數(shù)；是指示函數(shù)，若取為1，否則為0.

在近似修正先驗的數(shù)值規(guī)劃過程中，需分解目標(biāo)狀態(tài)的標(biāo)量和矢量，選擇標(biāo)量分量構(gòu)成凸的目標(biāo)函數(shù).實際應(yīng)用中，可通過最小二乘和多視圖對極幾何約束等測量手段將目標(biāo)初始位置X0定位在可行域Λ(Zk)內(nèi). 從幾何上看，可行域的中心處觀測誤差應(yīng)為最小.因此，構(gòu)建目標(biāo)函數(shù)為

在可行域內(nèi)，通過backtracking 回溯算法求解真實目標(biāo)狀態(tài)分布的均值，即

其中，s和dk分別是迭代步長和搜索方向，該數(shù)值規(guī)劃優(yōu)化方法的詳細步驟請參考本文作者之前的工作［19］. 根據(jù)該統(tǒng)計中心計算可行域的第二階矩，即方差.

綜上所述，在第k時刻，可行域的中心或修正先驗概率密度函數(shù)記為其中，下標(biāo)“1”用以區(qū)分原始先驗概率. 相應(yīng)地，原始先驗概率分布記為這里請注意，可通過擴維方法獲取目標(biāo)狀態(tài)矢量集［3］. 為了覆蓋多似然區(qū)域，通過加權(quán)p1( ·)和p0( ·)修正目標(biāo)狀態(tài)先驗為

其中，α是權(quán)重因子，用來衡量不同區(qū)域似然先驗對目標(biāo)狀態(tài)估計的影響. 其求解過程將在下一節(jié)詳細描述.

3.2 并行Gauss-Hermite 生成重要性分布

根據(jù)Gauss-Hermite 積分法則［5，6］，可通過有限個一維高斯分布迭代近似多維高斯分布的積分，即

其中，I表示單位矩陣是積分點及其權(quán)值；l是Hermite 多項式的特征根個數(shù)；是方差的Cholesky解，即

3.2.1 預(yù)測

在第k-1 時刻，從修正先驗分布選取2l-1個nX維積分點為

狀態(tài)預(yù)測的均值和協(xié)方差為

測量預(yù)測均值為

3.2.2 更新

計算測量的預(yù)測協(xié)方差和測量的互協(xié)方差為

計算卡爾曼濾波增益為

濾波輸出均值為

濾波輸出協(xié)方差為

至此，采用最大相關(guān)信息熵度量本文2個先驗分布之間的相似性. 構(gòu)建最大似然目標(biāo)函數(shù)，以求解權(quán)重因子α，即

其中，gi( ·)是第i個似然分布的核函數(shù)，定義為

其中，exp( ·)是以自然常數(shù)e為底的指數(shù)函數(shù).

3.3 重要性采樣更新

對應(yīng)式（13）的多模先驗分布，根據(jù)多模重要性分布采樣重要性樣本，即

計算樣本權(quán)重為

歸一化重要性樣本的權(quán)值為

根據(jù)蒙特卡羅法則，通過狄拉克函數(shù)δ( ·)加權(quán)近似目標(biāo)狀態(tài)的后驗概率密度函數(shù)為

在已知目標(biāo)狀態(tài)后驗分布的情況下，加權(quán)重要性樣本估計真實目標(biāo)狀態(tài)分布的均值和方差為

4 仿真實驗與分析

本節(jié)研究以下2 個例子. 第一個例子是一維單變量非平穩(wěn)增長模型（Univariate Nonstationary Growth Model，UNGM），將SIQF 算法與文獻［13］和文獻［19］中的擴展卡爾曼濾波（EKF）、無跡卡爾曼濾波（UKF）、傳統(tǒng)粒子濾波（Generical Particle Filter，GPF）、擴展卡爾曼粒子濾波和無跡卡爾曼濾波算法比較；第二個例子采用我國某通用民航公司實測民航航跡作為目標(biāo)數(shù)據(jù)，觀測傳感器為雷達. 由于測量環(huán)境干擾等不確定因素，部分觀測量之間時間間隔達數(shù)十秒，故稱之為稀疏觀測. 在該稀疏觀測場景中，EKF 和UKF 跟蹤誤差過大甚至發(fā)散［19］. 因此，將SIQF 算法與文獻［11］、文獻［16］和文獻［18］中的交互多模型擴展卡爾曼濾波（Interacting Multiple Model Extended Kalman Filter，IMMEKF）、交互多模型無跡卡爾曼濾波（Interacting Multiple Model Unscented Kalman Filter，IMMUKF）以及基于Rao-Blackwell 的多模型粒子濾波（Multiple Model Rao-Blackwell Particle Filter，MMRBPF）算法比較.

4.1 一維單變量非平穩(wěn)增長模型

一維單變量非平穩(wěn)增長模型的狀態(tài)方程和觀測方程模型分別為

其中，系數(shù)為已知常數(shù)，a=0.5，b=2.5，c=8，φ1=0.2，φ2=0.5；狀態(tài)初值X0在[0，1]范圍內(nèi)平均取值；υk是服從參數(shù)為[2，3]的伽馬分布的過程噪聲；ek是均值為0且方差為0.5的觀測噪聲.

構(gòu)建合適的目標(biāo)重要性函數(shù)是重要性采樣算法的關(guān)鍵步驟，它與實際目標(biāo)分布的匹配程度直接影響整體濾波性能和樣本數(shù)目及計算復(fù)雜度.

圖2～5 分別是采樣10，200，300，500 個樣本的結(jié)果. 從濾波誤差趨勢的定性比較分析可知：相比EKF和UKF 算法，GPF 顯示出處理非高斯噪聲的優(yōu)勢.EPF 和UPF 算法分別由EKF 和UKF 算法的非線性估計生成重要性函數(shù).Gauss-Hermite 方法估計精度高于EKF，且狀態(tài)方差有較重的尾部，因此SIQF濾波性能優(yōu)于UPF.

圖2 10個樣本

圖3 200個樣本

圖4 300個樣本

圖5是樣本增多到500 個，單獨比較基于蒙特卡羅采樣的數(shù)值方法的濾波性能. 隨著時間推進，在狀態(tài)演化過程中，EPF，UPF 和GPF 濾波算法的估計誤差在第30 s 后明顯增大，出現(xiàn)這個現(xiàn)象的主要原因是重要性函數(shù)與實際目標(biāo)的真實分布失配引起樣本多樣性衰減和權(quán)值退化. 相比UPF 濾波算法，SIQF 算法在無需犧牲計算復(fù)雜度的情況下將平均誤差減小了63%. 這主要是因為SIQF 融合截斷先驗和后驗反饋構(gòu)建覆蓋多似然的目標(biāo)重要性函數(shù)，從根本上增強了重要性樣本的多樣性和預(yù)測準(zhǔn)確性.

圖5 500個樣本

4.2 稀疏觀測機動目標(biāo)跟蹤

圖6 是我國某通用民航公司廣播式自動相關(guān)監(jiān)視（Automatic Dependent Surveillance Broadcast，ADS-B）系統(tǒng)監(jiān)測的某民航飛機航跡信息. 其中，黑色曲線是由全球定位系統(tǒng)（Golobal Positional System，GPS）接收的航跡，紅色曲線是監(jiān)測雷達系統(tǒng)的觀測數(shù)據(jù). 目標(biāo)狀態(tài)變量由位置、速度和轉(zhuǎn)彎率組成，表示為Xk=，其中ωk為第k時刻目標(biāo)的轉(zhuǎn)彎率.目標(biāo)狀態(tài)轉(zhuǎn)移矩陣和過程噪聲分別為

圖6 目標(biāo)GPS定位和觀測航跡

其中，T為雷達的采樣時間間隔，取為1 s；過程噪聲標(biāo)準(zhǔn)差為σv=0.1 km/s2，σω=0.1 rad/s2. 第k時刻的雷達測距rk和測角θk可由下式計算：

其中，(Xs，Ys)是雷達的位置坐標(biāo)，觀測噪聲標(biāo)準(zhǔn)差分別為σr=0.1 km和σθ=3 mrad.

圖7、圖8 和圖9 分別給出了IMMEKF，IMMUKF，MMRBPF 和本文提出的SIQF 濾波算法在二維空間、X軸和Y軸運動方向上跟蹤機動目標(biāo)的位置均方根誤差. 從該4 種濾波算法的均方根誤差趨勢的定性比較可知：相比傳統(tǒng)交互多模型濾波算法和MMRBPF 濾波算法，SIQF 算法的跟蹤效果在濾波精度方面表現(xiàn)出較大優(yōu)勢. 這主要是因為SIQF 算法采用最大相關(guān)信息熵準(zhǔn)則衡量截斷先驗和后驗反饋約束信息，在多模型參數(shù)下構(gòu)建的混合高斯目標(biāo)重要性函數(shù)覆蓋了多似然信息，有效減小了目標(biāo)重要性函數(shù)和實際目標(biāo)狀態(tài)的真實分布之間的偏差，從而提高了系統(tǒng)狀態(tài)在動態(tài)演化過程中對隨機噪聲的抗干擾能力.

圖7 二維空間的均方根誤差

圖8 X軸方向的均方根誤差

圖9 Y軸方向的均方根誤差

另外，為了分析SIQF算法的計算復(fù)雜度. 表1統(tǒng)計了IMMEKF，IMMUKF，MMRBPF 和SIQF 這4 種濾波器執(zhí)行100 輪蒙特卡羅實驗的誤差均值和所需的平均執(zhí)行時間. 從這2 個參數(shù)的定量比較可知：相比傳統(tǒng)的交互多模型濾波算法，SIQF 濾波算法的蒙特卡羅實驗執(zhí)行時間有所增加，主要原因在于近似有界似然和構(gòu)建重要性函數(shù)過程中的迭代約束優(yōu)化. 相比MMRBPF 濾波算法，SIQF 濾波算法取得與之相當(dāng)?shù)臑V波精度，而計算時間降低了一個數(shù)量級. 這主要是因為引入最大相關(guān)信息熵調(diào)制多模似然和目標(biāo)真實分布之間的匹配程度，而不是依靠增多模型個數(shù)以表征目標(biāo)機動行為，從而避免了多模型之間不必要的競爭. 同時，引入特征輔助變量加強了該算法學(xué)習(xí)有效似然先驗知識的自適應(yīng)性，在理論上實現(xiàn)了狀態(tài)空間降維.

表1 100輪蒙特卡羅輪實驗的濾波誤差均值平均執(zhí)行時間

5 結(jié)束語

本文分析了非線性濾波中觀測似然、重要性函數(shù)與目標(biāo)真實分布失配的根本原因，為了消減根源于模型歧義和預(yù)測有偏的2種模糊性，根據(jù)軟約束理論修正有界似然，引入了最大相關(guān)信息熵準(zhǔn)則構(gòu)建覆蓋多似然模態(tài)的貝葉斯序貫重要性積分濾波算法. 仿真實驗證明，該算法充分融合Gauss-Hermite 的高精度和重要性采樣算法的靈活性，提高了貝葉斯濾波器抗擊多峰似然和非高斯噪聲干擾的能力. 結(jié)合本文作者的課題研究，后續(xù)工作將基于空間目標(biāo)星圖和光度曲線信息，研究分布式貝葉斯重要性積分算法跟蹤和反演空間目標(biāo)姿態(tài).