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

張建峰 沈 軍 張昊平

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

0 引 言

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

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

1 相關(guān)工作

1.1 編碼孔徑光譜成像系統(tǒng)數(shù)學模型

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

圖1 CASSI成像系統(tǒng)抽象模型

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

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

(1)

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

(2)

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

g=Hf+ω

(3)

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

(4)

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

(5)

f=Ψθ

(6)

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

1.2 分裂Bregman迭代算法

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

考慮一個約束優(yōu)化問題如下：

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

(7)

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

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

循環(huán)執(zhí)行：

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

t←t+1

直到滿足停止迭代條件。

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

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

(8)

(9)

(10)

x=Dα

(11)

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

3 自適應分裂Bregman迭代的圖像重構(gòu)算法

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

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

(12)

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

(13)

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

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

(14)

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

(15)

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

(16)

3.1 x求解問題

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

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

(17)

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

(18)

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

(19)

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

(20)

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

3.2 字典更新

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

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

2) 重復執(zhí)行：

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

(3)k=k+1;

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

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

3.3 稀疏表示系數(shù)α求解問題

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

s.t.‖α‖0≤L

(21)

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

3.4 算法流程

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

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

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

3) 循環(huán)執(zhí)行下列步驟：

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

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

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

(4)t←t+1；

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

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

4 實 驗

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

圖2 原始光譜圖像

圖3 光譜圖像壓縮測量值

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

(22)

(23)

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

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

表1 重構(gòu)各個波段圖像的峰值信噪比 dB

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

圖4 不同程度噪聲影響下重構(gòu)圖像質(zhì)量

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

(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

利用結(jié)構(gòu)相似性(SSIM)指標進一步衡量圖像質(zhì)量，計算公式為：

(24)

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

表2 重構(gòu)光譜圖像的MSSIM值

5 結(jié) 語

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