核磁共振图像重构的即插即用算法

李金城，谢朝阳，李金兰，张 涛，邹 健

长江大学 信息与数学学院，湖北 荆州 434023

自20世纪70年代早期开始，核磁共振成像（magnetic resonance imaging，MRI）因其无创性且能有效地描述软组织病变等优点，彻底改变了放射学和医学，成为广泛应用于临床诊断的方法。限制MRI应用的主要因素是其数据采集过程比其他技术慢。在不降低重构质量的前提下，减少数据采集量一直是MRI成像研究的重点[1-3]。在过去的几十年里，基于稀疏先验的压缩感知（compressed sensing，CS）利用图像在某些变换域（如傅里叶变换、小波变换、梯度变换）中的稀疏性，精确地从少量随机测量数据中重构稀疏信号。压缩感知核磁共振成像（CS-MRI）使得从高度欠采样的K空间数据重构核磁共振图像成为可能，并成为多种重构方法的基础[4-5]。在MRI模型中，获取的K空间数据y∈Rm通常被表示为：

其中，x∈Rn为待重构图像；ηe∈Rm为测量噪声或干扰；A∈Rm×n(m＜n)是压缩感知采样矩阵。通常A=P×F,P∈Rm×n为欠采样模板，如随机采样模板、径向采样模板、笛卡尔采样模板等；F∈Rn×n为稀疏变换算子，如小波变换、傅里叶变换等。此时获取的K空间数据y要比观测核磁共振图像x的规模小得多。从低维y中重构唯一的高维图像x可以通过求解如下优化问题得到：

MRI中最常用的正则函数是L1正则，即φ(x)=L1范数的凸性可以保证问题（2）得到高效的求解。但是L1正则得到的解是真实信号的一个有偏估计，往往会低估x的高振幅分量[6]。非凸正则通常比L1正则更能有效地促进稀疏和更准确地估计x的高振幅分量。常用的非凸正则包括：光滑剪切绝对偏差（smoothly clipped absolute deviation，SCAD）正则，Lq范数φ(x)=‖x ‖q(0＜q＜1)等[7]。然而这些非凸正则，会使得问题（2）的目标函数是非凸的，使其陷入次优的局部最小值，增加了求解难度。Selesnick等人通过引入L1范数的Moreau包络和最小最大凹罚函数，构造了凸-非凸（convex-nonconvex，CNC）稀疏正则[6，8-9]，其表达式如下：

MRI重构的另一个重要问题是重构算法。近似点梯度下降（proximal gradient descent，PGD）是一种高效的求解非光滑问题（2）的一阶算法，其基本原理是利用不同策略来处理目标函数的光滑和非光滑部分，对光滑部分使用梯度迭代，对非光滑部分利用近似点算子求解[11-12]。近年来，随着深度学习的兴起，将传统优化算法与深度神经网络结合的数据驱动算法得到快速发展。即插即用（plug-and-play，PnP）算法将用合适的去噪器来替代传统稀疏优化算法中的近似点算子[13-14]，其优点是当没有足够的数据进行端到端训练时，传统或者预先训练好的去噪器可以作为即插即用框架下的模块化部件灵活插入使用。常用的去噪器包括传统的手工设计去噪器，如全变差（total variation，TV）、三维块匹配协同滤波（block-matching 3D，BM3D）[15]等。近年来，深度神经网络去噪器可以从大量自然图像中学习到更为丰富的图像先验，显示了优异的去噪性能[13，16]。基于深度学习的即插即用算法在图像重构、图像去噪、语音去噪和运动设备等方向得到了广泛的应用[17-19]。

本文利用即插即用近似点梯度下降算法来求解带凸-非凸正则的核磁共振重构模型，其主要贡献如下：

（1）给出了凸-非凸正则的近似点算子的迭代解形式，解决了采用凸-非凸正则的优化问题求解过程中较为困难的部分，从而提高计算效率。

（2）提出求解基于凸-非凸正则的核磁共振重构模型的近似点梯度下降算法，进而利用去噪器替换近似点算子，得到求解基于凸-非凸正则的核磁共振重构模型的即插即用的近似点梯度下降算法。其中即插即用算法可以灵活地引入更多先进的去噪器，可以有效提升去噪效果。另外在MRI重构中，采用凸-非凸正则能够更精确地重构图像在频域中的高频分量，进而能减少重构图像边缘轮廓和细小纹理等细节信息的丢失。

1 相关工作

定义1对于任意的α＞0，对于给定的适当、闭、凸函数g(x)，其近似点算子定义为：

其中，xi为x的第i个元素。从软阈值的定义可以看出，当xi的绝对值小于某一阈值α时，软阈值函数会将xi的值变为0，因此软阈值函数可理解为一个去噪器，参数α控制其去噪强度。

由近似点算子的定义，可得到近似点梯度下降算法的一般形式。考虑复合优化问题：

其中，f(x)凸且可微的，g(x)凸但不可微，且g(x)的近似点算子容易求得。

近似点梯度下降算法的基本迭代形式如下：

其中，∇f(x)是f(x)的梯度函数，α是步长，g的近似点算子proxαg如式（4）所定义。

2 凸-非凸稀疏正则的近似点算子

凸-非凸稀疏正则φb(x)的近似点算子定义如下：

由于φb(x)形式较为复杂，此时proxλφb(y)不再具有封闭解，本文将给出其迭代解。

证明 由近似点算子的定义可得Sb(x)的最优值点为：

代入Sb(x)中，则有：

容易求得Sb(x)的梯度如下：

命题1得证。

命题2 φb(x)的近似点算子为：

证明 将式（8）中的目标函数记为：

由近似点梯度下降算法，可得：

取步长μ=1，命题2得证。

3 即插即用近似点梯度下降算法

3.1 近似点梯度下降算法

本节将考虑运用近似点梯度下降算法求解基于凸-非凸稀疏正则的核磁共振图像重构模型（2）。令式（2）中,结合命题2，求解基于凸-非凸稀疏正则的近似点梯度下降（proximal gradient descent algorithm for convex-nonconvex sparse regularization，PGD-CNC）算法的主要步骤如下：

步骤1

步骤2

3.2 即插即用近似点梯度下降算法

此时若选择两个适当的去噪器Hσ1、Hσ2分别替换近似点算子、proxαλ‖·‖1，其中σ1、σ2为去噪参数，值越大去噪强度越大，则式（19）可改写为：

综合式（18）、（20），可得到求解基于凸-非凸稀疏正则的即插即用近似点梯度下降（plug-and-play proximal gradient descent algorithm for convex-nonconvex sparse regularization，PnP-PGD-CNC）算法如算法1所示。

算法1 MRI图像重构的PnP-PGD-CNC算法

即插即用的近似点梯度下降算法与近似点梯度下降算法的区别：前者可将预先训练好的去噪器直接插入到传统近似点梯度下降算法中的一个迭代步中。当去噪器满足一些简单条件，如非扩张、收缩时，可以保证算法的收敛性[17，19]。

算法1中计算复杂度较高的是z(k+1)的迭代，涉及到多个矩阵相乘，但在实际实验中都可以进行快速计算。由于P是欠采样矩阵，则PTP为矩阵元素间的乘积。若F为傅里叶变换，则有FT=F-1,F和FT在实验中可通过快速傅里叶变换（fast Fourier transform，FFT）和逆快速傅里叶变换（inverse fast Fourier transform，IFFT）实现。

4 数值实验与结果分析

4.1 实验设置

压缩感知核磁共振成像通过下采样获取更少的数据来加速核磁共振成像。本实验选取了三种不同的核磁共振图像：大脑、脑部血管和上半身。所有图像灰度值取值范围为0～1，尺寸设置为256×256像素；ηe为测量噪声且满足ηe～N(0,σeIk)，噪声级别为σe=15/255；P是欠采样模板，本实验应用了采样率均为30%的欠采样随机采样模板、径向采样模板和笛卡尔采样模板；F为傅里叶算子。待重构图像和采样模板如图1所示。

为了验证PnP-PGD-CNC算法的有效性，实验将ADMM-Net[20]、BM3D-MRI[21]、基于L1正则的近似点梯度下降算法（PGD-L1）[17]、基于L1正则的即插即用近似点梯度下降算法（PnP-PGD-L1）[17]以及本文提出的基于凸-非凸稀疏正则的近似点梯度下降算法（PGD-CNC）、基于凸-非凸稀疏正则的即插即用的近似点梯度下降算法（PnP-PGD-CNC）的核磁共振重构图像进行对比。其中ADMM-Net[20]是将传统的交替方向乘子法（alternating direction method of multipliers，ADMM）展开到深度神经网络中；BM3D-MRI[21]是使用非局部的BM3D图像模型进行MRI重构；其他四种算法的主要迭代步骤如表1所示。

表1中PGD-L1算法为正则项φ(x)=‖x‖1时，使用PGD算法得出的迭代解，其解由软阈值函数给出[17]；PnP-PGD-L1算法为文献[17]中提到的，由上述PGD算法与PnP算法相结合得出的迭代解；PGD-CNC算法即由式（18）、（19）给出；PnP-PGD-CNC算法即本文所提出的算法（算法1）。

表1 算法的迭代式Table 1 Iterative formula of algorithms

在PnP-PGD-L1和PnP-PGD-CNC算法中选取BM3D[15]和五种神经网络去噪器作为PnP框架下不同的去噪器。表2简要罗列了各种去噪器的特点。

表2 PnP框架下的去噪器Table 2 Denoiser under PnP framework

特别地，文献[25]中提到DnCNN需要为每个噪声水平分别学习一个模型，其他的去噪器通过单个模型可以处理各个噪声水平的噪声图片。因此，在PnP-PGDCNC算法中选用DnCNN去噪器时，去噪参数σ1=5σe/3=25/255,σ2=σe/3=5/255,在其余算法中σ1=σ2=σe=15/255。本实验中，参数设置均是在大量的实验下选择图像恢复效果较好时的参数，迭代次数均设置为各算法在50次内效果最好时的迭代次数。五种神经网络去噪器由多种不同的自然数据集训练得到[22-25]。实验操作都是在Intel®CoreTMi7-8750HM、2.20 GHz CPU、16 GB内存，NVIDIA GeForce GTX 1060 6 GB显卡的PC上进行，编程环境使用Python3.8，深度学习框架使用Pytorch。本文的主要实验代码、参数设置和结果可以从网站链接https：//github.com/zj15001/PNP-PGD-CNC下载。

本文实验结果由两方面进行呈现，第一方面是视觉效果图像对比，第二方面是数值结果比较。

为了从视觉效果上直观的评价重构效果，本文分别给出了表1中不同算法的重构效果图，以及重构后效果图与待重构核磁共振图像的差值图像。为方便更好地显示差异，选取图像的两个局部细节进行放大显示，并将差值图像的灰度值取值范围设置为0～0.2。

数值结果比较采用峰值信噪比(PSNR)的数值大小来作为评价核磁共振图像重构效果的重要指标，其定义如下：

4.2 实验结果分析

首先，采用随机采样模板对大脑图像进行采样后加噪声，再运用表1中的四种算法进行图像重构，具体重构图像和差值图像如图2所示。其中，每张子图的第一排为重构图像，第二排为重构图像与待重构图像的差值图像，每排下方的两个小图则是将两个局部细节进行放大显示；PnP-PGD-L1和PnP-PGD-CNC算法的重构结果中显示了部分去噪器（BM3D、FFDNet、IRCNN、DRUNet）的重构图像和差值图像。

由图2中各算法的重构图像及其差值图像的对比可以看出，PGD-L1算法的重构效果最差，图像分辨率较低，保留了较少的边缘和细节信息。PnP-PGD-L1和PGD-CNC算法相较PGD-L1算法都有不同程度的性能提升。本文提出的PnP-PGD-CNC算法整体要比其他算法的重构效果更好，其中效果最好的是采用了DRUNet去噪器的图2（j）。其左上角方框中的边缘线条和纹路更加清晰，右下角方框中的棱角部分更加分明，保留了更好的边缘轮廓信息，且差值图的颜色明显变深，表明了重构图像与原始图像的一致性较高。从数值结果来看，PGD-L1算法的重构图像（图2（a））峰值信噪比为23.78 dB，而在PnP-PGD-CNC框架下使用DRUNet去噪器的重构图像（图2（j））峰值信噪比为31.60 dB，显然PnP-PGD-CNC的视觉效果和数值结果更好。

最后，表3给出对于不同图像使用不同算法进行图像重构的数值比较结果。实验数据表明，对于本文给出的三种不同采样模板，PnP-PGD-L1算法和PGD-CNC算法的重构效果都要明显高于PGD-L1算法，分别突出了即插即用算法和凸-非凸稀疏正则的有效性。同时对比PnP-PGD-L1算法和PnP-PGD-CNC算法中使用不同去噪器时，神经网络去噪器的重构效果要优于其他去噪器，能够得到更高的峰值信噪比值，其中峰值信噪比值最高的是选用DRUNet去噪器。当PnP-PGD-CNC算法选用DRUNet去噪器时，较PGD-L1算法峰值信噪比值平均提升了6.26 dB。同时从PSNR值的大小也可以看出PnP-PGD-CNC算法的重构效果明显优于ADMMNet和BM3D-MRI算法。

表3 图像重构的峰值信噪比值Table 3 Peak signal-to-noise ratio of image reconstruction 单位：dB

5 结束语

本文提出了一种基于凸-非凸稀疏正则和即插即用近似点梯度下降的核磁共振图像重构算法，主要有以下特点：

（1）该算法是基于凸-非凸稀疏正则的，可以保证目标函数的整体凸性，并且能够准确地估计高振幅分量。本文给出了该正则项的近似点算子的迭代解形式，这为之后基于凸-非凸稀疏正则的优化问题的计算提供了便利。

（2）该算法中融入的即插即用算法，是将传统优化算法与深度学习相结合，与传统图像重构方法不同，它可以通过去噪器隐式定义即插即用先验。该框架具有足够的灵活性，可以与各种最先进的去噪器兼容。同时可以从理论上保证算法的收敛性，体现出其便捷、高效等优点。

为验证本文所提算法的有效性，与另外几种算法的重构结果进行比较。从视觉效果图像和数值(PSNR)结果比较都可以得出：本文所提算法对边缘轮廓和细小纹理的处理具有较为显著的性能提升，并且可以得到更高的峰值信噪比。

进一步的工作研究主要分为三方面：一是将继续探究本文所提出算法的收敛性理论证明；二是考虑设计更为高效的基于深度学习的去噪器，以达到更好的去噪效果；三是将凸-非凸稀疏正则和即插即用算法融入到其他的凸优化算法进行图像重构的研究。