孫亞楠,呂可嘉,張 瑞

(1.西北大學(xué) 醫(yī)學(xué)大數(shù)據(jù)研究中心， 陜西 西安 710127；2.西安交通大學(xué) 醫(yī)學(xué)部， 陜西 西安 710061)

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

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

圖1 一個(gè)正常的心動(dòng)周期Fig.1 Anormal cardiac cycle

20世紀(jì)50年代起,R峰的自動(dòng)檢測(cè)技術(shù)逐漸興起,眾多R峰自動(dòng)檢測(cè)方法被提出。這些方法大多基于差分運(yùn)算[3]、數(shù)字濾波器[4-7]、經(jīng)驗(yàn)?zāi)Ｊ椒纸?EMD)[8]、幾何模板匹配[9]和神經(jīng)網(wǎng)絡(luò)[10]等?；诓罘诌\(yùn)算的檢測(cè)算法,主要利用R峰具有較大幅值和斜率的特點(diǎn),對(duì)去噪后的ECG信號(hào)進(jìn)行差分運(yùn)算并與設(shè)定的閾值進(jìn)行比較,從而確定R峰[13]。這類(lèi)算法容易理解、操作簡(jiǎn)單、運(yùn)算速度快,但其抗噪性能較差從而導(dǎo)致R峰檢測(cè)的準(zhǔn)確率較低?；趲缀文０迤ヅ涞臋z測(cè)算法,主要利用相鄰心動(dòng)周期內(nèi)心電信號(hào)具有較低差異性的特點(diǎn),對(duì)任意心動(dòng)周期內(nèi)的心電信號(hào)與QRS波群模板的相關(guān)性進(jìn)行度量, 進(jìn)而依據(jù)相關(guān)性大小篩選出待檢QRS波群, 最終將振幅極大值確定為R峰[11]。 該算法的局限性在于計(jì)算復(fù)雜、 心電模板難以建立。 基于數(shù)字濾波器的檢測(cè)算法主要采用具有固定帶寬的濾波器將非QRS波群頻率占優(yōu)帶內(nèi)心電信號(hào)成分濾除,僅保留QRS波群頻率占優(yōu)帶內(nèi)的有效成分并結(jié)合決策規(guī)則完成R峰檢測(cè)。 該算法運(yùn)行速度快, 但其檢測(cè)性能與濾波器帶寬的選擇和集成器移動(dòng)窗口的大小密切相關(guān), 采用不同帶寬的數(shù)字濾波器對(duì)同一心電信號(hào)進(jìn)行R峰檢測(cè), 檢測(cè)結(jié)果有顯著差異。 基于經(jīng)驗(yàn)?zāi)Ｊ椒纸獾臋z測(cè)算法, 能夠?qū)⑿碾娦盘?hào)自適應(yīng)地分解成若干個(gè)本征模態(tài)函數(shù), 根據(jù)QRS波群的頻率分布范圍選擇幾個(gè)恰當(dāng)?shù)谋菊髂B(tài)函數(shù)進(jìn)行心電信號(hào)的重構(gòu), 進(jìn)而在重構(gòu)信號(hào)上進(jìn)行模極值對(duì)的檢測(cè)并將正模極值確定為R峰。 該算法不需要對(duì)心電信號(hào)進(jìn)行預(yù)處理, 很好地避免了參數(shù)選擇問(wèn)題, 但在含噪情況下選擇恰當(dāng)?shù)谋菊髂B(tài)函數(shù)(IMF)進(jìn)行ECG重構(gòu)極其困難且EMD分解十分耗時(shí), 從而導(dǎo)致該類(lèi)算法在臨床中難以應(yīng)用。

本文提出了一種新的基于小波變換的R峰自動(dòng)檢測(cè)方法。首先采用基于小波包分解的去基線漂移方法以及改進(jìn)的小波閾值法對(duì)原始心電信號(hào)進(jìn)行去噪處理;繼而,對(duì)去噪后的心電信號(hào)進(jìn)行小波變換,在恰當(dāng)?shù)念l率子帶上檢測(cè)出候選R峰,并將其對(duì)應(yīng)位置投影于心電信號(hào);最終,根據(jù)候選R峰在心電信號(hào)上的位置信息進(jìn)行精確篩選,完成R峰檢測(cè)。為了進(jìn)一步說(shuō)明所提R峰自動(dòng)檢測(cè)方法的有效性,本文選用MIT-BIH心律失常數(shù)據(jù)庫(kù)進(jìn)行驗(yàn)證并采用錯(cuò)誤檢測(cè)率(DER)、真陽(yáng)性率(+P)、敏感度(Se)、準(zhǔn)確率(Acc)4個(gè)指標(biāo)來(lái)評(píng)估該算法的性能。

1 小波變換與小波包分解

1.1 小波變換

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

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

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

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

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

1.2 小波包分解

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

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

(1)

其中

gk=(-1)kh1-k。

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

2 一種新的R峰自動(dòng)檢測(cè)方法

2.1 心電信號(hào)的去噪方法

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

實(shí)際應(yīng)用中最常用的去噪方法是由Donoho等人提出的軟閾值去噪法[14],其閾值函數(shù)為

(2)

(3)

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

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

(4)

(5)

(6)

2.2 R峰的自動(dòng)檢測(cè)方法

本小節(jié)介紹一種新的基于小波變換的R峰自動(dòng)檢測(cè)方法。首先,將去噪后的心電信號(hào)劃分為若干等長(zhǎng)的心電信號(hào)片段;其次,對(duì)每個(gè)心電信號(hào)片段進(jìn)行離散小波變換,并選擇恰當(dāng)?shù)念l率子帶信號(hào);再次在該子帶信號(hào)上采用自適應(yīng)閾值法確定出候選R峰;最后,根據(jù)候選R峰在心電信號(hào)上的位置信息進(jìn)行進(jìn)一步篩選,以完成最終的R峰檢測(cè)。具體步驟見(jiàn)算法2。

或

2)確定動(dòng)態(tài)閾值THR。

分別定義集合TUR的最小幅度平均值與最大幅度平均值為Min與Max,其中Min為T(mén)UR中前m1個(gè)最小元素的平均值,Max為T(mén)UR中前m2個(gè)最大元素的平均值。令

THR=d*(Max-Min)

(7)

3)確定候選R峰位置。

在TUR中選出大于動(dòng)態(tài)閾值THR的所有點(diǎn),記為

并稱(chēng)其為候選R峰位置序列。

步驟4在Si上完成R峰的自動(dòng)檢測(cè)。

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

2)按照下列循環(huán)對(duì)Si上的候選R峰位置進(jìn)行微調(diào),仍記為Ri:

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

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

End

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

首先,計(jì)算Si上的RR間期序列

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

(8)

(9)

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

將所有保留的R峰位置序列記為R={R1,R2,…,Rl*},其中l(wèi)*為最終檢測(cè)到的R峰個(gè)數(shù)。

3 數(shù)值實(shí)驗(yàn)結(jié)果及分析

3.1 心電數(shù)據(jù)集

本文所使用的心電數(shù)據(jù)來(lái)自美國(guó)麻省理工大學(xué)和Beth Israel 醫(yī)院合作建立的MIT-BIH心律失常數(shù)據(jù)庫(kù)[17]。該數(shù)據(jù)庫(kù)一共包含48個(gè)時(shí)長(zhǎng)約為30分鐘的雙導(dǎo)聯(lián)心電記錄,采樣率為360Hz,分辨率為11bit,且每個(gè)心電記錄中均含有一種或多種心律失常心電。本文從中隨機(jī)選取20個(gè)心電記錄,并采用每個(gè)記錄中第一個(gè)導(dǎo)聯(lián)的心電信號(hào)來(lái)驗(yàn)證所提R峰自動(dòng)檢測(cè)方法的有效性。文中所有實(shí)驗(yàn)均Matlar 8.5.0中運(yùn)行。

3.2 結(jié)果與分析

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

圖4 心電信號(hào)的去噪效果Fig.4 Denoising performance for the ECG signal

在R峰檢測(cè)過(guò)程中,取Q=180,則分割后的每個(gè)心電片段長(zhǎng)度為T(mén)=3 600。在DWT中,采用db4小波基函數(shù)并取N=4。由于心電信號(hào)的能量主要集中在0～4Hz內(nèi),故在低頻子帶信號(hào)A2上進(jìn)行候選R峰的確定。根據(jù)QRS波群寬度的參考值為0.06s～0.12s,實(shí)驗(yàn)中選取參數(shù)d=0.3,v=45。

表1 R峰檢測(cè)方法的性能評(píng)估Tab.1 Performanceof the proposed R peak detection method

本文采用敏感度Se、錯(cuò)誤檢測(cè)率DER、真陽(yáng)性率+P和準(zhǔn)確率Acc作為評(píng)估算法性能的度量指標(biāo),分別表示如下[18]:

(10)

(11)

(12)

(13)

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

表1為本文所提R峰自動(dòng)檢測(cè)方法應(yīng)用于20個(gè)心電信號(hào)的數(shù)值實(shí)驗(yàn)結(jié)果。從表1可看出本文所提方法的真陽(yáng)性率、敏感度和準(zhǔn)確率的均值分別達(dá)到99.67%，99.95%，99.62%。除心電記錄104，111，116，208和232外,其余檢測(cè)準(zhǔn)確率均在99.77%以上。

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

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

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

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

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

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

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

4 結(jié) 論

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

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

[2] 前田如矢. 一學(xué)就會(huì)心電圖[M].中國(guó):華夏出版社,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] 楊福生. 小波變換的工程分析與應(yīng)用[M].北京:科學(xué)出版社,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.