基于自适应阈值函数的四阶PDE去噪模型

Fourth-Order PDE Denoising Model Based on Adaptive Threshold Function

摘要: 针对偏微分方程在图像去噪过程中存在的阶梯效应，本文以YK模型为基础，引入差分曲率，将阈值和归一化差分曲率相结合形成自适应阈值函数，从而得到一个基于自适应阈值函数的新图像去噪模型。为验证模型的可行性，首先采用有限差分法进行数值计算，其次对Lena图像进行仿真实验，最后根据主观法和客观评价法对各模型去噪结果进行对比分析，结果表明文中提出的新模型去噪效果优于其他模型。

Abstract: Aiming at the step effect of partial differential equations in the process of image denoising, based on the YK model, this paper introduces differential curvature, and combines threshold and normalized differential curvature to form an adaptive threshold function, thereby obtaining a new image denoising model based on adaptive threshold function. In order to verify the feasibility of the model, the finite difference method is used for numerical calculation, and then the Lena image is simulated. Finally, the denoising results of each model are compared and analyzed according to the subjective method and the objective evaluation method, the results show that the denoising effect of the new model proposed in this paper is better than other models.

1. 引言

随着科学信息技术的飞速发展，数字图像处理技术已广泛应用到信息、工程、地球等科学领域。而数字图像在获取和传递的过程中易受到噪音污染，为获得质量更高的图像，国内外学者对图像去噪进行了诸多研究。目前图像去噪的方法有很多，主要分为变换域和空间域两大类去噪方法，其中空间域去噪方法包括均值滤波法、中值滤波法、非局部均值法及基于偏微分方程的方法等。随着数字图像处理技术的提高，基于数学理论的偏微分方程去噪方法已成为图像去噪领域的研究热点。

1990年Perona和Malik [1] 提出最经典的各向异性扩散方程，该模型的梯度下降流为：

$\frac{\partial u}{\partial t}=-\nabla ·\left(c\left(\left|\nabla u\right|\right)\nabla u\right)$， (1-1)

其中u表示图像的灰度值函数， $\nabla ·$ 表示散度， $\nabla$ 表示梯度算子， $c\left(\left|\nabla u\right|\right)$ 表示扩散系数。由于该模型在数学理论上来讲是一个病态模型，所有在去除噪音的同时也易产生阶梯效应。为克服PM模型出现的问题，You和Kaveh [2] 将梯度算子 $\nabla$ 改为拉普拉斯算子 $\Delta$，从而得到YK模型：

$\frac{\partial u}{\partial t}=-\Delta \left(c\left(\left|\Delta u\right|\right)\Delta u\right)$， (1-2)

虽然YK模型具有良好的去噪效果，但在去噪的过程中易出现过渡平滑和孤立的黑白斑点，后续学者纷纷对这两个模型进行研究，提出了许多改进的模型 [3] [4] [5]。

其中2014年Min等人 [6] 针对PM模型存在的问题，将梯度模扩展为灰度函数模，该模型的公式如下：

$\frac{\partial u}{\partial t}=-\Delta \left({\left(\frac{\left|u\right|}{k}\right)}^{-n}\Delta u\right)$， (1-3)

该模型利用灰度函数膜对图像进行去噪，根据不同的n对图像去噪得到不同的降噪效果。但该模型去噪过于模糊化，虽具有一定的去噪效果，但某些细节特征方面处理不佳。2019年闵莉花等人 [7] 利用图像灰度值和阈值函数提出的新去噪模型：

$\frac{\partial u}{\partial t}=-\Delta \left(\exp \left(-u/k\right)\Delta u\right)$， (1-4)

该模型虽然具有较好的去噪效果，但存在阈值k人工选择和图片不清晰的问题，原因可能是k值过大导致图像模糊化。

通过对以往模型的梳理，不难知道：阈值k是一个常数，需要人工给定，但在实验过程中很难确定合适的k。若k选择过大，易造成图像的过渡平滑；若k选择过小，导致图像去噪不彻底。因此本文引入差分曲率，对阈值k进行自适应改进。

本文其它部分构成如下：第二节介绍了什么是差分曲率，通过引入差分曲率，构建一种自适应阈值函数的PDE去噪模型；第三节先介绍去噪效果的评价方法，通过仿真实验得到各模型的去噪效果和评价指标，对比分析出各模型的去噪优劣；第四节对我们的成果进行了总结与展望。

2. 模型的改进

2.1. 差分曲率算子

差分曲率是由Chen等人 [8] (2010)首次提出的边缘检测算子，在图像去噪中主要用于区分边缘、平坦和独立噪声点区域。对应的表达式为：

$D=\left|\left|u_{\eta \eta }\right|-\left|u_{\xi \xi }\right|\right|$， (2-1)

其中 $u_{\eta \eta }$ 、 $u_{\xi \xi }$ 分别表示为灰度值函数图像u在梯度方向和垂直于梯度方向的二阶导数。设梯度方向和垂直于梯度方向的方向向量分别为：

$\eta =\left[u_x,u_y\right]/\sqrt{u_x^2+u_y^2}$， $\xi =\left[-u_y,u_x\right]/\sqrt{u_x^2+u_y^2}$ (2-2)

u对梯度方向的偏导为：

$u_{\eta }=\nabla u\cdot \eta =\left(u_x,u_y\right)\cdot \left[u_x,u_y\right]/\sqrt{u_x^2+u_y^2}=\left(u_x^2+u_y^2\right)/\left|\nabla u\right|$， (2-3)

其中 $\left|\nabla u\right|=\sqrt{u_x^2+u_y^2}$，从而得到u对梯度方向的二阶偏导：

$u_{\eta \eta }=\nabla u_{\eta }\cdot \eta =\left(\frac{u_x{u}_{xx}+u_y{u}_{xy}}{\sqrt{u_x^2+u_y^2}},\frac{u_y{u}_{yy}+u_x{u}_{xy}}{\sqrt{u_x^2+u_y^2}}\right)\cdot \frac{\left[u_x,u_y\right]}{\sqrt{u_x^2+u_y^2}}$， (2-4)

即：

$u_{\eta \eta }=\frac{u_{xx}{u}_x^2+2u_x{u}_y{u}_{xy}+u_{yy}{u}_y^2}{u_x^2+u_y^2}$， (2-5)

同理可得u对垂直于梯度方向的偏导和二阶导为：

$u_{\xi }=\left(u_x^2-u_y^2\right)/\left|\nabla u\right|$， $u_{\xi \xi }=\frac{u_{xx}{u}_y^2-2u_x{u}_y{u}_{xy}+u_{yy}{u}_x^2}{u_x^2+u_y^2}$， (2-6)

其中对x、y的偏导分别表示为：

$u_x=\frac{1}{2}\left[u\left(x+1,y\right)-u\left(x,y\right)\right]$， $u_y=\frac{1}{2}\left[u\left(x,y+1\right)-u\left(x,y\right)\right]$， (2-7)

对x、y的二阶偏导分别表示如下：

$u_{xx}=u\left(x+1,y\right)+u\left(x-1,y\right)-2u\left(x,y\right)$， (2-8)

$u_{yy}=u\left(x,y+1\right)+u\left(x,y-1\right)-2u\left(x,y\right)$， (2-9)

$u_{xy}=\frac{1}{4}\left[u\left(x+1,y+1\right)+u\left(x-1,y-1\right)-u\left(x-1,y+1\right)-u\left(x+1,y-1\right)\right]$. (2-10)

由公式(2-1)可以得到差分曲率算子D的以下特征：

1) 在边缘区域， $\left|u_{\eta \eta }\right|$ 较大， $\left|u_{\xi \xi }\right|$ 较小，所以D值较大；

2) 在平坦区域， $\left|u_{\eta \eta }\right|$ 和 $\left|u_{\xi \xi }\right|$ 都较小，所以D值较小；

3) 在独立噪声点区域， $\left|u_{\eta \eta }\right|$ 和 $\left|u_{\xi \xi }\right|$ 都较大，且几乎相等，所以D较小。

因此根据差分曲率D值大小，可以很好地区分边缘和平坦区域、独立噪声点。

2.2. 阈值函数的构建

针对前人提出的模型，可以看出阈值k需要人工给定，k值过大过小均会导致图像去噪效果不好。因此本小节从k入手，引入差分曲率，将阈值和归一化差分曲率相结合，得到了一个自适应阈值函数。

首先引入差分曲率D，对差分曲率进行归一化处理：

$d=\frac{D-\min D}{\max D-\min D}$， (2-11)

其次将阈值参数k与归一化差分曲率d相结合，形成一个自适应阈值函数：

$k^{\prime }=\frac{k_0}{d^{\beta }+1},\beta >0$， (2-12)

其中 $k_{\text{0}}$ 为初始阈值，现取 $k_{\text{0}}>k$，这里k为公式(1-4)中的值。通过对自适应阈值函数的分析可以得到如下特征：

1) 在平坦区域： $\left|u_{\eta \eta }\right|$ 较大， $\left|u_{\xi \xi }\right|$ 较大，从而 $D\to \text{0}$， $d\to 0$，有 $k^{\prime }>k$，因此增大了扩散强度；

2) 独立噪声区域： $\left|u_{\eta \eta }\right|$ 较小， $\left|u_{\xi \xi }\right|$ 较小，从而 $D\to \text{0}$， $d\to 0$，有 $k^{\prime }>k$，因此增大了扩散强度；

3) 在边缘区域： $\left|u_{\eta \eta }\right|$ 较大， $\left|u_{\xi \xi }\right|$ 较小，从而D较大， $d\to \text{1}$，当 $d>{\left(\frac{k_0-k}{k}\right)}^{\frac{1}{\beta }}$ 时，有 $k^{\prime }<k$，因此降低了扩散强度。

基于此，图像去噪能力得到提升的同时，边缘细节也得到了较好的保留。最后本文以YK模型为基础，引入差分曲率，将阈值和归一化差分曲率结合形成自适应阈值函数，提出一个基于自适应阈值函数的新图像去噪模型：

$\frac{\partial u}{\partial t}=-\Delta \left(\exp \left(-u/k^{\prime }\right)\Delta u\right)$. (2-13)

3. 仿真实验及结果分析

3.1. 去噪效果评价方法

本文的评价图像质量方法主要有两种：主观评价方法和客观评价方法。主观评价方法 [9] 主要依据降噪图像的视觉效果对降噪图像进行评价分析。客观评价方法 [10] 主要依据数学模型客观地对图像质量进行定量计算。本文采用的图像质量评价指标有信噪比(SNR)、峰值信噪比(PSNR)和均方根误差(RMSE)。其中PSNR反映的是图像处理前后数据变化的统计评价，其值越大，图像失真越少；SNR越大，去噪效果越好；RMSE越小，去噪能力越强。设图像大小为 $M\times N$， $u\left(x,y\right)$ 和 $U\left(x,y\right)$ 分别表示原图像和处理后的图像，则各评价指标的公式如下：

$\text{PSNR}=10\lg \left(\frac{{255}^2}{\frac{1}{MN}{\displaystyle \underset{i=1}{\overset{M}{\sum }}{\displaystyle \underset{j=1}{\overset{N}{\sum }}{\left[U\left(i,j\right)-u\left(i,j\right)\right]}^2}}}\right)$， (3-1)

和

$\text{SNR}=10\lg \left(\frac{{\displaystyle \underset{i=1}{\overset{M}{\sum }}{\displaystyle \underset{j=1}{\overset{N}{\sum }}{\left[u\left(i,j\right)-\bar{u}\left(i,j\right)\right]}^2}}}{{\displaystyle \underset{i=1}{\overset{M}{\sum }}{\displaystyle \underset{j=1}{\overset{N}{\sum }}{\left[U\left(i,j\right)-u\left(i,j\right)\right]}^2}}}\right)$， $\bar{u}\left(i,j\right)=\frac{1}{MN}{\displaystyle \underset{i=1}{\overset{M}{\sum }}{\displaystyle \underset{j=1}{\overset{N}{\sum }}u\left(i,j\right)}}$， (3-2)

$\text{RMSE}=\sqrt{\frac{{\displaystyle \underset{i=1}{\overset{M}{\sum }}{\displaystyle \underset{j=1}{\overset{N}{\sum }}{\left[U\left(i,j\right)-u\left(i,j\right)\right]}^2}}}{MN}}$. (3-3)

3.2. 数值求解

本小节使用有限差分法进行求解。u在梯度方向和垂直于梯度方向二阶导数表达如下：

${\left(u_{\eta \eta }\right)}_{i,j}^n=\frac{{\left(u_{xx}\right)}_{i,j}^n{\left({\left(u_x\right)}_{i,j}^n\right)}^2+2{\left(u_x\right)}_{i,j}^n{\left(u_y\right)}_{i,j}^n{\left(u_{xy}\right)}_{i,j}^n+{\left(u_{yy}\right)}_{i,j}^n{\left({\left(u_y\right)}_{i,j}^n\right)}^2}{{\left({\left(u_x\right)}_{i,j}^n\right)}^2+{\left({\left(u_y\right)}_{i,j}^n\right)}^2}$， (3-4)

${\left(u_{\xi \xi }\right)}_{i,j}^n=\frac{{\left(u_{xx}\right)}_{i,j}^n{\left({\left(u_y\right)}_{i,j}^n\right)}^2-2{\left(u_x\right)}_{i,j}^n{\left(u_y\right)}_{i,j}^n{\left(u_{xy}\right)}_{i,j}^n+{\left(u_{yy}\right)}_{i,j}^n{\left({\left(u_x\right)}_{i,j}^n\right)}^2}{{\left({\left(u_x\right)}_{i,j}^n\right)}^2+{\left({\left(u_y\right)}_{i,j}^n\right)}^2}$， (3-5)

从而得到差分曲率：

$D_{i,j}^n=\left|\left|{\left(u_{\eta \eta }\right)}_{i,j}^n\right|-\left|{\left(u_{\xi \xi }\right)}_{i,j}^n\right|\right|$， (3-6)

对差分曲率进行归一化处理：

$d_{i,j}^n=\frac{D_{i,j}^n-\min \left(D_{i,j}^n\right)}{\max \left(D_{i,j}^n\right)+\min \left(D_{i,j}^n\right)}$， (3-7)

故得到迭代后的自适应阈值函数：

${\left(k^{\prime }\right)}_{i,j}^n={\left(k_0\right)}_{i,j}^n/\left({\left(d_{i,j}^n\right)}^{\beta }+1\right),\beta >0$， (3-8)

现设

$g_{i,j}^n=\exp \left(-u_{i,j}^n/{\left(k^{\prime }\right)}_{i,j}^n\right)u_{i,j}^n$， (3-9)

采用中心差分格式计算上式的拉普拉斯算子：

$\Delta g_{i,j}^n=g_{i+1,j}^n+g_{i-1,j}^n+g_{i,j+1}^n+g_{i,j-1}^n-4g_{i,j}^n$， (3-10)

从而得到离散后的迭代方程：

$u_{i,j}^{n+1}-u_{i,j}^n=-\Delta t\cdot \Delta g_{i,j}^n$. (3-11)

3.3. 实验结果与分析

为验证新模型的可行性，本章使用Lena (图像大小：256 × 256)进行仿真实验。首先对Lena图像添加均值为0、方差为0.003的高斯噪声；其次固定时间步长 $\Delta t=0.01$，迭代次数 $n=100$ ；接着将提出的新模型与PM模型、YK模型、Min模型、式(1-4)模型分别进行降噪处理，最后根据主观视觉评价和客观评价指标对其去噪效果进行对比分析。经实验验证： $k=100$ 时，本章提出的新模型去噪效果达到最佳。模型处理后得到的图像视觉效果如图(图1)，各模型对应的客观评价指标对比分析如表(表1)。

Figure 1. Comparison of image denoising effects (mean 0, variance 0.003)

图1. 图像去噪效果对比(均值为0、方差为0.003)

从上图可以看出：PM模型对噪声图像处理效果不佳，且还出现了“阶梯效应”现象，这是因为PM模型的阈值参数具有一定的取值范围，当阈值不在对应的范围之内则去噪效果不好；YK模型虽然对该噪声图像具有一定的去噪效果，但某边缘区域比如头发、帽子和手臂等部分出现了很多的白色斑点，且从整体来看去噪后图像相对模糊；后面三个模型都具有较好的去噪效果，但Min模型和公式(1-4)模型对应的去噪后图像在某些细节方面比如帽子和右下角区域处理得不理想，而本文提出的新模型对图像的细节部分具有一定的保留性，比如红椭圆中白色的点亮度逐渐加强，且该模型去噪后的图像要比前面四个模型更加清晰特别是红框中呈现的图像区域部分。因此从视觉效果方面相比，本文模型在具有良好去噪性能的同时，也能较好的保护图像内部信息结构。

从以上数据来分析：PM模型的PSNR值和SNR值均低于噪声图像，RMSE值明显高于噪声图像，由此表明PM模型对含噪图像的去噪效果不好；其他模型的PSNR值和SNR值均高于噪声图像且RMSE值低于噪声图像，也说明了这四个模型对噪声图像具有一定的去噪效果。从PSNR值或者是SNR值来看，本文模型的PSNR值与公式(1-4)模型相比提高了1.27 dB，比YK模型升高了1.41 dB，本章模型的SNR值同Min模型和PM模型相比较，分别高出1.32 dB、14.92 dB；从RMSE值来看，新模型RMSE值低于PM模型41.53，和公式(1-4)模型相比降了1.43，同时也比Min模型和YK模型低了1.49、1.60。从而可以得出：本文模型与其他模型相比，其信噪比和峰值信噪比均得到了提高，均值方差有所下降，因此从客观指标数据来看本文提出的模型比其他模型更优。

Table 1. Comparison of objective indicators of image denoising effect

表1. 图像去噪效果客观指标对比

4. 结论

本文通过归纳总结以往的经典模型，引入差分曲率，将归一化的差分曲率与阈值函数相结合，提出一种基于自适应阈值函数的PDE模型。该模型的优越性在于阈值函数实现了自适应选择，去噪效果更好。通过对含噪图像进行仿真实验，各模型去噪结果与分析表明：本文提出的模型去噪性能优于其他模型。但我们的实验工作过于简单，只针对固定迭代次数的情况下进行实验，并没有选择最佳迭代次数；且该模型只进行了仿真实验并没有模型应用，因此接下来的工作需要对以上存在的问题进一步进行完善。