心脏病由于存在发病突然, 死亡率高的特点, 对人们生命健康造成了极大的威胁。因此寻找一种无创伤的、较为准确的方法诊断预测心脏类疾病是人们的迫切要求, 也是一项生物医电学科亟待解决的课题。目前临床心脏类疾病的诊断主要是分析心电图的特征, 一方面依赖于专业医生的诊断, 另一方面采用机器对采集到的心电信号进行分析, 而机器识别的方法目前较多的是线性分析方法, 如频谱分析^[1]。随着对心脏系统研究的深入, 科研人员发现心脏系统是一个典型的非线性系统。因此单纯采用线性分析方法无法有效描述心脏的运动状态, 也无法及时检测出心脏运动的非正常状态, 可能会延误心脏异常的诊断及预警。因此越来越多的科研人员试图采用非线性的分析方法研究心电信号的特征, 力图实现心脏疾病的早期准确诊断。如Babloyantz等对心脏电活动的非线性现象进行了深入的研究, 提出心脏电活动中也存在混沌特征^[2]。Goldberger等从临床分析方面得出, 心电信号具有分形的性质^[3]。刘雁等利用R-R Lorenz图的非线性特征分析心率变异^[4]。宋晋忠等针对心肌缺血的心电信号特征参数进行分析, 分析出其分形维数参数具有明显的差异^[5]。

从以上研究可以看出, 心电信号具有强的非线性特征, 而分形理论中的单重和多重分形特征作为分析非线性系统的一类重要方法, 用于心电信号的研究却较少。基于此背景, 本文以动态心电信号为研究对象分析其单重及多重分形特征。

1 相关理论基础 1.1 单重分形

单重分形特征又称自相似性, 是指某种过程的特征从不同的空间尺度或时间尺度来看都是相同的或表现为相似的。时间序列的单重分形特征通常是通过分析自相关函数和Hurst参数实现的。

严格意义上, 窦性心律的心电信号, 是一个随机过程, 其自相似性不是精确的自相似性, 只能考虑其统计的自相似特性。统计自相似是指一个随机变量在连续的尺度范围内测量所得的概率分布函数相互成比例关系, 又称为简单自相似。

1) 自相关函数

如果随机过程x_i存在期望值E(x_i)=μ和方差var(x_i), 对于∀i∈I(I=1, 2, 3…), x称为二阶过程, 则自相关函数就是k的函数^[6]

(1)

2) Hurst参数及R/S分析法

Hurst参数(简称H)作为度量时间序列的长相关性参数, 通常采用R/S分析法(又称为重标极差分析法)根据时间序列按照长相关特性的幂指数规律进行计算^[7]。

对于时间序列x_k(k=1, 2, 3…N), 将其等分成长度为n(≥3)的m(取N/n的整数部分)个连续的子序列X_α, 子序列的元素、均值、标准差分别记为x_{k, α}, e_{k, α}和S_k, _α(α=1, 2, …, m), (2)式则为时间序列x_k在时间增量长度为n上的重标极差

(2)

对(2)式两边取对数, 可得:

(3)

式中，R, S, n及C分别表示重新标度的极差、标准差、时间增量区间长度(样本观察值数目)、常数。根据(3)式即可计算出Hurst参数H。

根据统计学知识, H总是大于0。

当0≤H < 0.5时, 表示系统做反持续性分形布朗运动。表示时间序列下一刻的增量与上一刻的增量是负相关的。即上一刻心电信号电压的趋势是增加的, 则下一刻心电信号电压的趋势将减少;

当H=0.5时, 表示时间序列没有自相似性。即下一刻心电信号与前一刻的心电信号无关;

当0.5 < H < 1时, 表示系统具有长期持续性并做无周期的循环。这时意味着系统的变化是一个持续的时间序列。即前一刻的心电信号对下一刻的心电信号是有影响的, 存在着“记忆效应”, 即自相关性, 也称为时间序列具有单重分形特征。而且H值越接近1, 自相似程度越高。H越接近0.5, 其变化趋势也越不稳定, 自相似性也越弱^[8]。

因此, H参数可以作为描述心电信号变化的自相似参数, 表征其突发特性。

1.2 多重分形

通过多重分形可以描述分形体不同层次的特征, 即其局部奇异特性。进而从系统的局部出发来研究其整体特征。

把系统分为N个小区间, 设第i个小区间的线度大小为ε(ε→0), 分形体在该小区间的几率为P_i(ε), 将子集内奇异指数相同的单元数定义为N(ε), 则P_i(ε)和N(ε)满足下面的幂函数的子集:

(4)

(5)

式中, α是奇异指数, 反映分形上各个小区间的奇异程度。f(α)称为多重分形谱, 表示相同子集的分维数。

统计物理给出了规则和不规则多重分形谱的计算方法^[9]。首先, 定义一个配分函数χ_q(ε)≡∑P_i(ε)^q, 若配分函数收敛, 则等式

(6)

成立, 即可得到

(7)

τ(q)称为质量指数。设N(P_i)是概率为P_i(ε)的小区间的数目, 则有

(8)

将(4)、(5)、(8)式代入(6)式可以得出多重分形谱f(α)

(9)

1.3 心电信号

所分析的数据来自MIT-BIH数据库, 主要选取正常窦性心律及3种常见的心脏疾病:心房纤维性颤动、充血性心力衰竭和恶性室性心律失常的心电信号进行分析。信号采样周期为125 Hz, 数据长度为1 min。心电信号通常由1组P-QRS-T波组成。若QRS波群中的Q, R, S波形幅值较小时, 则分别由q, r, s表示^[10]。

一般正常心电信号在0.05~100 Hz频率范围内, 而90%的心电信号频谱能量集中在0.25 Hz到35 Hz之间, 其中频率最高的QRS综合波群频谱能量主要集中在2 Hz至20 Hz之间, 在12 Hz左右能量最大^[10]。由于所采集到的心电信号会受到高频波干扰, 为保证数据结果的准确性, 首先对原始信号进行滤波处理, 这里采用IIR巴特沃兹低通数字滤波器对数据进行滤波, 截至频率为45 Hz^[11]。图 1为经过滤波后各心电信号的时序图(每类选取1例, 但具有代表性)。从图 1可以看出, 不同类心电信号波形形状差异较大, 心率也相差较大。其中图 1a)的心率为79次/min。图 1b)的心率为97次/min, 呈QS波形。图 1c)的心率为104次/min, 呈RS或qRS波形。图 1d)的心率为74次/min, 波形较为复杂, 出现冠状T波和rsR型的特性。

本文从单重分形和多重分形的角度进一步分析这4类心电信号的非线性动力学特征, 旨在探索从分形的角度对心电信号进行定量分析的方法, 目的是实现对心脏疾病的早期无创准确的诊断。

2 心电信号的单重分形特征

这4类心电信号在一个心动周期内的自相关函数图(每类选取1例, 但具有代表性)如图 2所示。

从图 2可以看出, 正常窦性心律自相关函数迅速衰减, 在k=3就已经接近0, 表明其自相关特性较弱。而病态心电信号的自相关函数相对衰减较慢, 表现出明显的重尾特性。进一步分析心电信号的Hurst参数, 每类选取10个样本, 分析结果如表 1所示。

为检测Hurst参数在判断心动周期的应用效果, 定义α, 描述心动周期样本Hurst参数的波动程度; 定义β, 描述同类心动周期Hurst参数的稳定性。设H为样本的平均值, M为样本总数, N为样本中落到H的10%内的个数, 则:

(10)

(11)

从表 1可以看出H均在0.5~1之间, 表示此类信号均具有单重分形特征, 即具有状态的持续性和自相似性。进一步研究可以看出, 各类心电信号的β值都大于等于90%, 说明同类型心电信号的Hurst参数十分接近, 可以采用Hurst参数分析其单重分形特征。此外, 正常窦性心律的H值接近0.5, 而病态心电信号的H值更接近1, 尤其是恶性室性心律失常。这表示正常窦性心律的单重分形特征较弱, 变化趋势较不稳定, 表现出随机特征, 这一分析结果与自相关函数一致。因此, 我们可以初步根据一个心电信号的Hurst值判断该心电信号的类型。

3 心电信号的多重分形特征

根据单重分形分析心电信号的特征能够较为容易地分离出正常窦性心律, 但是对于非正常窦性心律, 由于存在个体差异, 不同类型的非正常窦性心律的Hurst参数存在交叠区域, 因此, 单纯依靠Hurst参数无法准确判断窦性心律类型, 需进一步分析各类心电信号的多重分形特征。

图 3为选取的4类心电信号的全部样本(每类10个)的多重分形谱图, 表 2为其多重分形谱参数。在表 2中, f_max=max(f(α_max), f(α_min))。为详细描述心电信号多重分形谱图的分类特征, 定义Δα和Δf(α)

(12)

(13)

定义K_Φ为Δf(α)出现正号的概率, 即f(α_max)大于f(α_min)的概率。设样本总数为N₁, 这里取10个样本, 出现Δf(α)为正的样本数为N₂, 则

(14)

对于心电信号的多重分形谱图, Δα值的大小反应心电信号电压概率分布均匀性, 即其混沌性的高低。这里Δα值越大, 说明电压分布越不均匀, 波动的强度越大, 混沌性也相应增强。结合表 2与图 3, 可以看出, 恶性室性心律失常的谱图Δα值相对较小, 与其他3种窦性心律相比较, 其他3种窦性心律电压分布不均匀, 波动量较大, 局部奇异特征较强。而Δf(α)的符号反映心电信号处于高值和低值的概率, 进而反映心动周期QRS波群中各波的强弱。当多重分形谱图向左钩时, 即Δf(α)>0时, 表明心电信号处于低值的机会大于高值, 反之亦然, 即多重分形谱的参数可以在一定程度上反映心电信号的变化范围和高低电压出现频率的变化。观察4种心电信号可知, 正常窦性心律的多重分形谱图均为向左钩, 恶性室性心律失常的谱图大部分也向左钩, 心房纤维性颤动的谱图大部分向右钩, 而充血性心力衰竭的谱图向左钩和向右钩的概率相同。

进一步, 发现同类型的心电信号的多重分形谱图出现了分类现象, 如正常窦性心律的谱图在向左钩的基础上可以分为2类, 而这2类与其他类型的窦性心律谱图没有重叠现象。而心房纤维性颤动、充血性心力衰竭、恶性室性心律失常可能由于病变原因不同, 其多重分形谱图并非都是向左钩或都是向右钩。为使图形特征更为明显, 我们根据图形的形状和α与f(α)的取值范围将每一类心电信号的样本的多重分形谱图做了分类, 分类的标准为:

1) Δf(α)的符号一致;

2) α_min及α_max值与中心值相差不超过0.01;

3) f_max值相差不超过0.2。

当同时满足以上3个条件时, 即认为心电信号的多重分形谱相似, 并将其划为一类。依据以上的分类标准, 将样本数据进行分类, 如表 2所示。图 4为相同类型相似样本根分类的多重分形谱图, 图 5为不同类型样本相似多重分形谱图。从图 4和图 5就可以较清晰地看出, 各类心电信号的多重分形谱图特点:

正常窦性心律(nsr)10个样本均是向左钩, 根据Δα和Δf(α)可以分为2类。而其他3种心电信号的向左钩的谱图均不处于这2类谱图中。

心房纤维性颤动(af)的谱图中, 向左钩有3个样本, 均不相似, 并且其中1个和充血性心力衰竭的向左钩样本相近(af2, chf10), 同时没有落入充血性心力衰竭的向左钩样本的相似类型中。而向右钩的7个样本, 5个相似, 还有2个样本与充血性心力衰竭相近, 但也没有落入其相似类型中。(af1, af4剩余2个与chf3未归类的并不相似)。已归类的af8, af10与已归类的chf2, chf8相近。

充血性心力衰竭(chf)的向右钩的5个样本中, 4个样本相似, 而向左钩的5个样本中3个样本特征相近, 其他2个样本, 其中一个(chf10)与af2相似, 另一个与其他样本均不相似。

恶性室性心律失常(mva)的向左钩的7个样本中, 4个样本十分接近, 还有3个样本与其余样本均不相似, 而向右钩的3个样本就十分接近了。

根据所设定的心电信号多重分形谱的分类划分标准, 所设定的4类心电信号的类型没有重叠。而能否根据心电信号的多重分形谱满足相对应的相似类型即可判定窦性心律的类型, 还需要海量数据的分析及判定标准的改进, 如中心值标准的设置方法的确定等。这也是我们未来要研究的内容之一。

通过分析关于临床对心脏疾病的研究, 引起心律失常的原因有多种, 如:风心病、冠心病或甲亢性心脏病等可能引发心房纤维性颤动, 冠心病、心瓣膜、高血压等可能引发充血性心力衰竭。而冠心病等心脑血管疾病也可能引发恶性恶性室性心律失常^[11]。由于同一病症如冠心病, 可能引起不同的心律失常, 而其心电信号的多重分形谱是否会有相似的特点也是我们未来研究的一项内容。

4 结论

论文以MIT-BIH数据库中的4类心电信号样本为研究对象, 分别从单重分形和多重分形的角度对心电信号时间序列的分形特征进行详细的分析。研究结果表明, 不同窦性心律的心电信号中存在自相似特征的同时, 也存在局部奇异特性。而且, 对于同一类型的心电信号多重分形谱出现分类现象。

无论是单重分形的自相关函数还是多重分形的谱参数, 对于同一类心电信号都存在相似的关联特征。尤其是多重分形谱参数的差异与心电信号的类型密切相关。因此采用分形理论可以对于心电信号的特征进行识别。而MIT-BIH数据库内的心电信号数据量有限, 本文仅对一些现有的数据进行分形特征分析, 而并没有结合发病的起因进行分析, 这便使得分析结果存在一定的局限性, 这也从心电信号的多重分形谱图中可以看出。后续工作可集中于海量数据, 针对每一类型的心电信号, 结合发病起因, 对其分形特征进一步分析, 为定量诊断提供一定的可靠的分析依据。为心脏疾病的无创诊断和预警, 提供一有效的分析手段。