曾春艷, 馬麗紅, 杜明輝
1.華南理工大學(xué)電子與信息學(xué)院,廣州510641
2.湖北工業(yè)大學(xué)電氣與電子工程學(xué)院,武漢430068
在壓縮感知(compressive sensing,CS)理論中,假設(shè)x∈RN為一個(gè)k稀疏信號(hào),Φ∈RM×N為測(cè)量矩陣,通過(guò)M個(gè)線性測(cè)量,可得到對(duì)x的一個(gè)觀測(cè)信號(hào)y=Φx,其中Φ可看作一個(gè)字典,它的每一列被稱作原子.如果k<<N且Φ滿足約束等距性質(zhì)(restricted isometry property,RIP),則存在精確重建算法使x從M=O(k lb(N))個(gè)測(cè)量值中實(shí)現(xiàn)無(wú)失真恢復(fù)[1].
兩種常用的CS重建方法分別基于l1范數(shù)最小和l0范數(shù)最小.前者利用線性規(guī)劃求解,有著很高的計(jì)算復(fù)雜度O(N3)[2];后者用迭代貪婪追蹤策略近似信號(hào)[3],其每次迭代包含最佳原子選擇、信號(hào)估計(jì)、殘差更新3步.匹配追蹤(matching pursuits,MP)[4]和正交匹配追蹤(orthogonal matching pursuit,OMP)[5]算法通過(guò)將殘差與最相關(guān)原子匹配來(lái)重建信號(hào),然而它們的解有時(shí)會(huì)被錯(cuò)誤原子干擾.為剪切錯(cuò)誤選擇的原子,在已知稀疏度下,可引入回溯策略.子空間追蹤(subspace pursuit,SP)[6]與壓縮感知匹配追蹤(compressed sensing matching pursuit,CoSaMP)[7]都是已知稀疏度k的回溯型重建算法,它們?cè)谠舆x擇環(huán)節(jié)分別保留k和2k個(gè)原子.另一種回溯算法是迭代硬閾值(iterative hard thresholding,IHT)[8]算法,它每次迭代保留近似信號(hào)幅值最大的k個(gè)元素xn=Hk(xn-1+ΦT(y-Φxn-1)).然而在實(shí)際應(yīng)用中稀疏度通常是未知的,為此回溯型自適應(yīng)正交匹配算法(backtracking-based adaptive orthogonal matching pursuit,BAOMP)[9]使用了一種基于閾值選擇的方法,它以固定的權(quán)重參數(shù)代替對(duì)稀疏度的估計(jì)來(lái)確定增減原子閾值,從而實(shí)現(xiàn)盲稀疏度下對(duì)原信號(hào)的重建[10].但BAOMP算法仍存在明顯缺陷:1)信號(hào)稀疏度k與閾值之間的關(guān)系沒有分析;2)經(jīng)過(guò)多次迭代后,搜索殘差會(huì)從相對(duì)小的值突然反彈到大的值,并且準(zhǔn)周期地反復(fù),降低了算法的收斂速度.
本文首先對(duì)殘差信號(hào)用類高斯分布建模,分析了BAOMP算法的閾值選擇方法與稀疏度已知算法使用稀疏度的一致性與差異點(diǎn),然后提出改進(jìn)的回溯型自適應(yīng)正交匹配算法(improved backtracking-based adaptiveorthogonal matching pursuit,IBAOMP),它采用80–20準(zhǔn)則判斷解信號(hào)的匹配程度,并在后續(xù)匹配步驟中引入可變步長(zhǎng)閾值,精細(xì)調(diào)整原子集容量,提高了選入原子的正確匹配率,從而避免了殘差信號(hào)的準(zhǔn)周期性失配.IBAOMP算法增加了2個(gè)參量:1)匹配準(zhǔn)確度閾值T,T決定新選原子集容量何時(shí)被調(diào)整;2)增量步長(zhǎng)S,S定義了閾值權(quán)重的調(diào)整步長(zhǎng).在第n次迭代時(shí),若上次迭代殘差模|rn-1|>T,則原子匹配策略與BAOMP算法相同;若|rn-1|≤T,則選入閾值權(quán)重μ1以步長(zhǎng)S增大到μ1+S,減少了新選原子集容量,使后續(xù)迭代步驟選擇的原子數(shù)目減少,從而能充分搜索出更接近原始信號(hào)的解.
因?yàn)樽罱K原子匹配集要包含k個(gè)原子,所以稀疏度k對(duì)精確重建十分關(guān)鍵,然而大多數(shù)情況下k為估計(jì)值.BAOMP算法通過(guò)設(shè)置閾值權(quán)重參數(shù)μ1和μ2來(lái)控制新選原子與淘汰原子數(shù)目,使重建信號(hào)近似k稀疏[9],選擇與淘汰原子策略如下:
1)新選原子集索引Cn
當(dāng)原子?i與殘差內(nèi)積滿足以下關(guān)系時(shí),原子索引i被選入集合Cn
式中,〈·〉表示內(nèi)積運(yùn)算,|·|表示取絕對(duì)值,i,j∈{1,2,···,N}.μ1度量了上次迭代殘差rn-1與原子?i內(nèi)積的最大值,較小的μ1將在每次搜索時(shí)選擇更多原子,因此加快了前期迭代的匹配速度,但會(huì)使迭代后期選入更多錯(cuò)誤原子.
2)淘汰原子集索引Γn
對(duì)新的候選原子索引集?n=Λn-1∪Cn,其中Λn-1為上次迭代近似信號(hào)支撐集,計(jì)算這些候選原子的近似系數(shù)
式中,“?”表示矩陣的偽逆.當(dāng)?n中對(duì)應(yīng)原子的近似系數(shù)小于淘汰閾值時(shí),該索引值加入淘汰原子索引集Γn,即
式中,i∈?,j∈Cn.μ2控制了淘汰原子數(shù)目.
下面本文通過(guò)對(duì)殘差信號(hào)建模分析BAOMP算法的閾值選擇方法與已知稀疏度k方法的關(guān)系,以及固定權(quán)重μ1導(dǎo)致的準(zhǔn)周期性失配現(xiàn)象.
1)稀疏度k與閾值的關(guān)系
BAOMP算法通過(guò)比較所有原子殘差內(nèi)積與μ1倍最大值來(lái)選擇新原子,從而用已選原子逼近上次迭代殘差.因?yàn)橛^測(cè)陣的元素通常是類高斯分布(如高斯隨機(jī)矩陣、部分傅里葉矩陣、Chirp感知矩陣),即使信號(hào)x是混合分布,殘差rn=y-也是類高斯分布.本文模擬超過(guò)300個(gè)殘差信號(hào)后確認(rèn)它們的分布如下:
式中,r為殘差r的分量,u為均值,δ為標(biāo)準(zhǔn)差,a為常數(shù).因此,搜索k個(gè)原子等價(jià)于通過(guò)搜尋一個(gè)正閾值t,使圖1中黑色區(qū)域與k成比例,從而逼近殘差,即
圖1 殘差信號(hào)幅值類高斯分布Figure 1 Gaussian-like distributed amplitudes of residuals
2)逼近稀疏度
對(duì)于類高斯分布,k可用μ1近似
式中,rmax是觀測(cè)信號(hào)或殘差的最大絕對(duì)值.k應(yīng)該由μ1和理論上可收縮的rmax近似得到,但實(shí)際上μ1為固定值,且rmax可能因錯(cuò)誤匹配而增大,于是對(duì)k的估計(jì)可能存在較大誤差.
3)準(zhǔn)周期性失配
盡管對(duì)k的估計(jì)可轉(zhuǎn)換成一個(gè)閾值尋找問(wèn)題,k依舊是未知的.當(dāng)殘差模值較小時(shí),固定的μ1會(huì)導(dǎo)致在如式(1)所示的全局搜索時(shí)選入錯(cuò)誤原子,并導(dǎo)致殘差模值劇烈反彈.圖2為圖像“Lena”在第34次和第35次迭代后的殘差,第34次迭代殘差的幅值范圍為[–1,+1],但一次錯(cuò)誤匹配迭代后,第35次迭代的殘差幅值范圍上升至[–71,+75].這兩次迭代的殘差分布如圖3所示.另一個(gè)高斯稀疏信號(hào)殘差模值|rn|準(zhǔn)周期性失配的例子如圖4所示,其中稀疏度k=12,觀測(cè)數(shù)目M=45,信號(hào)長(zhǎng)度N=256,第6次迭代后殘差2范數(shù)|rn|下降到3.6,但由于錯(cuò)誤匹配而導(dǎo)致第7次迭代后顯著上升至20.4.在極端情況下,當(dāng)達(dá)到最大迭代數(shù)時(shí),仍不能精確重建.
圖2 Lena圖像用BAOMP算法錯(cuò)誤匹配后殘差反彈Figure 2 Rebounded residuals of Lena after a mismatch by BAOMP
圖3 Lena圖像用BAOMP算法錯(cuò)誤匹配前后殘差幅值分布直方圖Figure 3 Histograms of amplitudes of residuals of Lena after a mismatch by BAOMP
為避免上述稀疏度k未知的BAOMP算法中固定的μ1引起的準(zhǔn)周期性失配現(xiàn)象,本文調(diào)整搜索方案去逼近稀疏度k.改進(jìn)的回溯型自適應(yīng)正交匹配算法(IBAOMP)定義了兩個(gè)參量:匹配準(zhǔn)確度閾值T和增量步長(zhǎng)S.與BAOMP中固定的μ1不同,T和S使μ1更加精準(zhǔn)地匹配殘差,當(dāng)殘差模值|rn|變小時(shí),選入的原子數(shù)目也隨之變少.具體調(diào)整如下:
圖4 高斯稀疏信號(hào)用BAOMP算法恢復(fù)時(shí)殘差范數(shù)隨迭代次數(shù)變化Figure 4 Residual norms versus the numbers of iterations of a Gaussian sparse signal by BAOMP
1)用可變?chǔ)?精確匹配
引入匹配準(zhǔn)確度閾值T,通過(guò)比較殘差模值|rn|確定μ1調(diào)整時(shí)刻.
當(dāng)|rn|>T時(shí),μ1保持不變來(lái)加速匹配;當(dāng)|rn|≤T時(shí),大多數(shù)正確原子已經(jīng)搜尋到,增加μ1來(lái)避免錯(cuò)誤投影與準(zhǔn)周期性失配.T的合理取值為
根據(jù)RIP原則,當(dāng)約束等距常數(shù)δ?1時(shí),||y||2≈||x||2,觀測(cè)信號(hào)y的模值可以代表x的能量.又根據(jù)80–20準(zhǔn)則,在大多數(shù)情況下,稀疏空間中不到20%的大幅值元素?fù)碛行盘?hào)80%的能量,因此式(8)是可行的.
2)增量步長(zhǎng)S控制
μ1以步長(zhǎng)S增加,如式(9)所示,S越大候選原子集容量越小
圖5顯示了用IBAOMP算法重建高斯稀疏信號(hào)時(shí)殘差模值|rn|隨迭代次數(shù)的變化,其中T=15,S=0.01,其他參數(shù)與圖4相同.作為對(duì)比,前3次類三角周期反彈并未被抑制.圖5中第18次迭代后殘差模值控制在|rn|<10-4,而圖4中,BAOMP算法中的殘差模在第18次迭代后又反彈了4次.
IBAOMP算法的實(shí)現(xiàn)步驟如下:
步驟1 迭代前期新選原子集索引.當(dāng)|rn|>T時(shí),與BAOMP算法相同產(chǎn)生新選原子集索引Cn,如式(1)所示.
圖5 高斯稀疏信號(hào)用IBAOMP算法恢復(fù)時(shí)殘差范數(shù)隨迭代次數(shù)變化Figure 5 Residual norms versus the numbers of iterations of a Gaussian sparse signal by IBAOMP
步驟2 迭代后期校正與步長(zhǎng)控制.當(dāng)|rn|≤T時(shí),校正μ1值至μ1+S,如式(9)所示.
步驟3 信號(hào)近似.產(chǎn)生新的候選原子索引集?n=Λn-1∪Cn,由等式(2)計(jì)算近似系數(shù)xΩn
步驟4 原子淘汰.由式(3)產(chǎn)生淘汰原子集索引Γn.
步驟5 更新支撐集索引Λn.更新近似系數(shù)xnΛ與殘差rn如下:
在IBAOMP算法的一次迭代中,計(jì)算N個(gè)內(nèi)積與鎖定Cn分別需要O(M N)和O(N)個(gè)標(biāo)準(zhǔn)乘,最小二乘計(jì)算和分別需要O(|?|N)和O(|Λ|N)個(gè)標(biāo)準(zhǔn)乘(其中|Λ|≤|?|<N).因此,與BAOMP算法相同,IBAOMP算法的一次迭代需要的計(jì)算量為O(M N)個(gè)標(biāo)準(zhǔn)乘.
為驗(yàn)證IBAOMP算法性能,進(jìn)行以下2組實(shí)驗(yàn):1)匹配準(zhǔn)確度閾值T和增量步長(zhǎng)S性能測(cè)試;2)比較IBAOMP算法與BAOMP、OMP、CoSaMP、St OMP、ROMP算法在不同稀疏度k和不同測(cè)量數(shù)M下精確重建概率.k稀疏信號(hào)為非零元素服從標(biāo)準(zhǔn)高斯分布的稀疏信號(hào),k=12,長(zhǎng)度N=256,觀測(cè)矩陣為歸一化高斯矩陣.可壓縮信號(hào)為3幅256×256圖像:頻率主要在低頻的“Peppers”圖像、低頻到高頻都有的“Lena”圖像、紋理豐富的高頻“Baboon”圖像,這3幅圖像的稀疏度依次由低到高,對(duì)于可壓縮信號(hào)采用部分傅里葉矩陣作為觀測(cè)矩陣.所有實(shí)驗(yàn)重復(fù)500次,在某一次重建實(shí)驗(yàn)中,對(duì)稀疏信號(hào)有‖x-‖2<10-4,對(duì)可壓縮信號(hào)有‖x-‖2<2時(shí),該次重建可看作精確重建.定義精確重建次數(shù)Nr與總次數(shù)N的比率為精確重建概率[11]
首先觀察當(dāng)S值固定,T均勻增大時(shí),匹配準(zhǔn)確度閾值T對(duì)信號(hào)精確重建性能的影響,T=0時(shí)的IBAOMP算法等同于BAOMP算法.然后選擇固定的T值,步長(zhǎng)S均勻增大時(shí),觀察S對(duì)信號(hào)精確重建性能的影響.實(shí)驗(yàn)首先用稀疏信號(hào)得到明確結(jié)果,然后用可壓縮信號(hào)測(cè)試實(shí)際應(yīng)用性能.
3.1.1 T和對(duì)S高斯稀疏信號(hào)重建性能的影響
在BAOMP算法中,μ1取0.4時(shí)精確重建效果最佳[9],參考這一結(jié)論將實(shí)驗(yàn)參數(shù)設(shè)置如下:S=0.01、M=45、μ1=0.4、μ2=0.6,T以步長(zhǎng)0.5由0增加到20,精確重建概率p和平均迭代結(jié)束時(shí)μ1的終值如圖6所示.
圖6 S=0.01時(shí)匹配準(zhǔn)確度閾值T對(duì)精確重建概率的影響Figure 6 Matching accuracy threshold T versus probability of exact reconstruction with S=0.01
從圖6中可以看出:隨著T的增加,μ1的終值增加,對(duì)應(yīng)的新增匹配原子集容量減少.T由0增加到20時(shí),IBAOMP算法的平均精確重建概率為82.4%,比BAOMP算法的69%高13.4%.當(dāng)T取19.5時(shí),IBAOMP算法的精確重建概率最高,達(dá)到89.6%,比BAOMP算法提高20.6%,說(shuō)明對(duì)于該高斯信號(hào)和步長(zhǎng)S的設(shè)定,T?取值使精確重建概率達(dá)到最高的19.5時(shí),μ1在合適的時(shí)刻被調(diào)整,最終能將精確重建概率提高到89.6%.
固定T=12、μ1=0.4,S以步長(zhǎng)0.001由0增加到0.02,實(shí)驗(yàn)結(jié)果如圖7所示.圖7表明S取較小值時(shí),后段迭代對(duì)μ1的調(diào)整較小,選取原子數(shù)目隨之平穩(wěn)減少,更適合該信號(hào)結(jié)構(gòu)特征,進(jìn)而有更高的精確重建概率.當(dāng)S=0.003時(shí),最終μ1=0.47已經(jīng)自適應(yīng)該信號(hào)結(jié)構(gòu),使得最終精確重建概率比BAOMP算法高出10.2%.
圖7 T=12時(shí)增量步長(zhǎng)S對(duì)精確重建概率的影響Figure 7 Incremental step S versus probability of exact reconstruction with T=12
3.1.2 T對(duì)可壓縮信號(hào)的影響
本文驗(yàn)證了重建圖像“Peppers”、“Lena”和“Baboon”時(shí)新增原子數(shù)目和淘汰原子數(shù)目隨迭代次數(shù)的變化.圖8為“Lena”圖像結(jié)果,其余兩幅圖結(jié)果類似,在此未給出.
由圖8(a)可看出,用BAOMP算法重建圖像,當(dāng)?shù)螖?shù)為35、66、97時(shí)新增原子數(shù)目發(fā)生劇烈變化,而此時(shí)對(duì)應(yīng)的淘汰原子數(shù)目(見圖8(b))也發(fā)生劇烈變化,表明這些迭代發(fā)生了錯(cuò)誤匹配.
圖9給出在參數(shù)T作用下對(duì)不同稀疏度信號(hào)重建時(shí)殘差模值隨迭代次數(shù)變化情況.圖9(a)和9(b)中“Peppers”與“Lena”圖像80%的能量由少量系數(shù)控制,IBAOMP算法明顯抑制了殘差模值的反彈.而對(duì)于圖9(c)中稀疏度較大的細(xì)節(jié)信號(hào)“Baboon”圖像,在選入錯(cuò)誤原子使得殘差模值反彈前,需要更多的迭代來(lái)達(dá)到初始匹配.雖然需要較長(zhǎng)的匹配過(guò)程,式(8)依舊有效,圖9中曲線拐點(diǎn)都出現(xiàn)在位置T=0.2norm(y).
參與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工具箱實(shí)現(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 不同測(cè)量數(shù)M下IBAOM P算法與其他重構(gòu)算法的精確重建性能比較
對(duì)k=12高斯稀疏信號(hào)用IBAOMP算法及其他典型貪婪算法重建結(jié)果如圖10所示.當(dāng)稀疏度保持不變時(shí),隨著測(cè)量數(shù)的增加,所有貪婪算法精確重建概率逐步上升.當(dāng)測(cè)量數(shù)M<75時(shí),由于測(cè)量數(shù)M過(guò)少,各類算法都不能以概率1完全重建原始信號(hào),但I(xiàn)BAOMP算法比BAOMP算法的精確重建概率高出26%;當(dāng)測(cè)量數(shù)M>75時(shí),這兩種算法都以概率1精確重建原信號(hào).CoSaMP和OMP算法在測(cè)量數(shù)M>80時(shí)的精確重建概率達(dá)到1,而StOMP和ROMP算法在測(cè)量數(shù)M=120時(shí)的精確重建概率分別為70.6%和98.2%.
圖10 不同測(cè)量數(shù)下IBAOMP算法與其他重構(gòu)算法比較Figure 10 Comparing IBAOMP with other algorithms with different numbers of measurement
3.2.2 不同稀疏度k下IBAOM P算法與其他重構(gòu)算法的精確重建性能比較
對(duì)M=148高斯稀疏信號(hào)用IBAOMP算法及其他典型貪婪算法重建結(jié)果如圖11所示.隨著稀疏度k的增加,需要的測(cè)量數(shù)更多,當(dāng)測(cè)量數(shù)M=148相對(duì)于稀疏度k不夠時(shí),精確重建概率開始下降.當(dāng)稀疏度k<70時(shí),IBAOMP算法以概率1精確重建;當(dāng)稀疏度k>70時(shí),IBAOMP算法比BAOMP算法的精確重建概率下降得更慢;當(dāng)稀疏度k=80時(shí),IBAOMP算法比BAOMP算法的精確重建概率高出17%.CoSaMP和OMP算法分別在稀疏度k>50和k>35時(shí)的精確重建概率迅速下降到零.StOMP算法在本實(shí)驗(yàn)中性能最差,當(dāng)稀疏度k為5時(shí),精確重建概率僅為75%.
圖11 不同稀疏度下IBAOMP算法與其他重構(gòu)算法比較Figure 11 Comparing IBAOMP with other algorithms with different sparsity levels
3.2.3 不同稀疏度圖像算法對(duì)比
分別用IBAOMP算法與BAOMP算法對(duì)稀疏度不同的圖像“Peppers”、“Lena”、“Baboon”進(jìn)行重建.為快速達(dá)到初步匹配,對(duì)于“Baboon”,取μ1=0.2;對(duì)于“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
本文提出的IBAOMP算法適用于壓縮感知中對(duì)稀疏度未知的信號(hào)重建.匹配準(zhǔn)確度閾值T觸發(fā)對(duì)μ1的調(diào)整,避免了準(zhǔn)周期性失配;增量步長(zhǎng)S調(diào)整了匹配原子集容量.實(shí)驗(yàn)表明IBAOMP算法解決了固定μ1可能帶來(lái)的原子錯(cuò)誤匹配及殘差模值增加的問(wèn)題,從而在計(jì)算復(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.