CN110292374B - 基于奇异谱分析和变分模态分解的心电信号去基线漂移方法 - Google Patents

Abstract

本发明公开了一种基于奇异谱分析和变分模态分解的心电信号去基线漂移方法，该方法将原始心电信号经过嵌入和奇异值分解得到奇异值和左右特征向量，通过检测右特征向量功率谱密度的峰值数及峰值所对应的频率将分解后的分量分成三组，对其中含有用信息的两组分量进行重构得到两个信号：第一部分去除基线漂移的心电信号ECG1和心电信号与基线漂移的混叠信号；为去除混叠信号中的基线漂移，采用变分模态分解对混叠信号进一步分解得到第二部分去除基线漂移的心电信号ECG2；最后，ECG1与ECG2相加得到最终去除基线漂移的心电信号。本发明结合奇异谱分析和变分模态分解去除心电信号中的基线漂移干扰结果准确且不依赖于使用者的先验知识。

Description

基于奇异谱分析和变分模态分解的心电信号去基线漂移方法

技术领域

本发明涉及心电信号处理领域，具体涉及一种基于奇异谱分析和变分模态分解的心电信号去基线漂移方法。

背景技术

心电信号是人类最早研究并用于医学临床的生物信号之一，它比其他信号更易于检测，并且具有较为直观的规律性，对心脏基本功能及其病理研究方面具有重要的参考价值。

纯净的心电信号有助于医学诊断分析，然而在心电信号的采集过程中，由于受人体自身因素或医疗器械的影响，不可避免的会带来各种噪声的干扰，如基线漂移、肌电干扰和工频干扰等。其中基线漂移主要是由于人体呼吸运动及采集电极滑动造成的，其特点是：频率超低，通常在0.05Hz～1Hz之间。基线漂移会使心电图整体上移，造成心电轨迹的严重失真，从而影响医学诊断。因此一个稳定有效的噪声去除方法是使用心电信号进行诊断分析的前提，消除基线漂移的影响，在临床医学的诊断中意义重大。

传统用于心电信号去除基线漂移干扰的方法可分为两类，一类是基于傅里叶变换和传统滤波器的线性滤波法，这类方法主要采用特定的函数将时间空间分离，选取固定的频带进行处理，并直接针对信号的全局信息，这种方式需要主观经验去选择频带，对生理信号这类非平稳信号，如何选择合适的频带是一个巨大的难题；另一类是基于非线性滤波结构的方法，这类方法主要包括：中值滤波、自适应滤波和数学形态滤波等，其特点是具有较好的鲁棒性，解决了线性滤波的固定截止频率的问题，但仍无法从理论上合理的给出相关参数的选择方案。

近些年来，小波变换、经验模态分解和变分模态分解被应用于心电信号基线漂移校正中。小波变换的应用难点是小波函数和分解层数的选取，这对处理结果有着较大的影响。经验模态分解法是一种基于数据驱动的自适应分解方法，被用于心电信号的基线漂移去除上，滤波效果有了较大的提升，但仍不可避免的引入了模态混叠的问题。变分模态分解的提出，改善了经验模函数分解中的模态混叠干扰问题，却同样出现了与非线性滤波结构算法的问题，难以从理论上合理的给出相关参数的选择方案。

发明内容

发明目的：为解决现有技术中存在的不足，本发明提供一种基于奇异谱分析和变分模态分解的心电信号去基线漂移方法，在去除基线漂移干扰的同时尽可能保留心电信息。

技术方案：本发明的一种基于奇异谱分析和变分模态分解的心电信号去基线漂移方法，包含以下步骤：

步骤1：获取心电信号数据；

步骤2：利用奇异谱分析法SSA从心电信号数据中提取出两个有用信号分量：第一部分去除基线漂移的心电信号ECG1和心电信号与基线漂移的混叠信号；

步骤3：利用变分模态分解法VMD将步骤2中获得的混叠信号分解成两个信号分量：第二部分去除基线漂移的心电信号ECG2和残余基线漂移；

步骤4：将步骤2得到的ECG1信号和步骤3中获得的ECG2信号相加，得到最终去除基线漂移后的心电信号。

进一步的，所述步骤2中心电信号进行SSA分解的方法按照以下步骤进行：

步骤2-1：嵌入

设给定的心电信号数据为一维时间序列x＝(x₁,x₂,…,x_N)，将其转化为轨迹矩阵X：

其中L为窗口长度，取值方法为：L≥f_s/f_i，f_s为心电信号数据的采样频率，f_i为所要提取信号的感兴趣频率；

步骤2-2：奇异值分解

由轨迹矩阵X生成矩阵S＝XX^T，其中(·)^T为矩阵的转置，计算矩阵S的特征值并按照降序排列：λ₁,λ₂,…,λ_L，则有λ₁≥λ₂≥…≥λ_L≥0，其对应的特征向量为u₁,u₂,…,u_L，令d＝arg max_i{λ_i＞0}，将轨迹矩阵X表示为d个初等矩阵的和：X＝X₁+X₂+…+X_d，其中

称为右特征向量，其中d＝arg max_i{λ_i＞0}中i的取值范围是：1-L，

中i的取值范围是：1-d；

步骤2-3：分组

将d个初等矩阵X_i分成不相交的三组：心电信号分量组、基线漂移分量组和心电信号与基线漂移混叠分量组，具体分组方法为：计算右特征向量ν_i的功率谱密度，利用峰值检测方法计算功率谱密度的峰值个数及其对应的频率，如果峰值个数大于1且对应的峰值频率同时分布于[0,1]和[1,f_s/2]区间，则该右特征向量ν_i所对应的初等矩阵X_i属于混叠分量组，否则，计算功率谱密度的最大峰值所对应的频率f_max，如果f_max≥1Hz，则该右特征向量ν_i所对应的初等矩阵X_i属于心电信号分量组，否则为基线漂移分量组；

将每个组内的初等矩阵分别相加得到三个合成矩阵：心电信号合成矩阵X_ECG、基线漂移合成矩阵X_BW和心电信号与基线漂移混叠合成矩阵X_mixing，轨迹矩阵X表示为：X＝X_ECG+X_BW+X_mixing；

步骤2-4：重构

将步骤2-3中分组得到的心电信号合成矩阵X_ECG和心电信号与基线漂移混叠合成矩阵X_mixing利用对角平均公式法重构为长度为N的一维时间序列，分别为心电信号ECG1和心电信号与基线漂移的混叠信号，其中对角平均法具体为：给定Y为L×K维矩阵，矩阵元素为y_ij，其中1≤i≤L，1≤j≤K，定义L^*＝min(L,K)，K^*＝max(L,K)，N＝L+K-1，

则重构时间序列G＝(g₀,g₁,…,g_N-1)可通过对角平均公式计算获得：

进一步的，所述步骤3中混叠信号进一步进行VMD分解包括如下步骤：

步骤3-1：计算混叠信号的功率谱密度，利用峰值检测方法计算功率谱密度的峰值个数，将其视为VMD分解的最优分解层数P；

步骤3-2：对混叠信号进行VMD分解，得到P个模态分量；

步骤3-3：计算每一个模态分量的功率谱密度，利用峰值检测方法计算功率谱密度的最大峰值所对应的频率f_max，如果f_max≥1Hz，其所对应的模态分量为心电信号分量，否则为基线漂移分量；

步骤3-4：将所有心电信号的模态分量合成为第二部分心电信号ECG2，所有基线漂移的模态分量合成为基线漂移。

进一步的，所述步骤2-3、步骤3-1和步骤3-3中的峰值检测方法为：给定一功率谱密度序列PSD＝(p₁,p₂,…,p_n)，如果p_i-p_i-1＞0，p_i-p_i+1＞0，且p_i＞th，则定义p_i点为峰值，其中i∈[1,n]，th＝min(PSD)+T*[max(PSD)-min(PSD)]，其中T＝0.01。

有益效果：本发明首先利用SSA方法对含基线漂移干扰的心电信号进行初次分解并对分解结果进行分类，去除信号中的基线漂移，提取出含有用信息的主心电信号和心电信号与基线漂移混叠信号，然后利用VMD方法对混叠信号进行二次分解，得到次心电信号和基线漂移，最后对两部分心电信号相加得到最终的去基线漂移结果，本发明具有以下的优点：

(1)改进SSA方法中因窗口长度L选取不当使得分解分量混叠的问题，容忍混叠分量的存在，通过识别算法将混叠分量进行分离并重构成混叠信号，进一步采用VMD方法分解。

(2)改进VMD方法直接用于基线漂移去除时最优分解层数难以选取的问题，本发明在步骤2中利用SSA方法已经将大部分心电信号分量和基线漂移分量进行了分离，在混叠分量中含有的心电信号和基线漂移较少且两者能量相当，便于利用检测功率谱密度的峰值的方法确定VMD的最优分解层数。

(3)提高了心电信号处理的精度，将SSA方法和VMD方法有机结合，在利用SSA去除基线漂移的基础上，进一步利用VMD去除残余的基线漂移，在去除基线漂移的同时尽可能保留心电信号信息。

(4)本发明中SSA方法和VMD方法的相关参数设置是一种数据驱动，自适应的设置方式，因此处理结果的精度不受先验知识的影响。

附图说明

图1为本发明基于奇异谱分析和变分模态分解的心电信号去基线漂移方法的整体流程图；

图2为本发明步骤2的流程图；

图3为本发明步骤3的流程图；

图4为本发明实施例中心电信号时域波形图(来自MIT/BIH心律失常数据库记录105号)；

图5为本发明实施例中加入了真实基线漂移(来自MIT/BIH心电噪声数据库)的心电信号时域波形图；

图6为本发明实施例中对心电信号进行奇异谱分析得到的第一部分心电信号和混叠信号时域波形图；

图7为本发明实施例中对混叠信号检测其功率谱密度函数的峰值个数的示意图；

图8为本发明实施例中对混叠信号进行VMD分解得到的分量及其对应的功率谱密度；

图9为本发明实施例中对混叠信号进行VMD分解得到的第二部分心电信号和残余基线漂移时域波形图；

图10为本发明实施例中去除基线漂移后的心电信号与原心电信号时域波形对比图。

具体实施方式

下面结合附图和实施例对本发明作进一步说明。

参照图1，本实施例所提供的基于奇异谱分析和变分模态分解的心电信号去基线漂移方法，按照以下步骤进行：

首先进入步骤1：获取心电信号数据，本实施例截取了5s来自MIT/BIT心律失常数据库的记录105号心电信号作为“干净”的心电信号x_ECG(t)，其波形如图4所示。为了验证本技术方案的有效性，基线漂移数据x_BW(t)来自MIT/BIT噪声数据库的真实基线漂移数据，同样截取了5s数据。将心电信号与基线漂移叠加得到含基线漂移的心电信号x(t)＝x_ECG(t)+x_BW(t)，即为本实施例中待处理的心电信号，其波形如图5所示。由于MIT/BIT数据库中信号采样频率f_s＝360Hz，故选取的信号采样点数为1800点；

然后进入步骤2：利用奇异谱分析法SSA从心电信号数据中提取出两个有用信号分量：第一部分去除基线漂移的心电信号ECG1和心电信号与基线漂移的混叠信号，参照图2，分解方法如下：

步骤2-1：嵌入

将待处理的1800点的心电信号数据x(t)转化为轨迹矩阵X：

其中L为窗口长度，取值方法为：L≥f_s/f_i，f_s为心电信号数据的采样频率，f_i为所要提取信号的感兴趣频率。

本实施例中，f_s＝360Hz，本发明的目的是去除基线漂移，其频率主要分布在1Hz以内，即感兴趣频率f_i＝1Hz。由L≥f_s/f_i可得窗口长度L≥360，本实施例取L＝361，基于以下两点考虑：一是本技术方案对心电信号首先采用SSA法进行初次分解，之后利用VMD法进行二次分解，在SSA中窗口长度L满足基本条件即可；二是L的增大会增加计算复杂度和计算机的存储，所以本实施例取L＝361；

步骤2-2：奇异值分解

由轨迹矩阵X生成矩阵S＝XX^T，其中(·)^T为矩阵的转置，计算矩阵S的特征值并按照降序排列：λ₁,λ₂,…,λ₃₆₁，则有λ₁≥λ₂≥…≥λ₃₆₁≥0，其对应的特征向量为u₁,u₂,…,u₃₆₁，arg max_i{λ_i＞0}＝361，则轨迹矩阵X可表示为361个初等矩阵的和：X＝X₁+X₂+…+X_d，其中

称为右特征向量，其中d＝arg max_i{λ_i＞0}中i的取值范围是：1-L，

中i的取值范围是：1-d；

步骤2-3：分组

将361个初等矩阵X_i分成不相交的三组：心电信号分量组、基线漂移分量组和心电信号与基线漂移混叠分量组，具体分组方法为：

首先，计算右特征向量ν_i的功率谱密度，采用Welch法估计功率谱密度，Welch法又称平均周期图法，该方法采用信号重叠分段、加窗函数和离散傅里叶变换等计算一个信号序列的自功率谱估计，可以减小估计值的随机起伏，有利于下一步对其进行峰值检测。本实施例对右特征向量ν_i的每个分量按照50％的重叠进行分段，加窗函数为汉明窗，对每段计算其周期图，然后求各个周期图的平均作为功率谱的估值；

其次，根据右特征向量ν_i的功率谱密度进行分组，利用峰值检测方法计算右特征向量ν_i的每个分量的功率谱密度的峰值个数及其对应的频率，本实施例中，右特征向量共361个，其中峰值个数大于1且对应的峰值频率同时分布于[0,1]和[1,f_s/2]区间的向量为第2、4、8、17、23、28、29、31、33、35、45、49、53、59个分量，这14个特征向量所对应的初等矩阵属于混叠分量组，否则，计算功率谱密度的最大峰值所对应的频率f_max，f_max小于1Hz的向量为第1个分量，这些特征向量所对应的初等矩阵属于基线漂移分量组，剩余向量为心电信号分量组。其中，峰值检测方法为：

给定一功率谱密度函数序列PSD＝(p₁,p₂,…,p_n)，如果p_i-p_i-1＞0，p_i-p_i+1＞0，且p_i＞th，则定义p_i点为峰值，其中i∈[1,n]，th＝min(PSD)+T*[max(PSD)-min(PSD)]，其中T＝0.01；

最后，将每个组内的初等矩阵分别相加得到三个合成矩阵：心电信号合成矩阵X_ECG、基线漂移合成矩阵X_BW和心电信号与基线漂移混叠合成矩阵X_mixing，轨迹矩阵X表示为：X＝X_ECG+X_BW+X_mixing；

步骤2-4：重构

将步骤2-3中分组得到的心电信号合成矩阵X_ECG和心电信号与基线漂移混叠合成矩阵X_mixing利用对角平均公式法重构为长度为1800点的一维时间序列，分别为心电信号ECG1和心电信号与基线漂移的混叠信号，其波形如图6所示，图6上半部分为ECG1信号，下半部分为混叠信号，由图可见，含基线漂移的心电信号经过SSA方法的初次分解和分组，已经得到了心电信号的大部分信息，但在混叠分量中仍含有残余的基线漂移干扰。本步骤中对角平均法具体为：

给定Y为L×K维矩阵，矩阵元素为y_ij，其中1≤i≤L，1≤j≤K，定义L^*＝min(L,K)，K^*＝max(L,K)，N＝L+K-1，

则重构时间序列G＝(g₀,g₁,...,g_N-1)可通过对角平均公式计算获得：

之后进入步骤3：利用VMD方法将步骤2中获得的混叠信号分解成两个信号分量：第二部分去除基线漂移的心电信号ECG2和残余基线漂移，参照图3，VMD分解包括以下步骤：

步骤3-1：计算混叠信号的功率谱密度函数，利用峰值检测方法计算功率谱密度函数的峰值个数，峰值检测结果如图7所示，本实施例中，共检测到了7个峰值，说明混叠信号中含有7个主频率成分，故对混叠信号进行VMD分解的最优层数为7。本步骤中功率谱密度计算方法和峰值检测方法与步骤2-3中所述的功率谱密度计算和峰值检测方法一致；

步骤3-2：对混叠信号进行VMD分解，得到7个模态分量，如图8左半部分所示。本实施例中，在应用VMD分解法时，除最优分解层数为7之外，其他参数均按照默认参数设置即可；

步骤3-3：计算每一个模态分量的功率谱密度，如图8右半部分所示。利用峰值检测方法计算功率谱密度的最大峰值所对应的频率f_max，所得结果如表1所示。

表1各模态分量功率谱密度最大峰值对应频率统计

模态1f<sub>max</sub>

模态2f<sub>max</sub>

模态3f<sub>max</sub>

模态4f<sub>max</sub>

模态5f<sub>max</sub>

模态6f<sub>max</sub>

模态7f<sub>max</sub>

0.6Hz

15.6Hz

4.2Hz

7.0Hz

11.4Hz

22.4Hz

2.8Hz

从表1可以看出，第2-7模态分量的功率谱密度最大峰值频率大于等于1Hz，故这6个分量为心电信号分量，第一个分量为基线漂移分量，本步骤中功率谱密度函数计算方法和峰值检测方法与步骤2-3中所述的功率谱密度计算和峰值检测方法一致；

步骤3-4：将所有心电信号的模态分量合成为第二部分心电信号ECG2，如图9上半部分所示，所有基线漂移的模态分量合成为基线漂移，如图9下半部分所示。

最后进入步骤4：将步骤2得到的ECG1信号和步骤3中获得的ECG2信号相加，则得到最终去除基线漂移后的心电信号，其波形如图10所示，实线为去除基线漂移后的心电信号，虚线为原“干净”的心电信号x_ECG(t)，从图中可以明显看出基线漂移信号x_BW(t)基本被完全去除，并完整的保留了真实心电信号x_ECG(t)的细节信息，说明本发明对去除心电信号的基线漂移的有效性。

综上所述，本发明将SSA分解法和VMD分解法进行了有效的结合，在对心电信号进行SSA分解去除基线漂移的基础上，通过VMD分解法进一步去除残余的基线漂移，在保留心电信号信息的同时尽可能的去除基线漂移干扰，提升了后续对心电信号进行分析和诊断的准确性，同时合理的给出了相关参数的选择方案，相比于其他去除基线漂移干扰方法，本发明具有良好的自适应性，处理结果的精度不受先验知识的影响，具有较强的适用性，值得推广。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，根据本发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护范围之内。

Claims (2)

1.一种基于奇异谱分析和变分模态分解的心电信号去基线漂移方法，其特征在于，包括以下步骤：

步骤1：获取心电信号数据；

步骤2：利用奇异谱分析法SSA从心电信号数据中提取出两个有用信号分量：第一部分去除基线漂移的心电信号ECG1和心电信号与基线漂移的混叠信号；

步骤3：利用变分模态分解法VMD将步骤2中获得的混叠信号分解成两个信号分量：第二部分去除基线漂移的心电信号ECG2和残余基线漂移；

步骤4：将步骤2得到的ECG1信号和步骤3中获得的ECG2信号相加，得到最终去除基线漂移后的心电信号；

步骤2中所述的心电信号进行SSA分解的方法按照以下步骤进行：

步骤2-1：嵌入

设给定的心电信号数据为一维时间序列x＝(x₁,x₂,…,x_N)，将其转化为轨迹矩阵X：

其中L为窗口长度，取值方法为：L≥f_s/f_i，f_s为心电信号数据的采样频率，f_i为所要提取信号的感兴趣频率；

步骤2-2：奇异值分解

由轨迹矩阵X生成矩阵S＝XX^T，其中(·)^T为矩阵的转置，计算矩阵S的特征值并按照降序排列：λ₁,λ₂,…,λ_L，则有λ₁≥λ₂≥…≥λ_L≥0，其对应的特征向量为u₁,u₂,…,u_L，令d＝argmax_i{λ_i＞0}，将轨迹矩阵X表示为d个初等矩阵的和：X＝X₁+X₂+…+X_d，其中

称为右特征向量，其中d＝argmax_i{λ_i＞0}中i的取值范围是：1-L，

中i的取值范围是：1-d；

步骤2-3：分组

将d个初等矩阵X_i分成不相交的三组：心电信号分量组、基线漂移分量组和心电信号与基线漂移混叠分量组，具体分组方法为：计算右特征向量ν_i的功率谱密度，利用峰值检测方法计算功率谱密度的峰值个数及其对应的频率，如果峰值个数大于1且对应的峰值频率同时分布于[0,1]和[1,f_s/2]区间，则该右特征向量ν_i所对应的初等矩阵X_i属于混叠分量组，否则，计算功率谱密度的最大峰值所对应的频率f_max，如果f_max≥1Hz，则该右特征向量ν_i所对应的初等矩阵X_i属于心电信号分量组，否则为基线漂移分量组；

将每个组内的初等矩阵分别相加得到三个合成矩阵：心电信号合成矩阵X_ECG、基线漂移合成矩阵X_BW和心电信号与基线漂移混叠合成矩阵X_mixing，轨迹矩阵X表示为：X＝X_ECG+X_BW+X_mixing；

步骤2-4：重构

将步骤2-3中分组得到的心电信号合成矩阵X_ECG和心电信号与基线漂移混叠合成矩阵X_mixing利用对角平均公式法重构为长度为N的一维时间序列，分别为心电信号ECG1和心电信号与基线漂移的混叠信号，其中对角平均法具体为：给定Y为L×K维矩阵，矩阵元素为y_ij，其中1≤i≤L，1≤j≤K，定义L^*＝min(L,K)，K^*＝max(L,K)，N＝L+K-1，

则重构时间序列G＝(g₀,g₁,…,g_N-1)可通过对角平均公式计算获得：

步骤3中所述的混叠信号进一步进行VMD分解包括如下步骤：

步骤3-1：计算混叠信号的功率谱密度，利用峰值检测方法计算功率谱密度的峰值个数，将其视为VMD分解的最优分解层数P；

步骤3-2：对混叠信号进行VMD分解，得到P个模态分量；

步骤3-3：计算每一个模态分量的功率谱密度，利用峰值检测方法计算功率谱密度的最大峰值所对应的频率f_max，如果f_max≥1Hz，其所对应的模态分量为心电信号分量，否则为基线漂移分量；

步骤3-4：将所有心电信号的模态分量合成为第二部分心电信号ECG2，所有基线漂移的模态分量合成为基线漂移。

2.根据权利要求1所述的基于奇异谱分析和变分模态分解的心电信号去基线漂移方法，其特征在于：所述的峰值检测方法为：给定一功率谱密度序列PSD＝(p₁,p₂,…,p_n)，如果p_i-p_i-1＞0，p_i-p_i+1＞0，且p_i＞th，则定义p_i点为峰值，其中i∈[1,n]，th＝min(PSD)+T*[max(PSD)-min(PSD)]，其中T＝0.01。

Images