石俊飛 林耀海 劉 璐
(1.西安電子科技大學(xué) 西安 710071;2.福建農(nóng)林大學(xué) 福州 350002)
極化合成孔徑雷達(dá)(SAR)圖像分類是極化SAR 圖像處理的重要任務(wù),對國防建設(shè),農(nóng)業(yè)發(fā)展都有很大的作用。最近,越來越多的學(xué)者開始關(guān)注極化SAR 圖像的分類。提出了很多極化SAR 數(shù)據(jù)模型和分類算法。例如,經(jīng)典的H/α 分類算法[1]是將極化SAR 數(shù)據(jù)分解為散射熵H 和散射角α,它們能夠反映地物的散射類型,通過H/α 平面將地物分為8 類。其他基于目標(biāo)分解的方法還有Freeman 分解[2]等。另外,一些數(shù)據(jù)模型被提出,從經(jīng)典的Wishart 分布[3]到后來發(fā)展的更高級的K 分布[4],G0 分布[5]和KummerU 分布[6]。相比于其他分布,KummerU 分布能夠更好的描述各種地物類型。
然而,基于這些分布的分類算法沒有考慮空間信息,使得同一地物容易受到斑點噪聲影響,得到椒鹽式的分類結(jié)果。馬爾科夫隨機場(MRF)[7]能夠很好的描述鄰域關(guān)系,因此,本文提出了基于KummerU 分布和MRF 的極化SAR 分類算法,該算法采用高級的KummerU 分布,能夠描述各類異質(zhì)地物。同時,使用MRF 模型加入鄰域信息,提高分類的區(qū)域一致性。期望最大化(EM)算法用來進行函數(shù)優(yōu)化。
極化SAR 數(shù)據(jù)通常用2 ×2 的S 矩陣表示,在Pauli 基下可以轉(zhuǎn)換為3 ×3 協(xié)方差矩陣C 和相干矩陣T。在相干斑一致性假設(shè)[8]下,S 滿足圓高斯分布,則C 矩陣服從Wishart 分布,定義如下:KummerU 分布退化Wishart 分布。因此,KummerU分布能描述各類型的數(shù)據(jù)模型。
其中,q 是通道數(shù)(一般情況下,q=3),n 是視數(shù),Σ 是平均協(xié)方差矩陣。
隨著雷達(dá)技術(shù)的發(fā)展,極化SAR 圖像分辨率越來越高,對于高分辨率極化SAR 圖像或異質(zhì)區(qū)域,如城區(qū)、森林等,相干斑一致性假設(shè)已經(jīng)難以滿足,因此,提出了非高斯的積模型[9],積模型假設(shè)協(xié)方差矩陣C 為兩個變量的積:一致的協(xié)方差矩陣Ch和紋理變量μ。即:
其中,Ch滿足Wishart 分布。
當(dāng)μ 為常數(shù)時,C 矩陣服從Wishart 分布,當(dāng)μ建模為gamma 分布時,C 矩陣服從K 分布,該分布能夠很好的描述森林等區(qū)域。當(dāng)μ 為逆Gamma 分布時,C 服從G0 分布,該分布適合極不勻質(zhì)區(qū)域,當(dāng)Gamma 中的參數(shù)變化時,G0 分布可以退化為K 和Wishart 分布。最近,一種新的分布被提取,該分布假設(shè)μ 服從Fisher 分布,定義如下:
其中L>0,M>0 為兩個形狀參數(shù),m 為尺度參數(shù)。Fisher 分布等于Gamma 和逆Gamma 分布的梅林卷積[10],因此,采用對數(shù)累積量的方法進行參數(shù)估計[11]。
根據(jù)Fisher 分布,得到C 矩陣服從KummerU 分布,公式如下:
KummerU 分布不僅能夠描述紋理結(jié)構(gòu),還能夠很好的描述勻質(zhì)區(qū)域。當(dāng)M→∞時,F(xiàn)isher 分布相當(dāng)于Gamma 分布,因此KummerU 分布退化為K 分布,當(dāng)L→∞時,F(xiàn)isher 分布相當(dāng)逆Gamma 分布,而KummerU 分布退化為G0 分布;當(dāng)L→∞且M→∞時,
對一幅二維圖像,假設(shè)觀測數(shù)據(jù)為y=,S 為位置集合,ys是在位置s 處的值,對應(yīng)的標(biāo)簽集合為,其中,cs ∈{1,…,L},L為類別個數(shù)。圖像分割可以看作求后驗概率。根據(jù)貝葉斯準(zhǔn)則,后驗概率可以表示為:
根據(jù)條件獨立性假設(shè),式(6)可以表示為:
其中,ηs為點s 的鄰域,cr為鄰域點r 對應(yīng)的類標(biāo)。
根據(jù)Gibbs 場和MRF 場的等價性[12],可以得到:
其中,U1(xs|cs)=-lnP(xs|cs)是觀測數(shù)據(jù)的能量函數(shù),是鄰域系統(tǒng)中可能基團Λ 的能量總和,因此,r∈ηs)的最大后驗概率就等于最小化如下能量函數(shù):
本文采用多層Ising模型來建模鄰域標(biāo)簽的分布,則平滑項能量函數(shù)可以寫為:
其中,β 是平滑項權(quán)重,它控制MRF 的平滑程度。δ(·)是Kronecker delta 函數(shù),定義為:
經(jīng)典的EM 算法用來進行能量函數(shù)優(yōu)化,H/α分類結(jié)果作為初始分割,因此,本文算法是一個無監(jiān)督的極化SAR 分類方法。具體步驟描述如下:
1)對極化數(shù)據(jù)進行7 ×7 的精致Lee 濾波;
2)采用H/α 分類算法進行初始分類,對每類估計KummerU 分布的紋理參數(shù)L,M;
3)EM 優(yōu)化迭代
4)用式(10)計算第j 個像素到第i 類的能量函數(shù)Uij,其中用式(5)計算數(shù)據(jù)項能量U1,用式(11)計算平滑項能量U2。
5)將每個j 像素賦給能量最小的類別,對每個像素重新分配類別。根據(jù)KummerU 分布,重新估計每類的紋理參數(shù)L,M,再回到4)。
迭代直到滿足停止條件,這里停止條件定義為迭代次數(shù)t≤iter。
為了驗證該算法的有效性,一個真實的極化SAR 圖像用來測試。該圖像是CONVAIR 衛(wèi)星拍攝的渥太華地區(qū)的全極化SAR 圖像,經(jīng)過10 視處理大小為222 ×342。對比算法為基于Wishart 的MRF方法[13]。實驗硬件條件為Intel(R)Core(TM)i3 CPU,4RAM。
圖1 給出了渥太華的全極化SAR 圖像偽彩圖,其中將Pauli 基作為RGB 三個通道顯示。從圖中可以看出該圖像中含有多種地物類型,左上角為城區(qū),城區(qū)內(nèi)部有道路,右下角有大片的裸地,其中有一些線目標(biāo)和樹木等點目標(biāo),各種地物尺度不一,形狀不同,因此,對其分類有一定的困難。
圖1(b)和(c)分別為基于Wishart 的MRF 方法和本文算法分類結(jié)果,從結(jié)果可以看出,本文算法在線目標(biāo)保持和邊界定位上能夠得到更好的結(jié)果。另外,在城區(qū)部分,本文算法能夠得到更好的結(jié)果。因為KummerU 分布能夠很好的描述異質(zhì)區(qū)。
圖1 渥太華地區(qū)極化SAR 圖像分類結(jié)果。
針對基于像素的分類結(jié)果對噪聲敏感的問題,本文提出了基于KummerU 和MRF 的極化SAR 分類方法。首先,對于復(fù)雜的地物類型,傳統(tǒng)的分布模型難以描述,因此,本文采用更高級的KummerU 模型對數(shù)據(jù)進行描述。另外,MRF 模型用來描述上下文關(guān)系,得到區(qū)域一致性更好的分類結(jié)果。
[1]Cloude S.R,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.
[2]L.Zhao,X.Zhou,Y.Jiang,and G.Kuang,Iterative classification of polarimetric SAR image based on the freeman decomposition and scattering entropy[C],in 1st Asian and Pacific Conference on Synthetic Aperture Radar,2007,473-476.
[3]J.-S.Lee and M.R.Grunes,Classification of multi-look polarimetric SAR data based on complex Wishart distribution[C],in National Telesystems Conference,1992,21-24.
[4]J.M.Beaulieu and R.Touzi,Segmentation of textured polarimetric SAR scenes by likelihood approximation[J],IEEE Transactions on Geoscience and Remote Sensing,2004,42:2063-2072.
[5]Niv.X,Ban,Y.An Adaptive Contextual SEM Algorithm for Urban Land Cover Mapping Using Multitemporal High-Resolution Polarimetric SAR Data[J],IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing,2012,5(4):1129-1139.
[6]Bombrun L,Vasile G,Gay M,Hierarchical Segmentation of Polarimetric SAR Images Using Heterogeneous Clutter Models[J],IEEE Transactions on Geoscience and Remote Sensing,2011,49:726-737.
[7]E.Rignot and R.Chellappa,Segmentation of polarimetric synthetic aperture radar data[J].IEEE Transactions on Image Processing,1992,1(3):281-300.
[8]LeeJS,Pottier E.Polarimetric radar imaging:from basics to applications[M].CRC press,2009.
[9]A.P.Doulgeris and T.Eltoft,Automated Non-Gaussian Clustering of Polarimetric SAR[C]//8th European Conference on Synthetic Aperture Radar(EUSAR),2010:1-4.
[10]Harant O,Bombrun L,Gay M,et al.Segmentation and Classification of Polarimetric SAR Data based on the KummerU Distribution[J].Proceedings of the Fourth International Workshop on Science & Applications of Sar Polarimetry & Polarimetric Interferometry Poiinsar,2009.
[11]S.N.Anfinsen and T.Eltoft,Application of the Matrix-Variate Mellin Transform to Analysis of Polarimetric Radar Images[J].IEEE Transactions on Geoscience and Remote Sensing,2011,49(6):2281-2295.
[12]HammersleyJM,Clifford P.Markov fields on finite graphs and lattices[J].1971.
[13]E.Rignot and R.Chellappa,Segmentation of polarimetric synthetic aperture radar data[J].IEEE Transactions on Image Processing,1992,1:281-300.