基于非平穩(wěn)系統(tǒng)辨識的心音包絡(luò)自適應(yīng)分割

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

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

0 概述

目前,心血管疾病已成為危害人類生命健康的嚴(yán)重疾病,心血管疾病死亡人數(shù)在全球疾病死亡人數(shù)中占據(jù)較大比例[1]。心音信號(Heart Sound Signal,HSS)分析是心血管疾病常用的輔助診斷方法之一,主要通過計算機輔助診斷分析系統(tǒng)的形式進(jìn)行應(yīng)用。心音信號的輔助診斷系統(tǒng)流程可劃分為分割與分類兩個主要步驟[2-3],其中分割是心音信號分析與診斷的基礎(chǔ)與前提,分割結(jié)果對整個系統(tǒng)的分析結(jié)果具有重要影響。然而,心音信號相當(dāng)敏感,極易受到儀器噪聲、環(huán)境噪聲、肺音、呼吸道音、腸鳴音及檢測者自身活動狀態(tài)等因素的影響。此外,各類病理原因也會導(dǎo)致心音中出現(xiàn)心雜音[2-3]。這些干擾因素的存在使得分割任務(wù)變得十分困難,且實際中采集到的心音信號通常包含多種噪聲,導(dǎo)致其可分析性急劇下降,使得原本較弱的心音信號變成一個較復(fù)雜的聲學(xué)信號。因此,心音信號的有效分割逐漸成為目前研究的熱點之一。

隨著心音信號研究的深入,各種分割方法相繼被提出,主流的心音信號分割方法可分為[4-6]基于包絡(luò)的方法、基于特征提取的方法、基于隱馬爾科夫模型(Hidden Markov Model,HMM)的方法和基于機器學(xué)習(xí)的方法4類。其中,基于包絡(luò)的方法因算法簡單、實時性強等特點[6-7]得到廣泛研究與應(yīng)用,特別是在實時檢測需求下成為首選方法[8],而其他3類分割方法的實時性相對較差[5,9-10]。本文借鑒包絡(luò)分割方法的思想,提出一種基于非平穩(wěn)性系統(tǒng)辨識的心音包絡(luò)提取方法,并將其應(yīng)用于心音閾值自適應(yīng)分割中。

1 基于包絡(luò)的分割方法

心音信號分割的首要任務(wù)是定位分割出基礎(chǔ)心音(Fundamental Heart Sound,FHS)信號s1(第一心音)和s2(第二心音),然后分割出收縮期(sys)和舒張期(dia),如圖1所示?；诎j(luò)的分割方法主要包括預(yù)處理、包絡(luò)提取、閾值分割3個核心部分。

圖1 心音時域波形圖Fig.1 Oscillogram of heart sound time domain

預(yù)處理主要包括預(yù)加重、消除趨勢項和降噪,其中:預(yù)加重為高頻提升,通過一階FIR濾波器即可完成[2];對于趨勢項,可采用最小二乘法擬合進(jìn)行消除[11];降噪中應(yīng)用較普遍的方法為基于小波變換的降噪方法,通過丟棄帶外小波系數(shù)、濾波帶內(nèi)小波系數(shù),實現(xiàn)噪聲抑制[12]。預(yù)處理后的心音信號特征更清晰,為包絡(luò)提取奠定了基礎(chǔ)。

心音信號包絡(luò)的提取方法主要包括香農(nóng)能量包絡(luò)法[3,5-6]、維奧拉積分包絡(luò)法[5]、希爾伯特變換包絡(luò)法[7,9]、樣條插值包絡(luò)法[6]以及經(jīng)驗?zāi)B(tài)分解(Empirical Mode Decomposition,EMD)包絡(luò)法[7,10]等。香農(nóng)能量包絡(luò)法在噪聲干擾較低的情況下能夠取得較好的效果,但若降噪后的噪聲干擾依舊較大,則提取的能量包絡(luò)特征會受到嚴(yán)重干擾。維奧拉積分包絡(luò)法能夠較好地提取出s1的特征,但對衰弱s2及心雜音的干擾信號包絡(luò)提取效果明顯下降。在希爾伯特變換法中,當(dāng)屬黃變換研究最深入,但該方法在通過極值點擬合包絡(luò)時通常存在一些固有缺陷,尤其是對于一些異常心音信號。樣條插值包絡(luò)法通過對極大值和極小值進(jìn)行插值獲得上包絡(luò)和下包絡(luò),可以直觀反映信號的幅值,但存在欠包絡(luò)的問題,且實際中欠包絡(luò)問題處理結(jié)果并不理想。EMD包絡(luò)法無需預(yù)先設(shè)定基函數(shù),自適應(yīng)能力較強,但在處理非平穩(wěn)信號時存在模態(tài)混疊及信號失真等問題。

分割過程主要包括閾值函數(shù)設(shè)計與分割成分判定兩部分。其中,閾值函數(shù)多數(shù)為雙閾值函數(shù)[7];任何一個心動周期中均包含第一心音、收縮期、第二心音、舒張期4個部分,因此通過閾值分割判定基礎(chǔ)心音后,根據(jù)時域特征即可區(qū)分出收縮期和舒張期[6-7,9]。

2 非平穩(wěn)性系統(tǒng)辨識原理

系統(tǒng)辨識主要是根據(jù)輸入與輸出時間序列來描述系統(tǒng)行為的過程,用于估計與建立控制領(lǐng)域內(nèi)的系統(tǒng)模型[13],其采用的非平穩(wěn)性系統(tǒng)辨識模型衍生于波束形成等信號處理領(lǐng)域[13-14]。文獻(xiàn)[14]提出一種估計傳遞函數(shù)(Transfer Functions,TFs)比率的方法,用于衡量不同信道間的傳遞差異。首先,假設(shè)傳遞函數(shù)比所期望信號變化的更加緩慢,進(jìn)一步假設(shè)噪聲信號與期望信號的幅值比是緩慢變化,將時間軸劃分為一系列分析間隔,在每一個分析間隔期間內(nèi)假定TFs和噪聲信號均近似為平穩(wěn)。然后,將該分析時間間隔劃分為幀,使得期望信號在每幀內(nèi)看作不變。以第K幀的信號為例,得到:

(1)

(2)

(3)

(4)

由于Um(t,ejw)和Z1(t,ejw)一般是相關(guān)的,因此無法通過計算式(4)來直接獲得無偏估計。文獻(xiàn)[15]指出采用最小二乘法可求得超定方程組Hm(ejw)(m=2,3,…,M為一組獨立方程)的無偏估計。方程式(5)的解為方程式(6)。

(5)

(6)

其中:〈〉表示取均值運算;Hm(ejw)為傳遞函數(shù)比率的估計結(jié)果,可用于衡量兩個通道信號間的相關(guān)性,若兩個通道信號相關(guān)性強,則Hm(ejw)值較大;若兩個通道信號相關(guān)性弱,則Hm(ejw)值較小。

3 基于非平穩(wěn)系統(tǒng)辨識的心音分割

3.1 基于系統(tǒng)辨識的包絡(luò)提取

雖然心音信號近似為準(zhǔn)周期性信號,但其本質(zhì)上是一類非平穩(wěn)信號,具備混沌特性[16]。為較好地辨識出基礎(chǔ)心音信號,需要提取更有效的特征包絡(luò),突出基礎(chǔ)心音特征,因此本文提出一種基于非平穩(wěn)系統(tǒng)辨識(Non-Stationary System Identification,NSSI)理論的包絡(luò)提取方法。

心音信號具有短時平穩(wěn)性,因此將時間段劃分得足夠小時,相鄰兩個時間段內(nèi)的心音信號具有較大的相似性。而對噪聲和部分雜音而言,其非平穩(wěn)性特征更加明顯,在適當(dāng)長度的時間段內(nèi),其相似性更低。根據(jù)非平穩(wěn)系統(tǒng)辨識原理,將原始心音段作為參考通道信號,將后移一個時間段的信號作為另一個對比通道,通過非平穩(wěn)系統(tǒng)辨識來求解其相關(guān)性,從而提取得到心音信號的特征包絡(luò),具體步驟如下:

1)對采集的心音信號做分幀處理。在此過程中,幀長wlen的選擇較為關(guān)鍵,其關(guān)系著心音與噪聲包絡(luò)特征的區(qū)分程度。文獻(xiàn)[17]指出,當(dāng)心音信號為10 ms～30 ms時具有短時平穩(wěn)性,因此選取的幀長應(yīng)在此范圍內(nèi)。此外,幀移的選擇也需注意,幀移過小會加大計算量,幀移過大會使部分對應(yīng)幀間的相關(guān)度偏低,使得特征包絡(luò)不明顯,實際處理中通常取其為幀長的一半[2,10]。本文按式(6)測試了多種類型心音信號下的相關(guān)度,并給出部分結(jié)果,如圖2所示。由于噪聲及心音信號都具備非平穩(wěn)性,且噪聲具有隨機性,心音具有一定的隨機變異性,因此圖2中相關(guān)度值波動較大屬于正?，F(xiàn)象。在選擇幀長時,要選擇s1和s2相關(guān)度值差距較小且能明顯區(qū)分于噪聲的幀長。圖2中用矩形框標(biāo)記了各情況下的可取幀長,為方便起見,選取普適性更強的16 ms幀長。

圖2 不同病況與幀長下的相關(guān)度測試結(jié)果Fig.2 Test results of the correlation at different patient conditions and frame lengths

2)求取參考信號x1與對比信號x2。將分幀處理后的信號去掉最后一幀,得到參考信號x1;將分幀處理后的信號去掉第一幀,得到對比信號x2。

(7)

(8)

(9)

(10)

其中,xcorr表示求相關(guān),FFT表示快速傅里葉變換。

4)按照式(11)給出的非平穩(wěn)系統(tǒng)相關(guān)性辨識結(jié)果求得包絡(luò)特征Hx。

(11)

5)每幀的信號Hx值個數(shù)與幀長相同,因此需要選擇一個代表性的值來表示該幀的包絡(luò)特征值,本文按照常規(guī)方法取升序排列50%的值,即取中值。由篩選的中值組成一組特征值構(gòu)成提取包絡(luò)。特征包絡(luò)Hx能夠較好地提取基礎(chǔ)心音信號特征,區(qū)分出基礎(chǔ)心音、噪聲和雜音。由于特征包絡(luò)的整體呈脈沖狀,若不做平滑處理,則會出現(xiàn)錯誤分割現(xiàn)象,因此需要對包絡(luò)作進(jìn)一步處理。本文首先采用最值濾波來平滑特征包絡(luò)Hx,即取窗內(nèi)最大值為平滑值,平滑窗長取為4(根據(jù)幀長確定);然后采用鋸齒展寬處理,得到Hx_s特征包絡(luò)。

如圖3所示,本文將降噪后的心音信號圖與所提取的包絡(luò)圖進(jìn)行對比,并用矩形框分別標(biāo)記時域干擾音及其對應(yīng)的特征包絡(luò)。由圖3可以看出,提取的包絡(luò)能較好地體現(xiàn)基礎(chǔ)心音特征并湮沒噪聲,而香農(nóng)能量包絡(luò)法和希爾伯特變換包絡(luò)法卻無法湮沒難以抑制的噪聲,在多組心音信號測試中均存在此類現(xiàn)象。由于良好的特征包絡(luò)將有助于分割出基礎(chǔ)心音信號,因此包絡(luò)的有效性與普適性將通過分割結(jié)果進(jìn)一步驗證。

圖3 特征包絡(luò)提取結(jié)果Fig.3 Extraction results of feature envelopes

3.2 自適應(yīng)閾值及分割點判定

在提取出包絡(luò)的基礎(chǔ)上進(jìn)行閾值選取分析。為避免人工設(shè)置閾值的不便,本文采用提取相關(guān)閾值參數(shù)來自適應(yīng)分配閾值。一般地,閾值是有用信號與噪聲的分界線,因此所以首先選擇能夠反映信號中噪聲程度的參數(shù),其次閾值被用于檢測特征包絡(luò)是否滿足要求,所以閾值還應(yīng)與提取的特征包絡(luò)相關(guān)。此外,由于類似基礎(chǔ)心音信號的雜音干擾或受試者心音異常,所以分割結(jié)果還需進(jìn)一步修正,以保證分割結(jié)果的正確性。整個過程具體如下:

1)選擇反映噪聲程度的參數(shù)——信噪比(Signal to Noise Ratio,SNR),但由于實際中的噪聲或者純凈信號,而通常是不可預(yù)知的,因此本文采用降噪的方式來估計純凈信號,獲取SNR。在包絡(luò)提取時已經(jīng)進(jìn)行小波降噪,但并不能將噪聲完全濾除,因此進(jìn)一步加大濾波參數(shù)為cora80%來抑制噪聲估計純凈心音,按式(12)計算SNR[12]。SNR估計值將用于提取閾值T1。

(12)

其中,s(n)為估計的純凈信號,d(n)為噪聲,N為信號長度。

(13)

其中,NH為包絡(luò)長度。

3)按照選擇的參數(shù)建立閾值函數(shù)。在包絡(luò)分割中,常用的閾值函數(shù)形式為雙門限閾值,建立雙閾值函數(shù)如式(14)和式(15)所示:

(14)

(15)

若SNR為0,則:

(16)

其中,α和β為常數(shù),文獻(xiàn)[18]指出閾值系數(shù)可通過遺傳算法進(jìn)行全局優(yōu)化,實驗得到α=1.4、β=1.6為最佳值。

4)根據(jù)閾值尋找分割點。分割點分為上分割點與下分割點[19],即沿時間軸取包絡(luò)與閾值T1的交點ui和包絡(luò)與閾值T2的交點rl的中點為分割點,當(dāng)ui>rl時為下分割點,當(dāng)ui

5)計算每組分割點間包絡(luò)所對應(yīng)的面積,對分割點進(jìn)行第一次篩檢[5,9]。由于所提取的基礎(chǔ)心音信號分割點間包絡(luò)所圍的面積要遠(yuǎn)大于噪聲信號分割點間的面積,因此可舍去所圍面積過小的分割點,排除噪聲包絡(luò)的干擾。

6)通過時域檢查波峰對分割點進(jìn)行第二次篩檢。提取每組分割點間的最高波峰,沿時間軸可構(gòu)成一組波峰時間點{pi|p1,)p2,…,pm},若波峰間的時間差滿足Δpi=pi+1-pi<0.1 s,則認(rèn)為pi+1或pi所對應(yīng)的分割點為錯誤分割點。如果采集的受試者心音存在類似基礎(chǔ)心音的心雜音、第三心音、第四心音,或者是奔馬律等異常類型心音,又或者是首位存在不完整的基礎(chǔ)心音信號,都極有可能導(dǎo)致分割點出錯。而由于目前的包絡(luò)分割方法多數(shù)未考慮此類情況,因此在實際應(yīng)用中容易加大分割誤差。本文根據(jù)心音信號的時域特征,提出相應(yīng)的誤分割點檢驗方法,具體步驟如下:

(1)識別不完整的基礎(chǔ)心音,去除單個分割點。由于用于分割的心音信號開始處和結(jié)束處可能存在不完整的基礎(chǔ)心音,會造成分割點判定結(jié)果存在單個的上分割點或下分割點,單個分割點會影響υi和νl匹配,同時會影響起始基礎(chǔ)心音信號的識別,因此需要檢測出單個的υi和νl,并將其從檢測結(jié)果中去除。在不完整的基礎(chǔ)心音檢測中,將第一個上分割點與第一個下分割點進(jìn)行比較,若第一個上分割點位置在第一個下分割點位置的后面,則判定ν1為不完整的下分割點,并利用類似方法檢測末尾的分割點,整個非完整基礎(chǔ)心音檢測及濾除過程可表示為:

(17)

其中,Qdel表示需要刪除的分割點,υend表示最后一個上分割點,νend表示最后一個下分割點。

(2)檢測心雜音或異常額外心音的分割點。實際中采集到的很多心雜音或異常額外心音與基礎(chǔ)心音非常類似,這些非基礎(chǔ)心音成分會對分割點的判定造成干擾。文獻(xiàn)[20]指出收縮期時長一般為160 ms～240 ms,在檢測誤檢分割點時,取收縮期時長最小值的一半(80 ms),由于基礎(chǔ)心音信號包絡(luò)值較非基礎(chǔ)心音信號包絡(luò)值更大,因此可采用如下判別方法:

(18)

當(dāng)誤檢點出現(xiàn)在舒張期的中間位置時,上述檢測方法難以有效判定,因此需要進(jìn)一步檢測。文獻(xiàn)[20]指出心動周期時長一般為600 ms～1 000 ms,取心動周期最短時長作為又一判定參量,采用如下檢測方法:

(19)

7)提取時域特征,識別心音成分。舒張期持續(xù)時長一般大于收縮期時長[18-19],故根據(jù)時域時長判定持續(xù)時間較長的成分為舒張期,較短的成分為收縮期[20-21],按順序處于舒張期與收縮期之間的成分為s1,則剩余成分為s2。

3.3 分割流程及評價指標(biāo)

根據(jù)上文所述方法完成包絡(luò)提取及分割任務(wù),具體流程如圖4所示。

圖4 自適應(yīng)閾值分割流程Fig.4 Procedure of adaptive threshold segmentation

由于本文方法主要面向復(fù)雜心音信號實時分割處理,難以通過直觀的分割定位結(jié)果對分割算法進(jìn)行有效性檢驗,因此本文采用文獻(xiàn)[22]提出的4類評價指標(biāo)作為檢驗指標(biāo):

(20)

(21)

(22)

PS=run_time

(23)

其中,TPFHS為基礎(chǔ)心音信號檢出正確率,FPFHS為非基礎(chǔ)心音信號檢入錯誤率,MPFHS為基礎(chǔ)心音未檢出率,TNFHS為正確檢入的基礎(chǔ)心音信號幀數(shù),FNFHS為錯誤檢入的非基礎(chǔ)心音信號幀數(shù),MNFHS為未檢入的基礎(chǔ)心音信號幀數(shù),PS為算法運行時間。TPFHS越高,FPFHS和TNFHS越低,代表分割精度越高;PS越小,代表實時性越好。

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

4.1 實驗數(shù)據(jù)及其處理

實驗數(shù)據(jù)由中國醫(yī)學(xué)科學(xué)院阜外醫(yī)院和南京郵電大學(xué)電子科學(xué)與工程學(xué)院提供。其中選用前面提供心音62條,后者提供心音43條,所有心音信號均降采樣至2 000 Hz。為模擬復(fù)雜環(huán)境,選用noiseX-92噪聲數(shù)據(jù)庫中的15類不同領(lǐng)域的噪聲按5 dB、10 dB、15 dB這3種信噪比加入心音信號中構(gòu)成待處理復(fù)雜心音,噪聲疊加不超過兩種,選用的心音信號中還包括17類異常心音信號。noiseX-92噪聲數(shù)據(jù)庫中的噪聲類型具體為白噪聲、粉紅噪聲、高頻通道噪聲、飛機駕駛艙噪聲1、飛機駕駛艙噪聲2、廠房噪聲1、廠房噪聲2、類語音噪聲、坦克噪聲、驅(qū)逐艦操作室噪音、驅(qū)逐艦機艙噪音、F16機艙噪聲、軍用車輛噪聲、機槍噪音和車內(nèi)噪聲。心音信號類型及數(shù)量如表1所示。

表1 心音信號類型及數(shù)量Table 1 Type and quantity of heart sound signal

105條心音信號經(jīng)處理后生成33 075條帶噪心音信號,共計516 285個心動周期,原始心音信號均在Cool Edit Pro下進(jìn)行手工標(biāo)注分割點。實驗軟件為Matlab2017a,計算機CPU為Intel?CoreTMi7-8700k@3.70 GHz,RAM為16.0 GB,硬盤為1 TB。

4.2 結(jié)果分析

本文方法分割結(jié)果如圖5所示,直觀反映了其能夠較精準(zhǔn)地分割出基礎(chǔ)心音信號,且不會將非基礎(chǔ)心音誤分割為基礎(chǔ)心音信號。本文方法與維奧拉積分包絡(luò)分割法[5]、改進(jìn)型希爾伯特-黃變換包絡(luò)雙閾值分割法[7]、基于心動周期估計的心音包絡(luò)分割法[9]、基于周期與波峰的自適應(yīng)閾值包絡(luò)分割法[10]和基于個性化高斯混合建模的包絡(luò)分割法[4]進(jìn)行對比,其中基于個性化高斯混合建模的包絡(luò)分割法采用在線訓(xùn)練方式,無需通過不同的訓(xùn)練集進(jìn)行訓(xùn)練,能夠避免對比差異。

圖5 基于本文方法的心音信號分割結(jié)果Fig.5 Segmentation results of heart sound signal based on the proposed method

表2給出了6種分割方法的評價指標(biāo)結(jié)果,可以看出,基于個性化高斯混合建模的分割法檢出正確率TPFHS最高為91.49%,其次為本文方法的89.21%,其余方法的檢出準(zhǔn)確率相對較低,這是因為本文方法可避免雜音和誤檢干擾。此外,可以發(fā)現(xiàn)其他4種基于包絡(luò)特征的分割方法中基礎(chǔ)心音錯誤檢入率FPFHS要高于未檢入率MPFHS,主要原因為在基礎(chǔ)心音末尾處噪聲存在拖尾,容易使非基礎(chǔ)心音的噪聲部分被誤判為基礎(chǔ)心音幀,而本文方法不存在此類問題。

表2 6種方法的分割性能比較Table 2 Comparison of segmentation performance of six methods

另外,在算法耗時PS指標(biāo)評測中,采用的心音信號段為4.5 s,即每5個心動周期為一段,獲取PS的平均值為最終結(jié)果。從表2可以看出,基于個性化高斯混合建模的包絡(luò)分割法的耗時為6.573 5 s,遠(yuǎn)高于其他分割方法,這是因為基于個性化高斯混合建模的包絡(luò)分割法采用在線訓(xùn)練方式,涉及到高斯混合建模、前向后向運算、Viterbi解碼、誤差修正等計算[4],復(fù)雜度遠(yuǎn)高于其他基于包絡(luò)的雙閾值分割方法;基于改進(jìn)型希爾伯特-黃變換包絡(luò)的分割法先進(jìn)行經(jīng)驗?zāi)B(tài)分解,然后做3次貝塞爾樣條插值運算,同時對端點效應(yīng)進(jìn)行優(yōu)化[7],其運算復(fù)雜度次于基于個性化高斯混合建模的包絡(luò)分割法;維奧拉積分包絡(luò)分割法在確定尺度與基礎(chǔ)心音的關(guān)系后直接提取包絡(luò)進(jìn)行分割[5],基于周期與波峰的自適應(yīng)閾值包絡(luò)分割法建立自相關(guān)估計周期與波峰的關(guān)系閾值后完成分割[10],本文方法在提取心音包絡(luò)后直接進(jìn)行分割與判別,復(fù)雜度均遠(yuǎn)低于前兩種方法。本文方法的耗時為0.105 2 s,能夠滿足實時分割要求,可為實時監(jiān)測的心音自動分析技術(shù)提供參考。

5 結(jié)束語

基于包絡(luò)的心音信號分割方法因算法簡單、時效性強等特點而得到廣泛應(yīng)用,其中的特征包絡(luò)提取與閾值處理算法設(shè)計為關(guān)鍵步驟。本文提出一種基于非平穩(wěn)系統(tǒng)辨識的心音特征包絡(luò)提取方法,采用遺傳算法設(shè)計自適應(yīng)雙閾值函數(shù)實現(xiàn)閾值分割。在近似基礎(chǔ)心音的心雜音處理上,根據(jù)時域特征進(jìn)行分割點篩查,識別并去除誤檢分割點和單獨分割點,從而得到正確的基礎(chǔ)心音分割點。實驗結(jié)果表明,本文方法相比維奧拉積分包絡(luò)分割法、改進(jìn)型希爾伯特-黃變換包絡(luò)雙閾值分割法、基于心動周期估計的心音包絡(luò)分割法和基于周期與波峰的自適應(yīng)閾值分割法分割精度更高,相比基于個性化高斯混合建模的包絡(luò)分割法實時性更強,提取的包絡(luò)能夠更好地反映基礎(chǔ)心音信號特征。但本文方法對于人體運動狀態(tài)下的心音信號分割精度不理想,后續(xù)將對此做進(jìn)一步研究。