基于廣義似然比的小波域SAR圖像相干斑抑制算法

侯建華，劉欣達，陳穩(wěn)，陳少波

(中南民族大學 電子信息工程學院，武漢430074)

摘要在聯(lián)合檢測與估計理論框架下推導(dǎo)出了Bayes萎縮函數(shù)表達式，并提出了一種基于廣義似然比的小波域SAR圖像去斑算法.該算法對含斑SAR圖像直接做冗余小波變換，求出小波系數(shù)所對應(yīng)的二值掩模；對相干斑噪聲和有用信號的似然條件概率分別建模為尺度指數(shù)分布和Gamma分布，根據(jù)二值掩模信息，采用最大似然估計得到兩種模型的參數(shù)并計算似然條件概率比.實驗結(jié)果表明：文中所給算法在有效濾除斑點噪聲的同時，也較好地保持了圖像的細節(jié)信息，在對人工加斑圖像和多幅實際SAR圖像的處理中獲得了令人滿意的結(jié)果.

關(guān)鍵詞聯(lián)合檢測與估計；SAR圖像去斑；小波變換；廣義似然比

收稿日期2015-01-20

作者簡介侯建華(1964-)，男，教授，博士，研究方向：圖像處理、模式識別，E-mail：hou8781@126.com

基金項目國家自然科學基金資助項目(61141010， 61201448)；湖北省自然科學基金資助項目(2012FFA113)

中圖分類號TN911.73文獻標識碼A

Generalized Likelihood Ratio Based SAR Image

Speckle Suppression Algorithm in Wavelet Domain

HouJianhua,LiuXinda,ChenWen,ChenShaobo

( College of Electronic Information Engineering, South-Central University for Nationalities, Wuhan 430074, China)

AbstractA Bayes shrinkage formula is derived under the framework of joint detection and estimation theory, and a wavelet SAR image despeckling algorithm is realized based on generalized likelihood ratio. Firstly, redundant wavelet transform is performed directly to the original speckled SAR images, and binary mask is obtained for each wavelet coefficient. We use scale exponential distribution and Gamma distribution, respectively, to model the likelihood conditional probability of speckle noise and useful signal. According to the mask, the parameters of the two modes are estimated by maximum likelihood estimation method, and thus the likelihood conditional probability ratio is calculated. Experiment results show that the proposed method can effectively filter the speckle noise, and at the same time preserve the image details as possible. Satisfactory results are achieved on both synthetically speckled images and real SAR images.

Keywordsjoint detection and estimation；SAR image despeckling；wavelet transform；generalized likelihood ratio

合成孔徑雷達(SAR)在衛(wèi)星遙感領(lǐng)域得到廣泛的應(yīng)用，但由于SAR的相干成像機理，SAR圖像所固有的相干斑噪聲嚴重降低了SAR圖像的可解譯程度，影響了后續(xù)的圖像處理任務(wù)，因此相干斑抑制方法研究成為SAR圖像預(yù)處理中一個非常重要的部分.SAR圖像去斑方法最初采用的是簡單空間域濾波方法，如均值濾波和中值濾波；從20世紀80年代開始到90年代初期，提出了空域自適應(yīng)濾波算法，典型代表為Lee濾波器[1]、Kuan濾波器[2]、Gamma MAP 濾波器[3,4]等；這些算法通過在圖像上取一個滑動窗,以窗內(nèi)所有像素作為濾波器的輸入值,基于窗口內(nèi)的局部統(tǒng)計特性進行濾波處理，獲得了優(yōu)于均值濾波和中值濾波的效果.

從20世紀90年代Donoho等人將小波分析應(yīng)用于圖像去噪開始[5],小波域去斑算法研究已成為SAR圖像去斑的主流方向[6-8].基于小波的SAR圖像去斑算法中研究最多的是小波閾值或萎縮方法，但這種方法沒有充分考慮斑點噪聲小波系數(shù)的統(tǒng)計特性；并且選取適當?shù)拈撝狄彩冀K是一個難題.

為了進一步提高算法的去斑性能，近年來，利用SAR圖像小波系數(shù)的統(tǒng)計分布特性,結(jié)合貝葉斯統(tǒng)計理論的SAR圖像去斑已成為該領(lǐng)域的一個研究熱點.Xie等提出了基于馬爾可夫隨機場的小波域貝葉斯斑點消除方法[9]；Achim等將真實信號所對應(yīng)的小波系數(shù)視為alpha-stable分布，斑點系數(shù)為高斯分布，利用貝葉斯最大后驗概率技術(shù)(MAP)對SAR圖像進行估計[10]；文獻[11]將有用信號、斑點噪聲的系數(shù)分別建模為拉普拉斯分布、瑞利分布，利用Bayesian MAP準則得到了解析的軟閾值去噪表達式.在基于小波統(tǒng)計模型的SAR去斑方法中，關(guān)鍵問題是如何對其中的有用信號和斑點噪聲信號建模，以及如何有效地對真實的信號做出估計，這些方面還有大量的工作有待于探索.

本文在文獻[12]的基礎(chǔ)上，實現(xiàn)了一種快速、有效的SAR圖像去斑算法.該方法沒有采用常規(guī)的對數(shù)變換，而是直接對SAR圖像做冗余小波分解，然后利用分類技術(shù)將小波系數(shù)劃分為兩類；其中斑點噪聲系數(shù)建模為尺度指數(shù)分布，有用信號的系數(shù)則建模為Gamma分布；根據(jù)Bayes最大似然估計得到這兩種模型的參數(shù)；最后基于聯(lián)合檢測與估計理論的廣義似然比(generalized likelihood ratio)框架，得到自適應(yīng)的萎縮函數(shù)表達式.

1基于廣義似然比的Bayes去斑理論

文獻[13]指出：在雷達信號處理、語音傳輸過程中，判斷系統(tǒng)觀察空間中的某一點是否代表信號，這個問題往往是不確定的.為此文獻[13]提出了信號的聯(lián)合檢測與估計理論，該理論已成功應(yīng)用于語音信號處理中.本文就是將聯(lián)合檢測與估計的方法應(yīng)用于小波域SAR圖像去斑.

1.1符號

為了清楚起見，先對符號做統(tǒng)一說明.設(shè)wk,,j,D表示原始SAR圖像在空間位置為k，分辨率尺度為j，圖像細節(jié)方向為D的小波系數(shù)，其去斑后相應(yīng)的估計值為yk,,j,D.下面在不至于混淆的情況下忽略尺度和方向標記j、D.

1.2Bayes萎縮函數(shù)

根據(jù)聯(lián)合檢測與估計理論[13,14]，設(shè)H1表示假設(shè)“小波系數(shù)對應(yīng)的是信號”，H0表示假設(shè)“小波系數(shù)對應(yīng)的是噪聲”.則真實信號的概率密度為：

p(y)=p(H0)·p(y|H0)+p(H1)·p(y|H1),

(1)

y的最小均方誤差估計(MMSE)為：

E(y|w)=p(H0|w) ·E(y|w,H0)+p(H1|w)·E(y|w,H1),

(2)

下面對(2)式化簡[14].首先，在H0假設(shè)成立的條件下，E(y|w,H0)≈0 ，則(2)式簡化為：

E(y|w)≈p(H1|w)·E(y|w,H1),

(3)

進一步，在H1假設(shè)成立的條件下，E(y|w,H1)≈y，則(3)式簡化為下面的萎縮函數(shù)形式：

(4)

在(4)式的具體實施中，用后驗概率p(Xk=1|Mk=mk)代替p(H1|w).考慮到位置k處的該后驗概率既與全局先驗有關(guān)(這里的全局先驗是針對一個子帶而言的)，也與位置k處的局部鄰域信息有關(guān)[14]，因此在計算后驗概率p(Xk=1|Mk=mk)過程中引入k點的鄰域信息，得到：

y=p(Xk=1|Mk=mk,X?k)·w.

(5)

下面推導(dǎo)后驗概率p(Xk=1|Mk=mk,X?k)·w的表達式.根據(jù)Bayes公式:

p(Xk=1|Mk=mk,X?k)=

(6)

上式推導(dǎo)中做了Mk=mk與X?k相互獨立的假設(shè)，其中:

(7)

(8)

ξk與ηk的乘積ξk·ηk稱為廣義似然比.ξk代表似然條件概率之比，其計算依賴于具體的條件概率模型，將在下面討論.ηk代表先驗信息，文獻[12]利用吉布斯隨機場(GRF)模型，提出了一種計算ηk的方法：

(9)

其中參數(shù)γ控制鄰域系數(shù)對算法的影響.最終得到的基于廣義似然比的萎縮函數(shù)表達式為：

(10)

2廣義似然比的計算

根據(jù)(10)式，ηk可以通過當前位置的鄰域信息得到，關(guān)鍵是似然條件概率比ξk的計算.其思想是首先估計每個位置所對應(yīng)的二值掩模；再對兩種假設(shè)下的小波系數(shù)的條件概率密度建模，結(jié)合掩模對模型參數(shù)進行估計，最后得到兩種假設(shè)下的條件概率密度，從而計算出ξk.

2.1掩模估計

圖像經(jīng)過冗余小波分解，每一個小波子帶都具有相同的尺寸.根據(jù)小波系數(shù)的尺度傳遞性，有用信號(如圖像邊緣)小波系數(shù)的幅度隨著分解尺度的增加而增加，而表示噪聲的小波系數(shù)正好系數(shù)則相反.利用該性質(zhì)，在尺度j的每個細節(jié)子帶wj,D中，按下式對二值掩模標記進行估計：

(11)

(12)

2.2似然函數(shù)估計

圖1是真實的SAR圖像(horsetrack)第1尺度HL1子帶中m0、m1的直方圖(實線)，即經(jīng)驗似然函數(shù)p(mk|0)和p(mk|1).

圖1　horsetrack圖像HL1子帶的經(jīng)驗似然函數(shù)估計 (實線)和分析模型曲線(虛線) Fig.1　Empirical likelihood estimation in HL1 subband of horsetrack image(solid line) and the analysis model curve (dashed line)

由圖1，文獻[12]提出用尺度指數(shù)分布、廣義Gamma分布分別近似斑點噪聲系數(shù)、有用信號的系數(shù)m0、m1的經(jīng)驗似然函數(shù)，從圖中可以看到，這兩種理論模型能夠較好的模擬實際直方圖分布.這兩種分布的概率密度函數(shù)分別為：

(13)

(14)

先假設(shè)模型參數(shù)a、b已知，于是廣義似然比為：

(15)

根據(jù)(10)式就可得到基于廣義似然比的萎縮函數(shù)表達式.

2.3模型參數(shù)估計

上面的兩種概率分布都只有一個參數(shù)，在給定觀測數(shù)據(jù)條件下利用最大似然估計(ML)理論可以很容易地求出各自的模型參數(shù).

對于尺度指數(shù)分布，令θ=1/a，則(13)式寫為：p(m,θ)=θ·exp(-θm)；

在給定觀測數(shù)據(jù)m0=(m1,m2,…,mn0)條件下，聯(lián)合概率密度為：

(16)

相應(yīng)的對數(shù)似然函數(shù)為：

(17)

(18)

同理，廣義Gamma分布模型參數(shù)b的參數(shù)的最大似然估計:

(19)

其中n0、n1分別是m0、m1的長度.

圖1同時給出了兩種分析模型對實際直方圖的擬合(虛線)，其中參數(shù)a、b按照(18)、(19)式計算.從圖1可見，尺度指數(shù)分布和廣義Gamma分布可以較好地表示m0、m1的條件概率密度函數(shù).

3去斑算法步驟

上述完整的小波域SAR圖像去斑算法步驟如下.

1) 對含斑圖像做N層冗余小波分解；

3) 對j=1,2，…,N-1尺度的各方向子帶：

iii) 按(18)、(19)式估計模型參數(shù)a、b；

iv) 由(15)式計算ξk，由(9)式計算ηk；

4) 做小波反變換，得到去斑后的恢復(fù)圖像.

具體應(yīng)用中，使用二次樣條小波，做4尺度冗余小波分解.

4實驗結(jié)果

4.1參與實驗的算法說明

我們選取了另外4種常用的算法與本文算法比較：

(1)傳統(tǒng)均值濾波；

(2)同態(tài)Wiener濾波[15]，該方法是在空域?qū)D像做對數(shù)變換后再進行Wiener濾波，最后做指數(shù)變換；

(3)增強型Gamma MAP濾波器，該方法被認為是目前空域自適應(yīng)濾波方法中最好的算法之一，在近年來的相干斑去噪文獻中通常被做為標準算法與其它方法進行比較[14]；

(4)BayesShrink軟閾值去噪算法[16]，記為BayesShrink-RWT；但與文獻[16]不同的是，這里要先對圖像做對數(shù)變換后再進行BayesShrink軟閾處理，最后做指數(shù)變換；并且常用冗余小波變換，它比基于正交小波的BayesShrink算法具有更好的去斑性能.

(5)本文方法，考慮鄰域信息，即根據(jù)(9)式計算ηk，其中γ取值為0.2，窗口為3×3.

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

首先在人工加斑圖像(即模擬SAR圖像)上對各種去斑方法進行比較，然后給出針對真實SAR圖像的實驗結(jié)果和分析.采用客觀評價與主觀評價相結(jié)合的方法綜合衡量算法的去斑性能，其中客觀評價指標為等效視數(shù)(ENL)；等效視數(shù)是衡量一幅SAR圖像相干斑噪聲相對強度的一種指標.等效視數(shù)越大,表明圖像上的相干斑越弱,可解譯性越好.等效視數(shù)可以由SAR圖像中一塊均勻區(qū)域I的均值平方與方差的比值計算，其定義為：

(20)

其中μ與σ分別是SAR圖像中一塊均勻區(qū)域內(nèi)的均值與標準差.

表1給出了在模擬Lena斑點圖像上的性能比較，圖2是各種算法的視覺效果比較.

表1不同去斑算法在模擬Lena斑點圖像上的性能比較

Tab.1Performance comparison of different despeckling

algorithm in Lena synthetically speckled image

去斑算法等效視數(shù)(ENL)斑點圖像19.12均值濾波161.76同態(tài)Wiener濾波151.21增強型Gammaf濾波161.76BayesShrink軟閾值256.47本文方法177.77

從表1可以看出，本文算法的等效視數(shù)僅比BayesShrink軟閾值方法低，高于其它三種算法.但從圖2給出的去斑后視覺效果考察，BayesShrink軟閾值方法處理后的圖像過平滑現(xiàn)象比較嚴重，其它三種算法處理后的圖像則還存在一些明顯的斑點噪聲；相比之下，從視覺效果上看，本文算法在去斑與保持圖像細節(jié)之間的擇中方面具有顯著的優(yōu)勢.

真實SAR圖像選自Sandia Natioanal Laboratories，分別是horse track (396×600)、agricultural field (350×600)、China Lake Airport (686×532).實驗結(jié)果見表2和圖3.

圖2　不同去噪算法在模擬Lena斑點圖像(SNR=11.75)上濾波后的視覺效果比較 Fig.2　Visual comparison of different despeckling algorithms on synthetically speckled image Lena (SNR=11.75)

去斑算法HorsetrackAgriculturalChinaLackAirport原始斑點圖像15.9413.8623.43均值濾波76.6659.0872.97同態(tài)Wiener濾波73.3051.1464.55增強型GammaMAP濾波76.6659.3371.13BayesShrink軟閾值去噪61.1545.7363.53本文算法121.85106.03118.11

圖3　本文算法對不同真實SAR圖像的去斑效果 Fig.3　Despecking effects of the proposed method on different real SAR images

從表2和圖3中可看出該算法能夠有效地濾除斑點噪聲，同時圖像的強、弱細節(jié)得以較好地保持.由此可見本文算法在對人工加斑圖像和多幅實際SAR圖像處理中都獲得了令人滿意的結(jié)果.

5結(jié)語

本文基于聯(lián)合檢測與估計的理論框架，設(shè)計并實現(xiàn)了一種基于廣義似然比的小波域SAR圖像去斑算法.該方法通過二值掩模將小波系數(shù)劃分為兩類：有用信號和斑點噪聲信號，利用實際SAR圖像小波系數(shù)的統(tǒng)計分布特性，分別對這兩類信號進行建模并求出模型參數(shù)，最終根據(jù)所推導(dǎo)的Bayes萎縮函數(shù)對SAR圖像做去斑處理.實驗結(jié)果表明，本文算法能有效地抑制斑點噪聲、同時又能較好的保持圖像的紋理細節(jié)和邊緣信息，能夠滿足后續(xù)圖像處理任務(wù)的要求.

參考文獻

[1]Lee J S. Refined filtering of image noise using local statistics[J]. Computer graphics and image processing, 1981, 15(4): 380-389.

[2]Kuan D T, Sawchuk A A, Strand T C, et al. Adaptive noise smoothing filter for images with signal-dependent noise[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1985 (2): 165-177.

[3]Kuan D T, Sawchuk A, Strand T C, et al. Adaptive restoration of images with speckle[J]. IEEE Transactions on Acoustics, Speech and Signal Processing, 1987, 35(3): 373-383.

[4]Lopes A, Nezry E, Touzi R, et al. Structure detection and statistical adaptive speckle filtering in SAR images[J]. International Journal of Remote Sensing, 1993, 14(9): 1735-1758.

[5]Donoho D L. Ideal spatial adaptation via wavelet shrinkage[J]. Biometrica, 1994, 81: 424-455.

[6]侯建華. 基于小波及其統(tǒng)計特性的圖像去噪方法研究 [D]. 武漢: 華中科技大學, 2007.

[7]Argenti F, Lapini A, Bianchi T, et al. A tutorial on speckle reduction in synthetic aperture radar images[J]. IEEE Geoscience and Remote Sensing Magazine, 2013, 1(3): 6-35.

[8]侯建華, 毛曉暉. 基于正態(tài)反高斯先驗?zāi)Ｐ偷男〔ㄈピ胨惴╗J]. 中南民族大學學報: 自然科學版, 2009 (3): 85-89.

[9]Xie H, Pierce L E, Ulaby F T. SAR speckle reduction using wavelet denoising and Markov random field modeling[J]. IEEE Transactions on Geoscience and Remote Sensing, 2002, 40(10): 2196-2212.

[10]Achim A, Tsakalides P, Bezerianos A. SAR image

denoising via Bayesian wavelet shrinkage based on heavy-tailed modeling[J]. IEEE Transactions on Geoscience and Remote Sensing, 2003, 41(8): 1773-1784.

[11]Jian-Hua Hou, Xiang-Ming Liu, Cheng-Yi Xiong,

Xiang He. Speckle reduction algorithm for synthetic aperture radar images based on Bayesian maximum a posteriori estimation in wavelet domain [J]. Optical Engineering, 2008, 47(5): 057004-1 — 057004-11.

[12]Pizurica A, Philips W, Lemahieu I, et al. Despeckling SAR images using wavelets and a new class of adaptive shrinkage estimators[C]//IEEE.Proceedings of International Conference on Image Processing. San Francisco: IEEE, 2001, 2: 233-236.

[13]Middleton D, Esposito R. Simultaneous optimum detection and estimation of signals in noise[J]. IEEE Transactions on Information Theory, 1968, 14(3): 434-444.

[14]Pizurica A. Image denoising using wavelets and spatial context modeling[D]. Ghent : Ghent University, 2002.

[15]Achim A, Bezerianos A, Tsakalides P. Novel Bayesian multiscale method for speckle removal in medical ultrasound images[J]. IEEE Transactions on, Medical Imaging, 2001, 20(8): 772-783.

[16]Chang S G, Yu B and Vetterli M. Adaptive wavelet thresholding for image denoising and compression[J]. IEEE Transactions on Image Processing, 2000, 9(9): 1532-1546.