经颅磁刺激导航个性化模型构建方法

张 芸 丁 辉 王广志

清华大学生物医学工程系，北京 100084

引言

经颅磁刺激(transcranial magnetic stimulation，TMS)是一种无创安全、操作简单的神经刺激技术[1]。作为神经调控的重要工具，TMS 广泛应用于神经精神类疾病的诊断与治疗。其中，TMS 线圈不准确、不稳定和难以重复定位，会导致神经系统疾病诊断结果存在变异性[2]，也会导致精神类疾病治疗存在异质性[3]。无框架立体定位系统与TMS 技术结合，可以很好地解决TMS 线圈定位的问题。这种导航TMS(navigated TMS，nTMS)的主要原理是基于个性化三维影像，制定个性化刺激计划，通过三维影像与物理空间注册，实现刺激计划到物理空间的映射，最终实现影像引导下刺激线圈的准确定位定向[4]。

影像引导是nTMS 的重要组成部分，它的前提是拥有患者的磁共振影像(magnetic resonance imaging，MRI)。在临床上，一些禁忌症，如体内存在电磁性、铁磁性或机械植入设备(包括心脏复律除颤器、心脏起搏器、金属碎片、血管夹、人工耳蜗或眼内植入物)，以及患有幽闭恐惧症或癫痫病等，造成部分患者无法获得MRI。同时，MRI 扫描费用会额外增加患者的经济负担，造成患者对nTMS 的接受度降低。这些情况都会导致nTMS 无法使用。

个性化近似影像已经被广泛应用在脑电溯源的研究中[5]。由于其中的近似影像构建需要佩戴脑电帽，且脑电溯源的关注精度与nTMS 不同，所以这些近似影像构建方法及精度结果在nTMS 中不完全适用，而针对nTMS 的近似影像研究较少。2012年，Carducci 等报道了基于近似影像nTMS 的定位可靠性[6]，但是没有公开近似影像的构建方法。2020年，Fleischmann 等对比发现，基于MRI 和基于近似影像引导的nTMS 评估运动皮层兴奋性的结果具有一致性，进一步说明基于近似影像nTMS 的可靠性，然而其中近似影像的构建需要借助外部工具箱[7]来实现，过程复杂且依赖于技术人员的经验，同时也没有针对近似影像的结构精度进行分析。综上分析，专门针对nTMS 的近似影像便捷构建方法、精度研究以及精度是否满足一些nTMS 定位需求都是未知的。

为了解决上述问题，本研究提出一种基于标准模板变形的方法，实现个性化模型的构建，以替代MRI 用于nTMS。首先获取患者头皮轮廓信息，其次重建成头皮表面，并在重建头皮和标准模板头皮上定位变形控制点，最后结合控制点配对信息，实现标准模板到患者头皮的变形，建立个性化模型。通过实验分析个性化模型的头皮精度、脑表面精度以及脑靶点精度，探讨模型在临床上的适用性。

1 材料和方法

图1为个性化模型构建流程。通过标准模板空间到患者头部空间的变形变换，实现标准模板头皮到患者头皮的匹配，从而建立个性化模型。首先，利用定位装置，获取患者头皮的轮廓点信息；其次，借助泊松表面重建方法，实现轮廓点到头皮表面的重建；接着，通过定位多条头皮曲线，实现重建头皮与标准模板头皮上控制点的定位；最后，基于控制点对信息结合薄板样条插值方法，实现标准模板到患者头皮的变形。

图1 个性化模型构建流程Fig.1 Individualized model construction flow chart

1.1 头皮轮廓点获取

TMS 过程往往不要求患者剃发。所以，大面积头皮轮廓点的获取是个难题。为解决这个问题，并结合后续轮廓点的处理需求，利用nTMS 光学定位装置及探针来拾取头皮轮廓点。如图2所示，患者佩戴跟踪头架，探针尖端轻触患者头皮，光学定位装置记录多个跟踪工具的空间位姿信息，包括定位装置坐标系下探针坐标系的实时位姿、 跟踪头架坐标系的实时位姿。由于探针尖端在探针坐标系下的位姿已知，头皮轮廓点在跟踪头架坐标系下的位姿(即探针尖端位姿)，可以计算如下:

图2 头皮轮廓点获取Fig.2 Diagram of scalp contour points obtaining

为了提高拾取效率，在nTMS 交互界面提示下完成整个头部轮廓点的拾取。为了保证完整性，头皮轮廓点需要覆盖到鼻子下方。同时，为了后续控制点定位，需要拾取3 个面部解剖点，包括左耳前端、右耳前端和鼻根。

1.2 头皮泊松重建

人工获取的头皮轮廓点分布不均匀，作为控制点来约束标准模板变形，存在标准模板头皮上对应控制点难定位的问题，容易造成变形过程中解剖结构的扭曲错位。为了解决这个问题，需要在患者头皮和标准模板头皮上定位统一的控制点，而控制点定位的前提是拥有平滑的头皮表面。本研究选用了融合全局和局部优点的泊松重建算法[8]，实现头皮轮廓点到头皮表面的重建。泊松重建的基本思想是建立表面采样点与表面指示函数的积分关系，结合划分块的方法，获得采样点的向量场，计算指示函数梯度场的逼近，从而建立泊松方程。通过拉普拉斯矩阵迭代实现泊松方程求解，并结合移动立方体法提取等值面，得到被测物体的表面。由于泊松重建算法较为成熟[8]，这里就不再赘述。

1.3 控制点定位

为了实现重建头皮和标准模板头皮上控制点的统一定位，笔者选用一种类似脑电电极定位的方法来确定控制点。脑电电极[9]均匀分布于整个头部，且电极位置具有一定的解剖学意义。尽管本研究的控制点定义方法类似国际脑电电极系统，但是分布更加稠密。控制点位置由头皮解剖点定义的多条曲线生成，解剖点包括左耳前端、右耳前端和鼻根以及由它们定义的枕骨隆凸和顶点。首先，如图3(a)所示，定义解剖点、矢状线和横断线。枕骨隆凸的定义是鼻根与左右耳前端中心连线与头皮的交点。顶点的定义是左右耳前端中心，沿鼻根—左耳前端—右耳前端平面法向量的射线与头皮的交点。矢状线的定义是头皮与鼻根—顶点—枕骨隆凸平面相交形成的曲线。定义距离鼻根和枕骨隆凸10%矢状线长度的点为鼻根辅助点和枕骨隆凸辅助点，基于这两个点定义平行于鼻根—左耳前端—枕骨隆凸平面的横断平面，该平面与头皮相交形成横断线。其次，如图3(b)所示，定义冠状线。在横断线上，以距离鼻根辅助点恒定长度放置控制点；在矢状线上，以距离鼻根辅助点恒定长度放置控制点。基于横断线上左右对称点和矢状线上对应点确定冠状平面，与头皮相交，形成冠状线。最后，如图3(c)所示，定义所有控制点。缩短矢状线、横断线和冠状线上控制点间隔，同时增加多条辅助曲线并忽略重叠点，将控制点数量增加至几百个。按照上述控制点定位方式，实现了重建头皮和标准模板头皮的控制点定位。

图3 控制点定义示意。(a)定义解剖点、矢状线、横断线；(b)定义冠状线；(c)定义所有控制点Fig.3 Diagram of control points definition. (a)Anatomical landmark points，sagittal lines and transverse lines definition； (b)Coronal lines definition； (c)All control points definition

1.4 标准模板变形

为了基于控制点实现标准模板到重建头皮的变形，选用了一种基于特征点匹配的非刚性图像配准算法——薄板样条插值算法(thin plate spline，TPS)[10]。TPS 方法是一种散乱数据插值变形方法，它可在标准模板x和重建头皮y之间建立一个光滑的三维位移插值函数f，从而实现三维网格的非刚性变形，且保持控制点对的严格重合。该三维插值函数由仿射变换部分和非线性变换部分构成，有

式中:Mx + t为全局仿射变换多项式，φ(x -xj)为局部非线性变换多项式；λj为控制点xj对应的权重系数；φ(x - xj)为TPS 基函数，定义

将m个控制点xj代入式(2)构成方程组，通过约束条件计算得到λj、M和t的值，最终得到标准模板到重建头皮的变形函数f(x)。

1.5 验证

在nTMS 中，个性化影像数据用于确定刺激计划。8 字形线圈(70 mm 直径Alpha 线圈)作为经颅磁刺激仪(Rapid2，Magstim，UK)的常用刺激工具，其有效刺激深度可到达头皮下方10 ～20 mm[11]，与头皮到脑表面距离接近[12]，主要实现大脑皮层的刺激。所以所关注的刺激计划包括基于功能信息确定脑靶点、基于脑靶点确定脑表面靶点、基于脑表面靶点确定头皮靶点。由本方法构建的个性化模型在用于nTMS 时，其头皮、脑表面以及脑靶点的精度决定了其是否能够适用于nTMS。通过实验对个性化模型的头皮、脑表面以及脑靶点精度进行验证，nTMS 定位装置选用光学追踪设备(Polaris Vicra，Northern Digital Inc，Canada)进行头皮轮廓点位置测量，其测量精度为0.25 mm。

1.5.1 头皮精度验证

个性化模型的构建过程涉及多个可控制参数，如头皮重建中轮廓点数量、标准模板变形中控制点数量等。为了量化这些参数对个性化模型精度的影响，笔者以头型模型为研究对象设计实验，探讨轮廓点数量对重建重复性精度和重建精度的影响，以及控制点数量对变形精度的影响。

为了研究轮廓点数量的影响，利用探针在头型模型表面拾取3 000 个点Pi(i=1，2，…，5)，重复该过程5 次。从Pi中随机采样n个点Pi，n(n=500，1 000，…，3 000)，由Pi，n重建得头皮Si，n，比较Si，n与Si+1，n的距离，计算距离平均值Di，n和标准差Stdi，n。由于重建头皮覆盖范围存在不完全重叠，所以只将距离小于2 mm 的点纳入计算范围。计算作为重建重复性精度和标准差。以头型模型的计算机断层影像(computed tomography，CT)头皮为标准，比较Si，n与CT 头皮的距离，计算距离平均值，作为重建精度，对5 次结果进行统计。同样，考虑二者不完全重叠，所以只将距离小于3 mm 的点纳入计算范围。其中，CT 头皮由CT 数据通过3D Slicer 软件[13]，手动提取。

为了研究控制点数量的影响，在头型模型CT头皮上定位m个控制点Pm(m=16，32，64，128，192，256)，基于Pm变形得头皮Sj，m(j=1，2，…，10)，重复该过程10 次。比较Sj，m与CT 头皮的距离，计算距离平均值Dj，m和标准差Stdj，m，由于变形头皮与CT头皮存在不完全重叠，所以只将控制点覆盖区域的点纳入计算范围。计算，作为变形精度和标准差。

为了评估基于合适的轮廓点数量和控制点数量构建的个性化模型精度，以头型模型和真实人体数据为研究对象，构建个性化模型；将重建头皮和个性化模型头皮分别与CT 头皮(MRI 头皮)对比，比较二者距离。其中，MRI 头皮由MRI 数据通过3D Slicer 软件[13]，手动提取。

1.5.2 脑表面精度验证

脑表面精度主要受标准模板结构特征影响。为了量化分析不同标准模板对脑表面精度的影响，以不同标准模板和真实人体数据为研究对象，构建个性化模型，并以人体MRI 脑表面为标准，评估个性化模型脑表面精度。为了深入探讨标准模板结构特征对脑表面精度的影响，设计实验对标准模板和人体MRI 结构特征进行量化分析。

为了评估个性化模型的脑表面精度，在人体MRI 头皮上采样轮廓点，基于标准模板构建个性化模型。个性化模型的脑表面为S，对应MRI 脑表面为SMRI，计算两者距离D和距离平均值。按照图3(b)中的冠状线对S进行三等份分块，得到近似额叶区、运动区和枕叶区。分别对每个区域的D统计距离平均值其中，标准模板包括了头皮和脑表面数据，它们是由MRI 数据，通过3D Slicer 软件[13]，手动提取。

为了量化分析标准模板结果特征，对标准模板、人体MRI 以及个性化模型的尺寸和头皮到脑表面的距离进行了量化计算。根据TPS 算法原理，全局仿射变换为主，局部非线性变换为辅。所以，标准模板变形过程可以近似为仿射变换过程，而仿射变换的参数由标准模板头皮与重建头皮的尺寸比例决定。同时，当变形后的标准模板头皮与重建头皮重合时，脑表面精度表现为头皮到脑表面距离的差别。综上分析，脑表面精度主要受标准模板头颅尺寸和头皮到脑表面距离的影响。所以，本研究对标准模板、人体MRI 和个体化模型的头颅尺寸以及头皮到脑表面距离(scalp to cortex distance，SCD)进行计算，计算方法如图4所示。由图3(a)中的横断线，对头皮进行截取而获得横断面，对其中头皮边界进行椭圆拟合，该椭圆的长轴与短轴分别作为头颅的长与宽，该横断面到顶点的距离作为高，SCD是指顶点到横断面间头部的头皮到脑表面距离的平均值。

图4 头颅尺寸计算Fig.4 Head size calculation diagram

1.5.3 脑靶点精度验证

不同TMS 应用场景需要不同脑靶点。当nTMS用于诊断神经系统疾病时，常常刺激手部肌肉对应脑区M1-HAND 来观察运动响应。当nTMS 用于精神类疾病治疗时，背外侧前额叶(dorsolateral prefrontal cortex，DLPFC)是常用靶区。为了探讨基于个性化模型的nTMS 在上述应用场景中的适用性，以真实人体数据为研究对象，构建个性化模型，并以人体MRI 脑靶点为标准，评估个性化模型的脑靶点精度。

为了量化分析脑靶点精度，以左侧M1-HAND脑区中心点和左侧DLPFC 脑区中心点为研究对象(以下分别简称为M1-HAND 和DLPFC)。计算个性化模型的M1-HAND 和DLPFC 与MRI 对应点的距离DMRI和DDLPFC，分别作为M1-HAND 精度和DLPFC 精度。M1-HAND 是根据解剖学知识[14]手动选定，DLPFC 是Talairach 坐标(-45，45，35)[15]对应点的脑表面最近点。脑电电极点F3 是一种简便且流行的DLPFC 替代靶点，根据控制点定位结果，确定F3 的脑表面最近点，并计算个性化模型的F3与MRI 的DLPFC 的距离DF3，以作为F3 精度。

2 结果

2.1 头皮精度

从真实人体三维影像中提取头部表面，并利用3D 打印技术制作头型模型。对该头型模型进行CT扫描，CT 体素大小为0.5 mm×0.5 mm×0.8 mm，尺寸为512×512×238。

轮廓点数量对重建重复性精度的影响如图5(a)所示，其中每个点为5 次测量的平均结果。随着轮廓点数量的增加，重复性精度呈现下降趋势，且均小于0.5 mm，说明点云重复性精度很好；标准差也随着轮廓点数量的增加呈现下降趋势，说明重建头皮变得越来越光滑。同时，随着轮廓点数量的增加，这种下降趋势变小。轮廓点数量对重建精度的影响如图5(b)所示，随着轮廓点数量的增加，误差呈现下降趋势；当轮廓点数量增加到2 000 时，这种下降趋势变小。因为轮廓点拾取过程的交互复杂性和耗时，所以选择2 000作为所需头皮轮廓点的最低数量。控制点数量对变形精度的影响如图5(c)所示，其中每组数据为10 次测量的平均结果。随着控制点数量的增加，变形误差有减小的趋势，在256 个点对时最小。由于控制点定位过程的自动化实现且耗时较短，所以选择256 个控制点。在头型模型表面和1 例真人头皮表面分别拾取2 000 个点云，并定位256 个控制点，重复5次，得到重建头皮和个性化模型头皮，如图6所示。对5 次结果统计，头型模型的重建头皮精度为(0.9±0.1)mm，个性化模型的头皮精度为(1.5±0.3)mm。真人的重建头皮精度为(2.3±0.2)mm，个性化模型的头皮精度为(2.8±0.5)mm。

图5 头皮精度影响因素。(a)重建重复性精度随轮廓点数量的变化；(b)重建精度随轮廓点数量的变化；(c)变形精度随控制点数量的变化Fig.5 Factors that affects scalp accuracy. (a)The relationship between reconstruction repeatability accuracy and contour points' number； (b)The relationship between reconstruction accuracy and contour points' number； (c)The relationship between warping accuracy and control points' number

图6 重建头皮与个性化模型头皮。(a)头型模型的头皮轮廓点；(b)头型模型的重建头皮精度；(c)头型模型的个性化模型头皮精度；(d)真人的头皮轮廓点；(e)真人的重建头皮精度；(f)真人的个性化模型头皮精度Fig.6 Reconstructed scalp and individualized model scalp. (a)Contour points of head model； (b)Reconstructed scalp accuracy of the head model； (c)Individualized model scalp accuracy of the head model；(d)Contour points of human head；(e)Reconstructed scalp accuracy of human head； (f)Individualized model scalp accuracy of human head

2.2 脑表面精度

共招募10 名健康志愿者(女性5 名)，年龄为25～32 岁，均为右利手，脑部均无器质性病变，均签署了知情同意书。影像数据为全脑T1 结构像，通过Philips 3.0T 磁共振成像平台获得，体素大小为0.9 mm×0.9 mm×1.0 mm，尺寸为256×256×180。在结构像中，脑组织部分由FreeSurfer (Harvard，USA)[16]分割得到。

选用3 套影像数据作为标准模板，分别为Colin27 ( MNI，Canada)[17]、 Chinese56 ( LONI，USA)[18]和Chinese1，三者均包括全脑结构像和脑组织结构像。Colin27 的两组影像体素大小均为1 mm×1 mm×1 mm，尺寸均为181×217×181；Chinese56 的两组影像体素大小均为0.5 mm×0.5 mm×0.5 mm，尺寸分别为442×512×369 和337×387×315。Chinese1 是一套人体MRI，影像采集参数同受试者的MRI 数据。

为了保证头皮精度准确，本研究选择了3 000个轮廓点，用256 个控制点来构建个性化模型。

图7(a)为基于不同标准模板的个性化模型的脑表面精度定量结果，同时对10 例数据进行统计:基于Chinese56 的个性化模型的脑表面精度平均值最差，为(4.4±1.1)mm；基于Colin27 的个性化模型的脑表面精度平均值居中，为(3.1±0.6)mm；基于Chinese1 的个性化模型的脑表面精度平均值最佳，为(2.5±0.5)mm。图7(b)～(d)分别为额叶区、运动区和枕叶区脑表面精度的定量结果。对10 例数据进行统计:基于Chinese56 的个性化模型脑表面精度较差的区域主要集中在额叶区和运动区，平均值分别为(5.6±1.6)mm 和(4.6±1.6)mm，基于Colin27 的个性化模型脑表面精度较差的区域主要集中在枕叶区，平均值为(3.4±0.9)mm。于枕叶区的中央纵裂附近，而基于Chinese1 的个性化模型的脑表面精度较好，且不存在中央纵裂处精度差的问题。

图7 个性化模型脑表面精度定量结果。(a)全脑精度；(b)额叶区精度；(c)运动区精度；(d)枕叶区精度Fig.7 Quantitative results of individualized models’ brain surface accuracy. (a)Whole brain accuracy； (b)Frontal cortex accuracy； (c)Motor cortex accuracy； (d)Occipital cortex accuracy

表1为头颅尺寸和SCD 结果。

表1 头颅尺寸和SCDTab.1 Head size and SCD

1)头颅尺寸。可见西方人模板Colin27(209 mm×171 mm×125 mm)比中国人模板Chinese56(201 mm×172 mm×101 mm)更长更高，中国人模板Chinese56(201 mm×172 mm×101 mm)和Chinese1(200 mm×167 mm×100 mm)与测试MRI 头颅尺寸的平均值(190 mm×167 mm×100 mm)接近。

图8 个性化模型脑表面精度定性结果Fig.8 Qualitative results of individualized models’ brain surface accuracy

2)SCD 结果。Colin27 的SCD(17.1 mm)最大；Chinese56 的SCD(10.1 mm)最小，接近测试MRI 的SCD 最小值(10.0 mm)；Chinese1 的SCD 为13.9 mm，接近测试MRI 的SCD 平均值(12.8 mm)。基于Colin27 的个性化模型的SCD 平均值为11.6 mm，基于Chinese56 的个性化模型的SCD 平均值为9.6 mm，基于Chinese1 的个性化模型的SCD 平均值为10.9 mm。Colin27 的SCD 虽然较大，但是其头颅尺寸也大，经过变形得到个性化模型的SCD 接近测试MRI 的平均值。Chinese56 的头颅尺寸与测试MRI平均值更接近，但是其SCD 却远远小于测试MRI 的平均值，经过变形后SCD 会更加小于测试MRI，这是基于Chinese56 的个性化模型脑表面精度较差的主要原因。

2.3 脑靶点精度

表2 个性化模型脑靶点精度Tab.2 Targets' accuracy of individualized models

3 讨论

本设计实现了一种基于标准模板变形的个性化模型构建方法，可以为无个体影像情况下的nTMS 提供个性化影像模型。因为用于约束标准模板变形的头皮轮廓点由nTMS 的定位装置采集得到，所以个性化模型处于定位装置坐标系下，根据该模型提供的靶点位置，可以直接实现影像引导的TMS 线圈定位，避免了空间注册流程，一定程度上缩短了nTMS 准备时间。

在个性化模型中，头皮和脑表面的准确是保证刺激计划准确制定的前提。本研究通过对比不同数量头皮轮廓点和变形控制点对个性化模型头皮精度的影响，结合耗时，为个性化模型的建立提供合理参数。通过拾取2 000 个轮廓点，定位256 个控制点，可以获得(1.5±0.3)mm 头型模型的头皮精度，(2.8±0.5)mm 的真人头皮精度。虽然探针接触拾取轮廓点的方式存在精度易受人为经验影响等问题，但是从实验结果可以看出，所得到的个性化模型的头皮精度是可以接受的。在后续的工作中，可以通过激光扫描等技术，实现非接触式轮廓点获取，提高轮廓点获取精度。不同标准模板的结构差异是脑表面精度优劣的重要影响因素，本研究通过对比不同标准模板对个性化模型脑表面精度的影响，结合头颅尺寸和SCD 计算结果，为个性化模型建立提供标准模板选择方案。基于常见标准模板Colin27 的个性化模型脑表面精度达(3.1±0.6)mm，误差分布均匀，可以用于实际应用。针对中国人标准模板是否更适合中国人，笔者做了尝试后发现，基于Chinese1 的个性化模型脑表面精度优于Colin27 的相关精度，也远远优于Chinese56 的相关精度。结合头颅尺寸和SCD 计算结果，发现Chinese56 结果较差的原因是其结构局限，表现为SCD 较小，远远偏离受试者的平均值。由于Chinese1 不是标准模板，所以后期应用时还可以尝试更多的中国人模板来提高精度。

脑靶点精度是个性化模型是否适用于临床应用的重要评价标准。笔者对手部肌肉对应脑区中心点M1-HAND 和精神类疾病治疗脑区中心点DLPFC 进行了精度测量，结果发现M1-HAND 的精度为9.2 mm。nTMS 用于神经系统疾病诊断时，会在影像模型脑表面以M1-HAND 为中心放置5×5 个间隔5 mm 的靶点阵列，影像引导线圈依次刺激这些靶点，并通过观察或记录肌肉响应作为诊断依据。9.2 mm 的精度意味着靶点阵列中心错位2 个靶点，但仍在靶点阵列范围内，所以是在可以接受的误差范围内。同时，个性化模型保证了影像引导下多靶点准确重复刺激，避免了刺激线圈在头部盲目移动刺激时造成的刺激效率低下等问题。实验结果还显示，DLPFC 的精度为8.9 mm，虽然没有MRI 提供的脑靶点准确，但优于目前临床普及的传统靶点F3 的精度(16.1 mm)。同时，个性化模型保证了影像引导下靶点的准确重复刺激，能够一定程度提高疾病治疗的效果。

4 结论

本研究提出个性化模型构建方法，为无个体影像的nTMS 提供了规划靶点的头部模型数据。nTMS 系统的定位装置拾取患者头皮轮廓点，而nTMS 系统自动实现个性化模型的构建，并基于该模型中的刺激靶点，实现刺激线圈定位定向。这种个性化模型的精度能够满足nTMS 神经系统疾病诊断中的定位需求，也能够满足精神类疾病治疗中的重复定位需求，进一步提高nTMS 在临床的普及率。