基于非对称高斯混合模型的遥感图像去噪算法

Remote Sensing Image Denoising Algorithm Based on Asymmetric Gaussian Mixture Model

摘要: 遥感图像在采集时，由于受到各种因素干扰，会严重影响图像的视觉效果，进而影响后续处理的准确性。因此，对遥感图像的噪声进行准确地建模是解决遥感图像噪声问题的关键。编码噪声分布最常选择是高斯分布、拉普拉斯分布和高斯混合分布，但它们总是与现实世界的遥感图像噪声不相容。考虑到遥感图像同时存在对称和非对称的噪声分布，本文在高斯混合分布基础上，引入了非对称参数，构建了一个基于非对称高斯混合分布模型(AMoG)的遥感图像去噪算法。该算法使用低秩矩阵分解将遥感图像近似为两个因子矩阵的乘积。对于模型的参数，使用了EM算法进行迭代更新。在合成数据集和真实数据集上的大量实验结果表明，该模型在PSNR、SSIM、FSIM、ERGA、SAM五种评价指标上均表现良好，表明了该算法在遥感图像去噪方面具有一定的优越性。

Abstract: Remote sensing images are often subject to various interferences during acquisition, which seriously affects the visual effect of the images, and then affects the accuracy of subsequent processing. Therefore, accurate modeling of the noise of remote sensing images is the key to solving the noise problem of remote sensing images. The most common choices for coded noise distributions are Gaussian, Laplace, and Gaussian mixtures, but they are always incompatible with real-world remote sensing image noise. Considering that there are both symmetrical and asymmetric noise distribution in remote sensing images, this paper introduces asymmetric parameters on the basis of Gaussian mixed distribution, and constructs a remote sensing image denoising algorithm based on asymmetric Gaussian mixed distribution model (AMoG). The algorithm uses low-rank matrix factorization to approximate the remote sensing image as the product of two-factor matrices. For the parameters of the model, the EM algorithm was used for iterative update. A large number of experimental results on synthetic datasets and real datasets show that the model performs well in five evaluation indexes: PSNR, SSIM, FSIM, ERGA and SAM, indicating that the algorithm has certain advantages in remote sensing image denoising.

1. 引言

近年来，基于机载和星载的遥感探测技术不断发展，为终端用户提供了丰富的光谱、空间、时间等信息。遥感图像已由几十个波段扩展至数百个波段，能够更全面地反映目标在目标区域的光谱特征。但是，由于遥感图像数据维度高、维数大、像元混合、采集过程中的光散射效应等原因，导致其出现高度的非线性与复杂性，给已有的数据处理与分析方法带来了巨大的挑战。为此，对遥感图像进行噪声去除已引起国内外学者的广泛重视。

图像去噪是指根据已获得的图像的先验信息，通过一个理想的算法或模型，将图像中的噪声降到最低，消除或消除，以尽可能多地保留原始图像的信息，提升图像的品质，为后续的处理提供帮助。经典的图像去噪方法主要包括基于滤波理论的去噪方法[1]、基于正则化理论的去噪方法[2] [3]、基于图像块相似性和稀疏性的去噪方法[4] [5]、基于稀疏表示和低秩理论的去噪方法[6] [7]、基于统计先验知识的去噪方法和噪声建模等。

噪声建模已广泛应用于遥感图像中，但是大多数遥感图像去噪方法主要考虑的是遥感图像先验谱和空间知识，只使用L₂损失项来校正重建与原始遥感图像之间的偏差。从统计理论的角度来看，利用这种损失项意味着遥感图像噪声遵循的是高斯分布。然而，在实际情况中，遥感图像噪声通常具有更为复杂的统计结构。这意味着这种简单的损失项过于主观，无法反映真实场景中非独立同分布的噪声情况。于是人们考虑能够超越高斯分布的更为复杂的去噪框架，近年来，低秩矩阵模型渐渐引起了关注。

利用低秩恢复或低秩矩阵逼近的方式进行图像处理是图像去噪的研究内容之一。在基于L₂范数的框架下Srebro等人基于期望最大化(EM) [8]算法构造了加权低秩逼近的计算策略[9]，计算结果表现良好。但该方法对初值较为敏感，同时每步迭代都需要进行奇异值分解，从而增大了计算代价。Buchanan等人提出了一种含有阻尼因子的牛顿法(DN) [10]，该方法基于目标函数的二阶泰勒展开，利用了二阶信息，有收敛速度快，计算精度高的优点。但是对于规模较大的数据，该方法每步迭代都需要计算一个庞大的Hessian矩阵，导致计算效率降低。为了解决这一问题，Okatani等人借助了Wiberg算法，提出了新的求解方法[11]，在L₂范数低秩矩阵分解问题中表现出了较好的计算效率。基于以上研究，Okatani等人提出Damped Wiberg (DW)算法[12]，进一步提升了算法的计算效率，同时改善了算法易收敛到局部最优的现象。

针对遥感图像噪声的去除，Meng等人提出将噪声建模为高斯混合(MoG) [13]，MoG建模能够处理具有未知统计分布的复杂遥感图像噪声。然而，该模型假设噪声相同和独立分布，这在很大程度上低估了非独立同分布的真实遥感图像背后的属性。Cao等人将噪声建模为一个更一般的非独立同分布指数功率(PMoEP) [14]混合分布，在真实噪声场景下实现了对遥感图像进行去噪。Chen等人首先观察到不同波段图像的噪声分布在统计上是不同的，但有一定的相关性，提出了非独立同分布高斯混合(NMoG)模型[15]来描述这一现象。这些工作取得了显著的进展，加深了人们对遥感图像去噪的认识。然而，MoG及其变体的成功本质上在于噪声普遍的近似特性，这使得MoG能够提供一种方法来模拟未知噪声。但以上这些方法却对精确地模拟复杂的遥感图像噪声是非常困难的。

这一困境促使人们继续探索。Xu等人揭示了合成噪声和真实世界的遥感图像噪声都表现出不对称的特性，提出了一种基于非对称拉普拉斯噪声的矩阵分解(BALMF)方法[16]，以提高去除混合噪声的性能。Wang等人将深度图像先验与非对称拉普拉斯分布相结合，提出了非对称拉普拉斯噪声建模深度图像先验(ALDIP)的方法[17]。然而，BALMF或ALDIP方法对于类拉普拉斯噪声是最优的，如脉冲、截止日期和条纹噪声，即稀疏噪声。

针对上述问题，本文假设每个波段的噪声服从非对称高斯混合分布，构建一种基于非对称高斯混合噪声的矩阵分解(AMoG)模型，以便更灵活地拟合混合噪声。该模型的参数使用EM算法进行迭代更新。通过对模拟数据和真实数据的一系列实验结果，表明了本文方法在去除遥感图像的混合噪声方面的优越性。

2. 相关工作

2.1. 非对称高斯混合分布

2.1.1. 高斯混合分布

高斯混合分布(GMM)是单高斯分布的延伸，就是用多个高斯概率密度函数精确地量化变量分布，是将变量分布分解为若干基于高斯概率密度函数分布的统计模型。它可以平滑地近似任何形状的分布。反过来解释就是，R个单高斯分布混合在一起，生成的模型就是高斯混合模型。

单一高斯分布由两个参数完全描述：均值$\mu$ ，方差$\sigma$ 。均值确定了分布的中心位置，标准差则决定了分布的离散程度。高斯分布的概率密度函数为：

$f\left(x;\mu ,\sigma \right)=\frac{1}{\sqrt{2\pi }\sigma }\exp \left[-\frac{{\left(x-\mu \right)}^2}{2{\sigma }^2}\right]$ (1)

高斯混合分布的概率密度函数可以表示为：

$f\left(x;\mu ,\sigma ,r\right)={\displaystyle \sum _{r=1}^R{\alpha }_{r}f\left(x;\mu ,\sigma \right)}$ (2)

其中，R为分布的个数；$0\le {\alpha }_r\le 1$ 为第r个高斯的权重，且${\displaystyle \sum _{r=1}^R{\alpha }_r}=1$ 。

2.1.2. 非对称高斯混合分布

基于高斯分布引申出的非对称高斯分布(AG)如下：

$f_{AG}\left(x;\mu ,\sigma ,\kappa \right)=\frac{2}{\sqrt{2\pi }\sigma \left(\kappa +1\right)}\times \{\begin{array}{l}\exp \left[-\frac{{\left(x-\mu \right)}^2}{2{\sigma }^2}\right],x\ge \mu \\ \exp \left[-\frac{{\left(x-\mu \right)}^2}{2{\kappa }^2{\sigma }^2}\right],x<\mu \end{array}$ (3)

其中，$\mu$ ，$\sigma$ 和$\kappa$ 分别是均值、方差和不对称参数。

非对称高斯混合分布(AMoG)是由多个单非对称高斯分布的模型通过线性组合而成，其联合概率密度函数如下：

$f_{AMoG}\left(x;\mu ,\sigma ,\kappa ,r\right)={\displaystyle \sum _{r=1}^R{\alpha }_{jr}{f}_{AG}\left(x;\mu ,\sigma ,\kappa \right)}$ (4)

因此，与高斯混合分布相比，非对称高斯混合分布可以更好地建模那些具有复杂的、非线性的、方向不同的数据分布。在实际的应用中，当数据的真实分布在某些方面偏离了标准的高斯形态，或者包含了明显的非线性结构时，非对称高斯混合模型能够提供更准确的拟合。这使得模型能够更精确地描述数据中的模态分布，从而更准确地进行分类、聚类或生成新的样本。

2.2. 低秩矩阵分解

基本的低秩矩阵分解模型[18]如下：

$\underset{U,V}{\min }{\Vert X-U{V}^{\text{T}}\Vert }_{L_P}+\tau \left(U,V\right)$ (5)

其中像素矩阵$X\in {ℝ}^{m\times n}$ 是待逼近矩阵，$U\in {ℝ}^{m\times r}$ ，$V\in {ℝ}^{n\times r}$ 是低维矩阵且$r<\min \left(m,n\right)$ ，$\tau \left(U,V\right)$ 为正则项，避免过度拟合。

该方法与凸松弛法不同，它不是直接优化矩阵变量$X\in {ℝ}^{m\times n}$ ，而是根据目标矩阵是低秩的信息将变量X用其双因子形式$U{V}^{\text{T}}$ 代替，然后直接求解所得到的非凸因子模型。

在很多实际应用中，人们所采集到的大规模的数据往往含有噪声、损毁或者误差等。但这些大规模数据有时候本身具有低秩特性，因此矩阵的低秩分解在对数据的修复或者去除噪声的影响这些方面得到了广泛的研究和应用。Xu等人[19]提出自适应分位数矩阵分解模型以恢复子空间。

超分辨率重建也是低秩矩阵分解的另一个应用场景。在这个场景中，多个低分辨率图像可以通过低秩矩阵分解及其扩展算法进行重构，从而得到一个更高分辨率的图像。Dian等人[20]结合了低秩和稀疏表示，提出了一种高光谱图像超分辨率方法，将融合问题定义为从低分辨率高光谱图像得到的光谱字典估计和两幅图像各自的回归系数。

此外，低秩矩阵分解还可以用于图像压缩，即将大量数据压缩成小量数据以减少存储空间。Haeffele等人[21]提出一种适用于大型数据集的矩阵分解技术，通过使用特定形式的正则化，实现对数据的高效压缩。

3. 基于非对称高斯混合模型的遥感图像去噪算法

3.1. 模型规划

本文将非对称高斯混合分布模型(AMoG)应用到遥感图像去噪中。给定一个波段为B，像素大小为$M\times N$ 的高光谱图像，用$Y\in {ℝ}^{MN\times B}$ 表示，它可以用两个因子矩阵($U\in {ℝ}^{MN\times R}$ 和$V\in {ℝ}^{B\times R}$ )加上噪声$N\in {ℝ}^{MN\times B}$ 来表示：

$Y=X+N=U{V}^{\text{T}}+N.$ (6)

等价于：

$y_{ij}=x_{ij}+n_{ij}=u_{i,:}{v}_{j,:}^{\text{T}}+n_{ij}.$ (7)

其中，$u_{i,:}$ 表示U的第i行，$v_{j,:}$ 表示V的第j行，$n_{ij}$ 是$y_{ij}$ 中的噪声。

由于高光谱图像的噪声是重尾和不对称的，不同的波段的统计特性是不同的，但具有一定的相关性。因此，所提出的模型可以利用非对称高斯混合模型进行建模。根据(3)式，将单一的非对称高斯模型的概率密度函数记作：$X~AG\left(X;\mu ,\sigma ,\kappa \right)$ 。于是，假设每个波段都被非对称高斯混合噪声破坏：

$n_{ij}~{\displaystyle \sum _{r=1}^R{\alpha }_{jr}A{G}_{jr}\left(n_{ij}|0,{\sigma }_{jr},{\kappa }_{jr}\right)}$ $\forall j=1,2,\cdots ,B.$ (8)

当不对称性参数$\kappa =1$ 时，非对称高斯分布退化为高斯分布的形式。这就是说高斯分布是非对称高斯分布的一个特例，非对称高斯分布相对于高斯分布具有更一般性。由于非对称高斯混合模型在非对称高斯模型的基础上，又可以拟合多个分布，这与自然图像数据更为一致，因此对于自然复杂图像的去噪效果会更好。

根据(7)式和(8)式，下列等式成立：

$y_{ij}~{\displaystyle \sum _{r=1}^R{\alpha }_{jr}A{G}_{jr}\left(y_{ij}|x_{ij},{\sigma }_{jr},{\kappa }_{jr}\right)}$ $\forall j=1,2,\cdots ,B.$ (9)

为了便于未知参数的估计，引入潜在变量${\gamma }_{ijr}$ ，其中${\gamma }_{ijr}\in \left\{0,1\right\}$ 且${\displaystyle \sum _{r=1}^R{\gamma }_{ijr}}=1$ ，这里${\gamma }_{ijr}=1$ 表示来自第r个非对称高斯分布，基于以上假设，$x_{ij}$ 的概率密度函数重新表示为：

$p\left(y_{ij}|x_{ij},{\alpha }_{jr},{\sigma }_{jr},{\kappa }_{jr}\right)={\displaystyle \prod _{r=1}^R{\left[{\alpha }_{jr}A{G}_{jr}\left(y_{ij}|x_{ij},{\sigma }_{jr},{\kappa }_{jr}\right)\right]}^{{\gamma }_{ijr}}}$ (10)

令$\Theta =\alpha ,\sigma ,\kappa$ ，其中，$\alpha ={\left\{{\alpha }_{jr}\right\}}_{B\times R}$ ，$\sigma ={\left\{{\sigma }_{jr}\right\}}_{B\times R}$ ，$\kappa ={\left\{{\kappa }_{jr}\right\}}_{B\times R}$ ，得到完整数据的对数似然函数：

$\begin{array}{c}L\left(Y;X,\Theta \right)={\displaystyle \sum _{j=1}^B{\displaystyle \sum _{i=1}^{MN}{\displaystyle \sum _{r=1}^R{\gamma }_{ijr}}}}\left[\log {\alpha }_{jr}+\log A{G}_{jr}\left(y_{ij}|x_{ij},{\sigma }_{jr},{\kappa }_{jr}\right)\right]\\ ={\displaystyle \sum _{j=1}^B{\displaystyle \sum _{i=1}^{MN}{\displaystyle \sum _{r=1}^R{\gamma }_{ijr}[\log {\alpha }_{jr}+\log \frac{2}{\sqrt{2\pi }}+\log \frac{1}{{\sigma }_{jr}}+\log \frac{1}{{\kappa }_{jr}+1}}}}-\frac{1}{2{\eta }_{ijr}{\sigma }_{jr}^2}{\left(y_{ij}-x_{ij}\right)}^2]\end{array}$ (11)

其中，当$y_{ij}\ge x_{ij}$ 时，${\eta }_{ijr}=1$ ；当$y_{ij}<x_{ij}$ 时，${\eta }_{ijr}={\kappa }_{jr}^2$ 。

3.2. 求解算法

在EM算法的框架下，对(11)式中的参数进行更新。

在对参数进行初始化时，采用K均值聚类算法对${\alpha }_{jr}$ 、$x_{ij}$ 、${\sigma }_{jr}$ 进行初始化。在对非对称参数${\kappa }_{jr}$ 进行初始化时，根据图像块的像素值，将其初始化为：

${\kappa }_{jr}=\sqrt{\frac{{\displaystyle \sum _{y_{ij}<x_{ij}}{\Vert y_{ij}-x_{ij}\Vert }_2}}{{\displaystyle \sum _{y_{ij}\ge x_{ij}}{\Vert y_{ij}-x_{ij}\Vert }_2}}}$ (12)

3.2.1. E步

在E步中，分别令$A{G}_{jr}\left(y_{ij}|x_{ij},{\sigma }_{jr},{\kappa }_{jr}\right)$ 在$y_{ij}\ge x_{ij}$ 和$y_{ij}<x_{ij}$ 的部分为$A{G}_{j{r}_1}\left(y_{ij}|x_{ij},{\sigma }_{jr},{\kappa }_{jr}\right)$ 和 $A{G}_{j{r}_2}\left(y_{ij}|x_{ij},{\sigma }_{jr},{\kappa }_{jr}\right)$ ，可以得到每个$y_{ij}$ 属于各个混合成分的后验概率：

$\widehat{{\gamma }_{ijr}^m}=E\left({\gamma }_{ijr}^m\right)=\frac{{\alpha }_{jr}A{G}_{j{r}_m}\left(y_{ij}|x_{ij},{\sigma }_{jr},{\kappa }_{jr}\right)}{{\displaystyle \sum _{r=1}^R{\alpha }_{jr}A{G}_{jr}\left(y_{ij}|x_{ij},{\sigma }_{jr},{\kappa }_{jr}\right)}},m=1,2$ (13)

3.2.2. M步

通过更新迭代(11)式的参数使其达到最大化，只保留函数中与参数有关的项，目标函数如下：

$l=\max L={\displaystyle \sum _{j=1}^B{\displaystyle \sum _{i=1}^{MN}{\displaystyle \sum _{r=1}^R{\gamma }_{ijr}[\log {\alpha }_{jr}+\log \frac{1}{{\sigma }_{jr}}+\log \frac{1}{{\kappa }_{jr}+1}}}}-\frac{1}{2{\eta }_{ijr}{\sigma }_{jr}^2}{\left(y_{ij}-x_{ij}\right)}^2]$ (14)

更新$\sigma$ 、$\alpha$ ：

${\sigma }_{jr}^2=\frac{{\displaystyle \sum _{y_{ij}\ge x_{ij}}{\gamma }_{ijr}^1\times {\left(y_{ij}-x_{ij}\right)}^2}+{\displaystyle \sum _{y_{ij}<x_{ij}}{\gamma }_{ijr}^2\times \frac{1}{{\kappa }_{jr}^2}{\left(y_{ij}-x_{ij}\right)}^2}}{{\displaystyle \sum _{y_{ij}\ge x_{ij}}{\gamma }_{ijr}^1}+{\displaystyle \sum _{y_{ij}<x_{ij}}{\gamma }_{ijr}^2}}$ (15)

${\alpha }_{jr}=\frac{{\displaystyle \sum _{y_{ij}\ge x_{ij}}{\gamma }_{ijr}^1}+{\displaystyle \sum _{y_{ij}<x_{ij}}{\gamma }_{ijr}^2}}{MN}$ (16)

更新$\kappa$ ：

关于${\kappa }_{jr}$ 的目标函数可以写为：

$l_1={\displaystyle \sum _{j=1}^B{\displaystyle \sum _{i=1}^{MN}{\displaystyle \sum _{r=1}^R{\gamma }_{ijr}\left[\log \frac{1}{{\kappa }_{jr}+1}-\frac{1}{2{\eta }_{ijr}{\sigma }_{jr}^2}{\left(y_{ij}-x_{ij}\right)}^2\right]}}}$ (17)

对上式中的${\kappa }_{jr}$ 求偏导：

$\frac{\partial l_1}{\partial {\kappa }_{jr}}={\gamma }_{ijr}\left[-\frac{1}{{\kappa }_{jr}+1}+\frac{1}{{\kappa }_{jr}^3{\sigma }_{jr}^2}{\left(y_{ij}-x_{ij}\right)}^2\right]$ (18)

由于$\frac{\partial l_1}{\partial {\kappa }_{jr}}$ 是非线性函数，所以使用Newton-Raphson法对${\kappa }_{jr}$ 进行估计，于是对上式的${\kappa }_{jr}$ 再次求偏导：

$\frac{\partial l_1}{\partial {\kappa }_{jr}\partial {\kappa }_{jr}}={\gamma }_{ijr}\left[\frac{1}{{\left({\kappa }_{jr}+1\right)}^2}-\frac{3}{{\kappa }_{jr}^4{\sigma }_{jr}^2}{\left(y_{ij}-x_{ij}\right)}^2\right]$ (19)

利用牛顿迭代关系式，假设迭代次数为t时${\kappa }_{jr}$ 的值为${\kappa }_{jr}^t$ ，在第$t+1$ 次迭代时有：

$\begin{array}{c}{\kappa }_{jr}^{t+1}={\kappa }_{jr}^t-\frac{{l^{\prime }}_1\left({\kappa }_{jr}^t\right)}{{l^{″}}_1\left({\kappa }_{jr}^t\right)}\\ ={\kappa }_{jr}^t-\frac{-\frac{1}{{\kappa }_{jr}^t+1}+\frac{1}{{\left({\kappa }_{jr}^t\right)}^3{\sigma }_{jr}^2}{\left(y_{ij}-x_{ij}\right)}^2}{\frac{1}{{\left({\kappa }_{jr}^t+1\right)}^2}-\frac{3}{{\left({\kappa }_{jr}^t\right)}^4{\sigma }_{jr}^2}{\left(y_{ij}-x_{ij}\right)}^2}\end{array}$ (20)

更新X：

根据(14)式，关于X的目标函数可以重写为：

$\begin{array}{c}{l}_2={\displaystyle \sum _{j=1}^B{\displaystyle \sum _{i=1}^{MN}{\displaystyle \sum _{r=1}^R{\gamma }_{ijr}\left[\frac{1}{2{\eta }_{ijr}{\sigma }_{jr}^2}{\left(y_{ij}-x_{ij}\right)}^2\right]}}}\\ ={\displaystyle \sum _{j=1}^B{\displaystyle \sum _{i=1}^{MN}\left({\displaystyle \sum _{r=1}^R{\gamma }_{ijr}\frac{1}{2{\eta }_{ijr}{\sigma }_{jr}^2}}\right){\left(y_{ij}-x_{ij}\right)}^2}}\\ ={\Vert W\odot \left(Y-X\right)\Vert }_F^2\end{array}$ (21)

其中，$\odot$ 表示元素乘积。对于W的元素，有：

$w_{ij}=\{\begin{array}{l}\sqrt{{\displaystyle \sum _{r=1}^R{\gamma }_{ijr}\frac{1}{2{\eta }_{ijr}{\sigma }_{jr}^2}}},i,j\in \Omega \\ 0,i,j\notin \Omega \end{array}$ (22)

其中，$w_{ij}$ 为W的元素。$\Omega$ 为是Y中非缺失项的集合。

最小化问题如下所示，可以使用DW算法来进行求解：

$\underset{U,V}{\min }{\Vert W\odot \left(Y-U{V}^{\text{T}}\right)\Vert }_F^2,\text{ s}\text{.t }X=U{V}^{\text{T}}$ (23)

算法1 基于AMoG的求解算法

输入：矩阵Y

初始化：$\Theta =\left\{\alpha ,\sigma ,\kappa \right\}$ ，因子矩阵$U,V$ ，混合分量数$R_{start}$ ，最大迭代$MaxIter$ ，容差Tol

1) 通过(13)更新$\gamma$ ；

2) 通过(15)更新$\sigma$ ；

3) 通过(16)更新$\alpha$ ；

4) 通过(20)更新$\kappa$ ；

5) 通过(23)更新$U,V$ ；

输出：参数$\Theta$ ，因子矩阵$U,V$ ，混合分量数$R_{final}$ 和后验概率$\gamma ={\left({\gamma }_{ijr}\right)}_{m\times n\times R_{final}}$

4. 实验

合成数据实验和真实场景内的数据实验是在一台拥有英特尔至强E5-2620v4CPU，2.1 GHz和31 GB内存的计算机上进行的，在Matlab R2019b中进行仿真。与MoG、PMoEP、NMoG、NGMeet、LRTFDFR [22]、BALMF六种方法进行对比。评估指标包括：峰值信噪比(PSNR)、结构相似性(SSIM)、特征相似性(FSIM)、相对无量纲全局误差(ERGAS)和光谱角制图(SAM)。

4.1. 合成数据实验

为了证明AMoG具有一定的优越性，在模拟实验中，选择了2个干净的遥感图像数据集进行合成数据实验。两个数据集分别是由反射光学系统成像光谱仪(ROSIS-03)捕获的帕维亚市中心(Pavia)数据集的子图像(PCC)，大小为200 × 200 × 80和子图像(PaviaU_A)大小为200 × 200 × 103。

表1模拟了不同类型的噪声对PCC和PaviaU_A的影响。考虑了六种情况，包括具有不同标准差$\sigma$ 的高斯噪声、具有不同强度p的脉冲噪声以及在所实现的高光谱图像的不同频带上具有不同比例$n_1$ 的死线噪声、不同比例$n_2$ 的条纹噪声的混合。

Table 1. Simulated noise setting

表1. 模拟噪声设置

情况	高斯噪声		脉冲噪声		死线噪声		条纹噪声
	波段	$\sigma$	波段	p	波段	$n_1$	波段	$n_2$
1	100%	0.05	--	--	--	--	--	--
2	100%	0.1	100%	0.01	--	--	--	--
3	100%	0.05	--	--	100%	3~10	--	--
4	100%	0.05	--	--	--	--	100%	3~10
5	100%	0.05	20%	0.01	20%	3~10	20%	10~30
6	100%	0.1	30%	0.01	30%	13~20	30%	20~50

表2和表3分别显示了PCC和PaviaU_A数据集上的6种情况的定量指标。最佳的结果用粗体字体突出显示。从此表中可以得到以下观察结果：

(1) 情况1只含有高斯噪声，虽然本文提出的AMoG不是专门为高斯噪声设计的，但在这种简单的理想情况下，它仍然可以获得与大多数竞争方法相比较的结果。

(2) 情况2~6包含了不同噪声类型的混合，使之与真实噪声更为相似。在这些情况下，AMoG实验结果表现较好。这可以合理地解释为，AMoG可以很好地模拟每个波段的重尾度和不对称性，且更适合对复杂的噪声进行建模。

图1和图2分别展示了情况4下各种去噪方法对PCC数据集的效果图(第60波段)以及情况5下各种去噪方法对PaviaU_A数据集的效果图(第103波段)。从图中可以看出所有方法都可以对噪声进行去除，但PMoEP和NMoG方法仍存在一些条纹噪声和截止日期噪声。MoG和LRTFDFR方法处理的结果过于平滑，细节信息不明显，可能会导致光谱失真。在被BALMF方法去噪的图像中，虽然能够显示出细节纹理，但会改变图像的对比度。相比之下，本文方法获得了更合理的视觉效果。

Table 2. Comparison of denoising effects of different algorithms on PCC datasets under different noise conditions

表2. 不同算法在不同噪声情况下对PCC数据集的去噪效果对比

情况	评价指标	Noisy	MoG	PMoEP	NMoG	LRTFDFR	BALMF	AMoG
1	PSNR	26.0150	39.2457	38.6649	39.2527	38.8334	38.0970	39.2595
	SSIM	0.7207	0.9784	0.9760	0.9745	0.9793	0.9720	0.9784
	FSIM	0.8628	0.9876	0.9865	0.9862	0.9858	0.9843	0.9878
	ERGAS	184.6058	41.3552	44.3153	51.8238	43.5660	47.2366	41.3155
	SAM	0.3375	0.0774	0.0830	0.0893	0.0633	0.0865	0.0771
2	PSNR	19.1585	31.8340	26.9529	33.6164	32.5079	32.8437	34.4109
	SSIM	0.4125	0.9087	0.7920	0.9360	0.9304	0.9138	0.9463
	FSIM	0.7308	0.9533	0.8913	0.9677	0.9546	0.9601	0.9785
	ERGAS	405.9134	98.6753	177.6485	77.5986	92.9199	87.6106	72.4463
	SAM	0.5241	0.1267	0.2079	0.0941	0.0946	0.1245	0.0889
3	PSNR	24.3388	36.1720	30.2209	38.0090	34.0096	37.8200	38.0139
	SSIM	0.6840	0.9658	0.9155	0.9767	0.9400	0.9707	0.9742
	FSIM	0.8456	0.9796	0.9496	0.9863	0.9554	0.9840	0.9860
	ERGAS	221.9761	63.6497	118.6285	47.5994	74.2050	49.1123	47.1643
	SAM	0.3823	0.1033	0.1548	0.0789	0.0792	0.0896	0.0731
4	PSNR	25.0288	36.1171	38.2768	37.8055	33.2894	37.9308	38.3930
	SSIM	0.6943	0.9666	0.9742	0.9702	0.9255	0.9708	0.9744
	FSIM	0.8488	0.9789	0.9849	0.9826	0.9427	0.9837	0.9851
	ERGAS	206.6823	63.8873	46.0881	50.4903	80.8614	48.2606	45.7388
	SAM	0.3659	0.1047	0.0846	0.0999	0.0844	0.0881	0.0844
5	PSNR	23.9304	35.9880	36.5621	37.3444	34.6567	37.9191	37.9237
	SSIM	0.6677	0.9612	0.9617	0.9639	0.9539	0.9711	0.9733
	FSIM	0.8415	0.9776	0.9796	0.9814	0.9681	0.9841	0.9851
	ERGAS	245.8752	69.1058	58.9979	54.9569	77.3999	48.6878	48.5230
	SAM	0.3950	0.0974	0.0963	0.0888	0.0877	0.0867	0.0742
6	PSNR	18.6694	32.6353	27.9668	32.1009	32.9470	32.4257	33.4051
	SSIM	0.3894	0.9132	0.8096	0.9181	0.9385	0.9078	0.9401
	FSIM	0.7155	0.9574	0.8981	0.9555	0.9590	0.9546	0.9680
	ERGAS	426.7157	87.2182	176.8810	98.7349	90.1368	91.1484	71.3959
	SAM	0.5542	0.1097	0.2879	0.1287	0.1118	0.1479	0.1034

Table 3. Comparison of denoising effects of different algorithms on PaviaU_A datasets under different noise conditions

表3. 不同算法在不同噪声情况下对PaviaU_A数据集的去噪效果对比

情况	评价指标	Noisy	MoG	PMoEP	NMoG	LRTFDFR	BALMF	AMoG
1	PSNR	26.0220	38.3396	38.3522	38.3051	37.8710	37.4874	38.3625
	SSIM	0.5808	0.9517	0.9518	0.9520	0.9503	0.9416	0.9523
	FSIM	0.8202	0.9800	0.9801	0.9801	0.9701	0.9766	0.9803
	ERGAS	199.2532	53.2327	53.2024	53.4001	53.1719	57.8561	53.1295
	SAM	0.2684	0.0655	0.0656	0.0658	0.0573	0.0699	0.0655
2	PSNR	19.1924	33.5029	27.2378	33.2650	33.3530	33.0357	33.7974
	SSIM	0.2952	0.8758	0.7089	0.8753	0.8940	0.8591	0.8772
	FSIM	0.6590	0.9540	0.8778	0.9541	0.9331	0.9481	0.9542
	ERGAS	434.8824	85.5763	205.1165	87.9298	87.2107	91.8153	84.8594
	SAM	0.4945	0.0965	0.1794	0.0968	0.0866	0.1012	0.0978
3	PSNR	24.5063	37.3400	36.2519	37.3345	34.5827	37.2348	37.3554
	SSIM	0.5489	0.9467	0.9416	0.9478	0.9108	0.9397	0.9484
	FSIM	0.8004	0.9776	0.9741	0.9780	0.9374	0.9760	0.9780
	ERGAS	232.3087	57.8956	63.8441	57.9977	74.5361	59.8328	58.6258
	SAM	0.3142	0.0685	0.0724	0.0689	0.0706	0.0714	0.0689
4	PSNR	25.0339	37.7050	37.7012	37.7659	34.0539	37.3287	37.9073
	SSIM	0.5540	0.9461	0.9466	0.9458	0.8954	0.9401	0.9490
	FSIM	0.8029	0.9772	0.9777	0.9775	0.9233	0.9760	0.9789
	ERGAS	223.7446	56.2408	56.3052	56.3277	78.4783	58.6367	55.1460
	SAM	0.2973	0.0686	0.0691	0.0689	0.0730	0.0708	0.0672
5	PSNR	24.0348	37.5971	33.9497	37.4438	37.4556	37.3090	37.8447
	SSIM	0.5339	0.9453	0.8969	0.9452	0.9482	0.9404	0.9489
	FSIM	0.7967	0.9775	0.9578	0.9779	0.9693	0.9761	0.9790
	ERGAS	261.2998	57.3987	89.6644	57.2622	55.5549	58.6010	55.1886
	SAM	0.3272	0.0687	0.1027	0.0698	0.0685	0.0709	0.0678
6	PSNR	19.2184	32.2288	26.7374	33.2121	33.0015	32.5792	33.4241
	SSIM	0.2864	0.8694	0.6923	0.8816	0.8846	0.8493	0.8872
	FSIM	0.6410	0.9485	0.8723	0.9518	0.9226	0.9436	0.9525
	ERGAS	431.9785	98.3936	210.1975	88.0881	90.3147	96.0397	88.0800
	SAM	0.5147	0.1125	0.1993	0.0981	0.0999	0.1071	0.0991

(a) 原始图像 (b) 噪声图像 (c) MoG (d) PMoEP

(e) NMoG (f) LRTFDFR (g) BALMF (h) AMoG

Figure 1. Comparison of PCC datasets by various denoising methods in Case 4 (Band 60)

图1. 情况4下各种去噪方法对PCC数据集的对比图(第60波段)

(a) 原始图像 (b) 噪声图像 (c) MoG (d) PMoEP

(e) NMoG (f) LRTFDFR (g) BALMF (h) AMoG

Figure 2. Comparison of the PaviaU_A dataset by various denoising methods in Case 5 (Band 103)

图2. 在情况5下各种去噪方法对PaviaU_A数据集的对比图(第103波段)

图3展示了PaviaU_A数据集在情况5下第103波段对应的误差图，观察到PMoEP方法细节保留最明显，这表明该方法去噪后仍存在部分噪声。MOG、NMoG、LRTFDFR方法仍存在少部分的高斯噪声和条纹噪声。相比于BALMF方法，本文方法对应的误差图更暗，更接近原始图像的误差图，这表明本文方法能够恢复更干净的结构和保存更多的局部细节。

(a) 原始图像 (b) 噪声图像 (c) MoG (d) PMoEP

(e) NMoG (f) LRTFDFR (g) BALMF (h) AMoG

Figure 3. PaviaU_A error plot of the dataset corresponding to Band 103 in Case 5

图3. PaviaU_A数据集在情况5下第103波段对应的误差图

(a) (b) (c) (d)

(e) (f) (g)

Figure 4. All algorithms in Case 5 PaviaU_A spectral curves on the pixels (188, 80) of the dataset. The horizontal axis is the band, and the vertical axis is the numeric value

图4. 所有算法在情况5中 PaviaU_A数据集的像素(188, 80)上的光谱曲线。横轴为波段，纵轴为数值

除了定量指标外，图4还显示了所有算法在情况5中PaviaU_A数据集的像素(188, 80)处的光谱曲线。地面真实值用棕色表示，去噪曲线用蓝色表示。观察到带噪图像受到强噪声的污染，PMoEP和BALMF方法不能很好地恢复真实谱。MoG、NMoG和LRTFDFR方法虽能够重建总体趋势，但仍存在剧烈的波动。作为比较，AGGD可以实现最接近真实谱的恢复，能够更好地拟合混合噪声。

4.2. 真实数据实验

真实数据实验的数据集选取的是由高光谱数字图像采集实验(HYDICE)得到的Urban数据集，大小为307 × 307 × 210，该数据集受到高斯噪声、截止日期噪声、条纹噪声、吸水率、大气和一些未知噪声的严重污染。本文选择了两个典型的频带来表示不同方法的恢复结果，图5显示了该数据中207波段的恢复结果。

(a) 噪声图像 (b) MoG (c) PMoEP (d) NMoG

(e) LRTFDFR (f) BALMF (g) AMoG

Figure 5. The denoising images of the Urban dataset under different denoising methods (Band 207)

图5. Urban数据集在不同去噪方法下的去噪图像(第207波段)

在图5的实验效果图中可以看到MoG和BALMF方法还存在大量噪声，PMoEP和LRTFDFR方法存在少量的条纹噪声，NMoG方法虽然能够很好地恢复出图像，但会改变图像的对比度。结果表明，本文提出的AMoG表现结果最好。

5. 结论

在传统的遥感去噪模型中，遥感图像噪声主要由对称分布建模，如高斯分布、拉普拉斯分布和高斯混合分布。但本文考虑到遥感图像噪声通常是不对称的，于是在高斯混合分布的基础上，将遥感图像噪声建模为非对称高斯混合分布，并提出了一个基于AMoG模型的遥感图像算法，从而灵活地拟合混合噪声。对于模型参数，使用EM算法进行求解。与具有代表性的遥感图像去噪算法相比，AMoG可以更好地对合成场景和真实场景的遥感图像进行去噪，表明了该算法具有一定的优越性。由于EM算法受初始值的影响，容易收敛至局部最优值，在后续工作中，可以考虑设计更优的算法推导求解模型参数。

基金项目

吉林省自然科学基金，NO. 20240101298JC；国家自然科学基金，NO. 12171054。