图像在获取和传输过程中, 不可避免地受到各种噪声的污染, 导致图像的真实信息被噪声所淹没, 降低图像的可读性与可应用性.因此我们有必要对含噪图像进行降噪处理, 使得复原图像能够满足图像分割、目标识别和边缘检测等后续工作的需要.此时, 描述含噪图像$f$的数学模型为

$ f=u+n $

(1.1)

其中$u$为清晰图像, $n$为方差为$\delta^2$的高斯白噪声.图像去噪已成为数字图像处理领域中一个非常重要的课题, 目前已有领域滤波、频域滤波和小波域滤波等多种去噪方法^[1].在$L^2$范数基础上构建的各向同性扩散模型, 由于其在去噪过程中的各项同性扩散特点, 使得复原图像中的边缘轮廓容易被模糊^[1].1992年Rudin、Osher和Fatemi^[2]在$L^1$范数基础上提出了全变差正则化(Total variation, TV)模型, 该模型在消除噪声的同时能够很好地保留图像的边缘轮廓, 但会将连续光滑信号变成分段的等值信号, 从而导致“分片常数”效应^{[3, 4]}, Song^[5]为克服TV模型容易产生“分片常数”效应的问题, 提出了广义TV模型, 但广义参数$p$的客观选取严重影响了含噪图像的复原质量.鉴于此, 张红英和黎芳等人^{[6, 7]}提出了自适应的广义TV模型, 该模型根据含噪图像的局部梯度信息自适应地选取广义参数$p$, 能够同时兼顾图像的平滑去噪与边缘保留.但根据局部梯度信息自适应选取的广义参数$p$对噪声有着很强的敏感性, 在图像受噪声污染较为严重时不能准确地得到广义参数$p$, 降低了图像复原质量^[8].

为提高含噪图像的复原质量, 本文先利用高斯滤波器对含噪图像作预处理, 然后构建8个边缘检测算子在预处理图像基础上计算广义参数$p$, 降低广义参数自适应选取过程中对噪声的敏感性, 以构建改进的广义全变差正则化图像去噪模型.该模型在消除噪声的同时能够保留图像的边缘轮廓等细节信息, 提高图像复原的质量.

受上述研究内容的启发, 文章提出如下的广义全变差正则化图像去噪模型(Generalized total variation, GTV):

$ \min\limits_{u\in BV(\Omega)\cap{L^2}(\Omega)}E(u)=\iint\limits_{\Omega}\frac{1}{p(x, y)}|\nabla u|^{p(x, y)}dxdy+\frac{\lambda}{2}\iint\limits_{\Omega}|u - f|^2dxdy, $

(2.1)

其中$\lambda>0$为正则化参数, 广义参数$p(x,y)\in[1,2]$随图像局部结构信息的变化而变化, $\Omega$为图像空间且$(x,y)\in\Omega$.

不失一般性, 后文讨论中将$p(x, y)$简记为$p$.当$p=1$时, 广义全变差正则化模(GTV)退化成基于$L^1$范数的全变差正则化去噪模型.此时的模型沿图像边缘轮廓方向扩散能力较弱, 能够保留图像的边缘轮廓信息, 但易将噪声当成边缘从而在图像的平滑区产生“分片常数”效应.当$p=2$时, GTV等效成基于$L^2$范数的各向同性扩散(调和)去噪模型, 在平滑区能够有效地去除噪声, 但因该模型朝各个方向的扩散能力相同, 容易模糊图像的边缘轮廓等细节信息.于是构建广义全变差正则化去噪模型的重点便转化成如何针对图像的局部结构信息自适应地选取广义参数$p$, 以达到在平滑噪声的同时能够保留图像的边缘轮廓等细节信息.

为降低广义参数$p$对噪声的敏感性, 文章先用一个高斯滤波器$G_\sigma$对含噪图像作预处理, 消除含噪图像平滑区的部分噪声, 降低将噪声当成图像假边缘的可能性.然后构建8个边缘检测算子$d_\theta$(大小为$6\times6$)在预处理图像的基础上计算广义参数$p$, 其中预处理图像$\hat{f}$为

$ \hat{f}=G_\sigma\otimes f, $

(2.2)

式中$\sigma>0$为高斯滤波器$G_\sigma$的标准差, $\otimes$为卷积运算符.此时的边缘检测算子$d_\theta$是由$d_0$经旋转角度$\theta\in\Theta$ ($\Theta=\{0, \pi/4, \pi/2, 3\pi/4, \pi, 5\pi/4, 3\pi/2, 7\pi/4\}$)后得到的, 其中边缘检测算子$d_0$为(见文献[9])

$   {d_0} = \frac{1}{{12}}\left[{\begin{array}{*{20}{c}} {}   {{{\rm{O}}_{\rm{1}}}}   {} \\ {{{\rm{O}}_{\rm{2}}}}   {\rm{M}}   {{{\rm{O}}_{\rm{2}}}} \\ {}   {{{\rm{O}}_{\rm{1}}}}   {} \\ \end{array}} \right], {\rm{M}} = \left[{\begin{array}{*{20}{c}} {\;\;1}   {\;\;2}   {\;\;2}   {\;\;1} \\ {-1}   {-2}   {-2}   { - 1} \\ \end{array}} \right], \nonumber\\    {{\rm{O}}_1} = \left[{\begin{array}{*{20}{c}} 0   0   0   0   0   0 \\ 0   0   0   0   0   0 \\ \end{array}} \right], {{\rm{O}}_2} = \left[{\begin{array}{*{20}{c}} 0 \\ 0 \\ \end{array}} \right]. $

(2.3)

根据对边缘检测算子$d_\theta$的定义, 广义全变差正则化去噪模型(2.1) 中的广义参数$p$为

$ p = 1 + \frac{1}{{1 + \frac{1}{2} \times \sqrt {\sum\limits_{\theta \in \Theta } {{{({d_\theta } \otimes \hat f)}^2}} } }} \in [1,2],~~~\Theta = \left\{ {0,\frac{\pi }{4},\frac{\pi }{2},\frac{{3\pi }}{4},\pi ,\frac{{5\pi }}{4},\frac{{3\pi }}{2},\frac{{7\pi }}{4}} \right\}. $

(2.4)

由于边缘检测算子$d_\theta$有着很强的对称性, 能够同时利用图像在8个方向上的灰度信息, 降低噪声对边缘检测结果的影响, 提高广义参数$p$计算的鲁棒性。为验证本文所定义的广义参数$p$的有效性, 以受标准差为10的高斯白噪声污染的Lena图像为例, 与文献[4, 5]中对$p$的计算结果进行了比较分析, 如图1所示

如图 1(a)所示, 清晰Lena图像中的边缘轮廓和光滑区域清晰可见, 文献[4, 5]在计算广义参数$p$时仅仅利用了预处理图像$\hat{f}$的局部梯度信息, 无法克服对噪声的敏感性, 如图 1(b)中计算的广义参数$p$出现了较大的误差.在图 1(c)中, 本文采用对称性较强的边缘检测算子$d_\theta$, 降低了广义参数$p$计算过程中对噪声的敏感性, 提高了广义参数$p$的计算精度.即在图像的平滑区域$\sum\limits_{\theta\in\Theta} {{{({d_\theta}\otimes\hat{f})}^2}}\to0$, 此时广义参数$p\to2$, GTV模型能够有效地平滑噪声; 在图像的边缘轮廓区域$\sum\limits_{\theta\in\Theta}{{{({d_\theta}\otimes\hat{f})}^2}}\to\infty$, 广义参数$p\to1$, 此时GTV模型沿图像边缘轮廓方向扩散能力较弱, 能够保留图像的边缘轮廓信息.在图像的其它区域, 广义参数$p$按照$\sum\limits_{\theta\in\Theta}{{{({d_\theta}\otimes\hat{f})}^2}}\to\infty$的大小自适应地变化, 同时兼顾图像的平滑去噪与边缘轮廓保留, 提高含噪图像的复原质量.

为了求出能量泛函(2.1) 的极小解, 本文先利用变分法推导该能量泛函对应的Euler-Lagrange偏微分方程(Partial differential equations, PDEs), 然后得到该PDEs对应的梯度下降流方程, 此时梯度下降流方程的稳定解即为对应能量泛函的极小解。

令$F(x, y, u, {u_x}, {u_y})=\frac{1}{p}|\nabla u|^p+\frac{\lambda}{2}|u-f|^2$, 此时对广义全变差正则化模型(2.1) 的最优解$u(x, y)$作一微扰动$v(x, y)$, 得$u(x, y)+v(x, y)$, 当$v(x, y)$和$v'(x, y)$足够小时, 由Taylor展开得到

$ F(x, y, u + v, {u_x} + {v_x}, {u_y} + {v_y}) = F(x, y, u, {u_x}, {u_y}) + \frac{{\partial F}}{{\partial u}}v + \frac{{\partial F}}{{\partial {u_x}}}{v_x} + \frac{{\partial F}}{{\partial {u_y}}}{v_y} + \cdots. $

(3.1)

联立模型(2.1) 和公式(3.1), 得到

$ E(u + v) \cong E(u)+\iint\limits_{\Omega}\left(\frac{{\partial F}}{{\partial u}}v + \frac{{\partial F}}{{\partial {u_x}}}{v_x} + \frac{{\partial F}}{{\partial {u_y}}}{v_y}\right)dxdy. $

(3.2)

根据分部积分法和Neumann边界条件, 得到

$ E(u + v) \cong E(u) + \iint\limits_{\Omega}{\left\{ {v \cdot \left[{\lambda (u-f)-\frac{d}{{dx}}\left( {\frac{{\partial F}}{{\partial {u_x}}}} \right)-\frac{d}{{dy}}\left( {\frac{{\partial F}}{{\partial {u_y}}}} \right)} \right]} \right\}}dxdy. $

(3.3)

当能量泛函$E(u)$达到极值, 对$u(x, y)$的任一足够小的微小扰动$v(x, y)$都应使得$E(u)$的极值不变, 则此时可得到能量泛函(2.1) 的Euler-Lagrange偏微分方程^[10]

$ \lambda (u - f) - \frac{d}{{dx}}\left( {\frac{{\partial F}}{{\partial {u_x}}}} \right) - \frac{d}{{dy}}\left( {\frac{{\partial F}}{{\partial {u_y}}}} \right) = 0. $

(3.4)

又因为$\partial{F}/\partial{u_x}=|\nabla u|^{p-2}{u_x}$, $\partial{F}/\partial{u_y}=|\nabla u|^{p-2}{u_y}$, 则得到最终的Euler-Lagrange偏微分方程为

$ \lambda (u - f) - \nabla \cdot \frac{{\nabla u\;\;\;}}{{|\nabla u{|^{2 - p}}}} = 0. $

(3.5)

此时能量泛函(2.1) 的极值问题归结为求解相应Euler-Lagrange偏微分方程问题, 但该非线性PDEs的求解非常困难.同时, 为抑制方程(3.5) 在图像平滑区域($\nabla u \approx0$)不稳定的问题, 引入一较小的常数$\beta>0$, 将方程(3.5) 改写成如下形式:

$ \lambda (u - f) - \nabla \cdot \frac{{\nabla u\;\;\;}}{{|\nabla u|_\beta ^{2 - p}}} = 0, $

(3.6)

其中$|\nabla u{|_\beta } = \sqrt {|\nabla u{|^2} + \beta }$.此时在该PDEs中引入一个时间辅助变量t, 将求解静态非线性PDEs问题转化成一个动态PDEs问题, 当该方程的迭代演化过程达到稳态时, 便可得到Euler-Lagrange偏微分方程的解.而该动态PDEs方程即为负梯度下降流方程, 此时图像去噪的迭代演化方程为

$ \left\{ \begin{array}{l} \frac{{\partial u}}{{\partial t}} = \nabla \cdot \frac{{\nabla u\;\;\;}}{{|\nabla u{|^{2 - p}}}} - \lambda (u - f), \\ {\left. u \right|_{t = 0}} = I, \\ {\left. {\frac{{\partial u}}{{\partial n}}} \right|_{\partial \Omega }} = 0, \end{array} \right. $

(3.7)

其中${\left. {\frac{{\partial u}}{{\partial n}}} \right|_{\partial \Omega }} = 0$为Neumann边界条件, $n$为图像$u$的区域边界.然后利用有限差分格式去离散上述的演化方程, 可得到对应的迭代求解格式^[1].利用散度的基本性质有

$ \begin{eqnarray}\label{eq:3.8} \nabla\cdot\frac{\nabla u}{|\nabla u|^{2 - p}}&=&(\nabla\cdot\nabla u)\cdot|\nabla u|^{p - 2}+\nabla(|\nabla u|^{p - 2})\cdot\nabla u\\ &=&|\nabla u|^{p - 4} \cdot\left[(p-1)(u_x^2{u_{xx} + u_y^2u_{yy}) + 2(p-2){u_x}{u_y}u_{xy} + u_x^2u_{yy} + u_y^2u_{xx}} \right].\nonumber \end{eqnarray} $

(3.8)

令空间步长为$h=1$, 时间步长为$\Delta t$, 变量梯度利用中心差分进行计算, 此时迭代演化方程(3.7) 的离散迭代求解格式为

$ {u_{k + 1}} = {u_k} + \Delta t \cdot \left[{\Phi ({u_k})-\lambda ({u_k}-f)} \right], $

(3.9)

其中

$ \Phi ({u_k}) = |\nabla {u_k}{|^{p -4}} \cdot \left[{(p-1)(u_{k, x}^2{u_{k, xx}} + u_{k, y}^2{u_{k, yy}}) + 2(p-2){u_{k, x}}{u_{k, y}}{u_{k, xy}} + u_{k, x}^2{u_{k, yy}} + u_{k, y}^2{u_{k, xx}}} \right], $

迭代初始值取${u_0}=f$.

为验证文中广义全变差正则化图像去噪模型(GTV)的有效性, 采用标准的Lena和Toys图像进行实验, 并与标准全变差正则化模型(Standard Total Variation, STV)和自适应全变差正则化模型(Adaptive Total Variation, ATV)的实验结果进行比较分析.实验仿真中, 输入图像的大小均为$256\times256$, 灰度级为256, 加性噪声的标准差分别为20和30.实验是在CPU Intel(R) 2.0GHz, 1G内存, Matlab 7.9.0的环境下进行, 综合采用传统的量化指标峰值信噪比(Peak signal-to-noise ratio, PSNR)和最近受关注的平均结构相似度^{[11, 12]}(Mean structural similarity, MSSIM)作为图像去噪质量的评价标准, MSSIM指标用以衡量复原图像与清晰图像在亮度、对比度和结构三方面的综合相似性, 其中PSNR和MSSIM的定义分别为

$ {\rm{PSNR}}(u, \bar u) = 10{\log _{10}}\frac{{{{255}^2} \times M \times N}}{{||u - \bar u||_2^2}}{\rm{dB}}, $

(4.1)

$ {\rm{MSSIM}}(u, \bar u) = \frac{1}{\Sigma }\sum\limits_{i = 1}^\Sigma {{{[l({u_i}, {{\bar u}_i})]}^\alpha } \cdot {{[c({u_i}, {{\bar u}_i})]}^\beta } \cdot {{[s({u_i}, {{\bar u}_i})]}^\gamma }}, $

(4.2)

其中$u$为清晰图像, $\bar{u}$为去噪图像, $M\times N$为图像大小, $\Sigma$为图像分块数目, ${\bar u_i}$与${u_i}$分别为去噪图像块与其对应的清晰图像块.本文取MSSIM指标中调整三个比较函数所占比重的参数$\alpha$, $\beta$与$\gamma$均为1, 广义全变差正则化模型(GTV)中的正则化参数为$\lambda=0.01$, 且离散迭代求解格式(3.8) 中的时间步长为$\Delta t=0.2$, 并满足停机准则${{||{u_k} - {u_{k - 1}}||_2^1} /{||{u_k}||_2^1}} < 5 \times {10^{-4}}$.此时的实验结果如图 2, 图 3和表1所示.

图 2和图 3分别为去噪前后的Lena和Toys图像, 表1为不同噪声水平(噪声标准差分别为20和30) 下的复原图像客观评价结果.由实验数据可知, 本文提出的广义全变差正则化模型(GTV)在不同噪声水平下的图像复原效果均优于其它两种方法, 人眼能直观地感受到图像复原质量的提高.在图 2(b)和图3(b)中, STV模型虽然在消除噪声的同时能够保留图像的边缘轮廓信息, 但在图像的平滑区有着明显的“分片常数”效应. ATV模型在一定程度上可以克服STV模型极易产生“分片常数”效应的不足, 但由于在其对应的广义参数$p$对噪声异常敏感, 使得在平滑区易将噪声当作边缘轮廓而保留, 导致复原图像出现了假边缘.而本文提出的GTV模型中的广义参数$p$对噪声有着极强的鲁棒性, 在消除图像平滑区噪声的同时, 能够很好地保留图像的边缘轮廓等结构信息, 提高了含噪图像的复原质量.GTV涉及8个方向检测算子的计算, 而ATV选取广义参数$p$的方式较为简单, 因此GTV在CPU计算时间上比ATV略有增加, 但较STV能够明显地提高图像去噪的计算效率.

文章在分析基于$L^1$范数的全变差正则化模型和基于$L^2$范数的调和模型的优点与不足的基础上, 提出了基于$L^p$($1\le p\le 2$)范数的广义全变差正则化图像去噪模型, 并引入边缘检测算子自适应选取广义参数$p$, 降低其对噪声的敏感性.实验结果表明, 本文提出的广义全变差正则化模型消噪效果明显, 并在消除噪声的同时能够保留图像的边缘轮廓等结构信息, 克服了全变差正则化模型中“分片常数”效应和调和模型中模糊边缘轮廓的不足, 且人眼能直观地感受到图像复原质量的提高.