基于高階統計量的心音特征提取與分類研究

朱燕瓊，潘家華，楊宏波，王威廉

(1. 云南大學信息學院，云南 昆明 650091；2. 云南省阜外心血管病醫(yī)院，云南 昆明 650102)

1 引言

心音信號是人體重要的生理信號[1]，它反映心臟及心血管系統機械運動狀況，含有心臟各部分及相互間作用的生理、病理信息[2]。心臟聽診是先天性心臟病(簡稱先心病)初診的主要手段，但需要醫(yī)生具有豐富的臨床經驗，這一過程具有一定的主觀性，缺乏客觀標準。高階統計量已經在其它生理信號的分析研究方面取得良好效果，ARMA(autoregressive moving average)模型廣泛用于對隨機信號的參數估計。本文采用高階統計量并結合ARMA模型對心音信號通過數字化分析手段，提取心音信號的特征進行分類，以達到機器輔助診斷的目的。

近年來對心音信號的研究，備受國內外學者的關注，心音信號預處理通常包括：去噪、提取包絡、分段定位，以及特征提取等[3，4]。國內利用高階統計量分析心音信號一般采用AR模型或MA模型進行參數估計[5]，國外對心音信號分析研究較國內早，L.Durand在1986年將傳統譜分析和現代譜分析方法應用于植入主動脈位置上的人工生物瓣膜音的分析，并將兩種譜分析的性能進行了比較[6]，Moukadem Ali運用高階統計量和時頻域特征對心臟壓力測試下心音的分類[7]，對嚴重的心動過速和快速性心律失常時，提取其心音的S1、S2的特征進行分析，但利用高階統計量ARMA模型分析心音信號尚未見報道。本研究對采集到的心音信號直接用高階統計量的ARMA模型和梅爾對數系數[8](MFSC)進行特征提取，并用三種K近鄰，決策樹，貝葉斯分類器進行分類研究和比較。

2 高階統計量的ARMA模型提取特征

2.1 高階統計量

高階統計量(high-order satistics，簡稱HOS)是指階數大于二階的統計量，如高階矩、高階累積量、高階矩譜，以及高階累積量譜。這些統計量對非高斯性、非最小相位、有色噪聲或非線性等必須考慮的問題有著特殊作用，能揭示常規(guī)功率譜分析所不能表現得重要信息[9]。對于高斯的信號，其高階累積量為0，因而可有效抑制高斯有色噪聲。心音信號也是一種非高斯、非線性、非平穩(wěn)的生理信號，傳統時頻域分析難以準確表征信號特征。高階統計量已在腦電信號、心電信號、肌電等生理信號分析上取得好的效果[10]，還廣泛應用于各個實際信號處理領域。例如，通信、聲納、雷達、生物工程、語音處理、地震信號處理、圖象處理、時延估計、系統辯識、自適應濾波以及陣列處理等領域都得到廣泛應用[9]。

r階矩和r階累積量的一般表達式，分別表示為式(1)、(2)

mk=mk1 k2…kn=mom(x1，x2，…，xk)

(1)

ck=ck1 k2…kn=cum(x1，x2，…，xk)

(2)

對隨機過程{x(n)}，根據高階累積量定義式(2)可知，{x(n)}的各階累積量為

C1x=E{x(n)}=0

(3)

C2x(τ)=E{x(n)x(n+τ)}=r(τ)

(4)

Ckx(τ1，τ2，…τn-1)=0，k≥3

(5)

由高階矩定義式(1)可推知，只有奇數階的高階矩才等于零，而偶數階的高階矩不等于零，即

(6)

高階統計量中，最簡單、最常用的是三階累積量和三階譜(又稱雙譜)，設高階累積量Ckx(τ1，τ2，…τn-1)絕對可和，即

(7)

則隨機信號{x(n)}的k階累積量譜Cnx(ω1，ω2，…，ωn-1)存在并連續(xù)，且定義為k階累積量的(k-1)維Fourier變換，即

Skx(ω1，ω2，…，ωn-1)=

(8)

式中k=3≤π，|ω2|≤π，|ω1+ω2|≤π。

據以上分析，本研究將高階統計量應用到心音信號的分析中，利用高階累積量抑制高斯有色噪聲，目的在于提取更能反映心音本質的特征，獲取更有效的分析結果。高階累積量具有相關函數無法比擬的優(yōu)點，能夠用來辨識非最小相位系統。

2.2 ARMA模型

心音信號被認為是一種隨機時間序列，為隨機信號建立參數模型是研究隨機信號的一種基本方法，ARMA模型是研究時間序列的重要方法，在自然科學和經濟學中得到了廣泛的應用。ARMA模型由自回歸模型與移動平均模型為基礎，是現代譜估計中的參數譜估計[11]。功率譜估計是數字信號處理的重要內容和重要應用點之一，主要研究信號在頻域中的各種特征，在有限數據中提取淹沒在噪聲中的有用信號。信號功率譜一般可以分為平譜、線譜和ARMA譜3種。平譜是白噪聲譜，線譜是一個或多個正弦信號功率譜，ARMA譜介于平譜和線譜兩種譜之間，是一個零極點模型，既有峰值又有谷點的譜，而AR模型則僅易于反映譜中的峰值，MA僅易于反映譜中的谷點。根據以上分析采用ARMA模型分析和計算心音信號的特征值。綜上所述，將高階統計量ARMA模型用于心音信號分析的優(yōu)點：

1)高階統計量ARMA模型不僅保留心音信號幅度信息，同時也保留了相位信息；

2)高階統計量ARMA模型具有更好的抗噪性，根據式(2)、(5)高斯信號的高階累積量C3，X(τ1，τ2)=0為，理論上能夠抑制高斯噪聲的影響，改善信噪比；

在觀測數據較短的情況下，要獲得高分辨率的雙譜估計，必須使用參數化雙譜估計方法，非參數化多譜估計是低分辨率的[15]。因此本實驗采用ARMA參數模型法對心音信號進行雙譜分析。設序列x(n)滿足三階平穩(wěn)，由式(3)給出ARMA(p，q)模型：

(9)

式中e(n)是具有有限的非零累積量的高階白噪聲。

2.3 心音信號提取特征

在matlab高階統計量工具箱中，q-slice(對角切片)又稱切片雙譜或1.5維譜，可以去除非線性耦合，而且具有抑制噪聲的能力。 殘差(residual)是模型的擬合值(fitted value)與實際觀測值之差。實驗中調用高階統計量工具箱中armaqs函數使用q-slice算法，調用armarts函數使用殘差時間序列估計ARAM模型參數，函數中的p和q個分量組成了心音信號的特征向量，經實驗對比，選取MA階數(函數中q)為5，AR階數(函數中p)為6時，累積量階數取默認值3，累積量的最大滯后數目為11，FFT的計算長度為256，在二維FFT中，f1和f2分別是FFT矩陣行和列相關的頻率向量，此時由頻率向量，相位對角切片分量和ARMA參數向量組成特征向量，參數個數為1046。預處理的部分分析圖如圖1-圖6所示。

圖1 房間隔缺損心音原信號

圖2 房間隔缺損信號的高階譜平面圖和三維圖

圖3 室缺信號的相位譜 三維圖和相位譜對角切片圖

圖4 正常心音原信號

圖5 正常信號的高階譜平面圖和三維圖

圖6 正常信號的相位譜 三維圖和相位譜對角切片圖

圖2、圖3、圖5、圖6是特征提取過程中，原始房間隔缺損心音信號和正常心音信號及對應高階統計量譜圖。心動周期中，因血流沖擊心室壁和大動脈等因素引起機械振動，產生心音。受環(huán)境噪聲、呼吸音等影響，從時域直接觀察4秒的異常心音信號和正常心音信號區(qū)別微小，但是觀察其對應高階統計量譜圖2和圖5，可以發(fā)現異常心音信號高階統計量譜圖2比正常心音信號高階統計量譜圖5中高頻成分要多，異常心音信號幅度數量級是10-4，正常心音信號幅度數量級是10-7。從相位三維圖3和對角切片圖6可以看出，室缺心音信號相位波動大，且相位起點更多從負相位到正相位?？梢姴捎酶唠A統計量易于對正常和異常心音進行區(qū)別性研究和分析，并提取更有效的特征。

3 梅爾頻率系數

MFCC中的梅爾刻度是一種基于人耳對等距的高音變化的感官判斷而制定的非線性頻率刻度，能較好地反映人耳對聲音的特點[12]。MFSC是省略離散余弦變換步驟的MFCC的特殊形式，直接對梅爾頻率系數取對數能量。MFSC相比較MFCC光譜能量能保持局部特性，在頻譜上更為平滑[13]。

人耳對聲音的感知是非線性的，以取對數的方式可以更好的描述。分幀只能反映本幀的心音信號特性，差分更能體現心音信號時域連續(xù)性，即常用一階差分(Δ)和二階差分(Δ-Δ)實現，如式(10)(11)所示。把得到的mfcc參數，一階差分，二階差分組成一維度矩陣，即M=幀數190*3，經過34(符號表示為N)個濾波器，即得到M×N的二維特征，在保證特征參數完整的前提下，將570*34(在這里M=570為移動的幀數，34為濾波器個數)二維的矩陣，按每階濾波器依次存放轉化為一維數據。最終得到19380個參數。

Δ=y(x+1)-y(x)

(10)

Δ-Δ=y(x+2)-2y(x+1)+y(x)

(11)

式中x為自變量。Δ是一階差分，Δ-Δ是二階差分。

圖7 心音信號的MFSC數據組織

MFSC處理的基本步驟：

1)預加重。將信號x(n)通過一個高通濾波器，其目的是補償高頻部分，突顯高頻的共振峰。

2)分幀、加窗。分幀是短時非平穩(wěn)的心音信號近似為平穩(wěn)信號，加窗可以克服頻譜泄露現象。分幀時信號交疊一般為50%，以提高時間分辨率。分幀的幀數M計算如式(12)所示

(12)

3)頻域變換。信號在時域上很難觀察其能量特性，通常把它轉換成頻域上的能量分布，能量分布來代表心音信號的特性。本文使用傅里葉變換，其計算公式如式(13)所示：

(13)

4)梅爾刻度轉換。由于頻域信號有很多冗余，梅爾濾波器組可以對頻域進行精簡，一個濾波器產生一個頻段值，每一個頻段用一個值表示。

5)能量取對數、差分。

4 分類識別模型

4.1 KNN

k近鄰(kNN，k-NearestNeighbor)是實踐中最成功的算法之一，它在理解、實踐上較簡單，但產生了令人難以置信的良好結果。它是一種通用分類器，在蛋白質分類，計算機視覺，計算幾何和圖形等方面都有應用。KNN是基于實例的學習(也稱為惰性學習)，其中函數僅是局部近似的，所有計算都會延遲到分類為止。在KNN中，給定測試樣本，基于某種距離度量找出訓練集中與其最靠近的k個訓練樣本，然后基于這k個“鄰居”的信息進行預測。任何給定實例均按其“ k”個最接近的鄰居按多數票進行分類，其中“ k”是一個重要參數，是正整數，通常很小，當k取不同值時，分類結果會有顯著不同。給定測試樣本x，若其最近鄰樣本為z，則最近鄰分類器出錯的概率就是x與z類別標記不同的概率[14]，如式(14)所示

(14)

4.2 決策樹

決策樹是一種樹形結構的分類器，通過順序詢問分類點的屬性決定分類點最終的類別。本文根據心音信號的MFSC特征和高階譜特征，構建一顆決策樹。分類是按照結點依次進行判斷，即可得到樣本所屬類別。隨著劃分過程不斷進行，希望決策樹的分支結點所包含的樣本盡可能屬于同一類別。即結點的“純度越高”， “信息熵”是度量樣本集和純度最常用的一種指標[14]。假定當前樣本集D中第k類樣本所占的比例為pk(k=1，2，…，|y|)，則D的信息熵定義如式(15)所示

(15)

4.3 樸素貝葉斯

樸素貝葉斯分類器是一個以貝葉斯定理為基礎的多分類的分類器。對于給定數據，首先基于特征的條件獨立性假設，學習輸入輸出的聯合概率分布，然后基于此模型，對給定的輸入，利用貝葉斯定理求出后驗概率最大的輸出，如式(16)。

(16)

式中P(c)是類先驗概率，P(c│x)是樣本x相對于類標記c的類條件概率或稱“似然”。

5 實驗設計與評估標準

5.1 實驗設計

本節(jié)重點對比討論本文的高階統計量特征提取和梅爾倒譜系數特征提取算法，主要研究內容是在matalb和tensorflow平臺上實現。研究所用的數據源為，從云南省阜外心血管病醫(yī)院臨床從和全國各地先心病篩查采集的心音信號，實驗取其中1000個心音樣本，800用于訓練，其中正常心音樣本600個，先心病心音樣本200個 ；200個樣本進行對比測試，該測試集中有正常心音樣本100個，先心病心音樣本100個。實驗心音信號時長為20s，采樣率為5kHz。

5.2 評估標準

二分類算法常采用靈敏度(符號記為：se)、特異度(符號記為：sp)和準確率(符號記為：acc)作為算法評價指標。根據金標準準則，表示正確分類的異常信號為真陽性(符號記為：TP)；錯誤分類的異常信號為假陽性(符號記為：FP)；表示正確分類的正常信號為真陰性(符號記為：TN)；表示錯誤分類的正常信號為假陰性(符號記為：FN)，靈敏度、特異度和準確率計算如式17‐式19所示。

(17)

(18)

(19)

6 實驗結果與分析討論

6.1 實驗結果

實驗中所用的1000例心音數據均來自經臨床確診的心音信號，其中訓練集800例，測試集200例。利用高階統計量的ARMA模型直接提取先心病心音信號的特征，用K近鄰，決策樹，貝葉斯分類器對其進行分類識別，并與用梅爾對數系數提取特征的方法進行了對比，結果如表1所示。從數據庫中隨機選取正常和異常各200列原始信號進行t檢驗，結果P<0.01，表明正常和異常數據具有顯著性差異。

表1 各類算法的比較

表中，分類結果sp為特異度，se為靈敏度，acc為準確率。

6.2 分析討論

從實驗結果可以看出K鄰近，決策樹，貝葉斯三種機器學習分類模型效果的差異性，貝葉斯分類在MFSC和高階統計量中均不理想， KNN沒有決策樹分類結果好，原因是本實驗中對于MFSC的特征來說k取值小了，但為保證實驗結果的可比性，k值根據訓練集綜合考慮決定。決策樹更適宜于這兩種特征提取后心音病例的分類，針對未去噪截取的原始心音信號，與MFSC相比高階統計量ARMA模型總體效果好。

7 結論

本文提出基于高階統計ARMA模型直接提取心音特征，省去了一般預處理中，對心音的降噪和分段定位過程，提高了分析效率，降低了特征參數數量，加快了運算速度。運用KNN，決策樹和貝葉斯三種不同的分類器進行分類判別評測，發(fā)現決策樹更適合ARMA和MFSC這兩種特征提取后心音病例的分類。實驗發(fā)現高階統計量ARMA模型與梅爾對數倒譜(MFSC)相比，針對未經過去噪，分段定位的原始心音信號，高階統計量ARMA模型在測試集上更有競爭力，計算速度更快。此方法預處理方法簡便，為機器輔助診斷和遠程聽診提供一種好的選擇。