磁共振扩散张量成像中扩散敏感梯度磁场方向分布方案的研究进展*

刘良友 高嵩 † 李莎 李兆同 夏一帆

1) (北京大学医学部医学技术研究院, 北京 100191)

2) (北京大学医学技术研究发展中心, 北京 100191)

磁共振扩散张量成像可以定量无创研究人体内水分子在三维空间中的各向异性扩散规律, 进而获取重要的病理及生理信息.为了得到水分子各向异性扩散信息, 需要按照一定的方案依次施加不同方向的扩散敏感梯度磁场, 测量水分子在这些方向上的扩散系数用以估算扩散张量.扩散张量成像测量结果的准确程度受梯度磁场方向分布方案的影响, 本文对扩散敏感梯度磁场方向分布方案进行综述, 包括完全随机方案、启发式方案、规则多面体式方案和数值优化方案等, 分析这些方案的优势与局限性, 并提出需进一步研究的问题.

综述

1 引 言

在均匀介质中水分子的随机扩散现象表现为各向同性, 其扩散程度可以用扩散系数D表示.在复杂的生物组织中, 水分子的扩散现象受组织结构的影响, 表现为各向异性扩散, 扩散程度用扩散张量D表示[1].扩散张量成像(diffusion tensor imaging,DTI)[2,3]技术是在标准磁共振成像(magnetic reso−nance imaging, MRI)脉冲序列中加入不同方向的扩散敏感梯度磁场(diffusion sensitive gradient,DSG), 测量不同方向上水分子的扩散情况, 最终获得人体内大量的生理和病理信息.DTI技术提供了一种在微观结构尺度上以非侵入性的方法来表征人体组织的结构和功能特性, 已经成为了解正常大脑组织结构[4,5]以及神经和精神性疾病病变过程的重要方式[6,7].

为了使DTI结果更加准确, 需要施加数量众多、线性无关且空间均匀分布的DSG.近年来为了进一步获取人体内复杂的微观结构, 提出了高角分辨率扩散成像(high angular resolution diffusion imaging, HARDI)[8]技术, 对DSG的数量和方向分布的均匀性提出了更高的要求.本文从DTI基本原理出发, 介绍了多种不同的DSG方向分布方案及特点, 同时提出DSG方向分布方案应需进一步优化的问题.

2 DTI基本原理

2.1 Bloch-Torrey方程

为了描述扩散现象与磁化强度M随时间的变化, Torrey在Bloch方程中加入扩散项得到关于扩散磁共振的方程[9]：

图1 Stejskal−Tanner序列, 其中两个同等的DSG脉冲置于180°RF脉冲两侧, 强度为G, 宽度为δ, 时间间隔为ΔFig.1.Stejskal−Tanner Scheme：Two diffusion sensitive gradients inserted before and after 180° RF refocusing pulse.G, amplitude;δ, duration of the DSG; Δ, time between the two sensitive gradient lobes.

其中B为磁场强度; γ为旋磁比; Mx, My, Mz是磁化矢量的3个分量; M0是平衡态下的磁化强度;T1, T2是纵向和横向弛豫时间.(1)式右端的前三项是经典的Bloch方程, 当施加 90°射频(radio frequency, RF)脉冲时, (1)式在横向平面上的磁化强度Mxy的解为

其中b为描述施加的DSG强度(G)、持续时间等因素的扩散梯度因子, 单位是s/mm2, 值越大表示扩散权重越大, 对水分子扩散更加敏感, 可以表示为

Stejskal−Tanner序列是利用水分子信号衰减的特性来测量水分子的扩散[10,11], 它是将一对大小和持续时间完全相同的DSG沿着某一方向置于180°脉冲的两侧(图1), 其信号衰减的大小可以由(1)—(3)式求解[12]：

其中S, S0表示施加和未施加DSG的回波信号值;Δ为两个脉冲之间的时间间隔; δ为脉冲的持续时间; D单位是mm2/s, 值越大表示水分子扩散运动的能力越强.

2.2 扩散张量成像

引入扩散张量D描述人体内水分子沿各向异性扩散规律[1], 即：

扩散张量矩阵为对称阵, 即Dij= Dji, i, j = x, y,z, 其中对角线元素Dxx, Dyy, Dzz和非对角线元素Dxy, Dxz, Dyz是扩散张量的6个独立分量, 表示水分子在不同方向上的扩散系数[1].

由(4)式可知, 对于各向异性扩散, 可以得到DTI图像在不同DSG下S与S0的关系：

其中gi= (gxigyigzi)是归一化矢量, 表示3个正交方向上的DSG强度, 等式左端为表观扩散系数(apparent diffusion coefficient, ADC), (6)式展开为

其中i代表DSG方向的数量, ADCi表示第i个方向的表观扩散系数, 张量D有6个独立分量, 因此需要在至少6个不共线的方向上施加DSG, 并代入方程(7)得到张量D.为了得到更加准确的结果, 并用于计算扩散张量D, 同时减少噪声对结果的影响, 通常会在众多的空间方向上施加DSG[13],而HARDI技术对DSG方向数量要求则更多.

通过对扩散张量矩阵对角化[14]可以得到扩散张量的特征值及特征矢量.特征值和特征矢量是旋转不变量, 且每个特征值对应一个特征矢量, 特征矢量之间相互垂直, 因此对于扩散的三维特性可以使用椭球体来描述, 椭球体的三个轴由扩散张量的三个特征值 λ1, λ2, λ3和特征矢量 ε1, ε2, ε3确定[15].同时可以利用特征值定义多种各向异性参数指标[16,17], 包括各向异性分数(fractional anisotropy,FA)、相对各向异性(relative anisotropy, RA)和容积比(volume ratio, VR)等, 还包括评价整体扩散效应指数的张量迹(trace)、平均扩散系数(medium diffusion, MD)等.

3 DSG方向分布方案

在DTI中, 通常在半个或整个球面上沿着至少6个线性无关的方向施加DSG, 并且要求DSG方向分布尽可能均匀, 才能获得较为准确的扩散张量D[18,19].一个空间均匀的DSG方向分布方案对于提高张量估算的精确性有着重要作用[20−22], 在一定程度上增加DSG的数量也有利于提高图像的信噪比(signal noise ratio, SNR)[23].HARDI技术对DSG的分布方案有更高要求, 但随着DSG数量的增多, 图像的采集时间也在不断增加, 因此大量研究提出了多种DSG方向分布方案[21,22].

3.1 随机分布方案

一般来说, 对于感兴趣区内组织结构的方向分布并不是已知的, 通常假设最佳的DSG方向是均匀分布在单位球面上[21], 且所有的DSG方向之间不共线.随机分布的DSG方向分布方案[24]即是在单位球表面随机分布一定数量的DSG方向, 按照随机方案产生的次序依次施加DSG, 最终获得随机分布方案的扩散信号.该方案虽然可以获得任意数量的DSG分布方向, 但是分布的均匀性较差(图2), 不能满足DSG方向分布尽可能均匀的要求, 并且随机分布的方案存在无法复现的问题, 因此随机分布是一种较差的DSG方向分布方案.

3.2 启发式分布方案

启发式DSG方向分布方案是基于MRI系统的x, y, z梯度方向来定义的一组基本方向, 对应着边长等于2, 中心点在(0, 0, 0)的立方体的三个面中心, 可以生成最多13个非共线的DSG方向(gx, gy, gz), 分别对应于面中心

边中心线

立方体对角线

图2 随机分布60个方向DSG方向分布方案Fig.2.DSG encoding scheme with random distribution in 60 directions.

启发式DSG方向分布方案包括：1)正交DSG方向分布方案, 由三个正交轴和三个边中心线组合成6个DSG分布方向, 即(G1+ G2), 称为金字塔锥体编码, 常用于化学屏蔽张量谱学中[25]; 2)倾斜双DSG方向分布方案, 由Basser和Pierpaoli[26]提出将所有的边中心线产生 6个DSG方向进行组合, 即(G2+ G3), 同样可以应用于化学屏蔽张量谱学中; 3)正交轴和四面体组合DSG方向分布方案, 由三个正交轴梯度和立方体对角线[27]梯度组合成7个方向, 即(G1+ G4), 不仅可以用于DTI[28]还可以用于热膨胀的张量测量[29];4)十面体DSG方向分布方案, Skare 和Nordell[30]将边中心线和立方体对角线组合成10个DSG方向, 即(G2+ G3+ G4); 5)完全启发式DSG方向分布方案, 将所有的面中心、边中心和对角线组合所产生的13个方向的DSG方向分布方案, 即(G1+ G2+ G3+ G4), 如图3所示.

图3 完全启发式13个方向DSG方向分布方案Fig.3.DSG encoding scheme with heuristic distribution in 13 directions.

启发式DSG方向分布方案由立方体对应的不同方向组合而成, 生成5种可供选择的子方案, 能够产生较为均匀的空间分布, 应用于张量测量的同时还在化学领域有较为广泛的应用, 但是启发式DSG方向分布方案所产生的方向数量最多仅为13个, 并且不能生成连续的方向数, 对于张量的精确估算来说, 该方案方向数量则略显不足.

3.3 几何多面体分布方案

利用高对称性的几何多面体生成DSG方向分布方案.几何多面体是由等边三角形、正方形或者正五边形组成的立体几何, 最常用的几何多面体包括四面体、立方体、八面体、十二面体和二十面体等.几何多面体的顶点数与面数之和比边数多出2个, 立方体的顶点与八面体的面彼此对应, 称为对偶, 十二面体与二十面体也具有类似的性质, 在这四种正几何多面体中, 从一个顶点通过中心绘制的线经过另一个顶点, 从一个面通过中心绘制的线经过另一个面, 从一条边通过中心绘制的线经过另一条边.因此它们的顶点、面、边可以定义一半的不共线方向, 如3, 4, 6, 10, 15个方向[31].而正四面体则是从一个顶点穿过中心经过另一个面的中心, 因此可以生成3或4个DSG方向, 称为四面体梯度方案[27], 且3和4个DSG方向可以组合生成7个DSG方向[32,33], 这就对应了四面体所有的顶点和边或者对应立方体、八面体的所有顶点和面, 并且可以调整其中一些方向的极性使其分布更加均匀[31].

立方体和八面体的十二条边可以定义6个均匀分布的DSG方向用于计算扩散张量值[21,26,34−36].同样立方体的六个面和十二条边可以组成9个DSG方向, 并且立方体的十二条边可以被二十面体的十二个顶点取代产生新的9个DSG方向, 后者的空间均匀性更好[37].13个DSG方向可以来自于立方体所有的面、顶点、边的组合, 其中任何一组的方向也可以由其他多面体DSG方向所取代或改变方向极性产生新的更加均匀的DSG方向分布.

二十面体包含的十二个顶点可以得到6个DSG方向分布, 二十个面可以得到10个DSG方向[21,31,36], 三十条边可以得到15个DSG方向[31,34].若6, 10, 15个梯度为单一方向, 则可以相互组合成为数量更多的分布均匀的DSG方案, 如6, 10,15, 16, 21, 25, 31等方向, 在这些方向组合中, 有的需要对一个或多个方向进行极性调整来获取最佳的DSG分布.对于一些其他数量的方向如12个方向, 则可以通过沿着x, y或z轴旋转90°来生成,或者通过删除二十面体15个方向中的立方体面产生的3个方向来获得.其中二十面体的31个方向如图4所示.

图4 二十面体31个方向DSG方向分布方案Fig.4.DSG encoding scheme with icosahedron distribution in 31 directions.

将二十面体每条边的中点作为新的顶点, 并将每个新的顶点相连, 来生成Ne= 5n2+ 1个方向,其中 n = 1, 2, 3, ··[38,39], 该方案可以产生 6, 21,46, 81, 126等方向数量.几何多面体能够生成较为均匀的DSG方向分布方案, 但是方向数量不多且不连续.

3.4 DISCOBALL分布方案

与DSG方向均匀的分布在球面上相比, 将DSG方向均匀的分布在圆环上更为简单.Stirnberg等[40]将三维的球面分布简化为二维的圆环分布,首先利用恒定的天顶角增量将球面分成均匀的片层, 然后将片层以恒定的方位角增量进行切割.每个片层的DSG方向的分布数量为：2Nsin(θ), 其中N为预定义的片层数, θ为天顶角增量, 同时将数量的结果四舍五入, 并且可以由

根据预定义的DSG方向数量Ntotal生成片层数来计算各片层上的DSG方向数, 最终得到整个球面的DSG方向的分布, 如图5所示.

该方案可以产生均匀分布的DSG方向, 但是子集分布的均匀性较差, 对于不完整扫描的情况容易对张量产生错误的估算结果.

图5 DISCOBALL 30个方向DSG方向分布方案Fig.5.DSG encoding scheme with DISCOBALL distribu−tion in 30 directions.

3.5 球面螺旋分布方案

基于单位球面恒定角速度采样的DSG方向分布方案会导致赤道位置DSG方向分布的密度较低, 而靠近两极处则密度较高, 对于尽可能均匀分布的空间DSG方向来说, 恒定角速度的DSG分布方案显然无法满足要求.Wong和Roos[41]提出球面螺旋的分布方案, 该方法是利用求解球面的方位角及顶角的速率变化生成恒定速率的螺旋路径,并扫过单位球面来产生均匀分布的DSG方向.DSG方向(gx, gy, gz)经归一化后可以写为

其中 N 为 DSG 方向数量, n = 1, 2, ··, N.此方法适用于DSG方向数量较多的分布方案, 可以提供覆盖单位球面的方向螺旋, 由于相邻点之间的距离相等, 因此DSG方向分布近似均匀[21,31], 如图6所示.该方案的不足之处是容易在两极附近产生空洞, 且DSG方向子集分布的均匀性较差.

图6 球面螺旋分布60个方向DSG方向分布方案Fig.6.DSG encoding scheme with spherical spiral distribu−tion in 60 directions.

3.6 Jones方案

Jones方案[20,21,31,42−44], 又称为静电力排斥或者最小作用力算法的DSG方向分布方案, 其建立在Conturo等[27]的平衡化学sp3杂化轨道中静电排斥模型的基础之上, 是目前应用范围最为广泛的DSG方向分布方案之一.假设建立一种模型,在该模型中DSG方向都是经过球体中心的线, 在线与球体表面相交处的两个点上放置单位点电荷,不断改变线的方向, 使产生的每个DSG方向都用一对点来表示, 在相反方向上的DSG同样可以用正方向上的DSG进行扩散衰减的测量.根据库仑定律, 一对点电荷之间的排斥力与点电荷之间距离的平方成反比, 因此, 用于在三维空间中均匀排列DSG方向的算法被用来优化这些正交梯度的方向,直到所有可能的电荷对之间的作用力总和最小.该方案生成的方向见图7(a).

该方法可以得到数量众多方向连续且分布均匀的DSG方向分布方案, 并且对于各向异性程度较高的组织, 众多不共线的DSG方向可以提供更为稳健的张量估算.但是由于初始方向的随机性,该方案无法复现, 单位点电荷之间最小静电力斥力算法也是一种迭代耗时的数值计算方法, 且很难得到一个简单的解析式, 其最大的不足是该方案的子集分布的均匀性不高, 无法满足一些临床要求.

图7 两种Jones 60个方向DSG方向分布方案 (a) Jones方案; (b)排序的Jones方案Fig.7.DSG encoding scheme with Jones (a) and Ordered Jones (b) in 60 directions.

Dubois等[22,45]提出把图像采集过程中位置接近的DSG方向进行分离, 以此来提高部分采集时张量估算的精确性, 称为加权的Jones方案, 主要方法是在能量公式中引入权重的概念：

其中Eij为两个不同方向i, j间相互作用的能量,ωij为相互作用权重, 在一个完全各向同性的方向分布方案中, ωij= 1.

Dubois等[22,45]提出三种不同的加权方案, 第一种产生DSG方向的子集, 使每个子集包含6个方向, 如果方向i, j在同一个子集中则ωij= 1, 否则, ωij= 0.6; 第二种方案是相邻子集间的权重为ωij= 0.8; 第三种方案是不产生子集, 但是权重因子随着序列的采集顺序不断下降, 如果|i — j| ≤10, 则 ωij= 1, 如果|i — j| ＞ 10, 则 ωij= |i — j|—α,α为常数值.

三种加权方案对应三种不同的采样方式, 第一种采样方式由不同方向的子集组成, 每个子集包含相同的方向数量(6个), 属于同一个子集的方向之间是完全相互作用的, 即权重因子等于1, 而属于不同子集的方向相互作用系数大小为0.6; 第二种采样方式原理与第一种方式类似, 不同的是方向在不同子集间的相互作用都是不相同的, 并且后一个子集应包括前一个子集中的方向, 各个子集的时间间隔越长, 相互作用越小, 权重因子也越小, 这就要求后一个子集必须在之前的集合基础上提供更加准确的空间方向信息; 第三种采样方式是完成一次完整的DSG方向采样, 不生成子集, 因此要求最少有6个DSG方向均匀分布, 在采样模型中,某一个方向与它时间上相距最近的10个方向完全相互作用, 与其他方向的相互作用随着时间的增加呈指数模型规律递减.该方案采用的是模拟退火法来最小化全局能量.

加权的Jones方案可以解决原Jones方案中子集方向分布不均匀的问题, 它的劣势在于虽然采用最小化加权能量公式产生DSG方向分布方案,但该方案却无法满足未加权的能量最小化, 因此其整体的DSG方向分布的均匀性不如Jones方案,另外由于能量公式中权重的变大会导致局部最小值数量的增加, 提高了求解全局最小值的难度.

Cook等[46]提出将所有的静电点集分为相同大小的子集, 并且使每个子集中的能量最小化以达到分布的各向同性, 具体来说就是最小化能量函数：

当 i, j在同一个子集中时, δij= 1, 否则 δij= 0.该方案的目的就是建立各向同性的子集, 利用每个子集来更加精确的估算扩散张量, 得到用于运动校正的张量信息.与未经排序的方向集相比, 提高了部分扫描的质量, 但是由于该方案是在独立处理每个子集方向, 导致对于前P个方向的优化无法达到最佳, 比如对于第一个子集和第二个子集的前一半组成的扫描序列可能会呈现各向异性较大的结果.

对于上述方法的不足, Cook等[43]再次提出一种新的方案来优化全局的采集顺序, 同时在不影响整体扫描质量的情况下, 提高了部分扫描的结果.该方案称为排序的Jones方案, 与上述方案不同的是, 子集的划分为嵌套式, 而后者将方向划分为非重叠的子集.目标是同时最小化所有子集的静电能量来求解最优排序：

EP为前P个方向子集的静电能量, 当P个方向各向同性分布时, 其静电能近似正比于P2, 因此归一化因子P—2使得每个子集的对目标函数具有相近的作用, 求解最优排序采用模拟退火算法, 方向如图7(b)所示[47].

4 总结与展望

由DTI原理可知, DSG分布越均匀, 对扩散张量的估算越精确, 同时对组织结构和纤维束走行分布的描述越真实.随机分布的DSG方向分布方案虽然可以产生数量较多且连续的DSG方向, 但DSG方向分布不均匀, 临床上很少采用这种方案;启发式方案弥补了随机方案方向分布不均匀的缺点, 但其产生的方向数量较少且方向数不连续, 应用范围较窄; 球面螺旋和DISCOBALL分布方案可以产生分布均匀且数量连续的DSG, 但是子集分布的均匀性较差, 应用较少; 几何多面体方案,能够产生均匀的DSG方向空间分布, 但是采样的数量和连续性受限; 目前应用最为广泛的是Jones DSG方向分布方案, 它能够产生任意数量且分布均匀的DSG方向, 并且加权和排序的Jones方案也解决了子集分布不均匀的问题, 可以应对临床采集过程中由于患者不配合或者不自主运动导致的数据集损坏的状况, 但是由于Jones方案的随机性, 造成梯度方案无法复现, 同时该方案采用迭代的数值计算方法, 导致计算过程复杂且耗时.

单纤维取向扩散张量模型假设每个体素中只能有一个主扩散方向, 但当纤维束出现交叉、分支、汇聚等情况时, 该模型则显得过于简单, 同时对于张量的估算及复杂的神经元微结构精度的判断变得不确定.HARDI[48−50]技术则沿着更多方向施加DSG, 因此它能够更准确解决纤维束交叉等问题, 其采样是分布于q空间的单球壳(single−shell)或者多球壳(multi−shell)上, 相较于 q空间的大量采样, HARDI减少了采样的数量, 降低了采集时间, 同时对b值的要求较低, 也得到了相对较好的SNR.但由于单球壳的HARDI (sHARDI)采样重建出的扩散概率分布函数不存在q空间的径向信息, 因此为了能够从径向信息中获取更多的扩散方向信息, 有研究提出多球壳HARDI (mHARDI)采样技术[51], 在每个球壳上DSG方向同样是均匀分布, 其在解决纤维束交叉等问题上的结果要优于sHARDI采样技术.

基于以上所有DSG方向分布方案所存在的不足和需要进一步探索的问题, 可考虑将黄金分割法应用于球面DSG方向分布方案, 并采用统计学方法验证该方案方向分布的均匀性, 进而得到DSG的均匀性对DTI结果的影响, 并且该方案既有望能够满足HARDI的要求, 同时在面对临床中可能碰到的数据集损坏问题时也可以得到较准确的张量估算结果.