应用于纳米刀的心脏有效不应期实时检测算法研究∗

董腾飞，潘 松，吴 松，郝珠珠，赵伟杰，杨 勇

(杭州电子科技大学自动化学院，杭州310018)

基于不可逆电穿孔技术(Irreversible Electroporation，IRE)的纳米刀在治疗肝癌、胰腺癌以及其他恶性肿瘤的场景下发挥着越来越强大的作用[1-2]。但是IRE治疗中释放的高压脉冲对心电信号的产生和传导会产生极大干扰，增大了引发心律失常等术中并发症的机率，故需要对患者的心电进行实时监测[3]，并使IRE治疗脉冲在心脏的有效不应期内释放，减小对心脏的正常工作过程的影响。故能否准确识别患者心脏的有效不应期，及时触发纳米刀释放高压脉冲，对降低术中并发症有着重要的意义。

由于在术中无法直接对心脏的有效不应期进行检测，本文根据体表心电图和有效不应期的对应关系，将心脏有效不应期的实时检测问题转化为R波实时检测问题。但由个体差异导致的心电波形变化，以及由高压脉冲刺激引起的骨骼肌强直收缩或心率失常等情况，会使采集到的心电波形与正常的心电波形产生较大的差异，从而给R波的准确识别带来困难；此外，通过检测心脏有效不应期来指导IRE治疗脉冲的释放时机，要求输出R波检测结果的实时性较高，算法的计算量较小。

目前，国内外常用的临床心电信号R波检测方法有硬件方法和软件方法两种。硬件方法响应速度快、实时性好，但是缺乏灵活性，鲁棒性较低。软件方法有时域分析法[4]、模板匹配法[5]、小波变换法[6]和神经网络法[7-9]等。时域分析法是应用最广、实现最简单、计算量最小的心电信号分析方法之一，但是变异的大P波、大T波，无法被滤波器滤除干净的噪声都有可能会降低检测结果的准确率，故该方法还有进一步提升的空间；小波变换是另一种应用较广、研究较深的方法，但是计算量大，在实际应用过程中软硬件设计复杂，实现困难；神经网络法多用于波形识别、病例诊断方向且计算量大。

以上几种方法并不能满足心脏有效不应期检测高准确率、高实时性、低计算量的要求。本文在时域分析法的基础上，从预处理步骤和R波识别步骤两个角度出发，以提高实时性和准确率、简化计算量为目标对算法加以改进创新，设计出了一套可以应用于IRE治疗的R波实时检测算法，从而达到识别延时低于15ms，且当识别到含有大量噪声的心电数据段时自动停止R波检测，提高算法鲁棒性的要求，进而通过固定延时实现了心脏有效不应期的实时检测定位，指导IRE治疗脉冲释放时机。

1 心电信号预处理

心电信号是微弱的低频信号，对外界干扰比较敏感[10]，故需要采用一定的预处理步骤滤除噪声，突出目标信号特征，简化计算量。这是准确分析心电信号的重要前提，本文采用的预处理流程如图1所示。

图1 预处理流程图

1.1 频率归一化处理

MIT-BIH心率失常数据库[11]的采样频率为360 Hz，QT数据库[12]的采样频率为250 Hz，采用频率不同会影响算法的预处理效果，故采用重采样的方法归一化为500 Hz。

1.2 滤波处理

由于本算法实时性要求较高，且只需要检测R波，故在保证检测效果的前提下，可以选用低阶IIR滤波器级联的方法进行滤波，降低有效信号频段通过滤波器的群延时。

常见心电信号噪声主要由肌电干扰、基线漂移、50 Hz工频干扰组成，本算法的滤波部分主要针对这几类噪声进行设计，滤波器参数如表1所示。

表1 滤波器参数

已知心电信号的主要能量集中在0～40(±7)Hz，QRS复波的主要能量在 6 Hz～18 Hz，使用MATLAB的fdatool工具分析得知，心电信号经过低通滤波器后，0～40 Hz的频率分量群延时约6 ms；经过高通滤波器后，5 Hz以上的频率分量群延时低于5 ms；经过带阻滤波器后，40 Hz以下的频率分量群延时低于3 ms。所以经过滤波处理后，QRS波段的群延时理论上不超过15 ms，且波形不会发生较大的非线性失真，满足了R波识别的实时性要求，为进一步的预处理做好了准备。

1.3 微分处理

本文使用微分函数来突出QRS波群在正常心电信号中斜率最大这一特征，将QRS波群从其他心电信号中区分出来。微分函数为:

式中:sig_f为经过滤波器预处理后的心电数据；sig_d为对sig_f进行微分处理后的数据。

1.4 低阈值归零处理

由于本文只需要检测R波，故可采用低阈值归零法对其进行处理，突出R波微分值较大的特征，同时滤除部分噪声干扰，简化算法需要处理的数据量，提高处理速度。低阈值归零处理的公式为:

式中:sig_z为低阈值归零处理后的数据；Thr_z为归零阈值。

1.5 平方处理

本文采取平方处理对信号进行非线性滤波，进一步突出R波微分值较大的特征，减弱其他波形或噪声的干扰。

式中:sig_s为sig_z平方后的结果。

1.6 移动窗口平均处理

对于正常心电信号而言，经过以上步骤的处理，R波分量已经非常突出，但当原始信号的信噪比较低或噪声的频段与心电信号的频段相重叠时，会导致滤波效果不佳，影响R波检测的准确率。

本文利用移动窗口平均法的特性将每个心跳周期分为有效信号区段和疑似噪声区段，并分别对其进行处理，达到提高检测结果阳性预测度的目的。移动窗口平均的运算公式为:

式中:sig_w为移动窗口平均后的数据；N为移动窗口的窗宽。当窗宽设计过大时，会导致噪声检测过于灵敏，将含有正常R波的有效信号区段误判为疑似噪声区段而造成漏检；当窗宽设计过小时，会使得噪声检测不够灵敏，将疑似噪声区段误判为有效信号区段而造成误判。正常人心电信号QRS波段最长不超过0.12 s，故窗宽时长选为0.15 s，采样频率为500 Hz时N的值为75。

2 R波检测算法

本文通过对MIT-BIH心率失常数据库和QT数据库的数据进行分析，设计了一套完整的R波实时检测方案，同时对部分变异心电信号和含有大量噪声的心电信号进行特殊处理，提高检测的阳性预测度。R波检测算法的整体框架图如图2所示。

图2 R波实时检测算法流程图

2.1 阈值初始化

当sig_f的长度超过5 s时，以1 s为间隔将其分为5段，保证每段中至少包含一个QRS波群，然后分别求最大值，再取中间大小的三个值求平均，记为SPK_f，根据SPK_f计算幅度阈值Thr_f，计算公式为:

微分阈值Thr_d1、Thr_d2和移动窗口平均阈值Thr_w的计算方法与上述过程相同，计算公式如下:

式中:SPK_d为对前5 s微分数据sig_d处理后的平均值；SPK_w为对前5 s移动窗口平均值sig_w处理后的平均值；a为比例系数，一般选取a＝0.035。

低阈值归零处理的归零阈值Thr_z计算公式如下:

式中:b为比例系数，一般选取b＝0.3。

2.2 噪声及室颤检测

本算法提出根据移动窗口平均后的数据sig_w与阈值Thr_w的关系，将每个心跳周期分为有效信号区段和疑似噪声区段。

第一阶段为有效信号区段，噪声标志位为0，当sig_w持续大于Thr_w的时间不超过150 ms时，进行R波检测，否则进入疑似噪声区段。程序运行框图如图3所示。

图3 有效信号区段程序流程图

第二阶段为疑似噪声区段，噪声标志位为1，当sig_w持续大于Thr_w的时间超过500 ms时，表明心电信号中含有较大噪声或患者出现心率过快、室颤等病症，此时为了防止误判暂停漏检机制，直到检测到下一个R波；如果持续时间超过5 s，则启动纠正机制；当sig_w持续小于Thr_w的时间超过150 ms时，返回第一阶段，开始R波检测。程序运行框图如图4所示。

如图5(a)、5(b)所示，本算法能够有效识别MIT-BIH心率失常数据库104号数据中含有噪声干扰的数据段和205号数据中含有异常心电波形的数据段，并自动停止检测，待心电信号恢复正常后，再重新开始检测。

图4 疑似噪声区程序流程图

图5 噪声数据

2.3 针对 QS、Qr、rS、rSr波形的处理

在一组数据中，若所有的心电波形都是如图6所示的 QS、Qr、rS、rSr′波，可以采用更换导联的方法显示向上的qRs波波形；若仅有部分心电波形是此类波形，则表明患者出现早搏或其他症状，需要重新判断是否适合进行IRE治疗，且根据这类心电波形来判断心脏有效不应期位置的准确率较低，故当检测到此类波形时，暂停R波检测。

算法实现如下:若存在某一时刻的心电数据同时满足式(10)和式(11)，且该时刻距离上一个R波大于200 ms，则判断出现此类心电波形。根据心脏有效不应期原则，同时为了避免将此类异常波形之后出现的大T波误判为R波，采取延时360 ms后再继续检测R波，且暂停漏检机制的策略，效果如图7所示。

图 6 QS、Qr、rS、rSr心电波形

图7 233号数据

2.4 平均RR间期

取最近8个RR间期的平均值为RR_mean，计算公式为:

式中:RR(i)为两个R波峰之间的采样点数。

2.5 R波检测

本算法在幅值、微分值两个维度设立动态阈值检测R波，方法如下:①如果连续四个点的微分值sig_d大于阈值Thr_d1，且该点距离相邻的R波超过200 ms，则表示检测到一个R波上升沿；②存在一点的微分值sig_d小于0，幅值sig_f大于阈值 Thr_f；③从满足条件1的点开始到找到一个满足条件2的点为止的这段时间内，任何一点的幅值和微分值都在上一个R波幅值和微分值的0.4倍到2倍之间。

当某一时刻的检测结果满足条件1，且50 ms内存在一点能够同时满足条件2和条件3，则判断找到一个新的R波波峰。其中，条件3是针对心电信号采集过程中突然出现的接触噪声或其他能够引起心电采集系统检测到信号突变的情况而设立的。

由于QRS波的特征参数会随着生理情况变化而变化，故采用动态阈值法提高算法鲁棒性。在检测到R波后同步更新阈值 Thr_f、Thr_d1、Thr_d2、Thr_w，动态阈值更新公式如下。

式中:sig_f_max、sig_d_max、sig_w_max 分别为检测到R波上升沿到检测到R波波峰这段时间内sig_f、sig_d和sig_w的最大值。

2.6 漏检机制及纠正机制

漏检机制:由于实时性要求较高，本文并未采用具有较大延时的修正和回检算法提高检出率，而是通过动态阈值法减少漏检。具体算法如下:

如果持续1.66倍RR间期没有检测到R波，则将微分阈值Thr_d1降低一半至Thr_d2，重复1.3.5节的算法继续进行检测，检测到R波后更新阈值Thr_f、Thr_d1、Thr_d2、Thr_w。 其中 Thr_f、Thr_d1、Thr_d2计算公式如下，Thr_w按照式(19)和式(20)进行更新。

纠正机制:如果持续5 s没有检测到R波且噪声标志位为1，则停止检测并提示原始信号信噪比过低；如果持续5 s没有检测到R波且噪声标志位为0，则暂停检测并自动重新初始化阈值，再利用新的阈值进行检测。

3 有效不应期检测

心脏的有效不应期是指使用3～5倍的阈刺激值对心肌细胞或组织进行刺激时，不会引起兴奋反应的一段时间，约200 ms～300 ms，对应于单细胞动作电位的0～2相和3相的前部，以心室肌为例相当于从体表心电图的QRS波开始一直持续到T波的前支。本文所述心脏有效不应期检测特指根据体表心电图检测心室肌有效不应期。同时，在心房肌、心室肌相对不应期刚开始的一段的时间内，应用较强的阈上刺激会有较大概率引发心房或心室颤动，这段期间被称为易颤期[13-14]。在对人体施加高于阈刺激值的电刺激时，应当避开心房肌和心室肌的易颤期。

我们可以根据体表心电图和有效不应期的对应关系，在检测到R波后延迟一段时间避开心房肌和心室肌的易颤期，准确地将IRE治疗脉冲的释放时间点控制在心脏有效不应期内。通过对算法延时、个体差异等因素的综合考虑，延时时间可以设置为50 ms～100 ms。

单细胞动作电位和心电图与不应期的对应关系如图8所示[15]，图中两块灰色区域分别表示心房肌易颤期和心室肌易颤期；a点表示R波波峰的位置；b点表示经过计算和固定延时后释放IRE治疗脉冲的位置。

图8 单细胞动作电位和心电图与不应期的对应关系

4 实验结果分析与讨论

4.1 数据库信号仿真实验

QT数据库标注了各个波形的起止点和峰值点，最适合用来进行波形定位准确性的评估，但部分标注的位置有较大误差，本文根据标记点位置的准确性，人工挑选出其中67组数据作为数据源进行仿真。

4.1.1 评价标准

4.1.1.1 准确性的评价标准

为了分析R波实时检测算法的定位准确性，本文采用如下几个参数:0257真阳性(TP):准确检测出的R波峰；0258假阳性(FP):非R波峰误判为R波峰；0259假阴性(FN):漏检的R波峰。

因为本算法要求检测延时低于15 ms，故取消了具有较大延时才能实现的R波修正和回检机制，且设计思路之一是通过过滤部分非标准R波心拍以及噪声较大的数据段以提高算法的鲁棒性，所以会出现部分R波漏检的情况，漏检数量由信号采样质量以及非标准R波心拍的数量决定，故不采用假阴性FN及其相关指标作为定量判断标准。

理论上本算法在R波峰后的15 ms内输出检测结果，但由于数据库官方标记位置有略微的误差，本文设定在检测结果之前20 ms，之后10 ms这一区间内存在官方R波标记点即为真阳性TP，反之为假阳性FP。

最终，本文通过计算阳性预测度(＋P)来定量评估算法，阳性预测度计算公式如下:

4.1.1.2 实时性的评价标准

在分析实时性时，本文采用平均延时(Td)来表示输出检测结果的延迟时间，并采用如下两个指标对该参数进行分析:

①均值:

4.1.2 仿真结果及分析

采用MATLAB进行仿真的结果如表2所示。

表2 QT数据库R波检测结果

续表2

根据表2计算平均阳性预测度、平均延时和方差的结果如表3所示。

表3 仿真结果分析表

对以上两表分析可知，仿真结果的平均阳性预测度为99.88%，最高为100%，最低为98.37%；平均延时为13.02 ms，最小为9.49 ms，最大为14.67 ms，达到了实时性和准确性的要求。

4.2 实测信号检测实验

将该算法移植到自主开发的心电信号采集分析系统，并采用SKX-2000G＋型心电信号模拟仪和人体实测心电作为数据源进行检测，实验结果如下:

4.2.1 使用心电信号模拟器作为数据源的测试结果

如下图所示，黑色实线代表心电波形；虚线代表不应期检测结果，虚线的上升沿代表检测到R波的时刻，下降沿代表理论上IRE治疗脉冲输出的时刻，该时刻位于心脏的有效不应期内。

①正常心电和窦性心电。可以100%的准确检出R波，实现不应期定位。

图9 正常心电

图10 窦性心电

②含有接触噪声。在含有大量噪声的时间段内，自动停止检测R波，噪声结束后再重新开始检测R波。

③叠加50 Hz工频干扰的心电信号。由图可见，50 Hz工频干扰被完美滤除，不影响R波的正常检测效果。

④高大T波。当出现高大T波时，能够准确检测出R波的波峰。

⑤倒尖角波形，模拟 QS、Qr、rS、rSr这类波形，可以完美滤除，不进行R波识别。

图11 含有大量噪声的心电

图12 叠加50 Hz工频干扰的心电

图13 大T波

图14 倒尖角波形

4.2.2 使用心电信号模拟器作为数据源的测试结果

利用心电信号采集分析系统对人体的心电信号进行采集识别，结果表明本算法不但可以有效的对正常人体心电进行识别，而且当心电信号含有较大的运动伪迹或其他噪声干扰时，同样可以起到较好的识别效果，具有较强的鲁棒性。

图15 正常心电

图16 含有运动伪迹

图17 含有大量噪声的心电

如上图所示，虚线的上升沿对应于R波下降沿顶端的位置，表明本算法可以在极小的延时内识别到R波，达到了实时性要求；折线的下降沿一般在S波上升沿过后的位置，表明本算法可以避开易颤期，准确的将IRE治疗脉冲释放的时间控制在心脏有效不应期内，达到了较好的效果。

5 总结

针对IRE治疗需要实时检测患者的心脏有效不应期这一需求，本文在时域分析法的基础上进行改进和创新，实现了复杂噪声条件下的R波实时检测功能和心脏有效不应期定位功能，平均阳性预测度达到99.84%，平均延迟在15 ms以内，且能够准确识别并排除部分变异心电信号和噪声干扰。对比于传统心电信号检测算法，在保证准确率的同时，实时性有了明显的提高，与杨衍菲等人针对除颤器设计的 QRS波斜率检测方法相比[16]，平均延时从17.3 ms降低到了13.02 ms。

本文算法测试环境与实际IRE治疗时的心电信号略有差异。首先，纳米刀消融的高频脉冲会为心电信号引入一些噪声。关于这一问题，本团队已经开展基于动物实验的相关研究，后续进展将在下一篇论文中呈现。其次，即使在精确识别心脏有效不应期的情况，外加全身麻醉及深度神经肌肉阻滞的情况下，纳米刀手术过程中仍有可能由于高压脉冲的作用而偶然引发自限性心律失常[17]。针对这一问题，算法设计检测到非正常心电信号后，会发送指令停止一段时间的检测，以暂停纳米刀的输出。当专业医疗人员观测到心电信号回归正常可以继续手术后，再次启用算法进行检测。

本文算法可用于指导临床纳米刀消融手术过程，降低对心脏正常生理活动的影响，提高手术安全性。下一步可以使用心电信号采集分析系统结合高频脉冲发生器展开相关动物实验，对本算法和心电信号采集分析系统做进一步的改进和优化。