許春冬 周 靜 應冬文, 侯雷靜 龍清華

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

引 言

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

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

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

1 心音成分及分割任務要求

心音一般是指由心肌收縮、心臟瓣膜關閉、血液撞擊心室壁和大動脈壁等活動引起的振動經(jīng)胸腔、胸壁等組織傳導至體表所產(chǎn)生的聲學信號[10-11]，是能夠反映人體心臟及心血管系統(tǒng)機械運動狀況的信號之一[12]。心臟結構如圖1所示，心臟通過循環(huán)的收縮活動來為人體不斷輸送血液，其活動狀態(tài)為一周期過程。對于正常人而言，心音信號近似于準周期信號，每個周期約在600 ms～1 s之間[2]，醫(yī)學上稱作心動周期。每個心動周期內(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ī)學上，常認為一個心動周期由s1，收縮期，s2和舒張期4個部分構成，且在任意一個心動周期中均存在，出現(xiàn)順序是固定循環(huán)的，如圖2所示。

圖1 心臟結構Fig.1 Structure of the heart

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

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

2 基于DHMM的心音分割

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

2.1 心音信號的標注

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

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

2.2 面向心音信號的HMM與DHMM

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.3 心動周期的提取

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

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

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

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

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

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

圖4 自相關提取心音周期Fig.4 Autocorrelation extraction heart sound cycle

2.4EEMD分解

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

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

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

2.5 心音信號的高斯建模

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

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

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

模型的約束關系可概括為[18]：

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

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

式中ρ為常數(shù)。

（2） 方差間的約束關系為

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

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

3 心音信號的分割流程

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

（1）輸入待分割的心音信號并作預處理，通過自相關包絡獲取心動周期均值T；

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

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

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

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

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

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

4 實驗結果及分析

實驗選用的低心率變異性數(shù)據(jù)來源于中國醫(yī)學科學院阜外醫(yī)院采集提供和“PhysioNet/Computing in Cardiology Challenge 2016”的數(shù)據(jù)集，并在心內(nèi)科從業(yè)人員的指導下對多組實驗數(shù)據(jù)進行了人工標記。實驗所用心音頻率為2 kHz,均為心率變異率較低的心音。在心率變異性較低的情況下，選用了4類干擾噪聲：脈沖噪聲、Babble類語音噪聲、白噪聲和粉紅噪聲。實驗中所用脈沖噪聲采樣率為1 kHz，Babble類語音噪聲及白噪聲采樣率為1.6 kHz，粉紅采樣率為44.1 kHz，并對4類噪聲作變采樣（上采樣或者降采樣）處理變?yōu)? kHz；幀長為20 ms，幀移為10 ms,窗函數(shù)為hamming窗,所有心音均已采用最小二乘法擬合進行趨勢項消除處理。實驗所用電腦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]進行了對比。

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

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

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

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

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

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

5 結束語

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

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