CN104599244A - 磁共振弥散张量成像的去噪方法和*** - Google Patents

Abstract

本发明提供了一种磁共振弥散张量成像的去噪方法和***，其方法基于磁共振弥散加权成像模型和采样噪声的高斯分布性质，利用平均弥散率的稀疏性，采用最大后验概率估计的方法直接由K空间数据获得每一个空间位置所对应的去噪后的弥散张量矩阵。本发明可以避免图像去噪的误差对弥散张量估计的影响，可以更有效地抑制弥散张量中的噪声，提高弥散张量的估计精度。

Description

磁共振弥散张量成像的去噪方法和***

技术领域

本发明涉及了磁共振成像技术，特别是涉及一种磁共振弥散张量成像的去噪方法和***。

背景技术

弥散张量成像(Diffusion Tensor Imaging DTI)是在弥散加权成像(Diffusion Weighted Imaging DWI)基础上发展起来的新的成像方法，其利用水分子的弥散各向异性进行成像，可无损的从微观领域评价组织结构的完整性，为疾病的预防、诊断及治疗提供更多的信息。但是，相比其他的磁共振成像技术，DTI需要较长的扫描时间，且信噪比较低。

为了提高弥散张量成像的信噪比，目前较为直接的方法是通过多次采样取平均和减小K空间采样区域的方法。这些方法在实际中有一定的应用，但会增加扫描时间和影响空间分辨率。另一种常见的方法是在获取K空间扫描数据后，首先重建出弥散加权图像，然后再运用信号处理的方法对图像进行去噪，最后，通过去噪后的图像计算弥散张量及各种弥散参数，用以生成弥散张量成像图。该方法应用较广，但图像去噪过程中的***误差有可能会进一步传递到后续的张量计算中，进而影响各种弥散参数图像的质量。

发明内容

基于此，有必要针对现有技术中的问题，提供一种磁共振弥散张量成像的去噪方法和***，其可以避免图像去噪的误差对弥散张量估计的影响，可以更有效地抑制弥散张量中的噪声，提高弥散张量的估计精度。

本发明提供了一种磁共振弥散张量成像的去噪方法，其包括：

图像数据获取步骤：获取磁共振弥散加权图像所对应的K空间数据；

去噪步骤：基于磁共振弥散加权成像模型和采样噪声的高斯分布性质，利用平均弥散率的稀疏性，采用最大后验概率估计的方法由所述K空间数据获得每一个空间位置所对应的去噪后的弥散张量矩阵；

弥散参数计算步骤：基于所述去噪后的弥散张量矩阵，获得弥散参数图。

在其中一个实施例中，所述去噪步骤包括：

基于磁共振弥散加权成像模型和采样噪声的高斯分布性质，利用平均弥散率的稀疏性，采用最大后验概率估计的方法构建去噪函数模型，所述去噪函数模型参见如下述公式(1)：

   公式(1)

其中，表示弥散张量矩阵的估计值，MD为平均弥散率；R(·)是作用于平均弥散率的稀疏约束函数，λ为相应的正则化参数；d_m为第m个弥散加权图像所对应的K空间数据；F表示傅立叶编码矩阵；表示第m个弥散加权图像，其中，I₀表示无弥散加权的参考图像，为第m个弥散加权图像的相位，b是弥散加权因子，g_m是第m个弥散加权图像所对应的弥散梯度向量g_m＝(g_xm,g_ym,g_zm)^T；

利用所述K空间数据，求解所述公式(1)，获得每一个空间位置所对应的去噪后的弥散张量矩阵。

在其中一个实施例中，所述稀疏约束函数约束平均弥散率的稀疏性。

在其中一个实施例中，所述稀疏约束函数利用L₁范函来约束平均弥散率在稀疏变换域内的稀疏性。

在其中一个实施例中，所述稀疏约束函数通过调用以下公式来计算获得：

$R\left(\mathrm{MD}\right)={\left|\right|\mathrm{\Ψ}\left(\frac{D_1+D_2+D_3}{3}\right)\left|\right|}_1$

其中，Ψ表示稀疏变换运算符，‖·‖₁表示取L₁范数，D₁、D₂、D₃分别为弥散张量矩阵主对角线的元素。

基于上述方法，本发明还提供了一种磁共振弥散张量成像的去噪***，其包括：

图像数据获取模块，用于获取磁共振弥散加权图像所对应的K空间数据；

去噪模块，用于基于磁共振弥散加权成像模型和采样噪声的高斯分布性质，利用平均弥散率的稀疏性，采用最大后验概率估计的方法由所述K空间数据获得每一个空间位置所对应的去噪后的弥散张量矩阵；及

弥散参数计算模块，用于基于所述去噪后的弥散张量矩阵，获得弥散参数图。

在其中一个实施例中，所述去噪模块包括：

模型构建单元，用于基于磁共振弥散加权成像模型和采样噪声的高斯分布性质，利用平均弥散率的稀疏性，采用最大后验概率估计的方法构建去噪函数模型，所述去噪函数模型参见如下述公式(1)：

   公式(1)

其中，表示弥散张量矩阵的估计值，MD为平均弥散率；R(·)是作用于平均弥散率的稀疏约束函数，λ为相应的正则化参数；d_m为第m个弥散加权图像所对应的K空间数据；F表示傅立叶编码矩阵；表示第m个弥散加权图像，其中，I₀表示无弥散加权的参考图像，为第m个弥散加权图像的相位，b是弥散加权因子，g_m是第m个弥散加权图像所对应的弥散梯度向量g_m＝(g_xm,g_ym,g_zm)^T；和

矩阵求解单元，用于利用所述K空间数据，求解所述公式(1)，获得每一个空间位置所对应的去噪后的弥散张量矩阵。

本发明提出利用平均弥散率的稀疏性，直接通过获取的K空间数据对弥散张量进行去噪成像的方法，该方法略过了传统的图像去噪的步骤，避免了图像去噪的误差对弥散张量估计的影响，可以更有效地抑制弥散张量中的噪声，提高弥散张量的估计精度。

附图说明

图1为本发明方法的一个实施例流程示意图；

图2为本发明方法的另一个实施例流程示意图；

图3为本发明***的一个实施例结构示意图；

图4为本发明***的另一个实施例结构示意图。

具体实施方式

本发明基于磁共振弥散加权成像技术，利用弥散张量模型和平均弥散率MD的固有特性(如稀疏性)，基于最大后验概率估计的理论框架，直接由采集的K空间数据获取去噪后的弥散张量，从而实现磁共振弥散张量成像的去噪。这里提到的磁共振弥散张量成像，是有别于弥散加权成像的一种新方法，可以用于描述大脑结构。举例来说，如果说核磁共振成像是追踪水分子中的氢原子，那么弥散张量成像便是依据水分子移动方向制图，弥散张量成像图(呈现方式与以前的图像不同)可以揭示脑瘤如何影响神经细胞连接，引导医疗人员进行大脑手术。它还可以揭示同中风、多发性硬化症、精神***症、阅读障碍有关的细微反常变化。

磁共振成像的物理机制可以表示为下述公式(2-1)：

d＝Fρ+n   公式(2-1)其中，d表示在磁共振仪上所采集的信号，F表示博立叶编码矩阵，ρ为磁共振重建的图像，n通常假定为复高斯白噪声。

在磁共振弥散加权成像中，第m个弥散加权图像ρ_m的成像模型可以写作下述公式(2-2)：

公式(2-2)其中，I₀表示无弥散加权的参考图像，为第m个弥散加权图像的相位，b是弥散加权因子(常量)，g_m是第m个弥散加权图像所对应的弥散梯度向量g_m＝(g_xm,g_ym,g_zm)^T。

对于弥散加权成像，只有一个弥散系数一个标量值来描述弥散属性，而在弥散张量成像中，弥散张量可以完全描述每一个方向上水分子的移动及在这些方向上水分子移动的相关性。而张量本质上就是一副三维空间的方向矢量图，可以用以显示脑图像中有方向性的白质纤维束内水分子移动的选择性。弥散张量通常由下述公式(2-3)矩阵表示，以下简称为弥散张量矩阵。

$D=\left(\begin{array}{ccc}{D}_1& D_4& D_5\\ D_4& D_2& D_6\\ D_5& D_6& D_3\end{array}\right)$    公式(2-3)

上述弥散张量矩阵D为对称矩阵，为了形象地表述弥散张量矩阵，可以进一步将弥散张量矩阵视为一个椭圆球体(ellipsoid)，其本征值代表了沿弥散椭球最大和最小轴的弥散系数。所以，弥散张量矩阵的三个本征值是最基本的旋转不变量(即不随弥散方向而改变)。它们是沿着三个坐标轴方向测量的主弥散系数。这三个坐标轴是组织固有的。每个本征值联系着一个本征向量，这个本征向量也是组织固有的。弥散张量矩阵的三个本征向量相互垂直，并构建了每个像素的局部参照纤维框架。在每个像素中，本征值从大到小排列：λ₁＝最大弥散系数，λ₂＝中级弥散系数，λ₃＝最低弥散系数。λ₁代表平行于纤维方向的弥散系数，λ₂和λ₃代表横向弥散系数。每个像素所对应的本征值λ₁、λ₂、λ₃以下统称为弥散张量特征值。

下面将结合附图详细说明本发明所提供的磁共振弥散张量成像的去噪方法和***的具体实施例。

如图1所示，本实施例提供的一种磁共振弥散张量成像的去噪方法，其包括以下步骤：

图像数据获取步骤110：获取磁共振弥散加权图像所对应的K空间数据；

去噪步骤120：基于磁共振弥散加权成像模型和采样噪声的高斯分布性质，利用平均弥散率的稀疏性，采用最大后验概率估计的方法由所述K空间数据获得每一个空间位置所对应的去噪后的弥散张量矩阵；

弥散参数计算步骤130：基于去噪后的弥散张量矩阵，获得弥散参数图。

基于上述实施例，上述步骤110中，所获得的K空间数据可以通过回波平面成像(echo planar imaging，EPI)技术和并行成像技术(parallel imaging)获得。当然本发明步骤110主要是为了获取K空间数据，而不限制所使用的采集技术。

基于上述实施例，如图2所示，上述去噪步骤120在具体实现时可以包括以下两个步骤。

步骤121，基于磁共振弥散加权成像模型(即上述公式(2-1)和(2-2))和采样噪声的高斯分布性质，利用平均弥散率MD的稀疏性，采用最大后验概率估计的方法构建去噪函数模型，此去噪函数模型参见如下述公式(1)：

   公式(1)

其中，表示弥散张量矩阵的估计值，MD为平均弥散率；R(·)是作用于平均弥散率的稀疏约束函数，λ为相应的正则化参数；d_m为第m个弥散加权图像所对应的K空间数据；F表示傅立叶编码矩阵；表示第m个弥散加权图像，其中，I₀表示无弥散加权的参考图像，为第m个弥散加权图像的相位，b是弥散加权因子，g_m是第m个弥散加权图像所对应的弥散梯度向量g_m＝(g_xm,g_ym,g_zm)^T；

步骤122，利用K空间数据，求解所述公式(1)，获得每一个空间位置所对应的去噪后的弥散张量矩阵。

上述步骤121中提到的MD，即平均弥散率(mean diffusivity，MD)，反映分子整体的弥散水平(平均椭球的大小)和弥散阻力的整体情况。MD只表示弥散的大小，而与弥散的方向无关。MD越大，组织内所含自由水分子则越多。计算公式如下述公式(2-4)所示：

MD＝(λ₁+λ₂+λ₃)/3   公式(2-4)其中，λ₁，λ₂，λ₃为弥散张量矩阵D的三个弥散张量矩阵特征值。

本发明基于采样噪声服从高斯分布，结合上述公式(2-1)和(2-2)，利用平均弥散率的固有特性(如稀疏性)，给出了上述公式(1)所述的去噪函数模型，用以表征弥散张量的最大后验概率估计。

而基于上述公式(2-4)，其中D₁、D₂、D₃分别为弥散张量矩阵主对角线的元素，于是上述公式(1)可以转变为下述公式(1-2)。

   公式(1-2)

针对公式(1-2)，在上述采用最大后验概率估计方法时，引入平均弥散率的稀疏约束。比如：上述稀疏约束函数约束平均弥散率的稀疏性，优选利用L₁范函来约束平均弥散率在稀疏变换域内的稀疏性。

在本发明的一个实施例中，所述稀疏约束函数通过调用以下公式(3)来计算获得。

$R\left(\mathrm{MD}\right)=R\left(\frac{D_1+D_2+D_3}{3}\right)={\left|\right|\mathrm{\Ψ}\left(\frac{D_1+D_2+D_3}{3}\right)\left|\right|}_1$ 公式(3)

即上述公式(1)中的作用于平均弥散率的稀疏约束函数R(MD)表示为上述公式(3)的内容，D₁、D₂、D₃分别为弥散张量矩阵主对角线的元素，运算符Ψ表示某稀疏变换，如小波变换等。‖·‖₁表示取L1范数。本实施例中利用L1范数‖·‖₁用来约束平均弥散率在稀疏变换域内的稀疏性。至此，弥散张量的去噪问题转变为结合公式(3)、利用数值最优化对上述公式(1)的求解问题，而对于上述公式(1)的常用求解方法有非线性共轭梯度法、模拟退火法、Bregman算法、FPC(Fixed-Point Continuation)算法、L1-magic算法、L1-LS算法、牛顿下降法、遗传算法等，在此不做详细说明。本发明不限于只采用L1范数，也可以采用其他范函来实现，具体可参见上述公式(3)进行类似设置。

基于上述实施例，上述弥散参数计算步骤130包括以下步骤：

根据去噪后的弥散张量矩阵，计算弥散参数图对应的弥散参数；这里提到的弥散参数主要是指以下提到的5个基本参数。

对于弥散张量成像技术，还需要根据上述弥散张量矩阵D计算各种弥散参数，如下述说明。

(1)MD，即平均弥散率(mean diffusivity MD)，MD反映分子整体的弥散水平(平均椭球的大小)和弥散阻力的整体情况。MD只表示弥散的大小，而与弥散的方向无关。MD越大，组织内所含自由水分子则越多。计算公式如下述公式(2-4)所示：

MD＝(λ₁+λ₂+λ₃)/3   公式(2-4)

(2)FA，即部分各向异性指数，是水分子各向异性成分占整个弥散张量的比例，它的变化范围从0～1。0代表弥散不受限制，比如脑脊液的FA值接近0；对于非常规则的具有方向性的组织，其FA值大于0，例如大脑白质纤维FA值接近1。FA值的计算公式如下下述公式(2-5)：

$\mathrm{FA}=\sqrt{3[{({\mathrm{\λ}}_1-\mathrm{MD})}^2+{({\mathrm{\λ}}_2-\mathrm{MD})}^2+{({\mathrm{\λ}}_3-\mathrm{MD})}^2]/2({{\mathrm{\λ}}_1}^2+{{\mathrm{\λ}}_2}^2+{{\mathrm{\λ}}_3}^2)}$    公式(2-5)

(3)RA，相对各向异性指数，是弥散张量的各向异性部分与弥散张量各向同性部分的比值，它的变化范围从0(各向同性弥散)到√2(无穷各向异性)。RA的计算公式为下述公式(2-6)：

$\mathrm{RA}=\sqrt{{({\mathrm{\λ}}_1-\mathrm{MD})}^2+{({\mathrm{\λ}}_2-\mathrm{MD})}^2+{({\mathrm{\λ}}_3-\mathrm{MD})}^2/3\left(\mathrm{MD}\right)}$    公式(2-6)

(4)VR，即容积比指数，是椭圆体与球体容积的比值。由于它的变化范围从1(即各向同性弥散)到0，所以，临床上更倾向于应用1/VR。VR的计算公式如下述公式(2-6)：

$\mathrm{RA}=({\mathrm{\λ}}_1\×{\mathrm{\λ}}_2\×{\mathrm{\λ}}_3)/{\mathrm{MD}}^3$    公式(2-6)

上述公式中，λ₁,λ₂,λ₃为弥散张量矩阵D的三个弥散张量矩阵特征值(前述已解释，请参见前述说明)。

(5)弥散张量轨迹(the trace of the diffusion tensor)，其是一个不变参数，Tr(D)＝D₁+D₂+D₃。其中，D₁、D₂、D₃为上述公式(2-3)所示矩阵的主对角线元素。

图1或图2为本发明一个实施例的流程示意图。应该理解的是，虽然图1或图2的流程图中的各个步骤按照箭头的指示依次显示，但是这些步骤并不是必然按照箭头指示的顺序依次执行。除非本文中有明确的说明，这些步骤的执行并没有严格的顺序限制，其可以以其他的顺序执行。而且，图1或图2中的至少一部分步骤可以包括多个子步骤或者多个阶段，这些子步骤或者阶段并不必然是在同一时刻执行完成，而是可以在不同的时刻执行，其执行顺序也不必然是依次进行，而是可以与其他步骤或者其他步骤的子步骤或者阶段结合实施的。以上各个实施例在具体说明中仅只针对相应步骤的实现方式进行了阐述，然后在逻辑不相矛盾的情况下，上述各个实施例是可以相互组合的而形成新的技术方案的，而该新的技术方案依然在本具体实施方式的公开范围内。

基于上述方法，本发明还提供了一种磁共振弥散张量成像的去噪***，如图3所示，该***包括：

图像数据获取模块210，用于获取磁共振弥散加权图像所对应的K空间数据；

去噪模块220，用于基于磁共振弥散加权成像模型和采样噪声的高斯分布性质，利用平均弥散率的稀疏性，采用最大后验概率估计的方法由所述K空间数据获得每一个空间位置所对应的去噪后的弥散张量矩阵；及

弥散参数计算模块230，用于基于去噪后的弥散张量矩阵，获得弥散参数图。

上述图像数据获取模块210主要用于图1所示方法中的步骤110，上述去噪模块220主要用于图1所示方法中的步骤120，上述弥散参数计算模块230主要用于图1所示方法中的步骤130。因此本实施例的***中各个功能模块的具体实现方式参见上述有关步骤110至步骤130的解释说明。

基于上述***，在本发明的一个实施例中，上述去噪模块220包括：

模型构建单元221，用于基于磁共振弥散加权成像模型和采样噪声的高斯分布性质，利用平均弥散率的稀疏性，采用最大后验概率估计的方法构建去噪函数模型；和

矩阵求解单元222，用于利用所述K空间数据，求解上述去噪函数模型，获得每一个空间位置所对应的去噪后的弥散张量矩阵。这里的去噪函数模型参见上述公式(1)。

上述模型构建单元221主要用于图2所示方法中的步骤121，上述矩阵求解单元222主要用于图2所示方法中的步骤122。因此本实施例中的各个功能模块的具体实现方式参见上述有关步骤121至步骤122的解释说明。

通过以上的实施方式的描述，本领域的技术人员可以清楚地了解到上述实施例方法可借助软件加必需的通用硬件平台的方式来实现，当然也可以通过硬件，但很多情况下前者是更佳的实施方式。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个非易失性计算机可读存储介质(如ROM、磁碟、光盘)中，包括若干指令用以使得一台终端设备(可以是手机，计算机，服务器，或者网络设备等)执行本发明各个实施例所述的***结构和方法。

综上所述，本发明运用平均弥散率的固有特性(如稀疏性)，直接通过采集的K空间数据获取去噪后的弥散张量矩阵，相比现有技术，其不需要通过增加采样次数或者减少K空间采样区域的方法，就可以提高弥散张量成像的信噪比，且还不需要通过对重建的弥散加权图像进行去噪处理来获取弥散张量矩阵及相关弥散参数，从而免除了图像去噪中的误差对弥散张量估计的影响，有效抑制了弥散张量中的噪声，提高了弥散张量估计的精度。

以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。

Claims (10)

1.一种磁共振弥散张量成像的去噪方法，其包括：

图像数据获取步骤：获取磁共振弥散加权图像所对应的K空间数据；

去噪步骤：基于磁共振弥散加权成像模型和采样噪声的高斯分布性质，利用平均弥散率的稀疏性，采用最大后验概率估计的方法由所述K空间数据获得每一个空间位置所对应的去噪后的弥散张量矩阵；

弥散参数计算步骤：基于所述去噪后的弥散张量矩阵，获得弥散参数图。

2.根据权利要求1所述的磁共振弥散张量成像的去噪方法，其特征在于，所述去噪步骤包括：

基于磁共振弥散加权成像模型和采样噪声的高斯分布性质，利用平均弥散率的稀疏性，采用最大后验概率估计的方法构建去噪函数模型，所述去噪函数模型参见如下述公式(1)：

       公式(1)

其中，表示弥散张量矩阵的估计值，MD为平均弥散率；R(□)是作用于平均弥散率的稀疏约束函数，λ为相应的正则化参数；d_m为第m个弥散加权图像所对应的K空间数据；F表示傅立叶编码矩阵；表示第m个弥散加权图像，其中，I₀表示无弥散加权的参考图像，为第m个弥散加权图像的相位，b是弥散加权因子，g_m是第m个弥散加权图像所对应的弥散梯度向量g_m＝(g_xm,g_ym,g_zm)^T；

利用所述K空间数据，求解所述公式(1)，获得每一个空间位置所对应的去噪后的弥散张量矩阵。

3.根据权利要求2所述的磁共振弥散张量成像的去噪方法，其特征在于，所述稀疏约束函数约束平均弥散率的稀疏性。

4.根据权利要求2所述的磁共振弥散张量成像的去噪方法，其特征在于，所述稀疏约束函数利用L1范函来约束平均弥散率在稀疏变换域内的稀疏性。

5.根据权利要求2所述的磁共振弥散张量成像的去噪方法，其特征在于，所述稀疏约束函数通过调用以下公式来计算获得：

$R\left(\mathrm{MD}\right)={\left|\right|\mathrm{\Ψ}\left(\frac{D_1+D_2+D_3}{3}\right)\left|\right|}_1$

其中，Ψ表示稀疏变换运算符，||□||₁表示取L₁范数，D₁、D₂、D₃分别为弥散张量矩阵主对角线的元素。

6.一种磁共振弥散张量成像的去噪***，其特征在于，所述***包括：

图像数据获取模块，用于获取磁共振弥散加权图像所对应的K空间数据；

去噪模块，用于基于磁共振弥散加权成像模型和采样噪声的高斯分布性质，利用平均弥散率的稀疏性，采用最大后验概率估计的方法由所述K空间数据获得每一个空间位置所对应的去噪后的弥散张量矩阵；及

弥散参数计算模块，用于基于所述去噪后的弥散张量矩阵，获得弥散参数图。

7.根据权利要求6所述的磁共振弥散张量成像的去噪***，其特征在于，所述去噪模块包括：

模型构建单元，用于基于磁共振弥散加权成像模型和采样噪声的高斯分布性质，利用平均弥散率的稀疏性，采用最大后验概率估计的方法构建去噪函数模型，所述去噪函数模型参见如下述公式(1)：

            公式(1)

其中，表示弥散张量矩阵的估计值，MD为平均弥散率；R(□)是作用于平均弥散率的稀疏约束函数，λ为相应的正则化参数；d_m为第m个弥散加权图像所对应的K空间数据；F表示傅立叶编码矩阵；表示第m个弥散加权图像，其中，I₀表示无弥散加权的参考图像，为第m个弥散加权图像的相位，b是弥散加权因子，g_m是第m个弥散加权图像所对应的弥散梯度向量g_m＝(g_xm,g_ym,g_zm)^T；和

矩阵求解单元，用于利用所述K空间数据，求解所述公式(1)，获得每一个空间位置所对应的去噪后的弥散张量矩阵。

8.根据权利要求7所述的磁共振弥散张量成像的去噪***，其特征在于，所述稀疏约束函数约束平均弥散率的稀疏性。

9.根据权利要求7所述的磁共振弥散张量成像的去噪***，其特征在于，所述稀疏约束函数利用L₁范函来约束平均弥散率在稀疏变换域内的稀疏性。

10.根据权利要求7所述的磁共振弥散张量成像的去噪***，其特征在于，

所述稀疏约束函数通过调用以下公式来计算获得：

$R\left(\mathrm{MD}\right)={\left|\right|\mathrm{\Ψ}\left(\frac{D_1+D_2+D_3}{3}\right)\left|\right|}_1$

其中，Ψ表示稀疏变换运算符，||□||₁表示取L₁范数，D₁、D₂、D₃分别为弥散张量矩阵主对角线的元素。