面向电子地图更新的水体自动提取方法研究

何正国，毛海亚，钟家晖，樊惠萍，苗立新，周淑芳，刘 烽，陈 婷

(1. 广州市城市规划自动化中心，广东 广州 510030；2. 二十一世纪空间技术应用股份有限公司，北京 100096)

水体提取的难点主要由光谱多样性、形态多样性和季节变化因素引起。水体包括面状、线状等不同的形态，同时随年份、季节波动，造成水体在不同区域和不同时相影像上形成迥异的特征，这都对特征选取、信息提取和质量控制等提出了较高要求，使得陆表水体信息提取技术更加复杂[1-4]。本文选取广州市作为研究区，基于高分辨率的北京二号遥感数据和面向对象的技术流程，提出一套基于全局最优阈值法的水体自动提取方法，旨在提高电子地图更新的自动化程度，缩短更新周期。

1 研究区与数据

1.1 研究区概况

广州市是我国的国家中心城市、广东省省会城市和珠三角城市群中心城市。国家“一带一路”建设已经进入纵深发展阶段，粤港澳大湾区城市群发展亦写入了2017年国务院《政府工作报告》，广州作为国家中心城市以及粤港澳大湾区核心门户城市之一，城市化进程越来越快。其空间数据的更新应与之相符，广州市电子地图数据作为典型的空间数据之一，应与城市发展同步，实现高频度的地图更新。

广州市目前正在开展“四标四实”工作，每天有大量的网格员采集数据，因此有必要对传统人工更新电子地图的方式进行革新，研究电子地图要素信息的自动提取方法。选取广州市南部作为研究区，见图1。该区域内水体类型丰富，既包括了天然的大型河流、湖泊，又包括了人工修筑的水库、坑塘，水体周边的地物环境也包含了建设用地、植被、裸地等不同类型，可以充分检验水体自动提取方法在不同地物构成下的适用性和有效性。

1.2 数据处理

由于城市建设步伐的加快，城市基础地理信息的变化也越来越快，公众对电子地图的现势性要求越来越高，常规的电子地图更新时效已经无法满足用户的需求。遥感技术对地物信息的快速、大面积获取能力，使其在电子地图的数据获取领域受到越来越多的重视。研究中采用了亚米级的北京二号卫星影像数据，用于水体的自动提取实验。

图1 广州研究区位置示意图

北京二号遥感卫星星座于2015年7月11日在印度的斯里赫里戈达岛(Sriharikota)卫星发射场成功发射。星座由3颗光学遥感卫星组成，均搭载有多光谱和全色成像仪，对全球任一地点的重访周期为1～2 d，在时间尺度上为水体的高密度动态监测提供了可能。每颗卫星在轨提供幅宽约24 km、0.8 m分辨率全色(450～650 nm)和3.2 m分辨率蓝(440～510nm)、绿(510～590 nm)、红(600～670 nm)、近红外(760～910 nm)多光谱图像，波段设置满足遥感水体信息提取的一般要求。

使用1景北京二号影像，成像时间为2018年1月12日，影像预处理主要包括系统几何校正、正射校正、波段融合等。结合地形图，通过选取地面同名控制点的方法对影像进行几何校正，纠正后各点误差在一个像元内。通过控制点采集及其优化、区域网平差模型构建及其解算，完成全色和多光谱影像的正射校正处理。采用Gram-Schmidt变换方法对校正后的全色、多光谱数据做融合处理，融合后的影像既保留了多光谱信息，又在原有基础上提高了其空间分辨率，得到亚米级的影像数据。

2 技术方法

2.1 面向对象分割

研究采用面向对象的技术流程，通过多尺度分割(multi-resolution segmentation)算法，使影像对象成为特征信息的载体[5]。该算法根据像元间的异质性，通过定义一定的阈值将影像分割成相对同质的对象，并以对象作为基本的分类单元。其综合利用了影像的光谱特征、几何特征，避免了对象破碎，提高了影像的分类精度[6-7]。

多尺度分割算法是一种区域生长算法。它基于影像对象的特征，针对不同土地覆盖类型的尺度效应，解决产生斑块与真实地物的边界拟合问题。分割对象的平均异质性最小化与像素的平均异质性最小化是影像分割的标准。当异质性大于等于给定的尺度时，区域就会停止生长，区域内像元将合并形成对象。异质性的计算公式为：

f=wcolor·[∑cwc·δc]+(1-wcolor)·

(1)

式中：wcolor为光谱异质性权重；wc为数据层权重；δc为数据层上光谱标准差；wcpt为紧致度异质性；l为对象长度(像素个数表示)；b为目标多边形最短边(像素个数表示)；n为像元个数。

对于多尺度分割，分割尺度的确定至关重要，是后续分类精度的重要保证。分割尺度越大，误分割图斑的数量会减少，但误分割的总面积会增大。为综合考虑两种因素，对误分割图斑个数和误分割图斑面积进行归一化处理，建立误分割指数模型，以期得到最佳分割尺度使得误分割指数达到最小值。误分割指数计算公式如式(2)所示。

ESI=n2+r2

(2)

式中：n为归一化后误分割图斑数；r为归一化后误分割图斑面积；ESI为误分割指数。

2.2 特征提取规则集

基于表观要素中的光谱特征和纹理特征等，构建水体异于其他地物类型的特征集。相较于传统的分类技术而言，面向对象的信息提取方法是基于对象的操作，能够充分利用提取对象的空间、质地和光谱特征，集合临近像元进行分割和分类，最终实现分类结果的矢量输出。

(1)波段均值(Mean_X)。是指每个对象内第X波段所有像元灰度值的平均。其中，红波段(Mean_R)与蓝波段(Mean_B)的灰度值均值，可以分别剔除大部分建设用地与蓝色屋顶的工厂板房。此外，部分水体表面被水生植被所覆盖，归一化植被指数(NDVI)取值范围与地表植被很接近，基于NDVI的提取混分严重，采用近红波段均值(Mean_NIR)可以扩大这部分水体与植被的差异，减少错提。亮度(Brightness)也属于波段均值的一种，是通过所有影像层计算出来的，该特征可以区分水体与周边地物，特别是与建筑物阴影的差异。

(2)归一化水体指数。水体因对入射能量(太阳光)具有强吸收性，所以在大部分遥感传感器的波长范围内，总体上呈现较弱的反射。Mcfeeters提出了归一化差异水体指数(NDWI)，旨在利用比值运算快速提取水体信息[8-9]。由于水体的反射从可见光到中红外波段逐渐减弱，在近红外和中红外波长范围内吸收性最强，几乎无反射，因此基于可见光波段和近红外波段的反差构成的NDWI可以突出影像中的水体信息。此外，由于植被在近红外波段的反射率一般最强，因此采用绿光波段与近红外波段的比值可以最大程度地抑制植被的信息，从而达到突出水体信息的目的。针对广州城市水体进行分析，多数水体对象的NDWI取值在0.5以上，可以比较明显地区别于植被与建设用地。

(3)归一化植被指数。比值型指数创建的基本原理就是在多光谱波段内，寻找出所要研究地类的最强反射波段和最弱反射波段，将强者置于分子，弱者置于分母。通过比值运算，进一步扩大二者的差距，使感兴趣的地物在所生成的指数影像上得到最大的亮度增强，而其他背景地物则受到普遍的抑制，从而达到突出感兴趣地物的目的。基于此原理，归一化植被指数是反映植物生长状况及植被覆盖度的最佳指示因子。研究区内水体对象的NDVI指数普遍处于-0.3以下，区别于大多数植被。

(4)灰度共生矩阵(GLCM)。灰度共生矩阵原理是：如果有n×n个灰度大小的影像，则该影像能够转化为n×n阶的灰度共生矩阵，在影像平面(x,y)中的某处(i,j)对应一个向量D=(Δx,Δy)。对于随意一个这样的向量和影像，可以计算由D向量偏移的成对出现的亮度值之间的联合概率密度P(i,j)。灰度共生矩阵的纹理信息都包含在移动方向为θ=0°、θ=45°、θ=90°和θ=135°的亮度值空间依赖矩阵中。能够用于提取纹理特征的纹理测度有很多统计值，而其中最常用的几个二阶概率纹理特征统计值主要是：均值(Mean)、方差(Variance)、同质性(Homogeneity)、对比度(Contrast)、非相似度(Dissimilarity)、角二阶矩(Angular Second Moment, ASM),以及熵(Entropy)和相关性(Correlation)，其中的每一个GLCM纹理测度都可以单独使用。

值得注意的是，在人工确定阈值的过程中，各单波段均值、亮度值及NDWI与NDVI的设定要依区域位置和时相的情况而定，不同区域、时相的影像差别较大。特别是在设定NDWI取值范围时，要根据研究区实际情况结合参考数据反复调试提取阈值的大小。这些虽然会达到理想的提取效果，但耗费较多人力与时间，在业务化实施中并非是最佳选择。

2.3 全局最优阈值算法

水体包括面状、线状等不同的形态，其提取难点主要由光谱多样性、形态多样性和季节变化因素引起。在采用水体指数、植被指数等进行水体信息提取过程中，指标阈值的确定是非常关键的。为了提高生产效率，同时针对本项目开展的业务化应用需求，采取了自动计算阈值的方法进行水体提取。

水体通常具有较大面积，呈现面状分布的特点，且内部特征相对均质，适用于根据水体指数、植被指数等设定阈值对影像进行前景与背景的划分。采用全局最优阈值的方法，其中心思想是基于迭代的方式将灰度图像用两个或两个以上正态分布的概率密度函数做近似表示，每次都取其中最显著的波峰来划分区域，然后依据各个区域的平均值选择合适阈值，重复过程直到该值收敛。方法实现如下。

遍历图像中各点读取像素，获得分割图像的最小灰度值T1以及最大灰度值T2，设定初始阈值T(0)，T(0)=(T1+T2)/2；根据阈值T(0)进行图像分割，将图像分为前景和背景两个区域。再次对分割后的图像进行遍历，计算出此时的前景T1和背景的平均灰度值T2，再次平均求出阈值T(1)。用gray(i,j)表示待处理图像像素(i,j)点的灰度值，N1表示图像前景像素的个数，N2表示图像背景像素的个数，计算如式(3)～式(5)。

(3)

(4)

T(1)=(T1+T2)/2

(5)

重复上述运算，直到T(k)=T(k+1)，迭代完成，此时用T(k)作为最终阈值对图像进行分割。

全局最优阈值法的水体提取技术首先根据影像的NDWI指数将其分为前景和背景值，然后利用全局最优阈值法自动获取NDWI、NDVI、Brightness指数的阈值，分别提取水体和剔除混分的植被、建设用地等要素。同时，考虑到在高分辨率遥感影像中，建筑物阴影是较为明显的干扰因素，且往往被错提为水体，因而在规则集中增加了纹理特征的灰度共生矩阵(GLCM)，用于进一步剔除建筑物阴影等干扰。经过对比实验，特别采用了纹理同质性指标(GLCM Homogeneity)作为对影像分布平滑性的测度，其计算公式如式(6)所示，匀质区域具有较高的特征值，反之则较低。

hom=∑∑p(i,j)/[1+(i-j)2]

(6)

3 结果与讨论

3.1 水体对象的多尺度分割

主要基于面向对象的方法进行水体信息提取，该方法根据像元间的异质性，通过定义一定的阈值将影像分割成相对同质的对象，并以对象作为基本的分类单元。其综合利用了影像的光谱特征、几何特征，避免了对象破碎，提高了影像的分类精度[10]。

先对研究区内的北京二号影像进行多尺度分割，分割后的特征值反映了每个对象范围内所有像元的平均值。当分割尺度过大时，水体边界与地面实体相差较大，边界不准确，对象内其他地物“杂质”过多、同质性差，为后期信息提取带来困难。而分割尺度过小时，虽然地物边界清晰，对象同质性较好，能够准确地反映地物空间分布状况，但会使分割过程及后期提取效率低下，混分的细碎斑块不易剔除，也难以形成规整的水体边界，不利于提取成果通过GIS方法等进行分类后处理(见图2)。经过反复对比，确定分割尺度为50(形状参数0.1，紧致度参数0.5)，以尽量将光谱、纹理特征接近的像元分割确定为同一对象，同时避免斑块过于细碎。经过多尺度分割后的影像对象产生了新的特征信息，相较于单个像元无疑含有更多的信息量。

本研究基于多尺度分割方法充分利用了多个尺度的信息，对分割结果进行约束优化，是目前一种较好的信息提取方法，应用也较广。多尺度分割是一种自下而上，采用异质性最小的区域合并策略的影像分割方法，分割时基于局部最佳相互适配原则，将每个像元逐步合并成较大的多边形对象。通过合理设置分割参数，使分割后的对象尽可能体现出研究区范围内不同构成要素间的差异，以便后续对样本特征参数的分析与选择。对水体进行多尺度分割时的参数选择应结合影像空间分辨率与地表覆盖的具体情况而定，分割结果既能满足同类要素的均一性，又能体现不同要素间的差异，并达到分割精度与效率的统一。

3.2 基于规则集的提取结果

获取了2018年覆盖广州研究区的北京二号融合影像，包含可见光、近红外共4个波段，空间分辨率0.8 m。在该区域内开展了水体(主要为河流、湖泊、坑塘等)的自动提取实验，确认方法及技术流程的可行性。共采用了NDWI、NDVI、Brightness、GLCM Homogeneity (band 1)与GLCM Homogeneity (band 2)等特征，均基于全局最优阈值法进行计算，确定各提取规则的特征值，见表1。首先充分利用NDWI指数在水体提取中的优势，基于迭代的方式将NDWI灰度图像的特征直方图用正态分布的概率密度函数做近似表示，取其中最显著的波峰将图像划分为前景和背景两个区域。其中，前景即为初步提取的水体，背景为其他地物类型。随后，对前景中的水体做进一步优化，不断剔除错提的其他地物。需要再次对分割后的影像进行遍历，分别依据NDWI、NDVI、Brightness、GLCM Homogeneity等每一类特征，计算出前景和背景的平均值，并由此获取各特征的合适阈值，重复过程直到该值收敛。最后获取的特征阈值，即是经过不断迭代优化后的取值。

注：(a)为北京二号融合影像，(b)、(c)、(d)所对应的分割尺度分别为30、50、70图2 不同分割尺度下的影像对象变化

表1 全局最优阈值法确定的提取规则阈值

水体的初步提取结果如图3所示，研究区范围内除了大型河流外，还存在较多的小型水体与坑塘，且零星分布。提取结果中的水体边界相对规整，后处理工作量小。除主要河流外，坑塘等小型水体往往受到季节、气象及人为因素的影响，年际变化较大，这也是南方水体普遍存在的明显特征。因而，使用自动提取方法可以较为准确、快速地获取其分布情况，有利于减少对此类小型水体的漏提。主要的后处理操作为排查与归并小图斑，剔除伪变化，同时也可以应对数据冗余。

图3 广州研究区水体自动提取结果

3.3 结果后处理与精度评价

在完成初步的信息提取后，还需要进一步开展结果的后处理工作，并主要针对于水体中面积较小的图斑对象，以有效减少错提和漏提这两类现象，见图4。首先将水体图斑进行合并，使空间位置上接壤的对象合并为一个较大的对象。这样一方面减小了后续操作的图斑数据量，提高运行效率；另一方面也便于通过设定面积阈值，剔除掉面积较小的孤立图斑。基于北京二号影像空间分辨率及研究区内水体分布的实际情况，确定面积阈值为4000像素，将小于该值的水体对象赋类为其他，可以对错提的阴影等要素进行剔除，见图4中(a)、(d)、(g)。同时，由于研究区内水运发达，水面船只数量繁多，见图4中(b)、(e)、(h)，农田水利设施广泛分布,见图4中(c)、(f)、(i)，考虑到水面船只、水利设施等对水体提取的干扰，将此类面积小于17 000像素且完全被水体图斑所包含的对象，也重新归类于水体，保证了水体图斑的完整，减少了漏提。

对水体自动提取的精度与效率进行分析，在研究区内生成300个随机点，通过人工判定方式对自动提取结果开展精度评价。经过混淆矩阵的计算，水体的生产者精度为64%，用户精度则达到了96%，Kappa系数为0.73。整个研究区面积达到了近300 km2，水体总面积约为98 km2，水体图斑总数超过了2800个，全流程的对象分割、信息提取可以在1 h之内完成。其中,除多尺度分割与GLCM纹理特征计算较为耗时以外，基于全局最优阈值法的各特征值获取均在3 min内完成，具有较高的计算效率。从最终提取水体对象的数量、面积及准确度而言，自动提取方法可以极大地提高工作效率，有利于快速获取水体要素分布情况，而人工参与只需要安排在后续的成果精细化修订环节即可。

注：(a)、(b)、(c)为北京二号融合影像，(d)、(e)、(f)为水体初步提取结果，(g)、(h)、(i)为后处理结果图4 针对错堤和漏堤的水体提取结果后处理

值得注意的是，基于高空间分辨率影像开展水体提取，在城市范围内往往存在建筑物等地物阴影的干扰问题。在本研究区内，部分高层建筑物、植被及山体的阴影会对相邻地物造成明显遮挡。如图4(a)所示，阴影遮挡下的地物光谱信息受到抑制，特别是植被、裸地等亮度、纹理等特征更接近于水体，导致在多尺度分割时就难以识别其准确边界，因而产生错提。如果在此基础上开展进一步的全覆盖地物类型提取，则会使误差继续累加。针对这一问题，可以考虑在研究区内提前圈定尽可能大的水体矢量范围，作为水体提取的目标区域，以便排除建成区内易混淆的高大建筑物阴影等。规则集运行中只提取目标区域内的水体，这样可以比较有效地规避建筑物及植被阴影的影响，减少水体错提。

4 结 论

本文以广州市部分水体为研究对象，基于北京二号高分辨率遥感影像，采用面向对象的技术流程，通过全局最优阈值算法实现了水体的自动化信息提取。研究表明，基于全自动的阈值确定方法，在无人工参与情况下，可以在较大面积的区域内实现水体要素的快速提取，应用效果较好。研究区内水体类型复杂，表观特征差异较大，针对电子地图更新的数据需求，通过面向对象的方法，提取的水体图斑边界规整，能够满足图4的精度要求。水体的提取作为电子地图数据更新的第一步，也为后续植被、建筑物、道路等信息提取与更新提供了重要参考，和传统的人工矢量化相比，有效提高了工作效率。提取结果的总体精度可以达到96%，Kappa系数为0.73。