馬曉雙，吳鵬海

安徽大學(xué)資源與環(huán)境工程學(xué)院，安徽 合肥 230601

全極化SAR(polarimetric SAR，PolSAR)系統(tǒng)作為傳統(tǒng)單極化SAR系統(tǒng)的一種更高級(jí)形式，能同時(shí)發(fā)射、接收水平及垂直的極化雷達(dá)波，能夠獲取地物更為豐富的后向散射信息，更利于地物場(chǎng)景的解譯，已在諸多行業(yè)領(lǐng)域得到廣泛應(yīng)用。特別是，我國(guó)于2016年發(fā)射了國(guó)內(nèi)首顆民用全極化SAR衛(wèi)星——高分3號(hào)，進(jìn)一步推動(dòng)了PolSAR的應(yīng)用前景。然而，由于成像機(jī)理的限制，使得PolSAR影像不可避免地存在相干斑噪聲。相干斑嚴(yán)重降低了SAR影像的質(zhì)量，也降低了利用SAR數(shù)據(jù)進(jìn)行地物解譯的能力[1-2]。因此，對(duì)于PolSAR數(shù)據(jù)的后續(xù)應(yīng)用而言，相干斑濾波是必不可少的預(yù)處理程序。

對(duì)于全極化SAR數(shù)據(jù)而言，噪聲不僅存在于極化協(xié)方差矩陣的對(duì)角線元素上，同時(shí)也存在于非對(duì)角線元素上。并且，對(duì)角線元素受到乘性噪聲的影響，而非對(duì)角線元素則受到乘性和加性混合的噪聲的影響[3]。這使得全極化SAR的去噪問題要比單極化SAR要復(fù)雜得多。早期的全極化SAR去噪的方法大都是基于文獻(xiàn)[4]提出的極化白化濾波器。之后，文獻(xiàn)[1]提出了一種基于線性最小均方誤差(linear minimum mean square error，LMMSE)原則的PolSAR去噪算法(增強(qiáng)Lee濾波)。這一方法極大地推動(dòng)了對(duì)極化SAR去噪的研究，也引領(lǐng)了之后一系列基于LMMSE濾波器的算法[5-6]。除此之外，基于二叉樹的算法[7]、各向異性擴(kuò)散算法[8-9]、變分正則化算法[10-11]、模擬退火算法[12]等也被應(yīng)用于極化SAR去噪。近年來，基于非局部均值(nonlocal means，NLM)思想的算法已成為了PolSAR影像去噪的研究熱點(diǎn)[13-15]，并取得了較好的去噪效果。NLM算法最先應(yīng)用于數(shù)字圖像處理[16]。其基本思想是認(rèn)為影像中往往都會(huì)存在重復(fù)出現(xiàn)的或相似的影像塊，因此像素的估計(jì)值可由它的非局部搜索窗口中的所有像素加權(quán)平均后得到，而其中每一個(gè)像素的權(quán)重由它的鄰域窗口與待處理像素鄰域窗口的差異性所獲得。差異越大，權(quán)重越小。

然而，傳統(tǒng)PolSAR非局部均值去噪算法在度量像素或影像塊極化相似度時(shí)基本利用到的都是含噪聲影像的某種統(tǒng)計(jì)特性，而沒有考慮到原始像素間的差異性。并且，許多PolSAR非局部均值算法僅僅考慮的是強(qiáng)度信息相似性，而未考慮包括極化相位差在內(nèi)的全極化信息，導(dǎo)致難以獲得穩(wěn)健準(zhǔn)確的相似性測(cè)度。有鑒于此，本文提出一種迭代優(yōu)化的PolSAR NLM算法。該算法的核心思想是：在每次迭代去噪過程中，同時(shí)考慮原始影像全極化噪聲統(tǒng)計(jì)特性和前一次迭代所得的相對(duì)低噪聲度影像的全極化信息，以完善影像塊間極化相似性的度量。試驗(yàn)部分的模擬PolSAR和真實(shí)PolSAR影像去噪結(jié)果都表明，本文提出的算法取得了較好的去噪效果，噪聲得到了有效的抑制，影像邊緣細(xì)節(jié)和極化散射特性都得到了較好的保持。

1 非局部均值及PolSAR相干斑統(tǒng)計(jì)特性

1.1 非局部均值去噪思想

NLM去噪算法濾波過程的通式可表示為

(1)

式中，f為原始噪聲影像；w代表一個(gè)權(quán)重函數(shù)。NLM算法認(rèn)為目標(biāo)像素的估計(jì)值由它的非局部搜索窗口中的像素加權(quán)平均后得到，而其中每一個(gè)像素的權(quán)重由該像素的鄰域窗口與待處理像素鄰域窗口的整體差異性所計(jì)算獲得。圖1展示了NLM的基本思想：對(duì)于以像素p為中心的鄰域塊而言，由于以q1為中心的鄰域塊的內(nèi)部結(jié)構(gòu)與之非常相似，所以在對(duì)p進(jìn)行去噪時(shí)，像素q1的權(quán)重w(p,q1)較大；而由于以q2和q3為中心的鄰域塊的內(nèi)部結(jié)構(gòu)與之相差較大，所以在對(duì)p進(jìn)行去噪時(shí)，像素q2和q3的權(quán)重很小?？梢钥吹剑捎贜LM算法不僅考慮了單一像素間影像值的接近性，還考慮到了兩個(gè)像素鄰域塊的結(jié)構(gòu)布局的相似性，這使其比起傳統(tǒng)的鄰域?yàn)V波方法更能取得穩(wěn)定有效的濾波效果。顯而易見，決定NLM算法濾波效果的一個(gè)重要因素就是影像塊間相似性的度量，本文的主要?jiǎng)?chuàng)新就是提出了一種迭代優(yōu)化的全極化影像塊相似性度量方式。

圖1 NLM算法的思路Fig.1 The diagram of the idea of nonlocal means

1.2 PolSAR相干斑統(tǒng)計(jì)特性

眾所周知，對(duì)于SAR的強(qiáng)度數(shù)據(jù)而言，其噪聲是一種服從Gamma分布的乘性噪聲[17]。然而，全極化數(shù)據(jù)的噪聲模型要復(fù)雜得多。對(duì)于噪聲完全發(fā)育的視數(shù)為L(zhǎng)的PolSAR影像而言，其極化協(xié)方差矩陣CL服從如下的復(fù)Wishart分布[18]

(2)

2 迭代優(yōu)化的PolSAR NLM去噪算法

要完善影像塊間極化相似性的度量，就必須要考慮到全極化信息。本文提出的PolSAR NLM算法的創(chuàng)新性是在每次迭代去噪過程中，同時(shí)考慮原始影像全極化噪聲統(tǒng)計(jì)特性和前一次迭代所得的相對(duì)低噪聲度影像的全極化信息。因此，如何度量原始影像塊全極化相似性及迭代所得影像塊全極化相似性是本文的主要關(guān)注點(diǎn)。

2.1 迭代優(yōu)化的PolSAR影像塊相似性度量

(3)

式中，m=n為影像的視數(shù)。對(duì)以上公式取自然對(duì)數(shù)并簡(jiǎn)化常數(shù)項(xiàng)，可得到

lnQ=6ln2+ln|X|+ln|Y|-2ln|X+Y|

(4)

式中，lnQ的取值范圍是(-∞,0]。像素越相似，則lnQ越大。

隨著影像迭代去噪過程的進(jìn)行，所得影像的噪聲程度會(huì)逐漸減弱。因此，在每次迭代的過程中，都可以利用上一次迭代所得影像的信息來完善像素的相似性測(cè)度。這一迭代優(yōu)化的相似性度量思想可以表示為

(5)

(6)

接下來可利用貝葉斯定理對(duì)式(6)右邊第1項(xiàng)進(jìn)行變化。貝葉斯公式如下

(7)

式中，P(B|A)代表?xiàng)l件A發(fā)生的前提下B發(fā)生的概率。據(jù)此，有

(8)

(9)

迭代相似性測(cè)度最終表示為

(10)

對(duì)于PolSAR數(shù)據(jù)而言，Wishart距離可表示為[18]

(11)

若考慮上述Wishart距離的對(duì)稱性，對(duì)其進(jìn)行調(diào)整即可得到如下的Kullback-Leibler散度距離度量公式[20]

(12)

最終可推得，式(12)所示的相似測(cè)度值仍為負(fù)值。值越大，代表像素間越相似。需要說明的是：本文提出的迭代去噪方法的主要思路是在每次迭代的過程中完善權(quán)重的計(jì)算，而迭代時(shí)仍然將像素的原始值同權(quán)重結(jié)合進(jìn)行加權(quán)平均，所以本質(zhì)上并不屬于“二次濾波”。

值得注意的是，由于式(4)中涉及矩陣行列式的對(duì)數(shù)，為了避免奇異值的出現(xiàn)，就應(yīng)當(dāng)保證|C|≠0，也即保證滿秩性。然而，當(dāng)影像視數(shù)小于協(xié)方差矩陣的維度時(shí)(即L<3)，影像中的一些像素的極化協(xié)方差矩陣將不再是滿秩的，此時(shí)就會(huì)導(dǎo)致數(shù)學(xué)計(jì)算的問題。為了解決這一問題，本文算法加入一種對(duì)待去噪影像進(jìn)行預(yù)處理的方法[22]，保證影像所有像素的協(xié)方差矩陣達(dá)到滿秩。這一協(xié)方差矩陣的預(yù)變換步驟可以表示為

(13)

式中，i和j代表協(xié)方差矩陣中元素的位置；γ=min(L/3,1)。這一過程的實(shí)質(zhì)就是保持協(xié)方差對(duì)角元素(即各通道的強(qiáng)度值)不變，而根據(jù)影像的視數(shù)對(duì)非對(duì)角元素進(jìn)行線性變換。研究表明[22]，以上步驟處理后的影像仍然滿足復(fù)Wishart分布，也即基于復(fù)Wishart分布的似然比檢驗(yàn)法依然適用于變換后的影像。這里需要說明的是，以上預(yù)處理過程只是用于計(jì)算每個(gè)像素的權(quán)重，而在加權(quán)平均過程中使用的仍然是原始的數(shù)據(jù)。因此，不會(huì)存在改變數(shù)據(jù)極化信息的弊端。

2.2 權(quán)重及濾波參數(shù)的設(shè)置

每次迭代濾波過程中，對(duì)于每個(gè)參考影像塊x而言，窗口內(nèi)其他影像塊y的權(quán)重為

(14)

式中，N為影像塊的像素?cái)?shù)；T在每次迭代時(shí)為一個(gè)固定參數(shù)，本文采用文獻(xiàn)[23]提出的“噪聲估算器”來確定T的取值：每一次迭代時(shí)，計(jì)算影像上所有相鄰像素的相似度，構(gòu)建由這些相似度所組成的直方圖，然后計(jì)算直方圖的90%分位數(shù)，即將此值當(dāng)作當(dāng)前時(shí)刻的參數(shù)T。最終，對(duì)于某個(gè)像素而言，其極化協(xié)方差矩陣的估計(jì)值為

(15)

式中，M為濾波窗口的大小。

本文提出的PolSAR NLM算法有3個(gè)參數(shù)需要預(yù)先設(shè)定，分別是迭代次數(shù)t、影像塊大小N和濾波窗口大小M。理論而言，隨著迭代次數(shù)增加，影像的去噪效果會(huì)越加優(yōu)化。但是，迭代次數(shù)增加也會(huì)帶來計(jì)算時(shí)間的增加。試驗(yàn)發(fā)現(xiàn)，對(duì)于大多數(shù)影像，當(dāng)?shù)螖?shù)達(dá)到3次時(shí)，就能達(dá)到優(yōu)良的去噪效果且繼續(xù)迭代后的效果提升非常有限。因此，本文提出的PolSAR NLM算法的迭代次數(shù)固定為3次。對(duì)于NLM算法而言，影像塊大小通常設(shè)置為3×3到7×7之間，而濾波窗口大小則通常設(shè)置為13×13到21×21之間。研究表明[14]，隨著影像塊和濾波窗口大小的增加，影像的噪聲會(huì)得到更明顯地抑制，但細(xì)節(jié)丟失會(huì)逐漸加重。因此，為了達(dá)到去噪與保持影像細(xì)節(jié)之間的折中，本文試驗(yàn)部分將算法的影像塊和濾波窗口大小分別設(shè)為5×5和17×17。

3 試驗(yàn)結(jié)果和分析

為了驗(yàn)證本文提出的算法的濾波效果，本文采用一幅模擬的PolSAR數(shù)據(jù)和兩幅真實(shí)的PolSAR影像進(jìn)行去噪的試驗(yàn)。作為對(duì)比，試驗(yàn)部分也展示了增強(qiáng)Lee濾波[1]、IDAN濾波[5]及Pretest NLM濾波算法[14]的去噪效果。為了達(dá)到去噪與影像細(xì)節(jié)保持的平衡，本文所有試驗(yàn)中增強(qiáng)Lee濾波的濾波窗口大小設(shè)為文獻(xiàn)[1]中所推薦的7×7；IDAN濾波的像素搜索窗口大小設(shè)置為15×15；Pretest NLM方法的影像塊和濾波窗口大小同本文方法設(shè)置一致。

3.1 模擬試驗(yàn)

模擬試驗(yàn)部分利用Monte Carlo方法[18]生成了一幅PolSAR影像。為了定量評(píng)價(jià)算法對(duì)噪聲的抑制及邊緣保持的效果，本文采用等效視數(shù)(equivalent number of looks，ENL)[24]和基于平均比的邊緣保持度(edge preservation degree based on ratio of the average,EPD-ROA)[25]這兩個(gè)指數(shù)。

ENL指數(shù)的定義為

ENL=1/σ2

(16)

式中，σ代表強(qiáng)度影像在同質(zhì)區(qū)的方差系數(shù)。ENL的估計(jì)方法一般都是人工選取影像內(nèi)幾塊大的同質(zhì)區(qū)域，利用以上公式求得每塊區(qū)域的等效視數(shù)后，取平均得到影像全局的等效視數(shù)。ENL越大，代表影像的噪聲水平越低。

EPD-ROA的定義如下

(17)

式中，m代表影像的大小；uD1(i)和uD2(i)分別表示去噪所得影像的某個(gè)方向上兩相鄰像素的強(qiáng)度值；類似地，uO1(i)和uO2(i)分別表示原始影像的某個(gè)方向上兩相鄰像素的強(qiáng)度值。EPD-ROA指數(shù)越接近于1，意味著濾波算法有著越好的邊緣保持能力。

(18)

最后對(duì)所有地物所有像素的ARB指數(shù)求均值即為參量θ最終的ARB指數(shù)。試驗(yàn)部分所采用的極化參數(shù)為Cloude極化分解參數(shù)(極化散射熵H，平均alpha角α，極化反熵A)[27]，原因是該分解模型能夠直觀地反映地物的物理散射特性。

圖2展示了幾種去噪算法的處理結(jié)果，表1為其定量評(píng)價(jià)結(jié)果。增強(qiáng)Lee濾波在較大程度上抑制了影像的噪聲。然而，仔細(xì)觀察可以發(fā)現(xiàn)，這一算法處理后影像呈現(xiàn)出一定的斑塊效應(yīng)，直接導(dǎo)致點(diǎn)目標(biāo)及邊緣受到輕度的模糊。產(chǎn)生這一問題的原因在于，該算法在濾波的過程中采用了一種選擇濾波窗口方向的策略，而非通常采用的方形窗口。與增強(qiáng)Lee濾波相比，雖然表1所示的ENL顯示出IDAN方法對(duì)噪聲的抑制能力有所提升，但其對(duì)于影像細(xì)節(jié)信息的模糊和弱化現(xiàn)象較為明顯。很顯然，與前面兩種傳統(tǒng)的局部濾波器相比，Pretest NLM濾波方法在噪聲抑制方面有了質(zhì)的提升，并且實(shí)現(xiàn)了對(duì)邊緣信息的優(yōu)良保持。然而，可以較明顯觀察到的是，Pretest NLM處理后的影像的點(diǎn)目標(biāo)模糊嚴(yán)重。相比之下，本文提出的去噪方法雖然ENL指數(shù)略微低于Pretest NLM算法，但影像的邊緣和細(xì)節(jié)信息都得到了較好的保持。就影像極化散射信息的保持性能而言，Pretest NLM方法和本文方法要明顯優(yōu)于其他兩種方法，并且這兩種方法對(duì)不同極化參數(shù)的估計(jì)精度各有優(yōu)劣。

圖2 模擬影像去噪結(jié)果Fig.2 Filtering results on the simulated image

表1 模擬影像的濾波定量評(píng)價(jià)結(jié)果Tab.1 Quantitative assessment on the simulated image

3.2 真實(shí)試驗(yàn)

除了模擬的PolSAR數(shù)據(jù)，兩幅真實(shí)的PolSAR影像也被用來驗(yàn)證去噪算法的優(yōu)劣。這兩幅影像分別為C波段和L波段的AirSAR傳感器所獲得的數(shù)據(jù)，視數(shù)為4。

圖3展示了幾種去噪算法對(duì)C波段AirSAR數(shù)據(jù)的去噪結(jié)果。可以看到，這幅影像場(chǎng)景較為復(fù)雜，既包括以表面散射為主的水體，也包含以體散射為主的森林以及散射機(jī)制較為復(fù)雜的城區(qū)。對(duì)于這類的影像，一個(gè)優(yōu)良的去噪方法應(yīng)當(dāng)在削弱噪聲的同時(shí)保持影像大部分的細(xì)節(jié)信息。很顯然，與模擬試驗(yàn)結(jié)果一樣，增強(qiáng)Lee濾波處理后的影像在同質(zhì)區(qū)存在著較明顯的斑塊效應(yīng)，部分邊緣和點(diǎn)目標(biāo)丟失；IDAN雖然更好地抑制了影像的噪聲水平，但是過平滑問題較為嚴(yán)重；Pretest NLM算法的噪聲抑制能力在4種算法當(dāng)中表現(xiàn)最佳，但也存在輕微的斑塊效應(yīng)，導(dǎo)致處理后影像城區(qū)部分一些紋理(主要是細(xì)小的道路)以及部分強(qiáng)散射點(diǎn)目標(biāo)丟失。相反，本文算法雖然噪聲去除能力略微遜于Pretest NLM，但細(xì)節(jié)保持能力有了較大的提升。值得注意的是，同Pretest NLM處理的結(jié)果相比，影像中機(jī)場(chǎng)跑道處的線狀地物(圖3(f)箭頭所指)隨著迭代的進(jìn)行反而變得更為模糊。出現(xiàn)這一問題的原因是：目視來看，該線狀目標(biāo)極化散射特性同機(jī)場(chǎng)周邊植被較為相近，因此初次濾波中大量權(quán)重略高的植被像素會(huì)參與對(duì)線狀地物的去噪。然而，初次濾波后這一目標(biāo)和植被像素的可區(qū)分性不僅沒有加大，目視來看反而變得更加接近。在這種情況下，迭代的后果就是這些像素的權(quán)重反而會(huì)提高，模糊問題不斷加深?？梢钥闯觯瑢?duì)于影像中的一些零星的細(xì)節(jié)目標(biāo)，如果本文提出的方法在初次去噪過程中不能較好地將它們與其他異質(zhì)地物“區(qū)分”開來，反而有可能會(huì)出現(xiàn)不斷過平滑的問題。因此，如何完善初始估計(jì)的精確性是本文方法未來值得改進(jìn)之處。

圖3 C波段AirSAR影像去噪結(jié)果Fig.3 Filtering results on the C band AirSAR image

圖4展示了幾種去噪算法對(duì)一景L波段PolSAR影像的去噪結(jié)果。與之前的試驗(yàn)結(jié)果相似，增強(qiáng)Lee濾波帶來的塊狀效應(yīng)制約了其濾波效果；IDAN能更好地抑制噪聲，但影像細(xì)節(jié)信息出現(xiàn)了嚴(yán)重的丟失；Pretest NLM算法不僅更好地抑制了噪聲，而且對(duì)于影像邊緣的保持有了一定的提升；相比之下，本文提出的算法雖然ENL指標(biāo)不如PolSAR NLM方法，但EPD-ROA指標(biāo)反映出本文算法在邊緣和細(xì)節(jié)保持方面相比于前者的優(yōu)勢(shì)性。顯然，本文提出的算法更好地實(shí)現(xiàn)了去噪和細(xì)節(jié)保持之間的平衡。表2為幾種算法對(duì)C波段和L波段影像濾波的評(píng)價(jià)結(jié)果。

圖4 L波段AirSAR影像去噪結(jié)果Fig.4 Filtering results on the L band AirSAR image

表2 真實(shí)影像的濾波定量評(píng)價(jià)結(jié)果Tab.2 Quantitative assessment on the real images

4 結(jié) 論

相干斑濾波是PolSAR數(shù)據(jù)預(yù)處理中重要的一環(huán)。傳統(tǒng)PolSAR非局部均值濾波算法在度量像素或影像塊相似度時(shí)基本利用到的都是含噪聲影像的統(tǒng)計(jì)特性，而沒有考慮到原始像素間的差異性，而且往往無法利用到全極化信息，因此難以獲得穩(wěn)健準(zhǔn)確的相似性測(cè)度。本文提出了一種迭代優(yōu)化的PolSAR非局部均值方法。該方法的創(chuàng)新之處在于：在每次迭代去噪過程中，同時(shí)考慮原始影像全極化噪聲統(tǒng)計(jì)特性和前一次迭代所得的相對(duì)低噪聲度影像的全極化信息，以完善影像塊間極化相似性的度量。模擬的PolSAR數(shù)據(jù)和不同波段的真實(shí)PolSAR影像的去噪試驗(yàn)都表明，本文提出的算法既能有效地抑制相干斑噪聲，又能較好地保持影像的邊緣、強(qiáng)散射點(diǎn)目標(biāo)、極化散射特性等細(xì)節(jié)信息。最后需要指出的是，相比較于傳統(tǒng)的局部濾波器，運(yùn)行效率相對(duì)較低是非局部濾波方法的一大弊端，特別是，本文方法還采取了迭代去噪的策略。因此，如何加速算法的運(yùn)行效率從而提升算法的實(shí)用性，是未來值得研究的方向之一。