基于Hilbert邊際能量的舒張期心音診斷算法*

任 靜， 趙治棟

(1.杭州電子科技大學 通信工程學院，浙江 杭州 310018； 2.浙江省杭電智慧城市研究中心，浙江 杭州 310018)

0 引 言

心音是一種重要的生理信號，是由心臟舒縮運動時大血管和心臟瓣膜等因素機械振動產生的，蘊含了豐富的生理病理信息。心臟產生病變會導致心臟內部結構發生變化，從而影響和改變心音信號。早期研究[1]表明冠脈狹窄會導致舒張期高頻心雜音增加，這為冠心病的無創智能診斷提供了理論依據。近年來，冠心病診斷技術成為研究熱點，國內外學者從時域、頻域、時頻域、非線性動力學等多個角度研究心音信號。其中在頻域能量方面，Gauthier D[2]采用快速傅里葉變換對舒張期心音信號進行分析，根據觀察正常人和冠心病患者的幅頻譜，以130 Hz作為固定頻率閾值，實驗證明冠心病患者的高頻能量占低頻能量的比值高于正常人。但由于個體心音幅頻譜存在差異，單一固定的閾值無法正確區分高低頻，另外傅里葉變換存在局限性。將傅里葉頻譜和邊際譜相比較可以發現，在傅里葉頻譜中，要求任意頻率都具有相同的幅值，容易出現虛假頻率，產生能量泄露[3]等問題，但是邊際譜涵蓋了信號整體頻率的幅值，任意頻率可能以不同或者相同的幅值出現在信號中。所以邊際譜能夠比較準確地反映信號的實際頻率分布情況[4]。

本文提出了基于希爾伯特—黃變換(Hilbert-Huang transform，HHT)[5]和能量譜分析相結合的方法，對舒張期心音信號進行HHT,根據不同能量譜的差異進行頻域分析，獲取多頻率心音特征，為冠心病的臨床診斷提供新思路。

1 HHT基本原理

HHT是由美籍華人黃鍔等人提出的一種時頻分析方法。主要分成2步：經驗模式分解(empirical mode decomposition，EMD)和Hilbert變換(Hilbert transform)。該方法非常適合對非平穩、非線性信號進行分析，因為其可得到極高的時頻分辨率。目前,廣泛應用于各科研領域。

1.1 EMD

EMD[6]適用于將非平穩信號分解成若干個平穩的具有一定物理意義的固有模態函數(intrinsic mode function，IMF)。從各IMF分量中可以提取能夠反映信號本質的特征信息，相比較從原始信號中直接提取特征的方法，該算法有利于減少特征信息之間的干涉或者耦合。另外，每個IMF分量需要滿足以下2個條件[7]：1)函數中所有的極值點與過零點的數目相等或者相差1；2)在函數上任意一點，局部極大值包絡與局部極小值包絡構成的均值包絡為零。EMD算法主要是通過不斷篩選和迭代獲取滿足條件的IMF分量，分解公式如下

(1)

式中ci(t),i=1,2,…,n為第i個IMF分量;rn(t)為殘余項，表示信號的趨勢;ci(t)依次包含了原始信號從高頻到低頻的成分，反映信號的動態特征;而rn(t)為信號的偏移量和穩態值。

1.2 Hilbert變換

基本原理如下：

1)對每一個IMF分量ci(t)進行Hilbert變換

(2)

2)由ci(t)和yi(t)構成的復信號zi(t)稱為ci(t)的解析信號

zi(t)=ci(t)+jyi(t)=ai(t)ejθi(t)

(3)

式中ai(t)為信號zi(t)的瞬時振幅；θi(t)為信號zi(t)的瞬時相位。

3)ai(t)，θi(t)及瞬時頻率fi(t)可按式(4)計算

(4)

4)Hilbert變換得到的頻率和幅值均為時間函數，可將幅值表示成頻率和時間的函數H(f,t)，即Hilbert幅值譜，簡稱Hilbert譜，記作

(5)

5)邊際譜是對整體信號每個頻率幅值的測量，表示同一頻率下幅值的累積[8]。將H(f,t)對時間積分可得到邊際譜，即

(6)

6)能量譜是對整體信號每個頻率能量的測量，表示同一頻率下能量的累積[9]。將H(f,t)的平方對時間積分可得到能量譜，即

(7)

1.3 主成分分析降維

本文采用主成分分析(principal component analysis，PCA)解決特征冗余問題[10]。PCA的基本原理是將n維特征映射到k維上(k≤n)，構造出k維全新的正交特征，也稱為主成分(principal component)。目的是簡化一些可以被其他特征替換的特征或者是能夠直接省略的特征。對于N×n的訓練數據，有N個樣本，每個樣本包含n個特征，先歸一化獲取參數，計算樣本特征的均值和協方差矩陣[11]，根據協方差矩陣獲取特征值和特征向量，特征向量根據特征值大小由大到小排序，選取前k列作為降維矩陣。對于測試集，采用訓練集歸一化后的參數進行歸一化，通過降維矩陣映射到低維空間上。

1.4 支持向量機分類

支持向量機(support vector machine，SVM)是一種有監督的分類器[12～14]，影響SVM分類精度的主要因素是核函數，本文采用徑向基核函數(radial basic function，RBF)，其滿足不同類別的樣本在特征空間的距離盡可能大，同一類別的樣本間距離盡可能小的需求[7]，并且其主要用于線性不可分的情況。另外，RBF需要確定的參數比較少，所以降低了函數的復雜程度，較易計算。核函數公式如下

(8)

分類過程：將數據區分成訓練集和測試集，對數據歸一化，應用RBF核函數，選取合適的參數σ和懲罰因子訓練訓練集，并使用測試集獲取分類精度。

2 心音信號特征分析

2.1 算法流程

算法流程如圖1所示。

圖1 算法流程

2.2 比值特征定義

取正常、異常舒張期心音各1例，如圖2所示。觀察其傅里葉譜和Hilbert邊際譜如圖3所示，可知傅里葉譜中出現豐富的虛假高頻成分，這是為了彌補信號的非平穩性而引入的，故采用傅里葉譜分析法描述信號的頻域特性是不適合的。分析正常異常舒張期心音的邊際能量譜如圖4所示，可知正常舒張期心音能量主要分布在0～100 Hz之間，包含很少的高頻心雜音；而異常舒張期心音含有較多的高頻心雜音，為冠脈狹窄導致心音異常產生的高頻分量。

圖2 正常和異常舒張期心音信號

圖3 正常異常舒張期心音的傅里葉譜和Hilbert邊際譜

圖4 正常異常舒張期心音的邊際能量譜

能量比值，定義為舒張期心音高頻能量占總能量的比值，由于大于500 Hz的能量太小可以忽略不計，所以，取500 Hz以內的能量作為總能量，高低頻能量的分界頻率為頻率閾值fth，則能量比值為

(9)

2.3 信號長度對特征的影響

分析正常和異常舒張期心音信號在不同長度時間序列下的能量比值的大小變化情況如圖5所示。觀察發現正常心音能量比值在100點時最大，之后下降趨于平穩，整體波動范圍較小，大約為0.14～0.25；異常心音在信號長度大于200點時，比值在0.5～0.6輕微波動。綜合考慮選取N=500作為樣本信號的長度。

圖5 不同數據長度對能量比值的影響

2.4 頻率閾值點序列選取

對36人的舒張期心音數據進行統計分析，其中正常和異常各18例，正常舒張期取360段，異常舒張期取345段，每段舒張期取500點。以個體為單位，提取能量比值均值和標準偏差，分析各頻率閾值下的箱線圖，如圖6所示，圖中每個頻率閾值對應的一組數據中左側為異常，右側為正常。結合能量比值均值和標準偏差的分布情況，選取80～160 Hz為頻率點序列，故將該頻率序列下的均值和標準偏差作為特征輸入分類器。

圖6 80～160 Hz頻率閾值下能量比值均值和標準偏差分布

3 實驗結果分析與討論

本文所用心音數據均由醫院冠脈造影技術確診，采樣頻率2 000 Hz，正常人和冠心病人各18例，采集4個位置心音(舍去不完整顯示的心音)，每個心音選取5段舒張期，正常心音360段舒張期，異常心音345段舒張期，每段舒張期心音取500點，以在頻率序列80～160 Hz下的能量比值的均值和標準偏差作為特征，通過PCA降維，設置訓練測試集不同分類方案，輸入分類器。實驗結果如圖7所示。

3.1 不同PCA特征值對分類結果的影響

采用PCA簡化樣本特征值，選取不同的參數，得到的分類準確率不同，如圖7(a)所示。本文共18個特征：9個頻率點下的特征均值和標準偏差。針對不同數量形式的特征集：1)以舒張期為單位，即特征集由每個舒張期的能量比值構成(k≤9)；2)以整個心音信號為單位，即特征集由每個信號所有舒張期比值的均值和標準偏差構成(k≤18)；3)以個人為單位，即特征集由每個人4個位置所有舒張期比值的均值和標準偏差構成(k≤18)。觀察圖7可知，k取到2時，后兩種特征集分類精度均達到最大值，且以個人為樣本單位的特征集能達到100 %分類準確度。

3.2 不同訓練集數目對分類結果的影響

根據不同的分類方案，即訓練集與測試集數目的比例不同，對分類結果有很大的影響，如圖7(b)所示?？芍柧毤贂r，準確度沒有參考價值；訓練集超過50 %以上，數量越多，分類準確度越高。

圖7 不同主成分數目和不同訓練集數目對分類結果的影響

3.3 不同分類器對分類結果的影響

針對本文中的小樣本[15]，采用不同的分類器，如K-最鄰近規則(KNN)，K均值(K-means)，樸素貝葉斯(Bayesian)，對原始數據進行分類，對比主成分分析和支持向量機相結合的分類方式(PCA+SVM)，如表1所示，可知該分類器分類精度高，分類時間短。

表1 不同分類器下的分類精度和分類時間

4 結束語

在HHT下，對臨床采集的心音數據的邊際能量譜進行分析，將不同頻率閾值下提取的高頻能量占總能量的比值作為特征集，采用PCA刪除冗余特征，通過SVM分類。結果表明異常舒張期心音的能量比值普遍大于正常心音，這與冠心病病理特征符合，即當冠脈堵塞時，舒張期心音會產生高頻分量，導致高頻部分能量增加。結合SVM分類結果，說明該方法可以正確區別兩類心音，有利于冠脈狹窄智能診斷的研究。

參考文獻:

[1] Dock W，Zoneraich S.A diastolic murmur arising in a stenosed coronary artery[J].Amer J Med,1967，42(4):617-619.

[2] Gauthier D,Akay Y M,Paden R G,et al.Spectral analysis of heart sounds associated with coronary occlusions[C]∥6th International Special Topic Conference on Information Technology Applications in Biomedicine,Tokyo,Japan,2007:49-52.

[3] 李建康,周宏月,宋向榮.車輛道路載荷譜的Hilbert邊際譜分析方法及應用[J].汽車技術,2013(10):55-59.

[4] 劉世金.Hilbert-Huang變換及其應用研究[J].高師理科學刊,2012(4):40-42.

[5] 王宏宇.基于希爾伯特—黃變換的語音識別特征提取方法研究[D].廣州：華南理工大學,2012:69.

[6] He Dhengyun，Yi Ding.Marginal Hilbert spectrum based on EMD reconstruction and its application in fault diagnosis of cooling tower[C]∥2013 the Third International Conference on Intelligent System Design and Engineering Applications(ISDEA),2013:926-929.

[7] 黃林洲,郭興明,丁曉蓉.EMD近似熵結合支持向量機的心音信號識別研究[J].振動與沖擊,2012(19):21-25.

[8] Le Shasha,Hu Weiping.Cough sound recognition based on Hilbert marginal spectrum[C]∥2013 the 6th International Congress on Image and Signal Processing(CISP),Hangzhou,China,2013:1346-1350.

[9] Qin Xiaoan,Zeng Xiangjun,Xu Yao,et al.HHT energy spectrum-based non-unit EHV line protection using traveling wave[C]∥The 2010 International Conference on Modelling,Identification and Control(ICMIC),Okayama,Japan,2010:466-470.

[10] Sankar D Sandeep Vara,Roy Lakshi Prosad.Principal component analysis(PCA) approach to segment primary components from pathologocal phonocardiogram[C]∥2014 International Conference on Communications and Signal Processing(ICCSP),2014:910-914.

[11] 胡順仁,李瑞平,包 明,等.基于主元分析的橋梁撓度傳感器故障診斷研究[J].傳感器與微系統,2014,33(6):9-12.

[12] 陶泳任,陳冠雄,沈海斌.便攜心電傳感器的心房肥大診斷研究[J].傳感器與微系統,2014,33(9):55-57.

[13] 江培海，黃啟俊，常 勝，等.基于SOC FPGA的心電檢測系統設計[J].傳感器與微系統,2016,35(2):74-77.

[14] 韓開旭，任偉建.基于改進Fisher核函數的支持向量機在推特數據庫情感分析中的應用 [J].自動化技術與應用,2015,34(11):30-35.

[15] 陳 曦,陳冠雄,沈海斌.基于SVM的ECG傳感器信號身份識別方法[J].傳感器與微系統,2014,33(10):40-42.