基于變分貝葉斯推理的高光譜圖像恢復

鄒長忠

(福州大學數(shù)學與計算機科學學院,福建 福州 350116)

0 引言

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

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

1 問題描述和估計模型

1.1 高光譜圖像降質模型

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

Y=X+V

(1)

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

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

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

(2)

其中

(3)

1.2 先驗分布

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

(4)

(5)

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

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

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

(6)

容易得到：

(7)

進一步可以展開為：

(8)

(9)

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

2 變分貝葉斯推理

2.1 變分貝葉斯推理

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

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

重復一、 二兩步,直到收斂或者達到停止條件.

2.2 算法

具體的算法描述如圖1.

3 測試數(shù)據(jù)與結果分析

提出的方法與NL,K-SVD,BM3D,Zhao和Rasti算法從衡量指標和視覺效果兩方面比較. 實驗結果看出本文的方法比這些去噪方法效果更好.

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

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

3.2 衡量指標和參數(shù)設置

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

3.3 實驗結果

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

圖2 Washington DC 恢復結果圖Fig.2 Results of several methods for Washington DC

圖3 Cave 恢復結果圖Fig.3 Results of several methods for Cave

圖4 Urban 恢復結果圖Fig.4 Results of several methods for Urban

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

表1 圖像恢復性能指標Tab.1 Compare of recovery of different images

4 結語

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

[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.