李淑芳 周衛(wèi)東 袁 琦 蔡冬梅

(山東大學(xué)信息科學(xué)與工程學(xué)院，濟南 250100)

引言

癲癇是以反復(fù)發(fā)作為特征的慢性神經(jīng)系統(tǒng)疾病或綜合征，其發(fā)作是腦內(nèi)神經(jīng)元陣發(fā)性異常電活動的臨床表現(xiàn)[1]。由于目前治療方法的局限，一些頑固性癲癇患者得不到有效的治療。若能在癲癇發(fā)作前預(yù)測其即將發(fā)生，即使是較短的時間，也可使患者或醫(yī)生能夠及時采取必要的預(yù)防保護措施，從而降低癲癇發(fā)作造成的損害。此外，對癲癇預(yù)測的研究也將有助于開發(fā)植入式癲癇控制芯片等新的治療手段，也有助于更清楚地認(rèn)識癲癇發(fā)作的機制。

越來越多的研究表明，癲癇發(fā)作過程不是突發(fā)的，而是一個隨時間演化的過程。在發(fā)作前，患者的某些生理指征已逐步發(fā)生改變，這一階段可能持續(xù)幾分鐘甚至幾小時，被稱為過渡階段或發(fā)作前期。由此，人們逐漸認(rèn)識到，癲癇發(fā)作是有可能被提前預(yù)知的。同時，隨著生理信號記錄手段及信號處理方法的不斷發(fā)展和完善，這些發(fā)作前的改變有可能被記錄并提取出來，從而實現(xiàn)癲癇發(fā)作的預(yù)測[2]。

目前，癲癇預(yù)測算法(seizure prediction algorithm，SPA)主要分為線性方法和非線性方法。線性方法應(yīng)用較早，主要是基于頻域分析。Salant等采用了基于參數(shù)模型的多變量頻譜估計的方法，檢測癲癇發(fā)作前期腦電(EEG)的變化，從而達到預(yù)測的目的[2];Litt在 2001年提出了累積能量法(accumulated energy)[3]。此后，Gigola 等用小波分析方法，計算累積能量的演變來進行癲癇預(yù)測[4]。由于大腦神經(jīng)電活動具有混沌特性，因此許多非線性方法逐漸被用于癲癇預(yù)測，如關(guān)聯(lián)維數(shù)[5]、相似性 指 數(shù)[6]、邊 際 可 測 性 方 法 (marginal predict ability)[7-8]、最大李雅普諾夫指數(shù)[9]等。

癲癇樣放電是癲癇發(fā)作的病理生理學(xué)基礎(chǔ)，臨床上常將棘波、尖波、棘慢復(fù)合波、尖慢復(fù)合波、多棘慢復(fù)合波等陣發(fā)性異常稱為癲癇樣放電。在臨床 EEG檢查中，最重要的是識別EEG中是否出現(xiàn)棘波和尖波。棘波時限為20～70 ms，多數(shù)波幅大于100 μV。在過去的研究中，棘波檢測已經(jīng)被用于癲癇腦電信號的檢測[10-11]，但還未有將棘波頻次變化用于癲癇發(fā)作預(yù)測的研究。為此，分析了癲癇發(fā)作前期腦電活動中棘波頻次的變化規(guī)律，提出并研究了一種新的基于腦電棘波頻次的癲癇預(yù)測算法。對21例癲癇病人含有87次發(fā)作的顱內(nèi)腦電數(shù)據(jù)進行實驗，評估了所提癲癇發(fā)作預(yù)測算法的預(yù)測性能。

1 方法

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

形態(tài)學(xué)濾波器是一種非線性濾波器，它基于信號的幾何特征，利用預(yù)先定義的結(jié)構(gòu)元素對信號進行匹配，提取出與其相似的信號特征。它可以將含有復(fù)雜成分的信號分解為具有物理意義的各個部分，使信號與背景分離，并保持其全局或局部的主要形態(tài)特征?；拘螒B(tài)學(xué)變換的基本運算包括腐蝕、膨脹、開運算和閉運算[12-14]，其中腐蝕和膨脹運算公式分別為

式中，f(n)和g(n)分別是定義在 F={0，1，…，N-1}和 G={0，1，…，M-1}上的離散函數(shù)，且 N?M。f(n)是輸入的時間序列，g(n)是結(jié)構(gòu)元素。f(n)關(guān)于結(jié)構(gòu)元素g(n)的形態(tài)開運算和閉運算分別定義為

Maragos提出了開、閉運算的級聯(lián)組合，構(gòu)造出形態(tài)開-閉和形態(tài)閉-開濾波器，可以同時去除信號中的正、負向脈沖。f(n)關(guān)于結(jié)構(gòu)元素g(n)的形態(tài)開-閉(OC)和閉-開(CO)濾波器分別定義為

采用開-閉、閉-開兩種濾波器的平均組合形式來解決幅度偏倚問題。信號f(n)經(jīng)開-閉、閉-開的平均組合濾波后，其輸出信號y(k)為

1.2 腦電棘波檢測

利用結(jié)構(gòu)元素對癲癇腦電信號進行形態(tài)開-閉、閉-開的平均組合濾波運算，從原始 EEG中檢測出棘波，具體步驟為：首先，對一段癲癇 EEG信號x(k)，利用式(5)和式(6)分別作形態(tài)開-閉和閉-開運算，去除 EEG中的瞬態(tài)信號;其次，利用式(7)得到EEG背景信號y(k);最后，將癲癇信號 x(k)與y(k)相減，得到棘波信號z(k)。

1.3 癲癇發(fā)作預(yù)測算法

對癲癇腦電數(shù)據(jù)分析發(fā)現(xiàn)，癲癇發(fā)作前棘波、多棘波、棘慢復(fù)合波等特征波會增加，為此提出利用棘波出現(xiàn)的頻次來預(yù)測癲癇發(fā)作，癲癇發(fā)作預(yù)測算法包括以下步驟。

步驟1：用帶通濾波器對原始 EEG進行濾波，以去掉肌電信號等高頻干擾，特別是50 Hz的工頻干擾。

步驟2：用移動窗技術(shù)對腦電數(shù)據(jù)逐段分析，每段長度取為5 s，移動步長與窗長相等。利用形態(tài)學(xué)濾波器檢測每段EEG信號的棘波，然后統(tǒng)計本段信號中的棘波數(shù)目。

步驟3：計算累計棘波

式中，num(i)代表第i段信號的棘波數(shù)目。

步驟4：計算棘波頻次SR(Spike Rate)

步驟5：使用平均濾波器平滑 SR，得到 SRm。

步驟6：如果SRm大于閾值，則發(fā)出預(yù)報。閾值的計算方法如下：

1)對患者j次癲癇發(fā)作的腦電信號，找出每次發(fā)作時的 SRm的最大值 max1，max2，…，maxj;

2)比較 max1，max2，…，maxj，找出這 j個值中的最小值min;

3)閾值thres=cons·min，其中 cons為一常數(shù)。

2 實驗及結(jié)果

2.1 腦電數(shù)據(jù)來源

實驗所用癲癇腦電數(shù)據(jù)來自德國Freiburg醫(yī)學(xué)院癲癇研究中心的腦電數(shù)據(jù)庫，數(shù)據(jù)全部是在難治性癲癇患者術(shù)前監(jiān)護中得到的顱內(nèi)腦電信號[15]。癲癇患者的腦電數(shù)據(jù)可分為4個階段，分別為發(fā)作前期(preictal)、發(fā)作期(ictal)、發(fā)作后期(postictal)以及發(fā)作間期(interictal)。數(shù)據(jù)記錄了21例癲癇患者6個導(dǎo)聯(lián)的腦電，總共有87次發(fā)作，509 h的發(fā)作間期，以及73 h的發(fā)作前期。在每一個患者的腦電記錄中，都至少包括50 min的發(fā)作前期數(shù)據(jù)以及足夠長的發(fā)作間期數(shù)據(jù)。對于每個患者，有3個顱內(nèi)電極位于或接近發(fā)作區(qū)。

2.2 結(jié)果及分析

對癲癇腦電信號進行形態(tài)學(xué)濾波，實現(xiàn)背景腦電和棘波分離。選取三角形結(jié)構(gòu)元素，有

由于棘波的周期在20～70 ms之間，腦電數(shù)據(jù)的采樣頻率為256 Hz，L取為5～18。如圖1顯示了一段含有棘波的癲癇發(fā)作時的腦電信號，對其做開-閉、閉-開平均組合濾波后的結(jié)果?？梢钥闯?，棘波可以很好地從癲癇發(fā)作信號中分離出來。

對德國Freiburg醫(yī)學(xué)院癲癇研究中心腦電數(shù)據(jù)庫的數(shù)據(jù)進行實驗。分析了87次發(fā)作的發(fā)作前期、發(fā)作期及發(fā)作后期EEG信號的棘波頻次變化，發(fā)現(xiàn)其中有51次在發(fā)作前幾分鐘甚至1 h出現(xiàn)上升趨勢。圖2表示平滑棘波頻次SRm的變化，當(dāng)SRm大于閾值時則報警。圖2中的水平實線為閾值，灰色垂直條代表癲癇發(fā)作期?？梢钥闯觯诎l(fā)作前期棘波數(shù)目有明顯增加。圖2(a)提前12 min預(yù)測到癲癇發(fā)作，圖2(b)提前4.5 min預(yù)測到癲癇發(fā)作，發(fā)出報警。

圖1 腦電信號棘波檢測。(a)癲癇腦電;(b)背景信號;(c)檢測出的棘波Fig.1 EEG spike detection.(a)Ictal EEG;(b)Background signal;(c)Spikes detected

圖2 平滑棘波頻次變化。(a)3號患者的平滑棘波頻次變化;(b)15號患者的平滑棘波頻次變化Fig.2 The changes of smoothed spike rate.(a)the change of smoothed spike rate for patient 3;(b)the change of smoothed spike rate for patient 15

在癲癇發(fā)作預(yù)測中，靈敏度、錯誤預(yù)測率和預(yù)測時間決定了預(yù)測方法的性能。表1列出了每個患者的靈敏度、錯誤預(yù)測率、最長預(yù)測時間、最短預(yù)測時間以及平均預(yù)測時間。從表中可知，本算法可以對21例癲癇患者中的20例做出癲癇預(yù)判。對全部的87次癲癇發(fā)作中準(zhǔn)確預(yù)測出65次，總預(yù)判率為74.7%，漏檢率為25.3%。對共計496 h的腦電數(shù)據(jù)，出現(xiàn) 55次誤判，誤判率為 0.111/h，特異性為88.9%。其中，有7名患者的靈敏度達到了1，6名患者的誤判率為0，3名患者的靈敏度達到1，同時誤判率為0。此外，最長預(yù)警時間可達到60 min，最短也有幾分鐘，平均預(yù)測時間為26.15 min。

在表1中，可見對2號、6號和8號患者的癲癇預(yù)測的靈敏度偏低。原因在于，雖然2號和6號患者發(fā)作期的SR顯著大于間歇期和發(fā)作前期的SR，但在發(fā)作前期沒有出現(xiàn)SR漸增的現(xiàn)象;而8號患者在發(fā)作期棘波沒有顯著增多，且發(fā)作前期沒有出現(xiàn)SR漸增的現(xiàn)象。

表1 本預(yù)測算法對21例癲癇患者的預(yù)測結(jié)果Tab.1 The prediction results of 21 patients using the proposed algorithm

表2給出本方法與相似指數(shù)、有效相關(guān)維、累積能量[15]和小波相干[16]等4種典型的癲癇預(yù)測方法所得到的發(fā)作預(yù)測性能的比較，實驗數(shù)據(jù)均來自德國Freiburg 21人癲癇腦電數(shù)據(jù)庫。從表2中可以看出，本方法具有更高的靈敏度，且其錯誤預(yù)測率低于相似指數(shù)、有效相關(guān)維和累積能量等3種方法的錯誤預(yù)測率。文獻[16]中用小波相干方法獲得的錯誤預(yù)測率較低，這是由于只采用了預(yù)測效果較好的15個患者的數(shù)據(jù)，而舍棄了其余6個患者的數(shù)據(jù)。

表2 本方法與其他方法的比較Tab.2 Comparison of our approach with other approaches

3 討論和結(jié)論

對癲癇發(fā)作進行預(yù)測和預(yù)警，對于建立新的治療方法、改善患者的生活質(zhì)量和降低癲癇損害，具有重要的意義。在臨床癲癇檢查中，最重要的是觀察EEG中是否出現(xiàn)棘波和尖波。在間歇期，大腦處于相對穩(wěn)定狀態(tài)，幾乎沒有癲癇樣放電;隨著癲癇發(fā)作的臨近，異常放電神經(jīng)元數(shù)目增加，棘波、尖波等異常波逐漸增多。對21例癲癇患者長程顱內(nèi)腦電進行分析，發(fā)現(xiàn)發(fā)作間期棘波數(shù)目較少，隨著發(fā)作的臨近，棘波出現(xiàn)頻繁，棘波頻次增多。這表明許多癲癇發(fā)作是一個隨時間演化的過程，可持續(xù)幾分鐘甚至幾小時。對這一過渡階段或發(fā)作前期的檢測，可用于癲癇發(fā)作的預(yù)測。將棘波頻次作為預(yù)測癲癇發(fā)作的特征，有較強的病理生理學(xué)基礎(chǔ)，且計算量小，并取得了較好的預(yù)測效果。本研究僅使用了顱內(nèi)腦電數(shù)據(jù)進行實驗，在以后的工作中將進一步研究頭皮腦電用于癲癇預(yù)測的性能，以及腦電在發(fā)作前期的改變狀況。

本研究提出并評估了一種基于棘波頻次預(yù)測癲癇發(fā)作的新算法。在實驗中，首先對數(shù)據(jù)進行分段，然后逐段采用形態(tài)學(xué)濾波器檢測棘波頻次(spike rate，SR)，最后用平滑后的 SR做出預(yù)判。實驗結(jié)果表明，所提出的算法有較高的敏感度和較低的誤判率。采用本方法對21例癲癇患者長程顱內(nèi)腦電進行癲癇預(yù)測，預(yù)測準(zhǔn)確率達到了74.7%，每小時錯誤預(yù)測次數(shù)僅為0.111次，平均預(yù)警時間為26.15 min。與相似指數(shù)、有效相關(guān)維、累積能量和小波相干4種癲癇預(yù)測方法相比，本方法具有更好的癲癇預(yù)測性能。

[1]劉曉燕.臨床腦電圖學(xué)[M].北京：人民衛(wèi)生出版社，2006：126-127.

[2]Salant Y，Gath I，Henriksen O.Prediction of epileptic seizures from two-channel EEG[J].Medical and Biological Engineering and Computing，1998，36(5)：549-556.

[3]Litt B，Esteller R，E chauz J，et al.Epileptic seizures may begin hours in advance of clinical on set：a report of five patients[J].Neuron，2001，30(1)：51-64.

[4]Gigola S，Ortiz F，D’Attellis CE，et al.Prediction of epileptic seizures using accumulated energy in a multiresolution framework[J].Journal of Neuroscience Methods，2004，138(1-2)：107-111.

[5]Lehnertz K，Elger CE.Can epileptic seizures be predicted?Evidence from nonlinear time series analysis of brain electrical activity[J].Physical Review Letters，1998，80(22)：5019-5023.

[6]Navarro V，Martinerie J，Le Van Quyen M，et al.Seizure anticipation in human neocortical partial epilepsy[J].Brain，2002，125(3)：640-655.

[7]Drury I，Smith B，Li Dingzhou，et al.Seizure prediction using scalp Electroencephalogram [J]. Experimental Neurology，2003，184(1)：9-18.

[8]Li Dingzhou，Zhou Weiping，Drury I，et al.Non-linear，noninvasive method for seizure anticipation in focal epilepsy[J].Mathematical Biosciences，2003，186(1)：63-77.

[9]Iasemidis LD，Shiau DS，Pardalos PM， et al. Long-term prospective on-line real-time seizure prediction[J].Clinical Neurophysiology，2005，116(3)：532-544.

[10]Dingle AA，Jones RD，Carroll GJ，et al.A multistage system to detect epileptiform activity in the EEG[J].IEEE Trans Biomed Eng，1993，40(12)：1260-1268.

[11]Wilson，Emerson R.Spike detection：a review and comparison of algorithms[J].Clinical Neurophysiology，2002，113(12)：1873-1881.

[12]王晶，徐光華，張四聰，等.基于形態(tài)學(xué)濾波器的棘波提取技術(shù)[J].中國生物醫(yī)學(xué)工程學(xué)報，2007，26(1)：69-73.

[13]陳永利，段會龍.數(shù)學(xué)形態(tài)學(xué)在心電信號處理中的應(yīng)用[J].中國醫(yī)療器械雜志，2006，30(6)：434-436.

[14]季虎，孫即祥，毛玲.基于小波變換與形態(tài)學(xué)運算的ECG自適應(yīng)濾波算法[J].信號處理，2006，3(22)：333-337.

[15]Thomas M，Matthias W，Richard A，et al.Comparison of three nonlinear seizure prediction methods by means of the seizure prediction characteristic[J].Physica D：nonlinear phenomena，2004，194(3-4)：357-368.

[16]Mirowski P，Madhavan D，Lecun Y，et al.Classification of Patterns of EEG Synchronization for Seizure Prediction [J].Clinical Neurophysiology，2009，120(11)：1927-1940.