張建峰 沈 軍 張昊平

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

0 引 言

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

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

1 相關(guān)工作

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

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

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

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

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

(1)

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

(2)

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

g=Hf+ω

(3)

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

(4)

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

(5)

f=Ψθ

(6)

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

1.2 分裂Bregman迭代算法

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

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

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

(7)

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

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

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

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

t←t+1

直到滿足停止迭代條件。

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

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

(8)

(9)

(10)

x=Dα

(11)

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

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

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

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

(12)

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

(13)

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

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)解。首先對(duì)x求偏導(dǎo)得到：

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

(17)

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

(18)

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

(19)

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

(20)

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

3.2 字典更新

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

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

2) 重復(fù)執(zhí)行：

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

(3)k=k+1;

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

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

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

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

s.t.‖α‖0≤L

(21)

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

3.4 算法流程

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

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

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

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

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

(2) 利用改進(jìn)的K-SVD算法訓(xùn)練學(xué)習(xí)得到與原始圖像相適應(yīng)的冗余字典D(t+1)；

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

(4)t←t+1；

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

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

4 實(shí) 驗(yàn)

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

圖2 原始光譜圖像

圖3 光譜圖像壓縮測(cè)量值

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

(22)

(23)

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

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

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

圖4展示了光譜圖像壓縮測(cè)量值在加入不同程度的噪聲下，光譜圖像的重構(gòu)效果。為了便于比較，這里對(duì)每個(gè)波段圖像的峰值信噪比進(jìn)行統(tǒng)計(jì)平均值然后比較?？梢钥闯觯帽疚姆椒ㄖ貥?gòu)的光譜圖像各個(gè)波段的峰值信噪比平均值總體在其他三種方法之上，但是在噪聲影響較大的情況下圖像重構(gòu)質(zhì)量提高較小，而當(dāng)噪聲影響較小時(shí)，圖像重構(gòu)質(zhì)量提高較大。

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

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

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

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

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

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

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

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

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

(24)

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

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

5 結(jié) 語

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