CN109416939A - 面向群体应用的经颅脑图谱生成方法、预测方法及其装置 - Google Patents

Abstract

一种经颅脑图谱生成方法和一种面向群体应用的经颅脑图谱预测方法及相应的经颅脑图谱预测装置。该经颅脑图谱生成方法包括如下步骤：在个体层面创建颅面坐标系；建立用于连接颅骨位置与大脑位置的经颅映射系统；使用马尔可夫链中的两步随机过程构建经颅脑图谱。经颅脑图谱将不可见的脑内图谱标号信息投射到可见的头皮上来，使得研究者或者医生可以使用这些脑结构和功能图谱信息，大大提升了脑图谱在经颅脑映射技术使用中的作用。

Description

面向群体应用的经颅脑图谱生成方法、预测方法及其装置

技术领域

本发明涉及一种经颅脑图谱生成方法，同时涉及一种面向群体应用的经颅脑图谱预测方法，还涉及相应的经颅脑图谱预测装置，属于认知神经科学技术领域。

背景技术

当前，认知神经科学所面临的一个重要课题是建立脑功能和脑结构之间的对应关系。以功能磁共振成像(fMRI)、功能近红外光谱(fNIRS)为代表的功能脑成像技术，可以使神经科学的研究者们以非侵入的方式对活体人脑的功能进行观测。

在现有的功能脑成像技术中，往往通过放置在颅骨表面的经颅成像装置，对颅内大脑的活动进行观测和干预，因此也称为经颅脑成像技术。经颅脑成像技术中存在两个空间的概念，其中一个是对经颅成像装置可见的颅骨表面空间，另一个是对经颅成像装置不可见的颅内脑空间。这两个空间的分离造成了经颅成像技术在功能数据定位和经颅成像装置正确放置两个方面的问题。解决这些问题的关键在于建立两个空间的对应关系。

经颅脑成像技术本身只能提供脑功能信息，而不能提供脑结构信息，并且只能将得到的脑功能信息定位到经颅成像装置在颅骨表面空间，而非人脑所在的颅内脑空间。这使得经颅脑成像技术无法在传统脑成像的定位框架下，实现经颅脑成像数据到标准脑空间(即MNI空间)的定位。这也意味着基于标准脑空间和脑图谱建立的传统脑成像研究的参照系统不能直接适用于经颅脑成像技术。

脑图谱(Brain Atlas，简称为BA)是脑成像研究中重要的标准参照系统。首先，脑图谱为脑成像研究提供了一个标准平台，使得基于不同被试大脑展开的不同研究的结果能够相互比较和印证，从而得到对大脑功能架构的全面认识。其次，脑图谱中提供的关于人脑的先验知识，是脑成像研究中划定感兴趣区域，划定脑网络节点的依据，也是解释和报告脑成像结果的依据。因此，将脑功能数据定位到脑图谱所在的标准脑空间，是脑成像研究分析的必要环节。但是，现有经颅脑成像技术并没有提供一个具备科学根据、并且对不同人群具备可预测性的通用脑图谱模型。实践中，可见的操作空间(即头皮表面，尤其是头皮上表面)和不可见的效用空间(即颅内脑空间)之间的分离仍然是有效应用经颅脑映射技术的最大挑战之一。

发明内容

本发明所要解决的首要技术问题在于提供一种经颅脑图谱生成方法。通过经颅脑图谱，可以将不可见的大脑内的图谱信息投射到可见的头皮表面上，使原本分离的操作空间和效用空间融合在一起。

本发明所要解决的另一个技术问题在于提供一种面向群体应用的经颅脑图谱预测方法。

本发明所要解决的又一个技术问题在于提供一种面向群体应用的经颅脑图谱预测装置。

为了实现上述目的，本发明采用下述的技术方案：

根据本发明实施例的第一方面，提供一种经颅脑图谱生成方法，包括如下步骤：

(1)在个体层面创建颅面坐标系；

(2)建立用于连接颅骨位置与大脑位置的经颅映射系统；

(3)使用马尔可夫链中的两步随机过程构建经颅脑图谱。

其中较优地，所述步骤(1)包括如下子步骤：

(11)在头皮表面识别五个颅骨标记Nz、Iz、AL、AR和Cz；

(12)将头皮表面和通过Nz、Cz和Iz的平面之间的相交曲线定义为颅骨赤道；

(13)给出头皮表面上的点p，经度曲线可以唯一地确定为头皮表面和通过AL、AR和p的平面之间的相交曲线，p'是颅骨赤道与经度曲线之间的交叉点；

(14)上头皮的任意点p由一对非负实数(p_e,p_l)唯一确定：

p_e＝L_NZ-p’/L_e,p_e∈[01]

p_l＝L_AL-p/L_AL-p-AR,p_l∈[01]

其中，L_Nz-p'是沿着颅骨赤道的Nz到p'的曲线长度，L_e是颅骨赤道的全长；L_AL-p是沿着全长为L_AL-p-AR的经度曲线从AL到p的曲线长度。

其中较优地，所述步骤(1)还包括步骤(15)：在标准半球上建立CPC空间；使用Hammer-Aitoff投影对带有CPC坐标的半球进行平面化，生成在扁平椭圆上呈现的CPC坐标系的地图。

其中较优地，所述步骤(2)包括如下子步骤：使用气球膨胀模型在个体空间中确定对应于给定的头皮表面任意点p的底层皮层位置c；

在所有皮层位置被空间标准化为MNI空间之后，聚合所有(p，c)对，生成确定性的个体经颅脑映射模型。

其中较优地，所述步骤(2)还包括如下步骤：

集成所有个体模型来生成群体层面的概率经颅脑映射模型：P(c|p)；

其中，p(p_e,p_l)∈CPC，c(x,y,z)∈C，C是MNI空间的子集。

其中较优地，所述步骤(2)中，所述概率经颅脑映射模型由下列公式生成：

其中，p_i是CPC₁₀₀网格位置，i＝1，2，...，Np；c_j是通过经颅脑映射技术获得的MNI空间中的皮质体素，j＝1，2，...，Nc；#是个体总数。

其中较优地，所述步骤(3)还包括如下步骤：

(31)给定点p(p_e，p_l)通过概率经颅映射P(c|p)被映射到MNI空间中的皮层位置c(x，y，z)；

(32)皮层位置c(x，y，z)被映射到标号空间L中的标号l。

其中较优地，在步骤(3)中，假设皮层域c上的点是脑中域b的点的子集，如果p和c离散化，P(l|p)由Chapman-Kolmogorov方程式表示：

其中，p_i是CPC空间中的离散位置，i＝1，2，...，Np；c_j是C在MNI空间中的离散位置，j＝1，2，...，Nc；l_k是特定脑图谱的标号，k＝1，2，...，N₁。

其中较优地，所述步骤(32)中使用LPBA40脑图谱、AAL标号图谱或Talairach图谱中的任意一种。

其中较优地，所述经颅脑图谱生成方法还包括步骤(4)：生成最大似然标号图和/或最大概率图。

根据本发明实施例的第二方面，提供一种面向群体应用的经颅脑图谱预测方法，包括如下步骤：

针对群体中的单独个体，当经颅设备的探头在给定坐标的头皮表面的任意位置进行刺激或记录时，通过上述的经颅脑图谱给出每个靶向脑区域被探及的概率。

其中较优地，所述经颅设备为经颅脑成像装置或经颅脑治疗装置中的任意一种。

根据本发明实施例的第三方面，提供一种面向群体应用的经颅脑图谱预测装置，用于实施上述的经颅脑图谱预测方法，其中：

所述经颅脑图谱预测装置为头盔或者头套形状，其上表面覆盖有上述的经颅脑图谱，下表面在使用时紧贴用户的头皮表面。

其中较优地，在所述经颅脑图谱预测装置的表面分布有若干小孔，用于供经颅设备的探头通过所述小孔接触用户的头皮表面。

其中较优地，所述小孔根据经颅脑图谱图像中的不同子区域排列成点阵形状，其排列密度根据所述子区域的面积呈反向变化。

本发明所提供的经颅脑图谱将不可见的脑内图谱标号信息投射到可见的头皮上来，使得研究者或者医生可以“直接”使用这些脑结构和功能图谱信息，大大提升了脑图谱在经颅脑映射技术使用中的作用。此外，经颅脑图谱将原本分离的操作空间和效用空间融合在一起，使得操作者就像在脑空间上操作一样，这将使得经颅脑映射技术中的定位等固有难题得到更加彻底的解决。

附图说明

图1(a)～图1(e)是经颅脑图谱的系列实施例示意图；

图2是从磁共振图像中识别颅骨标记的示意图；

图3是CPC坐标系的示意图；

图4给出了单个CPC坐标点的概率经颅脑映射示意图；

图5是TBA_LPBA的示意图；

图6是群体层面TBA_LPBA对个体的预测性结果示意图；

图7为本发明所提供的经颅脑图谱预测装置的一个实施例示意图。

具体实施方式

下面结合附图和具体实施例对本发明的技术内容做进一步的详细说明。

前已述及，对经颅成像装置可见的颅骨表面空间和对经颅成像装置不可见的颅内脑空间的分离是有效应用经颅脑映射技术的最大挑战之一。为了解决两个空间的对应问题，本发明首先提出了经颅脑图谱(Transcranial Brain Atlas，简称为TBA)的概念。经颅脑图谱是建立在头皮表面上的脑图谱，它将不可见的大脑内的图谱信息投射到可见的头皮表面(尤其是头皮上表面)上，使研究者或者医生可以直接使用这些涉及脑结构和脑功能的图谱信息。

具体地说，本发明首先显式地构建一个标准的颅骨坐标系统，用来实现对不同个体颅骨表面空间的定量描述。然后，依据颅-脑对应关系在人群上一致性的假设，建立一个标准颅骨表面空间到脑图谱所在的标准脑空间之间的对应关系。最后，本发明从颅骨表面空间-标准脑空间，以及脑图谱提供的标准脑空间-脑分区标号空间这两部分的对应关系中，求解出颅骨表面空间-脑分区标号空间之间的对应关系。作为结果，本发明将对应的脑标准空间和脑图谱中的信息反向呈现到颅骨坐标系统，从而形成一种新型的“经颅脑图谱”。经颅脑图谱的一个重要性质是能够仅仅通过颅骨位置信息，直接推断对应脑分区标号的信息，因此可以粗略地理解为一个建立在头皮表面上的脑图谱。经颅脑图谱本质上是一个建立在坐标化的头皮表面上的大脑功能地图，即在一个坐标化的大脑空间中，让传统脑图谱对应每个脑部位置及其功能或解剖学标号，从而描绘出经颅脑映射技术可及的各个皮层位置及其对应的图谱标号，并将它们明确呈现在作为可见操作空间的头皮表面上。经颅脑图谱可以将传统脑图谱中的先验脑分区信息，在人群水平颅-脑对应关系意义下映射到用于经颅脑成像装置放置的颅骨空间中。因此，可以看作是传统脑图谱在经颅成像技术领域的扩展。在经颅脑图谱的框架下，经颅数据在脑空间的定位可以转化为对经颅成像装置在颅骨空间的定位，这使得经颅脑映射技术的实时定位成为可能。此外，由于经颅脑图谱的标号信息是在颅骨空间显示的，这个特性非常有利于经颅脑图谱叠加到个体头皮表面进行显示，从而以直观的方式指导经颅成像装置在被试颅骨表面的放置。因此，经颅脑图谱的建立，解决了经颅脑映射技术中操作空间与效用空间分离的矛盾。

图1(a)～图1(e)是经颅脑图谱的系列实施例示意图。如图1(a)所示，假设在头皮表面上定义了统一的颅面坐标系，它描述了可以设置经颅脑映射技术的探头的所有可能的位置，并且将探头设置在给定的颅骨位置(图中黑点所在位置)。对于特定的个体(例如图1(b)中的受试者sub 1)，设置在其头皮表面的探头将能够探及特定皮层位置/脑区域(图中黄点所在位置)。但在群体层面，考虑到个体之间的解剖结构的差异，这种颅-脑对应关系可能不是确定的。如图1(c)所示，在标准化到标准脑空间之后，我们可以捕获这种概率对应关系的空间分布(即黑色圆圈内的彩色区域)。通过给定从大脑图谱所获得的解剖信息，这种概率对应关系能够提供如何从颅骨位置访问每个大脑区域的群体层面的概率，并且将其作为先验知识。经颅脑图谱本质上是为颅面坐标系定义的整个大脑空间映射该先验知识。具体在图1所示的系列实施例中，如果仅针对每个颅骨位置考虑最可能探及的脑分区标号及其关联的概率，则如图1(d)所示的最大似然标号图(MLM)和如图1(e)所示的最大概率图(MPM)可以作为经颅脑映射技术的探头排列的有用指导。

在本发明的实施例中，经颅脑图谱的生成过程主要包括三个步骤：1)在个体层面创建颅面坐标系；2)建立用于连接颅骨位置与大脑位置的经颅映射系统；3)使用马尔可夫链中的两步随机过程构建经颅脑图谱。

下面对经颅脑图谱的详细生成过程进行说明。

1)在个体层面创建颅面坐标系(Cranial surface coordinate system)

颅面坐标系需要满足两个基本要求：第一，它应为个体头皮表面提供一对一的映射；第二，为了便于进行群体层面的研究，对于颅面坐标系中的每个位置，应该使来自不同个体的下层皮层位置在神经解剖学上基本一致。

CPC坐标系的基本思想是在头表面上构造一个类似于经纬线的坐标系统。与地理经纬线系统不同的是，CPC坐标系对“经纬度”的确定是利用两次表面比例测量的方式确定的。在本发明的实施例中，通过如下三个步骤在个体的头皮表面建立被称为连续比例坐标空间(简称为CPC空间)的二维比例坐标系：

11)在个体空间(navie space)的头皮表面识别至少五个来源于10-20系统的颅骨标记Nz、Iz、AL、AR和Cz(参见图3(a))。在磁共振图像中识别颅骨标记的示例参见图2，其中，Iz是人类头骨的外枕骨隆起，斜方肌附着于其上；AL和AR是左耳和右耳前点，被识别为耳屏的峰值区域；Nz被确定为鼻梁上部的凹陷位置；Cz被确定为颅表面测地线AL-Cz-AR和Nz-Cz-Iz的交点，并分别平分两条颅表面测地线。

12)将颅骨赤道被定义为头皮表面和通过Nz、Cz和Iz的平面之间的相交曲线(即图3(a)中的曲线1)；

13)给出头皮表面上的点p，经度曲线(即图3(a)中的曲线2)可以唯一地确定为头皮表面和通过AL、AR和p的平面之间的相交曲线，p'是颅骨赤道与经度曲线之间的交叉点。

在以上三步定义的基础上，上头皮(高于由Nz、Iz、AL、AR点指定的曲线)的任意点p都可以由一对非负实数(p_e,p_l)唯一确定：

p_e＝L_Nz-p’/L_e,p_e∈[01] (1)

p_l＝L_AL-p/L_AL-p-AR,p_l∈[01] (2)

其中，L_Nz-p'是沿着颅骨赤道的Nz到p'的曲线长度，L_e是颅骨赤道的全长(从Nz到Iz)；L_AL-p是沿着全长为L_AL-p-AR的经度曲线从AL到p的曲线长度。如图3(b)所示，作为任意点的p点的表面位置由p'和p分别占这两条曲线的比例唯一地表示，计算公式参见公式(1)和公式(2)。

图3(b)是在头皮表面上建立的二维比例坐标系(简称为CPC坐标系)的示意图。其中，该二维比例坐标系为头皮表面上的任意点p向CPC空间提供一对一的映射。基于颅骨标记(Nz，Iz，AL，AR和Cz)和CPC坐标系(与尺度和形状成比例)定义的比例关系所确立的对象间对应关系，可以在个体层面的头皮表面上建立合理的解剖学对应关系。为了从单个视角可视化整个头皮表面，如图3(c)所示，在标准半球上建立了一个特殊的CPC空间。然后，使用现有的Hammer-Aitoff投影对带有CPC坐标系的半球进行平面化，生成在扁平椭圆上呈现的CPC坐标系的地图，申请人将其称为BNU地图(Beijing Normal University Map)(参见图3(d))，其实质为标准CPC坐标系的二维投影图像。在BNU地图的基础上，任何与头皮表面相关的脑功能数据都可以在该地图中呈现，从而能够在不同项目、人群、实验室甚至不同成像模式之间实现有效的比较。

2)建立用于连接颅骨位置与大脑位置的经颅脑映射模型(Transcranial brainmapping，简称为TBM)

一旦在个体层面的头皮表面建立了CPC空间，使用成熟的气球膨胀模型(Okamoto&Dan，2005)，可以在个体空间(例如，个体3D MRI图像)中确定对应于给定的头皮表面任意点p的底层皮层位置c。在所有皮层位置被空间归一化为标准脑空间(即MNI空间)之后，聚合所有(p，c)对，可以生成确定性的个体经颅脑映射模型。然后，通过集成所有的个体经颅脑映射模型来生成群体层面的概率经颅脑映射模型：

P(c|p) (3)

p(p_e,p_l)∈CPC，c(x,y,z)∈C，C是标准脑空间的子集，并包含经颅脑映射技术相关的所有皮层位置。当从给定坐标的头皮表面上任意点p(p_e，p_l)开始刺激或记录时，概率经颅脑映射模型给出每个靶向皮层位置c(x，y，z)的概率。

图4给出了单个CPC坐标所对应的概率经颅脑映射模型。当给出一对CPC坐标，例如P＝(0.4，0.6)，在每个个体层面的头皮表面上可以确定对应的点B，其所对应的皮层投射点C的位置使用成熟的气球膨胀模型在个体层面上的磁共振图像空间中识别出来。在图4中，下半部分的表格显示了任意点P＝(0.4，0.6)对应在概率经颅映射模型的不同位置的相应概率。

从颅骨空间一点到标号空间一个标号的映射，可以看做是一个两步映射。首先，从颅骨空间S映射到脑空间B，再从脑空间B映射到标号空间L。由于两步映射都是概率映射，因此，这个过程也可以看做一个两步随机过程。由于脑空间和脑分区标号的对应关系是由人脑自身的结构规律决定的，我们假定对于脑坐标空间中任意一点，其对应到各个脑区标号的概率是确定的，而与之前颅骨坐标空间到脑坐标空间的对应路径无关。因此，这个两步随机过程存在马尔可夫(Markov)性质。

3)使用马尔可夫链中的两步随机过程构建经颅脑图谱(Transcranial brainatlas)

本领域技术人员知晓，脑图谱是在概率框架中构建的。例如，传统的脑图谱(例如MNI图谱)描述的基本关系是条件概率：

P(l|b) (4)

其中，b(x，y，z)∈B，B是标准脑空间的子集，并且包含图谱的脑模板中的所有可能的脑组织点；l∈L，L包含所有可能的图谱标号，每个图谱标号都表示脑图谱中的特定脑区。对于给定的一对l和b，P(l|b)表示人脑中的位置b处出现图谱标号l可能性。

相应地，经颅脑图谱描述的基本关系也是条件概率：

P(l|p) (5)

其中p(p_e，p_l)∈CPC和l∈L。

在本发明的一个实施例中，可以通过马尔可夫链中的两步随机过程构建经颅脑图谱。具体地说，第一步：作为输入的给定点p(p_e，p_l)通过概率经颅映射P(c|p)被映射到标准脑空间中的皮层位置c(x，y，z)。第二步：特定的c(x，y，z)被映射到标号空间L中的标号l。研究者或者医生可以单独地通过皮层位置c(x，y，z)来预测图谱标号l，而无需考虑头皮表面上的位置p(p_e，p_l)，这里提供了公式(6)：

P(l|c,p)＝P(l|c) (6)

因此，两步随机过程使用了本领域技术人员常用的马尔科夫链。其中，假设皮层(域c)上的点是脑中(域b)的点的子集，P(l|p)可以使用公式(3)中的P(c|p)和公式(4)中的P(l|b)来计算。具体而言，如果p和c离散化，那么Chapman-Kolmogorov方程式表示：

其中，p_i是CPC空间中的离散位置(即任意点，下同)，i＝1，2，...，Np；c_j是C在标准脑空间中的离散位置，j＝1，2，...，Nc；l_k是特定脑图谱的图谱标号，k＝1，2，...，N₁。

需要说明的是，上述步骤所构建的经颅脑图谱是一种概率图谱，即当经颅设备(包括但不限于经颅脑成像装置如fMRI、fNIRS或经颅脑治疗装置如rTMS)的探头在给定坐标的头皮表面的任意位置p刺激或记录时，通过经颅脑图谱可以给出由l标记的每个靶向脑区域被探及的概率。它将不可见的脑内的图谱标号信息投射到可见的头皮上来，使得研究者或者医生可以直接使用这些脑结构和功能图谱信息，大大提升了脑图谱在经颅脑映射技术中的作用。

以上从整体上介绍了本发明所构建的经颅脑图谱生成方法。在实际应用中，每个人的颅骨形状、大脑结构存在明显的个体差异。上述采集个体特征构建的经颅脑图谱如何扩展到群体层面，从而囊括尽可能多的大脑内部结构共性特征，屏蔽尽可能少的大脑内部结构个体差异，就成为本领域研究者或者医生亟待解决的另一个技术课题。

为此，在本发明的一个实施例中，使用114个参与者的MRI数据集，提供了基于传统脑图谱的经颅脑图谱构建示例，并对构建的经颅脑图谱进行了验证，确认了本发明所提供的经颅脑图谱的可重复性和可预测性。

在上述两步随机过程的第二步中，本发明实施例考虑使用现有脑成像技术中三种通用的脑图谱实现。第一种是LPBA图谱(LONIProbabilistic Brain Atlas)，P_LPBA(l|b)，其中大脑皮质结构在标准脑空间基于宏观解剖分解被进行了标记(Shattuck，2008)。第二种是AAL标号图谱(Automated Anatomical Labeling Atlas)，P_AAL(l|b)，其中每个半球中的感兴趣区的宏观解剖结构在标准脑空间中的单个受试者的脑中被标记(Tzourio-Mazoyer等，2002)。最后一种是Talairach图谱，P_BA(l|b)，其中包括47个Brodmann区域(BA)标号，通过使用Brodmann方案的细胞类型在单个受试者的大脑中被定义(Talairach和Tournoux，1988；Lancaster等，2000)。

下面，首先对基于LPBA图谱的群体经颅脑图谱的生成过程及基于LPBA图谱的经颅脑图谱应用结果进行介绍，并对经颅脑图谱的可重复性和可预测性进行说明。

我们使用SLIM数据库(Southwest University Longitudinal ImagingMultimodal)(http://fcon_1000.projects.nitrc.org/indi/retro/southwestuni_qiu_index.html)中的114名青年参与者(年龄：18-24岁，63女性和51男性)的结构MRI(sMRI)数据集。使用磁化设备的快速采集梯度回波序列(MPRAGE)获取高分辨率3DT1加权结构图像(TR/TE＝1900ms/2.52ms，FA＝9°，FOV＝256×256mm²；切片＝176；厚度＝1.0mm；体素大小＝1×1×1mm³)。

对于每个参与者，使用MRIcron软件在个体层面上的3DMRI图像中利用视觉识别出Nz、AL、AR和Iz四个颅骨标记(识别颅骨标记的示意参见图2所示)，然后提取头皮表面和皮质表面。具体地说，使用SPM12(Wellcome Trust Centre for Neuroimaging,London,UK.http://www.fil.ion.ucl.ac.uk/spm)中的统一分割算法，将每个个体的3D MRI图像分割成六个组织图像：灰质、白质、脑脊液(CSF)、骨、软组织和空气/背景。生成脑图像(灰质+白质)和头图像(灰质+白质+CSF+骨+软组织)。在使用Gaussian kernel和FWHM[3,3,3]对脑图像和头图像平滑处理，以及使用阈值0.5进行二值化处理后，将SPM12中的表面提取算法应用于二值化图像，以提取参与者的头皮表面点云(粉红色)和皮层表面点云(灰色)(如图1b所示)。

在确定四个颅骨标记Nz、AL、AR和Iz的基础上，使用Jurcak等人提供的迭代算法(参见Jurcak,V.,D.Tsuzuki,and I.Dan,10/20,10/10,and 10/5systems revisited:their validity as relative head-surface-based positioning systems.Neuroimage,2007.34(4):p.1600-11.)将标记Cz标记在个体头皮表面点云中。然后，群体层面的经颅脑映射模型由公式(8)生成：

其中，p_i是CPC中的网格位置，i＝1，2，...，Np；c_j是通过经颅脑映射技术获得的标准脑空间中的皮质体素，j＝1，2，...，Nc。#是个体总数。例如，在图4中给出了头皮表面任意点p(0.4，0.6)的经颅脑映射模型114列表，并且，图4中灰色脑模型上的色块表示了给定颅面坐标点p＝(0.4，0.6)对应皮层投影点的空间频率分布P(c|p)。

考虑到公式(8)所需的离散形式，本发明实施例通过将p_e和p_l的整个范围分别均匀分割为100个分段将连续CPC空间离散化，生成一个均匀的网格，命名为CPC₁₀₀(如图3(b)所示)。CPC₁₀₀中两个相邻点之间的距离约小于3.5毫米，空间分辨率与大多数经颅脑映射技术兼容。

对于CPC₁₀₀中的每个点p，例如图4(c)中的(0.4，0.6)，相应的头皮位置s(x，y，z)(图4中的粉红点)由公式(1)和(2)从个体头皮点云中确定。然后，将成熟的气球膨胀模型应用于头皮位置s(x，y，z)以确定一个对应的皮层位置(图4中的黄点)。接着，将个体空间中的皮层位置(图4中的黄点)空间标准化到标准脑空间中以获得c。聚合所有(s，c)对生成个体层面上的经颅脑映射模型，即获得从CPC₁₀₀网格到标准脑空间的映射。

在实际计算中，我们对皮层投影点的MNI坐标进行了调整。第一，根据图谱图像的分辨率对所有投影点的MNI坐标进行了空间重采样。第二，对少数皮层投影点的MNI坐标进行了修正。由于配准误差，少数的皮层投影点标准化到标准脑空间后并不在脑图谱的范围，实际计算中，我们对这些偏离的点进行了修正，用最近邻搜索(Nearest NeighborSearching)将偏离的投影点限定到脑图谱的空间范围内。

作为结果，我们可以得到CPC坐标点到图谱体素的映射矩阵。例如，对LPBA图谱，共有136020个灰质体素。经过颅-脑投影和空间标准化，每个被试者的颅脑映射可以表示为一个9801×136020的二值矩阵。在群体层面上对该二值矩阵中的点求频率，得到矩阵Matrix_PB，Matrix_PB中的行表示估计得到的条件概率的分布律。对于AAL图谱，Matrix_PB的大小是9801×185355。对于BA图谱，Matrix_PB的大小是9801×403482。

根据脑图谱图像提供的信息，我们将图谱信息转化为体素到分区标号的映射矩阵。其中，对于AAL图谱和BA图谱，脑图谱图像的文件中为显式给出概率分布，因此，我们定义：

以此将分区编号转换为概率形式。

对于LPBA图谱，脑图谱图像提供的信息中直接给出了各个标号在各个空间位置出现的概率。因此，我们直接逐个索引脑图谱图像中的脑区域体素，并记录各个标号在该体素的概率值。

对整个标准脑空间中的位置，存在该位置不属于任何标号区域的情况，这个事件的对立事件是该位置至少被标号为某一脑区。我们显式地将“非脑”作为一种特殊的脑区标号，并定义其概率为：

使条件概率P(L|B)满足概率性质：

我们将上述得到条件概率分布记为矩阵形式Matrix_BL。对于AAL图谱，Matrix_BL大小为136020×121。对于BA图谱，Matrix_BL大小为185355×47，对于LPBA图谱，Matrix_BL大小为1336020×57。

根据上述步骤得到的Matrix_PB和Matrix_BL求解颅骨坐标到脑区标号的概率映射矩阵Matrix_PL，其中

Matrix_PL＝Matrix_PB·Matrix_BL (11)

对于LPBA图谱，Matrix_PL大小为9801×57；对于AAL图谱，Matrix_PL大小为9801×121；对于BA图谱，Matrix_PL大小为9801×47。Matrix_PL中的行，表示给定一个CPC颅骨坐标位置映射各个分区标号的映射概率；Matrix_PL中的列，表示CPC坐标映射到给定脑区标号l_k的条件概率P(L＝l_k|P)，这个概率表示了各个CPC坐标点对应到l_k的可能性。

确定两个半球以及四个主叶在CPC空间的定位是展示经颅脑图谱的第一步。构建主叶级经颅脑图谱并将其呈现在图5(a)和5b中。从图5(b)中的最大似然标号图(MLM)可以看出，在CPC空间中保留了脑空间中的四叶之间的原始空间拓扑结构，如双边对称性。图5(a)中的最大概率图(MPM)表明经颅脑图谱标号在人群上较高的一致性(中值(Median)＝0.946)。此外，图5(a)中明显的狭窄边界粗略地对应于划分脑叶的脑沟结构。

在本发明实施例中，具有35个子区域的经颅脑图谱(经颅脑图谱_LPBA)在BNU地图(图5(d)所示)以及三个不同视图中的立体定位头皮表面(图5(e))上呈现，其具有从LPBA图谱继承的颜色编码方案。原始LPBA图谱中有56个子区域，但在经颅脑图谱中21个子区域不可见。不可见的子区域主要位于大脑内侧和腹侧部分，经颅脑映射技术不可及。此外，经颅脑图谱中可见但小的子区域可能对应于原始脑图谱中的大结构。例如，楔前叶就是这样的子区域，其绝大部分位于中央纵裂内，因此经颅脑图谱和相应脑图谱中相同标号区域的尺寸不一定具有可比性。图5(c)中的经颅脑图谱_LPBA的MPM在每个子区域显示高的一致性(高达98％)，而低的一致性仅发生在边界附近。

综上所述，本发明首先提供了一个基于两步马尔可夫链模型的概率框架作为经颅脑图谱的理论基础。其中，第一步是从CPC空间中的头皮位置到MNI空间中的下划线皮层位置的颅骨-皮层映射(crainio-cortical mapping)。第二步是使用传统的大脑图谱构建经颅脑图谱，实质是一个从MNI空间的皮层位置到图谱标号空间的映射。需要说明的是，本发明利用概率框架提供了一个可扩展的经颅脑图谱模型，上述第二步所使用的大脑图谱可以由任何其他脑图谱替代。虽然，在本发明中仅提供了三种基于宏观解剖学的脑图谱(BA图谱，AAL图谱和LPBA图谱)用于构建经颅脑图谱，但类似的功能图谱、连接图谱以及其他图谱都可以被用来为特定应用提供功能性的经颅脑图谱。

经颅脑图谱的有效性主要体现在两个主要方面。首先，经颅脑图谱的构建是通过对人群部分个体的采样实现对人群水平解剖信息的估计，因此，经颅脑图谱在人群不同采样上构建的结果应该是一致的。其次，经颅脑图谱最终需要利用群体水平的解剖知识，实现对个体经颅数据的定位和导航，因此，基于经颅脑图谱实现的群体-个体可预测性也是评价经颅脑图谱有效性的另一个重要指标。

为了验证我们构建的经颅脑图谱，本发明实施例中以两种方式量化了它的性能。首先，我们测量了经颅脑图谱的可重复性。对于经颅脑图谱，高可重复性意味着从相同群体的不同样本中构建的群体经颅脑图谱是相似的。第二个有效性指标是可预测性。高可预测性意味着我们可以在较高的置信水平下以一个群体经颅脑图谱替代一个个体经颅脑图谱。在缺少个体结构像的情况下，这意味着实验者可以利用群体经颅脑图谱预测个体受试者或者患者头上的探头放置。

为了评估群体经颅脑图谱的可重复性，本发明实施例中将114名参与者的结构性sMRI数据集随机分为两组(GA和GB)，每组有57名参与者。GA和GB分别用于构建经颅脑图谱57A和经颅脑图谱57B。通过估计经颅脑图谱57A和经颅脑图谱57B之间的一致性来评估群体经颅脑图谱的可重复性。我们使用DICE指数来评价每个标号区域l的标号一致性。我们使用经颅脑图谱57A和经颅脑图谱57B的最大似然标记(l*_57A(p)和l*_57B(p))来计算DICE指数：

其中，||给出某一标号区域内的面积，X_lA和X_lB分别是l*_57A(p)和l*_57B(p))中标号l的区域。DICE指数的范围为0到1，较大的值表示较高的一致性。

计算得到的DICE值(Median＝0.95)请参见表1。

表1经颅脑图谱_LPBA的分区域DICE值

为了定量研究预测准确性，将114名参与者的sMRI数据随机分为构建组(GC，92名参与者)和测试组(GT，22名参与者)。在构建组上构建的经颅脑图谱用于预测测试组中每个参与者的个体经颅脑图谱。也就是说，将来自群体经颅脑图谱的最大似然标号与来自个体经颅脑图谱的相应标号进行比较，对每个CPC点p计算出预测精度的正确率，得到预测精度图。

图6(a)显示的任意一个测试组的参加者的个体经颅脑图谱和图6(b)所示的预测得到的经颅脑图谱(即经颅脑图谱92)非常相似。汇总所有22个单独预测性能的群体层面预测精度图(图6(c))显示出总体高预测精度(中值＝0.96)。在图6(d)中，用黄色编码的区域是精确度高于90％的区域，而用红色编码的区域预测精度低于90％，从中可以看出，精确度偏低的区域位于边界附近。为了评估构建的群体经颅脑图谱对于其他种族参与者的可预测性，一个由24名白人参与者(年龄：23.43±4.6；17名男性和7名女性)组成的群体被用作另一个测试组，得到了一个略低的预测精度(中值＝0.92)。

下面继续对经颅脑图谱_AAL和经颅脑图谱_BA的预测结果进行简单介绍。

经颅脑图谱_AAL中有76个子区域，其中原始AAL图谱中的120个子区域(每个半球为60个)中有44个子区域是不可见的。MPM显示出标号的整体高度一致性(中值＝0.86)，且仅在较小的子区域和接近边界处发生低度一致性，特别是在枕叶和靠近侧面沟槽处。区域性DICE值相当高(中值＝0.86)。群体层面预测精度的中值高达0.91，表明总体高预测精度。为了评估构建的经颅脑图谱对于其他参与者的可预测性，一组24名白人参与者(年龄：23.43±4.6；17名男性和7名女性)被用作另一个测试组，结果表明预测精度略低(中值＝0.88)。

经颅脑图谱_BA中有27个子区域，其中原始Talairach图谱中的50个子区域有23个不可见。MPM同样表现出总体高的标记一致性(中值＝0.87)，但在中央沟和侧龈附近的长而窄的子区域中表现为低度一致性。区域性DICE值相当高(中值＝0.92)。群体层面预测精度的中值高达0.91，表明总体高预测精度。为了评估构建的经颅脑图谱114对于其他参与者的可预测性，一组24名白人参与者(年龄：23.43±4.6；17名男性和7名女性)被用作另一个测试组，结果表明预测精度略低(中值＝0.88)。

上述研究可以得到一个重要结论：尽管存在个体差异，不同被试对象的同一颅骨地标点对应的皮层投影点，具有脑沟回或脑图谱分区级别的对应性。这种对应性一方面保证了基于颅骨位置的经颅成像装置放置，能够确保不同人之间沟回级别的一致性；另一方面，经过验证的10-20地标点的皮层对应位置，可以用于在不依赖个体解剖或脑成像信息的条件下，预测经颅成像装置可探及的皮层位置。

需要说明的是，对于一些经颅映射技术如fNIRS，很少进行额外的sMRI扫描。幸运的是，群体经颅脑图谱可以在这里用作现有知识，以显示特定目标的最可能的放置位置，以及每个头皮位置最可能的解剖标号。前已述及，我们验证了群体经颅脑图谱的可预测性。实验证明在我们的数据中，中值预测精度高于0.9，错误主要发生在标号区域的边界附近。对于其他更重视准确性的经颅映射技术(如TMS)，sMRI扫描更为常见。当sMRI数据可用时，可以基于个体经颅模型进一步构建个体经颅脑图谱，从而提供更有保证的标记精度。

在上述面向群体应用的经颅脑图谱生成方法和预测方法的基础上，本发明进一步提供一种可以在临床治疗中应用的经颅脑图谱预测装置。在本发明的一个实施例中，该经颅脑图谱预测装置可以制作成头盔或者头套形状，以覆盖较为完整的头皮表面。通过本发明获得的群体经颅脑图谱图像可以采用印刷或者喷绘方式覆盖在经颅脑图谱预测装置的上表面，该经颅脑图谱预测装置的下表面在实际使用时紧贴用户的头皮表面。

参见图7所示，作为上述实施例的一个变型例，该经颅脑图谱预测装置可以简化为类似泳帽的印有群体经颅脑图谱图像的头套。该头套由棉布或者化纤材料制成，成本低廉，便于大规模推广使用。在实践中，当利用上述经颅脑图谱生成方法获得某个用户的个体经颅脑图谱图像时，也可以将该图像直接喷绘或者打印在空白的头套上，从而在头套上呈现出反映用户个体特征的经颅脑图谱图像。用户在进行相应的脑科临床治疗时，可以随身携带该头套帮助医生进行准确定位。

在经颅脑图谱预测装置的表面上分布有若干小孔，以便让经颅成像装置或者经颅治疗装置的探头通过上述小孔接触用户的头皮表面，从而实现相应的颅内大脑观测或者治疗操作。上述供探头穿过的小孔优选根据群体经颅脑图谱图像中的不同子区域排列成点阵形状，其中对于面积较小的子区域，小孔的排列密度可以大一些(即排列更加密集)；对于面积较大的子区域，小孔的排列密度可以小一些(即排列相对疏松)。这样可以保证探头在不同子区域操作时，无论区域面积大小都能找到足够的操作位置。

相比传统脑图谱，本发明所提供的经颅脑图谱能够有效解决经颅脑成像研究中的定位问题。第一，颅骨坐标系统的建立，能够对实现对整个颅骨表面空间，也就是经颅成像装置的放置空间进行精确描述，保证了放置位置在不同个体之间的可重复性。第二，标准颅骨空间到标准脑空间的对应关系，使得研究者或者医生可以直接从经颅成像装置在颅骨的放置位置获得对应的MNI空间坐标，因此解决了在无MRI结构像情况下经颅数据的定位问题。第三，反向显示在被试颅骨表面的脑标号信息，便于以直观的方式指导经颅成像装置的准确放置。

本发明所提供的经颅脑图谱将不可见的脑内图谱标号信息投射到可见的头皮上来，使得研究者或者医生可以“直接”使用这些脑结构和功能图谱信息，大大提升了脑图谱在经颅脑映射技术使用中的作用。

以上对本发明所提供的面向群体应用的经颅脑图谱生成方法、预测方法及其装置进行了详细的说明。对本领域的一般技术人员而言，在不背离本发明实质精神的前提下对它所做的任何显而易见的改动，都将构成对本发明专利权的侵犯，将承担相应的法律责任。

Claims (15)

1.一种经颅脑图谱生成方法，其特征在于包括如下步骤：

(1)在个体层面创建颅面坐标系；

(2)建立用于连接颅骨位置与大脑位置的经颅映射系统；

(3)使用马尔可夫链中的两步随机过程构建经颅脑图谱。

2.如权利要求1所述的经颅脑图谱生成方法，其特征在于所述步骤(1)包括如下子步骤：

(11)在头皮表面识别五个颅骨标记Nz、Iz、AL、AR和Cz；

(12)将头皮表面和通过Nz、Cz和Iz的平面之间的相交曲线定义为颅骨赤道；

(13)给出头皮表面上的点p，经度曲线可以唯一地确定为头皮表面和通过AL、AR和p的平面之间的相交曲线，p'是颅骨赤道与经度曲线之间的交叉点；

(14)上头皮的任意点p由一对非负实数(p_e,p_l)唯一确定：

p_e＝L_NZ-p’/L_e,p_e∈[01]

p_l＝L_AL-p/L_AL-p-AR,p_l∈[01]

其中，L_Nz-p'是沿着颅骨赤道的Nz到p'的曲线长度，L_e是颅骨赤道的全长；L_AL-p是沿着全长为L_AL-p-AR的经度曲线从AL到p的曲线长度。

3.如权利要求2所述的经颅脑图谱生成方法，其特征在于所述步骤(1)还包括步骤(15)：在标准半球上建立CPC空间；使用Hammer-Aitoff投影对带有CPC坐标的半球进行平面化，生成在扁平椭圆上呈现的CPC坐标系的地图。

4.如权利要求1所述的经颅脑图谱生成方法，其特征在于所述步骤(2)包括如下子步骤：使用气球膨胀模型在个体空间中确定对应于给定的头皮表面任意点p的底层皮层位置c；在所有皮层位置被空间标准化为MNI空间之后，聚合所有(p，c)对，生成确定性的个体经颅脑映射模型。

5.如权利要求4所述的经颅脑图谱生成方法，其特征在于所述步骤(2)还包括如下步骤：

集成所有个体模型来生成群体层面的概率经颅脑映射模型：P(c|p)；

其中，p(p_e,p_l)∈CPC，c(x,y,z)∈C，C是MNI空间的子集。

6.如权利要求5所述的经颅脑图谱生成方法，其特征在于所述步骤(2)中，所述概率经颅脑映射模型由下列公式生成：

其中，p_i是CPC₁₀₀网格位置，i＝1，2，...，Np；c_j是通过经颅脑映射技术获得的MNI空间中的皮质体素，j＝1，2，...，Nc；#是个体总数。

7.如权利要求5所述的经颅脑图谱生成方法，其特征在于所述步骤(3)还包括如下步骤：

(31)给定点p(p_e，p_l)通过概率经颅映射P(c|p)被映射到MNI空间中的皮层位置c(x，y，z)；

(32)皮层位置c(x，y，z)被映射到标号空间L中的标号l。

8.如权利要求7所述的经颅脑图谱生成方法，其特征在于在步骤(3)中，假设皮层域c上的点是脑中域b的点的子集，如果p和c离散化，P(l|p)由Chapman-Kolmogorov方程式表示：

其中，p_i是CPC空间中的离散位置，i＝1，2，...，Np；c_j是C在MNI空间中的离散位置，j＝1，2，...，Nc；l_k是特定脑图谱的标号，k＝1，2，...，N₁。

9.如权利要求7所述的经颅脑图谱生成方法，其特征在于所述步骤(32)中使用LPBA40脑图谱、AAL标号图谱或Talairach图谱中的任意一种。

10.如权利要求1所述的经颅脑图谱生成方法，其特征在于还包括步骤(4)：生成最大似然标号图和/或最大概率图。

11.一种面向群体应用的经颅脑图谱预测方法，其特征在于包括如下步骤：

针对群体中的单独个体，当经颅设备的探头在给定坐标的头皮表面的任意位置进行刺激或记录时，通过权利要求1～10中任意一项所述的经颅脑图谱给出每个靶向脑区域被探及的概率。

12.如权利要求1所述的经颅脑图谱预测方法，其特征在于：

所述经颅设备为经颅脑成像装置或经颅脑治疗装置中的任意一种。

13.一种面向群体应用的经颅脑图谱预测装置，用于实施权利要求11或12所述的经颅脑图谱预测方法，其特征在于：

所述经颅脑图谱预测装置为头盔或者头套形状，其上表面覆盖有权利要求1～10中任意一项所述的经颅脑图谱，下表面在使用时紧贴用户的头皮表面。

14.如权利要求13所述的经颅脑图谱预测装置，其特征在于：

在所述经颅脑图谱预测装置的表面分布有若干小孔，用于供经颅设备的探头通过所述小孔接触用户的头皮表面。

15.如权利要求13所述的经颅脑图谱预测装置，其特征在于：

所述小孔根据经颅脑图谱图像中的不同子区域排列成点阵形状，其排列密度根据所述子区域的面积呈反向变化。