基于功率譜密度與卷積神經(jīng)網(wǎng)絡的心音分類

許春冬，辛鵬麗，周 靜，應冬文，2

1.江西理工大學 信息工程學院，江西 贛州341000

2.中國科學院大學 電子電氣與通信工程學院，北京100049

心血管疾病是全球患病率及疾病致死率的主要因素之一，在影響人類社會活動的同時，更嚴重危害著生命健康[1]。針對心血管疾病診斷的難題，出現(xiàn)了一系列輔助診斷方法，如心電分析、心音分析、CT檢測等[2]。其中，心音信號是由心臟進行周期性泵血運動產(chǎn)生的振動，經(jīng)人體組織傳導于體表形成的一類聲學信號，能夠反映心房、心室、瓣膜、動脈血管壁等結構的工作狀況，這也就使得通過心音來分析心血管疾病是可行的[3]。心音聽診是心音信號分析中常用的方法，被廣泛應用于臨床診斷中。然而，精確的聽診往往依賴于醫(yī)師的長期訓練經(jīng)驗，而這一過程又耗時耗力，因此需要一個用于心音分析的計算機輔助工具來協(xié)助診斷心血管疾病[3-4]。相比于其他輔助診斷方法，心音信號分析具備操作簡便、完全無損、成本低廉、便于遠程診斷、效率高等特點，且可檢測出初期的心血管疾病，能夠起到早檢測早治療的關鍵性作用[5]。

心音分析系統(tǒng)中，對正常與異常心音信號的分類方法研究最多，主要用于輔助醫(yī)師完成初步診斷，實現(xiàn)精確檢測，以節(jié)省醫(yī)療資源和成本[5-6]。在心音信號正常與否的判定問題中，研究者們提出了一系列方法。Gupta等人[7]提出了結合同態(tài)濾波和K 均值聚類的心音信號的小波分析分割方法，將分割結果形成特征向量，在兩類異常和一類正常心音分類中準確率達到97%，但當異常心音種類增多時分類精度明顯下降；Dokur 等人[3]提出通過離散小波變換加窗的方法完成心音的分割與特征提取，采用增量自組織映射網(wǎng)絡進行分類，但在處理噪聲環(huán)境復雜的心音時其分割及特征識別結果并不理想；Wang 等人[8]提出采用25～400 Hz 的巴特沃斯濾波器結合基于邏輯回歸的隱半馬爾科夫模型提高分割精度，最終獲得74.9%的分類結果；Zhang等人[9]提出了一種基于縮放圖譜和張量分解的特征提取方法，并通過支持向量機完成分類，在數(shù)據(jù)集C上取得了90%的分類結果。上述研究存在一個共同的問題，即在提取特征前要先有效地完成分割。然而，分割問題一直是心音信號分析中極具挑戰(zhàn)性的難題，到目前為止還沒有一種完美的分割方案。因此一些學者提出了直接提取整個心音段或心動周期的特征進行識別，從而規(guī)避因分割問題給系統(tǒng)帶來的誤差[10]。郭興明等人[11]提出了基于經(jīng)驗模式分解和關聯(lián)維數(shù)的心音特征分類方法，但心音信號具備多樣性，僅僅依靠關聯(lián)系數(shù)無法把握心音信號的本質(zhì)特征；汪晶[10]提出了多閾值融合心音遞歸圖和局部線性嵌入算法降維小波系數(shù)自相關特征的特征提取算法，但其時耗性高，且處理復雜病理心音效果不理想；Deng等人[12]提出將子帶包絡與自相關特征融合獲取擴散圖統(tǒng)一特征的分類方法，但在噪聲干擾較強時容易提取到噪聲信號；Singh-Μiller 等人[13]提出通過特征譜圖的均值、方差等特征來完成分類任務，同樣容易受噪聲干擾。

基于上述分類問題中遇到的難題，本文提出了一種基于功率譜密度（Power Spectral Density，PSD）時頻特征的正常與異常心音分類方法，通過有效的預處理和特征提取過程，采用卷積神經(jīng)網(wǎng)絡（Convolutional Neural Network，CNN）分類器，在Challenge 2016 數(shù)據(jù)集中獲得了良好的評價指標。

1 算法框架

本文提出了一種基于功率譜密度時頻特征的正常與異常心音分類方法，主要包括預處理、檢測心動周期并分段、提取PSD時頻特征矩陣、特征矩陣填補、構建分類器、分類器訓練與測試等，如圖1所示。首先，采用小波降噪的方法抑制心音中噪聲的干擾；然后，采用自相關法檢測心動周期，并按相關包絡分段出第一個心動周期，并對剩余部分循環(huán)提取，避免長時采集心音的自相關包絡帶來的心動周期誤差；獲取心動周期集后，即可提取心動周期的PSD時頻特征，并通過雙線性插值法將周期略短的心動周期PSD 特征進行補齊處理。根據(jù)Challenge 2016 的訓練集與測試集數(shù)據(jù)，獲取訓練特征集和測試特征集，送入構建的CNN 分類器進行迭代訓練，并通過測試進行泛化能力檢驗，從而獲得有效的正常與異常心音分類模型。

圖1 算法原理圖

2 PSD時頻特征的提取

針對心音分類，學者們提出了大量基于時域特征的分析方法[4]。然而，與時域分析相比，頻域分析往往更具優(yōu)勢，因為特定信號具有較為穩(wěn)定的頻帶，而時域卻呈現(xiàn)較強的非平穩(wěn)性[3-4]。頻域分析是一種較好的標準化分析工具，在許多臨床和生理現(xiàn)象中都得到了應用[4，9]。正常心音與異常心音在頻譜的分布上存在一定的差異，主要源于三尖瓣、二尖瓣、主動脈、心肌等結構工作狀況不同。在正常心音中能夠聽到清晰的第一心音（s1，10～140 Hz）和第二心音（s2，10～200 Hz），在部分兒童心音或異常心音中可聽到第三心音（s3，10～70 Hz），第四心音（s4，10～70 Hz）主要出現(xiàn)在異常心音中[14]。文獻[15]指出，異常心音在收縮期、舒張期會產(chǎn)生大量心雜音，且頻帶范圍較正常心音更高，部分異常心音頻率可達600 Hz左右。功率譜密度能夠聯(lián)合時頻域特征，較好地描述心音信號的時頻分布特征，故本文將通過功率譜密度分析對心音特征進行量化。

2.1 心音信號的預處理

心音信號來源于心臟泵血的機械振動，極易受到各類噪聲的干擾。因此，在對心音信號進行分析前，需進行降噪處理。小波降噪是心音分析中應用最廣的降噪方法，本文采用了小波降噪的方法，對心音信號進行了降噪處理。小波降噪的原理是通過小波系數(shù)來區(qū)分各尺度內(nèi)的噪聲與信號，其關鍵在于閾值函數(shù)的構建，本課題組對小波降噪進行了深入的研究，具體算法請參考文獻[16]。圖2 給出了隨機選取的四條心音記錄（均來自Challenge 2016 數(shù)據(jù)庫[14]）采用小波降噪方法降噪后的結果，可以看出該方法有效地抑制了噪聲，提升了心音信號的可分析性。

2.2 心動周期的檢測

圖2 降噪結果

心音信號具有顯著的周期性特征，其正常與否的特點通常在心動周期中能夠得到完整的體現(xiàn)。在功率譜特征的提取中，以心動周期為分析對象，獲取心臟每個循環(huán)供血過程中的細節(jié)特征，來判斷心臟工作的正常與異常。準確檢測心動周期的長度是檢測心音心動周期的前提。文獻[5]指出，通過自相關的方法能夠較為準確地獲取到心動周期的長度。然而，在長時采集心音段的自相關包絡中，由于后期的平滑處理，隨著時間的延長，無法保留心動周期的變異性特征，甚至使得用于區(qū)分心動周期的包絡特征變得模糊，如相關包絡變得矮小不利于波峰提取，如圖3所示。通過將心音信號劃分為一定長度的分析窗可以避免此類問題。本文采用自相關法檢測心動周期T[4-5]，并根據(jù)T 將心音信號分段為心動周期，其流程如下：

（1）心動周期長度一般在0.6～1.2 s 之間，將待處理的心音以4 s為時長取出心音的第一個時窗。

（2）按照式（1）計算時窗內(nèi)心音的自相關：

其中，mT為時窗內(nèi)第mT幀，mT=0,1,…,N,N 為時窗內(nèi)總的幀數(shù)，F(xiàn)r(lenT) 為短時平均能量，F(xiàn)rame 為幀長。在獲取心動周期時取自相關大于等于0部分即可，并做中值濾波處理，平滑包絡曲線。圖3展示了一條預處理后的心音段自相關包絡的提取結果。

（3）提取自相關包絡的波峰，記0 處的第一個波峰為P1，以0.5～1.5 s 采樣點為區(qū)間搜索下一個最高波峰P2，則P1～P2的時段即為第一個心動周期T1。

（4）跳過P2前的時段，以P2為新的波峰起點，重復步驟（3），搜索下一個波峰P3，獲取下一個新的周期T2。

（5）從心音中去除已檢測的前兩個心動周期部分，獲取新的時窗，并重復步驟（2）、（3）、（4），直到需要被檢測的剩余心音不足一個時窗的長度。

2.3 功率譜密度時頻分布矩陣的提取

圖3 心動周期的檢測

常用的功率譜估計方法主要有周期圖法、間接估計法、Bartlett估計法、Welch估計法四類[17]。周期圖法又被稱作直接估計法，其過程為先將待估計的信號x(n)進行快速傅里葉變換（Fast Fourier Transformation，F(xiàn)FT），然后取幅值平方并做平均得到估計值；Bartlett估計法是在周期圖法的基礎上按照分段求周期圖再取平均的方式進行優(yōu)化改進，改善了譜曲線的不穩(wěn)定和分辨率問題；Welch 估計法則是在Bartlett 估計法的基礎上增加了加窗分幀和幀移，避免譜估計為負并減小方差[18]。本文采用Welch估計法來估計心音信號的PSD，其步驟如下：

（1）根據(jù)fs=2 000 Hz 的采樣率，取窗長為win=30 ms，一是根據(jù)短時平穩(wěn)性使每一窗心音的PSD 具備平穩(wěn)性，二是使得心音PSD在時間域上具備較高的分辨率。幀移取inc=15 ms，使得PSD在時間域上具備較好的平滑性，按照窗長對信號進行分幀：

其中，m 為幀內(nèi)時域樣點，l 表示第l 幀，M 為總幀數(shù)，d1為窗函數(shù)。

（2）通過Welch 算法估計心音幀的功率譜密度PSD(w)：

其中，d2表示數(shù)據(jù)窗口，U 為歸一化因子。

（3）取對數(shù)獲得功率譜密度- -- ---PSD(w)：

（4）對功率譜密度進行正向歸一化處理：

上述過程中，窗函數(shù)的選擇對功率譜存在一定的影響，故測量了矩形窗（Boxcar）、海明窗（Hamming）、布萊克曼窗（Blackman）三類窗函數(shù)下功率譜密度的結果。圖4 給出了三類窗函數(shù)下一組心音信號PSD 均值的估計結果?？梢钥闯觯啾菳oxcar 窗和Blackman 窗，Hamming 窗能夠有效抑制低頻成分，突出基礎心音頻帶，更為適用于心音信號的PSD估計。采用Hamming窗提取的正向歸一化PSD 特征如圖5 所示。因此，d1和d2為Hamming窗。

圖4 不同窗函數(shù)下正向歸一化PSD均值的估計

2.4 雙線性插值法填補功率譜密度特征矩陣

第2.2節(jié)采用自相關法檢測心動周期并將心音信號分段為心動周期集，但心音具備多樣性，且心動周期間也具備心率變異性，這就使得分段出的心動周期長度并不一致，即每個心動周期幀數(shù)不一致。因此，心動周期集所求取出的功率譜密度特征矩陣大小也是不同的，而分類器輸入端的大小一般是固定的，故要求功率譜密度矩陣的大小一致。文獻[9]指出，可通過雙線性插值法來填補心動周期頻譜圖，本文采用此方法填補心動周期功率譜密度矩陣，使得各心動周期功率譜密度矩陣大小一致。

圖5 正向歸一化PSD特征

雙線性插值法心動周期功率譜填補的目的不是要在時間刻度上延長獲取更多的矩陣元素，而是在原始時間刻度內(nèi)增加功率譜密度幀數(shù)，從而提升時間維度的分辨率，這樣可以保持心動周期內(nèi)的頻譜分布基本保持不變。心音功率譜密度矩陣的雙線性插值法填補流程主要由以下四個步驟組成[9]：

（2）第i 個時間尺度點和第j 個頻率尺度點的插值因子為：

其中，round 表示四舍五入取整運算。

第i 個時間尺度點和第j 個頻率尺度點插值因子的絕對值偏差為：

（4）重復步驟（2）、（3），直至完成所有插值點的計算。

在心音信號的功率譜密度時頻矩陣填補中，只需填補時間尺度點，無需填充頻率尺度點。這是因為頻率尺度只與短時傅里葉變換的窗長有關，固定窗長其頻率尺度則一致；而時間尺度與心動周期長度相關，這恰是需要填補的。圖6給出了a0001中某一心動周期的填補結果，可看出時間軸長度不變，但增加了幀數(shù)。

圖6 雙線性法填充特征矩陣

3 構建CNN分類器與數(shù)據(jù)集

用于心音的分類器主要包括支持向量機（Support Vector Μachine，SVΜ）分類器、反向傳播神經(jīng)（Back Propagation，BP）網(wǎng)絡分類器、自組織映射（Self-Organizing Μapping，SOΜ）網(wǎng)絡分類器、隱馬爾科夫模型（Hidden Μarkov Μodel，HΜΜ）分類器、循環(huán)神經(jīng)網(wǎng)絡（Recurrent Neural Network，RNN）分類器、長短時記憶網(wǎng)絡（Long Short-Term Μemory Network，LSTΜ）分類器和卷積神經(jīng)網(wǎng)絡（CNN）分類器。其中，SVΜ 分類器、BP神經(jīng)網(wǎng)絡分類器不具備自主優(yōu)化特征的性能[9]；HΜΜ分類器多借助于高斯混合模型進行分類，當正常與異常心音時域分布類似時容易出現(xiàn)分類誤判[5，8]；RNN 分類器處理時間序列效果較好，但存在梯度消失和爆炸的問題；LSTΜ網(wǎng)絡在RNN基礎上，增加了門單元，解決了梯度消失和爆炸的問題；CNN 分類器相比前幾類分類器模型，能夠較好地處理多維數(shù)據(jù)，挖掘局部結構特點[19-20]。CNN 已經(jīng)在圖像處理領域取得了顯著成績，而心音PSD 時頻矩陣與圖像相同為二維矩陣，因此CNN 分類器更為適用。

CNN 是一種融入卷積結構的深度前饋神經(jīng)網(wǎng)絡，主要由輸入層、卷積層、激勵層、池化層、全連接層、輸出層組成[19，21]。其中卷積層、激勵層與池化層可進行交替連接，能夠有效地降低數(shù)據(jù)維度，提高模型的容錯能力。在CNN 分類器構建過程中，通過在有限范圍內(nèi)進行隨機搜索，能夠確定較為合適的卷積層數(shù)和卷積核大小。構建的CNN分類器的結構如圖7所示。

實驗中將99×41 的特征矩陣去除頻率尺度點為0、時間尺度點為0 所對應的行與列，構成98×40 的輸入矩陣，通過分類器并進行二分類。構建的CNN 分類器共計兩次卷積運算，第一卷積層的卷積核大小為5×5，數(shù)量為32，卷積水平步長為1，最大池化層為2×2，選擇ReLU函數(shù)為激勵函數(shù)，則第一次卷積運算后每個卷積核中的數(shù)據(jù)被壓縮為94×36 的特征圖，池化后特征圖大小為47×18；第二卷積層卷積核大小為2×3，數(shù)量為64，卷積水平步長為1，最大池化層為1×4，激活函數(shù)仍選擇ReLU，則第二次卷積運算后每個卷積核中的數(shù)據(jù)被壓縮為46×16 的特征圖，池化后特征圖大小為23×8；共三層全連接層，第一全連接層為1×1 024，第二全連接層為1×512，第三全連接層為1×256；輸出層為1×2，輸出結果為01 表示正常，輸出結果為10 表示異常。部分訓練參數(shù)設置如表1 所示，其中B 為下采樣層的map 的權值，Dropout為防止過擬合的trick。

表1 網(wǎng)絡參數(shù)設置

圖7 CNN分類器網(wǎng)絡結構

通過訓練數(shù)據(jù)集對CNN 分類器進行迭代訓練，并通過分類正確率和損失函數(shù)兩類指標對分類器進行評價。分類器為二分類器，損失函數(shù)采用常規(guī)的交叉熵損失函數(shù)即可。正確率ACC 和損失函數(shù)LOSS 評價指標的定義如式（19）、（20）所示：

其中，TP 為正常心音分類正確數(shù)，TN 為異常樣本分類正確數(shù)，F(xiàn)N 為正常心音分類錯誤數(shù)，F(xiàn)P 為異常樣本分類錯誤數(shù)。 N 為樣本總數(shù)，Osamp為第samp 個樣本真實標記值，O?samp為第samp 個樣本預測值。

需要進一步說明的是，TP、TN、FP、FN 是通過心音記錄心動周期的分類，完成心音記錄的判別獲得的。其中設置了容錯機制，即心音記錄中，異常心動周期數(shù)與心動周期總數(shù)比例小于10%時，判定心音記錄為正常，否則判定為異常。

選用Challenge 2016 數(shù)據(jù)集[14]構建訓練集與測試集。整個數(shù)據(jù)集由a、b、c、d、e 這5 個子數(shù)據(jù)庫構成，分別采集于1 000多位不同年齡、不同性別、不同身體狀況的受試者，共計3 126 條心音記錄。按照官方給出的數(shù)據(jù)構建測試集，從5個子數(shù)據(jù)庫中選用301條心音，除去測試集，訓練集共計2 825條心音。

4 實驗結果與分析

按照預處理、周期檢測、PSD 時頻特征提取的流程獲取PSD時頻特征數(shù)據(jù)集。其中訓練集2 825條心音中檢測出84 762個心動周期，測試集301條心音中檢測出7 326個心動周期，最終構成98×40×84 762的訓練集，98×40×7 326的測試集，均按列存儲，生成.mat的數(shù)據(jù)文件。

圖8 訓練結果

分類器在訓練集作用下經(jīng)過100次迭代訓練，其結果如圖8所示，訓練集與測試集在分類器中獲得的結果如表2 所示。圖8 訓練結果顯示，迭代訓練60 次左右ACC 和LOSS 值均基本收斂，收斂后ACC 值為88.70%左右，LOSS 值為1.084 4左右。訓練完成后，測試集在網(wǎng)絡中的ACC 為84.72%，LOSS 值為1.185 7。

表2 訓練集與測試集結果

此外，靈敏度（Sensitivity，Se）和特異性（Specificity，Sp）[20，23]被廣泛用于評價分類系統(tǒng)的性能。靈敏度和特異性計算方法如式（21）和式（22）所示：

其中，WA1、WA2、WN1、WN2為權系數(shù)，分別為官方標記的正常情況下可分析性好壞的心動周期數(shù)量比例以及異常情況下可分析性好壞的心動周期數(shù)量比例；N表示標記為正常的心音，A 表示標記為異常的心音，n表示分類輸出為正常的心音，q 表示分類輸出不確定的心音，a 表示分類輸出為異常的心音，下標1 表示心音可分析性良好，2表示心音可分析性較差。權系數(shù)取值如表3所示。

表3 心動周期比例權系數(shù)

本文算法在訓練集和測試集中獲得的靈敏度與特異性結果如表4所示。其中訓練集靈敏度為0.834 8，特異性為0.959 8，整體得分為0.897 2；測試集靈敏度為0.776 3，特異性達到0.946 3，整體得分為0.861 3。

表4 靈敏度與特異性

本文算法同自相關&擴散圖特征&SVΜ[12]、分割后多維特征&非線性徑向基SVΜ[8]、光譜圖&SVΜ[9]、分割提取324 維特征&BPANN[22]、熱譜圖&CNN[19]、AdaBoost&CNN[23]（官方網(wǎng)站公布數(shù)據(jù)）算法進行了對比。文獻[12]采用小波分解提取四層分解中第四細節(jié)系數(shù)層和第二細節(jié)系數(shù)層的系數(shù)，通過香農(nóng)包絡自相關特征，并通過擴散圖映射到對應空間，通過幾何空間中特征的歐氏距離采用SVΜ 進行分類，但小波分解層次選擇中存在有用信息丟失的風險；文獻[8]與文獻[22]方法相似，均是通過被動提取多維特征后采用分類器進行分類，這往往又與特征的好壞有關；文獻[9]將傅里葉變換后的頻譜特征通過偏最小二乘法進行降維，然后通過非線性SVΜ進行分類，但非線性SVΜ并不適用于偏最小二乘降維后的數(shù)據(jù)分類[24]；文獻[19]將Μel 頻率倒譜系數(shù)（Μel Frequency Cepstral Coefficient，ΜF(xiàn)CC）特征轉換為熱譜圖，采用SVΜ對其進行分類，然而ΜF(xiàn)CC是基于人耳的聽覺特性，用于心音信號特征的表征能力往往有限[25]；文獻[23]提取了心動周期4個成分25～400 Hz中的9個頻帶和13維的ΜF(xiàn)CC特征，按照25～45 Hz、45～80 Hz、80～200 Hz、200～400 Hz集合特征訓練CNN分類器，并采用AdaBoost 算法將4 個CNN 組合構成強分類器，但分類的結果仍取決于提取的特征，且沒有計算4個頻帶間的過渡信息；本文算法提取心動周期的功率譜密度特征，構成二維時頻特征矩陣，采用CNN分類器主動提取特征，降低了被動提取特征的風險，使得特征更適應分類器，獲得更強的分類性能，事實上這在許多圖像分類領域中已經(jīng)得到了體現(xiàn)。

上述實驗的SVΜ 模型中，除了文獻[8]選用了非線性徑向基核函數(shù)外，其余均選擇Sigmoid核函數(shù)，參數(shù)σ和懲罰因子c 的取值為（2，2），詳細的模型描述請參考具體文獻。

對比實驗結果如表5所示。實驗結果表明，在SVΜ分類器下，PSD時頻特征最終獲得了0.830 0的得分，較自相關與擴散圖特征（0.754 9）、光譜圖降維特征（0.784 1）、分割提取特征（0.782 6[8]、0.828 6[22]）的方法得分更高，這也進一步說明PSD 時頻特征矩陣能夠用于心音分類。此外，與分割提取324 維特征&BPANN、熱譜圖&CNN、AdaBoost&CNN 等算法相比，本文算法獲得了較好的Sp 得分（0.946 3），Se 得分（0.776 3）相對較低，但與熱譜圖&CNN算法相比提升了0.048 5，且獲得了較好的總體評分0.861 3。從表5 中的對比結果可以看出，相比SVΜ 分類器而言，PSD 時頻特征在CNN 分類器下具備更好的分類效果。從算法整體評價結果來看，在測試集中取得了不錯的效果，這也進一步說明本文提出的分類算法具備可行性與有效性，即所構建的CNN 網(wǎng)絡能夠有效區(qū)分正常心音與異常心音的PSD特征矩陣。

表5 不同算法下測試集評價指標的對比

5 總結與討論

本文提出了一種基于功率譜時頻特征和CNN的正常與異常心音分類算法，無需進行心音分割，且具備較好的分類效果。處理中，首先進行小波降噪預處理，提高數(shù)據(jù)集的可分析性；然后，通過自相關法檢測出有效的心動周期，并采用Welch法提取功率譜特征；接著，采用雙線性插值法填補各心動周期，保持幀數(shù)一致。針對于特征矩陣數(shù)據(jù)集，構建了CNN分類器，并采用Challenge 2016數(shù)據(jù)集獲取訓練集和測試集，并進行訓練和測試。實驗結果表明，本文算法在評價體系中分類正確率和特異性等指標較好，但在靈敏度上需進一步提升。同時，本文算法無需進行心音分割，避免了分割帶來的誤差和增加算法復雜度，且本文算法可為現(xiàn)代心音信號分析與處理提供技術性參考。此外，建立在準確心血管疾病心音記錄數(shù)據(jù)集的基礎上，本文算法可進行推廣，應用于具體異常類型的心血管疾病篩選與檢測。