CN116327167A - 一种用于脑部微组织结构测量的扩散弛豫谱成像方法 - Google Patents

Abstract

本发明公开了一种用于脑部微组织结构测量的扩散弛豫谱成像方法，包括：获取目标对象在多回波时间和多b值下的扩散磁共振数据；建立多组织脑微结构的扩散弛豫谱成像模型；通过基于球面平均技术的多参数联合优化方法获得多组织脑微结构的扩散弛豫时间和扩散敏感系数谱；将得到的弛豫时间、脑组织体积分数以及分布函数进行表征得到脑组织微观结构参数。本发明可提供组织成分特定的非T2加权的神经突内分数和自由水分数，解决了传统成分分数的TE依赖性。

Description

一种用于脑部微组织结构测量的扩散弛豫谱成像方法

技术领域

本发明属于磁共振成像和医学图像处理技术领域，具体涉及一种用于脑部微组织结构测量的扩散弛豫谱成像方法。

背景技术

脑部微组织结构测量的目的是为了实现对神经元微结构的特异性，以产生无创成像的独特临床生物标记物。由于生物组织十分复杂，微观结构建模通常需要基于特定的假设对组织结构进行模拟，以得到其感兴趣的目标指标。对于微观参数特异性成像，实现此目的的一种有效方法是增加数据的测量维度，以更好的区分组织微成分，如利用扩散和T1、T2弛豫和波谱成像之间的协同作用。相关进展包括弛豫扩散相关成像技术，已被应用于动物和离体大脑中。最近，这种多维度成像方法与微观结构建模相结合，以改善对人脑成像的组织微结构参数的估计。但是，多维度成像在提供特异性参数的同时也进一步增加了数据采集以及模型拟合的困难。

另外，在现有的扩散磁共振数据采集方案下，微观结构参数对扩散信号变化不敏感，限制了模型参数的特异性。其中，常见的多成分模型只考虑了不同成分间的扩散特性差异，而未考虑弛豫时间等其他特性的差异，而扩散磁共振数据也经常以单个TE进行采集。这导致从模型得出的定量指标会受到来自T1和T2加权的干扰。当成份之间存在T1或T2加权差异时，模型不能将成份之间的弛豫和扩散特性区分开来。在组织中神经突内和神经突外成份之间，T2的差异相比来说要小得多，但也被证明会导致在白质中一些扩散模型的导出参数具有回波时间依赖性这将减少模型衍生的定量指标的特异性，及研究结果的可解释性，特别是对于研究伴随着T2弛豫时间大幅改变的大脑发育和成熟的过程。

综上，脑组织微结构的特征在于扩散率和横向弛豫率之间的异质性。扩散磁共振成像及其高阶模型能够无创探测脑组织微结构。其中，受限、受阻和自由成分的水分子之间的信号差异可以通过高阶扩散模型来表征，从而估计体素内的细胞体、轴突纤维和自由成分的相对比例，从而用于检测大脑的发育、退化和疾病驱动的脑组织微结构改变。然而，组织成分的准确表征不仅受组织依赖的扩散系数的影响，也受到横向弛豫率的影响。明确考虑弛豫-扩散的耦合关系可能会显著改善脑组织微结构的表征。具有单一回波时间的标准扩散磁共振成像技术主要提供有关扩散率的信息，而弛豫扩散磁共振成像涉及多个回波时间和多个扩散加权强度，能够探测弛豫和扩散率之间的组织特征的耦合关系，因此提出一种用于脑部微组织结构测量的扩散弛豫谱成像方法具有重要意义。

发明内容

为了突破现有技术的挑战，本发明提供一种扩散弛豫谱成像方法，探测弛豫-扩散的耦合关系，并定量测量脑组织微结构参数，解决了传统成分分数的TE依赖性。

实现本发明目的的技术解决方案为：

一种用于脑部微组织结构测量的扩散弛豫谱成像方法，包括以下步骤：

步骤1、通过磁共振设备对脑部进行扩散加权成像扫描，根据不同b值和不同回波时间的设定，进行采样，得到多b值多回波时间的扩散磁共振成像图像；

步骤2、基于球面平均扩散弛豫谱成像模型对多b值多回波时间的扩散磁共振成像图像进行拟合得到具有脑组织特异性和b值依赖的弛豫时间，以及脑组织体积分数；

步骤3、基于扩散弛豫谱成像模型和步骤2中得到的弛豫时间对多b值多回波时间的扩散磁共振图像进行拟合得到无T2加权的具有脑组织特异性的纤维方向分布函数；

步骤4、将步骤2得到的弛豫时间、脑组织体积分数以及步骤3得到的分布函数进行表征得到脑组织微观结构参数。

与现有技术相比，本发明的有益效果为：本发明可提供组织成分特定的非T2加权的神经突内分数和自由水分数，解决了传统成分分数的TE依赖性；同时还提供了仅靠弛豫方法无法区分的神经突内、外T2弛豫时间；同时，我们通过理论推导和实验验证，显示出这些特异性参数可以从临床环境下的扫描设备中予以估计实现；本发明提高了组织参数特异性与定量特性，将为基于扩散磁共振成像的研究带来新的见解。

附图说明

图1是本发明的用于脑部微组织结构测量的扩散弛豫谱成像方法的流程图；

图2为根据本发明的方法在真实人脑数据(TE＝75,85,95,105,115,125,135ms,b＝400,800,1600,3200s/mm²)上的效果图，其中图(a)为健康样本的效果图，图(b)为神经胶质瘤患者样本的效果图；

图3为根据本发明的方法在真实人脑数据上计算的弛豫时间与b值的关系分析图；

图4为根据本发明的方法在不同神经胶质瘤患者上的效果图，(a)显示了白质病变中受限脑组织的弛豫时间效果图；(b)为传统单TE方法对神经突形态的估计图；(c)为根据本发明的方法对神经突形态的估计图；(d)为根据本发明的方法计算的轴突半径效果图；

图5为根据本发明的方法在健康人脑数据上计算的纤维方向分布，分别显示了在明确考虑弛豫和不考虑弛豫情况下的纤维方向分布的视觉比较。

具体实施方式

下面将结合附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

实施例1

本实施例提出提出一种扩散弛豫谱成像方法，获取目标对象在第一回波时间下的多b值(扩散敏感系数)扩散磁共振数据；在扩散时间和扫描重复时间均保持不变，对不同回波时间(TE)重复采集上述多b值扩散磁共振数据；根据采集的多次回波的多b值扩散磁共振数据，提出了一个统一的策略来估计(1)具有脑组织特异性的T2弛豫时间；(2)无T2加权的多尺度脑微结构组织参数；(3)无T2加权的多尺度纤维方向分布函数。

结合图1，对所述的扩散弛豫谱成像方法具体阐述为：

一种用于脑部微组织结构测量的扩散弛豫谱成像方法，包括以下步骤：

步骤1、通过磁共振设备对受试者脑部进行扩散加权成像扫描，根据不同b值和不同回波时间的设定，进行采样，得到多b值多回波时间的扩散磁共振成像图像；

步骤2、基于球面平均扩散弛豫谱成像模型对多b值多回波时间的扩散磁共振成像图像进行拟合得到具有脑组织特异性和b值依赖的弛豫时间，以及脑组织体积分数；

步骤3、基于扩散弛豫谱成像模型和步骤2中得到的弛豫时间对多b值多回波时间的扩散磁共振图像进行拟合得到无T2加权的具有脑组织特异性的纤维方向分布函数；

步骤4、将步骤2得到的弛豫时间、脑组织体积分数以及步骤3得到的分布函数进行表征得到脑组织微观结构参数。

进一步地，所述的步骤1中多b值多回波时间的扩散磁共振成像数据通过如下步骤获得：

步骤1.1、设定多b值梯度表，随着b值增加，梯度方向数量以2倍数量级增加；

步骤1.2、在扩散时间和扫描重复时间均保持固定的情况下，在受试者中进行多b值单次回波时间(TE)的扩散磁共振成像，获得各个b值对应的全采样磁共振图像；

步骤1.3、等步长增加回波时间(TE)，重复执行步骤1.2，并获取多b值多回波时间的扩散磁共振成像数据。

进一步地，所述球面平均扩散弛豫谱成像模型是扩散弛豫谱成像模型的球面平均化结果。

进一步地，所述的步骤2中球面平均扩散弛豫谱成像模型排除了纤维方向分布对于弛豫时间的影响。

进一步地，所述的步骤2中球面平均扩散弛豫谱成像模型将大脑的微观组织环境表征为脑组织内受限、受阻和自由扩散部分，并根据一个几何曲度来区分受限和受阻部分。

进一步地，所述的步骤2中球面平均扩散弛豫谱成像模型中每一个扩散部分通过一系列自旋包的球面平均信号的加权组合，并计算其对应的体积分数，从而能够探测跨越一系列扩散尺度的弛豫扩散耦合。

进一步地，所述的步骤2中球面平均扩散弛豫谱成像模型是步骤3中扩散弛豫谱成像模型的球面平均化结果；当步骤2求解得到弛豫时间后，步骤3将产生了具有组织特异性的多尺度纤维方向分布函数。

进一步地，所述的步骤4中脑组织微观结构参数通过如下步骤获得：

步骤4.1、无T2弛豫影响的方向不变性微观结构参数通过步骤2计算得到的具有组织特异性的体积分数来计算，包括轴突形态指数、轴突内体积分数、神经密度、微观各向异性指数、轴突的轴向扩散率、轴突的径向扩散率等；

步骤4.2、无T2弛豫影响的方向可变性微观结构参数通过步骤3计算得到的具有组织特异性的纤维方向分布函数来计算，包括定向受限扩散率、定向受阻扩散率等。

本发明方法称为扩散弛豫谱成像(RDSI)，该方法可以直接估计非T2加权的脑微结构组织分数比例和具有脑组织特异性的T2弛豫时间，具有以下特点：

本发明方法中的球面平均扩散弛豫谱成像模型排除了纤维方向分布对于弛豫时间的影响；

本发明方法中的球面平均扩散弛豫谱成像模型中每一个扩散部分通过一系列自旋包的球面平均信号的加权组合，并计算其对应的体积分数，从而能够探测跨越一系列扩散尺度的弛豫扩散耦合；

本发明方法中的扩散弛豫谱模型揭示了具有组织特异性的多尺度纤维方向分布；

利用本发明计算得到的具有组织特异性的体积分数，能够揭示无T2弛豫影响的方向不变性微观结构参数，包括轴突形态指数、轴突内体积分数、神经密度、微观各向异性指数、轴突的轴向扩散率、轴突的径向扩散率等；

利用本发明计算得到的具有组织特异性的纤维方向分布函数，能够揭示无T2弛豫影响的方向可变性微观结构参数，包括定向受限扩散率、定向受阻扩散率等。

实施例2

基于实施例1，对扩散弛豫谱成像模型和球面平均扩散弛豫谱成像模型具体设计包括：

所述扩散弛豫谱成像模型具体为：

扩散弛豫谱成像模型根据特定回波时间τ，扩散梯度方向g和b值采集的扩散磁共振信号S(τ，b，g)建模为

并且扩展为多组织模型为：

其中，S_r(b，h)，S_h(b，g)和S_f(b)分别表示脑组织内受限、受阻和自由扩散部分的扩散信号；梯度强度b依赖的表观弛豫率r(b)＝1/T₂(b)通过重复采集多个不同回波时间的扩散信号来估计，并表示为：

其中R(b，g，D_r)，R(b，g，D_h)，R(b，g，D_f)是与表观扩散系数D_r，D_h，D_f相关的具有脑组织特异性的响应函数，基于多尺度的响应函数与表观扩散系数，产生了具有组织特异性的多尺度纤维方向分布函数，在球谐空间中表示为球谐系数f(D_r)，f(D_h)，f(D_f)；算子A将球谐系数映射回纤维方向分布函数。

为了有效求解上模型，本发明通过球面平均技术，将上述方程简化，排除了纤维方向分布对于T2弛豫时间计算的影响，从而分步求解上述模型，所述所述球面平均扩散弛豫谱成像模型具体为：

其中w(D_r)、w(D_h)、w(D_f)是对应于脑组织内受限、受阻和自由扩散部分的体积分数，k(b，D_r)、k(b，D_h)、k(b，D_f)是响应函数R(b，g，D_r)，R(b，g，D_h)，R(b，g，D_f)的球面平均值，通过以下方程计算

本发明中，表观扩散系数D_r，D_h，D_f可以通过平行扩散系数λ_||和垂直扩散系数λ_⊥进行参数化表示在对应的集合脑组织内受限Λ_r，受阻Λ_h，自由扩散Λ_f中，并根据一个几何曲度φ来区分受限Λ_r和受阻Λ_h部分。因此，球面平均信号可以看做是自一系列自旋包的球面平均信号的加权组合，这使得我们能够探测跨越一系列扩散尺度的弛豫扩散耦合。所述平均扩散系数λ_||和垂直扩散系数λ_⊥的范围为[0e-33e-3]，所述几何曲度φ为π²/4。

为了求解该模型，本发明首先求解弛豫模式下的球面平均系数，然后从弛豫模式下的球面平均系数中分解弛豫项并求解弛豫率。最后使用估计出的弛豫率来求解纤维方向分布函数，具体如下：

我们将上述模型改写为

其中

是多次回波多b值的球面平均信号，K是响应函数的球面平均矩阵，E是弛豫矩阵，W是体积分数矩阵，/>

是Kronecker乘算子。为了求解X，我们将上述方程转换为以下问题求解

当求解得到X之后，E和W可以通过最小化以下约束非线性多变量问题来进行分解

当求解得到E之后，纤维方向分布函数的估计模型可以重写为一个严格凸的二次规划问题来求解

当求解得到W之后，轴突形态指数(包括平均轴突半径(Mean NR)、内部偏差指数(Std.NR)和相对轴突半径(Cov.NR))可以通过以下方式计算：

其中∈＞0是脉冲尺度，并且仅依赖扩散梯度的脉冲宽度δ和扩散时间Δ，

所述的受试者为人或实验动物，实验动物为鼠或猴或狗或猪。

通过理论推导和实验验证得到图2～图5，图2为根据本发明的方法在真实人脑数据(TE＝75,85,95,105,115,125,135ms,b＝400,800,1600,3200s/mm²)上的效果图，左边为健康样本，右边为神经胶质瘤患者样本，图3为根据本发明的方法在真实人脑数据上计算的弛豫时间与b值的关系分析图，显示出了在较高b值下，受限的脑组织中更多的体素的弛豫时间在100-200ms之间，尤其是高级别神经胶质瘤，图4为根据本发明的方法在不同神经胶质瘤患者上的效果图，(a)显示了白质病变中受限脑组织的弛豫时间，表明神经胶质瘤的弛豫时间比正常白质组织长；(b-c)显示本发明(c)相比于传统单TE方法(b)对神经突形态的估计更加敏感；(c-d)为根据本发明的方法计算的轴突半径，显示出了本发明对肿瘤范围的描绘和肿瘤内微环境的标准非常敏感，图5为根据本发明的方法在真是健康人脑数据上计算的纤维方向分布，分别显示了在明确考虑弛豫和不考虑弛豫情况下的纤维方向分布的视觉比较，通过图2-图5，显示出这些特异性参数可以从临床环境下的扫描设备中予以估计实现，验证本发明解决了传统成分分数的TE依赖性；同时还提供了仅靠弛豫方法无法区分的神经突内、外T2弛豫时间。

Claims (10)

1.一种用于脑部微组织结构测量的扩散弛豫谱成像方法，其特征在于，包括以下步骤：

步骤1、通过磁共振设备对脑部进行扩散加权成像扫描，根据不同b值和不同回波时间的设定，进行采样，得到多b值多回波时间的扩散磁共振成像图像；

步骤2、基于球面平均扩散弛豫谱成像模型对多b值多回波时间的扩散磁共振成像图像进行拟合得到具有脑组织特异性和b值依赖的弛豫时间，以及脑组织体积分数；

步骤3、基于扩散弛豫谱成像模型和步骤2中得到的弛豫时间对多b值多回波时间的扩散磁共振图像进行拟合得到无T2加权的具有脑组织特异性的纤维方向分布函数；

步骤4、将步骤2得到的弛豫时间、脑组织体积分数以及步骤3得到的分布函数进行表征得到脑组织微观结构参数。

2.根据权利要求1所述的一种用于脑部微组织结构测量的扩散弛豫谱成像方法，其特征在于，所述步骤1具体包括：

步骤1.1、设定多b值梯度表，随着b值增加，梯度方向数量以2倍数量级增加；

步骤1.2、在扩散时间和扫描重复时间均保持固定的情况下，进行多b值单次回波时间的扩散磁共振成像，获得各个b值对应的全采样磁共振图像；

步骤1.3、等步长增加回波时间，重复执行步骤1.2，获取多b值多回波时间的扩散磁共振成像数据。

3.根据权利要求1所述的一种用于脑部微组织结构测量的扩散弛豫谱成像方法，其特征在于，所述球面平均扩散弛豫谱成像模型是扩散弛豫谱成像模型的球面平均化结果。

4.根据权利要求3所述的一种用于脑部微组织结构测量的扩散弛豫谱成像方法，其特征在于，所述球面平均扩散弛豫谱成像模型排除纤维方向分布对于弛豫时间的影响，将大脑的微观组织环境表征为脑组织内受限、受阻和自由扩散部分，并根据一个几何曲度来区分受限和受阻部分；每一个扩散部分通过一系列自旋包的球面平均信号的加权组合，并计算其对应的体积分数，探测跨越一系列扩散尺度的弛豫扩散耦合。

5.根据权利要求4所述的一种用于脑部微组织结构测量的扩散弛豫谱成像方法，其特征在于，所述扩散弛豫谱成像模型具体为：

扩散弛豫谱成像模型根据特定回波时间τ，扩散梯度方向g和b值采集的扩散磁共振信号S(τ，b，g)建模为

并且扩展为多组织模型为：

其中，S_r(b，g)，S_h(b，g)和S_f(b)分别表示脑组织内受限、受阻和自由扩散部分的扩散信号；梯度强度b依赖的表观弛豫率r(b)＝1/T₂(b)通过重复采集多个不同回波时间的扩散信号来估计，并表示为：

其中R(b，g，D_r)，R(b，g，D_h)，R(b，g，D_f)是与表观扩散系数D_r，D_h，D_f相关的具有脑组织特异性的响应函数，基于多尺度的响应函数与表观扩散系数，产生了具有组织特异性的多尺度纤维方向分布函数，在球谐空间中表示为球谐系数f(D_r)，f(D_h)，f(D_f)；算子

将球谐系数映射回纤维方向分布函数。

6.根据权利要求5所述的一种用于脑部微组织结构测量的扩散弛豫谱成像方法，其特征在于，所述球面平均扩散弛豫谱成像模型具体为：

其中w(D_r)，w(D_h)，w(D_f)是对应于脑组织内受限、受阻和自由扩散部分的体积分数，k(b，D_r)，k(b，D_h)，k(b，D_f)是响应函数R(b，g，D_r)，R(b，g，D_h)，R(b，g，D_f)的球面平均值，具体通过以下方程计算

7.根据权利要求6所述的一种用于脑部微组织结构测量的扩散弛豫谱成像方法，其特征在于，所述表观扩散系数D_r，D_h，D_f通过平行扩散系数λ_||和垂直扩散系数λ_⊥进行参数化表示在对应的集合脑组织内受限Λ_r，受阻Λ_h，自由扩散Λ_f中，并根据一个几何曲度φ来区分受限Λ_r和受阻Λ_h部分。

8.根据权利要求7所述的一种用于脑部微组织结构测量的扩散弛豫谱成像方法，其特征在于，所述平均扩散系数λ_||和垂直扩散系数λ_⊥的范围为[0e-33e-3]，所述几何曲度φ为π²/4。

9.根据权利要求6所述的一种用于脑部微组织结构测量的扩散弛豫谱成像方法，其特征在于，所述扩散弛豫谱成像模型和球面平均扩散弛豫谱成像模型通过以下方式求解：

首先求解扩散弛豫谱成像模型的球面平均系数，然后从球面平均系数中分解弛豫项并求解弛豫率，最后使用估计出的弛豫率来求解纤维方向分布函数，具体如下：

球面平均扩散弛豫谱成像模型模型方程改写为：

其中

是多次回波多b值的球面平均信号，K是响应函数的球面平均矩阵，E是弛豫矩阵，W是体积分数矩阵，/>

是Kronecker乘算子，为了求解X，将上述方程转换为以下问题求解：

当求解得到X之后，E和W可以通过最小化以下约束非线性多变量问题来进行分解：

当求解得到E之后，纤维方向分布函数的估计模型重写为一个严格凸的二次规划问题来求解：

10.根据权利要求9所述的一种用于脑部微组织结构测量的扩散弛豫谱成像方法，其特征在于，所述步骤4具体包括：

步骤4.1、无T2弛豫影响的方向不变性微观结构参数通过步骤2计算得到的具有组织特异性的体积分数来计算，包括轴突形态指数、轴突内体积分数、神经密度、微观各向异性指数、轴突的轴向扩散率和轴突的径向扩散率；

步骤4.2、无T2弛豫影响的方向可变性微观结构参数通过步骤3计算得到的具有组织特异性的纤维方向分布函数来计算；

所述轴突形态指数包括平均轴突半径Mean NR、内部偏差指数Std.NR和相对轴突半径Cov.NR，具体为：

其中∈＞0是脉冲尺度，并且仅依赖扩散梯度的脉冲宽度δ和扩散时间Δ，

Std.NR＝std((∈(δ，Δ)w_Drλ_||λ_⊥)^1/4)

Images