CN110680308B - 基于改进emd与阈值法融合的心电信号去噪方法 - Google Patents

Abstract

本发明提出了一种基于改进EMD与阈值法融合的心电信号去噪方法，属于信号滤波技术领域。本发明通过叠加不同权重系数的白噪声解决模态混叠问题，通过最小二乘支持向量机的方法解决端点问题，通过保形样条插值的方法构造信号上下包络线，利用保形分段法来构造具有二阶逼近精度、分段少、运算量小的三次样条插值，该方法可以抑制包络拟合过冲/欠冲的问题，通过分解出的IMF的正交性与能量性质，提出IMF分量“筛分”终止的判据，保证了EMD分解的正交性与完备性，通过互信息的原则判定筛分出的IMF信号含有噪声的多少，来决定是否对其进行滤波处理，增加了EMD算法的快速性；改进了阈值函数，该阈值函数结合了软硬阈值的优点，对含有噪声的IMF进行滤波处理。

Description

基于改进EMD与阈值法融合的心电信号去噪方法

技术领域

本发明涉及一种心电信号去噪的方法，特别是基于改进EMD与阈值法融合进行心电信号去噪方法，属于信号滤波技术领域。

背景技术

心电信号是一种典型的非平稳、微弱生物电信号，广泛应用于各种心脏疾病的诊断与治疗。心电信号中常伴有非常严重的高频、低频噪声，且噪声频带常与心电信号频带有重叠，滤波预处理较困难。

传统的生物医学信号处理主要是以傅立叶理论为基础的。傅立叶信号处理技术在信号频谱分析方面以及与其相关联的数据压缩、信号检测、滤波等信号处理领域几乎无可替代。但傅氏变换的积分区间是由正无穷到负无穷的，它无法得到信号在某一段时间内的频谱含量。而小波变换由于其优良的时频分析特性和处理非平稳随机信号的能力，成为了处理心电等生物医学信号的一种行之有效的方法。同样，EMD由于其在分析非线性和非平稳性信号时所表现出的良好的适应性也己经开始被应用到了生物医学的处理领域。如心电图信号分析、血压信号去噪、心跳信号分析等，都已得到成功应用。

经验模态分解(Empirical modedecomposition，EMD)将信号分解成有限个本征模函数(Intrinsic mode function，IMF)之和。EMD分解充分考虑了信号本身的局部尺度特征，这样得到的每个IMF分量表示了原信号的一种尺度特征，包含了原信号真实的物理信息。EMD作为一种新的自适应信号时频处理方法，在机械故障诊断、特征提取、地球物理探测、医学分析等方面都有了广泛的应用，并且EMD方法也已扩展到二维信号处理领域。在图像边缘检测少、纹理分析、图像融合、图像压缩、图像滤波等领域都得到了很好的应用，这些都说明了EMD的有效性。

现如今的EMD方法还存在一些缺陷还需作进一步的研究。如：EMD快速算法的研究、EMD模态混叠问题、端点效应、信号的包络拟合问题以及IMF分量“筛分”终止的判据研究。方法中包络曲线拟合的好坏直接影响到的分解结果，原采用三次样条函数插值的方法是以极值点为已知点完成的，因为极值点分布的不均匀性，导致了其是非均匀点的拟合，会引起极值欠冲或过冲的问题，造成分解的较大误差。模态混叠问题是EMD使用过程中经常会遇到的一个问题，其主要由间断事件和信号的相互作用两方面原因引起。

本发明针对上述存在的缺陷，提出一种基于改进EMD与阈值法融合的心电信号去噪方法，致力于提高心电信号经验模态分解的质量并提高去除噪声的效果。

发明内容

针对现有技术中存在的不足，本发明的目的在于提出一种改进EMD与阈值法融合技术，以完成对心电信号的分解并去除心电信号的噪声。本发明通过叠加不同权重系数的白噪声解决模态混叠问题，通过最小二乘支持向量机的方法解决端点问题，通过保形样条插值的方法构造信号上下包络线，利用保形分段法来构造具有二阶逼近精度、分段少、运算量小的三次样条插值，该方法可以抑制包络拟合过冲/欠冲的问题，避免传统插值方法导致的插值误差随着分解过程的持续进行而出现的误差不断累积。通过分解出的IMF的正交性与能量性质，提出IMF分量“筛分”终止的判据，保证了EMD分解的正交性与完备性。通过互信息的原则判定筛分出的IMF信号含有噪声的多少，来决定是否对其进行滤波处理，这大大增加了EMD算法的快速性。并提出了一种改进的阈值函数，该阈值函数结合了软硬阈值的优点，对含有噪声的IMF进行滤波处理。

本发明旨在解决使用EMD分解时，造成的模态混叠问题，提出一种基于改进EMD与阈值法融合的心电信号去噪方法。

一种基于改进EMD与阈值法融合的心电信号去噪方法，包括如下步骤：

步骤1：对心电信号进行改进的EMD分解；

步骤1.1：对心电信号加入白噪声的方法可以表示为：

其中，

代表对心电信号第a次加入正系数白噪声后的信号；

代表对心电信号第a次加入与正系数相对应的负系数白噪声后的信号；y(t)代表心电信号；n(t)代表均值为零，单位方差的白噪声；ω_a为权重系数。以下步骤1.2～1.5是加入正系数白噪声的处理方法，加入负系数白噪声的处理方法与此类似，仅需将上标“+”改为“-”，含义表示与之对应的加入负系数白噪声的情况。

步骤1.2：为避免信号插值拟合时出现两端发散现象，即端点效应。通过利用最小二乘支持向量机延拓信号两端；对步骤1.1处理后的信号进行采样，首先，对

进行采样，信号采样序列为

N为采样点数。训练样本集为A＝{(h₁,g₁),(h₂,g₂),···,(h_l,g_l)}，其中：h_i为输入向量，

g_i为输出向量，

1≤i≤l；l为训练样本集个数；测试样本集B＝{(h_N-S+1,g_N-S+1),(h_N-S+2,g_N-S+2),···,(h_N,g_N)}，其中：1≤S为测试样本数；

通过测试样本集，最小二乘支持向量机的输出是

其中，

为延拓后的第一个信号，再将

作为原始数据新的边界点，以同样的方法得到第2个数据序列延拓值

以此类推，根据需要延拓数据点的个数可以得到全部延拓序列

其中，M为向右延拓的信号点数，对于给定数据序列向左延拓的方法与向右延拓的方法相同，并记向左延拓的序列为：

最终延拓后的序列记为

步骤1.3：通过保形样条插值的方法构造信号的上下包络线。

步骤1.3.1：通过比较步骤1.2得到的序列Z(t_i)中某点与其左右邻近点的大小关系来判断该点是否为极值点。具体操作如下：

当判断端点..与

时：

若

则点

为极大值点。

若

则点

为极小值点。

若

则点

为极大值点。

若

则点

为极小值点。

当判断除上述端点以外的其它点时：

若

-M+1≤i≤N+M-1，则点

为极大值点。

若

-M+1≤i≤N+M-1，则点

为极小值点。

通过上述方法可得序列

的极大值点序列为：

{x_max(t₀),x_max(t₁),...,x_max(t_b)}，其中b为极大值点个数；

极小值点序列为：

{x_min(t₀),x_min(t₁),...,x_min(t_c)}，其中c为极大值点个数。

步骤1.3.2：通过极大值点构造上包络线。

在每两个极大值点间插入两个数值点，定义如下：

其中，t_ij为位于时间t_i-1与t_i之间的第j个插入点。

插入数值点后的序列为：

{x_max(t₀),x_max(t₁₁),x_max(t₁₂),x_max(t₁),x_max(t₂₁),x_max(t₂₂),x_max(t₃),...,x_max(t_b)}

将上述序列作为T-B样条曲线的控制顶点，上包络线即T-B样条曲线可表达如下：

其中，j表示拟合的第j段T-B样条曲线；U表示上包络线；h_ij(t)为T-B样条基函数，当拟合t₀到t₁段曲线时，h_ij(t)的表达式如下：

其中：

拟合t₁:t₂、t₂:t₃、···、t_b-1:t_b段曲线时，h_ij(t)的表达式与上述类似。

步骤1.3.3：通过极小值点构造下包络线，构造方法与构造上包络线的方法类似，下包络线的表达如下：

其中，D表示下包络线；

步骤1.4：求上下包络线的均值：

其中，m_a(t)为第a次加入正系数白噪声后所得到的上下包络线的均值；

计算信号

与m_a(t)的差值：

如果C(t)不满足IMF定义截止条件，则重复步骤1.3；否则提取C(t)作为

为第a次加入正系数白噪声经第i次EMD分解的IMF；剩余量

步骤1.5：将步骤1.4的剩余量r(t)作为一个新的信号，即r(t)＝y(t)，再经步骤1.1～步骤1.4进行筛选来获得下一个更低频率的IMF，直到筛选出来的IMF分量满足如下规则，则停止筛分。

IMF分量“筛分”终止的停止准则为：

其中：IMF_i(t)为第i次EMD分解的IMF，n为IMF_i(t)的长度；

当ZSD＜θ时，停止筛分，作为优选，θ＝0.03。

步骤1.6：经上述分解后的信号可表达为如下公式：

其中：d为心电信号经EMD分解后的IMF的个数，

为第a次加入正系数白噪声经第i次EMD分解的IMF，r_d(t)为进行第d次分解时的残差。

与步骤1.1～步骤1.5相似，计算对心电信号加入负系数白噪声后的EMD分解，分解后的信号可表达为如下公式：

其中，

为第a次加入负系数白噪声经第i次EMD分解的IMF。

步骤1.7：对步骤1.6得到加入正系数及负系数白噪声后的IMF进行累加求平均，可得：

其中：

步骤2：通过排列互信息判断每个IMF与原信号的线性关系，对由步骤1得到的IMF进行有效信号判定，公式如下：

QI(IMF^(k),y)＝S(IMF^(k))+S(y)-S(IMF^(k),y) (13)

其中：S(IMF^(k))和S(y)分别为IMF^(k)(t)和y(t)的排列熵，S(IMF^(k))计算公式如下：

其中p(·)表示排列(·)的联合概率密度，计算公式如下：

其中i＝1,2,···,n-(d-1)τ，n为时间序列的长度，d为给定的嵌入维数，τ为时间延迟；#表示集合内元素数，

为IMF^(k)的一个状态向量，对应的排列是(π_i)；排列(π_i)表示时间序列IMF^(k)(t)映射到d维空间后每个状态向量所对应的一个排列顺序；

假设时间序列{IMF^(k)(t)}映射到d维相空间，定义其联合状态向量为：

其状态向量用以表征第i时刻序列的轨迹，通过此定义可以得到映射到状态空间的轨迹状态矩阵：

通过比较每个状态矩阵的行向量的近邻值的大小关系得一个排列顺序，即可得轨迹矩阵的排列矩阵：

S(y)与S(IMF^(k))的计算类似。

S(IMF^(k),y)表示IMF^(k)(t)和y(t)的联合排列熵，计算公式如下：

其中p(π_i,π_j)表示排列(π_i,π_j)的联合概率密度，计算公式如下：

其中下标i'＝1,2,···,n-(d-1)τ，

一对状态空间轨迹矩阵，其对应的排列是(π_i,π_j)；排列(π_i,π_j)表示IMF^(k)(t)和y(t)时间序列映射到d维空间后每个状态向量所对应的一个排列顺序，其表示方式如下：假设时间序列{IMF^(k)(t)}和{y(t)}映射到d维相空间，定义其联合状态矩阵为：

其状态向量用以表征第i时刻序列的轨迹，通过此定义可以得到映射到状态空间的轨迹状矩阵：

通过比较每个状态矩阵的行向量的近邻值的大小关系得一个排列顺序，即可得轨迹矩阵的排列矩阵：

根据每个IMF与原信号的互信息QI值判断每个IMF受噪声干扰的程度，其噪声量越多，QI值越小。将明显QI值小的IMF选出进行去噪处理。

步骤3：对步骤2选出来需要去噪的IMF进行阈值去噪处理，阈值公式如下：

其中：可调参数

可调参数α为正数，可调参数m为正数；阈值λ₁的计算公式如下：

阈值

σ为各个IMF分量中噪声标准差，σ＝median(IMF_i)/0.6745，median(·)为求中位数；

阈值λ₂的计算公式如下：

步骤4：步骤2剩下的IMF与步骤3进行处理后的信号进行信号重构，计算公式如下：

其中，y'(t)表示去噪后的信号。

至此，由步骤1到步骤4，完成了一种基于改进EMD与阈值法融合的心电信号去噪的过程。

有益效果

一种基于改进EMD与阈值法融合的心电信号去噪方法，相比于现有技术，本方法有如下有益效果：

1.本方法采用分别多次加入不同权重系数的正负白噪声的方法有效地解决了EMD中由于两个信号的相互作用导致的模态混叠问题，提高了算法的频率分辨率(这里所说的频率分辨率是指区分两个相邻频率分量的能力)；

2.本方法使用最小二乘支持向量机的方法对信号进行延拓并加窗，避免了利用三次样条插值方法对信号进行拟合，导致分解时出现的失真现象；

3.本方法使用保形样条插值的方法构造信号上下包络线，该方法利用保形分段法来构造具有二阶逼近精度、分段少、运算量小的三次样条插值方法来抑制包络拟合过冲/欠冲问题，避免传统插值方法导致的插值误差随着分解过程的持续进行而出现误差不断累积，造成严重误差。

4.本方法直接采用阈值法对IMF去噪，该方法有效的降低了运算时间。并改进了阈值函数，该阈值函数更好地保留有用信号的细节和边缘特性，在实现信号降噪的基础上，更好地提出了有效信号，保证信号的保幅度，提高了去噪效果，降低了运算时间。

5.本方法流程简单，易于实现，可以用于心电信号分析领域内相关软件的设计工作。

附图说明

图1为本发明的步骤示意图；

图2为本发明的流程图；

图3为本发明EMD分解的详细示意图；

具体实施方式

下面根据附图和实例对本发明进行详细说明，但本发明的具体实施方式不仅于此。

本实施例中的数据来源于中国人民解放军总医院心血管内科的患者采样数据，采用心电图仪对患者进行心电信号采集，信号时长5s，采样频率360Hz，获取患者在日常自由活动(非剧烈运动)情境下的心电信号数据。

本实施例阐述了将本发明“基于改进EMD与阈值法融合的心电信号去噪方法”应用于心电信号去噪场景下的流程。

图1为本方法步骤示意图，图2为具体的流程图，从图中可以看出，先采集心电信号，然后进行如下步骤：

步骤A：对采集到的信号进行改进的EMD分解；图3为本发明EMD分解的详细示意图，具体如下：

步骤A1：对采集到的信号y(t)加入正ω₁倍的均值为零，单位方差的白噪声n(t)；

具体到本实施例，白噪声幅值为心电信号标准差的0.05倍，即ω₁＝0.5，记加入噪声后的信号为x₁(t)；

步骤A2：利用最小二乘支持向量机对步骤A.1得到的x₁(t)信号进行延拓；

对x₁(t)进行采样得采样序列{x(1),x(2),x(3),···,x(N)}，N为采样点数，训练样本集为B＝{(h₁,g₁),(h₂,g₂),···,(h_l,g_l)}。

具体到本实施例，采样点数为2000，训练样本为100，可得其训练样本：输入样本为100×1900维数据。

向右延拓的第一个预测值为x(N+1)；再将x(N+1)作为原始数据新的边界点，以同样的方法得到第2个数据序列延拓值x(N+2)。以此类推，根据需要延拓数据点的个数可以得到全部延拓序列x(N+1)，x(N+2)…，x(N+M)，对于给定数据序列向前延拓的方法与向后延拓的方法相同。

由信号的两端点处向两侧各延拓100个采样点，即M＝100得扩展序列为：{x(1),x(2),x(3),···,x(2200)}。

步骤A3：利用保形样条插值的方法构造信号包络线；

具体到本实施例，通过比较步骤A.2得到的序列中的每一点与其左右邻近点的大小关系来判断该点是否为极值点，根据比较规则，比较后的极大值序列为{x_max(t₀),x_max(t₁),...,x_max(t_b)}，其中b为极大值点个数；通过公式(2)再两个极大值点间插入两个值作为T-B样条的控制顶点，利用公式(3)构造出信号的上包络线，记为

构造下包络线的方法与上包络线的类似，记为

步骤A4：利用公式(6)计算上下包络线的均值m_a(t)；

步骤A5：利用公式(7)计算步骤A.2得到的信号x₁'(t)和均值m_a(t)的差值；

根据公式(8)判断C(t)是否满足IMF的定义截至条件，即满足当ZSD＜0.3时，则提取C(t)作为IMF；否则利用公式

计算剩余量，并使剩余量为原始信号，重复步骤A1—步骤A5，直到剩余分量为单调函数时停止筛分。

步骤A6：改变步骤A.1的权重因子ω，并进行m次步骤A.1～步骤A.5的筛分；

具体到本实施例，白噪声幅值可为心电信号标准差的0.05倍，0.06和0.04倍，加入共100次；

步骤A7：利用公式(12)可得信号经EMD分解的所有IMF分量；

步骤B：通过排列互信息对得到的每个IMF进行有效信号判定；

步骤B1：求时间序列的排列熵；

具体到本实施例，嵌入维数为1001，时间延迟为1。根据排列概率密度函数的香农熵计算出时间序列IMF^(k)(t)和y(t)的排列熵S(IMF^(k))和S(y)；利用公式(18)计算出联合时间序列(IMF^(k)(t),y(t))的联合排列熵S(IMF^(k),y)。

步骤B2：利用公式(13)计算排列互信息；

步骤B3：根据每个IMF与原信号的互信息QI值判断每个IMF受噪声干扰的程度，其噪声量越多，QI值越小。将明显QI值小的IMF选出进行去噪处理。

步骤C：对步骤B选出来需要去噪的IMF进行阈值去噪处理；

具体到本实施例，将选出的IMF利用公式(23)对其进行去噪处理。

步骤D：步骤B剩下的IMF与步骤C进行处理后的信号进行信号重构。

通过公式(24)对信号进行重构。

至此，从步骤A到步骤D完成了本实施例一种基于改进EMD与阈值法融合的心电信号去噪方法。

Claims (9)

1.一种基于改进EMD与阈值法融合的心电信号去噪方法，包括步骤：

A.对采集到的信号进行改进的EMD分解，具体地：

A1.对采集到的信号y(t)加入ω_a倍的均值为零、单位方差的正系数白噪声n(t)，得到信号

下标a表示第a次加入白噪声；

A2.用最小二乘支持向量机对信号

两端进行延拓，得到延拓序列；

A3.用保形样条插值的方法构造所述延拓序列的上包络线

及下包络线

j表示第j段T-B样条拟合曲线，所述上包络线

由极大值点、以及在每两个极大值点间插入两个数值点构成T-B样条曲线，所述下包络线

由极小值点、以及在每两个极小值点间插入两个数值点构成T-B样条曲线；其中，所述在每两个极大值点间插入两个数值点具体为：定义所述数值点如下：

其中，t_ij为位于时间t_i-1与t_i之间的第j个插入点，b为极大值点个数，插入数值点后的序列为：

{x_max(t₀),x_max(t₁₁),x_max(t₁₂),x_max(t₁),x_max(t₂₁),x_max(t₂₂),x_max(t₃),...,x_max(t_b)}，

将上述序列作为T-B样条曲线的控制顶点，上包络线即T-B样条曲线表达如下：

其中，j表示拟合的第j段T-B样条曲线；U表示上包络线，h_ij(t)为T-B样条基函数；在每两个极小值点间插入两个数值点插入两个数值点构成T-B样条曲线的构造方法与构造上包络线的方法类似，下包络线的表达如下：

其中，D表示下包络线；

A4.计算第a次加入白噪声后所得到上下包络线的均值

以及信号

与m_a(t)的差值

如果C(t)不满足IMF定义的截止条件，则重复步骤A3；否则提取C(t)作为

表示第a次加入正系数白噪声经第i次EMD分解的IMF；计算剩余量

A5.将步骤A4的剩余量r(t)作为一个新的信号，即r(t)＝y(t)，再经步骤A1～A4进行筛选来获得下一个更低频率的IMF，直到筛选出来的IMF分量满足ZSD＜θ，则停止筛分，其中，θ为设定的阈值，ZSD计算式为：

IMF_i(t)为第i次EMD分解的IMF，n为IMF_i(t)的长度；

A6.经A1～A5分解后的信号表示为

其中，d为心电信号经EMD分解后的IMF的个数，

为第a次加入正系数白噪声经第i次EMD分解的IMF，r_d(t)为进行第d次分解时的残差；

与步骤A1～A5相似，对心电信号加入负系数白噪声后的信号进行EMD分解，分解后的信号表示为

其中，

为第a次加入负系数白噪声经第i次EMD分解的IMF；

A7.对步骤A6得到的加入正系数及负系数白噪声后IMF进行累加求平均，得

B.通过排列互信息判断每个IMF与原信号的线性关系，对由步骤A得到的IMF进行有效信号判定；

C.对步骤B选出来需要去噪的IMF进行阈值去噪处理，所述阈值计算式为：

其中，可调参数m为正数，可调参数

α为正数，

σ为各个IMF分量中噪声标准差，σ＝median(IMF_i(t))/0.6745，median(·)为求中位数；

N为采样点数；

D.步骤B剩下的IMF与步骤C进行处理后的信号进行信号重构。

2.如权利要求1所述的心电信号去噪方法，其特征在于，所述A2包括：对

进行采样得到采样序列为

N为采样点数；通过测试样本集，所述最小二乘支持向量机的输出是

其中，

为延拓后的第一个信号；再将

作为原始数据新的边界点，得到第2个数据序列延拓值

以此类推，根据需要延拓数据点的个数得到全部延拓序列

其中，M为向右延拓的信号点数；对于给定数据序列向左延拓得到

最终延拓序列为

3.如权利要求1或2所述的心电信号去噪方法，其特征在于，所述A3包括：比较所述延拓序列中某点与其左右邻近点的大小关系来判断该点是否为极值点；通过极大值点构造上包络线

通过极小值点构造下包络线

4.如权利要求3所述的心电信号去噪方法，其特征在于，所述判断该点是否为极值点的方法：对于端点

与

若

则点

为极大值点，否则，点

为极小值点；若

则点

为极大值点，否则，点

为极小值点；

对于除上述端点以外的其它点：若

则点

为极大值点；若

则点

为极小值点，其中，-M+1≤i≤N+M-1。

5.如权利要求1所述的心电信号去噪方法，其特征在于，所述步骤B包括：

B1.计算时间序列IMF^(k)(t)和y(t)的排列熵S(IMF^(k))和S(y)，其中IMF^(k)(t)表示第k次EMD分解的IMF时间序列；

B2.计算排列互信息；

B3.根据每个IMF与原信号的互信息QI值判断每个IMF受噪声干扰的程度，将QI值小的IMF选出进行去噪处理。

6.如权利要求1或5所述的心电信号去噪方法，其特征在于，所述排列互信息的计算式为：QI(IMF^(k),y)＝S(IMF^(k))+S(y)-S(IMF^(k),y)，其中，S(IMF^(k))和S(y)分别为IMF^(k)(t)和y(t)的排列熵，S(IMF^(k),y)为IMF^(k)(t)和y(t)的联合排列熵。

7.如权利要求6所述的心电信号去噪方法，其特征在于，所述S(IMF^(k))的排列熵的计算式为：

其中p(·)表示排列(·)的联合概率密度，排列(π_i)表示时间序列IMF^(k)(t)映射到给定嵌入维数空间后每个状态向量所对应的一个排列顺序，d为心电信号经EMD分解后的IMF的个数。

8.如权利要求6所述的心电信号去噪方法，其特征在于，所述IMF^(k)(t)和y(t)的联合排列熵S(IMF^(k),y)的计算式为

其中，排列(π_i,π_j)表示IMF^(k)(t)和y(t)时间序列映射到d维空间后每个状态向量所对应的一个排列顺序，联合概率密度

#表示集合{·}内元素数，下标i'＝1,2,···,n-(d-1)τ，(IMF_i' ^(k),Y_i')是一对状态空间轨迹矩阵，其对应的排列是(π_i,π_j)。

9.如权利要求1所述的心电信号去噪方法，其特征在于，所述步骤D中所述信号重构计算式为

其中，y'(t)表示去噪后的信号。

Images