畢思文,陳 浩,帥 通,李 娜
(1.中國(guó)電子科技集團(tuán)公司航天信息應(yīng)用技術(shù)重點(diǎn)實(shí)驗(yàn)室,河北 石家莊 050081;2.中國(guó)科學(xué)院遙感與數(shù)字地球研究所,北京100101)
數(shù)字圖像在成像和傳輸過(guò)程中,經(jīng)常會(huì)受到各種各樣的噪聲干擾[1]。為了使后期的圖像理解和圖像分割等操作更加有效,需要對(duì)受到噪聲干擾的圖像進(jìn)行處理。
在圖像去噪方面,人們提出了許多圖像去噪算法,文獻(xiàn)[2] 提出基于Cycle Spinning的圖像自適應(yīng)閾值去噪方法,該方法提高了去噪圖像峰值信噪比(PSNR),降低了均方誤差(MSE),去噪圖像清晰,獲得較好的視覺(jué)效果。 文獻(xiàn)[3]通過(guò)平穩(wěn)小波變換對(duì)圖像進(jìn)行小波分解,對(duì)于子圖像的高頻區(qū)域進(jìn)行閾值分割和雙邊濾波,利用平穩(wěn)小波更好的冗余性和平穩(wěn)不變性,更好地去除了SAR圖像的相干斑噪聲,實(shí)驗(yàn)表明這種改進(jìn)的去噪方法對(duì)SAR圖像的相干斑噪聲有很好的抑制效果。文獻(xiàn)[4]提出了一種新方法來(lái)濾除圖像中的云霧,結(jié)果表明,該算法優(yōu)于傳統(tǒng)濾波方法,且對(duì)于薄云圖像效果更佳。 文獻(xiàn)[5]在經(jīng)典小波閾值去噪算法的基礎(chǔ)上改進(jìn)了閾值函數(shù),提出了一種新的小波閾值去噪算法,實(shí)驗(yàn)結(jié)果表明該算法提高了信號(hào)特征的可分離性,具有較高的實(shí)用價(jià)值。
小波分析是為了彌補(bǔ)短時(shí)傅里葉變換的不足而發(fā)展起來(lái)的一門(mén)應(yīng)用數(shù)學(xué)學(xué)科[6]?;谛〔ǚ治龅娜ピ敕椒ㄗ钤缬蒑allat提出,他在1992年建立了小波變換快速算法,并將其運(yùn)用在信號(hào)和圖像的分解和重構(gòu)中[7]。1999年,Kingsbury提出了雙樹(shù)復(fù)小波變換(Dual Tree Complex Wavelet Transform,DTCWT),具有平移不變性,提供了6個(gè)方向的信息,因而具有較好的方向性和精確的相空間信息[8]。在DTCWT提出以后,有很多學(xué)者在基于雙樹(shù)復(fù)小波的去噪方法方面做了大量研究。文獻(xiàn)[9]提出了一種基于雙樹(shù)復(fù)小波變換和形態(tài)濾波的去噪算法。文獻(xiàn)[10]提出了一種基于非下采樣雙樹(shù)復(fù)小波域的圖像去噪算法,實(shí)驗(yàn)表明該算法比經(jīng)典算法提高了一定的峰值信噪比,且有良好的視覺(jué)效果,較好地保持了圖像中的紋理特征。文獻(xiàn)[11]綜合考慮空域?yàn)V波和變換域?yàn)V波的優(yōu)點(diǎn),提出了一種基于DTCWT和自適應(yīng)窗的圖像去噪算法,將雙樹(shù)復(fù)小波變化和自適應(yīng)橢圓窗口濾波相結(jié)合,考慮了小波分解的不同子帶具有不同的能量方向,以橢圓方向窗作為領(lǐng)域,獲得了比較好的去噪效果。但文獻(xiàn)[11]中使用的橢圓窗作為鄰域,系數(shù)估計(jì)時(shí)均需要采用橢圓模板進(jìn)行匹配,所以算法運(yùn)行時(shí)間較長(zhǎng)、復(fù)雜度略高。
本文借鑒多量子位量子疊加態(tài)的測(cè)量坍縮原理,根據(jù)雙樹(shù)復(fù)小波較好的方向特性,以坍縮后的狀態(tài)作為鄰域計(jì)算小波系數(shù)方差,利用雙樹(shù)復(fù)小波提供的方向信息和量子疊加態(tài)的測(cè)量坍縮原理運(yùn)用到圖像去噪中進(jìn)一步去噪。實(shí)驗(yàn)結(jié)果表明,本文提出的圖像去噪算法與文獻(xiàn)[11]的圖像去噪方法相比,去噪性能得到改善,運(yùn)行時(shí)間明顯提升。
DTCWT變換可以通過(guò)兩對(duì)濾波器組同時(shí)作用在輸入數(shù)據(jù)上來(lái)實(shí)現(xiàn)。復(fù)小波可以表示為:
ψ(t)=ψr(t)+jψi(t),
(1)
式中,Ψr(t)表示復(fù)小波的實(shí)部;Ψi(t)表示復(fù)小波的虛部。Ψr(t),Ψi(t)都是實(shí)函數(shù),因此,DTCWT可以表示為2個(gè)獨(dú)立的實(shí)小波變換,包含了2個(gè)平行的小波樹(shù):樹(shù)a和樹(shù)b。DTCWT分解示意圖如圖1所示,樹(shù)a和樹(shù)b的疊加濾波器組分別表示復(fù)數(shù)小波變換的實(shí)部和虛部,↓2表示隔點(diǎn)采樣[12]。為了保證濾波器的沖擊響應(yīng)對(duì)應(yīng)于復(fù)小波變換系數(shù)的實(shí)部和虛部,采用2棵樹(shù)的濾波器長(zhǎng)度分別為奇數(shù)和偶數(shù)且是線性相位。基本原理就是利用一對(duì)實(shí)濾波樹(shù),同時(shí)對(duì)輸入信號(hào)進(jìn)行分解,產(chǎn)生小波系數(shù)的實(shí)部和虛部。
圖1 DT-CWT分解示意
二維的雙樹(shù)復(fù)小波實(shí)部與虛部小波系數(shù)能提取±15°,±45°,±75°六個(gè)方向的高頻信息,相對(duì)離散小波變換(Discrete Wavelet Transform,DWT)具有近似平移不變性、多方向選擇性、更高定位精度和計(jì)算效率等優(yōu)點(diǎn)[13]。相比之下,DWT在每個(gè)尺度上有3個(gè)小波子帶,只能反映出水平數(shù)豎直和對(duì)角方向。二維DT-CWT在空間域和在2-D平面內(nèi)(理想化的)表示出具有6個(gè)不同的角度DT-CWT小波如圖2所示。
圖2 DT-CWT的脈沖響應(yīng)及在二維平面等效
若|Ψ1〉是Hilbert空間中的一個(gè)矢量,|Ψ2〉是另外一個(gè)矢量,由態(tài)疊加原理可知:
|Ψ〉=c1|Ψ1〉+c2|Ψ2〉,
(2)
也是Hilbert空間中的一個(gè)矢量。其中c1,c2分別是狀態(tài)|Ψ1〉|Ψ2〉的概率幅。且滿足歸一化條件:
c12+c22=1。
(3)
所以若量子系統(tǒng)處在|Ψ1〉和|Ψ2〉描述態(tài)中,則式(2)的線性疊加態(tài)|Ψ〉也是該系統(tǒng)的一個(gè)可能態(tài),這就是量子力學(xué)的態(tài)疊加原理[14]。量子比特[15](qubit)是量子信息理論中的一個(gè)重要概念,對(duì)一個(gè)具有2個(gè)基態(tài)的雙態(tài)量子系統(tǒng)[16]。若將2個(gè)基態(tài)分別記為|0〉和|1〉,記量子比特為:
|Ψq〉=a|0〉+b|1〉。
(4)
|Ψ〉= |Ψ(1)〉 |Ψ(2)〉…|Ψ(n)〉=
(5)
式中,|i〉表示第i個(gè)基態(tài);ωi為基態(tài)|i〉的概率幅,滿足歸一化條件:
(6)
由量子力學(xué)第三假設(shè)可知,設(shè)測(cè)量算子由{Mm}描述,測(cè)量前量子系統(tǒng)的最新?tīng)顟B(tài)是|Ψ〉,則測(cè)量后系統(tǒng)的狀態(tài)為[17]:
(7)
為了更好地理解量子測(cè)量坍縮原理,舉例如下:對(duì)于一個(gè)4×4的疊加態(tài)結(jié)構(gòu)元素(歸一化以后),若另取一個(gè)4×4的陣列,將陣列中邊緣的最大的2個(gè)灰度值取為1,其余取值為0,則得到{Mm},對(duì)應(yīng)|iM〉=|0001000000001000〉,若用此測(cè)量算子Mm=|0001000000001000〉〈00010000000010000|對(duì)鄰域圖像進(jìn)行測(cè)量,則鄰域圖像將坍縮到基態(tài)|iM〉,如圖3所示。
圖3 測(cè)量前和測(cè)量坍縮以后
由于圖像中不同點(diǎn)之間領(lǐng)域的特征不一樣,所以在不同移動(dòng)點(diǎn)所得到的測(cè)量算子和測(cè)量后的坍縮態(tài)也不同[18]。以上就是雙樹(shù)復(fù)小波變換以后,計(jì)算+45°方向鄰域,其他方向類(lèi)似。
假設(shè)原始圖像被均值為0,方差為σ2的高斯白噪聲污染,當(dāng)在小波域通過(guò)小波變換系數(shù)進(jìn)行去噪后采用以下模型:
Y(i,j)=X(i,j)+N(i,j),
(8)
式中,Y(i,j)為含噪小波系數(shù);X(i,j)為待估計(jì)的干凈小波系數(shù);N(i,j)為噪聲小波系數(shù)。采用MAP預(yù)估器可表示為:
(9)
(10)
(11)
將式(10)帶入式(9),并對(duì)其進(jìn)行求導(dǎo),令其導(dǎo)數(shù)為零可得:
(12)
式中,M為鄰域N中系數(shù)的個(gè)數(shù)。若已知圖像信號(hào)的方差分布,那么由最大后驗(yàn)概率估計(jì)(Maximum A Posterior,MAP)可得:
(13)
式中,fσ(σ2)為圖像信號(hào)的方差分布。由式(13)推導(dǎo)可得:
(14)
根據(jù)以上算法模型,本文去噪算法步驟如下:
① 對(duì)圖像進(jìn)行雙樹(shù)復(fù)小波變換;
② 利用量子態(tài)疊加的測(cè)量坍縮原理,以4×4的陣列測(cè)量分解層;
③ 以坍縮后的鄰域用式(8)和式(11)計(jì)算λ;
④ 重復(fù)步驟②,然后利用式(13)估計(jì)雙樹(shù)復(fù)小波系數(shù)的方差;
⑤ 利用式(8)得到去噪以后的系數(shù);
⑥ 雙樹(shù)復(fù)小波反變換,得到去噪后的圖像。
為了驗(yàn)證本文所提算法的有效性,實(shí)驗(yàn)中選用標(biāo)準(zhǔn)的512×512大小的8位灰度圖像作為實(shí)驗(yàn)對(duì)象。實(shí)驗(yàn)中,假定用均值為0、標(biāo)準(zhǔn)差分別為10,15,20高斯白噪聲污染,對(duì)測(cè)試圖像進(jìn)行4層的下采樣雙樹(shù)復(fù)小波分解。在文獻(xiàn)[19]中指出,對(duì)包含了高斯噪聲的圖像進(jìn)行3級(jí)小波分析,并對(duì)高斯噪聲小波系數(shù)的分布情況進(jìn)行了統(tǒng)計(jì)分析,發(fā)現(xiàn)大約有96%的噪聲系數(shù)分布在最外2層的細(xì)尺度子帶中[20]。因此,對(duì)第二、第三和第四分解層,采用利用量子態(tài)疊加測(cè)量坍縮原理,選取4×4的陣列中四周最大的2個(gè)值取值為1,其他的取為0,以此作為雙樹(shù)復(fù)小波分解以后的主方向,并且在此方向上作為系數(shù)估計(jì)的鄰域。將本文所提算法和文獻(xiàn)[11]做比較(PSNR 和運(yùn)行時(shí)間)。
實(shí)驗(yàn)效果圖和文獻(xiàn)[11]的數(shù)據(jù)對(duì)比結(jié)果如圖4和圖5所示。
圖4 2種去噪算法結(jié)果比較(Lena)
圖5 2種算法去噪結(jié)果比較(Barbara)
表 1 2種不同去噪方法得到的PSNR值比較
算法LenaBarbara1015202510152025平均運(yùn)行時(shí)間/s文獻(xiàn)[11]算法30.528.627.325.633.230.228.826.27.6本文算法30.828.627.525.833.330.428.926.62.1
通過(guò)表1可以看出,利用本文所提算法得到的PSNR值比文獻(xiàn)[11]要略高,從處理以后的效果圖來(lái)看,本文所提算法有相同的抑制噪聲的能力,但由于本文所提算法將量子態(tài)疊加的測(cè)量坍縮原理與雙樹(shù)復(fù)小波變換相結(jié)合,將坍縮狀態(tài)作為計(jì)算子帶的方向窗,既利用了雙樹(shù)復(fù)小波變換的方向特性,同時(shí)避免了文獻(xiàn)[11]中每次系數(shù)估計(jì)時(shí)都要進(jìn)行模板匹配的操作,因此,在運(yùn)行時(shí)間上比文獻(xiàn)[11]有很大的提高。
本文將雙樹(shù)復(fù)小波變換和量子態(tài)疊加的測(cè)量坍縮原理相結(jié)合,利用雙樹(shù)復(fù)小波變換較好的方向特性和疊加態(tài)結(jié)構(gòu)元素在不同位置坍縮為不同大小和形狀的屬性,避免了方形結(jié)構(gòu)元素所不具備的方向性特性和橢圓結(jié)構(gòu)元素帶來(lái)的復(fù)雜性。通過(guò)實(shí)驗(yàn)發(fā)現(xiàn),將量子力學(xué)中的理論運(yùn)用到圖像去噪中,大大提升了運(yùn)行速度并驗(yàn)證了本文算法的有效性。該算法雖然在運(yùn)行時(shí)間上相對(duì)于其他方法有較高提升,但是在評(píng)價(jià)指標(biāo)PSNR數(shù)據(jù)上的提高并不明顯,需要對(duì)算法進(jìn)行優(yōu)化改進(jìn),仍有較大的提升空間。