一种基于多尺度熵的心脏疾病检测方法与流程

本发明涉及到生物医学工程领域，尤其涉及到一种基于多尺度熵的心脏疾病检测方法。

背景技术：

在生活节奏日益加快的今天，心血管疾病一直是一个影响人们健康的重要因素。根据有关组织的调查，心血管疾病的死亡率多年来均超过癌症，居所有疾病之首。其中在我国，心血管疾病致死占农村居民总致死原因的44.6％，占城市居民总致死原因的42.52％。有报告称在我国患有心血管疾病的国民已接近3亿人。随着我国人口老龄化的加剧、城镇化进程的加快和人民生活水平的改变(如不合理饮食造成的肥胖或高血压)，我国心血管疾病患者的人数正在并仍将在未来若干年里快速增长。心血管疾病一旦发作，造成的后果十分严重。如果能够尽早发现、尽早确诊、尽早治疗，能够大大降低心血管疾病死亡率。

心电图记录了心脏在跳动过程中心脏电信号变化的曲线，是心脏电生理活动在人体体表的综合表现。心电图包含了丰富的心脏功能和病理信息，能够直观准确地反映心脏的电活动特性和表现心脏的工作状态，是心脏从兴奋发生到传播在到恢复整个过程的客观评价标准。通过分析患者心电图来评估患者的心脏功能和心血管疾病是一种非常方便和实用方法，在日常的诊断中这种方法也发挥了无可替代的作用。但仅仅通过人工分析心电信号来诊断心脏疾病，不仅加大了医生的负担，对患者来说也并不方便。而且人工分析有场所和时间的限制，但心血管疾病的发生是随时随地的，所以需要一种能应用在心电检测设备上的自动检测心脏疾病的方法。

技术实现要素：

为了解决背景技术种提到的问题，本发明提出了一种基于多尺度熵的心脏疾病检测方法。

本发明的技术方案是一种基于多尺度熵的心脏疾病检测方法，其特征在于，包括以下步骤：

步骤1：首先将原始心电信号放入带通滤波器滤去部分噪声，再使用信号微分、平方手段放大得到r波特征被放大的信号，将r波特征被放大的信号用动态阈值调整方法标记r波位置，获得心电信号的rr间期序列；

步骤2：根据心电信号的rr间期序列进行经验模态分解，为了抑制端点效应，需要先对信号进行延拓，再通过构建信号的上下包络线来对心电信号进行分解得到imf分量，根据imf分量得到健康人和心脏疾病患者心电信号的本征函数信号；

步骤3：通imf分量计算健康人和心脏疾病患者心电信号的本征函数信号的多尺度熵，利用支持向量机的分类功能对健康人和心脏疾病患者心电信号的本征函数信号进行分类，区分出正常人和心脏疾病患者的心电信号；

作为优选，步骤1中所述原始心电信号为zi(i∈1,2,...,n)，n为信号采样点的数量；

带通滤波由低通滤波器和高通滤波器叠加组成，滤波的频率范围为[fh-fl]；

其中，低通滤波器的传递函数为h(z)l,低通滤波器的截止频率为fl,即通过低通滤波器将频率高于fl的信号滤除，低通滤波器的增益为al，低通滤波器的过滤处理延迟dl个单位；

低通滤波器的传递函数为:

高通滤波器的传递函数为h(z)h,高通滤波器的截止频率为fh,即通过高通滤波器将频率低于fh的信号滤除，低通滤波器的增益为ah，低通滤波器的过滤处理延迟dh个单位；

高通滤波器的传递函数为:

该滤波器具有整数系数，在实际工作中可以极大减少处理器的工作负担；采用低通和高通滤波器结合的方式，配合了心电信号的采样频率；

原始信号zi经过低通函数h(zi)l和高通滤波器h(zi)h后，得到了滤去了噪声的心电信号zl,i(i∈1,2,...,n)，将zl,i(i∈1,2,...,n)微分，平方和移动窗口积分后得到r波特征被放大的信号z2,i(i∈1,2,...,n)，将z2,i通过动态阈值调整法标记其r波峰；

步骤1所述的将r波特征被放大的信号用动态阈值调整方法标记r波位置为：

整个检测方法有两套阈值，采用对半检测，spki是qrs复波的波峰幅值，npki是非qrs波的波峰幅值，thresholdi1是应用的第一组阈值，thresholdi2是应用的第二组阈值。spki，npki，thresholdi1和thresholdi2在最开始有一个运行估计值，如果所检出的峰值大于thresholdi1，该峰值被标记为信号峰值spki，否则该峰值被标记为噪声信号npki，阈值大小由下式更新：

spki＝0.125*peaki+0.875*spki(如果peaki是信号峰值)

npki＝0.125*peaki+0.875*npki(如果peaki是噪声峰值)

thresholdi1＝npki+0.25(spki-npki)

thresholdi2＝0.5thresholdi1

式中peak表示总体峰值。

采用双阈值技术来检测r波波峰位置，thresholdi1应用于滤波后的信号z1,i，thresholdi2应用于移动窗口积分产生的信号z2,i，当波峰在z1,i和z2,i都被标记为信号峰值时，它被确认为r波；假设检测到m个r波，得到r波序列rj(j∈1,2,...,m)，每两个相邻r波波峰位置相减得到序列的一个rr间期，rj所有的rr间期构成了信号的rr间期序列xi(i∈1,2,...,n)；

作为优选，步骤2所述心电信号的rr间期序列进行经验模态分解为：

首先要找到信号所有的极大值和极小值，对信号进行微分，得到微分后的信号tn(n∈1,2,...,n)。

tn＝xi-xi-1,i∈2,3,...,n

然后将信号分为两列，第一列d1＝tn1，n1∈(1，2，...，n-1)，第二列d2＝tn2，n2∈(2，3，...，n)；创建函数

得到新序列c1＝f1(d1*d2)，c2＝f1(d1)，c3＝f2(d1)；从而得到序列中的最大值和最小值；

indminj＝{x|x≠0,x∈c1&c2,j∈1,2,...,n}

indmaxj＝{x|x≠0,x∈c1&c3,j∈1,2,...,n}

步骤2中所述的用极大值和极小值点拟合包络线的方法使用三次样条插值法，极大值和极小值来自indminj和indmaxj，三次样条插值就是使用三次多项式去描述两个相邻点间的线段，三次样条插值后形成了信号的上包络线e+(d)和下包络线e-(d)；

将来自indminj和indmaxj的极大值和极小值放到平面直角坐标系中，构成坐标(d0，b0),……(dn，bn)，dn(n∈1,2,...,n,)表示该极大值或极小值在xi中的排序，是第几个点，bn表示对应点的值，步长可以表示为

hi＝di+1-di(i＝0,1,…,n-1)

根据数据节点和首位的端点条件可以得到矩阵方程

解矩阵方程，得到样条曲线的系数为

ai＝bi

在每个子区间xi≤x≤xi+1中，创建方程

gi(x)＝ai+gi(d-di)+li(d-di)²+fi(d-di)³

n+1个极大值或极小值构成了n个子区间，每个子区间的曲线方程联合起来构成了信号上包络线和下包络线；

步骤2中所述的分解步骤如下，利用信号的上包络线e+(d)和下包络线e-(d)计算均值包络线：

将原信号减去m1(d)得到一个去低频的新信号

h¹1(d)＝x(d)-m1(d)

此时信号h¹1(d)一般不满足在整个信号上极值点的个数和过零点的个数相差不大于且在任意点处上下包络的均值为零这两个条件，则重复上述步骤，假定k次后满足；则原信号x(d)的一阶本征模态函数分量为：

用原信号x(d)减去c1(d)，得到一个去高频成分的新信号r1(d)，有：

r1(d)＝x(d)-c1(d)

对r1(d)重复得到c1(d)的过程，得到第二个imf分量c2(d)，如此反复进行，直到第n阶imf分量cn(d)或其余量rn(d)小于预设值；或当残余分量rn(d)是单调函数或常量时，emd分解过程停止。c1(d)，c2(d)，…,cn(d)即为步骤2中所述最后得到健康人和心脏疾病患者心电信号的本征函数信号；

作为优选，步骤3中所述计算分解后信号的多尺度熵为：

将步骤2得到的imf分量cn(d)记为imp(p∈1,2,...,n)，先进性尺度τ(τ∈1,2,...,n)的粗粒化，构造出新向量m^*表示比较向量得长度，一般取常数1或2。

阈值r

r＝0.15*std({y^(τ)})

std表示计算序列粗粒化序列的标准差；

步骤3中由imp和im^*p重构得到的两个相量yq和y^*q，有定义

max表示取所有元素中的最大值,d表示两矢量之间的距离，由对应元素的最大差值决定；

统计满足以下条件的向量个数：

求对所有q值得平均值，即

令m^*＝m^*+1，重复上述步骤得到其中

在尺度τ下多尺度熵值mse(τ)如下

作为优选，步骤3中所述支持向量机：

首先将rr间期数据通过步骤2和步骤3处理后得到的多尺度熵添加标记位δ，健康人多尺度熵添加标记为δ＝0记为msen，心脏疾病患者多尺度熵添加标记位δ＝1记为msei，取msen和msei作为svm的训练集合记为u；

svm就是通过选择训练集合u里的数据，构建一个超平面保证分类精度的同时，能够使超平面两侧的空白区域也就是分类间隔最大化，设该超平面的方程为

(w^t·x)+b^*＝0

可以计算出其分类间隔为

所以整个分类机的目标函数可以写为：

约束条件为：

δi(w^tui+b)≥1,i＝1,...,n

引入lagrange函数：

其中a＞0为lagrange乘数，约束最优化问题的解由lagrange函数的鞍点决定，并且最优化问题的解在鞍点处满足对w和b的偏导为0，将该二次规划问题转化为相应的对偶问题：

解得最优解a^*＝(a^*1,a^*2,...,a^*l)^t；

计算最优权值向量w^*和最优偏置b^*，分别为：

式中，下标因此得到最优分类超平面(w^t·x)+b^*＝0；

判定条件为：

当信号标记位δ＝0时，信号为健康心电信号。当信号标记位δ＝1时，信号为心脏疾病患者心电信号。

本发明的优点为：

使用emd分解对心电信号进行处理和使用多尺度熵来衡量心脏的健康程度。

emd分解不需要事先选定基函数，而是从本身的尺度特征出发自适应的产生合适的模态函数imf，这些imf分量从高频到低频逐次分布，能够很好的反映出信号在任何时间局部的频率特性。

多尺度熵是一种衡量时间序列复杂性的方法，已有研究表明，相关生理信号复杂度出现降低是病理动力学的一个普遍特征，但具体应用到心电信号的分析中还是少数。

本发明通过emd分解和多尺度熵分析不仅可以及时检测出心脏的健康情况，还可以帮助研究人员在某种程度上了解疾病的原理。

附图说明

图1：本发明方法流程图；

图2：pan-tompkins算法流程图；

图3：emd算法流程图；

图4：多尺度熵算法流程图。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合附图及实施例对本发明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不用于限定本发明。

本发明具体实施方式中原始心电信号来自mit-bih心电数据库。

下面结合图1至图4论述本发明实施例。本发明实施例的具体步骤包括以下步骤：

步骤1：首先将原始心电信号放入带通滤波器滤去部分噪声，再使用信号微分、平方手段放大得到r波特征被放大的信号，将r波特征被放大的信号用动态阈值调整方法标记r波位置，获得心电信号的rr间期序列；

步骤1中所述原始心电信号为zi(i∈1,2,...,n)，n为信号采样点的数量；

带通滤波由低通滤波器和高通滤波器叠加组成，滤波的频率范围为[fh-fl]，一般地，可以取fh＝5hz，fl＝15hz；

其中，低通滤波器的传递函数为h(z)l,低通滤波器的截止频率为fl,即通过低通滤波器将频率高于fl的信号滤除，低通滤波器的增益为al＝36db，低通滤波器的过滤处理延迟dl＝6samples；

低通滤波器的传递函数为:

高通滤波器的传递函数为h(z)h,高通滤波器的截止频率为fh,即通过高通滤波器将频率低于fh的信号滤除，高通滤波器的增益为al＝32db，高通滤波器的过滤处理延迟dh＝16samples；

高通滤波器的传递函数为:

该滤波器具有整数系数，在实际工作中可以极大减少处理器的工作负担；采用低通和高通滤波器结合的方式，配合了心电信号的采样频率；

原始信号zi经过低通函数h(zi)l和高通滤波器h(zi)h后，得到了滤去了噪声的心电信号zl,i(i∈1,2,...,n)，将zl,i(i∈1,2,…,n)微分，平方和移动窗口积分后得到r波特征被放大的信号z2,i(i∈1,2,...,n)，将z2,i通过动态阈值调整法标记其r波峰；

步骤1所述的将r波特征被放大的信号用动态阈值调整方法标记r波位置为：

整个检测方法有两套阈值，采用对半检测，spki是qrs复波的波峰幅值，npki是非qrs波的波峰幅值，thresholdi1是应用的第一组阈值，thresholdi2是应用的第二组阈值。spki，npki，thresholdi1和thresholdi2在最开始有一个运行估计值，如果所检出的峰值大于thresholdi1，该峰值被标记为信号峰值spki，否则该峰值被标记为噪声信号npki，阈值大小由下式更新：

spki＝0.125*peaki+0.875*spki(如果peaki是信号峰值)

npki＝0.125*peaki+0.875*npki(如果peaki是噪声峰值)

thresholdi1＝npki+0.25(spki-npki)

thresholdi2＝0.5thresholdi1

式中peak表示总体峰值。

采用双阈值技术来检测r波波峰位置，thresholdi1应用于滤波后的信号z1,i，thresholdi2应用于移动窗口积分产生的信号z2,i，当波峰在z1,i和z2,i都被标记为信号峰值时，它被确认为r波；假设检测到m个r波，得到r波序列rj(j∈1,2,...,m)，每两个相邻r波波峰位置相减得到序列的一个rr间期，rj所有的rr间期构成了信号的rr间期序列xi(i∈1,2,...,n)；

步骤2：根据心电信号的rr间期序列进行经验模态分解，为了抑制端点效应，需要先对信号进行延拓，再通过构建信号的上下包络线来对心电信号进行分解得到imf分量，根据imf分量得到健康人和心脏疾病患者心电信号的本征函数信号；

步骤2所述心电信号的rr间期序列进行经验模态分解为：

首先要找到信号所有的极大值和极小值，对信号进行微分，得到微分后的信号tn(n∈1,2,...,n)。

tn＝xi-xi-1,i∈2,3,...,n

然后将信号分为两列，第一列d1＝tn1，n1∈(1，2，...，n-1)，第二列d2＝tn2，n2∈(2，3，...，n)；创建函数

得到新序列c1＝f1(d1*d2)，c2＝f1(d1)，c3＝f2(d1)；从而得到序列中的最大值和最小值；

indminj＝{x|x≠0,x∈c1&c2,j∈1,2,...,n}

indmaxj＝{x|x≠0,x∈c1&c3,j∈1,2,...,n}

步骤2中所述的用极大值和极小值点拟合包络线的方法使用三次样条插值法，极大值和极小值来自indminj和indmaxj，三次样条插值就是使用三次多项式去描述两个相邻点间的线段，三次样条插值后形成了信号的上包络线e+(d)和下包络线e-(d)；

将来自indminj和indmaxj的极大值和极小值放到平面直角坐标系中，构成坐标(d0，b0),……(dn，bn)，dn(n∈1,2,...,n,)表示该极大值或极小值在xi中的排序，是第几个点，bn表示对应点的值，步长可以表示为

hi＝di+1-di(i＝0,1,…,n-1)

根据数据节点和首位的端点条件可以得到矩阵方程

解矩阵方程，得到样条曲线的系数为

ai＝bi

在每个子区间xi≤x≤xi+1中，创建方程

gi(x)＝ai+gi(d-di)+li(d-di)²+fi(d-di)³

n+1个极大值或极小值构成了n个子区间，每个子区间的曲线方程联合起来构成了信号上包络线和下包络线；

步骤2中所述的分解步骤如下，利用信号的上包络线e+(d)和下包络线e-(d)计算均值包络线：

将原信号减去m1(d)得到一个去低频的新信号

h¹1(d)＝x(d)-m1(d)

此时信号h¹1(d)一般不满足在整个信号上极值点的个数和过零点的个数相差不大于且在任意点处上下包络的均值为零这两个条件，则重复上述步骤，假定k次后满足；则原信号x(d)的一阶本征模态函数分量为：

用原信号x(d)减去c1(d)，得到一个去高频成分的新信号r1(d)，有：

r1(d)＝x(d)-c1(d)

对r1(d)重复得到c1(d)的过程，得到第二个imf分量c2(d)，如此反复进行，直到第n阶imf分量cn(d)或其余量rn(d)小于预设值；或当残余分量rn(d)是单调函数或常量时，emd分解过程停止。c1(d)，c2(d)，…,cn(d)即为步骤2中所述最后得到健康人和心脏疾病患者心电信号的本征函数信号；

步骤3：通imf分量计算健康人和心脏疾病患者心电信号的本征函数信号的多尺度熵，利用支持向量机的分类功能对健康人和心脏疾病患者心电信号的本征函数信号进行分类，区分出正常人和心脏疾病患者的心电信号；

步骤3中所述计算分解后信号的多尺度熵为：

将步骤2得到的imf分量cn(d)记为imp(p∈1,2,...,n)，先进性尺度τ(τ∈1,2,...,n)的粗粒化，构造出新向量m*表示比较向量得长度，一般取常数1或2；

阈值r

r＝0.15*std({y^(τ)})

std表示计算序列粗粒化序列的标准差；

步骤3中由imp和im^*p重构得到的两个相量yq和y^*q，有定义

max表示取所有元素中的最大值,d表示两矢量之间的距离，由对应元素的最大差值决定；

统计满足以下条件的向量个数：

求对所有q值得平均值，即

令m^*＝m^*+1，重复上述步骤得到其中

在尺度τ下多尺度熵值mse(τ)如下

步骤3中所述支持向量机：

首先将rr间期数据通过步骤2和步骤3处理后得到的多尺度熵添加标记位δ，健康人多尺度熵添加标记为δ＝0记为msen，心脏疾病患者多尺度熵添加标记位δ＝1记为msei，取msen和msei作为svm的训练集合记为u；

svm就是通过选择训练集合u里的数据，构建一个超平面保证分类精度的同时，能够使超平面两侧的空白区域也就是分类间隔最大化，设该超平面的方程为

(w^t·x)+b^*＝0

可以计算出其分类间隔为

所以整个分类机的目标函数可以写为：

约束条件为：

δi(w^tui+b)≥1,i＝1,...,n

引入lagrange函数：

其中a＞0为lagrange乘数，约束最优化问题的解由lagrange函数的鞍点决定，并且最优化问题的解在鞍点处满足对w和b的偏导为0，将该二次规划问题转化为相应的对偶问题：

解得最优解a^*＝(a^*1,a^*2,...,a^*l)^t；

计算最优权值向量w^*和最优偏置b^*，分别为：

式中，下标因此得到最优分类超平面(w^t·x)+b^*＝0；

判定条件为：

当信号标记位δ＝0时，信号为健康心电信号。当信号标记位δ＝1时，信号为心脏疾病患者心电信号。

应当理解的是，本说明书未详细阐述的部分均属于现有技术。

应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。