CN102496159A - 一种联合fmri和meg的脑区因果连接检测方法 - Google Patents

Abstract

本发明涉及一种联合功能磁共振成像(FMRI)和脑磁图描记术(MEG)的脑区因果连接检测方法，该方法包括以下步骤：首先，对经过数据预处理的FMRI图像进行激活区坐标的提取；然后，基于提取的FMRI激活区坐标，在经过预处理的MEG数据上提取对应位置区域的脑区时间序列；最后，依据提取的MEG脑区时间序列，计算脑区间的因果连接强度和方向，利用有向网络图显示显著的连接。本发明所述方法是一种有效的联合FMRI和MEG两种成像模态的脑区因果连接检测方法，能够比仅利用FMRI图像的检测方法更完整准确地检测脑区因果连接。

Description

一种联合FMRI和MEG的脑区因果连接检测方法

技术领域

本发明属于图像处理领域，具体涉及一种利用脑功能成像技术的脑区因果连接检测方法，尤其涉及一种联合功能磁共振成像(FMRI)和脑磁图描记术(magneto encephalography，MEG)进行脑区因果连接强度和方向的检测的方法。

背景技术

随着科学技术的进步，无创性脑功能成像技术日新月异，为我们用多种成像技术进行脑的认知和研究提供了广阔的舞台。单一成像模态，由于成像机理决定了每一种成像技术都具有相应的优点和不足，因此没有一种单一的适合所有实验研究或临床应用的成像方法。目前，功能磁共振成像(Functional Magnetic Resonance Imaging，FMRI)技术以其高空间分辨率、非侵入式等特点在神经疾病诊断治疗方面得到了广泛应用。FMRI技术是通过测量由神经活动引起的脑血流和脑血氧等成分变化而造成的磁共振信号变化来反应人脑内部活动的。由于在代谢水平上测量磁共振信号的变化相较神经元活动具有一定的时间延迟，因此，FMRI技术仅能达到秒级时间分辨率精度。

近年来，脑磁图描记术(Magneto encephalography，MEG)，又称脑磁图仪，逐渐在神经影像领域得到应用。MEG技术可以实时测量脑内神经元之间信号传递时突触后电位变化所产生的磁场信号，因此能够准确测量毫秒级神经元电活动，具有非常强的时间定位能力。但是由于该技术难以准确定位颅内信号源，其空间分辨率较低。

因此，如果将FMRI技术与MEG技术的优势结合起来，我们就能达到同时提高空间分辨率和时间分辨率的目的。

脑功能一般遵循两个基本组织原则：功能集成化和功能特异化。在大空间刻度上，一个复杂的脑功能可能会由许多功能特异的脑区通过相互作用(集成)来完成；同时，某一特异性脑区也会对许多不同的刺激任务进行表示或加工，通过精细空间刻度上不同的分布式脑活动来对外部不同刺激进行表示。人脑是一个复杂的系统，在不同状态时期各个脑区之间的相互作用会发生变化，反映在脑区间的功能连接强度和方向上。格兰杰因果分析(Granger causality analysis，GCA)可以用于研究区域间的因果连接关系，为研究揭示人脑信息加工机制提供支持。

发明内容

本发明的目的在于提供一种同时具备高空间分辨率以及高时间分辨率优势的脑区因果连接强度及方向的检测方法。

为实现上述目的，本发明提出一种联合FMRI和MEG的脑区因果连接检测方法，包括以下步骤：

步骤Sa，对扫描脑区获取的FMRI图像进行预处理，并从预处理后的FMRI图像中提取激活区坐标；

步骤Sb，对扫描脑区获取的MEG数据进行预处理，并基于提取的FMRI激活区坐标，提取与该激活区坐标对应的脑区的脑区时间序列；

步骤Sc，依据提取的MEG脑区时间序列，计算脑区间的因果连接强度和方向，利用有向网络图显示显著的连接。

本发明所述的联合FMRI和MEG的脑区因果连接检测方法，能够有效地利用FMRI高空间分辨率和MEG高时间分辨率的优势。真实试验数据结果表明，本发明比传统的仅利用FMRI数据更完整准确地检测脑区因果连接，为脑功能数据分析和脑区因果连接检测方法提供了一种新思路，在基于多模态医学成像技术的脑认知加工研究中具有一定的应用价值。

附图说明

图1为本发明提供的联合FMRI和MEG的脑区因果连接检测的方法流程图。

图2是实例中使用本方法检测的脑区因果连接强度图。

图3是实例中仅使用单一模态检测的脑区因果连接强度图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。

本发明主要考虑利用功能磁共振成像(FMRI)技术高空间分辨率的优势，定位大脑FMRI图像上的关键神经元位置信息，继而在MEG技术获得上述神经元的高分辨率时间序列，从而在尽可能不丢失大脑信号活动时空二维信息的前提下，利用格兰杰因果分析方法进行脑区因果连接强度及方向的精确检测。

图1是本发明所涉及的一种脑区因果连接强度及方向的检测方法的一个具体实施例，其具体涉及联合利用功能磁共振成像(FMRI)和脑磁图描记术(MEG)的脑区因果连接强度及方向的检测方法。该实施例的具体实施步骤如下：

步骤Sa，对扫描脑区获取的FMRI图像进行预处理，并从预处理后的FMRI图像中提取激活区坐标；

1.对FMRI图像进行预处理

由于在磁共振成像的扫描过程中，被试的头动、图像内层间扫描时间的差别以及外加磁场的不均匀性等各种因素均会影响采集到的FMRI图像的质量，因此我们需要在保留脑功能图像细节的同时，通过使用脑功能磁共振图像与标准模板进行仿射配准变换方式的预处理，提高脑功能图像的信噪比。所谓仿射配准变换方式是指进行对原始图像缩放，旋转，平移后，根据原始图像和标准模板间的共有特征进行几何配准。

所述预处理步骤包括：切片扫描时间对齐、图像序列对齐、配准、标准化、平滑、滤波等。在本实施例中，对FMRI图像依次进行图像序列对齐、配准、标准化、空间平滑滤波。所属技术领域的技术人员也可以根据实际需要采用其他预处理方式，或者上述预处理方式的任意组合。

2.提取FMRI图像的激活区坐标

一般来说，提取激活区坐标的方法有基于一般线性模型(GeneralLinear Model，GLM)的方法，模式分类的方法和基于独立成分分析(Independent Component Analysis，ICA)的方法。其中基于ICA的脑功能激活成分提取模型如下：设A∈R^N×M为信号混合矩阵；X为脑功能信号的N维观测向量，C为由M(N≥M)维统计独立的信号分量构成的信号向量，其中包括与激活相关的脑功能信号分量。ICA分析法可以在混合矩阵A未知的情况下，寻找线性映射w，从观测的脑功能信号中提取不能被直接观测的原始信号c(≈y)：y＝w^Tx＝wAc。

本发明采用基于数据驱动的ICA方法，首先使用基于最小描述长度(Minimum Description Length，MDL)准则估计独立成分数，继而采用基于互信息最大化的算法计算解混矩阵，分离源信号，获得包含所需激活脑区的独立成分及其Talairach坐标(Talairach坐标是国际公认的一个标准化的脑图谱坐标系)。

步骤Sb，对扫描脑区获取的MEG数据进行预处理，并且基于步骤Sa提取的FMRI图像的激活区坐标，提取与该激活区坐标对应的脑区的时间维度信息(即脑区时间序列)；

1.MEG数据进行预处理

由于MEG数据在采集过程中容易受到环境噪声的影响，我们有必要对其进行预处理，以提高数据的信噪比；此外，由于MEG技术的高时间分辨率特性，采集的数据包含了大量的时间维度信息。因此，为保证计算精度与运算时间的平衡，预处理步骤包括：数据带通滤波，合适的重复下采样频率等。在本实施例中，对MEG数据依次进行0.5Hz-48Hz的带通滤波和300Hz的下采样频率。所属技术领域的技术人员也可以根据实际需要采用其他预处理方式，或者上述预处理方式的任意组合。

2.提取与激活区坐标对应的脑区的时间维度信息

如何从获得的MEG数据得到某个局部脑区的时间序列信号，就其本质来说，其实是一个脑磁逆问题，即利用扫描的磁场数据来推断脑内电流源的分布。已有研究证明了导体外的磁场数据无法唯一确定导体内的电流分布，因此脑磁逆问题是不适定的，我们必须在满足条件的解集中通过施加一定限制条件找出合理的解。目前，脑磁逆问题的求解方法主要分为两大类：一类基于图像重建技术，其基本思想是将电流分布区域离散成网格，对网络节点处磁源进行重建，重建的磁源参数主要是电流的强度。由于磁场方程中位置参数是非线性的，电流强度参数是线性的，重建时人为地限定了源的位置，因此磁场方程退化为线性方程。另一类为偶极子定位法，该类方法不限定磁源的位置，通常的求解过程是先建立一个目标函数，然后通过调整源的参数使目标函数达到极小，即通过执行一个非线性优化过程来获得磁源参数的解。

本发明的实施例采用基于图像重建技术的最小范数解(MinimumNorm Solution，MNS)求解方法。假设磁场方程组为D＝LS，其中，D为测得的磁场信号值，S代表节点处电流，L为已知的增益矩阵，取决于网格节点及磁场探测点的位置与头的电导率分布。MNS方法的限制条件是找到满足方程组解中范数最小的S^*，即

S^*＝min{||S||}，其中，||·||一般取Frobenius 2范数，即||S||＝(S^TS)^1/2，其中^T表示矩阵的转置。具体求解S^*时可以采用基于Moore-Penrose广义逆L⁺的方法，即

S^*＝L⁺D

对于广义逆L⁺的求解，本发明采用基于奇异值分解(Singular ValueDecomposition，SVD)的方法。设L为m×n阶实矩阵(m＜n)，其中，m为MEG扫描通道数，n为网格节点数。则对L进行奇异值分解可得

L＝U∑V^T

其中U＝(u₁，u₂，…，u_m)为m×m阶正交矩阵，V＝(v₁，v₂，…，v_n)为n×n阶正交矩阵，∑为m×n阶对角矩阵，即∑＝diag(δ₁，δ₂，…，δ_P)，其中P＝m，且δ₁≥δ₂≥…≥δ_P≥0，若矩阵L的秩rank(L)＝r，则上述序列中只有前r个值大于零，即δ₁≥δ₂≥…≥δ_r＞δ_r+1＝…＝δ_P＝0

那么奇异值分解可表述为：

$L={\mathrm{U\ΣV}}^T=\underset{i=1}{\overset{r}{\mathrm{\Σ}}}{\mathrm{\δ}}_i{u}_i{v_i}^T,$ 其中，u_i(i＝1，2，…，m)称为奇异值分解的左矢量，v_i(i＝1，2，…，n)称为奇异值分解的右矢量。

由此可求得Moore-Penrose广义逆L⁺为：

$L^{+}=\underset{i=1}{\overset{r}{\mathrm{\Σ}}}\frac{1}{{\mathrm{\δ}}_i}{v}_i{u}_i^T$

由此可根据S^*＝L⁺D求出MEG数据所隐含的源空间信息。从而根据ICA分析方法获得的局部激活值最大点体素的坐标，获得该局部脑区的时间序列信息。

步骤Sc，依据提取的MEG脑区时间序列，检测脑区间的因果连接强度和方向，利用有向网络图显示显著的连接。

目前，检测脑功能因果连接强度和方向的方法主要有两大类：采用动态因果模型(Dynamic Causal Modeling，DCM)和格兰杰因果分析(Granger Causality Analysis，GCA)。两者不同的是，DCM需要预先选定相互作用的区域，并假设这些区域的任意两个之间存在影响，这种预先假设的模型在验证一些有关大脑系统之间的假设时是起作用的，但如果对模型的指定出现偏差，就会导致错误的结论。而GCA则无需事先假设两点之间存在解剖结构的连接性，能更加完整和鲁棒地反应大脑之间的相互连接性。

本发明的实施例利用格兰杰因果分析方法检测研究脑区间的连接强度和方向。按照格兰杰因果的定义，若欲判断X是否引起Y，则考察Y的当前值在多大程度上可以由Y的过去值解释，然后考察加入X的滞后值是否能改善解释程度，如果X的滞后值有助于改善对Y的解释程度，则认为X是Y的格兰杰原因，Y是X的格兰杰结果。格兰杰因果关系检验假设了关于X和Y每一变量的预测的信息全部包含在这些变量的时间序列中。其步骤如下：

1.将当前的y对其所有的滞后项(y_t-1，y_t-2，…，y_t-q)做有约束回归，即

$y_t=\underset{i=1}{\overset{q}{\mathrm{\Σ}}}{\mathrm{\α}}_i{y}_{t-i}+u_t$

其中，y为一个激活区的时间序列，u_t为零均值白噪声，q为模型阶数，通过Akaike信息准则(Akaike Information Criterion，AIC)计算得到，

AIC(q)＝2log(det(∑))+2qm²/N

$\mathrm{\Σ}=R\left(0\right)+\underset{i=1}{\overset{q}{\mathrm{\Σ}}}A\left(i\right)R\left(i\right)$

其中X(t)＝(x₁(t)，x₂(t)，…，x_N(t))^T，x_i(t)(i＝1，…，N)为第i个激活区的时间序列，R(n)＝X(t)X^T(t+n)是激活区时间序列矩阵X(t)步长为n的协方差矩阵，∑为噪声协方差阵，m为激活区数目，N为时间点数目。

由此可得有约束残差平方和：

${\mathrm{RSS}}_R=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{{\hat{u}}_i}^2$

2.在上一步骤的回归式中加上滞后项x，做无约束回归，即

$y_t=\underset{i=1}{\overset{q}{\mathrm{\Σ}}}{\mathrm{\α}}_i{y}_{t-i}+\underset{j=1}{\overset{q}{\mathrm{\Σ}}}{\mathrm{\β}}_j{x}_{t-j}+e_t,$ 其中x，y为两个激活区的时间序列，e_t为零均值白噪声。

由此可得无约束残差平方和：

${\mathrm{RSS}}_{\mathrm{UR}}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{{\hat{e}}_i}^2$

3.构造F统计量如下

$F=\frac{({\mathrm{RSS}}_R-{\mathrm{RSS}}_{\mathrm{UR}})/q}{{\mathrm{RSS}}_{\mathrm{UR}}/N-2q-1}$

根据该统计量，我们即可在选定的显著性水平α下进行假设检验。若统计量F值大于自由度为q和N-2q-1的F分布对应的F_α值，则拒绝原假设(x_t不是y_t+1的格兰杰原因)，接受备择假设，即x_t是y_t+1的格兰杰原因。最后利用有向连接图显示显著的连接。

本发明所述的联合FMRI和MEG的因果连接强度的检测方法，可通过分析真实采集的实验数据得以说明：

(1)实验数据采集及分析

在实验中采用真实数据集作测试，一位健康男性被试参与了无刺激任务的静息状态FMRI扫描。

FMRI扫描采用T2^*加权梯度回波平面成像序列获取血氧水平依赖(Blood Oxygen Level Dependent，BOLD)数据。针对FMRI图像的预处理，我们采用统计参数图软件包(SPM5，http://www.fil.ion.ucl.ac.uk/spm/)对图像进行处理。之后利用GIFT软件(http://icatb.sourceforge.net/)进行激活区提取，分离出属于默认网络的成分并获得局部激活极大值点的Talairach坐标：左脑后扣带回(Left Posterior Cingulate Cortex，LPCC，-2，-40，24)，右脑后扣带回(Right Posterior Cingulate Cortex，RPCC，2，-44，21)，左脑内侧前额叶皮质(Left Medial Posterior Frontal Cortex，LMPFC，-4，42，-11)，右脑内侧前额叶皮质(Right Medial Posterior FrontalCortex，RMPFC，4，46，-11)。

脑磁图数据采集的采样频率设置为600Hz，利用加拿大CTF系统公司的CTF-151系统自带的综合三阶梯度的噪声消除技术对MEG信号进行降噪处理，对采集到的数据进行下采样(300Hz)和带通滤波(0.5-48Hz)。在此基础上进行脑磁图逆问题求解，利用前述4个极值点坐标提取MEG数据对应位置的时间序列。随后用本发明所述方法(方法A)结果对比仅用FMRI单一模态数据(方法B)的因果连接结果，比较相同激活区在两种不同模态下因果连接强度和方向的有向网络图。

(2)数据分析结果

利用上述两种方法(A和B)在实验数据上进行分析，其结果分别如表1和表2、图2和图3所示。其中，图2和图3所示箭头表示因果连接方向，对应表格中每行所对应的激活区指向每列对应的激活区，连接强度用线的宽度表示。通过分析结果，我们可知，用方法A检测的脑区因果连接显著强于方法B的结果。仅以LMPFC和LPCC两个激活脑区之间的关系为例，方法A检测到的由LMPFC指向LPCC的连接强度显著高于方法B的结果；同时，方法A检测到了由LPCC指向LMPFC的微弱连接，而方法B则检测不到该连接方向。

以上实验结果说明，本发明所述的联合FMRI和MEG的脑区因果连接检测方法，可以有效地提高脑区因果连接的检测性能，更完整、准确地反应大脑激活区的因果连接。

表1本发明方法(A)的激活脑区之间的因果连接

表2仅利用FMRI单一模态(B)的激活脑区之间的因果连接

以上所述，仅为本发明中的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉该技术的人在本发明所揭露的技术范围内，可理解想到的变换或替换，都应涵盖在本发明的包含范围之内，因此，本发明的保护范围应该以权利要求书的保护范围为准。

Claims (7)

1.一种联合FMRI和MEG的脑区因果连接检测方法，其特征在于，包括以下步骤：

步骤Sa，对扫描脑区获取的FMRI图像进行预处理，并从预处理后的FMRI图像中提取激活区坐标；

步骤Sb，对扫描脑区获取的MEG数据进行预处理，并基于提取的FMRI激活区坐标，提取与该激活区坐标对应的脑区的脑区时间序列；

步骤Sc，依据提取的MEG脑区时间序列，计算脑区间的因果连接强度和方向，利用有向网络图显示显著的连接。

2.如权利要求1所述的方法，其特征在于，步骤Sa中所述对扫描脑区获取的FMRI图像进行预处理，至少包括：

对扫描脑区获取的FMRI图像进行图像序列对齐、配准、标准化和空间平滑滤波。

3.如权利要求1所述的方法，其特征在于，步骤Sa中所述从预处理后的FMRI图像中提取激活区坐标，是采用基于数据驱动的ICA方法实现的。

4.如权利要求3所述的方法，其特征在于，所述采用基于数据驱动的ICA方法从预处理后的FMRI图像中提取激活区坐标，包括：

首先使用基于最小描述长度准则估计独立成分数，继而采用基于互信息最大化的算法计算解混矩阵，分离源信号，获得包含所需激活脑区的独立成分及其Talairach坐标。

5.如权利要求1所述的方法，其特征在于，步骤Sb中所述对扫描脑区获取的MEG数据进行预处理，至少包括：对扫描脑区获取的MEG数据进行数据带通滤波和下采样频率。

6.如权利要求1所述的方法，其特征在于，步骤Sb中所述基于提取的FMRI激活区坐标，提取与该激活区坐标对应的脑区的脑区时间序列，包括：

采用基于图像重建技术的最小范数解方法，对经过预处理的MEG数据进行脑磁逆问题求解，获得对应FMRI图像激活区坐标的脑区时间序列。

7.如权利要求1所述的方法，其特征在于，步骤Sc中所述依据提取的MEG脑区时间序列计算脑区间的因果连接强度和方向，是利用从MEG数据中提取的脑区时间序列，使用格兰杰因果方法分析时间序列之间的因果连接强度及方向。

Images