张雁峰，范西岸，尹志益，蒋铁钢

(1.广东工业大学 机电工程学院，广州 510006； 2.广东科技学院 机电工程学院, 广东 东莞 523083)(*通信作者电子邮箱1648421565@qq.com)

0 引言

在信号处理领域，传统的奈奎斯特采样方法的采样频率须高于实际信号最高频率的2倍，如此高密度的采样方式给数据量日益增加的信号处理过程增加了压力。压缩感知(Compressed Sensing, CS)，是由Donoho[1]提出来的一种亚采样信息获取理论。跟传统的奈奎斯特采样理论不同，CS采样频率可以远小于实际信号最高频率的两倍，即只需少量采样便可精确地重构原信号，从而大大降低了信号的存储和运输的代价。目前压缩传感理论已经在信息论、医疗成像、模式识别、雷达探测、地质勘探、图像压缩、图像超分辨率重建等领域受到高度关注和应用[2]。信号重构算法关系着信号重构速度和质量，一直是CS的研究热点。

对于长度为N的信号f，可以表示为一组正交基向量Ψ={ψi|i=1,2,…，N}的线性组合，即:

f=Ψx

(1)

其中:Ψ∈RN×N称为稀疏基;x∈RN×1为f在Ψ变换域中的系数向量，如果x中只有s(s≪N)个元素不为零(或远大于零，而其他元素接近于零)，则称f是s-稀疏的。

用一个M×N(M

y=ΦΨx=Ax

(2)

其中A∈RM×N称为传感矩阵。文献[3]证明了当传感矩阵A满足s阶有限等距性质(Restricted Isometry Property, RIP)时，可以通过求解式(3)以很高的概率精确求得x：

(3)

s. t.y=Ax

求解式(3)的过程称为CS信号重构。贪婪算法作为一种CS重构算法，由于具有操作简单、计算复杂度低的特点而被广泛应用。正交匹配追踪(Orthogonal Matching Pursuit, OMP)算法是一种典型的贪婪算法[4]，OMP算法每次迭代挑选一个原子，当采样数量较大或稀疏度较大时，算法复杂度急剧增加。文献[5]使用共轭梯度法代替OMP算法中每次迭代所需的最小二乘运算，在稀疏度较大的情况下优势明显。文献[6-8]在OMP算法的基础上，每次迭代挑选多个原子，一定程度上减少了找到正确支撑集的时间，但是需要选择合适的参数。文献[9]提出的迭代硬阈值(Iterative Hard Thresholding, IHT)算法是一种基于l0范数的凸优化算法，算法复杂度低，但是容易受传感矩阵缩放影响导致不稳定。文献[10]使用利普希茨收敛条件将IHT算法中的固定步长转化为一个下降序列，提出正则化 IHT(Normalized IHT, NIHT)算法，算法复杂度稍有增加，但是迭代次数大幅减少。IHT/NIHT的缺点是有可能导致某些支撑被重复选择，从而导致其找到正确的支撑集需要较多的迭代次数[11]。文献[12]提出共轭梯度 IHT(Conjugate Gradient IHT, CGIHT)算法, 该算法使用共轭梯度法逼近最优解，减少了迭代次数，但是在实际过程中，由于误差的累积和传感矩阵的列原子不完全正交，导致前后迭代中寻优方向不共轭，当迭代到最优解附近时还是可能出现类似梯度下降法的振荡现象，导致算法收敛速度变慢。

文献[13-15]中算法引入了回溯的思想，降低原子被重复选择的概率，计算复杂度相对较低且重构精度高。其中文献[13]中的子空间追踪(Subspace Pursuit, SP)算法和文献[14]中的基于回溯的IHT算法(Backtracking-based IHT algorithm, BIHT)在每次迭代中取当前迭代的s个支撑与前一次迭代的s个支撑组成大小为2s的候选集，并将观测向量重新投影到候选集对应的矩阵列所张成的空间，从投影值中选出s个支撑作为新的支撑集。文献[15]算与SP算法/BIHT采用类似的回溯过程，不同的是其每次迭代从大小为3s的候选集中选出s个元素作为支撑集。SP算法和文献[15]算法都是以OMP算法为基础，且都采用回溯思想，两种算法的重构性能相当。BIHT以IHT算法为基础，通过定点迭代法求解凸优化问题式(4)来求稀疏解，与SP算法和文献[15]算法采用最小二乘求解相比，BIHT迭代次数更多。

s. t. ‖x‖0≤s

(4)

本文以BIHT为基础，保留其每次迭代过程中的回溯过程，用梯度下降法和共轭梯度法结合的方式代替BIHT中单纯使用梯度下降法作为寻优方法，提出了一种基于回溯的共轭梯度迭代硬阈值算法(Backtracking-based Conjugate Gradient Iterative Hard Thresholding algorithm, BCGIHT)。通过一维信号和二维图像的重构仿真实验，与NIHT算法、BIHT、CGIHT算法、文献[15]的算法作对比，验证了BCGIHT的性能。

1 基于回溯的共轭梯度迭代硬阈值算法

1.1 算法描述

BCGIHT主要策略是：一方面，结合BIHT的回溯思想，使得迭代过程中支撑集被快速准确找到，确保算法重构概率和重构精度；另一方面，采用梯度下降法和共轭梯度法相结合的方式代替迭代硬阈值类算法所使用的单纯的梯度下降法，进一步提高算法的收敛速度。

首先在每次迭代中采用非线性算子Hs(α)获得系数向量xn中绝对值最大的s个支撑，并取前一步迭代的支撑集与当前s个支撑的并集作为候选集，将观测向量y在候选集对应的矩阵列所张成的空间重新做一次正交投影得到中间系数向量an，再筛选出an中绝对值最大的s个支撑作为新的支撑集。由于每次迭代中回溯了前一步迭代中xn-1的支撑，减少了某些支撑在第n-1次迭代中被设为0而在第n次迭代中重新被选择的次数，使得支撑被反复选择的概率大大降低，优化了支撑集的选择，从而可以快速找到正确支撑集。

优化方法中梯度下降法具有线性收敛速度，共轭梯度法具有超线性收敛速度。由于误差的累积及实际过程中传感矩阵的列原子不完全正交，导致在迭代过程中，前后寻优方向并不共轭，因此使用单纯的共轭梯度法作为寻优方法效果不好。BCGIHT在BIHT单纯使用梯度下降法的基础上引入共轭梯度法加速算法的收敛。以前后迭代所得支撑集是否相等作为准则来确定寻优方向：若前后支撑集相等，即Γn=Γn-1，说明由前后支撑集所对应的观测矩阵的列组成的矩阵相同，即AΓn=AΓn-1，因此AΓnTAΓn=AΓn-1TAΓn-1，这时可以继续用共轭梯度方向进行下一步寻优；而当前后迭代中的支撑集不相等时，上述两个对称正定矩阵AΓnTAΓn和AΓn-1TAΓn-1不再相等，前后寻优方向dn和dn-1不再是共轭的，因此调整寻优方向为梯度下降方向。相较于单纯地使用梯度方向或者共轭方向的方法，寻优方向更替策略能在每次迭代中自动调整寻优方向，从而加速算法收敛。

1.2 算法流程

以下流程中:Hs(α)是一个非线性算子，表示取α中绝对值最大的s个元素;μ为迭代步长;Γn表示第n次迭代的索引集合;〈·,·〉表示求内积;supp(·)表示向量非零元素索引;ε表示迭代停止误差阈值。BCGIHT算法具体流程如下：

步骤1 初始化x0=0，残差r0=y，初始梯度g0=ATy，初始寻优方向d0=g0，初始支撑集Γ0=supp(Hs(g0))，μ0=〈g0,g0〉/〈Ag0,Ag0〉，x1=Hs(x0+μ0d0)，Γ1=supp (x1)。

步骤2n≥1时，更新梯度方向gn=AT(y-Axn)，判断是否Γn=Γn-1，如果是，转步骤4，否则，转步骤3。

步骤3 梯度下降法。计算步长μn=〈gn Γ n,gn Γ n〉/〈AΓngnΓ n,AΓ ngnΓn〉，计算an=Hs(xn+μngn)。

步骤4 共轭梯度法。计算方向因子βn=〈gn Γ n,gnΓ n〉/〈g(n-1)Γ n-1,g(n-1)Γ n-1〉,更新dn=gn+βndn-1、μn=〈dnΓ n,gnΓ n〉/〈AΓ ndnΓ n,AΓndnΓn〉，计算an=Hs(xn+μndn)。

步骤6 更新残差rn+1=y-Axn+1。

步骤7 判断是否‖rn+1‖2<ε，若是，迭代停止；否则，n=n+1，返回步骤2继续迭代。

2 算法性能分析

1)算法收敛性分析。BCGIHT采用梯度与共轭梯度交替的方式作为寻优方向，由IHT算法和共轭梯度法的收敛性证明，可知在每次迭代中，残差‖rn‖2<‖rn-1‖2，因此BCGIHT算法至少可以收敛到局部最优解。

表1 不同重构算法所需的迭代次数Tab. 1 Iterative number required for different algorithms

3 仿真结果与分析

为了验证BCGIHT的性能，本文使用一维信号和图像进行一系列仿真实验，并将该算法与具有代表性的NIHT算法、BIHT、CGIHT算法、文献[15]算法作仿真对比。所有的实验都是在Intel i5-7500、8 GB DDR4内存和Matlab 2014b上完成。

3.1 一维信号重构

本次实验使用的测试信号长度为N=256，使用M×N的随机高斯矩阵作为采样矩阵，采样数量M=128，每次实验重复做1 000次，以平均值作为实验结果。重构成功率是指对同一实验进行多次重复实验，其重构均方根误差在一定的误差阈值λ内的概率，本次实验将误差阈值设为λ=E-12。在不同稀疏度下，NIHT、CGIHT、BIHT、BCGIHT和文献[15]算法的重构成功率对比结果如图1所示。

从图1可以看出,当采样次数固定，稀疏度越来越大时，所有算法的重构成功率都会变小，这是由于稀疏度s变大，导致由s个大系数组成的支撑的不确定性增加导致的;NIHT算法在稀疏度比较小时，重构成功率能保持较高水平，但是当s>35时，重构成功率急剧下降;BIHT和CGIHT相较于NIHT稍好，BCGIHT相对于BIHT和CGIHT有较大提高，在s=45时还能保持较高的重构成功率;文献[15]算法在s>40时重构成功率急剧下降，这是由于稀疏度较大时，候选支撑集元素之间相互影响导致。

图1 不同稀疏度下不同算法重构成功率对比Fig. 1 Reconstruction success rate comparison of different algorithms with different sparsity

一维信号在不同稀疏度情况下不同算法重构时间对比如图2所示。由图2可知，随着稀疏度的增加，五种算法的重构时间都趋于增加，且文献[15]算法重构时间增加得越来越明显；BCGIHT重构时间少于BIHT和文献[15]算法25%以上,BCGIHT重构时间介于BIHT和CGIHT算法之间;CGIHT算法所用时间最少，这是因为本实验设置的一维信号是严格稀疏信号，稀疏基为单位矩阵，不存在原子之间的相互干扰，AΓnTAΓn为严格对称正定矩阵，因此CGIHT算法能以很少的迭代次数达到迭代停止条件。

图2 不同稀疏度下不同算法重构时间对比Fig. 2 Reconstruction time comparison of different algorithms with different sparsity

3.2 图像重构

本次实验采用常用的测试图像Pepper(像素为512×512)验证BCGIHT重构效果。首先，对原始图像进行小波变换，然后利用高斯随机矩阵作为观测矩阵进行压缩采样，最后使用不同的算法重构图像，并比较分析重构结果。

在采样率为0.5(即M=0.5×512=256)、稀疏度为s=M/4=64时,不同算法重构出的Pepper图像的直观效果如图3所示。其中图3(b)～(f)的峰值信噪比(Peak Signal-to-Noise Ratio, PSNR)分别为33.05 dB、36.34 dB、31.58 dB、36.15 dB、34.58 dB。

图3 不同算法的重构结果Fig. 3 Reconstruction results of different algorithms

为了比较在不同的采样率下不同算法重构质量和重构时间，实验设置了三组不同的采样率，分别为0.3、0.5和0.7。用PSNR反映重构质量，其计算式为：

PSNR=10lg((N-1)2/MSE);

(5)

每种采样率下不同的算法分别进行100次蒙特卡洛实验，计算平均值作为结果。不同采样率下各重构算法的重构PSNR和重构时间对比如表2所示。从表2中可知，在采样率为0.3和0.5时，BIHT和BCGIHT的PSNR略高于文献[15]算法;在采样率为0.7时，三种算法PSNR相当,相较而言，NIHT算法和CGIHT算法PSNR较小。在采样率小于等于0.5时，BCGIHT重构时间最短；CGIHT算法和BCGIHT的重构时间均小于BIHT，这是因为BIHT中有最小二乘计算，因此重构时间相对较长，而BCGIHT虽然也需最小二乘计算，但是因为寻优方向改进了，因此其迭代次数减少了，所以其重构时间不仅没有增加，反而减少了，且相较于BIHT，对于不同的采样率，BCGIHT重构时间减少50%以上；文献[15]算法重构时间介于BIHT和BCGIHT之间，符合理论分析。

在原图上加上不同方差的高斯白噪声时，不同算法重构带噪图像的重构精度对比如表3表示。对比表2和表3可知，相比于原图重构，CGIHT算法PSNR减少幅度较小，且噪声方差越大时表现越明显，这是因为CGIHT算法没有使用回溯方法挑选支撑集，不存在候选集原子之间相互干扰的情况;NIHT算法、BIHT、BCGIHT和文献[15]算法抗噪性能相当;噪声方差较大时，BCGIHT抗噪性能略好。

表2 不同采样率下不同算法的重构PSNR和重构时间Tab. 2 Reconstruction PSNR and time for different algorithms at different sample rates

表3 采样率为0.5时噪声环境下不同算法重构PSNRTab. 3 Reconstruction PSNR of different algorithms with Gaussian noise at the sample rate of 0.5

4 结语

本文在压缩感知BIHT的研究基础上，针对BIHT迭代次数多且重构时间长的问题，提出BCGIHT并对其进行了收敛性和算法复杂度分析:一方面，用回溯思想确保正确支撑集被快速找到;另一方面，采用梯度与共轭梯度结合使用的策略来逼近最优解。一维信号和二维图像的仿真实验结果表明BCGIHT在同精度条件下与BIHT及类似算法相比，迭代次数更少，重构时间远低于BIHT，兼顾了算法重构精度和重构效率,且在不同强度噪声环境下的重构性能与BIHT及同类算法相当，具有较好的抗噪性。