楊虹 黃遠(yuǎn)輝 鐘宬 邵曉鵬

(西安電子科技大學(xué)技術(shù)物理學(xué)院，光電成像與圖像處理實(shí)驗(yàn)室，西安 710071)

1 引言

隨著遙感技術(shù)的發(fā)展，出現(xiàn)了多種光譜成像儀，例如擺掃成像儀、推掃成像儀和可調(diào)諧濾波成像儀[1-3]。這些成像儀的設(shè)計(jì)理念相對(duì)簡(jiǎn)單，但由于大都僅獲得數(shù)據(jù)立方體的一維或二維子集，因而需要對(duì)其它維度進(jìn)行實(shí)時(shí)掃描以此獲取完整的數(shù)據(jù)立方體; 再者，對(duì)于非相干光源，它們的光子收集率都較低，導(dǎo)致信噪比也偏低。一些最新發(fā)展起來(lái)的光譜成像儀，雖然光子收集率的問(wèn)題在很大程度上得到解決，但又引入了諸如“失蹤三角錐”等一些新問(wèn)題[4]。

上述這些傳統(tǒng)的光譜成像儀均存在一個(gè)很重要的特點(diǎn)，那就是獲得的觀測(cè)值的總數(shù)大于或等于重建的數(shù)據(jù)立方體總元素個(gè)數(shù)。與此不同，Brady等人2006年首次提出了壓縮編碼孔徑光譜成像的概念，這種成像方法基于壓縮感知理論，得到的觀測(cè)值總數(shù)遠(yuǎn)小于數(shù)據(jù)立方體的元素總數(shù)[5]。

壓縮編碼孔徑成像是近年來(lái)基于壓縮感知理論新興的一種成像技術(shù)，它在可見光成像方面的應(yīng)用研究相對(duì)較多，但在多光譜成像領(lǐng)域的研究相對(duì)較少。在保證圖像品質(zhì)的前提下由于壓縮編碼孔徑成像技術(shù)能大大減少采樣的數(shù)據(jù)量，使傳統(tǒng)遙感光譜成像面臨的最普遍的數(shù)據(jù)量龐大的問(wèn)題得到解決，所以本文提出將壓縮編碼孔徑成像應(yīng)用于遙感成像。然而，現(xiàn)有的壓縮編碼孔徑光譜成像系統(tǒng)僅能用于靜態(tài)成像，而遙感成像是個(gè)動(dòng)態(tài)成像的過(guò)程，因 此要想將壓縮編碼孔徑應(yīng)用于遙感光譜成像中，必須分析衛(wèi)星平臺(tái)運(yùn)動(dòng)對(duì)成像過(guò)程的影響，進(jìn)而克服由于衛(wèi)星平臺(tái)運(yùn)動(dòng)導(dǎo)致的圖像品質(zhì)的下降。

2 壓縮編碼孔徑光譜成像系統(tǒng)原理

2.1 壓縮編碼孔徑原理簡(jiǎn)介

壓縮編碼孔徑的靈感很大程度上來(lái)源于壓縮感知這個(gè)新興的領(lǐng)域[6-8]。壓縮感知也被稱為壓縮采樣或稀疏采樣，是一種尋找欠定線性系統(tǒng)的稀疏解，利用稀疏的或可壓縮的信號(hào)進(jìn)行信號(hào)重建的方法，是近年來(lái)由D. Donoho，E. Candès，J. Romberg和T. Tao等人提出的一種全新的采樣理論[9-11]。該理論指出: 只要信號(hào)是可壓縮的或在某個(gè)域上是稀疏的，就可以用一個(gè)與變換基不相關(guān)的測(cè)量矩陣將變換所得高維信號(hào)投影到一個(gè)低維空間上，使實(shí)際采樣率遠(yuǎn)低于 Nyquist采樣率，然后通過(guò)求解一個(gè)優(yōu)化問(wèn)題從這些少量的投影中以高概率重構(gòu)出原信號(hào)，這也證明了這樣的投影包含了重構(gòu)信號(hào)的足夠信息[11]。

編碼孔徑成像過(guò)程分為兩步: 第一步是編碼過(guò)程，光線透過(guò)編碼板每個(gè)小孔后在探測(cè)器上都形成投影圖像，這些投影圖像疊加在一起，在二維平面上形成退化的二維分布信號(hào)，即編碼圖像; 第二步是解碼過(guò)程，對(duì)探測(cè)器采集到的圖像數(shù)據(jù)進(jìn)行濾波和重建，復(fù)原出原始目標(biāo)圖像[12]，其基本過(guò)程如圖1所示:

圖1 編碼孔徑成像過(guò)程Fig.1 Imaging procedure of Compressive Coded Aperture

2.2 壓縮編碼孔徑快照光譜成像系統(tǒng)

2008年，Wagadarikar等人改進(jìn)了他們?cè)?007年提出的單分光元件編碼孔孔徑光譜成像儀[13]，其結(jié)構(gòu)如圖2所示。這是目前世界上發(fā)展相對(duì)較成熟的一套編碼孔徑光譜成像系統(tǒng)，由物鏡、編碼孔徑、帶通濾波器、F/8中繼透鏡、雙阿米西棱鏡和單色CCD探測(cè)器構(gòu)成。物鏡將場(chǎng)景成像于編碼孔徑平面，經(jīng)編碼孔徑編碼再經(jīng)過(guò)帶通濾波器，中繼鏡頭和雙阿米西棱鏡到達(dá)探測(cè)器平面。

圖2 單分光元件編碼孔徑快照光譜成像儀光路Fig.2 Optical path of Single Disperser Coded Aperture Snapshot Spectral Imager (SD-CASSI)

3 遙感衛(wèi)星平臺(tái)運(yùn)動(dòng)建模

隨著科技的進(jìn)步，人們對(duì)于光譜圖像品質(zhì)的要求越來(lái)越高，消除衛(wèi)星平臺(tái)運(yùn)動(dòng)引起的光譜混疊對(duì)提高光譜圖像品質(zhì)起著十分重要的作用。

3.1 衛(wèi)星平臺(tái)運(yùn)動(dòng)姿態(tài)及其對(duì)應(yīng)像移的計(jì)算

衛(wèi)星平臺(tái)的運(yùn)動(dòng)在空間上通常有3種運(yùn)動(dòng)姿態(tài): 俯仰、側(cè)滾和偏航。唐秋艷等人對(duì)衛(wèi)星運(yùn)動(dòng)姿態(tài)的改變對(duì)圖像造成的影響進(jìn)行研究[14]，并且引入了平均摻雜比的概念，架起了衛(wèi)星平臺(tái)參數(shù)與光譜成像之間的橋梁。在該項(xiàng)研究中，若將探測(cè)器陣列上一像素(m,n)處的光譜密度記為Bm(v)，其它像元的混疊光譜統(tǒng)一記作Bm+1(v)，混疊系數(shù)為iη，則(m,n)處的光譜密度可記為

由于平臺(tái)的穩(wěn)定系數(shù)一般是10-3(°)/s ～10-1(°)/s，即平臺(tái)運(yùn)動(dòng)引起的混疊通常都不會(huì)超過(guò)一個(gè)像元間隔，所以可以將Bm+1(v)直接視為相鄰像元的光譜密度，于是式(1)可改寫為

由光譜混疊模型中混疊光譜各成分的比例等于各地物面積之比可知[14-15]，平均混疊權(quán)重系數(shù)由像移路徑包圍的總面積決定。

俯仰會(huì)在探測(cè)器的x方向產(chǎn)生像移，像移為Δx(t)=ftan[Δα(t)]; 側(cè)滾在探測(cè)器的y方向產(chǎn)生像移，由圖可得Δy(t)=ftan[Δβ(t)]; 偏航在水平和垂直方向均產(chǎn)生像移，但由于垂直方向像移較小，可近似為0，其像移主要集中在水平方向，且位移大小與偏航角度相關(guān)Δx(t)=y·Δθ(t)。其中,f為系統(tǒng)焦距; Δα(t)，Δβ(t)為t時(shí)刻的振動(dòng)角度，由于衛(wèi)星平臺(tái)的穩(wěn)定度較高，俯仰和側(cè)滾的像移可分別近似表示成Δx(t)=f·Δα(t)，Δy(t)=f·Δβ(t)。

3.2 平均混疊權(quán)重系數(shù)矩陣的生成

在實(shí)際成像過(guò)程中，由于運(yùn)動(dòng)造成的像移往往是上述3種像移的疊加且衛(wèi)星的穩(wěn)定性導(dǎo)致像移不會(huì)超出一個(gè)像素范圍，因此混疊光譜可簡(jiǎn)化為原始光譜與其周圍8個(gè)相鄰像素光譜的疊加[15]，如圖3所示:

圖3 8鄰域像素光譜混疊示意Fig.3 Spectrum aliasing from the 8 neighbor pixels of the center pixel

觀察可得，探測(cè)器上的像移有4種可能的情況，可描述為相鄰8個(gè)光譜中的3個(gè)參與混疊，即光譜2,5,9參與混疊，光譜4,5,8參與混疊，光譜3,4,7參與混疊，光譜2,3,6參與混疊。

假設(shè)衛(wèi)星平臺(tái)的運(yùn)動(dòng)是一個(gè)隨機(jī)過(guò)程，即在一個(gè)運(yùn)動(dòng)周期內(nèi)4種光譜混疊情況等概率分布，且每種混疊情況的平均像移取最大像移量與最小相移量的平均值，設(shè)混疊光譜像元2～9所占平均混疊權(quán)重系數(shù)分別對(duì)應(yīng)，由圖3可得:

其中:

式(2)可改寫為:

其中，B(m,n)(v)表示不存在光譜混疊時(shí)的核心像素處的光譜信息;B(m,n)′(v)表示存在光譜混疊后的核心像素處的光譜信息:η(p,q)，B(p,q)(v)分別表示核心像素8個(gè)鄰域像素對(duì)應(yīng)的混疊權(quán)重系數(shù)和光譜信息，將式(8)代入式(7)，即可獲得每個(gè)像素存在混疊后的光譜。

4 仿真算法與結(jié)果分析

4.1 結(jié)合運(yùn)動(dòng)模糊模型的TwIST重建算法

在重建過(guò)程中，采用了Bioucas-Dias and Figueiredo提出的兩步迭代收縮閾值算法(Two-step Iterative Shrinkage/Thresholding , TwIST)[16]。此算法結(jié)合了迭代重權(quán)值收縮算法(Iterative Re-weighted Shrinkage ,IRS)和迭代收縮閾值算法(Iterative Shrinkage/Thresholding , IST)的優(yōu)點(diǎn)，既能解決方程嚴(yán)重欠定的問(wèn)題又具有較好的去噪效果。

TwIST算法解的形式為:

其中:

Γλ表示算子操作;α，β表示權(quán)重系數(shù);K表示編碼矩陣。

TwIST算法的流程大致可分為以下4部分:

1) 獲取參數(shù)β，α，和τ;

2) 賦初始值x0;

3) 迭代計(jì)算:

4) 當(dāng)估計(jì)值與原圖的差值小于某一閾值時(shí)，結(jié)束迭代，得到的x即為所求。

在本文中，此算法將求數(shù)據(jù)立方體的解等效為求如下這個(gè)非線性無(wú)約束最小化問(wèn)題的最優(yōu)解的過(guò)程:

式中g(shù)表示理想圖像;f表示真實(shí)圖像;H表示模糊矩陣;Φ(f)表示正則項(xiàng);τ表示正則系數(shù)。

其中，采用全變分的方法進(jìn)行正則化，即所謂的TV范數(shù):

對(duì)于連續(xù)和離散的情況，TV范數(shù)分別有如下形式:

4.2 基于CASSI系統(tǒng)的仿真試驗(yàn)設(shè)計(jì)

整個(gè)試驗(yàn)采用MATLAB進(jìn)行仿真和重建，其流程大致如下: 首先用編碼孔徑模板卷積選取多光譜原圖，然后再用運(yùn)動(dòng)模糊算子卷積編碼后的圖像，以此來(lái)模擬實(shí)際遙感編碼孔徑的成像過(guò)程; 在獲得編碼和加入模糊的圖像后，采用TwIST算法對(duì)圖像進(jìn)行重建，并對(duì)重建結(jié)果進(jìn)行了分析和總結(jié)。

在仿真過(guò)程中，采用了由Brady等人用單分光元件編碼孔徑快照光譜成像儀所拍攝的光譜數(shù)為 33的多光譜圖像作為源圖像，它對(duì)應(yīng)的成像場(chǎng)景的彩色圖像如圖4[13]。

在衛(wèi)星平臺(tái)參數(shù)方面，采用了美國(guó)宇航局戈達(dá)德航天中心提供的LAND-SAT-4衛(wèi)星在軌測(cè)量的振動(dòng)數(shù)據(jù)。衛(wèi)星平臺(tái)的振動(dòng)主要來(lái)源于太陽(yáng)能電池陣列驅(qū)動(dòng)所產(chǎn)生的頻率為1Hz、振幅為100μrad的振動(dòng)，以及衛(wèi)星反作用輪基波的二次諧波所產(chǎn)生的頻率為100Hz、振幅為4μrad和頻率為200Hz、振幅為0.6μrad的振動(dòng)，衛(wèi)星振動(dòng)可認(rèn)為是這3個(gè)諧波振動(dòng)的合成[15]:

圖4 成像場(chǎng)景的可見光圖像Fig.4 Visible image of the imaging scene

其中φ1,φ2,φ3是隨機(jī)生成的初始相位角，取值范圍為[0, 2π]。在實(shí)際仿真中將這3個(gè)初始相位角取為[φ1,φ2,φ3]=[0, π/6, π/3]，其周期遠(yuǎn)大于光譜儀的曝光時(shí)間。

在試驗(yàn)中，假設(shè)在一個(gè)振動(dòng)周期內(nèi)，平臺(tái)抖動(dòng)引起的光譜混疊等概率分布于4種混疊方式中，而每一塊像元的混疊面積取4種混疊情況下，該像元在一個(gè)振動(dòng)周期內(nèi)可能出現(xiàn)的最大混疊面積與最小混疊面積的平均值的求和，然后與原始像元面積求比值，即得到該像元對(duì)原始像元的混疊權(quán)重系數(shù)，核心處的原始像元混疊權(quán)重系數(shù)可通過(guò)用1減去其它8個(gè)鄰近像素的混疊權(quán)重系數(shù)的和得到，最終生成混疊權(quán)重系數(shù)矩陣，在焦距與像元邊長(zhǎng)約為10∶1時(shí)，混疊權(quán)重系數(shù)矩陣如圖5和圖6所示。

圖5 俯仰與側(cè)滾綜合作用的混疊權(quán)重系數(shù)矩陣Fig.5 Aliasing coefficient matrix of pitch and roll combined

圖6 偏航產(chǎn)生的混疊權(quán)重系數(shù)矩陣Fig.6 Aliasing coefficient matrix of yaw

比較兩個(gè)權(quán)重系數(shù)矩陣，不難發(fā)現(xiàn)第二個(gè)矩陣元素遠(yuǎn)小于第一個(gè)矩陣，可忽略。因此在實(shí)際操作中，僅采用由俯仰和側(cè)滾產(chǎn)生的混疊權(quán)重系數(shù)矩陣作為圖像的模糊算子，用此矩陣對(duì)探測(cè)器上獲得數(shù)據(jù)立方體進(jìn)行模糊操作。

在此需強(qiáng)調(diào)一點(diǎn)，當(dāng)平臺(tái)運(yùn)動(dòng)引起的像移超出一個(gè)像素范圍時(shí)，光譜混疊的范圍就不再局限于原始像素的8鄰域，而會(huì)向外擴(kuò)展到更大的范圍，此時(shí)需要重新獲得相應(yīng)的更大尺寸的混疊系數(shù)矩陣。

4.3 結(jié)論分析

實(shí)驗(yàn)結(jié)果表明，衛(wèi)星平臺(tái)的運(yùn)動(dòng)會(huì)導(dǎo)致嚴(yán)重的光譜混疊，若不針對(duì)其進(jìn)行去模糊操作，復(fù)原的圖像不能去除光譜混疊，與原圖差距很大; 而在重建中加入模糊算子進(jìn)行去模糊操作后，得到的圖像光譜混疊情況去除良好，分辨率得到很大提升。

圖7 場(chǎng)景在33個(gè)波段上的原始光譜圖像Fig.7 33 original spectral images of the scene

圖8 加入運(yùn)動(dòng)模型后重建出的圖像Fig.8 Recovered images applied the motion model

圖7為圖4場(chǎng)景下的33通道多光譜圖像，圖8為加入運(yùn)動(dòng)模糊后的多光譜圖像。對(duì)比兩幅圖，很明顯地看出與原圖相比，加入運(yùn)動(dòng)模型后的圖像在483nm, 505～537nm, 637nm和650nm波段處圖像強(qiáng)度提升很多，即發(fā)生了嚴(yán)重的光譜混疊。

圖9為在原圖7的基礎(chǔ)上加入運(yùn)動(dòng)模型，并在圖像重建時(shí)引入相應(yīng)模糊算子進(jìn)行去模糊操作后得到的圖像。為了更好地體現(xiàn)去模糊重建的效果，截取了第21波段圖像的局部圖進(jìn)行對(duì)比說(shuō)明。

從圖10中不難看出，加入運(yùn)動(dòng)模糊后的圖像在邊界處(香蕉與蘋果的分界處和蘋果與背景的邊界處)產(chǎn)生了嚴(yán)重的偽輪廓，并且丟失了很多細(xì)節(jié)信息(例如蘋果的柄在(b)中是看不見的)。經(jīng)過(guò)去模糊重建后，成功消除了邊界區(qū)域的偽輪廓，并且恢復(fù)了如蘋果柄這樣的細(xì)節(jié)信息。

圖9 進(jìn)行去模糊重建后的圖像Fig.9 Recovered image with deblured

圖10 第21譜段局部對(duì)比圖Fig.10 Local contrast of the 21st spectrum

圖11為含去模糊操作的重建結(jié)果和不含去模糊操作的重建結(jié)果的PSNR(峰值信噪比)對(duì)比圖。紅線表示去模糊過(guò)的圖像的PSNR值，藍(lán)線表示未去模糊的圖像的PSNR值?？梢姡ツ：膱D像PSNR值普遍遠(yuǎn)高于未去模糊的圖像。

圖12為原圖、含去模糊操作的重建結(jié)果和不含去模糊操作的重建圖像的平均梯度對(duì)比圖?？梢姡ツ：亟ǖ膱D像與原圖更一致，更接近真實(shí)值，未去模糊的圖像則相對(duì)有很大偏離。

圖11 重建結(jié)果PSNR對(duì)比Fig.11 PSNR of recovered images

圖12 重建結(jié)果平均梯度對(duì)比Fig.12 Average gradient of recovered images

通過(guò)以上對(duì)比，可以看出在多光譜遙感中，由于平臺(tái)抖動(dòng)引起的光譜混疊是不可忽視的，而重建過(guò)程中是否包含相應(yīng)的去模糊操作對(duì)復(fù)原圖像的品質(zhì)有至關(guān)重要的影響，因此對(duì)平臺(tái)抖動(dòng)進(jìn)行建模，再根據(jù)抖動(dòng)模型得到相應(yīng)模糊算子進(jìn)行去模糊重建對(duì)提高多光譜圖像品質(zhì)是十分有意義的。

5 結(jié)束語(yǔ)

本文通過(guò)對(duì)衛(wèi)星平臺(tái)的運(yùn)動(dòng)建模，獲得了對(duì)應(yīng)的模糊算子，進(jìn)而用此算子對(duì)多光譜圖像進(jìn)行模糊處理，來(lái)模擬編碼孔徑遙感成像的過(guò)程; 然后再將此模糊算子應(yīng)用于圖像重建過(guò)程中，獲得了很好的重建效果。除此之外，還進(jìn)行了對(duì)比試驗(yàn)——重建時(shí)不加入模糊算子，證明了模糊算子對(duì)重建圖像像質(zhì)的影響。通過(guò)整個(gè)仿真和重建過(guò)程，可知衛(wèi)星遙感平臺(tái)的抖動(dòng)對(duì)光譜成像過(guò)程具有不容忽視的影響，獲得準(zhǔn)確的運(yùn)動(dòng)模型將對(duì)后期的圖像重建工作十分重要。

壓縮編碼孔徑多光譜成像是一項(xiàng)基于壓縮感知(CS)理論發(fā)展起來(lái)的新技術(shù)，其在遙感領(lǐng)域的應(yīng)用仍有很長(zhǎng)的路要走。本文僅對(duì)將編碼孔徑運(yùn)用到遙感多光譜成像所遇到的第一道難題——嚴(yán)重的光譜混疊提出了初步的解決方案，對(duì)于文中所提及的混疊權(quán)重系數(shù)矩陣與焦距和像元尺寸的關(guān)系、重建的優(yōu)化算法及系統(tǒng)模型優(yōu)化等問(wèn)題，仍需進(jìn)一步研究。

References)

[1] Green R O, Eastwood M L, Sarture C M, et al. Imaging Spectroscopy and the Airborne VisibleInfrared Imaging Spectrometer(AVIRIS) [J]. Remote Sens. Environ, 1998, 65(3): 227–248.

[2] Herrala E, Okkonen J T, Hyvarinen T S, et al. Imaging Spectrometer for Process Industry Applications[J].Proc. SPIE, 1994,2248(33): 33–40.

[3] Morris H, Hoyt C, Treado P.Imaging Spectrometers for Fluorescence and Raman Microscopy: Acousto-optic and Liquid Crystal Tunable Filters[J]. Appl. Spec- trosc, 1994, 48(7): 857–866.

[4] Descour M R, Volin C E, Dereniak E L, et al.Demonstration of a High-speed Nonscanning Imaging Spectrometer[J]. Opt. Lett.,1997, 22(16): 1271–1273.

[5] Brady D J, Gehm M E.Compressive Imaging Spectrometers Using Coded Apertures[J]. Proc. SPIE, 2006, 6246: 62460A.

[6] Candès E, Romberg J, Tao T.Robust Uncertainty Principles: Exact Signal Reconstruction from Highly Incomplete Frequency Information[J]. IEEE Trans. Inf. Theory, 2006, 52(2): 489–509.

[7] Candès E, Tao T. Near Optimal Signal Recovery From Random Projections: Universal Encoding Strategies[J]. IEEE Trans. Inf.Theory, 2006, 52(12): 5406–5425.

[8] Donoho D. Compressed Sensing[J]. IEEE Trans. Inf. Theory, 2006, 52(4): 1289–1306.

[9] Donoho D L. Compressed Sensing[J]. IEEE Transactions on Information Theory, 2006, 52(4): 1289–1306.

[10] Candès E, Romberg J, Tao T. Robust Uncertainty Principles: Exact Signal Reconstruction from Highly Incomplete Frequency Information[J]. IEEE Transactions on Information Theory, 2006, 52(2): 489–509.

[11] Candès E. Compressive Sampling[C].In: Proceedings of International Congress of Mathematicians. Zürich, Switzerland: European Mathematical Society Publishing House, 2006: 1433–1452.

[12] 許曉峰.光學(xué)壓縮編碼成像及其復(fù)原算法研究[D]. 西安電子科技大學(xué)碩士研究生學(xué)位論文, 2011.XU Xiaofeng. Optical Compressive Coded Aperture Imaging and Research on the Recovery Algorithm [D] Master thesis of Xidian Univ., 2011. (in Chinese)

[13] Wagadarikar A A, Pitsianis N P, SUN Xiaobai,et al. Spectral Image Estimation for Coded Aperture Snapshot Spectral Imagers[J]. Proc. SPIE, 2008, 7076: 707602.

[14] 唐秋艷, 唐義, 曹瑋亮, 等. 衛(wèi)星平臺(tái)復(fù)雜振動(dòng)引起的光譜成像退化仿真研究[J]. 物理學(xué)報(bào),2012, 61(7): 070202-1-070202-9.TANG Qiuyan, TANG Yi, CAO Weiliang, et al. Simulation of Imaging Spectrometers Degraded by Satellite Vibrations[J].Acta Physica Sinica, 2012, 61(7): 070202-1-070202-9.(in Chinese)

[15] 王曉燕, 唐義, 唐秋艷, 等. 色散型推掃成像光譜儀在軌動(dòng)態(tài)成像建模與分析[J]. 北京理工大學(xué)學(xué)報(bào), 2012, 32(3):291–295.WANG Xiaoyan, TANG Yi, TANG Qiuyan, et al. Dynamic Imaging Modeling and Analysis of Dispersive Spectrometer with Pushbroom Imaging on Orbit [J] Journal of Beijing Institute of Technology, 2012, 32(3): 291–295. (in Chinese)

[16] Bioucas-Dias J M, Figueiredo M A. A New TwIST: Two-step Iterative Shrinkage /Thresholding Algorithms for Image Restoration[J].IEEE Trans Image Process, 2007, 16(12): 2992–3004.