李光廷 禹衛(wèi)東

①(中國科學院電子學研究所 北京 100190)

②(中國科學院研究生院 北京 100049)

1 引言

SAR圖像相干斑的存在給SAR圖像的理解與解譯帶來了極大的困難，相干斑抑制研究也一直是SAR圖像處理領域的重要課題。研究SAR圖像相干斑抑制的目的是在平滑斑點噪聲的同時較好地保持圖像的細節(jié)信息，至今已有許多經典算法取得了較好的濾波效果[1]。在實際中常用的濾波算法有Lee濾波，Kuan濾波，Frost濾波，Gamma MAP濾波及其它們的改進算法等的空域濾波，還有以小波變換為主的變換域濾波。空域濾波通常對圖像中的同質區(qū)域有著較好的去斑效果，而邊緣等細節(jié)區(qū)域的保持則不夠理想。在小波域濾波中，多種具有良好2維線奇異性刻畫能力的二代小波備受關注，但小波域濾波普遍存在偽吉布斯效應明顯、閾值選擇困難等問題，其中一類非下采樣 Contourlet變換(NonSubsampled Contourlet Transform,NSCT)近幾年比較受重視，由于它具有平移不變性可以大大降低偽吉布斯效應在相干斑抑制方面已經得到了應用[2,3]。由上面分析可見，相干斑抑制研究的關鍵在于如何充分利用像元周圍的結構信息，從而在平滑噪聲的同時來保持圖像細節(jié)。

本文研究的 Bilateral濾波器[4]是由 Tomasi和Manducci從直覺意義上提出一類空域濾波器，自提出以來許多學者對它進行了理論研究，完善了它的數學基礎[5,6]。由于Bilateral濾波器在濾波過程中同時利用中心像素與鄰域像素的空間相似性與灰度相似性來計算鄰域像素的權值，可以在平滑噪聲的同時較好地保持圖像細節(jié)，因此在圖像處理方面得到了越來越廣泛應用[7-9]。近年來，該算法被引入到了SAR圖像相干斑抑制的研究中[10]。由于Bilateral濾波器在計算圖像空間相似性與灰度相似性時均采用高斯函數，而高斯函數的方差系數選擇比較困難，文獻[9,10]對這一問題進行了研究。文獻[9]發(fā)現最優(yōu)的距離方差系數與圖像的噪聲標準差呈正比，而最優(yōu)灰度方差系數對圖像噪聲并不敏感；文獻[10]提出了一種利用等效視數及邊緣保持指數分兩步進行兩個方差系數選擇的方法。文獻[10]的方法取得了較好的相干斑抑制效果，但在進行方差系數選擇時計算量大，且最終選用恒定的方差系數。

由于相干斑的乘性特性，分析發(fā)現進行SAR圖像濾波時兩個方差系數設為定值是不合理的。針對這個問題，本文提出了一種自適應Bilateral濾波器；在進行空間相似性度量時，根據圖像的局域變差系數來選擇局域方差系數；在進行灰度相似性度量時，利用SAR圖像的似然概率函數作為度量模型；分別從兩方面解決了兩個方差系數選擇困難的問題。

2 Bilateral濾波器

Bilateral濾波器通過對窗口內不同像素設定不同的權重來進行中心像素值的估計，不同像素的權重由像素間的空間距離與灰度距離共同決定。給定輸入圖像f，則輸出h為

式中分母為歸一化系數。其中x=(x1,x2)代表中心像素的坐標，ξ=(ξ1,ξ2)是像素x的一個相鄰像素，c(ξ,x)和s(f(ξ),f(x))為兩像素的空間相似性與灰度相似性，兩個常用的相似性度量模型均為高斯函數

σd,σr分別為空間方差系數與灰度方差系數。

方差系數是高斯函數的形狀參數。對空間相似性來說，由式(2)可知，σd越大，則c(ξ,x)隨著空間距離增加而衰減的速度越慢，相反則越快。傳統(tǒng)的Bilateral算法在濾波過程中采用恒定的σd，這存在如下問題：若σd較大，則鄰域像素權值較大容易帶來細節(jié)的模糊；若σd較小，則鄰域像素的權值較小使得同質區(qū)域的去噪能力又太差。因此，最優(yōu)σd的選擇與區(qū)域的同質性有關，對同質性較強的區(qū)域應該選擇較大的σd來增強算法的噪聲平滑能力，對邊緣等同質性較弱的區(qū)域則應該選擇較小的σd來進行細節(jié)的保持。

對SAR圖像的灰度相似性來說，用高斯函數作為度量模型是不合理的。具體分析如下：多視SAR幅度圖像的相干斑通常用均值為1且服從均方根伽瑪分布的乘性噪聲來建模[11]，設g為真實灰度值，n為噪聲，則f=gn。對L視幅度SAR圖像，噪聲的概率模型為

用高斯函數度量兩像素的灰度相似性時，相似性僅由像素的灰度差決定。以灰度差為50的兩組像素來說明，由式(3)得s(60,10)=s(60,110)，而由式(4)得p(60/10)?p(60/110)，因此，用高斯函數來度量SAR圖像的灰度相似性是不合理的。分析還發(fā)現，若通過對數變換將乘性噪聲變?yōu)榧有栽肼?，由于變換后的噪聲仍然是非高斯的，所以用高斯函數來度量對數變換后 SAR圖像的灰度相似性也是不精確的。

綜上可得，當直接用Bilateral濾波器進行SAR圖像相干斑抑制時存在空間方差系數選擇困難和灰度相似性模型不合理的問題。針對這兩個問題，本文提出了一種自適應Bilateral濾波算法。

3 自適應Bilateral濾波器

3.1 基于SAR圖像變差系數的空間方差系數選擇

如上一節(jié)分析，用 Bilateral濾波器進行 SAR圖像相干斑抑制時，空間方差系數σd的選擇應該與區(qū)域的同質性有關。因此，我們可以建立σd與區(qū)域同質性的對應關系，從而根據區(qū)域的同質性來選擇σd。變差系數是 SAR圖像同質性測量最常用的指標，本文仍然用變差系數來進行區(qū)域的同質性測量。變差系數CV的計算公式為

式中var(f(x)),E[f(x)]是窗口內像素的方差和均值。

Lopes等人[12]根據兩個標準差Cu與Cmax將圖像分為3類區(qū)域：第1類滿足CV≤Cu，為均勻區(qū)域，是濾波過程中需要平滑的區(qū)域；第2類滿足Cu＜CV＜Cmax，為含有一定紋理信息的區(qū)域，是過度區(qū)域；第3類滿足CV≥Cmax，為非均勻區(qū)域，大多為圖像的細節(jié)，是需要保持的區(qū)域。其中

由對方差系數的分析知，方差系數σd與變差系數CV應成反比關系，本文所選的模型為

式中A,kd與Cd為待定參數，以濾波窗長N=7為例，參數確定過程如圖1所示。

由式(2)知，c(ξ,x)隨著兩像素空間距離的增加而減小，以0.5作為c(ξ,x)的閾值，3個未知參數按如下思路來確定：對第1類區(qū)域，當/2時，c(ξ,x)≥0.5，并且當CV=Cu時等號成立，如圖1(a),1(b)中菱形所注；對第3類區(qū)域，當時，c(ξ,x)≤0.5，并且當CV=Cmax時等號成立，如圖1(a),1(b)中方形所注。將上述等式成立時的條件代入式(2)得

圖1 模型參數的確定

解式(7)可得σd(Cu)與σd(Cmax)。結合式(6)所示模型的對稱性得

將σd(Cu)與σd(Cmax)代入式(6)并結合式(8)可求解3個未知參數，從而便確定了式(6)。

3.2 基于SAR圖像似然概率函數的灰度相似性度量

如第2節(jié)分析，用高斯函數作為SAR圖像的灰度相似性模型是不合理的。根據式(4)及相干斑的乘性模型可以得到f對g的似然概率函數[13]。本文用該似然函數進行鄰域像素與中心像素的灰度相似性度量，由此得到新的灰度相似性模型為

對于L視幅度SAR圖像，分析式(9)可以發(fā)現，s(f(ξ),f(x))由f(x)和f(ξ)/f(x)共同決定。但對當前像素x來說，f(x)是定值，式(1)中分母的歸一化操作消除了f(x)對s(f(ξ),f(x))的影響，如圖2(a)所示。因此，s(f(ξ),f(x))只與f(ξ)/f(x)有關，這與相干斑的乘性性質是一致的。

為了說明用似然函數來度量SAR圖像灰度相似性的合理性，仿真產生了灰度值為50和150的兩幅4視幅度SAR圖像，并分別對它們的概率密度進行了統(tǒng)計，如圖2(b)中‘3-﹡’所示；分別用似然函數與高斯函數(σr=60)來度量圖像的灰度相似性，如圖2(b)中‘1-﹡,2-﹡’所示。由圖可以看出，用似然函數度量的圖像灰度相似性更接近仿真圖像的概率密度分布，因此用似然函數作為SAR圖像的灰度相似性模型更合理。

3.3 自適應Bilateral濾波器的實現

根據前面的分析，自適應Bilateral濾波器的實現步驟可總結為

步驟1 輸入圖像，設置窗長N。

步驟2 若圖像視數L未知，計算等效視數ENL代替視數L。

圖2 似然函數進行灰度相似性度量結果

步驟3 計算局域變差系數CV，結合式(6)~式(8)確定σd與CV對應關系。

步驟4 執(zhí)行濾波，逐像素窗口操作，得到像素估計值。

(1)根據式(6)，由CV(x)確定σd(x)；

(2)將σd(x)代入式(2)求得各像素空間相似性，由式(9)求得像素灰度相似性；

(3)將式(2)，式(9)計算結果代入式(1)，求得像素估計值。

步驟5 將步驟2~步驟4重復3~5次，得到最終濾波結果。

4 實驗與結果分析

4.1 評價指標

在實驗結果的評價方面，這里用以下兩個指標進行評價：

(1)等效視數ENL反映圖像的噪聲平滑能力，在同質區(qū)域中，ENL越大，則相干斑抑制效果越好。設所選圖像的同質區(qū)域為Ihom，則ENL計算公式為[10]

(2)在細節(jié)保持評價方面，這里引入一新指標：細節(jié)保持系數(Detail-Preservation Index,DPI)，計算如下：

其中D為圖像的第3類區(qū)域，即滿足CV≥Cmax的區(qū)域，f與h分別為濾波前后圖像。DPI_M越接近1,DPI_V越小，則細節(jié)保持越好。

4.2 實驗結果

為了驗證提出方法的有效性，實驗用3幅SAR圖像作為測試圖像：(1)圖3(a)，仿真4視SAR圖像；(2)圖4(a)是Volgograd地區(qū)一幅6視TerraSAR-X幅度圖像；(3)圖5(a)，Flevol and地區(qū)的一幅4視AirSAR幅度圖像，3個圖像大小均為256×256。在實驗中，本文方法分別與精致Lee濾波(RE_Lee)，文獻[10]中的基本Bilateral濾波(BF)，基于NSCT局域高斯模型的MAP濾波[2](NSCT_GMAP)和基于NSCT雙變量模型的濾波算法[3](SNSCTBI)進行了比較。在參數設置方面，本文選用5×5窗長迭代5次，比較算法各參數設置為：RE_Lee選用7×7的窗，BF選用11×11的窗，NSCT_GMAP與SNSCTBI均采用‘cd 9-7’小波5層分解，高頻方向數均為[4 8 8 16]。各方法的濾波結果如圖3~圖5所示，濾波效果的評價如表1所示，各圖(a)中矩形框區(qū)域為選定用于進行ENL計算的區(qū)域。

圖3 仿真SAR圖像濾波前后比較

觀察圖3中各方法的濾波結果可以看出：RE_Lee濾波相比BF的噪聲平滑能力更強，但同質區(qū)域呈現出了具有一定方向性的斑塊，這是由濾波過程中的方向計算不準確引起的；NSCT_GMAP濾波與SNSCTBI濾波相比RE_Lee濾波具有更強的噪聲平滑能力，但邊緣細節(jié)出現了模糊，且在同質區(qū)域呈現了明顯的偽吉布斯效應；本文方法的同質區(qū)最平滑，邊界清晰。表1對各濾波結果與真值間的均方誤差(MMSE)進行了比較，結果顯示本文方法的處理結果具有最小的均方誤差值，這說明本文方法的濾波結果更接近于加噪聲前的真實圖像。

觀察圖4中各方法的濾波結果，尤其注意圖像中部的一系列點目標，本文方法的濾波結果中目標清晰，而其它方法均有一定程度的模糊，這是由于本文方法具有自動調節(jié)窗口大小的能力，在變差系數大的區(qū)域窗長很小，可以起到細節(jié)保持的作用；觀察圖5的濾波結果可以發(fā)現，本文方法的濾波結果在接近區(qū)域交界線的地方，仍然具有很強的噪聲平滑能力，這里因為本文方法可以根據像素的相關性對窗口內不同區(qū)域像素賦不同的權重，因此在接近區(qū)域交界線的地方仍然有盡可能多的同質區(qū)域像素起到噪聲平滑的作用。

實驗發(fā)現本文方法的濾波結果會存在一定的暗點或者暗斑，這是由于以當前像素作為真值來計算灰度相似性引起的，本文對濾波結果用簡單的排序濾波器去最小值解決了該問題。

比較表1中的各項指標我們可以發(fā)現：相比其它方法，本文方法的濾波結果在同質區(qū)域具有最大的等效視數，且衡量細節(jié)保持的DPI_M更接近于1,DPI_V更接近于0；這一結果也證明本文方法在噪聲平滑與細節(jié)保持上均優(yōu)于其它同類算法。

圖4 Volgograd地區(qū)SAR圖像濾波結果

圖5 Flevoland地區(qū)SAR圖像濾波結果

5 結論

SAR圖像的相干斑抑制一直是SAR圖像處理方面的研究熱點。由于Bilateral濾波器在進行圖像去噪時可以同時利用像素間的空間相似性與灰度相似性，近年來被引入到了SAR圖像相干斑抑制中。分析發(fā)現，用Bilateral濾波器進行相干斑抑制時存在空間方差系數選擇困難與灰度相似性模型不合理的問題，針對這兩個問題，本文提出了一種自適應Bilateral濾波器。該方法根據SAR圖像的局域變差系數來確定空間方差系數，從而調整高斯窗的形狀，起到自適應調節(jié)窗長的作用；同時根據SAR圖像的似然概率函數來度量窗內各像素與中心像素的灰度相似性，并根據相似性調節(jié)各像素的權重，起到方向選擇的作用；從這個意義上講，本文方法能夠根據圖像的結構信息自動調節(jié)濾波器的窗長與方向，是一種具有自適應能力的 SAR圖像相干斑抑制算法。仿真與實測SAR圖像的去斑實驗表明，自適應Bilateral濾波器在去噪能力與細節(jié)保持方面均優(yōu)于其它同類算法，是一種有效的相干斑抑制方法。

表1 濾波效果比較

[1]Amirmazlaghani M and Amindavar H.Two novel Bayesian multiscale approaches for speckle suppression in SAR images[J].IEEE Transactions on Geosciences and Remote Sensing,2010,48(7): 2980-2993.

[2]鳳宏曉,侯彪,焦李成,等.基于非下采樣 Contourlet域局部高斯模型和MAP的SAR圖像相干斑抑制[J].電子學報,2010,38(4): 811-816.Feng Hong-xiao,Hou Biao,Jiao Li-cheng,et al..SAR image despeckling based on local Gaussian model and MAP in NSCT domain[J].Acta Electronica Sinica,2010,38(4):811-816.

[3]賈建,陳莉.基于雙變量模型和非下采樣 Contourlet變換的SAR圖像相干斑抑制[J].電子與信息學報,2011,33(5):1088-1094.Jia Jian and Chen Li.SAR image despeckling based on bivariate threshold function in NSCT domain[J].Journal of Electronics＆Information Technology,2011,33(5):1088-1094.

[4]Tomasi C and Manduchi R.Bilateral filtering for grey and color images[C].Proceedings of IEEE International Conference on Computer Vision,New Delhi,India,1998:839-846.

[5]Barash D.A fundamental relationship between bilateral filtering,adaptive smoothing,and the nonlinear diffusion equation[J].IEEE Transactions on Pattern Analysis and Machine Intelligent,2002,24(6): 844-847.

[6]Elad M.On the origin of the bilateral filter and ways to improve it[J].IEEE Transactions on Image Processing,2002,11(10): 1141-1151.

[7]Han J W,Kim J H,Cheon S H,et al..A novel image interpolation method using the Bilateral filter[J].IEEE Transactions on Consumer Electronics,2010,56(1): 175-181.

[8]Lin C H,Tsai J S,and Chiu C T.Switching Bilateral filter with a texture/noise detector for universal noise removal[J].IEEE Transactions on Image Processing,2010,19(9):2307-2320.

[9]Zhang M and Gunturk B K.Multiresolution bilateral filtering for image denoising[J].IEEETransactions on Image Processing,2008,17(12): 2324-2333.

[10]Zhang W G,Liu F,and Jiao L C.SAR image despeckling via bilateral filtering[J].Electronics Letters,2009,45(15):781-783.

[11]Oliver C and Quegan S.Understanding Synthetic Aperture Radar Images[M].Boston.Artech House,1998: 75-120.

[12]Lopes A,Touzi R,and Nezry E.Adaptive speckle filters and scene heterogeneity[J].IEEE Transactions on Geosciences and Remote Sensing,1990,28(6): 992-1000.

[13]Molina D E,Gleich D,and Datcu M.Gibbs random field models for model-based despeckling of SAR images[J].IEEE Geoscience and Remote Sensing Letters,2010,7(1): 73-77.