基于KummerU和MRF的極化SAR分類算法研究

石俊飛 林耀海 劉 璐

(1.西安電子科技大學(xué) 西安 710071;2.福建農(nóng)林大學(xué) 福州 350002)

0 引言

極化合成孔徑雷達(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)化。

1 KummerU 模型

極化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→∞時，

2 本文算法

2．1 Markov Random Field 模型

對一幅二維圖像，假設(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ù)，定義為:

2．2 算法描述

經(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。

3 實驗結(jié)果和分析

為了驗證該算法的有效性，一個真實的極化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é)果。

4 結(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.