原子集校正及步长可控的稀疏度未知CS重构

曾春艳， 马丽红， 杜明辉

1.华南理工大学电子与信息学院，广州510641

2.湖北工业大学电气与电子工程学院，武汉430068

在压缩感知(compressive sensing,CS)理论中，假设x∈RN为一个k稀疏信号，Φ∈RM×N为测量矩阵，通过M个线性测量，可得到对x的一个观测信号y=Φx，其中Φ可看作一个字典，它的每一列被称作原子.如果k＜＜N且Φ满足约束等距性质(restricted isometry property,RIP)，则存在精确重建算法使x从M=O(k lb(N))个测量值中实现无失真恢复[1].

两种常用的CS重建方法分别基于l1范数最小和l0范数最小.前者利用线性规划求解，有着很高的计算复杂度O(N3)[2]；后者用迭代贪婪追踪策略近似信号[3]，其每次迭代包含最佳原子选择、信号估计、残差更新3步.匹配追踪(matching pursuits,MP)[4]和正交匹配追踪(orthogonal matching pursuit,OMP)[5]算法通过将残差与最相关原子匹配来重建信号，然而它们的解有时会被错误原子干扰.为剪切错误选择的原子，在已知稀疏度下，可引入回溯策略.子空间追踪(subspace pursuit,SP)[6]与压缩感知匹配追踪(compressed sensing matching pursuit,CoSaMP)[7]都是已知稀疏度k的回溯型重建算法，它们在原子选择环节分别保留k和2k个原子.另一种回溯算法是迭代硬阈值(iterative hard thresholding,IHT)[8]算法，它每次迭代保留近似信号幅值最大的k个元素xn=Hk(xn-1+ΦT(y-Φxn-1)).然而在实际应用中稀疏度通常是未知的，为此回溯型自适应正交匹配算法(backtracking-based adaptive orthogonal matching pursuit,BAOMP)[9]使用了一种基于阈值选择的方法，它以固定的权重参数代替对稀疏度的估计来确定增减原子阈值，从而实现盲稀疏度下对原信号的重建[10].但BAOMP算法仍存在明显缺陷：1)信号稀疏度k与阈值之间的关系没有分析；2)经过多次迭代后，搜索残差会从相对小的值突然反弹到大的值，并且准周期地反复，降低了算法的收敛速度.

本文首先对残差信号用类高斯分布建模，分析了BAOMP算法的阈值选择方法与稀疏度已知算法使用稀疏度的一致性与差异点，然后提出改进的回溯型自适应正交匹配算法(improved backtracking-based adaptiveorthogonal matching pursuit,IBAOMP)，它采用80–20准则判断解信号的匹配程度，并在后续匹配步骤中引入可变步长阈值，精细调整原子集容量，提高了选入原子的正确匹配率，从而避免了残差信号的准周期性失配.IBAOMP算法增加了2个参量：1)匹配准确度阈值T，T决定新选原子集容量何时被调整；2)增量步长S，S定义了阈值权重的调整步长.在第n次迭代时，若上次迭代残差模|rn-1|＞T，则原子匹配策略与BAOMP算法相同；若|rn-1|≤T，则选入阈值权重µ1以步长S增大到µ1+S，减少了新选原子集容量，使后续迭代步骤选择的原子数目减少，从而能充分搜索出更接近原始信号的解.

1 BAOM P算法中原子匹配与准周期性失配

因为最终原子匹配集要包含k个原子，所以稀疏度k对精确重建十分关键，然而大多数情况下k为估计值.BAOMP算法通过设置阈值权重参数µ1和µ2来控制新选原子与淘汰原子数目，使重建信号近似k稀疏[9]，选择与淘汰原子策略如下：

1)新选原子集索引Cn

当原子ϕi与残差内积满足以下关系时，原子索引i被选入集合Cn

式中，〈·〉表示内积运算，|·|表示取绝对值，i,j∈{1,2,···,N}.µ1度量了上次迭代残差rn-1与原子ϕi内积的最大值，较小的µ1将在每次搜索时选择更多原子，因此加快了前期迭代的匹配速度，但会使迭代后期选入更多错误原子.

2)淘汰原子集索引Γn

对新的候选原子索引集Ωn=Λn-1∪Cn，其中Λn-1为上次迭代近似信号支撑集，计算这些候选原子的近似系数

式中，“†”表示矩阵的伪逆.当Ωn中对应原子的近似系数小于淘汰阈值时，该索引值加入淘汰原子索引集Γn，即

式中，i∈Ω，j∈Cn.µ2控制了淘汰原子数目.

下面本文通过对残差信号建模分析BAOMP算法的阈值选择方法与已知稀疏度k方法的关系，以及固定权重µ1导致的准周期性失配现象.

1)稀疏度k与阈值的关系

BAOMP算法通过比较所有原子残差内积与µ1倍最大值来选择新原子，从而用已选原子逼近上次迭代残差.因为观测阵的元素通常是类高斯分布(如高斯随机矩阵、部分傅里叶矩阵、Chirp感知矩阵)，即使信号x是混合分布，残差rn=y-也是类高斯分布.本文模拟超过300个残差信号后确认它们的分布如下：

式中，r为残差r的分量，u为均值，δ为标准差，a为常数.因此，搜索k个原子等价于通过搜寻一个正阈值t，使图1中黑色区域与k成比例，从而逼近残差，即

图1 残差信号幅值类高斯分布Figure 1 Gaussian-like distributed amplitudes of residuals

2)逼近稀疏度

对于类高斯分布，k可用µ1近似

式中，rmax是观测信号或残差的最大绝对值.k应该由µ1和理论上可收缩的rmax近似得到，但实际上µ1为固定值，且rmax可能因错误匹配而增大，于是对k的估计可能存在较大误差.

3)准周期性失配

尽管对k的估计可转换成一个阈值寻找问题，k依旧是未知的.当残差模值较小时，固定的µ1会导致在如式(1)所示的全局搜索时选入错误原子，并导致残差模值剧烈反弹.图2为图像“Lena”在第34次和第35次迭代后的残差，第34次迭代残差的幅值范围为[–1,+1]，但一次错误匹配迭代后，第35次迭代的残差幅值范围上升至[–71,+75].这两次迭代的残差分布如图3所示.另一个高斯稀疏信号残差模值|rn|准周期性失配的例子如图4所示，其中稀疏度k=12，观测数目M=45，信号长度N=256，第6次迭代后残差2范数|rn|下降到3.6，但由于错误匹配而导致第7次迭代后显著上升至20.4.在极端情况下，当达到最大迭代数时，仍不能精确重建.

图2 Lena图像用BAOMP算法错误匹配后残差反弹Figure 2 Rebounded residuals of Lena after a mismatch by BAOMP

图3 Lena图像用BAOMP算法错误匹配前后残差幅值分布直方图Figure 3 Histograms of amplitudes of residuals of Lena after a mismatch by BAOMP

2 原子集校正及容量控制

为避免上述稀疏度k未知的BAOMP算法中固定的µ1引起的准周期性失配现象，本文调整搜索方案去逼近稀疏度k.改进的回溯型自适应正交匹配算法(IBAOMP)定义了两个参量：匹配准确度阈值T和增量步长S.与BAOMP中固定的µ1不同，T和S使µ1更加精准地匹配残差，当残差模值|rn|变小时，选入的原子数目也随之变少.具体调整如下：

图4 高斯稀疏信号用BAOMP算法恢复时残差范数随迭代次数变化Figure 4 Residual norms versus the numbers of iterations of a Gaussian sparse signal by BAOMP

1)用可变µ1精确匹配

引入匹配准确度阈值T，通过比较残差模值|rn|确定µ1调整时刻.

当|rn|＞T时，µ1保持不变来加速匹配；当|rn|≤T时，大多数正确原子已经搜寻到，增加µ1来避免错误投影与准周期性失配.T的合理取值为

根据RIP原则，当约束等距常数δ≪1时，||y||2≈||x||2，观测信号y的模值可以代表x的能量.又根据80–20准则，在大多数情况下，稀疏空间中不到20%的大幅值元素拥有信号80%的能量，因此式(8)是可行的.

2)增量步长S控制

µ1以步长S增加，如式(9)所示，S越大候选原子集容量越小

图5显示了用IBAOMP算法重建高斯稀疏信号时残差模值|rn|随迭代次数的变化，其中T=15，S=0.01，其他参数与图4相同.作为对比，前3次类三角周期反弹并未被抑制.图5中第18次迭代后残差模值控制在|rn|＜10-4，而图4中，BAOMP算法中的残差模在第18次迭代后又反弹了4次.

IBAOMP算法的实现步骤如下：

步骤1 迭代前期新选原子集索引.当|rn|＞T时，与BAOMP算法相同产生新选原子集索引Cn，如式(1)所示.

图5 高斯稀疏信号用IBAOMP算法恢复时残差范数随迭代次数变化Figure 5 Residual norms versus the numbers of iterations of a Gaussian sparse signal by IBAOMP

步骤2 迭代后期校正与步长控制.当|rn|≤T时，校正µ1值至µ1+S，如式(9)所示.

步骤3 信号近似.产生新的候选原子索引集Ωn=Λn-1∪Cn，由等式(2)计算近似系数xΩn

步骤4 原子淘汰.由式(3)产生淘汰原子集索引Γn.

步骤5 更新支撑集索引Λn.更新近似系数xnΛ与残差rn如下：

在IBAOMP算法的一次迭代中，计算N个内积与锁定Cn分别需要O(M N)和O(N)个标准乘，最小二乘计算和分别需要O(|Ω|N)和O(|Λ|N)个标准乘(其中|Λ|≤|Ω|＜N).因此，与BAOMP算法相同，IBAOMP算法的一次迭代需要的计算量为O(M N)个标准乘.

3 仿真实验及结果分析

为验证IBAOMP算法性能，进行以下2组实验：1)匹配准确度阈值T和增量步长S性能测试；2)比较IBAOMP算法与BAOMP、OMP、CoSaMP、St OMP、ROMP算法在不同稀疏度k和不同测量数M下精确重建概率.k稀疏信号为非零元素服从标准高斯分布的稀疏信号，k=12，长度N=256，观测矩阵为归一化高斯矩阵.可压缩信号为3幅256×256图像：频率主要在低频的“Peppers”图像、低频到高频都有的“Lena”图像、纹理丰富的高频“Baboon”图像，这3幅图像的稀疏度依次由低到高，对于可压缩信号采用部分傅里叶矩阵作为观测矩阵.所有实验重复500次，在某一次重建实验中，对稀疏信号有‖x-‖2＜10-4，对可压缩信号有‖x-‖2＜2时，该次重建可看作精确重建.定义精确重建次数Nr与总次数N的比率为精确重建概率[11]

3.1 参数T和S对精确重建结果的影响

首先观察当S值固定，T均匀增大时，匹配准确度阈值T对信号精确重建性能的影响，T=0时的IBAOMP算法等同于BAOMP算法.然后选择固定的T值，步长S均匀增大时，观察S对信号精确重建性能的影响.实验首先用稀疏信号得到明确结果，然后用可压缩信号测试实际应用性能.

3.1.1 T和对S高斯稀疏信号重建性能的影响

在BAOMP算法中，µ1取0.4时精确重建效果最佳[9]，参考这一结论将实验参数设置如下：S=0.01、M=45、µ1=0.4、µ2=0.6，T以步长0.5由0增加到20，精确重建概率p和平均迭代结束时µ1的终值如图6所示.

图6 S=0.01时匹配准确度阈值T对精确重建概率的影响Figure 6 Matching accuracy threshold T versus probability of exact reconstruction with S=0.01

从图6中可以看出：随着T的增加，µ1的终值增加，对应的新增匹配原子集容量减少.T由0增加到20时，IBAOMP算法的平均精确重建概率为82.4%，比BAOMP算法的69%高13.4%.当T取19.5时，IBAOMP算法的精确重建概率最高，达到89.6%，比BAOMP算法提高20.6%，说明对于该高斯信号和步长S的设定，Tœ取值使精确重建概率达到最高的19.5时，µ1在合适的时刻被调整，最终能将精确重建概率提高到89.6%.

固定T=12、µ1=0.4，S以步长0.001由0增加到0.02，实验结果如图7所示.图7表明S取较小值时，后段迭代对µ1的调整较小，选取原子数目随之平稳减少，更适合该信号结构特征，进而有更高的精确重建概率.当S=0.003时，最终µ1=0.47已经自适应该信号结构，使得最终精确重建概率比BAOMP算法高出10.2%.

图7 T=12时增量步长S对精确重建概率的影响Figure 7 Incremental step S versus probability of exact reconstruction with T=12

3.1.2 T对可压缩信号的影响

本文验证了重建图像“Peppers”、“Lena”和“Baboon”时新增原子数目和淘汰原子数目随迭代次数的变化.图8为“Lena”图像结果，其余两幅图结果类似，在此未给出.

由图8(a)可看出，用BAOMP算法重建图像，当迭代次数为35、66、97时新增原子数目发生剧烈变化，而此时对应的淘汰原子数目(见图8(b))也发生剧烈变化，表明这些迭代发生了错误匹配.

图9给出在参数T作用下对不同稀疏度信号重建时残差模值随迭代次数变化情况.图9(a)和9(b)中“Peppers”与“Lena”图像80%的能量由少量系数控制，IBAOMP算法明显抑制了残差模值的反弹.而对于图9(c)中稀疏度较大的细节信号“Baboon”图像，在选入错误原子使得残差模值反弹前，需要更多的迭代来达到初始匹配.虽然需要较长的匹配过程，式(8)依旧有效，图9中曲线拐点都出现在位置T=0.2norm(y).

3.2 IBAOM P算法与其他贪婪算法的比较

参与IBAOMP算法比较的有BAOMP、OMP、CoSaMP、StOMP和ROMP算法.BAOMP中µ1=0.4，µ2=0.6，IBAOMP中T=0.2norm(y)，S=0.01.OMP和St OMP使用SparseLab工具箱实现，OMP算法最大迭代数为k，StOMP算法最大迭代数为10，误差容限为10-5，CoSaMP和ROMP使用作者附录中给出的程序，最大迭代数为8k，误差容限10-4.

图8 Lena图像每次迭代新增及淘汰原子数目Figure 8 Numbers of adding and deleting atoms in each iteration of Lena

图9 不同稀疏度图像残差模值反弹Figure 9 Rebounded residuals norms of pictures with different sparsity levels

3.2.1 不同测量数M下IBAOM P算法与其他重构算法的精确重建性能比较

对k=12高斯稀疏信号用IBAOMP算法及其他典型贪婪算法重建结果如图10所示.当稀疏度保持不变时，随着测量数的增加，所有贪婪算法精确重建概率逐步上升.当测量数M＜75时，由于测量数M过少，各类算法都不能以概率1完全重建原始信号，但IBAOMP算法比BAOMP算法的精确重建概率高出26%；当测量数M＞75时，这两种算法都以概率1精确重建原信号.CoSaMP和OMP算法在测量数M＞80时的精确重建概率达到1，而StOMP和ROMP算法在测量数M=120时的精确重建概率分别为70.6%和98.2%.

图10 不同测量数下IBAOMP算法与其他重构算法比较Figure 10 Comparing IBAOMP with other algorithms with different numbers of measurement

3.2.2 不同稀疏度k下IBAOM P算法与其他重构算法的精确重建性能比较

对M=148高斯稀疏信号用IBAOMP算法及其他典型贪婪算法重建结果如图11所示.随着稀疏度k的增加，需要的测量数更多，当测量数M=148相对于稀疏度k不够时，精确重建概率开始下降.当稀疏度k＜70时，IBAOMP算法以概率1精确重建；当稀疏度k＞70时，IBAOMP算法比BAOMP算法的精确重建概率下降得更慢；当稀疏度k=80时，IBAOMP算法比BAOMP算法的精确重建概率高出17%.CoSaMP和OMP算法分别在稀疏度k＞50和k＞35时的精确重建概率迅速下降到零.StOMP算法在本实验中性能最差，当稀疏度k为5时，精确重建概率仅为75%.

图11 不同稀疏度下IBAOMP算法与其他重构算法比较Figure 11 Comparing IBAOMP with other algorithms with different sparsity levels

3.2.3 不同稀疏度图像算法对比

分别用IBAOMP算法与BAOMP算法对稀疏度不同的图像“Peppers”、“Lena”、“Baboon”进行重建.为快速达到初步匹配，对于“Baboon”，取µ1=0.2；对于“Peppers”和“Lena”，取µ1=0.3.重建结果如表1所示，可见IBAOMP算法明显优于BAOMP算法.

表1 不同稀疏度图像精确重建概率比较Table 1 Comparing probabilities of exact reconstruction of pictures with different sparsity levels

4 结语

本文提出的IBAOMP算法适用于压缩感知中对稀疏度未知的信号重建.匹配准确度阈值T触发对µ1的调整，避免了准周期性失配；增量步长S调整了匹配原子集容量.实验表明IBAOMP算法解决了固定µ1可能带来的原子错误匹配及残差模值增加的问题，从而在计算复杂度不变的前提下提高了精确重建概率.

[1]GIRYES R,ELAD M.RIP-based near-oracle performance guarantees for SP,CoSaMP,and IHT[J].IEEE Transactions on Signal Processing,2012,60(3):1465-1468.

[2]YANG Jingyu,PENG Yigang,XU Wenli,DAI Qionghai.Ways to sparse representation:a comparative study[J].Tsinghua Science and Technology,2009,14(4):434-443.

[3]TROPP J A,WRIGHT S J.Computational methods for sparse solution of linear inverse problems[J].Proceedings of the IEEE,2010,98(6):948-958.

[4]MALLAT S G,ZHANG Z F.Matching pursuits with time-frequency dictionaries[J].IEEE Transactions on Signal Processing,1993,41(12):3397-3415.

[5]TROPPJ A,GILBERT A C.Signal recovery from random measurements via orthogonal matching pursuit[J].IEEE Transactions on Information Theory,2007,53(12):4655-4666.

[6]VARADARAJAN B,KHUDANPUR S,TRAN T D.Stepwise optimal subspace pursuit for improving sparse recovery[J].IEEE Signal Processing Letters,2011,18(1):27-30.

[7]NEEDELL D,TROPP J A.CoSaMP:iterative signal recovery from incomplete and inaccurate samples[J].Applied and Computational Harmonic Analysis,2009,26(3):301-321.

[8]BLUMENSATHT,DAVIESM E.Iterative hard thresholding for compressed sensing[J].Applied and Computational Harmonic Analysis,2009,27(3):265-274.

[9]HUANG H L,MAKUR A.Backtracking-based matching pursuit method for sparse signal reconstruction[J].IEEE Signal Processing Letters,2011,18(7):391-394.

[10]GURBUZ A C,PILANCI M,ARIKAN O.Expectation maximization based matching pursuit[C]//IEEE International Conference on Acoustics,Speech and Signal Processing,2012:3313-3316.

[11]DAI W,MILENKOVIC O.Subspace pursuit for compressive sensing signal reconstruction [J].IEEE Transactions on Information Theory,2009,55(5):2230-2249.