李明飛, 袁梓豪, 趙琳琳, 孫曉潔
(1.中國航天科技集團有限公司量子工程研究中心,北京100094;2.北京航天控制儀器研究所,北京100039)
壓縮感知(壓縮采樣)(Compressive Sampling,CS)理論由 Candès等于 2006 年正式提出[1-2], 美國Rice大學在2008年利用單像素成像實驗驗證了其有效性。CS理論是對稀疏或可壓縮的信號進行少量非適應性的線性測量,可以近似100%地恢復出原始信號,可突破Shannon采樣定理的限制,同時還具有超靈敏探測等優勢[3]。目前,用于單像素成像壓縮感知算法的測量次數僅為Nyquist采樣極限的1%左右,即可重建出目標圖像[4]。壓縮感知算法是非線性迭代算法,且重建效果與所選稀疏基密切相關,故成像計算時間和穩定性同時受到極大挑戰。而強度關聯算法具有硬件要求低、成像速度快的特點,但其需要大于Nyquist采樣極限的測量次數,故采樣時間較長。
基于Walsh變換(Walsh Transform,WT)的方法提供了快速成像方案,該方案存在快速變換算法且硬件易于實現。采用二值矩陣作為調制矩陣,易于在高速數字微鏡器件(Digital Micro-mirror Device,DMD)上實現。調制矩陣生成速度快且無需存儲,Walsh矩陣作為測量矩陣具有正交特性,可完美重建圖像。基于離散余弦變換(Discrete Cosine Transform,DCT)的方法也提供了實現快速單像素成像的方案,除了具有快速變換的特點和正交性的優點外,DCT基是連續的,更適用于圖像壓縮,這是也聯合圖像專家組(Joint Photographic Experts Group,JPEG)采用DCT作為圖像壓縮方案的原因。當然,DCT相比于WT的重要區別還在于DCT的連續性,在用于單像素成像時不利于實現快速光學調制,如DCT無法直接使用DMD實現20kHz的二值調制。文獻[5]報道的抖動二值化方法使得DCT這類連續測量基實現高速調制測量成為現實,研究發現二值化后的DCT基與Walsh基非常相似,然而在單像素成像中二者誰更占優仍未有人研究。因此,在單像素成像中Walsh基和二值化后的DCT基各自優化排序后,哪一種成像信噪比更高和重建速度更快仍是值得研究的課題。
針對上述問題,本文開展了對基于Walsh基和二值化后的DCT基的單像素成像數值仿真實驗對比研究。對DCT矩陣進行了抖動二值化,分別對比了在不同采樣率條件下的基于兩種變換的單像素成像信噪比和成像時間,對比研究了相同采樣率下測量基三種不同排序時的單像素成像信噪比。
Walsh矩陣和Hadamard矩陣具有相同的矩陣元素,區別僅在于構造矩陣的行和列排列不同。Hadamard矩陣起源較早,應用極其廣泛[6-7]。Hadamard矩陣(簡稱H矩陣)具有如下特征:1)矩陣產生速度快;2)矩陣元僅有1元素。矩陣任意兩行(或列)正交, 有

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

式(2)中,b(k,l,p,q)為對 2 取模運算, 其表達式為

式(3)中,g0(p)=pn-1,g1(p)=pn-1+pn-2,g2(p)=pn-2+pn-3, …,gn-1(p)=p1+p0。 下標i=0, 1,2, …,n-1分別對應ki、li、pi、qi二進制表示時相應的位。由于Walsh矩陣元為-1和1兩種元素,需采用差分探測的方式,即先將矩陣元-1變為1,1變為0,構成D-1矩陣;將矩陣元-1變為0,1不變,產生D+1矩陣,通過利用D+1-D-1的方式來實現Walsh矩陣完整的測量。理論上全測量時,若圖像大小為N×N,則全采樣的測量次數為2×N×N,當采用互補調制[8]的方式時,實際測量次數需要N×N次。
二維DCT與離散Fourier變換(Discrete Fourier Transform,DFT)的數學表達式非常接近,實際上DCT就是DFT只截取余弦變換部分生成。二維DCT的數學表達式如下

其中,M和N為物體的行數和列數,0≤m≤M-1, 0≤n≤N-1。
DCT矩陣和Walsh矩陣的測量基相似,但重要區別在于DCT是具有灰度的測量基。若采用DMD為調制器,按8bit灰度投影計算,DCT基測量速度是Walsh基測量速度的1/8,相應的單像素成像速度也將是Walsh基方案的1/8,全采樣時二者信噪比相同。
經上述分析,從實際參考價值上看,編碼同為灰度編碼或同為二值編碼的比較才能具有實際意義,直接使用DCT的灰度測量基與Walsh基比較則意義不大。
張子 邦 等[5,8]引 入 了 抖 動 算 法, 將 快 速Fourier變換(Fast Fourier Transform, FFT)有灰度的正弦圖案轉換成二值正弦圖案,提高了基于FFT測量基的調制速度,使得單像素成像速度大幅提升。本文引入了抖動算法,將DCT基進行二值化,記為DCTb。不同于文獻[5]和文獻[8]之處在于,本文的DCT抖動二值化不對測量基進行上采樣或插值,即抖動前后的矩陣均與原始矩陣大小相等,這樣做的優勢在于在實際應用中更加接近真實情況。
前文提到的抖動算法的核心思想是基于誤差擴散原理,首先對當前像素值進行甄別,采用閾值法將當前像素值二值化后輸出,然后將輸入和輸出的像素值之差按一定的權重傳播到若干相鄰未處理的區域。誤差擴散原理可表達為

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

針對二值化DCTb矩陣,矩陣元的取值范圍為[-1,1],故要實現DCTb基調制,需要采用差分探測方案, 將式(9)中的 “0” 變為 “-1”, 與Walsh基相同。全測量時,若圖像大小為N×N,則全采樣的測量次數為2×N×N次。
在文獻[9]~文獻[14]中均提出了Walsh基壓縮測量的方案,并對壓縮率與圖像信噪比的關系進行了研究。上述研究的本質是對Walsh基的順序進行重新排列從而選取測量系統較大的值來近似重建圖像,在犧牲了一部分信噪比的前提下提升了重建速度。WT和DCT的壓縮采樣與傳統的CS理論不同,WT和DCT的壓縮采樣是有損壓縮,而CS理論可實現無損壓縮。然而在實現效果和實用性上,WT和DCT測量基的壓縮測量在圖像探測、目標識別等方面具有優勢,其計算時間復雜度為O[Mlog2(N)], 算法效率高于 CS 理論的倍數接近其迭代次數,M為壓縮測量次數。當圖像越大,所需測量次數M越大,快速變換算法相對于CS理論的優勢越明顯。
WT和DCT的測量基均為正交基,可通過快速變換實現單像素成像,兩者的區別在于:Walsh矩陣可通過1bit整數來表達,如在DMD的±12°兩個方向實現光學測量;DCT測量基由連續的余弦函數演化而來,表達時需要用空間光調制器進行量化,一般采用8bit整數來表達,如用DMD多次調制來產生灰度階的光強矩陣。在全采樣條件下,WT和DCT的測量基均得到了理想的成像效果,且WT測量基的測量速度是DCT的8倍。然而在壓縮采樣時,DCT測量基因其灰度的連續性優勢,壓縮效率高于WT,故DCT成為JPEG圖像壓縮的標準。
然而,針對二值化后的DCT測量基,也可采用1bit調制矩陣元進行成像,此時究竟哪一組測量基表現更佳、更適用于單像素成像尚無明確結論。為得出具有指導意義的結論,本文參考文獻[13]中的方法,對不同測量基的能量集中度、成像信噪比進行了對比分析。圖像質量評價標準選擇常用的兩種方法,即峰值信噪比(Peak Signal to Noise Ratio,PSNR)和結構相似度(Structural Similarity, SSIM)。
采用Lena圖像作為測試對象。圖1為Lena圖像及相應的WT和DCT權重,即經相應變換后的系數絕對值,圖像和權重矩陣均為128×128(像素)。
圖1中,WT系數和DCT系數均進行了取絕對值操作,并歸一化到[0,255]區間,用于顯示和對比不同變換的權重。
在單像素成像數值仿真對比實驗中,首先對Walsh基和二值化后的DCT基進行了權重排序,排序方法為:1)利用 STL-10數據庫[15]選取其中80000張灰度圖像,線性插值到128×128(像素);2)對每一幅圖像k分別作WT和DCT,得到變換系數圖yk;3)對yk取絕對值,記為|yk|; 4)再分別對所有變換系數累加求和,即5)按累加結果對IB中每1個系數權值IB(i)由大到小進行排序,最終將獲得通過圖像訓練得出的測量基順序。訓練得到的WT和DCT測量基順序與原始 測量基順序的對應關系如圖2所示。

圖1 不同測量基在圖像變換后對應的權重Fig.1 Corresponding weights of Lena after WT and DCT

圖2 訓練排序與原序對比Fig.2 Comparison between the trained order and the original order
圖2中,WT訓練的測量基排序與DCT訓練的測量基排序不同,反映出同一組圖片在不同測量基下的權重系數不同,這一點在圖1中也有體現。據此可以推測DCT和WT測量基對圖像的信息獲取效率將會有所差別,但哪種方式獲取的效率更高還不能明確得出結論,需要進一步比較。
圖2中的排序在不同壓縮采樣下對應的采樣系數譜圖如圖3所示。由于DCT和WT對圖像特征的提取方式不同,故對應最優的采樣排序策略也不相同。因此,本文考慮了更為公平的比較。通過80000張圖像的訓練后,相應的DCT和WT系數將各自訓練出從大到小的序列,且各自相應的訓練序列對于相應測量基均為最優排序。此時進行成像結果比較,顯然更加客觀。

圖3 二值化DCT和WT在訓練排序后的權重譜圖Fig.3 Sampling weights spectral diagram of binary DCT and WT after the trained order
從計算效率和能量采集效率兩個角度對比DCT和WT兩種方法的優勢,定義采樣率符號為SR, 圖 3(a)~圖3(h)均為不同 SR 下得到的測量系數譜圖。 其中, 圖3(a)~圖3(d)為用二值化的DCT測量基按DCT訓練的不同排序采樣SR(5%、15%、25%和50%)下對應的測量系數譜圖,計算得出了相應的圖像重建時間tDCT和能量集中度 ECR。圖 3(e)~圖 3(h)為 WT 測量基按 WT 訓練的不同排序采樣SR(5%、15%、25%和50%)下對應的測量系數譜圖,同樣計算得出了相應的圖像重建時間tWT和能量集中度ECR。對比兩組數據可知,在不同采樣率下,無論DCT還是WT,重建時間基本保持不變,而ECR隨采樣率的增加而增大,符合理論預測;在相同采樣率下,DCT重建時間在0.8ms數量級,小于WT的1.4ms;采樣率為5%時,能量集中度DCT優于WT,隨著采樣率的增加,能量集中度WT優于DCT。
將圖3中的權重譜圖按對應的DCT和WT進行重建, 得到的結果如圖4(a)~圖4(h)所示。 其中,圖4(a)~圖4(d)為 DCT 測量基按 DCT 訓練的不同排序采樣SR(5%、15%、25%和50%)下對應的重建圖像、圖像結構相似度SSIM和峰值信噪比PSNR,圖 4(e)~圖4(h)為 WT 測量基按 WT 訓練的不同排序采樣SR(5%、15%、25%和50%)下對應的重建圖像、圖像結構相似度SSIM和峰值信噪比PSNR。

圖4 二值化DCT和WT在訓練排序后壓縮測量單像素成像結果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的增加也呈現上升趨勢,但在SR=50%時出現拐點,與視覺判斷和PSNR趨勢不符。WT重建圖像的SSIM和PSNR均隨采樣率SR的增加而增大,圖像質量變好。在低采樣率條件下,SR為5%和15%時,DCT重建圖像的PSNR和SSIM均優于同樣采樣率下的WT;然而在SR為25%和50%時,情況發生了變化,WT重建圖像的PSNR和SSIM均優于同樣采樣率下的DCT,視覺效果與圖像評價值判斷結果相符。參考圖3的能量集中度ECR,單從Lena圖像的分析可知,在低采樣率條件下,DCT方法占優;而在高采樣率條件下,WT方法明顯優于DCT。
擴大測試圖像數據集,用統計的方法對重建圖像的PSNR、SSIM和ECR進行量化對比分析,使得研究結論具有普適性。測試方法和重建過程與Lena圖像研究方法完全相同,不同之處在于,測試圖像樣本數為隨機從STL-10圖庫中抽取的500張圖像,抽取圖像與訓練集中的80000幅圖像無關。具體的實驗步驟為:1)對500幅圖像進行插值, 從 96×96(像素)變為 128×128(像素); 2)利用DCT和WT分別對每1幅圖像進行測量,測量基按訓練權重排序,每增加1次測量,進行1次圖像恢復,并進行ECR、PSNR和SSIM的計算,得到相應SR下的數值;3)重復對每1幅圖像進行操作,并增加測量次數直到全采樣,共進行16384次測量。按上述方法得到的數據如圖5、圖6所示。

圖5 WT和DCT重建圖像對應的SSIM和ECR與測量次數的關系Fig.5 Relationship between SSIM,ECR corresponding to WT and DCT reconstruction image and measurement times

圖6 WT和DCT重建圖像對應的PSNR與測量次數的關系Fig.6 Relationship between PSNR corresponding to WT and DCT reconstruction image and measurement times
圖5表明,抖動二值化DCT重建圖像的SSIM在測量次數為2458次左右、對應的采樣率SR約為15%(2458/16384)時,其測量值優于同樣條件下的WT。而WT的SSIM在采樣率SR>15%(測量次數>2458次)時,其效果優于二值化的DCT成像結果。對于ECR,二值化的DCT在低采樣率條件下也有相同趨勢,即DCT的ECR優于WT;而在高采樣率條件下,WT相比DCT更加具有優勢。值得注意的是,二值化DCT方法的SSIM并不隨測量次數的增加而增大,而是先增大后減小;而WT方法則隨著測量次數的增加而增大。
進一步結合圖像的PSNR對二值化DCT和WT重建圖像進行分析,數據結果如圖6所示。基于二值化DCT在壓縮采樣時重建圖像的PSNR在采樣率較低時(接近10%,對應的測量次數為1638),DCT測量值優于同樣條件下的WT;隨著采樣率的增加,二值化DCT方法的PSNR先上升后下降,而WT方法的PSNR則隨著采樣率的增加而增大。通過與圖5對比發現,PSNR、SSIM和ECR給出了一致的結果,即抖動二值化DCT在低采樣率條件下的重建圖像相比WT具有一定優勢,PSNR在10%采樣率(測量次數為1638)下均可達20dB;當采樣率增大到50%(測量次數8192)左右,重建圖像的SSIM可優于0.9,PSNR優于30dB,此時WT方法顯著優于二值化DCT方法。
值得特別說明的是,二值化DCT方法重建圖像的SSIM和PSNR曲線均隨測量次數的增加先上升后下降,這說明隨著測量次數的持續增加,重建圖像的效果反而在下降,似乎與理論不符。分析其原因,由于對DCT測量基進行了抖動二值化,引入了量化誤差,測量增加后編碼空間網格變密,二值化的量化近似誤差隨之增大。因此,測量基的誤差導致了重建誤差,相關討論可參考文獻[16]~文獻[18]。
綜上所述,采樣率小于10%時,二值化DCT方法重建結果優于WT方法;采樣率在10%~20%時,WT方法的PSNR與二值化的DCT方法接近;在采樣率大于20%時,WT方法的PSNR優于二值化DCT方法;相同采樣率下,DCT圖像重建時間略優于WT,但二者在同一數量級;在采樣率低時,DCT方法的ECR優于WT方法,表明對于平滑圖像,DCT方法更加有效。
本文開展了基于WT和DCT測量基的單像素成像對比研究。針對抖動二值化DCT矩陣和WT方法,分別在不同采樣率條件下對單像素成像的PSNR和SSIM、圖像重建時間t和能量集中度ECR進行了研究。通過數值仿真對80000幅圖像訓練權重進行了排序,對壓縮采樣方案進行了優化,詳細分析了Lena圖像在抖動二值化DCT測量基和WT測量基的權重譜圖和重建結果,并從STL-10圖庫中隨機抽取500幅圖像進行了量化分析。
在研究方法上,本文提出了圖像訓練權重排序方法,從新角度研究了測量基排序后壓縮采樣與單像素成像圖像質量的關系。通過對DCT圖像進行抖動二值化,使其可在DMD上快速實現,并在此基礎上與WT方法進行了對比研究,研究方法具有指導意義。通過大量數值仿真,可以得到:采樣率小于10%時,二值化DCT方法重建結果優于WT方法;采樣率在10%~20%時,WT方法與二值化DCT方法效果接近;采樣率大于20%時,WT方法重建圖像質量優于二值化DCT方法。在相同采樣率下,DCT方法的圖像重建時間略優于WT方法,二者在同一數量級。該結論對實現快速、高PSNR單像素成像的測量基選取具有較好的參考價值。下一步,針對WT和DCT的壓縮采樣單像素成像仍有大量工作值得研究,包括針對不同圖像先分類再進行優化排序、排序效果在實驗系統的對比驗證等,研究成果將進一步推進單像素成像技術的實用化。