基于OGS-HL的遥感图像混合噪声去除算法

Hybrid Noise Removal Algorithm of Remote Sensing Image Based on OGS-HL

摘要: 遥感图像在成像过程中容易受到混合噪声污染，包括高斯噪声、条纹噪声和脉冲噪声等。这些混合噪声降低了遥感图像的质量，限制了其后续应用。为了解决这一问题，首先，通过对遥感图像的梯度值进行统计分析，发现遥感图像的空间梯度值是符合重尾分布，因此，设计了关于遥感图像空间梯度值的OGS-HL正则项来去除混合噪声模型，该正则项不仅可以减少全变分带来的阶梯效应，而且还可以对图像的梯度值进行合理的稀疏表示；其次，针对条纹噪声，考虑其具有低秩性且使用核范数来约束，而稀疏噪声则具有全局稀疏分布，并且采用L₁范数来约束；最后，采用交替方向乘子法和Majorization-Minimization算法来求解所提出的模型。通过与现有的算法进行比较，结果表明我们提出的算法在去除高水平混合噪声方面具有良好的效果。

Abstract: Remote sensing images are easily contaminated by mixed noise in the imaging process, including Gaussian noise, fringe noise and impulse noise. These mixed noises degrade the quality of remote sensing images and limit their subsequent applications. In order to solve this problem, firstly, through the statistical analysis of the gradient value of the remote sensing image, it is found that the spatial gradient value of the remote sensing image conforms to the heavy-tailed distribution. Therefore, the OGS-HL regularization term about the spatial gradient value of the remote sensing image is designed to remove the mixed noise model. Moreover, the gradient value of the image can be reasonably sparse represented. Secondly, the stripe noise was considered to have low rank and was constrained by the nuclear norm, while the sparse noise had global sparse distribution and was constrained by the L₁ norm. Finally, the alternating direction method of multipliers and the Majorization-Minimization algorithm are used to solve the proposed model. By comparing with the existing algorithms, the results show that our proposed algorithm has a good effect in removing high-level mixed noise.

1. 引言

遥感图像包含了丰富的空间和光谱信息，因此广泛应用于城市规划、农业和地球观测等等。但是遥感图像在其形成、传输和记录的过程中容易受到各种各样的噪声污染，这些噪声使得遥感图像的空间纹理细节信息变得模糊不清，导致遥感图像难以获得较好的视觉效果，所以去除混合噪声是遥感图像能够得以后续应用的关键步骤。

在过去的十几年里，国内外学者在图像去噪和复原方面提出了许多的方法，大多数图像复原方法都是基于优化方法建立的。而这些去噪和复原方法依赖于对图像的先验信息假设和噪声分量的统计分析，从而获得有效地信噪分离。

自从全变分方法提出来之后，由于其具有保留图像的边缘和细节纹理信息等优点，又可以在去除噪声的同时保持图像的平滑性。受到许多学者的广泛使用。1992年，Rudin等人提出了经典的全变分方法 [1] 。2008年，Wang等人在全变分的思想上提出了一种新的半二次交替最小化算法 [2] ，不仅适用于L₁范数，而且适用于L₂范数，并提高了全变分算法的计算速度。但是继续随着深入研究，逐渐发现了全变分的缺点，即容易产生阶梯效应。2015年，He等人提出了一种基于全变分正则化低秩矩阵分解方法 [3] 。另外，2022年，由Yin等人提出了一种新的全变分模型，解决了全变分的缺点，该模型利用L₀范数作为数据保真度来去除脉冲噪声，同时利用重叠组稀疏性全变分(Overlapping group sparse total variation, OGS-TV)作为正则项来减少阶梯效应 [4] 。2022年，杨俊慈等人提出了分数阶全变分泊松图像去噪的快速线性化增广拉格朗日模型 [5] 。

使用全变分正则项约束是探索遥感图像空间信息的另一种有效方法，并与光谱信息相结合进行去噪。2015年，Aggarwal等人提出了基于空谱全变分的图像去噪模型 [6] ，并提出了SSTV正则项(Spatio-Spectral Total Variation, SSTV)，SSTV是通过计算光谱差分和空间差分得出的，在没有权重的情况下，该方法可以恢复理想的图像，但是当得到的观测图像被高水平噪声污染时，就会产生类似噪声的伪影。2016年，Chang等人提出各向异性光谱–空间全变分模型用于遥感图像去除条带 [7] ，并提出了ASSTV正则项(Anisotropic spectral-spatial total variation, ASSTV)，ASSTV考虑了空间和光谱相关性，但是会产生光谱的过度拟合。2020年，Takeyama等人提出了基于空谱混合正则项的图像复原方法 [8] ，并提出了HSSTV正则项(Hybrid Spatio-Spectral Total Variation, HSSTV)，为了减少光谱的过度拟合，在光谱差分上施加了权重约束。2020年，陈育群等人提出了在高阶交叠组稀疏正则项的基础上，并结合图像的二阶梯度信息来建立了高阶交叠组稀疏正则项的图像恢复模型 [9] 。

现有的大多数关于混合噪声去除模型的研究没有把条纹噪声作为一个独立的分量来考虑，而这些研究学者是将它归为稀疏噪声来处理，导致了令人不满意的结果。因此我们提出了一种去除混合噪声算法。

所提出的算法贡献可以总结如下：

1) 由于以往的混合噪声去除方法是将条纹噪声归为稀疏噪声，并没有考虑条纹噪声独特的结构，所以本文将条纹噪声作为一个独立的分量，利用其低秩特性，从而获得更好的恢复结果。

2) 设计了关于遥感图像的OGS-HL正则项的混合噪声去除模型。

3) 采用交替方向乘子法和Majorization-Minimization算法求解所提出的模型，从仿真和真实实验的结果来看，与现有的去噪方法相比，该模型在去除高水平混合噪声方面具有良好的效果。

本文的其余部分安排如下。第2节介绍了低秩约束方法研究工作和OGS-HL正则项。提出的去除混合噪声模型及优化求解在第3节中进行描述。第4节是进行实验设置，并对仿真和真实实验进行了描述和分析，对本文涉及的参数进行了分析。最后在第5节中得出结论。

2. 相关工作

2.1. 关于使用低秩约束方法的研究进展

对于条纹噪声的低秩先验假设是一种常见的图像噪声处理方式，利用信号中条纹噪声的低秩性质来降低噪声的影响，提高图像的恢复质量的同时会丢失图像的部分信息。2013年，Lu等人提出了一种新颖的图正则项低秩表示(Low-rank representation, LRR)去条带算法 [10] ，首先使用LRR技术进行去条带，其次，在目标函数中引入图正则项，建立优化模型，从实验结果表明，该方法能够去除条带噪声和获得清晰的图像。2016年，Chang等人提出了一种基于低秩的单张图像分解模型(Low-rank based single-image decomposition model, LRSID) [11] ，实现了干净图像与条纹噪声的完美分离，为了进一步利用了遥感图像的光谱信息，并将二维图像分解方法扩展到三维上，该三维图像分解方法在去除条纹噪声方面获得了良好的实验结果。2020年，石建飞等人提出了结合非局部先验性与加权核范数最小化的图像去噪算法 [12] 。2021年，Zhang等人提出了一种基于双低秩矩阵分解(Double low-rank matrix decomposition, DLR)的图像去除混合噪声方法 [13] ，提出了新的图像退化问题，并针对遥感图像和条纹噪声使用核范数约束，然后建立恢复模型，从实验的结果来看，该模型能够去除大多数混合噪声污染的情况。

2.2. 重叠组稀疏性的超拉普拉斯先验

目前大部分研究工作集中在处理高斯噪声下的图像恢复问题，而高斯分布不能很好地描述图像中尾部较重的脉冲噪声。事实上，带有重尾分布的噪声通常出现在低频大气信号中，并且可以使用超拉普拉斯分布来近似描述，而非高斯或拉普拉斯分布。

此外，在自然界中的图像梯度值是符合这种重尾分布特性。因此，我们将这种分布推广到遥感图像上，并使用超拉普拉斯分布来刻画四组不同尺寸的清晰遥感图像是否也符合重尾分布。通过研究这些重尾分布的特征，我们可以更好地理解和处理遥感图像中的噪声，并进一步改进图像恢复算法的性能。

如图1所示，在图1(a)中，我们测试了Indian Pines数据集1，并选择了该图像尺寸大小为145 × 145 × 220，我们发现当p = 0.7时，该图像的空间和光谱的梯度值曲线是接近于超拉普拉斯曲线。我们又测试了在p = 0.7时低维度的遥感图像也是否符合这种情况，如图1(b)、图1(c)所示，我们选择Pansharpening数据集进行实验2，这两个图像的尺寸为256 × 256 × 8，发现这两个遥感数据集的空间梯度值曲线是接近于超拉普拉斯曲线，而它的光谱域梯度值分布成一条曲折的曲线，不符合超拉普拉斯曲线。我们又增加遥感图像的尺寸进行实验，选择了Pavia University scene数据集3且该图像尺寸大小为610 × 340 × 103进行实验，如图1(d)所示，发现当p = 0.9时，该数据集的空间梯度值曲线是接近于超拉普拉斯曲线，而光谱域的梯度值是不符合。

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

Figure 1. The gradient results of different remote sensing images were statistically analyzed respectively (a) Indian Pines (b) Pansharpening data set 8 (c) Pansharpening data set 14 (d) Pavia University scene

图1. 对不同遥感图像的梯度结果分别进行统计分析 (a) Indian Pines (b) Pansharpening数据集第8个 (c) Pansharpening数据集第14个 (d) Pavia University scene

近年来，研究学者对超拉普拉斯(Hyper-Laplacian, HL)先验表现出广泛的兴趣，因为它能够比拉普拉斯或高斯先验更准确地描述真实图像梯度值的重尾分布情况。超拉普拉斯先验仅能约束图像梯度值是符合重尾分布，而无法揭示图像中像素之间的相关性，也不能减少全变分带来的阶梯效应。因此，我们采用了一种新的正则项，即重叠组稀疏性超拉普拉斯先验(Overlapping group sparse hyper-laplacian, OGS-HL) [14] 。该超拉普拉斯先验能够促进图像梯度值的自然稀疏性，并通过利用重叠组稀疏性的结构信息来缓解恢复图像中的伪影，同时在不同程度上增强像素结构内部的相关性，从而有助于更准确地恢复图像。因此，我们将介绍重叠组稀疏超拉普拉斯先验。

${\tilde{V}}_{\left(i,j\right),K}=\left[\begin{array}{cccc}{v}_{i-m_1,j-m_1}& v_{i-m_1,j-m_1+1}& \cdots & v_{i-m_1,j+m_2}\\ v_{i-m_1+1,j-m_1}& v_{i-m_1+1,j-m_1+1}& \cdots & v_{i-m_1+1,j+m_2}\\ ⋮& ⋮& \ddots & ⋮\\ v_{i+m_2,j-m_1}& v_{i+m_2,j-m_1+1}& \cdots & v_{i+m_2,j+m_2}\end{array}\right]\in {ℝ}^{K\times K}$ (1)

其中 $m_1=\left[\frac{K-1}{2}\right]$ 和 $m_2=\left[\frac{K}{2}\right]$ ，表示以像素 $\left(i,j\right)$ 为中心，大小为K × K (K为选择的像素块大小)。

重叠组稀疏性超拉普拉斯先验(OGS-HL)定义如下：

${\Phi }_{OH}\left(f\right)={\phi }_{OH}\left(D_{x}f\right)+{\phi }_{OH}\left(D_{y}f\right)$ (2)

其中 ${\phi }_{OH}\left(v\right)=\underset{i,j=1}{\overset{n}{{\displaystyle \sum }}}{\Vert {\left|v_{\left(i,j\right),K}\right|}^q\Vert }_2=\underset{i,j=1}{\overset{n}{{\displaystyle \sum }}}\sqrt{\underset{k_1,k_2=-m_1}{\overset{m_2}{{\displaystyle \sum }}}{\left|v_{i+k_1,j+k_2}\right|}^{2q}}$ ， $q\in \left(0,1\right)$ ， $D_x$ 表示水平差分算子， $D_y$ 表示垂直差

分算子。又由于公式(2)式中的超拉普拉斯是非凸的，引入到目标函数中，使得我们建立的优化模型求解变得缓慢且困难。所以我们使用Majorization-Minimization (MM)算法进行求解，它是一种迭代算法，用于解决公式(2)复杂的优化问题。关于Majorization-Minimization算法的定义如下：

$\underset{v\in {ℝ}^{n^2}}{\min }\left\{R\left(v\right)=\frac{1}{2}{\Vert v-v_0\Vert }_2^2+\beta {\phi }_{OH}\left(v\right)\right\}$ (3)

关于上式(3)最后的求解结果如下：

$v^{\left(k+1\right)}={\left[I_{n^2}+\beta q\text{Λ}{\left(v^{\left(k\right)}\right)}^{\text{T}}\text{Λ}\left(v^{\left(k\right)}\right)T\left(v^{\left(k\right)}\right)\right]}^{-1}{v}_0$ (4)

其中 $I_{n^2}$ 表示为 $n^2$ 大小的单位矩阵， $\text{Λ}\left(v\right)$ 表示矩阵 $v$ 可以通过2D卷积运算轻松计算得出， $T\left(v\right)=diag\left({\left|v\right|}^{2q-2}\right)$ 。对于上述公式(1-4)具体由来，有兴趣的读者可以参考文献 [14] 中的详细过程。

对于使用超拉普拉斯先验方法，国内外研究学者还结合了其它方法。2018年，余义斌等人利用模糊图像及卷积核来联合生成先验，同时利用卷积核的谱特性先验和清晰图像梯度域的超拉普拉斯先验联合建立图像盲去模糊模型 [15] ，比单独使用卷积核先验来建模更合理，该模型能够恢复出清晰的图像，且收敛速度更快。2019年，Cheng等人提出了一种新的基于反射和光照先验的非凸Retinex变分模型 [16] ，采用超拉普拉斯先验来表示反射率梯度，使用混合超拉普拉斯算子和Tikhonov来表示光照梯度。

3. 遥感图像恢复模型的建立及优化求解

一般而言，在遥感图像成像过程中，观测到的遥感图像被各种混合噪声破坏，其中就包括高斯噪声、脉冲噪声、死像素、死线和条纹噪声等等。在本篇论文中，我们考虑如下遥感图像观测模型：

$Y=X+S+B+N$ (5)

其中Y表示观测到的遥感图像，X表示干净的遥感图像，S表示条纹噪声(死线)，B表示稀疏噪声(脉冲噪声、死像素)，N表示高斯噪声， $Y、X、S、B、N\in R^{m\times n\times h}$ ，其中m表示遥感图像的宽度，n表示遥感图像的高度，h表示遥感图像光谱带的数量。我们假设公式(5)中不同的噪声是独立分布的。我们方法目的是从公式(5)中恢复出干净的遥感图像X。

正如我们在前面所分析的，通过在条纹噪声上的大量研究结果，在条纹噪声图像上，它是具有重复的结构，可以描述为低秩特征且使用核范数进行约束。为了从观测的遥感图像Y上恢复出干净的图像X，除了使用OGS-HL正则项外，还对稀疏噪声分别使用了不同的范数来约束。

在遥感图像恢复模型中，我们使用𝐹范数作为高斯噪声的数据保真项，并采用重叠组稀疏性的超拉普拉斯先验来约束遥感图像的空间梯度值。此外，考虑到遥感图像光谱域的局部连续性，我们还引入了L₁范数进行约束。针对条纹噪声和稀疏噪声分别采用核范数和全局稀疏分布的L₁范数进行约束。

$\begin{array}{l}\underset{X,S,B}{\min }\frac{1}{2}{\Vert Y-X-S-B\Vert }_F^2+{\lambda }_1{\phi }_{OH}\left(D_{x}X\right)+{\lambda }_2{\phi }_{OH}\left(D_{y}X\right)\\ +{\lambda }_3{\Vert D_{z}X\Vert }_1+{\lambda }_4{\displaystyle {\sum }_{b=1}^h{\Vert S_b\Vert }_{\ast }}+{\lambda }_5{\Vert B\Vert }_1\\ \text{s}\text{.t}rank\left(S_b\right)\le r_1\end{array}$ (6)

其中公式(6)中的条纹噪声 $S_b$ 从 $b=1$ 到h个波段的求和。由于 $D_{x}X,D_{y}X,D_{z}X$ 是不可微的和非线性的，所以引入辅助变量x₁、x₂和x₃进行替换。

$\begin{array}{l}\underset{X,S,B}{\min }\frac{1}{2}{\Vert Y-X-S-B\Vert }_F^2+{\lambda }_1{\phi }_{OH}\left(x_1\right)+{\lambda }_2{\phi }_{OH}\left(x_2\right)\\ +{\lambda }_3{\Vert x_3\Vert }_1+{\lambda }_4{\displaystyle {\sum }_{b=1}^h{\Vert S_b\Vert }_{\ast }}+{\lambda }_5{\Vert B\Vert }_1\\ \text{s}\text{.t}rank\left(S_b\right)\le r_1{D}_{x}X=x_1,D_{y}X=x_2,D_{z}X=x_3\end{array}$ (7)

对于公式(7)约束形式的最小化问题的求解，我们采用交替方向乘子法，接下来，我们将介绍整个优化过程，并根据增广拉格朗日乘子法来求解优化问题(7)，它的增广拉格朗日函数如下：

$\begin{array}{c}{\mathcal{L}}_{\delta }\left(X,S,B,M,W\right)=\frac{1}{2}{\Vert Y-X-S-B\Vert }_F^2+{\lambda }_1{\phi }_{OH}\left(x_1\right)+{\lambda }_2{\phi }_{OH}\left(x_2\right)+{\lambda }_3{\Vert x_3\Vert }_1\\ +{\lambda }_4{\displaystyle {\sum }_{b=1}^h{\Vert S_b\Vert }_{\ast }}+{\lambda }_5{\Vert B\Vert }_1+\frac{{\delta }_1}{2}{\Vert D_{x}X-x_1+\frac{w_1}{{\delta }_1}\Vert }_F^2\\ +\frac{{\delta }_2}{2}{\Vert D_{y}X-x_2+\frac{w_2}{{\delta }_2}\Vert }_F^2+\frac{{\delta }_3}{2}{\Vert D_{z}X-x_3+\frac{w_3}{{\delta }_3}\Vert }_F^2\end{array}$ (8)

其中 $W=\left(w_1,w_2,w_3\right)$ 为拉格朗日乘子， $\delta =\left({\delta }_1,{\delta }_2,{\delta }_3\right)$ 为惩罚参数， $M=\left(x_1,x_2,x_3\right)$ 。

3.1. M子问题优化求解

$\begin{array}{l}{x}_1=\arg {\min }_{x_1}\left\{{\lambda }_1{\phi }_{OH}\left(x_1\right)+\frac{{\delta }_1}{2}{\Vert D_{x}X-x_1+\frac{w_1}{{\delta }_1}\Vert }_F^2\right\}\\ x_2=\arg {\min }_{x_2}\left\{{\lambda }_2{\phi }_{OH}\left(x_2\right)+\frac{{\delta }_2}{2}{\Vert D_{y}X-x_2+\frac{w_2}{{\delta }_2}\Vert }_F^2\right\}\end{array}$ (9)

对于公式(9)的求解，我们可以用Majorization-Minimization算法求解上面问题的解，求解的结果如下：

$\begin{array}{l}{x}_1^{\left(k+1\right)}={\left[I_{n^2}+\frac{{\lambda }_1}{{\delta }_1}q\text{Λ}{\left(x_1^{\left(k\right)}\right)}^{\text{T}}\text{Λ}\left(x_1^{\left(k\right)}\right)T\left(x_1^{\left(k\right)}\right)\right]}^{-1}\left(D_x{X}^{\left(k\right)}+\frac{w_1^{\left(k\right)}}{{\delta }_1}\right)\\ x_2^{\left(k+1\right)}={\left[I_{n^2}+\frac{{\lambda }_2}{{\delta }_2}q\text{Λ}{\left(x_2^{\left(k\right)}\right)}^{\text{T}}\text{Λ}\left(x_2^{\left(k\right)}\right)T\left(x_2^{\left(k\right)}\right)\right]}^{-1}\left(D_y{X}^{\left(k\right)}+\frac{w_2^{\left(k\right)}}{{\delta }_2}\right)\end{array}$ (10)

接着求解x₃子问题，与x₃相关的子问题如下：

$x_3=\arg {\min }_{x_3}\left\{{\lambda }_3{\Vert x_3\Vert }_1+\frac{{\delta }_3}{2}{\Vert D_{z}X-x_3+\frac{W_3}{{\delta }_3}\Vert }_F^2\right\}$ (11)

我们对(11)式使用软阈值进行求解，能够得出如下形式：

$x_3^{\left(k+1\right)}=soft\left(D_z{X}^k+\frac{w_3^{\left(k\right)}}{{\delta }_3},\frac{{\lambda }_3}{{\delta }_3}\right)$ (12)

其中 $soft\left(a,b\right)=\{\begin{array}{ll}a-b,\hfill & a>b\hfill \\ a+b,\hfill & a<b\hfill \\ 0,\hfill & \text{other}\hfill \end{array}$

3.2. X子问题优化求解

给定的观测遥感图像Y、条纹噪声S和稀疏噪声B，干净的遥感图像X优化方程如下：

$\begin{array}{c}X=\arg {\min }_X\frac{1}{2}{\Vert Y-X-S-B\Vert }_F^2+\frac{{\delta }_1}{2}{\Vert D_{x}X-x_1+\frac{w_1}{{\delta }_1}\Vert }_F^2\\ +\frac{{\delta }_2}{2}{\Vert D_{y}X-x_2+\frac{w_2}{{\delta }_2}\Vert }_F^2+\frac{{\delta }_3}{2}{\Vert D_{z}X-x_3+\frac{w_3}{{\delta }_3}\Vert }_F^2\end{array}$ (13)

快速傅里叶变换可以有效地解决公式(13)，可以得出X子问题的解，结果如下式子：

$\begin{array}{c}{X}^{\left(k+1\right)}={\mathcal{F}}^{-1}\left(\mathcal{F}\left(\frac{Q}{R}\right)\right)\end{array}$ (14)

其中 $Q={\delta }_1{D}_x^{\text{T}}\left(x_1^{\left(k+1\right)}-\frac{w_1^{\left(k\right)}}{{\delta }_1}\right)+{\delta }_2{D}_y^{\text{T}}\left(x_2^{\left(k+1\right)}-\frac{w_2^{\left(k\right)}}{{\delta }_2}\right)+{\delta }_3{D}_z^{\text{T}}\left(x_3^{\left(k+1\right)}-\frac{w_3^{\left(k\right)}}{{\delta }_3}\right)+Y-S^{\left(k\right)}-B^{\left(k\right)}$ 和 $R=1+{\delta }_1{D}_x^{\text{T}}{D}_x+{\delta }_2{D}_y^{\text{T}}{D}_y+{\delta }_3{D}_z^{\text{T}}{D}_z$ ， $\mathcal{F}$ 为傅里叶变换， ${\mathcal{F}}^{-1}$ 为傅里叶逆变换。

3.3. S子问题优化求解

给定的观测遥感图像Y、干净的遥感图像X和稀疏噪声B，条纹噪声S优化方程如下：

$S=\arg {\min }_S\left\{{\lambda }_4{\displaystyle {\sum }_{b=1}^h{\Vert S_b\Vert }_{\ast }}+\frac{1}{2}{\Vert Y-X-S-B\Vert }_F^2\right\}$ (15)

条纹噪声S的更新是一个低秩矩阵逼近问题，可以通过对矩阵奇异值分解进行软阈值操作来轻松解决，可以得出条纹噪声S子问题的解，如下形式：

$\{\begin{array}{l}S=U\left(shin{k}_{-}{L}_{\ast }\left(\underset{r_1}{{\displaystyle \sum }},{\lambda }_4\right)\right)V^{\ast }\\ shin{k}_{-}{L}_{\ast }\left(\underset{r_1}{{\displaystyle \sum }},{\lambda }_4\right)=diag\left\{max\left(\underset{r_1,ii}{{\displaystyle \sum }}-{\lambda }_4,0\right)\right\}\\ SVD\left(Y-X-B\right)=U\sum V^{\ast }\\ \underset{r_1,ii}{{\displaystyle \sum }}=\underset{r_1}{{\displaystyle \sum }}=diag\left({\delta }_i\right),1<i<r_1\end{array}$ (16)

3.4. B子问题优化求解

给定的观测遥感图像Y、干净的遥感图像X和条纹噪声S，与稀疏噪声B子问题相关的子问题优化方程如下：

$B=\arg {\min }_B\left\{{\lambda }_5{B}_1+\frac{1}{2}{\Vert Y-X-S-B\Vert }_F^2\right\}$ (17)

同样我们对(17)式使用软阈值方法进行求解，可以得出B子问题的解，与(12)式是一致的。

$B^{\left(k+1\right)}=soft\left(Y-X^{\left(k+1\right)}-S^{\left(k+1\right)},{\lambda }_5\right)$ (18)

3.5. W子问题优化求解

最后，拉格朗日乘子W可以通过如下式子并行更新。

$\begin{array}{l}{w}_1^{\left(k+1\right)}=w_1^{\left(k\right)}+{\delta }_1\left(D_x{X}^{\left(k+1\right)}-x_1^{\left(k+1\right)}\right)\\ w_2^{\left(k+1\right)}=w_2^{\left(k\right)}+{\delta }_2\left(D_y{X}^{\left(k+1\right)}-x_2^{\left(k+1\right)}\right)\\ w_3^{\left(k+1\right)}=w_3^{\left(k\right)}+{\delta }_3\left(D_z{X}^{\left(k+1\right)}-x_3^{\left(k+1\right)}\right)\end{array}$ (19)

算法的流程总结在算法1中。

4. 数值实验

4.1. 实验的数据集

在仿真实验中，我们使用了两个遥感图像数据集，即Indian Pines数据集和Washington DC Mall数据集⁴。在Indian Pines数据集上，选择了该图像尺寸大小分别为145 × 145 × 50和145 × 145 × 10，50个波段指的是从该数据集的第一个波段到第五十个波段，10个波段同样如此。而Washington DC Mall数据集，选择了该图像尺寸大小为256 × 256 × 10。对于真实实验，采用大小为307 × 307 × 210的HYDICE Urban Dataset数据集5。并选择了201~210波段进行计算。在实验之前，将遥感图像的像素值归一化为[0, 1]范围内。

4.2. 实验的指标及设置

在仿真实验中，我们使用平均峰值信噪比(MPSNR)和平均结构相似性(MSSIM)作为图像恢复结果的定量评价性指标。MPSNR和MSSIM是图像处理和计算机视觉领域中应用最为广泛的两个质量评价性指标。MPSNR根据均方误差衡量复原图像的质量，MSSIM依据复原图像与参考图像之间的相似度。MPSNR和MSSIM的值越高，恢复的图像效果越好。在实验结果中，我们将最好的实验结果进行了加粗显示。

在真实的实验中，由于缺乏参考图像，我们没有使用评价性指标来衡量真实实验恢复图像的质量，通过视觉来评论恢复的效果。

为了验证所提出算法的有效性。选择了与我们模型中比较近似的算法进行比较，并选择如下几种算法，例如：LRMR [17] 、LRTV [3] 、LRTDGS [18] 、LRSID [11] 、ASSTV [7] 、L1HyMixDe [19] 、FGSLR_2/3 [20] 、LMHTV [21] 。对于LRSID算法，它是基于单波段图像条纹噪声的去除算法，在文献 [11] 中，它将二维图像条纹噪声的去除算法扩展到三维图像上，即LRMID算法，在后面所有的实验中，我们使用的都是LRMID算法进行比较。而LRTDGS算法，它是基于张量的基础上建立的恢复模型，使用张量方法恢复的图像能够加强图像的空间和光谱的相关性，避免恢复之后图像信息的丢失。在仿真实验中，我们根据该算法所在的文章中所提供参数范围进行相应的设置，将以上几种算法的参数调至最优，并在该参数范围的基础上扩展了一下参数的范围。在真实的实验中，我们通过大量仿真实验结果的参数来调节真实实验的各个算法的参数。

在所有的仿真实验中，对于条纹噪声，主要是考虑周期性垂直条纹噪声，通过逐个波段添加该噪声，条纹的密度设置为50%，即在该波段上选择空间尺寸的50%来添加条纹数量。条纹的周期设置为10，即添加的条纹与条纹之间的间隔距离为10。条纹强度是变化的，即在该波段上添加条纹数值大小(记作S)。在添加的条纹噪声的基础上，将高斯噪声与条纹噪声进行相加，通过逐个波段添加高斯噪声，添加的是零均值，它的标准差是发生变化的(记作G)，利用相同标准差乘一个matlab中自带的randn函数来产生高斯噪声，而randn函数根据选择二维图像尺寸大小来生成高斯分布。在添加的前两种噪声的基础上，再与脉冲噪声进行加和，并使用matlab中自带的imnoise函数给图像逐个波段添加脉冲噪声，添加的百分比是变化的(记作B)。例如：G = 0.05、S = 10和B = 0.05，添加高斯噪声标准差水平为0.05、条纹强度为10和脉冲噪声百分比水平为0.05。

由于其它方法的遥感图像退化问题只考虑两种噪声，而LRSID只考虑条纹噪声的去除，并且我们的方法考虑了三种噪声，所以我们分别设置以下三种实验情况进行分析与讨论。

例1：针对Indian Pines数据集进行实验，选择了图像大小为145 × 145 × 50。主要做了两个实验，第一个考虑高斯噪声和条纹噪声，它的高斯噪声水平G = 0.05和条纹强度S = 10。第二个考虑高斯噪声、条纹噪声和脉冲噪声，它的高斯噪声水平G = 0.025、条纹强度S = 10和脉冲噪声B = 0.05。

例2：主要针对Washington DC Mall数据集进行实验，主要是考虑高斯噪声、条纹噪声和脉冲噪声。并设置高斯噪声水平G = 0.025、条纹噪声强度和脉冲噪声水平 $B=\left[0.05,0.1,0.15\right]$ 。

通过对不同噪声水平下大量反复参数的计算实验，并将参数调至最优。最后将本文涉及的参数范围确定如下，对于 ${\lambda }_1$ 、 ${\lambda }_2$ 、 ${\lambda }_3\in \left[{10}^{-1},{10}^{-2},{10}^{-3},{10}^{-4}\right]$ ， ${\lambda }_4\in \left[{10}^{-1},{10}^{-2},{10}^{-3}\right]$ ， ${\lambda }_5\in \left[{10}^{-2},{10}^{-3}\right]$ ， $\delta \in \left[{10}^{-1},{10}^{-2}\right]$ ， $r_1\le 1$ ，对于公式(10)的OGS-HL中 $q\in \left(0,1\right)$ ， $K\in \left(0,10\right)$ ， $Nit\in \left(1,10\right)$ ， $Nit$ 是OGS-HL正则项迭代的次数。终止准则 和算法迭代次数分别设为10⁻⁵和200。

通过大量的仿真和真实的遥感数据实验来验证基于重叠组稀疏性的超拉普拉斯先验的遥感图像混合噪声去除方法的有效性，该所有的数值实验在Matlab R2021a (笔记本Intel(R) Core(TM) i5-8250U CPU @1.60GHz 1.80 GHz，8 GB内存)上运行。

4.3. 实验结果

我们测试了遥感图像在不同波段上的去除混合噪声的性能，图2显示了在去除两种噪声和三种噪声后Indian Pines数据集上不同波段上的PSNR和SSIM值，并比较六种算法的性能。在图2(a)、图2(c)中显然我们提出的算法在某一些波段上的PSNR值是远远高于其它五种算法。此外，在图2(b)、图2(d)中，我们对于曲线变化较为显著的区域进行了放大观察。通过观测发现，在图像某些波段上，我们提出的算法SSIM值均高于其它五种算法。这些实验证明了我们所提出算法在某些波段上具有很好的效果与优越性能。

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

Figure 2. The peak signal-to-noise ratio (PSNR) and structural similarity (SSIM) values of Indian Pines data set under different mixed noise levels were analyzed, (a, b) G = 0.05 and S = 10 (c, d) G = 0.025, B = 0.05 and S = 10

图2. 对Indian Pines数据集在不同混合噪声水平下各波段的峰值信噪比(PSNR)和结构相似度(SSIM)值进行了分析，(a, b) G = 0.05和S = 10 (c, d) G = 0.025、B = 0.05和S = 10

我们又在Washington DC Mall数据集和Indian Pines数据集上添加了G = 0.075和S = 10的混合噪声。图3和图4分别显示了Washington DC Mall数据集和Indian Pines数据集在去除高斯噪声和条纹噪声的实验结果。通过对比两个数据集的实验结果，发现六种算法在去除两种混合噪声方面有所不同。对于Washington DC Mall数据集，所提出的方法与LRTV、ASSTV和LRSID等方法都有良好的视觉效果，而LRMR和LRTDGS方法在图像恢复时有较多残留的条纹印记。另外，对于Indian Pines数据集，LRMR方法存在大量条纹印记和稀疏噪声，LRTV和LRTDGS方法有少量条纹印记，ASSTV和LRSID方法则有大量稀疏噪声，相比之下，所提出的方法在去除这两种噪声方面表现出良好的视觉效果。综上，所提出的方法与LRTV、ASSTV和LRSID等方法在去除高斯噪声和条纹噪声方面效果良好，但LRMR和LRTDGS方法在某些情况下存在局限性。

为了充分证明我们提出算法的性能，我们在Washington DC Mall数据集上分别添加三种不同混合噪声水平并计算出MPSNR和MSSIM值。表1给出了三种不同混合噪声水平下，针对Washington DC Mall数据集的所有比较方法的两个定量评价性指标值。实验结果表明，随着脉冲噪声的添加，在G = 0.025、S = 10和B = 0.05时，我们的算法恢复的MPSNR值达到了33.7682，优于其它五种方法。同时，在G = 0.025、S = 50和B = 0.15时，我们算法计算出来的数值也更好。从各个方法恢复结果来看，当添加脉冲噪声水平增大时，其它方法在去除高水平噪声方面显示不足之处，而我们的算法则具有此优点。

Figure 3. The Washington DC Mall dataset removes two kinds of noise resulting images, (a) clean image (b) noisy image (G = 0.075, S = 10) (c) LRMR (d) LRTV (e) LRSID (f) ASSTV (g) LRTDGS (h) OUR

图3. Washington DC Mall数据集去除两种噪声结果图像，(a) 干净图像 (b) 噪声图像(G = 0.075, S = 10) (c) LRMR (d) LRTV (e) LRSID (f) ASSTV (g) LRTDGS (h) OUR

Figure 4. The Indian Pines dataset removed two noise result images, (a) clean image (145 × 145 × 10), (b) noisy image (G = 0.075, S = 10), (c) LRMR (d) LRTV (e) LRSID (f) ASSTV (g) LRTDGS (h) OUR

图4. Indian Pines数据集去除两种噪声结果图像，(a)干净图像(145 × 145 × 10) (b) 噪声图像(G = 0.075, S = 10) (c) LRMR (d) LRTV (e) LRSID (f) ASSTV (g) LRTDGS (h) OUR

Table 1. Six algorithms recover two evaluative index values of the results under three different mixed noise levels

表1. 六种算法在三种不同混合噪声水平下恢复结果的两个评价性指标值

最后，我们又继续在Washington DC Mall数据集和Indian Pines数据集上添加G = 0.025、S = 10和B = 0.1混合噪声。图5和图6分别显示了Washington DC Mall数据集和Indian Pines数据集在去除高斯噪声、脉冲噪声和条纹噪声方面的结果。六种算法在去除三种混合噪声方面有不同表现。对于Washington DC Mall数据集，LRMR和LRTV方法仍有部分条纹噪声残留，而LRMR方法还有稀疏噪声，LRSID方法去除了条纹噪声，但存在大量稀疏噪声，影响了图像的视觉效果，ASSTV方法有少许稀疏噪声，相比之下，我们提出的算法有显著的视觉效果。而对于Indian Pines数据集，LRMR方法仍有大量条纹印记和稀疏噪声，LRTV和LRTDGS方法有少量条纹印记残留，LRSID方法去除了条纹噪声，但处理稀疏噪声不佳，导致局部模糊，ASSTV方法中也有一些稀疏噪声。相比之下，我们提出的算法具有良好的视觉效果。综上，我们提出的算法在去除高斯噪声、脉冲噪声和条纹噪声具有良好的效果。

另外我们验证了所提出的算法在单一噪声方面的去除能力，而LRSID是去除条纹噪声，并与最新提出的算法进行比较，例如：L1HyMixDe、FGSLR_2/3和LMHTV。所以我们在Indian Pines数据集上分别添加条纹强度 $m=0.3$ 和不同条纹密度 $r=\left[0.3,0.2\right]$ 进行实验。如表2所示，通过两组非周期垂直条纹噪声实验，发现我们提出的算法在条纹噪声去除方面也表现良好，且这两个评价性指标值优于其它四种算法，并且发现FGSLR_2/3算法在Washington DC Mall数据集上去除单一的条纹噪声方面也是具有很好的实验结果。

Figure 5. Washington DC Mall data set removes three kinds of noise resulting images, (a) Clean image (b) Noisy image (G = 0.025, S = 10, B = 0.1) (c) LRMR (d) LRTV (e) LRSID (f) ASSTV (g) LRTDGS (h) OUR

图5. Washington DC Mall数据集去除三种噪声结果图像，(a) 干净图像 (b) 噪声图像(G = 0.025, S = 10, B = 0.1) (c) LRMR (d) LRTV (e) LRSID (f) ASSTV (g) LRTDGS (h) OUR

Figure 6. Indian Pines data set removes three noise resulting images, (a) Clean image (145 × 145 × 10) (b) Noisy image (G = 0.025, S = 10, B = 0.1) (c) LRMR (d) LRTV (e) LRSID (f) ASSTV (g) LRTDGS (h) OUR

图6. Indian Pines数据集去除三种噪声结果图像，(a) 干净图像(145 × 145 × 10) (b) 噪声图像(G = 0.025, S = 10, B = 0.1) (c) LRMR (d) LRTV (e) LRSID (f) ASSTV (g) LRTDGS (h) OUR

我们在HYDICE Urban数据集上评估所提出算法在真实遥感图像上的性能。图7展示了波段206的Urban数据集的去噪结果。从视觉效果来看，LRTV、LMHTV和LRMR方法恢复后的图像颜色较暗，且LRMR方法留下大量水平和垂直方向条纹印记，而LRTV方法只有少许垂直条纹印记。LRSID方法恢复的图像上残留大部分稀疏噪声。LRTDGS恢复后图像颜色也显得暗淡。ASSTV方法能够有效去除条纹和稀疏噪声，但在处理细节时部分表现不佳。相比其它五种算法，我们提出的算法整体效果更好。

Table 2. MPSNR and MSSIM values of remote sensing images under different aperiodic fringe noise levels by five algorithms

表2. 五种算法在不同的非周期条纹噪声水平下的遥感图像的MPSNR和MSSIM值

Figure 7. Different methods for denoising the 206th band of HYDICE Urban dataset, (a) noise image (b) LRMR (c) LRTV (d) LRSID (e) ASSTV (f) LRTDGS (g) LMHTV (h) OUR

图7. 不同的方法对HYDICE Urban数据集的第206波段去噪结果，(a) 噪声图像(b) LRMR (c) LRTV (d) LRSID (e) ASSTV (f) LRTDGS (g) LMHTV (h) OUR

5. 结论

我们针对遥感图像被高水平混合噪声污染问题，提出了一种基于OGS-HL正则项的遥感图像混合噪声去除方法。首先通过对遥感图像的空间梯度值进行统计分析并使用重叠组稀疏性的超拉普拉斯先验来约束；其次对遥感图像退化问题中各种类型噪声分量进行分析，再根据条纹噪声独特的结构和特性，对每一张条纹噪声进行低秩约束，建立了关于遥感图像去除混合噪声恢复模型，由于该模型是非凸的，所以采用交替方向乘子法和MM算法有效地求解了所提出模型的最优解；最后，进行了仿真和真实的遥感数据实验，从视觉和定量的评估角度清晰地显示出所提出算法的优越性能，特别在混合噪声水平越大的情况下依然能够恢复出理想的数值结果。然而，所提出的模型还有改进的余地，例如：引入L₀范数替换L₁范数来约束光谱梯度值。

NOTES

¹https://www.ehu.eus/ccwintco/index.php?title = Hyperspectral_Remote_Sensing_Scenes#Indian_Pines

²https://github.com/liangjiandeng/PanCollection

³https://www.ehu.eus/ccwintco/index.php/Hyperspectral_Remote_Sensing_Scenes#Pavia_University_scene

⁴https://engineering.purdue.edu/~biehl/MultiSpec/hyperspectral.html

⁵http://www.escience.cn/people/feiyunZHU/DatasetGT.html