閃爍噪聲下的自適應無跡卡爾曼算法

陳 政，劉 康

（昆明理工大學信息工程與自動化學院，云南 昆明 650500）

0 引言

卡爾曼最初提出的濾波基本理論只適用于線性系統，并且要求量測也必須是線性的[1]。卡爾曼算法可以看成是一種貝葉斯方法，貝葉斯方法通過對模型引入先驗分布，并利用似然函數求解后驗分布的方法來對模型中的不確定性進行建模[2]。在線性估計算法的發展歷程中，卡爾曼濾波扮演了重要的角色，但該算法實現的前提是系統模型已知，與實際工業中的狀態估計問題不相符[3]；同時，真實環境中衛星雷達受到的非高斯閃爍噪聲進一步加重了后驗概率密度的非高斯特性[4]，于是，傳統的卡爾曼濾波已經不再適用[5][6]。由此帶來了改變后的使用到非線性系統觀測中的濾波理論，如無跡卡爾曼濾波器（UKF），可用于導航、目標狀態估計、飛行器軌道確定等領域并獲得了良好效果，但噪聲統計特性未知時，UKF濾波精度下降[7]。

全球衛星導航系統在許多領域都發揮出了重要的作用，并且仍然在不斷地改進與創新[8]。近年來，隨著我國空間技術的發展，在軌運行的衛星越來越多[9]。衛星偵察范圍廣，可對目標進行長時間、大范圍的連續偵察和監視，獲取情報的時效性強，是現代偵察中不可或缺的重要手段[10]。人造衛星在空中飛行，當自身攜帶的能源耗盡時便被回收。基于此，設計出來一個復雜計算度低、精確度高的非線性動態系統的濾波算法尤為重要，本文結合EM算法和無跡卡爾曼濾波算法，提出一種新的方法，仿真結果表明，本文提出的方法濾波精度較高，能夠有效的跟蹤目標。

1 問題描述

1.1 閃爍噪聲

隨著計算機技術的不斷發展，利用計算機來實現目標追蹤功能成為目前計算機領域中最熱門的課題之一[11]。目標跟蹤是計算機視覺領域的一個重要研究方向，被廣泛應用于智能監控、目標識別、交通監視等方面[12]。

閃爍噪聲分布跟高斯分布的主要區別在于尾部較長，其中心區域則類似高斯分布[2]。本文采用混合高斯模型來模擬閃爍噪聲，其概率密度如下：

其中，N(x; μi,Pi)表示均值為μi，方差為 Pi的高斯分布，閃爍強度 ε ∈ [ 0,1]。

1.2 運動模型

目標跟蹤的只要目的就是估計移動目標的狀態軌跡[13]。考慮一般的目標跟蹤問題，在笛卡爾坐標系中，目標做勻速運動的離散時間狀態系統方程和觀測方程分別為：

其中，kx表示目標在K時刻的狀態向量，yk表示K時刻狀態的觀測值，F為狀態轉移矩陣，L為過程噪聲驅動矩陣，在一般的雷達目標跟蹤中，()·h為有界的非線性函數，過程噪聲1-kw 是均值為零方差為Q的高斯白噪聲，實際情況中，Q經常是未知的，觀測噪聲kv是均值為零方差為R的高斯白早噪聲，實際情況中，kv經常是閃爍噪聲。

2 算法描述

2.1 無跡卡爾曼濾波算法

考慮如下含有加性高斯噪聲的離散非線性系統[14-15]：

其中， k ≥ 0 為離散時間變量， x ∈ Rn為狀態向量，z ∈ Rm為輸出向量。非線性函數 fn∈Rm，h ∈ Rn。過程噪聲 wk-1和量測噪聲 vk分別為n維和m維的高斯白噪聲，并具有以下統計特性：

在很多情況下，系統噪聲都是非零均值。因此，給出噪聲均值非零情況下的 UKF 算法。由于系統過程噪聲和量測噪聲為非零均值的高斯噪聲且其均值分別為q和r。則系統噪聲可寫為如下形式：

顯然，kμ和ηk為互不相關的零均值高斯白噪聲，且其方差分別為Q和R。

此時, 離散非線性系統 (4) 可以被等價表述為以下形式：

則UKF 算法的計算步驟如下：

步驟1. 初始化系統變量

步驟2. 時間更新

根據確定的采樣策略，得到服從均值 xk-1、方差 pk-1的Sigma 點 { χi,k-1},i = 0,...,2n，為：

Sigma點的對應權值為：

其中，α描述了 Sigma 點的散布程度,通常取一小的正值（如 0.01），k通常取0；β用來描述 x的分布信息（Gauss情況下β的最優值為 2）；()i表示矩陣 ( n + λ )Pk-1的平方根矩陣的第k-1列； w(m)(i = 0,1,2,...,2n)為求一階統計特性時

i

的權系數 w(c)(i = 0,1,2,...,2n)為求二階統計特性時

i的權系數.計算 Sigma 點 {χi,k?1}, i = 0,··· ,2n 沿非線性函數 f ( .)的傳播結果， 并考慮非零均值噪聲的作用。

從而得到隨機向量沿非線性函數 (.)f 傳播的后驗均值和協方差：

步驟3. 量測更新

其對應權值與式 (9) 的權值相同。

計算 Sigma 點｛χi,k|k-1｝, i = 0,··· ,2n 沿非線性函數 h( .)的傳播結果，并考慮非零均值噪聲的作用。

得到隨機向量沿非線性函數 (.)h 傳播的后驗均值、自協方差和互協方差：

步驟4. 濾波更新

根據量測數據kz得到最小方差估計結果為：

步驟4 完成后，則返回步驟2，開始下一時刻的濾波過程。

2.2 EM算法在目標跟蹤中的應用

首先給出一般形式的線性狀態空間模型的參數估計問題：

其中F為狀態轉移矩陣，H為觀測矩陣， wt，vt為過程噪聲和觀測噪聲，它們分別服從以下分布：

對運動目標模型的參數估計即在知道tx和ty的情況下（將tx和ty合起來看做完全數據）對目標的運動參數F、H、Q、R進行估計[16]。設參數的F、H、Q、R的對數似然函數為：

多維高斯分布的概率密度函數：

且 xt+1～ N(F xt,Q)， yt～ N(H xt,R),即完全數據的似然也服從高斯分布，將這兩個式子帶入式（26）中，得：

將上式展開，我們就可以得到對數似然函數的最終表達式：

由EM算法的基本原理可知，需要對（28）式進行最大化處理，即：

化簡后得：

EM 算法是一種迭代算法，因此我們可以先假定F,H,Q,R的值，得到tx和ty的值后進行平滑濾波（E步），得到新的tx和ty的值后再利用式（30-33）估計出參數F,H,Q,R新的值（M步），然后再回到E步，這樣不斷的重復此過程直到參數收斂[17]。

3 仿真實驗

3.1 仿真場景

為了驗證本文提出的基于EM算法的改進無跡卡爾曼濾波算法（下面簡稱EMUKF）的濾波效果，在過程噪聲統計特性未知和帶有閃爍噪聲的情況下，分別采用UKF算法和EMUKF算法進行仿真。二維情況下，設雷達觀測站位于坐標原點，目標做近似勻速直線運動，運動方程描述如下：

向量 x = [ x x˙ y y˙]′，其中x和x˙分別表示X軸上的位置分量和速度分量，y和y˙分別表示Y軸上的位置分量和速度分量，目標的初始狀態為 x0=[50 1 50 - 1 ]′， zk是觀測向量， rk和θk分別是目標到觀測站的距離和角度。仿真中過程噪聲是均值為零的高斯白噪聲，觀測噪聲是閃爍噪聲，過程噪聲的統計特性都是未知。仿真中過程噪聲統計特性如下：

仿真中先驗過程噪聲 Qk=diag ｛ 0. 01 , 0.01｝。仿真時間為100 s，采樣周期T=1s。

3.2 仿真結果

針對過程噪聲統計特性未知，過程噪聲是一般高斯白噪聲和觀測噪聲是閃爍噪聲這兩種情況分別做 100次蒙特卡洛仿真。仿真實驗的均方根誤差（Root mean square error, RMSE）定義為：

其中N為蒙特卡洛仿真次數。

圖1為仿真中真實軌跡和兩種算法濾波軌跡的對比，圖2、3為仿真下兩種算法的均方根誤差曲線比較。表1是仿真下兩種濾波算法均方根誤差的均值和方差的統計數據。

圖1 三條不同軌跡Fig.1 Three different tracks

圖2 位置均方根誤差Fig.2 Root mean square error about position

通過上述仿真結果可以看出閃爍噪聲和過程噪聲統計特性未知會造成濾波效果差甚至發散的情況。本文提出的算法在閃爍噪聲下比標準的無跡卡爾曼濾波具有更良好的濾波效果，其均方根誤差也明顯小于標準的無跡卡爾曼濾波算法。

圖3 速度均方根誤差Fig.3 Root mean square error about velocity

表1 算法性能比較Tab.1 Algorithm performance comparison

4 結論

針對雷達閃爍噪聲情況下，過程噪聲統計特性未知的非線性目標跟蹤問題，采用無跡卡爾曼濾波的同時運用EM算法對過程噪聲進行參數估計，獲得較準確的參數后明顯的提高了濾波的精度和穩定性。特別是在閃爍噪聲的情況下，本文提出的自適應無跡卡爾曼濾波算法要明顯優于標準的無跡卡爾曼濾波，表現出很強的抗噪聲性能。

[1] 石章松, 劉忠等編著, 目標跟蹤與數據融合理論及方法[M]. 北京: 國防工業出版社, 2010. 7.

[2] 侯禹騰. 數據增廣求解貝葉斯Logistic 回歸模型的方法研究[J]. 軟件, 2014, 35(7): 109-115.

[3] 王寫, 劉妹琴, 張森林.卡爾曼濾波器、H∞濾波器和魯棒混合Kalman/H∞濾波器的比較[J].新型工業化, 2012, 2(7):36-45.

[4] I. Bilik; J. Tabrikian. Target tracking in glint noise environment using nonlinear non-Gaussian Kalman filter. 2006 IEEE Conference on Radar. 2006: 6 pp.

[5] Sorenson H W. Kalman filtering: Theory and Application.New York: IEEE, 1985.

[6] F. Daum. Nonlinear filters: beyond the Kalman filter. IEEE Aerospace and Electronic Systems Magazine.2005, 20(8):57 –69.

[7] 秦永元, 張洪鉞, 汪淑華等編著, 卡爾曼濾波與組合導航原理(第3版)[M]. 西安:西北工業大學出版社, 2015. 6.

[8] 陳楊毅, 陳凌宇.基于最優幾何的改進型多星座衛星選星算法[J]. 軟件, 2014,35(1): 26-31.

[9] 楊冬, 孫劍偉. 基于令牌桶算法的衛星數據地面傳輸流量控制方法研究[J]. 軟件, 2016, 37(3): 99-103.

[10] 張志恒, 尹路明, 王茂磊. 采用組合優化策略的電子偵察衛星規劃方法[J]. 軟件, 2014,35(4): 143-149

[11] 邱亞欽, 劉渭濱, 戶磊, 等.動態背景下的運動目標檢測方法[J].新型工業化, 2011,1(9): 22-29.

[12] 王江, 付強, 全權, 等.基于Kalman 濾波和直方圖匹配的雙目視覺跟蹤[J].新型工業化, 2013,3(2): 23-33.

[13] [美] Anton J。Haug著, 王欣, 于曉譯 Bayesian Estimation and Tracking A Practical Guide(貝葉斯估計與跟蹤實用指南)[M]. 北京: 國防工業出版社, 2014.

[14] 武飛, 柳炳利. Kalman 濾波技術在地球化學數據處理中的應用[J]. 軟件, 2013,34(9): 70-74.

[15] 韓崇昭, 朱洪艷, 段戰勝等著, 多源信息融合(第二版)[M].北京: 清華大學出版社, 2010.9.

[16] Ming Lei; Chongzhao Han; Panzhi Liu. Expectation Maximization (EM) algorithm-based nonlinear target tracking with adaptive state transition matrix and noise covariance. 2007 10th International Conference on Information Fusion. 2007:1–7.

[17] A.P. Dempster; N.M.Laird; D.B.Rubin. Maximum Likelihood from Incomplete Data via the EM Algorithm. Journal of the Royal Statistical Society. SeriesB (Methodological), 1977,39(1): 1-38.