（海軍航空工程學(xué)院 電子信息工程系，山東 煙臺(tái) 264001）

0 引言

擴(kuò)展卡爾曼濾波（EKF）是在工程實(shí)踐中應(yīng)用較為廣泛的濾波方法。該方法利用泰勒級(jí)數(shù)展開，將非線性濾波問(wèn)題轉(zhuǎn)化為線性濾波問(wèn)題。但是，EKF濾波在進(jìn)行濾波前假定過(guò)程噪聲和量測(cè)噪聲為互不相關(guān)零均值白噪聲，且噪聲的統(tǒng)計(jì)特性是已知的。噪聲的不準(zhǔn)確將會(huì)導(dǎo)致EKF方法給出的均方誤差偏離真實(shí)情況，有時(shí)不能準(zhǔn)確反映狀態(tài)估計(jì)精度。因此，在先驗(yàn)信息不完全的情況下，采用與真實(shí)值相差較多的觀測(cè)噪聲標(biāo)準(zhǔn)差進(jìn)行EKF 濾波，得到的狀態(tài)濾波方差有時(shí)無(wú)法反映狀態(tài)濾波精度，狀態(tài)濾波結(jié)果并不總是可靠的。同時(shí)在實(shí)際應(yīng)用中，動(dòng)態(tài)系統(tǒng)隨時(shí)可能受到各種外界干擾，因此觀測(cè)噪聲的統(tǒng)計(jì)特性隨時(shí)會(huì)發(fā)生變化，僅僅依靠有限的先驗(yàn)信息對(duì)觀測(cè)噪聲進(jìn)行描述是不可靠的，需要一種能夠?qū)τ^測(cè)噪聲進(jìn)行實(shí)時(shí)估計(jì)的方法[1-3]。

1 EMD方法

經(jīng)驗(yàn)?zāi)B(tài)分解方法（EMD）從本質(zhì)上講是對(duì)一個(gè)信號(hào)（或其導(dǎo)數(shù)，視所需的分解精度而定）進(jìn)行平穩(wěn)化處理，其結(jié)果是將信號(hào)中不同尺度的波動(dòng)或趨勢(shì)逐級(jí)分解開來(lái)，產(chǎn)生一系列具有不同特征尺度的數(shù)據(jù)序列，每一個(gè)序列被稱為一個(gè)本征模態(tài)函數(shù)（Intrinsic Mode Function，IMF）。最低頻率的IMF分量通常情況下代表原始信號(hào)的趨勢(shì)或均值。作為一種運(yùn)用，EMD方法可以有效地提取一個(gè)數(shù)據(jù)序列的趨勢(shì)或去掉該數(shù)據(jù)序列的均值。測(cè)試結(jié)果表明，EMD方法是目前提取數(shù)據(jù)序列趨勢(shì)或均值的最好方法[4]，EMD方法的另一目的是為了進(jìn)一步對(duì)各IMF分量進(jìn)行Hilbert變換，獲得信號(hào)的瞬時(shí)頻率。

通過(guò)對(duì)信號(hào)進(jìn)行平穩(wěn)化處理，可以分解出n個(gè)IMF分量ci(t)和趨勢(shì)項(xiàng)rn(t)。即

式中，每個(gè)IMF分量都是平穩(wěn)信號(hào)。

EMD方法（平穩(wěn)化過(guò)程）基本思想是：假如一個(gè)原始信號(hào)序列的極大值或極小值數(shù)目比上跨零點(diǎn)（或下跨零點(diǎn)）的數(shù)目多出兩個(gè)（或兩個(gè)以上），則該數(shù)據(jù)就需要進(jìn)行平穩(wěn)化處理。具體方法是：找出原始數(shù)據(jù)序列的所有極大值點(diǎn)并將其用樣條插值函數(shù)插值成為原始數(shù)據(jù)序列的上包絡(luò)線；找出原始數(shù)據(jù)序列的所有極小值點(diǎn)并用樣條插值函數(shù)插值成為原始數(shù)據(jù)序列的下包絡(luò)線；上、下包絡(luò)線的均值為原始數(shù)據(jù)序列的平均包絡(luò)線；將原始數(shù)據(jù)序列減去該平均包絡(luò)后即可得到一個(gè)去掉高頻的新數(shù)據(jù)序列。

在頻域上，信號(hào)經(jīng)EMD 分解的各個(gè)IMF的瞬時(shí)頻率之間有如下關(guān)系：第1個(gè)IMF 含最高瞬時(shí)頻率成份，第i（其中i≥2）個(gè)IMF的瞬時(shí)頻率幾乎處處是第i+1個(gè)IMF的瞬時(shí)頻率的兩倍[4]。因此，信號(hào)經(jīng)EMD 分解出的各個(gè)IMF可以看做是對(duì)原信號(hào)進(jìn)行帶通濾波的結(jié)果[4-8]。通過(guò)EMD方法可實(shí)現(xiàn)時(shí)空尺度濾波，即對(duì)一個(gè)能分解出n個(gè)IMF的信號(hào)其低通濾波為：

由于噪聲頻率一般較高，而各個(gè)IMF的頻率幾乎是按2的負(fù)冪次方的形式遞減[5]，所以各個(gè)IMF所包含的噪聲其強(qiáng)度也越來(lái)越弱。因此，低頻IMF幾乎就是期望信號(hào)的低頻分量。

2 EMD-EKF方法

考慮如下非線性模型

式中：假定過(guò)程噪聲和量測(cè)噪聲均為加性零均值白噪聲。

擴(kuò)展卡爾曼濾波的方法一般流程如下[1-3，9-10]：

預(yù)測(cè)：

狀態(tài)的一步預(yù)測(cè)

協(xié)方差的一步預(yù)測(cè)

更新：

量測(cè)預(yù)測(cè)值

協(xié)方差

增益

狀態(tài)更新方程

協(xié)方差更新方程

其中，I為與協(xié)方差同維的單位矩陣。

經(jīng)驗(yàn)?zāi)B(tài)分解方法具體步驟如下[7-8，11-12]：

1）計(jì)算出原始信號(hào)Z(t) 所有的極值點(diǎn)。

2）采用三次樣條插值算法，求解所有的極大值點(diǎn)構(gòu)成的上包絡(luò)線和所有的極小值點(diǎn)構(gòu)成的下包絡(luò)線，分別記為u1(t)和 v1(t)；在這一步驟的計(jì)算中，為了抑制邊緣效應(yīng)，需使用鏡像延拓的方法來(lái)延拓極值點(diǎn)。

3）記上、下包絡(luò)線的均值為：

并記信號(hào)與上、下包絡(luò)線的均值之差為：

4）判斷 h1(1)(t)是否滿足IMF的兩條性質(zhì)。若滿足，則 h1(1)(t)為IMF。否則，記 h1(1)(t)為Z (t)，重復(fù)步驟1）至步驟4），假定k次篩選后得到的h1(1)(t)滿足IMF定義，則Z (t)的第一階IMF分量c1(t)=h1(1)(t)。

5）記 Z1(t)=Z(t) ? c1(t)為新的待分析信號(hào)。重復(fù)步驟1）至步驟5），以得到第二個(gè)IMF，記為c2(t)，余項(xiàng) Z2(t)=Z1(t) ?c2(t)。

重復(fù)上述步驟，直至滿足停止條件，分解結(jié)束。如此，可得到余項(xiàng)Zm(t)。

經(jīng)驗(yàn)?zāi)B(tài)分解方法能夠分離信號(hào)和噪聲，因此可以用來(lái)估計(jì)含噪聲信號(hào)中噪聲的標(biāo)準(zhǔn)差。在觀測(cè)噪聲未知的情況下，可以首先選取一段長(zhǎng)度為L(zhǎng)的一段信號(hào)，即相當(dāng)于一個(gè)觀測(cè)區(qū)間，在這個(gè)觀測(cè)區(qū)間上進(jìn)行經(jīng)驗(yàn)?zāi)B(tài)分解，分離出其中的高頻噪聲，估算出觀測(cè)噪聲的標(biāo)準(zhǔn)差。根據(jù)這個(gè)標(biāo)準(zhǔn)差結(jié)合EMD 濾波，進(jìn)行狀態(tài)估計(jì)。

本文的方法與常用的EKF方法的主要區(qū)別在于量測(cè)噪聲協(xié)方差矩陣不是預(yù)先根據(jù)經(jīng)驗(yàn)值確定的，而是根據(jù)實(shí)時(shí)觀測(cè)信號(hào)值估計(jì)出來(lái)的。這就保證了計(jì)算用的量測(cè)噪聲協(xié)方差矩陣和實(shí)際系統(tǒng)觀測(cè)噪聲協(xié)方差基本一致，即使實(shí)際系統(tǒng)由于受到某種干擾，噪聲協(xié)方差發(fā)生了變化，該方法也能計(jì)算出量測(cè)噪聲的參數(shù)。

3 仿真結(jié)果

無(wú)源定位是EKF 濾波的一個(gè)重要應(yīng)用領(lǐng)域，下面就單站無(wú)源定位中的跟蹤定位問(wèn)題進(jìn)行仿真。

在實(shí)際的系統(tǒng)中，各個(gè)觀測(cè)量隨時(shí)都有可能受到各種外界干擾，導(dǎo)致觀測(cè)噪聲發(fā)生變化，假定運(yùn)動(dòng)單站通過(guò)方位角進(jìn)行無(wú)源定位，觀測(cè)時(shí)間間隔T=2s，連續(xù)觀測(cè)100 s。傳感器初始位置為(?5 000,0)，在y軸上做勻速直線運(yùn)動(dòng)v=4m/s，在x軸和y軸方向上的擾動(dòng) rx=0.1、r y=0.1。目標(biāo)初始位置為(0,0)，在y軸上做勻速直線運(yùn)動(dòng)v=3m/s，系統(tǒng)噪聲Q=16。定義相對(duì)位置誤差（Relative Position Error，RPE）來(lái)描述方法的性能：

式中：(xtrue,ytrue)為真實(shí)相對(duì)位置；(x,y)為相對(duì)位置的估計(jì)值。

圖1是在變?cè)肼晽l件下進(jìn)行了仿真，觀測(cè)噪聲在0～50 s為3 mrad，在50～100 s為5 mrad，圖1共對(duì)兩種情況進(jìn)行了仿真。情況一為本文方法，即首先用EMD 估計(jì)量測(cè)噪聲協(xié)方差，再進(jìn)行EKF 濾波，情況二為采用通過(guò)先驗(yàn)信息確定的量測(cè)噪聲方差進(jìn)行EKF 濾波。

圖1 量測(cè)噪聲為0～50 s為3 mrad，50～100 s為5 mrad

從圖1和表1中可以看出本文方法的濾波效果明顯好于使用經(jīng)驗(yàn)值的方法，這是因?yàn)椴捎肊MD方法可以實(shí)時(shí)確定量測(cè)噪聲方差，因此在此基礎(chǔ)上進(jìn)行的濾波效果更好。其中，進(jìn)行200次蒙特卡羅仿真，使用本文方法耗時(shí)15.64 s，使用經(jīng)驗(yàn)值方法耗時(shí)1.34 s。

表1 EMD-EKF方法和EKF方法平均相對(duì)誤差的比較

4 結(jié)束語(yǔ)

不準(zhǔn)確的量測(cè)噪聲方差統(tǒng)計(jì)特性往往會(huì)導(dǎo)致EKF 濾波性能惡化。本文針對(duì)這一問(wèn)題，采用EMD方法實(shí)時(shí)估計(jì)觀測(cè)噪聲序列的噪聲標(biāo)準(zhǔn)差，并利用估計(jì)量測(cè)噪聲序列的噪聲標(biāo)準(zhǔn)差進(jìn)行EKF 濾波。仿真結(jié)果表明，本文方法可以實(shí)現(xiàn)測(cè)量量測(cè)序列噪聲的變化，準(zhǔn)確估計(jì)出量測(cè)噪聲標(biāo)準(zhǔn)差，從而避免了因量測(cè)噪聲的不準(zhǔn)確而引起的EKF 性能下降。

[1]何友,修建娟,張晶煒,等.雷達(dá)數(shù)據(jù)處理及應(yīng)用[M].北京:電子工業(yè)出版社,2008:42-46.

[2]PEARSON J B,STEAR E B.Kalman filter applications in airborne radar tracking[J].IEEE Transactions on Aerospace and Electronic systems,1972(10):319-329.

[3]JOSEPH J,LAVILOLA JR.A comparison of unscented and extended kalman filtering for estimating quaternion motion[C]//The proceeding of the 2003 American Control Conference.2003:2435-2440.

[4]鄧擁軍,王偉,錢成春,等.EMD方法及Hilbert變換中邊界問(wèn)題的處理[J].科學(xué)通報(bào),2001,46(3):257-263.

[5]PATRICK FLANDRIN,GABRIEL BILLING,PAULO GONCALVES.Empirical mode decomposition as a filter bank[J].IEEE Signal Processing Letters,2004,11(2):112-114.

[6]陳克安,馬遠(yuǎn)良.自適應(yīng)有源噪聲控制——原理、算法及實(shí)現(xiàn)[M].西安:西北工業(yè)大學(xué)出版社,1993:133-140.

[7]楊世錫,胡勁松,吳昭同,等.旋轉(zhuǎn)機(jī)械振動(dòng)信號(hào)基于EMD的希爾伯特變換和小波變換時(shí)頻分析比較[J].中國(guó)電機(jī)工程學(xué)報(bào),2003,23(6):119-121.

[8]譚善文,秦樹人,湯寶平.Hilbert-Huang變換的濾波特性及其應(yīng)用[J].重慶大學(xué)學(xué)報(bào),2004,27(2):9-12.

[9]管旭軍,芮國(guó)勝.基于UKF的單站無(wú)源定位算法[J].電光與控制,2004,11(1):34-36.

[10]芮國(guó)勝,管旭軍.一種機(jī)動(dòng)目標(biāo)無(wú)源定位的新方法[J].宇航學(xué)報(bào),2005,26（增刊）:121-125.

[11]劉慧婷,張旻,程家興.基于多項(xiàng)式擬合算法的EMD端點(diǎn)問(wèn)題的處理[J].計(jì)算機(jī)工程與應(yīng)用,2004,40(16):84-86.

[12]杜陳艷,張榆鋒,楊平,等.經(jīng)驗(yàn)?zāi)B(tài)分解邊緣效應(yīng)抑制方法綜述[J].儀器儀表學(xué)報(bào),2009,30(1):55-59.