本发明属于神经信号处理研究领域，涉及脑功能网络建模方法，特别涉及一种静息态同步EEG-fMRI的脑功能网络建模方法。
背景技术：
人类所有的高级认知功能如思维、情感和意识依赖大脑，大脑是一个极其复杂的系统。大脑虽然只占人体2％的质量，却消耗了20％的能量。在没有任务的时候，大脑也会进行活动，人们将大脑这种不由外部刺激引起，不受被试意愿控制，自发产生神经活动的状态称作“静息态”。研究发现大脑在任务态比静息态增加的能量消耗通常不超过5％，而占大脑总能耗60％以上的能量用于自发神经活动，可见静息态在人脑功能中占据着重要地位。
脑功能网络是多个大脑区域的集合，这些脑区在静息态和任务态下都会进行同步活动。随着各种无创脑功能成像技术的出现，人类已经可以在微观和宏观不同尺度水平上探测大脑的活动。脑电图(Electroencephalogram，EEG)可以检测到脑活动在毫秒级的变化情况，具有很高的时间分辨率，而功能磁共振成像(functional Magnetic Resonance Imaging，fMRI)能以毫米级的精度定位同步振荡发生的位置，具有很高的空间分辨率。目前已经出现了多种基于EEG和fMRI的单模态脑功能网络建模方法。
静息态fMRI已被应用于多种神经疾病的研究中，如癫痫、阿尔兹海默症、精神分裂症、注意力不集中症、抑郁症等。对于阿尔兹海默症等疾病，fMRI研究结果与早期的研究相一致，但对于精神分裂症等疾病，研究结果的再现性较差。为了更深入的研究这些疾病，需要理解大脑中自发振荡的潜在机理，这可以通过多模态研究方法实现。此外，利用fMRI对静息态脑功能网络建模，由于没有外界刺激，要对所建网络进行解释非常困难，将EEG特征与fMRI信号结合，可以更好地解释脑功能网络代表的生理意义，并且对于还不清楚的EEG特征也能找到fMRI信号的支持证据。因此，越来越多的研究人员开始关注同步EEG-fMRI的研究。
同步EEG-fMRI的数据分析方法有比较方法、非对称融合方法和对称融合方法。比较方法仅仅对EEG和fMRI数据进行比较分析，没有进行模态间的数据融合；非对称融合方法假设EEG和fMRI测得的神经活动完全相同，这种假设缺乏一定的生理基础；目前已有的对称融合方法大多只关注事件相关电位与fMRI数据的融合，尚未见到对静息态同步EEG-fMRI的对称融合方法的研究。
基于以上研究背景，本发明提出了一种静息态同步EEG-fMRI的脑功能网络建模方法，该方法对同步数据进行对称融合，并利用融合结果进行脑功能网络建模，为有效整合两种模态的信息提供了一种新的思路和方案。
技术实现要素：
本发明提出了一种静息态同步EEG-fMRI的脑功能网络建模，该方法能够有效融合静息态下两种模态的信息，包含以下具体步骤：
S1.EEG信号预处理，提取带限能量(Band Limited Power，BLP)信号并构建回归项。同步采集EEG、fMRI数据时，EEG的采样率远远大于fMRI，大部分的EEG信号找不到时间上对应的fMRI信号，通常将EEG信号的BLP变化与一个表示血氧水平依赖(Blood Oxygenation Level Dependent，BOLD)在刺激呈现后变化的核函数进行卷积，用得到的信号作为拟合fMRI数据的回归项，通过比较回归项与BOLD信号的关系来研究EEG与fMRI之间的关系。本发明首先对同步记录的EEG信号通过预处理得到比较纯净的EEG信号，随后用独立成分分析(Independent Component Anaylysis，ICA)将信号分解为相互独立的源信号，最后计算源的BLP信号并构建回归项。
S2.fMRI信号预处理，提取各脑区的BOLD信号。原始fMRI数据中含有多种噪声，这些噪声可能是由于机器本身产生，也可能是被试的运动或生理活动产生，它们会对分析造成很大的干扰，需要通过预处理来提高原始数据的信噪比。由于fMRI数据中体素数目远远大于采集点数，要对fMRI数据进行时域分析，首先需要对其进行空间降维。本发明利用标准解剖模板将大脑灰质分成不同的脑区，将每个脑区中所有体素时间序列的平均作为该脑区产生的信号。
S3.对fMRI提取的各脑区BOLD信号和由EEG得到的回归项分别进行主成分分析(Principal Component Analysis，PCA)。在fMRI分析中，由于分割后的相邻脑区会有同步活动，使得相邻脑区计算得到的平均时间序列会有一定的相关性；由于大脑的网络同步活动，相隔较远脑区的时间序列之间也会存在较高的相关性。在EEG分析中，用每个独立成分的BLP信号构造的回归项之间也存在相关性。信号之间的这种相关会给之后的统计分析带来许多障碍。基于以上原因，首先要对fMRI中提取的各脑区BOLD信号和由EEG得到的回归项分别进行PCA，消除原信号之间的相关性和信息重叠。另外，静息态下大脑消耗的能量大部分用于自发神经活动，因此可以认为自发神经活动产生的信号是整个大脑信号的主成分，通过主成分分析的方法将它们提取出来，同时也去除了一些噪声成分。
S4.对fMRI提取的各脑区BOLD信号和由EEG得到的回归项的主成分进行典型相关分析(Canonical Correlation Analysis,CCA)。进行该步骤基于以下三点原因：(1)一种EEG模式可能与多个脑区的BOLD信号相关，一个脑区也可能产生不同模式的EEG信号，这表明，EEG和fMRI之间不是简单的一对一的关系，需要寻找一种方法，研究多个EEG模式与多个BOLD信号之间在总体上的关系。(2)同步颅内EEG-fMRI的比较方法研究表明，局部场电势的带限能量变化与局部BOLD信号是相关的。这表明，虽然EEG和fMRI两种成像方式基于不同的成像原理，但它们记录的自发神经活动信息有一定程度的重叠，因此可以通过比较两个模态信号之间的相关关系，判断它们是不是同一源产生的信号。(3)静息态大脑同时存在多个功能网络，它们分别代表不同的生理活动，且这些功能网络的激活时间过程互不相关。
S5.静息态同步EEG-fMRI的脑功能网络建模。将典型相关变量与EEG构建的回归项求相关，相关系数越大，表示该EEG模式与该大脑活动模式越相关，就越可能是由该脑活动产生的EEG信号，相关系数较大的EEG模式一起组成该脑活动所产生的EEG信号。用典型相关变量与脑区平均时间序列求相关，相关系数越大，表示该脑区与该大脑活动模式越相关，越有可能是参与该脑活动的脑区，相关系数较大的脑区一起构成该脑活动的功能网络。
附图说明
图1为本发明所述的建模方法流程图；
图2为利用EEG构建回归项的流程图；
图3为EEG回归项和BOLD信号分别进行主成分分析的示意图；
图4为BOLD信号(a)和EEG回归项(b)的主成分能量和累计贡献率；
图5为对EEG回归项和BOLD信号的主成分进行典型相关分析的示意图；
图6为由EEG回归项和BOLD信号的主成分求得的典型相关系数的直方图；
图7为脑功能网络建模结果。
具体实施方式
下面结合实例及附图对本发明作进一步详细的描述,但本发明的实施方式不限于此。
如图1所示为一种静息态同步EEG-fMRI的脑功能网络建模方法的具体流程图，包含以下部分：
一、EEG信号预处理，提取带限能量信号并构建回归项
1、EEG信号预处理
预处理包含以下六个步骤：(1)使用自适应伪迹相减法去除梯度伪迹和心电伪迹；(2)将EEG数据降采样到200Hz；(3)为了移除线性漂移与低频干扰，对降采样后的数据进行1Hz的高通滤波；(4)通过陷波滤波去除50Hz的工频噪声；(5)将所有电极信号重新参考到全导的平均信号；(6)去除坏数据段。通过以上处理得到比较纯净的EEG信号。
2、EEG信号独立成分分析
(1)采用信息最大化ICA算法将61导EEG信号分解为61个源信号；(2)经过预处理的EEG中仍然会存在一些残留噪声，使用EEGLAB扩展包ADJUST对独立成分进行分析，识别并去除噪声成分。
3、计算带限能量信号
通过短时傅里叶变换计算EEG的带限能量信号。对每个独立成分使用一个无重叠窗计算短时傅里叶变换，为了与fMRI每采集一个图像所用的时间匹配，窗的大小选择1.94s。之后将EEG分为delta(1～4Hz)，theta(4～8Hz)，alpha(8～13Hz)，beta(13～30Hz)，gamma(30～50Hz)五个频带，对每个频段内所有频率点在同一时刻的能量值求平均，得到五个频段的BLP变化，如图2所示。
4、构建回归项
通常神经脉冲电活动和与之相应的BOLD变化有4～8秒的延迟，将上一步得到的BLP信号作为一个自变量，通过与标准血氧响应函数(Hemodynamic Response Function，HRF)卷积，构建回归项，如图2所示。
二、fMRI信号预处理，提取各脑区的BOLD信号
1、fMRI信号预处理
在对fMRI数据进行分析之前，首先要进行一些必要的预处理，具体步骤如下：
(1)格式转换。将采集到的DICOM格式的文件转换为NIFTI格式，得到666个以img和hdr为后缀的文件，每个img后缀文件是由32张不同层的大脑图片组成。
(2)去除最开始采集的数据。fMRI设备刚启动时采集的数据不稳定，所以去除前十个时间点采集的图像。
(3)时间层校正。每扫描一个全脑图像需要1.94秒，该图像由32层不同时刻记录的图像组成，需要进行时间层校正，使得各层图像都是在相同时刻记录的。本实验中，每个图像层数为32，TR为1.94s，扫描顺序为降序，参考层为第16层。
(4)头动校正。本实验fMRI扫描需要连续进行大约半小时，被试不可避免的会有头部运动。此外，被试的自发生理活动如呼吸和心跳也可能导致头部运动，所以需要进行头动校正，减少头动对数据的影响。设定当头在某个方向的平移大于0.5毫米或旋转角度超过1度时，认为被试在这一时刻头动较大，去除这一时刻数据。
(5)空间标准化。由于不同被试脑部结构和大小存在一定的差异，为了方便后续的分析，需要将每次扫描的大脑数据标准化。空间标准化是对大脑整体进行形变和调整，即从原始空间中估计到标准空间上，调整后统一标准化到标准脑上。本实验中，使用EPI模板进行空间标准化。
(6)去线性漂移。由于扫描过程中机器温度升高或者被试的适应，随着时间变化会存在一个线性趋势，因此要去除这种线性漂移。
(7)滤波。人体的一些自发生理活动，如心跳、呼吸等会产生噪声，这些噪声的频率通常比大脑自发神经活动要高，本实验选择0.01～0.08Hz的带通滤波器将它们去除。
2、提取脑区平均时间序列
由于fMRI原始数据量特别大，经过预处理后，一个大脑的体素有61*73*61个，达到了百万级，需要采用合适的方法对原始数据进行空间降维。一种常用的方法就是将大脑灰质划分成不同的脑区，提取各个脑区的平均时间序列。本发明使用解剖自动标记(Anatomical Automatic Labeling，AAL)图谱把大脑灰质分为90个区域，之后对每个脑区中所有体素的时间序列求平均，得到该脑区的平均时间序列。
三、对fMRI提取的各脑区BOLD信号和由EEG得到的回归项分别进行主成分分析
经过前面的预处理步骤，分别得到了BOLD信号和EEG回归项的数据矩阵：X1和X2。其中X1为624*90的矩阵，行代表时间，列代表脑区；X2为624*225的矩阵，行代表时间，列由55个独立成分的5个频带串联得到的。
对X1和X2分别进行PCA，如图3所示。计算X1和X2各自的相关系数矩阵、特征值和特征向量，并计算每个特征值的累计贡献率。如图4(a)所示，fMRI前15个主成分的能量累计贡献率达到了87％，这15个主成分构成一个624*15的矩阵P1；如图4(b)所示，EEG前20个主成分的能量累计贡献率达到了70％，这20个主成分构成一个624*20的矩阵P2。此时，P1和P2分别表示从fMRI和EEG中提取的静息态大脑活动的主要特征。
四、对fMRI提取的各脑区BOLD信号和由EEG得到的回归项的主成分进行典型相关分析
对得到的fMRI的主成分P1和EEG的主成分P2进行CCA，如图5所示，得到15对典型相关变量S1和S2，它们都为624*15的矩阵，两个矩阵的对应列为一对典型相关变量。每对典型相关变量的典型相关系数如图6所示，设定阈值为0.5，大于该阈值的认为是典型相关的，得到七对最相关的典型相关变量。S1的每个典型相关变量表示一个脑活动对应的fMRI时间序列，S2的每个典型相关变量表示一个脑活动对应的EEG时间序列，且S1和S2的对应列具有较高的相关系数，认为它们是由同一个脑活动产生的信号。
五、静息态同步EEG-fMRI的脑功能网络建模
将S1的每一列与X1的每一列求相关，典型相关变量与各脑区BOLD信号的相关系数直方图如图7(a)所示，系数越大表示脑区与大脑活动越相关，越有可能是参与脑活动的脑区，相关系数较大的脑区一起构成该脑功能网络如图7(c)所示；将S2中的每一列与X2的每一列求相关，典型相关变量与各EEG回归项的相关系数直方图如图7(b)所示，系数越大表示该EEG模式与该大脑活动越相关，越可能是由该脑活动产生的EEG信号，相关系数较大的EEG模式一起组成该脑活动所产生的EEG信号，对图7(b)中每个频带内的相关系数求均值并取绝对值，得到典型相关变量与EEG各频带的相关度，如图7(d)所示。
利用本发明提出的建模方法得到了七个主要的脑功能网络，并且得到各个网络对应的EEG信号与脑区，这些脑网络与现有的研究结果具有很高的一致性，证明了该脑功能网络建模方法的有效性，为静息态脑功能研究提供了一种新的途径。