CN108830876A - 一种水域轮廓提取及面积检测方法 - Google Patents

Abstract

本发明公开一种水域轮廓提取及面积检测方法，步骤为：确定设定区域不同时间段卫星遥感图像；对图像进行锐化处理；批量进行兴趣区裁剪，并标定水域中心点；选取卫星影像进行处理，根据水体特征可计算得到得到每个像素点的NDWT值并获得水体判别阈；判定其水陆分类；若本次识别无陆地，则本次识别所有区域相连环形区以内均为水域；以原有半径为准，增长一个半径单位，若遇重合岸线则停止，非重合岸线部分继续扩展；若某次所有增长区内均无水域出现，则检测完毕；根据之间获得的增长区区域个数，结合增长区设定面积，测定图像中水域的总面积。本发明具有普适性的研究方法，能够大幅度减少水域轮廓提取中人工的部分，使水域轮廓提取更便捷与客观。

Description

一种水域轮廓提取及面积检测方法

技术领域

本发明属于数字图像处理与模式识别技术领域，特别涉及一种水域轮廓提取及面积检测方法。

背景技术

遥感技术是一种借助于卫星、飞机等进行空间探测的技术，具有探测范围广、快速、动态监测、经济客观等一系列优点。我国水域发展方面一度借助遥感技术进行探测。

随着遥感技术的成熟，水域成为了目前遥感领域的一大重点研究对象。水域轮廓能反应水体的多种信息，其体现了自然的发展方向，人类文明程度，同时也对改善附近城市气候有着十分重大的意义。在多重因素的影响下，如自然环境变化和人为改造，湖泊的形状在时空上一直发生在着变化，因此，从现有的光学遥感图像中自动提取水域轮廓线，是目前更新区域地图与地形数据的重要举措。同时，检测并计算水域面积的技术，可以为做水域未来面积预测提供必须的条件，人们通过面积预测可以预知洪涝灾害从而避免损失且可以更好地进行水资源保护。故提取水域的轮廓并测算面积对人类发展与开辟新领域都有着极大的影响。

从遥感图像上自动地提取地物要素一直是摄影测量与遥感领域的重要研究方向之一。作为一类重要的地物要素，水体的识别和提取，依然是通过作业人员目视判读、人工采集。虽然通过目视判读来识别影像中的水体能够满足一定的精度要求，但是在众多外界因素的影响下，极易造成错误判断，对最终结果产生恶性影响，耗时耗力，劳动强度大、生产效率低。

综合以上需求以及现有条件的不足，目前迫切的需要结合高科技信息化的手段，提出一种全自动不需人工干预的一种水域轮廓提取及面积检测方法。

发明内容

发明目的：针对现有技术中存在的问题，本发明提供一种具有普适性的研究方法，能够大幅度减少水域轮廓提取中人工的部分，使水域轮廓提取更便捷与客观的水域轮廓提取及面积检测方法。

技术方案：为解决上述技术问题，本发明提供一种水域轮廓提取及面积检测方法，包括如下步骤：

(1)确定设定区域不同时间段卫星遥感图像；

(2)对图像进行锐化处理；

(3)批量进行兴趣区裁剪，并标定水域中心点；

(4)选取卫星影像进行处理，根据水体特征可计算得到得到每个像素点的NDWT值，根据NDWT值可计算获得水体判别阈；

(5)通过将NDWT值与水体判别阈值进行比对，判定其水陆分类；

(6)若本次识别无陆地，则本次识别所有区域相连环形区，判别环形区以内均为水域；以原有半径为准，增长一个半径单位，若遇重合岸线则停止，非重合岸线部分继续扩展；

(7)若某次所有增长区内均无水域出现，则判定检测完毕，此时将其增长区区域内轮廓相接得到水域轮廓；

(8)根据之间获得的增长区区域个数，结合增长区设定面积，测定图像中水域的总面积。

进一步的，所述步骤(1)中所述卫星遥感图像为由Landsat系列卫星采集的光学影像。

进一步的，所述步骤(4)中NDWI值的表达式为：

NDWI＝(K(G)-K(R))/(K(G)+K(R))

K(G)是基于绿波段的归一化比值指数；

K(R)是基于近红外波段的归一化比值指数。

进一步的，所述步骤(4)中根据NDWT值计算获得水体判别阈的具体步骤如下：根据每一增长区的NDWI值，计算得到水体判别阈值；其中增长区区域的确定以半径单位做模拟长和模拟宽，以此长宽来确定增长区区域，圈定覆盖的若干相接的增长区后，测算出增长区区域内水体NDWI值，并计算得到水体判别阈值；

所述半径单位为以水域中心点为基点，以水域中心点到最近的额岸线间距离1％的像素点个数作为半径单位，围绕水域中心点做圆。

进一步的，所述步骤(5)中通过将NDWT值与水体判别阈值进行比对，判定其水陆分类的具体步骤如下：分别将所述图像每一次增加的水体增长区中NDWI值与水体识别阈值做对比处理：

若增长区的NDWI值较所述水体判别阈值大，则可判定其为水域；

若增长区的NDWI值较所述水体判别阈值小，则可判定其为陆地；

若本次识别无陆地，则本次识别所有区域相连环形区，判别环形区以内均为水域；以原有半径为准，增长一个半径单位；此时将原有水域范围设定为初值，将增长一个半径单位后的范围设定为复值。

进一步的，所述步骤(2)中对图像进行锐化处理的具体步骤如下：选用累积分布函数为对图像中像素个数多的灰度级进行展宽，而对图像中像素个数少的灰度进行压缩，使像素值均匀分布，完成对256个灰度级中各灰度数目的统计；在统计完成后计算灰度分布密度和直方图分布，最后对累计分布取整并对灰度值进行映射；其中，映射公式为：

其中，式中n表示像素的总和，而G为可能的灰度级总数。

与现有技术相比，本发明的优点在于：

本发明通过采用神经网络算法完成水域岸线轮廓提取与水域面积测算的方法，包括：确定设定区域的卫星遥感图像，将遥感图像锐化，计算得到遥感图像中感兴趣区域并进行像素识别。针对每一区域计算NDWI值，并由NDWI值计算得到阈值，用于判定其为水体或陆地。在图中人工识别并选取水域中心，利用每个像素点的水体指数与所述水体阈值进行对比，可鉴定其是否为水域，再由水域中心点为原点，以原有半径加一个半径单位向像素值扩展方向作圆，推算一周，继续比较NDWI值。以一个半径长度作为岸线停止的一个标准度量单位，重合岸线处停止比较，实际重合岸线的长度应大于设定岸线停止长度的整数倍。重合部分外多于此长度的部分，继续向前推进。最终确定岸线轮廓，提取轮廓，计算面积。以往人工标记得方法误差较大，应用上述方法提取得到得水域岸线结果具有更高的客观性，且可以全自动进行，不需人工干预。

附图说明

图1为本发明的流程图；

图2为实施例中洪泽湖处于丰水期时的卫星遥感影像示意图；

图3为实施例中多布扎湖处于丰水期时的卫星遥感影像示意图；

图4为实施例中洪泽湖丰水期时在湖泊中心标定水域中心点示意图；

图5为实施例中多布扎湖丰水期时在湖泊中心标定水域中心点的示意图；

图6为实施例中以洪泽湖丰水期为例的一个半径长度作为岸线停止标准度量单位的示意图；

图7为实施例中以洪泽湖丰水期为例的一个半径单位示意图；

图8为实施例中以洪泽湖丰水期为例的一个增长区示意图；

图9为实施例中洪泽湖丰水期时湖泊轮廓某一初值示意图；

图10为实施例中洪泽湖丰水期时湖泊轮廓在上图基础上增长一个半径单位的复值示意图；

图11为实施例中洪泽湖丰水期遇重合岸线停止时，非重合岸线(水域)部分继续向外扩展所得示意图；

图12为实施例中洪泽湖枯水期遇重合岸线停止时，非重合岸线(水域)部分继续向外扩展所得示意图；

图13为实施例中洪泽湖湖泊丰水期洪泽湖及附近水域轮廓提取示意图；

图14为实施例中洪泽湖湖泊丰水期洪泽湖主体轮廓提取示意图；

图15为实施例中洪泽湖湖泊枯水期洪泽湖及附近水域轮廓提取示意图；

图16为实施例中洪泽湖湖泊枯水期洪泽湖主体轮廓提取示意图。

具体实施方式

下面结合附图和具体实施方式，进一步阐明本发明。

一种水域轮廓提取及面积检测方法，参见图1所示，该方法包括：

S101、获取设定地理区域的卫星遥感图像；

具体的，本发明实施例用于在从卫星拍摄获得的遥感图像中识别水域与陆地，提取水域轮廓。卫星遥感图像具有较高的宏观性和实时性，具体在实施本发明实施例技术方案时，可根据需求选择特定日期，特定水域或特定地理区域的卫星遥感图像。

例如，就丰水期选择洪泽湖2000年5月的卫星影像(如图2所示)以及多布扎湖2000年7月的卫星影像(如图3所示)；

S102、利用直方图均衡化方法对图像进行锐化处理；

由于部分原始图像的像素占有很多的灰度级，而且分布均匀，具有较高的对比度和多变的灰度色调，导致图像局部对比度不高。为了改善图像的视觉效果，需对图像进行增强处理：即通过一定的手段对原图附加一些信息或变换数据，有选择地突出图像中感兴趣的特征或者抑制图像中某些不需要的特征，使图像与视觉相应特性相匹配。

具体的，选用累积分布函数为对图像中像素个数多的灰度级进行展宽，而对图像中像素个数少的灰度进行压缩，使像素值均匀分布，完成对256个灰度级中各灰度数目的统计。在统计完成后计算灰度分布密度和直方图分布，最后对累计分布取整并对灰度值进行映射，从而拓展像元取值的动态范围，提高了对比度和灰度色调的变化，使图像更加清晰。

其中，映射公式为：

其中，式中n表示像素的总和，而G为可能的灰度级总数。

S103、批量进行兴趣区裁剪，并标定水域中心点(如图4、图5所示)；

具体的，在卫星遥感影像中，覆盖面积较大区域较广，若目标水域过小，为使所观测图像更清晰更具说服力，应采用批量裁剪方式对目标进行裁剪，以裁剪后图像作为研究对象进行后续步骤。

例如，洪泽湖在条带号120，行编号37的卫星遥感影像上占据面积过小，不便于观测，则沿水域预留合适部分的陆地进行边缘分割得到批量裁剪图像。同样，多布扎湖位于条带号139，行编号40的卫星遥感影像上，其上包含多布扎湖同时包含多庆错湖，不便于对多布扎水域的识别与观测，故进行裁剪。

需要说明的是，由于研究对象不同，湖泊占据卫星影像大小不同，是否需要S103裁剪步骤依具体情况决定；对于湖泊面积过大的情况，也可采取拼接不同遥感影像的方式进行处理。

S104、选取卫星影像进行处理，根据水体特征可计算得到每个像素点的NDWI值，根据NDWI值可计算获得水体判别阈值；

具体的，从卫星遥感图像中可以获得每个像素点的基本数值，包括基于绿波段和基于近红外波段的归一化比值指数。

归一化水体指数计算公式：NDWI＝(K(G)-K(R))/(K(G)+K(R))

K(G)是基于绿波段的归一化比值指数。

K(R)是基于近红外波段的归一化比值指数。

S105、通过将NDWI值与水体判别阈值进行对比，判定其水陆分类；

具体的，在由于在遥感图像中水体和陆地的归一化指数有着明显的区别，一般情况下，水体归一化指数大于陆地的归一化指数，故可由此判别所检测区域是否为水域。可设定一个阈值用于判定检测区域是否为水体，通过讲像素点的NDWI值与其进行对比，确定此区域的水陆类型。

例如，以洪泽湖与多布扎湖为例，提取阈值如下：

湖畔归一化水体指数阈值

S106、若本次识别无陆地，则本次识别所有区域相连环形区，判别环形区以内均为水域。以原有半径为准，增长一个半径单位，若遇重合岸线则停止，非重合岸线部分继续向外扩展并重复(5)和(6)步骤。

具体的，首先设定湖泊中心点为识别起始点与参考点，从湖泊中心点向外扩展，根据计算得到的卫星遥感图像中各个像素点的归一化水体指数，计算水体阈值。在上次研究区域中所有NDWI值范围内取最小值作为水体归一化指数识别阈值A，并选择本次研究区域内出现频次最高的水体归一化指数，并命名为水体归一化指数识别阈值B，将阈值A与阈值B作比较，选择较小项作为最终的水体归一化指数阈值，用于后续步骤的比较。

进一步的，本次研究区域应为以中心点为圆心的圆形剔除掉上次研究区域，形成环形兴趣区。若本次识别未检测到陆地，则应继续扩大识别范围，在原有半径的基础上，增长一个半径单位，增长后区域剔除掉所有已研究过区域后，形成一个新的环形兴趣区。重复以上步骤，继续测算NDWI值并进行判定，若仍未发现陆地则继续扩大，检测到陆地则停止。

需要说明的是，在检测到某一处出现陆地时，由于陆地并不呈现规则的环形，故并不代表其他方向检测也将接近岸线。将出现陆地的一侧重合岸线确定出后，可以停止对重合岸线半径扩展，但非重合岸线部分应继续进行检测，直到出现线陆地为止。

S107、若某次所有增长区内均无水域出现，则判定检测完毕，此时将其增长区区域内轮廓相接得到水域轮廓。

具体的，仍利用水陆阈值进行判定，若均无水域，则认定已寻找到水域轮廓，此时将前期所有重合岸线用平滑的线段相连接，得到最终提取结果。

S108、根据之前获得的增长区区域个数，结合增长区设定面积，测定图像中水域的总面积。

具体的，在前期应对每一步的增长区进行数量统计，在检测到所有陆地时对增长区个数进行求和，最终得到增长区的总数。由于前期设定对于同一工程内每个增长区的面积必须相同，故最终根据单个增长区面积与总增长区个数可以便捷地得到水域地面积。

对于结果的精度可以利用Fleiss Kappa系数进行评价：Fleiss Kappa系数是检验实验标注结果数据一致性比较重要的参数，其可以用于多份标注结果的一致性检测。通过计算Pj和Pi可以对计算结果进行较为客观的评价：

对于Pj的计算，为每一列结果相乘除以任务总数，对应于Kappa系数中每个分类的随机一致概率，然后由于Fleiss Kappa没有一个参照的标注值，因此这里计算Pe理论一致性的时候，将每个Pj进行平方求和。

而对于Pi的计算，是对每一个标注任务进行实际一致性的计算，也就是要计算：在所有的排列中(这里用组合也可以，但是包含顺序比较方便计算)，具有一致性的排列有多少个，具体计算代码可借助python完成。

本发明提出一种利用神经网络算法完成水域岸线轮廓提取与水域面积测算的方法，包括：根据研究面积的条带号以及行编号确定设定区域的卫星遥感图像，将遥感图像锐化，计算得到遥感图像中感兴趣区域并进行像素识别。针对每一区域计算NDWI值，并由NDWI值计算得到阈值，用于判定其为水体或陆地。在图中人工识别并选取水域中心，利用每个像素点的水体指数与所述水体阈值进行对比，可鉴定其是否为水域，再由水域中心点为原点，以原有半径加一个半径单位向像素值扩展方向作圆，推算一周，继续比较NDWI值。

将一个半径单位作为岸线停止(如图6所示)的一个标准度量单位，重合岸线处停止比较，实际重合岸线的长度应大于一个半径单位的整数倍。重合部分外多于此长度的部分，继续向前推进。最终确定岸线轮廓，提取轮廓，计算面积。以往人工标记得方法误差较大，应用上述方法提取得到得水域岸线结果具有更高的客观性，且可以全自动进行，不需人工干预。

根据每一增长区的NDWI值，计算得到水体判别阈值；

此时将原有水域范围设定为初值，将增长一个半径单位后的范围设定为复值。

具体的，分别将所述图像每一次增加的水体增长区中NDWI值与水体识别阈值做对比处理，若：

若增长区的NDWI值较所述水体判别阈值大，则可判定其为水域；

若增长区的NDWI值较所述水体判别阈值小，则可判定其为陆地；

若本次识别无陆地，则本次识别所有区域相连环形区，判别环形区以内均为水域。以原有半径为准，增长一个半径单位。此时将原有水域范围设定为初值，将增长一个半径单位后的范围设定为复值(如图7所示)。

以水域中心点为基点，以固定像素点为半径单位，围绕水域中心点做圆，包括：

水域中心点由操作人员人为标定，便于识别兴趣区。但水域中心点的标定仅在将陆地错标为水域时会对结果产生恶性影响，若中心点标记不准确，即存在较大偏差，只要中心点仍标于水域，由于设定了则半径单位,故不会对结果造成影响。

具体的，所述圈定覆盖的若干相接的增长区，包括：

沿不同方向将圆周上长度分为若干相连的半径单位，以半径单位做模拟长和模拟宽，以此长宽确定一个增长区区域(如图8所示)。

具体的，在每次判定后识别新增加的增长区个数，并对增长区进行计数；

在每次判定后识别新增加的增长区个数，并对增长区进行计数；

要求每个增长区的面积恒等，便于在全部水域岸线识别结束后依据增长区的个数对水域面积进行计算。

以一个半径单位作为岸线停止的标准度量单位，包括：

每次扩大半径更新检测范围时都对水体增长区域进行计数(如图9、图10所示)图9和图10表示其所检测面积逐步增长的过程，最终检测完毕时根据水体增长区的单个面积及增长区的个数进行加和计算，得到水域的总面积。

具体的，根据水域轮廓形态等，结合如前所述半径单位，若检测中增长区出现陆地，图11、12表示检测中遇到岸线的情况，并认定其为重合岸线，则在重合岸线处停止扩展与比较，将以重合岸线作为此区域水陆分界线，并认定以重合岸线为始点向背离水域中心点向外扩展的部分视作陆地。实际重合岸线的长度应大于一个半径单位的整数倍。重合部分外多于此长度的部分，继续向前推进。

所述水域边界鉴定完毕时，需进行水域轮廓的连接和提取：

具体的，沿所识别到的所有增长区内轮廓进行逐次连接，所有方向选取最后一次扩展的增长区作为用于确定边界的增长区。根据所有连接线勾勒得到水域矢量边界，作为图像水域轮廓提取的最终结果，并自动生成缩略图，图13表示丰水期时对洪泽湖及周边进行提取获得的轮廓；图14表示丰水期时提取的洪泽湖主体轮廓；图15表示枯水期时对洪泽湖及周边进行提取获得的轮廓；图16表示枯水期时提取的洪泽湖主体轮廓。

根据增长区，进行面积计算：每次扩大半径更新检测范围时都对水体增长区域进行计数，最终检测完毕时根据水体增长区的单个面积及增长区的个数进行加和计算，得到水域的总面积。

对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。

Claims (6)

1.一种水域轮廓提取及面积检测方法，其特征在于，包括如下步骤：

(1)确定设定区域不同时间段卫星遥感图像；

(2)对图像进行锐化处理；

(3)批量进行兴趣区裁剪，并标定水域中心点；

(4)选取卫星影像进行处理，根据水体特征可计算得到得到每个像素点的NDWT值，根据NDWT值可计算获得水体判别阈；

(5)通过将NDWT值与水体判别阈值进行比对，判定其水陆分类；

(6)若本次识别无陆地，则本次识别所有区域相连环形区，判别环形区以内均为水域；以原有半径为准，增长一个半径单位，若遇重合岸线则停止，非重合岸线部分继续扩展；

(7)若某次所有增长区内均无水域出现，则判定检测完毕，此时将其增长区区域内轮廓相接得到水域轮廓；

(8)根据之间获得的增长区区域个数，结合增长区设定面积，测定图像中水域的总面积。

2.根据权利要求1所述的一种水域轮廓提取及面积检测方法，其特征在于，所述步骤(1)中所述卫星遥感图像为由Landsat系列卫星采集的光学影像。

3.根据权利要求1所述的一种水域轮廓提取及面积检测方法，其特征在于，所述步骤(4)中NDWI值的表达式为：

NDWI＝(K(G)-K(R))/(K(G)+K(R))

K(G)是基于绿波段的归一化比值指数；

K(R)是基于近红外波段的归一化比值指数。

4.根据权利要求3所述的一种水域轮廓提取及面积检测方法，其特征在于，所述步骤(4)中根据NDWT值计算获得水体判别阈的具体步骤如下：根据每一增长区的NDWI值，计算得到水体判别阈值；其中增长区区域的确定以半径单位做模拟长和模拟宽，以此长宽来确定增长区区域，圈定覆盖的若干相接的增长区后，测算出增长区区域内水体NDWI值，并计算得到水体判别阈值；

所述半径单位为以水域中心点为基点，以水域中心点到最近的额岸线间距离1％的像素点个数作为半径单位，围绕水域中心点做圆。

5.根据权利要求1所述的一种水域轮廓提取及面积检测方法，其特征在于，所述步骤(5)中通过将NDWT值与水体判别阈值进行比对，判定其水陆分类的具体步骤如下：分别将所述图像每一次增加的水体增长区中NDWI值与水体识别阈值做对比处理：

若增长区的NDWI值较所述水体判别阈值大，则可判定其为水域；

若增长区的NDWI值较所述水体判别阈值小，则可判定其为陆地；

若本次识别无陆地，则本次识别所有区域相连环形区，判别环形区以内均为水域；以原有半径为准，增长一个半径单位；此时将原有水域范围设定为初值，将增长一个半径单位后的范围设定为复值。

6.根据权利要求1所述的一种水域轮廓提取及面积检测方法，其特征在于，所述步骤(2)中对图像进行锐化处理的具体步骤如下：选用累积分布函数为对图像中像素个数多的灰度级进行展宽，而对图像中像素个数少的灰度进行压缩，使像素值均匀分布，完成对256个灰度级中各灰度数目的统计；在统计完成后计算灰度分布密度和直方图分布，最后对累计分布取整并对灰度值进行映射；其中，映射公式为：

其中，式中n表示像素的总和，而G为可能的灰度级总数。