数据采集与处理  2018, Vol. 33 Issue (4): 646-653   PDF    
改进EEMD算法在心电信号去噪中的应用
潘广贞 , 王凤 , 孙艳青     
中北大学软件学院, 太原, 030051
摘要: 集合经验模态分解(Ensemble empirical mode decomposition,EEMD)方法在去除心电信号噪声时,噪声本征模态函数(Intrinsic mode function,IMF)分量难以选择且将噪声分量直接去掉会导致信号失真。针对上述问题,提出了一种基于EEMD的自适应阈值算法。首先对含噪心电图(Electrocardiogram,ECG)数据进行EEMD分解,得到IMF,根据马氏距离进行信号IMF分量和噪声IMF分量的判定,然后通过果蝇优化算法确定噪声IMF的阈值,将经过阈值去噪的新的分量和剩余分量重构得到去噪后的ECG。最后,使用MIT-BIH数据库中的心电数据进行实验,实验结果表明,该方法在去噪同时能够较好地保留信号细节。
关键词: 集合经验模态分解    马氏距离    果蝇算法    心电信号    去噪    
Application of Improved EEMD Algorithm in ECG Signal Denoising
Pan Guangzhen, Wang Feng, Sun Yanqing     
School of Software Engineering, North University of China, Taiyuan, 030051, China
Abstract: In order to solve the problems that the intrinsic mode function(IMF) components are difficult to select and the noise components are always eliminated directly when removing the noise of the electrocardiogram(ECG) signal by using the ensemble empirical mode decomposition(EEMD) method, an adaptive thresholding algorithm based on EEMD is proposed. Firstly, the noisy ECG signal is decomposed to obtain the IMFs by the EEMD method, and then the noise IMFs and the siginal IMFs are judged according to the Mahalanobis distance. After that, the thresholding of the niose IMF is determined using the fruit fly optimization algorithm(FOA). The denoised ECG signals are reconstructed by the new IMFs and the rest of IMFs after thresholding denoising. Finally, the method is applied to ECG data in MIT-BIH database. The experimental results indicate that the method can preserve the signal details while denoising.
Key words: EEMD    Mahalanobis distance    FOA    ECG signal    denoising    
引言

心电图(Electrocardiogram,ECG)是医学领域一种重要的疾病诊断工具,是判断个人健康的重要依据。ECG数据采集于人体体表,在采集过程不可避免会受到噪声影响,基线漂移、工频和肌电是主要的噪声来源。这给心脏疾病诊断和分析带来巨大的困扰:基线漂移导致ST段偏离基线会被误诊为心肌梗死、冠状动脉供血不足等疾病;工频干扰对P、T波段的影响易被判别为高、低血钾症或心肌缺血、冠心病、高血压等疾病;肌电干扰会掩盖ECG心跳中的细节,从而弱化某些心脏疾病特征。噪声干扰不仅影响医生对心脏疾病的判断,也会对计算机辅助诊断过程中的特征提取及疾病自动识别造成困扰。因此消除掺杂在ECG心跳中的噪声干扰显得尤为重要[1]

ECG常用消噪方式有形态滤波、维纳滤波、卡尔曼滤波、经验模态(Emprical mode decomposition,EMD)以及小波阈值滤波等[2-5]。小波变换具有时频多分辨功能和较好的数据性,是处理诸如ECG等生物学数据的良好工具,但它对信号分解不具有自适应性[6]。EMD算法根据数据自身特点进行分解,但其带来的模态混叠效应无法避免[7]。集合经验模态分解(Ensemble empirical mode decomposition, EEMD)在EMD方法基础上稍作改进,有效提升分解效率[8]。EMD算法和EEMD算法在ECG去噪领域取得较好效果[9],艾延廷等人[10]将马氏距离与EEMD相结合有效区分出噪声IMF分量,但直接舍去噪声IMF分量会损失一部分信号;Nguyen等人[11]结合遗传算法较好的全局搜索性能对噪声IMF分量进行自适应阈值去噪,但遗传算法容易陷入局部最优。

针对传统EEMD算法在去除ECG噪声时存在的信、噪分量难于区分和去噪阈值难以确定的问题,本文对EEMD算法进行改进。含噪ECG数据经EEMD分解后得到一系列本征模态函数(Intrinsic mode functions, IMFs),针对噪声层和信号层难以区分问题,引入马氏距离;为更好地处理噪声IMF分量,通过果蝇算法计算每个噪声IMF最佳阈值,对其中每个分量进行去噪。最后,采用来自MIT-BIH的ECG数据进行对比实验。

1 基本原理 1.1 EEMD算法

EEMD算法在处理生物信号方面取得良好效果,不仅继承EMD算法自适应分解信号的优点,且有效避免EMD算法的模态混叠。该方法通过对原始ECG执行多次EMD分解,并在每次分解时加入白噪声,分解效果与分解次数呈正比。这些IMF分量可以用于频谱分析,IMF的频率随着其指数的增加而降低[8, 10]

EEMD方法的具体过程[12]

(1) 向原始数据x(t)中加入白噪声w(t),得到

$ y\left( t \right) = x\left( t \right) + w\left( t \right) $

(2) 对信号y(t)进行EMD分解得到

$ y\left( t \right) = \sum\limits_{i = 1}^N {{\rm{im}}{{\rm{f}}_j}\left( t \right) + {r_N}\left( t \right)} $

(3) 重复上述步骤(Huang[12]建议的分解次数为100次),得到

$ {y_i}\left( t \right) = \sum\limits_{i = 1}^N {{\rm{im}}{{\rm{f}}_{ij}}\left( t \right) + {r_{iN}}\left( t \right)} $

(4) 对上述的结果求平均,得到最终的IMF分量

$ {\rm{im}}{{\rm{f}}_j}\left( t \right) = \frac{1}{n}\sum\limits_{i = 1}^n {{\rm{im}}{{\rm{f}}_{ij}}\left( t \right)} $

最终EEMD的分解结果为

$ x\left( t \right) = \sum\limits_j {{\rm{im}}{{\rm{f}}_j}\left( t \right) + {r_N}\left( t \right)} $

利用EEMD算法去除ECG噪声步骤如图 1所示。

图 1 EEMD算法框图 Fig. 1 Block diagram of EEMD algorithm

1.2 马氏距离

马氏距离(Mahalanbis distance, MD)能够计算两个位置样本集的相似度,它对异常数值的敏感性使得它适合作为相似度测量工具,可以使用该度量工具判断两个一维信号概率密度函数(Probability density function, PDF)之间的相似性[13]

XY是从均值为μ,协方差矩阵为Σ的总体G中抽取的两个样品,则XY两点之间的马氏距离为

$ d_m^2\left( {X,Y} \right) = \left( {X - Y} \right){\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}^{ - 1}}\left( {X - Y} \right) $
1.3 果蝇优化算法

果蝇优化算法(Fly optimization algorithm, FOA)是一种新的寻优方法,该算法基于果蝇的觅食行为找到全局最优解,克服了其他寻优方式容易求得局部最优解缺陷,有更好的全局寻优性。由于该方法具有适应性强、简单便于实施等特点,使其得到广泛的应用。果蝇的视觉和嗅觉优于其他物种,在觅食过程中向食物气味浓度最高的方向飞去,在飞行过程中飞向觅食能力最强的果蝇个体。根据果蝇群体觅食过程,FOA算法通过反复迭代求得最佳解[14-15],其寻优步骤如下:

(1) 初始化参数:初始化整个种群的位置范围(LR)、活动范围(FR)、群体大小、迭代次数上限。可以通过以下等式获得初始群体位置(x0, y0)。

$ {x_0} = {\rm{rand}}\left( {{\rm{LR}}} \right) $
$ {y_0} = {\rm{rand}}\left( {{\rm{LR}}} \right) $

其中rand为随机数生成函数。

(2) 给予每个个体在觅食过程中确定飞行方向和计算距离的能力,其位置计算为

$ {x_i} = {x_0} + {\rm{rand}}\left( {{\rm{FR}}} \right) $
$ {y_i} = {y_0} + {\rm{rand}}\left( {{\rm{FR}}} \right) $

(3) 计算果蝇位置味道浓度:求个体在当前点气味浓度判断数(Si)和距离(Disti)。

$ {S_i} = 1/{\rm{Dis}}{{\rm{t}}_i} $
$ {\rm{Dis}}{{\rm{t}}_i} = \sqrt {x_i^2 + y_i^2} $

(4) 通过Si和函数Functioni求解每个个体在当前点的smelli。然后确定具有最佳气味浓度的果蝇。

$ {\rm{smel}}{{\rm{l}}_i} = {\rm{Function}}\left( {{S_i}} \right) $
$ \left[ {{\rm{bestmell}}\;{\rm{bestindex}}} \right] = \max \left( {{\rm{smel}}{{\rm{l}}_i}} \right) $

其中bestindex表示具有bestsmell的个体序号。

(5) 记录此时的bestsmell和序号bestindex的个体位置坐标,让剩余个体向该最佳点飞去。

(6) 重复步骤(1—5),判断每次得到的bestsmell,超出循环次数后,停止迭代,从而得到最优解bestsmell。

(7) 进入迭代优化,重复上述所有步骤,对每次迭代所得到的味道浓度进行比较分析。

2 本文算法

传统EEMD算法在处理ECG时,人为区分信、噪IMF分量会造成一定误差,对区分出的噪声IMF直接舍弃,而噪声IMF中往往还具有一定信息量,直接丢弃会导致信号失真,所以对噪声IMF的正确判断和合理处理是本算法提升去噪效果的关键。

2.1 噪声IMF的确定

EEMD算法将ECG分解为多个IMF,其中包括少数噪声IMF以及信号IMF。为区分出噪声IMF,采用基于PDF和MD的方法来判断所有IMFs中噪声IMF和信号IMF分界点。马氏距离的计算规则为

$ d\left( i \right) = {\rm{MD}}\left( {{\rm{PDF}}\left( {x\left( t \right)} \right),{\rm{PDF}}\left( {{\rm{im}}{{\rm{f}}_i}\left( t \right)} \right)} \right) $ (1)

将有用IMF和含噪IMF之间边界值设为γ,将该值定义为最后一个噪声IMF所对应的数值,根据PDF间的马氏距离得到。因此马氏距离拐点(即距离骤减点)前的IMF分量即为选定的噪声分量,边界值γ定义如下

$ \gamma = \mathop {\arg \max }\limits_{1 \le i \le N} \left\{ {d\left( i \right)} \right\} $ (2)
2.2 噪声IMF阈值的选取

识别出边界值后,就可以区分出含噪IMF和信号IMF。iγ的IMF为含噪数据,其余i>γ的IMF为不含噪数据。信号IMF予以保留,每个噪声IMF经过自适应阈值去噪处理。为了判断每一个噪声IMF分量阈值,使用FOA对阈值进行全局寻优。

重构的去噪信号如下

$ \hat x\left( t \right) = \sum\limits_{i = 1}^{r - 1} {\overline {{\rm{im}}{{\rm{f}}_i}} \left( t \right)} + \sum\limits_{i = \gamma }^N {{\rm{im}}{{\rm{f}}_i}\left( t \right)} + {r_N}\left( t \right) $ (3)

式中:$ {\overline {{\rm{imf}}} _i}\left( t \right)$为噪声IMF分量。本文对噪声IMF分量采用的阈值函数为

$ \overline {{\rm{im}}{{\rm{f}}_i}} \left( t \right) = \left\{ {\begin{array}{*{20}{c}} {{\mathop{\rm sgn}} \left( {{\rm{im}}{{\rm{f}}_i}\left( t \right)} \right)\left( {\left| {{\rm{im}}{{\rm{f}}_i}\left( t \right)} \right| - {T_i}} \right)}&{\left| {{\rm{im}}{{\rm{f}}_i}\left( t \right)} \right| > {T_i}}\\ 0&{\left| {{\rm{im}}{{\rm{f}}_i}\left( t \right)} \right| \le {T_i}} \end{array}} \right.,i = 1,2, \cdots ,\gamma $ (4)

式中:Ti为依赖于IMF的通用阈值,对噪声IMF分量的阈值表达式为

$ {T_i} = C\sqrt {{E_i}2{\rm{In}}N} $ (5)

式中:C为阈值系数,是可以通过实验确定的常数;N为信号长度。

i个噪声IMF分量的能量Ei计算如下

$ \begin{array}{*{20}{c}} {{E_i} = \frac{{{E_1}}}{\beta }{\rho ^{ - i}}}&{i = 2, \cdots ,\gamma } \end{array} $ (6)

式中:βρ为与筛选迭代次数有关的参数;E1为第1个IMF分量的能量。

由式(5,6)得到每个噪声IMF的阈值为

$ {T_i} = C\sqrt {\frac{{{E_1}}}{\beta }{\rho ^{ - i}}2{\rm{In}}N} $ (7)

为了计算噪声IMF分量中的相关阈值,首先要确定式(7)中的ρβC。为了解决这个问题,本文通过果蝇优化算法寻求ρβC这3个参数的最优解,从而计算每个噪声IMF分量的最佳阈值。

2.3 基本流程

本文提出算法的基本步骤为:

(1) 对原始ECG进行EEMD分解,向其中添加均值为零、方差恒定的白噪声,则原始ECG被分解为N个IMF和1个余项r(n)。

(2) 计算各IMF分量PDF与原始数据PDF之间马氏距离di(i=1, …, N), 确定边界值γ,则认为前γ个IMF为含噪IMF;第γ+1~N个IMF为信号IMF。

(3) 使用果蝇优化阈值对含噪IMF自适应去噪,同时保留信号IMF。

(4) 对信号IMF分量和去噪后的噪声IMF分量求和,重构去噪后心电信号。

上述步骤如图 2所示。

图 2 ECG去噪算法框图 Fig. 2 Block diagram of ECG denoising algorithm

3 实验仿真与分析

本文使用的实验平台Matlab10.0搭建在Windows 7上,ECG数据源取自公开数据库MIT-BIH。实验取其中Arrthythmia Database数据库第100号记录进行仿真实验,选取采样点为1 000。实验使用改进算法对加入噪声的ECG进行去噪处理,并与EEMD算法和小波阈值法对ECG处理结果进行直观效果对比和客观参数对比。

3.1 数据处理

向MIT-BIH/100号ECG中添SNB=10 dB噪声,如图 3所示。

图 3 原始ECG信号和加噪ECG信号 Fig. 3 Original ECG signal and noised signal

对含噪ECG心跳进行EEMD分解,IMF分量图如图 4所示。图中未见模态混叠,各IMF分量中QRS波群集中,在几个高频IMF分量中清晰观察到噪声存在,可见分解效果比较理想。

图 4 EEMD分解的IMF分量图 Fig. 4 Intrinsic mode component diagram of EEMD decomposition

各IMF分量和原始ECG数据PDF间的马氏距离如图 5所示,曲线在第4个IMF处发生骤减,通过该图可识别出边γ=3界值,这表示第一、第二和第三IMF分量是噪声分量,其余IMF是信号分量。

图 5 马氏距离 Fig. 5 Mahalanobis distance

图 6显示了使用小波阈值、EEMD、EEMD-GA和本文算法对仿真ECG的处理结果。对比图 6可知,这4种方法对ECG数据的去噪均有一定的贡献。其中,小波阈值处理过后的ECG数据光滑且低平,说明去噪强度过大,携带的信息量丢失明显;EEMD算法处理后的数据中明显可见噪声,且含有部分毛刺,说明去噪强度较弱;EEMD-GA算法[11]去噪效果优于EEMD算法,但还是存在噪声;本文算法能够有效地去除ECG信号中的噪声,同时信号中的细节得到了保留,与原始信号更为接近,去噪效果优于小波阈值法和EEMD方法。但本文算法由于进行了参数寻优,在耗时上要略高于其他算法。

图 6 算法去噪效果对比图 Fig. 6 Denoising effect comparison of different algorithms

3.2 评估指标

采用均方误差MSE和信噪比SNR这两个指标进行进一步分析,MSE和SNR能够量化去噪效果,直接反映不同方法对数据的处理能力,具体表达式如下

$ {\rm{SNR}} = 10{\log _{10}}\left( {\sum\limits_{t = 1}^N {{x^2}\left( t \right)} /\sum\limits_{t = 1}^N {{{\left( {x\left( t \right) - \hat x\left( t \right)} \right)}^2}} } \right) $ (8)
$ {\rm{MSE}} = \frac{1}{N}\sum\limits_{t = 1}^N {{{\left( {x\left( t \right) - \hat x\left( t \right)} \right)}^2}} $ (9)

式中:x(t)为加噪ECG数据,$ \widehat x\left( t \right)$为重构ECG数据,N为采样点数。

图 78分别显示不同输入信噪比下ECG数据去噪后的SNR值和MSE值。由于实验结果中小波阈值去噪算法的各项参数与其他3种算法差距较大,所以图中仅显示剩余3种算法的参数对比。

图 7 输入不同SNRs各方法去噪后SNR Fig. 7 Output SNR of different denoising algorithms under different input SNRs

图 8 输入不同SNRs各方法去噪后MSE Fig. 8 MSE of different denoising algorithms under different input SNRs

图 78参数对比显示,在输入信噪比相同的情况下,本文算法的输出SNR较其他两种方法有小幅提升,且均方差最小。

4 结束语

本文将EEMD与果蝇优化算法结合,提出一种自适应阈值算法,从噪声IMF选择和处理两个方面提升算法效率,解决EEMD方法人为区分信、噪IMF的弊端,以及对噪声IMF分量阈值难以确定的问题。该算法采用马氏距离区分出经EEMD分解得到的IMF中噪声IMF,然后用经过果蝇优化的阈值对噪声IMF进行去噪。实验结果显示本文算法具有较好鲁棒性,在去除心电信号噪声上取得理想效果,可以应用到其他生物信号相似处理上。但本文算法也存在计算复杂度高、计算量较大的问题,还需要进一步改进,这也是下一步的研究内容。

参考文献
[1]
殷俊鹏, 田应洪, 赖宗声, 等. 基于小波域数字滤波的心电信号BW去噪算法[J]. 计算机工程, 2013, 39(3): 267-271.
Yin Junpeng, Tian Yinghong, Lai Zongsheng, et al. BW algorithm for ECG signals based on wavelet domain digital filtering[J]. Computer Engineering, 2013, 39(3): 267-271.
[2]
王争, 何宏, 谭永红. 基于遗传优化函数曲线的小波阈值法心电信号除噪[J]. 计算机应用, 2014, 34(9): 2600-2603.
Wang Zheng, He Hong, Tan Yonghong. Wavelet threshold method for ECG signal denoising based on genetic optimization function curve[J]. Computer Applications, 2014, 34(9): 2600-2603. DOI:10.3969/j.issn.1001-3695.2014.09.009
[3]
张磊磊. 基于EMD的心电信号去噪方法研究及实现验证[D]. 重庆: 重庆邮电大学, 2016.
Zhang Leilei. ECG signal denoising method and verification based on EMD[D]. Chongqing: Chongqing University of Posts and Telecommunications, 2016.
[4]
欧阳波, 程栋, 王玲. 改进小波阈值算法在心电信号去噪中的应用[J]. 计算机工程与应用, 2015, 51(4): 213-217.
Ou Yangbo, Cheng Dong, Wang Ling. Application of improved wavelet threshold algorithm in ECG signal denoising[J]. Computer Engineering and Applications, 2015, 51(4): 213-217. DOI:10.3778/j.issn.1002-8331.1306-0078
[5]
田絮资, 杨建, 黄力宏. 心电信号去噪的数学形态学滤波器[J]. 计算机工程与应用, 2012, 48(2): 124-126.
Tian Xuzi, Yang Jian, Huang Lihong. Mathematical morphological filtering of ECG signal denoising[J]. Computer Engineering and Applications, 2012, 48(2): 124-126. DOI:10.3778/j.issn.1002-8331.2012.02.035
[6]
Sanyal A, Baral A, Lahiri A. Application of framelet transform in filtering baseline drift from ECG signals[J]. Procedia Technology, 2016, 4(4): 862-866.
[7]
杨向林, 严洪, 许志, 等. 基于Hilbert-Huang变换的ECG消噪[J]. 电子学报, 2011, 39(4): 819-824.
Yang Xianglin, Yan Hong, Xu Zhi, et al. ECG denoising based on hilbert-huang transform[J]. Acta Electronic Journal, 2011, 39(4): 819-824.
[8]
Mariyappa N, Sengottuvel S, Patel R, et al. Denoising of multichannel MCG data by the combination of EEMD and ICA and its effect on the pseudo current density maps[J]. Biomedical Signal Processing & Control, 2015, 18: 204-213.
[9]
吴其前, 张雄伟. 基于EEMD域统计模型的话音激活检测算法[J]. 数据采集与处理, 2012, 27(1): 51-56.
Wu Qiqian, Zhang Xiongwei. Voice activity detection algorithm based on EEMD domain statistical model[J]. Data Acquisition and Processing, 2012, 27(1): 51-56. DOI:10.3969/j.issn.1004-9037.2012.01.008
[10]
艾延廷, 冯研研, 周海仑. 小波变换和EEMD-马氏距离的轴承故障诊断[J]. 噪声与振动控制, 2015, 35(1): 235-239.
Ai Yanting, Feng Yanyan, Zhou Hailun. Wavelet transform and EEMD-Markov distance bearing fault diagnosis[J]. Noise and Vibration Control, 2015, 35(1): 235-239. DOI:10.3969/j.issn.1006-1355.2015.01.050
[11]
Nguyen P, Kim J M. Adaptive ECG denoising using genetic algorithm-based thresholding and ensemble empirical mode decomposition[J]. Information Sciences, 2016, 10: 499-511.
[12]
宋美. 基于集合经验模态分解和小波收缩算法的自适应心电信号去噪问题研究[J]. 生物数学学报, 2015(4): 629-638.
Song Mei. Research on adaptive ECG signal denoising problem based on aggregate empirical mode decomposition and wavelet shrinkage algorithm[J]. Journal of Biomathematics, 2015(4): 629-638.
[13]
张翔, 王士同. 一种基于马氏距离的可能性聚类方法[J]. 数据采集与处理, 2011, 26(1): 101-105.
Zhang Xiang, Wang Shitong. A possibility clustering method based on mahalanobis distance[J]. Data Acquisition and Processing, 2011, 26(1): 101-105. DOI:10.3969/j.issn.1004-9037.2011.01.019
[14]
Wu Lijuan, Cao Guohua. Seasonal SVR with FOA algorithm for single-step and multi-step ahead forecasting in monthly inbound tourist flow[J]. Knowledge-Based Systems, 2016, 15(110): 157-166.
[15]
李大中, 张坤, 赵杰, 等. 果蝇优化小波阈值超声检测信号去噪[J]. 中国测试, 2016, 42(7): 88-92.
Li Dazhong, Zhang Kun, Zhao Jie, et al. Drosophila optimized wavelet threshold ultrasonic detection signal denoising[J]. China Test, 2016, 42(7): 88-92.