牛頓迫近迭代算法在圖像恢復(fù)中的應(yīng)用

李旭超 劉 燕 李玉葉

1(赤峰學(xué)院計(jì)算機(jī)與信息工程學(xué)院 內(nèi)蒙古 赤峰 024000) 2(赤峰學(xué)院數(shù)學(xué)與統(tǒng)計(jì)學(xué)院 內(nèi)蒙古 赤峰 024000)

牛頓迫近迭代算法在圖像恢復(fù)中的應(yīng)用

李旭超1劉 燕1李玉葉2

1(赤峰學(xué)院計(jì)算機(jī)與信息工程學(xué)院 內(nèi)蒙古 赤峰 024000)2(赤峰學(xué)院數(shù)學(xué)與統(tǒng)計(jì)學(xué)院 內(nèi)蒙古 赤峰 024000)

由光滑與非光滑函數(shù)構(gòu)成的混合目標(biāo)函數(shù)，傳統(tǒng)的一階優(yōu)化算法，由于光滑函數(shù)一階逼近的欠準(zhǔn)確性和搜索步長的限制，很難獲得目標(biāo)函數(shù)的高精度解。針對(duì)此問題，提出二階牛頓迫近算子分裂迭代算法。對(duì)光滑函數(shù)進(jìn)行泰勒展開，獲得目標(biāo)函數(shù)的二階轉(zhuǎn)化模型，將轉(zhuǎn)化模型分解為牛頓迭代子問題和迫近迭代子問題；給出牛頓迭代子問題的搜索方向和最優(yōu)搜索步長；對(duì)算法的收斂特性進(jìn)行分析。利用被系統(tǒng)和噪聲退化的圖像進(jìn)行恢復(fù)實(shí)驗(yàn)，結(jié)果表明，該方法比現(xiàn)有方法峰值信噪比最高提高約2 dB，結(jié)構(gòu)相似測(cè)度提高約3%。

非光滑特性 牛頓迫近算法 迭代收斂 圖像恢復(fù)

0 引 言

在工業(yè)探傷、X射線斷層掃描成像及雷達(dá)成像等反問題中，受成像系統(tǒng)非適定(ill-posed)特性和噪聲的影響，很難獲得高精度有效解[1]。為解決此問題，常用光滑函數(shù)和非光滑函數(shù)的組合描述成像系統(tǒng)[2]，建立目標(biāo)優(yōu)化函數(shù)，表達(dá)式為：

(1)

式中：inf表示下確界，x*表示目標(biāo)函數(shù)的最優(yōu)值，s(x)是光滑函數(shù)，體現(xiàn)成像系統(tǒng)的統(tǒng)計(jì)特性，使采樣與成像系統(tǒng)相吻合；n(x)是非光滑函數(shù)，體現(xiàn)解的奇異特性，常用某一函數(shù)空間來描述，如一階、二階有界變差TV(total variation)函數(shù)空間[3]，變指數(shù)函數(shù)空間等[4]。但式(1)整體的非光滑特性，使得優(yōu)化算法設(shè)計(jì)十分困難。針對(duì)此問題，國內(nèi)外學(xué)者主要從兩個(gè)方面進(jìn)行研究。一是對(duì)式(1)進(jìn)行整體處理。如對(duì)n(x)進(jìn)行光滑化，使得式(1)整體可微分，利用經(jīng)典優(yōu)化理論，設(shè)計(jì)牛頓迭代算法、擬牛頓迭代算法和投影迭代算法等[5]。但整體處理獲得的海森矩陣規(guī)模往往比較大，若不具有特殊結(jié)構(gòu)，造成迭代算法收斂較慢，而且需要進(jìn)行光滑化和泰勒展開兩次逼近，很難獲得高精度有效解。二是分別利用光滑函數(shù)、非光滑函數(shù)的特性，對(duì)式(1)進(jìn)行分裂處理[6]。文獻(xiàn)[7]提出分裂兩步軟閾值迭代算法，但算法的收斂速度較慢。針對(duì)此問題，文獻(xiàn)[8]利用光滑函數(shù)的梯度形成一階迭代子問題，利用非光滑函數(shù)形成軟閾值迭代子問題，二者交替形成迭代閾值收縮算法，但算法收斂速度是次線性的。為提高算法的運(yùn)算速度，文獻(xiàn)[9]采用Nesterov多步加速策略[10]，提出收斂速度是二階的快速迭代軟閾值算法，但算法是有條件地收斂，依賴于系統(tǒng)矩陣的特征值，并且容易造成逼近解產(chǎn)生階梯效應(yīng)。

鑒于此，本文主要做以下工作。首先回顧軟閾值迭代算法，在此基礎(chǔ)上，對(duì)光滑函數(shù)和非光滑函數(shù)進(jìn)行分裂，形成牛頓迭代子問題和迫近迭代子問題，兩個(gè)子問題交替迭代形成牛頓迫近迭代算法，并給出步長更新準(zhǔn)則；然后分析算法的收斂特性；最后利用被系統(tǒng)和噪聲退化的圖像進(jìn)行恢復(fù)實(shí)驗(yàn)。在恢復(fù)性能定量評(píng)價(jià)上，根據(jù)峰值信噪比(PSNR)和結(jié)構(gòu)相似測(cè)度(SSIM)[11]，驗(yàn)證算法的有效性。

1 快速迭代軟閾值算法

1.1 目標(biāo)函數(shù)的一階逼近模型

利用Lipschitz常數(shù)，對(duì)光滑函數(shù)s(x)進(jìn)行一階逼近，表達(dá)式為：

(2)

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

(3)

經(jīng)整理，獲得式(3)最小化表達(dá)式為：

(4)

1.2 一階轉(zhuǎn)化模型進(jìn)行分裂

利用擬合項(xiàng)的光滑特性，將光滑函數(shù)轉(zhuǎn)化為最小二乘的形式，但由于模型整體的非光滑特性，利用經(jīng)典牛頓迭代算法無法直接求解。但隨著非光滑優(yōu)化理論的發(fā)展[12]，特別是迫近算子的提出，為非光滑函數(shù)模型的求解注入新的活力。將轉(zhuǎn)化模型式(4)表示成迫近算子的標(biāo)準(zhǔn)型[13]，表達(dá)式為：

(5)

利用光滑函數(shù)的梯度，xk形成一階不動(dòng)點(diǎn)迭代子問題，表達(dá)式為：

(6)

利用非光滑函數(shù)的次微分特性[13]，式(5)形成迫近迭代子問題，表達(dá)式為：

(7)

1.3 快速迭代軟閾值算法

利用光滑函數(shù)的梯度設(shè)計(jì)一階不動(dòng)點(diǎn)迭代子問題，利用非光滑函數(shù)的次微分特性設(shè)計(jì)迫近迭代子問題，二者交替形成快速軟閾值迭代算法。算法的具體步驟為：

初始化，計(jì)算L，y1=x0，t1=1，步長μ，最大迭代次數(shù)N；

循環(huán)迭代：

(8)

(9)

(10)

(11)

當(dāng)達(dá)到最大迭代次數(shù)或目標(biāo)函數(shù)連續(xù)兩次迭代之差小于給定值，循環(huán)終止，輸出恢復(fù)圖像。將快速迭代軟閾值算法應(yīng)用于圖像恢復(fù)，取得較快的收斂速度，但圖像恢復(fù)質(zhì)量不太理想，容易產(chǎn)生階梯效應(yīng)，且算法收斂步長受成像系統(tǒng)特征值的制約。

2 牛頓迫近迭代算法

2.1 目標(biāo)函數(shù)的二階逼近模型

從式(2)可知，用光滑函數(shù)一階梯度的Lipschitz常數(shù)作為校正項(xiàng)的系數(shù)，對(duì)光滑函數(shù)s(x)進(jìn)行逼近。但隨著迭代次數(shù)的增加，光滑函數(shù)的梯度隨之發(fā)生變化，導(dǎo)致Lipschitz常數(shù)L也隨之發(fā)生變化。也就是說，L是迭代次數(shù)的函數(shù)。因此，用常數(shù)L無法準(zhǔn)確逼近光滑函數(shù)，導(dǎo)致無法獲得高精度逼近解。為解決此問題，對(duì)光滑函數(shù)進(jìn)行泰勒展開，用隨迭代次數(shù)變化而變化的海森矩陣取代Lipschitz常數(shù)，表達(dá)式為：

(12)

式中：s(x，yk)表示對(duì)s(x)在yk處進(jìn)行泰勒展開，o(x-yk)表示高階逼近項(xiàng)，Hk表示海森矩陣，k表示迭代次數(shù)。式(1)的轉(zhuǎn)化模型表達(dá)式為：

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

(13)

忽略高階項(xiàng)和常數(shù)項(xiàng)，經(jīng)整理，獲得式(13)最小化表達(dá)式：

(14)

2.2 二階轉(zhuǎn)化模型進(jìn)行分裂

利用迫近算子，將式(14)表示成標(biāo)準(zhǔn)型，表達(dá)式為：

(15)

式中：xk形成二階牛頓迭代子問題，表達(dá)式為：

(16)

而式(15)形成迫近算子迭代子問題，表達(dá)式為：

(17)

為獲得模型有效解，下面對(duì)式(16)進(jìn)行算法設(shè)計(jì)。

2.3 牛頓迭代算法搜索方向和步長的確定

(18)

為使步長可變，算法具有自適應(yīng)性，式(16)可以表示為：

xk+1=yk+αkdk

(19)

(20)

而式(20)是迫近迭代子問題、牛頓迭代子問題的目標(biāo)解校正差，由式(18)可知，Hkdk的本質(zhì)是梯度，利用Karush-Kuhn-Tucke(KKT)條件[12]，則目標(biāo)函數(shù)獲得最優(yōu)解的條件表達(dá)式為：

(21)

(22)

而式(22)恰好是目標(biāo)函數(shù)式(1)取得最優(yōu)值的充分必要條件。由式(14)可知，在計(jì)算牛頓迭代子問題時(shí)，需要計(jì)算海森矩陣的逆矩陣，若其規(guī)模較大，計(jì)算其逆矩陣比較耗時(shí)，而且其特征值的幅值可能較小，具有奇異特性，為便于計(jì)算且避免矩陣的奇異特性，需對(duì)海森矩陣進(jìn)行逼近。若對(duì)稱正定矩陣Σk是光滑函數(shù)s(x)的海森矩陣Hk的逼近矩陣，那么Σk滿足Secant方程的三個(gè)條件[2]，即：

(23)

由BFGS原理[14]，海森矩陣的逼近矩陣Σk的更新表達(dá)式為：

(24)

由式(16)和式(19)，牛頓迭代子問題可以表示為：

(25)

為使式(25)中的αk具有自適應(yīng)性，步長更新表達(dá)式為：

(26)

2.4 牛頓迭代算法的搜索步長是最優(yōu)的

定理1牛頓迫近迭代算法的步長式(26)是最優(yōu)的。

證明：

(27)

由于目標(biāo)函數(shù)是凸函數(shù)，所以：

(28)

(29)

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

(30)

因?yàn)椋繕?biāo)函數(shù)的解唯一，由KKT條件式(21)，則有：

(31)

(32)

由次微分的定義，則有：

(33)

將式(33)代入式(32)，則有：

(34)

對(duì)式(30)進(jìn)一步放大，將式(34)代入式(30)，則有：

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

(35)

令O(akλk)=-akλk-ln(1-akλk)并代入式(35)，為使能量泛函最大限度地下降，式(35)關(guān)于步長ak取得極值，對(duì)其求偏導(dǎo)數(shù)，則有式(26)。證畢。

2.5 牛頓迫近迭代算法

初始化，最大迭代次數(shù)N，迭代殘差ε，正則項(xiàng)參數(shù)γ。

循環(huán)迭代：

(2) 計(jì)算牛頓迭代子問題式(25)獲得xk。

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

當(dāng)達(dá)到最大迭代次數(shù)或目標(biāo)函數(shù)連續(xù)兩次迭代之差小于ε，循環(huán)終止，輸出恢復(fù)圖像，否則，返回(1)。

3 牛頓迫近迭代算法的收斂性

證明：假設(shè)式(1)存在最優(yōu)解，由式(35)可知，當(dāng)步長取值為式(26)時(shí)，與步長相關(guān)函數(shù)取到極值的表達(dá)式為：

(36)

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

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

(37)

(38)

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

(39)

(40)

證畢。

4 被系統(tǒng)和噪聲降質(zhì)圖像的恢復(fù)模型

4.1 模糊圖像恢復(fù)形式化表示及算法流程圖

假設(shè)理想圖像為x，受系統(tǒng)A和泊松噪聲而降質(zhì)，經(jīng)采樣獲得的模糊圖像為η，為擬合成像過程，式(1)中的s(x)用Kullback-Leibler光滑函數(shù)[15]來描述，表達(dá)式為：

s(x)=Ax-ηlnAx

(41)

式中：A也稱為點(diǎn)擴(kuò)散函數(shù)[16]。

對(duì)于式(1)中的非光滑項(xiàng)，為體現(xiàn)圖像的奇異特性，用TV函數(shù)空間的半范數(shù)來描述，表達(dá)式為：

n(x)=γ|x|TV

(42)

式中：γ為正則項(xiàng)參數(shù)，|·|表示半范數(shù)。

將式(41)、式(42)代入式(1)，模糊圖像恢復(fù)的形式化表達(dá)式為：

(43)

由牛頓迫近迭代算法，則模糊圖像恢復(fù)算法流程圖，如圖1所示。

圖1 模糊圖像恢復(fù)算法流程圖

4.2 圖像恢復(fù)實(shí)驗(yàn)方法

對(duì)于式(41)中的A，模擬成像系統(tǒng)受大氣擾動(dòng)使圖像模糊，表達(dá)式為：

(44)

式中：φ為成像系統(tǒng)的孔徑，φ為大氣擾動(dòng)的相位，F(xiàn)-1為逆傅里葉變換，PSF(·,·)表示點(diǎn)擴(kuò)散函數(shù)，圖2為0.1≤φ≤0.3，φ=0.6時(shí)產(chǎn)生的點(diǎn)擴(kuò)散函數(shù)，z1、z2表示橫、縱坐標(biāo)。

圖2 點(diǎn)擴(kuò)散函數(shù)(200×200)

為驗(yàn)證本文算法的有效性，采用快速迭代軟閾值算法(算法1)、Matlab系統(tǒng)中的deconvlucy函數(shù)實(shí)現(xiàn)快速交替迭代期望最大值算法(算法2)與本文算法進(jìn)行對(duì)比實(shí)驗(yàn)，根據(jù)PSNR和SSIM的數(shù)值定量評(píng)價(jià)不同算法的恢復(fù)性能[17]。實(shí)驗(yàn)選用200×200的原始圖像，灰度級(jí)為0～255，如圖3所示。實(shí)驗(yàn)中，正則參數(shù)設(shè)置為γ=2.7×10-5，迭代殘差ε=10-8，最大迭代次數(shù)N=103。實(shí)驗(yàn)硬件配置：Genuine Intel(R)CPU T2300@ 1.66 GHz，1.66 GB內(nèi)存；軟件配置：Windows 8操作系統(tǒng)，Matlab R2011a。

圖3 原始圖像

4.3 實(shí)驗(yàn)結(jié)果及分析

從視覺效果來看，圖4(a)為降質(zhì)的Barbara圖像，可以看出，圖像非常模糊，細(xì)節(jié)信息丟失嚴(yán)重。而用算法1恢復(fù)的視覺效果較差，圖像產(chǎn)生明顯的階梯效應(yīng)，如圖4(b)。用算法2恢復(fù)圖像的視覺效果好于算法1，但圖像的細(xì)節(jié)信息恢復(fù)模糊，如圖4(c)。圖4(d)為用本文算法進(jìn)行恢復(fù)，圖像的細(xì)節(jié)信息恢復(fù)好于其他算法。

圖4 不同算法恢復(fù)Barbara圖像性能對(duì)比

圖5(a)為模糊的斑馬圖像，斑馬條紋比較模糊，無法分辨出“草地”的細(xì)節(jié)。圖5(b)為用算法1進(jìn)行恢復(fù)，產(chǎn)生非常明顯的階梯效應(yīng)。圖5(c)為用算法2恢復(fù)，恢復(fù)效果明顯好于算法1，但“草地”細(xì)節(jié)比較模糊。圖5(d)為用本文算法進(jìn)行恢復(fù)，斑馬的條紋和“草地”細(xì)節(jié)信息恢復(fù)比較理想。

圖5 不同算法恢復(fù)斑馬圖像性能對(duì)比

圖6(a)為模糊的鸚鵡圖像，“羽毛”細(xì)節(jié)信息比較模糊。圖6(b)為用算法1進(jìn)行恢復(fù)，產(chǎn)生大塊的階梯效應(yīng)，“羽毛”細(xì)節(jié)信息幾乎全部丟失。用算法2恢復(fù)的效果明顯好于算法1，如圖6(c)所示。圖6(d)為用本文算法進(jìn)行恢復(fù)，“羽毛”細(xì)節(jié)信息恢復(fù)比較理想，但在平穩(wěn)區(qū)域產(chǎn)生條紋，這是由于平穩(wěn)與非平穩(wěn)區(qū)域過渡處梯度幅值變化較大，造成海森矩陣逼近不準(zhǔn)確所致。

圖6 不同算法恢復(fù)鸚鵡圖像性能對(duì)比

圖7(a)為降質(zhì)的Haifa圖像。圖7(b)用算法1恢復(fù)，產(chǎn)生明顯的階梯效應(yīng)，細(xì)節(jié)信息嚴(yán)重被抹殺。圖7(c)為用算法2恢復(fù)，視覺效果好于算法1。圖7(d)為用本文算法進(jìn)行恢復(fù)，圖像恢復(fù)效果比較理想。

圖7 不同算法恢復(fù)Haifa圖像性能對(duì)比

從定量評(píng)價(jià)指標(biāo)來看，對(duì)于Barbara圖像，斑馬圖像和Haifa圖像，用本文算法獲得的PSNR和SSIM值高于其他兩種方法，如表1、表2和表4所示，說明本文算法恢復(fù)效果最好。但對(duì)于鸚鵡圖像，用算法2獲得的SSIM數(shù)值高于本文算法，但本文算法獲得的PSNR數(shù)值高于算法2，如表3所示，這也與恢復(fù)圖像的視覺效果相吻合。

表1 不同算法恢復(fù)Barbara性能對(duì)比

表2 不同算法恢復(fù)斑馬性能對(duì)比

表3 不同算法恢復(fù)鸚鵡性能對(duì)比

表4 不同算法恢復(fù)Haifa性能對(duì)比

5 結(jié) 語

從定性方面來看，對(duì)優(yōu)化函數(shù)中光滑部分和非光滑部分進(jìn)行分裂，通過降低海森矩陣規(guī)模，設(shè)計(jì)二階牛頓迫近迭代算法，這種算法具有二階收斂特性。并利用目標(biāo)函數(shù)光滑和非光滑部分的特性，推導(dǎo)出牛頓步的搜索方向和最優(yōu)步長更新準(zhǔn)則，并證明算法的收斂特性。從定量結(jié)果來看，該算法與一階迭代算法和同類分裂迭代算法進(jìn)行圖像恢復(fù)質(zhì)量對(duì)比，結(jié)果表明，該算法恢復(fù)圖像的視覺效果優(yōu)于其他算法，定量評(píng)價(jià)指標(biāo)PSNR、SSIM的數(shù)值高于一階迭代算法和同類分裂迭代算法。

[1] 韓波,李莉.非線性不適定問題的求解方法及其應(yīng)用[M].北京:科學(xué)出版社,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] 李旭超.能量泛函正則化模型在圖像恢復(fù)中的應(yīng)用[M].北京:電子工業(yè)出版社,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.凸優(yōu)化理論[M].北京:清華大學(xué)出版社,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視頻去噪算法的自適應(yīng)優(yōu)化方法[J].計(jì)算機(jī)應(yīng)用與軟件,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。國家自然科學(xué)基金項(xiàng)目(11402039)；2016年度內(nèi)蒙古自治區(qū)科技廳自然科學(xué)基金項(xiàng)目(2016MS0 602)；2016年度內(nèi)蒙古自治區(qū)高等學(xué)?？茖W(xué)研究項(xiàng)目(NJZY16254)。李旭超，教授，主研領(lǐng)域：圖像恢復(fù)。劉燕，教授。李玉葉，副教授。

TP391

A

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