河北工业大学，理学院，天津

收稿日期：2023年2月8日；录用日期：2023年3月4日；发布日期：2023年3月13日

摘要

由于图像复原具有不适定性，寻找一个有意义的图像先验仍然是图像处理上的重要挑战。在本篇论文中，提出了一个基于超拉普拉斯先验的交叠组稀疏模型用于去除脉冲噪声，且对于数据保真项采取Lp (0 < p < 1)范数。超拉普拉斯先验更能拟合自然图像的先验，交叠组稀疏是在图像去噪和图像去模糊中应用广泛。虽然是一项简单而广泛的研究，但该正则化在成像科学中仍有需求。为解决非凸非光滑最小化问题，我们采用交替乘子方向法(ADMM)作为主要算法框架。对于子问题分别使用FFT算法和MM算法。通过数值实验验证，所提出的模型在PSNR和SSIM值上是有效的。

关键词

超拉普拉斯先验，交叠组稀疏，脉冲噪声，ADMM

Removing Impulse Noise Method Based on Hyper Laplacian Prior on Overlapping Group Sparsity

Ping Bai

School of Science, Hebei University of Technology, Tianjin

Received: Feb. 8^th, 2023; accepted: Mar. 4^th, 2023; published: Mar. 13^th, 2023

ABSTRACT

Finding a meaningful image prior remains an important challenge in image processing due to the uncomfortable nature of image restoration. In this thesis, an overlapping group sparsity model based on a hyper-Laplacian prior is proposed for impulse noise removal, and the data fidelity term is taken to be the Lp (0 < p < 1) norm. The hyper-Laplacian prior is a better fit to the natural image prior, and overlapping group sparsity is widely used in image denoising and image deblurring. Although it is a simple and extensive study, the regularization is still in demand in imaging science. To solve the non-convex non-smooth minimization problem, the alternating multiplier direction method (ADMM) is used as the main algorithmic framework. For the subproblems, the FFT algorithm and the MM algorithm are used respectively. The proposed model is validated by numerical experiments for PSNR and SSIM values.

Keywords:Hyper-Laplacian Prior, Overlapping Group Sparsity, Impulse Noise, ADMM

Copyright © 2023 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

在景物成像过程中，由于目标的高速运动、散射、成像系统畸变和噪声干扰，致使最后形成的图像存在种种恶化，称之为“退化”。图像中常见的一种噪声类型是脉冲噪声。这类噪声通常由于成像过程中的错误切换、传输错误、电磁干扰和成像传感器的噪声而破坏图像。两种最常见的脉冲噪声是盐椒噪声和随机值噪声 [1] 。因此，恢复脉冲噪声污染的图像是图像处理的研究热点问题 [2] [3] 。

传统的去除脉冲噪声的图像恢复问题通常被建模为一个形式的线性逆问题

$\underset{u}{\min }\frac{1}{2}{\Vert Ku-f\Vert }_1+\lambda \psi \left(u\right),$ (1)

其中，第一项为数据保真项， $\psi \left(u\right)$ 为正则项， $\lambda$ 为正则项参数，用于平衡数据保真项和正则项。数据保真度项取决于噪声模型的噪声分布类型，而正则化项对图像的期望先验属性建模。

对于正则项，不同的选择有不一样的去噪模型。当 $\psi \left(u\right)={\Vert \nabla u\Vert }_1$ 该模型被称为TVL1 [4] 。然而，TV正则会产生阶梯效应 [5] 。为此，引入了一些修正的模型像高阶全变分(HOTVL1) [6] ，广义全变分正则(TGVL1) [7] 。

基于交叠组稀疏度(OGS) [8] 的模型用于自然图像去噪，并通过数值实验证明该模型比一些基于梯度的模型具有更好的鲁棒性。受OGS的启发，Selesnick和Chen [9] 考虑了组稀疏全变分正则模型，并将其应用于一维信号去噪问题。在文献 [10] 的基础上，Liu等人 [10] 将具有交叠组稀疏全变分的概念应用于脉冲噪声下的图像恢复，提出了一种脉冲噪声下图像去噪和去模糊的模型。然而，这个模型仍然存在轻微的阶梯效应。

L1正则项会产生过度惩罚的解，因此，使用L1范数的正则项或数据保真度项都受到这种影响。在这一点上，在图像重构中的正则化，由于非凸正则泛函的边缘保持能力，所以优先选择非凸正则泛函。而在数据保真项中，非凸泛函由于其对大误差和离群值的稳健性而被选。

在 [11] [12] [13] [14] [15] 中，涉及非凸泛函的数据保真度已被用于脉冲噪声下的图像恢复，这些泛函包括平滑限幅绝对偏差(SCAD)、指数型(ET)、Geman函数、LOG函数、Lp (0 < p < 1)范数和L0范数。然而，在这些工作中仍然使用TV正则化作为正则项，因此，阶梯效应仍然可能出现在恢复的图像中。此外，在 [16] 中也使用了非凸数据保真度并将其用于稀疏信号恢复和矩阵补齐。对于正则化的选择，非凸正则项已被证明可以产生清晰的边缘保持图像 [17] 。最近，文献 [18] 表明对于在超拉普拉斯先验上的交叠组稀疏正则具有鲁棒性。

基于超拉普拉斯先验可以很好地逼近自然图像梯度的重尾分布，交叠组稀疏可以通过引入额外的结构信息来弥补阶梯现象，我们引入了超拉普拉斯上的交叠组稀疏(OGS-HL)先验，用于脉冲噪声去噪，并且对于数据保真项采取Lp (0 < p < 1)范数。由于模型是非凸非光滑优化问题，ADMM方法非常适用。对于子问题采取FFT算法和MM算法进行求解，最后通过数值实验证明了图像去噪和去模糊的有效性。

2. 超拉普拉斯先验上的交叠组稀疏

Selesnic等人 [9] 定义了一个向量 $x\in R^n$ 的K点组： $x_{i,K}=\left[x\left(i\right),\cdots ,x\left(i+k-1\right)\right]\in R^K$ ， $x_{i,K}$ 可以看作是从索引i开始的K连续样本的一个块。对于一维情况，交叠组稀疏正则化项可以表示为：

$\psi \left(x\right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}{\Vert x_{i,K}\Vert }_2},$

其中，K表示组大小。

Liu等人 [10] 将交叠组稀疏正则从一维信号扩展到二维图像。一个图像矩阵 $f\in R^{n\times n}$ 的K × K点组表示为：

${\tilde{f}}_{i,j,K,K}=\left[\begin{array}{cccc}{f}_{i-a_1,j-a_1}& f_{i-a_1,j-a_1+1}& \cdots & f_{i-a_1,j+a_2}\\ f_{i-a_1+1,j-a_1}& f_{i-a_1+1,j-a_1+1}& \cdots & f_{i-a_1+1,j+a_2}\\ ⋮& ⋮& \ddots & ⋮\\ f_{i+a_1,j-a_1}& f_{i+a_1,j-a_1+1}& \cdots & f_{i+a_1,j+a_2}\end{array}\right]\in R^{K\times K},$

其中， $a_1=\lfloor \left(K-1\right)/2\rfloor$ ， $a_2=\lfloor K/2\rfloor$ ，且 $\lfloor \cdot \rfloor$ 表示小于或等于的最大整数。设 $f_{i,j,K,K}$ 是通过排列 ${\tilde{f}}_{i,j,K,K}$ 所有元素得到的向量。二维数组的交叠组稀疏正则化项为：

$\psi \left(f\right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}{\displaystyle \underset{j=1}{\overset{n}{\sum }}{\Vert f_{i,j,K,K}\Vert }_2}}.$

超拉普拉斯先验可以为图像梯度的重尾分布提供一个良性的近似，因此引起了人们的广泛兴趣。Jon等人 [18] 进一步研究了图像梯度的交叠组稀疏性，其稀疏性用Lq (0 < q < 1)伪范数测量。对于同样的K×K点组，定义超拉普拉斯先验交叠组稀疏正则化项为：

${\psi }_{\text{OH}}\left(f\right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}{\displaystyle \underset{j=1}{\overset{n}{\sum }}{\Vert {\left|f_{i,j,K,K}\right|}^q\Vert }_2}}={\displaystyle \underset{i=1}{\overset{n}{\sum }}{\displaystyle \underset{j=1}{\overset{n}{\sum }}\sqrt{{\displaystyle \underset{k_1=-a_1}{\overset{a_2}{\sum }}{\displaystyle \underset{k_2=-a_1}{\overset{a_2}{\sum }}{\left|f_{i+k_1,j+k_2}\right|}^{2q}}}}}}.$

其中，向量的 ${\left|.\right|}^q$ 是其元素是相应元素的绝对值的q次幂的向量。当K = 1时，上面式子是典型的超拉普拉斯先验；当q的取值为1时便是典型的交叠组稀疏正则模型，即 ${\psi }_{\text{OH}}\left(f\right)=\psi \left(f\right)$ 。

3. 提出的模型

在这一部分中，我们首先介绍了所提出的模型，然后在ADMM框架下对其进行求解。

脉冲噪声在现实生活中的已经存在，并且在生产生活中很容易产生脉冲噪声。我们提出了一个针对脉冲噪声的模型，该模型为：

$\underset{u}{\min }\lambda {\Vert Ku-f\Vert }_p^p+{\varphi }_{\text{OH}}\left(\nabla u\right)+I\left(u\right)$ (2)

其中， $0<p<\text{1}$ ， $\lambda >0$ ， $\nabla =\left(D_1,D_2\right)$ 。

为了解决所提出的模型，使用ADMM算法框架。通过引入辅助变量r，x₁，x₂，那么原始的复杂优化问题可以转化为以下的约束优化问题：

$\begin{array}{l}\underset{u}{\min }\lambda {\Vert r\Vert }_p^p+{\displaystyle \underset{i=1}{\overset{2}{\sum }}{\varphi }_{\text{OH}}}\left(v\right)+I\left(w\right)\\ \text{s}.\text{t}r=Ku-f,v=\nabla u,w=u.\end{array}$

那么，该约束优化问题增广拉格朗日方程可以表示为：

$\begin{array}{c}{L}_{\text{A}}\left(r,u,v,w,z\right)=\lambda {\Vert r\Vert }_p^p+{\varphi }_{\text{OH}}\left(v\right)+I\left(w\right)-z_v^{\text{T}}\left(v-\nabla u\right)+\frac{{\rho }_v}{2}{\Vert v-\nabla u\Vert }_2^2-z_w^{\text{T}}\left(w-u\right)\\ +\frac{{\rho }_w}{2}{\Vert w-u\Vert }_2^2-z_r^{\text{T}}\left(r-Ku+f\right)+\frac{{\rho }_r}{2}{\Vert r-Ku+f\Vert }_2^2\end{array}$

其中， $z_v$ ， $z_w$ 和 $z_r$ 是对偶变量(拉格朗日乘子)，变量 ${\rho }_1,{\rho }_2,{\rho }_3>0$ 为罚参数。

r子问题：r子问题即为求解如下极小化问题：

$r^{k+1}=\underset{r}{\arg \min }\frac{{\rho }_r}{2}{\Vert r-\left(K{u}^k-f+\frac{z_3^k}{{\rho }_3}\right)\Vert }_2^2+\lambda {\Vert r\Vert }_p^p$ (3)

令 $\bar{r}=K{u}^k-f+z_3^k/{\rho }_3$ ，(3)可以写成一个Lp (0 < p < 1)正则化问题：

$r^{k+1}=\underset{r}{\arg \min }\frac{1}{2}{\Vert r-\bar{r}\Vert }_2^2+\tau {\Vert r\Vert }_p^p,$

其中参数 $\tau$ 为正则化参数。求解非凸Lp正则化问题，可以采用IRL1算法 [19] 。定义 $g_i=\tau \cdot p/{\left(\left|r_i^k\right|+\epsilon \right)}^{1-p}$ 为权重， $r^{k+1}$ 的迭代为：

$r^{k+1}=\underset{r}{\arg \min }\frac{1}{2}{\Vert r-\bar{r}\Vert }_2^2+{\displaystyle \underset{i}{\sum }{g}_i}\left|r_i\right|,$

即为： $r^{k+1}=\text{sign}\left(\bar{r}\right)\odot \max \left(\left|\bar{r}\right|-g_i,0\right)$ .

u子问题：u子问题可以转化为以下极小化问题：

$\begin{array}{l}{u}^{k+1}=\underset{r}{\arg \min }\frac{{\rho }_r}{2}{\Vert r-Ku+f\Vert }_2^2-z_r^{\text{T}}\left(r-Ku+f\right)+\frac{{\rho }_v}{2}{\Vert v-\nabla u\Vert }_2^2\\ -z_v^{\text{T}}\left(v-\nabla u\right)+\frac{{\rho }_2}{2}{\Vert z-u\Vert }_2^2-z_2^{\text{T}}\left(z-u\right).\end{array}$

假定图像u符合周期边界条件，矩阵 ${\nabla }^{\text{T}}\nabla$ 和 $K^{\text{T}}K$ 是带循环块的块循环矩阵，可以通过2D快速傅立叶变换对角化。因此，u子问题的求解可以通过一个快速傅立叶变换(FFT)和一个逆FFT运算在傅立叶域中有效地求解。故此u子问题可以转化为：

$u^{k+1}={\mathcal{F}}^{-1}\left(rhs/lhs\right),$ (4)

其中，

$\begin{array}{l}lhs={\left|\mathcal{F}\left[K\right]\right|}^2+\frac{{\rho }_v}{{\rho }_r}{\left|\mathcal{F}[\nabla ]\right|}^2+\frac{{\rho }_z}{{\rho }_r},\\ rhs=F\left[K^{\text{T}}\left(r^{k+1}+f-\frac{z_3^k}{{\rho }_3}\right)+\frac{1}{{\rho }_1}{D}_1^{\text{T}}\left({\rho }_1{x}_1^k-z_1\right)+\frac{1}{{\rho }_1}{D}_2^{\text{T}}\left({\rho }_1{x}_2^k-z_1\right)+\frac{1}{{\rho }_r}\left({\rho }_z{z}^k-{\mu }_z\right)\right]\end{array}$

v子问题：v子问题可以转化为以下最小化问题：

$v^{k+1}=\underset{v}{\arg \min }\frac{{\rho }_v}{2}{\Vert v-\left(\nabla u^{k+1}+\frac{z_v^k}{{\rho }_v}\right)\Vert }_2^2+\frac{1}{{\rho }_v}{\varphi }_{\text{OH}}\left(v\right).$ (5)

为了解决这个子问题，采用了一种基于MM迭代方法，详见 [19] 。

w子问题：求解w子问题即为求解如下极小化问题：

$w^{k+1}=\underset{w}{\arg \min }\frac{{\rho }_w}{2}{\Vert w-\left(u^{k+1}+\frac{z_w^k}{{\rho }_w}\right)\Vert }_2^2+I\left(w\right)$ (6)

对于上述问题的求解可以使用投影法：

$w^{k+1}=P_{\Omega }\left[u^k+\frac{z_w^k}{{\rho }_w}\right]=\min \left(255,\max \left(u^{k+1}+\frac{z_w^k}{{\rho }_w},0\right)\right).$

拉格朗日乘子的更新：

$\begin{array}{l}{z}_1^{k+1}=z_1^k-{\rho }_v\left(v^{k+1}-\nabla u^{k+1}\right),\\ z_2^{k+1}=z_2^k-{\rho }_w\left(w^{k+1}-u^{k+1}\right),\\ z_3^{k+1}=z_3^k-{\rho }_3\left(r^{k+1}-K{u}^{k+1}+f\right).\end{array}$ (7)

因此，求解u的ADMM框架可以总结为算法1。

算法1.

4. 数值实验

针对带有脉冲噪声的图像问题进行了数值实验，并且与现有的其他算法进行了比较，通过数值实验证明了提出的方法的有效性。对于实验中的测试的图像，像素水平从256 × 256到640 × 640，具体的图像见图1。所有的实验结果用MATLAB 2019a在具备Intel(R) Core (TM) i5 2.7G处理器和4 GB的安装内存的Dell电脑上完成的。

本篇实验结果通过峰值信噪比(PSNR)和结构相似性(SSIM)两个指标来衡量恢复图像的质量。一般情况下，PSNR值越高意味着所恢复的图像质量越好，SSIM值小于等于1并且SSIM指数越接近1恢复的图像质量越好。

图1. 实验中的测试图像

PSNR被定义为：

$\text{PSNR}=10\log \frac{{\max }_{u,\hat{u}}^2}{1/N^2{\displaystyle {\sum }_{i,j=1}^N{\left[u_{i,j}-{\stackrel{⌢}{u}}_{i,j}\right]}^2}},$

其中， ${\max }_{u,\hat{u}}$ 是原始图像u和恢复图像 $\hat{u}$ 中的最大可能的像素值。其次，SSIM可以被定义为：

$\text{SSIM}=\frac{\left(2\bar{u}\overline{\stackrel{^}{u}}+C_1\right)\left(2{\delta }_{u,\hat{u}}+C_2\right)}{\left({\bar{u}}^2+{\overline{\stackrel{^}{u}}}^2+C_1\right)\left({\delta }_u^2+{\delta }_{\hat{u}}^2+C_2\right)},$

其中，

${\delta }_u^2=\frac{1}{N^2-1}{\displaystyle \underset{i,j=1}{\overset{N}{\sum }}{\left(u_{i,j}-\bar{u}\right)}^2},{\delta }_{\hat{u}}^2=\frac{1}{N^2-1}{\displaystyle \underset{i,j=1}{\overset{N}{\sum }}{\left({\hat{u}}_{i,j}-\overline{\stackrel{^}{u}}\right)}^2},{\delta }_{u,\hat{u}}=\frac{1}{N^2-1}{\displaystyle \underset{i,j=1}{\overset{N}{\sum }}\left(u_{i,j}-\bar{u}\right)\left({\hat{u}}_{i,j}-\overline{\stackrel{^}{u}}\right)}$

这里， $\bar{u}$ 表示u的平均值， $C_1={255}^2/{10}^4,C_2=9\cdot {255}^2/{10}^4$ 。一般来说，在实验中，我们算法的迭代终止准则为：

$\frac{\Vert u_{k+1}-u_k\Vert }{\Vert u_k\Vert }\le 1\times {10}^{-3},$

其中， $u_{k+1}$ 和 $u_k$ 分别表示在当前迭代和先前迭代的恢复图像。

4.1. 参数选择

首先，讨论了r子问题中的p参数。图2显示了该参数在去除高斯核为7 × 7噪声水平为0.3的椒盐噪声时PSNR和SSIM的影响。在图2中，观察到p在0.4之后趋于稳定。为了方便，取参数p = 0.5进行后续实验。

其次，讨论x_i (i = 1, 2)子问题中的q参数，图3显示了该参数在去除高斯核为7 × 7噪声水平为0.3的椒盐噪声时PSNR和SSIM的影响。对于参数q值的取值，在图3中可以看到当q接近1时，PSNR和SSIM逐渐增加，并且在q = 0.7，q = 0.8和q = 0.9时变得稳定。因此，在后续的实验中，简单的取q = 0.8值。

为了选择最优的内迭代次数N，我们选择测试图像中两个图像的“peppers”和“house”，图像“pepper”先被7 × 7高斯核模糊然后增加30%的椒盐噪声，图像“house”先被9 × 9平均核模糊然后增加40%的椒盐噪声。固定q = 0.8，p = 0.5，K = 3，λ = 120。在去模糊实验中，MM算法中组大小的值是变化的，其他参数保持不变。内部迭代次数N和PSNR、SSIM和CPU时间的值可见表1。从表1中，我们可以看到内迭代次数最优选择为5。

为了选择最优的交叠组稀疏正则项中的组大小K，我们选择测试图像中的“peppers”，“house”，“lighthouse”三个图像，这三个图像被先被7 × 7高斯核模糊然后增加30%的椒盐噪声。固定p = 0.5，q = 0.8，N = 5，λ = 120。MM迭代中的组大小值是变化的，其他参数保持不变。组大小K和图像的PSNR与SSIM 的关系详细见图4，且从图4中可以看出选择组大小为K = 3。

λ的值取决于图像和噪声水平。对于高斯核下的去模糊，λ的典型值在[75, 150]之间，而对于平均模糊情况，则在[35, 130]之间，在我们的实验中，取值为120。

最后是关于 ${\rho }_1,{\rho }_2,{\rho }_3$ 的选择，选择正确的 ${\rho }_1,{\rho }_2,{\rho }_3$ 值对算法的运行速度至关重要。对于非凸问题我们知道， ${\rho }_1,{\rho }_2,{\rho }_3$ 的值不能太小，如果太小，那么算法会不收敛，而如果太大，那么算法的运行速度非常地慢。我们根据经验设置了 ${\rho }_1=2.5,{\rho }_2=5,{\rho }_3=5$ 。

表1. 内迭代的选择

图2. p与PSNR、SSIM关系

图3. q与PSNR、SSIM关系

图4. 组大小与PSNR、SSIM关系

4.2. 图像去噪

为了说明该方法的优越性，我们将其与现有的模型进行了比较，即TVL1模型 [4] ，HOTVL1模型 [6] ，OGSL1模型 [10] ，HNHOTV-OGS模型 [16] 。在数值实验中，我们测试了实验中给出的6张图片。对于上述所提及的模型都采用ADMM算法求解，且上述的模型，取其相应论文中的各项参数，以达到最优的PSNR和SSIM值。图5、图6分别显示了TVL1模型，HOTVL1模型，OGSL1模型，HNHOTV-OGS模型的去噪效果和所提出的模型的去模糊结果。(a)~(g)是图像的整体效果，(a′)~(g′)是图像的局部效果。结果汇总在表2和表3中。从这些表中，可以观察到：首先，对于椒盐噪声(SP)和随机值噪声(RV)，所提出的方法几乎总是优于其他三种方法且PSNR和SSIM值比其他的四种方法高不少。因此，可以说，用该方法恢复的图像在结构上更接近于原始的干净图像。从表2和表3中的结果中可以观察到，当破坏为50%及以上(对于随机值的脉冲噪声为50%及以上)时，非凸数据保真度项(HNHOTV-OGS模型和所提出的模型)在抑制脉冲噪声方面的有效性。

图5. 显示了TVL1，HOTVL1，OGSL1，HNHOTV-OGS和proposed的去除增加5 × 5高斯核模糊和50%的椒盐噪声的去模糊结果

图6. 显示了TVL1，HOTVL1，OGSL1，HNHOTV-OGS和proposed在去除增加9 × 9平均核模糊和30%的随机值噪声的去模糊结果

表2. 现有的方法和所提出的方法对于增加5 × 5高斯核模糊和脉冲噪声的图像去模糊的PSNR值和SSIM值

表3. 现有的方法和所提出的方法对于增加9 × 9平均核模糊和脉冲噪声的图像去模糊的PSNR值和SSIM值

5. 结论

在这篇论文中，基于超拉普拉斯先验可以很好地逼近自然图像梯度的重尾分布，交叠组稀疏可以通过引入额外的结构信息来弥补阶梯现象，我们引入了超拉普拉斯上的交叠组稀疏先验，用于脉冲噪声去噪。数据保真项由于L1范数会产生过度惩罚的解，因此并不采L1范数而采取Lp (0 < p < 1)范数。由于该模型是非凸非光滑优化问题，ADMM方法非常适用。最后通过数值实验证明了图像去噪和去模糊的有效性。

文章引用

白 萍. 基于超拉普拉斯先验上的交叠组稀疏去除脉冲噪声方法
Removing Impulse Noise Method Based on Hyper Laplacian Prior on Overlapping Group Sparsity[J]. 应用数学进展, 2023, 12(03): 879-891. https://doi.org/10.12677/AAM.2023.123090

参考文献