薛 峭，趙書河

(南京大學(xué)地理信息科學(xué)系，南京 210093)

基于最小核值相似區(qū)算法的高分辨率遙感圖像分割方法

薛 峭，趙書河

(南京大學(xué)地理信息科學(xué)系，南京 210093)

采用最小核值相似區(qū)算法(Smallest Univalue Segment Assimilating Nucleus，SUSAN)計算QuickBird圖像的梯度，并采用標(biāo)記控制的分水嶺變換(Watershed Transform，WT)算法分割圖像，取得了較好的結(jié)果。SUSAN方法能有效地檢測圖像梯度，對噪聲不敏感;梯度值范圍明確，不因圖像而改變，為后續(xù)處理相關(guān)參數(shù)的選擇提供了便利;亮度閾值容易確定，模板半徑可選，具有很大的靈活性;適合于采用WT的遙感圖像的分割。采用基于SUSAN梯度和NDVI的標(biāo)記圖像，利用形態(tài)學(xué)灰度圖像重建方法修改梯度圖像，能夠有效地抑制梯度圖像中大量的局部灰度極小值，提高WT圖像分割的精度。

最小核值相似區(qū)算法(SUSAN);分水嶺變換(WT);QuickBird圖像;圖像分割

0 引言

高分辨率遙感圖像中地物景觀的結(jié)構(gòu)、形狀、紋理和細(xì)節(jié)等信息非常突出，基于像元分類的傳統(tǒng)算法難以從高空間分辨率遙感數(shù)據(jù)中提取所需的信息［1］，因此面向?qū)ο蟮母叻直媛蔬b感圖像處理方法得到不斷發(fā)展［2］。圖像分割是高分辨率遙感圖像信息提取和面向?qū)ο蠓诸惖闹匾h(huán)節(jié)［3］。

實現(xiàn)圖像分割的方法包括基于邊界的方法和基于區(qū)域的方法，前者檢測局部變化，后者則檢測像元和區(qū)域的相似性［4］。分水嶺變換(Watershed Transform，WT)屬于第二種方法。基于數(shù)學(xué)形態(tài)學(xué)原理的WT方法被廣泛地應(yīng)用于灰度圖像分割中。最初，WT方法被用來分割氣泡圖像和金相顯微鏡圖像，從圖像中提取出氣泡和顯微鏡下的目標(biāo)［5］。將圖像看作地表，灰度值為其高程值，在每個局部最低點位置打一個洞并且讓水從洞中涌出，從低到高淹沒整個地形;當(dāng)處在不同匯水盆地中的水要聚合到一起時，修建大壩阻止其聚合，大壩就對應(yīng)圖像的分割線［6］。分割線是閉合的，與圖像灰度的等值線一致，閉合的分水嶺線包圍的區(qū)域就是圖像分割中的目標(biāo)。

直接對圖像進(jìn)行WT分割將得到2個不好的結(jié)果:①很多分水嶺線沒有位于目標(biāo)的邊緣。這樣，即使匯水盆地大小合理，也不能保證分水嶺線包圍的是完整的目標(biāo)(理想的情況是目標(biāo)的邊緣高程值即灰度值高于非邊緣的灰度值，分水嶺線位于高處即形成目標(biāo)邊緣);②由于過分割，區(qū)域破碎，目標(biāo)內(nèi)部存在大量高程極小值，WT方法會形成眾多面積微小的匯水盆地。因此，在WT分割中常使用梯度圖像來保證目標(biāo)與其邊緣的高程差異，并且采取措施抑制梯度圖像目標(biāo)內(nèi)部大量高程極小值，以消除微小的匯水盆地。WT圖像分割的性能依賴于計算圖像梯度的方法和抑制梯度圖像極小值的方法［3］。典型的梯度檢測算子都用核算子來計算正交方向上的坡度值的強度，然后計算出總的梯度強度值［4］，但上述計算容易受到噪聲的影響。Canny梯度檢測屬于線性濾波梯度檢測算法，檢測步驟包括高斯平滑去噪、垂直和水平方向上的高斯算子邊緣檢測、細(xì)化邊緣的非極大值抑制和雙閾值處理(高閾值得到初步的梯度，低閾值保證初步梯度的連續(xù)性［7］)，其結(jié)果容易受到高斯濾波器標(biāo)準(zhǔn)差和高、低2個閾值的影響。1997年Stephen等［8］提出用最小核值相似區(qū)算法(Smallest Univalue Segment Assimilating Nucleus，SUSAN)檢測圖像的邊緣，其實質(zhì)是根據(jù)像元周圍局部區(qū)域內(nèi)與該像元亮度值相近的像元個數(shù)來判斷該點是否為邊緣。不用對圖像求導(dǎo)，可直接對圖像灰度進(jìn)行運算，效率高，定位準(zhǔn)確，具有積分特性，對噪聲不敏感［9－10］。

本文采用SUSAN方法計算QuickBird圖像梯度，并采用基于梯度和NDVI提取的標(biāo)記圖像對梯度進(jìn)行重建，抑制大量的梯度圖像的灰度極小值;最后進(jìn)行圖像的WT分割，取得較好的分割效果。

1 實驗數(shù)據(jù)

所用實驗數(shù)據(jù)為南京地區(qū)的QuickBird圖像，包括藍(lán)(450～520 nm)、綠(520～660 nm)、紅(630～690 nm)和近紅外(760～900 nm)4個多光譜波段，空間分辨率為2.44 m;全色波段(450～900 nm)，空間分辨率為0.61 m。本文截取的數(shù)據(jù)是512像元×512像元的全色圖像及與其對應(yīng)的多光譜圖像。典型地物包括河流、建筑物和植被等(圖1)。

圖1 實驗區(qū)QuickBird全色圖像Fig.1 QuickBird panchromatic image of the test site

2 SUSAN梯度計算

2.1 SUSAN 算法原理

SUSAN算法的原理如圖2所示。

圖2 SUSAN原理圖Fig.2 Principle of SUSAN

圖2中給出幾個有代表性的位置，暗色和周邊代表灰度圖像中的2類地物。

用如圖3所示的一個圓形模板遍歷圖像，模板半徑是3，圓形模版的中心位置稱作核(nucleus)，模板除去中心核外包含的區(qū)域大小m是36個像元(圖3中的暗色像元)。

圖3 SUSAN模板(r=3，m=36)Fig.3 SUSAN matrix(r=3，m=36)

若模板內(nèi)其他像元的灰度值與模板中心像元的灰度值之差小于一定的亮度差閾值t，就認(rèn)為二者相似。由相似的像元組成的區(qū)域稱為任意像元的核值相似區(qū)(USAN)［10］。USAN包含了重要的圖像結(jié)構(gòu)信息:在非邊緣處，USAN值較大(最大可為模板所含像元個數(shù));越靠近邊緣，USAN值越小。在圖2中的1處，模板位于角點，USAN值為1/4模板大小;當(dāng)模板的中心位于邊緣時，USAN值為1/2模板大小(如圖2中的3處);在圖2中的5處，USAN值則小于1/4模板大小。圖像中每一個像元對應(yīng)于一個USAN值，因此通過USAN值就能檢測出圖像的邊緣。用上述模板遍歷圖像，在任意像元位置，計算該像元的USAN值為

式中，u為像元的 USAN值，u(r0，r)=t為亮度差閾值，表示模板中心核像元亮度值與其他像元亮度值的差別，差別小則認(rèn)為二者相似。Stephen等［8］討論了t值的選取與圖像灰度范圍之間的關(guān)系，對于8 bit灰度圖像建議t取25。實驗表明，對于QuickBird圖像而言，t取值在25～40之間較合適，本文取t=30。

實際采用 u(r0，r)計算像元間的相似度，函數(shù)值在0和1之間。遍歷圖像完成計算后，得到一個新的浮點型灰度圖像，稱作該圖像的USAN圖，它反映每個像元的同質(zhì)區(qū)的大小(像元灰度最大值約為36，最小值接近0)，值大的像元為非邊緣，值小的像元為邊緣。

為了得到適合WT分割的梯度圖像，對USAN圖像進(jìn)一步做2步處理:①進(jìn)行灰度值求反，用該圖像的最大值減去每個像元的值作為該像元的值，使得目標(biāo)邊緣處高程值大，以適合WT的圖像分割;②將值小于此時圖像最大值的1/4倍的像元賦值為0，以初步抑制弱小梯度。得到的梯度圖像被稱作SUSAN梯度圖像。

2.2 SUSAN梯度圖像獲取

取研究區(qū)全色圖像，將4個多光譜波段的空間分辨率重采樣成與全色波段一樣大小。分別計算上述5幅圖像中每個像元的SUSAN梯度值，取最大值作為該位置上的最終梯度值。圖4為模板半徑r=5時得到的梯度圖像。大多數(shù)地物的內(nèi)部與邊緣有較好的可分性，但河流旁的樹木區(qū)域內(nèi)部與邊緣的差異不夠明顯，需要進(jìn)一步的處理。

圖4 SUSAN梯度圖(r=5)Fig.4 SUSAN gradients(r=5)

SUSAN梯度的重要特點是它的取值范圍只與模板半徑r有關(guān)，為后續(xù)處理相關(guān)參數(shù)的選擇提供便利;同時r可靈活選擇，以適應(yīng)不同的需求。下文將用到這個特點。

3 WT圖像分割

常用基于標(biāo)記控制的WT方法解決過分割問題［11］，前提是要知道被分割目標(biāo)的大致位置［12］。在梯度圖像中，目標(biāo)的位置在梯度值較小處，梯度值較大處則是目標(biāo)的邊緣。然而梯度較小的區(qū)域中仍存在較多的梯度極小值，會造成分割的破碎。反映圖像中目標(biāo)大概位置的另一種圖像叫做標(biāo)記圖像。典型的標(biāo)記圖像是目標(biāo)處為1、非目標(biāo)處為0且與原圖像大小相同的二值圖像。根據(jù)梯度圖像和標(biāo)記圖像，利用形態(tài)學(xué)灰度重建原理進(jìn)行梯度修改，可得到新的梯度圖像［13］。此時，對應(yīng)于標(biāo)記位置的梯度變得均一，因此分水嶺線不會破壞標(biāo)記范圍內(nèi)的目標(biāo)。標(biāo)記的作用是確定目標(biāo)的大致位置(以便于對梯度圖像進(jìn)行重建)，與目標(biāo)越接近的標(biāo)記越好。

3.1 標(biāo)記提取

圖像的梯度特征和光譜特征常用于提取標(biāo)記圖像。擴(kuò)展最小變換(Extended－minima Transform，EM)是常見的基于梯度的標(biāo)記提取方法［14］，它把目標(biāo)的大致位置處標(biāo)記為1，其余地方標(biāo)記為0。梯度圖像G經(jīng)過高度閾值為h的擴(kuò)展最小變換運算為

式中，E為輸出的與G大小相同的二值圖像;h為大于零的高程值。

若局部極小區(qū)域的最大高程差大于h，則在輸出圖像對應(yīng)的局部極小區(qū)位置上輸出1，否則為0。在得到的標(biāo)記圖像中，值為1處則對應(yīng)于目標(biāo)的大致位置。

同時采用基于光譜特征的標(biāo)記圖像提取方法。植被內(nèi)部與其邊緣梯度值的差異較小，但植被的光譜特征明顯，可采用NDVI來確定植被的范圍，利用NDVI提取植被的標(biāo)記圖像。將滿足NDVI＞0.2的區(qū)域作為植被的標(biāo)記，與擴(kuò)展最小變換提取的標(biāo)記做“或”運算，滿足其中之一者即可作為最終的標(biāo)記(最終的標(biāo)記圖像如圖5)。

圖5 標(biāo)記圖像Fig.5 Marker image

利用灰度圖像重建方法對梯度圖像修改之后，梯度圖像中的標(biāo)記處變成了均勻的極小區(qū)域(如圖6所示)?；叶戎亟刹捎枚喾N計算方法［13］，本文采用了形態(tài)學(xué)灰度圖像重建方法。

圖6 修改后的梯度圖像Fig.6 Modified gradient image

采用SUSAN梯度，再獲取基于梯度和NDVI的標(biāo)記圖像，用標(biāo)記控制的WT分割方法對QuickBird圖像進(jìn)行分割(分割結(jié)果如圖7)。整個過程在MATLAB R2007b軟件環(huán)境下實現(xiàn)。

圖7 分割結(jié)果Fig.7 Segmentation result

3.2 精度分析與參數(shù)選擇

由圖7可以看出，大多地物都能被較好地分割開來，但個別建筑物被分成了幾個區(qū)域，這與它們的復(fù)雜結(jié)構(gòu)有關(guān)。EM中對高度閾值h的選擇與圖像的梯度大小有關(guān)，梯度范圍與模板半徑r有關(guān)(模板半徑為 r，梯度范圍約是0～πr2，最大值 m約為πr2)。當(dāng)h滿足 h/m 在 0.1～0.4之間時，能取得較好的分割結(jié)果。本文中 r=5，h=20，h/m≈0.25，區(qū)域個數(shù)=259。隨著r的增加，梯度的寬度也會增加。取r=5，當(dāng)采用基于NDVI的標(biāo)記提取方法時，植被區(qū)域被選擇成標(biāo)記，成為灰度極小區(qū)域;區(qū)域邊緣的外側(cè)仍然存在著梯度，即“山谷”的周圍存在包圍該區(qū)域的“山脊”，成為區(qū)域的分水嶺線。此時r取值稍大更合適，而在其他圖像的分割中可選擇較小半徑的模板以提高計算速度。

對分割精度的評價采用區(qū)域個數(shù)與區(qū)域標(biāo)準(zhǔn)差的均值來衡量——區(qū)域個數(shù)反映分割結(jié)果是否破碎，區(qū)域標(biāo)準(zhǔn)差的均值則反映區(qū)域內(nèi)部是否均一［3］。區(qū)域平均標(biāo)準(zhǔn)差的計算公式為

式中，S為區(qū)域平均標(biāo)準(zhǔn)差;N為區(qū)域總數(shù);ni為第 i個區(qū)域的像元個數(shù);Ri為第 i個區(qū)域;f(x，y)為第i個區(qū)域內(nèi)像元灰度值。

當(dāng)r=5，h取h/m在0.1～0.4之間等間隔變化的100個值，可分別得到100個區(qū)域個數(shù)和區(qū)域平均標(biāo)準(zhǔn)差(圖8)。

圖8 h/m對分割區(qū)域數(shù)目和區(qū)域平均標(biāo)準(zhǔn)差的影響Fig.8 Influence of h/m on region number and average region standard deviation

隨著h/m的增加，區(qū)域個數(shù)減少，區(qū)域平均標(biāo)準(zhǔn)差則呈增大趨勢。當(dāng)h/m取值偏大時，區(qū)域個數(shù)偏少，整體的分割結(jié)果不是最好，但可以突出某些目標(biāo)。

4 結(jié)論

(1)本文采用最小核值相似區(qū)算法(SUSAN)計算QuickBird圖像的梯度，并采用標(biāo)記控制的分水嶺變換(WT)算法分割圖像，取得了較好的結(jié)果。

(2)SUSAN方法能有效地檢測圖像梯度，對噪聲不敏感;梯度值范圍明確，不因圖像而改變，為后續(xù)處理相關(guān)參數(shù)的選擇提供了便利;亮度閾值容易確定，模板半徑可選，具有很大靈活性，適合于采用WT的遙感圖像分割。

(3)采用基于SUSAN梯度和NDVI的標(biāo)記圖像，利用形態(tài)學(xué)灰度圖像重建方法修改梯度圖像，能夠有效地抑制梯度圖像中大量的局部灰度極小值，提高WT圖像分割的精度。

(4)標(biāo)記控制的WT圖像分割的關(guān)鍵是標(biāo)記圖像的精度。本文采用的是基于SUSAN梯度和NDVI的標(biāo)記圖像;今后的研究將進(jìn)一步考慮分割區(qū)域的其他特征，以獲取更加精確的標(biāo)記圖像。

［1］趙書河.高分辨率遙感數(shù)據(jù)處理方法實驗研究［J］.地學(xué)前緣，2006，13(3):60 －68.

［2］Baatz M，Sch?pe A.Object－ oriented and Multi－ scale Image Analysis in Semantic Networks［C］//Proceedings of the 2nd International Symposium on Operationalization of Remote Sensing，1999.

［3］肖鵬峰，馮學(xué)智，趙書河，等.基于相位一致的高分辨率遙感圖像分割方法［J］.測繪學(xué)報，2007，36(2):146 －151.

［4］Neoh H S，Hazanchuk A.Adaptive Edge Detection for Real－ time Video Processing using FPGAs［C］.GSPx 2004 Conference，2004.

［5］Beucher S，Lantuéjoul C.Use of Watersheds in Contour Detection［C］//Proc of International Workshop on Image Processing，Realtime Edgeand Motion Detection/Estimation，Rennes，1979:1－12.

［6］Gonzalez R C.Digital Image Processing［M］.2nd ed.New York:Prentice Hall，2002.

［7］Canny J.A Computational Approach to Edge Detection［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，1986，8(6):679－698.

［8］Stephen M S，Brady J M.SUSAN—A New Approach to Low Level Image Processing［J］.International Journal of Computer Vision，1997，23(1):45 －78.

［9］Marusic B，Skocir P，Juriji T，et al.Application of the SUSAN Filter to Wavelet Video － coding Artifact Removal［J］.Int J Electronics and Communications(AEU)，2006(60):56 －64.

［10］Mao W G，Zhang T W，Wang L.Detection of Text in Images Using SUSAN Edge Detector［J］.Journal of Harbin Institute of Technology(new series)，2005，12(1):34 －40.

［11］Meyer F，Beucher S.Morphological Segmentation［J］.Journal of Visual Communication and Image Representation，1990，1(1):21－46.

［12］Beucher S.The Watershed Transformation Applied to Image Segmentation［C］//Proceedings of the 10th Pfefferkorn Conference on Signal and Image Processing in Microscopy and Microanalysis，1991:299－314

［13］Vincent L.Morphological Grayscale Reconstruction in Image Analysis:Applications and EfficientAlgorithm ［J］.IEEE Transactions on Image Processing，1993，2(2):176 －201.

［14］Soille P.Morphological Image Analysis:Principles and Applications［M］.Berlin:Springer，1999.

Segmentation of the High Spatial Resolution Remotely Sensed Imagery Based on SUSAN

XUE Qiao，ZHAO Shu－h(huán)e
(Department of Geographic Information Science，Nanjing University，Nanjing 210093，China)

The SUSAN(Smallest Univalue Segment Assimilating Nucleus)method is used to detect gradient features from QuickBird imagery，and then the imagery is segmented using marker－ controlled WT(Watershed Transform)，and the segmentation result is satisfactory.The SUSAN method detects gradients well.It is not sensitive to noise and the values of the gradients are in a definite range and do not change with images，which offers convenience in selecting parameters in the later processes.The method is flexible because it is easy to choose the illumination threshold and the size of SUSAN matrix is not fixed.Based on the marker derived from both SUSAN gradients and NDVI，the gradients are modified using morphological grayscale reconstruction method，which efficiently constrains much local minima of the gradients and improves the segmentation precision.

Smallest Univalue Segment Assimilating Nucleus(SUSAN);Watershed Transform(WT);QuickBird imagery;Image segmentation

TP 751.1

A

1001－070X(2011)04－0037－05

2010－12－31;

2011－02－13

國家自然科學(xué)基金項目(編號:40501047)資助。

薛 峭(1987－)，男，碩士研究生，主要從事遙感與地理信息系統(tǒng)方面的研究。

趙書河，E－mail:zhshe@163.com。

(責(zé)任編輯:劉心季)