王 晨， 高美鳳

(江南大學(xué)輕工過程先進(jìn)控制教育部重點(diǎn)實(shí)驗(yàn)室，江蘇無錫214122)

信號(hào)通過抽樣由模擬信號(hào)轉(zhuǎn)變?yōu)閿?shù)字信號(hào)。為了預(yù)防信號(hào)畸變，傳統(tǒng)的抽樣以奈奎斯特采用定理為理論基礎(chǔ)。該定理指出:為了保證能夠完全恢復(fù)原信號(hào)，采樣速率必須大于等于信號(hào)帶寬的兩倍，否則信號(hào)頻譜會(huì)出現(xiàn)混疊，不能準(zhǔn)確地恢復(fù)原信號(hào)。隨著現(xiàn)代信息量日益增長(zhǎng)的需求，該定理則表現(xiàn)出越來越明顯的限制。一方面，隨著信號(hào)帶寬的增加，要求采樣速率和處理速度也要相應(yīng)提高，這樣會(huì)提高產(chǎn)品實(shí)現(xiàn)難度和成本;另一方面，在傳統(tǒng)的圖像或信號(hào)傳輸過程中，有效的數(shù)據(jù)往往只有一小部分，大部分采集數(shù)據(jù)會(huì)被舍棄，這樣造成了資源浪費(fèi)而且效率低下。

近些年來，一種新的采樣理論——壓縮感知［1］被廣泛關(guān)注。它利用其他變換空間描述信號(hào)，使得在保證信息不損失的情況下，用遠(yuǎn)低于采樣定理要求的速率在采樣信號(hào)的同時(shí)，又可完全恢復(fù)信號(hào)，即將對(duì)信號(hào)的采樣轉(zhuǎn)變?yōu)閷?duì)信息的采樣。壓縮感知突破了奈奎斯特采樣定理的限制，如果信號(hào)在某變換域中是稀疏的，并且其稀疏變換與觀測(cè)系統(tǒng)不相關(guān)，那么對(duì)信號(hào)無失真的采樣可通過優(yōu)化求解一個(gè)稀疏變換與觀測(cè)矩陣構(gòu)成的欠定方程來完成，因而可以從低分辨率觀測(cè)中恢復(fù)信號(hào)，大大提高了信息處理速度。隨著壓縮傳感的深入研究和廣泛應(yīng)用，信息處理領(lǐng)域正不斷得到較快發(fā)展［2］。

1 壓縮感知

壓縮感知理論指出，當(dāng)信號(hào)在某組基或者字典下是可壓縮或者稀疏的，那么就可以用一個(gè)與該基或字典不相關(guān)的觀測(cè)矩陣將信號(hào)線性投影為低維觀測(cè)值，其保留了重建信號(hào)所需要的重要信息，通過進(jìn)一步求解稀疏最優(yōu)化問題就能從線性觀測(cè)中精確地重構(gòu)原始信號(hào)。

壓縮感知理論框架如圖1所示:第1步為信號(hào)稀疏表示問題;第2步為信號(hào)的線性投影問題，設(shè)計(jì)一個(gè)與正交變換基Φ不相關(guān)的M×N維的觀測(cè)矩陣Ψ，保證原始信號(hào)從維降維到M維時(shí)重要信息不丟失;第3步為信號(hào)重構(gòu)問題，設(shè)計(jì)重構(gòu)算法保證從M維線性觀測(cè)向量中能夠重建原始信號(hào)。

圖1 CS理論框架Fig.1 Graph of CS theoretic framework

壓縮感知是利用信號(hào)具有稀疏性表示的先驗(yàn)知識(shí)，由觀測(cè)獲得少量的觀測(cè)值來進(jìn)行信號(hào)重構(gòu)?？紤]長(zhǎng)度為N的信號(hào)，將其看作RN中N×1的列向量，RN中的任意信號(hào)都可由一組N×1的矢量基或，那么信號(hào)X可以表示為

其中S=(s1，…，sN)T為信號(hào)在Φ域中的表示，若S的非零分量的個(gè)數(shù)為K(M?N)，稱S為K-稀疏的。衡量信號(hào)稀疏度常用的指標(biāo)有:(1)0范數(shù)，即信號(hào)中非零元素的個(gè)數(shù);(2)稀疏因子，即信號(hào)中非零元素的個(gè)數(shù)與元素總數(shù)的比值;(3)非線性逼近誤差，即用正交基Ψ表示信號(hào)時(shí)的稀疏程度或分解系數(shù)的能量集中程度。常用的稀疏基有DCT基、小波基、Gabor基及冗余字典等。

測(cè)量矩陣的設(shè)計(jì)是壓縮感知理論的一個(gè)重要前提和條件，既要減少壓縮測(cè)量次數(shù)又要確保信號(hào)的重建精度。Candès［3］等人指出，如果要精確重構(gòu)K階稀疏信號(hào) X，測(cè)量次數(shù) M需要滿足 M=O(Kln(N))，并且測(cè)量矩陣Ψ必須滿足約束等距性質(zhì)(Restrcited Isometry Property，RIP)。高斯隨機(jī)矩陣的優(yōu)點(diǎn)是它幾乎與任意稀疏信號(hào)都不相關(guān)，因此，當(dāng)測(cè)量矩陣是高斯隨機(jī)矩陣時(shí)，測(cè)量矩陣能以較大概率滿足RIP，因此可以通過選擇M×N的高斯測(cè)量矩陣得到所需的觀測(cè)矩陣

給定測(cè)量矩陣ΨM×N(M?N)，對(duì)信號(hào)X觀測(cè)過程如下:

式中，Y為經(jīng)過觀測(cè)后獲得的觀測(cè)值，包含了重構(gòu)原始信號(hào)的足夠信息;=為M×N矩陣，稱為傳感矩陣。由方程(2)看出，壓縮感知將信號(hào)X從N維降到M維觀測(cè)信號(hào)Y，由于觀測(cè)數(shù)量M遠(yuǎn)小于信號(hào)長(zhǎng)度N，所以方程(2)是個(gè)欠定線性方程組，存在無窮多解;同時(shí)由于S為K稀疏的，則可以通過系數(shù)分解算法求解式(2)的逆問題得到稀釋系數(shù)S，再通過式(1)重構(gòu)得到信號(hào)X。Tao［2］等指出當(dāng)傳感矩陣Θ滿足約束等距性質(zhì)RIP時(shí)，式(2)中的S可以通過L0的問題求解:

顯然L0為NP難的非凸優(yōu)化問題，Donoho，Chen和Saundrers指出求解更加簡(jiǎn)單的L1優(yōu)化問題會(huì)得到相同的解:

這一簡(jiǎn)單的改變使問題變成了凸優(yōu)化問題，可簡(jiǎn)化為線性規(guī)劃問題求解。目前有多種求解該問題的算法，典型重構(gòu)算法分為兩類:(1)貪婪追蹤算法:匹配追蹤(MP)［4］、正交匹配追蹤 (OMP)［5］、分段OMP(Stomp)［6］、正則化 OMP(Romp)［7］;(2)凸松弛法:基追蹤(BP)［8］、內(nèi)點(diǎn)法［9］、迭代閾值法［10］、共軛梯度投影法［11］等。

2 基于單層小波變換改進(jìn)的壓縮感知圖像重構(gòu)

2.1 單層小波變換稀疏分解

壓縮感知系統(tǒng)中，圖像的稀疏表示是實(shí)現(xiàn)壓縮感知的前提條件，目前普遍使用的的稀疏變換主要有正交變換稀疏表示和冗余字典稀疏分解［12］，其中正交變換構(gòu)造簡(jiǎn)單而且實(shí)現(xiàn)速度快，其中離散小波變換(DWT)進(jìn)行稀疏變換用的最為廣泛。在原有的壓縮感知圖像處理中，將原圖像首先進(jìn)行某種變換如DWT，然后構(gòu)造測(cè)量矩陣Ψ，利用該矩陣對(duì)小波變換后的稀疏系數(shù)進(jìn)行測(cè)量，得到M×N的測(cè)量系數(shù)，最后利用恢復(fù)算法通過得到的測(cè)量系數(shù)恢復(fù)原圖像。

小波分解將原圖像分解為低頻子帶和高頻子帶，高頻子帶可以認(rèn)為是稀疏的，而低頻子帶是原圖像在不同尺度下的逼近信號(hào)，不能認(rèn)為是稀疏的。如果將高低頻一起測(cè)量就破壞了低頻分量系數(shù)間的相關(guān)性，從而重構(gòu)效果不佳。在壓縮感知重構(gòu)算法中小波分解層數(shù)對(duì)圖像重構(gòu)有著重大的影響，分解層數(shù)越多重構(gòu)效果越好，大小為256×256，一般分解層數(shù)在4層以上。

文獻(xiàn)［13］提出了基于單層小波變換的壓縮感知改進(jìn)算法。根據(jù)小波分解圖像高低頻的特點(diǎn)，利用小波變換對(duì)原圖像進(jìn)行單層變換，對(duì)第一層高頻子帶進(jìn)行測(cè)量，對(duì)低頻逼近子帶直接觀測(cè)，這樣可以減少重構(gòu)圖像所需的數(shù)據(jù)量，加強(qiáng)重構(gòu)的效果。實(shí)現(xiàn)步驟如下:

1)對(duì)原圖像進(jìn)行一層小波分解得到{LL1，HL1，LH1，HH1}小波子帶系數(shù)。

2)構(gòu)造測(cè)量矩陣 Ψ 對(duì) LH1，HL1，HH1進(jìn)行測(cè)量，低頻LL1直接觀測(cè)。

3)利用OMP分別對(duì)高頻稀疏進(jìn)行重構(gòu)并與低頻子帶系數(shù)一起進(jìn)行小波反變換得到恢復(fù)圖像。

2.2 加權(quán)矩陣對(duì)小波分解稀疏增強(qiáng)

信號(hào)的稀疏程度決定了信號(hào)的重構(gòu)質(zhì)量和精度，因此對(duì)信號(hào)進(jìn)行有效的稀疏表示意義重大，它直接關(guān)系到壓縮感知的重構(gòu)精度。簡(jiǎn)單的稀疏正交變換雖然速度很快，但信號(hào)在此變換域下的稀疏性不高。因?yàn)楫?dāng)信號(hào)特征與稀疏變換的原子特征一致時(shí)能夠有效地得到信號(hào)精確的表示，但由于固定基不一定能靈活表示自然信號(hào)，因而在此稀疏域下不能夠足夠稀疏，則大大影響了壓縮感知的性能。對(duì)于文中所用的簡(jiǎn)單正交變換離散小波變換而言，它僅僅沿著圖像的水平和垂直兩個(gè)方向進(jìn)行變換。在圖像灰度漸變區(qū)域?qū)τ趫D像的稀疏表示效果很好，經(jīng)量化后高頻子帶產(chǎn)生大量的零系數(shù);但在非平滑區(qū)域，高頻子帶卻保留有較多的大幅度系數(shù)，滿足不了稀疏性的要求。因此，文中在此基礎(chǔ)上提出了一種增強(qiáng)其稀疏度的改進(jìn)方法。

二維離散小波對(duì)圖像進(jìn)行單層變換后保留低頻系數(shù)，高頻的稀疏系數(shù)S中除了有L個(gè)大系數(shù)，其他S～L個(gè)系數(shù)并不全為零，而有些是接近于零，不是絕對(duì)稀疏即不能足夠稀疏。為了使得高頻系數(shù)更加接近于絕對(duì)稀疏，文中提出引入加權(quán)矩陣的改進(jìn)方法。引入加權(quán)矩陣A，首先A是對(duì)角矩陣，且矩陣中系數(shù)

其中閾值系數(shù)λ為稀疏系數(shù)中選取的某一個(gè)系數(shù)的絕對(duì)值，加權(quán)系數(shù)設(shè)置的原則為使得稀疏系數(shù)中的大系數(shù)盡量保持不變而小系數(shù)盡量趨向于零［14］，即ai∝Si，加權(quán)系數(shù)矩陣左乘S得0，所以說當(dāng)p→+∞時(shí)，S中原本接近零的小系數(shù)右乘加權(quán)矩陣后變?yōu)榱?，因此，S'比S稀疏性更好，即引入加權(quán)矩陣后稀疏系數(shù)更加接近于絕對(duì)稀疏。

在對(duì)3個(gè)高頻子帶的觀測(cè)過程中引入加權(quán)矩陣增強(qiáng)其稀疏度，則壓縮感知的數(shù)學(xué)模型為

式中A為引入的對(duì)角加權(quán)矩陣，其中加權(quán)系數(shù){a1a2…an}在主對(duì)角線上而其他位置為0，相應(yīng)的壓縮感知重構(gòu)數(shù)學(xué)模型則變?yōu)?/p>

可將式(6)變換為

其中B=A－1，S'=AS。B矩陣也為對(duì)角矩陣B=diag{b1，b2，…，bn}且bi=，由此可知式(7)的最優(yōu)解即為式(8)的最優(yōu)解，即原圖像可以由重構(gòu)S'得到。

由此引入加權(quán)矩陣后可以使信號(hào)分解后的稀疏性進(jìn)一步加強(qiáng)，即將L個(gè)大系數(shù)保留，而剩余小系數(shù)變得更小，從而更接近理想的絕對(duì)稀疏，所以引入加權(quán)矩陣后3個(gè)高頻子帶的系數(shù)矩陣的稀疏性更好。

2.3 基于單層小波變換稀疏分解改進(jìn)與重構(gòu)

利用單層小波稀疏分解和加權(quán)系數(shù)矩陣的特點(diǎn)，通過單層小波變換獲得圖像高低頻子帶，保留低頻部分，對(duì)高頻系數(shù)進(jìn)行加權(quán)處理。

由上文提到重構(gòu)問題為NP難問題時(shí)，傳感矩陣必須滿足RIP特性，應(yīng)用加權(quán)系數(shù)矩陣后，傳感矩陣ΨΦ =Θ變?yōu)榱甩乏礎(chǔ)－1=Λ。由于加權(quán)系數(shù)矩陣僅僅使Θ中某些列向量的模的大小發(fā)生改變，而向量?jī)?nèi)部結(jié)構(gòu)并未發(fā)生改變，稀疏度也未發(fā)生改變。所以，如果原先稀疏矩陣與觀測(cè)矩陣不相關(guān)，則左乘加權(quán)系數(shù)矩陣后同樣滿足不相關(guān)性，因此Λ滿足RIP特性，可以重構(gòu)得到原圖像的恢復(fù)值。在重構(gòu)端由S'恢復(fù)到S時(shí)，則左乘矩陣B此時(shí)有

即由B和S'可以無損地恢復(fù)S。矩陣B:

其中μ為S'第L個(gè)最大值。

利用OMP算法恢復(fù)S'，最后將4個(gè)高低頻系數(shù)利用小波反變換恢復(fù)原圖像。

門限矩陣對(duì)OMP算法改進(jìn)后的算法流程如下:

輸入:傳感矩陣Λ =ΘA－1，觀測(cè)向量Y=ΛS'，稀疏度為K;

輸出:目標(biāo)信號(hào)S'的K稀疏逼近;

初始化:殘差r0=Y，索引集Γ0=?，t=1;

循環(huán)步驟1)～5):

5)判斷是否滿足t＞K，若滿足，則停止迭代;若不滿足，則執(zhí)行步驟1)。

由式(11)及上述流程可知:

從而有rt=由上述流程可以看出，門限矩陣加入前后目標(biāo)信號(hào)的稀疏逼近并沒有發(fā)生改變，又由于S'比S稀疏度更好，從而稀疏逼近與S'之間的誤差比與S之間的誤差小，所以重構(gòu)算法的精度得到了一定的改進(jìn)。

具體算法步驟如下:

1)對(duì)原始圖像進(jìn)行DWT，獲得高低頻4個(gè)小波子帶系數(shù)成分。

2)用加權(quán)矩陣A分別對(duì)LH1，HL1，HH1進(jìn)行處理，得到3個(gè)子帶的稀疏系數(shù)矩陣。

3)選擇合適的M值，構(gòu)造服從高斯分布的測(cè)量矩陣Ψ，分別對(duì)3個(gè)新的稀疏系數(shù)矩陣進(jìn)行測(cè)量，得到3個(gè)子帶的測(cè)量稀疏值矩陣，低頻LL1子帶系數(shù)保持不變。

5)將3個(gè)重構(gòu)高頻稀疏矩陣與LL1一起進(jìn)行小波反變換得到恢復(fù)的圖像。

3 實(shí)驗(yàn)及結(jié)果分析

實(shí)驗(yàn)基于PC機(jī)平臺(tái)(CPU主頻2.5 GHz，內(nèi)存2 GB)，并用Matlab7.0實(shí)現(xiàn)程序仿真。為了取得良好的實(shí)驗(yàn)結(jié)果，分別選擇了256×256 lena和256×256 Cameraman作為實(shí)驗(yàn)圖像，選取近似對(duì)稱小波Sym8對(duì)圖像進(jìn)行單層分解，測(cè)量矩陣為服從(0，1/N)的隨機(jī)高斯矩陣，利用OMP算法進(jìn)行恢復(fù)。根據(jù)單層小波分解高頻系數(shù)的特點(diǎn)，每個(gè)分量稀疏分解后的系數(shù)大系數(shù)L即為加權(quán)系數(shù)矩陣的閾值，則單層高頻水平分量的加權(quán)系數(shù)矩陣的閾值系數(shù)為第3 200個(gè)絕對(duì)值最大值，高頻垂直分量閾值系數(shù)為第800個(gè)絕對(duì)值最大值，高頻對(duì)角分量閾值系數(shù)為第190個(gè)絕對(duì)值最大值，選取參數(shù)p=20。

由于單層分解后高頻子帶均為128×128，所以測(cè)量矩陣行數(shù)滿足0＜M≤128，選取M=120進(jìn)行文中算法和單層小波變換壓縮感知算法對(duì)兩幅圖像進(jìn)行重構(gòu)，重構(gòu)結(jié)果如圖2，3所示。分別取M=40，50，…，120采用文中方法和僅僅采用單層小波變換的壓縮感知方法進(jìn)行仿真;選取M=100，110，…，250，采用傳統(tǒng)的算法進(jìn)行同樣仿真。

圖2 改進(jìn)算法對(duì)比Lena圖像Fig.2 Restoring algrithms contrast of lena image

圖像恢復(fù)質(zhì)量的評(píng)判標(biāo)準(zhǔn)為圖像的峰值信噪比，信噪比越高，圖像的恢復(fù)質(zhì)量越好。實(shí)驗(yàn)結(jié)果如圖4，5所示，其中橫軸為經(jīng)過測(cè)量后的系數(shù)點(diǎn)個(gè)數(shù)，縱軸為恢復(fù)圖像的峰值信噪比(PSNR)，PSNR為多次測(cè)量取平均值?？梢钥闯鑫闹刑岢龅母倪M(jìn)算法比原有算法有較明顯的提高。

4 結(jié)語

介紹了壓縮感知算法的原理和數(shù)學(xué)模型，在已有的單層小波分解壓縮感知算法的基礎(chǔ)上運(yùn)用加權(quán)矩陣對(duì)算法進(jìn)行了改進(jìn)，利用加權(quán)系數(shù)矩陣使信號(hào)的稀疏分解進(jìn)一步加強(qiáng)，接近更加理想的絕對(duì)稀疏;重構(gòu)時(shí)采用OMP算法，仿真效果驗(yàn)證了改進(jìn)方法的有效性。

［1］Donoho D.Compressed sensing［J］.IEEE Transactions on Information Theory，2006，52(4):1289-1306.

［2］LI De-lai，ZHANG Qiong，SHEN Hai-hong，et al.The application of compressive sensing based on wavelet in the reconstruction of ultrasonic medical image［C］//Second International Conference on Mechanic Automation and Control Engineering.Jinan，China:MACE，2011:5350-5353.

［3］Candes E，Romberg J，TAO T.Robust uncertainty principles:exact signal reconstruction from highly incomplete frequency information［J］.IEEE Transactions on Information Theory，2006，52(4):489-509.

［4］Mallat S G，ZHANG Zhi-feng.Matching pursuits with time-frequency dictionaries［J］.IEEE Transactions on Signal Processing，1993，41(12):3397-3415.

［5］Troppp J 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］Donoho D，Tsaig Y，Drori I.Sparse solution of underdetermined linear equations by stagewise orthogonal matching pursuit［J］.IEEE Transactions on Information Theory，2012，58(2):1094-1121.

［7］Needell，Deanna，Vershynin，et al.Signal recovery from incomplete and inaccurate measurement via regularized orthogonal matching pursuit［J］.IEEE Journal on Selected Topics in Signal Processing，2010，4(2):310-316.

［8］Hen S S B，Donoho D L，Saunders M A.Atomic decomposition by basis pursuit［J］.Siam Journal on Scientific Computing，1998，20(1):33-61.

［9］Malioutov D M，Cetin M，Willsky A S.Homotopy continuation for sparse signal representation［C］//2005 IEEE International Conference on Acoustics，Speech，and Signal Processing.Philadelphia:ICASSP，2005:733-736.

［10］Nowak R D，F(xiàn)igueiredo M A T.Sparse reconstruction by separable approximation［J］.IEEE Transactions on Signal Processing，2009，57(7):2479-249.

［11］Harmany Z，Thompson D，Willett R，et al.Gradient projection for linearly constrained convex optimization in sparse signal recovery［C］//2010 17th IEEE International Conference on Image Processing 2010.Hong Kong:ICIP，2010:3361-3364.

［12］Elad M.Sparse，Redundant Representations:From Theory to Applications in Signal and Image Processing［M］.New York:Springer，2010.

［13］岑翼剛，陳曉方，岑麗輝，等.基于單層小波變換的壓縮感知圖像處理［J］.通信學(xué)報(bào)，2010，31(8A):52-55.CEN Yi-gang，CHEN Xiao-fang，CEN Li-hui，et al.Compressed sensing based on the single layer wavelet transform for image processing［J］.Journal on Communications，2010，31(8A):52-55.(in Chinese)

［14］趙玉娟，鄭寶玉.壓縮感知中稀疏分解和重構(gòu)精度改進(jìn)的一種方法［J］.信號(hào)處理，2012，28(5):631-636.ZHAO Yu-juan，ZHENG Bao-yu.An improved method for sparse decomposition and resconstruction in compressed sensing［J］.Singal Processing，2012，28(5):631-636.(in Chinese)