王棟,李擴,劉曉芳,閆相國,王剛

(1.西安交通大學(xué)生物醫(yī)學(xué)信息工程教育部重點實驗室,710049,西安;2.西安交通大學(xué)第一附屬醫(yī)院神經(jīng)外科,710061,西安)

癲癇是一種常見的慢性腦部疾病,其發(fā)作表現(xiàn)為大腦神經(jīng)元過度同步化放電,并導(dǎo)致短暫性中樞神經(jīng)系統(tǒng)功能失常。腦電圖(EEG)是一種非侵入性的測量腦電活動的工具,其中蘊含大量的腦功能信息,對大腦疾病的診斷具有很高的價值,也可以用作癲癇發(fā)作的檢測[1]。在多數(shù)情況下,癲癇發(fā)作無法在短期內(nèi)預(yù)測,因此需要在較長的時間內(nèi)連續(xù)記錄腦電。長程腦電圖監(jiān)測可以有效地提供關(guān)于大腦活動、癲癇發(fā)作頻率等的信息,有助于臨床醫(yī)生對癲癇的診斷和治療。長程腦電經(jīng)過專家的目測檢查,可以確定癲癇發(fā)作的情況,但是這是一項相當(dāng)乏味而耗時的任務(wù),繁重的工作量容易引起醫(yī)護(hù)人員疲勞而導(dǎo)致人工檢測結(jié)果的準(zhǔn)確性下降[2],并且這很依賴于臨床醫(yī)生的自身經(jīng)驗和主觀判斷,不同的醫(yī)護(hù)人員所作的診斷有可能不一致,結(jié)果的可重復(fù)性較低。除此之外,如果在癲癇發(fā)作后較短的時間內(nèi)檢測到,也利于醫(yī)護(hù)人員采取相應(yīng)的措施來減少對患者的傷害。因此,計算機自動檢測技術(shù)一直是癲癇腦電圖監(jiān)測工作中的迫切需求。

近幾十年來,出現(xiàn)了各種各樣的癲癇自動檢測的方法,其中特征提取方法主要包括時域、頻域、時頻域和非線性的方法[3]。多方研究表明,腦電波信號表現(xiàn)為一種非線性、非平穩(wěn)的隨機過程,因此關(guān)聯(lián)維[4]、Lyapunov指數(shù)[5]、熵[6]等度量特征已經(jīng)被引入到信號分析和特征提取中,小波變換也被廣泛地應(yīng)用到了該領(lǐng)域[7-8]。除了單一通道的腦電,越來越多的人開始關(guān)注多通道腦電之間的相互關(guān)系,關(guān)注腦區(qū)各部分的信息交流[9-10]。有向傳遞函數(shù)(DTF)是通過腦電來計算腦區(qū)信息交流的方法之一[11]?，F(xiàn)有的大多數(shù)的癲癇發(fā)作檢測針對的都是低頻段,忽略頭皮腦電的高頻成分。

本文首先利用小波分解將腦電信號中的高頻γ波段提取出來,利用該波段腦電信號的流出信息強度作為特征來進(jìn)行癲癇發(fā)作檢測。結(jié)果表明,癲癇發(fā)作前后,高頻段的腦電流出信息強度差異很大,利用長程腦電的高頻成分對局灶性癲癇患者癲癇發(fā)作檢測有著很高的靈敏性和特異性。

1 數(shù)據(jù)和方法

1.1 癲癇患者和腦電數(shù)據(jù)

本文長程腦電數(shù)據(jù)是從西安交通大學(xué)第一附屬醫(yī)院10例臨床癲癇患者身上采集的。數(shù)據(jù)采集儀器是NIHON KOHDEN公司的EEG-1100腦電圖機,使用19導(dǎo)電極,電極位置按照國際10—20標(biāo)準(zhǔn),以CZ電極為參考電極,數(shù)據(jù)的通帶截止頻率為0.5～60 Hz,采樣頻率為200 Hz。癲癇發(fā)作的開始點和終止點由西安交通大學(xué)第一附屬醫(yī)院的兩位具有臨床經(jīng)驗的醫(yī)生通過視覺檢測和評估患者的實際臨床記錄和表現(xiàn)進(jìn)行了標(biāo)注。我們又人工剔除了原始腦電數(shù)據(jù)中的肌電等偽跡腦電。本文獲得的已采集腦電數(shù)據(jù)的總時長為121.8 h,包含了44次發(fā)作,其中癲癇發(fā)作期的平均時間為98.7 s。表1給出了每位患者的相關(guān)腦電數(shù)據(jù)信息。

1.2 癲癇檢測算法

癲癇腦電數(shù)據(jù)分為發(fā)作期腦電和非發(fā)作期腦電,發(fā)作檢測的目的是在長程腦電中識別發(fā)作期腦電。已有研究表明,癲癇發(fā)作與腦電的高頻成分有關(guān)系,癲癇發(fā)作時會有高頻腦電成分出現(xiàn)。在癲癇發(fā)作期和非發(fā)作期,頭皮腦電的高頻成分在不同電極間的信息流有可能是不同的,所以本文利用DTF算法來計算高頻腦電在不同腦區(qū)的信息流通情況。信息流是有方向的,即流出信息和流入信息,本文將各個電極的流出信息強度作為特征,用于癲癇檢測。本文利用高頻腦電成分進(jìn)行癲癇發(fā)作檢測的流程如圖1所示。首先將每位病人19通道的腦電數(shù)據(jù)進(jìn)行分割,分為發(fā)作腦電數(shù)據(jù)和非發(fā)作期腦電數(shù)據(jù)。對每個通道的腦電在滑動的時間窗內(nèi)用小波分解的方式提取高頻成分(γ波段),然后對19通道的高頻腦電建立多變量自回歸模型(MVAR)。利用有向傳遞函數(shù)來提取該頻段上各通道的信息交流特征,求出各通道的流出信息強度用于特征降維,最后經(jīng)SVM分類,通過五重交叉驗證得到檢測結(jié)果。

1.2.1 高頻腦電信號的提取 為了使用小波變換的方式提取出高頻波段的信號,根據(jù)小波分解的原理,當(dāng)信號的頻率為256 Hz時,對小波進(jìn)行兩層分解,即可得到大于32 Hz的腦電信號。小波分解中小波的選取非常重要,Db小波群對非平穩(wěn)信號(比如腦電信號)具有更好的分解效果[7-8],而db4小波更適合分析癲癇腦電信號[12]。因為原始腦電數(shù)據(jù)的采樣頻率為200 Hz,所以先進(jìn)行了一個256 Hz的重采樣,然后選擇db4小波對各個通道的腦電信號進(jìn)行兩層的小波分解,就能得到高頻腦電信號(>32 Hz)。

1.2.2 基于有向傳遞函數(shù)的特征提取 由于采集的腦電信號最小頻率為0.5 Hz,為了保證能達(dá)到最低分辨率,對腦電信號的分割需要不少于2 s[10]。本文對各個頻段的腦電信號進(jìn)行2 s無重疊的分割,利用有向傳遞函數(shù)進(jìn)行特征提取的第一步需要建立多變量自回歸模型(MVAR),然后在求得MAVR模型系數(shù)的基礎(chǔ)之上,利用有向傳遞函數(shù)(DTF)進(jìn)行特征值提取,該特征能夠反映不同導(dǎo)聯(lián)腦電信號之間的信息流動情況。

定義t時刻的N導(dǎo)聯(lián)(癲癇腦電的導(dǎo)聯(lián)數(shù)N=19)腦電信號為

X(t)=[X1(t),X2(t),…,XN(t)]T

(1)

式中:Xn(t)(n=1,2,…,N)是第n個通道的腦電信號。通過建立多通道自回歸模型,腦電序列可以表示如下

(2)

式中:p為MVAR模型的階數(shù);Ar為N×N的系數(shù)矩陣,r=1,2,…,p;E(t)為估計誤差,理想情況下是均值為0的非相關(guān)白噪聲。模型的階數(shù)p可以通過Akaike信息準(zhǔn)則(AIC)來確定[13]

F(p)=ln|∑(p)|+2pN2/M

(3)

式中:∑(p)為殘差的協(xié)方差矩陣;N為通道數(shù);M是全部用于計算的數(shù)據(jù)點數(shù)。選擇最小的F對應(yīng)的階數(shù),系數(shù)矩陣Ar的估計可以用Levison Wiggins Robinson(LWR)算法[14]求得。

對式(2)兩邊作傅里葉變換,得到

(4)

(5)

傳遞矩陣定義為

(6)

進(jìn)而可以獲得在頻率f上從導(dǎo)聯(lián)j到導(dǎo)聯(lián)i的信息流強度為

Dij(f)=

(7)

式中:Hij(f)為矩陣H(f)的第i行第j列元素;hi(f)為矩陣的第i列;Dij(f)表示在頻率f時,從導(dǎo)聯(lián)j到導(dǎo)聯(lián)i的信息流的強度和方向。

1.2.4 支持向量機分類 支持向量機(SVM)是一種新的機器學(xué)習(xí)方法,基本思想是:對于線性不可分樣本,經(jīng)非線性變換,將其映射到另一個高維空間中,在變換后的空間中尋找一個最優(yōu)的分界面(超平面),使之線性可分[15]。支持向量機對于二分類問題有著明顯的優(yōu)勢,本文選擇徑向基核函數(shù)用于發(fā)作期腦電和非發(fā)作期腦電的分類,采用五重交叉驗證的方法得到最終結(jié)果:將每一個病人的腦電數(shù)據(jù)平均分為5組,每次用4組數(shù)據(jù)做訓(xùn)練,一組數(shù)據(jù)用來測試,重復(fù)5次,保證每組都只測試1次,將5次的結(jié)果進(jìn)行平均作為最終測試結(jié)果。

2 結(jié) 果

算法的性能是根據(jù)每個患者臨床信息的分析統(tǒng)計結(jié)果的平均值來評估的,以避免由于單個患者有多次癲癇發(fā)作而引起的可能的偏差。通過比較本文算法對EEG片段做出的癲癇發(fā)作、非發(fā)作標(biāo)記以及醫(yī)生對此做出的標(biāo)記的差別來計算算法的性能。其中,真陽性TP是算法判斷正確的癲癇發(fā)作期片段,假陽性FP是算法判斷錯誤的癲癇發(fā)作期片段,真陰性TN是算法判斷正確的癲癇非發(fā)作期片段,假陰性FN是算法判斷錯誤的癲癇非發(fā)作期片段。算法的性能可以由以下指標(biāo)進(jìn)行評估[14,16],這些指標(biāo)的定義如下。

(1)正確率

%

(8)

(2)選擇性

%

(9)

(3)靈敏性

%

(10)

(4)特異性

%

(11)

(5)平均檢出率

ADR=0.5(SPE+SEN)×100%

(12)

本文獲得的已采集腦電數(shù)據(jù)的總時長為121.8 h,包含了44次發(fā)作,其中癲癇發(fā)作期的平均時間長度為98.7 s。以2 s時間窗作為分割,訓(xùn)練樣本的平均非發(fā)作期片段數(shù)為11 959,平均發(fā)作片段數(shù)為157,測試樣本的平均非發(fā)作期片段數(shù)為3 313,平均發(fā)作期片段數(shù)為39。表2為利用γ波段進(jìn)行癲癇發(fā)作檢測的結(jié)果,10位病人的平均正確率為98.4%,平均選擇性為60.7%,平均敏感性為93.4%,平均特異性為98.4%,平均檢出率為95.9%。

表2 基于高頻腦電的局灶性癲癇患者 癲癇發(fā)作檢測結(jié)果

3 討 論

3.1 基于DTF不同波段癲癇檢測結(jié)果對比

圖2 基于有向傳遞函數(shù)(DTF)算法的不同波段癲癇檢測結(jié)果

為了討論不同子頻帶的特征對癲癇檢測結(jié)果的影響,本文用不同子頻帶的特征進(jìn)行SVM分類。首先利用5層小波分解得到4種典型的腦電信號波段,即δ波(0.5～4 Hz)、θ波(4～8 Hz)、α波(8～16 Hz)、β波(16～32 Hz),然后采用DTF算法求出這4個子波上的流出信息特征,進(jìn)行癲癇發(fā)作檢測。將所有頻帶的檢測結(jié)果進(jìn)行比較,如圖2所示。通過單因素方差分析,全頻帶分類結(jié)果的5項指標(biāo)與γ波段分類結(jié)果的5項指標(biāo)無顯著性差異(P>0.05)。低頻段的4個子波段分類結(jié)果的5項指標(biāo)與γ波段分類結(jié)果的5項指標(biāo)有著顯著性差異(P<0.05)。這說明了不同子波段對于癲癇腦電的檢測都提供了不同量的特征;另一方面,說明了癲癇發(fā)作與非發(fā)作相比,γ波段對于癲癇腦電的檢測提供更多的特征,即γ波段的信息流強度變化顯著。

3.2 基于部分有向相干(PDC)不同波段癲癇檢測結(jié)果對比

為探究癲癇發(fā)作前后γ波段的信息流動變化是否顯著,本文采用了另一種腦功能連接的PDC算法。相較于DTF算法,PDC算法只會發(fā)現(xiàn)一些直接的因果關(guān)系,比如A、B、C 3者之間的關(guān)系,用PDC可能只能發(fā)現(xiàn)A到C的有向連接,而用DTF可能會發(fā)現(xiàn)A到B、B到C和A到C都有因果性聯(lián)系。因此,本文采用PDC算法進(jìn)行腦功能連接,求出腦電各個頻段上的流出信息強度特征,用以癲癇發(fā)作的檢測,檢測結(jié)果如圖3所示。單因素方差分析的結(jié)果與DTF算法檢測結(jié)果基本一致,其中全頻帶分類結(jié)果的5項指標(biāo)與γ波段分類結(jié)果的5項指標(biāo)無顯著性差異(P>0.05),其他4個子波段分類結(jié)果的5項指標(biāo)與γ波段分類結(jié)果的5項指標(biāo)有著顯著性差異(P<0.05)。

圖3 基于部分有向相干(PDC)算法的不同波段癲癇檢測結(jié)果

3.3 癲癇發(fā)作前后流出信息強度變化

為了探究癲癇發(fā)作前后流出信息強度的變化情況,尋找DTF和PDC算法在γ波段癲癇發(fā)作檢測效果優(yōu)于其他子波段的原因,將癲癇發(fā)作前后各個波段上的各通道的流出信息強度進(jìn)行了對比,并選取了一位病人(患者2),給出了在一次癲癇發(fā)作前后10 s的不同波段上的腦區(qū)流出信息強度的腦地形圖,如圖4所示。癲癇發(fā)作的前后10 s,δ波、θ波、α波、β波、γ波的流出信息強度變化不同。癲癇發(fā)作前10 s,各頻段的流出信息強度也不相同,δ波和θ波流出信息強度主要散布在右側(cè)前額葉、右側(cè)顳葉和枕葉;α波流出信息強度主要散布后側(cè)顳葉和枕葉;β波流出信息強度集中在額葉、頂葉和枕葉;γ波流出信息強度都很小,全頻帶流出信息強度也都很小,少量集中在枕葉。癲癇發(fā)作第10 s,δ波流出信息量主要在后額葉和右側(cè)后顳葉;θ波流出信息強度主要集中在后額葉、頂葉和左側(cè)前額葉;α波、β波、γ波和全頻帶流出信息強度集中的腦區(qū)位置很相似,主要在后額葉和頂葉。全頻帶的流出信息強度變化和γ波段變化很相似,而且它們在癲癇發(fā)作前后流出信息強度變化差異很大。特別在γ波,在癲癇發(fā)作之前,各個腦區(qū)的流出信息強度都比較弱,發(fā)作開始后,出現(xiàn)明顯的信息流在某些腦區(qū)匯聚和加強。這也可能是圖2中用γ和全頻段檢測結(jié)果的5項指標(biāo)無顯著性差異且比其他子頻段的檢測結(jié)果好的原因。

Ayala等用顱內(nèi)腦電(iEEG)的γ波段的功率譜進(jìn)行癲癇檢測取得了不錯的效果,并指出癲癇發(fā)作與γ波的功率譜有一定關(guān)系[17]。Lu等用頭皮腦電的高頻部分對癲癇病灶進(jìn)行定位,也說明了癲癇的發(fā)作與高頻腦電活動有關(guān)[18]。本文實驗結(jié)果顯示:γ波在癲癇發(fā)作前后流出信息強度變化顯著,說明了癲癇發(fā)作與γ波甚至更高的高頻振蕩有著密切的關(guān)系。對于局灶性癲癇患者,在癲癇發(fā)作時高頻成分的變化,導(dǎo)致了高頻波段的信息流在某一腦區(qū)聚集和增強,而且相對規(guī)律,這種特征可以被捕捉到并用于癲癇發(fā)作檢測。相對于其他子波段,在癲癇發(fā)作前后,信息流的聚集和變化情況表現(xiàn)為不規(guī)律狀態(tài),故沒有用高頻成分檢測效果好。

3.4 滑動時間窗對檢測結(jié)果的影響

本文分別在滑動時間窗為2、3和5 s時進(jìn)行了癲癇發(fā)作檢測。當(dāng)窗口為3 s時,檢測結(jié)果的平均正確率為98.1%,平均選擇性為55.2%,平均敏感性為93.4%,平均特異性為98.2%,平均檢出率為95.7%。當(dāng)窗口為5 s時,檢測結(jié)果的平均正確率為98.4%,平均選擇性為61.2%,平均敏感性為96.6%,平均特異性為98.5%,平均檢出率為97.6%。結(jié)果顯示,檢測的各項指標(biāo)無顯著性差異(P>0.05),所以我們認(rèn)為該方法與時間窗無關(guān)。為了減少運算復(fù)雜度和在癲癇發(fā)作后最短的時間檢測到,本文采用2 s的時間窗。

圖4 病人2的癲癇發(fā)作前后各個波段的腦電流出信息強度腦地形圖

3.5 與其他檢測方法比較

圖5 不同檢測方法的結(jié)果對比

選擇AR模型[19]和近似熵[20]兩種現(xiàn)存的癲癇發(fā)作檢測方法在相同數(shù)據(jù)集上進(jìn)行實驗。計算近似熵時時間窗為2 s,然后將每一導(dǎo)的信號求近似熵,將19導(dǎo)得到的特征進(jìn)行組合。對于AR模型方法,時間窗口也是2 s,設(shè)置的階數(shù)為3,并將19導(dǎo)得到的特征進(jìn)行組合。對比結(jié)果如圖5所示,AR模型方法檢測結(jié)果的平均正確率為72.6%,平均選擇性為4.9%,平均敏感性為51.2%,平均特異性為72.9%,平均檢出率為62.1%。近似熵方法檢測結(jié)果的平均正確率為93.8%,平均選擇性為23.9%,平均敏感性為85.5%,平均特異性為93.8%,平均檢出率為89.6%。本文方法檢測結(jié)果的5項指標(biāo)與AR模型和近似熵方法的檢測結(jié)果相比有顯著的提升(P<0.05)。

4 結(jié) 論

本文利用長時程頭皮腦電的高頻成分對局灶性癲癇患者進(jìn)行癲癇發(fā)作檢測,將多通道的腦電信號通過小波變換提取γ波段,然后在移動的時間窗內(nèi)建立MVAR模型,采用有向傳遞函數(shù)算法求出γ波段上不同導(dǎo)聯(lián)的流出信息強度,最后經(jīng)過SVM分類器進(jìn)行分類,識別出發(fā)作期腦電和非發(fā)作期腦電。通過五重交叉驗證得到的平均正確率為98.4%,平均選擇性為60.7%,平均敏感性為93.4%,平均特異性為98.4%,平均檢出率為95.9%。通過和其他波段腦電信號癲癇發(fā)作檢測結(jié)果比較,發(fā)現(xiàn)γ波段的檢測效果最好。另外,通過對癲癇發(fā)作前后不同腦區(qū)在不同子頻段上的流出信息強度進(jìn)行對比與分析,發(fā)現(xiàn)癲癇發(fā)作時在γ波段上信息流在某一腦區(qū)集中和增強,而且相對規(guī)律,這種特征能被本文方法捕捉并可以實現(xiàn)癲癇發(fā)作的檢測。對大腦高頻成分活動的研究可能對癲癇發(fā)作起因、發(fā)作檢測、發(fā)作預(yù)測和病灶定位有著重要意義。

[1] SUBASI A. Epileptic seizure detection using dynamic wavelet network [J]. Expert Systems with Applications, 2005, 29(2): 343-355.

[2] SHARMA R, PACHORI R B. Classification of epileptic seizures in EEG signals based on phase space representation of intrinsic mode functions [J]. Expert Systems with Applications, 2015, 42(3): 1106-1117.

[3] ACHARYA U R, SREE S V, SWAPNA G, et al. Automated EEG analysis of epilepsy: a review [J]. Knowledge-Based Systems, 2013, 45(3): 147-165.

[4] GENG S, ZHOU W. Nonlinear feature comparison of EEG using correlation dimension and approximate entropy [C]∥2010 3rd International Conference on Biomedical Engineering and Informatics. Piscataway, NJ, USA: IEEE, 2010: 978-981.

[5] GULER N, UBEYLI E, GULER I. Recurrent neural networks employing Lyapunov exponents for EEG signals classification [J]. Expert Systems with Applications, 2005, 29(3): 506-514.

[7] KHAN Y U, GOTMAN J. Wavelet based automatic seizure detection in intracerebral electroencephalogram [J]. Clinical Neurophysiology, 2003, 114(5): 898-908.

[8] 馬東華, 鄭旭媛, 王真. 基于形態(tài)成分分析的癲癇腦電棘波檢 [J]. 生物醫(yī)學(xué)工程雜志, 2013, 30(4): 710-713. MA Donghua, ZHENG Xuyuan, WANG Zhen. Detection of epileptic spike wave in EEG signal Based on morphological component analysis [J]. Journal of Biomedical Engineering, 2013, 30(4): 710-713.

[9] VAN MIERLO P, PAPADOPOULOU M, CARRETTE E, et al. Functional brain connectivity from EEG in epilepsy: seizure prediction and epileptogenic focus localization [J]. Progress in Neurobiology, 2014, 121: 19-35.

[10]WANG G, SUN Z, TAO R, et al. Epileptic seizure detection based on partial directed coherence analysis [J]. Journal of Biomedical and Health Informatics, 2016, 20(3): 873-879.

[11]KAMINSKI M J, BLINOWSKA K J. A new method of the description of the information flow in the brain structures [J]. Biological Cybernetics, 1991, 65(3): 203-210.

[12]ADELI H, ZHOU Z, DADMEHR N. Analysis of EEG records in an epileptic patient using wavelet transform [J]. Journal of Neuroscience Methods, 2003, 123(1): 69-87.

[13]AKAIKE H. A new look at the statistical model identification [J]. IEEE Transactions on Automatic Control, 1974, 19(6): 716-723.

[15]吳烜, 李京. 基于支持向量機的測厚儀CS值電壓漂移故障判定及處理 [J]. 計算技術(shù)與自動化, 2014, 33(1): 42-45. WU Xuan, LI Jing. Thickness gauge CS value voltage drift down approach based on support vector machines algorithm determines [J]. Computing Technology and Automation, 2014, 33(1): 42-45.

[16]AARABI A, GREBE R, WALLOIS F. A multistage knowledge-based system for EEG seizure detection in newborn infants [J]. Clinical Neurophysiology, 2007, 118(12): 2781-2797.

[17]AYALA M, CABRERIZO M, JAYAKAR P, et al. Subdural EEG classification into seizure and nonseizure files using neural networks in the gamma frequency band [J]. Journal of Clinical Neurophysiology, 2011, 28(1): 20-29.

[18]LU Y, WORRELL G A, ZHANG H C, et al. Noninvasive imaging of the high frequency brain activity in focal epilepsy patients [J]. IEEE Transactions on Biomedical Engineering, 2014, 61(6): 1660-1667.

[19]CHISCI L, MAVINO A, PERFERI G, et al. Real-time epileptic seizure prediction using AR models and support vector machines [J]. IEEE Transactions on Biomedical Engineering, 2010, 57(5): 1124-1132.

[20]YENTES J M, HUNT N, SCHMID K K, et al. The appropriate use of approximate entropy and sample entropy with short data sets [J]. Annals of Biomedical Engineering, 2013, 41(2): 349-65.