劉繼忠，鄭恩濤，賀艷濤，付珊珊，趙鵬

(南昌大學(xué)機(jī)電工程學(xué)院，江西南昌330031)

0 引言

近年來，壓縮感知（Compressed Sensing）在觀測矩陣、稀疏基的選擇以及重構(gòu)算法等三個方面已經(jīng)做了很多的研究，其中以重構(gòu)算法的優(yōu)化為最多，它的應(yīng)用已經(jīng)滲透到雷達(dá)、信號處理、圖像傳輸?shù)雀鞔箢I(lǐng)域.不論何種應(yīng)用，凸優(yōu)化算法一直以其重構(gòu)精度高而為大家所認(rèn)可，凸優(yōu)化中BP算法最為經(jīng)典，然而美中不足的是它的計算復(fù)雜度太高了，計算量大[1]，使得它在圖像傳輸這一領(lǐng)域并不被大家所青睞.2007年，Lu Gan提出了分塊壓縮感知(Block Compressed Sensing，BCS)，該方法不僅能解決大尺度圖像實時傳輸?shù)膯栴}，使計算復(fù)雜度大大減少，而且在觀測時每一個圖像塊的測量矩陣遠(yuǎn)小于整幅圖像的測量矩陣，便于存儲[2].

針對以上情況，本文提出基于分塊圖像的BP壓縮感知重構(gòu)算法（由于本文算法較多，以下未經(jīng)說明默認(rèn)本文所提算法為該算法）.該算法既能夠承襲BP算法重構(gòu)精度高的優(yōu)點，又融合了圖像分塊所帶來的計算便利，傳輸壓力低的特點.此外，還把比較成熟的應(yīng)用于腦電、心電領(lǐng)域的過完備字典稀疏表示法[3]移植到二維圖像壓縮重構(gòu)中，并且與傳統(tǒng)的不分塊圖像的BP、COSAMP、OMP、算法進(jìn)行對比得出：本文所提算法不僅重構(gòu)精度高，而且運行時間快.

1 基本原理

1.1 壓縮感知理論模型

CS理論對于稀疏信號用一個與變換基不相關(guān)的觀測矩陣將高維信號投影到一個低維空間上，然后通過求解一個等效優(yōu)化問題從這些投影中重構(gòu)出原信號[4?6].該模型具體過程如下：首先將信號s投影到觀測矩陣Φ上得到觀測向量y，即

式中y是M×1的觀測向量，Φ是M×N(M?N)的觀測矩陣[7]，信號s是長度為N的離散實值信號或者在基ψ=[ψ1,ψ2,...,ψN]上是稀疏的信號，即

式中α是N×1的稀疏向量，且僅有k?N個非零系數(shù)，k是稀疏度；ψ是N×N的稀疏基矩陣.將式(2)代入式(1)得

式中A為傳感矩陣.

因觀測維數(shù)M遠(yuǎn)小于信號維數(shù)N，所以無法直接從y的M個觀測值中解出信號s[8].根據(jù)s的稀疏性，可以通過求解式(3)的逆問題解得α，然后將α代入式(2)，求得原始信號s[3,9?12].

1.2 基空間追蹤（BP）算法

基追蹤算法從完備（過完備）原子庫中尋求信號最稀疏的表示,也就是利用最小化l1范數(shù)方法把信號的稀疏表示問題等價為有約束的極值問題,進(jìn)而將其轉(zhuǎn)化成求解線性規(guī)劃問題.具體過程如下：對于給定的信號

求s最簡單的表示形式,即求解式(5)：

其中原子庫ψri是線性獨立的，但是不再滿足正交性，并且滿足信號自適應(yīng)的選取.通常情況下集合ri是未知的，它依賴于待處理信號s[13,14].

1.3 分塊壓縮感知

分塊CS算法通過將圖像進(jìn)行分割，得到大小均為B×B的不重疊的各塊圖像，然后選取相同的觀測矩陣分別對其進(jìn)行觀測.具體過程如下：假設(shè)si是原始輸入圖像S的所分的第i個小圖像子塊，那么就可以將它們與觀測集合的關(guān)系表示如下式：

其中ΦB是大小為MB×B2的觀測矩陣，那么對整個圖像來說，觀測矩陣Φ是一個塊對角矩陣，Φ可表示如下[15]：

1.4 基于分塊圖像的BP重構(gòu)算法

本文所提算法中，首先將原始圖像（N×N）均勻分塊，對每個子塊圖像（B×B）的數(shù)據(jù)轉(zhuǎn)化成一列矩陣后進(jìn)行離散余弦變換，然后采用相同的觀測矩陣ΦB單獨對每個圖像子塊所轉(zhuǎn)換后的列矩陣進(jìn)行觀測采樣并重構(gòu)，最后采用BP算法單獨循環(huán)重構(gòu)出每個小塊圖像所對應(yīng)的列向量，并把列向量的數(shù)據(jù)轉(zhuǎn)化成塊矩陣，然后按照分塊順序重新組合得到新生成的圖像.此方法不僅可以減少運算復(fù)雜度，而且不需要傳輸所有的觀測數(shù)據(jù)，減少了傳感器部分的存儲量，從而也加快了重構(gòu)速度.

具體算法流程如下：

Step 1：將像素為N×N圖像S均勻分塊處理，設(shè)置各分塊圖像si的大小均為B×B，其中i=1,2···N/B；

Step 2：把所分大小為B×B的塊矩陣si轉(zhuǎn)化為si=B2?1的列矩陣;

Step 3：選取合適的M值，生成M×(B2)維的高斯隨機(jī)測量矩陣，并對其每列矩陣單獨循環(huán)做歸一化處理得到觀測矩陣ΦB，對si進(jìn)行測量得到測量矩陣yi；

Step 4：選用離散余弦變換基B2?B2維的ψ，并求出傳感矩陣A=ΦB?ψ,;

Step 5：把Step3、4中得出的測量矩陣yi和傳感矩陣A=ΦB?ψ,作為BP算法的輸入，初始化稀疏度為測量矩陣yi的維度的z倍，初始化稀疏系數(shù)αinitial=A0?yi，然后利用基于線性規(guī)劃的凸優(yōu)化問題求解L1范數(shù)問題，從而獲得稀疏解α；

Step 6：對重構(gòu)后的頻域信號通過稀疏表示的逆變換，即稀疏解α通過公式（2）變換到時域得到si，把列矩陣轉(zhuǎn)化為塊矩陣并進(jìn)行保存；

Step 7：循環(huán)執(zhí)行Step 2ΛStep 6,直到對所有小塊圖像si壓縮重構(gòu)完畢，進(jìn)而重新按照順序組合成整幅圖像.

1.5 其他算法

基于分塊圖像的COSAMP算法只是在分塊圖像的BP重構(gòu)算法中的step 5中把BP算法更換為COSAMP算法，其他初始化參數(shù)和執(zhí)行步驟如出一轍.

基于過完備字典稀疏表示的BP算法過程如下：首先把二維灰度圖像s對應(yīng)的（N?N）矩陣轉(zhuǎn)化成8?((N?N)/8)扁距陣s1；然后用隨機(jī)高斯矩陣歸一化后得到測量矩陣PHI并對s1進(jìn)行測量得到y(tǒng)1，接著運用奇異值分解的方法得到稀疏字典D[24]，傳感矩陣A由公式A=PHI?D得到；最后把A和y1作為BP算法的輸入求得稀疏系數(shù)向量α1，再由公式s1=D?α逆向得出重構(gòu)圖像的扁距陣形式，進(jìn)而逆向轉(zhuǎn)化成（N?N）的圖像矩陣.基于過完備字典稀疏表示的OMP算法只是把基于過完備字典稀疏表示的BP算法中的BP算法換成OMP算法，其他的過程不變.

2.工作條件：(1)相關(guān)的論文與著作以及當(dāng)代畫家畫作的圖片資料。(2)購買一些高清繪畫圖冊，以便更好的了解繪畫技能。(3)咨詢相關(guān)專業(yè)的老師，以及通過自己導(dǎo)師的指點改進(jìn)自己的不足。

2 實驗結(jié)果

2.1 評價標(biāo)準(zhǔn)和實驗環(huán)境

為了驗證上述的分析，本文采用客觀評價標(biāo)準(zhǔn)PSNR(單位為dB)來評價各種算法重構(gòu)圖像的質(zhì)量，采用運行時間t(單位為秒)來評價各種算法重構(gòu)算法的復(fù)雜度.同一算法在同等條件下，運行時間不是固定的而是有波動的，因此本文在相同硬件配置條件下，對同一個算法運行50次,取50次運行時間的平均值作為一個客觀評價標(biāo)準(zhǔn)[16].

定義圖像的壓縮比cr，重構(gòu)圖像的峰值信噪比PSNR的表達(dá)式如下:

其中M為壓縮后矩陣的行數(shù)，N為壓縮前矩陣的行數(shù)，為每段信號，為重構(gòu)的每段信號，L為原始信號的總段數(shù)[17].本文實驗都在win8.1系統(tǒng)下的Matlab2015a環(huán)境下完成，硬件配置為Inter(R)Core(TM)i7-5500U CPU@240GHz 240GHz,內(nèi)存為8.00GB.系統(tǒng)類型：64位操作系統(tǒng).

2.2 仿真實驗結(jié)果與分析

采用離散余弦基為正交基，觀測矩陣為服從高斯分布的隨機(jī)觀測矩陣，并對觀測矩陣的每一列逐一進(jìn)行歸一化處理.分別用以下算法：基追蹤算法（BP）、壓縮采樣匹配追蹤算法（COSAMP）[6]、基于圖像分塊的基追蹤算法（BP）、基于圖像分塊的壓縮采樣匹配追蹤算法（COSAMP）、正交匹配追蹤算法（OMP）等對二維灰度圖像（8位）進(jìn)行壓縮重構(gòu)；另外，還有基于過完備字典的BP算法和基于過完備字典的OMP算法，初始化稀疏度K= floor(1/4N)，N為壓縮前矩陣的行數(shù).本文算法如未明確說明均以DCT基做稀疏基表示.首先對標(biāo)準(zhǔn)圖片lena（256×256）圖，以采樣率cr=0.7時為例，以上7種算法重構(gòu)效果如圖1所示.

圖1 7種算法重構(gòu)效果圖

分別在采樣率為0.3，0.4，0.5，0.6，0.7時對以上幾種算法進(jìn)行仿真測試，結(jié)果分析如圖2.

圖2 7種算法與采樣率（cr）的關(guān)系

由圖2可知：總體而言， 隨著采樣率的增加，PSNR越來越大，對應(yīng)的重構(gòu)精度越來越高；總體上在不同算法之間，分塊-BP算法>分塊-COSAMP算法>KSVD-BP算法>BP算法>KSVDOMP>OMP>COSAMP算法.此外，當(dāng)cr=0.3～0.4時，前半部分的分塊-COSAMP算法>分塊-BP算法，以過完備字典為稀疏基的BP壓縮感知重構(gòu)算法比以DCT為稀疏基BP壓縮感知算法重構(gòu)精度高，分塊BP算法不但繼承了BP算法全局搜索的優(yōu)點，也繼承了分塊所帶來的增加稀疏基個數(shù)的優(yōu)點，總體效果最好.

由圖3可知：總體運行時間隨著采樣率的增加而增加，各算法之間運算時間大小關(guān)系，OMP分塊-BP算法，但其重構(gòu)精度優(yōu)勢并不明顯.結(jié)合運行時間，綜合圖2、3結(jié)果得出：本文所提算法不僅重構(gòu)精度比較高，而且運行時間比較短.

圖3 7種算法運行時間與采樣率（cr）的關(guān)系

此外，本文算法也對peppers.bmp(256×256)、camera.bmp(256×256)、baboon.bmp(256×256)、Plane.bmp(512×512)、barbara.jpg(512×512)、SAR1.jpg(512×512)等圖片進(jìn)行測試，設(shè)定cr=0.7.結(jié)果如圖5所示.

圖4 重構(gòu)效果圖

圖4中A～F左側(cè)均為原圖，右側(cè)均為重構(gòu)后的圖，由此可以看出該算法適用于（256×256）、（512×512）等大小不同、質(zhì)地紋理不同的圖片.

實驗a:以選擇最優(yōu)分塊大小為目的，設(shè)置稀疏度為K= floor(1/4N)（其中N為壓縮采樣后信號的維度），分別對lena（64×64）、lena（128×128）、lena（256×256）、peppers.bmp(256×256)、Plane.bmp(512×512))這些圖片，在采樣率為0.3，0.4，0.5，0.6，0.7時，分別設(shè)置分塊大小4×4、8×8、16×16、32×32、64×64進(jìn)行實驗.由于當(dāng)分塊大小為128×128、256×256時，從原理的角度講所分塊太大，使得觀測矩陣ΦB維數(shù)增多，計算量大大增加，從而失去了分塊的意義，故沒有給出具體運算結(jié)果.實驗結(jié)果如表1.

表1 在不同采樣率下不同分塊大小的重構(gòu)結(jié)果

由表1中A、B、C、D對應(yīng)實驗結(jié)果可知：對于lena（64×64）、lena（128×128）、lena（256×256）、Plane.bmp(512×512)等大小、質(zhì)地紋理不同的圖片，可觀察數(shù)據(jù)的大致規(guī)律分為cr≤0.4和cr≥0.5兩種情況來分析，在cr≤0.4時，PSNR隨著分塊大小從4*4到64*64時逐漸增大，其運算時間總體上越來越長，且在分塊為8*8以后迅速增長（其中只有l(wèi)ena（256×256）、Plane.bmp(512×512))這兩種圖片在cr=0.3時，t4?4略>t8?8，可忽略不計）.結(jié)合人類視覺需求≥25dB以上數(shù)據(jù)都滿足，考慮到圖像信號傳輸?shù)膶崟r性，運算時間越短越好，故選擇分塊為4*4時為最佳.然而4*4和8*8時運算時間接近，但8*8時對應(yīng)的PSNR要高出2～5dB不等，故在cr≤0.4分塊大小8*8為最優(yōu)；在cr≥0.5時，PSNR隨著分塊大小從4*4到32*32時逐漸增大，從32*32～64*64時再變小，t則相應(yīng)的在4*4～8*8時先變小(變化不大)，再由8*8～64*64時變長，與cr≤0.4時同理可得分塊大小8*8為最優(yōu).綜上所述：對于不同大小、質(zhì)地、紋理的圖片，在進(jìn)行圖像實時傳輸時，分塊大小為8*8時最優(yōu).

表2 不同稀疏度下對應(yīng)不同圖片的重構(gòu)結(jié)果

實驗b:以選擇最優(yōu)稀疏度為目的，按照a組的結(jié)果不妨設(shè)置分塊大小為8×8，分別采樣率為0.3，0.4，0.5，0.6，0.7時的情況下，分別對lena（64×64）、lena（128×128）、lena（256×256）、peppers.bmp(256×256)、camera.bmp(256×256)、baboon.bmp(256×256)、Plane.bmp(512×512)、barbara.jpg(512×512)、SAR1.jpg(512×512)等圖片，分別設(shè)置稀疏度為K= floor(x?N)（x=0.1,0.15,0.2···0.5；其中N為壓縮采樣后信號的維度）進(jìn)行實驗.（當(dāng)cr=0.4，cr=0.5，cr=0.6，cr=0.7，得出的結(jié)果過程和cr=0.3時相同，這里只列出了cr=0.3時的結(jié)果），實驗結(jié)果如表3.

由表2選擇出最高的PSNR和較小t對應(yīng)的稀疏度系數(shù)x可得表3.

表3 不同圖片在不同采樣率下的最優(yōu)稀疏度

表3中把表2中省略的cr=0.4，cr=0.5，cr=0.6，cr=0.7的最優(yōu)解也顯示了出來.由表3可知：對于不同大小、質(zhì)地、紋理的圖片，其最優(yōu)的稀疏度是不同的.總體上，最優(yōu)稀疏度隨著采樣率的增加而減小，設(shè)置稀疏度為K= floor(0.2～0.4?N)（其中N為壓縮采樣后信號的維度）為最優(yōu).

3 結(jié)論

本文針對BP壓縮感知重構(gòu)算法在圖像傳輸領(lǐng)域重構(gòu)精度不夠高、運行時間一般的特點提出基于分塊圖像的BP壓縮感知重構(gòu)算法.該算法不但繼承了圖像分塊降低了計算復(fù)雜度的特點，同時也承接了基追蹤（BP）算法全局搜索能力強(qiáng)、重構(gòu)相對精度高的優(yōu)點.并提出分塊大小為8?8、稀疏度為K= floor((0.2～0.4)?N)為最優(yōu).由基于過完備稀疏基的表示方法，可以使得稀疏表示更為精確，硬件也相對容易實現(xiàn)這一優(yōu)點，作者也曾試圖把本文所提算法DCT稀疏表示改為過完備字典稀疏（KSVD）表示，然而KSVD稀疏表示時運算復(fù)雜度大為提升，從理論上KSVD稀疏表示以接受扁矩陣為前提，而本文所分塊皆為正方矩陣，所以是否可以擴(kuò)展正矩陣為扁矩陣這一問題值得討論.另外也可考慮用稀疏二進(jìn)制隨機(jī)觀測矩陣代替本文的隨機(jī)高斯觀測矩陣.