蘇騰飛， 張圣微， 李洪玉

(內(nèi)蒙古農(nóng)業(yè)大學(xué)水利與土木建筑工程學(xué)院,呼和浩特 010018)

基于可變尺度Mean-Shift的農(nóng)田高分遙感影像分割算法

蘇騰飛， 張圣微， 李洪玉

(內(nèi)蒙古農(nóng)業(yè)大學(xué)水利與土木建筑工程學(xué)院,呼和浩特 010018)

為了提升農(nóng)田高分遙感影像(high spatial resolution remote sensing image,HRI)的信息提取效果，提出了一種新的農(nóng)田HRI分割算法。傳統(tǒng)的Mean-Shift(MS)HRI分割算法僅利用全局或單一的尺度參數(shù)； 而常規(guī)可變尺度MS算法在尺度參數(shù)估算中也只考慮光譜信息。這些都導(dǎo)致其分割結(jié)果難以完整地展現(xiàn)不同尺度的農(nóng)田區(qū)域。針對(duì)該問(wèn)題，在MS算法的基礎(chǔ)上進(jìn)行了改進(jìn)： 第一，提出了一種局部可變尺度參數(shù)的估計(jì)方法； 第二，提出了利用局部可變尺度進(jìn)行MS濾波的模型函數(shù)。該改進(jìn)算法主要包含3步： ①為了全面考慮不同波段的響應(yīng)變化，在MS濾波核函數(shù)中采用了對(duì)角化的尺度參數(shù)矩陣，并將其與采樣點(diǎn)密度估計(jì)模型相結(jié)合，導(dǎo)出了一種可變尺度MS濾波的迭代函數(shù)； ②為了提高算法的自動(dòng)化程度，利用局部光譜變化與邊界強(qiáng)度信息，提出了一種新的局部尺度參數(shù)估算方法； ③將MS濾波結(jié)果輸入到基于分形網(wǎng)絡(luò)演化方法(fractal net evolution approach,FNEA)的空間聚類算法中，得到最終的分割結(jié)果。利用RapidEye與OrbView3的2景HRI進(jìn)行了算法驗(yàn)證。實(shí)驗(yàn)結(jié)果表明，所提出的改進(jìn)算法能夠優(yōu)化農(nóng)田HRI分割的精度。

可變尺度； Mean-Shift； 農(nóng)田分割； 高分遙感影像

0 引言

近年來(lái)，一種新興的遙感影像解譯方式——面向地理對(duì)象的影像分析(geographical object based image analysis,GEOBIA)[1]在農(nóng)田遙感分類研究中逐漸增多[2-4]，這主要?dú)w因于2個(gè)方面： ①農(nóng)田遙感監(jiān)測(cè)需要高空間分辨率遙感影像(high spatial resolution remote sensing image,HRI)，以準(zhǔn)確提取作物信息。很多地區(qū)(如中國(guó))的農(nóng)田面積較小[5]，導(dǎo)致中低空間分辨率的遙感影像易產(chǎn)生混合像素光譜，HRI可有效地解決這一問(wèn)題，很多研究都表明GEOBIA適用于HRI的解譯[1-2]； ②在GEOBIA中，除農(nóng)田的光譜特征外，形狀和背景特征也能被有機(jī)地融合到作物種類判別算法中，從而顯著提高了作物識(shí)別精度。例如一些地區(qū)使用了圈形噴灌機(jī)，使農(nóng)田呈規(guī)則的圓形或扇形，這可被GEOBIA用來(lái)區(qū)分農(nóng)田和其他地物?？梢?jiàn)，利用GEOBIA和HRI開(kāi)展農(nóng)田分類研究具有重要意義。

GEOBIA首先就是影像分割，其目的是獲取地物的原始圖斑； 然后通過(guò)特征選擇和影像分類，賦予原始圖斑以真實(shí)地物的語(yǔ)義信息。不準(zhǔn)確的分割結(jié)果會(huì)導(dǎo)致后續(xù)分類的嚴(yán)重誤判，因此影像分割是GEOBIA的關(guān)鍵步驟。雖然目前已有大量的影像分割模型與軟件，但如何提高農(nóng)田HRI的分割精度依然是具有挑戰(zhàn)性的課題。這是因?yàn)椋?①HRI中的農(nóng)田區(qū)域內(nèi)部經(jīng)常出現(xiàn)較大的光譜變化； ②在同一地區(qū)，不同農(nóng)田的形狀、大小存在很大的差異； ③一些相鄰農(nóng)田的邊界十分模糊。以上因素都給傳統(tǒng)的影像分割帶來(lái)了一定困難。

本文在Mean-Shift(MS)遙感影像分割算法的基礎(chǔ)上進(jìn)行了改進(jìn)，使之更好地適用于農(nóng)田HRI的分割。之所以選用MS算法，主要有2個(gè)原因： ①M(fèi)S在GEOBIA中應(yīng)用廣泛[6]，具有抗噪性好[7]、能處理較大尺度影像[8]等優(yōu)點(diǎn)； ②基于MS的HRI分割算法可分為2步，第一步是邊界保持平滑濾波[9]； 第二步是空間聚類。其中第一步可以在不使邊界模糊的前提下有效地平滑農(nóng)田內(nèi)部的光譜值，從而提高農(nóng)田HRI的分割精度。眾多學(xué)者對(duì)基于MS的遙感影像分割算法進(jìn)行了較為深入的研究。Wang等[10]將一種監(jiān)督學(xué)習(xí)的邊界提取算法與MS濾波相結(jié)合，以優(yōu)化較大尺度HRI的分割效果。Ponti[11]利用MS算法處理從遙感影像中提取的植被指數(shù)，實(shí)現(xiàn)了MS在農(nóng)田遙感影像分割中的應(yīng)用。Lang等[7]提出了一種可變尺度參數(shù)的MS濾波算法，并將其用于抑制極化SAR影像的斑點(diǎn)噪聲。Banerjee等[12]利用k最近鄰(k-nearest neighbor,kNN)算法自動(dòng)估計(jì)MS的尺度參數(shù)，以提高M(jìn)S分割算法的自動(dòng)化程度。Huang等[13]發(fā)展了一套自適應(yīng)尺度參數(shù)選擇的MS分割算法，用以分割高光譜影像。Comaniciu等[14]提出了一種數(shù)據(jù)驅(qū)動(dòng)的可變尺度MS分割算法。Ming等[6]則利用空間統(tǒng)計(jì)方法來(lái)估計(jì)可變尺度MS分割算法中的尺度參數(shù)?？v觀上述研究可以發(fā)現(xiàn)，大多MS分割算法僅采用了全局或單一的尺度參數(shù)； 一些可變尺度的MS算法在其尺度參數(shù)估計(jì)中也僅考慮了光譜變化信息，而忽略了邊界等空間結(jié)構(gòu)信息。在實(shí)際農(nóng)田的HRI中，不同農(nóng)田的面積和田間光譜變化具有較大差異，這都會(huì)給常規(guī)MS分割策略帶來(lái)困難。

為了提高M(jìn)S分割算法在農(nóng)田HRI分割中的性能，本文提出了一種可變尺度的MS分割算法，其過(guò)程主要分為3步： ①將對(duì)角化的尺度參數(shù)矩陣與采樣點(diǎn)密度估計(jì)模型相結(jié)合，導(dǎo)出可變尺度MS濾波的迭代函數(shù)； ②利用局部光譜變化信息與矢量梯度邊界強(qiáng)度，提出一種新的局部尺度參數(shù)估算方法； ③將MS濾波結(jié)果輸入到基于分形網(wǎng)絡(luò)演化方法(fractal net evolution approach,FNEA)[15]的區(qū)域生長(zhǎng)算法中，得到農(nóng)田HRI的分割結(jié)果。并利用中國(guó)和美國(guó)的2景HRI進(jìn)行算法驗(yàn)證，以證明該算法可顯著優(yōu)化農(nóng)田HRI的分割精度。

1 算法原理

1.1 可變尺度MS濾波

MS最早由Fukunaga等[16]提出，用于多變量特征空間的密度中心估計(jì)。Comaniciu等[9]將MS運(yùn)用到影像特征空間分析中，發(fā)展了基于MS的影像分割算法。該算法主要包含 MS濾波和基于空間聚類的影像分割2個(gè)步驟。其中第一步的目的是獲取每個(gè)像素點(diǎn)的“mode”[9]，即像素點(diǎn)所屬類塊的密度中心。這實(shí)際上是對(duì)原始影像進(jìn)行了平滑濾波，但又與普通的低通平滑濾波不同——由于采用了基于核函數(shù)的密度梯度估計(jì)，MS濾波可以有效地保持區(qū)域邊界信息，因此該方法也被稱為“邊界保持平滑濾波”，該步驟對(duì)提高分割精度非常關(guān)鍵。

MS濾波的一般步驟是： ①對(duì)于影像中任意一個(gè)像素x初始化，選代次數(shù)j=1并且y1=x； ②當(dāng)jt，則第j+1次迭代濾波結(jié)果為

(1)

式中：t為收斂條件閾值(1個(gè)很小的正數(shù))；N()為濾波窗口內(nèi)的像素集合；KH()為核函數(shù)(多采用高斯函數(shù)形式)；H為KH()的尺度參數(shù)(亦稱其為“帶寬參數(shù)”)； ③得到濾波結(jié)果zi=yi,j+1。

通常，尺度參數(shù)有2種表達(dá)形式，其中最簡(jiǎn)單的形式[9]為

H1=h2I,

(2)

h=hs=hr，

(3)

式中hs和hr分別為空間和光譜尺度參數(shù)。這一形式的優(yōu)點(diǎn)是只需給出1個(gè)尺度參數(shù)h，即可完成濾波計(jì)算。在實(shí)際計(jì)算中，hs與hr可以不相等，此時(shí)的核函數(shù)為

(4)

式中：k()為核函數(shù)K()的輪廓函數(shù)[9]；C1為歸一化常數(shù)。

對(duì)于多光譜遙感影像，固定的hr很難適用于大多數(shù)實(shí)際情況，因?yàn)椴煌ǘ螌?duì)同一地物的響應(yīng)可能會(huì)有較大差別。所以第二種尺度參數(shù)的形式[9]更接近實(shí)際情況，即

(5)

式中：p為波段數(shù)； diag()表示對(duì)角矩陣。為了簡(jiǎn)便，將式(5)寫成向量形式，即

(6)

然后在式(6)的基礎(chǔ)上，進(jìn)而提出改變了的核函數(shù)形式，即

(7)

式中：C2為歸一化常數(shù)；hr，b為波段b的尺度參數(shù)。將式(7)帶入式(1)，可得基于H2的MS濾波迭代函數(shù)，即

(8)

雖然式(8)考慮了不同波段尺度參數(shù)的差異，但仍難以滿足實(shí)際情況的要求。在HRI成像過(guò)程中，由于氣候條件和人為干擾等因素，同一波段對(duì)不同地點(diǎn)的同種地物會(huì)存在光譜響應(yīng)上的差異。因此，需要在H2的基礎(chǔ)上考慮局部尺度參數(shù)。1.2節(jié)將詳細(xì)敘述提出的局部尺度參數(shù)估計(jì)方法，此處只給出局部尺度參數(shù)(即可變尺度的濾波迭代函數(shù))。

可變尺度的MS濾波算法有2種模型，分別是氣球和采樣點(diǎn)估計(jì)函數(shù)，在文獻(xiàn)[14]中均有詳細(xì)介紹，此處不再贅述。因?yàn)槲墨I(xiàn)[16]指出第一種模式的效果較差，本文采用第二種模式。在文獻(xiàn)[14]中式(17)的基礎(chǔ)上，將本文修改的式(6)代入其中。由于本文是按照式(7)的假設(shè)來(lái)進(jìn)行推導(dǎo)的，即認(rèn)為空間信息與光譜信息是相互獨(dú)立的，因此在密度中心估計(jì)中像素的位置與各波段光譜信息是連乘的形式。仿照式(8)，將其中隨窗口中心變化的尺度參數(shù)替換為隨任意像素位置變化的尺度參數(shù)形式，可得到本文最終的MS迭代函數(shù)，即

(9)

1.2 局部尺度參數(shù)估計(jì)

在農(nóng)田HRI中，要將各個(gè)農(nóng)田完整地分割出來(lái)是較為困難的。一方面，由于灌溉條件和耕作習(xí)慣等因素，即便在同一片農(nóng)田中作物之間的生長(zhǎng)速率也存在差異，這在HRI中反映為農(nóng)田內(nèi)部的光譜變化，會(huì)導(dǎo)致同一片農(nóng)田被錯(cuò)誤地分為若干區(qū)域； 另一方面，農(nóng)田之間經(jīng)常被較窄的小道或灌渠分隔，在HRI中以模糊邊界的形式顯示出來(lái)，這可能會(huì)導(dǎo)致一些光譜特征相近的農(nóng)田被分割為同一片區(qū)域。綜上可見(jiàn)，局部光譜變化和邊界強(qiáng)度信息對(duì)確定MS濾波參數(shù)是非常重要的。

前人在設(shè)計(jì)局部尺度參數(shù)的估計(jì)模型中，基本上只考慮了局部光譜變化信息[6,14]。從定性上來(lái)說(shuō)，較大的尺度參數(shù)易引起過(guò)平滑(即邊界信息的丟失)； 反之，較小的尺度參數(shù)難以對(duì)光譜變化較大的農(nóng)田區(qū)域進(jìn)行較好的平滑。鑒于此，在提出的尺度估計(jì)模型中，對(duì)邊界強(qiáng)度較高的區(qū)域，其尺度參數(shù)的估計(jì)結(jié)果應(yīng)適當(dāng)偏??； 而對(duì)于變化較大的農(nóng)田內(nèi)部，其尺度估計(jì)結(jié)果應(yīng)適當(dāng)偏大。綜合考慮以上因素，局部尺度參數(shù)應(yīng)與邊界強(qiáng)度成反比，而與農(nóng)田內(nèi)部的光譜變化成正比。為了避免因局部光譜變化較大而導(dǎo)致其被錯(cuò)誤地當(dāng)做邊界信息，在本文模型中，采取了單調(diào)遞減指數(shù)函數(shù)形式將邊界強(qiáng)度信息與局部光譜信息結(jié)合到一起，使其更適合農(nóng)田HRI的分割。其計(jì)算公式為

hr,edge(x,y,b)=exp[-En(x,y)]hr(x,y,b) ,

(10)

式中：En(x,y)為像素(x,y)的歸一化邊界強(qiáng)度值；hr(x,y,b)為像素(x,y)在波段b的尺度參數(shù)，即

(11)

(12)

其中：v為像素的光譜強(qiáng)度值；W(x,y,b)為以像素(x,y)為中心的窗口所包含的波段b子影像；I(b)為波段b的整景影像；c為影像最大灰度差異常數(shù)； 10與ln 2為經(jīng)驗(yàn)常數(shù)，均通過(guò)大量實(shí)驗(yàn)得出。由式(11)可以看出，對(duì)于變化較大的局部區(qū)域，其尺度參數(shù)較小，從而防止過(guò)度平滑。

HRI通常具有多個(gè)波段，因此采用了一種基于矢量梯度的邊界強(qiáng)度信息提取方法[17]。該方法需要首先計(jì)算出2個(gè)正交方向的梯度影像，本文采用高斯差分濾波來(lái)實(shí)現(xiàn)這一步驟，得到的梯度矢量為

(13)

式中?gh,1和?gv,1分別為波段1水平和豎直方向的梯度值。邊界強(qiáng)度值即為2×2矩陣VT·V的最大特征值λ。在實(shí)際計(jì)算中需要采用歸一化的邊界強(qiáng)度，即

(14)

從式(10)可以看出，對(duì)于影像中灰度變化較大或邊界強(qiáng)度較高的部位，所估計(jì)的局部尺度參數(shù)較小，這可以避免過(guò)大尺度可能導(dǎo)致的過(guò)度平滑。另外，對(duì)于內(nèi)部灰度變化較大的區(qū)域，由于其邊界強(qiáng)度通常較小(農(nóng)田內(nèi)部的灰度變化通常較為平緩)，因此其估計(jì)的尺度參數(shù)不會(huì)特別小，這有利于區(qū)域內(nèi)部的平滑。對(duì)于邊界較為模糊的部位，式(10)第一項(xiàng)依然可以使其局部尺度較低，從而降低過(guò)度平滑的概率。

1.3 基于FNEA的遙感影像分割

本文的主要目的是利用MS濾波算法提高農(nóng)田HRI的分割效果，因此在空間濾波階段采用了一種主流的遙感影像分割算法——FNEA來(lái)實(shí)現(xiàn)整體尺度的HRI分割。FNEA是一種基于區(qū)域生長(zhǎng)的空間聚類算法，在初始階段，各個(gè)像素被當(dāng)作單獨(dú)的斑塊； 然后通過(guò)迭代合并相似的斑塊，使均一的區(qū)域被分割出來(lái)。FNEA利用均一性變化量來(lái)衡量2個(gè)圖斑適于合并的程度，而均一性變化量的計(jì)算則考慮了形狀和光譜信息[14]。在算法執(zhí)行之前，F(xiàn)NEA需要設(shè)置形狀異質(zhì)性權(quán)重fs、緊湊異質(zhì)性權(quán)重fc與斑塊大小閾值fsize這3個(gè)參數(shù)； 前2個(gè)參數(shù)分別影響分割過(guò)程中光譜和形狀信息在均一性變化量計(jì)算中的比重，fsize則決定分割結(jié)果中圖斑的平均大小。一般而言，fs與fc分別被默認(rèn)為0.1和0.5。由于本文主要考察MS濾波結(jié)果對(duì)農(nóng)田HRI的分割效果，因此在實(shí)驗(yàn)中FNEA的fs與fc均采用了默認(rèn)設(shè)置。圖1示出了本文算法的流程。

圖1 本文算法流程

FNEA的實(shí)現(xiàn)參考了文獻(xiàn)[18]，利用基于紅黑搜索樹(shù)的全局最優(yōu)合并來(lái)實(shí)現(xiàn)快速、準(zhǔn)確的HRI分割。本文算法基于C++語(yǔ)言完成，實(shí)驗(yàn)平臺(tái)為Windows7操作系統(tǒng)，4Gb RAM，CPU為Intel Core I5-4200M。

2 算法實(shí)驗(yàn)

2.1 實(shí)驗(yàn)數(shù)據(jù)

本文實(shí)驗(yàn)共采用了2景農(nóng)田HRI，分別為RapidEye和OrbView3影像(以下簡(jiǎn)稱為S1和S2)，其參數(shù)見(jiàn)表1。

表1 實(shí)驗(yàn)數(shù)據(jù)信息

上述2種數(shù)據(jù)均包含藍(lán)光、綠光、紅光和近紅外波段，而RapidEye還提供紅邊波段(690～740 nm)。為便于對(duì)比分割效果，對(duì)S1和S2影像均采用了綠(Band2)、紅(Band3)、近紅外(Band4)波段組合。S1和S2影像Band4(R)Band3(G)Band2(B)假彩色合成結(jié)果如圖2所示。

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

圖2實(shí)驗(yàn)用Band4(R)Band3(G)Band2(B)假彩色合成影像

Fig.2FalsecolorimagescomposedofBand4(R)Band3(G)Band2(B)adoptedforexperiment

局部影像的位置分別見(jiàn)圖2中的白色和藍(lán)色方框，下文中分別將S1中的白色和藍(lán)色方框內(nèi)的局部影像簡(jiǎn)稱為S1a和S1b，將S2中的白色和藍(lán)色方框內(nèi)的局部影像簡(jiǎn)稱為S2a和S2b。在圖2中，S1影像顯示了中國(guó)山東地區(qū)的農(nóng)田，其田塊的平均面積較小，很多由田間小道和灌渠形成的邊界比較模糊，并且大多相鄰的田塊具有相似的光譜特征。相反，S2影像顯示了美國(guó)加利福尼亞地區(qū)的農(nóng)田，其田塊平均面積較大，形狀大多為規(guī)則的矩形，而且許多田塊之間的光譜變化較大。以上因素都會(huì)為準(zhǔn)確分割影像帶來(lái)一定的困難，因此，對(duì)S1和S2影像的分割實(shí)驗(yàn)可較為全面地驗(yàn)證本文算法對(duì)農(nóng)田HRI分割的優(yōu)化效果。

2.2 對(duì)比實(shí)驗(yàn)與結(jié)果分析

2.2.1 實(shí)驗(yàn)設(shè)置

為了驗(yàn)證本文算法的性能，采用了與其他3種算法進(jìn)行對(duì)比的實(shí)驗(yàn)，包括基于可變尺度的MS分割算法[14](簡(jiǎn)稱“VMS”)以及2種較新的MS分割算法。VMS采用的迭代濾波計(jì)算式與本文中的式(7)相同，其尺度參數(shù)的估計(jì)僅考慮了光譜信息。另2種較新的算法中，一種算法在文獻(xiàn)[12]中被提出，首先采用kNN進(jìn)行全局最優(yōu)光譜尺度參數(shù)hr的估計(jì)，然后利用估計(jì)的hr進(jìn)行MS濾波； 另一種算法是邊界加權(quán)(edge weighted, EW)的MS濾波[10]。

VMS,kNN和EW算法與本文算法均需要給出空間尺度參數(shù)hs，該參數(shù)也被當(dāng)作濾波窗口的尺度。經(jīng)過(guò)多次實(shí)驗(yàn)，將4種算法的hs均設(shè)置為5，一些相關(guān)研究[7-8,10]也都采用了相同的設(shè)置。除hs外，EW還需設(shè)置hr。在本文實(shí)驗(yàn)中發(fā)現(xiàn)，EW采用kNN估計(jì)的hr可以得到滿意的效果。另外，EW的輸入還包括邊界強(qiáng)度影像，由1.2節(jié)中的矢量梯度方法提供。

為便于對(duì)比，對(duì)4種MS濾波算法均采用FNEA完成最終的HRI分割。而對(duì)于S1和S2影像，其fsize需要分別通過(guò)測(cè)試和調(diào)節(jié)給出最佳設(shè)置。經(jīng)多次實(shí)驗(yàn)，發(fā)現(xiàn)S1和S2影像的最優(yōu)fsize分別為2 000和35 000。S2影像的參數(shù)更大的原因是因?yàn)镾2影像中農(nóng)田的平均尺度遠(yuǎn)大于S1影像。

為了定量地給出4種算法對(duì)農(nóng)田HRI分割的優(yōu)化效果，采用了一種定量的遙感影像分割評(píng)價(jià)方法[18-20]。該方法首先給出分割結(jié)果的精確率(P)和召回率(R)，分別反映過(guò)分割和亞分割錯(cuò)誤，并且數(shù)值范圍均為[0,1]。當(dāng)P和R均接近1時(shí)，說(shuō)明分割結(jié)果是精確的；P接近1而R接近0時(shí)，表明存在過(guò)分割錯(cuò)誤； 反之，P接近0而R接近1時(shí)，說(shuō)明存在亞分割錯(cuò)誤。再利用該方法根據(jù)P和R計(jì)算出F，以全面地反映分割精度。同樣，F(xiàn)的范圍也是[0,1]，其值越接近1說(shuō)明越精確。另外，為了更清晰地展現(xiàn)4種MS濾波算法對(duì)分割效果的提升，還與原始影像作為FNEA輸入時(shí)的分割結(jié)果進(jìn)行了對(duì)比。

由于S1和S2影像的尺度較大，觀察全景影像的分割結(jié)果不便于進(jìn)行細(xì)致、可靠的分割效果評(píng)價(jià)。因此，在S1和S2影像中分別選取了2個(gè)具有代表性的局部影像S1a,S1b和S2a,S2b進(jìn)行更為詳細(xì)的分割結(jié)果評(píng)價(jià)。

2.2.2 實(shí)驗(yàn)結(jié)果

圖3和圖4為本文算法的中間結(jié)果，其中圖3為矢量梯度方法提取的邊界強(qiáng)度。

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

圖3矢量梯度邊界強(qiáng)度

Fig.3Vectorgradientedgestrength

由圖3可以看出，大部分明顯的邊界都被提取出來(lái)了。雖然一些較弱的邊界沒(méi)有被明顯地顯示出來(lái)，但其邊界強(qiáng)度值依然高于區(qū)域內(nèi)部。圖4是各個(gè)波段的局部尺度影像。

(a) S1Band4影像 (b) S1Band3影像(c) S1Band2影像

(d) S2Band4影像 (e) S2Band3影像(f) S2Band2影像

圖4局部尺度參數(shù)估計(jì)結(jié)果

Fig.4Resultsoflocalscaleparameterestimation

從圖4可以看出，邊界位置基本具有較小的尺度參數(shù)； 而對(duì)于農(nóng)田內(nèi)部，尤其是灰度變化較小的區(qū)域，其尺度參數(shù)的數(shù)值較大。

對(duì)于濾波結(jié)果，kNN為S1和S2影像估計(jì)的最佳全局hr分別為6.27和28.65，前者小于后者主要是因?yàn)镾1影像中田塊的平均面積較小。

圖5所示本文算法對(duì)S1和S2影像的最終分割結(jié)果。

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

圖5本文算法整體分割結(jié)果

Fig.5Wholesegmentationresultsproducedbyproposedmethod

仔細(xì)觀察圖5可以看出，S1影像中主要的區(qū)域被較好地分割出來(lái)了； 而S2影像中較大的區(qū)域存在過(guò)分割現(xiàn)象。通過(guò)與S2影像原圖對(duì)比可以發(fā)現(xiàn)，S2影像中大部分被過(guò)分割的農(nóng)田，其內(nèi)部的灰度變化較大，且其中存在樹(shù)木和水塘等小尺度的地物，這些都是導(dǎo)致過(guò)分割的主要原因。因傳統(tǒng)MS，kNN與EW算法的最終分割結(jié)果與本文算法較為相似，為節(jié)省篇幅，未顯示其整景分割結(jié)果。

表2給出了對(duì)4種算法分割結(jié)果的定量評(píng)價(jià)。

表2 4種算法分割實(shí)驗(yàn)結(jié)果定量評(píng)價(jià)

從表2可見(jiàn)，對(duì)于2景實(shí)驗(yàn)影像，本文算法的F均最高，這說(shuō)明本文算法對(duì)農(nóng)田HRI分割優(yōu)化的效果較為顯著。值得一提的是，4種算法的P,R和F相差并不懸殊，這在圖6和圖7局部影像的分割效果對(duì)比中也可以得到印證。

仔細(xì)觀察圖6第一行可發(fā)現(xiàn)，本文算法的濾波結(jié)果對(duì)較大農(nóng)田內(nèi)部的平滑效果最佳(例如圖6(a)右下部顏色較深的大塊田地)。對(duì)于S1a和S1b影像，本文算法都得到了較好的分割結(jié)果。例如，S1a影像右上部的幾塊農(nóng)田之間有模糊的小道，在圖6第二行中，只有本文算法較好地區(qū)分了這些田地(圖6(e))； 其他算法在該區(qū)域普遍存在亞分割錯(cuò)誤(即把相鄰的農(nóng)田合并到了一起)。對(duì)于S1b影像，經(jīng)目視對(duì)比可見(jiàn)，4種算法的分割結(jié)果差異較小。kNN在S1b影像的分割中得到了較好的效果，其中個(gè)別田塊(例如圖6(o)下部的矩形農(nóng)田)被較為完整地分割出來(lái)。但在S1b影像右下部的細(xì)長(zhǎng)農(nóng)田分割中，本文算法將其中1個(gè)農(nóng)田錯(cuò)誤地分為2個(gè)小區(qū)域(圖6(m))； 除此以外，本文算法在S1b影像其他區(qū)域的分割中都取得了較為不錯(cuò)的效果。

(a) S1a影像本文濾波結(jié)果(b) S1a影像VMS濾波結(jié)果(c) S1a影像kNN濾波結(jié)果(d) S1a影像EW濾波結(jié)果

(e) (a)的分割結(jié)果 (f) (b)的分割結(jié)果 (g) (c)的分割結(jié)果 (h) (d)的分割結(jié)果

(i) S1b影像本文濾波結(jié)果 (j) S1b影像VMS濾波結(jié)果(k) S1b影像kNN濾波結(jié)果(l) S1b影像EW濾波結(jié)果

(m) (i)的分割結(jié)果 (n) (j)的分割結(jié)果 (o) (k)的分割結(jié)果 (p) (l)的分割結(jié)果

圖6S1局部影像的實(shí)驗(yàn)結(jié)果

Fig.6TestresultsofS1subsets

與圖6不同的是，圖7中的田塊具有更大的面積，但田塊內(nèi)部的光譜變化也更為顯著。一些農(nóng)田(特別是S2b影像中的田塊)具有較大的光譜變化，以至于4種算法都產(chǎn)生了過(guò)分割錯(cuò)誤。在S2a影像的分割結(jié)果中，EW的過(guò)分割錯(cuò)誤最為明顯； 在S2b影像中，4種算法均有明顯的過(guò)分割錯(cuò)誤。

(a) S2a影像本文濾波結(jié)果 (b) S2a影像VMS濾波結(jié)果(c) S2a影像kNN濾波結(jié)果(d) S2a影像EW濾波結(jié)果

(e) (a)的分割結(jié)果 (f) (b)的分割結(jié)果 (g) (c)的分割結(jié)果 (h) (d)的分割結(jié)果

(i) S2b影像本文濾波結(jié)果 (j) S2b影像VMS濾波結(jié)果(k) S2b影像kNN濾波結(jié)果(l) S2b影像EW濾波結(jié)果

(m) (i)的分割結(jié)果 (n) (j)的分割結(jié)果 (o) (k)的分割結(jié)果 (p) (l)的分割結(jié)果

圖7S2局部影像的實(shí)驗(yàn)結(jié)果

Fig.7TestresultsofS2subsets

對(duì)圖7進(jìn)一步觀察可以看出，本文算法、VMS與kNN的分割結(jié)果相差較小。然而表2中的定量評(píng)價(jià)結(jié)果顯示，本文算法優(yōu)于其他3種算法，這主要是因?yàn)閳D7只顯示了一部分面積較小農(nóng)田的分割結(jié)果，而表2中的定量評(píng)價(jià)結(jié)果是基于整景分割結(jié)果計(jì)算出來(lái)的。一些面積較大的田塊被傳統(tǒng)MS算法和kNN誤分割為若干區(qū)域，則是導(dǎo)致其F較低的主要原因。S2影像的分割實(shí)驗(yàn)可以更可靠地說(shuō)明，本文算法對(duì)農(nóng)田HRI分割具有較好的優(yōu)化效果。

3 結(jié)論

本文發(fā)展了一種基于可變尺度Mean-Shift(MS)的農(nóng)田高分遙感影像(HRI)分割算法，其主要貢獻(xiàn)在于將局部灰度變化與邊界強(qiáng)度信息加入到了局部尺度參數(shù)估計(jì)中，從而較為顯著地提高了MS濾波的邊界保持能力與區(qū)域平滑效果。利用2景較大尺度的農(nóng)田HRI進(jìn)行了實(shí)驗(yàn)驗(yàn)證，結(jié)果表明本文算法具有如下優(yōu)勢(shì)：

1)較高程度的自動(dòng)化。由于本文提出的局部尺度參數(shù)估計(jì)方法不需要預(yù)先設(shè)置任何參數(shù)，因而大大方便了用戶的操作。

2)優(yōu)越的邊界保持平滑濾波。實(shí)驗(yàn)結(jié)果說(shuō)明本文算法對(duì)農(nóng)田HRI分割的優(yōu)化效果要好于其他MS濾波算法。

本文算法的運(yùn)算速度相對(duì)較慢，對(duì)RapidEye影像進(jìn)行分割用時(shí)近2 min。在未來(lái)可以考慮通過(guò)優(yōu)化編程代碼或簡(jiǎn)化收斂條件來(lái)提升計(jì)算效率。在下一步工作中，也可以考慮本文算法在城市HRI或高分SAR影像分割中的應(yīng)用。

[1] Blaschke T,Hay G J,Kelly M,et al.Geographic object-based image analysis-towards a new paradigm[J].ISPRS Journal of Photogrammetry and Remote Sensing,2014,87:180-191.

[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] Vieira A,Formaggio A R,Rennó C D,et al.Object based image analysis and data mining applied to a remotely sensed Landsat time-series to map sugarcane over large areas[J].Remote Sensing of Environment,2012,123:553-562.

[5] 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.

[6] Ming D P,Li J,Wang J Y,et al.Scale parameter selection by spatial statistics for GEOBIA:Using mean-shift based multi-scale segmentation as an example[J].ISPRS Journal of Photogrammetry and Remote Sensing,2015,106:28-41.

[7] Lang F K,Yang J,Li D R,et al.Mean-shift-based speckle filtering of polarimetric SAR data[J].IEEE Transactions on Geoscience and Remote Sensing,2014,52(7):4440-4454.

[8] Michel J,Youssefi D,Grizonnet M.Stable mean-shift algorithm and its application to the segmentation of arbitrarily large remote sensing images[J].IEEE Transactions on Geoscience and Remote Sensing,2015,53(2):952-964.

[9] 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.

[10]Wang L G,Liu G Y,Dai Q L.Optimization of segmentation algorithms through mean-shift filtering preprocessing[J].IEEE Geoscience and Remote Sensing Letters,2014,11(3):622-626.

[11]Ponti M P.Segmentation of low-cost remote sensing images combining vegetation indices and mean shift[J].IEEE Geoscience and Remote Sensing Letters,2013,10(1):67-70.

[12]Banerjee B,Varma S,Buddhiraju K M,et al.Unsupervised multi-spectral satellite image segmentation combining modified mean-shift and a new minimum spanning tree based clustering technique[J].IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing,2014,7(3):888-894.

[13]Huang X,Zhang L P.An adaptive mean-shift analysis approach for object extraction and classification from urban hyperspectral imagery[J].IEEE Transactions on Geoscience and Remote Sensing,2008,46(12):4173-4185.

[14]Comaniciu D,Ramesh V,Meer P.The variable bandwidth mean shift and data-driven scale selection[C]//Proceedings of the 8th International Conference of Computer Vision.Vancouver:IEEE,2001:438-445.

[15]Baatz M,Sch?pe A.Multiresolution segmentation:An optimization approach for high quality multi-scale image segmentation[M]//Strobl J,Baschke T,Griesebner G.Angewandte Geographische Informations Verarbeitung XII.Karlsruhe:Wichmann Verlag,2000:12-23.

[16]Fukunaga K,Hostetler L.The estimation of the gradient of a density function,with applications in pattern recognition[J].IEEE Transactions on Information Theory,1975,21(1):32-40.

[17]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.

[18]蘇騰飛,李洪玉.一種兩階段區(qū)域生長(zhǎng)的遙感圖像分割算法[J].遙感技術(shù)與應(yīng)用,2015,30(3):476-485. Su T F,Li H Y.A two stage region growing method for remote sensing image segmentation[J].Remote Sensing Technology and Application,2015,30(3):476-485.

[19]Jiang H B,Su Y Y,Jiao Q S,et al.Typical geologic disaster surveying in Wenchuan 8.0 earthquake zone using high resolution ground LiDAR and UAV remote sensing[C]//Proceedings of SPIE 9262,Lidar Remote Sensing for Environmental Monitoring XIV.Beijing,China:SPIE,2014:926219.

[20]Yi L N,Zhang G F,Wu Z C.A scale-synthesis method for high spatial resolution remote sensing image segmentation[J].IEEE Transactions on Geoscience and Remote Sensing,2012,50(10):4062-4070.

(責(zé)任編輯：陳理)

VariablescaleMean-Shiftbasedmethodforcroplandsegmentationfromhighspatialresolutionremotesensingimages

SU Tengfei， ZHANG Shengwei， LI Hongyu

(CollegeofWaterConservancyandCivilEngineering，InnerMongolianAgriculturalUniversity，Hohhot010018，China)

In order to improve the effect of information extraction from high spatial resolution remote sensing images (HRI) of cropland, the authors put forward a new HRI segmentation algorithm. Due to the fact that the traditional Mean-Shift (MS) segmentation method only uses a global and single scale, and that some variable bandwidth MS only considers spectral information in their scale estimation process, and croplands with various sizes could be hardly extracted in one segmentation result, the authors improved a MS based approach to tackle this problem. The main consideration lies in two aspects: ① A local variable scale parameter estimation method is proposed; ② The model function for local variable scale is established for MS filtering. The proposed approach mainly consists of 3 parts: ① With the objective of comprehensively considering the response variation of different bands, the diagonal scale parameter matrix is adopted in the kernel function of MS filtering, and it is combined with sample point estimation model to derive the iterative function for variable scale MS filtering; ② For the purpose of increasing automation of the proposed method, local spectral variation and edge strength information are utilized to design a new local scale parameter estimation method; ③ For obtaining the final segmentation, the filtering result is used as input for the fractal net evolution approach (FNEA) which is a spatial clustering method. Two scenes of HRI acquired by RapidEye and OrbView3 were employed for experiment, and the results show that the proposed method can optimize the accuracy of cropland HRI segmentation.

variable scale； Mean-Shift； cropland segmentation； high spatial resolution remote sensing image

10.6046/gtzyyg.2017.03.06

蘇騰飛,張圣微,李洪玉.基于可變尺度Mean-Shift的農(nóng)田高分遙感影像分割算法[J].國(guó)土資源遙感,2017,29(3)：41-50.(Su T F,Zhang S W,Li H Y.Variable scale Mean-Shift based method for cropland segmentation from high spatial resolution remote sensing images[J].Remote Sensing for Land and Resources,2017,29(3)：41-50.)

2015-08-17；

2015-10-21

國(guó)家自然科學(xué)基金項(xiàng)目“科爾沁沙地典型生態(tài)系統(tǒng)水熱通量傳輸機(jī)理及其與植被耦合關(guān)系試驗(yàn)和模擬研究”(編號(hào)： 51569017)、“內(nèi)蒙古典型草原水文過(guò)程及其擾動(dòng)與觸發(fā)草地退化的水文臨界條件實(shí)驗(yàn)與模擬研究”(編號(hào)： 51269014)、內(nèi)蒙古自然科學(xué)基金項(xiàng)目“半干旱區(qū)沙地典型生態(tài)系統(tǒng)水熱通量傳輸機(jī)理研究”(編號(hào)： 2015MS0514)和中國(guó)博士后科學(xué)基金面上資助項(xiàng)目“西部地區(qū)博士后人才資助計(jì)劃”(編號(hào)： 2015M572630XB)共同資助。

蘇騰飛(1987-),男,碩士,實(shí)驗(yàn)師,主要從事面向?qū)ο蟮倪b感圖像分析算法方面的研究。Email： stf1987@126.com。

張圣微(1979-),男,博士,教授,主要從事定量遙感、生態(tài)水文及氣候變化等方面的研究。Email： zsw_imau@163.com。

TP 751.1

： A

： 1001-070X(2017)03-0041-10