牛頓迫近迭代算法在圖像恢復中的應用

李旭超 劉 燕 李玉葉

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