賀光碩,盧國梁*,尚偉
(1.山東大學(xué)機械工程學(xué)院,濟南 250061;2.山東大學(xué)高效潔凈機械制造教育部重點實驗室,濟南 250061;3.山東大學(xué)機械工程國家級實驗教學(xué)示范中心,濟南 250061;4.山東大學(xué)第二醫(yī)院,濟南 250033)
腦電圖(electroencephalogram,EEG)包含豐富的生理和疾病信息,在臨床醫(yī)學(xué)領(lǐng)域具有廣泛的應(yīng)用[1]。EEG信號的分析和處理能夠為某些神經(jīng)性疾病如癲癇[2]、急性腦血管病[3]、腦梗死[4]等疾病的早期預(yù)警和診斷提供相關(guān)依據(jù)。其主要技術(shù)手段是對腦電信號的狀態(tài)變化進(jìn)行監(jiān)測,即實時檢測腦電信號由正常狀態(tài)變?yōu)楫惓5倪^程[5]。
腦電信號動態(tài)變化的異常檢測是近幾十年來在腦電信號處理領(lǐng)域的重要研究課題,其流程主要包括:信號預(yù)處理、動態(tài)建模、決策制定[6]。其中關(guān)鍵的一步是對腦電信號進(jìn)行動態(tài)建模,并提取一個能夠表征腦電信號動態(tài)特性的綜合性指標(biāo)。然而,由其他生理信號和外界噪聲引起腦電信號的非平穩(wěn)、非線性和高變性等特點極大地影響了綜合性指標(biāo)的提取過程[7]。因此傳統(tǒng)的線性指標(biāo)如均值(mean)、方差(variance)、均方根(root mean square,RMS)和峰度(kurtosis)等極易受腦電信號高度復(fù)雜性的影響而產(chǎn)生各種誤報[8]。隨著腦電信號處理方法的研究,非線性指標(biāo)的發(fā)展克服了傳統(tǒng)指標(biāo)的不足,受到研究者的青睞[9-10]。圖模型(graph model,GM)作為一種譜分析算法,是一種獨特的非線性指標(biāo),能夠捕捉數(shù)據(jù)中隱藏的結(jié)構(gòu)和拓?fù)湫畔ⅲ驯挥糜谀X電信號分析[11]。Zhang等[12]利用圖模型的結(jié)構(gòu)信息來表示腦電傳感器的空間特征,然后進(jìn)行運動意圖檢測。Song等[13]使用圖模型建模并提取多通道腦電信號的特征用來進(jìn)行腦電信號的情緒識別。Haddad等[14]根據(jù)圖建模描述了每個病人的發(fā)作特征,從而預(yù)測癲癇發(fā)作。
雖然動態(tài)建模方法發(fā)展迅速,但一個不可回避且尚未完全解決的問題是如何檢測給定腦電信號中的微弱的異常變化,目前的方法在腦電信號微弱異常變化檢測方面效果不佳。檢測微弱變化能夠提前預(yù)測和診斷癲癇發(fā)作等神經(jīng)疾病,從而為醫(yī)生的治療做好提前準(zhǔn)備,在臨床應(yīng)用中具有重大意義[10]。
在腦電信號的處理過程中,通常把腦電信號發(fā)生異常變化的頻段分為δ、θ、α、β和 γ波,不同的頻率段包含不同的疾病信息和原始信號信息[15]。雖然現(xiàn)有的GM動態(tài)建模方法對噪聲具有很強的魯棒性,但在處理腦電信號的微弱異常變化時仍有很大的局限性。
針對以上問題,現(xiàn)提出一種基于多頻段圖模型的腦電信號微弱異常檢測方法。采用GM對δ、θ、α、β和 γ 5個頻段分別進(jìn)行建模來增強的腦電微弱變化信息,并提取各個頻段腦電信號的特征,采用自適應(yīng)權(quán)重融合算法融合頻段特征,生成綜合性指標(biāo)用來描述腦電信號的動態(tài)變化,以期解決腦電信號中微弱異常變化的檢測準(zhǔn)確率低的問題。
所用數(shù)據(jù)庫選自公開的CHB-MIT頭皮EEG數(shù)據(jù)庫[16](網(wǎng)址:https://physionet.org/content /chbmit/1.0.0/)和來自山東大學(xué)第二醫(yī)院神經(jīng)內(nèi)科的EEG數(shù)據(jù)庫[11,17]。
1.1.1 CHB-MIT頭皮EEG數(shù)據(jù)庫
該數(shù)據(jù)庫由患有難治性癲癇病的兒科受試者的腦電圖數(shù)據(jù)記錄組成,收集了22名受試者共24個案例,采樣率為256 Hz。每個案例(chb01、chb02等)包含來自一名受試者的9~42個連續(xù).edf文件并且大多數(shù)文件的時間跨度為1 h。如圖1所示,每個案例詳細(xì)地標(biāo)注了癲癇發(fā)作的時間和次數(shù),方便了實驗時對異常變化點的標(biāo)注。在每種案例中,隨機選擇一個含有癲癇發(fā)作的.edf文件作為測試數(shù)據(jù),共24個數(shù)據(jù),選擇第16個通道數(shù)據(jù)作為實驗所需的測試數(shù)據(jù)。
圖1 原始腦電信號癲癇發(fā)作時間(以chb01_03為例)Fig.1 Raw EEG signal seizure time(take chb01_03 as an example)
1.1.2 山東大學(xué)第二醫(yī)院神經(jīng)內(nèi)科EEG數(shù)據(jù)庫
該數(shù)據(jù)庫的所有EEG數(shù)據(jù)均來自癲癇患者,且所有受試者均已獲得書面知情同意。該實驗完全遵守《赫爾辛基宣言》和當(dāng)?shù)叵嚓P(guān)法律法規(guī)。實驗所用的癲癇檢測系統(tǒng)用于監(jiān)測夜間睡眠中的癲癇[17],同時收集受試者的EEG數(shù)據(jù)和相應(yīng)的監(jiān)控視頻。
原始腦電信號以1 024 Hz的采樣率記錄,并在進(jìn)一步分析之前降采樣為512 Hz。該數(shù)據(jù)集包含32通道的腦電信號,隨機選取其中的C3通道為測試數(shù)據(jù),共包含55個數(shù)據(jù)序列。由于癲癇發(fā)作時腦電信號的變化非常微弱,很難用肉眼直接判斷是否變化。因此,具有豐富臨床經(jīng)驗的神經(jīng)科醫(yī)生同時觀看監(jiān)控視頻和所有通道的數(shù)據(jù),來標(biāo)記EEG的異常變化位置。
癲癇等神經(jīng)源性疾病的頻率成分主要出現(xiàn)在δ、θ、α、β和 γ頻段[15,18]。在神經(jīng)性疾病發(fā)作時,這些頻率的幅值會相應(yīng)改變。因此,根據(jù)以上5個頻段的頻率定義一個0.1~70 Hz的帶通濾波器來捕獲該頻率的腦電信號。此外,使用特定的陷波濾波器來消除50 Hz的電力線的干擾。
根據(jù)帶通濾波器的截止頻率和腦電信號實際異常變化情況,腦電信號各個波段頻率范圍為:δ波(0.1~4 Hz)、θ波(4~8 Hz)、α波(8~12 Hz)、β波(12~30 Hz)、γ波(31~70 Hz)。對于一段給定的濾波后的腦電信號,分別進(jìn)行5個頻段的提取,如圖2所示??梢钥闯觯X電信號的5個頻段在異常變化檢測中具有不同表現(xiàn),如α、β和 γ頻段在異常變化位置反應(yīng)明顯,而δ和θ頻段則無明顯變化,因此所提算法對于腦電信號的微弱異常變化檢測具有重要意義。
圖2 原始腦電信號5個頻段的濾波信號Fig.2 Filtered signals of the five frequency bands of the raw EEG signal
基于多頻段圖模型方法,提出了腦電信號異常檢測框架,步驟如下。
步驟1構(gòu)造圖模型:對腦電信號的5個頻段(即δ、θ、α、β和 γ波)分別進(jìn)行GM動態(tài)建模。
步驟2相似性分?jǐn)?shù):對每個頻段的GM 運用距離函數(shù)得到能夠量化GM之間關(guān)系的相似性分?jǐn)?shù)。
步驟3自適應(yīng)權(quán)重融合算法:利用自適應(yīng)權(quán)重融合算法融合所有頻段的相似性分?jǐn)?shù)并最終得到一項綜合性指標(biāo)用來描述腦電信號的動態(tài)變化。
步驟4假設(shè)檢驗:最后采用假設(shè)檢驗來檢測腦電信號是否存在微弱的異常變化。
針對一段數(shù)據(jù),若檢測到異常變化,該框架將向用戶報告一個異常變化警報,然后繼續(xù)一個新的檢測過程。若沒有檢測到變化,該框架就會一直運行下去,直到數(shù)據(jù)運行結(jié)束。
對5個頻段中每個頻段分別進(jìn)行構(gòu)造圖模型,在1.2節(jié)收集的濾波后的腦電信號中,假定δ波的腦電信號為Xδ={Xδ(t)},t=1,2,…,N,其中Xδ(t)為t時刻的觀測值。
構(gòu)造圖模型的過程如下。
步驟1利用離散短時傅里葉變換(discrete short time Fourier transform,DSTFT)將信號從時域轉(zhuǎn)換到頻域。當(dāng)窗口長度固定時,DSTFT的計算公式為
(1)
式(1)中:φ*(·)為窗口函數(shù);p為DSTFT得到的頻率數(shù)據(jù);i=1,2,…,N/2;j為虛數(shù);N為窗口函數(shù)移動的時間步長,N=512;m為第m個時間步長。
步驟2計算功率譜,通過式(2)可以得到一系列周期圖P。
(2)
步驟3周期圖P(i,m)用于定義時間t處的頻率分辨率,即Fm={f(i)},i=1,2,…,N/2。其中,f(i)為頻率,F(xiàn)m為頻率f(i)的集合。將每個頻率分辨率視為一個節(jié)點,計算每兩個節(jié)點之間的歐氏距離di,j作為邊的權(quán)重。
(3)
(4)
對圖2中信號的5個頻段進(jìn)行以上步驟分析,得到相似性分?jǐn)?shù)如圖3所示??梢钥闯?,θ波的相似性分?jǐn)?shù)在異常變化位置沒有明顯變化,其他波段相似性分?jǐn)?shù)變化較為明顯,主要原因是腦電信號的變化信息隱藏在不同的頻率成分中。因此,為了充分利用這些重要的微弱變化信息來實現(xiàn)對EEG信號的自適應(yīng)檢測,需要一種自適應(yīng)權(quán)重融合算法對所有的相似性分?jǐn)?shù)進(jìn)行加權(quán)融合。
圖3 相似性分?jǐn)?shù)計算示例Fig.3 Example of calculation of the similarity scores
自適應(yīng)權(quán)重通過依據(jù)各個分量的重要性自適應(yīng)地計算權(quán)重,將所有相似性分?jǐn)?shù)組合在一起。計算過程如下。
(5)
步驟2計算方差,其計算公式為
(6)
步驟3使用方差自適應(yīng)地計算每個分量的權(quán)重wk,其計算公式為
(7)
步驟4所有頻段的相似性分?jǐn)?shù)可表示為綜合性指標(biāo)Qi,其計算公式為
(8)
基于綜合性指標(biāo){Q1,Q2,…,Qm,…,QM},采用常零假設(shè)檢驗進(jìn)行變化決策。對于第m段,假設(shè)檢驗如下:H0:沒發(fā)生變化,Qm∈??;H1:發(fā)生變化,Qm?Α,其中,Α=[μm-1-3σm-1,μm-1+3σm-1]為由常用的3σ準(zhǔn)則定義的置信區(qū)間,其中μm-1和σm-1分別為前m-1個片段綜合性指標(biāo)的平均值和標(biāo)準(zhǔn)差。
圖4給出了檢測結(jié)果的示例??梢钥闯?,使用的假設(shè)檢驗可以成功地檢測出異常發(fā)生的時間,證明了所提出的綜合性指標(biāo)的有效性。
圓圈表示測試方法產(chǎn)生的異常變化點圖4 檢測結(jié)果示例Fig.4 Example of detection result
兩個不同數(shù)據(jù)集的實驗結(jié)果,與傳統(tǒng)線性指標(biāo)如均值(mean)、方差(variance)、均方根(root mean square,RMS)和峰度(kurtosis)以及典型的非線性指標(biāo)經(jīng)驗?zāi)B(tài)分解(empirical mode decomposition,EMD)、短時傅里葉變換(short-time fourier transform,STFT)、連續(xù)小波變換(continuous wavelet transform,CWT)和離散小波變換(discrete wavelet transform,DWT)進(jìn)行了比較。采用常用指標(biāo)查準(zhǔn)率Precision、查全率Recall和F分?jǐn)?shù)Fscore來評價以上方法的檢測性能,定義為
(9)
(10)
(11)
式中:TP為真正例,表示正確的檢測結(jié)果;FP為假正例,表示錯誤的檢測結(jié)果;FN為假負(fù)例,表示沒有檢測結(jié)果。
顯然,高查準(zhǔn)率表示高檢測精度,高查全率表示高檢測靈敏度,F(xiàn)分?jǐn)?shù)對檢測的方法進(jìn)行綜合評價。
由于CHB-MIT頭皮EEG數(shù)據(jù)庫腦電信號的數(shù)據(jù)量過于龐大,只截取癲癇發(fā)作開始時間點的前60 s及后60 s作為實驗數(shù)據(jù)。圖5展示了所提方法用于測試CHB-MIT頭皮EEG數(shù)據(jù)庫實驗結(jié)果的示例。
顯然,所提方法綜合性指標(biāo)能夠及時檢測到腦電信號的微弱變化。由圖5可知,DWT有誤報點,EMD和STFT未能成功檢測到異常變化點,其余方法均表現(xiàn)良好的檢測結(jié)果,但是并不能通過此例證明其他方法比綜合性指標(biāo)方法好。如表1所示,綜合性指標(biāo)在查準(zhǔn)率、查全率和F分?jǐn)?shù)之間取得的效果最好,表現(xiàn)出了該方法的優(yōu)越性。盡管方差在查全率方面效果為100%,但是它的查準(zhǔn)率和F分?jǐn)?shù)沒有綜合性指標(biāo)的性能要好,這也并不能說明方差在檢測性能上好于綜合性指標(biāo)。
表1 CHB-MIT數(shù)據(jù)庫檢測結(jié)果進(jìn)行了3個指標(biāo)的比較Table 1 Comprehensive comparison in CHB-MIT scalp EEG database with three indicators
圓圈表示測試方法產(chǎn)生的異常變化點圖5 綜合性指標(biāo)與其他方法的比較結(jié)果Fig.5 Comparison results of comprehensive indicators and other methods
圖6展示了所提方法用于測試山東大學(xué)第二醫(yī)院神經(jīng)內(nèi)科EEG數(shù)據(jù)庫實驗結(jié)果的兩個示例。本實驗中測試數(shù)據(jù)的變化類型主要是幅值變化和頻率變化。圖6中測試數(shù)據(jù)1的異常變化主要是由幅值變化引起的。可以看出,除了均值、峰度和經(jīng)驗?zāi)B(tài)分解方法產(chǎn)生了誤報點,其余方法都能及時檢測到腦電數(shù)據(jù)的微弱異常變化。
檢測由頻率變化引起的腦電數(shù)據(jù)異常更具有挑戰(zhàn)性。如圖6(b)所示,很難從原始腦電信號中用肉眼直觀地識別腦電信號的微弱異常變化。經(jīng)過實驗,只有綜合性指標(biāo)能夠及時檢測到腦電信號微弱的異常變化,其他方法均產(chǎn)生延遲或誤報點。表2中的綜合比較結(jié)果驗證了這一點,只有綜合性指標(biāo)在查準(zhǔn)率、查全率和F分?jǐn)?shù)之間取得了最好的效果,并且查全率的結(jié)果為100%,這意味著該方法能夠檢測出所有潛在的腦電信號異常變化。腦電信號的非線性、非平穩(wěn)性和高變性使其極容易受外界影響,極大地增加了信號的異常檢測難度,所以綜合性指標(biāo)在查準(zhǔn)率方面未能達(dá)到百分之百也是有情可原?;谝陨辖Y(jié)果,綜合性指標(biāo)具有很高的檢測精度和靈敏度,并且能夠準(zhǔn)確反映腦電信號的動態(tài)變化。
圓圈表示測試方法產(chǎn)生的異常變化點圖6 山東大學(xué)數(shù)據(jù)集實驗結(jié)果比較Fig.6 Comparison of experimental results of Shandong University data set
表2 山東大學(xué)數(shù)據(jù)集檢測結(jié)果中三個指標(biāo)的比較Table 2 Comprehensive comparison in Shandong University data set with three indicators
在設(shè)計一個實用的變化檢測系統(tǒng)時,計算效率是一個需要考慮的重要問題。綜合性指標(biāo)與其他比較方法的檢測步驟基本一樣,均包括數(shù)據(jù)建模、指標(biāo)提取和決策3個階段。表3給出了綜合性指標(biāo)與其他方法的計算復(fù)雜度的比較結(jié)果。結(jié)果表明,在所有的方法中,決策的計算復(fù)雜度是相同的,但在其他兩個階段的計算復(fù)雜度卻有很大的不同。顯然,綜合性指標(biāo)的總計算復(fù)雜度是最大的,這意味著隨著腦電數(shù)據(jù)量的增加,該方法所消耗的時間也會大量增加。因此,在將來的工作中會進(jìn)一步提高綜合性指標(biāo)的計算效率,增加其實用價值。
表3 計算復(fù)雜度的比較Table 3 Comparison of the computational complexity
提出了一種基于多頻段圖模型的腦電信號微弱異常變化檢測新方法。該方法充分利用這五個頻段中腦電的微弱變化信息和GM的優(yōu)勢,對腦電信號δ、θ、α、β和 γ波的每個頻段單獨進(jìn)行圖模型建模,并通過自適應(yīng)權(quán)重融合算法融合各圖模型的相似性分?jǐn)?shù)形成綜合性指標(biāo)進(jìn)行假設(shè)檢驗。兩個實驗結(jié)果表明,所提方法查全率結(jié)果均為100%,表明所提方法能夠檢測所有潛在的腦電信號異常變化,因此該結(jié)果對腦電信號的微弱異常變化的檢測及實際應(yīng)用具有重要意義。由于所提方法的計算復(fù)雜度較高,下一步會對所提方法進(jìn)一步改進(jìn)以提高其計算效率,增加其應(yīng)用價值。