从模糊的、有噪声的或者部分观测图像做图像恢复可以看成一个反问题的应用, 数学表达式可以写成如下形式

$ \begin{equation} y=Ax+ \epsilon, \label{eq1} \end{equation} $

(1.1)

其中$x$, $y \in R^{m \times n}$分别代表着原始图像与观察到的图像, $\epsilon$是方差为$\sigma^{2}$的高斯噪声. $A$是从$R^{m \times n}$映射到$R^{m \times n}$的线性算子, 在图像的去卷积过程中, 取$A$为卷积算子.在反问题中, $A$是不可逆或者有很大条件数的病态的算子.由于噪声的存在, 完全恢复图像是不可能的, 所以只能使用正则化方法寻找近似解.修改后的模型为

$ \begin{equation} \underset{\text{x}}{\min} \frac{1}{2} \left\| y - Ax \right\|_{2}^{2} + \alpha R(x), \label{eq2} \end{equation} $

(1.2)

其中$\alpha R(x)$就是包含着解的先验信息的正则项, $\alpha$是正则化参数.

我们在正则项中使用$ R(x) = | Dx | $, 其中$D$取局部差分算子.这个正则项是非光滑非可分的, 一个自然的想法是使用光滑化的$l_{1}$范数逼近, 基于全变分(TV)的图像恢复经常使用这种方法.而TV有利于保持图像的锐利和物体的边界, 而且在恢复图像过程中保住分片光滑的物体^{[1, 2]}.但是同时TV也有一些缺点, 比如会在平坦区域产生阶梯现象造成图像细节丢失等.

在图像恢复过程中另一个重要性质就是图像的稀疏表示.最好的稀疏性由$L_{0}$范数生成, 但$L_{0}$范数导致了NP困难问题, 所以我们使用了$L_{1}$范数做逼近.利用小波框架表示图像则有利于保证稀疏性的存在^{[3, 4]}.本文将TV和小波框架两种正则化模型结合在一起以保持各自的优点.在正则化模型的使用过程中, 用ADMM算法以达到快速稳定的目标.在之后的数值实验中, 将看到新方法提高了原有的TV和框架模型的效果^{[5, 6]}.

TV去噪可以看做最有效和最难计算的几种方法之一.为了去除它的缺点, 近年已经发展出了原对偶方法, $FTVd$和增广拉格朗日方法等^[7].在这里将TV方法应用于原问题.

对于图像$x \in R^{m\times n}$, 用符号$\nabla$表示离散阶梯算子从$R^{m\times n}$映射到$R^{m \times n} \times R^{m \times n}$.即$\nabla x = (\nabla_{1} x , \nabla_{2} x)$, 在这里

$ \begin{eqnarray} (\nabla_{1} x)_{i, j} = \left\{ \begin{aligned} x_{i+1, j} - x_{i, j}, ~~ &i=1, \cdots, M-1;j=1, \cdots, N, \\ x_{1, j} - x_{M, j}, ~~ &i=M;j=1, \cdots, N; \end{aligned} \right. \end{eqnarray} $

(2.1)

$ \begin{eqnarray} (\nabla_{2} x)_{i, j} = \left\{ \begin{aligned} x_{i, j+1} - x_{i, j}, ~~ &i=1, \cdots, M;j=1, \cdots, N-1, \\ x_{i, 1} - x_{i, N}, ~~ &i=1, \cdots, M;j=N. \end{aligned} \right. \end{eqnarray} $

(2.2)

由于使用的是$L_{1}$范数, 定义$\left\| x \right\|_{TV} = \left\| \nabla x \right\|_{1} = \left| \nabla_{1} x \right| + \left| \nabla_{2} x \right|$, 称为各向异性, 那么在应用ADMM算法时会有两个正则化项.此时问题变成

$ \begin{equation} \underset{x}{\min} \lambda \left| \nabla_{1} x \right| + \lambda \left| \nabla_{2} x \right| + \frac{1}{2} \left\| y - Ax \right\|_{2}^{2}, \label{eq7} \end{equation} $

(2.3)

这里取$\lambda > 0$作为正则化参数, 用来平衡三项的最小化取值.把该最小化问题变成约束最优化问题

$ \begin{eqnarray} \underset{x, u}{\min} \frac{1}{2} \left\| y - Ax \right\|_{2}^{2} + \lambda \left\| u \right\|_{1}~~~~ {\rm s.t. }\quad u = \nabla x, \label{eq12} \end{eqnarray} $

(2.4)

再得到增广拉格朗日函数

$ \begin{equation} L(x, u;p_{2}) = \frac{1}{2} \left\| y - Ax \right\|_{2}^{2} + \lambda \left\| u \right\|_{1} + \frac{a}{2}\left\| u - \nabla x \right\|_{2}^{2} + \langle p_{2}, u - \nabla x \rangle, \label{eq14} \end{equation} $

(2.5)

这里$\thinspace p_{2} = (p_{2, 1}, p_{2, 2}) \in R^{M \times N} \times R^{M \times N}$为拉格朗日乘子, 作用是防止$\lambda > 0$在计算中趋向于无穷大.

类似的, 取框架分解算子$W$, 那么可以把原反问题写成最小化问题

$ \begin{equation} \underset{x}{\min} \mu \left\| Wx \right\|_{1} + \frac{1}{2} \left\| y - Ax \right\|_{2}^{2}, \end{equation} $

(2.6)

这里取$\mu > 0$为正则化参数.把该最小化问题变成约束最优化问题

$ \begin{eqnarray} \underset{x, v}{\min} \frac{1}{2} \left\| y - Ax \right\|_{2}^{2} + \mu \left\| v \right\|_{1}~~~~ {\rm s.t.} \quad v = Wx, \end{eqnarray} $

(2.7)

再得到增广拉格朗日函数

$ \begin{equation} L(x, v;p_{1}) = \frac{1}{2} \left\| y - Ax \right\|_{2}^{2} + \mu \left\| v \right\|_{1} + \frac{b}{2}\left\| v - Wx \right\|_{2}^{2} + \langle p_{1}, v - Wx \rangle, \end{equation} $

(2.8)

这里$p_{1} \in R^{M \times N}$为拉格朗日乘子, 作用是防止$\mu >0$在计算中趋向于无穷大.

结合两个模型, 在用小波框架作为第一步求出一个最优解之后, 把该最优解作为一个初始值带入到第二个模型算法之中做第二次运算, 以求出更好的解.于是得到新的方法如下

1.1:输入初始值$ y, A, \mu > 0, b > 0, x^{0} = y, v^{0} = Wx^{0}, p_{1}^{0} = 0 $;

1.2:对于$k = 0, 1, \cdots, n, $计算$v^{k+1} = {\rm Shrink}(Wx^{k} - \frac{p_{1}^{k}}{b}, \frac{\mu}{b}), $

$ x^{k+1} = (A^{T}A + bI)^{-1}(A^{T}y + bW^{T}v^{k+1} + W^{T}p_{1}^{k} ), p_{1}^{k + 1} = p_{1}^{k} - \delta_{b} (Wx^{k+1} - v^{k+1}); $

1.3:输出$ x^{n} $;

2.1:输入新的初始值$ y, A, \lambda > 0, a > 0, x^{n} = y, u^{0} = \nabla x^{n}, p_{2}^{0} = 0 $;

2.2:对于$k = n+1, \cdots, m, $计算$u^{k+1}=(u_{1}^{k+1}, u_{2}^{k+1}) $,

$ u_{1}^{k+1} = {\rm Shrink}(\nabla_{1} x^{k} - \frac{p_{2, 1}^{k}}{a}, \frac{\lambda}{a}), u_{2}^{k+1} = {\rm Shrink}(\nabla_{2} x^{k} - \frac{p_{2, 2}^{k}}{a}, \frac{\lambda}{a}), \\ x^{k+1} = (A^{T}A + a\nabla^{T}\nabla)^{-1}(A^{T}y + b\nabla^{T}u^{k+1} + \nabla^{T}p_{2}^{k} ), p_{2}^{k + 1} = p_{2}^{k} - \delta_{a} (\nabla x^{k+1} - u^{k+1}); $

2.3:达到停机准则后输出$ \hat{x} = x^{k + 1} $.

选定的停机准则为

$ \frac{\left\| x^{k+1} - x^{k} \right\|_{2}}{\max({\left\| x^{k} \right\|_{2}}, 1)} < \epsilon, $

其中$ \epsilon > 0 $是一个给定的值.

本文选用峰值信噪比$ PSNR = 10 \cdot \log{\frac{255^{2}}{MSE}} $指标来衡量图像的去噪效果, 其中$MSE$表示去噪前后图像的均方误差. $PSNR$值越大, 则说明去噪效果越好.选取B样条小波作为低通滤波器, 即$h_{0}=[\frac{1}{4}, \frac{1}{2}, \frac{1}{4}]$.两个高通滤波器分别为$h_{1}=\frac{\sqrt{2}}{4}[1, 0, -1]$, $h_{2}=\frac{1}{4}[-1, 2, -1]$, 三个滤波器组成实验所用框架算子.在去模糊实验中, $\delta_{a}=\delta_{b}=0.05$, $\mu=0.0002$, $\lambda=0.005$, $a=b=0.5$; 作为对照, TV部分$\lambda=0.0003$, $\delta_{a}=0.5$和$a=1$; 框架对照部分$\mu=0.0003$, $\delta_{b}=0.5$和$b=1$.实验结果在图 1中显示, 峰值信噪比见表 1.

在去噪声实验中, $\delta_{a}=0.5$, $\delta_{b}=1$, $\mu=5\times e^{-3}$, $\lambda=5\times e^{-3}$, $a=b=1$; 作为对照, TV部分$\lambda=0.01$, $\delta_{a}=1$和$a=1$; 框架对照部分$\mu=0.01$, $\delta_{b}=1.5$和$b=1$.实验结果在图 2中显示, 峰值信噪比见表 2.

在带噪声图像恢复后的背景和草地的细节中, 能够很明显地看到本文方法的优势在于保持了更多的细节信息.从PSNR的数据中, 可以看到从模糊图像和带噪声图像的恢复中, 本文的方法得到了更好的恢复结果.

针对图像从模糊和噪声中恢复的问题, 本文结合了小波框架和TV模型, 提出了新的$L_{1}$正则化模型, 并且应用ADMM算法来结算这个新的模型.数值实验的结果显示本文的方法在图像恢复过程中保留原始图像的边缘信息和细节上具有优势.于是本文的方法尤其适合恢复带有很多纹理细节的图像信息.在实验得到的PSNR数值结果佐证了这个判断.而这种将两种模型结合的思想可以应用在更多的场景之中.在今后的研究中可以应用同样的思想, 将不同的算子组合来恢复具有不同的几何信息和不同噪声的图像.