程凡強，朱永貴

(中國傳媒大學數(shù)據(jù)科學與智能媒體學院，北京 100024)

1 引言

線性逆問題在天體物理學、光學、信號/圖像處理和統(tǒng)計學等學科建設(shè)中應(yīng)用廣泛，反問題的跨學科性質(zhì)通過大量的文獻得以證明并迅速發(fā)展。本文使用快速迭代收縮閾值算法(FISTA)結(jié)合小波變換得到一種新的算法，使用該算法對彩色圖像進行去噪處理。對于這一類問題，在研究之前，需要從一個基本的線性逆問題入手，但此問題需要離散化以后才能進行，下面就對該問題進行離散化。這個離散的線性系統(tǒng)的形式為

Ax=b+w

(1.1)

當A∈Rm×n和b∈Rn已知時，w是未知噪聲(擾動)向量，x是待估計的未知真信號/圖像。例如，在圖像模糊問題中，b∈Rn代表模糊圖像，x∈Rn是未知真圖像，假設(shè)其大小與b相同，即m=n。b和x都是通過列字典序的方式形成的。在這些應(yīng)用中，矩陣A代表了模糊算子，由觀測到的模糊噪聲圖像b估計真圖像x的問題稱為圖像去模糊問題。

1.1 背景

解決問題(1.1)的經(jīng)典方法是最小二乘法(LS)[1]，用以下公式使數(shù)據(jù)誤差達到最小化，即

當m=n(如圖像處理應(yīng)用中的特殊情況)以及A是非奇異矩陣時，LS只是輔助解A-1b。在圖像去模糊中，A經(jīng)常是病態(tài)的[2]，在這種情況下，LS的解通常有一個巨大的范數(shù)，因此導(dǎo)致所得到的解毫無意義。為了克服這個困難，需要用正則化方法使原問題變成“良態(tài)”問題，從而使解穩(wěn)定、唯一。正則化的基本思想是用一個“鄰近的”條件良好的問題來代替原來的條件不好的問題，則該問題的解就近似于所需要的解。其中最開始流行的一種正則化是由Tikhonov提出的，后來被命名為Tikhonov正則化[3]，又稱為L2正則化方法，其實質(zhì)就是加了二范數(shù)的平方，即

(1.2)

上述最小化問題中的第二項是控制解的范數(shù)的正則化項。正則化參數(shù)λ>0在測量的保真度和噪聲靈敏度之間提供了一個權(quán)衡。L最常見選擇是接近一階或二階導(dǎo)數(shù)算子的恒等式或矩陣；可以參見文獻[4，5，6]。

后來，在研究圖像或信號處理中，人們又發(fā)現(xiàn)了另一種正則化項，叫作L1正則化，它的形式是

(1.3)

1.2 貢獻

凸優(yōu)化問題(1.3)可通過內(nèi)點法求解[11]。然而，在大多數(shù)應(yīng)用中，問題幾乎都是大規(guī)模的密集的矩陣數(shù)據(jù)，這基本排除了使用復(fù)雜的內(nèi)點法。因為這些問題的存在，推動了人們對求解(1.3)的更簡單的算法搜索。經(jīng)過多年探索，人們找到了一種公認的效果較好的方法，那就是迭代收縮閾值算法，其中每次迭代都涉及矩陣向量乘法，并且涉及和之后是收縮軟閾值步驟；請參見文獻[12，13，14，15]。

2 基本事實

在給出這個優(yōu)化算法之前，先給出一些基于梯度的基本事實。

2.1 梯度方法

考慮連續(xù)可微函數(shù)f∈Rn→R的無約束最小化問題，如下

(U) min{f(x)：x∈Rn}

解決(U)的最簡單方法之一是梯度法，其生成序列{xk}是由

x0∈Rn，xk=xk-1-tk▽f(xk-1)

(2.1)

其中tk>0是一個合適的步長。梯度迭代(參見文獻[16，17])可以被看作線性化函數(shù)f在xk-1處的近似正則化[18]，并且等價為

2.2 迭代收縮閾值算法(ISTA)

在求解擁有L1正則化項的問題時(這是這類問題的一種特殊情況)，最常用的方法就是ISTA方法，它被看作是梯度法的擴展，因此，ISTA方法所采用的思想與梯度法基本一致。在每次迭代時都用到矩陣向量乘法，然后涉及和之后就是軟閾值步驟。

下面給出ISTA的一般步驟是

xk+1=Γλt(xk-2tAT(Axk-b))

(2.2)

其中Γ是一個合適的步長，且Γα：Rn→Rn是由下式定義的運算符

Γα(x)i=(|xi|-α)+sgn(xi)Γ

(2.3)

在眾多的優(yōu)化文獻中，ISTA的收斂性分析已經(jīng)在不同的文獻中得到了很好的研究，包括各種修改與改進；可以參見文獻[7，9，19]，它的重點是確定序列{xk}收斂到(1.3)解的條件。ISTA的優(yōu)點在于它的簡單性；但卻忽視了本身收斂的緩慢性。

這就是該算法的核心部分，那么接下來給出完整的ISTA算法：

ISTA算法:輸入:L:=L(f) - ▽f的Lipschitz常數(shù)第零步:給定初始值x0∈Rn.第k(k≥1)步:計算xk=pL(xk-1)(2.4)

其中pL(·)是更一般的ISTA方法。

在下一節(jié)中，將給出本文所做的具體貢獻，在原有算法的基礎(chǔ)上進行改進，對彩色圖像進行處理，得到的效果很是不錯。

3 彩色圖像處理

在給出彩色圖像處理的程序與效果圖之前，我們還需要知道以下基本事實。

ISTA方法的出現(xiàn)，對于求解線性逆問題是一個極大的突破，很好的解決了“病態(tài)”問題，能夠使這類問題都有唯一解。然而，美中不足的就是收斂速度確實太慢，所以后人又在此基礎(chǔ)上進行了改進，發(fā)現(xiàn)并提出了很多種方法，但是應(yīng)用較多的或者是最常用方法就是FISTA。

當正則化項為零時，模型(1.2)或(1.3)和ISTA就簡化為梯度法。這是一種“最優(yōu)”的一階方法。值得注意的是，文獻[20]中開發(fā)的方法在每次迭代時不需要多一個梯度評估(即與梯度方法相同)，只是進行了一次巧妙的選擇，找到了一個易于計算的附加點而已。

在本節(jié)中，將會呈現(xiàn)新方法，既然作為ISTA的改進版，它必須保持原有的計算簡單性，而且還有比原來更快的收斂速度，顯然這些都已經(jīng)滿足。

3.1 小波去燥模型

在各種各樣的模型中，小波去燥模型是研究的重中之重，這種模型的離散化形式如下

(3.1)

其中B∈Rm×n是觀測到的模糊噪聲圖像，W：Rm×n→Rm×n是小波變換，X是原始圖像(“真”圖像)，λ是正則化參數(shù)。

處理L1范數(shù)正則化準則的基本原理是大多數(shù)圖像在小波域中具有稀疏表示。在去噪領(lǐng)域中，小波理論受到了許多學者的重視，他們應(yīng)用小波進行去噪并獲得了非常好的效果。這樣就不得不說一下小波變換所具備的特點：

(1)低熵性。小波系數(shù)的稀疏分布，使得圖象變換后的熵降低，從而更好的刻畫細節(jié)特征；

(2)多分辨率。由于采用了多分辨率的方法，所以可以非常好的刻畫信號的非平穩(wěn)特征，如邊緣、尖點、斷點等；

(3)去相關(guān)性。因為小波變換可以對信號進行去相關(guān)，且噪聲在變換后有白化趨勢，所以小波域比時域更有利于去噪；

(4)選基靈活性。由于小波變換可以靈活選擇變換基，從而對不同應(yīng)用場合，對不同的研究對象，可以選用不同的小波母函數(shù)，以獲得最佳的效果。

由于其特點，在后來的去噪模式中，逐步取代了傅里葉變換的地位。小波去噪的方法主要包括三個基本的步驟：第一，對含噪聲信號/圖像進行小波變換；第二，對變換得到的小波系數(shù)進行某種處理，以去除其中包含的噪聲；第三，對處理后的小波系數(shù)進行小波逆變換，得到去噪后的信號。

注：小波去噪方法的不同之處集中在第一步。

3.2 快速迭代收縮閾值算法(FISTA)

當求解問題(3.1)時，采用ISTA方法是個較為不錯的方法，只是迭代較為緩慢，因此，在此基礎(chǔ)上進行了改進，找到一個加速因子，把它應(yīng)用到ISTA方法中，得到的這個新的方法在計算量上并沒有增加多少，但是收斂速度成幾何倍增加，效果也令人滿意。

下面就給出FISTA算法：

FISTA算法:輸入:L:=L(f) - ▽f的Lipschitz常數(shù)第零步:給定初始值y1=x0∈Rn,t1=1.第k(k≥1)步:計算xk=pL(yk)(3.2)tk+1=1+1+4t2k2(3.3)yk+1=xk+(tk-1tk+1)(xk-xk-1)(3.4)

上述算法與ISTA的主要區(qū)別在于迭代收縮算子pL(·)不僅只用前一個點xk-1，而是使用前兩個點，即xk-2，xk-1，這樣，就能得到這個算法中的點yk。顯然這兩個方法的主要計算工作量并沒有發(fā)生變化，即在運算符pL(·)中。而(3.3)和(3.4)中tk+1，yk+1對所需的額外計算就顯得微不足道了。

4 算法代碼及實驗結(jié)果

在本節(jié)中，將給出基于小波去噪的彩色圖像FISTA的算法及對應(yīng)實驗結(jié)果。為了更好的表現(xiàn)所給代碼的去噪效果，我們把原圖一起進行比較，比較結(jié)果如下：

clc

close all

clear all

addpath(‘utils’)；

X=

double(imread(‘images/woman1.jpg’))；%double(imread(‘images/Lena.bmp’))；%double(imread(‘images/cameraman.pgm’))/255；

X=X/255；

subplot(1，3，1)

imshow(X，[ ])，title(‘Original Image’)；

sigma=0.1；

Bobs=X+sigma*randn(size(X))；

B=Bobs；

subplot(1，3，2)

imshow(B，[ ])，title(‘Noised Imgae’)

iwav=@(x)Wavedb1Phi(x，1)；% Inverse Wavelet Transform

wav=@(x)Wavedb1Phi(x，0)；% Wavelet Transform

pars.MAXITER=20；

alpha_1=12；%1e-4；

X_out=

denoise_wavelet_FISTA(B，wav，iwav，1/alpha_1，pars)；%denoise_wavelet_FISTA(B，wav，iwav，1e-4，pars)；

u=X_out；

subplot(1，3，3)

imshow(u，[])，title(‘Denoised Results’)；

Rel=norm(u(：)-X(：))*100/norm(X(：))

SNRValue=SNRZHU1(X，u)

SNRValue1=snr(u，X)

圖1 Rel=15.0587

圖2 SNRValue=16.4443

圖3 SNRValue1=12.5203

從上面的圖片1-3中不難看出，去噪效果還是很明顯的，而且很好的保留了原圖像的細節(jié)，圖1為原始圖像、圖2加噪圖像和圖3去噪圖像，上述結(jié)果分別為去噪圖像與原始圖像的差、去噪圖像與原始圖像的信噪比(兩種方法)，一般來說，信噪比指的是放大器的輸出信號與同時輸出的噪聲的比。信噪比越大，說明混在信號里的噪聲越小，聲音回放的音質(zhì)量越高，否則相反。

5 總結(jié)

在圖像處理應(yīng)用中，我們不能直接作用于圖像，或者說通常意義下的圖片，因此需要將其模型化處理，這樣就與線性反問題聯(lián)系到一起了。在解決這類問題的時候，用傳統(tǒng)的方法(即最小二乘法)解決起來非常麻煩，而且經(jīng)常求不出解，這就是我們所說的“病態(tài)”問題，為了得到最優(yōu)解，我們需要將方程進行處理，然后得到我們需要的結(jié)果。目前較為流行的方法就是加正則項，使原來的“病態(tài)”問題變成“良態(tài)”問題，這樣就能求解了，而且解也是唯一的。當然正則化項是多種多樣的，如本文提到的L1、L2正則化項，還有未涉及到的TV正則化項，等，這些都是根據(jù)不同情況加上的正則化項。

本文解決的問題是在小波變換下對含有L1正則化項的反問題，通過FISTA算法處理彩色圖像，從主觀來說，基于小波去噪的彩色圖像FISTA算法修復(fù)效果更好，具有更強的競爭力。這種方法的出現(xiàn)使得圖像處理變得更加簡單、高效、清晰。