基于多元高斯混合模型的图像去噪算法

Image Denoising Algorithm Based on Multivariate Gaussian Mixture Model

摘要: 高光谱图像(Hyperspectral Images, HSIs)在遥感和医学成像等领域具有广泛的应用，但在采集过程中容易受到各种噪声的干扰。尽管目前已有多种去噪方法应用于高光谱图像处理，但这些方法在应对复杂光谱特征和复杂或非均匀噪声分布时仍面临挑战，且部分方法可能导致图像细节的丢失，降低光谱数据的真实性和有效性。为克服这些不足，文章提出了一种基于多元高斯混合模型的高光谱图像去噪方法。该方法将噪声建模为多元高斯混合模型，通过多模态特征表示图像的复杂光谱结构，以更好地适应不同光谱分布。我们采用变分贝叶斯(Variational Bayes, VB)方法进行参数估计，从而改善了传统期望最大化(EM)算法易于陷入局部最优的局限性，提高了参数估计的稳定性和模型收敛效率。实验结果表明，本文方法在多个数据集和多种噪声情况下均表现出优异的去噪效果，还能更好地保持光谱特性和图像结构的一致性，验证了其在高光谱图像去噪任务中的有效性。

Abstract: Hyperspectral Images (HSIs) have widespread applications in fields such as remote sensing and medical imaging, but they are often subject to various types of noise during the acquisition process. Despite the availability of numerous denoising methods for HSI processing, these methods still face challenges in handling complex spectral features and non-uniform noise distributions, which may lead to the loss of image details and compromise the authenticity and effectiveness of spectral data. To address these shortcomings, this paper proposes a hyperspectral image denoising method based on a multivariate Gaussian mixture model. The proposed method models noise using a multivariate Gaussian mixture model, employing multimodal features to represent the complex spectral structure of images, thereby better adapting to diverse spectral distributions. We adopt the Variational Bayes (VB) method for parameter estimation to overcome the limitations of the traditional Expectation-Maximization (EM) algorithm, which tends to get stuck in local optima, thus improving the stability of parameter estimation and convergence efficiency. Experimental results demonstrate that the proposed method achieves superior denoising performance across multiple datasets and various noise levels, while also better preserving spectral characteristics and image structure, validating its effectiveness for HSI denoising tasks.

1. 引言

高光谱图像具有丰富的光谱信息，在遥感、农业监测、地质勘探[1]和医学成像等领域被广泛应用。与传统成像方式不同，高光谱成像可以在连续多个光谱波段上捕捉目标的光谱特征，因此在物质分类、成分识别和目标检测等任务中具有显著优势。然而，由于成像系统的限制以及外界环境因素的影响，高光谱图像在获取过程中经常受到噪声干扰。噪声不仅会导致光谱信息的失真，还会影响后续数据分析和处理的准确性，增加分类和识别任务的难度。因此有效的去噪方法在高光谱数据处理中至关重要。

在高光谱图像去噪领域，构建有效的图像先验模型是一个重要的挑战。图像先验模型通过对图像特定结构的假设，指导去噪算法在保留图像重要细节的同时去除噪声。近年来，研究者基于高光谱图像的低秩性和稀疏性特征，提出了大量的去噪方法[2]-[5]。这些方法可大致分为基于矩阵的去噪方法和基于张量的去噪方法。

基于矩阵的低秩分解方法是假设高光谱图像数据位于一个低维子空间中，具有较强的光谱相关性。这类方法将高维图像数据矩阵化并展开，通过对图像矩阵执行低秩分解，去除噪声并恢复原始图像。例如，基于核范数的低秩分解方法将秩最小化问题转化为核范数最小化问题，通过惩罚核范数来逼近低秩解[6]-[9]。然而，这些方法仅能捕捉光谱维度上的相关性，忽略了高光谱图像的多维空间信息，对空间维度上的复杂结构缺乏充分的建模能力，在复杂噪声干扰下的去噪效果往往不理想。为了同时考虑光谱和空间维度的相关性，张量分解方法应运而生。与矩阵方法不同，张量方法将高光谱图像视作三维结构，通过张量分解将其建模为低秩张量以去噪。张量分解可以充分利用高光谱数据的多维结构信息，近年来的研究表明此类方法在保持高光谱图像的多线性结构方面有显著优势。例如，Tucker分解[10]和CP分解[11]通过不同的张量秩模型对图像的内在结构建模，Renard等人[12]提出了LR张量近似(LRTA)模型通过采用Tucker分解进行HSI去噪。Liu等人[13]设计了一个PARAFAC模型，利用平行因子分析对HSI去噪，Peng等人[14]提出了一种可分解的非局部张量字典学习(TDL)模型，该模型充分考虑了空间上的非局部相似性和光谱间的全局相关性。最近的方法试图结合不同低阶的优点张量模型，如Kronecker表示基于KBR的张量稀疏性度量[15]。在这些方法中，TDL和KBR在更一般的噪声HSI情况下取得了最先进的性能。

除了图像先验模型，噪声建模也是影响去噪算法性能的关键因素。传统的去噪方法通常假设噪声服从独立同分布(i.i.d.)的高斯[16]或拉普拉斯分布[17]，这种假设适用于简单噪声情况。然而，在实际应用中，噪声往往是多成分混合的，且具有更复杂的统计分布特性。为更贴近实际噪声分布，近年来，高斯混合模型(Gaussian Mixture Model, GMM)被引入到噪声建模中[18]，凭借其强大的连续概率密度近似能力，显著提升了去噪效果。基于GMM的噪声模型进一步发展为指数幂混合模型(MoEP) [19]，通过增强分布的灵活性，提高了去噪算法在复杂噪声场景中的适应性。然而，传统的GMM和MoEP模型仍普遍采用i.i.d.假设，这使得它们对于非独立同分布噪声时表现有限。

高光谱图像中的噪声分布在不同波段上可能表现出较强的不均匀性。针对这一特点，最近的一些研究通过非独立同分布(non-i.i.d.)假设增强了噪声建模的灵活性，如non-i.i.d. MOG [20]以及狄利克雷过程高斯混合[21] [22]，从而产生更好的噪声拟合能力和更高的去噪精度。Ma等[23]提出NMoG-Tucker将图像先验和噪声先验整合到一个全贝叶斯生成模型中，并设计了一种高效的变分贝叶斯算法。

基于上述研究，本文假设每个波段的噪声服从非独立同分布(non-i.i.d.)的多元高斯混合分布，以更灵活地适应复杂的噪声分布。我们采用变分贝叶斯方法对模型参数进行迭代更新，以提高参数估计的稳定性和收敛效率。通过在模拟数据和真实数据上的一系列实验，结果表明本文方法在去除高光谱图像的复杂混合噪声方面表现出了显著的优越性。

2. 相关工作

2.1. 高斯混合模型

2.1.1. 单高斯分布

单高斯分布(也称为正态分布)是最常见的概率分布之一，用于描述单一随机变量的分布。其概率密度函数(PDF)为：

$p\left(x\right)=\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp \left(-\frac{{\left(x-\mu \right)}^2}{2{\sigma }^2}\right)$ (1)

其中，$\mu$ 是均值，${\sigma }^2$ 是方差。高斯分布的形状是对称的钟形曲线，均值决定其中心位置，方差则控制分布的宽度。

2.1.2. 多元高斯分布

多元高斯分布是单高斯分布的扩展，用于描述多个随机变量的联合分布。设$x\in {ℝ}^D$ 为D维随机变量，其概率密度函数为：

$p\left(x\right)=\frac{1}{{\left(2\pi \right)}^{D/2}{\left|\text{Σ}\right|}^{1/2}}\exp \left(-\frac{1}{2}{\left(x-\mu \right)}^{\text{T}}{\text{Σ}}^{-1}\left(x-\mu \right)\right)$ (2)

其中，$\mu \in {ℝ}^D$ 是均值向量，$\Sigma \in {ℝ}^{D\times D}$ 是协方差矩阵，$\left|\Sigma \right|$ 是协方差矩阵的行列式，${\Sigma }^{-1}$ 是协方差矩阵的逆。

2.1.3. 高斯混合模型

高斯混合模型(GMM)是一种由多个高斯分布组成的概率模型，通常用于表示复杂的、多模态的数据分布。GMM的概率密度函数为：

$p\left(x\right)={\displaystyle \sum _{k=1}^K{\pi }_k\mathcal{N}\left(x|{\mu }_k,{\text{Σ}}_k\right)}$ (3)

其中，${\pi }_k$ 是第k个高斯成分的混合权重，满足${\displaystyle \sum _{k=1}^K{\pi }_k=1}$ ，$\mathcal{N}\left(x|{\mu }_k,{\text{Σ}}_k\right)$ 是第k个高斯分布的概率密度函数，具有均值${\mu }_k$ 和协方差矩阵${\Sigma }_k$ ，K是高斯分布的成分数。

2.2. 变分贝叶斯

贝叶斯推断是一种利用贝叶斯定理来估计后验分布的方法。给定观测数据D和模型参数$\theta$ ，贝叶斯推断的目标是计算后验分布：

$p\left(\theta |D\right)=\frac{p\left(D|\theta \right)p\left(\theta \right)}{p\left(D\right)}$ (4)

其中，$p\left(D|\theta \right)$ 是似然函数，$p\left(\theta \right)$ 是先验分布，$p\left(D\right)$ 是边际似然。但在实际应用中，由于后验分布通常难以直接计算，因此需要使用近似方法。

变分推断方法通过引入一个变分分布$q\left(\theta \right)$ ，并最小化变分分布与真实后验分布之间的Kullback-Leibler (KL)散度来近似后验分布，即：

$\text{KL}\left(q\left(\theta \right)\left|\right|p\left(\theta |D\right)\right)=-{\displaystyle \int q\left(\theta \right)\log \left(\frac{p\left(D|\theta \right)p\left(\theta \right)}{q\left(\theta \right)}\right)\text{d}\theta }$ (5)

3. 基于多元高斯混合模型的图像去噪算法

3.1. 模型规划

将一个给定的HSI表示为一个矩阵$Y\in {ℝ}^{MN\times B}$ ，其中MN和B分别表示HSI的空间维度和光谱维度。数据矩阵Y可以被建模为：

$Y=X+E$ (6)

其中，$X\in {ℝ}^{MN\times B}$ 表示干净的高光谱图像，$E\in {ℝ}^{MN\times B}$ 表示嵌入的噪声。

结合低秩矩阵分解(Low-Rank Matrix Factorization, LRMF)方法，数据矩阵Y可以进一步建模为：

$Y=U{V}^{\text{T}}+E$ (7)

其中，$U\in {ℝ}^{MN\times R}$ 和$V\in {ℝ}^{B\times R}$ 分别表示空间和光谱特征的低秩基矩阵，R是所选的秩。

在光谱维度上，高光谱图像的每个像素都包含数十甚至数百个波段的数据。这些波段间往往具有较强的相关性，为了深入理解各波段间的相关性，本文计算了各波段之间的相关系数，并通过热力图进行了可视化。如图1，50个不同波段噪声的相关系数热力图所示，选取Shanghai的真实数据集，在这个数据集里对随机选取的50个不同波段噪声绘制了相关系数的热力图，图中黄色部分和深蓝色部分表示波段间是具有强相关性的。由此可见，某些波段间所承载的光谱信息具有显著的相似性。因此，在处理和对HSI去噪时，需要一个能够同时考虑这些高维特征及波段间相关性的模型。

为此，本文选择了多元高斯混合模型对高光谱图像中的噪声进行建模。将每个像素点的多波段光谱信息视为一个高维变量，通过混合多个高斯分布对这些高维数据进行建模。

将$e_{n\cdot }$ 表示为位于噪声矩阵E的第n行的元素，如上所述，位于第n行的噪声分布可以建模为：

$p\left(e_{n\cdot }\right)={\displaystyle \sum _{n=1}^N{\displaystyle \sum _{k=1}^K{\pi }_{k}N\left(e_{n\cdot }|{\mu }_k,{\Lambda }_k^{-1}\right)}}$ (8)

其中，${\pi }_k$ 是混合系数，有${\pi }_k\ge 0$ 且${\displaystyle {\sum }_{k=1}^K{\pi }_k=1}$ ，K为高斯分量数，${\mu }_k$ 和${\Lambda }_k$ ，分别为第k个高斯分量的平均值和精度矩阵。

Figure 1. Heatmap of correlation coefficients for noise across 50 different bands

图1. 50个不同波段噪声的相关系数热力图

接下来，我们引入一个独立的高斯-Wishart先验分布，控制每个高斯分布的均值和精度，形式为：

$\begin{array}{c}p\left(\mu ,\Lambda \right)=p\left(\mu |\Lambda \right)p\left(\Lambda \right)\\ ={\displaystyle \prod _{k}N\left({\mu }_k|m_0,{\left({\beta }_0{\Lambda }_k\right)}^{-1}\right)W\left({\Lambda }_k|W_0,{\nu }_0\right)}\end{array}$ (9)

其中，$W(\cdot )$ 表示Wishart分布。这样，就可以对不同波段噪声分布之间的相关性进行合理的编码。对于每个观测$e_{n\cdot }$ ，我们有一个对应的潜在变量$z_n$ ，它是⼀个“1-of-K”的二值向量，元素为$z_{nk}$ 。我们将观测数据集记作$E=\left\{e_1,\cdots ,e_n\right\}$ ，潜在变量$Z=\left\{z_1,\cdots ,z_N\right\}$ ，再结合式(6)，我们可以等价地将式(8)重写为：

$\begin{array}{l}{y}_n~{\displaystyle \prod _{n=1}^N{\displaystyle \prod _{k=1}^{K}N{\left(y_n|x_n+{\mu }_k,{\text{Λ}}_k^{-1}\right)}^{z_{nk}}}}\\ z_{nk}~\text{Multinomial}\left(z_{nk}|{\pi }_k\right)\\ {\pi }_k=\text{Dir}\left({\pi }_k|{\alpha }_0\right)\end{array}$ (10)

其中，$\text{Multinomial}(\cdot )$ 和$\text{Dir}(\cdot )$ 分别表示多项式分布和狄利克雷分布。

对于大多数确定性的LRMF模型，矩阵X的秩r是固定的。通过将问题建模到特定的生成模型中，可以自适应地从数据中学习秩。具体来说，假设U和V的列是由高斯分布生成的。对于$l=1,2,\cdots R$ ，其中R是一个大于真实秩r的预设值，那么

$u_{\cdot l}~N\left(u_{\cdot l}|0,{\gamma }_l^{-1}{I}_N\right),v_{\cdot l}~N\left(v_{\cdot l}|0,{\gamma }_l^{-1}{I}_B\right)$ (11)

其中，$u_{\cdot l}$ 和$v_{\cdot l}$ 分别为u和v的第l列。$I_N\left(I_B\right)$ 表示$N\times N\left(B\times B\right)$ 单位矩阵。${\gamma }_l$ 为$u_{\cdot l}$ 和$v_{\cdot l}$ 的精度，先验如下：

${\gamma }_l~\text{Gam}\left({\gamma }_l|{\xi }_0,{\delta }_0\right)$ (12)

U，V的每一个列对$u_{\cdot l}$ 和$v_{\cdot l}$ 具有相同的稀疏性轮廓，具有通用精度变量${\gamma }_{l\cdot }$ 的特征。我们已经证明，这种建模可以使某些${\gamma }_{ls}$ 的精度值较高，从而能够自动进行X的低秩估计。结合式(7)~(12)，我们可以构建完整的贝叶斯模型。该模型的图形模型表示如图2所示。目标是推断所有相关变量的后验：

$p\left(U,V,Z,\mu ,\Lambda ,\pi ,\gamma |Y\right)$ (13)

其中$\mu =\left\{{\mu }_1,\cdots ,{\mu }_K\right\}$ 和$\Lambda =\left\{{\Lambda }_1,\cdots ,{\Lambda }_K\right\}$ 分别表示每个单高斯分布参数的均值和精度，$Z={\left\{z_{nk}\right\}}_{N\times K}$ ，$\pi =\left\{{\pi }_1,\cdots ,{\pi }_K\right\}$ ，$\gamma =\left({\gamma }_1,\cdots ,{\gamma }_R\right)$ 。

Figure 2. Graphical model

图2. 图形化模型

3.2. 变分推断

我们使用变分贝叶斯(VB)方法来进行后验推理。VB的目标是使用变分分布$q\left(\theta \right)$ 来近似真实的后验$p\left(\theta |D\right)$ ，其中$\theta$ 表示参数集，D表示观测数据。为了实现这一目标，我们需要解决以下优化问题：

$\underset{q\in ℂ}{\min }\text{KL}\left(q\left(\theta \right)\left|\right|p\left(\theta |D\right)\right)=-{\displaystyle \int q\left(\theta \right)\ln \left\{\frac{p\left(\theta |D\right)}{q\left(\theta \right)}\right\}\text{d}\theta }$ (14)

其中，$\text{KL}\left(q\left|\right|p\right)$ 表示两个分布q和p之间的KL散度，$ℂ$ 表示概率密度的约束集，使最小化易于处理。假设C是关于一些不相交群的分布族：$q\left(\theta \right)={\displaystyle \prod _i{q}_i\left({\theta }_i\right)}$ ，则可以得到闭型最优解：

$q_i^{\ast }\left({\theta }_i\right)=\frac{\exp \left\{{\langle \ln p\left(\theta ,D\right)\rangle }_{\theta \backslash {\theta }_i}\right\}}{{\displaystyle \int \exp \left\{{\langle \ln p\left(\theta ,D\right)\rangle }_{\theta \backslash {\theta }_i}\right\}}}$ (15)

其中，$\langle \cdot \rangle$ 表示期望，而$\theta \backslash {\theta }_i$ 则表示从集合$\theta$ 中去除了${\theta }_i$ 后的集合。

基于这样的表示方法，我们可以进一步推断出式(13)中涉及的所有因子分解分布。假设后验分布式(13)的近似具有如下的因式分解形式：

$q\left(U,V,Z,\mu ,\text{Λ,}\pi \text{,}\gamma \right)={\displaystyle \prod _{n}q\left(u_{n\cdot }\right)}{\displaystyle \prod _{j}q\left(v_{j\cdot }\right)}{\displaystyle \prod _{nk}q\left(z_{nk}\right)}\times {\displaystyle \prod _{k}q\left({\mu }_k,{\Lambda }_k\right)q\left({\pi }_k\right)}{\displaystyle \prod _{l}q\left({\gamma }_l\right)}$ (16)

其中，$u_{n\cdot }$ 和$v_{j\cdot }$ 分别表示矩阵U和V的第n行和第j行。根据公式(15)，我们可以为公式(16)中的每个组成部分推导出具体的推断方程。

1) 噪声分量估计：

首先列出式(16)中涉及的噪声分量的更新方程。各波段噪声的均值和精度的后验分布更新如下：

$q^{\ast }\left({\mu }_k,{\text{Λ}}_k\right)=N\left({\mu }_k|m_k,{\left({\beta }_k{\text{Λ}}_k\right)}^{-1}\right)W\left({\text{Λ}}_k|W_k,{\nu }_k\right)$ (17)

其中，上式中的参数可以通过以下方式从数据中计算出来：

${\beta }_k={\beta }_0+N_k\text{, }{\overline{c}}_k=\frac{1}{N_k}{\displaystyle \sum _{n=1}^N{r}_{nk}\left(y_n-u_{n\cdot }{V}^{\text{T}}\right)}$(18)

$m_k=\frac{1}{{\beta }_k}\left({\beta }_0{m}_0+N_k{\overline{c}}_k\right)$ (19)

$S_k=\frac{1}{N_k}{\displaystyle \sum _{n=1}^N{r}_{nk}}{\left(y_n-u_{n\cdot }{V}^{\text{T}}-{\overline{c}}_k\right)}^{\text{T}}\left(y_n-u_{n\cdot }{V}^{\text{T}}-{\overline{c}}_k\right)$(20)

$W_k^{-1}=W_0^{-1}+N_k{S}_k+\frac{{\beta }_0{N}_k}{{\beta }_k}{\left({\overline{c}}_k-m_0\right)}^{\text{T}}\left({\overline{c}}_k-m_0\right)$(21)

$v_k=v_0+N_k$ (22)

潜变量Z的更新方程为：

$q^{\ast }\left(Z\right)={\displaystyle \prod _{n=1}^N{\displaystyle \prod _{k=1}^K{r}_{nk}^{z_{nk}}}}$ (23)

其中所涉及的参数是由以下公式计算得出：

$\begin{array}{l}{r}_{nk}=\frac{{\rho }_{nk}}{{\displaystyle {\sum }_{j=1}^K{\rho }_{nj}}}\\ \ln {\rho }_{nk}=\langle \ln {\pi }_k\rangle +\frac{1}{2}\langle \ln \left|{\text{Λ}}_k\right|\rangle -\frac{D}{2}\ln \left(2\pi \right)\\ -\frac{1}{2}\langle \left(y_n-u_{n\cdot }{V}^{\text{T}}-{\mu }_k\right){\text{Λ}}_k{\left(y_n-u_{n\cdot }{V}^{\text{T}}-{\mu }_k\right)}^{\text{T}}\rangle \end{array}$ (24)

不难确定上式(24)中涉及到的三个期望的一般表达式为：

$\langle \ln {\pi }_k\rangle =\psi \left({\alpha }_k\right)-\psi \left({\displaystyle \sum _{i=1}^K{\alpha }_i}\right)$ (25)

$\langle \ln \left|{\Lambda }_k\right|\rangle ={\displaystyle {\sum }_{i=1}^D\psi \left(\frac{{\nu }_k+1-i}{2}\right)}+D\ln 2+\ln \left|W_k\right|$ (26)

$\begin{array}{l}\langle \left(y_n-u_{n\cdot }{V}^{\text{T}}-{\mu }_k\right){\text{Λ}}_k{\left(y_n-u_{n\cdot }{V}^{\text{T}}-{\mu }_k\right)}^{\text{T}}\rangle \\ =D{\beta }_k{}^{-1}+v_k\left(y_n-u_{n\cdot }{V}^{\text{T}}-m_k\right)W_k{\left(y_n-u_{n\cdot }{V}^{\text{T}}-m_k\right)}^{\text{T}}\end{array}$ (27)

同样，属于第k个高斯分量的混合比例${\pi }_k$ 的更新方程可以写为

$q^{\ast }\left(\pi \right)={\displaystyle \prod _{n=1}^K{\pi }_k^{\alpha -1}}$ (28)

其中，${\alpha }_k={\alpha }_0+N_k$ ，$N_k={\displaystyle \sum _{n=1}^N{r}_{nk}}$ 。

2) 低秩成分的估计：

在完成噪声成分的更新之后，我们接着估计低秩成分$u_{n\cdot }$ (对于$n=1,\cdots ,N$ )和$v_{j\cdot }$ (对于$j=1,\cdots ,B$ )的后验分布：

$q^{\ast }\left(u_{n\cdot }\right)=N\left(u_{n\cdot }|{\mu }_{u_{n\cdot }},{\Sigma }_{u_{n\cdot }}\right)$ (29)

其中，

$\begin{array}{l}{\mu }_{u_{n\cdot }}={\displaystyle \sum _k\langle z_{nk}\rangle \left(y_n-\langle {\mu }_k\rangle \right)\langle {\Lambda }_k\rangle \langle V\rangle {\Sigma }_{u_{n\cdot }}}\\ {\Sigma }_{u_{n\cdot }}={\left({\displaystyle \sum _k\langle z_{nk}\rangle }\langle V^{\text{T}}{\Lambda }_{k}V\rangle +\langle \Gamma \rangle \right)}^{-1}\end{array}$ (30)

同理可以推断出V的后验分布，在这一过程中，我们采用逐行推断的方法进行计算，即：

$q^{\ast }\left(v_{j\cdot }\right)=N\left(v_{j\cdot }|{\mu }_{v_{j\cdot }},{\Sigma }_{v_{j\cdot }}\right)$ (31)

其中，

$\begin{array}{l}{\mu }_{v_{j\cdot }}=\left[{\displaystyle \sum _{n,k}\langle z_{nk}\rangle {\langle {\Lambda }_k\rangle }_{jj}}\left(y_{nj}-\langle {\mu }_{kj}\rangle \right)\langle u_{n\cdot }\rangle \right]{\Sigma }_{v_{j\cdot }}\\ {\Sigma }_{v_{j\cdot }}={\left({\displaystyle \sum _{n,k}\langle z_{nk}\rangle {\langle {\Lambda }_k\rangle }_{jj}\langle u_{n\cdot }^{\text{T}}{u}_{n\cdot }\rangle }+\langle \Gamma \rangle \right)}^{-1}\end{array}$ (32)

其中$\Gamma =\text{diag}\left(\langle \gamma \rangle \right)$ ，对于控制矩阵U和V秩的${\gamma }_l$ ，其后验分布服从伽马分布，即：

$q^{\ast }\left({\gamma }_l\right)=\text{Gam}\left({\gamma }_l|{\xi }_l,{\delta }_l\right)$ (33)

其中，

$\begin{array}{l}{\xi }_l=\frac{N+B}{2}+{\xi }_0\\ {\delta }_l=\frac{1}{2}{\displaystyle \sum _i\langle u_{il}^2\rangle }+\frac{1}{2}{\displaystyle \sum _j\langle v_{jl}^2\rangle }+{\delta }_0\end{array}$ (34)

Table 1. Image denoising based on multivariate gaussian mixture model

表1. 基于多元高斯混合模型的图像去噪

算法1. 基于多元高斯混合模型的图像去噪

输入：矩阵Y

初始化：k，超参数$\left\{m_0,{\beta }_0,v_0,W_0,{\alpha }_0,{\xi }_0,{\delta }_0\right\}$ ，因子矩阵$U,V$ ，混合分量数R

令迭代次数t = 0，执行以下操作直到收敛

1) 通过式(24)~(28)更新噪声分量$Z^t$ ，${\pi }^t$ 的近似后验

2) 通过式(17)~(22)噪声分量${\mu }^t$ ，${\Lambda }^t$ 的近似后验

3) 通过式(30)和(32)更新低秩分量$U^t$ ，$V^t$ 的近似后验

4) 通过式(33)和(34)更新噪声分量${\gamma }^t$ 中参数的近似后验

5) $t=t+1$

输出：因子矩阵$U,V$

完整的算法流程总结在表1。

4. 实验

为了验证所提出方法的有效性和去噪性能，我们在合成数据和真实高光谱图像(HSI)上进行了广泛的实验评估，并与BM4D [24]、TDL、KBR三种方法进行对比。在实验过程中，所有对比方法的参数均设置为默认值或通过手动调整，以确保得到最优的去噪效果。所有实验均在Intel(R) Core(TM) i7-1065G7 CPU @ 1.30GHz(1.50 GHz)、16.0 GB内存的64位Windows 10操作系统上，使用MATLAB 2021a进行了一系列实验。在实验评估中，我们采用了峰值信噪比(PSNR)和结构相似性(SSIM)作为评价指标，用于量化去噪效果的保真度和结构保持效果。通过与现有方法的全面对比，验证了所提出方法在不同噪声条件下的稳健性和去噪性能。

4.1. 合成HSI去噪实验

在本实验中，采用了两种HSI：华盛顿特区商场(Washington DC Mall)，尺寸为1208 × 307 × 191和帕维亚大学数据(Pavia)，尺寸为310 × 250 × 87。在裁剪HSI的主要部分并删除一些明显的视觉污染光谱通道后，华盛顿特区购物中心和帕维亚大学数据分别调整为200 × 200 × 160和310 × 250 × 87，将HSI的灰度值归一化为[0, 1]。现实世界的HSI通常被几种不同类型的噪声污染，包括最常见的高斯噪声、脉冲噪声、死线和条纹。为了模拟这些真实的HSI噪声场景，我们在原始的HSI数据中添加了6种噪声，见表2，模拟噪声设置。

Table 2. Simulated noise settings

表2. 模拟噪声设置

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

表2的模拟噪声设置模拟了不同类型的噪声对WDC数据集和Pavia数据集的影响。实验共考虑了六种情况，包括具有不同标准差$\sigma$ 的高斯噪声、具有不同强度p的脉冲噪声，以及在实现的高光谱图像的不同波段上，不同比例的死线噪声(比例$n_1$ )和条纹噪声(比例$n_2$ )的混合。

表3为不同算法在不同噪声情况下对WDC数据集的去噪效果对比，表4为不同算法在不同噪声情况下对Pavia数据集的去噪效果对比，它们分别显示了WDC和Pavia数据集上6种情况的定量指标。最佳的结果用粗体字体突出显示。从结果中可以明显看出，除了情况1外，所提出方法在其他情况下均表现出显著的优势。这是因为在独立同分布高斯噪声的情况下，传统方法的噪声假设更符合基本条件，在其模型中利用了更多有效的高光谱图像先验信息，因此往往表现相对更好。而在这种情况下，虽然本文提出的方法不是专门为高斯噪声设计的，但在这种简单的理想情况下，它仍然可以获得与大多数竞争方法相比较的结果。

在更为复杂但更贴近实际应用的非独立同分布噪声的情况2~6中，我们所提出方法的优势尤为明显。其在不同噪声情况下对WDC和Pavia数据集的去噪效果显著优于其他对比方法，且MPSNR (均方峰值信噪比)和MSSIM (平均结构相似性)均取得了较高的分数，表明其在去噪能力和保留图像细节方面的卓越性能。这可以通过其更好的噪声拟合能力来解释，即我们的方法能够更准确地提取嵌入在高光谱图像中的噪声，从而自然地提高了高光谱图像的恢复效果。

Table 3. Comparison of denoising performance on the WDC dataset under different noise conditions for various algorithms

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

	noisy	BM4D	TDL	KBR	ours
i.i.d. Gaussian Noise
MPSNR	26.02	36.66	40.15	38.75	39.28
MSSIM	0.763	0.973	0.988	0.985	0.985
Non-i.i.d. Gaussian Noise
MPSNR	27.36	34.90	30.04	36.75	37.21
MSSIM	0.806	0.949	0.880	0.974	0.977
Gaussian + Stripe Noise
MPSNR	26.3543	32.55	28.02	34.99	36.54
MSSIM	0.783	0.918	0.844	0.963	0.977
Gaussian + Deadline Noise
MPSNR	26.16	32.49	27.82	35.40	35.71
MSSIM	0.778	0.926	0.851	0.960	0.971
Gaussian + Impulse Noise
MPSNR	26.05	32.63	27.92	35.63	35.94
MSSIM	0.775	0.920	0.851	0.968	0.972
Mixture Noise
MPSNR	24.58	30.04	25.22	32.53	33.01
MSSIM	0.723	0.866	0.747	0.924	0.955

Table 4. Comparison of denoising performance on the Pavia dataset under different noise conditions for various algorithms

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

	noisy	BM4D	TDL	KBR	ours
i.i.d. Gaussian Noise
MPSNR	26.02	37.73	40.94	40.06	38.32
MSSIM	0.570	0.944	0.977	0.975	0.953
Non-i.i.d. Gaussian Noise
MPSNR	29.34	39.51	33.17	40.03	40.58
MSSIM	0.705	0.952	0.848	0.965	0.975
Gaussian + Stripe Noise
MPSNR	28.26	33.63	28.99	35.08	37.62
MSSIM	0.697	0.848	0.730	0.897	0.943
Gaussian + Deadline Noise
MPSNR	27.81	35.17	29.77	36.57	37.94
MSSIM	0.669	0.922	0.776	0.943	0.953
Gaussian + Impulse Noise
MPSNR	27.21	35.29	29.92	37.06	37.17
MSSIM	0.650	0.914	0.767	0.944	0.945
Mixture Noise
MPSNR	25.08	29.68	26.02	32.40	34.65
MSSIM	0.568	0.746	0.609	0.829	0.899

为了更直观地展示各方法在去噪效果上的差异，我们将去噪后的图像可视化。图3为情况4下各种去噪方法对WDC数据集的对比图(第79波段)，是我们选取WDC数据集在情况4下的第79波段不同方法去噪后的对比图像。图4为情况4下各种去噪方法对Pavia数据集的对比图(第79波段)，是Pavia数据集在情况4下的第79波段不同方法去噪后的对比图像。

从图中可以看到，各方法在噪声去除及细节保留方面的差异明显。BM4D方法对噪声有一定的抑制效果，但在细节恢复上仍显不足，图像边缘部分仍存在较多残余噪声，且细节纹理较为模糊。TDL方法在一定程度上改善了细节的恢复能力，相比BM4D更能保留图像的空间结构，但在复杂区域的去噪效果仍不理想。KBR方法在去噪和细节保持上表现更优，能够有效减少噪声并保留更多细节，然而在某些区域噪声仍有残留，尤其是在图像高频细节部分，表现略显不足。相比之下，我们的方法在去除噪声的同时，显著保持了图像的细节和结构，能够更好地还原原始图像的纹理和边缘信息，表现出更强的去噪效果和细节保留能力。因此，我们的方法在复杂噪声下展现出了较好的稳健性和适应性，相比其他方法更能实现高质量的遥感图像去噪。

Figure 3. Comparison of various denoising methods on the WDC dataset under condition 4 (band 79)

图3. 情况4下各种去噪方法对WDC数据集的对比图(第79波段)

Figure 4. Comparison of various denoising methods on the Pavia dataset under condition 4 (band 79)

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

4.2. 真实数据实验

真实数据实验的数据集选取的是由高光谱数字图像采集实验(HYDICE)得到的Urban数据集，大小为307 × 307 × 210，该数据集受到高斯噪声、截止日期噪声、条纹噪声、吸水率、大气和一些未知噪声的严重污染。其中一些波段受到大气和水的严重污染。为了全面地评估所提方法在复杂噪声环境下的有效性和鲁棒性，我们在实验中使用了全部的210个波段数据，未对污染波段进行任何去除，以更加接近真实应用场景。

在对比分析中，我们特别选取了两个典型的波段(如207波段和第103波段)作为代表，对比了不同去噪算法的效果。这些波段的噪声特征具有典型性，能够更好地反映各个方法在复杂噪声下的表现。图5为Urban数据集在不同去噪方法下的去噪图像(第207波段)，图6为Urban数据集在不同去噪方法下的去噪图像(第103波段)，它们分别展示了各方法在Urban数据集第207波段和103波段的去噪效果。

从图中可以看出，BM4D和TDL在去噪后依然存在较为明显的噪声残留，尤其是在高噪声区域，图像的细节和纹理信息丢失严重，整体图像质量较差。BM4D基于4D块匹配和硬阈值滤波，能够在某种程度上去除噪声，但在面对复杂多样的噪声类型时仍显得力不从心，去噪效果受限。TDL则利用了张量字典学习方法，通过对空间上的非局部相似性和光谱维度的全局相关性进行建模，表现相对BM4D有所改善，但噪声残留依然明显，且在纹理恢复上有所欠缺。KBR方法在去噪与细节保留上相较于BM4D和TDL有进一步的改善。其采用了一种结合Tucker和CP分解的Kronecker基表示来度量张量的稀疏性，能够更好地保留多维结构信息。在实验结果中可以看到，KBR去除了更多的噪声，并且细节保持得较为完整。然而，KBR依然存在一定的噪声残留，且在部分纹理信息的还原方面表现不足，导致图像的部分区域依旧存在较多模糊。我们的方法去噪效果相对最佳，噪声去除彻底，同时能很好地还原图像细节和纹理信息，整体图像质量更高。

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

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

Figure 6. The denoising images of the Urban dataset under different denoising methods (band 103)

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

5. 结论

本文针对高光谱图像中常见的复杂噪声问题，提出了一种基于多元高斯混合模型的去噪方法，并结合变分贝叶斯(VB)框架来进行参数估计和模型推导。通过引入多模态特征和稀疏表示的方法，我们能够更好地捕捉图像的复杂光谱结构，并适应多种噪声类型及其不均匀性。实验结果表明，相较于现有的经典去噪方法，本文提出的方法在多种噪声环境下均具有优异的去噪性能，能够在有效去除噪声的同时保留更多的光谱细节和空间结构。在后续工作中，我们计划进一步优化去噪模型，增强其对超参数设定的自动适应能力，以提升模型的实际应用效率。

基金项目

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