CN113786191B - 一种基于多通道奇异谱分析的心冲击信号去噪方法 - Google Patents

Abstract

一种基于多通道奇异谱分析的心冲击信号去噪方法属人体生理信号处理技术领域，本发明通过计算STA/LTA比值来获取心跳模板信号，利用模板与心冲击信号的相关性，将心冲击信号分割重建为一个空间结构具备连续性、且信号成分的秩已知的矩阵，使用多通道奇异谱分析分解矩阵，利用信号成分的低秩性与结构连续的相似性，直接选择与秩数相等数量的奇异值重构矩阵，使信号成分与噪声分离，再逆变换为一维心冲击信号。本发明能有效滤除心冲击信号监测系统在采集信号过程中受到的噪声干扰，在有效压制噪声的同时，还可恢复信号中波形失真、缺失的部分，进一步提高信号质量，同时避免重构信号时不易选择哪些奇异值属于信号成分的问题。

Description

一种基于多通道奇异谱分析的心冲击信号去噪方法

技术领域

本发明属于人体生理信号处理技术领域，具体涉及一种基于多通道奇异谱分析的心冲击信号去噪方法。

背景技术

随着国内社会老龄化，对老年人的睡眠健康状态的关注愈发强烈。老年人往往会患有呼吸疾病、心脏疾病、或其它一些慢性病，如呼吸障碍疾病将导致睡眠期间呼吸变浅或停止，容易导致生命危险。因此对一些独居老人来说，如果能在睡眠过程中实现有效的对突发疾病，心率变异等的异常状态检测，将具有十分重要的意义。心冲击信号(BCG)是一种无接触的心脏机械活动的监测方法，能无扰动、非创伤、连续记录人体生理信息，在人体睡眠健康检测方面有广阔的应用前景。在平静的卧姿睡眠状态下，从人体胸部位置采集的BCG信号一般不会包含大的体动伪影与外源干扰，但会包含着呼吸运动，还有源于人体本身机体，例如肌肉或内脏的微弱震动、乃至于人体与衣物或床褥的细小摩擦，都会在信号传递中造成的干扰，导致BCG心跳信号的形变，局部信号难以辨认以至于不能有效地提取心跳，进而造成无法准确的计算心率与疾病检测，因此需要实现BCG信号的噪声压制。

目前，常见的心冲击信号噪声抑制方法有应用数字滤波器进行低、高通滤波过滤噪声；利用小波变换进行多尺度的分析来得到心跳成分；利用经验模态分解、变分模态分解分析，区分信号分量与噪声分量；当采集多道BCG信号时，还可采用独立分量分析来降噪。在上述方法中，如数字滤波器滤波有时会导致滤波结果有些许变形，而且当噪声干扰与心跳信号同频就会存在噪声残留；小波分解、变分模态分解、经验模态分解等方法虽然能够在一定程度上分离开呼吸、心跳与噪声，但需要人为调整小波包，或是分解模态的个数、惩罚项等等参数，使用起来较为复杂，对噪声的压制效果不够理想，而且通常不能恢复变形的信号结构。上述的这些方法仍然是基于一维时序的处理方式，没有考虑到心跳是周期震荡的，心跳信号中各个心跳周期在一定程度上具备相似性，利用这一性质进行心冲击信号的去噪可以取得更佳的效果。

基于上述技术背景，说明已有的去噪方式仍存在缺陷，需要研究并提出更好的BCG信号去噪方法。

发明内容

本发明的目的在于提供一种基于多通道奇异谱分析的心冲击信号去噪方法。将一维的BCG信号重构为一个二维信号矩阵，利用有效心跳信号成分在信号矩阵中的低秩性与心跳信号的相似性，实现对心冲击信号中噪声的压制。

本发明的一种基于多通道奇异谱分析的心冲击信号去噪方法，包括下列步骤：

1.1采集心冲击(Ballistocardiogram，BCG)信号并去除呼吸信号：

用床垫式心冲击信号监测设备采集受测者的心冲击信号，取时长为T的一段BCG信号；

计算BCG信号的均值，并从BCG信号中去除均值以去除直流分量；

利用变分模态分解将去除均值的BCG信号分解为K个模态分量，搜索中心频率处于呼吸信号频带范围内的模态分量，作为呼吸信号分量，将所有呼吸信号分量相加重构得到呼吸信号，从BCG信号中减去呼吸信号得到心跳信号s(t)；

1.2构造心跳信号矩阵，包括下列步骤：

1.2.1计算s(t)的短、长时窗信号幅度平方的均值(STA/LTA)的比值r(t)；

计算r(t)频谱的最大峰值对应的频率，记为峰值频率f，将f的倒数与采样频率相乘得到长度P；

搜索r(t)中所有同时小于前、后相邻采样点幅值的点，作为局部最小点，搜索局部最小点中所有同时小于前、后相邻局部最小点幅值的点，作为各个心跳周期交界点；

以第一个心跳周期交界点为起点，以长度P对心跳信号s(t)平均分割，分割后的所有信号段做算术平均作为心跳模板信号G(t)；

1.2.2将心跳模板信号G(t)与心跳信号s(t)在首端对齐，滑动模板，步长为1，逐个窗口计算模板与心跳信号的相似系数，得到一条相似系数曲线；

相似系数曲线的波动会呈现周期性，与心跳信号的周期性一致，搜索相似系数曲线中所有同时大于前、后相邻相似系数幅值的点，作为局部相似点；

计算所有局部相似点的绝对中位差，如果一个局部相似点的值不大于且不小于3倍的绝对中位差，则保留；

计算所有局部相似点位置序列的差分值，如果一个局部相似点的前、后两个差分值都与长度P相近，则保留；

将保留下的局部相似点作为构造心跳信号矩阵的信号片段标志点，这些标志点会对应信号中的各个心跳周期；

1.2.3以每个构造心跳信号矩阵的信号片段标志点为中心，截取长度为P的心跳信号段，记作心跳信号段{S_i}，i＝1,2…N，同时记录下每条信号段在s(t)中的起终点位置；

将第i个截取出心跳信号段作为心跳信号矩阵S的第i列，把一维时序心跳信号转换为二维心跳信号矩阵S＝[S₁,S₂,…,S_N],S∈R^P×N；

1.3去除心跳信号矩阵中的噪声，包括下列步骤：

1.3.1以P/2为窗长，将S的每条列向量S_i都构造成相应的汉克尔矩阵M_i,

其中P₁＝P/2,N₁＝P-P₁+1，M_i中所有逆对角线上的元素都相等，得到汉克尔矩阵序列M＝(M₁,M₂,M₃…M_N)；

以N/2为窗长，将M构造成块汉克尔矩阵(Hankel Block Matrix)

其中P₂＝N/2,N₂＝N-P₂+1；

1.3.2对块汉克尔矩阵H做奇异值分解，得到按数值大小降序排列的奇异值及相应的特征向量；

选择前k个奇异值及其特征向量作为信号成分，一个干净的心跳信号矩阵由I个复平面波叠加形成,干净心跳信号矩阵对应的块汉克尔矩阵的秩同样是I，噪声会导致块汉克尔矩阵的秩增加，因此选择前k＝I个奇异值及相应的特征向量作为信号成分，剩余的奇异值及相应的特征向量作为噪声成分；

1.3.3选择前k个奇异值及相应的特征向量作为信号成分重构块汉克尔矩阵，得到去后噪的矩阵H′；

对去噪后的矩阵H′依次进行两次逆对角线平均化(汉克尔矩阵逆变换)，得到去噪后的心跳信号矩阵S′；

1.4根据步骤1.2.3中记录的每条信号段的起终点位置，将S′中对应的列向量取出，按位置放还并组成一维信号

若信号段间存在重叠，则对重叠部分做平均处理，得到去噪后的心跳信号

如果一个二维矩阵可认为是由I个复平面波叠加构成，那么其相应的块汉克尔矩阵的秩也会是I。二维矩阵中有效数据成分呈现出低秩性，受到噪声与干扰影响后，会导致块汉克尔矩阵的秩增加，因此可以通过做降秩处理来重构矩阵，以达到去噪的目的。假设一个干净的心跳信号矩阵由I个复平面波叠加形成，那么含噪的矩阵S的相关块汉克尔矩阵H的信号成分秩就应当为I，选择前I个奇异值与相应的特征向量做为信号成分来重构矩阵，即完成了对心跳信号矩阵的去噪。此外由于心跳信号的各个周期是相似的，那么在矩阵S中，尤其是各心跳周期的IJK波部分会主要组成一个横向的W形状波纹，呈现出空间结构的连续相似特征，多通道奇异谱分解可以有效利用这种相似性。

本发明所解决的技术问题：通常获取心跳模板信号需要外源信息辅助标记，如ECG信号，但实际家庭应用中并不适合采取这种接触式测量方式，因此本发明通过计算STA/LTA比值来获取心跳模板信号，能利用BCG信号本身的波动提取各心跳周期的信息；使用多通道奇异谱分析心跳信号矩阵，直接选择与秩数相等数量的奇异值与相应的特征向量重构矩阵，从而区分信号与噪声，避免重构信号时不易选择哪些奇异值属于信号成分的问题；由于心跳的各周期信号是相似的，奇异谱分解得到的是各周期共同的特征，所以在压制噪声的同时，还可以恢复信号中波形失真、缺失的部分，进一步提高信号质量。

本发明的有益效果在于：心冲击信号的采集过程中易受到噪声干扰，使信号的主要形态特征失真，影响心率计算、疾病识别的准确率，而且当前常见的去噪方法对噪声的压制效果不够理想，通常不能恢复变形的信号。针对以上问题，本发明提出一种基于多通道奇异谱分析的心冲击信号去噪方法，考虑到每个周期的心跳信号在一定程度上具备相似性，基于这一特性将每个心跳周期提取出来重建为一个波形空间结构存在连续性，且主要信号成分的秩已知的矩阵，而噪声在时空间上不具备相似性，因此通过多通道奇异谱分解，选择与秩数量相同的奇异值重构矩阵，通过降秩从而达到了去噪的目的。本发明能有效滤除心冲击信号监测系统在采集信号过程中受到的噪声干扰，并且奇异谱分解的出的信号成分包含的是各心跳周期整体共有的信息，相似的成分可以更好的帮助恢复结构失真部分；选择哪些奇异值与特征向量重构信号可做为已知的，不用额外设定选取规则，避免了因难以抉择哪些奇异值代表信号成分，而导致出现噪声残留或信号损失的问题。

附图说明

图1为一种基于多通道奇异谱分析的心冲击信号去噪方法的流程图；

图2为一个周期的理想心冲击信号；

图3为模拟的干净BCG信号；

图4为从实际信号中估计出的噪声；

图5为添加了噪声后的模拟BCG信号；

图6为获得的心跳模板信号；

图7为将一维信号分割并重建二维矩阵的示意图；

图8为模拟含噪数据的信号矩阵；

图9为应用本发明对模拟数据去噪后的结果；

图10为应用SSA对模拟数据去噪后的结果；

图11为一段实测BCG信号；

图12为应用本发明对实测BCG信号去噪后的结果；

图13为应用SSA方法对实测BCG信号去噪后的结果。

具体实施方式

为了使本发明所要解决的技术问题、技术方案及有益效果更加清楚明白，以下结合附图及实施例，对本发明进行详细说明。

本发明的处理流程如图1所示，首先获取BCG心跳模板信号；然后计算模板与BCG信号的相关性来构造信号矩阵，图中黑色圆形标记是构造心跳信号矩阵的信号片段标志点；接着进行信号矩阵的去噪，最后恢复为一维信号。

本发明的一种基于多通道奇异谱分析的心冲击信号去噪方法，包括下列步骤：

1.1采集心冲击(Ballistocardiogram，BCG)信号并去除呼吸信号：

用床垫式心冲击信号监测设备采集受测者的心冲击信号，取时长为T的一段BCG信号；

计算BCG信号的均值，并从BCG信号中去除均值以去除直流分量；

利用变分模态分解将去除均值的BCG信号分解为K个模态分量，搜索中心频率处于呼吸信号频带范围内的模态分量，作为呼吸信号分量，将所有呼吸信号分量相加重构得到呼吸信号，从BCG信号中减去呼吸信号得到心跳信号s(t)；

1.2构造心跳信号矩阵，包括下列步骤：

1.2.1计算s(t)的短、长时窗信号幅度平方的均值(STA/LTA)的比值r(t)；

计算r(t)频谱的最大峰值对应的频率，记为峰值频率f，将f的倒数与采样频率相乘得到长度P；

搜索r(t)中所有同时小于前、后相邻采样点幅值的点，作为局部最小点，搜索局部最小点中所有同时小于前、后相邻局部最小点幅值的点，作为各个心跳周期交界点；

以第一个心跳周期交界点为起点，以长度P对心跳信号s(t)平均分割，分割后的所有信号段做算术平均作为心跳模板信号G(t)；

1.2.2将心跳模板信号G(t)与心跳信号s(t)在首端对齐，滑动模板，步长为1，逐个窗口计算模板与心跳信号的相似系数，得到一条相似系数曲线；

相似系数曲线的波动会呈现周期性，与心跳信号的周期性一致，搜索相似系数曲线中所有同时大于前、后相邻相似系数幅值的点，作为局部相似点；

计算所有局部相似点的绝对中位差，如果一个局部相似点的值不大于且不小于3倍的绝对中位差，则保留；

计算所有局部相似点位置序列的差分值，如果一个局部相似点的前、后两个差分值都与长度P相近，则保留；

将保留下的局部相似点作为构造心跳信号矩阵的信号片段标志点，这些标志点会对应信号中的各个心跳周期；

1.2.3以每个构造心跳信号矩阵的信号片段标志点为中心，截取长度为P的心跳信号段，记作心跳信号段{S_i}，i＝1,2…N，同时记录下每条信号段在s(t)中的起终点位置；

将第i个截取出心跳信号段作为心跳信号矩阵S的第i列，把一维时序心跳信号转换为二维心跳信号矩阵S＝[S₁,S₂,…,S_N],S∈R^P×N；

1.3去除心跳信号矩阵中的噪声，包括下列步骤：

1.3.1以P/2为窗长，将S的每条列向量S_i都构造成相应的汉克尔矩阵M_i,

其中P₁＝P/2,N₁＝P-P₁+1，M_i中所有逆对角线上的元素都相等，得到汉克尔矩阵序列M＝(M₁,M₂,M₃…M_N)；

以N/2为窗长，将M构造成块汉克尔矩阵(Hankel Block Matrix)

其中P₂＝N/2,N₂＝N-P₂+1；

1.3.2对块汉克尔矩阵H做奇异值分解，得到按数值大小降序排列的奇异值及相应的特征向量；

选择前k个奇异值及其特征向量作为信号成分，一个干净的心跳信号矩阵由I个复平面波叠加形成,干净心跳信号矩阵对应的块汉克尔矩阵的秩同样是I，噪声会导致块汉克尔矩阵的秩增加，因此选择前k＝I个奇异值及相应的特征向量作为信号成分，剩余的奇异值及相应的特征向量作为噪声成分；

1.3.3选择前k个奇异值及相应的特征向量作为信号成分重构块汉克尔矩阵，得到去后噪的矩阵H′；

对去噪后的矩阵H′依次进行两次逆对角线平均化(汉克尔矩阵逆变换)，得到去噪后的心跳信号矩阵S′；

1.4根据步骤1.2.3中记录的每条信号段的起终点位置，将S′中对应的列向量取出，按位置放还并组成一维信号

若信号段间存在重叠，则对重叠部分做平均处理，得到去噪后的心跳信号

实施例：

下面通过对模拟数据和实测数据的实验来说明本发明的效果，并结合附图进一步说明本发明的实施流程，以便本领域的技术人员更好地理解本发明。

一个周期的理想心冲击信号形态如图2所示，主要由H、I、J、K、L、M、N等7个波组成，其中IJK峰群幅度普遍较大，构成了BCG信号的主要形态结构，被称为W型组合波，也是本发明要去噪并恢复的主要波形结构。图3所示是基于高斯核函数微分方程模拟出BCG信号。与图2类似，每个心跳都包含了七个子波，每个子波对应一个高斯函数，改变各个高斯函数的幅度，中心坐标位置与方差，即对应着改变每个子波的幅度，相位与持续时间。组合七个子波从而模拟出的一段心率为65次/min的BCG信号。图4是从实际BCG信号中估计出的噪声，可以看到在4s左右噪声幅度较大。将噪声数据添加到模拟BCG信号中，得到含噪的模拟数据如图5所示，信噪比为4.98dB。含噪信号在3-5s的数据上，已经无法识别出心跳规则的W形组合波，波形结构被破坏。

分别采用本发明的方法与一维SSA去噪方法处理模拟含噪数据，来验证本发明的噪声压制能力。对于本发明方法，首先获得心跳的模板信号如图6所示，模板信号依然是含噪的，但从整体趋势能看出符合心跳信号的形态。图7是将一维时序BCG信号重建为信号矩阵的具体过程示意图。图中实线为BCG信号，粗体虚线为心跳模板信号与BCG信号的相关系数曲线，圆圈为构造心跳信号矩阵的信号片段标志点，记录这些点的位置，同时以这些标记点为中心，在前后共取P个采样点对BCG时序信号进行分割，同时将各段作为列向量逐次排列作为信号矩阵。图8是根据模拟含噪数据重构成的信号矩阵，可以看出各心跳尽管受到了噪声干扰，但IJK峰群在仍呈现出一条横向的W状一致性纹路，呈现出空间结构的连续相似性，也证明心跳周期是相似的，而噪声在时空间上没有相似性。图9中虚线是将信号矩阵去噪并重构回一维序列的去噪结果，实线是干净的模拟信号。图10则是利用一维SSA去噪后的结果。在去噪后，SSA方法的信噪比提升到11.9dB，虽然能适当滤除噪声，但是造成部分有效信号的损失，在信号受干扰较严重3-5s处，L、M、N等幅度较小的子波并没有得到良好的恢复。而本发明的结果中，信噪比提升到13.7dB，去噪后曲线变得平滑清晰，3-5s处心跳波形恢复良好，噪声被有效的抑制，证明了本方法的有效性。

采用实测的BCG信号验证本发明的有效性。图11为一段实测BCG信号，信号第4s前后的数据受到干扰，信号结构发生了形变。图12为应用本发明方法的去噪结果，图13是应用SSA的去噪结果。对比实验结果图12和图13，整体来看两种方法都有一定的噪声抑制能力，含噪信号的一些细小毛刺、伪峰得到了抑制。但SSA方法没能恢复第4s前的紊乱波形，而本发明将紊乱心跳的W形状凸显出来，变得容易辨别。可以证明本发明在对心冲击信号进行噪声抑制与结构恢复上可以取得较好的效果。

Claims (1)

1.一种基于多通道奇异谱分析的心冲击信号去噪方法，其特征在于，包括下列步骤：

1.1采集心冲击信号，即采集BCG信号并去除呼吸信号：

用床垫式心冲击信号监测设备采集受测者的心冲击信号，取时长为T的一段BCG信号；

计算BCG信号的均值，并从BCG信号中去除均值以去除直流分量；

利用变分模态分解将去除均值的BCG信号分解为K个模态分量，搜索中心频率处于呼吸信号频带范围内的模态分量，作为呼吸信号分量，将所有呼吸信号分量相加重构得到呼吸信号，从BCG信号中减去呼吸信号得到心跳信号s(t)；

1.2构造心跳信号矩阵，包括下列步骤：

1.2.1计算s(t)的短、长时窗信号幅度平方的均值的比值r(t)；

计算r(t)频谱的最大峰值对应的频率，记为峰值频率f，将f的倒数与采样频率相乘得到长度P；

搜索r(t)中所有同时小于前、后相邻采样点幅值的点，作为局部最小点，搜索局部最小点中所有同时小于前、后相邻局部最小点幅值的点，作为各个心跳周期交界点；

以第一个心跳周期交界点为起点，以长度P对心跳信号s(t)平均分割，分割后的所有信号段做算术平均作为心跳模板信号G(t)；

1.2.2将心跳模板信号G(t)与心跳信号s(t)在首端对齐，滑动模板，步长为1，逐个窗口计算模板与心跳信号的相似系数，得到一条相似系数曲线；

相似系数曲线的波动会呈现周期性，与心跳信号的周期性一致，搜索相似系数曲线中所有同时大于前、后相邻相似系数幅值的点，作为局部相似点；

计算所有局部相似点的绝对中位差，如果一个局部相似点的值不大于且不小于3倍的绝对中位差，则保留；

计算所有局部相似点位置序列的差分值，如果一个局部相似点的前、后两个差分值都与长度P相近，则保留；

将保留下的局部相似点作为构造心跳信号矩阵的信号片段标志点，这些标志点会对应信号中的各个心跳周期；

1.2.3以每个构造心跳信号矩阵的信号片段标志点为中心，截取长度为P的心跳信号段，记作心跳信号段{S_i}，i＝1,2…N，同时记录下每条信号段在s(t)中的起终点位置；

将第i个截取出心跳信号段作为心跳信号矩阵S的第i列，把一维时序心跳信号转换为二维心跳信号矩阵S＝[S₁,S₂,…,S_N],S∈R^P×N；

1.3去除心跳信号矩阵中的噪声，包括下列步骤：

1.3.1以P/2为窗长，将S的每条列向量S_i都构造成相应的汉克尔矩阵M_i,

其中P₁＝P/2,N₁＝P-P₁+1，M_i中所有逆对角线上的元素都相等，得到汉克尔矩阵序列M＝(M₁,M₂,M₃…M_N)；

以N/2为窗长，将M构造成块汉克尔矩阵H，

其中P₂＝N/2,N₂＝N-P₂+1；

1.3.2对块汉克尔矩阵H做奇异值分解，得到按数值大小降序排列的奇异值及相应的特征向量；

选择前k个奇异值及其特征向量作为信号成分，一个干净的心跳信号矩阵由I个复平面波叠加形成,干净心跳信号矩阵对应的块汉克尔矩阵的秩同样是I，噪声会导致块汉克尔矩阵的秩增加，因此选择前k＝I个奇异值及相应的特征向量作为信号成分，剩余的奇异值及相应的特征向量作为噪声成分；

1.3.3选择前k个奇异值及相应的特征向量作为信号成分重构块汉克尔矩阵，得到去后噪的矩阵H′；

对去噪后的矩阵H′依次进行两次逆对角线平均化，得到去噪后的心跳信号矩阵S′；

1.4根据步骤1.2.3中记录的每条信号段的起终点位置，将S′中对应的列向量取出，按位置放还并组成一维信号

若信号段间存在重叠，则对重叠部分做平均处理，得到去噪后的心跳信号

Images