陳 軍,杜培軍,2,譚 琨
(1.中國(guó)礦業(yè)大學(xué)國(guó)土環(huán)境與災(zāi)害監(jiān)測(cè)國(guó)家測(cè)繪地理信息局重點(diǎn)實(shí)驗(yàn)室,徐州 221116;2.南京大學(xué)江蘇省地理信息技術(shù)重點(diǎn)實(shí)驗(yàn)室,南京 210023)
全極化合成孔徑雷達(dá)(polarimetric synthetic aperture Radar,PolSAR)系統(tǒng)是在傳統(tǒng)的單極化、雙極化和多極化雷達(dá)系統(tǒng)基礎(chǔ)上發(fā)展起來(lái)的新型主動(dòng)微波成像系統(tǒng),如何通過(guò)全極化SAR系統(tǒng)提供的極化信息對(duì)地物進(jìn)行分類已成為全極化SAR圖像應(yīng)用的一個(gè)重要方面[1]。傳統(tǒng)的單極化和多極化SAR圖像處理通過(guò)提取紋理信息、極化信息等進(jìn)行分類,然而,無(wú)論單極化還是多極化SAR系統(tǒng)獲取的都是部分極化信息;而全極化SAR系統(tǒng)包含同極化、交叉極化等不同的極化信息,可以全面反映地物目標(biāo)的物理性質(zhì)。
國(guó)內(nèi)外學(xué)者在全極化SAR圖像極化特征提取方面開(kāi)展了大量研究[2-7]。利用極化分解獲得的特征可以對(duì)全極化 SAR 圖像進(jìn)行分類[8-15]。Lee等[8]首先用H/Alpha分解法提取的特征對(duì)圖像進(jìn)行初始聚類,得到8個(gè)聚類中心;然后用描述多視協(xié)方差矩陣的復(fù)Wishart迭代分類器對(duì)圖像進(jìn)行分類(簡(jiǎn)稱 H/Alpha- Wishart分類器)。Pottier等[9]提出了H/Alpha/A-Wishart分類器,在H/Alpha分解的基礎(chǔ)上加入A特征,將圖像聚成16類,再對(duì)圖像進(jìn)行復(fù)Wishart迭代分類。Lee等[10]進(jìn)一步結(jié)合Freeman分解法對(duì)上述方法進(jìn)行改進(jìn),提出Freeman-H/Alpha-Wishart分類器,該分類器通過(guò)Freeman分解提取體散射、表面散射和二次散射信息,然后根據(jù)H/Alpha/A-Wishart分類器進(jìn)行分類,分類中要求各類別必須保持純凈性(即各類別內(nèi)部各像元必須基于相同的散射性質(zhì)),必須是體散射、表面散射或二次散射中的一種,該分類器已經(jīng)隱含參考Freeman分解進(jìn)行分類,但并沒(méi)有直接利用Freeman分解得到的特征,而是簡(jiǎn)單地要求分類結(jié)果必須具有相同的散射性質(zhì),因此還具有一定局限性。
無(wú)論H/Alpha-Wishart分類器、Freeman-H/Alpha-Wishart分類器或者 H/Alpha-SPAN[11]分類器,其分類結(jié)果均明顯依賴于H/Alpha初始聚類結(jié)果,導(dǎo)致最終分類精度不高。針對(duì)這一問(wèn)題,本文提出一種改進(jìn)的分類方法,該方法充分利用了全極化SAR圖像的表面散射、體散射、二面角散射、螺旋體散射和線散射等極化特征,在沒(méi)有人工干預(yù)的情況下能夠得到較好的聚類結(jié)果。
1.1.1 H/Alpha/A 分解與初始聚類
H/Alpha/A分解主要是提取全極化SAR圖像的散射熵H、散射角α(即Alpha,下同)和反熵A信息。其基本原理如下:
極化相干矩陣T可以表示為3個(gè)獨(dú)立的相干矩陣之和[2],即
式中:λi為特征值(i=1,2,3);ei為矩陣T的第i個(gè)復(fù)特征向量;為ei的共軛矩陣。
在此基礎(chǔ)上,Cloude等[2]定義了熵H,用于表示地物散射的隨機(jī)程度,即
式中Pi為各特征值所占的百分比,即
散射角α代表二面角散射和表面散射的平均散射程度,即
反熵A表示優(yōu)勢(shì)散射以外的較弱散射分量間的關(guān)系[9],即
根據(jù)H和α的不同取值,即可確定圖像每個(gè)像元在H/Alpha平面(圖1)中的位置,從而將像元?dú)w入8個(gè)不同的類別;對(duì)每個(gè)類別取平均操作,得到初始聚類結(jié)果。加入A后,設(shè)置1個(gè)閾值,可以得到16個(gè)不同的類別。
圖1 H/Alpha 分類圖[2]Fig.1 Classification image of H/A lpha[2]
1.1.2 基于復(fù)Wishart分布的極大似然分類
以1.1.1節(jié)初始聚類結(jié)果作為中心,進(jìn)行基于復(fù)Wishart分布的極大似然分類。其原理[8]如下:
極化協(xié)方差矩陣的概率密度分布函數(shù)服從n個(gè)自由度的復(fù)Wishart分布,其分布函數(shù)PT為
其中
式中:n為視數(shù);q代表不同的極化方式,在互易性情況下,q=3;Tr為矩陣的跡;C為極化協(xié)方差矩陣;∑=E(C);K為歸一化因子;Γ(·)為Gamma函數(shù)。對(duì)矩陣進(jìn)行線性變換可以得到極化相干矩陣T。因此,極化相干矩陣的概率密度分布函數(shù)也服從n個(gè)自由度的復(fù)Wishart分布。
進(jìn)一步對(duì)圖像進(jìn)行開(kāi)窗操作,計(jì)算窗口內(nèi)所有像元的平均極化相干矩陣;然后分別求該矩陣與8個(gè)類別中心的距離[8],即
式中:∑i為第i個(gè)類別中心的極化相干矩陣;T為窗口內(nèi)所有像元極化相干矩陣的均值矩陣。比較8個(gè)距離,將像元?dú)w為距離最小的那一類。
H/Alpha的初始聚類過(guò)程是根據(jù)圖1內(nèi)H和α的不同取值進(jìn)行初始類別劃分的。而實(shí)際上,H和α取值和不同的地物類型并沒(méi)有嚴(yán)格的對(duì)應(yīng)關(guān)系[8]。另外,該聚類在算法實(shí)現(xiàn)上也是采用近似的形式,沒(méi)有完全按照?qǐng)D1內(nèi)的曲線進(jìn)行邊界劃分,因此,由圖1得到的初始聚類結(jié)果存在一定的偏差。對(duì)于基于初始聚類中心的分類算法,初始聚類中心的位置直接影響最終分類結(jié)果,若該位置選擇不當(dāng),必然會(huì)降低分類精度。因此,對(duì)極化特征分解和初始聚類過(guò)程進(jìn)行改進(jìn),有望提高全極化SAR圖像的分類精度。
本研究從2個(gè)方面優(yōu)化初始聚類效果:①引入MCSM分解對(duì)全極化SAR圖像進(jìn)行極化分解,提取表面散射、二次散射、體散射、螺旋體散射和線散射等更為全面的基本散射信息作為特征,解決僅使用H,α和A特征進(jìn)行初始聚類的不足;②使用聚類效果更好的ISODATA聚類算法代替圖1中的簡(jiǎn)單聚類算法,以提高初始聚類精度。
MCSM分解的極化協(xié)方差陣[6]為
式中:fs,fd,fv,fh和fw分別為表面(surface)散射、二次(double)散射、體(volume)散射、螺旋體(helix)散射和線(wire)散射系數(shù);α,β,γ和ρ為未知數(shù);α*,β*,γ*和 ρ*分別為 α,β,γ 和 ρ的共軛。求解式(9),即可得到各散射系數(shù)及α和β的值;γ和ρ可以通過(guò)下式求解[6],即
式中SHH,SVV,SHV和SVV分別表示極化方式為HH,VV,HV和VV的散射矩陣。
根據(jù)各散射成分的系數(shù),即可計(jì)算各散射功率的值,計(jì)算公式[6]為
式中:Ps為表面散射功率;Pd為二次散射功率;Pv為體散射功率;Ph為螺旋體散射功率;Pw為線散射功率。
對(duì)極化協(xié)方差矩陣進(jìn)行線性變換可以得到極化相干矩陣。因此,極化相干矩陣同樣可以通過(guò)MCSM分解得到以上各散射成分的系數(shù)。從當(dāng)前極化目標(biāo)分解算法的發(fā)展和評(píng)述來(lái)看,MCSM分解的方法能充分提取出各種自然和人工地物的各種散射成分的信息,將這些信息作為特征向量可以較全面地反映圖像的物理性質(zhì)[16]。
ISODATA是一種運(yùn)行速度較快、可以自動(dòng)合并和分裂聚類結(jié)果的非監(jiān)督分類或聚類方法。與K-means聚類算法相比,ISODATA可以進(jìn)行類別的合并和分裂,靈活性較好,聚類結(jié)果更為合理。大量實(shí)驗(yàn)分析表明,ISODATA聚類具有算法實(shí)現(xiàn)簡(jiǎn)單、運(yùn)算速度快、聚類效果較好等特點(diǎn)。
傳統(tǒng)的ISODATA聚類算法多采用歐式距離計(jì)算樣本與聚類中心的距離。但考慮到極化數(shù)據(jù)不同于光學(xué)數(shù)據(jù)的概率分布,本文對(duì)比歐式距離和馬氏距離后,采用馬氏距離進(jìn)行距離計(jì)算,即
式中樣本si與中心mj分別為
改進(jìn)MCSM-Wishart分類器的核心思想是:①引入MCSM分解對(duì)圖像進(jìn)行極化分解,提取表面散射、二次散射、體散射、螺旋體散射和線散射等更為全面的基本散射信息作為特征;②使用聚類結(jié)果更好的ISODATA算法。因此,MCSM-Wishart分類器的改進(jìn)主要包括4個(gè)部分:圖像預(yù)處理、極化特征的提取、初始聚類結(jié)果的生成和復(fù)Wishart分布的迭代分類。具體分類步驟如下:
1)數(shù)據(jù)預(yù)處理。讀入全極化SAR圖像,進(jìn)行多視處理,按照通用方法將多視比設(shè)為7∶1。設(shè)置濾波窗口大小為7像元×7像元,進(jìn)行Refined Lee濾波,得到圖像的極化相干矩陣T。
2)MCSM分解。設(shè)置的窗口大小與步驟1)保持一致,對(duì)極化相干矩陣T進(jìn)行MCSM分解,提取表面散射、二次散射、體散射、螺旋體散射和線散射等多種散射成分,將其組合成特征向量,用以表征極化信息。
3)ISODATA初始聚類。設(shè)置ISODATA聚類算法的初始條件(最小、最大類別等),對(duì)步驟2)產(chǎn)生的特征向量進(jìn)行迭代聚類,生成初始聚類結(jié)果。
4)計(jì)算聚類中心。在各聚類內(nèi)部進(jìn)行平均操作,計(jì)算出初始聚類中心。
5)計(jì)算距離并聚類。選擇步驟4)生成的聚類中心,設(shè)置循環(huán)條件,進(jìn)行基于復(fù)Wishart分布的迭代分類;根據(jù)式(8)計(jì)算圖像內(nèi)所有像元到各聚類中心的距離,將像元?dú)w入距離最小的類別。
6)判斷是否滿足循環(huán)終止條件,若不滿足,返回步驟4)繼續(xù)進(jìn)行計(jì)算;否則,完成分類。
為驗(yàn)證本文算法,選擇南京溧水和鹽城濱海研究區(qū)2景ALOS的1.1級(jí)PALSAR(phased array L-band SAR)圖像(全極化方式)的子圖像為數(shù)據(jù)源進(jìn)行實(shí)驗(yàn)。2景圖像的獲取時(shí)間分別為2011年3月22日和2009年4月9日。其中,研究區(qū)1(南京溧水)圖像大小為400像元×400像元,研究區(qū)2(鹽城濱海)圖像大小為500像元×500像元。研究區(qū)內(nèi)包括水體、建筑物、耕地、林地和道路等地物類型。實(shí)驗(yàn)結(jié)果采用目視解譯數(shù)據(jù)和Google Earth圖像,通過(guò)綜合對(duì)比的方法確定測(cè)試樣例,進(jìn)行精度評(píng)價(jià)。研究區(qū)1的Google Earth圖像成像時(shí)間為2013年4月10日;研究區(qū)2的TM圖像成像時(shí)間為2010年4月29日。2個(gè)研究區(qū)均屬縣級(jí)行政單位下的鄉(xiāng)鎮(zhèn)地區(qū),全極化SAR圖像、TM圖像和Google Earth圖像在獲取時(shí)間間隔內(nèi)土地利用變化較小,因此雖然在時(shí)相上有一定差距,但仍具有參考價(jià)值。
實(shí)驗(yàn)采用的處理軟件為POLSARPro(the polarimetric SAR data processing and educational tool)軟件,數(shù)據(jù)預(yù)處理包括多視處理、濾波和極化特征分解等,傳統(tǒng) H/Alpha/A-Wishart分類算法及監(jiān)督-Wishart分類算法采用該軟件自帶的功能,改進(jìn)算法采用C語(yǔ)言實(shí)現(xiàn)。
研究區(qū)1和研究區(qū)2的Pauli基分解[1](極化校準(zhǔn))及假彩色合成圖像如圖2所示。
圖2 研究區(qū)Pauli基分解及假彩色合成圖像Fig.2 False color com posite images and Pau li decomposition of two study areas
首先進(jìn)行圖像預(yù)處理工作,包括多視、濾波和極化特征分解。對(duì)圖像進(jìn)行多視處理,多視比設(shè)為7∶1;對(duì)多視處理后的圖像進(jìn)行窗口大小為7像元×7像元的Refined Lee濾波處理;通過(guò)MCSM分解提取出各種散射成分。然后進(jìn)行初始聚類,ISODATA算法設(shè)置的初始條件為:最小類別為8,最大類別為12,類內(nèi)最小個(gè)數(shù)為2,類間最小距離為4,最大迭代次數(shù)為20,最大類標(biāo)準(zhǔn)差為1。通過(guò)聚類結(jié)果計(jì)算聚類中心,最后進(jìn)行復(fù)Wishart分類,得到分類結(jié)果。
作為對(duì)比,本文還分別進(jìn)行了H/Alpha-Wishart和H/Alpha/A-Wishart算法的非監(jiān)督分類實(shí)驗(yàn),并選擇訓(xùn)練樣本進(jìn)行了監(jiān)督-Wishart分類實(shí)驗(yàn)。其他分類器如 Freeman-H/Alpha-Wishart,H/Alpha-SPAN等的分類精度與H/Alpha/A-Wishart分類器的分類精度相似。為便于研究,僅以H/Alpha-Wishart和H/Alpha/A-Wishart算法進(jìn)行實(shí)驗(yàn)。為便于比較分類結(jié)果,結(jié)合訓(xùn)練樣本和監(jiān)督分類結(jié)果對(duì)非監(jiān)督聚類后的結(jié)果進(jìn)行了類別合并(圖3和圖4)。
圖3 研究區(qū)1不同算法分類結(jié)果Fig.3 Classification results by different algorithms for study area one
從圖3(a)可以看出,在研究區(qū)1,原H/Alpha-Wishart算法將水體和耕地混分。原因在于,當(dāng)水體面積較少時(shí),H/Alpha分解提取出的特征不足以反映地物的極化性質(zhì),因此影響了最終的精度。雖然對(duì)H/Alpha/A-Wishart算法有所改善,但整體精度仍低于改進(jìn)后的MCSM-Wishart算法精度。通過(guò)對(duì)比圖3(b)與圖3(c)可以看出,在研究區(qū)1,H/Alpha/A-Wishart算法存在耕地和林地混分、林地和建筑物混分的情況,而MCSM-Wishart算法則較少存在以上混分現(xiàn)象。另外,MCSM-Wishart算法提取出的水體也多于H/Alpha/A-Wishart算法,這是因?yàn)?,MCSM算法提取出了更多的極化特征用于分類。
圖4 研究區(qū)2不同方法分類結(jié)果Fig.4 C lassification resu lts by different algorithm s for study area two
對(duì)比圖4(a)(b)(c)可以看出,在研究區(qū)2,H/Alpha-Wishart和H/Alpha/A-Wishart算法均不同程度地存在建筑和道路混分、道路和水體混分的情況;而MCSM-Wishart算法則能較好地區(qū)分出建筑、道路和水體,圖4下部的水體中混有道路的情況得到了改善,圖4右上部的建筑物中混有道路的情況得到改進(jìn)。主要原因是MCSM方法提取出了螺旋體散射、線散射等較多的能反映城市地區(qū)地物類型的特征,因此能較好地對(duì)城市及周邊的不同人造地物進(jìn)行類別劃分。
為比較各算法的分類精度,結(jié)合Google Earth圖像、TM圖像對(duì)2個(gè)研究區(qū)進(jìn)行人工采樣,手工選擇了各種地物類型的測(cè)試樣本如表1所示。
表1 測(cè)試樣本Tab.1 Test sam p les (像元)
4種算法的分類精度見(jiàn)表2。
表2 4種算法的分類精度Tab.2 Classification results of four kinds of algrithm
從表2可以看出,無(wú)論研究區(qū)1還是研究區(qū)2,MCSM-Wishart算法的分類精度均有較大提升,分別比 H/Alpha/A-Wishart算法提高了4.22%和11.81%;Kappa 系數(shù)則分別提高了0.05 和0.16。
與監(jiān)督-Wishart算法相比,MCSM-Wishart算法精度稍差,說(shuō)明極化特征提取方法還有需要完善之處。4種算法的執(zhí)行效率見(jiàn)表3。
表3 4種算法的執(zhí)行效率Tab.3 Efficiency of four kinds of algorithm(s)
從表3中4種算法的運(yùn)行時(shí)間來(lái)看,MCSMWishart算法效率高于H/Alpha/A-Wishart算法。雖然監(jiān)督-Wishart算法需要的時(shí)間更短,但這個(gè)時(shí)間并未包含選擇訓(xùn)練樣本的時(shí)間、也沒(méi)有考慮訓(xùn)練樣本選擇的實(shí)際困難。因此總的來(lái)看,MCSMWishart算法是一種有效的非監(jiān)督分類算法。
綜上所述,無(wú)論總體精度、Kappa系數(shù)還是運(yùn)行效率,MCSM-Wishart算法均較 H/Alpha-Wishart算法和H/Alpha/A-Wishart算法有了明顯的提高。
1)經(jīng)典的全極化SAR圖像非監(jiān)督分類主要以H/Alpha/A分解提取的極化特征為基礎(chǔ)進(jìn)行圖像分類,但總體分類精度往往較低,不能滿足應(yīng)用要求。造成這種情況的原因是多方面的,如該分解法并不能完全提取反映全極化SAR圖像的全部特征,選擇的分類器并不能很好地對(duì)圖像進(jìn)行分類等。
2)引入MCSM算法是提高全極化SAR圖像分類準(zhǔn)確性的有效途徑之一。本文在傳統(tǒng)的極化信息提取基礎(chǔ)上,引入多極化特征提取的MCSM算法,全面提取出反映城市地物類型的極化特征;然后通過(guò)改進(jìn)的聚類算法進(jìn)行初始聚類,以此為基礎(chǔ)對(duì)原H/Alpha/A-Wishart分類器進(jìn)行改進(jìn)。
3)以ALOS的PALSAR圖像為遙感數(shù)據(jù)源進(jìn)行綜合分類實(shí)驗(yàn)的結(jié)果表明,改進(jìn)的MCSM-Wishart分類算法較經(jīng)典的H/Alpha/A-Wishart分類器在分類精度、Kappa系數(shù)、算法效率等指標(biāo)上均有很大的提高。因此,提取全極化SAR圖像的線散射等極化特征,然后與二次散射、表面散射和螺旋體散射等經(jīng)典的極化特征相結(jié)合,輔以速度快、精度高的ISODATA聚類算法,可以有效地提高全極化SAR圖像分類的精度。
4)下一步的研究將充分考慮從全極化SAR圖像中提取其他的特征(如紋理和空間幾何信息等),并將其與極化特征相結(jié)合,對(duì)H/Alpha/A-Wishart分類器作進(jìn)一步改進(jìn)。
志謝:感謝歐洲空間局提供開(kāi)源軟件 POLSARPro。
[1] Cloude SR,Pottier E.A review of target decomposition theorems in Radar polarimetry[J].IEEE Transactions on Geoscience and Remote Sensing,1996,34(2):498-518.
[2] Cloude SR,Pottier E.An entropy based classification scheme for land applications of polarimetric SAR[J].IEEE Transactions on Geoscience and Remote Sensing,1997,35(1):68-78.
[3] Freeman A,Durden S L.A three-component scatteringmodel for polarimetric SAR data[J].IEEE Transactions on Geoscience and Remote Sensing,1998,36(3):963-973.
[4] Yamaguchi Y,Moriyama T,Ishido M,etal.Four-component scatteringmodel for polarimetric SAR image decomposition[J].IEEE Transactions on Geoscience and Remote Sensing,2005,43(8):1699-1706.
[5] Yamaguchi Y,Sato A,BoernerW M,et al.Four-component scattering power decomposition with rotation of coherency matrix[J].IEEE Transactions on Geoscience and Remote Sensing,2011,49(6):2251-2258.
[6] Zhang LM,Zou B,Cai H J,et al.Multiple-component scattering model for polarimetric SAR image decomposition[J].IEEE Geoscience and Remote Sensing Letters,2008,5(4):603-607.
[7] 王 慶,曾琪明,廖靜娟.基于極化分解的極化特征參數(shù)提取與應(yīng)用[J].國(guó)土資源遙感,2012,24(3):103-110.doi:10.6046/gtzyyg.2012.03.19.Wang Q,Zeng Q M,Liao J J.Extraction and application of polarimetric characteristic parameters based on polarimetric decomposition[J].Remote Sensing for Land and Resources,2012,24(3):103-110.doi:10.6046/gtzyyg.2012.03.19.
[8] Lee JS,GrunesM R,Ainsworth T L,etal.Unsupervised classification using polarimetric decomposition and the complex Wishart classifier[J].IEEE Transactions on Geoscience and Remote Sensing,1999,37(5):2249-2258.
[9] Pottier E,Lee JS.Application of the H/A/α polarimetric decomposition theorem for unsupervised classification of fully polarimetric SAR data based on the Wishart distribution[J].Proceeding of the Committee on Earth Observing Satellites SARWorkshop.Toulouse:NASA,1999:335-340.
[10] Lee JS,Grunes M R,Pottier E,et al.Unsupervised terrain classification preserving polarimetric scattering characteristics[J].IEEE Transactions on Geoscience and Remote Sensing,2004,42(4):722-731.
[11] Cao F,Hong W,Wu Y R.An improved Cloude-Pottier decomposition using H/α/SPAN and complex Wishart classifier for polarimetric SAR classification[C]//International Conference on Radar.Shanghai:IEEE,2006:1-4.
[12] 楊 磊,劉 偉,王志剛.加權(quán)全極化SAR圖像非監(jiān)督Wishart分類方法[J].電子與信息學(xué)報(bào),2008,30(12):2827-2830.Yang L,Liu W,Wang ZG.Weighted-based unsupervised Wishart classification of fully polarimetric SAR image[J].Journal of Electronics and Information Technology,2008,30(12):2827-2830.
[13] 郎豐鎧,楊 杰,趙伶俐,等.基于Freeman散射熵和各向異性度的極化SAR影像分類算法研究[J].測(cè)繪學(xué)報(bào),2012,41(4):556-562.
Lang F K,Yang J,Zhao L L,et al.Polarimetric SAR data classification with freeman entropy and anisotropy analysis[J].Acta Geodaetica et Cartographica Sinica,2012,41(4):556-562.
[14] 趙伶俐,楊 杰,李平湘,等.極化SAR影像弱散射地物統(tǒng)計(jì)分
類[J].遙感學(xué)報(bào),2013,17(2):312-310.Zhao L L,Yang J,Li P X,et al.Statistical classification of weak backscattering scatterers of PolSAR image[J].Journal of Remote Sensing,2013,17(2):312-319.
[15] Paladini R,Martorella M,Berizzi F.Classification ofman-made targets via invariant coherency-matrix eigenvector decomposition of polarimetric SAR/ISAR images[J].IEEE Transactions on Geoscience and Remote Sensing,2011,49(8):3022-3034.
[16] Zhang LM,Lu D,TangW Y.A GS-based built-up area detection method using polarimetric SAR images[C]//2012 IEEE International Geoscience and Remote Sensing Symposium.Munich:IEEE,2012:5911-5914.