CN104123701B - 基于莫罗包络平滑l1/全变差的图像脉冲噪声去除方法 - Google Patents

Abstract

基于莫罗包络平滑l₁/全变差的图像脉冲噪声去除方法，建立去除脉冲噪声的数学模型；对数学模型进行全变差范数的平滑化处理；对平滑化处理后的模型进行求解，得到去噪后的图像；对模型求解算法进行收敛速度分析，确定拥有加速步骤的求解算法流程。本发明基于莫罗包络平滑l₁/全变差的图像脉冲噪声去除方法，利用莫罗包络对全变差范数进行平滑化处理，平滑后的全变差范数为可微凸函数，是原函数的可微紧下界，具有迭代形式的解析解，也是原函数的解，其梯度很容易计算；由于求解的外部循环属于收缩算子，提出了加速策略，给出了加速版的算法，大大提高了收敛速度。

Description

基于莫罗包络平滑l1/全变差的图像脉冲噪声去除方法

技术领域

本发明属于信号处理技术领域，具体涉及一种基于莫罗包络平滑l₁/全变差的图像脉冲噪声去除方法。

背景技术

图像在获取、传输、存储及处理过程中，由于成像质量、传输介质、存储器件和处理手段的影响，不可避免地出现一些随机、离散的亮点或暗点，称之为脉冲噪声。脉冲噪声属于非全局噪声，去除比较困难。例如电荷耦合器件(CCD)因感光面积太小、热稳定性不足、长时间曝光等原因会产生脉冲噪声。为了去除脉冲噪声，一方面需要改善硬件质量，这需要较高的成本；另一方面可以从算法入手，设计合适的算法可以大幅度提高图像质量。

目前，去除图像脉冲噪声的基本方法是中值滤波和均值滤波。均值滤波是指对像素某个邻域内的所有像素值求平均值，用以代替原像素值。中值滤波则是选择像素某个邻域内的所有像素值的中间值，用中间值来代替原像素值。这两种方法只在脉冲噪声密度较小时有效，且容易造成局部模糊。维纳滤波法也可用来去除脉冲噪声，但它只涉及图像信息中的二阶统计特性，不涉及更重要高阶统计特性，因此其去噪效果并不理想。

发明内容

本发明的目的是提供一种基于莫罗包络平滑l₁/全变差的图像脉冲噪声去除方法，该方法能够去除较高密度的脉冲噪声，修复前后图像对比度和形 态不变，不易产生局部模糊。

本发明所采用的技术方案是，基于莫罗包络平滑l₁/全变差的图像脉冲噪声去除方法，具体包括以下步骤：

步骤1，建立去除脉冲噪声的数学模型；

步骤2，对步骤1建立的数学模型中的全变差范数进行平滑化处理；

步骤3，对步骤2平滑化处理后的模型进行求解，得到去噪后的图像；

步骤4，对步骤3模型求解算法进行收敛速度分析，确定拥有加速步骤的求解算法流程。

本发明的特点还在于，

步骤1具体为：

含脉冲噪声图像的数学模型为

b＝x+w (1)

其中，矩阵x、w和b分别表示尺寸为n×n的无噪声图像、脉冲噪声和退化图像；

去除脉冲噪声的数学模型为

min||x||_TV

s.t||x-b||₀≤ε (2)

为便于表达，设r＝n²，将x看成向量，R为n×n的一阶差分Toeplitz方阵即

$R=\left[\begin{array}{cccc}1& & & \\ -1& 1& & \\ & ...& ...& \\ & & -1& 1\end{array}\right]---\left(3\right)$

令 $C=\left[\begin{array}{c}{I}_{n\×n}\⊗R\\ R\⊗I_{n\×n}\end{array}\right],$ I_{n × n}为单位矩阵，表示矩阵的Kronecker积，易知||x||_TV＝Φ(Cx)，因此最小化||x||_TV就是降低噪声的影响；

由于||x-b||₀失去了幅度特征，难以求导，公式(2)无法求解，又因l₁-范数对l₀-范数具有良好的近似特性，故选择l₁-范数作为保真项，公式(2)转化为：

minΦ(Cx)

s.t||x-b||₁≤ε (4)

正则化后可得

minF(x)＝λ||x-b||₁+Φ(Cx)

s.tx＝{x_i,j,0≤x_i,j≤1} (5)。

步骤2中平滑化处理的具体步骤为：

定义任意凸函数f(x)的莫罗包络为

${\mathrm{env}}_{\mathrm{af}}\left(x\right)=\underset{y}{\inf }[\frac{{\left|\right|y-x\left|\right|}_2^2}{2\mathrm{\α}}+f\left(y\right)]---\left(6\right)$

由凸分析知，任何凸函数的莫罗包络都是可微凸函数，且inf env_{α f}(x)＝f(x)，即凸函数的莫罗包络是其本身的紧下界，因此 同时，当α趋于0时，env_{α f}(x)也趋于f(x)；env_{α f}(x)的梯度为▽env_{α f}(x)＝α^-1[1-prox_{α f}(x)]，其中

为构造全变差范数的莫罗包络，需利用公式(5)的增广拉格朗日(augmented Lagrange)形式，即引入中间变量q，使公式(5)转化为

${\min }_F(q,x)=\mathrm{\λ}{\left|\right|x-b\left|\right|}_1+\mathrm{\Φ}\left(q\right)+\frac{1}{2\mathrm{\α}}{\left|\right|a-\mathrm{Cx}\left|\right|}_2^2$

s.tx＝{x_i,j,0≤x_i,j≤1} (7)

将(7)式中替换为Φ(Cx)的莫罗包络，可得

minF_env(Cx,x)＝λ||x-b||₁+env_αΦ(Cx)

s.tx＝{x_i,j,0≤x_i,j≤1} (8)。

步骤3，求解的具体算法流程为：

输入：含噪图像b，常数m、λ和α；

初始化：令x⁰＝b；

迭代： $x^k=b+{\mathrm{prox}}_{\frac{\mathrm{\λ}{\left|\right|\·\left|\right|}_1}{m}}[x^{k-1}-\frac{C^T}{\mathrm{m\α}}(I-{\mathrm{prox}}_{\mathrm{\α\Φ}})\left({\mathrm{Cx}}^k\right)-b],$

其中 ${\mathrm{prox}}_{m^{-1}\mathrm{\λ}{\left|\right|\·\left|\right|}_1}\left(x\right)=\max \left(\right|x|-m^{-1})\×\mathrm{sign}\left(x\right),{\mathrm{prox}}_{\mathrm{\α}{\left|\right|\·\left|\right|}_2}(\·)=\max ({\left|\right|\·\left|\right|}_2-\mathrm{\α})\×\frac{(\·)}{{\left|\right|\·\left|\right|}_2};$

输出：x^k+1。

步骤4，收敛速度分析如下：

每循环一次，内部循环l次prox_αΦ(·)，则总的收敛速度为O[l(1+s)/k]。其中k为的迭代次数，l为prox_αΦ(·)的迭代次数，s为小于1的常数；

拥有加速步骤的求解算法流程为：

输入：含噪图像b，常数m、λ和α；

初始化：令y⁰＝b；t⁰＝1；

迭代：

$x^k=b+{\mathrm{prox}}_{\frac{\mathrm{\λ}{\left|\right|\·\left|\right|}_1}{m}}[y^{k-1}-\frac{C^T}{\mathrm{m\α}}(I-{\mathrm{prox}}_{\mathrm{\α\Φ}})\left({\mathrm{Cy}}^k\right)-b];$

$t^{k+1}=\frac{1}{2}(\sqrt{1+4{\left(t^k\right)}^2}+1),y^{k+1}=x^k+\frac{t^k-1}{t^{k+1}}(x^k-x^{k-1});$

输出：x^k+1。

本发明的有益效果是，本发明基于莫罗包络平滑l₁/全变差的图像脉冲噪声去除方法，利用莫罗包络对全变差范数进行平滑化处理，平滑后的全变差范数为可微凸函数，是原函数的可微紧下界，具有迭代形式的解析解，也是原函数的解，其梯度很容易计算；由于求解的外部循环属于收缩算子(shrinkage operator)，提出了加速策略，给出了加速版的算法，大大提高了收敛速度。

附图说明

图1为1024×1024标准lena图像；

图2为图1添加密度为10％的脉冲噪声的lena图像；

图3为采用本发明方法对图2进行恢复的lena图像；

图4为图3的脸部放大图；

图5为1024×1024标准cameraman图像；

图6为图5加密度为10％的脉冲噪声的cameraman图像；

图7为采用本发明方法对图6进行恢复的cameraman图像；

图8为图7的局部放大图；

图9为图1添加密度为20％的脉冲噪声的lena图像；

图10为采用本发明方法对图9进行恢复的lena图像；

图11为图10的脸部放大图；

图12为图1添加密度为30％的脉冲噪声的lena图像；

图13为采用本发明方法对图12进行恢复的lena图像；

图14为图13的脸部放大图；

图15为四种不同算法对添加不同密度脉冲噪声的1024×1024标准lena 图像进行去噪处理的峰值信噪比；

图16为四种不同算法对添加不同密度脉冲噪声的1024×1024标准cameraman图像进行去噪处理的峰值信噪比。

具体实施方式

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

本发明基于莫罗包络平滑l₁/全变差的图像脉冲噪声去除方法，具体包括以下步骤：

步骤1，建立去除脉冲噪声的数学模型；

含脉冲噪声图像的数学模型为

b＝x+w (1)

其中，矩阵x、w和b分别表示尺寸为n×n的无噪声图像、脉冲噪声和退化图像；

去除脉冲噪声的数学模型为

min||x||_TV

s.t||x-b||₀≤ε (2)

为便于表达，设r＝n²，将x看成向量，R为n×n的一阶差分Toeplitz方阵即

$R=\left[\begin{array}{cccc}1& & & \\ -1& 1& & \\ & ...& ...& \\ & & -1& 1\end{array}\right]---\left(3\right)$

令 $C=\left[\begin{array}{c}{I}_{n\×n}\⊗R\\ R\⊗I_{n\×n}\end{array}\right],$ I_{n × n}为单位矩阵，表示矩阵的Kronecker积，易知||x||_TV＝Φ(Cx)，因此最小化||x||_TV就是降低噪声的影响；

由于||x-b||₀失去了幅度特征，难以求导，公式(2)无法求解，又因l₁-范数对l₀-范数具有良好的近似特性，故选择l₁-范数作为保真项，公式(2)转化为：

minΦ(Cx)

s.t||x-b||₁≤ε (4)

正则化后可得

minF(x)＝λ||x-b||₁+Φ(Cx)

s.tx＝{x_i,j,0≤x_i,j≤1} (5)；

步骤2，对步骤1建立的数学模型中的全变差范数的进行平滑化处理；

平滑化处理的具体步骤为：

定义任意凸函数f(x)的莫罗包络为

${\mathrm{env}}_{\mathrm{af}}\left(x\right)=\underset{y}{\inf }[\frac{{\left|\right|y-x\left|\right|}_2^2}{2\mathrm{\α}}+f\left(y\right)]---\left(6\right)$

由凸分析知，任何凸函数的莫罗包络都是可微凸函数，且inf env_{α f}(x)＝f(x)，即凸函数的莫罗包络是其本身的紧下界，因此 同时，当α趋于0时，env_{α f}(x)也趋于f(x)；env_{α f}(x)的梯度为▽env_{α f}(x)＝α^-1[1-prox_{α f}(x)]，其中

为构造全变差范数的莫罗包络，需利用公式(5)的增广拉格朗日(augmented Lagrange)形式，即引入中间变量q，使公式(5)转化为

${\min }_F(q,x)=\mathrm{\λ}{\left|\right|x-b\left|\right|}_1+\mathrm{\Φ}\left(q\right)+\frac{1}{2\mathrm{\α}}{\left|\right|a-\mathrm{Cx}\left|\right|}_2^2$

s.tx＝{x_i,j,0≤x_i,j≤1} (7)

将(7)式中替换为Φ(Cx)的莫罗包络，可得

minF_env(Cx,x)＝λ||x-b||₁+env_αΦ(Cx)

s.tx＝{x_i,j,0≤x_i,j≤1} (8)；

步骤3，对步骤2平滑化处理后的模型进行求解，得到去噪后的图像；

求解的具体算法流程为：

输入：含噪图像b，常数m、λ和α；

初始化：令x⁰＝b；

迭代： $x^k=b+{\mathrm{prox}}_{\frac{\mathrm{\λ}{\left|\right|\·\left|\right|}_1}{m}}[x^{k-1}-\frac{C^T}{\mathrm{m\α}}(I-{\mathrm{prox}}_{\mathrm{\α\Φ}})\left({\mathrm{Cx}}^k\right)-b],$

其中 ${\mathrm{prox}}_{m^{-1}\mathrm{\λ}{\left|\right|\·\left|\right|}_1}\left(x\right)=\max \left(\right|x|-m^{-1})\×\mathrm{sign}\left(x\right),{\mathrm{prox}}_{\mathrm{\α}{\left|\right|\·\left|\right|}_2}(\·)=\max ({\left|\right|\·\left|\right|}_2-\mathrm{\α})\×\frac{(\·)}{{\left|\right|\·\left|\right|}_2};$

输出：x^k+1；

步骤4，对步骤3模型求解算法进行收敛速度分析，确定拥有加速步骤的求解算法流程；

收敛速度分析如下：

每循环一次，内部循环l次prox_αΦ(·)，则总的收敛速度为O[l(1+s)/k]。其中k为的迭代次数，l为prox_αΦ(·)的迭代次数，s为小于1的常数；

拥有加速步骤的求解算法流程为：

输入：含噪图像b，常数m、λ和α；

初始化：令y⁰＝b；t⁰＝1；

迭代：

$x^k=b+{\mathrm{prox}}_{\frac{\mathrm{\λ}{\left|\right|\·\left|\right|}_1}{m}}[y^{k-1}-\frac{C^T}{\mathrm{m\α}}(I-{\mathrm{prox}}_{\mathrm{\α\Φ}})\left({\mathrm{Cy}}^k\right)-b];$

$t^{k+1}=\frac{1}{2}(\sqrt{1+4{\left(t^k\right)}^2}+1),y^{k+1}=x^k+\frac{t^k-1}{t^{k+1}}(x^k-x^{k-1});$

输出：x^k+1。

通过以下仿真数据和图像进一步说明本发明的优点。

1、仿真条件

仿真中采用1024×1024的标准lena图像和cameraman图像作为仿真对象。

2、仿真内容

仿真1，分别对1024×1024的标准lena和cameraman图像添加密度为10％的脉冲噪声，用本发明方法中拥有加速步骤的算法对其进行去噪处理，参见图1至图8。

仿真2，对1024×1024的标准lena图像添加密度为20％的脉冲噪声，用本发明方法中拥有加速步骤的算法对其进行去噪处理，参见图9至图11。

仿真3，对1024×1024的标准lena图像添加密度为30％的脉冲噪声，用本发明方法中拥有加速步骤的算法对其进行去噪处理，参见图12至图14。

仿真4，分别用本发明方法中拥有加速步骤的算法、非解耦形式的l₁/TV方法、中值滤波法、均值滤波法对添加不同密度脉冲噪声的1024×1024标准lena图像进行去噪处理。处理结果参见图15。

仿真5，分别用本发明方法中拥有加速步骤的算法、非解耦形式的l₁/TV方法、中值滤波法、均值滤波法对添加不同密度脉冲噪声的1024×1024标准cameraman图像进行去噪处理。处理结果参见图16。

3、仿真结果分析

根据仿真1，计算得修复后lena图像PSNR(峰值信噪比)达到36.21dB， 修复图像与原图像几乎完全一致，脸部特征和头发及帽沿等细节和纹理保持得很好，整幅图像轮廓鲜明，噪声去除得非常干净，没有产生局部模糊现象。修复后的cameraman图像PSNR达到34.34dB，修复图像与原图像也几乎完全一致，脸部特征和相机及手套等细节保持得很好，噪声也去除得非常干净，没有产生局部模糊现象。

根据仿真2，计算得修复后图像的PSNR为34.56dB。从图10可以看出，修复效果依然很理想，图像的整体轮廓未受影响，纹理和细节得也到了很好的保持，修复图像与原图像几乎完全一致，没有噪声残留，几乎没有产生局部模糊现象。

根据仿真3，计算得修复后图像的PSNR为32.71dB。相比仿真1和2的结果，图13中出现了局部模糊现象，但并不是很明显，说明即使在高密度噪声环境下本发明方法也能够取得较好的恢复效果。产生局部模糊是因为当脉冲噪声密度较大时，其能量较大，此时最小化(5)式会迫使全变差范数值低于理想值，使原本不需平滑的区域被平滑以满足最小化(5)式。所以，局部模糊产生的区域为像素值梯度较大的区域，如边缘、纹理细节处。

根据仿真4和仿真5，由图15和16可以看出，四种方法中本发明方法拥有最高的峰值信噪比，并且随着脉冲噪声密度的增加，峰值信噪比的下降最为平缓，表现最为稳定。非解耦形式的l₁/TV方法稍逊于本发明方法，峰值信噪比大约低0.5-1dB，下降速度也很平缓。中值滤波法很不稳定，当脉冲噪声密度较小时，其峰值信噪比与前两种方法差不多，甚至还略高一点，但随着噪声密度的增大，中值滤波法的峰值信噪比迅速下降，当噪声密度高于40％时，其性能甚至差于均值滤波法。

Claims (1)

1.基于莫罗包络平滑l₁/全变差的图像脉冲噪声去除方法，其特征在于，具体包括以下步骤：

步骤1，建立去除脉冲噪声的数学模型；

步骤1具体为：

含脉冲噪声图像的数学模型为

b＝x+w (1)

其中，矩阵x、w和b分别表示尺寸为n×n的无噪声图像、脉冲噪声和退化图像；

去除脉冲噪声的数学模型为

min||x||_TV

s.t||x-b||₀≤ε (2)

为便于表达，设r＝n²，将x看成向量，R为n×n的一阶差分Toeplitz方阵即

$R=\left[\begin{array}{cccc}1& & & \\ -1& 1& & \\ & ...& ...& \\ & & -1& 1\end{array}\right]---\left(3\right)$

令I_n×n为单位矩阵，表示矩阵的Kronecker积，易知||x||_TV＝Φ(Cx)，因此最小化||x||_TV就是降低噪声的影响；

由于||x-b||₀失去了幅度特征，难以求导，公式(2)无法求解，又因l₁-范数对l₀-范数具有良好的近似特性，故选择l₁-范数作为保真项，公式(2)转化为：

minΦ(Cx)

s.t||x-b||₁≤ε (4)

正则化后可得

min F(x)＝λ||x-b||₁+Φ(Cx)

s.t x＝{x_i,j,0≤x_i,j≤1} (5)

步骤2，对步骤1建立的数学模型中的全变差范数进行平滑化处理；

步骤2中平滑化处理的具体步骤为：

定义任意凸函数f(x)的莫罗包络为

${\mathrm{env}}_{\mathrm{\α}f}\left(x\right)=\underset{y}{inf}\[\frac{\left|\right|y-x|{|}_2^2}{2\mathrm{\α}}+f\left(y\right)\]---\left(6\right)$

由凸分析知，任何凸函数的莫罗包络都是可微凸函数，且inf env_αf(x)＝f(x)，即凸函数的莫罗包络是其本身的紧下界，因此同时，当α趋于0时，env_αf(x)也趋于f(x)；env_αf(x)的梯度为其中

为构造全变差范数的莫罗包络，需利用公式(5)的增广拉格朗日形式，即引入中间变量q，使公式(5)转化为

$minF(q,x)=\mathrm{\λ}\left|\right|x-b|{|}_1+\mathrm{\Φ}\left(q\right)+\frac{1}{2\mathrm{\α}}||q-Cx|{|}_2^2$

s.t x＝{x_i,j,0≤x_i,j≤1} (7)

将(7)式中替换为Φ(Cx)的莫罗包络，可得

min F_env(Cx,x)＝λ||x-b||₁+env_αΦ(Cx)

s.t x＝{x_i,j,0≤x_i,j≤1} (8)；

步骤3，对步骤2平滑化处理后的模型进行求解，得到去噪后的图像；

步骤3，求解的具体算法流程为：

输入：含噪图像b，常数m、λ和α；

初始化：令x⁰＝b；

迭代：

其中

输出：x^k+1；

步骤4，对步骤3模型求解算法进行收敛速度分析，确定拥有加速步骤的求解算法流程；

步骤4，收敛速度分析如下：

每循环一次，内部循环l次prox_αΦ(·)，则总的收敛速度为O[l(1+s)/k]，其中k为的迭代次数，l为prox_αΦ(·)的迭代次数，s为小于1的常数；

拥有加速步骤的求解算法流程为：

输入：含噪图像b，常数m、λ和α；

初始化：令y⁰＝b；t⁰＝1；

迭代：

$x^k=b+{\mathrm{prox}}_{\frac{\mathrm{\λ}\left|\right|\·|{|}_1}{m}}\[y^{k-1}-\frac{C^T}{m\mathrm{\α}}(I-{\mathrm{prox}}_{\mathrm{\α}\mathrm{\Φ}})\left({\mathrm{Cy}}^k\right)-b\];$

$\begin{array}{cc}{t}^{k+1}=\frac{1}{2}(\sqrt{1+4{\left(t^k\right)}^2}+1)& y^{k+1}=x^k+\frac{t^k-1}{t^{k+1}}(x^k-x^{k-1})\end{array};$

输出：x^k+1。