一种基于最小二乘优化的压缩感知图像快速重建方法

摘要

本发明公开了一种基于最小二乘优化的压缩感知图像快速重建方法，基于最小二乘方法实现信号的优化重构，利用新定义的整体相关性度量参数-整体相关度，选择针对图像信号的最相关原子，减少迭代次数,引入分块重构理论并重新设计分块大小和测量矩阵,降低了重构操作的计算规模；基于最小二乘方法实现信号的优化重构,保证了重构精度和收敛速度。实验结果表明,FBWRFI算法也可以显著降低信号的重构时间，并使随信号增大而高速增长的重构时间的增长趋势变为线性,证明了算法的有效性。

说明书

技术领域

本发明涉及一种压缩感知图像重建方法，具体地涉及一种面向图像的、基 于最小二乘优化的压缩感知采样快速重建方法。

背景技术

基于稀疏表示理论和泛函分析-逼近理论的压缩感知(compressedsensingor compressivesampling,CS)是一种新型的采样方法，在对可压缩信号采样时能够突破著 名的Shannon-Nyquist采样定理的限制，以远低于两倍信号最大带宽的采样率对信号采 样，从而减少了采集的数据量。以其用少量采样值表示全长的信号的优点，压缩感知在 出现之初就吸引了信号、通信、电子信息、统计理论、编解码理论和计算机等众多领域 的研究热情，被认为是信息科学近年来最重大的研究成果。压缩感知基本原理如下：

假设一维离散信号s∈R^N×1是可压缩的(这是前提条件)，压缩感知方法的信号采样 是利用特别设计的观测矩阵Φ∈R^M×N(这里M＜＜N)将可压缩信号s从N维投影到M 维，即：y＝Φs；这里的y∈R^M×1是信号s的压缩采样，而采样y的长度M与信号s的长 度N与之比M:N就是在压缩感知框架下的采样率，该采样率一定小于1。由于y＝Φs是 欠定方程组，无法直接求解，压缩感知方法从近似重构的角度来考虑求解问题。

考虑信号s可压缩的前提和稀疏表示理论，对于可压缩信号s一定存在某稀疏变换基 Ψ∈R^N×N，使得s在该变换基上的表示θ∈R^N×1是稀疏的，即：

(1)s＝Ψθ，

式中θ的大部分元素为0或接近于0；通常θ中非0元素的个数K(K＜＜M＜＜N)被 称为θ的稀疏度。从而有y＝Φs＝ΦΨθ＝Aθ；这里的A＝ΦΨ被称为信息算子或感知 矩阵。y＝Aθ可以通过0-范数优化问题求解

(2)min||θ||₀s.t.y＝Aθ。

而根据Donoho、Candès等人研究提出的压缩感知理论，在感知矩阵A具有有限等距性 质(RestrictedIsometryProperty,RIP)的条件下，0-范数优化可以转化为1-范数优化问 题：

(3)min||θ||₁s.t.y＝Aθ，

其中，A具有RIP性质的等价条件是观测矩阵Φ与稀疏变换基Ψ不相干，这一点在设计 观测矩阵Φ时已得到保证。

由于θ是稀疏的，通过(2)和(3)式可以得到其精确或近似逼近解然后利用(1) 式得到原始信号s的精确或近似逼近解常用压缩感知算法主要包括针对(2)式的贪 婪算法和针对(3)式的凸松弛算法。最经典的压缩感知算法就是凸松弛算法中的基追踪 (BasisPursuit,BP)算法，所需的观测值数目少(即采样率低)、计算精度高；匹配追 踪(MatchingPursuit,MP)、正交匹配追踪(OrthogonalMP,OMP)及其改进算法都是 贪婪算法的代表，收敛速度快、计算复杂度较低，实际应用较多。

压缩感知方法的信号重建是近似优化重构，非常适合于处理图像信号，目前其已在 磁共振、X-射线扫描、雷达成像、遥感图像处理、超频谱图像分析等方面有了很多应用， 并已开发出单像素相机。压缩感知算法的提出都是以一维信号为处理对象的，对于二维 信号一般以列或行向量为重建单位分别重建，当每一个列或行信号都被重建后，整个二 维信号也就重建成功了，图像信号的重建通常也是按列处理的。我们把图像信号S记为 [s₁,s₂,…,s_n]，其中s_i(i＝1,2,…,n)为信号的第i列，在处理中每一个s_i分别进行重构。

压缩感知方法具有采样数目少的优点，但其应用仍然有很多难题等待解决，较长的 信号重构时间是其中的一个关键。如BP算法的计算复杂度高达O(N³)，收敛速度较快 的OMP算法的计算复杂度也达到了O(NK²)。在配置为Intel(R)Core(TM)2DuoCPU E46002.4GHz、3GB内存、Windows7、MatlabR2011b的个人计算机上的重构不同 大小信号时的重构时间如表1和2所示。

表1OMP和BP算法重构一维信号时，信号大小每增加1倍重构时间的变化

表2OMP和BP算法重构二维信号时，信号大小每增加1倍重构时间的变化

压缩感知算法重构信号所需的时间较长，并随信号增大以远大于信号增幅的速率增 长，而且这个速率随信号的增大还在不断增长。对比表1和表2可以发现，二维信号的 重构时间远远大于一维信号，且其增长幅度也较大；如大小为1024×1的信号重构时间是 512×1的5.55倍，而大小为1024×1024的信号重构时间是512×512的11.76倍。

我们在压缩感知算法及其应用的研究中，利用日渐普及的高性能计算设备和云计算 平台加速压缩感知算法，以提高算法的执行速度，验证了压缩感知方法在物联网数据采 集中的应用，实现了压缩感知算法的并行化、多/多核CPU加速、GPU加速和云平台加 速等，加速后的算法执行时间如图1和表3所示。

表3OMP算法重构信号时，不同数量计算资源的云加速效果(编程语言Python，云平台)

从图1和表3中可知，压缩感知算法的并行化及加速虽然能显著减少算法的执行时 间，但不能改变算法随信号增大而以极高的速率增长的趋势。压缩感知算法的重构时间 高速增长，且增长趋势不断扩大，使重构算法的执行时间无法估计，不利于在云平台加 速时的自动资源分配。

发明内容

针对上述技术问题，本发明目的是：提供一种基于最小二乘优化的压缩感 知图像快速重建方法FBWRFI，它采用最小二乘优化方法实现图像的重建，以保证收敛速 度和重构精度；以每一个原子和二维残差之间的整体相关度测量取代原子与一维残差之 间的相关度测量，以减少迭代次数、降低计算复杂度；引入块压缩感知理论并重新设计 快大小和测量矩阵，以降低优化迭代中的参数规模。该方法能够显著降低从压缩的采样 中重建原始图像的时间，并使图像的重构时间由随信号增大而以极高速率增加变为线性 增长，提高了重构时间的可预估性，有利于使用云资源加速时的资源智能分配。

本发明的技术方案是：

一种基于最小二乘优化的压缩感知图像快速重建方法，其特征在于，包括以下步 骤：

S01：把原始图像信号S:N×N基于压缩感知方法的压缩采样Y:M×N，根据块测量 矩阵Φ_B:M_B×B划分为个互不重叠的块Y_i:M_B×B，其中是基于压缩感 知方法进行压缩采样时的采样率，

S02：利用WRFI重构每个块S_t:B×B，压缩采样为Y_t:M_B×B，稀疏系数记Θ_t:B×B， 其过程如下：

S02-01：计算块感知矩阵A_B＝Φ_BΨ_B:M_B×B中每一个原子α_i(i＝1,2,…,B)与当前二 维残差R(k-1):M_B×B的整体相关度其中k表示当前是第k次 迭代，α_i是A_B的第i列，二维残差R(k-1)的初始值为Y_i并在每次迭代后更新，R(k-1)_j是R(k-1)的第j列，Ψ_B:B×B是稀疏变换基且Ψ_B和Φ_B线性无关；

S02-02：通过对比ρ_i(i＝1,2,…,B)，在序号不包含在原子序号集合Λ_k-1的原子中选出 与当前残差最大相关原子的序号λ_k＝max_{i＝1,2,…,N}(ρ_i)；其中Λ_k-1是前面k-1次迭代中选 择的最相关原子的序号组成的集合，λ_k表示第k次迭代选出的最相关原子的序号；

S02-03：把λ_k并入集合Λ_k-1，把原子并入集合A_k-1得到：Λ_k＝[Λ_k-1,λ_k]， 其中是A_B的第λ_k列，A_k-1是前面k-1中得到的最相关原子的集合；

S02-04：基于A_k和最小二乘优化方法获取本次迭代的近似逼近解Θ_tk；Θ_tk是第k次 迭代得到的稀疏系数Θ_t的近似逼近；

S02-05：判断迭代条件：达到设定的重构精度或设定的最大迭代次数K，则结束重 建块S_t的迭代计算并跳转到S02-07；Θ_tk是S_t的稀疏表示系数Θ_t的近似逼近解

S02-06：更新参数R(k)＝Y-A_kΘ_tk和k＝k+1，跳转至S02-01继续迭代；这里k是 新的迭代次数，R(k)是新的二维残差；

S02-07：利用S_t＝Ψ_BΘ_t得到S_t的近似逼近解

S02-08：判断是否每个块S_t都已经重建：即当t≥n则跳至S03，否则跳至S02重构 下一个块S_t+1；

S03：将所有的重构信号拼接得到完整的重构信号

优选的，在重构过程之前还包括采样过程，采样过程包括以下步骤：

S11：把原始图像信号S:N×N分割为个大小为B×B的块，若原始信号 边缘的分块大小不足B×B，将不足的部分的元素填充为0；

S12：对于每一个分块S_i:B×B(i＝1,2,…,n)，使用相同的稀疏变换基Ψ_B:B×B进行 稀疏变换S_i＝Ψ_BΘ_i，其中Θ_i:B×B是S_i的稀疏系数；

S13：对于每一个分块S_i:B×B(i＝1,2,…,n)，使用相同的测量矩阵Φ_B:M_B×B进行 投影变换Y_i＝Φ_BS_i，其中Y_i:M_B×B是S_i的压缩采样，就是采样率或压缩比，而 Y_i＝Φ_BS_i＝Φ_BΨ_BΘ_i可记为Y_i＝A_BΘ_i，A_B就是压缩感知方法的感知因子；

S14：把每一个Y_i(i＝1,2,…,n)拼接起来组成Y就是原始图像信号 S的压缩采样，就是采样率或压缩比。

优选的，所述块测量矩阵Φ_B的大小为M_B×B，分块大小大于等于256×256，在 分块时当原始信号边缘的分块大小不足256×256，将不足的部分的元素填充为0；当重 构完成后，将填充部分的值去掉。

与现有技术相比，本发明的优点是：

本发明采用最小二乘优化方法实现图像的重建，以原子和二维残差之间 的整体相关度测量取代原子和一维残差之间的相关度测量，减少了迭代次 数、降低计算复杂度，引入块压缩感知理论降低优化迭代中的参数规模和计 算规模，并使重构时间随信号增大呈线性增长，提高了重构时间的可预估性。

附图说明

下面结合附图及实施例对本发明作进一步描述：

图1为OMP算法重构信号时，多CPU并行的加速效果图；

图2为本发明基于最小二乘优化的压缩感知图像快速重建方法的整体 重构设计的流程图；

图3为本发明基于最小二乘优化的压缩感知图像快速重建方法的原理 图；

图4为使用本发明方法重构不同图像的效果图；

图5为在不同采样率下重构图像的效果图；

图6为重构不同大小图像的效果图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明了，下面结合具体实施 方式并参照附图，对本发明进一步详细说明。应该理解，这些描述只是示例 性的，而并非要限制本发明的范围。此外，在以下说明中，省略了对公知结 构和技术的描述，以避免不必要地混淆本发明的概念。

实施例：

本发明主要包括两部分：整体重构设计WRFI和基于分块压缩感知的 快速重建方法FBWRFI，它们之间的关系是：WRFI是以原子和二维残差之 间的整体相关性测量取代原子和一维残差之间的相关性测量后，面向图像整 体的、较少迭代次数的快速重建方法；FBWRFI是对需要重建的图像分块后， 每一个图像块均采用WRFI方法进行重建，在保证较少迭代次数的前提下降 低迭代中的参数规模和计算量，以获取更快的重建速度和可预估的重构时 间。

WRFI在重构过程中将二维信号作为一个整体来处理，而不是按列分别处理。为了 测量二维信号残差与测量矩阵各列(原子)之间的相关度，就需要定义新的相关度测量 参数，称之为二维信号的整体相关度(whole-correlation)。

定义1(整体相关度).向量α∈R^m×1和矩阵Φ∈R^m×n的整体相关度(whole-correlation) 为其中ρ为整体相关度，是是矩阵Φ的第i列。

对原始信号S:N×N(为了简单，假设S的行和列相同)，设压缩采样Y:M×N， 测量矩阵Φ:M×N，稀疏变换基Ψ:N×N，记A＝ΦΨ:M×N且有S＝AΘ，Θ:N×N是 稀疏的。WRFI的主要思想是：在每次迭代时(第t次迭代)，首先计算感知矩阵A中的 每一个原子(A的每一列)α_i:M×1(i＝1,2,…,N)与当前二维残差R(t-1):M×N(第t-1 次迭代时更新的残差，初始值R(0)＝Y)的整体相关度ρ_i，并选出最大相关度 这里最相关原子的序号是λ_t，最相关原子是然后将和A_t-1合并为这里的A_t-1是前t-1次迭代选出的最相关原子组成的矩阵；接着基 于A_t和最小二乘优化方法获取本次迭代的近似估计：

当近似估计达到设定的精度或迭代计数器k达到设定的最大迭代次数K时，迭 代结束，本次迭代的近似估计就是Θ的重建结果；否则更新残差和增量t，并继续迭 代。WRFI的算法流程如图2所示。

在本发明方法中，对信号S:N×N，其稀疏变换基Ψ:N×N和测量矩阵 Φ:M×N是线性无关，这是前提条件。在S＝ΨΘ和Y＝ΦS＝ΦΨΘ＝AΘ中， Y:M×N、Θ:N×N、A＝ΦΨ:M×N；对于S的每一列s_i:N×1都有s_i＝Ψθ_i(i＝1,2,…N)和y_i＝Φs_i＝ΦΨθ_i＝Aθ_i(i＝1,2,…,N)，因Ψ和Φ是线性无关的，根据压缩 感知理论所有的s_i(i＝1,2,…,N)都可以被高概率地重构，从而作为整体的S也就可以被高 概率精确重构，也即本方法理论上是可行的。

WRFI通过将选择原子用于二维信号整体重建降低了原子选择和重构操作中的迭代 次数，但并不能改变重构时间的随信号增大而带来的高增长率。这主要是因为随着信号 规模的增加，算法中使用的稀疏变换基、观测矩阵及其它一些中间参数的规模也随之增 大，从而带来了优化算法中每一次迭代的计算量呈指数级增加。为了降低重构算法所需 的存储和计算的规模，将信号的大小保持在一个相对较小的程度是必须的。在面向图像 的二维整体重构得到WRFI的基础上，引入了块压缩感知(blockcompressedsensing)， 并最终形成了基于最小二乘优化的快速图像重建压缩感知方法FBWRFI方法。

FBWRFI首先把原始图像S:N×N被分割成若干个大小为B×B的小块(图像)。基 于压缩感知理论，图像S的稀疏变换就是对每一个块图像S_i:B×B使用相同的稀疏变换 基Ψ_B:B×B进行操作：S_i＝Ψ_BΘ_i，Θ_i:B×B就是块图像S_i的稀疏系数；图像S采样就 是使用相同的测量矩阵Φ_B:M_B×B对每一个块图像S_i进行投影变换Y_i＝Φ_BS_i；若记 A_B＝Φ_BΨ_B:M_B×B，则Y_i＝A_BΘ_i，其中M_B:B可表示采样率/压缩比。因为设计测量矩 阵Φ_B:M_B×B时总是M_B＜＜B，采样率也就小于1，实现了压缩采样。对于整个信号S来 说，稀疏变换基Ψ:N×N和观测矩阵Φ:M×N(M＝M_B×k，)分别如(5) 式所示。

FBWRFI方法的(整体重构设计WRFI)原理与过程如图3所示。FBWRFI方法中 设计的块测量矩阵Φ_B的大小为M_B×B，综合考虑重构时间和重构精度后把分块大小定 义为256×256，这个分块大小是可以调整的，调整的原则基于实际应用中对重构图像精 度的需求，分块越大则重建图像的精度越高。在分块时，原始信号边缘各个分块大小可 能不足B×B，需要将不足部分的元素值填充为0；当重构完成后，再将这部分的值去掉 即可。

本发明的具体实现步骤，包括：

输入：

原始信号S:N×N，压缩采样Y:M×N，块测量矩阵Φ_B:M_B×256，块稀疏变换基 Ψ_B:B×B，块稀疏度K；

输出：

重构信号$\hat{S}:N\times N.$

算法过程：

(1)FBWRFI的初始化

1)把压缩采样Y根据Φ_B分块Y_i:M_B×256(i＝1,2,…,n)，为块数量；

2)t＝1，t为块计数器；

3)存放重构信号的参数

4)最大块数n；

(2)使用WRFI重构每一个块图像S_t:256×256(其压缩采样为Y_t，t＝1,2,…,n)

1)WRFI的初始化

a)残差；R₀＝Y_t

b)已选出原子序号集合

c)已选择的原子组成的矩阵

d)迭代次数计数器k＝1；

2)获取块图像S_t稀疏表示系数的Θ_t:256×256近似解的迭代过程

a)计算块感知矩阵A_B＝Φ_BΨ_B:M_B×256中每一原子α_i(i＝1,2,…,256)与当前二维残 差R(k-1):M_B×256的整体相关度R(k-1)_j是R(k-1)的第j列；

b)通过对比ρ_i(i＝1,2,…,256)，在序号不包含在集合Λ_k-1的原子中选出与当前残差 最大相关的原子序号λ_k：λ_k＝max_{i＝1,2,…,256}(ρ_i)；

c)把λ_k并入集合Λ_k-1、把原子并入集合A_k-1得到：Λ_k＝[Λ_k-1,λ_k]，

d)利用A_k和最小二乘优化方法得到本次迭代的近似逼近解Θ_tk；

e)判断迭代结束条件：Θ_tk达到设定重构精度或迭代计数器k达到最大迭代次数K， 则结束块S_t的重建迭代并跳转至步骤(2)-5)，这时最后一次迭代中得到的Θ_tk就是S_t的稀 疏表示系数Θ_t的近似逼近解否则更新参数R(k)＝Y-A_kΘ_tk和k＝k+1，并返回步骤 (2)-2)-a)继续迭代；

5)利用S_t＝Ψ_BΘ_t得到S_t的近似逼近解

(3)保存重建的块并使t指向下一个需要重构的块

1)把第t块重构信号存入

2)t＝t+1；

(4)判断循环结束条件

1)如果t＞n，退出循环；

2)否则，返回步骤(2)继续循环；

(5)得到完整的重构信号

计算复杂度和空间复杂度

对于图像S:N×N，稀疏变换基Ψ:N×N，测量矩阵Φ:M×N，稀疏表示系数 Θ:N×N，压缩采样Y:M×N，如果采用按列重构的方法，OMP的计算复杂度为 $O(N\times N(K_1^2+K_2^2+...+K_N^2\left)\right)=O\left(N^2\right(K_1^2+K_2^2+...+K_N^2\left)\right),$假设这里的K_i(i＝1,2,…,N) 相等或接近，并设K＝max_{i＝1,2,…,N}(K_i)，那么可以记为 O(N²K²)。

不同于OMP方法，WRFI利用原子与二维信号之间相关度的测量，在每次迭代时， 选出针对当前二维残差R(t-1):M×N的最大相关原子，并利用该原子更新矩阵A_t及直接 计算二维系数Θ的近似估计Θ_t:N×N，最终可以直接得到二维的最佳近似逼近解 而不是信号S某列s_i:N×1的重构信号从而将迭代次数减少为OMP 算法的1/N，即O(NK²)，有效地降低了计算复杂度。

FBWRFI方法进一步通过分块重构的方式降低了优化操作中的参数规模，降低了优 化方法的计算量。它设计每个分块大小B＝256，这样一个分块的计算复杂度为O(256K′²)， K′是分块的最大稀疏度。对图像S:N×N，其分块数位其计算复 杂度可表示为O(256nK′²)；且相对于OMP和WRFI方法的优化重建过程中使用的矩阵 和Φ:M×N，FBWRFI在优化重建中使用Ψ_B:256×256和Φ:M_B×256。也就是说，OMP 和WRFI的优化迭代中以N×N矩阵为计算对象，随着图像的增大，其计算时间的增长 率将以指数级增加；而FBWRFI的优化迭代中以256×256矩阵为计算对象，随着图像 的增加，只是相应增加了几个基本运算的循环，其计算时间的增长率理论上可以实现线 性增加。FBWRFI与现有压缩感知方法的计算复杂度对比如表4所示。

表4算法计算复杂度对比

FBWRFI方法效果

重建不同图像，本组实验验证FBWRFI方法对不同图像的作用效果，实验结果表6 和图4证明了算法的有效性(这里的测试图像大小为512×512)。

表6FBWRFI对不同图像的重构时间和重构质量

信号长度 重构时间(s) PSNR Lenna 16.35 27.07 Barbara 16.10 24.49 Peppers 16.74 27.47 Boat 16.48 24.89 Fingerprint 17.07 18.20 MRI 17.35 27.23

不同采样率下重建图像效果，本组实验验证不同采样率下FBWRFI方法重建图像的 效果，实验结果表7和图5证明了算法的有效性(这里的测试图像大小为512×512)。

表7不同采样率下FBWRFI重建图像的时间和质量

不同大小图像的重建效果，本组实验验证FBWRFI方法重建大同大小的图像的效 果，实验结果表8和图6证明了算法的有效性(这里的测试图像大小为256×256、 512×512和1024×1024)。

表8FBWRFI和OMP重建不同大小图像的时间和质量对比

快速重构算法的研究一直都是压缩感知理论研究的热点，我们提出面向图像重建的 快速压缩感知方法FBWRFI方法，它采用最小二乘优化方法实现图像的重建，基于图像 的整体相关度实现原子和二维残差之间的相关度测量以减少迭代次数、降低计算复杂度， 引入块压缩感知降低优化迭代中的参数规模以降低计算量，并使重构时间随信号增大而 呈线性增长，提高了重构时间的可预估性。

应当理解的是，本发明的上述具体实施方式仅仅用于示例性说明或解释 本发明的原理，而不构成对本发明的限制。因此，在不偏离本发明的精神和 范围的情况下所做的任何修改、等同替换、改进等，均应包含在本发明的保 护范围之内。此外，本发明所附权利要求旨在涵盖落入所附权利要求范围和 边界、或者这种范围和边界的等同形式内的全部变化和修改例。