腦磁圖腦功能連接網絡癲癇棘波識別方法研究

張航宇，李 彬，尹春麗，劉 凱，王玉平，張軍鵬

1.四川大學 電氣信息學院 自動化系，成都610065

2.首都醫科大學 宣武醫院 神經內科，北京100053

1 引言

癲癇是大腦神經元異常放電，導致大腦功能障礙的一種慢性腦部疾患，以腦部神經元過度放電所致的突然反復和短暫的中樞神經系統功能失常為特征[1]。在我國，癲癇的患病率約為3‰ 至6‰ ，難治性癲癇約占20%[1-2]，目前全世界患有癲癇疾病的人數約占總人口數的1%，且每年以0.02%～0.05%的速度增長[3]。

自發功能連接（spontaneous Functional Connectivity，sFC）已成為認知神經科學的一個主要工具[3]。自發功能連接是基于對腦區間的自發波動信號之間的統計關系的分析，一些研究提出它能反映腦區之間的深層神經互動[4]。作為神經科學研究的一個工具，自發功能連接能夠幫助理解脫離任務環境下的大腦是如何工作的[5]。自發功能連接作為一個重要且廣泛使用的工具，它可以檢測典型大腦的功能連接并記錄神經發育[6]，還可以描繪神經學和精神病學失調病癥的異常腦活動的特性，這些病癥包括中風、帕金森癥、阿爾茲海默癥、癲癇和自閉癥。

自發波動由Biswal 基于功能磁共振數據在血液氧合中第一次觀察到，隨后同樣的網絡在大腦電生理學記錄腦電圖和腦磁圖中也被觀察到[7]。很多研究出于高空間分辨率的原因使用功能磁共振成像，但是它的低時間分辨率（0.5～2 Hz）和血管決定的對比度使它不能直接研究被認為是腦區間信息交流機制的高頻神經皮層活動[8]。除此之外，神經活動的時間尺度是典型的遠遠快于fMRI的記錄能力的[9]。

腦磁圖（Magnetoencephalography，MEG）是一種非侵入式的測量神經集束樹狀活動所引發磁場的方法，它的采樣率超過500 Hz，能夠直接用來研究快速神經活動[10]。為了獲取多樣的神經生理學相關頻段，本研究使用MEG，這是由于它的時間分辨率高。有些研究已經證明了使用皮層區域間相位關系作為功能連接量度的有效性。特別的，相位鎖定值是衡量兩個時間序列間相位同步性的方法，它以前被用于分析MEG 的靜息態連接，來建立腦功能連接網絡。腦磁圖檢查是輔助癲癇診斷的一個非常重要的手段。癲癇在腦磁圖中的主要表現是棘波和尖波，有時對它們不作區分，統稱為癲癇瞬變現象[11]或棘波。棘波和尖波突出于背景活動，波幅較高，周期為20～200 ms。

腦磁信號的分析主要是對大腦異常活動的檢測分析，這些工作目前都是由醫療工作者根據經驗通過對患者腦磁圖的視覺檢測完成的。視覺檢測有許多不利的因素，一般一次檢測要記錄60 min 的數據，從很長的數據中找出棘波非常耗費精力，對分析者的判斷力有很高要求，并且醫生在工作中要接觸大量癲癇病例，在繁重的工作任務下分類結果的準確性無法保證。而且，不同的專家對同一記錄的判斷結果也不盡相同。因此，對異常大腦電生理信號的自動檢測就顯得十分重要[12]。

不論選用何種檢測標準，通常要求棘波檢測結果有較高的正確率，較低的漏檢率和誤檢率。本文提出一種基于Desikan-Killiany 腦區劃分的相位鎖定值腦功能連接網絡全頻段的機器學習棘波檢測方法，可對大腦棘波放電狀態進行客觀的標定，有很大的潛力應用于腦磁圖癲癇棘波的自動識別與標記的輔助工具，降低醫生的勞動強度，提升檢測正確率并降低漏檢率和誤檢率。

2 數據與方法

2.1 數據采集

本文使用的MEG數據來自于首都醫科大學宣武醫院腦磁圖中心做檢查的20名被確診為島葉及島蓋癲癇的患者，均在手術前進行了MEG檢查，其中男11例，女9例；年齡15～52歲，平均28.7歲。

MRI 應用1.5T 或3.0T 進行標準的MRI 掃描，包括橫斷位SE 序列T1W1 及TSE 序列T2W1（層厚5 mm）；垂直于右側海馬長軸的斜冠狀位和平行于海馬長軸的橫斷位液體衰減反轉恢復（FLAIR）序列（層厚5 mm）。

MEG 應用芬蘭Elekta Neuromag Oy 公司生產的NM20215A-G 型306 信道全頭型生物磁儀對患者進行檢測。患者于安靜狀態下于磁屏蔽室記錄60 min自發腦磁圖數據。帶通濾波為0.10～330 Hz，采樣率1 000 Hz。腦磁源成像（MSI）采用標準程序306導全頭型腦磁圖系統進行采集，對發作間期癲癇樣波進行標記、離線分析。腦磁波與病人MRI 圖像進行融合，在MRI 上自動生成MSI圖像。

腦磁圖中心的醫生根據經驗為每位病例挑選出三段含棘波的數據區段，每段各10 s，并標出棘波峰值的時間點；作為對比同時提供每位病例三段正常靜息態的數據區段，每段各10 s，為了保證本次選出的棘波數據與正常靜息態數據的正確性與典型性，特意請不同的醫生進行了交叉檢查確認，確保了本次研究采用的兩種數據均為典型棘波態與典型正常靜息態的腦磁圖數據。

2.2 數據預處理

獲取了醫生用傳統方法挑選出的全部20例病例的MEG數據后，采用基于數學軟件矩陣實驗室（MATLAB）平臺下的MEG 分析軟件Brainstorm（brain recordings analysis tool）進行數據預處理。首先，對MEG 進行濾波得到0.1～500 Hz頻段數據；然后，去除數據中的偽跡，主要包括眼電、眼動等對MEG的干擾。

所有數據均嚴格按照上述方法依次進行處理后，從醫生提供的10 s 有棘波數據段中截取含棘波峰值的數據段2 000 ms，從10 s無棘波數據段中間截取相同長度的數據段；最后共獲得2 000 ms時長的有棘波數據段與無棘波數據段各60段，這些數據是進一步分析的基礎。

2.3 腦功能網絡建立

所有數據均通過基于Uhuntu 平臺的工作站，調用FreeSurfer 軟件包，進行自動化分析處理。分析處理過程分為體積處理流程與表面處理流程兩大部分。體積處理流程包括圖像灰度標準化、不均勻磁場的校正、到Talairaeh 空間的配準、非腦組織的去除和白質（White Matter，WM）與灰質（Grey Matter，GM）的分割等。表面處理流程包括對白質曲面進行三維重建，從白質曲面出發，沿著灰質梯度方向向外膨脹得到灰質外表面曲面。定義灰質曲面和白質曲面間的距離為腦皮層厚度，采用T-average算法計算皮層厚度。灰質外表面向外膨脹，得膨脹曲面。膨脹曲面經過球狀形變后與模板進行高維配準，依據Desikan-Killiany圖譜對皮層進行自動分區。Desikan-Killiany 圖譜將全腦分為70 個腦區（左右半腦各35 個），其中胼胝體沒有灰質厚度，因此最后得到68個腦區（左右半腦各34個）的皮層厚度。對全部數據用Brainstorm軟件中自帶的Dipole方法進行求源處理，接著對每位病例每個波段求源過的數據按照Desikan-Killiany圖譜進行下采樣。

腦網絡最基本和最關鍵的兩個組成元素是節點和邊：節點代表了大腦的各個腦區，邊則能反映各個腦區之間的聯系。Desikan-Killiany模板所劃分的68個腦區被定義為腦網絡中的節點，如圖1 所示，每個腦區是一個節點。

圖1 全腦68個腦區分區圖

腦功能網絡使用相位鎖定值（Phase Locking Value，PLV）[13-14]構建。相位鎖定值是對兩個時間序列信號在一段時間內保持恒定相位分離的傾向的度量。對于電生理學記錄的時間序列之間靜息態相位鎖定的測量，可以表征在給定頻率下這些時間序列之間的相位差隨時間變化的特性。PLV的計算公式如下：

其中，N 是采樣點的個數，θ1和θ2是采樣點n 處的瞬時相位值，相位鎖定值為復數，其模的范圍為0到1，0代表隨機的相位關系，1代表固定的相位關系。相位鎖定是一種無方向的量度方式，因此是對稱的。全腦相位鎖定網絡是在每一對皮層腦區時間序列測量結果之間計算出來的，由于相位鎖定的對稱性，只需要在每一對時間序列測量結果之間計算一個相位鎖定值。

下采樣之后得到每個源的數據，之后對頻率范圍為0.1～500 Hz的全部兩組數據中的所有個體數據段分7個頻段分別求出相位鎖定值（PLV）矩陣。全部7個頻段從F1低頻到F7高頻的對應頻率范圍如表1所示。

表1 頻段對應表

為了有代表性，隨機選取單一頻段下相位鎖定值矩陣作為典型展示，在此選取3號病例alpha頻段下2 000 ms時長數據段的相位鎖定值矩陣展示于圖2中。

圖2 相位鎖定值矩陣

3 模型設計

本文使用了三種傳統的分類算法，分別是線性邏輯回歸分類、基于線性核函數和徑向基核函數的支持向量機分類與高斯樸素貝葉斯分類算法。以上模型均使用了scikit-learn框架下的函數實現。

3.1 線性分類原理

本文使用Sigmoid函數作為線性邏輯回歸（Logistic-Regression）分類器函數。Sigmoid函數是使用范圍最廣的一類分類函數，該函數定義為：

其中，y 是樣本被判別為正例的概率，1-y 為樣本被判別為負例的概率。為保證模型的訓練精度，損失函數引入L1正則化項，其中設定正則化系數倒數C=1.0，模型預測值為二分類的概率值，其損失函數為：

計算損失函數取最小值時的參數(w,b)，即能得到最優的分類模型。

3.2 支持向量機分類原理

支持向量機在高維度或無窮維度空間中，構建一個超平面或者一系列的超平面，借助超平面去實現分割，使正負樣本之間的空白最大。該模型對非線性關系分類有較好的效果。基本原理如下。

給定訓練向量xi∈RP,i=1,2,…,n 和一個標簽向量y ∈{1,-1}n。在樣本空間存在一個超平面wTx+b=0 可以實現對樣本進行有效分類，則有：

使得式（4）等號成立的距離超平面最近的幾個點構成支持向量，正反方向支持向量到超平面的距離之和為使得該距離最小的參數即為最優參數。用作為式目標函數，按照拉格朗日乘子法進行優化，簡化為對偶問題：

當原始空間超平面不能對數據進行有效分隔時，需要將數據映射到高維空間，即。在對偶問題中，尋找核函數使得，以便在較低維度空間實現內積計算。本文采用線性核函數（kernel：linear）和徑向基核函數（kernel：rbf）。

3.3 樸素貝葉斯分類

樸素貝葉斯分類（Naive Bayesian Classifier，NBC）由于計算高效、精確度高，并具有堅實的理論基礎而得到廣泛的應用。一般情況下在貝葉斯分類中的所有屬性都直接或間接地發揮作用，即所有的屬性都參與分類，而不是一個或幾個屬性決定分類[15]。貝葉斯分類器的分類原理是通過某標簽的先驗概率，利用貝葉斯公式計算出其后驗概率，選擇具有最大后驗概率的類作為該對象所屬的類。簡單地假設任意兩對特征之間相互獨立，假設每個樣本為獨立同分布，根據中心極限定理，標準化之后的每個特征的數據滿足高斯分布的假設，即P(xi|y)～Gaussian，根據貝葉斯定理，則有：

根據每對特征之間相互獨立的假設，從而有：

根據式（7）即可計算y 的后驗分布。

3.4 原始特征構建

在研究中，先使用作為原始特征空間的PLV復數矩陣數據本身（plv）、該矩陣數據的實部（plv_real）、虛部（plv_imag）、幅角（plv_angle）和模（plv_abs）分別進行了測試。結果表明，在使用的所有模型中，模數據表現最為優異。以樸素貝葉斯分類器模型為例，采用五項數據分類的準確率如表2所示。

表2 基于樸素貝葉斯分類器模型的準確率分布表

經測試，在本文使用的其他四種分類器模型中以上五種數據均有類似表現，因此本文所提到的數據均為PLV復數矩陣的模（plv_abs）數據，選定該原始特征空間復數矩陣的模數據作為原始特征數據集。

在實驗中使用以上四種分類器對plv_abs數據分別進行學習分類，實驗結果如表3所示。

表3 基于原始數據的不同模型分類準確率和AUC結果比較

4 實驗結果與分析

4.1 實驗流程

首先根據120×15 946的原始特征數據集作為輸入，然后對輸入數據集做了標準化，最后對標準化后的數據集應用卡方檢驗、F檢驗與迭代特征消除三種方法進行了特征提取；應用四種分類器分別對上述數據進行分類實驗，對全部實驗結果進行了分析與比較。原始特征數據集構建之后的處理流程如圖3所示。

圖3 實驗流程圖

4.2 實驗結果

為定性地比較各種不同的方法，本文對比了四種分類器的受試者工作特征曲線（Receiver Operating Characteristic Curve，ROC）圖。ROC 曲線是反映真陽性和假陽性的綜合指標。本文的四種分類器對結果的預測均為一個概率值，其中有棘波的標注為正例，正常無棘波狀態標注為負例。按照判別為正例的概率值對測試結果進行排序，再以不同概率值作為正反例判別的閾值，并計算相應閾值下的真陽性率和假陽性率，即可繪制得到分類器的ROC曲線。一般根據對分類任務的不同要求，根據曲線對分類器性能的評判存在不同的標準，且本文的分類器ROC曲線存在交叉，因此本文按照通用的AUC（Area Under the ROC Curve）方式對分類器性能進行評價。AUC是ROC曲線與軸所包絡部分的面積，面積越大，表明該分類器的分類性能越好。

4.2.1 原始特征數據集的實驗比較

使用plv_abs 數據作為原始特征數據集輸入，在對四種個分類模型進行訓練后，實驗結果如表3與圖4所示。

圖4 原始數據的ROC圖

從表3與圖4中可以看出使用樸素貝葉斯分類器模型準確率為0.875，比其他三種模型都要高，其AUC 值為0.914，也是所有分類器中第二高的，說明對原始數據而言使用線性核函數的支持向量機分類器性能最好。

4.2.2 特征標準化的實驗比較

對原始特征數據集的各列特征向量數據均按照公式（8）進行了標準化處理：

其中，xi表示第i 個特征向量，σi表示該列特征的標準差，μi表示第i 個特征的均值，xstd表示標準化之后的特征向量。

使用標準化后數據作為輸入，在對四種模型進行訓練后，測試結果見表4和圖5。

在對原始訓練集的各列數據進行標準化處理后，線性邏輯回歸分類器的準確率從0.771提升為0.917，使用線性核函數的支持向量機分類器的準確率從0.833提升為0.917，而徑向基核函數支持向量機分類器的準確率從0.500 提升為0.938，說明提升最高，樸素貝葉斯分類器的準確率也從0.875小幅提升為0.896，說明對原始訓練集的各列數據進行標準化處理有利于提升準確率，并且四種分類器就準確率而言性能相差不大，結合作為分類器性能評價標準的AUC 值，徑向基核函數支持向量機分類器的AUC 值為0.951，在四種分類器中值最大，故綜合上述結果，使用標準化后的特征數據集作為輸入，選用徑向基核函數支持向量機分類器進行棘波自動分類的效果最好，成功率為93.8%。

表4 基于標準化數據的不同模型分類準確率和AUC結果比較

圖5 原始數據標準化的ROC圖

4.2.3 特征提取

在現實的人類識別中，根據經驗知識能夠輕易識別棘波的存在及相關影響，很明顯未能短時間速實現對15 946 維數據的觀測和判斷，因此推斷，數據的維度中必然存在大量的冗余，因此在適當容忍模型性能降低的基礎上進行了特征選擇，以降低模型的復雜度。

（1）基于單變量的特征選擇

基于單變量的統計測試，通過過濾的方式選擇特征。本文采用卡方檢驗和F檢驗方法，逐一對單個特征與實驗標簽進行相關性檢驗，并根據相關性排序選擇部分特征。基于對提取后的特征多次測試，判定選擇32個特征，即整體特征數的0.2%時就能達到較高的分類準確率又能大幅度降低特征數量。F 檢驗和卡方檢驗所提取的32個特征并不一致。

（2）基于模型的循環特征消除（RFE）

本文采用Logistic回歸分類模型，并使用了L1進行正則化，選擇了32個特征。

4.2.4 特征選擇后數據的實驗比較

分別使用3 種特征選擇方法篩選后的32 個特征向量作為輸入，在對4 種分類模型進行訓練后，實驗結果如表5和圖6～8所示。

表5 不同特征選擇方式下分類器準確率比較

圖6 基于卡方檢驗提取特征的學習分類ROC比較

圖7 基于F檢驗提取特征的學習分類ROC比較

圖8 基于迭代特征消除后特征的學習分類ROC比較

從表5 與表6 中可以看出，經過迭代特征消除方法進行特征選擇后，在四種分類器下的分類準確率與分類器的AUC 值均低于其余兩種選擇方法，說明此種特征選擇方法的效果不好，故重點關注卡方檢驗和F檢驗特征提取選擇方法。

與特征標準化后的實驗結果比較，線性邏輯回歸分類器的準確率從0.917 降低為0.833 與0.833，線性核函數支持向量機分類器的準確率從0.917 降低為0.896 與0.792，而徑向基核函數支持向量機分類器的準確率從0.938降低為0.833與0.854，說明經過特征選擇后，使用上述三種分類器分類，準確率會有所降低，而對于樸素貝葉斯分類器，使用卡方檢驗和F檢驗特征提取后的數據作為輸入，其準確率從0.896 提升為0.917，分類器AUC 值從0.912 提升到0.951，比其余三種分類器均要高，說明針對卡方檢驗和F 檢驗特征提取，樸素貝葉斯分類器性能最優，分類效果最好。

表6 不同特征選擇方式下各模型的AUC比較

5 結果討論

5.1 模型選擇

在本文的研究中，通過醫生的精心挑選得到正常靜息態與棘波態腦磁圖樣本數據，數據被適當地預處理后，分七個頻段經相位同步的PLV算法計算，得到68個腦區之間的PLV復數值矩陣作為特征空間，取該特征空間的模值數據作為原始特征數據集，然后對原始特征數據集做了標準化，最后對標準化后的數據集應用卡方檢驗、F 檢驗與迭代特征消除三種方法進行特征提取，最終構成了五種用于輸入分類器的數據集；本研究中，應用了機器學習的方法，包括線性邏輯回歸、線性核函數支持向量機、徑向基核函數支持向量機與樸素貝葉斯分類算法，構建了四種對應的分類器，將五種輸入數據集依次使用四種分類器進行分類實驗，得到相應的分類準確率，對分類器本身使用繪制受試者工作特征曲線（ROC）圖并計算對應的AUC值方式對分類器性能進行評價，最終得到，在應用標準化后特征數據集不進行特征選擇提取時，徑向基核函數支持向量機分類器的準確率最高，為93.8%，AUC 值為0.951 所以性能最優；應用卡方檢驗、F 檢驗進行特征提取后的特征數量大為減少，對應此種精簡特征數據集時，樸素貝葉斯分類器的準確率雖與完整特征數據集下最高分類準確率相比稍有降低，也可達到91.7%，AUC值為0.951，故性能與完整特征數據集下的徑向基核函數支持向量機分類器相當。

5.2 基于特征選擇的腦區識別

根據特征選擇的結果，32個特征基本能夠實現對數據的較好分類。但基于F 檢驗和卡方檢驗提取的特征存在較大不同。按照68×68 矩形區域的上三角部逐行編號，兩類檢驗的提取結果如表7。

從表7 可以得到，卡方檢驗和F 檢驗提取特征的分布上存在18個重疊特征，因此可以斷定，這18個特征對判別是否是癲癇癥患者具有重要影響。

表7 卡方檢驗和F檢驗提取特征編號的分布

6 結論和建議

對本文列出的18 個關鍵特征位置還原計算，得到所處的位置分別為：（1，20）（1，22）（1，60）（9，13）（12，20）（15，31）（15，65）（18，33）（18，48）（18，52）（25，42）（28，38）（30，61）（30，62）（31，38）（32，54）（55，66）（59，65），18個特征全部處于delta頻段。因此，可以斷定，決定癲癇的腦電信號可以完全由delta頻段的信號完成判別。建議在后續醫療實踐中，對以上18 個特征所對應的18 對大腦皮層區域之間的信號傳遞加以重點關注，下一步從臨床實踐上對以上腦區開展研究。