网站地图
-
在获取图像的过程中, 获取的图像不可避免地受到噪声的干扰。图像的噪声可以理解为影响人的视觉器官或传感器对所接收图像源信息进行理解或分析的各种因素。它对图像的输入、采集和处理的各个环节, 以及输出结果的全过程都有影响。因此消除噪声已成为图像处理中极重要的问题。数字图像的噪声主要来源于图像的获取和传输过程。常见的噪声有:白噪声、冲击噪声、量化噪声、椒盐噪声。本文中考虑白噪声的一个特例——高斯噪声的消除。
最近几年涌现了许多数字图像处理的方法, 半光滑牛顿法(semi-smooth Newton method, SSN)[1-6]也是其中一种, 即采用光滑函数来逼近不可微目标函数。随着对半光滑问题研究的不断深入, 该方法得到迅速的发展, 是最优化领域中极其活跃的研究方向之一。利用半光滑牛顿法求解最优化问题时, 由于线性方程组系数矩阵可能是病态的, 常采用Tikhovon正规化方法[7]进行求解。目前采用的SSN在求解非对称线性方程组的计算时间长、收敛速度慢。文中将继续探究怎样利用半光滑牛顿法快速求解非对称线性方程组。
本文中首先给出图像去噪模型, 其次介绍最优化问题的最优化条件和求解该条件的结合广义最小残差法(generalized minimum residual method, GMRES)[7-9]的半光滑牛顿法, 并分析其局部超线性收敛性, 然后利用数值结果验证该算法的性能并与其它算法作比较, 最后得出结论。
-
在数学上, 去噪模型表示为:
$ \mathit{\boldsymbol{y}} = \mathit{\boldsymbol{x}} + \mathit{\boldsymbol{\eta }} $
(1) 式中, x是原始数据或无噪声数据, η是服从高斯分布的随机噪声, y是观察到的数据。由于观察到的图像是带有噪声的, 所以要解决的问题是怎样由观察到的数据y估计原始数据x。
图像去噪的本质是一个病态的反问题[10]。所谓病态即不适定, 反问题的不适定性对相应问题的数值解法将产生本质的影响。正则化理论[11-12]常用来处理这种病态问题。(1)式的噪声为高斯白噪声, 原始数据可以通过求解下面的最小二乘正则化问题来估计:
$ \mathop {\min }\limits_\mathit{\boldsymbol{x}} \frac{1}{2}\left\| {\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{x}}} \right\|_2^2 + \lambda {\left\| {\mathit{\boldsymbol{Ax}}} \right\|_1} $
(2) 式中, 第1项为数据拟合项, 反映了观察图像对原始图像的逼近程度; 第2项为正则项, λ为正则化参量, 用来平衡近似解的逼近程度和平滑性; A为差分矩阵; ‖·‖2表示2-范数, ‖·‖1表示1-范数。
由于1-范数的不可微性, 求解最小化问题(2)式的稳定性和有效性将会变得困难。
定义共轭函数。设函数f:Rn→R, 定义函数f*Rn→R。则有f*(y)=$\mathop {\sup }\limits_{\mathit{\boldsymbol{x}} \in {\rm{dom}}f} $(yTx-f(x)), 称此函数为函数f的共轭函数[13]。其中R为实数空间, Rn表示n维实数空间, dom f为函数f的定义域, T表示转置。
利用共轭函数定义[13-14], 知道‖Ax‖1可表示为$\mathop {\max }\limits_{{{\left\| z \right\|}_\infty } \le 1} $zTAx。从而(2)式转化为:
$ \mathop {\max }\limits_{{{\left\| \mathit{\boldsymbol{z}} \right\|}_\infty } \le 1} \mathop {\min }\limits_\mathit{\boldsymbol{x}} \frac{1}{2}\left\| {\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{x}}} \right\|_2^2 + \lambda {\mathit{\boldsymbol{z}}^{\rm{T}}}\mathit{\boldsymbol{Ax}} $
(3) 式中, z是一个对偶变量, ‖·‖∞表示无穷范数。
由于(2)式中1-范数的不可微性和高斯白噪声的随机性, 想要稳定的解出x是困难的, 而(3)式对变量x是可导的。与最优化问题(2)式相比, 最大最小化问题[15](3)式可以转换成一个带约束的最小二乘问题。(3)式关于变量x的解为:
$ \mathit{\boldsymbol{x}} = \mathit{\boldsymbol{y}} - \lambda {\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{z}} $
(4) 将(4)式代回(3)式, 则对偶变量z是如下优化问题的解:
$ \mathop {\min }\limits_{{{\left\| \mathit{\boldsymbol{z}} \right\|}_\infty } \le 1} \frac{1}{2}\left\| {{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{z}} - \frac{\mathit{\boldsymbol{y}}}{\lambda }} \right\|_2^2 $
(5) 反过来, 只需求出z, 然后代入(4)式即可求解出x。传统上, 基于投影的梯度下降法、半隐性梯度法等用来求解(5)式。但是这些方法需要的迭代次数多, 收敛速度慢。下一节中将运用具有超线性收敛速度的半光滑牛顿法来解(5)式。
-
在这一节中, 首先导出最优化问题(5)式的最优化条件, 然后介绍最优化条件的半光滑性, 最后利用结合GMRES的半光滑牛顿法对其求解以及分析该算法的优越性——局部超线性收敛性。
-
引入Lagrange乘子, 将有约束问题的(5)式转化为无约束问题:
$ \mathop {\min }\limits_{\mathit{\boldsymbol{z}},\mathit{\boldsymbol{\alpha }}} \frac{1}{2}\left\| {{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{z}} - \frac{\mathit{\boldsymbol{y}}}{\lambda }} \right\|_2^2 + \sum\limits_{i = 0}^n {{\mathit{\boldsymbol{\alpha }}_i}\left( {1 - \mathit{\boldsymbol{z}}_i^2} \right)} $
(6) 由此可得最优化条件:
$ \left\{ \begin{array}{l} \mathit{\boldsymbol{A}}\left( {{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{z}} - \frac{\mathit{\boldsymbol{y}}}{\lambda }} \right) + 2\mathit{\boldsymbol{\alpha z}} = 0\\ {\mathit{\boldsymbol{\alpha }}_i}\left( {1 - \mathit{\boldsymbol{z}}_i^2} \right) = 0 \end{array} \right. $
(7) 式中, , z=(z1, z2, …, zn)T, i, j=1, 2, 3, …, n。
由于αi≥0, 1-zi2≥0, (7)式的第2个式子也称为互补问题[14]。
-
定义F-B函数。对于互补问题, 利用由:
$ \varphi \left( {a,b} \right) = \sqrt {{a^2} + {b^2}} - a - b $
(8) 定义的Fischer-Burmeister(F-B)函数[13]φ:R2→R来进行定式。其中a≥0, b≥0。
利用函数φ定义ψi:Rn→R(i=1, 2, 3, …, n)函数:ψi(x)=φ(xi, Fi(x)), 并令ψi(x)=(ψ1(x), …, ψn(x))T, 则关于向量值函数ψ:Rn→Rn的方程ψ(x)=0与互补问题等价。但ψ在满足xi=Fi(x)=0的x处不可微, 称为非光滑方程组。这也是大家熟知的非线性互补问题(nonlinear complementarity problem, NCP), 本文中采用半光滑牛顿法求解这类非光滑方程组。下面引入次微分及牛顿导数的定义。
定义次微分。记φ可微的点的全体构成的集合为Sφ, 定义:
$ \begin{array}{*{20}{c}} {{\partial _{\rm{B}}}\varphi \left( \mathit{\boldsymbol{x}} \right) = \left\{ {\mathop {\lim }\limits_{k \to \infty } \nabla \varphi \left( {{\mathit{\boldsymbol{x}}^{\left( k \right)}}} \right)\left| {\mathop {\lim }\limits_{k \to \infty } {\mathit{\boldsymbol{x}}^{\left( k \right)}} = \mathit{\boldsymbol{x}},} \right.} \right.}\\ {\left. {{\mathit{\boldsymbol{x}}^{\left( k \right)}} \subseteq {S_\varphi }} \right\}} \end{array} $
(9) 该式称为φ在x处的次微分或Bouligand次微分(简记为B次微分)[16]。
定义牛顿/弱可微。假定f:Rn→Rm, x∈int dom f, 函数f在x处牛顿可微的定义是, 存在矩阵Df(z)∈Rm×n满足:
$ \mathop {\lim }\limits_{\mathit{\boldsymbol{z}} \in {\rm{dom}}\;f,\mathit{\boldsymbol{z}} \ne \mathit{\boldsymbol{x}},z \to x} \frac{{{{\left\| {f\left( \mathit{\boldsymbol{z}} \right) - f\left( \mathit{\boldsymbol{x}} \right) - \mathit{\boldsymbol{D}}f\left( \mathit{\boldsymbol{z}} \right)\left( {\mathit{\boldsymbol{z}} - \mathit{\boldsymbol{x}}} \right)} \right\|}_2}}}{{{{\left\| {\mathit{\boldsymbol{z}} - \mathit{\boldsymbol{x}}} \right\|}_2}}} = 0 $
(10) 则称f:Rn→Rm在点x是牛顿(或弱)可微的, 且将Df(z)称为f在x处的牛顿导数[7](或弱函数)。
众所周知, 局部Lipschitz连续函数均是牛顿可微的(Rademacher定理[16])。它在x处的Clarke次微分与B次微分之间的关系式∂φ(x)=co(∂Bφ(x))成立。其中co(∂Bφ(x))表示∂Bφ(x)的凸包。
由于F-B函数是局部Lipschitz连续函数, 则其Clarke次微分可由下式给出:
$ \begin{array}{*{20}{c}} {\partial \varphi \left( {a,b} \right) = }\\ {\left\{ \begin{array}{l} {\left[ {1 - \frac{a}{{\sqrt {{a^2} + {b^2}} }},1 - \frac{b}{{\sqrt {{a^2} + {b^2}} }}} \right]^{\rm{T}}},\left( {\left( {a,b} \right) \ne \left( {0,0} \right)} \right)\\ {\left( {1 - \xi ,1 - \eta } \right)^{\rm{T}}}\left| {{\xi ^2} + {\eta ^2} \le 1,\left( {\left( {a,b} \right) = \left( {0,0} \right)} \right)} \right. \end{array} \right.} \end{array} $
(11) 由F-B函数的定义, 令F-B函数为:
$ \begin{array}{*{20}{c}} {f\left( {{\mathit{\boldsymbol{\alpha }}_i},{\mathit{\boldsymbol{z}}_i}} \right) = \sqrt {\mathit{\boldsymbol{\alpha }}_i^2 + {{\left( {1 - \mathit{\boldsymbol{z}}_i^2} \right)}^2}} - }\\ {{\mathit{\boldsymbol{\alpha }}_i} - \left( {1 - \mathit{\boldsymbol{z}}_i^2} \right) = 0} \end{array} $
(12) 由于(12)式是牛顿可微的和局部Lipschitz连续的, 则Clarke次微分为:
$ \begin{array}{*{20}{c}} {\partial f\left( {{\mathit{\boldsymbol{\alpha }}_i},{\mathit{\boldsymbol{z}}_i}} \right) = }\\ {\left\{ \begin{array}{l} {\left[ {\frac{{ - 2{\mathit{\boldsymbol{z}}_i}\left( {1 - \mathit{\boldsymbol{z}}_i^2} \right)}}{{\sqrt {\mathit{\boldsymbol{\alpha }}_i^2 + {{\left( {1 - \mathit{\boldsymbol{z}}_i^2} \right)}^2}} }} + 2{\mathit{\boldsymbol{z}}_i},\frac{{{\mathit{\boldsymbol{\alpha }}_i}}}{{\sqrt {\mathit{\boldsymbol{\alpha }}_i^2 + {{\left( {1 - \mathit{\boldsymbol{z}}_i^2} \right)}^2}} }} - 1} \right]^{\rm{T}}},\\ \;\;\;\;\left( {\left( {{\mathit{\boldsymbol{\alpha }}_i},1 - {\mathit{\boldsymbol{z}}_i}} \right) \ne \left( {0,0} \right)} \right)\\ {\left( {2{\mathit{\boldsymbol{z}}_i} + \xi ,1 - \eta } \right)^{\rm{T}}}\left| {{\xi ^2} + {\eta ^2} \le 1} \right.,\\ \;\;\;\;\left( {\left( {{\mathit{\boldsymbol{\alpha }}_i},1 - {\mathit{\boldsymbol{z}}_i}} \right) = \left( {0,0} \right)} \right) \end{array} \right.} \end{array} $
(13) 根据F-B函数的定义可将最优化条件(7)式变形为:
$ \left\{ \begin{array}{l} \mathit{\boldsymbol{A}}\left( {{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{z}} - \mathit{\boldsymbol{y}}{\lambda ^{ - 1}}} \right) + 2\mathit{\boldsymbol{\alpha }}z = 0\\ \sqrt {\mathit{\boldsymbol{\alpha }}_i^2 + {{\left( {1 - \mathit{\boldsymbol{z}}_i^2} \right)}^2}} - {\mathit{\boldsymbol{\alpha }}_i} - \left( {1 - {\mathit{\boldsymbol{z}}_i}^2} \right) = 0 \end{array} \right. $
(14) 运用Taylor展式得:
$ \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{P}}&\mathit{\boldsymbol{Q}}\\ \mathit{\boldsymbol{H}}&\mathit{\boldsymbol{K}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\delta _\mathit{\boldsymbol{z}}}}\\ {{\delta _\mathit{\boldsymbol{\alpha }}}} \end{array}} \right] = - \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{C}}_1}}\\ {{\mathit{\boldsymbol{C}}_2}} \end{array}} \right] $
(15) 式中, P=AAT+2α, Q=2z, C1=A(ATz-yλ-1)+2αz, ${\mathit{\boldsymbol{C}}_2} = \sqrt {\mathit{\boldsymbol{\alpha }}_i^2 + {{\left( {1 - z_i^2} \right)}^2}} - {\mathit{\boldsymbol{\alpha }}_i} - \left( {z_i^2} \right)$; δz和δα表示迭代步长。
当(αi, 1-zi)≠(0, 0)时, H是以H1=-2zi(1-zi2)/$\sqrt {\mathit{\boldsymbol{\alpha }}_i^2 + {{\left( {1 - z_i^2} \right)}^2}} + 2{z_i}$为对角线元素的对角矩阵, K是以K1=$\sqrt {\mathit{\boldsymbol{\alpha }}_i^2 + {{\left( {1 - z_i^2} \right)}^2}} + 2{z_i}$为对角线元素的对角矩阵。当(αi, 1-zi)=(0, 0)时, H是以H2=2zi+ξ为对角线元素的对角矩阵, K是以K2=1-η为对角线元素的对角矩阵, 且ξ2+η2≤1。
-
传统在解线性方程组(15)式时[3], 利用了预处理共轭梯度法(preconditioned conjugate gradients method, PCG)[17-18]。这种方法是处理大型病态方程组的有效数值迭代算法, 对矩阵进行预处理使条件数大大减小, 从而减小其病态性。但由于PCG方法是一种数值迭代方法, 主要求解对称的线性方程组, 而(15)式的系数矩阵为非对称矩阵, 如果应用PCG, 那么必须先转化为正规方程组, 增加了计算量, 影响SSN最终的计算时间。所以接下来, 介绍求解该非对称线性方程组的另一种方法——GMRES, 并利用数值实验对比整体效果。
广义最小残差法于1986年由SAAD和SCHULTZ提出[7], 最吸引人的地方在于算法的高效性和稳定性。每一步只需要进行少量(1个或2个)矩阵向量的乘积, 与解同样非线性方程组的其它方法相比具有存贮量少的特点, 并且能充分利用矩阵的稀疏结构, 计算过程不再需要单独对矩阵进行存储, 从而大大节约了内存, 且易于进行计算。
GMRES算法的具体步骤如下:
(1) 选初始迭代向量x(0), 计算残量r(0)=b-Gx(0)和初始正交向量v(1)=r(0)/‖r(0)‖2。
(2) 对j=1, 2, …; i=1, 2, …, j。
$ \left\{ \begin{array}{l} {\mathit{\boldsymbol{h}}^{\left( {i,j} \right)}} = \left( {\mathit{\boldsymbol{G}}{\mathit{\boldsymbol{v}}^{\left( j \right)}},{\mathit{\boldsymbol{v}}^{\left( i \right)}}} \right),{{\mathit{\boldsymbol{\hat v}}}^{\left( {j + 1} \right)}} = \mathit{\boldsymbol{G}}{\mathit{\boldsymbol{v}}^{\left( j \right)}} - \sum\limits_{i = 1}^j {{\mathit{\boldsymbol{h}}^{\left( {i,j} \right)}}{\mathit{\boldsymbol{v}}^{\left( i \right)}}} \\ {\mathit{\boldsymbol{h}}^{\left( {j + 1,i} \right)}} = {\left\| {{{\mathit{\boldsymbol{\hat v}}}^{\left( {j + 1} \right)}}} \right\|_2},{\mathit{\boldsymbol{v}}^{\left( {j + 1} \right)}} = {{\mathit{\boldsymbol{\hat v}}}^{\left( {j + 1} \right)}}/{\mathit{\boldsymbol{h}}^{\left( {j + 1,i} \right)}} \end{array} \right. $
(16) (3) 形成近似解x(k)=x(0)+u(k), 其中u(k)为最小二乘问题的解:
$ \begin{array}{*{20}{c}} {{{\left\| {{\mathit{\boldsymbol{r}}^{\left( k \right)}}} \right\|}_2} = {{\left\| {\mathit{\boldsymbol{b}} - \mathit{\boldsymbol{G}}\left( {{\mathit{\boldsymbol{x}}^{\left( 0 \right)}} + {\mathit{\boldsymbol{u}}^{\left( k \right)}}} \right)} \right\|}_2} = }\\ {{{\left\| {{\mathit{\boldsymbol{r}}^{\left( 0 \right)}} - \mathit{\boldsymbol{G}}{\mathit{\boldsymbol{u}}^{\left( k \right)}}} \right\|}_2} = \min \left\{ {{{\left\| {{\mathit{\boldsymbol{r}}^{\left( 0 \right)}} - \mathit{\boldsymbol{Gu}}} \right\|}_2}} \right\}} \end{array} $
(17) 利用此方法对(15)式求解δz和δα, 然后代入迭代公式:
$ \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{z}}^{\left( {k + 1} \right)}}}\\ {{\mathit{\boldsymbol{\alpha }}^{\left( {k + 1} \right)}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{z}}^{\left( k \right)}}}\\ {{\mathit{\boldsymbol{\alpha }}^{\left( k \right)}}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{\delta _\mathit{\boldsymbol{z}}}}\\ {{\delta _\mathit{\boldsymbol{\alpha }}}} \end{array}} \right] $
(18) 解得最优值z, 最后代回(4)式可求得最优解x。
结合GMRES的半光滑牛顿法见下。其中, 函数为:x=fSSN(y, A, λ)。输入为:y, A, λ。
(1) 初始化z(0), α(0)。
(2) for k=1:c, 其中c为最大迭代次数。
(3) 用半光滑牛顿法求解(14)式的第2个式子:
$ \partial f\left( {{\mathit{\boldsymbol{\alpha }}_i},{\mathit{\boldsymbol{z}}_i}} \right) = \left\{ \begin{array}{l} {\left( {{\mathit{\boldsymbol{H}}_1},{\mathit{\boldsymbol{K}}_1}} \right)^{\rm{T}}},\left( {\left( {{\mathit{\boldsymbol{\alpha }}_i},1 - {\mathit{\boldsymbol{z}}_i}} \right) \ne \left( {0,0} \right)} \right)\\ {\left( {{\mathit{\boldsymbol{H}}_2},{\mathit{\boldsymbol{K}}_2}} \right)^{\rm{T}}},\left( {\left( {{\mathit{\boldsymbol{\alpha }}_i},1 - {\mathit{\boldsymbol{z}}_i}} \right) = \left( {0,0} \right)} \right) \end{array} \right. $
(19) 式中, 当(αi, 1-zi)≠(0, 0)时, H1是以-2zi(1-zi2)/$\sqrt {\mathit{\boldsymbol{\alpha }}_i^2 + {{\left( {1 - z_i^2} \right)}^2}} + 2{z_i}$为对角线元素的对角矩阵, K1是以${\mathit{\boldsymbol{\alpha }}_i}/\sqrt {\mathit{\boldsymbol{\alpha }}_i^2 + {{\left( {1 - z_i^2} \right)}^2}} - 1$为对角线元素的对角矩阵; 当(αi, 1-zi)=(0, 0)时, H1是以2zi+ξ为对角线元素的对角矩阵, K1是以1-1-η为对角线元素的对角矩阵, 且ξ2+η2≤1。
(4) 通过求解(18)式更新z和α。
(5) 当迭代次数达到c或相对误差l < 10-4, 结束。
(6) 代回(4)式:x=y-λATz。
(7) 输出x。
-
在本小节中将阐述结合GMRES的半光滑牛顿法的优越性——局部超线性收敛性。这里, 局部意味着只有当选择的初始点y(0)=(z(0), α(0))充分接近真实解y†=(z†, α†)时, 结论才成立, 其中变量z是变量x的对偶变量。
由于方程组(14)式的第1式是光滑函数, 所以现在只考虑(14)式的第2式:
$ \begin{array}{*{20}{c}} {F\left( \mathit{\boldsymbol{y}} \right) = F\left( {\mathit{\boldsymbol{z}},\mathit{\boldsymbol{\alpha }}} \right) = \sqrt {{\mathit{\boldsymbol{\alpha }}^2} + {{\left( {1 - {\mathit{\boldsymbol{z}}^2}} \right)}^2}} - }\\ {\mathit{\boldsymbol{\alpha }} - \left( {1 - {\mathit{\boldsymbol{z}}^2}} \right) = 0} \end{array} $
(20) 式中, F为Lipschitz连续且半光滑函数。求解(13)式的广义牛顿法迭代公式如下:
$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{y}}^{\left( {k + 1} \right)}} = {\mathit{\boldsymbol{y}}^{\left( k \right)}} - \mathit{\boldsymbol{V}}_k^{ - 1}F\left( {{\mathit{\boldsymbol{y}}^{\left( k \right)}}} \right),}\\ {\left( {{\mathit{\boldsymbol{V}}_k} \in {\partial _{\rm{B}}}F\left( {{\mathit{\boldsymbol{y}}^{\left( k \right)}}} \right)} \right)} \end{array} $
(21) 设y†为(20)式的解, 若任意的V∈∂F(y†)为非奇异的, 参考文献[19]中证明了由(21)式给出的迭代公式是适定的, 并具有超线性收敛性质, 即:
$ \left\| {{\mathit{\boldsymbol{y}}^{\left( {k + 1} \right)}} - {\mathit{\boldsymbol{y}}^\dagger }} \right\| = o\left( {\left\| {{\mathit{\boldsymbol{y}}^{\left( k \right)}} - {\mathit{\boldsymbol{y}}^\dagger }} \right\|} \right) $
(22) 通过第2节中的讨论知道, 结合GMRES的半光滑牛顿法在求解非奇异非对称线性方程组问题时, 有着收敛速度快、占用存储空间小的优点, 且算法是局部超线性收敛的。但该算法只适用于高斯白噪声的去噪问题, 其它噪声类型的情况需要另做推广。
-
本文中讨论运用结合GMRES的半光滑牛顿法对含有随机高斯白噪声图像去噪, 并与结合PCG的半光滑牛顿法、ADMM[20-22]方法进行比较。虽然ADMM算法出现比较早, 但是在最近几年才应用到图像处理中且效果比较好。作者做了大量对比实验, 程序编写在MATLAB R2007b中进行。采用定量指标:峰值信噪比(peak signal-to-noise ratio, PSNR)和均方误差(mean squared error, MSE)对图像质量进行评价。其定义分别如下:
$ \left\{ \begin{array}{l} {R_{{\rm{PSNR}}}} = 10\lg \frac{{{{255}^2}}}{{\sum\limits_i {\sum\limits_j {\frac{{{{\left( {{\mathit{\boldsymbol{J}}_{ij}} - {\mathit{\boldsymbol{I}}_{ij}}} \right)}^2}}}{{MN}}} } }}\\ {R_{{\rm{MSE}}}} = \frac{1}{{MN}}\sum\limits_i {\sum\limits_j {{{\left( {{\mathit{\boldsymbol{J}}_{ij}} - {\mathit{\boldsymbol{I}}_{ij}}} \right)}^2}} } \end{array} \right. $
(23) 式中, M和N为图像尺寸大小, I和J分别是原始无噪图像和恢复后的图像, Iij和Jij表示图像像素。RPSNR值越大或RMSE值越小, 恢复效果越好。在理论和实践的指导下, 一般使用原始图像x的连续迭代值之间的差‖x(k+1)-x(k)‖22来度量算法的收敛性。另外引入相对误差l作为算法的停止条件。其定义如下:
$ l = \frac{{\left\| {{\mathit{\boldsymbol{x}}^{\left( {k + 1} \right)}} - {\mathit{\boldsymbol{x}}^{\left( k \right)}}} \right\|_2^2}}{{\left\| {{\mathit{\boldsymbol{x}}^{\left( {k + 1} \right)}}} \right\|_2^2}} $
(24) 设定当l≤10-4或者迭代次数达到最大迭代次数c时停止迭代。
-
先给出差分矩阵A和信号x:
$ \left\{ \begin{array}{l} \mathit{\boldsymbol{A}} = {\left( {\begin{array}{*{20}{c}} { - 1}&1&{}&{}&{}\\ {}&{ - 1}&1&{}&{}\\ {}&{}& \ddots&\ddots &{}\\ {}&{}&{}&{ - 1}&1 \end{array}} \right)_{\left( {n - 1} \right) \times n}}\\ {x_i} = \left\{ \begin{array}{l} 5{x_i} - 0.5,\left( {0.1 < {x_i} \le 0.3} \right)\\ 1.20,\left( {0.4 < {x_i} \le 0.42} \right)\\ 1.20,\left( {0.45 < {x_i} \le 0.48} \right)\\ 1.0,\left( {0.85 < {x_i} \le 0.95} \right)\\ 0,\left( {{\rm{otherwise}}} \right) \end{array} \right. \end{array} \right. $
(25) 式中, x=(x1, x2, …, xn), i=1, 2, …n。去噪模型为(1)式, 这里噪声数据是服从n(0, 1)正态分布的随机高斯白噪声。选择不同的正则化参量λ, 在解非对称线性方程组时比较分别结合GMRES, PCG的半光滑牛顿法和ADMM这3种算法, 观察定量指标RMSE的变化情况。
图 1a中给出了高斯噪声的方差δ=8的信号图形的真实数据与噪声数据。图 1b、图 1c和图 1d表示在高斯噪声的方差δ=8、正则化参量λ=50的情况下, 分别结合GMRES, PCG的半光滑牛顿法和ADMM这3种算法的对比去噪图形。图 1e表示均方误差RMSE随迭代时间t的变化曲线。图 1a、图 1b、图 1c和图 1d的横坐标v均表示离散点的坐标, 无单位。图 1a纵坐标w表示观察到的离散点所对应的值, 图 1b、图 1c和图 1d的纵坐标w分别表示用不同方法恢复离散点所对应的值, 均无单位。表 1中给出当正则化参量λ=50时, 随着参量β的变化, 3种算法的迭代次数k、迭代时间t和均方误差RMSE的变化数据。
Figure 1. Experiment results of 1-D signal(δ=8, λ=50)
表 1中, k1, k2和k3分别表示结合GMRES, PCG的SSN和ADMM算法恢复信号的迭代次数; t1, t2和t3分别表示结合GMRES, PCG的SSN和ADMM算法恢复信号的迭代时间; m1, m2和m3分别表示结合GMRES, PCG的SSN和ADMM算法恢复信号的RMSE。
Table 1. Experiment parameters when λ=50
β k1 t1 M1 k2 t2 M2 k3 t3 M3 2.2 29 2.1485 0.0662×105 15 19.9191 0.0992×105 1000 5.2583 0.1394×105 2.3 49 2.3740 0.0501×105 17 21.2814 0.0702×105 1000 5.4383 0.1225×105 2.4 86 3.1606 0.0484×105 20 25.9716 0.0581×105 1000 4.3907 0.1349×105 2.5 70 4.1595 0.0347×105 15 17.9194 0.0495×105 1000 4.1653 0.1121×105 2.6 44 1.5466 0.0450×105 15 16.9880 0.0579×105 1000 4.0729 0.1467×105 2.7 54 3.3145 0.0470×105 15 18.7179 0.0663×105 1000 4.0212 0.1565×105 2.8 45 3.0259 0.0437×105 17 19.6219 0.0774×105 1000 4.2353 0.1620×105 2.9 62 2.2909 0.0472×105 16 19.7418 0.0861×105 1000 4.1886 0.1748×105 3.0 45 3.0073 0.0509×105 24 28.2398 0.0743×105 1000 4.9195 0.1854×105 3.1 56 2.1859 0.0591×105 22 28.3739 0.0820×105 1000 5.0466 0.1949×105 3.2 80 2.7877 0.0421×105 18 22.8552 0.0565×105 1000 3.9656 0.1915×105 3.3 44 3.2971 0.0514×105 15 17.4596 0.0793×105 1000 4.0851 0.1971×105 3.4 84 5.0000 0.0511×105 15 19.5113 0.0866×105 1000 5.6617 0.2348×105 3.5 73 2.7443 0.0432×105 18 23.2715 0.0594×105 1000 4.4383 0.2249×105 由表 1看出:ADMM算法运行1000次没有达到收敛标准, 随着参量β不同, 运行时间可能不同。虽然结合GMRES的SSN迭代次数不是最少, 但迭代时间比结合PCG的SSN和ADMM的少很多, 而且它的RMSE比结合PCG的SSN和ADMM的RMSE小很多, 效果很好。综合图 1和表 1发现, 结合GMRES的SSN和结合PCG的SSN都具有局部超线性收敛性, ADMM只有线性收敛性。
-
选取图像处理中常用的Cameraman图像作为测试目标, 并对其加入标准方差δ=12的高斯白噪声。Cameraman图像中不但含有像素跳跃区域(相机支架)、图像渐变区域(天空), 而且还含有图像震荡区域(草坪)。实验中取正则参量λ=8, 图 2中给出了结合GMRES的半光滑牛顿法和结合PCG的半光滑牛顿法去噪结果。其中图 2a为原始图像; 图 2b是含噪图像; 图 2c是结合GMRES的半光滑牛顿法的去噪结果; 图 2d是结合PCG的半光滑牛顿法的去噪结果; 图 2e和图 2f分别是定量指标RPSNR和RMSE随时间t的变化曲线图。
Figure 2. Experimental results of 2-D image(δ=12, λ=8)
从图 2可以看出, 无论从RPSNR, RMSE还是时间t方面作比较, 结合GMRES的半光滑牛顿法均优于结合PCG的半光滑牛顿法。
-
利用半光滑牛顿法处理全变差噪声消除问题。结合GMRES的半光滑牛顿法具有局部超线性收敛性的优越性, 不仅能够有效的提高恢复信号和图像的质量, 而且与结合PCG的半光滑牛顿法、ADMM相比时, 收敛速度快且稳定。
全变差噪声消除问题的半光滑牛顿法
Semi-smooth Newton method for total variation noise removal
-
摘要: 为了达到全变差噪声消除的图像去噪目的，将去噪问题转换为优化问题。采用了结合广义最小残差法的半光滑牛顿法来解决相关优化问题，求解非对称线性方程组，进行了理论分析和实验验证，取得了将该方法与其它方法应用于1维信号、2维图像去噪实验的大量可行数据。结果表明，结合广义最小残差法的半光滑牛顿法的收敛速度比结合预处理共轭梯度法的半光滑牛顿法和交替方向乘子法更快，而且能够有效地消除噪声。Abstract: In order to remove the noise of image based on total variation, the denoising problem was converted to optimization problem. Semi-smooth Newton method incorporated by generalized minimum residual method was used to solve the associated optimization problem and non-symmetric linear equations. After theoretical analysis and experimental verification, a great deal of feasible data of removal noise experiment for 1-D signal and 2-D image were obtained by different methods. The results show that semi-smooth Newton method incorporated by generalized minimum residual method converges faster than that incorporated by preconditioned conjugate gradients method and alternating direction method of multipliers algorithm. The proposed method can remove the noise of image effectively.
-
Figure 1. Experiment results of 1-D signal(δ=8, λ=50)
a—real data and noise data b—SSN incorporated by GMRES c—SSN incorporated by PCG d—alternating direction method of multipliers e—RMSE variation curve vs. iteration time t
Figure 2. Experimental results of 2-D image(δ=12, λ=8)
a—original image b—noise image c—recover image of SSN incorporated by GMRES d—recover image of SSN incorporated by PCG e—RPSNR variation curve vs. iteration time t f—RMSE variation curve vs. iteration time t
Table 1. Experiment parameters when λ=50
β k1 t1 M1 k2 t2 M2 k3 t3 M3 2.2 29 2.1485 0.0662×105 15 19.9191 0.0992×105 1000 5.2583 0.1394×105 2.3 49 2.3740 0.0501×105 17 21.2814 0.0702×105 1000 5.4383 0.1225×105 2.4 86 3.1606 0.0484×105 20 25.9716 0.0581×105 1000 4.3907 0.1349×105 2.5 70 4.1595 0.0347×105 15 17.9194 0.0495×105 1000 4.1653 0.1121×105 2.6 44 1.5466 0.0450×105 15 16.9880 0.0579×105 1000 4.0729 0.1467×105 2.7 54 3.3145 0.0470×105 15 18.7179 0.0663×105 1000 4.0212 0.1565×105 2.8 45 3.0259 0.0437×105 17 19.6219 0.0774×105 1000 4.2353 0.1620×105 2.9 62 2.2909 0.0472×105 16 19.7418 0.0861×105 1000 4.1886 0.1748×105 3.0 45 3.0073 0.0509×105 24 28.2398 0.0743×105 1000 4.9195 0.1854×105 3.1 56 2.1859 0.0591×105 22 28.3739 0.0820×105 1000 5.0466 0.1949×105 3.2 80 2.7877 0.0421×105 18 22.8552 0.0565×105 1000 3.9656 0.1915×105 3.3 44 3.2971 0.0514×105 15 17.4596 0.0793×105 1000 4.0851 0.1971×105 3.4 84 5.0000 0.0511×105 15 19.5113 0.0866×105 1000 5.6617 0.2348×105 3.5 73 2.7443 0.0432×105 18 23.2715 0.0594×105 1000 4.4383 0.2249×105 
下载: 导出CSV
-
[1] MICHAEL K, QI L Q, YANG Y F, et al. On semi-smooth Newton methods for total variation minimization[J]. Manufactured in the Netherlands, 2007, 27(3):265-276. [2] PANG Zh F, LÜ J Ch. The modified semi-smooth Newton algorithm based on the ROF mode[J]. Pure Mathematics, 2011, 1(4):26-29(in Chinese). [3] YANG M. Semi-smooth Newton method for solving nonlinear complementarity problems[D]. Xi'an: Xi'an Electronic and Science University, 2009: 9-20(in Chinese). [4] CHEN X J, NASHED Z H, QI L Q. Smoothing methods and semismooth methods for nondifferentiable operator equations[J]. SIAM Journal Numerical Analysis, 2000, 38(4):1200-1216. [5] CHEN M Y. Application of semi-smooth Newton method for solving elastic wave parameters inversion in porous media[D]. Dalian: Dalian Maritime University, 2013: 1-23(in Chinese). [6] GRIESSE R, LOREN D A. A semismooth Newton method for Tikhonov functionals with sparsity constraints[J]. Inverse Problem, 2008, 24(3):1-22. SAAD Y, SCHULTZ M. GMRES:a generalized minimal residual algorithm for solving nonsymmetric linear systems.SIAM Journal of Scientific and Statistical Computing, 1986, 7(3):856-869. [7] GUAN P Y, LI C G, JING H F. A weighted GMRES method based on polynomial preconditioning generalized minimal residual method[J]. Journal of Engineering Mathematics, 2014, 31(5):697-706(in Chinese). [8] MEI D. GMRES algorithm based on solution space decomposition and its application in image processing[J]. Computer and Digital Engineering, 2009, 37(12):139-143(in Chinese). [9] WANG Q, DAI H. Regularization GMERR method for solving discrete ill-posed problems[J]. Mathematica Numerica Sinica, 2013, 35(2):195-204(in Chinese). [10] WANG M J, PAN Z Y. A regularization method for solving nonlinear ill-posed problems[J]. Journal of Xi'an Polytechnic University, 2013, 27(5):675-679(in Chinese). [11] FUHRY M. A new tikhovon regularization method[J].Numerical Algorithms, 2012, 59(3):433-445. doi: 10.1007/s11075-011-9498-x [12] CHAMBOLLE A, POCK T. A first-order primal-dual algorithm for convex problems with applications to imaging[J]. Journal of Mathematical Imaging and Vision, 2011, 40(1):120-145. doi: 10.1007/s10851-010-0251-1 [13] WANG S N, XU J, HUANG X L. Convex optimization[M]. Beijing:Tsinghua University Press, 2013:85-90(in Chinese). [14] CAI X J, HAN D, XU L L. An improved first-order primal-dual algorithm with a new correction step[J]. Springer Science+Business Media, 2013, 57(4):1419-1428. [15] STEPHEN B, LIEVEN V. Convex optimization[M]. Cambridge, UK:World Book Publishing Company, 2004:90-92. [16] ZHANG T B, XIE J, ZHU J J. The solution of the problem of the solution of the conjugate gradient method and its application in GPS[J]. Surveying and Mapping Engineering, 2010, 19(4):60-63(in Chinese). [17] SUN Y J, SUN Q. Large complex linear equations with two conjugate gradient method[J]. Computer Engineering and Appliction, 2007, 43(36):19-20(in Chinese). [18] LI W R. Semismooth newton algorithm for NCP with a non-monotone line search[D]. Tianjin: Tianjin University, 2008: 16-20(in Chinese). [19] DANIELE G, LANURE BF, GILLES A. ADMM algorithm for demosaicking deblurring denoising[J/OL].(2013-10-18)[2016-01-13].https://hal.archives-ouvertes.fr/hal-00874610/file/admcolorpaper-5.pdf. [20] JOAO F C M, JOAO M F X, PEDRO M Q A. ADMM for consensus on colored networks[J]. Decision and Control, 2012, 8267(1):5116-5121. [21] YASAR K A, ORHAN A. ADMM based mainlobe power constrained phase-only sidelobe supperssion[C]//2014 IEEE 22nd Signal Processing and Communications Applications Conference (SIU). New York, USA: IEEE, 2014: 610-614. [22] YASAR K A, ORHAN A. ADMM based mainlobe power constrained phase-only sidelobe supperssion[C]//2014 IEEE 22nd Signal Processing and Communications Applications Conference (SIU).New York, USA: IEEE, 2014: 610-614. -
点击查看大图
图(2) / 表(1)
计量
- 文章访问数: 4909
- HTML全文浏览量: 3016
- PDF下载量: 416
- 被引次数: 0
