一种心电信号心率变异性特征提取方法

1.本发明涉及心电信号分析领域，具体地说是一种心电信号心率变异性特征提取方法。


背景技术：

2.时间序列分析是统计学的分支之一，分析方法主要是时域法、频域法和非线性分析法。为了克服以往对原始时间序列单尺度分析的局限性，新的时间序列分析方法应运而生，多尺度时间序列算法被提出。时间序列的多尺度是通过粗粒化过程来完成的，其主要思想是将原时间序列通过粗粒化求出多个时间尺度的时间序列，从而对其潜在的动力学机制进行更加深入的探究。如今，对时间序列的多尺度分析得到了广泛的使用，且已被证明可以很好的将单变量时间序列转化为多尺度时间序列并提取特征，然而，如何提取一个随尺度变化的特征这一问题还有待解决。
3.近年来，随着人们物质生活水平的提高及人口老龄化进程的加快，全球心血管疾病发病率呈逐年上升趋势。有研究表明，心血管疾病已成为世界范围内死亡的首因，全球每年有1300万余人死于心血管疾病。因此，精准的诊断出心脏疾病信号显得尤为重要。研究发现，人的心脏系统是一个极其复杂的非线性系统，心电信号(ecg)是记录心电图机采集的心脏在每一个心动周期产生电活动的变化图，不同生理和病理条件下的人群，心脏各部位形体的形态不同，ecg信号自然也存在很大的差别。此外，某些心脏疾病信号在病发前与正常人的心电信号无明显差别，单从心电图上很难准确区分。


技术实现要素：

4.本发明的目的就是提供一种心电信号心率变异性特征提取方法，该方法可以较为准确地对心电信号心率变异性特征进行提取，实现对心电信号更全面、准确的分析。
5.本发明是这样实现的：一种心电信号心率变异性特征提取方法，包括如下步骤：
6.s1、获取原始心电信号，采用pan-tompkins方法对原始心电信号进行r特征波检测，提取心率变异性信号rrm，1≤m≤k，k为正整数；
7.s2、在心率变异性信号rrm中顺序选取n个采样点构成一组hrv时间序列片段{rri：1≤i≤n}；n≤k；
8.s3、对步骤s2中所形成的hrv时间序列片段{rri：1≤i≤n}进行50尺度庞加莱散点图分析，得到hrv时间序列片段在50个尺度上的庞加莱散点分布图，并提取每一尺度的特征参数sd1s，1≤s≤50；
9.s4、提取hrv时间序列片段的阈值特征w，阈值特征w的计算公式为：
10.w＝z－y111.其中，z为第s
x
尺度下的sd1值，y1为阈值；
12.阈值y1的计算公式为：
[0013][0014]
其中，h和d分别表示第s
x
尺度下健康个体和心脏疾病个体sd1的均值；
[0015]sx
的求解过程如下：
[0016]
①
选取若干健康个体和若干心脏疾病个体的原始心电信号，并依据步骤s1获得心率变异性信号，依据步骤s2构建hrv时间序列片段数据集，依据步骤s3对hrv时间序列片段数据集中的每一组hrv时间序列片段均在50个尺度上提取特征参数sd1s；
[0017]
②
分别绘制健康个体和心脏疾病个体hrv时间序列片段各尺度下sd1的误差棒图；
[0018]
③
对健康个体和心脏疾病个体的误差棒图进行多尺度分析，寻找使得健康个体和心脏疾病个体误差棒图重合部分最小的尺度因子，该尺度因子即为s
x
；
[0019]
上述h和d即可在求解s
x
的过程中获得；
[0020]
s5：若阈值特征w大于0，则判定该数据为健康个体对应的数据，若小于0，则该数据为心脏疾病个体对应的数据。
[0021]
本发明还提供了另一种并列的技术方案，该并列的技术方案中的前三步骤与上述s1-s3均相同，s4和s5如下：
[0022]
s4、选取前十个尺度的sd1计算hrv时间序列片段的尺度变化特征sv，尺度变化特征sv的计算公式为：
[0023][0024]
其中，n＝9；
[0025]
s5：若尺度变化特征sv大于0，则判定该数据为健康个体对应的数据，若小于0，则该数据为心脏疾病个体对应的数据。
[0026]
在第一个技术方案中，优选的，尺度因子s
x
等于39。
[0027]
上述两个技术方案中，优选的，步骤s2中n为3000。
[0028]
上述两个技术方案中，优选的，步骤s3具体如下：
[0029]
首先，对hrv时间序列片段进行粗粒化计算，各尺度的粗粒化时间序列计算公式为：
[0030][0031]
其中，s为尺度因子，1≤s≤50；则每个尺度下的时间序列长度等于原始时间序列长度n除以尺度因子s，即ns＝n/s；则第s尺度下的时间序列为
[0032]
接着，依次求出1-50个尺度下时间序列的庞加莱散点图的短轴sd1，即{sd11，sd12，
…
，sd1
50
}；
[0033]
各尺度下的庞加莱散点图的短轴sd1计算公式为：
[0034]
[0035]
其中，其中，为δrr的均值，ns为第s尺度时间序列的长度，即ns＝n/s。
[0036]
上述两个技术方案中，优选的，步骤s1具体如下：
[0037]
①
采用级联低通和高通滤波器组成的带通滤波器去除原始心电信号中的噪声；
[0038]
②
滤波后，采用具有传递函数的五点微分函数，放大qrs波的斜率信息，以便区分qrs波群；
[0039]
③
对信号逐点平方；
[0040]
④
对信号进行移动窗口积分；
[0041]
⑤
调整阈值，检测出qrs波群的位置，并得到r波。
[0042]
本发明提供了一种心电信号心率变异性特征提取方法，通过多尺度庞加莱散点图分析建立了hrv时间序列片段的多尺度特征，通过研究不同尺度下的特征变化特性，建立了表征心电信号心率变异性的有效特征：阈值特征和尺度变化特征，提高了特征的有效性，实现了心电信号更加全面的表示和更系统的分析。
附图说明
[0043]
图1是本发明心电信号心率变异性特征提取方法的整体流程图。
[0044]
图2是尺度粗粒化原理示意图。
[0045]
图3是本发明实施例中健康个体和心脏疾病个体50个尺度的误差棒图。
[0046]
图4是本发明实施例中尺度变化特征sv和阈值特征w的特征分布图。
具体实施方式
[0047]
为使本发明的目的、技术方案和优点更加清楚，下面结合实例和附图对本发明做进一步详细说明。
[0048]
如图1所示，本发明提出的心电信号心率变异性特征提取方法的步骤如下：
[0049]
s1、获取原始心电信号，采用pan-tompkins方法对原始心电信号进行r特征波检测，提取心率变异性(hrv)信号{rrm:1≤m≤k}，k为正整数。
[0050]
原始健康个体和心脏疾病个体的心电信号分别来源于physionet(https://www.physionet.org/)中的normal sinus rhythm database和sudden cardiac death holterdatabase，其中包含健康个体18例，心脏疾病个体23例，心脏疾病个体中只有20例用于分析。
[0051]
利用pan-tompkins方法定位r波并提取hrv信号{rrm:1≤m≤k}，具体步骤为：
[0052]
①
采用级联低通和高通滤波器组成的带通滤波器去除原始心电信号中的噪声。
[0053]
②
滤波后，采用具有传递函数的五点微分函数，放大qrs波的斜率信息，以便区分qrs波群。
[0054]
③
对信号逐点平方。
[0055]
④
对信号进行移动窗口积分。
[0056]
⑤
调整阈值，检测出qrs波群的位置，并得到r波。
[0057]
s2、在hrv信号中顺序选取3000个采样点构成一组hrv时间序列片段{rri：1≤i≤3000}，按照同样的选取方法，构建hrv时间序列片段数据集。
[0058]
通过分析上面提到的数据库，共构建健康个体hrv时间序列片段566组，心脏疾病个体hrv时间序列片段255组。
[0059]
s3、对每组hrv时间序列片段{rri：1≤i≤3000}进行50尺度庞加莱散点图分析，得到hrv时间序列片段在50个尺度上的庞加莱散点分布图，并提取每一尺度的特征参数{sd1s：1≤s≤50}，具体步骤如下：
[0060]
首先求出每组hrv时间序列片段各尺度下的时间序列，此过程是通过时间序列的粗粒化完成的，粗粒化原理示意图如图2所示，其中，
·
代表的是原始时间序列的各值，yj是粗粒化后s尺度的时间序列各值。则s尺度下的时间序列各值的计算公式为：
[0061][0062]
其中，n是原始hrv序列长度，本实施例中n为3000，s为尺度因子，则每个尺度下的时间序列长度等于原始时间序列长度n除以尺度因子s，即ns＝n/s。则第s尺度的时间序列为
[0063]
依次求出1-50个尺度下时间序列的庞加莱散点图的短轴sd1，即{sd11，sd12，
…
，sd1
50
}。
[0064]
各尺度下的庞加莱散点图的短轴sd1计算公式为：
[0065][0066]
其中，δrri为δrr中第i个元素(1≤i≤ns)，为δrr均值，ns为第s尺度时间序列的长度，即ns＝n/s，n＝3000。
[0067]
s4、提取hrv时间序列片段的第一特征：阈值特征w。具体步骤为：
[0068]
首先，分别绘制健康个体和心脏疾病个体各尺度下sd1的误差棒图。误差棒图是一种统计学方法，由均值和标准差组成，如图3所示，图中
□
和
○
分别代表健康个体和心脏疾病个体各尺度下sd1的均值，在均值上下分布的两条横线代表标准差，误差棒图的绘制是通过计算机算法完成的。
[0069]
然后对健康个体和心脏疾病个体的误差棒图(如图3所示)进行多尺度分析寻找使得健康个体和心脏疾病个体误差棒图重合部分最小的尺度因子s
x
，经计算可知，在第39尺度下重合部分最小，计算过程是由计算机算法完成。
[0070]
在尺度s
x
下设立阈值y1，阈值y1的计算公式为：
[0071][0072]
其中，h和d分别表示第s
x
尺度下健康个体和心脏疾病个体sd1的均值；
[0073]
则阈值特征w的计算公式为：
[0074]
w＝z－y1ꢀꢀꢀ
(4)
[0075]
其中，z为s
x
尺度的sd1值，y1为阈值。
[0076]
若公式(4)计算出阈值特征w的值大于0，则判定该数据为健康数据，若小于0，则该数据为心脏疾病数据。
[0077]
s5、提取hrv时间序列片段的第二特征：尺度变化特征sv。
[0078]
选取前十个尺度的sd1计算尺度变化特征sv；
[0079]
尺度变化特征sv的计算公式为：
[0080][0081]
其中，n＝9，sd1s为第s尺度下的sd1值。
[0082]
若公式(5)计算出尺度变化特征sv的值大于0，则判定该数据为健康数据，若小于0，则该数据为心脏疾病数据。
[0083]
s6、特征有效性分析。
[0084]
由步骤2可知，所有健康个体和心脏疾病个体共构建数据821组，其中，健康数据566组，心脏疾病数据255组。表1和表2分别是采用阈值特征w和尺度变化特征sv识别健康和心脏疾病数据的结果。
[0085]
表1采用阈值特征w识别健康个体和心脏疾病个体的结果
[0086][0087]
表2采用尺度变化特征sv识别健康个体和心脏疾病个体的结果
[0088][0089]
由表1和表2可知，健康和心脏疾病数据共821组，采用阈值特征w和尺度变化特征sv分别准确识别出健康和心脏疾病数据639组和709组，能够正确识别出大部分数据，识别准确率高，能够说明本发明提出的特征的有效性。
[0090]
图4为阈值特征w和尺度变化特征sv的特征分布图，其中，横坐标是尺度变化特征sv的值，纵坐标是阈值特征w的值，由图中可知心脏疾病数据大致分布在图的左下方，健康数据在右上方，区分明显，能够进一步说明本发明提出的特征的有效性。