朱紅運，龐建國，先治文

(中國人民解放軍63729部隊，山西 太原 030000)

1 引言

在工程測試中，由于測量設(shè)備本身誤差、隨機誤差和各種環(huán)境因素等干擾，使得測量數(shù)據(jù)存在噪聲，部分測量數(shù)據(jù)還會包含有嚴重偏離目標真實值的數(shù)據(jù)，通常將這些嚴重偏離目標真實值的數(shù)據(jù)稱為“野值”。根據(jù)是否連續(xù)，野值可分為孤立型野值和連續(xù)型野值[1]，孤立型野值表現(xiàn)形式是孤立的點，而連續(xù)型野值則成片出現(xiàn)，也稱為斑點型野值。在對測量數(shù)據(jù)降噪時，野值會給降噪結(jié)果帶來很大的誤差，甚至會使信號嚴重失真。因此，為獲取有效測量數(shù)據(jù)，需采用魯棒性強、濾波精度高的降噪方法，對原始數(shù)據(jù)進行降噪，并將野值剔除或進行必要的修正。

卡爾曼濾波作為線性高斯系統(tǒng)的最優(yōu)濾波算法，具有理論基礎(chǔ)完備、計算簡便等優(yōu)點，已廣泛應(yīng)用于測量數(shù)據(jù)的降噪處理[2-3]。傳統(tǒng)的卡爾曼濾波魯棒性相對較差[4]，為進一步提高其魯棒性和濾波精度，許多學者對卡爾曼濾波算法進行了改進。文獻[5-7]在傳統(tǒng)卡爾曼濾波算法的基礎(chǔ)上，通過引入調(diào)節(jié)因子，設(shè)計了一種自適應(yīng)卡爾曼濾波算法，提高了卡爾曼濾波器的抗野值干擾能力，但上述方法在選取調(diào)節(jié)因子時還需要依靠經(jīng)驗；文獻[8]提出了一種能動態(tài)調(diào)整抗野值門限的自適應(yīng)卡爾曼濾波算法，該方法可根據(jù)信號變化率自適應(yīng)調(diào)整抗野值門限，有效降低野值的影響；當野值連續(xù)出現(xiàn)時，雖然該方法門限能夠自適應(yīng)調(diào)整，具有一定的抗連續(xù)型野值干擾能力，但由于相鄰野值的變化率是隨機的，在抗野值門限調(diào)整時必然會存在一定的誤差；文獻[9]提出了一種基于最小二乘和卡爾曼濾波的跟蹤定位算法，其首先采用最小二乘法對目標點進行估算，而后使用卡爾曼濾波算法對估算值進一步修正，實現(xiàn)了高精度跟蹤定位，但該方法需要對每個觀測值進行擬合，隨著數(shù)據(jù)增多，擬合誤差會逐漸增大，對連續(xù)型野值的抗干擾能力有限。

鑒于此，為有效降低野值對濾波的影響，提高卡爾曼濾波精度，本文提出了一種基于卡爾曼和最小二乘的抗野值降噪方法，該方法在進行卡爾曼濾波過程中，采用最小二乘擬合方法對野值進行修正，并通過引入調(diào)節(jié)因子降低連續(xù)型野值擬合誤差的影響，提高其抗連續(xù)型野值干擾能力，最后通過仿真驗證了該方法的有效性。

2 卡爾曼濾波算法

卡爾曼濾波算法基本原理是求動態(tài)系統(tǒng)狀態(tài)序列的線性最小方差誤差估計，利用動態(tài)的狀態(tài)方程和觀測方程來描述系統(tǒng)。其在濾波過程中，首先建立描述隨機動態(tài)變量隨時間變化的先驗?zāi)Ｐ停笤趯崟r觀測隨機變量基礎(chǔ)上，根據(jù)前一次估計值即可推導(dǎo)出當前目標狀態(tài)的最優(yōu)估計值。

首先，使用系統(tǒng)的狀態(tài)方程和觀測方程來描述需要測量的系統(tǒng)

X(k)=A·X(k-1)+B·U(k)+w(k)

(1)

Z(k)=H·X(k)+y(k)

(2)

式中，X(k)、X(k-1)分別為k與k-1時刻系統(tǒng)的狀態(tài)；A、B分別為系統(tǒng)狀態(tài)轉(zhuǎn)移方程和控制方程；w(k)為協(xié)方差為Q的系統(tǒng)噪聲；Z(k)為k時刻系統(tǒng)的觀測值，H為觀測狀態(tài)轉(zhuǎn)移矩陣，y(k)為協(xié)方差為R的觀測噪聲。

預(yù)測方程為

X(k|k-1)=AX(k-1|k-1)+BU(k)

(3)

P(k|k-1)=AP(k-1|k-1)AT+Q

(4)

式中，X(k|k-1)為k-1時刻對k時刻系統(tǒng)狀態(tài)的預(yù)測；P為誤差協(xié)方差，P(k|k-1)為k-1時刻對k時刻誤差協(xié)方差的預(yù)測；AT為A的轉(zhuǎn)置。

而后對狀態(tài)方程進行修正，具體如下

Kg(k)=P(k|k-1)HT/(HP(k|k-1)HT+R)

(5)

e(k)=Z(k)-HX(k|k-1)

(6)

X(k|k)=X(k|k-1)+Kg(k)e(k)

(7)

P(k|k)=(I-Kg(k)H)P(k|k-1)

(8)

式中，Kg(k)為k時刻卡爾曼增益，e(k)為k時刻殘差序列，X(k|k)為濾波后的狀態(tài)值，至此完成一次完整的卡爾曼濾波。

當濾波達到穩(wěn)定狀態(tài)時，e(k)為零均值服從正態(tài)分布的隨機序列，此時濾波精度較高。假設(shè)k時刻前觀測值均正常，在k時刻出現(xiàn)野值，則此時觀測方程可表示為

Z′(k)=Z(k)+Δ

(9)

式中，Z(k)與Z′(k)分別為疊加野值前后的觀測值，Δ為野值分量。

在野值存在時，殘差序列可表示為

e′(k)=Z(k)-HX(k|k-1)+Δ

(10)

由式(7)可知，在求濾波后狀態(tài)值時，野值分量被放大為原來的Kg(k)倍引入到狀態(tài)估計中，特別是當野值連續(xù)出現(xiàn)時，殘差得不到及時修正，會嚴重影響濾波精度，因而，在使用卡爾曼濾波器進行濾波時，必須對野值進行剔除或修正。

3 卡爾曼濾波與最小二乘擬合的聯(lián)合應(yīng)用

3.1 野值判別

由于異常數(shù)據(jù)會在殘差中體現(xiàn)出來，因而通過判別殘差是否超過了合理范圍即可對野值進行辨識。對于一個服從正態(tài)分布的序列，根據(jù)高斯理論，其測量值落在[-3σ，3σ]以外的概率不超過0.3%(σ為序列的標準差)，可認為落在該區(qū)域以外的值為異常值，這就是萊特準則判別方法，也稱為3σ方法。因此可以以3σ為野值判斷閾值，殘差中落于[-3σ，3σ]以外的觀測值即可認為是野值。

當判出野值時，需對野值進行適當?shù)奶幚恚壳俺Ｓ玫姆椒ㄊ侵苯邮褂妙A(yù)測值代替野值[7]，但該方法未對預(yù)測值進行修正，從而降低了濾波精度，特別是當野值連續(xù)出現(xiàn)時，濾波效果會受到嚴重影響。為此，本文采用最小二乘擬合方法對野值進行修正，以提高濾波精度。

3.2 基于最小二乘擬合的野值修正

設(shè)存在n+1個點(xk，yk)，k=0，1，…，n，可求一個m次的最小二乘擬合多項式

Pm(x)=a0+a1x+a2x2+…+amxm

(11)

為求得上式，首先構(gòu)造滿足點(xk，yk)的正交基函數(shù){Qj(x)(j=0，1，…，m)}；其中多項式Qj(x)函數(shù)次數(shù)不大于m，則Pm(x)可表示為

Pm(x)=q0Q0(x)+q1Q1(x)+…+qmQm(x)

(12)

式中

(13)

此時，正交基函數(shù){Qj(x)(j=0，1，…，m)}的遞推表達式為

(14)

式中

(15)

(16)

將m個連續(xù)已知的點坐標代入到擬合多項式中，由最小二乘擬合原理得出擬合系數(shù)，求得擬合多項式，而后即可利用該擬合多項式外推出下一點的數(shù)值。

設(shè)k時刻系統(tǒng)的觀測值Z(k)為野值，X(k|k-1)為卡爾曼濾波器對k時刻系統(tǒng)狀態(tài)的預(yù)測值，X′(k)為外推擬合值，此時用擬合值代替野值，殘差可更新為

e(k)=X′(k)-HX(k|k-1)

(17)

而后將其代入卡爾曼濾波器進行濾波計算，完成單個野值的剔除。

受擬合多項式次數(shù)限制，當野值連續(xù)出現(xiàn)時，擬合誤差會逐漸積累，濾波準確度會隨之降低，為保證濾波精度，提高卡爾曼濾波抗連續(xù)型野值干擾能力，引入調(diào)節(jié)因子λ，此時，濾波后狀態(tài)X(k|k)可表示為

X(k|k)=X(k|k-1)+λKg(k)e(k)

(18)

由式(18)可知，隨著野值連續(xù)個數(shù)的增加，λ逐漸減小，即當野值連續(xù)出現(xiàn)時，外推擬合值對預(yù)測值的修正作用逐漸降低，由此可有效抑制擬合誤差帶來的影響，提高濾波精度。

綜上，本文提出的基于卡爾曼濾波和最小二乘擬合的抗野值降噪方法具體步驟可表述為：

步驟1：對原始信號數(shù)據(jù)進行卡爾曼濾波計算，求得k時刻殘差；

步驟2：采用萊特準則對k時刻殘差進行判別，判斷該時刻觀測值是否為野值；

步驟3：若步驟2中k時刻觀測值為正常值，則采用該觀測值對卡爾曼預(yù)測值進行修正，而后求得最終濾波結(jié)果；

步驟4：若步驟2中k時刻觀測值為野值，則采用最小二乘方法得出該時刻擬合值；

步驟5：使用步驟4中擬合值代替野值，更新k時刻殘差；

步驟6：設(shè)計調(diào)節(jié)因子λ；

步驟7：使用步驟6設(shè)計的調(diào)節(jié)因子對步驟5中更新后的殘差進行修正，而后由式(18)求得最終濾波結(jié)果。

4 仿真與分析

為檢驗所提出算法的性能，首先采用傳統(tǒng)卡爾曼濾波算法、文獻[8]算法、文獻[9]算法及本文算法分別對含孤立型野值的含噪仿真信號進行降噪，原始含噪信號及不同算法降噪結(jié)果如圖1所示。含噪信號信噪比為10dB，其表達式為y(t)=s(t)+n(t)+δ(t)；式中，y(t)為含噪信號，s(t)為幅值為1的原始正弦波信號，n(t)為服從正態(tài)分布的噪聲序列，δ(t)為孤立型野值，野值分別位于信號第300、500、700及1000個采樣點處。

圖1 不同算法對含孤立型野值含噪信號的降噪結(jié)果

由圖1可以看出，上述4種算法對噪聲均有一定的抑制作用，但傳統(tǒng)卡爾曼濾波算法抗野值干擾能力相對較弱，降噪后信號在野值點處存在一定的失真現(xiàn)象。

為進一步更好的比較各算法的性能，引入噪聲抑制率(noise suppression ration,NSR)和信號失真率(signal distortion ration,SDR)兩個指標參數(shù)進行評估。兩指標表達式分別為

(19)

(20)

式中，y′(n)為降噪后信號，s(n)為原始信號，y(n)為含噪信號。

由式(19)～式(20)可知，在采用上述指標評定時，噪聲抑制率越大、信號失真率越小，則表明算法的降噪性能越好。不同算法降噪后噪聲抑制率(NSR)和信號失真率(SDR)如表1所示。

表1 不同算法降噪后NSR與SDR

由表1可知，與傳統(tǒng)卡爾曼濾波算法相比，文獻[8]、文獻[9]算法對噪聲抑制能力更強，但由于文獻[9]在進行濾波時，用最小二乘法擬合值代替了觀測值，隨著數(shù)據(jù)增多，擬合誤差會逐漸增大，導(dǎo)致該方法降噪后信號失真率相比文獻[8]更大；本文算法在采用最小二乘擬合時，僅對野值進行了修正，能夠有效消除野值的影響，具有最好的降噪性能。

為檢驗所提算法抗連續(xù)型野值干擾能力，在含噪仿真信號中加入連續(xù)型野值，而后分別采用上述算法對含噪信號進行降噪處理，原始含噪信號及不同算法降噪結(jié)果如圖2所示。含噪信號信噪比為10dB，其表達式為y(t)=s(t)+n(t)+δ(t)；式中，y(t)為含噪信號，s(t)為幅值為1的原始正弦波信號，n(t)為服從正態(tài)分布的噪聲序列，δ(t)為連續(xù)型野值(連續(xù)6個野值)，連續(xù)型野值分別在信號第300、500、700及1000個采樣點處開始出現(xiàn)。

圖2 不同算法對含連續(xù)型野值含噪信號的降噪結(jié)果

由圖2可以看出，當野值連續(xù)出現(xiàn)時，傳統(tǒng)卡爾曼濾波算法無法有效降低野值的影響，降噪后信號中仍存在突變點，其余算法降噪后信號在野值對應(yīng)位置也出現(xiàn)了不同程度的失真。

不同算法降噪后噪聲抑制率(NSR)和信號失真率(SDR)如表2所示。

表2 不同算法降噪后NSR與SDR

由表2可知，文獻[8]、文獻[9]算法降噪性能均優(yōu)于傳統(tǒng)卡爾曼濾波算法，同樣由于文獻[9]在進行濾波時，用最小二乘法擬合值代替了觀測值，隨著數(shù)據(jù)增多，擬合誤差會逐漸增大，致使該方法降噪后信號失真率相比文獻[8]更大；當野值連續(xù)出現(xiàn)時，會對文獻[8]設(shè)計的自適應(yīng)抗野值門限造成一定的誤差，導(dǎo)致該方法的噪聲抑制率低于文獻[9]算法；在對連續(xù)型野值處理時，本文算法通過引入調(diào)節(jié)因子，有效降低了擬合誤差帶來的影響，提高了濾波精度，因而具有最強的抗連續(xù)型野值干擾能力。

為進一步驗證所提算法的魯棒性，采用該算法分別對信噪比為10dB、8dB、6dB的含孤立型和連續(xù)型野值的含噪信號進行降噪，結(jié)果分別如表3、表4所示。

表3 含孤立型野值含噪信號的降噪結(jié)果

表4 含連續(xù)型野值含噪信號的降噪結(jié)果

由表3、表4可知，無論野值為孤立型還是連續(xù)型，隨著信噪比減小，噪聲抑制率均逐漸增大，表明本文算法具有較強的魯棒性，在強噪聲影響下仍具有很好的降噪性能及抗野值干擾能力；此外，隨著噪聲強度增強，由于受噪聲影響，降噪后信號失真率緩慢增大，屬正?，F(xiàn)象，降噪結(jié)果表明，當噪聲較強時，信號失真率仍可以滿足要求。

5 結(jié)論

本文提出了一種基于卡爾曼和最小二乘的抗野值降噪方法，該方法在進行卡爾曼濾波過程中，使用最小二乘擬合方法對野值進行修正，有效降低了野值對降噪結(jié)果的影響。由實驗結(jié)果和理論分析可知，在進行卡爾曼濾波時，使用最小二乘擬合值代替測量信號中的野值可有效提高卡爾曼濾波算法性能，降低野值對降噪結(jié)果的影響；當野值連續(xù)出現(xiàn)時，逐漸減弱野值擬合結(jié)果對卡爾曼預(yù)測值的修正作用，可進一步提高卡爾曼濾波算法抗連續(xù)型野值的干擾能力。