基于改進經(jīng)驗模態(tài)分解域內(nèi)心動物理特征識別模式分量的心電信號重建*

牛曉東 盧莉蓉 王鑒 韓星程 郭樹言 王黎明?

1) (中北大學,信息探測與處理山西省重點實驗室,太原 030051)

2) (長治醫(yī)學院物理教研室,長治 046000)

3) (長治醫(yī)學院生物醫(yī)學工程系,長治 046000)

心電圖(electrocardiogram,ECG)診斷心臟疾病的嚴格標準,要求有效地消除噪聲并準確地重建ECG信號.經(jīng)驗模式分解(empirical mode decomposition,EMD)方法重建ECG 信號中,模式混疊及重建采用模式分量的識別以經(jīng)驗為基礎,導致重建ECG 信號準確度降低,且方法不具有自適應和通用性.本文首先基于積分均值定理提出一種改進的EMD 方法—積分均值模式分解(integral mean mode decomposition,IMMD)方法,經(jīng)5000 個高斯白噪聲樣本的蒙特卡羅法驗證,IMMD 方法比EMD 具有更優(yōu)多分辨率分析能力,能夠有效地緩解模式混疊.其次,基于ECG 信號內(nèi)固有心動物理特征量識別重建ECG 信號所采用的模式分量,具有現(xiàn)實物理意義,因此,方法具有自適應和通用性.經(jīng)驗證,提出方法重建47 例ECG 信號與原ECG 信號的相關系數(shù)中: 31 例優(yōu)于變分模式分解方法; 33 例優(yōu)于Haar 小波軟閾值法; 42 例優(yōu)于集總經(jīng)驗模式分解方法; 45 例優(yōu)于EMD 方法.相關系數(shù)均值為0.8904,方差為0.0071,表現(xiàn)穩(wěn)定且最優(yōu).

1 引 言

心電圖(electrocardiogram,ECG)記錄的是心臟生理電活動信號,ECG 檢查一般通過ECG 中P,QRS,T 特征波波形、時長及幅值識別診斷各種心臟疾病.例如P 波形狀高尖,且振幅 ≥ 0.25 mV判斷為右房肥大; P,R 間期延長 > 0.21 s 存在房室傳導阻滯等[1].因此,ECG 信號輕微失真,也可能導致心臟疾病誤診.然而,ECG 信號是一種非線性、非平穩(wěn)弱信號,極易被噪聲污染.噪聲一般包括兩類: 50 Hz (在一些國家為60 Hz)電源線工頻及由呼吸、電極或人體運動引起的基線漂移(baseline wander,BW)等窄帶噪聲; 主要由肌肉收縮引起的肌電干擾(muscle artifact,MA)等寬帶噪聲[2].兩種噪聲常常疊加引起ECG 信號畸變,嚴重影響ECG 診斷評估心臟病準確性,甚至導致診斷失效.嚴格的ECG 檢查標準量以及ECG 信號的噪聲特點,要求有效消除噪聲后重建的ECG 信號具有高準確度.

窄帶噪聲頻譜與ECG 信號頻譜不重疊,常見線性和非線性信號處理方法都可有效校正[3?5],而寬帶噪聲全頻域污染ECG 信號,導致嚴重特征波失真,成為ECG 信號消噪重點.迄今為止,已有多種非線性信號處理方法應用于消除寬帶噪聲并重建ECG 信號,但各種方法均存在缺點,重建準確度不令人滿意.其中,自適應濾波需要額外地采集噪聲相關參考信號[6]; 奇異值分解(singular value decomposition,SVD)[7]和小波閾值法[8,9]中,閾值過小消噪不理想,過大將導致波形嚴重畸變,并且閾值不能自適應噪聲水平的波動[10].

經(jīng)驗模式分解(empirical mode decomposition,EMD)是一種自適應且數(shù)據(jù)驅動的非線性信號處理方法[11],可以克服SVD 和小波閾值法缺乏自適應的缺點,成為近年來生理信號[12]、醫(yī)學圖像[13]、特別是ECG 信號重建的研究熱點.然而,EMD 域內(nèi)重建ECG 信號存在兩個問題: 1)瞬時模態(tài)分量(intrinsic mode function,IMF)之間的模式混疊,導致ECG 信號重建不準確.集總經(jīng)驗模式分解(ensemble empirical mode decomposition,EEMD)[14]以及變分模式分解(variational mode decomposition,VMD)[15]能夠緩解模式混疊,并已應用于ECG 信號消噪和重建[16,17].然而,嚴格意義上,對每一個ECG 信號,EEMD 都需要預先多次測試并選擇合適的輔助噪聲幅度和數(shù)量;VMD 同樣需要預先多次測試并選擇合適的分解層數(shù)、平衡參數(shù)、噪聲松弛等參數(shù),因此,兩種方法都不屬于自適應且不具有通用性.2) 重建ECG 信號采用的模式分量基于工作者經(jīng)驗識別[18,19],不具有通用性.采用IMF 的香農(nóng)熵、平均值、方差等統(tǒng)計特征量的閾值識別IMF,閾值也基于經(jīng)驗設置[20].

極值域模式分解(extremum field mean mode decomposition,EMMD)[21]是一種改進的EMD 方法.EMMD 采 用 積 分 均 值 定 理(integral mean value theorem,IMVT)求信號均值,能夠增強分解信號能力,但EMMD 算法分解信號存在過分解問題.為了在EMD 域內(nèi)進一步提高重建ECG 信號準確度,同時保證重建方法具有自適應和通用性,本文提出一種重建ECG 信號新方法.新方法首先改進EMMD,保留IMVT 均值方法并解決算法導致的過分解問題.為了區(qū)別于EMMD,本文將改進的EMMD 稱為積分均值模式分解(integral mean mode decomposition,IMMD).經(jīng)驗證可以發(fā)現(xiàn),IMMD 比EMD 具有更優(yōu)的多分辨率分解能力,并且可以有效地緩解IMF 之間的模式混疊.其次,心電活動或心肌細胞的除極、復極過程對應一個完整的心動周期,因此,ECG 信號呈現(xiàn)心動的周期性波動,其IMF 分量同樣表現(xiàn)出心動物理周期特征.本文采用心動物理周期特征識別屬于ECG信號的IMF,識別方法符合ECG 信號的物理本質(zhì),并保證重建ECG 信號方法的自適應及通用性.

2 方法理論

2.1 IMMD 方法

設xi(t)(ti≤t≤ti+1,i=1,2,··· ,n ?1) 為信號x(t) 相鄰兩個極值x(ti)和x(ti+1) 之間的局部,IMVT 求信號局部xi(t) 的均值為

mi(t)固定于局部中點.基于(1)式求出所有局部均值點后,三次樣條插值所有均值點構建均值曲線m(t).

原型模式函數(shù)(proto mode function,PMF)[22]即為

上述過程稱為一次均值篩選[11],PMF 記為PMF1.PMF1作為一個新信號,迭代重復均值篩選過程k(k= 2,3,···)次,得到PMFk.當PMFk滿足柯西篩選停止準則[11]:

篩選停止,PMFk即為IMF1.將剩余分量r1(t)=x(t) – IMF1 作為新信號,重復上述所有過程得到IMF2.同理,可得其余IMF 分量及剩余分量r(t).

除局部均值外,EMMD 方法利用相鄰兩個局部均值mi(t)和mi+1(t)加權求取二者之間極值點ti+1處均值,即:

其中,

然而,m(ti+1)攜帶數(shù)據(jù)信息是mi(t)和mi+1(t)所攜信息的重復.并且均值點數(shù)增加,極易導致樣條插值構造均值曲線的不平滑,最終造成信號過分解.因此,IMMD 方法僅采用局部均值經(jīng)三次樣條插值構建均值曲線.

2.2 IMMD 方法的多分辨率分解信號能力

文獻[23]研究了EMD 分解5000 個長度N=512,均值為0,方差為1 的高斯白噪聲數(shù)據(jù)樣本,得到如下結論: EMD 類似于小波變換Mallat 算法,等效于一個恒定品質(zhì)因子Q的二分(或二進)濾波器組,具有多分辨率分解信號能力.本節(jié)將采用同樣方法,分析IMMD 分解高斯白噪聲特性,以此說明方法分解信號的性能.

IMMD 分解每個高斯白噪聲樣本,得到至少10 個IMF.對所有樣本相應階數(shù)IMF 求平均傅里葉功率譜,結果見圖1(a).IMF1 相當于一個高通濾波器,其余IMF 等效于一組重疊帶通濾波器,且后一IMF 等效帶通濾波器中心頻率約是前一 IMF等效帶通濾波器中心頻率的2/3.IMMD 方法分解高斯白噪聲等效于濾波器組,且具有三分特性.

除IMF1 外,后一IMF 等效帶通濾波器頻帶寬度大致是前一IMF 的2/3,所以等效帶通濾波器組具有恒Q性質(zhì).因此,IMFs 的功率譜存在自相似性[23]:

圖1 IMMD 和EMD 方法分解高斯白噪聲的等效濾波器組特性 (a) IMMD (實線)和EMD (虛線)的IMF 分量的平 均 功 率 譜; (b) 基 于(5)式,IMMD (實 線)和EMD (虛線)的IMF 分量的平均功率譜坍縮重合Fig.1.Equivalent filter banks of IMMD and EMD decomposing Gauss white noises: (a) Averaged power spectra of IMFs from IMMD (solid curves) and EMD (dotted curves);(b) collapse and coincidence of the average power spectrum of IMFs from IMMD (solid curves) and EMD (dotted curves) based on Eq.(5).

其中Sk′(f) 是第k′個IMF 的 功率 譜(Sk(f) 是 第k個IMF 的功率譜),k′ >k≥2 ;ρ是常數(shù),對于IMMD,ρ=3/2 ,對于EMD,ρ=2 .基于(5)式對平均功率譜標準化,所有帶通IMF 的功率譜坍縮并重合為一條曲線.對于EMD,帶通濾波器混疊較多低頻成分,關于中心頻率對稱度差,窄帶特性差,需要適當調(diào)整ρ值為2.01[24].對于IMMD,等效帶通濾波器關于中心頻率對稱度及窄帶特性明顯優(yōu)于EMD,兩套等效帶通濾波器組的平均功率譜標準化結果如圖1(b)所示.

IMMD 分解高斯白噪聲等效濾波器組具有三分特性,等效濾波器個數(shù)多于EMD (圖1),濾波器組窄帶特性及恒Q性質(zhì)優(yōu)于EMD,因此,具有更好的多分辨率分析能力.除上述5000 個高斯白噪聲樣本數(shù)據(jù)的蒙特卡羅方法實驗外,本課題組基于大量數(shù)據(jù)實驗,驗證得到: 相比EMD,IMMD 具有更好的多分辨率分析能力,可以有效地緩解IMF分量之間模式混疊.

2.3 心動物理特征識別IMF

ECG 信號是記錄人體心臟電活動信號,具有心動物理周期特性,其IMF 分量將呈現(xiàn)心動周期或心率(heart rate,HR)特征.

任一個IMF 分量可表達為[25]

其中φ(t) 為IMF 頻率調(diào)制,a(t) 為IMF 振幅調(diào)制.心動周期特征通過兩種模式作用于IMF 分量:1) 屬于ECG 信號的低階IMF,心動周期特征作用于調(diào)幅a(t) .此時,IMF 分量的包絡具有心動周期特性,包絡的頻譜中最大幅值對應頻率等于HR.2) 屬于ECG 信號的高階IMF,心動周期特征作用于調(diào)頻φ(t) .此時,IMF 具有周期性心動的諧波特性,頻譜中最大幅值對應頻率為HR 的整數(shù)倍(通常為1—3 倍).噪聲屬性IMF 分量不能表現(xiàn)上述兩種心動周期特性.

2.4 IMMD 域內(nèi)心動物理特征識別模式分量重建ECG 方法

1) 含噪ECG 信號經(jīng)IMMD 分解為一套IMFs分量集合;

2) 求出每一個IMF 的上包絡功率譜,若譜中最大幅值對應頻率等于HR,則IMF 被識別為ECG信號分量;

3) 對步驟2 中未被識別為ECG 信號分量的IMF,求出每一個IMF 功率譜,若譜中最大幅值對應頻率等于HR 整數(shù)倍,則IMF 也被識別為ECG信號分量;

4) 采用上述所有被識別為ECG 信號分量的IMF重建ECG 信號.

3 實驗驗證

3.1 105 ECG 信號[26]實驗

選取MIT-BIH 心律失常數(shù)據(jù)庫中105 ECG、壓力測試數(shù)據(jù)庫中BW 和MA 噪聲[27],疊加構成含噪105 ECG 信號,采樣頻率為360 Hz,采樣長度為3.6×103,HR 等于1.4 Hz,如圖2 所示.可見,原105 ECG 信號的噪聲污染程度較低,便于本節(jié)下述實驗中的定量分析.

圖3 含噪105 ECG 信號的IMF (藍色曲線)及其包絡(紅色曲線)Fig.3.IMF envelopes (red curves) and IMFs (blue curves) of noisy No.105 ECG signal.

圖2 實驗采用信號波形 (a) 105 ECG 信號; (b) BW 信號; (c) MA 信號; (d)含噪105 ECG 信號Fig.2.Waveform of signals used in experiment: (a) No.105 ECG signal; (b) BW signal; (c) MA signal; (d) the noisy No.105 ECG signal.

圖4 含噪105 ECG 信號IMFs 的功率譜(藍色曲線),及IMFs 包絡的功率譜(紅色曲線)Fig.4.Power spectra of envelopes of IMFs (red curves) and of IMFs (blue curves) of noisy No.105 ECG signal.

IMMD 分解含噪105 ECG 信號得到IMF 分量及其包絡如圖3 所示,對應頻譜如圖4 所示.IMF1和IMF2 及其包絡的頻域中,幅值最高點處頻率不等于1.4 Hz,因此沒有表現(xiàn)出心動周期特征作用IMF 的兩種模式(詳見2.3 節(jié)),被識別為MA 噪聲高頻分量; IMF3—8 包絡的頻域中,幅值最高點對應頻率等于心率1.4 Hz,符合心動周期特征作用IMF 的第一種模式; IMF9—12 的頻域中,幅值最高點對應頻率分別等于心率的3,3,2,1 倍,識別為ECG 信號諧波分量,符合心動周期特征作用IMF 的第二種模式; IMF13—r 沒有表現(xiàn)出心動周期特征作用IMF 的兩種模式,統(tǒng)一歸類為低頻噪聲.產(chǎn)生的原因包括: MA 噪聲低頻分量,BW 噪聲分量,IMMD 方法端點效應以及原105 ECG 信號的非零均值.經(jīng)上述心動物理周期特征識別,IMF3—12 為ECG 信號分量,IMF1 和IMF2 為高頻噪聲,IMF13—r 為低頻噪聲.

為了顯示所提出方法的能力,將本文方法同近年來常用VMD、小波軟閾值法、EEMD 以及EMD重建105 ECG 信號進行對比,利用這五種方法重建的105 ECG 信號如圖5 所示.由圖5 可見,五種方法都很好地消除了BW 噪聲.由于MA 噪聲寬頻特性,五種方法重建ECG 信號中仍然存在少量MA 噪聲分量,特別是Harr 小波軟閾值方法.EMD方法重建ECG 信號畸變最為嚴重,其次是EEMD.實際中,EEMD 對含噪105 ECG 信號的多次重建結果之間都有輕微不同,這是由輔助白噪聲的隨機性引起的[28].另外,VMD 方法中特征R 波的峰值失真比本文提出方法嚴重,例如,圖5(b)中第10個R 波波峰峰值損失14.5%,本文方法為3.2%.

采用重建ECG 與原ECG 信號的相關系數(shù)R定量描述重建準確度,信噪比(signal-to-noise ratio,SNR)及均方誤差(mean square error,MSE)定量描述消噪能力,五種方法重建105 ECG 信號R,SNR,MSE 的值如表1 所列.可 見,提出方法的R 和SNR 值最大,MSE 值最小,重建105 ECG 信號和消噪能力最優(yōu).

3.2 47 例含噪ECG 信號實驗

為了進一步驗證提出方法的有效性以及普遍適用性,選取MIT-BIH 心律失常數(shù)據(jù)庫中其余ECG 信號(心律失常ECG 信號對應類型見表2),分別疊加BW 和MA 噪聲,構成47 例含噪ECG信號(除去周期性極差的232 號),并采用上述五種方法重建47 例原ECG 信號.

圖5 原105 ECG 信號(藍色點虛線)與由五種方法重建的105 ECG 信號(紅色實線) (a) 本文方法; (b) VMD; (c) Haar 小波軟閾值; (d) EMMD; (e) EMDFig.5.Original No.105 ECG signal (blue dotted curves) and the No.105 ECG signals (red solid curves) reconstructed by 5 methods: (a) The proposed method; (b) VMD; (c) Haar wavelet with soft threshold; (d) EEMD; (e) EMD.

表1 五種方法重建105 ECG 信號的特征量值Table 1.Characteristic values of No.105 ECG signals reconstructed by 5 methods.

表2 實驗采用的部分ECG 信號對應心律失常類型Table 2.Type of arrhythmia corresponding to some ECG signals used in the experiment.

由于大多數(shù)原ECG 信號本身含有一定量的低頻和高頻噪聲,嚴重影響SNR 和MSE 值的準確度,因此,僅采用相關系數(shù)R對五種方法進行對比評估.利用五種方法重建5 組47 例ECG 信號與原ECG 信號的相關系數(shù)柱狀圖如圖6 所示.其中,本文提出方法有31 例R值優(yōu)于VMD,33 例R值優(yōu)于Haar 小波,42 例R值優(yōu)于EEMD,45 例R值優(yōu)于EMD.

采用五種方法重建47 個ECG 信號的相關系數(shù)R值的統(tǒng)計盒形圖如圖7 所示,平均值及方差見表3.表3 中,本文提出方法對應相關系數(shù)均值最高,方差僅大于小波閾值.由圖7 可見,本文提出方法明顯優(yōu)于其余四種方法.本文提出方法重建ECG 信號最為穩(wěn)定、準確,其次為VMD,Haar小波閾值法,EEMD,EMD.

表3 五種方法重建47 個ECG 信號的R 值的均值與方差Table 3.Means and variances of R values of 47 ECG signals reconstructed by 5 methods.

4 討論與結論

相比經(jīng)驗性包絡相減法求均值(EMD,EEMD方法內(nèi)),本文提出的IMMD 方法采用IMVT 求均值,具有嚴格的均值數(shù)學基礎,因此,IMF 能夠從信號中被準確地篩選剝離.經(jīng)5000 個高斯白噪聲樣本的蒙特卡羅法統(tǒng)計驗證,IMMD 方法具有比EMD,EEMD (同EMD 方 法 一 樣,EEMD 具 有恒Q二分濾波器組特性[29])更優(yōu)的多分辨率分解信號能力,可以有效緩解IMF 之間的模式混疊.

圖6 采用五種方法重建47 個ECG 信號的R 值柱狀圖Fig.6.Histogram of R values of 47 ECG signals reconstructed by 5 methods.

圖7 采用五種方法重建47個ECG信號的相關系數(shù)R 值的統(tǒng)計盒子圖Fig.7.Five box-plots of R values of 47 ECG signals reconstructed by 5 methods.

另外,ECG 信號具有心動周期、心率等波動物理特征,ECG 信號的IMF 分量同樣具有心動周期或HR 的特性.因此,本文方法中使用心動周期或HR 識別ECG 信號的IMF 分量,符合ECG 的物理本質(zhì)特性.EMD,EEMD 和VMD 域內(nèi)ECG 重建一般采用QRS 特征波經(jīng)驗性識別IMF,適用于時域中具有顯著QRS 波動的IMF 識別(例如圖3中IMF4).如果存在: 1) 時域中QRS 波動不顯著但屬于ECG 信號的低階IMF (例如圖3 中IMF3);2) 時域完全沒有QRS 特征但屬于ECG 信號的高階IMF (例如圖3 中IMF12),該方法識別錯誤.小波閾值法中固定閾值不能自適應小波系數(shù)不同局部內(nèi)噪聲水平波動,因此,重建ECG 信號可能存在局部消噪不理想(例如圖5(c)).所以,本文提出方法,能夠比上述四種方法進一步提高ECG信號重建準確度,且方法具有自適應和通用性.經(jīng)47 例ECG 信號實驗驗證,提出的方法重建ECG與原ECG 信號相似度平均值為0.8904,方差為0.0071,且有31 例相似度值優(yōu)于VMD,33 例相似度值優(yōu)于Haar 小波,42 例相似度值優(yōu)于EEMD,45 例相似度值優(yōu)于EMD.本文提出方法重建ECG 信號能力體現(xiàn)了物理特征或現(xiàn)象在生物電信號處理中具有重要作用.

本方法中,建議通過含噪ECG 信號IMF 分量包絡頻譜圖,經(jīng)驗確定HR.HR 可以通過嚴格的RR 波間隔得到,但需要大量算法實現(xiàn).實際上,人心率一般為1—1.7 Hz (60—100 次/min).對于采集良好的ECG 信號,其IMF 分量包絡的頻譜圖中,1—2 Hz 之間的幅度最大值對應頻率即為HR,并且它常常出現(xiàn)在多個IMF 分量包絡的頻譜中,容易辨別(如圖4 所示).本文選取的47 個ECG信號,都是基于上述經(jīng)驗方法確定HR.

提出的方法具有一定局限.對于極其特殊的232 ECG 信號(竇性心動過緩、一級房室傳導阻滯和頻繁異位心房運動,ECG 信號停頓持續(xù)長達6 s)重建,本文提出方法失效.關于心動周期特性極差的ECG 信號重建,以及從物理現(xiàn)象本質(zhì)探索ECG信號處理,是本課題組下一步的工作.