计算机应用   2016, Vol. 36 Issue (12): 3411-3417  DOI: 10.11772/j.issn.1001-9081.2016.12.3411
0

引用本文 

汪浩然, 夏克文, 任苗苗, 李绰. 结合PCA及字典学习的高光谱图像自适应去噪方法[J]. 计算机应用, 2016, 36(12): 3411-3417.DOI: 10.11772/j.issn.1001-9081.2016.12.3411.
WANG Haoran, XIA Kewen, REN Miaomiao, LI Chuo. Adaptive Denoising of Hyperspectral Remote Sensing Image Based on Principal Component Analysis and Dictionary Learning with Noise Estimation[J]. JOURNAL OF COMPUTER APPLICATIONS, 2016, 36(12): 3411-3417. DOI: 10.11772/j.issn.1001-9081.2016.12.3411.

基金项目

国家自然科学基金资助项目(51208168);天津市自然科学基金资助项目(13JCYBJC37700);河北省自然科学基金资助项目(E2016202341);大学生创新创业训练计划项目(河北省重点)(201510080051)

通信作者

夏克文(1956-),男,湖南邵阳人,教授,博士,主要研究方向:石油测井数据挖掘、智能天线953833825@qq.com

作者简介

汪浩然(1990-),男,天津人,硕士研究生,主要研究方向:高光谱图像处理、压缩感知;
任苗苗(1992-),女,河北石家庄人,硕士研究生,主要研究方向:合成孔径雷达图像处理、压缩感知;
李绰(1995-),女,河北石家庄人,主要研究方向:粒计算、粗糙集

文章历史

收稿日期:2016-05-23
修回日期:2016-07-04
结合PCA及字典学习的高光谱图像自适应去噪方法
汪浩然1, 夏克文1, 任苗苗2,3, 李绰1    
1. 河北工业大学 电子与信息工程学院, 天津 300401 ;
2. 中国科学院 电子学研究所, 北京 100190 ;
3. 中国科学院大学, 北京 100190
摘要: 高光谱图像各波段图像噪声分布复杂,传统去噪方法难以达到理想效果。针对这一问题,在主成分分析(PCA)的基础上,结合噪声估计和字典学习,提出一种新的高光谱去噪方法。首先,对原始高光谱数据进行主成分变换得到一组主成分图像并根据能量比重将其划分为清晰图像组和含噪图像组;然后,根据任一波段图像的信息,利用奇异值分解(SVD)对图像进行噪声估计,再将得到的噪声估计方法与K-SVD字典学习去噪算法结合,提出一种具备自适应噪声估计特性的字典学习去噪算法,并将其应用于信息量较小的含噪图像组进行去噪处理;最后,按各主成分图像对应的信息量比例进行加权融合得到最终的去噪图像。通过对模拟与实际高光谱遥感图像的实验表明,与PCA、PCA-Bish、PCA-Contourlet三种去噪方法相比,所提方法去噪后图像的峰值信噪比(PSNR)可以提升1~3 dB,且具有更多的细节信息和更好的视觉效果。
关键词: 高光谱遥感    主成分分析    噪声估计    奇异值分解    字典学习    
Adaptive Denoising of Hyperspectral Remote Sensing Image Based on Principal Component Analysis and Dictionary Learning with Noise Estimation
WANG Haoran1, XIA Kewen1, REN Miaomiao2,3, LI Chuo1     
1. School of Electronics and Information Engineering, Hebei University of Technology, Tianjin 300401, China ;
2. Institute of Electronics, Chinese Academy of Sciences, Beijing 100190, China ;
3. University of Chinese Academy of Sciences, Beijing 100190, China
Abstract: The distributed state of noise existing among different bands of hyperspectral remote sensing image is complex, so the traditional denoising methods are hard to achieve the desired effect. In order to solve this problem, based on Principal Component Analysis (PCA), a novel denoising method for hyperspectral data was proposed combining with noise estimation and dictionary learning. Firstly, a group of the principal component images were achieved from the original hyperspectral data by using the PCA transform, which were divided into clear image group and noisy image group according to the corresponding energy. Then, according to any band image from noisy hyperspectral data, the noise standard deviation of the image was estimated via a noise estimation method based on Singular Value Decomposition (SVD). Meanwhile, combining this noise estimation method with denoising method via K-SVD dictionary learning, a new dictionary learning denoising method with adaptive noise estimation characteristics was proposed and applied to denoise those images from noisy image group with low energy where noise mainly existed. Finally, the final denoising image was obtained by weighted fusion according to the corresponding energy of each principal component image. The experimental results on simulated and real hyperspectral remote sensing data show that, compared with PCA, PCA-Bish and PCA-Contourlet, the Peak Signal-to-Noise Ratio (PSNR) of the image denoised by the proposed algorithm is improved by 1-3 dB, and more detailed information and better visual effect of the denoised image by the proposed method are achieved.
Key words: hyperspectral remote sensing    Principal Component Analysis (PCA)    noise estimation    Singular Value Decomposition (SVD)    dictionary learning    
0 引言

近些年来,高光谱遥感图处理技术得到了广泛关注,并已被成功应用于农业、林业、环境监测、考古探索、军事侦察等多个领域中[1]。高光谱遥感图像可以看作是由一维光谱信息和二维空间信息联合组成的三维数据块。然而,由于太阳光经地面反射后穿过大气时会发生散射,散射光到达传感器就会对传感器接收的图像造成噪声污染;同时,高光谱图像的电磁波在传播过程中会受到多重因素的影响,也会导致高光谱在成像时遭到严重的噪声干扰。这些噪声污染必然会对后续的解混[2]和分类[3]等环节造成严重影响。因此,高光谱图像去噪技术在高光谱遥感领域,乃至整个遥感学科中都占据着重要地位。

与单幅遥感图像、全色图像和多光谱图像相比,高光谱图像的噪声去除问题更加复杂。一方面因素单个波段图像所含的光谱频率范围较窄,导致各波段图像间有着很强的相关性;另一方面因素是高光谱图像的波段数较多,每个波段图像所受噪声的影响程度不同。上述特性使得传统的遥感图像去噪方法并不适用于高光谱图像的噪声去除,因此针对高光谱遥感图像的去噪技术研究是很有必要的。

传统的高光谱遥感图像去噪方法基本可以分为两大类:1) 在空间域或频率域进行噪声去除。如Yuan等[4]提出的光谱-空间域自适应全变分去噪模型,Rasti等[5]提出了基于稀疏降秩回归小波的高光谱图像去噪方法, 此类方法充分考虑了高光谱数据的二维空间信息,但并没能充分利用高光谱数据各波段间的光谱信息。2) 用特征降维方法进行去噪。如Stephan等[6]利用主成分分析(Principal Component Analysis, PCA)进行高光谱图像的消噪处理, Green等[7]提出最小噪声分离(Minimum Noise Fraction, MNF)理论, 该类方法充分考虑了高光谱图像各波段间的相关性,但由于缺乏对单张图像本身细节(即二维空间信息)的考虑,导致去噪效果也不甚理想。

考虑到以上两类去噪方法的缺陷,近些年出现了将特征降维技术和空间域去噪技术结合使用的去噪方法。此类方法同时考虑了高光谱图像在空间域的细节信息和多个波段间的相关性,比起前两类方法能够实现更好的去噪效果。Chen等[8]提出了利用PCA对数据进行降维,保留能量高的主成分,而对能量低的主成分,利用二维的二元小波阈值收缩方法进行去噪。之后,Chen等[9]又提出了结合PCA和块匹配4D(Block-Matching 4D, BM4D)滤波器的高光谱图像消噪方法。常威威等[10]利用PCA变换在高维数据特征降维中的高效性和Contourlet变换良好的去噪能力,提出了基于PCA和Contourlet变换的高光谱遥感图像消噪方法,得到了较好的去噪效果。上述方法均在PCA的基础上结合了多尺度几何变换技术应用于高光谱图像去噪,但应用多尺度几何变换对图像进行稀疏表示时得到的字典具有特定的解析公式,能够隐含地快速实现;但其对应的稀疏字典结构是固定的,没有充分考虑图像本身的属性,这也导致其在对图像内容特殊的高光谱遥感图像去噪时,自适应性受到局限。

因此,针对上述缺陷,本文提出一种结合特征提取方法主成分分析和字典学习并具有噪声估计功能的高光谱图像去噪方法。一方面,充分利用主成分分析对高维数据高效的信息提取能力,使图像的主要信息得以保留;另一方面,对每个噪声成分图像进行自适应的噪声强度估计和字典学习,再结合稀疏表示方法对其进行去噪。

1 自适应噪声估计特性的K-SVD去噪算法 1.1 基于K-SVD的图像去噪模型

稀疏表示目前广泛应用于图像去噪领域,而如何构造合适的字典是进行稀疏表示的重要环节。字典构造方法主要可以分为两大类:一类是固定字典构造方法,常用的如离散余弦变换(Discrete Cosine Transform, DCT)字典、小波字典等。此类方法的缺点在于,使用同一种稀疏基对细节特征差别较大的图像进行稀疏表示时效果差距较大。第二类方法是设计自适应性的冗余字典。此类方法弥补了第一类方法的不足,能够根据图像自身的特点,自适应性地选择稀疏基,从而实现对各类图像的精确稀疏逼近。其中,Aharon等[11]提出的K-奇异值分解(K-Singular Value Decomposition, K-SVD)算法是目前最具代表性、应用最广泛的自适应字典学习算法。以K-SVD算法为基础,Elad等[12]又提出了利用K-SVD算法进行图像去噪的去噪模型。

假设Y=X+V,其中,Y为大小为N×N的含噪图像,X为待重构图像,V为均值为0、标准差为σ的高斯白噪声信号。基于K-SVD算法的图像去噪模型所对应的能量形式为:

$\begin{array}{l} {\rm{\{ }}{\mathop {\bf{\alpha }}\limits^ \wedge _{{\bf{ij}}}}{\bf{,}}\mathop {\bf{D}}\limits^ \wedge {\bf{,}}\mathop {\bf{X}}\limits^ \wedge {\rm{\} }}{\bf{ = }}{\rm{arg}}\mathop {{\rm{min}}}\limits_{{{\bf{\alpha }}_{{\bf{ij}}}}{\bf{,D,X}}} \lambda \left\| {{\bf{X}} - {\bf{Y}}} \right\|_{\rm{2}}^{\rm{2}}\\ + \sum\limits_{ij} {{\mu _{ij}}{{\left\| {{{\bf{\alpha }}_{ij}}} \right\|}_0}} + \sum\limits_{ij} {\left\| {{\bf{D}}{{\bf{\alpha }}_{ij}} - {{\bf{R}}_{ij}}{\bf{X}}} \right\|_2^2} \end{array}$ (1)

模型(1) 中,式中的第1项代表重构图像X和含噪图像Y之间的误差;第2项和第3项为确保X的图像块RijX的稀疏表示具有确定误差界的先验条件。其中,RijX代表图像X被分割为i×j个大小为$\sqrt n \times \sqrt n $图像块中的第(i, j)个图像块。D∈Rn×k为稀疏字典,k是字典原子的个数。αij为RijX在字典D上的稀疏表示系数。 $\mathop {{{\bf{\alpha }}_{{\rm{ij}}}}}\limits^ \wedge ,\mathop {\bf{D}}\limits^ \wedge ,\mathop {\bf{X}}\limits^ \wedge $分别代表第(i, j)个图像块对应的稀疏表示系数、学习得到的字典和去噪后的图像。λ对应图像X和Y之间的保真度,与噪声强度成反比,原文中取λ=30/σ。 μij为稀疏度‖α‖0的权重。

1.2 基于奇异值分解的噪声估计方法

文献[11]中提出的利用K-SVD构造自适应字典的去噪方法已经证明对各类图像具有很强的普适性和很好的去噪效果。但该算法的使用有个前提,即图像的噪声强度是已知的。在实际应用中,需要先对未知图像的噪声强度进行估计,才能继续进行去噪操作, 因此本文从噪声估计这一角度出发,利用SVD方法对图像进行噪声估计研究。

在大多数情况下,高光谱图像噪声均可以用高斯分布来进行建模。本文算法也假设图像的噪声为加性高斯白噪声(Additive White Gaussian Noise,AWGN),即噪声幅度的概率密度函数为:

$f(x) = {e^{ - {x^2}/2{\sigma ^2}}}/\left( {\sqrt {2\pi } \sigma } \right)$ (2)

式中,σ为噪声幅度的标准差。

1.2.1 奇异值分解和基于奇异值分解的噪声估计

奇异值分解(SVD)在数据压缩、图像处理和模式识别等多个领域都得到了广泛的应用。对图像矩阵进行奇异值分解:设XN为大小为m×n、秩为r、概率密度函数满足式(2) 的高斯白噪声图像矩阵,则XN的奇异值分解为:

${{\bf{X}}_N} = {{\bf{U}}_N} \times {{\bf{S}}_N} \times {{\bf{V}}_N}^T$ (3)

式中,SN为XN的奇异值矩阵,由XN XN T的特征值的平方根按降序排列组成。使用si(i=0, 1, …, r)表示SN中的每个奇异值。

同时,由SVD的性质,易知:

$\sigma = \sum\limits_{j = 1}^r {{{\bf{S}}_N}(j)} $ (4)

式中:用L表示选用的尾部奇异值的个数,用函数 PL 表示L个尾部奇异值的平均值,自变量为噪声标准差σ,即:

${{\mathop{\rm P}\nolimits} _L}(\sigma ) = \frac{1}{L}\sum\limits_{j = r - L + 1}^r {{s_N}(j)} ,1 \le L \le r$ (5)

文献[13]的研究表明:当L>>1时(通常取L≥r/4,PL (σ)与σ大致呈线性关系,即:

${{\mathop{\rm P}\nolimits} _L}(\sigma ) = k\sigma ,L > > 1$ (6)

式中,k为线性函数的斜率。

在上述结论基础上可作进一步推广。对含AWGN的图像X,由于噪声和原始图像的加性关系,可进一步将X的PL 分解为PLS和PLN两部分,二者分别代表原图像和噪声对PL 的贡献。又由于原图像的信息量是确定的,根据SVD的性质可知,PLS的取值应为一未知常数,即:

${P_{LS}} = {\rm{ }}b{\rm{;}}b为常数$ (7)

根据式(6) ~(7) 可知:

${{\mathop{\rm P}\nolimits} _L} = {{\mathop{\rm P}\nolimits} _{L{\mathop{\rm S}\nolimits} }} + {{\mathop{\rm P}\nolimits} _{L{\mathop{\rm N}\nolimits} }} = ks + b$ (8)

由此,可以得出一个结论:对含噪图像进行SVD,其尾部L个奇异值的平均值(即PL )与所含AWGN的标准差(即σ)之间大致呈线性关系。要根据式(8) 进行噪声估计,就要同时确定噪声标准差σ和参数k、b。可以通过对含噪图像分别加入2组不相关的已知噪声强度的AWGN,从而构建一个3元3次方程组来确定σ、k、b的值。

现假设图像中含有加性高斯白噪声N,噪声标准差是σ。若对含噪图像再加入与N不相关,强度为σ1的加性高斯白噪声N1。可以证明:加噪后图像所含噪声Nsum1仍是AWGN,且其标准差${\sigma _{sum1}} = \sqrt {{\sigma ^2} + {\sigma _1}^2} $

证明

$\begin{array}{l} {\sigma _{sum1}} = \sqrt {{\mathop{\rm E}\nolimits} \left[ {{{(N + {N_1})}^2}} \right] - {{\mathop{\rm E}\nolimits} ^2}(N + {N_1})} \\ \;\;\;\;\;\;\;\; = \sqrt {{\mathop{\rm E}\nolimits} ({N^2}) + {\mathop{\rm E}\nolimits} ({N_1}^2) - {{\mathop{\rm E}\nolimits} ^2}(N) - {{\mathop{\rm E}\nolimits} ^2}({N_1})} \\ \;\;\;\;\;\;\;\; = \sqrt {[{\mathop{\rm E}\nolimits} ({N^2}){\rm{ - }}{{\mathop{\rm E}\nolimits} ^2}(N)] + [{\mathop{\rm E}\nolimits} ({N_1}^2) - {{\mathop{\rm E}\nolimits} ^2}({N_1})]} \\ \;\;\;\;\;\;\;\; = \sqrt {{\sigma ^2} + {\sigma _1}^2} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \end{array}$ (9)

由式(9) ,可以得到

${{\mathop{\rm P}\nolimits} _{1L}} = k\sqrt {{\sigma ^2} + {\sigma _1}^2} + b$ (10)

同理易知, 若对含噪图像加入与N不相关、强度为σ2的加性高斯白噪声N2,可得

${{\mathop{\rm P}\nolimits} _{2L}} = k\sqrt {{\sigma ^2} + {\sigma _2}^2} + b$ (11)

将式(8) 、(10) 、(11) 联立,使用Matlab解方程组可以得到σ的估计值$\mathop \sigma \limits^ \wedge $

$\begin{array}{l} \mathop \sigma \limits^ \wedge = \frac{{( - {{\mathop{\rm P}\nolimits} _L}^2 + {\sigma _1}^2{{\mathop{\rm P}\nolimits} _L}^2 + 2{\sigma _2}^2{{\mathop{\rm P}\nolimits} _L}{{\mathop{\rm P}\nolimits} _{L1}} - 2{\sigma _1}^2{{\mathop{\rm P}\nolimits} _L}{{\mathop{\rm P}\nolimits} _{L2}} - {\sigma _2}^2{{\mathop{\rm P}\nolimits} _{L1}}^2 + {\sigma _1}^2{{\mathop{\rm P}\nolimits} _{L2}}^2)}}{{\sqrt {({\sigma _2}^2{{\mathop{\rm P}\nolimits} _L}^2 + {\sigma _1}^2{{\mathop{\rm P}\nolimits} _{L2}}^2 - {\sigma _2}^2{{\mathop{\rm P}\nolimits} _{L1}})} }} \times \\ \frac{1}{{2\sqrt {({{\mathop{\rm P}\nolimits} _L}{{\mathop{\rm P}\nolimits} _{L2}}^2 - {{\mathop{\rm P}\nolimits} _{L1}}{{\mathop{\rm P}\nolimits} _{L2}}^2 + {{\mathop{\rm P}\nolimits} _{L1}}{P_L}^2 + {{\mathop{\rm P}\nolimits} _{L1}}^2{{\mathop{\rm P}\nolimits} _{L2}} - {{\mathop{\rm P}\nolimits} _{L2}}{{\mathop{\rm P}\nolimits} _L}^2 - {{\mathop{\rm P}\nolimits} _L}{{\mathop{\rm P}\nolimits} _{L1}}^2)} }} \end{array}$ (12)

至此,得到了对含噪图像X进行噪声估计的方法,具体步骤如下:

1) 对含噪图像X进行SVD;

2) 确定合适的L取值,计算得到X的奇异值矩阵尾部L个奇异值的均值PL

3) 对X加入标准差为σ1、均值为零的AWGN,得到图像X1,再对X1进行SVD并计算得到X1的奇异值矩阵尾部L个奇异值的均值P1L

4) 对X加入标准差为σ2,均值为零的AWGN,得到图像X2;对X2进行SVD,并计算得到X2的奇异值矩阵尾部L个奇异值的均值P2L

5) 将上述计算结果代入式(12) ,计算得到X的噪声标准差的估计值$\mathop \sigma \limits^ \wedge $

1.2.2 具有自适应噪声估计特性的K-SVD去噪算法

结合1.2.1节提出的图像噪声估计方法,本节提出一种具有自适应噪声估计特性的K-SVD去噪算法。为便于后文描述,该方法简记为自适应KSVD去噪算法(Adaptive K-SVD, AKSVD),用于求解式(1) 中的KSVD去噪模型。AKSVD算法描述如下:

输入 含噪图像X,图像块的维数n,字典的原子数k,迭代次数J,噪声因子C。

输出 去噪后图像$\mathop {\bf{X}}\limits^ \wedge $

过程如下:

步骤1 初始化。令$\mathop {\bf{X}}\limits^ \wedge = {\bf{X}}$,初始字典$\mathop {\bf{D}}\limits^ \wedge $为DCT字典。

步骤2 噪声估计。利用1.2.1 节提出的噪声估计算法对进行噪声估计,得到预估噪声标准差σ,计算得到图像保真度λ=30/σ。

步骤3 迭代。先求解图像块RijX对应的稀疏系数$\overset{\wedge }{\mathop{{{\mathbf{\alpha }}_{\text{ij}}}}}\,$,再根据由$\overset{\wedge }{\mathop{{{\mathbf{\alpha }}_{\text{ij}}}}}\,$组成的稀疏矩阵$\mathop {\bf{\alpha }}\limits^ \wedge $更新稀疏字典$\mathop {\bf{D}}\limits^ \wedge $。即每次迭代包含2个步骤,迭代持续J次:

1) 固定字典${\hat{D}}$,通过OMP算法[14]求解$\mathop {\bf{\alpha }}\limits^ \wedge \forall i,j,\;{\mathop {\bf{\alpha }}\limits^ \wedge _{ ij}} = \mathop {\arg \min }\limits_{{{\bf{\alpha }}_{ij}}} {\left\| {{{\bf{\alpha }}_{ij}}} \right\|_0},\;(\left\| {{{\bf{R}}_{ij}}\mathop {\bf{x}}\limits^ \wedge - \mathop {\bf{D}}\limits^ \wedge {{\bf{\alpha }}_{ij}}} \right\|_2^2 \le n(C{\sigma ^2}))$

2) 固定$\mathop {{{\bf{\alpha }}_{{\rm{ij}}}}}\limits^ \wedge $,逐个对$\mathop {\bf{D}}\limits^ \wedge $中的原子${{\bf{d}}_m}(m = 1,2, \cdots k)$进行更新:

① 找出满足${\omega _m} = \{ [i,j]|{\mathop {\bf{\alpha }}\limits^ \wedge _i}_j(m) \ne 0\} $ 的图像块;

② 对每个满足[i, j]∈ωm的图像块,计算对应的残差 ${\bf{e}}_{ij}^m = {{\bf{R}}_{ij}}{{\bf{X}}_{ij}} - \sum\limits_{l \ne m} {{{\bf{d}}_m}{{\bf{\alpha }}_{ij}}(l)} $,并得到对应的残差矩阵Em

③ 对Em进行SVD,得到更新的$\mathop {{{\bf{d}}_m}}\limits^ \wedge $$\overset{\wedge }{\mathop{{{\mathbf{\alpha }}_{\text{ij}}}}}\,$

步骤4 代入求解。将之前得到的$\lambda ,\mathop {{{\bf{\alpha }}_{ij}}}\limits^ \wedge \mathop {\bf{D}}\limits^ \wedge $代入式(13) ,求得$\mathop {\rm{X}}\limits^ \wedge $

$\mathop {\bf{X}}\limits^ \wedge = {(\lambda {\bf{I}} + \sum\limits_{ij} {{\bf{R}}_{ij}^T{{\bf{R}}_{ij}}} )^{ - 1}}(\lambda {\bf{Y}} + \sum\limits_{ij} {{\bf{R}}_{ij}^T{\bf{D}}{{\mathop {\bf{\alpha }}\limits^ \wedge }_{ij}}} )$ (13)
2 结合PCA和AKSVD的高光谱图像去噪方法 2.1 主成分分析

主成分分析(PCA)是一种常用的数据降维方法,主要是通过对数据的协方差矩阵进行特征值分解,特征值分解后得到的特征向量对应数据的主成分,相应的特征值对应数据在各个主成分上的权重。对高光谱图像进行PCA变换后,大部分信息集中在前几个主成分中,可以通过保留前几个主成分达到去噪的目的。

设原始数据X=(x1, x2, …, xm),其中xi是N维矢量,则X的PCA变换为:

${\bf{Y}} = {{\bf{U}}^T}({\bf{X}} - \overline {\bf{X}} )$ (14)

式中,$\overline {\bf{X}} = (\overline {{{\bf{x}}_{\rm{1}}}} ,\overline {{{\bf{x}}_{\rm{2}}}} , \cdots \overline {{{\bf{x}}_m}} ), = \frac{1}{N}\sum\limits_{j = 1}^N {{{\bf{x}}_{ij}}} $为各矢量均值;而U=(u1, u2, …, ud),是由X的协方差矩阵$\sum {_{\bf{x}}} = {\mathop{\rm E}\nolimits} \{ ({\bf{X}} - \overline {\bf{X}} ){({\bf{X}} - \overline {\bf{X}} )^T}\} $ 求得的按降序排列的特征值λi对应的特征向量ui组成的。d为所取主成分的个数。

对应的PCA逆变换为:

${{\bf{X}}^{\rm{*}}} = {\bf{UY}}$ (15)

式中,X*为使用X的前d个主成分分量进行重构得到的重构数据。对于含噪数据X,舍弃后(N-d)个主成分分量能够起到一定的去噪效果。

2.2 结合PCA和AKSVD的高光谱图像去噪方法

如前文所述,本文所提出去噪方法的基本原理是将PCA方法对高维数据的关键信息提取能力和AKSVD算法对单张图像的噪声估计和去噪能力结合起来,对高光谱图像进行去噪处理。本文将PCA和AKSVD算法结合起来,设计出一种进的高光谱图像去噪算法,设包含N个波段的含噪高光谱数据为X,算法流程如图 1所示。

图 1 本文去噪算法的流程

1) 对X进行PCA操作得到其对应的PCA变换矩阵TPCA,提供备选的主成分分量以便后文对X中的单幅波段图像(记为Xj(1≤j<N))进行预去噪,且将TPCA的N个成分分量所对应的能量比重记为EIV(j)(1≤j≤N)。

2) 根据EIV,对TPCA进行分类。传统的PCA去噪方法是选取几个较大的主成分分量用于PCA变换去噪,但该做法会导致图像损失较多的细节信息。本文根据相邻成分对应的EIV之比对TPCA进行分组处理,将TPCA分为2部分,分别为:清晰成分分量(记为TPCA(high))和含噪成分分量(记为TPCA(low))。经多次对比实验发现,当1≤i<N且相邻成分分量的EIV临界比值取EIV(i)/EIV(i+1) =1.1时,清晰图像和含噪图像的区分度最好(即利用得到的TPCA(high)对图像进行PCA变换和PCA逆变换后得到的图像信噪比达到最大)。因此,分类判别准则如式(16) 所示:

$\left\{ \begin{array}{l} {{\bf{T}}_{PCA(high)}}(i) = {{\bf{T}}_{PCA}}(i),\;\;\;\;\;\;if\;\frac{{EIV(i)}}{{EIV(i + 1)}} > 1.1\\ {{\bf{T}}_{PCA(low)}}(i - k) = {{\bf{T}}_{{\rm{PCA}}}}(i),\;if\;\frac{{EIV(i)}}{{EIV(i + 1)}} \le 1.1 \end{array} \right.$ (16)

式中:TPCA(high)的成分数量为k;TPCA(low)的成分数量为N-k。

3) 分别对Xj(1≤j≤N))逐一使用TPCA(high)和TPCA(low)进行PCA变换和PCA逆变换(IPCA),则Xj能够得到k张图像组成的清晰图像组Xj(k)和N-k张图像组成的含噪图像组Xj(N-k)

4) 使用AKSVD算法对Xj(N-k)(1≤j≤N)进行去噪处理。

5) 根据原本各成分分量对应的EIV,将Xj对应的清晰图像组Xj(k)和去噪后的Xj(N-k)进行加权融合得到最终去噪图像Xj ′,最后将N幅去噪后的波段图像依次排列即可得到完成去噪的高光谱数据X′。

3 仿真实验与结果分析

为验证本文提出的去噪方法的有效性,本章选择将本文所提出的去噪方法和以下3种去噪方法进行对比:1) 文献[6]提出的直接利用PCA去噪方法(简记为PCA方法);2) 文献[8]提出的结合PCA和二维的二元小波阈值收缩的去噪方法(简记为PCA-Bish方法);3) 文献[10]提出的基于PCA和Contourlet变换的去噪方法(简记为PCA-Contourlet方法)。为充分检验以上4种方法对高光谱遥感图像的去噪效果,针对模拟高光谱数据和真实高光谱数据,本文设计了2组对比实验。为便于对比,本文进行KSVD去噪时所选原子块的大小均为4×4。

3.1 模拟数据对比实验

本节选用的模拟高光谱图像数据由AVIRIS sensor获取的真实地物数据集Indian Pines参考美国地质调查局数字光谱实验室所提供的光谱数据库“splib06a”合成得到。Indian Pines数据集涵盖了16类确定类别的地物,本文从光谱数据库“splib06a”中找出与这16类地物及未标记类别的地物近似的光谱数据,按照如下方式对无噪的模拟高光谱图像进行构建[15]:首先,分别对每一类地物按照选定的近似光谱数据赋值;然后,将整个数据的数值范围线性拉伸到[0,255]后截取前64个波段,得到大小为145×145×64的模拟数据。

假定噪声为AGWN,为充分衡量4种方法对模拟高光谱数据的去噪效果,本文分别考虑了3种含噪情况:1) 随机选取32个波段加噪且各波段噪声强度相同,噪声标准差分别为σ={10, 15, 20};2) 对所有波段加噪且噪声强度在[3, 12]内呈均匀分布;3) 对所有波段加噪且每个波段噪声强度取为该波段最大像素值的8%。为便于去噪性能对比,本节所用4种方法在进行PCA变换时均保留前80%的主成分分量。

除视觉效果外,本文还使用信噪比(Signal-to-Noise Ratio, SNR)、峰值信噪比(Peak Signal-to-Noise Ratio, PSNR)和能较好反映人眼视觉感受的图像相似度(Structural Similarity Index Measurement, SSIM)三个指标对不同方法的性能进行比较。对于大小为M×N的高光谱图像,其中SNR、PSNR和SSIM指标定义分别为:

$SNR = 10lg\frac{{\sum\limits_{i,j,p} {{\bf{X}}{{(i,j)}^2}} }}{{\sum\limits_{i,j} {{{\left( {\overline {\bf{X}} (i,j) - {\bf{X}}(i,j)} \right)}^2}} }}$ (17)
$PSNR = 10lg\left[ {\frac{{{{255}^2}}}{{\sum\limits_{i,j} {\frac{{{{(\overline {\bf{X}} (i,j) - {\bf{X}}(i,j))}^2}}}{{M \times N}}} }}} \right]$ (18)
$SSIM = \frac{{(2{u_{\bf{X}}}{u_{\overline {\bf{X}} }} + {C_1})(2{\sigma _{\bf{X}}}{\sigma _{\overline {\bf{X}} }} + {C_2})}}{{({u_{\bf{X}}}^2 + {u_{\overline {\bf{X}} }}^2 + {C_1})({\sigma _{\bf{X}}}^2 + {\sigma _{\overline {\bf{X}} }}^2 + {C_2})}}$ (19)

式中:i=1, 2, …, M; j=1, 2, …, N;X表示不含噪声的参考数据;X表示含噪声的数据或者去噪后的数据;uX表示图像X灰度的均值;σX表示图像X灰度的标准差;C1=(255K1)2和C2=(255K2)2为常数,;本文中取K1=K2=0.01。

图 2中显示了对模拟数据Indian Pines、含噪声数据(σ=20) 和4种方法去噪后数据中随机选取的2个波段图像。图中第1、2行图像分别对应模拟数据的第11和第44波段图像。从第一行图像可以看出,图 2(c)列对应的PCA方法去噪后的图像比较模糊,边缘细节不清晰;相比PCA方法,图 2(d)图 2(e)列对应的PCA-Bish和PCA-Contourlet方法去噪后图像的边缘信息得到了较好保留,但图像的细节信息(右侧中部色块),尤其是对比度方面损失较为严重,这主要是由于Bish小波阈值收缩方法和Contourlet方法去噪时使用固定的稀疏基,由于经过PCA变换产生的含噪图像组的图像特征比较复杂,固定的稀疏基难以准确地对图像进行稀疏表示,造成去噪过程中图像失真较为严重,以致融合后的图像细节信息丢失。由图 2(f)列对应的图像可以看出,相比前3种方法,本文提出的PCA-AKSVD方法既能较好地去除噪声影响,同时尽可能多地保留了原图像的细节信息。

由于高光谱图像波段众多,仅凭某个波段图像的视觉效果评定,并不足以说明方法的有效性,为此需要进行全部波段的数据分析。表 1统计了在第1) 种噪声情况下含噪高光谱图像及不同方法去噪后的高光谱图像的三项评价指标的平均值和标准差。通过对表中数据的观察,可以看出SNR、PSNR和SSIM三者的变化趋势基本一致,由此对4种方法的去噪效果进行分析:1) 在噪声的标准差很低时,本文方法并无明显优势。这主要是因为在噪声很小时PCA变换产生的含噪图像组保留了绝大多数图像信息(即其中所含的噪声几乎可以忽略),除PCA外的三种方法对其进行去噪处理时会造成图像细节信息的丢失,而PCA方法直接使用主成分图像,不存在这一问题。2) 当噪声的标准差逐渐增大时,本文方法去噪性能的优势越发明显。例如,当σ=15时,本文方法的SNR和PSNR的平均值已明显优于其他三种方法;而当σ=20时,本文方法的各项指标均达到最优。如本文方法的PSNR的平均值和标准差(Standard deviation, Std)分别优于排名第二的PCA-Contourlet方法达到1.5426 dB和43.6%。由此可见:在高光谱数据各波段受到较为严重的同等强度加性噪声污染时,本文方法的去噪能力较其他三种方法有着明显优势。

图 2 不同算法对模拟数据Indian Pines第11和第44波段去噪后图像与原图像
表 1 4种去噪方法对模拟数据Indian Pines的去噪性能比较(各波段噪声强度相同)

针对噪声强度随波段变化的第2) 和3) 种噪声情况,表 2展示了在第2) 和第3) 种噪声情况下不同方法去噪后的高光谱图像的三项评价指标的平均值和标准差。对于噪声强度随波段变化且覆盖全部波段的噪声情况2) 和3) ,观察表 2中数据,可以得到:a)在SNR和PSNR指标上,4种方法中按平均值和标准差排名去噪效果从好到差依次是:PCA-AKSVD、PCA-Contourlet、PCA-Bish、PCA。其中,针对噪声情况2) ,本文方法去噪得到的高光谱各波段图像的SNR平均值和标准差分别优于次优的PCA-Contourlet方法达到2.0132 dB和55.4%。b)在SSIM指标上,4种方法中按标准差排名去噪效果从好到差依次是:PCA-AKSVD、PCA-Contourlet、PCA、PCA-Bish。其中,针对噪声情况3) ,本文方法去噪后得到的高光谱各波段图像的SSIM平均值分别比PCA和PCA-Bish方法高出了1.55%和1.19%。

由于噪声强度随波长变化,为进一步比较不同方法在每个波段上的去噪性能,图 3显示了分别在第2) 和第3) 种噪声情况下不同去噪方法去噪后数据在各个波段上的PSNR和SSIM之差。为便于观察,图 3中还增加了一条黑色水平线,表示4种方法对各波段去噪后评价指标的平均值。从图 3(a)可以看出,4种方法去噪得到的高光谱各波段图像的PSNR变化趋势非常稳定,从高到低依次是:PCA-AKSVD、PCA-Contourlet、PCA-Bish、PCA,且只有本文方法的PSNR值始终高于平均值。观察图 3(b),可以看出:PCA方法去噪的性能偏低,在绝大多数波段SSIM值低于4种方法在该波段的平均值;PCA-Bish方法较PCA方法效果有所提升,部分波段达到甚至超过了4种方法在该波段达到的平均值(如10~20波段)。PCA-Contourlet方法在大多数波段都取得了较好的去噪效果,整体比较稳定且大多数波段优于前两种方法。而本文方法去噪效果较PCA-Contourlet方法又有所提升,尤其在前20个波段提升效果显著。

综合来看,针对含加性噪声的模拟高光谱图像去噪时,对比三种同类方法,本文方法不仅在图像视觉效果方面更好,而且对绝大多数波段的去噪效果和稳定性上有着明显优势。

表 2 4种去噪方法对模拟数据Indian Pines的去噪性能比较(各波段噪声强度随波段变化)
图 3 4种不同去噪方法对模拟数据Indian_Pines去噪后结果比较
3.2 真实数据对比实验

本节选用的高光谱图像数据参考由AVIRIS sensor获取的KSC(Kennedy Space Center)图像。KSC数据是由AVIRIS sensor在美国佛罗里达州的肯尼迪航天中心采集,光谱分辨率约为10 nm,包含176个波段。本文截取大小为256×256×128的高光谱数据用于验证不同方法对含噪真实高光谱数据的去噪性能。去噪前将所有数据的数值范围线性拉伸到[0,255], 该数据前20个波段中含有较强的噪声。为便于去噪性能对比,4种方法在选择主成分分量时均保留前70%的主成分分量。

图 4显示了对真实高光谱数据KSC使用4种方法去噪前后从数据中随机选取的某两个波段图像。

图 4 实际数据KSC经不同算法去噪前后的第3和第14波段图像

图 4可以看出,PCA方法去噪后的图像的右下和中上部出现了明显的失真现象,在4种方法中效果最差。而其余三种方法去噪后的图像的纹理细节均比较清晰,靠目视很难分辨优劣。为更好地评价4种方法的去噪效果,除去前文使用的客观评价指标,本节选用了一个更能反应图像视觉效果的指标:图像清晰度(DeFinition,DF),对于尺寸为M×N的图像定义为:

$DF = \sum\limits_{i = 1}^M {\sum\limits_{j = 1}^N {[{\bf{X}}(i + 2,j) - {\bf{X}}(i,j)]} } $ (20)

其中,X(i, j)表示图像X对应像素点(i, j)的灰度值。

表 3图 5分别显示了采用4种不同方法去噪后KSC数据各波段图像DF的统计数据和DF差值的变化曲线。结合表 3中数据和图 5,可以看出:整体来看,PCA方法的去噪效果最差,其DF值的平均值和标准差均明显次于其他三种方法;在数据的前100波段之内,PCA方法得到的DF值始终最低。PCA-Bish方法和PCA-Contourlet方法得到的DF值的平均值和标准差非常接近,因此在整体性能上相差不大, 但PCA-Bish方法的稳定性明显不足。例如,在前80个波段中,PCA-Bish方法得到的DF值为最优,但在80波段之后突然下滑,尤其是在100波段之后变为了最低值。相比之下,PCA-Contourlet方法得到的DF值虽然始终没能达到最优,但稳定性明显胜过PCA-Bish方法。与PCA-Contourlet方法相比,本文方法得到的DF值的平均值和标准差的优势分别为2.64%和48.1%;除去在前80个波段仅次于PCA-Bish方法,在其余绝大多数波段,本文方法的DF值均为最优值。同时,本文方法得到的DF值曲线是4种方法中唯一始终位于平均值曲线上方的。由以上分析不难得出:对真实高光谱数据KSC去噪时,本文方法在绝大多数波段的去噪效果均优于其余三种方法,且具有很强的稳定性。

算法的实现是通过Matlab 7.0平台编程实现的,所用系统为Window 7,硬件条件为Core i5-2540M@2.5 GHz、1.5 GB内存。表 4列出了4种方法运算时间的比较结果。

表 3 4种方法对真实数据KSC去噪后各波段图像清晰度的结果
图 5 不同算法去噪后KSC图像1~128波段数据的结果
表 4 不同算法对KSC图像数据进行去噪的运算时间

本文方法运算时间要长于其余3种方法,这主要是因为,对尺寸为M×N的图像,AKSVD算法要对其进行3次SVD以及(M×N)/(4×4) 次的KSVD稀疏分解,这一过程会增加一定的时间开销。但结合去噪结果来看,本文方法可以视为是用时间代价换取了去噪效果的提升,而且通过对AKSVD方法所用原子块尺寸的调整,可以实现时间开销和去噪效果之间的有效置换,这一特点也提高了本文方法在高光谱遥感图像处理的工程应用中的适用性。

4 结语

由于高光谱图像各波段噪声强度分布的复杂性,传统的特征降维类去噪方法并不能满足高光谱遥感数据的去噪要求。本文先将SVD的噪声估计能力以及KSVD字典学习去噪算法的消噪能力结合起来,提出了一种具有自适应噪声估计特性的字典学习去噪算法AKSVD,再将AKSVD算法与特征降维方法PCA所具备的关键信息提取能力结合起来,提出了一种具备自适应噪声估计特性的高光谱图像去噪方法PCA-AKSVD方法。

与传统的去噪方法不同,本文方法不只利用高光谱图像整体数据信息,还通过噪声估计和字典学习实现了对单幅图像结构信息的充分利用。通过对各类对比实验表明:与同类方法相比,本文方法在消除噪声的同时,也最大限度保留了原始数据的细节信息,显著提高了去噪后高光谱图像的画面质量。稍显不足的是,本文是从去噪效果的角度对高光谱遥感图像去噪方法进行设计,而在时间开销上则考虑得不足,如何降低去噪方法的时间复杂度将在今后的研究中进一步探讨和研究。

参考文献
[1] MANOLAKIS D, TRUSLOW E, PIEPER M, et al. Detection algorithms in hyperspectral imaging systems:an overview of practical algorithms[J]. IEEE Signal Processing Magazine, 2014, 31 (1) : 24-33. doi: 10.1109/MSP.2013.2278915
[2] GUERRA R, SANTOS L, LOPEZ S, et al. A new fast algorithm for linearly unmixing hyperspectral images[J]. IEEE Transactions on Geoscience & Remote Sensing, 2015, 53 (12) : 6752-6765.
[3] LI H, XIAO G, XIA T, et al. Hyperspectral image classification using functional data analysis[J]. IEEE Transactions on Cybernetics, 2014, 44 (9) : 1544-1555. doi: 10.1109/TCYB.2013.2289331
[4] YUAN Q, ZHANG L, SHEN H. Hyperspectral image denoising employing a spectral-spatial adaptive total variation model[J]. IEEE Transactions on Geoscience & Remote Sensing, 2012, 50 (10) : 3660-3677.
[5] RASTI B, SVEINSSON J R, ULFARSSON M O. Design of personal rapid transit circulator for major activity center:Hacienda Business Park, Pleasanton, California[J]. Transportation Research Record Journal of the Transportation Research Board, 2007, 43 (1) : 104-113.
[6] STEPHAN K, HIBBITTS C A, HOFMANN H, et al. Reduction of instrument-dependent noise in hyperspectral image data using the principal component analysis:applications to Galileo NIMS data[J]. Planetary & Space Science, 2008, 56 (3/4) : 406-419.
[7] GREEN A A, BERMAN M, SWITZER P, et al. A transformation for ordering multispectral data in terms of image quality with implications for noise removal[J]. IEEE Transactions on Geoscience & Remote Sensing, 1988, 26 (1) : 65-74.
[8] CHEN G, QIAN S E. Denoising of hyperspectral imagery using principal component analysis and wavelet shrinkage[J]. IEEE Transactions on Geoscience & Remote Sensing, 2011, 49 (3) : 973-980.
[9] CHEN G, BUI T D, QUACH K G, et al. Denoising hyperspectral imagery using principal component analysis and block-matching 4D filtering[J]. Canadian Journal of Remote Sensing, 2014, 40 (1) : 60-66. doi: 10.1080/07038992.2014.917582
[10] 常威威, 郭雷, 刘坤, 等. 基于Contourlet变换和主成分分析的高光谱数据噪声消除方法[J]. 电子与信息学报, 2009, 31 (12) : 2892-2896. ( CHANG W W, GUO L, LIU K, et al. Denoising of hyperspectral data based on Contourlet transform and principal component analysis[J]. Journal of Electronics & Information Technology, 2009, 31 (12) : 2892-2896. )
[11] AHARON M, ELAD M, BRUCKSTEIN A K. K-SVD:an algorithm for designing overcomplete dictionaries for sparse representation[J]. IEEE Transactions on Signal Processing, 2006, 54 (11) : 4311-4322. doi: 10.1109/TSP.2006.881199
[12] ELAD M, AHARON M. Image denoising via sparse and redundant representations over learned dictionaries[J]. IEEE Transactions on Image Processing, 2006, 15 (12) : 3736-3745. doi: 10.1109/TIP.2006.881969
[13] LIU W, LIN W S. Additive white Gaussian noise level estimation in SVD domain for images[J]. IEEE Transactions on Image Processing, 2013, 22 (3) : 872-883. doi: 10.1109/TIP.2012.2219544
[14] TROPP J A, GILBERT A C. Signal recovery from random measurements via orthogonal matching pursuit[J]. IEEE Transactions on Information Theory, 2007, 53 (12) : 4655-4666. doi: 10.1109/TIT.2007.909108
[15] 霍雷刚, 冯象初. 基于主成分分析和字典学习的高光谱遥感图像去噪方法[J]. 电子与信息学报, 2014, 36 (11) : 2723-2729. ( HUO L G, FENG X C. Denoising of hyperspectral remote sensing image based on principal component analysis and dictionary learning[J]. Journal of Electronics & Information Technology, 2014, 36 (11) : 2723-2729. )