融合分数阶全变分和重叠组稀疏的遥感图像复原算法

Remote Sensing Image Restoration Algorithm Combining Fractional-Order Total Variation and Overlapping Group Sparse

摘要: 脉冲噪声的随机性和高对比度导致其在遥感图像中难以预测和定位，为了去除脉冲噪声，本文提出了一种融合分数阶全变分先验和重叠组稀疏先验的遥感图像复原算法。该模型采用l₀范数作为数据保真项以避免l₁范数的过度惩罚，利用重叠组稀疏先验来消除阶梯效应，同时分数阶全变分先验能够更有效地保留图像中的边缘和纹理信息。我们使用优化最小化算法和交替方向乘子法来进行求解，并与L0-OGSTV、HNHOTV、L0-TV三种算法进行对比，实验结果表明，本文所提出的算法在峰值信噪比和结构相似度上均优于其他几种算法。

Abstract: The randomness and high contrast of impulse noise cause it to be difficult to predict and localize in remote sensing images. In order to remove the impulse noise, this paper proposes a remote sensing image restoration algorithm that integrates fractional-order total variational prior and overlapping group sparse prior. The model adopts the l₀ norm as the data fidelity term to avoid the over-penalization of the l₁ norm, and utilizes the overlapping group sparse prior to eliminating the staircase effect, while the fractional-order total variation prior can retain the edge and texture information in the image more effectively. We use the majorization-minimization algorithm and the alternating direction multiplier method to solve the problem and compare it with the three algorithms, L0-OGSTV, HNHOTV, and L0-TV, and the experimental results show that the algorithm proposed in this paper outperforms the other algorithms in terms of the peak signal-to-noise ratio and the structural similarity.

1. 引言

遥感图像的成像过程具有高度复杂性，在曝光期间，由于成像系统与目标的相对运动，引起光轴偏移、像面抖动、目标物体的移动等因素，往往会导致所获取的图像受到模糊和噪声的影响。这些噪声和模糊严重影响了图像的质量和后续的分析应用。为了提升遥感图像的质量，去除模糊噪声是遥感图像领域的一项重要技术。图像退化的过程可以描述为原始图像和点扩散函数(PSF)的卷积，并加上噪声。图像去模糊可分为盲去卷积和非盲去卷积两类，在PSF已知的情况下来复原图像，称为非盲去卷积，反之则为盲去卷积[1]。全变分(TV)模型[2]是图像复原中流行的正则化方法之一，它能够有效地恢复出图像的边缘信息，但其考虑的是像素点的一阶梯度，往往会产生阶梯效应，降低了图像的对比度。一些研究者在TV模型上进行了改进，2010年，Bredies等人[3]在TV模型的基础上，提出了广义全变分(TGV)模型来去除噪声，其协调了图像的低阶导和高阶导，同时能够有效保护边缘信息和减轻阶梯效应。2010年，Pu等人[4]提出了另一种新的变分模型——分数阶全变分(FOTV)。2018年，Yang等人[5]提出了一种基于FOTV和自适应估计两个正则化参数的非盲图像去模糊模型，避免了选择正则化参数困难的问题，同时更好地恢复了图像的纹理细节。2020年，Kazemi Golbaghi等人[6]提出了基于整数阶和分数阶总变差的混合图像去噪方法。2022年，Zhu等人[7]提出了一种带有盒约束的新的变分模型，用于恢复受脉冲噪声污染的图像。所提出的模型由分数阶全变分正则化和$l_p$ ($0<p<1$ )保真项组成。此外，该模型具有保持锐利边缘和去除块效应的优点。在文献[8]中，R Parvaz利用分数阶全变分和桢变换，改进了用于处理脉冲噪声问题的图像恢复非凸模型。

近年来，重叠组稀疏模型受到广泛关注，相比于TV正则项，重叠组稀疏模型采用图像梯度的结构信息作为稀疏性的度量标准，从而改进了图像梯度的稀疏性，并且在恢复图像边缘的同时能够消除阶梯伪影。2021年Tarmizi Adam [9]提出了高阶非凸全变分与重叠群稀疏度相结合的算法来脉冲噪声去除。该模型对数据保真项和二阶全变差正则化项使用非凸$l_p$ 范数($0<p<1$ )，并结合了重叠组稀疏正则化器。在文献[10]中，Li等人提出了带有重叠组稀疏先验的$l_{2,p}$ 正则化全变差用于脉冲噪声图像修复，所提出的先验继承了全变差正则化器的优势，同时促进组稀疏性。2024年，Li等人[11]提出了一种具有重叠组稀疏性的广义非凸非光滑四向全变分图像复原算法，其有效缓解了与TV正则化相关的阶梯伪影,并通过利用图像像素的特定领域信息来提高恢复质量。在脉冲噪声条件下，$l_1$ 范数在稀疏信号和图像恢复方面具有很好的性能。然而，$l_1$ 范数往往会过度惩罚信号[12]。此外，$l_1$ 范数对脉冲噪声的离群值特征不具有鲁棒性[13] [14]。一些研究者使用$l_0$ 范数作为数据保真项，其实验结果表明，该方法是优于$l_1$ 范数的。

基于上述原因，本文提出一种融合分数阶全变分和重叠组稀疏的遥感图像复原算法，用于去除脉冲噪声。该模型采用$l_0$ 范数为数据保真项，FOTV和重叠组稀疏先验作为正则化项，避免了$l_1$ 范数正则化可能产生的过度惩罚问题，并能够利用图像梯度的结构信息，保留图像细节和边缘，同时减少传统全变差方法容易引入的阶梯伪影。

2. 预备知识

2.1. 分数阶全变分正则化

在图像处理中，通常使用Grünwald-Letnikov分数阶导数，分数阶定义为

${\nabla }^{\alpha }f\left(x\right)=\underset{h\to 0}{\lim }\frac{{\displaystyle \sum _{k\ge 0}{\left(-1\right)}^k{C}_k^{\alpha }f\left(x-kh\right)}}{h^{\alpha }},\alpha >0,$

其中权重$C_k^{\alpha }$ 为

$C_k^{\alpha }=\frac{\Gamma \left(\alpha +1\right)}{\Gamma \left(k+1\right)\Gamma \left(\alpha +1-k\right)},$ (1)

式(1)中$\Gamma (\cdot )$ 为伽马函数，对于指定的微分阶$\alpha$ 可取任意正实数，若$\alpha =1$ ，则为全变分模型，随着展开项数k的增加，$C_k^{\alpha }$ 趋近于0。当$h=0$ 时，分数阶梯度算子可以用前K项分数阶差分近似为以下形式：

${\nabla }^{\alpha }f\left(x\right)={\displaystyle \sum _{k=0}^{K-1}{\left(-1\right)}^k{C}_k^{\alpha }f}\left(x-k\right).$

给定一个图像，推广为二维离散化形式，其被离散为一个矩形网格$\left\{\left(x_i,y_j\right):1\le i\le m,1\le j\le n\right\}$ 。由此图像可以表示为欧几里得空间$R^{m\times n}$ 中的一个矩阵，记为$u_{i,j}=u\left(x_i,y_j\right)$ 。离散分数阶梯度定义为：

${\nabla }^{\alpha }u={\left[D_1^{\alpha }u,D_2^{\alpha }u\right]}^{\text{T}},$

其中沿x轴和y轴的分数阶导数$D_1^{\alpha }u,D_2^{\alpha }u\in R^{m\times n}$ 近似为

${\left(D_1^{\alpha }u\right)}_{i,j}={\displaystyle \sum _{k=0}^{K-1}{\left(-1\right)}^k}{C}_k^{\alpha }{u}_{i-k,j},$

${\left(D_2^{\alpha }u\right)}_{i,j}={\displaystyle \sum _{k=0}^{K-1}{\left(-1\right)}^k}{C}_k^{\alpha }{u}_{i,j-k},$

其中，K是在每个像素处近似分数阶导数的相邻像素的个数。将u的离散分数阶全变分定义为

${\Vert {\nabla }^{\alpha }u\Vert }_1={\displaystyle \sum _{i,j}\left(\left|{\left(D_1^{\alpha }u\right)}_{i,j}\right|+\left|{\left(D_2^{\alpha }u\right)}_{i,j}\right|\right)},$

其中${\Vert \cdot \Vert }_1$ 表示$l_1$ 范数，根据${\left({\nabla }^{\alpha }\right)}^{\ast }={\left(-1\right)}^{\alpha }di{v}^{\alpha }$ ，对于$p=\left(p^{\left(1\right)},p^{\left(2\right)}\right)\in R^{m\times n}\times R^{m\times n}$ ，离散分数阶散度 $di{v}^{\alpha }p\in R^{m\times n}$ 由[15]给出：

${\left(di{v}^{\alpha }p\right)}_{i,j}={\left(-1\right)}^{\alpha }{\displaystyle \sum _{k=0}^{K-1}{\left(-1\right)}^k{C}_k^{\alpha }}\left(p_{i+k,j}^{\left(1\right)}+p_{i,j+k}^{\left(2\right)}\right).$

2.2. 重叠组稀疏正则化器

对于二维图像$g\in {ℝ}^{n\times n}$ ，定义[16]一个$N\times N$ 的像素组$\tilde{g}\left(i,j\right)$ ：

${\tilde{g}}_{\left(i,j\right),N}=\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 {ℝ}^{N\times N},$

其中$m_1=\lfloor \frac{N-1}{2}\rfloor$ ，$m_2=\lfloor \frac{N}{2}\rfloor$ ，N表示组大小，$\lfloor \cdot \rfloor$ 表示下取整函数，即将任何实数转换为小于或等于它的最接近的整数。${\tilde{g}}_{\left(i,j\right),N}$ 可以看作是以g的$\left(i,j\right)$ 为中心的连续$N\times N$ 采样的正方形块。设$g_{\left(i,j\right),N}$ 是通过叠加${\tilde{g}}_{\left(i,j\right),N}$ 的N列而得到的向量，即$g_{\left(i,j\right),N}={\tilde{g}}_{\left(i,j\right),N}(:)$ 。因此，二维图像$g\in {ℝ}^{n\times n}$ 的重叠组稀疏正则化算子定义为

$\phi \left(g\right)={\displaystyle \sum _{i,j=1}{\Vert g_{\left(i,j\right),N}\Vert }_2}.$ (2)

3. 本文算法

3.1. 本文模型

图像的纹理和细节信息具有非局部相似性。根据图1，可以看出某一点的分数阶差分具有非局部特性，即该点的差分是邻域内若干个点的信息共同决定的。因此，分数阶差分在描述图像纹理信息时更加符合实际情况。另外，分数阶差分对高频信息的增强程度较低，能够有效防止将边缘信息误认为噪声点而去除，从而在一定程度上保留了图像的边缘和轮廓。而重叠群稀疏先验从全局角度出发，让图像的边缘轮廓得以保留。通过上述分析，本文设计了分数阶先验和重叠组稀疏先验相结合的遥感图像复原模型：

$\underset{0\le u\le 1}{\min }{\Vert O\odot \left(Hu-f\right)\Vert }_0+\lambda \phi \left(\nabla u\right)+{\Vert {\nabla }^{\alpha }u\Vert }_1$ (3)

其中H表示模糊核，$\lambda$ 为平衡参数，$\odot$ 表示哈达玛乘积，$O\in {\left\{0,1\right\}}^n$ 由使用者决定，当$O_i=0$ 时，表示位置i处的像素值为离群值，当$O_i=1$ 时，此时表示位置i处的像素值为潜在的离群值。当含噪图像$f_i$ 中的像素为$u_{\max }$ 或$u_{\min }$ 时，$O_i=0$ 。当$f_i$ 取其他值时，$O_i=1$ 。$\phi (\cdot )$ 是由式(2)定义的组稀疏正则化算子，$1\le \alpha \le 2$ 是分数阶差分阶数，需要说明的是：

(1) 当$0<\alpha <1$ 时，部分低频分量会非线性增加，而部分高频分量会非线性减少。这使得纹理和噪声的频率更加接近，增加了区分噪声和纹理的难度，导致二者信息同时被去除。

(2) 当$\alpha >1$ 时，部分低频分量会非线性减少，而部分高频分量会非线性增加。此时纹理和噪声的频率差异会被放大，使得区分纹理和噪声变得更容易。虽然$\alpha$ 的值越大，图像纹理信息的保留效果越好，但此时也会将图像中的一些边缘轮廓信息误判为噪声点，从而去除，导致恢复出的图片存在模糊现象，影响图像的视觉效果，因此本文设定$1\le \alpha \le 2$ 。

(a) 整数阶差分 (b) 分数阶差分

Figure 1. The difference between integer-order and fractional-order derivatives

图1. 整数阶差分和分数阶差分的差异

3.2. 模型求解

由于式(3)中涉及$l_0$ 范数最小化问题，根据文献[12]中的引理，我们首先将(3)式等价地表述为具有平衡约束的数学规划问题，如(4)式所示：

$\begin{array}{l}\underset{0\le u\le 1}{\min }\langle 1,1-s\rangle +\lambda \phi \left(\nabla u\right)+{\Vert {\nabla }^{\alpha }u\Vert }_1\\ \text{s}\text{.t}.s\odot \left|O\odot \left(Hu-f\right)\right|=0.\end{array}$ (4)

针对(4)式的求解，我们在ADMM框架下，引入辅助变量，得到以下等价的约束最小化问题：

$\underset{0\le u\le 1}{\min }\langle 1,1-s\rangle +\lambda \phi \left(\nabla u\right)+{\Vert {\nabla }^{\alpha }u\Vert }_1$

$\text{s}\text{.t}\text{.}x=\nabla u,y=Hu-f,x_1={\nabla }^{\alpha }u,s\odot \left|O\odot y\right|=s\odot O\odot \left|y\right|=0.$

其增广拉格朗日形式为

$\begin{array}{l}ℒ\left(u,s,x,y,x_1,{\theta }_1,{\theta }_2,{\theta }_3,{\theta }_4\right)\\ =\langle 1,1-s\rangle +\lambda \phi \left(x\right)+{\Vert x_1\Vert }_1+\langle \nabla u-x,{\theta }_1\rangle +\frac{{\beta }_1}{2}{\Vert \nabla u-x\Vert }_2^2+\langle Hu-f-y,{\theta }_2\rangle +\frac{{\beta }_2}{2}{\Vert Hu-f-y\Vert }_2^2\\ +\langle s\odot O\odot \left|y\right|,{\theta }_3\rangle +\frac{{\beta }_3}{2}{\Vert s\odot O\odot \left|y\right|\Vert }_2^2+\langle {\nabla }^{\alpha }u-x_1,{\theta }_4\rangle +\frac{{\beta }_4}{2}{\Vert {\nabla }^{\alpha }u-x_1\Vert }_2^2,\end{array}$

其中，${\beta }_1,{\beta }_2,{\beta }_3,{\beta }_4$ 为惩罚参数，${\theta }_1,{\theta }_2,{\theta }_3,{\theta }_4$ 为拉格朗日乘数。

下面，针对各个子问题进行求解：

(1) u子问题

$u^{k+1}=\underset{u}{\arg \min }\frac{{\beta }_1}{2}{\Vert \nabla u-x^k+\frac{{\theta }_1^k}{{\beta }_1}\Vert }_2^2+\frac{{\beta }_2}{2}{\Vert Hu-f-y^k+\frac{{\theta }_2^k}{{\beta }_2}\Vert }_2^2+\frac{{\beta }_4}{2}{\Vert {\nabla }^{\alpha }u-x_1^k+\frac{{\theta }_4^k}{{\beta }_4}\Vert }_2^2$ (5)

求解(5)式等价于求解下面关于u的线性方程组

$\left({\beta }_1{\nabla }^{\text{T}}\nabla +{\beta }_2{H}^{\text{T}}H+{\beta }_4{\left({\nabla }^{\alpha }\right)}^{*}{\nabla }^{\alpha }\right)u={\nabla }^{\text{T}}\left({\beta }_1{x}^k-{\theta }_1^k\right)+H^{\text{T}}\left({\beta }_{2}f+{\beta }_2{y}^k-{\theta }_2^k\right)+{\left({\nabla }^{\alpha }\right)}^{*}\left({\beta }_4{x}_1^k-{\theta }_4^k\right)$

在周期边界条件下，由于${\left({\nabla }^{\alpha }\right)}^{*}{\nabla }^{\alpha }$ ，$H^{\text{T}}H$ 和${\nabla }^{\text{T}}\nabla$ 有各自的循环变换矩阵，因此我们可以利用快速傅里叶变换和其逆变换对()式进行求解，得到u的闭合解为

$u^{k+1}={ℱ}^{-1}\left(\frac{ℱ\left({\nabla }^{\text{T}}\right)\odot ℱ\left({\beta }_1{x}^k-{\theta }_1^k\right)+ℱ\left(H^{\text{T}}\left({\beta }_{2}f+{\beta }_2{y}^k-{\theta }_2^k\right)\right)+ℱ\left({\left({\nabla }^{\alpha }\right)}^{*}\left({\beta }_4{x}_1^k-{\theta }_4^k\right)\right)}{{\beta }_2ℱ\left(H^{\text{T}}H\right)+{\beta }_1ℱ\left({\nabla }^{\text{T}}\right)\odot ℱ(\nabla )+{\beta }_4ℱ\left({\left({\nabla }^{\alpha }\right)}^{*}{\nabla }^{\alpha }\right)}\right)$ (6)

其中$ℱ$ 表示快速傅里叶变换，${ℱ}^{-1}$ 表示快速傅里叶逆变换。

(2) s子问题

$s^{k+1}=\underset{s}{\arg \min }\left({\theta }_3^k\odot O\odot \left|y^k\right|-1\right)s+\frac{1}{2}{\beta }_{3}O\odot y^k\odot y^k\odot s^2$ (7)

(7)式通过投影计算可以得到其解的形式为：

$s^{k+1}=\min \left(1,\max \left(0,-\frac{{\theta }_3^k\odot O\odot \left|y^k\right|-1}{{\beta }_{3}O\odot y^k\odot y^k}\right)\right)$ (8)

(3) x子问题

$x^{k+1}=\underset{x}{\arg \min }\frac{{\beta }_1}{2}{\Vert x-\left(\nabla u^{k+1}+\frac{{\theta }_1^k}{{\beta }_1}\right)\Vert }_2^2+\lambda \phi \left(x\right)$ (9)

(9)式是重叠组稀疏问题，由于此类问题难以直接进行求解，我们采用优化最小化(MM)算法对其进行求解，可以得到x子问题的解为

$x^{k+1}={\left[I+\frac{\lambda }{{\beta }_1}\Lambda {\left(x^k\right)}^{\text{T}}\left(x^k\right)\right]}^{-1}{x}_0,$ (10)

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

${\left[\Lambda \left(x\right)\right]}_{l,l}=\sqrt{{\displaystyle \sum _{i,j=-m_1}^{m_2}{\left[{\displaystyle \sum _{k_1,k_2=-m_1}^{m_2}{\left|x_{r-i+k_1,t-j+k_2}\right|}^2}\right]}^{-\frac{1}{2}}}}.$

(4) y子问题

$y^{k+1}=\underset{y}{\arg \min }\frac{{\beta }_2}{2}{\Vert y-\left(H{u}^{k+1}-f+\frac{{\theta }_2^k}{{\beta }_2}\right)\Vert }_2^2+\frac{{\beta }_3}{2}{\Vert s^{k+1}\odot O\odot \left|y\right|+\frac{{\theta }_3^k}{{\beta }_3}\Vert }_2^2$ (11)

将(11)式中的常数项舍弃，利用四维收缩算子计算(11)式，得到y子问题的解为

$y^{k+1}=\frac{H{u}^{k+1}-f+\frac{{\theta }_2^k}{{\beta }_2}}{\left|H{u}^{k+1}-f+\frac{{\theta }_2^k}{{\beta }_2}\right|}\odot \max \left(\frac{{\beta }_2\left|H{u}^{k+1}-f+\frac{{\theta }_2^k}{{\beta }_2}\right|-s^{k+1}\odot O\odot {\theta }_3^k}{{\beta }_2+{\beta }_3{s}^{k+1}\odot s^{k+1}\odot O},0\right).$ (12)

(5) $x_1$ 子问题

$x_1^{k+1}=\underset{x_1}{\arg \min }{\Vert x_1\Vert }_1+\frac{{\beta }_4}{2}{\Vert {\nabla }^{\alpha }{u}^{k+1}-x_1+\frac{{\theta }_4^k}{{\beta }_4}\Vert }_2^2$ (13)

利用一维收缩算子得到$x_1$ 子问题封闭形式的解为

$x_1^{k+1}=\text{shrink}\left({\nabla }^{\alpha }{u}^{k+1}+\frac{{\theta }_4^k}{{\beta }_4},\frac{1}{{\beta }_4}\right)$ (14)

其中

$\text{shrink}\left(v,\gamma \right)=\mathrm{sgn}\left(v\right)\odot \max \left\{\left|v\right|-\gamma ,0\right\}.$

(6) 更新拉格朗日乘子

$\begin{array}{l}{\theta }_1^{k+1}={\theta }_1^k+{\beta }_1\left(\nabla u^{k+1}-x^{k+1}\right)\\ {\theta }_2^{k+1}={\theta }_2^k+{\beta }_2\left(H{u}^{k+1}-f-y^{k+1}\right)\\ {\theta }_3^{k+1}={\theta }_3^k+{\beta }_3\left(s^{k+1}\odot O\odot \left|y^{k+1}\right|\right)\\ {\theta }_4^{k+1}={\theta }_4^k+{\beta }_4\left({\nabla }^{\alpha }{u}^{k+1}-x_1^{k+1}\right)\end{array}$ (15)

算法流程
步骤1：	Input $K,\lambda ,\alpha ,{\beta }_1,{\beta }_2,{\beta }_3,{\beta }_4$ ，初始化$u^0=f,{\theta }_1^0,{\theta }_2^0,{\theta }_3^0,{\theta }_4^0$ 。
步骤2：	给定迭代终止条件$tol=1\times {10}^{-4}$ ，给定迭代次数$k=1,2,\cdots ,\max$ ，迭代步骤如下： 根据(6)式更新u； 根据(8)式更新s； 根据(10)式更新x； 根据(12)式更新y； 根据(14)式更新$x_1$ ； 根据(15)式更新${\theta }_1,{\theta }_2,{\theta }_3,{\theta }_4$ ；
步骤3：	如果$\Vert u^{k+1}-u^k\Vert /\Vert u^k\Vert \le tol$ 迭代终止，否则继续迭代更新。 输出图像u。

4. 实验

在本节中，我们对本文所提出的模型进行脉冲噪声下遥感图像去模糊实验。本文实验是在处理器为Intel(R) Core(TM) i7-8750H CPU @ 2.20 GHz 2.21 GHz，运行内存为8.00 GB的64位Win10操作系统上进行的，所有仿真实验结果均是在MATLAB (2022a)环境中计算得出。为了验证本文所提出的算法对脉冲噪声遥感图像恢复的有效性，我们在开放的遥感卫星图像数据集中选取六幅大小为256 × 256不同场景的图片，分别为高尔夫球场、高速公路、海港、十字路口、住宅和天桥，如图2所示。

实验选用高斯模糊、均匀模糊和运动模糊，噪声均选用椒盐噪声模拟模糊噪声图像。每种模糊核选取如下尺寸，高斯模糊核宽度尺寸为7，滤波器的标准差为5。平均模糊核大小为7 × 7。运动模糊核的运动角度为8，运动的长度30，每种模糊情况分别添加密度为10%、20%和30%的椒盐噪声。每种模糊环境下的复原效果，我们分别选取两种场景进行展示，本文算法与对比实验算法的图像复原结果的对比图如图3~8所示。我们选择峰值信噪(PSNR)和结构相似度(SSIM)作为对比实验的客观评价指标，各算法的PSNR与SSIM的数值结果在表1、表3和表5中展示，运行时间结果在表2、表4和表6中展示。

Image1 Image2 Image3

Image4 Image5 Image6

Figure 2. Clear images

图2. 清晰图像

图3和图4展示了在10%椒盐噪声和高斯模糊下不同模型的去噪去模糊的视觉效果。在图3中，我们能够看到L0-OGSTV模型整体上的细节恢复效果不如本文算法，HNHOTV模型在右侧黑色区域边缘出现伪影，L0-TV模型在白色区域边缘处出现了伪影。在图4中，L0-TV和L0-OGSTV模型公路线在恢复效果上明显暗一些，HNHOTV模型在车头位置出现了伪影块。表1展示了在10%椒盐噪声和高斯模糊下不同模型的去噪去模糊的数值结果，表2展示了不同模型的运行时间。从表1中，可以看出本文算法的PSNR值和SSIM值均高于其他对比算法。

Figure 3. Restoration results of Gaussian blurred remote sensing image 1 with 10% impulse noise under different models

图3. 含有10%脉冲噪声的高斯模糊遥感图像1在不同模型下的复原结果

Figure 4. Restoration results of Gaussian blurred remote sensing image 2 with 10% impulse noise under different models

图4. 含有10%脉冲噪声的高斯模糊遥感图像2在不同模型下的复原结果

Table 1. Objective evaluation index of different models for recovering Gaussian blurred remote sensing images with 10% impulse noise

表1. 不同模型复原含有10%脉冲噪声的高斯模糊遥感图像的客观评价指标

image	L0-OGSTV		HNHOTV		L0-TV		Ours
	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM
1	31.8754	0.7925	31.8078	0.8033	30.9055	0.7464	31.9755	0.7998
2	28.3056	0.8357	28.6899	0.8426	26.6970	0.7792	28.8597	0.8455
3	30.7373	0.8450	30.5809	0.8455	28.86339	0.7866	30.8699	0.8500
4	30.4347	0.8637	30.5419	0.8623	28.5900	0.8134	30.6762	0.8670
5	31.9727	0.9494	31.7591	0.9406	29.4760	0.9253	32.0681	0.9480
6	31.4039	0.8600	31.4731	0.8593	27.7605	0.8017	31.5908	0.8644

Table 2. Runtimes for different models to recover Gaussian blurred remote sensing images with 10% impulse noise

表2. 不同模型复原含有10%脉冲噪声的高斯模糊遥感图像的运行时间

image	运行时间/s
	L0-OGSTV	HNHOTV	L0-TV	Ours
1	6.36	22.71	9.97	9.71
2	7.68	17.33	11.12	12.19
3	8.39	24.04	10.97	12.51
4	9.48	22.88	9.25	13.30
5	7.87	22.28	11.58	11.68
6	8.03	25.06	10.56	11.47

Figure 5. Restoration results of average blurred remote sensing image 3 with 20% impulse noise under different models

图5. 含有20%脉冲噪声的均匀模糊遥感图像3在不同模型下的复原结果

Figure 6. Restoration results of average blurred remote sensing image 4 with 20% impulse noise under different models

图6. 含有20%脉冲噪声的均匀模糊遥感图像4在不同模型下的复原结果

Table 3. Objective evaluation index of different models for recovering average blurred remote sensing images with 20% impulse noise

表 3. 不同模型复原含有20%脉冲噪声的均匀模糊遥感图像的客观评价指标

image	L0-OGSTV		HNHOTV		L0-TV		Ours
	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM
1	31.9709	0.7992	31.7406	0.8027	30.9898	0.7517	32.0352	0.8055
2	29.5968	0.8558	29.8082	0.8567	27.2664	0.7900	29.9403	0.8607
3	30.7822	0.8485	30.4126	0.8438	28.9657	0.7923	30.8464	0.8518
4	30.6301	0.8655	30.5869	0.8631	28.5704	0.8168	30.8226	0.8679
5	32.1153	0.9495	31.6951	0.9423	29.6930	0.9271	32.1809	0.9491
6	31.6260	0.8626	31.4814	0.8575	27.9299	0.8065	31.6925	0.8647

Table 4. Runtimes for different models to recover average blurred remote sensing images with 20% impulse noise

表4. 不同模型复原含有20%脉冲噪声的均匀模糊遥感图像的运行时间

image	运行时间/s
	L0-OGSTV	HNHOTV	L0-TV	Ours
1	7.24	22.73	10.94	9.88
2	10.24	22.67	13.74	14.61
3	9.25	23.00	10.28	13.51
4	10.17	22.65	10.72	16.03
5	9.72	22.17	10.08	11.70
6	8.31	22.72	10.51	11.69

图5和图6展示了在20%椒盐噪声和均匀模糊下不同模型的去噪去模糊的视觉效果。在图5中，我们能够看到L0-OGSTV模型减弱了图像的细节，HNHOTV模型加深了整个图像的纹路，L0-TV模型在恢复船的轮廓上没有达到很好的效果。在图6中，L0-TV在左面路面处出现伪影，HNHOTV在车处加深了道路，L0-OGSTV和HNHOTV模型在车的周围出现了伪影。表3展示了在20%椒盐噪声和均匀模糊下不同模型的去噪去模糊的数值结果，表4展示了不同模型的运行时间。从表3中，可以看出本文算法的PSNR值和SSIM值均高于其他对比算法。

图7和图8展示了在30%椒盐噪声和运动模糊下不同模型的去噪去模糊的视觉效果。在图7中，HNHOTV模型在右侧路面处出现伪影，观察住宅部分，能够看出本文模型在刻画纹理和细节区域时包含了更多的信息。在图8中，L0-OGSTV模型在恢复汽车边缘上没有达到很好的效果，HNHOTV模型在右上方马路边缘上有伪影，L0-TV模型在路灯光影的恢复上不够清晰。表5展示了在30%椒盐噪声和运动模糊下不同模型的去噪去模糊的数值结果，表6展示了不同模型的运行时间。从表5中，可以看出本文算法的PSNR值和SSIM值均高于其他对比算法。

Figure 7. Restoration results of motion blurred remote sensing image 5 with 30% impulse noise under different models

图7. 含有30%脉冲噪声的运动模糊遥感图像5在不同模型下的复原结果

Figure 8. Restoration results of motion blurred remote sensing image 6 with 30% impulse noise under different models

图 8. 含有30%脉冲噪声的运动模糊遥感图像6在不同模型下的复原结果

Table 5. Objective evaluation index of different models for recovering motion blurred remote sensing images with 30% impulse noise

表5. 不同模型复原含有30%脉冲噪声的运动模糊遥感图像的客观评价指标

image	L0-OGSTV		HNHOTV		L0-TV		Ours
	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM	PSNR	SSIM
1	34.0503	0.8852	32.0193	0.8366	33.2008	0.8601	34.0725	0.8833
2	33.9206	0.9349	30.9689	0.8876	32.3845	0.9178	34.0107	0.9367
3	32.4904	0.9016	29.8946	0.8576	31.2968	0.8780	32.5267	0.9019
4	34.4604	0.9403	31.6514	0.9056	32.5843	0.9262	34.5829	0.9408
5	35.0099	0.9666	30.2360	0.9161	33.4368	0.9634	35.1280	0.9693
6	34.8918	0.9290	31.0223	0.8682	29.6619	0.9014	34.9741	0.9306

Table 6. Runtimes for different models to recover motion blurred remote sensing images with 30% impulse noise

表6. 不同模型复原含有30%脉冲噪声的运动模糊遥感图像的运行时间

image	运行时间/s
	L0-OGSTV	HNHOTV	L0-TV	Ours
1	6.52	23.12	10.55	9.46
2	8.11	22.46	9.88	9.95
3	8.02	22.38	10.62	10.61
4	8.11	22.38	9.44	10.34
5	7.48	21.85	10.34	10.46
6	8.50	23.49	9.47	10.71

通过主观上视觉效果的对比，相比于其他模型，本文模型无论是在边缘轮廓的保留还是对复杂纹理细节的恢复上都具有更好的效果。表中的数值结果也客观地说明了本文模型在复原含有脉冲噪声的模糊遥感图像上可以达到更好的效果，更能接近原始清晰图像。

5. 结论

本文对开放的遥感卫星图像数据集中不同场景的图像进行图像复原研究，设计了分数阶全变分先验与重叠组稀疏先验相结合的遥感图像复原算法。该算法采用$l_0$ 范数为数据保真项，FOTV和重叠组稀疏先验作为正则化项，避免了$l_1$ 范数正则化可能产生的过度惩罚问题，并能够利用图像梯度的结构信息，保留图像细节和边缘，同时减少传统全变分方法容易引入的阶梯伪影，有效复原了脉冲噪声下的高斯模糊、均匀模糊以及运动模糊的遥感图像。