基于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]。