級(jí)聯(lián)自適應(yīng)局部投影降噪方法

徐禮勝, 崔慧穎, 吳俊鼎, 王仲怡

(1. 東北大學(xué) 醫(yī)學(xué)與生物信息工程學(xué)院, 遼寧 沈陽(yáng) 110169; 2. 沈陽(yáng)東軟智能醫(yī)療科技研究院有限公司, 遼寧 沈陽(yáng) 110167)

信號(hào)在采集過(guò)程中易受不同程度的噪聲干擾,這些噪聲破壞信號(hào)結(jié)構(gòu),降低了識(shí)別信號(hào)特征的準(zhǔn)確性[1];因此,噪聲抑制是信號(hào)研究的基礎(chǔ)[2].當(dāng)信號(hào)頻譜與噪聲頻譜存在明顯差別時(shí),線(xiàn)性濾波器效果顯著[3-4].而產(chǎn)生自混沌系統(tǒng)的信號(hào)通常具有寬頻譜的特性,信號(hào)頻譜與噪聲頻譜重疊時(shí),傳統(tǒng)線(xiàn)性濾波方法會(huì)導(dǎo)致信號(hào)失真.

基于非線(xiàn)性動(dòng)力學(xué)理論的局部投影算法(local projection, LP)為寬頻譜信號(hào)的降噪提供了一種新的選擇.期望信號(hào)與噪聲在相空間的不同維度子空間中具有不同的結(jié)構(gòu),局部投影算法則據(jù)此區(qū)分兩者.依據(jù)Takens嵌入定理[5]重構(gòu)信號(hào)得到的嵌入空間即為相空間,相空間中的光滑非線(xiàn)性流形反映系統(tǒng)狀態(tài),因此系統(tǒng)產(chǎn)生的期望信號(hào)在相空間中的軌跡位于或非?？拷饣蔷€(xiàn)性流形[6],而噪聲的相空間軌跡因噪聲并不具備任何來(lái)源約束,因此分散在流形四周,甚至可能位于流形上.LP算法在重構(gòu)相空間中為每個(gè)相點(diǎn)劃分鄰域,對(duì)鄰域內(nèi)的非線(xiàn)性流形進(jìn)行局部線(xiàn)性逼近.依據(jù)特征值分解理論區(qū)分信號(hào)空間和噪聲空間,通過(guò)向噪聲子空間投影濾去噪聲成分.該算法目前已被成功應(yīng)用于呼吸音信號(hào)濾波[7]、核磁共振激光數(shù)據(jù)降噪[8]、腦機(jī)接口腦電圖增強(qiáng)[9]、語(yǔ)音信號(hào)分析[10-11]等.

自適應(yīng)鄰域選取是局部投影算法的研究重點(diǎn)之一.Przybya等[12]使用K-means聚類(lèi)算法形成鄰域,取得了較好的降噪效果.Chelidze[13]重構(gòu)了信號(hào)的短時(shí)軌跡鏈?zhǔn)?并利用平滑正交分解識(shí)別局部子空間,降噪效果優(yōu)于基于傳統(tǒng)正交分解的降噪方法.經(jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decomposition,EMD)作為一種完全由數(shù)據(jù)驅(qū)動(dòng)的方法,不需要任何先驗(yàn)知識(shí)即可自適應(yīng)地將復(fù)雜數(shù)據(jù)分解為有限個(gè)本征模態(tài)函數(shù),適用于非線(xiàn)性時(shí)間序列的處理[14-16].因此,本文提出基于經(jīng)驗(yàn)?zāi)B(tài)分解的鄰域選取策略,根據(jù)自相關(guān)函數(shù)確定分解得到的噪聲分量,再通過(guò)評(píng)估噪聲分量在原始信號(hào)中的占比確定鄰域半徑.此外,目前有許多結(jié)合小波理論對(duì)信號(hào)進(jìn)行局部投影降噪處理的研究[17-19],Wan等[20]通過(guò)將局部投影算法應(yīng)用于腦電圖的小波變換,去除了噪聲干擾.小波閾值去噪方法理論成熟,降噪效果明顯,但其在恢復(fù)信號(hào)非線(xiàn)性流形結(jié)構(gòu)方面性能不佳.LP算法通過(guò)將噪聲分量投影到低維子空間中以減少其與非線(xiàn)性流形的偏差,能夠最大程度地恢復(fù)信號(hào)的真實(shí)流形結(jié)構(gòu).綜上,本文提出局部投影算法與小波閾值去噪級(jí)聯(lián)的降噪方法.

1 方 法

1.1 局部投影降噪算法

假設(shè)x(t)(t=1,2,…,l)是動(dòng)力系統(tǒng)輸出的一組觀測(cè)時(shí)間序列.Takens嵌入定理已經(jīng)證明,當(dāng)嵌入維數(shù)m>2D2(D2為系統(tǒng)的關(guān)聯(lián)維數(shù))時(shí),從x(t)中每隔時(shí)間τ取出一個(gè)觀察值作為時(shí)間序列的狀態(tài)測(cè)量值,即可重構(gòu)出與原動(dòng)力系統(tǒng)具有相同拓?fù)涮匦缘膍維相空間.重構(gòu)得到的信號(hào)為

X(t)=[x(t),x(t+τ),…,x(t+(m-1)τ)]T.

(1)

其中,t=1,2,…,l-(m-1)τ.

理論上,當(dāng)嵌入維數(shù)大于或等于m時(shí),系統(tǒng)的動(dòng)力學(xué)特性即可得到充分展現(xiàn)[21].但是,更大的嵌入維數(shù)也意味著更大的計(jì)算成本.因此,選取一個(gè)大于2D2的最小嵌入維數(shù)是必要的.基于嵌入窗思想的C-C算法[22]以計(jì)算量小、結(jié)果準(zhǔn)確的特點(diǎn)成為目前最常用的計(jì)算嵌入維數(shù)和時(shí)間延遲的方法.

在重構(gòu)的相空間中,原始動(dòng)力系統(tǒng)產(chǎn)生的混沌吸引子僅存在于維數(shù)為m0的低維子空間,m0維子空間的正交補(bǔ)空間叫做零子空間,其維數(shù)為Q=m-m0.信號(hào)無(wú)噪聲的情況下,零子空間中無(wú)元素;當(dāng)信號(hào)含噪聲時(shí),噪聲因其自身的無(wú)方向性存在于任何子空間.顯然,此時(shí)零子空間中的元素全部由噪聲產(chǎn)生.因此,把鄰域內(nèi)的所有相點(diǎn)投影到m0維零子空間中,即可將低維動(dòng)力學(xué)成分從高維空間中分離出來(lái),使混沌吸引子流形回歸到真實(shí)形態(tài).

1.2 基于經(jīng)驗(yàn)?zāi)B(tài)分解的鄰域選取策略

鄰域半徑是保證局部投影算法降噪質(zhì)量的重要參數(shù)[23].半徑太小會(huì)導(dǎo)致部分相點(diǎn)丟失,使鄰域內(nèi)的期望信號(hào)被噪聲淹沒(méi);半徑太大則會(huì)包含過(guò)多無(wú)用相點(diǎn),導(dǎo)致局部線(xiàn)性逼近的結(jié)果不準(zhǔn)確.為避免憑經(jīng)驗(yàn)選取鄰域半徑的盲目性,本文利用經(jīng)驗(yàn)?zāi)B(tài)分解方法估計(jì)觀測(cè)時(shí)間序列的噪聲水平,據(jù)此選取半徑進(jìn)行鄰域的劃分.

經(jīng)驗(yàn)?zāi)B(tài)分解方法認(rèn)為任何時(shí)間序列都是由多個(gè)具有不同特征尺度的本征模態(tài)函數(shù)(intrinsic mode function,IMF)組成的,即

(2)

式中:x(t)為原始信號(hào);imfi(t)(i=1,2,…,n)為分解得到的一組IMF;rn(t)為殘余分量.IMF分量的頻率從高到低代表了x(t)的不同成分,而噪聲通常表現(xiàn)出高頻特性.確定代表噪聲的IMF分量后,即可評(píng)估信號(hào)的噪聲水平,從而選取合適的鄰域范圍.

自相關(guān)函數(shù)是反映信號(hào)與其自身在不同時(shí)間尺度上相似程度的統(tǒng)計(jì)量,隨機(jī)噪聲在不同時(shí)間尺度上的隨機(jī)性使其與自身表現(xiàn)為弱相關(guān)性,而一般信號(hào)則表現(xiàn)為較強(qiáng)相關(guān)性,因此可用自相關(guān)函數(shù)區(qū)分隨機(jī)噪聲和一般信號(hào).

(3)

其中,M是重構(gòu)噪聲分量的絕對(duì)中值偏差，

(4)

定義鄰域半徑為

(5)

其中,σdata是實(shí)測(cè)時(shí)間序列的標(biāo)準(zhǔn)差,即根據(jù)估得的噪聲標(biāo)準(zhǔn)差與原時(shí)間序列的標(biāo)準(zhǔn)差之比確定鄰域半徑.

1.3 級(jí)聯(lián)自適應(yīng)局部投影降噪算法

局部投影降噪算法對(duì)原始信號(hào)的非線(xiàn)性流形有明顯的修復(fù)作用,但對(duì)高頻噪聲的去除效果不如頻域降噪方法[24].如圖1所示,脈搏信號(hào)在經(jīng)過(guò)兩次局部投影降噪處理后仍殘留部分噪聲.

圖1 含噪脈搏信號(hào)和兩次局部投影降噪后信號(hào)

傳統(tǒng)局部投影方法通過(guò)迭代的方式去除殘留噪聲,實(shí)際上,2～3次迭代之后信號(hào)的信噪比提升量會(huì)逐漸收斂.如圖2所示,持續(xù)迭代使信號(hào)平滑的代價(jià)是幅值的持續(xù)衰減,這種丟失信息的降噪方法顯然是不可取的.

本文方法只進(jìn)行1～2次局部投影處理,這樣可以盡可能多地保留信號(hào)非線(xiàn)性特征,并且信噪比得到了大幅提升,而小波閾值方法正適用于中等信噪比的信號(hào)降噪處理.雖然小波閾值方法能夠較好地抑制高頻噪聲,但恢復(fù)混沌吸引子流形結(jié)構(gòu)的能力有限[25].因此,本文綜合兩種方法的優(yōu)點(diǎn),1～2次局部投影后再通過(guò)小波閾值方法去除殘留的少量噪聲成分即可達(dá)到既保留信號(hào)特征又平滑噪聲的目的,本文的降噪方法流程如圖3所示,具體步驟如下:

1) 相空間重構(gòu).觀測(cè)時(shí)間序列在m維重構(gòu)相空間中可表示為

X=[X(1),X(2),…,X(N)],N=l-(m-1)τ.

(6)

2)鄰域劃分.對(duì)相空間中的任意相點(diǎn)X(i),定義其鄰域?yàn)?/p>

Ui={X(j)‖X(j)-X(i)‖

(7)

其中:X(j)為滿(mǎn)足條件的鄰點(diǎn);‖·‖表示歐幾里得距離;r為鄰域半徑,由前述基于經(jīng)驗(yàn)?zāi)B(tài)分解方法得到.

3)局部投影.定義協(xié)方差矩陣：

(8)

圖2 原始信號(hào)與三次局部投影降噪后的信號(hào)對(duì)比

計(jì)算C(i)的特征向量和對(duì)應(yīng)的特征值.理論上,特征向量對(duì)應(yīng)的方向可以衡量信號(hào)與噪聲的不同方向.較大特征值對(duì)應(yīng)的特征向量所在的方向代表信號(hào)子空間;反之,較小的Q個(gè)特征值λ1,λ2,…,λQ對(duì)應(yīng)的特征向量eq(q=1,…,Q)所在的方向代表了由噪聲組成的零子空間[8].引入權(quán)矩陣R對(duì)坐標(biāo)進(jìn)行修正以避免邊界效應(yīng),定義r11=rmm=1 000,其余元素全為1[26].各個(gè)相點(diǎn)在m0維子空間中的投影即為噪聲分量,

(9)

其中,為避免各個(gè)相點(diǎn)首尾坐標(biāo)畸變帶來(lái)的誤差,修正后的相點(diǎn)可表示為

(10)

4) 更新所有相點(diǎn),對(duì)修正后的所有相點(diǎn)進(jìn)行反重構(gòu).

5) 對(duì)反重構(gòu)得到的信號(hào)進(jìn)行小波閾值降噪處理,抑制高頻噪聲.

圖3 級(jí)聯(lián)降噪方法框圖

2 結(jié)果與討論

2.1 Lorenz系統(tǒng)時(shí)間序列實(shí)驗(yàn)

本節(jié)以L(fǎng)orenz系統(tǒng)產(chǎn)生的混沌時(shí)間序列檢驗(yàn)所提算法的有效性,Lorenz系統(tǒng)的動(dòng)力學(xué)方程為

(11)

選取σ=10,R=28,b=8/3,設(shè)定積分步長(zhǎng)h=0.003,此時(shí)系統(tǒng)表現(xiàn)為混沌狀態(tài)，采用四階龍格-庫(kù)塔法求解方程.取x分量上的4 000個(gè)數(shù)據(jù)點(diǎn)作為實(shí)驗(yàn)信號(hào),以添加噪聲后信噪比為20 dB的信號(hào)為例,首先對(duì)仿真信號(hào)進(jìn)行經(jīng)驗(yàn)?zāi)B(tài)分解,并計(jì)算鄰域半徑用于局部投影降噪.圖4分別為原始信號(hào)、含噪信號(hào)、小波閾值方法降噪信號(hào)和本文方法降噪信號(hào)的相空間重構(gòu)圖.圖4b顯示,信噪比為20 dB時(shí)吸引子流形已經(jīng)完全被噪聲破壞,無(wú)法分辨.圖4c展示了小波閾值去噪方法的處理結(jié)果,高頻噪聲雖然得到了抑制但吸引子流形的結(jié)構(gòu)較為混亂.當(dāng)嵌入維數(shù)m=11,時(shí)間延遲τ=1,鄰域半徑r=0.592 8,并選用“db4”小波基時(shí),本文方法的降噪結(jié)果如圖4d所示,信號(hào)的相空間軌跡得到了基本恢復(fù),并且流形結(jié)構(gòu)更加有序,顯著改善了單獨(dú)使用小波閾值方法的不足.特別地,對(duì)于吸引子流形中軌跡密集的位置,本文方法的處理結(jié)果具有更清晰的細(xì)節(jié)特征,更加接近原始信號(hào),這也說(shuō)明利用本文方法選取的鄰域可以提高逼近非線(xiàn)性流形的準(zhǔn)確率.

為定量評(píng)估算法性能,向Lorenz信號(hào)中添加不同信噪比的噪聲,分別計(jì)算小波閾值降噪與本文方法降噪后信號(hào)的信噪比(signal to noise ratio,SNR)和均方誤差(mean square error,MSE).通常來(lái)說(shuō),SNR越大MSE越小的降噪后的信號(hào)越接近原始信號(hào).結(jié)果如表1所示,相較于小波閾值去噪方法,本文方法在中高強(qiáng)度的噪聲范圍內(nèi)能獲得更高的SNR和更低的MSE,說(shuō)明該方法在濾除噪聲的同時(shí)能夠最大程度地接近原始信號(hào).

表1 不同噪聲強(qiáng)度下Lorenz信號(hào)降噪后的信噪比和均方誤差

2.2 脈搏信號(hào)實(shí)驗(yàn)

分別采集橈動(dòng)脈、肱動(dòng)脈、頸動(dòng)脈三處的脈搏波.為對(duì)比評(píng)估本文算法在脈搏信號(hào)降噪中的性能,將三處實(shí)測(cè)信號(hào)降噪后得到的信號(hào)作為干凈信號(hào),采用SNR和MSE作為評(píng)價(jià)指標(biāo),分別使用小波閾值去噪方法、本文方法對(duì)添加不同強(qiáng)度噪聲的三處脈搏信號(hào)進(jìn)行降噪處理.圖5結(jié)果顯示本文方法整體優(yōu)于小波閾值去噪方法,三處脈搏信號(hào)的實(shí)驗(yàn)結(jié)果均顯示本文方法能夠獲得更大的信噪比和更小的均方誤差.噪聲能量較大時(shí),信號(hào)幾乎完全被噪聲覆蓋,兩種方法性能差別較小,都不能很好地恢復(fù)信號(hào).噪聲能量相對(duì)較小時(shí),信號(hào)結(jié)構(gòu)受噪聲影響較小,局部投影方法恢復(fù)信號(hào)結(jié)構(gòu)的性能發(fā)揮程度有限,而小波閾值去噪方法對(duì)此類(lèi)噪聲敏感性更強(qiáng).因此,隨著SNR的提高,兩種方法在降噪效果上的差別逐漸變小.中等噪聲強(qiáng)度下,本文方法的優(yōu)勢(shì)更加明顯.此時(shí),含噪信號(hào)經(jīng)過(guò)局部投影處理,信號(hào)基本結(jié)構(gòu)得以恢復(fù),為小波閾值去噪方法平滑噪聲的處理提供了前期基礎(chǔ).相比單獨(dú)使用小波閾值去噪方法,本文方法明顯提高了信噪比,降低了均方誤差.其中,本文方法在頸動(dòng)脈脈搏信號(hào)的降噪處理中優(yōu)勢(shì)更加明顯.頸動(dòng)脈脈搏波波形面積大,主波下降速度慢,潮波、重搏波位置更高,信號(hào)波幅變化相對(duì)較小,主波持續(xù)時(shí)間更長(zhǎng),與橈動(dòng)脈、肱動(dòng)脈脈搏波形態(tài)差別較大.“db4”小波基無(wú)法很好地匹配該波形,致使小波閾值方法性能有所下降.本文方法在不同部位的脈搏信號(hào)表現(xiàn)出穩(wěn)定的降噪效果.

圖4 Lorenz信號(hào)降噪前后相空間重構(gòu)圖

圖5 兩種方法在不同位置脈搏波上的降噪結(jié)果對(duì)比

圖6展示了對(duì)實(shí)測(cè)脈搏信號(hào)分別使用小波閾值去噪方法、本文方法的相空間重構(gòu)結(jié)果和時(shí)域圖像,其中時(shí)間延遲τ=1、嵌入維數(shù)m=6、鄰域半徑r=1.369 7,選用“db4”小波基.原始信號(hào)如圖6a所示,從時(shí)域圖像中尚可觀察到信號(hào)的部分周期特征,但相空間軌跡已無(wú)法分辨.對(duì)比圖6b與圖6c,小波閾值去噪方法雖然使信號(hào)的重構(gòu)軌跡變得平滑,但結(jié)構(gòu)雜亂,軌跡堆疊現(xiàn)象明顯,時(shí)域圖像中各周期波形形態(tài)差異較大;而本文方法在恢復(fù)信號(hào)相空間基本結(jié)構(gòu)的基礎(chǔ)上去除了殘留的高頻噪聲,得到的信號(hào)曲線(xiàn)光滑,軌跡更加收斂,混沌吸引子呈現(xiàn)出有規(guī)律的幾何結(jié)構(gòu),且信號(hào)時(shí)域波形規(guī)整有序.實(shí)測(cè)脈搏信號(hào)降噪結(jié)果說(shuō)明本方法恢復(fù)信號(hào)原始特征的能力明顯優(yōu)于小波方法.

圖6 實(shí)測(cè)脈搏信號(hào)降噪前后時(shí)域和相空間重構(gòu)圖像

2.3 心電信號(hào)實(shí)驗(yàn)

為驗(yàn)證改進(jìn)方法在生理信號(hào)降噪方面的普適性,本節(jié)選擇實(shí)測(cè)心電信號(hào)(electrocardiogram signal, ECG)作為驗(yàn)證對(duì)象.心電信號(hào)在采集過(guò)程中易受噪聲干擾,噪聲會(huì)對(duì)心電信號(hào)QRS波的檢測(cè)產(chǎn)生不良影響,降低計(jì)算機(jī)自動(dòng)診斷的準(zhǔn)確性[27].心電信號(hào)中的噪聲主要有基線(xiàn)漂移、工頻干擾和肌電干擾[28].其中,肌電干擾的頻譜范圍較寬,與ECG頻譜存在重疊,傳統(tǒng)方法很難將此類(lèi)噪聲與心電信號(hào)完全分離.

圖7a中原始心電信號(hào)受到噪聲干擾.P波、Q波、S波和T波被覆蓋,其對(duì)應(yīng)的相空間軌跡雜亂重疊,無(wú)法分辨.小波方法對(duì)心電信號(hào)中的噪聲有一定的抑制作用,圖7b中的信號(hào)波形得到了大致恢復(fù),但仍存在明顯瑕疵,不同周期中的P波、Q波、S波和T波恢復(fù)程度差異較大,更換不同小波基所得結(jié)果并無(wú)明顯差別,說(shuō)明小波方法的結(jié)果對(duì)小波基的選取依賴(lài)性較大.利用本文方法對(duì)心電信號(hào)作降噪處理,取時(shí)間延遲τ=1、嵌入維數(shù)m=9、鄰域半徑r=0.257 1,選用“db4”小波基.由圖7c可知,本文方法對(duì)心電信號(hào)中的噪聲有明顯抑制作用且降噪效果穩(wěn)定可靠,并不受信號(hào)形態(tài)影響,降噪后的心電信號(hào)完全平滑,形態(tài)穩(wěn)定,各個(gè)周期波峰清晰.相空間軌跡結(jié)構(gòu)穩(wěn)定,且幅值較小的波所對(duì)應(yīng)的軌跡同樣得到了基本恢復(fù),規(guī)整有序.

圖7 實(shí)測(cè)心電信號(hào)降噪前后時(shí)域和相空間重構(gòu)圖像

3 結(jié) 語(yǔ)

具有非線(xiàn)性、非平穩(wěn)性特征的微弱信號(hào)易受噪聲干擾,而傳統(tǒng)線(xiàn)性濾波器在期望信號(hào)頻帶與噪聲頻帶發(fā)生重疊時(shí)性能不佳.針對(duì)此類(lèi)問(wèn)題,本文提出了一種基于EMD的噪聲自適應(yīng)鄰域選取策略和級(jí)聯(lián)局部投影降噪方法.通過(guò)向Lorenz信號(hào)中添加不同程度高斯白噪聲驗(yàn)證提出方法的有效性,不同降噪方法下含噪Lorenz信號(hào)去噪前后的SNR和MSE結(jié)果表明,本方法能夠有效去除強(qiáng)高斯白噪聲,并在恢復(fù)信號(hào)非線(xiàn)性動(dòng)力學(xué)特性方面更有效.心血管解剖結(jié)構(gòu)的分形特點(diǎn)及生物電傳導(dǎo)系統(tǒng)的分形結(jié)構(gòu)使與之密切相關(guān)的生理信號(hào)表現(xiàn)出混沌現(xiàn)象與非線(xiàn)性動(dòng)力學(xué)行為.因此,本文量化分析了所提方法在脈搏信號(hào)和心電信號(hào)中的降噪效果.脈搏信號(hào)仿真實(shí)驗(yàn)的定量結(jié)果表明,與小波閾值降噪方法相比,本文方法抑制噪聲的能力有明顯提高.實(shí)測(cè)心電信號(hào)降噪結(jié)果也證實(shí)本文方法消除噪聲干擾的能力顯著優(yōu)于小波方法.降噪前后的信號(hào)相空間重構(gòu)結(jié)果表明本文方法能夠在中、高強(qiáng)度噪聲下恢復(fù)信號(hào)的非線(xiàn)性結(jié)構(gòu),保留信號(hào)特征.