一种卫星遥感影像的城市水体提取方法与流程

本发明属于水体遥感影像技术领域，主要涉及一种卫星遥感影像的城市水体提取方法。

背景技术：

城市是人类社会高度发展的体现，城市水体作为城市生态系统中重要的因素，在维持城市生态系统稳定性上具有十分重要的作用。城市水体的改变将会对生活产生巨大的变化，可能会引发一些灾害，如城市水源滞留、水资源短缺，甚至引发一些和人类健康生活相关的疾病。因此，了解和掌握城市水体分布以及水域面积的变化已经成为了人们日益关注的焦点。

近年来，随着遥感技术的应用与发展，遥感影像在自然资源调查、动态监测、自然地表水源规划等方面发挥着着越来越重要的角色，利用遥感技术进行地表的监测也得到了越来越多的科研工作者的关注；遥感影像可以以一个不同的视角去观察地球表面的地物，实时的监测地表的变化。在水体提取技术中，通过遥感数据及时准确地获取城市水体信息成为目前主流的水体提取方式。到目前为止，诸多学者提出了大量遥感影像的水体提取方法，但大部分都是基于中低分辨率遥感影像，由于影像分辨率较低，面积较小的水体未能有效的提取；尤其是对于城市地区水体的提取，城市地区水体面积大小参差不齐，存在诸多小面积人工湖泊和细小的河流。因此，在城市水体提取中应更多的采用高分辨率遥感影像，以资源3号卫星为例，资源3号卫星遥感影像有着5.8m分辨率和5.2km的宽幅，它为城市水体提取提供了理想的多光谱影像数据，资源3号具体参数如下表1：

表1资源三号卫星影像参数

水体是遥感影像中很常见的地物类别，也是重要的基础地理信息之一，其动态信息的快速获取，对水资源调查、水利规划、环境监测与保护等事业都有着十分明显的实用价值和科学意义。对此，诸多学者很早就对此展开了研究，提出了许多有效的水体信息自动提取模型。大致可以分为4类：(a)单波段和多波段阈值分割法(single-band or multiple-band threshold method)，(b)水体指数法(Water indices)，(c)线性分解模型(linear un-mixing model)，(d)监督与非监督分类方法(supervised or unsupervised classification method)。除此之外，还有一些其他的方法，如：基于数字高程模型的水体提取技术，基于微波遥感(Microwave Remote Sensing)影像的水体提取技术，面向对象(Object Oriented)技术的水体提取技术等。但这些方法并不常用，总得来说，水体指数法由于其模型简单、方便，且精度较高，在实际中最为常用。

然而随着遥感影像分辨率的提高，大部分高分辨率遥感影像(如；WorldView-2,IKONOS,RapidEye,and资源3号)没有像Landsat TM/ETM+/OLI这么多可以利用的波段用于水体的提取，因此，MNDWI(改进的归一化差异水体指数法)和AWEI(自动水体指数)将无法使用，因为大部分的高分辨率遥感影像只有4个波段(蓝色、绿色、红色、近红外波段)，缺少MNDWI/AWEI计算所需的短红外波段(SWIR)。因此，在使用NDWI(归一化差异水体指数)对高分辨率影像进行水体提取的时候就将会产生一些问题，如阴影无法去除的问题，尤其是城市地区高层建筑物的阴影，在高分辨率影像上表现尤为突出。在高分辨率遥感影像上，城市区域的高层建筑物阴影与水体很难区分，尽管目前有相关学者对这一方面进行了研究，如基于面向对象的技术，通过计算高层建筑物阴影区域的纹理特性来对高层建筑物阴影进行检测；虽然可以达到预期效果，但由于纹理的描述与计算相对复杂且耗时较长，所以从计算时间上考虑该方法并不是一个理想的阴影检测方法。也有基于SVM特征训练进行的阴影检测，以期达到去除水体检测中高层建筑物阴影的影响。SVM是一种分类精度较高的方法，但SVM训练需要花去较多的时间，尤其是当训练样本数目较多且样本特征向量维度较高时。若通过采用高分辨率遥感影像阴影检测方法(morphological shadow index,MSI)与NDWI相结合的方式对WorldView-2高分辨率影像进行水体提取，以期提高水体检测精度；尽管该方法原理简单，但由于该方法以NDWI方法作为水体提取的基础，其检测精度并不会很高，尤其是周围植被茂密的细小面积水域，水体的光谱特性将受到严重的污染，水体光谱特征表现极其不稳定，同时城区水体具有悬浮泥沙含量高、水体富营养化严重、受各种污染物污染较大等特点，使得城区水体与自然界中未受污染的水体表现不同的光学特性。

技术实现要素：

针对如何从高分辨率卫星遥感影像上，进行城市区域的水体提取问题，特别是有效的区分高层建筑物阴影与水体和提高水体提取的精度问题，本发明提出了一种卫星遥感影像的城市水体提取方法(Automatic urban water extraction method,AUWEM)。

本发明的具体技术方案如下：

一种卫星遥感影像的城市水体提取方法，包括以下步骤：

步骤1：遥感影像数据的预处理，即对遥感影像数据进行正射校正和大气校正；

步骤2：预处理后的遥感影像数据，包括蓝色band1数据、绿色band2数据、红色band3数据和近红外band4数据，选取预处理后的遥感影像数据中的蓝色band1数据代替归一化差异水体指数NDWI的计算公式中的绿色band2数据，获得新的归一化差异水体指数NNDWI1，新的归一化差异水体指数NNDWI1的计算公式为：

此计算公式即为NNDWI1指数模型，利用此模型通过阈值分割即得到NNDWI1的阈值分割结果，也即为NNDWI1水体提取结果；

步骤3：对预处理后的遥感影像数据中包括的四个波段数据，即蓝色band1数据、绿色band2数据、红色band3数据和近红外band4数据进行PCA变换，并将PCA变换后的第一主成分分量Component1替代归一化差异水体指数NDWI的计算公式中的绿色band2数据，获得另一个新的归一化差异水体指数NNDWI2，即：

其中，Component1表示PCA变换的第一主成分分量，此计算公式即为NNDWI2指数模型，利用此模型通过阈值分割即得到NNDWI2的阈值分割结果，也即为NNDWI2的水体提取结果；

步骤4：将步骤2得到的NNDWI1的阈值分割结果和步骤3中得到的NNDWI2的阈值分割结果进行叠加，将得到的结果定义为新的归一化差异水体指数NNDWI的阈值分割结果，即NNDWI1的阈值分割结果与NNDWI2的阈值分割结果进行叠加，其计算公式为：

NNDWI＝(segmentation_NNDWI1)∪(segmentation_NNDWI2)

式中segmentation_NNDWI1表示NNDWI1的阈值分割结果，segmentation_NNDWI2表示NNDWI2的阈值分割结果，此计算公式即为NNDWI指数模型，利用此指数模型得到NNDWI的水体提取结果；

步骤5：对预处理后的遥感影像数据中的近红外band4数据进行阈值分割，得到近红外band4数据的阈值分割结果；

步骤6：对NNDWI的水体提取结果中的大面积水体对象和小面积对象进行分割，NNDWI的水体提取结果中，像素个数大于设定阈值的为大面积水体对象，像素个数小于等于设定阈值的为小面积对象；

步骤7：对步骤6中得到的小面积对象进行数学形态学膨胀处理，得到膨胀后的小面积对象，将步骤5得到的近红外band4数据的阈值分割结果作为约束条件，即采用膨胀后的小面积对象和近红外band4数据的阈值分割结果求交集的方式对膨胀后的小面积对象进行约束，约束的数学表达式为：

component2＝(dilate_component)∩(segmentation_band4)

式中，dilate_component表示膨胀后的小面积对象，segmentation_band4表示近红外band4数据的阈值分割结果，component2表示约束后的小面积对象；

步骤8：对步骤7得到的约束后的小面积对象进行阴影检测与去除，得到小面积水体对象；

阴影检测与去除，是指对每个小面积对象中的每个像元进行波谱关系的描述，并判断该像元是否满足阴影像元的条件，记录并统计每个小面积对象中阴影像元的个数，当一个小面积对象中阴影像元所占比例大于阈值T时，把该小面积对象判定为建筑物阴影对象，小于等于阈值T时小面积对象则判定为小面积水体对象，阴影像元所占比例即为小面积对象中阴影像元的个数与该小面积对象中总像元个数的比值，区分小面积对象中小面积水体对象和阴影对象的函数表达式为：

式中，n表示某一小面积对象中总像元个数，m为该小面积对象中阴影像元的个数；

阴影像元的条件，是指满足阴影像元的波谱大小关系，即满足以下三个不等式条件：

步骤9：将步骤6中得到的大面积水体对象与步骤8中得到的小面积水体对象进行叠加，即将步骤6中得到的大面积水体对象和步骤8中得到的小面积水体对象求并集，得到卫星遥感影像的城市水体提取结果。

本发明的有益效果如下：

(1)新的归一化差异水体指数NNDWI1和NNDWI2可以有效的对水体进行初始提取，提高后续水体的提取精度。

(2)在阴影对象获取方面，为了更加准确的获取阴影对象，对小面积对象进行膨胀处理；同时为了限制在真实地表阴影区域，对膨胀的对象采用近红外band4数据的阈值分割结果进行约束。

(3)为避免阴影对象检测过程中时间消耗，对小面积对象采用基于波谱特性的描述，减少传统中采用纹理特征进行描述的时间消耗，提高方法的计算效率。

(4)AUWEM方法的分类精度要高于NDWI的分类精度和最大似然法的分类精度。在5个实验区的实验结果中，AUWEM平均Kappa系数约为93.0％，NDWI的平均Kappa系数约为86.2％，最大似然法分类精度介于两者之间，平均Kappa系数约为88.6％；同时AUWEM漏检与虚警总错误率也要少于NDWI的分类结果和最大似然法的分类结果，AUWEM平均漏检与虚警总错误率约为11.9％，最大似然法平均漏检与虚警总错误率约为18.2％，NDWI平均漏检与虚警总错误率约为22.1％％。除此之外，AUWEM具有更高的水体边缘检测精度。

附图说明

图1是本发明具体实施方式中的NNDWI的水体提取结果，其中图(a)是从资源3号卫星获得的遥感影像，图(b)是NNDWI1水体提取结果，图(c)是NNDWI2水体提取结果，图(d)NNDWI水体提取结果；

图2是本发明具体实施方式中的膨胀约束过程示意图；

图3(a)～(e)是本发明具体实施方式中的阴影区域不同波谱特性曲线；

图4是本发明具体实施方式中的卫星遥感影像的城市水体提取方法流程图；

图5是本发明具体实施方式中的不同方法在5个实验区的实验结果；

图6(a)～(f)是本发明具体实施方式中的不同方法在5个实验区的6个指标直方图；

图7是本发明具体实施方式中的水体边缘精度评估中待评估区域获取示意图；

图8(a)～(c)是本发明具体实施方式中的不同方法在5个实验区水体边缘检测精度比较。

具体实施方式

下面结合附图和具体实施方式对本发明进行详细说明。

一种卫星遥感影像的城市水体提取方法，以资源3号卫星遥感影像数据为例，包括以下步骤：

步骤1：遥感影像数据的预处理，即对遥感影像数据进行正射校正和大气校正；对实验区域的影像利用RPC+30DEM进行无控制点正射校正，采用Feyisa G L et al.的FLAASH(Fast Line-of-Sight Atmospheric Analysis of Spectral Hypercubes)大气校正模型进行影像的大气校正，均在ENVI5.2软件中完成。其中，资源3号FLAASH大气校正需要的定标系数可通过网站http://www.cresda.com/CN/Downloads/dbcs/index.shtml进行下载，光谱响应函数可通过网站http://www.cresda.com/CN/Downloads/gpxyhs/index.shtml进行下载。

步骤2：预处理后的遥感影像数据，包括蓝色band1数据、绿色band2数据、红色band3数据和近红外band4数据，选取预处理后的遥感影像数据中的蓝色band1数据代替归一化差异水体指数NDWI(Normalized Difference Water Index)的计算公式中的绿色band2数据，获得新的归一化差异水体指数NNDWI1(New Normalized Difference Water Index 1)，新的归一化差异水体指数NNDWI1的计算公式为：

此表达式即为NNDWI1指数模型。利用此模型在编程软件Matlab 2014a中通过阈值分割得到NNDWI1阈值分割结果，也即为NNDWI1水体提取结果。

步骤3：对预处理后的遥感影像数据中包括的四个波段数据，即蓝色band1数据、绿色band2数据、红色band3数据和近红外band4数据进行PCA变换，并将PCA变换后的第一主成分分量Component1替代归一化差异水体指数NDWI的计算公式中的绿色band2数据，获得另一个新的归一化差异水体指数NNDWI2(New Normalized Difference Water Index 2)，即：

其中，Component1表示PCA变换的第一主成分分量。

此表达式即为NNDWI2指数模型。利用此模型在编程软件Matlab 2014a中通过阈值分割得到NNDWI2阈值分割结果，也即为NNDWI2水体提取结果。

步骤4：将步骤2得到的NNDWI1的阈值分割结果和步骤3中得到的NNDWI2的阈值分割结果进行叠加，将得到的结果定义为新的归一化差异水体指数NNDWI(New Normalized Difference Water Index)的阈值分割结果，即NNDWI1的阈值分割结果与NNDWI2的阈值分割结果进行叠加，其计算公式为：

NNDWI＝(segmentation_NNDWI1)∪(segmentation_NNDWI2)

式中segmentation_NNDWI1表示NNDWI1指数影像的阈值分割结果，segmentation_NNDWI2表示NNDWI2指数影像的阈值分割结果。此计算公式即为NNDWI指数模型，利用此模型在编程软件Matlab 2014a中得到NNDWI水体提取结果。

通过实验分析，实验中NNDWI1的阈值分割结果所采用的阈值与NNDWI2的阈值分割结果所采用的阈值都设为0，可获得较为理想的实验结果。NNDWI的水体提取结果综合了NNDWI1的水体提取结果和NNDWI2的水体提取结果，避免了单一指数对水体的漏检现象发生。如图1(a)～(d)所示。NNDWI1对含泥沙高的水体检测效果较为明显，而NNDWI2对受植被光谱信息干扰严重的水体较为敏感，因此，实际水体提取中结合NNDWI1的水体提取结果和NNDWI2的水体提取结果，整合生成新的水体提取结果，提高后续水体提取精度。

步骤5：对预处理后的遥感影像数据中的近红外band4数据进行阈值分割，得到近红外band4数据的阈值分割结果。

通过分析大量NNDWI水体提取结果的影像资料发现，除部分城区面积较小的人工池塘与湖泊之外，在一般情况下，城市高层建筑物阴影对象的面积一般要小于水体对象的面积。所以在实际的分析中，只对面积相对较小的地物进行检测，因为在这部分面积相对较小的对象中包括几乎所有的阴影对象，同时也包括小面积水体对象。

步骤6：对NNDWI的水体提取结果中的大面积水体对象和小面积对象进行分割，NNDWI的水体提取结果中，像素个数大于设定阈值的为大面积水体对象，像素个数小于等于设定阈值的为小面积对象；

小面积对象获取模型，可表示为：

其中，t为设定的阈值，其值取最大阴影对象的像素个数，component表示NNDWI水体指数提取结果中的离散对象，area(component)表示NNDWI水体指数提取结果中的离散对象的面积，area(component)>t为大面积水体对象，area(component)≤t为小面积水体或者阴影对象。

步骤7：对步骤6中得到的小面积对象进行数学形态学膨胀处理，得到膨胀后的小面积对象，将步骤5得到的近红外band4数据的阈值分割结果作为约束条件，即采用膨胀后的小面积对象和近红外band4数据的阈值分割结果求交集的方式对膨胀后的小面积对象进行约束，约束的数学表达式为：

component2＝(dilate_component)∩(segmentation_band4)

式中，dilate_component表示膨胀后的小面积对象，segmentation_band4表示近红外band4数据的阈值分割结果，component2表示约束后的小面积对象；

该膨胀后的小面积对象更加完整的包括了小面积水体对象与阴影的阴影像素，对膨胀后的小面积对象进行约束也保持了小面积对象的膨胀操作限制在地表真实阴影区域的范围内。

以步骤6得到的小面积对象为建筑物阴影对象为例的膨胀约束过程示意图如图2所示，从资源3号假彩色影像中，分别获取建筑物阴影对象和Band4影像数据阈值分割结果，将得到的建筑物阴影对象进行数学形态学膨胀处理，获得建筑物阴影对象的膨胀结果，将建筑物阴影对象的膨胀结果与Band4影像数据阈值分割结果进行求交集处理，即得到在Band4影像数据阈值分割结果约束下的膨胀后的建筑物阴影对象。

步骤8：对步骤7得到的约束后的小面积对象进行阴影检测与去除，得到小面积水体对象；

在NNDWI水体提取结果中，基本上只包含小面积水体对象与阴影，所以只需要对小面积水体对象与阴影的特征进行研究与分析，找到适合区分小面积水体对象与阴影的特征。在实验中发现，尽管纹理特征可以很好的描述小面积水体对象与阴影，但由于地物的纹理特征，如：灰度共生矩阵，其计算相对复杂、耗时较长，不太适于用于小面积水体对象与阴影的区分，所以在实验中采用地物的光谱特征描述小面积水体对象和阴影的像素，以此作为依据区分小面积水体对象与阴影。

通过分析大量的水体与阴影的波谱特性曲线，得到水体像元的波谱关系满足不等式：

band2>band4

而阴影的像元波谱曲线较为复杂，实验中分析总结出了以下5种波谱特性曲线，如图3(a)～(e)所示。

根据上述波谱曲线各波段对应的大小关系，实验中总结出阴影像元的波谱关系满足以下三个不等式条件：

由此可知，阴影像元的条件，是指阴影像元的波谱关系满足以上三个不等式条件。

阴影检测与去除，是指对每个小面积对象中的每个像元进行波谱关系的描述，并判断该像元是否满足阴影像元的条件，记录并统计每个小面积对象中阴影像元的个数，当一个小面积对象中阴影像元所占比例大于阈值T时，把该小面积对象判定为建筑物阴影对象，小于等于阈值T时小面积对象则判定为小面积水体对象，阴影像元所占比例即为小面积对象中阴影像元的个数与该小面积对象中总像元个数的比值，区分小面积对象中小面积水体对象和阴影对象的函数表达式为：

式中，n表示某一小面积对象中总像元个数，m为该小面积对象中阴影像元的个数。阈值T是通过实验统计得来的，对资源3号遥感影像数据的阴影像元进行统计发现，当T取0.5时可以很好的区分水体与阴影对象。

步骤9：将步骤6中得到的大面积水体对象与步骤8中得到的小面积水体对象进行叠加，即将步骤6中得到的大面积水体对象和步骤8中得到的小面积水体对象求并集，得到卫星遥感影像的城市水体提取结果。

一种卫星遥感影像的城市水体提取方法的总体流程图如图4所示。

为验证方法的有效性，分别采用NDWI方法水体提取结果和最大似然法(MaxLike)水体提取结果进行对比实验。在中国境内选取了5处不同地区且具有不同周边环境的影像用于实验，它们包括湖泊和河流，分别位于中国地区的北京市、武汉市、苏州市、广州市，其中武汉市选用了两个不同覆盖区域的影像。资源3号影像详细信息描述如表2所示，实验中详细的参数的设置情况如表3所示，其中band4的数值首先按下式归化到0-255取值范围，然后再选取阈值进行分割的。

表2实验数据详细信息

表3不同实验地区阈值设定(其中T1，T2，T3，分别为NNDWI1、NNDWI2、band4的分割阈值)

为便于不同方法分类结果的目视判读与分析，对正确分类的水体像元赋予灰色，正确分类的非水体像元赋予黑色，错误分类的像元赋予白色，实验结果如图5所示。从图5的分类结果目视判读可以发现，本发明提出的AUWEM的水体提取分类精度要好于NDWI的水体提取分类精度和最大似然法的水体提取分类精度。AUWEM对水体边缘混合像元可以很好的分类(结合Beijing、Wuhan_1和Wuhan_2地区的水体分类结果)、对细小的池塘水体检测性能好于NDWI和最大似然法(结合Suzhou地区水体分类结果)、对房屋建筑物的阴影可以很好的去除(结合Suzhou和Wuhan_2地区水体分类结果)。

不同实验区域三种方法水体提取分类精度比较统计结果如表4所示，从表4的统计结果中发现，AUWEM水体提取分类精度要高于NDWI和最大似然法。AUWEM在这5个实验区域的分类精度是最高，平均Kappa系数达93.0％；而NDWI的分类精度最低，平均Kappa系数约为86.2％；最大似然法分类精度介于两者之间，平均Kappa系数约为88.6％。

表4不同实验区域三种方法的精度统计

为了更加详细的评估AUWEM对水体的提取精度，采用生产者精度、用户精度、Kappa系数、漏检率、虚警率、总错误率6个指标来描述方法的水体提取精度。三种不同方法在5个实验地区的6个指标直方图如图6(a)～(f)所示，从各个指标的直方图可以发现，AUWEM的水体提取分类精度要高于NDWI的水体提取分类精度和最大似然法的水体提取分类精度。AUWEM在水体提取的虚警率方面，除在Suzhou地区达到9.1％左右，其他的实验区域都是低于5％；在水体漏检率方面，5个地区漏检率都要明显少于NDWI和最大似然法。当水体提取的虚警率与水体漏检率都低时，总错误率自然也就是最低。从图6(a)～(f)中可以看出AUWEM总错误率最低，其次是最大似然法，总错误率最高的是NDWI，其中AUWEM平均总错误率约为11.9％，最大似然法平均总错误率约为18.2％，NDWI平均总错误率约为22.1％％。

对于水体提取分类生产者精度，AUWEM水体提取分类生产者精度最高，平均约为91.6％；最大似然法水体提取分类生产者精度次之，平均约为84.8％；NDWI的水体提取分类生产者精度最低，平均约为81.6％。在水体提取用户精度方面，最大似然法的水体提取用户精度要大于的AUWEM和NDWI，平均精度用户高达96.6％；AUWEM方法次之，平均精度用户约为96.4％；最差的仍然属于NDWI方法，平均精度用户约为95.7％。

为了更加客观评价三种方法提取的水体的边缘检测精度，设计以下方法来进行精度评价，方法具体实施描述如下：

首先利用参考影像采用Canny算子获取三种方法提取的水体的边界线；

对获取的边界线进行数学形态学的膨胀处理，建立以边界线为中心半径为4个像元的缓冲区域；

然后对缓冲区域的像元进行判断，假设缓冲区域总像元值为N，正确分类的像元数目的N_R，漏检像元数目为No，虚警数目为Nc，则：

其中，A+Eo+Ec＝100％。A表示边缘像元正确划分类别的比例，在这里称它为边缘像元正确分类精度；Eo表示边缘像元漏检的比例，称它为边缘像元漏检率；Ec表示边缘像元虚警的比例，称它为边缘像元虚警率。将经过数学形态学膨胀处理的水体的边界线与从资源3号卫星获取的遥感影像进行叠加，得到边缘精度评估中待评估区域。边缘精度评估中待评估区域获取示意图如图7所示。

根据上述方法统计了实验区域的水体边缘检测精度，分别统计了边缘像元虚警率(Commission Error)、边缘像元漏检率(Omission Error)和边缘像元正确分类精度(Accuracy of edge detection)，并对实验结果进行比较，如图8(a)～(c)所示。从图8的比较结果可以很清晰的发现：AUWEM方法边缘像元正确分类精度要好于NDWI和最大似然法。其中AUWEM方法水体边缘像元正确分类精度最高达93.7691％(Guangzhou地区)，最小精度为79.5798％(Wuhan_2地区)；NDWI水体边缘像元正确分类精度最高达84.0917％(Suzhou地区)，最小精度为69.8310％(Beijing地区)；最大似然法水体边缘像元正确分类精度最高达85.8149％(Guangzhou地区)，最小精度为69.7974％(Wuhan_2地区)。