基于非平稳系统辨识的心音包络自适应分割

许春冬,周 静,应冬文,2,龙清华

(1.江西理工大学 信息工程学院,江西 赣州 341000;2.中国科学院声学研究所 语言声学与内容理解重点实验室,北京 100190)

0 概述

目前,心血管疾病已成为危害人类生命健康的严重疾病,心血管疾病死亡人数在全球疾病死亡人数中占据较大比例[1]。心音信号(Heart Sound Signal,HSS)分析是心血管疾病常用的辅助诊断方法之一,主要通过计算机辅助诊断分析系统的形式进行应用。心音信号的辅助诊断系统流程可划分为分割与分类两个主要步骤[2-3],其中分割是心音信号分析与诊断的基础与前提,分割结果对整个系统的分析结果具有重要影响。然而,心音信号相当敏感,极易受到仪器噪声、环境噪声、肺音、呼吸道音、肠鸣音及检测者自身活动状态等因素的影响。此外,各类病理原因也会导致心音中出现心杂音[2-3]。这些干扰因素的存在使得分割任务变得十分困难,且实际中采集到的心音信号通常包含多种噪声,导致其可分析性急剧下降,使得原本较弱的心音信号变成一个较复杂的声学信号。因此,心音信号的有效分割逐渐成为目前研究的热点之一。

随着心音信号研究的深入,各种分割方法相继被提出,主流的心音信号分割方法可分为[4-6]基于包络的方法、基于特征提取的方法、基于隐马尔科夫模型(Hidden Markov Model,HMM)的方法和基于机器学习的方法4类。其中,基于包络的方法因算法简单、实时性强等特点[6-7]得到广泛研究与应用,特别是在实时检测需求下成为首选方法[8],而其他3类分割方法的实时性相对较差[5,9-10]。本文借鉴包络分割方法的思想,提出一种基于非平稳性系统辨识的心音包络提取方法,并将其应用于心音阈值自适应分割中。

1 基于包络的分割方法

心音信号分割的首要任务是定位分割出基础心音(Fundamental Heart Sound,FHS)信号s1(第一心音)和s2(第二心音),然后分割出收缩期(sys)和舒张期(dia),如图1所示。基于包络的分割方法主要包括预处理、包络提取、阈值分割3个核心部分。

图1 心音时域波形图Fig.1 Oscillogram of heart sound time domain

预处理主要包括预加重、消除趋势项和降噪,其中:预加重为高频提升,通过一阶FIR滤波器即可完成[2];对于趋势项,可采用最小二乘法拟合进行消除[11];降噪中应用较普遍的方法为基于小波变换的降噪方法,通过丢弃带外小波系数、滤波带内小波系数,实现噪声抑制[12]。预处理后的心音信号特征更清晰,为包络提取奠定了基础。

心音信号包络的提取方法主要包括香农能量包络法[3,5-6]、维奥拉积分包络法[5]、希尔伯特变换包络法[7,9]、样条插值包络法[6]以及经验模态分解(Empirical Mode Decomposition,EMD)包络法[7,10]等。香农能量包络法在噪声干扰较低的情况下能够取得较好的效果,但若降噪后的噪声干扰依旧较大,则提取的能量包络特征会受到严重干扰。维奥拉积分包络法能够较好地提取出s1的特征,但对衰弱s2及心杂音的干扰信号包络提取效果明显下降。在希尔伯特变换法中,当属黄变换研究最深入,但该方法在通过极值点拟合包络时通常存在一些固有缺陷,尤其是对于一些异常心音信号。样条插值包络法通过对极大值和极小值进行插值获得上包络和下包络,可以直观反映信号的幅值,但存在欠包络的问题,且实际中欠包络问题处理结果并不理想。EMD包络法无需预先设定基函数,自适应能力较强,但在处理非平稳信号时存在模态混叠及信号失真等问题。

分割过程主要包括阈值函数设计与分割成分判定两部分。其中,阈值函数多数为双阈值函数[7];任何一个心动周期中均包含第一心音、收缩期、第二心音、舒张期4个部分,因此通过阈值分割判定基础心音后,根据时域特征即可区分出收缩期和舒张期[6-7,9]。

2 非平稳性系统辨识原理

系统辨识主要是根据输入与输出时间序列来描述系统行为的过程,用于估计与建立控制领域内的系统模型[13],其采用的非平稳性系统辨识模型衍生于波束形成等信号处理领域[13-14]。文献[14]提出一种估计传递函数(Transfer Functions,TFs)比率的方法,用于衡量不同信道间的传递差异。首先,假设传递函数比所期望信号变化的更加缓慢,进一步假设噪声信号与期望信号的幅值比是缓慢变化,将时间轴划分为一系列分析间隔,在每一个分析间隔期间内假定TFs和噪声信号均近似为平稳。然后,将该分析时间间隔划分为帧,使得期望信号在每帧内看作不变。以第K帧的信号为例,得到:

(1)

(2)

(3)

(4)

由于Um(t,ejw)和Z1(t,ejw)一般是相关的,因此无法通过计算式(4)来直接获得无偏估计。文献[15]指出采用最小二乘法可求得超定方程组Hm(ejw)(m=2,3,…,M为一组独立方程)的无偏估计。方程式(5)的解为方程式(6)。

(5)

(6)

其中:〈〉表示取均值运算;Hm(ejw)为传递函数比率的估计结果,可用于衡量两个通道信号间的相关性,若两个通道信号相关性强,则Hm(ejw)值较大;若两个通道信号相关性弱,则Hm(ejw)值较小。

3 基于非平稳系统辨识的心音分割

3.1 基于系统辨识的包络提取

虽然心音信号近似为准周期性信号,但其本质上是一类非平稳信号,具备混沌特性[16]。为较好地辨识出基础心音信号,需要提取更有效的特征包络,突出基础心音特征,因此本文提出一种基于非平稳系统辨识(Non-Stationary System Identification,NSSI)理论的包络提取方法。

心音信号具有短时平稳性,因此将时间段划分得足够小时,相邻两个时间段内的心音信号具有较大的相似性。而对噪声和部分杂音而言,其非平稳性特征更加明显,在适当长度的时间段内,其相似性更低。根据非平稳系统辨识原理,将原始心音段作为参考通道信号,将后移一个时间段的信号作为另一个对比通道,通过非平稳系统辨识来求解其相关性,从而提取得到心音信号的特征包络,具体步骤如下:

1)对采集的心音信号做分帧处理。在此过程中,帧长wlen的选择较为关键,其关系着心音与噪声包络特征的区分程度。文献[17]指出,当心音信号为10 ms～30 ms时具有短时平稳性,因此选取的帧长应在此范围内。此外,帧移的选择也需注意,帧移过小会加大计算量,帧移过大会使部分对应帧间的相关度偏低,使得特征包络不明显,实际处理中通常取其为帧长的一半[2,10]。本文按式(6)测试了多种类型心音信号下的相关度,并给出部分结果,如图2所示。由于噪声及心音信号都具备非平稳性,且噪声具有随机性,心音具有一定的随机变异性,因此图2中相关度值波动较大属于正常现象。在选择帧长时,要选择s1和s2相关度值差距较小且能明显区分于噪声的帧长。图2中用矩形框标记了各情况下的可取帧长,为方便起见,选取普适性更强的16 ms帧长。

图2 不同病况与帧长下的相关度测试结果Fig.2 Test results of the correlation at different patient conditions and frame lengths

2)求取参考信号x1与对比信号x2。将分帧处理后的信号去掉最后一帧,得到参考信号x1;将分帧处理后的信号去掉第一帧,得到对比信号x2。

(7)

(8)

(9)

(10)

其中,xcorr表示求相关,FFT表示快速傅里叶变换。

4)按照式(11)给出的非平稳系统相关性辨识结果求得包络特征Hx。

(11)

5)每帧的信号Hx值个数与帧长相同,因此需要选择一个代表性的值来表示该帧的包络特征值,本文按照常规方法取升序排列50%的值,即取中值。由筛选的中值组成一组特征值构成提取包络。特征包络Hx能够较好地提取基础心音信号特征,区分出基础心音、噪声和杂音。由于特征包络的整体呈脉冲状,若不做平滑处理,则会出现错误分割现象,因此需要对包络作进一步处理。本文首先采用最值滤波来平滑特征包络Hx,即取窗内最大值为平滑值,平滑窗长取为4(根据帧长确定);然后采用锯齿展宽处理,得到Hx_s特征包络。

如图3所示,本文将降噪后的心音信号图与所提取的包络图进行对比,并用矩形框分别标记时域干扰音及其对应的特征包络。由图3可以看出,提取的包络能较好地体现基础心音特征并湮没噪声,而香农能量包络法和希尔伯特变换包络法却无法湮没难以抑制的噪声,在多组心音信号测试中均存在此类现象。由于良好的特征包络将有助于分割出基础心音信号,因此包络的有效性与普适性将通过分割结果进一步验证。

图3 特征包络提取结果Fig.3 Extraction results of feature envelopes

3.2 自适应阈值及分割点判定

在提取出包络的基础上进行阈值选取分析。为避免人工设置阈值的不便,本文采用提取相关阈值参数来自适应分配阈值。一般地,阈值是有用信号与噪声的分界线,因此所以首先选择能够反映信号中噪声程度的参数,其次阈值被用于检测特征包络是否满足要求,所以阈值还应与提取的特征包络相关。此外,由于类似基础心音信号的杂音干扰或受试者心音异常,所以分割结果还需进一步修正,以保证分割结果的正确性。整个过程具体如下:

1)选择反映噪声程度的参数——信噪比(Signal to Noise Ratio,SNR),但由于实际中的噪声或者纯净信号,而通常是不可预知的,因此本文采用降噪的方式来估计纯净信号,获取SNR。在包络提取时已经进行小波降噪,但并不能将噪声完全滤除,因此进一步加大滤波参数为cora80%来抑制噪声估计纯净心音,按式(12)计算SNR[12]。SNR估计值将用于提取阈值T1。

(12)

其中,s(n)为估计的纯净信号,d(n)为噪声,N为信号长度。

(13)

其中,NH为包络长度。

3)按照选择的参数建立阈值函数。在包络分割中,常用的阈值函数形式为双门限阈值,建立双阈值函数如式(14)和式(15)所示:

(14)

(15)

若SNR为0,则:

(16)

其中,α和β为常数,文献[18]指出阈值系数可通过遗传算法进行全局优化,实验得到α=1.4、β=1.6为最佳值。

4)根据阈值寻找分割点。分割点分为上分割点与下分割点[19],即沿时间轴取包络与阈值T1的交点ui和包络与阈值T2的交点rl的中点为分割点,当ui>rl时为下分割点,当ui

5)计算每组分割点间包络所对应的面积,对分割点进行第一次筛检[5,9]。由于所提取的基础心音信号分割点间包络所围的面积要远大于噪声信号分割点间的面积,因此可舍去所围面积过小的分割点,排除噪声包络的干扰。

6)通过时域检查波峰对分割点进行第二次筛检。提取每组分割点间的最高波峰,沿时间轴可构成一组波峰时间点{pi|p1,)p2,…,pm},若波峰间的时间差满足Δpi=pi+1-pi<0.1 s,则认为pi+1或pi所对应的分割点为错误分割点。如果采集的受试者心音存在类似基础心音的心杂音、第三心音、第四心音,或者是奔马律等异常类型心音,又或者是首位存在不完整的基础心音信号,都极有可能导致分割点出错。而由于目前的包络分割方法多数未考虑此类情况,因此在实际应用中容易加大分割误差。本文根据心音信号的时域特征,提出相应的误分割点检验方法,具体步骤如下:

(1)识别不完整的基础心音,去除单个分割点。由于用于分割的心音信号开始处和结束处可能存在不完整的基础心音,会造成分割点判定结果存在单个的上分割点或下分割点,单个分割点会影响υi和νl匹配,同时会影响起始基础心音信号的识别,因此需要检测出单个的υi和νl,并将其从检测结果中去除。在不完整的基础心音检测中,将第一个上分割点与第一个下分割点进行比较,若第一个上分割点位置在第一个下分割点位置的后面,则判定ν1为不完整的下分割点,并利用类似方法检测末尾的分割点,整个非完整基础心音检测及滤除过程可表示为:

(17)

其中,Qdel表示需要删除的分割点,υend表示最后一个上分割点,νend表示最后一个下分割点。

(2)检测心杂音或异常额外心音的分割点。实际中采集到的很多心杂音或异常额外心音与基础心音非常类似,这些非基础心音成分会对分割点的判定造成干扰。文献[20]指出收缩期时长一般为160 ms～240 ms,在检测误检分割点时,取收缩期时长最小值的一半(80 ms),由于基础心音信号包络值较非基础心音信号包络值更大,因此可采用如下判别方法:

(18)

当误检点出现在舒张期的中间位置时,上述检测方法难以有效判定,因此需要进一步检测。文献[20]指出心动周期时长一般为600 ms～1 000 ms,取心动周期最短时长作为又一判定参量,采用如下检测方法:

(19)

7)提取时域特征,识别心音成分。舒张期持续时长一般大于收缩期时长[18-19],故根据时域时长判定持续时间较长的成分为舒张期,较短的成分为收缩期[20-21],按顺序处于舒张期与收缩期之间的成分为s1,则剩余成分为s2。

3.3 分割流程及评价指标

根据上文所述方法完成包络提取及分割任务,具体流程如图4所示。

图4 自适应阈值分割流程Fig.4 Procedure of adaptive threshold segmentation

由于本文方法主要面向复杂心音信号实时分割处理,难以通过直观的分割定位结果对分割算法进行有效性检验,因此本文采用文献[22]提出的4类评价指标作为检验指标:

(20)

(21)

(22)

PS=run_time

(23)

其中,TPFHS为基础心音信号检出正确率,FPFHS为非基础心音信号检入错误率,MPFHS为基础心音未检出率,TNFHS为正确检入的基础心音信号帧数,FNFHS为错误检入的非基础心音信号帧数,MNFHS为未检入的基础心音信号帧数,PS为算法运行时间。TPFHS越高,FPFHS和TNFHS越低,代表分割精度越高;PS越小,代表实时性越好。

4 实验结果与分析

4.1 实验数据及其处理

实验数据由中国医学科学院阜外医院和南京邮电大学电子科学与工程学院提供。其中选用前面提供心音62条,后者提供心音43条,所有心音信号均降采样至2 000 Hz。为模拟复杂环境,选用noiseX-92噪声数据库中的15类不同领域的噪声按5 dB、10 dB、15 dB这3种信噪比加入心音信号中构成待处理复杂心音,噪声叠加不超过两种,选用的心音信号中还包括17类异常心音信号。noiseX-92噪声数据库中的噪声类型具体为白噪声、粉红噪声、高频通道噪声、飞机驾驶舱噪声1、飞机驾驶舱噪声2、厂房噪声1、厂房噪声2、类语音噪声、坦克噪声、驱逐舰操作室噪音、驱逐舰机舱噪音、F16机舱噪声、军用车辆噪声、机枪噪音和车内噪声。心音信号类型及数量如表1所示。

表1 心音信号类型及数量Table 1 Type and quantity of heart sound signal

105条心音信号经处理后生成33 075条带噪心音信号,共计516 285个心动周期,原始心音信号均在Cool Edit Pro下进行手工标注分割点。实验软件为Matlab2017a,计算机CPU为Intel®CoreTMi7-8700k@3.70 GHz,RAM为16.0 GB,硬盘为1 TB。

4.2 结果分析

本文方法分割结果如图5所示,直观反映了其能够较精准地分割出基础心音信号,且不会将非基础心音误分割为基础心音信号。本文方法与维奥拉积分包络分割法[5]、改进型希尔伯特-黄变换包络双阈值分割法[7]、基于心动周期估计的心音包络分割法[9]、基于周期与波峰的自适应阈值包络分割法[10]和基于个性化高斯混合建模的包络分割法[4]进行对比,其中基于个性化高斯混合建模的包络分割法采用在线训练方式,无需通过不同的训练集进行训练,能够避免对比差异。

图5 基于本文方法的心音信号分割结果Fig.5 Segmentation results of heart sound signal based on the proposed method

表2给出了6种分割方法的评价指标结果,可以看出,基于个性化高斯混合建模的分割法检出正确率TPFHS最高为91.49%,其次为本文方法的89.21%,其余方法的检出准确率相对较低,这是因为本文方法可避免杂音和误检干扰。此外,可以发现其他4种基于包络特征的分割方法中基础心音错误检入率FPFHS要高于未检入率MPFHS,主要原因为在基础心音末尾处噪声存在拖尾,容易使非基础心音的噪声部分被误判为基础心音帧,而本文方法不存在此类问题。

表2 6种方法的分割性能比较Table 2 Comparison of segmentation performance of six methods

另外,在算法耗时PS指标评测中,采用的心音信号段为4.5 s,即每5个心动周期为一段,获取PS的平均值为最终结果。从表2可以看出,基于个性化高斯混合建模的包络分割法的耗时为6.573 5 s,远高于其他分割方法,这是因为基于个性化高斯混合建模的包络分割法采用在线训练方式,涉及到高斯混合建模、前向后向运算、Viterbi解码、误差修正等计算[4],复杂度远高于其他基于包络的双阈值分割方法;基于改进型希尔伯特-黄变换包络的分割法先进行经验模态分解,然后做3次贝塞尔样条插值运算,同时对端点效应进行优化[7],其运算复杂度次于基于个性化高斯混合建模的包络分割法;维奥拉积分包络分割法在确定尺度与基础心音的关系后直接提取包络进行分割[5],基于周期与波峰的自适应阈值包络分割法建立自相关估计周期与波峰的关系阈值后完成分割[10],本文方法在提取心音包络后直接进行分割与判别,复杂度均远低于前两种方法。本文方法的耗时为0.105 2 s,能够满足实时分割要求,可为实时监测的心音自动分析技术提供参考。

5 结束语

基于包络的心音信号分割方法因算法简单、时效性强等特点而得到广泛应用,其中的特征包络提取与阈值处理算法设计为关键步骤。本文提出一种基于非平稳系统辨识的心音特征包络提取方法,采用遗传算法设计自适应双阈值函数实现阈值分割。在近似基础心音的心杂音处理上,根据时域特征进行分割点筛查,识别并去除误检分割点和单独分割点,从而得到正确的基础心音分割点。实验结果表明,本文方法相比维奥拉积分包络分割法、改进型希尔伯特-黄变换包络双阈值分割法、基于心动周期估计的心音包络分割法和基于周期与波峰的自适应阈值分割法分割精度更高,相比基于个性化高斯混合建模的包络分割法实时性更强,提取的包络能够更好地反映基础心音信号特征。但本文方法对于人体运动状态下的心音信号分割精度不理想,后续将对此做进一步研究。