基于递归残差U型网络的快速磁共振成像方法   0    0

[0001] 本发明属于磁共振成像领域，涉及一种基于递归残差U型网络的快速磁共振成像方法。

[0002] 磁共振成像(Magnetic Resonance Imaging，MRI)技术因其具有较高的软组织分辨率，对人体没有电离辐射损伤等优点在临床应用中非常受欢迎。但是由于MRI的成像速度较慢，成像过程中被检者的生理性运动等往往会造成成像伪影，难以满足实时成像的要求，因此如何加快MRI的成像速度是MRI理论和技术研究的热点之一。减少MRI的k空间数据采集总量是一种有效的加速方法，但是数据采集量的大量减少，会带来图像质量的显著下降，尽管研究者可通过多种重建算法，例如，压缩感知(Compressed Sensing，CS)技术等，提高欠采样图像重建的质量，但是往往需要较长的重建时间，难以满足实时重建的临床需求。

[0003] 2015年，由于卷积神经网络具有离线训练、快速成像的优点，研究者开始将它应用到快速磁共振成像研究中。卷积神经网络利用大量的数据对网络进行训练，生成参数优化后的卷积神经网络，再将需要重建的欠采样数据利用训练好的网络进行快速高质量重建，已成为这两年MRI快速成像领域的研究热点。

[0004] 目前已申请的基于卷积神经网络的快速MRI成像方面的专利有：磁共振成像方法和系统(申请号：CN201310633874.6)，基于深度网络模型来估计k空间中已采样点和未采样点之间的映射关系，从而估计出完整的k空间数据来重建磁共振图像；一种基于深度卷积神经网络的快速磁共振成像方法及装置(申请号：CN201580001261.8)，提出基于深度卷积神经网络来实现快速磁共振成像；一种磁共振成像方法及装置(申请号：CN201710236330.4)，提出基于多对比度的欠采和全采MRI图像来训练网络；一种基于卷积神经网络的多对比度磁共振图像重建方法(申请号：CN201711354454.9)，提出利用多对比度MRI图像来训练卷积神经网络；一种磁共振成像方法和系统(申请号：CN201611006925.2)，基于深度学习方法提高CS‑MRI图像重建的质量和速度；基于机器学习的并行磁共振成像GRAPPA方法(申请号：CN201210288373.4)，提供一种基于机器学习的并行磁共振成像GRAPPA方法；基于深度卷积网的一维部分傅里叶并行磁共振成像方法(申请号：CN201710416357.1)和一种基于深度学习和凸集投影的磁共振重建方法(申请号：CN201810306848.5)，这两个专利将深度学习方法应用于磁共振并行成像领域。目前已申请的基于U型网络的快速MRI成像方面的专利有：
一种基于双通路U型卷积神经网络的医学图像分割方法(申请号：CN201810203917.X)，主要用于医学图像分割。基于深度学习进行薄层磁共振图像重建的方法(申请号：
CN201810435306.8)，将厚层磁共振图像用三维U型网络生成薄层磁共振图像，提高了磁共振图像中灰质、白质和全部脑容积的估计准确率。目前已申请的基于递归残差神经网络的快速MRI成像方面的专利有：一种基于增强递归残差网络的磁共振图像超分辨重建方法(申请号：CN201810251558.5)，以反复递归的残差模块为基本单元建立递归残差网络，获得磁共振超分辨重建效果；一种基于递归残差注意力网络的图像超分辨率重建方法(申请号：
CN201810600397.6)，加入残差注意力网络和递归网络，增加层数的同时而不增加新的参数提高超分辨率重建的精度。目前还未能查询到任何基于递归残差U型网络的快速磁共振成像方法的授权发明专利或申请。

[0005] 目前国内外已发表的基于深度学习卷积神经网络的快速磁共振成像方面的文章有：Wang  S等人提出基于卷积神经网络进行快速磁共振图像重建(Wang  S,et al.Accelerating magnetic resonance imaging via deep learning,in Proc.IEEE 13th Int.Conf.Biomedical Imaging,pp.514–517,2016.)。Yu S等人提出基于生成对抗网络的深度学习方法加速CS‑MRI重建(Yu S,Dong H,Yang G,et al.Deep de‑aliasing for fast compressive sensing MRI.arXiv preprint arXiv:1705.07137,2017.)。Yang Y等人提出在Generic‑ADMM‑Net网络的非线性转换层中加入广义运算算子构成Complex‑ADMM‑Net实现图像重建(Yang Y,et al.ADMM‑Net:A deep learning approach for compressive sensing MRI.arXiv:1705.06869v1,2017.)。Lee D等人提出深度伪影学习网络用于CS‑MRI并行成像(Lee D,Yoo J,Ye J C.Deep artifact learning for compressed sensing and parallel MRI.arXiv preprint arXiv:1703.01120,2017.)，通过幅值网络和相位网络直接估计混叠伪影，可将欠采样重建的图像减去估计的混叠伪影来获得无混叠图像。Hammernik K等人提出一种深度变分网络并加速基于并行成像的MRI重建(Hammernik K et al.Learning a variational network for reconstruction of accelerated MRI data,Magn.Reson.Med.vol.79,no.6,pp.3055‑3071,2017.)。目前已发表的基于U型网络的MRI快速成像方面的文章有：Jin K H等人(JinK H,et al.Deep convolutional neural network for inverse problems in imaging.IEEE Transactions on Image Processing,2017,26(9):4509‑4522.)提出基于滤波反投影的深度学习网络结构解决成像中的逆问题，其中的基础网络是U型网络结构。2018年，Yang G等人提出了以U型网络为基础的生成对抗网络DAGAN，用于CS‑MRI成像(Yang G,et al.Dagan:Deep de‑aliasing generative adversarial networks for fast compressed  sensing  MRI reconstruction.IEEE Transactions on Medical Imaging,2018,37(6):1310‑1321.)。
Hyun C M等人(Hyun C M,Kim H P,Lee S M,et al.Deep learning for undersampled MRI reconstruction.Physics in medicine and biology,10.1088/1361‑6560/aac71a，
2018.)为U型网络MRI快速成像提供了数学理论依据。Corey‑Zumar将U型网络应用于脑部和前列腺图像重建(https://github.com/Corey‑Zumar/MRIReconstruction/tree/master/sub‑mrine)。目前已发表的基于递归残差U型网络的文章有：Md Z A等人使用递归残差U型网络进行视网膜图像中的血管分割、皮肤癌分割和肺病变分割(Md Z A,Mahmudul H,Chris Y.Recurrent Residual Convolutional NeuralNetwork based on U‑Net(R2U‑Net)forMedical Image Segmentation,arXiv:1802.06955,2018.)。Md Z A等人将R2U‑Net首次应用于公开数据集中的核分割，清楚地证明了R2U‑Net模型对核分割任务的鲁棒性(Md Z A,Chris Y,Tarek M T,and Vijayan K A.Nuclei Segmentation with Recurrent Residual Convolutional Neural Networks based U‑Net(R2U‑Net),IEEE National Aerospace and Electronics Conference,arXiv:1811.03447,2018.)。

[0006] 以上发表的基于深度学习卷积神经网络的快速MRI成像方面的文章或申请的发明专利，主要是基于卷积神经网络或U型网络的深度学习方法来进行快速MRI成像，或者是基于递归残差卷积神经网络的快速MRI成像，或者是基于递归残差U型网络的医学图像分割，没有出现过任何基于递归残差U型网络的快速MRI成像的专利或文章。

[0007] 本发明针对现有的基于U型网络的快速MRI成像在成像质量上的不足，对现有的U型网络模块进行改进，加入递归残差卷积神经网络模块，在增加U型网络层数的同时不增加新的参数，能在重建时恢复出更多图像的高频信息。

[0008] 本发明包含三个步骤：数据处理、模型训练、图像重建。

[0009] 1.数据处理

[0010] 数据处理包含两个步骤：数据采样、填零重建。

[0011] 1‑1.数据采样

[0012] 全采样k空间数据为kref(xk,yk)，其中xk表示k空间频率编码(Frequency Encoding,FE)方向，yk表示相位编码(Phase Encoding,PE)方向。全采样k空间数据kref(xk,yk)由全采样参考图像Iref(x,y)经离散傅里叶变换(DFT)得到。

[0013] kref(xk,yk)＝DFT(Iref(x,y))                 (1)

[0014] 对于k空间数据的模拟欠采样，使用1维均匀欠采样模板mask，在PE方向隔n行进行数据采集，在PE方向中心的4％行进行采集，在FE方向全采集。用欠采样模板和全采样k空间数据点乘表示采样过程。

[0015] ku(xk,yk)＝mask*kref(xk,yk)                 (2)

[0016] 其中，ku(xk,yk)表示欠采样k空间数据；*表示点乘；mask为全采样k空间数据相同大小的矩阵，采集的点用1表示，未采集的点用0表示。

[0017]

[0018] 1‑2.填零重建

[0019] 将得到的欠采样k空间数据进行离散傅里叶反变换，生成填零重建图像，即欠采样的图像Iu(x,y)。全采样参考图像Iref(x,y)和欠采样图像Iu(x,y)组成一对训练数据。

[0020] 2.网络训练

[0021] 网络训练包含两个步骤：构建网络、训练网络。

[0022] 2‑1.构建网络

[0023] 构建网络包含2个步骤：U型网络、递归残差模块。

[0024] 2‑1‑1.U型网络

[0025] U型网络包含降采样和升采样两部分，采样包含卷积(Convolution,Conv)、批标准化(Batch Normalization，BN)、激活函数(Activation Function,AF)、池化(Pooling)、转置卷积(Transposed Convolution)五个基本单元。

[0026] 卷积公式如下：

[0027] Cl＝Wl*Cl‑1+bl                     (4)

[0028] 其中，Cl表示第l层输出；Wl表示第l层的卷积核，其大小为wh×wh×nl，nl表示通道数，wh×wh表示单通道卷积核尺寸；Cl‑1表示第l层输入；bl表示第l层偏置量。

[0029] 批标准化公式如下：

[0030]

[0031]

[0032]

[0033]

[0034] 其中，μ，ρ分别是批数据的均值和方差；T为批数据大小；是批标准化输出；γ、β为经验参数。

[0035] 激活函数公式如下：

[0036]

[0037] 其中，σ为激活函数。

[0038] 池化公式如下：

[0039] C1＝pool(Cl)                       (10)

[0040] 其中，pool为池化。

[0041] 转置卷积公式如下：

[0042] Cl＝WlT*Cl‑1+bl                    (11)

[0043] 其中，WlT为转置卷积核。

[0044] 降采样包含4层，每层包含1个递归残差模块和1个池化模块。

[0045] 降采样的最后一层输出经过一次递归残差模块处理，得到升采样的输入。

[0046] 升采样包含4层，每层包含1个转置卷积模块和1个递归残差模块。在升采样中，每层的转置卷积输出Ca会和同层降采样递归残差计算输出Cb进行合并(concat)，得到升采样递归残差模块输入Ccat，合并公式如下：

[0047] Ccat＝concat(Ca+Cb)                   (12)

[0048] 2‑1‑2.递归残差模块

[0049] 递归残差模块包含2个递归(Recurrent)计算和1个残差(Residual)计算，递归计算由卷积模块组成每个卷积模块包含卷积、批标准化和激活函数三个单元，递归计算公式如下：

[0050]

[0051]

[0052] 其中，t表示递归次数；BN表示批标准化；Cls表示残差计算的输出。

[0053] 残差计算公式如下：

[0054] Cl‑1＝Wl‑2*Cl‑2+bl‑2                  (15)

[0055] Cls＝recus(Cl‑1)                     (16)

[0056] Clo＝Cl‑1+Cls                      (17)

[0057] 其中，recus表示递归计算；Clo表示递归残差模块输出。

[0058] 2‑2.网络训练

[0059] 网络训练包含3个步骤：损失函数、优化算法、循环计算。

[0060] 2‑2‑1.损失函数

[0061] 使用均方根误差函数(Mean Squared Error，MSE)作为反向传播的损失函数，求得输出层的损失值loss。对数据集 损失函数计算公式如下：

[0062]

[0063] 其中，角标i表示批数据的第i个矩阵；i＝(1,2…T)，RR_Unet表示递归残差U型网络；θ表示训练的网络参数。

[0064] 2‑2‑2.优化算法

[0065] 优化算法使用随机梯度下降法(Stochastic gradient descent，SGD)算法流程如下：

[0066]

[0067] θ←θ‑lr×g                     (20)

[0068] 其中，g表示第i批数据梯度均值； 表示参数梯度；lr表示学习率。

[0069] 2‑2‑3.循环计算

[0070] 以差值dif作为损失值和损失阈值的判断条件：

[0071]

[0072] 其中，τ表示损失阈值。

[0073] 计算差值dif，若差值不小于零，则执行优化算法；若差值小于零，则迭代结束。执行n次循环计算，直到求得最优网络参数θ。

[0074] 3.图像重建

[0075] 用训练好的递归残差U型网络对测试图像Itu(x,y)进行在线重建，得到重建预测图像Irecon(x,y)，重建算法如下：

[0076] Irecon(x,y)＝RR_Unet(Itu(x,y),θbest)              (22)

[0077] 其中，θbest表示训练好的网络最优参数。

[0078] 将重建预测图像进行离散傅里叶变换得到k空间预测数据krecon(xk,yk)，用实际采集到的k空间数据ktu(xk,yk)去替换krecon(xk,yk)中相应位置的数据，最后的重建图像用Iout(x,y)表示：

[0079] Iout(x,y)＝IDFT(krecon(xk,yk)(1‑mask)+ktu(xk,yk))       (23)[0080] 其中，1表示与采样mask相同的1矩阵。

[0081] 本发明采用的递归残差U型网络的快速MRI成像，在降采样和升采样的每一层都使用了递归残差模块，递归残差模块不仅能解决网络反向传播带来的梯度爆炸和梯度消失的问题，还能在增加网络层数的同时不增加新的网络参数，在训练和重建的过程中有更好的细节表现。同时采用随机梯度下降法(SGD)，能加快算法收敛速度，更快训练出较优模型。此外，采用规则欠采样方式进行数据采集，实现方式简单。本发明具有以下特点：

[0082] 1、在U型网络每个采样层加入递归残差模块，解决网络反向传播带来的梯度爆炸和梯度消失的问题，增加网络层数的同时不增加新的网络参数。

[0083] 2、在优化器选择上使用了随机梯度下降法，加快算法收敛速度，更快训练出较优模型。

[0084] 3、在采样中使用1维均匀欠采样，能加快采集速度，更易于硬件实现。

[0088] 本发明包含三个步骤：数据处理、模型训练、图像重建。

[0089] 1.数据处理

[0090] 数据处理包含两个步骤：数据采样、填零重建。

[0091] 1‑1.数据采样

[0092] 全采样k空间数据为kref(xk,yk)(如图1(a)所示)，其中xk表示k空间频率编码(Frequency Encoding,FE)方向，yk表示相位编码(Phase Encoding,PE)方向。全采样k空间数据kref(xk,yk)由全采样参考图像Iref(x,y)(如图1(b)所示)经离散傅里叶变换(DFT)得到。

[0093] kref(xk,yk)＝DFT(Iref(x,y))                  (1)

[0094] 对于k空间数据的模拟欠采样，使用1维均匀欠采样模板mask(如图1(c)所示)，在PE方向隔n行进行数据采集，在PE方向中心的4％行进行采集，在FE方向全采集。用欠采样模板和全采样k空间数据点乘表示采样过程。

[0095] ku(xk,yk)＝mask*kref(xk,yk)                (2)

[0096] 其中，ku(xk,yk)表示欠采样k空间数据(如图1(d)所示)；*表示点乘；mask为全采样k空间数据相同大小的矩阵，采集的点用1表示，未采集的点用0表示。

[0097]

[0098] 1‑2.填零重建

[0099] 将得到的欠采样k空间数据进行离散傅里叶反变换，生成填零重建图像，即欠采样的图像Iu(x,y)(如图1(e)所示)。全采样参考图像Iref(x,y)和欠采样图像Iu(x,y)组成一对训练数据。

[0100] 2.网络训练

[0101] 网络训练包含两个步骤：构建网络、训练网络。

[0102] 2‑1.构建网络

[0103] 构建网络包含2个步骤：U型网络、递归残差模块。

[0104] 2‑1‑1.U型网络

[0105] U型网络(如图2所示)包含降采样和升采样两部分，采样包含卷积(Convolution,Conv)、批标准化(Batch Normalization，BN)、激活函数(Activation Function,AF)、池化(Pooling)、转置卷积(Transposed Convolution)五个基本单元。

[0106] 卷积公式如下：

[0107] Cl＝Wl*Cl‑1+bl                       (4)

[0108] 其中，Cl表示第l层输出；Wl表示第l层的卷积核，其大小为wh×wh×nl，nl表示通道数，wh×wh表示单通道卷积核尺寸；Cl‑1表示第l层输入；bl表示第l层偏置量。

[0109] 批标准化公式如下：

[0110]

[0111]

[0112]

[0113]

[0114] 其中，μ，ρ分别是批数据的均值和方差；T为批数据大小；是批标准化输出；γ、β为经验参数。

[0115] 激活函数公式如下：

[0116]

[0117] 其中，σ为激活函数。

[0118] 池化公式如下：

[0119] Cl＝pool(Cl)                      (10)

[0120] 其中，pool为池化。

[0121] 转置卷积公式如下：

[0122] Cl＝WlT*Cl‑1+bl                    (11)

[0123] 其中，WlT为转置卷积核。

[0124] 降采样包含4层，每层包含1个递归残差模块和1个池化模块。

[0125] 降采样的最后一层输出经过一次递归残差模块处理，得到升采样的输入。

[0126] 升采样包含4层，每层包含1个转置卷积模块和1个递归残差模块。在升采样中，每层的转置卷积输出Ca会和同层降采样递归残差计算输出Cb进行合并(concat)，得到升采样递归残差模块输入Ccat，合并公式如下：

[0127] Ccat＝concat(Ca+Cb)                 (12)

[0128] 2‑1‑2.递归残差模块

[0129] 递归残差模块包含2个递归(Recurrent)计算和1个残差(Residual)计算(如图2所示)，递归计算由卷积模块组成每个卷积模块包含卷积、批标准化和激活函数三个单元，递归计算公式如下：

[0130]

[0131]

[0132] 其中，t表示递归次数；BN表示批标准化；Cls表示残差计算的输出。

[0133] 残差计算公式如下：

[0134] Cl‑1＝Wl‑2*Cl‑2+bl‑2                   (15)

[0135] Cls＝recus(Cl‑1)                     (16)

[0136] Clo＝Cl‑1+Cls                      (17)

[0137] 其中，recus表示递归计算；Clo表示递归残差模块输出。

[0138] 2‑2.网络训练

[0139] 网络训练包含3个步骤：损失函数、优化算法、循环计算。

[0140] 2‑2‑1.损失函数

[0141] 使用均方根误差函数(Mean Squared Error，MSE)作为反向传播的损失函数，求得输出层的损失值loss。对数据集 损失函数计算公式如下：

[0142]

[0143] 其中，角标i表示批数据的第i个矩阵；i＝(1,2…T)，RR_Unet表示递归残差U型网络；θ表示训练的网络参数。

[0144] 2‑2‑2.优化算法

[0145] 优化算法使用随机梯度下降法(Stochastic gradient descent，SGD)算法流程如下：

[0146]

[0147] θ←θ‑lr×g                       (20)

[0148] 其中，g表示第i批数据梯度均值； 表示参数梯度；lr表示学习率。

[0149] 2‑2‑3.循环计算

[0150] 以差值dif作为损失值和损失阈值的判断条件：

[0151]

[0152] 其中，τ表示损失阈值。

[0153] 计算差值dif，若差值不小于零，则执行优化算法；若差值小于零，则迭代结束。执行n次循环计算，直到求得最优网络参数θ。

[0154] 3.图像重建

[0155] 用训练好的递归残差U型网络对测试图像Itu(x,y)进行在线重建，得到重建预测图像Irecon(x,y)，重建算法如下：

[0156] Irecon(x,y)＝RR_Unet(Itu(x,y),θbest)              (22)

[0157] 其中，θbest表示训练好的网络最优参数。

[0158] 将重建预测图像进行离散傅里叶变换得到k空间预测数据krecon(xk,yk)，用实际采集到的k空间数据ktu(xk,yk)去替换krecon(xk,yk)中相应位置的数据，最后的重建图像用Iout(x,y)表示：

[0159] Iout(x,y)＝IDFT(krecon(xk,yk)(1‑mask)+ktu(xk,yk))         (23)[0160] 其中，1表示与采样mask相同的1矩阵。

[0161] 以下结合人体头部的MRI数据，对基于递归残差U型网络快速磁共振成像方法进行实例说明，如图1所示，(a)为全采样k空间数据kref(xk,yk)，经过离散傅里叶变换得到全采样图像，(b)为全采样图像数据Iref(x,y)，(c)为1维均匀欠采样模板mask，采样方式为每隔4行采集一行数据，并采集中间16行主要数据，数据尺度大小均为256*256，对(a)进行欠采样生成欠采样k空间数据，(d)为欠采样k空间数据ku(xk,yk)，(d)经过离散傅里叶变换得到欠采样图像Iu(x,y)，本发明采集3200对图像数据，其中3000对数据用于递归残差U型网络的训练，200对数据用于重建对比。网络结构如图2所示，包含4层降采样，4层升采样和每个采样层的一个递归残差模块。递归残差模块包含1次残差计算和2次递归计算，递归计算的递归次数为2。实验环境为INTEL I5‑4460 16G内存,NVIDIA GTX1080 8G显存,Windows10,Python3.6.5,Pytorch1.0.1。基于U型网络训练共2.5小时，重建时间约为1.2秒，基于递归残差U型网络训练共4小时，重建时间约1.6秒。图3显示了各重建方法的结果图，(a)是原始图像，(b)是填零重建图像，(c)是基于U型网络的重建图像，(d)是基于递归残差U型网络的重建图像，(e)、(f)、(g)分别是(b)、(c)、(d)与原始图像(a)的差值图。从图中可以直观的看出填零重建的误差最大，基于U型网络的重建图像(c)能较好的还原出原始图像，基于递归残差U型网络的重建图像(d)能在基于U型网络的基础上进一步提高图像的细节信息，误差更小。总相对误差(total relative error,TRE)公式如下：

[0162]

[0163] 填零重建的TRE值为1.6215，基于U型网络的重建TRE值为0.9707，而本发明的重建TRE值为0.5351。

[0164] 通过本发明可以看出，相同欠采样的条件下，在加入递归残差模块的U型网络能比U型网络在增加网络层数的同时没有增加网络参数。相比U型网络，虽然本专利提出的方法训练时间有所增加，但重建图像质量提升显著，能恢复出更多的图像细节。

[0085] 图1是采用本发明进行数据采集的示意图；

[0086] 图2是本发明的网络结构图；

[0087] 图3是图像重建实例的结果对比图；