fMRI数据分析处理原理及方法

来源： 整理文件的时候翻到的，来源已经找不到了囧感觉写得还是不错，贴在这里保存。

近年来，血氧水平依赖性磁共振脑功能成像(Blood oxygenation level-dependent functional magnetic resonance imaging, BOLD-fMRI)技术得到极快的发展，除了与扫描硬件、扫描技术的进步有关外，更得力于以图形图像等计算机科学为核心的相关学科的支持：图像数据的后处理技术成为fMRI中的关键环节

一、功能图像数据的性质

功能磁共振数据包括解剖（结构）像和功能像两类。解剖像采用高分辨的T1、T2及FSPGR三维成像方式。功能像的处理是fMRI数据处理的关键。因为脑皮层活动瞬息变化，相应要求足够快的成像序列对某一个刺激任务造成的皮层活动进行记录，并且要有对脑血氧代谢的产物——脱氧血红蛋白产生的T2*缩短效应敏感，EPI(Echo planar Imaging)、FLASH(Fast Low Angle Shot)等序列可以满足这两个条件，现在大都采用EPI序列采集fMRI功能像。

EPI于频率编码上采用一系列反向梯度，通过一次激发产生建成一幅MR图像的所有信号，基于小角度激发的GRE-EPI(Gradient echo- Echo planar Imaging)技术，在很短的TR时间内得到一系列（数幅至数十幅）图像。每次采集得到的图像组成一个脑体积(Volume)，相应要求在fMRI实验组块(Epoch/block Paradigm)设计时，每个组块的时间必须为TR时间的整数倍。实际的血流动力相应是一个缓慢的过程，任务激发后信号经过一个小的下降期开始上升，4-8秒达到高峰然后缓慢下降，11-14秒恢复。在事件相关设计(Event-related Paradigm)时，如果不考虑两（次）任务间的相互作用，需要保证间隔时间大于一次响应时间 。但也有研究显示短的刺激间隔时间对统计结果并无多大影响 。（见图1）。

EPI序列以极快的采集速度，在一个数分钟的实验(Session)中，产生数百至数千幅图像，几十个不同时间的脑体积成为EPI图像的时间序列(Time-series Image)。快速以牺牲图像的分辨率为代价，典型的EPI图像采集矩阵为64×64，提高采集矩阵会延长采样时间并且导致更严重的图像几何变形。除此之外，EPI序列图像对外在磁场环境的影响十分敏感，微弱的BOLD信号会伴有大量的干扰成分。较突出的问题有：

1. 扫描过程中的头部运动的影响。虽然可以采取各种物理方法加以限制，但头部的运动还是难以完全消除，其副作用远不止于功能像与结构像叠加融合时的不匹配。头部微小的运动会使激活体素位置改变而造成真正功能信号的改变，场强为1.5 Tesla时，BOLD信号本身只有0.5-2.0％，但通常两个相邻体素的信号差都大于10％，大脑边缘的甚至达到70％。而且头部运动可能是激活相关的规律性运动，将导致激活区的完全错误，严重影响实验的结果。

2. 易感性伪影。由于梯度磁场的高速切换产生的MR设备导体表面强度的涡流，人体头部组织磁敏感性的差别，尤其是靠近副鼻窦等含有空气的空腔时，导致局部磁场不均匀，将使重建的EPI图像在相位编码的方向上产生几何变形，使功能区不准确。

3. 扫描设备和生理运动产生的噪声干扰，多属高频噪声。生理运动包括呼吸、心跳等，特别是这些运动与任务相关时，对BOLD信号的检出影响更大。同时，由于BOLD效应是血流调节，激活区域信号的改变速率有限，生理自发活动会引起热噪声和高时间频率的波动、扫描硬件的不稳可以产生低频漂移。

4. N/2伪影或鬼影(N/2 ghosts)，由于不准确的采集时序和不均匀的静磁场，k空间交替的回波呈献一定的相位差，以方向相反频率读出梯度交替MR信号奇、偶回波的EPI序列，信号经傅立叶变换重建后出现沿相位编码方向的成对假影。是EPI图像质量受损的最大原因。

EPI图像数据的大量、低分辨率及干扰因素严重等特点决定了（1）务必除去与BOLD信号相关的干扰信号，提高信噪比。（2）大量时间序列的四维EPI图像，要通过转化为三维的形式表现出来。（3）通过合适的算法把真正的BOLD信号提取出来。（4）低分辩率的功能像要与高分辨率的解剖像叠加融合，或配准到已知的空间解剖结构中加以表现。是fMRI数据处理和分析的主要任务。可分为数据的处理、分析和结果的呈示（见图2）。

二、功能图像数据的处理

1. 校正(Re-alignment)。
头部运动的校正是一个理想的单体 (Subject)单模态(Modality)配准，常基于刚体运动模型，迭代计算平移、旋转参数，使参考图像（通常为时间序列的第一幅）与后续序列图像之间的不匹配程度最小化，实现所有时间序列图像的配准。三维空间校正选用三个方向的平移与三个坐标轴的旋转6个参数对头部刚体模型进行描述；三维配准时还需要考虑每个脑体积中（TR时间内）头部运动的影响，以二维配准方法分别校正每一幅图像。Friston强调了基于自动回归移动平均模型(Autoregistration Moving Average ,ARMA)的重要性，可以消除对象自旋激励历史中的运动影响。

此外还必需注意在EPI多层采集过程中，同脑体积中每层采集时间的轻微差异（数十毫秒）。在组块设计实验时，由于每个任务组块时间持续时间较长（数秒至数十秒），可以不考虑这些时间差异；但在事件相关设计时，任务激发的时间性要求高，就必需对每层采集的不同时间差异进行校正，保证组成每个脑体积的数十层图像在相同时间内完成。常采用Sinc 法插值。

通常每一个实验采集数百至数千幅图像，大量的数据使校正过程非常耗时，某些机器附带商业软件为了提高处理速度，达到实时效果，而舍弃此步骤。快速运动校正算法的开发对实时成像（Real-time imaging）十分有意义 。

2. 配准(Registration)。低分辨率的EPI功能图像经常需要叠加在高分辨率的解剖图像上进行功能区的辨认，通过配准功能激活映射图和解剖图像实现。因Ghost效应及磁敏感效应导致EPI图像的几何及强度变形，需要对变形的图像进行反卷积(Unwrapping)校正。这是一个单体多模态配准。J. Asbnrner等提出联合头部尺寸和形态的贝叶斯最大后估计量(Bayesian Maximum a posterior estimator, MAP)方法，利用中间图像实现多种类型的功能数据和解剖数据的精确配准。但如果进行数据的空间归一化，这些变形也都可解决。

3. 归一化(Normalize)。将检测的功能激活区准确地映射到高分辩率的解剖结构图上是fMRI可视化的关键，功能激活映射图根本不含任何解剖信息，无法和解剖图配准，但功能映射图和功能图像可共享同样坐标系统，故可以先把功能图像与解剖图像配准，将得到的变换应用于功能映射图与解剖像之间。把空间校正产生的平均图像(Mean image)或配准好的解剖图像与预先设计好的标准解剖空间的模板图像(Template image)的卷积参数应用于每一个断层图像(Slice image)。这样就可以保证不同样本、不同模态的图像数据在相同的坐标系统进行评价。对于单样本分析，可以不归一化到标准空间，而是到单独创建的模板上；对于脑占位或梗塞等脑结构明显受损的样本图像，务必不能归一化到正常的模板上，自动算法的线性和非线性转换过程中会抹除所有受损部位的特有信息，使归一失败，对于这样的样本，除了用单独创建的模板外，还可以采用有偿函数遮盖（Cost-function Masking）技术对病变部位进行遮盖处理，然后再归一化到标准空间中以资比较 。归一化的本质是一个多体多模态配准。

Talairach and Tournoux系统是最经典的标准解剖系统 ，数据来自于实体解剖，Talairach and Tournoux系统和Brodmann's分区之间的对应关系现在已颇为详知，文献资料十分丰富。加大拿McGill 大学Montreal Neurological Institute建立的MNI系统，采用305例正常人的MR脑扫描，经过映射到Talairach and Tournoux获得，如著名的软件SPM99，标准模板即采用MNI系统。MNI系统尚无与Brodmann's分区的对应信息资料，MNI系统脑模较Talairach and Tournoux系统稍大，虽然有的使用者把二者对等使用，但最好采取一定的方法进行坐标点互换 。

由于全局的脑血流改变以及扫描硬件不稳定，时间序列图像的平均图像的平均信号强度随时间发生与功能活动无关的改变，使得每次刺激的响应不在同一水平，减少了统计检测功能激活信息效果。需要调整每一副图像使其平均值等于全局的平均值，即时间序列的归一化。

4. 平滑(Smooth)。

对于硬件不稳及生理运动产生的干扰信号，可以通过平滑消除：空间平滑减小MR图像随机噪声、提高信噪比与功能激活数据的检测能力。通过将fMRI数据与一个三维高斯函数进行卷积积分形成一个滤波器，滤波器的平滑范围可用高斯核(Gaussian kernel)的全宽半高(FWHM)来表示。理论上高斯核应该与反应区的尺度一样，但要保证高斯核一定要大于一个体素的尺度，否则将造成数据再采样，使内在分辩下降。信噪比较低时，采用较宽的滤波器，检测到的激活区覆盖较大的范围。多样本对比的样本间分析时， FWHM也要大一些（8mm），以使各样本数据能够投射到共同的功能解剖像上，减少样本间差异。滤波器虽然可以有效地滤掉特定频率的噪声，也会牺牲一部分频率相当的真正BOLD信号。

对于时间序列信号的低频漂移，可以采用与BOLD信号波形相似的滤波器（FWHM＝2.8mm），对每个体素的时间序列进行时间平滑；如果用短TR采集功能像，可用频带抑制或最小均方适应滤波器去除与呼吸心跳相关的生理噪声。提高反应体素时间过程的信噪比，增加统计检测信号的能力。

此外，虽然真正的BOLD信号主要源于激活脑组织的毛细血管中的血氧代谢的贡献，但由于大血管的流空流入效应，在非激活区也有大量的脱氧血红蛋白流入，造成信号增高，称为“流入性伪影”，出现在较多引流静脉的皮层区域。低场强机器的伪信号更严重，提高场强可以减少这种大血管效应，SE-EPI序列也可以减少流入效应，对于单层EPI成像，通过增加射频翻转脉冲的作用时间可以限制血流敏感性。但多层EPI则无法满足每层足够长的翻转脉冲时间。有学者通过加权各种组织的统计参数图T对比来减少其影响 。

三、功能数据的分析

在数据预处理之后，采用适当的算法把真正的代表激活的象素提取出来，即功能数据的分析。在最先的fMRI研究中 ，仅采用图像相减的简单方法来演示任务依赖的脑区域。这是基于心理学的Pure insertion假说的认知相减(Cognition Subtraction)原理，用任务状态的图像减去控制状态的图像，差值图像高的灰度值反应的就是任务产生的有效活动区。这种方法对判别活动和非活动体素的阈值设置太过草率，并且对运动相关的效应以及其它未知原因干扰特别敏感。更可靠的脑图是采用参数和非参数检验的方法。
常用的有零假设t检验，基于每个体素计算，加权平均信号差异，t值大于设定的阈值（如p=0.05）的体素认为是激活，常以伪彩的形式表现出来。相关系数法脑图中，每个体素都与线性交叉相关系数r值有关，此系数表示时间序列的体素信号强度与参考函数的相关性，测定的是时间序列过程中体素的灰阶值与期望的氧代谢反应间的关系，相关性大于设定阈值的体素认为是激活。此外还有F检验，z检验等。从统计的观点来看，这些参数检验可以认为是广义线性模型(General Linear Model ,GLM)的特例。GLM是由K.J. Friston和其同事用作PET数据处理时开发的一个标准的统计工具，可以将所有感兴趣和非感兴趣的因素都包含于设计矩阵，如果能够充分考虑时间序列间的时间空间自相关，可以用于fMRI数据的分析。

假设检验时，首先构建关于某一统计量的统计参数映射图，计算每个体素反应的时间过程与参考函数之间的线性相关系，根据检验的显著性水平确定一个阈值，对零假设进行检验，通过阈值化统计参数映射图判别激活与非激活。构建统计映射参数图时，重要的是每个体素的灰度水平时间过程与期望的血流动力函数相似程度。通常选用血流动力相应的脉冲函数与一个理想的on-off函数的卷积积分作为参考函数，所以准确地建立血流动力相应模型十分重要。目前已提出多种建立血流动力相应模型的方法。如Bandettini的傅立叶频谱分析技术、Bullmore的正、余弦波的线性组合拟合实验数据等。

上述方法最主要的问题是如何选择阈值分隔统计参数映射图，来确定激活与非激活体素。确定合适的未校正的单象素显著性阈值非常困难，低的阈值可以增加激活检出的敏感性，但将非激活区作为激活区的可能性增大，增加检测结果的假阳性率。并且GLM方法对每个体素进行假设检验，变成多假设(Multiple Comparisons)检验，总体脑体素的检验将导致更多错误率。为了控制假阳性，常用Boferroni法校正，但过于保守，导致检出率下降。采用高斯随机理论(Gaussian Random Field)，可以保证体素以上水平多假设时假阳性的发生，需要采用相同高斯核对图像进行平滑，以保证数据逼进高斯分布。对这样的数据，K.J. Friston提出检验的等级理论，可分为集合(Set)水平、聚类(Cluster)水平和体素(Voxel)水平，虽然增强了统计能力，但降低了空间分布特性。也可认为真正激活的体素相邻聚类的超过阈值的可能性也较大，用联合强度阈值与聚类尺寸阈值分隔统计参数映射图法，在降低假阳性发生的同时又不降低统计能力。蒙特卡罗仿真技术不需很多假设，但较耗时。

前面提到的体素依赖方法只适用于时间参数已明确知道的任务设计的实验数据分析，对于未知刺激任务时间的实验，如睡眠、癫痫放电等自发生理活动的数据分析时，将无法应用。这类实验的数据可采用主成分分析(Principal Component Analysis, PCA)和独立成分分析(Independent Component Analysis, ICA)等多变量分析，将fMRI数据分解成正交的空间成分或具有不同时间过程的独立的成分，提取包含于时间序列图像中的功能信息，不需要任何血流动力学响应的时间过程数据及皮层幅度的先验假设，其实验设计也就无需依赖任何实验模型（如组块或事件相关）。 故体素依赖的单变量方法又称模式驱动(paradigm driven)，相应多变量分析称数据驱动(data driven)分析模式。PCA通过检测随实验条件变化的开始一部分空间特征模式的时间形式，确定与反应有关的功能系统的分布特征，侧重于描述功能系统的分布而不是定位，用于探索各功能区之间的相互联系。ICA通过提取一系列空间独立的空间模型，相比PCA更侧重空间定位，最适合于探索一个新假说模型的发生而非已知假设的检验。如fMRI对药物作用、睡眠、饥饿感的中枢机制研究等 ，近来有把时间聚类分析(Temporal Clustering Analysis) 用于无EEG联合的癫痫灶定位研究中 。PCA和ICA的缺点是对于大部分的不同成分的数据相关性难以给出一个生理解释。

四、功能磁共振数据可视化方法

fMRI数据经过处理和分析，以直观的形式表现出来，以方便结果观察和引用。除了解剖像与映射参数图叠加外，还可采用大脑皮层重建，提供关于大脑皮层表面解剖结构和几何特性，依此对反应的功能区进行皮层定位。对标准T1解剖像进行灰白质及脑脊液成分分隔(Segment)，行皮层解剖重建。在功能磁共振的视网膜脑图(Retinotopic map)技术中 ，把枕叶的脑沟脑回结构展开，在平面图像上进行评价。有许多重建方法，如基于体素方、2D轮廓重建方法等，重要的是保持皮层的解剖拓扑结构。

现在普遍采用Brodmann's 分区对脑功能区定位，由于脑皮层结构的特异性，除了初级运动及感觉皮层区域较恒定外，其它功能区与解剖关系之间变异普遍存在，除非对神经解剖以及Talairach and Tournoux系统很熟悉，一般都难以在肉眼下对反应区进行功能定位。如前所述，把个体脑图归一化标准脑结构之后，就可以方便地对反应区坐标点按Brodmann's分区进行确认，也有专业的软件自动处理 。

以上简单介绍了fMRI数据处理与分析的原理及方法。这些步骤的实现均靠软件根据不同算法完成。专业软件多种多样，但方法和步骤都基本相同。国际较为通用的功能影像软件有综合的处理分析软件，如英国伦敦大学神经影像科学系Wellcome实验室的SPM(Statistical Parametric Mapping)系列软件以及MCW AFNI(Medical College of Wisconsin Analysis of Functional Neuroimages)软件等。以及专项功能处理的软件，如图像浏览、格式转换的MRIcro软件

(http://www.psychology.nottingham.ac.uk/staff/cr1/mricro.html)、运动校正的AIR软件(http://bishopw.loni.ucla.edu/AIR3/)等。很多是开放的免费软件，可以在相关的网站上查询下载。