沈 伊,孫 靜,楊宏波,王威廉*

(1. 云南大學(xué)信息學(xué)院,云南 昆明 650500;2. 云南省阜外心血管病醫(yī)院,云南 昆明 650102)

1 引言

心音是人體重要的生理信號(hào)之一,其包含的生理和病理信息是臨床心臟疾病診斷重要依據(jù)。例如:對先天性心臟病(congenialheartdisease,CHD)的初診和篩查主要依靠心臟聽診。CHD是造成新生兒和嬰幼兒死亡的主要原因,在全國多地均位居新生兒出生缺陷的首位[1]。心臟聽診,對醫(yī)生的經(jīng)驗(yàn)和主觀性有極大的依賴,具有局限性、不穩(wěn)定性和不準(zhǔn)確性[2]。因此,用現(xiàn)代信息技術(shù)分析研究心音和輔助診斷技術(shù)具有重要意義。

近年來,國內(nèi)外學(xué)者對基于計(jì)算機(jī)的心音信號(hào)分析已有了顯著成果。主流分析流程包括以下三個(gè)步驟:預(yù)處理、特征提取和分類識(shí)別。其中,心音信號(hào)的預(yù)處理是心音研究的基礎(chǔ)和前提,主要包含去噪、包絡(luò)提取和分段定位三個(gè)部分。典型的心音分割研究可歸納為以下幾類:參考心電的分割[3,4]、基于時(shí)域特征的分割[5-7]、基于神經(jīng)網(wǎng)絡(luò)的分割[8]、基于概率模型的分割[9,10]。

Giordano N等[3]和王楨楨[4]通過心電的R波、T波與心音存在的對應(yīng)關(guān)系對心音進(jìn)行分割,前者對24例正常心音分割準(zhǔn)確率達(dá)到99.2%,后者對正常和異常心音進(jìn)行分割,平均準(zhǔn)確率達(dá)到98.5%。Mecheri等[5]通過計(jì)算各心動(dòng)周期之間的相關(guān)系數(shù)矩陣進(jìn)行相似性識(shí)別對心音進(jìn)行分割;Chen等[6]利用改進(jìn)型希爾伯特包絡(luò)并通過雙門限法對心音信號(hào)進(jìn)行分割;Zhao Zhan等[7]預(yù)先估計(jì)出心動(dòng)周期,提出一種基于心動(dòng)周期估計(jì)的心音分割,可以正確識(shí)別出80%以上的心動(dòng)周期。Yao Chen等[8]提出了一個(gè)持續(xù)時(shí)間長短記憶網(wǎng)絡(luò)結(jié)合持續(xù)時(shí)間特征的方法對心音進(jìn)行分割,準(zhǔn)確率達(dá)到96.36%。侯雷靜等[9]提出一種面向心音分割的個(gè)性化高斯混合建模方法(PMM)對480個(gè)心音信號(hào)進(jìn)行周期檢測與狀態(tài)分割,準(zhǔn)確率達(dá)到93.4%;Oliveira[10]使用邏輯回歸函數(shù)作為先驗(yàn)概率建立HMM模型進(jìn)行心音狀態(tài)分割,平均準(zhǔn)確率達(dá)93%。

上述研究中,參考心電的心音分割算法是利用心電圖(electrocardiogram,ECG)與心音圖(phonocardiogram,PCG)的對照關(guān)系對心音進(jìn)行分割。該方法需要對心音、心電進(jìn)行同步采集和同步處理,采集儀器和采集過程較為復(fù)雜,尤其將采集儀器放置在新生兒或嬰兒身上時(shí)比較困難。此外,在左心室或右心室肥大的兒童中,心臟軸左偏或右偏會(huì)導(dǎo)致心電圖信號(hào)異常,使心音分割復(fù)雜化[11]。基于神經(jīng)網(wǎng)絡(luò)的心音分割算法主要是將已標(biāo)記的心音用神經(jīng)網(wǎng)絡(luò)進(jìn)行訓(xùn)練,從而得到心音的狀態(tài)序列。該方法分割準(zhǔn)確率高,但通常需要大量數(shù)據(jù)、復(fù)雜的網(wǎng)絡(luò)結(jié)構(gòu)以及計(jì)算機(jī)強(qiáng)大的計(jì)算能力,且不利于算法部署到移動(dòng)設(shè)備。基于概率模型的心音分割算法使用心音特征參數(shù)作為模型的觀測序列,然后對心音狀態(tài)進(jìn)行分割。該方法性能很大程度受限于轉(zhuǎn)移概率矩陣參數(shù),在處理不同樣本時(shí),分割性能會(huì)發(fā)生變化,沒有較好的普適性。

針對上述問題,考慮到基于時(shí)域特征的分割能保留心音時(shí)間序列中的峰值信息,具有運(yùn)算速度快、普適性較好以及算法部署至移動(dòng)設(shè)備較為簡單的優(yōu)勢[12],本文提出一種包絡(luò)結(jié)合突變點(diǎn)檢測與峰值點(diǎn)搜索的心音分割算法,目的是在保證準(zhǔn)確分割的前提下能快速分割心音信號(hào),且具有較強(qiáng)的魯棒性,能適應(yīng)不同聽診部位導(dǎo)致的數(shù)據(jù)差異性問題,便于對心音信號(hào)進(jìn)行特征提取,提高分類識(shí)別的準(zhǔn)確率。

2 心音組成及算法框架

2.1 心音的組成

心音包括第一心音(S1)、第二心音(S2)、第三心音(S3)和第四心音(S4)四個(gè)成分。其中,S3、S4一般不是研究的重點(diǎn)[13]。一個(gè)完整的心動(dòng)周期由S1、心臟舒張期、S2和心臟收縮期構(gòu)成。通常情況下,心臟舒張期要長于心臟收縮期。心音與心電的對應(yīng)關(guān)系如圖1所示,ECG信號(hào)R峰的位置對應(yīng)PCG信號(hào)中S1的起始,S2的位置出現(xiàn)在T波末端附近[14]。

圖1 心音與心電對應(yīng)關(guān)系圖

2.2 算法框架

本文提出的心音分割算法系統(tǒng)框架如圖2所示。主要包括心音預(yù)處理、突變點(diǎn)檢測、峰值搜索和心音分割四個(gè)部分。其中,突變點(diǎn)指的是S1、S2與收縮期、舒張期的分界點(diǎn)。

圖2 心音分割算法系統(tǒng)框架

算法具體步驟如下:第一步,對心音信號(hào)進(jìn)行預(yù)處理。主要包括以下幾個(gè)過程:將心音信號(hào)采樣頻率統(tǒng)一至1000Hz,之后對心音信號(hào)進(jìn)行去噪和幅值歸一化處理。第二步,定位突變點(diǎn)。對心音短時(shí)能量與頻譜質(zhì)心進(jìn)行閾值判斷,區(qū)分出心音信號(hào)振動(dòng)部分與非振動(dòng)部分,其分界點(diǎn)即為心音信號(hào)突變點(diǎn)。第三步,峰值搜索。根據(jù)突變點(diǎn)確定大小區(qū)間,用其代替?zhèn)鹘y(tǒng)初始閾值對心音包絡(luò)極大值進(jìn)行前、后向搜索,在濾除偽峰的同時(shí)通過峰值點(diǎn)補(bǔ)償算法提高有效峰值的檢出率。第四步,心音分割。利用峰峰值之間的距離定位S1和S2,并以突變點(diǎn)為分割點(diǎn),將心音信號(hào)分割成逐個(gè)完整的心動(dòng)周期。

3 方法

3.1 預(yù)處理

心音屬于人體微弱信號(hào),由于人體系統(tǒng)和檢測環(huán)境的復(fù)雜性,采集心音時(shí)往往會(huì)伴隨各種環(huán)境噪聲或電子信號(hào)的干擾。去噪的目的是濾除心音數(shù)據(jù)中的雜音,保留有效和關(guān)鍵的信息。近年來,正交小波分解在主頻率范圍與噪聲頻率沒有交集的非平穩(wěn)隨機(jī)信號(hào)的去噪效果上表現(xiàn)出色。原因是正交小波分解算法可將信號(hào)分解為多個(gè)頻段,并根據(jù)實(shí)際需要對各頻段成分進(jìn)行交換和重構(gòu),將有效信號(hào)從環(huán)境噪聲中分離并提升信號(hào)有效成分的信噪比[15]。雙正交小波是正交小波概念的推廣。由于正交小波缺乏對稱性,相應(yīng)的濾波器不能保持線性相位,而雙正交小波可解決正交小波線性相位和正交性要求的矛盾[16]。

首先將心音的采樣頻率統(tǒng)一至1000Hz,降低數(shù)據(jù)大小,提高運(yùn)算速度。然后對心音信號(hào)進(jìn)行去噪,選擇bior5.5小波將心音信號(hào)分為五層,選取適合該段心音的閾值函數(shù)得到小波系數(shù),濾除噪聲后對信號(hào)進(jìn)行重構(gòu),得到去噪后的心音信號(hào)。原始信號(hào)去噪時(shí)域波形如圖3所示。圖3(a)為未經(jīng)處理的原始心音,圖中A、B為收縮期、舒張期雜音。

圖3 小波去噪對比圖

由圖3(b)可以看出,原始心音信號(hào)使用db6小波進(jìn)行五層分解后仍有噪聲與信號(hào)重疊部分,而采用雙正交bior5.5小波進(jìn)行五層分解后,心音波形更加清晰,噪聲基本濾除,如圖3(c)所示。

3.2 突變點(diǎn)檢測

根據(jù)心音產(chǎn)生機(jī)制,S1是由心室收縮開始時(shí)二尖瓣和三尖瓣關(guān)閉,瓣葉突然緊張引起振動(dòng)而產(chǎn)生;S2是由心室舒張開始時(shí)主動(dòng)脈瓣和肺動(dòng)脈瓣突然關(guān)閉引起的瓣膜振動(dòng)所產(chǎn)生[17]。因此相較于舒張期和收縮期,S1和S2的振動(dòng)更加明顯,平均能量更高。在語音信號(hào)處理中,短時(shí)能量(Short-term energy)可以反映信號(hào)的時(shí)域特征,頻譜質(zhì)心(Spectral Centroid)可以用來計(jì)算頻譜的位置[18]。結(jié)合兩者對心音信號(hào)進(jìn)行分析,可有效的區(qū)分出心音信號(hào)中的振動(dòng)部分與非振動(dòng)部分,尋找出心音信號(hào)的突變點(diǎn)。突變點(diǎn)檢測流程框圖如圖4所示。

圖4 突變點(diǎn)檢測流程

突變點(diǎn)檢測算法具體步驟如下:首先,對歸一化的去噪心音信號(hào)進(jìn)行分幀處理。計(jì)算每幀信號(hào)的短時(shí)能量和頻譜質(zhì)心。其中,短時(shí)能量表示一幀信號(hào)樣點(diǎn)值的平方和,頻譜質(zhì)心是基于信號(hào)能量分布的頻率一階矩,公式如式(1)和式(2)所示。設(shè)M1和M2分別為短時(shí)能量和頻譜質(zhì)心直方圖第一大、第二大局部最若某一段信號(hào)的短時(shí)能量與頻譜質(zhì)心同時(shí)小于其閾值,則認(rèn)為該段信號(hào)沒有發(fā)生振動(dòng)。反之,若某一段信號(hào)的短時(shí)能量與頻譜質(zhì)心同時(shí)大于其閾值,則認(rèn)為該段信號(hào)發(fā)生了振動(dòng)。并將同時(shí)大于短時(shí)能量與頻譜質(zhì)心閾值的信號(hào)段(也即兩者的交集部分)的左右兩個(gè)端點(diǎn)定義為該段心音。信號(hào)的突變點(diǎn)。突變點(diǎn)位置如圖5(c)中“*”號(hào)所示。

圖5 心音信號(hào)短時(shí)能量、頻譜質(zhì)心及信號(hào)突變點(diǎn)

(1)

式中,E表短時(shí)能量,S為幀信號(hào),i為第i幀。

(2)

式中,C表示頻譜質(zhì)心,Xi(k),k=1,…,N是第i幀信號(hào)的離散傅里葉變換,N為幀長度。

(3)

其中,W為手動(dòng)設(shè)置的閾值,本文取W=5。

3.3 峰值搜索

目前峰值搜索算法有:簡單比較尋峰法、導(dǎo)數(shù)尋峰法、協(xié)方差尋峰法等。其中,簡單比較尋峰法計(jì)算簡單,適合查找孤立的峰值,但搜索時(shí)容易受到連續(xù)背景信息和噪聲的影響;導(dǎo)數(shù)尋峰法適合簡單的峰形,但檢出率不高且適應(yīng)能力較差;協(xié)方差尋峰法主要優(yōu)點(diǎn)是分辨重峰的能力強(qiáng),適合尋找弱峰,但計(jì)算過程復(fù)雜、運(yùn)算速度較慢。由于去噪過程無法保證將信號(hào)完全從環(huán)境雜音中剝離,容易導(dǎo)致偽峰的產(chǎn)生,并且上述算法在包絡(luò)的尋峰過程中均無法保證偽峰的全部濾除,不利于后續(xù)對峰值的定位與心音成分的識(shí)別。因此,本文提出一種導(dǎo)數(shù)尋峰法與閾值搜索結(jié)合的尋峰算法,并采用峰值補(bǔ)償增加有效峰值的檢出,最大程度地減少偽峰對尋峰過程的干擾。具體步驟如下:

1)首先對去噪信號(hào)進(jìn)行維奧拉積分包絡(luò)的提取,公式如式(4)所示,然后對包絡(luò)進(jìn)行歸一化處理,公式如式(5)所示。雖然維奧拉積分包絡(luò)具有簡單、處理速度快和實(shí)時(shí)性強(qiáng)的優(yōu)勢[19],但包絡(luò)存在毛刺尖峰,因此用三階Savitzky-Golay(S-G)濾波器對窗口內(nèi)的信號(hào)包絡(luò)進(jìn)行三階多項(xiàng)式擬合,得到平滑后的維奧拉積分包絡(luò)。

(4)

(5)

式中,Et為心音信號(hào)的維奧拉積分包絡(luò),E為平滑后的維奧拉積分包絡(luò)。

2)根據(jù)突變點(diǎn)確定峰值搜索的大小區(qū)間。公式如式(6)和式(7)所示。

(6)

(7)

式中,TH1為小區(qū)間,TH2為大區(qū)間。

3)通過二階導(dǎo)數(shù)法求取包絡(luò)所有極大值并記錄局部極大值索引表R1。利用大小區(qū)間對R1進(jìn)行正、反向搜索,并計(jì)算正反搜索結(jié)果的峰值點(diǎn)距離平均均方差,將平均均方差較小的一個(gè)搜索結(jié)果記為R2。根據(jù)式計(jì)算閾值TH3,并將其作為小區(qū)間搜索結(jié)果的判斷閾值,篩選出第三次峰值索引結(jié)果集R3。

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

(8)

式中,Env(R2)為結(jié)果集R2對應(yīng)的幅值,Median(Env(R2))表示取結(jié)果集R2對應(yīng)的幅值的中位數(shù)。

為了減少有效信息的丟失,提高S1、S2定位識(shí)別的準(zhǔn)確率,增加偽峰消除的步驟。利用S1、S2的最大持續(xù)時(shí)間結(jié)合其幅值高于收縮期與舒張期幅值的特點(diǎn)對R3結(jié)果集進(jìn)行偽峰消除,增加檢出正確率,得到最終定位結(jié)果集Rend。峰值點(diǎn)補(bǔ)償算法如圖6所示。

圖6 峰值點(diǎn)補(bǔ)償算法

圖7為峰值搜索及補(bǔ)償算法效果圖。圖7(a)中由于心音收縮期雜音較強(qiáng)產(chǎn)生了偽峰,經(jīng)過峰值定位及補(bǔ)償算法后,增加了有效峰值的檢出,如圖7(b)所示。

圖7 偽峰去除結(jié)果圖

PCG信號(hào)中S1和S2的正確識(shí)別定位難度較大的原因有很多:在時(shí)域中,其持續(xù)時(shí)間和形態(tài)相似;在頻域中,其頻譜成分互相重疊(S1從10到140Hz,S2從10到200 Hz)[20]。此外,噪聲與正常心音和異常心音在形態(tài)學(xué)上的相似性使得在時(shí)域上識(shí)別它們非常困難。通常情況下,心音的舒張期比收縮期長,可以以此作為判定S1、S2的依據(jù)。根據(jù)突變點(diǎn)作為分割點(diǎn),可以分割出心音的S1、收縮期、S2與舒張期。圖8可從心電與心音內(nèi)在關(guān)系的角度客觀評價(jià)分割正常心音、異常心音的準(zhǔn)確性。虛線表示心音的四個(gè)狀態(tài),其中,狀態(tài)“1”表示第一心音(S1)、狀態(tài)“2”表示收縮期、狀態(tài)“3”表示第二心音(S2)、狀態(tài)“4”表示舒張期。從圖8可以看出分割的四個(gè)狀態(tài)與心電圖之間略有偏差,但由于該偏差在誤差范圍之內(nèi),可忽略不計(jì)。

圖8 基于本文算法的正常、異常心音分割結(jié)果

4 實(shí)驗(yàn)結(jié)果與分析

4.1 實(shí)驗(yàn)環(huán)境與數(shù)據(jù)說明

實(shí)驗(yàn)所用電腦CPU型號(hào)為AMD Ryzen7 4800HS,RAM為16.0GB,軟件環(huán)境為MATLAB2019b。本文實(shí)驗(yàn)數(shù)據(jù)從PhysioNet/Computing in Cardioloy Challenge選取了1000例心音信號(hào),并通過同步采集的ECG信號(hào)對S1、S2的位置進(jìn)行了標(biāo)記,標(biāo)記結(jié)果通過多種R峰與T波探測器綜合給出[21]。所選取的心音中包含539例正常心音,461例異常心音。

4.2 評價(jià)指標(biāo)

為評估該心音分割算法性能以及其在測試集中準(zhǔn)確定位S1、S2的能力,本文采用ECG-PCG對照法對其進(jìn)行評估。由于心音分割不是分類問題,因此不存在真正的負(fù)樣本,因此:

1)若算法定位出的S1起始位置在ECG R峰的100ms內(nèi),則將其認(rèn)為分段準(zhǔn)確(TP)。

2)若算法定位出的S2中心位置在ECG T波末端的100ms內(nèi),則將其認(rèn)為分段準(zhǔn)確(TP)。

3)其余情況均認(rèn)為是錯(cuò)誤分段(FP)。

使用Se：靈敏度,ACC：準(zhǔn)確率,F1：評估精度分?jǐn)?shù)(綜合評價(jià)指標(biāo)),SN：每分鐘內(nèi)分割心音的個(gè)數(shù),四種評估標(biāo)準(zhǔn),分別由式(9)-(11)定義。

(9)

(10)

(11)

4.3 實(shí)驗(yàn)結(jié)果

對正常心音、異常心音分割的性能比較如表1所示。本文算法分割性能與引言中提到的其中幾種分割算法進(jìn)行對比的結(jié)果如表2所示。對比試驗(yàn)計(jì)算SN指標(biāo)時(shí),統(tǒng)一采用4.5s的心音段,計(jì)算一分鐘內(nèi)可分割心音的總個(gè)數(shù)。

表1 正常、異常心音分割的性能比較

表2 5種分割方法的性能比較

4.4 實(shí)驗(yàn)結(jié)果分析

通過表1中的準(zhǔn)確率評價(jià)指標(biāo)可看出,本文算法對于正常心音的分割準(zhǔn)確率較高于異常心音。通過分析,這是由于異常心音中病理性特征的存在,往往使異常心音更為復(fù)雜,導(dǎo)致對異常心音的突變點(diǎn)檢測發(fā)生錯(cuò)誤、偽峰去除不完整或S1、S2難以定位。

通過表2可看出,本文算法在準(zhǔn)確度、精度和F1指標(biāo)要優(yōu)于方法1-3。實(shí)驗(yàn)表明,所提算法通過心音振動(dòng)段與非振動(dòng)段進(jìn)行突變點(diǎn)檢測并根據(jù)突變點(diǎn)計(jì)算自適應(yīng)閾值,對心音更具有普適性,并用峰值搜索對偽峰進(jìn)行了濾除并進(jìn)行了峰值補(bǔ)償操作,增加了有效峰值的檢出率,使分割準(zhǔn)確率有了提升。通過SN指標(biāo)可看出,雖然方法4具有較高的準(zhǔn)確率和F1,但由于其包含高斯混合建模、前后向運(yùn)算和修正誤差等計(jì)算過程,導(dǎo)致該分割法相較于其它分割算法計(jì)算過程復(fù)雜,耗時(shí)較長。由于方法2需要先進(jìn)行經(jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decomposition,EMD)找到固有模態(tài)函數(shù)(intrinsicmodefunction,IMF)并進(jìn)行插值運(yùn)算,因此其運(yùn)算復(fù)雜度僅次于方法4。本文算法在尋找到心音突變點(diǎn)和心音有效峰值后可根據(jù)心音時(shí)域特征直接進(jìn)行分割,復(fù)雜度遠(yuǎn)遠(yuǎn)低于以上兩種算法。且相較于方法4,本文算法更容易部署到移動(dòng)端,更有利于先心病機(jī)器輔助診斷系統(tǒng)的普及。

5 結(jié)束語

基于時(shí)域特征的心音分割方法算法具有計(jì)算過程簡單、運(yùn)算速度快的特點(diǎn)。其中突變點(diǎn)檢測與峰值定位的過程為關(guān)鍵步驟。本文提出一種利用短時(shí)能量和頻譜質(zhì)心結(jié)合的突變點(diǎn)檢測算法和基于維奧拉積分包絡(luò)的自適應(yīng)閾值尋峰算法,根據(jù)心音信號(hào)時(shí)域特點(diǎn)進(jìn)行峰值篩查與峰值補(bǔ)償,識(shí)別并去除偽峰,從而得到正確的S1、S2定位識(shí)別,實(shí)現(xiàn)了心音信號(hào)的快速分割。相比個(gè)性化高斯混合建模的包絡(luò)分割法,本文算法在保證準(zhǔn)確率的情況下更加快速且具有可移植性。相較于傳統(tǒng)閾值選取分段算法,本文算法解決了其適用性低、自適應(yīng)弱和依賴調(diào)節(jié)系數(shù)等問題,采用大小區(qū)間搜索與峰值補(bǔ)償算法得到有效的心音峰值點(diǎn),保證了分段結(jié)果的準(zhǔn)確性。但本文算法對于環(huán)境噪聲過大且無法將心音從噪聲剝離的信號(hào)分割效果尚不理想,有待進(jìn)一步研究,以達(dá)到更好的分割效果。