基于希尔伯特变换和自适应阈值的R波检测算法

郭田雨，严荣国,2，方旭晨,徐玉玲，陶争屹

(1.上海理工大学健康工程与科学学院，上海 200093； 2.上海介入医疗器械工程技术研究中心，上海 200093; 3.上海杨浦区市东医院，上海 200090)

0 引 言

我国城乡居民因心血管疾病死亡的人数大约占全球的三分之一，这个比例还在不断上升[1]，而心电图(electrocardiogram, ECG)是通过人体特定位置的电极测量大量心肌细胞产生的电位信号，涵括了心脏的许多生理和病理信息，其检测和分析有助于心血管疾病的辅助治疗。R波的定位是对心率监测[2]、心律失常检测[3]、心跳分类[4]等心电图相关分析应用中的一个重要步骤，所以快速、准确地识别出信号中具有较大幅值和斜率的R波在心电图分析中具有重大意义和价值。检测R波的传统算法[5]主要包括差分阈值算法[6]、小波变换算法[7]、经验模态分解法[8]以及人工神经网络类法[9]。这些方法在特定检测环境下具有较高检测率，但在应用到临床或者处理强噪声和异常心电信号时有一定的局限性。近几年，更多研究者逐渐将深度神经网络算法运用到R波检测[10]中，然而，模型训练需要大量的样本，样本的质量和数量会影响模型的泛化能力，很难满足实际临床应用。

希尔伯特变换最早是由Bolton和Westphal应用在心电图分析中[11]，至今希尔伯特变换已与多种方法结合用来识别心电图中的特征信号。本文提出一种基于希尔伯特变换和自适应阈值的R波检测算法,可尽可能同时保证实时性和准确性，该方法可有效滤除心电噪声并快速准确定位R波，实现容易，并且计算不复杂，具有较好的鲁棒性和实时性。本文详细介绍该算法流程，并使用PhysioNet数据库中的MIT-BIH心率失常标准数据库、QT健康数据库、NST噪音压力测试数据库、European ST-T数据库以及临床采集的心电数据对该算法的性能进行验证。分析和测试结果表明，该算法能够在保证检测准确度的前提下快速完成R波检测任务。

1 希尔伯特变换原理

对一个连续时间信号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 心电信号预处理

人体内外环境很容易影响微弱的心电信号的采集，采集得到的心电信号往往伴随着工频干扰(50 Hz)、基线偏移(0.05 Hz～1 Hz)、肌电干扰(60 Hz～300 Hz)和运动伪迹(3 Hz～14 Hz)等噪声[13]，会影响后续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所示，分别统计了共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间期绝对不规则[24]，I型传导阻滞的心电信号特征表现为PR间期延长至大于200 ms，故2种异常心电叠加的信号检测准确率较低[25]。

表5 临床数据检测结果

为了评估本文算法的计算复杂度以及实时性，本文在同一实验环境下，即CPU为AMD Ryzen 5 3500U with Radeon Vega Mobile Gfx@2.10 GHz，内存为8 GB，操作系统为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算法短，表明本文算法具有良好的鲁棒性和实时性。