CN107945129B - 一种mri图像重构方法 - Google Patents

Abstract

本发明公开一种MRI图像重构方法；本发明方法综合考虑换变换域稀疏性与梯度域稀疏性，结合非均匀傅立叶变换建立压缩感知模型，并以交替方向乘子法为基础构建最小优化求解算法，得到MRI重构图像，与基于回溯线搜索的共轭梯度法和快速迭代收缩阈值法相比，本发明方法在各个采样率下均获得了更高质量的重构图像，且算法执行时间更少。

Description

一种MRI图像重构方法

技术领域

本发明属于磁共振成像领域，特别涉及一种MRI图像重构技术。

背景技术

作为二十世纪医学成像最重要的进展之一，磁共振成像(Magnetic ResonanceImaging,MRI)已经成为临床医学检查的重要手段之一，为临床诊断提供了非常有价值的信息。但是，MRI的成像原理就决定了其较长的扫描时间，过慢的成像速度，影响了MRI的进一步推广和应用。

为了提高MRI成像速度，研究者们提出了多种加速采样方法，如并行采样、钥孔成像等。这些加速技术都是在传统的奈奎斯特采样定理框架下进行的，受到采样频率的限制，成像速度提升有限。文献“Michael Lustig,David Donoho,Hohn M Pauly.Sparse MRI:theapplication of compressed sensing for rapid MR imaging[J].Magnetic Resonancein Medicine,2007,58(6):1182–95”第一次将压缩感知(Compressed Sensing,CS)理论引入到MRI成像中(CS-MRI)。这是一个充分利用信号稀疏性或可压缩性的全新信号采集、编解码理论。MRI医学图像本身具有可压缩性，且MRI采集到的数据是空间频率编码(即k空间数据)，这使得压缩感知理论非常适合处理MRI数据。基于压缩感知理论的磁共振成像能够以远低于奈奎斯特采样频率进行数据采样，并从这些欠采样数据中实现图像恢复，从而缩短磁共振扫描时间，提高成像速度。

CS-MRI利用了图像稀疏性这一先验信息，从不完全的采样数据中依据稀疏性原则重构原始图像。重构图像的过程在数学上来讲就是最小l₀范数优化问题。但直接求解该问题数值解极不稳定，是一个NP-Hard问题，且抗噪能力较差。因此，研究者们通过引入l_p范数、全变分(Total Variation,TV)等正则项使该问题变为一个更容易处理的凸优化问题。然后根据不同的正则化模型，采用多种优化算法求最稀疏解，如共轭梯度法(ConjugateGradient,CG)、变量分裂法(Variable Splitting,VS)、迭代硬阈值法(Iterative HardThresholding,IHT)、分裂Bregman迭代法(Splitting Bregman Iteration,SBI)，以及快速迭代收缩阈值法(Fast Iterative Thrinkage/Thresholding Algorithm,FISTA)等。这些算法大多基于迭代过程进行求解，各个算法的执行效率和重构图像质量也有差异。

交替方向乘子法(Alternating Direction Method of Multipliers,ADMM)也是一种迭代优化算法，它综合了变量分裂法、拉格朗日乘子法和交替方向最小化算法的优点，是一种既可以求解约束性的又可以求解非约束的最小优化问题非常有效的方法。ADMM算法在图像处理方面最初主要用于图像降噪，后来开始应用于MRI压缩成像。针对不同的CS-MRI重构模型，算法的求解方法也不相同。如文献“Runge V M,Richter J K,Heverhagen JT.Speed in Clinical Magnetic Resonance[J].Investigative Radiology,2017,52(1):1-17.”针对基于二阶全变分(second-order TV)模型的MR图像重构问题利用变量分裂法将目标问题分为多个子问题，然后采用迭代收缩算法分别求解；而文献“Yang Zhenzhen,YangZhen.Fast linearized alternating direction method of multipliers for theaugmented l1-regularized problem[J].Signal Image&Video Processing,2015,9(7):1601-1612.”针对l₁正则化模型提出一种基于线性展开的快速ADMM方法(作者称为FLADMM)来求解，文献“Chen Shanshan,Du Hongwei,Wu Linna,et al.Compressed sensing MRIvia fast linearized preconditioned alternating direction method ofmultipliers[J].Biomedical Engineering Online,2017,16(1):53.”将这种快速算法用于求解添加全变分平方项(quadratic term)的正则化模型。这些算法要么采用更高阶次的正则化模型，要么增加更多的正则化项，使得算法复杂度增加。

一般的MRI数据采集采用均匀采样，对这种方法采样得到的数据只需要做一个快速傅立叶逆变换即可重构原始图像。但这种方法采集的数据量大，采样时间长，且重建图像会产生较强的混叠伪影。现在常用非均匀采样，如螺旋形采样、径向采样等，以减少采样数据量，提高采样速度。

发明内容

本发明为解决上述技术问题，提出了一种MRI图像重构方法。该方法引入辅助变量，将原始模型中的优化问题分解为多个更容易求解的子优化问题，再交替求解这些问题，以得到原始问题的最稀疏解，从而高概率重构原始MRI图像。

本发明采用的技术方案为：一种MRI图像重构方法，包括：

S1、采用非均匀傅立叶变换对MRI信号进行插值处理；

S2、采用l₁范数和全变分混合正则化项建立压缩感知模型；

其中，m表示MRI图像，S表示MRI信号，

表示经欠采样所得的测量数据，m表示由

恢复得到的图像；α和β为正则化参数，表示平衡两个正则化项在目标函数中所占的比重；Φ为感知矩阵；

S3、基于ADMM重构MRI图像；具体为：由步骤S2得到的优化问题为：

其中，‖·‖₁表示求l₁范数，

表示求l₂范数，

表示图像的空间域的有限差分，D为TV算子，ρ1、ρ₂是添加的惩罚项的惩罚参数，u₁、u₂是拉格朗日乘子的函数。

进一步地，如步骤S3所述，通过交替最小化方法求解，分别优化变量m、z、u₁以及u₂，得到重构MRI图像，具体为：

A1、通过固定变量z，u₁和u₂，得到求解变量m第i+1次的值的子优化问题：

A2、通过固定变量m，u₁和u₂，得到求解变量z第i+1次的值的子优化问题：

A3、根据下式利用得到的mⁱ⁺¹和zⁱ⁺¹更新变量u₁和u₂：

最后，通过循环执行步骤A1、A2和A3，直到满足迭代终止条件，得到重建图像信号m。

更进一步地，所述迭代终止条件为：

|SSIM(mⁱ)-SSIM(m^i-1)|＜ε₃；

其中，SSIM(·)表示求SSIM值运算，ε₃表示迭代终止阈值。

本发明的有益效果：本发明方法综合考虑MRI图像在变换域和梯度域下的稀疏性，并针对目前MRI系统常用的非均匀采样，使用非均匀傅立叶变换实现插值处理，然后在此基础上构建出MRI图像重构模型，最后采用交替方向乘子法求解该模型；与基于回溯线搜索的共轭梯度法和快速迭代收缩阈值法相比，本发明方法在各个采样率下均获得了更高质量的重构图像，且算法执行时间更少。

附图说明

图1为本发明实施例提供的方案流程图；

图2为本发明实施例提供的乐高积木实验中采用CG、FISTA和ADMM三种算法在采样率同为10％、20％、30％下的重构图像；

图3为本发明实施例提供的经CG算法、FISTA算法和本文ADMM方法重构的乐高积木图像的平均SSIM值和算法平均执行时间对比图。

具体实施方式

为便于本领域技术人员理解本发明的技术内容，下面结合附图对本发明内容进一步阐释。

如图1所示为本发明的方案流程图，本发明的技术方案为：一种MRI图像重构方法，具体包括如下步骤：

S1、采用非均匀傅立叶变换对MRI信号进行插值处理；

S2、采用l₁范数和全变分混合正则化项建立压缩感知模型；

MR信号s和MR图像m之间的关系可以用式(1)表示：

s＝∫_Rm·e^-2πik·rdr (1)

可见，MR获取的信号s是m在空间频域的傅立叶域的数据表示，即k空间数据。从而可知，MRI图像重构过程，就是在稀疏域内求最稀疏解的过程，即为l₀最小优化问题。依据压缩感知理论，只要信号s是稀疏的，并且其测量矩阵满足有限等距特性(RestrictedIsometry Property,RIP)，就可以用l₁范数代替求解l₀范数的最小优化问题精确重构原始的MRI图像，即：

其中，||·||表示欧几里得范数，

表示经欠采样所得的测量数据，m表示由

恢复得到的图像，Φ为感知矩阵。式(2)中的模型利用了MRI图像本身具有稀疏性这一先验信息，稀疏变换采用恒等变换。

TV正则化在图像处理中被大量使用，它可以很好的保留图像边缘信息。在压缩感知模型中，即使采用其它正则化方法，通常也将TV作为一个惩罚项包含在其中。因此，本发明综合l₁和TV正则化项的最优化模型为：

式(3)中综合了l₁正则化和TV正则化项，分别对应MRI图像本身的稀疏性和梯度域的稀疏性，称为L₁TV模型。式(3)中α和β为正则化参数，平衡两个正则化项在目标函数中所占的比重。

此外，因为非均匀采样(如径向采样)一般中心采样密度高，而边缘采样密度低，所以需要对采样密度进行补偿。所以，模型中Φ可以描述如下：

Φ＝U{D_CF[F_NUFFT(·)]}

其中，F_NUFFT为非均匀傅立叶变换，D_CF为密度补偿函数(density compensationfunction,DCF)，U为欠采样矩阵。测量矩阵主要包含以下几个过程：

首先，经过非均匀傅立叶变换得到变换域的数据；

然后，利用补偿函数进行采样密度补偿；

最后，经过欠采样得到测量数据。

在本实施例中采用经典的维诺图(Voronoi Diagram)法生成DCF密度补偿函数。

S3、基于ADMM重构MRI图像；

考虑式(3)所示的优化问题，其中全变分即为空间梯度域的l₁范数，即：

其中D为TV算子。引入辅助变量

则优化问题式(3)可转换如下的约束优化问题。

上式(4)的增广拉格朗日表示式如下：

其中，λ₁、λ₂为拉格朗日乘子，ρ₁、ρ₂惩罚参数。令u₁＝λ₁/ρ₁，u₂＝λ₂/ρ₂，式(5)经过简单的代数运算可转换为：

直接求解式(6)比较困难，引入交替最小化方法，并采用Gauss-Seidel思想，分别优化m，z，u₁和u₂变量。每次迭代时只对一个变量进行优化，而固定其它几个变量，就可以把式(6)的最小优化问题转换为更容易处理的几个子问题；具体为：

首先，固定变量z，u₁和u₂为第i次的值，变量m第i+1次的值可以通过如下子优化问题得到：

同理，固定变量m，u₁和u₂为第i次的值，变量z第i+1次的值可以通过如下子优化问题得到：

然后，利用得到的mⁱ⁺¹和zⁱ⁺¹更新变量u₁和u₂：

由于ADMM算法在解决这类最小优化问题时具有全局收敛性，所以循环执行式(7)-(10)，直到满足迭代终止条件，就可以得到重建图像信号m。

在本实施例中式(8)的优化问题可以直接采用熟知的软阈值方法求解，得：

其中，soft(·,·)表示软阈值函数，其定义为：

soft(x,μ)＝sgn(x)max(|x|-μ,0)

其中，sgn(·)为符号函数。

式(7)包含了l₁范数和两个l₂范数的平方项，直接求解比较困难，在本实施例中利用泰勒公式展开，求它的近似解。令：

显然，函数f(m)连续可微，且函数f(m)在点mⁱ处的梯度值为gⁱ＝▽f(mⁱ)，即：

则式(7)可变换为：

其中以γⁱ为元素的对角矩阵γ近似表示海森矩阵(Hesian matrix)▽²f(x)。显然，问题(14)依然可以使用软阈值方法求解，即：

式(3)所表示的优化问题，是要求取目标函数值最小时满足约束条件的变量m的值m^*。因此，最直接的想法就是求得的变量m越接近m^*越好，也就可以设计迭代终止条件为：

‖m-m^*‖₂＜ε₁ (16)

其中，ε₁＞0为允许的误差。

由于m^*是未知的，式(16)实际上是不可能实现的。所以在实际实现的时候，常常用目标函数或者变量的相对变化值作为结束条件，即：

其中，xⁱ为第i次的目标函数值或者变量值。

由于本发明所述算法最后的变量是重构的MRI图像数据，而最后判断图像质量的时候，本发明采用图像的结构相似度(Structural Similarity Image Metric,SSIM)作为衡量重构图像质量好坏的指标。因此，本实施例以图像的SSIM值的变化作为式(3)所表示的优化问题的迭代终止的条件，即：

|SSIM(mⁱ)-SSIM(m^i-1)|＜ε₃ (18)

其中，SSIM(·)表示求SSIM值运算，ε₃表示迭代终止阈值，是一个大于0的比较小的值，两次迭代后的SSIM值相差小于该值就结束迭代。

为了验证前述算法在非均匀MRI图像重构应用中的有效性，以一个乐高积木的MRI径向采样数据为实验对象，将本发明方法和基于回溯线搜索的共轭梯度法和当前流行的FISTA法进行比较，比较三种算法的重构图像质量和算法执行时间。

图2所示为乐高积木实验中采用CG、FISTA和ADMM三种算法在采样率同为10％、20％、30％下的重构图像。图中第一行至第三行对应的采样率分别为10％、20％和30％，第一列为经CG算法重构所得图像，第二列为经FISTA算法重构所得图像，第三列为经ADMM算法重构所得图像。图像下方所标示的SSIM值是该采样率下进行10次随机采样，然后使用相应重构算法分别进行MRI图像重构所得图像的SSIM平均值。

从图中可以看出，采样率越高，采样数据从原始图像中获取的信息也就越多，重构图像的SSIM值也越大，图像质量越高。同时，针对相同的采样率，采用本文ADMM方法重构的图像比采用CG算法和FISTA算法得到的图像噪声更少，质量更高。图中SSIM数据显示，即使采样率为20％，经过本文算法所重构图像的SSIM值也达到了0.8以上，肉眼上看和MRI全采样重构图比较接近，当采样率达到30％时就更接近了。

图3显示了采样率在5％到30％之间，以1％为步进的各个采样率下分别经CG算法、FISTA算法和本文ADMM方法重构的乐高积木图像的平均SSIM值和算法平均执行时间。图中粗实线示的是平均SSIM值，细虚线表示的是算法平均执行时间。从图中可以直观的看出，在各个采样率下，采用非均匀MRI压缩图像重构的ADMM算法恢复得到的图像质量明显更高。同时也可以看出，CG算法所消耗的执行时间最长，本文算法在整个采样范围内执行时间几乎不变，比FISTA算法的执行时间略少，且采样率越低，时间差越大。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的权利要求范围之内。

Claims (3)

1.一种磁共振成像图像重构方法，包括：

S1、采用非均匀傅立叶变换对磁共振成像信号进行插值处理；

S2、采用l₁范数和全变分混合正则化项建立压缩感知模型；

其中，‖m‖_TV表示全变分正则化项，m表示磁共振成像图像，S表示磁共振成像信号，

表示经欠采样所得的测量数据，m表示由

恢复得到的图像；α和β为正则化参数，表示平衡两个正则化项在目标函数中所占的比重；Φ为感知矩阵；

S3、基于交替方向乘子法重构磁共振成像图像；具体为：由步骤S2得到的优化问题为：

其中，‖·‖₁表示求l₁范数，

表示求l₂范数，z＝Dm表示图像的空间域的有限差分，D为全变分算子，ρ₁、ρ₂是添加的惩罚项的惩罚参数，u₁、u₂是拉格朗日乘子的函数；

通过交替最小化方法求解，分别优化变量m、z、u₁以及u₂，得到重构磁共振成像图像。

2.根据权利要求1所述的磁共振成像图像重构方法，其特征在于，所述步骤S3所述分别优化m，z，u₁和u₂变量，通过交替最小化方法求解得到重构磁共振成像图像，具体为：

A1、通过固定变量z，u₁和u₂，得到求解变量m第i+1次的值的子优化问题：

通过固定变量m，u₁和u₂，得到求解变量z第i+1次的值的子优化问题：

A2、根据下式利用得到的mⁱ⁺¹和zⁱ⁺¹更新变量u₁和u₂：

最后，通过循环执行步骤A1、A2，直到满足迭代终止条件，得到重建图像信号m。

3.根据权利要求2所述的磁共振成像图像重构方法，其特征在于，

所述迭代终止条件为：

|SSIM(mⁱ)-SSIM(m^i-1)|＜ε₃；

其中，SSIM(·)表示求结构相似度值运算，ε₃表示迭代终止阈值。

Images