基于分水嶺變換的高分辨率遙感影像分割

繆曉冬 陳曾敏

(南通市測繪院有限公司, 江蘇 南通 226000)

0 引言

遙感技術(shù)能在不接觸目標(biāo)的情況下獲取其幾何、物理、位置等信息,如何更好地從遙感影像中獲取信息是近幾十年來遙感研究的主要挑戰(zhàn)[1]。傳統(tǒng)的基于像元的遙感影像處理方法都是建立在遙感影像光譜信息極其豐富、地物間光譜差異較為明顯的基礎(chǔ)上的。而隨著傳感器技術(shù)的發(fā)展,遙感影像的空間分辨率變得越來越高,高分辨影像的分析存在著一些問題,即地物內(nèi)部像元的差異可能非常大,基于像元的影像分析方法無法充分利用高分辨率遙感影像上的信息,所以對于高分辨率遙感影像,傳統(tǒng)的基于像元的分析方法不再適用,而基于對象的高分辨率遙感影像分析方法更加適合[2-3]。

基于對象的影像分析方法通常包括兩個步驟：一是影像分割,二是特征提取和分類。其中最關(guān)鍵、最基本的步驟是影像分割。影像分割得到勻質(zhì)的、不重疊的圖像對象,然后根據(jù)各個對象的屬性(光譜、紋理、形狀等)進(jìn)行分類對象的特征提取,分類的準(zhǔn)確性主要取決于影像分割的質(zhì)量。分水嶺變換是一種基于區(qū)域的方法[4],分水嶺變換具有局部一致性特性,分割結(jié)果保留了更加精細(xì)的細(xì)節(jié),分水嶺變換最大的優(yōu)勢是可以提取出完整的邊緣,不需要邊緣連接的過程。但是,根據(jù)分水嶺變換的定義,影像上只要有一個局部極小值點(diǎn),就會產(chǎn)生一個分割區(qū)域,整個影像上局部極小值點(diǎn)非常多,再加上噪聲影響,不管是將其應(yīng)用在原始影像圖上,還是影像梯度圖上,都會產(chǎn)生嚴(yán)重的過分割現(xiàn)象。

Giuseppe Masi[5],Giuseppe Scarpa[6],Raffaele Gaetano[7],Giovanni Poggi[8]和Luisa Verdoliva[9]5人逐步提出了一種改進(jìn)的分水嶺算法,從使用形態(tài)學(xué)標(biāo)記的邊界標(biāo)志填充算法(edge mark fill,EMF),到進(jìn)一步引入光譜標(biāo)記的EMF+算法,最后到能夠處理同時包含全色數(shù)據(jù)和多光譜數(shù)據(jù)的高分辨率遙感影像的高分遙感邊界標(biāo)志填充算法(multispectral remote-sensing edge mark fill,MR-EMF),三種算法主要基于標(biāo)記控制來改進(jìn)原始分水嶺算法的過分割現(xiàn)象。栗敏光等人[10]采用形態(tài)學(xué)混合開閉重構(gòu)濾波對原始影像進(jìn)行平滑處理,以此來抑制原始分水嶺算法的過分割現(xiàn)象;楊軍等人[11]提出了一種基于巴氏系數(shù)和灰度共生矩陣的區(qū)域合并方法來改進(jìn)分水嶺算法的分割結(jié)果;閆鵬飛等[12]針對傳統(tǒng)分水嶺分割算法易受噪聲干擾且分割尺度參數(shù)難以自適應(yīng)選擇的問題,提出了一種尺度自適應(yīng)的高分辨遙感影像分水嶺分割方法;閆利等[13]采用鄰域自適應(yīng)閾值、區(qū)域合并方法對種子點(diǎn)選取過程進(jìn)行優(yōu)化,提出基于多尺度標(biāo)記和邊緣約束的多尺度分割算法。但是在許多學(xué)者的研究中,算法不夠完整,缺乏預(yù)處理及后處理步驟,本文針對分水嶺變換的過分割問題,提出一種改進(jìn)的分水嶺變換的高分辨率遙感影像分割算法,來補(bǔ)充基礎(chǔ)算法缺失的預(yù)處理及后處理步驟,從而達(dá)到抑制原始分水嶺變換過分割現(xiàn)象、提高分割精度的目的。

1 算法流程

本文主要以基于邊緣距離變換的形態(tài)學(xué)標(biāo)記及光譜標(biāo)記控制的分水嶺算法為基礎(chǔ),結(jié)合形態(tài)學(xué)混合開閉濾波操作、基于巴氏系數(shù)和灰度共生矩陣的區(qū)域合并方法,來補(bǔ)充基礎(chǔ)算法缺失的預(yù)處理及后處理步驟,研究了算法各步驟對分割結(jié)果的影響,并根據(jù)基于對象匹配的定量分析方法,對本文算法與易康軟件(eCognition)中的多尺度分割算法進(jìn)行定量對比實(shí)驗(yàn)分析,算法流程圖如圖1所示。

1.1 形態(tài)學(xué)混合開閉濾波

分水嶺變化直接作用于原影像或者梯度圖像上時,由于原始高分辨率遙感影像像元灰度值往往很復(fù)雜并且還時常伴有噪聲,會造成分割的輪廓偏移、產(chǎn)生“偽區(qū)域”以及導(dǎo)致過分割。為了消除噪聲并削弱同類區(qū)域內(nèi)部細(xì)節(jié),要對原始影像進(jìn)行濾波平滑處理。而與傳統(tǒng)的線性濾波方式不同,基于數(shù)學(xué)形態(tài)學(xué)的非線性濾波方法在消除噪聲、削弱細(xì)節(jié)的同時能保留對象的幾何形狀特征,更適合與分水嶺分割相結(jié)合。用于混合開閉重構(gòu)操作的結(jié)構(gòu)元素的形狀和尺寸對影像去噪和平滑也有影響,這里使用圓形結(jié)構(gòu)元素,為了與之后的形態(tài)學(xué)標(biāo)記提取使用的圓形結(jié)構(gòu)元素一致,其半徑需小于影像中最小分割對象的最小外接圓的半徑,防止損失最小待分割地物對象。

1.2 基于邊緣的標(biāo)記控制的分水嶺分割

分水嶺變換最早來自地貌學(xué)基本概念,將灰度圖像視作地形圖,每個點(diǎn)的灰度值作為它的高程,其基本思想為：在地形圖上的每個區(qū)域極小值打一個洞,讓整個區(qū)域從上往下勻速降入水中,直至其被完全淹沒,這期間,當(dāng)不同匯水盆地上升的水匯集時,在匯集地建立水壩來阻止匯集,最終建立的水壩即為分水嶺,本文調(diào)用矩陣實(shí)驗(yàn)室(matrix laboratory，MATLAB)中watershed函數(shù)實(shí)現(xiàn)分水嶺變換。

1.2.1Canny邊緣檢測算法

1.2.1.1消除噪音

傳統(tǒng)Canny邊緣檢測使用高斯濾波器對圖像進(jìn)行平滑去噪處理,(2k+1)×(2k+1)尺寸的高斯核Hij構(gòu)造如下：

(3)

式中,σ為標(biāo)準(zhǔn)差;k為不小于0的常數(shù)。高斯核的尺寸會影響邊緣檢測的效果,尺寸越大,對噪聲的靈敏度越低,邊緣定位誤差也會增加;同時,標(biāo)準(zhǔn)差σ也會影響邊緣檢測效果。

1.2.1.2計(jì)算梯度

使用索伯算子(Sobel)得到水平方向和垂直方向的一階導(dǎo)數(shù)值(Gx,Gy),梯度和梯度方向計(jì)算如下：

式中,GX和Gy分別為垂直方向和水平方向的索伯濾波器;G為梯度;θ為梯度方向。將連續(xù)的梯度方向劃分為四個離散的方向：水平、垂直以及兩個對角線方向,為下一步非最大化抑制做準(zhǔn)備。

1.2.1.3非最大化抑制

此步驟起到細(xì)化邊緣的作用,對每個像素進(jìn)行掃描,檢查像素是否是在其梯度方向上附近的局部最大值,若是,則保留為邊緣;若不是,則賦梯度值為0去除。

1.2.1.4雙閾值滯后檢測

設(shè)置一高一低兩個梯度閾值,對邊緣像素進(jìn)行篩選：若梯度值大于高閾值,則為邊緣;若梯度值小于小閾值,則舍去;若梯度值介于兩者之間,判斷其是否與大于高閾值地邊緣像素相連,若是則保留為邊緣,若不是則舍去。

1.2.2形態(tài)學(xué)標(biāo)記方法

1.2.2.1邊緣檢測

對原始圖像進(jìn)行Canny邊緣檢測,得到邊緣圖。

1.2.2.2距離變換

將邊緣圖(二值圖,1代表邊緣,0代表非邊緣)的像素值替換該像素到最近邊緣的距離,得到距離變換圖,并取負(fù)得到最終進(jìn)行分水嶺分割的數(shù)字高程模型(digital elevation model,DEM),此時邊緣的像素值為極大值0,距離邊緣越遠(yuǎn),像素值越小。

1.2.2.3尋找極小值標(biāo)記種子

尋找DEM的局部極小值,設(shè)置為起始標(biāo)記種子。

1.2.2.4標(biāo)記膨脹

對種子點(diǎn)進(jìn)行形態(tài)學(xué)膨脹操作,膨脹的結(jié)構(gòu)元素為圓形,其半徑為種子點(diǎn)到最近邊緣的距離再減去一個縮進(jìn)距離Δ(默認(rèn)值為3個像素),當(dāng)兩個種子點(diǎn)膨脹后相交時,將其歸并為一個標(biāo)記。

1.2.3光譜信息提取與形態(tài)學(xué)標(biāo)記合并算法

形態(tài)學(xué)標(biāo)記的生成完全基于原始影像邊緣和形態(tài)學(xué)準(zhǔn)則,而忽略了影像的光譜信息,因此利用光譜信息提取光譜標(biāo)記,并最后將其與形態(tài)學(xué)標(biāo)記合并會取得更好的分割效果。同形態(tài)學(xué)標(biāo)記的提取一樣,為了保留邊緣,將邊緣像素及其周圍像素剔除,在剩余區(qū)域提取光譜標(biāo)記,并根據(jù)光譜標(biāo)記連通區(qū)域的光譜方差對較大方差區(qū)域進(jìn)行分割。

1.2.3.1生成光譜區(qū)域

對邊緣進(jìn)行膨脹操作,再取反,使得生成光譜標(biāo)記的區(qū)域不包含邊緣及其相鄰區(qū)域。

1.2.3.2提取初始光譜標(biāo)記

通過分水嶺擬合標(biāo)記侵蝕(watershed marked erosion,WME)使光譜標(biāo)記的每個連通區(qū)域與基于邊緣的分水嶺分割保持一致。

1.2.3.3 樹型馬爾科夫隨機(jī)場(tree shape Markov random field,TS-MRF)

統(tǒng)計(jì)各個光譜標(biāo)記的連通區(qū)域的光譜均方差,對均方差最大的1%(默認(rèn)值)個連通區(qū)域進(jìn)行TS-MRF分割為兩類。

1.2.3.4分水嶺擬合標(biāo)記侵蝕

再對1.2.3.3產(chǎn)生的光譜標(biāo)記使用WME,得到最后的光譜標(biāo)記。1.2.3.1和1.2.3.2都是為了使得最后分割圖保留邊緣信息,1.2.3.3以均方差作為判定依據(jù),對均方差較大的光譜標(biāo)記連通區(qū)域使用TS-MRF分割為兩類,以此更好地將不同地物分割出來。

1.3 基于巴氏系數(shù)和灰度共生矩陣的區(qū)域合并

巴氏系數(shù)是對兩個離散統(tǒng)計(jì)樣本重疊度的近似計(jì)算,在灰度影像上即灰度分布直方圖的相似度計(jì)算,區(qū)域1與區(qū)域2之間的巴氏系數(shù)B(R1,R2)定義如下：

(6)

式中，L為影像的灰度級；p(i)為灰度為i的像素?cái)?shù)量占總像素的比例,巴氏系數(shù)越接近1,兩區(qū)域越相似,同類地物的概率越大。

灰度共生矩陣是對灰度值之間空間關(guān)系的統(tǒng)計(jì),用來描述影像的紋理特征,其矩陣維度一般為影像的灰度級,每個灰度共生矩陣元素記為：Pd(i,j),其中d表示在方向,即統(tǒng)計(jì)在該方向上的灰度值空間關(guān)系,(i,j)為一對灰度值,Pd(i,j)具體的值定義為在方向d上灰度值對(i,j)出現(xiàn)的次數(shù),d包括方向和步長,根據(jù)灰度共生矩陣,可以定義出14個紋理特征標(biāo)量,相關(guān)研究分析發(fā)現(xiàn)只有角二階矩(angular second moment, ASM)、熵(entropy, ENT)、對比度(contrast, CON)和相關(guān)性(correlation, COR)4個紋理特征與人類視覺感知特性有明確的對應(yīng)關(guān)系[15],因此本文采用這4個特征量來描述紋理信息,作為相鄰區(qū)域合并的一條準(zhǔn)則。角二階矩反映圖像灰度分布均勻程度和紋理粗細(xì)程度,熵反映紋理的復(fù)雜程度或非均勻程度,對比度反映了圖像的清晰度和紋理溝紋的深淺程度,相關(guān)性反映的是紋理的方向性。計(jì)算這些紋理特征標(biāo)量時需將灰度共生矩陣歸一化(每個元素除以矩陣所有元素之和)。

2 實(shí)驗(yàn)分析

對上述完整的算法流程進(jìn)行實(shí)驗(yàn)分析,定性、定量地給出其性能評價,并將其與基于距離變換的分水嶺分割算法以及商業(yè)軟件eCognition Developer中使用的多尺度分割技術(shù)進(jìn)行比較。

2.1 數(shù)據(jù)介紹及預(yù)處理

使用公開數(shù)據(jù)集Zurich_dataset_v1.0進(jìn)行實(shí)驗(yàn)分析。該數(shù)據(jù)集獲取自2002年8月份瑞士蘇黎世地區(qū)的20幅QuickBird影像,包括紅、綠、藍(lán)及近紅外四個波段,并全色銳化至0.6 m的分辨率。同時,還包括20張相應(yīng)的地物類別標(biāo)注圖,共注釋了8類地物：道路、樹木、裸土、鐵路、建筑物、草地、水體和泳池,但并不是每個像素都被標(biāo)注。本文選取第十二幅(782×622)進(jìn)行分析,記為zh12,如圖2所示。

圖2 zh12實(shí)驗(yàn)原始影像

圖像的預(yù)處理主要是針對公開數(shù)據(jù)集Zurich_dataset_v1.0提供的地物類別標(biāo)注圖,地物分類圖中同一類別的像元共有一個標(biāo)注,而影像分割與分類不同,需要將同一類別的不同對象賦予不同的標(biāo)注。因此對分類圖進(jìn)行處理,使得每一類的不同對象分開,賦予不同的標(biāo)注。

2.2 參數(shù)設(shè)置及定量精度評價方法

使用基于區(qū)域的評價標(biāo)準(zhǔn),其主要概念為：將算法所得分割圖中的每個對象與真實(shí)分割圖中的每個對象進(jìn)行匹配,并根據(jù)區(qū)域判定標(biāo)準(zhǔn)將真實(shí)分割圖的每個對象分為四類。

2.2.1正確分割(CS)

若真實(shí)分割對象(G)與算法分割對象(S)的重疊區(qū)域大于各自的75%(α),則將真實(shí)分割對象記為正確分割：

(11)

2.2.2過分割(OS)

若真實(shí)分割對象被多個算法分割對象所分割,且滿足下式,則將真實(shí)分割對象即為過分割。

(12)

2.2.3欠分割(US)

若算法分割對象被多個真實(shí)分割對象所分割,且滿足下式,則將真實(shí)分割對象記為欠分割。

(13)

2.2.4漏分割(MEs)

除去上述三類真是分割對象后,剩下的即為漏分割。

2.3 結(jié)果分析

下面選取本文算法、基于距離變換的分水嶺算法、eCognition-30和eCognition-120得到的分割結(jié)果進(jìn)行展示,如圖3所示。

(a)本文算法

統(tǒng)計(jì)分割對象尺寸與像元數(shù)的分布情況,如圖4所示,橫坐標(biāo)x表示所有尺寸處于(2x-1,2x]分割對象,縱坐標(biāo)y表示這些分割對象的像元數(shù)量之和。可以看出,本文算法得到的分割結(jié)果尺寸分布比較均勻,大尺寸、小尺寸對象都有,最后還有個明顯的最大值,這是因?yàn)樵跋褡髠?cè)的水體被很好地分割出來;而基于距離變換的分水嶺算法得到的分割結(jié)果的尺寸都偏小,大值很少;而eCognition軟件的分割結(jié)果的尺寸分布與尺度參數(shù)相關(guān)性很大,都是分布在一個峰值尺寸左右,尺度參數(shù)越大,峰值尺寸越大。

圖4 分割對象尺寸與該尺寸像元數(shù)關(guān)系

根據(jù)對象匹配原則,計(jì)算各個算法的基于區(qū)域的精度評價得分,如表1所示,各項(xiàng)得分是根據(jù)被標(biāo)記為CS、OS、US和MEs的真實(shí)分割區(qū)域面積和與所有參與對象匹配的分割對象面積和之比,并且不是所有參與對象匹配的分割對象都有對應(yīng)的真實(shí)分割對象,第六列為各算法得到的分割對象數(shù)量。

從表1可以看出,本文算法的CS(正確分割)得分是最高的,過分割得分處于中間水平,欠、漏分割很少,基于距離變換的分水嶺分割算法、eCognition-20、eCognition-30和eCognition-50過分割現(xiàn)象很嚴(yán)重,eCognition軟件的后三項(xiàng)分割得分以及分割數(shù)和尺度參數(shù)相關(guān)性很高,尺度參數(shù)越大,OS得分越小、US得分越大、MEs得分越大,分割數(shù)量越少,即過分割現(xiàn)象越少、欠分割、漏分割現(xiàn)象越多,而CS得分呈現(xiàn)先上升再下降的情況,最優(yōu)尺度參數(shù)在80附近。

表1 各分割算法的對象匹配分?jǐn)?shù)

選取的四個算法所得分割結(jié)果的對象匹配效果從圖5可以看出,本文提出的算法很好地把原影像左側(cè)的水體分割出來,其余三種算法均對水體得到過分割結(jié)果。并且,尺度參數(shù)為120的eCognition算法的漏分割現(xiàn)象很嚴(yán)重。

圖5 各算法分割結(jié)果的分割匹配圖

3 結(jié)束語

本文以基于邊緣距離變換的形態(tài)學(xué)標(biāo)記及光譜標(biāo)記控制的分水嶺分割為基礎(chǔ),加入了形態(tài)學(xué)混合開閉濾波預(yù)處理操作以及基于巴氏系數(shù)及灰度共生矩陣的區(qū)域合并分割后操作,在保留了邊緣信息的情況下,充分利用了影像的光譜信息以及紋理信息,大大改善了原始分水嶺分割的過分割現(xiàn)象,提高了建筑物這種具有明顯邊界的對象的分割精度,以及植被區(qū)域的分割效果也顯著增強(qiáng)。

根據(jù)基于對象匹配的定量評價方法,在公開數(shù)據(jù)集及自制標(biāo)注數(shù)據(jù)上將本文算法與eCognition軟件中的多尺度分割算法進(jìn)行對比,展現(xiàn)出了本文算法的優(yōu)勢：①總體分割視覺效果更好,最終分割對象數(shù)量適中;②正確分割的定量得分較高;③初始參數(shù)設(shè)置較少(eCognition軟件中的多尺度分割算法需要耗時的試錯過程來尋找最佳分割尺度參數(shù))。