韓 萍, 韓賓賓

(中國民航大學天津市智能信號與圖像處理重點實驗室, 天津 300300)

0 引 言

極化合成孔徑雷達(polarimetric synthetic aperture radar,PolSAR)是一種先進的雷達成像系統(tǒng),它具有全天時、全天候、高分辨率的特點。由于系統(tǒng)固有缺陷,PolSAR圖像存在嚴重的相干斑噪聲,它不僅降低了圖像的視覺效果,而且也降低了后續(xù)分類、目標檢測和識別等處理的精度[1-3]。因此,PolSAR圖像相干斑抑制研究是圖像解譯預處理的主要內(nèi)容之一[4-6]。

由于通道串擾,PolSAR圖像的相干斑抑制更加復雜。早期具有代表性的濾波算法有極化白化濾波[7]、最優(yōu)加權(quán)濾波[8]和極化矢量濾波[9],這些算法的主要問題是引入了通道串擾,目標的極化特性沒有得到很好保持;文獻[10] 提出PolSAR圖像的相干斑濾波準則,該準則指出,為避免通道串擾和保持目標的極化特征,應盡可能地選取與濾波像素相近的鄰域像素進行同等程度的濾波?；谠摐蕜t,文獻[10] 同時給出基于邊緣配準窗和最小均方誤差濾波器的精致Lee濾波算法,目前很多濾波算法都在該算法的基礎上發(fā)展而來[11-17],如基于散射模型的PolSAR濾波算法、基于多級分類的濾波算法以及自適應增強Lee濾波算法等。這些算法主要對濾波像素篩選方法進行了改進,可以分為3類:①基于模板窗的鄰域像素篩選,這類算法的優(yōu)點是篩選速度快,缺點是模板篩選隱含假設模板內(nèi)部像素為同質(zhì)像素,這種假設在紋理豐富區(qū)域可能不成立;②基于分類的鄰域像素篩選,這類算法利用未濾波數(shù)據(jù)進行分類,然后以分類結(jié)果為依據(jù)進行同質(zhì)像素篩選,缺點是分類結(jié)果可靠性比較低;③基于相似性度量的鄰域像素篩選,這類算法通過對鄰域像素與濾波像素的相似性進行度量,選擇相似度比較大的鄰域像素作為同質(zhì)像素,能否綜合利用PolSAR圖像的極化信息來度量像素相似性是這類算法的關(guān)鍵。

基于分析,本文提出基于典型散射差異指數(shù)(typical scattering difference index, TSDI)的PolSAR圖像Lee濾波算法。該算法首先引入TSDI度量,TSDI能夠綜合利用PolSAR圖像3個極化通道數(shù)據(jù)度量濾波像素與鄰域像素的相似程度,然后基于該度量采用自適應閾值法篩選出濾波像素的同質(zhì)像素,之后再用篩選的同質(zhì)像素進行Lee濾波。對于篩選閾值,首先利用Parzen窗估計同質(zhì)區(qū)域TSDI的概率分布,然后根據(jù)估計結(jié)果計算閾值。因為算法的同質(zhì)像素篩選過程綜合利用了PolSAR數(shù)據(jù)極化特性,因而同質(zhì)像素與濾波像素更加相近,濾波結(jié)果相對于精致Lee濾波結(jié)果相干斑抑制更加徹底,同時圖像邊緣清晰度和目標極化特性保持更加完好。

1 極化相干矩陣

對于單站雷達,在滿足互易性假設下,極化散射信息可以表示成如下復數(shù)矩陣形式

(1)

式中,SXX和SYY分別為X通道和Y通道發(fā)射接收極化散射系數(shù)；SXY為X通道發(fā)射Y通道接收極化散射系數(shù)。所謂互易性,即X通道發(fā)射Y通道接收回波系數(shù)等于Y通道發(fā)射X通道接收回波系數(shù)。為抑制相干斑,PolSAR數(shù)據(jù)通常會進行多視處理:

(2)

式中,T為常用的PolSAR圖像相干矩陣；Tij為相干矩陣T元素；L表示PolSAR圖像的視數(shù)。

文獻[18-19]提出散射相似性參數(shù)的概念,利用PolSAR目標的相干矩陣與典型散射體的Pauli矢量計算而得

(3)

結(jié)合式(3),容易得到任意相干矩陣與3種典型散射體的相似性參數(shù)

(4)

式中,rs為奇次散射相似性參數(shù);rd為偶次散射相似性參數(shù);rv為體散射相似性參數(shù);y為目標散射功率。由式(4)可以看出,T11可以解釋為目標的奇次散射能量;T22為偶次散射能量;T33為體散射能量。

2 精致Lee濾波

精致Lee濾波主要包括以下兩個步驟:

步驟1利用模板篩選同質(zhì)像素:Lee設計的8個邊緣模板如圖1所示。按文獻[10] 步驟,利用功率圖計算邊緣方向,并將模板內(nèi)部白色像素作為同質(zhì)像素參與到模板中心像素濾波中。

圖1 精致Lee濾波的8種邊緣模板Fig.1 8 edge-directed template of Refine-Lee filter

(5)

3 本文方法

算法首先引入TSDI,該指數(shù)能夠全面地利用PolSAR數(shù)據(jù)極化信息進行像素篩選,因而篩選的鄰域像素比精致Lee濾波更加相近中心像素,濾波后圖像邊緣和目標極化信息保持更加完好。本文所提算法的基本思路為:①利用濾波像素與鄰域像素之間的TSDI,采用自適應閾值法選出濾波像素的同質(zhì)鄰域像素;②利用同質(zhì)鄰域像素進行Lee濾波。算法流程圖如圖2所示。

圖2 算法流程圖Fig.2 Flow chart of algorithm

3.1 典型散射差異指數(shù)TSDI

在局部區(qū)域內(nèi),如果兩個像素的奇次散射能量、偶次散射能量和體散射能量都相近,則可以認為這兩個像素屬于相同的散射體?；谠撛瓌t,引入TSDI概念:

(6)

(1)g≥1,當且僅當兩個目標的3種散射能量都相等時,典型散射差異指數(shù)g=1;

(2) 兩目標相干矩陣的差異性越大,典型散射差異指數(shù)g值越大;

(3) 典型散射差異指數(shù)g只與兩目標相干矩陣對角元素的相對差異有關(guān)。

3.2 基于Parzen窗的TSDI分布估計

為確定同質(zhì)像素篩選的判決閾值,對同質(zhì)區(qū)域的TSDI概率分布進行統(tǒng)計,并用Parzen窗概率密度估計方法[20-21]擬合概率分布。Parzen窗概率密度估計方法在很多領域都有廣泛應用,它以數(shù)據(jù)樣本為中心利用窗函數(shù)對概率密度函數(shù)進行插值擬合。當數(shù)據(jù)樣本足夠大時,該方法能夠給出相當準確的概率密度估計。

假設樣本集合X={x1,x2,…,xM},Parzen窗概率密度函數(shù)為

(7)

式中,ke(·)表示Parzen窗的核函數(shù),可以是矩形函數(shù)、高斯函數(shù)等；h表示核函數(shù)的寬度；M表示樣本個數(shù)。本文使用最常用的Gauss窗函數(shù)和文獻[22] 給出的最優(yōu)窗寬度,即

(8)

式中,σ表示樣本標準差。

考慮到TSDI取值范圍比較大,為減少估計誤差,先對TSDI取自然對數(shù)再用Parzen窗進行估計

x=ln(g)

(9)

得到同質(zhì)區(qū)域TSDI的概率密度p(x)后,由給定檢測率Pd可確定對應的篩選閾值th：

(10)

利用篩選閾值,在中心像素鄰域內(nèi)進行如下方式的同質(zhì)像素篩選：

(11)

式中,δk,l=1表示窗內(nèi)(k,l)位置像素為當前濾波像素的同質(zhì)像素,參與Lee濾波;相反,δk,l=0表示(k,l)位置像素為非同質(zhì)像素,不參與濾波。

4 實驗及分析

為了驗證本文算法的有效性,分別利用美國AIRSAR系統(tǒng)和UAVSAR系統(tǒng)采集的多幅圖像進行實驗,由于篇幅所限本文只給出San Francisco地區(qū)和Gulf Coast地區(qū)的部分實驗結(jié)果。AIRSAR系統(tǒng)采集的San Francisco地區(qū)PolSAR數(shù)據(jù)如圖3(a)所示,數(shù)據(jù)大小為900×700,數(shù)據(jù)包含了海洋、森林、城區(qū)等典型地物。UAVSAR系統(tǒng)采集的Gulf Coast地區(qū)PolSAR數(shù)據(jù)如圖3(b)所示,數(shù)據(jù)大小為401×675,數(shù)據(jù)場景為局部機場區(qū)域,既有勻質(zhì)的道路、水體和植被,又有紋理豐富的城區(qū),還有一些明顯的點、線目標。

圖3 實驗數(shù)據(jù)PauliRGB圖Fig.3 PauliRGB image of experiment data

圖3(a)中圓角矩形區(qū)域是海洋,以表面散射為主,總像素為22 000,圖3(b)圓角矩形區(qū)域是植被,以體散射為主,總像素數(shù)為21 343。這些區(qū)域的散射特性分布比較均勻,適合用來估計TSDI的概率分布。由于數(shù)據(jù)樣本比較大,可以把平均相干矩陣作為目標真實相干矩陣。區(qū)域內(nèi)相干矩陣與平均相干矩陣的自然對數(shù)TSDI概率分布和Parzen窗估計結(jié)果如圖4所示?；鸩駰U圖為真實TSDI概率分布,實曲線為擬合結(jié)果,由圖4可知,Parzen窗較好地擬合了TSDI概率分布。χ2檢驗是一種通用的假設檢驗方法,如果χ2匹配系數(shù)小于卡方臨界值,則接受估計結(jié)果作為真實結(jié)果。San Francisco地區(qū)和Gulf Coast地區(qū)的χ2匹配系數(shù)分別為210.9和229.9,χ2臨界值為254.13,可以接受兩種估計分布為真實概率分布。比較圖4(a)和圖4(b)可知,雖然San Francisco和Gulf Coast的采集系統(tǒng)不同,目標散射類型也不同,但TSDI概率分布仍一致,說明TSDI具有較高的魯棒性。[0.5,1.0]范圍內(nèi)的TSDI和對應檢測概率如表1所示。由表1可知，當ln(th)>0.5即th>1.65時兩個地區(qū)都能達到90%以上的檢測概率;當ln(th)>1.0即th>2.72時能達到98%以上,在[1.65,2.72]之間隨機選取閾值也能取得不錯的篩選結(jié)果。針對不同PolSAR圖像可以采用兩種閾值選取方式:①當PolSAR圖像存在大片勻質(zhì)區(qū)域時可以用Parzen窗估計TSDI概率分布,再由檢測概率確定篩選閾值;②當PolSAR圖像不存在大片勻質(zhì)區(qū)域時,可在[1.65,2.72]隨機選取閾值。

圖4 均勻區(qū)域典型散射差異指數(shù)分布與擬合Fig.4 Distribution and fitting result of typical scattering difference index in homogeneous region

對數(shù)典型散射差異指數(shù)ln(g)0.50.60.70.80.91.0典型散射差異指數(shù)g1.651.822.012.232.462.72SanFrancisco檢驗概率Pd0.920.940.970.980.990.99GulfCoast檢驗概率Pd0.950.960.970.980.980.99

4.1 視覺效果評價

為驗證改進算法的有效性,首先從視覺角度對濾波結(jié)果進行比較。PauliRGB圖除了能夠表示目標的強度信息,也能表示目標的散射信息和極化信息,在PauliRGB圖中進行視覺分析能更全面地表現(xiàn)各種算法的濾波效果。

San Francisco地區(qū)的濾波結(jié)果如圖5所示,其中圖5(a)為原始PauliRGB圖,本文對圖中一些典型的散射目標進行了標注,1區(qū)為山林,2區(qū)為海洋,3區(qū)為城區(qū),4區(qū)為金門大橋,5區(qū)為高爾夫球場,6區(qū)為跑馬場,7區(qū)為城區(qū)主干道。實驗分別選取了5×5,7×7,9×9,11×11這4種窗,對精致Lee濾波和本文算法進行比較。由于文獻[10] 只給出7×7窗的濾波結(jié)果,為便于比較本文也只給出窗口為7×7的濾波結(jié)果。圖5(b)為精致Lee濾波后的PauliRGB圖;圖5(c)為本文方法濾波后的PauliRGB圖。從圖中可以看出,兩種濾波算法都較好地抑制了相干斑,海洋和跑馬場區(qū)域的視覺效果有所改善;但是,精致Lee濾波算法在紋理信息豐富的區(qū)域產(chǎn)生了過度平滑,比如高爾夫球場的紋理模糊嚴重,城區(qū)的塊狀效應非常明顯,且城區(qū)的道路出現(xiàn)了斷斷續(xù)續(xù)的情況,甚至一些小的道路消失不見。而本文算法則在這方面表現(xiàn)較好,圖中城區(qū)、高爾夫球場的紋理和城區(qū)道路都清晰可見。Gulf Coast地區(qū)的濾波結(jié)果如圖6所示。從整體上看,本文算法的視覺效果更加自然,同質(zhì)區(qū)域非常平整,點、線目標清晰,城區(qū)的清晰度明顯改善。

4.2 定量評價

采用常用的PolSAR圖像濾波評價指標等效視數(shù)(equivalent number of looks,ENL)和邊緣保持指數(shù)(edge preserved index,EPI)進行定量分析。ENL越大表明圖像相干斑越小,算法抑制效果越好;EPI越接近1,表明算法目標保持效果越好。圖7給出了計算各項指標的采樣區(qū)域;表2和表3分別給出了San Francisco和Gulf Coast地區(qū)的計算結(jié)果。

圖5 San Francisco地區(qū)濾波效果Fig.5 Filter effect in San Francisco region

圖6 Gulf Coast地區(qū)濾波效果Fig.6 Filter effect in Gulf Coast region

圖7 ENL和EPI分析的數(shù)據(jù)采樣區(qū)域Fig.7 Sampling regions for ENL and EPI analysis

圖7中選擇A區(qū)、B區(qū)計算ENL,其中San Francisco地區(qū)的A區(qū)和B區(qū)分別為海洋和森林,Gulf Coast地區(qū)的A區(qū)和B區(qū)分別為灌木林和草坪,這些區(qū)域的散射特性比較均勻單一,可以看成同質(zhì)區(qū)域;選擇C區(qū)、D區(qū)計算EPI,San Francisco地區(qū)的C區(qū)和D區(qū)為山頂建筑物和城區(qū)邊界,Gulf Coast地區(qū)的C區(qū)和D區(qū)為機場周邊建筑,這些區(qū)域邊緣豐富。如表2和表3所示,濾波前A區(qū)和B區(qū)的ENL都比較低,相干斑嚴重;濾波后兩個區(qū)域的ENL都有所增加,而且兩種方法都符合一個比較普遍的認知:窗口越大,相干斑抑制越好,ENL越大。對比來看,本文方法的ENL要高于精致Lee濾波算法,這說明新方法的相干斑抑制能力更強。相對應的,濾波后兩種濾波算法都損失了一定的邊緣信息,從表1看出,San Francisco地區(qū)的EPI,除9×9窗外新方法都高于精致Lee濾波;Gulf Coast地區(qū)的EPI,新方法在各種窗口下都高于精致Lee濾波,這說明新方法的圖像邊緣保持能力更好。而且,隨著窗口增大本文方法的ENL和EPI都有增大的趨勢,這是因為窗口越大同質(zhì)區(qū)篩選可選的鄰域像素越多。但是,隨著窗口增大本文算法的運算量也迅速增大,如以w×w窗對大小m×n的PolSAR圖像進行濾波,同質(zhì)區(qū)篩選需要進行6mnw2次除法運算和5mnw2次加法運算,運算量與窗口尺寸的平方成正比。因此，本文算法適合應用在對圖像實時處理要求不高的場合。綜合考慮濾波效果和運算量,9×9和11×11是比較理想的濾波窗口。

表2 San Francisco地區(qū)的ENL和EPI定量分析表

表3 Gulf Coast地區(qū)的ENL和EPI定量分析表

4.3 極化特性保持評價

為比較兩種算法的極化特性保持能力,圖8和圖9給出了7×7窗口的點目標同極化響應和交叉極化響應圖。采樣點的位置如圖7箭頭所示,San Francisco地區(qū)的采樣點為海洋中的艦船,坐標為(260,332);Gulf Coast地區(qū)的采樣點為停機坪上的飛機,坐標為(170,282),它們都與周圍目標存在較大的散射差異。從圖8和圖9中容易看出本文算法的極化保持能力要好于精致Lee濾波。進一步的分析表明,這是因為精致Lee濾波篩選出的鄰域像素與中心像素存在較大差異所致,而本文方法篩選的鄰域像素與中心像素差異較小。圖10給出了San Francisco地區(qū)采樣點的濾波像素篩選情況,白色表示參與濾波的像素。精致Lee算法選擇中心點下方所有像素參與濾波,這些像素中混雜了很多海洋目標,而本文算法只選用與中心像素比較相近的10個像素參與濾波。

圖8 San Francisco地區(qū)圖中(260,332)的目標同極化和交叉極化響應圖Fig.8 Co-polarization and cross-polarization signatures of(260,332) pixel in San Francisco region

圖9 Gulf Coast地區(qū)圖中(170,282)的目標同極化和交叉極化響應圖Fig.9 Co-polarization and cross-polarization signatures of (170,282)pixel in Gulf Coast region

圖10 San Francisco地區(qū)極化特性分析采樣點篩選的鄰域濾波像素Fig.10 Selected neighborhood pixels of sampling point for polarimetric characteristic analysis in San Francisco

5 結(jié) 論

針對濾波后PolSAR圖像的邊緣和目標的極化特性保持不夠理想的問題,本文提出基于TSDI的PolSAR圖像Lee濾波算法。算法利用濾波像素和鄰域像素之間的TSDI,采用自適應閾值法選出濾波像素的同質(zhì)像素進行Lee濾波。利用實測數(shù)據(jù)進行實驗,結(jié)果表明本文濾波效果優(yōu)于精致Lee濾波。算法在濾波過程中需要逐個對圖像像素進行同質(zhì)像素篩選,其運算量與鄰域窗口平方成正比,運算量相對于精致Lee濾波較大,如何減少運算量是下一步要研究的重點。

[1] GU J, YANG J, ZHANG H, et al. Speckle filtering in polarimetric SAR data based on the subspace decomposition[J]. IEEE Trans.on Geoscience & Remote Sensing, 2004, 42(8): 1635-1641.

[2] KERSTEN P R, LEE J S, AINSWORTH T L. Unsupervised classification of polarimetric synthetic aperture radar images using fuzzy clustering and EM clustering[J]. IEEE Trans.on Geoscience & Remote Sensing, 2005, 43(3): 519-527.

[3] VASILE G, TROUVE E, LEE J S, et al. Intensity-driven adaptive neighborhood technique for polarimetric and interferometric SAR parameters estimation[J]. IEEE Trans.on Geoscience & Remote Sensing, 2006, 44(6): 1609-1621.

[4] 周曉光, 匡綱要, 萬建偉. 多極化SAR圖像斑點抑制綜述[J]. 中國圖象圖形學報, 2008, 13(3): 377-385.

ZHOU X G, KUANG G Y, WAN J W. A review of polarimetric SAR speckle reduction[J]. Journal of Image & Graphics, 2008, 13(3): 377-385.

[6] LEE J S, AINSWORTH T L, WANG Y. On polarimetric SAR speckle filtering[C]∥Proc.of the IEEE Geoscience and Remote Sensing Symposium, 2012: 111-114.

[7] NOVAK L M, BURL M C. Optimal speckle reduction in polarimetric SAR imagery[J]. IEEE Trans.on Aerospace & Electronic Systems, 1990, 26(2): 293-305.

[8] LEE J, HOPPEL K, MANGO S A. Unsupervised estimation of speckle noise in radar images[J]. International Journal of Imaging Systems & Technology, 2010, 4(4): 298-305.

[9] LIN Q, ALLEBACH J P. Combating speckle in SAR images: vector filtering and sequential classification based on a multiplicative noise model[J]. IEEE Trans.on Geoscience & Remote Sensing, 1990, 28(28): 647- 653.

[10] LEE J S, GRUNES M R, GRANDI G D. Polarimetric SAR speckle filtering and its impact on classification[C]∥Proc.of the Geoscience and Remote Sensing, 1997: 1038-1040.

[11] LEE J S, GRUNES M R, SCHULER D L, et al. Scattering-model-based speckle filtering of polarimetric SAR data[J]. IEEE Trans.on Geoscience & Remote Sensing, 2006, 44(1): 176-187.

[12] 韓萍,董菲,吳仁彪.基于多級分類的極化SAR圖像斑點抑制[J]. 系統(tǒng)工程與電子技術(shù),2012,34(7): 1360-1365.

HAN P, DONG F, WU R B. PolSAR image speckle reduction based on multi-stage classification[J]. Systems Engineering and Electronics, 2012, 34(7): 1360-1365.

[13] LANG F, YANG J, LI D. Adaptive-window polarimetric SAR image speckle filtering based on a homogeneity measurement[J]. IEEE Trans.on Geoscience & Remote Sensing, 2015, 53(10): 5435-5446.

[14] LEE J S, AINSWORTH T L, WANG Y, et al. Polarimetric SAR speckle filtering and the extended sigma filter[J].IEEE Trans.on Geoscience & Remote Sensing,2015,53(3):1150-1160.

[15] 李平湘, 馬曉雙, 沈煥鋒. 基于局部同質(zhì)區(qū)最大化搜索的極化SAR影像最小均方誤差濾波方法[J]. 測繪地理信息, 2016, 41(3): 1-7.

LI P X, MA X S, SHEN H F. A minimum-mean-square error filter based on the search of the maximal local homogenous area for PolSAR image despeckling[J]. Journal of Geomatics, 2016, 41(3): 1-7.

[16] 楊學志,陳靖,周芳,等.基于同質(zhì)像素預選擇的極化SAR圖像非局部均值濾波[J].電子與信息學報,2015,37(12):2991-2999.

YANG X Z, CHEN J, ZHOU F, et al. Polarimetric SAR image despeckling using non local means filter based on homogeneous pixels preselection[J]. Journal of Electronics & Information Technology, 2015, 37(12): 2991-2999.

[17] 郎豐鎧, 楊杰, 李德仁. 極化SAR圖像自適應增強Lee濾波算法[J]. 測繪學報, 2014, 43(7): 690-697.

LANG F K, YANG J, LI D R. An adaptive enhanced Lee speckle filter for polarimetric SAR image[J]. Acta Geodaetica et Cartographica Sinica, 2014, 43(7): 690-697.

[18] YANG J, PENG Y N, LIN S M. Similarity between two scattering matrices[J].Electronics Letters, 2001,37(3):193-194.

[19] AN W, ZHANG W, YANG J, et al. Similarity between two targets and its application to polarimetric target detection for sea area[C]∥Proc.of the 24th Progress in Electromagnetics Research Symposium, 2008:515-520.

[20] 張宏稷,楊健,李延,等.基于條件熵和Parzen窗的極化SAR艦船檢測[J].清華大學學報,2012(12):1693-1697.

ZHANG H J, YANG J, LI Y, et al. Ship detection in polarimetric SAR images based on the conditional entropy and Parzen windows[J]. Journal of Tsinghua University, 2012(12): 1693-1697.

[21] DUDA R O, HART P E, STORK D G. Pattern classification[M].2nd ed. New York:Wiley, 2001.

[22] WASSERMAN L. All of nonparametric statistics (springer texts in statistics)[M]. New York: Springer-Verlag, 2006.