牛顿迫近迭代算法在图像恢复中的应用

李旭超 刘 燕 李玉叶

1(赤峰学院计算机与信息工程学院 内蒙古 赤峰 024000) 2(赤峰学院数学与统计学院 内蒙古 赤峰 024000)

牛顿迫近迭代算法在图像恢复中的应用

李旭超1刘 燕1李玉叶2

1(赤峰学院计算机与信息工程学院 内蒙古 赤峰 024000)2(赤峰学院数学与统计学院 内蒙古 赤峰 024000)

由光滑与非光滑函数构成的混合目标函数，传统的一阶优化算法，由于光滑函数一阶逼近的欠准确性和搜索步长的限制，很难获得目标函数的高精度解。针对此问题，提出二阶牛顿迫近算子分裂迭代算法。对光滑函数进行泰勒展开，获得目标函数的二阶转化模型，将转化模型分解为牛顿迭代子问题和迫近迭代子问题；给出牛顿迭代子问题的搜索方向和最优搜索步长；对算法的收敛特性进行分析。利用被系统和噪声退化的图像进行恢复实验，结果表明，该方法比现有方法峰值信噪比最高提高约2 dB，结构相似测度提高约3%。

非光滑特性 牛顿迫近算法 迭代收敛 图像恢复

0 引 言

在工业探伤、X射线断层扫描成像及雷达成像等反问题中，受成像系统非适定(ill-posed)特性和噪声的影响，很难获得高精度有效解[1]。为解决此问题，常用光滑函数和非光滑函数的组合描述成像系统[2]，建立目标优化函数，表达式为：

(1)

式中：inf表示下确界，x*表示目标函数的最优值，s(x)是光滑函数，体现成像系统的统计特性，使采样与成像系统相吻合；n(x)是非光滑函数，体现解的奇异特性，常用某一函数空间来描述，如一阶、二阶有界变差TV(total variation)函数空间[3]，变指数函数空间等[4]。但式(1)整体的非光滑特性，使得优化算法设计十分困难。针对此问题，国内外学者主要从两个方面进行研究。一是对式(1)进行整体处理。如对n(x)进行光滑化，使得式(1)整体可微分，利用经典优化理论，设计牛顿迭代算法、拟牛顿迭代算法和投影迭代算法等[5]。但整体处理获得的海森矩阵规模往往比较大，若不具有特殊结构，造成迭代算法收敛较慢，而且需要进行光滑化和泰勒展开两次逼近，很难获得高精度有效解。二是分别利用光滑函数、非光滑函数的特性，对式(1)进行分裂处理[6]。文献[7]提出分裂两步软阈值迭代算法，但算法的收敛速度较慢。针对此问题，文献[8]利用光滑函数的梯度形成一阶迭代子问题，利用非光滑函数形成软阈值迭代子问题，二者交替形成迭代阈值收缩算法，但算法收敛速度是次线性的。为提高算法的运算速度，文献[9]采用Nesterov多步加速策略[10]，提出收敛速度是二阶的快速迭代软阈值算法，但算法是有条件地收敛，依赖于系统矩阵的特征值，并且容易造成逼近解产生阶梯效应。

鉴于此，本文主要做以下工作。首先回顾软阈值迭代算法，在此基础上，对光滑函数和非光滑函数进行分裂，形成牛顿迭代子问题和迫近迭代子问题，两个子问题交替迭代形成牛顿迫近迭代算法，并给出步长更新准则；然后分析算法的收敛特性；最后利用被系统和噪声退化的图像进行恢复实验。在恢复性能定量评价上，根据峰值信噪比(PSNR)和结构相似测度(SSIM)[11]，验证算法的有效性。

1 快速迭代软阈值算法

1.1 目标函数的一阶逼近模型

利用Lipschitz常数，对光滑函数s(x)进行一阶逼近，表达式为：

(2)

QL(x,yk)=s(x,yk)+n(x)

(3)

经整理，获得式(3)最小化表达式为：

(4)

1.2 一阶转化模型进行分裂

利用拟合项的光滑特性，将光滑函数转化为最小二乘的形式，但由于模型整体的非光滑特性，利用经典牛顿迭代算法无法直接求解。但随着非光滑优化理论的发展[12]，特别是迫近算子的提出，为非光滑函数模型的求解注入新的活力。将转化模型式(4)表示成迫近算子的标准型[13]，表达式为：

(5)

利用光滑函数的梯度，xk形成一阶不动点迭代子问题，表达式为：

(6)

利用非光滑函数的次微分特性[13]，式(5)形成迫近迭代子问题，表达式为：

(7)

1.3 快速迭代软阈值算法

利用光滑函数的梯度设计一阶不动点迭代子问题，利用非光滑函数的次微分特性设计迫近迭代子问题，二者交替形成快速软阈值迭代算法。算法的具体步骤为：

初始化，计算L，y1=x0，t1=1，步长μ，最大迭代次数N；

循环迭代：

(8)

(9)

(10)

(11)

当达到最大迭代次数或目标函数连续两次迭代之差小于给定值，循环终止，输出恢复图像。将快速迭代软阈值算法应用于图像恢复，取得较快的收敛速度，但图像恢复质量不太理想，容易产生阶梯效应，且算法收敛步长受成像系统特征值的制约。

2 牛顿迫近迭代算法

2.1 目标函数的二阶逼近模型

从式(2)可知，用光滑函数一阶梯度的Lipschitz常数作为校正项的系数，对光滑函数s(x)进行逼近。但随着迭代次数的增加，光滑函数的梯度随之发生变化，导致Lipschitz常数L也随之发生变化。也就是说，L是迭代次数的函数。因此，用常数L无法准确逼近光滑函数，导致无法获得高精度逼近解。为解决此问题，对光滑函数进行泰勒展开，用随迭代次数变化而变化的海森矩阵取代Lipschitz常数，表达式为：

(12)

式中：s(x，yk)表示对s(x)在yk处进行泰勒展开，o(x-yk)表示高阶逼近项，Hk表示海森矩阵，k表示迭代次数。式(1)的转化模型表达式为：

QHk(x,yk)=s(x，yk)+n(x)

(13)

忽略高阶项和常数项，经整理，获得式(13)最小化表达式：

(14)

2.2 二阶转化模型进行分裂

利用迫近算子，将式(14)表示成标准型，表达式为：

(15)

式中：xk形成二阶牛顿迭代子问题，表达式为：

(16)

而式(15)形成迫近算子迭代子问题，表达式为：

(17)

为获得模型有效解，下面对式(16)进行算法设计。

2.3 牛顿迭代算法搜索方向和步长的确定

(18)

为使步长可变，算法具有自适应性，式(16)可以表示为：

xk+1=yk+αkdk

(19)

(20)

而式(20)是迫近迭代子问题、牛顿迭代子问题的目标解校正差，由式(18)可知，Hkdk的本质是梯度，利用Karush-Kuhn-Tucke(KKT)条件[12]，则目标函数获得最优解的条件表达式为：

(21)

(22)

而式(22)恰好是目标函数式(1)取得最优值的充分必要条件。由式(14)可知，在计算牛顿迭代子问题时，需要计算海森矩阵的逆矩阵，若其规模较大，计算其逆矩阵比较耗时，而且其特征值的幅值可能较小，具有奇异特性，为便于计算且避免矩阵的奇异特性，需对海森矩阵进行逼近。若对称正定矩阵Σk是光滑函数s(x)的海森矩阵Hk的逼近矩阵，那么Σk满足Secant方程的三个条件[2]，即：

(23)

由BFGS原理[14]，海森矩阵的逼近矩阵Σk的更新表达式为：

(24)

由式(16)和式(19)，牛顿迭代子问题可以表示为：

(25)

为使式(25)中的αk具有自适应性，步长更新表达式为：

(26)

2.4 牛顿迭代算法的搜索步长是最优的

定理1牛顿迫近迭代算法的步长式(26)是最优的。

证明：

(27)

由于目标函数是凸函数，所以：

(28)

(29)

O((xk+1-xk)Hk(xk+1-xk)(xk+1-xk))+αk·

(30)

因为，目标函数的解唯一，由KKT条件式(21)，则有：

(31)

(32)

由次微分的定义，则有：

(33)

将式(33)代入式(32)，则有：

(34)

对式(30)进一步放大，将式(34)代入式(30)，则有：

Φ(xk+1)≤Φ(xk)-[αkβk-O(αkλk)]

(35)

令O(akλk)=-akλk-ln(1-akλk)并代入式(35)，为使能量泛函最大限度地下降，式(35)关于步长ak取得极值，对其求偏导数，则有式(26)。证毕。

2.5 牛顿迫近迭代算法

初始化，最大迭代次数N，迭代残差ε，正则项参数γ。

循环迭代：

(2) 计算牛顿迭代子问题式(25)获得xk。

(4) 用式(26)更新步长，用式(19)更新解。

当达到最大迭代次数或目标函数连续两次迭代之差小于ε，循环终止，输出恢复图像，否则，返回(1)。

3 牛顿迫近迭代算法的收敛性

证明：假设式(1)存在最优解，由式(35)可知，当步长取值为式(26)时，与步长相关函数取到极值的表达式为：

(36)

因此，式(35)可以表示为：

Φ(xk+1)≤Φ(xk)-Mk

(37)

(38)

Φ(xk+1)≤Φ(x1)-τM

(39)

(40)

证毕。

4 被系统和噪声降质图像的恢复模型

4.1 模糊图像恢复形式化表示及算法流程图

假设理想图像为x，受系统A和泊松噪声而降质，经采样获得的模糊图像为η，为拟合成像过程，式(1)中的s(x)用Kullback-Leibler光滑函数[15]来描述，表达式为：

s(x)=Ax-ηlnAx

(41)

式中：A也称为点扩散函数[16]。

对于式(1)中的非光滑项，为体现图像的奇异特性，用TV函数空间的半范数来描述，表达式为：

n(x)=γ|x|TV

(42)

式中：γ为正则项参数，|·|表示半范数。

将式(41)、式(42)代入式(1)，模糊图像恢复的形式化表达式为：

(43)

由牛顿迫近迭代算法，则模糊图像恢复算法流程图，如图1所示。

图1 模糊图像恢复算法流程图

4.2 图像恢复实验方法

对于式(41)中的A，模拟成像系统受大气扰动使图像模糊，表达式为：

(44)

式中：φ为成像系统的孔径，φ为大气扰动的相位，F-1为逆傅里叶变换，PSF(·,·)表示点扩散函数，图2为0.1≤φ≤0.3，φ=0.6时产生的点扩散函数，z1、z2表示横、纵坐标。

图2 点扩散函数(200×200)

为验证本文算法的有效性，采用快速迭代软阈值算法(算法1)、Matlab系统中的deconvlucy函数实现快速交替迭代期望最大值算法(算法2)与本文算法进行对比实验，根据PSNR和SSIM的数值定量评价不同算法的恢复性能[17]。实验选用200×200的原始图像，灰度级为0～255，如图3所示。实验中，正则参数设置为γ=2.7×10-5，迭代残差ε=10-8，最大迭代次数N=103。实验硬件配置：Genuine Intel(R)CPU T2300@ 1.66 GHz，1.66 GB内存；软件配置：Windows 8操作系统，Matlab R2011a。

图3 原始图像

4.3 实验结果及分析

从视觉效果来看，图4(a)为降质的Barbara图像，可以看出，图像非常模糊，细节信息丢失严重。而用算法1恢复的视觉效果较差，图像产生明显的阶梯效应，如图4(b)。用算法2恢复图像的视觉效果好于算法1，但图像的细节信息恢复模糊，如图4(c)。图4(d)为用本文算法进行恢复，图像的细节信息恢复好于其他算法。

图4 不同算法恢复Barbara图像性能对比

图5(a)为模糊的斑马图像，斑马条纹比较模糊，无法分辨出“草地”的细节。图5(b)为用算法1进行恢复，产生非常明显的阶梯效应。图5(c)为用算法2恢复，恢复效果明显好于算法1，但“草地”细节比较模糊。图5(d)为用本文算法进行恢复，斑马的条纹和“草地”细节信息恢复比较理想。

图5 不同算法恢复斑马图像性能对比

图6(a)为模糊的鹦鹉图像，“羽毛”细节信息比较模糊。图6(b)为用算法1进行恢复，产生大块的阶梯效应，“羽毛”细节信息几乎全部丢失。用算法2恢复的效果明显好于算法1，如图6(c)所示。图6(d)为用本文算法进行恢复，“羽毛”细节信息恢复比较理想，但在平稳区域产生条纹，这是由于平稳与非平稳区域过渡处梯度幅值变化较大，造成海森矩阵逼近不准确所致。

图6 不同算法恢复鹦鹉图像性能对比

图7(a)为降质的Haifa图像。图7(b)用算法1恢复，产生明显的阶梯效应，细节信息严重被抹杀。图7(c)为用算法2恢复，视觉效果好于算法1。图7(d)为用本文算法进行恢复，图像恢复效果比较理想。

图7 不同算法恢复Haifa图像性能对比

从定量评价指标来看，对于Barbara图像，斑马图像和Haifa图像，用本文算法获得的PSNR和SSIM值高于其他两种方法，如表1、表2和表4所示，说明本文算法恢复效果最好。但对于鹦鹉图像，用算法2获得的SSIM数值高于本文算法，但本文算法获得的PSNR数值高于算法2，如表3所示，这也与恢复图像的视觉效果相吻合。

表1 不同算法恢复Barbara性能对比

表2 不同算法恢复斑马性能对比

表3 不同算法恢复鹦鹉性能对比

表4 不同算法恢复Haifa性能对比

5 结 语

从定性方面来看，对优化函数中光滑部分和非光滑部分进行分裂，通过降低海森矩阵规模，设计二阶牛顿迫近迭代算法，这种算法具有二阶收敛特性。并利用目标函数光滑和非光滑部分的特性，推导出牛顿步的搜索方向和最优步长更新准则，并证明算法的收敛特性。从定量结果来看，该算法与一阶迭代算法和同类分裂迭代算法进行图像恢复质量对比，结果表明，该算法恢复图像的视觉效果优于其他算法，定量评价指标PSNR、SSIM的数值高于一阶迭代算法和同类分裂迭代算法。

[1] 韩波,李莉.非线性不适定问题的求解方法及其应用[M].北京:科学出版社,2011.

[2] Vogel C R.Computational Methods for Inverse Problems[M].Philadelphia,Pennsylvania,USA:Society for Industrial and Applied Mathematics,2002.

[3] Fan Qibin,Jiang Dandan,Jiao Yuling.A multi-parameter regularization model for image restoration[J].Signal Processing,2015,114(9):131-142.

[4] 李旭超.能量泛函正则化模型在图像恢复中的应用[M].北京:电子工业出版社,2014.

[5] Landi G,Piccolomini E L.NPTool:a MATLAB software for nonnegative image restoration with Newton projection methods[J].Numerical Algorithms,2014,62(3):487-504.

[6] Deng L J,Guo H,Huang T Z.A fast image recovery algorithm based on splitting deblurring and denoising[J].Journal of Computational & Applied Mathematics,2015,287(C):88-97.

[7] Bioucas-Dias J M,Figueiredo M A T.A new TwIST:two-step iterative shrinkage/thresholding algorithms for image restoration[J].IEEE Transactions on image processing,2007,16(12):2992-3004.

[8] Beck A,Teboulle M.A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems[J].SIAM Journal on Imaging Sciences,2009,2(1):183-202.

[9] Beck A,Teboulle M.A fast dual proximal gradient algorithm for convex minimization and applications[J].Operations Research Letters,2014,42(1):1-6.

[10] Becker S,Bobin J,Candès E J.NESTA:A Fast and Accurate First-order Method for Sparse Recovery[J].SIAM Journal on Imaging Sciences,2009,4(1):1-39.

[11] Wang H,Ho A T S,Li S.A novel image restoration scheme based on structured side information and its application to image watermarking[J].Signal Processing Image Communication,2014,29(7):773-787.

[12] Bertsekas D P.凸优化理论[M].北京:清华大学出版社,2011.

[13] Villa S,Salzo S,Baldassarre L,et al.Accelerated and Inexact Forward-Backward Algorithms[J].SIAM Journal on Optimization,2013,23(23):1607-1633.

[14] Narushima Y,Yabe H,Ford J A.A Three-Term Conjugate Gradient Method with Sufficient Descent Property for Unconstrained Optimization[J].SIAM Journal on Optimization,2011,21(1):212-230.

[15] Harmany Z T,Marcia R F,Willett R M.This is SPIRAL-TAP:Sparse Poisson Intensity Reconstruction ALgorithms-theory and practice[J].IEEE Transactions on Image Processing A Publication of the IEEE Signal Processing Society,2012,21(3):1084-1096.

[16] Cho H M,Cho H S,Kim K S,et al.Experimental study on the application of a compressed-sensing-based deblurring method in x-ray nondestructive testing and its image performance[J].NDT &E International,2015,75(1):1-7.

[17] 李小利,杨晓梅.基于RPCA视频去噪算法的自适应优化方法[J].计算机应用与软件,2016,33(9):215-220.

APPLICATIONOFNEWTONPROXIMALITERATIVEALGORITHMFORIMAGERESTORATION

Li Xuchao1Liu Yan1Li Yuye2

1(CollegeofComputerandInformationEngineering,ChifengUniversity,Chifeng024000,InnerMongolia,China) >2(CollegeofMathematicsandStatistics,ChifengUniversity,Chifeng024000,InnerMongolia,China)

The mixture object function is composed of smooth and non-smooth function. The traditional first-order optimization algorithm is limited by the first-order approximation of the smooth function and the search step. And it is difficult to obtain a high-precision solution of the objective function. Therefore, we propose a second order Newton proximal operator splitting iterative algorithm. Firstly, Taylor expansion of the smoothing function was used to obtain the two order transformation model of the objective function. The transformed model was decomposed into Newton iterative subproblem and proximal iterative subproblem. Then, the search direction and the optimization search step length of Newton iterative sub-problem were given. Finally, the convergence property was analyzed. Taking advantage of image blurred by system and noise for restoration, we perform the recovery experiments. The results show the PSNR (peak signal to noise ratio) of the proposed method is about 2 dB higher than other methods, and the SSIM (structural similarity index measure) is improved by about 3%.

Non-smooth property Newton proximal algorithm Iterative convergence Image restoration

2017-04-01。国家自然科学基金项目(11402039)；2016年度内蒙古自治区科技厅自然科学基金项目(2016MS0 602)；2016年度内蒙古自治区高等学校科学研究项目(NJZY16254)。李旭超，教授，主研领域：图像恢复。刘燕，教授。李玉叶，副教授。

TP391

A

10.3969/j.issn.1000-386x.2017.11.038