一種多閾值融合心音遞歸圖的設(shè)計與應(yīng)用

成謝鋒， 汪 晶， 王 悅

(1. 南京郵電大學(xué) 電子科學(xué)與工程學(xué)院，南京 210003；2. 江蘇省射頻集成與微組裝工程實驗室，南京 210003)

心音信號是人體重要的生理信號之一，它包含了關(guān)于心臟各個部分如心房、心室、心血管以及各個瓣膜功能狀態(tài)的生理信息，具有普遍性、穩(wěn)定性、獨(dú)特性和可采集性等生物特征[1]。課題組在心音信號的基礎(chǔ)上構(gòu)建了心血管健康評估系統(tǒng)(進(jìn)入了研究生電子設(shè)計大賽全國總決賽并獲獎)-“生命心衣”，為了達(dá)到快而高效的效果，本文從二維的角度對心音部分進(jìn)行分析。

心音信號是典型的非線性信號，并且具有混沌特性，而遞歸圖是分析時間序列周期性、混沌性的一種重要的非線性分析方法[2]。但在利用遞歸圖進(jìn)行分析時存在兩個關(guān)鍵性難題：第一，遞歸圖適合對較短數(shù)據(jù)的處理，一旦數(shù)據(jù)太長，就會明顯增加數(shù)據(jù)的處理時間，心音是準(zhǔn)周期性信號，一般對心音信號的分析至少是2～3個周期，甚至更長,因此，必須對心音信號進(jìn)行前期處理；第二，遞歸圖對閾值參數(shù)的要求很高，如果閾值選擇過小，遞歸圖中將幾乎沒有遞歸點(diǎn)，也難以提供有用的遞歸結(jié)構(gòu)特征；若閾值選擇過大，某些遞歸結(jié)構(gòu)特征將會被掩蓋，從而出現(xiàn)錯誤的遞歸特征。對于遞歸圖閾值的選擇，很多學(xué)者進(jìn)行了相關(guān)的研究工作[3-5]，但這些方法對于心音特征的表征不一定有效，因此，尋找一種適用于心音信號分析的閾值至關(guān)重要。

為了解決上述存在的兩個難點(diǎn)，本文提出了相應(yīng)的解決方法。首先，為了解決數(shù)據(jù)長度的問題，前期需要對心音信號進(jìn)行等分等長處理，既選出心音信號所在頻段,又保留了心音信號適當(dāng)?shù)臄?shù)據(jù)長度,并得到用于遞歸圖處理的等分等長心音系數(shù)。其次，針對遞歸圖閾值的選取問題，本文提出了一種自適應(yīng)閾值的獲取方法，并且考慮到心音信號的復(fù)雜性，通過單個閾值獲取的遞歸圖在表征心音信號的特征時將出現(xiàn)困難,因此，在自適應(yīng)閾值的基礎(chǔ)上構(gòu)建了一種多閾值融合心音遞歸圖，這種心音遞歸圖能最佳的表現(xiàn)心音特征。通過從多閾值融合心音遞歸圖中提取D′2/S′2值和灰度共生矩(Gray-Level Co-occurrence-Matrix, GLCM)的4種特征：能量、對比度、相關(guān)性和熵，并將這些特征作為支持向量機(jī)(Support Vector Machine, SVM)的輸入進(jìn)行分析。

1 分析的前提

圖1中，Ⅰ為正常心音的波形； 波形中的T1為第一心音S1的持續(xù)時間；T2為第二心音S2的持續(xù)時間；T12為S1～S2的持續(xù)時間；T21為S2～S1的持續(xù)時間。Ⅱ為Ⅰ對應(yīng)的遞歸圖。為了明確兩者之間的關(guān)系，并考慮到Ⅰ和Ⅱ的橫坐標(biāo)點(diǎn)數(shù)不一致，因此，將Ⅱ中縱橫坐標(biāo)擴(kuò)展成與Ⅰ的橫坐標(biāo)長度一致，然后將Ⅱ中橫坐標(biāo)方向上的矩形圖塊左右邊界線與Ⅰ對照分析,并用直線標(biāo)記，逆時針旋轉(zhuǎn)90°，即得到Ⅲ。結(jié)合Ⅲ，并對Ⅰ和Ⅱ中的圖形同時從左至右進(jìn)行對比，Ⅰ中第一次S1持續(xù)的時間與Ⅱ中第一次出現(xiàn)矩形圖塊左側(cè)的空白間隔相等，Ⅰ中第一次S1～S2的持續(xù)時間與Ⅱ中第一次出現(xiàn)矩形圖塊的寬度相等，Ⅰ中第一次S2持續(xù)的時間與Ⅱ中第一次出現(xiàn)矩形圖塊右側(cè)的空白間隔相等，Ⅰ中第一次S2～S1的持續(xù)時間與Ⅱ中第二次出現(xiàn)矩形圖塊的寬度相等。以此類推，這種關(guān)系表現(xiàn)在Ⅲ中，即從上往下的直線依次對應(yīng)Ⅱ中從上至下矩形圖塊的左右邊界線。當(dāng)只考慮Ⅱ中對角方向上的情況時，可以得到Ⅳ中的對角化圖形。將Ⅰ與Ⅳ進(jìn)行對比可知，Ⅳ小方塊的邊長T12′與Ⅰ中的T12相等，Ⅳ大方塊的邊長T21′與Ⅰ中的T21相等，因此對于時域波形上T12/T12可以通過遞歸圖中T12′×T12′/T21′×T21′得到， 時域上S1/S2可以通過Ⅳ中T12′/T21′得到。由于一維波形上T12/T12和S1/S2在二維遞歸圖上得到了放大，因此利用遞歸圖對心音進(jìn)行分析是可行的。

圖1 心音與遞歸圖的關(guān)系Fig.1 The relationship between heart sound and recurrence plot

在一個心動周期內(nèi)，可以分為收縮期(systolic,S)和舒張期(diastolic,D)兩部分，其中收縮期持續(xù)時間為S1的開始到S2開始前， 舒張期持續(xù)時間為S2的開始到下一個心動周期S1開始前， 將T1，T2，T12，T21與心音的收縮期和舒張期聯(lián)系起來，可以得到S=T1+T12，D=T2+T21， 為了便于后續(xù)分析令S″=T12,D″=T21。

因此，本文設(shè)計的一種多閾值融合心音遞歸圖和分類識別方法的原理圖，如圖2所示。

圖2 系統(tǒng)框架圖Fig.2 The frame diagram of system

2 分析方法

2.1 獲取等分等長心音系數(shù)

利用RP對心音信號進(jìn)行處理時存在的第一個關(guān)鍵難題，即：RP適合對較短數(shù)據(jù)的處理，數(shù)據(jù)過長會明顯增加數(shù)據(jù)的處理時間，如何降低數(shù)據(jù)的長度。

定義1對信號s(t),通過某種變換降低數(shù)據(jù)的長度,獲取信號所在的頻段，消除其它頻段雜音,稱為等分處理;再對所選信號進(jìn)行多層等長度分析,凸顯信號的時頻細(xì)節(jié)特征, 稱為等長處理;這種既保留原始信號s(t)的所有特征信息同時又降低信號長度的方法，稱為等分等長處理。

設(shè)對信號s(t)等分p段，ε(t)為單位階躍信號,bk(k=1,2,…,p)為等分處理系數(shù)，則經(jīng)過等分處理后，任意一段信號可表示為

(1)

任意選擇等分處理后的一段信號進(jìn)行等長處理, 得到q個等長信號，則等長處理后的任意一個信號可表示為

ski(t)=aisk(t)i=1,2,…,q

(2)

式中：ai為等長處理系數(shù)。

定義2設(shè)Hc(t)是時序心音信號經(jīng)過等分等長處理得到的，在保留原始心音信號的所有特征信息的同時降低了心音信號的長度，則Hc(t)稱為等分等長心音系數(shù)。

等分等長心音系數(shù)Hc(t)具有如下特性：

(1)Hc(t)不僅保留了心音信號的周期性，能量集中性，相對穩(wěn)定性等特征，并且頻段為心音信號的主要頻段范圍，凸顯了心音的時頻特征；

(2)Hc(t)的數(shù)據(jù)長度相對于原始心音信號的長度明顯降低了，同時又避免了心音信號長度過短的情況，很好地滿足了遞歸圖的處理要求。

等分等長心音系數(shù)Hc(t)獲取步驟：

步驟1任意選擇待分析心音信號中4～6 s長度的信號；

步驟2為了獲取心音信號所在的頻段，采用心音小波對原始心音信號進(jìn)行m層小波分解，使得第m層的頻帶低于500 Hz, 并令第m層的低頻信號為Hl(t), 高頻信號為Hh(t)；

步驟3為了凸顯心音的時頻特征，對等分處理后的Hl(t)進(jìn)行d層小波包分解， 得到u=2d個信號， 即Hl1(t)，Hl2(t)，…，Hlu(t)， 選擇經(jīng)過等分等長處理后所有長度一致的， 得到2k+1 個信號，即等分等長心音系數(shù)Hc(t)。

2.2 多閾值融合心音遞歸圖的獲取方法

利用RP對心音信號進(jìn)行處理時存在的第二個關(guān)鍵難題，即：RP對閾值參數(shù)的要求很高，閾值的設(shè)置直接影響RP對心音信號內(nèi)部結(jié)構(gòu)的呈現(xiàn)情況,如何選擇合適的閾值。

2.2.1 一種自適應(yīng)閾值θ的獲取方法

設(shè){x1,x2,x3,…,xn}為時序心音信號， 嵌入維數(shù)為m， 延時時間為τ， 則相空間重構(gòu)后的矩陣為

(3)

令N=n-(m-1)τ，對相空間矩陣任意兩列之間求距離，得到距離矩陣為

(4)

距離矩陣反映了在相空間中，每一個向量在相空間的位置狀態(tài)與其他向量之間的距離變化情況，因此從距離矩陣出發(fā)選擇合適的閾值更具有合理性。令Yi={di1,di2,…,dii,…,diN}， 其中i=1,2,…,N，為了反映距離矩陣D中每一個向量Yi的統(tǒng)計分布，對其求均值，即為

(5)

對距離矩陣中每一個向量Yi進(jìn)行均值分析后， 得到了均值向量，即:E={E1,E2,…,Ei,…,EN}T。 為了便于數(shù)據(jù)分析，只保留小數(shù)點(diǎn)后四位，同時去掉均值向量E中重復(fù)的元素，可得到變換后的互異均值向量，即Et={E1,E2,…,Ei,…,Ep}T,p

θ=Ei|max(W)

(6)

2.2.2 基于自適應(yīng)閾值θ的多閾值融合心音遞歸圖

本文提出的閾值選擇方法是選取權(quán)重最大的元素對應(yīng)的均值，當(dāng)權(quán)重最大的元素出現(xiàn)重復(fù)現(xiàn)象時，閾值選擇將會出現(xiàn)困難。因此，本文在自適應(yīng)閾值的基礎(chǔ)上提出了一種多閾值融合心音遞歸圖。

定義3設(shè)R為等分等長心音系數(shù)經(jīng)過多個閾值處理得到的遞歸矩陣的融合，則稱之為多閾值融合遞歸矩陣。

為了便于對多個閾值的分析，將Et改寫為Et1={E1,E2,…,Ei,…,Ep}T，W改寫為Wt1={w1,w2,…,wi,…,wp}。首先進(jìn)行閾值分類，選擇Wt1中最大元素在Et1中對應(yīng)的值作為第一類閾值e1f={ε11,ε12,ε13,…,ε1f}，其中f

選取不同個數(shù)的閾值進(jìn)行融合處理得到的多閾值融合遞歸圖在表征心音特征上是不同的。灰度共生矩陣的灰度級數(shù)一般取8級或16級[7]，考慮到運(yùn)算效率，因此，本文選擇了前8個閾值進(jìn)行分析，并按照只選擇前兩個，只選擇前三個,只選擇前4個，只選擇前5個，只選擇前6個，只選擇前7個，只選擇前8個，即7種情況進(jìn)行分析。對數(shù)據(jù)庫中的心音采用這7種情況進(jìn)行多閾值融合處理，依次對得到的多閾值融合心音遞歸圖提取D′2/S′2值和灰度共生矩特征，并按照十折交叉驗證[8]的方法對文中兩個數(shù)據(jù)庫的normal和abnormal情況進(jìn)行分析，結(jié)果如圖3所示。

圖3 7種閾值選擇方法的識別率分析Fig.3 The recognition rate analysis of seven threshold selection methods

根據(jù)文獻(xiàn)[9]中的互信息法和Cao算法計算得到嵌入維數(shù)和延時時間參數(shù)，設(shè)定為m=3，τ=1，同時選擇三個閾值進(jìn)行多閾值融合處理。圖4為對一例正常心音進(jìn)行等分等長處理后得到的其中一個等分等長心音系數(shù)Hc1(t)進(jìn)行多閾值遞歸圖融合的分析。

圖4 等分等長心音系數(shù)的遞歸圖融合原理圖Fig.4 The schematic diagram of recurrence plot fusion on equal part and equal length heart sound coefficient

圖4(a)為歸一化心音經(jīng)過等分等長處理后得到的等分等長心音系數(shù)Hc1(t)； 圖4(b)為閾值為e1f時， 對應(yīng)的遞歸圖R1；圖4(c)為閾值為e2v時，對應(yīng)的遞歸圖R2； 圖4(d)為閾值為e3z時，對應(yīng)的遞歸圖R3； 圖4(e)為三個閾值得到的遞歸圖進(jìn)行融合后的多閾值融合遞歸圖。從圖4(e)中可以看出，經(jīng)過多閾值處理得到的多閾值融合遞歸圖不僅保留了心音特征T1，T2，S″和D″，同時S″和D″的長度更加精確了。

2.3 特征提取

2.3.1 提取多閾值融合心音遞歸圖的D′2/S′2值

將得到的多閾值融合心音遞歸圖進(jìn)行對角化處理，則D′2/S′2值的獲取步驟為：

步驟1由于采集的心音信號不都是周期的整數(shù)倍，得到的S″和D″的個數(shù)不都相等，首先統(tǒng)計多閾值融合心音遞歸圖中S″和D″的個數(shù)，若兩者相等不做舍去操作，若兩者不等，舍棄較多的個數(shù)，使兩者保持相等；

步驟2考慮到心音是非平穩(wěn)準(zhǔn)周期信號，隨機(jī)選擇其中一個S″和D″會造成較大的誤差性，因此，分別對步驟1中的S″和D″進(jìn)行求和處理，得到S′=ΣS″，D′=ΣD″；

步驟4分別對S′和D′進(jìn)行平方處理，得到S′2和D′2，并獲取最終的特征值value=D′2/S′2。

2.3.2 提取多閾值融合心音遞歸圖的GCLM

為了獲取多閾值融合心音遞歸圖中不同灰度級之間的空間相關(guān)性，采用灰度共生矩陣提取紋理特征。Haralick等定義了14種描述灰度共生矩陣[10]的特征參數(shù)，本文選取最常用的4種提取多閾值融合心音遞歸圖的紋理特征。

(1) 二階矩

(7)

(2) 對比度

(8)

(3) 相關(guān)性

(9)

式中：μ1，μ2，σ1，σ2分別為均值和標(biāo)準(zhǔn)差。

(4) 熵

(10)

通過對式(4)的分析可知，距離矩陣是關(guān)于斜率為-1的對角線對稱的，即得到的遞歸圖是關(guān)于135°角對稱的，因此在計算灰度共生矩陣時，只用考慮0°，45°，135°三個方向的紋理特征。

3 數(shù)據(jù)分析

3.1 數(shù)據(jù)來源

文中的數(shù)據(jù)庫來源于課題組心音數(shù)據(jù)庫和文獻(xiàn)[11]中的心音數(shù)據(jù)庫，由于遞歸圖可以將D/S值放大，因此將兩個數(shù)據(jù)庫中的異常心音分為收縮期雜音(Systolic Murmur, SM)和舒張期雜音(Diastolic Murmur, DM)。由于數(shù)據(jù)的來源不同，使用的采樣頻率有差異，為了便于數(shù)據(jù)分析，統(tǒng)一將心音信號轉(zhuǎn)換成采樣頻率為2 000 Hz的形式。

3.2 結(jié)果分析

3.2.1 不同閾值選擇方法的對比分析

為了驗證本文提出的自適應(yīng)閾值獲取方法的有效性，分別與文獻(xiàn)[12]中標(biāo)準(zhǔn)差的10%，文獻(xiàn)[13]中標(biāo)準(zhǔn)差的15%，文獻(xiàn)[14]中標(biāo)準(zhǔn)差的25%和文獻(xiàn)[15]中標(biāo)準(zhǔn)差的5倍4種閾值選擇方法進(jìn)行對比分析。并采用十折交叉驗證方法對normal與abnormal的分類以及SM與DM的分類進(jìn)行分析，結(jié)果如表1所示。

從表1可以看出，對于文中的兩個數(shù)據(jù)庫，不管是對正常心音與異常心音的分類，還是對收縮期雜音與舒張期雜音的分類，本文提出的方法較其他4種方法識別效果更好，即本文提出的閾值設(shè)置方法是有效的。

表1 5種閾值選擇方法的識別率對比分析

3.2.2 不同遞歸圖選取方法的對比分析

為了驗證本文提出的多閾值融合心音遞歸圖的可行性，選擇3種方法獲取遞歸圖，方法1選擇一個閾值獲取遞歸圖；方法2選擇三個閾值獲取遞歸圖，不進(jìn)行遞歸圖融合處理；方法3將三個閾值的遞歸圖進(jìn)行融合處理。同樣地，采用十折交叉驗證方法對normal與abnormal的分類以及SM與DM的分類進(jìn)行分析，并采用文獻(xiàn)[16]中的評價指標(biāo)，最終的結(jié)果如表2所示。

表2 識別結(jié)果

從表2可以看出，采用多閾值融合遞歸圖不僅可以有效地實現(xiàn)正常心音與病理心音的分類,還能較好地實現(xiàn)異常心音的分類。其中文獻(xiàn)[11]中數(shù)據(jù)庫的整體識別效果比課題組數(shù)據(jù)庫的整體識別效果更好，原因在于文獻(xiàn)[11]中的心音信號比較干凈，幾乎沒有背景噪聲。

3.2.3 不同算法的對比分析

為了進(jìn)一步驗證本文提出方法的有效性，利用3種算法與本文提出的算法進(jìn)行對比分析，并按照十折交叉驗證方法對normal和abnormal情況進(jìn)行分析，結(jié)果如表3所示。從表中數(shù)據(jù)可知，本文的識別效果較好，即本文的分析方法是有效的。

4 算法應(yīng)用

課題組自主設(shè)計了心血管健康評估系統(tǒng)-“生命心衣”，其中包括心音傳感器，心電傳感器，脈搏傳感器，血氧傳感器。首先利用自主設(shè)計的傳感器實現(xiàn)4路信號的實時采集，并對采集后的信號利用LABVIEW平臺進(jìn)行前期預(yù)處理，包括數(shù)據(jù)存儲和心音信號的時頻分析，然后對采集的數(shù)據(jù)進(jìn)行心音，心電，脈搏，血氧分析，通過提取它們各自的特征參數(shù)構(gòu)建心血管評估指標(biāo)，最后利用得到的多個指標(biāo)構(gòu)建心血管評估模型和預(yù)測模型。系統(tǒng)框架圖如圖5所示。

圖5 生命心衣系統(tǒng)框圖Fig.5 The system diagram of life heart clothing

通過采用本文所述方法對采集的心音信號提取特征，能夠有效減少多路信號之間的影響，并且不用對心音信號進(jìn)行分段處理，就能實現(xiàn)心音的正常與異常分類，提高了計算的實時性，并與同步采集的心電、脈搏、血氧信息一起，在線分析人體血糖(GLU)，平均動脈壓(PM)、心率(HR)、心臟搏血量(SV)、脈搏波傳導(dǎo)時間(PPT)、脈搏傳導(dǎo)速度(PWV)等多個生理參數(shù)，通過對這些生理參數(shù)數(shù)據(jù)的融合處理，構(gòu)建心血管評估模型，并通過對測試者長期測試數(shù)據(jù)的分析，構(gòu)建預(yù)測模型，可預(yù)估測試者未來心血管健康的發(fā)展趨勢。

5 結(jié) 論

針對課題組構(gòu)建的心血管健康評估系統(tǒng)中對心音的分析，采用本文提出的分析方法，可以得到有效地評估效果。因此，本文設(shè)計的一種多閾值融合心音遞歸圖，并提取其D′2/S′2值和灰度共生矩特征進(jìn)行分類識別為心音的分析提供了一種新途徑。

(1) 采用等分等長處理方法進(jìn)行前期處理，避免了對S1，S2進(jìn)行分段操作，同時獲取的等分等長心音系數(shù)，很好地滿足了遞歸圖的處理要求，也是獲取較高識別率的基礎(chǔ)。

(2) 遞歸圖對閾值的要求較高，本文設(shè)計自適應(yīng)閾值的獲取方法，簡單且容易獲取，并且相對于其他閾值選擇方法，識別效果提高了10%～30%。

(3) 多閾值心音遞歸圖的最佳閾值個數(shù)為3，相對于其他閾值個數(shù)，識別效果最好，并且將三個閾值得到的遞歸圖進(jìn)行融合，信號的內(nèi)部結(jié)構(gòu)得到了有效地呈現(xiàn)，遞歸圖的紋理特征更加豐富。

(4) 通過提取多閾值融合遞歸圖的D′2/S′2值和灰度共生矩特征參數(shù)，識別效果可以達(dá)到90%～100%，并且在相同條件下，相對于其他特征提取算法，識別率提高了5%～15%。