基于希爾伯特變換和自適應閾值的R波檢測算法

郭田雨，嚴榮國,2，方旭晨,徐玉玲，陶爭屹

(1.上海理工大學健康工程與科學學院，上海 200093； 2.上海介入醫(yī)療器械工程技術研究中心，上海 200093; 3.上海楊浦區(qū)市東醫(yī)院，上海 200090)

0 引 言

我國城鄉(xiāng)居民因心血管疾病死亡的人數大約占全球的三分之一，這個比例還在不斷上升[1]，而心電圖(electrocardiogram, ECG)是通過人體特定位置的電極測量大量心肌細胞產生的電位信號，涵括了心臟的許多生理和病理信息，其檢測和分析有助于心血管疾病的輔助治療。R波的定位是對心率監(jiān)測[2]、心律失常檢測[3]、心跳分類[4]等心電圖相關分析應用中的一個重要步驟，所以快速、準確地識別出信號中具有較大幅值和斜率的R波在心電圖分析中具有重大意義和價值。檢測R波的傳統(tǒng)算法[5]主要包括差分閾值算法[6]、小波變換算法[7]、經驗模態(tài)分解法[8]以及人工神經網絡類法[9]。這些方法在特定檢測環(huán)境下具有較高檢測率，但在應用到臨床或者處理強噪聲和異常心電信號時有一定的局限性。近幾年，更多研究者逐漸將深度神經網絡算法運用到R波檢測[10]中，然而，模型訓練需要大量的樣本，樣本的質量和數量會影響模型的泛化能力，很難滿足實際臨床應用。

希爾伯特變換最早是由Bolton和Westphal應用在心電圖分析中[11]，至今希爾伯特變換已與多種方法結合用來識別心電圖中的特征信號。本文提出一種基于希爾伯特變換和自適應閾值的R波檢測算法,可盡可能同時保證實時性和準確性，該方法可有效濾除心電噪聲并快速準確定位R波，實現容易，并且計算不復雜，具有較好的魯棒性和實時性。本文詳細介紹該算法流程，并使用PhysioNet數據庫中的MIT-BIH心率失常標準數據庫、QT健康數據庫、NST噪音壓力測試數據庫、European ST-T數據庫以及臨床采集的心電數據對該算法的性能進行驗證。分析和測試結果表明，該算法能夠在保證檢測準確度的前提下快速完成R波檢測任務。

1 希爾伯特變換原理

對一個連續(xù)時間信號x(t)經過Hilbert變換后如式(1)所示：

(1)

原信號的解析信號：

(2)

記原信號x(t)的頻譜為Sin(ω)，則解析信號y(t)對應的頻譜Yin(ω)為：

(3)

從解析信號中可以得到原信號x(t)的包絡B(t)(瞬時振幅)為：

(4)

B(t)在復平面上的瞬時相位為：

(5)

Hilbert變換根據其物理意義可知，它是一個奇函數[12]，即Hilbert變換后的極值點對應于原始信號的過零點[13]。因此，可以用一階差分和Hilbert變換來提取預處理后的信號包絡定位R波。如圖1所示，選取200個采樣點，分別繪制原始信號、差分信號(經過濾波之后)、希爾伯特變換后信號，可以看到在標注線x=17采樣點的位置，與3條曲線交點的Y值分別對應原始信號的R峰、差分信號的過零點、希爾伯特變換包絡的極值點。

圖1 Hilbert變換和原始心電信號的對應關系

2 算法實現

本文提出的R波檢測算法首先對原始心電信號進行濾波、一階差分、幅值歸一化等預處理，P波和T波以及基線漂移等噪聲被抑制，同時突出R波的能量，處理后的信號數據有利于Hilbert包絡提??；然后，使用Hilbert變換提取包絡，對包絡信號進行平滑、增強其幅值處理，提高下一步R波定位的準確性；最后，結合自適應雙閾值方法對R波進行定位和檢測。本文的算法流程如圖2所示。

圖2 R波檢測算法流程

2.1 心電信號預處理

人體內外環(huán)境很容易影響微弱的心電信號的采集，采集得到的心電信號往往伴隨著工頻干擾(50 Hz)、基線偏移(0.05 Hz～1 Hz)、肌電干擾(60 Hz～300 Hz)和運動偽跡(3 Hz～14 Hz)等噪聲[13]，會影響后續(xù)R波檢測效果，所以需要經過一定的預處理來濾除噪聲。因為QRS復合波的能量主要集中在5 Hz～20 Hz，所以本文采用通頻帶為4 Hz～22 Hz的40階FIR帶通濾波器。本文選取MIT-BIH心律失常數據庫116號記錄進行處理(116號記錄包含嚴重漂移噪聲和小幅值R波)，濾波后的信號如圖3(b)所示，可以看到P波、T波被明顯削弱，R峰突出。

R波是心電信號最陡的部分，一階差分運算可以增強其斜率信息，使R波更加突出。差分后的極大值對的過零點對應原來信號的R峰點。一階差分的計算公式如式(6)所示：

d(n)=filter_ECG(n+1)-filter_ECG(n)

(6)

其中，filter_ECG(n)指濾波后的心電信號數據；d(n)為一階差分運算后的數據。

采用幅值歸一化處理差分后的數據，將幅值范圍擴展到-1～1,有利于下一步包絡提取，幅值歸一化公式如式(7)所示：

(7)

預處理前后過程的心電信號波形圖如圖3所示。由圖3可以看出，噪聲、P波和T波經過預處理后得到了有效抑制。

2.2 希爾伯特包絡分析

經過差分和幅值歸一化后，此時信號的過零點對應R峰，再對其做Hilbert變換求包絡，包絡信號的極大值點即為R峰點，在原始心電信號經過預處理后進行Hilbert變換求包絡處理，可以增強中、高等強度的信號能量。

變換結果如圖4所示。在圖4(a)中能夠看出，通過希爾伯特變換后，會得到小而排列緊密的波峰，需要提取其包絡線，使能量更加集中。由圖4(b)可以看出經過Hilbert變換包絡過程得到了較為平滑的包絡。在處理含噪比較嚴重或者P波、T波對R波的影響較大的心電信號時，還需要窗口大小為36個采樣點(當處理數據庫采樣頻率為250 Hz時，設置窗口大小為25個采樣點)的移動均值濾波處理使包絡更加平滑，因為MIT-BIH心律失常數據庫中的心電數據是以360個/s采樣點進行采樣，即頻率為360 Hz，而成年人的QRS波群的時間一般為0.06 s～0.11 s[14]，濾波后的波形如圖4(c)所示。

(a) 原始心電信號

(a) Hilbert變換結果

2.3 R波增強

經過上述步驟后，原始心電信號的幅值越來越小，不利于R波定位和檢測。本文采用滑動窗口積分，大小為17個采樣點來增大絕對振幅，進一步平滑波形。滑動窗口積分公式為：

z(n)=[y(n)+y(n+1)+…+y(n+17)]

(8)

其中，y(n)是移動均值濾波后數據，z(n)是滑動窗口積分后數據。

滑動窗口積分處理波形圖如圖5所示。

圖5 滑動窗口積分

2.4 R波定位

閾值的設定最為關鍵，也是決定檢測算法是否準確的因素之一。ECG信號的形狀和幅值會由于不同受試者的生理差異以及同一受試者不同時間的影響而變化，因此，采用自適應方法來對閾值進行設定。本文采用雙閾值，通過比較波峰大小與上、下限閾值來更新閾值[15]，檢測步驟如下：

Step1從第8個數據點開始，利用極值點判斷當前點與其前后共15個數據點的峰值大小，當y(i)>y(i-j)并且y(i)>y(i+j)(1≤j≤7)成立則篩選出候選波峰。

Step2設定初始雙閾值，選取濾波后數據前3 s內最大峰值的4倍作為高閾值Thr1，選取窗口積分后數據前3 s內平均值的2.5倍作為低閾值Thr0，依次比較候選波峰與初始雙閾值的大小。

Step2.1若候選波峰大于高閾值，則認為檢測到一個R波，直接存儲在R_buf1中，但說明閾值偏低程度很大，更新高低閾值為:

Thr1=0.7×mean(AMP_buf1)

(9)

Thr0=0.25×mean(AMP_buf1)

(10)

Step2.2若候選波峰在高低閾值之間，則認為檢測到一個R波，直接存儲在R_buf1中，閾值更新如下：

(11)

Thr0=0.4×peak

(12)

其中，peak指當前檢測到的峰值，AMP_buf1指存儲當前峰值前后共8個特征峰值。

Step2.3若候選波峰小于低閾值，則認為沒有檢測到R波。

Step3為了避免閾值變化較大，設定了閾值的下限，并通過多次實驗驗證，高閾值下限Thr1_lim=4.5，低閾值下限Thr0_lim=0.35，若檢測的當前R波峰值peak小于2個閾值下限時，自適應閾值分別為2個下限值。

Step4R波回檢和漏檢。

(13)

其中，R_buf(end)為檢測到的最后一個R峰值，mean_Interval為最后檢測到的8個RR間期的平均間期。從第2個R峰開始依次檢測，如果RR間期小于0.4平均間期，則認為R波可能存在多檢的情況，將當前R波的峰值與前、后R波的峰值的絕對值進行比較，保留較大的值作為正確的R波，如果RR間期大于平均間期的1.66倍，則認為可能存在漏檢的R波，那么在當前R峰右邊0.06 s到后一個R峰左邊0.06 s之間，尋找幅值比前后8個平均幅值大的點。

上述自適應閾值方法定位的R波檢測結果如圖6所示，圖6(a)是R波在特征信號中的定位結果，圖6(b)是R波在原始信號中的定位結果。

(a) R波在特征信號中定位結果

3 實驗測試結果和分析

3.1 評估指標

本文以敏感度(Se)、真陽性率(+P)、準確率(Acc)和檢測誤差(DER)作為評價算法性能的指標，計算公式如下：

(14)

(15)

(16)

(17)

其中，TP(True Positive)表示真實R波被正確檢出的個數；FP(False Positive)表示非R波被判定為R波的個數；FN(False Negative)表示真實R波被漏檢的個數；TB代表每條心電記錄標記的R波總個數。DER值越小，檢測效果越好。

3.2 心電信號的檢測結果

本文使用了4個具有不同采樣頻率和信噪比的數據庫，包括MITDB(Fs=360 Hz)，QTDB(Fs=250 Hz)，NSTDB(Fs=360 Hz)，EDB(Fs=250 Hz)。MIT-BIH心律失常數據庫驗證結果如表1所示，分別統(tǒng)計了共48條心電數據的正確檢出、誤檢和漏檢數量，最終得到總體準確率為99.13%的檢測結果。

表1 MITDB的R波檢測結果

從表1可以看出，對于大多數心電記錄敏感度、真陽性率、準確率檢測效果較好，本文方法在這幾個編號(104，105，108，116，201，207)的記錄檢測效果不佳。MIT-BIH數據庫中包括在其他樣本數據庫中不太常見但臨床上明顯的心律失常心電信號，表2是該數據庫中具有代表性的8條異常心電信號。本文將這8條記錄與其他文獻中所提方法檢測結果進行對比。這些記錄在常用方法中檢測效果也不理想。對比結果如表3所示。

表2 典型異常心電信號

表3 不同算法在MIT-BIH數據庫性能評估

表4是本文算法和其他文獻算法在QTDB、NSTDB、EDB這3個數據庫中對R波檢測結果對比，表3、表4中結果證明本文算法的敏感度Se、真陽性率+P較高，檢測誤差DER較低。這說明結合Hilbert和自適應雙閾值的方法可以有效提高算法的魯棒性和檢測準確性。

表4 不同算法在QTDB、NSTDB、EDB中性能評估

表5是使用本文算法對來自臨床的一位70歲女性患有房顫和I型傳導阻滯病的病人10組心電數據檢測結果。從表5可以得到10組臨床數據檢測平均準確率為97.51%，在測量過程中患者多次說話和移動身體，所以噪聲非常嚴重。房顫患者的心電信號特征表現為P波消失和RR間期絕對不規(guī)則[24]，I型傳導阻滯的心電信號特征表現為PR間期延長至大于200 ms，故2種異常心電疊加的信號檢測準確率較低[25]。

表5 臨床數據檢測結果

為了評估本文算法的計算復雜度以及實時性，本文在同一實驗環(huán)境下，即CPU為AMD Ryzen 5 3500U with Radeon Vega Mobile Gfx@2.10 GHz，內存為8 GB，操作系統(tǒng)為64位Windows10，MATLAB版本為2019a，將本文算法與Pan and Tompkins算法[26]進行對比實驗。Pan and Tompkins算法處理一條時長約30 min的心電數據平均時間為17.02 s，而本文算法約為2.85 s，大大縮短了計算所耗費的時間。圖7描述了MITDB中每條記錄的時間消耗。

圖7 處理一條記錄的時間(MITDB)

4 結束語

本文提出了一種基于希爾伯特變換和自適應雙閾值的R波檢測算法。該方法首先通過預處理消除噪聲，突出R波能量，之后通過Hilbert進一步抑制了P波和T波，并對變換后的信號應用自適應雙閾值對其R波快速準確定位。本文對MITDB、QTDB、NSTDB、EDB這4個數據庫中的心電數據進行處理，并將結果與其他的方法進行了比較，結果表明本文算法具有較高的敏感度、真陽性率和較低的檢測誤差，并且計算時間比Pan and Tompkins算法短，表明本文算法具有良好的魯棒性和實時性。