曾春艷， 馬麗紅， 杜明輝

1.華南理工大學電子與信息學院，廣州510641

2.湖北工業(yè)大學電氣與電子工程學院，武漢430068

在壓縮感知(compressive sensing,CS)理論中，假設(shè)x∈RN為一個k稀疏信號，Φ∈RM×N為測量矩陣，通過M個線性測量，可得到對x的一個觀測信號y=Φx，其中Φ可看作一個字典，它的每一列被稱作原子.如果k＜＜N且Φ滿足約束等距性質(zhì)(restricted isometry property,RIP)，則存在精確重建算法使x從M=O(k lb(N))個測量值中實現(xiàn)無失真恢復(fù)[1].

兩種常用的CS重建方法分別基于l1范數(shù)最小和l0范數(shù)最小.前者利用線性規(guī)劃求解，有著很高的計算復(fù)雜度O(N3)[2]；后者用迭代貪婪追蹤策略近似信號[3]，其每次迭代包含最佳原子選擇、信號估計、殘差更新3步.匹配追蹤(matching pursuits,MP)[4]和正交匹配追蹤(orthogonal matching pursuit,OMP)[5]算法通過將殘差與最相關(guān)原子匹配來重建信號，然而它們的解有時會被錯誤原子干擾.為剪切錯誤選擇的原子，在已知稀疏度下，可引入回溯策略.子空間追蹤(subspace pursuit,SP)[6]與壓縮感知匹配追蹤(compressed sensing matching pursuit,CoSaMP)[7]都是已知稀疏度k的回溯型重建算法，它們在原子選擇環(huán)節(jié)分別保留k和2k個原子.另一種回溯算法是迭代硬閾值(iterative hard thresholding,IHT)[8]算法，它每次迭代保留近似信號幅值最大的k個元素xn=Hk(xn-1+ΦT(y-Φxn-1)).然而在實際應(yīng)用中稀疏度通常是未知的，為此回溯型自適應(yīng)正交匹配算法(backtracking-based adaptive orthogonal matching pursuit,BAOMP)[9]使用了一種基于閾值選擇的方法，它以固定的權(quán)重參數(shù)代替對稀疏度的估計來確定增減原子閾值，從而實現(xiàn)盲稀疏度下對原信號的重建[10].但BAOMP算法仍存在明顯缺陷：1)信號稀疏度k與閾值之間的關(guān)系沒有分析；2)經(jīng)過多次迭代后，搜索殘差會從相對小的值突然反彈到大的值，并且準周期地反復(fù)，降低了算法的收斂速度.

本文首先對殘差信號用類高斯分布建模，分析了BAOMP算法的閾值選擇方法與稀疏度已知算法使用稀疏度的一致性與差異點，然后提出改進的回溯型自適應(yīng)正交匹配算法(improved backtracking-based adaptiveorthogonal matching pursuit,IBAOMP)，它采用80–20準則判斷解信號的匹配程度，并在后續(xù)匹配步驟中引入可變步長閾值，精細調(diào)整原子集容量，提高了選入原子的正確匹配率，從而避免了殘差信號的準周期性失配.IBAOMP算法增加了2個參量：1)匹配準確度閾值T，T決定新選原子集容量何時被調(diào)整；2)增量步長S，S定義了閾值權(quán)重的調(diào)整步長.在第n次迭代時，若上次迭代殘差模|rn-1|＞T，則原子匹配策略與BAOMP算法相同；若|rn-1|≤T，則選入閾值權(quán)重μ1以步長S增大到μ1+S，減少了新選原子集容量，使后續(xù)迭代步驟選擇的原子數(shù)目減少，從而能充分搜索出更接近原始信號的解.

1 BAOM P算法中原子匹配與準周期性失配

因為最終原子匹配集要包含k個原子，所以稀疏度k對精確重建十分關(guān)鍵，然而大多數(shù)情況下k為估計值.BAOMP算法通過設(shè)置閾值權(quán)重參數(shù)μ1和μ2來控制新選原子與淘汰原子數(shù)目，使重建信號近似k稀疏[9]，選擇與淘汰原子策略如下：

1)新選原子集索引Cn

當原子?i與殘差內(nèi)積滿足以下關(guān)系時，原子索引i被選入集合Cn

式中，〈·〉表示內(nèi)積運算，|·|表示取絕對值，i,j∈{1,2,···,N}.μ1度量了上次迭代殘差rn-1與原子?i內(nèi)積的最大值，較小的μ1將在每次搜索時選擇更多原子，因此加快了前期迭代的匹配速度，但會使迭代后期選入更多錯誤原子.

2)淘汰原子集索引Γn

對新的候選原子索引集?n=Λn-1∪Cn，其中Λn-1為上次迭代近似信號支撐集，計算這些候選原子的近似系數(shù)

式中，“?”表示矩陣的偽逆.當?n中對應(yīng)原子的近似系數(shù)小于淘汰閾值時，該索引值加入淘汰原子索引集Γn，即

式中，i∈?，j∈Cn.μ2控制了淘汰原子數(shù)目.

下面本文通過對殘差信號建模分析BAOMP算法的閾值選擇方法與已知稀疏度k方法的關(guān)系，以及固定權(quán)重μ1導致的準周期性失配現(xiàn)象.

1)稀疏度k與閾值的關(guān)系

BAOMP算法通過比較所有原子殘差內(nèi)積與μ1倍最大值來選擇新原子，從而用已選原子逼近上次迭代殘差.因為觀測陣的元素通常是類高斯分布(如高斯隨機矩陣、部分傅里葉矩陣、Chirp感知矩陣)，即使信號x是混合分布，殘差rn=y-也是類高斯分布.本文模擬超過300個殘差信號后確認它們的分布如下：

式中，r為殘差r的分量，u為均值，δ為標準差，a為常數(shù).因此，搜索k個原子等價于通過搜尋一個正閾值t，使圖1中黑色區(qū)域與k成比例，從而逼近殘差，即

圖1 殘差信號幅值類高斯分布Figure 1 Gaussian-like distributed amplitudes of residuals

2)逼近稀疏度

對于類高斯分布，k可用μ1近似

式中，rmax是觀測信號或殘差的最大絕對值.k應(yīng)該由μ1和理論上可收縮的rmax近似得到，但實際上μ1為固定值，且rmax可能因錯誤匹配而增大，于是對k的估計可能存在較大誤差.

3)準周期性失配

盡管對k的估計可轉(zhuǎn)換成一個閾值尋找問題，k依舊是未知的.當殘差模值較小時，固定的μ1會導致在如式(1)所示的全局搜索時選入錯誤原子，并導致殘差模值劇烈反彈.圖2為圖像“Lena”在第34次和第35次迭代后的殘差，第34次迭代殘差的幅值范圍為[–1,+1]，但一次錯誤匹配迭代后，第35次迭代的殘差幅值范圍上升至[–71,+75].這兩次迭代的殘差分布如圖3所示.另一個高斯稀疏信號殘差模值|rn|準周期性失配的例子如圖4所示，其中稀疏度k=12，觀測數(shù)目M=45，信號長度N=256，第6次迭代后殘差2范數(shù)|rn|下降到3.6，但由于錯誤匹配而導致第7次迭代后顯著上升至20.4.在極端情況下，當達到最大迭代數(shù)時，仍不能精確重建.

圖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引起的準周期性失配現(xiàn)象，本文調(diào)整搜索方案去逼近稀疏度k.改進的回溯型自適應(yīng)正交匹配算法(IBAOMP)定義了兩個參量：匹配準確度閾值T和增量步長S.與BAOMP中固定的μ1不同，T和S使μ1更加精準地匹配殘差，當殘差模值|rn|變小時，選入的原子數(shù)目也隨之變少.具體調(diào)整如下：

圖4 高斯稀疏信號用BAOMP算法恢復(fù)時殘差范數(shù)隨迭代次數(shù)變化Figure 4 Residual norms versus the numbers of iterations of a Gaussian sparse signal by BAOMP

1)用可變μ1精確匹配

引入匹配準確度閾值T，通過比較殘差模值|rn|確定μ1調(diào)整時刻.

當|rn|＞T時，μ1保持不變來加速匹配；當|rn|≤T時，大多數(shù)正確原子已經(jīng)搜尋到，增加μ1來避免錯誤投影與準周期性失配.T的合理取值為

根據(jù)RIP原則，當約束等距常數(shù)δ?1時，||y||2≈||x||2，觀測信號y的模值可以代表x的能量.又根據(jù)80–20準則，在大多數(shù)情況下，稀疏空間中不到20%的大幅值元素擁有信號80%的能量，因此式(8)是可行的.

2)增量步長S控制

μ1以步長S增加，如式(9)所示，S越大候選原子集容量越小

圖5顯示了用IBAOMP算法重建高斯稀疏信號時殘差模值|rn|隨迭代次數(shù)的變化，其中T=15，S=0.01，其他參數(shù)與圖4相同.作為對比，前3次類三角周期反彈并未被抑制.圖5中第18次迭代后殘差模值控制在|rn|＜10-4，而圖4中，BAOMP算法中的殘差模在第18次迭代后又反彈了4次.

IBAOMP算法的實現(xiàn)步驟如下：

步驟1 迭代前期新選原子集索引.當|rn|＞T時，與BAOMP算法相同產(chǎn)生新選原子集索引Cn，如式(1)所示.

圖5 高斯稀疏信號用IBAOMP算法恢復(fù)時殘差范數(shù)隨迭代次數(shù)變化Figure 5 Residual norms versus the numbers of iterations of a Gaussian sparse signal by IBAOMP

步驟2 迭代后期校正與步長控制.當|rn|≤T時，校正μ1值至μ1+S，如式(9)所示.

步驟3 信號近似.產(chǎn)生新的候選原子索引集?n=Λn-1∪Cn，由等式(2)計算近似系數(shù)xΩn

步驟4 原子淘汰.由式(3)產(chǎn)生淘汰原子集索引Γn.

步驟5 更新支撐集索引Λn.更新近似系數(shù)xnΛ與殘差rn如下：

在IBAOMP算法的一次迭代中，計算N個內(nèi)積與鎖定Cn分別需要O(M N)和O(N)個標準乘，最小二乘計算和分別需要O(|?|N)和O(|Λ|N)個標準乘(其中|Λ|≤|?|＜N).因此，與BAOMP算法相同，IBAOMP算法的一次迭代需要的計算量為O(M N)個標準乘.

3 仿真實驗及結(jié)果分析

為驗證IBAOMP算法性能，進行以下2組實驗：1)匹配準確度閾值T和增量步長S性能測試；2)比較IBAOMP算法與BAOMP、OMP、CoSaMP、St OMP、ROMP算法在不同稀疏度k和不同測量數(shù)M下精確重建概率.k稀疏信號為非零元素服從標準高斯分布的稀疏信號，k=12，長度N=256，觀測矩陣為歸一化高斯矩陣.可壓縮信號為3幅256×256圖像：頻率主要在低頻的“Peppers”圖像、低頻到高頻都有的“Lena”圖像、紋理豐富的高頻“Baboon”圖像，這3幅圖像的稀疏度依次由低到高，對于可壓縮信號采用部分傅里葉矩陣作為觀測矩陣.所有實驗重復(fù)500次，在某一次重建實驗中，對稀疏信號有‖x-‖2＜10-4，對可壓縮信號有‖x-‖2＜2時，該次重建可看作精確重建.定義精確重建次數(shù)Nr與總次數(shù)N的比率為精確重建概率[11]

3.1 參數(shù)T和S對精確重建結(jié)果的影響

首先觀察當S值固定，T均勻增大時，匹配準確度閾值T對信號精確重建性能的影響，T=0時的IBAOMP算法等同于BAOMP算法.然后選擇固定的T值，步長S均勻增大時，觀察S對信號精確重建性能的影響.實驗首先用稀疏信號得到明確結(jié)果，然后用可壓縮信號測試實際應(yīng)用性能.

3.1.1 T和對S高斯稀疏信號重建性能的影響

在BAOMP算法中，μ1取0.4時精確重建效果最佳[9]，參考這一結(jié)論將實驗參數(shù)設(shè)置如下：S=0.01、M=45、μ1=0.4、μ2=0.6，T以步長0.5由0增加到20，精確重建概率p和平均迭代結(jié)束時μ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的終值增加，對應(yīng)的新增匹配原子集容量減少.T由0增加到20時，IBAOMP算法的平均精確重建概率為82.4%，比BAOMP算法的69%高13.4%.當T取19.5時，IBAOMP算法的精確重建概率最高，達到89.6%，比BAOMP算法提高20.6%，說明對于該高斯信號和步長S的設(shè)定，T?取值使精確重建概率達到最高的19.5時，μ1在合適的時刻被調(diào)整，最終能將精確重建概率提高到89.6%.

固定T=12、μ1=0.4，S以步長0.001由0增加到0.02，實驗結(jié)果如圖7所示.圖7表明S取較小值時，后段迭代對μ1的調(diào)整較小，選取原子數(shù)目隨之平穩(wěn)減少，更適合該信號結(jié)構(gòu)特征，進而有更高的精確重建概率.當S=0.003時，最終μ1=0.47已經(jīng)自適應(yīng)該信號結(jié)構(gòu)，使得最終精確重建概率比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”時新增原子數(shù)目和淘汰原子數(shù)目隨迭代次數(shù)的變化.圖8為“Lena”圖像結(jié)果，其余兩幅圖結(jié)果類似，在此未給出.

由圖8(a)可看出，用BAOMP算法重建圖像，當?shù)螖?shù)為35、66、97時新增原子數(shù)目發(fā)生劇烈變化，而此時對應(yīng)的淘汰原子數(shù)目(見圖8(b))也發(fā)生劇烈變化，表明這些迭代發(fā)生了錯誤匹配.

圖9給出在參數(shù)T作用下對不同稀疏度信號重建時殘差模值隨迭代次數(shù)變化情況.圖9(a)和9(b)中“Peppers”與“Lena”圖像80%的能量由少量系數(shù)控制，IBAOMP算法明顯抑制了殘差模值的反彈.而對于圖9(c)中稀疏度較大的細節(jié)信號“Baboon”圖像，在選入錯誤原子使得殘差模值反彈前，需要更多的迭代來達到初始匹配.雖然需要較長的匹配過程，式(8)依舊有效，圖9中曲線拐點都出現(xiàn)在位置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工具箱實現(xiàn)，OMP算法最大迭代數(shù)為k，StOMP算法最大迭代數(shù)為10，誤差容限為10-5，CoSaMP和ROMP使用作者附錄中給出的程序，最大迭代數(shù)為8k，誤差容限10-4.

圖8 Lena圖像每次迭代新增及淘汰原子數(shù)目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 不同測量數(shù)M下IBAOM P算法與其他重構(gòu)算法的精確重建性能比較

對k=12高斯稀疏信號用IBAOMP算法及其他典型貪婪算法重建結(jié)果如圖10所示.當稀疏度保持不變時，隨著測量數(shù)的增加，所有貪婪算法精確重建概率逐步上升.當測量數(shù)M＜75時，由于測量數(shù)M過少，各類算法都不能以概率1完全重建原始信號，但IBAOMP算法比BAOMP算法的精確重建概率高出26%；當測量數(shù)M＞75時，這兩種算法都以概率1精確重建原信號.CoSaMP和OMP算法在測量數(shù)M＞80時的精確重建概率達到1，而StOMP和ROMP算法在測量數(shù)M=120時的精確重建概率分別為70.6%和98.2%.

圖10 不同測量數(shù)下IBAOMP算法與其他重構(gòu)算法比較Figure 10 Comparing IBAOMP with other algorithms with different numbers of measurement

3.2.2 不同稀疏度k下IBAOM P算法與其他重構(gòu)算法的精確重建性能比較

對M=148高斯稀疏信號用IBAOMP算法及其他典型貪婪算法重建結(jié)果如圖11所示.隨著稀疏度k的增加，需要的測量數(shù)更多，當測量數(shù)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算法與其他重構(gòu)算法比較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.重建結(jié)果如表1所示，可見IBAOMP算法明顯優(yōu)于BAOMP算法.

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

4 結(jié)語

本文提出的IBAOMP算法適用于壓縮感知中對稀疏度未知的信號重建.匹配準確度閾值T觸發(fā)對μ1的調(diào)整，避免了準周期性失配；增量步長S調(diào)整了匹配原子集容量.實驗表明IBAOMP算法解決了固定μ1可能帶來的原子錯誤匹配及殘差模值增加的問題，從而在計算復(fù)雜度不變的前提下提高了精確重建概率.

[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.