基于卡爾曼和最小二乘的抗野值降噪方法研究

朱紅運(yùn)，龐建國，先治文

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

1 引言

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

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

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

2 卡爾曼濾波算法

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

首先，使用系統(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時(shí)刻系統(tǒng)的狀態(tài)；A、B分別為系統(tǒng)狀態(tài)轉(zhuǎn)移方程和控制方程；w(k)為協(xié)方差為Q的系統(tǒng)噪聲；Z(k)為k時(shí)刻系統(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時(shí)刻對(duì)k時(shí)刻系統(tǒng)狀態(tài)的預(yù)測；P為誤差協(xié)方差，P(k|k-1)為k-1時(shí)刻對(duì)k時(shí)刻誤差協(xié)方差的預(yù)測；AT為A的轉(zhuǎn)置。

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

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時(shí)刻卡爾曼增益，e(k)為k時(shí)刻殘差序列，X(k|k)為濾波后的狀態(tài)值，至此完成一次完整的卡爾曼濾波。

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

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

(9)

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

在野值存在時(shí)，殘差序列可表示為

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

(10)

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

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

3.1 野值判別

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

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

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

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

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

(11)

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

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

(12)

式中

(13)

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

(14)

式中

(15)

(16)

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

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

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

(17)

而后將其代入卡爾曼濾波器進(jìn)行濾波計(jì)算，完成單個(gè)野值的剔除。

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

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

(18)

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

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

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

步驟2：采用萊特準(zhǔn)則對(duì)k時(shí)刻殘差進(jìn)行判別，判斷該時(shí)刻觀測值是否為野值；

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

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

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

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

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

4 仿真與分析

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

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

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

為進(jìn)一步更好的比較各算法的性能，引入噪聲抑制率(noise suppression ration,NSR)和信號(hào)失真率(signal distortion ration,SDR)兩個(gè)指標(biāo)參數(shù)進(jìn)行評(píng)估。兩指標(biāo)表達(dá)式分別為

(19)

(20)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5 結(jié)論

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