基于分水岭变换的高分辨率遥感影像分割

缪晓冬 陈曾敏

(南通市测绘院有限公司, 江苏 南通 226000)

0 引言

遥感技术能在不接触目标的情况下获取其几何、物理、位置等信息,如何更好地从遥感影像中获取信息是近几十年来遥感研究的主要挑战[1]。传统的基于像元的遥感影像处理方法都是建立在遥感影像光谱信息极其丰富、地物间光谱差异较为明显的基础上的。而随着传感器技术的发展,遥感影像的空间分辨率变得越来越高,高分辨影像的分析存在着一些问题,即地物内部像元的差异可能非常大,基于像元的影像分析方法无法充分利用高分辨率遥感影像上的信息,所以对于高分辨率遥感影像,传统的基于像元的分析方法不再适用,而基于对象的高分辨率遥感影像分析方法更加适合[2-3]。

基于对象的影像分析方法通常包括两个步骤：一是影像分割,二是特征提取和分类。其中最关键、最基本的步骤是影像分割。影像分割得到匀质的、不重叠的图像对象,然后根据各个对象的属性(光谱、纹理、形状等)进行分类对象的特征提取,分类的准确性主要取决于影像分割的质量。分水岭变换是一种基于区域的方法[4],分水岭变换具有局部一致性特性,分割结果保留了更加精细的细节,分水岭变换最大的优势是可以提取出完整的边缘,不需要边缘连接的过程。但是,根据分水岭变换的定义,影像上只要有一个局部极小值点,就会产生一个分割区域,整个影像上局部极小值点非常多,再加上噪声影响,不管是将其应用在原始影像图上,还是影像梯度图上,都会产生严重的过分割现象。

Giuseppe Masi[5],Giuseppe Scarpa[6],Raffaele Gaetano[7],Giovanni Poggi[8]和Luisa Verdoliva[9]5人逐步提出了一种改进的分水岭算法,从使用形态学标记的边界标志填充算法(edge mark fill,EMF),到进一步引入光谱标记的EMF+算法,最后到能够处理同时包含全色数据和多光谱数据的高分辨率遥感影像的高分遥感边界标志填充算法(multispectral remote-sensing edge mark fill,MR-EMF),三种算法主要基于标记控制来改进原始分水岭算法的过分割现象。栗敏光等人[10]采用形态学混合开闭重构滤波对原始影像进行平滑处理,以此来抑制原始分水岭算法的过分割现象;杨军等人[11]提出了一种基于巴氏系数和灰度共生矩阵的区域合并方法来改进分水岭算法的分割结果;闫鹏飞等[12]针对传统分水岭分割算法易受噪声干扰且分割尺度参数难以自适应选择的问题,提出了一种尺度自适应的高分辨遥感影像分水岭分割方法;闫利等[13]采用邻域自适应阈值、区域合并方法对种子点选取过程进行优化,提出基于多尺度标记和边缘约束的多尺度分割算法。但是在许多学者的研究中,算法不够完整,缺乏预处理及后处理步骤,本文针对分水岭变换的过分割问题,提出一种改进的分水岭变换的高分辨率遥感影像分割算法,来补充基础算法缺失的预处理及后处理步骤,从而达到抑制原始分水岭变换过分割现象、提高分割精度的目的。

1 算法流程

本文主要以基于边缘距离变换的形态学标记及光谱标记控制的分水岭算法为基础,结合形态学混合开闭滤波操作、基于巴氏系数和灰度共生矩阵的区域合并方法,来补充基础算法缺失的预处理及后处理步骤,研究了算法各步骤对分割结果的影响,并根据基于对象匹配的定量分析方法,对本文算法与易康软件(eCognition)中的多尺度分割算法进行定量对比实验分析,算法流程图如图1所示。

1.1 形态学混合开闭滤波

分水岭变化直接作用于原影像或者梯度图像上时,由于原始高分辨率遥感影像像元灰度值往往很复杂并且还时常伴有噪声,会造成分割的轮廓偏移、产生“伪区域”以及导致过分割。为了消除噪声并削弱同类区域内部细节,要对原始影像进行滤波平滑处理。而与传统的线性滤波方式不同,基于数学形态学的非线性滤波方法在消除噪声、削弱细节的同时能保留对象的几何形状特征,更适合与分水岭分割相结合。用于混合开闭重构操作的结构元素的形状和尺寸对影像去噪和平滑也有影响,这里使用圆形结构元素,为了与之后的形态学标记提取使用的圆形结构元素一致,其半径需小于影像中最小分割对象的最小外接圆的半径,防止损失最小待分割地物对象。

1.2 基于边缘的标记控制的分水岭分割

分水岭变换最早来自地貌学基本概念,将灰度图像视作地形图,每个点的灰度值作为它的高程,其基本思想为：在地形图上的每个区域极小值打一个洞,让整个区域从上往下匀速降入水中,直至其被完全淹没,这期间,当不同汇水盆地上升的水汇集时,在汇集地建立水坝来阻止汇集,最终建立的水坝即为分水岭,本文调用矩阵实验室(matrix laboratory，MATLAB)中watershed函数实现分水岭变换。

1.2.1Canny边缘检测算法

1.2.1.1消除噪音

传统Canny边缘检测使用高斯滤波器对图像进行平滑去噪处理,(2k+1)×(2k+1)尺寸的高斯核Hij构造如下：

(3)

式中,σ为标准差;k为不小于0的常数。高斯核的尺寸会影响边缘检测的效果,尺寸越大,对噪声的灵敏度越低,边缘定位误差也会增加;同时,标准差σ也会影响边缘检测效果。

1.2.1.2计算梯度

使用索伯算子(Sobel)得到水平方向和垂直方向的一阶导数值(Gx,Gy),梯度和梯度方向计算如下：

式中,GX和Gy分别为垂直方向和水平方向的索伯滤波器;G为梯度;θ为梯度方向。将连续的梯度方向划分为四个离散的方向：水平、垂直以及两个对角线方向,为下一步非最大化抑制做准备。

1.2.1.3非最大化抑制

此步骤起到细化边缘的作用,对每个像素进行扫描,检查像素是否是在其梯度方向上附近的局部最大值,若是,则保留为边缘;若不是,则赋梯度值为0去除。

1.2.1.4双阈值滞后检测

设置一高一低两个梯度阈值,对边缘像素进行筛选：若梯度值大于高阈值,则为边缘;若梯度值小于小阈值,则舍去;若梯度值介于两者之间,判断其是否与大于高阈值地边缘像素相连,若是则保留为边缘,若不是则舍去。

1.2.2形态学标记方法

1.2.2.1边缘检测

对原始图像进行Canny边缘检测,得到边缘图。

1.2.2.2距离变换

将边缘图(二值图,1代表边缘,0代表非边缘)的像素值替换该像素到最近边缘的距离,得到距离变换图,并取负得到最终进行分水岭分割的数字高程模型(digital elevation model,DEM),此时边缘的像素值为极大值0,距离边缘越远,像素值越小。

1.2.2.3寻找极小值标记种子

寻找DEM的局部极小值,设置为起始标记种子。

1.2.2.4标记膨胀

对种子点进行形态学膨胀操作,膨胀的结构元素为圆形,其半径为种子点到最近边缘的距离再减去一个缩进距离Δ(默认值为3个像素),当两个种子点膨胀后相交时,将其归并为一个标记。

1.2.3光谱信息提取与形态学标记合并算法

形态学标记的生成完全基于原始影像边缘和形态学准则,而忽略了影像的光谱信息,因此利用光谱信息提取光谱标记,并最后将其与形态学标记合并会取得更好的分割效果。同形态学标记的提取一样,为了保留边缘,将边缘像素及其周围像素剔除,在剩余区域提取光谱标记,并根据光谱标记连通区域的光谱方差对较大方差区域进行分割。

1.2.3.1生成光谱区域

对边缘进行膨胀操作,再取反,使得生成光谱标记的区域不包含边缘及其相邻区域。

1.2.3.2提取初始光谱标记

通过分水岭拟合标记侵蚀(watershed marked erosion,WME)使光谱标记的每个连通区域与基于边缘的分水岭分割保持一致。

1.2.3.3 树型马尔科夫随机场(tree shape Markov random field,TS-MRF)

统计各个光谱标记的连通区域的光谱均方差,对均方差最大的1%(默认值)个连通区域进行TS-MRF分割为两类。

1.2.3.4分水岭拟合标记侵蚀

再对1.2.3.3产生的光谱标记使用WME,得到最后的光谱标记。1.2.3.1和1.2.3.2都是为了使得最后分割图保留边缘信息,1.2.3.3以均方差作为判定依据,对均方差较大的光谱标记连通区域使用TS-MRF分割为两类,以此更好地将不同地物分割出来。

1.3 基于巴氏系数和灰度共生矩阵的区域合并

巴氏系数是对两个离散统计样本重叠度的近似计算,在灰度影像上即灰度分布直方图的相似度计算,区域1与区域2之间的巴氏系数B(R1,R2)定义如下：

(6)

式中，L为影像的灰度级；p(i)为灰度为i的像素数量占总像素的比例,巴氏系数越接近1,两区域越相似,同类地物的概率越大。

灰度共生矩阵是对灰度值之间空间关系的统计,用来描述影像的纹理特征,其矩阵维度一般为影像的灰度级,每个灰度共生矩阵元素记为：Pd(i,j),其中d表示在方向,即统计在该方向上的灰度值空间关系,(i,j)为一对灰度值,Pd(i,j)具体的值定义为在方向d上灰度值对(i,j)出现的次数,d包括方向和步长,根据灰度共生矩阵,可以定义出14个纹理特征标量,相关研究分析发现只有角二阶矩(angular second moment, ASM)、熵(entropy, ENT)、对比度(contrast, CON)和相关性(correlation, COR)4个纹理特征与人类视觉感知特性有明确的对应关系[15],因此本文采用这4个特征量来描述纹理信息,作为相邻区域合并的一条准则。角二阶矩反映图像灰度分布均匀程度和纹理粗细程度,熵反映纹理的复杂程度或非均匀程度,对比度反映了图像的清晰度和纹理沟纹的深浅程度,相关性反映的是纹理的方向性。计算这些纹理特征标量时需将灰度共生矩阵归一化(每个元素除以矩阵所有元素之和)。

2 实验分析

对上述完整的算法流程进行实验分析,定性、定量地给出其性能评价,并将其与基于距离变换的分水岭分割算法以及商业软件eCognition Developer中使用的多尺度分割技术进行比较。

2.1 数据介绍及预处理

使用公开数据集Zurich_dataset_v1.0进行实验分析。该数据集获取自2002年8月份瑞士苏黎世地区的20幅QuickBird影像,包括红、绿、蓝及近红外四个波段,并全色锐化至0.6 m的分辨率。同时,还包括20张相应的地物类别标注图,共注释了8类地物：道路、树木、裸土、铁路、建筑物、草地、水体和泳池,但并不是每个像素都被标注。本文选取第十二幅(782×622)进行分析,记为zh12,如图2所示。

图2 zh12实验原始影像

图像的预处理主要是针对公开数据集Zurich_dataset_v1.0提供的地物类别标注图,地物分类图中同一类别的像元共有一个标注,而影像分割与分类不同,需要将同一类别的不同对象赋予不同的标注。因此对分类图进行处理,使得每一类的不同对象分开,赋予不同的标注。

2.2 参数设置及定量精度评价方法

使用基于区域的评价标准,其主要概念为：将算法所得分割图中的每个对象与真实分割图中的每个对象进行匹配,并根据区域判定标准将真实分割图的每个对象分为四类。

2.2.1正确分割(CS)

若真实分割对象(G)与算法分割对象(S)的重叠区域大于各自的75%(α),则将真实分割对象记为正确分割：

(11)

2.2.2过分割(OS)

若真实分割对象被多个算法分割对象所分割,且满足下式,则将真实分割对象即为过分割。

(12)

2.2.3欠分割(US)

若算法分割对象被多个真实分割对象所分割,且满足下式,则将真实分割对象记为欠分割。

(13)

2.2.4漏分割(MEs)

除去上述三类真是分割对象后,剩下的即为漏分割。

2.3 结果分析

下面选取本文算法、基于距离变换的分水岭算法、eCognition-30和eCognition-120得到的分割结果进行展示,如图3所示。

(a)本文算法

统计分割对象尺寸与像元数的分布情况,如图4所示,横坐标x表示所有尺寸处于(2x-1,2x]分割对象,纵坐标y表示这些分割对象的像元数量之和。可以看出,本文算法得到的分割结果尺寸分布比较均匀,大尺寸、小尺寸对象都有,最后还有个明显的最大值,这是因为原影像左侧的水体被很好地分割出来;而基于距离变换的分水岭算法得到的分割结果的尺寸都偏小,大值很少;而eCognition软件的分割结果的尺寸分布与尺度参数相关性很大,都是分布在一个峰值尺寸左右,尺度参数越大,峰值尺寸越大。

图4 分割对象尺寸与该尺寸像元数关系

根据对象匹配原则,计算各个算法的基于区域的精度评价得分,如表1所示,各项得分是根据被标记为CS、OS、US和MEs的真实分割区域面积和与所有参与对象匹配的分割对象面积和之比,并且不是所有参与对象匹配的分割对象都有对应的真实分割对象,第六列为各算法得到的分割对象数量。

从表1可以看出,本文算法的CS(正确分割)得分是最高的,过分割得分处于中间水平,欠、漏分割很少,基于距离变换的分水岭分割算法、eCognition-20、eCognition-30和eCognition-50过分割现象很严重,eCognition软件的后三项分割得分以及分割数和尺度参数相关性很高,尺度参数越大,OS得分越小、US得分越大、MEs得分越大,分割数量越少,即过分割现象越少、欠分割、漏分割现象越多,而CS得分呈现先上升再下降的情况,最优尺度参数在80附近。

表1 各分割算法的对象匹配分数

选取的四个算法所得分割结果的对象匹配效果从图5可以看出,本文提出的算法很好地把原影像左侧的水体分割出来,其余三种算法均对水体得到过分割结果。并且,尺度参数为120的eCognition算法的漏分割现象很严重。

图5 各算法分割结果的分割匹配图

3 结束语

本文以基于边缘距离变换的形态学标记及光谱标记控制的分水岭分割为基础,加入了形态学混合开闭滤波预处理操作以及基于巴氏系数及灰度共生矩阵的区域合并分割后操作,在保留了边缘信息的情况下,充分利用了影像的光谱信息以及纹理信息,大大改善了原始分水岭分割的过分割现象,提高了建筑物这种具有明显边界的对象的分割精度,以及植被区域的分割效果也显著增强。

根据基于对象匹配的定量评价方法,在公开数据集及自制标注数据上将本文算法与eCognition软件中的多尺度分割算法进行对比,展现出了本文算法的优势：①总体分割视觉效果更好,最终分割对象数量适中;②正确分割的定量得分较高;③初始参数设置较少(eCognition软件中的多尺度分割算法需要耗时的试错过程来寻找最佳分割尺度参数)。