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

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

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

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

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

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

1 相關(guān)理論

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.1 aEEG 的獲取

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

通過(guò)棘波信號(hào)Z(k)可以計(jì)算棘波率（spike rate,SR），表示如下：

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

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

3 實(shí)驗(yàn)結(jié)果與分析

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

為了評(píng)估本文所提出方法的有效性和優(yōu)越性，本文通過(guò)計(jì)算混淆矩陣（confusion matrix,CM）來(lái)體現(xiàn)該算法得到的檢測(cè)結(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）是評(píng)估癲癇發(fā)作檢測(cè)算法性能的3 個(gè)重要指標(biāo)，可由CM 計(jì)算得到，其中AC為檢測(cè)正確的EEG 片段的數(shù)量與所有EEG 片段數(shù)量的比值，可通過(guò)式（17）計(jì)算得到。SE 表示癲癇發(fā)作檢測(cè)算法正確識(shí)別的發(fā)作次數(shù)與實(shí)際發(fā)作事件數(shù)之比，它衡量了本文所提出的方法對(duì)癲癇發(fā)作的識(shí)別能力，由式（18）計(jì)算得到。SP 代表實(shí)際為癲癇發(fā)作期并被檢測(cè)為發(fā)作間期的 EEG 片段占所有發(fā)作間期EEG 片段比例，反映了分類器檢驗(yàn)負(fù)樣本的能力，由式（19）表示。

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

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

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

3.2 結(jié)果分析

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)語(yǔ)

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