李斯卉, 呂可嘉, 潘 敏, 張 瑞

(1.西北大學 醫(yī)學大數(shù)據(jù)研究中心，陜西 西安 710127； 2.西安交通大學 醫(yī)學部,陜西 西安 710061)

癲癇是最常見的神經(jīng)功能障礙之一,世界范圍內(nèi)的患病率接近0.9%。其最主要的臨床癥狀表現(xiàn)為癲癇性發(fā)作,通常是由大腦內(nèi)大量神經(jīng)細胞群異常超同步放電所引起。腦電圖(Electroencephalogram,EEG)通過追蹤和記錄腦電波來呈現(xiàn)大腦的放電活動,是目前用于檢測癲癇性發(fā)作等神經(jīng)功能障礙異常腦電模式的有效手段。然而,傳統(tǒng)的癲癇診斷往往需要專業(yè)醫(yī)師依據(jù)經(jīng)驗對長時程腦電圖通過視覺上的檢查來完成,這一過程不僅耗時、主觀性強,而且大量噪聲的存在也使得這一工作極具挑戰(zhàn)。因此,近幾十年來,越來越多的研究者開始致力于開展癲癇性發(fā)作自動檢測的研究。

神經(jīng)學家的近期研究結(jié)果表明,人的大腦可以看作是一個非線性動力系統(tǒng)。基于此,各種建立在非線性方法基礎(chǔ)上的癲癇性發(fā)作自動檢測方法被相繼提出。從分析系統(tǒng)混沌性的角度出發(fā),Beyli等人計算腦電信號的Lyapunov指數(shù)作為特征,并應(yīng)用多層感知器神經(jīng)網(wǎng)絡(luò)(MLPNN)來實現(xiàn)癲癇性發(fā)作的自動檢測[26]?；谙到y(tǒng)復(fù)雜性分析,Zhang 等人結(jié)合分形維數(shù)和梯度增強提出了一種癲癇自動檢測方法[28],Song 等人提出了一種基于模糊條件Renyi熵的癲癇性發(fā)作自動檢測方法[24]。除此之外, Zhou 等人提出了一種基于小波變換的癲癇性發(fā)作檢測方法,該方法提取殘差和波動指數(shù)作為特征且應(yīng)用貝葉斯線性判別分析(BLDA)作為分類器[29]。Li等人在文獻[16] 中定義了腦電信號三維相空間表示的核半徑,并以此作為所提取的特征結(jié)合超限學習機和支撐向量機完成發(fā)作的自動檢測[16]。在文獻[7] 中,Chen等人從EEG子帶信號中提取了近似熵、遞歸量化分析和樣本熵三個非線性度量,并將其融合特征作為支撐向量機和超極限學習機的輸入以完成癲癇性發(fā)作的自動檢測。在非線性相互依賴性分析方面,Song 等人提出了基于龐加萊散點圖的特征,通過刻畫延遲龐加萊散點圖的分散程度和分布一致性并結(jié)合ELM完成自動檢測[23]。為了分析系統(tǒng)的高階統(tǒng)計特性,Acharya等人用高階譜(HOS)結(jié)合SVM來區(qū)分正常、發(fā)作期和發(fā)作間期的EEG信號[1]。此外,用模糊相似性指數(shù)[14]和動力學相似性指數(shù)[17]對腦電信號的非線性時空演變進行分析,且被用于分析發(fā)掘從非癲癇狀態(tài)到癲癇發(fā)作狀態(tài)之間的過渡。在文獻[22] 中,Song等人提取了基于馬氏相似度的特征并進一步結(jié)合樣本熵進行特征融合,最終采用超限學習機完成癲癇檢測。

本文提出了一種新的腦電信號相似性的分析方法,并進一步探討了其在癲癇性發(fā)作自動檢測中的應(yīng)用。首先根據(jù)時間延遲方法計算腦電信號的軌跡矩陣(TMS);其次,通過對軌跡矩陣進行奇異值分解(SVD)定義腦電信號的特征表示向量(FRV);進而,基于信號的特征表示向量計算待測腦電與模版腦電之間的最大余弦相似度(MCS),并將其定義為所提取的癲癇腦電特征;最后,將所提取特征輸入到分類器中完成癲癇性發(fā)作的自動檢測。

1 方 法

1.1 相空間和軌跡矩陣

動力系統(tǒng)可以通過一組狀態(tài)以及狀態(tài)間的轉(zhuǎn)換規(guī)則來描述[13]。系統(tǒng)的狀態(tài)通常以向量形式表示,而狀態(tài)向量的集合則定義為系統(tǒng)的相空間。我們可以通過評估其相空間來研究系統(tǒng)的動力學特性[15],然而在大多數(shù)情況下,一個系統(tǒng)的相空間很難直接獲取。基于此,通過所記錄的系統(tǒng)輸出信號來重建該系統(tǒng)相空間的不同重建方法被相繼提出,其中時間延遲法是應(yīng)用最為廣泛的一類[25]。

給定時間序列s={s1,s2,…,sN},其相空間可以由下式表達:

A={A1,…,Aj,…,AL}T

(1)

Aj={sj,sj+τ…,sj+(m-1)τ}T。

(2)

其中m表示嵌入維度,τ為時間延遲,L=N-(m-1)τ。由于重構(gòu)的相空間A以矩陣形式表示,因此A也稱為信號s的軌跡矩陣(trajectory matrix,TM)。需要強調(diào)的是,參數(shù)m和τ的選擇非常重要。本文采用互信息法(MI)確定參數(shù)τ的取值[4],并根據(jù)參考文獻[5] 中的方法確定m的取值。

1.2 EEG信號的特征表示向量

奇異值分解是最具代表性的一類矩陣分解方法,在機器學習的算法中占有重要地位[20]。

給定EEG信號s={s1,s2,…,sN}。根據(jù)式(1)和(2)可得到s的軌跡矩陣

則矩陣A的奇異值分解可表示為

(3)

注1除SVD分解外,常用的矩陣分解方法還包括LU分解,Cholesky分解,QR分解和LDLT分解等。其中Cholesky分解和LDLT分解均假設(shè)矩陣是正定的,QR分解假設(shè)矩陣是方陣且非奇異的。因此,針對本文的研究對象EEG信號,LU分解也可用來獲得EEG信號的特征表示向量。數(shù)值實驗結(jié)果表明,奇異值分解在完成癲癇性發(fā)作檢測方面優(yōu)于LU分解(見表4)。因此,本文所提方法中采用奇異值分解。

1.3 EEG信號間的最大余弦相似度

余弦相似性(cosine similarity,CS)通過在內(nèi)積空間中計算兩個非零向量夾角的余弦值來度量其相似性。給定兩個非零向量x=(x1,x2,…,xn)T和y=(y1,y2,…,yn)T,則x與y之間的余弦相似性定義為:

(4)

其中θ為x與y之間的角度,x·y表示x與y的內(nèi)積,‖·‖表示向量長度。在此基礎(chǔ)上,可進一步將其擴展為度量兩個矩陣之間的相似性。具體地,給定兩個具有相同階數(shù)q×n的矩陣X=(x1,x2,…,xn)q×n和Y=(y1,y2,…,yn)q×n,其中x1,x2,…,xn和y1,y2,…,yn分別是X和Y的列向量。對于每對列向量xi和yj,根據(jù)式(6)計算它們之間的余弦相似性CS(xi,yj),則有

MX,Y=

定義2設(shè)X和Y是兩個同階矩陣,定義

(5)

為X和Y之間的行最大余弦相似度(MCS-R),其中

(6)

同理,定義

(7)

為X和Y之間的列最大余弦相似度(MCS-C),其中

(8)

注2在定義2中,行最大余弦相似度MCS-R和列最大余弦相似度MCS-C顯然滿足MCS(r)(Y,X)=MCS(c)(X,Y)。

注3在定義2中,如果兩個矩陣不滿足同階但具有相同的行數(shù)(或相同的列數(shù)),則可同理定義它們之間的行(列)最大余弦相似度。

結(jié)合定義1和定義2,我們進一步刻畫兩個EEG信號之間的相似性如下。

MCS(s,t)=MCS(r)(Us,Ut)

(9)

為s和t的最大余弦相似度,且

MCS(t,s)=MCS(c)(Us,Ut)

(10)

為t和s的最大余弦相似度。

1.4 基于EEG信號相似性分析的癲癇性發(fā)作自動檢測方法

在上述工作的基礎(chǔ)上,本小節(jié)提出一種新的基于最大余弦相似度的癲癇腦電特征,并結(jié)合有效的分類器以完成癲癇性發(fā)作的自動檢測。該方法總結(jié)為以下TFM-SZ算法。

算法Ⅰ(TFM-SZ算法) 令S={s1,s2,…,sΓ}為腦電信號集合,其中si={si,1,si,2,…,si,N}表示第i個腦電信號,Γ表示S中腦電信號的總數(shù)。

步驟2構(gòu)建參考腦電信號。

1)從S中隨機選取p個發(fā)作間期腦電信號,記為{si1,si2,…,sip},其中p小于S中的發(fā)作間期腦電信號總數(shù);

步驟4計算參考腦電信號與每個待測腦電信號之間的最大余弦相似度。

對于i=1:Γ

End

步驟5定義癲癇腦電特征(TFM-F)為

(11)

步驟6將提取的特征F輸送到某個分類器中。

TFM-SZ算法的流程圖如圖1所示。

圖1 TFM-SZ算法的流程圖Fig.1 The flow chart of TFM-SZ Algorithm

2 數(shù)值試驗結(jié)果與分析

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

本文數(shù)值實驗中使用的腦電數(shù)據(jù)來自兩個數(shù)據(jù)庫:Bonn數(shù)據(jù)庫和CHB-MIT數(shù)據(jù)庫。Bonn數(shù)據(jù)庫的詳細信息如表1所示,圖2展示數(shù)據(jù)集D和數(shù)據(jù)集E的兩個腦電片段(單通道)。CHB-MIT數(shù)據(jù)庫采集自波士頓兒童醫(yī)院,采集自22名患有難治性癲癇的受試者(17名女性,年齡為1.5～19歲;5名男性,年齡為3～22歲)。共包含664個腦電圖記錄,采樣率為256Hz,分辨率為16bit,大多數(shù)腦電信號為23通道(個別記錄為24或26通道)。圖3展示了來自CHB-MIT數(shù)據(jù)庫中“記錄chb03”的一個23通道的腦電片段。

表1 BONN數(shù)據(jù)庫的詳細信息Tab.1 Details of BONN database

注：①BONN數(shù)據(jù)庫取自德國BONN大學癲癇系[2,24];②每個數(shù)據(jù)集包括100個單通道腦電圖,采樣率為173.61Hz,每個腦電圖段持續(xù)23.6s;③所有EGG都沒有造假。

圖2 BONN數(shù)據(jù)庫的EEG片段Fig.2 Sample EEG recordings from Boon database

2.2 結(jié)果與分析

在數(shù)值實驗中,我們從Bonn數(shù)據(jù)庫的數(shù)據(jù)集D中隨機選擇p=30個發(fā)作間期腦電信號,再從每個信號中隨機選取持續(xù)時長為10s的腦電片段(即d=173.6×10=1 736)。在特征值TFM-F的計算過程中,設(shè)定嵌入維數(shù)m=5和延遲時間τ=1。在超限學習機(extreme learning machine,ELM)[10]中,采用可加性隱節(jié)點G(a,b,x)=g(ax+b),隱節(jié)點參數(shù)(a,b)服從均勻分布隨機選取,通過十倍交叉驗證選擇隱節(jié)點個數(shù)為80。在支持向量機(support vector machine,SVM)[9]中,使用最新的Libsvm軟件包3.22版本,核函數(shù)選取徑向基函數(shù)。根據(jù)網(wǎng)格搜索法,正則參數(shù)C和寬度g取為(2-4,2)。所有的實驗均在Matlab 9.1.0中運行。

圖3 CHB-MIT數(shù)據(jù)庫的EEGFig.3 Sample EEG recordings from CHB-MIT databse

2.2.1 Bonn數(shù)據(jù)庫上的性能驗證 本文提取了一個新的基于相似性分析的癲癇腦電特征TFM-F,并將其用于實現(xiàn)發(fā)作期腦電、發(fā)作間期腦電和正常腦電的分類與檢測。

首先，驗證所提取特征TFM-F分別結(jié)合ELM和SVM在完成癲癇性發(fā)作自動檢測中的性能。在執(zhí)行SVM和ELM的過程中均進行50次測試,每次測試隨機生成個數(shù)均衡的訓練集和測試集[24]。表2列出了性能的比較結(jié)果,包括分類準確率(ACC)、標準差(Std)、訓練時間(training time)和參數(shù)選擇時間(parameter selection time)。從表2可看出,ELM相比SVM,分類準確率略高一點,然而其標準差更小,這表明ELM比SVM的分類性能更加穩(wěn)定。同時,ELM的訓練時間僅為SVM的1/430。由于ELM具有較好的分類性能,因此在進一步的數(shù)值實驗中我們均采用ELM作為分類器。

其次，比較兩種矩陣分解方法以及5種相似性度量在特征提取中的性能。表3,4中結(jié)果可展示所提方法的優(yōu)良性能,包括敏感性(SEN),特異性(SPE),精度(PRE)和準確率(ACC)[22]。從表3可以看出,由“TFM(SVD)+ELM” 方法得到的敏感性,特異性,精度和準確率均遠遠高于通過 “TFM(LU)+ELM”方法所得到的值。此外從表4可知,與其他4種相似性度量方法相比,“TFM(Cosine)+ELM”的性能表現(xiàn)最好。上述結(jié)果表明,本文特征提取算法中所采用的SVD分解和余弦相似性最為有效。

表2 “TFM-F+ELM”和“TFM-F+SVM”的性能比較Tab.2 Performance comparison between "TFM-F+ELM" and "TFM-F+SVM" on BONN database

注:運行時間=訓練時間+參數(shù)選擇時間。

表3 特征提取中兩種矩陣分解LU和SVD的性能比較Tab.3 Performance comparisonbetween LU and SVD in extraction the feature TFM-F

表4 特征提取中5種相似計算的性能比較Tab.4 Performance comparison between consine similarity and other four similarity measures in extracting the feature TFM-F

然后對本文所提取特征TFM-F與模糊相似性指數(shù)(FSI)和樣本熵(SE)的性能進行系統(tǒng)比較。從表5可以看出,與FSI和SE相比,本文所提取特征TFM-F的性能有明顯提高。具體來說, TFM-F的漏檢率和誤報率均達到0.15%,與 FSI的結(jié)果相比分別減少了近55倍和45倍。此外,FSI所耗費的訓練時間是本文所提方法的44倍左右。另一方面,雖然SE所耗費的訓練時間只有TFM-F的三分之一,但其漏檢率和誤報率分別比TFM-F升高了近3倍和20倍。此外,本文采用ROC曲線更為直觀地比較上述3個特征的性能(見圖4),顯然TFM-F的AUC(ROC曲線下方包圍的面積)最大,表明其檢測性能最好。

最后將本文所提出的自動檢測方法TFM-SZ和已有6種自動檢測方法進行比較(所有結(jié)果均采用相同數(shù)據(jù)集)。從表7可以看出,本文方法將準確率從97.53%提高至99.85%。

圖4 3個特征TFM-F,FSI和SE的ROC曲線Fig.4 ROC curves for three features TFM-F, FSI and SE

特征ACC/%SEN/%SPE/%PRE/%TFM-F+ELM99.8599.8599.8699.85FSI+ELM93.6491.7593.1092.03SE+ELM96.9499.6294.2994.51

2.2.2 CHB-MIT數(shù)據(jù)庫上的性能驗證 本小節(jié)采用CHB-MIT數(shù)據(jù)庫中的頭皮腦電數(shù)據(jù),驗證所提癲癇性發(fā)作自動檢測方法TFM-SZ的性能。這一數(shù)據(jù)庫中的腦電信號含有大量偽跡并具有嚴重的數(shù)據(jù)不平衡現(xiàn)象,因此會給自動檢測的有效實現(xiàn)帶來更大困難。然而,這樣的挑戰(zhàn)在臨床應(yīng)用中卻具有重大的現(xiàn)實意義。本文選取了來自CHB-MIT數(shù)據(jù)庫中5名患者的腦電數(shù)據(jù)進行數(shù)值實驗。表6中詳細列出了上述數(shù)據(jù)的詳細信息,包括記錄時間、發(fā)作次數(shù)和發(fā)作持續(xù)時間。采用小波軟閾值法去除偽跡,其中選用4階Daubechies小波作為小波母函數(shù)進行5層分解。通過重疊時長為3s的滑動窗口,把所有去噪后的腦電信號分割成5s長的腦電片段并將這些片段隨機分為個數(shù)相同的訓練集和測試集?？紤]到CHB-MIT數(shù)據(jù)庫中發(fā)作數(shù)據(jù)和非發(fā)作數(shù)據(jù)的嚴重不平衡,本文采用敏感性和誤報率(FDR)來評估算法性能,計算公式如下:

敏感性=發(fā)作檢測正確的次數(shù)/發(fā)作的總次數(shù),

誤報率=發(fā)作檢測錯誤的次數(shù)/記錄總時長。

表6中展示的實驗結(jié)果包括檢測正確的次數(shù)、檢測錯誤的次數(shù),敏感性和誤報率?？梢钥闯?只有患者“chb03”出現(xiàn)了一次錯誤檢測,其誤報率為0.0625。此外,在患者“chb08”(女,3.5歲)8h的腦電圖記錄中共出現(xiàn)了兩次發(fā)作,但本文所提方法僅成功地檢測出一次發(fā)作,敏感性為50%,這一結(jié)果很可能是由該患者年齡過小所造成的。

表6 所提自動檢測方法TFM-SZ在CHB-MIT數(shù)據(jù)庫上的性能評估Tab.6 Performance evaluation of the proposed method TFM-SZ on CHB-MIT database

表7 TFM-SZ的建議方法和BONN數(shù)據(jù)庫上的其他現(xiàn)有方法的性能比較Tab.7 Performance comparison among the proposed method TFM-SZ and other existing methods on BOON database

3 結(jié) 論

本文提出了一種新的腦電信號相似性的分析方法,進而在這一方法的基礎(chǔ)上定義了待測腦電與模版腦電之間的最大余弦相似度為癲癇腦電特征,并將其應(yīng)用于癲癇性發(fā)作的自動檢測。文中從三個方面驗證了所提癲癇檢測方法TFM-SZ的性能:①所提取特征TFM-F的性能評估(Bonn數(shù)據(jù)庫),包括兩個矩陣分解方法LU 與SVD的性能比較、余弦相似度與其他四種相似性度量方法的比較、以及3個特征間的性能比較;②所提出癲癇性發(fā)作自動檢測方法TFM-SZ的性能評估(Bonn數(shù)據(jù)庫和CHB-MIT數(shù)據(jù)庫),包括兩個分類器ELM和SVM的性能比較,所提方法TFM-SZ與其他6種已有方法的性能比較,以及TFM-SZ針對原始采集的頭皮腦電數(shù)據(jù)上的檢測性能驗證。所有數(shù)值實驗結(jié)果表明,本文所提出的腦電信號相似性分析方法可以成功應(yīng)用于癲癇性發(fā)作的自動檢測。