梁慶真，郭興明，袁志會

(重慶大學 生物工程學院，生物流變科學與技術教育部重點實驗室，重慶 400044)

心音信號是一種非線性非平穩(wěn)的微振動信號，在采集的過程中不可避免的會引入干擾噪聲，使得部分有用信息淹沒在噪聲中[1]，因此在對其進行進一步分析之前，心音去噪顯得尤為重要。目前，小波去噪在心音去噪中已經(jīng)得到廣泛的應用，它克服了傅里葉分析不能對非平穩(wěn)信號進行有效分析的缺點，但是心音信號的頻帶與噪聲頻帶存在部分或全部重疊，所以采用傳統(tǒng)的小波分解重構法對心音進行線性去噪時，會相應的去除掉與噪聲頻帶相重疊的部分心音，在這個過程中信號本身不可避免的會損失掉部分非線性特征信息。采用非線性小波變換函數(shù)[2]的離散小波閾值去噪法(Discrete Wavelet Threshold,DWT)能較好的抑制噪聲，且能得到更精確的重構信號[3]。基于相空間重構技術的局部投影降噪法(Local Projection,LP)[4-5]，它將一維的時間序列拓展到高維相空間，利用系統(tǒng)信號和噪聲信號在高維相空間上分布特性的差異來濾除噪聲信號，在濾除噪聲的同時可以有效保留微弱特征信號成分[6]。本文綜合兩者的優(yōu)點，提出噪聲水平自適應估計的局部投影與離散小波閾值相結合的去噪方法，該方法既能得到精確的重構信號又能保留信號的非線性特征，采用該方法對Lorenz序列以及實測心音進行降噪研究，用信噪比和均方誤差來檢驗該方法的去噪效果，并提取去噪前后信號的最大Lyapunov指數(shù)來說明該方法對信號非線性特征的影響。同時，為了體現(xiàn)本文方法的優(yōu)勢，還將與局部投影和離散小波閾值去噪進行對比分析。

1 基本原理

1.1 噪聲水平自適應估計的局部投影降噪

Zn=[Zn,Zn+τ,Zn+2τ,…,Zn+(m-1)τ]

(1)

在理想的無噪聲(Wn=0)的情況下，有如下的相空間映射關系：

(2)

當系統(tǒng)存在噪聲(Wn≠0)時，式(2)可表示為：

(3)

式(2)在Zn的鄰域Nn內可線性化展開為：

(4)

(5)

Nn,

k=[1,n)∪(n,L-(m-1)τ]}

(6)

同理，式(3)在Zn的鄰域內線性化展開為：

(7)

其中:ηn為噪聲總和。式(7)表明在m維的相空間里，存在一個超平面(由鄰域內相點局部線性化得到)，相點與超平面的偏差代表了噪聲部分。在這個m維的相空間里，存在著一個吸引子空間和零子空間(由Q個相互正交的向量組成，可記為P=[P1,P2,…,PQ])，由系統(tǒng)動力學行為所決定的吸引子只被局限在吸引子空間中，當系統(tǒng)受到噪聲污染時，整個相空間中都會有分量，必定是噪聲使得零子空間產生了分量，局部投影降噪算法就是通過向吸引子空間投影的方法來濾除噪聲成分。投影后的修正值可表示為：

(8)

1.2 離散小波閾值降噪

帶噪聲的信號模型可表示為：

(9)

其中:f(i)為真實信號，y(i)為含噪信號，e(i)是噪聲。一般情況下，噪聲是方差為σ2的高斯噪聲，服從N(0,σ2)。通常噪聲不可知，由SURE程序來估計閾值，定義如下：

(10)

離散閾值小波的去噪方法[10]為：

(1) 對含噪信號y(i)做小波變換得小波系數(shù)Wj,k。

(2) 對Wj,k進行閾值處理,得到估計小波系數(shù)，

2 局部投影和離散小波閾值降噪算法

局部投影降噪是將時間序列重構的相空間劃分成不同的子空間，特征值較大那部分是背景信號(相對較強的信號)，特征值較小的部分是微弱信號和噪聲，故可以多次采用局部投影降噪結合小波去噪算法來提取完整信號，先取較大的m，用局部投影降噪去除隨機噪聲部分，保留了背景信號和微弱信號，然后取較小的m，再次用局部投影降噪算法，這時特征值較大的信號為背景信號，特征值較小的即為微弱信號，但是微弱信號中還包含一些噪聲，而離散小波閾值去噪可以有效的保留信號信息，所以再用離散小波閾值降噪將此噪聲去除，即可得到較完整的心音信號。算法步驟如下：

(11)

再計算Cn，Cn=(R·An)T(R·An),R為m階對角權重矩陣，用來抑制相點兩端元素畸變，保留穩(wěn)定的中間元素，其中R11=Rmm=1 000。

(3) 計算矩陣Cn的特征值Λ和特征向量P，其中Λ=diag(λ1,λ2,…,λm)，且(λ1≥λ2≥…≥λm)，P=[P1,P2,…,Pm]，對于特征向量P，[P1,P2,…,Pm0]且(m0

(7) 修正去噪后的相點，并更新所有相點，對更新后的相點進行加權反重構處理。

算法分離噪聲的流程圖如圖1。

圖1 分離噪聲的流程圖

3 實驗及結果分析

3.1 數(shù)值仿真分析

為驗證本文算法的有效性，以Lorenz系統(tǒng)產生的混沌吸引子為例進行分析。

Lorenz系統(tǒng)的混沌動力學方程為：

x′=σ(y-x)y′=x(R-z)-yz′=xy-bz(12)

其中:σ=10,R=28,b=8/3。設定積分步長h=0.003，取x方向的12 000個點作為研究對象，添加SNR為4.142 6高斯白噪聲，重構嵌入維數(shù)m=11，時延τ=20。這時特征值較大的表征吸引子整體流行的信號，特征值較小的表征精細結構的信號，二者結合白噪聲構成了該Lorenz信號。分別用離散小波閾值(DWT)、局部投影(LP)以及本文算法對其進行去噪。DWT法選取了去噪效果較好的db5小波，并采用了SURE閾值估計。

圖2 Lorenz系統(tǒng)降噪效果圖

圖2中的(a)、(b)、(c)、(d)、(e)分別為原始信號、加噪信號、DWT去噪信號、LP去噪信號和本文方法去噪信號?？梢钥闯?，直接用DWT去噪效果不是很理想，LP法已基本保留了特征值較大和較小的信號成分，已將大的隨機噪聲去除，但是去噪后的曲線在細微處依然不是很平滑，說明特征值較小的微弱信號中還含有部分隨機噪聲，影響了Lorenz的精細結構；本文方法又用DWT法將微弱信號中隨機噪聲去除，去噪后的吸引子曲線較為平滑，有著較為精細的結構，與原始信號吸引子大致相同，對噪聲有明顯的抑制作用。

為了檢驗該算法的去噪效果，選擇了兩個評價指標——信噪比(SNR)和均方誤差(MSE)。去噪后SNR越高，且MSE越小，表明越接近原始信號，去噪效果就越好。在不同信噪比的噪聲下，分別計算三種方法降噪后的信噪比和均方誤差，如表1所示。

表1 不同方法降噪后的信噪比和均方誤差

由表1可知，在SNR不同的噪聲中，本文方法較DWT法和LP算法都有較高的SNR和較低的MSE，說明該方法能有效的濾除噪聲，且能較大程度的逼近原始信號。

最大Lyapunov指數(shù)是診斷和描述非線性混沌系統(tǒng)的重要參數(shù)[11]，最能表征信號的非線性混沌特征。本文采用了改進的Rosenstein方法[12-13]計算最大Lyapunov指數(shù)，它對重構時延，維數(shù)和噪聲水平具有一定的魯棒性，算出的最大Lyapunov指數(shù)較為準確。提取加噪前后以及去噪后該系統(tǒng)的最大Lyapunov指數(shù)如表2所示。

表2 不同方法降噪后的最大Lyapunov指數(shù)

由表2可知，加噪后的最大Lyapunov指數(shù)明顯偏離實際值，導致有效信息丟失，說明噪聲對最大Lyapunov指數(shù)這一特性參數(shù)有影響。DWT和LP算法雖然能使這一參數(shù)朝著原始信號逼近，但是存在較大的誤差，不能作為原始信號特征值的近似，本文方法使該參數(shù)在較小的誤差范圍內等同于未加噪時的原始信號，證實該去噪方法能有效保留信號的非線性特征。

3.2 實測信號分析

心音信號由運動心力檢測儀(ECCM，專利號01256971.2)采集，采樣頻率為11 025 Hz。正常心音(Normal Heart Sound，NHS)、二尖瓣狹窄(Mitral Stenosis，MS)、主動脈關閉不全(Aortic Incompetence，AI)三類心音的采樣的時間序列長度分別為7 300、6 000、4 800。LP法中由互信息法和CAO算法求得三類心音的延遲時間τ和嵌入維數(shù)m見表3。

表3 不同類型心音信號的延遲時間τ和嵌入維數(shù)m

限于篇幅，本文僅列出正常心音信號的相空間降噪圖。對一例含噪的正常心音進行相空間重構，分別用DWT、LP和本文方法進行降噪的效果如圖3。

圖4為不同降噪算法的去噪時序圖，圖中的?、?、?、?分別代表原始信號、DWT去噪、LP去噪和本文方法去噪后的信號。

圖3 各種降噪算法的相空間去噪效果圖

圖4 三種類型心音的降噪時序圖

圖3～4中可以看出，DWT法去噪后相空間曲線流形依然雜亂，LP法降噪效果優(yōu)于DWT法，已經(jīng)趨于平滑；從時序圖來看，在病理信號上兩者都保留了體現(xiàn)信號病理特征的高頻部分，保證了心音信號的完整性，但是LP法去噪后的信號還有些許毛刺；本文方法降噪后的心音信號相空間圖曲線較為平滑，也比較清晰，去噪后的時序圖幾乎沒有毛刺，也保留了心音信號的精細結構，說明該方法既能抑制噪聲又能得到較為精確的信號，從而表明該方法的有效性。

提取去噪后信號的最大Lyapunov指數(shù)如表4所示(表中不含噪信號的最大Lyapunov指數(shù)是相關學者研究的數(shù)據(jù)[14])?？梢钥闯觯N類型心音的最大Lyapunov指數(shù)均為正值，表明心音信號中存在混沌現(xiàn)象[15]，其中正常心音的最大Lyapunov指數(shù)最小，說明病理心音較正常心音混沌性強，正常心音具有較強的周期性。并且DWT法和LP法得出的最大Lyapunov指數(shù)已經(jīng)明顯偏離了不含噪信號的最大Lyapunov指數(shù)值，不能作為不含噪信號非線性特征的近似，而本文算法得出的最大Lyapunov指數(shù)在較小的誤差范圍內等同于不含噪信號，表明該方法有效保留了信號的非線性特征。

表4 不同方法降噪后的最大Lyapunov指數(shù)

4 結 論

綜合局部投影和離散小波去噪的優(yōu)點，本文提出了噪聲自適應估計的局部投影與離散小波閾值相結合的降噪方法，并提取信號的最大Lyapunov指數(shù)這一表征信號非線性的特征量，通過Lorenz序列以及實測心音這種微振動生理信號的實驗分析，證實了該方法能有效濾除噪聲且能在誤差范圍內逼近不含噪信號，得到精確的重構信號，保留了更多的非線性特征信息，對比局部投影和離散小波閾值去噪算法，表明了該方法的優(yōu)越性，為心音信號的后續(xù)處理奠定了基礎。

[1] 韓陽.基于心音信號分析的心臟病診斷方法研究[D].蘭州：蘭州理工大學,2012.

[2] 熊新兵,陳亞光.非線性小波變換及其應用[J].中南民族大學學報(自然科學版),2005,24(2):43-45.

XIONG Xin-bing, CHEN Ya-guang. Nonlinear wavelet transform and its application[J]. Journal of south-central university for nationalities(Natural Science Edition),2005,24(2): 43-45.

[3] Emeterio J L S, RodriguezHernandez M A. Wavelet denoising of ultrasonic a-scans for detection of weak signals[C]. //IEEE on 19th International Conference on Systems, Signals and Image Processing, 2012:48-51.

[4] 張來斌,李峰,段禮祥.基于噪聲水平自適應估計的往復壓縮機振動信號的局部投影降噪方法[J].振動與沖擊,2010,29(1):53-57.

ZHANG Lai-bin, LI Feng, DUAN Li-xiang. Local projection denoise method for reciprocating compressor vibration signals based on adaptive estimation of noise level[J].Journal of Vibration and Shock, 2010,29(1):53-57.

[5] 陽建宏,徐金梧,楊德斌,等.鄰域自適應選取的局部投影非線性降噪方法[J].振動與沖擊,2006,25(4):64-67.

YANG Jian-hong,XU Jin-wu, YANG De-bin et al. Nonlinear denoise method by local projection with Adaptive neighborhood selection [J].Journal of Vibration and Shock, 2006,25(4):64-67.

[6] Miller N B. Reflections from gradual transition sound absorbers[J]. Journal of the Acoustical Society of America，1958,30(10).

[7] Cao L. Practical method for determining the minimum embedding dimension of a scalar time series. Physica D: Nonlinear Phenomena[J]. 1997,110(1-2): 43-50.

[8] Fraser A M, Swinney H L. Independent coordinates for strange attractors from mutual information[J]. Physical Review A. 1986,33(2): 1134-1140.

[9] 楊志安,王光瑞,陳式剛. 用等間距分格子法計算互信息函數(shù)確定延遲時間[J]. 計算物理,1995,12(4):442-447.

YANG Zhi-an, WANG Guang-rui, CHEN Shi-gang. Spacing divided grid method to calculate the mutual information function to determine the delay time[J].Chinese Journal of Computational Physics, 1995, 12(4):442-447.

[10] 劉希靈,李夕兵,洪亮,等.基于離散小波變換的巖石SHPB測試[J].爆炸與沖擊，2009,29(1):67-71.

LIU Xi-ling, LI Xi-bing, HONG Liang, et al. SHPB test of rock based on discrete wavelet transform[J]. Journal of Explosion and Shock Waves, 2009, 29(1): 67-71.

[11] 呂金虎,陸君安,陳士華. 混沌時間序列分析及其應用[M]. 武漢:武漢大學出版社,2005.

[12] 梁勇,孟橋,陸佶人. Lyapunov指數(shù)的算法改進與加權預測[J].聲學技術,2006,25(5):463-467.

LIANG Yong, MENG Qiao, LU Ji-ren. Improved algorithm of Lyapunov exponent and weighted prediction[J]. Technical Acoustics, 2006,25(5):463-467.

[13] Kantz H, Schreiber T. Nonlinear time series analysis[M]. Cambridge: Cambridge Univ. Press, 1977：65-74.

[14] 丁曉蓉.基于混沌理論的心音信號非線性動力學分析[D].重慶:重慶大學,2012.

[15] 姜桂仁. 混沌時序的特征量分析及相空間重構[D]. 鎮(zhèn)江:江蘇大學,2005.