侯建華，陳 穩(wěn)，劉欣達，陳少波( 中南民族大學 智能無線通信湖北省重點實驗室，武漢 430074 )

?

基于異質性分類的小波域SAR圖像去斑

侯建華，陳 穩(wěn)，劉欣達，陳少波
( 中南民族大學 智能無線通信湖北省重點實驗室，武漢 430074 )

摘要：提出了一種非同態(tài)濾波框架下的小波域SAR圖像相干斑抑制算法。將小波域中的后向散射信號和斑點噪聲分別建模為正態(tài)逆高斯分布、高斯分布，在貝葉斯最大后驗準則下推導出信號估計的表達式；為了提高模型參數的估計精度，引入多尺度局部變差系數作為異質性測度，并提出用對數正態(tài)分布對該測度進行擬合；基于異質性測度的統(tǒng)計分布特點，本文對小波子帶中的系數進行分類，用累積量的估計方法計算每一類小波系數的模型參數。實驗結果表明，與傳統(tǒng)的同態(tài)濾波框架下的同類算法，以及與未采用分類技術的非同態(tài)濾波同類算法相比較，本文算法在主、客觀性能評價上均有一定的優(yōu)勢，取得了較好的去斑效果。本文提出的基于異質性測度對小波系數分類的思想為SAR圖像去斑算法的研究提供了一種新的途徑。

關鍵詞：SAR圖像去斑；正態(tài)逆高斯分布；異質性測度；分類

0 引 言

由于合成孔徑雷達(Synthetic Aperture Radar, SAR)圖像固有的相干成像特點，在獲取SAR圖像時不可避免地會引入相干斑，嚴重地阻礙了后續(xù)的圖像處理工作，如圖像分割、圖像分類、目標識別和邊緣檢測等等。因此，相干斑抑制是SAR圖像預處理中一個重要的研究課題。

近二十年來，小波變換因其所具有的諸多良好特性，在SAR圖像去斑處理中得到了廣泛的應用。由于相干斑噪聲的乘性特征，需要首先將其變換為加性噪聲的形式，為此有兩種方法：1) 直接對含斑圖像取對數變換；2) 在形式上將含斑圖像分解為信號與噪聲相加，但該噪聲實際與信號相關。前者稱為同態(tài)濾波方法，后者稱為非同態(tài)濾波。同態(tài)濾波方法的一個優(yōu)勢在于可以借用成熟的去除加性噪聲的理論和技術，因此早期的SAR圖像去斑算法大多屬于此類方法[1-2]，但由于對數變換的非線性，這類方法會造成同質區(qū)域后向散射系數的損失。近年來，非同態(tài)濾波算法受到了更多研究者的關注[3-5]。

異質性測量與分析在SAR圖像處理領域得到了成功的應用，異質性反映了圖像紋理信息變化程度，主要用來描述同一目標區(qū)域內不同場景之間的差別[6]。目前異質性測量方法主要包括：基于變差系數[7]、基于算術與幾何均值比[8]、基于信息論的方法[9]，這些方法已被應用于小波域SAR圖像去斑，且其效果還是比較顯著的。例如，文獻[10]提出了多尺度局部變差系數的概念；文獻[5]提出了以“紋理”能量作為異質性測度，對小波系數進行分類的思想。與[5]類似，文獻[11]也以異質性測度為依據對小波系數進行分類，并根據異質性對閾值做自適應調整。此外，文獻[6]在空間域從信息論的角度出發(fā)，提出了一種新的定量化測量SAR圖像異質性參數的指標，并用于斑點噪聲抑制。

本文提出了一種將小波域貝葉斯去噪與異質性測量相結合的SAR圖像去斑算法。采用非同態(tài)濾波框架，對小波域后向散射信號和斑點噪聲的統(tǒng)計分布建模，求出Bayesian MAP估計表達式；將文獻[10]提出的多尺度局部變差系數作為異質性測度對小波系數分類，分別計算不同類別小波系數所對應的模型參數，從而提高模型參數的估計精度，最終改善去斑算法的性能。

1 小波域貝葉斯SAR圖像去斑

1.1 Bayesian 最大后驗概率估計

假設斑點噪聲完全發(fā)育，則SAR圖像可以表示為以下乘性噪聲模型：

式中：Y為噪聲污染的觀測信號，即含斑圖像；X為待估計的真實無噪信號，Z為與X相互獨立的相干斑噪聲，對上式做平穩(wěn)小波變換(SWT)：

式中：wY、wX及wN分別對應Y、X和N=X( Z -1)的小波系數。根據最大后驗概率估計(MAP)準則，在給定觀測數據wY條件下，使得后驗概率最大的wX即為真實信號的估計：

與自然圖像類似，真實的SAR圖像一般是由大部分的平滑區(qū)域加小部分的異質區(qū)域(邊緣和紋理)構成，平滑區(qū)域小波系數趨近于零，而異質區(qū)域小波系數幅值偏大，從直方圖分布上看，各細節(jié)子帶中的小波系數表現(xiàn)出高峰度、尖峰型和重拖尾的非高斯性，因此高斯分布已經不能很好的描述這種分布特性。正態(tài)逆高斯分布(Normal Inverse Gaussian, NIG)是由具有不同均值的高斯分布和反高斯分布構成的混合模型，在理論上，其參數選取的靈活性有助于克服傳統(tǒng)先驗模型的缺陷[12]。近年來，已經有文獻報道將NIG應用于乘性斑點噪聲的處理[4，13]。NIG概率分布密度函數(pdf)為

其中：K1為第二類一階修正貝塞爾函數，α參數控制分布陡峭度，δ為尺度參數。

假設后向散射信號小波系數wX的先驗分布PwX( wX)為正態(tài)逆高斯分布，如式(4)所示，斑點噪聲的小波系數wN的概率密度函數為

將式(4)、式(5)代入式(3)，參考Hyvarinen的方法[14]，得到：

1.2 參數估計

其中：Zσ為乘性噪聲標準差，h為小波分解濾波器組中的高通濾波器。

對于正態(tài)逆高斯模型，我們采用文獻[16]提出的基于二階累積量的參數估計法。由于該方法又要用到待估計量wX(即無噪圖像)，為此，與文獻[11]類似，我們先用frost濾波器處理原始SAR圖像，得到參考圖像，再對其作平穩(wěn)小波變換，令該信號小波系數為，用此估計正態(tài)逆高斯模型的參數。

根據文獻[17]，NIG模型的形參α、尺度參數δ的估計表達式分別為

其中：k2、k4分別代表二階、四階累積量。k2和k4的統(tǒng)計無偏估計表達式為[18]

其中：M代表小波系數樣本數，m2及m4分別為樣本小波系數的2階和4階中心矩。

2 異質性測量與分析

2.1 異質性測度及統(tǒng)計分布

SAR圖像可以分為同質區(qū)域和紋理區(qū)域，理想的去斑效果應該是在平滑同質區(qū)域的同時保護紋理細節(jié)信息，但抑制噪聲與保持圖像細節(jié)兩者一般是難以兼得的。一個自然的思路就是將像素劃分到不同類別的區(qū)域、再采取不同的處理策略[5，11]。異質性反映了SAR圖像中不同場景之間圖像紋理信息的變化程度，本文選擇Argenti等提出的多尺度局部變差系數(MLCV)作為異質度測量指標，以此為依據對小波系數分類，再按類別估計正態(tài)逆高斯模型參數，提高模型參數的估計精度。異質度測度定義如下[10]：

圖1 Horsetrack圖像第一層子帶MLCV經驗直方圖以及對數正態(tài)分布圖Fig.1 Histogram of MLCV in the first level of Horsetrack and logarithmic normal distribution pdf

我們統(tǒng)計真實SAR圖像異質性測度直方圖，發(fā)現(xiàn)MLCV呈不對稱性分布，先迅速遞增到最大值后再逐漸減小，稱該最大值所對應的c(k)為模態(tài)T。圖1給出了Horsetrack圖像第一層子帶MLCV統(tǒng)計直方圖。進一步，考慮到該直方圖的非對稱性和重拖尾特性，我們提出用對數正態(tài)分布對MLCV直方圖進行擬合，見圖1中的實線部分。從圖中可以看出，對數正態(tài)分布曲線能夠較好的擬合實際的異質性測度直方圖，并且由于對數正態(tài)分布模型簡單，參數估計也很容易得到，因此是一個理想的異質性測度分布理論模型。

2.2 小波系數分類

研究表明，模態(tài)T可以作為不同異質性區(qū)域分類的依據，當c(k)小于T，其對應位置處像素的異質度偏低；相反則像素異質度則偏高。本文提出將每個小波子帶分類為同質區(qū)域、若干個異質區(qū)域，以及孤立散射點區(qū)域，每個系數的分類方法如下：

對于同質區(qū)域、以及n個異質區(qū)域中的像素，我們假設每類像素分別服從不同參數的正態(tài)逆高斯分布，分別估計每一類的NIG參數。對于強異質類區(qū)域，相干斑未完全發(fā)育，不滿足乘性噪聲模型(1)，在本文去斑處理中，其值保持不變。

圖2顯示的是Horsetrack圖像第一層子帶經過異質性分類，共分了四類，an分別取1、1.5，4。每個子圖中的白色區(qū)域代表本類像素集合，黑色區(qū)域則為其他類像素集合，參考原始真實圖像可以發(fā)現(xiàn)，(a)圖中白色區(qū)域對應著原始真實圖像中的同質區(qū)域，而(b)、(c)、(d)圖中白色區(qū)域則反映了不同異質度的紋理區(qū)域。因此，作為異質性測度指標，MLCV為小波圖像分類提供了另外一種有效的手段。

圖2 Horsetrack圖像第一層子帶基于異質性分類圖Fig.2 The classification maps of Horsetrack for sub-band in the first level

3 實驗結果與分析

為了驗證本文算法的有效性，分別對人工加斑圖像(即模擬SAR圖像)、實際SAR圖像進行相干斑抑制實驗，主要從濾波效果，邊緣保護和量化評估等方面分析討論。參與的比較實驗有：1) 經典的空域ehgamma MAP濾波；2) BayesShrink軟閾值去噪算法[19]; 3) 基于NIG的同態(tài)濾波[13]，以下簡稱為“homomorphic-NIG”；4) 基于NIG的非同態(tài)濾波，該算法與本文算法的區(qū)別在于未采用異質性分類技術，以下簡稱 “Non-homomorphic-NIG”；實驗參數設置如下：ehgamma MAP空域濾波算法的局部窗口尺寸選為5×5；四種基于小波的算法選擇4尺度小波分解；本文算法在對小波細節(jié)子帶進行區(qū)域劃分時，通過實驗我們選取a1=1，a3=1.5，a3=4經驗值將每層細節(jié)子帶分為同質區(qū)域，兩類異質區(qū)域及孤立散射區(qū)域四類。

3.1 人工加斑圖像

采用對無噪圖像人工加斑的方法獲得含斑圖像，由于有真實的參考圖像(即無噪圖像)，我們可以獲得客觀的性能評價指標。對512×512 Lena圖像人為加入乘性斑點噪聲，產生一幅SNR為13.17的含斑圖像。

我們選取四種常用的相干斑抑制指標對上述五種算法做定量評價，包括：信噪比(SNR)，均值(MEAN),等效視數(ENL)和邊緣保持指數(EPI)。經上述五種算法處理后圖像的視覺效果如圖3所示。從圖3可以看出，經過ehgamma濾波、BayesShrink算法處理后的圖像，還存留了少數的斑點噪聲；基于NIG的非同態(tài)濾波方法則有所改善。相干斑抑制能力最好的是本文方法和基于NIG的同態(tài)濾波方法，與前者(基于NIG的同態(tài)濾波方法)相比，本文算法利用了異質性分類技術，兼顧了相干斑抑制與紋理細節(jié)的保護，因此處理過的圖像更加清晰，得到了滿意的去斑效果。

去斑性能實驗結果如表1。從表1可以看出，本文算法在上述四個評價指標上均優(yōu)于其它去斑算法，在去斑能力和保護圖像細節(jié)的折中方面取得了最好的效果。此外，基于NIG的同態(tài)濾波算法沒有采用基于異質性的分類技術，但利用統(tǒng)計模型較好的描述了小波系數的分布特性，該算法綜合評價指標僅次于本文算法，也取得了較好的效果。

圖3 五種算法對仿真SAR圖像濾波后的對比Fig.3 Visual comparison of various speckle reduction methods on synthetically speckled image

表1 五種算法對仿真SAR圖像去斑性能比較Table 1 Performance comparison of 5 methods on synthetically speckled image

3.2 真實SAR圖像

對于真實SAR圖像，由于無法獲得真實的參考圖像，因此采用等效視數和均值作為客觀評價指標。3幅真實SAR圖像選自Sandia National Laboratories，分別為Horsetrack、Agricultural、Pipeline。圖4給出了五種算法對Horsetrack圖像處理后的視覺效果對比；與前面人工加斑圖像實驗中的分析類似，從相干斑抑制能力考察，本文方法和基于NIG的同態(tài)濾波方法優(yōu)于其他的方法；但通過觀察可以發(fā)現(xiàn)，基于NIG的同態(tài)濾波方法處理后圖像(圖4(d))在右上角的紋理區(qū)域，相對于本文方法(圖4(f))有一定程度的邊緣損失，同時本文方法處理后圖像的清晰度更高。

圖4 五種算法對Horsetrack圖像處理后的視覺效果Fig.4 Visual comparison of various speckle reduction methods on real SAR image “Horsetrack”

表2給出了實驗結果，從中可看出：1) 本文方法的等效視數均較大于其他四種算法，表現(xiàn)出良好的相干斑抑制能力。2) 從均值考慮，本文方法與原始圖像的均值最為接近；本文方法處理后圖像的均值與基于NIG非同態(tài)濾波幾乎相同，明顯高于其他算法。

表2 五種算法對真實SAR圖像去斑性能比較Table 2 Performance comparison of five methods for real SAR image

4 結 論

本文提出了一種新的SAR圖像去斑算法，其本質上屬于基于小波統(tǒng)計模型的Bayesian方法。與以往算法不同之處在于：1) 采用了非同態(tài)濾波框架，將正態(tài)逆高斯函數和高斯函數作為后向散射信號和斑點噪聲的先驗分布；2) 引入多尺度局部變差系數(MLCV)作為異質度測量指標，研究了MLCV的統(tǒng)計分布特性，提出用對數正態(tài)分布對MLCV進行建模；3) 基于異質性測度，將子帶中的小波系數劃分到具有不同異質性程度的類別，再按照各類別分別估計正態(tài)逆高斯模型參數。實驗結果表明，該方法不僅利用了小波系數的統(tǒng)計先驗知識，還較好的兼顧了邊緣細節(jié)信息的保護；其中分類估計的方法對提高算法的性能起到了重要的作用。

參考文獻：

[1] BHUIYAN M I H，AHMAD M O，SWAMY M N S. Spatially adaptive wavelet-based method using the Cauchy prior for denoising the SAR images [J]. IEEE Transaction on Circuits and Systems for Video Technology(S1051-8215)，2007，17(4)： 500-507.

[2] HOU Jianhua，LIU Xiangming，XIONG Chengyi，et al. Speckle reduction algorithm for synthetic aperture radar images based on Bayesian maximum a posteriori estimation in wavelet domain [J]. Optical Engineering(S0097-3286)，2008，47(5)： 057004-1-057004-11.

[3] FOUCHER S，BéNIéG B，BOUCHER J M. Multiscale MAP filtering of SAR images [J]. IEEE Transaction on Image Processing(S1057-7149)，2001，10(1)：49-60.

[4] SOLB? S，ELTOFT T. Γ-WMAP：a statistical speckle filter operating in the wavelet domain [J]. International Journal of Remote Sensing(S0143-1161)，2004，25(5)：1019-1036.

[5] BIANCHI T，ARGENTI F，ALPARONE L. Segmentation-based MAP despeckling of SAR images in the undecimated wavelet domain [J]. IEEE Transaction on Geoscience and Remote Sensing(S0196-2892)，2008，46(9)：2728-2742.

[6] CHEN Jie，ZHU Jing，LUO Chun-sheng，et al. A novel speckle filter for SAR images based on information-theoretic heterogeneity measurements [J]. Chinese Journal of Aeronautics(S1000-9361)，2009，22(5)：528-534.

[7] LOPES A，TOUZI R，NEZRY E. Adaptive speckle filters and scene heterogeneity [J]. IEEE Transaction on Geoscience and Remote Sensing(S0196-2892)，1990，28(6)：992-1000.

[8] MARTIO BEAUCHEMIN，KEITH P B THOMSON，GEOFFREY EDWARDS. The ratio of the arithmetric to the geometric mean：A first-order statistical test for multilook SAR image homogeneity [J]. IEEE Transaction on Geoscience and Remote Sensing(S0196-2892)，1996，34(2)：604-606.

[9] AIAZZI B，ALPARONE L，BARONTI S. Information-theoretic heterogeneity measurement for SAR imagery [J]. IEEETransaction on Geoscience and Remote Sensing(S0196-2892)，2004，43(3)：619-624.

[10] ARGENTI F，ALPARONE L. Speckle removal from SAR images in the undecimated wavelet domain [J]. IEEE Transaction on Geoscience and Remote Sensing(S0196-2892)，2002，40(11)：2363-2374.

[11] LI Hengchao，HONG Wen，WU Yiyong，et al. Bayesian wavelet shrinkage with heterogeneity-adaptive threshold for SAR image despeckling based on generalized gamma distribution [J]. IEEE Transaction on Geoscience and Remote Sensing (S0196-2892)，2013，51(4)：2388-2402.

[12] BARNDORFF-NIELSEN O E. Normal inverse Gaussian distributions and stochastic volatility modeling [J]. Scandinavian Journal of Statistics(S0303-6898)，1997，24(1)：1-13.

[13] BHUIYAN M I H，AHMAD M O，SWAMY M N S. Wavelet-based image denoising with the normal inverse Gaussian prior and linear MMSE estimator [J]. IET Image Processing(S1751-9659)，2008，2(4)：203-217.

[14] HYV?RINEN A. Sparse code shrinkage：Denoising of nongaussian data by maximum likelihood estimation [J]. Neural Computation(S0899-7667)，1999，11(7)：1739-1768.

[15] ARGENTI F，BIANCHI T，ALPARONE L. Multiresolution MAP despeckling of SAR images based on locally adaptive generalized Gaussian pdf modeling [J]. IEEE Transaction on Image Processing(S1057-7149)，2006，15(11)：3385-3399.

[16] NICOLAS J M. Introduction to second kind statistic：Application of log-moments and log-cumulants to SAR image law analysis [J]. Traitement Du Signal(S0765-0019)，2002，19(3)：139-167.

[17] BHUIYAN M I H，AHMAD M O，SWAMY M N S. Wavelet-based despeckling of medical ultrasound images with the symmetric normal inverse Gaussian prior [C]// IEEE International Conference on Acoustics，Speech and Signal Processing，Acoustic，F(xiàn)rance，April 15-20，2007，1：721-724.

[18] FADILI J M，BOUBCHIR L. Analytical form for a Bayesian wavelet estimator of images using the Bessel K form densities [J]. IEEE Transaction on Image Processing(S1057-7149)，2005，14(2)：231-240.

[19] CHANG S G，YU Bin，VETTERLI M. Adaptive wavelet thresholding for image denoising and compression [J]. IEEE Transaction on Image Processing(S1057-7149)，2000，9(9)：1532-1546.

Wavelet SAR Image Despeckling Based on Heterogeneity Classification

HOU Jianhua，CHEN Wen，LIU Xinda，CHEN Shaobo
( Hubei Key Laboratory of Intelligent Wireless Communications, South-Central University for Nationalities, Wuhan 430074, China )

Abstract:A Bayesian wavelet speckle reduction algorithm for SAR image is developed under the non-homomorphic framework. We use Normal Inverse Gaussian (NIG) function for modeling backscattered signal in wavelet domain, and Gaussian function for speckle noise (i.e. signal-dependent noise). The estimation formula of noise-free signal is derived by Bayesian maximum a posteriori (MAP) criterion. With regarding to estimation of model parameters, we introduce Multiscale Local Coefficient of Variation (MLCV) as heterogeneity measure, the histogram of which can be well fitted by logarithmic normal distribution. Based on heterogeneity measure, each coefficient in wavelet sub-band is classified into one of several different heterogeneity scenes, and NIG model parameters are computed in each class through cumulants estimation method. Experiment results show that, compared with its counterpart algorithm in homomorphic framework and its counterpart algorithm in non-homomorphic framework without heterogeneity based classification, our method has obvious advantage in terms of both subjective and objective evaluation, and has obtained satisfactory de-speckled image. A classification method of wavelet coefficients is proposed by heterogeneity measure, which could provide a new means for the research of SAR image despeckling.

Key words:SAR images despeckleing; Normal inverse Gaussian distributiom; heterogeneity measure; classification

作者簡介：侯建華(1964-)，男(漢族)，湖北武漢人。博士，教授，主要從事圖像處理和計算機視覺方面的研究。E-mail：hou878l@126.com。

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

收稿日期：2015-03-27； 收到修改稿日期：2015-06-07

文章編號：1003-501X(2016)02-0055-07

中圖分類號：TP391

文獻標志碼：A

doi：10.3969/j.issn.1003-501X.2016.02.010