李明飛， 袁梓豪， 趙琳琳， 孫曉潔

（1.中國航天科技集團(tuán)有限公司量子工程研究中心，北京100094；2.北京航天控制儀器研究所，北京100039）

0 引言

壓縮感知（壓縮采樣）（Compressive Sampling，CS）理論由 Candès等于 2006 年正式提出［1-2］， 美國Rice大學(xué)在2008年利用單像素成像實驗驗證了其有效性。CS理論是對稀疏或可壓縮的信號進(jìn)行少量非適應(yīng)性的線性測量，可以近似100%地恢復(fù)出原始信號，可突破Shannon采樣定理的限制，同時還具有超靈敏探測等優(yōu)勢［3］。目前，用于單像素成像壓縮感知算法的測量次數(shù)僅為Nyquist采樣極限的1%左右，即可重建出目標(biāo)圖像［4］。壓縮感知算法是非線性迭代算法，且重建效果與所選稀疏基密切相關(guān)，故成像計算時間和穩(wěn)定性同時受到極大挑戰(zhàn)。而強(qiáng)度關(guān)聯(lián)算法具有硬件要求低、成像速度快的特點，但其需要大于Nyquist采樣極限的測量次數(shù)，故采樣時間較長。

基于Walsh變換（Walsh Transform，WT）的方法提供了快速成像方案，該方案存在快速變換算法且硬件易于實現(xiàn)。采用二值矩陣作為調(diào)制矩陣，易于在高速數(shù)字微鏡器件（Digital Micro-mirror Device，DMD）上實現(xiàn)。調(diào)制矩陣生成速度快且無需存儲，Walsh矩陣作為測量矩陣具有正交特性，可完美重建圖像?；陔x散余弦變換（Discrete Cosine Transform，DCT）的方法也提供了實現(xiàn)快速單像素成像的方案，除了具有快速變換的特點和正交性的優(yōu)點外，DCT基是連續(xù)的，更適用于圖像壓縮，這是也聯(lián)合圖像專家組（Joint Photographic Experts Group，JPEG）采用DCT作為圖像壓縮方案的原因。當(dāng)然，DCT相比于WT的重要區(qū)別還在于DCT的連續(xù)性，在用于單像素成像時不利于實現(xiàn)快速光學(xué)調(diào)制，如DCT無法直接使用DMD實現(xiàn)20kHz的二值調(diào)制。文獻(xiàn)［5］報道的抖動二值化方法使得DCT這類連續(xù)測量基實現(xiàn)高速調(diào)制測量成為現(xiàn)實，研究發(fā)現(xiàn)二值化后的DCT基與Walsh基非常相似，然而在單像素成像中二者誰更占優(yōu)仍未有人研究。因此，在單像素成像中Walsh基和二值化后的DCT基各自優(yōu)化排序后，哪一種成像信噪比更高和重建速度更快仍是值得研究的課題。

針對上述問題，本文開展了對基于Walsh基和二值化后的DCT基的單像素成像數(shù)值仿真實驗對比研究。對DCT矩陣進(jìn)行了抖動二值化，分別對比了在不同采樣率條件下的基于兩種變換的單像素成像信噪比和成像時間，對比研究了相同采樣率下測量基三種不同排序時的單像素成像信噪比。

1 WT和DCT與測量基排序

1.1 WT與測量基

Walsh矩陣和Hadamard矩陣具有相同的矩陣元素，區(qū)別僅在于構(gòu)造矩陣的行和列排列不同。Hadamard矩陣起源較早，應(yīng)用極其廣泛［6-7］。Hadamard矩陣（簡稱H矩陣）具有如下特征：1）矩陣產(chǎn)生速度快；2）矩陣元僅有1元素。矩陣任意兩行（或列）正交， 有

若K為矩陣階數(shù)，則K需等于H矩陣的行數(shù)k或列數(shù)l，H（k，l）為Hadamard矩陣元。H矩陣正交歸一，有HHT＝KI，I為單位矩陣，且滿足H＝H-1。對H矩陣按逆序Gray碼排序后可得到Walsh矩陣（簡稱W矩陣），設(shè)圖像矩陣元為T（m，n），其二維WT的表達(dá)式如下

式（2）中，b（k，l，p，q）為對 2 取模運算， 其表達(dá)式為

式（3）中，g0（p）＝pn-1，g1（p）＝pn-1＋pn-2，g2（p）＝pn-2＋pn-3， …，gn-1（p）＝p1＋p0。 下標(biāo)i＝0， 1，2， …，n-1分別對應(yīng)ki、li、pi、qi二進(jìn)制表示時相應(yīng)的位。由于Walsh矩陣元為-1和1兩種元素，需采用差分探測的方式，即先將矩陣元-1變?yōu)?，1變?yōu)?，構(gòu)成D-1矩陣；將矩陣元-1變?yōu)?，1不變，產(chǎn)生D＋1矩陣，通過利用D＋1-D-1的方式來實現(xiàn)Walsh矩陣完整的測量。理論上全測量時，若圖像大小為N×N，則全采樣的測量次數(shù)為2×N×N，當(dāng)采用互補(bǔ)調(diào)制［8］的方式時，實際測量次數(shù)需要N×N次。

1.2 DCT與測量基

二維DCT與離散Fourier變換（Discrete Fourier Transform，DFT）的數(shù)學(xué)表達(dá)式非常接近，實際上DCT就是DFT只截取余弦變換部分生成。二維DCT的數(shù)學(xué)表達(dá)式如下

其中，M和N為物體的行數(shù)和列數(shù)，0≤m≤M-1， 0≤n≤N-1。

DCT矩陣和Walsh矩陣的測量基相似，但重要區(qū)別在于DCT是具有灰度的測量基。若采用DMD為調(diào)制器，按8bit灰度投影計算，DCT基測量速度是Walsh基測量速度的1／8，相應(yīng)的單像素成像速度也將是Walsh基方案的1／8，全采樣時二者信噪比相同。

經(jīng)上述分析，從實際參考價值上看，編碼同為灰度編碼或同為二值編碼的比較才能具有實際意義，直接使用DCT的灰度測量基與Walsh基比較則意義不大。

1.3 DCT基二值化-抖動算法

張子 邦 等［5，8］引 入 了 抖 動 算 法， 將 快 速Fourier變換（Fast Fourier Transform， FFT）有灰度的正弦圖案轉(zhuǎn)換成二值正弦圖案，提高了基于FFT測量基的調(diào)制速度，使得單像素成像速度大幅提升。本文引入了抖動算法，將DCT基進(jìn)行二值化，記為DCTb。不同于文獻(xiàn)［5］和文獻(xiàn)［8］之處在于，本文的DCT抖動二值化不對測量基進(jìn)行上采樣或插值，即抖動前后的矩陣均與原始矩陣大小相等，這樣做的優(yōu)勢在于在實際應(yīng)用中更加接近真實情況。

前文提到的抖動算法的核心思想是基于誤差擴(kuò)散原理，首先對當(dāng)前像素值進(jìn)行甄別，采用閾值法將當(dāng)前像素值二值化后輸出，然后將輸入和輸出的像素值之差按一定的權(quán)重傳播到若干相鄰未處理的區(qū)域。誤差擴(kuò)散原理可表達(dá)為

式（8）中，G（m，n）為輸入的 DCT 圖案，G*（m，n）為像素點（m，n）鄰域像素加入的量化誤差擴(kuò)散值與輸入圖案灰度之和，量化誤差error（m，n）擴(kuò)散到鄰域像素的二維權(quán)重函數(shù)（也稱為核函數(shù)）為w（k，l）。T為閾值，歸一化后的灰度圖像一般取T＝0.5，像素值G*（m，n）可被二值化為抖動的圖案b（m，n）， 其表達(dá)式為

針對二值化DCTb矩陣，矩陣元的取值范圍為［-1，1］，故要實現(xiàn)DCTb基調(diào)制，需要采用差分探測方案， 將式（9）中的 “0” 變?yōu)?“-1”， 與Walsh基相同。全測量時，若圖像大小為N×N，則全采樣的測量次數(shù)為2×N×N次。

1.4 測量基排序?qū)崿F(xiàn)壓縮采樣

在文獻(xiàn)［9］～文獻(xiàn)［14］中均提出了Walsh基壓縮測量的方案，并對壓縮率與圖像信噪比的關(guān)系進(jìn)行了研究。上述研究的本質(zhì)是對Walsh基的順序進(jìn)行重新排列從而選取測量系統(tǒng)較大的值來近似重建圖像，在犧牲了一部分信噪比的前提下提升了重建速度。WT和DCT的壓縮采樣與傳統(tǒng)的CS理論不同，WT和DCT的壓縮采樣是有損壓縮，而CS理論可實現(xiàn)無損壓縮。然而在實現(xiàn)效果和實用性上，WT和DCT測量基的壓縮測量在圖像探測、目標(biāo)識別等方面具有優(yōu)勢，其計算時間復(fù)雜度為O［Mlog2（N）］， 算法效率高于 CS 理論的倍數(shù)接近其迭代次數(shù)，M為壓縮測量次數(shù)。當(dāng)圖像越大，所需測量次數(shù)M越大，快速變換算法相對于CS理論的優(yōu)勢越明顯。

1.5 不同測量基對比分析

WT和DCT的測量基均為正交基，可通過快速變換實現(xiàn)單像素成像，兩者的區(qū)別在于：Walsh矩陣可通過1bit整數(shù)來表達(dá)，如在DMD的±12°兩個方向?qū)崿F(xiàn)光學(xué)測量；DCT測量基由連續(xù)的余弦函數(shù)演化而來，表達(dá)時需要用空間光調(diào)制器進(jìn)行量化，一般采用8bit整數(shù)來表達(dá)，如用DMD多次調(diào)制來產(chǎn)生灰度階的光強(qiáng)矩陣。在全采樣條件下，WT和DCT的測量基均得到了理想的成像效果，且WT測量基的測量速度是DCT的8倍。然而在壓縮采樣時，DCT測量基因其灰度的連續(xù)性優(yōu)勢，壓縮效率高于WT，故DCT成為JPEG圖像壓縮的標(biāo)準(zhǔn)。

然而，針對二值化后的DCT測量基，也可采用1bit調(diào)制矩陣元進(jìn)行成像，此時究竟哪一組測量基表現(xiàn)更佳、更適用于單像素成像尚無明確結(jié)論。為得出具有指導(dǎo)意義的結(jié)論，本文參考文獻(xiàn)［13］中的方法，對不同測量基的能量集中度、成像信噪比進(jìn)行了對比分析。圖像質(zhì)量評價標(biāo)準(zhǔn)選擇常用的兩種方法，即峰值信噪比（Peak Signal to Noise Ratio，PSNR）和結(jié)構(gòu)相似度（Structural Similarity， SSIM）。

2 實驗效果對比

2.1 數(shù)值仿真方法

采用Lena圖像作為測試對象。圖1為Lena圖像及相應(yīng)的WT和DCT權(quán)重，即經(jīng)相應(yīng)變換后的系數(shù)絕對值，圖像和權(quán)重矩陣均為128×128（像素）。

圖1中，WT系數(shù)和DCT系數(shù)均進(jìn)行了取絕對值操作，并歸一化到［0，255］區(qū)間，用于顯示和對比不同變換的權(quán)重。

在單像素成像數(shù)值仿真對比實驗中，首先對Walsh基和二值化后的DCT基進(jìn)行了權(quán)重排序，排序方法為：1）利用 STL-10數(shù)據(jù)庫［15］選取其中80000張灰度圖像，線性插值到128×128（像素）；2）對每一幅圖像k分別作WT和DCT，得到變換系數(shù)圖yk；3）對yk取絕對值，記為｜yk｜； 4）再分別對所有變換系數(shù)累加求和，即5）按累加結(jié)果對IB中每1個系數(shù)權(quán)值IB（i）由大到小進(jìn)行排序，最終將獲得通過圖像訓(xùn)練得出的測量基順序。訓(xùn)練得到的WT和DCT測量基順序與原始 測量基順序的對應(yīng)關(guān)系如圖2所示。

圖1 不同測量基在圖像變換后對應(yīng)的權(quán)重Fig.1 Corresponding weights of Lena after WT and DCT

圖2 訓(xùn)練排序與原序?qū)Ρ菷ig.2 Comparison between the trained order and the original order

圖2中，WT訓(xùn)練的測量基排序與DCT訓(xùn)練的測量基排序不同，反映出同一組圖片在不同測量基下的權(quán)重系數(shù)不同，這一點在圖1中也有體現(xiàn)。據(jù)此可以推測DCT和WT測量基對圖像的信息獲取效率將會有所差別，但哪種方式獲取的效率更高還不能明確得出結(jié)論，需要進(jìn)一步比較。

2.2 重建時間與能量集中度對比

圖2中的排序在不同壓縮采樣下對應(yīng)的采樣系數(shù)譜圖如圖3所示。由于DCT和WT對圖像特征的提取方式不同，故對應(yīng)最優(yōu)的采樣排序策略也不相同。因此，本文考慮了更為公平的比較。通過80000張圖像的訓(xùn)練后，相應(yīng)的DCT和WT系數(shù)將各自訓(xùn)練出從大到小的序列，且各自相應(yīng)的訓(xùn)練序列對于相應(yīng)測量基均為最優(yōu)排序。此時進(jìn)行成像結(jié)果比較，顯然更加客觀。

圖3 二值化DCT和WT在訓(xùn)練排序后的權(quán)重譜圖Fig.3 Sampling weights spectral diagram of binary DCT and WT after the trained order

從計算效率和能量采集效率兩個角度對比DCT和WT兩種方法的優(yōu)勢，定義采樣率符號為SR， 圖 3（a）～圖3（h）均為不同 SR 下得到的測量系數(shù)譜圖。 其中， 圖3（a）～圖3（d）為用二值化的DCT測量基按DCT訓(xùn)練的不同排序采樣SR（5%、15%、25%和50%）下對應(yīng)的測量系數(shù)譜圖，計算得出了相應(yīng)的圖像重建時間tDCT和能量集中度 ECR。圖 3（e）～圖 3（h）為 WT 測量基按 WT 訓(xùn)練的不同排序采樣SR（5%、15%、25%和50%）下對應(yīng)的測量系數(shù)譜圖，同樣計算得出了相應(yīng)的圖像重建時間tWT和能量集中度ECR。對比兩組數(shù)據(jù)可知，在不同采樣率下，無論DCT還是WT，重建時間基本保持不變，而ECR隨采樣率的增加而增大，符合理論預(yù)測；在相同采樣率下，DCT重建時間在0.8ms數(shù)量級，小于WT的1.4ms；采樣率為5%時，能量集中度DCT優(yōu)于WT，隨著采樣率的增加，能量集中度WT優(yōu)于DCT。

2.3 圖像PSNR與SSIM對比

將圖3中的權(quán)重譜圖按對應(yīng)的DCT和WT進(jìn)行重建， 得到的結(jié)果如圖4（a）～圖4（h）所示。 其中，圖4（a）～圖4（d）為 DCT 測量基按 DCT 訓(xùn)練的不同排序采樣SR（5%、15%、25%和50%）下對應(yīng)的重建圖像、圖像結(jié)構(gòu)相似度SSIM和峰值信噪比PSNR，圖 4（e）～圖4（h）為 WT 測量基按 WT 訓(xùn)練的不同排序采樣SR（5%、15%、25%和50%）下對應(yīng)的重建圖像、圖像結(jié)構(gòu)相似度SSIM和峰值信噪比PSNR。

圖4 二值化DCT和WT在訓(xùn)練排序后壓縮測量單像素成像結(jié)果Fig.4 Compress and measure single-pixel imaging results of binary DCT and WT after the trained order

由圖4可知，無論DCT還是WT，重建圖像的PSNR均隨采樣率SR的增加而不斷增大，圖像視覺效果也變得越來越好；DCT重建圖像的SSIM隨采樣率SR的增加也呈現(xiàn)上升趨勢，但在SR＝50%時出現(xiàn)拐點，與視覺判斷和PSNR趨勢不符。WT重建圖像的SSIM和PSNR均隨采樣率SR的增加而增大，圖像質(zhì)量變好。在低采樣率條件下，SR為5%和15%時，DCT重建圖像的PSNR和SSIM均優(yōu)于同樣采樣率下的WT；然而在SR為25%和50%時，情況發(fā)生了變化，WT重建圖像的PSNR和SSIM均優(yōu)于同樣采樣率下的DCT，視覺效果與圖像評價值判斷結(jié)果相符。參考圖3的能量集中度ECR，單從Lena圖像的分析可知，在低采樣率條件下，DCT方法占優(yōu)；而在高采樣率條件下，WT方法明顯優(yōu)于DCT。

2.4 重建圖像評價值統(tǒng)計分析

擴(kuò)大測試圖像數(shù)據(jù)集，用統(tǒng)計的方法對重建圖像的PSNR、SSIM和ECR進(jìn)行量化對比分析，使得研究結(jié)論具有普適性。測試方法和重建過程與Lena圖像研究方法完全相同，不同之處在于，測試圖像樣本數(shù)為隨機(jī)從STL-10圖庫中抽取的500張圖像，抽取圖像與訓(xùn)練集中的80000幅圖像無關(guān)。具體的實驗步驟為：1）對500幅圖像進(jìn)行插值， 從 96×96（像素）變?yōu)?128×128（像素）； 2）利用DCT和WT分別對每1幅圖像進(jìn)行測量，測量基按訓(xùn)練權(quán)重排序，每增加1次測量，進(jìn)行1次圖像恢復(fù)，并進(jìn)行ECR、PSNR和SSIM的計算，得到相應(yīng)SR下的數(shù)值；3）重復(fù)對每1幅圖像進(jìn)行操作，并增加測量次數(shù)直到全采樣，共進(jìn)行16384次測量。按上述方法得到的數(shù)據(jù)如圖5、圖6所示。

圖5 WT和DCT重建圖像對應(yīng)的SSIM和ECR與測量次數(shù)的關(guān)系Fig.5 Relationship between SSIM，ECR corresponding to WT and DCT reconstruction image and measurement times

圖6 WT和DCT重建圖像對應(yīng)的PSNR與測量次數(shù)的關(guān)系Fig.6 Relationship between PSNR corresponding to WT and DCT reconstruction image and measurement times

圖5表明，抖動二值化DCT重建圖像的SSIM在測量次數(shù)為2458次左右、對應(yīng)的采樣率SR約為15%（2458／16384）時，其測量值優(yōu)于同樣條件下的WT。而WT的SSIM在采樣率SR＞15%（測量次數(shù)＞2458次）時，其效果優(yōu)于二值化的DCT成像結(jié)果。對于ECR，二值化的DCT在低采樣率條件下也有相同趨勢，即DCT的ECR優(yōu)于WT；而在高采樣率條件下，WT相比DCT更加具有優(yōu)勢。值得注意的是，二值化DCT方法的SSIM并不隨測量次數(shù)的增加而增大，而是先增大后減小；而WT方法則隨著測量次數(shù)的增加而增大。

進(jìn)一步結(jié)合圖像的PSNR對二值化DCT和WT重建圖像進(jìn)行分析，數(shù)據(jù)結(jié)果如圖6所示。基于二值化DCT在壓縮采樣時重建圖像的PSNR在采樣率較低時（接近10%，對應(yīng)的測量次數(shù)為1638），DCT測量值優(yōu)于同樣條件下的WT；隨著采樣率的增加，二值化DCT方法的PSNR先上升后下降，而WT方法的PSNR則隨著采樣率的增加而增大。通過與圖5對比發(fā)現(xiàn)，PSNR、SSIM和ECR給出了一致的結(jié)果，即抖動二值化DCT在低采樣率條件下的重建圖像相比WT具有一定優(yōu)勢，PSNR在10%采樣率（測量次數(shù)為1638）下均可達(dá)20dB；當(dāng)采樣率增大到50%（測量次數(shù)8192）左右，重建圖像的SSIM可優(yōu)于0.9，PSNR優(yōu)于30dB，此時WT方法顯著優(yōu)于二值化DCT方法。

值得特別說明的是，二值化DCT方法重建圖像的SSIM和PSNR曲線均隨測量次數(shù)的增加先上升后下降，這說明隨著測量次數(shù)的持續(xù)增加，重建圖像的效果反而在下降，似乎與理論不符。分析其原因，由于對DCT測量基進(jìn)行了抖動二值化，引入了量化誤差，測量增加后編碼空間網(wǎng)格變密，二值化的量化近似誤差隨之增大。因此，測量基的誤差導(dǎo)致了重建誤差，相關(guān)討論可參考文獻(xiàn)［16］～文獻(xiàn)［18］。

綜上所述，采樣率小于10%時，二值化DCT方法重建結(jié)果優(yōu)于WT方法；采樣率在10%～20%時，WT方法的PSNR與二值化的DCT方法接近；在采樣率大于20%時，WT方法的PSNR優(yōu)于二值化DCT方法；相同采樣率下，DCT圖像重建時間略優(yōu)于WT，但二者在同一數(shù)量級；在采樣率低時，DCT方法的ECR優(yōu)于WT方法，表明對于平滑圖像，DCT方法更加有效。

3 結(jié)論

本文開展了基于WT和DCT測量基的單像素成像對比研究。針對抖動二值化DCT矩陣和WT方法，分別在不同采樣率條件下對單像素成像的PSNR和SSIM、圖像重建時間t和能量集中度ECR進(jìn)行了研究。通過數(shù)值仿真對80000幅圖像訓(xùn)練權(quán)重進(jìn)行了排序，對壓縮采樣方案進(jìn)行了優(yōu)化，詳細(xì)分析了Lena圖像在抖動二值化DCT測量基和WT測量基的權(quán)重譜圖和重建結(jié)果，并從STL-10圖庫中隨機(jī)抽取500幅圖像進(jìn)行了量化分析。

在研究方法上，本文提出了圖像訓(xùn)練權(quán)重排序方法，從新角度研究了測量基排序后壓縮采樣與單像素成像圖像質(zhì)量的關(guān)系。通過對DCT圖像進(jìn)行抖動二值化，使其可在DMD上快速實現(xiàn)，并在此基礎(chǔ)上與WT方法進(jìn)行了對比研究，研究方法具有指導(dǎo)意義。通過大量數(shù)值仿真，可以得到：采樣率小于10%時，二值化DCT方法重建結(jié)果優(yōu)于WT方法；采樣率在10%～20%時，WT方法與二值化DCT方法效果接近；采樣率大于20%時，WT方法重建圖像質(zhì)量優(yōu)于二值化DCT方法。在相同采樣率下，DCT方法的圖像重建時間略優(yōu)于WT方法，二者在同一數(shù)量級。該結(jié)論對實現(xiàn)快速、高PSNR單像素成像的測量基選取具有較好的參考價值。下一步，針對WT和DCT的壓縮采樣單像素成像仍有大量工作值得研究，包括針對不同圖像先分類再進(jìn)行優(yōu)化排序、排序效果在實驗系統(tǒng)的對比驗證等，研究成果將進(jìn)一步推進(jìn)單像素成像技術(shù)的實用化。