基于Ginzburg-Landau泛函二值圖像復原的快速算法

孫寶宸，李 君，常慧賓

（天津師范大學數學科學學院，天津 300387）

在圖像的采集和傳輸等過程中，噪聲污染和信息缺失會導致圖像質量降低.圖像噪聲有多種類型[1]，如高斯噪聲、泊松噪聲、椒鹽噪聲等.噪聲的多樣性使得圖像去噪成為一項富有挑戰且被廣泛關注的研究課題，相關學者提出了多種經典的去噪方法[2-3].

1992年，Rudin等提出的一種以全變差作為正則項的凸ROF模型[2]，可以通過求解其歐拉-拉格朗日方程得到全局最優解，且解是唯一的，之后，一系列快速求解的一階優化算法被相繼提出[4].針對具體噪聲的統計性質，可以通過保留全變差正則項改進不同的保真項度量來提高復原圖像的質量.ROF模型具有良好的性質，但是全變差項的不可微性使算法設計變得困難，為解決這一問題，一系列光滑化技術被相繼提出.文獻[5]首次提出使用一類基于雙阱勢函數（doublewell potential function）的修正模型（簡稱DP模型）處理圖像銳化問題，之后該模型也被應用于圖像增強[6]、修復[7-8]和去噪[9]等問題.在文獻[10-11]中，基于非局部建模的DP模型首次用于圖像分割與數據分類問題.DP模型是逐點可微但非凸的，因此對模型求解具有雙重困難：一是直接求解歐拉-拉格朗日方程只能得到局部最優解而無法保證達到全局最優解；二是DP模型的一階變分是非線性的，很難設計有效的優化算法和梯度流迭代格式.對于全隱梯度流格式，文獻[12]證明了若時間步長滿足一定的約束條件，則可以保證能量穩定；對于顯隱格式，在相場模型的求解中也得到了類似的結論[13].

近年來，一類保能量衰減的顯隱迭代格式被應用于相場模型的計算中，如，乘子方法、IEQ格式和SAV格式[13].這類迭代格式通過引入輔助變量改寫目標能量泛函，將非線性梯度流方程轉化為多變量線性系統，但在實際計算中沒有增加需要求解的變量，并能保證改寫的等價能量泛函隨迭代進行是衰減的.此外，交替方向乘子算法在圖像處理的反問題求解中取得了巨大的成功，但目前未見其應用于求解DP模型.

本文首先引入基于Ginzburg-Landau泛函松弛的雙阱勢模型和2種經典算法，進而通過構造特殊的優化問題，建立了Lie格式與MBO格式的聯系，對MBO格式提出了一種從算子分裂優化技術角度的理解方式；然后通過引入輔助變量設計了一類新的算子分裂算法，其子優化問題具有閉形式解；最后通過二值圖像復原實驗評估各種算法的有效性.

1 雙阱勢模型和經典算法

1.1 雙阱勢模型

ROF模型[2-3]的連續形式定義為

其中：λ為模型參數；Du為u的分布導數；f為污染圖像；u為待復原的圖像；Ω為二維歐氏空間的一個有界閉集.全變差函數空間為

顯然E1ROF（u）關于u是下凸的（若無特殊說明，本文所指的凸性均為下凸性）.對于上述優化問題直接設計迭代格式求解歐拉-拉格朗日方程[4]，即式（1），獲得唯一的最優解.

其中：ε為正常數；W（u）為雙阱勢函數（如W（u）=0.25（u2-1）2，它有2個全局極小值）.文獻[16]證明了當ε→0+時，GL（u）會Γ收斂[3]到u的全變差.

基于上述結果將BV（Ω）松弛到H1（Ω）可建立雙阱勢模型（DP模型）如下

1.2 MBO格式（全隱格式）

文獻[17]提出了一種有效的算子分裂格式（記為MBO）.采用MBO格式進行變量更新需要2步迭代：第1步利用向后Euler差分格式求解反應擴散方程ut=Δu；第2步使用1/2水平集進行閾值化.由于DP模型中W（u）項含有2個駐點0和1，所以第2步迭代可表達為

將MBO應用于求解DP模型得到一種修正的MBO格

1.3 凸分裂格式

對于非凸優化問題，若通過能量恒等變換可將非凸的能量泛函轉化為凸能量與凹能量之和的形式，則凸分裂格式可視為一種簡單有效的數值求解算法，但它通常只有一階收斂速度.文獻[7]將DP模型分解為EDP（u）=E1（u）+E2（u），其中E1（u）=‖Δu‖2+μ‖u（u-1）‖2，為構造凸能量，在不改變目標能量的前提下取E（2u）=λ‖f-u‖2，對E（iu）（i=1、2）設計凸分裂格式.將E（iu）分解為E（iu）=E（i1u）-E（i2u），其中：

文獻[7]證明了上述凸分裂格式（convex splitting with 2 mixtures，CS2M）是無條件能量穩定的，即對任意的正整數n，有EDP（un+1）≤EDP（un）.

2 快速迭代求解方法

本節設計一些高效求解DP模型的迭代格式.首先將界面問題的新算法推廣到圖像復原問題的求解中，然后設計一類基于交替方向乘子方法（alternating direction method of multipliers，ADMM）的新算法.本文算法中均假設u滿足Neumann邊界條件，即0，其中：Γ為區域Ω的邊界，n為單位外法向量.

2.1 相場算法的推廣

利用梯度流方程求解能量泛函F的極小化問題，其一般框架為[13]

2.1.1 EQCN格式

文獻[21-22]對一類具有特殊結構的目標泛函設計了二階的顯隱差分格式，其核心是通過引進輔助變量將能量泛函的非平方項轉化為二次形式，進一步設計能量穩定的時間、空間全離散數值格式，本文考慮其時間半離散格式.定義輔助變量V=u（u-1），將DP模型進行變量分解，則式（3）轉化為

式（11）的梯度流AC方程組為

文獻[22]通過引入外插和Crank-Nicolson（CN）中心插商構建了迭代格式較為簡單且具有二階收斂速度的EQCN格式.記

利用式（14）進行變量更新，最后使用式（4）進行閾值化.EQCN格式可以保證能量穩定且具有二階收斂速度[23].

2.1.2 SAVCN格式

文獻[13]中提出了具有一階收斂速度的向前Euler格式、高階CN和BDF時間半離散格式，本文主要考慮SAVCN格式，其變量更新滿足

結合式（15）和式（16）對線性系統進行迭代求解，最后使用式（4）進行閾值化.SAVCN格式可以保證能量穩定且具有二階收斂速度[13].

2.2 交替方向乘子方法

交替方向乘子方法[24-26]可以處理非凸非光滑的復雜約束優化問題.考慮如下優化問題

上述框架可以自然地推廣到無窮維空間上.基于此，本文運用算子分裂技術，通過引入輔助變量將無約束優化問題轉化為約束優化問題，設計了2種迭代格式，其子優化問題具有閉形式解.

2.2.1 兩變量交替方向乘子算法

引入輔助變量v使得u=v，將W（u）進行算子分裂，原無約束優化問題可修改為約束優化問題：

通過交替優化變量和乘子可以得到求解DP問題的新算子分裂格式，其更新框架為

由于其目標泛函是凸可微的，因此可由一階最優性條件

同u子問題的求解過程一致，可得

乘子Λ依照式（18c）進行更新.滿足終止條件后使用式（4）對u進行閾值化.上述優化算法簡記為ADMM2V，其算法流程見表1.

表1 ADMM2V算法流程Tab.1 Algorithm flow of ADMM2V

2.2.2 三變量交替方向乘子算法

類比EQCN格式的解耦思想，本文設計了一種三變量交替方向乘子算法.引入2個輔助變量w和v將W（u）進行算子分裂，無約束優化問題可以轉化為約束優化問題：

其中：r1和r2為正的增廣拉格朗日系數，Λ1和Λ2為乘子.通過交替優化變量和乘子可以得到求解DP問題的新算子分裂格式，其更新框架為

對于乘子Λ1和Λ2按照式（21d）進行更新，滿足終止條件后使用式（4）對u進行閾值化.上述優化算法簡記為ADMM3V，其算法流程見表2.

表2 ADMM3V算法流程Tab.2 Algorithm flow of ADMM3V

3 數值實驗

由于更新格式（12）和（15）考慮的是Neumann邊界條件，因此本文對離散圖像添加鏡像邊界，依據文獻[27]，對連續的梯度算子離散得Δu（i，j）=（Δux（i，j），Δuy（i，j）），令Nx和Ny分別為離散圖像矩陣的行數和列數，其中：

其中：p=（p1，p2），p1=Δux，p2=Δuy.外迭代終止條件設為迭代1 000次或變量u的相對殘差小于10-4，即‖un+1-un‖/‖un‖＜10-4.對子問題的線性方程組使用共軛梯度法迭代求解，并設內迭代終止條件為迭代10次或相對誤差小于10-10.

3.1 圖像去噪

本文主要測試了3種加性噪聲（高斯噪聲、椒鹽噪聲和泊松噪聲）的復原問題.盡管針對不同種類的噪聲，可以設計不同的保真度量以獲得更好的去噪效果，但對于不可微的保真項會使基于梯度信息的SAVCN和EQCN格式失效，因此在圖像去噪實驗中統一采用L2模數作為數據保真項.采用信噪比（SNR）評價去噪效果，

其中：u為去噪后圖像，Img為無噪聲污染圖像.分別采用指紋、條形碼和二維碼圖像（Nx=Ny=256）進行高斯去噪、椒鹽去噪和泊松去噪實驗，計算mMBO、CS2M、EQCN、SAVCN、ADMM2V和ADMM3V等6種算法的SNR.基于2.2中對mMBO格式的解析，只討論其他5種算法的數值收斂性和能量衰減性.

對二值化的指紋圖像添加服從正態分布N（0，σ2）的高斯噪聲，取σ=0.2、0.5和0.8的噪聲圖像，見圖1，各算法復原圖像的SNR見表3，6種算法SNR的最大值用黑體標出，下同.

圖1 指紋真實圖像和高斯噪聲圖像Fig.1 Real image and Gaussian noisy images of fingerprint

由表3可見，除mMBO外，其他5種算法的SNR基本相同，均高于mMBO.圖2給出了變量u的相對誤差和總能量EDP（u）的變化曲線.由圖2（a）—（c）可見，5種算法的收斂速度沒有顯著差異.由圖2（d）—（f）可見，各算法都保證了能量衰減性，但能量衰減的速度不一致.當σ=0.2時，EQCN和SAVCN能量下降得最快；當σ=0.5時，各算法的能量降速基本一致，當σ=0.8時，ADMM能量下降得最快.

表3 高斯去噪實驗的SNRTab.3 SNR of image denoising experiments of Gaussian noisydB

圖2 高斯去噪實驗相對誤差與總能量的變化情況Fig.2 Changes of relative residuals and total energy in Gaussian denoising experiment

對二值化的條形碼圖像添加強度為ρ的椒鹽噪聲，取ρ=0.05、0.10和0.20的噪聲圖像，見圖3，各算法復原圖像的SNR見表4.由表4可見，除mMBO外，其他5種算法的SNR基本相同，均高于mMBO.圖4給出了變量u的相對誤差和總能量EDP（u）的變化曲線.由圖4（a）—（c）可見，5種算法的收斂速度沒有顯著差異.由圖4（d）—（f）可見，各算法都保證了能量衰減性，對于ρ=0.05、0.10和0.20，均為ADMM能量下降得最快.

圖3 條形碼真實圖像與椒鹽噪聲圖像Fig.3 Real image and salt and pepper noisy images of barcode

表4 椒鹽去噪實驗的SNRTab.4 SNR of image denoising experiments of salt and pepper noisy dB

圖4 椒鹽去噪實驗相對誤差與總能量的變化情況Fig.4 Changes of relative residuals and total energy in salt and pepper denoising experiment

對二值化的二維碼圖像逐點添加泊松噪聲，服從Poisson（ηf（i，j））/η分布，取η=4、8和16的噪聲圖像，見圖5，各算法復原圖像的SNR見表5.由表5可見，除mMBO外，其他5種算法的SNR均較高，ADMM2V的SNR最高.圖6給出了變量u的相對誤差和總能量EDP（u）的變化曲線.由圖6（a）—（c）可見，5種算法的收斂速度沒有顯著差異，但在優先保證EDP（u）充分下降的前提下，2種ADMM算法的收斂速度較快.由圖6（d）—（f）可見，各算法都保證了能量衰減性，對于η=4、8和16，均為ADMM能量下降得最快.

圖6 泊松去噪實驗相對誤差與總能量的變化情況Fig.6 Changes of relative residuals and total energy in Poisson denoising experiment

表5 泊松去噪實驗的SNRTab.5 SNR of image denoising experiments of Poisson noisy dB

圖5 二維碼真實圖像與泊松噪聲圖像Fig.5 Real image and Poisson noisy images of QR code

3.2 圖像修復

對于圖像修復問題需要將DP模型中的λ∈R+修改為

其中：Ωb為待修復部分；λ0為一個正常數.考慮簡單結構和復雜結構2類圖像的修復問題，以修復后圖像與真實圖像的均方誤差（MSE）評價修復效果，

在簡單結構圖像修復實驗中，對二值化的矩形拼接圖（Nx=Ny=128）設置不同比例的區域信息缺失，缺失比例分別設為30%、50%和70%，見圖7，各算法修復圖像的MSE見表6，6種算法MSE的最小值用黑體標出，下同.

圖7 矩形拼接圖的真實圖像與待修復圖像Fig.7 Real image of rectangular and images to be inpainting

由表6可見，SAVCN算法的修復效果最優.圖8給出了變量u的相對誤差和總能量EDP（u）的變化曲線.由圖8（a）—（c）可見，在迭代初期，2種ADMM算法收斂較快，而隨著迭代次數增多，SAVCN格式的收斂速度超過其他算法.由圖8（d）—（f）可見，各算法都保證了能量衰減性，但CS2M的總能量隨迭代次數增多先上升后下降，最終達到穩態.

表6 矩形拼接圖像修復實驗數值結果（MSE）Tab.6 Numerical results（MSE）of rectangular image inpainting experiments

圖8 矩形拼接圖像修復數值實驗相對誤差與總能量的變化情況Fig.8 Changes of relative residuals and total energy in image inpainting numerical experiments for rectangular

在復雜結構圖像修復實驗中，對條形碼和二維碼圖像（Nx=Ny=256）分別設置不同的局部區域信息缺失，見圖9，各算法修復圖像的MSE見表7.

圖9 條形碼和二維碼的真實圖像與待修復圖像Fig.9 Real images of barcode and QR code and images to be inpainting

由表7可見，SAVCN算法的修復效果最優.圖10給出了變量u的相對誤差和總能量EDP（u）的變化曲線.由圖10（a）和（b）可見，對于條形碼圖像，CS2M收斂最快，對于二維碼圖像，2種ADMM算法和CS2M收斂速度基本一致，高于其他算法.由圖10（c）和（d）可見，各算法均可保證目標總能量最終達到穩態.

圖10 條形碼和二維碼圖像修復數值實驗相對誤差與總能量變化情況Fig.10 Changes of relative residuals and total energy in image inpainting numerical experiments for barcode and QR code

表7 條形碼和二維碼圖像修復實驗數值結果（MSE）Tab.7 Numerical results（MSE）of barcode and QR code images inpainting experiments

對于以上圖像去噪和修復實驗，圖11給出了ADMM和SAVCN算法中超參數對圖像復原效果的影響.由圖11（a）可見，ADMM2V的參數r2L對最終的結果數值影響是敏感的.由圖11（b）和（c）可見，SAVCN的參數τ和ADMM3V的參數r1、r2對最終的結果數值影響均不敏感.

圖11 算法參數對圖像復原結果的影響Fig.11 Influences of algorithm parameters on image restoration results

4 結論

本文利用DP模型處理二值圖像復原問題，主要考慮了三種成熟的半隱迭代格式和兩種新的ADMM算法.數值算例的結果表明，本文所提出的ADMM方法去噪實驗中可以恢復到信噪比更高的結果且保持能量穩定，在修復實驗中，SAVCN格式最終修復所得結果更好，但其收斂速度較慢

鑒于ADMM算法在圖像去噪實驗中表現出穩定的優勢，因此在后續的研究中，將對ADMM算法的收斂性進行分析.除此之外，考慮到SAVCN格式在圖像修復中表現良好，從而將在后續研究中設計更為穩定的能量衰減迭代格式.