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

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

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

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

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

1 方 法

1.1 局部投影降噪算法

假設x(t)(t=1,2,…,l)是動力系統輸出的一組觀測時間序列.Takens嵌入定理已經證明,當嵌入維數m>2D2(D2為系統的關聯維數)時,從x(t)中每隔時間τ取出一個觀察值作為時間序列的狀態(tài)測量值,即可重構出與原動力系統具有相同拓撲特性的m維相空間.重構得到的信號為

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

(1)

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

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

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

1.2 基于經驗模態(tài)分解的鄰域選取策略

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

經驗模態(tài)分解方法認為任何時間序列都是由多個具有不同特征尺度的本征模態(tài)函數(intrinsic mode function,IMF)組成的,即

(2)

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

自相關函數是反映信號與其自身在不同時間尺度上相似程度的統計量,隨機噪聲在不同時間尺度上的隨機性使其與自身表現為弱相關性,而一般信號則表現為較強相關性,因此可用自相關函數區(qū)分隨機噪聲和一般信號.

(3)

其中,M是重構噪聲分量的絕對中值偏差，

(4)

定義鄰域半徑為

(5)

其中,σdata是實測時間序列的標準差,即根據估得的噪聲標準差與原時間序列的標準差之比確定鄰域半徑.

1.3 級聯自適應局部投影降噪算法

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

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

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

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

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

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

(6)

2)鄰域劃分.對相空間中的任意相點X(i),定義其鄰域為

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

(7)

其中:X(j)為滿足條件的鄰點;‖·‖表示歐幾里得距離;r為鄰域半徑,由前述基于經驗模態(tài)分解方法得到.

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

(8)

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

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

(9)

其中,為避免各個相點首尾坐標畸變帶來的誤差,修正后的相點可表示為

(10)

4) 更新所有相點,對修正后的所有相點進行反重構.

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

圖3 級聯降噪方法框圖

2 結果與討論

2.1 Lorenz系統時間序列實驗

本節(jié)以Lorenz系統產生的混沌時間序列檢驗所提算法的有效性,Lorenz系統的動力學方程為

(11)

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

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

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

2.2 脈搏信號實驗

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

圖4 Lorenz信號降噪前后相空間重構圖

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

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

圖6 實測脈搏信號降噪前后時域和相空間重構圖像

2.3 心電信號實驗

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

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

圖7 實測心電信號降噪前后時域和相空間重構圖像

3 結 語

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