基于变分贝叶斯推理的高光谱图像恢复

邹长忠

(福州大学数学与计算机科学学院,福建 福州 350116)

0 引言

高光谱遥感图像是一类由一两百个波段组成的光谱图像,覆盖的谱范围很广,一般为400～2 500 nm,因为有很高光谱分辨率,这类图像得到广泛应用[1-2]. 因为传感器本身性能问题,遥感图像传感器得到的高光谱图像都带有噪声[3-4]. 不同于一般的自然图像,高光谱图像由于其本身的特性,需要处理的光谱数量之多,这些都给去噪等问题带来新的困难. 目前经典的去噪方法,有基于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分别表示高光谱图像波段数和每个波段对应的像素个数；Y∈RN×n表示观测到的降质图像；V∈RN×n表示高斯噪声. 我们的问题是,由观测到的图像Y, 如何有效地估计X. 本文通过贝叶斯最大后验估计方法来估计需要恢复的图像.

由观测的图像Y估计X,首先建立最大似然函数max {P(Y|X)},最大化似然函数等价于最小化负对数似然函数,即： min{-logP(Y|X)}. 模型中考虑的是高斯噪声,所以噪声服从高斯分布,即V～N(0,Rv),其中Rv是噪声方差,进一步假设方差已知,那么似然函数服从高斯分布为：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,则似然函数具体表示为：

(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,假设协方差矩阵Rz中的每个ri服从伽马分布,即ri～Gamma(α0,β0),取伽马分布有两个原因,一是该分布可以比较好地模拟方差的变化,另一个原因是该分布是共轭先验分布,以便于后面的概率推导. 进一步假设zi之间独立,那么Z的先验概率具体为：

(4)

(5)

相应的联合分布函数表达为：P(Z,Rz,Y,W)=P(Y|WZ)P(Z|Rz)P(Rz) .

为方便计算,根据公式(2)～(4),进一步求其联合分布函数的对数形式：

根据贝叶斯估计原则,需要最大后验密度函数,目标函数表示为：

(6)

容易得到：

(7)

进一步可以展开为：

(8)

(9)

式(9)是似然函数及先验概率和超参数先验概率函数的组合,对其取对数形式对应到优化问题就是反映高斯噪声的数据保真项和正则项的组合. 本文采用变分贝叶斯方法推导.

2 变分贝叶斯推理

2.1 变分贝叶斯推理

第一步,求Z的后验密度函数：

那么根据伽马分布性质,计算其均值为：

重复一、 二两步,直到收敛或者达到停止条件.

2.2 算法

具体的算法描述如图1.

3 测试数据与结果分析

提出的方法与NL,K-SVD,BM3D,Zhao和Rasti算法从衡量指标和视觉效果两方面比较. 实验结果看出本文的方法比这些去噪方法效果更好.

3.1 测试数据

三个实际的高光谱图像被测试. 第一个数据是拍摄于华盛顿(WashingtonDCMall),由AVIRIS[20]传感器在1996年拍摄,该数据有191个谱波段,覆盖范围为0.4～2.5μm, 选取每个波段构成256×256个像素,使用全部波段. 第二个数据是Cave数据库[21],该数据是日常对象的高光谱图像,共有32个场景,每个场景有512×512×31个大小,选取其中一个场景名为oil_painting来做实验. 第三个数据是HYDICEUrban数据库,该数据拍摄于Urban地区,在去除水蒸气波段后,使用的测试数据大小为128×128×163.

3.2 衡量指标和参数设置

衡量指标有：PSNR,RMSE,SAM[22]. 其中PSNR衡量恢复的图像的信噪比,该性能指标值越大越好；RMSE衡量恢复的图像的最小均方误差,越小越好；SAM是恢复的图像的谱扭曲程度,越小越好. 所有的实验都是在Intel(R)Core(TM)i7-3537UCPU@2.50GHzand8GBRAM机子上,用MatlabR2012b编程实现. 实验添加的噪声方差为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方法去噪效果不均匀,造成局部噪声比较大,所以存在图像局部模糊现象,如图2和图4所示. 提出的方法结果图,从噪声和谱扭曲上看明显优于其他三个方法,效果非常接近真实图像. 表 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

综上,由这些指标和图像实验结果图,可以看出本文的方法明显优于其他三个方法. 分析其主要原因是处理恢复问题,把高光谱图像作为一个整体来建模,另外小波基作用下能很好地表示因子的稀疏性,构建小波稀疏因子的先验信息能很好地模拟图像在变换域后的稀疏特性,所以本文的方法能有效去除噪声,同时还能保持原来光谱间的关系,重构的结果非常地接近真实图像.

表1 图像恢复性能指标Tab.1 Compare of recovery of different images

4 结语

提出一种基于变分贝叶斯推理的高光谱图像恢复方法,该方法将高光谱图像恢复问题描述成一个病态反问题. 为建立优化估计模型,首先建立高光谱图像含噪的观测模型,然后提出用于描述高斯噪声的二范数保真项,为约束该病态问题的解空间,进一步提出基于小波变换后因子的稀疏先验分布作为正则项,从而建立最大后验联合分布函数. 该联合优化模型包括用于描述高斯分布的数据保真项和描述小波因子稀疏特性的0均值高斯分布,其中小波因子分布的方差作为未知的超参数估计,从而具备更好的解空间搜索能力. 针对该优化模型的解析表达式求解困难问题,最后采用变分贝叶斯方法推导,与采样方法相比,变分贝叶斯能快速地收敛到极小值. 为验证算法的有效性,利用实际的三幅公认高光谱图像做实验,实验结果从效果衡量指标和实验结果视觉图两方面验证本文提出的方法优于目前经典的恢复方法.

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