一種新的心電信號R峰自動檢測方法

孫亞楠,呂可嘉,張 瑞

(1.西北大學 醫學大數據研究中心， 陜西 西安 710127；2.西安交通大學 醫學部， 陜西 西安 710061)

心電圖(Electrocardiogram,ECG)是通過心電圖儀從體表記錄的能夠反映心臟興奮產生、傳導和恢復過程中所產生的生物電位變化情況的曲線[1]。一次正常的心跳(即一個心動周期),通常由P波、QRS波群和T波組成,其持續時間大約為0.8s(如圖1所示)。最早出現的波形是由心房興奮產生的P波(即反映心房除極過程),通常表現為向上的圓鈍平滑小波,正常的持續時間為0.06s～0.12s;其次是由心室興奮所產生的QRS波群,為尖銳的上下波(即由向下的Q波,向上的R波,及緊跟的S波共同構成),表示心室除極過程和心房復極過程的復合,正常的持續時間為0.06s～0.1s;最后出現的是由心室興奮恢復而產生的緩慢T波,代表心室復極過程,正常持續時間為0.1s～0.25s[2]。

R波具有較高的振幅和較大的斜率,是ECG信號中最重要的組成部分,蘊含著重要的生理病理信息。R波波峰(R峰)檢測是心率變異性分析、心血管疾病診斷、生物識別和ECG編碼系統的首要步驟。只有準確檢測到R峰后,才能正確定位出P波和T波,進一步計算心率、RR間期、PR間期、ST 段偏移等心電圖特征,從而完成心電圖分析等工作。因此,為了實現海量心電數據的實時判讀,開展心電信號R峰自動檢測的研究具有極為重要的臨床意義。

圖1 一個正常的心動周期Fig.1 Anormal cardiac cycle

20世紀50年代起,R峰的自動檢測技術逐漸興起,眾多R峰自動檢測方法被提出。這些方法大多基于差分運算[3]、數字濾波器[4-7]、經驗模式分解(EMD)[8]、幾何模板匹配[9]和神經網絡[10]等。基于差分運算的檢測算法,主要利用R峰具有較大幅值和斜率的特點,對去噪后的ECG信號進行差分運算并與設定的閾值進行比較,從而確定R峰[13]。這類算法容易理解、操作簡單、運算速度快,但其抗噪性能較差從而導致R峰檢測的準確率較低。基于幾何模板匹配的檢測算法,主要利用相鄰心動周期內心電信號具有較低差異性的特點,對任意心動周期內的心電信號與QRS波群模板的相關性進行度量, 進而依據相關性大小篩選出待檢QRS波群, 最終將振幅極大值確定為R峰[11]。 該算法的局限性在于計算復雜、 心電模板難以建立。 基于數字濾波器的檢測算法主要采用具有固定帶寬的濾波器將非QRS波群頻率占優帶內心電信號成分濾除,僅保留QRS波群頻率占優帶內的有效成分并結合決策規則完成R峰檢測。 該算法運行速度快, 但其檢測性能與濾波器帶寬的選擇和集成器移動窗口的大小密切相關, 采用不同帶寬的數字濾波器對同一心電信號進行R峰檢測, 檢測結果有顯著差異。 基于經驗模式分解的檢測算法, 能夠將心電信號自適應地分解成若干個本征模態函數, 根據QRS波群的頻率分布范圍選擇幾個恰當的本征模態函數進行心電信號的重構, 進而在重構信號上進行模極值對的檢測并將正模極值確定為R峰。 該算法不需要對心電信號進行預處理, 很好地避免了參數選擇問題, 但在含噪情況下選擇恰當的本征模態函數(IMF)進行ECG重構極其困難且EMD分解十分耗時, 從而導致該類算法在臨床中難以應用。

本文提出了一種新的基于小波變換的R峰自動檢測方法。首先采用基于小波包分解的去基線漂移方法以及改進的小波閾值法對原始心電信號進行去噪處理;繼而,對去噪后的心電信號進行小波變換,在恰當的頻率子帶上檢測出候選R峰,并將其對應位置投影于心電信號;最終,根據候選R峰在心電信號上的位置信息進行精確篩選,完成R峰檢測。為了進一步說明所提R峰自動檢測方法的有效性,本文選用MIT-BIH心律失常數據庫進行驗證并采用錯誤檢測率(DER)、真陽性率(+P)、敏感度(Se)、準確率(Acc)4個指標來評估該算法的性能。

1 小波變換與小波包分解

1.1 小波變換

小波變換是一種進行信號分析的時間-尺度方法,它具有多尺度分辨率分析的特點,在時域和頻域上同時具有表征信號局部特征的能力。通過伸縮和平移變換來實現心電信號在多個不同尺度上的細化,將心電信號的時域及頻域的局部特征同時反映在各個不同的頻率帶上,突出心電信號的不同特征。小波變換可分為連續小波變換(Continuous wavelet transform, CWT)和離散小波變換(Discrete wavelet transform, DWT)[11-13]。

由于連續小波變換中的基函數具有較強的相關性,會導致變換后得到的小波系數產生信息冗余,從而影響信號分析的質量,故本文僅采用離散小波變換對心電信號進行分析。給定母小波ψ(t),采用尺度因子aj=2j和位移因子bj,k=2jk對ψ(t)進行伸縮和平移變換,得到一組小波基函數序列

在此基礎上，信號S(t)的離散小波變換(DWT)定義為

給定信號S,其離散小波變換過程如圖2所示。設采樣率為2L,則由奈奎斯特采樣定律可知,第一層分解后可得L/2～L高頻子帶信號D1(細節部分)及0～L/2低頻子帶信號A1(近似部分),其中D1由信號S經過高通濾波器濾波和二元下采樣得到,A1由信號S經過低通濾波器濾波和二元下采樣得到;接下來,對A1進行第二層分解,可得L/22～L/2高頻子帶信號D2和0～L/22低頻子帶信號A2;依次類推。

圖2 離散小波變換的示意圖Fig.2 The diagram of discrete wavelet transformation

1.2 小波包分解

小波包分解是一種能夠很好地對信號進行時頻局部化分析的信號處理方法。相比小波變換而言,它能夠同時對信號的低頻以及高頻成分進行表征,為信號中未被細化的高頻部分提供更精細的分解,從而發掘高頻部分蘊含的重要信息,更加充分地對信號進行分析和研究。

給定一個信號S={s(1),s(2),…,s(n)}。令{H,G}為一對正交鏡像濾波器,其中H={hk}k∈Z為低通濾波器且G={gk}k∈Z為高通濾波器,則S的小波包分解系數可表示為[11]

(1)

其中

gk=(-1)kh1-k。

圖3 小波包分解示意圖Fig.3 The diagram of wavelet packet decomposition

2 一種新的R峰自動檢測方法

2.1 心電信號的去噪方法

ECG信號在采集過程中會受到多種噪聲的干擾,主要包括基線漂移、工頻干擾和肌電偽跡[2]。這些背景噪聲的存在會大大降低R峰檢測的可靠性。因此,為了高效地實現R峰的自動檢測,需要對ECG信號進行去噪處理。本小節將介紹一種新的基于小波包分解與改進的小波閾值法相結合的心電信號去噪方法。

實際應用中最常用的去噪方法是由Donoho等人提出的軟閾值去噪法[14],其閾值函數為

(2)

(3)

根據所定義的閾值函數(3),本文提出一種新的心電信號去噪方法(見算法1)。

算法1(心電信號去噪方法) 給定原始心電信號S={s(1),s(2),…,s(n)},其中n表示信號長度。

(4)

(5)

(6)

2.2 R峰的自動檢測方法

本小節介紹一種新的基于小波變換的R峰自動檢測方法。首先,將去噪后的心電信號劃分為若干等長的心電信號片段;其次,對每個心電信號片段進行離散小波變換,并選擇恰當的頻率子帶信號;再次在該子帶信號上采用自適應閾值法確定出候選R峰;最后,根據候選R峰在心電信號上的位置信息進行進一步篩選,以完成最終的R峰檢測。具體步驟見算法2。

或

2)確定動態閾值THR。

分別定義集合TUR的最小幅度平均值與最大幅度平均值為Min與Max,其中Min為TUR中前m1個最小元素的平均值,Max為TUR中前m2個最大元素的平均值。令

THR=d*(Max-Min)

(7)

3)確定候選R峰位置。

在TUR中選出大于動態閾值THR的所有點,記為

并稱其為候選R峰位置序列。

步驟4在Si上完成R峰的自動檢測。

1)采用小波變換逆變換將r映射到Si上,得到Si上的候選R峰位置序列,記為

2)按照下列循環對Si上的候選R峰位置進行微調,仍記為Ri:

①對Ri(p)進行半徑為v的加窗處理;

②將Ri(p)更新為[Ri(p)-v,Ri(p)+v]上最大值所對應的位置。

End

3)確定Si上最終的R峰位置。

首先,計算Si上的RR間期序列

其次,分別計算每一RR間期的左、右近鄰平均值為:

(8)

(9)

最后,按如下規則確定最終的R峰位置:若RR(p)

將所有保留的R峰位置序列記為R={R1,R2,…,Rl*},其中l*為最終檢測到的R峰個數。

3 數值實驗結果及分析

3.1 心電數據集

本文所使用的心電數據來自美國麻省理工大學和Beth Israel 醫院合作建立的MIT-BIH心律失常數據庫[17]。該數據庫一共包含48個時長約為30分鐘的雙導聯心電記錄,采樣率為360Hz,分辨率為11bit,且每個心電記錄中均含有一種或多種心律失常心電。本文從中隨機選取20個心電記錄,并采用每個記錄中第一個導聯的心電信號來驗證所提R峰自動檢測方法的有效性。文中所有實驗均Matlar 8.5.0中運行。

3.2 結果與分析

在去噪過程中,采用db6小波基函數對心電信號進行M=10層小波包分解。圖4展示了“記錄116”與“記錄208”兩段原始心電信號及其去噪后的心電信號。從圖中可以看出,本文所提出的去噪方法有效去除了基線漂移及工頻干擾等噪聲。

圖4 心電信號的去噪效果Fig.4 Denoising performance for the ECG signal

在R峰檢測過程中,取Q=180,則分割后的每個心電片段長度為T=3 600。在DWT中,采用db4小波基函數并取N=4。由于心電信號的能量主要集中在0～4Hz內,故在低頻子帶信號A2上進行候選R峰的確定。根據QRS波群寬度的參考值為0.06s～0.12s,實驗中選取參數d=0.3,v=45。

表1 R峰檢測方法的性能評估Tab.1 Performanceof the proposed R peak detection method

本文采用敏感度Se、錯誤檢測率DER、真陽性率+P和準確率Acc作為評估算法性能的度量指標,分別表示如下[18]:

(10)

(11)

(12)

(13)

其中TP(true positive,TP)表示正確檢測真實R峰的個數;FP(false positive,FP)表示將非R峰判定為R峰的個數; FN(false negative,FN)表示漏檢真實R峰的個數。

表1為本文所提R峰自動檢測方法應用于20個心電信號的數值實驗結果。從表1可看出本文所提方法的真陽性率、敏感度和準確率的均值分別達到99.67%，99.95%，99.62%。除心電記錄104，111，116，208和232外,其余檢測準確率均在99.77%以上。

由于心電記錄104，111和208中含有嚴重的基線漂移和波形突變現象,且受噪聲干擾和人工干預的影響較大,部分未被去除的噪聲尖波在檢測過程中被誤檢為R峰;心電記錄232中存在長達6s的停頓間歇,檢測過程中容易將停頓間歇內的毛刺誤檢為R峰,從而造成誤檢個數增多;心電記錄116和208中含有眾多的小振幅QRS波群,同時208中存在大范圍室性早搏現象,這些現象均會使得漏檢個數增多。圖5～10分別展示了上述5條心電記錄片段(時長15s)的R峰自動檢測結果。圖中紅色星號表示檢測出的R峰,綠色點代表真實的R峰位置。

圖5 含有嚴重肌肉噪聲的ECG片段的R峰自動檢測結果(記錄104)Fig.5 R peak detection result for the ECG signal with serious muscle noise (Record 104)

圖6 含有嚴重基線漂移和突變的ECG片段的R峰自動檢測結果(記錄111)Fig.6 R peak detection result for the ECG signal with serious baseline wander and abrupt changes (Record 111)

圖7 含有QRS波群形態突變的ECG片段的R峰自動檢測結果(記錄116)Fig.7 R peak detection result for the ECG signal with sudden changes in QRS amplitudes(Record 116)

圖8 包含小QRS波群的ECG片段的R峰自動檢測結果(記錄116)Fig.8 R peak detection result for the ECG signal with small QRS complexes (Record 116)

圖9 具有室性早搏和小QRS波群的ECG片段的R峰自動檢測結果(記錄208)Fig.9 R peak detection result for the ECG signal with wide premature ventricular contractions (PVCs), and small QRS complexes (Record 208)

圖10 包含長達6s停頓的ECG片段的R峰自動檢測結果(記錄232)Fig.10 R peak detection result for the ECG signal include long pauses up to 6s in duration (Record 232)

4 結 論

本文提出了一種新的基于小波變換的R峰自動檢測方法。首先,在小波包分解和小波閾值去噪法的基礎上,本文結合基線漂移、工頻干擾和肌電偽跡的頻率占優特點,提出了一種新的心電信號去噪方法,并采用該方法對原始心電信號進行去噪處理。其次,對去噪后的心電信號進行小波變換,采用自適應閾值法在恰當的頻率子帶上檢測出后候選R峰,并經逆變換將其位置映射到原始心電信號。最后,根據RR間期的局部變化趨勢進行R峰的精細篩選,以完成最終的R峰檢測。本文選用MIT-BIH心律失常數據庫中的20個心電記錄對所提R峰自動檢測方法進行性能評估,采用錯誤檢測率(DER)、真陽性率(+P)、敏感度(Se)、準確率(Acc)作為評價指標。實驗結果表明,本文所提方法具有較高的檢測性能且計算速度快,能夠為臨床輔助診斷提供一定的指導依據。

[1] 李中健，李世鋒，申繼紅，等.心電圖學系列講座(三)——心電圖一般知識[J]. 中國全科醫學,2014,17(03):360-362.

[2] 前田如矢. 一學就會心電圖[M].中國:華夏出版社,2010.

[3] YEH Y, WANG W. QRS complexes detection for ECG signal：The difference operation method[J].Computer Methods and Programs in Biomedicine, 2008, 91(3): 245-254.

[4] PAN J, TOMPKINS W J. A real-time QRS detection algorithm[J]. IEEE Transactions on Biomedical Engineering, 1985, 32(3): 230-236.

[5] HAMILTON P S, TOMPKINS W J. Quantitative investigation of QRS detection rules using the MIT/BIH arrhythmia database[J]. IEEE Transactions on Biomedical Engineering, 1986, 33(12): 1157-1165.

[6] ARZENO N M, DENG Z, POON C, et al. Analysis of first-derivative based QRS detection algorithms[J]. IEEE Transactions on Biomedical Engineering, 2008, 55(2): 478-484.

[7] BENITEZ D, GAYDECKI P, ZAIDI A, et al. The use of the Hilbert transform in ECG signal analysis[J]. Computers in Biology and Medicine, 2001, 31(5): 399-406.

[8] HONGYAN X, MINSONG H. A new QRS detection algorithm based on empirical mode decomposition[C]∥International Conference on Bioinformatics and Biomedical Engineering, 2008: 693-696.

[9] SUAREZ K, SILVA J, BERTHOUMIEU Y, et al. ECG beat detection using a geometrical matching approach[J]. IEEE Transactions on Biomedical Engineering, 2007, 54(4): 641-650.

[10] ABIBULLAEV B, SEO H D. A new QRS detection method using wavelets and artificial neural networks[J]. Journal of Medical Systems, 2011, 35(4): 683-691.

[11] 楊福生. 小波變換的工程分析與應用[M].北京:科學出版社,1999,1-145.

[12] GOMES P, SOARES F, CORREIA J H, et al. ECG data-acquisition and classification system by using wavelet-domain hidden markov models[C]∥International Conference of the Ieee Engineering in Medicine and Biology Society, 2010: 4670-4673.

[13] ADDISON P S. Wavelet transforms and the ECG: a review[J]. Physiological Measurement, 2005, 26(5).

[14] DONOHO D L. De-noising by soft-thresholding[J]. IEEE Transactions on Information Theory, 1995, 41(3): 613-627.

[15] DONOHO D L, JOHNSTONE I M. Adapting to unknown smoothness via wavelet shrinkage[J]. Journal of the American Statistical Association, 1995, 90(432): 1200-1224.

[16] ZHANG Q, WANG L, GAO P, et al. An innovative wavelet threshold denoising method for environmental drift of fiber optic gyro[J]. Mathematical Problems in Engineering, 2016: 1-8.

[17] GOLDBERGER A L, AMARAL L A, GLASS L, et al. PhysioBank, PhysioToolkit, and PhysioNet components of a new research resource for complex Physiologic signals[J]. Circulation, 2000, 101(23): 215-220.

[18] MANIKANDAN M S, SOMAN K P. A novel method for detecting R-peaks in electrocardiogram (ECG) signal[J].Biomedical Signal Processing and Control, 2012, 7(2): 118-128.