摘要:心率变异性的复杂波动反映了心脏的自主调节功能. 本文提出了一种新的心率变异性度量方法—ICBN方法, 该方法通过改进的自适应噪声完备集合经验模态分解方法对心率变异性信号进行分解, 得到多个模态分量, 计算每个模态分量的bubble熵得到熵值向量, 把该向量映射成复杂网络, 通过计算网络的特征参数, 对心率变异性在不同时频尺度状态下的非线性特征之间的耦合关系进行度量. 首先, 采用时域、频域和ICBN分析方法对29名充血性心力衰竭病人和29名正常窦性心律对象的心率变异性进行分析, 结果表明: 时域指标三角指数HRVTi, 频域指标LF/HF, 网络层级加权值WB, 平均点权值PW, 特征路径长度CL具有统计学差异; 基于网络层级加权值WB, 特征路径长度CL, 频域指标LF/HF和Fisher判别方法的识别模型对充血性心力衰竭病人的识别正确率达到89.66%. 然后, 又对43名房颤心律失常患者和43名正常窦性心律对象的心率变异性进行分析, 结果表明: 时域指标SDNN, pNN50, RMSSD, 频域指标LF/HF, 网络层级加权值WB, 平均点权值PW具有统计学差异; 时域指标pNN50, RMSSD, 频域指标LF/HF和网络层级加权值WB, 平均点权值PW作为特征向量, Fisher判别方法作为分类器, 对房颤心律失常患者的识别正确率达到91.86%. 综合以上实验结果可知, 本文为心率变异性的度量研究提供了一种新的思路.
关键词: 心率变异性/
改进的自适应噪声完备集合经验模态分解(ICEEMDAN)/
bubble熵/
复杂网络

English Abstract

全文HTML

--> --> -->

心率变异性(heart rate variability, HRV)是指逐次心跳间期的波动, 常用来评估自主神经系统(autonomic nervous system, ANS)对心脏活动的调节作用, 是评价自主神经活动的重要无创指标, HRV的分析能够了解ANS在心脏发病过程中的作用^[1-3].
传统的HRV分析方法集中于时域和频域分析^[4-6]. 时域分析法主要是利用相邻正常心跳间隔的标准差(standard deviation of normal to normal, SDNN)、相邻正常心跳间隔差值的均方根(root mean square of the successive normal, RMSSD)、相邻正常心跳间隔差值超过50 ms的比例(count divided by the total number of all NN intervals, pNN50)、相邻正常心跳间隔的总个数除以相邻正常心跳间隔直方图的高度(HRV triangular index, HRVTi)等指标对HRV进行评价. 频域分析法是把HRV分解为不同的频率成分, 并将其相对强度定量为功率, 提供了各种频率成分的功率谱测定, 分为甚低频段(very low frequency, VLF)、低频段(low frequency, LF)、高频段(high frequency, HF)、总功率谱(total power, TP)等指标. 然而, 时域指标不能表达出HRV信号的时变特性, 所以对自主神经系统的反映十分有限, 而频域分析只能给出全局的频率信息, 缺少局部和不同频率间的耦合信息. 时域和频域的分析方法都是对HRV信号做线性分析, 很难透过HRV信号的非线性本质^[7,8]. 因此, 通过非线性分析方法对HRV进行分析得到了广泛的研究, 包括Poincaré散点图^[9]、去趋势波动分析^[10]和近似熵^[11]、样本熵^[12]、排列熵^[13]、bubble熵^[14]等熵分析法. 熵分析法是HRV分析中比较有代表性的一类方法, 因HRV的无序波动中, 包含了心脏活动中有序的内在动力学特性, 而熵测度能够度量心脏活动的这种有序程度, 所以在HRV信号的研究中, 熵分析法得到了广泛的应用. 但对HRV信号直接进行熵测度分析忽略了信号内部的波动细节. 为了增强对HRV信号内部细节的观察, 文献[15]把经验模态分解(empirical mode decomposition, EMD)方法与熵分析方法相结合, 但该方法忽略了残差信号, 同时没有考虑不同本征模态函数(intrinsic mode function, IMF)之间的耦合关系, 造成了信息量的丢失, 而且EMD分解算法存在模态混叠和端点效应等不足.
目前, 对HRV信号的分析研究中, 并没有对不同时频尺度下的HRV信号之间的耦合关系进行分析. 本文根据非线性和非平稳的HRV信号在不同时刻具有不同频率成分的特点, 提出了一种基于改进的自适应噪声完备集合经验模态分解(improved complete ensemble EMD with adaptive noise, ICEEMDAN)和bubble熵(bubble entropy)的复杂网络(complex network, CN)映射的HRV度量方法—ICBN方法, 该方法主要研究不同时频尺度空间下的HRV信号之间的耦合关系. 本研究先通过ICEEMDAN对HRV信号进行时频尺度化, 然后将得到的不同时频段间的IMF分别计算bubble熵值, 构建熵值向量, 根据有限穿越水平可视图法(limited penetrable horizontal visibility graph, LPHVG)^[16,17]将熵值向量映射成复杂网络, 通过网络参数特征对不同时频尺度空间化下的IMF间的耦合关系进行度量, 从而实现对HRV信号的度量.

本文的实验方法结构如图1所示. 首先, 对带有人工校准的心跳标注和时间戳的文件进行HRV提取, 并进行片段分割; 然后, 通过时域分析、频域分析、ICBN方法对HRV信号进行分析, 提取时域、频域指标和耦合特征指标; 最后, 通过统计分析和交叉校验对各特征的分类性能进行对比分析.
图 1 HRV分析框架图
Figure1. Framework of the HRV analysis.

2

-->

2.1.研究对象

实验数据取自PhysioNet提供的数据集^[18]. 在这里共使用了4个数据集. 第一个数据集是chf2, 包含29个患有充血性心力衰竭(congestive heart failure, CHF)病人的长时心电信号(electrocardiogram, ECG), 患者年龄从34到79岁. 第二个数据集是nsr2, 包含54个正常窦性心律(normal sinus rhythm, NSR)对象的长时ECG信号, 其中男性30名, 年龄从28.5到76岁, 女性24名, 年龄从58到73岁. 第三个数据集是nsr, 包含18个正常窦性心律的长时ECG信号, 包括男性5名, 年龄从26到45岁, 女性13名, 年龄从20到50岁. 第四个数据集是ltaf, 包含84个房颤心律失常(atrial fibrillation, AF)患者的长时ECG信号. 原始ECG信号的采样频率是128 Hz, 每个信号都带有人工校准的心跳标注, 从而可以获得用于数据分析的RR间期序列, 即HRV信号(图2). 为了排除昼夜节律对信号分析的影响, 文献[19]采用上午清醒时间段的信号进行分析, 本文采用同样的方法进行实验.
图 2 HRV信号获得示意图
Figure2. Schematic of obtaining HRV signal.

2 -->

2.2.时域、频域分析方法

时域分析的统计特征如表1所示, 包括SDNN, RMSSD, pNN50, HRVTi这4个时域特征指标.

指标	单位	描述与定义
SDNN	${\rm{ms}}$	相邻正常心跳间隔的标准差${\rm{SDNN}} = {\sqrt {\dfrac{1}{N}\displaystyle \sum \limits_{i = 1}^N \left( {{\rm{RR}}{s_i} - \dfrac{1}{N}\mathop \sum \limits_{i = 1}^N {\rm{RR}}{s_i}} \right)} ^{}}$
RMSSD	ms	相邻正常心跳间隔差值的平方和均值的均方根${\rm{RMSSD}} = \sqrt {\dfrac{1}{{N - 1}}\displaystyle \sum \limits_{i = 1}^{N - 1} {{\left( {{\rm{RR}}{s_{i + 1}} - {\rm{RR}}{s_i}} \right)}^2}} $
pNN50	%	相邻正常心跳间隔差值超过50毫秒的比例${\rm{PNN}}50 = \dfrac{{{\rm{num}}\left[ {\left( {{\rm{RR}}{s_{i + 1}} - {\rm{RR}}{s_i}} \right) > 50{\rm{ ms}}} \right]}}{{N - 1}}$
HRVTi	—	相邻正常心跳间隔的总个数除以相邻正常心跳间隔直方图的高度

表1时域分析统计特征
Table1.Statistical features in time domain.

本文采用Lomb-Scamble周期图法计算功率谱密度, 相对基于FFT的方法具有更高的估计精度^[20]. 本文计算的频域特征包括TP (≤ 0.4 Hz), VLF (0.003—0.040 Hz), LF (0.04—0.15 Hz)和 HF(0.15—0.40 Hz)范围的功率谱, LF和HF的比值LF/HF, 5个频域特征指标. 具体的描述与定义如表2所示.

指标	单位	描述与定义	频率范围
Total power	${\rm{m}}{{\rm{s}}^2}$	所有频率范围的功率谱总和	≤ 0.4 Hz
VLF	${\rm{m}}{{\rm{s}}^2}$	甚低频范围内的功率谱	0.003—0.040 Hz
LF	${\rm{m}}{{\rm{s}}^2}$	低频范围内的功率谱	0.04—0.15 Hz
HF	${\rm{m}}{{\rm{s}}^2}$	高频范围内的功率谱	0.15—0.40 Hz
LF/HF	%	LF$\left[ {{\rm{m}}{{\rm{s}}^2}} \right]$ 与HF$\left[ {{\rm{m}}{{\rm{s}}^2}} \right]$的比值	—

表2频域分析统计特征
Table2.Statistical features in frequency domain.

2 -->

2.3.ICBN分析方法

本文提出的ICBN分析方法的实现过程如图3所示. 首先, 读取带有人工校准标注的数据文件并提取HRV信号, 根据标注文件剔除异常点; 然后, 通过ICEEMDAN对提取的HRV信号在不同时频尺度下进行分解, 得到不同时频尺度空间下的IMF; 再通过计算每个IMF的bubble熵值, 构建bubble熵向量, 对熵向量进行网络的映射; 最后, 通过对网络的特征进行计算, 完成对HRV信号在不同时频尺度下的耦合特征分析.
图 3 ICBN分析方法实现过程框图
Figure3. Implementation process of ICBN analysis method.

3 -->

2.3.1.RRs提取

获取带有人工标注的ECG信号文件和人工标注的注释文件, 提取每个对象上午清醒时间段1 h长的HRV信号, 定义为$x\left( i \right), 1 \leqslant i \leqslant N,\;i \in Z$, 其中$N$表示信号长度.
3 -->

2.3.2.ICEEMDAN分解

ICEEMDAN基于EMD方法提出, 是针对非线性和非稳态数据的自适应分析方法, 它改善了EMD方法的模态混叠效应, 减少了IMF中的噪声成分, 同时丰富了IMF的物理意义^[21]. 对于输入HRV信号$x$, ICEEMDAN对其分解的具体实现过程如下.
1)对于$j = 1, 2, \cdots, J$, 通过EMD分解${x^{\left( j \right)}} = x + {\beta _0}{E_1}\left( {{w^{\left( j \right)}}} \right)$, 得到第一个残差:

${r_1} = \left\langle {M\left( {{x^{\left( j \right)}}} \right)} \right\rangle .$

2)在第一阶段$k = 1$, 计算第一个模态分量:

${\rm{IM}}{{\rm{F}}_1} = x - {r_1}.$

3)对于$j = 1, 2, \cdots, J$, 通过EMD分解${r_1} + {\beta _0}{E_2}\left( {{w^{\left( j \right)}}} \right)$, 得到第二个模态分量为

${\rm{IM}}{{\rm{F}}_2} = {r_1} - {r_2} = {r_1} - \left\langle {M\left( {{r_1} + {\beta _1}{E_2}\left( {{w^{\left( j \right)}}} \right)} \right)} \right\rangle .$

4)对于$k = 3, \cdots, K$, 计算第$k$个残差:

${r_k} = \left\langle {M\left( {{r_{k - 1}} + {\beta _{k - 1}}{E_k}\left( {{w^{\left( j \right)}}} \right)} \right)} \right\rangle .$

5)计算第$k$个模态分量:

${\rm{IM}}{{\rm{F}}_k} = {r_{k - 1}} - {r_k}.$

6)对下一个$k$回到步骤4).
重复4)到6), 直到获取的残差不能继续被EMD分解, 则输入信号被分解为

$x = \sum\limits_{k = 1}^K {{\rm{IM}}{{\rm{F}}_k}}.$

以上公式中${E_k}\left( \cdot \right)$为产生第$k$个模态分量的算子; $M\left( \cdot \right)$是局部均值算子; ${w^{\left( i \right)}}$为零均值单位方差白噪声; 系数${\beta _k} = {\varepsilon _k}{\rm{std}}\left( {{r_k}} \right)$允许在每一次迭代中选择信噪比, 其中${\rm{std}}\left( \cdot \right)$为标准偏差算子. 研究表明^[21], 噪声幅值${\varepsilon _0} = 0.2$时, ICEEMDAN具有最好的分解效果, 分解次数$J$在取值50和200时具有相似的分解结果, 所以本文设置${\varepsilon _0} = 0.2$, $J = 50$.
3 -->

2.3.3.计算IMFs的bubble熵

Bubble熵是一种更加稳定和精确的度量方法, 不需要设置容限值且度量结果不依赖于嵌入维数, 是基于条件排列熵和Renyi排列熵而提出的. 相对于样本熵和近似熵, bubble熵具有更好的区分能力^[14]. 对于输入的HRV信号$x:\left\{ {x\left( i \right):1 \leqslant i \leqslant N} \right\}$, 通过ICEEMDAN在不同时频尺度上分解出不同的模态分量${\rm{IMFs}}$, 定义为${\rm{IMFs}} = [{\rm{IM}}{{\rm{F}}_1}, {\rm{IM}}{{\rm{F}}_2}, \cdots, $${\rm{IM}}{{\rm{F}}_K}] $, 其中每个模态分量定义为${\rm{IM}}{{\rm{F}}_k}:\{{\rm im}{f_k}\left( i \right): $$1 \leqslant i \leqslant N, 1 \leqslant k \leqslant K \}$, 对${\rm{IM}}{{\rm{F}}_k}$进行计算bubble熵的过程如下.
将${\rm{IM}}{{\rm{F}}_k}$嵌入一个$m$维空间, 对于每一个$im{f_k}\left( i \right)$, 进行相空间重构, 得到

$\begin{split}{\rm{IM}}{{\rm{F}}_k}\left( i \right) = & \left[{\rm{im}}{f_k}\left( i \right), {\rm{im}}{f_k}\left( {i + L} \right), \cdots,\right. \\ & \left.{\rm{im}}{f_k}\left( {i + \left( {m - 1} \right)L} \right) \right], \end{split}$

其中$L$和$m$分别为延迟因子和嵌入维数. 将${\rm{IM}}{{\rm{F}}_k}\left( i \right)$的$m$个重构分量$\left[{\rm{im}}{f_k}\left( i \right), {\rm{im}}{f_k}\left( {i + L} \right), \cdots,\right.$$\left.{\rm{im}}{f_k}\left( {i + \left( {m - 1} \right)L} \right)\right] $使用冒泡排序法, 按照升序进行重新排列, 统计每一个重构分量的必要交换次数${n_j}\left( {j = 1, 2, \cdots, m} \right)$, 在bubble熵的计算中, $L = 1$, 所以$m$维矢量的个数为$N - m + 1$, 计算其对应次数出现的概率$p_j^{} = \dfrac{{{n_j}}}{{N - m + 1}}$. 嵌入维数m根据文献[14]设置为15.
Renyi排列熵可根据如下公式计算:

${H_2} = {\rm{ - }}\log \left( {\sum\limits_{i = 1}^n {p_j^2} } \right).$

条件排列熵可根据如下公式计算:

${\rm{cpe}} = {\rm{PE}}\left( {m + 1} \right) - {\rm{PE}}\left( m \right),$

其中, ${\rm{PE}}\left( m \right)$代表m维空间排列熵, ${\rm{PE}}\left( {m + 1} \right)$代表m+1维空间排列熵.
条件Renyi排列熵为

${\rm{crpe}} = \left( {H_2^{m + 1} - H_2^m} \right)/\log \left( {m + 1} \right).$

根据公式(8), (9)和(10)得bubble熵的计算公式为

${\rm{be}} = \left( {H_{{\rm{swaps}}}^{m + 1} - H_{{\rm{swaps}}}^m} \right)/\log \left( {m + 1/m - 1} \right), $

其中$1/\log \left( {m + 1/m - 1} \right)$是归一化因子.
对于模态分量${\rm{IMFs}} = \left[ {{\rm{IM}}{{\rm{F}}_1}, {\rm{IM}}{{\rm{F}}_2}, \cdots, {\rm{IM}}{{\rm{F}}_K}} \right]$, 得到bubble熵值向量${\rm{BE}} = \left[ {{\rm{b}}{{\rm{e}}_1}, {\rm{b}}{{\rm{e}}_2}, \cdots, {\rm{b}}{{\rm{e}}_K}} \right]$.
3 -->

2.3.4.Bubble熵网络的映射

LPHVG具有相对更强的连接性, 进而在较小的时间尺度上可以更多地反应信号的内在信息和特征, 同时不影响对序列长程波动趋势的影响^[16,17], 因此本文选择LPHVG方法构建bubble熵网络.
对于时间序列$x\left( i \right),1 \leqslant i \leqslant N,i \in Z$, 通过ICEEMDAN在不同时频尺度上分解出IMFs, ${\rm{IMF}}s = \left[ {{\rm{im}}{f_1}, {\rm{im}}{f_2}, \cdots, {\rm{im}}{f_K}} \right]$, 对每一个IMF进行非线性分析计算出bubble熵值, 得到向量${\rm{BE}} = \left[ {{\rm{b}}{{\rm{e}}_1}, {\rm{b}}{{\rm{e}}_2}, \cdots, {\rm{b}}{{\rm{e}}_K}} \right]$, 根据LPHVG, 将${\rm{BE}}$向量映射成网络$G = \left( {V, \xi } \right)$. 通过对网络$G = \left( {V, \xi } \right)$的特征参数进行计算, 完成对IMFs之间的耦合关系分析, 其中节点集合为$V = \left\{ {\left. {{v_1},{v_2},\cdots,{v_K}} \right\}} \right.$, $\xi \subseteq $$\left\{ {\left. {\left( {{v_\delta }, {v_\vartheta }} \right):{v_\delta }, {v_\vartheta } \in V, {v_\delta } \ne {v_\vartheta }} \right\}} \right.$是边所构成的集合, 且$\left( {{v_\delta }, {v_\vartheta }} \right) \in \xi \Leftrightarrow \left( {{v_\vartheta }, {v_\delta }} \right) \in \xi $, $\left| V \right|$表示图$G$的阶. 定义$W = \left[ {{w_{\delta \vartheta }}} \right]$为图$G$的加权邻接矩阵, 其中${w_{\delta \vartheta }} \ne 0 \Leftrightarrow \left( {{v_\delta }, {v_\vartheta }} \right) \in \xi $, 2个节点${v_\delta }$和${v_\vartheta }$的权值分别为${\rm{b}}{{\rm{e}}_\delta }$和${\rm{b}}{{\rm{e}}_\vartheta }$, 定义节点${v_\delta }$和${v_\vartheta }$之间边的权值为

${w_{\delta \vartheta }} = \left\{ \begin{aligned}&\left| {\left( {{\rm{b}}{{\rm{e}}_\delta } - {\rm{b}}{{\rm{e}}_\vartheta }} \right)} \right|, \;{\text{节点}}{v_\delta }{\text{和}}{v_\vartheta }{\text{相连}},\\& 0,\quad \quad \quad \qquad \;\;{\text{节点}}{v_\delta }{\text{和}}{v_\vartheta }{\text{不相连}}.\end{aligned} \right.$

节点${v_\delta }$的邻居集合定义为${T_\delta } =\{{v_\vartheta } \in V: $${w_{\delta \vartheta }} \ne 0 \} =\left\{ {\left. {{v_\vartheta } \in V:\left( {{v_\delta }, {v_\vartheta }} \right) \in \xi } \right\}} \right.$, 节点${v_\delta }$的度定义为${d_\delta }$, 则${d_\delta } = \left| {{T_\delta }} \right|$, 图$G$的度矩阵定义为${{{D}}_G} = \left[ {{d_1}, {d_2}, \cdots, {d_K}} \right] = \left[ {\left| {{T_1}} \right|, \left| {{T_2}} \right|, \cdots, \left| {{T_K}} \right|} \right]$, 根据节点之间的边权值, 节点${v_\delta }$的点权值${S_\delta }$定义为与它关联的边权值之和, 即${s_\delta } = \sum\limits_{\vartheta = 1}^K {{a_{\delta \vartheta }}{w_{\delta \vartheta }}} $, 图$G$的点权矩阵定义为${{{S}}_G} = \left[ {{s_1}, {s_2}, \cdots, {s_K}} \right]$.
3 -->

2.3.5.计算网络特征指标

在将bubble熵值向量映射为加权网络后, 对5个复杂网络参数进行考察, 分别为平均集聚系数${\rm{MC}}$, 特征路径长度${\rm{CL}}$, 网络拓扑熵${\rm{TE}}$, 网络层级加权值${\rm{WB}}$, 平均点权值${\rm{PW}}$.
平均集聚系数${\rm{MC}}$:　平均集聚系数是指在网络中与同一个节点连接的2个节点之间相互连接的平均概率, 该系数用来刻画网络的局域结构性质^[22].

${\rm{MC}} = \frac{1}{K}\sum\limits_{\delta = 1}^K {\frac{{2\left| {{T_\delta }} \right|}}{{{\kappa _\delta }\left( {{\kappa _\delta } - 1} \right)}}}, $

其中, 与节点${v_\delta }$直接相连接的节点数为$\left| {{T_\delta }} \right|$, ${\kappa _\delta }$为所有与节点${v_\delta }$相连接的节点之间直接相连接的边数.
特征路径长度${\rm{CL}}$:　特征路径长度表示所有节点之间距离的平均值, 它描述了网络中节点间的平均分离程度, 即网络有多小^[22].

${\rm{CL}} = \frac{2}{{K\left( {K - 1} \right)}}\sum\limits_{\delta = 1}^K {\sum\limits_{\vartheta = \delta + 1}^K {{w_{\delta \vartheta }}} } .$

网络标准拓扑熵TE:　网络拓扑熵可以更简洁的度量网络的序状态, 对网络拓扑熵进行归一化得到网络标准拓扑熵, 其定义为^[23]

${\rm{TE}} = \frac{{ - 2 \displaystyle\sum\limits_{\delta = 1}^K {{\gamma _\delta } \cdot \ln {\gamma _\delta }} - \ln 4\left( {K - 1} \right)}}{{2\ln K - \ln 4\left( {K - 1} \right)}},$

其中${\gamma _\delta } = {{{d_\delta }}\bigg/ {\displaystyle\sum\limits_{\delta = 1}^K {{d_\delta }} }}$为节点重要程度, 显然$0 \leqslant {\rm{TE}} \leqslant 1$.
针对IMFs的bubble熵网络的实际物理意义, 本文在构建的bubble熵网络的基础上提出2个网络特征, 分别是网络层级加权值${\rm{WB}}$和平均点权值${\rm{PW}}$.
网络层级加权值WB:　网络层级加权值以不同权值, 对不同层级地位的网络节点的权值总量进行计算, 图$G$中节点${v_\delta }$的$\delta $表示IMF的层级位置, 每个节点${v_\delta }$的权值为${\rm{b}}{{\rm{e}}_\delta }$, 则网络层级加权值为

${\rm{WB}} = \sum\limits_{\delta = 1}^K {{\rm{b}}{{\rm{e}}_\delta }} \cdot \delta .$

平均点权值PW:　用来刻画网络中点权值的平均分布情况, 定义为图$G$的点权向量${{{S}}_G} = \left[ {{s_1}, {s_2}, \cdots, {s_K}} \right]$除以度向量${{{D}}_G} = \left[ {{d_1}, {d_2}, \cdots, {d_K}} \right]$.

$\begin{split}& {\rm{PW}} = \frac{{{{{S}}_G}}}{{{{{D}}_G}}} = \\ & \frac{{\left[ {\displaystyle\sum\limits_{\vartheta = 1}^K {{a_{1\vartheta }}{w_{1\vartheta }}}, \sum\limits_{\vartheta = 1}^K {{a_{2\vartheta }}{w_{2\vartheta }}, \cdots, \sum\limits_{\vartheta = 1}^K {{a_{K\vartheta }}{w_{K\vartheta }}} } } \right]}}{{\left[ {\left| {{T_1}} \right|, \left| {{T_2}} \right|, ..., \left| {{T_K}} \right|} \right]}}.\end{split}$

采用上述的时域、频域分析方法和本文提出的ICBN分析方法, 分别对不同人群的HRV信号进行分析, 然后通过T检验对各类特征分别进行差异分析, 并取$p < 0.001$的特征作为最终的极显著差异特征进行分类, 最后通过Fisher判别模型对实验对象进行识别. 为了得到更加可靠、稳定的结果, 采用留一法交叉检验实验的平均值作为最后的识别结果.
2

-->

3.1.应用于健康人和CHF患者的HRV

首先对如下两组HRV信号进行分析, 即CHF组和NSR1组, 其中CHF组采用chf2数据集的29例CHF患者; NSR1组采用nsr2数据集的前29例健康人.
3 -->

3.1.1.NSR1和CHF的HRV信号差异性分析

HRV信号时域、频域特征和ICBN分析方法计算的IMF间耦合特征的均值和标准差如表3所示. 从NSR1和CHF这2类人群的对比中可以看出, WB, PW, CL, HRVTi, LF/HF具有极显著统计学差异($p < 0.001$); TE, MC指标具有显著统计学差异($p < 0.01$). 其中, 耦合特征中WB, PW和CL具有极显著的统计学差异($p < 0.001$); 时频域指标中HRVTi和LF/HF具有极显著的统计学差异($p < 0.001$).

指标		NSR1(mean$ \pm $SD)	CHF(mean$ \pm $SD)	标准误差差值的95%置信区间		$p$
				下限	上限
ICBN	WB	21.853$ \pm $1.479	27.835$ \pm $7.741	3.050	8.914	0***
	PW	0.563$ \pm $0.051	0.455$ \pm $0.103	–0.151	–0.066	0***
	TE	0.956$ \pm $0.019	0.937$ \pm $0.037	–0.034	–0.005	0.009**
	CL	1.135$ \pm $0.133	0.954$ \pm $0.194	–0.268	–0.094	0***
	MC	0.684$ \pm $0.035	0.705$ \pm $0.026	–0.038	–0.006	0.009**
时域	SDNN	81.507$ \pm $38.566	59.535$ \pm $44.76	–43.951	0.007	0.05
	pNN50	11.476$ \pm $14.676	10.772$ \pm $14.110	–8.277	6.870	0.853
	RMSSD	51.172$ \pm $54.895	60.307$ \pm $58.497	–20.707	38.976	0.542
	HRVTi	6.886$ \pm $2.452	4.093$ \pm $1.494	–3.861	–1.725	0***
频域	TP	1.809$ \pm $4.909	1.443$ \pm $4.052	–2.734	2.001	0.758
	VLF	0.0003$ \pm $0.448	1.387$ \pm $6.214	–1.265	3.370	0.367
	LF	0.212$ \pm $0.333	0.119$ \pm $0.316	–0.263	0.078	0.281
	HF	1.597$ \pm $4.608	1.322$ \pm $3.745	–2.483	1.934	0.804
	LF/HF	0.288$ \pm $0.184	0.108$ \pm $0.083	–0.255	–0.105	0***
注: *, **, ***分别代表$p < 0.05$$p < 0.01$, $p < 0.001$.

表3NSR1和CHF患者不同分析方法下的结果
Table3.Statistical analysis results of HRV index under different analysis methods.

图4是具有极显著差异的特征的均值和标准差. 从图4中可以看出, 在NSR1组中, PW, CL, HRVTi, LF/HF的均值比CHF组的均值高, PW的均值比CHF组的均值低. ICBN方法中的指标WB, PW和CL, 在CHF对象中表现出了更高的标准差, 而时域指标HRVTi和频域指标LF/HF, 在NSR1对象中表现出了更高的标准差.
图 4 NSR1和CHF这2组对象具有极显著差异的指标的均值与标准差
Figure4. Mean and variance of the indicators with very significant differences between the two groups of NSR1 and CHF.

3 -->

3.1.2.CHF病人的识别

通过Fisher判别模型对NSR1和CHF这2组HRV信号进行分类, 实现对不同指标分类性能的验证. 分类结果通过正确率(accuracy, Acc)、灵敏度(sensitivity, Sen)、特异度(specificity, Spe)、ROC曲线下的面积(area under curve, AUC)来进行评价. 实验结果如表4所列, 通过单独的每个特征进行分类时, WB, PW, CL, HRVTi和LF/HF这5个特征之中, WB表现出了更高的Acc, Sen和AUC, 分别是79.3%, 90.48%和81.72%, LF/HF表现出了更高的Spe, 为80.95%.

指标	TP	TN	FP	FN	Acc	Sen	Spe	AUC
WB	19	27	10	2	79.31	90.48	72.97	81.72
PW	19	24	10	5	74.14	79.17	70.59	75.53
CL	20	23	9	6	74.14	76.92	71.88	74.40
HRVTi	23	19	6	10	72.41	69.70	76.00	72.85
LF/HF	25	17	4	12	72.41	75.68	80.95	78.32
注: TP, 被判定为CHF病人的数量; TN, 被判定为NSR1对象的数量; FP, NSR1对象被判定为CHF病人的数量; FN, CHF病人被判定为NSR1对象的数量; 正确率${\rm{Acc}} = \dfrac{{{\rm{TP + TN}}}}{{{\rm{TP \!+\! FP \!+\! TN \!+\! FN}}}} \times 100\% $; 灵敏度${\rm{Sen}} = \dfrac{{T{\rm{P}}}}{{T{\rm{P}} + {\rm{FN}}}} \times 100\% $; 特异度${\rm{Spe}} = \dfrac{{{\rm{TN}}}}{{{\rm{FP \!+\! TN}}}} \times 100\% $; ${\rm{AUC}} = \dfrac{1}{2}\left( {\dfrac{{{\rm{TP}}}}{{{\rm{TP + FN}}}}{\rm{ + }}\dfrac{{{\rm{TN}}}}{{{\rm{TN + FP}}}}} \right) \times 100\% $.

表4不同特征的CHF识别性能对比
Table4.Performance comparisons of different indices for CHF recognition

本文通过穷举法, 将具有极显著差异的5个指标WB, PW, CL, HRVTi和LF/HF中的所有2, 3, 4, 5个指标的组合作为特征, 通过Fisher判别模型对NSR1和CHF病人进行分类, 将分类正确率最高的5种组合列出(具体见表5). 其中, WB&CL&LF/HF表现出了最高的Acc和Spe, 分别是89.66%和92.59%, WB&PW表现出了最高的Sen和AUC, 分别是100%和90.28%.

指标	TP	TN	FP	FN	Acc	Sen	Spe	AUC
WB&CL&LF/HF	25	27	4	2	89.66	92.59	87.1	89.85
WB&PW&CL&HRVTi&LF/HF	24	27	5	2	87.93	92.31	84.38	88.35
WB&PW&CL&LF/HF	24	27	5	2	87.93	92.31	84.38	88.35
WB&PW	22	29	7	0	87.93	100	80.56	90.28
WB&CL&HRVTi&LF/HF	25	25	4	4	86.20	86.21	86.21	86.21
WB&PW&HRVTi&LF/HF	24	26	5	3	86.20	88.89	83.87	86.38
WB&CL&HRVTi	25	25	4	4	86.20	86.21	86.21	86.21
WB&PW&HRVTi	24	26	5	3	86.20	88.89	83.87	86.38

表5不同特征组合的CHF识别性能对比
Table5.Performance comparisons of different indices for CHF recognition.

2 -->

3.2.应用于健康人和AF患者的HRV

对如下2组HRV信号进行分析: AF组采用ltaf数据集的前43例患者, 为了和AF组形成公平对比, 并和上一组实验进行区分, NSR2组采用nsr数据集的18例健康人和nsr2数据集的后25例健康人, 共43个实验对象.
3 -->

3.2.1.NSR2和AF的HRV信号差异性分析

NSR2组和AF组的HRV信号时域、频域和ICBN指标如表6所列. 从NSR2和AF这2类人群的对比可知, WB, PW, SDNN, pNN50, RMSSD和LF/HF具有极显著统计学差异 (p < 0.001);CL, HRVTi, TP, LF和HF指标具有显著统计学差异(p < 0.01); TE和 MC指标具有统计学差异($p < 0.0{\rm{5}}$). 其中, 耦合特征中WB和PW具有极显著的统计学差异(p < 0.001); 时频域指标中SDNN, pNN50, RMSSD和LF/HF具有极显著的统计学差异($p < 0.001$).

指标		NSR2(mean$ \pm $SD)	AF(mean$ \pm $SD)	标准误差差值的95%置信区间		$p$
				下限	上限
ICBN	WB	21.483$ \pm $1.367	24.243$ \pm $3.105	1.731	3.789	0***
	PW	0.567$ \pm $0.074	0.454$ \pm $0.090	–0.148	–0.077	0***
	TE	0.941$ \pm $0.050	0.960$ \pm $0.010	0.004	0.035	0.013*
	CL	1.113$ \pm $0.146	0.999$ \pm $0.170	–0.183	–0.047	0.001**
	MC	0.687$ \pm $0.043	0.664$ \pm $0.032	–0.024	–0.008	0.04*
时域	SDNN	74.698$ \pm $26.193	139.016$ \pm $62.480	10.331	43.773	0***
	pNN50	10.123$ \pm $9.610	45.495$ \pm $30.687	4.904	25.620	0***
	RMSSD	36.402$ \pm $19.003	170.926$ \pm $97.980	15.220	104.256	0***
	HRVTi	7.735$ \pm $3.210	6.049$ \pm $2.488	–2.918	–0.455	0.008**
频域	TP	0.615$ \pm $0.612	13.493$ \pm $19.369	7.000	18.754	0.001**
	VLF	0.0002$ \pm $0.0003	0.017$ \pm $0.063	–2.256	36.147	0.083
	LF	0.157$ \pm $0.154	1.204$ \pm $1.620	0.553	1.540	0.002**
	HF	0.458$ \pm $0.515	12.272$ \pm $17.908	6.381	17.247	0.001**
	LF/HF	0.515$ \pm $0.419	0.126$ \pm $0.059	–0.519	–0.259	0***
注: *, **, ***分别代表$p < 0.05$, $p < 0.01$, $p < 0.001$.

表6NSR2和AF患者在不同分析方法下的结果
Table6.Statistical analysis results of HRV index under different analysis methods

图5是具有极显著差异特征的均值和标准差. 从图5中可以看出, 具有极显著差异的指标中, 在NSR2组中PW和LF/HF的均值比AF组的均值高, WB, SDNN, pNN50和RMSSD的均值比AF组的均值低. ICBN方法中的指标WB和PW, 时域指标SDNN, pNN50和RMSSD在AF对象中表现出了更高的标准差, 频域指标LF/HF在NSR2对象中表现出了更高的标准差.
图 5 NSR2和AF这2组对象具有极显著差异的指标的均值与标准差
Figure5. Mean and variance of the indicators with very significant differences between the two groups of NSR2 and AF.

3 -->

3.2.2.AF患者的识别

通过Fisher判别模型对NSR2和AF这2组HRV信号进行分类, 实现对不同指标分类性能的验证. 实验结果如表7所示, 通过单独的每个特征进行分类时, WB, PW, SDNN, pNN50, RMSSD和LF/HF这6个特征之中, PW表现出了最高的Acc, Sen和AUC, 分别是83.72%, 96.77%和86.57%, LF/HF表现出了最高的Spe, 是96.00%.

指标	TP	TN	FP	FN	Acc	Sen	Spe	AUC
WB	30	37	13	6	77.91	83.33	74.00	78.67
PW	30	42	13	1	83.72	96.77	76.36	86.57
SDNN	31	34	12	9	75.58	77.50	73.91	75.71
pNN50	28	39	15	4	77.91	87.50	72.22	79.86
RMSSD	29	33	14	10	72.09	74.36	70.21	72.29
LFHF	42	24	1	19	76.74	68.85	96.00	82.43

表7不同特征的AF识别性能对比
Table7.Performance comparisons of different indices for AF recognition.

本文通过穷举法, 将具有极显著差异的6个指标的所有2, 3, 4, 5, 6个指标的57种组合作为特征, 通过Fisher判别模型对NSR2和AF病人进行分类, 将分类正确率最高的5种组合列出(表8). 其中, WB&PW&pNN50&RMSSD&LFHF表现出了最高的Acc, Spe和AUC, 分别是91.86%, 89.13%和92.07%, WB&PW&SDNN&pNN50&RMSSD表现出了最高的Sen, 为97.3%.

指标	TP	TN	FP	FN	Acc	Sen	Spe	AUC
WB&PW&pNN50&RMSSD&LFHF	38	41	5	2	91.86	95.00	89.13	92.07
WB&PW&SDNN&LFHF	38	40	5	3	90.70	92.68	88.89	90.79
WB&PW&RMSSD&LFHF	38	40	5	3	90.70	92.68	88.89	90.79
WB&PW&SDNN&RMSSD	37	41	6	2	90.70	94.87	87.23	91.05
WB&PW&SDNN&pNN50&RMSSD	36	42	7	1	90.70	97.30	85.71	91.51

表8不同特征的AF识别性能对比
Table8.Performance comparisons of different indices for AF recognition.

本文首先研究比较了NSR1组对象和CHF患者HRV信号时域、频域指标和ICBN指标之间的差异, 然后又对另一NSR2组对象和AF患者进行了同样的分析.
从表3可知, 对于CHF和NSR1这2组实验对象, 所有的时域和频域指标中, HRVTi和LF/HF表现出了极显著的统计学差异($p < 0.001$). CHF与NSR1相比, 时域指标HRVTi显著减小, 表明充血性心力衰竭病人心脏活动发生变化, 交感神经活动增强, 副交感神经活动下降, 频域指标LF/HF的下降反应了交感神经和副交感神经系统活动的均衡状态发生变化, 这与文献[4,6]的研究结论相一致, 因此本文选择HRVTi和LF/HF作为时域和频域的识别指标.
从图4可知, ICBN分析方法中, CHF与NSR1相比, WB指标均值增大, PW指标和CL指标均值减小. WB指标反映的是不同层级模态分量的波动特性的加权值, CHF患者的HRV信号波动特性相对NSR1组的健康人更强, 信号频率成分更复杂, 所以CHF患者的WB均值更大. PW指标反映的是模态分量的波动特性在整个网络中的分布情况, 而NSR1组对象的HRV信号频率成分更简单, 根据LPHVG构建的网络节点的度值较小, CHF患者的HRV信号频率成分复杂, 构建的网络边数更多, 节点的度值更大, 所以CHF患者具有更小的PW值. 而CL值反映的是构建的网络的所有节点之间最短路径长度的平均值, 是网络的全局特征, 因为CHF患者的网络具有更多的边, 所有节点之间有更多的路径和更短的路径, 所以CHF组具有更小的CL值. ICBN方法是对HRV在不同时频尺度空间下的非线性特征之间的关系进行分析, 反映了充血性心力衰竭病人交感神经和副交感神经之间复杂的耦合关系. ICBN方法中CHF对象的WB指标、PW指标和CL指标标准差均大于NSR1组, 因标准差是对数据离散程度的度量, 而ICBN方法是通过度量每个IMF的复杂性来实现对HRV信号的度量, WB指标、PW指标和CL指标标准差越大表示CHF组病人的HRV信号之间差异大, 自主神经系统的调控越活跃, 而NSR1组实验对象之间的差异则较小, 实验对象的自主神经系统活动更有规律. 本文的研究初步表明, ICBN方法能够对CHF病人和健康人的HRV信号进行度量并进行区分, 但其背后的生理学机制有待深入研究.
从表4和表5可知, 在ICBN方法中, WB, PW和CL相对于时域指标HRVTi和频域指标LF/HF具有更高的识别正确率, WB, CL和LF/HF这3个指标组合的情况下具有最高的分类正确率. 实验结果表明, ICBN方法提高了分类的正确率, WB, PW和CL指标, 对于CHF病人的筛选具有相对于时域HRVTi指标和频域LF/HF指标更好的分类效果, 通过表5还可以看出不是指标越多分类效果越好. 因为指标之间并不相互独立, 过多的指标会导致分类器出现过拟合现象, 使模型不稳定, 分类性能下降. 分类正确率最高的组合中均包含本文提出的ICBN方法中的指标, 表明了本文提出的ICBN方法对HRV信号的度量具有一定的应用价值.
从表6可知, 对于AF和NSR2这2组实验对象, 所有的时域和频域指标中, SDNN, pNN50, RMSSD和LF/HF表现出了极显著的统计学差异($p < 0.001$), 这与文献[24]的研究结论相一致, 所以本文选择SDNN, pNN50, RMSSD和LF/HF作为时域和频域指标.
从图5可知, ICBN方法中, AF与NSR2相比, WB指标均值增大, PW指标均值减小, 与CHF患者和NSR1的情况相同. WB指标反映的是不同层级模态分量的波动特性的加权值, 因疾病原因, AF患者的HRV信号波动特性更强, 信号频率成分更复杂, 所以AF患者的WB均值更大. PW指标是网络单位度值下的点权值, 而NSR2组的HRV信号频率成分相对AF组简单, 根据LPHVG构建的网络节点的度值小, AF患者的HRV信号频率成分复杂, 构建的网络边数更多, 网络中度值更大, 所以AF患者具有更小的PW值. AF组的WB和PW指标的标准差均大于NSR2组, 因标准差是对数据离散程度的度量, 而ICBN方法是通过度量每个IMF的复杂性来实现对HRV信号的度量的, WB和PW指标标准差越大表示AF组病人的HRV信号之间差异越大, 而NSR2组实验对象之间的差异则较小, HRV信号更有规律.
从表7和表8可知, 在ICBN方法中, WB和PW指标相对于时域指标SDNN, pNN50和RMSSD, 频域指标LF/HF具有更高的识别正确率. 对于AF病人的筛选, 在组合所有极显著差异指标后, WB, PW, pNN50, RMSSD和LFHF组合表现出了最高的分类正确率, 且分类正确率最高的组合中均包含本文提出的ICBN方法中的指标. 实验结果表明, ICBN方法提高了AF病人分类的正确率, 对HRV信号的度量具有一定的应用价值. 还可以看出, 因为指标之间并非完全独立, 过多的指标组合会导致分类器出现过拟合现象, 反而会降低分类正确率.
ICBN分析方法首先通过ICEEMDAN方法对HRV信号进行分解, 因为健康人和病人(CHF与AF患者)的HRV信号在频率成分上存在一定的差异, ICEEMDAN将分解出不同数量和波动特性的IMFs, 因HRV信号内在的生理信息不能单纯地通过信号分解进行分离, 因此在不同IMFs之间存在着耦合的生理信息, 所以通过bubble熵对IMFs的无序特性进行度量, 然后将得到的bubble熵值通过LPHVG方法构建成网络, 通过对网络特征的分析实现不同IMF之间耦合生理信息的分析. 因为健康人和病人将构建出具有不同网络特征的网络, 通过网络特征WB, PW, TP, CL和MC对网络进行度量, 实现对HRV信号的度量. ICBN方法中包含了5个指标, 即WB, PW, CL, MC和TE. WB指标反映的是不同层级IMF的波动特性的加权值, CHF和AF患者的HRV信号频率成分相对于健康人更复杂, 因此AF和CHF患者的WB值更大. PW指标反映的是网络的单位度值下的点权值, 健康人的HRV信号规律性强, 频率成分简单, 得到的IMF数量少, 因此LPHVG构建的网络、节点和边的数量少, 节点的度值小. 相反, CHF和AF患者的HRV信号频率成分复杂, 构建的网络具有更大的度值, 因此具有更小的PW值. CL值反映的是构建的网络的所有节点之间路径长度的平均值, 是网络的全局特征, 如果构建的网络的节点之间有更多的路径和更短的路径, 则具有更小的CL值. MC指标反应的是网络的平均聚集程度, 在节点数相差不大的情况下, 如果构建的网络具有更多的边, 则聚集程度更高. TE反映的是网络的序特征, 网络的结构越规则, TE值越小. WB和PW这2个网络特征度量指标的计算过程是对HRV信号不同时频尺度下IMFs的bubble熵值直接进行应用, 而CL, MC和TE更多的是对网络结构特征的度量, 对模态分量的bubble熵值是一种间接的应用, 所以WB和PW相对于CL, MC和TE具有更好的分类效果.

本文首先将提出的ICBN方法和传统的时域、频域分析方法分别应用于NSR1组和CHF组实验对象, 结果表明, ICBN方法相对于时域和频域指标具有更高的识别正确率, 并且本文提出的WB和CL指标与频域指标LF/HF作为特征向量, 通过Fisher判别方法对CHF病人识别正确率达89.66%. 然后将ICBN方法、时域和频域分析方法应用于NSR2组和AF组实验对象, 结果表明ICBN方法中的指标相对于时域和频域指标具有更高的识别正确率. 将时域指标pNN50, RMSSD, 频域指标LF/HF, WB和PW指标作为特征向量, Fisher判别方法对AF患者的识别正确率达到91.86%. 综合以上实验结果可知, 本文提出的ICBN方法具有一定的心率变异性度量能力. ICBN分析方法在不同时频尺度空间下对非线性和非平稳的HRV的复杂波动特性进行分析, 是评价交感神经和副交感神经系统调节的潜在指标. ICBN分析方法的提出, 为心脏自主调节的复杂波动分析和HRV信号的度量分析提供了新的思路.