CN104504657B - 磁共振弥散张量去噪方法和装置 - Google Patents

Abstract

本发明提供了一种磁共振弥散张量去噪方法和装置，其中，该方法包括：采集K空间数据；根据采集的K空间数据，求取弥散张量的最大似然估计值；利用弥散张量和弥散参数的稀疏性，将求取的最大似然估计值作为初始值，计算弥散张量的最大后验概率估计值；根据求取的所述最大后验概率估计值，计算弥散系数。本发明解决了现有技术中无法在不增加扫描时间和不影响空间分辨率的同时，实现弥散张量的高精度去噪的技术问题，达到了有效抑制弥散张量中的噪声，提高弥散张量估计精度的技术效果。

Description

磁共振弥散张量去噪方法和装置

技术领域

本发明涉及磁共振成像技术领域，特别涉及一种磁共振弥散张量去噪方法和装置。

背景技术

磁共振弥散张量成像(Diffusion Tensor Imaging，简称为DTI)是利用组织中水分子自由热运动的各向异性原理，探测组织的微观结构，以达到研究人体功能的目的。

目前，DTI是唯一一种可以活体显示脑白质纤维束的无创性成像方法。然而，DTI的信噪比比较低，这在很大程度上限制了DTI在临床上的广泛应用。

磁共振弥散张量成像的物理机制可以表示为：

d＝Fρ+n (公式1)

其中，d表示在磁共振仪上所采集的信号，F表示傅立叶编码矩阵，ρ表示弥散加权图像，n通常假定为复高斯白噪声。

在弥散张量成像中，第m个弥散加权图像ρ_m可以表示为：

其中，I₀∈R^N×1表示无弥散加权的参考图像，表示第m个弥散加权图像的相位，b表示弥散加权因子(常量)，g_m表示第m个弥散加权图像所对应的弥散梯度向量g_m＝(g_xm,g_ym,g_zm)^T，D表示弥散张量，具体的D可以表示为：

其中，每个D_i∈R^N×1，N表示磁共振图像像素点的个数。

具体的，弥散张量成像的过程，主要包括以下步骤：

步骤1：由采集信号d重建图像ρ；

步骤2：根据模型张量(即上述公式2)估计每一空间位置所对应的张量矩阵D；

步骤3：根据下述公式4和公式5，由张量矩阵D计算得到各种弥散参数，例如：MD、FA等：

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

其中，λ₁,λ₂,λ₃为矩阵D的特征值。

为了提高弥散张量成像的信噪比，目前较为直接的方法是通过多次采样取平均或者是减小K空间的采样区域。虽然这些方法在实际应用中可以提高信噪比，但是这些方法会增加扫描时间或者影响空间分辨率。

目前，还有一种常用的提高弥散张量成像质量的方法。该方法在获取K空间的扫描数据后，首先重建出弥散加权图像。然后运用信号处理的方法对图像进行去噪。最后，通过去噪后的图像计算弥散张量及各种弥散参数。即，在处理的过程中需要先重建图像，然后对图像进行去噪。可是在图像去噪过程中产生的***误差往往会进一步传递到后续的张量计算中，从而影响弥散张量的精度。

针对如何在不增加扫描时间和不影响空间分辨率的同时，实现弥散张量的高精度去噪的技术问题，目前尚未提出有效的解决方案。

发明内容

本发明实施例提供了一种磁共振弥散张量去噪方法，以解决现有技术中无法在不增加扫描时间和不影响空间分辨率的同时，实现弥散张量的高精度去噪的技术问题，该方法包括：

采集K空间数据；

根据采集的K空间数据，求取弥散张量的最大似然估计值；

利用弥散张量和弥散参数的稀疏性，以求取的最大似然估计值作为初始值，计算弥散张量的最大后验概率估计值；

根据求取的最大后验概率估计值，计算弥散系数。

在一个实施例中，根据采集的K空间数据，求取弥散张量的最大似然估计值，包括：

按照以下公式计算弥散张量的最大似然估计值：

其中，表示弥散张量的最大似然估计值，d_m表示第m个弥散加权图像所对应的K空间数据，F表示傅立叶编码矩阵，I₀表示无弥散加权的参考图像，表示第m个弥散加权图像的相位，b表示弥散加权因子，g_m＝(g_xm,g_ym,g_zm)^T表示第m个弥散梯度向量。

在一个实施例中，利用弥散张量和弥散参数的稀疏性，将求取的最大似然估计值作为初始值，计算弥散张量的最大后验概率估计值，包括：

按照以下公式计算弥散张量的最大后验概率估计值：

其中，R(·)表示作用于弥散张量D的惩罚函数，G(·)表示作用于弥散参数FA的惩罚函数，γ₁、γ₂表示正则化参数、D₁,D₂,D₃,D₄,D₅,D₆代表不同方向的6个弥散张量。

在一个实施例中，R(·)约束每个弥散张量的稀疏性，或者R(·)约束弥散张量的联合稀疏性。

在一个实施例中，在R(·)约束每个弥散张量稀疏性的情况下，

在R(·)约束弥散张量的联合稀疏性的情况下：

其中，Ψ为运算符，表示稀疏变换，ψ_i表示Ψ的第i列，||·||₁表示求L1范数，||·||₂表示求L2范数。

本发明实施例还提供了一种磁共振弥散张量去噪装置，以解决现有技术中无法在不增加扫描时间和不影响空间分辨率的同时，实现弥散张量的高精度去噪的技术问题，该装置包括：

K空间数据采集模块，用于采集K空间数据；

最大似然估计求取模块，用于根据采集的K空间数据，求取弥散张量的最大似然估计值；

最大后验概率估计模块，用于利用弥散张量和弥散参数的稀疏性，将求取的最大似然估计值作为初始值，计算弥散张量的最大后验概率估计值；

弥散系数计算模块，用于根据求取的所述最大后验概率估计值，计算弥散系数。

在一个实施例中，所述最大似然估计求取模块，具体用于按照以下公式计算弥散张量的最大似然估计值：

其中，表示弥散张量的最大似然估计值，d_m表示第m个弥散加权图像所对应的K空间数据，F表示傅立叶编码矩阵，I₀表示无弥散加权的参考图像，表示第m个弥散加权图像的相位，b表示弥散加权因子，g_m＝(g_xm,g_ym,g_zm)^T表示第m个梯度向量。

在一个实施例中，所述最大后验概率估计模块，具体用于按照以下公式计算弥散张量的最大后验概率估计值：

其中，R(·)表示作用于弥散张量D的惩罚函数，G(·)表示作用于弥散参数FA的惩罚函数，γ₁、γ₂表示正则化参数、D₁,D₂,D₃,D₄,D₅,D₆代表不同方向的6个弥散张量。

在一个实施例中，R(·)约束每个弥散张量的稀疏性，或者R(·)约束弥散张量的联合稀疏性。

在一个实施例中，在R(·)约束每个弥散张量的稀疏性的情况下，

在R(·)约束弥散张量的联合稀疏性的情况下：

其中，Ψ为运算符，表示稀疏变换，ψ_i表示Ψ的第i列，||·||₁表示求L1范数，||·||₂表示求L2范数。

在本发明实施例中，在采集到K空间数据后，直接根据K空间数据确定弥散张量的最大似然估计值，然后依据求取的最大似然估计值，结合弥散张量、弥散参数的稀疏性计算弥散张量的最大后验概率估计值，从而确定去噪后的弥散张量。这有别于直接对弥散加权图像进行去噪的方法，避免了因为图像去噪而引入的***误差，解决了现有技术中无法在不增加扫描时间和不影响空间分辨率的同时，实现弥散张量的高精度去噪的技术问题，有效地抑制了弥散张量中的噪声。进一步的，因为通过弥散张量、弥散参数的稀疏性来进行去噪，从而更有效地提高了弥散张量的估计精度。

附图说明

此处所说明的附图用来提供对本发明的进一步理解，构成本申请的一部分，并不构成对本发明的限定。在附图中：

图1是根据本发明实施例的磁共振弥散张量去噪方法的方法流程图；

图2是根据本发明实施例的磁共振弥散张量去噪装置的结构框图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，下面结合实施方式和附图，对本发明做进一步详细说明。在此，本发明的示意性实施方式及其说明用于解释本发明，但并不作为对本发明的限定。

发明人考虑到，在现有的弥散张量成像过程中，需要先得到弥散加权图像，再对弥散加权图像进行去噪。然而这样会引入***误差，从而导致最终的弥散张量估计出现误差。因此，在本例中提出了一种磁共振弥散张量去噪方法。该方法直接通过获取的K空间数据对弥散张量进行去噪，无需求取弥散加权图像，也不需要对弥散加权图像进行去噪。如图1所示，该方法包括以下步骤：

步骤101：采集K空间数据；

步骤102：根据采集的K空间数据，求取弥散张量的最大似然估计值；

步骤103：利用弥散张量和弥散参数的稀疏性，以求取的最大似然估计值作为初始值，计算弥散张量的最大后验概率估计值；

步骤104：根据求取的最大后验概率估计值，计算弥散系数，其中，弥散系数可以包括：MD、FA。

在上述实施例中，在采集到K空间数据后，直接根据K空间数据确定弥散张量的最大似然估计值。然后以求取的最大似然估计值作为初值，结合弥散张量、弥散参数的稀疏性计算弥散张量的最大后验概率估计值，进而确定去噪后的弥散张量。这不同于对弥散加权图像进行去噪的方法，避免了因为图像去噪而引入的***误差，解决了现有技术中无法在不增加扫描时间和不影响空间分辨率的同时，实现弥散张量的高精度去噪的技术问题，有效抑制了弥散张量中的噪声。进一步的，因为通过弥散张量、弥散参数的稀疏性来进行去噪，从而更有效地提高了弥散张量的估计精度。

具体的，假设采样噪声服从高斯分布，结合上述公式1和公式2，弥散张量的最大似然估计值可以写作：

其中，表示弥散张量的最大似然估计值，d_m表示第m个弥散加权图像所对应的K空间数据，F表示傅立叶编码矩阵，I₀表示无弥散加权的参考图像，表示第m个弥散加权图像的相位，b表示弥散加权因子，g_m＝(g_xm,g_ym,g_zm)^T表示第m个弥散加权图像所对应的弥散梯度向量。

可以按照以下公式计算弥散张量的最大后验概率估计值：

其中：

因此，上述公式可以表示为：

其中，R(·)表示作用于弥散张量D的惩罚函数，G(·)表示作用于弥散参数的惩罚函数，γ₁、γ₂表示正则化参数、D₁,D₂,D₃,D₄,D₅,D₆代表不同方向的6个弥散张量。假定D和FA都是稀疏的，那么可以令：

惩罚函数G(·)＝||Ψ(·)||₁分别来约束6个弥散张量：(D₁,D₂,D₃,D₄,D₅,D₆)和弥散参数FA的稀疏性，

或者，令：

来约束弥散张量的联合稀疏性；

其中，运算符Ψ表示某种稀疏变换(例如：小波变换等)，ψ_i表示Ψ的第i列，||·||₁表示求L1范数，用于约束弥散张量在稀疏变换域内的稀疏性，||·||₂表示求L2范数。

这样，弥散张量的去噪问题等价于上述的最优化求解问题。具体的可以采用非线性共轭梯度法进行求解。

即，在本例中，利用弥散张量模型和弥散张量、弥散参数FA的固有特性(如稀疏性)，直接通过获取的K空间数据对弥散张量进行去噪，该方法略过了传统的图像去噪的步骤，避免了图像去噪的误差对张量估计的影响，可以更有效地抑制弥散张量中的噪声，提高弥散张量的估计精度。

基于同一发明构思，本发明实施例中还提供了一种磁共振弥散张量去噪装置，如下面的实施例所述。由于磁共振弥散张量去噪装置解决问题的原理与磁共振弥散张量去噪方法相似，因此磁共振弥散张量去噪装置的实施可以参见磁共振弥散张量去噪方法的实施，重复之处不再赘述。以下所使用的，术语“单元”或者“模块”可以实现预定功能的软件和/或硬件的组合。尽管以下实施例所描述的装置较佳地以软件来实现，但是硬件，或者软件和硬件的组合的实现也是可能并被构想的。图2是本发明实施例的磁共振弥散张量去噪装置的一种结构框图，如图2所示，包括：K空间数据采集模块201、最大似然估计求取模块202、最大后验概率估计模块203和弥散系数计算模块204，下面对该结构进行说明。

K空间数据采集模块201，用于采集K空间数据；

最大似然估计求取模块202，用于根据采集的K空间数据，求取弥散张量的最大似然估计值；

最大后验概率估计模块203，用于利用弥散张量和弥散参数的稀疏性，将求取的最大似然估计值作为初始值，计算弥散张量的最大后验概率估计值；

弥散系数计算模块204，用于根据求取的所述最大后验概率估计值，计算弥散系数。

在一个实施例中，最大似然估计求取模块202，具体用于按照以下公式计算弥散张量的最大似然估计值：

其中，表示弥散张量的最大似然估计值，d_m表示第m个弥散加权图像所对应的K空间数据，F表示傅立叶编码矩阵，I₀表示无弥散加权的参考图像，表示第m个弥散加权图像的相位，b表示弥散加权因子，g_m＝(g_xm,g_ym,g_zm)^T表示第m个梯度向量。

在一个实施例中，最大后验概率估计模块203，具体用于按照以下公式计算弥散张量的最大后验概率估计值：

其中，R(·)表示作用于弥散张量D的惩罚函数，G(·)表示作用于弥散参数FA的惩罚函数，γ₁、γ₂表示正则化参数、D₁,D₂,D₃,D₄,D₅,D₆代表不同方向的6个弥散张量。

在一个实施例中，R(·)约束每个弥散张量的稀疏性，或者R(·)约束弥散张量的联合稀疏性。

在一个实施例中，在R(·)约束每个弥散张量的稀疏性的情况下，

在R(·)约束弥散张量的联合稀疏性的情况下：

其中，Ψ为运算符，表示稀疏变换，ψ_i表示Ψ的第i列，||·||₁表示求L1范数，||·||₂表示求L2范数。

在另外一个实施例中，还提供了一种软件，该软件用于执行上述实施例及优选实施方式中描述的技术方案。

在另外一个实施例中，还提供了一种存储介质，该存储介质中存储有上述软件，该存储介质包括但不限于：光盘、软盘、硬盘、可擦写存储器等。

从以上的描述中，可以看出，本发明实施例实现了如下技术效果：在采集到K空间数据后，直接根据K空间数据确定弥散张量的最大似然估计值，然后依据求取的最大似然估计值，结合弥散张量、弥散参数的稀疏性计算弥散张量的最大后验概率估计值，从而确定去噪后的弥散张量。这有别于直接对弥散加权图像进行去噪的方法，避免了因为图像去噪而引入的***误差，解决了现有技术中无法在不增加扫描时间和不影响空间分辨率的同时，实现弥散张量的高精度去噪的技术问题，有效抑制了弥散张量中的噪声。进一步的，因为通过弥散张量、弥散参数的稀疏性来进行去噪，从而更有效地提高了弥散张量的估计精度。

显然，本领域的技术人员应该明白，上述的本发明实施例的各模块或各步骤可以用通用的计算装置来实现，它们可以集中在单个的计算装置上，或者分布在多个计算装置所组成的网络上，可选地，它们可以用计算装置可执行的程序代码来实现，从而，可以将它们存储在存储装置中由计算装置来执行，并且在某些情况下，可以以不同于此处的顺序执行所示出或描述的步骤，或者将它们分别制作成各个集成电路模块，或者将它们中的多个模块或步骤制作成单个集成电路模块来实现。这样，本发明实施例不限制于任何特定的硬件和软件结合。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明实施例可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种磁共振弥散张量去噪方法，其特征在于，包括：

采集K空间数据；

根据采集的K空间数据，求取弥散张量的最大似然估计值；

利用弥散张量和弥散参数的稀疏性，将求取的最大似然估计值作为初始值，计算弥散张量的最大后验概率估计值；

根据求取的所述最大后验概率估计值，计算弥散系数；

其中，利用弥散张量和弥散参数的稀疏性，将求取的最大似然估计值作为初始值，计算弥散张量的最大后验概率估计值，包括：

按照以下公式计算弥散张量的最大后验概率估计值：

其中，R(·)表示作用于弥散张量D的惩罚函数，G(·)表示作用于弥散参数FA的惩罚函数，γ₁、γ₂表示正则化参数、D₁,D₂,D₃,D₄,D₅,D₆代表不同方向的6个弥散张量。

2.如权利要求1所述的方法，其特征在于，根据采集的K空间数据，求取弥散张量的最大似然估计值，包括：

按照以下公式计算弥散张量的最大似然估计值：

其中，表示弥散张量的最大似然估计值，d_m表示第m个弥散加权图像所对应的K空间数据，F表示傅立叶编码矩阵，I₀表示无弥散加权的参考图像，表示第m个弥散加权图像的相位，b表示弥散加权因子，g_m＝(g_xm,g_ym,g_zm)^T表示第m个梯度向量。

3.如权利要求2所述的方法，其特征在于，R(·)约束每个弥散张量的稀疏性，或者R(·)约束弥散张量的联合稀疏性。

4.如权利要求3所述的方法，其特征在于：

在R(·)约束每个弥散张量的稀疏性的情况下，

在R(·)约束弥散张量的联合稀疏性的情况下：

其中，Ψ为运算符，表示稀疏变换，ψ_i表示Ψ的第i列，||·||₁表示求L1范数，||·||₂表示求L2范数。

5.一种磁共振弥散张量去噪装置，其特征在于，包括：

K空间数据采集模块，用于采集K空间数据；

最大似然估计求取模块，用于根据采集的K空间数据，求取弥散张量的最大似然估计值；

最大后验概率估计模块，用于利用弥散张量和弥散参数的稀疏性，将求取的最大似然估计值作为初始值，计算弥散张量的最大后验概率估计值；

弥散系数计算模块，用于根据求取的所述最大后验概率估计值，计算弥散系数；

其中，所述最大后验概率估计模块，具体用于按照以下公式计算弥散张量的最大后验概率估计值：

其中，R(·)表示作用于弥散张量D的惩罚函数，G(·)表示作用于弥散参数FA的惩罚函数，γ₁、γ₂表示正则化参数、D₁,D₂,D₃,D₄,D₅,D₆代表不同方向的6个弥散张量。

6.如权利要求5所述的装置，其特征在于，所述最大似然估计求取模块，具体用于按照以下公式计算弥散张量的最大似然估计值：

其中，表示弥散张量的最大似然估计值，d_m表示第m个弥散加权图像所对应的K空间数据，F表示傅立叶编码矩阵，I₀表示无弥散加权的参考图像，表示第m个弥散加权图像的相位，b表示弥散加权因子，g_m＝(g_xm,g_ym,g_zm)^T表示第m个梯度向量。

7.如权利要求6所述的装置，其特征在于，R(·)约束每个弥散张量的稀疏性，或者R(·)约束弥散张量的联合稀疏性。

8.如权利要求7所述的装置，其特征在于：

在R(·)约束每个弥散张量的稀疏性的情况下，

在R(·)约束弥散张量的联合稀疏性的情况下：

其中，Ψ为运算符，表示稀疏变换，ψ_i表示Ψ的第i列，||·||₁表示求L1范数，||·||₂表示求L2范数。