蔡 俊，張翔風(fēng)

(安徽理工大學(xué) 電氣與信息工程學(xué)院，安徽 淮南 232001)

隨著社會(huì)的高速發(fā)展，人們的醫(yī)療健康意識(shí)逐漸提高，對(duì)心血管疾病的關(guān)注日益增多，心電信號(hào)中包含了關(guān)于診斷心臟疾病的大量醫(yī)學(xué)信息，所以心電信號(hào)研究成為現(xiàn)下科研中的一大熱門課題。而心電信號(hào)去噪則是信號(hào)預(yù)處理的一個(gè)關(guān)鍵步驟。心電信號(hào)是一種微弱的、非線性、不平穩(wěn)的信號(hào)，幅度一般在0.01～5mV，頻率分布大約在0.05～100Hz[1]，其中信號(hào)能量主要集中在0.25～35Hz左右[2]。心電信號(hào)在采集過程中易受到包括基線漂移、肌電干擾和50Hz工頻干擾的噪聲影響。這些干擾信號(hào)和心電信號(hào)的頻率接近，會(huì)造成心電信號(hào)的畸變，對(duì)后續(xù)心電信號(hào)特征波檢測、特征提取、分類等處理產(chǎn)生干擾，從而影響心臟疾病的正確診斷。

近年來，有許多學(xué)者提出了一些心電信號(hào)去噪算法，例如傳統(tǒng)濾波法、小波閾值去噪法、經(jīng)驗(yàn)?zāi)B(tài)分解(Empirical Mode Decomposition，EMD)、集合經(jīng)驗(yàn)?zāi)B(tài)分解(Ensemble Empirical mode Decomposition，EEMD)、互補(bǔ)集合經(jīng)驗(yàn)?zāi)B(tài)分解(Complementary Ensemble Empirical Mode Decomposition，CEEMD)。張勇[3]等利用中值濾波器濾除噪聲，提高了整體的運(yùn)算速度，但其幅頻特性和過渡帶的影響很容易造成波形失真。檀雪等[4]提出一種改進(jìn)小波閾值函數(shù)去噪方法，克服了原有的軟、硬閾值函數(shù)去噪的缺點(diǎn)，可以更好地保持心電信號(hào)的原有特征。然而，其參數(shù)選擇困難，很有可能導(dǎo)致攜帶信號(hào)信息的小波系數(shù)被錯(cuò)誤剔除，從而影響了去噪的效果。盧麗蓉等[5]采用基于心動(dòng)周期和經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)的方法來進(jìn)行心電信號(hào)降噪，雖然去噪的結(jié)果優(yōu)于普通經(jīng)驗(yàn)?zāi)B(tài)分解法，但無法解決EMD算法存在的模態(tài)混疊問題。趙偉等[6]應(yīng)用MCA-EEMD算法去除心電信號(hào)中的工頻干擾，EEMD雖然有效抑制了模態(tài)混疊，但是迭代次數(shù)和白噪聲幅值的設(shè)置，需要依賴人為經(jīng)驗(yàn)，如果數(shù)值設(shè)置不當(dāng)，那么就無法解決模態(tài)混疊問題，同時(shí)也會(huì)出現(xiàn)重構(gòu)誤差。董立超等[7]應(yīng)用CEEMD-小波包算法去除心音信號(hào)的噪聲，CEEMD算法加入噪聲集合的次數(shù)少，計(jì)算效率高，但CEEMD的不足之處分解到的固有模態(tài)函數(shù)(Intrinsic Mode Function，IMF)分量中仍然會(huì)有白噪聲殘留，而且算法自適應(yīng)差。因?yàn)镋CG信號(hào)經(jīng)過經(jīng)驗(yàn)?zāi)B(tài)分解后，噪聲集中在高頻分量中，所以傳統(tǒng)的經(jīng)驗(yàn)?zāi)B(tài)分解去噪時(shí)，會(huì)直接丟棄少量高頻分量，而這也會(huì)將高頻分量中包含的有用信息一并丟失，導(dǎo)致重構(gòu)信號(hào)的誤差過大，不利于后續(xù)的分析處理。

CEEMDAN是在EMD和EEMD基礎(chǔ)上提出的一種算法，該方法通過在分解的每一階段添加自適應(yīng)的白噪聲，計(jì)算唯一殘差信號(hào)，來得到IMF分量[8]。不僅有效地抑制了模態(tài)混疊的情況，而且分解之后得到的IMF分量中包含的殘余噪聲更少，可以降低篩選次數(shù)，減少了重構(gòu)誤差，重構(gòu)能力和完備程度大大增高。

本文針對(duì)心電信號(hào)的處理中存在的噪聲干擾問題，引入一種CEEMDAN-小波包聯(lián)合去噪辦法。通過導(dǎo)入心電信號(hào)，對(duì)其進(jìn)行CEEMDAN分解，得到IMF分量之后，繪制頻譜圖。通過相關(guān)系數(shù)的變化規(guī)律來篩選出包含噪聲的高頻IMF分量，再根據(jù)心電信號(hào)的頻率分布范圍來驗(yàn)證篩選是否正確。包含基線漂移的低頻IMF分量，則根據(jù)心電信號(hào)的特征波形來確定。對(duì)包含噪聲的高頻IMF分量進(jìn)行小波包分解，將含有基線漂移的低頻IMF分量利用中值濾波進(jìn)行噪聲消除，最后，將去除過噪聲的IMF分量和不含噪的IMF分量進(jìn)行再次重構(gòu)得到高信噪比的ECG信號(hào)。通過結(jié)合小波包分析和CEEMDAN，可以綜合兩者的優(yōu)點(diǎn)，既能更好地處理非線性非平穩(wěn)信號(hào)，噪聲信號(hào)可以被更精細(xì)地去除掉，又能在去噪后較好地保留心電信號(hào)的完整性。

1 聯(lián)合去噪算法原理

1.1 CEEMDAN分解算法

經(jīng)驗(yàn)?zāi)B(tài)分解算法是希爾伯特-黃變換(Hilbert-Huang Transform，HHT)的核心算法，可以將含噪信號(hào)分解為多個(gè)保留局部信號(hào)的數(shù)據(jù)序列，即IMF分量，適用于處理非平穩(wěn)信號(hào)[9]。完全自適應(yīng)噪聲集合經(jīng)驗(yàn)?zāi)B(tài)分解(Complete EEMD with Adaptive Noise，CEEMDAN)[10]是TORRES等提出的改進(jìn)算法。

設(shè)δi是滿足標(biāo)準(zhǔn)正態(tài)分布的高斯白噪聲信號(hào)，i=1,2,…,n，a為系數(shù)，取值為[1，2]，λ為幅值，CEEMDAN算法原理如下。

(1)

(2)

(3)用原信號(hào)序列Y(t)減去第一個(gè)本征模態(tài)分量M1(t)，構(gòu)造殘差e1(t)如式(3)所示。

e1(t)=Y(t)-M1(t)

(3)

(4)

(5)繼續(xù)減去第二個(gè)本征模態(tài)分量M2(t)，構(gòu)造殘差e2(t)如式(5)所示。

e2(t)=e1(t)-M2(t)

(5)

(6)重復(fù)上述步驟，直到殘差信號(hào)ek(t)不能再繼續(xù)分解，最終可以得到K個(gè)本征模態(tài)分量。最終得到的重構(gòu)信號(hào)Y(t)如式(6)所示。

(6)

1.2 小波包去噪

小波變換是通過挑選快速衰減的小波基，進(jìn)行伸縮、平移來表示信號(hào)的過程，常用于分析非平穩(wěn)信號(hào)[11]。一般的小波變換只對(duì)上一層的低頻信號(hào)進(jìn)行分解，而高頻細(xì)節(jié)部分則不再分析。在此基礎(chǔ)上，發(fā)展出了小波包分解，可以對(duì)小波變換中未加處理的高頻信號(hào)進(jìn)行分解，進(jìn)一步去除高頻分量中包含的信息[12]。

小波包變換的去噪步驟如下。

(1)在確定小波基和分解層數(shù)的情況下，分解原始信號(hào)。

(2)選擇合適的閾值λ，對(duì)分解得到的全部小波包系數(shù)進(jìn)行去噪，將小于該閾值的系數(shù)當(dāng)作噪聲分量加以去除。

(3)疊加重構(gòu)去噪過的小波包系數(shù)，得到去噪信號(hào)。

常見的小波去噪閾值函數(shù)主要包括硬閾值函數(shù)和軟閾值函數(shù)[13]。硬閾值函數(shù)如式(7)所示。

(7)

軟閾值函數(shù)如式(8)所示。

(8)

其中sgn(.)是符號(hào)函數(shù)，x為小波包系數(shù)，S為處理之后的小波包系數(shù)，λ為選取的閾值。

小波包分解的去噪效果會(huì)受到閾值以及閾值函數(shù)選取的影響。從均方差數(shù)據(jù)來考量，硬閾值去噪比軟閾值去噪效果好，但會(huì)使處理后的信號(hào)附加抖動(dòng)，平滑性不好[14]。軟閾值函數(shù)處理過的小波包系數(shù)，連續(xù)性較好，不存在震蕩，原始信號(hào)的平滑性被很好保留。綜合比較，本文選取軟閾值函數(shù)來對(duì)含噪ECG進(jìn)行處理。

1.3 中值濾波去噪

因?yàn)樾碾娦盘?hào)中包含的基線漂移變化緩慢、幅值較大，實(shí)驗(yàn)過程中發(fā)現(xiàn)使用小波包閾值去噪算法對(duì)其處理效果不好，所以改為用中值濾波器來進(jìn)行去噪[15]。中值濾波器在處理非平穩(wěn)信號(hào)中效果較為顯著，該方法既可以濾除噪聲，又可以保留原始信號(hào)的邊緣特征，能夠屏蔽異常值。該算法的原理是，把數(shù)字序列中一個(gè)點(diǎn)的值，用其所在的一個(gè)鄰域內(nèi)所有點(diǎn)的中值來取代，從而來消除噪聲點(diǎn)，即選擇n個(gè)采樣點(diǎn)為一個(gè)窗口，取中位數(shù)，窗口每次移動(dòng)一個(gè)采樣點(diǎn)，最后做窗口疊加運(yùn)算[16]。

1.4 聯(lián)合去噪流程

聯(lián)合去噪算法流程如圖1所示。結(jié)合CEEMDAN和小波包算法，通過導(dǎo)入含噪心電信號(hào)，對(duì)其進(jìn)行CEEMDAN分解。得到數(shù)個(gè)IMF分量后，計(jì)算各模態(tài)分量與原始心電信號(hào)的相關(guān)系數(shù)，同時(shí)通過傅里葉變換，得出各個(gè)IMF分量的頻率分布范圍。之后按照相關(guān)系數(shù)來篩選含噪IMF分量，并用頻譜圖加以核驗(yàn)。使用小波包閾值去噪和中值濾波法對(duì)包含噪聲的IMF分量進(jìn)行降噪處理，之后重構(gòu)經(jīng)過處理的IMF分量與不包含噪聲的IMF分量，得到純凈的心電信號(hào)。

圖1 聯(lián)合去噪流程圖

2 信號(hào)仿真分析

為了驗(yàn)證本文方法的有效性，生成模擬的潔凈心電信號(hào)與模擬噪聲疊加作為實(shí)驗(yàn)對(duì)象。該模擬信號(hào)結(jié)合了心電信號(hào)特征波形的特點(diǎn)，幅值不超過2.5mV。

2.1 仿真ECG信號(hào)構(gòu)建

心電圖特征波形包括P波、QRS波群、T波、U波。P波的時(shí)限為0.12s，高度為0.25mV，正常的P波與R波之間的間期在0.12～0.20s。QRS波群會(huì)由年齡不同而導(dǎo)致持續(xù)時(shí)間不同，但健康的成年人一般不會(huì)超過0.11s。S-T段在0.12～0.14s之間，電壓在0.05mV到0.1mV之間。T波的幅值范圍是0.1～0.8mV，持續(xù)時(shí)間最長為0.25s。正常U波的持續(xù)時(shí)間通常為0.16～0.25s，且幅值<0.1mV。Q-T間期在心電信號(hào)中持續(xù)時(shí)間最長，可達(dá)0.44s左右。根據(jù)上述的信息，模擬計(jì)算出心電特征值為0.25mV的P波、0.35mV的T波、0.035mV的U波、1.6mV的R波、R-R間隔為0.83s、PR間隔為0.16s的10s理想心電信號(hào)，如圖2(a)所示。

肌電干擾的頻率分布一般在5～2000Hz，屬于高頻噪聲，基線漂移一般是低頻擾動(dòng)，頻率小于1Hz，因此不易與心電信號(hào)進(jìn)行剝離[17]。在此數(shù)據(jù)基礎(chǔ)上，在模擬的純凈心電信號(hào)中添加頻率為0.15Hz，幅度為模擬信號(hào)15%的正弦波信號(hào)，作為基線漂移的模擬信號(hào)；添加頻率為50Hz的正弦信號(hào)作為工頻干擾的模擬信號(hào)；添加20dB的隨機(jī)噪聲作為肌電干擾的模擬信號(hào)，以此疊加構(gòu)成含有噪聲的心電信號(hào)，如圖2(b)所示。

(a)理想ECG信號(hào)

2.2 仿真ECG信號(hào)降噪處理

通過CEEMDAN分解理想的加噪心電信號(hào)，得到一定數(shù)量的IMF分量。如果只靠經(jīng)驗(yàn)和分量圖來挑選含噪分量，這種挑選方式較為主觀，且后續(xù)實(shí)驗(yàn)耗時(shí)耗力。相關(guān)系數(shù)表示了兩個(gè)信號(hào)之間的相似程度，數(shù)值越大，則相關(guān)性越強(qiáng)，所以引入相關(guān)系數(shù)來挑選含噪分量，可以提高算法效率，更為客觀地挑選出含噪的高頻分量。模擬心電信號(hào)被CEEMDAN分解成9個(gè)IMF分量，如圖3所示。

圖3 加噪模擬信號(hào)分解結(jié)果

可以計(jì)算出各個(gè)IMF分量與加噪ECG信號(hào)的相關(guān)系數(shù)pk[18]如式(9)所示。

(9)

ck是分解得到的IMF分量，s是加噪的理想ECG信號(hào)。

計(jì)算出的各分量的相關(guān)系數(shù)如表1所示。相關(guān)系數(shù)小于0.1的分量被定義為偽分量，即含有較少原始信息的分量，在相關(guān)系數(shù)大于0.1的分量當(dāng)中，出現(xiàn)的第一個(gè)局部極大值之前的IMF分量被定義為噪聲主導(dǎo)分量[19]。由表1所示，IMF1的相關(guān)系數(shù)為0.0451小于0.1，所以是偽分量，IMF4為局部最大值，相關(guān)系數(shù)為0.7540，在IMF4之前的模態(tài)分量即噪聲主導(dǎo)分量，包含大量高頻噪聲，所以可以得出IMF1～I(xiàn)MF3包含高頻噪聲。由圖3可以得出，添加的干擾信號(hào)主要集中在IMF1～I(xiàn)MF3分量中，低頻分量IMF7～I(xiàn)MF9的頻率逐漸減小，看不出明顯的心電信號(hào)特征波，主要包含基線漂移。

表1 模擬信號(hào)各分量的相關(guān)系數(shù)

心電信號(hào)分解所用的小波基，不僅要考慮其正交性、緊支撐性、光滑性，同時(shí)也要考慮分解層數(shù)。分解層數(shù)過多，也會(huì)導(dǎo)致信號(hào)中的原始信息被濾除。經(jīng)過多次實(shí)驗(yàn)，最終選擇coif5小波基對(duì)IMF1～I(xiàn)MF3分量進(jìn)行三層小波包分解，之后對(duì)分解得到的小波包系數(shù)進(jìn)行軟閾值去噪，信號(hào)重構(gòu)，得到去噪后的高頻IMF分量。實(shí)驗(yàn)采用中值濾波算法消除IMF7～I(xiàn)MF9的基線漂移，最后重構(gòu)疊加降噪后的分量與不含噪分量，組成降噪后的心電信號(hào)。

圖4為加入模擬噪聲的理想ECG信號(hào)去噪結(jié)果。由圖2(b)可以看出，原始理想ECG信號(hào)中加入噪聲干擾之后，圖形出現(xiàn)許多毛刺，這些干擾掩蓋了理想ECG信號(hào)中的有用的生理信息。由圖4可知，含噪ECG信號(hào)經(jīng)過CEEMDAN-小波包聯(lián)合去噪之后，不僅保留了原始的生理特征信息，同時(shí)剔除掉了多余的噪聲干擾，使得信號(hào)曲線更為光滑，具有較好的去噪效果。

圖4 理想ECG信號(hào)降噪結(jié)果

3 實(shí)驗(yàn)分析

3.1 數(shù)據(jù)處理

本文通過對(duì)MIT-BIH心律異常數(shù)據(jù)庫的118號(hào)數(shù)據(jù)進(jìn)行分析，來進(jìn)一步對(duì)本文算法進(jìn)行有效性驗(yàn)證。實(shí)驗(yàn)樣本采用含有采集過程中的基線漂移、肌電干擾以及工頻干擾信號(hào)的118號(hào)數(shù)據(jù)，采樣點(diǎn)數(shù)為3600個(gè)，采樣頻率為360Hz。對(duì)118號(hào)數(shù)據(jù)經(jīng)過CEEMDAN分解得到12個(gè)IMF分量。對(duì)12個(gè)IMF分量進(jìn)行傅里葉變換得到的頻譜圖，如圖5所示。

圖5 118號(hào)數(shù)據(jù)本征模態(tài)分量的頻譜圖

心電118號(hào)數(shù)據(jù)各模態(tài)分量的相關(guān)系數(shù)如表2所示。由表2可以得出，相關(guān)系數(shù)的變化趨勢為先增高再降低，IMF5的相關(guān)系數(shù)最大，為0.6548。根據(jù)上述提到的定義，所以IMF1～I(xiàn)MF4為含噪的高頻分量，主要包含肌電干擾和工頻干擾，需要進(jìn)行小波包閾值去噪。

表2 118號(hào)ECG數(shù)據(jù)各分量相關(guān)系數(shù)

觀察圖5中的前4個(gè)IMF分量可以看出，分量的頻率超過35Hz，超出了大部分心電信號(hào)的能量分布范圍，再次驗(yàn)證IMF1～I(xiàn)MF4是包含噪聲的高頻分量。IMF9～I(xiàn)MF12分量的心電信號(hào)的特征波形不明顯，主要是正弦波形的疊加，是包括基線漂移干擾的低頻分量，所以對(duì)其進(jìn)行中值濾波處理。

之后對(duì)118號(hào)ECG數(shù)據(jù)分別進(jìn)行小波包去噪、EMD-小波包去噪、EEMD-小波包去噪、CEEMD-小波包去噪這四種不同算法進(jìn)行去噪。得到結(jié)果如圖6所示。

(a)118號(hào)ECG數(shù)據(jù)

由圖6可以看出五種算法所得到的結(jié)果，皆濾除了一部分噪聲。但是圖6(c)小波包去噪得到的波形存在較大失真，屬于過度濾除，很可能會(huì)丟失重要的醫(yī)學(xué)信息，例如QRS波群和T、P波的幅值、間隔以及頻率信息會(huì)因此而改變，整體去噪效果不佳。傳統(tǒng)小波包去噪算法在重構(gòu)時(shí)，根據(jù)頻率分段來選擇丟棄幾個(gè)節(jié)點(diǎn)系數(shù)，雖然圖形效果較好，但是被丟棄的節(jié)點(diǎn)系數(shù)中可能包含有用的高頻信息，這對(duì)后續(xù)心電信號(hào)的處理是不利的。

EMD-小波包算法、EEMD-小波包算法以及CEEMD-小波包算法，在去噪處理之后，局部依然保留了許多毛刺，這會(huì)導(dǎo)致有效的生理信息被淹沒。而圖6(b)所示，經(jīng)過本文算法所得到的心電信號(hào)波形則濾除了更多的噪聲，最大限度地保持了信號(hào)的原始信息，圖形更為明晰，去噪效果更佳。

3.2 去噪算法評(píng)估

為了驗(yàn)證本文算法在處理心電信號(hào)上的有效性，采用信噪比(Signal-Noise ratio，SNR)和均方根誤差(Root Mean Square Error，RMSE)兩個(gè)指標(biāo)來進(jìn)行評(píng)估[20]。信噪比是指信號(hào)中有效信號(hào)和噪聲的比值，數(shù)值越大越好。均方根誤差是用來描述經(jīng)過去噪重構(gòu)之后的信號(hào)和原信號(hào)之間的誤差，數(shù)值越小越好。信噪比和均方根誤差的定義如式(10)和式(11)所示。

(10)

(11)

其中，xi是原始信號(hào)，yi是去噪之后的信號(hào)。

除了118號(hào)數(shù)據(jù)，本文還引用了MIT-BIH心律失常數(shù)據(jù)庫的105、109號(hào)ECG數(shù)據(jù)，這三組數(shù)據(jù)都包含著多種噪聲干擾。實(shí)驗(yàn)分別應(yīng)用本文算法以及前人應(yīng)用過的4種算法去噪，對(duì)算法之間的指標(biāo)進(jìn)行比較。由表3可以看出，處理105號(hào)數(shù)據(jù)時(shí)，小波包算法、EMD-小波包算法、EEMD-小波包算法、CEEMD-小波包算法的信噪比分別為22.6010、23.2240、27.4238、27.4305，呈現(xiàn)遞增趨勢，均方根誤差分別為0.0181、0.0163、0.0104、0.0103，呈現(xiàn)遞減趨勢，證明聯(lián)合算法的降噪效果要優(yōu)于小波包算法，EEMD和CEEMD算法優(yōu)于EMD算法。而本文算法的信噪比達(dá)到了31.6446，均方根誤差降為了0.0096，將這幾組數(shù)據(jù)進(jìn)行對(duì)比，可以得出，本文算法的信噪比更大，均方根誤差更小，有更好的信噪分離效果。

表3 不同信號(hào)應(yīng)用不同算法去噪的指標(biāo)對(duì)比表

4 結(jié)論

針對(duì)現(xiàn)有心電信號(hào)去噪算法所存在的不足，在CEEMDAN算法的基礎(chǔ)上，提出了一種CEEMDAN-小波包聯(lián)合去噪的算法。為了驗(yàn)證算法的可行性，本文首先根據(jù)一個(gè)完整心電信號(hào)周期所包含的幅值、頻率、間隔等信息，生成純凈的模擬心電信號(hào)。在此基礎(chǔ)上加入噪聲信號(hào)，得到含噪心電信號(hào)，對(duì)此仿真信號(hào)應(yīng)用本文算法進(jìn)行去噪，得到的去噪效果較好。之后引入MIT-BIH心律失常數(shù)據(jù)庫的105、109、118號(hào)數(shù)據(jù)，并使用前人應(yīng)用過的四種算法和本文算法去噪，來進(jìn)行算法之間的對(duì)比。

該算法利用CEEMDAN算法解決了EMD算法存在的模態(tài)混疊、EEMD算法重構(gòu)誤差大、CEEMD算法自適應(yīng)性差等問題，同時(shí)利用小波包去除了心電信號(hào)包含的多種噪聲干擾。從波形圖結(jié)果可以得出，通過CEEMDAN-小波包算法可以在消除噪聲分量的同時(shí)，有效保留了原始信號(hào)中的信息。而且該算法與小波包單獨(dú)降噪、EMD-小波包聯(lián)合降噪、EEMD-小波包聯(lián)合降、CEEMD-小波包聯(lián)合降噪方法相比，得到的信噪比最大，均方根誤差最小，證明了本文算法優(yōu)于其他算法，而且為以后的心電信號(hào)預(yù)處理提供了新思路，具有比較廣闊的應(yīng)用前景。