针对先天性心脏病相关肺动脉高压听诊特征不明显，已有的机器辅助诊断算法相对复杂等问题，提出一种基于第二心音信号高频分量统计特征的分析方法。首先，采用端点检测自适应分割方法提取第二心音。其次，使用离散小波变换分解出高频分量，并提取该分量的赫斯特（Hurst）指数、勒佩尔-齐夫（Lempel-Ziv）信息和样本熵等统计特征。最后，使用这些特征训练极端梯度提升算法（XGBoost）分类器，在三分类中准确率达到了80.45%。该方法无需进行降噪处理，特征提取速度快，且只需三个特征即可实现较好的多分类效果，有望用于先天性心脏病相关肺动脉高压早期筛查。
引用本文： 杨炫锴, 孙静, 杨宏波, 郭涛, 潘家华, 王威廉. 基于第二心音统计特征的先天性心脏病相关肺动脉高压诊断方法. 生物医学工程学杂志, 2024, 41(1): 41-50. doi: 10.7507/1001-5515.202304037 复制
0 引言
肺动脉高压（pulmonary arterial hypertension，PAH）是一种潜在致命性心血管疾病，其特征为肺动脉压力升高，当压力超过一定阈值时，会使得右心房扩大，最终导致右心衰竭[1]。然而，PAH的症状通常在疾病晚期才显现，早期筛查很容易遗漏。PAH的病因复杂，其中约40%是由先天性心脏病（congenital heart disease，CHD）引起的[2]。CHD是一种胎儿时期发育异常导致的心血管畸形，心脏缺陷导致肺循环的压力与容量超负荷，最终形成PAH。CHD相关PAH（CHD-associated PAH，CHD-PAH）由于同时夹杂着两种症状，致死率更高。有创介入式的右心导管术是诊断PAH的金标准，但存在较大的风险。目前无创诊断最有价值的评估方法是超声心动图检查，但在偏远地区，由于缺乏设备和技术人员支持，且该设备不便携带，难以在大规模筛查中推广。
心音是一种机械振动，是心脏血流动力学和心血管系统相互作用的结果。由于心音包含了大量有关心血管系统生理状况的信息，因此可以用于评估心脏的健康状况[3]。计算机智能辅助诊断（computer-aided diagnosis，CAD）通过提取心音信号的病理特征，进行机器学习训练并构建模型用于疾病辅助诊断，能有效配合超声心动图，提高筛查效率，主要步骤分为信号预处理、特征提取、构建分类器等。
预处理阶段主要为降噪与分割。降噪能抑制采集中的呼吸音、杂音等。Dwivedi等[4] 结合稳态小波变换与集合经验模态分解（ensemble empirical mode decomposition，EEMD），来抑制环境音和噪声，以提高后续分类准确率（accuracy，Acc）。Simonneau等[5]研究表明，第一心音（first heart sound，S1）和第二心音（second heart sound，S2）与PAH的严重程度相关。其中PAH患者的S2病理特征更明显，会出现异常亢进与分裂，有助于提高PAH的诊断Acc [6]。但从心音中人工提取S2，困难且复杂[7]，而现有的心音分割算法则没有提取S2，转而通过分割出单个心动周期来研究病理特征，如基于多尺度的希尔伯特包络分割[8]和持续时间隐马尔可夫分割等[9]。因此，如果能高效准确地分割出S2，将是目前有助于促进PAH诊断的一个研究方向。
特征提取的质量会影响分类的结果，Kaddoura等[10]研究心音信号的声学特征，估计心音信号的频谱对数并进行傅里叶逆变换，以创建倒谱特征对PAH与正常心音进行二分类。Alqudah等[11]提取心音信号中S2的能量比作为特征，来区分PAH患者和健康受试者。Ge等[12]提取心动周期和S2的时频域特征等150个特征用于后续分类。
分类阶段，当数据集充足时神经网络能有效提升分类Acc。Wang等[7]基于长短时记忆方法建立了一种小波散射卷积网络，提取受试者的小波散射系数作为特征，来对正常心音与CHD二分类。对于CHD-PAH这样的小样本数据，经典的机器学习能够有效根据信号特征对心音进行分类，不易受小数据量的影响而产生过拟合。已有研究报道，随机森林（random forest，RF）[9]和K近邻法（k-nearest neighbor，KNN）[13]在心音二分类上有一定效果。
目前大多数研究都是针对PAH和正常心音的二分类[5.10-11]，CHD-PAH的多分类问题仍有待改进，比如在预处理阶段，降噪算法有可能滤除病理信息且耗时，达不到实际应用的要求。在心音分割中，大多数分割算法依赖于同步采集心电信号来对心音信号进行辅助分割，对信号采集要求高，只通过心音信号对S2部分进行分割是目前研究的难点。特征提取阶段，目前的研究集中于提取心音的声学和时频域特征，如梅尔倒谱系数（Mel-scale frequency cepstral coefficients，MFCC）[14]和梅尔频谱系数（Mel-frequency spectral coefficient，MFSC）[15]等。然而，声学特征是针对语音信号设计的，而心音信号是由血液流动冲击瓣膜所产生的信号，与语音信号存在一定的差异，并且在听诊中，PAH和CHD-PAH仅通过听觉是很难区分，因为它们之间的声学与时频特征相似度高，分类效果不理想。此外，对于CHD-PAH的多分类诊断依赖于大量特征[12]，特征提取耗时长，不利于大规模的早期筛查任务。在分类器阶段，由于CHD-PAH的数据量少，神经网络在少量数据上容易出现过拟合的问题，传统机器学习更适用于小样本数据分类。
鉴于上述现状，本文旨在建立一个快速有效的CHD-PAH诊断算法。在预处理阶段，不进行降噪，使用自适应的双阈值分割法得到心音S2部分，再提取S2高频分量中的统计特征，包括赫斯特（Hurst）指数、勒佩尔-齐夫（Lempel-Ziv）复杂度（Lempel-Ziv complexity，LZC）和样本熵（sample entropy，SampEn）。与传统方法不同的是，本文算法专注于研究心音的高频特性而非声学相关特性，这些高频分量被用来描述短而突然的变化特征杂音，能更有效地凸显出PAH患者S2亢进、分裂的特征。最后，采用贝叶斯优化的极端梯度提升算法（extreme gradient boosting，XGBoost）对正常、CHD和CHD-PAH进行三分类诊断。
1 方法
本文算法通过以下几个步骤对CHD-PAH进行诊断：① 使用一种自适应分割方法来截取出心音信号中S2部分。② 对S2部分进行离散小波变换（discrete wavelet transformation，DWT）分解，提取出高频分量。③ 使用Hurst指数、LZC和样本熵进行复杂性度量，提取出S2高频分量的统计特征。④ 特征通过XGboost分类器进行CHD-PAH的三分类诊断。CHD-PAH分类流程如图1所示。
图1
CHD-PAH流程图
Figure1.
CHD-PAH flow chart
1.1 预处理
心脏跳动是一个近似于周期性的过程，通常分为S1、收缩期、S2和舒张期四个部分。在听诊PAH患者时，S2中肺部成分加重是其确诊的主要特征之一[6]。如图2所示，局部方框圈出部分波形即为S2，PAH患者表现为S2部分出现亢进与异常分裂。相比于整个心动周期，仅使用S2部分即能有效地凸显出PAH的特征[13]。并且，将原始信号分割为单独的S2部分后，还能提高特征提取速度。
图2
CHD-PAH的S2部分出现异常
Figure2.
Anomalies in the S2 part of CHD-PAH
本文采用了一种阈值检测算法分割S2，如图1所示。首先对心音信号进行短时傅里叶变换（short-time Fourier transform，STFT），并使用汉明窗进行分帧，如式(1)所示：
 |
其中， t为时间，f为频率，j为虚数，π为圆周率，x(τ) 为t = τ时的心音信号，ω(τ − t)为窗函数。心音 x(τ)被分帧后，提取每一帧的短时能量Ei与频谱扩散Si，如式(2)～式(3)所示：
 |
 |
其中，xi(n)，n = 1, ···, N表示心音信号的第i帧，N为帧数长度，fi为第i帧的频率，u1为频谱中值，si为第i帧的频谱值， b1与b2为频谱边界值。所得Ei与Si，如图3所示。
图3
短时能量与频谱扩散
Figure3.
Short-time energy and spectral diffusion
计算短时能量与频谱扩散对应的直方图与直方图阈值TE与TS，阈值计算如式(4)所示：
 |
其中，M1、M2为直方图中前两个最大值的位置，W为设定值。根据阈值对序列进行判断并标记。其中，满足Ei ≥ TE的为标记点1（flag1），满足Si ≥ TS的为标记点2（flag2），这些标记被认为是心音S1或S2分量，称之为待定分量。如图4所示，根据阈值判断并标记序列，黑色波形是原始心音信号的频率质心，绿色波形是待定分量，即被标记为心音S1或心音S2的频率质心，蓝色星号代表S1或S2频率质心之间的距离，即舒张期或收缩期。
图4
心音待定分量标记
Figure4.
Suspected components of the heart sound waveform
计算这些待定分量之间的距离（distance，Dist）（以符号Dist表示），如式(5)所示：
 |
其中，{flags}表示标记点集合，flagx_i与flagx_i+1分别为两个距离最近的标记点， Fs表示采样率。一个心动周期中S1与S2之间的间隔通常小于58 ms[16]，因此实验中设定了一个58 ms的判定距离。如果待定分量之间的距离小于判定距离，则这些分量会被判定为心音S1或S2。将序列映射到原始心音信号中，根据心音的生理特征，S2与S1之间的信号被定义为舒张期，具有最大的时间间隔[17]。如图5所示，红色虚线为状态序列，根据状态序列的长度进行选取，最长的一段为舒张期，将舒张期前的分量标记为S2，从而实现了对心音信号中S2部分的精准截取。
图5
心音分割标记图
Figure5.
Heart sound segmentation marker chart
1.2 特征提取
心音是由心脏在收缩和开合时血流撞击产生的振动信号组成[18]。心血管疾病导致心脏结构变化，使得心音出现病理特征。本文选取心音S2部分的长记忆性、信息内容和随机性作为统计特征用于分类。
1.2.1 高频分量的分解
DWT将信号分解成多个不同频率的子带，其中高频子带包含信号中的瞬时变化和细节信息[19]。对于心音信号，高频分量包含着心脏收缩、瓣膜开合等短暂且快速的变化，因此选择高频分量可以很好地描述心音信号的动态特性。此外，直接提取高频分量还可以避免低频噪声的影响[20]，去除降噪步骤，进一步提高算法的速度。DWT需要两个基本小波来分析信号：父子波φ(t)和母子波ψ(t)。其中，φ(t)用于提取低频部分，ψ(t)用于提取高频部分。尺度j下的低频振荡v和高频振荡ω的离散形式分别由Sj和Dj表示，如式(6)～式(7)所示：
 |
 |
其中，k表示求和的下界，N表示上界，即从k = 1开始取值，一直取到N。多贝西小波四号（Daubechies-4）用于对S2部分进行三级分解，分解后将312.5 Hz以上的分量认为是高频分量。在尺度j处，Sj表示信号的总体趋势，Dj表示其详细波动。后续的统计特征都是从高频分量Dj提取的。
1.2.2 统计特征提取
统计特征主要包括Hurst指数、LZC与样本熵。这些特征分别用于获取心音高频分量中的长期记忆性、信息内容和规律性[21]，能有效地描述心音信号S2部分的病理特征。
Hurst指数即有偏的随机游走，用于衡量时间序列的长期记忆性。心音信号本质上是一种非平稳时间序列。Hurst值大于0.5，表示为一个具有持续性并且极度相关的时间序列；Hurst值等于0.5，表示为一个不相关的时间序列；Hurst值小于0.5，表示为一个非持续不相关的时间序列，本文用重标极差分析法计算Hurst指数。
给定的时间序列{I1, I2, ···, Iq}，有q个分量，对于任意正整数υ≥1(υ=1, 2, ···, n)，计算其平均值e(υ)与方差S(υ)，如式(8)、式(9)所示:
 |
 |
对q个分量，计算其累计差Z(q)与极差R(υ)，如式(10)、式(11)所示:
 |
 |
引入无量纲比，通过方差S(υ)与极差R(υ)来计算重标极差Rs(υ)，如式(12)所示:
 |
通过计算对应的对数方程，并求取对应的斜率即可得到Hurst系数，其中b是一个常数，H是Hurst指数，如式(13)所示:
 |
LZC可量化时间序列的随机性和复杂性，适用于衡量非线性、非平稳信号的信息内容。LZC对信号中小扰动和噪声不敏感，因此在生理检测中可用于确定心音信号的基本模式与复杂性，正常心音通常以重复模式为特征，而异常心音则以缺乏重复模式为特征。
对于一个给定的信号{x1, x2, …, xN}，将其转化为二进制，转化公式如式(14)所示:
 |
阈值Td为信号{x1, x2, …, xN}平均值。得到yi以后，用c(n)作为yi的复杂度计数器来计算yi中不同字串的数目。c(n)存在着上限，通过b(n)进行归一化得到LZC值C，如式(15)、式(16)所示:
 |
 |
样本熵，是一种衡量时间序列信号复杂性的方法；熵值越低，表示信号具有更高的规律性和可预测性。在心音信号分析中，高的样本熵值，可能表明心音的不规则性或变异性程度较高，并且相较于其它熵测量方法，样本熵的优势在于其稳定性，抗噪、抗干扰能力强。
定义一个长度为m的滑动窗口，从时间序列数据{s1, s2, ···, sN}中提取出所有长度为m的子序列{Sm(1), Sm(2), ···, Sm(N − m + 1)}，其中Sm(i) = si, si + 1, ···, si + m − 1，对于每个子序列，计算其它子序列与该子序列距离小于等于r的数量，表示为Cm(r)。其中，阈值r为原始信号的标准差与自定义权值w的乘积，如式(17)所示:
 |
对于长度为m+1的子序列，计算其它子序列与该子序列距离小于等于r的数量，表示为Cm+1(r)，计算样本熵（以符号SampEn表示），如式(18)所示:
 |
1.3 分类识别
XGBoost，是一种高效的机器学习算法，它采用多个决策树模型的集成来解决二分类和多分类问题。该算法通过构建一系列的决策树模型，并进行结合，每个决策树都基于前面决策树的残差误差构建，因此，能够不断提高分类效果。XGBoost的数学描述如式(19)所示:
 |
其中，t=1, 2, …, T，T表示分类回归树的总数，a=1, 2, …, N，Za为输入，ma为实际值，
为预测值，ft表示第t棵分类回归树。再设置XGBoost的目标函数obj(θ)，如式(20)所示：
 |
其中，
为损失函数，用于衡量实际值ma与预测值
的偏差值。
为正则项，对t个分类回归树的复杂度进行求和。
XGBoost不断优化后，得到第k棵决策树fk的模型fk(Za)，第k个决策树的目标函数H(k)，如式(21)所示：
 |
贝叶斯优化算法，是一种高效的超参数优化算法，它利用先验信息不断优化目标函数，在短时间内找到更好的超参数组合。该算法具有快速找到最优解、强鲁棒性和通用性的优点[22]。在预处理与特征提取后，本文采用贝叶斯优化算法来获取合适的分类模型参数。
2 实验与分析
2.1 数据与说明
实验所用心音数据集采集于云南省阜外心血管病医院，所用采集设备为实验室自行研发的心音信号采集设备，传感器为The ONE心音传感器（Thinklabs Medical LLC Inc.，美国）。招募的志愿者年龄在6个月～18岁，未成年受试者监护人及成年受试者均已签署知情同意书。每个志愿者按5个心脏听诊位置各采集20 s心音信号，采样频率为5 000 Hz。志愿者都经过超声心动仪检测，CHD-PAH患者还经过临床医生在术前与术中检测患者的压力进行确认，以保证数据的可靠性。本次研究所用的心音数据的采集通过了云南大学人体研究材料伦理委员会及云南省阜外心血管病医院医学伦理委员会的批准并授权使用（CHSRE2021008、IRB2020-BG-028）。
从前述采集的心音数据集中选取345名志愿者的1 725例心音数据，志愿者年龄限制18岁以下，这是为了保证训练模型在青少年早期筛查中的有效性。其中，正常、CHD与CHD-PAH心音比例为1:1:1，以保证样本的平衡性。CHD心音包括室间隔缺损、房间隔缺损、动脉导管未闭等常见疾病。将数据按7:2:1的比例划分为训练集、测试集、验证集。
2.2 实验环境
实验所用配置：中央处理器i5-11400H@2.70GHz（Intel Inc.，美国)，独立显卡NVIDIA GeForce RTX3060（Nvidia Inc.，美国）。信号处理软件Matlab2021b（MathWorks Inc.，美国），深度学习框架TensorFlow 2.0（Google Inc.，美国），编程语言Python 3.7（Python Software Foundation Inc.，美国）
2.3 评价指标
使用灵敏度（sensitivity，Se）、特异性（specificity，Sp）、精准率（precision，Pr）、Acc等指标用于评价分类器的性能，如式(22) ～ 式(25)所示：
 |
 |
 |
 |
其中，真阳性（true positive，TP）表示正类别中正确分类样本的数目；真阴性（true negative，TN）表示负类别中正确分类样本的数目；假阴性（false negative，FN）表示负类别中错误分类样本的数目； 假阳性（false positive，FP）表示正类别中错误分类样本的数目。
F1分数，是一种兼顾了Pr与召回率（recall, Re）的分类模型评价指标，如式(26)、式(27)所示：
 |
 |
修正Acc（modified Acc，MAcc）是根据复杂生理信号研究资源网（research resource for complex physiologic signals，PhysioNet）心音挑战赛中所提出的，对Sp和Se指标添加不定权系数0.5，如式(28)所示：
 |
2.4 实验方案设计
本文根据参数设置对结果的影响进行讨论，贝叶斯参数寻优算法对结果的影响如图6所示。为了凸显本文算法速度的优势，与小波降噪、EEMD去噪进行对比，计算处理1 725例心音的运行时间与分类结果，结果如表1所示。
图6
贝叶斯参数优化算法寻优图
Figure6.
Optimization by Bayesian parameter optimization algorithm
表1
预处理运行时间与分类效果
Table1.
Pre-treatment running time and classification results
为了验证统计特征对分类的提升效果，本文还进行了特征消融实验，分别使用三个特征单独进行分类，来区分正常心音、CHD心音与CHD-PAH心音，分类结果总结如表2所示。与文献[9, 13-15]中的声学特征MFCC、MFSC进行对比，其中MFCC采用与人耳结构相似的三角滤波器组，MFSC更注重对信号在高低频内的特征提取。分类模型选择KNN与RF。KNN是一种基于距离度量的非参数分类方法，适用于小规模数据集；而RF通过构建多个决策树来进行分类，能与XGBoost形成对比，分类结果如表3所示。
表2
消融实验分类结果
Table2.
Classification results of ablation experiments
表3
声学特征对比实验分类评价指标
Table3.
Classification evaluation metrics for acoustic features comparison experiment
为了验证本文算法的有效性，选取近两年文献[23-24]提及的心音分类算法，以本文数据集进行训练，对比分类结果。文献[23]对心音进行改进的自适应噪声完全EEMD（the improved complete EEMD with adaptive noise，ICEEMDAN），提取分解后本征模态函数（intrinsic mode function，IMF）的多尺度样本熵，通过逻辑回归（logistic regression，LR)进行分类。文献[24]中对心音进行自适应噪声完备EEMD（complete EEMD with adaptive noise，CEEMDAN），提取IMF的排列熵与支持向量机（support vector machine，SVM）结合后实现分类，分类结果总结如表4所示。
表4
分类算法对比实验分类评价指标
Table4.
Classification evaluation metrics for comparison of classification algorithms
所有实验中，机器学习都通过贝叶斯优化算法进行参数寻优，XGBoost中决策树的数量设为150，学习率设为0.1，最大深度设为5。KNN中k值设为7。RF中决策树的数量设为300，最大深度设置为6。
2.5 实验结果与分析
2.5.1 贝叶斯参数优化实验
如图6所示，对XGBoost的参数进行了贝叶斯优化算法寻优，其中学习率的范围设定在0.1～1.0之间，最大深度设定在3～8，决策树的数量设定在10～600。图6中展示了10组寻优出来的参数以及对应的Acc，参数的顺序为学习率、最大深度、决策树的数量。贝叶斯优化算法在短时间对参数进行寻优，获得了分类Acc较高的一组参数，即学习率0.1，最大深度5，决策树数量150，Acc达到了80.45%。
2.5.2 预处理实验
如表1所示，为本文算法（S2分割+高频分量）与文献[4]、文献[20]中去噪算法的对比实验，小波降噪与EEMD降噪是心音分类中有效的预处理算法。可以看出，由于直接提取高频分量，省略了降噪步骤，本文模型在运行时间上远低于小波降噪与EEMD降噪算法。同时，小波降噪与EEMD降噪过程中可能将一些微弱的病理信息去除，导致分类效果不佳。本文算法直接去除低频分量，虽然会导致低频中部分特征被舍弃，但高频分量中的特征能保证Acc达到80%以上，实现了快速有效的诊断。
2.5.3 消融实验
消融实验结果如表2所示，Hurst指数、LZC和样本熵三个特征都有效地提取了心音高频分量的部分特征，在三分类实验中具有一定的分类效果，其中LZC与样本熵的Acc都达到了70%以上，但同时使用三个特征进行分类的Acc达到了80%以上，各项指标也远高于单独使用，验证了三个特征组合后对分类效果的提升。
2.5.4 三分类实验
三分类实验结果如表3所示，对已经分割出S2部分的心音进行特征提取与分类，将不同的特征和分类器进行组合，得到分类结果与相关指标。基于高频统计特征的XGBoost模型分类效果优于其它算法，表现为：① 高频统计特征相比于MFCC与MFSC在Acc 、MAcc和F1三个评价指标上更优。MFCC的核心是构建一个非均匀分布的频域滤波器组，其特点是低频区域分布密集、高频区域分布稀疏，忽略了信号的相位信息和高频细节，因此在PAH与CHD-PAH的分类中效果不佳。通过对比Acc值，证实了高频统计特征能更好地提取CHD-PAH的病理特征。通过对比MAcc值，表明该特征不易受环境噪声影响。通过对比F1值，得出使用该特征进行分类得到的假阴性更少。② 声学特征MFCC、MFSC与本文特征在XGBoost模型上各项评价指标都高于RF与KNN，体现了XGBoost在小样本多分类任务上的优越性。
2.5.5 对比分类实验
分类算法对比实验结果如表4所示，将文献[23-24]心音分类算法用于本文的CHD-PAH数据集上进行三分类实验，并与本文算法进行比较。
从结果上看，三种算法采用的特征与分类器在三分类上都有一定效果，但本文算法在Acc、Pr、Re和 F1 值上都优于其它方法。① 对比Acc值，证实了分割S2来进行后续特征提取，能更好地提供病理信息，提高了筛查中的检出率。② 通过对比MAcc值，表明直接提取高频分量，去除低频部分的噪声，可提高算法抗噪声能力。③ 通过对比F1值，得出使用S2的高频统计特征相较多尺度样本熵与排列熵，可提取到的病理信息更丰富，得到更少的假阴性病例。④ 采用多个决策树模型的集成，XGBoost在小样本的分类效果也要优于SVM与LR。
综上所述，基于心音S2部分高频分量的统计特征，有效提取了CHD-PAH的病理特征，结合XGBoost分类模型，在三分类综合评价指标上相比现有算法有所提升。
3 结论
本研究的主要目的是设计一个计算机辅助诊断系统，以区分正常、CHD与CHD-PAH心音。本文通过自适应分割出心动周期中S2部分，能更好地凸显CHD-PAH的病理特征，再利用S2的高频分量统计特征来训练贝叶斯优化算法调整的XGBoost模型。本研究的主要特点在于省略了预处理中降噪步骤，通过DWT将原始信号缩短并提取了其高频分量的主要病理特征，这使得特征提取更为快速，整体算法所需时间与资源量低，在三分类上Acc达到80.45%，为CHD-PAH快速辅助诊断提供一种可行的思路。
重要声明
利益冲突声明：本文全体作者均声明不存在利益冲突。
作者贡献声明：杨炫锴负责实验构思与设计、数据采集、数据分析和论文写作；孙静负责实验指导、论文写作指导和数据收集；杨宏波负责医学指导和数据收集；郭涛负责医学指导；潘家华负责医学指导和数据审查；王威廉负责理论指导、实验设计指导、论文整体审查和写作指导。
伦理声明：本文所使用的心音数据得到云南大学人体研究材料伦理委员会以及云南省阜外心血管病医院伦理委员会批准采集并授权使用（CHSRE2021008、IRB2020-BG-028）。
0 引言
肺动脉高压（pulmonary arterial hypertension，PAH）是一种潜在致命性心血管疾病，其特征为肺动脉压力升高，当压力超过一定阈值时，会使得右心房扩大，最终导致右心衰竭[1]。然而，PAH的症状通常在疾病晚期才显现，早期筛查很容易遗漏。PAH的病因复杂，其中约40%是由先天性心脏病（congenital heart disease，CHD）引起的[2]。CHD是一种胎儿时期发育异常导致的心血管畸形，心脏缺陷导致肺循环的压力与容量超负荷，最终形成PAH。CHD相关PAH（CHD-associated PAH，CHD-PAH）由于同时夹杂着两种症状，致死率更高。有创介入式的右心导管术是诊断PAH的金标准，但存在较大的风险。目前无创诊断最有价值的评估方法是超声心动图检查，但在偏远地区，由于缺乏设备和技术人员支持，且该设备不便携带，难以在大规模筛查中推广。
心音是一种机械振动，是心脏血流动力学和心血管系统相互作用的结果。由于心音包含了大量有关心血管系统生理状况的信息，因此可以用于评估心脏的健康状况[3]。计算机智能辅助诊断（computer-aided diagnosis，CAD）通过提取心音信号的病理特征，进行机器学习训练并构建模型用于疾病辅助诊断，能有效配合超声心动图，提高筛查效率，主要步骤分为信号预处理、特征提取、构建分类器等。
预处理阶段主要为降噪与分割。降噪能抑制采集中的呼吸音、杂音等。Dwivedi等[4] 结合稳态小波变换与集合经验模态分解（ensemble empirical mode decomposition，EEMD），来抑制环境音和噪声，以提高后续分类准确率（accuracy，Acc）。Simonneau等[5]研究表明，第一心音（first heart sound，S1）和第二心音（second heart sound，S2）与PAH的严重程度相关。其中PAH患者的S2病理特征更明显，会出现异常亢进与分裂，有助于提高PAH的诊断Acc [6]。但从心音中人工提取S2，困难且复杂[7]，而现有的心音分割算法则没有提取S2，转而通过分割出单个心动周期来研究病理特征，如基于多尺度的希尔伯特包络分割[8]和持续时间隐马尔可夫分割等[9]。因此，如果能高效准确地分割出S2，将是目前有助于促进PAH诊断的一个研究方向。
特征提取的质量会影响分类的结果，Kaddoura等[10]研究心音信号的声学特征，估计心音信号的频谱对数并进行傅里叶逆变换，以创建倒谱特征对PAH与正常心音进行二分类。Alqudah等[11]提取心音信号中S2的能量比作为特征，来区分PAH患者和健康受试者。Ge等[12]提取心动周期和S2的时频域特征等150个特征用于后续分类。
分类阶段，当数据集充足时神经网络能有效提升分类Acc。Wang等[7]基于长短时记忆方法建立了一种小波散射卷积网络，提取受试者的小波散射系数作为特征，来对正常心音与CHD二分类。对于CHD-PAH这样的小样本数据，经典的机器学习能够有效根据信号特征对心音进行分类，不易受小数据量的影响而产生过拟合。已有研究报道，随机森林（random forest，RF）[9]和K近邻法（k-nearest neighbor，KNN）[13]在心音二分类上有一定效果。
目前大多数研究都是针对PAH和正常心音的二分类[5.10-11]，CHD-PAH的多分类问题仍有待改进，比如在预处理阶段，降噪算法有可能滤除病理信息且耗时，达不到实际应用的要求。在心音分割中，大多数分割算法依赖于同步采集心电信号来对心音信号进行辅助分割，对信号采集要求高，只通过心音信号对S2部分进行分割是目前研究的难点。特征提取阶段，目前的研究集中于提取心音的声学和时频域特征，如梅尔倒谱系数（Mel-scale frequency cepstral coefficients，MFCC）[14]和梅尔频谱系数（Mel-frequency spectral coefficient，MFSC）[15]等。然而，声学特征是针对语音信号设计的，而心音信号是由血液流动冲击瓣膜所产生的信号，与语音信号存在一定的差异，并且在听诊中，PAH和CHD-PAH仅通过听觉是很难区分，因为它们之间的声学与时频特征相似度高，分类效果不理想。此外，对于CHD-PAH的多分类诊断依赖于大量特征[12]，特征提取耗时长，不利于大规模的早期筛查任务。在分类器阶段，由于CHD-PAH的数据量少，神经网络在少量数据上容易出现过拟合的问题，传统机器学习更适用于小样本数据分类。
鉴于上述现状，本文旨在建立一个快速有效的CHD-PAH诊断算法。在预处理阶段，不进行降噪，使用自适应的双阈值分割法得到心音S2部分，再提取S2高频分量中的统计特征，包括赫斯特（Hurst）指数、勒佩尔-齐夫（Lempel-Ziv）复杂度（Lempel-Ziv complexity，LZC）和样本熵（sample entropy，SampEn）。与传统方法不同的是，本文算法专注于研究心音的高频特性而非声学相关特性，这些高频分量被用来描述短而突然的变化特征杂音，能更有效地凸显出PAH患者S2亢进、分裂的特征。最后，采用贝叶斯优化的极端梯度提升算法（extreme gradient boosting，XGBoost）对正常、CHD和CHD-PAH进行三分类诊断。
1 方法
本文算法通过以下几个步骤对CHD-PAH进行诊断：① 使用一种自适应分割方法来截取出心音信号中S2部分。② 对S2部分进行离散小波变换（discrete wavelet transformation，DWT）分解，提取出高频分量。③ 使用Hurst指数、LZC和样本熵进行复杂性度量，提取出S2高频分量的统计特征。④ 特征通过XGboost分类器进行CHD-PAH的三分类诊断。CHD-PAH分类流程如图1所示。
图1
CHD-PAH流程图
Figure1.
CHD-PAH flow chart
1.1 预处理
心脏跳动是一个近似于周期性的过程，通常分为S1、收缩期、S2和舒张期四个部分。在听诊PAH患者时，S2中肺部成分加重是其确诊的主要特征之一[6]。如图2所示，局部方框圈出部分波形即为S2，PAH患者表现为S2部分出现亢进与异常分裂。相比于整个心动周期，仅使用S2部分即能有效地凸显出PAH的特征[13]。并且，将原始信号分割为单独的S2部分后，还能提高特征提取速度。
图2
CHD-PAH的S2部分出现异常
Figure2.
Anomalies in the S2 part of CHD-PAH
本文采用了一种阈值检测算法分割S2，如图1所示。首先对心音信号进行短时傅里叶变换（short-time Fourier transform，STFT），并使用汉明窗进行分帧，如式(1)所示：
|
其中， t为时间，f为频率，j为虚数，π为圆周率，x(τ) 为t = τ时的心音信号，ω(τ − t)为窗函数。心音 x(τ)被分帧后，提取每一帧的短时能量Ei与频谱扩散Si，如式(2)～式(3)所示：
|
|
其中，xi(n)，n = 1, ···, N表示心音信号的第i帧，N为帧数长度，fi为第i帧的频率，u1为频谱中值，si为第i帧的频谱值， b1与b2为频谱边界值。所得Ei与Si，如图3所示。
图3
短时能量与频谱扩散
Figure3.
Short-time energy and spectral diffusion
计算短时能量与频谱扩散对应的直方图与直方图阈值TE与TS，阈值计算如式(4)所示：
|
其中，M1、M2为直方图中前两个最大值的位置，W为设定值。根据阈值对序列进行判断并标记。其中，满足Ei ≥ TE的为标记点1（flag1），满足Si ≥ TS的为标记点2（flag2），这些标记被认为是心音S1或S2分量，称之为待定分量。如图4所示，根据阈值判断并标记序列，黑色波形是原始心音信号的频率质心，绿色波形是待定分量，即被标记为心音S1或心音S2的频率质心，蓝色星号代表S1或S2频率质心之间的距离，即舒张期或收缩期。
图4
心音待定分量标记
Figure4.
Suspected components of the heart sound waveform
计算这些待定分量之间的距离（distance，Dist）（以符号Dist表示），如式(5)所示：
|
其中，{flags}表示标记点集合，flagx_i与flagx_i+1分别为两个距离最近的标记点， Fs表示采样率。一个心动周期中S1与S2之间的间隔通常小于58 ms[16]，因此实验中设定了一个58 ms的判定距离。如果待定分量之间的距离小于判定距离，则这些分量会被判定为心音S1或S2。将序列映射到原始心音信号中，根据心音的生理特征，S2与S1之间的信号被定义为舒张期，具有最大的时间间隔[17]。如图5所示，红色虚线为状态序列，根据状态序列的长度进行选取，最长的一段为舒张期，将舒张期前的分量标记为S2，从而实现了对心音信号中S2部分的精准截取。
图5
心音分割标记图
Figure5.
Heart sound segmentation marker chart
1.2 特征提取
心音是由心脏在收缩和开合时血流撞击产生的振动信号组成[18]。心血管疾病导致心脏结构变化，使得心音出现病理特征。本文选取心音S2部分的长记忆性、信息内容和随机性作为统计特征用于分类。
1.2.1 高频分量的分解
DWT将信号分解成多个不同频率的子带，其中高频子带包含信号中的瞬时变化和细节信息[19]。对于心音信号，高频分量包含着心脏收缩、瓣膜开合等短暂且快速的变化，因此选择高频分量可以很好地描述心音信号的动态特性。此外，直接提取高频分量还可以避免低频噪声的影响[20]，去除降噪步骤，进一步提高算法的速度。DWT需要两个基本小波来分析信号：父子波φ(t)和母子波ψ(t)。其中，φ(t)用于提取低频部分，ψ(t)用于提取高频部分。尺度j下的低频振荡v和高频振荡ω的离散形式分别由Sj和Dj表示，如式(6)～式(7)所示：
|
|
其中，k表示求和的下界，N表示上界，即从k = 1开始取值，一直取到N。多贝西小波四号（Daubechies-4）用于对S2部分进行三级分解，分解后将312.5 Hz以上的分量认为是高频分量。在尺度j处，Sj表示信号的总体趋势，Dj表示其详细波动。后续的统计特征都是从高频分量Dj提取的。
1.2.2 统计特征提取
统计特征主要包括Hurst指数、LZC与样本熵。这些特征分别用于获取心音高频分量中的长期记忆性、信息内容和规律性[21]，能有效地描述心音信号S2部分的病理特征。
Hurst指数即有偏的随机游走，用于衡量时间序列的长期记忆性。心音信号本质上是一种非平稳时间序列。Hurst值大于0.5，表示为一个具有持续性并且极度相关的时间序列；Hurst值等于0.5，表示为一个不相关的时间序列；Hurst值小于0.5，表示为一个非持续不相关的时间序列，本文用重标极差分析法计算Hurst指数。
给定的时间序列{I1, I2, ···, Iq}，有q个分量，对于任意正整数υ≥1(υ=1, 2, ···, n)，计算其平均值e(υ)与方差S(υ)，如式(8)、式(9)所示:
|
|
对q个分量，计算其累计差Z(q)与极差R(υ)，如式(10)、式(11)所示:
|
|
引入无量纲比，通过方差S(υ)与极差R(υ)来计算重标极差Rs(υ)，如式(12)所示:
|
通过计算对应的对数方程，并求取对应的斜率即可得到Hurst系数，其中b是一个常数，H是Hurst指数，如式(13)所示:
|
LZC可量化时间序列的随机性和复杂性，适用于衡量非线性、非平稳信号的信息内容。LZC对信号中小扰动和噪声不敏感，因此在生理检测中可用于确定心音信号的基本模式与复杂性，正常心音通常以重复模式为特征，而异常心音则以缺乏重复模式为特征。
对于一个给定的信号{x1, x2, …, xN}，将其转化为二进制，转化公式如式(14)所示:
|
阈值Td为信号{x1, x2, …, xN}平均值。得到yi以后，用c(n)作为yi的复杂度计数器来计算yi中不同字串的数目。c(n)存在着上限，通过b(n)进行归一化得到LZC值C，如式(15)、式(16)所示:
|
|
样本熵，是一种衡量时间序列信号复杂性的方法；熵值越低，表示信号具有更高的规律性和可预测性。在心音信号分析中，高的样本熵值，可能表明心音的不规则性或变异性程度较高，并且相较于其它熵测量方法，样本熵的优势在于其稳定性，抗噪、抗干扰能力强。
定义一个长度为m的滑动窗口，从时间序列数据{s1, s2, ···, sN}中提取出所有长度为m的子序列{Sm(1), Sm(2), ···, Sm(N − m + 1)}，其中Sm(i) = si, si + 1, ···, si + m − 1，对于每个子序列，计算其它子序列与该子序列距离小于等于r的数量，表示为Cm(r)。其中，阈值r为原始信号的标准差与自定义权值w的乘积，如式(17)所示:
|
对于长度为m+1的子序列，计算其它子序列与该子序列距离小于等于r的数量，表示为Cm+1(r)，计算样本熵（以符号SampEn表示），如式(18)所示:
|
1.3 分类识别
XGBoost，是一种高效的机器学习算法，它采用多个决策树模型的集成来解决二分类和多分类问题。该算法通过构建一系列的决策树模型，并进行结合，每个决策树都基于前面决策树的残差误差构建，因此，能够不断提高分类效果。XGBoost的数学描述如式(19)所示:
|
其中，t=1, 2, …, T，T表示分类回归树的总数，a=1, 2, …, N，Za为输入，ma为实际值，为预测值，ft表示第t棵分类回归树。再设置XGBoost的目标函数obj(θ)，如式(20)所示：
|
其中，为损失函数，用于衡量实际值ma与预测值的偏差值。为正则项，对t个分类回归树的复杂度进行求和。
XGBoost不断优化后，得到第k棵决策树fk的模型fk(Za)，第k个决策树的目标函数H(k)，如式(21)所示：
|
贝叶斯优化算法，是一种高效的超参数优化算法，它利用先验信息不断优化目标函数，在短时间内找到更好的超参数组合。该算法具有快速找到最优解、强鲁棒性和通用性的优点[22]。在预处理与特征提取后，本文采用贝叶斯优化算法来获取合适的分类模型参数。
2 实验与分析
2.1 数据与说明
实验所用心音数据集采集于云南省阜外心血管病医院，所用采集设备为实验室自行研发的心音信号采集设备，传感器为The ONE心音传感器（Thinklabs Medical LLC Inc.，美国）。招募的志愿者年龄在6个月～18岁，未成年受试者监护人及成年受试者均已签署知情同意书。每个志愿者按5个心脏听诊位置各采集20 s心音信号，采样频率为5 000 Hz。志愿者都经过超声心动仪检测，CHD-PAH患者还经过临床医生在术前与术中检测患者的压力进行确认，以保证数据的可靠性。本次研究所用的心音数据的采集通过了云南大学人体研究材料伦理委员会及云南省阜外心血管病医院医学伦理委员会的批准并授权使用（CHSRE2021008、IRB2020-BG-028）。
从前述采集的心音数据集中选取345名志愿者的1 725例心音数据，志愿者年龄限制18岁以下，这是为了保证训练模型在青少年早期筛查中的有效性。其中，正常、CHD与CHD-PAH心音比例为1:1:1，以保证样本的平衡性。CHD心音包括室间隔缺损、房间隔缺损、动脉导管未闭等常见疾病。将数据按7:2:1的比例划分为训练集、测试集、验证集。
2.2 实验环境
实验所用配置：中央处理器i5-11400H@2.70GHz（Intel Inc.，美国)，独立显卡NVIDIA GeForce RTX3060（Nvidia Inc.，美国）。信号处理软件Matlab2021b（MathWorks Inc.，美国），深度学习框架TensorFlow 2.0（Google Inc.，美国），编程语言Python 3.7（Python Software Foundation Inc.，美国）
2.3 评价指标
使用灵敏度（sensitivity，Se）、特异性（specificity，Sp）、精准率（precision，Pr）、Acc等指标用于评价分类器的性能，如式(22) ～ 式(25)所示：
|
|
|
|
其中，真阳性（true positive，TP）表示正类别中正确分类样本的数目；真阴性（true negative，TN）表示负类别中正确分类样本的数目；假阴性（false negative，FN）表示负类别中错误分类样本的数目； 假阳性（false positive，FP）表示正类别中错误分类样本的数目。
F1分数，是一种兼顾了Pr与召回率（recall, Re）的分类模型评价指标，如式(26)、式(27)所示：
|
|
修正Acc（modified Acc，MAcc）是根据复杂生理信号研究资源网（research resource for complex physiologic signals，PhysioNet）心音挑战赛中所提出的，对Sp和Se指标添加不定权系数0.5，如式(28)所示：
|
2.4 实验方案设计
本文根据参数设置对结果的影响进行讨论，贝叶斯参数寻优算法对结果的影响如图6所示。为了凸显本文算法速度的优势，与小波降噪、EEMD去噪进行对比，计算处理1 725例心音的运行时间与分类结果，结果如表1所示。
图6
贝叶斯参数优化算法寻优图
Figure6.
Optimization by Bayesian parameter optimization algorithm
表1
预处理运行时间与分类效果
Table1.
Pre-treatment running time and classification results
为了验证统计特征对分类的提升效果，本文还进行了特征消融实验，分别使用三个特征单独进行分类，来区分正常心音、CHD心音与CHD-PAH心音，分类结果总结如表2所示。与文献[9, 13-15]中的声学特征MFCC、MFSC进行对比，其中MFCC采用与人耳结构相似的三角滤波器组，MFSC更注重对信号在高低频内的特征提取。分类模型选择KNN与RF。KNN是一种基于距离度量的非参数分类方法，适用于小规模数据集；而RF通过构建多个决策树来进行分类，能与XGBoost形成对比，分类结果如表3所示。
表2
消融实验分类结果
Table2.
Classification results of ablation experiments
表3
声学特征对比实验分类评价指标
Table3.
Classification evaluation metrics for acoustic features comparison experiment
为了验证本文算法的有效性，选取近两年文献[23-24]提及的心音分类算法，以本文数据集进行训练，对比分类结果。文献[23]对心音进行改进的自适应噪声完全EEMD（the improved complete EEMD with adaptive noise，ICEEMDAN），提取分解后本征模态函数（intrinsic mode function，IMF）的多尺度样本熵，通过逻辑回归（logistic regression，LR)进行分类。文献[24]中对心音进行自适应噪声完备EEMD（complete EEMD with adaptive noise，CEEMDAN），提取IMF的排列熵与支持向量机（support vector machine，SVM）结合后实现分类，分类结果总结如表4所示。
表4
分类算法对比实验分类评价指标
Table4.
Classification evaluation metrics for comparison of classification algorithms
所有实验中，机器学习都通过贝叶斯优化算法进行参数寻优，XGBoost中决策树的数量设为150，学习率设为0.1，最大深度设为5。KNN中k值设为7。RF中决策树的数量设为300，最大深度设置为6。
2.5 实验结果与分析
2.5.1 贝叶斯参数优化实验
如图6所示，对XGBoost的参数进行了贝叶斯优化算法寻优，其中学习率的范围设定在0.1～1.0之间，最大深度设定在3～8，决策树的数量设定在10～600。图6中展示了10组寻优出来的参数以及对应的Acc，参数的顺序为学习率、最大深度、决策树的数量。贝叶斯优化算法在短时间对参数进行寻优，获得了分类Acc较高的一组参数，即学习率0.1，最大深度5，决策树数量150，Acc达到了80.45%。
2.5.2 预处理实验
如表1所示，为本文算法（S2分割+高频分量）与文献[4]、文献[20]中去噪算法的对比实验，小波降噪与EEMD降噪是心音分类中有效的预处理算法。可以看出，由于直接提取高频分量，省略了降噪步骤，本文模型在运行时间上远低于小波降噪与EEMD降噪算法。同时，小波降噪与EEMD降噪过程中可能将一些微弱的病理信息去除，导致分类效果不佳。本文算法直接去除低频分量，虽然会导致低频中部分特征被舍弃，但高频分量中的特征能保证Acc达到80%以上，实现了快速有效的诊断。
2.5.3 消融实验
消融实验结果如表2所示，Hurst指数、LZC和样本熵三个特征都有效地提取了心音高频分量的部分特征，在三分类实验中具有一定的分类效果，其中LZC与样本熵的Acc都达到了70%以上，但同时使用三个特征进行分类的Acc达到了80%以上，各项指标也远高于单独使用，验证了三个特征组合后对分类效果的提升。
2.5.4 三分类实验
三分类实验结果如表3所示，对已经分割出S2部分的心音进行特征提取与分类，将不同的特征和分类器进行组合，得到分类结果与相关指标。基于高频统计特征的XGBoost模型分类效果优于其它算法，表现为：① 高频统计特征相比于MFCC与MFSC在Acc 、MAcc和F1三个评价指标上更优。MFCC的核心是构建一个非均匀分布的频域滤波器组，其特点是低频区域分布密集、高频区域分布稀疏，忽略了信号的相位信息和高频细节，因此在PAH与CHD-PAH的分类中效果不佳。通过对比Acc值，证实了高频统计特征能更好地提取CHD-PAH的病理特征。通过对比MAcc值，表明该特征不易受环境噪声影响。通过对比F1值，得出使用该特征进行分类得到的假阴性更少。② 声学特征MFCC、MFSC与本文特征在XGBoost模型上各项评价指标都高于RF与KNN，体现了XGBoost在小样本多分类任务上的优越性。
2.5.5 对比分类实验
分类算法对比实验结果如表4所示，将文献[23-24]心音分类算法用于本文的CHD-PAH数据集上进行三分类实验，并与本文算法进行比较。
从结果上看，三种算法采用的特征与分类器在三分类上都有一定效果，但本文算法在Acc、Pr、Re和 F1 值上都优于其它方法。① 对比Acc值，证实了分割S2来进行后续特征提取，能更好地提供病理信息，提高了筛查中的检出率。② 通过对比MAcc值，表明直接提取高频分量，去除低频部分的噪声，可提高算法抗噪声能力。③ 通过对比F1值，得出使用S2的高频统计特征相较多尺度样本熵与排列熵，可提取到的病理信息更丰富，得到更少的假阴性病例。④ 采用多个决策树模型的集成，XGBoost在小样本的分类效果也要优于SVM与LR。
综上所述，基于心音S2部分高频分量的统计特征，有效提取了CHD-PAH的病理特征，结合XGBoost分类模型，在三分类综合评价指标上相比现有算法有所提升。
3 结论
本研究的主要目的是设计一个计算机辅助诊断系统，以区分正常、CHD与CHD-PAH心音。本文通过自适应分割出心动周期中S2部分，能更好地凸显CHD-PAH的病理特征，再利用S2的高频分量统计特征来训练贝叶斯优化算法调整的XGBoost模型。本研究的主要特点在于省略了预处理中降噪步骤，通过DWT将原始信号缩短并提取了其高频分量的主要病理特征，这使得特征提取更为快速，整体算法所需时间与资源量低，在三分类上Acc达到80.45%，为CHD-PAH快速辅助诊断提供一种可行的思路。
重要声明
利益冲突声明：本文全体作者均声明不存在利益冲突。
作者贡献声明：杨炫锴负责实验构思与设计、数据采集、数据分析和论文写作；孙静负责实验指导、论文写作指导和数据收集；杨宏波负责医学指导和数据收集；郭涛负责医学指导；潘家华负责医学指导和数据审查；王威廉负责理论指导、实验设计指导、论文整体审查和写作指导。
伦理声明：本文所使用的心音数据得到云南大学人体研究材料伦理委员会以及云南省阜外心血管病医院伦理委员会批准采集并授权使用（CHSRE2021008、IRB2020-BG-028）。
