CN104473660A - 一种基于子带能量包络自相关特征的异常心音识别方法 - Google Patents

Abstract

基于子带能量包络自相关特征的异常心音识别方法，涉及一种异常心音识别方法，解决心音识别中存在依赖心音分段、不能处理变长信号、无法适应现实噪声环境的特征提取和分类识别的问题。将输入的心音信号进行能量标准化处理后降采样并进行带通滤波；对近似子带能量包络信号和细节子带能量包络信号分别计算其自相关序列，并分别截取两个序列的前M个值作为输入心音信号的近似子带能量包络自相关特征和细节子带能量包络自相关特征；依据能量包络自相关特征分别构造作用其上的两个由R^M空间到R^e空间的散射映射，降维后进行融合得到心音信号的能量包络自相关特征；测试数据进行特征提取并输入到分类模型中进行分类。避免了心音分段过程，而且提高了噪声环境下的鲁棒性。

Description

一种基于子带能量包络自相关特征的异常心音识别方法

技术领域

本发明涉及一种异常心音识别方法，涉及生物信号识别技术及智能信息处理领域。

背景技术

据世界卫生组织统计的结果，心血管疾病已经成为当今威胁人类生命健康的头号杀手。2004年大约有1700万人死于心血管疾病，占全部疾病的29％，其中720万人死于心脏病，而且由心脏导致的患病和死亡的人数和比例还在逐年增加。在我国，依据“中国心血管病报告2012”公布的结果：我国心血管病患病率处于持续上升态势，估计心血管病的患者人数为2.9亿，即每10个成年人中就有2人患心血管病，每年约有350万人死于心血管病，占总死亡原因的41％，居各种疾病之首。因此，及时发现心血管疾病的患病症状、掌握病情状态对保证人类的生命健康具有极其重要的意义。

人体是一个异常复杂的有机整体，当人体患病时会产生各种生理信号，通过解读这些生理信号能够及早地发现所患疾病，并能够对病情实行有效地监控。对心血管疾病而言，心音(Heart Sound)就是人体最为重要的生理信号之一。心音是心脏及心血管系统机械运动状况的反映，是在心动周期中由于心肌收缩舒张、瓣膜启闭以及血流冲击心室壁和大动脉等引起的一种机械振动。在医学研究中，心音图(Phonocardiogram，PCG)是评估心脏功能状态的一种基本的方法。与心电图(Electrocardiogram，ECG)方法相比，心脏或心血管疾病所表现出的病理现象会较早地体现在心音信号上，在患者尚未感觉到痛感和ECG异常等症状之前，心音就可以为早期诊断提供可靠的信息。基于心音信号的传统医疗手段是心音听诊和心音图，然而它们在临床应用中皆有不足之处。听诊需要医生使用听诊器依靠耳听的方式进行诊断，然而心音的主要频率成分是在40～150Hz之间，而人耳听觉系统自身频率敏感的范围是1000～3000Hz，这意味着这种通过人耳听诊方式会损失大量的信息，从而影响对疾病的诊断。心音图通过将心音信号进行图像化弥补了听诊的不足，但需要更为专业的设备，也增加了诊疗的费用。除设备和听诊方式自身的不足外，更为重要的是，以上两种方法均需要依赖医生自身的专业知识和经验才能做出诊断，其诊断过程是主观的，缺少定量的分析方法，因而难以保证诊断结果的可靠性，往往还需要借助更为复杂、费用更高的其他诊疗方法进行确诊。此外，对普通人或需要经常进行病情监控的患者来说，由于其既没有专业方面的知识，也无法购置昂贵的医疗设备，上述两种方法显然都难以实现家庭听诊。

心音是声音信号的一种，其蕴含着有关人体健康状况的重要信息，通过提取这种信息并进行有效的识别能够实现客观的数字化心音听诊，从而能为患者提供可靠的诊断结果。就心音识别而言，心音分段与心音分类是其主要研究内容。心音分段主要是对心音信号中出现的第一心音信号S1和第二心音信号S2的检测，并进而确定心动周期等重要生理参数，同时也是接下来心音分类任务的预处理。心音分类需进一步从异常心音信号中识别出是何种心脏杂音(如心缩期杂音、心舒期杂音)以及其他干扰声音等，并进而分辨出相应的疾病类型。在现实的应用环境下，采集到的心音信号通常会伴有多种干扰噪声，如听诊器与衣服或皮肤的摩擦声、背景音乐声、人的说话声等，因此这需要心音识别方法具有较强的鲁棒性。为此，欧盟支持的PASCAL(Pattern Analysis,Statistical Modelling andComputational Learning)计划将噪声环境下的心音分段和分类作为2012年度的挑战任务。然而，噪声环境下、结构复杂、非平稳的心音信号识别是一个新兴的、较难解决的问题，目前尚处于研究的起步阶段，缺少系统性的研究方法。

1)现有的心音识别/分类方法大多需要先对输入心音信号进行分段处理，依据分段的结果进行特征提取和分类识别，因而其识别性能在很大程度上依赖于心音分段的精度。然而，严重心脏病患者的心音中通常伴有显著的心杂音(murmur)，而且在现实中的心音采集(听诊)过程中通常伴有噪声(如听诊器与皮肤或衣服的摩擦音、语音、音乐及其他环境噪声等)，这些情况都会大大降低心音分段的精度，从而明显降低心音识别系统的性能。

2)现有的心音识别/分类方法通常以固定时长来处理心音信号，以便进行特征提取和分类识别。然而在现实应用中所采集的心音信号通常具有较大的变动范围，这限定了固定长度识别方法的使用。

发明内容

要解决的技术问题：

为解决上述心音识别中存在依赖心音分段、不能够处理变长输入信号、无法适应于现实噪声环境的特征提取和分类识别的问题，进而提供一种基于子带能量包络自相关特征的异常心音识别方法。

本发明为解决上述技术问题采取的技术方案是：

一种基于子带能量包络自相关特征的异常心音识别方法，由以下步骤实现：

步骤一、将输入的心音信号进行能量标准化处理，然后降采样为2kHz，并对其采用6阶巴特沃斯滤带通(20～900Hz)波器对降采样后的心音信号进行带通滤波以滤除截止频率以外的声音和噪声，得到预处理后的信号；

步骤二、依据心音和心杂音信号的频率分布特点，采用4阶Daubechies小波对步骤一中预处理后的信号进行4层分解，对分解所得到的小波近似子带系数a₄和细节子带系数d₂分别进行分帧处理并提取均值香农能量包络，从而得到近似子带能量包络信号E_a和细节子带能量包络信号E_d；

步骤三、对步骤二中所得到的近似子带能量包络信号E_a和细节子带能量包络信号E_d分别计算其自相关序列，并分别截取两个序列的前M个值作为输入心音信号的近似子带能量包络自相关特征r_a∈R^M和细节子带能量包络自相关特征r_d∈R^M,其中M的取值范围为:200≤M≤300；

步骤四、依据步骤三中所得到的能量包络自相关特征r_a和r_d分别构造作用其上的两个由R^M空间到R^e空间的散射映射(diffusion maps)ψ_a:R^M→R^e和ψ_d:R^M→R^e，两种能量包络自相关特征r_a和r_d经各自的散射映射变换后的结果为：ψ_a(r_a)∈R^e和ψ_d(r_d)∈R^e，其中e<M，从而实现降维；然后，再将ψ_a(r_a)和ψ_d(r_d)以向量拼接的方式融合在一起，最终得到心音信号的能量包络自相关特征z；

步骤五、依据步骤一至四所述的特征提取过程，对训练集数据中L个心音信号进行特征提取得到L个心音能量包络自相关特征并结合每个心音信号的类别标识信息y_i得到训练集然后，以支撑向量机(SVM)为分类器训练分类模型，从而完成模型训练过程；

步骤六、依据步骤一至四所述的特征提取过程，对测试数据中J个心音信号进行特征提取得到J个心音能量包络自相关特征将输入到由步骤五中训练得到分类模型中进行分类。

步骤三中，依据下式计算近似子带能量包络自相关特征r_a和细节子带能量包络自相关特征r_d

$r_a\left(m\right)=\underset{n=0}{\overset{N-m-1}{\mathrm{\&Sigma;}}}{E}_a(n+m)E_a\left(n\right),m\&GreaterEqual;1$

式中：

当x取“a”时，E_x(n)表示近似子带能量包络信号，当x取“d”时，E_x(n)表示细节子带能量包络信号，n为能量包络信号的时间索引，m为时间索引的偏移量。

步骤四中散射影射ψ_a:R^M→R^e和ψ_d:R^M→R^e的构造过程为：在构造散射影射ψ_a和ψ_d时，除所采用的特征集合不同外，其构造过程完全相同，因此统一表示为散射影射ψ，并以表示由L个心音信号中提取出的能量包络自相关特征，r_i为近似子带能量包络自相关特征r_a或细节子带能量包络自相关特征r_d；散射影射ψ的构造过程包含以下四个步骤：

步骤四(一)、构造特征集上的转移矩阵P；

步骤四(二)、对转移矩阵P进行特征值分解以得到散射影射ψ；

步骤四(三)、对散射影射ψ进行扩散以使其可以应用于新数据。

步骤四(一)中使用如下方法构造特征集上的转移矩阵P：采用高斯核函数度量特征集中任意两个特征r_i和r_j间的相似关系w(r_i,r_j)：

$w(r_i,r_j)=\exp \left\{\frac{{\left|\right|r_i-r_j\left|\right|}^2}{{\mathrm{\&delta;}}^2}\right\}$

其中：δ为核宽；则特征集上的转移矩阵P定义为

P＝D^-1W

其中：矩阵W中的元素定义为W_ij＝w(r_i,r_j)，对角阵D的对角元素定义为

步骤四(二)中使用如下方法对转移矩阵P进行分解并得到散射影射ψ：

对转移矩阵P进行特征值分解：

P＝USV^T

其中：U＝[u₁,...,u_K]和V＝[v₁,...,v_K]分别为转移矩阵P的右特征向量和左特征向量，S是由P的特征值所构成的对角阵列，即S＝diag(λ₁,...,λ_K)；则对特征集中的任意特征r∈Ω，其散射映射定义为：ψ(r)＝[λ₁u₁(r),...,λ_Ku_K(r)]。

步骤四(三)中使用如下方法对散射影射ψ进行扩散以使其可以应用于新数据：对任意特征其散射映射由下式得到

ψ(r^new)＝p^newU

其中：U为对转移矩阵P进行特征分解时所得到的右特征向量， $w_i^{\mathrm{new}}=\exp \{-\frac{{\left|\right|r^{\mathrm{new}}-r_i\left|\right|}^2}{{\mathrm{\&delta;}}^i}\},r_i\&Element;\mathrm{\&Omega;},d^{\mathrm{new}}=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{w}_i^{\mathrm{new}}.$

步骤四中不同能量包络自相关特征r_a和r_d的散射映射ψ_a(r_a)和ψ_d(r_d)按下式进行融合得到能量包络自相关特征z：

$z=\left[\begin{array}{c}{\mathrm{\&psi;}}_a\left(r_a\right)\\ {\mathrm{\&psi;}}_d\left(r_d\right)\end{array}\right]\&Element;R^{2e}.$

本发明具有以及有益效果：本发明通过使用不同小波子带能量包络的自相关特征来表示整体心音信号信息，不仅避免了心音分段过程，而且提高了噪声环境下的鲁棒性。通过使用散射映射将两种自相关特征进行融合，提高了系统的识别性能。本发明方法不依赖心音分段、能够处理变长输入信号、适应于现实噪声环境的特征提取和分类识别方法，是一种无需进行心音分段处理、不受固定时长限制、可应用于现实噪声环境下的异常心音识别方法。

附图说明

图1为本发明所提出的心音识别方法的流程框图。

具体实施方式

具体实施方式一：本实施方式所述的基于子带能量包络自相关特征的异常心音识别方法由以下步骤实现：

步骤一、心音预处理

将输入的心音信号进行能量标准化处理，然后降采样为2kHz，并对其采用6阶巴特沃斯滤带通(20～900Hz)波器对降采样后的心音信号进行带通滤波以滤除截止频率以外的其他声音和噪声；注：假定输入的心音信号采样频率高于2kHz，预处理后的输出信号采样频率近似为2kHz；

步骤二、提取平均香农能量包络

(1)对预处理后采样频率为2kHz的心音信号采用4阶Daubechies小波进行4层分解得到：一个近似系数序列a₄和四个细节系数序列d₄,d₃,d₂,d₁，它们分别对应着信号中的频率范围分别为0～125Hz，125～250Hz，250～500Hz，500～1000Hz和1000～2000Hz的信号成分；

(2)依据心音和心杂音信号的频率分布特点，分别选取小波子带近似系数a₄和小波子带细节系数d₂计算其平均香农能量包络：

$E_a=-\frac{1}{M_a}\underset{i=1}{\overset{M_a}{\mathrm{\&Sigma;}}}{a}_4{\left(i\right)}^2{\log a}_4{\left(i\right)}^2---\left(1\right)$

$E_d=-\frac{1}{M_d}\underset{i=1}{\overset{M_d}{\mathrm{\&Sigma;}}}{d}_4{\left(i\right)}^2{\log d}_4{\left(i\right)}^2---\left(2\right)$

其中：M_a和M_d分别为作用于两个子带系数的窗口长度，E_a和E_d是所对应的子带平均香农能量包络；

步骤三、计算自相关特征

对步骤二中所得到的能量包络信号E_a和E_d依下式分别计算其单边自相关函数：

$r_x\left(m\right)=\underset{n=0}{\overset{N-m-1}{\mathrm{\&Sigma;}}}{E}_x(n+m)E_x\left(n\right),m\&GreaterEqual;1---\left(3\right)$

其中：当x取“a”时，E_x(n)表示近似子带能量包络信号，当x取“d”时，E_x(n)表示细节子带能量包络信号，n为能量包络信号的时间索引，m为时间索引的偏移量。

步骤四、子带特征融合

为了将所得到的子带特征使用适当的分类器进行分类，需要对其进行融合处理，为此需要为每个子带特征分别构造散射映射，然后依据映射后的结果进行特征融合。子带特征融合过程由以下步骤构成：

(1)构造特征集上的转移矩阵P。令为由L个训练集心音得到的子带特征集合，采用高斯核函数度量特征集中任意两个特征r_i和r_j间的相似关系w(r_i,r_j)：

$w(r_i,r_j)=\exp \left\{\frac{{\left|\right|r_i-r_j\left|\right|}^2}{{\mathrm{\&delta;}}^2}\right\}---\left(4\right)$

其中：δ为通过核宽。则特征集上的转移矩阵P定义为

P＝D^-1W   (5)

其中：矩阵W中的元素定义为W_ij＝w(r_i,r_j)，对角阵D的对角元素定义为

(2)计算散射映射：对转移矩阵P进行特征值分解：

P＝USV^T   (6)

其中：U＝[u₁,...,u_K]和V＝[v₁,...,v_K]分别为转移矩阵P的左特征向量和右特征向量，S是由P的特征值所构成的对角阵列，即S＝diag(λ₁,...,λ_K)。通过将特征向量的参数化可得到任意x∈Ω的散射映射为：

ψ(x)＝[λ₁u₁(x),...,λ_Ku_K(x)]   (7)

(3)扩展散射映射：对任意特征其散射映射由下式得到

ψ(r^new)＝p^newU   (8)

其中：U为对转移矩阵P进行特征分解时所得到的右特征向量， $w_i^{\mathrm{new}}=\exp \{-\frac{{\left|\right|r^{\mathrm{new}}-r_i\left|\right|}^2}{{\mathrm{\&delta;}}^i}\},r_i\&Element;\mathrm{\&Omega;},d^{\mathrm{new}}=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{w}_i^{\mathrm{new}};$

(4)子带特征融合：依据上述步骤(1)(2)(3)对训练集数据所得到的自相关特征r_a和r_d分别构造两种特征集上的散射映射ψ_a和ψ_d，然后将映射后的特征合并为最终特征：

$z=\left[\begin{array}{c}{\mathrm{\&psi;}}_a\left(r_a\right)\\ {\mathrm{\&psi;}}_d\left(r_d\right)\end{array}\right]\&Element;R^{2e}---\left(9\right)$

步骤五、依据训练数据的特征集合训练支撑向量机(SVM)模型，并使用此模型对新输入的心音信号进行分类；

步骤六、依据步骤一至四所述的特征提取过程，对测试数据中J个心音信号进行特征提取得到J个心音能量包络自相关特征将输入到由步骤五中训练得到分类模型中进行分类。

实验结果：

为验证本发明所提出方法的有效性并对其性能进行评估，实验采用欧盟支持的PASCAL挑战的公开数据库作为测试数据。测试数据包含A和B两个数据库，其中数据库A采用数字听诊器采集，包含正常心音、心杂音、额外心音和其他声音(语音、音乐等)四个类别。数据库B采用苹果手机采集，包含正常心音、心杂音和早搏心音三个类别。实验以挑战赛所公布的伦敦大学学院和波尔图大学的实验结果为参考基线(其中波尔图大学的识别算法取得了本次比赛的第一名)，以挑战赛所提供的评估程序做性能评估方法。对数据库A和数据库B的识别结果如表1和表2所示。实验结果表明：本发明所提出识别方法在两个数据库上几乎全部超过了参加比赛算法的识别结果，这充分证明了本发明方法的有效性。

表1 数据库A的识别结果

表2 数据库B的识别结果

Claims (8)

1.一种基于子带能量包络自相关特征的异常心音识别方法，其特征在于它由以下步骤实现：

步骤一、将输入的心音信号进行能量标准化处理，然后降采样为2kHz，并对其采用6阶巴特沃斯滤带通波器对降采样后的心音信号进行带通滤波以滤除截止频率以外的声音和噪声，得到预处理后的信号；

步骤二、依据心音和心杂音信号的频率分布特点，采用4阶Daubechies小波对步骤一中预处理后的信号进行4层分解，对分解所得到的小波近似子带系数a₄和细节子带系数d₂分别进行分帧处理并提取均值香农能量包络，从而得到近似子带能量包络信号E_a和细节子带能量包络信号E_d；

步骤三、对步骤二中所得到的近似子带能量包络信号E_a和细节子带能量包络信号E_d分别计算其自相关序列，并分别截取两个序列的前M个值作为输入心音信号的近似子带能量包络自相关特征r_a∈R^M和细节子带能量包络自相关特征r_d∈R^M,其中M的取值范围为:200≤M≤300；

步骤四、依据步骤三中所得到的能量包络自相关特征r_a和r_d分别构造作用其上的两个由R^M空间到R^e空间的散射映射ψ_a:R^M→R^e和ψ_d:R^M→R^e，两种能量包络自相关特征r_a和r_d经各自的散射映射变换后的结果为：ψ_a(r_a)∈R^e和ψ_d(r_d)∈R^e，其中e<M，从而实现降维；然后，再将ψ_a(r_a)和ψ_d(r_d)以向量拼接的方式融合在一起，最终得到心音信号的能量包络自相关特征z；

步骤五、依据步骤一至四所述的特征提取过程，对训练集数据中L个心音信号进行特征提取得到L个心音能量包络自相关特征并结合每个心音信号的类别标识信息y_i得到训练集然后，以支撑向量机为分类器训练分类模型，从而完成模型训练过程；

步骤六、依据步骤一至四所述的特征提取过程，对测试数据中J个心音信号进行特征提取得到J个心音能量包络自相关特征将输入到由步骤五中训练得到分类模型中进行分类。

2.根据权利要求1所述的一种基于子带能量包络自相关特征的异常心音识别方法，其特征在于步骤三中，依据下式计算近似子带能量包络自相关特征r_a和细节子带能量包络自相关特征r_d

$r_a\left(m\right)=\underset{n=0}{\overset{N-m-1}{\mathrm{\&Sigma;}}}{E}_a(n+m)E_a\left(n\right),m\&GreaterEqual;1$

式中：

当x取“a”时，E_x(n)表示近似子带能量包络信号，当x取“d”时，E_x(n)表示细节子带能量包络信号，n为能量包络信号的时间索引，m为时间索引的偏移量。

3.根据权利要求1或2所述的一种基于子带能量包络自相关特征的异常心音识别方法，其特征在于步骤四中散射影射ψ_a:R^M→R^e和ψ_d:R^M→R^e的构造过程为：在构造散射影射ψ_a和ψ_d时，除所采用的特征集合不同外，其构造过程完全相同，因此统一表示为散射影射ψ，并以表示由L个心音信号中提取出的能量包络自相关特征，r_i为近似子带能量包络自相关特征r_a或细节子带能量包络自相关特征r_d；散射影射ψ的构造过程包含以下四个步骤：

步骤四(一)、构造特征集上的转移矩阵P；

步骤四(二)、对转移矩阵P进行特征值分解以得到散射影射ψ；

步骤四(三)、对散射影射ψ进行扩散以使其可以应用于新数据。

4.根据权利要求3所述的一种基于子带能量包络自相关特征的异常心音识别方法，其特征在于步骤四(一)中使用如下方法构造特征集上的转移矩阵P：采用高斯核函数度量特征集中任意两个特征r_i和r_j间的相似关系w(r_i,r_j)：

$w(r_i,r_j)=\exp \{-\frac{{\left|\right|r_i-r_j\left|\right|}^2}{{\mathrm{\&delta;}}^2}\}$

其中：δ为核宽；则特征集上的转移矩阵P定义为

P＝D^-1W

其中：矩阵W中的元素定义为W_ij＝w(r_i,r_j)，对角阵D的对角元素定义为

5.根据权利要求4所述的一种基于子带能量包络自相关特征的异常心音识别方法，其特征在于步骤四(二)中使用如下方法对转移矩阵P进行分解并得到散射影射ψ：

对转移矩阵P进行特征值分解：

P＝USV^T

其中：U＝[u₁,...,u_K]和V＝[v₁,...,v_K]分别为转移矩阵P的右特征向量和左特征向量，S是由P的特征值所构成的对角阵列，即S＝diag(λ₁,...,λ_K)；则对特征集中的任意特征r∈Ω，其散射映射定义为：ψ(r)＝[λ₁u₁(r),...,λ_Ku_K(r)]。

6.根据权利要求5所述的一种基于子带能量包络自相关特征的异常心音识别方法，其特征在于步骤四(三)中使用如下方法对散射影射ψ进行扩散以使其可以应用于新数据：对任意特征其散射映射由下式得到

ψ(r^new)＝p^newU

其中：U为对转移矩阵P进行特征分解时所得到的右特征向量， $w_i^{\mathrm{new}}=\exp \{-\frac{{\left|\right|r^{\mathrm{new}}-r_i\left|\right|}^2}{{\mathrm{\&delta;}}^2}\},r_i\&Element;\mathrm{\&Omega;},d^{\mathrm{new}}=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{w}_i^{\mathrm{new}}.$

7.根据权利要求1、2、4、5、或6所述的一种基于子带能量包络自相关特征的异常心音识别方法，其特征在于步骤四不同能量包络自相关特征r_a和r_d的散射映射ψ_a(r_a)和ψ_d(r_d)按下式进行融合得到能量包络自相关特征z：

$z=\left[\begin{array}{c}{\mathrm{\&psi;}}_a\left(r_a\right)\\ {\mathrm{\&psi;}}_d\left(r_d\right)\end{array}\right]\&Element;R^{2e}.$

8.根据权利要求3所述的一种基于子带能量包络自相关特征的异常心音识别方法，其特征在于步骤四不同能量包络自相关特征r_a和r_d的散射映射ψ_a(r_a)和ψ_d(r_d)按下式进行融合得到能量包络自相关特征z：

$z=\left[\begin{array}{c}{\mathrm{\&psi;}}_a\left(r_a\right)\\ {\mathrm{\&psi;}}_d\left(r_d\right)\end{array}\right]\&Element;R^{2e}.$