自適應分裂Bregman迭代的編碼孔徑光譜圖像重構方法

張建峰 沈 軍 張昊平

1(哈爾濱工程大學信息與通信工程學院 黑龍江 哈爾濱 150001)2(上海無線電設備研究所 上海 200090)

0 引 言

光譜成像技術通過獲取目標的二維空間信息和一維光譜信息來組成圖譜結合的數據立方體[1]，它能夠幫助我們識別場景中出現的不同物體，具有廣闊的應用前景。傳統的光譜成像技術，在探測器一次測量中只能獲得數據立方體的一個切面，需要反復測量才能獲得目標圖像，采樣過程不僅花費大量時間，而且測量值占大量存儲空間，越來越不能滿足現實測量的要求。編碼孔徑光譜成像技術[2-5]是結合壓縮感知理論的新型光譜成像方法，很好地打破了傳統光譜成像技術的局限性，該系統利用二維探測器捕捉的是一個光譜圖像數據立方體的多通道投影。因此，利用編碼孔徑光譜成像系統在很大程度上減少了采樣數據，提高了圖像采集效率。另一方面，編碼孔徑光譜成像系統采集到的圖像是高度壓縮的[6]，因此從二維測量值中準確高質量重構出三維的光譜圖像立方體成為了一個重要的挑戰。為了從壓縮測量值中更好地重構出光譜圖像，人們做了大量的研究。廣泛應用的是將該重構問題轉化為l1范數最小化問題[3,7-8]，利用相應的重構算法進行求解。Arguello等[7]利用梯度稀疏投影稀疏重建(GPSR)[9]算法，它使用二維小波變換矩陣與一維余弦變換基的克羅內克積構成稀疏基代入重構框架進行求解。Wagadarikar等[3]則利用兩步迭代收縮閾值(TwIST)框架進行求解，它主要將l1范數替換成全變分模型作為正則化項[10]，這種方法能快速地恢復出原光譜圖像。文獻[11]則采用一種近似消息傳遞框架(AMP)，引入自適應維納濾波器重構光譜圖像。除了像文獻[7]中使用三維解析字典外，采用字典學習的方式提高光譜圖像重構質量也成為一種有效的方式。文獻[12]分別將解析字典和K-SVD離線訓練字典代入重構框架，結果顯示采用K-SVD離線訓練的字典重構的光譜圖像質量得到有效提高。文獻[13]則在壓縮投影和光譜圖像間建立神經網絡，通過深度學習的方式重構光譜圖像。但是構建大量與目標光譜圖像相匹配的訓練集十分困難，并且學習過程會耗費大量時間。

上述算法雖然在一定程度上實現光譜圖像重建，但都沒能充分利用光譜圖像在特定冗余字典下的稀疏特性。壓縮感知理論指出，信號在冗余字典下表示越稀疏，則重構的質量越好[14-15]。可見構造合適的冗余字典對于圖像重構至關重要，因此本文提出了一種自適應分裂Bregman迭代的編碼孔徑光譜圖像重構方法。該算法以分塊測量方式構建的圖像重構框架為基礎[16-17]，很好地將分裂Bregman迭代算法[18]與字典學習算法相結合。在迭代的過程中利用前一次圖像的估計值構造訓練樣本，并結合改進的K-SVD算法自適應學習更新字典，此時得到的字典能對光譜圖像實現更有效的稀疏表示，多次迭代之后得到待重構圖像。相比之下文獻[12]則采用其他光譜圖像樣本學習得到字典，其字典原子不一定與待重構圖像相適應。值得注意的是改進的K-SVD算法，其特點是利用K-SVD算法[19]訓練好的字典再次對樣本進行稀疏表示，提取出誤差較大的原子組成新的矩陣并進行一次SVD運算，將變換后前若干奇異值向量加入到原字典當中，進一步使該冗余字典與待重構圖像相適應。仿真實驗將本文算法與廣泛使用的GPSR、TwIST以及OMP重構光譜圖像算法進行了比較，通過視覺比較重構光譜圖像的細節以及利用峰值信噪比(PSNR)和結構相似性(SSIM)指標作出對比，證明了本文算法的有效性。

1 相關工作

1.1 編碼孔徑光譜成像系統數學模型

本文主要以編碼孔徑快照光譜成像系統(CASSI)模型為基礎進行研究。CASSI的成像過程主要包括：首先物鏡將目標光源聚焦到編碼孔徑模板上，實現了空間信息采樣；再通過色散元器件對不同譜段的投影進行偏移實現編碼；然后利用二維的焦平面陣列有探測器測得混疊的圖譜信息；最后根據混疊的信息重構出目標圖像。圖像測量模型如圖1所示。

圖1 CASSI成像系統抽象模型

假設F0和F1都為M×N×L數據立方體，分別表示原始光譜圖像和通過編碼孔徑模板后的圖像。

F1(x,y,λ)=T(x,y)×F0(x,y,λ)

(1)

式中：x、y表示空間坐標，λ表示光譜值;T表示模板，一般用服從伯努利分布的二值矩陣表示，″0″表示不透光，″1″表示透光。原始數據立方體中的所有高光譜片段都使用同樣的空間模板來采樣。采用光色散器件將光譜圖像各個波段的圖像沿色散方向彼此錯開實現編碼，最終在探測器陣列上成像。成像圖案表示如下：

(2)

式中：s(λ)代表光譜沿色散方向的移動函數。式(2)中假設沿x軸色散。以參考波長λ為中心，根據棱鏡色散系數決定不同波長光譜側向移動空間距離大小，參考波長以0度入射時假定無偏移產生。假設目標場景的空間維度為M×N，同時光譜維度為L，則總的像素個數為M×N×L。假設比鄰的光譜間側向色散距離為1，在探測器陣列接收過程中光強沿光譜方向疊加，最后得到一張測量圖像，其維度為M×(N+L-1)。為了便于計算，將整個過程表示為：

g=Hf+ω

(3)

式中：g是G的向量形式，長度為M×(N+L-1)；f=vec([f0,f1,…,fL-1]) 是F的向量表示，fk為第k個波段的圖像數據；w表示測量過程中的噪聲；H為觀測矩陣，對于CASSI來說，H是由編碼孔徑形式和色散元器件共同決定的[16]，表示整個采樣和色散過程。

(4)

假設f在冗余字典Ψ下的表示是稀疏的，則可以采用壓縮感知計算方法根據測量值重構出圖譜信息，表示為：

(5)

f=Ψθ

(6)

式中:θ為某冗余字典Ψ下的表示系數;K為稀疏性約束因子。

1.2 分裂Bregman迭代算法

為了方便后續部分的討論，在這里概要介紹經典的分裂Bregman算法，它是由Goldstein和Osher提出的用于解決l1范數正則化優化問題。文獻[18]證明了該算法的收斂性，在迭代的過程中不斷逼近最優解。該算法在提高圖像處理質量的同時保持高效的計算效率，在圖像處理領域應用廣泛[20]。

考慮一個約束優化問題如下：

minu∈RN,ν∈RMf(u)+g(ν) s.t.u=Gν

(7)

式中：G∈RN×M，并且f:RN→R,g:RN→R是凸函數。SBI算法對該問題求解的具體步驟如下：

初始化：設置t、b0、u0、ν0的初始值，參數μ>0。

循環執行：

b(t+1)←b(t)-(u(t+1)-Gν(t+1))

t←t+1

直到滿足停止迭代條件。

2 分塊測量編碼孔徑光譜成像數學模型

本文主要采用重疊分塊測量的方式構造數學模型。分塊的好處是可以很好地降低計算復雜度[16-17]，為之后構造訓練樣本集提供前提條件。分塊即在空間維進行分塊測量，光譜維保持不變，進行獨立的測量和恢復，最后合成整個目標圖像。其過程可以表達為：

(8)

(9)

(10)

x=Dα

(11)

式中：D為稀疏基或冗余字典；α為x在D下的稀疏表示系數。測量矩陣H的構造方式不變，只是矩陣的規模變小。

3 自適應分裂Bregman迭代的圖像重構算法

在這一部分主要介紹了對原有數學模型進行變換，利用分裂Bregman迭代求解出相應的圖像塊；以及每次迭代過程中將求解的圖像塊構造訓練樣本集，結合改進的K-SVD算法自適應學習更新字典；多次迭代后逐漸逼近求解出待重構圖像。

考慮到式(10)中采用的l0范數是非凸函數且是NP難問題，一般采用l1范數來最佳凸逼近l0范數。則上述模型可以改為如下優化問題：

(12)

通過借鑒分裂Bregman算法將式(12)分裂成兩個子問題進行最小化求解，這樣做的好處是對于子問題的求解難度要遠低于對原問題求解，并且經過反復迭代逐漸逼近最優解。

(13)

利用x(t+1)構造訓練樣本集r(t+1)=x(t+1)-b(t)，利用改進K-SVD算法進行字典更新，即對下式進行求解：

s.t.‖αk‖≤l,?k

(14)

在更新后的字典下求解稀疏表示系數為：

(15)

b(t+1)=b(t)-(x(t+1)-D(t+1)α(t+1))

(16)

3.1 x求解問題

式(13)可采用最小二乘法得出最優解。首先對x求偏導得到：

HT(y-Hx)+μ(x-D(t)α(t)-b(t))=0

(17)

進而估計出x(t+1)如下：

(18)

在每次迭代過程中字典D和系數α都會進行更新，使得每次估計值不斷接近目標圖像，最終求解出目標圖像。考慮到HTH不具有對角結構，故求逆(HTH+μI)求逆過程計算代價很大，在每次循環的過程中引入泰勒公式二次項變換為：

(19)

此時將式(19)代入式(17)同樣采用最小二乘法解得：

(20)

式中：I為單位矩陣。這樣求逆的部分變成對角陣，降低了計算的復雜性。

3.2 字典更新

在成像的過程中冗余字典D的構造至關重要，它直接影響到圖像稀疏表示，進而影響圖像重構質量。本文在分裂Bregman迭代求解的基礎上采用一種改進的K-SVD算法進行自適應訓練更新冗余字典。主要是利用迭代過程中待重構圖像的近似值構造樣本訓練集結合K-SVD字典學習算法更新字典；之后再次對樣本逐列進行稀疏表示，提取誤差較大的原子組成新的矩陣，進行一次SVD運算，將前若干列左奇異值向量加入到字典當中。這樣做的好處是豐富了字典原子的特征，能夠對原始圖像進行更好的稀疏表示。具體實現步驟如下：

1) 初始化：設置閾值ε，訓練樣本集r,空矩陣w,k=0，利用K-SVD算法得到字典D。

2) 重復執行：

(1) 對rk在字典D下進行稀疏表示，得到系數αk；

(3)k=k+1;

(4) 直到遍歷整個樣本集。

3) 對w進行一次SVD變換，提取前若干列左奇異值向量加入到字典D當中。

3.3 稀疏表示系數α求解問題

對于稀疏系數的求解問題就是對式(14)進行求解，該問題可以等價于求解為下列問題：

s.t.‖α‖0≤L

(21)

對α的求解就變成了常規的稀疏編碼問題，該式可以通過正交匹配追蹤算法(OMP)進行求解，其中L表示稀疏系數α的稀疏度。

3.4 算法流程

根據以上分析，給出具體的算法流程如下：

1) 輸入：分塊取樣算子R、測量矩陣H、μ、δ、ε、tol以及分塊測量值構成的矩陣y。

2) 初始化：冗余字典D0、x(0)、b(0),本文中D0為隨機抽取同類型其他光譜圖像塊組成的字典；

3) 循環執行下列步驟：

(1) 利用式(9)解出x(t+1)、D(t)、α(t)分別表示上一次更新得到的字典以及相應的系數,x(t)為上一次估計的圖像值,y表示原始測量值；

(2) 利用改進的K-SVD算法訓練學習得到與原始圖像相適應的冗余字典D(t+1)；

(3) 利用OMP算法求解式(21),利用更新后的字典再求獲得稀疏表示系數α(t+1)；

(4)t←t+1；

(4) 滿足循環終止條件,輸出x；

(5) 根據式(9)解出目標圖像。

4 實 驗

實驗通過仿真CASSI系統，以及利用本文方法重構光譜圖像，比較其成像效果。實驗中采用的原始光譜圖像來自于哥倫比亞大學光譜圖像數據庫和ICVL數據集，如圖2所示截取圖像大小為256×256×8，包括8個波段，波長范圍為440～650 nm。通過經典的CASSI系統仿真得到的測量圖如圖3所示。光譜圖像重構即由測量值恢復出原始各個波段的圖像。

圖2 原始光譜圖像

圖3 光譜圖像壓縮測量值

本文采用分塊采樣的方式，使用6×6的窗口進行重疊取塊，移動步長為4，分塊可以有效降低計算復雜度。仿真得到測量值矩陣，然后代入式(12)求解框架中，利用提出的自適應分裂Bregman算法進行求解，重構出目標光譜圖像。該框架中的測量矩陣H是由采樣模板T以式(4)構成，而仿真中采樣模板T主要為服從伯努利分布的隨機矩陣，仿照大量文獻中分布參數設置為0.5。實驗中利用本文算法進行光譜圖像重構，設置初始參數μ=τ=0.05，ε=tol=103。最后與文獻[7]采用GPSR算法、文獻[3]采用的TwIST算法以及壓縮感知中常用的OMP算法重構的光譜圖像進行比較。比較的方法從主觀和客觀上考慮，主觀評價是圖像的成像效果，客觀評價首先以常用的峰值信噪比來進行衡量。考慮到光譜圖像的特殊性，在計算時分別求出每個波段對應的PSNR值，計算如下式：

(22)

(23)

式中:λ對應第λ個波段的圖像，I和K分別表示重構后的圖像以及原始圖像，MAXIλ表示第λ個波段上圖像的最大值。

表1展示了重構的各個波段圖像的峰值信噪比。可以看出，利用本文方法重構的光譜圖像峰值信噪比明顯高于其他算法，而GPSR算法和OMP算法總體上好于TwIST算法。其中：本文算法與GPSR算法和OMP算法求解框架都是以稀疏系數的l1范數作為正則化條件，引入稀疏基或冗余字典；而TwIST是以圖像空間信息平滑性作為正則化條件進行約束。可以看出通過利用光譜圖像在某一冗余字典下的稀疏性作為約束更加符合解決編碼孔徑光譜圖像重構的問題。本文在分裂Bregman算法迭代求解過程中利用改進的K-SVD自適應訓練更新字典，使得冗余字典能夠對待重構圖像進行更好的稀疏表示，充分利用光譜圖像在特定冗余字典下具有稀疏性這一特征，進而提高圖像的質量。

表1 重構各個波段圖像的峰值信噪比 dB

圖4展示了光譜圖像壓縮測量值在加入不同程度的噪聲下，光譜圖像的重構效果。為了便于比較，這里對每個波段圖像的峰值信噪比進行統計平均值然后比較。可以看出，利用本文方法重構的光譜圖像各個波段的峰值信噪比平均值總體在其他三種方法之上，但是在噪聲影響較大的情況下圖像重構質量提高較小，而當噪聲影響較小時，圖像重構質量提高較大。

圖4 不同程度噪聲影響下重構圖像質量

圖5-圖7是利用三組光譜圖像進行測試的結果，分別選取其中一個波段的圖像并截取一部分進行顯示，通過觀察圖5刻度尺部分、圖6建筑物邊緣以及圖7的路燈部分清晰可見，利用OMP和GPSR重構出的光譜圖像在邊緣部分不平滑，出現很多相互交錯的噪聲，總體效果要好于TwIST算法恢復出的圖像；相比之下TwIST恢復出的圖像雖然很平滑但是圖像細節部分丟失嚴重，整體較為模糊；而利用本文方法恢復出的圖像細節更加清晰，而且圖像邊緣部分也恢復得很好。

(a) OMP算法 (b) GPSR算法

(c) TwIST算法 (d) 本文方法圖5 測試圖像1

(a) OMP算法 (b) GPSR算法

(c) TwIST算法 (d) 本文方法圖6 測試圖像2

(a) OMP算法 (b) GPSR算法

(a) TwIST算法 (b) 本文方法圖7 測試圖像3

利用結構相似性(SSIM)指標進一步衡量圖像質量，計算公式為：

(24)

式中：μ表示圖像的均值；σ表示圖像方差或協方差；c1和c2表示常數。這里的圖像均指某一特定波段上的圖像。同樣為了方便比較，采用平均結構相似性(MSSIM)，即對每個波段圖像的SSIM求統計平均。表2展示了上述三幅測試圖像重構以后的MSSIM值，MSSIM值越大則表示重構圖像與原圖像越接近，可以看出本文方法與其他三種算法相比具有一定的優勢。

表2 重構光譜圖像的MSSIM值

5 結 語

為了解決編碼孔徑光譜成像技術中，利用傳統方法造成光譜圖像重構質量不高等問題，本文提出了一種自適應分裂Bregman迭代的圖像重構算法。該算法很好地將分裂Bregman迭代算法與字典學習相結合，在每次迭代的過程中自適應訓練更新字典。充分利用光譜圖像在特定冗余字典下具有稀疏性這一特征，提高圖像重構質量。實驗通過與傳統光譜圖像重構過程中廣泛使用的GPSR、OMP以及TwIST算法進行比較，結果表明本文算法性能更好。同時，與上述三種方法相比本文算法也存在重構時間長的問題，下一步將進一步改善字典學習的過程，縮短重構時間。