譚平，蔡自興

(中南大學(xué) 信息科學(xué)與工程學(xué)院，湖南 長沙，410083)

卡爾曼濾波是最常用的最優(yōu)狀態(tài)估計方法，在滿足線性和標準噪聲(零均值、高斯分布)的條件下，通過遞推形式得到最小二乘意義下的最優(yōu)狀態(tài)估計[1]。而基于小波多尺度分解的正交變換能夠把單一尺度的狀態(tài)估計轉(zhuǎn)化到不同尺度上的狀態(tài)估計，得到比原來單一尺度上好的處理效果[2]，特別針對 1/f類分形信號，小波變換可以使得不同尺度下的小波系數(shù)序列的自相關(guān)和互相關(guān)函數(shù)衰減很快[3-4]。為此，研究者們將小波多尺度分析和卡爾曼方法結(jié)合，提出了多尺度卡爾曼濾波方法。多尺度卡爾曼濾波能充分發(fā)揮小波變換與卡爾曼濾波的優(yōu)點，取長補短，在信號去噪[5-6]、信號檢測[7]、信息融合[8-9]、目標跟蹤等方面效果明顯。現(xiàn)有的多尺度卡爾曼濾波通常采用分段方式，先對狀態(tài)序列進行小波多尺度分解，接著在各個尺度上進行卡爾曼濾波，然后通過小波重構(gòu)得到狀態(tài)估計值。但是，這種方式不滿足在線實時處理，在進行多尺度分解時需要得到t時刻后的數(shù)據(jù)才能進行分解[10-11]，即在t時刻，不能同時得到不同尺度上的信息，而且存在較大的延時，特別是分解尺度越多，其延時越大。為了實時準確的得到各個尺度上的信息，本文作者提出基于無抽取Haar算法[12]的實時卡爾曼濾波，該方法能在t時刻獲得各個尺度信息，降低了實時響應(yīng)時間，另結(jié)合各小波閾值去噪方法，在信號處理方面表現(xiàn)出良好去噪性能。

1 多尺度卡爾曼濾波方法

1.1 多尺度動態(tài)系統(tǒng)建模

當(dāng)采樣時刻為t(t≥0)，在某一尺度j上建立動態(tài)系統(tǒng)：

系統(tǒng)隨機噪聲序列和觀測隨機噪聲序列滿足：

1.2 多尺度卡爾曼濾波方法

多尺度卡爾曼濾波方法(Multi-scale Kalman Filter，MKF)是在通過小波變換后，在各個尺度上對信號進行卡爾曼濾波[13]，圖1所示為多尺度卡爾曼濾波的結(jié)構(gòu)，其實現(xiàn)過程如圖2所示。在計算時刻t的多尺度信號時，需要使用(t-2J-1)到(t+2J-1-1)(J為最大尺度)時刻的數(shù)據(jù)，因此其信號處理延時了2J-1-1個時間步長。當(dāng)實時性要求較高時，將無法滿足系統(tǒng)要求。

圖1 多尺度Kalman濾波結(jié)構(gòu)Fig.1 Structure of multi-scale Kalman filter

從圖2可知這種方法存在以下2 個缺點：(1) 在t時刻計算得到的小波系數(shù)，不能為下次小波系數(shù)的計算所使用，需要重復(fù)計算，效率低；(2) 由于小波的對稱性，為了計算t 時刻小波系數(shù)，需要延時2J-1-1個時間步長。圖2 中延時為7 個時間步長，因此，實時性較差。

圖2 t時刻小波變換中小波系數(shù)所使用數(shù)據(jù)示意圖Fig.2 Input signal used to calculate wavelet coefficient at t time step

2 多尺度實時卡爾曼濾波算法

2.1 多孔算法

考慮時間序列信號{c0,t}，可定義為在t時刻函數(shù)f(x)與尺度函數(shù)φ(x)的內(nèi)積，即

尺度函數(shù)滿足方程：

這里h是與尺度函數(shù)相關(guān)的離散低通濾波。由定義可知：信號通過低通濾波后與相鄰尺度信號非常接近，而2層之間的距離通過因子2從一層增加到另一層。

在時刻t和尺度j下，逼近系數(shù){cj,t}是函數(shù)f(x)與尺度函數(shù)內(nèi)積：

它也可以通過卷積計算得到：

因此，連續(xù)2個層之間的差可以表示為：

也可獨立表達為：

這里，小波函數(shù)定義為：

原始信號的重構(gòu)可以表示為：

2.2 基于無抽取Haar算法的Kalman濾波

為了減少邊界影響，同時提高變換的實時性，Renaud等[12]建議選用濾波器系數(shù)盡量窄的小波函數(shù)。并針對局部階躍效應(yīng)提出了無抽取Haar算法，該算法使用Haar小波基的低通濾波器h(k)={1/2,1/2}替換標準多孔算法，這里h是非對稱的，考慮到第1層小波分解，可以使用原始信號與h進行卷積，則：

通過式(10)和(11)變換可知：在任意時刻 t，都不需要使用t時刻以后的數(shù)據(jù)來計算小波系數(shù)，因此，增強了實時性。顯然該算法由于不需要平移信號，因此，信號＜x1，x2，…，xt＞在尺度j下的小波系數(shù)嚴格等價于信號＜x1，x2，…，xN＞的前t個小波系數(shù)，其中N＞t。該方法在實際中計算容易，只需要加減運算和右移運算即可完成，而且由于數(shù)據(jù)是進行有規(guī)律地更新，前面計算得到的小波系數(shù)可以為后面的計算使用，不需要重新計算，提高計算效率。

圖 3所示為計算各個尺度逼近系數(shù)所用到的數(shù)據(jù)，對應(yīng)的小波系數(shù)可以通過式(11)所得。

結(jié)合上面所提的快速計算方法，采用卡爾曼方法對各個尺度信息進行濾波，能有效提高收斂速度。為了提高精度，在小波分解后，在各個尺度上進行小波去噪，采用軟閾值方法[14]，小波系數(shù)的去噪公式如下：

圖3 t時刻小波變換中系數(shù)c所使用數(shù)據(jù)示意圖Fig.3 Input signal are used to calculate approximation coefficient c at t time step in wavelet transform

圖4 所示為基于無抽取Haar算法的實時多尺度卡爾曼濾波算法(Real time multiscale Kalman filter，RMKF)結(jié)構(gòu)圖。該算法的處理步驟如下(輸入為t時刻采樣數(shù)據(jù)xt；輸出為t時刻濾波數(shù)據(jù)yt)：(1) 采用式(10)和(11)進行小波分解，令 c0,t=xt；(2) 利用式(13)在各個尺度上進行閾值去噪；(3) 在各個尺度上進行Kalman濾波；(4) 利用式(9)進行小波重構(gòu)后得到濾波信號。

圖4 多尺度實時卡爾曼濾波算法結(jié)構(gòu)圖Fig.4 Structure of multi-scale real time Kalman filter

顯然，該算法能夠在t時刻處理完數(shù)據(jù)，具有較強的實時性，其計算復(fù)雜度為 S=J×(4+Skf)(其中，J為最大尺度；Skf為卡爾曼濾波計算復(fù)雜度)。

3 實驗結(jié)果及討論

3.1 傳感器去噪性能比較

為了驗證該方法的效果，通過三軸低精度加速度和一個精度較高的 GPS/INS中慣性測量單元(Inertial Measurement Unit, IMU)進行驗證，其中高精度數(shù)據(jù)作為參照對象。實驗在中南大學(xué)自主改裝的智能獵豹越野車上進行。

首先實時采集數(shù)據(jù)，然后分別采用小波閾值、普通卡爾曼濾波，MKF和MRKF 4種方法對低精度數(shù)據(jù)進行去噪處理，將結(jié)果與高精度的IMU加速度數(shù)據(jù)進行比較，這里采用均方誤差(Mean square error，MSE)對性能進行評價：

其中：xd為算法處理結(jié)果，xh為高精度傳感器采集數(shù)據(jù)。圖6所示為在設(shè)定狀態(tài)轉(zhuǎn)移誤差為q=0.02時，4種算法對x軸數(shù)據(jù)的處理結(jié)果。為了驗證算法的有效性，實驗重復(fù)20次，結(jié)果取平均值，如表1所示。

圖6 4種算法對x方向加速度的處理結(jié)果Fig.6 Denoise result of four algorithms to x acceleration

表1 在q=0.02時，4種算法的均方誤差Table 1 MSEs of four algorithms when q=0.02

從圖 6可以看出：(1) 小波閾值去噪處理算法對勻速情況下能有較好的處理，但是在加速度變換特別是智能車顛簸時，處理后數(shù)據(jù)噪聲依然很大，與實際數(shù)據(jù)依然相差很大。(2) 在狀態(tài)轉(zhuǎn)移誤差估計準確時，Kalman濾波、MKF和MRKF對大噪聲污染數(shù)據(jù)都有較好的處理能力，處理后的數(shù)據(jù)基本接近高精度數(shù)據(jù)。

為了驗證實時性能，在采樣速率為20 Hz，共有4個尺度的情況下，統(tǒng)計4種算法的平均延時，結(jié)果如表2所示?？梢姡篕alman延時最小，MRKF其次，小波去噪和MKF濾波時間較長，其原因是這2種算法本身需要等待2J-1-1個時間步長。

表2 4種算法平均延時Table 2 Average delay of four algorithms s

3.2 狀態(tài)轉(zhuǎn)移誤差q對濾波結(jié)果的影響

對于一個動態(tài)系統(tǒng)，其系統(tǒng)狀態(tài)轉(zhuǎn)移誤差通常比較難以估計，為了進一步比較Kalman濾波和RMKF的去噪能力，對系統(tǒng)方程中的狀態(tài)轉(zhuǎn)移誤差q進行分析。設(shè)定狀態(tài)轉(zhuǎn)移誤差q從0.001到0.700變化，對每個不同的q計算2種方法的均方誤差，為了保證計算的統(tǒng)計特性，實驗重復(fù)20次，取平均值得到結(jié)果如圖7所示。由圖7可知：(1) 存在最小估計誤差，即誤差估計有下限；(2) 當(dāng)狀態(tài)轉(zhuǎn)移誤差模型準確時，卡爾曼濾波是最優(yōu)估計，q=0.01時，其均方誤差最??；而q小于該值后，估計誤差均方誤差將顯著增大，因此，實際中不能將狀態(tài)轉(zhuǎn)移誤差設(shè)得太小。(3) 當(dāng)狀態(tài)誤差模型不精確且q＞0.01時，多尺度卡爾曼濾波后的均方誤差比普通Kalman濾波的小，其估計更精確。

圖7 3種算法的均方誤差隨系統(tǒng)狀態(tài)誤差估計變換Fig.7 MSEs of three algorithms with system state error estimated changes

因此，對于一個狀態(tài)轉(zhuǎn)移誤差不能確定的系統(tǒng)，采用RMKF方法能有效提高系統(tǒng)的性能。

4 結(jié)論

(1) 提出了基于無抽取Haar算法的卡爾曼濾波實時處理方法，通過采用簡單的加減、移位運算在t時刻完成多尺度變換，克服了基于分段式方法計算量大、實時性差的問題。該方法提高了算法實時性，并且有效減少重復(fù)運算，提高運算效率。

(2) 該方法結(jié)合小波閾值去噪，對強背景噪聲數(shù)據(jù)具有較強去噪性能，并通過實驗得到了驗證。

[1] 徐寧壽. 隨機信號估計與系統(tǒng)控制[M]. 北京:北京工業(yè)大學(xué)出版社, 2001: 135-141.XU Ning-shou. Random signal estimation and system control[M]. Beijing: Beijing Industrial University Press, 2001:135-141.

[2] 文成林. 多尺度估計理論及其應(yīng)用[M]. 北京: 清華大學(xué)出版社, 2002: 115-121.WEN Cheng-lin. Multiscale estimation theory and its application[M]. Beijing: Tsinghua University Press, 2002:115-121.

[3] CHEN Bor-sen, HOU Wen-sheng. Deconvolution filter design for fractal signal transmission systems: A multi-scale Kalman filter bank approach[J]. IEEE Transactions on Signal Processing,1997, 45(5): 1359-1364.

[4] Hirchuren G A, Attellis C E. Estimation of fractal Brownian motion with multiresolution Kalman filter banks[J]. IEEE Transactions on Signal Processing, 1999, 47(5): 1431-1433.

[5] 柯熙政, 任亞飛. 多尺度多傳感器融合算法在微機電陀螺數(shù)據(jù)處理中的應(yīng)用[J]. 兵工學(xué)報, 2009, 30(7): 994-998.KE Xi-zheng, REN Ya-fei. The application of multi-scale sensor fusion algorithm to MEMS gyroscope data processing[J]. Acta Armamentarii, 2009, 30(7): 994-998.

[6] 覃方君, 許江寧, 等. 基于小波卡爾曼濾波的加速度計降噪方法[J]. 武漢理工大學(xué)學(xué)報: 交通科學(xué)與工程版, 2009, 33(1):49-52.QIN Fang-jun, XU Jiang-ning. A wavelet based kalman filter for accelerometers de-noising[J]. Journal of Wuhan University of Technology: Transportation Science & Engineering, 2009, 33(1):49-52.

[7] 尹世榮, 王蔚然. 多尺度卡爾曼濾波探測氣體濃度[J]. 激光與紅外, 2005, 35(9): 679-681.YIN Shi-rong, WANG Wei-ran. Multiscale kalman filtering gas concentration detecting[J]. Laser & Infrared, 2005, 35(9):679-681.

[8] ZHAO Tong-zhou, WANG Yan-li, WANG Hai-hui. Image Fusion Based on Multi-scale Kalman Filtering[C]//2010 Second International Workshop on Education Technology and Computer Science(ETCS). Wuhan, 2010, 3: 207-215.

[9] 崔培玲, 王桂增, 潘泉. 一類動態(tài)多尺度系統(tǒng)融合估計算法的分析[J]. 控制理論與應(yīng)用, 2007, 24(1): 84-89.CUI Pei-ling, WANG Gui-zeng, PAN Quan. Analyzing of the fusion estimation algorithm of a class of dynamic multiscale system[J]. Control Theory & Applications, 2007, 24(1): 84-89.

[10] 時偉, 吳美平, 薛祖瑞. 基于小波卡爾曼混合濾波的激光陀螺信號處理[J]. 兵工自動化, 2005, 24(4): 57-58.SHI Wei, WU Mei-ping, XUE Zu-rui. Signal processing for laser gyro based on wavelet-kalman mixed filtering[J]. Ordnance Industry Automation, 2005, 24(4): 57-58.

[11] 丁寧, 周新志. 基于改進多孔算法的時間序列預(yù)測[J]. 系統(tǒng)仿真學(xué)報, 2007, 19(17): 4082-4085.DING Ning, ZHOU Xin-zhi. Time series forecasting based on improved á trous algorithm [J]. Journal of System Simulation,2007, 19(17): 4082-4085.

[12] Renaud O, Starck J L, Murtagh F. Wavelet-based combined signal filtering and prediction[J]. IEEE Transactions on Systems,Man, and Cybernetics, 2005, 36(6): 1241-1251.

[13] 李強, 王其申. 1/f分形噪聲的一種多尺度Kalman濾波方法[J].量子電子學(xué)報, 2007, 24(1): 7-12.LI Qiang, WANG Qi-shen. 1/f fractal noise attenuation using multiscale Kalman filter[J]. Chinese Journal of Quantum Electronics, 2007, 24(1): 7-12.

[14] Donoho D L. Denoising by soft-thresholding[J]. IEEE Transactions on Information Theory, 1995, 41(3): 613-627.