基于模糊遞歸和最優(yōu)硬閾值的局部投影降噪算法

王 東, 崔忠馬, 陳文東, 舒 勤,*

(1. 四川大學(xué)電氣工程學(xué)院, 四川 成都 610065; 2. 北京遙感設(shè)備研究所, 北京 100084)

0 引 言

實(shí)際測得的信號(hào)易受噪聲干擾,從而破壞信號(hào)的結(jié)構(gòu)特征,因此降噪是后續(xù)各種信號(hào)處理的基礎(chǔ)[1-2]。傳統(tǒng)信號(hào)降噪方法基于線性系統(tǒng)的假設(shè),認(rèn)為信號(hào)和噪聲的頻譜不完全重合。但是，由混沌系統(tǒng)產(chǎn)生的信號(hào)往往具有寬頻譜、內(nèi)在偽隨機(jī)等混沌特性,信號(hào)和噪聲頻譜重合,使得傳統(tǒng)線性濾波方法失效[3]。

Cawley等[4]和Sauer等[5]提出的基于混沌動(dòng)力學(xué)理論的局部投影降噪算法為具有寬頻譜特性的混沌信號(hào)降噪提供了一種新方法。該方法首先根據(jù)Takens定理重構(gòu)信號(hào),得到一個(gè)與原始動(dòng)力系統(tǒng)微分同胚的相空間[6]。期望信號(hào)的吸引子被限制在一個(gè)低維流形上,而噪聲的吸引子分散在流形周圍。局部投影降噪算法根據(jù)信號(hào)與噪聲在相空間中局部動(dòng)力學(xué)特性與局部幾何特性的不同，區(qū)分信號(hào)子空間和噪聲子空間,利用幾何投影去除噪聲分量,再將相空間反重構(gòu)為時(shí)間信號(hào),從而達(dá)到降噪的目的。在實(shí)際情況中，許多信號(hào)都有混沌特性,目前該算法成功應(yīng)用的場景包括人類語音信號(hào)處理[7]、振動(dòng)信號(hào)分析[8]，還包括生物信號(hào)處理,如腦電信號(hào)處理[9]、心電信號(hào)處理[10]、呼吸聲音信號(hào)處理[11],以及激光數(shù)據(jù)處理[12]、故障檢測[13]等。

鄰域選取和子空間劃分是局部投影算法的兩個(gè)研究重點(diǎn)。鄰域選取的方法很多,但都存在不足。Cawley等[4]和Sauer等[5]直接指定鄰點(diǎn)數(shù)目,這種方式易受人為因素影響。Kantz等[12]利用遞歸圖估計(jì)鄰域半徑,效果優(yōu)于原始方法,但計(jì)算復(fù)雜。馮飛龍等[14]利用小波分解估計(jì)初始鄰域半徑,再進(jìn)行鄰點(diǎn)搜索,徐禮勝等[15]采用經(jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decomposition, EMD)估計(jì)鄰域半徑[15],但這些方法都需要預(yù)先估計(jì)噪聲水平。Przybya等[16]和Kotas等[17]利用K-means聚類算法確定鄰域,但聚類數(shù)難以確定。對于子空間劃分,現(xiàn)有方法直接指定子空間的維數(shù),或根據(jù)特征值的大小對子空間進(jìn)行劃分。Chelidze等[18]通過重構(gòu)信號(hào)的短時(shí)軌跡,利用平滑正交分解識(shí)別子空間。但考慮到混沌吸引子本身具有分?jǐn)?shù)維的性質(zhì)[19],系統(tǒng)局部鄰域的動(dòng)力學(xué)特性以及每個(gè)鄰域內(nèi)的噪聲分量占據(jù)的空間也不盡相同,且在實(shí)際情況中可能并不知道原始動(dòng)力系統(tǒng)的相空間維數(shù),所以每個(gè)局部鄰域應(yīng)該進(jìn)行不同的子空間劃分。

針對以上問題,本文提出基于模糊遞歸圖與最優(yōu)硬閾值準(zhǔn)則的局部投影降噪算法。首先，根據(jù)模糊遞歸圖對鄰域進(jìn)行選擇。為避免計(jì)算鄰域協(xié)方差矩陣,直接將鄰域矩陣進(jìn)行奇異值(singular value decomposition, SVD)分解；然后，根據(jù)最優(yōu)硬閾值對局部鄰域的信號(hào)子空間和噪聲子空間進(jìn)行劃分,避免人為因素的影響；最后，針對高斯白噪聲,采用本文所提方法分別對仿真Lorenz信號(hào)與實(shí)測含噪心電圖(electrocardiogram, ECG)信號(hào)進(jìn)行仿真研究,并與其他局部投影算法以及其他ECG信號(hào)降噪方法進(jìn)行對比,仿真結(jié)果驗(yàn)證了本文方法的有效性。

1 局部投影算法原理

1.1 相空間重構(gòu)

根據(jù)Takens嵌入定理,對于無限長、無噪的d維混沌吸引子的標(biāo)量時(shí)間信號(hào){x(t)},在拓?fù)洳蛔兊囊饬x下可以找到一個(gè)m維的嵌入相空間。其中,m≥2d+1。而對于有限長、含噪的信號(hào)，可采用坐標(biāo)延遲重構(gòu)[20]。

設(shè)x1,x2,…,xN為一長度為N的含高斯白噪聲的單變量時(shí)間信號(hào),選定一個(gè)時(shí)間延遲τ和嵌入維數(shù)m,構(gòu)造如下的相空間矢量:

Xi=[xi,xi+τ,…,xi+(m-1)τ]T,i=1,2,…,M

(1)

式中:M=N-(m-1)τ。τ由平均互信息量法[21]確定,即

(2)

式中：yi=xi+τ。取I(τ)第一個(gè)極小值點(diǎn)對應(yīng)的時(shí)間作為時(shí)間延遲。m由Cao[22]所提的方法確定,即

(3)

1.2 局部投影算法

基于相空間重構(gòu)局部投影降噪算法,利用信號(hào)與噪聲在相空間中軌線的動(dòng)力學(xué)特性與幾何特性的不同,保留信號(hào)分量,抑制噪聲分量,最大程度地恢復(fù)原始信號(hào)的吸引子流形。

對于相點(diǎn)Xn,將其在鄰域內(nèi)線性化展開:

(4)

(5)

(6)

第n個(gè)相點(diǎn)Xn的鄰域加權(quán)矩陣為

(7)

對式(7)進(jìn)行SVD分解,得到:

(8)

式中:Vn=[a1,a2,…,aK]由信號(hào)子空間與噪聲子空間構(gòu)成。將SVD分解得到的奇異值按從小到大進(jìn)行排列，有σ1≥σ2≥…≥σK-1≥σK,其中大的奇異值對應(yīng)信號(hào)子空間,小的奇異值對應(yīng)噪聲子空間。根據(jù)Vn得到由噪聲引起的分量,減去第n個(gè)相點(diǎn)在噪聲子空間中的投影,即得到修正后的相點(diǎn):

(9)

為了進(jìn)一步提高降噪效果,避免局部線性化產(chǎn)生較大誤差,本文采用Moore等[23]提出的質(zhì)心修正方法。修正后的質(zhì)心為

(10)

2 改進(jìn)的局部投影算法

2.1 鄰域選擇

在局部投影降噪算法中,鄰域的選擇十分重要。鄰域選擇得過小,會(huì)損失有效信息,受噪聲干擾嚴(yán)重;鄰域選擇得過大,會(huì)使得線性逼近的效果不好。

遞歸圖[24]可以用來提取時(shí)間信號(hào)中的相關(guān)信息。但是，遞歸圖的缺點(diǎn)是,動(dòng)力系統(tǒng)的遞歸模式可視化對相似閾值的選取十分敏感[25]。Pham等[26]在遞歸圖基礎(chǔ)上提出了模糊遞歸圖。相比于傳統(tǒng)遞歸圖,模糊遞歸圖不需要選取相似閾值,且模糊遞歸圖以灰度圖的形式展示,能夠?yàn)槟Ｊ椒治鎏峁└S富的信息。本文采用模糊遞歸圖來對鄰域進(jìn)行選擇。

(11)

式中：μ(Xi,Xj)∈[0,1]表示Xi和Xj之間的一種模糊相似性的度量,其具有如下性質(zhì)。

(1) 自反性:

μ(Xi,Xi)=1

(12)

(2) 對稱性:

μ(Xi,Vj)=μ(Vj,Xi)

(13)

(3) 傳遞性:

μ(Xi,Xj)=max[min{μ(Xi,Vk),μ(Xj,Vk)}]

(14)

式中：i=1,2,…,M;k=1,2,…，c；c為聚類數(shù),1

μ(Xi,Xj)的值根據(jù)模糊C均值(fuzzy C-Mean, FCM)聚類算法[27]計(jì)算得到。FCM算法通過最小化如下模糊目標(biāo)函數(shù)實(shí)現(xiàn):

(15)

式中：ω∈[1,+∞)為模糊度參數(shù);U=[μij](i=1,2,…,M;j=1,2,…,c)是劃分矩陣;V=(V1,V2,…,Vc)是聚類中心矩陣;Vj表示第j個(gè)聚類中心;d(Xi,Vj)表示某一范數(shù),本文采用歐式范數(shù)。上述模糊目標(biāo)函數(shù)滿足

(16)

為了得到最優(yōu)的U和V,通過迭代過程數(shù)值求解目標(biāo)函數(shù)的最小值,迭代過程如下:

(17)

(18)

當(dāng)滿足‖Ut-Ut+1‖<ε時(shí),停止迭代。其中,t為迭代次數(shù),ε為給定的精度水平。

模糊遞歸圖具有對稱性,可以看作是相空間狀態(tài)矢量之間的一種模糊關(guān)系。模糊遞歸圖用灰度圖表示,灰度值代表了狀態(tài)矢量對之間的模糊關(guān)系,這與遞歸圖是相互兼容的。灰度值越小,則表示這兩個(gè)狀態(tài)矢量越相似。一個(gè)灰度值為0的像素點(diǎn)代表了兩個(gè)狀態(tài)矢量完全相似,即代表動(dòng)力系統(tǒng)中一個(gè)100%的遞歸事件[28]。

2.2 信號(hào)與噪聲子空間劃分

在進(jìn)行投影修正之前,需要根據(jù)SVD分解得到的奇異值對信號(hào)子空間與噪聲子空間進(jìn)行劃分。本文利用Gavish等[29]提出的最優(yōu)硬閾值準(zhǔn)則對信號(hào)子空間與噪聲子空間進(jìn)行劃分,不需要預(yù)先估計(jì)噪聲水平。最優(yōu)硬閾值γ為

γ=ω(β)σmed

(19)

式中:σmed為奇異值的中位數(shù);ω(β)=λ(β)/μβ,λ(β)為

(20)

對于前述的鄰域加權(quán)矩陣Zn∈Rm×n;當(dāng)m=n時(shí),β=1;當(dāng)m>n時(shí),β=n/m;當(dāng)m

(21)

可以通過數(shù)值方法求解得到式(21)積分方程中μβ的值,從而求得閾值γ。

根據(jù)最優(yōu)硬閾值準(zhǔn)則,將大于等于閾值γ的奇異值所對應(yīng)的奇異向量所形成的子空間作為信號(hào)子空間,小于閾值γ的奇異值所對應(yīng)的奇異向量所張成的子空間作為噪聲子空間,由此實(shí)現(xiàn)了子空間的劃分。

2.3 改進(jìn)算法的基本步驟

基于以上分析,本文所提的局部投影降噪算法步驟如算法1所示。

算法1 改進(jìn)的局部投影降噪算法步驟輸入:含噪信號(hào)x(t)輸出:降噪信號(hào)x^(t)開始 相空間重構(gòu):利用式(2)計(jì)算平均互信息量并取第一個(gè)局部極小值對應(yīng)的時(shí)間作為時(shí)間延遲τ;利用式(3)計(jì)算E1(m),取停止變化時(shí)對應(yīng)的維數(shù)作為嵌入維數(shù)m;對x(t)進(jìn)行相空間重構(gòu),得到式(1)表示的相空間X。 循環(huán)(i=1∶M) 1. 選定參考相點(diǎn)Xi。 2. 鄰域選擇:根據(jù)式(11)計(jì)算重構(gòu)相空間的模糊遞歸圖,得到第i個(gè)參考相點(diǎn)的鄰域。 3. 計(jì)算鄰域質(zhì)心和鄰域矩陣:由式(10)和式(7)分別計(jì)算第i個(gè)參考相點(diǎn)的鄰域質(zhì)心和鄰域矩陣。 4. SVD分解:根據(jù)式(8)對鄰域矩陣進(jìn)行SVD分解,得到奇異值與右奇異向量。 5. 子空間劃分:根據(jù)式(19)的最優(yōu)硬閾值準(zhǔn)則對信號(hào)子空間與噪聲子空間進(jìn)行劃分。 6. 投影修正:根據(jù)式(9)進(jìn)行投影修正。 結(jié)束 反重構(gòu):將相空間恢復(fù)為時(shí)間信號(hào)x^(t),恢復(fù)方式采用式(22)[30],以減小誤差。結(jié)束

為減小反重構(gòu)所產(chǎn)生的誤差，進(jìn)行如下操作：

(22)

為達(dá)到較好的效果，需將以上步驟重復(fù)迭代幾次。

3 仿真實(shí)例與工程應(yīng)用

3.1 Lorenz混沌信號(hào)仿真

本文首先采用Lorenz混沌系統(tǒng)信號(hào)進(jìn)行仿真。Lorenz混沌系統(tǒng)由如下偏微分方程組描述：

(23)

當(dāng)參數(shù)取σ=10, r=28, b=8/3時(shí),系統(tǒng)呈現(xiàn)出混沌特性。采用四階Runge-Kutta方法數(shù)值求解上述偏微分方程組,初始值取x0=10, y0=10, z0=10,積分步長為0.01。取X分量中的2 000個(gè)數(shù)據(jù)點(diǎn)進(jìn)行仿真。對信號(hào)添加高斯白噪聲,添加噪聲后的信噪比為8 dB,然后使用本文所提的局部投影方法進(jìn)行降噪處理。降噪效果如圖1所示。

由圖1可知,本文方法具有較好的降噪效果,在噪聲較強(qiáng)的情況下也能夠恢復(fù)信號(hào)。由圖1可以看出，降噪后，信號(hào)比較光滑,看不到噪聲的影響,同時(shí)基本保持了Lorenz混沌吸引子的幾何形狀。

3.2 參數(shù)對算法的影響

本小節(jié)討論鄰域選擇和閾值這兩個(gè)參數(shù)對算法的影響。首先是鄰域選擇對算法的影響。Lorenz信號(hào)添加噪聲后的信噪比為8 dB。選取不同鄰域?qū)π盘?hào)進(jìn)行降噪,每個(gè)鄰域仿真200次，然后取輸出信噪比的平均值作為最終輸出信噪比。由圖2可知,當(dāng)鄰域選取合適時(shí),降噪后信噪比達(dá)到最大。鄰域選擇過小或過大都會(huì)使得信噪比下降。這是因?yàn)猷徲蜻x擇過小,信號(hào)受噪聲影響明顯;而鄰域選擇過大,對于相空間軌線的逐段線性逼近效果差。

圖1 Lorenz混沌信號(hào)降噪前后時(shí)域波形圖和相空間圖Fig.1 Time-domain waveform and phase space diagram of Lorenz chaotic signal before and after denoising

圖2 鄰域選擇對算法的影響Fig.2 Effect on the algorithm of neighborhood selection

仍然選取上述信號(hào),在計(jì)算得到的閾值上疊加一個(gè)隨機(jī)誤差后再進(jìn)行降噪處理。由圖3可以看出,在有誤差和無誤差時(shí)，其降噪效果基本相同,這表明算法具有較好的魯棒性。

圖3 閾值誤差對算法的影響Fig.3 Effect of threshold error on the algorithm

3.3 不同局部投影方法對比

為了進(jìn)一步研究本文所提方法的降噪效果,將其他局部投影降噪算法與本文所提方法進(jìn)行對比。選取降噪后的信噪比、均方誤差、復(fù)雜度和耗時(shí)作為衡量降噪效果的評價(jià)指標(biāo)。信噪比的計(jì)算公式為

(24)

(25)

首先,在信噪比為8 dB的情況下,計(jì)算混沌信號(hào)重要的特征量:李雅普諾夫指數(shù)、關(guān)聯(lián)維數(shù),然后再統(tǒng)計(jì)每個(gè)方法計(jì)算的耗時(shí)。每種方法設(shè)定相同參數(shù),迭代8次,仿真200次,然后計(jì)算平均耗時(shí)。由表1可以看出,由本文方法降噪后的信號(hào)計(jì)算得到的李雅普諾夫指數(shù)和關(guān)聯(lián)維數(shù)最接近原始Lorenz信號(hào)的對應(yīng)值,這表明Lorenz混沌吸引子中的確定性結(jié)構(gòu)得到了較好的保留。另外,由表1可以看出，本文方法較為耗時(shí)，這是因?yàn)楸疚姆椒ㄔ谟?jì)算過程中涉及了模糊遞歸圖的求解以及數(shù)值求解積分方程。但是本文方法在犧牲計(jì)算時(shí)間的情況下,具有更好的降噪效果。

表1 不同方法降噪后信號(hào)的李雅普諾夫指數(shù)、關(guān)聯(lián)維數(shù)以及計(jì)算耗時(shí)(SNR=8 dB)

同樣選定Lorenz混沌信號(hào),添加0～20 dB的高斯白噪聲,然后分別用這些方法進(jìn)行降噪處理。為減小偶然誤差,每種方法仿真200次,最后取200次的平均值作為最終的結(jié)果。

不同降噪方法的降噪后輸出信噪比和均方誤差的結(jié)果對比分別如圖4和圖5所示。基于小波的局部投影與基于EMD的局部投影的降噪效果接近,因?yàn)檫@兩種方法均是采用預(yù)先估計(jì)噪聲水平方式確定鄰域的。平滑子空間局部投影通過平滑正交分解識(shí)別子空間,效果優(yōu)于標(biāo)準(zhǔn)局部投影方法。遞歸局部投影采用遞歸圖計(jì)算最優(yōu)鄰域半徑,K均值局部投影利用K均值聚類方法選取鄰域,能夠?qū)崿F(xiàn)更好的去噪效果。而本文方法同時(shí)考慮了鄰域選取、子空間劃分以及質(zhì)心修正問題,相比上述方法,本文所提方法具有最高的輸出信噪比和最低的均方誤差,這說明方法較好地保留了信號(hào)時(shí)域波形的形狀,具有較好的降噪性能。

圖4 不同局部投影降噪方法的輸出信噪比Fig.4 Output signal to noise ratio of different local projective denoising methods

圖5 不同局部投影降噪方法的均方誤差Fig.5 Mean square error of different local projective denoising methods

綜上所述,本文所提的基于模糊遞歸圖和最優(yōu)硬閾值準(zhǔn)則的局部投影方法能夠有效濾除混沌信號(hào)中的噪聲。

3.4 實(shí)測ECG信號(hào)仿真

ECG可以記錄心臟的生理活動(dòng),在醫(yī)療臨床領(lǐng)域應(yīng)用廣泛。但是，ECG信號(hào)在儀器測量過程中容易被噪聲污染,從而影響各種生理特征的檢測、提取和識(shí)別,容易造成誤診,因而ECG信號(hào)降噪對實(shí)際臨床診斷具有重要的意義和價(jià)值。目前，常見的ECG信號(hào)去噪方法有奇異譜分析法、EMD方法、小波閾值方法等[31]。

ECG信號(hào)屬于非平穩(wěn)非線性信號(hào),研究表明ECG信號(hào)具有復(fù)雜的動(dòng)力學(xué)特性,表現(xiàn)出一定的混沌特性,因此非線性動(dòng)力系統(tǒng)的方法作為一種研究心臟病的有效工具,被廣泛應(yīng)用于ECG信號(hào)的研究當(dāng)中。

本文選取MIT-BIH(Massachusettes Institute of Technology-Boston’s Beth Israel Hospital)噪聲壓力測試數(shù)據(jù)庫中編號(hào)為118e12的實(shí)測數(shù)據(jù),其采樣頻率為f=360 Hz,采集到的信號(hào)含有基線漂移、噪聲等影響。

首先，利用最小二乘擬合去掉信號(hào)內(nèi)的基線漂移；然后，再使用本文方法、奇異譜分析方法、EMD方法和小波閾值方法分別對ECG信號(hào)進(jìn)行處理及對比。處理結(jié)果如圖6所示。

由圖6對比去噪前后的相圖,可以發(fā)現(xiàn)噪聲大大減少，本文方法去噪后的相圖表現(xiàn)出了較為清晰的吸引子結(jié)構(gòu),說明原來含有隨機(jī)噪聲的ECG信號(hào)在經(jīng)過降噪處理后，其確定性成分得以加強(qiáng),展現(xiàn)出了ECG信號(hào)原來的特征。對比本文方法,奇異譜分析(singular spectrum analysis, SSA)方法、EMD方法和小波閾值(wavelet threshold, WT)方法雖然也有較好的去噪效果,但時(shí)域波形還是不夠光滑,而且相空間軌線細(xì)節(jié)部分呈現(xiàn)出雜亂堆疊的現(xiàn)象。

圖6 ECG信號(hào)降噪前后波形圖和相空間圖Fig.6 Waveform and phase space diagram of ECG signal before and after denoising

為進(jìn)一步說明降噪效果,計(jì)算如圖7所示的原始信號(hào)和降噪信號(hào)在嵌入維數(shù)不同時(shí)的偽最近鄰點(diǎn)比例。當(dāng)偽最近鄰點(diǎn)比例等于0或者小于某個(gè)值而保持不變,可認(rèn)為不含偽最近鄰點(diǎn);當(dāng)時(shí)間信號(hào)受到噪聲影響時(shí),混沌吸引子軌線受到影響,偽最近鄰點(diǎn)的比例也會(huì)受到影響。由圖7可知,降噪前偽最近鄰點(diǎn)的比例在嵌入維數(shù)為4時(shí)降到較低,但沒有降至0。降噪后,當(dāng)嵌入維數(shù)為7時(shí)，偽最近鄰點(diǎn)的比例就降至0。這進(jìn)一步說明了應(yīng)用本文所提方法對實(shí)測ECG信號(hào)可進(jìn)行有效的降噪處理。

圖7 ECG降噪前后的虛假最近鄰點(diǎn)比例Fig.7 Proportion of false nearest neighbor points before and after ECG noise reduction

4 結(jié) 論

本文提出的基于模糊遞歸圖和最優(yōu)硬閾值的局部投影算法,主要針對原始局部投影算法中的鄰域選擇問題和子空間劃分問題。在相空間重構(gòu)的基礎(chǔ)上,首先使用模糊遞歸圖對鄰域進(jìn)行選擇,同時(shí)為避免求取原始方法中的協(xié)方差矩陣,本文直接使用SVD分解,然后采用最優(yōu)硬閾值準(zhǔn)則對信號(hào)子空間和噪聲子空間進(jìn)行劃分。通過對Lorenz混沌信號(hào)進(jìn)行仿真,并與其他局部投影算法相比較,表明該方法具有一定的優(yōu)越性。最后，將其應(yīng)用于實(shí)測含噪的ECG信號(hào),有效地降低了原始信號(hào)中的噪聲,使信號(hào)的特征更加明顯。同時(shí)，將本文所提算法與其他類型的ECG信號(hào)降噪方法進(jìn)行了對比,證明了本文方法的有效性,為后續(xù)的信號(hào)處理奠定了基礎(chǔ)。