-
电磁波的能量在不同的物质中有不同的衰减,磁共振成像(magnetic resonance imaging, MRI)利用这一特性,通过外加梯度磁场检测电磁波,从而可知构成物体原子核的位置和种类,绘制出物体内部影像[1-2]。MRI在医学成像上应用非常广泛,已经成为临床医学检查的重要手段之一,为临床诊断提供了非常有价值的信息。
但是,MRI扫描时间过长、成像较慢给病人带来了额外的痛苦,而且由于被检者的器官运动,如呼吸、吞咽、心脏自然跳动等,都会造成图像模糊和对比度失真等情况,无法满足动态实时高精度成像要求,影响MRI的进一步推广应用[3-4]。为了提高MRI成像速度,研究者们提出了多种$k$空间加速成像方法,如部分傅里叶成像、钥孔成像、并行成像等技术[5],但这些方法对成像时间的减少有限,而且会额外增加硬件成本,或者使得成像质量降低。
文献[6]在前人的基础上系统性地提出了压缩感知(compressed sensing, CS)理论。这是一个充分利用信号稀疏性或可压缩性的全新信号采集、编解码理论。基于压缩感知理论的磁共振成像可以以远低于奈奎斯特采样的频率进行采样,并从这些欠采样数据中实现图像重构。该方法可以显著降低采样量,并缩短磁共振扫描时间,提高成像速度。
依据压缩感知理论,只要信号本身是稀疏的或者在某个变换域下是稀疏的,那么就可以用一个与变换基不相关的观测矩阵将变换所得高维信号投影到一个低维空间上。然后通过求解一个最小化优化问题,从这些少量的投影中重构出原始信号。因此,压缩感知在数学上来讲就是最小${\ell _0}$范数优化问题。但直接求解该问题数值极不稳定,是一个NP-Hard问题。针对这一问题,研究者们提出了一系列寻找次优解的优化算法,如正交匹配追踪算法(orthogonal matching pursuit, OMP)[7]、迭代硬阈值法(iterative hard thresholding, IHT)[8]、共轭梯度法(conjugate gradient, CG)[9-10]、分裂Bregman迭代法(splitting bregman-based, SBB)[11]、曲波变换(curvelet transform, CT)[12]等。各种重构算法在运算复杂度、重构图像质量等方面各有优劣。
本文提出一种基于预测线搜索的共轭梯度MRI压缩成像方法。采用该方法可以减少CG算法中的线搜索次数,从而提高算法的执行效率。
-
MR信号$s({\mathit{\boldsymbol{k}}})$和MR图像$m({\mathit{\boldsymbol{r}}})$之间的关系可表示为:
$$ s({\mathit{\boldsymbol{k}}}) = \int {m({\mathit{\boldsymbol{r}}})} {{\rm{e}}^{ - 2\pi {\rm{i}}kr}}{\rm{d}}{\mathit{\boldsymbol{r}}} $$ (1) 可见,MR获取的信号$s({\mathit{\boldsymbol{k}}})$是$m({\mathit{\boldsymbol{r}}})$在空间频域的傅里叶域的数据表示,即$k$空间数据。
依据压缩感知理论,只要信号$s({\mathit{\boldsymbol{k}}})$是稀疏的,并且其测量矩阵满足有限等距特性(restricted isometry property, RIP),就可以通过求解${\ell _1}$范数的最优化问题精确重构原始的MRI图像。在压缩感知模型中,将全变分(total variation, TV)作为一个惩罚函数包含在其中[9],可以更好地保留图像的边界信息。通常$k$空间测得的数据是有噪声的,可以采用基追踪降噪法[13]对含有噪声的感知数据进行处理。这样,基于${\ell _1}$范数的CS-MRI模型如式(2)所示,然后通过求解这个无约束最小化优化问题重建原始信号$m({\mathit{\boldsymbol{r}}})$,该信号满足如下关系:
$$ \begin{array}{l} \mathop {\min }\limits_{\hat m({\mathit{\boldsymbol{r}}})} {\lambda _1}{\left\| {\hat m({\mathit{\boldsymbol{r}}})} \right\|_1} + {\lambda _2}{\left\| {\hat m({\mathit{\boldsymbol{r}}})} \right\|_{{\rm{TV}}}} + \\ \frac{1}{2}\left\| {\mathit{\Phi} \hat m({\mathit{\boldsymbol{r}}}) - \bar s({\mathit{\boldsymbol{k}}})} \right\|_2^2 \end{array} $$ (2) 式中,$\left\| \cdot \right\|$表示欧几里得范数;$\mathit{\Phi} $表示一种与$k$空间欠采样模式对应的傅里叶变换;$\bar s({\mathit{\boldsymbol{k}}})$是欠采样所得的测量数据;$\hat m({\mathit{\boldsymbol{r}}})$是根据信号$\bar s({\mathit{\boldsymbol{k}}})$恢复的被测物体的图像;${\left\| {\hat m({\mathit{\boldsymbol{r}}})} \right\|_1}$为$\hat m({\mathit{\boldsymbol{r}}})$的${\ell _1}$范数,体现图像在变换域的稀疏性;${\left\| {\hat m({\mathit{\boldsymbol{r}}})} \right\|_{{\rm{TV}}}}$为TV算子,体现图像在梯度域上的稀疏性;常数${\lambda _1}$和${\lambda _2}$控制${\ell _1}$范数和TV算子在目标函数中的权重。模型中稀疏变换采用恒等变换。
-
共轭梯度法是求解线性和非线性无约束最小化问题的有效方法之一。考虑式(2)的优化问题,令:
$$\begin{array}{l} f(\hat m({\mathit{\boldsymbol{r}}})) = \mathop \lambda \nolimits_1 {\left\| {\hat m({\mathit{\boldsymbol{r}}})} \right\|_1} + \mathop \lambda \nolimits_2 {\left\| {\hat m({\mathit{\boldsymbol{r}}})} \right\|_{{\rm{TV}}}} + \\ \frac{1}{2}\left\| {\mathit{\Phi} \hat m({\mathit{\boldsymbol{r}}}) - \bar s({\mathit{\boldsymbol{k}}})} \right\|_2^2 \end{array}$$ (3) 若$f(\hat m({\mathit{\boldsymbol{r}}}))$是连续可微的,采用共轭梯度法求解此最优化问题的迭代公式为:
$$ \hat m{({\mathit{\boldsymbol{r}}})_{k + 1}} = \hat m{({\mathit{\boldsymbol{r}}})_k} + {\alpha _k}{d_k} $$ (4) $${d_k} = \left\{ {\begin{array}{*{20}{c}} { - {g_k}}&{k = 1}\\ { - {g_k} + {\beta _k}{d_{k - 1}}}&{k \ge 2} \end{array}} \right.$$ (5) 式中,${g_k} = \nabla f(\hat m{({\mathit{\boldsymbol{r}}})_k})$为函数$f(\hat m({\mathit{\boldsymbol{r}}}))$在$\hat m{({\mathit{\boldsymbol{r}}})_k}$处的梯度;${\alpha _k}$为步长;${d_k}$为搜索方向;${\beta _k}$是一个标量,为方向调控参数。
CG方法需要求取目标函数的梯度,这就要求函数$f(\hat m({\mathit{\boldsymbol{r}}}))$连续可微。式(3)中${\ell _1}$范数等于矩阵中各元素$\hat m{({\mathit{\boldsymbol{r}}})_i}$绝对值的和,而绝对值函数是不连续的,所以在实际算法实现时,用一个连续函数$\left| {\hat m{{({\mathit{\boldsymbol{r}}})}_i}} \right| \approx $ $\sqrt {\hat m({\mathit{\boldsymbol{r}}})_i^* * \hat m{{({\mathit{\boldsymbol{r}}})}_i} + \mu } $代替[9],$\mu \in [{10^{ - 15}},{10^{ - 6}}]$,有:
$$ \frac{{d\left| {\hat m{{({\mathit{\boldsymbol{r}}})}_i}} \right|}}{{d\hat m{{({\mathit{\boldsymbol{r}}})}_i}}} \approx \frac{{\hat m{{({\mathit{\boldsymbol{r}}})}_i}}}{{\sqrt {\hat m({\mathit{\boldsymbol{r}}})_i^* * \hat m{{({\mathit{\boldsymbol{r}}})}_i} + \mu } }} $$ (6) TV算子是图像梯度值的${\ell _1}$范数,在计算时也作类似处理。记W1和W2为两个对角矩阵,且其对角线元素分别为${({W_1})_i} = \sqrt {\hat m({\mathit{\boldsymbol{r}}})_i^* * \hat m{{({\mathit{\boldsymbol{r}}})}_i} + \mu } $和${({W_2})_i} = \sqrt {(\Delta \hat m({\mathit{\boldsymbol{r}}}))_i^* * (\Delta \hat m{{({\mathit{\boldsymbol{r}}})}_i}) + \mu } $,其中${( \cdot )^{\rm{*}}}$表示共轭,$\Delta \hat m({\mathit{\boldsymbol{r}}})$表示图像梯度值[14]。这样,目标函数$f(\hat m({\mathit{\boldsymbol{r}}}))$的梯度可表示为:
$$ \begin{array}{l} \nabla f(\hat m({\mathit{\boldsymbol{r}}})) \approx {\lambda _1}W_1^{ - 1}\hat m({\mathit{\boldsymbol{r}}}) + {\lambda _2}{\Delta ^*}W_2^{ - 1}\Delta \hat m({\mathit{\boldsymbol{r}}}) + \\ {\mathit{\Phi} ^*}(\mathit{\Phi} \hat m({\mathit{\boldsymbol{r}}}) - \bar s({\mathit{\boldsymbol{k}}})) \end{array} $$ (7) -
共轭梯度法中${\beta _k}$控制着搜索方向,不同的${\beta _k}$对应不同的共轭梯度法[15]。文献[9-10]采用的${\beta _k}$如下所示:
$$ \beta _k^{{\rm{FR}}} = \frac{{{{\left\| {{g_k}} \right\|}^2}}}{{{{\left\| {{g_{k - 1}}} \right\|}^2}}} $$ (8) 基于式(8)的共轭梯度法常称为FR-CG法[16]。FR-CG方法运算简单,也有较好的数值效果,但为了保证算法的收敛性,需要步长${\alpha _k}$满足强Wolfe条件[17]。后来一些研究者又对这些公式进行了改进,提出了一些改进方法[18]和新的方法[19],这些方法改善了原有方法的依赖条件,也有较好的数值效果和收敛性,但增加了算法的复杂度。经过对多种CG方法的分析比较,本文采用DY-CG方法[17]来实现MRI图像的重构,其方向调控参数为:
$$\beta _k^{{\rm{DY}}} = \frac{{{{\left\| {{g_k}} \right\|}^2}}}{{{\mathit{\boldsymbol{d}}}_{k - 1}^{\rm{T}}({g_k} - {g_{k - 1}})}} $$ (9) 式中不仅包含了目标函数的梯度值,还包含了上一次迭代方向${{\mathit{\boldsymbol{d}}}_{k - 1}}$,使得该方法具有更好的性能,它只需要步长${\alpha _k}$满足标准Wolfe线搜索条件就具有全局收敛性。
-
共轭梯度法中的步长${\alpha _k}$一般通过线搜索确定。文献[9]采用了回溯线搜索(backtracking line search, BLS)方法,该方法思路为:先确定一个足够大的初始值,然后从该初始值开始逐渐减小,直到满足Wolfe条件,即:
$$ \begin{array}{l} f(\hat m{({\mathit{\boldsymbol{r}}})_k} + {\alpha _k}{{\mathit{\boldsymbol{d}}}_k}) \le f(\hat m{({\mathit{\boldsymbol{r}}})_k}) + {c_1}{\alpha _k}{\mathit{\boldsymbol{d}}}_k^{\rm{T}}\nabla f(\hat m{({\mathit{\boldsymbol{r}}})_k})\\ \;\;\;\;\;{\mathit{\boldsymbol{d}}}_k^{\rm{T}}\nabla f(\hat m{({\mathit{\boldsymbol{r}}})_k} + {\alpha _k}{{\mathit{\boldsymbol{d}}}_k}) \ge {c_2}{\mathit{\boldsymbol{d}}}_k^{\rm{T}}\nabla f(\hat m{({\mathit{\boldsymbol{r}}})_k}) \end{array} $$ 式中,${c_1}$、${c_2}$为常数,且${\rm{0}} < {c_1} < {c_2} < 1$。
CG方法求解最优化问题是一个迭代过程,每一次迭代都需要进行一次线搜索来获取合适的步长,且随着迭代次数的增加,步长值呈下降趋势。如果每次迭代中都将线搜索初始值设置为一个相同的最大值,会增加线搜索的时间;如果选择太小,可能得不到合适的步长值。因此,在迭代过程中需要调整线搜索初始步长的大小。
文献[9]设计一个系数$\beta $来修正步长初始值,即:
$$ {\alpha _{(k + 1)}}_0 = \left\{ {\begin{array}{*{20}{l}} {{\alpha _{(k + 1)}}_0\beta }&{i > 2}\\ {{\alpha _{k0}}}&{i = 1,2}\\ {{\alpha _{(k + 1)}}_0/\beta }&{i = 0} \end{array}} \right. $$ (10) 式中,$0 < \beta < 1$;$i$为第$k$次迭代时的线搜索次数, $i \geqslant 0$;${\alpha _{k0}}$和${\alpha _{(k{\rm{ + }}1)0}}$分别为第$k$次迭代和第k+1次迭代中线搜索初始值。每次迭代开始时都设置$i = 0$,若${\alpha _{k0}}$满足Wolfe条件,则$i$就保持为0。
但这种步长初值选择的方法只简单使用了一个系数,在确定第$k{\rm{ + }}1$次线搜索初始值的时候只利用了第$k$次迭代的线搜索初始值,而不管上一次线搜索后得到实际值。这样,不管第$k$次迭代线搜索初值和实际值之间的差值是多少,在进行第$k{\rm{ + }}1$次迭代时线搜索初始值缩小或放大的倍数不变。而这个差值越大,意味着搜索次数会越多,搜索所消耗的时间越长。因此,可以考虑利用这个差值,来预测第$k{\rm{ + }}1$次迭代时线搜索的初始步长值,如式(11)所示:
$$ {\alpha _{(k{\rm{ + }}1)0}} = {\alpha _{k0}} + \beta ({\alpha _k} - {\alpha _{k0}}) $$ (11) 式中,$0 < \beta < 1$;${\alpha _k}$为第$k$次迭代最后采用的实际步长值。从式(11)可以看出,若${\alpha _k} - {\alpha _{k0}}$越大,则${\alpha _{(k{\rm{ + }}1)0}}$相对于${\alpha _{k0}}$的变化也越大。这种初值确定法实际上用上一次迭代过程中的步长初始值和实际值来预测下一次迭代时采用的线搜索初始值,称这种方法为预测线搜索法(prediction line search, PLS),基于PLS的DY-CG算法具体步骤可描述为:
1) 设置最大线搜索次数M,最大迭代次数N,参数${\lambda _1}$、${\lambda _2}$、${c_1}$、${c_2}$和$\beta $的值;
2) 设置$\hat m{({\mathit{\boldsymbol{r}}})_1}{\rm{ = }}{\mathit{\Phi} ^{\rm{H}}}\bar s({\mathit{\boldsymbol{k}}})$,${g_1} = \nabla f(\hat m{({\mathit{\boldsymbol{r}}})_1})$,如果${g_1} = 0$,返回$\hat m{({\mathit{\boldsymbol{r}}})_1}$,结束;
3) 设置$k = 1$,${d_1} = - {g_1}$,${\alpha _k} = {\alpha _{k0}} = 1$,$i = 0$;
4) 计算Wolfe条件,如果成立则执行步骤6);
5) 如果$i < M$,则$i = i + 1$,${\alpha _k} = {\alpha _k}\beta $,转到步骤4),否则给出错误提示,结束;
6) 计算$\hat m{({\mathit{\boldsymbol{r}}})_{k + 1}} = \hat m{({\mathit{\boldsymbol{r}}})_k} + {\alpha _k}{{\mathit{\boldsymbol{d}}}_k}$,${g_{k + 1}} = $ $\nabla f(\hat m{({\mathit{\boldsymbol{r}}})_{k + 1}})$,如果${g_{k + 1}} = 0$或$k \geqslant N$,返回$\hat m{({\mathit{\boldsymbol{r}}})_{k + 1}}$,结束;
7) 用式(8)计算$\beta _{k + 1}^{{\rm{DY}}}$,用式(5)计算${{\mathit{\boldsymbol{d}}}_{k + 1}}$,${\alpha _{k + 1}} = {\alpha _{(k + 1)0}} = {\alpha _{k0}} + \beta ({\alpha _k} - {\alpha _{k0}})$,$i = 0$,$k = k + 1$,转到步骤4)。
-
为了证明前述方法的有效性,本文采用Shepp- Logan图[20]作为试验对象进行随机欠采样与重构分析。Shepp-Logan图像由10个椭圆构成人的大脑模拟图,其原始图像如1a所示。分析时采用几个不同的采样率对数据进行欠采样,然后将数据进行分析,并将分析结果与文献[9]中零填充法和文献[10]中基于FR-CG的压缩感知方法所得结果进行比较,以分析算法的有效性。
实验中采用MATLAB的phantom函数生成512×512的Shepp-Logan图像,然后将该图像做二维FFT变换,对得到的$k$空间数据进行欠采样。由于$k$空间数据对重构所得图像的贡献是不均匀的,中心部分主要提供了图像的能量,边缘部分则提供了图像的细节[9]。所以在欠采样时采用变密度随机欠采样方法进行采样。图 1b是采样率为0.3时的采样矩阵图,可见中心采样密度高,周围采样密度低。
重构图像的质量采用重构图像与原始图像的结构相似性指标SSIM[21] (structural similarity)进行衡量。两个图像${I_1}$和${I_2}$的SSIM定义式为:
$$ {\rm{SSIM}}({I_1},{I_2}) = \frac{{(2{\mu _1}{\mu _2} + {c_1})(2{\sigma _{12}} + {c_2})}}{{(\mu _1^2{\rm{ + }}\mu _2^2 + {c_1})(\sigma _1^2 + \sigma _2^2 + {c_2})}}$$ (12) 式中,${\mu _1}$、${\mu _2}$分别为图像${I_1}$和${I_2}$的平均值;${\sigma _1}$、${\sigma _2}$分别为图像${I_1}$和${I_2}$的方差;${\sigma _{12}}$为图像${I_1}$和${I_2}$的协方差。两幅一模一样的图像,其SSIM的值等于1。计算重构图像和原始图像的SSIM,所得值越大就表明重构图像的质量越好。
本实验中的所有算法均在MATLAB 2016a中完成,计算机CPU为Intel i5 2.5 GHz,内存为4 GB。
-
实验中分别用0.1、0.2和0.3的采样率进行采样,重构操作中${\lambda _1}$取值0.01,${\lambda _2}$取值0.05,最大迭代次数设置为25次,最大线搜索次数为150次,PLS法中的修正因子$\beta $取值0.7。采样率为0.3时图像重构过程中线搜索参量如图 2所示,最后的重构图像如图 3所示。
图 2显示的是采样率为0.3时两种线搜索参数方法的实际参量,包括线搜索步长初始值、实际值和每次从步长初始值经过搜索得到实际步长值的线搜索次数。图 2a显示的是采用BLS方法的参数值,图 2b是采用PLS方法所得到参量值。从图中可以看到两种线搜索方法在每次迭代时使用的步长几乎一样,但BLS方法的初始值选取比实际值大得多,所以需要经过多次迭代才能得到合适的步长值,这样所用的搜索时间也较长。而PLS方法采用的线搜索初始值更接近最终的实际值,所以线搜索次数少,甚至初值就是最后的实际值,如图 2b所示的第5次和第6次迭代。表 1显示了两种算法执行时间,可见在几种采样率下,基于PLS法的图像重构时间均少于基于BLS法的重构时间。
表 1 重构算法执行时间
s 采样率 0.3 0.2 0.1 BLS法 5.942 7 5.903 0 5.884 2 PLS法 5.101 2 5.079 5 5.115 5 图 3中3行图像是在采样率分别为0.3、0.2和0.1下进行采样得到的。每行包括了3个图像,第一个为零填充后重构得到的图像,第二个为采用FR-CG方法重构所得的图像,第三个为采用DY-CG方法重构所得的图像。
从图中可以看出,随着采样率的降低,采样的数据中保留原始图像的信息也就越少,重构图像3a、3d、3g的SSIM值依次减少,图像质量依次降低。同时,也可以从图中直观地看到,采用DY-CG方法重构的图像比采用FR-CG方法重构的图像质量高。
An Improved Conjugate Gradient Algorithm for Compressed MRI Imaging
-
摘要: 为了缩短磁共振成像系统的扫描时间,压缩感知方法利用欠采样数据和非线性恢复算法实现系统的实时或准实时成像需求。通过联合考虑MRI图像在变换域和梯度域下的稀疏性,提出了一种基于预测线搜索方法的共轭梯度算法来重建磁共振图像。针对共轭梯度算法中线搜索次数过多和运行时间过长问题,采用基于预测的方法来优化搜索步长值,以此缩短算法执行时间和减少线搜索次数。仿真实验利用磁共振图像的10%、20%和30%的下采样数据进行图像重建,结果显示基于该预测线搜索方法的压缩成像算法执行时间少于回溯线搜索法的执行时间,重构图像质量优于零填充法和FR共轭梯度法,验证了该算法的有效性。Abstract: In order to shorten the scanning time of magnetic resonance imaging (MRI), a combination of sub-nyquist sampling data and nonlinear reconstruction algorithms is adopted in the compressed sensing method to realize the real-time or quasi real-time imaging requirement. Considering the sparseness of MRI images in the transform domain and the gradient domain, a conjugate gradient algorithm is proposed to reconstruct the magnetic resonance image using the prediction line search method. The conventional conjugate gradient algorithm suffers a long computation time because of too many attempts in the line search process. We address this problem by optimizing the searching step size with a prediction approach, which significantly reduces the line search cost and accelerates convergence. Simulation results, with under-sampling rate at 10%, 20% and 30%, show that the prediction line search method achieves the best image reconstruction resolution in comparison to the zero filling method and the FR conjugate gradient method, while its time cost is less than the backtracking line search method, which illustrate the effectiveness of the proposed algorithm.
-
表 1 重构算法执行时间
s 采样率 0.3 0.2 0.1 BLS法 5.942 7 5.903 0 5.884 2 PLS法 5.101 2 5.079 5 5.115 5 -
[1] BROOKES M J, VRBA J, MULLINGER K J, et al. Source localisation in concurrent EEG/fMRI:Applications at 7T[J]. Neuroimage, 2009, 45(2):440-452. doi: 10.1016/j.neuroimage.2008.10.047 [2] ZHANG Y, PETERSON B, DONG Z. A support-based reconstruction for Sense MRI[J]. Sensors, 2013, 13(4):4029-4040. doi: 10.3390/s130404029 [3] 王水花, 张煜东.压缩感知磁共振成像技术综述[J].中国医学物理学杂志, 2015, 32(2):158-162. doi: 10.3969/j.issn.1005-202X.2015.02.002 WANG Shui-hua, ZHANG Yu-dong. Survey on compressed sensing magnetic resonance imaging technique[J]. Chinese Journal of Medical Physics, 2015, 32(2):158-162. doi: 10.3969/j.issn.1005-202X.2015.02.002 [4] RUNGE V M, RICHTER J K, HEVERHAGEN J T. Speed in clinical magnetic resonance[J]. Investigative Radiology, 2017, 52(1):1-17 http://d.old.wanfangdata.com.cn/OAPaper/oai_pubmedcentral.nih.gov_3381802 [5] 翁卓, 谢国喜, 刘新, 等.基于k空间加速采集的磁共振成像技术[J].中国生物医学工程学报, 2010, 29(5):785-790. doi: 10.3969/j.issn.0258-8021.2010.05.023 WENG Zhuo, XIE Guo-xi, LIU Xin, et al. Development of fast magnetic resonance imaging techniques based on k-space accelerated collection[J]. Chinese Journal of Biomedical Engineering, 2010, 29(5):785-790. doi: 10.3969/j.issn.0258-8021.2010.05.023 [6] DAVID L D. Compressed sensing[J]. IEEE Transactions on Information Theory, 2006, 52(4):1289-1306. doi: 10.1109/TIT.2006.871582 [7] SOUSSEN C, GRIBONVAL R, IDIER J, et al. Joint K-step analysis of orthogonal matching pursuit and orthogonal least squares[J]. IEEE TransInform Theory, 2013, 59(5):3158-3174. doi: 10.1109/TIT.2013.2238606 [8] BLUMENSATH T, DAVIES M E. Iterative hard thresholding for compressed sensing[J]. Applied and Computational Harmonic Analysis, 2009, 27(3):265-274. http://d.old.wanfangdata.com.cn/OAPaper/oai_arXiv.org_0805.0510 [9] LUSTIG M, DONOHO D, PAULY J M. Sparse MRI:the application of compressed sensing for rapid MR imaging[J]. Magnetic Resonance in Medicine, 2007, 58(6):1182-1195. doi: 10.1002/(ISSN)1522-2594 [10] YANG Xia-han, ZHENG Y R, Yang M, et al. Compressed sensing with non-uniform fast fourier transform for radial Ultra-short Echo Time (UTE) MRI[C]//2015 IEEE 12th International Symposium on Biomedical Imaging (ISBI). New York, USA: IEEE, 2015, 1: 918-921. [11] BI Dong-jie, MA Lan, XIE X, et al. A splitting bregman-based compressed sensing approach for radial UTE MRI[J]. IEEE Transaction on Applied Superconductivity, 2016, 26(7):1-5. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=9a00cbcd530d1ceee3a6e55d933cbd1c [12] JAING Xiao-ping, DING Hao, ZHANG Hua, et al. Study on compressed sensing reconstruction algorithm of medical image based on curvelet transform of image block[J]. Neurocomputing, 2017, 220:191-198. doi: 10.1016/j.neucom.2016.04.062 [13] CHEN S S, DAVID L D, SAUNDERS M A. Atomic decomposition by basis pursuit[J]. SIAM Journal on Scientific Computing, 2001, 43(1):129-159. doi: 10.1137-S1064827596304010/ [14] RUDIN L, OSHER S, FATEMI E. Non-linear total variation noise removal algorithm[J]. Physica D, 1992, 60:259-268. doi: 10.1016/0167-2789(92)90242-F [15] WILLIAN W H, ZHANG Hong-chao. A survery of nonlinear conjugate gradient methods[J]. Pacific Journal of Optimization, 2006, 2(1):35-58. [16] FLETCHER R, REEVES C. Function minimization by conjugate gradients[J]. Computer Journal, 1964, 7(2):149-154. doi: 10.1093-comjnl-7.2.149/ [17] DAI Y H, YUAN Y. A nonlinear conjugate gradient method with a strong global convergence property[J]. SIAM Journal on Optimization, 1999, 10(1):177-182. doi: 10.1137/S1052623497318992 [18] LU Ai-guo, LIU Hong-wei, ZHENG X, et al. A variant spectral-type FR conjugate gradient method and its global convergence[J]. Applied Mathematics and Computation, 2011, 217(12):5547-5552. doi: 10.1016/j.amc.2010.12.028 [19] FATEMI M. A new efficient conjugate gradient method for unconstrained optiomization[J]. Journal of Computational & Applied Mathematics, 2016, 300:207-216. doi: 10.1023/A%3A1012930416777 [20] SHEPP L A, LOGAN B F. The fourier reconstruction of a head section[J]. IEEE Transactions on Nuclear Science, 1974, 21(3):21-43. doi: 10.1109/TNS.1974.6499235 [21] WANG Z, BOVIK A C, SHEIKH H R, et al. Image quality assessment:from error visibility to structural similarity[J]. IEEE Transactions on Image Process, 2004, 13(4):600-612. doi: 10.1089-fpd.2009.0394/