陳 政，劉 康

（昆明理工大學(xué)信息工程與自動(dòng)化學(xué)院，云南 昆明 650500）

0 引言

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

全球衛(wèi)星導(dǎo)航系統(tǒng)在許多領(lǐng)域都發(fā)揮出了重要的作用，并且仍然在不斷地改進(jìn)與創(chuàng)新[8]。近年來(lái)，隨著我國(guó)空間技術(shù)的發(fā)展，在軌運(yùn)行的衛(wèi)星越來(lái)越多[9]。衛(wèi)星偵察范圍廣，可對(duì)目標(biāo)進(jìn)行長(zhǎng)時(shí)間、大范圍的連續(xù)偵察和監(jiān)視，獲取情報(bào)的時(shí)效性強(qiáng)，是現(xiàn)代偵察中不可或缺的重要手段[10]。人造衛(wèi)星在空中飛行，當(dāng)自身攜帶的能源耗盡時(shí)便被回收?；诖?，設(shè)計(jì)出來(lái)一個(gè)復(fù)雜計(jì)算度低、精確度高的非線性動(dòng)態(tài)系統(tǒng)的濾波算法尤為重要，本文結(jié)合EM算法和無(wú)跡卡爾曼濾波算法，提出一種新的方法，仿真結(jié)果表明，本文提出的方法濾波精度較高，能夠有效的跟蹤目標(biāo)。

1 問(wèn)題描述

1.1 閃爍噪聲

隨著計(jì)算機(jī)技術(shù)的不斷發(fā)展，利用計(jì)算機(jī)來(lái)實(shí)現(xiàn)目標(biāo)追蹤功能成為目前計(jì)算機(jī)領(lǐng)域中最熱門(mén)的課題之一[11]。目標(biāo)跟蹤是計(jì)算機(jī)視覺(jué)領(lǐng)域的一個(gè)重要研究方向，被廣泛應(yīng)用于智能監(jiān)控、目標(biāo)識(shí)別、交通監(jiān)視等方面[12]。

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

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

1.2 運(yùn)動(dòng)模型

目標(biāo)跟蹤的只要目的就是估計(jì)移動(dòng)目標(biāo)的狀態(tài)軌跡[13]?？紤]一般的目標(biāo)跟蹤問(wèn)題，在笛卡爾坐標(biāo)系中，目標(biāo)做勻速運(yùn)動(dòng)的離散時(shí)間狀態(tài)系統(tǒng)方程和觀測(cè)方程分別為：

其中，kx表示目標(biāo)在K時(shí)刻的狀態(tài)向量，yk表示K時(shí)刻狀態(tài)的觀測(cè)值，F(xiàn)為狀態(tài)轉(zhuǎn)移矩陣，L為過(guò)程噪聲驅(qū)動(dòng)矩陣，在一般的雷達(dá)目標(biāo)跟蹤中，()·h為有界的非線性函數(shù)，過(guò)程噪聲1-kw 是均值為零方差為Q的高斯白噪聲，實(shí)際情況中，Q經(jīng)常是未知的，觀測(cè)噪聲kv是均值為零方差為R的高斯白早噪聲，實(shí)際情況中，kv經(jīng)常是閃爍噪聲。

2 算法描述

2.1 無(wú)跡卡爾曼濾波算法

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

其中， k ≥ 0 為離散時(shí)間變量， x ∈ Rn為狀態(tài)向量，z ∈ Rm為輸出向量。非線性函數(shù) fn∈Rm，h ∈ Rn。過(guò)程噪聲 wk-1和量測(cè)噪聲 vk分別為n維和m維的高斯白噪聲，并具有以下統(tǒng)計(jì)特性：

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

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

此時(shí), 離散非線性系統(tǒng) (4) 可以被等價(jià)表述為以下形式：

則UKF 算法的計(jì)算步驟如下：

步驟1. 初始化系統(tǒng)變量

步驟2. 時(shí)間更新

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

Sigma點(diǎn)的對(duì)應(yīng)權(quán)值為：

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

i

的權(quán)系數(shù) w(c)(i = 0,1,2,...,2n)為求二階統(tǒng)計(jì)特性時(shí)

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

從而得到隨機(jī)向量沿非線性函數(shù) (.)f 傳播的后驗(yàn)均值和協(xié)方差：

步驟3. 量測(cè)更新

其對(duì)應(yīng)權(quán)值與式 (9) 的權(quán)值相同。

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

得到隨機(jī)向量沿非線性函數(shù) (.)h 傳播的后驗(yàn)均值、自協(xié)方差和互協(xié)方差：

步驟4. 濾波更新

根據(jù)量測(cè)數(shù)據(jù)kz得到最小方差估計(jì)結(jié)果為：

步驟4 完成后，則返回步驟2，開(kāi)始下一時(shí)刻的濾波過(guò)程。

2.2 EM算法在目標(biāo)跟蹤中的應(yīng)用

首先給出一般形式的線性狀態(tài)空間模型的參數(shù)估計(jì)問(wèn)題：

其中F為狀態(tài)轉(zhuǎn)移矩陣，H為觀測(cè)矩陣， wt，vt為過(guò)程噪聲和觀測(cè)噪聲，它們分別服從以下分布：

對(duì)運(yùn)動(dòng)目標(biāo)模型的參數(shù)估計(jì)即在知道tx和ty的情況下（將tx和ty合起來(lái)看做完全數(shù)據(jù)）對(duì)目標(biāo)的運(yùn)動(dòng)參數(shù)F、H、Q、R進(jìn)行估計(jì)[16]。設(shè)參數(shù)的F、H、Q、R的對(duì)數(shù)似然函數(shù)為：

多維高斯分布的概率密度函數(shù)：

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

將上式展開(kāi)，我們就可以得到對(duì)數(shù)似然函數(shù)的最終表達(dá)式：

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

化簡(jiǎn)后得：

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

3 仿真實(shí)驗(yàn)

3.1 仿真場(chǎng)景

為了驗(yàn)證本文提出的基于EM算法的改進(jìn)無(wú)跡卡爾曼濾波算法（下面簡(jiǎn)稱(chēng)EMUKF）的濾波效果，在過(guò)程噪聲統(tǒng)計(jì)特性未知和帶有閃爍噪聲的情況下，分別采用UKF算法和EMUKF算法進(jìn)行仿真。二維情況下，設(shè)雷達(dá)觀測(cè)站位于坐標(biāo)原點(diǎn)，目標(biāo)做近似勻速直線運(yùn)動(dòng)，運(yùn)動(dòng)方程描述如下：

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

仿真中先驗(yàn)過(guò)程噪聲 Qk=diag ｛ 0. 01 , 0.01｝。仿真時(shí)間為100 s，采樣周期T=1s。

3.2 仿真結(jié)果

針對(duì)過(guò)程噪聲統(tǒng)計(jì)特性未知，過(guò)程噪聲是一般高斯白噪聲和觀測(cè)噪聲是閃爍噪聲這兩種情況分別做 100次蒙特卡洛仿真。仿真實(shí)驗(yàn)的均方根誤差（Root mean square error, RMSE）定義為：

其中N為蒙特卡洛仿真次數(shù)。

圖1為仿真中真實(shí)軌跡和兩種算法濾波軌跡的對(duì)比，圖2、3為仿真下兩種算法的均方根誤差曲線比較。表1是仿真下兩種濾波算法均方根誤差的均值和方差的統(tǒng)計(jì)數(shù)據(jù)。

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

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

通過(guò)上述仿真結(jié)果可以看出閃爍噪聲和過(guò)程噪聲統(tǒng)計(jì)特性未知會(huì)造成濾波效果差甚至發(fā)散的情況。本文提出的算法在閃爍噪聲下比標(biāo)準(zhǔn)的無(wú)跡卡爾曼濾波具有更良好的濾波效果，其均方根誤差也明顯小于標(biāo)準(zhǔn)的無(wú)跡卡爾曼濾波算法。

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

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

4 結(jié)論

針對(duì)雷達(dá)閃爍噪聲情況下，過(guò)程噪聲統(tǒng)計(jì)特性未知的非線性目標(biāo)跟蹤問(wèn)題，采用無(wú)跡卡爾曼濾波的同時(shí)運(yùn)用EM算法對(duì)過(guò)程噪聲進(jìn)行參數(shù)估計(jì)，獲得較準(zhǔn)確的參數(shù)后明顯的提高了濾波的精度和穩(wěn)定性。特別是在閃爍噪聲的情況下，本文提出的自適應(yīng)無(wú)跡卡爾曼濾波算法要明顯優(yōu)于標(biāo)準(zhǔn)的無(wú)跡卡爾曼濾波，表現(xiàn)出很強(qiáng)的抗噪聲性能。

[1] 石章松, 劉忠等編著, 目標(biāo)跟蹤與數(shù)據(jù)融合理論及方法[M]. 北京: 國(guó)防工業(yè)出版社, 2010. 7.

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

[3] 王寫(xiě), 劉妹琴, 張森林.卡爾曼濾波器、H∞濾波器和魯棒混合Kalman/H∞濾波器的比較[J].新型工業(yè)化, 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] 秦永元, 張洪鉞, 汪淑華等編著, 卡爾曼濾波與組合導(dǎo)航原理(第3版)[M]. 西安:西北工業(yè)大學(xué)出版社, 2015. 6.

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

[9] 楊冬, 孫劍偉. 基于令牌桶算法的衛(wèi)星數(shù)據(jù)地面?zhèn)鬏斄髁靠刂品椒ㄑ芯縖J]. 軟件, 2016, 37(3): 99-103.

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

[11] 邱亞欽, 劉渭濱, 戶(hù)磊, 等.動(dòng)態(tài)背景下的運(yùn)動(dòng)目標(biāo)檢測(cè)方法[J].新型工業(yè)化, 2011,1(9): 22-29.

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

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

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

[15] 韓崇昭, 朱洪艷, 段戰(zhàn)勝等著, 多源信息融合(第二版)[M].北京: 清華大學(xué)出版社, 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.