0 引言

近年来随着医学影像技术的迅速发展，在磁共振成像技术基础上发展起来的功能磁共振成像(functional Magnetic Resonance Imaging, fMRI)^[1-3]技术，由于其能够无创伤性地对脑功能进行准确的定位, 并且具有较好的可重复性和可行性，因此得到了广泛的关注和研究。fMRI主要是基于核磁共振成像信号的血氧水平依赖性(Blood Oxygenation Level Dependence, BOLD)^[4]，通过改变局部脱氧血红蛋白和氧合血红蛋白的相对含量引起磁共振信号的变化，从而可以有效地检测到大脑皮层不同功能区域的激活信息，获得准确的空间定位，有助于对脑的基础研究和临床治疗^[5-7]。因此，对fMRI数据进行准确的分析是实现脑功能精确定位的基础和前提。

常见的一种分析fMRI时间序列的方法是采用一般线性模型(General Linear Model, GLM)^[8-9]。简单地说，它是将观测信号通过单个或者多个解释变量的线性组合来构成的模型，这种方法给研究人员进行数据分析提供了极大的灵活性和丰富的选择性。然而，即使在最简单的fMRI实验模型中，回归量之间也往往存在一定的关联，即共线性^[10]。例如在一个fMRI实验中，反馈总是紧随刺激2 s之后出现，由于血流动力学反应的模糊效应，刺激和反馈的回归量将会高度相关。这会影响后续模型中回归量的参数估计，从而得到不确定甚至错误的结论。因此，消除模型回归量的共线性成为一个亟待解决的问题。

Andrade等^[11]针对fMRI数据分析中回归模型协变量相关的问题，提出移除协变量之间冗余的部分。具体地说是从其中一个回归量减去与之相关且乘以一个特定系数的回归量，但是最终所得到的参数估计结果是不明确的。之后，Erdeniz等^[12]针对决策变量：奖励预测误差(Reward Prediction Error, RPE)和奖励结果(Reward Outcomes, RO)之间内在的共线性问题，提出采用正交化的方法来解决。利用统计参数图(Statistical Parametric Mapping, SPM)软件自带的正交化方法对RPE和RO信号进行正交化处理，得到它们对应的表示合成BOLD信号的脑区激活。最近，Mumford等^[13]也指出正交化可以解决fMRI模型中的回归量共线性问题，并阐述了正交化可能对一般线性模型回归量参数估计的影响，但没有给出具体的正交化方法，只是通过均值中心化来说明正交化的作用，这对于一些复杂的模型不太适用。此外，上述研究主要是从理论上介绍，缺乏真实的数据实验。

基于上述研究存在的不足，本文具体提出了一种正交化方法，通过将fMRI时间序列的一般线性模型中共线的回归量正交，消除其中一个回归量对另一个回归量在结果变量中的影响。当执行正交化后，共线的部分将分配给特定的回归量，从而提高该回归量的显著性，其结果也更符合真实情况。最后，本文分别采用一些合成数据和当前一个流行的fMRI数据分析软件包——脑功能磁共振图像软件包(Functional magnetic resonance imaging of the brain Software Library, FSL)进行实验，实验结果表明，正交化可以消除一般线性模型的共线性，并提高感兴趣回归量的显著性。

1 基于fMRI数据的一般线性模型

在对fMRI时间序列进行分析时，通常采用一般线性模型，它假设体素k上同一任务的时间序列或不同任务序列的实验数据是一些未知参数β_i的线性组合:

$ {{Y}_{i}}^{k}={{x}_{i1}}{{\beta }_{1}}^{k}+{{x}_{i2}}{{\beta }_{2}}^{k}+\ldots +{{x}_{im}}{{\beta }_{m}}^{k}+{{\varepsilon }_{i}}^{k};\ \ \ i=1, 2, \ldots, n $

(1)

式中：m是未知参数个数；n是实验测量次数；x_ij是与任务或时间有关，但与具体脑区(像素)无关的已知参数，它组成的矩阵X通常又称为设计矩阵；ε_i^k为体素k处的误差，这里假设它们之间相互独立，并服从平均值是0、标准偏差均为σ_k²的正态分布N(0，σ_k²)，写成矩阵的形式，可简写为：

$ {\mathit{\boldsymbol{Y}}^k} = \mathit{\boldsymbol{X}}{\mathit{\boldsymbol{\beta }}^k} + {\mathit{\boldsymbol{\varepsilon }}^k} $

(2)

其中：Y^k是数据组成的列向量；β^k是未知参数组成的列向量；ε^k是误差项组成的列向量。经过这样的转换后，原本对Y^k作统计分析，现在改为拟合出β^k后，得到了许多关于的图像，然后再对它们进行统计分析，脑功能激活图实际上就是根据对参数的统计推断而得到的。一般线性模型主要用来解决两个问题：一是对模型的未知参数进行估计；二是对估计得到的参数作假设检验，以此来判断回归量的显著性。下文对此进行具体的介绍。

1.1 模型参数的估计

建立一般线性模型的方程(2) 后，通常采用最小二乘法来求解参数β的估计值。$ \hat{\beta } $满足等式：

$ {S_e} = \min \left\{ {{{\left( {\mathit{\boldsymbol{Y}} - \mathit{\boldsymbol{X}}\hat \beta } \right)}^{\rm{T}}}\left( {\mathit{\boldsymbol{Y}} - \mathit{\boldsymbol{X}}\hat \beta } \right)} \right\} $

(3)

令$ \frac{\partial {{S}_{e}}}{\partial \hat{\beta }}=0 $，得到关于$ \hat{\beta } $的正规方程：

$ {\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{X}}\hat \beta = {\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{Y}} $

(4)

如果X满秩，那么X^TX可逆，则参数β的估计为：

$ \hat \beta = {\left( {{\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{X}}} \right)^{ - 1}}{\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{Y}} $

(5)

其误差为：

$ {\rm{Var}}\left\{ {\hat \beta } \right\} = {\sigma _k}^2{\left( {{\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{X}}} \right)^{ - 1}} $

(6)

其中：σ_k²是Y_i^k的方差，一般不能直接通过实验测量到，但是可以将其作为方程(1) 中的另一个待求参数，采用最大似然法求出它的估计值为：

$ {\hat \sigma _k}^2 = {S_e}/\left( {n - p} \right) $

(7)

式中：S_e是方程(2) 最小二乘法拟合后得到的残差；n是X的行数；p是X的列数，也即X的秩。

1.2 假设检验

在建立模型并估计出其参数后，需要对参数作统计推断。如果方程(2) 的设计矩阵X满秩，那么参数的估计值符合多维正态分布：$ \hat{\beta } $~N(β，σ²(X^TX)^-1)，它们的线性组合也服从一维正态分布：c^Tβ~N(c^Tβ，σ²c^T(X^TX)^-1c)。根据Fisher定律，$ \hat{\beta } $和$ {{\hat{\sigma }}^{2}} $是相互独立的。因此，模型参数组合量c^Tβ将服从学生分布：

$ \frac{{{\mathit{\boldsymbol{c}}^{\rm{T}}}\hat \beta - {\mathit{\boldsymbol{c}}^{\rm{T}}}\beta }}{{\sqrt {{{\mathit{\hat \sigma }}^2}{\mathit{\boldsymbol{c}}^{\rm{T}}}{{\left( {{\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{X}}} \right)}^{ - 1}}\mathit{\boldsymbol{c}}} }}\tilde{\ }{\mathit{t}_{n - p}} $

(8)

其中，t_n-p是自由度为n-p的t分布。例如，对于假设H：c^Tβ=d，可以通过式(9) 计算t值或p值来检验。

$ t = \frac{{{\mathit{\boldsymbol{c}}^{\rm{T}}}\hat \beta - d}}{{\sqrt {{{\mathit{\hat \sigma }}^2}{\mathit{\boldsymbol{c}}^{\rm{T}}}{{\left( {{\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{X}}} \right)}^{ - 1}}\mathit{\boldsymbol{c}}} }} $

(9)

一般显著性检验的原假设形式为H₀：c^Tβ=0。当检验某个变量的显著性时，首先计算其具体观测值的统计量t值，然后再通过比较选定的显著性水平和计算得到的t值的大小来判断是否显著。若t > t_α(n-p)(其中α为显著性水平)，则拒绝原假设，表明该变量对结果有显著效应；若t≤t_α(n-p)，则接受原假设，此时β=0，表明该变量对结果没有影响。

对于X不满秩的情况，X^TX无法直接求出，此时需要用求伪逆的方法计算，但上述方法仍然适用。

2 共线性的影响

对fMRI时间序列的一般线性模型进行分析时，模型的回归量之间往往存在一定的共线性，这会对后面的参数估计产生不良的影响，甚至导致出现错误的统计推断。例如，如果模型中两个解释变量X₁和X₂具有共线性，那么它们中的一个变量就可以由另外一个变量表征。这时X₁和X₂的参数并不反映各自与结果变量之间的结构关系，而是反映它们对结果变量的共同影响，所以各自的参数已失去了原来的意义，甚至可能会表现出反常的现象，比如估计参数结果本来应该是正的，但实际结果却是负的。其次，对于存在共线性的模型，其参数估计量方差会变大，根据式(8) 计算可知其显著性检验统计量t变小，从而检验接受原假设H₀：c^Tβ=0的可能性增大，这样会使本来影响很大的因素误判为不显著，使模型失去可靠性。此外，由于参数估计量的方差变大，因而对样本值的反应十分敏感，即当样本观测值稍有变化时，模型参数就有很大差异，致使模型难以应用。因此，消除模型共线性对于正确理解fMRI数据参数估计的解释有重要的影响。

3 正交化

为了消除模型中回归量共线性的影响，可以将其中一个回归量与其对应共线的回归量进行正交。为了简化模型，本文主要讨论两个回归量共线的情况。对于两个共线的回归量X′₁和X′₂，如图 1(a)所示，当检验模型中回归量X′₁的显著性时，首先必须有效地消除回归量X′₂对X′₁在整个模型中的影响，即将共线回归量X′₂关于X′₁正交，如图 1(b)所示。

通过计算推导得到规范后正交化的公式为：

$ \mathit{\boldsymbol{X}}{{'}_{2}}^{\bot }=\mathit{\boldsymbol{X}}{{'}_{2}}-\mathit{\boldsymbol{X}}{{'}_{1}}\mathit{\boldsymbol{X}}{{'}_{2}}\mathit{\boldsymbol{X}}{{'}_{1}} $

(10)

经过正交化后，原来两个共线的回归量X′₁和X′₂之间不再相关，即X′₁和X′₂之间相关系数为0，而且此时相应的估计参数也会发生变化。对于回归量X′₂，其参数估计值保持不变，而对于X′₁，由于两个回归量X′₁和X′₂之间共享的部分都已经分配给它，其参数估计将会发生改变。

对于一般线性模型，它具有这样的一个基本性质：只有仅属于该回归量的部分才会影响该回归量的参数估计值。利用维恩图表示上述正交化的过程，如图 2所示。两个共线的回归量X₁和X₂，当检验模型回归量的显著性时，仅属于X₁的部分决定估计参数β₁，仅属于X₂的部分决定估计参数β₂。当把X₂关于X₁正交化后，X₁和X₂之间共享的部分全部分配给X₁，此时等效于处理两个独立的部分X₁和X₂，如图 2(c)所示，对于X₂来说，正交化前后模型仅属于X₂的部分是相同的，因此其参数估计值不变，而对于X₁，由于正交化后X₁和X₂之间共享的部分全部分配给X₁，即属于X₁的部分发生变化，其参数估计值也会随之改变。

4 实验结果与分析

实验环境设置如下：电脑选用i5-4460，3.2 GHz处理器，12 GB内存，Linux发行版CentOS6.5系统，编程语言为Matlab。为了验证本文提出的正交化方法的有效性，分别采用合成数据以及真实数据来进行实验。

4.1 合成数据

合成数据主要通过建模一个刺激过程和主体的反应来验证正交化的作用。由于刺激和响应这两个事件发生得极为相近，因此一般线性模型中刺激和响应的回归量将会高度相关。

图 3所示为一个共线性的例子(纵坐标仅表示数值大小，无量纲单位)，根据第一章假设检验相关理论知识，分别计算出模型的检验统计量t值，如表 1所示。其中带三角形的实线表示刺激，带正方形的实线表示模拟主体刺激2 s后的响应。图 3(a)为创建的合成BOLD信号，它包含刺激和响应相关的效应。将图 3(a)中的刺激和响应分别与血流动力学响应函数进行卷积，得到卷积后的信号曲线，如图 3(b)所示。从图 3(b)中可以看出，刺激和响应之间存在严重的共线性。因此尽管信号与刺激和响应都相关，但根据表 1可知，实际上只有刺激的检验统计量t值显著，这显然是不符合真实情况的。

当对信号中的响应量感兴趣时，就要消除刺激量的共线性影响，如图 3(c)所示，它表示刺激量与响应量正交后的图形。由于只是对刺激量进行正交化，响应量的曲线没有变化，它的解释也不发生改变，但是显然刺激量曲线和响应量不再相关。此外，由表 1可以看出，正交后响应量的统计量t值增大，而刺激量的t值保持不变，究其原因主要是正交后，两个共线回归量之间共享的部分已经全部分配给响应量，其单独表示信号的部分较正交之前增多，而对于刺激量，由于正交化前后仅属于刺激量的部分不变，其统计量t值也保持不变，但此时它的解释发生改变。

图 3(d)表示与图 3(c)相同的情况，仅仅只改变刺激量和响应量正交的顺序，即将响应量关于刺激量正交，此时刺激量的检验统计量t增大，而响应量则不变。此外，由表 1可知，正交化后刺激和响应回归量的显著性t值与单独检验其中一种回归量显著性t值基本相同，说明正交化后的结果更加符合真实情况。

4.2 真实数据

实验数据由Siemens Trio Tim 3.0T磁共振机器扫描获得，采用标准正交头线圈，先进行常规扫描，扫描序列包括T1加权像(T1 Weighted Imaging, T1WI)、T2加权像(T2 Weighted Imaging, T2WI)、磁共振血管造影(Magnetic Resonance Angiography, MRA)、解剖像(对脑组织的容积扫描，层数为160，层厚度l mm)。功能序列采用平面回波扫描序列，基本的扫描参数有：TR(Repetition Time)为重复时间；TE(Echo Time)为回波时间；Slices为层数；FOV(Field Of View)为扫描视野，根据需要扫描的范围确定其值。相应扫描的具体参数取值为：TR=2 000 ms，TE=50 ms，Slices =35，FOV为280 mm×280 mm，这些都是自下而上间隔扫描。功能扫描分为静息、左手指运动、右手指运动三部分完成。首先进行静息扫描，扫描过程中受试者平躺于扫描床上，不做任何任务，闭眼休息。之后为运动扫描，运动任务为左、右手指运动，实验范式采用组块(Block)设计，为静息—左手指运动—右手指运—静息—左手指运动—右手指运动，每个Block持续80 s，一共4个Block。在扫描的过程中，受试者闭上双眼，肘部伸直，双手平放于扫描台上，配合操作者的口令完成静息—运动交替的过程，整个扫描过程中受试者头部保持不动。

FSL^[14-15](http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/)是一种常见的处理fMRI数据的软件，它能够允许用户指定回归量进行正交化。利用FSL对采集的图像进行处理，经过一系列的头动校正、与标准模板配准、标准化、高斯平滑等预处理过程以及统计分析后，最终得到与左、右手指运动对应的脑激活区域。由于左、右手指运动任务彼此之间存在一定的相关性，因此，当检测左手指运动任务的脑激活区域时，将右手指运动任务关于其进行正交化；同理，检测右手指运动的脑激活区域时，则将左手指运动任务与其正交，然后再分别与之前未作正交化的结果进行对比分析。

左、右手指运动任务及分别正交化后的实验结果如图 4~6所示，表 2~4分别为对应图 4~6左、右手指运动任务主要激活区域的定位。

根据实验结果图 4可知，左手指运动任务激活区域主要对应大脑的右半部分，右手指运动激活区域主要对应大脑的左半部分。在实验过程中，左、右手指运动任务之间相互影响，当只对其中一种任务的激活区域感兴趣时，就需要消除另一种任务对它的影响，即通过正交化来消除二者之间的相关性。图 5、6中的黑色圆圈部分分别表示执行正交化后与未做正交化脑激活区域的不同。从表 2~3可以看出，将左手指运动任务关于右手指运动任务正交化后，右手指运动的激活区域与之前相比发生变化，激活体素数增多，而左手指运动的激活体素数保持不变。这是因为经过正交化后，左、右手指运动任务变成两个独立的回归量，原来由左、右手指运动共同表示的激活区域现在全部由右手指运动表示。虽然左手指运动回归量的估计参数保持不变，但它现在表示减去两者共线之后的部分，其解释发生改变。

同样地根据表 2、4可以知道，当检验右手指运动激活的区域时，改变回归量正交化的顺序，此时左手指运动的激活体素数发生变化，而右手指运动的激活体素数保持不变。正交化只是改变其中一个任务回归量的参数估计，而对另一个回归量的参数估计没有影响，但值得注意的是，此时模型中该回归量的参数解释发生改变。

5 结语

共线性的出现会降低基于功能磁共振成像的一般线性模型结果的效力，甚至产生不可靠的参数估计。实验结果和数据分析表明，本文提出的正交化方法能够有效地消除共线性的影响，得到更加准确的结果。通过正交分解共线的回归量，并把它们之间共享的部分分配给感兴趣的回归量，从而提高该回归量的显著性。然而，本文研究中还存在一些需要进一步深入探讨的地方，譬如正交化的顺序会影响显著性检验的结果，因此，对于两个以上的回归量如何采取合适的方法消除共线性等需进一步研究。