CN106067165B - 基于聚类化稀疏随机场的高光谱图像去噪方法 - Google Patents

Abstract

本发明公开了一种基于聚类化稀疏随机场的高光谱图像去噪方法，用于解决现有高光谱图像去噪方法去噪性能差的技术问题。技术方案是根据光谱特性将高光谱图像中的像素聚成若干个类别，利用马尔科夫随机场模型构造图像的先验模型。先验中定义的类内结构稀疏性势能函数和图结构势能函数，分别约束了类内像素光谱维上的相关性和空间维上的相似性。并且，采用正则化回归模型，联合表征了先验学习模型和去噪模型。由于联合考虑了高光谱图像光谱维上的相关性和空间维上的相似性，去噪性能得以提升。在CAVE数据集上的去噪实验表明，当噪声图像的信噪比为16.5分贝时，去噪结果获得的峰值信噪比达到了33.3分贝。

Description

基于聚类化稀疏随机场的高光谱图像去噪方法

技术领域

本发明涉及一种高光谱图像去噪方法，特别是涉及一种基于聚类化稀疏随机场的高光谱图像去噪方法。

背景技术

高光谱图像包含的光谱信息，反映了成像场景在不同电磁波段下的反射比，能够对场景进行更为全面和可靠的描述。因此，高光谱图像在目标探测、医疗诊断和农业监控等诸多方面发挥出了极大的优势。然而，在实际成像过程中，高光谱图像往往会不可避免地受到噪声的污染，导致图像质量下降，影响了高光谱图像在诸多应用上的性能发挥。因此，高光谱图像去噪成为了其应用中的关键步骤。

文献“Maggioni M,Katkovnik V,Egiazarian K,et al.Nonlocal transform-domain filter for volumetric data denoising and reconstruction[J].ImageProcessing,IEEE Transactions on,2013,Vol.22(1),p119-133.”公开了一种有效的高光谱图像去噪算法。该算法将相似的三维图像块聚集成数据组，通过同时挖掘组内的局部相似性和不同组之间的非局部相似性，并且利用变换域的滤波技术，实现高光谱图像的去噪。然而，该算法忽略了高光谱图像中光谱维上存在的强相关性，去噪性能有限。

发明内容

为了克服现有高光谱图像去噪方法去噪性能差的不足，本发明提供一种基于聚类化稀疏随机场的高光谱图像去噪方法。该方法根据光谱特性将高光谱图像中的像素聚成若干个类别，利用马尔科夫随机场模型构造图像的先验模型。先验中定义的类内结构稀疏性势能函数和图结构势能函数，分别约束了类内像素光谱维上的相关性和空间维上的相似性。并且，采用正则化回归模型，联合表征了先验学习模型和去噪模型。由于联合考虑了高光谱图像光谱维上的相关性和空间为维上的相似性，去噪性能得以提升。

本发明解决其技术问题所采用的技术方案是：一种基于聚类化稀疏随机场的高光谱图像去噪方法，其特点是包括以下步骤：

步骤一、建立噪声图像的观测模型F＝X+N，其中表示噪声图像，表示噪声；假设噪声N为高斯白噪声，服从矩阵正太分布其中，Σ_n＝diag(λ)表示以λ为对角线元素的对角矩阵，表示行间协方差矩阵；为列向量，控制不同波段的噪声程度；I为对应尺寸的单位矩阵，表示列间协方差矩阵；假设干净的高光谱图像X在光谱字典上稀疏表示为X＝ΦY，表示稀疏表示系数矩阵，n_d表示字典中原子的个数；结合噪声分布、观测模型以及图像的稀疏表示模型，获得似然函数，

其中，表示矩阵Q的加权迹范数，Q＝ΦY-F；

步骤二、根据光谱特性将X中的像素聚为K个类别，表示第k个类别中所有像素的集合，n_k表示该类别中像素的个数，k＝1,...,K；X_k在字典Φ下的稀疏表示系数矩阵表示为针对X的稀疏表示系数矩阵Y设计聚类化稀疏随机场先验；首先，给出Y的马尔科夫随机场模型，如下

其中，Γ_k＝diag(γ_k)，Σ_k＝diag(η_k)；表示Y_k自我表示的权重矩阵，其对角线元素全为0；为结构稀疏性势能函数，为图结构势能函数；其次，在结构稀疏性势能函数中，为了描述Y_k内在的结构稀疏性，进一步对γ_k引入如下的伽马分布，

其中，接着，在图结构势能函数中，为了描述类内光谱自我表示误差的稀疏性，对η_k引入如下的伽马分布，

其中，最后，为了避免先验学习过程中的过拟合现象，对W_k中的每一列引入独立同分布的正太分布先验，如下

其中，表示W_k中的第i个列向量，∈为给定参数，∈＝20；表示均值为0，协方差矩阵为∈^-1I的正太分布；为了便于表示，将上述具有级联结构的先验模型，统一表示为如下的先验形式，

其中，变量表示先验模型中的所有参数，p_csf(Y|Θ)则表示聚类化稀疏随机场先验；

步骤三、将先验学习模型和去噪模型进行联合建模；首先，利用最大化后验概率估计从噪声图像中估计先验模型参数Θ和噪声程度参数λ，如下

然而，式子(7)无法求解；为此，对先验模型中的进行如下近似，

其中，M_k＝Y′_kW_k，Y′_k表示在前一次迭代中与X_k对应的稀疏表示系数矩阵；将式子(8)代入到式子(7)中，通过积分，并引入-2log运算，式子(7)简化为：

其中，tr(·)表示迹范数，F_k为与X_k对应的局部噪声图像；||·||_F表示弗罗贝尼乌斯范数；进一步，将已知的数学关系，代入式子(9)中，得到联合先验学习模型和去噪模型的正则化回归模型，如下

步骤四、采用交替最小化方法将式子(10)分解为若干个子问题进行迭代求解，直至收敛；在每个子问题中仅优化一个变量而固定其他所有变量；具体步骤如下：

①初始化λ，为对应长度的全1向量，利用现有去噪算法的结果初始化X，计数变量初始化t＝0；

②将X中的每一列看作一个样本，进行主成分分解，用得到的正交基构成字典Φ，并利用K-均值聚类方法依据光谱特性将X中的像素聚为K＝30个类别；

③固定X，得到关于权重矩阵W_k的优化子问题，如下

求解得到更新所有的W_k，k＝1,...,K；

④固定γ_k,η_k和λ，根据上一次迭代中得到的X_k，计算M_k＝Φ^TX_kW_k，然后得到关于Y的优化子问题，如下

求解得到更新所有的Y_k，并根据聚类关系重新构成Y，然后更新高光谱图像X＝ΦY；

⑤固定Y_k,η_k和λ，得到关于γ_k的子问题，如下

求解得到γ_k中第j个元素γ_jk的更新形式，如下

其中，z_j为向量中的第j个元素，为向量中的第j个元素；此处diag(·)表示取矩阵对角线元素组成列向量；表示中的第j个元素；根据关系式(14)，更新γ_k中的所有元素；

⑥固定γ_k，得到的子问题

求解得到中第j个元素更新中的所有元素；

⑦固定Y_k,M_k,ν_k和λ，得到关于η_k的优化子问题，如下

求解得到η_k中第j个元素η_jk的更新形式，如下

其中，α_j为向量中的第j个元素，为向中的第j个元素；diag(·)表示取矩阵对角线元素组成列向量；ν_jk表示ν_k中的第j个元素；根据关系式(17)，更新η_k中的所有元素；

⑧固定η_k，得到关于ν_k的优化子问题，如下

求解得到ν_k中第j个元素ν_jk＝2/η_jk，更新ν_k中的所有元素；

⑨固定Y_k,γ_k和η_k，得到关于λ的优化子问题，如下

求解得到λ中第j个元素λ_jk的更新形式，如下

其中，β_jk是向量中的第j个元素，是向量中的第j个元素；diag(·)表示取矩阵对角线元素组成列向量；利用式子(20)更新λ中的所有元素；

⑩计数变量t加1，如果t<4，循环执行步骤②～⑩；否则，退出程序。

本发明的有益效果是：该方法根据光谱特性将高光谱图像中的像素聚成若干个类别，利用马尔科夫随机场模型构造图像的先验模型。先验中定义的类内结构稀疏性势能函数和图结构势能函数，分别约束了类内像素光谱维上的相关性和空间维上的相似性。并且，采用正则化回归模型，联合表征了先验学习模型和去噪模型。由于联合考虑了高光谱图像光谱维上的相关性和空间为维上的相似性，去噪性能得以提升。在CAVE数据集上的去噪实验表明，当噪声图像的信噪比为16.5分贝时，去噪结果获得的峰值信噪比达到了33.3分贝。

下面结合具体实施方式对本发明作详细说明。

具体实施方式

本发明基于聚类化稀疏随机场的高光谱图像去噪方法具体步骤如下：

对于包含n_b个波段，每个波段空间上包含n_r行和n_c列的高光谱图像，将其每一个波段拉伸成为一个行向量，重新组成一个二维矩阵，(n_p＝n_r×n_c表示一个波段上的像素数目)。X中的任一行对应特定波段下所有像素值；每一列则对应某个像素的光谱。本发明设计的图像去噪方法，具体包含以下四个步骤：

1、建立噪声图像的成像模型。

本发明主要处理与信号独立的加性噪声，对应的观测模型为F＝X+N，其中表示噪声图像，表示噪声。假设噪声N为高斯白噪声，服从矩阵正太分布其中，Σ_n＝diag(λ)表示以λ为对角线元素的对角矩阵，表示行间协方差矩阵。为列向量，控制不同波段的噪声程度。I为对应尺寸的单位矩阵，表示列间协方差矩阵。假设干净的高光谱图像X可以在光谱字典上稀疏表示为X＝ΦY，表示稀疏表示系数矩阵，n_d表示字典中原子的个数。结合噪声分布，观测模型，以及图像的稀疏表示模型，可以获得如下的似然函数，

其中，表示矩阵Q的加权迹范数，此处Q＝ΦY-F。

2、建立聚类化稀疏随机场的高光谱图像先验模型。

根据光谱特性将X中的像素聚为K个类别，表示第k个类别中所有像素的集合，n_k表示该类别中像素的个数，k＝1,...,K。X_k在字典Φ下的稀疏表示系数矩阵则可以表示为本发明针对X的稀疏表示系数矩阵Y设计了聚类化稀疏随机场先验。首先，给出了Y的马尔科夫随机场模型，如下

其中，Γ_k＝diag(γ_k)，Σ_k＝diag(η_k)。表示Y_k自我表示的权重矩阵，其对角线元素全为0。为结构稀疏性势能函数，为图结构势能函数。其次，在结构稀疏性势能函数中，为了描述Y_k内在的结构稀疏性，进一步对γ_k引入如下的伽马分布，

其中，接着，在图结构势能函数中，为了描述类内光谱自我表示误差的稀疏性，对η_k引入如下的伽马分布，

其中，最后，为了避免先验学习过程中的过拟合现象，对W_k中的每一列引入独立同分布的正太分布先验，如下

其中，表示W_k中的第i个列向量，∈为给定参数，本发明中∈＝20。表示均值为0，协方差矩阵为∈^-1I的正太分布。为了便于表示，将上述具有级联结构的先验模型，统一表示为如下的先验形式，

其中，变量表示先验模型中的所有参数，p_csf(Y|Θ)则表示本发明提出的聚类化稀疏随机场先验。

3、建立去噪模型。

本发明提出了一种统一的正则化回归模型，将先验学习模型和去噪模型进行联合建模。首先，利用最大化后验概率估计从噪声图像中估计先验模型参数Θ和噪声程度参数λ，如下

然而，式子(7)无法求解。为此，对先验模型中的进行如下近似，

其中，M_k＝Y′_kW_k，Y′_k表示在前一次迭代中与X_k对应的稀疏表示系数矩阵。将式子(8)代入到式子(7)中，通过积分，并引入-2log运算，式子(7)可简化为：

其中，tr(·)表示迹范数，F_k为与X_k对应的局部噪声图像。||·||_F表示弗罗贝尼乌斯范数。进一步，将已知的数学关系，代入式子(9)中，得到联合先验学习模型和去噪模型的正则化回归模型，如下

4、模型求解。

本发明采交替最小化方法将式子(10)分解为若干个子问题进行迭代求解，直至收敛。在每个子问题中仅优化一个变量而固定其他所有变量。具体步骤如下：

①初始化λ，为对应长度的全1向量，利用现有去噪算法的结果初始化X，计数变量初始化t＝0；

②将X中的每一列看作一个样本，进行主成分分解，用得到的正交基构成字典Φ，并利用K-均值聚类方法依据光谱特性将X中的像素聚为K＝30个类别；

③固定X，得到关于权重矩阵W_k的优化子问题，如下

求解得到利用该规则更新所有的W_k，k＝1,...,K；

④固定γ_k,η_k和λ，根据上一次迭代中得到的X_k，计算M_k＝Φ^TX_kW_k，然后得到关于Y的优化子问题，如下

求解得到利用该关系更新所有的Y_k，并根据聚类关系重新构成Y，然后更新高光谱图像X＝ΦY；

⑤固定Y_k,η_k和λ，得到关于γ_k的子问题，如下

求解得到γ_k中第j个元素γ_jk的更新形式，如下

其中，z_j为向量中的第j个元素，为向量中的第j个元素。此处diag(·)表示取矩阵对角线元素组成列向量。表示中的第j个元素。根据关系(14)，更新γ_k中的所有元素；

⑥固定γ_k，得到的子问题

求解得到中第j个元素利用该关系更新中的所有元素；

⑦固定Y_k,M_k,ν_k和λ，得到关于η_k的优化子问题，如下

求解得到η_k中第j个元素η_jk的更新形式，如下

其中，α_j为向量中的第j个元素，为向中的第j个元素。此处diag(·)表示取矩阵对角线元素组成列向量。ν_jk表示ν_k中的第j个元素。根据关系式(17)，更新η_k中的所有元素；

⑧固定η_k，得到关于ν_k的优化子问题，如下

求解得到ν_k中第j个元素ν_jk＝2/η_jk，利用该关系更新ν_k中的所有元素；

⑨固定Y_k,γ_k和η_k，得到关于λ的优化子问题，如下

求解得到λ中第j个元素λ_jk的更新形式，如下

其中，β_jk是向量中的第j个元素，是向量中的第j个元素。此处diag(·)表示取矩阵对角线元素组成列向量。利用式子(20)更新λ中的所有元素；

⑩计数变量t加1，如果t<4，循环执行步骤②～⑩；否则，退出程序。

总之，本发明采用的图像先验同时描述了高光谱图像光谱维的相关性和空间维上的相似性，较现有的去噪方法，更加全面地表征了高光谱图像的三维特点，并且采用的联合去噪模型能够进行自适应的先验学习和图像去噪，最终实现了去噪性能的有效提升。在CAVE数据集上的去噪实验表明，当噪声图像的信噪比为16.5分贝时，本发明去噪结果获得的峰值信噪比达到了33.3分贝。

Claims (1)

1.一种基于聚类化稀疏随机场的高光谱图像去噪方法，其特征在于包括以下步骤：

步骤一、建立噪声图像的观测模型F＝X+N，其中表示噪声图像，表示噪声；假设噪声N为高斯白噪声，服从矩阵正太分布其中，Σ_n＝diag(λ)表示以λ为对角线元素的对角矩阵，表示行间协方差矩阵；为列向量，控制不同波段的噪声程度；I为对应尺寸的单位矩阵，表示列间协方差矩阵；假设干净的高光谱图像X在光谱字典上稀疏表示为X＝ΦY，表示稀疏表示系数矩阵，n_d表示字典中原子的个数；结合噪声分布、观测模型以及图像的稀疏表示模型，获得似然函数，

其中，表示矩阵Q的加权迹范数，Q＝ΦY-F；

步骤二、根据光谱特性将X中的像素聚为K个类别，表示第k个类别中所有像素的集合，n_k表示该类别中像素的个数，k＝1,...,K；X_k在字典Φ下的稀疏表示系数矩阵表示为针对X的稀疏表示系数矩阵Y设计聚类化稀疏随机场先验；首先，给出Y的马尔科夫随机场模型，如下

其中，表示Y_k自我表示的权重矩阵，其对角线元素全为0；为结构稀疏性势能函数，为图结构势能函数；其次，在结构稀疏性势能函数中，为了描述Y_k内在的结构稀疏性，进一步对γ_k引入如下的伽马分布，

其中，接着，在图结构势能函数中，为了描述类内光谱自我表示误差的稀疏性，对η_k引入如下的伽马分布，

其中，最后，为了避免先验学习过程中的过拟合现象，对W_k中的每一列引入独立同分布的正太分布先验，如下

其中，表示W_k中的第i个列向量，∈为给定参数，∈＝20；表示均值为0，协方差矩阵为∈^-1I的正太分布；为了便于表示，将上述具有级联结构的先验模型，统一表示为如下的先验形式，

其中，变量表示先验模型中的所有参数，p_csf(Y|Θ)则表示聚类化稀疏随机场先验；

步骤三、将先验学习模型和去噪模型进行联合建模；首先，利用最大化后验概率估计从噪声图像中估计先验模型参数Θ和噪声程度参数λ，如下

然而，式子(7)无法求解；为此，对先验模型中的进行如下近似，

其中，M_k＝Y′_kW_k，Y′_k表示在前一次迭代中与X_k对应的稀疏表示系数矩阵；将式子(8)代入到式子(7)中，通过积分，并引入-2log运算，式子(7)简化为：

其中，tr(·)表示迹范数，F_k为与X_k对应的局部噪声图像；||·||_F表示弗罗贝尼乌斯范数；进一步，将已知的数学关系，中，得到联合先验学习模型和去噪模型的正则化回归模型，如下

步骤四、采用交替最小化方法将式子(10)分解为若干个子问题进行迭代求解，直至收敛；在每个子问题中仅优化一个变量而固定其它所有变量；具体步骤如下：

①初始化λ，为对应长度的全1向量，利用现有去噪算法的结果初始化X，计数变量初始化t＝0；

②将X中的每一列看作一个样本，进行主成分分解，用得到的正交基构成字典Φ，并利用K-均值聚类方法依据光谱特性将X中的像素聚为K＝30个类别；

③固定X，得到关于权重矩阵W_k的优化子问题，如下

求解得到更新所有的W_k，k＝1,...,K；

④固定γ_k,η_k和λ，根据上一次迭代中得到的X_k，计算M_k＝Φ^TX_kW_k，然后得到关于Y的优化子问题，如下

求解得到更新所有的Y_k，并根据聚类关系重新构成Y，然后更新高光谱图像X＝ΦY；

⑤固定Y_k,η_k和λ，得到关于γ_k的子问题，如下

求解得到γ_k中第j个元素γ_jk的更新形式，如下

其中，z_j为向量中的第j个元素，为向量中的第j个元素；此处diag(·)表示取矩阵对角线元素组成列向量；表示中的第j个元素；根据关系式(14)，更新γ_k中的所有元素；

⑥固定γ_k，得到的子问题

求解得到中第j个元素更新中的所有元素；

⑦固定Y_k,M_k,ν_k和λ，得到关于η_k的优化子问题，如下

求解得到η_k中第j个元素η_jk的更新形式，如下

其中，α_j为向量中的第j个元素，为向量中的第j个元素；diag(·)表示取矩阵对角线元素组成列向量；ν_jk表示ν_k中的第j个元素；根据关系式(17)，更新η_k中的所有元素；

⑧固定η_k，得到关于ν_k的优化子问题，如下

求解得到ν_k中第j个元素ν_jk＝2/η_jk，更新ν_k中的所有元素；

⑨固定Y_k,γ_k和η_k，得到关于λ的优化子问题，如下

求解得到λ中第j个元素λ_jk的更新形式，如下

其中，β_jk是向量中的第j个元素，是向量中的第j个元素；diag(·)表示取矩阵对角线元素组成列向量；利用式子(20)更新λ中的所有元素；

⑩计数变量t加1，如果t＜4，循环执行步骤②～⑩；否则，退出程序。