CN111724318A - 基于混合高阶偏微分方程模型的图像去噪方法 - Google Patents

Abstract

本发明公开了一种基于混合高阶偏微分方程模型的图像去噪方法，包括如下步骤：对输入的含噪图像进行二维离散傅里叶变换，将图像从空域转换到频域，使待处理图像变换为离散化数字图像；采用有限差分法对离散化数字图像进行求解：计算每个离散化数字图像的具有二阶精度的中心差分值，并根据计算的中心差分值计算扩展图像的梯度模值和扩散张量，找到图像的边缘；根据混合高阶偏微分方程模型的边缘检测函数自适应检测扩展后图像的边缘；进行多次迭代直到结束，得到扩散图像的二维离散傅里叶变换图，即为去噪结果图像。所述方法既能够保护图像边缘信息又能够抑制区域内部的阶梯效应，图像处理效果好。

Description

基于混合高阶偏微分方程模型的图像去噪方法

技术领域

本发明涉及图像处理方法技术领域，尤其涉及一种基于混合高阶偏微分方程模型的图像去噪方法。

背景技术

数字图像是获取信息的主要来源，已经广泛应用于各类学科领域，但是图像采集信息时会受到环境背景光、CCD热噪声、读出噪声、A/D转换噪声、量化噪声的干扰，使得图像像质差、信噪比降低。因此，图像去噪是图像处理和计算机视觉中的首要问题。近年来，对噪声信号的恢复主要采用非线性手段，偏微分方程(PDE)方法是近年来的研究热点，具有很强的局域自适应性。基于Fourier变换的方法完全没有局域性，所以只能用于平稳信号的处理，不适于图像处理。基于小波变换的方法具有很好的时频域双重定域性，但是其采用的分离变量方法使得自适应能力受到限制，即便是脊小波和轮廓小波等自适应能力依旧不足。PDE是建立在连续模型之上的，仅依赖于某个像素点的邻域内，具有超强的局域自适应能力。PDE方法核心是利用图像的一个重要局部特征-梯度模值，将图像的滤波过程与图像边缘的检测过程结合起来，赋予其较线性滤波更优良的性能。

为了达到去噪并同时保护边缘的要求，Perona和Malik于1990年首先提出了扩散PDE模型，即PM模型。模型机理是传导系数依赖于图像梯度：在图像比较平坦的区域，传导系数能自动增大，使得区域中较小的不规则噪声被平滑；而在图像的边缘附近，传导系数能自动减小，使得边缘不受影响。虽然PM模型在图像处理中潜力巨大，但是数学研究表明其初值问题是病态的。Catte F等将PM模型进行了修改，引入了Gaussian函数，提出了正则化PM方程(亦称为CLMC模型)克服了PM方程的病态性质，并通过数学证明其方程是一个完全的适定问题。二阶PDE使用关于梯度模值的递减函数作为能量函数的积分。尽管这种PDE能够在去噪和保护边缘之间实现良好的折中，但往往会导致图像出现阶梯效应。

为解决这一问题，一些学者提出了高阶偏微分方程，尤其是四阶PDE法：2000年You和Kaveh提出的PDE(简称为YK)试图使成本函数最小化，从而避免了阶梯效应，其中成本函数是图像强度函数的拉普拉斯绝对值的递增函数。2002年Gilboa和Soche提出了自适应的向前向后(FABD)扩散模型，非线性扩散系数根据图像特征进行局部调整，可将扩散过程从前向模式转换为后向模式，去除局部噪声并将边缘锐化。数值实验证明：四阶模型能够避免阶梯效应。但是，图片的边缘可能会受到过度平滑的影响，使得边缘模糊。此外，在数学分析中，这些四阶非线性PDE很少开发。为了解决阶梯效应和边缘模糊问题，许多作者提出了一些混合高阶正则化模型。其中2016年，Dong和Chen提出了一个用于去噪的统一变分模型(DC模型)，在目标函数的正则化项中使用了不同分数阶导数的组合，在保存图像纹理和消除阶梯效应方面取得了很好的效果，但是并没有解决所建模型解的存在性问题。

发明内容

本发明所要解决的技术问题是如何提供一种既能够保护图像边缘信息又能够抑制区域内部的阶梯效应的图像去噪方法。

为解决上述技术问题，本发明所采取的技术方案是：一种基于混合高阶偏微分方程模型的图像去噪方法，其特征在于包括如下步骤：

对输入的含噪图像进行二维离散傅里叶变换，将图像从空域转换到频域，使待处理图像变换为离散化数字图像；

采用有限差分法对离散化数字图像进行求解：分别计算每个离散化数字图像在某点沿x方向和y方向的一阶差分

值以及二阶差分

值，并根据计算的各阶差分值计算扩展图像的一阶梯度模值

二阶梯度模值

和混合高阶偏微分方程模型所用的边缘检测函数

根据混合高阶偏微分方程模型的边缘检测函数自适应检测扩展后图像的边缘；

计算模型中的参数：

通过

进行迭代，其中γ是极小的正数保证有意义，

如果n＝I,计算直到结束得到扩散图像

的二维离散傅里叶变换图，即为去噪结果图像，迭代停止，否则令n＝I+1继续迭代；

d_1x、d_1y、d_2x、d_2y没有明确定义，就是为了方便设定的参数，具体是所建模型的两项内容

和

式中：含噪图像为u₀，并设置u_n＝u₀，α＝1,2，ω1,ω2∈{0,1,2,...,m-1}分别是对应于空间变量x和y的DFT频域变量，F是二维离散傅里叶变换，

表示的是对u_n进行二维离散傅里叶变换。

进一步的技术方案在于：假定原始离散图像u为m×m像素，含噪图像为u₀，并设置u_n＝u₀，用具有二阶精度的中心差分来近似计算对空间x和y的高阶偏导数，则

和

的近似值如下：

式中：α＝1,2，ω1,ω2∈{0,1,2,...,m-1}分别是对应于空间变量x和y的DFT频域变量，F是二维离散傅里叶变换，F^-1是F的逆变换，

表示的是对u_n进行二维离散傅里叶变换；

使用上述公式分别计算

的值。

进一步的技术方案在于，所述混合高阶偏微分方程模型的边缘检测函数的构造方法如下：

模型的特殊形式是整数阶导数，相关的泛函问题是：

式中，u为图像灰度函数，

表示u的梯度矢量模值即

用来检测边缘：当图像存在边缘时，梯度值较大，在平滑部分，梯度值较小，对平滑区域有抑制作用；u₀:Ω→R是采集到的含噪图像，

是具有Lipschitz边界的有界域；λ是大于零的参数；

D²u表示图像u的二阶导数，即

是具有Lipschitz边界的有界域，并且μ＞0,λ＞0是平衡目标函数中的参数；

h₁(x,y)是一个边界检测器，并由：

其中K₁和K₂是两个正常数；

考虑最小化问题：

式(7)对应的欧拉-拉格朗日方程如下：

是梯度向量，方向导数变化最快的方向，

是以梯度向量为速度场的散度，边界条件为：

使用最速下降法，导出相关的热流，得到最小化问题：

提出混合高阶偏微分方程：

进行参数选取后修改方程为：

式中，h(x,y)是边界检测器，定义为：

式中，ε为很小的一个正数，目的是确保边缘检测器小于1。

采用上述技术方案所产生的有益效果在于：所述方法使用混合高阶偏微分方程模型的边缘检测函数对图像进行处理，混合高阶偏微分方程模型将二阶项和四阶项进行结合，在区域内部四阶项起作用，平滑的同时抑制阶梯效应，在图像边缘处二阶项起作用，去噪的同时保护边缘信息不丢失。通过弱解存在性和数值计算验证了模型的有效性，并将其应用于图像去噪中，实验结果表明：通过所述方法去噪得到的图像在主观效果和客观指标上均优于传统去噪模型，较未改进的高阶去噪模型，在评价指标：SNR、PSNR和MSSIM上分别提高了1.89％、1.31％、和1.49％。

附图说明

下面结合附图和具体实施方式对本发明作进一步详细的说明。

图1a-图1c是本发明实施例中三种模型对lena图像进行去噪效果的结果图；

图2是本发明实施例中不同去噪方法的PSNR值随标准差的变化曲线图；

图3a-图3i为lena通过不同模型去噪后局部放大的结果图(图3a-图3i分别为原图、加噪图、双边滤波、小波全变分去噪、PM模型去噪、YK模型去噪、AFBD模型去噪、DC模型去噪以及本申请所述方法中使用的模型的去噪结果图)；

图4a-图4i为Cameraman通过不同模型去噪后局部放大的结果图(图4a-图4i分别为原图、加噪图、双边滤波、小波全变分去噪、PM模型去噪、YK模型去噪、AFBD模型去噪、DC模型去噪以及本申请所述方法中使用的模型的去噪结果图)；

图5a-图5i为Barbara的纹理图(图5a-图5i分别为原图、加噪图、双边滤波、小波全变分去噪、PM模型去噪、YK模型去噪、AFBD模型去噪、DC模型去噪以及本申请所述方法中使用的模型的去噪结果图)；

图6是本发明实施例所述方法的流程图。

具体实施方式

下面结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

在下面的描述中阐述了很多具体细节以便于充分理解本发明，但是本发明还可以采用其他不同于在此描述的其它方式来实施，本领域技术人员可以在不违背本发明内涵的情况下做类似推广，因此本发明不受下面公开的具体实施例的限制。

如图6所示，本发明实施例公开了一种基于混合高阶偏微分方程模型的图像去噪方法，包括如下步骤：

对输入的含噪图像进行二维离散傅里叶变换，将图像从空域转换到频域，使待处理图像变换为离散化数字图像；

采用有限差分法对离散化数字图像进行求解：计算每个离散化数字图像的具有二阶精度的中心差分值，并根据计算的中心差分值计算扩展图像的梯度模值和扩散张量，找到图像的边缘；

根据混合高阶偏微分方程模型的边缘检测函数自适应检测扩展后图像的边缘；

进行多次迭代直到结束，得到扩散图像的二维离散傅里叶变换图，即为去噪结果图像。

下面结合具体内容对上述方法进行详细说明：

模型的提出：

Perona-Malik类型的PDE模型与变分能量最小化问题有很强的联系，考虑图像还原的最小化问题：

式中，u为图像灰度函数，

表示u的梯度矢量模值即

用来检测边缘：当图像存在边缘时，梯度值较大，在平滑部分，梯度值较小。u₀:Ω→R是采集到的含噪图像，

是具有Lipschitz边界的有界域；λ是大于零的参数。式子前半部分为正则项，后半部分为保真项。与函数E(u)相关的形式梯度流由下式给出：

式中，

是梯度向量，方向导数变化最快的方向，

是以梯度向量为速度场的散度。扩散方程用于去噪的原理即：比较图像上的每一个点，若某一点的值小于周围点则增加它，否则就减小它，进而建立下列形式的自适应前后扩散方程(简称AFBD)：

这是从以下最小化问题得出的：

式中

是边缘检测器，由下式给出：

式中K是一个自然数；G_σ表示二维高斯核G_σ＝(2πσ)^-1exp(-(x²+y²)/2σ²)。*代表卷积。实验结果表明上述模型通过使用梯度相关性克服了热方程众所周知的边缘拖尾效应。但是，在有噪声的情况下得到的图像也显示出阶梯效应。

现有技术提出了一个变分模型(简称DC模型)，该模型使用了两个不同的分数阶导数的组合。在此，此模型的特殊形式是整数阶导数，相关的泛函问题是：

D²u表示图像u的二阶导数，即：

是具有Lipschitz边界的有界域，并且μ＞0,λ＞0是平衡目标函数中的参数。h₁(x,y)是一个边界检测器，并由：

其中K₁和K₂是两个正常数。上述模型的主要优点是在处理纹理方面的优势，以及消除阶梯效应。但是，并没有分析模型的解是否存在，是否是唯一存在的问题。

考虑最小化问题：

式(7)对应的欧拉-拉格朗日方程如下：

边界条件为：

使用最速下降法，能够导出相关的热流，得到最小化问题：

基于上述分析，本申请提出混合高阶偏微分方程：

进行参数选取后修改方程为：

式中，h(x,y)是边界检测器，定义为：

式中，ε为很小的一个正数，目的是确保边缘检测器小于1。为了简洁起见，省略了图像保真项，因为该值比较小在模型数学分析中影响可以忽略。

所提模型在图像区域内部，

和|D²u|都比较小，边界检测函数就近似等于1，这时，四阶项起主要作用，抑制阶梯效应。在图像的边界处，

和|D²u|都比较大，这时边界检测函数的值接近于0，模型的二阶项起主要作用，从而能保存图像的边缘。如图1a-图1c所示，比较了二阶模型、四阶模型和本申请所提模型对lena脸部区域图像的去噪效果：(a)为二阶模型去噪效果，图像的边缘得以保留但是出现了阶梯效应；(b)为四阶模型去噪效果，图像没有出现阶梯效应但是图像十分模糊；(c)为本申请模型去噪效果，结合了二阶和四阶的优点，既保留了图像的边缘信息又避免了阶梯效应。

弱解存在性分析：

为方程(12)添加初边值条件：

u(x,y,0)＝u₀(x,y)(x,y)∈Ω (14)

其中n＝(n₁,n₂)为边界

上的法向量。证明式(12)、(14)和(15)的弱解存在性，分为三部分进行：构造逼近解、估计逼近解和对逼近解求极限。

定义：若u∈L²(0,T；L²(Ω))∩BPV¹∩BPV²∩W^1,∞(0,T；[L²∩W^2,1(Ω)]^*) (16)

则称函数u是式(12)、(14)和(15)的弱解。并且满足：

其中，BPV^k代表Banach空间，Banach空间中的常微分方程如下：

取ε＝ε_k→0时，{u_k}为式(17)的解。取

则有：

||z_1,k||L^∞(Ω×(0,T)；R²)||≤1,||z_2,k||L^∞(Ω×(0,T)；R²)||≤1 (19)

固定测试函数w，求式(17)与2(u_k-w)的内积并对时间积分，得：

其中，

因为

所以有：

得知：

u_k→u在L^∞(0,T；L²(Ω))中弱*收敛，

在L2(0,T；(Hm(Ω))*)中弱收敛，

z_1,k→z₁,在L^∞(Ω×(0,T))中弱*收敛，

z_2,k→z₂,在L^∞(Ω×(0,T))中弱*收敛

故，初值条件式(14)成立。

在C([0,T]；L^∞(Ω))中强收敛，

D²(G_σ*u_k)＝G_σ*D²u_k→G_σ*D²u＝D²(G_σ*u)在C([0,T]；L^∞(Ω))中强收敛

有

对式(18)取测试函数v:

得：

左右两边同时取极限即为式(17)。故测试函数v和w满足条件。

数值计算和实验结果

本申请采用有限差分法求偏微分方程定解问题的解析解，为了得到离散化数字图像，将时间[0，T]等间隔细分为t_n＝nΔt,Δt＝T/L,n＝0,...,L，其中L是一个正整数，假定原始离散图像u为m×m像素。

用具有二阶精度的中心差分来近似对空间x的高阶偏导数，则

和

的近似值(参见公式(13))如下：

式中：α＝1,2，ω1,ω2∈{0,1,2,...,m-1}分别是对应于空间变量x和y的DFT频域变量。F是二维离散傅里叶变换，F^-1是F的逆变换，

表示的是对u进行二维离散傅里叶变换；

算子共轭由如下公式进行计算：

式中conj(·)是复共轭。

当图像尺寸m为奇数时，

是实数；当m为偶数时，

包含复数项。但是，当m够大时，复数项很小可以忽略，所以，当m为偶数时，可以通过扩展初始图像的方法消除复数项：

其中，

为u₀的扩展图像，大小为(m+1)×(m+1)，是奇数。

鉴于式(12)中

和|D²u|作为分子出现，为使式子有意义，使用

和

进行数值计算，其中γ为极小的正数，本申请取其值为10^-6。以下是本申请所提的基于混合高阶偏微分方程模型应用于图像去噪中的具体流程：

1)对输入的含噪图像进行二维离散傅里叶变换，将图像从空域转换到频域；

2)待处理图像已经是离散化数字图像，符合有限差分法的条件，所以采用有限差分法求解。分别计算每个离散化数字图像在某点沿x方向和y方向的一阶差分

值以及二阶差分

值，并根据计算的各阶差分值计算扩展图像的一阶梯度模值

二阶梯度模值

和混合高阶偏微分方程模型所用的边缘检测函数为

3)根据边缘检测函数

自适应检测图像边缘，在边缘附近

都足够大，边缘检测函数接近于零，扩散几乎被停止。

4)进行多次迭代直到结束，得到扩散图像的二维离散傅里叶变换图，即为去噪结果图像。

实验结果：

为了证明本申请所提混合高阶模型去噪和保持边缘的效果优势，用传统双边滤波和小波变换方法与偏微分方程方法：PM模型、YK模型、AFBD模型(式(3))、DC模型(式(6))、本申请模型(式(12))进行了去噪性能对比，使用三个评价标准对去噪后的图像进行了判定。其中信噪比(SNR)和峰值信噪比(PSNR)用来衡量图像去噪的客观指标，平均结构相似度(MSSIM)用来衡量图像去噪的视觉保真度。实验参数选取如下：Δt＝1/16,K₁＝0.02,K₂＝0.02,μ＝1/5,ε＝10^-6，高斯核为5×5。各模型中的参数选取与文献保持一致。

图2为Barbara图像加入不同噪声后，采用不同模型进行去噪得到的PSNR值变化曲线。

图3a-图3i、图4a-图4i以及图5a-图5i呈现了加入σ＝25的噪声后，不同模型的去噪结果。图3a-图3i为lena通过不同模型去噪后局部放大的结果，以便更好观察细节信息。从图看出，双边滤波结果还是存在大量的噪声；小波全变分去噪过渡平滑丢失了边缘信息；PM模型能够很好的保留图像的边缘信息，但是图像会出现阶梯效应和斑点；YK模型消除了分块现象，图像更加自然，但是依旧存在孤立的斑点且理论上不适定；AFBD模型也能够很好的保留边缘，减小阶梯效应，但是还会出现轻微的阶梯效应；DC模型和本申请所提模型能够兼顾两个问题，Lena脸部没有出现阶梯效应和斑点，并且很好的保留了边缘信息。图4a-图4i为Cameraman通过不同模型去噪后局部放大的结果，PM模型和AFBD模型依然会使图像产生阶梯效应，而DC模型和本申请模型没有产生接替效应并将边界纹理很好的保留下来。图5a-图5i是Barbara的纹理图，同前者一致，本申请模型能够在去噪的同时将纹理信息保留。通过图像的视觉对比，可以看出本申请所提模型在去噪和保存纹理两方面的优势。

为更清晰地展现去噪效果对比，列出不同模型去噪的SNR、PSNR、MSSIM值结果，如表1所示。从表中数据可知本申请模型具有最高的SNR、PSNR和MSSIM值，较其他模型有很大的优越性。

综上，本申请提出了一种混合高阶偏微分方程模型，该模型将二阶项和四阶项进行结合，在区域内部四阶项起作用，平滑的同时抑制阶梯效应，在图像边缘处二阶项起作用，去噪的同时保护边缘信息不丢失。通过弱解存在性和数值计算验证了模型的有效性，并将其应用于图像去噪中，实验结果表明：本申请所述去噪模型得到的图像在主观效果和客观指标上均优于传统去噪模型，较未改进的高阶去噪模型，在评价指标：SNR、PSNR和MSSIM上分别提高了1.89％、1.31％、和1.49％。

Claims (3)

1.一种基于混合高阶偏微分方程模型的图像去噪方法，其特征在于包括如下步骤：

对输入的含噪图像u₀进行二维离散傅里叶变换，将图像从空域转换到频域，使待处理图像变换为离散化数字图像；

采用有限差分法对离散化数字图像进行求解：分别计算每个离散化数字图像在某点沿x方向和y方向的一阶差分

值以及二阶差分

值，并根据计算的各阶差分值计算扩展图像的一阶梯度模值

二阶梯度模值

和混合高阶偏微分方程模型所用的边缘检测函数

根据混合高阶偏微分方程模型的边缘检测函数自适应检测扩展后图像的边缘；

计算混合高阶偏微分方程模型中的参数：

通过

进行迭代，其中γ是极小的正数保证有意义，

如果n＝I,计算直到结束得到扩散图像

的二维离散傅里叶变换图，即为去噪结果图像，迭代停止，否则令n＝I+1继续迭代；

将时间[0，T]等间隔细分为t_n＝nΔt,Δt＝T/L,n＝0,...,L，其中L是一个正整数，d_1x、d_1y、d_2x、d_2y无明确定义，具体是所建混合高阶偏微分方程模型中的两项内容

和

含噪图像为u₀，并设置u_n＝u₀，α＝1,2，ω1,ω2∈{0,1,2,...,m-1}分别是对应于空间变量x和y的DFT频域变量，F是二维离散傅里叶变换，

表示的是对u_n进行二维离散傅里叶变换。

2.如权利要求1所述的基于混合高阶偏微分方程模型的图像去噪方法，其特征在于：

假定原始离散图像u为m×m像素，含噪图像为u₀，并设置u_n＝u₀，用具有二阶精度的中心差分来近似计算对空间x和y的高阶偏导数，则

和

的近似值如下：

式中：α＝1,2，ω1,ω2∈{0,1,2,...,m-1}分别是对应于空间变量x和y的DFT频域变量，F是二维离散傅里叶变换，F^-1是F的逆变换，

表示的是对u_n进行二维离散傅里叶变换；

使用上述公式分别计算

的值。

3.如权利要求1所述的基于混合高阶偏微分方程模型的图像去噪方法，其特征在于，所述混合高阶偏微分方程模型的边缘检测函数的构造方法如下：

模型的特殊形式是整数阶导数，相关的泛函问题是：

式中，u为图像灰度函数，

表示u的梯度矢量模值即

用来检测边缘：当图像存在边缘时，梯度值较大，在平滑部分，梯度值较小，对平滑区域有抑制作用；u₀:Ω→R是采集到的含噪图像，

是具有Lipschitz边界的有界域；λ是大于零的参数；

D²u表示图像u的二阶导数，即

是具有Lipschitz边界的有界域，并且μ＞0,λ＞0是平衡目标函数中的参数；

h₁(x,y)是一个边界检测器，并由：

其中K₁和K₂是两个正常数；

考虑最小化问题：

式(7)对应的欧拉-拉格朗日方程如下：

是梯度向量，方向导数变化最快的方向，

是以梯度向量为速度场的散度，边界条件为：

使用最速下降法，导出相关的热流，得到最小化问题：

提出混合高阶偏微分方程：

进行参数选取后修改方程为：

式中，h(x,y)是边界检测器，定义为：

式中，ε为很小的一个正数，目的是确保边缘检测器小于1。

Images