技术领域
本发明涉及一种心电信号分析方法，具体涉及一种基于形态学的T波交替散点图法的 心电信号分析方法。属于生物医学信号处理技术领域。
背景技术
T波交替(TWaveAlternans，TWA)是在规整的心律时，一种体现在体表心电图上的 T波形态、极性和振幅逐拍交替变化的心电变异现象。虽然T波交替的产生机理与室性心 律失常的关系尚未有一致的定论，但大量的临床和动物实验研究表明，T波交替与室性心律 失常、室速、室颤、心源性猝死等有直接的联系，微伏级T波交替已成为预测室性心律失 常的重要指标。
T波交替是一种非稳态的变异性现象。按照T波交替检测统计方法的不同，T波交替检 测被分为短时傅里叶变换法、符号变换法和非线性法。频谱分析法(SM)是短时傅里叶变 换法中最成熟一种方法，散点图法是非线性法的一种。随着数字信号处理技术和计算机技 术的快速发展，微伏级的T波交替被检出的精确度越来越高。使用数字信号处理技术检测 体表电路图中是否存在微伏级的T波交替现象，已作为预测室性心律失常、心源性猝死风 险的独立指标，并成为一项新的无创心电检测方法。
文献《利用Poincare散点图法检测t波交替的指标研究》(李斐、赵捷等，现代生物医 学进展，VOL.12，NO.14MAY2012)和《基于Poincare散点图与谱分析法的T波交替检测》 (李斐，山东师范大学，2012年)主要是根据散点图形态和T波交替幅值变化的特点，利 用散点到x+y＝0直线的距离均值作为定量检测指标来判断T波交替的存在与否。由于该方 法受T波幅值变化影响较大，对于幅值差别较大的心电，会出现误检，降低了检测结果的 准确性。
文献《心率的Poincare散点图量化指标》(生物医学工程学杂志，2000：17(4)，P433-436) 只提出了基于散点图的四个参数：长轴、短轴、面积和角度，并没有将矢量角度指数应用 于实际心电进行检测，也没验证其可靠性。
中国专利申请201210064226公开了一种心电图中的T波交替检测的方法，采用相关分 析法检测T波交替的存在与否，其缺点是相关分析法的频率分辨性较差，如当呼吸调制的 振幅较大时，对相关分析法结果的影响要比谱分析法大得多。
中国专利申请201110418860.3公开了一种基于模型的动态心电图T波交替定量分析方 法，采用T波的解析函数即维格纳分布(Winger-VilleDistribution，WVD)对心电信号进 行定量分析，而WVD容易产生交叉干扰项，加之其运算量较大，已超过了现有硬件的处 理能力。
中国专利申请200410003212.1公开了一种用于改善T波交替的测量的方法和系统，利 用T波交替是每二次搏动自动重复并已经连接到基本心脏不稳定性的心电图的ST波段或T 波的规则的或搏动对搏动的变化的特点，将某些搏动从T波交替测量中排除出去，该方法 只是提取出T波交替信号，而并没有对T波交替信号做进一步的检测，缺乏对搏动变率可 靠性的检验。
发明内容
本发明的目的是为克服上述现有技术的不足，提供一种基于形态学的T波交替散点图 法的心电信号分析方法。针对谱分析法检测T波交替的计算复杂、检测每个心拍T波交替 幅值的困难等缺点，提出一种基于形态学的T波交替散点图法的研究，为国内相关临床研 究开展提供一种新的T波交替检测有效量化指标。本发明的技术方案是通过T波窗口分析 法，对选自美国MIT/BIH标准心律失常数据库和欧洲ST-T心电数据库的128个连续心拍的 T波进行采样，利用所得采样点作一次差分散点图，用散点图的矢量角度指数对T波交替 进行定量分析，得出判断T波交替存在与否的矢量角度指数取值范围。将散点图法的矢量 角度指数和谱分析法的T波交替幅值对T波交替的检测结果进行比较，判断矢量角度指数 检测T波交替的可靠性。
为实现上述目的，本发明采用下述技术方案：
一种基于形态学的T波交替散点图法的心电信号分析方法，包括以下步骤：
(1)选取心电信号并进行心电信号预处理；
(2)获取采样点：采用T波窗口分析法，对128个心拍进行采样，每一个心拍采用T 波窗口法选定7个采样点，就得到检测数据共有128*7个采样点；
(3)绘制散点图法：将128个相邻心拍的采样点作一次差分，绘制一次差分的T波散点 图；
(4)确定矢量角度指数取值范围，以判断T波交替存在与否；
(5)绘制功率谱曲线图：利用快速傅立叶变换将采样点的幅度转化成能量谱并进行叠 加平均得到其功率谱曲线图；
(6)计算T波交替幅度值(VTWA)：
对心电信号中N个连续ST-T心动周期的一组采样点进行FFT快速傅里叶变换分析，其 功率谱的计算公式如下：
P l + 1 = 1 N Σ i = 1 m | Σ n = 0 N - 1 b i W N n l | ]]>
其中，0≤l≤127，N＝128，m＝7，bi为第i个样点的值。
设定0.46～0.49周期/心拍内的噪声视为背景噪声，其平均频率值用Sn表示；0.5周期/心拍 处的频率值用S0.5表示，T波交替的功率值记作STWA，STWA的计算公式如下：
STWA＝S0.5-Sn
规定：当STWA<0时，说明0.5周期/心拍处的频率值小于背景噪声频率值，S0.5完全被背景噪声覆盖，则T波交替幅度值VTWA＝0，判定T波交替不存在；当STWA>0，T波交替幅度值为判定T波交替存在；
通过观察0.5cycles/beat处是否存在一明显的谱线，即为T波交替的标识谱线，可判断 在心动周期内是否存在T波交替：若线谱不明显，则心跳周期中不存在T波交替；若有明 显线谱，则存在T波交替。
所述步骤(3)、(4)、(6)均通过MATLAB7.0软件实现。
所述步骤(1)中的心电信号是从美国MIT/BIH标准心律失常数据库和欧洲ST-T心电 数据库中获得。
步骤(1)中所述心电信号预处理是指去除基线漂移、工频干扰及肌电干扰的随机干扰， 获得具有明显特征的心电信号。
步骤(4)中所述判断T波交替存在与否的散点图矢量角度指数取值范围的确定，包含 以下步骤：
41)根据散点聚集区域的不同形态以及散点图的各种定量参数值，可得到所需的不同的 信息；通过计算散点图的矢量角度指数(VAI)来判断T波交替的存在与否，其计算公式为：
V A I = Σ i = 1 N | θ i - 45 | N , ]]>
式中，N＝128，xi为第i点的横坐标，θi为第i点到坐标原点的连线与横坐标之间的夹 角度，即θi＝arctan[(xi+2-xi+1)/(xi+1-xi)]；
用软件仿真时，将45度化成弧度制取为0.7854，VAI的大小表示相邻T波间的幅值差 沿45°线的散开程度；
42)规定：当矢量角度指数0.9rad≤VAI≤1.03rad时，则在心电信号中，含有T波交替； 当VAI＜0.9rad或VAI＞1.03rad时，则在心电信号中，不含T波交替。
最后，进行了仿真结果分析：首先用MATLAB7.0软件和对谱分析法和散点图法所得 矢量角度指数和T交替幅值的数据进行对比分析，检验散点图中矢量角度指数的可靠性； 再用Kappa一致性检验进行分析谱分析法与散点图法各个指标检测结果的一致性，检验散 点图中矢量角度指数的可靠性。
所述仿真结果分析，具体包含以下步骤：
a)选取美国MIT/BIH标准心律失常数据库和欧洲ST-T心电数据库的数据，并连续取 128个周期心电信号进行软件仿真检测。MIT/BIH标准心律失常数据库原始抽样频率为360 Hz，而欧洲ST-T心电数据库的频率为250Hz，故为分析数据方便，对两组数据库的数据进 行再抽样，抽样频率为200Hz；
b)用CurveFittingTool工具对仿真所得的VTWA和VAI的离散数据进行二次曲线拟合，VTWA和VAI之间关系式如：
VTWA＝4300*VAI2-7185*VAI+3001
VTWA和VAI二次曲线的拟合度为0.867。VTWA和VAI之间的相关系数，两者之间的相 关系数γ＝0.8601。由于VT波交替和VAI之间的互相关系数γ＝0.8601，这表明由散点图法 得到的VAI与非常成熟的谱分析法得到的VT波交替具有很强的相关性，同时这也说明矢 量角度VAI能够很好应用到T波交替的检测中。
本发明的优点是：
1.本发明采用矢量角度指数，即散点图每一点到原点的连线与X轴构成的角度Ai与45 度之差的和的平均值，来判断T波交替的存在与否，故该方法与T波幅值变化率有关，进 而减小了因T波幅值变化而使检测准确率下降的可能性。
2.本发明将矢量角度指数应用于实际心电进行检测，并找出了判断T波交替存在与否的 矢量角度指数合适阈值。
3.本发明采用非线性方法中的散点图法，非线性法的抗干扰能力较强，当存在异常心搏、 急剧变化、人为干扰以及“大尾巴”的生理噪声时，非线性方法的检测结果比较理想，鲁 棒性很好。
4.本发明采用用散点图的矢量角度指数对T波交替进行定量分析，计算比较简单，易于 硬件操作。
5.本发明不仅找出了判断T波交替存在与否的矢量角度指数合适阈值，还对其可靠性作 出评估。
6.具有直观形态且计算简单，由于是差值作图，抗干扰性较强，对信号没有高质量要求， 不仅可以测得相邻T波交替幅值，而且包含时域信息。
附图说明
图1为不存在和存在T波交替散点图(a)存在T波交替散点图(b)不存在T波交替的 散点图；
图2为不存在T波交替交替的功率谱图和存在T波交替交替的功率谱图(a)不存在T 波交替的功率谱图(b)存在T波交替交替的功率谱图；
图3为VTWA和VAI的曲线拟合图。
具体实施方式
下面结合附图和实施例对本发明进行进一步的阐述，应该说明的是，下述说明仅是为 了解释本发明，并不对其内容进行限定。
一种基于形态学的T波交替散点图法的心电信号分析方法，包括以下步骤：
(1)选取心电信号并进行心电信号预处理：从美国MIT/BIH标准心律失常数据库和欧 洲ST-T心电数据库中获得心电信号，进行去除基线漂移、工频干扰及肌电干扰的随机干扰， 获得具有明显特征的心电信号。用简单整系数法去掉50Hz工频干扰，利用零相位数字滤波 器去除心电信号的基线漂移以及采用bior2.2小波函数进行小波变换阈值去噪法去除肌电干 扰，进而获得具有明显特征的心电信号。
(2)获取采样点：
我们采用T波窗口分析法，以便于更准确地检测微伏TWA(MTWA)。由RR间期的大小来选择T波窗口起点及窗宽，T波起点按经验选取见表1，窗宽取
表1
RR间期(s) RR<＝0.6 0.6<RR<1.1 RR>＝1.1 T波起点(距R峰) 60ms 100ms 150ms
根据QT和RR间期的关系式，计算出窗口的终点，如公式所示:
Q T = k R R , ( k = 0.39 ± 0.05 ) ]]>
其中，Q表示Q波起点，近似在R峰前50ms处；T表示T波起点；RR表示间期。k是常 数，一般取k＝0.4。
采用T波窗口分析法，对128个心拍进行采样，每一个心拍采用T波窗口法选定7个采样点， 就得到检测数据点集合，共有128*7个采样点；本发明采用T波窗口分析法，以便于更准确地 检测微伏T波交替(MT波交替)，即选取L＝128个连续的心跳周期，每一个周期采用T波窗口法 选定一组m＝7个采样点。
在一个心跳周期内，将T波分析窗口平均分配，每两个采样数据间的间距为ID。采样点 bi由如下公式得到，公式为：
bi+1＝b1+i*ID(i＝1，2…6)
利用QT和bi+1的公式，若每一个心跳周期采样点数为7个，且连续选取L个心跳周期， 可得到7*L个采样点，进而我们得到一组信号采样点序列P＝{x1,x2,...,xi,...}，取序列中的采样 点的差值组成新的序列Q＝{x2-x1,x3-x2,...,xi-xi-1,...}。
(3)绘制散点图：将128个相邻心拍的采样点作一次差分，绘制一次差分的T波散点图； 得出(xi+2-xi+1)和(xi+1-xi)的散点图，其中i是识别心跳序数指数。根据散点聚集区域的不同 形态以及散点图的各种定量参数值，可得到所需的不同信息。
(4)确定矢量角度指数取值范围，以判断T波交替存在与否的；
矢量角度指数(VAI)的计算公式为：
V A I = Σ i = 1 N | θ i - 45 | N , ]]>
式中，N＝128，xi为第i点的横坐标，θi为第i点到坐标原点的连线与横坐标之间的夹 角度，即θi＝arctan[(xi+2-xi+1)/(xi+1-xi)]。用软件仿真时，将45度化 成弧度制取为0.7854。VAI的大小表示相邻T波间的幅值差沿45°线的散开程度。当矢量角 度指数0.9rad≤VAI≤1.03rad时，则在心电信号中，含有T波交替；当VAI＜0.9rad或VAI ＞1.03rad时，则在心电信号中，不含T波交替。如图1不存在T波交替散点图和存在T 波交替散点图(a)存在T波交替散点图(b)不存在T波交替的散点图。
(5)绘制功率谱曲线图：利用快速傅立叶变换将采样点的幅度转化成能量谱并进行叠 加平均得到其功率谱曲线图；将采样点的能量谱对应叠加法，以提高信噪比，从而能够更 加明显地显现出T波功率谱图中的T波交替。谱分析法是一种趋近于标准化的T波交替频 域检测方法，也是目前应用最广泛的方法。
(6)计算T波交替幅度值(VTWA)：
对心电信号中N个连续ST-T心动周期的一组采样点进行FFT快速傅里叶变换分析，其 功率谱的计算公式如下：
P l + 1 = 1 N Σ i = 1 m | Σ n = 0 N - 1 b i W N n l | , ]]>
其中，0≤l≤127，N＝128，m＝7，bi为第i个样点的值。
设定0.46～0.49周期/心拍内的噪声视为背景噪声，其平均频率值用Sn表示；0.5周期/心拍 处的频率值用S0.5表示，T波交替的功率值记作STWA，STWA的计算公式如下：
STWA＝S0.5-Sn
规定：当STWA<0时，说明0.5周期/心拍处的频率值小于背景噪声频率值，S0.5完全被背景噪声覆盖，则T波交替幅度值VTWA＝0，判定T波交替不存在；当STWA>0，T波交替幅度值为判定T波交替存在。
通过观察0.5cycles/beat处是否存在一明显的谱线(即为T波交替的标识谱线)，可判断 在心动周期内是否存在T波交替。若线谱不明显，则心跳周期中不存在T波交替；若有明 显线谱，则存在T波交替。以100信号的00:17:13～00:19：00时间范围内的功率谱为例， 如图2(a)在0.5cycles/beat处不存在一明显的谱线，故在心跳周期中不存在T波交替交替。 以e0105信号的00：05：20～00：07：40的谱为例，如图2(b)在0.5cycles/beat处存在 一明显的谱线，故在心跳周期中存在T波交替。
仿真结果分析：
本发明矢量角度指数(VAI)和T波交替幅值(VTWA)进行比较。在MATLAB7.0软件的 环境下，用CurveFittingTool工具对实验中VTWA和VAI的离散数据进行二次曲线拟合，VTWA和 VAI之间关系式如：
VTWA＝4300*VAI2-7185*VAI+3001
VTWA和VAI的二次曲线拟合如图3，该坐标系中纵坐标为VAI，横坐标为VTWA，其二 次曲线的拟合度为0.867。VTWA和VAI之间的相关系数，两者之间的相关系数γ＝0.8601。 由于VTWA和VAI之间的互相关系数γ＝0.8601，这表明由散点图法得到的VAI与由已经非 常成熟的谱分析法得到的VT波交替具有很强的相关性，同时这也说明矢量角度VAI能够 很好应用到T波交替的检测中。
为了从定性的角度验证对于散点图和谱分析两种方法在检验T波交替时的一致性关系， 本发明选用Kappa一致性检验。用两种方法分别对N个观察对象逐个判断其属于C类别中 的哪一类，其判断结果常以C×C列联表的形式表示，如表2：
表2
Kappa系数的计算公式为：
K＝Po-Pe/1-Pe
其中，Po＝∑Aii/N，Pe＝∑ai*bi，ai＝Ai/N，bi＝Bi/N。P0-Pe为实际一致率，1-Pe为非机 遇一致率，Aii为CC联表中的主对角线上的实际值，N表示总例数，Ai、Bi表示第i行、 第i列的边缘值，ai、bi表示第i行、第i列的边缘概率。目前，对一致性强度的参考指标 划分有细微差异，Landies和Koch将Kappa系数划分为六个区段代表不同强弱的一致性， 如表3：
表3
Kappa系数K <0 0.0～0.20 0.20～0.40 0.40～0.60 0.60～0.80 0.80～1.00 一致性程度 极差 微弱 弱 中度 高度 极强
本发明采用更为细致的六级划分来判断。由于实际资料计算的K值只是一个两本的统 计量，存在抽样误差，故计算所得K值是否来自K值为“0”的总体即离散数据的一致程度是 由于机遇造成的，这需要假设检验即u检验，检验公式为：
S e ( K ) = P 0 ( 1 - P 0 ) N ( 1 - P e ) u = K S e ( K ) ]]>
式中，u为标准正态分位数，Se(K)为K的标准误。用Kappa系数的计算原理定性判断 T波交替的存在与否，可将不存在定性为阴性，用“-”表示；将存在定义为阳性，用“+” 表示，就得到了一个两分类，故可用一个2*2联表的形式表示检测结果。通过统计学中的 Kappa一致性检验，进行分析谱分析法与散点图法各个指标检测结果的一致性，设H0两种 方法对应的结果无关，H1两种方法对应的结果有关系。谱分析法指标与各个散点图法指标 的联表如表4所示：
表4
P0＝0.8222，Pe＝0.5238，K＝0.6266，参考评价原则，K∈(0.4,0.75)一致性好，同时，μ为 9.2419，大于95％标准正态分布位数1.96，故P<0.05，拒绝H0，接受H1，可认为两种检 测方法具有一致性。
以上两种评价方法表明VAI与VTWA具有很强的相关性，矢量角度VAI能够很好应用到T 波交替的检测中。
上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明保护范围的 限制，在本发明的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的 各种修改或变形仍在本发明的保护范围以内。