CN113180680B - 一种改进的基于奇异谱分析的心电信号降噪方法 - Google Patents

Abstract

本发明涉及一种改进的基于奇异谱分析的心电信号降噪方法。具体流程包括：对心电信号选择合适的窗口长度构造轨迹矩阵并进行奇异值分解；计算每一个奇异值所对应的信号重构分量并对第一个重构分量采用低阶多项式拟合的方法去除趋势项；采用陷波器去除各重构分量中的工频干扰；通过计算重构分量间的相关系数的绝对值并与阈值比较选择出有效的信号分量；最后将有效信号分量叠加便可得到降噪后的心电信号。本发明方法由数据驱动，无需参考输入，适用于高斯白噪声、基线漂移、肌电干扰、电力线等噪声的去除，可提高后续的信号分析效果。本方法不仅信号失真程度小，而且在较低信噪比情况下也能较好地降噪，可推广至单路或多路其他生理信号的降噪。

Description

一种改进的基于奇异谱分析的心电信号降噪方法

技术领域

本发明涉及心电信号的降噪处理，具体涉及一种改进的基于奇异谱分析的心电信号降噪方法。

背景技术

体表心电图是通过体表电极和采集设备记录到的标准化导联上的心脏电活动的波形，包含着丰富的生理和病理信息，被认为是人体的重要生理信号。心电图能反映受试者的心脏各部位的生理状况，是评价心脏功能和预测心脏疾病的重要依据之一。

心电信号本身非常微弱，通常情况下最大幅值不超过5mV，在采集过程中极易受到各种干扰，如人体呼吸运动引起的基线漂移、肌肉收缩与颤动引起的肌电干扰、环境噪声、电力线噪声等等。这些干扰和噪声的存在使得信号的波形发生变化，影响研究者对于波形的形态、心跳间隔等特征的获取，导致误诊或者使得信号失去临床价值。因此对心电信号进行降噪的预处理极为重要，有利于为后续的进一步信号分析提供基础，从而降低仪器的误警率。

近年来，研究人员一直致力于开发鲁棒性高的降噪技术，常用的降噪方法包括数字滤波法、自适应滤波、基于小波变换的降噪方法等，都有一定的局限性。如，数字滤波方法无法直接去除频谱与心电信号有重叠的噪声，而且可能改变心电信号的局部特征波甚至引起信号失真；自适应滤波通常需要一个高质量的参考输入信号，这种参考信号通常难以获取；小波基函数的选择对于基于小波变换的降噪方法有很大的影响。

奇异谱分析是一种独立于信号模型由数据驱动的时间序列分解方法，它将时间序列分解为趋势、振荡和噪声分量，它们本质上独立且可加。基于奇异谱分析的方法已在气候学、生物医学图像处理等方面得到有效应用。但在心电信号降噪领域的应用还有待更深入的探索，如：对于基线漂移、工频干扰等噪声难以去除；对于如何选择用于叠加的重构分量的个数也是一个难题。

发明内容

本发明的目的在于提出了一种改进的基于奇异谱分析的心电信号降噪方法。本发明方法能够对与心电信号频谱有重叠的噪声进行去除，且能高效地去除多种噪声包括基线漂移、电力线干扰、肌电干扰等。本发明方法不依赖于函数或者模型，由数据驱动，无需参考输入。该方法不仅不会改变局部特征波，使得信号失真程度小，对有效信息损害小，而且在较低信噪比情况下也能较好地实现降噪。

本发明提出了一种改进的基于奇异谱分析的心电信号降噪方法，具体步骤如下：

(1)读取一段采样频率为f_s、长度为N的单导联心电信号数据，得到一组一维心电信号片段；

(2)对于步骤(1)得到的一维心电信号片段X＝[x₁,x₂,…,x_N]，设定窗口长度为L(L≤N/2)，将数据长度为N的原始单导联心电信号映射成K个L滞后向量，并用这些滞后向量形成一个轨迹矩阵，该轨迹矩阵表示为X_Traj：

其中，K＝N–L+1，且K＞＞L；

对轨迹矩阵X_Traj做奇异值分解，将得到的L个奇异值按大小降序排列，即：

λ₁≥λ₂≥…≥λ_L

轨迹矩阵X_Traj可以写成上述L个奇异值λ₁，λ₂，…，λ_L对应的分量矩阵Z_i相加的形式，i为1，2至L：

(3)对于步骤(2)得到的分量矩阵Z_i，对它们所在逆对角线上的数值累加求平均，得到信号重构分量C_i(i为1，2至L)，可以表示为：

(4)仅针对步骤(3)得到的信号重构分量C_i中的第一个重构分量C₁，用低阶多项式拟合C₁得到f_{_}C₁，再用C₁减去f_{_}C₁，得到去除非线性趋势后的新的第一个重构分量C₁’；

(5)对步骤(4)中得到的新的第一个重构分量C₁’和信号重构分量C_i(i为2，3至L)，采用截止频率为50Hz的陷波滤波器去除所有分量中的工频干扰，得到滤波后的重构分量FC_i(i为1，2至L)；

(6)对经过步骤(5)得到的滤波后的重构分量FC_i，计算它们中任意两个滤波后的重构分量FC_p和FC_q(p,q均为1，2至L)之间的相关系数ρ，得到大小为L×L的相关系数矩阵；其中，相关系数ρ的计算公式为：

设循环次数为w(初始值为1，且w≤L)，分别将第w个滤波后的重构分量FC_w与第一个滤波后的重构分量FC₁、第二个滤波后的重构分量FC₂之间的相关系数的绝对值记为CORR_1w和CORR_2w；然后设定经验阈值为t(0<t<1)，若相关系数CORR_1w＞t或CORR_2w＞t，则w＝w+1，并继续该循环，否则循环结束；当循环结束时，设置用于叠加的重构分量个数r为r＝w–1；

(7)将步骤(6)选择的前r个滤波后的重构分量(即前r个奇异值λ₁，λ₂，…，λ_r所对应的重构分量FC₁，FC₂，…，FC_r)叠加，得到降噪后的心电信号片段Y＝[y₁,y₂,…,y_N]，即：

本发明中，步骤(4)仅针对步骤(3)得到的信号重构分量C_i中的第一个重构分量C₁，采用减去低阶多项式拟合f_{_}C₁的方法得到去除非线性趋势后的新的第一个重构分量C₁’，该步骤可去除信号中的基线漂移噪声。

本发明中，步骤(5)中对步骤(4)中得到的新的第一个重构分量C₁’和信号重构分量C_i(i为2，3至L)，进行截止频率为50Hz的陷波滤波，该步骤可去除信号中的工频干扰。

本发明中，步骤(6)中使用相关系数阈值法来确定要选择的重构分量个数：计算任意两个滤波后的重构分量FC_p和FC_q(p,q均为1，2至L)之间的相关系数矩阵，并计算所有滤波后的重构分量FC_w分别与第一个滤波后的重构分量FC₁、第二个滤波后的重构分量FC₂之间的相关系数的绝对值，将这两个值与设定的阈值进行循环比较(循环次数w的初始值为1，且w≤L)；当且仅当两者都小于等于阈值时，循环结束；将循环结束时的w值减去1，作为后续用于叠加的重构分量的个数。

本发明具有以下有益效果：

1.本发明方法仅由数据驱动，无需依赖模型函数或者其他参考信号，适用于高斯白噪声、基线漂移、肌电干扰、电力线等噪声的去除，可提高后续的信号分析效果；

2.本发明方法通过多项式拟合的方法去除第一个重构分量C₁中的非线性趋势，使得原始心电信号序列中的基线漂移得到有效去除，最后得到的降噪后信号序列中不含或只含非常少量的基线漂移噪声；

3.本发明方法通过让所有重构分量C_i都经过一个截止频率为50Hz的陷波滤波器，尽可能地去除原始信号中的工频干扰，将奇异谱分析方法与滤波方法有效结合，有利于噪声的进一步去除；

4.本发明方法通过滤波后的重构分量之间的相关系数阈值法来实现滤波后的重构分量个数的选择，在低信噪比情况下也能较好地实现分量个数的选择，对提高降噪算法的鲁棒性有一定意义；

5.本发明方法不仅信号失真程度小，而且在较低信噪比情况下也能较好地降噪，可推广至单路或多路其他生理信号的降噪。

附图说明

图1是对心电信号(ECG)应用本发明方法进行降噪的总体流程图。

图2是本实施例中原始干净信号、添加10dB高斯白噪声后的带噪信号、降噪后信号的对比图。其中(a)为10s原始干净信号片段的示意图，(b)为干净信号片段添加10dB高斯白噪声后得到的带噪信号的示意图，(c)为采用本发明方法对带噪信号进行降噪后的信号示意图。

图3是本实施例中带噪信号、降噪后信号与噪声的对比图。其中，(a)为带噪信号与降噪后信号的对比示意图，(b)为本实施例残差(去除掉的噪声)的示意图。

具体实施方式

下面结合附图和实施例对本发明方法及其应用做进一步详细说明。实施例仅用于对本发明进行说明，而非限定条件。在本发明技术方案的基础上，凡根据本发明原理对该实施例的各种变形或修正，均不应排除在本发明的保护范围之外。

实施例1：将本发明的一种改进的基于奇异谱分析的心电信号降噪方法应用于添加10dB高斯白噪声的心电信号(ECG)，其流程如图1所示；本实施例的数据来自公开的数据库MIT-BIH Arrhythmia Database；工作流程如下：

(1)选择MIT-BIH Arrhythmia Database数据库中标号为100的记录中的某10s片段，其采样率为360Hz，使用滤波等方法对其进行预处理，使原始信号尽可能地干净；得到的干净ECG信号如图2(a)所示；对得到的原始干净ECG信号添加10dB的高斯白噪声，添加噪声后得到的带噪信号如图2(b)所示；

(2)对步骤(1)中得到的10s带噪信号片段，设定窗口长度L为90，将数据长度为3600的原始心电信号映射成3511个90-滞后向量，并用这些滞后向量形成一个轨迹矩阵，表示为X_Traj：

对轨迹矩阵X_Traj做奇异值分解，将得到的90个奇异值按大小降序排列，即：

λ₁≥λ₂≥…≥λ₉₀

轨迹矩阵X_Traj可以写成上述90个奇异值λ₁，λ₂，…，λ₉₀对应的分量矩阵Z_i相加的形式，i为1，2至90：

(3)对于步骤(2)得到的分量矩阵Z_i，对它们所在逆对角线上的数值累加求平均，得到信号重构分量C_i(i为1，2至90)，可以表示为：

(4)仅针对对步骤(3)得到的重构分量C_i中的第一个重构分量C₁，用低阶多项式拟合C₁得到f_{_}C₁，再用C₁减去f_{_}C₁，得到去除非线性趋势后的新的重构分量C₁’；

(5)对步骤(4)中得到的新的第一个重构分量C₁’和信号重构分量C_i(i为2，3至90)，采用截止频率为50Hz的陷波滤波器去除所有分量中的工频干扰，得到滤波后的重构分量FC_i(i为1，2至90)；

(6)经过步骤(5)得到的滤波后的重构分量FC_i，计算它们中任意两个滤波后的重构分量FC_p和FC_q(p,q均为1，2至90)之间的相关系数，得到大小为90×90的相关系数矩阵；设循环次数为w(初始值为1，且w≤90)，分别将第w个滤波后的重构分量FC_w与第一个滤波后的重构分量FC₁、第二个滤波后的重构分量FC₂之间的相关系数的绝对值记为CORR_1w和CORR_2w；然后设定经验阈值为t＝0.0015，若相关系数CORR_1w＞t或CORR_2w＞t，则w＝w+1，并继续该循环，否则循环结束；当循环结束时，设置用于叠加的重构分量个数r为r＝w–1；这里得到的r＝14；

(7)将步骤(6)选择的前14个滤波后的重构分量叠加，得到降噪后的心电信号片段Y＝[y₁,y₂,…,y₃₆₀₀]，即：

降噪后心电信号如图2(c)所示；图2展示了原始干净信号、带噪信号、降噪后信号的对比图，可见噪声被明显地去除，且降噪后信号未有明显失真，有效信息保留完整；

表1展示了对原始干净心电信号分别添加10dB高斯白噪声、基线漂移、肌电干扰、电力线噪声后应用本发明方法进行降噪后得到的评价指标结果；

评价指标分别为：

(1)降噪后信噪比SNR_OUT：

(2)均方误差MSE：

(3)相关系数CORR：

表1

10dB噪声	SNR<sub>OUT</sub>	MSE	CORR
				高斯白噪声	17.4148	2.3148e-04	0.9909
基线漂移	20.3121	1.1879e-04	0.9954
				肌电干扰	13.8507	5.2591e-04	0.9804
电力线噪声	33.0667	6.2996e-06	0.9998

Claims (1)

1.一种改进的基于奇异谱分析的心电信号降噪方法，其特征在于具体步骤如下：

(1)读取一段采样频率为f_s、长度为N的单导联心电信号数据，得到一组一维心电信号片段；

(2)对于步骤(1)得到的一维心电信号片段X＝[x₁,x₂,…,x_N]，设定窗口长度为L，L≤N/2，将数据长度为N的原始单导联心电信号映射成K个L滞后向量，并用这些滞后向量形成一个轨迹矩阵，该轨迹矩阵表示为X_Traj：

其中，K＝N–L+1，且K＞＞L；

对轨迹矩阵X_Traj做奇异值分解，将得到的L个奇异值按大小降序排列，即：

λ₁≥λ₂≥…≥λ_L

轨迹矩阵X_Traj写成上述L个奇异值λ₁，λ₂，…，λ_L对应的分量矩阵Z_i相加的形式，其中i为1，2至L：

(3)对于步骤(2)得到的分量矩阵Z_i，对它们所在逆对角线上的数值累加求平均，得到信号重构分量C_i，其中i为1，2至L，表示为：

(4)仅针对步骤(3)得到的信号重构分量C_i中的第一个重构分量C₁，用低阶多项式拟合C₁得到f_C₁，再用C₁减去f_C₁，得到去除非线性趋势后的新的第一个重构分量C₁’；

(5)对步骤(4)中得到的新的第一个重构分量C₁’和信号重构分量C_i，i为2，3至L，采用截止频率为50Hz的陷波滤波器去除所有分量中的工频干扰，得到滤波后的重构分量FC_i，i为1，2至L；

(6)对经过步骤(5)得到的滤波后的重构分量FC_i，计算它们中任意两个滤波后的重构分量FC_p和FC_q之间的相关系数ρ，p,q均为1，2至L，得到大小为L×L的相关系数矩阵；其中，相关系数ρ的计算公式为：

设循环次数为w，初始值为1，且w≤L，分别将第w个滤波后的重构分量FC_w与第一个滤波后的重构分量FC₁、第二个滤波后的重构分量FC₂之间的相关系数的绝对值记为CORR_1w和CORR_2w；然后设定经验阈值为t，0<t<1，若相关系数CORR_1w＞t或CORR_2w＞t，则w＝w+1，并继续该循环，否则循环结束；当循环结束时，设置用于叠加的重构分量个数r为r＝w–1；

(7)将步骤(6)选择的前r个滤波后的重构分量，即前r个奇异值λ₁，λ₂，…，λ_r所对应的重构分量FC₁，FC₂，…，FC_r叠加，得到降噪后的心电信号片段Y＝[y₁,y₂,…,y_N]，即：

Images