基于集合經驗模態分解和信號結構分析的心電信號R波識別算法

林金朝 李必祿 李國權* 黃正文 龐 宇

①(重慶郵電大學通信與信息工程學院 重慶 400065)

②(布魯內爾大學電子與計算機工程系 倫敦 UB8 3PH)

③(光電信息感測與傳輸技術重點實驗室 重慶 400065)

1 引言

心電信號(ElectroCardioGram,ECG)實時記錄著心臟的活動狀況，攜帶著豐富的心律以及病理信息[1]。心電信號由波、段、間期組成，各自攜帶著對應的臨床信息，是心血管疾病診斷的重要依據。正常的心電信號由P波，QRS復合波，T波組成，不同的波對應著心臟不同的電活動過程[2]。R波是QRS復合波中特征最明顯的波形，是確定心電信號各波段的重要參考，是ECG自動分析的重要前提[3]。

針對心電信號中R波的識別，硬件方式主要通過峰值電壓檢測器結合電壓比較器和單穩電路對R波進行檢測[4]，但會受到器件特性等因素的影響。目前文獻以及應用主要采用數字信號處理等軟件的方式進行R波的檢測，檢測思路主要有實時檢測方法和模型匹配方法。實時檢測方法根據心電信號中R波峰值和斜率的特征對R波進行檢測；模型檢測方法通過學習的方式獲取R波的模板，然后計算模板與心電信號的相關性來實現R波的檢測。Pan等人[5,6]提出一種基于導數濾波器的R波檢測算法。該算法利用導數濾波器對帶通濾波后的心電信號進行微分以獲取R波的斜率信息，對得到的斜率信息通過自適應雙閾值的方法以實現R波的識別。該算法是一種通過閾值進行處理并且具有低復雜度的實時識別方法。由于特征波頻率的變化，基于濾波的方法檢測性能受到較大影響[7]。吳建等人[8]提出一種基于差分閾值和模板匹配的方法對心電信號的R波進行識別。該方法利用1階差分獲取R波初始模板，通過滑動窗口的方式將初始模板與窗口內的信號進行匹配，尋找相似度最大的匹配信號來定位R波，成功完成一次R波檢測后更新模板繼續下一次匹配。該方法閾值的設定取決于上一次檢測到的R波，當相鄰R波出現幅值相差較大的情況，漏檢現象比較明顯。Merah等人[9]利用平穩小波變換對心電信號進行多尺度分解，然后對各細節分量進行能量，頻率以及相關性分析，根據選取出的合適細節分量的局部極值點信息檢測R波。孫亞楠等人[10]利用改進的小波閾值法自適應地在恰當的頻率子帶上提取出R波的候選集，然后根據RR間期的局部變化趨勢對R波候選集進行篩選進而實現R波的識別。季虎等人[11,12]對心電信號采取小波變換進行多尺度分解，通過模極大值的方式在分解尺度上完成R波的識別。但基于小波變換的方法太過依賴小波基的選擇，小波基函數影響著各細節分量中的能量分布，從而影響R波識別的準確率[13]。行鴻彥等人[14]利用經驗模態分解(Empirical Mode Decomposition,EMD)將心電信號分解成一系列本征模態函數(Intrinsic Mode Function,IMF)，然后結合軟閾值的方法進行預處理，再利用模極大值和R波特征點對應的關系實現R波的檢測。但EMD存在的模態混疊問題影響R波對應特征點的識別。文獻[15]提出一種基于有限狀態機的自適應閾值的R波檢測算法。該算法在預處理階段去除噪聲的同時初始化R峰值閾值，然后利用有限狀態機根據信號的變化趨勢和之前檢測到的R峰值對閾值自適應修正，根據自適應的閾值判定R波的位置。但該方法需要進行復雜的計算[16]。文獻[17]將人工神經網絡(Artificial Neural Network,ANN)引入心電R波的識別。該方法提出一種基于sigmoidal徑向基函數優化非線性自適應白化濾波器對心電信號進行預處理以抑制噪聲同時增強QRS波，然后根據決策邏輯確定R波的位置。R波識別準確率可達99.91%。但ANN在訓練階段需要大量心電信號的先驗信息，運算量大，內存資源消耗高，難以用于實時檢測[13,18]。然而，在心電信號在采集和傳輸的過程中受到工頻干擾，肌電干擾以及基線漂移的影響，掩蓋了QRS復合波的部分特征，嚴重影響了R波的定位精度。大部分方法在對R波進行識別之前，都會對心電信號進行預處理以排除噪聲的干擾，排除噪聲干擾的同時可能會丟失心電信號的部分有用信息。因此，Safari等人[19]提出一種基于集合經驗模態分解(Ensemble Empirical Mode Decomposition,EEMD)和獨立成分分析(Independent Component Analysis,ICA)的方法對帶噪聲的心電信號直接進行R波的識別，但識別的靈敏度和準確率有待提高。

本文針對帶噪聲的心電信號提出一種基于EEMD和信號結構分析高準確率的R波識別算法。首先利用EEMD將帶噪聲的心電信號分解成一系列本征模態分量，然后對分解后的各模態分量作ICA分析以提取出R波特征最明顯的成分；最后對R波最明顯的成分進行結構分析，從而實現對R波的準確定位。

2 集合經驗模態分解

為了解決經驗模態分解(Emp ir ical Mod e Decomp osition,EMD)存在的模態混疊問題，Wu等人[20]在原始信號上加入若干次輔助白噪聲，把每一次構造的信噪混合體進行EMD分解，最后對各模態分量取平均值獲取逼近的真實模態。EEMD具體步驟為

(1)在原始信號x(t)上分別加入N次不同大小的輔助白噪聲構成不同的信噪混合體

其中，Xi(t)(i=1,2,···,N)為第i次加入輔助白噪聲后的信號，k i為 第i次 加入的白噪聲與原始信號x(t)幅值標準差比值，σx為 原始信號標準差，n(t)為歸一化白噪聲。

(2)對 信 噪 混 合 體X i(t)(i=1,2,···,N)進 行EMD分解，得到m個IMF分量cij(t)(j=1,2,···,M)和一個殘余量ri(t)，可以表示為

(3)對每次EMD分解對應的IMF分量和殘余量求平均值

3 本文算法

針對心電信號R波的識別，由于采集和傳輸過程中各種噪聲的干擾，大部分算法都會在R波識別之前進行預處理，排除各種噪聲對R波識別的干擾。在預處理過程中，各種預處理算法都會在一定程度上破壞心電信號的有用成分，并且預處理過程會增加整個R波識別的處理時間，難以應用于實時處理。針對帶噪聲的心電信號，本文提出一種無需預處理直接對R波進行識別的算法。首先利用EEMD將帶噪聲的心電信號分解成一系列IMF和殘余分量。EEMD分解后的各模態分量可看作是心電信號不同成分，噪聲以及R波的線性組合。即R波作為一個源信號分布在各個IMF中[19]。然后通過ICA算法從各模態分量中提取出R波。最后將ICA分離出的R波分量進行結構分析實現R波的識別。

3.1 EEMD分解參數的確定

EEMD對信號進行分解之前需要確定兩個十分重要的參數：添加的輔助白噪聲的大小k以及集合平均次數N。通常情況下，這兩個參數的設置是根據經驗進行確定，這大大降低了算法的自適應性和適用范圍。為了有效避免模態混疊，保護信號中的有用成分，使得信號的分解結果最優化，本文總結了k和N選取的具體過程。

(1)輸入待處理信號x(t)，并計算出該信號幅值標準差σx。

(2)將x(t)通過EMD分解成一系列IMF和殘余量，選取IMF1作為高頻分量，并計算出IMF1幅值標準差σx。

(3)根據EEMD中加入白噪聲的可依據準則[22]：0

3.2 本文算法應用流程

(1)通過對EEMD分解后的IMF做ICA分析，從IMF中提取出R波的源信號，通過對ICA提取的R波分量進行結構分析實現對R的識別。本文算法流程圖如圖1所示，具體實現步驟如下：輸入帶噪聲的心電信號，其中y(n)中混有工頻干擾、肌電干擾以及基線漂移3種噪聲，根據上文確定的添加白噪聲的大小k和集合平均的次數N對y(n)進行EEMD分解。

圖1 本文算法流程圖

(2)獨立成分分析是一種將多元(多維)統計數據分解為多個統計獨立且非高斯的分量線性和的統計方法。ICA是一個不斷迭代優化的過程，使得從觀測信號分離出的各個獨立分量最大程度接近各信號源。假設存在N個隨機信號z1,z2,···,z n，可由n個相互獨立的非高斯信號s1,s2,···,s n線性表示，即：z i=b i1s1+b i2s2+···+b in s n。其矩陣表達式為

式中，B為混合矩陣，S為獨立源向量，Z為觀測信號向量。

ICA分析就是在源信號S中各分量未知且相互獨立和B未知的條件下，由觀測信號Z估計分離矩陣W，用輸出信號Y估計S，即

為了從觀測信號X中分離出各個獨立源信號，本文采用基于最大負熵的FastICA算法。該方法以負熵為搜索目標，依次提取各獨立源信號。該算法詳細推導過程見文獻[22]。利用EEMD分解后的IMF構造ICA算法中的觀測信號矩陣Z

并利用基于最大負熵的FastICA算法從觀測信號矩陣中提取出R波的源信號H(n)。

(a)滑動窗口內的特征值滿足正態分布。

(b)R峰是滑動窗口內的局部極大值。

(c)在滑動窗口內隨機選擇的特征為R峰的概率小于或等于32%。

(6)若特征H(i)滿足以上條件并且超出1σ界限，則該特征是一個離群值。利用該特征構造R峰值包絡節點

4 仿真結果及分析

為了驗證本文算法的R波識別效果，本文所使用的心電信號數據來自QT數據庫。該數據庫中的數據主要從現有的心電數據庫中選擇，包括美國麻省理工大學與Beth Israel醫院合作建立的MIT-BIH心律失常數據庫，歐洲心臟病學會ST-T數據庫以及Beth Israel醫院收集的其他幾個心電數據庫[24]。該數據庫一共包含105個時長為15 min雙通道的動態心電記錄，采樣頻率為250 Hz，分辨率為11 bit，并且每個心電記錄均含有一種或多種噪聲干擾，同時伴隨著一種或多種心律失常。在每個記錄中，節拍都是由專家在一個小間隔內使用交互式圖形顯示手動確定。該數據庫的心電信號具有現實變化的各種QRS形態[25]。本文隨機選取12個心電記錄，并采用每個記錄中的第1通道的數據x(n)作為實驗心電信號源。在x(n)上疊加5 d B高斯白噪聲構造帶噪聲的心電信號y(n)( 其中y(n)混有工頻干擾、肌電干擾以及基線漂移3種噪聲，工頻干擾和基線漂移來自x(n)，肌電干擾采用5 d B的高斯白噪聲模擬)。為了驗證本文算法的實際應用效果，分別通過Pan-Tomkins算法，EEMD-ICA算法和本文算法分別對y(n)中的R波進行識別，利用仿真圖和評價指標對R波的識別效果進行定性分析和定量分析來評估本文算法的性能。

為了對R波的識別效果進行定量分析，引入靈敏度Sen，陽性準確率+P和準確率Acc這3個指標對仿真結果進行評估[23]。

(1)靈敏度

式中，TP為正確識別到的R波的個數，FN漏檢的R波的個數。

(2)陽性準確率

式中，TP為正確識別到的R波的個數，FP誤檢的R波的個數。

(3)準確率：

4.1 算法有效性驗證

為了驗證本文算法的有效性，通過Pan-Tomkins算法，EEMD-ICA算法和本文算法分別對心電信號中的R波進行識別，利用仿真圖和評價指標對R波的識別效果進行定性分析和定量分析來評估本文算法的性能。

(1)定性分析。圖2展示了信號sel223即x(n)和加入5 d B高斯白噪聲的sel223信號即y(n)。從圖2中可以看出，sel223信號中自帶工頻干擾和基線漂移，本文采用5 d B高斯白噪聲模擬肌電干擾并疊加在x(n)上構造信噪混合體。分別采用Pan-Tomkins算法，EEMD-ICA算法和本文算法對y(n)中的R波進行識別，為了方便更清楚地從仿真圖中看到R波的檢測情況，本文將在y(n)中識別到的R波的位置標記在來自Q T 數據庫的信號x(n)上。在利用EEMD-ICA算法和本文算法對y(n)中的R波進行識別之前，需要對EEMD中添加的輔助白噪聲的大小k以及集合平均次數N進行確定，由于參數的確定是根據待處理的信號自身特征自適應確定的，不同信號的k和N的值是不一樣的。

圖2 原始sel223信號和添加5 dB高斯白噪聲的sel223信號

從圖3中可以看出，Pan-Tomkins算法對混有工頻干擾，基線漂移以及肌電干擾的心電信號y(n)進行R波檢測的時候，在第85個采樣點處，由于該處R波的幅值與附近采樣點的R波幅值相比較小從而發生了漏檢現象；在第8571和9650個采樣點處，由于在8166～10600之間心電信號發生失常，導致Pan-Tomkins算法在這兩個位置發生誤檢現象。從圖4中可以看出，EEMD-ICA算法在對y(n)進行R波檢測的時候，在第85個采樣點處發生漏檢現象；在第8571個采樣點處發生誤檢現象。EEMD-ICA算法通過ICA算法對EEMD分解的各模態分量進行盲源分離，從各模態分量中分離出R波的源信號，起到了增強R波特征的作用，但是該方法對R波的檢測方法仍采取了Pan-Tomkins算法，雖然對誤檢現象有所改善，但改善效果并不理想。從圖5中可以看出，本文算法正確檢測出y(n)中所有的R波，消除了前兩種算法在在第85個采樣點處發生的漏檢現象，修正了前兩種算法在在第8571個采樣點處發生的誤檢。

圖3 Pan-Tomkins算法對y (n)的R波檢測結果

圖4 EEMD-ICA算法對y (n)的R波檢測結果

圖5 本文算法對y (n)的R波檢測結果

(2)定量分析。本文從QT數據庫隨機選取12個心電記錄，并采用每個記錄中的第1通道的數據x(n)作 為實驗心電信號源。在x(n)上疊加5 d B高斯白噪聲構造帶噪聲的心電信號。針對帶噪聲心電信號y(n)中R波的識別，本文通過靈敏度，陽性準確率，準確率這3個指標來定量評估本文算法的性能，結果如表1所示。

從表1中可以看出，本文提出的算法對疊加5 dB高斯白噪聲的12組心電信號進行R波的識別，在15411個R波中漏檢10個R波，誤檢10個R波，錯檢20個，具有很好的識別效果；此外，識別的靈敏度可達99.94%，陽性準確率可達99.94%，準確率達到99.87%。由圖3和圖4可知，EEMD-ICA算法對帶高斯白噪聲的sel223信號片段的識別效果略優于Pan-Tommkins算法，但隨著測試數據的增多，通過評價指標可得到EEMD-ICA算法性能不如Pan-Tomkins算法。

表1 本文算法R波識別性能評估

本文分別通過Pan-Tomkins算法、EEMD-ICA算法以及帶預處理的識別算法[26]對帶噪聲的心電信號y(n)進行R波識別，識別性能統計結果如表2所示。與帶預處理算法比較，由于存在噪聲干擾，該算法需要首先進行預處理操作，通過EMD閾值法對噪聲進行濾除后再進行識別。從表2可以看出，本文算法識別性能優于該帶預處理算法，并且耗時遠遠小于該算法，表現出更優異的識別性能。

表2 3種R波識別算法性能對比

Pan-Tomkins算法中的導數濾波導致心電信號丟失前后的相關信息，從而使得漏檢現象十分明顯。與Pan-Tomkins算法進行比較，本文算法充分利用了信號中的有效信息，不存在前面算法的相關信息丟失現象，在耗時方面雖弱于Pan-Tomkins算法，但對帶噪心電信號R波識別的漏檢數和誤檢數分別減少了85個和15個，漏檢現象得到明顯改善，具有更好的識別性能；并且在識別靈敏度、陽性準確率以及準確率方面均有較大的提升。由于實際應用中R波準確定位對疾病的診斷尤為重要，盡管本文算法耗時相對較長，但具有更強的適用性，更能夠體現算法的實際應用效果。與EEMD-ICA算法進行比較，本文算法在識別靈敏度和準確率方面有著較為明顯的提升，對帶噪聲心電信號R波識別的漏檢數和誤檢數分別減少了134個和30個，識別性能得到顯著提升并且耗時略較少。

4.2 算法穩定性分析

為了驗證本文算法的穩定性，利用其去識別病變明顯的心電信號中的R波。本文分別對最容易出現錯檢的長停頓的心電信號片段和T波高大的心電信號片段進行R波的識別。R波的識別結果分別見圖6和圖7。

從圖6和圖7中可以看出，本文算法針對最容易出現錯檢的長停頓的心電信號片段和T波高大的心電信號片段的R波具有良好的識別效果。

圖6 長停頓心電信號片段R波識別效果

圖7 T波高大心電信號片段R波識別效果

5 結束語

本文針對帶噪聲的心電信號提出一種無需預處理過程直接識別R波的算法。首先利用EEMD將帶噪聲的心電信號分解成一系列IMF和殘余分量。EEMD分解后的各模態分量可看作是心電信號不同成分，噪聲以及R波的線性組合。即R波作為一個源信號分布在各個IMF中[15]。然后通過ICA算法從各模態分量中提取出R波。最后將ICA分離出的R波分量通過滑動窗口獲取R波的包絡線獲得初步檢測的R波，并對初步檢測到的R波根據優化標準進行不斷迭代優化，最終實現R的識別。利用ICA對存在于IMF中的目標特征進行增強的思想，可以推廣到適合多通道數據處理的MEMD[27]和NA-MEMD[28]之中。本文選用QT數據庫的12組心電記錄，并對其疊加5 dB的高斯白噪聲構造帶噪聲的心電信號作為本文算法的待識別信號，通過定量分析和定性分析對本文算法的R波識別效果進行衡量。此外，本文針對最容易出現錯檢的長停頓的心電信號片段和T波高大的心電信號片段進行識別來驗證本文算法的穩定性，仿真結果表明，本文算法對帶噪聲的心電信號以及異常心電信號都具有較高的識別性能，對心電信號R波的實時檢測具有十分重要的意義。