基于深度学习的MR医学图像重建方法研究

Study on Deep Learning-Based MR Medical Image Reconstruction Method

摘要: 磁共振(Magnetic Resonance, MR)成像是一种先进的成像方式，其特殊的成像方法导致采集数据的时间较长，容易受噪声的干扰；通过忽略部分数据(欠采样)可以提高采集速度，但会使图像分辨率丧失或产生图像伪影，导致图像质量下降。近年来，即插即用图像复原可以灵活地处理不同的图像复原任务，被广泛应用在自然图像中。受此启发，本文提出了一种基于深度学习的MR图像重建框架，通过训练卷积神经网络模型来建立深度去噪先验，并将其嵌入到交替迭代优化算法中，进而隐式地对图像先验进行建模。以解决欠采样MR图像中的噪声和模糊问题，从而重建出高质量图像。和其他先进算法相比，我们的方法在MR图像重建上具有更好的竞争力。

Abstract: Magnetic Resonance (MR) imaging is an advanced imaging method, its special imaging method leads to a long time to collect data, easy to be disturbed by noise; by ignoring part of the data (undersampling), the acquisition speed can be improved, but the image resolution will be lost or image artifacts will be generated, resulting in image quality degradation. In recent years, Plug-and-Play image restoration has been widely used in natural images because it can flexibly handle different image restoration tasks. Inspired by this, this paper proposes a deep learning-based MR image reconstruction framework, which establishes a deep denoising prior by training a convolutional neural network model and embedding it into an alternating iterative optimization algorithm to implicitly model the image prior. To solve the problem of noise and blur in undersampled MR image and reconstruct high-quality images. Compared with other advanced algorithms, our method has better competitiveness in MR image reconstruction.

1. 引言

在医学成像中，图像在传输和采集过程中容易受到不同类型和不同程度的噪声影响 [1] [2] 。MRI是一种先进的成像方式，提供了优越的软组织对比度，但由于采集时间较长，除了受到噪声的影响，患者的不经意移动等还可能造成图像模糊，从而降低了图像的质量，导致人类或机器产生错误的决策，使得诊断、分析和治疗过程变得困难 [3] [4] 。因此，加快图像采集，消除图像重建过程中产生的噪声和伪影变得至关重要。

近年来，深度学习已经被证明是一个强大的应用工具 [5] [6] ，卷积神经网络(Convolutional Neural Network, CNN)和生成对抗网络(Generative Adversarial Networks, GAN)在许多机器学习任务中取得了最先进的性能 [7] [8] 。卷积神经网络由于能有效地提取局部图像特征而获得了较高的性能，在医学图像重建中被广泛应用 [9] [10] ，能够很好地消除噪声和伪影。

深度即插即用迭代算法(Plug-and-Play, P&P)的主要思想是借助可变分裂算法，例如交替方向乘子法(Alternating Direction Method of Multipliers, ADMM) [11] 和半二次方分裂(Half Quadratic Splitting, HQS) [12] ，可以分别处理数据项和先验项 [13] 。最近，关于P&P图像恢复的工作表明 [14] ，去噪模型可以隐式地作为基于模型的方法的图像先验来解决许多逆问题。当通过具有大建模能力的深度CNN网络有区别地学习去噪器时，这种特性对于P&P图像恢复具有相当大的优势 [15] 。然而，虽然更深和更大的CNN模型正迅速得到普及，但是现有的P&P图像重建方法由于缺乏合适的降噪模型先验而阻碍了其性能。随后，Zhang等人 [16] 将深度学习模型引入作为去噪先验，在自然图像重建中得到很好地应用。

受上述研究启发，本文提出了一种基于深度学习的MR图像重建框架，该框架通过训练高度灵活和有效的CNN模型来建立深度去噪先验，引入P&P的交替迭代优化算法和自动调整机制，以解决MR图像重建中的去噪、去模糊问题，重建出高质量的图像，在速度和性能优势上具有很好的优势。

2. 相关工作

P&P图像重建算法在解决传统的自然图像逆问题中被广泛利用 [17] [18] [19] [20] 。Danielyan等人 [21] 利用Nash均衡推导出一种迭代解耦去模糊BM3D (IDDBM3D)方法用于图像去模糊。Venkatakrishnan等人 [15] 提出了第一个将去噪模型视为“P&P先验”的工作。在此之前，Zoran等人 [19] 也提到了类似的即插即用思想，其中，采用HQS算法被用于图像去噪、去模糊；Chambolle等人 [22] 使用了ADMM和HQS的替代方法，即先验对耦算法来解耦数据项和先验项；Teodoro等人 [23] 将特定高斯混合模型(GMM)去噪器插入到ADMM中，以求解图像去模糊和压缩成像问题。尽管这些方法在性能上有一定的突破，但速度上还有待提高。P&P图像重建算法可以通过隐式地定义去噪先验，当去噪先验为深度学习模型时，基于模型的方法将获得速度和性能优势 [24] 。

随着网络设计和基于梯度的优化算法等深度学习技术的发展，基于CNN的去噪模型在效率等方面都表现出了良好的性能 [25] [26] 。在此之后，一系列基于CNN的即插即用的重建工作被提出 [27] [28] [29] 。例如，Romano等人 [27] 提出了用TNRD去噪器能够对图像去模糊；Tirer等人 [28] 提出了使用IRCNN去噪器的迭代去噪和反向投影，用于去噪和去模糊。即插即用的图像重建可以利用深度学习模型来做先验去噪，现有的方法 [30] [31] [32] [33] 模型普遍使用DnCNN或IRCNN。DnCNN作为去噪先验往往需要大量的迭代才能获得令人满意的性能，这将增加计算负担。IRCNN去噪模型可以处理广泛的噪声水平，但它对特定的噪声水平及模糊图像处理上鲁棒性不高。目前，最适合即插即用重建的去噪器是FFDNet [18] ，它可以通过将噪声水平图作为输入来处理广泛的噪声水平，但缺乏提高即插即用重建性能的有效性。Zhang等人 [16] 利用使用的U-Net [34] 和ResNet [35] 进行架构设计来改进FFDNet，在自然图像处理上有着显著成绩，但并未应用在医学图像处理上。

基于此，本文在文献 [16] 基础上加以改进，提出一种用于MR医学图像重建的深度学习模型，加深网络层数提高模型的鲁棒性，并加入一个自动调整机制来自动设置参数，通过实验证明了其与特定任务技术方法相比在图像去噪和去模糊方面更具有的竞争力，在速度和性能上具有很好的优势。

3. 方法理论

目前，U-Net已成为医学图像分析领域研究最为广泛的模型，ResNet能通过短跳跃连接解决梯度消失的问题，结合U-Net和ResNet的类似想法可以在其他作品中找到 [36] [37] ，故，我们将残差块集成到U-Net中，以含有噪声的欠采样模糊图作为输入，以实现有效的去噪先验建模。整个网络结构主要由CNN和P&P方法两部份组成。在3.1节对主干网络的CNN模型进行概述；3.2节对即插即用方法进行简要介绍；3.3节为本文所提的方法的详细描述。

3.1. 主干网络

本文所提出的网络架构如图1所示。我们的模型有能力通过单一模型处理各种噪声水平，由六层组成，每层在2 × 2步长卷积(SConv)下采样和2 × 2转置卷积(TConv)上采样操作之间都有一个连接。从第一层到第六层，每一层的通道数分别为16、32、64、128、256和512个。每层都采用了连续的残差块，每个残差块包含一个Conv激活函数，由于对于有偏差的网络，偏差的幅度将远远大于滤波器，这可能会损害其通用性，故我们在所有Conv、SConv和TConv层中都没有使用偏置项。

Figure 1. Deep learning-based MR medical image reconstruction network architecture

图1. 基于深度学习的MR医学图像重建网络架构

3.2. 即插即用方法

在文献 [28] 中详细介绍了即插即用的方法，一般情况下，图像重建的目的是从脏图中恢复出干净的图像，如公式(1)所示。

$\begin{array}{l}Y=\tau \left(x\right)+n\\ \tau \left(x\right)=Hx\end{array}$ (1)

其中，x为未知的原始图像，y为观测后的图像，H是一个退化矩阵，n为加性高斯白噪声(AWGN)， $\tau \left(x\right)$ 通过指定不同的退化操作可以相应的得到不同的图像重建任务。在这里，我们的操作是欠采样导致图像模糊，故，H是一个稀疏矩阵。几乎任何恢复x的方法都涉及代价函数(由保真度和惩罚项组成)，它由期望得到最小化解。其代价函数表示为：

$f\left(\tilde{x}\right)=\frac{1}{2{\sigma }_e^2}{\Vert y-H\tilde{x}\Vert }_2^2+s\left(\tilde{x}\right)$ (2)

其中， $\tilde{x}$ 为优化变量， ${\Vert .\Vert }_2$ 为欧氏范数。然而，利用变量分裂，P&P方法将最小化问题描述为：

${\min }_{\tilde{x},\tilde{v}}l\left(\tilde{x}\right)+\beta s\left(\tilde{v}\right)\text{s}\text{.t}\text{.}\tilde{x}=\tilde{v}$ (3)

其中， $l\left(\tilde{x}\right)\triangleq \frac{1}{2{\sigma }_e^2}{\Vert y-H\tilde{x}\Vert }_2^2$ 是公式(2)中的保真度项，而β是一个正参数，使得代价函数变得复杂。通过ADMM [11] 构造一个增广拉格朗日可以解决该问题，即：

$\begin{array}{c}{L}_{\lambda }=l\left(\tilde{x}\right)+\beta s\left(\tilde{v}\right)+u^{T\left(\tilde{x}-\tilde{v}\right)}+\frac{\lambda }{2}{\Vert \tilde{x}-\tilde{v}\Vert }_2^2\\ =l\left(\tilde{x}\right)+\beta s\left(\tilde{v}\right)+\frac{\lambda }{2}{\Vert \tilde{x}-\tilde{v}+\tilde{u}\Vert }_2^2-\frac{\lambda }{2}{\Vert \tilde{u}\Vert }_2^2\end{array}$ (4)

其中，u是对偶变量， $\tilde{u}\triangleq \frac{1}{\lambda }u$ 是缩放对偶变量， $\lambda$ 是ADMM的惩罚参数，ADMM算法在以下三个步骤中进行迭代直到收敛。

$\begin{array}{l}{x}_k^{\vee }=\underset{\tilde{x}}{\arg \min }{L}_{\lambda }\left(\tilde{x},v_{k-1}^{\vee },u_{k-1}^{\vee }\right)\\ v_k^{\vee }=\underset{\tilde{v}}{\arg \min }{L}_{\lambda }\left(x_k^{\vee },\tilde{v},u_{k-1}^{\vee }\right)\\ u_k^{\vee }=u_{k-1}^{\vee }+\left(x_k^{\vee }-v_k^{\vee }\right)\end{array}$ (5)

结合公式(4)和(5)，则有：

$\begin{array}{l}{x}_k^{\vee }=\underset{\tilde{x}}{\arg \min }l\left(\tilde{x}\right)+\frac{\lambda }{2}{\Vert \tilde{x}-\left(v_{k-1}^{\vee }-u_{k-1}^{\vee }\right)\Vert }_2^2\\ v_k^{\vee }=\underset{\tilde{v}}{\arg \min }\frac{\lambda }{2\beta }{\Vert \left(x_k^{\vee }+u_{k-1}^{\vee }\right)-\tilde{v}\Vert }_2^2+s\left(\tilde{v}\right)\\ u_k^{\vee }=u_{k-1}^{\vee }+\left(x_k^{\vee }-v_k^{\vee }\right)\end{array}$ (6)

s(x)是通过选择 $D\left(·;\sigma \right)$ 来隐式定义的，故，P&P方法不需要知道或显式地指定先验函数s(x)，详细的算法描述见表1所示。P&P方法需要大量的迭代，这通常要求P&P能够收敛到一个好的解决方案，以及需要设计合适的 $\beta$ 和 $\lambda$ 参数值。

Table 1. Plug-and-Play algorithm

表1. 即插即用算法

3.3. 我们的方法

众所周知，图像重建是一个逆问题，在此之前，还需要采用正则化来约束空间解 [38] 。从贝叶斯角度来看，通过求解最大后验估计(Maximum APosteriori, MAP)可以得到 $\hat{x}$。

$\hat{x}=\underset{x}{\arg \max }\log p\left(y|x\right)+\log p\left(x\right)$ (7)

其中， $\log p\left(y|x\right)$ 表示观测值y的对数似然值， $\log p\left(x\right)$ 表示干净图像x的先验，并且与退化图像y无关。公式(7)可以重新表述为：

$\hat{x}=\underset{x}{\arg \min }\frac{1}{2{\sigma }^2}{\Vert y-\tau \left(x\right)\Vert }^2+\lambda R\left(x\right)$ (8)

该公式最小化了由一个数据项 $\frac{1}{2{\sigma }^2}{\Vert y-\tau \left(x\right)\Vert }^2$ 和一个正则化参数为 $\lambda$ 的正则化项 $\lambda R\left(x\right)$ 组成的能量函数。为了解耦前面所述的数据项和先验项，我们引入了一个辅助变量z，得到一个约束优化公式：

$\hat{x}=\underset{x}{\arg \min }\frac{1}{2{\sigma }^2}{\Vert y-\tau \left(x\right)\Vert }^2+\lambda R\left(x\right)\text{s}\text{.t}\text{.}z=x$ (9)

通过最小化处理可以得到：

$L_u\left(x,y\right)=\frac{1}{2{\sigma }^2}{\Vert y-\tau \left(x\right)\Vert }^2+\lambda R\left(z\right)+\frac{u}{2}{\Vert z-x\Vert }^2$ (10)

其中，µ是一个惩罚因子。可以通过迭代解决x和z的子问题，同时保持其余的变量固定。

$\{\begin{array}{l}{x}_k=\underset{x}{\arg \min }{\Vert y-\tau \left(x\right)\Vert }^2+\mu {\sigma }^2{\Vert x-z_{k-1}\Vert }^2(\text{a})\\ z_k=\underset{z}{\arg \min }\mathrm{a}\frac{1}{2{\left(\sqrt{\lambda /\mu }\right)}^2}{\Vert z-x_k\Vert }^2+R\left(z\right)(\text{b})\end{array}$ (11)

公式(11)将数据项和先验项解耦为两个独立的子问题。具体地说，公式(11a)的目的是找到 $z_{k-1}$ 的近端点，通常有一个快速的闭型解，取决于 $\tau$ ；而从贝叶斯角度来看，(11b)的子问题对应于噪声水平为 $\sqrt{\lambda /\mu }$ 的高斯噪声 $x_k$。因此，任何高斯去噪器都可以插入到交替迭代中来求解公式(8)。为了解决这个问题，我们将公式(11b)整理为：

$z_k=Denoiser\left(x_k,\sqrt{\lambda /\mu }\right)$ (12)

故，我们利用CNN的高灵活性、GPU的高效率和强大的深度网络建模能力等优点，学习单一的CNN去噪器来替换(12)。通过加入先验去噪CNN模型，引入自动调整机制设置合适的参数值，能够有效解决P&P存在的问题。对具有均匀模糊(或图像反褶积)的模糊图像进行去模糊的退化模型通常表示为：

$y=x\otimes k+n$ (13)

其中， $x\otimes k$ 表示干净图像与模糊核之间的二维卷积。假设在圆形边界条件下进行卷积，则得到公式11(a)的快速解：

$x_k=F^{-1}\left(\frac{\overline{F\left(k\right)}F\left(y\right)+{\alpha }_{k}F\left(z_k-1\right)}{\overline{F\left(k\right)}F\left(k\right)+{\alpha }_k}\right)$ (14)

其中， $F(.)$ 和 $F^{-\text{1}}(.)$ 表示快速傅里叶变化(Fast Fourier Transform, FFT)和逆傅里叶变化， $\overline{F(.)}$ 表示 $F(.)$ 的复共轭。可以注意到，模糊内核k只涉及到公式(14)，明确地处理了模糊的失真。具体详细算法描述见表2所示。

Table 2. CNN-based algorithm

表2. 基于CNN的算法

4. 实验分析

4.1. 数据集与实验细节

训练数据：我们使用来自MGH-USCHCP19公共数据库 [39] 的大脑图像作为训练数据集，采用T1加权三维MRI采集协议MPRAGE，重复时间TR = 2530 ms，回声时间TE = 1.15 ms，反转时间TI = 1100 ms，带宽BW = 651 Hz Px-1在西门子Skyra3TMRI平台(西门子医疗解决方案，德国埃尔兰根)。使用35名受试者的冠状面t1加权切片生成一个6125张图像数据集(每个受试者产生175张切片)。对于每张图像，中央的256 × 256像素被裁剪并下采样到128 × 128。

测试数据：第一次评估实验中使用的测试数据取自训练对象之外的另一个受试者。测试数据还包括t2加权MRI采集协议SPACE，TR = 3200 ms，TE = 561 ms，FOV = 224 mm × 244 mm，BW = 744 Hz Px-1，在西门子Skyra3TMRI平台上，并将其调整为128 × 128来测试模型。实验表明，我们的模型受益于图像样本量的增加和大批量的批处理。

实验细节：我们首先对干净图像进行下采样处理，然后添加噪声水平为 $\sigma$ 的AWGN来合成模糊图像。实验通过WinSCP和PuTTY连接GPU服务器，使用NVIDIA GTX 2080Ti GPU(24G)在公共Pytorch框架上实现。其中，学习率为5e−4，迭代200个epoch，噪声率为25。并通过峰值信噪比(Peak Signal to Noise Ratio, PSNR)作为模型的评价指标，用于客观地分析模型的鲁棒性，公式如下。

$\text{PSNR}=10{\log }_{10}\left(\frac{{\text{MAX}}_1}{\sqrt{\frac{1}{mn}{\displaystyle {\sum }_{i=0}^{m-1}{\displaystyle {\sum }_{j=0}^{n-1}{\left[I\left(i,j\right)-K\left(i,j\right)\right]}^2}}}}\right)$ (15)

其中，MAX_I是表示图像点颜色的最大数值，I表示大小为m × n的干净单色图像，K为噪声近似值，峰值信噪比越高，则算法性能越好。

4.2. 实验结果与分析

为了评估本文方法的有效性，首先选择了一位受试者的水平面T1加权切片进行测试，分别对噪声图像、欠采样图像及含有噪声的欠采样图像进行重建，如图2所示。其中，噪声水平设置为25，欠采样图像通过模糊算子进行模拟(丢失40%的像素)。

Figure 2. Examples of reconstruction results of this method

图2. 本文方法的重建结果示例

最左边为原始的干净图像，其次上半部分从左到右分别是噪声图像、欠采样模糊图像及欠采样噪声图像，下半部分本文所提出的方法重建出的图像。可以看出，在处理同时含有噪声和模糊的图像时，其峰值信噪比(PSNR)达到33.19 db，能够很好地处理较为复杂的情况。

此外，本文选择了受试者的冠状面、失状面及水平面的T1加权切片(如图3所示)和水平面及矢状面的T2加权切片(如图4所示)进行测试，与先进的方法(P&P [16] 和IRCNN [30] )进行比较。“our”表示本文提出的方法，粗体表示最优结果。

如图3所示，第一行为T1加权平面图，本文方法的PSNR能够达到33.19 db，比P&P增益了0.51 db，比IRCNN增益了0.45 db；第二行为矢量图，其PSNR比P&P增益了0.39 db，比IRCNN模型增益了0.12 db；第三行为冠状图，相比之下，IRCNN的PSNR为最优结果，比本文方法增益了0.08 db。

Figure 3. Reconstruction results of T1-weighted images with different methods

图3. 不同方法下对T1加权图像的重结果

为了更好地分析我们模型的适用性，使用T2加权切面图像进行测试，如图4所示。从第一行的平面图中我们可以看到，本文方法的PSNR比P&P方法增益了0.92 db，比IRCNN算法增益了0.49 db；从第二行的冠状图中可以看出，本文的PSNR为31.23 db，比P&P方法增益了0.05 db，比IRCNN增益了0.18 db。

Figure 4. Reconstruction results of T2-weighted images with different methods

图4. 不同方法下对T2加权图像的重建结果

综上所述，无论哪种切片模式下，本文的方法都具有良好的竞争力。尽管在T1加权冠状图的测试上IRCNN的峰值信噪比比本文方法增益了0.08 db，但无论在模型架构还是速度上，本文方法都具有竞争力。通过引入自动调整机制，能在较少的调整参数情况下进行训练。和IRCNN方法相比，本文方法在速度和性能上占据优势。

5. 结论

在本文中，我们介绍了基于CNN的MR医学图像重建算法，结合交替迭代优化算法，能够在减少迭代次数的情况下重建出高质量图像，有效地减少了训练时间。引入自动调整机制使得我们的方法比P&P需要调整的参数少。通过实验证明，我们的方法在去除噪声、模糊等方面与其他算法相比，更具有竞争力，在没有特定任务训练的情况下，它比基于学习的方法更灵活。

基金项目

国家自然科学基金(11963003，12273007)；SKA专项资助(2020SKA0110300)；贵州省科技厅(黔科合基础-ZK[2022]一般143)；国家重点研发计划(2018YFA0404602，2018YFA0404603)；贵州省教育厅青年科技人才成长项目(黔教合KY字[2018] 119号)；贵州大学引进人才科研基金(贵大人基合字(2018) 60号)；贵州大学培育项目(贵大培育[2020] 76号)。

NOTES

^*通讯作者。