鄒長(zhǎng)忠

(福州大學(xué)數(shù)學(xué)與計(jì)算機(jī)科學(xué)學(xué)院,福建 福州 350116)

0 引言

高光譜遙感圖像是一類由一兩百個(gè)波段組成的光譜圖像,覆蓋的譜范圍很廣,一般為400～2 500 nm,因?yàn)橛泻芨吖庾V分辨率,這類圖像得到廣泛應(yīng)用[1-2]. 因?yàn)閭鞲衅鞅旧硇阅軉栴},遙感圖像傳感器得到的高光譜圖像都帶有噪聲[3-4]. 不同于一般的自然圖像,高光譜圖像由于其本身的特性,需要處理的光譜數(shù)量之多,這些都給去噪等問題帶來新的困難. 目前經(jīng)典的去噪方法,有基于TV 正則方法[5]； 利用非局部相似性方法來模擬圖形的非局部方法 (NL)[6-8]. 基于圖像的稀疏表示方法,如Elad等[9-10]提出的的K-SVD方法,該方法基本思想是通過學(xué)習(xí)到的字典和OMP方法[11]求解稀疏因子來估計(jì)去噪后的圖像,另一個(gè)圖像恢復(fù)算法是BM3D[12],該方法充分探索光譜圖像的三維表示,通過非局部相似性方法來去除噪聲. 當(dāng)然這些方法直接應(yīng)用在處理高光譜圖像去噪問題上,通過對(duì)高光譜圖像的各波段獨(dú)立處理,但不能充分考慮高光譜圖像內(nèi)在的結(jié)構(gòu)關(guān)系. 考慮高光譜波段之間特性的一些去噪方法也有相關(guān)的文獻(xiàn)報(bào)道,如Chen等[13]提出基于主成分分析和小波收縮方法的遙感圖像去噪,Yuan等[14]提出基于空-譜自適應(yīng)全變分模型去除光譜圖像噪聲. Zhao等[15]提出基于稀疏表示和低秩約束的高光譜圖像去噪方法,Rasti等[16]提出基于小波域的光譜正則方法解決高光譜去噪問題.

本文構(gòu)建高光譜圖像的表示模型,充分挖掘圖像空—譜特性,建立圖像恢復(fù)模型,達(dá)到較好的去噪效果. 首先,建立高光譜圖像的降質(zhì)表示模型,利用圖像在變換域后的稀疏特性,通過小波基矩陣作為變換矩陣,構(gòu)建均值為0的高斯分布作為小波因子的先驗(yàn)分布,然后建立最大后驗(yàn)估計(jì)模型,最后提出變分貝葉斯推導(dǎo)方法. 通過實(shí)際的高光譜遙感圖像實(shí)驗(yàn),提出的方法與K-SVD[10]算法 、 非局部算法NL[6]、 BM3D[12]算法、 Zhao等[15]和Rasti等[16]算法等作比較.

1 問題描述和估計(jì)模型

1.1 高光譜圖像降質(zhì)模型

為整體描述一個(gè)高光譜圖像,將每個(gè)波段表示成一個(gè)向量,作為這個(gè)高光譜圖像的行向量,實(shí)際的觀測(cè)圖像受高斯噪聲影響,從而可描述圖像降質(zhì)模型為：

Y=X+V

(1)

式中：X∈RN×n表示真實(shí)的圖像,其中N,n分別表示高光譜圖像波段數(shù)和每個(gè)波段對(duì)應(yīng)的像素個(gè)數(shù)；Y∈RN×n表示觀測(cè)到的降質(zhì)圖像；V∈RN×n表示高斯噪聲. 我們的問題是,由觀測(cè)到的圖像Y, 如何有效地估計(jì)X. 本文通過貝葉斯最大后驗(yàn)估計(jì)方法來估計(jì)需要恢復(fù)的圖像.

由觀測(cè)的圖像Y估計(jì)X,首先建立最大似然函數(shù)max {P(Y|X)},最大化似然函數(shù)等價(jià)于最小化負(fù)對(duì)數(shù)似然函數(shù),即： min{-logP(Y|X)}. 模型中考慮的是高斯噪聲,所以噪聲服從高斯分布,即V～N(0,Rv),其中Rv是噪聲方差,進(jìn)一步假設(shè)方差已知,那么似然函數(shù)服從高斯分布為：P(Y|X)～N(X,Rv).

近來,稀疏表示[18]已經(jīng)成為非常有效地表示圖像結(jié)構(gòu)特性的一種先驗(yàn)表示. 本文采用小波基[19]作為變換矩陣,關(guān)于小波基的構(gòu)造在實(shí)驗(yàn)部分會(huì)加以說明,通過引入小波基矩陣W,X可以表示為X=WZ,原問題轉(zhuǎn)化為Y=WZ+V. 在估計(jì)好Z后,通過計(jì)算X=WZ, 即可得到最后估計(jì)的圖像. 以列向量形式來表示,記Y=[y1,y2, …,yn],X=[x1,x2, …,xn],V=[v1,v2, …,vn],那么對(duì)于每個(gè)列向量yi,有yi=Wzi+vi,則似然函數(shù)具體表示為：

(2)

其中

(3)

1.2 先驗(yàn)分布

由于圖像在小波變換后具有稀疏性,所以Z的大部分因子取值很小,或者接近于0,為了較好地刻畫Z的稀疏性,假設(shè)Z=[z1,z2, …,zn]中的每個(gè)zi服從均值為0,方差為Rz的高斯分布,即：P(zi|Rz)～N(0,Rz). 其中定義Rz=diag(r1-1,r2-1,…,rn-1),每個(gè)ri用于描述相應(yīng)波段的方差,此外為了更好地估計(jì)Z,假設(shè)協(xié)方差矩陣Rz中的每個(gè)ri服從伽馬分布,即ri～Gamma(α0,β0),取伽馬分布有兩個(gè)原因,一是該分布可以比較好地模擬方差的變化,另一個(gè)原因是該分布是共軛先驗(yàn)分布,以便于后面的概率推導(dǎo). 進(jìn)一步假設(shè)zi之間獨(dú)立,那么Z的先驗(yàn)概率具體為：

(4)

(5)

相應(yīng)的聯(lián)合分布函數(shù)表達(dá)為：P(Z,Rz,Y,W)=P(Y|WZ)P(Z|Rz)P(Rz) .

為方便計(jì)算,根據(jù)公式(2)～(4),進(jìn)一步求其聯(lián)合分布函數(shù)的對(duì)數(shù)形式：

根據(jù)貝葉斯估計(jì)原則,需要最大后驗(yàn)密度函數(shù),目標(biāo)函數(shù)表示為：

(6)

容易得到：

(7)

進(jìn)一步可以展開為：

(8)

(9)

式(9)是似然函數(shù)及先驗(yàn)概率和超參數(shù)先驗(yàn)概率函數(shù)的組合,對(duì)其取對(duì)數(shù)形式對(duì)應(yīng)到優(yōu)化問題就是反映高斯噪聲的數(shù)據(jù)保真項(xiàng)和正則項(xiàng)的組合. 本文采用變分貝葉斯方法推導(dǎo).

2 變分貝葉斯推理

2.1 變分貝葉斯推理

第一步,求Z的后驗(yàn)密度函數(shù)：

那么根據(jù)伽馬分布性質(zhì),計(jì)算其均值為：

重復(fù)一、 二兩步,直到收斂或者達(dá)到停止條件.

2.2 算法

具體的算法描述如圖1.

3 測(cè)試數(shù)據(jù)與結(jié)果分析

提出的方法與NL,K-SVD,BM3D,Zhao和Rasti算法從衡量指標(biāo)和視覺效果兩方面比較. 實(shí)驗(yàn)結(jié)果看出本文的方法比這些去噪方法效果更好.

3.1 測(cè)試數(shù)據(jù)

三個(gè)實(shí)際的高光譜圖像被測(cè)試. 第一個(gè)數(shù)據(jù)是拍攝于華盛頓(WashingtonDCMall),由AVIRIS[20]傳感器在1996年拍攝,該數(shù)據(jù)有191個(gè)譜波段,覆蓋范圍為0.4～2.5μm, 選取每個(gè)波段構(gòu)成256×256個(gè)像素,使用全部波段. 第二個(gè)數(shù)據(jù)是Cave數(shù)據(jù)庫(kù)[21],該數(shù)據(jù)是日常對(duì)象的高光譜圖像,共有32個(gè)場(chǎng)景,每個(gè)場(chǎng)景有512×512×31個(gè)大小,選取其中一個(gè)場(chǎng)景名為oil_painting來做實(shí)驗(yàn). 第三個(gè)數(shù)據(jù)是HYDICEUrban數(shù)據(jù)庫(kù),該數(shù)據(jù)拍攝于Urban地區(qū),在去除水蒸氣波段后,使用的測(cè)試數(shù)據(jù)大小為128×128×163.

3.2 衡量指標(biāo)和參數(shù)設(shè)置

衡量指標(biāo)有：PSNR,RMSE,SAM[22]. 其中PSNR衡量恢復(fù)的圖像的信噪比,該性能指標(biāo)值越大越好；RMSE衡量恢復(fù)的圖像的最小均方誤差,越小越好；SAM是恢復(fù)的圖像的譜扭曲程度,越小越好. 所有的實(shí)驗(yàn)都是在Intel(R)Core(TM)i7-3537UCPU@2.50GHzand8GBRAM機(jī)子上,用MatlabR2012b編程實(shí)現(xiàn). 實(shí)驗(yàn)添加的噪聲方差為0.04和0.06. 本文采用db1小波構(gòu)造小波基矩陣W.

3.3 實(shí)驗(yàn)結(jié)果

因?yàn)椴豢赡芙o出所有波段的結(jié)果圖,采用偽彩色方法來表示去噪結(jié)果. 通過選擇其中的三個(gè)波段,給出噪聲方差為0.06的結(jié)果圖. 選擇顯示W(wǎng)ashington DC結(jié)果圖的波段是50、 40、 25 ,Cave的波段是26、 14、 6,Urban的波段是50、 30、 10. 結(jié)果圖分別顯示在圖2～4. 從這些效果圖可看出NL和K-SVD方法恢復(fù)的圖像對(duì)光譜扭曲比較大,特別是圖3恢復(fù)圖像的左半部分顏色很明顯存在與真實(shí)圖像有較大差異,這是因?yàn)镹L和K-SVD方法都是孤立地對(duì)每個(gè)波段進(jìn)行去噪,這樣沒有考慮波段之間的相關(guān)性,在去噪時(shí)不能充分地利用圖像光譜之間的相似特性,也就是在求解圖像去噪反問題時(shí),不能有效地利用先驗(yàn)信息對(duì)最后的解加以約束,從而使得恢復(fù)出來的圖像失去了譜間的相關(guān)性,造成光譜失真嚴(yán)重. 而相對(duì)BM3D,Zhao和Rasti方法,這些方法處理時(shí)因?yàn)榭紤]圖像三維信息,所以譜扭曲相對(duì)小,但去噪效果一般,這些恢復(fù)圖像仍然存在不少噪聲,此外Zhao和Rasti方法去噪效果不均勻,造成局部噪聲比較大,所以存在圖像局部模糊現(xiàn)象,如圖2和圖4所示. 提出的方法結(jié)果圖,從噪聲和譜扭曲上看明顯優(yōu)于其他三個(gè)方法,效果非常接近真實(shí)圖像. 表 1給出了恢復(fù)的量化指標(biāo). 從表1看出本文的方法恢復(fù)的圖像結(jié)果噪聲比較小,另外衡量恢復(fù)圖像的最小均方誤差RMSE也對(duì)應(yīng)很小,針對(duì)高光譜圖像,另一個(gè)很重要的指標(biāo)SAM,用以表示圖像的譜扭曲程度,提出方法使得最后圖像的譜扭曲小,這個(gè)與實(shí)驗(yàn)結(jié)果圖一致,與我們的模型考慮高光譜圖像譜間關(guān)系一致.

圖2 Washington DC 恢復(fù)結(jié)果圖Fig.2 Results of several methods for Washington DC

圖3 Cave 恢復(fù)結(jié)果圖Fig.3 Results of several methods for Cave

圖4 Urban 恢復(fù)結(jié)果圖Fig.4 Results of several methods for Urban

綜上,由這些指標(biāo)和圖像實(shí)驗(yàn)結(jié)果圖,可以看出本文的方法明顯優(yōu)于其他三個(gè)方法. 分析其主要原因是處理恢復(fù)問題,把高光譜圖像作為一個(gè)整體來建模,另外小波基作用下能很好地表示因子的稀疏性,構(gòu)建小波稀疏因子的先驗(yàn)信息能很好地模擬圖像在變換域后的稀疏特性,所以本文的方法能有效去除噪聲,同時(shí)還能保持原來光譜間的關(guān)系,重構(gòu)的結(jié)果非常地接近真實(shí)圖像.

表1 圖像恢復(fù)性能指標(biāo)Tab.1 Compare of recovery of different images

4 結(jié)語(yǔ)

提出一種基于變分貝葉斯推理的高光譜圖像恢復(fù)方法,該方法將高光譜圖像恢復(fù)問題描述成一個(gè)病態(tài)反問題. 為建立優(yōu)化估計(jì)模型,首先建立高光譜圖像含噪的觀測(cè)模型,然后提出用于描述高斯噪聲的二范數(shù)保真項(xiàng),為約束該病態(tài)問題的解空間,進(jìn)一步提出基于小波變換后因子的稀疏先驗(yàn)分布作為正則項(xiàng),從而建立最大后驗(yàn)聯(lián)合分布函數(shù). 該聯(lián)合優(yōu)化模型包括用于描述高斯分布的數(shù)據(jù)保真項(xiàng)和描述小波因子稀疏特性的0均值高斯分布,其中小波因子分布的方差作為未知的超參數(shù)估計(jì),從而具備更好的解空間搜索能力. 針對(duì)該優(yōu)化模型的解析表達(dá)式求解困難問題,最后采用變分貝葉斯方法推導(dǎo),與采樣方法相比,變分貝葉斯能快速地收斂到極小值. 為驗(yàn)證算法的有效性,利用實(shí)際的三幅公認(rèn)高光譜圖像做實(shí)驗(yàn),實(shí)驗(yàn)結(jié)果從效果衡量指標(biāo)和實(shí)驗(yàn)結(jié)果視覺圖兩方面驗(yàn)證本文提出的方法優(yōu)于目前經(jīng)典的恢復(fù)方法.

[1] MANOLAKIS D,SHAW G. Detection algorithms for hyperspectral imaging applications[J]. IEEE Signal Processing Magazine. 2002,19(1): 29-43.

[2] CHANG C I. Hyperspectral data exploitation: theory and application[M]. New Jersey: Wiley,2007.

[3] AIAZZI B,ALPARONE L, BARDUCCI A,etal. Information theoretic assessment of sampled hyperspectral imagers[J]. IEEE Transactions Geoscience and Remote Sensing,2001,39(7): 1447-1458.

[4] KEREKES J,BAUM J. Full-spectrum spectral imaging system analytical model[J]. IEEE Transactions Geoscience and Remote Sensing,2005,5(2): 571-580.

[5] RUDIN L,OSHER S,FATEMI E. Nonlinear total variation based noise removal algorithms[J]. Physica D: Nonlinear Phenomena,1992,60(1): 259-268.

[6] KINDERMANN S,OSHER S,JONES P. Deblurring and denoising of images by nonlocal functionals[J]. SIAM Multiscale Modeling and Simulation,2005,4(4): 1091-1115.

[7] BOULANGER J,KERVRANN C,BOUTHEMY P,etal. Patch-based nonlocal functional for denoising fluorescence microscopy image sequences[J]. IEEE Transactions Medical Imaging,2010, 29(2): 442-454.

[8] DURAN J,BUADES A,COLL B,etal. A nonlocal variational model for pansharpening image fusion[J]. SIAM Journal on Imaging Sciences,2014,7(2): 761-796.

[9] STARCK J L,ELAD M,DONOHO D L. Image decomposition via the combination of sparse representations and a variational approach[J]. IEEE Transactions on Image Processing,2005,14(10): 1570-1582.

[10] AHARON M,ELAD M,BRUCKSTEIN A. K-SVD: an algorithm for designing of overcomplete dictionaries for sparse representation[J]. IEEE Transactions Signal Processing, 2006,54(11): 4311-4322.

[11] TROPP J. Greed is good: algorithmic results for sparse approximation[J]. IEEE Transactions on Information Theory,2004,50(10): 2231-2242.

[12] DABOV K,FOI A,KATKOVNIK V,etal. Image denoising by sparse 3d transform-domain collaborative filtering[J]. IEEE Transactions Image Processing,2007,16(8): 2080-2094.

[13] CHEN G,QIAN S. Denoising of hyperspectral imagery using principal component analysis and wavelet shrinkage[J]. IEEE Transactions on Geoscience and Remote Sensing,2011,49(3): 973-980.

[14] YUAN Q,ZHANG L,SHEN H. Hyperspectral image denoising employing a spectral-spatial adaptive total variation model[J]. IEEE Transactions on Geoscience and Remote Sensing,2012,50(10): 3660-3677.

[15] ZHAO Y,YANG J. Hyperspectral image denoising via sparse representation and low-rank constraint[J]. IEEE Transactions on Geoscience and Remote Sensing. 2015,53(1): 296-308.

[16] RASTI B,SVEINSSON J R,Ulfarsson M,etal. Hyperspectral image denoising using first order spectral roughness penalty in wavelet do-main[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing,2014,7(6): 2458-2467.

[17] BIOUCAS J,PLAZA A,DOBIGEON,etal. Hyperspectral unmixing overview: geometrical,statistical,and sparse regression-based approaches[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing,2012,5(2): 354-379.

[18] DAUBECHIES I,DEFRIESE M,MOL C. An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[J]. Commun Pure Appl Math ,2004,57(11): 1413-1457.

[19] DONOHO D,JOHNSTONE I. Ideal spatial adaptation by wavelet shrinkage[J].Biometrika,1994,81(3): 425-455.

[20] VANE G,GREEN R O,CHRIEN T G,etal. The airborne visible/infrared imaging spectrometer(AVIRIS)[J]. Remote Sens Environ Jun,1993,44(2): 127-143.

[21]YASUMA F,MITSUNAGE T,ISO D,etal. Generalized assorted pixel camera: post-capture control of resolution,dynamic range and spectrum[J]. IEEE Transactions on Image Processing,2010,19(9): 2241-2253.

[22] WALD L. Quality of high resolution synthesised images: is there a simple criterion[C]// Proc 3rd Conf on Fusion of Earth Data. SEE/URISCA: Sophia Antipolis,2000: 99-103.