CN112001981A - 基于广义涅斯捷罗夫加速共轭梯度算法的压缩采样mr影像重建方法 - Google Patents

Abstract

一种基于广义涅斯捷罗夫加速共轭梯度算法的压缩采样MR影像重建方法，其特征是包括如下步骤：步骤(1)：重构原始磁共振图像：在基于压缩感知的磁共振成像中，欠采样磁共振图像重构的基本模型为：min||X||₀ s.t.Y＝ΦX，其中X是重建的MR图像，是Y在k空间的欠采样MR图像数据，Φ是观察矩阵；Φ表达为：Φ＝μF，其中F是一个二维傅里叶变换，μ代表一个欠采样模式；在满足约束等距性的条件下，可以实现图像的精确重建；步骤(2)：通过范数正则化找到最稀疏的解：步骤(3)：引入TV模型和改进涅斯捷罗夫加速算法，提出了广义涅斯捷罗夫加速的全变分共轭梯度算法进行CS‑MRI重建。本发明进一步提高收敛速度，获得更好的MR图像重建效果。

Description

基于广义涅斯捷罗夫加速共轭梯度算法的压缩采样MR影像重 建方法

技术领域

本发明涉及计算机信息处理技术，具体是一种基于广义涅斯捷罗夫加速共轭梯度算法的压缩采样MR影像重建方法。

背景技术

磁共振成像(MRI)技术属无创伤、无射线检查，具有对比分辨率高、多方位、多参数采集及功能成像等优势，它提供的软组织对比度的能力是大多数其他成像方式所无法比拟的，可以使用非侵入式检查方式检查和量化组织的代谢和生理特征，从而提供关于病理的有价值的信息，目前已经成为一种重要的临床检查方式。然而，MRI技术成像速度相对较慢，数据采集时间过长；MR图像容易受到与运动相关伪影的影响，相对较低的时间分辨率限制了会随呼吸移动的身体部位(如腹部和心脏)的成像，MRI的缺点已成为其在实际应用中的主要障碍。

近年来，Donoho和Candes等人在前人的研究基础上系统的提出了压缩感知(Compressed Sensing，CS)理论，根据该理论：如果待重建图像在某个变换域是稀疏的，就可以通过低于或远低于奈奎斯特采样标准的方式对其进行数据采样并高概率重构该图像。

磁共振图像由于其自身在特定变换域中的稀疏性可以应用压缩感知理论，在MRI应用过程中，根据压缩感知理论,图像重建所需的数据来自于对k空间数据进行欠采样,选取合适的欠采样方式可以用少量的数据准确恢复出所需磁共振图像。

MRI重构方法主要有有最基本梯度下降系列算法，以及衍生而来的阈值收缩系列算法和分裂布雷格曼迭代。其中阈值收缩系列算法采用收缩算子来求解L1最优化问题提高了迭代效率，而收缩算子的本质也是梯度下降。分裂布雷格曼迭代则讲多约束问题分裂为L1和L2最优化问题，分别可以用收缩算子和高斯-塞舌尔迭代来交替求解，可有效求解多约束重构问题。然而，作为不少算法的基础，梯度下降法仍然被很多学者研究和应用，尤其在信号去噪、信号重构、机器学习、深度学习等领域。因此对梯度下降算法的持续改进很有必要。

为了在早期迭代中快速收敛以及进一步提高图像重建的分辨率，本文，以共轭梯度算法为基础，引入TV模型和改进涅斯捷罗夫加速算法，提出了广义涅斯捷罗夫加速的全变分共轭梯度算法。

GNACG_TV使用涅斯捷罗夫加速来优化梯度下降法方法，这种方法通常被称为NAG，但是这种形式的加速因子是不受控制的，存在以下问题：收敛曲线不一定是最优的。加速因子不受控制，只与迭代次数有关。为了解决这个问题，提出了广义涅斯捷罗夫加速度。通过引入目标函数的Frobenius范数作为参数，使其不仅与迭代次数有关，而且与迭代过程有关，从而保证了迭代过程的收敛性。

发明内容

本发明提供一种基于广义涅斯捷罗夫加速共轭梯度算法的压缩采样MR影像重建方法。进一步提高收敛速度，获得更好的MR图像重建效果。

为了加快磁共振成像(MRI)的扫描速度，提高磁共振(MR)图像重建的质量，提出了一种基于压缩感知的快速MRI技术。涅斯捷罗夫的加速梯度下降法(NAG)算法使用涅斯捷罗夫加速来优化梯度下降法(GD)方法。然而，这种形式的加速因子使用固定的迭代曲线更新，不能适应不同的迭代过程。提出了广义涅斯捷罗夫加速度概念。结合总变分模型，提出了一种基于总变分的广义涅斯捷罗夫加速共轭梯度算法。将加速度因子进行了广义扩展，将目标函数的Frobenius范数作为参数引入，使加速度因子不仅与迭代次数有关，而且与迭代过程有关，保证了迭代过程的收敛性。在不同采样率下对腹部、头部和脚踝的三幅图像进行了实验，结果表明，GNACG_TV算法在去噪性能、均方误差(MSE)、峰值信噪比(PSNR)和结构相似度(SSIM)性能上均优于共轭梯度算法(CG)、基于总变分的共轭梯度算法(CG_TV)和基于总变分的涅斯捷罗夫加速共轭梯度算法(NACG_TV)。通过与定性和定量分析结果的比较，表明该方法比其他3种方法能更好地重建欠采样图像。GNACG_TV可以进一步提高基于涅斯捷罗夫加速算法的收敛速度，获得更好的重建性能。

本发明的技术方案如下：

一种基于广义涅斯捷罗夫加速共轭梯度算法的压缩采样MR影像重建方法，其特征是包括如下步骤：

在基于磁共振的压缩感知成像技术中，欠采样磁共振图像重建的基本模型是：

min||X||₀ s.t.Y＝ΦX (1)

其中X是重建的MR图像，Y是在k空间的欠采样MR图像数据，Φ是观察矩阵；

Φ被表达为：

Φ＝μF (2)

其中F是一个二维傅里叶变换，μ代表一个欠采样模式；在满足约束等距性的条件下，可以实现图像的精确重建；

上述方程是一个欠定问题，通过范数正则化得到最优解。L0-norm范数优化是一个NP-hard问题；常用的方法是将其转化为最优凸逼近L1-norm范数优化，其中X表示为：

X＝Ψ^Tθ (3)

其中Ψ＝[ψ₁,ψ₂,...,ψ_N]∈R^N×N和θ是投影系数矩阵,

其中ΦΨ^T是编码矩阵，且

A＝ΦΨ^T (4)

考虑到成像过程中噪声的影响，将上述模型转化为:

min||θ||₁ s.t.||Aθ-Y||₂＜ε (5)

其中，阈值参数ε控制重建图像对测量数据的保真度，阈值参数ε的值通常设置为低于期望噪声水平；

将L1-norm范数转换为更简单的L2-norm范数，MR图像可以通过优化下列问题重建：

其中λ是正则化参数，用于平衡保真项和正则项的比例；通过求解方程(6)可以得到最优解；并将最优解引入到方程(3)中重建磁共振的原始影像；

基于全变分算法的广义涅斯捷罗夫加速共轭梯度：

在CS-MRI图像重建模型中加入了TV正则项，并在此基础上建立了图像重建模型，方程式(6)可更改为：

引入整体变分的正则化可以更好地保持图像的边缘信息，利用正则化参数λ₁,λ₂调整权重，选择合适的正则化参数值可以达到理想的重建效果；

采用一种简单的CS-MR图像重建模型，即在方程式(7)中设置稀疏正则项的权重系数λ₁至0，方程式(7)可简化为:

TV算法的梯度公式如下:

被表示为：

其中，x_i,j是图像中像素的值，v_i,j就是TV算法在x_i,j位置上的梯度；

CS-MR图像重建可以通过在现有共轭梯度法的基础上优化上述模型并提出广义涅斯捷罗夫加速算法的概念来解决；

在CG算法中，第k+1次迭代的公式是:

x_k+1＝x_k+a_kd_k (11)

其中x_k是当前迭代点,a_k步长，也称为学习率，d_k是搜索方向；

在每次迭代过程中，CG方法每次迭代的迭代点是在前一次迭代过程中更新的函数位置。根据方程式(9)，步长和搜索方向的确定是为了找到一个合适的功能点位置，因此确定下一次迭代的合适点对于共轭梯度法来说也是一个重要步骤；

结合涅斯捷罗夫加速算法的概念，利用当前迭代点和前一个迭代点的位置，修改共轭梯度法以更新方程式(9)为：

x_k+1＝z_k+1+a_kd_k (12)

z_k+1为：

t_k+1为：

n_k为：

||f(θ)_k||_F和||f(θ)_k-1||_F为目标函数迭代k和k-1次的Frobenius范数；h是一个可调整的权重，k的取值范围为k＞3；

随着迭代次数的增加，n_k≈2，

与来自NAG算法的t_k+1一致；

用GNACG_TV算法重建MR图像：

压缩感知理论可以从变换域中的稀疏图像重建图像；本发明通过螺旋桨采样轨迹获取观测数据，利用GNACG_TV重建图像；采样轨迹如图1所示。GNACG_TV的完整算法流程如下：

基于GNACG_TV算法的CS-MR图像重建算法：

输入：

Y为K空间中欠采样的MR图像数据

ψ为多贝西小波基础

初始化:

k＝1,n₁＝0.2,t₁＝1,z₁＝0,λ＝0.01,h＝5

当k＜max_iter时,进入循环

k＞3

通过方程式(15)计算n_k；

通过方程式(14)计算t_k+1；

通过方程式(13)计算z_k+1；

通过方程式(12)计算x_k+1并更新搜索方向；

k＝k+1；

结束，得到输出：重建的MR图像X＝Ψ^Tθ。

本发明的有益效果如下：

经过定量和定性分析表明，与CG、CG_TV和NACG_TV算法相比，本发明的一种基于广义涅斯捷罗夫加速共轭梯度算法的压缩采样MR影像重建方法具有较好的重建效率和稳定性，在均方误差、峰值信噪比和结构相似性有较高的性能，在去噪性能方面有较好的鲁棒性。结果表明，该方法比其他3种方法能更好地重建欠采样MR图像，并能进一步提高基于涅斯捷罗夫加速度的收敛速度，获得更好的重建效果。

附图说明

图1为螺旋桨取样轨迹示意图。

图2为本发明的测试用MR影响示意图。

图3为在20％的采样率采用不同算法进行重建，局部细节和误差图像的的腹部MR对比图像。

图4为对于20％采样率的头部图像，采用不同的算法重建局部细节和误差的对比图像。

图5为对于20％采样率的踝部图像，采用不同的算法重建局部细节和误差的对比图像。

图6为在20％的采样率下，用MSE、PSNR和SSIM算法重建的腹部MR图像。

图7为在20％的采样率下，用MSE、PSNR和SSIM算法重建的头部MR图像。

图8为在20％的采样率下，用MSE、PSNR和SSIM算法重建的踝部MR图像。

图9为在40％的采样率采用不同算法进行重建，局部细节和误差图像的的腹部MR对比图像。

图10为不同的采样率下，用MSE、PSNR和SSIM算法重建的腹部MR图像。

图11为在不同的算法和噪声水平下样本比率为20％的腹部MR图像去噪效果对比图。

具体实施方式

下面结合附图对本发明做进一步描述。

如图1至图11，在基于磁共振的压缩感知成像技术中，欠采样磁共振图像重建的基本模型是：

min||X||₀ s.t.Y＝ΦX (1)其中X是重建的MR图像，Y是在k空间的欠采样MR图像数据，Φ是观察矩阵；

Φ被表达为：

Φ＝μF (2)

其中F是一个二维傅里叶变换，μ代表一个欠采样模式；在满足约束等距性的条件下，可以实现图像的精确重建；

上述方程是一个欠定问题，通过范数正则化得到最优解。L0-norm范数优化是一个NP-hard问题；常用的方法是将其转化为最优凸逼近L1-norm范数优化，其中X表示为：

X＝Ψ^Tθ (3)

其中Ψ＝[ψ₁,ψ₂,...,ψ_N]∈R^N×N和θ是投影系数矩阵,

其中ΦΨ^T是编码矩阵，且

A＝ΦΨ^T (4)考虑到成像过程中噪声的影响，将上述模型转化为:

min||θ||₁ s.t.||Aθ-Y||₂＜ε (5)

其中，阈值参数ε控制重建图像对测量数据的保真度，阈值参数ε的值通常设置为低于期望噪声水平。

将L1-norm范数转换为更简单的L2-norm范数，MR图像可以通过优化下列问题重建：

其中λ是正则化参数，用于平衡保真项和正则项的比例。通过求解方程(6)可以得到最优解；并将最优解引入到方程(3)中重建磁共振的原始影像。

基于全变分算法的广义涅斯捷罗夫加速共轭梯度：

在此基础上，引入了TV模型的概念，提出了一种广义涅斯捷罗夫加速算法用于CS-MRI重建共轭梯度法，具体的解决方法如下:

在CS-MRI图像重建模型中加入了TV正则项，并在此基础上建立了图像重建模型，方程式(6)可更改为：

引入整体变分的正则化可以更好地保持图像的边缘信息，利用正则化参数λ₁,λ₂调整权重，选择合适的正则化参数值可以达到理想的重建效果；

采用一种简单的CS-MR图像重建模型，即在方程式(7)中设置稀疏正则项的权重系数λ₁至0，方程式(7)可简化为:

TV算法的梯度公式如下:

被表示为：

其中，x_i,j是图像中像素的值，v_i,j就是TV算法在x_i,j位置上的梯度；

CS-MR图像重建可以通过在现有共轭梯度法的基础上优化上述模型并提出广义涅斯捷罗夫加速算法的概念来解决；

在CG算法中，第k+1次迭代的公式是:

x_k+1＝x_k+a_kd_k (11)

其中x_k是当前迭代点,a_k步长，也称为学习率，d_k是搜索方向；

在每次迭代过程中，CG方法每次迭代的迭代点是在前一次迭代过程中更新的函数位置。根据方程式(9)，步长和搜索方向的确定是为了找到一个合适的功能点位置，因此确定下一次迭代的合适点对于共轭梯度法来说也是一个重要步骤；

结合涅斯捷罗夫加速算法的概念，利用当前迭代点和前一个迭代点的位置，修改共轭梯度法以更新方程式(9)为：

x_k+1＝z_k+1+a_kd_k (12)

z_k+1为：

t_k+1为：

n_k为：

||f(θ)_k||_F和||f(θ)_k-1||_F为目标函数迭代k和k-1次的Frobenius范数；h是一个可调整的权重，k的取值范围为k＞3；

随着迭代次数的增加，n_k≈2，

与来自NAG算法的t_k+1一致；

用GNACG_TV算法重建MR图像：

压缩感知理论可以从变换域中的稀疏图像重建图像。本发明通过螺旋桨采样轨迹获取观测数据，利用GNACG_TV重建图像。采样轨迹如图1所示。GNACG_TV的完整算法流程如下：

基于GNACG_TV算法的CS-MR图像重建算法：

输入：

Y为K空间中欠采样的MR图像数据

ψ为多贝西小波基础

初始化:

k＝1,n₁＝0.2,t₁＝1,z₁＝0,λ＝0.01,h＝5

当k＜max_iter时，进入循环

k＞3

通过方程式(15)计算n_k；

通过方程式(14)计算t_k+1；

通过方程式(13)计算z_k+1；

通过方程式(12)计算x_k+1并更新搜索方向；

k＝k+1；

结束，得到输出：重建的MR图像X＝Ψ^Tθ。

表1.在20％的采样率重建图像的性能对比表

如图2所示，MR图像的腹部，头部和踝部选择作为实验图像，以测试重建性能的CS-MR图像。MR的腹部图像是由南京鼓楼医院的飞利浦ingenia 1.5t核磁共振扫描仪拍摄的。采用江苏省人民医院的联影1.5t umr 560和3.0t umr 770磁共振成像仪采集头部和踝部MR图像。

腹部图像矩阵尺寸＝448*448，头部图像矩阵尺寸＝336*384，脚踝图像矩阵尺寸＝480*480。为了控制不同矩阵大小的影响，在算法处理前对图像进行重采样，三幅图像重采样后的矩阵大小为512*512。本发明采用MSE、PSNR和SSIM三个评价参数来对比评价本发明的MR图像重建的效果。

MSE:这个参数反映了重建图像和原始图像之间的差异程度，定义为:

PSNR:这个参数反映了最大信号量与噪声强度的比值，定义为:

SSIM:该参数用于计算原始图像和重建图像之间的相似性，定义为:

其中

和

分别是参考图像和原始图像之间的亮度、对比度和结构相关函数；α、β和γ分别表示亮度、对比度和结构相关函数的权重；通常，这三个参数的乘积是1。

本发明的图像重建也得到了改善。实施例采用20％的采样率选取样本，并与现有的共轭梯度算法(CG)、基于总变分的共轭梯度算法(CG_TV)和基于总变分的涅斯捷罗夫加速共轭梯度算法(NACG_TV)进行比较。图3-图5中的(a)-(d)显示了使用CG、CG_TV、NACG_TV和GNACG_TV的梯度下降法算法重建的MR图像。图3-图5中(e)-(h)显示上述四种图像重建算法的局部细节。图3-图5中(i)-(l)显示了上述四种重建图像算法的重建误差。为了便于观测，对每幅误差图像进行了5倍增强，图6-8显示了在20％的采样率下使用不同算法进行20次迭代后，MSE、PSNR和SSIM评价参数的收敛曲线。表1列出了采用CG、CG_TV、NACG_TV和GNACG_TV的腹部、头部和踝部MRI数据的MSE、PSNR和SSIM的最终收敛值，采样率为20％。为了验证该算法在不同采样率下的重建效果，比较了腹部MR图像在其他采样率下的重建效果。图9显示重建图像，局部细节图像和重建误差图像的人腹部在40％的采样率。图10显示了腹部重建图像质量评价指标在10％、30％和40％采样率下的差异。由此可见，本发明的算法具有明显优势。

莱斯噪声被加入到磁共振图像中，噪声级由标准差测量。以腹部图像为例，采样率为20％，图11显示了消噪效果在莱斯噪声水平s＝5～25，以及在s＝5～10和条件下的细节图像。

图3-5中的(a)-(d)分别显示了CGCG_TV NACG_TV和GNACG_TV算法在20％的采样率下对人体腹部、头部和踝关节的重建图像。图3-5中的(e)-(h)显示重建图像的局部细节。与其他三种算法相比，GNACG_TV算法可以减少吉布斯伪影，显示更多的局部细节和更清晰的轮廓。图3-5中的(i)-(l)显示重建图像的误差。该算法可以产生较少的误差点，并能更精确地重建图像。图9表明，在较高的采样率下，该算法可以提供更丰富、更准确的重建图像细节。图6-8表明，随着迭代次数的增加，四种算法的最小均方误差逐渐减小，峰值信噪比和峰值信噪比逐渐增大。但是，四种算法的收敛速度存在明显差异，本发明的算法收敛速度明显快于其他三种算法。与其他三种算法相比，本发明的算法具有更好的性能。图10显示了同一图像在不同采样率下的评价参数。该算法具有更好的性能。

表1显示，本发明的算法对于不同数据源和人体不同部位的MR图像具有良好的稳定性和重建性能。图11显示了不同莱斯噪声水平下的去噪效果。本发明的算法的去噪能力明显优于其他三种算法，具有较好的鲁棒性。

本发明的利用图像质量(MSE)、峰值信噪比(PSNR)、结构相似度(SSIM)和去噪效果对GNACG_TV算法进行了定性和定量分析。将该算法与CG、CG_TV和NACG_TV算法进行了比较，证明了该算法的有效性。

以上显示和描述了本发明的基本原理、主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的仅为本发明的优选例，并不用来限制本发明，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。

Claims (4)

1.基于广义涅斯捷罗夫加速共轭梯度算法的压缩采样MR影像重建方法，其特征是包括如下步骤：

在基于磁共振的压缩感知成像技术中，欠采样磁共振图像重建的基本模型是：

min||X||₀ s.t.Y＝ΦX (1)

其中X是重建的MR图像，是Y在k空间的欠采样MR图像数据，Φ是观察矩阵；

Φ被表达为：

Φ＝μF (2)

其中F是一个二维傅里叶变换，μ代表一个欠采样模式；在满足约束等距性的条件下，可以实现图像的精确重建；

上述方程是一个欠定问题，通过范数正则化得到最优解。L0-norm范数优化是一个NP-hard问题；常用的方法是将其转化为最优凸逼近L1-norm范数优化，其中X表示为：

X＝Ψ^Tθ (3)

其中Ψ＝[ψ₁，ψ₂，...，ψ_N]∈R^N×N和θ是投影系数矩阵，

其中ΦΨ^T是编码矩阵，且

A＝ΦΨ^T (4)

考虑到成像过程中噪声的影响，将上述模型转化为：

min||θ||₁ s.t.||Aθ-Y||₂＜ε (5)

其中，阈值参数ε控制重建图像对测量数据的保真度，阈值参数ε的值通常设置为低于期望噪声水平；

将L1-norm范数转换为更简单的L2-norm范数，MR图像可以通过优化下列问题重建：

其中λ是正则化参数，用于平衡保真项和正则项的比例；通过求解方程(6)可以得到最优解；并将最优解引入到方程(3)中重建磁共振的原始影像。

2.根据权利要求1所述的基于广义涅斯捷罗夫加速共轭梯度算法的压缩采样MR影像重建方法，其特征是所述方法还包括如下步骤：

基于全变分算法的广义涅斯捷罗夫加速共轭梯度：

在CS-MRI图像重建模型中加入了TV正则项，并在此基础上建立了图像重建模型，方程式(6)可更改为：

引入整体变分的正则化可以更好地保持图像的边缘信息，利用正则化参数λ₁，λ₂调整权重，选择合适的正则化参数值可以达到理想的重建效果；

采用一种简单的CS-M R图像重建模型，即在方程式(7)中设置稀疏正则项的权重系数λ₁至0，方程式(7)可简化为：

TV算法的梯度公式如下：

被表示为：

其中，x_i，j是图像中像素的值，v_i，j就是TV算法在x_i，j位置上的梯度；

CS-MR图像重建可以通过在现有共轭梯度法的基础上优化上述模型并提出广义涅斯捷罗夫加速算法的概念来解决；

在CG算法中，第k+1次迭代的公式是：

x_k+1＝x_k+a_kd_k (11)

其中x_k是当前迭代点，a_k步长，也称为学习率，d_k是搜索方向；

在每次迭代过程中，CG方法每次迭代的迭代点是在前一次迭代过程中更新的函数位置。根据方程式(9)，步长和搜索方向的确定是为了找到一个合适的功能点位置。

3.根据权利要求1所述的基于广义涅斯捷罗夫加速共轭梯度算法的压缩采样MR影像重建方法，其特征是所述方法还包括确定下一次迭代的合适点的步骤：

结合涅斯捷罗夫加速算法的概念，利用当前迭代点和前一个迭代点的位置，修改共轭梯度法以更新方程式(9)为：

x_k+1＝z_k+1+a_kd_k (12)

z_k+1为：

t_k+1为：

n_k为：

||f(θ)_k||_F和||f(θ)_k-1||_F为目标函数迭代k和k-1次的Frobenius范数；h是一个可调整的权重，k的取值范围为k＞3；

随着迭代次数的增加，n_k≈2，

与来自NAG算法的t_k+1一致；

用GNACG_TV算法重建MR图像：

压缩感知理论可以从变换域中的稀疏图像重建图像；通过螺旋桨采样轨迹获取观测数据，利用GNACG_TV重建图像。

4.根据权利要求1所述的基于广义涅斯捷罗夫加速共轭梯度算法的压缩采样MR影像重建方法，其特征是所述方法还包括GNACG_TV的算法如下：

基于GNACG_TV算法的CS-MR图像重建算法：

输入：

Y为K空间中欠采样的MR图像数据

ψ为多贝西小波基础

初始化：

k＝1，n₁＝0.2，t₁＝1，z₁＝0，λ＝0.01，h＝5

当k＜max_iter进入循环

k＞3

通过方程式(15)计算n_k；

通过方程式(14)计算t_k+1；

通过方程式(13)计算z_k+1；

通过方程式(12)计算x_k+1并更新搜索方向；

k＝k+1；

结束，得到输出：重建的MR图像X＝Ψ^Tθ。

Images