基于DHMM的低心率變異性心音的分割方法

許春冬 周 靜 應(yīng)冬文, 侯雷靜 龍清華

（1.江西理工大學(xué)信息工程學(xué)院，贛州，341000；2.中國科學(xué)院聲學(xué)研究所語言聲學(xué)與內(nèi)容理解重點實驗室，北京，100190）

引 言

隨著社會生活的不斷復(fù)雜化，心血管疾病的患病率逐年上升，已經(jīng)成為全球慢性病的主要致死原因之一[1]。在心血管類疾病的診斷中，常用的方法包括血壓檢測[1]、心電分析[2]、心臟超聲檢查[1-2]、冠狀動脈CT[1]、心肌灌注同位素圖像法[3]和心音分析[1-3]等。其中心音分析不僅可以檢測出受試者是否存在異常，還可通過其作出病況預(yù)測，且方便、快速、經(jīng)濟(jì)，是一種完全無損的檢測方法。心音分析已經(jīng)得到了廣泛應(yīng)用，并成為心血管疾病主要的輔助診斷方法之一。

心音信號由心臟機(jī)械工作產(chǎn)生，其蘊(yùn)含著豐富的生理與病理信息。心音分割是心音自動分析的主要步驟之一。心音周期性分析及分割的根本目的是通過提取各心音成分上的病理特征，使后續(xù)的心音分類及診斷目標(biāo)性更強(qiáng)、準(zhǔn)確度更高。現(xiàn)實中，心音極易受到臨床環(huán)境噪音、體內(nèi)雜音、信道差異、以及心率波動等因素的影響，這使得分割過程變得十分困難。

現(xiàn)階段主要的分割算法包括基于包絡(luò)的分割方法[4]、基于特征提取分析的分割方法[5]、基于機(jī)器學(xué)習(xí)的分割方法[6]、基于隱馬爾可夫模型（Hidden Markov model,HMM）的分割方法[7-9]。其中基于HMM的方法能夠通過觀察值序列和隱含狀態(tài)序列間的關(guān)系反映出的隨機(jī)過程判斷出存在的合理狀態(tài)序列，這種特點使HMM在心音分割處理中得到了一定深度的研究。HMM在語音活動性檢測中已經(jīng)得到了深入的研究，但在心音定位分割及分類領(lǐng)域，從模型建立和分割精度要求角度來看，還存在一定的研究空間?，F(xiàn)階段基于HMM的分割算法多是通過大量訓(xùn)練集進(jìn)行的非監(jiān)督學(xué)習(xí)，容易造成模型泛化而無法保證提取的信號特征的與信號的高度匹配，會丟失部分細(xì)節(jié)特征，從而使得對目標(biāo)信號的分割精度下降。因此，結(jié)合HMM在心音分割中的研究現(xiàn)狀，本文采用基于時間相關(guān)性的隱馬爾可夫模型（Duration-dependent hidden Markov model,DHMM）探求精度更高的分割算法。

1 心音成分及分割任務(wù)要求

心音一般是指由心肌收縮、心臟瓣膜關(guān)閉、血液撞擊心室壁和大動脈壁等活動引起的振動經(jīng)胸腔、胸壁等組織傳導(dǎo)至體表所產(chǎn)生的聲學(xué)信號[10-11]，是能夠反映人體心臟及心血管系統(tǒng)機(jī)械運(yùn)動狀況的信號之一[12]。心臟結(jié)構(gòu)如圖1所示，心臟通過循環(huán)的收縮活動來為人體不斷輸送血液，其活動狀態(tài)為一周期過程。對于正常人而言，心音信號近似于準(zhǔn)周期信號，每個周期約在600 ms～1 s之間[2]，醫(yī)學(xué)上稱作心動周期。每個心動周期內(nèi)可包含4個心音成分[1,3,7]：第1心音（s1）、第2心音（s2）、第3心音（s3）和第4心音（s4）。一般s1和s2極易被聽到，而s3和s4主要出現(xiàn)在特定人群的心音中，如s3主要出現(xiàn)在兒童心音信號中，s4主要出現(xiàn)在一些病理心音信號中。在醫(yī)學(xué)上，常認(rèn)為一個心動周期由s1，收縮期，s2和舒張期4個部分構(gòu)成，且在任意一個心動周期中均存在，出現(xiàn)順序是固定循環(huán)的，如圖2所示。

圖1 心臟結(jié)構(gòu)Fig.1 Structure of the heart

醫(yī)學(xué)上將s1和s2稱作基礎(chǔ)心音（Foundamental heart sound,FHS）[2,5,8]。心音分割的主要任務(wù)是精準(zhǔn)有效地將各心動周期劃分為s1，收縮期，s2及舒張期，以達(dá)到有效提取各心音成分特征參數(shù)的目的。此外，心音分割可為后續(xù)的模式識別（分類）工作提供定位基準(zhǔn)，是實現(xiàn)瓣膜病、心率衰竭和冠心病等疾病無創(chuàng)、快捷診斷的前提與基礎(chǔ)。心動周期分割的有效與否，將直接影響到整個心音分析系統(tǒng)的精準(zhǔn)性，因此，在心音自動分析處理中，分割步驟極為重要。

圖2 心音成分Fig.2 Composition of heart sound

2 基于DHMM的心音分割

在心音信號的分割中，一些基于HMM的方法取得了顯著的效果[7-9]?；贖MM的方法一般是通過大量人工標(biāo)注的訓(xùn)練集提取特征參數(shù)，然后根據(jù)分割狀態(tài)建立相應(yīng)的模型，從而完成分割任務(wù)。

2.1 心音信號的標(biāo)注

心音分割方法中，使用心電圖與心音信號同步采集的方法，為心音分割提供了有效的參考，是目前精度最高的分割方法之一[2,12]。但在臨床診斷中，該方法由于操作復(fù)雜，增加了醫(yī)護(hù)人員的操作難度，臨床應(yīng)用可行性較差。同時，為了便于心音分析在常規(guī)體檢項目中應(yīng)用，對實時性的要求也較高。故該方法只在對精度有較高要求的情況下得到了一定的應(yīng)用。

部分研究者提出通過參考心電信號標(biāo)注訓(xùn)練集訓(xùn)練模型的方法來提高分割精度。然而，該方法在擴(kuò)大更新標(biāo)注訓(xùn)練集時耗時耗力，顯然不可行。只通過采集心音信號標(biāo)注的方法，雖然精度略低于參考心電信號的方法，但其臨床實用性強(qiáng)，操作簡單，得到了廣大研究者和醫(yī)護(hù)人員的認(rèn)可，成為現(xiàn)階段主要的研究方向。

2.2 面向心音信號的HMM與DHMM

HMM是一種關(guān)于時間序列的統(tǒng)計模型，它在滿足馬爾可夫性的基礎(chǔ)上，描述了隱含的隨機(jī)序列轉(zhuǎn)移及生成對應(yīng)的可觀測隨機(jī)序列的過程，是一類雙離散隨機(jī)序列[13]模型。一個簡單的HMM過程的描述如圖3所示。隱含狀態(tài)隨機(jī)序列轉(zhuǎn)移過程描述的是有限個可知狀態(tài)的隨機(jī)轉(zhuǎn)移過程，其轉(zhuǎn)移狀態(tài)序列是不可知的；可觀測隨機(jī)序列描述的是隱含狀態(tài)與可觀測序列間的一種關(guān)聯(lián)統(tǒng)計關(guān)系，可觀測序列是可知的[14]。HMM隱含序列轉(zhuǎn)移到下一時刻的狀態(tài)僅與隱含序列當(dāng)前時刻的狀態(tài)相關(guān)，而與觀測隨機(jī)序列和隱含序列過去時刻的狀態(tài)無關(guān)。因此，HMM符合心音序列的統(tǒng)計模型，能夠較好地描述心音序列的短時平穩(wěn)性及整體的非平穩(wěn)性。

圖3HMM過程Fig.3 The process of HMM

一個基本的HMM中的關(guān)鍵參數(shù)主要包括有限數(shù)目的隱含狀態(tài)集合S及隱含狀態(tài)序列St，各隱含狀態(tài)所有的可觀察值集合VOS及與隱含序列對應(yīng)的觀察序列Vt，各隱含狀態(tài)的初始概率分布πc，各隱含狀態(tài)間的轉(zhuǎn)移概率矩陣δ，以及各隱含狀態(tài)可能的觀察值對應(yīng)的概率分布矩陣Ψ。其中πc和δ決定了隱含狀態(tài)序列，Ψ決定了觀察值序列，因此往往將一個基本的HMM模型概括為

在利用HMM分割心音信號時，其有限狀態(tài)集合可表示為

則其轉(zhuǎn)移矩陣可表示如表1，其中ST表示跳轉(zhuǎn)至本狀態(tài)的概率，SN表示跳轉(zhuǎn)至下一狀態(tài)的概率，STi+SNi=1，STi＞SNi,1 ≤i≤ 4。

由于HMM的狀態(tài)轉(zhuǎn)移存在兩種情況，即本狀態(tài)至本狀態(tài)或者本狀態(tài)至下一狀態(tài)，容易出現(xiàn)個別點誤判現(xiàn)象，且轉(zhuǎn)移狀態(tài)數(shù)目較多，概率匹配相對復(fù)雜。因此考慮將HMM模型進(jìn)行優(yōu)化，以尋求更好的性能?？紤]采用基于時間相關(guān)性的隱馬爾可夫模型(Duration-dependent hidden Markov model,DHMM)[9]來對分割過程進(jìn)行優(yōu)化。

表1 心音的HMM狀態(tài)轉(zhuǎn)移矩陣Tab.1 HMM state transition matrix of heart sounds

表2 心音的DHMM狀態(tài)轉(zhuǎn)移矩陣Tab.2 DHMM state transition matrix of heart sounds

DHMM與HMM的不同之處在于對狀態(tài)轉(zhuǎn)移加入了時間關(guān)系的條件限制，即經(jīng)過一段相關(guān)聯(lián)時間的累積，當(dāng)各狀態(tài)跳轉(zhuǎn)時間的聯(lián)合轉(zhuǎn)移概率滿足條件時才判斷為狀態(tài)轉(zhuǎn)移，不滿足則認(rèn)為仍是當(dāng)前狀態(tài)[9]。同時，僅將轉(zhuǎn)移狀態(tài)劃分為基礎(chǔ)心音和非基礎(chǔ)心音兩類，無需劃分為s1，收縮期，s2和舒張期4類，模型更加精簡化。由于心音信號總是按照s1，收縮期，s2和舒張期的成分循環(huán)出現(xiàn)，且舒張期持續(xù)時長一般大于收縮期持續(xù)時長，故能較容易地在區(qū)分基礎(chǔ)心音與非基礎(chǔ)心音的基礎(chǔ)上將4類成分區(qū)分開來。因此在DHMM下的有限狀態(tài)集合可表示為

其狀態(tài)轉(zhuǎn)移矩陣可表示為表2。

在DHMM中，前向概率δt(j)[15]可表示為

即i時刻隱含值Sti為j狀態(tài)sj，i+1時刻隱含值St(i+1)不為sj，即發(fā)生跳轉(zhuǎn)。則Sti為sj狀態(tài)最后時間的概率可表示為

式中：pj(d)表示狀態(tài)時長d下的分布概率密度函數(shù)[15-16]；δt-d(i)表示t-d為i狀態(tài)下的最后時刻的概率；aij表示從i狀態(tài)跳轉(zhuǎn)至j狀態(tài)的概率。則可將前向概率表示為

其中d和i的最大化值將用于后續(xù)的回溯算法，詳細(xì)可參考文獻(xiàn)[15,16]。

2.3 心動周期的提取

心音信號的各心動周期之間的差異較小，且醫(yī)學(xué)上往往是在機(jī)體活動穩(wěn)定的狀態(tài)下在進(jìn)行采集，以保證心率的變異性較低。在理想的環(huán)境下采集到的正常心音信號一般呈準(zhǔn)周期特性，即各心動周期的持續(xù)時長基本一致。

心音信號具有短時平穩(wěn)性，研究中普遍認(rèn)為在10～30 ms內(nèi)心音信號可視作平穩(wěn)[5]。各心動周期變異性較低的情況下，通過自相關(guān)法能夠較好地提取出心動周期的長度，自相關(guān)方法的特征包絡(luò)提取可表達(dá)為

式中：E(n)為待分割的心音信號的第n幀信號短時平均能量，r為信號的總長度，F(xiàn)L為幀長。求取的自相關(guān)包絡(luò)如圖4所示。在圖4中，每隔一段時間會出現(xiàn)一個峰值，每兩個峰值間對應(yīng)的時移恰好是一個心動周期持續(xù)時間的長度。在心動周期的提取中，可采用如下方法：

（1）采用中值濾波方法將包絡(luò)進(jìn)行平滑，去除每兩個峰值點間存在干擾的極大值點；

（2）提取相關(guān)函數(shù)包絡(luò)的極大值點；

（3）對提出的極大值點進(jìn)行降序排列，取出前兩個極值點的位置信息；

圖4 自相關(guān)提取心音周期Fig.4 Autocorrelation extraction heart sound cycle

2.4EEMD分解

根據(jù)2.3節(jié)中提取的心動周期T，將心音信號劃分為若干個心動周期，并對各心音周期進(jìn)行集合經(jīng)驗?zāi)B(tài)分解（Ensemble empirical mode decomposition,EEMD）[17],得到多個本征模態(tài)函數(shù)（Intrinsic mode function,IMF）分量，每個IMF分量中都可能包含有用心音信號。這里采用IMF分量與心音信號的相關(guān)系數(shù)差值來衡量心音成分含量的多少，剔除部分無用的IMF分量，在抑制噪聲的同時以盡量多的含心音信號的IMF分量重構(gòu)來表征心音信號，其可分析性更高。相關(guān)系數(shù)的定義如下

式中：N為分量總數(shù)，xk表示第k個心動周期為第k個心動周期的均值；fi表示第i個IMF分量表示各IMF分量的均值。相關(guān)系數(shù)差值可表示為

在相關(guān)系數(shù)的判決中，求得的互相關(guān)系數(shù)差值ΔRki的第一個局部最小值fk，則fk+1可視作模態(tài)混疊分量。這里將前k個IMF分量視為噪聲分量，其中心有效的音信號分量極少，因此可直接丟棄，將剩余分量進(jìn)行重構(gòu)即得到各心動周期的模態(tài)本征信號。

2.5 心音信號的高斯建模

基礎(chǔ)心音和非基礎(chǔ)心音存在某種約束關(guān)系[18-19]，本文分析了心動周期內(nèi)的短時平均能量分布及其直方分布，建立如圖5所示的高斯混合模型（Gaussian mixture model,GMM）。

圖5（a）為一心動周期內(nèi)短時平均能量分布圖，圖5（b）為噪聲干擾下的心音信號時域波形圖，與圖5(a)相對應(yīng)；圖5(c)為根據(jù)短時能量分布特征進(jìn)行的高斯擬合過程。由于非基礎(chǔ)心音信號持續(xù)時間較長，出現(xiàn)次數(shù)較多，且短時平均能量較小，集中在0附近，故擬合的高斯波形狹窄而高，即高斯模型方差較?。换A(chǔ)心音信號持續(xù)時間較短，總的出現(xiàn)次數(shù)較小，短時平均能量分布較為分散，故擬合的高斯波形寬廣而矮，即方差較大。4類常見噪聲的干擾均符合上述情況，因此可根據(jù)直方分布情況建立GMM模型，即式中：σn為標(biāo)準(zhǔn)差，ηn為期望值，n=1或2。當(dāng)n=1時，表示為非基礎(chǔ)心音；當(dāng)n=2時，表示為基礎(chǔ)心音。

圖5 心動周期的高斯建模Fig.5 Gaussian modeling of cardiac cycle

模型的約束關(guān)系可概括為[18]：

前瞻性研究。本研究數(shù)據(jù)分析采用SPSS 22.0軟件包，采用重復(fù)測量的方差分析、卡方檢驗進(jìn)行數(shù)據(jù)比較。以P＜0.05為差異有統(tǒng)計學(xué)意義。

（1） 基礎(chǔ)心音和非基礎(chǔ)心音狀態(tài)間的約束關(guān)系為式中：τ＞0，用于防止誤判和基礎(chǔ)心音泄漏，其估計為

式中ρ為常數(shù)。

（2） 方差間的約束關(guān)系為

非基礎(chǔ)心音段將基礎(chǔ)心音方差設(shè)置為σ1，通過約束條件（1，2）使得模型區(qū)分兩種狀態(tài)的能力得以保持。

通過約束關(guān)系建立合適的GMM模型，由DHMM模型訓(xùn)練，來完成基礎(chǔ)心音與非基礎(chǔ)心音的判別，此時模型可表示為

3 心音信號的分割流程

本算法是基于準(zhǔn)周期心音信號前提提出的，分割處理的心音信號心率變異性應(yīng)較小。整個心音分割流程可概括為以下6個步驟，流程圖如圖6所示。

（1）輸入待分割的心音信號并作預(yù)處理，通過自相關(guān)包絡(luò)獲取心動周期均值T；

（2）通過EEMD分解計算每個心動周期的IMF分量，計算相關(guān)系數(shù)R(xk,fi)并用其差值ΔRki第1個極小值后的分量重構(gòu)信號；

（3）通過約束關(guān)系訓(xùn)練出適用于區(qū)分基礎(chǔ)心音和非基礎(chǔ)心音的GMM模型，求得期望ηn及方差σn；

（4）計算DHMM的初始概率分布πc，狀態(tài)轉(zhuǎn)移概率矩陣δ，以及各狀態(tài)可能的觀察值對應(yīng)的概率分布矩陣Ψ，結(jié)合GMM模型，生成分割模型λ={πc,δ,Ψ,ηn,σn}；

（5）對狀態(tài)跳轉(zhuǎn)的判斷及修正，得到基礎(chǔ)心音及非基礎(chǔ)心音的區(qū)分結(jié)果；

（6）通過心音信號的時域特征識別出收縮期及舒張期，區(qū)分出s1，收縮期，s2和舒張期，得到心動周期內(nèi)的分割結(jié)果。

圖6 心音分割算法流程圖Fig.6 Flow chart of heart sound segmentation algorithm

4 實驗結(jié)果及分析

實驗選用的低心率變異性數(shù)據(jù)來源于中國醫(yī)學(xué)科學(xué)院阜外醫(yī)院采集提供和“PhysioNet/Computing in Cardiology Challenge 2016”的數(shù)據(jù)集，并在心內(nèi)科從業(yè)人員的指導(dǎo)下對多組實驗數(shù)據(jù)進(jìn)行了人工標(biāo)記。實驗所用心音頻率為2 kHz,均為心率變異率較低的心音。在心率變異性較低的情況下，選用了4類干擾噪聲：脈沖噪聲、Babble類語音噪聲、白噪聲和粉紅噪聲。實驗中所用脈沖噪聲采樣率為1 kHz，Babble類語音噪聲及白噪聲采樣率為1.6 kHz，粉紅采樣率為44.1 kHz，并對4類噪聲作變采樣（上采樣或者降采樣）處理變?yōu)? kHz；幀長為20 ms，幀移為10 ms,窗函數(shù)為hamming窗,所有心音均已采用最小二乘法擬合進(jìn)行趨勢項消除處理。實驗所用電腦CPU型號為Intel(R)Core(TM)i78700kCPU-@3.70 GHz，RAM為16.0 GB，硬盤為1 T，仿真軟件為MATLAB2017a。為了檢驗分割算法的有效性及可靠性，驗證了4類噪聲在不同信噪比下對算法分割精度的干擾，并與基于邏輯回歸的隱半馬爾可夫模型(Logistic regression hidden semi-Markov model,LRHSMM)的算法[8]和基于Hilbert的算法[4]進(jìn)行了對比。

實驗中，共處理238條樣本信號，約8 000個心動周期。由于部分心音信號并不是太理想，故采用了4個評價指標(biāo)[20]對算法進(jìn)行評價：算法處理耗時（PS）、基礎(chǔ)心音檢出正確率（TPS12）、基礎(chǔ)心音錯誤檢入率（FPS12）和基礎(chǔ)心音未檢出率（MPS12）。計算方法為

式中：TNS12為基礎(chǔ)心音正確檢出幀數(shù)，F(xiàn)NS12為基礎(chǔ)心音錯誤檢入幀數(shù)，MNS12為基礎(chǔ)心音未檢出幀數(shù)。計算中，以人工標(biāo)注的數(shù)據(jù)為參考標(biāo)準(zhǔn)。

表3中給出了測試白噪聲干擾，SNR分別為15，10，5和0 dB下3種算法的基礎(chǔ)心音檢出正確率，圖7中給出了不同信噪比4類噪聲下的檢出正確率。

表3 白噪聲下的評價指標(biāo)對比Tab.3 Comparison of evaluation indicators under white noise

圖7 不同信噪比4類噪聲干擾下的檢出正確率Fig.7 Detection accuracy rate under four types of noise interference with different SNRs

圖7中反映出，Hilbert算法在信噪比較高的情況下檢出正確率較高，可達(dá)92%以上，具有較好的分割性能，但隨信噪比降低，其檢出正確率也明顯下降,且波動較大，抗噪性能相對較弱。在信噪比較高時，LRHSMM算法檢出正確率可達(dá)94%左右，略高于Hilbert方法檢出正確率，但隨信噪比降低LRHSMM算法的檢出正確率明顯高于Hilbert方法，且檢出正確率波動較小，抗噪性能較好。本文提出的DHMM方法檢出正確率可達(dá)94%以上，略高于LRHSMM方法。這是由于本文算法在低心率變異下通過GMM模型進(jìn)行約束能更好地匹配心音信號特征，降低了泛化效果對分割精度的干擾，且不受受試者間心動周期差異的影響。

在PS指標(biāo)中，取平均作為PS評價值。Hilbert方法每條心音檢出耗時在6.7 s左右，耗時短，實時性好，但復(fù)雜聲學(xué)環(huán)境下分割精度相對較低；LRHSMM方法每條心音檢出耗時在30 s左右，運(yùn)算時間明顯高于Hilbert算法，但其分割精度也高于Hilbert算法；本文DHMM方法每條心音檢出耗時在25 s左右，運(yùn)算時長較LRHSMM短，且LRHSMM方法通過大量訓(xùn)練集進(jìn)行訓(xùn)練，耗時耗力，而本文算法采用在線訓(xùn)練方法更為簡便，其分割精度也較高。

圖8給出了原始心動周期與4類噪聲干擾下的心動周期在DHMM中分割結(jié)果的對比。從圖8（ad）分割結(jié)果與原始心動周期分割結(jié)果可以直觀地看出，本文算法在各類噪聲干擾下依然具備良好的分割性能。

5 結(jié)束語

在心率變異性較低的心音信號分割中，DHMM較好地結(jié)合了心音信號的準(zhǔn)周期性、短時平穩(wěn)性及其主要成分在時域下的分布特點。DHMM簡化了HMM中的狀態(tài)轉(zhuǎn)移過程，由s1，收縮期，s2和舒張期間的轉(zhuǎn)換直接變?yōu)榛A(chǔ)心音及非基礎(chǔ)心音間的轉(zhuǎn)換。通過自相關(guān)函數(shù)可較準(zhǔn)確地提取出心音周期，采用EEMD提取模態(tài)本征信號表征心音信號，分析心音短時平均能量分布直方圖可建立GMM模型以區(qū)分基礎(chǔ)心音和非基礎(chǔ)心音，當(dāng)基礎(chǔ)心音的判斷滿足關(guān)聯(lián)時間要求時方可發(fā)生狀態(tài)跳轉(zhuǎn)。該過程較好地描述了心音信號的狀態(tài)轉(zhuǎn)移過程，且在低心率變異情況下，本文提出的方法能更好地匹配心音信號特征，較LRHSMM及Hilbert方法更適用于心率變異性較低的心音定位分割。后續(xù)工作可對如何提高算法對心音變異容忍度進(jìn)行研究，以提高算法的適應(yīng)能力，擴(kuò)大算法的適用范圍。

圖8 SNR=10 dB時4類噪聲干擾下的分割效果Fig.8 Segmentation effect under four types of noise interference with SNR=10 dB