基于多頻段圖模型的微弱腦電信號異常檢測

賀光碩，盧國梁*，尚偉

(1.山東大學機械工程學院，濟南 250061;2.山東大學高效潔凈機械制造教育部重點實驗室，濟南 250061；3.山東大學機械工程國家級實驗教學示范中心，濟南 250061；4.山東大學第二醫(yī)院，濟南 250033)

腦電圖(electroencephalogram，EEG)包含豐富的生理和疾病信息，在臨床醫(yī)學領域具有廣泛的應用[1]。EEG信號的分析和處理能夠為某些神經(jīng)性疾病如癲癇[2]、急性腦血管病[3]、腦梗死[4]等疾病的早期預警和診斷提供相關依據(jù)。其主要技術手段是對腦電信號的狀態(tài)變化進行監(jiān)測，即實時檢測腦電信號由正常狀態(tài)變?yōu)楫惓５倪^程[5]。

腦電信號動態(tài)變化的異常檢測是近幾十年來在腦電信號處理領域的重要研究課題，其流程主要包括：信號預處理、動態(tài)建模、決策制定[6]。其中關鍵的一步是對腦電信號進行動態(tài)建模，并提取一個能夠表征腦電信號動態(tài)特性的綜合性指標。然而，由其他生理信號和外界噪聲引起腦電信號的非平穩(wěn)、非線性和高變性等特點極大地影響了綜合性指標的提取過程[7]。因此傳統(tǒng)的線性指標如均值(mean)、方差(variance)、均方根(root mean square,RMS)和峰度(kurtosis)等極易受腦電信號高度復雜性的影響而產(chǎn)生各種誤報[8]。隨著腦電信號處理方法的研究，非線性指標的發(fā)展克服了傳統(tǒng)指標的不足，受到研究者的青睞[9-10]。圖模型(graph model，GM)作為一種譜分析算法，是一種獨特的非線性指標，能夠捕捉數(shù)據(jù)中隱藏的結構和拓撲信息，已被用于腦電信號分析[11]。Zhang等[12]利用圖模型的結構信息來表示腦電傳感器的空間特征，然后進行運動意圖檢測。Song等[13]使用圖模型建模并提取多通道腦電信號的特征用來進行腦電信號的情緒識別。Haddad等[14]根據(jù)圖建模描述了每個病人的發(fā)作特征，從而預測癲癇發(fā)作。

雖然動態(tài)建模方法發(fā)展迅速，但一個不可回避且尚未完全解決的問題是如何檢測給定腦電信號中的微弱的異常變化，目前的方法在腦電信號微弱異常變化檢測方面效果不佳。檢測微弱變化能夠提前預測和診斷癲癇發(fā)作等神經(jīng)疾病，從而為醫(yī)生的治療做好提前準備，在臨床應用中具有重大意義[10]。

在腦電信號的處理過程中，通常把腦電信號發(fā)生異常變化的頻段分為δ、θ、α、β和 γ波，不同的頻率段包含不同的疾病信息和原始信號信息[15]。雖然現(xiàn)有的GM動態(tài)建模方法對噪聲具有很強的魯棒性，但在處理腦電信號的微弱異常變化時仍有很大的局限性。

針對以上問題，現(xiàn)提出一種基于多頻段圖模型的腦電信號微弱異常檢測方法。采用GM對δ、θ、α、β和 γ 5個頻段分別進行建模來增強的腦電微弱變化信息，并提取各個頻段腦電信號的特征，采用自適應權重融合算法融合頻段特征，生成綜合性指標用來描述腦電信號的動態(tài)變化，以期解決腦電信號中微弱異常變化的檢測準確率低的問題。

1 實驗數(shù)據(jù)來源

1.1 腦電數(shù)據(jù)庫的選取

所用數(shù)據(jù)庫選自公開的CHB-MIT頭皮EEG數(shù)據(jù)庫[16](網(wǎng)址：https://physionet.org/content /chbmit/1.0.0/)和來自山東大學第二醫(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所示，每個案例詳細地標注了癲癇發(fā)作的時間和次數(shù)，方便了實驗時對異常變化點的標注。在每種案例中，隨機選擇一個含有癲癇發(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 山東大學第二醫(yī)院神經(jīng)內(nèi)科EEG數(shù)據(jù)庫

該數(shù)據(jù)庫的所有EEG數(shù)據(jù)均來自癲癇患者，且所有受試者均已獲得書面知情同意。該實驗完全遵守《赫爾辛基宣言》和當?shù)叵嚓P法律法規(guī)。實驗所用的癲癇檢測系統(tǒng)用于監(jiān)測夜間睡眠中的癲癇[17]，同時收集受試者的EEG數(shù)據(jù)和相應的監(jiān)控視頻。

原始腦電信號以1 024 Hz的采樣率記錄，并在進一步分析之前降采樣為512 Hz。該數(shù)據(jù)集包含32通道的腦電信號，隨機選取其中的C3通道為測試數(shù)據(jù)，共包含55個數(shù)據(jù)序列。由于癲癇發(fā)作時腦電信號的變化非常微弱，很難用肉眼直接判斷是否變化。因此，具有豐富臨床經(jīng)驗的神經(jīng)科醫(yī)生同時觀看監(jiān)控視頻和所有通道的數(shù)據(jù)，來標記EEG的異常變化位置。

1.2 數(shù)據(jù)的預處理

癲癇等神經(jīng)源性疾病的頻率成分主要出現(xiàn)在δ、θ、α、β和 γ頻段[15,18]。在神經(jīng)性疾病發(fā)作時，這些頻率的幅值會相應改變。因此，根據(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)。對于一段給定的濾波后的腦電信號，分別進行5個頻段的提取，如圖2所示?？梢钥闯?，腦電信號的5個頻段在異常變化檢測中具有不同表現(xiàn)，如α、β和 γ頻段在異常變化位置反應明顯，而δ和θ頻段則無明顯變化，因此所提算法對于腦電信號的微弱異常變化檢測具有重要意義。

圖2 原始腦電信號5個頻段的濾波信號Fig.2 Filtered signals of the five frequency bands of the raw EEG signal

2 實驗方法

基于多頻段圖模型方法，提出了腦電信號異常檢測框架，步驟如下。

步驟1構造圖模型：對腦電信號的5個頻段(即δ、θ、α、β和 γ波)分別進行GM動態(tài)建模。

步驟2相似性分數(shù)：對每個頻段的GM 運用距離函數(shù)得到能夠量化GM之間關系的相似性分數(shù)。

步驟3自適應權重融合算法：利用自適應權重融合算法融合所有頻段的相似性分數(shù)并最終得到一項綜合性指標用來描述腦電信號的動態(tài)變化。

步驟4假設檢驗：最后采用假設檢驗來檢測腦電信號是否存在微弱的異常變化。

針對一段數(shù)據(jù)，若檢測到異常變化，該框架將向用戶報告一個異常變化警報，然后繼續(xù)一個新的檢測過程。若沒有檢測到變化，該框架就會一直運行下去，直到數(shù)據(jù)運行結束。

2.1 構造圖模型

對5個頻段中每個頻段分別進行構造圖模型，在1.2節(jié)收集的濾波后的腦電信號中，假定δ波的腦電信號為Xδ={Xδ(t)},t=1,2,…,N，其中Xδ(t)為t時刻的觀測值。

構造圖模型的過程如下。

步驟1利用離散短時傅里葉變換(discrete short time Fourier transform,DSTFT)將信號從時域轉換到頻域。當窗口長度固定時，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作為邊的權重。

(3)

2.2 相似性分數(shù)

(4)

對圖2中信號的5個頻段進行以上步驟分析，得到相似性分數(shù)如圖3所示?？梢钥闯?，θ波的相似性分數(shù)在異常變化位置沒有明顯變化，其他波段相似性分數(shù)變化較為明顯，主要原因是腦電信號的變化信息隱藏在不同的頻率成分中。因此，為了充分利用這些重要的微弱變化信息來實現(xiàn)對EEG信號的自適應檢測，需要一種自適應權重融合算法對所有的相似性分數(shù)進行加權融合。

圖3 相似性分數(shù)計算示例Fig.3 Example of calculation of the similarity scores

2.3 自適應權重融合算法

自適應權重通過依據(jù)各個分量的重要性自適應地計算權重，將所有相似性分數(shù)組合在一起。計算過程如下。

(5)

步驟2計算方差，其計算公式為

(6)

步驟3使用方差自適應地計算每個分量的權重wk，其計算公式為

(7)

步驟4所有頻段的相似性分數(shù)可表示為綜合性指標Qi，其計算公式為

(8)

2.4 假設檢驗

基于綜合性指標{Q1,Q2,…,Qm,…,QM}，采用常零假設檢驗進行變化決策。對于第m段，假設檢驗如下：H0:沒發(fā)生變化，Qm∈??；H1:發(fā)生變化，Qm?Α，其中，Α=[μm-1-3σm-1,μm-1+3σm-1]為由常用的3σ準則定義的置信區(qū)間，其中μm-1和σm-1分別為前m-1個片段綜合性指標的平均值和標準差。

圖4給出了檢測結果的示例?？梢钥闯觯褂玫募僭O檢驗可以成功地檢測出異常發(fā)生的時間，證明了所提出的綜合性指標的有效性。

圓圈表示測試方法產(chǎn)生的異常變化點圖4 檢測結果示例Fig.4 Example of detection result

3 實驗結果及討論

兩個不同數(shù)據(jù)集的實驗結果，與傳統(tǒng)線性指標如均值(mean)、方差(variance)、均方根(root mean square,RMS)和峰度(kurtosis)以及典型的非線性指標經(jīng)驗模態(tài)分解(empirical mode decomposition,EMD)、短時傅里葉變換(short-time fourier transform,STFT)、連續(xù)小波變換(continuous wavelet transform,CWT)和離散小波變換(discrete wavelet transform,DWT)進行了比較。采用常用指標查準率Precision、查全率Recall和F分數(shù)Fscore來評價以上方法的檢測性能，定義為

(9)

(10)

(11)

式中：TP為真正例，表示正確的檢測結果；FP為假正例，表示錯誤的檢測結果；FN為假負例，表示沒有檢測結果。

顯然，高查準率表示高檢測精度，高查全率表示高檢測靈敏度，F(xiàn)分數(shù)對檢測的方法進行綜合評價。

3.1 CHB-MIT頭皮EEG數(shù)據(jù)庫

由于CHB-MIT頭皮EEG數(shù)據(jù)庫腦電信號的數(shù)據(jù)量過于龐大，只截取癲癇發(fā)作開始時間點的前60 s及后60 s作為實驗數(shù)據(jù)。圖5展示了所提方法用于測試CHB-MIT頭皮EEG數(shù)據(jù)庫實驗結果的示例。

顯然，所提方法綜合性指標能夠及時檢測到腦電信號的微弱變化。由圖5可知，DWT有誤報點，EMD和STFT未能成功檢測到異常變化點，其余方法均表現(xiàn)良好的檢測結果,但是并不能通過此例證明其他方法比綜合性指標方法好。如表1所示，綜合性指標在查準率、查全率和F分數(shù)之間取得的效果最好，表現(xiàn)出了該方法的優(yōu)越性。盡管方差在查全率方面效果為100%，但是它的查準率和F分數(shù)沒有綜合性指標的性能要好，這也并不能說明方差在檢測性能上好于綜合性指標。

表1 CHB-MIT數(shù)據(jù)庫檢測結果進行了3個指標的比較Table 1 Comprehensive comparison in CHB-MIT scalp EEG database with three indicators

圓圈表示測試方法產(chǎn)生的異常變化點圖5 綜合性指標與其他方法的比較結果Fig.5 Comparison results of comprehensive indicators and other methods

3.2 山東大學第二醫(yī)院神經(jīng)內(nèi)科EEG數(shù)據(jù)庫

圖6展示了所提方法用于測試山東大學第二醫(yī)院神經(jīng)內(nèi)科EEG數(shù)據(jù)庫實驗結果的兩個示例。本實驗中測試數(shù)據(jù)的變化類型主要是幅值變化和頻率變化。圖6中測試數(shù)據(jù)1的異常變化主要是由幅值變化引起的?？梢钥闯觯司?、峰度和經(jīng)驗模態(tài)分解方法產(chǎn)生了誤報點，其余方法都能及時檢測到腦電數(shù)據(jù)的微弱異常變化。

檢測由頻率變化引起的腦電數(shù)據(jù)異常更具有挑戰(zhàn)性。如圖6(b)所示，很難從原始腦電信號中用肉眼直觀地識別腦電信號的微弱異常變化。經(jīng)過實驗，只有綜合性指標能夠及時檢測到腦電信號微弱的異常變化，其他方法均產(chǎn)生延遲或誤報點。表2中的綜合比較結果驗證了這一點，只有綜合性指標在查準率、查全率和F分數(shù)之間取得了最好的效果，并且查全率的結果為100%，這意味著該方法能夠檢測出所有潛在的腦電信號異常變化。腦電信號的非線性、非平穩(wěn)性和高變性使其極容易受外界影響，極大地增加了信號的異常檢測難度，所以綜合性指標在查準率方面未能達到百分之百也是有情可原?；谝陨辖Y果，綜合性指標具有很高的檢測精度和靈敏度，并且能夠準確反映腦電信號的動態(tài)變化。

圓圈表示測試方法產(chǎn)生的異常變化點圖6 山東大學數(shù)據(jù)集實驗結果比較Fig.6 Comparison of experimental results of Shandong University data set

表2 山東大學數(shù)據(jù)集檢測結果中三個指標的比較Table 2 Comprehensive comparison in Shandong University data set with three indicators

3.3 計算效率

在設計一個實用的變化檢測系統(tǒng)時，計算效率是一個需要考慮的重要問題。綜合性指標與其他比較方法的檢測步驟基本一樣，均包括數(shù)據(jù)建模、指標提取和決策3個階段。表3給出了綜合性指標與其他方法的計算復雜度的比較結果。結果表明，在所有的方法中，決策的計算復雜度是相同的，但在其他兩個階段的計算復雜度卻有很大的不同。顯然，綜合性指標的總計算復雜度是最大的，這意味著隨著腦電數(shù)據(jù)量的增加，該方法所消耗的時間也會大量增加。因此，在將來的工作中會進一步提高綜合性指標的計算效率，增加其實用價值。

表3 計算復雜度的比較Table 3 Comparison of the computational complexity

4 結論

提出了一種基于多頻段圖模型的腦電信號微弱異常變化檢測新方法。該方法充分利用這五個頻段中腦電的微弱變化信息和GM的優(yōu)勢，對腦電信號δ、θ、α、β和 γ波的每個頻段單獨進行圖模型建模，并通過自適應權重融合算法融合各圖模型的相似性分數(shù)形成綜合性指標進行假設檢驗。兩個實驗結果表明，所提方法查全率結果均為100%，表明所提方法能夠檢測所有潛在的腦電信號異常變化，因此該結果對腦電信號的微弱異常變化的檢測及實際應用具有重要意義。由于所提方法的計算復雜度較高，下一步會對所提方法進一步改進以提高其計算效率，增加其應用價值。