基于最小核值相似区算法的高分辨率遥感图像分割方法

薛 峭，赵书河

(南京大学地理信息科学系，南京 210093)

基于最小核值相似区算法的高分辨率遥感图像分割方法

薛 峭，赵书河

(南京大学地理信息科学系，南京 210093)

采用最小核值相似区算法(Smallest Univalue Segment Assimilating Nucleus，SUSAN)计算QuickBird图像的梯度，并采用标记控制的分水岭变换(Watershed Transform，WT)算法分割图像，取得了较好的结果。SUSAN方法能有效地检测图像梯度，对噪声不敏感;梯度值范围明确，不因图像而改变，为后续处理相关参数的选择提供了便利;亮度阈值容易确定，模板半径可选，具有很大的灵活性;适合于采用WT的遥感图像的分割。采用基于SUSAN梯度和NDVI的标记图像，利用形态学灰度图像重建方法修改梯度图像，能够有效地抑制梯度图像中大量的局部灰度极小值，提高WT图像分割的精度。

最小核值相似区算法(SUSAN);分水岭变换(WT);QuickBird图像;图像分割

0 引言

高分辨率遥感图像中地物景观的结构、形状、纹理和细节等信息非常突出，基于像元分类的传统算法难以从高空间分辨率遥感数据中提取所需的信息［1］，因此面向对象的高分辨率遥感图像处理方法得到不断发展［2］。图像分割是高分辨率遥感图像信息提取和面向对象分类的重要环节［3］。

实现图像分割的方法包括基于边界的方法和基于区域的方法，前者检测局部变化，后者则检测像元和区域的相似性［4］。分水岭变换(Watershed Transform，WT)属于第二种方法。基于数学形态学原理的WT方法被广泛地应用于灰度图像分割中。最初，WT方法被用来分割气泡图像和金相显微镜图像，从图像中提取出气泡和显微镜下的目标［5］。将图像看作地表，灰度值为其高程值，在每个局部最低点位置打一个洞并且让水从洞中涌出，从低到高淹没整个地形;当处在不同汇水盆地中的水要聚合到一起时，修建大坝阻止其聚合，大坝就对应图像的分割线［6］。分割线是闭合的，与图像灰度的等值线一致，闭合的分水岭线包围的区域就是图像分割中的目标。

直接对图像进行WT分割将得到2个不好的结果:①很多分水岭线没有位于目标的边缘。这样，即使汇水盆地大小合理，也不能保证分水岭线包围的是完整的目标(理想的情况是目标的边缘高程值即灰度值高于非边缘的灰度值，分水岭线位于高处即形成目标边缘);②由于过分割，区域破碎，目标内部存在大量高程极小值，WT方法会形成众多面积微小的汇水盆地。因此，在WT分割中常使用梯度图像来保证目标与其边缘的高程差异，并且采取措施抑制梯度图像目标内部大量高程极小值，以消除微小的汇水盆地。WT图像分割的性能依赖于计算图像梯度的方法和抑制梯度图像极小值的方法［3］。典型的梯度检测算子都用核算子来计算正交方向上的坡度值的强度，然后计算出总的梯度强度值［4］，但上述计算容易受到噪声的影响。Canny梯度检测属于线性滤波梯度检测算法，检测步骤包括高斯平滑去噪、垂直和水平方向上的高斯算子边缘检测、细化边缘的非极大值抑制和双阈值处理(高阈值得到初步的梯度，低阈值保证初步梯度的连续性［7］)，其结果容易受到高斯滤波器标准差和高、低2个阈值的影响。1997年Stephen等［8］提出用最小核值相似区算法(Smallest Univalue Segment Assimilating Nucleus，SUSAN)检测图像的边缘，其实质是根据像元周围局部区域内与该像元亮度值相近的像元个数来判断该点是否为边缘。不用对图像求导，可直接对图像灰度进行运算，效率高，定位准确，具有积分特性，对噪声不敏感［9－10］。

本文采用SUSAN方法计算QuickBird图像梯度，并采用基于梯度和NDVI提取的标记图像对梯度进行重建，抑制大量的梯度图像的灰度极小值;最后进行图像的WT分割，取得较好的分割效果。

1 实验数据

所用实验数据为南京地区的QuickBird图像，包括蓝(450～520 nm)、绿(520～660 nm)、红(630～690 nm)和近红外(760～900 nm)4个多光谱波段，空间分辨率为2.44 m;全色波段(450～900 nm)，空间分辨率为0.61 m。本文截取的数据是512像元×512像元的全色图像及与其对应的多光谱图像。典型地物包括河流、建筑物和植被等(图1)。

图1 实验区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)，模板除去中心核外包含的区域大小m是36个像元(图3中的暗色像元)。

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

若模板内其他像元的灰度值与模板中心像元的灰度值之差小于一定的亮度差阈值t，就认为二者相似。由相似的像元组成的区域称为任意像元的核值相似区(USAN)［10］。USAN包含了重要的图像结构信息:在非边缘处，USAN值较大(最大可为模板所含像元个数);越靠近边缘，USAN值越小。在图2中的1处，模板位于角点，USAN值为1/4模板大小;当模板的中心位于边缘时，USAN值为1/2模板大小(如图2中的3处);在图2中的5处，USAN值则小于1/4模板大小。图像中每一个像元对应于一个USAN值，因此通过USAN值就能检测出图像的边缘。用上述模板遍历图像，在任意像元位置，计算该像元的USAN值为

式中，u为像元的 USAN值，u(r0，r)=t为亮度差阈值，表示模板中心核像元亮度值与其他像元亮度值的差别，差别小则认为二者相似。Stephen等［8］讨论了t值的选取与图像灰度范围之间的关系，对于8 bit灰度图像建议t取25。实验表明，对于QuickBird图像而言，t取值在25～40之间较合适，本文取t=30。

实际采用 u(r0，r)计算像元间的相似度，函数值在0和1之间。遍历图像完成计算后，得到一个新的浮点型灰度图像，称作该图像的USAN图，它反映每个像元的同质区的大小(像元灰度最大值约为36，最小值接近0)，值大的像元为非边缘，值小的像元为边缘。

为了得到适合WT分割的梯度图像，对USAN图像进一步做2步处理:①进行灰度值求反，用该图像的最大值减去每个像元的值作为该像元的值，使得目标边缘处高程值大，以适合WT的图像分割;②将值小于此时图像最大值的1/4倍的像元赋值为0，以初步抑制弱小梯度。得到的梯度图像被称作SUSAN梯度图像。

2.2 SUSAN梯度图像获取

取研究区全色图像，将4个多光谱波段的空间分辨率重采样成与全色波段一样大小。分别计算上述5幅图像中每个像元的SUSAN梯度值，取最大值作为该位置上的最终梯度值。图4为模板半径r=5时得到的梯度图像。大多数地物的内部与边缘有较好的可分性，但河流旁的树木区域内部与边缘的差异不够明显，需要进一步的处理。

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

SUSAN梯度的重要特点是它的取值范围只与模板半径r有关，为后续处理相关参数的选择提供便利;同时r可灵活选择，以适应不同的需求。下文将用到这个特点。

3 WT图像分割

常用基于标记控制的WT方法解决过分割问题［11］，前提是要知道被分割目标的大致位置［12］。在梯度图像中，目标的位置在梯度值较小处，梯度值较大处则是目标的边缘。然而梯度较小的区域中仍存在较多的梯度极小值，会造成分割的破碎。反映图像中目标大概位置的另一种图像叫做标记图像。典型的标记图像是目标处为1、非目标处为0且与原图像大小相同的二值图像。根据梯度图像和标记图像，利用形态学灰度重建原理进行梯度修改，可得到新的梯度图像［13］。此时，对应于标记位置的梯度变得均一，因此分水岭线不会破坏标记范围内的目标。标记的作用是确定目标的大致位置(以便于对梯度图像进行重建)，与目标越接近的标记越好。

3.1 标记提取

图像的梯度特征和光谱特征常用于提取标记图像。扩展最小变换(Extended－minima Transform，EM)是常见的基于梯度的标记提取方法［14］，它把目标的大致位置处标记为1，其余地方标记为0。梯度图像G经过高度阈值为h的扩展最小变换运算为

式中，E为输出的与G大小相同的二值图像;h为大于零的高程值。

若局部极小区域的最大高程差大于h，则在输出图像对应的局部极小区位置上输出1，否则为0。在得到的标记图像中，值为1处则对应于目标的大致位置。

同时采用基于光谱特征的标记图像提取方法。植被内部与其边缘梯度值的差异较小，但植被的光谱特征明显，可采用NDVI来确定植被的范围，利用NDVI提取植被的标记图像。将满足NDVI＞0.2的区域作为植被的标记，与扩展最小变换提取的标记做“或”运算，满足其中之一者即可作为最终的标记(最终的标记图像如图5)。

图5 标记图像Fig.5 Marker image

利用灰度图像重建方法对梯度图像修改之后，梯度图像中的标记处变成了均匀的极小区域(如图6所示)。灰度重建可采用多种计算方法［13］，本文采用了形态学灰度图像重建方法。

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

采用SUSAN梯度，再获取基于梯度和NDVI的标记图像，用标记控制的WT分割方法对QuickBird图像进行分割(分割结果如图7)。整个过程在MATLAB R2007b软件环境下实现。

图7 分割结果Fig.7 Segmentation result

3.2 精度分析与参数选择

由图7可以看出，大多地物都能被较好地分割开来，但个别建筑物被分成了几个区域，这与它们的复杂结构有关。EM中对高度阈值h的选择与图像的梯度大小有关，梯度范围与模板半径r有关(模板半径为 r，梯度范围约是0～πr2，最大值 m约为πr2)。当h满足 h/m 在 0.1～0.4之间时，能取得较好的分割结果。本文中 r=5，h=20，h/m≈0.25，区域个数=259。随着r的增加，梯度的宽度也会增加。取r=5，当采用基于NDVI的标记提取方法时，植被区域被选择成标记，成为灰度极小区域;区域边缘的外侧仍然存在着梯度，即“山谷”的周围存在包围该区域的“山脊”，成为区域的分水岭线。此时r取值稍大更合适，而在其他图像的分割中可选择较小半径的模板以提高计算速度。

对分割精度的评价采用区域个数与区域标准差的均值来衡量——区域个数反映分割结果是否破碎，区域标准差的均值则反映区域内部是否均一［3］。区域平均标准差的计算公式为

式中，S为区域平均标准差;N为区域总数;ni为第 i个区域的像元个数;Ri为第 i个区域;f(x，y)为第i个区域内像元灰度值。

当r=5，h取h/m在0.1～0.4之间等间隔变化的100个值，可分别得到100个区域个数和区域平均标准差(图8)。

图8 h/m对分割区域数目和区域平均标准差的影响Fig.8 Influence of h/m on region number and average region standard deviation

随着h/m的增加，区域个数减少，区域平均标准差则呈增大趋势。当h/m取值偏大时，区域个数偏少，整体的分割结果不是最好，但可以突出某些目标。

4 结论

(1)本文采用最小核值相似区算法(SUSAN)计算QuickBird图像的梯度，并采用标记控制的分水岭变换(WT)算法分割图像，取得了较好的结果。

(2)SUSAN方法能有效地检测图像梯度，对噪声不敏感;梯度值范围明确，不因图像而改变，为后续处理相关参数的选择提供了便利;亮度阈值容易确定，模板半径可选，具有很大灵活性，适合于采用WT的遥感图像分割。

(3)采用基于SUSAN梯度和NDVI的标记图像，利用形态学灰度图像重建方法修改梯度图像，能够有效地抑制梯度图像中大量的局部灰度极小值，提高WT图像分割的精度。

(4)标记控制的WT图像分割的关键是标记图像的精度。本文采用的是基于SUSAN梯度和NDVI的标记图像;今后的研究将进一步考虑分割区域的其他特征，以获取更加精确的标记图像。

［1］赵书河.高分辨率遥感数据处理方法实验研究［J］.地学前缘，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］肖鹏峰，冯学智，赵书河，等.基于相位一致的高分辨率遥感图像分割方法［J］.测绘学报，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－he
(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

国家自然科学基金项目(编号:40501047)资助。

薛 峭(1987－)，男，硕士研究生，主要从事遥感与地理信息系统方面的研究。

赵书河，E－mail:zhshe@163.com。

(责任编辑:刘心季)