优化VMD与NLM结合的信号去噪

尹佳璠，陈小奇，李世林，邵 磊，李 季，刘宏利+

(1.天津理工大学 电气电子工程学院 复杂系统控制理论及应用重点实验室，天津 300384；2.天津喜诺生物医药有限公司 仪器研发部，天津 300457)

0 引 言

人体心电信号是判断人类心脏疾病的重要依据。且心电信号具有不平稳、随机性、非线性等特点。在采集过程中又易被噪声干扰，例如周围环境和机器内部电场磁场导致的工频干扰，因此，需要抑制ECG信号中噪声干扰。

针对ECG信号微弱、不平稳的特性，经验模态分解(EMD)无疑是最合适且应用最广泛的，但是，经验模态分解在进行信号重构时会存在模态混叠的缺陷。后续也有研究人员通过添加白噪声的方法优化经验模态分解。文献[1]利用改进的小波阈值去噪方法与EEMD结合，并利用过零率选取合适的固有模态分量。文献[2]提出了一种变分模态分解的参数选择的方法，文献[3]运用CEEMDAN与小波阈值结合的去噪方法对ECG信号进行去噪处理。但这些方法的不足之处在于噪声分量和信号分量分离的不完全以及去噪效果不够理想等缺点。研究人员提出了一种全新的变分模态分解(variational mode decomposition，VMD)，解决了模态混叠的问题。但是，这种分解方法需要通过确定合适的参数才能达到较好的分解效果。传统的参数选取方法主要通过经验选取，参数选取非常不准确。

因此，本文设计了一种参数优化的变分模态分解算法，将原始信号分解，引入样本熵的概念计算分解后信号分量中含噪声的比例，将噪声分量与信号分量区分开，对噪声分量进行非局部均值去噪(non-local means，NLM)，得到的去噪后的信号分量与低阶分量重构，得到去噪后的ECG信号。通过对实际信号进行多种方法分解重构并与原信号对比，验证了参数优化的变分模态分解后的效果最好。

1 变分模态分解原理

变分模态分解(variational mode decomposition，VMD)通过迭代法搜寻变分模型中的最优解，确定模态分量及其中心频率和带宽，变分模态分解本质上是基于维纳滤波、Hilbert变换以及频率混合推广得来的[4]。运用这种方法可以将信号分解为K个模态分量，每个分量的中心频率为ωK，将固有模态分量的带宽问题变成具有约束的变分问题。采用这种不同于经验模态分解的处理方式，不直接从已有的信号中分离，因此有效避免了类EMD算法中不可避免的模态混叠和端点效应。变分模态分解算法步骤基本可以总结为构造变分模型和模型的求解两步，理论推导如下。

1.1 构造变分模型

对原始信号x(t)分解得到K个模态分量uK(t)，进行Hilbert变换，得到解析信号

(1)

对每个解析后的信号与预估中心频率e-jωKt混合，然后将频谱调制到频带上

(2)

计算以上信号梯度的平方L2范数得到约束条件

(3)

(4)

其中：{uK}={u1,…,uK}是分解后得到的K个IMF分量；{ωK}={ω1,…,ωK}为每个分量的中心频率。δ(t)为单位脉冲函数；j为虚数单位；*表示卷积运算；∂K为计算偏导数。

1.2 变分模型求解

对于式(3)、式(4)得到的约束条件，引入惩罚因子α和Lagrange算子λ(t)，然后利用交替方向的乘子算法求出增广函数的鞍点，然后通过迭代的方式求出各个参数的最优解，最后逆变换得到各个模态分量。惩罚因子可以保证信号收到噪声干扰时进行分解重构后的信号具有较高的精度，Lagrange算子λ(t)可以保证约束的严格性[5]，将两个参数有带入后，得到

(5)

运用交替方向乘子算法(alternate direction method of multiplers, ADMM)，交替更新{uK}，{ωK}，λ

(6)

运用傅里叶变换，假设在已知其它变量的情况下，一次求解一个变量。因此，在n+1次循环后可得

(7)

(8)

(9)

2 变分模态参数优化方法

2.1 K-L散度

使 K-L散度又称相对熵，可以判断两个概率分布的相似程度，等价于两个概率分布的信息熵差值，公式如下

(10)

式中：p(xi)为真实数据的概率分布，q(xi)为理论概率分布。两组数据越相近则相对熵越小，两组数据分布差别越大，则相对熵越大[6]。

2.2 变分模态分解的参数优化方法

根据变分模态分解方法和理论推导，变分模态分解算法虽然克服了传统经验模态分解及其改进方法的缺点，但分解前需要设定几个参数值才能得到最好的分解效果，其中分解层数K和惩罚因子α对分解结果的影响最大。通常情况下，进行变分模态分解前需要人为选取K，α的值，但是人为选取无法使变分模态分解达到最好的效果。因此，简单直接选取参数是急需解决的问题。文献[5]提出了一种运用遗传算法优化的变分模态分解参数。文献[7]运用粒子群算法优化变分模态分解中参数的选取。文献[4]提出了改进变分模型并通过粒子群算法选取参数的优化变分模态分解算法。但无论是基于遗传算法还是粒子群算法的优化，除了这两种算法都容易陷入局部最优的缺陷以外，其运算所需要的时间也极长，这和算法本身的十分复杂有关。因此，这几种方法无论是从计算复杂度运算时间，还是可操作性上都不够理想。

为了寻找一种思路简单又便于运行的算法，本文根据文献[8]提出的思想，结合稀疏字典学习的思路。改进了停止条件，通过固定分解层数，优化惩罚因子，取得最优参数。

2.3 参数优化流程

参数优化流程如图1所示，具体操作方法主要分为两步：首先确定惩罚因子α，设定模态K的取值范围，K的取值范围设定为[2,20]，初始值为2，对原始信号进行变分模态分解，再通过计算确定每一次分解所得的本征模态分量(IMF)的相对熵，获得相对熵最小值时所对应的K的最优解。再根据得到最优的K值，确定惩罚因子α的最优解。设定α的范围为[100,2000]，步长为50。再寻找最小相对熵所对应的α的值，是为最优的α值。此时，得到的K，α值都是最优的。变分模态分解的参数便确定了下来。这个方法简单且易于执行，用时少，运行逻辑简单。在代码编写和数据处理上更加方便。

图1 参数优化流程

3 样本熵

由于经过分解后的信号分量包含了信号主导分量和噪声主导分量，本文通过对信号分量求样本熵，判断信号分量的中噪声含量的大小，样本熵的计算原理如下：

(1)对一个信号序列{x(i),i=1,2,…,N},假定一个m维向量序列

Xm(i)={x(i),x(i+1),…，x(i+m-1)},1≤i≤N-m+1

(11)

(2)将向量Xm(i),Xm(j)中对应向量的绝对差值设为d[Xm(i),Xm(j)]，即

d[Xm(i),Xm(j)]=
maxk=0,……,m-1(|x(i+k)-x(j+k)|)

(12)

(3)给定一个参数r，对于每个Xm(i)，确定d[Xm(i),Xm(j)]

(13)

(4)定义B(m)(r)为

(14)

(5)令m=m+1，重复以上述4个步骤，得到B(m+1)(r);

(6)样本熵的定义为，在同一参数r上，有

(15)

当N为有限值时

(16)

4 非局部均值去噪

非局部均值去噪算法(non-local mean)能够有效去除噪声并保留信号的真实性。该算法是利用了图像的自相似性，即设置一个搜索窗口，在这个搜索窗口内寻找尽可能多的相似值，然后运用加权平均值得到真实的值并同时去除了噪声干扰[9]。通过样本熵判断出噪声主导的信号分量并对这些分量进行非局部均值(non-local mean)去噪处理得到去噪分量，在去除噪声干扰的同时，保证了信号的细节。运算原理如图2所示。

图2 一维NLM原理

(17)

(18)

式中：B是以i点为中心的相似区域，L为领域B中包含点的个数，λ是滤波器参数。

5 仿真结果

为了验证本文方法的去噪效果，本文首先对仿真信号进行去噪实验。仿真信号采用MATLAB中自带的Bumps信号进行实验，采样点数为1000。Bumps信号与ECG信号相似，都是非平稳非线性的连续信号，具有代表性，可以验证本文算法的优越性。经过前文中提出的方法，确定了仿真信号的最优分解参数为：K=7，α=250。向仿真信号中添加信噪比为20 dB的高斯白噪声，分别对比传统的小波阈值去噪、文献[1]提出的EEMD-小波阈值去噪、文献[3]提出的CEEMDAN小波阈值去噪以及文献[9]的传统NLM去噪效果。由于篇幅限制，本文只截取了最终的去噪结果作为对比，去噪效果如图3所示。

图3 Bumps信号去噪

图3提供了Bumps信号在加入白噪声信号后的去噪图像，对比图3中各个方法去噪后的结果能够发现，本文提出的方法在信号保真度和去噪效果上都较好，有效抑制了高频部分的噪声和毛刺，并且完全还原了峰值部分的原始信号状态，基线部分的噪声也去除了大部分，但还残留有少量的噪声。相比之下，图3(b)中小波阈值去噪效果最差，曲线峰值上的噪声虽然处理干净但是也损失了大量的有效信号，基线部分的噪声去噪去除效果更差，甚至改变了原始信号。图3(c)中EEMD与小波阈值结合去噪方法同样在基线部分的去噪表现不好，有大量的毛刺起伏，基线部分不够平稳和圆滑。图3(d)中CEEMDAN与小波阈值结合的去噪方法相比较EEMD-小波阈值去噪效果好了很多，但仔细观察会发现毛刺现象比较严重，而波峰顶部去噪后也有明显的信号失真，真实性不如VMD-NLM去噪。图3(e)中传统的NLM去噪方法虽然抑制基线部分毛刺比较出色，但是明显严重损害了信号的真实性，幅值出现了较大的改变。因此，其它文献中改进的几种方法去噪效果都不如本文算法。为了量化去噪效果，使得去噪效果更加的直观。本文通过计算信号的信噪比(SNR)和均方误差(MSE)两个指标判断去噪情况。

为了验证去噪效果的准确性以及合理性，降低偶然性，本文进行了重复性实验，通过向仿真信号中分别加入信噪比为5 dB、10 dB、15 dB、20 dB的噪声信号，将各种方法去噪后的信噪比和均方误差列出成表。根据表1可以直观得出，本文算法在Bumps仿真信号的去噪过程中，均能够得到最高的信噪比数值和最低的均方误差数值。在不同的噪声干扰下，去噪情况都优于其它算法。本文实验环境为matlab2017a，win10-64位系统，8 G运行内存。详细结果见表1。

表1 去噪效果比较

6 实际心电信号去噪

为了验证本文去噪算法在实际采集的信号中的效果，验证算法的可行性。本文从麻省理工大学MIT-BIH心律失常数据库中提取了噪声较多的112 m实际ECG信号为实验信号，该数据库共享了医院采集的众多心电信号，在国际上具有极高的权威性。为了便于观察，本文截取了1000个采样点，采样频率为360 Hz进行实验验证。

为了验证本文提出的优化参数的变分模态分解算法与其它分解算法相比所具有的优越性，对ECG信号分别使用参数优化的变分模态分解算法、经验模态分解算法(EMD)、EEMD和CEEMDAN算法分解重构后比较，对比重构后信号与原始信号的相似程度。根据前文提出的方法，寻得最优参数为K=17,α=100，由于篇幅限制，只给出了各个方法重构后得信号与原始信号的相关系数对比表见表2。

表2 重构效果

据表2提供的重构信号和原始信号的相关系数，可直观得知，经验模态分解(EMD)算法重构后的效果相对而言最差，而经过改进的EEMD和CEEMDAN算法效果逐渐变好，但是经过参数优化后的变分模态分解算法效果最好，本文算法的重构信号与原始信号的相关系数能达到100%，各处的细节都有较高的还原度。经过以上步骤能够验证经过本文优化过的分解算法能够完美还原原始信号。而其它的几种常用的分解方法虽然效果也不错，但是仍然比本文优化的算法稍逊一筹。

然后对含有噪声的ECG信号进行优化参数的变分模态分解(VMD)去噪，如图4所示，再对每个分量求取样本熵如图5所示，本文选取的样本熵参数为：嵌入维数m=1，相似度阈值r=0.15*A，A为每个模态分量的标准差。根据样本熵与信号噪声的关系，样本熵数值越大，则信号中噪声含量越多。样本熵数值越小，有效信号含量越多的原理。有效信号主要集中在第一、二，十一，十二，十三、十四分量，余下的分量均为噪声主导的分量。

图4 VMD分量

图5 样本熵曲线

对噪声主导分量进行非局部均值去噪(non-local means)，根据文献中提供的NLM相关参数的确定方法，滤波参数h=0.5*σ，σ为噪声的标准方差。经过本文的多次实验，确定了合适的非局部均值去噪参数。NLM参数设定为：搜索窗口t=5，邻域窗口f=5，滤波参数h=200。将去噪后的分量与有效信号主导分量重构得到去噪信号如图6所示。

图6 去噪结果

为了直观展示本文算法与第5节中文献中提出的算法在去噪效果上的差异。根据传统的判断信号噪声含量的指标主要时信噪比(SNR)和均方误差(MSE)

(19)

(20)

其中，si为原信号，fi为去噪信号，N为信号长度。

通过图7与表3所给出的图形和数据，可以看出本文提出的去噪方法无论在图像表现上还是量化的数值比较上都比其它算法更为优秀，信号保真度较高，在去除噪声的同时保证了有效信号不流失。与本文算法相比，图7(a)中小波阈值去噪在实际心电信号上的去噪表现一般，噪声去除不够完整，且改变了信号的完整性，曲线也显得十分扭曲。图7(b)中EEMD与小波阈值结合的去噪方法在低频信号处去噪效果不够完美，在S波后的第一个波峰上的处理不好，毛刺较多，这一部分基本没有去噪，但在波谷的噪声信号去除的不错曲线比较圆滑。图7(c)中CEEMDAN与小波阈值结合的去噪方法改变了信号的幅值，使信号出现了一些震荡，同时在Q波波谷的去噪效果一般，保留了些许毛刺噪声。图7(d)中非局部均值去噪效果不错，但是信号失真严重，有效信号流失较多，并且出现了漂移现象。因此，通过对比几种不同的算法在编号102 m心电信号上的去噪效果，根据表3所示的信噪比与均方误差的数值，可以明显发现，本文提出的改进参数的变分模态分解算法，无论是在信噪比还是均方误差上都要优于其余算法。在去除了毛刺和噪声的同时，保证了曲线的平滑性而心电信号中最关键的Q、R、S波收到的损害最小，保证了心电信号的真实性和可靠性。

表3 去噪结果

图7 其它算法去噪结果

为了进一步验证本文算法的普适性并且减少偶然性，提高准确性和可信度，本文从MIT-BIH心率失常数据库中分别选取了100 m、102 m、103 m、109 m、116 m、113 m等6种不同的心电信号进行重复性去噪实验，分别列出几种算法的去噪结果。由于篇幅限制去噪后的信噪比和均方误差见表4(由于表格长度限制，信噪比精确到小数点后3位，均方误差精确到小数点后4位)。

表4 不同信号去噪效果对比

7 结束语

本文提出了一种参数优化的变分模态分解，通过确定一个变量搜索另一个变量的方式，利用相对熵作为停止条件，寻找最优的参数。根据优化过的变分模态分解将实际ECG信号分解得到多个模态分量，再通过计算各个模态分量的样本熵，确定噪声主导分量和信号主导分量，对噪声主导的分量进行非局部均值去噪，最后将信号分量和去噪分量重构，得到去噪后的信号。对比分析EEMD与小波阈值结合的去噪方法、CEEMDAN与小波阈值结合的去噪方法以及普通非局部均值滤波等几种先进算法，本文提出的变分模态分解和非局部均值结合的去噪方法能够自适应不同的信号，能得到最好的信噪比和均方误差，在最大程度上保存信号的真实性。