eemd优缺点_改进EEMD算法在心电信号去噪中的应用

引言

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

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

EEMD方法的具体过程[为

(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[建议的分解次数为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 EEMD算法框图 Fig. 1 Block diagram of EEMD algorithm

1.2 马氏距离

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

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

$

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算法通过反复迭代求得最佳解

(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 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 原始ECG信号和加噪ECG信号 Fig. 3 Original ECG signal and noised signal

对含噪ECG心跳进行EEMD分解，IMF分量图如

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

各IMF分量和原始ECG数据PDF间的马氏距离如γ=3界值，这表示第一、第二和第三IMF分量是噪声分量，其余IMF是信号分量。

图 5 马氏距离 Fig. 5 Mahalanobis distance

[去噪效果优于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为采样点数。

图 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

4 结束语

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