陳少波，侯建華，張 華

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

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

陳少波，侯建華，張 華

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

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

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

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

已知N視SAR圖像的強(qiáng)度服從具有2N個(gè)自由度的Gamma分布[11]。

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

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

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

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

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

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

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

2 Perona-Malik模型分析

2.1 P-M擴(kuò)散模型簡(jiǎn)介

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

式（7）中：

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

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

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

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

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

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

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

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

2.3 梯度閾值的非線性衰減

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

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

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

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

2.4 擴(kuò)散最優(yōu)停止時(shí)間

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

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

滿足時(shí)迭代過程中止。

3 新的斑點(diǎn)噪聲抑制算法

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

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

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

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

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

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

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

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

圖3 仿真實(shí)驗(yàn)1濾波效果

圖4 仿真實(shí)驗(yàn)2濾波效果

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

其中N，S，E，W是North，South，East和West的簡(jiǎn)寫，τ為時(shí)間步長(zhǎng)，符號(hào)?并不是梯度算子，而是表示差分，即

4 實(shí)驗(yàn)結(jié)果與分析

4.1 實(shí)驗(yàn)結(jié)果評(píng)價(jià)指數(shù)

國(guó)外公共數(shù)據(jù)庫(kù)MSTAR數(shù)據(jù)是公開評(píng)價(jià)SAR算法性能的標(biāo)準(zhǔn)數(shù)據(jù)。下面通過對(duì)MSTAR數(shù)據(jù)的實(shí)驗(yàn)來驗(yàn)證本文方法的相干斑抑制效果。實(shí)驗(yàn)中，分別采用Lee濾波、SRAD濾波和本文方法進(jìn)行處理。

在對(duì)處理結(jié)果進(jìn)行定量評(píng)價(jià)時(shí)，采用圖像均值μ、等效視數(shù)（ENL）和邊界保持指數(shù)（EPI）幾個(gè)評(píng)價(jià)指標(biāo)。

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

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

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

4.2 實(shí)驗(yàn)結(jié)果

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

表1 實(shí)驗(yàn)結(jié)果1

表2 實(shí)驗(yàn)結(jié)果2

4.3 參數(shù)分析

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

圖5 ε和γ對(duì)ENL的影響曲線示意圖

圖6 ε和γ對(duì)EPI的影響曲線示意圖

5 總結(jié)

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

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

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

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.

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

陳少波（1980—），男，博士，講師，研究方向：圖像去噪、字符識(shí)別；侯建華（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