一種改進的全極化SAR圖像MCSM－Wishart非監(jiān)督分類方法

陳 軍，杜培軍，2，譚 琨

(1.中國礦業(yè)大學國土環(huán)境與災害監(jiān)測國家測繪地理信息局重點實驗室，徐州 221116;2.南京大學江蘇省地理信息技術(shù)重點實驗室，南京 210023)

0 引言

全極化合成孔徑雷達(polarimetric synthetic aperture Radar，PolSAR)系統(tǒng)是在傳統(tǒng)的單極化、雙極化和多極化雷達系統(tǒng)基礎上發(fā)展起來的新型主動微波成像系統(tǒng)，如何通過全極化SAR系統(tǒng)提供的極化信息對地物進行分類已成為全極化SAR圖像應用的一個重要方面［1］。傳統(tǒng)的單極化和多極化SAR圖像處理通過提取紋理信息、極化信息等進行分類，然而，無論單極化還是多極化SAR系統(tǒng)獲取的都是部分極化信息;而全極化SAR系統(tǒng)包含同極化、交叉極化等不同的極化信息，可以全面反映地物目標的物理性質(zhì)。

國內(nèi)外學者在全極化SAR圖像極化特征提取方面開展了大量研究［2－7］。利用極化分解獲得的特征可以對全極化 SAR 圖像進行分類［8－15］。Lee等［8］首先用H/Alpha分解法提取的特征對圖像進行初始聚類，得到8個聚類中心;然后用描述多視協(xié)方差矩陣的復Wishart迭代分類器對圖像進行分類(簡稱 H/Alpha－ Wishart分類器)。Pottier等［9］提出了H/Alpha/A－Wishart分類器，在H/Alpha分解的基礎上加入A特征，將圖像聚成16類，再對圖像進行復Wishart迭代分類。Lee等［10］進一步結(jié)合Freeman分解法對上述方法進行改進，提出Freeman－H/Alpha－Wishart分類器，該分類器通過Freeman分解提取體散射、表面散射和二次散射信息，然后根據(jù)H/Alpha/A－Wishart分類器進行分類，分類中要求各類別必須保持純凈性(即各類別內(nèi)部各像元必須基于相同的散射性質(zhì))，必須是體散射、表面散射或二次散射中的一種，該分類器已經(jīng)隱含參考Freeman分解進行分類，但并沒有直接利用Freeman分解得到的特征，而是簡單地要求分類結(jié)果必須具有相同的散射性質(zhì)，因此還具有一定局限性。

無論H/Alpha－Wishart分類器、Freeman－H/Alpha－Wishart分類器或者 H/Alpha－SPAN［11］分類器，其分類結(jié)果均明顯依賴于H/Alpha初始聚類結(jié)果，導致最終分類精度不高。針對這一問題，本文提出一種改進的分類方法，該方法充分利用了全極化SAR圖像的表面散射、體散射、二面角散射、螺旋體散射和線散射等極化特征，在沒有人工干預的情況下能夠得到較好的聚類結(jié)果。

1 H/Alpha/A－Wishart分類方法及改進思路

1.1 分類方法

1.1.1 H/Alpha/A 分解與初始聚類

H/Alpha/A分解主要是提取全極化SAR圖像的散射熵H、散射角α(即Alpha，下同)和反熵A信息。其基本原理如下:

極化相干矩陣T可以表示為3個獨立的相干矩陣之和［2］，即

式中:λi為特征值(i=1，2，3);ei為矩陣T的第i個復特征向量;為ei的共軛矩陣。

在此基礎上，Cloude等［2］定義了熵H，用于表示地物散射的隨機程度，即

式中Pi為各特征值所占的百分比，即

散射角α代表二面角散射和表面散射的平均散射程度，即

反熵A表示優(yōu)勢散射以外的較弱散射分量間的關(guān)系［9］，即

根據(jù)H和α的不同取值，即可確定圖像每個像元在H/Alpha平面(圖1)中的位置，從而將像元歸入8個不同的類別;對每個類別取平均操作，得到初始聚類結(jié)果。加入A后，設置1個閾值，可以得到16個不同的類別。

圖1 H/Alpha 分類圖［2］Fig.1 Classification image of H/A lpha［2］

1.1.2 基于復Wishart分布的極大似然分類

以1.1.1節(jié)初始聚類結(jié)果作為中心，進行基于復Wishart分布的極大似然分類。其原理［8］如下:

極化協(xié)方差矩陣的概率密度分布函數(shù)服從n個自由度的復Wishart分布，其分布函數(shù)PT為

其中

式中:n為視數(shù);q代表不同的極化方式，在互易性情況下，q=3;Tr為矩陣的跡;C為極化協(xié)方差矩陣;∑=E(C);K為歸一化因子;Γ(·)為Gamma函數(shù)。對矩陣進行線性變換可以得到極化相干矩陣T。因此，極化相干矩陣的概率密度分布函數(shù)也服從n個自由度的復Wishart分布。

進一步對圖像進行開窗操作，計算窗口內(nèi)所有像元的平均極化相干矩陣;然后分別求該矩陣與8個類別中心的距離［8］，即

式中:∑i為第i個類別中心的極化相干矩陣;T為窗口內(nèi)所有像元極化相干矩陣的均值矩陣。比較8個距離，將像元歸為距離最小的那一類。

1.2 分類器存在的問題及改進思路

H/Alpha的初始聚類過程是根據(jù)圖1內(nèi)H和α的不同取值進行初始類別劃分的。而實際上，H和α取值和不同的地物類型并沒有嚴格的對應關(guān)系［8］。另外，該聚類在算法實現(xiàn)上也是采用近似的形式，沒有完全按照圖1內(nèi)的曲線進行邊界劃分，因此，由圖1得到的初始聚類結(jié)果存在一定的偏差。對于基于初始聚類中心的分類算法，初始聚類中心的位置直接影響最終分類結(jié)果，若該位置選擇不當，必然會降低分類精度。因此，對極化特征分解和初始聚類過程進行改進，有望提高全極化SAR圖像的分類精度。

2 改進的MCSM－Wishart分類算法

本研究從2個方面優(yōu)化初始聚類效果:①引入MCSM分解對全極化SAR圖像進行極化分解，提取表面散射、二次散射、體散射、螺旋體散射和線散射等更為全面的基本散射信息作為特征，解決僅使用H，α和A特征進行初始聚類的不足;②使用聚類效果更好的ISODATA聚類算法代替圖1中的簡單聚類算法，以提高初始聚類精度。

2.1 MCSM 分解

MCSM分解的極化協(xié)方差陣［6］為

式中:fs，fd，fv，fh和fw分別為表面(surface)散射、二次(double)散射、體(volume)散射、螺旋體(helix)散射和線(wire)散射系數(shù);α，β，γ和ρ為未知數(shù);α*，β*，γ*和 ρ*分別為 α，β，γ 和 ρ的共軛。求解式(9)，即可得到各散射系數(shù)及α和β的值;γ和ρ可以通過下式求解［6］，即

式中SHH，SVV，SHV和SVV分別表示極化方式為HH，VV，HV和VV的散射矩陣。

根據(jù)各散射成分的系數(shù)，即可計算各散射功率的值，計算公式［6］為

式中:Ps為表面散射功率;Pd為二次散射功率;Pv為體散射功率;Ph為螺旋體散射功率;Pw為線散射功率。

對極化協(xié)方差矩陣進行線性變換可以得到極化相干矩陣。因此，極化相干矩陣同樣可以通過MCSM分解得到以上各散射成分的系數(shù)。從當前極化目標分解算法的發(fā)展和評述來看，MCSM分解的方法能充分提取出各種自然和人工地物的各種散射成分的信息，將這些信息作為特征向量可以較全面地反映圖像的物理性質(zhì)［16］。

2.2 ISODATA 聚類

ISODATA是一種運行速度較快、可以自動合并和分裂聚類結(jié)果的非監(jiān)督分類或聚類方法。與K－means聚類算法相比，ISODATA可以進行類別的合并和分裂，靈活性較好，聚類結(jié)果更為合理。大量實驗分析表明，ISODATA聚類具有算法實現(xiàn)簡單、運算速度快、聚類效果較好等特點。

傳統(tǒng)的ISODATA聚類算法多采用歐式距離計算樣本與聚類中心的距離。但考慮到極化數(shù)據(jù)不同于光學數(shù)據(jù)的概率分布，本文對比歐式距離和馬氏距離后，采用馬氏距離進行距離計算，即

式中樣本si與中心mj分別為

2.3 MCSM－Wishart分類器改進

改進MCSM－Wishart分類器的核心思想是:①引入MCSM分解對圖像進行極化分解，提取表面散射、二次散射、體散射、螺旋體散射和線散射等更為全面的基本散射信息作為特征;②使用聚類結(jié)果更好的ISODATA算法。因此，MCSM－Wishart分類器的改進主要包括4個部分:圖像預處理、極化特征的提取、初始聚類結(jié)果的生成和復Wishart分布的迭代分類。具體分類步驟如下:

1)數(shù)據(jù)預處理。讀入全極化SAR圖像，進行多視處理，按照通用方法將多視比設為7∶1。設置濾波窗口大小為7像元×7像元，進行Refined Lee濾波，得到圖像的極化相干矩陣T。

2)MCSM分解。設置的窗口大小與步驟1)保持一致，對極化相干矩陣T進行MCSM分解，提取表面散射、二次散射、體散射、螺旋體散射和線散射等多種散射成分，將其組合成特征向量，用以表征極化信息。

3)ISODATA初始聚類。設置ISODATA聚類算法的初始條件(最小、最大類別等)，對步驟2)產(chǎn)生的特征向量進行迭代聚類，生成初始聚類結(jié)果。

4)計算聚類中心。在各聚類內(nèi)部進行平均操作，計算出初始聚類中心。

5)計算距離并聚類。選擇步驟4)生成的聚類中心，設置循環(huán)條件，進行基于復Wishart分布的迭代分類;根據(jù)式(8)計算圖像內(nèi)所有像元到各聚類中心的距離，將像元歸入距離最小的類別。

6)判斷是否滿足循環(huán)終止條件，若不滿足，返回步驟4)繼續(xù)進行計算;否則，完成分類。

3 實驗與分析

為驗證本文算法，選擇南京溧水和鹽城濱海研究區(qū)2景ALOS的1.1級PALSAR(phased array L－band SAR)圖像(全極化方式)的子圖像為數(shù)據(jù)源進行實驗。2景圖像的獲取時間分別為2011年3月22日和2009年4月9日。其中，研究區(qū)1(南京溧水)圖像大小為400像元×400像元，研究區(qū)2(鹽城濱海)圖像大小為500像元×500像元。研究區(qū)內(nèi)包括水體、建筑物、耕地、林地和道路等地物類型。實驗結(jié)果采用目視解譯數(shù)據(jù)和Google Earth圖像，通過綜合對比的方法確定測試樣例，進行精度評價。研究區(qū)1的Google Earth圖像成像時間為2013年4月10日;研究區(qū)2的TM圖像成像時間為2010年4月29日。2個研究區(qū)均屬縣級行政單位下的鄉(xiāng)鎮(zhèn)地區(qū)，全極化SAR圖像、TM圖像和Google Earth圖像在獲取時間間隔內(nèi)土地利用變化較小，因此雖然在時相上有一定差距，但仍具有參考價值。

實驗采用的處理軟件為POLSARPro(the polarimetric SAR data processing and educational tool)軟件，數(shù)據(jù)預處理包括多視處理、濾波和極化特征分解等，傳統(tǒng) H/Alpha/A－Wishart分類算法及監(jiān)督－Wishart分類算法采用該軟件自帶的功能，改進算法采用C語言實現(xiàn)。

研究區(qū)1和研究區(qū)2的Pauli基分解［1］(極化校準)及假彩色合成圖像如圖2所示。

圖2 研究區(qū)Pauli基分解及假彩色合成圖像Fig.2 False color com posite images and Pau li decomposition of two study areas

首先進行圖像預處理工作，包括多視、濾波和極化特征分解。對圖像進行多視處理，多視比設為7∶1;對多視處理后的圖像進行窗口大小為7像元×7像元的Refined Lee濾波處理;通過MCSM分解提取出各種散射成分。然后進行初始聚類，ISODATA算法設置的初始條件為:最小類別為8，最大類別為12，類內(nèi)最小個數(shù)為2，類間最小距離為4，最大迭代次數(shù)為20，最大類標準差為1。通過聚類結(jié)果計算聚類中心，最后進行復Wishart分類，得到分類結(jié)果。

作為對比，本文還分別進行了H/Alpha－Wishart和H/Alpha/A－Wishart算法的非監(jiān)督分類實驗，并選擇訓練樣本進行了監(jiān)督－Wishart分類實驗。其他分類器如 Freeman－H/Alpha－Wishart，H/Alpha－SPAN等的分類精度與H/Alpha/A－Wishart分類器的分類精度相似。為便于研究，僅以H/Alpha－Wishart和H/Alpha/A－Wishart算法進行實驗。為便于比較分類結(jié)果，結(jié)合訓練樣本和監(jiān)督分類結(jié)果對非監(jiān)督聚類后的結(jié)果進行了類別合并(圖3和圖4)。

圖3 研究區(qū)1不同算法分類結(jié)果Fig.3 Classification results by different algorithms for study area one

從圖3(a)可以看出，在研究區(qū)1，原H/Alpha－Wishart算法將水體和耕地混分。原因在于，當水體面積較少時，H/Alpha分解提取出的特征不足以反映地物的極化性質(zhì)，因此影響了最終的精度。雖然對H/Alpha/A－Wishart算法有所改善，但整體精度仍低于改進后的MCSM－Wishart算法精度。通過對比圖3(b)與圖3(c)可以看出，在研究區(qū)1，H/Alpha/A－Wishart算法存在耕地和林地混分、林地和建筑物混分的情況，而MCSM－Wishart算法則較少存在以上混分現(xiàn)象。另外，MCSM－Wishart算法提取出的水體也多于H/Alpha/A－Wishart算法，這是因為，MCSM算法提取出了更多的極化特征用于分類。

圖4 研究區(qū)2不同方法分類結(jié)果Fig.4 C lassification resu lts by different algorithm s for study area two

對比圖4(a)(b)(c)可以看出，在研究區(qū)2，H/Alpha－Wishart和H/Alpha/A－Wishart算法均不同程度地存在建筑和道路混分、道路和水體混分的情況;而MCSM－Wishart算法則能較好地區(qū)分出建筑、道路和水體，圖4下部的水體中混有道路的情況得到了改善，圖4右上部的建筑物中混有道路的情況得到改進。主要原因是MCSM方法提取出了螺旋體散射、線散射等較多的能反映城市地區(qū)地物類型的特征，因此能較好地對城市及周邊的不同人造地物進行類別劃分。

為比較各算法的分類精度，結(jié)合Google Earth圖像、TM圖像對2個研究區(qū)進行人工采樣，手工選擇了各種地物類型的測試樣本如表1所示。

表1 測試樣本Tab.1 Test sam p les (像元)

4種算法的分類精度見表2。

表2 4種算法的分類精度Tab.2 Classification results of four kinds of algrithm

從表2可以看出，無論研究區(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算法精度稍差，說明極化特征提取方法還有需要完善之處。4種算法的執(zhí)行效率見表3。

表3 4種算法的執(zhí)行效率Tab.3 Efficiency of four kinds of algorithm(s)

從表3中4種算法的運行時間來看，MCSMWishart算法效率高于H/Alpha/A－Wishart算法。雖然監(jiān)督－Wishart算法需要的時間更短，但這個時間并未包含選擇訓練樣本的時間、也沒有考慮訓練樣本選擇的實際困難。因此總的來看，MCSMWishart算法是一種有效的非監(jiān)督分類算法。

綜上所述，無論總體精度、Kappa系數(shù)還是運行效率，MCSM－Wishart算法均較 H/Alpha－Wishart算法和H/Alpha/A－Wishart算法有了明顯的提高。

4 結(jié)論

1)經(jīng)典的全極化SAR圖像非監(jiān)督分類主要以H/Alpha/A分解提取的極化特征為基礎進行圖像分類，但總體分類精度往往較低，不能滿足應用要求。造成這種情況的原因是多方面的，如該分解法并不能完全提取反映全極化SAR圖像的全部特征，選擇的分類器并不能很好地對圖像進行分類等。

2)引入MCSM算法是提高全極化SAR圖像分類準確性的有效途徑之一。本文在傳統(tǒng)的極化信息提取基礎上，引入多極化特征提取的MCSM算法，全面提取出反映城市地物類型的極化特征;然后通過改進的聚類算法進行初始聚類，以此為基礎對原H/Alpha/A－Wishart分類器進行改進。

3)以ALOS的PALSAR圖像為遙感數(shù)據(jù)源進行綜合分類實驗的結(jié)果表明，改進的MCSM－Wishart分類算法較經(jīng)典的H/Alpha/A－Wishart分類器在分類精度、Kappa系數(shù)、算法效率等指標上均有很大的提高。因此，提取全極化SAR圖像的線散射等極化特征，然后與二次散射、表面散射和螺旋體散射等經(jīng)典的極化特征相結(jié)合，輔以速度快、精度高的ISODATA聚類算法，可以有效地提高全極化SAR圖像分類的精度。

4)下一步的研究將充分考慮從全極化SAR圖像中提取其他的特征(如紋理和空間幾何信息等)，并將其與極化特征相結(jié)合，對H/Alpha/A－Wishart分類器作進一步改進。

志謝:感謝歐洲空間局提供開源軟件 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ù)提取與應用［J］.國土資源遙感，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］.電子與信息學報，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］.測繪學報，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］.遙感學報，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.