許春冬,周 靜,應冬文,2,龍清華

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

0 概述

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

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

1 基于包絡的分割方法

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

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

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

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

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

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

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

(1)

(2)

(3)

(4)

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

(5)

(6)

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

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

3.1 基于系統辨識的包絡提取

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

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

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

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

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

(7)

(8)

(9)

(10)

其中,xcorr表示求相關,FFT表示快速傅里葉變換。

4)按照式(11)給出的非平穩(wěn)系統相關性辨識結果求得包絡特征Hx。

(11)

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

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

圖3 特征包絡提取結果Fig.3 Extraction results of feature envelopes

3.2 自適應閾值及分割點判定

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

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

(12)

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

(13)

其中,NH為包絡長度。

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

(14)

(15)

若SNR為0,則:

(16)

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

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

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

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

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

(17)

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

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

(18)

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

(19)

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

3.3 分割流程及評價指標

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

圖4 自適應閾值分割流程Fig.4 Procedure of adaptive threshold segmentation

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

(20)

(21)

(22)

PS=run_time

(23)

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

4 實驗結果與分析

4.1 實驗數據及其處理

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

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

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

4.2 結果分析

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

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

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

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

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

5 結束語

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