基于排列熵的室性心律失常检测算法

耿向南

南京医科大学第一附属医院，江苏 南京，210029

0 前言

室性心律失常是一种很危险的心脏病，恶性室性心律失常甚至可以使病人心脏性猝死。准确地检测出心脏病病人发生室性心律失常的时刻，对心脏病病人的病情诊断和治疗都起着很大的作用。常见室性心律失常的检测方法有相关维的分析[1]、利用近似熵[2]、多重分形分析[3]、频域分析[4]和小波分析[5]等。大多数对于室性心律失常研究的重点都放在了区分VF(室性纤颤)和VT(室性心动过速)上。但如果一个心电自动分析系统不能精确快速地定位何时出现室性心律失常，即使分类率再高，实际应用价值也不会太大。Yinhe cao等人[6]提出了利用排列熵判断心电时间序列的不规则性，从而检测出发生病变的心电信号。zhang daoming等人[7]在其基础上利用排列熵移动窗口检测室性心律失常发生的时间，获得了较好的效果。但是，因为窗宽选择固定，没有考虑到覆盖整个心律节拍，因此，检出的室性心律失常发生的时间并不十分精确，不能准确地定位到某个节拍，且文中只对室性心律失常何时出现作了研究，并未对何时结束给出说明。针对上述文献中的不足，本文提出了一种利用检出R峰作为辅助，综合排列熵和RR间期的改变作为判定规则，精确定位室性心律失常的算法。心电数据来源于MIT-BIH心律失常和恶性室性心律失常数据库。结果表明，该算法取得了很高的检出率，可以精确地定位室性心律失常发生的节拍，并且当R峰的检测存在一定程度的漏检和误检时，也能利用R峰漏检回溯很好地定位室性心律失常发生的时间。

1 方法

1.1 基于Hilbert变换[8]的R峰检测

R峰的检出是为了近似算出心律节拍的起始点和结束点，从而得到室性心律失常发生的时间。对于受噪声污染的心电信号，检测R峰需要通过一系列的变换把R峰突出。对心电信号进行预处理，滤掉基线漂移、高频噪声，降低P、T波的影响。首先对心电信号进行小波分解，根据心电信号的频率分布特点进行重建，对重建后的信号一阶差分后进行Hilbert变换，将得到的结果与预处理后的信号综合考虑，根据阈值判定R峰。

1.1.1 基本原理

对于一个时域信号x(t)，Hilbert变换的定义如下：

对信号进行Hilbert变换后，信号的拐点对应于Hilbert变换信号与横轴的零交叉点，信号在连续的两个正和负的极值点之间的零交叉点对应于Hilbert变换信号的峰值点。

1.1.2 算法实现

利用Mallat算法对心电信号进行6层小波分解，选择的是双正交样条小波（bior4.4），它由具有最小支撑范围且容易在计算机上实现的B样条函数构造。QRS波群的能量主要集中在17 Hz左右，对于采样频率为360 Hz的信号，R峰的能量主要在第三、四、五尺度上，这里结合第三、四、五尺度的心电信号作为滤波后的心电信号。这样可以滤去基线漂移、高频干扰和大部分P、T波的能量。将预处理后的心电信号进行归一化处理。

计算一阶差分：

对 进行Hilbert变换：

归一化处理：

x(n)的峰值点对应着y(n)连续的两个正和负的极值点之间的零交叉点，y(n)在连续的两个正和负的极值点之间的零交叉点对应于Hilbert变换信号的峰值点。即x(n)的峰值对应着R峰位置。综合考虑 和 ，得到：

式中v(k)表示测量噪声。系统(4)中x(k)和d(k)同样被当做某种可预测的信号，这里假设d(k)是一个随机信号，通过以下模型描述：

选取一段心电信号，以D(n)的幅值判定R峰，取得了很好的效果，见图1。

图1 以D(n)幅值判定R峰

1.2 基于排列熵的室性心律失常的检测

1.2.1 基本原理

排列熵是一种基于相空间重构的熵算法，主要是计算相空间内点的符号序列出现概率，然后以Shannon信息熵[9]的形式计算出熵值。排列熵的变化反映并放大了时间序列微小的细节变化，可以很好地检测出复杂系统的动力学突变。

对于时间序列{x(i),i=1,2……}，以延迟1将其嵌入m维相空间：X(i)=[x(i),x(i+1),……x(i+(m-1))],对于相空间中任意一点中的m个元素进行升序排列： [x(i+(j1-1)), x(i+(j2-1))≤……≤x(i+(jm-1))]，如果存在x(i+(ji1-1))=x(i+(ji2-1))的情况，就按j值的大小排序，这样，相空间中任意一点都可以得到一组序列符号(j1,j2,……jm),这是m!种排列情况中的一种，计算每一种排列出现的概率Pj，则时间序列的k种不同符号序列的排列熵就可以按照Shannon信息熵的形式定义为：

对Hp(m)用ln(m!)进行归一化处理，则

1.2.2 算法实现

排列熵的变化可以反应出动力系统混沌特性的微小改变，一段正常的心电信号的混沌特性与心律失常的心电信号混沌特性必然不一样，利用排列熵来检测出室性心律失常可以得到相当好的效果。

（1）利用本文提出的方法检出心电信号每个心律 节 拍 的 R峰 的 位 置， 记 作R(1),R(2),……R(n)。 以Q(1)=(R(1)+R(2))/2作为第一个心律节拍的起始位置，Q(2)=(R(2)+R(3))/2作为第一个心律节拍的结束位置，则Q(1)到Q(2) 间的所有采样点近似作为第一个完整的心律节拍。依次计算Q(i)=(R(i)+R(i+1))/2，i

（2）以8个心律节拍作为窗宽。从Q(1)开始，对Q(1)到Q(9)间的所有采样点利用G-P算法[10-11]计算出适合的嵌入维，进行排列熵的计算，计作PE(1)。以一个心律节拍为移动窗口步长，对Q(2)到Q(10)间的所有采样点用G-P算法计算出适合的嵌入维，进行排列熵的计算。计作PE(2)，重复该步骤，依次得到Q(3)到Q(3+N)，Q(4)到Q(12)……Q(i)到Q(i+8)间的排列熵。分别计作PE(3)、PE(4)……PE(i)。以Q(1)到V(10)这一段采样点为例，若Q(1)到Q(9)的心律都为正常心率，Q(9)到Q(10)为室性心律失常心律，那么Q(1)到Q(9)和Q(2)到Q(10)的这两段时间序列的混沌特性有所不同，计算排列熵可以得出PE(1)>PE(2)，即移动窗口中新加入了室性心率失常的采样点，排列熵会有一定程度的降低。降低的大小取决于发生室性心律失常的混乱程度和采样点数。

窗口在移动时，若舍去的一个心律节拍为正常心率，新加入的一个心律节拍若为室性心律异常节拍，则排列熵必然会降低，若舍去和新加入的心律节拍均为室性心律异常节拍，则会影响对新加入心律节拍的判断。需要避免舍去的心律节拍为室性心律失常节拍的情况，所以要对窗口的大小做一个自适应的改变。当窗口移动时，若舍弃的一个心律节拍是已检出的室性心律失常节拍，则将窗口大小改为9个节拍，即暂不舍弃室性心律节拍，新加入一个心律节拍，得到PE'(i+1)，然后将窗口调为8个节拍，计算排列熵即PE (i+1)，用PE' (i+1)取代PE(i+1)与PE(i)进行比较才能准确地反应窗口新加入的节拍的混沌特性。这样得到了PE和PE’ 交叉的序列。

（3）综合考虑RR间期与排列熵的关系，得到：

式中c为常数，本文中取15。mean(RR)为序列的平均RR间期。

当PR(i)的值＜阈值时，那么就可以认为在第Q(i+8)～Q(i+9)段的节拍发生了室性心律失常。当发生室性心律失常时的RR间期较小时，新加入的室性心律失常的节拍的采样点不足以使得排列熵有较大下降时，RR_DI(i)的取值便会较小。而发生室性心律异常时的RR间期较大时，所加入的室性心律失常的节拍的采样点足够使得排列熵有较大下降。正常心率的排列熵的波动很小，RR间期的波动也不大，则RR_DI(i)的值也不会＜阈值。

（4）漏检回溯。R峰的检测会存在一定的漏检，需要对算法进行漏检回溯。当新加入窗口节拍Q(i+8)～Q(i+9)时，若R(i+9)-R(i+8)>1.4mean(RR)，则认为Q(i+8)～Q(i+9)间可能存在R峰漏检。取Q'(i+9)=(Q(i+8)+Q(i+9))/2，计算Q(i)～Q(i+9)间的排列熵PE''与Q(i)～Q(i+8)间的排列熵PE(i)做比较，若PE''-PE<-0.1，则认为在Q(i+8)～Q(i+9)间存在漏检，并且漏检的节拍为室性心律失常节拍，发生的位置定在Q'(i+9)，但会与实际的R峰点位置有较小的出入。

2 实验和结果

为了验证本文提出的方法，从MIT-BIH心律失常及恶性室性心律失常数据库中提取信号，检测发生室性心律失常的节拍。

取其中一段包含室性心动过速（VT）的信号，取c=15，计算式PR(i)。以-0.06为阈值进行检测结果，见图2。可以看出在R峰被正确检出的情况下每个发生VT的节拍都被正确地检出。

图2 室性心律失常检测结果

R峰的检测会存在漏检，以图2的一段心电信号为例，假使发生漏检在图中标出（miss）相应位置；在进行漏检回溯时，检测到了漏检的室性心律失常节拍在图中标注“+”。检测结果，见图3。

从图3可见，在漏检回溯时将漏检的室性心律失常节拍检出，但是因为没有检出R峰具体位置的缘故，检出的室性心律失常的位置相对有些偏差。

通过阈值划分，可以很明显地区分正常心律节拍和室性心律失常节拍，与数据库中专家的注解对比，所有标注“V”的节拍全部检出。

本文采用的方法相比有关文献[7]中以固定的2000点采样点为移动窗口宽度，100点采样点为移动步长和-0.01的斜率为检出标准的方法相比，虽然计算相对较为复杂，但是获得了更高的精度，并且把室性心律失常发生的时间定位到了每个节拍，获得了理想的效果。

3 小结

精确检出室性心律失常是区分VT和VF的先决条件，针对目前对室性心律失常的研究，提出了一种利用检测R峰辅助检出室性心律失常的算法，算法首先对心电信号进行小波变换，然后对一阶差分进行Hilbert变换，可以很精软公司推出的Windows平台应用程序开发环境，除了支持Microsoft SQL Server，它还支持IBM DB2和Oracle数据库。在解决问题之前，首先确定服务器端安装的Oracle是32 bit版本还是64 bit版本，Visual studio2010是32位的，只能调用32位（X86）程序来运行。安装时需要把客户端的版本尽量与服务器的版本一致。

（2）Sql*plus（与Oracle进行交互的客户端工具）中出现中文乱码或其他字符乱码。首先查看服务端字符集[12]：一种方法查select userenv（‘language’）from dual；另一种方法查select *from NLS_DATABASE_PAEAMETERS。假设查询出的服务端字符集为AMERICAN_AMERICA.US7ASCII，此时需要设置客户端字符集，方法一见上文（1.2 Oracle客户端字符集的注册表修改）；方法二是修改环境变量，右键，选择我的电脑/属性/高级系统/高级/环境变量/新建用户变量，变量名写入NLS_LANG，变量值写入AMERICAN_AMERICA.US7ASCII。

[1]梁志刚,李坤成.医学图像存储与传输系统[J].中国医学影像技术,2003,19(3):365-366.

[2]言伟强,刘鹏程,高文清,等.PACS/RIS在放射科医生日常工作中的应用[J].医疗设备信息,2005,20(12):53-55.

[3]Steckel RJ.The current applications of PACS to radiology practice[J].Radiology,1994,190(3):50-52.

[4]罗敏,彭承琳,王小林,等.医学图像存储与传输系统的初步应用与效益分析[J].中华放射学杂志,2003,37(2):150-155.

[5]谈永奇,候道勇.基于Oracle的数据复制技术在HIS中的应用[J].医疗卫生装备,2010,31(10):66-67.

[6]滕永昌.Oracle 10g数据库系统原理[M].北京:机械工业出版社,2006:491-509.

[7]金重午,章士正.放射科RIS/PACS构建的技术应用探讨[J].中国医学计算机成像杂志,2006,12(4):286-288.

[8]蒋昆,罗二平,漆家学等.HIS与PACS软件融合方案的设计与实现[J].医疗卫生装备,2007,28(7):37-38.

[9]陆伟.基于Web技术的医生工作站调阅LIS/RIS/PACS图文报告的接口设计与实现[J].中国医疗设备,2012,27(3):38-40.

[10]辜丽川.PACS与HIS/RIS集成网点的研究与设计[J].医疗卫生装备,2005,26(11):35.

[11]景宇宁.完美世界—Visual Studio 2010新理念选析[J].科技情报开发与经济,2010,20(19):104-105.

[12]龚维荣,周顺平,万波.浅谈Oracle数据库基于索引的SQL语句优化方法[J].计算机工程与应用,2003,39(5):196-199.