基于單像素相機重構矩陣優化的影像采集和重構方法

程 濤，劉玉安

（1.廣西科技大學汽車與交通學院，廣西 柳州 545006；2.武漢大學測繪遙感信息工程國家重點實驗室，湖北 武漢 430079）

0 引言

壓縮感知理論表明，如果信號是稀疏的就能以遠低于奈奎斯特（Nyquist）采樣頻率的采樣率采集信號，并能以高概率精確重建信號。但是，影像等信號一般都不是稀疏的，而是可壓縮的。因此，需要對影像等可壓縮信號作稀疏變換以適應壓縮感知的要求。可壓縮信號的壓縮感知模型如式（1）所示。

式（1）中：y是測量數據，y∈RM；Φ 是測量矩陣，Φ∈RM×N；Ψ 是稀疏變換基，Ψ ∈RN×N；x是可壓縮信號，x∈RN，x＝ΨTα；α為x的稀疏變換域系數，α∈RN，α＝Ψx；重構矩陣ΠR＝ΦΨT。

壓縮感知的數據采集處理過程可分為測量和重構兩個階段。在測量階段通過測量矩陣Φ以y＝Φx采集測量數據y；在重構階段，求解式（1）得到變換域系數，并通過＝取得重構信號。

對于可壓縮信號，如果測量矩陣和稀疏變換基的不相關性好，那么重構矩陣的信號重構能力就強［2］。在大家尋找與稀疏變換基不相關性好的測量矩陣時，Elad［3］卻提出了截然不同的測量矩陣設計思路。Elad［3］和Duarte［4］等在稀疏變換基已確定的情況下，以重構矩陣各列相關性最小化或平均化為目標優化重構矩陣；然后，通過稀疏變換基和優化后的重構矩陣（簡稱優化矩陣）解算出與之相應的測量矩陣。盡管優化矩陣有更好的重構能力，但是測量矩陣不再是原測量矩陣。即使原測量矩陣是便于電路設計的確定性矩陣或便于硬件實現的單像素相機采用的0－1隨機矩陣，測量矩陣也不會再保持原測量矩陣的優良性質。Elad等人的重構矩陣優化算法盡管能夠事前確定稀疏變換基，卻無法事前確定測量矩陣，也就無法改變當前測量矩陣難以硬件實現和設計的窘境，因而應用價值有限。

稀疏信號無需作稀疏變換，所以這時測量矩陣就是重構矩陣。文獻［2，5—6］通過測量矩陣的行向量正交化以及列向量單位化來優化測量矩陣。優化后的測量矩陣盡管能改善信號重構效果，但是和Elad的方法一樣，優化后的測量矩陣性質存在不確定性和未知性，也存在難以硬件實現的問題，不適合工業化應用。

文獻［1］針對稀疏信號提出了測量矩陣和測量數據的事后優化算法和模型。實現了在測量階段采用事前確定的測量矩陣（0－1稀疏矩陣）采集測量數據，在重構階段采用優化后的測量矩陣（高斯矩陣）重構稀疏信號，取得了很好的效果，并給出了相關的理論分析和證明。

如能將文獻［1］的測量矩陣分離算法推廣應用于可壓縮信號，就可克服Elad等人重構矩陣優化算法的缺點。在稀疏變換基確定的情況下，實現以事前確定的測量矩陣采集數據，以事后優化的重構矩陣重構信號。從而初步解決當前測量矩陣難以硬件實現和設計的難題，并且保持高的信號重構能力。

離 散 余 弦 變 換 （Discrete cosine transform，DCT）編碼是影像／視頻壓縮采用的主流技術［7］。單像素相機的0－1隨機矩陣不但易于硬件實現而且所需存儲空間小運算速度快。因此，以離散余弦變換矩陣作為式（1）的稀疏變換基，以單像素相機的0－1隨機矩陣作為式（1）的測量矩陣。為便于與同行研究成果作對比和硬件設計，測量矩陣的大小為128×256，各行向量包含32個隨機均勻分布的1。

盡管單像素相機測量矩陣（0－1隨機矩陣）和稀疏變換基（DCT矩陣）的不相關性很差，但是理論和實驗證明，基于單像素相機重構矩陣優化的方法依然能取得好的信號重構效果。從而回避了Donoho等人［2］提出的“測量矩陣和稀疏變換基不相關性好”的測量矩陣設計原則。本文針對壓縮感知在遙感中的真實應用情境，提出了無需把二維影像轉化成一維信號的基于線陣推掃和重構矩陣優化的單像素相機影像采集和重構方法。

1 重構矩陣的優化算法和近似矩陣的計算公式

文獻［1］的測量矩陣優化算法，只適用于稀疏信號，可使優化后的測量矩陣具備行模相近、列模相近、各行不相關性好和各列不相關性好以及服從高斯分布5個特征。影像等可壓縮信號作稀疏變換后，可近似視為稀疏信號。因此，在文獻［1］的測量矩陣優化算法中引入稀疏變換基，即可變成適用于可壓縮信號的重構矩陣優化方法，如算法1所示。

算法1：

輸入：測量矩陣Φ和稀疏變換基Ψ初始化：ΠR＝Φ ΨT，設定迭代次數i的初始值為0，最大值imax默認為100迭代：在第i次迭代執行以下步驟

1）以Jarque－Bera檢驗計算ΠR各列和各行服從高斯分布的行數Jri和列數Jci；

2）計算ΠR各列向量間的相關系數，取出其絕對值的最大值μcimax；計算各行向量間的相關系數，取出其絕對值的最大值μrimax；

3）計算ΠR各行向量的模，取出其最大值normrimax和最小值 normri min；

4）正交規范化ΠR各行向量，單位化ΠR各列向量；

5）使i＝i＋1，判斷i＜imax，如果是退出迭代，否則返回執行步驟1；

輸出：優化后的重構矩陣ΠO。

因為算法1也是以行變換為主的算法，所以只要在文獻［1］的近似矩陣計算公式中引入稀疏變換基，把測量矩陣換成針對可壓縮信號的重構矩陣就可得到重構矩陣的近似矩陣計算公式，如式（2）所示。

表1 3類矩陣的相關參數Tab.1 Parameters of 3matrices

表1列出了測量矩陣采用單像素相機的0－1隨機矩陣，稀疏變換基采用DCT矩陣的算法1迭代100步時重構矩陣ΠR、優化矩陣ΠO和近似矩陣ΠT的各種統計學參數（表1中各行數據，除“Jarque－Bera檢驗”外，“／”號左邊的數據表示最小值，“／”號右邊的數據表示最大值；“Jarque－Bera檢驗”行，“／”號左邊的數據表示符合檢驗的列數；“／”號右邊的數據表示符合檢驗的行數）。μcmax表示列間相關系數絕對值的最大值；μrmax表示行間相關系數絕對值的最大值。因為采用的測量矩陣是各行都包含32個1的0－1隨機矩陣，而且DCT矩陣是正交矩陣，所以重構矩陣的行模都為

2 重構矩陣優化算法的迭代過程分析

根據文獻［1］的理論和實驗分析，只要優化矩陣和近似矩陣同時具備行模相近、列模相近、各行不相關性好和各列不相關性好就會對高斯稀疏信號具有很好的重構能力，如果優化矩陣和近似矩陣各行和各列幾乎都滿足高斯分布就會同時具備高斯矩陣對各類稀疏信號的普適性。因此本文以此為判據分析圖1—圖4所表示的算法1對重構矩陣的迭代優化過程和相應的近似矩陣的性質變化。圖1和圖4中，由于迭代0處的y值過大，故只能部分顯示，0處的y值可由表1查得。

圖1 矩陣各行模的極值與迭代數的關系Fig.1 Row norms extreme value of matrix vs.iteration

由圖2（a）可見優化矩陣的各行列相關系數絕對值的最大值在迭代過程中迅速收斂，在第13次迭代后數值不再變化，曲線變成直線。由表1可見，優化矩陣的行相關系數為0，即各行完全正交；優化矩陣的列相關系數也由0.839減小到0.364。由圖2（b）可見近似矩陣的各列相關系數絕對值的最大值在迭代過程中也迅速變小收斂，行相關系數最大值到第70次迭代后才基本穩定。由表1可見，近似矩陣的行相關系數為0.02，各行幾乎完全正交；近似矩陣的列相關系數也由0.839減小到0.279。

圖2 矩陣各行列相關系數絕對值的最大值與迭代的關系Fig.2 Maximum of absolute value of correlation coefficient of matrix rows and columns vs.iteration

由圖3（a）可見優化矩陣的各行列在第60次迭代后基本都服從高斯分布（第100次迭代后各行列服從高斯分布的數量各為247和123，如表1所示），說明重構矩陣已演變為高斯矩陣。圖3（b）和（a）的曲線構型非常相似，也是在第60次迭代后各行列基本都服從高斯分布（第100次迭代后各行列服從高斯分布的數量各為246和123，如表1所示）。由圖4可見近似矩陣的列模極值收斂到與1很接近的值。由于算法1每次迭代的最后一步都是“單位化ΠR各列向量”，所以優化矩陣各列模的值都為1，如表1所示。

由表1和圖1－圖4可見，優化矩陣具備各行不相關、各列不相關性好、各行模相等和各列模相等4個優點，而且優化矩陣各行和各列幾乎都服從高斯分布。基于算法1和式（2）的近似矩陣很好地繼承了優化矩陣的性質。因此，文獻［1］的測量矩陣分離算法可推廣應用于適用于可壓縮信號的重構矩陣。

圖3 矩陣服從高斯分布行列數與迭代的關系Fig.3 Number of rows and columns which follows the gauss distribution vs.iteration

圖4 近似矩陣各列模極值與迭代數的關系Fig.4 coumn norm extreme value of approximate matrix vs.iteration

3 重構矩陣優化算法的理論分析

進一步通過Kolmogorov－Smirnov檢驗可以發現，優化矩陣不僅服從高斯分布，而且服從N（0，1／M）分布。因此，優化矩陣各行模的平方服從χ2分布，即各行模的數學期望為。由圖1（a）可見，第16次迭代后，行模極值都收斂于，與其數學期望值相同。

盡管測量矩陣是各行都包含32個隨機均勻分布的1的0－1隨機矩陣，但依然可視為各元素滿足獨立同分布性質。由表1可見，由測量矩陣和DCT矩陣（稀疏變換基）乘積所得的重構矩陣，通過Jarque－Bera檢驗服從高斯分布的行列數分別為121和247。因此，重構矩陣可視為高斯矩陣。

算法1中行正交化操作盡管涉及多個行的運算，根據線性變換不變性原理可知正交化后的各行和各列依然符合高斯分布。根據大數定律可知，在多次迭代過程中，各行和各列的高斯分布越來越接近于樣本均值為0和樣本方差為1／M的理想狀態，服從N（0，1／M）分布的行列數分別為120和243。因此，優化矩陣各元素可視為來自服從標準高斯分布的樣本總體。近似矩陣與優化矩陣的性質相近，近似矩陣各元素也可視為來自服從標準高斯分布的樣本總體，服從N（0，1／M）分布的行列數分別為121和245。

對于不同的重構矩陣，算法1中迭代次數最大值的選取要滿足如下條件：能夠使圖1、圖2和圖4中優化矩陣和近似矩陣的曲線隨迭代收斂穩定，滿足行模相近、列模相近、各行不相關性好和各列不相關性好4個優點；能夠使圖3中優化矩陣和近似矩陣服從高斯分布的行列數大于矩陣行列數的95%（即使隨機生成的高斯矩陣也無法保證95%的行列數服從高斯分布）。

4 基于稀疏信號重構實驗的重構矩陣優化算法的驗證

盡管ΠR和ΠT是適用于可壓縮信號的包含稀疏變換基的矩陣，但為了驗證算法1和式（2）的有效性，分別以重構矩陣ΠR、優化矩陣ΠO和重構矩陣ΠT直接測量并重構高斯稀疏信號和0－1稀疏信號，檢驗3類矩陣的性能，實驗結果如圖5所示（采用Matlab在普通臺式機上做仿真實驗，CPU主頻2 GHz，內存1GB）。圖5以正交匹配追蹤算法（orthogonal matching pursuit，OMP）對每個稀疏度K值重復試驗500次，然后計算稀疏信號的精確重構概率。為更清晰地表達實驗結果，圖5只表示各曲線重構概率為1的最后一個點和重構概率為0的第一個點，以及兩點之間的曲線部分。

由圖5可見，優化矩陣和近似矩陣的重構效果遠遠好于重構矩陣，整條曲線都位于重構矩陣的右側。優化矩陣和近似矩陣的重構效果基本一致，兩條曲線幾乎重合。說明算法1和式（2）都是非常有效可行的。因為0－1隨機矩陣行列不相關性差［8－9］，DCT矩陣也不是高斯矩陣，所以兩者構成的重構矩陣的行列不相關性不好（如表1所示），由圖5可見重構矩陣的重構效果極差。

圖5 類矩陣的重構概率與稀疏度的關系Fig.5 Prob.of exact recovery vs.the sparsity by 3matrices

5 基于重構矩陣優化算法的單像素相機影像采集重構方法的驗證

基于算法1和式（2）的近似矩陣計算方法可使近似矩陣很好地繼承優化矩陣的性質。因此，能以過渡矩陣Τ實現重構矩陣ΠR和測量數據y的事后優化。在測量階段采用易于硬件實現的0－1隨機矩陣，以y＝Φx采集數據；在重構階段以過渡矩陣Τ將測量數據y和重構矩陣ΠR完善優化：yT＝Τy，ΠT＝ΤΠR，并以式（3）重構可壓縮信號x的稀疏變換域系數α。最后，通過就能得到重構的可壓縮信號

文獻［10］等的壓縮感知影像處理多是從原始影像中抽取16×16的子塊，然后逐列首尾相接成向量。由于子塊包含區域有限，均勻性好，因此能獲得較好的影像重構效果。但是這樣的方法并不利于影像的采集和測量矩陣的設計。在遙感等海量數據采集環境中多是采用線陣推掃數據采集方法。單像素相機采用的最小全變分法（Total Variation，TV）計算復雜度大，難以應用于海量數據處理環境。因此影像重構實驗采用OMP算法。

圖5只是驗證了理想情況下重構矩陣、優化矩陣和近似矩陣的性能。為進一步驗證真實情況下算法1、式（2）和式（3）的可用性和有效性，本文采用Lena圖逐列作線陣推掃式影像采集和重構實驗［11］。

近似矩陣、重構矩陣和優化矩陣的影像采集和重構過程如下。

近似矩陣ΠT：首先，以測量矩陣Φ（0－1隨機矩陣）通過y＝Φx對Lena圖逐列掃描采集；其次，對測量數據y和重構矩陣ΠR事后優化：yT＝Τy，ΠT＝ΤΠR；然后，以OMP算法逐列求解min‖α‖0s.t.yT＝ΠTα，得到；最后，解算＝即可得到Lena圖的重構影像。

重構矩陣ΠR：首先，以測量矩陣Φ（0－1隨機矩陣）通過y＝Φx對Lena圖逐列掃描采集；然后，以OMP算法逐列求解min‖α‖0s.t.y＝ΠRα，得到；最后，解算＝。

優化矩陣ΠO：首先，對影像逐列作稀疏變換α＝Ψx，取得α；其次，通過y＝ΠOα逐列取得Lena圖的測量數據y；然后，以OMP算法逐列求解min‖α‖0s.t.y ＝ ΠOα，得到；最后，解算＝（優化矩陣ΠO無法應用于真實的影像采集作業，只是為了與重構矩陣和近似矩陣作影像重構能力的對比分析，而虛擬了優化矩陣的影像采集）。

3類矩陣（重構矩陣、優化矩陣和近似矩陣）的重構效果如圖6所示（從上到下分別是Lena、Boats和意大利撒丁島Mulargia湖圖；從左至右分別是重構矩陣、優化矩陣和近似矩陣的重構影像，最右邊的是原始影像）。Mulargia湖的遙感影像由Landsat－5衛星在波段4拍攝于1996年7月。表2是3類矩陣重構圖的信噪比（Signal to Noise Ratio，SNR）和峰值信噪比（Peak Signal to Noise Ratio，PSNR）。

由表2可見，重構矩陣的SNR和PSNR最小；優化矩陣的SNR和PSNR比重構矩陣分別約提高0.5～1.5dB；近似矩陣的SNR和PSNR與優化矩陣基本相當，變化不大。這進一步說明算法1、式（2）和式（3）都是有效可行的。OMP算法并非是一個高效的算法，但是它是其他各類匹配追蹤算法的基礎，如采用其他匹配追蹤算法有可能進一步改善影像的重構效果。因此，工業化應用中可以采用基于近似矩陣和線陣推掃模式的單像素相機影像采集和重構方案。

表2 3類矩陣重構Lena、Boats和Mulargia湖圖的信噪比和峰值信噪比Tab.2 SNR and PSNR of reconstructed Lena、Boats and Mulargia lake of 3matrices

圖6 3類矩陣的重構效果圖Fig.6 Reconstruction effect drawing of 3matrices

6 結論

在稀疏變換基（DCT矩陣）和單像素相機測量矩陣（0－1隨機矩陣）已確定的情況下，本文提出了以事前確定的測量矩陣（0－1隨機矩陣）采集數據，以事后優化的重構矩陣（近似矩陣）重構信號的方法。既初步解決了當前測量矩陣難以硬件實現和設計的難題，又改善提高了信號的重構效果。同時回避了“測量矩陣和稀疏變換基不相關性好”的測量矩陣設計原則。實驗和理論分析表明無需把二維影像轉化成一維信號的基于線陣推掃和重構矩陣優化的單像素相機影像采集和重構方法是可行的。擬將本文的測量矩陣轉變為0－1稀疏循環矩陣作進一步的深入研究。

［1］程濤，朱國賓，劉玉安.基于0－1稀疏循環矩陣的測量矩陣分離研究 ［J］.光學學報，2013，33（2）：1－6.

［2］DONOHO D L.Compressed sensing ［J］.Information Theory，IEEE Transactions on，2006，52（4）：289－306.

［3］ELAD M.Optimized Projections for Compressed Sensing［J］.Signal Processing，IEEE Transactions on，2007，55（12）：695－702.

［4］DUARTE－CARVAJALINO J M，SAPIRO G.Learning to Sense Sparse Signals：Simultaneous Sensing Matrix and Sparsifying Dictionary Optimization［J］.Image Processing，IEEE Transactions on，2009，18（7）：395－408.

［5］DONOHO D L，TSAIG Y，DRORI I，et al.Sparse Solution of Underdetermined Systems of Linear Equations by Stagewise Orthogonal Matching Pursuit［J］.Information Theory，IEEE Transactions on，2012，58（2）：94－121.

［6］CANDES E J，TAO T.Near－Optimal Signal Recovery From Random Projections：Universal Encoding Strategies［J］.Information Theory，IEEE Transactions on，2006，52（12）：6－25.

［7］潘榕，劉昱，侯正信，等.基于局部DCT系數的圖像壓縮感知編碼與重構 ［J］.自動化學報，2011，37（6）：74－81.

［8］DUARTE M F，DAVENPORT M A，TAKHAR D，et al.Single－Pixel Imaging via Compressive Sampling ［J］.Signal Processing Magazine，IEEE，2008，25（2）：83－91.

［9］AKCAKAYA M，JINSOO P，TAROKH V.A Coding Theory Approach to Noisy Compressive Sensing Using Low Density Frames ［J］.Signal Processing，IEEE Transactions on，2011，59（11）：69－79.

［10］劉亞新，趙瑞珍，胡紹海，等.用于壓縮感知信號重建的正則化自適應匹配追蹤算法 ［J］.電子與信息學報，2010，32（11）：3－7.

［11］程濤，朱國賓，劉玉安.基于二維壓縮感知的定向遙感和變化檢測研究 ［J］.紅外與毫米波學報，2013，32（5）：1－6.