基于改进重叠组稀疏的遥感图像去噪算法

Remote Sensing Image Denoising Algorithm Based on Improved Overlap Group Sparsity

摘要: 脉冲噪声对遥感图像的质量有着显著的负面影响，它会破坏图像的连续性，降低图像的可视性和信息的准确性，从而影响遥感图像的应用效果。本文通过融合自适应中值滤波技术和组稀疏模型，设计了基于改进重叠组稀疏的模型，以有效地消除遥感图像的脉冲噪声，并消除梯度伪影现象。由于本文所提出的模型是非凸问题，我们利用最大–最小化(MM)方法和交替方向乘子法(ADMM)对模型进行求解。实验结果表明，本文提出的模型在峰值信噪比(PSNR)和结构相似度(SSIM)方面优于其他四种算法。

Abstract: Impulse noise has a significant negative effect on the quality of remote sensing image. It can destroy the continuity of image, reduce the visibility and sourcing circumstances of image, and affect the application effect of remote sensing image. By fusing adaptive median filter and group sparse model, a new model based on improved overlap group sparse model is designed to eliminate impulse noise and gradient artifacts in remote sensing images. Since the model presented in this paper is non-convex, we use the maximum-minimum (MM) method and the alternating direction multiplier (ADMM) method to solve the model. Experimental results show that the proposed model outperforms the other four algorithms in Peak signal-to-noise ratio (PSNR) and structural similarity (SSIM).

1. 引言

遥感图像的复原可以被视为一个线性逆问题，这是成像科学中的核心问题之一，对于众多图像处理应用而言至关重要。其中，去噪作为遥感图像复原的一个重要环节，旨在从受噪声影响的观测图像中估计并还原出原始清晰图像。

数学上，图像去噪的模型可以表示为：

$b=u+n,$ (1)

其中$b\in R^{n\times 1}$ 为降级后的图像，$u\in R^{n\times 1}$ 为期望的原始干净图像，$n\in R^{n\times 1}$ 为加性噪声。

遥感图像在采集和传输过程中会受到脉冲噪声的影响，脉冲噪声可以分为固定值脉冲噪声和随机值脉冲噪声(RVIN)两种类型，其中固定值脉冲噪声也被称为椒盐噪声(SPN)。在图像处理阶段，可以使用滤波算法、增强算法等技术来降低噪声并提高图像质量。以往，学者们已经提出了自适应中值滤波器(AMF) [1]来处理遥感图像中的脉冲噪声，然而AMF对细节丰富的图像可能产生不良影响。所以Erkan (2019) [2]，Zhang and Li (2014) [3]，Thanh (2020) [4]，Wang (2016) [5]，Khan (2018) [6]，Singh (2020) [7]，Mojica Vargas (2018) [8]，Enginoglu (2019) [9]，和Li (2014) [10]等人提出了一些基于平均滤波器的算法，虽然基于平均滤波器的算法简单易懂、高效，但同时也会导致边缘信息损失、细节信息丢失。最近，Mohd Rafi Lone等人提出的NNFM方法[11]利用距离较近的像素具有更高的相关性这一特点来去除脉冲噪声。虽然该方法通常能够较好地保留图像的细节和边缘信息，但该方法涉及的参数较多，如果参数设置不当可能会导致处理效果不佳或引入新的噪声。

但现在，学者们开始运用正则化技术构造优化模型，以更好地平衡去噪和细节保留之间的关系。最初，由Rudin等人将总变分模型(TV) [12]应用于去噪、去模糊、分割和超分辨率等问题。TV模型能够获得清晰的恢复结果并保留图像边缘，但复原的图像容易受到梯度伪影的影响。针对上述问题，Liu等人提出了重叠组稀疏(OGS)全变分图像恢复模型[13]。文献[14]中将超拉普拉斯先验与OGS结合使用，得到更稀疏的图像表示约束模型来有效地降低梯度伪影的影响。Selesnick和Chen [15]提出了用于一维TV去噪的OGS惩罚函数，数值实验表明该方法有效地减少了梯度伪影。综上所述，OGS正则化器已成功应用于混合图像恢复模型，特别是在消除梯度伪影方面[16]-[18]，取得了显著的成功。

在图像去噪中，$l_2$ 范数保真度被广泛用于恢复被加性高斯噪声损坏的图像[19]。然而，$l_2$ 范数对异常值敏感，当存在异常值时，图像恢复效果不理想。而基于贝叶斯统计原理，$l_1$ 范数保真度比$l_2$ 范数更适合恢复被脉冲噪声破坏的图像[20]。先前的研究[21]-[25]表明$l_1$ 范数在稀疏信号和脉冲噪声下的图像恢复中表现良好。在TV正则化器中使用$l_1$ 范数也有助于最小化凸目标函数[16]。尽管如此，$l_1$ 范数可能会产生过度惩罚的解[26]-[28]，在捕获脉冲噪声的异常值特征方面不够稳健。最近，文献[26]提出了一种基于$l_0$ $l_0$ 全变差($l_0$ -TV)脉冲噪声去除方法，该方法使用$l_0$ 范数作为数据保真度，并使用近端交替方向乘子法(PADMM)来恢复图像。结果表明，该方法的性能优于基于$l_1$ 范数的方法。文献[29]提出了一种重叠群稀疏全变分方法，该方法将数据保真度范数与重叠群稀疏全变分(OGSTV)相结合，用于恢复脉冲噪声下的模糊图像。

现有的去噪算法大多是基于模糊逻辑，这无疑增加了计算复杂度。所以，本文通过融合自适应中值滤波技术和组稀疏模型，设计了基于改进重叠组稀疏的模型，以有效地消除遥感图像的脉冲噪声。该模型首先利用改进的自适应中值滤波去除遥感图像中大部分脉冲噪声，并利用梯度域导向滤波保留图像边缘信息；利用自适应中值滤波计算噪声掩模矩阵，并利用该噪声掩模矩阵提高重叠组稀疏模型对脉冲噪声去除的能力，并消除图像中的梯度伪影。数值结果表明，我们的方法与其他四种先进的方法相比是非常有效和有竞争力的。

2. 预备知识

2.1. 改进的自适应中值滤波算法

在本节中，首先我们提供一些基本概念，以便在图像的滤波窗口中找到无噪声像素，然后给出了该方法及其实现的算法，如算法1。

我们设$u:={\left[u_{ij}\right]}_{m\times n}$ 为图像矩阵(IM)，其中$u_{ij}$ 是无符号整数，且$0\le u_{ij}\le 255$ 。如果$u_{ij}=0$ 或者$u_{ij}=255$ ，则$u_{ij}$ 称为u的噪声像素，否则，则为u的干净像素；如果存在某些i和j，使得$u_{ij}$ 是u的带噪声的像素，则称u为噪声图像矩阵(NIM)。如果设$u:={\left[u_{ij}\right]}_{m\times n}$ 为NIM，那么$C:={\left[c_{ij}\right]}_{m\times n}$ 称为u的二进制矩阵表示为：

$c_{ij}=\{\begin{array}{l}0，c_{ij}是u的噪声像素\\ 1，其他\end{array}，$ (2)

如果我们设$u:={\left[u_{ij}\right]}_{m\times n}$ 和$t\in \left\{1,2,\cdots ,\min \left\{m,n\right\}\right\}$ ，则${\left[{\overline{u}}_{rs}\right]}_{(m+2t)\times (n+2t)}$ 称为u的t对称矩阵，u的t对称矩阵用${\overline{u}}_{t-sym}$ 或者${\overline{u}}_t$ 表示，其定义如下：

$\left[\begin{array}{cccccccccc}{u}_{tt}& \dots & u_{t1}& u_{t1}& u_{t2}& \dots & u_{tn}& u_{tn}& \dots & u_{t(n-t+1)}\\ ⋮& \ddots & ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& \ddots & ⋮\\ u_{1t}& \cdots & u_{11}& u_{11}& u_{12}& \cdots & u_{1n}& u_{1n}& \cdots & u_{1(n-t+1)}\\ u_{1t}& \cdots & u_{11}& u_{11}& u_{12}& \cdots & u_{1n}& u_{1n}& \cdots & u_{1(n-t+1)}\\ u_{2t}& \cdots & u_{21}& u_{21}& u_{22}& \cdots & u_{2n}& u_{2n}& \cdots & u_{2(n-t+1)}\\ u_{3t}& \cdots & u_{31}& u_{31}& u_{32}& \cdots & u_{3n}& u_{3n}& \cdots & u_{3(n-t+1)}\\ ⋮& \ddots & ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& \ddots & ⋮\\ u_{mt}& \cdots & u_{m1}& u_{m1}& u_{m2}& \cdots & u_{mn}& u_{mn}& \cdots & u_{m(n-t+1)}\\ u_{mt}& \cdots & u_{m1}& u_{m1}& u_{m2}& \cdots & u_{mn}& u_{mn}& \cdots & u_{m(n-t+1)}\\ ⋮& \ddots & ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& \ddots & ⋮\\ u_{(m-t+1)t}& \cdots & u_{(m-t+1)1}& u_{(m-t+1)1}& u_{(m-t+1)2}& \cdots & u_{(m-t+1)n}& u_{(m-t+1)n}& \cdots & u_{(m-t+1)(n-t+1)}\end{array}\right]$ ，(3)

如果设$u:={\left[u_{ij}\right]}_{m\times n}$ 和$k\in \left\{1,2,\cdots ,t\right\}$ ，那么${\overline{u}}_t$ 的k近似矩阵$u_{ij}^k$ 定义为：

${\left[\begin{array}{ccc}{\overline{u}}_{(i+t-k)(j+t-k）}& \cdots & {\overline{u}}_{(i+t-k)(j+t+k）}\\ ⋮& {\overline{u}}_{(i+t)(j+t）}& ⋮\\ {\overline{u}}_{(i+t+k)(j+t-k）}& \cdots & {\overline{u}}_{(i+t+k)(j+t+k）}\end{array}\right]}_{(2k+1)\times (2k+1)},$ (4)

其中，由$u_{ij}^k$ 的所有正数项组成的非递减矩阵${\widehat{u}}_{ij}^k$ 称为$u_{ij}^k$ 的正则项矩阵(REM)。所有元素为零的矩阵称为零矩阵，记为[0]。最后我们计算中值$m=\text{median}\left({\widehat{u}}_{ij}^k\right)$，用于替换噪声像素。

算法1. 改进的自适应中值滤波算法

	步骤1：输入$u:={[u_{ij}]}_{m\times n}$ ，其中$\min \left\{m,n\right\}\ge 5.$
	步骤2：将变量$u$ 从uint8形式转换为double形式
	步骤3：对于$t$ 从5到1 获得关于$u$ 的二进制矩阵$C:={[c_{ij}]}_{m\times n}$ ，获得${\overline{u}}_t$ 和${\overline{C}}_t$ 对于所有$i$ 和$j$ 如果$c_{ij}=0$ 对于$k$ 从1到$t$ 如果$C_{ij}^k\ne 0$ 获得$u_{ij}^k$ 和${\widehat{u}}_{ij}^k$ $u_{ij}\leftarrow median({\widehat{u}}_{ij}^k)$ 结束循环 End If End For End If End For End For 步骤4：将$u$ 从double形式转换为uint8形式

本文采用改进的自适应脉冲噪声检测算法计算噪声掩膜矩阵$M$ ，该掩膜矩阵$M$ 的大小与图像的尺寸相同，其分量元素表示对应位置是否为噪声，噪点位置对应的分量设置为1，以便在后续处理中提高重叠组稀疏模型对脉冲噪声去除的能力，对于非噪声位置，其分量设置为0，以保留图像的信息。

2.2. 重叠组稀疏正则化器

为了能够更好地在去噪过程中消除TV先验所产生的块状伪影，文献[14]提出了一种重叠组稀疏正则化器。对于给定的向量$x\in R^n$ ，其$K$ 块组表示为：

$x_{i,k}=\left[x\left(i\right),\cdots ,x\left(i+K-1\right)\right]\in R^K,$ (5)

其中，$x_{i,K}$ 可以被视为$x$ 中从索引$i$ 开始的大小为$K$ 的连续样本块，常用的组稀疏正则化器[14] [29]定义为：

$\phi \left(x\right)={\displaystyle \sum _i{\left[{\displaystyle \sum _{k=0}^{K-1}{\left|x\left(i+k\right)\right|}^2}\right]}^{\frac{1}{2}}}.$ (6)

对于二维的情况，图像$u\in R^{n\times 1}$ 的$K\times K$ 块由[12]给出：

$\tilde{u}{\left(i,j\right)}_K=\left[\begin{array}{cccc}u\left(i-m_1,j-m_1\right)& u\left(i-m_1,j-m_1+1\right)& \cdots & u\left(i-m_1,j+m_2\right)\\ u\left(i-m_1+1,j-m_1\right)& u\left(i-m_1+1,j-m_1+1\right)& \cdots & u\left(i-m_1+1,j+m_2\right)\\ ⋮& ⋮& \ddots & ⋮\\ u\left(i+m_2,j-m_1\right)& u\left(i+m_2,j-m_1+1\right)& \cdots & u\left(i+m_2,j+m_2\right)\end{array}\right]\in R^{K\times K},$ (7)

式中$m_1=\left[\frac{K-1}{2}\right]$ ，$m_2=\left[\frac{K}{2}\right]$ ，其中$\left|x\right|$ 表示不大于$x$ 的最大正整数。$\tilde{u}{\left(i,j\right)}_K$ 的中心坐标是$\left(i,j\right)$ ，设$u{\left(i,j\right)}_K$ 为矩阵$\tilde{u}{\left(i,j\right)}_K$ 的$K$ 列叠加得到的向量，即$u{\left(i,j\right)}_K=\tilde{u}{\left(i,j\right)}_K(:)$ 。然后，二维重叠群稀疏正则化器可以表示为：

$\phi \left(u\right)={\displaystyle \sum _{i,j=1}{\Vert u{\left(i,j\right)}_K\Vert }_2}.$ (8)

在线性逆问题的背景下，OGS被用作正则化器以模拟信号在原始形式或其它领域中的分组或聚类行为，并已成功减少基于总变差的信号和图像恢复中的梯度伪影。

3. 本文模型及求解

3.1. 模型建立

本文提出一种基于改进重叠组稀疏的模型来去除脉冲噪声，模型如下：

$\underset{0\le u\le 1}{\min }{\Vert o\odot \left(u-b\right)\Vert }_0+\lambda \phi \left(\nabla u\right),$ (9)

其中，$\lambda >0$ 是控制OGS正则器的正则化参数(8)，u是经过改进的自适应中值滤波处理后并采用梯度域导向滤波进行增强的图像。为了便于计算和验证结果，我们只考虑位于$\left[0,1\right]$ 范围内的所有图像。

在本文模型中的数据保真度项使用$l_0$ 范数项，众所周知，$l_0$ 范数可以用来计算向量中非零元素的数量，从而实现稀疏性。而保真度项${\Vert o\odot \left(u-b\right)\Vert }_0$ 表示降级后的图像与需要的原始干净图像之间的误差。由于$l_0$ 范数计算向量中的非零元素，这使得$o\odot \left(u-b\right)$ 稀疏，并且非常适合脉冲噪声的稀疏性质。因此，本文的模型对脉冲噪声的异常值特征具有更强的鲁棒性。这就是在我们提出的模型中采用$l_0$ 范数作为稳健脉冲噪声恢复的数据保真度的主要原因。

那么问题(9)可以重新表述为：

$\begin{array}{l}\underset{0\le u,v\le 1}{\min }\langle 1,1-v\rangle +\lambda \phi \left(\nabla u\right),\\ s.t.v\odot \left|o\odot \left(u-b\right)\right|=0.\end{array}$ (10)

我们利用变量拆分，将问题重新表述为以下约束优化问题：

$\begin{array}{l}\underset{0\le u,v\le 1}{\min }\langle 1,1-v\rangle +\lambda \phi \left(x\right),\\ s.t.\nabla u=x,u-b=y,v\odot \left|o\odot y\right|=v\odot o\odot \left|y\right|=0.\end{array}$ (11)

3.2. 模型求解

为了解决(11)式问题，我们采用ADMM算法，那么(11)的增广拉格朗日函数为：

$\begin{array}{l}{ℒ}_{\mathcal{A}}\left(u,v,x,y,\pi ,{\pi }_z,{\pi }_0\right)=\langle 1,1-v\rangle +\lambda \phi \left(x\right)+\langle \nabla u-x,\pi \rangle +\frac{{\beta }_1}{2}{\Vert \nabla u-x\Vert }^2\\ \text{+}\langle u-b-y,{\pi }_z\rangle +\frac{{\beta }_2}{2}{\Vert u-b-y\Vert }^2\text{+}\langle v\odot o\odot \left|y\right|,{\pi }_0\rangle +\frac{{\beta }_3}{2}{\Vert v\odot o\odot \left|y\right|\Vert }^2,\end{array}$ (12)

其中，变量$\pi$ ，${\pi }_z$ 和${\pi }_0$ 是与(11)的约束相关的拉格朗日乘子。${\beta }_1$ ，${\beta }_2$ 和${\beta }_3>0$ 是二次罚分项${\Vert \nabla u-x\Vert }^2$ ，${\Vert u-b-y\Vert }^2$ 和${\Vert v\odot o\odot \left|y\right|\Vert }^2$ 所对应的惩罚参数。

基于ADMM方案，我们可以采用交替的方式解决以下问题。在此过程中，我们会深入探讨每个问题的解决方案，并在后文呈现完整的算法流程，即算法2，以清晰展示整个求解过程。

$\{\begin{array}{l}{u}^{k+1}=\arg {\min }_u\langle \nabla u-x^k,{\pi }_k\rangle +\frac{{\beta }_1}{2}{\Vert \nabla u-x^k\Vert }^2\text{ +}\langle u-b-y^k,{\pi }_z^k\rangle +\frac{{\beta }_2}{2}{\Vert u-b-y^k\Vert }^2,\\ v^{k+1}=\arg {\min }_v\langle v,1-v\rangle +\langle v\odot o\odot \left|y^k\right|,{\pi }_0^k\rangle \text{ +}\frac{{\beta }_3}{2}{\Vert v\odot o\odot \left|y^k\right|\Vert }^2,\\ x^{k+1}=\arg {\min }_x\lambda \phi \left(x\right)+\langle \nabla u^{k+1}-x,{\pi }^k\rangle +\frac{{\beta }_1}{2}{\Vert \nabla u^{k+1}-x\Vert }^2,\\ y^{k+1}=\arg {\min }_y\langle u^{k+1}-b-y,{\pi }_z^k\rangle +\frac{{\beta }_2}{2}{\Vert u^{k+1}-b-y\Vert }^2\\ \text{ +}\langle v^{k+1}\odot o\odot \left|y\right|,{\pi }_0^k\rangle +\frac{{\beta }_3}{2}{\Vert v^{k+1}\odot o\odot \left|y\right|\Vert }^2,\\ {\pi }^{k+1}={\pi }^k+{\beta }_1\left(\nabla u^{k+1}-x^{k+1}\right),\\ {\pi }_z^{k+1}={\pi }_z^k+{\beta }_2\left(u^{k+1}-b-y^{k+1}\right),\\ {\pi }_0^{k+1}={\pi }_0^k+{\beta }_3\left(v^{k+1}\odot o\odot \left|y^{k+1}\right|\right).\end{array}$ (13)

式(13)中的u的子问题等价于：

$u^{k+1}=\underset{u}{\arg \min }\frac{{\beta }_1}{2}{\Vert \nabla u-x^k+\frac{{\pi }^k}{{\beta }_1}\Vert }^2+\frac{{\beta }_2}{2}{\Vert u-b-y^k+\frac{{\pi }_z^k}{{\beta }_2}\Vert }^2,$ (14)

根据一阶最优性条件，这个问题要求我们解线性方程组：

$\left({\beta }_1{\nabla }^{\text{T}}\nabla +{\beta }_2\right)u={\nabla }^{\text{T}}\left({\beta }_1{x}^k-{\pi }^k\right)+\left({\beta }_{2}b+{\beta }_2{y}^k-{\pi }_z^k\right),$ (15)

在u的周期边界条件下，矩阵${\nabla }^{\text{T}}\nabla$ 具有块循环和循环块(BCCB)的结构，因此可以用二维离散傅里叶变换(FFT)对角化，即利用二维FFT和二维IFFT可以有效地求解u：

$u^{k+1}={ℱ}^{-1}\left(\frac{ℱ\left({\nabla }^{\text{T}}\right)\odot ℱ\left({\beta }_1{x}^k-{\pi }^k\right)+ℱ\left({\beta }_{2}b+{\beta }_2{y}^k-{\pi }_z^k\right)}{{\beta }_2ℱ+{\beta }_1ℱ\left({\nabla }^{\text{T}}\right)\odot ℱ(\nabla )}\right).$ (16)

对于求解$v$ 的子问题，可以得到下式：

$v^{k+1}=\arg \min \frac{1}{2}{\beta }_{3}o\odot y^k\odot y^k\odot v^2+\left({\pi }_0^k\odot o\odot \left|y^k\right|-1\right)v.$ (17)

因此，解$v^{k+1}$ 我们通过投影运算得到：

$v^{k+1}=\min \left(1,\max \left(0,-\frac{{\pi }_0^k\odot o\odot \left|y^k\right|-1}{{\beta }_{3}o\odot y^k\odot y^k}\right)\right).$ (18)

这个子问题是凸集上的投影，它能确保恢复图像的像素值在0到1之间。

对于求解式(13)中的x的子问题，其最小化函数可以表示为：

$\begin{array}{l}{x}^{k+1}=\underset{x}{\arg \min }\frac{{\beta }_1}{2}{\Vert x-\left(\nabla u^{k+1}+\frac{{\pi }^k}{{\beta }_1}\right)\Vert }^2+\lambda \phi \left(x\right),\\ =\underset{x}{\arg \min }\frac{1}{2}{\Vert x-\left(\nabla u^{k+1}+\frac{{\pi }^k}{{\beta }_1}\right)\Vert }^2+\frac{\lambda }{{\beta }_1}\phi \left(x\right).\end{array}$ (19)

不难看出，此问题我们可以使用MM方法[30]来解决，则：

$x^{k+1}=\left(I+\mu \Lambda {\left(x^k\right)}^{\text{T}}\Lambda {\left(x^k\right)}^{-1}\right)x_0,k=0,1,\cdots ,$ (20)

然后我们利用掩膜矩阵$M$ 进行更新，当迭代次数$k$ 小于3时，

$x^{k+1}=M\odot \left(I+\mu \Lambda {\left(x^k\right)}^{\text{T}}\Lambda {\left(x^k\right)}^{-1}\right)x_0+\nabla u^{k+1}\odot \left(1-M\right),k=0,1\cdots ,$ (21)

其中，$x_0=\nabla u^{k+1}+\frac{{\pi }^k}{{\beta }_1}$ 。

求解式(13)中的y的子问题等价于：

$y^{k+1}=\underset{y}{\arg \min }\frac{{\beta }_2}{2}{\Vert y-\left(u^{k+1}-b+\frac{{\pi }_z^k}{{\beta }_2}\right)\Vert }^2+\frac{{\beta }_3}{2}{\Vert v^{k+1}\odot o\odot \left|y\right|+\frac{{\pi }_0^k}{{\beta }_3}\Vert }^2,$ (22)

通过展开(22)并去掉常数项，我们可以将(22)重写为：

$y^{k+1}=\underset{y}{\text{argmin}}\frac{1}{2}{\Vert y-\frac{{\beta }_2\left(u^{k+1}-b+\frac{{\pi }_z^k}{{\beta }_2}\right)}{{\beta }_2+{\beta }_3{\left(v^{k+1}\odot o\right)}^2}\Vert }^2+\frac{v^{k+1}\odot o\odot {\pi }_0^k}{{\beta }_2+{\beta }_3{\left(v^{k+1}\odot o\right)}^2}\odot \left|y\right|.$ (23)

其最小值可以通过以下四维收缩算子计算得到：

$y^{k+1}=\frac{\frac{{\beta }_2\left(u^{k+1}-b+\frac{{\pi }_z^k}{{\beta }_2}\right)}{{\beta }_2+{\beta }_3{\left(v^{k+1}\odot o\right)}^2}}{\left|\frac{{\beta }_2\left(u^{k+1}-b+\frac{{\pi }_z^k}{{\beta }_2}\right)}{{\beta }_2+{\beta }_3{\left(v^{k+1}\odot o\right)}^2}\right|}\odot \max \left(\left|\frac{{\beta }_2\left(u^{k+1}-b+\frac{{\pi }_z^k}{{\beta }_2}\right)}{{\beta }_2+{\beta }_3{\left(v^{k+1}\odot o\right)}^2}\right|-\frac{v^{k+1}\odot o\odot {\pi }_0^k}{{\beta }_2+{\beta }_3{\left(v^{k+1}\odot o\right)}^2},0\right),$ (24)

化简(24)式可得到更新的$y^{k+1}$ ：

$y^{k+1}=\frac{u^{k+1}-b+\frac{{\pi }_z^k}{{\beta }_2}}{\left|u^{k+1}-b+\frac{{\pi }_z^k}{{\beta }_2}\right|}\odot \max \left(\frac{{\beta }_2\left|u^{k+1}-b+\frac{{\pi }_z^k}{{\beta }_2}\right|-v^{k+1}\odot o\odot {\pi }_0^k}{{\beta }_2+{\beta }_3{v}^{k+1}\odot v^{k+1}\odot o},0\right).$ (25)

最后，拉格朗日乘子$\pi$ 、${\pi }_z$ 、${\pi }_0$ 根据式(13)中的方程进行更新。

算法2. 本文算法

输入：正则化参数$\lambda$ ，组块大小$K$ ，惩罚参数${\beta }_1,{\beta }_2,{\beta }_3>0$ ，迭代次数$k$

初始化：初始图像$u^0=b$ ，$k=0$ ，拉格朗日乘子$\pi ,{\pi }_z,{\pi }_0$ 预先计算：矩阵 输出：恢复图像$u$ Fork = 0到迭代次数 Do

根据(16)式计算$u^{k+1}$ ，根据(18)式计算$v^{k+1}$ 根据(21)式计算$x^{k+1}$ ，根据(25)式计算$y^{k+1}$ 根据(13)式计算${\pi }^{k+1}$ ${\pi }_z^{k+1}$ ${\pi }_0^{\text{k+}1}$

If $\frac{\left|\right|u_{k+1}-u_k\left|\right|}{\left|\right|u_k\left|\right|}\le 1\times {10}^{-4}$ then Break $k=k+1$

End Do

4. 数值实验

本文实验部分全部使用Matlab2021a软件进行，电脑操作系统为Windows11，硬件平台为Intel(R) Core(TM) i9-12900H CPU @2.50GHz处理器、32GB RAM。我们将本文的方法与AMF [1]、NNFM [10]、l₀-TV [25]以及l₀-OGSTV [28]方法进行比较，并给出了数值结果，以说明所提出的算法对遥感图像恢复的有效性，所使用的测试图像均来自AID数据集如图1所示。

Figure 1. Test images of different scenarios

图1. 不同场景的测试图像.

4.1. 参数设置

为了获得可接受的图像恢复质量，我们首先设置这些参数，即组块大小K，内部迭代次数N和正则化参数$\lambda$ 。

为了找到最优的组块大小K，我们进行了脉冲噪声水平为30%的去噪实验，使用了“Iamge01”，“Iamge03”和“Iamge05”三个图像来获得K的最佳选择。当我们使其他参数保持不变时，令K的值发生变化，结果显示当K > 5时，PSNR和SSIM值会缓慢降低；当K = 5时，三幅图像的PSNR和SSIM值都最大。因此，在实验中我们将组块大小K设置为5。我们使用K = 5时的实验去除“Iamge01”和“Iamge05”两幅图像中 30%的脉冲噪声，以确定重叠组稀疏子问题的内部迭代次数N。当我们改变N时，且其他参数保持不变，增加内部迭代的次数会导致所提出算法的CPU时间增加，而恢复图像的PSNR和SSIM值几乎相同。因此，在实验中，我们设置N = 5来平衡CPU时间和恢复图像质量。最后一个需要适当调整的参数是正则化参数$\lambda$ ，它在噪声去除中起着重要的作用，该值取决于遥感图像本身以及噪声水平。通常，对于纯脉冲噪声的去除，其值在[0.10, 1.10]之间。

4.2. 实验结果及分析

在本文实验中，我们用不同程度的脉冲噪声破坏了九幅600 × 600的遥感图像。测试的具体噪声水平分别为10%、20%、30%和40%、50%、60%，这些值表示图像中被噪声损坏的像素的百分比。所提方法的所有参数均按前面的讨论得到。我们通过人工选择$l_0$ -TV和$l_0$ -OGSTV的正则化参数，以便得到最佳的PSNR和SSIM值，并且本文采用定量评价的方式，获得了五种算法的PSNR和SSIM值，见表1~9。

Figure 2. Comparison of the effects of different methods for removing pulse noise (ND = 20%) (Image 01)

图2. 不同方法去除脉冲噪声(ND = 20%)的效果对比图(Image 01)

Figure 3. Comparison of the effects of different methods for removing pulse noise (ND = 40%) (Image 07)

图3. 不同方法去除脉冲噪声(ND = 40%)的效果对比图(Image 07)

Figure 4. Comparison of the effects of different methods for removing pulse noise (ND = 60%) (Image 06)

图4. 不同方法去除脉冲噪声(ND = 60%)的效果对比图(Image 06)

Table 1. Evaluation index results of different algorithms at different density levels (10%~60%) (Image 01)

表1. 不同算法在不同密度水平(10%~60%)下的评价指标结果(Image 01)

ND (%)	AMF		NNFM		$l_0$ -TV		$l_0$ -OGSTV		OURS
	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM
10	38.4469	0.9795	37.0155	0.9795	39.1986	0.9773	40.2316	0.9807	40.5041	0.9812
20	36.3300	0.9692	34.7247	0.9635	37.1661	0.9650	38.3233	0.9697	38.4669	0.9701
30	34.2040	0.9533	33.2774	0.9479	35.5964	0.9514	36.8031	0.9581	36.9756	0.9587
40	32.3750	0.9313	32.0492	0.9317	34.2540	0.9361	35.5110	0.9452	35.6428	0.9459
50	30.6287	0.9025	31.1041	0.9145	32.8005	0.9162	34.1228	0.9291	34.2589	0.9299
60	28.1191	0.8524	30.0880	0.8933	31.2787	0.8894	32.6681	0.9087	32.7858	0.9094

Table 2. Evaluation index results of different algorithms at different density levels (10%~60%) (Image 02)

表2. 不同算法在不同密度水平(10%~60%)下的评价指标结果(Image 02)

ND (%)	AMF		NNFM		$l_0$ -TV		$l_0$ -OGSTV		OURS
	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM
10	34.7434	0.9609	40.7337	0.9889	35.8879	0.9675	36.8458	0.9724	36.8946	0.9723
20	33.1834	0.9489	31.3989	0.9350	33.9403	0.9507	34.8905	0.9577	34.9791	0.9580
30	31.3401	0.9297	30.3571	0.9270	32.5329	0.9349	33.6546	0.9449	34.8405	0.9449
40	29.6708	0.9020	29.1564	0.9039	31.4210	0.9179	32.4779	0.9300	32.4656	0.9301
50	28.1027	0.8663	28.1759	0.8791	30.2811	0.8970	31.3247	0.9117	31.2956	0.9118
60	25.9320	0.8085	27.2608	0.8519	29.0201	0.8699	30.1303	0.8904	30.0525	0.8898

Table 3. Evaluation index results of different algorithms at different density levels (10%~60%) (Image 03)

表3. 不同算法在不同密度水平(10%~60%)下的评价指标结果(Image 03)

ND (%)	AMF		NNFM		$l_0$ -TV		$l_0$ -OGSTV		OURS
	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM
10	38.1864	0.9831	36.7964	0.9858	37.4577	0.9805	40.0167	0.9863	40.2065	0.9865
20	35.4793	0.9740	33.9972	0.9721	34.8694	0.9676	37.5283	0.9768	37.7201	0.9770
30	33.1512	0.9605	32.3487	0.9587	33.2388	0.9558	35.8514	0.9679	36.0107	0.9683
40	31.0538	0.9416	30.9422	0.9441	32.1069	0.9426	34.4933	0.9575	34.6443	0.9580
50	29.1681	0.9169	29.6984	0.9277	30.8585	0.9271	33.1024	0.9454	33.2076	0.9460
60	26.8051	0.8689	28.5994	0.9086	29.3191	0.9019	31.4362	0.9261	31.5489	0.9272

通过观察图2可以看出，图像在脉冲噪声浓度不大的情况下，四种不同的算法都可以去除噪声，并且可以较好地保存图像细节。但是观察图2(b)可以看出，AMF处理后的图像会变得模糊，而相比于图2(b)，不难看出图2(c)经过NNFM处理后，可以更好地解决去噪过程中图像出现模糊的问题，但是图像飞机场跑道以及边界恢复效果不是很好，仍然有噪声存在。我们可以看出图2(d)和图2(e)的$l_0$ -TV方法和$l_0$ -OGSTV对于图像飞机场跑道部分恢复效果有明显提高，但是图像边界部分仍然存在伪影。而我们提出的方法基本上可以去除掉噪声和消除伪影，且能更好地恢复边界细节，并且在PSNR和SSIM值上也高于其它四种算法，见表1。

通过对商业街道图像的边缘区域进行视觉检测，如图3(b)和图3(c)所示，可以看出本文提出的方法在恢复体育场图像细节方面优于AMF和NNFM方法。值得注意的是，AMF和NNFM方法的处理会使图像部分区域仍然存在噪声，而通过观察图3(d)和图3(e)可以看出$l_0$ -TV和$l_0$ -OGSTV方法可以更好地去除噪声，但是在边缘区域仍然存在伪影。此外，图3(f)清楚地表明，本文提出的方法有效地去除了脉冲噪声，同时保持了图像边缘纹理，且在评价指标的数值上也有所体现，见表7。

Table 4. Evaluation index results of different algorithms at different density levels (10%~60%) (Image 04)

表4. 不同算法在不同密度水平(10%~60%)下的评价指标结果(Image 04)

ND (%)	AMF		NNFM		$l_0$ -TV		$l_0$ -OGSTV		OURS
	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM
10	41.3157	0.9884	39.9195	0.9897	42.6721	0.9901	44.3479	0.9926	44.5195	0.9928
20	38.8005	0.9825	37.1438	0.9799	40.4333	0.9839	42.0974	0.9877	42.2395	0.9879
30	36.4474	0.9732	35.5556	0.9709	38.5943	0.9774	40.3680	0.9825	40.4638	0.9827
40	34.3230	0.9599	34.3858	0.9621	36.9858	0.9691	38.9577	0.9767	39.0587	0.9770
50	32.3132	0.9412	33.4305	0.9520	35.5324	0.9586	37.3699	0.9684	37.4109	0.9686
60	29.0613	0.9023	32.4219	0.9399	33.8233	0.9429	35.6010	0.9568	35.6275	0.9570

Table 5. Evaluation index results of different algorithms at different density levels (10%~60%) (Image 05)

表5. 不同算法在不同密度水平(10%~60%)下的评价指标结果(Image 05)

ND (%)	AMF		NNFM		$l_0$ -TV		$l_0$ -OGSTV		OURS
	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM
10	33.5803	0.9745	31.2754	0.9764	36.7674	0.9905	40.2322	0.9951	40.6434	0.9953
20	30.9291	0.9618	29.0455	0.9750	32.8990	0.9781	36.1621	0.9880	36.3551	0.9882
30	28.6557	0.9420	27.4577	0.9363	30.6116	0.9619	33.5725	0.9776	33.6674	0.9779
40	26.6557	0.9156	26.2591	0.9156	28.7185	0.9414	31.2130	0.9624	31.1873	0.9623
50	25.7917	0.8795	25.2685	0.8927	27.1551	0.9145	29.3214	0.9415	29.1802	0.9405
60	23.0405	0.8208	24.2825	0.8658	25.5741	0.8765	27.5299	0.9131	27.1971	0.9099

Table 6. Evaluation index results of different algorithms at different density levels (10%~60%) (Image 06)

表6. 不同算法在不同密度水平(10%~60%)下的评价指标结果(Image 06)

ND (%)	AMF		NNFM		$l_0$ -TV		$l_0$ -OGSTV		OURS
	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM
10	30.2816	0.9228	26.8233	0.9035	36.1704	0.9853	37.6296	0.9891	37.7927	0.9893
20	28.7402	0.9033	25.7886	0.8699	32.5357	0.9637	34.0084	0.9732	34.1159	0.9734
30	26.7880	0.8691	24.8674	0.8348	30.0446	0.9340	31.4417	0.9495	31.5436	0.9499
40	25.3026	0.8250	24.0212	0.7977	28.2199	0.8943	29.4656	0.9149	29.5419	0.9154
50	23.7812	0.7682	23.2583	0.7568	26.5512	0.8426	27.7201	0.8685	27.7858	0.8692
60	22.2113	0.6922	22.4616	0.7102	25.0505	0.7761	26.1940	0.8114	26.2189	0.8119

通过比较图4所示的教堂图像中房顶区域的视觉质量，可以明显看出，本研究提出的方法在捕捉房屋边缘细节方面表现出色，优于其他方法。AMF和NNFM方法在去除噪声的过程中仍存在伪影，并且还留有部分噪声未去除。而本文方法在保留教堂边界的细节方面表现出优越的能力，并且梯度伪影抑制能力更为明显，且在数值结果上优于其它四种算法，见表6。

Table 7. Evaluation index results of different algorithms at different density levels (10%~60%) (Image 07)

表7. 不同算法在不同密度水平(10%~60%)下的评价指标结果(Image 07)

ND (%)	AMF		NNFM		$l_0$ -TV		$l_0$ -OGSTV		OURS
	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM
10	38.6161	0.9809	35.3577	0.9751	40.6558	0.9846	42.1039	0.9883	42.2142	0.9885
20	36.0397	0.9691	33.6355	0.9594	38.3553	0.9754	39.8262	0.9808	39.9405	0.9811
30	33.7753	0.9517	32.4033	0.9431	36.5690	0.9653	38.0976	0.9726	38.1587	0.9729
40	31.7791	0.9272	31.2754	0.9265	35.0286	0.9527	36.5913	0.9632	36.6395	0.9635
50	30.0133	0.8954	30.3078	0.9076	33.4074	0.9359	35.0240	0.9505	35.1016	0.9511
60	27.3885	0.8421	29.2087	0.8844	31.5860	0.9108	33.3808	0.9335	33.4590	0.9343

Table 8. Evaluation index results of different algorithms at different density levels (10%~60%) (Image 08)

表8. 不同算法在不同密度水平(10%~60%)下的评价指标结果(Image 08)

ND (%)	AMF		NNFM		$l_0$ -TV		$l_0$ -OGSTV		OURS
	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM
10	42.4882	0.9907	36.8547	0.9860	43.5801	0.9894	45.1743	0.9919	45.6352	0.9921
20	39.4442	0.9849	35.3957	0.9777	41.4466	0.9840	43.0504	0.9875	43.5852	0.9884
30	37.0059	0.9754	34.1668	0.9678	39.8384	0.9781	41.6387	0.9831	42.0620	0.9841
40	32.6700	0.9630	33.1096	0.9571	38.1855	0.9713	40.4359	0.9788	40.7819	0.9796
50	32.6720	0.9443	32.1961	0.9482	36.8176	0.9629	39.1907	0.9730	39.3935	0.9739
60	29.6331	0.9068	31.5237	0.9345	34.8508	0.9506	37.3188	0.9650	37.7388	0.9662

Table 9. Evaluation index results of different algorithms at different density levels (10%~60%) (Image 09)

表9. 不同算法在不同密度水平(10%~60%)下的评价指标结果(Image 09)

ND (%)	AMF		NNFM		$l_0$ -TV		$l_0$ -OGSTV		OURS
	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM
10	33.6372	0.9521	28.8733	0.9516	35.8022	0.9663	37.0898	0.9724	37.1080	0.9725
20	31.3748	0.9355	27.7039	0.9237	33.4953	0.9467	34.9315	0.9561	34.9842	0.9563
30	29.2126	0.9085	26.7600	0.8965	31.7315	0.9248	33.1468	0.9378	33.2321	0.9380
40	27.3701	0.8728	25.8908	0.8666	30.1865	0.8996	31.6362	0.9170	31.7218	0.9174
50	25.8212	0.8273	25.0977	0.8348	28.7745	0.8682	30.2326	0.8917	30.3473	0.8921
60	23.8153	0.7569	24.1709	0.7964	27.2301	0.8269	28.6359	0.8586	28.7093	0.8592

5. 结论

本文通过融合自适应中值滤波技术和组稀疏模型，设计了基于改进重叠组稀疏的模型，以有效地消除遥感图像的脉冲噪声。该模型首先利用改进的自适应中值滤波去除遥感图像中大部分脉冲噪声，并利用梯度域导向滤波保留图像边缘信息；利用自适应中值滤波计算噪声掩模矩阵，并利用该噪声掩模矩阵提高重叠组稀疏模型对脉冲噪声去除的能力，并消除图像中的梯度伪影。本文方法在有效地保留图像边缘的同时减少了梯度伪影，并且数值实验结果表明，我们提出的模型在PSNR和SSIM数值方面优于其它四种算法。

基金项目

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