基于EMD的自适应ECG信号基线漂移消除方法研究

郭树言，王黎明，牛晓东，张岩，李集照，赵越

中北大学信息探测与处理山西省重点实验室，山西太原030051

前言

医学领域在采集心电（ECG）信号时，往往会由于各种原因叠加不同类型的噪声，如高频电噪声、肌电噪声、由于病人呼吸和活动产生的基线漂移（BW）噪声，这些噪声的存在有时会严重干扰后续对ECG信号QRS 复合波的检测，而QRS 复合波形状的一些具体指标是医生诊断某些疾病的关键，因而在QRS复合波检测前对采样得到的原始ECG 信号进行滤波是必不可少的［1-3］。然而，ECG 信号是一种典型的非线性非平稳的信号，使用传统时频域信号处理方法（如短时傅里叶变换、小波变换）都有频谱展宽的问题，且原始ECG 信号叠加的各种噪声的频域范围很宽，这样就导致传统时频域信号处理方法在消除ECG信号各种噪声时失真严重或效果不理想。

近年来，EMD 方法由于对非线性非平稳信号良好的处理效果逐渐被国内外学者应用在ECG 信号去噪领域，尤其是在BW 去除方面提出了许多优秀的方法，如EMD 方法与低通滤波器相结合［4］，EMD 方法与形态学滤波相结合［5］，EMD 方法与BW 特点相结合进行经验滤波［1，6］等。上述方法都是通过EMD 分解产生ⅠMF 分量，再对单独的ⅠMF 分量进行处理或检验。在Weng 等［4］研究中，通过对高阶ⅠMF 分量进行低通滤波，得到对应ⅠMF的低频成分作为BW 的部分信号。但低通滤波器截止频率的选择没有自适应性，而且ⅠMF 分量本质上也是非线性非平稳信号，直接进行频域低通滤波而不考虑EMD方法中最重要的瞬时频率的观点，必然造成对单一的ⅠMF 分量，在有的时间范围内对信号过处理，而在有的时间范围内滤波不彻底的问题；而文献［5］提出的方法中为判断每个ⅠMF 分量是否是零均值的，引入固定的阈值，从而判断ⅠMF 分量中是否含有BW 噪声的成分；这样，由于固定阈值造成算法处理数据时自适应性的缺失不可避免。文献［6］采用与文献［5］相同的BW 去除算法，自适应性不强。本文在详细研究EMD 分解产生的ⅠMF 分量物理意义的基础上，提出一种自适应的BW 去除算法，利用EMD 自适应性和ⅠMF 频率筛选性质筛选出BW 噪声，并对算法性能进行定性和定量分析。

1 经验模式分解

1.1 EMD方法的物理意义

EMD 方法是1998年黄锷教授等人提出的，它本来是为了将非平稳非线性信号分解为一系列的固有模式函数（ⅠMF）。这样对ⅠMF 分量进行希尔伯特变换之后得到的瞬时频率才是有意义的。因为认为每个ⅠMF分量在任意时刻只有1个瞬时频率，而原始信号是所有ⅠMF 分量的叠加（将剩余分量算为1 个ⅠMF），在原始信号的任意时间点上有多个瞬时频率，所以对原始信号直接求希尔伯特变换得到的瞬时频率是没有物理意义的［7-14］。从这一意义上说，如果将原始信号理解为多个信号源产生信号的叠加，那每一个ⅠMF分量就代表了1个信号源产生的信号，虽然每一个ⅠMF 分量也是非线性非平稳的，但其瞬时频率和幅度的变化一定有一个过程，即瞬时频率和幅度的变化是连续的，而不是突变的。这样，每一个ⅠMF分量的瞬时频率总在一个范围中，而高阶ⅠMF的瞬时频率所处的范围总体上总是低于比它低阶的ⅠMF 分量的。如果组成一个信号的信号源发射信号的频率总体上的高低已知，我们可以通过ⅠMF 分量的频率范围来确定这个ⅠMF对应的信号源。

1.2 EMD过程描述

EMD 方法的核心概念是时间尺度，它被定义为信号的一个极值点到与它相邻的极值点的时间范围（如果是数字信号，定义为数字点数）。如果是原始信号，这个时间尺度就被视为信号中最高频分量（第一个ⅠMF 分量）的频率指标，而其他ⅠMF 分量合起来被视为是第一个ⅠMF 分量的趋势，因为ⅠMF 分量是零均值的，EMD 方法中直接将原始信号上下包络求均值作为趋势量，然后原始信号减去趋势量得到第一个ⅠMF 分量［7，12，15-18］。因为得到的ⅠMF 分量一般不满足零均值条件，需要将上面的过程迭代多次，具体过程如下：

（1）对一个信号x(t)，先找出他的所有极大值点和极小值点，然后利用三次样条插值函数求得信号极大值和极小值的包络，并求两个包络的均值，得到一条均值曲线m1(t)，令h1(t)=x(t)-m1(t)，得到h1(t)；

（2）对得到的h1(t)进行步骤（1）操作，得到h2(t)，并重复10 次，得到第一个ⅠMF 分量e1(t)，并令r1=x(t)-e1(t)，r1为对应于第一个ⅠMF分量的剩余分量；

（3）对r1(t)进行步骤（1）、（2），得到e2(t)，e3(t),…,em(t)和r(t)，r(t)为剩余分量；

结果，EMD将信号x(t)分解成如下形式：

其中，k为ⅠMF 分量的阶数，而m为总的ⅠMF 分量个数，r(t)为剩余分量［19-20］。

2 基于EMD的消除ECG信号BW的算法

2.1 ⅠMF分量频率指标的确定

由于ECG 信号的BW 噪声主要由患者的呼吸或运动引起，相对ECG 信号来说是一种低频噪声，所以由以上理论可知，通过EMD 分解，低频的BW 噪声被分解到高阶ⅠMF分量中。

前面提到每一个ⅠMF 分量都有一个连续的频率范围，而且高阶的ⅠMF分量频率范围低于低阶的ⅠMF分量（可能会有重叠）。如果将每一个ⅠMF 分量相邻两个极大值点或极小值点中间的时间尺度视为其局部的频率指标，定义为Sk，k为第k个时间尺度，求出这些时间尺度的平均值，定义为每一个ⅠMF 分量的“周期”，并将其作为ⅠMF 分量的频率指标，用Tn来表示，n为ⅠMF阶数。有：

其中N为第n个ⅠMF分量中Sk的个数。

这里必须指出，理论上高阶ⅠMF 分量的Tn小于低阶ⅠMF 分量，但有必要进行验证，分别使用GSTA（一组表征年平均全球表面温度异常的数据）［8］和MⅠT-BⅠH 心率失常数据库的100 信号（MLⅠⅠ导联，取3 600 个数据点）进行验证。由表1 可以看出无论GSTA数据还是MⅠT-BⅠH 100数据的ⅠMF的“周期”都是随着ⅠMF 分量阶数的升高而增大的，说明这里定义的“周期”在实际使用中与理论一致，可以反映ⅠMF分量的频率指标。

2.2 基于EMD的消除ECG信号BW的算法

基于第2.1 节给出的ⅠMF 分量的“周期”Tn，得到一种自适应的ECG 信号的BW 去除算法。假设一种情况，1 个ECG 信号中只有1 个心率周期的数据，通过EMD 分解之后最高阶的ⅠMF 分量的“周期”理应是1个心率周期的长度，拓展到1个ECG 信号中有多个心率周期的数据，当通过EMD 分解之后产生的ⅠMF 分量中，有理由认为包含ECG 信号的ⅠMF 分量的“周期”不大于心率周期，即EMD 分解产生的“周期”大于ECG 信号周期的ⅠMF 分量都认为是BW 噪声的分解。

这里值得一提的是，通过对理想ECG 信号模型进行EMD 分解，发现ECG 本身存在直流分量，即理想ECG 信号模型最后的剩余分量（最后一个ⅠMF 分量）（图1）。我们有理由认为，医学上测量的ECG 信号存在直流分量，同时通过EMD 分解，ECG 信号的直流分量被分解到最后一个ⅠMF 分量中。所以在基于EMD 的ECG 信号BW 去除算法中，直接将最后一个ⅠMF 分量去除是不合适的。图1 中ⅠMF4 和ⅠMF5中明显存在由于端点效应造成的扰动，在ⅠMF6 中边界效应造成的扰动向内移动，在之后的ⅠMF 分量都是接近零的一条直线，直到最后的ⅠMF12，明显偏离零值，在0.12 附近，可以视为模拟标准ECG 信号的“直流部分”。

鉴于上述分析，我们不能直接将最后一个ⅠMF分量舍弃，由于ECG信号是类周期信号，将一段ECG信号分成一个一个的周期，每个周期单独进行EMD分解，每个周期分解得到的最后一个ⅠMF 分量是几乎一样的，这样理论上一段ECG 信号经过EMD 分解得到的最后一个ⅠMF 分量应该近似是一条直线。而实际上，EMD 分解ECG 信号得到的最后一个ⅠMF 分量本身就近似一条直线，同时由于BW 来源于患者呼吸或在测量信号过程中的一些动作，EMD 分解后不可能存在近乎直流的分量，所以EMD 分解ECG 信号的最后一个ⅠMF分量就是来源于ECG信号本身的。

于是，本文基于EMD 的消除ECG 信号BW 的算法具体步骤如下：

（1）使用EMD 方法对信号x(t)进行处理，得到n个ⅠMF分量；

（2）计算每个ⅠMF 分量的“周期”Tn，如式（2）所示；

（3）计算x(t)的平均R-R 周期，作为x(t)的心率周期t；

（4）将大于t的Tn对应的ⅠMF 分量去除，但保留最后一个ⅠMF 分量，然后将剩余的ⅠMF 分量叠加，得到去除基线漂移噪声的心电信号xdelBW(t)。

3 仿真结果分析与应用

3.1 定性分析

为直观看到本文算法的滤波效果，选用美国麻省理工学院MⅠT-BⅠH 心率异常数据库中的100 信号和103 信号（都采用MLⅠⅠ导联，取65 536 个数据点），用本文算法进行处理效果如图2 和图3 所示。为方便观察滤波效果，这里没有将最后一个ⅠMF 分量加入滤波后的信号，因为最后的ⅠMF 分量是原始ECG信号的直流成分，去除之后只是将滤波后的信号向上平移，与原始信号交错开，方便观察。

由图2 和图3 可以看出经过本算法处理的ECG信号明显平稳了许多，那些明显的非ECG 的趋势信号被滤除，而ECG信号的细节得以保留。

3.2 定量分析

为进一步验证本文算法的性能，需要进一步的定量分析。因为MⅠT-BⅠH 心率异常数据库中的数据本身含有的BW 噪声未知，许多研究都会对数据库中的数据进行一些处理，以得到近似没有BW 噪声的信号，但这些处理显然无法完全得到干净的信号，致使后续的实验和研究不能让人完全信服。所以本文采用ECGSYN［9］生成的模拟理想ECG 信号x(t)（图4a），其中加入的BW为：bw= 0.1sin(( 2π/4096)t)），并加入低频的正弦波作为模拟的BW 噪声bw(t)，得到带有BW 噪声的ECG 信号X(t（)图4b），这样就提前已知了信号和噪声的全部信息。

将合成的信号用EMD方法和本文算法进行处理后，得到信号X′(t)和噪声信号bw′(t)。通过x(t)与X′(t)相关系数，bw(t)与bw′(t)的相关系数，基线矫正率3 个指标来判断本文去除BW 算法的性能，并与文献［1］中的相应算法进行比较。其中，相关系数越高，说明信号还原越精确；基线矫正率越低，滤出的BW越完整。

图1 EMD分解模拟标准ECG信号Fig.1 Decomposition of analog standard electrocardiogram(ECG)signals by empirical mode decomposition(EMD)

图2 MIT-BIH 100原始数据及去除基线漂移后的数据Fig.2 MIT-BIH 100 raw data and the data after baseline wander removal

相关系数：

基线矫正率：

其中，分别用ρnn′表示n(t)与n′(t)的相关系数，用ρxX′表示x(t)与X′(t)相关系数，BR表示BW率。

由表2 可知，本文算法分离出的ECG 信号X′(t)和BW 信号n′(t)与原始的理想ECG 信号x(t)和BW信号n(t)几乎完全一致，相较于文献［1］的算法有更佳的滤除BW的效果。

4 结束语

本文在前人研究基础上，针对基于EMD 方法滤除ECG 信号BW 算法存在的不足，提出一种更加贴合EMD方法物理意义的自适应的滤除ECG信号BW的算法。在提出ⅠMF 分量平均周期的概念，以及阐述它与ECG 信号关系的基础上，通过排除不包含ECG 信号信息的ⅠMF 分量来滤除ECG 信号中的BW噪声。本文首先使用ECG 数据验证了提出的ⅠMF 分量平均周期概念的实用性，之后通过本算法处理大量ECG 数据直观体现滤除BW 噪声的效果，最后，用理想ECG 信号叠加正弦信号合成的信号进行定量分析。经实验验证，本文提出的方法由于贴合EMD 的物理意义和其自适应性，相较于其他基于EMD 算法的ECG 信号BW 去除算法有很大的优势，不论是信号的相关系数还是BW 噪声的相关系数都有所提高。不过，EMD的端点效应的影响依然较大，是算法得到信号的误差的主要来源。

图3 MIT-BIH 103原始数据及去除基线漂移后的数据Fig.3 MIT-BIH 103 raw data and the data after baseline drift removal

图4 理想心电及叠加的模拟基线漂移信号Fig.4 Ideal ECG and superimposed analog baseline drift signal

表2 文献［1］算法与本文算法性能指标的比较Tab.2 Comparison of performance indexes between the algorithm in literature［1］and the proposed algorithm