张 毅,张彦峰
(武汉大学遥感信息工程学院,湖北 武汉 430079)
基于光谱-空间特征的遥感影像水体提取方法
张毅,张彦峰
(武汉大学遥感信息工程学院,湖北武汉 430079)
地表水体是人类生存和社会发展不可取代的基本要素,使用遥感影像提取地表水体相比传统实地调研的方式,具有高效、省时、省力的特点。传统的水体提取方法一般仅使用一种或少量类型特征的结合,在处理高分辨率遥感影像时容易出现误提取。针对以上问题,提出一种结合多种光谱和空间特征的自动化水体提取方法。通过两幅高分一号影像上的定性和定量实验结果表明,该方法能够比传统方法更稳定更精确的提取影像中的水体信息。
遥感;水体提取;影像分类;光谱特征;空间特征
地表水体是人类生存和社会发展不可取代的基本要素。地表水体的研究对气候模型,农业稳定性,河流变化,湿地监管,地表水体调研和管理,环境监测等领域起着重要的作用[1,2]。准确地绘制出地表水体的空间分布是水文相关领域研究的前提,也是制定相关政策非常关键的一步[3]。传统的以实地调研为主的水体提取方式需要耗费大量的人力物力资源,且耗时较长难以获得及时的信息。遥感是一种利用非接触手段获取信息的方式,通过从遥感影像中解读水体信息,大大降低了获得地表水体信息的成本。
目前的国内外研究中,遥感影像提取水体的主要研究主要集中在特征提取和水体提取算法两个方面。合适的特征对正确的从遥感影像提取水体至关重要,目前在水体提取中研究的侧重点为水体特征的研究主要包括:
1)光谱特征:遥感影像获取的原始信息进行辐射定标、大气校正后的结果,是水体提取过程中最基本的特征,大部分水体提取过程中均会使用光谱特征。总体上来说,水体在蓝绿两波段上呈现反射特性,在近红外波段上呈现吸收的特性。根据水体和背景像素的光谱反射率差异,得到基于规则的方法得到影像中的水体,或使用光谱反射率作为特征训练分类器如支持向量机(SVM)[4],对影像进行二值分类得到水体和非水体信息。
2)水体指数:通过分析不同地物对遥感影像波段的响应情况对波段进行组合、代数运算,从而增强水体信息,并抑制其他地物信息。常见的水体指数有:NDWI[5],MNDWI[6],AEWI[7],HRWI[4]等。
3)纹理特征:使用纹理特征或纹理特征与光谱的结合,提取水体信息。文献[8]在水体提取问题上提出了方向方差算子和纹理表相结合的多纹理特征结合的水体提取算法,在高分辨率灰度影像的实验上得到了很好的提取结果。文献[9]比较了仅使用颜色特征和颜色特征加灰度共生矩阵角二阶矩特征的提取效果,实验表明,增加纹理特征后,提取结果有了很大的提高。
4)颜色特征:对影像的红绿蓝波段反射率进行颜色空间变化如HIS变换、CIELAB变换等,利用变化后的颜色空间特征作为水体提取的特征。文献[10]研究了HIS变换后的色调(H)、饱和度(S)和亮度(I),对高分辨率遥感影像进行水体提取。文献[9]通过CIELAB变换,使用亮度分量(L)进行分割,并使用颜色分量(a和b)和亮度分量的纹理特征构成特征向量提取航空影像上的水体。
5)形状特征:通过分割的方式得到对象,提取对象的形状特征,利用水体和非水体对象在形状特征上的差异,提取影像中的水体对象。文献[11]使用紧致度(Compactness)、主要点(Criticalpoints)和对称性(Symmetry)三个形状特征成功区分了影像中的阴影和水体。
目前研究中,一般水体提取方法仅使用以上一类或两类特征,然而,由于地表水体的复杂性,以上一种或两种特征的结合并不能很好的从影像中提取水体信息。首先,“同物异谱”导致同样的水体可能具有不同的光谱、水体指数、颜色特征,而“同谱异物”现象使得水体和背景信息在这些特征上存在混淆。其次,一般的纹理特征在计算时需要进行取窗口的操作,在水体边缘上的纹理信息会产出模糊的现象,无法完全反映水体的特点。最后,自然界中的水体往往具有不确定的形状,因此,很多情况下形状特征无法准确区分水体和背景地物。针对以上问题,本文提出一种基于光谱-空间特征的遥感影像水体提取算法,结合了光谱、水体指数、像元长度指数和面积特征,多种类型特征的结合更好的描述了遥感影像上水体的特征,鲁棒地提取遥感影像中的水体信息。
本文方法使用了光谱、水体指数、像元长度特征和面积特征,通过结合光谱和空间特征,一方面减轻了“同物异谱”和“同谱异物”问题对水体和非水体像素的影响,另一方面引入的像元长度特征综合影像中像素的局部结构特征,有效的反映了像元邻域内的相似性,且避免了求取纹理特征时滑动窗口操作可能产生的问题。本文方法总体上可以分为两个部分,第一部分使用双阈值分割去除影像中的大部分非水体像素,第二部分将第一步得到的水体待选像素聚集成水体待选对象,并使用面积特征和水体指数去除非水体对象,最后对得到的水体对象进行区域生长,得到精确的水体边界信息。方法的整体流程如图1所示。
图1 水体提取方法整体流程图
1.1辐射定标
首先对遥感影像进行辐射定标预处理,根据公式(1)将卫星各载荷的通道观测技术值DN转换为卫星载荷入瞳处等效表观辐射亮度:
式中:Gain为定标斜距;DN为卫星载荷观测值;Bias为定标截距,单位为;以上参数均由实际卫星数据给出;
再根据如下公式,得到大气顶端的表观辐射量:
其中,L为上一步得到的入瞳处的表观辐射量,d为年积日即影像获取时间为一年中的第几天,θ为太阳高度角,以上数据可从影像的元数据获得,ESUN为太阳表观辐射率均值,不同的卫星传感器会给定相应的取值。
1.2基于光谱-空间特征的双阈值分割
本步骤的目标为结合光谱和像元的局部结构特征通过两步阈值分割处理,去除不符合要求的非水体像素。水体指数从光谱的角度反映了水体和非水体像素的差异,其取值范围为[-1,1],水体的NDWI取值一般为正,且取值越大该像素为水体的可能性也越大,水体指数的计算公式如下所示:
其中,green为遥感影像上的绿波段,NIR为近红外波段的辐射量;
其次,本文引入像元长度指数 (PixelLength Index,PLI),衡量水体内部的相似性。像元长度指数是像元结构特征集合[13](StructuralFeatureSet,SFS)中的一种特征,该特征考虑了相邻像元的相似性,综合了像素间的上下文信息,具有减少同质性区域的噪声和光谱变化的特性。像元长度特征反映了像元在其邻域内相似性的最大值。由于同一水体对象内部的光谱差异较小,且不同水体对象之间的形态差异较大,使用像元长度指数既可以反映影像内部水体之间的相似性,也减小了水体对象形态差异造成的像元结构差异。
在获得影像的像元长度指数时,为了加大水体之间的相似性,以及水体和非水体之间的差异性,本文使用拉伸后的水体指数 SNDWI(Stretched NDWI)特征计算像元长度指数。由公式得到水体指数特征(NDWI),并根据公式(4)进行拉伸。
使用SNDWI得到的光谱相似性测度为下式:
其中,phi表示第i个方向线上中心像素和邻域像素之间的异质性,Psncen为中心像素的SNDWI特征,Psncen为邻域像素的SNDWI特征。
对预处理后的影像数据提取水体指数和像元长度指数,通过双阈值分割,得到水体待选像素,具体步骤如下所示:首先,对影像的SNDWI特征计算得到的像元长度特征使用给定阈值进行二值化,从像元的局部相似性,排除了影像中大量的非水体像素。获得上步中通过阈值分割的像素的NDWI特征值,使用大津法对其进行阈值分割,利用水体和非水体像素在水体指数上的差异进一步排除影像中的非水体像素。
1.3基于双阈值分割结果的水体对象提取
经过上一步的处理,得到了通过光谱和空间特征约束的像素的二值结果。首先,根据二值图上像素的连通性,将8邻域内相邻的连接像素聚集为对象,生成水体待选对象。
根据聚集的水体对象的面积进行阈值分割,给定面积阈值,若对象的面积大于给定的阈值则认为该对象为大面积水体待选对象,否则将其作为小面积水体待选对象。给定一个较大的面积阈值,如10000,则得到的大面积水体待选对象为面积很大的平滑区域,且具有较高的NDWI值。由于遥感影像上存在面积大且水体指数高的弱纹理非水体地物概率非常小,故接受大面积水体待选对象为大面积水体对象。而由于小面积的对象中含有背景地物的概率较高,故对小面积水体待选对象将进行进一步的处理以排除背景地物。
获得小面积水体待选对象覆盖范围内所有像素的NDWI特征,根据其直方图形状使用峰谷的方法(peak-and-valleymethod)获得阈值再一次进行阈值分割。为了减小分割时的漏检率,因此选择最左边的谷点对应的NDWI值作为分割的阈值。判断直方图上峰谷的方式为:
1)如果pn≥pn+1&pn≥pn-1,则为直方图上的峰值;
2)如果pn≥pn+1&pn≥pn-1,则为直方图上的谷值;
其中,pn表示直方图上n位置对应的频率。
然而,由上述条件得到的结果往往会包含较多的伪峰和伪谷,因此,通过对直方图进行核函数估计(Kerneldensityestimation,KDE)得到平滑后的概率密度曲线(直方图),从平滑后的结果中提取出峰值和谷值。核函数密度估计是非参数的估计随机变量概率密度函数的一种方法,常被从来实现数据平滑的目的。其计算公式如下所示:
其中,Kh(·)为核函数,本文使用高斯核函数。图2给出了经过第一步处理后剩余像素的NDWI直方图,以及对该直方图进行核函数估计的结果。从图中可以看出,蓝色部分的直方图中存在伪峰值,而平滑后的概率密度曲线上不存在伪峰值,因此,使用平滑后的结果可以得到正确的阈值。
图2 影像NDWI直方图核函数估计结果
如果概率密度估计得到的结果中不包含谷值点,则说明经过第一步的处理后,剩下的小面积水体待选对象中不包含非水体像素,则不对其进行进一步的阈值分割。
由以上两步可以看出,本文方法旨在对大面积的水体使用较小的NDWI阈值进行分割,而对面积较小的水体使用较高的NDWI阈值分割,并且在使用PLI进行阈值分割时去除了大量的与水体具有类似NDWI取值的背景像素。
由于水体边缘像素的PLI和NDWI特征比水体中心像素小,故在上述处理中可能被错误的去除,因此,在得到水体对象后,通过区域生长的方式得到更精确的水体边界。将所有水体边缘像素加入种子点堆栈,若种子点边缘像素与其邻域像素之间异质性测度小于给定阈值,则认为该邻域像素为水体,并将其加入种子点堆栈,重复进行上述操作直到所有种子点均被处理完成。这里使用的异质性测度为光谱角,因为该特征综合考虑了影像的所有光谱特征,计算公式如(7)所示,给定的阈值为thsa。
其中,SA(x,y)表示光谱向量之间的光谱角,xi和yi分别表示种子点像素和邻域像素在i波段上的光谱反射率,n表示影像的波段个数。对于具有可见光和近红外四个波段的高分辨率遥感影像来说,x和y均为四波段波谱反射率组成的四维向量。
本文使用两幅高分一号卫星的8米分辨率多光谱影像进行水体提取实验。高分一号卫星是我国高分辨率对地观测系统重大专项的第一个卫星,具有重要的战略意义和地位。其8m分辨率的多光谱卫星有近红外和可见光四个波段,重访周期为4d,具有较高空间分辨率和较短的回访周期的特点,因此,具有较高的应用价值。本文使用两幅高分一号数据进行实验,通过对比人工勾选的水体提取结果定性和定量的比较了本文提取方法与使用光谱特征的支持向量机(典型的水体提取方法)的提取结果。
本文使用如下方法进行定量的精度评点:对比参考影像使用混淆矩阵统计,得到正确提取的水体像素(TP),错误提取的水体像素(FP),正确提取的非水体像素(TN),错误提取的非水体像素(FN),四类的像素个数。使用统计结果得到水体的生产者精度PA=TP/(TP+FN),生产者精度反映了水体的漏检率,使用者精度UA=TP/(TP+FP),使用者精度反映了水体的误检率,总体精度OA=(TP+TN)/N,其中N为影像中的像素总体个数和kappa系数(如式(8)所示)。
其中,
图3(a)显示了本文实验中所用的第一幅影像(影像一),图3(b)显示了相应的参考水体掩膜。影像一相幅大小为4548*4500,经度范围为120.0874度到 120.5829度,纬度范围为 40.0612度到40.4397度。从影像中可以看出大面积水体的光谱特征并不统一,其大部分面积呈蓝色,水体左下角颜色变深,并且在水体内部可以看到条纹状的差异。影像中的背景地物主要由黄色的裸地、深色山体和亮建筑物组成。
图3 实验影像一(a为伪彩色显示结果,b为参考水体掩膜)
图4 影像一的提取结果(a为SVM结果,b为本文方法结果)
表1 影像一精度评定结果
图4给出了使用SVM和本文方法提取的结果,表1给出了影像一提取结果的定量评价。从提取结果中可以看出,两种方法均较好的提取出了遥感影像中的水体,本文方法得到的大面积水体更为完整,小面积的误提取更少。从定量的结果可以看出,本文方法在四个指标上均优于SVM的提取结果。其中,使用者精度提高了约3%,说明本文方法极大的减小了误检情况的发生。
图5 实验影像二(a为伪彩色显示结果,b为参考水体掩膜)
图5(a)显示了本文实验中所用的第二幅影像(影像二),图5(b)为其相应的参考水体掩膜。影像二相幅大小为4548*4500,经度范围为120.1701度到120.6267度,纬度范围为39.1023度到30.4952度。从影像中可以看出,大面积的水体本身的光谱发生了较大的变化,左下方主要呈现蓝色,而右上方呈现绿色,说明右上方的水体在红波段的反射率较高。不仅如此,影像中存在小面积的暗水体。影像二的背景信息复杂,既包括亮建筑物,也包括暗建筑物,还存在较多的植被、裸地以及阴影信息。
图6给出了SVM和本文方法对影像二的提取结果,表2给出了定量的精度评价结果。与人工勾选的参考影像相比可以看出,SVM和本文方法均有效提取出了影像中的主要水体信息,且在小面积水体上均存在一定的误检问题,但整体上本文方法得到的结果更为完整,且误检的像素更少。从定量评价的角度可以看出,本文的生产者精度比SVM高约10%,使用者精度高约13.5%,整体精度高2.3%,因此,使用本文方法较好的提高了水体的提取精度。
图6 影像二的提取结果(a为SVM结果,b为本文方法结果)
表2 影像二精度评定结果
从总体上了来说,本文方法的提取精度均高于SVM的提取精度,多种特征的结合以及像元长度指数的引入,有效的提高了本文方法在大面积水体提取结果上的完整度,且有效的去除了小面积的噪声。不仅如此,SVM为监督分类方法,使用该方法进行水体提取时,需要人工选择训练样本,并通过多次实验选取合适的参数。而本文方法为自动提取方法,计算速度快。
针对目前遥感影像水体提取中使用的特征存在的多种问题,本文提出了一种基于光谱-空间特征的遥感影像水体提取算法。本文方法结合像元长度指数、光谱指数以及对象的面积特征,从光谱和空间特征的角度区分影像中的水体和非水体像素。通过两幅高分一号数据,定性和定量的比较了本文方法的有效性。从总体上来说,本文方法能够高效自动化的提取遥感影像上的水体。
[1] Du Z,Linghu B,Ling F,et al.Estimating surface water area changes using time-series Land sat data in the Qingjiang River Basin,China[J].Journal of Applied Remote Sensing,2012,6 (1):063609-063609.
[2] Sun F,Sun W,Chen J,et al.Comparison and improvement of methods for identifying water bodies in remotely sensed imagery[J].International journal of remote sensing,2012,33 (21):6854-6875.
[3] Gessner M O,Hinkelmann R,Nützmann G,et al.Urban water interfaces[J].Journal of Hydrology,2014,514:226-232.
[4] Yao F,Wang C,Dong D,et al.High-Resolution Mapping of Urban Surface Water Using ZY-3 Multi-Spectral Imagery[J]. Remote Sensing,2015,7(9):12336-12355.
[5] McFeeters S K.The use of the Normalized Difference Water Index(NDWI)in the delineation of open water features[J].International journal of remote sensing,1996,17(7):1425-1432.
[6] Xu H.Modification of normalised difference water index (NDWI)to enhance open water features in remotely sensed imagery[J].International Journal of Remote Sensing,2006,27 (14):3025-3033.
[7] Feyisa G L,Meilby H,Fensholt R,et al.Automated Water Extraction Index:A new technique for surface water mapping using Landsat imagery[J].Remote Sensing of Environment,2014,140:23-35.
[8] Wang H,Pan L,Zheng H.Multi-texture-model for water extraction based on remote sensing image[C]//Image and Signal Processing,2008.CISP'08.Congress on.IEEE,3:710-714.
[9] ZhaoM,ZhangY,ShangH,et al.Water Area Extraction from RGB Aerophotograph Based on Chromatic and Textural Analysis[C]//GEO Processing 2011,The Third International Conference on Advanced Geographic Information Systems,Applications,and Services:46-52.
[10] 柳稼航.基于视觉特征的高分辨率光学遥感影像目标识别与提取技术研究[D].上海交通大学,2011.
[11] Li B,Zhang H,Xu F.Water extraction in high resolution remote sensing image based on hierarchical spectrum and shape features[C]//IOP Conference Series:Earth and Environmental Science.IOP Publishing,2014,17(1):012123.
[12] 黄昕.高分辨率遥感影像多尺度纹理,形状特征提取与面向对象分类研究[D].武汉:武汉大学,2009.
P407.8