ECG ×AI: 机器/深度学习的ECG应用入门（3）

Aiwiscal

QRS心拍定位: 解决识别对象问题
1.引言
上一节我们解决了数据来源问题,成功将ECG数据读入了Matlab中.而我们现在还不能马上使用高大上的机器/深度学习算法,原因在于一个问题:在这个问题中,对于我们的机器/深度学习来说,直接处理的对象是什么?
是我们在上一节说的一条"记录"吗?对于大多数研究/论文来说,并不是所谓的记录,而是ECG的基本单位:心拍,就是在最开始的时候介绍的,包含了P波，QRS波，T波的信号段。那为什么不是“记录”呢？我猜测的原因可能如下：
1）一条充分长的记录中包含的信息太多，可能包含了很多种类型的心拍，直接全局处理不利于精细的识别。
2）一条充分长的记录数据量偏大，若一次性全局投入后端机器/深度学习算法，不利于复杂而又大量的迭代。
3）目前可用的开源数据库数据量偏少，若直接以一条记录为一个算法处理单位，则总体样本量大大下降，不利用机器学习，尤其是深度学习的应用。以MIT-BIH数据库为例，若以记录为单位处理，则只有48个样本；若以心拍为单位处理，样本量可达10W+。
暂时看不懂也没关系，总之记住，我们在这里直接处理的是心拍，而非一整条记录。
那么问题就来了，我们读到Matlab中的ECG数据都是以记录形式存在的，包含了形态各异，数量不等的心拍，例如下图给出了100，114，119信号的一部分，我们怎么样才能定位提取这些心拍呢？


这就引出了ECG算法领域的经典问题：QRS波定位检测问题。通过检测QRS波的位置，可以大体定位心拍所在的位置。

2.基本思想1）为什么是QRS波？
       直观来看，各类型心拍的各波段中，QRS波往往是最大、最明显的、最尖锐的，本身就利于检测；而在数学上，称为QRS波的“奇异性”，表现为斜率的突变，存在不可导的点，即波峰（或波谷）点。利用该性质可配合多种处理手法，如差分法，小波变换法等等进行有效检测。
2）基本步骤
       预处理+自适应阈值。预处理的目的在于两点：
       (1) 消除噪声和其他的杂波（P波，T波）等。
       (2) 同时使波形模式变得更为单一，QRS波的变化更为突出。
       预处理完成后，配合自适应阈值完成QRS波的定位。自适应的意思是能够适应不同类型的记录。自适应机制越好的算法，能够在各种类型的记录下良好地区分目标与非目标，从而达到很高的准确率和可用度。自适应阈值的设计原则，个人总结如下：
       (1) 自适应阈值变换机制中的经验参数需固定。在自适应阈值的变化机制中，有些参数可能需要我们根据经验提前设定好。这没有问题。但是一旦确定好之后检测性能时，对于不同的记录，这些经验参数不可以再重新设定。
       (2) 自适应阈值要跟随波形实时地，稳健地变化。自适应阈值，并不是至给一整条记录计算好一个阈值后就一成不变了，而是要跟随波形实时地变化，例如当检测到的QRS波振幅呈总体上升态势时，自适应阈值也要相应地提高，反之则要相应降低；自适应阈值的变化往往不可过于剧烈，可通过设定上下限的方式使其变化更加稳健。
      (3) 自适应阈值要有针对性地防错检、漏检的机制。例如可以回溯之前的波形，再次进行检测防止漏检；利用“不应期”（当一个QRS波出现后，短时间（~0.24s）内不会再出现第二个）的含义，在检测到一个QRS波后短时间内不再检出QRS波，防止错检。当然并不是说有了这些机制就完全可以避免错检漏检，而是有效地减少。
      以上几条关键点是本人自己总结的，可能有所疏漏，仅供参考。另外这些所谓目的和原则对于初学者来说都过于抽象，而且往往在阈值的设计中交织在一起。下面一部分我们将结合具体算法来进行讲解。

3.示例：QRS波定位算法设计与评估
在QRS波识别领域，不得不提的就是经典的Pan-Tompkins算法[1]和Cuiwei Li等人发表的小波变换法[2]。已经有大神实现并开源了Pan-Tompkins算法的Matlab代码，已经放到了本人github上：pan_tomkin.m。不过这里不打算以这两个算法为例讲解，因为这两个算法对初学者来说比较复杂。下面是以本人原创的一个QRS波定位算法进行分析。
1）预处理。这一部分的总体步骤图如下：

（1）我们对心电信号进行40阶FIR带通滤波，通带为15~25Hz，大致为QRS波所在频段。关于滤波器的设计可以使用Matlab的Fdatool工具，得到相应的滤波器系数，然后套用filter内建函数完成滤波。以数据库的103号信号为例，滤波前后的波形示意如下：

通过对比可以看到，滤波后P波，T波被明显削弱，体现了上一部分中预处理目的的第（1）点。

（2）对滤波后的波形“双斜率”预处理。所谓“双斜率”处理是参考自另一篇文献[3]。基本思想是分别在一个点的左右两侧的某个区间内寻找最大平均斜率与最小平均斜率，然后分别用左侧最大斜率减去右侧最小斜率，用右侧最大斜率减去左侧最小斜率，再求取两者中的最大者。过程有些复杂，但其实其基本动机就是利用QRS波两侧较陡的性质，只有QRS波这样的尖峰在经过上述处理时才会有很大的响应。在这里，我们设定这个寻找斜率的区间是左右两侧0.015s~0.060s处，此为经验参数。经过处理，可对比前后波形：

可以看到，双斜率处理后，波形模式更为单一，体现了预处理目的的第（2）点。但是，波形出现了双峰现象，一定程度上不利于精准检测。因此，在双斜率处理后继续低通滤波（截止频率5Hz，经验参数），使得波形更光滑：

可以看到，低通滤波后的波形变得非常光滑，杂波基本消失，模式也非常单一，基本上已经达到了进一步阈值处理的要求。
（3）滑动窗口积分。仔细观察前几步，每一步得到的波形由于滤波或是求斜率的因素，其幅度值越来越小，而过小的幅值其实不利于检测。这里我们利用滑动窗口积分，使得绝对振幅增大，并使波形进一步光滑，滑动窗口宽度设为17个采样点，为经验参数：



积分后，波形幅值明显增大（注意纵坐标），体现了预处理目的的第（2）点。至此，预处理阶段完成，经过处理后，原始的心电信号变成了一个个模式单一的波峰组成的信号，每个波峰对应一个QRS波。相比原信号，预处理后的信号更易定位检测，可以说达到了预处理的目的。

2)自适应阈值设计
        与上述所说的自适应阈值设计原则相契合，这个算法的思路是：
        （1）自适应阈值要跟随信号实时变化。问题是如何决策此时阈值的变化。这里我采用了一种常用的策略：双阈值。一高一低，当某个波峰超过低阈值时，我们认为检测到了一个QRS波，然后通过比较波峰振幅与高低阈值的关系，调整阈值。
        （2）为了保证变化的稳健，阈值需要根据前面已经检测到的正确波峰振幅变化，另外设定下限，以确保不会过高或过低。
通过以上两点，我设计自适应机制如下（THR1为高阈值，THR2为低阈值）：



其中，peak指检测到的峰值，peak_buffer指存储了当前峰值之前8个峰值，mean(.)则是求均值。THR1_lim与THR2_lim为两个经验常数，代表两个阈值变化的下限，在这里分别取值为0.3和0.23。
具体来说，当当前波峰高于高阈值时，我们倾向于认为此时阈值偏低，阈值变化为之前检测到的8个波峰平均值的0.7倍，0.25倍（均为经验系数），如果相比最近的8个峰值，阈值确实偏低程度很大，则就算乘上了一个小于1的系数（例如0.7，0.25），自适应阈值还是会增大，如果偏低程度不大，则不一定增大。之所以乘上系数还是因为要保持变化的稳健，保证只有阈值明显偏低时才会明显增大；当当前波峰在高低阈值之间，我们倾向于调低高阈值，而低阈值要与高阈值保持距离，因此也倾向于调低一些，不过我们倾向于认为此时的波形较为异常（因为较小，在两阈值之间），不宜根据前面的平均峰值变化，因此此时阈值变化与当前波峰的值相关。具体变化机制可以参考上面的公式，体现实时变化原则。另外重复一下，当波峰低于低阈值时，认为是噪声。
         （3）防漏检与错检。其实本身双阈值的设计一定程度上就可以防漏检。因为一高一低的阈值，相比只有一个阈值，可以捕捉到更多层次的波峰。另外在防错检方面，我采用了“不应期”的方式，当两个波峰过于靠近时，只取较大的波峰。这个“过近”指的是距离低于0.24s，即不应期的长度；双阈值的变化存在下限，也防止了一些噪声被错检为QRS波。
另外，双阈值的初始值是被提前设定好的经验系数，与其他经验系数一样，一旦确定就不会因记录的不同而变化。体现了上述自适应阈值设计的第（1）原则。

matlab的旋律

以上阈值的操作，均在原信号预处理后的特征信号上进行，检测到的波峰都认为对应于原信号中的一个QRS波。另外预处理中的两次滤波都造成了信号的延迟，我们检测到特征信号的波峰后，可减去延迟，即可得到原信号中的QRS波的位置。

说了这么多确实有些抽象，还是看一下直观的图：

上图这是以100号信号为例。由于本身信号比较稳定，因此阈值的变化主要在信号开始。可以看到初始时阈值明显偏低，自适应阈值很快向上变化找到了合适的值。下面再看一条变化比较大的信号119：

可以看出双阈值可以很好地跟随信号变化，从而正确检出QRS波的位置。注意,这里检测的"正确"并非是指精准地落在QRS波的波峰或波谷中,尽管精确的人工标注往往就是在波峰或波谷上.美国EC38标准提供了一个评判标准:与人工标准差距在150ms内,即可认为定位成功.事实上,这个标准相当宽松,以至于现在的QRS波定位算法的准确率普遍在99%以上(以MIT-BIH数据库为测试数据), 而一些出色的算法甚至达到了99.9%以上.
以上原创算法已开源,放在本人的github上: DS_detect.m, 已封装为函数。以上的算法可能理解起来相当抽象，建议结合代码（有注释）体会。


3）算法评估
上一部分中曾说过，目前QRS波定位检测算法的“准确率”已经普遍在99%以上。但这么说其实不准确，因为衡量QRS波检测算法的指标并非我们印象中的准确率，而是有两个常用指标：敏感度（SE）和正预测率（P+），定义如下：

其中TP代表检测正确的心拍个数，FN代表漏检的心拍个数，FP代表错检的心拍个数。我们尝试对整个MIT-BIH数据库跑一遍，统计算法的性能参数。这一部分不是算法的核心，所以在这里不再讲具体的实现步骤，不过已经开源了测试代码S_test.m,为Matlab脚本文件。在运行之前注意数据库路径的设置。注意在代码的第164行，调用了上面的算法DS_detect。

DS_detect和DS_test最好在一个文件夹中，或者将这些文件放到Matlab可以检索到的目录中。以上都没有问题后，直接运行DS_test,等待其完成。当然也可以把DS_detect改为其他算法，只要搞清楚其输出哪个是QRS波的索引就可以。DS_test程序完成后，在工作空间中可以找到名为INFO的变量，双击打开：

从左到右，第一列：记录编号，第二列：人工标记心拍数，第三列：算法正确检出心拍数，第四列：算法错检心拍数，第五列：算法漏检心拍数，第六列：Se，第七列：P+。看最后一行：

可以看出算法共错检648个心拍，漏检380个心拍，总体Se达到了99.65%，P+达到了99.41%，汇总如下：

在我的笔记本电脑上，该算法处理一条30min，采样率360Hz的记录需要14s左右。各项指标虽达不到顶尖，但也不算差，算是有一定的可用度了。接下来，我们将使用该算法将我们需要的心拍截取出来。

4.心拍截取这里其实有一个矛盾。无论我们的QRS波定位算法有多厉害，都很难100%将人工标注过的所有心拍定位准确，必然有错检和漏检。这就导致我们无法100%利用好人工标记的心拍。这里我们采用一个折衷的手段：忽略漏检的心拍，而错检和正确检出的心拍都截取出来。由于QRS定位算法的、性能已经够好了，所以截取出来的心拍绝大多数是正确检出的心拍，只有少部分为错检心拍。而我们接下来所要使用的机器/深度学习算法都具有一定的鲁棒性，较少的错检心拍数不会对结果产生太大影响。至于心拍的类型（标签），我们以距离该心拍最近的人工标记的类型为准。
另外，我们截取心拍其实不是精确截取，即从P波起点到T波终点，而是采取一种较为粗略的手段：以QRS波的位置为基准，分别向前向后包括若干点，然后将这一段数据点截取出来作为心拍。如下图所示为什么不精确截取呢？是因为检测P波和T波及其起始点或终点非常困难，目前也没有一个性能足够好的算法满足应用需求。所以，大部分文献采用的是这种粗略截取的方法，只不过前后包括的点数有区别：

这里我们向左包含100个点，向右包含150个点，即截取的每个心拍长度为250个点（约0.7s）这一部分的实现代码也已开源：SegBeat.m 为脚本文件。同样，数据库路径，代码文件存放位置要正确。该代码与DS_test.m有类似之处，也调用了DS_detect算法进行QRS检测（第164行，同样也可以置换为其他QRS波定位算法）。至于其他实现细节，作为初学者可不必了解，因为确实很绕，并且也并非是我们的主线任务。运行完成后（时间可能会比较长），在工作空间可以找到4个变量，Nb,Lb,Rb,Vb,为我们接下来的目标类型心拍数据集合，分别为“正常（N）”，“左束支阻滞（LBBB）”，“右束支阻滞（RBBB）”，“室性早搏（PVC）”，示例如下（可能横坐标表示有些不准确，请忽略。。。。）：

实际上，同一类型的心拍间变化非常大，并非只是示例的四个心拍那么简单，这也就是我们为什么后面要使用高大上的机器/深度学习算法来分辨。我们以.mat的形式将这四个变量保存（N_dat.mat,L_dat.mat,R_dat.mat,V_dat.mat），它们各自包含了一种类型的心拍，是我们接下来使用机器/深度学习算法要进行分类识别的目标。至此，我们所需要的识别对象已经成功准备好了。

5.小结

本节内容很多，也很抽象。实际上，QRS波检测定位算法的设计很大程度都要靠经验，只不过个人认为有了上面说的几条原则的指导下，方向可以更明确。相关代码已经开源,大家可以结合代码来加深理解,多上手,多尝试.把目标心拍截取之后,我们接下来就可以进行机器/深度学习的应用了.

** 相关代码文件下载：https://github.com/Aiwiscal/ECG-ML-DL-Algorithm-Matlab

** 推荐文献：[1]  Pan J, Tompkins W J. A real-time QRS detection algorithm[J]. IEEE transactions on biomedical engineering, 1985 (3): 230-236.

                   [2]   Li C, Zheng C, Tai C. Detection of ECG characteristic points using wavelet transforms[J]. IEEE Transactions on Biomedical  Engineering, 1995, 42(1): 21-28.

                        [3] Wang Y, Deepu C J, Lian Y. A computationally efficient QRS detection algorithm for wearable ECG sensors[C]//Engineering in Medicine and Biology Society, EMBC, 2011 Annual International Conference of the IEEE. IEEE, 2011: 5641-5644.