基于fMRI功能網絡和貝葉斯矩陣分解的腦電源成像方法

劉 柯 楊 東 鄧 欣

(重慶郵電大學計算機科學與技術學院 重慶 400065)

1 引言

腦電(ElectroEncephaloGraphy, EEG)通過在頭皮利用電極檢測皮層神經元同步活動產生的電信號，以探測皮層神經活動，具有無創性、時間分辨率高(毫秒級)等特點。通過頭皮EEG信號估計皮層神經電活動稱為腦電源成像(EEG Source Imaging,ESI)。ESI在認知神經領域和臨床應用中具有重要的作用，皮層活動無創定位可用于診斷與大腦相關的生理、心理功能異常，例如癲癇、精神分裂和腦腫瘤等[1]。

由于從頭皮傳感器到皮層神經活動的映射不唯一，存在無窮多組皮層源信號滿足測量的頭皮腦電圖。因此ESI是具有無窮解高度欠定的逆問題[2]。為重構皮層源信號，常常將大腦皮層分割成若干小三角塊，每個三角塊為一個偶極子或源，源信號和EEG信號之間的關系利用一個線性方程表示。源信號重構通過求解一個線性逆問題實現[3]。然而，由于偶極子數量遠超頭皮傳感器，該線性方程高度欠定，需要先驗假設約束解空間，以獲得唯一解[1]。

一種經典的ESI算法是基于L2范數約束的最小范數估計(Minimum Norm Estimate, MNE)[4]。MNE在所有可能的源中選擇一個能量最小的(以L2范數度量)作為源信號估計。由于皮層表面的源信號更接近頭皮傳感器，更容易被傳感器檢測到，因此MNE偏向于表面源。為彌補深度偏差，主要有兩種方法。第1種方法是加權最小范數估計(weighted Minimum Norm Estimate, wMNE)[5]，wMNE使用導聯矩陣列的范數對不同位置的偶極子進行加權。第2種方法是在MNE的解基礎上對估計源信號的方差信息進行歸一化，代表算法有動態統計參數成像[6]和標準化低分辨率腦電磁成像[7]。盡管基于L2范數約束的ESI算法計算簡便，但其結果都過于彌散，覆蓋了皮層的大部分區域。

為估計皮層活動的位置和尺寸信息，多重稀疏先驗(Multiple Sparse Prior, MSP)[8]在經驗貝葉斯框架下，對全腦進行等間隔采樣以構建空間先驗協方差，并根據EEG數據自動選擇有效的空間先驗。文獻[9]在層級貝葉斯框架下，利用數據驅動分塊(Data Driven Parcelization, DDP)技術將皮層分割為若干個互不重合的團塊以構建空間協方差。為消除EEG信號中頭動、眼動等影響，文獻[10]利用L1范數表示EEG信號擬合殘差，同時利用結構化稀疏約束，提出了L1范數殘差與結構化稀疏源成像(L1-norm Residual and Structured Sparsity based Source Image, L1R-SSSI)算法，獲得局部平滑和全局稀疏的源信號估計。功能磁共振成像(functional Magnetic Resonance Imaging, fMRI)可以精確定位大腦活動，但時間分辨率不足[11]。越來越多的研究試圖融合EEG-fMRI以克服單模態腦功能信號的不足[3,11]。本文主要討論基于fMRI空間約束的EEG源成像的非對稱融合。其中fwMNE(fMRI-weighted Minimum-Norm Estimate)[12]是最常見的基于fMRI約束的ESI算法。fwMNE利用fMRI激活圖約束EEG源活動。然而在EEG與fMRI耦合關系還不明確的情況下，將fMRI的激活區域作為真實源具有很大風險。參數經驗貝葉斯(Parametric Empirical Bayesian, PEB)放寬了這一嚴苛的假設[13,14]。PEB將fMRI信息作為先驗信息引入，并通過由EEG數據確定的超參數衡量其相對權重。文獻[14]利用獨立成分分析(Independent Component Analysis,ICA)提取fMRI數據中的腦功能網絡并構建協方差先驗，提出了網絡源成像(NEtwork based SOurce Imaging, NESOI)。這些研究表明fMRI信息能有效提高EEG源信號重構性能。

以上方法都是基于空間信息對源進行約束。為利用EEG信號中的時域信號，文獻[15]假設所有的源活動是時間基函數(Temporal Basis Functions,TBFs)的線性組合，并通過變分貝葉斯(Variational Bayesian, VB)期望最大化算法實現模型推斷。文獻[16]將空間建模與源信號的時間建模相結合，在空域和時域中分別使用L1范數正則化和L2范數正則化，并利用時間基函數得到了空間上稀疏時間上連續的源估計。文獻[17]提出基于貝葉斯框架的時空基函數全數據驅動源成像，對彌散源有更精確的定位。這些研究表明時域信息有利于ESI的精準重構。

為同時利用fMRI的空間信息和EEG信號的時域信息，本文在研究[14,17]基礎上，提出了一個更廣義的基于fMRI腦網絡和時空約束的EEG源重構算法(Functional Network based Spatio-Temporal Constrains Source Imaging, FN-STCSI)。對于時域約束，FN-STCSI假設源信號是若干時間基函數的線性組合。以往利用TBFs的研究需要事先將TBFs確定下來。但事先設定的TBFs會明顯地影響源成像的質量，確定TBFs的模態和數量依然是一個挑戰。針對這一問題，本文在經驗貝葉斯概率框架下，根據矩陣分解的思想，將源信號S分解為若干個TBFs的線性組合：S=W Φ。通過數據自驅動的方式，時間基函數矩陣Φ和相應的權重矩陣W同時從EEG數據中學習得到。對于空間約束，在經驗貝葉斯概率空間下，FN-STCSI假設W的先驗協方差是多個空間協方差分量的加權和。首先，FNSTCSI利用獨立主成分分析提取fMRI中的腦功能網絡構建先驗協方差；其次，為體現腦活動空間連續局部同質的特性，FN-STCSI引入MSP和DDP作為空間協方差。利用變分貝葉斯推斷技術，通過最大化自由能來逼近模型參數的后驗概率分布，利用數據自驅動的方式獲得不同fMRI空間先驗和EEG空間先驗對源信號的貢獻，實現EEG-fMRI時空信息融合，以準確估計皮層源活動的位置、尺寸和時間序列活動。

本文主要貢獻如下：(1) 在經驗貝葉斯框架下，提出了一種具有時空約束的基于fMRI功能網絡信息的 ESI方法。(2) 基于fMRI 獨立主成分網絡和 EEG 記錄構建空間協方差，使用自相關決策(Automatic Relevance Determination, ARD)自動選擇與大腦活動相關的協方差分量。

本文的結構如下：第2節介紹貝葉斯時空正向模型、該算法的具體實現以及相關源成像算法，第3節介紹仿真設計和性能指標，第4節將提出的算法應用于模擬和真實EEG數據集進行分析，第5節對本文的研究成果進行討論，第6節對本文進行總結。

2 方法

2.1 FN-STCSI算法

根據正向建模，EEG信號與腦源信號之間的關系[3]可表示為

根據式(7)，FN-STCSI的全概率模型為

對應的概率圖如圖1所示，其中，方框變量為已知變量，圓圈變量為未知變量。變量間箭頭表示依賴關系。

2.2 先驗約束

許多神經影像數據分析的結果已表明皮層神經活動可表示為功能同質的皮質團塊。文獻[9]提出基于皮層表面的DDP方法將皮層表面劃分為若干互不重疊的團塊。DDP利用EEG數據來指導空間聚類，首先使用多變量源預定位對腦源進行預定位，將得到的腦源作為種子點，使用區域生長算法迭代種子點周圍的區域，直到達到指定空間大小。這樣整個大腦被分割成若干皮質團塊，最終得到基于EEG信息構建的空間協方差，本文將此方法的協方差記為C2。

為利用fMRI信號的空間信息，FN-STCSI利用空間獨立主成分分析將fMRI信號分解為若干獨立主成分，對fMRI獨立成分與EEG腦源空間進行配準，將fMRI體素的Z分數賦值給離該體素最近的EEG源空間格點，使得fMRI獨立成分轉變為EEG源空間的強度矩陣Wd×ks，d為EEG源格點數，ks為獨立成分個數。將每個fMRI獨立成分定義為腦功能網絡。將強度矩陣W二值化為激活矩陣U，若Wij的Z分數大于等于3.0，則Uij=1.0，否則Uij=0。基于腦功能網絡的協方差定義為[14]

本文同時考慮了以上3種協方差成分，協方差成分集C= {C1,C2,C3}。每個空間先驗的貢獻由經驗貝葉斯框架下的EEG信號自動確定，因此FNSTCSI可以融合來自EEG和fMRI信號的空間信息，提高腦源估計的性能。

本文使用臺式機(i7-8700U CPU 3.2 GHz和8 GB RAM)進行數值仿真實驗。EEG傳感器的數量為70，源的數量為8196，基于EEG信息的DDP協方差成分數量大約為 350。MSP的協方差成分數量為256×3=768，基于fMRI信息的協方差成分數量為3～22 ，如果設置初始TBFs的數量5，FN-STCSI在大約110次迭代后收斂，耗時大約20 s。

3 數值仿真實驗設計與性能評價指標

實驗采用開源的多模態人臉識別數據集(https://openneuro.org/datasets/ds000117/)，該數據集共16個被試參與了對熟悉、陌生和干擾面孔的面部識別任務。EEG數據利用Elekta VectorView系統(70個電極，鼻參考)采集。MRI數據通過3T Siemens TIM Trio采集，其中包括1個1 mm×1 mm×1 mm的T1加權結構磁共振成像以及若干3 mm×3 mm×4 mm的T2加權功能磁共振成像。EEG和fMRI數據利用SPM12(Statistical Parametric Mapping)進行預處理。數據集介紹以及具體處理過程參照文獻[20]。

3.1 EEG數據預處理

本文采用被試15的EEG數據，對BIDS(Brain Imaging Data Structure)格式的原始EEG數據進行格式轉換，采樣率保持為1100 Hz，進行6～40 Hz的帶通濾波，然后提取熟悉面孔識別任務下每個試次—200～600 ms (0 ms代表刺激開始)的EEG數據進行分析。利用眼電圖通道對每個試次進行眼電偽跡去除，最后將所有試次EEG信號的平均應用于仿真EEG數據構造和實驗數據分析。

3.2 正演模型的建立

采用被試15的個體頭模型，將被試的MRI與標準MNI(Montreal Neurological Institute)空間配齊，得到標準化的頭模型。將皮層表面分割為8196個網格，每個網格作為一個源信號，源信號方向垂直于皮質表面。利用SPM12，采用邊界元模型，基于70個EEG電極分布位置，得到導聯矩陣L ∈R70×8196。

3.3 fMRI功能網絡的提取

利用SPM12對fMRI數據進行預處理，包括切片間圖像采集時間校正，頭動校正，以及MNI空間配準。所有被試fMRI預處理完成后，采用GIFT(Group ICA of fMRI Toolbox)提取功能網絡。通過主成分分析，最終得到25個獨立成分。每個空間獨立成分進行Z變換，Z分數大于3的體素設定為激活體素。

3.4 模擬EEG數據的構造

3.5 性能評價指標與對比算法

為驗證FN-STCSI的有效性，本文將以下5個算法作為對比方法：(1)L1R-SSSI[10]；(2)wMNE[5]；(3)fwMNE[12]；(4)MSP[8]；(5)NESOI[14]。其中fwMNE和NESOI為EEG-fMRI皮層源成像算法，L1R-SSSI, wMNE和MSP為EEG源成像算法。同時，通過蒙特卡羅數值仿真，分別比較了在不同SNR, SNIR，不同有效fMRI先驗個數和無效fMRI先驗個數下，各個算法的性能(每種情況分別進行50次蒙特卡羅仿真)。

本文利用4個指標定量評估不同算法的成像性能：(1)受試者工作特征曲線下的面積(Area Under the receiver operating characteristic Curve, AUC)，用于評估重建源的檢測靈敏度和特異性；(2)空間彌散度(Spatial Dispersion, SD)，用于描述重建源的空間模糊度和彌散度；(3)定位誤差(Distance of Localization Error, DLE)，用于描述重構源相對于模擬源間的誤差距離。(4)相對均方誤差(Relative Mean Square Error, RMSE)，用于描述重構源與模擬源之間的相對均方誤差，表示重構波形的準確性。對于不同的源成像算法，其AUC越大，SD,DLE, RMSE越小，代表成像質量越好。

為評估結果統計顯著性，本文先采用ANOVA分析。如果ANOVA分析統計結果顯著，再利用Bonferroni校正的配對t檢驗對FN-STCSI和其他對比方法兩兩進行比較，以驗證FN-STCSI得到的結果是否顯著更優。由于L1R-SSSI是基于單個時間點的算法，L1R-SSSI選擇EEG幅值最大的時刻進行源重構，其余算法對0～600 ms整個時間段進行源重構。對于成像可視化，本文顯示指定時間段或單個時間點重構源的絕對值，成像閾值由Otsu方法[21]確定。

4 實驗結果

4.1 數值仿真實驗結果

4.1.1 有效fMRI先驗影響

為研究fMRI與EEG在不同耦合程度下的成像結果，本節測試在總fMRI先驗為3的情況下，不同有效fMRI先驗個數下的性能(有效fMRI先驗為來自生成EEG數據的D)，固定SNR=5 dB, SNIR=0 dB。圖3是fwMNE, NESOI和FN-STCSI在有效fMRI先驗個數從0到3時的性能評價指標。隨著有效fMRI先驗個數的增加，這3個方法的成像性能均能有效提升，表現為逐漸增大的AUC(p ＜0.05)，逐漸減小的SD(p ＜0.05)和DLE(p ＜0.05)值。當有效fMRI先驗個數為3時，此時fMRI先驗與EEG仿真源完全吻合，fwMNE能夠近乎完美地重構皮層源信號，但這種情況在實際中基本不可能發生。

4.1.2 無效fMRI先驗影響

考慮到實際情況，從fMRI提取的大部分功能網絡與EEG源無關。為驗證無關fMRI網絡對源成像性能的影響，本文保持有效fMRI先驗為2個，而無效fMRI先驗個數從5增加到20。圖4表明，fwMNE隨無效fMRI先驗的增加，性能顯著降低，表現為顯著降低的AUC(p ＜0.05)和逐漸升高的SD(p ＜0.05),DLE(p ＜0.05)，證明fwMNE高度依賴性和fMRI先驗的有效性。得益于ARD機制，FN-STCSI和NESOI根據EEG數據判斷每個fMRI先驗對源信號的貢獻，從而去除冗余的協方差成分。隨著無效fMRI先驗個數的增加，這兩個算法的性能指標沒有明顯變化，表明其對無效fMRI先驗具有較高的魯棒性。這在實際應用中非常重要，因為事先并不知道哪些功能網絡是有效的。

4.1.3 不同信噪比下源信號重構性能比較

考慮到EEG與fMRI不完全耦合這一普遍情況，并且事先并不知道哪些功能網絡為有效fMRI先驗，哪些是無效fMRI先驗，本文選取2個有效fMRI先驗，同時將剩余的無關功能網絡也納入先驗中(2個有fMRI先驗，22個無效fMRI先驗)，并固定SNIR=5 dB。

如圖5所示，測量噪聲對所有算法均有明顯影響。隨著信噪比的增加，所有方法的AUC(p ＜0.05)逐漸增加，SD(p ＜0.05)逐漸降低。在所有SNR條件下，相對于基于經驗貝葉斯框架的MSP, NESOI和FN-STCSI算法，基于L2范數約束的算法(wMNE,fwMNE)具有更大的SD(p ＜0.05)和DLE(p ＜0.05)，表明產生了更模糊的源信號估計。相對于僅利用EEG信息的MSP算法，同時利用EEG-fMRI信息的FN-STCSI和NESOI具有更好的成像性能。相對于僅使用了空間約束的MSP和NESOI，FN-STCSI具有最大的AUC(p ＜0.05)和最小的SD(p ＜0.05)。

圖6為不同SNIR條件下的算法性能。隨著腦源噪聲逐漸降低，MSP, NESOI和FN-STCSI性能有所提高，表現為逐漸增加的AUC(p ＜0.05)值，逐漸減小的SD(p ＜0.05)和DLE(p ＜0.05)值。而L1RSSSI, wMNE和fwMNE對腦源噪聲并不敏感，各性能指標無明顯變化。經不同SNR和SNIR條件下的對比，證明FN-STCSI對傳感器噪聲和腦源噪聲具有強魯棒性。

圖7是SNR=10 dB和SNIR=5 dB條件下成像的一個例子，其中，總fMRI先驗個數為25，其中3個為有效fMRI先驗，剩余22個為無效fMRI先驗。SNR為10 dB, SNIR為5 dB。L1R-SSSI在頂葉區出現了許多不存在的“偽源”，wMNE和fwMNE的成像結果則過于彌散和模糊，MSP和NESOI沒能準確重構源的尺寸。FN-STCSI得到了最高的AUC，最低的SD和DLE，說明FN-STCSI能準確地重構源的位置和尺寸。

4.2 真實數據檢驗

將FN-STCSI應用與真實的人臉識別數據集中，預處理過程同3.1節—3.3節。圖8(a)是被試15熟悉面部識別任務平均后產生的EEG波形。對于基于單個時間點進行源成像算法L1R-SSSI，選取刺激后177 ms的腦電信號進行源重建，其余算法將利用0到600 ms間EEG信號進行源重建，結果如圖8(b)所示，除L1R-SSSI外，wMNE, fwMNE,MSP, NESOI和FN-STCSI均定位到了雙側梭形皮質，這與先前的研究一致。但wMNE和fwMNE過于彌散。對于人臉識別數據的分析，本文得到了符合神經生理學研究的結果，也與文獻[14,19]的結果一致。

5 討論

由于ESI問題是高度欠定的，需要生理或物理先驗信息約束解空間以獲得唯一解。基于L2范數約束的MNE和fwMNE對空間范圍不敏感，往往產生過度模糊和彌散的估計[22]。如圖5和圖6的高SD和高DLE已經說明了這一點。為提高重構源的空間分辨率，一些研究提出基于稀疏約束的算法(如基于L1范數的方法)，但原始源空間中的稀疏懲罰會產生過度聚焦的估計，無法恢復源的空間范圍。為重構彌散源，文獻[9]提出的算法提出了協方差分量，這些協方差分量是基于DDP和局部空間平滑模型推導的，假設每個包內的函數是同質的，得到了比傳統的基于L1和L2范數方法更好的源重構。另外，基于EEG-fMRI時空分辨率互補的特性，文獻[14]利用fMRI的空間信息對EEG進行約束，以提高皮層神經活動重構性能。但這些彌散源成像方法未考慮源的潛在時間相關性，可以利用這些相關性來進一步增強源成像性能。

本文將時域空域信息結合，更精確地重建了腦源的位置、范圍和時間過程。本文使用ICA提取fMRI功能網絡，并將提取的功能網絡全部納入協方差先驗，但提取的功能網絡并非全部與大腦活動相關。在fMRI先驗有效的情況下，即使是fwMNE都有非常好的效果(如圖3所示)，但隨著無效fMRI先驗增多，成像性能也隨之下降(如圖4所示)。實際應用中，本文不知道某個fMRI功能網絡是否有效。事實上，只有少數功能網絡與大腦活動相關。對于fMRI先驗的有效性，本文采取ARD策略，根據EEG數據自動決策某fMRI先驗的有效性，從而融合有效fMRI先驗信息并過濾無效fMRI先驗。在信噪比對比實驗中，由于存在有效fMRI功能網絡，FN-STCSI與NESOI各方面指標均優于沒有利用fMRI功能網絡的L1R-SSSI、MNE和MSP算法，證明利用fMRI信息能有效提高源成像的性能。由于FN-STCSI利用了EEG的時域信息，在空間先驗相同的條件性能下優于NESOI，證明時域約束的加入能有效提高成像性能。

本文主要進行任務態的EEG/fMRI人臉識別分析，而利用靜息態fMRI得到的功能網絡，可使用FN-STCSI分析靜息態EEG/fMRI數據。另外FNSTCSI假設測量噪聲服從高斯分布，不能很好地擬合EEG測量過程中存在的頭動、眼動等異常干擾。為獲得更魯棒的ESI成像算法，下一步工作將利用Laplace先驗或廣義高斯先驗分布表示測量噪聲，以進一步克服EEG信號中的異常干擾。此外，鑒于深度神經網絡的深層結構和非線性具有的強大的特征表示能力，其也被逐漸應用于ESI問題求解[23,24]。文獻[24]提出了深度腦神經網絡(Deep-BraiNNet)實現了魯棒的稀疏皮層源信號時空活動重構。DeepBraiNNet利用長短期記憶的循環神經網絡表示皮層源信號的時空活動，數值仿真和在運動想象實驗數據的結果表明，其能有效克服測量噪聲的影響，準確估計皮層神經活動。未來也將進一步利用深度神經網絡模型表示皮層彌散源信號的時空關系，并融合fMRI信息，實現數據和模型驅動的EEG-fMRI源成像算法。

6 結束語

本文提出了一種新的貝葉斯算法FN-STCSI。通過將源分解為若干TBFs的線性組合，有效地利用的EEG中的時域信息；利用ICA從fMRI提取的功能網絡作為協方差，有效地利用了fMRI的空間信息，從而實現了高時間分辨率EEG與高空間分辨率fMRI的有效融合，有效地解決了彌散源成像的問題。仿真結果表明，FN-STCSI算法的性能優于L1R-SSSI, wMNE, fwMNE, MSP和NESOI等基準算法。FN-STCSI在探測靈敏度、源范圍、位置和幅度誤差方面獲得了更精確的估計。對于真實的實驗數據，FN-STCSI也提供了更有意義的神經生理學結果。未來工作將應用FN-STCSI到解碼腦機接口的EEG信號[25]以及分析各種精神疾病背后的皮層腦網絡[26,27]。