CN115105088A - 一种改进的基于小波域稀疏特性的心电信号去噪方法 - Google Patents

Abstract

一种改进的基于小波域稀疏特性的心电信号去噪方法，通过对所有不同尺度的小波系数进行优化计算，避免了小波阈值去噪过程中将较小的小波系数舍弃的问题，从而，能够更好地保留心电信号的细节特征。能够提高降噪后心电信号的组稀疏特性，并减轻基于全变分方法求解过程中出现的阶梯伪影，保持了原始心电信号的波形特征。充分利用了心电信号在时域和频域的组稀疏特性，且该方法不仅适用于非重叠的组稀疏信号，当信号相邻分组出现重叠时，优化变量之间是相互耦合的，该方法仍然适用。通过选取参数化的非凸惩罚函数，并限定正则化参数的区间，保证了总代价函数的严格凸性，且该优化算法计算效率高、收敛速度快。

Description

一种改进的基于小波域稀疏特性的心电信号去噪方法

技术领域

本发明涉及心电信号降噪技术领域，具体涉及一种改进的基于小波域稀疏特性的心电信号去噪方。

背景技术

心电图(Electrocardiogram，ECG)是人体最重要的生物电信号之一，它从各个方面反映了心脏电活动等可靠特征，因此，它是诊断和治疗心血管疾病的无创且有效的临床工具。然而，在心电信号的采集过程中，常常会受到各种噪声和伪影的干扰，这将对后续的心电信号的特征提取、心拍分类和心律失常检测等产生重要影响。在这些噪声中，基线漂移和工频干扰噪声集中在一些较窄的频段，通过选择适当的滤波器可以很容易地去除这些噪声。肌电干扰噪声频谱较宽，其频谱特性与白噪声相似，且在0.01Hz～100Hz范围内与心电信号频谱重叠。由于频谱特性的重叠性，在去除该类噪声的同时，很容易破坏心电信号的重要波形特性，从而，降低医生诊断的准确性。

心电信号的去噪算法有维纳滤波、经验模态分解、小波变换和基于稀疏特性的去噪方法。维纳滤波器是最小均方误差意义上的最优线性滤波方法，可应用于时域的因果系统，或者频域的非因果系统，然而，由于心电信号的非平稳特性，维纳滤波不能很好地抑制心电信号中的噪声。经验模态分解方法去噪时通常会丢弃一些包含噪声的初始固有模态，这将导致重构心电信号的波形失真，特别是QRS波。基于小波变换的去噪方法是目前常用的具有多分辨率特性的信号分析方法，它将噪声信号分解成若干级，根据小波阈值技术去除不同尺度的噪声成分，然后通过小波逆变换进行无噪声信号重构。然而，由于小波阈值函数的不连续性，小波阈值去噪方法往往会引入新的噪声尖峰。基于稀疏特性的去噪方法，如全变分方法能够保持信号的边缘信息，但是，基于全变分的去噪方法常常会导致信号在不连续点处的峰值欠估计，如阶梯状的伪影等。

发明内容

本发明为了克服以上技术的不足，提供了一种改进的基于小波域稀疏特性的心电信号去噪方，通过融合非凸的正则化和改进的全变分思想，实现了不同尺度小波系数的独立优化，在去噪的同时能够更好地保留心电信号的波形特征。

本发明克服其技术问题所采用的技术方案是：

一种改进的基于小波域稀疏特性的心电信号去噪方法，包括如下步骤：

a)建立如y(m)＝x(m)+ε(m)

的心电信号数学模型，其中m表示第m个采样点，

N为采样点总数，y为含有噪声的心电信号，y∈R^N，R^N为N维的实数空间，x为干净的心电信号，x∈R^N，ε为高斯白噪声，ε∈R^N，对干净的心电信号x进行小波变换，得到小波系数w；

b)建立关于小波系数w的凸优化问题，通过公式

计算得到凸优化问题的最优解w^*，式中||·||₂为L2范数，F(w)为非凸成本函数，arg表示F(w)的最小值对应的变量，W为小波变换矩阵，D为一阶微分矩阵，D∈R^(N-1)×N，

λ_j、α均为正则化参数，a_j为参数，λ_j＞0、α＞0，

为参数化的反正切惩罚函数，

atan(·)为反正切函数，

为w_j中起始下标为n且由K_j个连续点组成的向量，

M_j为w_j的长度，K_j为w_j的分组大小，w_j为第j个尺度的小波系数，j＝1,2,...,J，J为小波尺度总数，

为向量DW^Tw起始下标为l且由

个连续点组成的向量，0≤l≤N-2，

为信号

的分组大小，

为通过降噪得到的干净的心电信号，

T为转置；c)基于交替方向乘子方法令u＝w，u为辅助变量，得到关于变量u和w的凸优化问题；

d)通过公式

迭代求解最小化问题，得到凸优化问题的最优解w^*，式中s为变量，s维度与w维度相同，s＝(Wy+ρ(u-d))/(1+ρ)，ρ为常数，ρ＞0，d为优化变量，

变量s的第j个子向量，通过公式

计算得到第k+1步迭代的小波系数

为

的第n个元素，

为w^(k+1)的第j个尺度的分量，式中

为v_j中起始下标为n-i且由K_j个连续点组成的向量，v_j为辅助变量，

为

的一阶导数，

e)通过公式

迭代求解最小化问题，得到凸优化问题的最优解u^*，f(·)为非线性函数，prox_f(W^T(w+d))为邻近算子，

z为辅助变量，g(z)的替代函数为G(z,p)，

式中

为与z无关的常数，

ξ＝Dp，

p为辅助变量，[Γ(ξ)]_l,l为Γ(ξ)第l行第l列的元素，当z≠p时，G(z,p)≥g(z)，当z＝p时，G(p,p)＝g(z)，

f)通过公式

计算得到辅助变量z第k+1步迭代的值z^{(k +1)}，

I为单位矩阵，δ＝Dz，δ^(k)为δ第k步的迭代值，k∈{1,2,...,K₀}，K₀为迭代步数的最大值，

为Γ(δ^(k))第l行第l列的元素，d(^k)为优化变量d的第k步迭代的值；

g)通过公式u^(k+1)＝w^(k+1)+d^(k)+W(z^(k+1)-W^T(w^(k+1)+d^(k)))计算得到变量u的第k+1步迭代的值u^(k+1)；

h)通过公式d^(k+1)＝d^(k)-(u^(k+1)-w^(k+1))计算得到优化变量d的第k+1步迭代的值；

i)给定常数c₀＞0，如果满足收敛条件F(w^(k+1))＜c₀，则执行步骤j)，如果不满足收敛条件则重复执行步骤d)至步骤h)直至满足收敛条件；

j)通过公式

计算得到降噪后的干净的心电信号

进一步的，步骤a)中干净的心电信号x的小波系数为w，w＝Wx，W为小波变换矩阵，W满足W^TW＝I，I为单位矩阵，T为矩阵转置，w_j,t为w的元素，j为小波尺度，t为时间。

进一步的，步骤c)包括如下步骤：

c-1)通过公式

建立优化问题，式中

本发明的有益效果是：通过对所有不同尺度的小波系数进行优化计算，避免了小波阈值去噪过程中将较小的小波系数舍弃的问题，从而，能够更好地保留心电信号的细节特征。能够提高降噪后心电信号的组稀疏特性，并减轻基于全变分方法求解过程中出现的阶梯伪影，保持了原始心电信号的波形特征。充分利用了心电信号在时域和频域的组稀疏特性，且该方法不仅适用于非重叠的组稀疏信号，当信号相邻分组出现重叠时，优化变量之间是相互耦合的，该方法仍然适用。通过选取参数化的非凸惩罚函数，并限定正则化参数的区间，保证了总代价函数的严格凸性，通过凸优化方法能够得到去噪问题的唯一解，且该优化算法计算效率高、收敛速度快。

具体实施方式

下面对本发明做进一步说明。

一种改进的基于小波域稀疏特性的心电信号去噪方法，包括如下步骤：

a)建立如y(m)＝x(m)+ε(m)

的心电信号数学模型，其中m表示第m个采样点，

N为采样点总数，y为含有噪声的心电信号，y∈R^N，R^N为N维的实数空间，x为干净的心电信号，x∈R^N，ε为高斯白噪声，ε∈R^N，对干净的心电信号x进行小波变换，得到小波系数w。心电信号去噪的目的是充分利用心电信号在小波域的组稀疏特性，通过优化方法从含有噪声的信号y中恢复出干净的心电信号x。

b)建立关于小波系数w的凸优化问题，通过公式

计算得到凸优化问题的最优解w^*，式中||·||₂为L2范数，F(w)为非凸成本函数，arg表示F(w)的最小值对应的变量，W为小波变换矩阵，D为一阶微分矩阵，

λ_j、α均为正则化参数，a_j为参数，

为参数化的反正切惩罚函数，

atan(·)为反正切函数，

为w_j中起始下标为n且由K_j个连续点组成的向量，

M_j为w_j的长度，K_j为w_j的分组大小，w_j为第j个尺度的小波系数，j＝1,2,...,J，J为小波尺度总数，

为向量DW^Tw起始下标为l且由

个连续点组成的向量，0≤l≤N-2，

为信号

的分组大小，

为通过降噪得到的干净的心电信号，

T为转置。

c)基于交替方向乘子方法令u＝w，u为辅助变量，得到关于变量u和w的凸优化问题。

d)选取增广拉格朗日函数为：

通过交替求解关于w和u的最小化问题，则关于变量u和w的凸优化问题可以通过如下的迭代算法求解得到，具体为：通过公式

迭代求解最小化问题，得到凸优化问题的最优解w^*，合并二次项之后，w^*可通过迭代求解如下的最小化问题得到

式中s为变量，s维度与w维度相同，s＝(Wy+ρ(u-d))/(1+ρ)，ρ为常数，ρ＞0，d为优化变量。进一步H(w)可表示为

s_j为变量s的第j个子向量。选取如下的替代函数

这里v_j为辅助变量，

的一阶导数，

为向量v_j中起始下标为n且由K_j个连续点组成的向量，C₁表示与w_j无关的常数，Q_j(w_j,v_j)满足条件，

根据优化最小算法，可以得到k+1步的迭代值

j＝1,2,...,J，

这里，

表示w_j(n)第k步的迭代值。具体的：通过公式

计算得到第k+1步迭代的小波系数

为

的第n个元素，

为w^(k+1)的第j个尺度的分量，式中

为v_j中起始下标为n-i且由K_j个连续点组成的向量，v_j为辅助变量，

为

的一阶导数，

e)求解最小化问题

选取非线性函数

可以得到

由于小波变换W是正交的，根据邻近算子的正交线性变换性质，可以通过公式

迭代求解最小化问题，得到凸优化问题的最优解u^*，f(·)为非线性函数，prox_f(W^T(w+d))为邻近算子，

z为辅助变量，g(z)的替代函数为G(z,p)，

式中

为与z无关的常数，

ξ＝Dp，

p为辅助变量，[Γ(ξ)]_l,l为Γ(ξ)第l行第l列的元素，当z≠p时，G(z,p)≥g(z)，当z＝p时，G(p,p)＝g(z)，

f)根据优化最小化算法求解步骤e)中的优化问题，即邻近算子prox_f(W^T(w+d))，具体的通过公式

计算得到辅助变量z第k+1步迭代的值z^(k+1)，

I为单位矩阵，δ＝Dz，δ^(k)为δ第k步的迭代值，k∈{1,2,...,K₀}，K₀为迭代步数的最大值，

第l行第l列的元素，d(^k)为优化变量d的第k步迭代的值。

g)根据矩阵逆引理，通过公式u^(k+1)＝w^(k+1)+d^(k)+W(z^(k+1)-W^T(w^(k+1)+d^(k)))计算得到变量u的第k+1步迭代的值u^(k+1)。

h)通过公式d^(k+1)＝d^(k)-(u^(k+1)-w^(k+1))计算得到优化变量d的第k+1步迭代的值。

i)给定常数c₀＞0，如果满足收敛条件F(w^(k+1))＜c₀，则执行步骤j)，如果不满足收敛条件则重复执行步骤d)至步骤h)直至满足收敛条件。这里F(w^(k+1))为与步骤b)相同的参数化的成本函数，具体的，

j)通过公式

计算得到降噪后的干净的心电信号

从而实现了心电信号去噪的目的。

与传统的小波阈值方法不同，本发明通过对所有不同尺度的小波系数进行优化计算，避免了小波阈值去噪过程中将较小的小波系数舍弃的问题，从而，能够更好地保留心电信号的细节特征。与标准的凸优化方法相比，本发明能够提高降噪后心电信号的组稀疏特性，并减轻基于全变分方法求解过程中出现的阶梯伪影，保持了原始心电信号的波形特征。本发明充分利用了心电信号在时域和频域的组稀疏特性，且该方法不仅适用于非重叠的组稀疏信号，当信号相邻分组出现重叠时，优化变量之间是相互耦合的，该方法仍然适用。通过选取参数化的非凸惩罚函数，并限定正则化参数的区间，保证了总代价函数的严格凸性，通过凸优化方法能够得到去噪问题的唯一解，且该优化算法计算效率高、收敛速度快。

实施例1：

步骤a)中干净的心电信号x的小波系数为w，w＝Wx，W为小波变换矩阵，W满足W^TW＝I，I为单位矩阵，T为矩阵转置，w_j,t为w的元素，j为小波尺度，t为时间。

实施例2：

步骤c)包括如下步骤：

c-1)通过公式

建立优化问题，式中

最后应说明的是：以上所述仅为本发明的优选实施例而已，并不用于限制本发明，尽管参照前述实施例对本发明进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (3)

1.一种改进的基于小波域稀疏特性的心电信号去噪方法，其特征在于，包括如下步骤：

a)建立如

的心电信号数学模型，其中m表示第m个采样点，

N为采样点总数，y为含有噪声的心电信号，y∈R^N，R^N为N维的实数空间，x为干净的心电信号，x∈R^N，ε为高斯白噪声，ε∈R^N，对干净的心电信号x进行小波变换，得到小波系数w；

b)建立关于小波系数w的凸优化问题，通过公式

计算得到凸优化问题的最优解w^*，式中||·||₂为L2范数，F(w)为非凸成本函数，arg表示F(w)的最小值对应的变量，W为小波变换矩阵，D为一阶微分矩阵，D∈R^(N-1)×N，

λ_j、α均为正则化参数，a_j为参数，λ_j＞0、α＞0，

为参数化的反正切惩罚函数，

atan(·)为反正切函数，

为w_j中起始下标为n且由K_j个连续点组成的向量，

M_j为w_j的长度，K_j为w_j的分组大小，w_j为第j个尺度的小波系数，j＝1,2,...,J，J为小波尺度总数，

为向量DW^Tw起始下标为l且由

个连续点组成的向量，0≤l≤N-2，

为信号

的分组大小，

为通过降噪得到的干净的心电信号，

T为转置；

c)基于交替方向乘子方法令u＝w，u为辅助变量，得到关于变量u和w的凸优化问题；

d)通过公式

迭代求解最小化问题，得到凸优化问题的最优解w^*，式中s为变量，s维度与w维度相同，s＝(Wy+ρ(u-d))/(1+ρ)，ρ为常数，ρ＞0，d为优化变量，

s_j为变量s的第j个子向量，通过公式

计算得到第k+1步迭代的小波系数

为

的第n个元素，

为w^(k+1)的第j个尺度的分量，式中

为v_j中起始下标为n-i且由K_j个连续点组成的向量，v_j为辅助变量，

为

的一阶导数，

e)通过公式

迭代求解最小化问题，得到凸优化问题的最优解u^*，f(·)为非线性函数，prox_f(W^T(w+d))为邻近算子，

z为辅助变量，g(z)的替代函数为G(z,p)，

式中

为与z无关的常数，

ξ＝Dp，

p为辅助变量，[Γ(ξ)]_l,l为Γ(ξ)第l行第l列的元素，当z≠p时，G(z,p)≥g(z)，当z＝p时，G(p,p)＝g(z)，

f)通过公式

计算得到辅助变量z第k+1步迭代的值z^(k+1)，

I为单位矩阵，δ＝Dz，δ^(k)为δ第k步的迭代值，k∈{1,2,...,K₀}，K₀为迭代步数的最大值，

[Γ(δ^(k))]_l,l为Γ(δ^(k))第l行第l列的元素，d(^k)为优化变量d的第k步迭代的值；

g)通过公式u^(k+1)＝w^(k+1)+d^(k)+W(z^(k+1)-W^T(w^(k+1)+d^(k)))计算得到变量u的第k+1步迭代的值u(^k+1)；

h)通过公式d^(k+1)＝d^(k)-(u^(k+1)-w^(k+1))计算得到优化变量d的第k+1步迭代的值；

i)给定常数c₀＞0，如果满足收敛条件F(w^(k+1))＜c₀，则执行步骤j)，如果不满足收敛条件则重复执行步骤d)至步骤h)直至满足收敛条件；

j)通过公式

计算得到降噪后的干净的心电信号

2.根据权利要求1所述的改进的基于小波域稀疏特性的心电信号去噪方法，其特征在于：步骤a)中干净的心电信号x的小波系数为w，w＝Wx，W为小波变换矩阵，W满足W^TW＝I，I为单位矩阵，T为矩阵转置，w_j,t为w的元素，j为小波尺度，t为时间。

3.根据权利要求1所述的改进的基于小波域稀疏特性的心电信号去噪方法，其特征在于，步骤c)包括如下步骤：

c-1)通过公式

建立优化问题，式中