基于aEEG尖峰和cEEG棘波提取的癲癇發(fā)作檢測算法

吳端坡，王紫萌，董 芳，馮 維，劉兆霆

（1.杭州電子科技大學(xué) 通信工程學(xué)院，浙江 杭州 310018；2.浙大城市學(xué)院 信息與電氣工程學(xué)院，浙江 杭州 310015；3.浙江環(huán)瑪信息科技有限公司，浙江 杭州 310052）

癲癇是大腦神經(jīng)元突發(fā)性異常放電導(dǎo)致短暫的大腦功能障礙的一種慢性疾病[1]。目前在全球有超過7 000萬人患有癲癇病，中國的癲癇患者有900萬之多[2]。目前癲癇疾病診斷的主要依據(jù)是癲癇癥狀學(xué)、連續(xù)腦電圖（continuous electroencephalography,cEEG）和振幅整合腦電圖（amplitude integrated electroencephalography,aEEG）的檢測結(jié)果[3-4]。但是目前的神經(jīng)科醫(yī)生在利用腦電數(shù)據(jù)進(jìn)行癲癇確診時仍采用腦電技師人工對采集的腦電數(shù)據(jù)進(jìn)行逐秒數(shù)據(jù)篩查，工作量大，而且需要依靠讀圖人的個人經(jīng)驗，診斷結(jié)果可能會產(chǎn)生誤判漏判[5]。因此癲癇發(fā)作自動檢測算法得到了廣泛的關(guān)注。

棘波是典型的癲癇特征波形，相對于背景波形，棘波波形尖銳，具有高幅和瞬變的特性，目前臨床上的癲癇檢查主要是通過人眼檢測識別腦電信號的棘波[6]。棘波檢測的方法有很多，Saltzberg 等[7]早在1967 年研究了一種基于棘波周期、幅值和二階導(dǎo)數(shù)等參數(shù)的棘波檢測方法。小波變換由于具有良好的時頻局部化性質(zhì)，在棘波檢測領(lǐng)域中已開展應(yīng)用。Kalayci 等[8]對EEG 信號進(jìn)行小波分解，得到的小波系數(shù)作為神經(jīng)網(wǎng)絡(luò)的輸入?yún)?shù)，但對背景信號的抑制效果較差，棘波檢測的效果不夠理想。

MATLAB 是信號處理與分析領(lǐng)域常用的數(shù)據(jù)分析平臺[9-10]，“MATLAB 與仿真”是我院一直開設(shè)的經(jīng)典課程。本文在上述研究的基礎(chǔ)上，根據(jù)癲癇發(fā)作時的臨床特征，并結(jié)合aEEG 和cEEG 的信號特征，通過MATLAB 仿真分析了癲癇發(fā)作時期棘波的數(shù)目與時間的關(guān)系，提出了一種基于aEEG 尖峰和cEEG 棘波提取的聯(lián)合癲癇檢測方法。

1 相關(guān)理論

1.1 數(shù)據(jù)描述與數(shù)據(jù)預(yù)處理

本文使用公開的CHB-MIT 癲癇腦電數(shù)據(jù)庫進(jìn)行聯(lián)合癲癇發(fā)作檢測研究[11]。該數(shù)據(jù)庫的數(shù)據(jù)來源于波士頓兒童醫(yī)院，包含了24 組數(shù)據(jù)（chb01-chb24），分別來自23 名患者（chb01 和chb21 來自同一名女性患者），其中包括6 名男性，年齡在3—22 歲；17 名女性，年齡在1.5—19 歲。如圖1 所示，所有的EEG 數(shù)據(jù)都是通過國際10-20 電極分布系統(tǒng)采集的[12]。該數(shù)據(jù)庫中EEG 信號的采樣頻率為256 Hz，分辨率為16 bits，所有的數(shù)據(jù)都已經(jīng)由專業(yè)的醫(yī)生進(jìn)行了發(fā)作開始和發(fā)作結(jié)束的標(biāo)記。數(shù)據(jù)庫中的數(shù)據(jù)大部分是23 導(dǎo)聯(lián)的，也有一部分是22 或24 導(dǎo)的，本文使用其中的19導(dǎo)數(shù)據(jù)進(jìn)行癲癇發(fā)作檢測的研究。圖2 和3 分別顯示了癲癇發(fā)作間期和發(fā)作期的腦電信號。

圖1 國際10-20 電極分布系統(tǒng)

圖2 癲癇發(fā)作間期腦電信號

圖3 癲癇發(fā)作期腦電信號

在EEG 信號采集過程中，由于設(shè)備本身、肌肉活動和眨眼等因素容易產(chǎn)生噪聲和偽跡，在本研究中，使用1～48 Hz 的帶通濾波器濾除干擾噪聲和偽跡。癲癇發(fā)作持續(xù)時間短，且容易反復(fù)，為了更好地提取癲癇發(fā)作的棘波特征，將連續(xù)的EEG 信號分為若干個包含了19 個通道的EEG 片段，長度5 s，有4 s 重疊，并對每個EEG 片段進(jìn)行分析處理。

1.2 形態(tài)學(xué)濾波器

形態(tài)學(xué)濾波器是一種非線性濾波器。它利用信號的幾何特征將信號與預(yù)先設(shè)定的結(jié)構(gòu)元素進(jìn)行匹配，提取出與結(jié)構(gòu)元素相似的信號特征[13-14]。形態(tài)學(xué)變換的基本操作包括腐蝕、膨脹、形態(tài)開運算和形態(tài)閉運算，其定義分別如式（1）—（4）所示。

其中，f(n)代表輸入的時間序列，g(n)代表結(jié)構(gòu)元素。f(n)和g(n)的序列長度分別為N和M，并且N＞＞M。符號?、⊕、?、·分別代表腐蝕、膨脹、形態(tài)學(xué)開操作和形態(tài)學(xué)閉操作。

形態(tài)學(xué)開操作和閉操作對信號處理的影響是不同的：形態(tài)學(xué)開操作可以平滑信號中的正向脈沖，形態(tài)學(xué)閉操作可以平滑信號中的負(fù)向脈沖。為了同時去除信號中的正、負(fù)脈沖，Maragos[15]提出了一種將開、閉運算級聯(lián)組合的方式，構(gòu)造出形態(tài)開-閉（open-closing,OC）和形態(tài)閉-開（close-opening,CO）濾波器。

形態(tài)OC 和形態(tài)CO 運算公式分別為：

2 癲癇發(fā)作檢測算法

算法主要包括數(shù)據(jù)預(yù)處理、振幅整合、基于aEEG尖峰提取的發(fā)作片段檢測和基于cEEG 棘波提取的癲癇發(fā)作檢測4 部分。算法流程如圖4 所示。

圖4 癲癇發(fā)作檢測算法整體流程圖

2.1 aEEG 的獲取

為了提高癲癇發(fā)作檢測算法的效率，本文將cEEG信號轉(zhuǎn)化為aEEG 信號，然后采用尖峰提取的方法對癲癇發(fā)作進(jìn)行檢測，確定癲癇發(fā)作的主要通道。原始cEEG 轉(zhuǎn)換成aEEG 的過程如下：

（1）非對稱的二次窄帶濾波。首先對原始cEEG數(shù)據(jù)進(jìn)行2～15 Hz 的二次窄帶非對稱濾波，剔除與癲癇發(fā)作聯(lián)系極小的極低幅極低頻分量的快波，進(jìn)而分析對檢測癲癇發(fā)作有意義的頻率范圍。

（2）幅值壓縮。如果將原始的cEEG 信號直接壓縮在時間軸上，一些較低的振幅分量會丟失，腦電信號中的癲癇病理特征缺失會影響發(fā)作檢測的準(zhǔn)確性。因此，aEEG 以半對數(shù)方式進(jìn)行振幅整合，縮小原始cEEG 信號的幅值變化范圍，使反映腦功能變化的aEEG信號能夠顯示在有限的圖紙上。具體步驟如下，首先取原始cEEG 信號幅值的絕對值，然后選取5 μV 作為閾值進(jìn)行判斷。如果幅值＜5 μV，那么該點的幅值保持不變，如果幅值在5 μV 以上，則執(zhí)行半對數(shù)運算。

（3）時間壓縮。對原始cEEG 信號的時間軸進(jìn)行高度壓縮可以使長時腦電信號反映在aEEG 較短的一段波形中，便于醫(yī)生從宏觀角度觀察病人的腦功能狀態(tài)。本文考慮到所用數(shù)據(jù)的采樣率為256 Hz，因此將壓縮比設(shè)置為256∶1。

（4）下包絡(luò)提取。aEEG 可以明顯地反映原始腦電信號的波動變化趨勢，為此，可以將aEEG 信號分割成5 s 長的片段，取每個片段的最小值，得到aEEG信號的下包絡(luò)，用于后續(xù)的棘波檢測。

2.2 基于aEEG 尖峰提取的發(fā)作片段檢測

對aEEG 下包絡(luò)進(jìn)行尖峰提取的第一步就是選擇一個合適的結(jié)構(gòu)元素，這對尖峰檢測的結(jié)果有很大的影響。為了去除cEEG 信號中的瞬時分量，并將背景腦電信號從原始信號中分離出來，本文選擇了能夠最大程度反映尖峰信號幾何特征的三角形結(jié)構(gòu)元素，如圖5 所示，公式如下：

其中，A代表結(jié)構(gòu)元素的中心高度，2L代表結(jié)構(gòu)元素的寬度。

圖5 三角形結(jié)構(gòu)元素結(jié)構(gòu)

因為aEEG 下包絡(luò)表示了cEEG 信號的最小波幅變化，因此為了提高尖峰提取的精度，進(jìn)一步抑制背景信號，本文使用了一種改良的形態(tài)學(xué)濾波方法，設(shè)置兩個結(jié)構(gòu)元素，根據(jù)數(shù)據(jù)集中所有數(shù)據(jù)的癲癇發(fā)作平均時間來確定結(jié)構(gòu)元素的寬度，經(jīng)統(tǒng)計，癲癇平均發(fā)作時間為46 s，即2L的值為46。此外，數(shù)據(jù)集中aEEG 信號下包絡(luò)的最小和最大波幅為結(jié)構(gòu)元素高度A的范圍，為[8.5,10]。因此對aEEG 下包絡(luò)進(jìn)行棘波提取所用的兩個結(jié)構(gòu)元素表示如下：

aEEG 下包絡(luò)信號P(k)經(jīng)形態(tài)學(xué)OC 操作和形態(tài)學(xué)CO 操作得到的結(jié)果如下所示：

由于開操作具有擴(kuò)展性而閉操作具有反擴(kuò)展性，因此單獨進(jìn)行OC 操作或CO 操作會使結(jié)果出現(xiàn)統(tǒng)計偏倚現(xiàn)象。針對以上現(xiàn)象，本文使用平均加權(quán)的方法消除統(tǒng)計偏倚。aEEG 下包絡(luò)信號P(k)經(jīng)過OC 操作和CO 操作的平均加權(quán)后，得到的輸出信號即為背景信號Q(k)：

將aEEG 下包絡(luò)信號與背景信號Q(k)相減得到尖峰信號S(k)：

最后使用閾值法判斷該aEEG 下包絡(luò)片段是否為癲癇發(fā)作，設(shè)發(fā)作閾值為：

最后采用投票法確定癲癇發(fā)作的位置，當(dāng)一個腦電片段中尖峰幅值＞發(fā)作閾值thresholda的通道數(shù)≥10，則說明該腦電片段為發(fā)作片段C。

2.3 基于cEEG 棘波檢測的癲癇發(fā)作檢測

對cEEG 的棘波提取步驟與aEEG 的尖峰提取步驟類似，最大的不同之處就是結(jié)構(gòu)元素的選擇，通常情況下，棘波的幅值在150 μV 以上，持續(xù)時間t為20～70 ms，cEEG 信號的采樣率fs為256 Hz，記2Lc為cEEG 棘波檢測中結(jié)構(gòu)元素的寬度，Lc的計算公式為，因此設(shè)Ac=200 μV，Lc的取值范圍為2～9，當(dāng)Lc=5 時，可以有效地利用CHB-MIT EEG 數(shù)據(jù)提取cEEG 的棘波。

設(shè)cEEG 腦電信號為C(k)。與aEEG 尖峰提取方法相同，背景信號棘波信號Z(k)=C(k)-Y(k)。

通過棘波信號Z(k)可以計算棘波率（spike rate,SR），表示如下：

其中，numj(i)表示第i個片段信號中第j個通道的棘波數(shù)量，I表示aEEG 下包絡(luò)片段的總數(shù)目。在本文中，cEEG 片段長度T為5 s。

根據(jù)式（14）和（15）計算cEEG 的發(fā)作閾值thresholdc和SR，采用投票法判斷是否有癲癇發(fā)作，當(dāng)19 通道腦電片段中SR＞thresholdc的通道數(shù)≥10，則說明該片段確實為癲癇發(fā)作片段。

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

3.1 性能評估指標(biāo)

為了評估本文所提出方法的有效性和優(yōu)越性，本文通過計算混淆矩陣（confusion matrix,CM）來體現(xiàn)該算法得到的檢測結(jié)果與專家標(biāo)記結(jié)果之間的差異。Smn表示被標(biāo)記為標(biāo)簽m并被歸類為標(biāo)簽n的EEG 片段數(shù)量。在CM 中，癲癇發(fā)作期和癲癇發(fā)作間期分別作為標(biāo)簽1 和標(biāo)簽2。CM 表示如下：

準(zhǔn)確性（accuracy,AC）、敏感性（sensitivity,SE）和特異性（specificity,SP）是評估癲癇發(fā)作檢測算法性能的3 個重要指標(biāo)，可由CM 計算得到，其中AC為檢測正確的EEG 片段的數(shù)量與所有EEG 片段數(shù)量的比值，可通過式（17）計算得到。SE 表示癲癇發(fā)作檢測算法正確識別的發(fā)作次數(shù)與實際發(fā)作事件數(shù)之比，它衡量了本文所提出的方法對癲癇發(fā)作的識別能力，由式（18）計算得到。SP 代表實際為癲癇發(fā)作期并被檢測為發(fā)作間期的 EEG 片段占所有發(fā)作間期EEG 片段比例，反映了分類器檢驗負(fù)樣本的能力，由式（19）表示。

其中，NUMtest表示通過aEEG 尖峰和cEEG 棘波提取的癲癇發(fā)作檢測（M3）方法檢測到的癲癇發(fā)作事件數(shù)目；NUMreal表示實際的癲癇發(fā)作事件數(shù)據(jù)。

此外，該方法還計算了癲癇發(fā)作檢測中應(yīng)用廣泛的基于癲癇事件的假陽性率（false positive ratio based on the event,FPRE），進(jìn)一步評價了本文所提出方法的性能和通用性。FPRE 表示1 h 內(nèi)的誤警次數(shù)，計算方法由式（20）表示。

其中，F(xiàn)P 表示檢測到癲癇發(fā)作間期的事件數(shù)，TD 表示EEG 信號記錄的總持續(xù)時間。事件代表從發(fā)作開始到發(fā)作結(jié)束的時間，可能包含多個EEG 片段。

3.2 結(jié)果分析

在CHB-MIT 數(shù)據(jù)集中，本文首先使用了aEEG的下包絡(luò)進(jìn)行癲癇發(fā)作片段的檢測，然后對該發(fā)作片段的cEEG 信號進(jìn)行棘波提取，通過閾值檢測和投票法實現(xiàn)癲癇發(fā)作區(qū)域的檢測。

圖6 給出了chb04 腦電數(shù)據(jù)不同階段的癲癇發(fā)作檢測結(jié)果。其中1 代表癲癇發(fā)作期，0 代表發(fā)作間期。圖6(a)是由臨床醫(yī)生標(biāo)記的2 起真實發(fā)作事件；圖6(b)是經(jīng)過aEEG 尖峰檢測得到的5 起發(fā)作事件；圖6(c)是由cEEG 棘波檢測得到的4 起發(fā)作事件；圖6(d)是經(jīng)過aEEG 尖峰檢測和cEEG 棘波檢測后得到的2 起發(fā)作事件。可以看出，本文提出的癲癇發(fā)作檢測方法得到的結(jié)果與醫(yī)生真實標(biāo)注的發(fā)作事件相吻合。

圖6 不同方法得到的癲癇發(fā)作檢測結(jié)果

為了評價該方法的通用性，本文使用AC、SP、SE 和FPRE 對所提出的檢測方法性能進(jìn)行了評價。根據(jù)每個患者的檢測結(jié)果列出表1。其中，chb06 病例在 實際中有10 次癲癇發(fā)作，但M3 方法未檢出癲癇發(fā)作。同樣地，chb14 和chb15 分別實際發(fā)作8 次和20 次，檢測到發(fā)作次數(shù)為0，這大大降低了SE 的整體性能。從表1 中可以看出，該方法可以針對不同患者進(jìn)行癲癇發(fā)作檢測，具有較強(qiáng)的廣泛性。

為了更清晰地表現(xiàn)本文所提出的癲癇發(fā)作檢測算法的優(yōu)越性，圖7 給出了不同檢測算法的性能指標(biāo)。可以看出，使用aEEG+cEEG 進(jìn)行癲癇發(fā)作檢測的AC和SP 指標(biāo)要高于僅使用aEEG 和僅使用cEEG 的癲癇發(fā)作檢測方法。

表2 是本文研究的方法與其他方法性能指標(biāo)的比較，所有方法均使用CHB-MIT 數(shù)據(jù)集進(jìn)行驗證。其中文獻(xiàn)[16]采用基于腦電圖信號動力學(xué)的癲癇發(fā)作檢測方法，結(jié)合主成分分析得到表2 中的指標(biāo)。文獻(xiàn)[17]提出了一種基于平穩(wěn)小波變換的非特異性癲癇發(fā)作檢測方法。文獻(xiàn)[18]為了解決長期腦電圖記錄中癲癇病發(fā)作的離線監(jiān)測問題，提出了一種基于稀疏表示和局部Gabor 二值模式的癲癇發(fā)作檢測算法。文獻(xiàn)[19]采用快速傅立葉變換和功率譜密度的方法實現(xiàn)癲癇發(fā)作智能檢測。

表1 24 組數(shù)據(jù)的檢測結(jié)果

圖7 不同方法的性能指標(biāo)對比

表2 本方法與其他方法的比較

由表2 可以看出，本方法的AC、SP 和FPRE 均高于其他方法，但SE 指標(biāo)略低。出現(xiàn)這種現(xiàn)象的原因有兩種，其一是因為癲癇種類多，低幅值的發(fā)作無法通過aEEG 識別到，例如患者chb14 在發(fā)作的時候幅值沒有明顯升高，所以無法通過aEEG 形成拱門；其二是本方法的SE 和其他方法中的計算方法不同，本文是從醫(yī)學(xué)角度出發(fā)，以事件為基礎(chǔ)的方式進(jìn)行計算，更符合醫(yī)學(xué)領(lǐng)域的性能評估方式。

文獻(xiàn)[17]并沒有對數(shù)據(jù)集中所有的數(shù)據(jù)進(jìn)行研究，因此文獻(xiàn)[17]中的FPRE 可能比實際值高。本方法對CHB-MIT 數(shù)據(jù)集中目前所有的數(shù)據(jù)都進(jìn)行了檢測分析，所以FPRE 的性能得到了提高。

4 結(jié)語

本文通過MATLAB 仿真平臺對CHB-MIT 數(shù)據(jù)庫中的24 組數(shù)據(jù)進(jìn)行發(fā)作檢測的研究，并對每組數(shù)據(jù)的SR 進(jìn)行了統(tǒng)計分析。可以得出，發(fā)作間期的棘波數(shù)最少，SR 最低，隨著時間的推移，距離發(fā)作期越近，SR逐漸增大，在發(fā)作結(jié)束后逐漸減少。本文提出并評估了基于aEEG 尖峰和cEEG 棘波提取的癲癇發(fā)作檢測算法，采用本方法對24 組數(shù)據(jù)進(jìn)行癲癇檢測，準(zhǔn)確率最高可達(dá)99.9%，特異性高達(dá)99%，每小時的誤檢次數(shù)最低為0 次。在以后的研究中還可以深入分析發(fā)作間期與癲癇棘波之間的關(guān)系，為癲癇發(fā)作預(yù)警技術(shù)提供理論與實驗支持。