CN101287410A - 脑功能分析方法及脑功能分析程序 - Google Patents

Abstract

本发明提供一种方法包括：获取脑功能数据和弥散张量数据(S10)，基于该弥散张量数据计算相邻体素间的联系度(S30)，基于该脑功能数据和相邻体素间联系度来形成数据估计值(S40)，对该数据估计值进行非参数的回归分析(S50)，基于分析结果形成和显示图像(S60，S70)。

Description

脑功能分析方法及脑功能分析程序

技术领域

本发明涉及一脑功能分析方法及脑功能分析程序，该脑功能分析方法及脑功能分析程序使用不同检测技术获得的数据分析脑功能，如：fMRI(BOLD成像)，PET(正电子发射断层成像)或类似技术。

背景技术

如今典型的非介入性脑功能检测方法有fMRI，PET，MEG(脑磁图描技术)及其类似技术，其中fMRI被认为具有很强的空间分辨能力，并被广泛应用。

fMRI是一种能通过反映脑不同物理量来确定脑激活区域的非常有效的脑功能检测手段(参见非专利文献1，2)。它除了能和检测脑结构的MRI解剖图像一样反映活体器官的质子密度，纵向弛豫时间T₁和横向弛豫时间T2以外，还能特别的检测出脑激活区域中血液流量的增加。现在已知，在脑激活区域的血液流量会增加，同时血液中血红蛋白的磁性也会因其带氧状态即含氧血红蛋白(氧合血红蛋白)和脱氧血红蛋白(去氧血红蛋白)而不同。人们认为fMRI在激活区域的信号(BOLD信号)会增加，因为在增加的动脉血液中，具有干扰磁场能力的去氧血红蛋白的量减少了。因此，我们可以使用fMRI，以脑工作时BOLD信号的变化为线索，确定进行该工作的脑区域(激活区域)。

fMRI检测得到的BOLD信号的时间序列数据的分析技术包括基于一种广泛的线性模型的SPM(统计参数绘图)(参见非专利文献1)，或是基本主要组份分析或独立组份分析的数据分析(参见非专利文献3)，或类似分析。这些分析技术的特点是其输出的结果中的BOLD信号的时间序列数据的每个三维像素(体素)都各自经过统计处理，最后形成图像从而确定脑的激活区域。

然而在上述所说的分析技术中存在着以下问题：在分析数据时，神经网络结构并未被考虑。脑中许多神经细胞通过神经突触构成复杂的网络，而根据近些年来的脑研究，我们认为脑中每个神经区域通过互相协作形成上述神经网络，进而作为一个整体发挥出更强的功能。例如：在脑完成一项工作时，其中的很多区域都出现激活现象。虽然上述所说的分析技术可以用于分析该激活现象，确定激活区域，但是却很难确定激活区域间的联系。

造成上述现象的原因有多个。因为fMRI是如上所述依据血液流量检测血液信号的，fMRI可以捕捉到具有较高血液流量的脑灰质的活动(神经细胞体)，但却很难捕捉到血液流量较低的脑白质的活动(神经细胞轴突或神经纤维)。

同时为了捕捉作为神经网络结构基础的神经纤维组的走向，一种DTI(弥散张量成像)方法在近些年来越来越受到人们的关注。该方法使用MRI检测身体器官内质子的弥散度，作为一种新的参数(参见非专利文献4)。当时使用常规的解剖性MRI时，神经纤维在T₁加权图像中是强信号，在T2加权图像中是弱信号。神经纤维在T₁加权图象中变成强信号是因为髓磷脂的存在。髓磷脂是由双分子层类脂物和大型蛋白质组成，并沿着神经纤维走向成形。因此才出现了各向异性，即沿着神经纤维走向上的质子的弥散常数很大，而其在垂直方向上的弥散常数很小。DTI是使用MPG(弥散梯度磁场，)检测弥散的各向异性，从而突出质子弥散情况的技术。该技术可以应用于分析ST(Stejskal-Tanner)脉冲序列检测获得的BOLD信号的时间序列数据的强度S′(l，m，k，i)，在SE(回波)脉冲的前后，STG(Stejskal-Tanner Gradient梯度)脉冲被加入用于检测弥散。

其中：l，m，k为正离散变量，代表三维体素的位置，它们分别代表体素在X，Y，Z方向上的值。另外，i为正离散变量，代表检测时间。BOLD信号强度S′(l，m，k，i)表达式如下：

$S^{\′}(l,m,k,i)={\mathrm{\ρ}}^{\′}(l,m,k,i)\exp \left(\text{-b}{\stackrel{\→}{G}}^{T}D(l,m,k)\stackrel{\→}{G}\right).\·\·\·\left(1\right)$

当弥散加权图像形成后，式1中的弥散张量将成为分析数据。该弥散张量如下所示：

$D(l,m,k)=\left(\begin{array}{ccc}{D}_{\mathrm{ll}}(l,m,k)& D_{\mathrm{lm}}(l,m,k)& D_{\mathrm{lk}}(l,m,k)\\ D_{\mathrm{ml}}(l,m,k)& D_{\mathrm{mm}}(l,m,k)& D_{\mathrm{mk}}(l,m,k)\\ D_{\mathrm{kl}}(l,m,k)& D_{\mathrm{km}}(l,m,k)& D_{\mathrm{kk}}(l,m,k)\end{array}\right)\·\·\·\left(2\right)$

其中：ρ′(l，m，k，i)代表当没有使用MPG时，BOLD信号强度(普通脑功能数据分析对象)，而b为MPG强度参数。请注意ρ′(l，m，k，i)如下式所示：

${\mathrm{\ρ}}^{\′}(l,m,k,i)\∝f\left(v\right)\·{\mathrm{\ξ}}^{\′}(l,m,k,i)\·(1-\exp (-\frac{T_R}{T_1}))\·\exp (-\frac{T_E}{T_2})\·\·\·\left(3\right)$

其中：f(v)代表流动速度；TR为反复时间；TE为回波时间；ξ′(l，m，k，i)为质子密度。

近些年来，在分析fMRI测得的BOLD信号ρ′(l，m，k，i)的时间序列数据后再使用弥散张量数据D(l，m，k)将目标脑区域连接起来的脑分析方法已经被广泛使用(参见非专利文献5)。

非专利文献1：“Human Brain Function：2nd-Ed.″，Richard S.J.Frackowiak，et al，ELSEVIER ACADEMIC PRESS，2004”

非专利文献2：“Image of Mind″，M.I.Posner and M.E.Raichle，WH Freeman&Co，1997”。

非专利文献3：“Independent Component Analysis：Theory andApplications″，T.W.Lee，Kluwer Acadmic，1988”。

非专利文献4：“Korede wakaru diffusion MRI″S.Aoki，O.Abe，Syuujyun sha，2002”。

非专利文献5：“Combined functional MRI and tractography todemonstrate the connectivity of the human primary motor cortex invivo″，Guye M，et al.，Neuroimage，Vol.19，pp.1349-1360，2003”。

发明内容

在上述的脑功能分析方法中，弥散张量D(l，m，k)并不是用于分析BOLD信号ρ′(l，m，k，i)时间序列数据本身的。因此其有效性值得怀疑。

本发明就上述问题提出了相应的解决方案，还提供了基于非介入检测方法(如：fMRI，PET或类似方法)测得的脑功能数据和能够确定神经纤维组走向的弥散张量上，分析脑激活区域间连接结构的脑功能分析方法和脑功能分析程序。

一方面，本发明提供了一种脑功能分析方法，该方法包括：一、在从体素到体素基础上，获得能够确定脑激活区域的脑功能数据，和能够确定脑中质子弥散程度的弥散张量。二、在弥散张量基础上形成相邻体素间联系度估计值。

另一方面，本发明还提供了一种脑功能分析程序，该程序使计算机获得可以确定脑激活区域的从体素到体素的脑功能数据；能够确定脑中质子弥散度的从体素到体素的弥散张量；在弥散张量基础上形成相邻体素间联系度估计值，并在该基础上分析脑数据。

附图说明

图1是本发明实施例1中脑功能分析仪器的构型示意图。

图2是典型的fMRI检测方法示意图。

图3是典型的fMRI图(二维切片图)，其中图3A是位于床上的头，图3B为组成头部二维横截面(二维切片图)的一个体素。

图4是二维切片图中脑功能信息体素所对应的时间序列数据和弥散张量表。

图5是图4中脑功能信息时间序列数据的一个例子，其中图5A为检测得到的实际数值，图5B为基于临界值，通过组合两个实际检测值获得的数值。

图6是基于图4所示脑功能信息时间序列数据和弥散张量数据所得立体构型，其中图6A是位于床上的头，图6B是由图4脑功能信息中的时间序列数据和弥散张量数据生成k-th二维切片图像Sk的步骤。

图7是由图6二维切片图像形成的三维头部图像，其中图7A是位于床上的头，图7B是通过收集二维切片图像Sk形成的三维头部图像。

图8是通过使用图1中脑功能分析仪器分析数据的步骤流程图。

图9显示了本发明数据预处理技术，其中图9A显示了在与一体素相邻的两个体素均被认为是灰质的区域中的数据预处理；图9B显示了在与一体素相邻的两个体素中只有一个被认为是灰质的区域中的数据预处理；图9C显示了在与一体素相邻的两个体素均不被认为是灰质的区域中的数据预处理。

图10是本发明实施例2中脑功能分析仪器的构型示意图。

图11是通过使用图10中脑功能分析仪器分析数据的步骤流程图。

图12是通过使用图10中数据平滑方法实现的平滑技术。

图13为实施例2修饰例中脑功能分析仪器构型示意图。

图14为使用图13中脑功能分析仪器进行数据分析的步骤流程图。

图15为本发明实施例3中脑功能分析仪器构型示意图。

图16为使用图15中脑功能分析仪器进行数据分析的步骤流程图。

图17为使用图15的聚集装置300的聚集技术，其中，图17A显示了聚集处理前联系度矢量值，图17B显示了聚集处理后联系度矢量值。

图18为实施例3修饰例中脑功能分析仪器构型示意图。

图19为使用图18中脑功能分析仪器进行数据分析的步骤流程图。

图20为本发明实施例4中脑功能分析仪器构型示意图。

图21为使用图20中脑功能分析仪器进行数据分析的步骤流程图。

图22为实施例4修饰例中脑功能分析仪器构型示意图。

图23为使用图22中脑功能分析仪器进行数据分析的步骤流程图。

图24显示了脑在声音刺激下进行简单的重复时，通过不同的脑功能分析方法所得到的结果，其中，图24A和24B显示了使用SPM的T-检测的结果(图24A中T-检测的临界值被修正，而图24B中T-检测的临界值并未被修正)。图24C为组合使用SPM和跟踪成像法的脑功能分析法所得到的分析结果。图24D为使用本发明实施例1中的脑功能分析法所得到的分析结果。

具体实施方式

本发明的特征在于通过非介入性检测设备，如MRI或类似设备在从体素到体素基础上获取脑功能数据及弥散张量数据，进而确定脑激活区域位置，根据弥散张量数据形成相邻体素间的联系度估计值，同时利用该估计值分析脑功能数据。在不同分析技术中，如何考虑上述所说的脑功能数据中相邻体素的联系度估计值，及如何使用该估计值可以有所不同。

因此，以下实施例将结合附图对本发明进行详细描述。

实施例一

图1为本发明实施例一中脑功能分析仪器构造示意图。本发明在实施例一中提出了一种脑功能分析仪器10，其带有以下配置：脑功能数据获取装置1，用于获得MRI设备50测得的脑功能数据中原始时间序列数据ρ′(l，m，k，i)(如：等式(3))；弥散张量数据获取装置2，用于获得MRI设备50测得的弥散张量数据D(l，m，k)(如：等式(2))；数据预处理装置3，用于预处理脑功能数据获取装置1获取的脑功能信息原始时间序列数据ρ′(l，m，k，i)；体素间联系度计算装置4，用于根据弥散张量数据获取法2获得的弥散张量数据D(l，m，k)来计算代表相邻体素联系度的联系度矢量(如：

)；数据估计值形成装置5，用于从数据预处理装置3预处理后的脑功能信息时间序列数据ρ′(l，m，k，i)和通过体素间联系度计算装置4得到的联系度矢量

中得出预定估计值Q；数据分析装置6，用于计算数据估计值Q的极值(在本实施例中为最小值)；图像生成装置7，用于生成数据分析装置6计算得到的数据图像；图像显示装置8，用于显示图像生成装置7生成的图像；存储装置9，其带有次储存，用于记录上述装置2到装置8中获得、计算、形成、分析和生成的各种数据和主存储，用于储存电脑可读的脑功能分析程序，用于执行脑功能分析的每一步骤(此后将详细描述)。CPU(中央处理器)根据存储装置9中储存的脑功能分析程序控制上述装置2至装置9的每一步骤，或是类似过程。以上所用符号的意义在以下内容中将描述。

在此，MRI设备50是由以下部件组成：一磁性组件，该组件带有能够产生静磁场，进而引起核磁共振的静磁场线圈，和能够发出高频共振波从而检测共振信号的高频线圈，和能够产生梯度磁场，进而在共振信号中加入位置信息编码的梯度磁场线圈，或类似部件，和用于控制激活这些线圈的系统，或类似系统，其中，根据操作人员的要求，MRI设备50产生不同的fMRI数据(与血液流量相关的脑功能信息的原始时间序列数据ρ′(l，m，k，i)，与神经纤维走向相关的弥散张量D(l，m，k)，或类似数据)，并将这些数据传输给脑功能分析仪器10。请注意，本发明并非一定需要通过MRI设备50获得脑功能信息的原始时间序列数据ρ′(l，m，k，i)，也可以通过其它非介入性检测设备获得该数据，如PET，通过它也可以获得多种有关脑功能的数据。

同时，对于图像显示装置8，有多种显示设备可用，例如：阴极射线管(CRT)显示器，TFT液晶显示器，等离子显示器，或类似设备，或者是各种打印机，例如：喷墨打印机，激光打印机，或类似设备。另外，存储装置9由，例如：RAM(随机存取存储器)，ROM(只读内存)或类似设备组成。另外，存储装置9的主存储和次存储是相互独立分开的，因此，主存储中的内容可以存放在例如光盘中，如硬盘，floppy软盘(注册商标)，或是只读CD-ROMs，磁带，存储芯片或者类似设备。

本实施例中，虽然脑功能分析仪器10使用的显示装置包括的主要形成装置10A，图像显示装置8和存储装置9是组合在一起的，但是图像显示装置8，或者图像显示装置8和存储装置9，也可以与主要形成装置10A分开来用，作为独立的图像显示单元和存储单元。在任何组合方法中，脑功能分析仪器10都是由电脑来完成。另外，CPU 100从存储装置9中读取脑功能分析程序，控制上述装置2到装置9中的每一个步骤。

这里所说的电脑是指一种能够根据预先设定的标准处理具有一定结构的输入信息，然后构建和输出处理结果，例如普通的电脑，超级电脑，主机，工作站，微型电脑，服务器，或类似设备都包括在其中。另外，也可以使用一种系统(例如分布式计算机系统)，其中包括两台或更多台电脑通过网络连接(如：内部网，局域网，广域网)，或由其组合组成的联系网络相互连接。

另外，为了增强对本发明的理解，图2至图7对本发明的脑功能分析流程进行了概括性的描述。图2显示了典型的fMRI检测方法。图3显示了典型的fMRI图像(二维切片图像)，图3A为位于床上的头，图3B显示了组成头部二维横截面(二维图像)的体素。

例如，在某个实施例中，我们让研究对象完成一项任务(如：让手指轻轻拍打)，与这项任务相关的脑区域位置通过fMRI检测被确定。在图2所示的例子中，一套程序在一次测量中被重复三次，该程序是由一段时间T内的任务和相同长度时间T内的休息组成。这里，水平线为时间轴，其中“ON”表示任务时间，“OFF”表示休息时间。通常在一个实验中，fMRI进行几十次检测。虽然在本实验中，fMRI进行了24次检测，但是任务开始和结束时的检测图像一般不被使用。这是考虑到检测时间的滞后性。所以，这就意味着有效的fMRI检测实际只进行了18次(相应的线为18条粗线)。ti(i＝1，2，...，24)表示各个实际检测时间。以下内容中，ti将只以i表示(i＝1，2，...，I)(分离的正整数)。

图3B是通过fMRI检测获得的二维切片图像Sk。这里的K表示实验中获得的二维切片图像数量。在常规fMRI检测中，我们一般可以获得20至30张二维切片图像(图片中，K＝20至30)。虽然对图3B中的二维切片图像Sk(k＝1，2...K)只是做了表面性的描述，但是它们却是由三维立体像素，又称体素(voxel)所组成的。虽然由64×64体素组成的二维切片图像作为一个整体显示在一张图片中，但是该二维切片图像也可以被分为不同数量体素。虽然二维切片图像Sk中的任何一个体素都可以通过两个分离的变量l和m表达为(l，m)，但是它也可以通过在预定规则下对其设定一个数进而只用一分离正变量j表示。从数学的角度，我们可以在预定规则下得到j和(l，m)一一对应关系(对于(Sk)来说，j≡(l，m))。在本实例中，由于l的上限L和m的上限M都是64(L＝64，M＝64)，所以导致了j＝1，2，...，4096。在以下内容中，这些表示方法都将被适当地使用。另外，图3B中的圆圈表示头部横截面的轮廓。

图4表示了通过脑功能数据获取装置1获得的，并经过数据预处理装置3预处理后的(将在后面描述)构成二维切片图像Sk的体素j的脑功能信息时间序列数据ρ′(j，k，i)(≡ρk(j，i))和通过弥撒张量数据获取装置2得到的D(j，k)(≡Dk(j))。如背景技术所述，当弥散加权图像通过DTI方法生成后，也可同时从ST脉冲序列图像中获得二维切片图像Sk的脑功能信息原始时间序列数据ρ′(l，m，k，i)(≡ρ′k(l，m，i))和弥散张量D(l，m，k)(≡Dk(l，m))。该ST脉冲序列图像是在SE脉冲前后加入STG脉冲用于检测弥散而获得的。ρk(j，i)代表对原始时间序列数据ρ′k(l，m，i)进行先设预定预处理后(将在后面描述)获得的脑功能信息时间序列数据。就是说，图中的ρk(j，i)代表二维切片图像Sk中体素j在某一时间点i上的脑功能信息。同时Dk(j)代表二维切片图像Sk中某一体素j的弥散张量信息。ρk(j，i)特指数量，Dk(j)在矩阵中由等式(2)代表。他们可以是如等式(4)中所示的值，如：

${\mathrm{\ρ}}_k(j,i)=5000,D_k\left(j\right)=\left(\begin{array}{ccc}1.0& 1.2& 1.1\\ 1.2& 0.8& 0.9\\ 1.1& 0.9& 0.7\end{array}\right).\·\·\·\left(4\right)$

由于弥散张量是如本实例所示的典型的对称张量(非对角线元素具有对称性)，其大体拥有六个部分。另外，图中的CLASS表示该任务是否被完成，“ON(＝1)”和“OFF(＝0)”分别代表“ON”和“OFF”，如图2所示。

图5是图4表示的脑功能信息时间序列数据ρk(j，i)的一个实例的表格。虽然图5(A)中脑功能信息ρk(j，i)的实际时间序列数据是连续数值，但是一旦“4000”被设定为激活临界值之后，该数据可在预处理时被进一步二元化，即体素值大于该临界值的处于激活“A”状态，而小于该临界值的处于非激活“I”状态，如图5(B)所示。

图6为图4中弥散张量Dk(j)和脑功能信息时间序列数据ρk(j，i)的立体构型图。其中图6A显示了置于床上的研究对象的头部。图6B显示了通过图4脑功能信息时间序列数据和弥散张量生成k-th二维切片图像Sk的步骤。图6B中体素j1和j2上的五边形标记表示每个体素内部空间(即体素的梯度场和张量场)，它们对应图4的脑信息时间序列数据和弥散张量。

例如：在背景技术所述的SPM中，可进行如下预处理：

(a)重排：下面的图像与上面的二维切片图像排列一致(如，S1)。这样可以补偿检测过程中头部移动带来的位置偏差，去除相应的错误信号。

(b)空间标准化：为了对多个研究对象收集数据并进行比较，我们将每个研究对象的数据调整为Talairach标准脑。

(c)平滑化：Gaussian型过滤器被用于过滤原始带噪音的时间序列数据。这样可以在不降低空间分辨率的情况下提高分析的灵敏度。

然后，对于通过MRI设备50获得的脑功能信息原始时间序列数据ρ′k(l，m，i)，我们使用多种不同的检测方法检测其中的每一个体素。而对于由多个研究对象组成的不同小组，我们在各组间进一步通过对象的个人数据检测其与任务相联系的脑激活区域。

另外，在SPM中，我们事先建立了一个激活预测模型，该模型与经过上述方法预处理后的功能数据的时间序列数据ρk(j，i)(＝ρk(l，m，i))的相配程度需要看不使用弥散张量Dk(l，m)的一般线性模型所得的估计参数。

同时，在本发明中，除了上述所说的预处理以外，我们还对使用通过弥散张量获取装置2获得的弥散张量Dk(j)(＝Dk(l，m))对表达脑功能信息的原始时间序列数据ρ′k(l，m，i)进行(d)型预处理。该预处理方法将在以下内容中详细描述。随后，对使用弥散张量Dk(j)(＝Dk(l，m))预处理的表达脑功能信息的原始时间序列数据ρ′k(l，m，i)进行数据处理，该处理方法将在以下内容中详细描述。由此产生某图像显示目标的二维图像ak(j)(＝ak(l，m))。

$\mathrm{CV}=\frac{1}{I}\underset{i=1}{\overset{I}{\mathrm{\Σ}}}{\{\mathrm{\φ}\left(i\right)-{\widehat{\mathrm{\φ}}}_{\left[i\right]}\}}^2.\·\·\·\left(28\right)$

除了i-th值外，

是通过计算回归系数得到的值，由此来预测使用回归系数得到的i-th的检测值。

接下来，在步骤S60中，图像生成装置7生成图片，其中拥有正回归系数(即正关联关系)的体素被设置为偏向白色，拥有负回归系数(即负关联关系)的体素被设置为偏向黑色。例如，基于上述计算得到的回归系数

设定回归系数为0(即无关联)的体素为灰色度标准。

另外，可以显示体素值频率在直方图中为前5％或更高的值为灰色图像，或者可以根据预设规则生成全色图像。

随后，在步骤S70中，图像显示装置8显示上述立体图像。

如上所述，在实施例一中，除了对根据MRI设备测得的反映脑功能信息时间序列数据所作的误差评估以外，估计值也包括基于MRI设备测得的弥散张量数据对空间方向上连续性的评估。估计值也进行了非参数线性回归分析，所得到的回归系数是用于图像的显示，从而使脑激活区域的位置可以被确定，同时还考虑到脑激活区域间的连接结构。

虽然在S20中进行了(d)型预处理，而不是本实验中普通的(a)到(c)型预处理，但是(d)型预处理也可以不用，这在以下实施例也是同样的。

虽然在实施例一中，等式(23)用于计算数据估计值，但是我们可以根据本发明目的采用的分析方法计算各种不同的估计值。因此，在以下几个对实施例一的修饰中，为了避免重复，我们只是描述其与实施例一不同的部分。

<实施例一的修饰例1>

在实施例一中，数据分析装置6在S50步中使用线性回归方程(等式(19))进行非参数型回归分析。但是，更严格的说，回归表达可以通过一般n维方程式进行，如：二次方程式，三次方程式，或类似方程式。在进行这样的非线性回归分析时，我们可以使用神经网络模型。典

值得注意的是，为了简化的目的，我们以二维切片图像作为一个例子，但实际上图7A和7B所示的头部立体图形数据ak(j)(＝a(j，k)＝a(n))是三维数据分析。使用弥散张量数据Dk(j)的数据处理也被用于那个实例中，以下将作详细介绍。这里的“n”是代表体素位置的独立正数变量，在预定规则下，立体构形中的所有体素都有一个相应的指定值。虽然代表二维切片图形的在X，Y坐标上的位置变量j和代表第三维方向(z方向)的变量k在图片中被分别说明，但是在数学上，它等同于上述所说的单一变量“n”(n≡(j，k))。请注意，使用单一变量“n”和在X，Y，Z轴上使用三个变量(j，l，k)来代表同一个三维体素在以下内容中都将被应用到(n≡(l，m，k)≡(j，k))。

以下我们将具体描述在以上假定条件下使用图1所示脑功能分析仪器10进行的数据分析方法。图8是使用图1所示脑功能分析仪器10分析数据的步骤流程图。

在步骤S10中，通过脑功能数据获取装置1获得MRI设备50检测得到的表达脑功能信息的原始时间序列数据ρ′(l，m，k，i)，同时通过弥散张量信息获取装置2获得MRI设备50检测得到的弥散张量D(l，m，k)。

下一步，通过数据预处理装置3进行各种预处理，例如位置补偿，去除噪音，或如类似上述SPM(a)到(c)的预处理方法，用于处理在上述S20中获得的表达脑信息的原始时间序列数据ρ′(l，m，k，i)。

在以上同一步骤中，为了对上述常规预处理后的表达脑功能信息时间序列数据ρ′(l，m，k，i)作进一步(d)型预处理(数据转化为神经信号)，可使用数据预处理装置3进行计算，例如使用以上获得的弥散张量D(l，m，k)在符合弥散张量D(l，m，k)的本征等式(5)的本征矢量(

α＝1，2，3)中计算最大本征值λM的相应值()。

$D(l,m,k){\stackrel{\&RightArrow;}{v}}_{\mathrm{\&alpha;}}(l,m,k)={\mathrm{\&lambda;}}_{\mathrm{\&alpha;}}{\stackrel{\&RightArrow;}{v}}_{\mathrm{\&alpha;}}(l,m,k),\&CenterDot;\&CenterDot;\&CenterDot;\left(5\right)$

为了简单起见，以下我们假设|λ1|≥|λ2|≥|λ3|，最大本征值λ1相应的本征矢量(

)为

另外，本征矢量的基准应当通过1标准化。

随后，在以上同一步骤中，数据预处理装置3，就某个体素(l，m，k)而言，首先确定一对应方向矢量的相邻体素，该相邻体素与体素(l，m，k)本征矢量

的内积的绝对值在与体素(l，m，k)相邻的26个体素(三维)中是最大的。

我们假设满足该条件的相邻体素是(l+1，m+1，k)。那么由体素(l，m，k)决定的该相邻体素的带有点对称关系的另一相邻体素就被确定了，在此称为体素(l-1，m-1，k)，如图9所示。因为体素(l+1，m+1，k)和(l-1，m-1，k)与体素(l，m，k)本征矢量

的内积是最大的，我们认为它们与体素(l，m，k)的联系是最强的。

如果是如图9(C)所示的情况，那么与体素(l，m，k)相关的体素(l+1，m+1，k)和(l-1，m-1，k)可以被认为是走方相同的神经纤维体素。

接下来在同一步骤中，数据预处理装置3根据以下三种不同情况改写体素(l，m，k)的表达脑功能信息的时间序列数据ρ′(l，m，k，i)。

(情况1)根据预先设定的脑解剖学数据(通过MRI设备50获得的加权图像T1)(参见图9(A))，体素(l+1，m+1，k)和(l-1，m-1，k)都被认为是灰质(神经元细胞体：图片中表示为G)。

在这种情况下，体素(l，m，k)的表达脑功能信息时间序列数据ρ′(l，m，k，i)的数值将根据下式改写：

$\mathrm{\&rho;}(l,m,k,i)=\frac{{\mathrm{\&rho;}}^{\&prime;}(l+1,m,+1,k,i)+{\mathrm{\&rho;}}^{\&prime;}\left(\text{l-1,m-1,k,i}\right)}{2}.\&CenterDot;\&CenterDot;\&CenterDot;\left(6\right)$

(情况2)根据预先设定的脑解剖学数据(通过MRI设备50获得的加权图像T1)(参见图9(B))，体素(l+1，m+1，k)(或者体素(l-1，m-1，k))被认为是灰质。

在这种情况下，体素(l，m，k)的表达脑功能信息时间序列数据ρ′(l，m，k，i)的数值将根据下式改写：

ρ(l，m，k，i)＝ρ′(l+1，m+1，k，i)(orρ′(l-1，m-1，k，i)).…    (7)

(情况3)根据预先设定的脑解剖学数据(通过MRI设备50获得的加权图像T1)(参见图9(C))，体素(l+1，m+1，k)和(l-1，m-1，k)都不被认为是灰质(即被认为白质(神经纤维：图片中表示为W))。

在这种情况下，体素(l，m，k)的表达脑功能信息时间序列数据ρ′(l，m，k，i)的数值将根据下式改写：

$\mathrm{\&rho;}(l,m,k,i)=\frac{{\mathrm{\&rho;}}^{\&prime;}(l+1,m,+1,k,i)+{\mathrm{\&rho;}}^{\&prime;}(l,m,k,i)+{\mathrm{\&rho;}}^{\&prime;}\left(\text{l-1,m-1,k,i}\right)}{3}.\&CenterDot;\&CenterDot;\&CenterDot;\left(8\right)$

所有的体素都将经过这种处理。因此，被认为是灰质的体素(灰质体素)中表达脑功能信息的时间序列数值可以如图9(A)和(B)下部图像所示，该数据值可被反映(传输)给相邻的被认为是白质的体素(白质体素)的表达脑功能信息的时间序列数值。另外，当被认为是白质的体素如图9(C)中下部图像所示彼此相邻时，这些体素中反映脑功能的数将被平滑处理。该预处理可以提高后续步骤的效果。

注意(d)型预处理并不局限于上述所说的实例。但是与灰质体素相邻并从灰质体素中接受反映脑功能信息的时间序列数值的白质体素可以按上述类似步骤进一步将数值传递给邻近的白质体素。此外，当灰质体素中反映脑功能信息的时间序列数值通过上述步骤传递到白质体素时，一加权数值可根据弥散张量D(l，m，k)进行传递。另外，在上述步骤中被改写的白质体素值和原始的白质体素值之间较大的值可以被设定为该白质体素的值。

接下来，在步骤S30中，基于上述所说的弥散张量，

$D(l,m,k)=\left(\begin{array}{ccc}{D}_{\mathrm{ll}}(l,m,k)& D_{\mathrm{lm}}(l,m,k)& D_{\mathrm{lk}}(l,m,k)\\ D_{\mathrm{lm}}(l,m,k)& D_{\mathrm{mm}}(l,m,k)& D_{\mathrm{mk}}(l,m,k)\\ D_{\mathrm{lk}}(l,m,k)& D_{\mathrm{mk}}(l,m,k)& D_{\mathrm{kk}}(l,m,k)\end{array}\right),\&CenterDot;\&CenterDot;\&CenterDot;\left(9\right)$

使用体素间联系度计算装置4，计算两相邻体素联系度矢量，

$\stackrel{\&RightArrow;}{C}(l,m,k)={(c_l(l,m,k),c_m(l,m,k),c_k(l,m,k))}^T\&CenterDot;\&CenterDot;\&CenterDot;\left(10\right)$

该联系度矢量可从下式所得：

$\stackrel{\&RightArrow;}{C}(l,m,k)=(1-\mathrm{\&alpha;})\stackrel{\&RightArrow;}{1}+\mathrm{\&alpha;}{\stackrel{\&RightArrow;}{D}}^{\&prime;}(l,m,k).\&CenterDot;\&CenterDot;\&CenterDot;\left(11\right)$

在此，等式(12)中的每个组份都可如等式(13)所定义，

${\stackrel{\&RightArrow;}{D}}^{\&prime;}(l,m,k)={({D^{\&prime;}}_l(l,m,k),{D^{\&prime;}}_m(l,m,k),{D^{\&prime;}}_k(l,m,k))}^T\&CenterDot;\&CenterDot;\&CenterDot;\left(12\right)$

${{\stackrel{\&RightArrow;}{D}}^{\&prime;}}_l(l,m,k)=\frac{D_l(l,m,k)+D_l(l+1,m,k)}{2},$

${{\stackrel{\&RightArrow;}{D}}^{\&prime;}}_m(l,m,k)=\frac{D_m(l,m,k)+D_m(l,m+1,k)}{2},\&CenterDot;\&CenterDot;\&CenterDot;\left(13\right)$

${{\stackrel{\&RightArrow;}{D}}^{\&prime;}}_k(l,m,k)=\frac{D_m(l,m,k)+D_m(l,m,k+1)}{2}.$

在以下内容中，这被称为相邻体素间平均弥散度矢量。标记α为满足不等式(14)的参数。

是满足等式(15)的常数矢量。

0≤α≤1，…    (14)

$\stackrel{\&RightArrow;}{1}={\left(\mathrm{1,1,1}\right)}^T.\&CenterDot;\&CenterDot;\&CenterDot;\left(15\right)$

在等式(13)中，等式(16)是代表质子弥散度的矢量(以下被称为弥散度矢量)，其中第一组份Dl(l，m，k)代表质子在体素(l，m，k)X轴方向上的弥散度，第二组份Dm(l，m，k)代表质子在体素(l，m，k)Y轴方向上的弥散度，第三组份Dk(l，m，k)代表质子在体素(l，m，k)Z轴方向上的弥散度。

$\stackrel{\&RightArrow;}{D}(l,m,k)={(D_l(l,m,k),D_m(l,m,k),D_k(l,m,k))}^T\&CenterDot;\&CenterDot;\&CenterDot;\left(16\right)$

这里所说的X，Y和Z方向假定和图7显示的方向相一致。在等式(9)中，弥散张量被假定为对称张量。

使用具有上述假设的弥散度矢量

通过等式(13)形成的等式(11)的相邻体素间联系度矢量

的每个组份Cl(l，m，k)，Cm(l，m，k)和Ck(l，m，k)分别代表体素(l，m，k)和体素(l+1，m，k)的联系度，体素(l，m，k)和体素(l，m+1，k)的联系度，体素(l，m，k)和体素(l，m，k+1)的联系度。

需要注意的是在本实施例中，虽然在倾斜方向上相邻的体素(如体素(l，m，k)和体素(l+1，m+1，k))的联系度没有被直接考虑，但是其已经通过垂直方向上的体素(体素(l+1，m，k)或者体素(l，m+1，k))被间接考虑了。显然，我们也可以设置成直接考虑倾斜方向上相邻的体素的联系度。

弥散度矢量

可以由以下等式计算：

(I)当被应用于椭圆形的模型中，

$D_l(l,m,k)=\sqrt{\frac{2\left|D\right|}{D_{\mathrm{mm}}(l,m,k)D_{\mathrm{kk}}(l,m,k)-D_{\mathrm{mk}}{(l,m,k)}^2}},$

$D_m(l,m,k)=\sqrt{\frac{2\left|D\right|}{D_{\mathrm{ll}}(l,m,k)D_{\mathrm{kk}}(l,m,k)-D_{\mathrm{lk}}{(l,m,k)}^2}},\&CenterDot;\&CenterDot;\&CenterDot;\left(17\right)$

$D_k(l,m,k)=\sqrt{\frac{2\left|D\right|}{D_{\mathrm{mm}}(l,m,k)D_{\mathrm{mm}}(l,m,k)-D_{\mathrm{lm}}{(l,m,k)}^2}}$

这里的|D|代表由等式(9)代表的矩阵的行列式。

(II)使用由等式(9)代表的矩阵的对角元素Dll(l，m，k)，Dmm(l，m，k)和Dkk(l，m，k)，可以得到以下等式：

D_l(l，m，k)＝|D_ll(l，m，k)|，D_m(l，m，k)＝|D_mm(l，m，k)|and D_k(l，m，k)＝|D_kk(l，m，k)|

…    (18)

(III)使用一代表弥散各向异性的指标，弥散各向异性典型代表形式包括下式所示的部分各向异性：

$\sqrt{\frac{3}{2}\&CenterDot;\frac{\underset{\mathrm{\&alpha;}=1}{\overset{3}{\mathrm{\&Sigma;}}}{({\mathrm{\&lambda;}}_{\mathrm{\&alpha;}}-D_{\mathrm{AV}})}^2}{\underset{\mathrm{\&alpha;}=1}{\overset{3}{\mathrm{\&Sigma;}}}{{\mathrm{\&lambda;}}_{\mathrm{\&alpha;}}}^2}}\&CenterDot;\&CenterDot;\&CenterDot;\left(19\right)$

和相对各向异性：

$\frac{1}{\sqrt{6}}\sqrt{\frac{\underset{\mathrm{\&alpha;}=1}{\overset{3}{\mathrm{\&Sigma;}}}{({\mathrm{\&lambda;}}_{\mathrm{\&alpha;}}-D_{\mathrm{AV}})}^2}{D_{\mathrm{AV}}}},\&CenterDot;\&CenterDot;\&CenterDot;\left(20\right)$

弥散度矢量

同样也可以被确定。

这里的DAV是

$D_{\mathrm{AV}}=\frac{\underset{\mathrm{\&alpha;}=1}{\overset{3}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_{\mathrm{\&alpha;}}}{3}.\&CenterDot;\&CenterDot;\&CenterDot;\left(21\right).$

(IV)另外，可根据预先设定的脑解剖数据(通过MRI设备50获得的T1加权图像)确定脑白质体素和灰质体素，然后上述涉及的(I)的弥散度用于白质体素，而将其它弥散度矢量(如：以上提到的(II)或(III))用于灰质体素。

如上所述，作为本发明中的弥散度矢量

任何可以通过MRI设备50获得的弥散张量D(l，m，k)确定神经纤维走向的矢量都可以被应用为本发明的弥散矢量。

另外，除了与α参数相关的不等式(14)以外，α参数也可以被设置为在白质体素和灰质体素间转换。比如，在灰质体素中拥有一个大的值，而在白质体素中拥有一个小的值。

此外，除了等式(11)的联系度矢量

以外，也可以使用以下等式的联系度矢量：

$\stackrel{\&RightArrow;}{C}(l,m,k)=\stackrel{\&RightArrow;}{1}+\mathrm{\&alpha;}{\stackrel{\&RightArrow;}{D}}^{\&prime;}(l,m,k)\&CenterDot;\&CenterDot;\&CenterDot;\left(22\right)$

另外，为了突出神经纤维的走向，可以提高等式(11)或者(22)中联系度矢量

的每个组分cl(l，m，k)，cm(l，m，k)和ck(l，m，k)至其N-th次方(N为2或2以上自然数)，从而得到一新的联系度矢量。

另外，通过将这些组分数值乘以a得到的值可能用于白质体素中，通过将这些组分数值乘以b得到的值可能用于灰质体素中。此外，只有这些组分中最大的组分才可以用上述方法突出。还有，只有拥有大于或者等于{cl(l，m，k)+cm(l，m，k)+ck(l，m，k)}/3数值的组分才可以用上述方法突出。

联系度矢量

可以用一定方式形成以使灰质体素只与相邻的六个体素相关联(六个方向)，而白质体素只与和本征等式(5)中最大本征值对应的本征矢量的内积为最大值的体素以及与该体素拥有点对称关系的体素(两个方向)相关联。

之后，在步骤S40，为了根据步骤S20中预处理后的反映脑功能信息的时间序列数据ρ(l，m，k，i)和步骤S40中计算得到的联系度矢量

对图5(A)所示数据进行非参数型回归分析，如果存在六个相邻体素，数据估计值形成装置5将形成如下数据估计值Q：

$Q=\frac{1}{I}\underset{i=1}{\overset{I}{\mathrm{\&Sigma;}}}{\{\mathrm{\&phi;}\left(i\right)-\widehat{\mathrm{\&phi;}}\left(i\right)\}}^2$

$+\mathrm{\&kappa;}[\underset{l=1}{\overset{L-1}{\mathrm{\&Sigma;}}}{C}_l(l,m,k){\{\hat{a}(l+1,m,k)-\hat{a}(l,m,k)\}}^2$

$+\underset{m=1}{\overset{M-1}{\mathrm{\&Sigma;}}}{C}_m(l,m,k){\{\hat{a}(l,m+1,k)-\hat{a}(l,m,k)\}}^2$

$+\underset{k=1}{\overset{K-1}{\mathrm{\&Sigma;}}}{C}_k(l,m,k){\{\hat{a}(l,m,k+1)-\hat{a}(l,m,k)\}}^2].\&CenterDot;\&CenterDot;\&CenterDot;\left(23\right)$

这里的标记i为时间，标记l，m和k分别为代表体素在X轴，Y轴，Z轴方向上位置的变量。请注意，实际中图7所示的立体构型体素中可能有体素不代表脑(神经)，与其相关的l，m，k需被去除。同时，

是图5(A)代表CLASS的变量的时间序列，它可以为1(＝ON)或0(＝OFF)。在本实施例的非参数型回归分析中，当图5(A)所示的反映脑功能信息时间序列数据ρ(l，m，k，i)和代表CLASS的变量的时间序列的为解释变量时，我们在这些变量间假设一个线性回归等式，其中为回归系数：

$\widehat{\mathrm{\&phi;}}\left(i\right)=\underset{l=1}{\overset{L}{\mathrm{\&Sigma;}}}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}\hat{a}(l,m,k)\mathrm{\&rho;}(l,m,k,i)\&CenterDot;\&CenterDot;\&CenterDot;\left(24\right)$

这里的

和分别表示预测结果，k是代表相邻体素间连续性的加权参数。

虽然如上述所说的，等式(23)通过相邻体素间的连续性推算出了估计值，但是该估计值也可以通过相邻体素间的平滑性推算出。

本实施例中，以非参数回归分析作为数据分析技术的基础在于，当解释后的变量

与解释变量ρ(l，m，k，i)相比，的数量是非常小的。在图5(A)所示的实例中出现了64×64×64＝262144个解释变量ρ(l，m，k，i)，该数目与体素的数量相同，然而解释后变量

的检测值只有18个。一般来说，解释后变量的检测值为10个，也可能超过100个。因为在这种情况下，无法进行常规的回归分析，本发明中我们使用了非参数回归分析(非专利文献6：“Spline Smoothing and NonparametricRegression″，R.L.Eubank，Marcel Dekker，Newyork，1988”)。

通常在非参数回归分析中，我们通过假定解释后变量

的连续性来确定使估计值Q′(等式(25))最小化的

$Q^{\&prime;}=\frac{1}{I}\underset{i=1}{\overset{I}{\mathrm{\&Sigma;}}}{\{\mathrm{\&phi;}\left(i\right)-\widehat{\mathrm{\&phi;}}\left(i\right)\}}^2+\mathrm{\&kappa;}\underset{i=1}{\overset{I-1}{\mathrm{\&Sigma;}}}{\{\widehat{\mathrm{\&phi;}}(i+1)-\widehat{\mathrm{\&phi;}}\left(i\right)\}}^2\&CenterDot;\&CenterDot;\&CenterDot;\left(25\right)$

这里右边第一项表示解释后变量在每个时间i的预测结果

的误差，第二项是用于估计解释后变量

在某个时间方向上的连续性。

然而如图5(A)所示，我们不能就解释后变量

假定连续性。因此，在本实施例中等式(23)表示的数据估计值Q是根据增强的非参数回归分析得出的，这不是对解释后变量

而是对解释变量ρ(l，m，k，i)设置的连续性约束(非专利文献7：“Nonparametric regression analysisof brain function images”，H.Tukimoto et al.，Institute of Electronics，Information and Communication Engineers D-II，Vol.J84，No.12，pp2623-2633，December，2001)。

在此，我们将对等式(23)进行解释。右边第一项表示解释后变量

在每个时间i与预测结果

的误差，这与等式(25)相同。

第二到第四项表示相邻体素对解释后变量

影响的连续性，该连续性可以通过相邻解释变量ρ(l，m，k，i)对解释后变量的影响进行估计，也就是等式(24)的回归系数的连续性。

另外，在本实施例中，由体素间联系度计算装置4计算得到的联系度矢量

的每一个组份Cl(l，m，k)，Cm(l，m，k)和Ck(l，m，k)分别在第二至四项中作为加权因子引入相邻体素间在X轴，Y轴，Z轴方向上的连续性中。

由于这些组分分别是平均弥散度矢量

的每个组份D′l(l，m，k)，D′m(l，m，k)和D′k(l，m，k)的函数，如等式(11)所示，因此等式(23)的数据估计值Q能反映神经纤维(神经元)的走向(质子的弥散各向异性)。

接下来，在步骤S50中，数据分析装置6将通过计算来确定在S40步中形成的数据估计值Q的极值。在本实施例中，我们计算获得了使等式(23)中数据估计值Q最小化的回归系数

这里为了简便起见，我们从等式(23)改写数据估计值Q″为：

$Q^{\&prime;\&prime;}=\frac{1}{I}\underset{i=1}{\overset{I}{\mathrm{\&Sigma;}}}{\{\mathrm{\&phi;}\left(i\right)-\widehat{\mathrm{\&phi;}}\left(i\right)\}}^2+\mathrm{\&kappa;}\underset{n=1}{\overset{L^3-1}{\mathrm{\&Sigma;}}}C\left(n\right){\{\hat{a}\left(n^{\&prime;}\right)-\hat{a}\left(n\right)\}}^2\&CenterDot;\&CenterDot;\&CenterDot;\left(26\right)$

其中，n是代表体素位置的独立的正数变量，所有具有三维构型的体素都被给予一预先指定的数值，从而被重新编排为一维，该变量与三维显示变量(l，m，k)一一对应。另外，为了简化目的，我们设置为L＝M＝K，因此体素的数量为L3。同时，当所有的体素都以这种方法在一维方向上重排后，系数C(n)代表体素n和相邻体素n′的联系度。该系数是根据等式(11)的连接矢量

特别确定的。

例如，当只限于图3(B)中的二维切片图像Sk时，与n＝130体素相邻的体素有四个：n＝127，129，131，257。例如：根据S30步计算的联系度矢量

体素n＝127和体素n＝130，体素n＝129和体素n＝130，体素n＝131和体素n＝130，体素n＝257和体素n＝130间的联系度分别为2，2，5，5。

在这里，联系度是由

所获得。

该实例中，通过使用最小平方法，回归系数

可以被确定如下：

$\stackrel{\&RightArrow;}{a}={(X^{T}X+I{\mathrm{\&kappa;}}_{O}C)}^{-1}{X}^T\stackrel{\&RightArrow;}{\mathrm{\&phi;}}\&CenterDot;\&CenterDot;\&CenterDot;\left(27\right)$

其中，X为拥有ρ(n，i)组分的L3xI矩阵，

是拥有组分

的I维矢量，

是拥有如下

组分的L3维矢量，κO是κ的优化值。κ由交叉确认方法确定，即κ使等式(28)最小化。

型的神经网络模型包括由输入层，中间层，和输出层组成的多层次神经网络。

在本修饰例的步骤S40中，数据估计值形成装置5形成估计值，其中基于体素间联系度计算装置4获得的联系度矢量

所计算得到的评估相邻体素间连续性的估计值被加入到平方误差值中，该平方误差是神经网络模型中反映正常学习的估计值。

随后，在步骤S50中，使用误差后续增值法或类似方法，通过数据分析装置6进行分析从而确定估计值的最小值，而后步骤S60中，图片生成装置7生成输入(脑功能信息的时间序列数据)和输出的敏感度(即输入值的改变对输出值的影响)，从而生成类似于实施例一的显示图像。然后，在S70步中，图片显示装置8显示上述生成的立体图像。

这使我们有可能获得与实施例一相似的效果。

<实施例一修饰例2>

在步骤S40中，数据估计值形成装置5使用实施例一中的反映脑信息的时间序列数据ρ(l，m，k，i)作为解释变量，

作为解释后变量，形成等式(23)的估计值Q，其中，在本修饰例中，每个体素(l，m，k)都有一个形成的估计值Q″′如下：

$Q^{\&prime;\&prime;\&prime;}=\frac{1}{I}\underset{i=1}{\overset{I}{\mathrm{\&Sigma;}}}{\{\mathrm{\&rho;}(l,m,k,i)-\widehat{\mathrm{\&rho;}}(l,m,k,i)\}}^2$

$+\frac{\mathrm{\&kappa;}}{I}\underset{i=1}{\overset{I}{\mathrm{\&Sigma;}}}[C_l(l,m,k){\{\widehat{\mathrm{\&rho;}}(l+1,m,k,i)-\widehat{\mathrm{\&rho;}}(l,m,k,i)\}}^2$

$+C_m(l,m,k){\{\widehat{\mathrm{\&rho;}}(l,m+1,k,i)-\widehat{\mathrm{\&rho;}}(l,m,k,i)\}}^2$

$+C_k(l,m,k){\{\widehat{\mathrm{\&rho;}}(l,m,k+1,i)-\widehat{\mathrm{\&rho;}}(l,m,k,i)\}}^2]\&CenterDot;\&CenterDot;\&CenterDot;\left(29\right)$

该估计值的形成是通过使用反比例关系，即使用

作为解释变量，ρ(l，m，k，i)作为解释后变量。

其中给出了下列限制条件：

$\widehat{\mathrm{\&rho;}}(l,m,k,i)=\underset{l=1}{\overset{L}{\mathrm{\&Sigma;}}}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}\hat{a}(l,m,k)\mathrm{\&phi;}\left(i\right)+b\&CenterDot;\&CenterDot;\&CenterDot;\left(30\right)$

随后在步骤S50中，使用交叉证实方法，通过数据分析装置6确定最佳值κ，从而确定上述估计值Q″′的最小值，使用概率性搜寻算法如：起源算法(Genetic Algoritm)，确定回归系数

随后在步骤S60中，根据体素的脑功能数据ρ(l，m，k)，通过图像生成装置7生成一个作为图像显示的数值，该数值通过不同分析所得结果被确定为重要的。然后，在S70中，通过图像显示装置8显示以上生成的立体图像。

这使我们有可能获得与实施例一相似的效果。

<实施例一修饰例3>

在步骤S40中，数据估计值形成装置5使用实施例一中的反映脑信息的时间序列数据ρ(l，m，k，i)作为解释变量，

作为解释后变量，形成等式(23)的估计值Q，其中，在本修饰例中，我们也可以使用

作为解释变量，ρ(l，m，k，i)作为解释后变量，使用多层神经网络模型进行非线性回归分析。请注意在本修饰例中，误差后续增值法不可以在S50中用于确定极值，因此，数据分析装置6使用交叉证实法确定最佳κ值，通过概率性搜寻算法如起源算法进行数据分析。随后在步骤S60中，根据体素的脑功能数据ρ(l，m，k)，通过图像生成装置7生成一个作为图像显示的数值，该数值通过上述不同分析所得结果被确定为重要的。然后，在S70中，通过图像显示装置8显示以上生成的立体图像。

这使我们有可能获得与实施例一相似的效果。

<实施例一修饰例4>

实施例一S50中使用等式(24)的线性回归等式，通过数据分析装置6进行非参数型回归分析，也可以使用分类分析技术，如：判别分析(非专利文献8：Story of multivariate analysis″，T.Arima，S.Ishimura，Tokyo Tosho，1987)，定量II(非专利文献9)，决策树(非专利文献10：″Data analysis by A.I.″，J.R.Quinlan，Toppan，1995)，支撑矢量机(非专利文献11：″Introduction to support vector machine″，NelloCristianini and John Shawe-Taylor，Kyoritsu shuppan，2005)或者类似技术作为分析技术。

在任何情况下，使用反映脑功能信息时间序列数据ρ(l，m，k，i)作为解释变量，使用任务和休息的时间序列数据作为外部标准，通过各种辨别方程式，决策树或类似方程式对数据进行分类。在本修饰例步骤S40中，通过数据估计值形成装置5形成估计值，其中于步骤S30中计算得到的联系度矢量

被加入到每个分析技术的普通估计值中。

随后在步骤S50中，使用每个分析技术，通过数据分析装置6进行数据分析。如果估计值出现上述改变，普通的解决方案将不适用，因此可以使用交叉证实法，通过数据分析装置6确定优化κ值，使用概率性搜寻算法进行数据分析。

随后在S60中，根据体素的脑功能数据ρ(l，m，k)，通过图像生成装置7生成一个用于图像显示数值，该数值通过上述不同分析结果被确定是重要的。然后，在步骤S70中，通过图像显示装置8显示以上生成的立体图像。

这使我们有可能获得与实施例一相似的效果。

<实施例一修饰例5>

当使用实施例一步骤S50中等式(24)的线性回归等式，通过数据分析装置6进行非参数回归分析，独立组份分析也可以被用作一种数据分析技术。(非专利文献12：″Detailed explanation of IndependentComponent Analysis Novel world of signal analysis″，Aapo Hyvarinenet al.，Tokyo Denki University Press 2005)。

在本修饰例中，数据估计值形成装置5在步骤S40中基于体素间联系度计算装置4获得的联系度矢量

形成估计值。

在本修饰例中，根据独立成份分析，分析反映脑功能信息时间序列数据ρ(l，m，k，i)中的“空间独立性”条件被更改为“基于弥散张量信息的联系度的空间独立性”。独立组份分析有几种特定技术，但是在使用交叉关联度最小化的技术时，交叉关联度其实并不是被最小化了，而是在联系度的基础上使其更接近估计值。结果，具有较大联系度的体素间的交叉关联度变成了“独立程度”，而具有较小联系度(或没有关联度)的体素间的交叉关联度被最小化，与联系度的估计值更相近，更确切的讲，例如：在评估方程中，基于联系度的估计值的加权因子被设置于平方总数的每一个值，该评估方程使关联矩阵中的非对角元素的平方总和最小化。当联系度越大，加权因子越小。

在随后的步骤S50中，结合上述形成的估计值，通过分析装置6进行独立组份分析。其中，通过使用交叉证实方法确定优化κ值，同时通过概率性搜寻算法如起源算法进行数据分析。

随后在步骤S60中，通过图像生成装置7生成上述获得的独立组分显示图像。然而，由于上述分析产生了许多独立组分，因此需要从所有的行矢量中提取与“任务/休息”关联值最高的行矢量。

在随后的步骤S70中，通过图像显示装置8显示上述形成的立体图像。

这使我们有可能获得与实施例一相似的效果。

请注意，除了主要的组分分析以外，作为从反映脑功能信息时间序列数据中提取预设主要组份的提取技术的因素分析和定量III，但不是上述的独立组份分析，也可以以上述相似的方式进行分析。

[实施例二]

图10是本发明实施例二中脑功能分析仪器构造示意图。实施例二中的脑功能分析仪器20装备有体素间联系度计算装置24，而不是体素间联系度计算装置4，还有实施例一中脑功能分析仪器10的数据估计值形成装置5和数据分析装置6，数据估计值形成装置25，数据平滑装置200和数据分析装置26。对于与实施例一中相同的配置，其符号也相同，因此之后将不做重复描述。请注意，虽然脑功能分析仪器20使用显示控制台模式，其中本实施例中的主要形成装置20A，图形显示装置8，存储装置9被组合在一起，但是也可以使用如下配置：将图形显示装置8，或者图形显示装置8和存储装置9，和主要形成装置20A分开，分别作为独立的图像显示单元和储存单元。

依照实施例一详细描述的方式，从弥散张量数据获取装置2获取的弥散张量数据D(l，m，k)(等式(9))中计算出平均弥散度矢量(等式(12))后，体素间联系度计算装置24计算联系度矢量

然后据根据平均弥散度矢量

的要素D′l(l，m，k)、D′m(l，m，k)和D′k(l，m，k)的大小由数据预处理装置3来预处理脑功能信息的时间序列数据ρ(l，m，k，i)，再对所得数据进行平滑处理，从而决定相邻体素间的平滑程度。

数据平滑装置200包括一个加权平均过滤器，使用体素间联系度计算装置24计算得到的联系度矢量

来平滑处理由数据预处理装置3预处理后的脑功能信息时间序列数据ρ(l，m，k，i)，然后使用加权平均过滤器平滑处理由数据预处理装置3预处理后的脑功能信息时间序列数据ρ(l，m，k，i)。

数据分析装置26对由数据平滑装置200利用如SPM或相似技术平滑处理后的数据进行数据分析。

图11是脑功能分析仪器20根据图10所示实施例二进行数据分析的流程图。

在步骤S100中，脑功能数据获取装置1获取由MRI设备50测得的脑功能信息原始时间序列数据ρ′(l，m，k，i)。同一步骤中，弥散张量信息获取装置2获取由MRI设备50测得的弥散张量数据D(l，m，k)。

接着，在步骤S110中，数据预处理装置3对在步骤S100中获得的脑功能信息原始时间序列数据ρ′(l，m，k，i)进行如实施例一所述的(a)到(d)型预处理，进而产生脑功能信息时间序列数据ρ(l，m，k，i)。

接着，在步骤S120中，体素间联系度计算装置24按照实施例一详细介绍的方法，根据步骤S100中获得的弥散张量数据D(l，m，k)初步计算平均弥散度矢量

接着，根据等式(31)计算联系度矢量

进而通过平滑预处理的脑功能信息时间序列数据ρ(l，m，k，i)确定相邻体素间的平滑性，该脑功能信息时间序列数据ρ(l，m，k，i)根据平均弥散度矢量

的要素D′l(l，m，k)、D′m(l，m，k)和D′k(l，m，k)的大小由数据预处理装置3进行预处理。

$\stackrel{\&RightArrow;}{W}(l,m,k)=(w_l(l,m,k),w_m(l,m,k),w_k(l,m,k))$

$\&equiv;(\frac{w{D^{\&prime;}}_l(l,m,k)}{\underset{\mathrm{\&alpha;}=l,m,k}{\mathrm{\&Sigma;}}{D^{\&prime;}}_{\mathrm{\&alpha;}}(l,m,k)},\frac{w{D^{\&prime;}}_m(l,m,k)}{\underset{\mathrm{\&alpha;}=l,m,k}{\mathrm{\&Sigma;}}{D^{\&prime;}}_{\mathrm{\&alpha;}}(l,m,k)},\frac{w{D^{\&prime;}}_k(l,m,k)}{\underset{\mathrm{\&alpha;}=l,m,k}{\mathrm{\&Sigma;}}{D^{\&prime;}}_{\mathrm{\&alpha;}}(l,m,k)}).\&CenterDot;\&CenterDot;\&CenterDot;\left(31\right)$

接着，在步骤S130中，数据平滑装置26根据等式(32)，使用步骤S120中计算获得的联系度矢量

各向异性地平滑在步骤S110中预处理后的脑功能信息时间序列数据ρ(l，m，k，i)。

$\stackrel{\&OverBar;}{\mathrm{\&rho;}}(l,m,k,i)=\frac{1}{w+W}[\{\mathrm{W\&rho;}(l,m,k,i)$

$+w_l\&CenterDot;\{\mathrm{\&rho;}(l-1,m,k,i)+\mathrm{\&rho;}(l+1,m,k,i)\}$

$+w_m\&CenterDot;\{\mathrm{\&rho;}(l,m-1,k,i)+\mathrm{\&rho;}(l,m+1,k,i)\}$

$+w_k\&CenterDot;\{\mathrm{\&rho;}(l,m,k-1,i)+\mathrm{\&rho;}(l,m,k+1,i)\}].\&CenterDot;\&CenterDot;\&CenterDot;\left(32\right)$

此处，w和W为加权系数，该加权系数经过标准化以满足以下等式(33)：

$\left\{\underset{\mathrm{\&alpha;}=l,m,k}{\mathrm{\&Sigma;}}{D^{\&prime;}}_{\mathrm{\&alpha;}}(l,m,k)\right\}+W=w+W=1.\&CenterDot;\&CenterDot;\&CenterDot;\left(33\right)$

图12描述了实施一使用的一数据平滑技术，同时列出了二维切片图像Sk在位置(l，m)的9个体素的脑功能信息时间序列数据ρk(l，m，i)。在该例子中，ρk(l，m，i)＝5，ρk(l-1，m，i)＝ρk(l-1，m-1，i)＝ρk(l，m-1，i)＝ρk(l，m+1，i)＝ρk(l+1，m，i)＝2，and ρk(l+1，m-1，i)＝ρk(l-1，m+1，i)＝ρk(l+1，m+1，i)＝3，此处，ρk(l，m-1，i)，ρk(l，m+1，i)，ρk(l-1，m，i)，ρk(l+1，m，i)用来平滑其中的ρk(l，m，i)。同时，假设一种满足D ′k(l，m，k)＝0的情况(也就是说在Z方向上不存在质子发散)，和D′l(l，m，k)＝0.8、D′m(l，m，k)＝0.2，w＝0.4，W＝0.6，在这种情况下，联系度矢量

根据等式(31)得到等式(34)所示结果，以使得位于(l，m，k)处的体素的脑功能信息时间序列数据ρ(l，m，k，i)＝5根据等式(36)平滑处理为等式(35)。

$\stackrel{\&RightArrow;}{W}(l,m,k)=(0.4\&times;\mathrm{0.8,0.4}\&times;\mathrm{0.2,0})=\left(\mathrm{0.32,0.08,0}\right)\&CenterDot;\&CenterDot;\&CenterDot;\left(34\right)$

$\stackrel{\&OverBar;}{\mathrm{\&rho;}}(l,m,k,i)=\frac{1}{1}(0.6\&times;5+0.32\&times;2\&times;2+0.08\&times;2\&times;2)=4.44\&CenterDot;\&CenterDot;\&CenterDot;\left(35\right)$

$\stackrel{\&RightArrow;}{S}(l,m,k)=(S_l(l,m,k),S_m(l,m,k),S_k(l,m,k))={\stackrel{\&RightArrow;}{D}}^{\&prime;}(l,m,k)-\stackrel{\&RightArrow;}{h}(l,m,k)\&CenterDot;\&CenterDot;\&CenterDot;\left(\text{36}\right)$

这一处理是针对二维切片图像Sk中所有体素的脑功能信息时间序列数据ρ(l，m，k，i)进行的，经过对该数据的平滑处理，使得所有方向(也就是平均弥散度矢量

的各要素D′l(l，m，k)、D′m(l m，k)和D′k(l，m，k))中拥有最大组份(本实施例中为X方向)的方向，即神经纤维的走向，的作用变大。

另外，如等式(5)所述，只有与最大本征值λM相应的本征矢量

具有相同方向的两个体素才可以得到平滑化。

接着，在步骤S140中，数据分析装置26用SPM或相似的分析技术对在步骤S140中经过平滑的脑功能信息的时间序列数据ρ(l，m，k，i)进行数据分析。

接着，在步骤S150中，图像生成装置7根据数据分析装置26分析的结果生成颜色图像，如灰色或全彩图。

其后，在步骤S160中，图像显示装置8立体显示上述生成的图像。

如上所述，在实施例二中，从MRI设备中获得的脑功能信息时间序列数据依据从MRI设备获得的弥散张量数据计算得到的联系度矢量进行平滑处理，在考虑脑激活区域之间连接结构的同时，确定脑激活区域。

需要指出的是，等式(31)的联系度矢量

或平滑处理等式(32)的运算法则都只是其中的一个例子，其它的实施方式也是可能的。

此外，上述分析清楚地表明在本实施例中的数据平滑处理技术并不依赖于某一特定的数据分析技术。

<实施例二修饰例>

在实施例二中，尽管步骤S140使用SPM或相似技术进行的数据分析在后，步骤S130依据弥散张量数据D(l，m，k)计算得到的联系度矢量

对脑功能信息时间序列数据ρ(l，m，k，i)进行平滑处理在前，颠倒这些步骤的顺序也可以达到本发明目的。

因而，我们将对上述实施例二的修饰例进行描述。为避免重复，我们在以下内容中只描述与实施例二不同的部分。

图13是一个符合实施例二的修饰例的脑功能分析仪器构造示意图。尽管本修饰例的脑分析仪器20′的配置与实施例二的仪器配置基本相似，但不同的是数据平滑装置200和数据分析装置26被互相调换了。因此，由脑功能分析仪器20′进行的数据分析程序与实施例二中的不同。注意，尽管脑功能分析仪器20′采用了显示控制台模式，其中主要形成装置20A′、图像显示装置8和存储装置9在本实例中是组合在一起的，其也可能采用另外一种构造，即图像显示装置8或图像显示装置8和图像存储装置9从主要形成装置20A′中分离出来分别作为独立的图像显示单元和存储单元。

图14是图13中脑功能分析仪器20′进行数据分析程序的流程图。在本修饰例中，数据分析装置26′首先在步骤S121中使用SPM或相似分析技术对数据预处理装置3预处理后的脑功能信息时间序列数据ρ(l，m，k，i)进行数据分析。接着，在步骤S131中，通过体素间联系度计算装置24，按照与实施例二中相似的方法计算如等式(31)所示的联系度矢量

接着，在步骤S141中，数据平滑装置200′利用等式(32)按照与实施例二中相似的方法对脑功能数据ρ(l，m，k)进行平滑处理，该数据由数据分析装置26′进行分析产生的。接着，在步骤S150中，图像生成装置7根据预设的处理方法生成有关上述平滑处理数据ρ(l，m，k)的图像颜色，如灰色或全彩图像。

此后，在步骤S160中，图像显示装置8立体显示上述产生的图像。

因此，本修饰例获得与实施例二相似的结果。

[实施例三]

图15是本发明实施例三的脑功能分析仪器的构造示意图。实施例三的脑功能分析仪器30包括如下配置：体素间联系度计算装置34，而不是体素间联系度计算装置4，数据估计值形成装置5，实施例一中脑功能分析仪器10的数据分析装置6，聚集装置300和数据分析装置36。与实施例一配置相同的，其符号也相同，下文将省略重复的内容。注意，尽管脑功能分析仪器30采用了显示控制台模式，其中主要形成装置30A、图像显示装置8和存储装置9在本实例中是组合在一起的，其也可能采取如下配置：图像输出装置8，或图像输出装置8和图像存储装置9从主要形成装置30A中分离出来分别作为独立的图像显示单元和存储单元。

通过弥散张量数据获取装置2，根据实施例一所述的方法获得的弥散数据D(l，m，k)(等式(9))计算出平均弥散度矢量

后，体素间联系度计算装置34根据平均弥散度矢量的各要素D′l(l，m，k)，D′m(l，m，k)和D′k(l，m，k)的大小计算出一个联系度矢量

作为聚集各体素的一个指数。

聚集装置300使用由体素间联系度计算装置34计算得出的联系度矢量

聚集体素。

数据分析装置36利用SPM或相似分析技术对数据预处理装置3预处理后的脑功能信息时间序列数据ρ(l，m，k，i)进行数据分析。该数据分析同时使用聚集装置300聚集获得的体素组(群)为单位。

图16是图15所示实施例三中脑功能分析仪器30进行数据分析的流程图。

在步骤S200中，脑功能数据获取装置1获取MRI设备50测得的脑功能信息原始时间序列数据ρ′(l，m，k，i)。同时，弥散张量信息获取装置2获取MRI设备50测得的弥散张量数据D(l，m，k)。

接着，在步骤S210中，数据预处理装置3对在步骤S200中获得的脑功能信息原始时间序列数据ρ′(l，m，k，i)进行如实施例一所述的(a)到(d)型预处理，进而产生脑功能信息时间序列数据ρ(l，m，k，i)。

接着，在步骤S220中，体素间联系度计算装置34根据在实施例一所述方式，首先从在步骤S200中获得的弥散张量数据D(l，m，k)中计算出平均弥散度矢量

接着，如等式(36)所示，计算出联系度矢量该矢量为根据平均弥散度矢量

的各要素D′l(l，m，k)，D′m(l，m，k)和D′k(l，m，k)进行体素聚集的指数。

$\stackrel{\&RightArrow;}{S}(l,m,k)=(S_l(l,m,k),S_m(l,m,k),S_k(l,m,k))={\stackrel{\&RightArrow;}{D}}^{\&prime;}(l,m,k)-\stackrel{\&RightArrow;}{h}(l,m,k)\&CenterDot;\&CenterDot;\&CenterDot;\left(36\right)$

此处，

是一个代表聚集临界值的常数矢量，表示为等式(37)：

$\stackrel{\&RightArrow;}{h}(l,m,k)=(h_x,h_y,h_z).\&CenterDot;\&CenterDot;\&CenterDot;\left(37\right)$

接着，在步骤S230中，通过聚集装置300聚集体素，其中在步骤S220中计算得到的联系度矢量的各要素Sl(l，m，k)，Sm(l，m，k)和Sk(l，m，k)的值分别在X轴、Y轴和Z轴方向变为正数。

图17为使用图15所示聚集装置300的聚集技术图解，并显示二维切片图像Sk的16个体素。图17(A)中相邻体素间的箭头值代表平均弥散度矢量在X轴或Y轴方向的值。如果等式(37)的常数矢量的每一要素的值分别设定为hl＝hm＝3和hk＝0，联系度矢量

的要素Sl(l，m，k)和Sm(l，m，k)的值也变为如图17(B)中所示。

此时，通过有次序的连接相邻体素，其中联系度矢量

的要素Sl(l，m，k)和Sm(l，m，k)在X轴和Y轴方向均为正值，则被虚线包围的体素将最终形成一束。

由于联系度矢量

变为平均弥散度矢量D′(l，m，k)的函数，所以上述聚集可以反映神经纤维的走向。

接着，在步骤S240中，数据分析装置36使用SPM或相似分析对在步骤S210中预处理后的脑功能信息时间序列数据ρ(l，m，k，i)进行数据分析，该数据分析覆盖了步骤S230中聚集的每一体素组(束)。在那种情况下则需要对一个聚集束或类似聚集等同物内的脑功能信息时间序列数据进行平均化处理。

接着，在步骤S250中，根据步骤S240分析后的数据

通过图像生成装置7，在预设处理规则下，生成颜色图像，如：灰色或全彩图。

之后，在步骤S260中，图像显示装置8显示步骤S250产生的立体图像。

如上所述，在实施例三中，由MRI设备获取的脑功能信息时间序列数据根据该设备获取的弥散张量数据被聚集成束，而后对各个束相应的脑功能信息时间序列数据进行分析。这样就能在确定脑激活功能区域的同时也考虑到脑激活功能区域间联系。

需要指出的是，等式(36)的联系度矢量

只是其中的一个例子，也可能使用其它的实施方式。

另外，如等式(5)所述，只有与最大本征值λM相应的本征矢量

具有相同方向的两个体素才可以被聚集在一起。

此外，上述分析清楚地表明了本实施例中的聚集处理技术并不依赖于一个具体的数据分析技术。

<实施例三的修饰例1>

虽然在实施例三中，首先使用根据步骤S230中获得的弥散张量数据D(l，m，k)计算联系度矢量

聚集体素，然后就每一个聚集得到的束在步骤S240使用SPM或类似技术对脑功能信息时间序列数据ρ(l，m，k，i)进行分析。但颠倒这些步骤的顺序也可以达到本发明目的。

我们将对上述实施例三中的修饰例进行阐述。为避免重复，下文只阐述那些与实施例三中不同的内容。

图18是实施例三修饰例中的脑功能分析仪器的结构示意图。虽然本修饰例的脑功能分析仪器30′的结构基本上与实施例三相同，但是它们区别在于聚集装置300和数据分析装置36被调换了。因此，由脑功能分析仪器30′执行的数据分析程序与实施例三不同。注意，尽管脑功能分析仪器30′采用了显示控制台模式，其中主要形成装置30A′、图像显示装置8和存储装置9在本实例中是组合在一起的，其也可能采用如下配置：图像显示装置8，或图像显示装置8和图像存储装置9从主要形成装置30A′中分离出来分别作为独立的图像显示单元和存储单元。

图19是图18所示脑功能分析仪器30′执行数据分析程序的流程图。本修饰例步骤S221中，数据分析装置36’首先使用SPM或相似技术对由数据预处理装置3预处理后的脑功能信息时间序列数据ρ(l，m，k，i)进行数据分析。接着，在步骤S231中，体素间联系度计算装置34以与实施例三中相似的方式计算等式(36)中的联系度矢量

接着，在步骤S241中，聚集装置300′采用与实施例三相似的方式根据由数据分析装置36′分析获得的脑功能数据ρ(l，m，k)对图17所示体素进行聚集。接着，在步骤S250中，图像生成装置7根据聚集的脑功能数据基于预设方法，生成颜色图像，如灰色或全彩图像。

此后，在步骤S260中，图像显示装置8显示上述生成的立体图像。

因此，本修饰例获得了与实施三相似的结果。

<实施例三修饰例2>

使用SPM或相似分析技术在对由数据预处理装置3预处理后的脑功能信息时间序列数据ρ(l，m，k，i)进行数据分析后，使用实施例三修饰例1中体素间联系度计算装置34计算出的联系度矢量

对分析后的脑功能数据ρ(l，m，k)进行聚集，但也可以考虑用分类分析技术作为数据分析方法。

因此本修饰例将描述使用决策树(分类分析技术的一种)作为数据分析技术的情况。比如，当使用决策树同时又使用解释变量(属性)作为脑功能信息时间序列数据ρ(l，m，k)，使用外部标准(任务)作为图4中任务和休息时间序列数据时，不必要的属性(体素)就可以被除去，同时只保留那些必需的属性(体素)。这样，通过使用联系度矢量

留下的体素就被聚集了。

比如，上述实例可以如下步骤进行：(1)选定某一特定体素，(2)对体素周围联系度不低于某一定值的体素进行聚集，(3)重复程序(2)一直到联接度低于该定值的体素全部被去除，(4)如果仍有其它体素存在，这些体素将被选出，然后回到处理程序(2)。

接着，在步骤S250中，图像生成装置7基于预设的方法生成上述聚集的脑功能数据

的颜色图像，如灰色或全彩色图像。

此后，在步骤S260中，图像显示方式8显示上述生成的立体图像。

因此，本修饰例获得了与实施三相似的结果。

[实施例四]

图20是实施例四中的脑功能分析仪器的结构示意图。实施例四的脑功能分析仪器40包括如下配置：体素间联系度计算装置44(而不是体素间联系度计算装置4)，数据估计值形成装置5，实施例一中脑功能分析仪器10的数据分析装置6，和数据测试装置400。与实施例一配置相同的，其符号也相同，下文将省略重复的内容。注意，尽管脑功能分析仪器40采用了显示控制台模式，其中主要形成装置40A、图像显示装置8和存储装置9在本实例中是组合在一起的，其也可能采取如下配置：图像显示装置8，或图像显示装置8和图像存储装置9从主要形成装置40A中分离出来分别作为独立的图像显示单元和存储单元。

体素间联系度计算装置44计算联系度矢量

通过对每个体素脑功能信息时间序列数据ρ(l，m，k，i)(其由数据预处理装置3作为等式(12)所示平均弥散度矢量

的函数从弥散张量数据获取装置2获得的弥散张量数据D(l，m，k)(等式(9))进行预处理)进行多种数据测试(比如，t测试、f测试、z测试)来调整参考值(测试值、重要性级别或类似数值)。

数据测试装置400根据体素间联系度计算装置44计算获得的联系度矢量

纠正参考值，在纠正过的参考值预处理后，再对每个体素脑功能数据ρ(l，m，k)进行数据测试。

图21是图20所示脑功能分析仪器40执行数据分析程序的流程图。

在步骤S300中，脑功能数据获取装置1获取由MRI设备50测得的脑功能信息原始时间序列数据ρ′(l，m，k，i)，同一步骤中，弥散张量信息获取装置2获取同样由MRI设备50测得的弥散张量数据D(l，m，k)。

接着，在步骤S310中，数据预处理装置3对在步骤S300中获得的脑功能信息原始时间序列数据ρ′(l，m，k，i)进行如实施例一所述的(a)到(d)型预处理，从而产生脑功能信息时间序列数据ρ(l，m，k，i)。

接着，在步骤S320中，通过体素间联系度计算装置44计算联系度矢量

$\stackrel{\&RightArrow;}{T}(l,m,k)=(T_l(l,m,k),T_m(l,m,k),T_k(l,m,k))\&CenterDot;\&CenterDot;\&CenterDot;\left(38\right)$

用于调整参考值(测试值、重要性级别或类似数值)。例如：对每一个体素的脑功能信息时间序列数据ρ(l，m，k，i)进行Z测试，该数据在步骤S310根据在该步骤获得的弥散张量数据作为等式(12)所示平均弥散度矢量的函数进行了预处理。

接着，在步骤S330中，通过数据测试装置400对每一个体素在步骤S310中预处理过的脑功能信息时间序列数据ρ(l，m，k，i)进行Z测试，进而计算一个Z分值。接着，同一步骤中，利用一个带高Z分值的体素(例如，以体素(l，m，k))作为参考点，参考点周围(26个方向)的一体素(如体素(l+1，m，k))带有的重要性级别值为a(l+1，m，k)，该值被降低为如等式(39)所示：

a′(l+1，m，k)＝{1-αT_i(l，m，k)}a(l+1，m，k)    …(39)

此处，α是一依赖于根据步骤S320中计算获得的联系度矢量

的加权系数。

接着，在步骤S340中，图像生成装置7对被在步骤S330中根据修正后的重要性级别被认定为在图像数据中具重要性的体素的脑功能数据ρ(l，m，k)生成一个值。

接着，在步骤S350中，图像显示装置8显示步骤S340中生成的立体图像数据。

注意，尽管测试重要性级别的值是根据上述弥散张量数据获取装置2获得的弥散张量数据D(l，m，k)计算得出等式(38)中的联系度矢量

来进行调整的，测试值也可以取代测试重要性级别使用该方法进行调整。

比如，在每一体素的任务图像和休息图像之间执行Z测试。任务图像是指执行某特定任务的图像(比如在图2中由粗线代表的图像1到3)，休眠图像是指未执行任务的图像(比如在图2中由粗线代表的图像5到7)。在这种情况下，体素(l，m，k)的Z分值(Z测试的测试值)通常由以下等式(40)代表：

$z(l,m,k)=\frac{|\mathrm{Mt}(l,m,k)-\mathrm{Mc}(l,m,k)|}{\sqrt{\mathrm{\&sigma;t}{(l,m,k)}^2+\mathrm{\&sigma;c}{(l,m,k)}^2}}\&CenterDot;\&CenterDot;\&CenterDot;\left(40\right)$

此处，Mt(l，m，k)代表任务图像的脑功能信息时间序列数据ρ(l，m，k，i)的平均值；Mc(l，m，k)代表休息图像的脑功能信息时间序列数据ρ(l，m，k，i)的平均值；σt(l，m，k)代表任务图像的脑功能信息时间序列数据ρ(l，m，k，i)的标准偏差，σc(l，m，k)代表休息图像的脑功能信息时间序列数据ρ(l，m，k，i)的标准偏差。在这一前提下，如果z＝0，可以说任务图像的脑功能信息时间序列数据ρ(l，m，k，i)的平均值与休息图像的脑功能信息时间序列数据ρ(l，m，k，i)的平均值没有区别。如果Z＝1，则意味着任务图像的脑功能信息时间序列数据ρ(l，m，k，i)的平均值与休息图像的脑功能信息时间序列数据ρ(l，m，k，i)的平均值由于时间序列数据的变化而变大了。从这一趋势中可见带高值Z分值的体素是与任务相关的激活区域。

正是因此，使用在步骤S320中计算获得的等式(38)的联系度矢量

等式(40)可改变成如下所示：

$z(l,m,k)=b(l,m,k)\frac{|\mathrm{Mt}(l,m,k)-\mathrm{Mc}(l,m,k)|}{\sqrt{\mathrm{\&sigma;t}{(l,m,k)}^2+\mathrm{\&sigma;c}{(l,m,k)}^2}},\&CenterDot;\&CenterDot;\&CenterDot;\left(41\right)$

b(l，m，k)在步骤S330中被修改为：

b′(l，m，k)＝{1+βT_l(l，m，k)}b(l，m，k)    …(42)

此处，β是一依赖于步骤S320通过等式(38)计算所得联系度矢量

的加权系数。

接着，在步骤S340中，图像生成装置7根据在所有体素的一般重要性级别基础上被认为是图像数据中重要的体素的脑功能数据ρ(l，m，k)生成一个值。然后，在步骤S350中，图像显示装置8显示步骤S340生成的立体图像数据。

如上所述，在实施例四中，在MRI设备获得的弥散张量的基础上，用于检测MRI设备获得的脑功能信息时间序列数据的参考值(测试值或重要性级别)在每个体素的基础上得到了调整。这样可以在考虑脑激活区域之间连接结构的同时，确定脑激活区域。

实施例四的修饰例

在实施例四中，当对每一体素的脑功能信息时间序列数据ρ(l，m，k，i)执行不同的数据测试(比如T测试、F测试或Z测试)时，数据测试的参考值(测试值、重要性级别或类似数值)被运用弥散张量数据获取装置2获得的弥散张量数据计算得出的联系度矢量

进行调整。为了实现本发明的目的，当对每一体素的脑功能信息时间序列数据ρ(l，m，k，i)使用不同的分析技术进行数据分析时，不论该分析技术是否具有统计学重要性，该数据检测(比如，T测试，F测试和Z测试)还是可以同样进行。

接着，在本修饰例中，我们考虑使用联系度矢量

调整参考值(测试值、重要性级别或类似数值)。

需要说明的是，为避免重复，仅对那些与实施例四不同的部分进行描述。

图22为实施例4修饰例中脑功能分析仪器构型示意图。本修饰例的脑功能分析仪器40′的配置不同如下：数据测试装置400′和数据分析装置46′取代了实施例四中的数据测试装置400。

数据分析装置46′为每一个体素对脑功能数据获取装置1获取的由数据预处理装置3预处理的脑功能信息时间序列数据ρ(l，m，k，i)进行数据分析，该分析使用如传统分类、回归和相关性的分析技术。

数据测试装置400′使用依据联系度矢量

的多种数据测试(比如，T测试，F测试和Z测试)来调整用于测试数据分析装置46′进行的分析技术是否具有统计学重要性的参考值。

图23是图22中脑功能分析仪器执行数据分析的流程图。

在步骤S360中，数据分析装置46′针对每一个体素的，对其由经脑功能数据获取装置1获取的经数据预处理装置3预处理过的脑功能信息时间序列数据ρ(l，m，k，i)进行数据分析，该分析使用如传统分类、回归和相关性的分析技术。

1.分类分析技术(判别分析、决策树、支撑矢量机或相类似技术)。

1-1：在判别分析的情况下，数据分类是通过多种判定函数来进行的，同时，使用解释性变量作为脑功能信息时间序列数据ρ(l，m，k，i)，外部标准作为任务和休息的时间序列数据。

1-2：在决策树的情况下，数据分类是用决策树来进行的，同时使用解释性变量(属性)作为任务和休息的时间序列数据，并通过分散脑功能信息时间序列数据ρ(l，m，k，i)形成几个类别，作为外部标准。

2.回归分析装置(使用多种函数(线性，非线性)进行回归分析)。

2-1：在线性回归分析情况下，数据分析是通过一个线性回归等式来进行的，同时使用解释性变量作为任务和休息，解释后变量作为脑功能信息时间序列数据ρ(l，m，k，i)。

3.相关性分析技术(一种利用估算模型(函数)和实际数据之间相关性的技术)。

数据分析通过利用预设的模型(比如，输入一种刺激则产生脑功能信息时间序列数据的函数)和脑功能信息的实际时间序列数据ρ(l，m，k，i)之间的相关性来进行的。

接着，从步骤S320到步骤S331(与实施例四相似)，数据测试装置400′对在步骤S360中获得的结果是否具有统计学重要性进行测试。同时，测试的参考值(测试值，重要性级别或类似数值)利用联系度矢量

(等式(38))进行校正。

例如，当对分类分析技术获得的结果进行f测试时，f测试将对在步骤S360中每一体素获得的结果进行测试，产生一个f值。接着，在同一步骤中，使用带高f值的体素(如体素(l，m，k))作为参考点，参考点周围(26个方向)的一体素(如体素(l+1，m，k))的重要性级别的a(l+1，m，k)值可根据步骤S230中计算出的联系度矢量

利用等式(39)来进行调整。

在回归分析技术的情况下，临界值(与重要性级别相对应)可以按上述方法调整，同样，在相关性分析技术的情况下，相关性的临界值可以按照上述方法进行调整。

接着，在步骤S340中，图像生成装置7根据在步骤S330中修正过的重要性级别而被认定为图像数据中重要体素的脑功能数据ρ(l，m，k)生成一个值。

接着，在步骤S350中，图像显示装置8立体显示步骤S340中生成的图像数据。

因此，本修饰例获得了与实施四相似的结果。

以上描述的实例仅仅是一些范例，本发明并不仅限于此。

比如，在每一个实例中，计算相邻体素联系度步骤可以紧随获取脑功能数据和弥散张量数据步骤。

另外，至于数据分析技术，使用子波分析、逻辑回归分析或相类似分析。

此外，也可以把实施例二或实施例三合并于实施例一。首先，当实施例二并入实施例一时，仅需引入一个配置，在该配置中，图1显示的脑功能分析仪器10另外带有与图10或13显示的数据平滑装置200或200′具有相同功能的数据平滑装置。接着，必要的是，在数据分析装置6确定一个回归系数以使实施例一步骤S50中的数据估计值Q最小化后，使用上述数据平滑装置对回归系数值进行各向异性地平滑处理。结果，生成图像的白质(神经纤维)的活跃度有可能比实施例一中获得的图像更突出。

同时，当实施例三并入实施例一时，仅需引入一个配置，其中，图1显示的脑功能分析仪器10进一步带有与图15或18显示的聚集装置300或300′具有相同功能的聚集装置。接着，必要的是，在数据分析装置6决定一个回归系数以使实施例一步骤S50中的数据估计值Q最小化后，使用上述聚集装置处理回归系数值以聚集各体素。结果，生成图像的白质(神经纤维)的活跃度有可能比实施例一中获得的图像更突出。

此外，从上述实例的结果可以看出，通过进一步使用弥散张量数据D(l，m，k)，也有可能扩展白质体素的激活区域。技术基本上与跟踪成像分析技术相似。

进一步讲，也可能以把上述技术(回归分析、相关性、分类、独立组份分析、测试或类似技术)仅仅应用于白质体素来捕获白质的激活。

同时，虽然我们在每个实施例的实验中都假定了整体设计，但是本发明也可以应用于与事件相关的fMRI。

另外，本发明适用于联系度分析。联系度分析包括多种技术，比如在SEM(结构等式模型)(协方差结构分析)的情况下，相关性和测试分析被应用在其计算中，但本发明在那种情况下也是适用的，平滑和聚集处理也同样适用。

神经元和神经纤维之间的连接信息可以从上述弥散张量数据中获得的，但如果连接信息(知识)可以从解剖结构上获知，那也是可以使用的。

除了上述描述的实例外，如果不偏离本发明的主题，其他各种实施方式也是可能的，本发明的范围仅限制于权利要求。

最后，利用本发明脑功能分析方法的脑功能分析效果将会被显示。图24显示了主体在声音刺激下执行简单重复任务时脑功能分析的结果。在实验中，重复其听到的来自耳机的声音被定义为任务，主体不思考或尽可能地不思考被定义为休息。任务的音量数(测量频率)设定为48，休息的音量数(测量频率)设定为60，未被用来分析的音量数(测量频率)设为12(参见图2)。另外，TR(重复时间与图2中音宽相对应)设定为5000毫秒。主体的数量为8人，在SPM和本发明的脑功能分析方法中都获得优异结果的主体的分析结果将在下文列出。

图24A和24B都显示了使用SPM的T测试的结果，其中图24A显示了T测试临界值被修正过的分析结果，图24B显示了T测试临界值没有被修正过的分析结果。图24C显示综合SPM和跟踪成像分析技术进行的分析结果，图24D显示了使用本发明实施例一的脑功能分析方法的分析结果。注意，位置补偿、标准化和平滑处理是作为SPM的预处理程序来进行的。滤波的半值宽度设定为9毫米。另外，dTV作为弥散张量分析的软件用于跟踪成像。

从解剖学上已知，Wernicke区(感官语言中枢)和Broca区(运动语言中枢)在简单重复过程中被激活，且二者由一束叫做弓形纤维束的神经纤维相连。该弓形纤维束在说话时被使用。

图24A显示，使用SPM时，当临界值被修正时，Wernicke区和Broca区的活动不能被检测到。图24B显示，当使用SPM且T测试的临界值未被修正时，二区域的活动可以被检测到，但连接二区域的弓形纤维束不能被检测到。与这些结果相比，如果利用结合SPM和跟踪成像的分析技术，如图24C所示，不仅韦尼克氏区(Wernicke Area)和卜洛柯区(Broca’s Area)还有弓形纤维束本身都可以被检测到。

同时，当使用本发明实施例一的脑功能分析方法时，除能够发现图24D中显示的二区域和连接二区域的弓形纤维外，还可以检测到弓形纤维束激活(在神经纤维中的脉冲传送形成弓形纤维束激活)。

此外，在结合SPM和跟踪成像的分析技术中，分析人员需事先确定神经纤维的一个跟踪起点和一个跟踪终点以用来跟踪神经纤维(在本例中为弓形纤维束)。因此，在神经纤维的跟踪起点和跟踪终点不能事先确定时(比如在某任务中，灰质的活动区域没有事先决定)，该技术将无法使用。与此相对应，根据本发明实施例一的脑功能分析方法，有可能检测出图24D显示的神经纤维的活性，尽管神经纤维的跟踪起点和跟踪终点不能事前确定。

工业应用

如上所述，根据本发明能够确定激活的脑功能区域，即使分析者并未事先确定神经纤维的起点和终点，仍有可能确定神经纤维(脑活动区域之间的连接结构)进而可能确定神经纤维的激活(神经纤维中的脉冲传输)，因此为医学诊断和治疗高度脑功能紊乱如痴呆症，失语证、精神分裂症或类似症状提供了极其有效的脑功能分析装置和脑功能分析程序。

Claims (18)

1.一种脑功能分析方法，其特征在于，包括：

第一步，从体素到体素的基础上，获取能够定位脑激活区域的脑功能数据和能够确定上述脑内质子弥散度的弥散张量数据，

第二步，在上述弥散张量数据的基础上形成相邻体素间联系度的估计值，

第三步，在上述相邻体素间联系度估计值的基础上分析上述脑功能数据。

2.如权利要求1所述的脑功能分析方法，其特征在于，上述第三步是一个回归分析，在该步骤中，将上述脑功能数据和任务的一者用作解释性变量，并将上述脑功能数据和上述任务的另一者用作解释后变量，计算估计值的最大值或最小值的一者，其中，将上述相邻体素间联系度估计值并入到上述脑功能数据每一体素的估计值中。

3.如权利要求1所述的脑功能分析方法，其特征在于，上述第三步是测试上述脑功能数据的方法，上述测试的参考值基于上述相邻体素间联系度进行调整。

4.如权利要求1所述的脑功能分析方法，其特征在于，上述第三步是对上述脑功能数据进行分类的技术，上述分类的参考值基于上述相邻体素间联系度估计值进行调整。

5.如权利要求1所述的脑功能分析方法，其特征在于，上述第三步是获取上述脑功能数据与预定模式的相关性的技术，上述相关性的参考值基于上述相邻体素间联系度估计值进行调整。

6.如权利要求1所述的脑功能分析方法，其特征在于，上述第三步是从上述脑功能数据中提取主要要素的技术，上述提取的参考值基于相邻体素间联系度估计值进行调整。

7.如权利要求1所述的脑功能分析方法，其特征在于，在上述第三步中，上述脑功能数据基于相邻体素间联系度估计值进行平滑。

8.如权利要求1所述的脑功能分析方法，其特征在于，在上述第三步中，上述脑功能数据基于相邻体素间联系度估计值进行聚集。

9.如权利要求1至8中的任一项所述的脑功能分析方法，其特征在于，还包括第四步，其进行下述的预处理，即，基于第一步中获得的弥散张量数据，将通过同一步骤中获得的脑功能数据确定的激活区域体素的脑功能数据的值传送到另一个体素的脑功能数据的值。

10.一种脑功能分析程序，其特征在于，使计算机作为下述装置而起作用，

脑功能数据获取装置，从体素到体素的基础上，获取能够确定脑激活区域的脑功能数据；

弥散张量数据获取装置，从体素到体素的基础上，获取能够确定上述脑中质子弥散度的弥散数据；

数据估计值形成装置，基于上述弥散张量数据形成相邻体素间联系度的估计值；和

数据分析装置，基于上述相邻体素间联系度估计值对脑功能数据进行分析。

11.如权利要求10所述的脑功能分析程序，其特征在于，上述数据分析装置是一个回归分析装置，在该步骤中，将上述脑功能数据和任务的一者用作解释性变量，并将上述脑功能数据和任务的另一者用作解释后变量，计算估计值的最大值或最小值的一者，其中，上述相邻体素间联系度估计值并入到上述脑功能数据每一体素的估计值中。

12.如权利要求10所述的脑功能分析程序，其特征在于，上述数据分析装置是用来检测上述脑功能数据的装置，上述检测的参考值基于上述相邻体素间联系度估计值进行调整。

13.如权利要求10所述的脑功能分析程序，其特征在于，上述数据分析装置是用来对上述脑功能数据进行分类的装置，上述分类的参考值基于上述相邻体素间联系度估计值进行调整。

14.如权利要求10所述的脑功能分析程序，其特征在于，上述脑功能分析装置是获取上述脑功能数据与预定模式的相关性的装置，上述相关性的参考值基于上述相邻体素间联系度估计值进行调整。

15.如权利要求10所述的脑功能分析程序，其特征在于，上述数据分析装置是从上述脑功能数据中提取主要要素的装置，上述提取的参考值基于上述相邻体素间联系度估计值进行调整。

16.如权利要求10所述的脑功能分析程序，其特征在于，上述数据分析基于上述相邻体素间联系度估计值对上述脑功能数据进行平滑。

17.如权利要求10所述的脑功能分析程序，其特征在于，上述数据分析装置基于上述相邻体素间联系度估计值对上述脑功能数据进行聚集。

18.如权利要求10至17中的任一项所述的脑功能分析程序，其特征在于，其使上述计算机作为数据预处理装置而起作用，基于第一步中获得的弥散张量数据，使来源于同一步骤中获得的脑功能数据的激活区域体素的脑功能数据的值，传送到另一个体素的脑功能数据的值。

Images