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.   集合经验模态分解

为了解决经验模态分解(Empirical Mode Decomposition, EMD)存在的模态混叠问题，Wu等人^[20]在原始信号上加入若干次辅助白噪声，把每一次构造的信噪混合体进行EMD分解，最后对各模态分量取平均值获取逼近的真实模态。EEMD具体步骤为

(1) 在原始信号$x(t)$上分别加入N次不同大小的辅助白噪声构成不同的信噪混合体

其中，${X_i}(t)$($i{\rm{ = }}1,2, ··· ,N$)为第$i$次加入辅助白噪声后的信号，${k_i}$为第$i$次加入的白噪声与原始信号$x(t)$幅值标准差比值，${\sigma _x}$为原始信号标准差，$n(t)$为归一化白噪声。

(2) 对信噪混合体${X_i}(t)$($i{\rm{ = }}1,2, ··· ,N$)进行EMD分解，得到$m$个IMF分量${c_{ij}}(t)$($j{\rm{ = }}1,2, ··· ,M$)和一个残余量${r_i}(t)$，可以表示为

(3) 对每次EMD分解对应的IMF分量和残余量求平均值

式中，${c_j}(t)$表示EEMD分解后得到的第$j$个IMF分量，${r_m}(t)$为EEMD分解后的残余量。

(4)最终得到的分解结果为

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)$，并计算出该信号幅值标准差$\sigma x$。

(2)将$x(t)$通过EMD分解成一系列IMF和残余量，选取IMF₁作为高频分量，并计算出IMF₁幅值标准差$\sigma x$。

(3)根据EEMD中加入白噪声的可依据准则^[22]：$0 < k < {A / 2}$，确定$k$的大小，其中$A{\rm{ = }}{{{\sigma _1}} / {{\sigma _x}}}$。一般情况下，$k{\rm{ = }}{A / 4}$就能有效避免信号在分解时出现的模态混叠问题。

(4)根据EEMD中集合平均次数$N$与白噪声大小$k$的关系：$e = {k / {\sqrt N }}$，确定$N$的大小，其中$e$为预设相对误差，通常设为1%即可满足要求^[21]。

3.2   本文算法应用流程

(1)通过对EEMD分解后的IMF做ICA分析，从IMF中提取出R波的源信号，通过对ICA提取的R波分量进行结构分析实现对R的识别。本文算法流程图如图1所示，具体实现步骤如下：输入带噪声的心电信号，其中$y(n)$中混有工频干扰、肌电干扰以及基线漂移3种噪声，根据上文确定的添加白噪声的大小$k$和集合平均的次数$N$对$y(n)$进行EEMD分解。

图 1  本文算法流程图

下载: 全尺寸图片 幻灯片

(2)独立成分分析是一种将多元(多维)统计数据分解为多个统计独立且非高斯的分量线性和的统计方法。ICA是一个不断迭代优化的过程，使得从观测信号分离出的各个独立分量最大程度接近各信号源。假设存在N个随机信号${z_1},{z_2}, ··· ,{z_n}$，可由n个相互独立的非高斯信号${s_1},{s_2}, ··· ,{s_n}$线性表示，即：${z_i} = {b_{i1}}{s_1} + {b_{i2}}{s_2} + ··· + {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)$。

(3) 从$H(n)$中提取出所有波峰的特征点，波峰集的 定义为

式中，${\rm{sgn}} (\cdot)$表示符号函数，$v$表示满足波峰集的特征值。

(4) 计算以时间步长$i$为中心，大小为$2\alpha + 1$的滑动窗口内特征值的平均值${\alpha _i}$和标准差${{\rm{sd}} _i}$

(5) 若特征$H(i)$是滑动窗口内的一个离群值，则假定它是构成R峰值包络线的一个节点，判断该特征是否是滑动窗口内的离群值，离群值需满足以下条件^[23]：

(a) 滑动窗口内的特征值满足正态分布。

(b) R峰是滑动窗口内的局部极大值。

(c) 在滑动窗口内随机选择的特征为R峰的概率小于或等于32%。

(6) 若特征$H(i)$满足以上条件并且超出$1\sigma$界限，则该特征是一个离群值。利用该特征构造R峰值包络节点

(7) 对${\rm{NT}}$按${t_i}$进行升序排列，对包络节点之间的值进行线性插值，使得每个特征$H(i)$在构造的R峰包络线${\rm{ET}} (i)$中都有对应的值。为了达到这一要求，滑动窗口内必须至少有一次心跳。因此在执行步骤(4)的时候，取$\alpha {\rm{ = }}{{\left( {{F_{\rm{s}}} - 1} \right)} / 2}$，${F_{\rm{s}}}$为心率。

(8) 根据构建的R峰值包络线确定R波的位置，但存在漏检和误检现象。需要对检测到的R波进行修正，对算法进行优化。搜索$Ri$内连续太近的波峰对，去除振幅最小的一个以排除误检。搜索$Ri$内连续太远的波峰对，将它们之间具有最高值的特征处判定为R波来“添加”一个R波，这样可以排除漏检现象。根据心率对R峰值太过接近和太过分离进行判断，判断标准分别为

其中，${R_i}(k)$表示第$i$次迭代搜索过程中R峰$k$，${R_i}(j)$表示第$i$次迭代搜索过程中R峰$j$，${F_{\rm{s}}}$表示所检测心电信号的心率，${{\rm{HR}} _{\max }}$表示最大心率，本文取${{\rm{HR}} _{\max }}{\rm{ = }}220\;{\rm{bpm}}$，${{\rm{HR}} _{\min }}$表示最小心率，本文取${{\rm{HR}} _{\min }}{\rm{ = }}30\;{\rm{bpm}}$。

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 dB高斯白噪声构造带噪声的心电信号$y(n)$(其中$y(n)$混有工频干扰、肌电干扰以及基线漂移3种噪声，工频干扰和基线漂移来自$x(n)$，肌电干扰采用5 dB的高斯白噪声模拟)。为了验证本文算法的实际应用效果，分别通过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 dB高斯白噪声的sel223信号即$y(n)$。从图2中可以看出，sel223信号中自带工频干扰和基线漂移，本文采用5 dB高斯白噪声模拟肌电干扰并叠加在$x(n)$上构造信噪混合体。分别采用Pan-Tomkins算法，EEMD-ICA算法和本文算法对$y(n)$中的R波进行识别，为了方便更清楚地从仿真图中看到R波的检测情况，本文将在$y(n)$中识别到的R波的位置标记在来自QT数据库的信号$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 dB高斯白噪声构造带噪声的心电信号。针对带噪声心电信号$y(n)$中R波的识别，本文通过灵敏度，阳性准确率，准确率这3个指标来定量评估本文算法的性能，结果如表1所示。

表 1  本文算法R波识别性能评估

ECG记录	R峰总数	漏检	误检	错检总数	灵敏度(%)	阳性准确率(%)	准确率(%)
sel100	1134	0	3	3	100.00	99.74	99.74
sel103	1048	0	0	0	100.00	100.00	100.00
sel116	1185	0	0	0	100.00	100.00	100.00
sel213	1642	0	1	1	100.00	99.94	99.94
sel221	1247	1	4	5	99.92	99.68	99.60
sel223	1309	3	2	5	99.77	99.85	99.62
sel230	1077	0	0	0	100.00	100.00	100.00
sel301	1351	2	0	2	99.85	100.00	99.85
sel310	2012	3	0	3	99.85	100.00	99.85
sel803	1026	0	0	0	100.00	100.00	100.00
sel853	1113	1	0	1	99.91	100.00	99.91
sel891	1267	0	0	0	100.00	100.00	100.00
合计	15411	10	10	20	99.94	99.94	99.87

下载: 导出CSV 

| 显示表格

从表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算法。

本文分别通过Pan-Tomkins算法、EEMD-ICA算法以及带预处理的识别算法^[26]对带噪声的心电信号$y(n)$进行R波识别，识别性能统计结果如表2所示。与带预处理算法比较，由于存在噪声干扰，该算法需要首先进行预处理操作，通过EMD阈值法对噪声进行滤除后再进行识别。从表2可以看出，本文算法识别性能优于该带预处理算法，并且耗时远远小于该算法，表现出更优异的识别性能。

表 2  3种R波识别算法性能对比

R波识别算法	R波总数	漏检	误检	错检总数	灵敏度(%)	阳性准确率(%)	准确率(%)	平均处理时间(s)
Pan-Tomkins算法[4]	15411	95	25	120	99.38	99.84	99.22	1.7194
EEMD-ICA算法[17]	15411	144	40	184	99.07	99.74	98.81	76.9896
本文算法	15411	10	10	20	99.94	99.94	99.87	76.9335
带预处理算法[26]	15411	13	140	153	99.91	99.09	99.02	114.607

下载: 导出CSV 

| 显示表格

Pan-Tomkins算法中的导数滤波导致心电信号丢失前后的相关信息，从而使得漏检现象十分明显。与Pan-Tomkins算法进行比较，本文算法充分利用了信号中的有效信息，不存在前面算法的相关信息丢失现象，在耗时方面虽弱于Pan-Tomkins算法，但对带噪心电信号R波识别的漏检数和误检数分别减少了85个和15个，漏检现象得到明显改善，具有更好的识别性能；并且在识别灵敏度、阳性准确率以及准确率方面均有较大的提升。由于实际应用中R波准确定位对疾病的诊断尤为重要，尽管本文算法耗时相对较长，但具有更强的适用性，更能够体现算法的实际应用效果。与EEMD-ICA算法进行比较，本文算法在识别灵敏度和准确率方面有着较为明显的提升，对带噪声心电信号R波识别的漏检数和误检数分别减少了134个和30个，识别性能得到显著提升并且耗时略较少。

4.2   算法稳定性分析

为了验证本文算法的稳定性，利用其去识别病变明显的心电信号中的R波。本文分别对最容易出现错检的长停顿的心电信号片段和T波高大的心电信号片段进行R波的识别。R波的识别结果分别见图6和图7。

图 6  长停顿心电信号片段R波识别效果

下载: 全尺寸图片 幻灯片

图 7  T波高大心电信号片段R波识别效果

下载: 全尺寸图片 幻灯片

从图6和图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波的实时检测具有十分重要的意义。