尹宝才 郭晓明 施云惠 丁文鹏

（北京工业大学城市交通学院，多媒体与智能软件技术北京市重点实验室，北京，100124）

引 言

目前研究表明，压缩感知理论可以在图像少量观测值下完成对该图像的重建［1-3］。通常，图像具有块结构特性并且可以在确定小波基下稀疏表示［4，5］。TV-Wavelet-L1（TVWL1）模型采用全变分正则化约束（Total-variation，TV）和小波L1范式正则化约束的组合稀疏表示方式，使得基于压缩感知理论的图像重建算法精确度得到提高［6，7］。

对应于图像的分解和重建过程，图像的稀疏表示方法分为分析模型和综合模型。例如：综合的稀疏表示方法是一个图像x∈Rd可以描述为x＝Dα，这里D∈Rd×n是变换矩阵，α是x的稀疏表示系数，‖α‖0＝k≪n。而分析稀疏表示方法是向量Ωx可以被稀疏化，Ω是一个线性算子，相对稀疏性可以定义为向量Ωx的零元素的数量。

当综合模型被广泛学习研究时，分析模型的稀疏表示方法却很少被关注［8］。但是，最近几年，分析稀疏表示方法得到了人们的广泛关注。基于TVWL1模型，相较于分析方式，综合稀疏表示方式的TV正则化约束更为复杂并且求解过程更加困难。所以，通常应用基于分析稀疏表示方式的TVWL1模型［9］来进行图像重建工作，该模型如下

对于简单正则化约束，例如单独的TV或L1范数约束，迭代收缩阈值算法（Iterative shrinkagethresholding algorithm，ISTA）和快速迭代收缩阈值算法都可以进行高效的图像重建［10］。但是它们不适用于TV和L1范数组合的正则化约束。并且，TV和L1范数是不平滑的，它使得求解问题（1）变得更加困难。

目前，一些算法［11，12］已经应用到求解TVWL1模型。组合分裂算法（Composite splitting algorithm，CSA）和快速组合分裂算法（Fast composite splitting algorithm，FCSA）［9］是两个高效的求解问题（1）的算法，它们应用于图像重建，并且有着出色的图像重建质量。这两个算法都是把原问题分解为两个子问题：L1正则化约束和TV正则化约束。然后，采用已有技术分别求解这两个子问题来完成图像重建。然而，在少量观测值的情况下，这两个算法不能够完成对MR图像的高质量重建。

本文提出了一种新的求解TVWL1模型的图像重建算法，该算法将问题（1）分解为几个子问题，并利用分析稀疏表示特性构建子问题的求解算法，再通过组合子问题的解来获得重建图像。实验结果表明，与现有的一些算法相比，本文提出的新算法可以显著提高重建图像的主客观质量。

1 TVWL1模型求解算法

1.1 近似梯度算法

给出一个连续的凸函数g（x）和任意标量ρ＞0，凸函数g的prox算子定义如下［10，13］

无约束优化可以分为两部分

式（3）中f（x）是平滑的凸函数，g（x）是包含prox算子的凸函数。近似梯度算法定义如下

▽表示梯度算子，k是迭代步长。

1.2 组合分裂算法

可用于求解问题（1）的CSA算法［9］是一种基于变量和算子分裂技术的方法。它把组合正则化问题（1）分解为两个子问题：（1）把变量x分解为两个变量｛xi｝i＝1，2；（2）对｛xi｝i＝1，2独立执行算子分裂求解；（3）获得由｛xi｝i＝1，2线性组合而成的重建结果x。FCSA是CSA加入加速过程的改进算法。算法1给出了用于MR图像重建的FCSA/CSA。

算法1 快速组合分裂算法

（用rk＋1＝xk替换步骤（6）和（7）为CSA）

为了更加高效地求解问题（1），关键是如何应用相对稀疏度限制去改善图像的重建质量。

2 基于分析稀疏表示方法的图像重建算法

在学习求解问题（1）之前，首先回顾一下分析稀疏表示模型中的贪心分析追踪算法［14］。

2.1 贪心分析追踪算法

GAP算法是求解下面分析模型的L1最小化问题

式中x采用分析稀疏表示方式。GAP的核心思想是找到Φx的零元素的位置，并且从Φ中移除与Φx非零元素对应的行，用Φ的子集去执行迭代计算。算法2给出了GAP算法的详细描述。

算法2 贪心分析追踪算法

表达式（6）解析式为

GAP算法需要一个大的内存空间去存储大尺度矩阵并且在计算大尺度矩阵的逆矩阵时，有着很高运算复杂度，这样就很难实现大尺度图像的重建。

2.2 基于分析稀疏表示算法

本文提出的算法把图像重建问题分解为两个子问题：由近似梯度算法求解的TV正则化约束问题和由贪心分析追踪算法求解的L1正则化约束问题。这个算法被命名为基于分析稀疏表示算法（Analysis-sparse-based algorithm，ASBA）。图像重建算法框架如算法3所示。

算法3 基于分析稀疏表示算法

ASBA是一种迭代算法，它分为近似梯度求解和分析追踪求解两个过程。TV正则化约束问题很容易通过类似于CSA的近似梯度算法求解。而本算法的核心部分是由分析追踪算法求解L1正则化约束问题，即求解问题（6）的过程。

实际上，问题（6）可以转换为下面的形式

式中：A＝RTR＋，A是正定对称矩阵。所以，式（7）可以通过共轭梯度算法（Conjugate gradient algorithm，CG）高效求解。于是，结合CG算法和GAP算法得到一个新的算法——基于共轭梯度的GAP算法（Conjugate-gradient-based GAP，CGGAP），它在算法4中详细给出。

算法4 基于共轭梯度GAP算法

算法4中，x0是由近似梯度过程计算得到初始值是重建图像；n是共轭梯度算法迭代次数；K是GAP算法迭代次数；共轭梯度算法（CG）在算法5中给出；其他参数在第一部分有详细描述。

算法5 共轭梯度算法

在算法4和算法5中，由于部分傅里叶变换矩阵R和小波变换矩阵Φ尺度很大并且非常稠密，与其相关的计算需要很大的内存才能顺利完成。所以，在程序实现过程中，采用把大矩阵分块存储技术来降低运算对内存的需求。

实验部分将展示重建图像的主客观质量。实验结果表明，本文提出的新算法的图像重建质量要明显好于以前提出的一些算法。

3 实验结果和分析

实验的软件环境是 MATLAB7.6，实验代码是基于文献［12］提供的程序编写，实验内容是观测和重建自然图像和MR图像。实验图像包括自然图像Lena和Boat，医学图像 Heart，Brain，Chest和Artery，它们的尺度均为128＊128。观测方法采用部分傅里叶观测。采样率定义为m/n。噪声模型是由 Matlab产生的高斯噪声，定义为σ＊randn（m，1），σ＝0.01。正则化参数α和β分别设置为0.001和0.035。为了实验的公平性，在对比不同算法重建结果时，每个算法的迭代次数设置为10。小波变换矩阵Φ设置为2D小波。实验将比较每个算法的主客观图像重建质量，客观评价标准使用峰值信噪比（Peak-signal-to-noise-ratio，PSNR）。实验算法是ASBA，CSA和FCSA。

3.1 ASBA图像重建效果

图1是在不同采样率下ASBA重建自然图像Lena效果图，图1（a～d）对应的重建图像PSNR分别为21.68，28.72，32.34和38.62dB。图2是在不同采样率下ASBA重建MR图像Brain效果图，图2（a～d）对应的重建图像PSNR分别为18.13，24.42，27.85和35.21dB。从图中可以看出ASBA可以很好地求解TVWL1模型，并且随着采样率的提高，ASBA重建图像的主观质量明显提高。

图1 不同采样率下ASBA重建Lena图像效果Fig.1 The results of the reconstructed Lena images by ASBA at different sample rates

图2 不同采样率下ASBA重建Brain图像效果Fig.2 The results of the reconstructed Brain images by ASBA at different sample rates

3.2 图像重建主观评价结果

通过求解问题（1）来重建图像，从视觉效果上对本文提出的新算法ASBA和已有算法CSA＼FCSA进行评估。图3，4给出了采样率在25%时，ASBA，FCSA和CSA的图像重建质量对比。从图中可以看出，无论是自然图像还是医学图像，本文提出的ASBA重建的图像主观质量要明显好于FCSA和CSA。

图3 采样率25%不同算法重建Lena图像Fig.3 Lena image is reconstructed by different algorithm at sampling ratio near 25%

图4 采样率25%不同算法重建Brain图像Fig.4 Brain image is reconstructed by different algorithm at sampling ratio near 25%

图3（a）是Lena图像原图，图3（b～d）分别是由ASBA，FCSA，CSA重建得到的Lena图像，它们对应的PSNR分别是32.35，28.06和28.01dB。图4（a）是Brain图像原图，图4（b～d）分别是由ASBA，FCSA，CSA重建得到的Brain图像，它们对应的PSNR分别是27.92，24.19和23.53dB。

3.3 图像重建客观评价结果

图5是不同采样率下不同算法重建图像的PSNR对比曲线图。从图中可以看出，在不同采样率下，ASBA重建图像的PSNR与FCSA和CSA相比，可以提升2～5dB。表1列出的实验结果显示ASBA在所有采样率下的图像重建质量最好。

图5中曲线图上对应的四个采样率点是11%，16%，25%和43%。图5（a，b）图中的3条曲线分别对应ASBA，FCSA和CSA三种算法重建图像Lena和Brain的PSNR-Sampling Ratio曲线。

图5 重建图像PSNR-Sampling Ratio曲线Fig.5 PSNR-Sampling Ratio curve of the reconstruction images

表1 不同采样率下不同算法重建图像实验结果Table 1 Experimental results for image reconstruction over different sampling ratios with different methods

4 结束语

本文提出了一个基于分析稀疏表示的图像重建算法。该算法将TVWL1模型分解为TV和L1正则化约束两个子问题并分别进行求解，从而达到高效图像重建的目的。实验部分，通过与CSA和FCSA重建图像的比较，可以看出新提出的算法可以高质量重建MR图像和自然图像。

［1］Candès E J，Tao T.Near-optimal signal recovery from random projections：Universal encoding strate-gies［J］.Information Theory，IEEE Transactions on，2006，52（12）：5406-5425.

［2］Candès E J，Romberg J，Tao T.Robust uncertainty principles：Exact signal reconstruction from highly incomplete frequency information［J］.Information Theory，IEEE Transactions on，2006，52（2）：489-509.

［3］Donoho D L.Compressed sensing［J］.Information Theory，IEEE Transactions on，2006，52（4）：1289-1306.

［4］Candès E，Romberg J.Sparsity and incoherence in compressive sampling［J］.Inverse Problems，2007，23（3）：969.

［5］Romberg J.Imaging via compressive sampling［J］.Signal Processing Magazine，IEEE，2008，25（2）：14-20.

［6］Lustig M，Donoho D L，Santos J M，et al.Compressed sensing MRI［J］.Signal Processing Magazine，IEEE，2008，25（2）：72-82.

［7］Yang J，Zhang Y，Yin W.A fast alternating direction method for TVL1-L2signal reconstruction from partial Fourier data［J］.Selected Topics in Signal Processing，IEEE，2010，4（2）：288-297.

［8］Elad M，Milanfar P，Rubinstein R.Analysis versus synthesis in signal priors［J］.Inverse Problems，2007，23（3）：947-968.

［9］Huang J，Zhang S，Metaxas D.Efficient MR image reconstruction for compressed MR imaging［J］.Med-ical Image Analysis，2011，15（5）：670-679.

［10］Beck A，Teboulle M.A fast iterative shrinkagethresholding algorithm for linear inverse problems［J］.SIAM Journal on Imaging Sciences，2009，2（1）：183-202.

［11］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.

［12］Ma S，Yin W，Zhang Y，et al.An efficient algorithm for compressed MR imaging using total variation and wavelets［C］∥Computer Vision and Pattern Recognition，CVPR 2008.［S.L.］：IEEE，2008：1-8.

［13］Beck A，Teboulle M.Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems［J］.Image Processing，IEEE Transactions on，2009，18（11）：2419-2434.

［14］Nam S，Davies M E，Elad M，et al.The cosparse analysis model and algorithms［J］.Applied and Computational Harmonic Analysis，2013，34（1）：30-56.