基于馬爾科夫判別譜聚類的極化SAR圖像分類方法

張向榮 于心源 唐 旭 侯 彪 焦李成

(西安電子科技大學人工智能學院，智能感知與圖像理解教育部重點實驗室，國際智能感知與計算聯合研究中心 西安 710071)

1 引言

極化SAR (PolSAR) 圖像地物分類在地質勘探、地形分析以及災害監測等方面具有廣泛的應用[1]。隨著理論水平的逐年提高，加上圖像處理領域的巨大需求，極化SAR圖像的分類問題成為了研究領域的熱點。

根據算法是否依賴于數據的先驗樣本，極化SAR圖像分類算法可以具體地分為有監督與無監督兩大類。無監督分類方法主要是利用極化數據的統計特征對同類型像素進行分類，不需要進行訓練，因此當訓練樣本嚴重不足的情況下，采用無監督分類方法的優勢十分明顯。且無監督方法的圖像分類過程簡單，充分利用了圖像的有效信息，適用范圍十分廣泛。無監督的極化SAR分類主要分為兩類：基于地物目標電磁散射特性和統計特性的分類法[2]以及基于聚類分析和圖像處理技術的分類法[3]。然而目前基于無監督的方法大多從極化數據的統計特性和散射特性角度出發，很少從圖像本身角度去考慮，不能充分利用極化數據的特征信息，無法全面地描述地面目標的物理屬性。因此，如何深度挖掘極化SAR圖像的特征信息，采用高效的處理方法提高分類精度，是當前極化SAR分類面臨的挑戰，同時也是本文重點關注的問題。

譜聚類在分析復雜的數據結構信息時，通過得到數據點的不同相似圖來預測聚類標簽，往往能夠顯示出其較強的聚類能力。將判別聚類方法應用到極化SAR圖像的分類中有兩個優點：首先，判別聚類可以將有監督的判別能力引入無監督分類問題中，現有的監督學習工具在弱監督任務和無監督任務中具有良好的重要性能；其次，判別聚類是一個通用框架，它允許將不同的判別損失函數或其他特定領域的約束合并到一個單一的損失函數中，適用于不同的應用程序，并且具有很強的靈活性[4]。

但是，現有的聚類方法只是簡單地將數據中的信息組合在一起，產生的噪聲會大大降低聚類的性能。為了解決上述問題，本文提出一種基于馬爾科夫的判別譜聚類方法(Markov Discriminative Spectral Clustering，MDSC)，具有低秩和稀疏分解的特點。在應用本方法時，首先需要構造一個原始的概率轉移矩陣，然后用其恢復一個真實的低秩概率轉移矩陣作為標準馬爾科夫譜聚類方法的關鍵輸入。為了能夠對極化SAR的數據信息進行多層次利用，本文在目標函數中引入了判別信息以提升聚類精度。

對于本方法中的目標函數的優化問題，在概率轉移矩陣上有一個低秩約束，同時在該矩陣的每一行上有一個概率單純形約束，本文提出了一種基于增廣拉格朗日乘子法的優化方法來解決這個有難度的優化問題。本文后續在各種實際數據集上進行實驗，結果表明，本方法具有較好的準確率，表現出了良好的分類性能。

2 基于馬爾科夫的判別譜聚類

2.1 標準馬爾科夫譜聚類框架

由隨機游走理論分析，當隨機游走到某個分類時，在該類中停留的概率較大，游走到其他類的概率較小。由此，Meila等人[5]提出，譜聚類可以用圖結構上的馬爾科夫隨機游動框架來描述。馬爾科夫鏈狀態簇是根據每個狀態到平穩態的距離進行聚類的，可以在拓撲圖上的隨機游走框架中連接，狀態轉移概率圖可以看作是一個有向圖。由此將譜聚類的圖譜理論應用于馬爾科夫鏈中，以達到聚類的目的。求解馬爾科夫隨機游動的轉移概率矩陣的特征值問題可以用來確定圖上的歸一化分割。本文所用的標準馬爾科夫譜聚類算法流程如下：

步驟1 計算所有數據點的相似度，構建相似度矩陣S；

步驟2 計算概率轉移矩陣P=D-1S以及它的平穩分布π=Pπ；

步驟3 構造拉普拉斯矩陣L=∏-1/2(∏P+PT∏)，其中∏表示對角元素為π(i)的對角矩陣；

步驟4 對L進行特征分解，得到前k個最小的特征向量；

步驟5 將這前k個特征向量作為矩陣的列向量構建特征矩陣；

步驟6 將特征矩陣的每一行作為數據點，利用k-均值算法對其進行聚類。

2.2 構造概率轉移矩陣

在馬爾科夫判別譜聚類方法中，最為關鍵的一步是如何構造一個精確的概率轉移矩陣。本文通過低秩和稀疏分解的方法，恢復真實的概率轉移矩陣，并將其用作標準馬爾科夫譜聚類方法的輸入，以獲得最終的聚類解決方案。

該方法的基本假設有兩個:

(1) 無向加權圖G的特征足以發現大部分聚類信息；

(2) 提取的特征可能會被噪聲破壞，即這些噪聲可能會導致一小部分數據點被分配到錯誤的類。

圖1 真實的概率轉移矩陣構造概圖Fig. 1 Real probability transfer matrix construction profile

在現實的譜聚類問題中，可以假設同一簇內任意兩點之間的轉移概率較高，而不同簇內兩點之間的轉移概率較低且近似為0，從而導致概率轉移矩陣的秩往往較低。綜上所述，根據這些觀察結果可以假設，反映潛在真實聚類信息的概率轉移矩陣往往是低秩的。

誤差矩陣E表示了P和之間的差異。根據假設，提取的特征足以識別大多數集群結構，所以P中的元素與中相應的元素只有一小部分顯著不同，可以說誤差矩陣E趨于稀疏。

綜上所述，在低秩稀疏假設下，可以將真實的概率轉移矩陣構造問題表示為，

由于原始概率轉移矩陣P的構建中只考慮了數據點之間的相似度，為了對數據信息進行多層次的充分利用，在目標函數中引入表示判別信息的判別損失函數Ec

其中，Ec表示極化 SAR 圖像分類的判別損失函數，1是全一向量，，表示非負平衡參數，跡范數是的秩在譜范數的單位球上的凸包絡，在實際問題中，最小化跡范數能夠得到理想的低秩結構[6]。表示的每一行都是一個概率分布，這強制保證了一定是一個概率轉移矩陣。

2.3 引入判別信息

判別聚類是一個將不同的判別損失函數或其他特定約束合并到一個損失函數中的通用框架。在進行聚類的過程中，可以通過分析不同類別的樣本信息，得到各自的特點與規律，進而構建出更加準確的判別準則對數據進行分類。本節將一種新的基于判別聚類的模型引入馬爾科夫譜聚類算法中，提高算法的信息利用率，改善分類精度。

本節根據SR模型[7]重新設計了判別損失函數，它結合了判別聚類項和正則項。前一項負責利用判別信息構建softmax損失函數，后一項負責降低由噪聲和異常值引起的過擬合。

綜上所述，在判別聚類的基礎上，將判別損失函數定義為

其中，W表示分類器矩陣，L(P，W|X)是softmax損失函數，R(W)表示正則項。引入softmax損失函數是為了解決由于不同類別的像素數不同而導致的樣本不平衡問題，該函數度量了分類器W與原始的概率轉移矩陣之間的一致性。用交叉熵來定義softmax損失函數L(P，W|X)

其中，e是自然常數，k表示類的個數。

引用正則項R(W)來降低由噪聲和異常值引起的過擬合。R(W)定義如下

判別損失函數Ec(P，W|X)雖然遵循SR模型的式，但本質上是不同的。在softmax分類方法中，訓練數據集中的ground-truth類標簽為常量。而在判別方法中，類標是在無監督算法下需要通過聚類得到的變量。

3 模型優化求解

在本節中將利用增廣拉格朗日乘子法[8]來解決上一節中低秩和概率單純形約束下的目標函數優化問題。與樸素拉格朗日方法相比，該方法提高了算法的魯棒性，并放寬了函數的強凸約束，使變換后的問題更易于求解。

下面根據 Xia 等人[9]提出的魯棒多視角譜聚類(Robust Multi-view Spectral Clustering， RMSC)方法，對，E，Ec進行更新求解。

首先假設Ec是已知的，式(2)對應的增廣拉格朗日函數為

其中，H是拉格朗日算子，μ>0是自適應懲罰參數。

(1) 求解E固定時，優化問題可以簡化為

利用奇異值閾值法[10]可以得到E的近似解

(3) 求解Ec：根據式(10)，由于可以通過迭代原始的概率轉移矩陣P得到，因此可以通過固定的P來最小化分類器W，并使用迭代優化算法來解決這個問題。

對Ec進行求導，梯度可計算為

根據這個求導式，使用L-BFGS優化算法[10]來最小化式。

綜上所述，真實的概率轉移矩陣構建算法框架如下：

步驟2 根據式(8)更新E；

步驟5 使用L-BFGS優化算法更新Ec;

得到真實的概率轉移矩陣后，即可按照表標準馬爾科夫譜聚類的算法流程，得到極化SAR圖像的最終分類結果，如圖2所示。

圖2 本文算法框架圖Fig. 2 Algorithm frame diagram

4 實驗結果與分析

4.1 實驗數據與環境說明

由于本文算法是在馬爾科夫譜聚類算法基礎上的改進模型，對目標函數的不足進行優化，并引入了判別信息。因此在實驗中選取3個對比算法如下：

(1) 聯合正則譜聚類(Co-Regularized spectral clustering， Co-Reg)[11]：譜聚類的共正則化方法，Kumar于2011年提出；

(2) 混合馬爾科夫鏈(Mixture of Markov Chains， MMC)[12]：Zhou和Burges于 2007年提出的混合馬爾科夫鏈方法，這是與本文所提基于馬爾科夫鏈的判別譜聚類算法最相關的方法；

(3) SR-MO算法[13]：Haixia Bi提出的無監督判別聚類方法，利用監督Softmax 邏輯回歸 (Softmax logistic Regression， SR)模型和大量特征進行無監督分類并在分類過程中采用了馬爾可夫隨機場優化算法(Markov random field Optimization， MO)，且考慮了空間關系。

本節在荷蘭Flevoland地區小農田和大農田、德國Oberpfaffenhofen地區和西安地區這4幅真實的極化SAR數據上進行實驗，以上4種數據分別來自不同的成像系統，包含不同的波段與數據類型，通過以上數據證明本文算法的有效性。仿真實驗均是在主頻2.50 GHz的Intel(R) Core(TM) i5-7300HQ CPU， 8 G的內存環境和Windows10操作系統中編程實現的。實驗結果均為MATLAB R2017a的軟件環境中進行10次實驗的平均值。

本文實驗主要用總體分類精度OA、平均分類精度AA以及Kappa系數作為分類的評價指標。

4.2 荷蘭Flevoland地區農田小圖實驗結果

本實驗數據為NASA在1989年使用AIRSAR系統獲得的荷蘭Flevoland地區L波段的農田小圖數據，該組圖像的大小為300×270。圖像主要包含裸土、馬鈴薯、甜菜、大麥、豌豆、小麥6種農作物。

Co-Reg， MMC以及SR-MO 3種不同對比算法和本文算法對Flevoland地區小農田圖的總體分類精度OA、平均分類精度AA以及Kappa系數如表1所示；分類結果圖如圖3所示。

在圖3中，圖3(a)是荷蘭Flevoland地區農田小圖的Pauli分解偽彩圖，圖3(b)是地物類標圖，圖3(c)-圖3(e)分別是Co-Reg， MMC以及SR-MO算法的分類結果圖，圖3(f)是本文算法的分類結果圖。

從總體分類精度和平均分類精度來看，本文方法均為最高且總體分類精度達92.43%，分別比Co-Reg，MMC和SR-MO算法高出8.43%， 5.35%和1.13%。雖然MMC在馬鈴薯這一類的分類效果最好，對豌豆、小麥等大多數地物類別也更加準確，但這也導致了將部分甜菜等地物錯分為馬鈴薯，影響分類精度。這也說明了與單純馬爾科夫譜聚類算法相比，本文提出的引入判別信息的馬爾科夫譜聚類算法更有優勢。并且在分類結果圖中，可以看出本文算法對不同的地物分類比較均衡，邊緣更清晰，孤立像素更少，顯示出算法的平滑效果，Kappa系數也是最高的，驗證了本文方法的有效性。

4.3 德國Oberpfaffenhofen地區實驗結果

本實驗數據為德國國家宇航中心DLR使用ESAR系統拍攝的德國Oberpfaffenhofen地區L波段極化SAR數據的局部，400×450，分辨率為3×2.2 m。該區域主要分為農田、居民區、林地、道路和其它地物5類。

Co-Reg， MMC和SR-MO 3種不同對比算法和本文算法對德國Oberpfaffenhofen地區的總體分類精度OA、平均分類精度AA以及Kappa系數如表2所示，分類的結果圖如圖4 所示。

在圖4中，圖4(a)是德國Oberpfaffenhofen地區的Pauli分解偽彩圖，圖4(b)是地物類標圖，圖4(c)-圖4(e)分別是Co-Reg， MMC以及SR-MO算法的分類結果圖，圖4(f)是本文算法的分類結果圖。

從圖4中能夠看出，圖4(c)的雜點最多，區域一致性最差。圖4(d)整體分類效果較好，但和本文算法相比，對邊緣像素點的分類效果不太理想。圖4(e)和圖4(f)相比，圖4(e)將林地和開放型區域錯分的像素點較多。圖4(d)錯分的雜點較多，將大部分林地區域錯分為開放型區域。對比算法的實驗結果圖中，圖4(e)分類效果最好。

由實驗結果表2可以得到本文算法的總體分類精度為79.74%，分別比Co-Reg， MMC和SR-MO方法高出6.11%， 5.03%和1.52%。Kappa系數和平均分類正確率也優于對比方法。這說明了本文算法對德國Oberpfaffenhofen地區的分類效果較好，地物之間的分界較為清晰，能夠識別出農田和道路，以及大部分的林地區域。同時，可以看到本文算法在處理除過道路之外的區域時，分類效果很好，且更加穩定，尤其在圖中間的農場區域，分類更為連貫平整，視覺效果好，雜點較少。Co-Reg方法對于道路的誤分現象較為嚴重，大部分道路沒有被識別。SR-MO方法能夠有效識別大部分的道路、郊區、林地，以及農田，卻將大部分的農田誤分為郊區，分類效果也不夠理想，這是由于僅僅基于判別信息進行分類，沒有反映極化數據的統計特征。不過這4類算法對中間道路的分類效果都不是很好，還有很大的改進完善空間，但相對來說本文算法的道路邊界更為清晰平滑，視覺效果更好。

表1 4種算法對Flevoland小農田圖的分類結果Tab. 1 Classification results of four algorithms for Flevoland small farmland map

圖3 荷蘭Flevoland地區農田小圖的偽彩圖、類標圖以及不同算法的分類結果圖Fig. 3 Pseudo-color map， class diagram and classification results of different algorithms for farmland maps in the Flevoland region of the Netherlands

表2 4種算法對德國Oberpfaffenhofen地區的分類結果Tab. 2 Classification results of four algorithms for the Oberpfaffenhofen region of Germany

4.4 西安地區實驗結果

本實驗數據為由加拿大太空署RADARSAT-2系統獲取的西安地區極化SAR圖像，該圖像大小為512×512，主要有河流、城區、植被3種地物分類。

Co-Reg， MMC以及SR-MO 3種不同對比算法和本文算法對西安地區的總體分類精度OA、平均分類精度AA以及Kappa系數如表3所示；分類的結果圖如圖5所示。

在圖5中，圖5(a)是西安地區的Pauli分解偽彩圖，圖5(b)是地物類標圖，圖5(c)-圖5(e)分別是Co-Reg， MMC以及SR-MO算法的分類結果圖，圖5(f)是本文算法的分類結果圖。

由表3可以看出本文算法的總體分類精度為85.03%，分別比Co-Reg， MMC和SR-MO高出11.03%， 8.33%和2.76%，尤其在城區和植被的分類中，本文算法均表現出了良好的性能。從圖5的視覺效果上分析，Co-Reg和MMC算法的分類效果較差，通常情況下某一區域內的樣本點應屬于同一類地物，但這兩幅結果圖的整個圖像充滿斑點點，區域內雜點過多。相比之下，SR-MO和本文算法的兩幅圖視覺效果很好，雜點較少，河流區域分類較好。相比與SR-MO算法，本文算法的城區部分分類較好，能較好地保持區域一致性。

圖4 德國Oberpfaffenhofen地區數據的偽彩圖、類標圖以及不同算法的分類結果圖Fig. 4 Pseudo-color map， class diagram and data classification results of different algorithms in the Oberpfaffenhofen region of Germany

表3 4種算法對西安地區的分類結果Tab. 3 Classification results of four algorithms for Xi'an area

4.5 荷蘭 Flevoland 地區大農田地區實驗結果

本節實驗美國NASA/JPL AIRSAR系統于1989獲得的Flevoland地區四視L波段的大圖數據，圖像大小為750×1024，分辨率為12.1×6.7 m。包含15類地物：蠶豆、油菜籽、裸地、土豆、甜菜、小麥2、豌豆、小麥3、苜蓿、大麥、小麥、草地、森林、水域和建筑物。設置Ns=15，K=9， 4種不同算法的總體分類精度OA、平均分類精度AA以及Kappa系數如表4所示，不同算法的結果圖如圖6所示。

從圖6中能夠看出，圖6(c)的雜點最多，區域一致性最差，油菜籽、甜菜和小麥等區域被大量誤分為水域，區域之間沒有明顯區分，分類效果較差。圖6(e)整體分類效果較好，但和本算法相比，對邊緣像素點的分類效果不太理想。和圖6(f)相比，圖6(d)與圖6(e)區域錯分的像素點較多，小麥、草地和建筑物等區域均有較多雜點。對比算法的實驗結果圖中，圖6(e)分類效果最好。

與對比算法相比，本文算法對Flevoland地區大農田的農作物分類的正確率都很高，穩定性很好，沒有偏差，尤其在油菜籽、土豆、苜蓿、森林這幾類上，分類正確率明顯優于其它3種算法。但是本章算法在處理草地和小麥區域時，最終結果并沒有達到理想狀態，還有進一步提升的空間。但是整體來說，本文算法對比于其他對比算法的分類效果最好，區域一致性較好，邊界更清晰，錯分點更少。

圖5 西安地區數據的偽彩圖、類標圖以及不同算法的分類結果圖Fig. 5 Pseudo-color map， class diagram and data classification results of different algorithms in Xi'an area

表4 4種算法對荷蘭 Flevoland 地區大農田圖的分類結果Tab. 4 Classification results of four algorithms for large farmland maps in the Flevoland region of the Netherlands

4.6 參數敏感性分析

在本文方法中有兩個權衡參數λ，β和正則項參數ξ。通常的做法是在無監督聚類中根據經驗設置參數。下面本節對荷蘭Flevoland小農田、德國Oberpfaffenhofen和西安地區數據集進行實驗，觀察不同值和對總體分類精度的影響，如圖7-圖10所示。可以觀察到：

綜上所述，本文算法在合適的區間內對其參數相對不敏感，區間內參數的變化對分類總精度影響較小，且參數的可調節的范圍較大。因此，本文所提方法有較好的參數穩定性。這使得本文算法易于使用，無需進行太多的權衡參數調優。

5 結束語

本文提出一種基于馬爾科夫的低秩稀疏的判別譜聚類方法。首先構造一個概率轉移矩陣用于恢復一個真實的低秩轉移概率矩陣作為標準馬爾科夫聚類方法的關鍵輸入。然后在目標函數中引入判別信息，達到對數據信息的充分利用。本文采用基于增廣拉格朗日乘子法的優化方法來求解低秩和概率單純形約束下的目標函數。通過應用4種典型的實驗數據，證明了本文算法在分類精度、參數敏感性等方面具有優勢，最終的分類效果也更好。

圖6 荷蘭 Flevoland 地區大農田數據的偽彩圖、類標圖以及不同算法的分類結果圖Fig. 6 Pseudo-color map， class diagram and classification results of different algorithms for large farmland data in the Flevoland region of the Netherlands

圖7 荷蘭小農田中不同和下的分類結果圖Fig. 7 Classification results of different and below in small Dutch farmland

圖8 德國地區中不同和下的分類結果Fig. 8 Classification results for different and below in the German region

圖9 西安地區中不同和下的分類結果圖Fig. 9 Classification results of different and subordinates in Xi'an area

圖10 不同正則項參數的分類結果圖Fig. 10 Classification results of different regular item parameter