光學遙感圖像變量分裂迭代快速復原算法

宋義剛，肖亮，韋志輝，2，黃麗麗

(1.南京理工大學 計算機科學與技術(shù)學院，江蘇 南京210094;2.南京理工大學 應(yīng)用數(shù)學系，江蘇 南京210094)

0 引言

光學遙感成像是對地觀測和深空探測的重要數(shù)據(jù)獲取途徑。然而成像過程中光學系統(tǒng)散焦、運動和大氣擾動等造成的模糊效應(yīng)以及成像系統(tǒng)光電因素引起的噪聲干擾等嚴重影響了圖像的視覺質(zhì)量和空間分辨率。如何去除模糊效應(yīng)和抑制噪聲干擾是光學遙感圖像復原需要解決的問題，也是光學遙感圖像超分辨的關(guān)鍵技術(shù)［1］。

目前圖像復原算法可分為頻域算法和空間域算法。Wiener 濾波算法簡單、便于實現(xiàn)、在噪聲水平不大的情形能夠獲得較好的復原效果，但隨著噪聲水平增大復原效果急劇下降并存在較多的“寄生波紋效應(yīng)”;約束最小二乘算法具有Wiener 濾波類似的頻域公式，能夠減少“寄生波紋效應(yīng)”但是不能很好的保持圖像的邊緣［2］。文獻［3］基于Contourlet變換的遙感圖像去噪新算法，為遙感圖像噪聲抑制提供了新的手段。文獻［4］研究提出了一種基于超分辨圖像處理的模糊圖像盲復原算法，其利用推定復原圖像的拉普拉斯算子獲得參數(shù)的估計，算法用于遙感圖像復原取得較好的效果。文獻［5］提出將多準則優(yōu)化理論引入到光學層析成像圖像重建取得很好的結(jié)果，對光學遙感圖像復原也具有指導意義。

光學遙感圖像復原是一個病態(tài)重建問題，為克服重建過程的病態(tài)性，人們提出基于正則化和偏微分方程的圖像復原方法，如文獻［7－8］提出基于各向異性擴散方程的圖像濾波算法，文獻［9］針對航天湍流降質(zhì)圖像給出極大似然估計正則化復原算法。Rudin－Osher 提出了圖像去噪的全變差(TV)圖像模型［10］。TV 模型將圖像建模為有界變差空間(BV 空間)的函數(shù)，能夠很好保持圖像的邊緣結(jié)構(gòu)，成為廣為認可的圖像復原和超分辨領(lǐng)域的圖像正則化模型［10－14］。但是基于TV 模型的復原算法還存在一些問題，包括:復原圖像會存在“階梯效應(yīng)”;如果采取最速下降法求解，算法的收斂速度較低，經(jīng)實驗表明需要迭代上百次才能達到收斂，并且每次迭代的復雜度非常高，降低了該算法的實用性。最近一些學者提出了一些快速算法，例如:針對TV 去噪模型的快速投影算法［11］，多重網(wǎng)格算法［12］，基于引入輔助變量的代價函數(shù)方法［13］等。

本文給出了一種新的變量分裂快速迭代快速復原算法。首先通過引入待去模糊圖像的輔助變量、圖像梯度的輔助變量，構(gòu)造逼近于原始TV 復原模型的代理代價函數(shù)。通過分析代理代價函數(shù)特殊結(jié)構(gòu)，充分利用循環(huán)周期卷積核與圖像的卷積等價于傅里葉域的乘積，并且由傅里葉變換不會改變歐式范數(shù)的性質(zhì)，設(shè)計了三階段解析迭代算法。算法具有如下優(yōu)點:1)迭代過程每個階段均是具有解析解的子問題求解，并且每次迭代僅需2 次快速傅里葉變換和1 次快速逆傅里葉變換，算法的時間復雜度僅為O(M2logM)，其中M2為圖像像素個數(shù)。2)迭代過程考慮了人眼對圖像平坦區(qū)域和邊緣區(qū)域噪聲的感知特性，通過估計迭代系統(tǒng)中自適應(yīng)的正則化參數(shù)，改善了圖像復原的效果，一定程度上克服了“階梯效應(yīng)”。

1 已有全變差正則化復原算法與問題分析

1.1 光學遙感成像降質(zhì)模型

下面討論一般遙感成像系統(tǒng)的降質(zhì)模型。記F(h)為點擴散函數(shù)h 的傅里葉變換，并令(ξ，η)∈為頻域坐標，遙感成像光學系統(tǒng)的調(diào)制傳遞函數(shù)(MTF)一般包含如下3 部分:

1)感光耦合器件(CCD)的調(diào)制傳遞函數(shù)。在焦平面上CCD 敏感像素是一個細小的［－ c/2，c/2］2窗口，其中c 為相鄰CCD 像元的距離，可以用一個矩形函數(shù)來描述，這樣可推導出CCD 的MTF為

可假設(shè)遙感成像光學系統(tǒng)為圓對稱系統(tǒng)，則總體MTF 是各部分MTF 的乘積:

假設(shè)u(x，y)為未降質(zhì)的圖像，從而降質(zhì)觀測圖像u0(x，y)為

或 u0(x，y)=u(x，y)?h(x，y)+n，

式中:(x，y)為空間坐標;(u，v)為頻率域坐標;F 和F－1分別表示傅里葉正反變換;n 表示噪聲;?表示卷積算子。

1.2 基本模型與問題分析

圖像復原實質(zhì)上從u0中估計合適的u，屬于Hardmard 意義下非適定數(shù)學反問題?；赗udin－Osher 連續(xù)全變差模型，人們提出了圖像正則化復原算法［11－13］，

式中:參數(shù)α 為正的拉格朗日乘子，連續(xù)全變差模型定義為

利用變分法和梯度下降法(GD)轉(zhuǎn)換為求解非線性方程［10］:

式中:H 為h 對應(yīng)的卷積算子，為避免|u| =0，需引入?yún)?shù)ε ＞0，將其替換為全變差圖像復原模型能夠減少“寄生波紋效應(yīng)”也能較好保持圖像邊緣，但是會導致圖像中出現(xiàn)“階梯效應(yīng)”。

另外雖然梯度下降法是目前求解非線性變分問題最直觀、最常用且實現(xiàn)比較簡單的算法，但算法存在如下缺點:1)收斂速度較慢;2)在求解過程中需引入?yún)?shù)ε ＞0，且因為參數(shù)ε 的引入迭代解不能嚴格收斂于(6)式的精確解;3)每次迭代都要計算一次H*H，由于線性空間不變的模糊算子的支集一般很大，當受到模糊影響的遙感圖像尺寸很大時去模糊的計算代價更高。文獻［13］針對該問題，通過引入輔助變量w 和w 與u 的懲罰項，構(gòu)造逼近于模型(2)式的代理代價函數(shù)，

顯然當參數(shù)α1足夠大時，能夠逼近于模型(2)式，并可采取交錯迭代法求解:

其好處是將原始TV 復原轉(zhuǎn)化為去模糊和TV去噪2 個子問題。而去模糊是二次凸泛函的最優(yōu)，可通過共扼梯度法快速求解;而TV 去噪可采取Chambolle 快速投影算法。從而避免了每次迭代都要計算一次H*H，且Chambolle 投影算法是TV 去噪的精確求解算法，實驗結(jié)果表明該算法收斂速度提高一個數(shù)量級。該算法對TV 去噪引起的“階梯效應(yīng)”仍沒有很好解決。

2 代理函數(shù)模型與快速優(yōu)化算法

2.1 代理函數(shù)模型建立

本文在離散化框架下建立基于全變差的代理代價函數(shù)圖像復原模型。不失一般性，考慮M ×M 大小的離散觀測圖像u0∈RM2，理想圖像u∈RM2和噪聲n∈RM2.光學系統(tǒng)的點擴散函數(shù)為h，則降質(zhì)模型為

定義離散梯度算子D∶RM2→RM2為

這里i=M，j=M 時采取循環(huán)邊界條件，wi，j表示向量w 的第(iM+j)個分量。則

因此原始基于TV 的圖像復原模型的離散化形式為

令w=h?u，則降質(zhì)模型

因此光學系統(tǒng)圖像的圖像復原問題等效于先去噪然后去模糊的2 個子問題迭代求解。根據(jù)TV 模型，本文構(gòu)建帶約束優(yōu)化模型

式中:‖·‖22為2－范數(shù)。利用Lagrange 乘子法將有約束問題轉(zhuǎn)化為無約束代價函數(shù)最優(yōu):

式中α1、α2為正的正則化參數(shù)。比較(12)式和(9)式，2 式之間的差別是(12)式新增了變量w 和相關(guān)的懲罰項‖w－u?h‖22;增加了另一個變量v 以及相關(guān)的懲罰項(‖Dxu－v1‖22+‖Dyu－v2‖22)。當α1、α2充分大時，(12)式與(9)式等價。

注意到(12)式中的代價函數(shù)包含4 項，第1 項為l1范數(shù)，而第2、第3 和第4 項為l2范數(shù)，它們都是凸的，因此代價函數(shù)是一個凸優(yōu)化問題。因此根據(jù)凸分析理論，(12)式存在唯一的最優(yōu)解。下面對新的代理代價函數(shù)最優(yōu)設(shè)計具體的迭代過程和快速算法。

2.2 子問題優(yōu)化求解

使用交替最小化算法求解(12)式。由一初值u(0)，利用該方法得到一迭代序列

使得

進一步考慮到:1)點擴散函數(shù)(PSF)與圖像的卷積等價于調(diào)制傳遞函數(shù)與圖像傅里葉變換的乘積;2)采取循環(huán)邊界條件的向前差分算子與圖像作用的過程可以表示為線性卷積過程，同樣可以轉(zhuǎn)化為向前差分算子形成的循環(huán)矩陣的傅里葉變換與圖像傅里葉變換的乘積;3)傅里葉變換屬于線性變換，不會改變歐式范數(shù)(Parseval 定理)的性質(zhì)。則記F(u)，F(xiàn)(u0)，F(xiàn)(w)，F(xiàn)(Dx)，F(xiàn)(Dy)，F(xiàn)(h)分別為u，u0，w，Dx，Dx，h 的傅里葉變換，并分別令α =因此最小化問題(13)～(15)式可重寫為

下面分別對3 個子問題的最優(yōu)解進行推導，并給出計算公式。

v－子問題(13)式:u(i－1)固定時，直接進行演算，該問題對應(yīng)的解為

w－子問題(17)式:當u(i－1)固定時，其目標代價函數(shù)是關(guān)于F(w)的二次泛函，其解可以顯式的表示為

由(20)式可以看出，F(xiàn)(w(i))可以看作是第i－1 次的校正頻譜F(u(i－1))與原始圖像頻譜的插值，相當于進行了一次頻譜校正。

u－子問題(18)式:將F(v(i))和F(w(i))固定，求(18)式的最小解可以顯示表示為

其中(20)式和(21)式“?”表示矩陣點積，“* ”表示復共扼。

進一步進行逆傅里葉變換可得

(22)式表明，u－子問題等效于圖像高頻信息與w－子問題求解得到的插值校正頻譜信息F(w(i))的頻譜重新校正過程。整個過程表明三階段解析迭代過程具有較好的頻譜外推特性。

2.3 參數(shù)選擇

從上節(jié)可以看出，構(gòu)建原始全變差圖像恢復模型的代理代價函數(shù)的好處是等效地將模型(9)式轉(zhuǎn)化為3 個簡單子問題的迭代求解過程，而3 個子問題最優(yōu)均具有顯式解析解。三階段解析迭代僅需要輸入降質(zhì)圖像u0，正則化參數(shù)α1，α2， (λ 或者等效參數(shù)以及迭代終止誤差參數(shù)。

由于本文算法中目標代價函數(shù)當α1，α2充分大時，才能很好地逼近于全變差復原模型。因此基于延拓方法的思想［14］，本文設(shè)置α1，α2的最大上限值α1max和α2max，在迭代過程中逐漸調(diào)整α1，α2的值達到最大上限。終止條件采取α1，α2的最大上限值α1≥α1max或α2≥α2max以及Re Diff= ‖u(i+1)－u(i)‖2/‖u(i+1)‖2＜ξstop進行控制。

基于全變差模型的圖像復原算法能夠較好地克服“寄生波紋效應(yīng)”，但是全變差模型容易產(chǎn)生“階梯效應(yīng)”，特別是參數(shù)λ 的選擇對“階梯效應(yīng)”的產(chǎn)生有影響;不恰當?shù)摩?容易導致平坦區(qū)域的“階梯效應(yīng)”非常明顯，從而降低復原圖像的信噪比。根據(jù)代理函數(shù)模型，λ 的作用是平衡正則化項和保真項‖w－u0‖22懲罰之間的權(quán)衡參數(shù)。按照人眼對圖像的感知特性，人眼對圖像刺激的響應(yīng)較少依賴于絕對亮度，而更依賴于對周圍亮度的局部對比度變化。Weber－Fechner 定律［15］表明在一個相當寬的變化范圍內(nèi)，主觀上剛可辨別亮度差異Δu 與背景的亮度u 之比稱為Weber 比，它是一個常數(shù)(約為2%).另一方面圖像平坦區(qū)域由于其本身方差很小且人眼對該區(qū)域的“階梯效應(yīng)”最為敏感，可認為方差的增加純粹是噪聲和“階梯效應(yīng)”引起的，因此應(yīng)增加正則化項的作用，即選取較小的λ;而邊緣區(qū)，由于其方差本身較大以及人眼的掩蓋效應(yīng)，使得人眼對該區(qū)域的“階梯效應(yīng)”不敏感，因此應(yīng)增加保真項的懲罰，選取較大的λ，紋理區(qū)域應(yīng)介于邊緣區(qū)和平坦區(qū)之間。λ 的選擇可通過圖像的局部方差和梯度信息進行。由于在算法迭代過程中由(19)式我們可得到圖像的梯度信息vki，j，因此本文構(gòu)造自適應(yīng)

式中σi，j采取魯棒性方差自動估計方法［16］，σi，j=1.486·median‖vki，j－median|vki，j|‖.

2.4 快速算法描述與復雜度分析

綜合3.2 和3.3 節(jié)，本文算法可歸結(jié)為如下步驟:

初始化:u0=u0，正則化參數(shù)α1，α2，λ，最大上限值α1max和α2max，迭代終止誤差ξstop.

解析迭代

while α1＜α1max或α2＜α2maxdo

1)固定u(i－1)，根據(jù)(19)式計算梯度信息vii，j;校正得到F(w(i));

3)空域vki，j的梯度信息進行FFT 得到F(v1(i))，

2)根據(jù)(20)式在傅里葉變換域進行線性頻譜F(v2(i))，連同步驟2)得到的F(w(i))代入(21)式計算F(u(i))，并進行快速IFFT 得到u(i);

4)根據(jù)(23)式進行參數(shù)更新。

終止條件判斷:

如果Re Diff =‖u(i+1)－u(i)‖2/‖u(i+1)‖2＞ξstop，α1= 2α1，α2= 2α1，按照(23)式自適應(yīng)調(diào)整(i，j)，重新迭代。否則輸出u(i)，算法結(jié)束。

本文算法屬于變量分離的子問題分裂算法，包含3 個子問題的求解。v－子問題對應(yīng)的解析求解vki，j的算法復雜度為O(M2);而w－子問題和u－子問題的過程中，F(xiàn)(Dx)，F(xiàn)(Dy)，F(xiàn)(h)*，F(xiàn)(h)在迭代過程中保持不變，可以預先計算;解析求解F(w(i))和F(u(i))需要2 次快速傅里葉變換和一次快速逆傅里葉變換，其復雜度為O(M2lgM2)=O(M2lgM).而算法本質(zhì)上是利用傅里葉變換對循環(huán)矩陣的對角化特性，與共扼梯度法一樣屬于有限步收斂算法。因此算法的整體復雜度為O(M2lgM).

關(guān)于本文算法的空間復雜度，由于預先存儲的F(Dx)，F(xiàn)(Dy)，F(xiàn)(h)*，F(xiàn)(h)等信息均與圖像的大小相同，因此空間復雜度為O(M2).

3 實驗結(jié)果與分析

文中所有算法均在MATLAB 7.1 平臺實現(xiàn)，計算環(huán)境為IBM W510 筆記本電腦、Intel 1.73 GHz CPU，3 G 內(nèi)存。為了驗證本文算法的有效性。用于定量評估去模糊算法的性能指標采取改進信噪比(ΔSNR)、峰值信噪比(PSNR)和相對誤差(ReErr).設(shè)u0，u ，u* 分別為M ×N 的降質(zhì)圖像、參考圖像和去模糊圖像，ΔSNR、PSNR 和ReErr 定義如下:

一般而言改進信噪比、峰值信噪比越大越好，而相對誤差越小越好。為了更好地評價算法的結(jié)構(gòu)保持特性，本文采取結(jié)構(gòu)相似度指標(SSIM)來評價圖像的恢復質(zhì)量［17］。SSIM 指標與人類視覺特性基本吻合，與人的主觀測評圖像質(zhì)量有較好的匹配性。實驗中對比算法包括:MATLAB 7.1 中Wiener濾波、約束最小二乘算法(RLS)，Lucy－Richardson 算法;同時我們和全變差交替子空間投影算法(TVAP)［13］和文獻［14］中FTVd 算法進行的對比。本文算法中初始化參數(shù)α1=α2=320，α1max=α2max=40 000，參數(shù)λ 初值為1.0，迭代過程由(23)式估計。TVAP、FTVd 和本文算法迭代終止條件均取為Re Diff=‖u(i+1)－u(i)‖2/‖u(i+1)‖2＜10－4.實驗中圖1(a)為一幅清晰的512 ×512 光學遙感圖像，遙感成像系統(tǒng)MTF 模擬曲線如圖1(b)所示。圖2(a)為根據(jù)降質(zhì)模型模擬生成的降質(zhì)圖像，其中加入噪聲為均值為0、方差為15 高斯噪聲。從圖2(b)－(f)所示復原圖像的視覺質(zhì)量來看，約束最小二乘(RLS)算法和Lucy－Richardson 算法對于“寄生波紋效應(yīng)”有所改善，但“階躍”邊緣和線狀邊緣保持效果不好。FTVd 算法和TVAP 的結(jié)果圖像能夠極大地減少“寄生波紋效應(yīng)”，但視覺上的“階梯效應(yīng)”較為明顯。本文算法的復原圖像(如圖2(f))很好地克服了兩方面的不利影響，并且復原圖像的細節(jié)較為清晰。

圖1 實驗圖像與成像系統(tǒng)MTFFig.1 The test image and MTF of the imaging system

從復原算法的圖像質(zhì)量客觀評價指標來看，本文算法雖然與TVAP 同樣利用了全變差模型，由于TVAP 基于Chambolle 投影［11］的精確求解算法，圖像質(zhì)量的客觀評價指標有所改善。而本文算法中通過自適應(yīng)調(diào)整參數(shù)(i，j)，圖像質(zhì)量的客觀評價指標得到進一步改善。表1給出了關(guān)于樣本圖像圖2(a)的各類指標的提升情況，可以看到本文算法的PSNR 比FTVd 提高0.6 dB，而本文算法的ΔSNR比FTVd 提高0.3 dB.另外，本文算法取得了較好的SSIM 指標，表明圖像的結(jié)構(gòu)保持性能很好。

從算法所花費的CPU 時間來看，Lucy－Richardson 算法和RLS 算法是經(jīng)典迭代的算法，時間較長，分別需要34.76 s 和7.11 s.TVAP 算法和FTVd、本文算法雖然都屬于迭代算法，但因為前者采取了Chambolle 投影算法，而FTVd 采用了變量分離算法，都是目前快速全變差算法，時間較短。而本文算法僅用了1.36 s 就達到優(yōu)于FTVd 算法的圖像質(zhì)量，因此本文算法得到了很好的圖像質(zhì)量－計算時間性價比。

圖2 不同算法的復原圖像Fig.2 Recovered images of different algorithms

表1 不同復原算法性能指標比較結(jié)果Tab.1 Quantitative comparison between the proposed algorithm and state－of－the art algorithms

4 結(jié)論

本文給出了基于全變差圖像復原模型的一種三階段解析迭代快速復原算法，該算法基于輔助變量的全變差代理代價函數(shù)模型，將圖像復原轉(zhuǎn)化為3個子優(yōu)化問題;并通過分析代理代價函數(shù)特殊的結(jié)構(gòu)，得到了算法時間復雜度僅為O(M2logM)的三階段解析迭代快速復原算法。算法考慮了人眼對圖像平坦區(qū)域和邊緣區(qū)域噪聲的感知特性，通過估計迭代系統(tǒng)中自適應(yīng)的正則化參數(shù)，較好克服了“階梯效應(yīng)”，改善了圖像復原的效果。實驗證明算法是一種適應(yīng)于非紋理圖像的優(yōu)良復原算法。需要指出，對于實際光學遙感圖像的復原，一般需要對模糊系統(tǒng)進行辨識，本文算法是假設(shè)系統(tǒng)點擴散函數(shù)已知的圖像非盲復原，對于融合模糊辨識的盲復原算法將是下一步的主要工作;另外聯(lián)合稀疏性和正則性［18］的圖像復原也是重要的研究方向。

References)

［1］ Mohammal D A.Super－resolution:a short review，a new method based on hidden Markov modeling of HR image and future challenges［J］.Computer Journal，2009，52(1):126－141.

［2］ 鄒謀炎.反卷積和信號復原［M］.北京:國防工業(yè)出版社，2001.ZOU Mou－yan.De－convolution and signal restoration［M］.Beijing:National Defense Industry Press，2001.(in Chinese)

［3］ 張晶晶，方勇華.基于Contourlet 變換的遙感圖像去噪新算法［J］.光學學報，2008，28(3):462－466.ZHANG Jing－jing，F(xiàn)ANG Yong－h(huán)ua.Novel denoising method for remote sensing image based on contourlet transform［J］.Acta Optica Sinica，2008，28(3):462－466.(in Chinese)

［4］ 劉揚陽，金偉其，蘇秉華，等.基于超分辨力圖像復原算法的模糊系統(tǒng)辨識［J］.光電子·激光，2005 ，16(2):213－216.LIU Yang－yang ，JIN Wei－qi ，SU Bing－h(huán)ua，et al.Identification of blurred system based on super－resolution reconstructing image schemes［J］.Journal of Optoelectronics Laser，2005，16(3):213－216.(in Chinese)

［5］ 孟靜，王加俊，黃賢武，等.一種光學層析圖像的多準則重建方法［J］.光學學報，2006，26(9):1340－1344.MENG Jing，WANG Jia－jun，HUANG Xian－wu，et al.Multi－criterion reconstruction method for optical tomography［J］.Acta Optica Sinica，2006，26(9):1340－1344.(in Chinese)

［6］ Shao W Z，Wei Z H.Edge－and－corner preserving regularization for image interpolation and reconstruction［J］.Image and Vision Computing，2008，26:1591－1606.

［7］ 王懷野，張科，李言俊.各向異性濾波在紅外圖像處理中的應(yīng)用［J］.紅外與毫米波學報，2005，24(2):109－112.WANG Huai－ye，ZHANG Ke，LI Yan－jun.Anisotropic Gaussian filtering for infrared image［J］.Journal of Infrared And Millimeter Waves，2005，24(2):109－112.(in Chinese)

［8］ 洪漢玉，張?zhí)煨?，余國?航天湍流降質(zhì)圖像的極大似然估計規(guī)整化復原算法［J］.紅外與毫米波學報，2005，24(2):130－134.HONG Han－yu，ZHANG Tian－xu，YU Guo－liang.Regularized restoration algorithm of astronautically turbulence degraded images using maximum－likehood estimation ［J］.Journal of Infrared and Millimeter Waves，2005，24(2):130－134.(in Chinese)

［9］ Rudin L，Osher S，F(xiàn)atemi E.Nonlinear total variation based noise removal algorithms［J］.Physica D，1992，60(1－4):259－268.

［10］ Marquina A，Osher S J.Image super－resolution by TV－regularization and Bregman iteration［J］.Journal of Scientific Computing，2008，37:367－382.

［11］ Chambolle A.An algorithm for total variation minimization and application［J］.Journal of Mathematical Imaging Vision，2004，20(1):89－97.

［12］ Chen K，Tai X C.A nonlinear multigrid method for total variation minimization from image restoration ［J］.Journal of Scientific Computing，2007，33(2):115－138.

［13］ Huang Y M，Micheal K NG，Wen Y W.A fast total variation minimization method for image restoration［J］.Muti－scale Model Simulation，2008，7(2):774－795.

［14］ Wang Y，Yang J，Yin W，et al.A new alternating minimization algorithm for total variation image reconstruction［J］.SIAM Journal on Imaging Sciences，2008，1(3):248－272.

［15］ Shen J H.On the foundations of vision modeling，I.Weber’s law and Weberized TV restoration［J］.Physica D:Nonlinear Phenomena，2003，175(3－4):241－251.

［16］ Black M，Sapiro G.Edge as outiers:anisotropic smoothing using local image statistics［C］∥Proceedings of the Scale－space Conference.Berlin:Springer－Verlag，1999:259－270.

［17］ Wang Z，Bovik A C，Sheikh H R，et al.Image quality assessment:from error visibility to structural similarity［J］.IEEE Transaction on Image Processing，2004，13(4):600－612.

［18］ Afonso M V，Bioucas－Dias J M，F(xiàn)igueiredo M A T.Fast image recovery using variable splitting and constrained optimization［J］.IEEE Transaction on Image Processing，2010，19(9):2345－2356.