基于VMD与包络峰值提取的R波检测算法

庞 宇，汪立宇，陈亚军

(重庆邮电大学 光电工程学院，重庆 400065)

0 引 言

心电(electrocardiogram，ECG)信号的波形是诊断心血管疾病的重要参考。传统上心电图的分析是由QRS波群开始的，它的能量高、振幅大，因此最容易被检测。而QRS波群的定位检测的第一步就是先检测R峰的位置，完成该处理后才能根据规律检测Q波和S波[1]。

对于ECG信号的R波检测，有差分阈值[2]、小波变换[3]等算法，均具有良好的检测定位效果。但是，传统的差分阈值法检测R波原理简单，对低频干扰抑制作用好，但容易受高频干扰影响从而降低识别率；小波变换检测计算量略大且对硬件平台要求高。此外还有香农能量变换法[4]、粒子群优化算法[5]、神经网络[6]的模式识别算法等，各种方法都在特定的检测环境下表现出不同的优势。经验模态分解(empirical mode decomposition，EMD)及相关理论的提出，为非平稳信号的处理找到了一条新途径，也适用于ECG的波形检测[7]，但算法本身存在模态混叠的缺点。Dragomiretskiy K等提出了变分模态分解算法(variational mode decomposition，VMD)[8]，该算法可有效消除EMD存在的模态混叠问题，在生物体征信号的分析处理领域具有很大的应用价值。

VMD是在EMD基础上，结合数学变分法的理论提出的一种信号分析方法，适用于非平稳信号的分析，在心电信号R波检测中具有较强的抗干扰能力与较高的识别率[9]。本文通过对预设尺度K值选取设定最优化算法，将ECG信号分解提取出特征化模态分量，然后对该分量进行包络峰值提取且设定复检规则，能有效排除背景噪声的影响，提高R波定位的准确率。

1 R波检测定位算法

1.1 VMD处理

VMD是一种处理非线性、非平稳信号的有效方法，可以将原始ECG信号分解成有限数量的二维本征模态函数(bi-dimensional intrinsic mode function，BIMF)[10]，当再现输入信号时，经分解得到的子信号具有特定的稀疏特性[11]。

VMD将EMD方法中本征模态函数定义为窄带BIMF，即

xk(t)=Ak(t)cos[φk(t)]

(1)

其中，φk(t)非递减，φ′k(t)≥0；包络线非负，即Ak(t)≥0；且Ak(t)与瞬时频率ωk(t)=φ′k(t)相对于相位φk(t)的变换是缓变的。

VMD方法将信号放入约束变分模型中求取最优解，确定限带本征模函数的中心频率及带宽，模型表达式为

(2)

(3)

其中，xk(t)、ωk分别表示信号的第k个模态分量及其对应的中心频率，{xk}={x1,x2,…,xK}、{ωk}={ω1,ω2,…,ωK}则分别表示K个模态分量及中心频率的集合。上述约束优化问题引用增广拉格朗日修正方程求解，表达式如下

(4)

式中：α——惩罚参数；λ——拉格朗日因子。

根据文献[10，11]中的方法，从谱域解得到的模态及对应的更新中心频率表示为

(5)

(6)

然后利用交替方向乘子算法[12]求该模型最优解，得到的结果即是信号分解后的K个限带模态分量。算法步骤如下：

其中τ为双上升时间步长，ε为收敛性判别标准的公差。

原始ECG信号包含有PQRST各个波段有效信息，以及工频、肌电、基线漂移、运动伪迹干扰等各种噪声，经过以上过程即可分解为从低频区到高频区的各个模态分量，对于R波的检测定位只需分析特定模态分量以排除噪声与其它无效信息的干扰。

1.2 K值选取

VMD技术中，K值需要预先设定，其大小直接影响了信号分解的结果。文献[13]指出，当K值最优时，BIMF末分量的中心频率将达到最大值。因此，可以通过观察中心频率随不同K值的变化，确定最优预设尺度。

这里设定确定K值的算法：初设一个较大k值，令K=k和K=k-1，若两次分解得到的BIMF末分量的中心频率相近，则说明设定的k值偏大；再取K=k-1和K=k-2，依次循环，当两次分解得到的最后一个BIMF分量中心频率差异较大时，可给出最优K值。通过设定一个阈值a(视信号频域而定)来判断中心频率是否相近。在此通过一个仿真信号进行分析。仿真信号如下

f=0.5cos(2*pi*2*t)+0.25cos(2*pi*24*t)+
0.2cos(2*pi*288*t)+η

(7)

式中：η为幅值标准差0.1的加性高斯白噪声。不同K值下，对信号f进行分解，各BIMF分量对应中心频率ω分布见表1。

表1 信号f在不同K值下的中心频率分布

表1指出，当K=4时，BIMF4的中心频率与K=5时BIMF5的基本相等。当K=5时，出现过分解，其中BIMF3和BIMF4的中心频率接近，所以K值取4。该实验为K值选取算法提供了依据，使用MIT-BIH数据库中的103心电信号仿真分析，如表2所示，K=4时，中心频率ω4几乎达到最大值，考虑到ECG信号90%的能量集中在0.05 Hz～45 Hz且P波(0 Hz～17 Hz)、QRS波群(0 Hz～33 Hz)与T波(0 Hz～9 Hz)等特征波形的频率分布，所以选择K=4。

表2 信号103在不同K值下的中心频率

1.3 R波定位

对于R波的定位，核心思想是：首先以变分模态分解为基础，确定预设尺度K值的自适应准则；K个不同的模态分量是原始信号在不同频域(从低频到高频)的表达，然后根据频域相关性提取出R波所在的有效分量，进行归一化并平方处理；最后通过平滑滤波获得包络峰值，标记为逻辑R波所在的索引位置。

逻辑上R波波峰仅仅可以初步缩小R波波峰所在位置的区域，不能做到准确定位。平滑滤波在峰值点的定位过程中存在时移偏差，需要附加峰值搜寻算法定位真正的R波波峰[14]，该算法需要在上一步归一化平方后、平滑处理前的第三模态信号区域峰值群(每个平滑包络峰值点的附近存在几个很陡峭的尖峰)中搜索最大幅值点，设定该复检规则可确保算法的准确性。可在提取到的平滑包络峰值前后特定窗口区选取N个采样点，最大的幅值点对应原始信号的索引位置即为真正R波峰值位置。故对应的R波定位表达式为

(8)

所以R波的后续定位算法如下：

(1)信号有效分量归一化平方处理；

(2)滑动平均滤波器对{N(BIMF)}2进行包络平滑处理；

(3)提取平滑包络的峰值得到R波的索引位置，标记为逻辑上的R波；

(4)在{N(BIMF)}2搜索包络峰值区域最大的幅值点，在平滑包络峰值的窗口Wt分别取N个采样点，取最大幅值点；

(5)寻找到的最大幅值点对应的即是R波在该模态分量中的索引位置，在原始信号中标记R波。

综上，本文的算法框架如图1所示。

图1 R波检测定位算法框架

2 信号仿真分析

为了验证所提出算法的有效性，本文首先选取MIT-BIH数据库中103信号做过程模拟。图2为103信号时域与频域图。

图2 103信号

1.2节有关K值选取的实验给出103信号的最优模态分解数为4，图3以及组图4为VMD分解给出的每个子信号与对应频谱。

图3 103信号VMD分解后的模态分量

图4 各模态分量对应频谱

分解后的BIMF1模态主要由输入信号的最低频率组成，BIMF3由信号的高频QRS区组成，T波和P波完全减弱，终模态BIMF4由心电其余高频成分或是50 Hz工频噪声组成，结合1.2节中所提到的ECG信号频率分布情况，提取BIMF3作为有效分量做R波定位分析。图5所示为BIMF3经过归一化再平方处理效果图。

图5 BIMF3归一化再平方处理

归一化平方处理能够突出信号幅值特征，放大大的值，缩小小的值。接着，将归一化平方后的BIMF3取平滑包络的峰值，定义为逻辑上的R波在模态分量上的索引位置，标记如图6所示。

图6 逻辑R波对应的索引位置标记

用滑动平均滤波器对归一化平方的BIMF3处理得到的平滑包络，其标记的逻辑峰值是将波形的每一处峰值群的几个尖峰平均处理得到的一个峰值集合，必然存在对于真正R波索引位置的微小时移，每一个标记的逻辑峰值附近存在不等的尖峰如图7所示。需要根据R波定位算法的复检规则，搜寻真正的最大幅值点标记。

图7 包络峰值的多个实际尖峰

一般人的心率绝对小于300 次/分钟，即在一定间隔内，只能出现一段QRS波，医学上称之为“不应期”。心电图学中，R峰的不应期在200 ms。MIT-BIH库中信号的采样率为360 Hz，以及QRS波群的时间宽度为100 ms，可设定在初始定位位置前后取50 ms窗口，依次取25个采样点[15]，标记最大幅值点，对应真正R波的定位位置标记如图8所示。

图8 103信号的R波定位

以上仿真实验复刻了算法的基本流程。信号103属于数据库中较为纯净的ECG信号，为了检验算法的鲁棒性，再选取库中几个典型的对于R波检测干扰较大的信号仿真验证，信号来源于108、109、113、209，分别对应特征为R波倒置、R波倒置加基线漂移、T波尖耸、强噪声。其定位效果如图9所示。

图9 具有不同特征的R波定位效果

3 结果与讨论

计算算法的灵敏度与准确率，需要多样本参考，本实验继续采用MIT-BIH库中样本仿真，主要针对5类典型ECG信号：正常窦性心律(NSR)、左束支传导阻滞(LBBB)、右束支传导阻滞(RBBB)、室性早搏(PVC)与房性早搏(APB)。

本文选用灵敏度(SEN)、阳性准确率(+P)以及准确率(Acc)评估算方法指标，引入真阳性(TP)、假阳性(FP)、真阴性(TN)、假阴性(FN)4个参数，分别代表正确检测出的R峰数量、误将尖峰噪声视为R峰的数量、正确判断为尖峰噪声数量、误将R峰视为尖峰噪声的数量，其计算式如下

(9)

(10)

(11)

考虑到真阴性(TN)指标，正确判断尖峰噪声数量的工作难以进行并且不是本文研究重点，所以可忽略该指标影响，然后结合文献[1，8，10，14]对评价指标的规定重新定义式(11)

(12)

样本信号采用30 min时程数据，对应的R波识别统计结果见表3。

表3 本文算法的R波识别统计结果

为进一步表明算法优势，将本文算法与参考文献中出现的算法比对灵敏度、阳性准确率及准确率指标(部分文献数据不足)，得到对照结果见表4。

4 结束语

对于原始ECG信号噪声强R波不便识别的问题，本文通过采用VMD技术对信号实现各特征频域的分解，排除噪声的干扰从而提取到R波所在的有效模态分量，对该分量进行分析能大大简化问题思路。接下来通过对分量信号归一化平方处理进一步滤除干扰，平滑滤波提取包络峰值初步确定R波位置，然后根据窗口采样方法校正R波索引位置从而使定位更加准确。借助MIT-BIH库中原始信号的支持，采用多样本仿真获取统计结果验证了算法的有效性。本文算法在实验仿真中表现出了高灵敏度高准确率的优势，但相对于传统算法计算量较大且过程相对繁琐，适合应用于对定位要求相对较高的领域。

表4 本文算法与相关文献算法结果对比