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

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

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

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

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

1 方 法

1.1 局部投影降噪算法

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

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

(1)

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

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

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

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

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

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

(2)

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

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

(3)

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

(4)

定義鄰域半徑為

(5)

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

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

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

圖1 含噪脈搏信號和兩次局部投影降噪后信號

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

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

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

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

(6)

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

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

(7)

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

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

(8)

圖2 原始信號與三次局部投影降噪后的信號對比

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

(9)

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

(10)

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

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

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

2 結(jié)果與討論

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

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

(11)

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié) 語

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