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

賀光碩，盧國梁*，尚偉

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

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

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

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

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

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

1 實驗數據來源

1.1 腦電數據庫的選取

所用數據庫選自公開的CHB-MIT頭皮EEG數據庫[16](網址：https://physionet.org/content /chbmit/1.0.0/)和來自山東大學第二醫院神經內科的EEG數據庫[11,17]。

1.1.1 CHB-MIT頭皮EEG數據庫

該數據庫由患有難治性癲癇病的兒科受試者的腦電圖數據記錄組成，收集了22名受試者共24個案例，采樣率為256 Hz。每個案例(chb01、chb02等)包含來自一名受試者的9～42個連續.edf文件并且大多數文件的時間跨度為1 h。如圖1所示，每個案例詳細地標注了癲癇發作的時間和次數，方便了實驗時對異常變化點的標注。在每種案例中，隨機選擇一個含有癲癇發作的.edf文件作為測試數據，共24個數據，選擇第16個通道數據作為實驗所需的測試數據。

圖1 原始腦電信號癲癇發作時間(以chb01_03為例)Fig.1 Raw EEG signal seizure time(take chb01_03 as an example)

1.1.2 山東大學第二醫院神經內科EEG數據庫

該數據庫的所有EEG數據均來自癲癇患者，且所有受試者均已獲得書面知情同意。該實驗完全遵守《赫爾辛基宣言》和當地相關法律法規。實驗所用的癲癇檢測系統用于監測夜間睡眠中的癲癇[17]，同時收集受試者的EEG數據和相應的監控視頻。

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

1.2 數據的預處理

癲癇等神經源性疾病的頻率成分主要出現在δ、θ、α、β和 γ頻段[15,18]。在神經性疾病發作時，這些頻率的幅值會相應改變。因此，根據以上5個頻段的頻率定義一個0.1～70 Hz的帶通濾波器來捕獲該頻率的腦電信號。此外，使用特定的陷波濾波器來消除50 Hz的電力線的干擾。

根據帶通濾波器的截止頻率和腦電信號實際異常變化情況，腦電信號各個波段頻率范圍為：δ波(0.1～4 Hz)、θ波(4～8 Hz)、α波(8～12 Hz)、β波(12～30 Hz)、γ波(31～70 Hz)。對于一段給定的濾波后的腦電信號，分別進行5個頻段的提取，如圖2所示。可以看出，腦電信號的5個頻段在異常變化檢測中具有不同表現，如α、β和 γ頻段在異常變化位置反應明顯，而δ和θ頻段則無明顯變化，因此所提算法對于腦電信號的微弱異常變化檢測具有重要意義。

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

2 實驗方法

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

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

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

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

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

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

2.1 構造圖模型

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

構造圖模型的過程如下。

步驟1利用離散短時傅里葉變換(discrete short time Fourier transform,DSTFT)將信號從時域轉換到頻域。當窗口長度固定時，DSTFT的計算公式為

(1)

式(1)中:φ*(·)為窗口函數；p為DSTFT得到的頻率數據；i=1,2,…,N/2；j為虛數；N為窗口函數移動的時間步長，N=512；m為第m個時間步長。

步驟2計算功率譜，通過式(2)可以得到一系列周期圖P。

(2)

步驟3周期圖P(i,m)用于定義時間t處的頻率分辨率，即Fm={f(i)},i=1,2,…,N/2。其中，f(i)為頻率，Fm為頻率f(i)的集合。將每個頻率分辨率視為一個節點，計算每兩個節點之間的歐氏距離di,j作為邊的權重。

(3)

2.2 相似性分數

(4)

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

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

2.3 自適應權重融合算法

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

(5)

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

(6)

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

(7)

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

(8)

2.4 假設檢驗

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

圖4給出了檢測結果的示例。可以看出，使用的假設檢驗可以成功地檢測出異常發生的時間，證明了所提出的綜合性指標的有效性。

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

3 實驗結果及討論

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

(9)

(10)

(11)

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

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

3.1 CHB-MIT頭皮EEG數據庫

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

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

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

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

3.2 山東大學第二醫院神經內科EEG數據庫

圖6展示了所提方法用于測試山東大學第二醫院神經內科EEG數據庫實驗結果的兩個示例。本實驗中測試數據的變化類型主要是幅值變化和頻率變化。圖6中測試數據1的異常變化主要是由幅值變化引起的。可以看出，除了均值、峰度和經驗模態分解方法產生了誤報點，其余方法都能及時檢測到腦電數據的微弱異常變化。

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

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

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

3.3 計算效率

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

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

4 結論

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