一种基于极值能量分解法的心率变异性信号分析方法
技术领域
本发明涉及一种心电图信号分析,尤其涉及一种基于极值能量分解法的心率变异性信号分析方法。
背景技术
生理信号是由生命体多个***相互作用产生的,不同***作用的时间和强度不同,导致生理信号具有时间和空间上的复杂性。心率变异性(HRV)是指测量连续心动周期之间的时间变异数,准确地说,应该是测量连续出现的正常P-P间期之间的差异的变异数。但由于P波不如R波明显或P波顶端有时宽钝,所以对心率变异性信号的研究通常用与P-P间期相等的R-R间期信号(RRI)来代替。研究表明,HRV可做为植物神经***活动的无创性检测指标,尤其在判断某些心血管疾病的预后方面有着重要意义。
现有技术中研究心脏的长时调节节律(<1Hz))常用心律变异性信号(HRV信号)作为分析对象。大量的研究表明,人体HRV信号具有长时相关性和非线性动力学复杂性,并且年龄和疾病会导致动力学复杂性降低。对心率变异性(Heart Rate Varibility,HRV)信号的研究常用的是RR间期(Interbeat Intervals,RRI)信号,即连续RRI信号R波之间的时间间隔信号。
研究HRV信号的能量改变最常用的方法是功率谱分析(PSD)。PSD通过傅里叶变换将HRV信号的功率转化成频率的函数,研究不同频域范围的功率大小,通常HRV频谱分为高频(HF)、低频(LF)和极低频(VLF)等几个频段。LF/HF比值有重要的临床价值。心脏疾病会引起HRV功率谱的改变,比如心衰和心肌梗死引起归一化HF增高、LF和VLF降低。然而,PSD不是一种基于数据驱动的方法,并且对频域的分段比较粗糙,导致细节缺失,分割也不够灵活。
因此,亟待解决上述问题。
发明内容
发明目的:本发明的目的是提供一种可采用较少数据即可直观反映心电图能量分布和信号波动的真实规律的基于极值能量分解法的心率变异性信号分析方法。
技术方案:为实现以上目的,本发明公开了一种基于极值能量分解法的心率变异性信号分析方法,包括如下步骤:
(1)、获取给定时间和给定采样频率下的未知状态的ECG信号,然后对ECG信号进行去噪预处理,从中提取RRI信号,得到未知状态的RRI信号x(t);
(2)、将RRI信号x(t)作为原始信号,求出原始信号的所有局部极值点,然后将原始信号的所有极大值点和所有极小值点采用样条曲线连起来分别形成上包络线emax和下包络线emin,得到上、下包络线的包络均值信号m(t)=(emax+emin)/2;
(3)、将原始信号x(t)减去包络均值信号m(t),得到h(t)=x(t)-m(t);然后判断h(t)是否满足极值模态函数的判定条件,如果不满足,将h(t)作为原始信号返回至步骤(3),直到hk(t)满足极值模态函数的判定条件,则记c1(t)=hk(t),作为第一个极值模态函数分量;
(4)、将原始信号x(t)减去第一个极值模态函数分量c1(t),得到余量r1(t)=x(t)-c1(t),然后判断hk(t)是否满足停止准则,如果不满足,将r1(t)作为新的原始序列x(t),返回至步骤(2)和(3);如果hk(t)满足停止准则但n<8时,返回步骤(1)重新获取原始信号;如果hk(t)满足停止准则且n≥8时,得到第2、3、…、n个极值模态函数分量及余量rn(t),于是将原始信号x(t)分解为n个极值模态函数分量和一个余量,即
(5)、对极值模态函数分量ci(t),i=1,2,…,n,进行频谱分析得到各极值模态函数分量的中心频率;
(6)、将原始信号x(t)分解得的n个极值模态函数分量,代表了原始信号不同频段的分量,然后计算其各个分量的能量
Ei=∫|ci(t)|2dt,i=1,2,…,n
将每一个能量值归一化,得到归一化的能量分布向量
pi=Ei/E,i=1,2,...,n
其中,
第一个分量p
1表示最高频段的能量,代表了信号在最高频段范围内能量分布的比例,最后一个分量p
n表示信号在最低频段范围内能量分布的比例;根据归一化的能量分布向量绘制归一化能量分布图,其中横坐标表示分量层次,纵坐标表示归一化的能量分布向量值,曲线表示平均值,误差棒表示标准差;
(7)、选取第二个分量p2至第七个分量p7,计算能量差异值EDV,EDV=(p2+p3+p4)-(p5+p6+p7),当EDV≤第一标准值时,则判定该RRI信号为正常心率变异性信号,当第一标准值<EDV<第二标准值时,则判定该RRI信号为疑似异常心率变异性信号,当EDV≥第二标准值时,则判定该RRI信号为异常心率变异性信号。
其中,所述原始信号x(t)所需最少数据量N=2n+1,其中n为分解出的极值模态函数分量的数量。
优选的,所述步骤(1)中去噪预处理的具体方法为:将ECG信号经过40Hz零相位FIR低通滤波器滤波消除高频噪声,然后经过中值滤波器去除基线漂移。
再者,所述步骤(3)中极值模态函数的判定条件为:(a)、在整个数据序列中,极值点的数量与过零点的数量相等或者相差一个;(b)、在任意时刻,上下包络线对于时间轴对称。
进一步,所述步骤(4)中hk(t)满足停止准则公式为:
优选的,所述步骤(7)中的第一标准值为-0.15,第二标准值为0.08。
有益效果:与现有技术相比,本发明具有以下显著优点:
本发明采用极值能量分解方法(Extremum Energy Decomposition,EED)分析RRI信号,将原始信号分解为多个分量,也就是极值分量函数,计算每一个分量的能量,得到其能量分布;本发明的可依据生RRI信号自身的波动规律将信号分解为从高频到低频的不同时间层次信号,对频段的分割较为细致;极值分解在所有层次上得到的数据长度相同,因而不会导致数据长度减小,从而使其可以用于短时间数据分析,即需要很少数据量即可分析得到准确结果;EED对于不同层次分量能量分析不容易受噪声的影响。
附图说明
图1为本发明中原始信号的示意图;
图2为本发明中原始信号求取包络线的示意图;
图3为本发明中原始信号的减去包络均值信号的示意图;
图4为本发明中得到第一个极值模态函数分量的示意图;
图5为本发明中极值能量分解法的流程示意图;
图6为本发明实施例1中RRI信号的EED分解示意图;
图7为本发明实施例1中RRI信号中归一化能量分布示意图。
具体实施方式
下面结合附图对本发明的技术方案作进一步说明。
如图1、图2、图3、图4和图5所示,本发明的一种基于极值能量分解法的心率变异性信号分析方法,包括如下步骤:
(1)、获取给定时间和给定采样频率下的未知状态的ECG信号,然后对ECG信号进行去噪预处理,从中提取RRI信号,得到未知状态的RRI信号x(t);其中去噪预处理的具体方法为:由于ECG能量主要集中在0~40Hz,将ECG信号经过40Hz零相位FIR低通滤波器滤波消除高频噪声,然后经过中值滤波器去除基线漂移;
(2)、将RRI信号x(t)作为原始信号,原始信号x(t)所需最少数据量N=2n+1,其中n为分解出的极值模态函数分量的数量;求出原始信号的所有局部极值点,然后将原始信号的所有极大值点和所有极小值点采用样条曲线连起来分别形成上包络线emax和下包络线emin,得到上、下包络线的包络均值信号m(t)=(emax+emin)/2;
(3)、将原始信号x(t)减去包络均值信号m(t),得到h(t)=x(t)-m(t);然后判断h(t)是否满足极值模态函数的判定条件,如果不满足,将h(t)作为原始信号返回至步骤(2),直到hk(t)满足极值模态函数的判定条件,则记c1(t)=hk(t),作为第一个极值模态函数分量;其中极值模态函数的判定条件为:(a)、在整个数据序列中,极值点的数量与过零点的数量相等或者相差一个;(b)、在任意时刻,上下包络线对于时间轴对称;
(4)、将原始信号x(t)减去第一个极值模态函数分量c1(t),得到余量r1(t)=x(t)-c1(t),然后判断hk(t)是否满足停止准则,如果不满足,将r1(t)作为新的原始序列x(t),返回至步骤(2)和(3);如果hk(t)满足停止准则但n<8时,返回步骤(1)重新获取原始信号;如果hk(t)满足停止准则且n≥8时,得到第2、3、…、n个极值模态函数分量及余量rn(t),于是将原始信号x(t)分解为n个极值模态函数分量和一个余量,即
其中hk(t)满足停止准则公式为:
ε表示筛选门限,取0.2~0.3之间;满足停止准则的极值模态分解则满足如下两个条件:(a)最后得到的极值模态函数分量c
n(t)或者余量r
n(t)小于预先设定的阈值;(b)残余信号r
n(t)成为单调信号,不能从中再提取出极值模态函数信号;
(5)、对极值模态函数分量ci(t),i=1,2,…,n,进行频谱分析得到各极值模态函数分量的中心频率;
(6)、将原始信号x(t)分解得的n个极值模态函数分量,代表了原始信号不同频段的分量,然后计算其各个分量的能量
Ei=∫|ci(t)|2dt,i=1,2,…,n
将每一个能量值归一化,得到归一化的能量分布向量
pi=Ei/E,i=1,2,...,n
其中,
第一个分量p
1表示最高频段的能量,代表了信号在最高频段范围内能量分布的比例,最后一个分量p
n表示信号在最低频段范围内能量分布的比例;根据归一化的能量分布向量绘制归一化能量分布图,其中横坐标表示分量层次,纵坐标表示归一化的能量分布向量值,曲线表示平均值,误差棒表示标准差;
(7)、选取第二个分量p2至第七个分量p7,计算能量差异值EDV,EDV=(p2+p3+p4)-(p5+p6+p7),当EDV≤第一标准值时,则判定该RRI信号为正常心率变异性信号,当第一标准值<EDV<第二标准值时,则判定该RRI信号为疑似异常心率变异性信号,当EDV≥第二标准值时,则判定该RRI信号为异常心率变异性信号;其中的第一标准值为-0.15,第二标准值为0.08。
本发明采用的极值能量分解方法(Extremum Energy Decomposition,EED),其是基于极值模态函数的概念的方法,极值模态函数是同时满足下面两个条件的具有单一频率的一类信号,两个条件为:
(a)、在整个数据序列中,极值点(包括极大值和极小值)的数量与过零点的数量必须相等或者最多相差一个;
(b)在任意时刻,局部极大值点形成的上包络线与局部极小值点形成的下包络线的均值为零,也就是说局部上下包络线对于时间轴局部对称;
上面两个条件,条件(a)类似于高斯正态平稳过程对于传统窄带的要求,条件(b)保证了由极值模态函数计算得到的瞬时频率有物理意义。
本发明的极值模态函数分解终止的标准选择要适中,条件太严格,会导致最后几个极值模态函数分量失去意义;条件太宽松,会导致有用的分量丢失;在实际应用中,也可以根据需求设定需要分解的极值模态函数分量层数,当满足分解层数时即终止计算。
实施例1
将EED分析方法用于分析健康人和CHF患者ECG在不同层次下的能量分布。
健康人的一种基于极值能量分解法的心率变异性信号分析方法,包括如下步骤:
(1)、从physionet的RR间期数据库nsr2db获取健康人群的ECG信号;其中数据包含54个健康人(年龄28~76,平均61),然后对ECG信号进行去噪预处理,从中提取RRI信号,得到RRI信号x(t);其中去噪预处理的具体方法为:由于ECG能量主要集中在0~40Hz,将ECG信号经过40Hz零相位FIR低通滤波器滤波消除高频噪声,然后经过中值滤波器去除基线漂移;
(2)、将RRI信号x(t)作为原始信号,原始信号x(t)所需最少数据量N=2n+1=29,其中n为分解出的极值模态函数分量的数量,n=8;求出原始信号的所有局部极值点,然后将原始信号的所有极大值点和所有极小值点采用样条曲线连起来分别形成上包络线emax和下包络线emin,得到上、下包络线的包络均值信号m(t)=(emax+emin)/2;
(3)、将原始信号x(t)减去包络均值信号m(t),得到h(t)=x(t)-m(t);然后判断h(t)是否满足极值模态函数的判定条件,如果不满足,将h(t)作为原始信号返回至步骤(2),直到hk(t)满足极值模态函数的判定条件,则记c1(t)=hk(t),作为第一个极值模态函数分量;其中极值模态函数的判定条件为:(a)、在整个数据序列中,极值点的数量与过零点的数量相等或者相差一个;(b)、在任意时刻,上下包络线对于时间轴对称;
(5)、将原始信号x(t)减去第一个极值模态函数分量c1(t),得到余量r1(t)=x(t)-c1(t),将r1(t)作为新的原始序列x(t),返回至步骤(2)和(3),得到第2、3、…、8个极值模态函数分量及余量r8(t),于是将原始信号x(t)分解为8个极值模态函数分量和一个余量,即
得到如图6所示的健康人的极值模态分解示意图,可以看到分量1具有最高的频率,信号在最短的时间尺度上波动,随着分量序号增加,频率逐渐降低;
(6)、对极值模态函数分量ci(t),i=1,2,…,8,进行频谱分析得到各极值模态函数分量的中心频率,得到频域分析结果图;
(7)、将原始信号x(t)分解得的8个极值模态函数分量,代表了原始信号不同频段的分量,然后计算其各个分量的能量
Ei=∫|ci(t)|2dt,i=1,2,…,8
将每一个能量值归一化,得到归一化的能量分布向量
pi=Ei/E,i=1,2,...,8
其中,
第一个分量p
1表示最高频段的能量,代表了信号在最高频段范围内能量分布的比例,最后一个分量p
n表示信号在最低频段范围内能量分布的比例。
(8)、选取第二个分量p2至第七个分量p7,计算能量差异值EDV,EDV=(p2+p3+p4)-(p5+p6+p7)=-0.2092±0.2940。
CHF患者的一种基于极值能量分解法的心率变异性信号分析方法,包括如下步骤:
(1)、从physionet的RR间期数据库chf2db数据库获取ECG信号,其中chfdb数据库包含29个充血性心衰(Congestive Heart Failure,CHF)患者(年龄34~79,平均55),然后对ECG信号进行去噪预处理,从中提取RRI信号,得到RRI信号x(t);其中去噪预处理的具体方法为:由于ECG能量主要集中在0~40Hz,将ECG信号经过40Hz零相位FIR低通滤波器滤波消除高频噪声,然后经过中值滤波器去除基线漂移;
(2)、将RRI信号x(t)作为原始信号,原始信号x(t)所需最少数据量N=2n+1=29,其中n为分解出的极值模态函数分量的数量,n=8;求出原始信号的所有局部极值点,然后将原始信号的所有极大值点和所有极小值点采用样条曲线连起来分别形成上包络线emax和下包络线emin,得到上、下包络线的包络均值信号m(t)=(emax+emin)/2;
(3)、将原始信号x(t)减去包络均值信号m(t),得到h(t)=x(t)-m(t);然后判断h(t)是否满足极值模态函数的判定条件,如果不满足,将h(t)作为原始信号返回至步骤(2),直到hk(t)满足极值模态函数的判定条件,则记c1(t)=hk(t),作为第一个极值模态函数分量;其中极值模态函数的判定条件为:(a)、在整个数据序列中,极值点的数量与过零点的数量相等或者相差一个;(b)、在任意时刻,上下包络线对于时间轴对称;
(4)、将原始信号x(t)减去第一个极值模态函数分量c1(t),得到余量r1(t)=x(t)-c1(t),将r1(t)作为新的原始序列x(t),返回至步骤(2)和(3),得到第2、3、…、8个极值模态函数分量及余量r8(t),于是将原始信号x(t)分解为8个极值模态函数分量和一个余量,即
(5)、对极值模态函数分量ci(t),i=1,2,…,8,进行频谱分析得到各极值模态函数分量的中心频率;
(6)、将原始信号x(t)分解得的8个极值模态函数分量,代表了原始信号不同频段的分量,然后计算其各个分量的能量
Ei=∫|ci(t)|2dt,i=1,2,…,8
将每一个能量值归一化,得到归一化的能量分布向量
pi=Ei/E,i=1,2,...,8
其中,
第一个分量p
1表示最高频段的能量,代表了信号在最高频段范围内能量分布的比例,最后一个分量p
n表示信号在最低频段范围内能量分布的比例;根据健康人和CHF患者的归一化的能量分布向量绘制归一化能量分布图,其中横坐标表示分量层次,纵坐标表示归一化的能量分布向量值,曲线表示平均值,误差棒表示标准差;
(7)、选取第二个分量p2至第七个分量p7,计算能量差异值EDV,EDV=(p2+p3+p4)-(p5+p6+p7)=0.2642±0.4070。
本发明为了进一步分析了健康人和CHF患者HRV信号不同层次分量的平均中心频率,得到表1的结果。
表1健康人和CHF患者HRV信号不同分量层次的平均中心频率
可以得出随着分量层次增加HRV的中心频率逐渐降低。
传统的功率谱密度(Power Spectral Density,PSD)方法对频域分割的典型方式为:HF(0.15~0.4Hz),LF(0.04~0.15Hz),VLF(0.0033~0.04Hz)。本发明EED方法对分量层次1,频率高于HF;层次2在HF的频率范围;层次3、4在LF的范围;层次5~7在VLF的频率范围;层次8低于VLF。另外可以看到相同层次时CHF患者的频率略高于健康人,反应了心脏疾病对HRV波动节律的影响,同样层次下CHF患者HRV波动更快。
本发明将二个群体的HRV信号分解得的极值模态函数分量Ci(t),计算获得归一化能量分布向量,并作EED曲线图,如图7所示。图7为健康人和CHF患者RR间期信号EED分析结果示意图,其中数据长度为10000点,曲线表示均值,误差棒表示标准差。曲线上方的符号*表示两组人群能量T检验p<0.01。在层次选择上,去除了层次1和高于7的层次,包括余量。层次1容易受到噪声的影响,导致能量较大的波动,引起结果标准差过大,并且其频率可高至几kHz以上,因而没有明确的生理意义;高于7的层次反应了信号的长时节律,非常容易受到外界环境的影响,并且其频率很低,生理意义不明。图7中在分量低层次(层次2、3)上,CHF患者的归一化能量值明显高于健康人,在分量高层次(层次>5)上,发生相反的变化,健康人能量高于CHF患者能量;健康人能量在层次2~5比较平稳,层次大于5时,能量缓慢上升,而CHF患者在层次2~4迅速下降,之后趋于稳定;在4个层次(2、3、6、7)上,二者能量有显著性区别(p<0.01)。为了对比,本发明增加了替代数据(Healthy Surrogate,CHF Surrogate)的EED分析结果,如图7所示,替代数据是通过将原始数据随机化的方法生成,替代数据的能量分布随着尺度增加单调下降,相比较CHF患者,在小时间尺度上能量更高,在长时间尺度上能量更低。
为了评估EED曲线在低层次和高层次分解上的能量分布差异,计算各人群的能量差异值EDV,高的EDV值表示RRI信号更高的分量低层次能量分布和更低的分量高层次能量分布。计算得到健康人、CHF患者以及他们替代数据的EDV值,如表2所示。
表2健康人和CHF患者的EDV值
*代表健康人和CHF患者T检验结果p<0.0001。
**代表替代数据与其原始数据T检验结果p<0.0001。
从表2可以看出,健康人和CHF患者EDV值具有显著的差别,CHF患者EDV值比健康人高很多。健康人的EDV值小于0,说明在分量高层次有更高的能量,预示着高层次分量有着更高的调节强度。两个人群的EDV值与他们的替代数据EDV值均有显著性差别,人体HRV的EDV明显小于随机序列。
心率变异性(HRV)是指测量连续心动周期之间的时间变异数,准确地说,应该是测量连续出现的正常P-P间期之间的差异的变异数。但由于P波不如R波明显或P波顶端有时宽钝,所以对心率变异性信号的研究通常用与P-P间期相等的R-R间期信号(RRI)来代替。研究表明,HRV可做为植物神经***活动的无创性检测指标,尤其在判断某些心血管疾病的预后方面有着重要意义。
在临床与实际应用中,本发明提出可采用512(n=8)个连续RR间期的短时心搏信号用于上述的EED分析是有效的,并有一定的数据余量;由上述研究表明,EED分解达到7个分量层次(n=7)已有很好的结果,则所需数据量最少可取N=2n+1=27+1=256。