CN103006215A - 基于局部平滑回归的脑功能区定位方法 - Google Patents

Abstract

一种基于局部平滑回归的脑功能区定位方法，对数据进行预处理并确定设计矩阵X；以体素v_i为球心、r为半径建立球形选区，提取球形选区中的所有体素的时间序列；根据球形选区内所有体素的时间序列和设计矩阵形成目标函数，并对目标函数进行优化；计算体素v_i的条件特异性效应；转向下一个体素v_i+1，然后重复步骤S2至步骤S4，直到对全脑每一个体素都进行过上述步骤为止；为全脑映射图设定阈值，从而得到和刺激条件相关的脑功能区定位图。基于单体素的回归和基于高斯平滑滤波的广义线性模型都可以视为本发明的特例。本发明可被整合进“探照灯”方法所使用的框架，在求得回归系数之后计算不同预测子系数之间的马氏距离。通过调节超参数α和β，获得不同程度的平滑效果。

Description

基于局部平滑回归的脑功能区定位方法

技术领域

本发明属于图像处理技术领域，具体涉及一种基于局部平滑回归的脑功能区定位方法。

背景技术

功能磁共振成像(functional Magnetic Resonance Imaging，fMRI)以其高时空分辨率，非侵入式等特点在神经疾病诊断治疗和认知神经科学研究等方面得到了广泛应用。fMRI一般指基于血氧水平依赖(bloodoxygen level-dependent，BOLD)的磁共振成像，它通过测量由神经活动引起的脑血流和脑血氧等成分变化而造成的磁共振信号变化来反应脑活动。脑是一个复杂的系统，在受到刺激条件或经历病变时脑的磁共振图像会发生相应的变化。利用脑功能区定位，可以找到某些刺激条件特异性的大脑激活区。

随着fMRI技术的发展，高分辨率fMRI得到越来越广泛的应用。高分辨率成像使得研究人员可以看清精细尺度的神经活动情况，一些以前存在争议的问题也有望得到解决。但是，现存的脑功能区定位方法并不适用于高分辨率成像数据的分析。传统的基于单体素的脑功能映射图的方法，如广义线性模型，依赖于空域高斯平滑滤波。空域平滑滤波会掩盖掉高分辨率成像特有的有价值的细节信息。然而，我们又不能简单摒弃空域平滑滤波。因为平滑滤波在增强功能对比度噪声比方面、提高统计假设的有效性方面均有至关重要的作用。假如不做空域平滑滤波，最终生成的脑映射图类似椒盐噪声，而不是通常见到的团状激活图。在设置一定的阈值之后，产生的是一些零散、细小的激活区，难以和噪声区分开来。另外，高分辨率fMRI是以损失功能对比度噪声比为代价获得高分辨率的。所以，简单地摒弃空域平滑滤波是不可取的。

为解决这一问题，一种被称为“探照灯”的方法出现了。这种方法直接使用在未经平滑滤波的数据上，可以保存数据的细节信息。“探照灯”方法的基本思想是：(1)“探照灯”方法考虑“探照灯”内的所有邻域体素，而传统的基于单体素的广义线性模型一次只考虑单一体素；(2)“探照灯”方法求解一个多变量多元线性回归问题，而传统的基于单体素的广义线性模型求解的是一个单变量多元回归问题；(3)“探照灯”方法采用马氏距离作为衡量不同刺激条件引起大脑活动差异的测度，而传统的基于单体素的广义线性模型采用欧氏距离衡量不同刺激条件引起的大脑活动差异。尽管“探照灯”方法相比较与之前的方法有了很大进步，但是并没有充分利用邻域体素的信息。它使用最小二乘法求解多变量多元线性回归问题，表面上看是同时利用了“探照灯”内的所有体素信息。实际上，这种无规则化的优化问题等价于单体素多元回归问题。也就是说，在没有任何约束的情况下，这种对多个体素进行联合回归的方法等价于对其中每一个体素进行单体素回归。

发明内容

本发明的目的是提供一种基于局部平滑回归的脑功能区定位方法。

为实现上述目的，一种基于局部平滑回归的脑功能区定位方法，包括：

S1对数据进行预处理并确定设计矩阵X；

S2以体素v_i为球心、r为半径建立球形选区，提取球形选区中所有体素的时间序列；

S3根据球形选区内所有体素的时间序列和设计矩阵形成目标函数，并对目标函数进行优化；

S4计算体素v_i的条件特异性效应；

S5转向下一个体素v_i+1，然后重复步骤S2至步骤S4，直到对全脑每一个体素都进行过上述步骤为止；

S6为全脑映射图设定阈值，从而得到和刺激条件相关的脑功能区定位图。

基于单体素的回归和基于高斯平滑滤波的广义线性模型都可以视为本发明的特例。本发明还可以被整合进“探照灯”方法所使用的框架，在求得回归系数之后计算不同预测子系数之间的马氏距离。通过调节超参数α和β，可以获得不同程度的平滑效果，从而实现更加灵活、精确的脑功能区定位。

附图说明

图1是本发明方法的流程图；

图2是采用不同方法进行脑功能区定位的结果对比图，(A)是用广义线性模型处理未经平滑滤波的数据的结果；(B)是用广义线性模型处理经过平滑滤波的数据的结果；(C)是“探照灯”方法的结果；(D)是本发明方法的结果；

图3是以超参数α和β为自变量、聚合度

为因变量绘制的等高线图，等高线刻度显示的是当前聚合度的值，(a)，(b)和(c)分别来自三个典型被试的数据；

图4是单纯使用“探照灯”方法和将本发明方法与“探照灯”方法框架相融合的效果对比图，(a)是在几个经典的感兴趣区域上，两种方法的检测率，(b)是在这几个经典的感兴趣区域上，两种方法探测到的感兴趣区域的体积。

具体实施方式

下面将结合附图详细描述本发明的功能区定位方法，应指出的是，所描述的实施例仅旨在便于对本发明的理解，而对其不起任何限定作用。

图1示出了本发明方法的流程图。

步骤1：数据预处理并确定设计矩阵；

为限制T1效应的影响，我们舍弃每一扫描阶段的前3张扫描图片，然后对留存的扫描图片进行时域校正、空域校正、去除不同扫描阶段的基线差异、进行高通滤波以去除扫描机器漂移和低频伪影。

设计矩阵X∈R^N×P，其中N为扫描的时间点数，P为回归子的个数。P个回归子分别代表对不同刺激条件下大脑响应的预期、头部运动校正参数、以及不同扫描阶段的基线差异。

步骤2：提取球形选区中所有体素的时间序列；

假设连续采集了N幅功能磁共振成像图片，每幅图片包含M个体素，体素v_i的时间序列为y_i∈R^N，i＝1，…，M。以体素v_i为中心、r为半径建立球形选区。体素v_i的时间序列y_i以及邻域体素v_j，j∈NN(i)的时间序列y_j∈R^N被同时用来估计v_i的基重b_i∈R^P，其中NN(i)表示v_i的近邻。

步骤3：根据设计矩阵和球形选区内所有体素的时间序列形成目标函数，并对目标函数进行优化；

其中第一项是常用的针对v_i的广义线性模型，第二项是针对球形选区中其他体素的正则化项。NN(i)表示v_i的邻域体素，y_i表示v_i的时间序列，y_j表示邻域体素v_j的时间序列，b_i表示v_i的基重(回归系数)，ξ_j∈R^P是v_j的容差，f_j＜1是v_j的高斯权重系数，α和β是超参数。目标函数实际上是一个凸二次优化问题，求解该问题可以得到：

$b_i={((\stackrel{\_}{f}\mathrm{\β}+1)X^{T}X-\stackrel{\‾}{f}\mathrm{\β}{X}^T\mathrm{XAX})}^{-1}{X}^T(y_i+\mathrm{\β}\stackrel{\‾}{y}-\mathrm{\βXA}\stackrel{\‾}{y})$

其中，A＝(X^TX+αI)^-1X^T

$\stackrel{\‾}{y}=\underset{\underset{j\≠1}{j\∈\mathrm{NN}\left(i\right)}}{\mathrm{\Σ}}{f}_j{y}_j$

$\stackrel{\‾}{f}=\underset{\underset{j\≠1}{j\∈\mathrm{NN}\left(i\right)}}{\mathrm{\Σ}}{f}_j$

调整超参数α和β可以获得不同程度的平滑效果。当设置β＝0，上述目标函数的结果变成：

b_i＝(X^TX)^-1X^Ty_i

此时，局部平滑回归方法退化为单体素回归方法，回归结果仅和v_i有关，而与邻域体素无关。当设置α＝+∞，β＝1时，目标函数的结果变成：

$b_i={\left((\stackrel{\‾}{f}+1)X^{T}X\right)}^{-1}{X}^T(y_i+\stackrel{\‾}{y})$

此时，局部平滑回归方法等价于先对数据进行高斯平滑滤波，再做广义线性回归。

步骤4：计算体素v_i的条件特异性效应；

假设c∈R^P是对比向量，表示了不同预测子之间的比较关系。那么c^Tb_i作为比较结果，可以用来表征体素v_i的条件特异性效应。假设B∈R^P×M表示全脑体素的回归系数矩阵，则c^TB表示条件特异性脑功能图。

或者采用马氏距离作为条件特异性效应的测度。假设Y_S∈R^N×L是球形选区中的时间序列矩阵，其中L是球形选区内的体素个数。B_S∈R^P×L是球形选区内的回归系数矩阵，那么马氏距离定义如下：

Δ²＝a∑^-1a^T

其中，a＝c^TB_S，∑＝E^TE，E＝Y_S-XB_S。

步骤5：转向下一个体素v_i+1，然后重复步骤2至步骤4，直到对全脑每一个体素都进行过上述步骤为止；

步骤6：为以上得到的全脑映射图设定合适的阈值，从而得到和刺激条件相关的脑功能区定位图；

首先将之前得到的条件特异性效应c^TB或者Δ²转换成z值分布图，也就是使其符合高斯分布。体素v_i在z值分布图中的值用z_i表示。设定阈值z_thresh，假如z_i＞z_thresh得到满足，则v_i是一个超阈值激活体素。所有相邻的激活体素组成了超阈值激活团。假定全脑共有V个激活体素和C个激活团，则聚合度

可以表示为：

$\stackrel{\‾}{V}=\frac{V}{C}.$

运行结果

为了验证本发明的方法，我们用3T的西门子磁共振扫描仪采集了一批高分辨率的功能磁共振成像数据。受试者共有10名(4名男性，6名女性，平均年龄22岁)，实验设计采用的是一个通常用于定位面孔加工脑区的传统实验。实验过程中，8个20秒的面孔图片呈现阶段和8个20秒的物体图片呈现阶段交替进行，受试者需要做的仅仅是被动观看图片。对每一个受试者，我们采集了288张3D回波平面成像(echo-planar-imaging，EPI)全脑图像，每张图像的空间分辨率是2mm×2mm×2.2mm。

依据步骤1中的数据预处理方法对采集到的数据进行预处理。特别注意，用于“探照灯”方法和局部平滑回归方法分析的数据不需要平滑滤波处理，而用于传统的广义线性模型分析的数据须先用6×6×6mm³的高斯核滤波。

采用不同方法进行脑功能区定位的结果如图2所示。可以看到，局部平滑回归方法找到了相对完整的与面孔加工相关的感兴趣区域，并且在感兴趣区域内，面孔、物体引起的激活差异更大。

以超参数α和β为自变量、聚合度

为因变量绘制的等高线图如图3所示。图3(a)、(b)和(c)分别显示了三个典型被试的等高线图。可以看到，针对不同被试，超参数α和β的最优组合不同。而传统的广义线性模型只考虑了α＝+∞，β＝1这一种情况，所以不一定能得到最佳脑功能区定位结果。

局部平滑回归方法是一种灵活的模型，在求得全脑共M个体素对应的回归系数[b₁，b₂，…，b_M]后，既可以直接比较不同刺激条件下回归系数的差值，也可以像“探照灯”方法那样计算它们的马氏距离。假如需要计算马氏距离，应按照步骤4所述执行。图4是单纯使用“探照灯”方法和将本发明方法与“探照灯”方法框架相融合的效果对比图。如图4(a)所示，将本发明方法与“探照灯”方法框架相融合，找到了全部10名受试者的双侧梭状回面孔区。如图4(b)所示，将本发明方法与“探照灯”方法框架相融合，找到的双侧枕叶面孔区体积远远大于单纯使用“探照灯”方法所找到的枕叶面孔区的体积。

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

Claims (7)

1.一种基于局部平滑回归的脑功能区定位方法，包括：

S1对数据进行预处理并确定设计矩阵X；

S2以体素v_i为球心、r为半径建立球形选区，提取球形选区中所有体素的时间序列；

S3根据球形选区内所有体素的时间序列和设计矩阵形成目标函数，并对目标函数进行优化；

S4计算体素v_i的条件特异性效应；

S5转向下一个体素v_i+1，然后重复步骤S2至步骤S4，直到对全脑每一个体素都进行过上述步骤为止；

S6为全脑映射图设定阈值，从而得到和刺激条件相关的脑功能区定位图。

2.如权利要求1所述的方法，其特征在于所述设计矩阵X按下式表示：

X∈R^N×P，其中N为扫描的时间点数，P为回归子的个数，P个回归子分别代表对不同刺激条件下大脑响应的预期、头部运动校正参数、以及不同扫描阶段的基线差异。

3.如权利要求1所述的方法，其特征在于，利用球形选区中邻域体素的时间序列来辅助估计球心体素vi的回归系数，结合设计矩阵X可以得到目标函数如下：

其中第一项是常用的针对v_i的广义线性模型，第二项是针对球形选区中其他体素的正则化项。NN(i)表示v_i的邻域体素，y_i表示v_i的时间序列，y_j表示邻域体素v_j的时间序列，b_i表示v_i的回归系数，ξ_j是v_j的容差，f_j是v_j的高斯权重系数，α和β是超参数。

4.如权利要求3所述的方法，其特征在于，调整超参数α和β获得不同程度的平滑效果。当设置β＝0，上述目标函数的结果变成：

b_i＝(X^TX)^-1X^Ty_i

此时，局部平滑回归方法退化为单体素回归方法，回归结果仅和v_i有关，而与邻域体素无关。当设置α＝+∞，β＝1时，目标函数的结果变成：

$b_i={\left((\stackrel{\‾}{f}+1)X^{T}X\right)}^{-1}{X}^T(y_i+\stackrel{\‾}{y})$

$\stackrel{\‾}{y}=\underset{\underset{j\≠1}{j\∈\mathrm{NN}\left(i\right)}}{\mathrm{\Σ}}{f}_j{y}_j$

$\stackrel{\‾}{f}=\underset{\underset{j\≠1}{j\∈\mathrm{NN}\left(i\right)}}{\mathrm{\Σ}}{f}_j$

此时，局部平滑回归方法等价于先对数据进行高斯平滑滤波，再做广义线性回归。

5.如权利要求3所述的方法，其特征在于，在求得全脑共M个体素对应的回归系数后，直接比较不同刺激条件下回归系数的差值或计算它们的马氏距离。

6.如权利要求5所述的方法，其特征在于，设定阈值z_thresh，假如z_i＞z_thresh得到满足，则v_i是一个超阈值激活体素，其中，z_i表示体素v_i在特异性效应z值分布图中的对应值。

7.如权利要求6所述的方法，其特征在于，所有相邻的激活体素组成了超阈值激活团。假定全脑共有V个激活体素和C个激活团，则聚合度

可以表示为：

$\stackrel{\‾}{V}=\frac{V}{C}.$

Images