基于超像素MRF的農田地區高分遙感影像分割

蘇騰飛， 張圣微， 李洪玉

(內蒙古農業大學水利與土木建筑工程學院，呼和浩特 010018)

0 引言

從遙感影像中快速、準確地提取農田信息，是農業遙感的基礎研究內容。近年來，很多學者都指出利用高空間分辨率遙感影像(high resolution remote sensing image，HRI)開展農田監測的必要性[1-3]，其中一個重要理由是： 很多地區，尤其是中國的一些區域，農田平均面積較小[4]； 如果利用中等空間分辨率的遙感影像，如空間分辨率為250 m的MODIS或30 m的Landsat數據，極易因像元混合導致錯誤的分割結果。因此，發展適用于HRI的農田信息提取算法是非常必要的。

HRI可以捕獲眾多地物的細節信息，使得很多地物在影像中以面目標的形式呈現，這就使HRI像素之間產生了較高的空間關聯性。如何有效利用像素之間的空間信息，是發展HRI解譯算法的關鍵。馬爾科夫隨機場(Markov random field，MRF)作為一種影像分割算法模型，在其Gibbs勢能函數中直接將像素的光譜和空間信息相關聯，以提升像素類別場在空間上的連續性，因而受到很多學者的青睞。

然而，基于MRF影像分割算法的求解通常具有較高的計算復雜度。主流的MRF模型求解都采用了Gibbs勢能函數最小化的策略[5-9]，通過逐像素的迭代計算，試圖求出勢能函數最低時所對應的分割結果。為了提高MRF模型的求解效率，很多學者將超像素的圖模型引入到MRF模型中，以最大限度地減少MRF的求解空間，從而提高計算效率。孫巍等[7]將簡單線性迭代聚類(simple linear iterated clustering，SLIC)超像素生成的結果用于條件隨機場的建模； 劉磊等[8]將高階勢能引入到超像素MRF算法中； Yu等[9]借助區域生長的策略，也發展了一套基于超像素的MRF分割算法； Qin等[10]將Yu的算法拓展到多波段影像的分割中； 蘇騰飛等[11]將全局最優合并的策略引入到Yu的算法中。

雖然超像素MRF模型具有求解效率高的優勢，但目前主流的超像素生成算法難以將超像素的邊界與影像中地物的邊界完全地吻合，而這類錯誤往往增大了MRF求解的不確定性。超像素算法的目的是將影像分割為大小近似相等、形狀大致相似的若干區域(或超像素)[12-14]。目前主流的超像素生成算法，例如SLIC[12]，Mean-Shift[13]，TurboPixel[14]等，雖然都可以取得較好的超像素分割結果，但是也難以保證超像素邊界與影像中地物的實際邊界完全符合，這類現象會導致超像素MRF算法在交互勢函數最小化的過程中發生錯誤。近年來雖然出現了較多基于超像素MRF算法的研究，但是都忽略了這一問題[7-11]。

在農田HRI的分割中，邊界信息是非常重要的。很多相鄰的農田地塊在HRI中具有十分相似的光譜特征(例如種植同類作物的農田)，而它們之間間隔的通常是灌渠或道路，在HRI中對這類地物的解譯主要依靠的是邊界信息[15]。因此，將超像素MRF算法應用于農田HRI分割時，如何最大限度地利用邊界信息，并且避免超像素生成算法引入的邊界不完全符合的問題，對于分割精度的提高將是非常關鍵的因素。

針對農田地區的HRI，提出了一種基于超像素MRF的HRI分割算法。為了避免超像素邊界不完全符合造成的誤差對分割精度產生的負面影響，提出一種超像素的高階鄰域模型，并重新定義了超像素MRF的交互勢函數，以全面利用邊界信息來提高分割精度。在本文模型的基礎上，還給出了一種交互勢函數調節參數的自動估計方法，以進一步提高精度。

1 算法原理

1.1 超像素MRF模型

遙感影像可被看作一個矩形格陣：S={s(i，j)|1≤i≤H，1≤j≤W}，其中i和j分別代表像素的列、行號，W和H分別為影像的總列、行數。在S的基礎上，定義2個隨機場X=(Xs)s∈S與Y=(Ys)s∈S。X為標號場，其元素x的值取自整數集合L={l|1≤l≤L′}；Y為觀測場，其元素y由像素光譜值表示。在基于超像素的MRF模型中，S被劃分為在空間上互不重疊的若干區域(超像素)：R={rt|1≤t≤M}，rt={(i，j)|l(i，j)=t，(i，j)∈S}。其中每個超像素rt包含若干像素，并且所有rt的并集等于S，M為總區域數。

在遙感影像分割問題中，已知的是影像的像素值，求解的是各個像素的類別。MRF將該問題建模，為求取P(x)P(y|x)最大化的問題，一般將其進行負對數變換，即

U(x|y)=-ln(P(x|y))=-ln(P(x))-ln(P(y|x))=U(x)+U(y|x)，

(1)

進而轉化為勢能函數最小化的問題，即

(2)

式中：P(x|y)是X在Y條件下的條件概率分布，也是MRF模型的目標函數；P(x)為X的先驗分布；P(y|x)是Y在X條件下的條件概率分布；U(x)為一階勢函數，反映了X與Y在概率分布上的關聯性；U(y|x)為交互勢函數，表示X與Y的元素在空間分布上的關系； β是2個勢能函數項的調節參數，需要根據具體影像來設置或估算。

一階勢函數一般采用高斯分布模型來建模[9-11]，其計算方法為： 對于某個區域i，其屬于類別c的勢能可由下式估計，即

(3)

在MRF迭代求解的初期，馬氏距離常因協方差矩陣的估計偏差而產生較大的錯誤。因此本文在具體實現中采用了范數距離取對數的方式，即

(4)

式中： ||·||表示2階范數；Z2同樣為歸一化常數，與Z1類似。取自然對數的目的是為了減少范數距離的數值變化范圍。

1.2 高階鄰域交互勢函數建模

傳統的交互勢函數主要利用的是邊界強度等影像的空間背景信息，引導求解過程向正確的方向收斂，并保證在空間上相鄰的像素或超像素具有一致的類別，因此該項又被稱為類邊界懲罰函數。基于超像素MRF的交互勢函數的一般形式為

(5)

式中：rs與rt表示2個相鄰的區域；g(·)表示2個區域公共邊界的邊界強度均值；xs與xt分別為2個區域的類別標號。

一般而言，傳統計算方法僅考慮了2個相鄰區域的公共邊界的像素[7-10]。假如超像素過分割的結果與實際地物的邊界完全符合(如圖1(a))，則傳統計算方法可以準確地反映2個區域之間的邊界強度信息。但在很多實際情況中，超像素生成算法需要兼顧超像素內部灰度上的一致性和超像素之間在形狀上的相似性，這使得超像素之間的邊界像素難以完美地與實際地物邊界符合(如圖1(b))，此時傳統計算方法難以準確反映2個區域的邊界強度信息，圖1(a)和(b)中的顏色反映的是地物的光譜值。

(a) 與實際地物邊界完全符合(b) 與實際地物邊界有誤差 (c)r1和r2的高階鄰域范圍

圖1超像素邊界示意圖

Fig.1Illustrationofsuper-pixels’boundary

為了解決以上問題，本文提出了一種高階鄰域的交互勢函數模型。該模型考慮了超像素之間高階鄰域像素的邊界強度值，以更準確地描述地物邊界信息，從而提高MRF的分割精度。該模型與式(5)類似，但增加了一個高階鄰域范圍參數hn，即

(6)

與式(5)不同的是，式(6)中計算方法考慮了2個超像素邊界像素的高階鄰域，如圖1(c)，淺灰色表明了hn為3時高階鄰域范圍內的像素，在此基礎上計算式可表示為

(7)

式中：D(r1，r2)是超像素r1與r2公共邊界的像素集合；N(p，hn)是像素p以hn為范圍的高階鄰域集合像素；p和m分別表示屬于集合D與N的某一像素；qm表示像素m處的邊界強度值，其范圍在(0，1)之間；d(r1，r2)表示超像素r1與r2公共邊界的長度，以像素為單位； 系數3為經驗常數，是根據多次實驗確定的。本文的邊界強度計算方法采用了文獻[10]的矢量梯度方法，該方法可以從多波段影像中提取單波段的邊界強度影像。

1.3 基于邊界懲罰的β參數估計方法

β參數的作用是調節一階勢函數與交互勢函數對MRF總勢能函數的貢獻。β可由用戶根據實際情況來設置，但此方法受主觀因素影響較大，且往往需要多次實驗，費時費力。β還可根據影像中的統計特征量來估計[9-11]，但該方法需要較多的類別樣本，而且對于一景影像的分割，該方法的β是固定的。根據式(5)和(6)，可知交互勢函數主要懲罰的是類別不同區域的邊界，而不同類別區域的邊界具有不同程度的邊界強度信息，因此在類別邊界懲罰中應考慮到具體類別的統計信息。基于該思路，本文提出了一種自動的β估計方法，即

β(xs,xt)=ln(‖μs,t‖)，

(8)

式中μs,t為類別xs與xt的光譜均值向量相減所得的向量。需要注意的是，其計算結果與式(4)中的對數范數距離類似，因此，式(8)在考慮到具體類別統計信息的同時，也在數值上使交互勢函數和一階勢函數更好地相互影響，從而更為精確地引導MRF分割算法的求解過程。

1.4 求解方法與算法流程

很多主流的MRF模型求解算法采用的是迭代條件模式(iterated conditional mode，ICM)[5-6]。該方法是在最大先驗概率(maximum a priori，MAP)的理論上建立的。ICM通過逐像素的勢能函數最小化迭代計算，試圖求出Gibbs勢能函數最低時所對應的分割結果。在本文超像素MRF模型的基礎上，給出基于ICM的求解方法及具體流程。

輸入包括原始HRI、單波段的邊界強度影像、超像素生成算法的參數、類別數目參數c、高階鄰域范圍參數hn和最大迭代次數Tmax。

具體流程如下：

1)利用超像素生成算法，對輸入的HRI進行過分割。

2)利用最近鄰域法，初始化各個超像素的類別。

3)設置迭代次數T=0，并進行迭代計算。①估計各個類別的光譜均值向量μ； ②根據式(8)計算出不同類別邊界的參數β； ③利用式(2)計算各個超像素的類別； ④T加1。

4)若T≤Tmax，返回步驟3。

5)結束并輸出結果。

由算法流程可知，在算法初始階段，需對各個超像素的類別進行初始化。該步驟可用非監督聚類與大多數投票的方法來實現[11]，即先利用K-Means等算法給出粗略的聚類，再將每一個超像素類別初始化為其中大多數像素所屬的類別。然而在實驗中，該方法易受到K-Means初始化的影響而使后續步驟產生較差的分割效果。

因此本文利用了一種簡單、高效的半監督策略來完成超像素類別的初始化，該方法包含3個步驟： ①用戶在原始影像中選取樣本，每一個類別僅選取5～10個具有代表性的像素作為樣本； ②利用樣本計算各個類別的光譜均值向量，作為每個類別的中心； ③根據各個類別的中心向量，利用最近鄰域法對各個超像素進行類別初始化。實驗表明，該方法可在用戶干預較少的前提下保證較高的精度。

步驟3和4實際上是利用ICM實現的。值得一提的是，步驟3中第③步，在估計超像素類別的計算時，將式(4)，(6)和(8)帶入式(2)，以估計本文超像素MRF模型的勢能函數最小時所對應的超像素類別。

2 實驗驗證

2.1 實驗數據

為了驗證本文算法的性能，共選取了2景農田地區的高空間分辨率多光譜遙感影像，以開展分割實驗。

本文進行實驗的計算機配置為： Windows7，4 G內存，CPU： Intel I5 4200 m(2.5 GHz)。2景影像分別由中國的資源三號衛星和美國的OrbView-3衛星獲取，以下對其分別簡稱為S1和S2，其基本信息見表1。

表1 實驗影像數據基本信息Tab.1 Basic information of the image data used for experiment

2景影像均包含近紅外、紅光、綠光和藍光4個波段。由于原始影像幅寬巨大，而本文實驗的主要目的是驗證算法對于農田地區的分割精度，因此在S1和S2中選取了具有代表性的子影像開展實驗，子影像大小均為400像素×400像素，采用近紅外(R)，紅光(G)，綠光(B)波段假彩色合成影像，如圖2所示。

(a) S1子影像(b) S2子影像

圖2實驗用的子影像

Fig.2Subsetsforexperiment

由圖2可見，S1和S2分別展現了不同地貌特征的農田地區。S1覆蓋了內蒙古錫林郭勒盟草原地區的農田，主要作物包括玉米和小麥，且當地的農業主要依賴雨水和河水灌溉。該區的農田平均面積較大，但農田內部的光譜變化也較為明顯。另外，該區的草場由于牧草收割活動而呈現出明顯的草場區域邊界。S2為寧夏中衛市的郊區區域。當地屬半干旱氣候，依賴黃河水灌溉，盛產小麥和葵花等經濟作物。相比S1，S2中農田的平均面積較小，但其農田內部的光譜均勻性更高。值得一提的是，S2中存在較多的休耕田地，且其內部光譜由于雜草的存在而呈現出較大變化。

為了開展算法精度的定量評價，以專家手動提取的分割結果作為基準，計算不同MRF算法分割結果的總精度(overall accuracy，OA)與Kappa系數進行定量對比分析。

2.2 參數hn對算法性能的影響

為了進一步分析本文算法性能(分割精度與計算速度)與高階鄰域參數hn的關系。研究了S1和S2分割實驗的運算時間、OA和Kappa隨hn的變化規律，如圖3所示。hn的范圍被設置為[1，15]。值得一提的是，當hn=1時，高階鄰域交互勢函數與普通交互勢函數的計算相同。

(a) S1子影像OA與hn(b) S2子影像OA與hn

(c) S1子影像Kappa與hn(d) S2子影像Kappa與hn

(e) S1子影像運算時間與hn(f) S2子影像運算時間與hn

圖3本文算法性能與hn的關系

Fig.3Relationshipbetweentheperformanceoftheproposedalgorithmandhn

在分割精度方面，S1和S2子影像的最高精度均是在hn>1時得到的，這較為有力地證明了本文提出的高階鄰域交互勢函數對于提升分割精度的作用。對于S1子影像，hn=3時精度達到最高，而S2子影像最高精度所對應的hn為6。這表明，對于不同的影像，其最佳hn的取值可能存在差異。另外一個較為明顯的規律是，隨著hn的增加，雖然S1和S2子影像的精度起初增加并達到最高值，但隨后精度呈下降趨勢，且在hn較大時S1和S2子影像的精度均低于hn=1時的精度。這在一定程度上說明，超像素邊界像素的高階鄰域在一定范圍內對提高分割精度具有積極作用，但是過大的高階鄰域范圍反而會誤導MRF的求解過程，以致精度降低。需要說明的是，本文高階鄰域交互勢函數提出的目的，是為了降低超像素與實際地物邊界不完全符合造成的誤差對MRF求解的負面影響。在此前提下，hn應近似等于超像素分割與實際地物邊界的誤差。因此，過高或過低的hn均不利于邊界強度等空間背景信息在MRF求解過程中起到積極作用。

在計算速度方面，S1和S2子影像均顯示出了同一規律： 運算時間隨hn的增大而變長。有趣的是，雖然S1和S2子影像的像素數目相同，超像素的數目也相近，但S2的運算時間明顯長于S1。這主要是由于： 根據式(6)，只有在相鄰超像素類別不同時，高階鄰域交互勢函數的計算才會進行； S1中的地物平均面積較大，因此相鄰超像素之間類別不同的情況較少，而S2中很多地物的面積較小，這使得相鄰超像素之間類別不同的情況更多，因此交互勢函數的計算次數更多，導致耗時更久。

2.3 算法對比分析

為了全面、客觀地驗證本文算法的性能，選取了另外3種超像素MRF算法，以開展對比實驗。算法分別為： ①無β參數自動估計的本文算法(以下簡稱M1)，該算法也采用了高階鄰域交互勢函數，且其hn的設置與本文算法相同； ②一種高階MRF影像分割算法[8](以下簡稱M2)； ③一種融合了MRF與區域生長的多波段影像分割算法[10](以下簡稱M3)。本文算法簡稱為M0。進行M0與M1對比的目的是驗證本文提出的β參數自動估計方法對分割精度的影響。

為了消除初始化對分割結果的影響，4種算法均采用SLIC生成的超像素，其結果見圖4(a)和(b)。另外，4種算法求解過程中所需的邊界強度影像，也都由同一方法計算[10]，見圖4(c)和(d)。此外，4種算法在分割S1和S2子影像時，類別數目參數c均分別設置為5和4； 分割S1和S2子影像時的hn分別被設置為3和6；Tmax均設置為50，以平衡分割精度與收斂速度。

(a) S1子影像超像素分割(b) S2子影像超像素分割(c) S1子影像邊界提取結果 (d) S2子影像邊界提取結果

圖4S1和S2子影像超像素分割與邊界強度提取結果

Fig.4Super-pixelsegmentationresultsandedgestrengthextractionresultsofS1andS2subsets

S1子影像4種不同MRF分割算法的分割結果如圖5所示。

(a) M0 (b) M1 (c) M2 (d) M3

圖5S1子影像4種MRF算法的分割結果

Fig.5S1subsetsegmentationresultsof4differentMRFalgorithms

觀察圖5(a)可以發現，本文算法M0得到了較好的分割結果，其中大部分地物內部的類別具有較好的類別一致性。雖然個別超像素被分為錯誤的類別，例如上方中央位置的農田(這主要是由于該處較大的光譜變化增加了分割難度)。M1的分割結果與M0相近，但M1未能精確分割S1子影像下方中央位置牧草場和收割后草場的交界區域(圖5(b))，這主要是由于2種地物的光譜較為接近，再加上該處超像素與實際地物邊界未能完全符合，導致交互勢函數難以正確引導分割結果，這也在一定程度上顯示了本文β自動估計方法的優越性。通過觀察圖5(c)和(d)，可以發現M2和M3將很多河流周邊的區域錯分為河流，導致精度較低。這主要是由于M2和M3的交互勢函數未能充分利用空間背景信息以及不同類別之間的β參數差異，導致分割效果欠佳。另外，M3的分割結果中很多玉米被錯分為小麥，除了其交互勢函數模型的欠缺外，其超像素合并的求解策略容易加大類間區分的難度，也是影響其精度的一個重要原因。

S2子影像4種算法的分割結果如圖6所示。

(a) M0 (b) M1 (c) M2 (d) M3

圖6S2子影像4種MRF算法的分割結果

Fig.6S2subsetsegmentationresultsof4differentMRFalgorithms

對比圖6(a)和(b)可以發現，M0與M1的分割結果在總體上較為接近，其主要差異在于玉米和葵花這2個類別的區分上。由圖2(b)可見，S2中的玉米和葵花具有較為相似的光譜特征，且這2個類別在空間上交錯分布，增加了這2個類別區分的難度。通過仔細對比M0與M1的結果，可以觀察到M0對于玉米和葵花的區分效果更佳，這一現象在S2子影像右下方的農田較為明顯，M1將其中的很多葵花錯分為玉米，而M0在該區的分割更為成功。盡管整體上M2的分割精度與M0十分接近，但圖6(c)顯示M2將一些玉米錯分為葵花。這在一定程度上說明，M2雖然可以較為有效地區分易混淆的類別，但在S2的分割上其效果依舊次于本文算法。對于M3，圖6(d)顯示其分割結果中很多休耕地和農田被錯分為城鎮，分割效果較差。

為了進一步定量評價S1與S2子影像4種算法的分割精度，各算法的OA和Kappa如表2所示。

表2 S1與S2子影像4種算法的分割精度定量評價Tab.2 Quantitative evaluation of segmentationaccuracy by 4 algorithms for S1 and S2 subsets

從表2中可以看出，S1子影像M0的分割精度最高，M3的精度最差，M0和M1的分割精度較為接近，這與圖5的觀察是一致的； S2子影像雖然整體精度都較高，但同樣是M0的分割精度最高，M3的精度最差，也定量地說明本文算法具有最佳的分割精度。以上實驗結果和分析都較為有力地證明了本文算法的優勢。對于農田區域的HRI，本文的超像素MRF模型可以更為準確地利用邊界強度等空間背景信息，從而引導MRF求解過程得到更為精確的分割結果。

3 結論

本文發展了一種超像素MRF的影像分割算法，以提升高空間分辨率遙感數據在農田地區的信息提取精度。本文算法將高階鄰域模型引入到MRF的交互勢函數中，并采用了逐類別的β參數估計方法。通過實驗驗證，得出結論如下：

1)將高階鄰域模型引入到超像素MRF影像分割算法中，可以抑制超像素與實際地物邊界不完全符合造成的誤差對MRF求解的負面影響，使MRF更充分地利用空間背景信息，從而提高分割精度。

2)提出β自動估計方法是基于范數距離計算的，相比于傳統的馬氏距離，該方法能夠更為精確地估計不同類別信息對β參數調節的貢獻，有利于提高分割精度。

通過對高階鄰域參數hn的分析表明，不同影像的最佳hn可能存在差異。因此在未來的工作中，發展最佳hn的自動估計方法是十分必要的。另外，本文僅考慮了農田地區的光學遙感影像； 如何將本文算法拓展到SAR影像或非農田區域的遙感影像，也是值得探索的方向。

[1] L?w F,Conrad C,Michel U.Decision fusion and non-parametric classifiers for land use mapping using multi-temporal RapidEye data[J].ISPRS Journal of Photogrammetry and Remote Sensing,2015,108:191-204.

[2] Kim H O,Yeom J M.Effect of red-edge and texture features for object-based paddy rice crop classification using RapidEye multi-spectral satellite image data[J].International Journal of Remote Sensing,2014,35(19):7046-7068.

[3] 楊閆君,占玉林,田慶久,等.基于GF-1/WFV NDVI時間序列數據的作物分類[J].農業工程學報,2015,31(24):155-161.

Yang Y J,Zhan Y L,Tian Q J,et al.Crop classification based on GF-1/WFV NDVI time series[J].Transactions of the Chinese Society of Agricultural Engineering,2015,31(24):155-161.

[4] Liu M W,Ozdogan M,Zhu X J.Crop type classification by simultaneous use of satellite images of different resolutions[J].IEEE Transactions on Geoscience and Remote Sensing,2014,52(6):3637-3649.

[5] Besag J.On the statistical analysis of dirty pictures[J].Journal of the Royal Statistical Society,Series B Methodological,1986,48(3):259-302.

[6] Kumar S,Hebert M.Discriminative random fields[J].International Journal of Computer Vision,2006,68(2):179-201.

[7] 孫 巍,郭 敏.基于SLIC與條件隨機場的圖像分割算法[J].計算機應用研究,2015,32(12):3817-3820,3824.

Sun W,Guo M.Image segmentation based on SLIC and conditional random field[J].Application Research of Computers,2015,32(12):3817-3820,3824.

[8] 劉 磊,石志國,宿浩茹,等.基于高階馬爾可夫隨機場的圖像分割[J].計算機研究與發展,2013,50(9):1933-1942.

Liu L,Shi Z G,Su H R,et al.Image segmentation based on higher order Markov random field[J].Journal of Computer Research and Development,2013,50(9):1933-1942.

[9] Yu Q Y,Clausi D A.IRGS:Image segmentation using edge penalties and region growing[J].IEEE Transactions on Pattern Analysis and Machine Intelligence,2008,30(12):2126-2139.

[10] Qin A K,Clausi D A.Multivariate image segmentation using semantic region growing with adaptive edge penalty[J].IEEE Transactions on Image Processing,2010,19(8):2157-2170.

[11] 蘇騰飛,李洪玉.一種結合GBM的MRF遙感圖像分割算法[J].內蒙古農業大學學報(自然科學版),2015,36(1):143-149.

Su T F,Li H Y.GBM-combined MRF method for remote sensing image segmentation[J].Journal of Inner Mongolia Agricultural University(Natural Science Edition),2015,36(1):143-149.

[12] Achanta R,Shaji A,Smith K,et al.SLIC superpixels compared to state-of-the-art superpixel methods[J].IEEE Transactions on Pattern Analysis and Machine Intelligence,2012,34(11):2274-2282.

[13] Comaniciu D,Meer P.Mean shift:A robust approach toward feature space analysis[J].IEEE Transactions on Pattern Analysis and Machine Intelligence,2002,24(5):603-619.

[14] Levinshtein A,Stere A,Kutulakos K N,et al.TurboPixels:Fast superpixels using geometric flows[J].IEEE Transactions on Pattern Analysis and Machine Intelligence,2009,31(12):2290-2297.

[15] Su T F,Li H Y,Zhang S W,et al.Image segmentation using mean shift for extracting croplands from high-resolution remote sensing imagery[J].Remote Sensing Letters,2015,6(12):952-961.