侯凤贞 黄晓林 宁新宝*

1(南京大学电子科学与工程学院生物医学电子研究所，南京 210093)

2(中国药科大学基础部信息管理教研室，南京 210009)

引言

心脏病是一种常见的多发慢性疾病，近年来，随着生活环境的变更、生存节奏的加快以及社会竞争和各方面压力的加剧，心脏病的患病率逐年升高。在我国，每年有几十万人死于心脏病，全世界更有1/3人口的死亡是由心脏病引起的。因此，如何有效地检测与评价心脏的功能状况，以便对心脏疾病进行准确的预报和诊断，是目前治疗心脏疾病的一个重要研究课题［1］。

在临床诊断心脏疾病时，作为一种无创检测手段，体表心跳频率变异信号(Heart Rate Variability，HRV)分析受到很大的重视。目前对于 HRV信号的分析，常采用时域分析法、频域分析法以及非线性动力学分析法［1－2］。

时域分析方法是研究HRV的最简单易实现的方法，主要是借用统计学的方法，对时间序列进行统计性分析。研究表明，在冠心病患者，严重的心肌病患者、充血性心脏病、心肌梗塞患者及老年患者中整个HRV时域分析指标的下降预示着术后状态较差和死亡威胁的增加［1］。时域分析已经广泛的应用来比较心肌梗死的治疗效果，估计冠状动脉受阻，预测心律失常等等［1］。

频域分析法即功率谱分析法，是将一定时间内的连续RR间期序列进行离散傅里叶变换得到以频率为横坐标、功率为纵坐标的功率谱曲线的方法。频域分析法提供了能量随频率变化分布的基本信息，能够较为有效地鉴别出交感神经和迷走神经各自对 HRV 的影响［2］。

近些年来，人们越来越清楚地认识到，心脏电活动与传统研究的随机信号和周期信号都有所不同，它由多种独立因素协调控制，犹如一个远离平衡点的多输入物理系统，富含复杂的非线性成分［1］。因此，人们开始广泛应用非线性的概念和方法对心脏动力学进行研究。众多的非线性分析方法，如幂律分析、熵分析等已经广泛地用于HRV信号的分析中;许多衡量HRV非线性特性的指标，如尺度指数、样本熵值、李雅普诺夫指数等被陆续提了出来［2－5］。另外，在对 HRV信号的分析中结合一些行之有效的研究手段，如多尺度分析、符号动力学分析、时间逆转分析，也都获得了一定的进展［6－8］。

但是，对照比较充盈性心衰(congestive heart failure，CHF)患者和正常人的各指标，我们发现很难用一个单一的指标来完全区分CHF患者和健康人。鉴于单一指标所表达出来的信息具有片面性，且各指标的使用都有其优缺点，联合多个指标进行分析有可能提高CHF的诊断正确率。作为生物控制论的一个显著成果，人工神经网络的工作原理和功能特点接近于人脑，不是按照给定的程序一步步地机械执行，而是根据自身适应环境、总结规律、完成运算、识别和控制工作。目前，人工神经网络在模式识别与分类领域中已经取得了较为显著的成果［9］。因此在本研究选取了多个从 HRV信号中提取的指标，联合起来作为人工神经网络的特征参数，利用误差反向传播(back propagation，BP)算法，建立起用于心衰诊断的智能模型，拟为充盈性心衰的诊断提供了一条无创、简捷、高准确率的补充途径。

1 实验数据与方法

1.1 实验数据

复杂生理信号研究资源库 PhysioNet是由NIBIB(National Institute of Biomedical Imaging and Bioengineering)和 NIGMS(NationalInstituteof General Medical Sciences)联合赞助的开放服务资源库［10］。本研究使用的数据来源于其中的充盈性心衰chfdb和chf2db数据库，以及正常窦性心率数据库nsrdb和nsr2db。其中包括72例健康人，平均年龄:(54.6±16.0)岁，44例 CHF患者，平均年龄(55.6±11.3)岁。几个数据库中样本采样时间都在18 ～24 h之间，为长时采样。但是在计算各HRV指标时，考虑到昼夜节律的影响，数据长度都选定为20000点(约5 h的采样时间)。

1.2 常用HRV指标的计算

1.2.1 时域指标的计算

对于数据长度为N的HRV序列{RRi:1≤i≤N}，借用统计学的方法进行统计性分析性时可以得到时域分析中的长时指标，常用的指标及其计算方法如下［1］:

· 全部正常窦性RR间期均值aRR，可依下式计算得到:

·全部正常窦性RR间期标准差SDNN，可依下式计算得到:

·全程相邻RR间期之差的均方根RMSSD，可依下式计算得到:

·每5 minRR间期平均值的标准差SDANN和每5 min的RR间期标准差的平均值SDNN index:本研究中将每个长度为N的HRV序列按每段500点划分成K=N/500小段，先计算每小段内的RR间期平均值，再计算这K个平均值的标准差得到 SDANN指标;而若先计算每小段内RR间期的标准差，再计算这K个标准差的平均值则得到SDNNindex指标。

1.2.2 频域指标的计算

按照国际上的标准，将HRV序列进行离散傅里叶变换，并对得到的功率谱依据频率的不同可以划分为4个频段［2］:0～0.003 Hz的超低频段(ultra low frequency，ULF)、0.003～0.04 Hz的极低频段(very low frequency，VLF)、0.04～0.15 Hz的低频段(low frequency，LF)和 0.15～0.4 Hz的高频段(high frequency，HF)。

对于长时HRV信号的频域分析，通常用的指标有上述4个频段的功率［1～2］，在本研究中我们分别记为 ULFP(ULF Power)、VLFP(VLF Power)、LFP(LF Power)、HFP(HF Power)。此外，所有小于0.4 Hz的频率成分的总功率TP(Total Power)也是一个常用指标［1］。

1.2.3 非线性指标的计算

在本研究中，对数据长度为 N的 HRV序列{RRi:1≤i≤N}，计算如下几个非线性指标，这些指标最近被报道能较成功地区分CHF患者和健康人:

·MSE4:是对RRi的增量序列做多尺度样本熵分析时，在尺度4下的熵值［11］。其计算方法可依如下步骤进行。

步骤1:获得增量序列

步骤2:对上述增量序列进行尺度为4的粗粒化:

y序列的长度为K=(N－1)/4。

步骤3:对序列y做m维嵌入，得到矢量序列:B(m)(i)=(yi，yi+1，…，yi+m－1)，1 ≤ i≤ K － m+1

步骤4:当两个矢量 B(m)(i)和 B(m)(j)的距离小于一个设定值r时，则认为两矢量是相同的。定义n(m)i是与B(m)(i)相同的矢量个数，则m维嵌入下的矢量相同的概率为:

步骤5:将嵌入维数增加到 m+1，同法计算C(m+1)(r).

步骤6:则y序列的样本熵值，也就是MSE4指标值为

·P8:是对 RRi做高维的时间不可逆性分析，在嵌入维为8时参数Pm的值［12］。其计算方法包括以下步骤。

步骤 1:对 RRi序列做 8维嵌入，得到矢量序列:

步骤2:计算该矢量序列在各个二维平面(RRi，RRi+n)(1≤n≤7)上的投影关于该平面主对角线RRi=RRi+n的不对性测度:

式中，ΔRRn=RRi+n－ RRi，N(ΔRRn)表示值小于 0的ΔRRn的个数，N(ΔRRn≠0)则表示值不等于0的ΔRRn的个数。

步骤3:所求的非线性指标值为

· POLVAR20:这是一个符号动力学参数。其计算方法包括以下步骤［13］。

步骤4:对 RRi的增量序列 ΔRRi(1≤i≤N－1)，按照绝对值小于20 ms符号化为0，否则为1的方式，符号化仅含0、1两个符号的序列Di(1≤i≤N－1);

步骤5:对Di序列做6维嵌入，得到矢量序列:

步骤6:统计矢量“000000”在矢量序列 B(i)中出现的次数n，即可得到所求的指标值:

1.3 HRV指标的评价

一般以受试者工作特征曲线(receiver operator characteristic curve，ROC)下的面积(area under the curve，AUC)为依据来判断指标的诊断价值［14］。如果定义灵敏度(sentivity，sen)为真阳性率，即被正确诊断为异常的人数与总异常人数的比值，定义特意度(specificity，spe)为真阴性率，即被正确诊断为正常的人数与总正常人数的比值，那么以(1－spe)为横轴，以sen为纵轴，将待检验的指标值按不同的阈值用最基本的线性感知分类器分类，分别计算 sen和spe，将结果绘制于图上，就可以得到 ROC曲线。一般AUC值应在1.0和0.5之间。在 AUC＞0.5的情况下，AUC越接近于 1，说明诊断效果越好。AUC在0.5～0.7时有较低准确性，AUC在 0.7～0.9时有一定准确性，AUC在 0.9以上时有较高准确性。AUC=0.5时，说明诊断方法完全不起作用，无诊断价值［14］。

1.4 BP网络的设计

三层的 BP网络(即一个输入层、一个隐层、一个输出层)可以拟合任意的非线性函数［15］。在本研究中，为了联合时域、频域、非线性动力学分析方法，选择常用的时域指标中AUC值最大的指标，频域指标中AUC值最大的指标以及3个非线性指标MSE4、P8、POLVAR20，共同作为 BP 网络的特征向量，因此输入层神经元即为5个。网络的输出层神经元个数K由需要分类的种数决定，本研究只研究健康人和CHF患者的区别，故选定输出神经元为一个。而为了不增加网络的复杂程度，确定网络隐层为一层，隐层神经元数定为7。隐层传递函数选择非线性激活函数tansig();输出层传递函数亦选择非线性激活函数tansig()。

1.5 BP网络的训练与评估

采用 Levenberg－Marquard 算法［16］来训练如上的BP网络;同时为了提高网络的泛化能力，采用及早停止法(early stopping)来控制训练过程［17］。在这种控制模式下，网络开始训练前需将全部的116个样本按照6∶2∶2的比例随机划分为训练(Train)样本、验证(Validation)样本、测试(Test)样本3个子集。使用训练样本进行学习，同时监测验证样本的误差。验证样本的误差在网络训练的初期是逐渐减小的，而当网络出现过拟合时，其误差便开始上升，此时就应停止训练。测试样本用于比较不同的网络模型，它的误差开始增加的时间应接近验证样本。

考虑到一次训练的结果可能并不具有代表性，因此，本研究中共进行了10000次的随机划分训练、验证、测试样本的训练。在每一次的训练结束后，借助matlab神经网络工具箱中的sim函数，将全部的116个样本的特征向量输入至训练好的网络进行仿真。对于健康人样本，期望网络输出值为-1，其仿真输出值不大于0时视为判断正确;对于CHF患者，期望输出为1，其仿真输出值大于0时视为判断正确。由于每次训练时的测试样本子集都是随机确定的，可能每次的测试子集都不相同，因此，本研究中使用全样本集的识别正确率，即116个样本中判断正确的样本百分数，而不是用测试子集的识别正确率来评估神经网络的建模结果。

2 结果

对116个HRV样本数据分别进行了时域、频域、非线性动力学分析，分别计算了上述各种常用指标及各指标在诊断CHF时的AUC值，结果如表1所示。从表1中选择出AUC值最大的时域指标SDNN、频域指标 LFP、三个非线性指标 MSE4、P8、POLVAR20共同作为 BP网络的输入向量，训练10000次后，对于全样本集的识别正确率平均为86.97%，标准差为8.86%，而最大识别正确率高达99.14%，即如图1所示，116个样本中仅有编号为100的样本被误判。而从图2则可以看出仅经过21步的迭代后，验证集的误差开始呈现上升趋势，网络停止训练;测试集的误差也在第22步迭代后开始增加，说明网络已经取得了较好的泛化能力。

表1 健康人组和CHF组常用HRV指标值(均值 ±标准差)及 AUC值Tab.1 The values(mean±standard errors)of HRV indicesfor both healthy and CHF group and the corresponding AUC values

3 讨论与结论

作为对比，选择表1中AUC值最大的MSE4指标，编程计算得出:仅仅采用该指标进行 CHF组和健康人组分类时，对于实验中的116个样本的最大正确率为85.34%。

图1 最大正确率时116个样本的仿真值Fig.1 The simulation results for the 116 samples when the maximum accuracy achieved

图2 最大正确率时网络训练过程的平方误差曲线Fig.2 The mean square error curve of network training when the maximum accuracy achieved

分别考察仅采用5个频域参数 ULFP、VLFP、LFP、HFP、TP，仅采用4个有一定诊断价值的时域参数 aRR、SDNN、SDANN、SDNNindex，仅采用三个非线性参数 MSE4、P8、POLVAR20，以及仅采用表 1 中AUC 值最大的 5个参数 MSE4、VLFP、LFP、HF、TP 作为BP网络的特征参数来对上述116个样本进行分类，同样地每种试验都进行了10000次的随机划分训练、验证、测试样本的训练，并逐次进行全部样本的仿真测试，结果如表2所示。同时，t检验的结果确证，联合时域指标SDNN、频域指标LTP、及3个非线性指标时网络的识别性能显著高于仅上述采用某一类指标时的识别性能(t检验的P值均小于0.05)。

表2 采用不同特征参数的BP网络用于CHF患者和健康人分类时的分类结果Tab.2 The classification results of BP neural network with different features used in CHF and healthy groups

上述对比实验的结果表明，联合时域、频域及非线性指标作为神经网络的特征参数来区分两类人群比仅使用某一个或者某一类指标来分类，能获得更高的识别正确率。因此可认为，对于HRV信号的分析，无论是时域分析方法、频域分析方法还是非线性动力学分析方法，所得到的每一个指标都只能片面地反映心脏状态的一个方面，只有联合线性(时域分析、频域分析)和非线性方法，使用尽可能多的指标才可能全面地揭示心脏的动力学特征，从而提高心脏疾病的诊断正确率。虽然本研究中仅仅对心衰的诊断做了一定的探索，但从探索结果来看，联合HRV信号分析的线性和非线性方法，使用多参数来进行心脏疾病的诊断，极有可能取得较高的诊断正确率。因此，将进一步充实临床数据，继续探讨本方法在不同心脏疾病诊断中的应用效果。同时，考虑到临床上采集长时心电数据的不方便性，基于短时HRV信号分析的智能化诊断系统也是未来的一个研究重点。

［1］黄晓林.心率变异性的分析方法研究［D］.南京:南京大学，2008.

［2］Task Force.Heart rate variability:standard of measurement，physiological interpretation，and clinical use ［J］.European Heart Journal，1996，17(3):354 －381.

［3］Ivanov P C，Chen Zhi，Hu Kun，et al.Multiscale aspects of cardiac control［J］.Physica A，2004，344:685 － 704.

［4］Peng Chungkang，Havlin S，Stanley H E，et al.Quantification of scaling exponents and crossover phenomena in nonstationary heartbeat time series［J］.Chaos，1995，5(1):82 － 87.

［5］Kobayashi M，Musha T.1/f fluctuation of heartbeat period［J］.IEEE Transactions on Biomedical Engineering，1982，29:456 －457.

［6］Costa M，Goldberger AL，Peng Chungkang.Multiscale entropy analysis of complex physiologic time series［J］.Physical Review Letters，2002，89:068102.

［7］Li Jin，Ning Xinbao.Dynamic complexity detection in shortterm physiological series using base－scale entropy［J］.Physical Review E，2006，73:052902.

［8］Hou Fengzhen，Zhuang Jianjun，Bian Chunhua，et al.Analysis of heartbeat asymmetry based on multi－scale time irreversibility test［J］.Physica A，2010，389:754 －760.

［9］王瑞，李杰平，卢志刚，等.基于人工神经网络遥感图像分类的应用研究［J］.科技情报开发与经济，2011，21(3):108－111.

［10］Goldberger AL，Amaral LAN，Glass L，et al.PhysioBank，PhysioToolkit，and PhysioNet:components of a new research resource for complex physiologic signals ［J］. Circulation，2000，101(23):215－220.

［11］Huang Xiaolin，Ning Xinbao，Wang Xinlong.Multiscale entropy of heart beat interval increment series and its application in diagnosis of congestive heart failure ［J］. Chinese Science Bulletin，2009，54(20):3784－3789.

［12］Hou Fengzhen，Ning Xinbao，Zhuang Jianjun，et al.Highdimensionaltime irreversibility analysisofhuman interbeat intervals［J］.Medical Engineering and Physics，2011，33(5):633－637.

［13］Wessel N，Riedl M，Kurths J. Is the normalheartrate“chaotic”due to respiration?［J］CHAOS，2009，19:028508.

［14］颜虹.医学统计学［M］.北京:人民卫生出版社，2005.217－226.

［15］宁新宝，黄晓林，徐寅林，等.神经网络分析方法用于心脏病诊断的研究［J］.生物医学工程学杂志，2004，21(2):284－287.

［16］赵弘，周瑞祥，林廷圻.基于 Levenberg－Marquardt算法的神经网络监督控制［J］.西安交通大学学报，2002，36(5):523－527.

［17］戴荣，李仲奎.三维地应力场BP反分析的改进［J］.岩石力学与工程学报，2005，24(1):83－88.