CN102508184B - 一种基于移动平均时间序列模型的脑功能激活区检测方法 - Google Patents

Abstract

本发明公开了一种基于移动平均时间序列模型的磁共振成像脑激活区检测方法，该方法包括以下步骤：使用移动平均算法，对于每个个体，逐体素计算个体(individual subject)体素信号强度时间序列模型；在此基础上使用EM-algorithm计算每个个体在组群中的权重值，逐体素构建组群(group subjects)体素信号强度时间序列动态变化模型；最后通过区域增长方法，度量组群体素间信号强度时间序列模型的活动相似性，确定脑功能激活区。在真实fMRI数据集上的实验说明，本发明所述方法是一种有效的脑功能激活区提取方法。

Description

一种基于移动平均时间序列模型的脑功能激活区检测方法

技术领域

本发明属于生物信息、图像识别技术领域，涉及一种基于时间序列状态分析的核磁共振成像(functional Magnetic Resonance Imaging，fMRI)图像后处理技术。

背景技术

功能磁共振成像(fMRI)技术是在传统磁共振成像(Magnetic ResonanceImaging，MRI)技术的基础上发展而来。fMRI技术基于血氧水平依赖(BloodOxygenation Level-Dependent，BOLD)，通过测量由脑血流或脑血氧变化而引起的磁共振信号变化来反映脑活动。由于其具有无创性、可重复性等优点，并可提供较高的图像空间分辨率，该技术已被广泛应用于神经科学、临床医学等领域的研究。作为一种非侵入性检测手段，fMRI可为一些重要脑功能区(如运动、语言等)定位。脑功能定位研究是神经科学研究发展极为迅速的领域之一，为人类的认知与思维活动及疾病的康复治疗与预测带来了新的前景。目前，脑功能定位研究已广泛应用于神经外科、认知功能相关研究、康复治疗与预测、社会心理等方面。随着基础研究与临床方面的应用深入，合理有效的分析fMRI数据并提取与外部刺激相关的脑激活区对于脑功能活动区的精确定位愈发重要。

现有fMRI脑功能激活区检测方法大多基于广义线性模型(GeneralLinear Model，GLM)方法。该模型本质上是一种模型驱动(Data Driven)的一元回归分析方法，一般先指定一个先验模型(即一个包含实验设计、人为增加的心跳、呼吸、头动等因素的模型)，将该模型分别与大脑中的每一个体素的时间序列进行匹配，确定每一个体素的“激活”状态。因此，该方法的应用前提是fMRI时间序列中“激活”体素的信号强度时间序列变化需符合实验设计中的指定的先验模型。然而，越来越多的影像学研究及临床应用表明，神经激活信号的时间变化模式并无准确先验。因此，基于GLM模型的先验假设过于单一，并不能够很好地表示fMRI数据。脑活动的分布模式、血液动力学模型等不能由GLM模型恰如其分地模拟。

发明内容

为了避免现有方法的不足，本发明设计一种基于移动平均时间序列模型的脑功能激活区检测方法，以克服现有方法依赖先验模型的缺点，达到有效提取脑激活区目的。

为实现上述目的，本发明提出的一种基于移动平均时间序列模型的磁共振成像脑激活区检测方法，其特征在于，该方法包括以下步骤：

步骤1，对采集到的原始功能磁共振成像(fMRI)图像进行预处理；

步骤2，对于每个被测试者逐体素计算每个被测试者的个体体素信号强度时间序列模型；

步骤3，对于组群，计算每个个体在组群中的权重值，并根据全脑个体体素信号强度时间序列模型，逐体素构建组群体素信号强度时间序列模型；

步骤4，根据全脑组群体素信号强度时间序列模型，应用区域增长方法实现脑激活区的提取。

本发明针对现有脑影像fMRI数据分析方法过于依赖先验模型这一问题，采用移动平均方法构建大脑中每个体素的时间序列状态变化模型，该时间序列动态变化模型不依赖与先验知识；通过区域增长方法度量各体素时间序列动态变化的活动相似性，可自适应地提取脑功能激活区。将该方法应用于fMRI影像数据集的结果亦表明，该方法是一种有效的脑功能激活区提取方法。

附图说明

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

图2是本发明中个体体素信号强度时间序列动态变化模型。

图3是本发明中组群体素信号强度时间序列动态变化模型。

图4是本发明中提取的脑功能区激活图。

具体实施方式

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

参照图1，本发明提出的一种基于移动平均时间序列模型的脑功能激活区检测方法，具体包括以下步骤：

步骤1、对采集到的原始fMRI图像进行预处理。

首先对输入的三维fMRI图像进行时间点校正和头动校正；再将校正后的图像进行标准化处理：采用SPM(Statistical Parametric Mapping，SPM，统计参数映射)软件自带的EPI(Echo Planar Image，平面回波成像)模板，将校正后的图像归一化到蒙特利尔神经科学研究所(Montreal NeurologicalInstitute，MNI)标准脑模板上，并把三维fMRI图像中的每个体素重新采样至3×3×3mm³；然后对标准化后的图像使用高斯核函数进行平滑，其中，高斯核函数的全宽半高值(FWHM)为6×6×6mm³；最后对平滑后的图像去除线性漂移，将扫描基线的漂移进行去除。

步骤2、对于每个被测试者(individual subject)，逐体素计算每个被测试者的个体(individual subject)体素信号强度时间序列模型。

步骤2进一步包括以下几个步骤：

步骤2.1，构建个体体素信号强度时间序列预测统计量；

对于每个体素，假设第i个个体体素的fMRI神经响应的时间序列分布状态X满足以下两种数据分布形式：基线状态X₀～N(θ₀，∑_i)，以及激活状态X₁～N(θ₁，∑_i)，两种状态分布的均值不同，分别对应于无刺激状态和刺激状态。其中，θ₀，θ₁为均值，∑_i为协方差。那么，实测的信号强度时间序列统计模型为：

$x_t^i=s_t^i+{\mathrm{\ϵ}}_t^i,i=1,...,m,t=1,...,n$

其中，t为时间序列取值的时间点，为第i个个体在时刻t测得的信号强度，为第i个个体在时刻t所对应的噪音项，为第i个个体对应于整个时间序列的噪音向量，服从均值为0，协方差为∑_i的正态分布，m为所有被测试者人数，n为时间点个数，表示第i个个体在t时刻的信号强度期望值，满足：

$s_t^i=\left\{\begin{array}{c}{\mathrm{\θ}}_0,t=1,...,\mathrm{\τ}\\ {\mathrm{\θ}}_1,t=\mathrm{\τ}+1,...,n\end{array}\right.,$

其中，τ为从基线状态到激活状态的状态改变所对应的时间点，在τ时刻，fMRI时间序列的状态发生改变，如由外界刺激条件引起大脑皮层神经活动的增强(或减弱)，进而使时间序列中的数据从均值为θ₀转变到均值为θ₁的分布状态。

利用移动平均法构建个体体素信号强度时间序列预测统计量满足：

$z_t^i=\mathrm{\λ}{x}_t^i+(1-\mathrm{\λ})z_{t-1}^i,\∀1\≤t\≤n,$

其中，λ为权重系数，0≤λ≤1，λ决定了该模型检测fMRI时间序列信号突然发生变化的能力，即时效性。λ越大，时效性越强，反之越弱。这里选择λ＝0.2。

步骤2.2，基于个体体素信号强度时间序列预测统计量，统计推断个体体素信号强度时间序列变化状态，完成个体体素信号强度时间序列模型的建立。

通过对数据分布形式加以判断，可以判断特定时刻时间序列信号的状态是否发生变化。为此做出如下假设：零假设(H₀)为只存在一种均值状态，均值μ为常数，即信号状态没有发生改变；备择假设(H₁)为存在两种状态(基线状态与激活状态)，均值μ对应不同值，即：

H₀：μ＝θ₀，t＝1，...，n；

H₁：μ＝θ₀，t＝1，...，τ；

    μ＝θ₁，t＝τ+1，...，n；

为了检验个体体素信号强度的状态是否发生变化，定义检验统计量T为：

$T_t^i=\frac{z_t^i-{\mathrm{\θ}}_0}{\sqrt{\mathrm{Var}\left(z_t^i\right)}},$

其中，为统计量在t时刻的方差。

检验统计量T服从自由度为d的t分布t(∑，df)，控制阈限(ControlLimits，即置信度区间)定义为：

${\mathrm{\θ}}_0\±t^{*}\sqrt{\mathrm{Var}\left(z_t^i\right)},$

其中，t^*为假阳性率(一般取值为95％)。一旦统计量超过该阈限值，则表明t时刻该体素处于“激活”状态。

由于检验统计量T服从自由度为d的t分布t(∑，df)，因此fMRI时间序列中每个时间点t都对应着一个概率p值。为遍历整个时间序列，需对p值进行校正。令为对应于所有时间点的检验统计量向量，服从多变量t分布，即：优选地，可采用蒙特-卡罗(Monte Carlointegration)方法在时间维度上对p值进行FWER(Family-wise Error Rrate)校正，以估计t^*值。具体来说，可由t(∑，df)分布随机生成2000个长度为n的特征向量，用T_k表示其中一个特征向量。对于T_k，可求得T_k各分量的绝对值最大，表示为：

$T_{\max }=\underset{k\∈\{1,\·\·\·,n\}}{\max }\left[\right|T_k\left|\right\}.$

因此，对应于随机生成的2000个特征向量，可求得2000个T_max值。将T_max按照由大到小的顺序排序：T_max，1＞T_max，2＞…＞T_max，2000，由于蒙特-卡罗方法的校正过程本质上是多重检验的过程，一般通过控制假阳性率(FDR，False Discovery Rate)来决定p值。实践过程中FDR的值一般选取0.05。因此t^*取为T_max，95，t^*即为校正后的p值。

对于每个个体i的任一体素，用表示第i个个体对应于整个时间序列的信号强度预测统计量，模拟了信号强度时间序列动态变化，通过与阈限值比较，可以判断某时刻该体素是否处于“激活”状态。

至此，通过计算个体体素信号强度时间序列预测统计量，并在此基础上估计信号强度时间序列变化状态，完成了个体体素信号强度时间序列模型的建立。

上述描述是针对单一体素举例，将上述步骤应用于全脑体素，可以得到全脑个体体素信号强度时间序列模型。

步骤3、对于组群(group subjects)，计算每个个体在组群中的权重值，并根据全脑个体体素信号强度时间序列模型，逐体素构建组群(groupsubjects)体素信号强度时间序列模型。

步骤3进一步包括以下几个步骤：

步骤3.1，构建组群体素信号强度时间序列预测统计量；

对于组群，考虑到不同被测试者之间的个体差异，实测的信号时间序列统计模型可写为：

$x_t^i=s_t^i+{\mathrm{\ϵ}}_t^i,i=1,\·\·\·,m$

$s_t^i=s_t^{\mathrm{pop}}+{\mathrm{\η}}_t,t=1,\·\·\·,n$

其中，为第i个个体在时刻t测得的信号强度；为第i个个体在时刻t所对应的噪音项， 表示第i个个体在时刻t的信号强度期望值；η_t为时刻t组群内不同被测试者个体差异引起的噪音项，这里∑_B为不同个体间的噪音协方差矩阵，m为个体个数，n为时间点个数，满足：

$s_t^{\mathrm{pop}}=\left\{\begin{array}{cc}{\mathrm{\θ}}_0,& t=\mathrm{1,2},\·\·\·,\mathrm{\τ}\\ {\mathrm{\θ}}_1,& t=\mathrm{\τ}+1,\·\·\·,n\end{array}\right.,$

将该模型写为矩阵形式，可得：

$\stackrel{{\→}^i}{x}=\stackrel{{\→}^i}{s}+\stackrel{{\→}^i}{\mathrm{\ϵ}}.$

$\stackrel{{\→}^i}{s}={\stackrel{\→}{s}}_{\mathrm{pop}}+\stackrel{\→}{\mathrm{\η}}$

由上面2个公式，该模型又可写为：

$\stackrel{{\→}^i}{x}={\stackrel{\→}{s}}_{\mathrm{pop}}+\stackrel{\→}{\mathrm{\η}}+\stackrel{{\→}^i}{\mathrm{\ϵ}}={\stackrel{\→}{s}}_{\mathrm{pop}}+\stackrel{{\→}^i}{\mathrm{\ζ}},$

其中，用表示第i个个体的总噪音项，令V_i＝∑_i+∑_B，为第i个个体的体素信号强度时间序列预测统计量。则第i个个体信号强度时间序列预测统计量的噪音协方差矩阵可写为：

$V_i^{*}=\mathrm{\Λ}{V}_i{\mathrm{\Λ}}^T={\mathrm{\Σ}}_i^{*}+{\mathrm{\Σ}}_B^{*}$

其中，Λ为：

λ为权重系数，0≤λ≤1，此处取λ＝0.2。第i个个体信号强度时间序列预测统计量的噪音协方差矩阵可看做由2部分组成，分别对应于个体噪音项的部分和对应于由组群内不同被测试者个体差异引起的噪音项部分。

组群体素信号强度时间序列预测统计量可由个体体素信号强度时间序列预测统计量的权重和表示，具有如下形式：

${\stackrel{\→}{z}}_{\mathrm{pop}}={\left[\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{V}_i^{{*}^{-1}}\right]}^{-1}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{V}_i^{{*}^{-1}}\stackrel{{\→}^i}{z},$

其中，为第i个个体的体素信号强度时间序列预测统计量。令上式可写为：

${\stackrel{\→}{z}}_{\mathrm{pop}}=D^{-1}D\stackrel{{\→}^i}{z},$

由此可知，组群体素信号强度时间序列预测统计量可看做由个体体素信号强度时间序列预测统计量权重组合而成，D为权重系数。

步骤3.2，计算每个个体在组群中的权重值D，优选地，采用期望最大化算法计算每个个体在组群中的权重值D；

对于组群内不同被试个体差异引起的噪音项∑_B而言，由于个体间的信号是相互独立的，因此不同个体间的噪音协方差为0，即：∑_B的非对角线元素都为0，因此∑_B可写为如下形式：

${\mathrm{\Σ}}_B=I_n{\mathrm{\σ}}_B^2,$

${\mathrm{\Σ}}_B^{*}=\mathrm{\α}{\mathrm{\Λ}}^T\mathrm{\Λ},$

其中，I_n为行数为n的单位矩阵；为第i个个体噪音协方差项由组群内不同被测试者个体差异引起的部分；为行数为n，列数为1服从的随机变量，的每个分量服从标准正态分布。α代表需要求解的部分，Λ^TΛ代表的已知部分。

因此为了估计权重值D，只需估计α。

优选地，采用期望最大化算法(EM-algorithm)估计α。根据期望最大化算法，定义下式：

$V_G^{*}=\mathrm{\α}{Q}_G+{\mathrm{\Σ}}_G,$

$Q_G=I_m\⊗{\mathrm{\Λ}}^T\mathrm{\Λ},$

其中，Q_G为权重系数λ对应的算法初始化参数项；为第m个个体体素信号强度时间序列预测统计量对应的噪音协方差；∑_G为个体体素信号强度时间序列预测统计量对应的算法初始化参数项；为组群个体体素信号强度时间序列预测统计量对应的噪音协方差矩阵；为行数为m的单位矩阵。定义G＝[I_n，I_n，…，I_n]^T是行数为mn、列数为n的矩阵，其中I_n为行数为n的单位矩阵。

令为一个行数为m，列数为n的矩阵，其中，中的第i个行向量对应于第i个个体的组群体素信号强度时间序列模型。

最大期望算法经过计算期望和最大化两个步骤交替进行计算：

1.计算期望(E步骤)： $V_G^{*}=\mathrm{\α}{Q}_G+{\mathrm{\Σ}}_G$

$C={\left(G^T{V}_G^{{*}^{-1}}G\right)}^{-1}$

2.最大化(M步骤)： $P=V_G^{{*}^{-1}}-V_G^{{*}^{-1}}\mathrm{GC}{G}^T{V}_G^{{*}^{-1}}$

$g=-\frac{1}{2}\mathrm{tr}\left({\mathrm{PQ}}_G\right)+\frac{1}{2}\mathrm{tr}\left(P^T{Q}_{G}P{\stackrel{\→}{z}}_G{\stackrel{\→}{z}}_G^T\right)$

$H=\frac{1}{2}\mathrm{tr}\left(P^T{Q}_{G}P{Q}_G\right)$

α＝α+H^-1g，

重复上面的E步骤和M步骤直到收敛。其中，为行数为m，列数为n的矩阵，第i个行向量为第i个个体的组群体素信号强度时间序列模型；C为数据集的期望值；P、g、H为使得步骤1期望最大化时的参数估计。

步骤3.3，基于组群体素信号强度时间序列预测统计量，统计推断组群体素信号强度时间序列变化状态。

为检测组群时间序列信号的状态是否发生变化，定义假设：零假设为只存在一种均值状态，均值μ为常数，即信号状态没有发生改变；备择假设为存在两种状态(基线状态与激活状态)，均值μ对应不同值，即：

H₀：μ＝θ₀，t＝1，...，n，

H₁：μ＝θ₀，t＝1，...，τ，

    μ＝θ₁，t＝τ+1，...，n，

定义检验统计量T为：

$T_t=\frac{z_{{\mathrm{pop}}_t}-{\mathrm{\θ}}_0}{\sqrt{\mathrm{Var}\left(z_{{\mathrm{pop}}_t}\right)}}$

控制阈限定义为：

${\mathrm{\θ}}_0\±t^{*}\sqrt{\mathrm{Var}\left(z_{{\mathrm{pop}}_t}\right)}$

其中，t^*为假阳性率(一般取值为95％)。一旦统计量超过该阈限值，则表明t时刻该体素处于“激活”状态。

同步骤2.2的处理方法，可采用蒙特-卡罗(Monte Carlo integration)方法在时间维度上对概率p进行FWER(Family-wise Error Rrate)校正，估计t^*值。可由协方差矩阵的对角线元素计算而得。一旦统计量超过该阈限值，则表明t时刻该体素处于“激活”状态。

步骤3.4，基于全脑个体体素信号强度时间序列模型、每个个体在组群中的权重值D以及组群体素信号强度时间序列变化状态，逐体素构建组群体素信号强度时间序列模型。

至此，基于步骤一得到的每个个体体素信号强度时间序列模型，通过EM算法求得每个个体在组群中的权重值，可推断组群体素信号强度时间序列模型。

上述说明针对单一体素举例，将上述步骤应用于全脑体素，可以得到全脑组群体素信号强度时间序列模型。

步骤4、根据全脑组群体素信号强度时间序列模型，应用区域增长(region growing)方法实现脑激活区的提取。

步骤4进一步包括以下几个步骤：

步骤4.1，度量不同体素间的活动相似性。

优选地，使用组群体素信号强度时间序列模型的皮尔逊(Pearson)相关系数作为区域一致性准则，来度量不同体素间的活动相似性。

步骤4.2，使用体素v₀作为种子点，区域增长算法每次选择同种子点活动最相似的邻域体素加入局部一致脑区；迭代进行区域增长算法，直到局部一致脑区达到指定局部区域大小后停止，局部区域内的体素信号强度时间序列动态变化基本一致。

步骤4.3，将上述步骤4.1和步骤4.2应用到全脑体素，根据脑激活区内体素信号强度时间序列动态变化基本一致这一原理，可实现脑激活区提取。

为了说明本发明所述方法的有效性，将该方法应用于针刺机理脑成像研究中。

1.数据采集及实验设计

本测试数据里来自16名被测试者(8男8女)，均为右利手，并符合针刺及fMRI实验要求。实验设计采用传统on-off组块(block)设计，三个刺激(on)组块的时间间隔分别为2分钟和1分钟。穴位选择临床中较为常用的镇痛穴足三里(ST36)，全部实验由专业针灸师操作。

2.数据处理

数据预处理与统计分析均采用统计参数图分析软件(StatisticalParametric Mapping 2，SPM2)。首先，对图像进行时间点校正并采用最小二乘法对功能数据进行头动校正(所有被试的头动平移未超过1mm，旋转角度未超过1°)；其次，对校正后的数据利用12参数的仿射变换将平均图像配准到MNI(Montreal Neurological Institute)标准模板，并将体素重切为3mm×3mm×3mm的大小；采用6mm FWHM(Full Width Half Maximum)高斯核函数对数据进行平滑。最后，用rest软件去除线性漂移，将扫描基线的漂移进行去除。

将本文所述的基于时间序列状态分析的方法应用于预处理后数据，结果如附图所示。图2为分别对应于16个被试的某一体素信号强度时间序列变化模型；图3为该体素组群信号强度时间序列变化模型；图4为应用该方法得到的脑区激活图(XY-轴位观，XZ-冠状位观，YZ-矢状位观)。实验结果说明，本发明所述的基于时间序列状态分析的脑激活区检测方法可有效提取脑激活区。

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

Claims (6)

1.一种基于移动平均时间序列模型的磁共振成像脑激活区检测方法，其特征在于，该方法包括以下步骤：

步骤1，对采集到的原始功能磁共振成像fMRI图像进行预处理；

步骤2，对于每个被测试者逐体素计算每个被测试者的个体体素信号强度时间序列模型；

步骤3，对于组群，计算每个个体在组群中的权重值，并根据全脑个体体素信号强度时间序列模型，逐体素构建组群体素信号强度时间序列模型；

步骤4，根据全脑组群体素信号强度时间序列模型，应用区域增长方法实现脑激活区的提取；

所述步骤2进一步包括以下步骤：

步骤2.1，使用移动平均算法来构建个体体素信号强度时间序列预测统计量；

步骤2.2，基于个体体素信号强度时间序列预测统计量，统计推断个体体素信号强度时间序列变化状态，完成个体体素信号强度时间序列模型的建立；

所述步骤3进一步包括以下步骤：

步骤3.1，构建组群体素信号强度时间序列预测统计量；

步骤3.2，计算每个个体在组群中的权重值D；

步骤3.3，基于组群体素信号强度时间序列预测统计量，统计推断组群体素信号强度时间序列变化状态；

步骤3.4，基于全脑个体体素信号强度时间序列模型、每个个体在组群中的权重值D以及组群体素信号强度时间序列变化状态，逐体素构建组群体素信号强度时间序列模型；

所述步骤4进一步包括以下步骤：

步骤4.1，度量不同体素间的活动相似性；

步骤4.2，使用体素v₀作为种子点，使用区域增长算法选择同种子点活动最相似的邻域体素加入局部一致脑区；迭代进行区域增长算法，直到局部一致脑区达到指定局部区域大小后停止；

步骤4.3，将上述步骤4.1和步骤4.2应用到全脑体素，根据脑激活区内体素信号强度时间序列动态变化基本一致这一原理，实现脑激活区的提取。

2.如权利要求1所述的方法，其特征在于，所述步骤1中对采集到的原始fMRI图像进行预处理进一步包括：对fMRI图像进行时间点校正、头动校正、标准化处理、重新采样、平滑以及去除线性漂移。

3.如权利要求2所述的方法，其特征在于，对fMRI图像进行标准化处理时采用统计参数映射软件自带的平面回波成像(EPI)模板对fMRI图像进行标准化处理。

4.如权利要求1所述的方法，其特征在于，所述步骤2.2进一步采用蒙特-卡罗方法在时间维度上对概率p进行FWER校正来统计推断个体体素信号强度时间序列变化状态，其中，概率p为fMRI时间序列中时间点t所对应的概率值。

5.如权利要求1所述的方法，其特征在于，所述步骤3.2进一步使用期望最大化算法来计算每个个体在组群中的权重值D。

6.如权利要求1所述的方法，其特征在于，所述步骤4.1进一步使用组群体素信号强度时间序列模型的皮尔逊相关系数作为区域一致性准则，度量不同体素间的活动相似性。