帶不確定混合噪聲系統的變分貝葉斯期望最大濾波算法

馬天力，張 揚，陳超波

（1.西安工業大學電子信息工程學院，西安 710021；2.陜西省自主系統與智能控制國際聯合研究中心，西安 710021）

狀態估計的目的是從受到噪聲污染或者干擾的觀測信號中，消除噪聲以及干擾的影響，準確獲得最優的系統狀態[1]。典型的濾波算法是采用卡爾曼濾波器(Kalman filter, KF)，其假設量測噪聲已知且服從單一高斯分布。但在實際過程中，量測噪聲不僅包含概率特性未知的隨機噪聲[3]，還可能存在邊界未知的有界噪聲[4,5]。例如GPS/INS組合導航系統中，GPS解算誤差屬于非高斯隨機噪聲，而INS誤差是已知其邊界的有界噪聲。

對于存在未知隨機噪聲系統的狀態估計，主要采用自適應濾波技術[6]，該方法是一種基于貝葉斯理論的遞歸估計處理策略。需要已知隨機變量的概率統計特性，并在遞歸計算噪聲參數的同時進行狀態估計。其包括狀態增廣法[7]、交互式多模型理論[8]以及序貫蒙特卡洛方法[9]。另一類方法是基于期望最大（Expectation Maximization, EM）理論的批處理技術[10]，從包含有隱變量的量測數據中計算最大似然參數。對于包含未知但有界（Unknown-But-Bounded,UBB）噪聲系統的狀態估計，由于噪聲統計特性未知，僅已知噪聲所在范圍的邊界信息，通常采用集員濾波（Set membership filter, SMF）[11]進行求解，與概率方法不同，其以包含系統真實狀態的外定界橢球集合為基礎，估計得到的是狀態變量的可行集而不是后驗概率密度，該可行解集由與系統模型、初始狀態和UBB噪聲假設相一致的所有狀態點構成。

上述方法僅針對存在某一種不確定噪聲進行計算，對于兩種或多種未知噪聲同時存在的混合不確定噪聲系統，一方面非高斯噪聲存在致使集員濾波器選取噪聲邊界困難，過小的噪聲邊界造成估計精度下降；另一方面，因集員噪聲統計特性為均勻分布，采用相應的貝葉斯濾波算法則無法獲得其解析解。解決混合噪聲條件下的狀態估計問題的核心是如何將兩類不確定噪聲進行統一描述。現有方法主要包括三類，一種是基于序貫蒙特卡洛采樣理論。Zou[12]提出箱粒子濾波算法（Box-Particle filter, BPF），該方法在粒子濾波基礎上結合區間分析理論對非高斯噪聲干擾下的狀態進行估計，但受到箱粒子數目的影響，越大的箱粒子數目，算法運行時間越長。一種是基于隨機集合理論[13-15]，將加性集員不確定納入廣義似然函數中。Handbeck[16]利用融合隨機集合提出統計與集合理論信息濾波器（Statistical and Set theory Information filter,SSI），傳統單一的濾波方法得到的值是SSI估計器的邊界情況，即當隨機誤差消失，其收斂于集合估計；當有界誤差為零，其收斂于貝葉斯估計。對于混合噪聲，其估計結果是統計狀態下具有不確定邊界的解集。胡斌[17]將隨機觀測集看作量化單元，根據量測集識別理論對于量化單元進行計算。另一種策略是非精確概率描述[18,19]，其用概率密度函數的閉型凸集對混合噪聲進行構建。Noack[20]將混合不確定噪聲表示為概率密度的集合，將系統狀態轉移與似然密度的信度需求分別擴展到狀態轉移與似然密度凸集的形式，提出信度狀態濾波器（Credal state filter, CSF）。Klumpp[21]對基于隨機集合的SSI濾波器和基于集合概率密度的CSF進行了比較，CSF估計結果相比SSI更保守。此外，Ginger[22]結合序貫蒙特卡洛理論提出非精確采樣概念，運用混合非精確采樣集合對非精確量測進行建模，并在貝葉斯框架下通過集員理論中有界支撐傳播非精確信息。但該方法對參數具有較強的敏感度。Henningsson[23]根據控制中觀測器的思想，用橢球包含混合噪聲中有界集合部分的估計誤差，并通過線性矩陣不等式得到最優濾波增益，該算法的性能受到一個權系數的影響，而該系數的取值取決于實際系統總體噪聲中隨機不確定和有界不確定的某種權重關系，但一般這種權重關系并不明確。江濤[24]提出混合噪聲的聯合濾波算法（Combined filter, CF），其將統計特性未知但有界的噪聲引入到卡爾曼濾波模型中，得到一組包含集合運算的卡爾曼濾波方程。其中UBB噪聲應用SMF的思想進行處理，而隨機噪聲應用卡爾曼濾波的思想進行處理。該算法問題在于統計噪聲為參數已知的高斯分布。

本文針對未知隨機不確定噪聲和有界不確定噪聲共同存在下的狀態估計問題，提出基于混合噪聲框架的狀態估計方法，將隨機和有界不確定噪聲運用未知參數的混合統計噪聲模型進行描述，將混合噪聲條件下的狀態估計問題描述為未知參數聯合統計噪聲條件下的參數與狀態的估計問題，并設計變分期望最大算法求解其狀態及噪聲參數。

1 系統模型

假設含有隨機噪聲和有界噪聲的不確定系統的狀態方程為

其中xk為k時刻系統狀態向量，xk?Rn；wk為零均值，方差為Qk的高斯過程噪聲，其中Qk為非負正定協方差矩陣。zk為k時刻的系統量測，vk+ek為混合測量噪聲，其中vk為零均值非高斯隨機量測噪聲，ek為UBB噪聲，ek∈Ek?Rm，Ek為邊界未知的有界集合。假設在上述模型中，ek= 0，即量測噪聲為典型的隨機模型，該模型可以用卡爾曼濾波器進行求解。

隨機、集員和混合不確定統計特性模型的概率密度函數表示如圖1所示。由于有界和隨機并沒有嚴格的劃分，兩者可以近似轉換。對于有界噪聲，若利用統計特性進行表述，則高斯白噪聲方差σ為有界噪聲邊界的1/3倍[24]。

圖1 隨機與集員混合模型Fig.1 Combined stochastic and set-membership model

基于該特性，有界不確定噪聲可近似描述為參數未知的高斯噪聲，并構建兩種不確定噪聲的混合高斯模型。

由于高斯混合項中含有隱變量、且每個高斯分布的均值、方差和權重系數均未知。若采用貝葉斯方法估計，則需考慮參數先驗分布。對于混合權重系數lα，其先驗為Dirichlet分布[10]。

其中λ0是濃度參數。噪聲高斯項均值和方差的先驗服從Gaussian和Wishart分布。

其中υ0為Wishart分布的自由度，?0為精度矩陣。

令 Θ =[α,μ, Σ,S]，本文的目的是在已知觀測集合Z={z1,z2,… ,zN}的情況下，求解最優系統狀態X={x0,x1,…xn}和參數集合Θ的最大后驗分布。

由狀態空間模型可知：

系統的似然函數為：

2 變分EM濾波算法

2.1 變分貝葉斯學習

對于未知參數Θ的估計，若系統量測Z已知，則根據貝葉斯規則：

一般情況下，邊緣概率密度函數p(Z)的積分形式復雜，難以求得其解析解，導致難以得到參數后驗概率密度函數的解析表達式p(X, Θ|Z)。針對這一問題，根據先驗分布滿足共軛指數分布族，變分貝葉斯（Variation Bayesian, VB）理論提出構建一個新的分布q(X,Θ)去逼近真實后驗分布p(X, Θ|Z)，無需給定其概率密度函數形式，表達形式與p(X, Θ|Z)一致，區別在于參數不同。故邊緣似然函數的log形式表示為：

其中，自由能量函數 F (q(X,Θ)是logp(Z)的下界函數[6]。K L(q(X, Θ) ||p(X, Θ |Z))為q(X,Θ)和p(X,Θ|Z)之間的Kullback-Leibler散度，其用來衡量兩個概率分布之間的差異。因KL散度非負，由式(11)可知，只需使KL散度最小，求最大化自由能量函數max F (q(X,Θ) )，即可獲得最佳近似分布。為了使模型參數可求，VB采用均值域逼近理論將多變量聯合概率分布近似逼近為各個變量邊緣分布的乘積[25]。q(X,Θ)的表達式如下：

對于每一個參數的近似分布q( Θi)對KL散度

求偏導，即可求得q(Θi)的通解表達式

2.2 變分EM估計策略

由于在聯合后驗概率分布計算過程中，所求參數相互耦合。故本文將EM理論中參數迭代計算后驗參數思想引入VB理論，設計參數耦合更新策略，通過反復迭代求解最大期望。在變分貝葉斯期望（Variational Bayesian expectation, VBE）步驟中，利用前向-后向遞推算法計算隱變量的后驗，并在變分貝葉斯最大化（Variational Bayesian maximization, VBM）步驟中對模型超參數進行更新。

A 變分貝葉斯期望步驟

根據文獻[25]，利用前向遞推算法求解隱變量的充分統計量，通過輸入隱馬爾可夫模型參數矩陣以及觀測序列計算該序列發生的概率。假設ak(xk) =p(xk|z1:k)為隱變量的后驗概率密度，結合公式(14)可知狀態xk的邊緣后驗為：

其中ζk=p(zk|zk-1)為歸一化因子。a(xk)中關于xk-1邊緣概率密度函數滿足高斯分布形式，其均值為，方差為經過簡化：

將上式帶入式(15)中，并對式(15)進行邊緣化處理，使其只含有隱狀態xk項，便能夠計算高斯分布的均值θk和方差Ωk，如下所示

與前項遞推濾波過程不同，后向濾波根據觀測量zk+1:N估計現有狀態xk的后驗，采用并行遞推的形式實現。假設bk(xk) =p(zk+1:N|xk)，xk滿足終止條件bN(xN) = 1。則：

與ak(xk)類似，b(xk)同樣滿足高斯分布的形式，則：

近似分布q(xk)的后驗概率密度計算采用前向-后向聯合濾波，如下所示：

則系統狀態后驗的充分統計量為：

B 變分貝葉斯最大化步驟

一旦獲得隱變量的最大后驗分布，在VBM步驟中需要計算系統模型的超參數。由于參數變量空間維度高導致其計算困難，所提算法利用近似邊緣分布去分別逼近其統計特性。根據共軛指數分布的性質[26]，高斯混合模型參數的近似邊緣分布定義為：

對于高斯混合模型系數αl的推理，假設其變分后驗分布與先驗分布形式相同，均滿足Dirichlet分布，根據式(13)(14)，對聯合概率密度函數的近似分布取對數后為：

由于混合系數的后驗概率密度形式與先驗相同，經過變分推斷：

同理，對于混合噪聲高斯項均值μl和方差Σl的聯合概率密度函數近似分布取對數：

其中Tr(·)表示矩陣的跡，根據式(6)，q(μl, Σl)的后驗分布為高斯威沙特分布的形式。首先考慮均值μl：

上式展開后，合并關于μl的同類項，得到新的高斯分布q(μl|Σl) =N(μl|ml,(βlΣl)-1)，其中：

與混合系數的計算類似，上式為遞推形式。由式(6)可知，lμ和Σl相互耦合，無法利用式(30)直接進行求解。因此，采用兩者聯合概率密度函數減去lμ的概率密度函數。

對上式中跡矩陣進行合并簡化，其中高斯混合項方差的后驗近似分布為，則：

至此，式(24)-(36)為VBEM算法的狀態及參數更新方程。式(24)(25)中狀態均值與方差的計算與超參數有關，式(33)(35)中超參數的更新與狀態均值、方差同樣相關，相互耦合。在VBEM中，需要迭代求解狀態與參數，直至達到收斂。該方法與EM算法類似，其均是尋求能量函數關于隱變量分布最大值的方法。不同之處在于，EM算法在期望最大化步驟中，并未賦予參數任何先驗信息，而是作為未知隨機變量進行計算。而VBEM算法考慮參數近似后驗分布q(Θ)，對參數的超參數進行更新計算。若將所求參數分布假設為Dirichlet函數，則VBEM對概率模型參數進行點估計求解，所提出的VBEM算法簡化為EM算法。

3 仿真驗證

為了驗證所提算法的有效性，分別采用卡爾曼濾波器（KF），聯合濾波（CF），箱粒子濾波算法（BPF）以及VBEM算法對GPS/INS松耦合系統開展仿真實驗。取東北天地理坐標系為導航坐標系，以INS誤差作為狀態變量，GPS和INS的位置差作為觀測變量，系統狀態方程為：

其中x=[φE,φN,φU,δVE,δVN,δVU,δL,δλ,δh,εE,εN,εU,?E,?N,?U]T，φE、φN和φU為姿態誤差角；δVE、δVN和δVU為東北天向的速度誤差；δL、δλ和δh為經緯度及高度誤差，εE、εN與εU為陀螺的一階馬爾科夫漂移誤差；?E、?N與 ?U為加速度計零點漂移誤差。F為系統狀態矩陣。

系統量測方程為：

其中z=[VEG-VEI,VNG-VNI,VUG-VUI,LG-LI,λG-λI,hG-hI]T；LG、LI、λG、λI、hG、hI分別為GPS和INS觀測的位置信息；VEG、VEI、VNG、VNI、VUG、VUI分別為GPS和INS獲得的速度分量，Hk=[03×3,I6×6,03×3]為觀測矩陣。

模擬車輛在二維空間運動，其行駛軌跡如圖2所示。車輛初始位置為北緯34.246 °，東經108.909 °，初始速度為10 m/s；初始姿態為方位角0 °，俯仰角為0 °，橫滾角為0 °。仿真時間為600 s。vk為服從Student-T分布的隨機噪聲序列，均值為0，自由度為1，協方差矩陣δL=δλ= 10-5(°)。ek=[ek,L,ek,λ]T為 量 測 UBB噪 聲 ，ek,L,ek,λ?[-1 ,1]× 1 0-4(°)。

圖2 車輛真實軌跡Fig.2 Trajectory of the vehicle

圖3為四種算法100次蒙特卡洛仿真實驗下的經緯度的誤差結果。從圖中可以看出，對于同時存在未知參數統計噪聲和未知邊界的有界噪聲的組合導航系統，變分EM算法的經緯度均方根誤差明顯低于其他三種算法。主要原因是VBEM將混合噪聲建模為未知參數的高斯混合模型，進行聯合估計。而KF算法僅針對量測噪聲參數已知的高斯分布模型。因此，結果存在較大偏差，估計性能不足。另外，CF算法與BPF算法均是采用區間分析策略，該策略需要已知噪聲的上、下界。由于噪聲模型有偏，不準確的噪聲邊界設定造成估計結果出現誤差。表1為四種算法的經緯度均方根誤差(RMSE)統計表。

圖3 不同算法的經緯度誤差Fig.3 The latitude and longitude errors from different method

表1 四種算法的經緯度誤差統計表Tab.1 Statistical table of latitude and longitude error

從表中可以看出，VBEM算法相比于KF，CF和BPF在緯度的平均均方根誤差分別下降了83.32%、80.81%和79.84%。在經度方面，VBEM算法相比于KF，CF和BPF平均均方根誤差分別下降了61.28%、60.97%、56.78%。

本文對算法的運算時間進行了分析，表2為四種算法的相對運算仿真時間對比，令卡爾曼濾波算法運算時間為單位1。由表2可知，本文所提算法雖然相比于KF，CF具有較長的運算時間，但明顯低于基于蒙特卡洛理論的BPF算法。

表2 四種算法的相對運行時間Tab.2 The running time of different method

為了驗證VBEM算法中相關參數對系統的影響，分別從初始量測噪聲均值以及精度矩陣對系統影響角度進行了分析。圖4為不同初始噪聲均值mL，mλ下的經緯度誤差。從圖中可以看出，在噪聲初始均值不相同的情況下，采用變分EM算法均可使系統收斂。初始噪聲均值越貼近真實值，其經緯度誤差越小。圖5為不同初始精度矩陣Λ下的濾波結果。可以看出，初始噪聲方差越接近真實值，經緯度的均方根誤差越小，主要是因為VBEM參數遞歸中需要引入初始噪聲參數進行計算，較大誤差的先驗信息會對濾波系統造成干擾。

圖4 不同初始噪聲均值下的經緯度誤差Fig.4 The latitude and longitude errors from different initial noise means

圖5 不同初始精度矩陣下的經緯度誤差Fig.5 The latitude and longitude errors from different initial precision

為了說明所提算法中遞歸循環次數對系統的影響，進行了相關驗證。采用不同迭代次數分別計算整個濾波過程的均方根誤差。從圖6中可以看出，所提算法的收斂效果受到迭代次數的影響。迭代次數越低，相對應濾波算法的誤差越大。迭代次數為20時，在遞歸計算過程中自由能量函數并沒有達到最大，造成估計結果有偏。

圖6 不同迭代次數下的經緯度誤差Fig.6 The latitude and longitude errors from different iteration times

3 結 論

針對混合不確定量測噪聲條件下的狀態估計問題，提出基于變分貝葉斯期望最大理論的魯棒濾波算法。由于隨機噪聲統計特性與有界噪聲邊界參數均未知，若僅采用基于統計信息或基于區間分析的策略，易造成誤差偏移。因此，本文將隨機不確定噪聲與未知有界噪聲運用未知參數的高斯混合模型進行統一表示，通過變分期望算法對模型隱變量進行推理，同時運用變分最大化算法更新模型參數。最終獲得狀態估計結果。仿真試驗結果表明，在含有混合不確定噪聲的情況下，所提算法具有較高的估計精度。鑒于變分EM具有較強的自適應性，未來考慮非線性系統模型以及具有復雜時變特性的隨機及有界噪聲模型。