图像分块压缩感知中的自适应粒重建算法*

李 然 孙艳歌,2 张清洁 刘宏兵

(1.信阳师范学院计算机与信息技术学院，信阳，464000； 2.北京交通大学计算机与信息技术学院，北京，100044)

引 言

奈奎斯特(Nyquist)频域采样定理是传统图像编码(例如JPEG和JPEG2000)的理论基础，其要求图像变换次数至少为图像总像素个数，即对图像实施全变换，这导致图像编码计算量较大，但大部分变换系数在编码过程中被丢弃，造成能耗浪费。由此可看出，具有高计算复杂度和低能耗利用率的传统图像编码不适合应用于“轻”采集点场合如：能耗受限的无线传感器网络节点。另外，全变换中稀少系数占据了大部分有用信息，若这些重要系数丢失，将极大衰退图像重建质量，因此，这也对无线通信中的容错机制提出挑战。压缩感知(Compressed sensing, CS)[1-2]以欠Nyquist速率变换信号，但仍可精确复原信号，这激发了图像压缩感知作为新图像编码方法出现[3]，将CS随机采样看作是图像变换，则其实现了在变换图像的同时以降维形式压缩图像，大大节约了编码成本，吸引了大量研究者的关注。

提高图像压缩感知的率失真性能是众多研究者致力研究的目标，目前主流方法是发掘可自适应图像非平稳统计特性的稀疏表示模型，以改善最小l1范数重建的收敛性能：Chen等[4]利用空间相关性建立多假设预测模型，复原更稀疏的残差改善重建质量；常侃等[5]改进了多假设模型，以分层重建方式恢复图像；Becker等[6]提出利用一阶分析法的 NESTA 算法，求解最小全变差(Total variation, TV)模型，确保稀疏分解的快速性和稳健性；Zhang等[7]利用非局部自相似性提出组稀疏模型，以图像子块构成群组，精简图像的稀疏表达式；王超等[8]提出贪婪迭代方法重建图像；Wu等[9]则引入局部自回归(Auto regressive, AR)模型追踪图像非平稳统计特性，提高稀疏表示基底解析图像的性能；郑海波等[10]提出自适应平均分组方法，通过改善测量方式获得重建图像质量提升。上述方法的率失真性能提升是以在重建端引入高计算复杂度为代价，往往随着图像维数增加，重建时间将急速上升，例如，Wu等[9]的算法重建一幅512×512的图像，需要约1个小时，这使之失去应用价值。相比于上述方法忽略计算复杂度而追求高重建质量，Gan[11]，Mun等[12]提出的平滑投影Landweber算法以低计算复杂度确保了良好的率失真性能，是性价比较好的算法之一，而该算法存在问题在于采用固定稀疏表示基底(例如，DCT基、小波基等)，无法在迭代过程中根据图像内容变化自适应更改基底构成，因此，其提高重建质量的潜力未被完全开发。主成分分析(Principal component analysis, PCA)[13]是可完全去冗余的最佳正交变换矩阵，文献[14]采用PCA在Landweber迭代过程不断学习更新稀疏表示基底，由于其可自适应于图像统计特性变化，因此，获得了一定的性能提升，然而，在PCA学习过程中，忽略了图像局部结构特性平稳，以全局统计特性更新基底，阻碍了PCA去相关的效果，影响基于Landweber迭代的重建算法性能提升。

针对Landweber迭代中，PCA未根据图像平稳局部结构特性更新稀疏表示基底的问题，本文提出自适应粒重建算法，借助粒计算(Granular computing, GrC)理论[15]，将图像分解为若干粒，任一粒是结构特性相似的图像子块集合，利用PCA学习出它们对应的最优表示基底，以获得高效的硬阈值收缩性能，提升图像重建质量。实验结果表明，本文提出算法可有效改善图像压缩感知的率失真性能，并确保了良好的主观视觉质量。

本文将探讨图像压缩感知中的自适应粒重建算法。首先，简述图像分块压缩感知与粒计算聚类两个背景知识；接着重点讨论提出的自适应粒重建算法，并分析仿真实验结果；最后对本文加以总结。

1 图像分块压缩感知和粒计算聚类

1.1 图像分块压缩感知

分块压缩感知(Block compressed sensing, BCS)[11]的框架如下所述。首先，尺寸为N = Ic×Ir的图像x被划分n个尺寸为B×B的块xi, i= 1,…,n, n=N/B2；接着，构造尺寸为B2×B2的高斯随机矩阵，由规范正交化生成正交矩阵Θ，随机提取MB(≪ B2)行，形成尺寸为MB×B2的测量矩阵ΦB，最后，对xi测量得到观测向量yi，表达式为

yi=ΦBxi

(1)

即原始图像的M个CS观测值，M = n×MB。整幅图像的总测量矩阵Φ是块对角矩阵，n个对角元素皆为ΦB，表达式为

(2)

BCS可采用最小线性均方误差准则(Minimummeansquareerror,MMSE)估计较好的初始图像，加速重构过程。图像块的MMSE初始解计算如下

(3)

(4)

式中ρ取0.9至1之间，凭经验ρ取0.95；分块尺寸B取32。

重建图像时，可建立各块观测向量yi与图像x的关系模型如下

(5)

式中E对应于图像分块排序的初等变换矩阵，Θ=Φ·E。利用式(5)可得l2-l1范数最小化重建模型如下

(6)

式中‖·‖2和‖·‖1分别为l2和l1范数，Ψ是稀疏表示基，λ为固定取值的正则化因子。式(6)可采用平滑投影Landweber算法以低计算复杂度求解，尤其是Ψ可在迭代过程中更新。因此，基于PCA的自适应稀疏表示基学习有望实现。

1.2 粒计算聚类

美国控制论专家Zadeh1979年提出GrC理论[15]，指出在现实生活中信息粒广泛存在，是对现实的一种抽象。GrC的3个重要概念是粒化形式、组织形式和因果关系[16]：粒化就是将研究对象的全体划分成各个部分；组织就是将部分通过一些特定算子融合成整体；因果关系则涉及原因与结果之间的联系。由此可看出，GrC与聚类具有天然联系，训练集可看作是研究对象全体，将样本视作单个粒，则粒化就是聚类，且组织实现了聚类的反过程，而因果关系描述了样本间存在的内在联系。不同于传统聚类方法(例如，K-means[17])，GrC聚类可灵活选择各种形状表示粒，粒间关系也可使用更加成熟的代数系统描述，因此具有更好的推广能力。

训练集的点群具有不规则形状，而GrC将粒表示为各种形状，例如，空间中的超菱形、超球、超正方形和超盒等，可较好地区分点群边界。粒的向量形式为G= (c, gr)，c是粒中心，gr是粒度。假设ci为包含在粒G中的点，则粒度为

(7)

式中‖·‖p为p范数，作为距离测度。不同距离测度表示不同形状的粒，例如，p=1表示对应的粒为超菱形粒，p=2表示对应的粒为超球粒，p=∞表示为超正方形粒。粒之间的算子包括合并算子∨与分解算子∧，合并算子∨就是将较小的两粒合并为较大的粒，分解算子∧就是将较大的粒分解为较小的粒，其中合并算子∨用于样本点的聚类。对于两粒G1= (c1, r1)和G2= (c2, r2)，其合并粒为

(8)

式中

图1 二维空间中球形粒合并Fig.1 Merge of spherical granules in two-dimensional space

图1显示了二维空间两个球形粒G1= (0.1, 0.2, 0.2)和G2= (0.25, 0.25, 0.2)的合并过程，它们的合并结果为粒G1∨G2= (0.175 0, 0.225 0, 0.279 1)。

根据粒向量化表示方法，粒间关系必然涉及向量间关系，而向量间偏序关系与粒间包含关系是不一致的，因此认为粒间具有模糊包含关系。为了度量粒合并时的模糊包含关系，引入评价函数v(G)，设计如下算子

(9)

式中v(G):G→R是粒空间到实数空间的映射，采用下述非线性正评价函数

(10)

针对聚类问题的训练集S，构造由粒集G、模糊包含关系μ(·,·)、合并算子∨组成的代数系统〈G,μ〉，通过粒间模糊包含程度，控制粒合并，逐步划分S至若干个粒，而每个粒就对应于一个类别。图2显示了粒表示分别为超球形与正方形时，GrC以模糊包含关系与粒度大小控制有条件的粒合并，实现了良好的聚类效果，尤其是对于包含度的调控十分有利于正确区分边界处训练样本的归属。对于图像信号来说， 在某空间邻域内， 各块在结构特征上有着天然联系，利

图2 二维空间球粒与正方形粒的GrC聚类结果Fig.2 GrC clustering results of spherical and square granules in two-dimensional space

用GrC聚类就可自适应于图像内容实现块的自动划分。各聚类中，块结构特征极其相似，更有利于寻找到专门适应于固定数据模式的基本波形实现稀疏表示。

2 自适应粒重建算法

2.1 平滑投影Landweber迭代

在BCS框架下，采用凸优化算法(例如，NESTA[6])仅能逐块独立重建，而由于图像局部结构特征非平稳变化，造成稀疏性分布不均匀，因此，各块复原水平不一致，出现块效应。另外，分块尺寸较小也制约了分块重建性能，造成重建图像出现大量噪声。为了克服上述问题，文献[11] 结合凸集投影、硬阈值收缩和维纳滤波，提出基于平滑投影Landweber迭代的CS图像重建算法，其可在低计算复杂度下确保良好的重建质量，具体的算法流程如表1所示。在MMSE线性估计初始解中，具有大量块效应和噪声，因此，在每次迭代中，首先，实施3×3的维纳滤波运算，抑制块效应以及平滑图像；然后，再投影滤波后的图像回凸集С = {xi:yi=ΦBxi, i=1,…,n}；接着，作基于稀疏表示基Ψ的图像变换，并对变换系数作硬阈值收缩，消除残留噪声；最后，投影回凸集С。重复上述流程，直到满足终止准则为止。硬阈值收缩过程中，由于噪声存在，迭代解不稀疏，通过硬阈值判决强制噪声分量为0，得到更纯净的稀疏表示系数。

表1 平滑投影Landweber迭代流程

2.2 基于粒的PCA硬阈值收缩

阈值收缩过程是否能够有效地消除噪声将决定平滑投影Landweber迭代的重建性能优劣，文献[11]采用重叠块离散余弦变换和非抽样小波变换作硬阈值去噪，但它们并不能较好地捕获点、线、边缘等方向特征，导致重建图像质量不理想。为了解决上述问题，文献[12]提出基于轮廓波、双树小波等方向小波的双收缩策略，获得了一定性能提升。上述方法均是构造各种特定小波使图像有用信息尽量多地集中于少数表示系数上，而删除大量近似为0的表示系数，从而消除无用噪声。文献[18]指出，采取自适应含噪图像内容的稀疏冗余字典作硬阈值收缩去噪，会比特定小波更易于表示图像重要结构，但K-SVD[19]等字典学习算法计算复杂度高，因此，不适合于要求低计算复杂度的平滑投影Landweber迭代。文献[13]提出若干PCA硬阈值收缩去噪方法，提取含噪图像中图像子块作为样本，利用PCA训练正交变换矩阵，去除样本间空间相关性，将图像的有用信号和噪声分离，从而保证了硬阈值收缩去噪的高效性。另外，PCA硬阈值收缩去噪的计算复杂度不高，十分适合于平滑投影Landweber迭代。由于图像仅在局部具有平稳统计特性，而在全局范围内结构特征变化较大，因此，将全局所有图像子块作为样本组成训练集，生成单一PCA矩阵，会降低PCA解析稀疏性的能力。更加符合实际的做法是，根据局部平稳特性，划分全体图像子块至若干子集，任一子集均是结构特征相近的子块集合。那么，利用PCA生成各子集的稀疏表示基底，可专属于某特定结构特征，保证子集内所有样本均有简洁的稀疏表达式。

根据上述思路，在利用PCA学习稀疏表示基底之前，先对图像子块进行聚类。K-means[17]是主流的聚类算法，但它需预先设定聚类个数，无法根据数据模式自适应分类，而图像非平稳统计特性导致图像子块的结构特征类别差异较大，那么，K-means十分不利于图像子块划分。GrC不需设定固定的聚类个数，而根据粒间模糊包含关系，通过粒度阈值调控粒合并，得到较为合理的聚类个数，因此本文选择自适应于数据模式的GrC聚类图像子块。基于粒的PCA硬阈值收缩具体步骤如下：

步骤2以pi作为最细的原子粒，构成训练集S，设置粒度控制参数ρ，利用合并算子∨对粒进行有条件合并，直至所有训练样本都包含在某个粒中，最终每个粒对应一个类别。S的聚类流程如表2所示，输出由L个粒组成的粒集GS= {G1,G2,…,GL}。

步骤3计算第l个粒Gl的d×d协方差矩阵Σl如下

(11)

(12)

式中|Gl|为粒Gl中包含的样本个数。

步骤4对协方差矩阵Σl作特征值分解，得特征值λ1≥…≥λd，及对应的规范化特征向量e1,e2, …,ed。

步骤5列排列规范化特征向量e1,e2, …,ed组成d×d正交变换矩阵Γl= [e1,e2,…,ed]，正交变换粒Gl内各图像子块作

(13)

步骤6作硬阈值收缩如下

(14)

(15)

(16)

步骤7反变换回像素域：

i

(17)

步骤8对粒集GS中所有粒作硬阈值收缩后，反变换回子块，并重新合并成整幅图像。由于子块互有重叠，对于重叠区域，取所有重叠像素的平均值。

表2 GrC聚类划分训练集S流程

由上述流程可知，对于每个粒均对应一个PCA矩阵，但是所有粒包含的样本总体仍保持不变，所以与学习全局PCA矩阵所需的计算复杂度相比，基于粒的PCA学习算法仍保持了低计算复杂度，仅是GrC聚类增加了额外的计算复杂度O(m·Wp2)。

3 实验结果与分析

采用5幅包含不同程度平滑、边缘和纹理细节的512×512标准测试图像Lenna，Barbara，Parrot，Leaves和House测试本文提出的算法性能，其参数设置如下：正则化因子λ为2.5，平滑投影Landweber迭代的终止容限ε为0.01，图像子块尺寸Wp为7，粒度阈值T为0.03。一方面，对比不同稀疏表示基底下硬阈值收缩性能，提出算法采用的是基于粒的PCA基(简称为GrC-PCA)，参考基底有DCT变换基、双树小波变换基(简称为DDWT)和全局PCA基(简称为GPCA)；另一方面，比较提出算法与多假设预测平滑Landweber算法(简称为MH_SPL)[4]、基于TV最小的NESTA算法(简称为NESTA_TV)[6]和组稀疏算法(简称为GSR)[7]之间的重建性能。在所有实验中，分块尺寸B为32，预设总测量率SR分别取0.2～0.6。评价客观性能的指标采用峰值信噪比(Peak signal-noise ratio, PSNR)，但考虑到随机变化的测量矩阵，实验中5次重建图像并计算PSNR值取其平均。运行实验的计算机配置为：主频3.40 GHz的酷睿i7 CPU，Windows 7 64位操作系统和Matlab 7.6仿真实验软件。

3.1 硬阈值收缩性能对比

采用DCT基、DDWT基和GPCA基作硬阈值收缩，与本文算法使用的GrC-PCA基对比5幅测试图像的平均率失真曲线，如图3所示。可看到在任何测量率下，GrC-PCA基均获得最高的PSNR值，尤其是优于DCT基1 dB左右。由于考虑了局部平稳结构特性，相比于GPCA基，GrC-PCA基也获得了0.2 dB左右的PSNR增益。图4显示了采用不同聚类方法下本文硬阈值收缩方案的性能对比，可看出若采用传统K-means实施图像子块聚类， 率失真性能有所衰退， 而采用GrC聚类方法， 在任何测量率下均获得了PSNR值提升，尤其是在测量率为0.2时，PSNR增益可达0.6 dB左右。

图3 不同稀疏表示基底下硬阈值收缩性能对比

Fig.3 Comparison of hard-thresholding shrinkage performances under the different sparse representation bases

图4 采用不同聚类方法硬阈值收缩性能对比

Fig.4 Comparison of hard-thresholding shrinkage performances using the different clustering methods

3.2 重建性能对比

为了评估本文整体算法的率失真性能，选择如下目前较流行的重建算法作参照，即MH_SPL[4]、NESTA_TV[6]和GSR[7]。表3列出了测量率由0.2～0.4取值时，5幅测试图像的平均PSNR值，可看到本文算法明显优于MH_SPL与NESTA_TV算法，分别获得了0.48 dB和3.27 dB的PSNR增益。然而，当对比于GSR算法时，本文算法平均衰退了0.64 dB，这是由于GSR采用的组稀疏模型不仅考虑了局部结构特征平稳特性，而且也考虑了图像子块的非局部相似性。因此，它比GrC_PCA基具有更好的解析稀疏能力。但是，如表4所示，GSR算法具有较高的计算复杂度，当在不同测量率下重建Lenna图像时，平均重建时间高达2 542.65 s，而本文算法仅需30.20 s，由此可知，若同时考虑重建计算复杂度，本文算法具有较好的性价比。与采用了DDWT基的MH_SPL算法相比，本文算法平均增加了18.54 s，这是由于PCA需在每次迭代中重新学习稀疏表示基底，因此增加了一定的计算复杂度，但也带来了PSNR值增益。图5显示了测量率为0.3时，不同算法重建Lenna图像的主观视觉质量对比，可看到NESTA_TV重建图像包含了大量块效应，MH_SPL算法消除了块效应，但却引入了噪声，而本文算法获得了与GSR算法相近的主观视觉质量。

表3不同算法重建5幅测试图像的平均PSNR值比较

Tab.3ComparisonofaveragePSNRvaluesoffivetestimagesusingthedifferentmethodsdB

重建算法测量率SR0．20．30．4平均MH＿SPL31．6233．5234．8233．32NESTA＿TV28．5630．8932．1530．53GSR32．5634．8235．9534．44提出算法32．1234．0835．2133．80

表4不同算法重建Lenna图像所消耗时间对比

Tab.4ComparisonofconsumedtimewhenusingthedifferentmethodstoreconstructLennaimages

重建算法测量率SR0．20．30．4平均MH＿SPL6．5912．5215．8811．66NESTA＿TV44．9844．6344．8144．81GSR2272．152519．232836．562542．65提出算法32．0530．6127．9530．20

图5 测量率为0.3时，不同算法重建Lenna图像的主观视觉质量对比Fig.5 Comparison of subjective visual quality using the different methods to reconstruct Lenna image when the measurement rate is 0.3

4 结束语

本文提出采用GrC聚类改善平滑投影Landweber迭代过程中的PCA硬阈值收缩性能。考虑到图像的局部平稳统计特性，首先，将图像分解成若干粒，由于GrC可依据数据模式自适应划分训练集，因此，粒内样本具有相似的结构特征；接着，利用PCA学习专属于各粒的稀疏表示基底，并对粒内子块作硬阈值收缩去噪；最后，合并各子块形成整幅图像。实验结果表明，本文算计可确保良好的率失真性能，重建图像主观视觉质量较好，且具有较低的计算复杂度。

尽管GrC聚类与平滑投影Landweber迭代的结合能有效改善图像压缩感知率失真性能，但是仍与目前高级别算法具有一定差距。GrC聚类与压缩感知测量策略的结合也具有提升率失真性能的潜力，因此，在下一步工作中将研究基于GrC聚类的压缩感知自适应测量方法，进一步提高图像压缩感知的率失真性能。

[1] Candes E J, Romberg T, Tao T. Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[J]. IEEE Transactions on Information Theory, 2006,52(2):489-509.

[2] Donoho D L. Compressed sensing[J]. IEEE Transactions on Information Theory, 2006,52(4):1289-1306.

[3] Deng C, Lin W, Lee B S, et al. Robust image coding based upon compressive sensing[J]. IEEE Transactions on Multimedia, 2012,14(2):278-290.

[4] Chen C, Tramel E W, Fowler J E. Compressed sensing recovery of images and video using multihypothesis predictions[C]//Proc of Conference Record of the 46th Asilomar Conference. Pracific Grove, CA: IEEE Signal Processing Society Press, 2011:1193-1198.

[5] 常侃,覃团发,唐振华.基于多重假设的视频压缩感知分层重建[J].数据采集与处理,2013,28(6):730-738.

Chang Kan, Qin Tuanfa, Tang Zhenhua. Multi-hypothesis-based hierarchical reconstruction for compressed video sensing[J]. Journal of Data Acquisition and Processing, 2013,28(6):730-738.

[6] Becker S, Bobin J, Candès E J. NESTA: A fast and accurate first-order method for sparse recovery[J]. SIAM Journal on Imaging Sciences, 2011,4(1):1-39.

[7] Zhang Jian, Zhao Debin, Xiong Ruiqin, et al. Group-based sparse representation for image restoration[J]. IEEE Transactions on Image Processing, 2014,23(8):3336-3351.

[8] 王超.基于压缩感知的贪婪迭代重构算法[J].数据采集与处理,2012,27(S2):298-303.

Wang Chao. Greedy iterative reconstruction algorithm based on compressive sensing[J]. Journal of Data Acquisition and Processing, 2012,27(S2):298-303.

[9] Wu Xiaolin, Dong Weisheng, Zhang Xiangjun, et al. Model-assisted adaptive recovery of compressed sensing with image applications[J]. IEEE Transactions on Image Processing, 2012,21(2):451-458.

[10] 郑海波,朱秀昌.基于分块压缩感知与平均分组的图像多描述编码[J].数据采集与处理,2014,29(5):764-769.

Zheng Haibo, Zhu Xiuchang. Image multi-description coding method based on equally grouping and block compressive sensing strategy[J]. Journal of Data Acquisition and Processing, 2014,29(5):764-769.

[11] Gan L. Block compressed sensing of natural images[C]//Proc of International Conference on Digital Signal Processing. Cardiff, UK: IEEE Signal Processing Society Press, 2007:403-406.

[12] Mun S, Fowler J E. Block compressed sensing of images using directional transforms[C]//Proc of International Conference on Image Processing. Cario, Egypt: IEEE Signal Processing Society Press, 2009:3021-3024.

[13] Deledalle C A, Salmon J, Dalalyan A. Image denoising with patch based PCA: Local versus global[C]//Proc of the 22nd British Machine Vision Conference. Dundee, British: Springer, 2011:1-10.

[14] 李然,干宗良,朱秀昌.基于PCA硬阈值收缩的平滑投影Landweber图像压缩感知重构[J].中国图象图形学报,2013,18(5):504-514.

Li Ran, Gan Zongliang, Zhu Xiuchang. Smoothed projected Landweber image compressed sensing reconstruction using hard thresholding based on principal components analysis[J]. Journal of Image and Graphics, 2013,18(5):504-514.

[15] Zadeh L A. Fuzzy sets and information granulation[M]. Washington: North Holland Publishing, 1979.

[16] Zadeh L A. Toward a theory of fuzzy information granulation and its centrality in human reasoning and fuzzy logic[J]. Fuzzy Sets and Systems, 1997,90(2):111-127.

[17] Bishop C M. Pattern recognition and machine learning[M]. New York: Springer, 2006:424-428.

[18] Elad M, Aharon M. Image denoising via sparse and redundant representation over learned dictionaries[J]. IEEE Transaction on Image Processing, 2006,15(12):3736-3745.

[19] Aharon M, Elad M, Bruckstein A. K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation[J]. IEEE Transaction on Signal Processing, 2006,54(11):4311-4322.