CN105980876A - 用于定量样品中的各向同性弥散和/或各向异性弥散的方法 - Google Patents

Abstract

根据本发明构思的一方面，提供了一种用于定量样品中的各向同性弥散和/或各向异性弥散的方法，所述方法包括：利用弥散编码磁梯度脉冲序列G_i＝1…m对所述样品进行弥散加权磁共振测量，其中每个磁梯度脉冲序列G_i被生成为使得用于所述磁梯度脉冲序列G_i的弥散编码张量b_i具有一至三个非零本征值，其中q_i(t)与公式成正比，τ是回波时间。所述方法还包括收集通过对所述样品进行所述测量得到的表示磁共振回波信号的数据，其中所述数据的至少一个子组表示利用引起各向异性弥散加权的一组磁梯度脉冲序列得到的回波信号，并且其中用于所述磁梯度脉冲序列组中的每个梯度脉冲序列的弥散编码张量具有三个非零本征值，所述三个本征值中的至少一个与另外两个本征值不同。所述方法还包括利用所述数据计算各向同性弥散程度和/或各向异性弥散程度。

Description

用于定量样品中的各向同性弥散和/或各向异性弥散的方法

技术领域

本发明构思涉及用于定量样品中的各向同性弥散和/或各向异性弥散的方法。

背景技术

各种各样的多孔材料(从溶致液晶[1]到脑组织[2])包含在介观尺度上具有不同大小、形状和取向度的各向异性孔。材料的完整表征要求对所有这些参数进行估计，但遗憾地是，当使用基于Stejskal-Tanner序列[3]与两个磁场梯度脉冲的常规弥散MRI法时，这些参数对检测到的MRI(磁共振成像)信号的影响糟糕地纠缠在一起。在下文中该序列可以称作单脉冲场梯度(sPFG)序列或实验。

在弥散MRI(dMRI)中，图像的每个体素(其通常可为毫米级的)包含关于水在微米尺度上的平移信息[15]。将sPFG应用于弥散张量成像(DTI)中，能够定量平均弥散率(MD，也即表观弥散系数，ADC)和弥散各向异性(部分各向异性，FA)。虽然基于sPFG的DTI测量对细胞结构中的变化非常敏感，但sPFG通常仅能对高度组织化的白质束提供稳健估计。在较为无序的组织中，它几乎不能对这类变化的性质提供洞察，从而导致常见的误判。例如，FA的变化被认为表示白质完整性，但是很多因素(细胞死亡、水肿、神经胶质增生、炎症、髓鞘形成中的变化、交叉纤维的连接性增加、细胞外水分或细胞内水分增加等等)都会引起FA的变化。测量的有限特异性(例如FA和MD)阻碍了我们将测量结果与神经病理学或与局部解剖变化(例如连接性差异)关联起来[24,25,26,27]。与sPFG不同，通过提供关于体素内细胞形状、大小和膜特性的分布信息，非常规的dMRI序列能够开始在大脑宏观尺度和微观尺度之间建立联系。

依据化学位移与弥散各向异性张量之间的形式模拟，已经表明固态NMR(核磁共振)技术(例如，“魔角旋转”)能够适用于弥散MRI[4]。在其最简单的形式中，q-向量的魔角旋转能够用于评估各向同性弥散率的分布，而不受各向异性的混杂影响。

WO 2013/165312公开了如何通过移相向量q(t)的连续或离散调制实现弥散加权回波信号衰减的各向同性弥散加权，从而例如通过利用魔角旋转使得对回波信号的各向异性影响最小化。WO 2013/165313公开了通过分析利用两种不同梯度调制方案所得到的回波衰减曲线定量微观弥散各向异性和/或平均弥散率的方法，其中一个梯度调制方案基于各向同性弥散加权，另一个梯度调制方案基于非各向同性弥散加权。WO 2013/165313公开了可以通过例如利用单脉冲梯度自旋回波(PGSE)实现非各向同性弥散加权。

尽管这些现有技术方法能实现对回波信号衰减的各向同性贡献和各向异性贡献的分离以及对(尤其是)微观部分各向异性的定量，但在一些情况下就用于引起弥散加权的梯度调制方案来说具有更大的自由度，且仍能够分析和定量诸如微观弥散各向异性和/或平均弥散的微观结构性质(例如为了利用弥散波谱进行组织表征的目的)是可取的。例如，各向同性弥散编码在某些情况下对硬件的转换速率和最大量值提出了较高要求，较为老式且低廉的设备很难满足这一要求。

发明内容

本发明构思的一个目的是提供用于定量样品中的各向同性弥散和/或各向异性弥散的方法，该方法不需要使用引起各向同性弥散编码的弥散编码磁梯度脉冲序列。通过以下发明内容可以理解其它目的。

根据本发明构思的一个方面，提供了用于定量样品中的各向同性弥散和/或各向异性弥散的方法，该方法包括：

利用弥散编码磁梯度脉冲序列G_i＝1…m对样品进行弥散加权磁共振测量，其中每个磁梯度脉冲序列G_i被生成为使得用于磁梯度脉冲序列G_i的弥散编码张量b_i具有一至三个非零本征值，其中q_i(t)与成正比，τ是回波时间。

该方法还包括收集通过对样品进行所述测量得到的表示磁共振回波信号的数据，其中所述数据的至少一个子组表示利用引起各向异性弥散加权的一组磁梯度脉冲序列得到的回波信号，并且其中用于所述磁梯度脉冲序列组中的每个梯度脉冲序列的弥散编码张量具有三个非零本征值，所述三个本征值中的至少一个与另外两个本征值不同。该方法还包括利用所述数据计算各向同性弥散程度和/或各向异性弥散程度。

除了别的以外，本发明的方法基于这样的认识：利用磁梯度脉冲序列(或更短的“脉冲序列”)进行弥散编码，使得弥散编码张量具有三个非零本征值，其中至少一个与另外两个不同，并引起各向异性弥散加权，从而使得控制样品材料中弥散各向异性对回波信号的影响成为可能。如下文将更详细地描述，该方法能够准确表征微观弥散特性(例如各向同性和各向异性弥散)，尤其是样品内部微区室的微观弥散特性，该内部微区室比磁共振测量的空间分辨率更小。此外，这种表征能够实现，而不需要如现有技术一样依赖于对各向同性弥散加权的使用(该情形可通过具有三个相等的非零本征值实现)。这可使得能在更大范围的设备上进行准确的弥散测量。

根据本发明方法，每个磁梯度脉冲序列G_i被生成为使得用于磁梯度脉冲序列G_i的弥散编码张量b_i具有1至3个非零本征值。换句话说，每个磁梯度脉冲序列G_i被生成为使得存在用于脉冲序列G_i的弥散编码张量表示b_i，其具有1至3个非零本征值。类似地，对于引起各向异性弥散加权的上述磁梯度脉冲序列组中的每个磁梯度脉冲序列来说，存在弥散编码张量表示，其具有3个非零本征值，其中至少一个与另外两个本征值不同。引起各向异性弥散加权的磁梯度脉冲序列组可以形成至少一个弥散编码磁梯度脉冲序列G_i＝1…m子组。用于引起各向异性弥散加权的每个所述梯度脉冲序列的弥散编码的所述至少一个本征值可以有利地与另外两个本征值中的任一个相差至少5％，并且甚至更优选地相差至少10％。这可以确保样品中有足够程度的各向异性弥散加权，有利于后续的计算并降低了硬件要求。

所述数据子组可表示从相同的样品部分获得的回波信号，该部分包括多个分体积(partial volume)，各分体积呈现不同程度的各向同性弥散或不同程度和/或取向的各向异性弥散，其中计算各向同性弥散程度和/或各向异性弥散程度可以包括针对所述分体积中的至少一个计算各向同性弥散程度的估计值和/或各向异性弥散程度的估计值。

特别地，所述部分可以具有与弥散加权磁共振测量的空间分辨率匹配的空间延伸度。因而，每一个分体积可以具有小于该空间分辨率的延伸度。这种分体积在下文中可称作“微观分体积”。因此，在上下文中，利用所述数据计算出的各向同性弥散程度和/或各向异性弥散程度可以称作样品的亚分辨率或“微观”分体积的各向同性弥散程度和/或各向异性弥散程度。

每个分体积的弥散可具有弥散张量表示D。换句话说，每个分体积可具有能通过相应的弥散张量D定义的弥散。因此，在所述部分中，可以用弥散张量D的分布(例如高斯分布)表示。

优选地，可以在样品上进行多个弥散加权磁共振测量。至少两个(优选地，多个)弥散编码磁梯度脉冲序列G_i＝1…m拥有具有一至三个非零本征值的张量表示b_i。至少两个(优选地，多个)弥散编码磁梯度脉冲序列彼此不同。

根据一个实施例，所述引起各向异性弥散加权的磁梯度脉冲序列组形成第一组磁梯度脉冲序列，并且所述数据子组形成表示由第一组磁梯度脉冲序列获得的第一回波衰减曲线的第一数据子组，其中所述数据还至少包括表示由引起各向同性或各向异性弥散加权的第二组磁梯度脉冲序列获得的第二回波衰减曲线的第二数据子组。因而，可以基于表示两个不同回波衰减曲线的回波信号对各向同性和/或各向异性弥散进行定量。应当指出的是，决定“第一”和“第二”仅应理解为相应磁梯度脉冲序列组和数据子组的标记。它们并没有必然地意指任何特定顺序，例如，第一组脉冲序列在第二组脉冲序列之前施加到样品上。实际上，它们可以按相反顺序或者甚至是任意交错的方式施加。

第一组中的每个脉冲序列可被生成为使得用于所述脉冲序列的弥散编码张量的第一本征值和第二本征值彼此相等。用于第一组的弥散编码张量的第三本征值与第一本征值和第二本征值不同(从而引起各向异性弥散加权)。此外，第二组中的每个脉冲序列可为使得用于所述脉冲序列的弥散编码张量的第一本征值和第二本征值彼此相等。采用弥散编码张量的不同大小的第一本征值和第二本征值转换为使样品材料中的弥散各向异性对所得回波信号的影响发生改变。因而，以这种方式生成磁梯度脉冲序列使关于各向同性和各向异性检测样品的弥散特性称为可能。

第一组磁梯度脉冲中的脉冲序列和第二组磁梯度脉冲中的脉冲序列可以具有不同的最大梯度幅值。关于弥散加权张量，其可以表示为用于第一(或第二)组中的脉冲序列的弥散编码张量的迹线，其在整个第一(或第二)组中有所变化。因此，弥散加权的强度可以是变化的。

根据一个实施例，对于第一组中的每个脉冲序列，存在能由下式定义的第一弥散编码张量不变量Δ_b,1：

${\mathrm{\Δ}}_{b,1}=\frac{1}{b}(b_{zz}^{PAS}-\frac{b_{yy}^{PAS}+b_{xx}^{PAS}}{2}),b=b_{xx}^{PAS}+b_{yy}^{PAS}+b_{zz}^{PAS}$

其中，b_xx ^PAS表示用于所述脉冲序列的弥散编码张量的第一本征值，b_yy ^PAS表示用于所述脉冲序列的弥散编码张量的第二本征值，而b_zz ^PAS表示用于所述脉冲序列的弥散编码张量的第三本征值，并且

其中第一组脉冲序列被生成为使得第一组中的脉冲序列的第一弥散编码张量不变量Δ_b,1彼此相等。通过以这种方式控制第一组弥散编码磁梯度脉冲序列的生成，由第一数据子组表示的第一回波衰减曲线可以表示使用具有相同程度的各向异性Δ_b,1的弥散编码张量获得的回波衰减曲线。

类似地，对于第二组中的每个脉冲序列，可存在能由下式定义的第二弥散编码张量不变量Δ_b,2：

${\mathrm{\Δ}}_{b,2}=\frac{1}{b}(b_{zz}^{PAS}-\frac{b_{yy}^{PAS}+b_{xx}^{PAS}}{2}),b=b_{xx}^{PAS}+b_{yy}^{PAS}+b_{zz}^{PAS}$

其中，b_xx ^PAS表示用于所述脉冲序列的弥散编码张量的第一本征值，b_yy ^PAS表示用于所述脉冲序列的弥散编码张量的第二本征值，而b_zz ^PAS表示用于所述脉冲序列的弥散编码张量的第三本征值，并且

其中第二组脉冲序列使得第二组中的脉冲序列的第二弥散编码张量不变量Δ_b,2彼此相等，并且Δ_b,2不同于Δ_b,1。

根据一个实施例，计算各向同性弥散程度和/或各向异性弥散程度包括：

通过分析利用第一组中的脉冲序列获得的第一回波信号与利用第二组中的脉冲序列获得的第二回波信号之间的变化、改变或差异，计算各向同性弥散程度和/或各向异性弥散程度。由于第一组中的脉冲序列和第二组中的脉冲序列可被生成为呈现不同程度的各向异性(即Δ_b,1和Δ_b,2)，因此能通过第一回波信号与第二回波信号之间的变化、改变或差异(例如，与幅值相关)估算各向同性弥散程度和/或各向异性弥散程度。为简化估算，第一回波信号和第二回波信号可以利用最大梯度幅值相等的梯度脉冲(换句话说，相等的弥散加权幅值b)获得。

在本发明方法的范围内，除了所述数据的第一子组和所述数据的第二子组凹，所述数据可以至少包括第三数据子组，其表示利用引起各向异性弥散加权的第三组磁梯度脉冲序列获得的第三回波衰减曲线，

其中用于第三组中的每个梯度脉冲序列的弥散编码张量具有3个非零本征值，其中第一本征值和第二本征值彼此相等且与第三本征值不同，并且

其中，对于第三组中的每个脉冲序列，存在能由下式定义的第三弥散编码张量不变量Δ_b,3：

${\mathrm{\Δ}}_{b,3}=\frac{1}{b}(b_{zz}^{PAS}-\frac{b_{yy}^{PAS}+b_{xx}^{PAS}}{2}),b=b_{xx}^{PAS}+b_{yy}^{PAS}+b_{zz}^{PAS}$

其中，b_xx ^PAS表示用于所述脉冲序列的弥散编码张量的第一本征值，b_yy ^PAS表示用于所述脉冲序列的弥散编码张量的第二本征值，而b_zz ^PAS表示用于所述脉冲序列的弥散编码张量的第三本征值，并且

其中第三组脉冲序列使得第三组中的脉冲序列的第三弥散编码张量不变量Δ_b,3彼此相等，并且Δ_b,3不同于Δ_b,2和Δ_b,1。沿其他恒定的编码张量各向异性线(例如Δ_b,3)获得表示其他回波衰减曲线的数据能够对样品的弥散特性进行延伸检测。

根据一个实施例，第一组中的每个脉冲序列使得Δ_b,1>0，第二组中的每个脉冲序列使得Δ_b,2＝0，从而引起各向同性弥散加权，而第三组中的每个脉冲序列使得Δ_b,3<0。这使得能够对样品的弥散特性的“形状”进行估算和分析。具体地讲，使得估算弥散主要是各向同性的(即球形的)、主要是单向的(即扁球形的)还是主要是平面的(即扁长的)变成可能。

根据一个实施例，该方法还包括基于表示所述回波信号的数据(根据上文其至少包括第一数据子组、第二数据子组和第三数据子组)，计算概率分布，该概率分布指示所述回波信号中的每一个与模型各向同性弥散参数D_iso和/或模型各向异性弥散参数Δ_D的多个不同值中每一个相关的概率。除了别的以外，这使得能够分析在测量样品过程中得到的测量结果，包括呈现不同程度的各向同性弥散和/或各向异性弥散的域。通过概率分布可以确定这类域的数量(即分量)。

所述概率分布可以通过确定方程组的数值解而计算得到，该方程组将通过所述数据表示的回波信号与核函数和所述概率分布的乘积关联起来。所述方程组尤其可以是线性方程组。

概率分布可以是联合概率分布p，而核函数可以是包括至少M x N个元素的矩阵K，对于弥散加权幅值b、弥散编码张量不变量Δ_b、模型各向同性弥散参数D_iso和模型各向异性弥散参数Δ_D的组合值，每个所述元素基于下述公式的积分：

其中A＝3bD_isoΔ_bΔ_D，

对于弥散加权幅值b、弥散编码张量不变量Δ_b、模型各向同性弥散参数D_iso和模型各向异性弥散参数Δ_D的不同组合值，可以计算得到矩阵元素。

根据一个实施例，该方法还包括：

将所述第一组磁梯度脉冲序列中的每个脉冲序列分多次施加到样品上，相对于固定的实验室系采用不同的梯度脉冲取向，并通过对针对不同取向得到的回波信号测量值求平均以形成所述第一数据子组。这可以称作“粉末平均(powder averaging)”，借此，当存在域取向的某些优先排列时，能够模拟随机域取向的影响。这种“粉末平均化”也可以在第二组磁梯度脉冲序列上进行。也就是说，通过将所述第二组磁梯度脉冲序列中的每个脉冲序列分多次施加到样品上，相对于固定的实验室系采用不同的梯度脉冲取向，并通过对针对不同取向得到的回波信号测量值求平均以形成所述第二数据子组。

根据一个实施例，每一个所述弥散编码磁梯度脉冲序列G_i形成(单独的)三重受激回波序列的一部分。当对包括具有相对较短的横向弛豫时间T₂的材料域的样品进行测量时，这是尤其有利的。

根据一个实施例，该方法还包括：

基于将回波信号E与弥散编码张量b和弥散张量D关联起来的函数展开形成方程组；

通过利用由所述数据表示的回波信号测量值和弥散编码张量b_i的至少一个子组的表示值确定方程组的解，计算平均弥散张量<D>和弥散张量协方差张量S；

通过将S投影到体积基E _体积上，计算协方差张量S的不变体积分量S_体积；

通过将S投影到剪切基E _剪切上，计算协方差张量S的不变剪切分量S_剪切；以及

利用不变体积分量S_体积和/或不变剪切分量S_剪切计算各向同性弥散程度和/或各向异性弥散程度。

该实施例使得能够对微观结构弥散特性(例如，就微观弥散各向异性方面)的差异进行定量，而不需要轴向对称的弥散编码张量。特别是，不变体积分量S_体积可以生成对在上述样品部分的不同分体积之间的各向同性弥散程度的变化的估算。不变剪切分量S_剪切可以生成对在上述样品部分的不同分体积之间的各向异性弥散方向的变化的估算。

计算得到的平均弥散张量<D>可以针对上述包括多个分体积的部分进行计算。<D>可以表示对所述部分的平均弥散张量的估算值。相似地，弥散张量协方差张量S可以表示对所述部分弥散张量分布的协方差估算值。

方程组可以是线性方程组，其中由所述数据表示的回波信号测量值可以用于构成线性方程组的常数，并且所述弥散编码张量的至少一个子组可以用于构成该线性方程组的参数。特别是，该方程组等同于函数E(b)＝<exp(–<b,D>)>的累积展开。

基于不变剪切分量S_剪切和平均弥散张量<D>的平方到剪切基E _剪切上的投影的总和，可以计算各向异性弥散程度。特别是，可以基于所述总和的方差计算各向异性弥散程度。

计算得到的各向异性弥散程度可以进一步基于平均弥散张量<D>的平方到体积基E _体积上的投影与所述总和之间的比值进行计算。特别是，各向异性弥散程度可基于所述比值计算为微观部分各向异性μFA的估算值。

通过计算协方差张量S的矩阵表示与体积基E _体积的矩阵表示之间的内积，可以计算S到体积基E _体积上的投影。

通过计算协方差张量S的矩阵表示与剪切基E _剪切的矩阵表示之间的内积，可以计算S到剪切基E _剪切上的投影。

通过计算<D>平方的矩阵表示与体积基E _体积的矩阵表示之间的内积，可以计算<D>平方到体积基E _体积上的投影。

特别是，微观部分各向异性μFA可以按下述公式计算：

其中

本发明方法和上文公开的本发明方法的实施例可为用于弥散MRI的方法，其中计算的各向同性弥散程度和/或各向异性弥散程度可以用作弥散MRI数据的体素的对比参数。

附图说明

通过下面结合附图对本发明构思的优选实施例的说明性和非限制性详细描述，将能够更好地理解本发明的上述及其他目的、特征和优点。在附图中，除非另有说明，否则类似的附图标记用于指示类似的元件。

图1示出梯度波形的一些代表性示例。

图2示出理论回波衰减信号S(b,Δ_b)与bD_iso的关系。

图3和图4示出多个实验结果。

具体实施方式

为便于理解本发明构思，下面将参考附图提供一些理论概念的讨论。

各向异性高斯弥散

高斯弥散过程的方向性记录在弥散张量D中[11]。在其主轴系统(PAS)中，张量是沿对角线的元素D_xx ^PAS、D_yy ^PAS和D_zz ^PAS。在本公开的上下文中，可以方便地利用各向同性值D_iso、各向异性Δ_D和不对称性η_D表征弥散张量：

$\begin{array}{c}{D}_{iso}=(D_{xx}^{PAS}+D_{yy}^{PAS}+D_{zz}^{PAS})/3\\ {\mathrm{\&Delta;}}_D=\frac{1}{3D_{iso}}(D_{zz}^{PAS}-\frac{D_{yy}^{PAS}-D_{xx}^{PAS}}{2})\\ {\mathrm{\&eta;}}_D=\frac{D_{yy}^{PAS}-D_{xx}^{PAS}}{2D_{iso}{\mathrm{\&Delta;}}_D}\end{array}---\left(1\right)$

应当指出的是，该形式让人想起在固态NMR中用于化学位移张量的公式(例如在[12]和[13]中那样)。

各元素根据惯例|D_zz ^PAS-D_iso|>|D_yy ^PAS-D_iso|>|D_xx ^PAS-D_iso|进行排序。对方程(1)中的数值因子进行选择以得到范围在-1/2≤Δ_D≤1和0≤η_D<1的参数。当η_D＝0时，张量是轴对称的。Δ_D的正值和负值分别对应于扁长形(prolate)和扁球形(oblate)的张量形状。用于对元素排序的惯例保证了Z轴是长球形和扁球形的张量的主对称轴。但是，应当注意，在不脱离本发明构思的范围的前提下，也可以使用其他惯例(例如，使用x轴或y轴作为对称轴)。

利用欧拉旋转矩阵R_z(α)、R_y(β)和R_z(γ)，PAS从实验室系的一般旋转得到下述的实验室系弥散张量zz-元素：

$D_{zz}(\mathrm{\&alpha;},\mathrm{\&beta;})=D_{zz}^{PAS}{\cos }^2\mathrm{\&beta;}+\frac{D_{yy}^{PAS}+D_{xx}^{PAS}}{2}{\sin }^2\mathrm{\&beta;}-\frac{D_{yy}^{PAS}-D_{xx}^{PAS}}{2}{\sin }^2\mathrm{\&beta;}cos2\mathrm{\&alpha;}---\left(2\right)$

利用方程(1)的关系式，该表达式可以重排为：

D_zz(α，β)＝D_iso[1+Δ_D(2P₂(cosβ)-η_Dsin²βcos2α)] (3)

其中，P₂(x)＝(3x²-1)/2是第二勒让德多项式(Legendre polynomial)。对于轴对称张量，η_D＝0，方程(3)简化为：

D_zz(β)＝D_iso[1+2Δ_DP₂(cosβ)] (4)

。在一致性检验中，注意到正如所预期的，将方程(1)代入方程(4)得到D_zz(0)＝D_zz ^PAS和D_zz(π/2)＝(D_xx ^PAS+D_yy ^PAS)/2。

弥散加权张量b

利用时变磁场梯度G^T(t)＝{G_x(t),G_y(t),G_z(t)}，用关于平移运动的信息编码NMR信号。瞬时移相向量q(t)通过下述时间积分得到：

$q\left(t\right)=\mathrm{\&gamma;}\underset{0}{\overset{t}{\&Integral;}}G\left(t^{\&prime;}\right){\mathrm{dt}}^{\&prime;}---\left(5\right)$

其中，γ是研究的核的磁旋比。在测量过程中，当自旋磁化是重相时，即q(τ)＝0时，可以在回波时间τ记录回波信号，，换句话说，在时间点t＝τ，此刻自旋磁化是重相的。假定为高斯弥散，信号振幅S可以写为：

S＝S₀exp(-b∶D) (6)

其中，S₀是在零梯度振幅的信号强度(即，非弥散加权回波信号)，b:D表示广义数积，定义为：

$b:D=\underset{i}{\&Sigma;}\underset{j}{\&Sigma;}{b}_{ij}{D}_{ij}---\left(7\right)$

弥散加权矩阵b通过以下公式得到：

$b=\underset{0}{\overset{\mathrm{\&tau;}}{\&Integral;}}q\left(t\right)q^T\left(t\right)dt---\left(8\right)$

与方程(1)类似，b-矩阵可以用总弥散加权b、各向异性Δ_b和不对称性η_b表征：

$\begin{array}{c}b=b_{xx}^{PAS}+b_{yy}^{PAS}+b_{zz}^{PAS}\\ {\mathrm{\&Delta;}}_b=\frac{1}{b}(b_{zz}^{PAS}-\frac{b_{yy}^{PAS}+b_{xx}^{PAS}}{2})\\ {\mathrm{\&eta;}}_b=\frac{3}{2}\frac{b_{yy}^{PAS}-b_{xx}^{PAS}}{{\mathrm{b\&Delta;}}_b}\end{array}---\left(9\right)$

可注意到，通过方程(5)到方程(8)中的b-矩阵的定义可以在关于弥散NMR和MRI的标准教科书中找到，参见例如Price的第4.4.1章[14]或Callaghan的第9.7.2章[15]。但是，使用固态NMR术语的b-矩阵表征是新的，并且如下文所示，不但简化了标记，还提供了用于设计测量协议和分析方法的框架。b-矩阵的元素根据二阶张量规则[16]旋转变换。因此，在下文中，b-矩阵(即，b)可因而被称作b-张量。

q-向量的可变角旋转

在球面坐标中，q-向量可以用其倾斜角ζ(t)、方位角Ψ(t)和幅值qF(t)定义，其中q是最大幅值，F(t)是归一化到区间0≤F(t)≤1的时间依赖性幅值。可以从下述关系式中得到笛卡尔分量：

$\begin{array}{c}{q}^T\left(t\right)=\{q_x\left(t\right),q_y\left(t\right),q_z\left(t\right)\}\\ =qF\left(t\right)\{\sin \mathrm{\&zeta;}\left(t\right)\cos \mathrm{\&psi;}\left(t\right),\sin \mathrm{\&zeta;}\left(t\right)\sin \mathrm{\&psi;}\left(t\right),\cos \mathrm{\&zeta;}\left(t\right)\}\end{array}---\left(10\right)$

将方程(10)代入方程(8)，利用标准三角关系得到下述用于b-张量元素的表达式：

$\begin{array}{c}{b}_{xx}=q^2\underset{0}{\overset{t_E}{\&Integral;}}{F}^2\left(t\right)\frac{1-{\cos }^2\mathrm{\&zeta;}\left(t\right)+{\sin }^2\mathrm{\&zeta;}\left(t\right)\cos \&lsqb;2\mathrm{\&psi;}\left(t\right)\&rsqb;}{2}dt\\ b_{yy}=q^2\underset{0}{\overset{t_E}{\&Integral;}}{F}^2\left(t\right)\frac{1-{\cos }^2\mathrm{\&zeta;}\left(t\right)-{\sin }^2\mathrm{\&zeta;}\left(t\right)\cos \&lsqb;2\mathrm{\&psi;}\left(t\right)\&rsqb;}{2}dt\\ b_{zz}=q^2\underset{0}{\overset{t_E}{\&Integral;}}{F}^2\left(t\right){\cos }^2\mathrm{\&zeta;}\left(t\right)dt\\ b_{xy}=b_{yx}=q^2\underset{0}{\overset{t_E}{\&Integral;}}{F}^2\left(t\right)\frac{{\sin }^2\mathrm{\&zeta;}\left(t\right)\sin \&lsqb;2\mathrm{\&psi;}\left(t\right)\&rsqb;}{2}dt\\ b_{xz}=b_{zx}=q^2\underset{0}{\overset{t_E}{\&Integral;}}{F}^2\left(t\right)\frac{\sin \&lsqb;2\mathrm{\&zeta;}\left(t\right)\&rsqb;\cos \mathrm{\&psi;}\left(t\right)}{2}dt\\ b_{yz}=b_{zy}=q^2\underset{0}{\overset{t_E}{\&Integral;}}{F}^2\left(t\right)\frac{\sin \&lsqb;2\mathrm{\&zeta;}\left(t\right)\&rsqb;\sin \mathrm{\&psi;}\left(t\right)}{2}dt\end{array}---\left(11\right)$

如Eriksson等人的文献[4]中所示，若ζ是常量，含Ψ(t)的所有术语为零，q-向量的轨迹具有至少三重对称，即

$\begin{array}{c}\mathrm{\&psi;}\left(t\right)=\mathrm{\&psi;}\left(0\right)+\mathrm{\&psi;}(t+\frac{n\mathrm{\&tau;}}{N})+\frac{2\mathrm{\&pi;}n}{N}\\ F\left(t\right)=F(t+\frac{n\mathrm{\&tau;}}{N})\\ N\&GreaterEqual;3,n=1,2,...,N\end{array}---\left(12\right)$

另一种将方程(11)中带Ψ(t)的术语置零的方法是当所述轨迹满足如下条件(参见参考文献[4])：

$\mathrm{\&psi;}\left(t\right)=\mathrm{\&psi;}\left(0\right)+\frac{2\mathrm{\&pi;}n}{t_d}\underset{0}{\overset{t}{\&Integral;}}F{\left(t^{\&prime;}\right)}^2{\mathrm{dt}}^{\&prime;}---\left(13\right)$

其中，n是非零整数，t_d是由下式得到的有效弥散时间：

$t_d=\underset{0}{\overset{\mathrm{\&tau;}}{\&Integral;}}F{\left(t\right)}^{2}dt---\left(14\right)$

从几何学来说，Ψ(t)和F(t)的这种调制对应于q-向量绕z轴以与F(t)²成正比的角速度dΨ(t)/dt自旋，同时沿循孔径2ζ的圆锥表面上的路径。

对方程(11)中的b-张量元素显式求值后得到：

$\begin{array}{c}{b}_{xx}=b_{yy}=q^2{t}_d\frac{1-{\cos }^2\mathrm{\&zeta;}}{2}\\ b_{zz}=q^2{t}_d{\cos }^2\mathrm{\&zeta;}\\ b_{xy}=b_{yx}=b_{xz}=b_{zx}=b_{xz}=b_{zx}=0\end{array}---\left(15\right)$

将其代入方程(9)生成：

$\begin{array}{c}b=q^2{t}_d\\ {\mathrm{\&Delta;}}_b=P_2\left(\cos \mathrm{\&zeta;}\right)\\ {\mathrm{\&eta;}}_b=0\end{array}---\left(16\right)$

假定q-向量的轨迹遵守方程(12)或(13)，b-张量是轴对称的，其中z轴为主对称轴，并具有通过经由在魔角ζ＝acos(1/3^1/2)(即，使用来自参考文献[4]的术语)的各向同性弥散加权将角ζ从常规的在ζ＝0时单向弥散加权改变为在ζ＝π/2时所谓的循环编码(即，使用来自参考文献[17]的术语)的可调整的各向异性。只要Ψ(t)和F(t)的时间调制保持相同，ζ的调整仅会影响Δ_b，而不会影响q、t_d或b的值。

梯度调制函数的数值优化

对于q-向量的给定调制F(t)，其时间依赖性取向通过角Ψ(t)(通过方程(13)的积分得到)和选定的ζ常数值得到。q-向量的笛卡尔分量可以用方程(10)计算得到，而梯度调制函数由以下导数得到：

$G\left(t\right)=\frac{1}{\mathrm{\&gamma;}}\frac{dq\left(t\right)}{dt}---\left(17\right)$

概念性简单调制函数是，在区间0≤t≤τ中F(t)＝1，否则F(t)＝0，对应于t＝0和τ时无限短且强的梯度脉冲。对于在具有有限梯度能力的MRI硬件上的实际执行，需要在F(t)＝0和1之间具有更少的急剧转变。[10]中描述了用于在临床MRI扫描仪上寻找用于q-MAS弥散加权的最佳梯度波形的一种可能的程序。简单来说，F(t)可以展开为：

$F\left(t\right)=\underset{m=1}{\overset{M}{\&Sigma;}}\frac{a_m}{2}\&lsqb;1-cos\left(\frac{2\mathrm{\&pi;}mt}{\mathrm{\&tau;}}\right)\&rsqb;---\left(18\right)$

其中，对系数a_m进行迭代优化以得到给定波形持续时间τ约束范围内的最大弥散加权和三个梯度通道中的每一个上的最大梯度幅度G_max。优化途径的详细内容在参考文献[10]中有所描述。表1中列出了对于圆锥角ζ＝0和ζ＝π/2联合优化的最终结果，得到一个轴向梯度调制函数和两个径向梯度调制函数，它们可以叠加以得到实验室系中的任意圆锥角和取向的q-向量调制。

表1中的系数a_m可以称作经数值优化的q-VAS梯度波形的系数。利用表1，通过选择回波时间值τ以及q-向量倾斜角ζ和幅值q可以得到明确的梯度调制函数，随后利用方程(18)计算归一化的q-幅值调制F(t)。利用方程(14)可以计算有效弥散时间t_d。利用方程(13)可以计算q-向量方位角Ψ(t)。利用方程(10)可以计算q-向量的笛卡尔分量。最后，利用方程(17)可以计算梯度向量的笛卡尔分量。

图1中示出了包括相应的q-向量轨迹和b-张量元素的梯度波形的一些代表性示例。由于q-向量的方位角Ψ(t)随时间变化，图1中例示的调制可以称作q-向量的可变角旋转(q-VAS)。如图1第1列示出，图1的1-4行对应于角ζ＝0°、35.3°、54.7°和90°。图1的第2列示出了从表1列出的系数得到的梯度调制函数G_x(t)、G_y(t)和G_z(t)(分别为点线、短划线和实线)。图1的第3列示出了相应的q-向量调制函数q_x(t)(点线)、q_y(t)(短划线)、q_z(t)(实线)和q(t)(点划线)。图1的第4列示出了相对于正x-、y-、z-轴的q-向量轨迹的3D曲线图(黑色线)。可以看出，第二和第三行中的q-向量轨迹位于具有孔径2ζ的圆锥表面上。第一行中的q-向量轨迹与z-轴对齐。第四行中的q-向量轨迹位于xy-平面上。图1的第5列示出了相应的b-张量各向异性Δ_b＝1、0.5、0和-0.5。

表1

有效弥散系数D _zz ^eff

将方程(15)代入方程(7)得到：

$b:D=q^2{t}_d\&lsqb;(D_{xx}+D_{yy})\frac{1-{\cos }^2\mathrm{\&zeta;}}{2}+D_{zz}{\cos }^2\mathrm{\&zeta;}\&rsqb;---\left(19\right)$

，使用方程(16)中的关系式以及弥散张量的迹线的旋转不变量D_iso＝(D_xx+D_yy+D_zz)/3，其能够改写为下式：

b∶D＝b[D_iso+Δ_b(D_zz-D_iso)] (20)

。假定弥散张量是轴对称的，代入方程(4)得到：

b∶D＝bD_iso[1+2Δ_bΔ_DP₂(cosβ)] (21)

根据以下公式，b-值后的因子可以解释为有效弥散系数D_zz ^eff(β)，其取决于在实验室系中根据下式的角β的b-张量的各向异性和弥散张量的取向：

$D_{zz}^{eff}\left(\mathrm{\&beta;}\right)=D_{iso}\&lsqb;1+2{\mathrm{\&Delta;}}_b{\mathrm{\&Delta;}}_D{P}_2\left(cos\mathrm{\&beta;}\right)\&rsqb;---\left(22\right)$

比较方程(4)和(22)表明，内在弥散各向异性的影响按照Δ_b值缩放，而该值取决于z-轴和旋转q-向量之间的角ζ。D_zz ^eff的值在D_zz ^eff(0)＝D_iso(1+2Δ_bΔ_D)和D_zz ^eff(π/2)＝D_iso(1-Δ_bΔ_D)的范围内。

粉末平均信号(powder-averaged signal)衰减和有效弥散分布

考虑了宏观样品，其由随机取向的、具有相同D_iso和Δ_D值的微观各向异性域的集合组成。在存在域取向的某些优先排列的情况下，通过对数据进行“粉末平均”(即针对q-轨迹的对称轴的一系列方向记录数据，随后对各方向的结果进行平均)，能够模拟随机域取向的影响。在对由磁梯度脉冲序列(由以z-轴作为对称轴(即η_D＝0)的q-向量表示)编码的宏观样品的测量中，每个域产生了可通过将方程(21)代入方程(6)计算得到的信号。

$\frac{S(b,{\mathrm{\&Delta;}}_b)}{S_0}=\exp \{-{\mathrm{bD}}_{iso}\&lsqb;1+2{\mathrm{\&Delta;}}_b{\mathrm{\&Delta;}}_D{P}_2\left(cos\mathrm{\&beta;}\right)\&rsqb;\}---\left(23\right)$

对来自所有域的贡献进行积分，得到：

$\frac{S(b,{\mathrm{\&Delta;}}_b)}{S_0}=\underset{0}{\overset{\mathrm{\&pi;}/2}{\&Integral;}}{P}_{\mathrm{\&beta;}}\left(\mathrm{\&beta;}\right)\exp \{-{\mathrm{bD}}_{iso}\&lsqb;1+2{\mathrm{\&Delta;}}_b{\mathrm{\&Delta;}}_D{P}_2\left(cos\mathrm{\&beta;}\right)\&rsqb;\}d\mathrm{\&beta;}---\left(24\right)$

其中，P_β(β)是角分布函数，在区间0≤β≤π/2内归一化。域取向的随机分布对应于P_β(β)＝sinβ，基于方程(24)的积分求值，其产生：

$\begin{array}{c}\frac{S(b,{\mathrm{\&Delta;}}_b)}{S_0}=\exp (-{\mathrm{bD}}_{iso})\&CenterDot;\exp \left(\frac{A}{3}\right)\&CenterDot;\frac{\sqrt{\mathrm{\&pi;}}}{2}\frac{\mathrm{\&gamma;}(1/2,A)}{\sqrt{A}}\\ A=3{\mathrm{bD}}_{iso}{\mathrm{\&Delta;}}_b{\mathrm{\&Delta;}}_D\end{array}---\left(25\right)$

在方程(25)中，γ(s,x)是低阶不完全伽马函数。该函数可以例如方便地在Matlab中利用“gammainc”函数进行数值求值。应当注意，虽然当自变量为负时，方程(25)中γ和平方根因子是假设的，但它们的比值仍为真实的且为正。基于γ(1/2,x)∝erf(x^1/2)，γ(s,x)因子还可以错误函数“erf”写入。

图2示出了根据方程(25)的理论信号S(b,Δ_b)与bD_iso和Δ_bΔ_D的关系，其中b是弥散加权幅值，D_iso是各向同性弥散率，Δ_b和Δ_D分别是弥散加权b和弥散张量D的各向异性。所述表面利用方程(25)计算，参数b、D_iso、Δ_b和Δ_D由方程(1)和(9)中的相应张量本征值限定。

方程(25)中的表达式提供了分析实验数据的基础。当Δ_b＝0，方程(25)简化为单指数衰减：

$\frac{S(b,{\mathrm{\&Delta;}}_b=0)}{S_0}=\exp (-{\mathrm{bD}}_{iso})---\left(26\right)$

从而提供一种简单的提取弥散张量各向同性部分的方法，而不存在各向异性的混杂影响。当在1/D_iso阶在某些有限值下b保持不变，Δ_D值可以从S的特征变量根据Δ_b进行确定。

根据下述公式，方程(25)中的多指数信号衰减可以解释为有效弥散率的P(D_zz ^eff)分布的拉普拉斯变换：

$\frac{S\left(b\right)}{S_0}=\underset{0}{\overset{\mathrm{\&infin;}}{\&Integral;}}P\left(D_{zz}^{eff}\right)\exp (-{\mathrm{bD}}_{zz}^{eff}){\mathrm{dD}}_{zz}^{eff}---\left(27\right)$

其中

$P\left(D_{zz}^{eff}\right)=\frac{1}{2\sqrt{3D_{iso}{\mathrm{\&Delta;}}_b{\mathrm{\&Delta;}}_D\&lsqb;D_{zz}^{eff}-D_{iso}(1-{\mathrm{\&Delta;}}_b{\mathrm{\&Delta;}}_D)\&rsqb;}}---\left(28\right)$

D_zz ^eff在从D_iso(1-Δ_bΔ_D)到D_iso(1+2Δ_bΔ_D)的范围内；否则，P(D_zz ^eff)＝0。该分布在D_zz ^eff＝D_iso(1-Δ_bΔ_D)有一奇点，对应于域取向β＝π/2。方程(28)与针对轴对称化学位移各向异性张量获得的“粉末图”NMR光谱类似。该分布的平均值是D_iso，而第二中心矩μ2和第三中心矩μ3分别是：

${\mathrm{\&mu;}}_2=\frac{4}{5}{\left(D_{iso}{\mathrm{\&Delta;}}_b{\mathrm{\&Delta;}}_D\right)}^2---\left(29\right)$

和

${\mathrm{\&mu;}}_3=\frac{16}{35}{\left(D_{iso}{\mathrm{\&Delta;}}_b{\mathrm{\&Delta;}}_D\right)}^3---\left(30\right)$

多域材料的信号强度

由具有不同的D_iso和Δ_D值的域的集合组成的材料产生粉末平均化的信号，其可以表达为积分变换式：

$\frac{S(b,{\mathrm{\&Delta;}}_b)}{S_0}=\underset{-1/2}{\overset{1}{\&Integral;}}\underset{0}{\overset{\mathrm{\&infin;}}{\&Integral;}}K(b,{\mathrm{\&Delta;}}_b,D_{iso},{\mathrm{\&Delta;}}_D)P(D_{iso},{\mathrm{\&Delta;}}_D){\mathrm{dD}}_{iso}{\mathrm{d\&Delta;}}_D---\left(31\right)$

其中P(D_iso,Δ_D)是D_iso和Δ_D的2D联合概率分布。内核K(b,Δ_b,D_iso,Δ_D)由方程(25)的右侧得到，并将2D“解析空间”(D_iso,Δ_D)映射到2D“采集空间”(b,Δ_b)。从实验数据组I(b,Δ_b)估算P(D_iso,Δ_D)可以认为是不适定性问题，并且可以得益于特殊方法以确保数值稳定性。基于用于解决NMR弥散中的类似问题的方法和弛豫相关法[18-20]，以及一些来自压缩感知[8,9]的其他启发，可以使用下文概述的程序。

方程(31)可以是离散的并且可以写成矩阵形式：

s＝Kp (32)

其中，s是针对(b,Δ_b)的M个组合测得的信号振幅的列向量，p是用于(D_iso,Δ_D)的N个离散对的所求概率列向量，K是针对MxN格的(b,Δ_b)和(D_iso,Δ_D)对计算的内核K(b,Δ_b,D_iso,Δ_D)的矩阵版本。方程(32)是一组线性方程，如果问题是超定的(即M>N)，则原则上可以通过直接矩阵反演法对该方程求解。遗憾的是，内核的“平滑”使得这种直接方法对实验噪声极其敏感，对于输入数据向量s而言微小的变化会引起解向量p的剧烈波动。假定材料由几个离散分量组成，可能有利的是寻找稀疏解，这表示大多数元素p是零。

基于上述考虑，与实验数据s一致的稀疏解p可以通过最小化下述函数进行估算：

$f\left(p\right)=\underset{i=1}{\overset{N}{\&Sigma;}}\underset{j=1}{\overset{M}{\&Sigma;}}|K_{ij}{p}_j-s_i{|}^2+\mathrm{\&lambda;}\underset{j=1}{\overset{M}{\&Sigma;}}|p_j|---\left(33\right)$

其中，第一项是最小二乘失配值，第二项是利用正则化参数λ加权的l₁-范数。结合元素p的非负约束，方程(33)可以用公式表示为二次规划问题：

f(p)＝p^TK^TKp-2s^TKp+s^Ts+λ1^Tp (34)

其能容易地利用例如Matlab的优化工具箱中的“quadprog”函数求解。

示例性实验

在该部分及之后的部分，将描述用于证明原理的实验的多个示例及它们的结果。根据这些示例，实验在带有洗涤剂气溶胶-OT和由H₂O和D₂O的等摩尔混合物组成的水的溶致液晶上进行。基于平衡相图(参见例如参考文献[21])，选择洗涤剂浓度以提供三种不同的液晶相：层状相(25重量％和75重量％)、双连续立方相(80重量％)和反六角相(85重量％)。首先称量样品加入10mL瓶中，充分混合，随后取400μL转移至5mm一次性NMR管中。通过记录小角度x-射线散射图和²H NMR光谱，单独检验相对称性。通过将装有25重量％气溶胶-OT的5mm NMR管插入装有癸醇的10mm NMR管中，制备具有两种不同弥散张量组分的样品。

在具有11.7T磁铁且配备了Bruker MIC-5微成像探针并且能在三个正交方向上提供最大3T/m的磁场梯度的Bruker Avance-II 500MHz光谱仪上进行弥散磁共振实验。q-VAS梯度调制通过在标准¹H回波脉冲序列中的180°RF脉冲的两侧上包括图1所示的波形来实现。在后半个自旋回波中记录的信号经傅里叶变换后得到高分辨光谱，从而能够将水¹H信号与源自洗涤剂的信号分开。为了研究具有短的横向弛豫时间T₂的系统，有利的是通过调整三重受激回波序列中的脉冲场梯度的方向来实现弥散加权各向异性Δ_b的变化。三个梯度方向具有方位角Ψ＝0°、120°和240°，以及倾斜角ζ，其通过方程(16)得到Δ_b。

利用三重受激回波研究75重量％、80重量％和85重量％气溶胶-OT/水样品，采用三对持续时间δ＝1ms以及前沿分离点Δ＝100ms的梯度脉冲实施。通过改变梯度脉冲的振幅和角ζ，对采集空间(b,Δ_b)矩形网格进行采样。最大梯度振幅接近1T/m，并且针对不同样品进行调整以得到大致相同的最大信号衰减。立方相样品和反六角相样品均得到可用三重受激序列检测到的足够窄的气溶胶-OT共振线。根据经验，发现水和气溶胶-OT共振线在25℃重叠，但是通过将温度升高至80℃，该重叠变得不显著。为保持一致，75重量％、80重量％和85重量％的样品都在80℃下进行研究。根据平衡相图，样品保持与在25℃时相同的液晶相，并且通过²H光谱测量得到确认。

利用在梯度调制持续时间τ＝140ms以及最大梯度振幅为0.090T/m的自旋回波版本，在25℃下研究了25重量％气溶胶-OT/水/癸醇样品。通过改变最大梯度振幅和ζ以Z字形图形对(b,Δ_b)空间进行采样。虽然能够通过傅里叶变换分开共振线，仍将它们共同记录下来以得到包含具有不同弥散行为的多分量的信号。实际上，只有水和癸醇的端甲基具有足够长的T₂，以在漫长的自旋回波序列下留存。

为了确保数据与随机分布的域取向相对应，如方程(25)要求的那样，重复进行采集并对31个不同的“圆锥取向”(即，q-向量轨迹的主对称轴的取向)求平均。这些方向根据静电排斥方案(参见参考文献[22]和[23])进行选择。

示例性实验的结果

图3示出了层状、双连续立方和反六角液晶相的实验数据。数据用方程(25)拟合，得到与已知微观结构一致的弥散各向异性Δ_D值。值得注意的是，Δ_D的记号能够根据Δ_b从信号特征变量提取，只要bD_iso是一阶或更高阶。

更详细地说，图3示出了针对(a)层状，(b)双连续立方和(c)反六角类型的AOT/水液晶的表示测量的信号衰减S(b,Δ_b)与弥散加权幅值b和各向异性Δ_b关系的数据。最上面一行示出了这些类型的示意性图示。这些几何形状表征数十纳米长度范围中相应的水区室几何结构。实心圆表示在(b,Δ_b)空间中的矩形网格上采样的实验数据点。网格图示出了利用初始信号强度S₀、各向同性弥散率D_iso和弥散各向异性Δ_D作为可调参数时方程(25)对实验数据的拟合。该拟合得到D_iso/10^-9m²s^-1＝3.53(层状)，2.37(立方)和1.22(反六角)，以及Δ_D＝-0.38(层状)，0.00(立方)和0.80(反六角)。

图4示出了液晶/癸醇样品的结果。(b,Δ_b)空间的Z字形采样便于在信号与b关系的2D图中显示数据。在这样的图中，具有不同的D_iso的多组分的存在可以视作下包络数据的曲率，而Δ_D的非零值导致在由采样图给定的频率处发生振荡。振荡的振幅与Δ_D的幅值有关，而与Δ_b＝1和-0.5时对应的局部极值之间的比值给出Δ_D的记号。

通过前述的l₁-正则化模型无关法将数据S(b,Δ_b)转换为概率分布P(D_iso,Δ_D)。得到的分布包含在(D_iso＝10^-10m²/s,Δ_D＝0)和(D_iso＝10^-9m²/s,Δ_D＝-0.5)的组分，分别与癸醇和水对应。一旦知道分布包含两个组分，通过两组分模型拟合能更精确地确定它们的坐标和幅值，其中来自每个组分的信号利用方程(25)描述。这种拟合的结果也在图4b中示出。所得的结果与已知的癸醇各向同性弥散和液晶相的层状对称一致。更详细地说，图4a示出了针对装有AOT/水层状液晶(内管)和癸醇(外管)的管套管样品的实验水及癸醇信号S(b,Δ_b)与弥散加权幅值b的关系。最上面示出了(b,Δ_b)-空间的Z字形采样图。图4b示出了与S(b,Δ_b)数据一致的概率密度P(D_iso,Δ_D)。轮廓线示出了l₁-正则化模型无关估计的结果，而十字线表示两组分模型拟合结果，给出D_iso/10^-9m²s^-1＝0.083(癸醇)和1.33(水)以及Δ_D＝0.00(癸醇)和-0.496(水)。图4a中的点表示实验数据，而黑线代表两组分模型拟合。

实施例描述

根据本发明构思，提供了一种用于定量样品中的各向同性弥散和/或各向异性弥散的方法。参考前述说明，各向同性弥散可例如通过弥散张量D的各向同性值D_iso进行定量(如上面结合方程1所定义)。各向异性弥散可以例如通过弥散张量D的各向异性Δ_D进行定量。

在所述方法中进行的各种计算可以例如利用可存储于非暂态计算机存储介质上或者内嵌到其上的一组软件指令实现。

所述方法可利用当前技术中的NMR光谱仪或MRI装置进行。如本领域中众所周知，这种装置可包括一个或多个用于控制装置运行的处理器，特别是生成磁梯度脉冲序列、采集回波信号以及对测量的信号进行采样和数字化，用于形成表示所得回波信号的数据。生成弥散编码磁梯度脉冲序列可利用软件指令进行，所述软件指令可存储于计算机可读介质上(例如非暂态计算机可读存储介质)并通过所述装置的一个或多个处理器执行。所述软件指令可以例如存储于装置的存储器的程序/控制段，所述装置的一个或多个处理器可访问所述程序/控制段。所收集的表示测量值的数据可存储于所述装置的、或计算机的或者可与所述装置相连接的类似装置的数据存储器中。所述方法的计算可利用软件指令实现，所述软件指令可存储于计算机可读介质上并通过所述装置的一个或多个处理器执行。但是，同样可能的是，在与NMR光谱仪或MRI装置分开的装置上进行计算，例如在计算机上。所述装置和所述计算机可以例如设置为经由通信网络(例如LAN/WLAN)或者经由一些其他串行或并行通信接口进行通信。应当进一步指出的是，不利用软件指令，所述方法的操作可在所述装置的专用电路中实施，例如在一个或多个集成电路中、在一个或多个专用集成电路(ASIC)或者现场可编程门阵列(FPGA)中实施，仅举几例。

所述方法包括利用可表示为G_i＝1…m的弥散编码磁梯度脉冲序列对样本进行弥散加权磁共振测量。参照前述说明，每个磁梯度脉冲序列G_i可由弥散编码张量b_i定义，或者具有弥散编码张量b_i形式的表示，弥散编码张量b_i具有三个非零本征值，其中(方程(5))，且q_i(t)为与(方程(8))成正比的移相向量。

为了得到回波信号，每个弥散编码磁梯度脉冲序列可利用一个或成像的磁梯度和任选地磁梯度校正梯度补充，正如本领域中众所周知的。因此，每个弥散加权磁共振测量可利用包括弥散编码磁梯度脉冲序列Gi的磁脉冲梯度序列、成像磁梯度序列和任选地，校正磁梯度序列进行。在一些情况下，这些序列可适时重叠。但是，即使在这种情况下，至少一部分序列可由具有上述特性的弥散编码张量b_i描述或表示。

所述方法还包括收集表示回波信号测量值的数据。数据收集可包括对从样本的目标部分接收的回波信号进行采样和数字化。与上述一致，每个回波衰减测量可为回波信号振幅S_i的形式，假定为高斯弥散，该回波信号振幅S_i具有如方程(6)所给出的与弥散编码张量b_i和D的依赖关系。测量可为自旋回波测量或受激回波测量，例如三重受激回波测量。

所述部分可如上所述包括多个“微观”分体积，所述分体积呈现不同的各向同性弥散程度(例如上文定义的Diso)或不同的各向异性弥散程度和/或取向(例如上文定义的Δ_D)。如本领域众所周知，NMR光谱仪或MRI装置的空间分辨率尤其受限于磁场强度、施加至样品的梯度脉冲序列幅值和转换速率。以下公开的数据分析使得能够估算所述部分内的微观分体积的弥散特性，即，超出了测量的常规分辨率的限制。为了辨别对应于所述部分的回波信号分量，可对来自样本的测量信号进行快速傅里叶变换，如在本领域内众所周知的那样，从而将来自样本的每个回波信号的光谱分量变换为样本的多个空间区域。

根据本发明的方法，表示回波信号测量值的数据的至少一个子组利用引起各向异性弥散加权的一组磁梯度脉冲序列得到，其中用于所述磁梯度脉冲序列组中的每个梯度脉冲序列的弥散编码张量具有三个非零本征值，三个本征值中的至少一个与另外两个本征值不同。参见方程(9)，每个脉冲序列的三个非零本征值可表示为b_xx ^PAS、b_yy ^PAS和b_zz ^PAS。类似地，所述组中的每个脉冲序列能够由总弥散加权b、各向异性Δ_b和不对称性η_b限定，遵循方程(9)中的定义。这些参数形成每个脉冲序列的不变参数。参数b表示梯度序列的弥散加权幅值。

现在将介绍一个实施例，其中回波信号利用生成的磁梯度脉冲序列G_i＝1…m得到，使得不对称参数η_b等于零。这意指对应的弥散加权编码张量的至少两个本征值相等(即，b_xx ^PAS＝b_yy ^PAS)。表示回波信号测量值的数据可按照以上所述收集。所述数据可包括多个不同的数据子组，每个子组利用不同的脉冲序列组得到。例如，所述数据可包括第一数据子组、第二数据子组和第三数据子组。第一数据子组、第二数据子组和第三数据子组可分别表示利用第一组(记为G_i＝1…j)、第二组(记为G_i＝j+1…k)和第三组(记为G_i＝k+1…m)磁梯度脉冲序列得到的回波信号测量值。尽管提及了这三个数据子组和脉冲序列组，但是所述方法也适用于更多或更少的这种数据子组和脉冲序列组。

第一组脉冲序列可分别生成为具有相同程度的各向异性Δ_b,1但不同的最大振幅b(根据方程(9)和(16)定义)。类似地，第二组脉冲序列可分别生成为具有相同程度的各向异性Δ_b,2但不同的最大振幅b。同样地，第三组脉冲序列可分别生成为具有相同程度的各向异性Δ_b,3但不同的最大振幅b。例如，脉冲序列可利用结合表1描述的方法生成。参照图2，可以理解，第一数据子组、第二数据子组和第三数据子组各自表示沿相应的恒定的各向异性线(即Δ_b,1和Δ_b,2)得到的相应回波衰减曲线。利用因此收集的数据，可计算出弥散张量D的D_iso和/或弥散张量D的Δ_D。

对于包括具有大致相似程度的各向异性和各向同性弥散的随机定向微观分体积的集合的样本，可通过利用方程(25)拟合数据计算用于表征每个分体积中的弥散的随机取向弥散张量D的D_iso和Δ_D(的估算值)。例如，如果利用Δ² _b＝0得到第二数据子组，那么弥散张量D的各向同性值D_iso可利用方程26和第二数据子组直接计算得出。

如果样品的分体积中存在弥散取向的某些优先排列，则可应用粉末平均化，其中第一组磁梯度脉冲序列、第二组磁梯度脉冲序列和第三组磁梯度脉冲序列中的每一组中的每个脉冲序列都可多次施加至样品，且相对于固定的实验室系采用不同的梯度脉冲取向。然后，第一、第二和第三数据子组可通过对针对不同取向得到的回波信号测量值取平均形成。之后，可以与针对包括随机取向域的样品相同的方式应用方程(25)和(26)。

对于包括呈现不同程度的各种同性和/或各向异性弥散的微观分体积的集合的样本(即，两组分或多组分材料)，所述数据(进行或不进行粉末平均化取决于是否存在取向的优先排列)可用于估算每个组分的D_iso和/或Δ。具体地讲，可采用结合方程(32-34)描述的数据分析方法。核矩阵可利用成对的模型参数(b,Δ_b)和模型参数(D_iso,Δ_D)计算出来，其中(b,Δ_b)的值对应于在测量期间使用的每个梯度脉冲序列的值，并且，例如基于特定样本的(D_iso,Δ_D)的可能数值范围的先验知识，选择(D_iso,Δ_D)的值以覆盖样品的目标区域。如果与弥散MRI结合使用，D_iso和/或Δ_D可用于生成表示样品的所述部分的体素的对照。可进行类似的计算，以生成表示样本的其他部分的体素的对照。

尽管在上文中描述了利用(例如)方程(25-26)计算D_iso和Δ_D，但应当注意，也可采用其他方法计算或估算参数。测量数据也可利用不同的模型函数进行拟合，所述模型函数将弥散加权信号与相关的弥散度量(例如D_iso或Δ_D)关联起来。作为方程(25)的替代形式一个示例，可利用根据扩散率的分布的时刻μ₂、μ₃等的方程(25)的展开形式。另一个示例是通过伽马分布函数取扩散率的分布的近似值。

在上文中，已经公开了各种实施例，其中弥散参数基于回波信号测量值进行计算，回波信号测量值利用轴对称的弥散编码张量得到。下文将公开本发明方法的不需要轴对称弥散编码张量的实施例。一个目的是提供能够定量微观结构中的差异的方法，例如在基于一部分样品的“微观”分体积的弥散特性的微观弥散各向异性方面。为便于理解以下实施例，现在将提供一些理论概念的讨论。

理论

考虑包括分体积集合(例如，“微观”分体积)的样品的一部分(将对其进行弥散NMR/MRI)，其中在每个分体积中所述弥散为高斯型并通过弥散张量D进行描述。在所述部分内的这些微环境中的弥散特性可利用高斯分布关于张量进行建模。张量D因此可称为具有期望值<D>的随机变量,其中<·>表示关于该部分中的分布的积分。然后D的协方差可通过利用根据以下公式的协方差的标准定义限定的4阶张量S得出：

$\underset{\&OverBar;}{S}=<D^{\&CircleTimes;2}>-<D{>}^{\&CircleTimes;2}---\left(35\right)$

其中为从D与其自身的外积得到的四阶张量。来自于一部分的弥散编码MR-信号E包括多个这种微观分体积，每个都具有高斯弥散，可通过以下公式估算：

E(b)＝<exp(-<b，D>)>， (36)

其中<·,·>为内积。为了便于理解，应当注意，方程(36)基于方程(6)，但是区别在于E表示归一化的回波信号强度(即S/S₀)并包括来自于所述部分的每一微观环境的信号分布。此外，为了简化符号，下文将采用内积取代方程(7)中所使用的广义上的数积。从方程(36)的累积展开为：

$E\left(b\right)\&ap;\exp (-<b_2<D>>+\frac{1}{2}<\underset{\&OverBar;}{B},\underset{\&OverBar;}{S}>)---\left(37\right)$

其中且S为所述部分内的张量的协方差(其在弥散MRI的情况下由体素表示)。为便于理解，下面提供累积展开的具体推导过程。

方程37中的近似值为累积展开，其中根据以下公式，当大约b＝0时展开e(b)＝logE(b)(其中b对应于总弥散加权，类似于方程9)：

$e\left(b\right)=log<\exp (-b<N,D>)>\&ap;f\left(0\right)+{\mathrm{bf}}^{\&prime;}\left(0\right)+\frac{1}{2}{b}^2{f}^{\&prime;\&prime;}\left(0\right)---\left(38\right)$

其中

$f^{\&prime;}\left(b\right)=\frac{E^{\&prime;}\left(b\right)}{E\left(b\right)}---\left(39\right)$

$f^{\&prime;\&prime;}\left(b\right)=\frac{E^{\&prime;\&prime;}\left(b\right)}{E\left(b\right)}-{\left(\frac{E^{\&prime;}\left(b\right)}{E\left(b\right)}\right)}^2---\left(40\right)$

对于b＝0，这些函数求值得到：

E(0)＝1 (41)

E′(0)＝-<<N，D>>＝-<N，<D>> (42)

E″(0)＝<<N，D>²>＝<<N，D>>＝<N，<D>> (43)

其中因此，

在弥散张量D的情形中，两个共同的不变量表示为平均弥散率(MD)和部分各向异性(FA)。MD可根据以下公式作为在各向同性基E_iso上的投影进行计算，

MD(D)＝<D，E_iso>， (44)

其中E_iso＝1/3I，即单位张量(identity tensor)的三分之一。这相当于结合方程(1)定义的参数D_iso。类似地，四阶协方差张量S可投影到各向同性基上，以得到旋转不变参数。但是，各向同性的4阶张量S具有两个各向同性的分量，其与力学领域中类似，可理解为4阶应力张量的体积模量和剪切模量(参见例如参考文献[28])。这两个基可通过以下公式定义：

和

其中I为单位张量。注意，E ₁＝E _体积和E ₂＝E _剪切是正交的，即

<E _i，E _j>＝δ_ij. (47)

其中如果i＝j，则δ_ij为1，否则为0。

正如在方程(44)中一样，类似地估算平均弥散率MD，4阶协方差张量可投影到其两个各向同性的基元素上。根据以下公式，投影至E _体积上得出平均弥散率V_MD(其也可以表示为S_体积)的协方差：

<S，E _体积>＝<MD(D)²>-MD(<D>)²＝V_MD· (48)

这由以下方程(49-53)推导出来。

$=<<D^{\&CircleTimes;2}>,E_{iso}^{\&CircleTimes;2}>-<<D{>}^{\&CircleTimes;2},E_{iso}^{\&CircleTimes;2}>=---\left(50\right)$

$=<<D^{\&CircleTimes;2},E_{iso}^{\&CircleTimes;2}>>-<<D>,E_{iso}{>}^2=---\left(51\right)$

＝<<D，E_iso>²>-<<D，E_iso>>²＝ (52)

＝<MD²>-<MD>²＝V_MD (53)

将S投影至E _剪切上得出与张量本征值的协方差相关的另一不变参数，<S，E _剪切>＝<V_λ(D)>-V_λ(<D>)＝ΔV_λ， (54)

(其也表示为S_剪切)，其中V_λ(·)得出弥散张量本征值(λ_i,i＝1,2,3)的协方差，

$=\frac{1}{3}{\mathrm{\&Sigma;}}_{i=1}^3{\mathrm{\&lambda;}}_i^2-\frac{1}{3}{\left({\mathrm{\&Sigma;}}_{i=1}^3{\mathrm{\&lambda;}}_i\right)}^2.---\left(56\right)$

这由以下方程(57-64)推导出来。考虑到映射不是S而是在E _剪切上的投影，

＝<D²，E_iso>-<D，E_iso>²＝ (59)

$=\frac{1}{3}{\mathrm{\&Sigma;}}_{i=1}^3{\mathrm{\&lambda;}}_i^2-{\left(\frac{1}{3}{\mathrm{\&Sigma;}}_{i=1}^3{\mathrm{\&lambda;}}_i\right)}^2=V_{\mathrm{\&lambda;}}\left(D\right)---\left(60\right)$

其中利用了以下关系：

$<D^2,E_{iso}>=\frac{1}{3}Tr\left(DD\right)=\frac{1}{3}<D,D>=<D^{\&CircleTimes;2},\frac{1}{3}\underset{\&OverBar;}{I}>---\left(61\right)$

基于上文，S的投影遵守：

＝<V_λ(D)>-V_λ(<D>) (64)

正如本发明人所理解，V_MD可解释为弥散张量的体积变化(即大小的变化)，ΔV_λ可解释为它们的剪切(即，方向之间的变化)。

在常规的弥散张量成像(DTI)中，通常考虑不变参数是部分各向异性(FA)。其由平均张量的本征值的归一化方差定义(参见例如参考文献[29])。

$FA(<D>)=\sqrt{\frac{3}{2}\frac{\frac{1}{3}{\mathrm{\&Sigma;}}_{i=1}^3{\mathrm{\&lambda;}}_i^2-{(\frac{1}{3}{\mathrm{\&Sigma;}}_{i=1}^3-{\mathrm{\&lambda;}}_i)}^2}{\frac{1}{3}{\mathrm{\&Sigma;}}_{i=1}^3{\mathrm{\&lambda;}}_i^2}}---\left(65\right)$

通过使用方程(44-48)和(54-56)，根据下式，FA可以表达为常规平均弥散张量(通过取与其自身的外积升至更高的阶)在体积基和剪切基上的投影：

因此，通过将平均弥散张量的方差(V_λ(<D>))替换为微观区室张量(即样品的所述部分的分体积的弥散张量)的本征值的平均方差(<V_λ(D)>)，根据下式可以得到微观FA(μFA)：

其中若微观张量共享同一组的本征值，则它们共享同一V_λ(D)值，在这种情况下，μFA将产生微观张量的精确FA。

通过利用方程(67)计算μFA，由于在μFA计算中，外积作用在局部弥散张量上，而不是全局平均张量上，该参数变得对取向分散不再敏感。这也使得实施和估算简单。

张量的沃伊特符号(Voigt notation)

出于实施目的，张量D和S有利地用沃伊特符号表示，这使得D可以表示为大小为6×1的列向量：

$d={\left(D_{xx}{D}_{yy}{D}_{zz}\sqrt{2}{D}_{yz}\sqrt{2}{D}_{xz}\sqrt{2}{D}_{xy}\right)}^T.---\left(68\right)$

现在第四阶(4阶)张量S可以用6×6方差-协方差矩阵表示，由于其根据下式由d定义：

S＝<dd^T>-<d><d>^T， (69)

。完全形式下，S由下式得出：

$\underset{\&OverBar;}{S}=\left(\begin{array}{cccccc}{\mathrm{\&Sigma;}}_{xxxx}& {\mathrm{\&Sigma;}}_{xxyy}& {\mathrm{\&Sigma;}}_{xxzz}& \sqrt{2}{\mathrm{\&Sigma;}}_{xxyz}& \sqrt{2}{\mathrm{\&Sigma;}}_{xxxz}& \sqrt{2}{\mathrm{\&Sigma;}}_{xxxy}\\ {\mathrm{\&Sigma;}}_{yyxx}& {\mathrm{\&Sigma;}}_{yyyy}& {\mathrm{\&Sigma;}}_{yyzz}& \sqrt{2}{\mathrm{\&Sigma;}}_{yyyz}& \sqrt{2}{\mathrm{\&Sigma;}}_{yyxz}& \sqrt{2}{\mathrm{\&Sigma;}}_{yyxy}\\ {\mathrm{\&Sigma;}}_{zzxx}& {\mathrm{\&Sigma;}}_{zzyy}& {\mathrm{\&Sigma;}}_{zzzz}& \sqrt{2}{\mathrm{\&Sigma;}}_{zzyz}& \sqrt{2}{\mathrm{\&Sigma;}}_{zzxz}& \sqrt{2}{\mathrm{\&Sigma;}}_{zzxy}\\ \sqrt{2}{\mathrm{\&Sigma;}}_{yzxx}& \sqrt{2}{\mathrm{\&Sigma;}}_{yzyy}& \sqrt{2}{\mathrm{\&Sigma;}}_{yzzz}& 2{\mathrm{\&Sigma;}}_{yzyz}& 2{\mathrm{\&Sigma;}}_{yzxz}& 2{\mathrm{\&Sigma;}}_{yzxy}\\ \sqrt{2}{\mathrm{\&Sigma;}}_{xzxx}& \sqrt{2}{\mathrm{\&Sigma;}}_{xzyy}& \sqrt{2}{\mathrm{\&Sigma;}}_{xzzz}& 2{\mathrm{\&Sigma;}}_{xzyz}& 2{\mathrm{\&Sigma;}}_{xzxz}& 2{\mathrm{\&Sigma;}}_{xzxy}\\ \sqrt{2}{\mathrm{\&Sigma;}}_{xyxx}& \sqrt{2}{\mathrm{\&Sigma;}}_{xyyy}& \sqrt{2}{\mathrm{\&Sigma;}}_{xyzz}& 2{\mathrm{\&Sigma;}}_{xyyz}& 2{\mathrm{\&Sigma;}}_{xyxz}& 2{\mathrm{\&Sigma;}}_{xyxy}\end{array}\right)---\left(70\right)$

四阶张量(例如B)或者可以用6×6矩阵(与方程(70)类似)表示，或者根据下式可以沃伊特符号由21×1列向量b表示：

$\underset{\&OverBar;}{b}=\left(\begin{array}{c}\sqrt{1}{b}_{xx}{b}_{xx}\\ \sqrt{1}{b}_{yy}{b}_{yy}\\ \sqrt{1}{b}_{zz}{b}_{zz}\\ \sqrt{2}{b}_{yy}{b}_{zz}\\ \sqrt{2}{b}_{xx}{b}_{zz}\\ \sqrt{2}{b}_{xx}{b}_{yy}\\ \sqrt{4}{b}_{xx}{b}_{yz}\\ \sqrt{4}{b}_{xx}{b}_{xz}\\ \sqrt{4}{b}_{xx}{b}_{xy}\\ \sqrt{4}{b}_{yy}{b}_{yz}\\ \sqrt{4}{b}_{yy}{b}_{xz}\\ \sqrt{4}{b}_{yy}{b}_{xy}\\ \sqrt{4}{b}_{zz}{b}_{yz}\\ \sqrt{4}{b}_{zz}{b}_{xz}\\ \sqrt{4}{b}_{zz}{b}_{xy}\\ \sqrt{4}{b}_{yz}{b}_{yz}\\ \sqrt{4}{b}_{xz}{b}_{xz}\\ \sqrt{4}{b}_{xy}{b}_{xy}\\ \sqrt{8}{b}_{yz}{b}_{xz}\\ \sqrt{8}{b}_{xz}{b}_{xy}\\ \sqrt{8}{b}_{xy}{b}_{yz}\end{array}\right)---\left(71\right)$

使用沃伊特符号时，体积基和剪切基可以表示为：

和

如上所述，这两个各向同性体积基和剪切基是利用正交设计，如能够从方程(45)和(46)看到的，将这些基函数加和可以得出结构简单的对角单位矩阵。

内积和外积

在实现过程中使用矩阵和向量表示张量的一个优点在于内积和外积易于在软件上实现。张量(例如，D)的外积可以根据下式计算：

$D^{\&CircleTimes;2}={\mathrm{dd}}^T---\left(74\right)$

因此，可以定义例如或者用沃伊特符号

内积用<·,·>表示，根据下式，其可以定义为按元素相乘，接着再相加：

<D，N>＝d^Tn (75)

或

<S，B>s ^T b (76)

根据下式，两个矩阵的内积还可以定义为：

<D，N>＝Tr(DN^T) (77)

实施例描述

考虑到上述内容，根据优选实施例，利用多个不同的弥散编码磁梯度脉冲序列在样品上进行多个弥散加权回波衰减测量，其中每个磁梯度脉冲序列G_i被生成为使得对于每个磁梯度脉冲序列G_i具有一个弥散编码张量b_i。多个脉冲序列可以包括其弥散编码张量具有1至3个(优选地，具有2至3个)非零本征值的脉冲序列的组合。上述关于在NMR分光仪或MRI装置上实施该方法的讨论也可以应用在本发明的实施例上。

通常情况下，通过时间依赖性梯度建立q-向量以穿过q-空间内的一任意路径。弥散编码张量的阶(即非零本征值的数目)取决于该路径，在sPFG的情况下为1，当第一和第二梯度脉冲沿正交方向施加时为2，在各向同性编码的情况下(比如三重PFG[4]或q-MAS[30])为3。例如，平面弥散编码张量(即在平面内旋转对称的编码)可以通过产生平面q-空间轨迹的一组时变梯度来实现。由于b-值是时间和q-值两者的函数，通过在小q-值时采用较低速度，能够通过改变在q-空间中的遍历速度保证不变的角b-值。在小q时，长弥散时间能建立相同的编码功率(b-值)，q-值越大，弥散时间越短。

利用脉冲序列G_i获得的表示回波衰减测量值的数据组{E₁,…E_m}可以被收集并排列在列向量表示中。b和<D>可以(使用沃伊特符号)排列在大小为6×1的列向量表示中，用b_V和<d>表示；类似地，四阶张量B和S可以排列在大小为21×1的列向量表示中，由b和s表示。利用方程(5)和(8)中的定义可以获得每个b_i的元素。根据下式，现在方程(37)中的内积可以通过简单矩阵运算表达：

$<b,<D>>=b_V^T<d>---\left(78\right)$

和

<B，S>＝b ^T s (79)

由于方程(37)是线性模型，可以利用伪逆矩阵对下述方程组求解，从而对<d>和s进行估算：

$\left(\begin{array}{c}log{E}_1\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ \log E_m\end{array}\right)=\left(\begin{array}{ccc}1& -b_{V1}^T& \frac{1}{2}{\underset{\&OverBar;}{b}}_1^T\\ \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;\\ 1& -b_{Vm}^T& \frac{1}{2}{\underset{\&OverBar;}{b}}_m^T\end{array}\right){\left(\begin{array}{ccc}{E}_0& <d>& \underset{\&OverBar;}{s}\end{array}\right)}^T---\left(80\right)$

基于方程(36)的累积展开，方程(80)形成线性方程组。数据{E₁,…E_m}形成该线性方程组的常数。脉冲序列G_i的弥散编码张量表示b_v1…_m及其四阶张量表示b₁…_m形成该线性方程组的参数。

在E₀、<d>和s中，该模型总共具有1+6+21＝28个自由参数。为了能够使用伪逆矩阵估算方程(80)的解，应当利用变化形状的测量张量获得数据，这样体积分量(b ^T _K e _体积)和剪切分量(b ^T _K e _剪切)之间的相关性小于1。这解释了为什么利用常规的sPFG不能分开S的两个各向同性分量。假定编码张量的体积分量和剪切分量不是完全相关的，若测量的数量超过28，可以执行伪逆矩阵。但是，如果仅求<d>和s的投影，则可以使用较少量的测量。

在上文中，本发明的构思主要参考有限数量的示例进行了描述。但是，本领域的技术人员能容易理解的是，在由所附权利要求书限定的本发明的构思的范围内，除了上述公开的那些以外的其他示例也是同样可能的。

参考文献列表

在上述公开中，成对括号“[]”之间的一个或多个数字指下述参考文献列表中相应编号的参考文献。

[1]P.T.Callaghan,O.Examination of the lamellar phase ofAerosol OT/water using pulsed field gradient nuclear magnetic resonance,J.Phys.Chem.,87(1983)1737-1744(P.T.Callaghan，O.利用脉冲场梯度核磁共振对气溶胶OT/水的层状相的检查，《物理化学杂志》，1983年第87卷第1737-1744页)。

[2]M.E.Moseley,J.Kucharczyk,H.S.Asgari,D.Norman,Anisotropy indiffusion-weighted MRI,Magn.Reson.Med.,19(1991)321-326(M.E.Moseley，J.Kucharczyk，H.S.Asgari，D.Norman，弥散加权MRI中的各向异性，《医学磁共振》，1991年第19卷第321-326页)。

[3]E.O.Stejskal,J.E.Tanner,Spin diffusion measurements:Spin echoes inthe presence of a time-dependent field gradient,J.Chem.Phys.,42(1965)288-292(E.O.Stejskal，J.E.Tanner，自旋弥散测量：时间依赖性场梯度存在下的自旋回波，《物理化学杂志》，1965年第42卷第288-292页)。

[4]S.Eriksson,S.Lasic,D.Topgaard,Isotropic diffusion weighting inPGSE NMR by magic-angle spinning of the q-vector,J.Magn.Reson.,226(2013)13-18(S.Eriksson，S.Lasic,D.Topgaard，利用q-向量魔角旋转的PGSE NMR各向同性弥散加权，《磁共振杂志》，2013年第226卷第13-18页)。

[5]L.Frydman,G.C.Chingas,Y.K.Lee,P.J.Grandinetti,M.A.Eastman,G.A.Barrall,A.Pines,Variable-angle correlation spectroscopy in solid-statenuclear magnetic resonance,J.Chem.Phys.,97(1992)4800-4808(L.Frydman，G.C.Chingas，Y.K.Lee，P.J.Grandinetti，M.A.Eastman，G.A.Barrall，A.Pines，固态核磁共振中的可变角关联光谱学，《物理化学杂志》，1992年第97卷第4800-4808页)。

[6]A.Bax,N.M.Szeverenyi,G.E.Maciel,Chemical shift anisotropy inpowdered solids studied by 2D FT NMR with flipping of the spinning axis,J.Magn.Reson.,55(1983)494-497(A.Bax，N.M.Szeverenyi,G.E.Maciel，利用带旋转轴翻转的2D FT NMR研究的粉末固体中的化学位移各向异性，《磁共振杂志》，1983年第55卷第494-497页)。

[7]Y.-Q.Song,L.Venkataramanan,M.D.Hürlimann,M.Flaum,P.Frulla,C.Straley,T1-T2correlation spectra obtained using a fast two-dimensionalLaplace inversion,J.Magn.Reson.,154(2002)261-268(Y.-Q.Song，L.Venkataramanan，M.D.Hürlimann，M.Flaum，P.Frulla，C.Straley，利用快速二维拉普拉斯反演获得的T1-T2关联光谱，《磁共振杂志》，2002年第154卷第261-268页)。

[8]M.Lustig,D.Donoho,J.M.Pauly,Sparse MRI:The application ofcompressed sensing for rapid MR imaging,Magn.Reson.Med.,58(2007)1182-1195(M.Lustig，D.Donoho，J.M.Pauly，稀疏MRI：用于快速MR成像的压缩感知的应用，《医学磁共振》，2007年第58卷第1182-1195页)。

[9]M.Mobli,M.W.Maciejewski,A.D.Schuyler,A.S.Stern,J.C.Hoch,Sparsesampling methods in multidimensional NMR,Phys.Chem.Chem.Phys.,14(2012)10835-10843(M.Mobli，M.W.Maciejewski，A.D.Schuyler，A.S.Stern，J.C.Hoch，多维NMR中的稀疏采样方法，《物理化学与化学物理》，2012年第14卷第10835-10843页)。

[10]D.Topgaard,Isotropic diffusion weighting in PGSE NMR:Numericaloptimization of the q-MAS PGSE sequence,Microporous Mesoporous Mater.,178(2013)60-63(D.Topgaard，PGSE NMR中的各向同性弥散加权：q-MAS PGSE序列的数值优化，《微孔与介孔材料》，2013年第178卷第60-63页)。

[11]P.J.Basser,J.Mattiello,D.Le Bihan,MR diffusion tensorspectroscopy and imaging,Biophys.J.,66(1994)259-267(P.J.Basser，J.Mattiello，D.Le Bihan，MR弥散张量光谱学和成像，《生物物理学杂志》，1994年第66卷第259-267页)。

[12]K.Schmidt-Rohr,H.W.Spiess,Multidimensional solid-state NMR andpolymers,Academic Press,San Diego,(1994)(K.Schmidt-Rohr，H.W.Spiess，多维固态NMR和聚合物，学术出版社，圣地亚哥，1994年)。

[13]M.Duer,Introduction to solid-state NMR spectroscopy,BlackwellPublishing Ltd,Oxford,(2004)(M.Duer，固态NMR光谱学简介，布莱克威尔出版公司，牛津，2004年)。

[14]W.S.Price,NMR studies of translational motion,CambridgeUniversity Press,Cambridge,(2009)(W.S.Price，平移运动的NMR研究，剑桥大学出版社，剑桥，2009年)。

[15]P.T.Callaghan,Translational dynamics&magnetic resonance,OxfordUniversity Press,Oxford,(2011)(P.T.Callaghan，平移动力学与磁共振，牛津大学出版社，牛津，2011年)。

[16]G.B.Arfken,H.-J.Weber,Mathematical methods for physicists,4thed.,Academic Press,San Diego,(1995)(G.B.Arfken，H.-J.Weber，《物理学家的数学方法》，第4版，学术出版社，圣地亚哥，1995年)。

[17]V.J.Wedeen,D.L.Rosene,R.Wang,G.Dai,F.Mortazavi,P.Hagmann,J.H.Kaas,W.-Y.I.Tseng,The geometric structure of the brain fiber pathways,Science,335(2012)1628-1634(V.J.Wedeen，D.L.Rosene，R.Wang，G.Dai，F.Mortazavi，P.Hagmann，J.H.Kaas，W.-Y.I.Tseng，脑纤维通路的几何结构，《科学》，2012年第335卷第1628-1634页)。

[18]K.P.Whittal,Analysis of large one-dimensional and two-dimensionalrelaxation data sets,J.Magn.Reson.A,110(1994)214-218(K.P.Whittal，大型一维和二维弛豫数据组分析，《磁共振杂志》，1994年第110卷第214-218页)。

[19]D.Bernin,D.Topgaard,NMR diffusion and relaxation correlationmethods:New insights in heterogeneous materials,Curr.Opin.Colloid InterfaceSci.,18(2013)166-172(D.Bernin，D.Topgaard，NMR弥散和弛豫关联方法：非均质材料的新见解，胶体与界面科学新见，2013年第18卷第166-172页)。

[20]J.Mitchell,T.C.Chandrasekera,L.F.Gladden,Numerical estimation ofrelaxation and diffusion distributions in two dimensions,Prog.Nucl.Magn.Reson.Spectrosc.,62(2012)34-50(J.Mitchell，T.C.Chandrasekera，L.F.Gladden，二维弛豫分布和弥散分布的数值估算，《核磁共振光谱学进展》，2012年第62卷第34-50页)。

[21]B.B.Lindman,K.Holmberg,B.Kronberg,Surfactants andpolymers in aqueous solution,John Wiley&Sons Ltd,(1998)(B.B.Lindman，K.Holmberg，B.Kronberg，水溶液中的表面活性剂和聚合物，约翰威立国际出版公司，1998年)。

[22]M.Bak,N.C.Nielsen,REPULSION,A novel approach to efficient powderaveraging in solid-state NMR,J.Magn.Reson.,125(1997)132-139(M.Bak，N.C.Nielsen，排斥：固态NMR中有效粉末平均化的新方法，《磁共振杂志》，1997年第125卷第132-139页)。

[23]D.K.Jones,M.A.Horsfield,A.Simmons,Optimal strategies formeasuring diffusion in anisotropic systems by magnetic resonance imaging,Magn.Reson.Med.,42(1999)515-525(D.K.Jones，M.A.Horsfield，A.Simmons，利用磁共振成像在各向异性系统中测量弥散的最佳策略，《医学磁共振》，1999年第42卷第515-525页)。

[24]Le Bihan,D.,J.B.H.:Diffusion MRI at 25:exploring brain tissuestructure and function(2012)(Le Bihan D.，J.B.H.，25处的弥散MRI：探索大脑组织结构和功能，2012年)。

[25]Assaf,Y.,P.O.:Diffusion MRI at 25:exploring brain tissuestructure and function(2008)(Assaf,Y.，P.O.，25处的弥散MRI：探索大脑组织结构和功能，2012年)。

[26]Jones,D.K.,K.T.R.T.R.:White matter integrity,fiber count,andother fallacies:the dos and don’t’s of diffusion mri(2013)(Jones,D.K.，K.T.R.，T.R.，白质完整性、纤维支数和其他谬误：弥散MRI的注意事项，2013年)。

[27]Vos,S.B.,J.D.K.J.B.V.M.A.L.A.:The influence of complex whitematter architecture on the mean diffusivity in diffusion tensor mri of thehuman brain(2012)(Vos，S.B.,J.D.K.，J.B.，V.M.A.，L.A.，人脑弥散张量MRI中复杂白质体系结构对平均弥散率的影响，2012年)。

[28]Moakher,M.:Fourth-order cartesian tensors:old and new facts,notions and applications.The Quart.Jour.of Mechanics and Applied Math.61(2)(2008)181-203(Moakher，M.，四阶笛卡尔张量：新旧事实、概念和应用，《力学与应用数学季刊》，2008年第61卷第2期第181-203页)。

[29]Basser,P.J.,Pierpaoli,C.:Microstructural and physiologicalfeatures of tissues elucidated by quantitative-diffusion-tensor mri.Journalof Magnetic Resonance,Series B 111(3)(1996)209-219(Basser,P.J.,Pierpaoli,C.，利用定量弥散张量MRI描述的组织微结构特征和生理学特征，《磁共振杂志B辑》，1996年第111卷第3期第209-219页)。

[30]Valette,J.,Giraudeau,C.,Marchadour,C.,Djemai,B.,Geffroy,F.,Ghaly,M.A.,Le Bihan,D.,Hantraye,P.,Lebon,V.,Lethimonnier,F.:A new sequence forsingle-shot diffusion-weighted nmr spectroscopy by the trace of the diffusiontensor.MRM 68(6)(2012)1705-1712(Valette，J.，Giraudeau，C.，Marchadour,C.，Djemai,B.，Geffroy,F.，Ghaly,M.A.，Le Bihan,D.，Hantraye,P.，Lebon,V.，Lethimonnier，用于利用弥散张量轨迹的单射击弥散加权NMR光谱的新序列，《MRM》，2012年第68卷第6期第1705-1712页)。

Claims (20)

1.一种用于定量样品中的各向同性弥散和/或各向异性弥散的方法，所述方法包括：

利用弥散编码磁梯度脉冲序列G_i＝1…m对所述样品进行弥散加权磁共振测量，其中每个磁梯度脉冲序列G_i被生成为使得用于所述磁梯度脉冲序列G_i的弥散编码张量b_i具有一至三个非零本征值，其中q_i(t)是与成正比的移相向量，并且τ是回波时间；

收集表示磁共振回波信号测量值的数据，所述数据的至少一个子组表示利用引起各向异性弥散加权的一组磁梯度脉冲序列得到的回波信号，其中用于所述磁梯度脉冲序列组中的每个梯度脉冲序列的弥散编码张量具有三个非零本征值，所述三个本征值中的至少一个与另外两个本征值不同；以及

利用所述数据计算各向同性弥散程度和/或各向异性弥散程度。

2.根据权利要求1所述的方法，其中所述数据子组表示来自一部分所述样品的回波信号，所述部分包括多个分体积，所述多个分体积具有不同程度的各向同性弥散或不同程度和/或不同取向的各向异性弥散，其中计算各向同性弥散程度和/或各向异性弥散程度包括针对所述分体积中的至少一个计算各向同性弥散程度的估计值和/或各向异性弥散程度的估计值。

3.根据权利要求1至2中任一项所述的方法，

其中所述磁梯度脉冲序列组形成第一组磁梯度脉冲序列，并且所述数据子组形成利用所述第一组磁梯度脉冲序列得到的表示第一回波衰减曲线的第一数据子组，并且

其中所述数据还包括至少第二数据子组，所述第二数据子组表示利用引起各向同性或各向异性弥散加权的第二组磁梯度脉冲序列得到的第二回波衰减曲线。

4.根据权利要求3所述的方法，

其中所述第一组中的每个脉冲序列使得用于所述脉冲序列的弥散编码张量的第一本征值和第二本征值彼此相等，并且

其中所述第二组中的每个脉冲序列使得用于所述脉冲序列的弥散编码张量的第一本征值和第二本征值彼此相等。

5.根据权利要求4所述的方法，

其中所述第一组磁梯度脉冲和所述第二组磁梯度脉冲序列中的脉冲序列具有不同的最大梯度幅值。

6.根据权利要求3至5中任一项所述的方法，

其中，对于所述第一组中的每个脉冲序列，存在能由下式定义的第一弥散编码张量不变量Δ_b,1：

${\mathrm{\&Delta;}}_{b,1}=\frac{1}{b}(b_{zz}^{PAS}-\frac{b_{yy}^{PAS}+b_{xx}^{PAS}}{2}),b=b_{xx}^{PAS}+b_{yy}^{PAS}+b_{zz}^{PAS}$

其中，b_xx ^PAS表示用于所述脉冲序列的弥散编码张量的第一本征值，b_yy ^PAS表示用于所述脉冲序列的弥散编码张量的第二本征值，而b_zz ^PAS表示用于所述脉冲序列的弥散编码张量的第三本征值，并且

其中所述第一组脉冲序列使得所述第一组中的所述脉冲序列的所述第一弥散编码张量不变量Δ_b,1彼此相等。

7.根据权利要求6所述的方法，

其中，对于所述第二组中的每个脉冲序列，存在能由下式定义的第二弥散编码张量不变量Δ_b,2：

${\mathrm{\&Delta;}}_{b,2}=\frac{1}{b}(b_{zz}^{PAS}-\frac{b_{yy}^{PAS}+b_{xx}^{PAS}}{2}),b=b_{xx}^{PAS}+b_{yy}^{PAS}+b_{zz}^{PAS}$

其中，b_xx ^PAS表示用于所述脉冲序列的弥散编码张量的第一本征值，b_yy ^PAS表示用于所述脉冲序列的弥散编码张量的第二本征值，而b_zz ^PAS表示用于所述脉冲序列的弥散编码张量的第三本征值，并且

其中所述第二组脉冲序列使得所述第二组中的所述脉冲序列的所述第二弥散编码张量不变量Δ_b,2彼此相等，并且Δ_b,2不同于Δ_b,1。

8.根据权利要求3至7中任一项所述的方法，其中计算各向同性弥散程度和/或各向异性弥散程度包括：

通过分析利用所述第一组中的脉冲序列获得的第一回波信号与利用所述第二组中的脉冲序列获得的第二回波信号之间的变化、改变或差异，来计算各向同性弥散程度和/或各向异性弥散程度。

9.根据权利要求3至8中任一项所述的方法，其中，除了所述数据的第一子组和所述数据的第二子组外，所述数据至少包括第三数据子组，所述第三数据子组利用引起各向异性弥散加权的第三组磁梯度脉冲序列而获得，

其中用于所述第三组中的每个梯度脉冲序列的弥散编码张量具有3个非零本征值，其中第一本征值和第二本征值彼此相等且与第三本征值不同，并且

其中，对于所述第三组中的每个脉冲序列，存在能由下式定义的第三弥散编码张量不变量Δ_b,3：

${\mathrm{\&Delta;}}_{b,3}=\frac{1}{b}(b_{zz}^{PAS}-\frac{b_{yy}^{PAS}+b_{xx}^{PAS}}{2}),b=b_{xx}^{PAS}+b_{yy}^{PAS}+b_{zz}^{PAS}$

其中b_xx ^PAS表示用于所述脉冲序列的弥散编码张量的第一本征值，b_yy ^PAS表示用于所述脉冲序列的弥散编码张量的第二本征值，而b_zz ^PAS表示用于所述脉冲序列的弥散编码张量的第三本征值，并且

其中所述第三组脉冲序列使得所述第三组中的所述脉冲序列的所述第三弥散编码张量不变量Δ_b,3彼此相等，并且Δ_b,3不同于Δ_b,2和Δ_b,1。

10.根据权利要求9所述的方法，其中所述第一组中的每个脉冲序列使得Δ_b,1>0，所述第二组中的每个脉冲序列使得Δ_b,2＝0，而所述第三组中的每个脉冲序列使得Δ_b,3<0。

11.根据权利要求3至10中任一项所述的方法，还包括：

基于表示所述回波信号的数据计算概率分布，所述概率分布指示所述回波信号中的每一个与模型各向同性弥散参数D_iso和/或模型各向异性弥散参数Δ_D的多个不同值中的每一个相关的概率。

12.根据权利要求11所述的方法，其中所述概率分布通过确定方程组的解而计算得到，所述方程组将由所述数值表示的所述回波信号与核函数和所述概率分布的乘积关联起来。

13.根据权利要求12所述的方法，其中所述概率分布为联合概率分布，并且所述核函数为包括至少M×N个元素的矩阵，对于弥散加权幅值b、弥散编码张量不变量Δ_b、所述模型各向同性弥散参数D_iso和所述模型各向异性弥散参数Δ_D的组合值，每个所述元素基于下述公式的积分：

其中A＝3bD_isoΔ_bΔ_D。

14.根据前述权利要求中任一项所述的方法，还包括：

将所述第一组磁梯度脉冲序列中的每个脉冲序列分多次施加到所述样品上，相对于固定的实验室系采用不同的梯度脉冲取向，并且

通过对针对所述不同取向得到的回波衰减测量值求平均以形成所述第一数据子组。

15.根据前述权利要求中任一项所述的方法，其中每一个所述弥散编码磁梯度脉冲序列G_i形成三重受激回波序列中的一部分。

16.根据权利要求1至3中任一项所述的方法，还包括：

基于将回波信号E与弥散编码张量b和弥散张量D关联起来的函数展开形成方程组；

通过利用由所述数据表示的回波信号测量值和弥散编码张量b_i的至少一个子组的表示值确定所述方程组的解，来计算平均弥散张量<D>和弥散张量协方差张量S；

通过将S投影到体积基E _体积上，计算协方差张量S的不变体积分量S_体积；

通过将S投影到剪切基E _剪切上，计算协方差张量S的不变剪切分量S_剪切；

利用所述不变体积分量S_体积和/或所述不变剪切分量S_剪切计算各向同性弥散程度和/或各向异性弥散程度。

17.根据权利要求16所述的方法，其中所述方程组等同于函数E(b)＝<exp(–<b,D>)>的累积展开。

18.根据权利要求16至17中任一项所述的方法，其中基于所述不变剪切分量S_剪切和所述平均弥散张量<D>的平方到所述剪切基E _剪切上的投影的总和，计算所述各向异性弥散程度。

19.根据权利要求18所述的方法，其中基于所述平均弥散张量<D>的平方到所述体积基E _体积上的投影与所述总和之间的比值，计算所述各向异性弥散程度。

20.根据权利要求19所述的方法，其中基于所述比值，所述各向异性弥散程度被计算为微观部分各向异性μFA的估算值。