一种基于白帽变换的细小水体信息提取方法

张 强,吴 波

(福州大学地理空间信息技术国家地方联合工程研究中心,空间数据挖掘与信息共享教育部重点实验室,福建 福州 350116)

0 引言

水体的储量和位置与人类水产活动息息相关,快速准确地获取地表水体信息对于区域水土保持监测和水资源的合理开发利用具有重要的现实意义. 由于遥感影像中水体光谱特性较为明显,且在空间分布上有较强的连通性,易于通过遥感影像快速准确地提取[1]. 相关研究人员从不同的角度提出各种水体信息提取方法,包括谱间分析法[2]、 影像分类法[3]、 水体指数法[4-5]、 特征变化法[6]等. 这类方法虽然在宽阔水体或大面积水域提取中取得较好效果,却难以准确提取河面较为狭窄的细小水体. 这是由于开阔水体的光谱较为纯净,其影像特征具有较好的区分度,因而容易识别. 而细小水体通常定义为影像中像元表观宽度小于或等于3个像元的狭窄河流[7],受遥感影像空间分辨率的限制,细小水体大都呈现混合像元,在影像中以混合光谱的形态出现,不再具有稳定的光谱特征,单纯依据光谱信息难以准确提取[1].

为分析细小水体与其它地类特征信息,图1给出6类典型地表覆盖类型(即开阔水体、 细小水体、 植被、 城区、 道路、 山体阴影)在归一化差异水体指数(normalized difference water index, NDWI)[4]和改进的归一化差异水体指数(modified normalized difference water index, MNDWI)[5]两种典型水体指数特征影像中的变化情况. 图中以[均值-标准差,均值+标准差]作为地类的集中分布区间,每类地物选取200个样本进行统计分析. 其中细小水体样本从细小河流中选取； 开阔水体样本从水体的中心区域选取. 其他4类分别为对应地类具有代表性的样本. 从图1可看出,开阔水体的NDWI与MNDWI特征值远高于陆地背景,易与其它地类区分. 而细小水体除了易与植被区分,其特征值与城区、 道路和山区阴影均存在不同程度的混淆. 这表明仅依赖光谱信息难以区分细小水体,需引入其它特征辅助提取细小水体信息[8].

图1 NDWI、 MNDWI特征值统计Fig.1 Statistical feature value of NDWI, MNDWI

基于以上认识,杨树文等[9]在初步提取水体的基础上,结合形态学膨胀和图像细化连接漏提河道. 该方法在漏提较少的河道,较准确地补全了细小水体,但其补全过程忽视了水体原有的纹理信息,导致河道间断较大时误提严重. 李艳华等[8]在光谱信息的基础上引入纹理、 形状等空间特征,通过获取最优分割尺度来分割影像,从而提取细小水体,但其需要人为设定较多的参数,因而不利于水体信息的快速提取,在大范围区域内难以适用. Jiang等[7]利用线性特征模型LFE增强特征影像中的线状水体,有效地提取了细小水体,但提取过程中的多个经验阈值降低了该方法的自动化程度. 为了更加完整、 高效地提取细小水体,本文以中等分辨率的Landsat影像为数据源,在MNDWI初步增强水体信息的基础上,提出一种基于形态学的细小水体特征的后处理方法,并结合最大类间方差[10]的自动阈值选取方法进行双阈值分割,解决现有方法中[2-6]细小水体难以完整提取的问题,实现细小水体准确、 快速的提取.

1 细小水体信息的提取方法

图2 细小水体提取流程 Fig.2 Flow chart of narrow water extraction

综合考虑细小水体在遥感影像中的光谱与纹理特征,并根据水体在空间分布上具有连通的特性,设计一种基于影像线性增强的细小水体提取方法,通过经验阈值与自动阈值结合的双阈值分割,实现细小水体的完整提取,基本的技术流程如图2所示,主要包括以下4个步骤.

1) 影像预处理. 为保证影像光谱信息的原始性,在获取Landsat ETM+/OLI原始数据的基础上,仅对其进行必要的大气辐射校正、 系统几何校正.

2) 细小水体特征增强. 先计算遥感影像的MNDWI,其公式为：

(1)

式中：DNGreen代表绿光波段的灰度值；DNMIR代表短波红外波段的灰度值. 图1表明MNDWI虽然在局部范围内增大了水体与陆地背景间的反差,但细小水体仍易与城区、 道路、 山体阴影等地物产生混淆. 借助线性形态学算子和白帽变换对细小水体进行增强,构建出一种结合形态学的细小水体指数MNWI(morphologicalnarrowwaterindex),其详细步骤总结如下.

Step 1 定义线性结构体. 由于细小水体具有典型的线状特征,方向性显著,而建筑物、 山体阴影具有面状特征. 本文的线性结构体定义为se=strelem(d,s), 其中d表示线性结构体的方向(分别为0°、 45°、 90°、 135°),s表示细小水体的尺度(即宽度在几个像元内的水体被定义为细小水体,在本文中s=3).

Step 2 形态学白帽重构(white top-hat). 为增强宽度小于结构体se的线状地物,对MNDWI影像进行形态学白帽重构,其公式如下所示：

WTH(d,s)=MNDWI-γMNDWI(se)

(2)

式中： γMNDWI(se)表示MNDWI特征影像在结构元素se作用下的开运算结果；WTH则表示在尺度s和方向d作用下的白帽重构结果. 经过此次重构,小于结构元素se的细小水体将高亮显示,大于结构元素se的开阔水体及背景地物将被显著抑制.

Step 3 构建形态学细小水体指数MNWI. WTH重构虽然有效地抑制了无关背景,但可能仍有少量小面积建筑会被增强. 由于建筑物的面状特性,其WTH在4个不同方向上有较高的灰度值,故利用WTH在各个方向上最大、 最小值的反差,抑制面状背景,突出线状细小水体. 据此原理构建MNWI公式如下：

MNWI=WTHmax-WTHmin

(3)

式中：WTHmax为4个方向上白帽重构的最大值；WTHmin为4个方向上白帽重构的最小值.

3) 双阈值分割. 首先在MNDWI选择一个较大的经验阈值,获取水体信息的初始分割. 为确保目标水体不受道路噪声干扰,本文基于多次统计实验,确认经验阈值Te= 0.2是一个较好的阈值. 然后,利用最大类间方差法(Otsu)在MNWI影像中自动确定另外一个阈值,获取MNWI影像的潜在水体区域. 如果潜在的水体区域与经验阈值确定的水体区域在空间上相连通,则将潜在水体区域判定为真实水体,否则该区域为非水体. 这样就初步提取了MNWI影像中的水体信息.

4) 影像后处理. 在初步提取的水体信息中,混有少量道路或桥梁地物,采用归一化建筑物指数(normalizeddifferencebuilt-upindex,NDBI)[11]来剔除这些误提的目标,利用NDBI>0.05 进一步剔除这些误提的地物. 此外,如果分割斑块的像元小于30,则将其剔除以减少小斑块噪声影响,获取最终的水体信息.

2 实验结果与分析

2.1 研究区数据

如图3所示,选取不同地形条件下的LandsatETM+/OLI影像进行细小水体提取实验.

图3 研究区的遥感影像Fig.3 Remote sensing image of study areas

影像1为2001年9月1日获取的Landsat7ETM+影像地表覆盖类型, 以落叶针叶林、 灌木林和城市建筑用地为主,地形起伏较为明显. 因此,研究区1以河北省潘家口水库为中心选取,涵盖滦河流域内的多条细小支流. 影像2为2014年12月13日获取的Landsat8OLI影像,其轨道号为119-42. 该区域属亚热带海洋性季风气候,平均降雨量1 400～2 000mm,是中国雨量最丰富的区域之一. 本文选取两个区域进一步分析细小水体的提取精度,其中研究区2位于福建省罗源县霍口溪流域,大小为1 000px×1 000px,该区域的植被覆盖类型以常绿针叶林为主,地形同样存在较大起伏,流域范围内涵盖黄浦溪、 后溪、 斌溪等支流. 研究区3位于闽江中游支流尤溪河流域,大小为1 000px×1 000px,该区域的地表覆盖类型以常绿针叶林、 城市建筑用地为主,北部地形较为平坦,南部地形较为起伏,流域范围内涵盖青印溪、 清溪、 文江溪等多条支流.

本文以高分辨率的谷歌地球影像为基础,对研究区范围内的细小水体进行手工数字化,生成参考水体数据,以确保实验验证过程中细小水体的准确性.

2.2 MNWI影像分析

利用MNWI对研究区1的MNDWI影像进行线性增强,增强前后分别如图4(a)、 4(c)所示. 分析特征影像可知,在MNDWI影像中,裸地和城区的亮度值与细小河流较为接近,在水体信息提取时易造成混淆； 而在线性增强后的MNWI影像中,细小河流的亮度明显增强,更容易与背景中的城区和裸地区分. 对4(a)、 4(c)中标注的样本区域进行放大,分别如图4(b)、 4(d)所示. 通过对比可以发现,MNDWI中细小水体的亮度差较大,这使得分割提取细小水体时较暗的部分容易出现漏提,而这一现象在MNWI中得到改善,细小水体整体亮度值较高,从而更有利于阈值分割.

图4 比较MNDWI与MNWI中的细小水体Fig.4 Comparison of narrow water in MNDWI and MNWI image

水体信息绘制该区域MNWI影像的ROC曲线[12], 如图5所示. 进一步分析可发现,ROC曲线前半段上升较快,当虚警率仅为5%时,查全率便可达到90%以上,而当查全率到达95%后,增速逐渐减缓,这表明真实的参考细小水体与MNWI影像呈高度正相关性,MNWI有较好的细小水体增强效果.

图5 样本区域MNWI影像的ROC和AUC曲线Fig.5 ROC and AUC curves of MNWI image in the focus area

图6 MNWI影像的特征值统计Fig.6 Statistical feature value of MNWI image

利用图1所选取的地物样本对MNWI影像进行特征统计,结果如图6所示. 与MNDWI影像中各种地物的特征值分布进行对比发现,MNWI影像中的细小水体得到显著增强,其特征值远高于其它地物,因而容易通过阈值分割出细小水体. 由于道路同样呈线状分布,它同样在增强后的MNWI影像中具有较大的特征值,使得难以通过单一阈值的方式将二者准确区分.

分别利用单阈值和双阈值分割提取样本区域内的水体信息,结果如图7所示. 定性分析4种水体分割结果可以发现,以单阈值20分割MNWI特征影像,虽较完整地提取了样本区域内的水体信息,但却误提了影像下侧的一段道路,并且提取结果中存在较严重的斑块噪声； 当单阈值提升至30时,误提的道路信息明显减少,但水体信息的完整度也随之下降； 最后将单阈值提升至40,虽然完全避免了道路的误提,但提取的河流信息却中出现较多的断点； 而双阈值分割则较好地避免了道路信息的干扰,且降噪效果较为明显.

图7 MNWI影像分割结果对比Fig.7 Comparison of MNWI image segmentation results表1 不同方法分割MNWI影像的精度评价Tab.1 Accuracy evaluation of different methods segment MNWI images (%)

此外,本文还基于2.1中的参考水体数据,对样本区的全部区域进行定量精度评价(如表1所示). 观察评价结果可以发现,随着单阈值的升高,细小水体提取结果的用户精度逐步提升,误提取现象不断减少； 但此时制图精度也在逐渐下降,细小水体中出现越来越多的漏提取现象. 然而使用双阈值分割方法可避免此类问题,不仅获取了较高的用户精度,降低了混淆地物的误提,而且制图精度高达99.1%,能够更为完整地提取细小水体. 因此本文使用双阈值分割方法进行细小水体提取.

2.3 细小水体提取结果与分析

为了评价细小水体提取结果,图8分别显示对本文方法和最优阈值分割MNDWI影像法的提取结果. 图中采用不同的颜色对准确提取、 漏提、 误提的细小水体分别进行标注. 结果表明本文方法在个研究区均较完整地提取了细小水体. 研究区1中的6条主要支流提取结果均较为完整； 而漏提水体主要集中在潘家口水库的东侧,漏提河道宽度小于10 m. 这些漏提的河道在Landsat影像中不足一个像元,本文的方法增强细小水体线性后仍难取得较好的效果. 同样,研究区2和3的细小水体整体提取效果较好,但研究区2在斌溪和霍口溪段出现两处明显漏提,主要是因为该区域河道普遍狭窄以及沿岸树阴遮挡较为严重. 由图8提取结果知,最优阈值法的细小水体漏提和间断都很严重,仅有少量支流提取较为完整,其中研究区2的水体漏提最为严重.

此外,为了对细小水体提取结果进行定量评价,本文以研究区的全部区域为评价对象,以制图、 用户精度为评价指标,对本文方法、 最优阈值法、 文献[9]中基于形态学膨胀的后处理方法(对最优阈值分割结果进行后处理)、 文献[7]分割LFE影像的提取方法进行精度评价,结果如表2所示. 其中由于目标水体和噪声道路均以线状形式在提取结果中出现,为使评价结果更具代表性,本文采用细小河流长度作为统计指标参与精度评价.

图8 3个研究区的多光谱影像及其水体提取结果Fig.8 Multi-spectral images and water extraction results of 3 study areas表2 研究区内不同细小水体提取方法的精度评价Tab.2 Accuracy evaluation of different narrow water extraction methods in the study areas (%)

对比4种不同方法的提取精度,可以发现最优阈值分割法的制图精度最低,且该方法在研究区2的制图精度明显低于其他区域,这表明该方法提取的细小水体完整性较差且在河道狭窄处精度下降严重. 文献[9]中的方法虽然提高了制图精度,并减少了细小水体的漏提,但在研究区2内的用户精度却下降近50%. 这可能是因为通过膨胀滤波和图像细化操作实现的断点连接,忽略了地类的真实纹理,难以捕捉细小水体的真实位置,较多的误提使得该方法用户精度较低. 文献[7]中基于LFE影像的提取方法与本文方法都取得了较好的制图精度,较为完整地提取了细小水体. 但本文方法自动化程度更高,在确定初始目标水体分割阈值的前提下可自动提取细小水体,并且误提率更低,具有更好的用户精度.

3 结语

现有的水体指数多对开阔水体提取有效,但细小水体的提取中容易存在漏提现象,影响水体的整体提取精度. 本文以MNDWI的水体指数特征影像为基础,利用形态学白帽变换,构建出一种细小水体指数MNWI,使得细小水体与背景更易区分. 实验结果表明,本文方法可较好地提取细小水体,与最优阈值分割等方法相比,水体提取更为连续,且有着更高的提取精度. 研究表明MNWI对不同宽度细小水体的增强效果有所差别,当宽度大于20 m时(即宽度在Landsat影像中大于2/3像素)有较好的增强效果,但当宽度小于15 m时增强效果显著下降.

尽管本文采用MNDWI构造形态学细小水体指数MNWI,本文的方法同样适用其他常用水体指数,如NDWI等. 本文主要以中等分辨率的Landsat影像为数据源进行方法设计,该方法是否适用于Spot-4或GF-1等中等分辨率影像、 Worldview或Quickbird等高分辨率影像,仍有待今后开展大量的相关验证工作.

[1] 朱长明, 骆剑承, 沈占锋, 等. DEM辅助下的河道细小线性水体自适应迭代提取[J]. 测绘学报, 2013, 42(2): 277-283.

[2] 都金康, 黄永胜, 冯学智, 等. SPOT卫星影像的水体提取方法及分类研究[J]. 遥感学报, 2001, 5(3): 214-219.

[3] 刘炜, 王聪华, 赵尔平,等. 基于面向对象分类的细小河流水体提取方法研究[J]. 农业机械学报, 2014, 45(7): 237-244.

[4] 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.

[5] 徐涵秋. 利用改进的归一化差异水体指数(MNDWI)提取水体信息的研究[J]. 遥感学报, 2005, 9(5): 589-595.

[6] 王刚, 李小曼, 田杰. 几种TM影像的水体自动提取方法比较[J]. 测绘科学, 2008, 33(3): 141-142.

[7] JIANG H, FENG M, ZHU Y,etal. An automated method for extracting rivers and lakes from Landsat imagery[J]. Remote Sensing, 2014, 6(6): 5067-5089.

[8] 李艳华, 丁建丽, 闫人华. 基于国产GF-1遥感影像的山区细小水体提取方法研究[J]. 资源科学, 2015, 37(2): 408-416.

[9] 杨树文, 薛重生, 刘涛, 等. 一种利用TM影像自动提取细小水体的方法[J]. 测绘学报, 2010, 39(6): 611-617.

[10] OTSU N A. Threshold selection method from gray - scale histograms[J]. IEEE Trans on SMC, 1979, 9(1): 62-66.

[11] 查勇, 倪绍祥, 杨山. 一种利用TM图像自动提取城镇用地信息的有效方法[J]. 遥感学报, 2003, 7(1): 37-40.

[12] WU B, HUANG B, FUNG T,etal. Projection of land use change patterns using kernel logistic regression.[J]. Photogrammetric Engineering & Remote Sensing, 2009, 75(8): 971-979.