心率变异性(Heart Rate Variability，简称HRV)，被定义为连续心跳间期的变化，是一种非常有意义的分析信号. 文献[1-2]指出其变化主要受到来自心脏和循环系统的自主神经调节. 准确有效地分析HRV时间序列可以帮助我们快速地诊断心脏的健康状况，例如在心律失常和心脏性猝死等方面的预测可以让我们及时作出适当的预防措施. 因此此类研究具有重要的临床应用意义.

HRV时间序列通常可以通过测量心电图(Electrocardiography，简称ECG)随时间变化的连续R波峰之间的时间差来获得. 目前，临床上针对HRV时间序列的分析主要有时域特征和频域特征的指标. 如时域特征中的间期均值(Mean of the RR Intervals，简称MEAN)、间期总体标准差(Standard Deviation of the RR Intervals，简称SDNN)、间期均值标准差(Standard Deviation of the Averages of the RR Intervals, 简称SDANN)和相邻间期差大于50 ms所占百分比(Percentage of Successive RR Intervals >50 ms，简称PNN50)等指标，以及频域特征中信号功率谱的总功率(Total Power，简称TP)、甚低频段功率(Power within Very Low Frequency Band，简称VLF)、低频段功率(Power within Low Frequency Band，简称LF)和高频段功率(Power within High Frequency Band，简称HF)等指标. 然而，基于时域和频域特征的传统分析方法主要分别基于统计学以及平稳模型的谱估计. 因此在分析HRV时间序列这种非稳态且非线性的复杂信号时，容易丢失信号的时间动态信息. 这样导致了以上分析方法在研究HRV时间序列的有用信息方面仍然十分有限. 为了解决此类问题，基于分形维度 ^[3-4]和基于相关维度^[5-6]的复杂度测量的HRV分析方法被相继提出，并推动了HRV与生理系统调节相关性的研究. 从信息领域的角度出发，近年来有学者提出了各种基于如近似熵、排列熵和样本熵等的方法^[7-10]用于研究HRV时间序列中的复杂度、随机性和规律性等特征信息，然而上述基于熵的分析方法仍然是在单一时间尺度特征上对HRV时间序列进行分析，因此并不能充分挖掘反映人体生理变化差异性的重要信息.

为了进一步区分出不同生理健康状况所反映的HRV时间序列，本文提出了一种基于互补式集合平均经验模态分解算法(Complementary Ensemble Empirical Mode Decomposition，简称CEEMD)^[11]和改进的排列熵算法(Modified Permutation Entropy，简称mPE)^[12]的新型HRV时间序列分析方法. 互补式集合平均经验模态分解算法是对传统经验模态分解算法^[13](Empirical Mode Decomposition，简称EMD)的改进，理论上可以用于分解任何类型的非平稳及非线性的数据，同时还能有效克服原始EMD算法在实际应用中所出现的模态混叠问题. 因此在本文所提的方法中，先使用互补式集合平均经验模态分解算法将HRV时间序列分解成具有原信号不同时间尺度特征的本征模态函数(Intrinsic Mode Functions，简称IMFs)，然后通过综合分析每个本征模态函数的改进排列熵来给出对不同HRV时间序列的综合评估. 使用该分析方法可以在不同的时间尺度特征上充分分析HRV时间序列的排列熵，弥补了原始排列熵方法^[14]以及改进的排列熵方法在分析HRV中时间尺度特征单一的不足. 最后基于MIT-BIH^[15]心率失常数据库的实验结果表明，与其他现有的先进方法相比，本文所提出的HRV时间序列分析方法在区分心率失常病人与无明显心率失常病人方面具有更显著的表现.

1 经典算法原理 1.1 互补式集合经验模态分解算法(CEEMD)

针对传统经验模态分解算法容易因间断信号或干扰信号而引起的模态混叠问题，文献[11]提出了一种引入了正负白噪声的互补式分解算法，相比于其他的改进算法如文献[16]中的集合平均经验模态分解算法(Ensemble Empirical Mode Decomposition，简称EEMD)，该算法在克服了信号分解过程中容易出现模态混叠问题的同时，还极大降低了因引入白噪声所导致的残差问题. 面对HRV时间序列的高度复杂性以及不稳定性，本文将使用互补式集合平均经验模态分解算法对HRV序列进行多个时间尺度的分解. 其具体的算法流程由以下步骤给出：

步骤1：加入一组正白噪声序列 $ \{ {n_i}(t)\} $ 以及所对应的负白噪声序列 $ \{ {\rm{ - }}{n_i}(t)\} $ 到原始给定的信号 $ S(t)$ 中，其中 $ i = 1,2,\cdots,N$ . 因此可表示为

$ r_i^ + (t) = S(t) + a{n_i}(t) $

(1)

和

$ r_i^ - (t) = S(t) - a{n_i}(t) \text{，}$

(2)

其中a和N分是一组所加白噪声的幅值和总数.

步骤2：使用经验模态分解算法^[13]将已加入了白噪声的信号 $ r_i^ + (t)$ 和 $ r_i^ - (t)$ 进行分解并得到两组本征模态函数序列集 $ \{ I_{ij}^ + (t)\} $ 和 $ \{ I_{ij}^ - (t)\} $ . 其中 $ I_{ij}^ + (t)$ 和 $ I_{ij}^ - (t)$ 分别表示第i个已加入正白噪声的信号所分解得到的第j个本征模态函数，以及第i个已加入负白噪声的信号所分解得到的第j个本征模态函数.

步骤3：原信号 $ S(t)$ 中的第j个真实的本征模态函数可以表征为

$ {I_j}(t) = \frac{1}{{2N}}\sum\limits_{i = 1}^N {\left[ {I_{ij}^ + (t) + I_{ij}^ - (t)} \right]} \text{。}$

(3)

步骤4：最终原信号 $ S(t)$ 可重构为

$ S(t) = \sum\limits_{j = 1}^p {{I_j}(t) + r(t)} \text{，}$

(4)

其中p为所有真实本征模态函数的总数， $ r(t)$ 是残差信号.

1.2 排列熵算法(PE)

排列熵是一种非常有效的描述信号复杂度和随机性的数学度量工具，同时具有计算简洁的特点. 对给定的信号序列 $ \{ x(n)\} $ ， $ i = 1,2,\cdots,N$ ，其对应排列熵的计算流程为

步骤1：对于时刻n，可构造m维的向量 $ { X}(n)$ ，其可表示为

$ { X}(n) = {\left[ {x(n),x(n + \tau ),\cdots ,x(n + (m - 1)\tau )} \right]^{\rm{T}}} \text{，}$

(5)

其中m可称为嵌入维数， $ \tau $ 为时延.

步骤2：得到向量 $ { X}(n)$ 后对其元素以升序的方式重新排序，使其满足

$ \begin{split}&x(n + ({j_1} - 1)\tau ) \leqslant x(n + ({j_2} - 1)\tau )\leqslant\cdots \leqslant\\ & x(n + ({j_m} - 1)\tau ). \end{split} $

(6)

其中 $ {j_k}$ ， $ k = 1,\cdots,m$ 是向量 ${ X}(n)$ 中元素的索引值. 并且当 $x(n + ({j_{k1}} - 1)\tau ) = x(n + ({j_{k2}} - 1)\tau )$ 时，则排列顺序根据 $ {j_k} $ 值的大小而定. 即需满足

$ x(n + ({j_{k1}} - 1)\tau ) \leqslant x(n + ({j_{k2}} - 1)\tau ),{j_{k1}} \leqslant {j_{k2}} \text{。}$

(7)

步骤3：因此对于不同的 $ n = 1,2,\cdots,k$ 和 $ k \leqslant m!$ ，可构造得到不同 $ { X}(n)$ 所对应的符号序列 $ { S}(n)$ ，其中 $ { S}(n)$ 可表示为

$ { S}(n) = {\left[ {{j_1},{j_2},\cdots,{j_m}} \right]^{\rm{T}}}. $

(8)

不难看出 $ { S}(n)$ 最多有m!种不同的排列组合. 同时计算每种 $ { S}(n)$ 所对应的概率并标记为 $ {P_1},{P_2},\cdots,{P_k}$ .

步骤4：最后序列 $ \{ x(n)\} $ 的排列熵可以定义为

$ {\rm{PE}}(m,\tau ) = - \sum\nolimits_v^k {{P_v}\ln {P_v}},$

(9)

同时其归一化的排列熵可以定义为

$ {\rm{P}}{{\rm{E}}_{{\rm{normalized}}}}(m,\tau ) = \frac{{{\rm{PE}}(m,\tau )}}{{\ln (m!)}} \text{。}$

(10)

2 多时间尺度特征的HRV序列排列熵

正常的心率起伏波动较小，呈现较稳定的周期性，当心率发生明显的变异时，其对应的HRV时间序列的复杂度会随之发生改变. 因此使用排列熵这一有效的计算指标可以用来研究HRV时间序列的复杂程度和随机性等信息，从而帮助本文研究不同生理健康状态下所反映的HRV时间序列的差异性.

传统的基于排列熵的HRV分析方法无法做到从多个时间尺度分析HRV时间序列，不能充分挖掘蕴藏在HRV中的重要信息特征. 因此在本文所提出的HRV分析方法中，互补式集合平均经验模态分解算法被用于将HRV时间序列自适应地分解为具有多时间尺度特征的本征模态函数IMFs. 得益于互补式集合平均经验模态分解算法在处理HRV时间序列时拥有的抗噪声、抗模态混叠的特性，其对HRV时间序列的分解结果如图1所示，分解得到的IMFs按频率从大到小分布，清晰地呈现不同的时间尺度特征，因此分析这些IMFs更有利于挖掘单一时间尺度下传统方法所不能充分得到的信息.

然而在较低频率或较大时间尺度的序列中，如IMF4到IMF10，则有很大的概率会出现等值的情况，若仍然按照原始排列熵的处理方法来忽略掉等值的状态，则可能导致无法精确地描述序列的复杂度和随机性等信息. 因此通过使用改进的排列熵算法在表征符号序列 $ S(n)$ 时，对于等值的情况如

$ x(n + ({j_{k1}} - 1)\tau ) = x(n + ({j_{k2}} - 1)\tau ),{j_{k1}} < {j_{k2}} \text{，}$

(11)

则使用相同的符号 $ {j_{k1}}$ . 在原来的基础上，与原始排列熵算法相似，在给定嵌入维数m和时延 $ \tau $ 后，可计算得到k种符号序列 $ { S}(n)$ 与每种符号序列 $ { S}(n)$ 所对应的概率，标记为 $ P_1^{'},P_2^{'},\cdots,P_k^{'}$ ，其中 $ k \leqslant {K_{{\rm{upper}}}}$ . 关于 $ {K_{{\rm{upper}}}}$ 的具体数值则与m的取值相关，其可由递归的方式求得. 本文在表1中列出了当m取值为3~7时的 $ {K_{{\rm{upper}}}}$ 值.

最后，改进的排列熵可定义为

$ {\rm{mPE}}(m,\tau ) = - \sum\nolimits_v^k {P_v^{'}\ln P_v^{'}} \text{，}$

(12)

其归一化的形式可定义为

$ {\rm{mP}}{{\rm{E}}_{{\rm{normalized}}}}(m,\tau ) = \frac{{{\rm{mPE}}(m,\tau )}}{{\ln ({K_{{\rm{upper}}}})}} \text{。}$

(13)

在本文的研究中，甚低频段的IMFs以及残差信号对于分析HRV时间序列并无明显的参考意义. 因此本文选取了分解所得IMFs中序号为1~8的前8个IMFs作为分析HRV时间序列的研究对象，并使用改进的排列熵算法计算每个IMF的排列熵. 最后将其所有排列熵的加权平均和作为新型的HRV时间序列的总体分析指标. 在下文中简称该指标为CEEMD-mPE. 同时，图2给出了所提HRV时间序列分析方法的一般性流程框架.

3 实验分析 3.1 数据库

为了验证本文所提出的HRV分析方法的有效性，本文使用了下载自PhysioNet^[17]的由美国麻省理工学院提供的心率失常数据库^[15](Massachusetts Institute of Technology-BethIsrael Hospital arrhythmia database，简称MIT-BIH). 该数据库由48个时长约30 min、采样率为360 Hz的ECG记录组成. 每个记录的心拍类型标注均由不少于两名心脏病专家的独立意见综合给出. 在AAMI^[18]标准中，所有记录中的心跳类型可分为正常心跳(Normal，简称N)、室上性异位搏动(Supraventricular Ectopic Beat，简称SVEB)、室性异位搏动(Ventricular Ectopic Beat，简称VEB)、融合心拍(Fusion，简称F)以及存在争议的未知心拍(Questioned，简称Q). 根据AAMI标准的指导意见，编号为102、104、107和217的记录因为包含了起博心跳，所以不被考虑. 同时为了便于研究对比，本文实验从其余44个记录中选择并分成两组记录，同时根据每个记录的心拍R峰的标注信息提取RR间期得到原始的HRV时间序列. 其中第1组包含了无明显心率失常的记录，每个记录的正常心拍数占比在97%以上，第2组包含了较明显心率失常的记录，每个记录正常心拍数在85%以下. 两组记录都满足争议疑问心拍数占比少于5%的条件，其中两组记录所包含的具体记录编号可分别由如下所示：

第1组记录(无明显心率失常)：100、101、103、109、111、112、113、115、117、121、122、123、212、230、234。

第2组记录(有明显心率失常)：106、119、200、201、203、208、213、221、222、223、228、232、233。

3.2 实验结果分析

本文实验通过使用单因素方差分析法(One-way analysis of variance，简称ANOVA)对各组样本的基于不同方法得到的计算指标进行了显著差异性分析. 其中原假设 $ {H_0}$ 为各组计算指标拥有同样的均值与方差，即各组不存在显著性的差异. 相反，则备择假设 $ {H_1}$ 为各组计算指标存在显著性差异. 其中ANOVA分析的两个主要的统计量为F值和p值，F值的定义为

$ F = \frac{{\sigma _{\rm O}^2}}{{\sigma _{\rm I}^2}} \text{，}$

(14)

$ \sigma _{\rm O}^2$ 和 $ \sigma _{\rm I}^2$ 分别代表组间变异和组内变异，当组间样本存在明显差异时，则通常所对应的F值远大于1. 同时在假设检验中，p值表示当拒绝原假设时发生错误的概率. 一般情况下，若p值小于0.05时，则认为发生该类错误的情况为小概率事件，因此选择接受备择假设 $ {H_1}$ 即组间存在显著性差异. 否则就无法选择拒绝原假设 $ {H_0}$ .

在本文的实验中，首先通过使用所提出的分析方法对提取自两组记录的HRV时间序列进行了分析. 参数的选择方面，本文使用了互补式集合平均经验模态分解算法分解HRV时间序列所得的序号为1~8的前8个IMFs，然后使用改进的排列熵算法计算每个IMF嵌入维数为 $ m = 5$ ，时延为 $ \tau = 5$ 的排列熵. 每个记录具体的计算结果由表2和表3给出. 同时由图3绘制了两组记录使用互补式集合平均经验模态分解算法分解后不同IMFs的mPE分布.

通过对比不难发现，由表2与表3的数据可知. 使用本文所提的分析方法，第1组记录计算所得的CEEMD-mPE指标均低于第2组记录的CEEMD-mPE指标. 两组记录的CEEMD-mPE总体平均差异为0.050 9，个体之间最大的差异为0.112. 在对比两组记录的不同IMFs间的mPE分布时发现，序号为4~7的IMFs间的mPE差异是相对较明显的，其中两组记录总体之间的平均最大差异达到0.139 4，个体间的最大差异可达0.326 6. 其具体的差异由图3可见.

最后本文比较了所提方法得到的CEEMD-mPEHRV时间序列分析指标与其他分别计算自原始排列熵算法，改进的排列熵算法以及样本熵算法的分析指标在两组记录上的显著性差异. 基于ANOVA和箱线图的分析统计结果分别由表4和图4可见，使用本文所提的分析方法在两组样本间的分析指标其F值远大于1且p值远小于0.05，即组间的CEEMD-mPE存在显著性差异，箱线图的统计结果也证实了CEEMD-mPE可以完全区分两组记录. 相反，基于其它3种方法得到的分析指标其F值均与1相近，其中基于排列熵方法的p值虽然较小于0.05，但由箱线图可见其仍然出现了较大程度的重叠与异常值，因此差异性有限且不明显. 而其余的两种方法则表明其p值均大于0.05且同样存在箱线图显示的较大程度的重叠与异常情况，所以组间的分析指标均不存在显著性差异.

4 结论

本文结合互补式集合平均经验模态分解算法与改进的排列熵算法提出了一种新型的HRV时间序列分析方法，并计算得到一种新型的用于分析HRV时间序列的指标CEEMD-mPE. 基于MIT-BIH心率失常数据库的实验结果表明该分析指标相较于由其他先进的方法计算得出的指标可以更好地区分心率失常病人，具有统计显著性的区分度. 同时，以后可作为一项重要的特征或指标，以用于挖掘更多蕴藏在HRV时间序列中的生理信息. 因此本文的研究具有重大的临床指导意义以及未来广泛的应用前景.