CN108765313A - 基于类内低秩结构表示的高光谱图像去噪方法 - Google Patents

Abstract

本发明提供了一种基于类内低秩结构表示的高光谱图像去噪方法。首先，采用光谱聚类算法将高光谱数据划分成多个类，实现对高光谱图像类内低秩特性的有效利用。然后，利用奇异值分解方法对每一类的类内低秩矩阵进行分解，并使用结构稀疏方法对奇异值进行建模。最后，利用优化方法自适应地进行模型求解。本发明方法有效利用了高光谱图像的局部相似性和非局部相似性，不需要输入噪声方差或其它先验信息，且可以适用于各种不同噪声情况，具有较好的去噪效果和良好的适应性。

Description

基于类内低秩结构表示的高光谱图像去噪方法

技术领域

本发明属高光谱图像的处理技术领域，具体涉及一种基于类内低秩结构表示的高光谱图像去噪方法。

背景技术

高光谱图像是一个三维数据块，丰富的空-谱信息使其被广泛应用于资源勘探、环境监测、生物医学等领域。然而，高光谱成像受大气、太阳角度、传感器混合效应等的影响，图像中不可避免地引进噪声，因此，去噪是高光谱图像处理领域的一个研究热点。如文献“Maggioni M,Katkovnik V,Egiazarian K,et al.Nonlocal Transform-Domain Filterfor Volumetric Data Denoising and Reconstruction[J].IEEE Transactions onImage Processing A Publication of the IEEE Signal Processing Society,2012,22(1):119-133.”针对高光谱图像这类三维数据，提出了四维块匹配算法。该算法将图像数据划分为许多有重叠区域的三维小图像块，并把相似图块作为一组进行处理。该算法同时利用了存在于图块内部的局部相似性和图块之间的非局部相似性。在变换域上，每组数据都具有高度的稀疏性，使得该算法能完成对图像信号和噪声的有效分离。四维块匹配算法需要图像噪声方差的先验知识，但在实际情况下，噪声往往是未知的。因此，需要人为地估计噪声方差作为算法输入，但由于噪声的多样性和估计的不准确性，该算法对噪声的适应性较差。

发明内容

为了克服现有技术的不足，并且进一步提高高光谱图像去噪方法对不同噪声的适应性，本发明提出了一种基于类内低秩结构表示的高光谱图像去噪方法。首先，采用光谱聚类算法将高光谱数据划分成多个类，实现对高光谱图像类内低秩特性的有效利用。然后，利用奇异值分解方法(SVD)对每一类的类内内低秩矩阵进行分解，并使用结构稀疏方法对奇异值进行建模。最后，利用优化方法自适应地进行模型求解。

一种基于类内低秩结构表示的高光谱图像去噪方法，其特征在于包括以下步骤：

步骤一：首先，将给定的三维高光谱图像数据转化成二维矩阵R，矩阵大小为n_b×n_P，其中，n_b为光谱波段个数，n_p为像元的个数；然后，采用K-means方法对矩阵R进行聚类处理，得到K个类Y₁,…,Y_K，K的取值范围为30～90，每一类Y_k所包含的像元数为n_k，k＝1,…,K。

步骤二：建立去噪模型如下：

其中，λ为噪声方差，||·||_F表示Frobenius范数，X＝[X₁,…,X_K]为去噪后图像，p(X_k|Θ_lk)为X_k的低秩先验，Θ_lk表示低秩先验p(X_k|Θ_lk)中的参数。利用奇异值分解的方法将低秩矩阵X_k表示为将公式(1)转化为：

其中，s_k为由S_k对角线上的奇异值组成的向量，即奇异值的噪声观测项；p(s_k|Θ_sk)表示奇异值向量s_k的稀疏先验，Θ_sk表示稀疏先验p(s_k|Θ_sk)中的参数；U_k和V_k通过奇异值分解Y_k得到。

步骤三：首先，引入辅助变量t_k，将公式(2)转化为：

其中，β为预定义的标量参数，β＝0.1；p(t_k|Θ_tk)表示辅助变量t_k的稀疏先验，Θ_tk表示稀疏先验p(t_k|Θ_tk)中的参数，||·||₂表示L₂范数。

然后，采用半二次分裂算法对公式(3)进行求解，得到t_k，k＝1,…,K，并得到去噪后的图像X，具体为：

步骤1：稀疏学习秩估计，固定{s_k}，{t_k}的估计表示为：

令s_k＝t_k+ε_k，ε_k是噪声腐蚀项，ε_k服从均值为0的高斯分布，s_k的似然函数表示为∝表示正比于，diag(·)表示matlab中的diag函数；使用重加权拉普拉斯先验对t_k中的结构信息进行表示，令t_k服从高斯分布，t_k的概率密度函数表示为γk为高斯分布中的方差向量，Γ_k＝diag(γ_k)；令γ_k服从伽马分布，γ_k的概率密度函数表示为κ_k为伽马分布中的形状参数向量，R_k＝min(n_b,n_k)；使用基于隐变量的贝叶斯稀疏学习方法对参数ε_k、γ_k和κ_k进行学习，得到如下稀疏先验：

其中，表示加权迹范数， 公式(5)转化为一个统一的变分框架，即：

利用交替最小化算法对公式(6)进行优化求解，具体步骤如下：

步骤a：初始化，记初始迭代次数为0，分别令ε_k、γ_k、κ_k为值为1的R_k维向量，即ε_k＝γ_k＝κ_k＝[1,...,1]^T，s_k初始化为对Yk进行奇异值分解得到的奇异值组成的向量；

步骤b：按照t_k＝Γ_k(∑_bk)^-1s_k更新t_k，其中Γ_k＝diag(γ_k)；

步骤c：按照更新ε_k，其中，Q＝t_k-sk，

步骤d：按照更新γ_k，其中，

步骤e：根据更新参数κ_k；

步骤f：迭代次数加1，如果迭代次数大于15或t_k的更新差异小于10^-6，则停止，此时得到最终的t_k；否则，回到步骤b。

步骤(2)：低秩矩阵修复，固定{t_k}，{s_k}的估计表示为：

根据将对{s_k}的求解转化为以下等价问题：

其中，T_k为对角矩阵，t_k为T_k对角线元素组成的向量。公式(8)的闭式解为：

其中，按公式(9)求X_k，k＝1,...,K，即获得去噪后图像X＝[X₁,…,X_K]。

本发明的有益效果是：由于采用K-means方法对原始图像进行聚类，实现了对高光谱图像局部相似性和非局部相似性的有效利用；利用奇异值分解方法对每一类的类内低秩矩阵进行分解，并使用重加权拉普拉斯先验对奇异值进行稀疏表示，进而将求解低秩矩阵秩的问题转化为求解矩阵奇异值稀疏性的问题。本发明方法不需要输入噪声方差或其它先验信息，且可以适用于各种不同噪声情况，具有较好的去噪效果和良好的适应性。

具体实施方式

本发明提供了一种基于类内低秩结构表示的高光谱图像去噪方法，其具体过程如下：

1、高光谱图像聚类

由于高光谱图像数据具有局部和非局部空间的相似性，所以可以采用K-means聚类方法对其进行聚类，划分成多个类。首先，对于给定的高光谱图像数据包括n_r行、n_c列和n_b个波段。为计算方便，先将三维高光谱图像数据转化成二维矩阵其中，n_p＝n_r×n_c为像元的个数；其次，采用K-means方法对矩阵R进行聚类处理，得到K个类Y₁,…,Y_K，每一类即为R中的一个同质区域，n_k为每一个类中的像素数，k＝1,...,K，K的取值范围为30～90，本实施例中K＝90，高光谱图像R可以重新表示为Y＝[Y₁,…,Y_K]。

2、建立去噪模型

聚类得到的每一类不仅存在明显的光谱间相似性，而且还存在光谱内相似性，因此每一类都存在低秩结构。对于第k类Y_k，其含噪观测模型可以低秩表示为：

Y_k＝X_k+N_k (10)

其中，X_k为第k类低秩矩阵，即Y_k的背景图像，N_k为第k类噪声项。当聚类数K为1时，上式即退化为对整个图像的低秩表示，即Y＝X+N，其中Y表示原始高光谱图像即含噪图像，X为对应背景图像的低秩矩阵，即待求的去噪后图像，N为噪声项。假定噪声服从零均值矩阵正态分布∑_n＝λI为噪声方差矩阵，λ表示噪声方差，指示噪声强度，I为单位矩阵。根据公式(10)的观测模型，得到如下的似然估计：

其中，p(Y|X,λ)表示Y的似然函数，∝表示正比于，||·||_F表示Frobenius范数。

给定噪声方差λ和X_k的低秩先验p(X_k|Θ_lk)，Θ_lk为低秩先验p(X_k|Θ_lk)中的参数。根据公式(11)，对低秩矩阵X的估计可以转化为最大后验估计：

其中，p(X|Y)表示X的后验概率。公式(12)等价于：

通过奇异值分解(SVD)的方法对低秩矩阵进行表示，可以将求解低秩矩阵秩的问题转化为求解矩阵奇异值稀疏性的问题。因此，根据SVD方法，X_k可以表示为公式(13)可以重新表示为：

其中，表示S_k对角线上奇异值组成的向量，即奇异值的噪声观测项，R_k＝min(n_b,n_k)。p(s_k|Θ_sk)表示奇异值向量s_k的稀疏先验，Θ_sk表示稀疏先验p(s_k|Θ_sk)中的参数；为简单起见，U_k和V_k通过SVD分解Y_k得到。公式(13)中低秩矩阵修复问题转化为公式(14)中的寻找合适的稀疏表示s_k的问题。

3、模型求解

对于一个复杂的p(s_k|Θ_lk)(例如混合高斯先验)，公式(14)将很难求解，因此，首先，引入辅助变量t_k，将公式(14)改写为：

其中，p(t_k|Θ_tk)表示t_k的稀疏先验，Θ_tk表示稀疏先验p(t_k|Θ_tk)中的参数，β是预定义的标量参数，β＝0.1。

然后，采用半二次分裂算法(half-quadratic splitting scheme)将公式(15)划分成两个简单的子问题(稀疏学习秩估计、低秩矩阵修复)，通过循环求解这两个子问题，即可获得去噪后的图像X。具体为：

步骤(1)：稀疏学习秩估计，固定{s_k}，{t_k}的估计可以表示为：

由于一个准确的矩阵X_k的秩估计取决于在合适的稀疏先验下对s_k及t_k稀疏性的估计。另外，在固定s_k的情况下，可通过变量s_k学习到一个针对特定数据的t_k的稀疏先验。公式(16)中t_k的求解问题相当于在已知稀疏先验的条件下，求解最大后验概率问题。

将奇异值的噪声观测项表示为s_k＝t_k+ε_k，ε_k是噪声腐蚀项，ε_k服从均值为0的高斯分布，则s_k的似然函数表示为：

其中，diag(·)表示matlab中的diag函数。本发明采用重加权的拉普拉斯先验来模拟t_k中的结构信息。

首先，令t_k服从高斯分布，t_k的概率密度函数表示为：

其中，γ_k为高斯分布中的方差向量，Γ_k＝diag(γ_k)；

其次，令γ_k服从伽马分布，γ_k的概率密度函数表示为：

其中，κ_k为伽马分布中的形状参数向量。

为了能够数据自适应地获取参数ε_k、γ_k和κ_k，使用基于隐变量的贝叶斯稀疏学习方法对上述参数进行学习，得到一个针对特定数据的稀疏先验：

其中，辅助变量t_k在积分过程中被消除， 表示加权迹范数，利用关系将上式转化为一个统一的变分框架：

接下来利用交替最小化算法对上式进行优化求解，具体步骤如下：

步骤a：初始化，记初始迭代次数为0，分别令ε_k、γ_k、κ_k为值为1的R_k维向量，即ε_k＝γ_k＝κ_k＝[1,...,1]^T，s_k初始化为奇异值分解Yk得到的奇异值组成的向量；

步骤b：更新t_k：

t_k＝Γ_k(∑_bk)^-1s_k (22)

其中，Γ_k＝diag(γ_k)；

步骤c：更新ε_k；

其中，Q＝t_k-s_k，

步骤d：更新γ_k：

其中，

步骤e：更新κ_k：

步骤f：迭代次数加1，如果迭代次数大于15或相邻两次迭代得到的t_k的差值小于10^-6，则停止，此时得到最终的t_k；否则，回到步骤b。

按以上过程求解，即可得到{t_k}，k＝1,...,K。

步骤(2)：低秩矩阵修复，固定{t_k}，{s_k}的估计可以表示为：

因为对{s_k}的求解可以转化为一个等价问题：

其中，T_k为对角矩阵，t_k为T_k对角线元素对应的向量。对于每个X_k上式的二次最优化问题的闭式解为：

其中，按照公式(28)即可求得去噪后图像X＝[X₁,…,X_k,…,X_K]，k＝1,...,K。

本发明提出的一种基于类内结构表示的高光谱图像去噪方法，不需要噪声方差或其它先验信息，且可以适用于各种不同噪声情况，具有较好的去噪效果和良好的适应性。在Washington DC Mall(https://engineering.purdue.edu/～biehl/MultiSpec/hyperspectral.html)数据集上，均一高斯噪声(各个波段施加相同强度的高斯噪声)情况下，利用本发明方法去噪后图像的PSNR(峰值信噪比)值比利用当前最优方法CMESSC(Structured Sparse Coding based Hyperspectral Imagery Denoising with Intra-cluster Filtering，IEEE Transactions on Geoscience&Remote Sensing,2017,PP(99):1-17.)去噪后图像的PSNR值高2.46dB；非均一高斯噪声(各个波段施加不同强度的高斯噪声)情况下，利用本发明方法去噪后图像的PSNR值比利用CMESSC方法去噪后图像的PSNR值高3.58dB。

Claims (1)

1.一种基于类内低秩结构表示的高光谱图像去噪方法，其特征在于包括以下步骤：

步骤一：首先，将给定的三维高光谱图像数据转化成二维矩阵R，矩阵大小为n_b×n_P，其中，n_b为光谱波段个数，n_p为像元的个数；然后，采用K-means方法对矩阵R进行聚类处理，得到K个类Y₁,…,Y_K，K的取值范围为30～90，每一类Y_k所包含的像元数为n_k，k＝1,…,K；

步骤二：建立去噪模型如下：

其中，λ为噪声方差，||·||_F表示Frobenius范数，X＝[X₁,…,X_K]为去噪后图像，p(X_k|Θ_lk)为X_k的低秩先验，Θ_lk表示低秩先验p(X_k|Θ_lk)中的参数；利用奇异值分解的方法将低秩矩阵X_k表示为将公式(1)转化为：

其中，s_k为由S_k对角线上的奇异值组成的向量，即奇异值的噪声观测项；p(s_k|Θ_sk)表示奇异值向量s_k的稀疏先验，Θ_sk表示稀疏先验p(s_k|Θ_sk)中的参数；U_k和V_k通过奇异值分解Y_k得到；

步骤三：首先，引入辅助变量t_k，将公式(2)转化为：

其中，β为预定义的标量参数，β＝0.1；p(t_k|Θ_tk)表示辅助变量t_k的稀疏先验，Θ_tk表示稀疏先验p(t_k|Θ_tk)中的参数，||·||₂表示L₂范数；

然后，采用半二次分裂算法对公式(3)进行求解，得到t_k，k＝1,…,K，并得到去噪后的图像X，具体为：

步骤1：稀疏学习秩估计，固定{s_k}，{t_k}的估计表示为：

令s_k＝t_k+ε_k，ε_k是噪声腐蚀项，ε_k服从均值为0的高斯分布，s_k的似然函数表示为∝表示正比于，diag(·)表示matlab中的diag函数；使用重加权拉普拉斯先验对t_k中的结构信息进行表示，令t_k服从高斯分布，t_k的概率密度函数表示为γ_k为高斯分布中的方差向量，Γ_k＝diag(γ_k)；令γ_k服从伽马分布，γ_k的概率密度函数表示为κ_k为伽马分布中的形状参数向量，R_k＝min(n_b,n_k)；使用基于隐变量的贝叶斯稀疏学习方法对参数ε_k、γ_k和κ_k进行学习，得到如下稀疏先验：

其中，表示加权迹范数， 公式(5)转化为一个统一的变分框架，即：

利用交替最小化算法对公式(6)进行优化求解，具体步骤如下：

步骤a：初始化，记初始迭代次数为0，分别令ε_k、γ_k、κ_k为值为1的R_k维向量，即ε_k＝γ_k＝κ_k＝[1,...,1]^T，s_k初始化为对Yk进行奇异值分解得到的奇异值组成的向量；

步骤b：按照t_k＝Γ_k(∑_bk)^-1s_k更新t_k，其中Γ_k＝diag(γ_k)；

步骤c：按照更新ε_k，其中，Q＝t_k-s_k，

步骤d：按照更新γ_k，其中，

步骤e：根据更新参数κ_k；

步骤f：迭代次数加1，如果迭代次数大于15或t_k的更新差异小于10^-6，则停止，此时得到最终的t_k；否则，回到步骤b；

步骤(2)：低秩矩阵修复，固定{t_k}，{s_k}的估计表示为：

根据将对{s_k}的求解转化为以下等价问题：

其中，T_k为对角矩阵，t_k为T_k对角线元素组成的向量；公式(8)的闭式解为：

其中，按公式(9)求X_k，k＝1,...,K，即获得去噪后图像X＝[X₁,…,X_K]。