心电相空间映射定量分析的一种新方法——棋盘映射

　　【摘要】 目的 鉴于传统的标量坐标延时相空间重构方法存在运算复杂、耗时大、稳定性差等缺点，本文提出了相空间定量分析的一种新方法——“棋盘”映射(tessellated map, TM)。方法 首先以逐拍心电波形为列矢量构造二维矩阵X，再通过KL变换得到X的主分量，然后由X的最大特征向量重构X。基于重构误差，从X中筛选出高SNR的正常心拍，得到二维矩阵Y。由模值归一化的Y及其衍生差分阵Z，得到“棋盘格”映射图。通过空白区域的比例、点密度、棋盘映射“位图”3个指标对TM进行定量评估。 结果 实验结果表明，TM能够对精神应激和非应激进行区分，能够判断心律失常，能够用于身份识别。结论 TM能够有效地从ECG中提取信息，值得深入研究。

　　【关键词】 心电图;相空间重构; KL变换; 棋盘映射

　　A Novel Quantitative Analysis Method Called Tessellated Map for ECG Mapped Phase Space.LI Yanjun, YAN Hong, WANG Zengli, YANG Xianglin.Space Medicine & Medical Engineering,2010,23(1):57～62

　　Abstract: Objective The traditional delay time phase space reconstruction method has the disadvantages of complicated operation, time consuming and unstableness. This paper presented a novel approach called Tessellated Map (TM) for ECG phase space quantitative analysis. Methods Firstly, 2D matrix X with each ECGcycle as a column vector was constructed. Secondly, the principal components of X based on KarhunenLoeve Transform (KLT) were gotten to reconstruct X from the first eigenvector of XXT. Then the high SNR ECGcycle form X based on reconstruction error was selected, to construct the 2D matrix Y. Afterwards, TM was gotten from united Y and Z (derived from the difference of Y).Three TMbased quantitative indexes, consisted of ratio of blank zone, density of points and TM Bitmap (TMB) were gotten. Results This approach could differentiate mental stress from nonstress, estimate arrhythmia,and could be used for personal identification. Conclusion TM is a useful way to obtain information from ECG, and has its value for further research.

　　Key words: ECG; phase space reconstruction; KarhunenLoeve transform; tessellated map

　　Address reprint requests to:YAN Hong.State Key Lab of Space Medicine Fundamentals and Application, China Astronaut Research and Training Center, Beijing 100094, China

　　对心脏混沌特性的研究是当今生物信息领域的研究热点之一，而ECG是了解心脏状态的最简易手段之一。对ECG混沌特性的研究，集中于ECG节律与波形分析两个方面，ECG相空间也分为基于RR间期的心率相空间与基于心电波形的波形相空间两类。

　　1 原理

　　1.1 标量坐标延时相空间重构

　　对ECG混沌特性的定量度量，最常用的指标是李氏指数(Lyapunov exponent, LE)，其经典算法是由Packard等[1]提出的由单变量时间序列通过延时坐标映射来重构相空间的方法。此方法需要给出嵌入维数m和延迟时间τ。τ有经验值、自相关法、互信息法等几种算法，m通常大于互关联维数即可。嵌入维m和延时τ与信号采样率Fs、导联[2]、ECG噪声水平、个体差异等诸多因素相联系。王振洲等[2]指出LE和导联有关联，噪声干扰可能是不同导联ECG的LE差别的原因之一。由于嵌入维m和延时τ选择的多样性，引起的差别更大。此外，受噪声干扰影响较小的Eckmann[3]算法的运算负荷非常重，当前文献报道的测试结果通常由1 min甚至10 s左右数据得到，因此，对长程ECG数据的LE指标的实时计算更加困难。

　　1.2 两导联心电信号的正交合成

　　黄玲等[4]探讨了一种融合两导联心电信号的心电正交合成图法，并基于QRS环最大向量的幅度与角度指标进行了定性和定量分析。其研究表明[4]，由各类病例心拍构成的正交合成图的形态较正常图的差别明显[4]。

　　1.3 单导联心电的“棋盘映射”

　　由于本文提出的相空间定量分析方法以类似“棋盘”的方格为基础，故本文称其为“棋盘”映射(tessellated mapping, TM)。概括讲，先以心电幅度及其差分值来重构相空间，进而计算“棋盘”映射的定量指标。该方法简单、快捷。本文实验表明，该方法能够从ECG中提取有用信息。

　　将ECG幅度归一化得到Am(图1c)，计算Am的差分(图1d)，将差分值归一化得到Diff_A，将(Am, Diff_A)映射到二维图上(图1a)，本文称之为基于幅度及幅度差分的相空间重构法。图1(a)有3个环，最大环与最小环是QRS波群构成的，中等环是由STT段构成的。如图1(b)，由于P波幅度一般远低于QRS波与T波，因此P环通常隐匿于原点附近而不易分辨。

　　2 方法

　　2.1 算法流程

　　1)心电信号S0除噪，得到S1;2)由S1构造心拍矩阵X;3)通过KL变换，从X中筛选出高信噪比的窦性心拍矩阵Y′;4)对Y′进行幅度归一化，得到矩阵Y;5)由Y构造其差分矩阵Z′，对Z′进行幅度归一化，得到矩阵Z;6)由矩阵Y与Z进行二维相空间重构，进行方格划分并进行定量计算。

　　2.2 心拍矩阵的构造

　　通过低通滤波去除ECG的高频干扰，用分段多项式拟合基漂的方法逐拍抑制基漂，然后定位R波，进而得到RR间期序列，从R波峰点向两侧扩展，提取长度为0.8*RR的各个心拍(其中可能包含极少数错检和漏检的心拍)。

　　每个心拍x有M个点(M=Fs*0.8*RR)，N个心拍。数据矩阵的构造如下：

　　X[M][N]=[x1 x2 … xN](1)

　　2.3 基于KL变换的心拍筛选

　　1) 构造相关方阵C：

　　C=E{XXT}(2)

　　2) 基于方阵C计算特征值{λk}及其对应的特征向量{vk}，特征值λk代表了信号在特征向量vk上的能量投影，选取前K个最大特征值λk及其对应的特征向量vk。

　　3) 通过内积运算得到KL变换系数a，a是标量。每个心拍以K个KL系数a来表达。

　　an,k= (3)

　　其中k=1,2,…K; n=1,2,...N; an,k是第n个心拍在第k个特征向量上的投影强度。

　　4) 用变换系数(公式 4)重构x∧n(对应于原xn)，由于舍弃了以噪声为主的次要信息，信号的信噪比通常得到改善。

　　x∧n=∑Kk=1vkan,k(4)

　　5) 第n个心拍的重构误差为：

　　Bstd(n)=x∧n-xn(5)

　　一般而言，Bstd(n)值较高的心拍有3类：心律失常心拍、低SNR的正常心拍、R波定位错误的“伪心拍”，从X中剔除这些心拍，得到只包 图1 ECG相空间重构图

　　Fig.1 Phase space map of ECG

　　(a) TM of one beat; (b) TM of P wave; (c) ECG waveform; (d) difference of ECG waveform

　　含高SNR的正常心拍集Y。对Y的每个列矢量进行幅度模值的归一化。

　　2.4 棋盘映射及其量化指标

　　2.4.1 差分矩阵

　　对Y的每个列矢量求差分，构成矩阵Z(Z的每个列矢量进行了模值归一化)。差分方法基于5点求得(公式 6)，其抗干扰性较强。

　　y′i=-2i-2-yi-1+yi+1+2yi+2(6)

　　2.4.2 棋盘映射

　　把向量(Y,Z)逐拍映射到“棋盘格”上，如图3，其横坐标为归一化的幅度，纵坐标代表归一化的幅度差分。ECG幅度的归一化矩阵Y与幅度差分的归一化矩阵Z的取值范围都为[-1,1]。把正方形的边长进行M等分后，则得到M2个正方形(图3);把(Y,Z)逐拍映射到“棋盘格”上并进行量化分析。

　　2.4.3 定量指标

　　设M2个正方形中，有K个方格内有点，其它M2K个方格中没有点。有点的方格，按点数目从多到少依次为N1>N2>…>NK，点数目的总数为Nt。本文提出以下两个标量指标：

　　空区比重：

　　Rt=(M2-K)/ M2(7)

　　点密度：

　　Di=Ni /Nt(8)

　　其中1≤i≤K。

　　由图5可以直观地看到，心拍波形的一致性越好，则空区比重Rt越大;反之，心拍波形的差异性越大，则Rt越小。

　　本文还提出了一个矢量指标，即棋盘位图(tessellated mapping bitmap, TMB)。对含有M2个单元的“棋盘格”Map，可建立与之对应的数组bool IsEmpty[M][M]，如果Map[i][j]中有象点，则令IsEmpty[i][j]=0，否则IsEmpty[i][j]=1;这样，IsEmpty就标记了“棋盘格”中的象点分布。IsEmpty可尝试用于心律失常的判定、基于心电的身份识别等。二值位图中，‘0’代表黑色，‘1’代表白色，因此由矢量IsEmpty绘制出的二值位图中，黑色像素代表“非空”棋盘格。

　　棋盘位图刻画了棋盘映射的形态细节，而标量指标则为棋盘映射的宏观量。

　　3 讨论

　　3.1 KLT的心拍筛选功能

　　以MIT心律失常数据库(http: //www.physionet.org/ ) 的100号、101号样本的MLII导联为例(图2结果都是基于12 s数据量的运算，KLT特征向量只选取协方差阵的最大特征向量)。图2(a)室性期前收缩‘V’的重构误差是正常窦性心拍‘N’的8倍以上，图2(c)房性期前收缩‘A’的重构误差是正常窦性心拍‘N’的2倍左右，图2(b)受噪声严重影响的未分类心拍‘Q’的重构误差是正常窦性心拍‘N’的10倍以上。可见，当矩阵X以‘N’型心拍为主，即心律失常心拍、受噪声严重污染心拍相对较少时，基于KLT主分量的心拍重构误差序列Bstd(n)中，较高Bstd值的主要是心律失常心拍、低SNR的正常心拍以及误检的“伪心拍”，从二维矩阵X中剔除Bstd值较高的心拍列矢量，仅保留Bstd值较低的心拍列矢量以构成新的二维矩阵Y，则基于Y重构的ECG相空间有更高的一致性。

　　3.2 心拍筛选对相空间重构结果的改善

　　由于体表心电图都不同程度地受到噪声干扰的影响，ECG预处理只能一定程度地抑制噪声。因此，使用所有心拍重构出的相空间有较强的随机性(图3a)。以高SNR的窦性心拍重构相空间时，映射效果有了明显的改善(图3b)。

　　3.3 方格参数对TM的影响

　　映射象点的横、纵坐标都在[-1,1]的范围内。方格的边长a变化时，得到的空区比重Rt和点密度Di都有所不同。显然，a越小，则“粗粒化”的精度就越高，Rt就相应越大，Di则相应越小，而且TMB的精度也越高。图4直观反映了所述趋势，方格尺度a改变时，TM参数虽然差异较大，但数据变化趋势基本一致。

　　3.4 TM与两导心电合成向量环的比较

　　两者的相同之处在于都降低了相空间重构方法的复杂度，且方法简单、稳定、可靠;都是从一维心电向两维图像的演化，且都可向三维及更高维欧氏空间推广。

　　不同点有：TM的最初定义是基于单导心电，而心电向量环的最初定义是基于两导心电;心电向量环的量化指标是两个标量，即向量的长度与角度，而TM同时定义了标量指标与矢量指标，通过TM逻辑位图对TM映射的细节进行刻画。因此，本文提出的“棋盘格”相空间映射，实际上刻画了ECG各个心拍的一致性程度。通过KL变换剔除异常心拍、错检心拍，使得后续分析真实可信;对逐拍ECG幅度和幅度差分的归一化，使得“棋盘格”定量分析更简单快捷。

　　4 应用

　　4.1 精神应激与非应激的区分

　　采用牛津仪器(Oxford Instrument)的Medilog AR12(5个电极,3导联)采集9位受试者的数据，Medilog AR12为16位AD转化、1 024 Hz采集率。ECG采集过程为：受试者先休息15 min，再施加30 min精神应激(高声朗读英文报刊)，最后休息15 min。

　　以15 min分段处理，得到样本No1的结果(图5)。显然，较之非应激，应激状态ECG波形的一致性减小。比较休息状态(图5a,d)与精神应激状态(图5b,c)，易发现休息状态与应激状态的主要差别在中环(STT段)上，应激状态时STT段环的弥散度更高;其次，应激状态QRS环的弥散度也比休息状态高。

　　以5 min分段，计算9位受试者ECG的TM空区比值Rt并进行方差检验，结果如图6。图6的一级横坐标代表样本序号(a～i)，次级横坐标为休息bW、应激W、休息pW 3个阶段。从图6可见，休息时的Rt一般比应激高，且多数样本存在显著性差异。

　　4.2 心律失常的判别

　　以MIT心律失常数据库的100号MLII导联(图7)与102号V5导联(图8)为例，发现‘N’( normal beat)与‘P’(paced beat)、‘f’( fusion of paced and normal beat)、‘V’( premature ventricular contraction )的TMB[图7(c),(f),(i)与图8(c),(f),(i),(l),M=20]有较大差异，‘N’与‘A’(atrial premature beat)的差别不大(表1)。

　　4.3 身份识别

　　该研究基于心电信号24位受试者(记为a～x)的身份识别。采用牛津仪器的Medilog AR12采集15 min静息状态的数据，以1 min分段，则每个受试者有15段数据。个体内部的TMB差别称为“组内差别”，个体之间的TMB差别称为“组间差别”。用ECG的TMB识别身份的前提是组间差别足够大，组内差别足够小，并且满足“组间的最小差别大于组内的最大差别”， 即组间最大相关系数<组内最小相关系数。如表2，用TMB的相关系数区分这24位受试者的最佳阈值Th在[0.5438, 0.5475]区间内，即大于Th则认为同组(为同一人)，小于Th则认为非同组(是不同的人)。由于r(d,v) =0.60273，即不能区分样本d与v，因此这24位受试者基于TMB(M=100)的相关系数的身份识别率为22/24=91.67%。

　　5 结论

　　相比LE计算，TM有几个明显优势：1)运算快捷。对于ECG波形的LE计算，时间消耗十分惊人，故通常LE都基于1 min时间长度以内的ECG计算，而TM则有明显的运算效率优势，ECG的时间长度几乎不受限制。2)较高的可信度。LE需要设定嵌入维数和延迟时间，较优参数的设定十分困难，不当的设定则可能造成错误，而且其计算末期还需要拟合斜率，而斜率稳定段的寻找有很大随机性和主观性，而TM指标则有较好的稳定性和一致性。

　　本文首次提出了“棋盘格”映射图思想，并提出两个标量指标——空区比重Rt和点密度Di，与一个向量指标——棋盘逻辑位图TMB。从分析结果看，Rt可以用来判别相空间的弥散程度，而TMB则刻画了相空间的分布细节，可用于应激状表2 基于TMB相关系数的身份识别 (n=24)

　　态的判断、心律失常的判定、身份识别等，且有较好的效果。本文提出的“棋盘格”映射已应用到心电波形变异性[5]的分析中。

　　【参考文献】

　　[1] Packard NH, Crutchfield JP, Farmer JD, et al. Geometry from a time series [J]. Physical Review Letters, 1980, 45(9):712716.

　　[2] WANG Zhenzhou, LI Zheng ,WEI Yixiang. A comparative study of lyapunov spectra of CAD patients’ and healthy people’s synchronous 12 Lead ECG signals [J]. Chinese Journalof Biomedical Engineering, 2004, 23(6): 504510.

　　[3] Eckmann JP, Kamphomt SO. Lyapunov exponents from time seties [J]. Physical Review A, 1986, 34:4971 4979.

　　[4] HUANG Ling, YAN Hong. Analysis and processing of twolead ECG signals based on orthogonal synthesizing [J]. Space Medicine & Medical Engineering, 2005, 18 (3): 177181.

　　[5] LI Yanjun,YAN Hong,WANG Zengli.Study of electrocardiowaveform variability[J].Science in China C,2009,39(12):11811187.