CN102366323A - 一种基于pca和gca的磁共振脑成像因果连接强度的检测方法 - Google Patents

Abstract

本发明公开了一种基于主成分分析PCA和格兰杰因果分析GCA的磁共振脑成像因果连接强度的检测方法，该方法包括以下步骤：对经过数据预处理的脑功能图像进行激活区域多体素时间序列的提取，得到激活区域多个体素的时间序列矩阵；使用PCA对每个激活区多个时间序列矩阵空间降维得到主要成分，平均主要成分的激活值得到一时间序列；建立所有激活区域的时间序列间的多变量自回归模型；计算各个时间序列间的偏相关系数；通过直接传递函数DTF方法计算dDTF值，得到激活脑区间的因果连接强度和方向；使用替代数据法统计检验连接强度的显著性，在有向网络图上将结果显示。在真实数据集上的实验说明，本发明所述方法是一种有效的磁共振脑成像因果连接强度的检测方法。

Description

一种基于PCA和GCA的磁共振脑成像因果连接强度的检测方法

技术领域

本发明属于图像处理领域，具体涉及一种磁共振脑成像因果连接强度的检测方法。尤其涉及使用主成分分析PCA和格兰杰因果分析GCA进行磁共振脑成像因果连接强度的检测。

背景技术

功能磁共振成像(functional Magnetic Resonance Imaging，fMRI)以其高时空分辨率，非侵入式等特点在神经疾病诊断治疗方面得到了广泛应用。fMRI一般指基于血氧水平依赖(blood oxygen level-dependent，BOLD)的磁共振成像，它通过测量由神经活动引起的脑血流和脑血氧等成分变化而造成的磁共振信号变化来反应脑活动。人脑是一个复杂的系统，在受到刺激条件或经历病变时各个脑区间存在的相互作用也在发生变化，反映在脑区间的功能连接强度和方向变化上。脑功能数据连接强度分析一般是根据测得的fMRI数据寻找其活动能显著区分不同试验条件(如，刺激条件和基线条件)的脑区，再检测脑区之间的相互连接强度。合理的脑功能数据连接强度分析方法需同时考虑到人脑功能一般组织原则以及fMRI数据本身一些特性。本发明为脑区间功能连接的分析和可视化提供一种有效途径，为神经疾病的早期诊断和临床治疗提供帮助。

脑功能一般遵循两个基本组织原则：功能集成化和功能特异化。在大空间刻度(规模)上，一个复杂的脑功能可能会由许多功能特异的脑区通过相互作用(集成)来完成；反过来，一个特异性脑区也会对许多不同的刺激任务进行表示或加工，通过精细空间刻度上不同的分布式脑活动来对外部不同刺激进行表示。另一方面，随着MRI技术的进步，fMRI图像空间分辨率一直在提高。传统fMRI数据在体素各维宽度约为4毫米时可获得较高的信噪比。当前，各维宽度为2毫米的体素已经可以在标准3T临床磁共振(MRI)机器上可靠地获得。随着超高场(≥7T)磁共振技术的使用，fMRI图像的空间分辨率正在向亚毫米的刻度上逼近。fMRI图像空间分辨率的提高，为我们在更精细的空间刻度上研究脑功能提供了可能。目前，传统的检测脑功能因果连接强度方法大多采用动态因果模型(dynamic causalmodeling，DCM)，它需要预先选定相互作用的区域，并假设这些区域的任意两个之间存在影响，这种预先假设的模型在验证一些有关大脑系统之间的假设时是起作用的，但如果出现了对模型的错误指定，如漏掉了一个间接作用或新起作用的区域，就会导致错误的结论。格兰杰因果分析(Granger causality analysis，GCA)检测方法可以用于研究区域间的因果关系，而不需要先验知识，克服了DCM的局限性，且其突出的优越性在于无需事先假设两点之间存在解剖结构的连接性。同时GCA可反映因果方向性，这为脑神经网络回路的研究提供了方便，这是普通脑功能连接强度检测方法所无法做到的。另外，以往对于同一脑激活区所包含多个体素的时间序列，采取简单平均处理得到一个平均的时间序列，这样会丢掉大量有用脑活动信号信息。本发明采用主成分分析，对多体素的时间序列进行空间上的降维，尽量消除各时间点上所有体素脑活动信号值所构成的向量之间的相关性，能够利用脑活动信号信息的主要成分，有效地提高磁共振脑功能图像因果连接强度的检测性能，更完整更准确地检测脑激活区之间的因果连接。

发明内容

为了克服已有技术的不足，本发明的目的在于设计一种提高磁共振脑成像因果连接强度的检测结果准确性及完整性的检测方法。

为实现上述目的，本发明提出一种基于主成分分析PCA和格兰杰因果分析GCA的磁共振脑成像因果连接强度的检测方法，包括以下步骤：

步骤Sa：对经过数据预处理的脑功能图像进行激活区域多体素时间序列的提取，得到激活区域多个体素的时间序列矩阵；

步骤Sb：使用PCA对每个激活区多个体素的时间序列矩阵进行空间降维得到保留大部分信息量的主要成分，再平均主要成分的激活值得到一个时间序列；

步骤Sc：建立所有激活区域的时间序列间的多变量自回归模型；

步骤Sd：依据所述多变量自回归模型，计算各个时间序列间的偏相关系数；

步骤Se：依据偏相关系数，通过直接传递函数DTF方法计算dDTF值，得到激活脑区间的因果连接强度和方向；

步骤Sf：使用替代数据法统计检验连接强度的显著性，并利用有向网络图显示显著的连接。

本发明针对磁共振脑成像因果连接强度的检测，通过激活区多体素时间序列矩阵提取、PCA降维、使用基于多变量自回归模型及传递函数方法的格兰杰因果分析计算激活区之间的因果连接强度，再用替代数据的方法检验连接强度的统计显著性，提高了因果连接强度检测的准确性和鲁棒性。真实数据试验结果证明该算法比传统采用基于简单平均激活区所有体素激活值的GCA方法能更完整更准确地检测到脑激活区之间的因果连接，为脑功能数据分析方法提供了一种新途径。

附图说明

图1是本发明算法计算流程示意图。

图2是根据本发明算法的实施例获得的因果连接强度图。

图3是使用基于简单平均体素方法获得的因果连接强度图。

具体实施方式

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

参照图1，本发明所涉及的一种磁共振脑成像因果连接强度的检测方法，尤其涉及使用主成分分析和格兰杰因果分析进行磁共振脑成像因果连接强度的检测的具体实施步骤如下：

步骤Sa，对经过数据预处理的脑功能图像进行激活区域多体素时间序列的提取，得到激活区域多个体素的时间序列矩阵；

1.数据预处理

由于磁共振扫描过程中各种各样的噪声的影响，个体自身存在尺度和位置上的差异，非常有必要在分析数据之前对数据做一定的预处理。在整个的实验的数据获取中，主要的噪声信息来源有：(1)物理头动；(2)图像内层间扫描时间差别；(3)外在磁场的不均匀性等。脑功能图像预处理是在保留脑功能图像细节的同时，使用脑功能图像与标准模板进行仿射配准变换方式的预处理，并提高脑功能图像的信噪比。常见的预处理步骤有：切片扫描时间对齐，图像序列对齐，联合配准，标准化(或称均一化)，空间平滑滤波和时间平滑滤波等。

2.独立成分分离

一般来说，分离独立成分的方法有基于广义线性模型(General LinearModels，GLM)的方法、模式分类方法和独立成分分析方法(Independentcomponent analysis，ICA)。其中，基于ICA的分离独立成分模型如下：令v为脑功能观测量(维度)，其由多个信源S经混合矩阵A组合而成(v＝AS)。现在的任务是：在S与A均未知的条件下，求取解混矩阵W，使得v通过W后所得输出y(y＝W·v)是S的最优逼近。

本发明采用独立成分分析方法来分离独立成分，分离步骤如下：

假设分离出的各分量互相独立，由中心极限定理可知，一定条件下相互独立的随机变量之和趋向于高斯分布，所以两个相互独立的随机变量之和比任何一个参与求和的随机变量更加接近于高斯分布，这样最后求出的非高斯性最大的组合随机变量，即距离高斯分布最远的随机变量y就可以认为是源信号v的一个独立成份，这种使距离高斯分布最远的随机变量之间相互独立的方法称为非高斯最大化方法。另外，还可用很多其它方法来定义这种接近程度，如基于信息论的最小化互信息法和最大似然估计法也是对ICA进行估计的方法。并由此发展出很多ICA算法，如最大(小)化一些相关的代价函数，基于随机梯度的自适应算法，还有现在常被采用的快速ICA算法以及Infomax算法等。在处理fMRI数据时，我们一般采用Infomax算法寻找解混矩阵，这样也就能有效地分离多个非高斯性最大的源信号。

Infomax算法是一种基于信息最大化传输的单层前反馈神经网络算法，其基本思想是利用了似然函数极大法估计中的一个准则，即估计随机变量v独立性的最大化等价于估计v的某一非线性函数独立性的最大化。它的主要特点是：引入了非线性函数f(·)，使得新向量y＝f(u)的熵极大化等价于u的各分量独立性的最大化，其中u为输入v各分量的线性组合。在该算法中选用Sigmoid函数作为f(·)，因为Sigmoid函数具有超高斯信号pdf(概率密度函数)的性质，所以能有效地分离多个超高斯分布的源信号。

Sigmoid函数的形式为：f(x)＝(1+e^-x)^-1。

所谓信息最大化传输(熵极大化)就是使输入v与输出y的互信息最大，互信息定义为：

I(y，v)＝H(y)-H(y/v)，

其中H(y)是输出信号y的边缘熵，H(y/v)为条件熵。

为求互信息I(y，v)的最大值，将等式两边对W求导

$\frac{\∂}{\∂W}I(y,v)=\frac{\∂}{\∂W}H\left(y\right),$

其中，条件熵H(y/v)不依赖于解混矩阵W。

输出熵H(y)定义为

$H\left(y\right)=-E\left[\ln f_y\left(y\right)\right]=-\underset{-\∞}{\overset{+\∞}{\∫}}f\left(y\right)\ln f_y\left(y\right)\mathrm{dy},$

其中，E表示求取数学期望，f_y(y)为y的概率密度，它可以由输入v的概率密度f_v(v)得到

$f_y\left(y\right)=\frac{f_v\left(v\right)}{\left|J\right|},$

其中，J为雅可比行列式。

根据以上两式，可得

H(y)＝E[ln|J|]-E[ln f_v(v)]，

其中，E[ln f_v(v)]不依赖于解混矩阵W。得到H(y)后，代入之前的等式

$\frac{\∂}{\∂W}I(y,v)=\frac{\∂}{\∂W}H\left(y\right),$

可得

$\frac{\∂H\left(y\right)}{\∂W}=\frac{\∂I(y,v)}{\∂W}=\frac{\∂}{\∂W}E\left[\ln \right|J\left|\right],$

为了实时分离各独立源，用瞬时值代替均值，可得

$\frac{\∂H\left(y\right)}{\∂W}=\frac{\∂I(y,v)}{\∂W}=\frac{\∂}{\∂W}\ln \left|J\right|,$

选用sigmoid函数作为非线性单元，即f(u)＝(1+e^-u)^-1，y＝f(u)，u＝Wv+w₀，其中w₀为设定的随机初始向量，u为输入v各分量的线性组合，将上式带入

I(y，v)＝H(y)-H(y/v)，

$\frac{\∂H\left(y\right)}{\∂W}=\frac{\∂I(y,v)}{\∂W}=\frac{\∂}{\∂W}\ln \left|J\right|,$

可得迭代增量ΔW：

ΔW ∝{(W^T)^-1+(1-2y)v^T}，

其中，∝表示正比例关系，我们用自然梯度代替常规梯度来避免矩阵求逆，将上式右乘W，得到：

ΔW＝μ[I+(1-2y)u^T]W，

其中，μ为矫正系数，取值在0.9～1之间。

判断Infomax算法是否收敛的性能函数为：

$E=\underset{i=1}{\overset{a}{\mathrm{\Σ}}}(\underset{j=1}{\overset{a}{\mathrm{\Σ}}}\frac{\left|p_{\mathrm{ij}}\right|}{{\max }_k\left|p_{\mathrm{ik}}\right|}-1)+\underset{j=1}{\overset{a}{\mathrm{\Σ}}}(\underset{i=1}{\overset{a}{\mathrm{\Σ}}}\frac{\left|p_{\mathrm{ij}}\right|}{{\max }_k\left|p_{\mathrm{kj}}\right|}-1),$

其中，P＝(p_ij)＝WA为分离矩阵与混合矩阵的乘积，称为全局转换矩阵，p_ij为矩阵元素，i＝1，2，...，a，j＝1，2，...，a，k＝1，2，...，a，a为样本总数。

其算法的主要步骤如下：

(1)初始化解混矩阵W(随机数)；

(2)取第i个样本矢量v_i，i＝1，2，…，a，a为样本总数；

1.计算解向量u_i＝Wv_i和网络输出y_i＝[1+exp(-u_i)]^-1，

2.计算权值增量， $\mathrm{\ΔW}=\mathrm{\μ}\frac{\∂H\left(y\right)}{\∂W}=\mathrm{\μ}[I+(1-2y)u^T]W,$ 其中y＝(y₁，y₂，…，y_a)^T，u＝(u₁，u₂，…，u_a)^T，T为矩阵转置符号。

3.更新权值W_t＝W_t-1+ΔW，其中t＝1，2，…，end，end为迭代次数；

(3)是否达到收敛条件E＜ε，ε为无穷小量，约等于0，如取ε＝10^-10，是则结束；否则回到第(2)步。

如果达到收敛的条件，那么就可以得到y＝(y₁，y₂，…，y_a)^T，从而得到分离出的独立成分y₁，y₂，…，y_a。

3.激活区多体素时间序列提取

依据ICA对每个被测试者分离出独立成分y₁，y₂，…，y_a后，在总体独立成分数据上进行单样本t检验，获得平均成分。对每个平均成分中每个体素的时间序列求出方差值，作为此体素的激活值(T值)。通过对比这些平均成分的不同的激活脑区与已知网络成分中的激活脑区，选择出我们感兴趣的成分，然后从挑选出的感兴趣的成分中确定出激活值(T值)最大点体素的坐标(此坐标为标准脑模板的坐标)，其中，每个体素的激活值(T值)为此体素在不同时间点上形成的时间序列的方差值，并以此坐标为中心作一个半径为3倍体素棱长(即采集图像的分辨率，一般为3mm×3mm×3mm)的球，该球中的体素均属于同一激活区。依据预处理后的脑功能图像数据，提取包含于此球体内部的各个体素在不同时间点上激活值的时间序列矩阵Y′(矩阵维数为D×N)，其中D为包含于球体内部的体素数目，N为时间点数。

步骤Sb，使用PCA对每个球体(激活区)中的多个体素的时间序列矩阵进行空间降维，得到保留大部分信息量的主要成分，再平均主要成分的激活值即可得到一个时间序列；

依据激活区内大部分的信息和能量包含于少数的主要成分的原则，对球体内所有激活体素的激活值进行协方差分析来降低体素的空间维数，即把激活体素的时间序列矩阵的空间维数从Y′(矩阵维数为D×N)降低到X′(矩阵维数为d×N)，其中，Y′为提取出的包含于球体内部的各个体素的时间序列矩阵，D为球体内激活体素个数，N为时间点个数，X′为降维后的时间序列矩阵，d为保留的主成分个数。根据矩阵分解原理，存在非奇异变换，使得X′＝(W′)^TY′，其中W′为D×d的矩阵，上标T表示矩阵转置。

1.计算球体内的所有激活体素的相关系数矩阵

为保证数据分析结果的准确性，使用总和标准化法标准化具有D个体素的时间序列矩阵Y′，即，对矩阵每行求和得到[s₁，…，s_D]^T，矩阵Y′每行各个元素除以其相应的行和s_i(i＝1，…，D)，得到新矩阵Y″，计算Y″的相关系数矩阵ρ＝(ρ_ij)(矩阵ρ的维数为D×D)，

${\mathrm{\ρ}}_{\mathrm{ij}}=\frac{\mathrm{Cov}(Y_i^{\′\′},Y_j^{\′\′})}{\sqrt{D\left(Y_i^{\′\′}\right)}\sqrt{D\left(Y_j^{\′\′}\right)}},$

其中，ρ_ij为相关系数矩阵元素，Cov()表示两时间序列的协方差，Y″_i，Y″_j分别表示矩阵Y″的第i，j行矩阵元素，D()表示每行矩阵元素的方差。

2.主成分排序

为了实现主成分排序，计算相关系数矩阵ρ的特征值，并将特征值按从大到小的顺序排序：λ₁≥λ₂≥…≥λ_D，其相应的特征向量为v₁，v₂，…，v_D，这样就完成了对成分的排序，排在越前面的特征向量在整个数据分析中承担的主要意义越大。

3.计算累计贡献率

贡献率表示所定义的主成分在整个数据分析中承担的主要意义占多大的比重，当取前d个主成分来代替原来全部变量时，累计贡献率的大小反应了这种取代的可靠性，累计贡献率越大，可靠性越大；反之，则可靠性越小。一般要求累计贡献率达到85％，即

$\frac{\underset{i=1}{\overset{d}{\mathrm{\Σ}}}{\mathrm{\λ}}_i}{\underset{i=1}{\overset{D}{\mathrm{\Σ}}}{\mathrm{\λ}}_i}\≥85\%.$

选择前d个特征值对应的特征向量构成矩阵W′＝[v₁，…，v_d]，再由X′＝(W′)^TY′获得由PCA降维后的主成分时间序列矩阵X′(维数为d×N)。平均该激活区的主成分时间序列矩阵X′的激活值就可得到该激活区的一个时间序列x_i(t)(i＝1，…，m)，其中，m为所有激活区的数目。

步骤Sc，建立所有激活区的时间序列间的多变量自回归模型(MVAR)，并确定模型的阶数p；

建立时间序列间的多变量自回归模型：

$X\left(t\right)-\underset{n=1}{\overset{p}{\mathrm{\Σ}}}A\left(n\right)X(t-n)=E\left(t\right),$

其中，X(t)＝(x₁(t)，x₂(t)，...，x_m(t))^T，x_i(t)(i＝1，…，m)为第i个激活区经PCA降维后的主成分的平均时间序列，m为所有激活区数目，p为模型阶数，A(n)(n＝1，…，p)为m×m的系数矩阵，E(t)为零均值白噪声。

模型的阶数p依靠Akaike information criterion(AIC)准则来决定，

AIC(p)＝2log(det(∑))+2pm²/N，

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

其中，AIC(p)表示AIC统计量，R(n)＝X(t)X^T(t+n)是X(t)步长为n的协方差矩阵，∑为噪声协方差阵，m为激活区数目，N为时间点数目。

随着模型的阶数逐渐增大，AIC(p)取最小的统计量时，此时的p即选求模型的阶。

步骤Sd，依据上一步骤确定的多变量自回归模型，计算各个时间序列间的偏相关系数；

1.将多变量自回归模型从时域转换到频域，以得到频率转移矩阵。

为了得到时间序列在频域上的因果关系，利用Fourier变换将步骤Sc中的时域多变量自回归模型变换到频域：

$X\left(f\right)[{\mathrm{\δ}}_{\mathrm{ij}}-\underset{n=1}{\overset{p}{\mathrm{\Σ}}}{a}_{\mathrm{ij}}\left(n\right)e^{-i2\mathrm{\πfn}}]=E\left(f\right),$

${\mathrm{\δ}}_{\mathrm{ij}}=\left\{\begin{array}{cc}1,& i=j\\ 0,& i\≠j.\end{array}\right.,$

其中，f表示频率，X(f)为X(t)在频域上的表示，δ_ij为狄拉克δ函数，i＝1，…，m，j＝1，…，m，m为所有激活区数目，a_ij(n)为矩阵A(n)的元素，E(f)为E(t)在频域上的表示。

令其中，a_ij(f)为矩阵A(f)的元素，A(f)为A(n)在频域上的表示，则有

X(f)＝A^-1(f)E(f)＝H(f)E(f)，

其中，H(f)＝A^-1(f)为频域转移矩阵，H(f)包含了时间序列间相互作用的所有信息，其矩阵元素h_ij(f)表示第j个激活区对第i个激活区的非标准化的DTF值，其中，DTF值为直接传递函数值，表征了任意两个激活区之间的因果关系的强度。

2.计算偏相关系数

为了获得以频率的函数表示的第i个激活区与第j个激活区时间序列间的关系，即描述两个时间序列相互关联程度的数字特征，计算交叉频谱S(f)：

S(f)＝H(f)VH^H(f)，

然后计算时间序列间的偏相关系数θ_ij(f)：

${\mathrm{\θ}}_{\mathrm{ij}}\left(f\right)=\frac{M_{\mathrm{ij}}^2\left(f\right)}{M_{\mathrm{ii}}\left(f\right)M_{\mathrm{jj}}\left(f\right)},$

其中，V表示零均值白噪声矩阵E(t)的方差，H(f)为频域转移矩阵，上标H表示共轭转置，θ_ij(f)表示第i个激活区与第j个激活区时间序列间的偏相关系数，M_ij(f)表示交叉频谱S(f)的代数余子式。

步骤Se，依据偏相关系数，通过直接传递函数(DTF)方法计算dDTF值，得到激活脑区间的因果连接强度和方向；

为得到激活脑区间的因果连接强度和方向，通过直接传递函数(DTF)方法计算dDTF值：

${\mathrm{dDTF}}_{\mathrm{ij}}=\underset{f}{\mathrm{\Σ}}{h}_{\mathrm{ij}}\left(f\right){\mathrm{\θ}}_{\mathrm{ij}}\left(f\right),$

其中，dDTF_ij值表示激活脑区j到激活脑区i的因果连接强度值，h_ij(f)表示频域转移矩阵H(f)的元素，θ_ij(f)表示时间序列间的偏相关系数。dDTF_ij值表示了非标准化的DTF值与偏相关系数的乘积在所有频率成分上的总和，强调了激活区之间直接的因果联系。

步骤Sf，使用替代数据法统计检验连接强度的显著性，并利用有向网络图显示显著的连接。

DTF函数(即频域转移矩阵H(f))与生成的时间序列之间有很高的非线性关系，不容易建立dDTF值的估计量的分布函数。为了检验由dDTF值表征的格兰杰因果关系的显著性，使用替代数据生成经验假设分布的方法进行统计检验。

为保证生成的替代数据与原有的时间序列具有相同的频谱，将原始的时间序列使用Fourier变换到频域，相位平均分布在(-π，π)，随机抽取在此范围的相位2500次并使之服从均匀分布，抽取一次即得到一条替代数据组成的新时间序列，然后再将在频域中抽取2500次所抽取出的时间序列变换回时域得到作为替代数据的时间序列。每次得到作为替代数据的时间序列之后重新计算激活区之间的dDTF值。由每个被测试者的所有激活区之间的可能因果连接产生一个假设的正态分布，对于每个被测试者的每个连接，用真实的dDTF值与它相关的假设正态分布作假设检验得到p值。p值是将观察结果认为有效即具有总体代表性的犯错概率，在这里即使真实的dDTF值与它相关的假设正态分布的均值相等，也会在p值的概率下不相等，即犯错的概率为p。为了得到组分析的结果，假设所有被测试者的同一个连接的p值是相互独立的，则统计量 $X^2=-2\underset{i=1}{\overset{K}{\mathrm{\Σ}}}{\log }_e\left(p_{\mathrm{ij}}\right),$ (j＝1，…，m)服从自由度为2K的卡方分布。

其中，K为被测试者数目，p_ij为第i个被测试者的第j个连接的p值；统计量X²大于查表值X²(2K)的连接则视为显著的连接，最后利用有向网络图显示显著的连接。

本发明所述的因果连接强度的检测方法，可通过真实的磁共振脑成像数据得以说明：

(1)真实数据实验

为展示本发明的效果，在实施方案中采用真实数据集作测试，共12个被测试者参与了实验，9男3女。实验中采用针刺刺激一个中医传统穴位丘墟，用于治疗听觉相关的疾病。实验过程由专业的针灸师操作，采用T2*加权梯度回波平面成像(Echo-Planar Imaging，EPI)序列获取针刺刺激后BOLD fMRI静息数据。

采用SPM5(http://www.fil.ion.ucl.ac.uk/spm/)对数据进行预处理，包括切片扫描时间对齐，图像序列对齐，联合配准，标准化(或称均一化)，空间平滑滤波。然后进行激活区提取，利用GIFT(http://icatb.sourceforge.net/)处理，包括估计独立成分，寻找解混矩阵，重建信号。分离出默认脑网络成分和听觉网络成分，并得到激活区，前脑岛(anterior insula，AI)和颞上回(superior temporal gyrus，STG)属于听觉网络成分；前扣带回(anteriorcingulate cortex，ACC)，颞下回(inferior temporal gyrus，ITG)和后扣带回(posterior cingulate cortex，PCC)属于默认脑网络成分。提取这五个激活区的时间序列，并用本发明所述方法(方法A)对比现有的基于简单平均激活区所有体素激活值的GCA方法(方法B)，获得激活区之间的因果连接强度和方向的有向网络图，并将激活区的因果连接强度作为算法性能的度量。

(2)实验结果

在真实实验数据集上两种方法的激活脑区之间的显著的因果连接强度及方向如表格1和表格2所示。如图2(方法A)和图3(方法B)所示箭头表示连接强度的方向，由表格中每列所对应的激活区指向每行所对应的激活区，连接强度由线条的宽度表示。方法A得到的因果连接强度显著高于方法B得到因果连接强度，如AI到ITG的因果连接强度0.021(方法A)＞0.012(方法B)，STG到AI的因果连接强度0.267(方法A)＞0.206(方法B)等。方法A中出现的STG到ACC的微弱因果连接强度，且方法A中出现的ITG与ACC之间的微弱的双向因果连接在方法B中均未被检测到。

以上实验结果说明，本发明所述的一种基于PCA和GCA的磁共振脑成像因果连接强度的检测方法，有效地提高了磁共振脑功能图像因果连接强度的检测性能，更完整更准确地检测脑激活区之间的因果连接。

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

表格2基于简单平均激活区所有体素激活值的GCA方法(B)的结果

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种基于主成分分析PCA和格兰杰因果分析GCA的磁共振脑成像因果连接强度的检测方法，其特征在于，包括以下步骤：

步骤Sa：对经过数据预处理的脑功能图像进行激活区域多体素时间序列的提取，得到激活区域多个体素的时间序列矩阵；

步骤Sb：使用PCA对每个激活区多个体素的时间序列矩阵进行空间降维得到保留大部分信息量的主要成分，再平均主要成分的激活值得到一个时间序列；

步骤Sc：建立所有激活区域的时间序列间的多变量自回归模型；

步骤Sd：依据所述多变量自回归模型，计算各个时间序列间的偏相关系数；

步骤Se：依据偏相关系数，通过直接传递函数DTF方法计算dDTF值，得到激活脑区间的因果连接强度和方向；

步骤Sf：使用替代数据法统计检验连接强度的显著性，并利用有向网络图显示显著的连接。

2.如权利要求1所述的检测方法，其特征在于，所述脑功能图像的数据预处理是在保留脑功能图像细节的同时，使用脑功能图像与标准模板进行仿射配准变换方式的预处理。

3.如权利要求1所述的检测方法，其特征在于，提取激活区域多体素时间序列的步骤为：

(1)利用独立成分分析ICA方法得到不同激活模式的独立成分；

(2)获得所述独立成分中激活值最大点体素的坐标，以此坐标为中心作一个半径为3倍体素棱长的球体；

(3)依据预处理后的脑功能图像数据，提取包含于此球体内部的各个体素在不同时间点上激活值的时间序列。

4.如权利要求3所述的检测方法，其特征在于，所述获得独立成分中激活值最大点体素的坐标的步骤具体包括：

(1)在总体独立成分数据上进行单样本t检验，获得平均成分；

(2)对每个平均成分中每个体素的时间序列求得方差值，作为此体素的激活值；

(3)从感兴趣的成分中确定出激活值最大点体素的坐标。

5.如权利要求1所述的检测方法，其特征在于，使用PCA对每个激活区多个体素的时间序列矩阵进行空间降维得到保留大部分信息量的主要成分，再平均主成分得到一个时间序列的步骤为：

(1)使用总和标准化法标准化具有D个体素的时间序列矩阵Y′，得到新矩阵Y″；

(2)计算Y″的相关系数矩阵ρ＝(ρ_ij)，矩阵ρ的维数为D×D，

${\mathrm{\ρ}}_{\mathrm{ij}}=\frac{\mathrm{Cov}(Y_i^{\′\′},Y_j^{\′\′})}{\sqrt{D\left(Y_i^{\′\′}\right)}\sqrt{D\left(Y_j^{\′\′}\right)}},$

其中，ρ_ij为相关系数矩阵元素，Cov()表示两时间序列的协方差，Y″_i，Y″_j分别表示矩阵Y″的第i，j行元素，D()表示每行矩阵元素的方差；

(3)计算相关系数矩阵ρ的特征值，并将特征值按从大到小的顺序排序：λ₁≥λ₂≥…≥λ_D，其相应的特征向量为v₁，v₂，…，v_D；

(4)计算累计贡献率达到85％的前d个特征值，其中累计贡献率按下式计算：

$\frac{\underset{i=1}{\overset{d}{\mathrm{\Σ}}}{\mathrm{\λ}}_i}{\underset{i=1}{\overset{D}{\mathrm{\Σ}}}{\mathrm{\λ}}_i}\≥85\%;$

(5)取所述前d个特征值对应的特征向量构成矩阵W′＝[v₁，…，v_d]；

(6)由非奇异变换X′＝(W′)^TY′获得降维后的时间序列矩阵X′；

(7)平均X′的激活值，得到一个时间序列x_i(t)(i＝1，…，m)，其中，m为所有激活区数目。

6.如权利要求1所述的检测方法，其特征在于，所有激活区的时间序列间的多变量自回归模型为：

$X\left(t\right)-\underset{n=1}{\overset{p}{\mathrm{\Σ}}}A\left(n\right)X(t-n)=E\left(t\right),$

其中，X(t)＝(x₁(t)，x₂(t)，...，x_m(t))^T，x_i(t)(i＝1，…，m)为第i个激活区经PCA降维后的主成分的平均时间序列，p为模型阶数，A(n)(n＝1，…，p)为m×m的系数矩阵，E(t)为零均值白噪声。

7.如权利要求6所述的检测方法，其特征在于，步骤Sd具体包括以下步骤：

(1)将步骤Sc中得到的多变量自回归模型从时域转换到频域，得到频率转移矩阵H(f)：

X(f)＝A^-1(f)E(f)＝H(f)E(f)，

其中，f表示频率，X(f)为X(t)在频域上的表示，A(f)为A(n)在频域上的表示，E(f)为E(t)在频域上的表示，H(f)＝A^-1(f)为频域转移矩阵；

(2)计算各个时间序列间的偏相关系数θ_ij(f)。

8.如权利要求7所述的检测方法，其特征在于，所述偏相关系数θ_ij(f)的计算包括：

(1)计算交叉频谱S(f)：

S(f)＝H(f)VH^H(f)，

其中，V表示零均值白噪声矩阵E(t)的方差，上标H表示共轭转置；

(2)计算偏相关系数θ_ij(f)：

${\mathrm{\θ}}_{\mathrm{ij}}\left(f\right)=\frac{M_{\mathrm{ij}}^2\left(f\right)}{M_{\mathrm{ii}}\left(f\right)M_{\mathrm{jj}}\left(f\right)},$

其中，θ_ij(f)表示第i个激活区与第j个激活区时间序列间的偏相关系数，M_ij(f)表示交叉频谱S(f)的代数余子式。

9.如权利要求1所述的检测方法，其特征在于，步骤Se中通过直接传递函数DTF方法计算dDTF值具体为：

${\mathrm{dDTF}}_{\mathrm{ij}}=\underset{f}{\mathrm{\Σ}}{h}_{\mathrm{ij}}\left(f\right){\mathrm{\θ}}_{\mathrm{ij}}\left(f\right),$

其中，dDTF_ij表示激活脑区j到激活脑区i的因果连接强度值，h_ij(f)表示频域转移矩阵H(f)的元素，θ_ij(f)表示时间序列间的偏相关系数。

10.如权利要求1所述的检测方法，其特征在于，步骤Sf中使用替代数据的方法进行统计检验连接强度的显著性，并利用有向网络图显示显著的连接。

Images