脑功能分析方法及脑功能分析程序的制作方法

专利名称：脑功能分析方法及脑功能分析程序的制作方法
技术领域：
，，，发明涉及一脑功能分、,方法,及脑功能分析程序，该脑功能分析方
歴RI (BOLD成像)，PET (正电子发射断层成像)或类似技术。
背景技术：
如今典型的非介入性脑功能检测方法有fMRI， PET， MEG(脑磁 图描技术)及其类似技术，其中歴RI被认为具有很强的空间分辨能力， 并被广泛应用。
fMRI是一种能通过反映脑不同物理量来确定脑激活区域的非常有 效的脑功能检测手段(参见非专利文献l， 2)。它除了能和检测脑结构 的MRI解剖图像一样反映活体器官的质子密度，纵向弛豫时间L和横 向弛豫时间T2以夕卜，还能特别的检测出脑激活区域中血液流量的增加。 现在已知，在脑激活区域的血液流量会增加，同时血液中血红蛋白的磁 性也会因其带氧状态即含氧血红蛋白(氧合血红蛋白)和脱氧血红蛋白 (去氧血红蛋白)而不同。人们认为fMRI在激活区域的信号(BOLD 信号)会增加，因为在增加的动脉血液中，具有干扰磁场能力的去氧血 红蛋白的量减少了。因此，我们可以使用fMRI，以脑工作时BOLD信 号的变化为线索，确定进行该工作的脑区域(激活区域)。
fMRI检测得到的BOLD信号的时间序列数据的分析技术包括基于 一种广泛的线性模型的SPM (统计参数绘图)(参见非专利文献1 )，或 是基本主要组份分析或独立组份分析的数据分析(参见非专利文献3 ), 或类似分析。这些分析技术的特点是其输出的结果中的BOLD信号的 时间序列数据的每个三维像素(体素)都各自经过统计处理，最后形成 图4象从而确定脑的激活区域。
然而在上述所说的分析技术中存在着以下问题在分析数据时，神经网络结构并未被考虑。脑中许多神经细胞通过神经突触构成复杂的网 络，而根据近些年来的脑研究，我们认为脑中每个神经区域通过互相协
作形成上述神经网络，进而作为一个整体发挥出更强的功能。例如在 脑完成一项工作时，其中的很多区域都出现激活现象。虽然上述所说的 分析技术可以用于分析该激活现象，确定激活区域，但是却很难确定激 活区域间的联系。
造成上述现象的原因有多个。因为fMRI是如上所述依据血液流量 检测血液信号的，fMRI可以捕捉到具有较高血液流量的脑灰质的活动 (神经细胞体)，但却很难捕捉到血液流量较低的脑白质的活动(神经 细胞轴突或神经纤维)。
同时为了捕捉作为神经网络结构基础的神经纤维组的走向， 一种 DTI (弥散张量成像)方法在近些年来越来越受到人们的关注。该方法 使用MRI检测身体器官内质子的弥散度，作为一种新的参数(参见非 专利文献4)。当时使用常规的解剖性MRI时，神经纤维在L加权图像 中是强信号，在T2加权图像中是弱信号。神经纤维在1\加权图象中变 成强信号是因为髓磷脂的存在。髓磷脂是由双分子层类脂物和大型蛋白 质组成，并沿着神经纤维走向成形。因此才出现了各向异性，即沿着神 经纤维走向上的质子的弥散常数很大，而其在垂直方向上的弥散常数很 小。DTI是使用MPG (弥散梯度磁场，)检测弥散的各向异性，从而突 出质子弥散情况的技术。该技术可以应用于分析ST ( Stejskal-Tanner) 脉冲序列检测获得的BOLD信号的时间序列数据的强度S' (1， m, k， i)，在SE (回波)脉冲的前后，STG ( Stejskal-Tanner Gradient梯度) 脉冲被加入用于检测弥散。
其中1， m， k为正离散变量，代表三维体素的位置，它们分别代 表体素在X, Y， Z方向上的值。另外，i为正离散变量，代表检测时间。 BOLD信号强度S' (1， m, k， i)表达式如下
<formula>formula see original document page 6</formula> (1)
当弥散加权图像形成后，式l中的弥散张量将成为分析数据。该弥
散张量如下所示<formula>formula see original document page 7</formula>
(2)
其中p' (1， m， k， i)代表当没有使用MPG时，BOLD信号强 度(普通脑功能数据分析对象)，而b为MPG强度参数。请注意p' (1， m， k， i)如下式所示
<formula>formula see original document page 7</formula> (3)
其中f(v)代表流动速度；TR为反复时间；TE为回波时间； (1， m， k， i)为质子密度。
近些年来，在分析fMRI测得的BOLD信号p' (1， m， k， i)的时 间序列数据后再使用弥散张量数据D (1， m， k)将目标脑区域连接起 来的脑分析方法已经被广泛使用(参见非专利文献5)。
非专利文献l: "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， W H Freeman & Co， 1997"。
非专利文献3: "Independent Component Analysis: Theory and Applications", T. W. Lee， Kluwer Acadmic， 1988"。
非专利文献4: "Korede wakaru diffusion MRI" S. Aoki， O. Abe， Syuujyun sha， 2002，，。
非专利文献5: "Combined functional MRI and tractography to demonstrate the connectivity of the human primary motor cortex in vivo", Guye M, et al" Neuroimage, Vol.19, pp.1349-1360, 2003，，。

发明内容
在上述的脑功能分析方法中，弥散张量D(l, m, k)并不是用于
7分析BOLD信号p ， ( 1， m， k, i)时间序列数据本身的。因此其有效 性值得怀疑。
本发明就上述问题提出了相应的解决方案，还提供了基于非介入检 测方法(如fMRI， PET或类似方法)测得的脑功能数据和能够确定 神经纤维组走向的弥散张量上，分析脑激活区域间连接结构的脑功能分 析方法和脑功能分析程序。
一方面，本发明提供了一种脑功能分析方法，该方法包括 一、在 从体素到体素基础上，获得能够确定脑激活区域的脑功能数据，和能够 确定脑中质子弥散程度的弥散张量。二、在弥散张量基础上形成相邻体 素间联系度估计值。
另一方面，本发明还提供了一种脑功能分析程序，该程序使计算机 获得可以确定脑激活区域的从体素到体素的脑功能数据；能够确定脑中 质子弥散度的从体素到体素的弥散张量；在弥散张量基础上形成相邻体 素间联系度估计值，并在该基础上分析脑数据。


图l是本发明实施例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显示了在与一体素相邻的两个体素均不被i^为是灰质的 区域中的数据预处理。
图IO是本发明实施例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测得的脑功能数据中原始时间序列 数据p'(l, m， k， i)(如等式(3));弥散张量数据获取装置2，用 于获得MRI设备50测得的弥散张量数据D(1， m, k)(如等式(2)); 数据预处理装置3，用于预处理脑功能数据获取装置1获取的脑功能信 息原始时间序列数据p'(1， m， k， i);体素间联系度计算装置4，用于 根据弥散张量数据获取法2获得的弥散张量数据D (1， m， k)来计算 代表相邻体素联系度的联系度矢量(如勿，附力)；数据估计值形成装 置5,用于从数据预处理装置3预处理后的脑功能信息时间序列数据p'
(1， m， k， i)和通过体素间联系度计算装置4得到的联系度矢量勿，附，" 中得出预定估计值Q;数据分析装置6,用于计算数据估计值Q的极值
(在本实施例中为最小值)；图像生成装置7，用于生成数据分析装置6计算得到的数据图像；图像显示装置8，用于显示图像生成装置7生成
的图像；存储装置9，其带有次储存，用于记录上述装置2到装置8中 获得、计算、形成、分析和生成的各种数据和主存储，用于储存电脑可 读的脑功能分析程序，用于执行脑功能分析的每一步骤(此后将详细描 述)。CPU (中央处理器)根据存储装置9中储存的脑功能分析程序控 制上述装置2至装置9的每一步骤，或是类似过程。以上所用符号的意 义在以下内容中将描述。
在此，MRI设备50是由以下部件组成 一磁性组件，该组件带有 能够产生静磁场，进而引起核磁共振的静磁场线围，和能够发出高频共 振波从而检测共振信号的高频线團，和能够产生梯度磁场，进而在共振 信号中加入位置信息编码的梯度磁场线團，或类似部件，和用于控制激 活这些线團的系统，或类似系统，其中，根据操作人员的要求，MRI 设备50产生不同的歴RI数据(与血液流量相关的脑功能信息的原始 时间序列数据p，(l， m， k， i)，与神经纤维走向相关的弥散张量D (1， m， k)，或类似数据)，并将这些数据传输给脑功能分析仪器10。请注 意，本发明并非一定需要通过MRI设备50获得脑功能信息的原始时间 序列数据p' (1， m， k， i)，也可以通过其它非介入性检测设备获得该 数据，如PET，通过它也可以获得多种有关脑功能的数据。
同时，对于图像显示装置8，有多种显示设备可用，例如阴极射 线管(CRT)显示器，TFT液晶显示器，等离子显示器，或类似i殳备， 或者是各种打印机，例如喷墨打印机，激光打印机，或类似设备。另 外，存储装置9由，例如RAM (随机存取存储器)，ROM (只读内存) 或类似设备组成。另外，存储装置9的主存储和次存储是相互独立分开 的，因此，主存储中的内容可以存放在例如光盘中，如硬盘，fl叩py软 盘(注册商标)，或是只读CD-ROMs,磁带，存储芯片或者类似设备。
本实施例中，虽然脑功能分析仪器10使用的显示装置包括的主要 形成装置IOA，图像显示装置8和存储装置9是组合在一起的，但是图 像显示装置8，或者图像显示装置8和存储装置9，也可以与主要形成 装置IOA分开来用，作为独立的图像显示单元和存储单元。在任何组合 方法中，脑功能分析仪器IO都是由电脑来完成。另外，CPU IOO从存 储装置9中读取脑功能分析程序，控制上述装置2到装置9中的每一个步骤。
这里所说的电脑是指一种能够根据预先设定的标准处理具有一定 结构的输入信息，然后构建和输出处理结果，例如普通的电脑，超级电 脑，主机，工作站，微型电脑，服务器，或类似设备都包括在其中。另 外，也可以使用一种系统(例如分布式计算机系统)，其中包括两台或
更多台电脑通过网络连接(如内部网，局域网，广域网)，或由其组 合组成的联系网络相互连接。
另外，为了增强对本发明的理解，图2至图7对本发明的脑功能分 析流程进行了概括性的描述。图2显示了典型的fMRI检测方法。图3 显示了典型的fMRI图像(二维切片图像)，图3A为位于床上的头，图 3B显示了组成头部二维横截面(二维图像)的体素。
例如，在某个实施例中，我们让研究对象完成一项任务(如让手 指轻轻拍打)，与这项任务相关的脑区域位置通过MRI检测被确定。 在图2所示的例子中， 一套程序在一次测量中被重复三次，该程序是由 一段时间T内的任务和相同长度时间T内的休息组成。这里，水平线 为时间轴，其中"ON"表示任务时间，"OFF，，表示休息时间。通常在一 个实验中，fMRI进行几十次检测。虽然在本实验中，fMRI进行了 24 次检测，但是任务开始和结束时的检测图像一般不被使用。这是考虑到 检测时间的滞后性。所以，这就意味着有效的fMRI检测实际只进行了 18次(相应的线为18条粗线)。ti ( i=l,2，...，24)表示各个实际检测 时间。以下内容中，ti将只以i表示(i=l，2，...，I)(分离的正整数)。
图3B是通过fMRI检测获得的二维切片图像Sk.这里的K表示实 验中获得的二维切片图像数量。在常规fMRI检测中，我们一般可以获 得20至30张二维切片图像(图片中，K^20至30)。虽然对图3B中的 二维切片图像Sk (k=l， 2... K)只是做了表面性的描述，但是它们却是由 三维立体像素，又称体素(voxel)所组成的。虽然由64x64体素组成的 二维切片图像作为一个整体显示在一张图片中，但是该二维切片图像也 可以被分为不同数量体素。虽然二维切片图像Sk中的任何一个体素都 可以通过两个分离的变量l和m表达为(1, m)，但是它也可以通过在 预定规则下对其设定一个数进而只用 一分离正变量j表示。从数学的角度，我们可以在预定规则下得到j和(l， m)--对应关系(对于(Sk)来说，
j三(l， m))。在本实例中，由于1的上限L和m的上限M都是64 (L=64,M=64)，所以导致了」'=1， 2， ...， 4096。在以下内容中，这 些表示方法都将被适当地使用。另外，图3B中的圆圏表示头部横截面 的轮廓。
图4表示了通过脑功能数据获取装置1获得的，并经过数据预处理 装置3预处理后的(将在后面描述)构成二维切片图像Sk的体素j的 脑功能信息时间序列数据p' (j， k， i)(Epk(j，i))和通过弥撒张量数据 获取装置2得到的D(j,k)(EDk(j))。如背景技术所述，当弥散加权图像 通过DTI方法生成后，也可同时从ST脉冲序列图像中获得二维切片图 像Sk的脑功能信息原始时间序列数据p'(l， m， k， i)(三p'k(l， m， i))和弥散 张量D(l， m, k)(三Dk (1， m))。该ST脉冲序列图像是在SE脉沖前后加入 STG脉冲用于检测弥散而获得的。pk(j, i)代表对原始时间序列数据p'k(l， m, i)进行先设预定预处理后(将在后面描述)获得的脑功能信息时间序 列数据。就是说，图中的pk(j, i)代表二维切片图像Sk中体素j在某一 时间点i上的脑功能信息。同时DkG)代表二维切片图像Sk中某一体素 j的弥散张量信息。pk(j，i)特指数量，Dk(j)在矩阵中由等式(2)代表。 他们可以是如等式(4)中所示的值，如
<formula>formula see original document page 13</formula>(4)
由于弥散张量是如本实例所示的典型的对称张量(非对角线元素具
有对称性)，其大体拥有六个部分。另外，图中的CLASS表示该任务是 否被完成，"ON (=1)，，和"OFF ^0)"分别代表"ON"和"OFF"，如图2 所示。
图5是图4表示的脑功能信息时间序列数据pk(j，i)的一个实例的表 格。虽然图5 ( A)中脑功能信息pk(j， i)的实际时间序列数据是连续数 值，但是一旦"4000"被设定为激活临界值之后，该数据可在预处理时被 进一步二元化，即体素值大于该临界值的处于激活"A，，状态，而小于该 临界值的处于非激活"I"状态，如图5 (B)所示。图6为图4中弥散张量Dk(j)和脑功能信息时间序列数据pk(j， i)的 立体构型图。其中图6A显示了置于床上的研究对象的头部。图6B显 示了通过图4脑功能信息时间序列数据和弥散张量生成k-th 二维切片 图像Sk的步骤。图6B中体素jl和j2上的五边形标记表示每个体素内 部空间(即体素的梯度场和张量场)，它们对应图4的脑信息时间序列 数据和弥散张量。
例如在背景技术所述的SPM中，可进行如下预处理
(a) 重排下面的图像与上面的二维切片图像排列一致(如，Sl)。 这样可以补偿检测过程中头部移动带来的位置偏差，去除相应的错误信 号。
(b) 空间标准化为了对多个研究对象收集数据并进行比较，我 们将每个研究对象的数据调整为Talairach标准脑。
(c) 平滑化Gaussian型过滤器被用于过滤原始带噪音的时间序 列数据。这样可以在不降低空间分辨率的情况下提高分析的灵敏度。
然后，对于通过MRI设备50获得的脑功能信息原始时间序列数据 p'k(l， m， i)，我们使用多种不同的检测方法检测其中的每一个体素。而 对于由多个研究对象组成的不同小组，我们在各组间进一步通过对象的 个人数据检测其与任务相联系的脑激活区域。
另外，在SPM中，我们事先建立了一个激活预测模型，该模型与 经过上述方法预处理后的功能数据的时间序列数据pk(j， i)(=pk(l， m， i)) 的相配程度需要看不使用弥散张量Dk(l， m)的一般线性模型所得的估计 参数。
同时，在本发明中，除了上述所说的预处理以外，我们还对使用通 过弥散张量获取装置2获得的弥散张量Dk(j) (=Dk(l， m))对表达脑功能 信息的原始时间序列数据p'k(l， m， i)进行(d )型预处理。该预处理方法 将在以下内容中详细描述。随后，对使用弥散张量Dk(j) (=Dk(l， m))预 处理的表达脑功能信息的原始时间序列数据p'k(l， m， i)进行数据处理， 该处理方法将在以下内容中详细描述。由此产生某图像显示目标的二维 图像ak(j)(-ak(l, m))。除了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" (nE(j, k))。请注意，使用单一变量"n，，和在X， Y， Z 轴上使用三个变量(j， 1， k)来代表同一个三维体素在以下内容中都将被 应用到(n三(l， m，k)E(j，k))。
以下我们将具体描述在以上假定条件下使用图l所示脑功能分析仪 器10进行的数据分析方法。图8是使用图1所示脑功能分析仪器10分 析数据的步骤流程图。
在步骤S10中，通过脑功能数据获取装置1获得MRI设备50检测 得到的表达脑功能信息的原始时间序列数据p，(l， m， k， i)，同时通过弥散 张量信息获取装置2获得MRI设备50检测得到的弥散张量D(l， m， k)。
下一步，通过数据预处理装置3进行各种预处理，例如位置补偿， 去除噪音，或如类似上述SPM(a)到(c)的预处理方法，用于处理在 上述S20中获得的表达脑信息的原始时间序列数据p'(l， m， k, i)。
在以上同一步骤中，为了对上述常规预处理后的表达脑功能信息时 间序列数据p'(l， m， k， i)作进一步(d)型预处理(数据转化为神经信 号)，可使用数据预处理装置3进行计算，例如使用以上获得的弥散张 量D(l, m， k)在符合弥散张量D(l, m， k)的本征等式(5)的本征矢量 (t(/，《a)， a=l，2，3)中计算最大本征值XM的相应值(^(z，w，")。
<formula>formula see original document page 16</formula>… (5)
为了简单起见，以下我们假设IXligi2gp31，最大本征值ll相应的 本征矢量(巧(z,附，")为4(/,附,A)。
另外，本征矢量&(z，"^)的基准应当通过l标准化。随后，在以上同一步骤中，数据预处理装置3，就某个体素(l,m,k) 而言，首先确定一对应方向矢量的相邻体素，该相邻体素与体素(l， m，k) 本征矢量^仏"^)的内积的绝对值在与体素(l， m， k)相邻的26个体素 (三维)中是最大的。
我们假设满足该条件的相邻体素是(1+1, m+l， k)。那么由体素(l, m， k)决定的该相邻体素的带有点对称关系的另一相邻体素就被确定了，在 此称为体素(l-l， m-l， k)，如图9所示。因为体素(1+1， m+l， k)和(1-1， m-l， k)与体素(l， m， k)本征矢量^仏"a)的内积是最大的，我们认为它们与体
素(l， m，k)的联系是最强的。
如果是如图9(C)所示的情况，那么与体素(l, m,k)相关的体素(l+l， m+l， k)和(1-1, m-l, k)可以被^人为是走方相同的神经纤维体素。
接下来在同一步骤中，数据预处理装置3根据以下三种不同情况改 写体素(l， m， k)的表达脑功能信息的时间序列数据p'(l， m， k, i)。
(情况1)根据预先设定的脑解剖学数据(通过MRI设备50获得 的加权图像Tl)(参见图9(A))，体素(1+1, m+l， k)和(l-l， m-l, k)都被认 为是灰质(神经元细胞体图片中表示为G)。
在这种情况下，体素(l，m,k)的表达脑功能信息时间序列数据p'(l， m, k， i)的数值将根据下式改写
<formula>formula see original document page 17</formula>… (6)
(情况2 )根据预先设定的脑解剖学数据(通过MRI设备50获得 的加权图像Tl)(参见图9(B))，体素(1+1， m+l， k)(或者体素(l-l, m誦l， k))被认为是灰质。
在这种情况下，体素(l,m,k)的表达脑功能信息时间序列数据p，(l， m, k， i)的数值将根据下式改写
<formula>formula see original document page 17</formula> (7)
(情况3 )根据预先设定的脑解剖学数据(通过MRI设备50获得 的加权图像Tl)(参见图9(C)),体素(1+1, m+l， k)和(l-l, m-l, k)都不被认为是灰质(即被认为白质(神经纤维图片中表示为W))。
在这种情况下，体素(l，m，k)的表达脑功能信息时间序列数据p'(l， m， k, i)的数值将根据下式改写
,,,.、//(/ + 1,附,+1,&，/) + //(/，附，^:，/) + / '(/ ——1，A，/) /0、 /w, A:") =-^-, … (A)
所有的体素都将经过这种处理。因此，被认为是灰质的体素(灰质
体素)中表达脑功能信息的时间序列数值可以如图9(A)和(B)下部图像 所示，该数据值可被反映(传输)给相邻的被认为是白质的体素(白质 体素)的表达脑功能信息的时间序列数值。另外，当被认为是白质的体 素如图9(C)中下部图像所示彼此相邻时，这些体素中反映脑功能的数将 被平滑处理。该预处理可以提高后续步骤的效果。
注意(d)型预处理并不局限于上述所说的实例。但是与灰质体素 相邻并从灰质体素中接受反映脑功能信息的时间序列数值的白质体素 可以按上述类似步骤进一步将数值传递给邻近的白质体素。此外，当灰 质体素中反映脑功能信息的时间序列数值通过上述步骤传递到白质体 素时， 一加权数值可根据弥散张量D(l， m， k)进行传递。另外，在上述 步骤中被改写的白质体素值和原始的白质体素值之间较大的值可以被 设定为该白质体素的值。
接下来，在步骤S30中，基于上述所说的弥散张量，
Z),J/,附，"Z)汰G，附，A)、 D(/, a)= Z),",(/，w,A:) Dmm(/，m，A:) 1^(/,m，"，… (9)
、D汰(/,/7a) a t(/，《a) a^(/，挑,"乂
使用体素间联系度计算装置4，计算两相邻体素联系度矢量，
5(/，m,A:)二(c,(/,w，A:), cm(/,w,W， q(/，w,A:)f … (10)
该联系度矢量可从下式所得
5(/，m,A:) = (l — )T + c^'(/, n,A:).… (11)
在此，等式(12)中的每个组份都可如等式(13)所定义，<formula>formula see original document page 19</formula>
在以下内容中，这被称为相邻体素间平均弥散度矢量。标记a为满 足不等式(14)的参数。T是满足等式(15)的常数矢量。<formula>formula see original document page 19</formula>
在等式(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轴方向 上的弥散度。<formula>formula see original document page 19</formula>
这里所说的X， Y和Z方向假定和图7显示的方向相一致。在等式 (9)中，弥散张量被假定为对称张量。
使用具有上述假设的弥散度矢量A^，W,通过等式(13)形成的等 式(11 )的相邻体素间联系度矢量5(^，"的每个组份Cl(l, m, k), Cm(l, m， k)和Ck(l, m, k)分别代表体素(l， m， k)和体素(l+l， m， k)的联系度，体素 (1， m， k)和体素(l， m+l， k)的联系度，体素(l， m， k)和体素(l， m， k+l)的联 系度。
需要注意的是在本实施例中，虽然在倾斜方向上相邻的体素(如体 素(l， m, k)和体素(1+1, m+l， k))的联系度没有被直接考虑，但是其已 经通过垂直方向上的体素(体素(1+1, m， k)或者体素(l, m+l， k))被间接 考虑了。显然，我们也可以设置成直接考虑倾斜方向上相邻的体素的联系度。
弥散度矢量力(/，附，"可以由以下等式计算 (I)当力(^，"被应用于椭圆形的模型中
<formula>formula see original document page 20</formula>
V (/, (/， m，" — Am (/,附，Ay 这里的IDI代表由等式(9)代表的矩阵的行列式。
(II)使用由等式(9)代表的矩阵的对角元素Dll(l， m， k)， Dmm(l， m, k)和Dkk(l， m, k)，可以得到以下等式
<formula>formula see original document page 20</formula>
(18)
(III)使用一代表弥散各向异性的指标，弥散各向异性典型代表形
式包括下式所示的部分各向异性
<formula>formula see original document page 20</formula>
和相对各向异性:
<formula>formula see original document page 20</formula>(20)
弥散度矢量力仏同样也可以被确定. 这里的DAV是(IV)另外，可根据预先设定的脑解剖数据(通过MRI设备50获 得的T1加权图像)确定脑白质体素和灰质体素，然后上述涉及的(I ) 的弥散度用于白质体素，而将其它弥散度矢量(如以上提到的(II)或 (III))用于灰质体素。
如上所述，作为本发明中的弥散度矢量力(^力，任何可以通过MRI 设备50获得的弥散张量D(l， m， k)确定神经纤维走向的矢量都可以被应 用为本发明的弥散矢量。
另外，除了与a参数相关的不等式(14)以外，a参数也可以被设 置为在白质体素和灰质体素间转换。比如，在灰质体素中拥有一个大的 值，而在白质体素中拥有一个小的值。
此外，除了等式(11)的联系度矢量勿，附，"以外，也可以使用以下 等式的联系度矢量
勿，m,A)" + a5'(/，附，A:)… (22)
另外，为了突出神经纤维的走向，可以提高等式(11)或者(22) 中联系度矢量^，w，"的每个组分cl(l， m, k)， cm(l, m, k)和ck(l, m， k)至 其N-th次方(N为2或2以上自然数)，从而得到一新的联系度矢量。
另外，通过将这些组分数值乘以a得到的值可能用于白质体素中， 通过将这些组分数值乘以b得到的值可能用于灰质体素中。此外，只有 这些组分中最大的组分才可以用上述方法突出。还有，只有拥有大于或 者等于(cl(l， m, k)+cm(1， m， k)+ck(1， m， k)}/3数值的组分才可以用上述 方法突出。
联系度矢量^^，"可以用一定方式形成以使灰质体素只与相邻的 六个体素相关联(六个方向)，而白质体素只与和本征等式(5)中最大 本征值对应的本征矢量的内积为最大值的体素以及与该体素拥有点对 称关系的体素(两个方向)相关联。之后，在步骤S40，为了根据步骤S20中预处理后的反映脑功能信 息的时间序列数据p(l， m, k， i)和步骤S40中计算得到的联系度矢量 勿，^a)，对图5(A)所示数据进行非参数型回归分析，如果存在六个相邻 体素，数据估计值形成装置5将形成如下数据估计值Q:
<formula>formula see original document page 22</formula> (23)
这里的标记i为时间，标记1， m和k分别为代表体素在X轴，Y轴， Z轴方向上位置的变量。请注意，实际中图7所示的立体构型体素中可 能有体素不代表脑(神经)，与其相关的l, m， k需被去除。同时，cp(i) 是图5(A)代表CLASS的变量的时间序列，它可以为l(-ON)或 O(-OFF)。在本实施例的非参数型回归分析中，当图5(A)所示的反映脑 功能信息时间序列数据p(l， m， k， i)和代表CLASS的变量的时间序列的 cp(i)为解释变量时，我们在这些变量间假设一个线性回归等式，其中
<formula>formula see original document page 22</formula>为回归系数
<formula>formula see original document page 22</formula>
这里的^和3分别表示预测结果，k是代表相邻体素间连续性的加权 参数。
虽然如上述所说的，等式(23)通过相邻体素间的连续性推算出了 估计值，但是该估计值也可以通过相邻体素间的平滑性推算出。
本实施例中，以非参数回归分析作为数据分析技术的基础在于，当 解释后的变量(p(i)与解释变量p(l， m， k， i)相比，q)(i)的数量是非常小的。 在图5( A)所示的实例中出现了 64x64x64=262144个解释变量p(l, m, k, i)，该数目与体素的数量相同，然而解释后变量cp(i)的检测值只有18个。 一般来说，解释后变量的检测值为10个，也可能超过100个。因为在这种情况下，无法进行常规的回归分析，本发明中我们使用了非参数回
归分才斤(非专利文献6 : "Spline Smoothing and Nonparametric Regression", R. L. Eubank, Marcel Dekker， Newyork， 1988，，)。
通常在非参数回归分析中，我们通过假定解释后变量(p(l， m， k， i) 的连续性来确定使估计值Q'(等式(25))最小化的勿)。
<formula>formula see original document page 23</formula> (25)
这里右边第 一项表示解释后变量q)(i)在每个时间i的预测结果勿)的 误差，第二项是用于估计解释后变量(p(i)在某个时间方向上的连续性。
然而如图5(A)所示，我们不能就解释后变量(()(i)假定连续性。因 此，在本实施例中等式(23 )表示的数据估计值Q是根据增强的非参数 回归分析得出的，这不是对解释后变量(p(i)而是对解释变量p(l， m， k， i) 设置的连续性约束(非专利文献7: "Nonparametric regression analysis of brain function images", H. Tukimoto et al" Institute of Electronics, Information and Communication Engineers D-II， Vol. J84， No.l2， pp2623-2633, December, 2001 )。
在此，我们将对等式(23)进行解释。右边第一项表示解释后变量 (p(i)在每个时间i与预测结果&)的误差，这与等式(25)相同。
第二到第四项表示相邻体素对解释后变量(p(i)影响的连续性，该连 续性可以通过相邻解释变量p(l， m， k， i)对解释后变量的影响进行估计， 也就是等式(24)的回归系数^^，"的连续性。
另夕卜，在本实施例中，由体素间联系度计算装置4计算得到的联系 度矢量5(/，附，"的每一个组份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"为
<formula>formula see original document page 24</formula>(26)
其中，n是代表体素位置的独立的正数变量，所有具有三维构型的 体素都被给予一预先指定的数值，从而被重新编排为一维，该变量与三 维显示变量(1， m, k)——对应。另外，为了简化目的，我们设置为 L=M=K，因此体素的数量为L3。同时，当所有的体素都以这种方法在 一维方向上重排后，系数C (n)代表体素n和相邻体素n'的联系度。该 系数是根据等式(11)的连接矢量^^，^特别确定的。
例如，当只限于图3(B)中的二维切片图像Sk时，与11=130体素相 邻的体素有四个n=127, 129， 131, 257。例如根据S30步计算的联系 度矢量勿，m，"，体素！!=127和体素n=130，体素n=129和体素n=130， 体素n=131和体素n=130，体素n=257和体素n=130间的联系度分别为 2， 2， 5, 5。
在这里，联系度是由
2讽127)-,0)}2， 2{5(129)-a(130)}2， 5讽131)-,0)}2， 5讽257)-S(130)}2所获得。
该实例中，通过使用最小平方法，回归系数^")可以被确定如下
<formula>formula see original document page 24</formula>(27)
其中，X为拥有p(n， i)组分的L3xl矩阵，^是拥有组分cp(i)的I维 矢量，S是拥有如下W")组分的L3维矢量，kO是k的优化值。k由交 叉确认方法确定，即k使等式(28)最小化。型的神经网络模型包括由输入层，中间层，和输出层组成的多层次神经 网络。
在本修饰例的步骤S40中，数据估计值形成装置5形成估计值，其 中基于体素间联系度计算装置4获得的联系度矢量勿，"^)所计算得到 的评估相邻体素间连续性的估计值被加入到平方误差值中，该平方误差 是神经网络模型中反映正常学习的估计值。
随后，在步骤S50中，使用误差后续增值法或类似方法，通过数据 分析装置6进行分析从而确定估计值的最小值，而后步骤S60中，图片 生成装置7生成输入(脑功能信息的时间序列数据)和输出的敏感度(即 输入值的改变对输出值的影响)，从而生成类似于实施例一的显示图像。 然后，在S70步中，图片显示装置8显示上述生成的立体图像。
这使我们有可能获得与实施例一相似的效果。
<实施例一修饰例2>
在步骤S40中，数据估计值形成装置5使用实施例一中的反映脑信 息的时间序列数据p(l，m，k，i)作为解释变量，cp(i)作为解释后变量，形 成等式(23)的估计值Q,其中，在本修饰例中，每个体素(l， m， k)都 有一个形成的估计值Q，"如下
g"=丄t {/ (/, m人0 -w人z.)}2
+丄t [c'仏w，"(々G+1附,'.)—々(人w， & 0}2
+ Cm (/, n, A:) 附+1, A:, Z) — /3(/，附人z')}2
+ Ci(/，w，it){》(/,w,A: + U)-y3(/,m,A:,f)}2] … (29)
该估计值的形成是通过使用反比例关系，即使用cp(i)作为解释变量， p(l, m， k, i)作为解释后变量。
其中给出了下列限制条件
》(/，m人/)二ZZJ^(/，附,A:)^(0 + 6 ... (30)
z=i m=i随后在步骤S50中，使用交叉证实方法，通过数据分析装置6确定 最佳值k，从而确定上述估计值q"'的最小值，使用概率性搜寻算法如 起源算法(Genetic Algoritm)，确定回归系数^，附，"。
随后在步骤S60中，根据体素的脑功能数据p(l， m, k),通过图像生 成装置7生成一个作为图像显示的数值，该数值通过不同分析所得结果 被确定为重要的。然后，在S70中，通过图像显示装置8显示以上生成 的立体图像。
这使我们有可能获得与实施例一相似的效果。
<实施例一修饰例3>
在步骤S40中，数据估计值形成装置5使用实施例一中的反映脑信 息的时间序列数据p(l，m，k, i)作为解释变量，(p(i)作为解释后变量，形 成等式(23)的估计值Q，其中，在本修饰例中，我们也可以使用cp(i) 作为解释变量，p(l，m，k，i)作为解释后变量，使用多层神经网络模型进 行非线性回归分析。请注意在本修饰例中，误差后续增值法不可以在 S50中用于确定极值，因此，数据分析装置6使用交叉证实法确定最佳 k值，通过概率性搜寻算法如起源算法进行数据分析。随后在步骤S60 中，根据体素的脑功能数据p(l， m， k)，通过图像生成装置7生成一个作 为图像显示的数值，该数值通过上述不同分析所得结果被确定为重要 的。然后，在S70中，通过图像显示装置8显示以上生成的立体图像。
这使我们有可能获得与实施例一相似的效果。
<实施例一<务饰例4>
实施例一S50中使用等式(24)的线性回归等式，通过数据分析装 置6进行非参数型回归分析，也可以使用分类分析技术，如判别分析 (^^专矛J文献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 )，支撑矢量机(非 专矛J文献 11 : "Introduction to support vector machine", Nello Cristianini and John Shawe画Taylor, Kyoritsu shuppan， 2005 )或者类似 技术作为分析技术。在任何情况下，使用反映脑功能信息时间序列数据p(l， m， k， i)作为 解释变量，使用任务和休息的时间序列数据作为外部标准，通过各种辨 别方程式，决策树或类似方程式对数据进行分类。在本修饰例步骤S40 中，通过数据估计值形成装置5形成估计值，其中于步骤S30中计算得 到的联系度矢量^，"^)被加入到每个分析技术的普通估计值中。
随后在步骤S50中，使用每个分析技术，通过数据分析装置6进行 数据分析。如果估计值出现上述改变，普通的解决方案将不适用，因此 可以使用交叉证实法，通过数据分析装置6确定优化K值，使用概率性 搜寻算法进行数据分析。
随后在S60中，根据体素的脑功能数据p(l， m， k)，通过图像生成装 置7生成一个用于图像显示数值，该数值通过上述不同分析结果被确定 是重要的。然后，在步骤S70中，通过图像显示装置8显示以上生成的 立体图像.
这使我们有可能获得与实施例一相似的效果。 <实施例一修饰例5>
当使用实施例一步骤S50中等式(24)的线性回归等式，通过数据 分析装置6进行非参数回归分析，独立组份分析也可以被用作一种数据 分才斤技术。(4^专利文献12: "Detailed explanation of Independent Component Analysis Novel world of signal analysis", Aapo Hyvarinen et al" Tokyo Denki University Press 2005 )。
在本修饰例中，数据估计值形成装置5在步骤S40中基于体素间联 系度计算装置4获得的联系度矢量勿，"^)形成估计值。
在本修饰例中，根据独立成份分析，分析反映脑功能信息时间序列 数据p(l, m， k, i)中的"空间独立性"条件被更改为"基于弥散张量信息的 联系度的空间独立性"。独立组份分析有几种特定技术，但是在使用交 叉关联度最小化的技术时，交叉关联度其实并不是被最小化了，而是在 联系度的基础上使其更接近估计值。结果，具有较大联系度的体素间的 交叉关联度变成了"独立程度"，而具有较小联系度(或没有关联度)的 体素间的交叉关联度被最小化，与联系度的估计值更相近，更确切的讲，例如在评估方程中，基于联系度的估计值的加权因子被设置于平方总 数的每一个值，该评估方程使关联矩阵中的非对角元素的平方总和最小 化。当联系度越大，加权因子越小。
在随后的步骤S50中，结合上述形成的估计值，通过分析装置6进 行独立组份分析。其中，通过使用交叉证实方法确定优化K值，同时通 过概率性搜寻算法如起源算法进行数据分析。
随后在步骤S60中，通过图像生成装置7生成上述获得的独立组分 显示图像。然而，由于上述分析产生了许多独立组分，因此需要从所有 的行矢量中提取与"任务/休息"关联值最高的行矢量。
在随后的步骤S70中，通过图像显示装置8显示上述形成的立体图像。
这使我们有可能获得与实施例一相似的效果。
请注意，除了主要的组分分析以外，作为从反映脑功能信息时间序 列数据中提取预设主要组份的提取技术的因素分析和定量III,但不是 上述的独立组份分析，也可以以上述相似的方式进行分析。
图20是实施例四中的脑功能分析仪器的结构示意图。实施例四的 脑功能分析仪器40包括如下配置体素间联系度计算装置44 (而不是 体素间联系度计算装置4)，数据估计值形成装置5，实施例一中脑功能 分析仪器10的数据分析装置6，和数据测试装置400。与实施例一配置 相同的，其符号也相同，下文将省略重复的内容。注意，尽管脑功能分 析仪器40采用了显示控制台模式，其中主要形成装置40A、图像显示 装置8和存储装置9在本实例中是组合在一起的，其也可能采取如下配 置图像显示装置8，或图像显示装置8和图像存储装置9从主要形成 装置40A中分离出来分别作为独立的图像显示单元和存储单元。
体素间联系度计算装置44计算联系度矢量^，^)，通过对每个体 素脑功能信息时间序列数据p(l， m， k， i)(其由数据预处理装置3作为等 式(12)所示平均弥散度矢量5'("w，"的函数从弥散张量数据获取装置 2获得的弥散张量数据D(l， m, k)(等式(9 ))进行预处理)进行多种 数据测试(比如，t测试、f测试、z测试)来调整参考值(测试值、重 要性级别或类似数值)。
数据测试装置400根据体素间联系度计算装置44计算获得的联系 度矢量勿，m，^)纠正参考值，在纠正过的参考值预处理后，再对每个体素
脑功能数据p(l， m, k)进行数据测试。
图21是图20所示脑功能分析仪器40执行数据分析程序的流程图。
在步骤S300中，脑功能数据获取装置1获取由MRI设备50测得 的脑功能信息原始时间序列数据p，(l，m,k，i),同一步骤中，弥散张量信 息获取装置2获取同样由MRI设备50测得的弥散张量数据D(l, m， k)。
接着，在步骤S310中，数据预处理装置3对在步骤S300中获得的 脑功能信息原始时间序列数据p'(l, m， k， i)进行如实施例一所述的(a) 到(d)型预处理，从而产生脑功能信息时间序列数据p(l，m，k，i)。
接着，在步骤S320中，通过体素间联系度计算装置44计算联系度矢量f(/,w,A:)，
<formula>formula see original document page 38</formula>(38)
用于调整参考值(测试值、重要性级别或类似数值)。例如对每 一个体素的脑功能信息时间序列数据p(i， m， k， i)进行Z测试，该数据在
步骤S310根据在该步骤获得的弥散张量数据作为等式(12)所示平均
弥散度矢量力'(^，^的函数进行了预处理。
接着，在步骤S330中，通过数据测试装置400对每一个体素在步 骤S310中预处理过的脑功能信息时间序列数据p(l， m， k， i)进行Z测试， 进而计算一个Z分值.接着，同一步骤中，利用一个带高Z分值的体 素(例如，以体素(l，m，k))作为参考点，参考点周围(26个方向)的 一体素(如体素(1+1， m， k))带有的重要性级别值为a(l+l， m， k)，该值 被降低为如等式(39)所示
<formula>formula see original document page 38</formula> (39)
此处，a是一依赖于根据步骤S320中计算获得的联系度矢量f (z，m，。 的加权系数。
接着，在步骤S340中，图像生成装置7对被在步骤S330中根据修 正后的重要性级别被认定为在图像数据中具重要性的体素的脑功能数 据p(l，m, k)生成一个值。
接着，在步骤S350中，图像显示装置8显示步骤S340中生成的立 体图像数据。
注意，尽管测试重要性级别的值是根据上述弥散张量数据获取装置 2获得的弥散张量数据D(l， m， k)计算得出等式(38)中的联系度矢量 f(/，m，^)来进行调整的，测试值也可以取代测试重要性级别使用该方法进
行调整。
比如，在每一体素的任务图像和休息图像之间执行Z测试。任务图 像是指执行某特定任务的图像(比如在图2中由粗线代表的图像1到3 )， 休眠图像是指未执行任务的图像(比如在图2中由粗线代表的图像5到7)。在这种情况下，体素(l， m, k)的Z分值(Z测试的测试值)通常由 以下等式(40)代表
此处，Mt(l， m， k)代表任务图像的脑功能信息时间序列数据p(l， m， k, i)的平均值；Mc(l，m，k)代表休息图像的脑功能信息时间序列数据p(l， m,k,i)的平均值；(yt(l，m,k)代表任务图像的脑功能信息时间序列数据 p(l，m，k,i)的标准偏差，oc(l，m，k)代表休息图像的脑功能信息时间序 列数据p(l，m，k, i)的标准偏差。在这一前提下，如果z一，可以说任务 图像的脑功能信息时间序列数据p(I， m， k， i)的平均值与休息图像的脑 功能信息时间序列数据p(l，m，k，i)的平均值没有区别。如果Z4，则意 味着任务图像的脑功能信息时间序列数据p(l, m， k， i)的平均值与休息 图像的脑功能信息时间序列数据p(l， m， k， i)的平均值由于时间序列数 据的变化而变大了。从这一趋势中可见带高值Z分值的体素是与任务相 关的激活区域。
正是因此，使用在步骤S320中计算获得的等式(38)的联系度矢 量JU^a;i，等式(40)可改变成如下所示
b(l， m, k)在步骤S330中被修改为
6'(/,附,/1) = {1 + / 7;(/，附,&)}6(/，附,^0… (42)
此处，p是一依赖于步骤S320通过等式(38)计算所得联系度矢
量f(/，^a)的加权系数。
接着，在步骤S340中，图像生成装置7根据在所有体素的一般重 要性级别基础上被认为是图像数据中重要的体素的脑功能数据p(l， m， k) 生成一个值。然后，在步骤S350中，图像显示装置8显示步骤S340生 成的立体图像数据。
如上所述，在实施例四中，在MRI设备获得的弥散张量的基础上，用于检测MRI设备获得的脑功能信息时间序列数据的参考值(测试值 或重要性级别)在每个体素的基础上得到了调整。这样可以在考虑脑激 活区域之间连接结构的同时，确定脑激活区域。
实施例四的修饰例
在实施例四中，当对每一体素的脑功能信息时间序列数据p(l，m，k, i)执行不同的数据测试(比如T测试、F测试或Z测试)时，数据测试 的参考值(测试值、重要性级别或类似数值)被运用弥散张量数据获取 装置2获得的弥散张量数据计算得出的联系度矢量勿，"^)进行调整。为 了实现本发明的目的，当对每一体素的脑功能信息时间序列数据p(l， m， k, i)使用不同的分析技术进行数据分析时，不论该分析技术是否具有统 计学重要性，该数据检测(比如，T测试，F测试和Z测试)还是可以 同样进行。
接着，在本修饰例中，我们考虑使用联系度矢量^，m，"调整参考 值(测试值、重要性级别或类似数值)。
需要说明的是，为避免重复，仅对那些与实施例四不同的部分进行 描述。
图22为实施例4修饰例中脑功能分析仪器构型示意图。本修饰例 的脑功能分析仪器40'的配置不同如下数据测试装置400'和数据分析 装置46'取代了实施例四中的数据测试装置400。
数据分析装置46'为每一个体素对脑功能数据获取装置1获取的由 数据预处理装置3预处理的脑功能信息时间序列数据p(l， m， k， i)进行数 据分析，该分析使用如传统分类、回归和相关性的分析技术。
数据测试装置400，使用依据联系度矢量勿，附，"的多种数据测试(比 如，T测试，F测试和Z测试)来调整用于测试数据分析装置46'进行 的分析技术是否具有统计学重要性的参考值。
图23是图22中脑功能分析仪器执行数据分析的流程图。
在步骤S360中，数据分析装置46，针对每一个体素的，对其由经脑 功能数据获取装置1获取的经数据预处理装置3预处理过的脑功能信息时间序列数据p(l, m， k， i)进行数据分析，该分析使用如传统分类、回归 和相关性的分析技术。
1. 分类分析技术(判别分析、决策树、支撑矢量机或相类似技术)。
1-1:在判别分析的情况下，数据分类是通过多种判定函数来进行 的，同时，使用解释性变量作为脑功能信息时间序列数据p(l，m，k，i)， 外部标准作为任务和休息的时间序列数据。
1- 2:在决策树的情况下，数据分类是用决策树来进行的，同时使 用解释性变量(属性)作为任务和休息的时间序列数据，并通过分散脑
功能信息时间序列数据p(l， m， k， i)形成几个类别，作为外部标准。
2. 回归分析装置(使用多种函数(线性，非线性)进行回归分析)。
2- 1:在线性回归分析情况下，数据分析是通过一个线性回归等式 来进行的，同时使用解释性变量作为任务和休息，解释后变量作为脑功 能信息时间序列数据p(l， m， k， i)。
3. 相关性分析技术(一种利用估算模型(函数)和实际数据之间 相关性的技术)。
数据分析通过利用预设的模型(比如，输入一种刺激则产生脑功能 信息时间序列数据的函数)和脑功能信息的实际时间序列数据p(l， m, k， i)之间的相关性来进行的。
接着，从步骤S320到步骤S331 (与实施例四相似)，数据测试装
置400，对在步骤S360中获得的结果是否具有统计学重要性进行测试。
同时，测试的参考值(测试值，重要性级别或类似数值)利用联系度矢 量f(/，m，A)(等式(38))进行校正。
例如，当对分类分析技术获得的结果进行f测试时，f测试将对在 步骤S360中每一体素获得的结果进行测试，产生一个f值。接着，在 同一步骤中，使用带高f值的体素(如体素(l,m，k))作为参考点，参考 点周围(26个方向)的一体素(如体素(1+1, m， k))的重要性级别的a(l+l， m，k)值可根据步骤S230中计算出的联系度矢量勿，附，。，利用等式(39) 来进行调整。在回归分析技术的情况下，临界值(与重要性级别相对应)可以按 上述方法调整，同样，在相关性分析技术的情况下，相关性的临界值可 以按照上述方法进行调整。
接着，在步骤S340中，图像生成装置7根据在步骤S330中修正过 的重要性级别而被认定为图像数据中重要体素的脑功能数据p(l， m， k) 生成一个值。
接着，在步骤S350中，图像显示装置8立体显示步骤S340中生成 的图像数据。
因此，本4务饰例获得了与实施四相似的结果。
以上描述的实例仅仅是一些范例，本发明并不仅限于此。
比如，在每一个实例中，计算相邻体素联系度步骤可以紧随获取脑 功能数据和弥散张量数据步骤。
另外，至于数据分析技术，使用子波分析、逻辑回归分析或相类似 分析。
此外，也可以把实施例二或实施例三合并于实施例一。首先，当实 施例二并入实施例一时，仅需引入一个配置，在该配置中，图l显示的 脑功能分析仪器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显示，当4吏用SPM且T测试的 临界值未被修正时，二区域的活动可以被检测到，但连接二区域的弓形 纤维束不能被检测到.与这些结果相比，如果利用结合SPM和跟踪成 4象的分析技术，如图24C所示，不仅韦尼克氏区(Wernicke Area)和 卜洛柯区(Broca，sArea)还有弓形纤维束本身都可以被检测到。
同时，当使用本发明实施例一的脑功能分析方法时，除能够发现图 24D中显示的二区域和连接二区域的弓形纤维外，还可以检测到弓形纤 维束激活(在神经纤维中的脉冲传送形成弓形纤维束激活)。
此外，在结合SPM和跟踪成像的分析技术中，分析人员需事先确 定神经纤维的一个跟踪起点和一个跟踪终点以用来跟踪神经纤维(在本 例中为弓形纤维束)。因此，在神经纤维的跟踪起点和跟踪终点不能事 先确定时(比如在某任务中，灰质的活动区域没有事先决定)，该技术 将无法使用。与此相对应，根据本发明实施例一的脑功能分析方法，有 可能检测出图24D显示的神经纤维的活性，尽管神经纤维的跟踪起点和 跟踪终点不能事前确定。
工业应用
如上所述，根据本发明能够确定激活的脑功能区域，即使分析者并 未事先确定神经纤维的起点和终点，仍有可能确定神经纤维(脑活动区 域之间的连接结构)进而可能确定神经纤维的激活(神经纤维中的脉冲 传输)，因此为医学诊断和治疗高度脑功能紊乱如痴呆症，失语证、精 神分裂症或类似症状提供了极其有效的脑功能分析装置和脑功能分析 程序。
权利要求
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中的任一项所述的脑功能分析程序，其特 征在于，其使上述计算机作为数据预处理装置而起作用，基于第一步中 获得的弥散张量数据，使来源于同一步骤中获得的脑功能数据的激活区 域体素的脑功能数据的值，传送到另一个体素的脑功能数据的值。
全文摘要
本发明提供一种方法包括获取脑功能数据和弥散张量数据(S10)，基于该弥散张量数据计算相邻体素间的联系度(S30)，基于该脑功能数据和相邻体素间联系度来形成数据估计值(S40)，对该数据估计值进行非参数的回归分析(S50)，基于分析结果形成和显示图像(S60，S70)。
文档编号A61B5/055GK101287410SQ20068003806
公开日2008年10月15日 申请日期2006年10月6日 优先权日2005年10月12日
发明者月本洋 申请人:学校法人东京电机大学