心电图是临床上医生诊断心脏疾病的一个常规而且高效的技术手段, 心电图各波与心肌动作电位有着密切关系. 心肌梗死根据病变程度不同, 在心电图上可以体现出不同的形态. 目前在心肌梗死的检测上, 多采用传统的信号处理方法, 如傅里叶变换, 小波变换等, 该传统方法对于波形质量要求高, 对于有噪声干扰和形态不好的波形效果不佳. 也有一些学者提出使用SVM^[1], KNN, 决策树等机器学习方法^[2], 这一方法前期依赖于传统方法, 需要提取相关特征, 对于非典型性心肌梗死的检测^[3], 由于特征不明显, 故准确率不高. 目前也有学者单纯使用卷积神网络诊断心肌梗死, 但该方法对于典型心肌梗死的检测准确率没有直接提取特征的效果好, 并且时间复杂度较大. 由于非典型心肌梗死波形复杂多变, 特征模糊不易提取, 所以BP神经网络在非典型心肌梗死的检测上表现一般.

本文针对上述问题, 综合考虑准确率和算法时间复杂度因素, 提出BP和CNN结合的方法. 该方法在非典型性和典型性心肌梗死的检测上都具有较好效果. 首先对ECG信号做前期处理, 主要包括滤波, 去基线, 差分等. 之后针对典型心肌梗死, 对各个波做准确定位和识别, 再根据心肌梗死的典型特征, 提取出相应的特征值, 采用BP神经网络训练; 对于非典型性心肌梗死, 由于特征模糊, 在对ECG信号前期处理后, 采用卷积神经网络训练^[4]. 算法整体流程如图1所示.

1 研究方法

本文主要从数据获取, 数据前期处理, 波形检测, 特征提取, 构建神经网络及识别等几个方面叙述.

1.1 数据获取

数据主要来源于MIT-BIH数据库和中国医科大学附属第一医院.

1) MIT-BIH是国际上公认的心电数据库, 由美国麻省理工大学提供. 在其ECG Database中包含PTB Diagnostic ECG Database数据库, 其中有268个病人数据, 心肌梗死的共有148人. 该数据库以其固定格式编码存储, 需要特定工具处理. 每个病人数据中包含.dat(ECG数据), .hea(病人详细信息)和.xyz(Frank导联数据)三种格式文件, 采样率为1000. 除了常规的12导联外, 还有vx, vy, vz三个Frank导联, 共15个导联.

2) 中国医科大学附属第一医院提供了从2008年到2015年间心电科的病人数据. 从中提取出所有的心电数据, 并对心电数据按医生诊断结果进行详细分类, 提取出所有患有心肌梗死的病人的心电数据, 并进一步核对数据. 其中, 诊断为心肌梗死(包含疑似心肌梗死)的病人481人. 该数据为xml格式, 需要从中提取信息, 文件包含病人姓名, 性别, 年龄, 诊断结果等信息. 采样率为500, 包含常规12导联数据.

1.2 ECG数据前期处理

对于MIT-BIH需要使用WFDB工具包根据不同数据格式提取相应数据. 对于医大一院数据, 需要从xml中根据关键字段提取ECG数据和诊断结果. 之后根据二者不同采样率和幅值设置不同的滤波系数等. 数据先经过带通滤波器去除一定的噪声干扰, 之后经过改进的均值滤波平滑后, 保证数据形态特征不发生较大改变. 然后, 将数据做去基线处理, 防止基线漂移. 最后将数据做放大处理, 以便后期提取特征值.

1.2.1 滤波

本文根据滤波效果, 波形特点和对时间复杂度的要求尝试了多种滤波方法, 最终比较了均值滤波和butterworth滤波.

(1) 均值滤波如下:

假如滤波器移动平均窗宽为5, 经过平滑处理返回值yy, 则yy为:

$ yy\left( 1 \right) = y\left( 1 \right) $

(1)

$ yy\left( 2 \right) = \left( {y\left( 1 \right) + y\left( 2 \right) + y\left( 3 \right)} \right)/3 $

(2)

$ yy\left( {\rm{n}} \right) = \frac{1}{5}\mathop \sum \nolimits_{i = n - 2}^{n + 2} {y_i}\left( {n \ge 3} \right) $

(3)

(2) Butterworth滤波如下:

$ {\left| {H\left( {jw} \right)} \right|^2} = \frac{1}{{1 + {{\left( {w/{w_c}} \right)}^{2n}}}} = \frac{1}{{1 + {\varepsilon ^2}{{\left( {w/{w_p}} \right)}^{2n}}}} $

(4)

其中, n是滤波器的阶数, w_c是截止频率, w_p是通频带边缘频率, $\dfrac{1}{{1 + {\varepsilon ^2}}}$ 是|H(w)|²在通频带边缘的数值.

如果令s=jw, 拉普拉斯变换在虚轴s=jω上的性质可以得到:

$ H\left( s \right)H\left( { - s} \right) = \frac{1}{{1 + {{\left( {s/j} \right)}^{2n}}}} $

(5)

由极点在单位圆上可以得到:

$ {s_k} = {e^{\frac{{j\pi }}{{2n}}\left( {2k + n - 1} \right)}},\;\;k = 1,2,3, \cdots $

(6)

因此传递函数为:

$ H\left( s \right) = \frac{1}{{\left( {s - {s_1}} \right)\left( {s - {s_2}} \right) \cdots \left( {s - {s_n}} \right)}} $

(7)

1到10阶的Butterworth多项式因子可以通过查表得到.

图2–图4为两种方法的对比结.

从图中可以看出, 经过均值滤波后, ECG波形基本没有发生形变, 只是随着平滑点数的增加, 幅值有所降低, 为此采用一种补偿方法, 将其补偿到原来的数值. 但是经过Butterworth滤波后, 除了波形幅值有所降低外, 在波形拐点处会有一定的形变, 尤其在S波会有一个明显向下的尖小的波, 因此导致波形在一定程度上有所失真.

1.2.2 去基线

由于存在基线漂移, 对提取特征值造成很大干扰, 尤其是心肌梗死判断中, 准确识别出ST段是否抬高或压低是及其关键的, 如果存在基线漂移, 就不能准确识别这一重要特征. 在一个周期中采用插值方法找到基线, 然后将所有数据减去基线, 得到去完基线的ECG信号.

1.3 QRS波群检测

理论上每一个完整节拍ECG信号都有P波, QRS波群, T波组成. 由于在一些特殊情况下, 可能出现P波消失以及T波不明显的情况, 但都会有QRS波群, 所以, 我们先检测QRS波群, 然后, 由QRS波群确定P波和T波. 文献[5]提出一种改进的DWT算法来做波形匹配. 首先, 对波形求一阶导数, 得到QRS波群的大概位置序列, 之后在对原始波形求二阶导, 找到更为精确的位置, 然后根据位置序列的索引范围从原始波形中找到QRS波群. 之后, 在经过滤波等处理后的波形中按照索引范围, 以及原始波形高度的阈值等多维特征共同确定QRS波群位置, 并对其修正. 图5–图7分别为滤波后原始数据, 一阶导数据和二阶导数据.

由图可以看到, 原始波形中QRS波群波峰的位置索引靠近一阶导峰值位置, 理论上在一阶导数值为0对应的位置上, 而一阶导为0的位置变化率最快, 因此与二阶导谷值对应索引更加接近, 故选择在二阶导谷值附近查找原始信号的QRS波峰. 由QRS波群进而定位P波和T波.

1.4 ECG特征提取

在准确识别了QRS波群的前提下, 需要对ECG波形做进一步的特征提取. 国际上公认的常规12导联是标准I, II, III, 加压单极肢体导联aVR, aVL, aVF和单极胸腔壁导联V1–V6^[6]. 12导联之间存在如下向量运算关系:

$ aVR + aVL + aVF = 0 $

(8)

$ II = I + III $

(9)

$ aVR = - 0.5*\left( {I + II} \right) $

(10)

$ aVL = I - 0.5*II $

(11)

$ aVF = II - 0.5*I $

(12)

每一个单极导联都由周期性的P, QRS, T波等组成. 个别还有U波, 但目前没有发现明确的生理学意义(U波异常可能与缺钾有关, 但不是绝对的), 我们暂不研究. 一个完整心动周期如图8所示.

其中, P波代表心房肌除极的电位变化, PR间期代表心房开始除极的时间, QRS波群代表心室肌除极的电位变化, ST段代表心室缓慢复极的过程, T波代表心室快速复极的电位变化. 结合医学知识和专家意见, 为了能通过心电图准确判断出心肌梗死或心律失常等疾病, 我们需要准确提取的特征有P波, QRS波群, T波的振幅, 宽度, 斜率等, 以及PQ间期, PR间期, ST段等. 对于心肌梗死, 我们着重提取Q波振幅和宽度, S-T段抬高或压低幅度, 并将其转换为正弦值, T波是否倒置.

1.5 心肌梗死特征

心肌梗死又叫心肌梗塞, 是冠状动脉闭塞, 血流中断, 使部分心肌因严重的持久性缺血而发生局部坏死. 在ECG上的表现大概可以分两类: 急性心肌梗死和非典型心肌梗死.

1.5.1 急性心肌梗死在ECG的表现

1) 异常Q波(Q波的深度大于等于R波高度的四分之一).

2) S-T段抬高或压低或T波倒置. (ST段抬高超过0.1 mv或压低超过0.05 mv)^[7].

1.5.2 非典型心肌梗死在ECG的表现

1) II III avF导联QRS波振幅降低且QRS波起始出现顿挫或切迹^[8].

2) R波振幅对比下降明显, 出现顿挫或切迹. (主要为V1–V3导联).

3) P波增宽, 粗钝, 切迹成W或M型.

4) T波高耸, T/R>1.

1.6 神经网络构建

针对心肌梗死程度不同, 我们将其分为4类. 采用one-hot编码方式, 如表1所示.

对于非数值型特征值, 例如T波是否倒置, 若T波倒置, 该特征值置1, 否则, 置–1. 最后综合所有特征值来确定心肌梗死的程度. 初步提取的特征有: 异常Q波(Q波振幅和宽度), ST段抬高, ST段压低, T波倒置, R波振幅等. 对于急性心肌梗死, 采用BP神经网络, 其中激活函数使用logsig, 其函数原型为:

$ {\rm{logsig}}\left( n \right) = \frac{1}{{1 + {e^{ - n}}}} $

(13)

训练函数采用trainlm, 其函数为L-M算法, 它同时具有梯度法和牛顿法的优点. 当λ很小时, 步长等于牛顿法步长, 当λ很大时, 步长约等于梯度下降法的步长. 算法如下:

考虑如下信赖域模型:

$ \min \;\;\;\left| {f\left( {{x^k}} \right) + A\left( {{x^k}} \right)\left( {x - {x^k}} \right)} \right|^2 $

(14)

$ {\rm{s}}.{\rm{t}}.\;\;\;{\left| {x - {x^k}} \right|^2} \le {h_k} $

(15)

其中, h_k为信赖域半径. 该方程的解可有如下方程得到:

$ \left( {A{{\left( {{x^k}} \right)}^{\rm{T}}}A\left( {{x^k}} \right) + {\lambda _k}I} \right)z = - A{\left( {{x^k}} \right)^{\rm{T}}}f\left( {{x^k}} \right) $

(16)

$ {x^{k + 1}} = {x^k} - {\left( {A{{\left( {{x^k}} \right)}^{\rm{T}}}A\left( {{x^k}} \right) + {\lambda _k}I} \right)^{ - 1}}A{\left( {{x^k}} \right)^{\rm{T}}}f\left( {{x^k}} \right) $

(17)

如果 $\left\| {{{\left( {A{{\left( {{x^k}} \right)}^{\rm{T}}}A\left( {{x^k}} \right)} \right)}^{ - 1}}A{{\left( {{x^k}} \right)}^{\rm{T}}}f\left( {{x^k}} \right)} \right\| \le {h_k}$ , 则λ_k=0, 否则λ_k>0.

由于 $\left( {A{{\left( {{x^k}} \right)}^{\rm{T}}}A\left( {{x^k}} \right) + {\lambda _k}I} \right)$ 正定, 从而产生的方向z是梯度下降方向. 提取特征值后, 将每一行特征值对应的病症作为标签, 由于数据量较少, 多次实验后采用5层BP神经网络进行训练效果最佳.

BP神经网络结构如图9.

对于非典型性心肌梗死, 由于判断条件复杂多变, 在ECG表现的特征不一, 医学上至今也未能给出明确的指标, 为了更充分的提取特征, 采用卷积神经网络:

微积分中卷积的表达式为:

$ s\left( t \right) = \mathop \smallint \nolimits x\left( {t - a} \right)w\left( a \right)da $

(18)

离散形式:

$ s\left( t \right) = \mathop \sum \nolimits_a x\left( {t - a} \right)w\left( a \right) $

(19)

用矩阵表示为(*表示卷积运算):

$ s\left( t \right) = \left( {X*W} \right)\left( t \right) $

(20)

如果是二维的卷积, 则表示式为:

$ s\left( {i,j} \right) = \mathop \sum \nolimits_m \mathop \sum \nolimits_n x\left( {i - m,j - n} \right)w\left( {m,n} \right) $

(21)

其中, W为卷积核, X为输入. 如果X是一个二维输入的矩阵, 而W也是一个二维的矩阵. 如果X是多维张量, 那么W也是一个多维的张量^[9].

本文卷积神经网络部分结构如图10所示^[10,11].

将波形按10 s一次采样, 经过平滑, 滤波, 去基线等前期处理后, 送入卷积神经网络^[12], 分别经过卷积层, dropout层, 池化层, BN层, 如此循环, 最后经过Flatten层, 全连接层输出结果, 本文采用1-D卷积网络^[13].

2 实验结果

为了最小化模型结构风险, 本文采用K折交叉验证的方式, 随机地将数据切分为K个大小相同但互不相交的子集, 然后利用K–1个子集训练数据, 剩余的子集测试模型, 如此重复选择, 选出K次测试中误差最小的模型^[14]. 为了防止过拟合或欠拟合, 在训练过程中加入L1, L2混合使用的正则化. 并且本文采用训练误差和预测误差随着训练次数变化的动态曲线, 当训练误差和预测误差都在降低时, 说明还处于欠拟合, 如果预测误差开始上升, 说明处于过拟合阶段, 应该停止训练. 图11为误差随训练次数变化的曲线图.

由图可以看到, 当训练次数为700次左右时, 虽然训练误差还在降低, 但预测误差已有上升趋势, 此时网络泛化误差最小, 网络训练效果最好.

本文分别使用MIT-BIH和医大一院数据进行预测. 其中, MIT-BIH数据中没用明确区分心肌梗死类型, 所以我们只预测两类. 医大一院根据医生诊断结果分为四类. 网络预测的统计结果如表2、表3所示.

由预测结果可以看出, 对于急性心梗和正常人群预测效果理想, 因为其特征差别较大, 易于区分. 但对于非典型性心肌梗死, 由于其特征不明显, 变化较多, 所以正确率有所降低.

测试相同PTB数据集, 与其他方法的效果对比见表4.

由表3可以看出, 本文算法(CNN&BP)相对于其他方法, 正确率有所提高,主要原因是在非典型性心肌梗死的检测上, 明显优于其他方法.

为了更好的评估模型, 本文采用精准率和召回率, F值衡量模型预测结果, 并做出RoC曲线图, 也有采用MAE评价模型. 定义如下:

TP: 预测为正, 实际为正;

FP: 预测为正, 实际为负;

TN: 预测为负, 实际为负;

FN: 预测为负, 实际为正.

精准率P为:

$ P = \frac{{{{TP}}}}{{{{TP}} + {{FP}}}} $

(22)

召回率R为:

$ R = \frac{{{{TP}}}}{{{{TP}} + {{FN}}}} $

(23)

F值为:

$ F = \frac{{2*{{TP}}}}{{2*{{TP}} + {{FP}} + {{FN}}}} $

(24)

由上述公式可知, 精准率描述预测为真, 并且预测结果对的概率. 召回率描述预测为真占所有真集的概率. 而F值是精准率和召回率的调和平均, 只有当二者都高时, F才会大. 以TPR为y轴, 以FPR为x轴, 我们就直接得到了RoC曲线,TPR和FPR公式如下:

$ {{TPR}} = \frac{{{{TP}}}}{{{{TP}} + {{FN}}}} $

(25)

$ {{FPR}} = \frac{{{{FP}}}}{{{{FP}} + {{TN}}}} $

(26)

从FPR和TPR的定义可以理解, TPR越高, FPR越小, 模型和算法就越高效, 或者说在RoC曲线下所围面积越大, 算法效果越好, 图12为几种算法的RoC曲线:

由上图RoC曲线可知, CNN&BP算法在精准率和召回率上都优于其他算法.

3 总结

本文首先对ECG信号做滤波,去基线等处理, 根据波形特点在准确识别QRS波群后, 提取相应特征. 针对心肌梗死不同的类型, 采用BP神经网络和卷积神经网络^[15]相结合的方法, 并对12导联做分析. 实验结果表明, 针对不同心肌梗死类型采用不同神经网络预测在准确率上提高2%, 但由于非典型心肌梗死特征多样且不明显, 所以对此仍有待改进和完善, 文献[16]提出可以根据12导联确定心肌梗死具体起源位置. 虽然心电图自动分析程序已经达到了较高的准确率, 但是由于疾病种类繁多, 临床表现不一, 给诊断带来一定困难, 仍不能完全取代医生, 只是起到一个辅助诊断的作用. 对于一些复杂情况, 仍需要结合临床诊断.