符号动力学在心率变异性分析中的参数选择*

宋爱玲 黄晓林 司峻峰 宁新宝

(南京大学电子科学与工程学院，生物医学电子工程研究所，近代声学教育部重点实验室，南京 210093)

符号动力学在心率变异性分析中的参数选择*

宋爱玲 黄晓林司峻峰 宁新宝

(南京大学电子科学与工程学院，生物医学电子工程研究所，近代声学教育部重点实验室，南京 210093)

(2010年4月30日收到;2010年5月24日收到修改稿)

时间序列的符号动力学信息熵 Hk因其计算简单快速，对数据量要求小，而被应用于心率变异性(heart rate variability，HRV)分析，然而符号化的参数选择至今却并未形成统一标准.HRV作为典型的生理信号，存在着极大的个体间差异和非平稳性，要获得稳健的一致性分析，在符号化过程中必须考虑符号化参数 α与序列本身均值、标准差的综合影响.文中，首先以仿真噪声序列为对象，考察了3个参数对于Hk的影响及三者相互之间的关联性，研究表明当满足特定关系时，Hk的曲线簇收敛于反映序列动力特性的Hk-up;随后在对15例心跳间隔序列的分析中，验证了Hk-up在消除个体间差异及减弱非平稳干扰影响两方面都优于α取固定值时的研究结果.

符号动力学，熵，心率变异性

PACS:05.45.Tp，87.19.Hh

1.引 言

健康人的逐拍心率或心跳间隔不是整齐划一的，而是存在着微小的波动［1，2］，俗称心率变异性(heart rate variability，HRV).HRV的存在绝不仅是外界刺激的结果，而是由心脏作为复杂的多输入的非线性系统的本质造成，即使是健康人在安静状态下也普遍存在［3］.目前所知的健康人的心率至少受3个因素——自动窦性节律、交感神经以及副交感神经的影响［1—3］，而去神经支配后的自动窦性节律是相对稳定的，因此，HRV主要反映了自主神经系统对心脏功能的调控.以往研究表明:自主神经系统与心血管疾病以及心源性猝死有着密切联系，例如由于交感张力的升高与副交感张力的降低，使室颤阈降低，易引发室速、室颤和猝死，而HRV的降低则可以作为急性心梗后死亡的强有力的独立的预测因子［4，5］.同时，对于许多其他与自主神经功能相关的非心脏性疾病，如糖尿病、睡眠呼吸停止、癫痫等，以及运动医学等方面，HRV也具有重要的临床诊断或研究价值［6—8］.

HRV的分析方法主要包括线性方法和非线性方法两大类［1，2］.其中如频域分析等传统的线性方法至今仍被广泛应用于临床诊断.20世纪90年代以后，随着对心率波动信号的非线性本质的认识的逐步深入［9］，HRV的非线性分析成为研究的热点［2，10］.尽管HRV的非线性分析目前在临床应用中还未形成普遍的统一的标准，但是作为心脏动力系统非线性本质的刻画，非线性方法仍然显示出其独特的重要性［11］.

在HRV的非线性分析中，有一类方法——与符号动力学相结合的熵分析方法，由于其对数据长度要求小，计算简单快速，而显得极具吸引力.在该类方法中，首先将时间序列，例如心跳间隔序列，在幅度域上进行符号化(粗粒化)编码，从而将模拟量的时间序列简化为由有限个符号构成的符号序列，然后考察某一固定长度下时间序列中各种子串模式的出现概率，并对其求信息熵.幅度域的符号化一方面极大地提高了计算速度，另一方面，如果方式选择得当，尽管会丢失一些细节信息，但却能在保存动力系统的本质特性的同时大大降低噪声的影响，因此选择恰当的符号化依据成为符号化熵分析方法中的关键.符号化依据的选择一般有3种方式: 1)以整个时间序列的均值或方差为基准作为符号划分的阈值［12，13］，或直接定义临界点［14］，例如Kurths等［12，13］提出的符号动力学信息熵;2)以时间序列的局域波动作为符号划分的依据，例如李锦等［15］提出的基本尺度熵;3)以时间序列中相邻两点的差值(一阶差分)作为符号化的依据，例如卞春华等［13，16］提出的符号序列熵.后两种方法中，无论是以局域波动作为符号划分的依据，还是以时间序列的一阶差分作为编码的依据，其结果都忽视了时间序列的整体信息或低频特性［17］，例如在对高斯白噪声的基本尺度熵计算中，其熵值远远达不到理论值，而不同类别信号间的差异也因为低频信息的丢失而大大减小，不利于信号的区分.对于这两种方法的具体讨论，可参考文献［13，15—17］.本文中主要讨论的是第一种方法，即符号动力学信息熵的参数影响及其选择依据.

2.符号动力学信息熵

符号动力学信息熵最早由 Kurths等［12，13］提出，其计算方法如下:对于数据长度为N的时间序列{xi}(1≤i≤N)，将其转换成符号序列{si}(1≤i≤N)，其中si∈{0，1，2，3}，具体的转换方法如下:

式中μ代表时间序列{xi}的均值，α则为一个特殊参数.

对于符号序列{si}，其中所有长度为m的子串的集合可表示为{Wj}(1≤j≤N-m+1)，Wj=(sj，sj+1，…，sj+m-1)，考察{Wj}的分布特性，即对出自符号表{0，1，2，3}的长度为m的所有可能的4m个符号串，统计其在{Wj}中出现的次数C(l)(1≤l≤4m)，并估计出概率

对其计算信息熵

从而得到时间序列{xi}的符号动力学信息熵.Hk反映了时间序列在符号化以后各种长度为m的子串的模式丰富程度和分布特性，出现的符号串模式越多、分布越分散，则熵值越高;符号串模式越少、分布越集中，则熵值越低.例如，对于时间上无相关的完全随机序列，如高斯白噪声序列，其长度为m的子串在字符表{0，1，2，3}中取值的可能组合共有4m种，且各种组合模式出现的概率基本一致，因而计算得到的熵值达到最大值log(4)，而对于有特定时序结构的序列，其子串的模式则会远远小于4m，同时在分布上也很可能出现不均衡现象，从而熵值小于log(4m).

关于计算Hk时参数的选择，子串的长度m的增加会导致计算量及必要数据长度的指数增加，而对于结果却无本质影响，因此，一般取到3就足够了［12，13］;时间序列长度N的选择，为使得Hk具有统计学意义，一般要求N4m［12］;而符号化过程中的参数选择则非常重要.从(1)式可以看出，序列{xi}的均值μ和参数α都会影响符号化的结果，另外，尽管(1)式中没有直接反映，但是可以推断，序列{xi}的标准差(standard deviation，SD)也会影响符号化的结果.事实上，这三个参数并非完全独立，而是相互关联的.在以往研究中，Kurths和 Wessel等只给出了在研究HRV时，α可以取0.1［12］或者取0.03到0.07之间［13］，然而，HRV作为典型的生理信号，个体间差异大，且同一个体也往往存在着极大的非平稳性，固定的参数选取会导致分析结果受到心跳间隔序列的均值、标准差等参数的影响，而无法获得稳健的一致性结果.因此，在下面的章节中，将进行详细的符号化参数分析.

3.基于仿真噪声序列的符号化参数分析

本节将以仿真噪声序列为例详细讨论3个参数的影响以及相互之间的关联性.

3.1.μ，α和SD对Hk的影响

首先以高斯白噪声序列为对象，考察μ，α和SD对于Hk计算的影响.为简化起见，事先将长度N= 1000的噪声序列进行了零均值和归一化标准差的预处理，然后人为地改变其μ或者SD，并计算不同α取值时的Hk，计算中，取m=3.结果如图1所示，其中图1(a)为固定SD=1时，Hk随μ和α变化的曲线，图1(b)为固定μ=10时，Hk随SD和α变化的曲线.

图1 高斯白噪声序列Hk随α，μ，SD的变化情况 (a)固定SD时，随μ变化的Hk-α曲线簇，(b)固定μ时，随SD变化的Hk-α曲线簇

从图1中可见，Hk的值随着时间序列的 μ，SD以及参数 α的变化而变化，即不同的序列均值、标准差和参数α都会影响计算出来的Hk;另一方面，对于多数α-Hk曲线，都存在着一个极大值，不同曲线的极大值基本相等，即收敛于一个固定值，这里定义为Hk-up，该值非常接近理论推导值 log(43)，这说明:在μ，SD，α三个参数的特定组合下，Hk能达到理论值，从而能够正确的反映时间序列的内在动力学特性;而参数取值不当时，Hk则可能会远远小于log(43)，从而不能真实地反映序列的复杂程度.

除了高斯序列，我们还对1/f噪声和布朗噪声做了相同的计算.计算结果显示，对于这两种噪声序列，其Hk随着μ，SD以及α的变化情况与图1类似，并且 α-Hk曲线簇也存在着收敛的极大值Hk-up——1/f噪声为3.5，布朗噪声为1.9.1/f噪声和布朗噪声都存在着时间相关性，两者的区别在于尺度因子的不同.对于存在时间相关的序列，其长度为m的子串会集中在相对有限的模式，因而其Hk理论上应小于无时间相关性的高斯序列，且尺度因子越大，相关性越强，其 Hk则越小，这与我们的计算结果是相符的.

图2 获得Hk-up时的μ-α关系曲线 (a)高斯白噪声，(b)1/f噪声，(c)布朗噪声

3.2.最佳参数设置

既然对于三种典型噪声序列，其α-Hk曲线簇的极大值都收敛于Hk-up，并且通过高斯序列已经验证了Hk-up与理论值非常接近，我们推断Hk-up可以作为时间序列复杂性的一种稳健的反映.下面仍然以上述三种仿真噪声序列为对象，考察要获得 Hk-up，μ，SD以及α之间必须满足的关系.简单起见，这里将SD固定为1，着重研究μ-α关系.固定μ时，SD-α关系可以类似方法获得，在此不赘述.

图2中分别是2000点的标准差为1的高斯序列(a)，1/f噪声序列(b)，布朗噪声序列(c)获得Hk-up时，μ-α关系双对数曲线.从图2中可见，三种序列与Hk-up对应的μ-α关系双对数曲线都呈现出非常明显的线性，因此，针对三种噪声各150例独立的序列，进行了μ-α关系双对数曲线的线性回归并作了F检验，详细结果列于表1.

表1 基于仿真噪声序列的μ-α双对数曲线线性回归

从表1中可看出，μ在2.5—100之间时，经 F检验证实，三种噪声μ-α关系的双对数曲线的线性关系显著，且线性拟合的斜率均近似为-1，截距也都在-0.4附近.三种噪声之间，无论是拟合的斜率还是截距，差异均不显著(t检验P>0.05)，这说明:对于这三种存在不同时序结构的序列，要获得反映其动力学复杂性的 Hk-up，当 SD=1时，其 μ-α关系遵循基本同等的规律，即满足

由(4)式可以看成在计算符号动力学信息熵Hk时参数选取的最佳方案.有了这个标准，按(1)式符号化时，时间序列单纯由均值、标准差的线性变化带来的影响基本被消除.

4.HRV的分析应用

在本节中我们研究的序列来自于Chaos在2008年发起的专题讨论“正常心率是混沌的吗?”所提供的数据库*，其中包括 5例正常窦性心率(normal sinus rhythm，NSR)(年龄30±10)，5例充盈性心衰(congestive heart failure，CHF)患者(年龄59±9)，5例房颤(atrial fibrillation，AF)患者(年龄不详)的长时(20—24 h)心跳间隔序列.生理信号的个体差异和非平稳性，一直是阻碍其分析获得稳健的一致性结果的一个重要原因.例如在HRV研究中，采集心电数据时不同个体之间所处状态的差异，以及单个个体在持续心电采集时心跳间隔序列的非平稳性，都会严重影响到分析结果.本节中，针对上述两种情况，分别研究了利用最佳参数设置后，符号动力学信息熵对于HRV分析的改善.

首先对所有15例数据的最初2000点，分别利用固定参数，取 α=0.05(仿照文献［13］)，α= 0.01，以及第3节得到的最佳参数设置方法，计算了各自的符号动力学信息熵，分别记为 Hk1，Hk2和Hk-up，结果如图3所示.图3中，横轴是15例数据的编号，其中1—5为NSR(n1nn—n5nn)，6—10对应CHF患者(c1nn—c5nn)，11—15则对应 AF患者(a1nn—a5nn);纵轴是符号动力学信息熵 Hk的计算值;图中‘’代表应用最佳参数方案获得的熵值Hk-up，而‘*’和‘+’则分别代表不考虑个体间均值差异而固定取α=0.05时的熵值Hk1和 α=0.01时的熵值Hk2.

图3 心跳间隔序列的Hk1，Hk2和Hk-up结果比较

由图3可见，三组样本的 Hk-up相对于Hk2都有明显的改变，而对于 ID为4以及7—15的序列，Hk-up与 Hk1也有显著的差别，特别是，采用 Hk-up时，AF患者完全与NSR和CHF区分开来，同时AF类别内的个体差异也明显减小.为进一步说明利用Hk-up更能反映序列的动力学本质而尽可能排除均值等线性因素的影响，图4中分别给出了CHF组(a)和AF组(b)的μ—Hk的散点图，图中实线代表拟合直线.从图4中可见，Hk-up与序列均值的线性相关性相对于Hk1，Hk2减小，这说明，采用标准化参数后，组内成员之间尽管采样时存在均值、标准差等状态的不同，Hk-up仍然能真实地反映动力系统的复杂性，从而获得一致性的结果.

图4 心跳间隔序列的μ—Hk散点图 (a)CHF组，(b)AF组

接下来，则通过对长时数据的分段分析来考察在同一个体中，Hk-up减弱均值波动影响的能力.对于每例数据，我们都将其按2000点一段分为31段，分别计算每段的均值、Hk1、Hk-up，再考察每例样本Hk1，Hk-up的波动范围及其与均值之间的相关性，对于线性相关显著(P<0.05)的则对其进行线性拟合，比较线性拟合的斜率，图5画出了具有代表性的a1nn数据的分段分析结果，其中图5(a)为熵值随分段时间的波动，其横轴代表分段的序号，纵轴代表符号动力学信息熵，图5(b)则为其μ—Hk的散点图，实线代表线性拟合.由图5可见，就 a1nn而言，Hk-up的波动范围明显小于 Hk1，且 Hk-up对于 μ的线性依赖性也显著降低.三组共15例样本的全部结果列于表2中，其中，Δ代表参数的值域，k代表线性拟合的斜率，NS代表线性相关不显著(P>0.05).

从表2中可以看到，就熵值的变化范围而言，除了n3nn，c1nn，c4nn(表中有下划线的数据)以外，其他样本中Hk-up波动都是小于 Hk1的;而从 μ—Hk的相关性分析及线性拟合斜率来看，除了n3nn(表中有下划线数据)的 Hk-up随 μ的变化比 Hk1大，n4nn，n5nn，c1nn，c4nn两种熵值与μ均无显著线性相关，其他样本都表现出 Hk-up随 μ的变化减小或者与 μ的线性相关性的削弱.这两个方面都说明，采用最佳参数设置后的Hk-up对于HRV的复杂性描述更稳定，更不易受均值波动的影响.而在三组样本中，以AF组的改善效果为最好，这是由于 α取固定值0.05，即是文献［13］在完全不考虑 AF病人的前提下针对正常范围的平均心率和标准差情况挑选过的，所以就NSR，CHF而言该α与最佳参数偏离不是特别大，而AF组其均值和标准差较NSR和CHF组有较大偏离，因而α取0.05的选择严重偏离其最佳参数，从而造成看起来Hk-up对于 AF组的改善效果最好.

图5 a1nn的长时心跳间隔序列分段分析 (a)Hk随分段时间的波动情况，(b)μ—Hk散点图

表2 长时心跳间隔序列分段分析的Hk波动范围及与μ的线性相关性分析

图6为三种噪声和三组HRV的Hk-up值的比较.

由图6可见，AF组可以完全与NSR和CHF组区分开来，而且与高斯噪声接近，而NSR和CHF组的Hk-up则处于1/f噪声和布朗噪声之间.这说明，AF患者的心跳间隔序列存在着严重的时间相关性缺失，其随机性与白噪声接近，而NSR和CHF患者的心跳间隔序列则存在着时间相关性，这种相关性比1/f噪声略强，而比布朗噪声又要弱，这与以往通过功率谱或去趋势波动(detrended fluctuation analysis，DFA)等方法得到的研究结果是一致的［18，19］.

图6 仿真噪声和HRV的Hk-up对照

5.讨论和结论

以往研究表明，心率波动序列作为确定性系统产生的非随机序列，仅由概率分布函数来描述是不够的，其排列的先后顺序也蕴含了大量的内在动力学特性［1］，符号动力学信息熵则能从一定程度上衡量这种时间前后的关联性，从而评价时间序列的复杂性.符号动力学信息熵计算简单快速，对数据量要求不大，但是符号化参数的选取一直以来并未形成统一标准，而HRV作为生理信号的一种，具有个体间差异大和非平稳等特点，因此，如何尽量消除个体差异使参数具有可比性，以及如何消除非平稳干扰的影响而获得一个稳健的一致性的分析结果，使得为符号化建立一个参数选取的最佳方案成为必须解决的问题.

本文对符号动力学信息熵Hk的符号化方法进行了参数分析:首先通过三种代表性的仿真噪声序列(高斯噪声、1/f噪声、布朗噪声)，研究了符号化参数α与序列本身均值、标准差的影响及相互关联性，结果表明当 α与均值、标准差满足特定关系((4)式)时，Hk会收敛于反映真实动力学特性的Hk-up;然后，针对三组共15例心跳间隔序列验证了利用最佳参数设置方案的 Hk-up无论是消除不同个体间状态差异，还是减弱同一个体在长时心电采集中的基础心率的影响，都优于固定α的Hk.最后，通过比较三组HRV与三种噪声的Hk-up，认为AF组心跳间隔序列在时间相关性上与高斯白噪声基本接近，而NSR组和CHF组则更接近1/f噪声，即NSR和CHF组的心跳间隔序列存在着类似1/f噪声的时间相关性，这与以往研究认为正常人 HRV包含1/f成分［18，19］的观点是相符的.

同时我们也看到，在NSR和CHF组的Hk-up分段分析结果中，Hk-up仍然存在着不可忽略的波动，既然基础心率等线性非平稳干扰带来的影响已经消除，那么 Hk-up的波动很可能是由心脏动力系统本身固有的低频特性所引起，2000点的心跳间隔序列不足以包含这种低频成分，这很可能就是导致NSR组和 CHF组之间仅由 Hk-up无法区分的重要原因.要详细研究这些低频成分，需要有大量的长时心跳间隔序列，并且要在采集过程中严格标定测试者各个时间段的状态，这将是我们进一步工作的方向.

感谢南京邮电大学地理和生物信息学院马千里博士对文中噪声的获得和处理提出的宝贵建议.

［1］Camm A J，Malik M，Bigger J T，Breithardt G，Cerutti S，Cohen R J，Coumel P，Fallen E L，Kennedy H L，Kleiger R E，Lombardi F，Malliani A，Moss A J，Rottman J N，Schmidt G，Schwartz P J，Singer D H 1996 Eur.Heart J.17 354

［2］Kleiger R E，Stein P K，Bigger J T 2005 A.N.E.10 88

［3］Ivanov P C，Chen Z，Hu K，Stanley H E 2004 Physica A 344 685

［4］Kleiger R E，Miller J P，Bigger J T，Moss A J 1987 Am.J. Cardiol.59 256

［5］Malik M，Farell T，Camm A J 1990 Am.J.Cardiol.66 1049

［6］Huikuri H V，Makikallio T H，Airaksinen K E J，Seppanen T，Puukka P，Raiha I J，Sourander L B 1998 Circulation 97 2031

［7］Stein P K，Kleiger R E 1999 Annu.Rev.Med.50 249

［8］Zhuang J J，Ning X B，Zou M，Sun B，Yang X 2008 Acta Phys. Sin.57 2805(in Chinese)［庄建军、宁新宝、邹 鸣、孙 飙、杨 希2008物理学报57 2805］

［9］Li C，Tang D K，Fang Y，Sun J T，Ding G H，Poon C S，Wu G Q 2009 Acta Phys.Sin.58 1348(in Chinese)［李 程、汤大侃、方 勇、孙锦涛、丁光宏、Poon C S、吴国强2009物理学报58 1348］

［10］Ning X B，Bian C H，Wang J，Chen Y 2006 Chin.Sci.Bull. 51 385

［11］Wessel N，Riedl M，Kurths J 2009 Chaos 19 028508

［12］Kurths J，Voss A，Saparin P，Witt A，Kleiner H J，Wessel N 1995 Chaos 5 88

［13］Wessel N，Ziehmann C，Kurths J，Meyerfeldt U，Schirdewan A，Voss A 2000 Phys.Rev.E 61 733

［14］Liu X F，Yu W L 2008 Acta Phys.Sin.57 2587(in Chinese)［刘小峰、俞文莉2008物理学报57 2587］

［15］Li J，Ning X B，Wu W，Ma X F 2005 Chin.Phys.14 2428

［16］Bian C H，Ma Q L，Si J F，Wu X H，Shao J，Ning X B，Wang D J 2009 Chin.Sci.Bull.54 4610

［17］Huang X L，Cui S Z，Ning X B，Bian C H 2009 Acta Phys. Sin.58 8160(in Chinese)［黄晓林、崔胜忠、宁新宝、卞春华2009物理学报58 8160］

［18］Kobayashi M，Musha T 1982 IEEE Trans.Biomed.Eng.29 456

［19］Peng C K，Havlin S，Stanley H E，Goldberger A L 1995 Chaos 5 82

PACS:05.45.Tp，87.19.Hh

Optimum parameters setting in symbolic dynamics of heart rate variability analysis*

Song Ai-Ling Huang Xiao-LinSi Jun-Feng Ning Xin-Bao
(Key Laboratory of Modern Acoustics of Ministry of Education，Institute of Biomedical Electronic Engineering，School of Electronic Science and Engineering，Nanjing University，Nanjing 210093，China)

30 April 2010;revised manuscript

24 May 2010)

The Shannon entropy Hkin symbolic dynamics analysis of time series has less data demand and can be complemented easily，and therefore is applied to heart rate variability(HRV)analysis.However，the criterion has not been established as to the parameter α，which is used during the transformation into symbols.Like all other physiological signals，great differences between individuals as well as nonstationarity are usual in HRV.Therefore，to achieve a robust and consistent analysis，the parameter α should be considered together with the mean and standard error of the series.In this paper，the integrative effects on Hkof the three parameters were studied firstly in simulation time series，and it was suggested that under certain conditions，the clusters of Hkapproach a convergent upper boundary named Hk-up，which reflects the intrinsic dynamics of the time series.Then，in the heart beat interval series analysis of 15 subjects，the abilities of Hk-upto alleviate the impacts brought by both difference between individuals and nonstationarity were testified.

symbolic dynamics，entropy，heart rate variability

*国家自然科学基金(批准号:60701002)资助的课题.

*http://www.physionet.org/challenge/chaos/

*Project supported by the National Natural Science Foundation of China(Grant No.60701002).