磁共振成像已经成为医学领域里非常具有影响力的无创诊断成像技术之一。然而，成像时间长一直制约着其进一步发展，通过物理手段提高成像速度已经达到极限^[1]。在不影响成像质量的前提下，减少采样时间，通过设计更好的重建算法来提高成像速度为研究者们提供了新思路。如何利用少量采样数据重建出高质量的磁共振图像成为近年来的热点问题。压缩传感^[2-5]具有可以从少量随机采样数据中恢复信号的性质，因此，可以应用于磁共振图像的重建问题，进而发展成为磁共振成像领域里主流的成像技术，如Lustig提出的压缩传感磁共振成像方法(compressed sensing MRI，CS-MRI^[6-8])。压缩传感要求待重建的信号在稀疏字典或稀疏变换(统称为稀疏字典)下具有稀疏表示，目前的稀疏字典主要分为2类：

第一类是全局稀疏字典，常用的有：小波变换^[6-8]、曲线小波变换^[9]，奇异值分解^[10]等。这类预定义字典能对图像的整体进行稀疏表示，然而，使用全局稀疏字典会使重建的图像丢失图像边缘^[7-11]等的一些比较精细的特征，如图 1(b)所示。该图是将小波变换作为全局稀疏字典的CS-MRI^[8]的重建结果，从图中可以看出边缘信息的缺失。

第二类是局部稀疏字典，这类字典以图像块为基本单位，通过学习得到的局部稀疏字典来求得所有图像块的稀疏表示^[12-21]。以局部稀疏表示为约束，重建的图像具有比较精细的图像细节，然而常常会丢失图像的整体结构信息，常见的方法有PBDW^[17]，PANO^[18]和DLMRI^[19]。图 1(c)是DLMRI在图像块的大小为10×10且块与块之间不重叠的情况下所得到的重建结果，从图中可以看出图像块之间不平滑的块状效应。

然而，独立使用全局或局部稀疏字典，会分别导致图像细节或图像整体结构信息的丢失。笔者提出一个联合利用局部和全局稀疏字典来表示图像稀疏结构的成像模型(global and local sparse MRI，GLSMRI)，既保证了图像的整体结构信息，也较好地保留了比较精细的图像细节，克服了单独使用全局或局部稀疏字典所带来的问题。该模型由2个子模型构成：1)局部稀疏模型，它通过学习一个局部稀疏字典来对图像块进行稀疏表示；2)全局稀疏模型，它通过预定义的稀疏字典来对图像进行稀疏表示。研究的GLSMRI在传统的压缩传感框架下进行求解，并且求解过程以迭代的方式进行，来逐步获得更好的重建结果。笔者使用大脑数据集对GLSMRI进行验证，并对模型中的几个重要参数进行讨论。

1 压缩感知磁共振成像方法回顾

设X代表待重建的2D磁共振图像，Y=F_uX。代表了k空间的传感过程，其中F_u是部分傅立叶变换，Y是采样的k空间数据。假设ψ代表全局稀疏字典或者从图像块中学习得到的局部稀疏字典，则X的稀疏表示可以定义为：α=ψX。基于压缩感知来从采样数据Y中重构磁共振图像的成像模型如下^[20]

$ {P_0}:\min {\left\| {\psi X} \right\|_1}{\text{s}}.{\text{t}}.\;\left\| {Y-{F_u}X} \right\|_2^2 \leqslant \varepsilon $

式中：ε代表了重建过程中允许的误差向量。在模型中l₁范数加强了ψX的稀疏性，误差约束保证了对k空间采样数据的保真度。

当ψ为小波变换(DWT)时，成像模型就变成了基于全局稀疏字典的CS-MRI算法。此外，当稀疏字典是利用字典学习从图像块中学习出来的时候，成像模型就变成了基于局部稀疏字典的DLMRI算法，成像模型如下

$ \begin{gathered} \mathop {\min }\limits_{X.D.\Lambda = \left( {aij\;|\forall i, j} \right)} \left\{ {\sum\limits_{ij} {\left\| {{R_{ij}}X-{D_{aij}}} \right\|_2^2 + \lambda \left\| {Y-{F_u}X} \right\|_2^2} } \right\}, \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\text{s}}.{\text{t}}.\;{\left\| {{a_{ij}}} \right\|_D} \leqslant {T_0}\;\forall \;i, j, \hfill \\ \end{gathered} $

式中：R_ij是图像块提取操作；D是局部稀疏字典；α_ij是每个图像块R_ijX在D上的稀疏表示; T₀代表了稀疏水平。

2 联合局部和全局稀疏表示的磁共振图像重建方法 2.1 成像模型

当前的磁共振成像方法单独利用图像的全局或者局部稀疏结构来进行磁共振图像的重构，常常导致丢失图像的精细特征或整体结构。GLSMRI联合利用了图像的全局稀疏表示和局部稀疏表示来对磁共振图像进行重构，其成像模型如下

$ \begin{gathered} {P_{0\;:\;}}\mathop {\min }\limits_{X.D.\Lambda = \left( {aij\;|\forall i, j} \right)} \left\{ {\lambda \left\| {Y-{F_u}X} \right\|_2^2 + {\lambda _L}\sum\limits_{ij} {\left\| {{R_{ij}}X-{D_{aij}}} \right\|_2^2 + {\lambda _G}{{\left\| \Psi \right\|}_1}} } \right\}, \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\text{s}}.{\text{t}}.\;{\left\| {{a_{ij}}} \right\|_D} \leqslant {T_0}\;\forall \;i, j, \hfill \\ \end{gathered} $

(1)

目标函数中第一项是k空间数据保真项，其主要作用是使得重建图像的k空间数据与采样得到的k空间数据尽可能接近；第二项是局部稀疏项；第三项是全局稀疏约束，其中ψ是全局稀疏字典。参数λ_L和λ_C用于平衡图像的局部和全局稀疏约束。显然，当λ_L=0或者λ_C=0时, GLSMRI分别退化成DLMRI或CS-MRI。

2.2 求解算法

下面对问题P₀进行求解，其求解过程可以分成2个步骤。首先，学习局部稀疏字典D，然后获得每个图像块的稀疏表示α_ij。其次，联合局部和全局稀疏约束，在压缩传感框架重建图像，详细描述如下。

2.2.1 图像的局部稀疏表示

P₀中图像的局部稀疏表示模型如下

$ \begin{gathered} {P_{1\;:\;}}\mathop {\min }\limits_{D.\Lambda = \left( {aij\;|\forall i, j} \right)} \sum\limits_{ij} {\left\| {{R_{ij}}X-{D_{{\alpha _{ij}}}}} \right\|_2^2}, \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\text{s}}.{\text{t}}.\;{\left\| {{a_{ij}}} \right\|_D} \leqslant {T_0}\;\forall \;i, j, \hfill \\ \end{gathered} $

(2)

笔者使用过完备字典训练算法K-SVD^[12]从初始图像中训练出局部稀疏字典D^*，当D^*确定后，利用正交匹配算法OMP^[22-23]算法确定每个图像块的稀疏表示α_ij^*。

2.2.2 局部和全局稀疏约束下的图像重建

在D^*和α_ij^*已知的情况下，GLSMRI的重建子问题可以由下面公式来描述

$ \begin{gathered} {P_{2\;}}:\mathop {\min }\limits_X f\left( x \right) = \left\| {Y-{F_u}X} \right\|_2^2 + {\lambda _L} \hfill \\ \sum\limits_{ij} {\left\| {{R_{ij}}X-{D^*}\alpha _{ij}^*} \right\|_2^2 + {\lambda _G}\left\| {\psi X} \right\|} \hfill \\ \end{gathered} 。 $

(3)

子问题P₂可以用非线性共轭梯度算法来求解^[24]。f(X)的梯度为

$ \begin{gathered} \nabla f\left( X \right) = 2{\lambda _L}\sum\limits_{ij} {\left( {R_{ij}^T{R_{ij}}X-R_{ij}^T{D^*}\alpha _{ij}^*} \right)} \hfill \\ + 2F_u^H\left( {{F_u}X-Y} \right) + {\lambda _G}\nabla \left( {{{\left\| {\psi X} \right\|}_1}} \right) \hfill \\ \end{gathered} 。 $

(4)

进一步化简为

$ \begin{gathered} \nabla f\left( X \right) = 2{\lambda _L}\beta \left( {X-{X^*}} \right) + 2F_u^H \hfill \\ \left( {{F_u}X-Y} \right) + {\lambda _G}\nabla \left( {{{\left\| {\psi X} \right\|}_1}} \right), \hfill \\ \end{gathered} $

(5)

其中

$ {X^*} = \frac{{\sum\limits_{ij} {R_{ij}^T{D^*}\alpha _{ij}^*} }}{c}, $

(6)

式中：X^*是由图像的局部稀疏表示α_ij^*构成的中间图像，具有较精细的图像细节，但图像块之间会出现块状效应。式6中c代表图像块对像素点的贡献度，可以通过c=n/r²进行计算，n是一个图像块中像素点的个数，r代表 2个相邻的图像块之间的重叠距离。

这2个步骤可以重复迭代进行来逐步改善重建图像的质量，整个算法的流程如图 2所示。

与DLMRI相比，GLSMRI做了2方面的改进

1) 在成像模型方面，DLMRI没有全局稀疏约束项，即λ_G=0其优化方程为2λ_Lβ(X-X^*)+2F_u^H(F_uX-Y)=0，并用最小二乘方法求解该方程。由于x^*是对重建图像的粗略估计，而最小二乘的解常常过拟合，导致不能很好地消除X^*中存在的块状效应。GLSMRI通过加入全局稀疏约束来抑制过拟合，从而有效地消除了X^*中的块状效应。

2) 在求解方法方面，DLMRI的最小二乘解等价于简单的k空间回填：对于采样点(i, j)，用X^*对应的k空间数据S=F^-1X^*和采样点数据Y的加权平均进行回填, 即{S(i, j)+vY(i, j)}/(1+v), 其中v=1/λ_L。对于未采样点(i, j)，简单地用S(i, j)进行回填。GLSMRI增加了全局稀疏约束，并使用非线性共轭梯度算法进行迭代求解，避免了简单k空间回填导致的过拟合，从而可以得到更好的重建结果。

为了验证观点，设计了如图 3所示的对比实验。相比于图 3(c)所示的DLMRI重建图像，图 3(b)中GLSMRI的重建图像无明显的块状效应，并且比图 3(a)中CS-MRI的重建图像有更精细的图像细节。

3 实验结果及讨论

实验数据是大小为512×512的大脑图像^[25]。采用了基于高斯随机分布的2种采样模式：沿相位编码方向的笛卡尔随机采样，如图 4(a)所示；二维随机采样，如图 4(b)所示。这2种采样模式在低频部分具有相对密集的采样数据。K-SVD算法的代码从^[19]获得，为了满足算法的求解需要，修改了CS-MRI中的非线性共轭梯度求解工具。所有的模拟和重建实验在MATLAB R2012b (math works, Natick, MA)下进行。

实验中，用零填充图像作为初始图像以获取初始局部稀疏字典。GLSMRI的几个重要参数设置如下：图像块的大小设定为6×6，字典大小设定为36，重叠距离设定为1，训练字典所需的图像块个数设定为5 000，用db4小波变换做为全局稀疏字典。为了对算法进行评价，将DLMRI作为对比方法并使用峰值信噪比(PSNR)来量化重建图像的质量，此外，还展示了GLSMRI中的全局稀疏约束是如何影响图像重建的。对于大脑图像，λ_L=1/300，λ_L/λ_G=50，即局部稀疏权重是全局稀疏权重的50倍。在字典学习阶段，从X中随机选择一部分图像块，然后从中训练得到D，此外，采用周期延拓的方式来处理边界的图像块^[19]。

在图 5中，尽管从图 5(b-c)看不出重建图像之间有明细的视觉差别，但是可以从图 5(g)看出GLSMRI在重建质量上超过了DLMRI，尤其在用箭头标注的图像纹理变化丰富的位置。图 5(e-f)揭示了GLSMRI相比于DLMRI可以恢复更多的图像细节，并且图 5(h)表明对于不同的下采样因子，GLSMRI仍然优于DLMRI。当采用二维随机采样时，两个方法的差别更为明显，从图 6(g)可以看出，GLSMRI可以更接近真实图像。图 5(d)和图 6(d)表明GLSMRI相比于DLMRI需要更少的迭代次数，具体表现为：GLSMRI经过6次迭代稳定下来，DLMRI经过11次迭代稳定下来。最后，从PSNR来看，GLSMRI相对于DLMRI有1.8 dB的提升，从误差图 6(e-f)来看，DLMRI的误差包含更多的结构性的信息。

在参数选择实验中，验证了GLSMRI中几个重要的参数以及它们的选择，具体策略为：改变其中一个参数的同时固定其余参数。对于实验数据，仍使用大脑图像，并将下采样因子设置成4。

对于重叠距离r，图 7(a)表明，当r增加时，PSNR随之下降，因此将r设置为1。为了选择λ_L与λ_G的最优组合，首先将λ_L固定为1/300，图 7(b)显示了PSNR随权衡局部稀疏和全局稀疏的比值λ_L/λ_G的变化曲线。该图显示：随着比值λ_L/λ_G从0增加到40，PSNR逐渐增大，当λ_L/λ_G从40继续增加时，曲线趋于稳定，当λ_L/λ_G趋于无穷时，PSNR逐渐下降。其中，λ_L/λ_G为0和λ_L/λ_G为无穷大这2种极端情况分别代表了CS-MRI和DLMRI。结果表明，λ_L/λ_G=50是一个合理的选择。最后，固定λ_L/λ_G=50，图 7(c)的PSNR随λ_L的变化曲线表明λ_L=1/300是一个合适的选择。

4 结语

笔者提出了一个新颖的联合利用局部稀疏和全局稀疏约束来捕捉磁共振图像的局部和整体结构信息的重建模型。分2步来解决该模型：首先，学习局部稀疏字典，然后对图像块进行稀疏编码；其次，联合利用局部稀疏和全局稀疏约束来进行重建。模拟实验表明，GLSMRI超过了现有的仅利用局部稀疏结构或全局稀疏结构的方法。几个重要参数的选取原理也被阐述。将来的工作将着眼于：开发更为鲁棒的局部稀疏字典和使用自适应的全局稀疏字典来改善重建质量。