陳少波，侯建華，張 華

中南民族大學(xué) 電信學(xué)院，武漢 430074

改進的P-M擴散方程在相干斑抑制中的應(yīng)用

陳少波，侯建華，張 華

中南民族大學(xué) 電信學(xué)院，武漢 430074

相干斑抑制是合成孔徑雷達圖像處理的一個非常重要的而又未完全解決的環(huán)節(jié)。目前主要有兩種處理手段：多視處理和濾波后處理。多視處理是以降低圖像的分辨率為代價的，因此人們更多地去研究濾波后處理技術(shù)。在過去的20年里，斑點噪聲的后處理抑制算法得到了深入研究，涌現(xiàn)出許多較為成熟的算法。這些算法主要分兩類：（1）基于局部統(tǒng)計特性的空域濾波[1-5]：Lee、Frost、Kuan和MAP，這類濾波算法大都基于均勻目標場景而提出來的；它們利用局部統(tǒng)計參數(shù)，能夠有效地濾除斑點噪聲。（2）基于多分辨率分析的算法，其典型思路為通過小波變換這種多尺度分解的方式將SAR圖像變換到小波域，通過對小波域的系數(shù)進行重新調(diào)整來達到去斑的目的[6]?；诜蔷€性擴散方程提出的濾波理論近年來在光學(xué)圖像去噪方面取得了很好的應(yīng)用。Yu Yongjian和Santiago Aja-Fernandez等人將各向異性擴散方程成功應(yīng)用于相干斑抑制：SRAD[7]和DPAD[8]。這類算法的思路是：基于局域統(tǒng)計特性的自適應(yīng)濾波器可以轉(zhuǎn)化成線性擴散模型，于是他們將局部統(tǒng)計特性引入擴散系數(shù)中，得到了相應(yīng)的去斑模型。SRAD方法與傳統(tǒng)的空間濾波相比有更加好的相干斑抑制效果，但是也存在不足之處：該方法對目標的判斷不準確，在迭代次數(shù)較多的情況下，會導(dǎo)致圖像嚴重模糊。Yu Yongjian曾經(jīng)在文獻[7]中指出直接將加性噪聲抑制的各向異性擴散方程用來抑制相干斑噪聲時，不但不能抑制噪聲，反而會增強相干斑；但可以先通過對圖像取對數(shù)，將相干斑噪聲轉(zhuǎn)換為加性噪聲，然后選取適合加性噪聲抑制的各向異性擴散方程來對轉(zhuǎn)化后的圖像進行擴散處理；目前關(guān)于各向異性擴散方程對加性噪聲的抑制算法的研究比較深入，因此思路將會更加寬廣一點。

本文首先分析了對數(shù)變化對于相干斑噪聲的影響，說明了適合去除加性噪聲的P-M擴散方程應(yīng)用到相干斑噪聲抑制的可行性。其次，為了使P-M擴散方程在有效去除經(jīng)過對數(shù)變換的相干斑噪聲的同時保持SAR圖像的邊緣特征，對P-M擴散方程進行了分析改進：（1）在分析穩(wěn)健統(tǒng)計學(xué)和各向異性擴散的聯(lián)系的基礎(chǔ)之上[9]，選擇了性能穩(wěn)定的擴散系數(shù)函數(shù)；（2）通過分析擴散方程的實際處理情況，對擴散過程中的關(guān)鍵參數(shù)梯度閾值進行了修改，使之更加符合實際的擴散過程，更加能保持圖像的邊緣特征[10]；（3）通過引進擴散停止判斷準則，使處理過程既準確又高效。最后給出了本文所提出的相干斑噪聲抑制算法具體的步驟，給出實驗結(jié)果，并對實驗結(jié)果進行了分析。

1 對數(shù)變換對相干斑噪聲的影響

已知N視SAR圖像的強度服從具有2N個自由度的Gamma分布[11]。

式中，I0為均值。Γ(N)為Gamma函數(shù)：

取對數(shù)變化，D=-ln(I)，則對數(shù)變換后相干斑噪聲模型的概率密度為：

式中，D0=-ln(I0)，D0僅改變分布的位置而并不影響其形狀。

為了簡化表示，可設(shè)D0=0，則

當N=1時，式（4）被稱為Fisher-Tippett分布，而式（1）退化為指數(shù)分布。

對式（4）中的exp(·)項泰勒展開，取前三項，得到：

式（6）表明，在獨立視數(shù)足夠大時（通常N＞3就可得到較滿意的近似），對數(shù)變化相干斑噪聲近似為高斯白噪聲，其方差是1/N。此時，式（1）也趨向高斯分布，但較式（6）慢很多。因此，把基于加性高斯噪聲的圖像濾波算法用于對數(shù)變化后的SAR圖像比直接用SAR圖像更加合適。

2 Perona-Malik模型分析

2.1 P-M擴散模型簡介

P-M[12]模型在熱傳導(dǎo)方程的基礎(chǔ)上加入了一個根據(jù)圖像的位置控制擴散速度的函數(shù)C(x，y，t)，得到了以下的各向異性擴散模型：

式（7）中：

最終，P-M模型的完整表達式為：

同時，Perona和Malik給出了g(s)的兩種形式：

2.2 擴散系數(shù)特性分析

按照P-M擴散方程的設(shè)計思想：擴散系數(shù)在光滑區(qū)域取值很大，而在邊緣區(qū)域則取值很小，這樣光滑區(qū)域的點就會因周圍點的平均作用而抑制掉其中所存在的隨機性噪聲，而邊緣區(qū)域的點則不會被其周圍的點平均掉，所以邊緣區(qū)域的灰度值基本不變，達到在去噪的同時保持邊緣信息的效果。由于Perona和Malik提出的擴散系數(shù)的擴散速度很大，在梯度閾值附近仍然存在擴散，故由Perona和Malik給出的擴散系數(shù)很可能在擴散的過程中消除一些不太明顯而又重要的細節(jié)邊緣和弱化高對比度的區(qū)域。Black M J在1998年給出了各向異性擴散和穩(wěn)健統(tǒng)計學(xué)的聯(lián)系、穩(wěn)健估計模型的推導(dǎo)過程和具體含義，并給出了一些具有穩(wěn)健統(tǒng)計特性的邊緣終止函數(shù)，其中比較有代表性的是Tukey提出的Biweigh函數(shù)（其中k為梯度閾值）：

下面將基于Tukey的Biweight擴散系數(shù)和Perona-Malik擴散系數(shù)（以式（7）為例）與函數(shù)圖像相結(jié)合，分析和比較一下這兩個擴散系數(shù)及其各自對應(yīng)的誤差模型和影響函數(shù)(φ(·))的特性。

從圖1和圖2中各自的邊緣終止函數(shù)和它們對應(yīng)的影響函數(shù)可以看到，基于Tukey函數(shù)的影響函數(shù)比Perona-Malik影響函數(shù)遞減程度更快，且當梯度超過某一閾值后，影響函數(shù)遞減為0，即模型不再進行擴散，從而避免了邊緣弱化的問題，因此使用基于Tukey函數(shù)的邊緣終止函數(shù)能達到更好地保護圖像邊緣和細節(jié)的效果。

圖1 Perona-Malik邊緣終止函數(shù)、對應(yīng)的影響函數(shù)和穩(wěn)健統(tǒng)計學(xué)誤差

圖2 基于Tukey函數(shù)的邊緣終止函數(shù)、對應(yīng)的影響函數(shù)和穩(wěn)健統(tǒng)計學(xué)誤差模型

2.3 梯度閾值的非線性衰減

在用P-M擴散方程對帶噪圖像作去噪處理的過程中，梯度閾值k對于圖像中邊緣保持有著非常關(guān)鍵的作用。從擴散過程的本質(zhì)來看，隨著擴散的不斷進行，圖像梯度的絕對值會越來越小。隨著梯度絕對值的減小，梯度閾值k也應(yīng)該相應(yīng)減小，從而保持圖像的邊緣不被丟失，因為邊緣保持的條件是| ?I|＞k。

針對這種情況，在仿真中使梯度閾值k隨時間變小，即擴散系數(shù)變?yōu)槿缦滦问剑?/p>

考慮擴散過程中前后圖像的關(guān)聯(lián)性和實現(xiàn)過程中的方便性，給出如下基于遞歸形式的梯度衰減形式：

其中可由γ來控制擴散的速度，隨著γ的不斷變大，退化過程越來越慢。分析上式，可以知道，隨著t的增加（即擴散次數(shù)的增加），梯度閾值k會不斷減小，達到了預(yù)期的效果。

2.4 擴散最優(yōu)停止時間

當P-M擴散方程應(yīng)用于圖像平滑的時候，由于模型中沒有考慮處理后圖像與原始圖像的逼近項，因此利用模型抑制就必須考慮最優(yōu)停止時間的問題。本文采用相對峰值信噪比來控制擴散的程度，其定義如下：

式中I(k)和I(k+1)分別是圖像經(jīng)過第k次和第k+1次迭代濾波后的結(jié)果。于是當?shù)兄箺l件

滿足時迭代過程中止。

3 新的斑點噪聲抑制算法

大量的仿真實驗證明，P-M模型對高斯白噪聲的處理效果是比較好的。根據(jù)SAR圖像的噪聲特點，結(jié)合前面的關(guān)于相干斑噪聲及其經(jīng)過對數(shù)變換后的統(tǒng)計特性的一些結(jié)論，本文選擇P-M擴散方程來去除SAR圖像的相干斑噪聲。綜合穩(wěn)健統(tǒng)計學(xué)知識和梯度閾值的非線性退化的思想，本文對P-M擴散方程進行了改進，建立了如下的擴散系數(shù)：

其中，k(t)由式（14）給出。

本文提出的斑點噪聲的濾波算法步驟描述如下：

步驟1對原始SAR圖像的加性噪聲進行抑制：通常原始SAR圖像中除了包含乘性的相干斑噪聲外，還包含有加性噪聲，亦即SAR圖像的數(shù)學(xué)模型為Y=FX+n。因此，在抑制乘性的相干斑噪聲之前，首先采用高斯濾波器來對加性噪聲n進行抑制。

步驟2對SAR圖像作對數(shù)變換：利用對數(shù)變換= 20lg(Y)將乘性噪聲模型轉(zhuǎn)換為加性噪聲模型。

步驟3利用改進的P-M方程（18）對對數(shù)域圖像的加性噪聲F~進行抑制，本步驟又分為如下幾步：

其中，g(?I，t)如式（17）。

（1）對對數(shù)域圖像作如下處理：，這樣中的加性噪聲′將在理論上趨近于高斯白噪聲。

圖3 仿真實驗1濾波效果

圖4 仿真實驗2濾波效果

（2）利用改進的P-M模型對進行噪聲抑制：在處理之前要確定高斯濾波函數(shù)的參數(shù)σ，擴散模型的迭代次數(shù)N，步長Dt，梯度閾值參數(shù)的初始值k(1)以及其衰減系數(shù)λ；利用式（9）對Y~進行噪聲抑制，其離散化形式如式（19）所示，利用式（15）和式（16）所表示的擴散停止準則來確定擴散的停止時間。

其中N，S，E，W是North，South，East和West的簡寫，τ為時間步長，符號?并不是梯度算子，而是表示差分，即

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

4.1 實驗結(jié)果評價指數(shù)

國外公共數(shù)據(jù)庫MSTAR數(shù)據(jù)是公開評價SAR算法性能的標準數(shù)據(jù)。下面通過對MSTAR數(shù)據(jù)的實驗來驗證本文方法的相干斑抑制效果。實驗中，分別采用Lee濾波、SRAD濾波和本文方法進行處理。

在對處理結(jié)果進行定量評價時，采用圖像均值μ、等效視數(shù)（ENL）和邊界保持指數(shù)（EPI）幾個評價指標。

均值是衡量圖像整體特征的指標，反映圖像的平均亮度，處理前后要求均值要基本保持。

等效視數(shù)是衡量相干斑的相對強度的一種，定義為：ENL=μ2/σ2，其中μ和σ2分別為圖像均勻區(qū)域的均值和方差。等效視數(shù)越高，表明相干斑抑制越好。

圖像的邊緣保持指數(shù)是衡量算法對圖像邊緣保持程度的重要指標，其定義為：ps(i，j)是濾波后圖像像元，pο(i，j)是原始圖像像元，ps(i，j)和pο(i，j)都位于邊緣區(qū)域，i是行數(shù)，j是列數(shù)。EPI的取值范圍為：0～1，EPI的值為1時圖像的邊緣得到完全保持，EPI的值為0時，圖像已經(jīng)變?yōu)橐粋€平面，沒有任何變化。

4.2 實驗結(jié)果

仿真結(jié)果見圖3和圖4；實驗數(shù)據(jù)見表1和表2。從仿真結(jié)果和實驗數(shù)據(jù)可以看出，Lee和SRAD在一定程度上降低了圖像的標準差，且等效視數(shù)較大；但是在保留邊緣細節(jié)方面效果不是很好；從視覺效果來說平滑效果也不理想。相比之下，本文提出的算法不僅可以較好地保持均值，大大降低標準差，對均勻區(qū)域的平滑效果較好，而且在保持和增強邊緣細節(jié)方面效果也很明顯，所得到的SAR圖像有良好的視覺效果。

表1 實驗結(jié)果1

表2 實驗結(jié)果2

4.3 參數(shù)分析

從前面的討論中可以知道，參數(shù)γ和ε的選取對模型的擴散行為有著重要的影響。由圖5知，隨著參數(shù)γ的變大，梯度閾值的非線性衰減越厲害，這樣就會對圖像的邊緣細節(jié)保護比較好，EPI會有所升高；但是由于梯度閾值的減小，很多噪聲點也被保留下來了，因此ENL會有所下降。由圖6知，隨著ε的變小，擴散停止時間會增加，這樣圖像會越來越光滑，ENL有所提高，而EPI值會有所下降。因此在實際的處理過程中要根據(jù)實際的應(yīng)用環(huán)境來選擇參數(shù)γ和ε，以達到更佳的應(yīng)用效果。

圖5 ε和γ對ENL的影響曲線示意圖

圖6 ε和γ對EPI的影響曲線示意圖

5 總結(jié)

本文提出的相干斑噪聲抑制算法是建立在P-M[11]擴散模型和對數(shù)變化對相干斑噪聲的理論推導(dǎo)之上的。通過透徹分析P-M模型的擴散行為，結(jié)合實際仿真環(huán)境，對擴散系數(shù)、梯度閾值、擴散停止時間的選擇等方面進行改進，并給出了詳細的算法步驟和實驗結(jié)果。通過理論分析與實驗結(jié)果說明該算法在去相干斑和保持邊緣等方面都有一定優(yōu)勢。

[1]Lee S J.Digital image enhancement and noise filtering by use of local statistics[J].IEEE Transactions on Pattern Analysis and Machine Intelligence，1980，2（2）：165-168.

[2]Frost S V.A model for radar images and its application to adaptive digital filtering of multiplicative noise[J].IEEE Trans on Pattern Analysis and Machine Intelligence，1982，4（2）：157-166.

[3]Kuan T D.Adaptive restoration of images with speckle[J]. IEEE Trans on ASSP，1987，35（3）：373-383.

[4]Lopes，Nezry E，Touzi R，et al.Structure detection and statistical adaptive speckle filtering in SAR images[J].Int J Remote Sensing，1993，14（9）：1735-1758.

[5]Lopes，Touzi R，Nezry E.Adaptive speckle filters and scene hoterogeneity[J].IEEE Transactions on Geoscience and Remoted Sensing，1990，28（6）：992-1000.

[6]Xie H，Pierce L E，Ulaby F T.SAR speckle reduction using wavelet denoising and Markovrandom field modeling[J].IEEE Trans on Geosci Remote Sens，2002，40（10）：2196-2212.

[7]Yu Y J，Action S T.Speckle reduction anisotropic diffusion[J]. IEEE Transactionson ImageProcessing，2002，11（11）：1260-1270.

[8]Aja-Fernandez S，Alberola-Lopez C.On the estimation of the coefficient of variation for anisotropic diffusion speckle filtering[J].IEEE Trans on Image Processing，2006，15（9）：2694-2701.

[9]Black M J，Sapiro G，Marimont D H，et al.Robust anisotropic diffusion[J].IEEE Transactions on Image Processing，1998，7（3）：421-432.

[10]王毅，張良培，李平湘.多光譜遙感圖像的自適應(yīng)各向異性擴散濾波[J].遙感學(xué)報，2005，9（6）：659-666.

[11]Ulaby F，Dobson M C.Handbook of radar scattering statistics for terrain[M].Norwood，MA：Artech House，1989.

[12]Perona P，Malik J.Scale-space and edge detection using anisotropic diffusion[J].IEEE Transactions on Pattern Analysis and Machine Intelligence，1990，12（7）：629-639.

CHEN Shaobo,HOU Jianhua,ZHANG Hua

College of Electronics and Information Engineering,South-Central University for Nationalities,Wuhan 430074,China

The filter based on P-M diffusion equation directly has been proven to be efficient to the additive noise,but it often fails when facing the multiplicative noise（such as the SAR image speckle）.A novel speckle reduction method based on improved P-M diffusion equation is presented.By discussing the influence of applying log transfer to speckle,the theoretical foundation of applying P-M diffusion equation to speckle reduction is provided.A new diffusion coefficients based on Tukey’s Biweight Estimator error norm by recurring to the relationship between robust statistics and P-M diffusion is built.Simultaneously,the nonlinear time-dependent cooling technique for gradient threshold is incorporated into the new diffusion coefficients.In the experiment, the method is proven good performance in reducing speckle noise and preserving edges and details at the same time.Moreover, the filtered SAR images look good.

speckle;log-transform;Perona-Malik（P-M）diffusion equation;robust statistics;gradient threshold

直接基于Perona-Malik擴散方程的濾波算法對于加性噪聲非常有效，但是對于乘性噪聲（如合成孔徑雷達（SAR）圖像相干斑噪聲）收效甚微。提出了一種基于改進的Perona-Malik擴散方程抑制SAR圖像相干斑噪聲的新算法。分析對數(shù)變化對相干斑噪聲的影響，為將P-M擴散方程應(yīng)用于相干斑噪聲抑制奠定了理論基礎(chǔ)；通過P-M擴散和穩(wěn)健統(tǒng)計學(xué)的聯(lián)系，建立了基于Biweight Estimator誤差模型的擴散系數(shù)；同時利用非線性衰減技術(shù)對梯度閾值的選擇改進。實驗表明，該方法不僅有效抑制了SAR圖像相干斑噪聲，較好地保持了細節(jié)和邊緣信息，而且視覺效果比較好。

相干斑；對數(shù)變換；P-M擴散方程；穩(wěn)健統(tǒng)計；梯度閾值

A

TN911.73

10.3778/j.issn.1002-8331.1303-0403

CHEN Shaobo,HOU Jianhua,ZHANG Hua.Application of improved P-M diffusion equation in speckle reducing. Computer Engineering and Applications,2013,49（18）：1-5.

國家自然科學(xué)基金（No.61201448）；湖北省自然科學(xué)基金（No.2011CHB043）；武漢市科技攻關(guān)項目（No.201212521825）；中南民族大學(xué)中央高?；究蒲袠I(yè)務(wù)費專項資金項目（No.CZY10001）。

陳少波（1980—），男，博士，講師，研究方向：圖像去噪、字符識別；侯建華（1964—），男，博士，教授；張華（1972—），男，博士，講師。E-mail：chenshaobo1980@qq.com

2013-03-26

2013-05-17

1002-8331（2013）18-0001-05

CNKI出版日期：2013-06-08 http://www.cnki.net/kcms/detail/11.2127.TP.20130608.1001.029.html