基于突变点检测与峰值搜索的心音分割算法

沈 伊,孙 静,杨宏波,王威廉*

(1. 云南大学信息学院,云南 昆明 650500;2. 云南省阜外心血管病医院,云南 昆明 650102)

1 引言

心音是人体重要的生理信号之一,其包含的生理和病理信息是临床心脏疾病诊断重要依据。例如:对先天性心脏病(congenialheartdisease,CHD)的初诊和筛查主要依靠心脏听诊。CHD是造成新生儿和婴幼儿死亡的主要原因,在全国多地均位居新生儿出生缺陷的首位[1]。心脏听诊,对医生的经验和主观性有极大的依赖,具有局限性、不稳定性和不准确性[2]。因此,用现代信息技术分析研究心音和辅助诊断技术具有重要意义。

近年来,国内外学者对基于计算机的心音信号分析已有了显著成果。主流分析流程包括以下三个步骤:预处理、特征提取和分类识别。其中,心音信号的预处理是心音研究的基础和前提,主要包含去噪、包络提取和分段定位三个部分。典型的心音分割研究可归纳为以下几类:参考心电的分割[3,4]、基于时域特征的分割[5-7]、基于神经网络的分割[8]、基于概率模型的分割[9,10]。

Giordano N等[3]和王桢桢[4]通过心电的R波、T波与心音存在的对应关系对心音进行分割,前者对24例正常心音分割准确率达到99.2%,后者对正常和异常心音进行分割,平均准确率达到98.5%。Mecheri等[5]通过计算各心动周期之间的相关系数矩阵进行相似性识别对心音进行分割;Chen等[6]利用改进型希尔伯特包络并通过双门限法对心音信号进行分割;Zhao Zhan等[7]预先估计出心动周期,提出一种基于心动周期估计的心音分割,可以正确识别出80%以上的心动周期。Yao Chen等[8]提出了一个持续时间长短记忆网络结合持续时间特征的方法对心音进行分割,准确率达到96.36%。侯雷静等[9]提出一种面向心音分割的个性化高斯混合建模方法(PMM)对480个心音信号进行周期检测与状态分割,准确率达到93.4%;Oliveira[10]使用逻辑回归函数作为先验概率建立HMM模型进行心音状态分割,平均准确率达93%。

上述研究中,参考心电的心音分割算法是利用心电图(electrocardiogram,ECG)与心音图(phonocardiogram,PCG)的对照关系对心音进行分割。该方法需要对心音、心电进行同步采集和同步处理,采集仪器和采集过程较为复杂,尤其将采集仪器放置在新生儿或婴儿身上时比较困难。此外,在左心室或右心室肥大的儿童中,心脏轴左偏或右偏会导致心电图信号异常,使心音分割复杂化[11]。基于神经网络的心音分割算法主要是将已标记的心音用神经网络进行训练,从而得到心音的状态序列。该方法分割准确率高,但通常需要大量数据、复杂的网络结构以及计算机强大的计算能力,且不利于算法部署到移动设备。基于概率模型的心音分割算法使用心音特征参数作为模型的观测序列,然后对心音状态进行分割。该方法性能很大程度受限于转移概率矩阵参数,在处理不同样本时,分割性能会发生变化,没有较好的普适性。

针对上述问题,考虑到基于时域特征的分割能保留心音时间序列中的峰值信息,具有运算速度快、普适性较好以及算法部署至移动设备较为简单的优势[12],本文提出一种包络结合突变点检测与峰值点搜索的心音分割算法,目的是在保证准确分割的前提下能快速分割心音信号,且具有较强的鲁棒性,能适应不同听诊部位导致的数据差异性问题,便于对心音信号进行特征提取,提高分类识别的准确率。

2 心音组成及算法框架

2.1 心音的组成

心音包括第一心音(S1)、第二心音(S2)、第三心音(S3)和第四心音(S4)四个成分。其中,S3、S4一般不是研究的重点[13]。一个完整的心动周期由S1、心脏舒张期、S2和心脏收缩期构成。通常情况下,心脏舒张期要长于心脏收缩期。心音与心电的对应关系如图1所示,ECG信号R峰的位置对应PCG信号中S1的起始,S2的位置出现在T波末端附近[14]。

图1 心音与心电对应关系图

2.2 算法框架

本文提出的心音分割算法系统框架如图2所示。主要包括心音预处理、突变点检测、峰值搜索和心音分割四个部分。其中,突变点指的是S1、S2与收缩期、舒张期的分界点。

图2 心音分割算法系统框架

算法具体步骤如下:第一步,对心音信号进行预处理。主要包括以下几个过程:将心音信号采样频率统一至1000Hz,之后对心音信号进行去噪和幅值归一化处理。第二步,定位突变点。对心音短时能量与频谱质心进行阈值判断,区分出心音信号振动部分与非振动部分,其分界点即为心音信号突变点。第三步,峰值搜索。根据突变点确定大小区间,用其代替传统初始阈值对心音包络极大值进行前、后向搜索,在滤除伪峰的同时通过峰值点补偿算法提高有效峰值的检出率。第四步,心音分割。利用峰峰值之间的距离定位S1和S2,并以突变点为分割点,将心音信号分割成逐个完整的心动周期。

3 方法

3.1 预处理

心音属于人体微弱信号,由于人体系统和检测环境的复杂性,采集心音时往往会伴随各种环境噪声或电子信号的干扰。去噪的目的是滤除心音数据中的杂音,保留有效和关键的信息。近年来,正交小波分解在主频率范围与噪声频率没有交集的非平稳随机信号的去噪效果上表现出色。原因是正交小波分解算法可将信号分解为多个频段,并根据实际需要对各频段成分进行交换和重构,将有效信号从环境噪声中分离并提升信号有效成分的信噪比[15]。双正交小波是正交小波概念的推广。由于正交小波缺乏对称性,相应的滤波器不能保持线性相位,而双正交小波可解决正交小波线性相位和正交性要求的矛盾[16]。

首先将心音的采样频率统一至1000Hz,降低数据大小,提高运算速度。然后对心音信号进行去噪,选择bior5.5小波将心音信号分为五层,选取适合该段心音的阈值函数得到小波系数,滤除噪声后对信号进行重构,得到去噪后的心音信号。原始信号去噪时域波形如图3所示。图3(a)为未经处理的原始心音,图中A、B为收缩期、舒张期杂音。

图3 小波去噪对比图

由图3(b)可以看出,原始心音信号使用db6小波进行五层分解后仍有噪声与信号重叠部分,而采用双正交bior5.5小波进行五层分解后,心音波形更加清晰,噪声基本滤除,如图3(c)所示。

3.2 突变点检测

根据心音产生机制,S1是由心室收缩开始时二尖瓣和三尖瓣关闭,瓣叶突然紧张引起振动而产生;S2是由心室舒张开始时主动脉瓣和肺动脉瓣突然关闭引起的瓣膜振动所产生[17]。因此相较于舒张期和收缩期,S1和S2的振动更加明显,平均能量更高。在语音信号处理中,短时能量(Short-term energy)可以反映信号的时域特征,频谱质心(Spectral Centroid)可以用来计算频谱的位置[18]。结合两者对心音信号进行分析,可有效的区分出心音信号中的振动部分与非振动部分,寻找出心音信号的突变点。突变点检测流程框图如图4所示。

图4 突变点检测流程

突变点检测算法具体步骤如下:首先,对归一化的去噪心音信号进行分帧处理。计算每帧信号的短时能量和频谱质心。其中,短时能量表示一帧信号样点值的平方和,频谱质心是基于信号能量分布的频率一阶矩,公式如式(1)和式(2)所示。设M1和M2分别为短时能量和频谱质心直方图第一大、第二大局部最若某一段信号的短时能量与频谱质心同时小于其阈值,则认为该段信号没有发生振动。反之,若某一段信号的短时能量与频谱质心同时大于其阈值,则认为该段信号发生了振动。并将同时大于短时能量与频谱质心阈值的信号段(也即两者的交集部分)的左右两个端点定义为该段心音。信号的突变点。突变点位置如图5(c)中“*”号所示。

图5 心音信号短时能量、频谱质心及信号突变点

(1)

式中,E表短时能量,S为帧信号,i为第i帧。

(2)

式中,C表示频谱质心,Xi(k),k=1,…,N是第i帧信号的离散傅里叶变换,N为帧长度。

(3)

其中,W为手动设置的阈值,本文取W=5。

3.3 峰值搜索

目前峰值搜索算法有:简单比较寻峰法、导数寻峰法、协方差寻峰法等。其中,简单比较寻峰法计算简单,适合查找孤立的峰值,但搜索时容易受到连续背景信息和噪声的影响;导数寻峰法适合简单的峰形,但检出率不高且适应能力较差;协方差寻峰法主要优点是分辨重峰的能力强,适合寻找弱峰,但计算过程复杂、运算速度较慢。由于去噪过程无法保证将信号完全从环境杂音中剥离,容易导致伪峰的产生,并且上述算法在包络的寻峰过程中均无法保证伪峰的全部滤除,不利于后续对峰值的定位与心音成分的识别。因此,本文提出一种导数寻峰法与阈值搜索结合的寻峰算法,并采用峰值补偿增加有效峰值的检出,最大程度地减少伪峰对寻峰过程的干扰。具体步骤如下:

1)首先对去噪信号进行维奥拉积分包络的提取,公式如式(4)所示,然后对包络进行归一化处理,公式如式(5)所示。虽然维奥拉积分包络具有简单、处理速度快和实时性强的优势[19],但包络存在毛刺尖峰,因此用三阶Savitzky-Golay(S-G)滤波器对窗口内的信号包络进行三阶多项式拟合,得到平滑后的维奥拉积分包络。

(4)

(5)

式中,Et为心音信号的维奥拉积分包络,E为平滑后的维奥拉积分包络。

2)根据突变点确定峰值搜索的大小区间。公式如式(6)和式(7)所示。

(6)

(7)

式中,TH1为小区间,TH2为大区间。

3)通过二阶导数法求取包络所有极大值并记录局部极大值索引表R1。利用大小区间对R1进行正、反向搜索,并计算正反搜索结果的峰值点距离平均均方差,将平均均方差较小的一个搜索结果记为R2。根据式计算阈值TH3,并将其作为小区间搜索结果的判断阈值,筛选出第三次峰值索引结果集R3。

TH3=Median(Env(R2))*10%

(8)

式中,Env(R2)为结果集R2对应的幅值,Median(Env(R2))表示取结果集R2对应的幅值的中位数。

为了减少有效信息的丢失,提高S1、S2定位识别的准确率,增加伪峰消除的步骤。利用S1、S2的最大持续时间结合其幅值高于收缩期与舒张期幅值的特点对R3结果集进行伪峰消除,增加检出正确率,得到最终定位结果集Rend。峰值点补偿算法如图6所示。

图6 峰值点补偿算法

图7为峰值搜索及补偿算法效果图。图7(a)中由于心音收缩期杂音较强产生了伪峰,经过峰值定位及补偿算法后,增加了有效峰值的检出,如图7(b)所示。

图7 伪峰去除结果图

PCG信号中S1和S2的正确识别定位难度较大的原因有很多:在时域中,其持续时间和形态相似;在频域中,其频谱成分互相重叠(S1从10到140Hz,S2从10到200 Hz)[20]。此外,噪声与正常心音和异常心音在形态学上的相似性使得在时域上识别它们非常困难。通常情况下,心音的舒张期比收缩期长,可以以此作为判定S1、S2的依据。根据突变点作为分割点,可以分割出心音的S1、收缩期、S2与舒张期。图8可从心电与心音内在关系的角度客观评价分割正常心音、异常心音的准确性。虚线表示心音的四个状态,其中,状态“1”表示第一心音(S1)、状态“2”表示收缩期、状态“3”表示第二心音(S2)、状态“4”表示舒张期。从图8可以看出分割的四个状态与心电图之间略有偏差,但由于该偏差在误差范围之内,可忽略不计。

图8 基于本文算法的正常、异常心音分割结果

4 实验结果与分析

4.1 实验环境与数据说明

实验所用电脑CPU型号为AMD Ryzen7 4800HS,RAM为16.0GB,软件环境为MATLAB2019b。本文实验数据从PhysioNet/Computing in Cardioloy Challenge选取了1000例心音信号,并通过同步采集的ECG信号对S1、S2的位置进行了标记,标记结果通过多种R峰与T波探测器综合给出[21]。所选取的心音中包含539例正常心音,461例异常心音。

4.2 评价指标

为评估该心音分割算法性能以及其在测试集中准确定位S1、S2的能力,本文采用ECG-PCG对照法对其进行评估。由于心音分割不是分类问题,因此不存在真正的负样本,因此:

1)若算法定位出的S1起始位置在ECG R峰的100ms内,则将其认为分段准确(TP)。

2)若算法定位出的S2中心位置在ECG T波末端的100ms内,则将其认为分段准确(TP)。

3)其余情况均认为是错误分段(FP)。

使用Se：灵敏度,ACC：准确率,F1：评估精度分数(综合评价指标),SN：每分钟内分割心音的个数,四种评估标准,分别由式(9)-(11)定义。

(9)

(10)

(11)

4.3 实验结果

对正常心音、异常心音分割的性能比较如表1所示。本文算法分割性能与引言中提到的其中几种分割算法进行对比的结果如表2所示。对比试验计算SN指标时,统一采用4.5s的心音段,计算一分钟内可分割心音的总个数。

表1 正常、异常心音分割的性能比较

表2 5种分割方法的性能比较

4.4 实验结果分析

通过表1中的准确率评价指标可看出,本文算法对于正常心音的分割准确率较高于异常心音。通过分析,这是由于异常心音中病理性特征的存在,往往使异常心音更为复杂,导致对异常心音的突变点检测发生错误、伪峰去除不完整或S1、S2难以定位。

通过表2可看出,本文算法在准确度、精度和F1指标要优于方法1-3。实验表明,所提算法通过心音振动段与非振动段进行突变点检测并根据突变点计算自适应阈值,对心音更具有普适性,并用峰值搜索对伪峰进行了滤除并进行了峰值补偿操作,增加了有效峰值的检出率,使分割准确率有了提升。通过SN指标可看出,虽然方法4具有较高的准确率和F1,但由于其包含高斯混合建模、前后向运算和修正误差等计算过程,导致该分割法相较于其它分割算法计算过程复杂,耗时较长。由于方法2需要先进行经验模态分解(empirical mode decomposition,EMD)找到固有模态函数(intrinsicmodefunction,IMF)并进行插值运算,因此其运算复杂度仅次于方法4。本文算法在寻找到心音突变点和心音有效峰值后可根据心音时域特征直接进行分割,复杂度远远低于以上两种算法。且相较于方法4,本文算法更容易部署到移动端,更有利于先心病机器辅助诊断系统的普及。

5 结束语

基于时域特征的心音分割方法算法具有计算过程简单、运算速度快的特点。其中突变点检测与峰值定位的过程为关键步骤。本文提出一种利用短时能量和频谱质心结合的突变点检测算法和基于维奥拉积分包络的自适应阈值寻峰算法,根据心音信号时域特点进行峰值筛查与峰值补偿,识别并去除伪峰,从而得到正确的S1、S2定位识别,实现了心音信号的快速分割。相比个性化高斯混合建模的包络分割法,本文算法在保证准确率的情况下更加快速且具有可移植性。相较于传统阈值选取分段算法,本文算法解决了其适用性低、自适应弱和依赖调节系数等问题,采用大小区间搜索与峰值补偿算法得到有效的心音峰值点,保证了分段结果的准确性。但本文算法对于环境噪声过大且无法将心音从噪声剥离的信号分割效果尚不理想,有待进一步研究,以达到更好的分割效果。