CN107240113A

CN107240113A - 一种基于空间剖面线的半自动水体范围提取方法

Info

Publication number: CN107240113A
Application number: CN201710370602.XA
Authority: CN
Inventors: 陈晓玲; 常好雪; 张展; 蔡晓斌; 王玉玺
Original assignee: Wuhan University WHU
Current assignee: Wuhan University WHU
Priority date: 2017-05-23
Filing date: 2017-05-23
Publication date: 2017-10-10
Anticipated expiration: 2037-05-23
Also published as: CN107240113B

Abstract

本发明公开了一种基于空间剖面线的半自动水体范围提取方法，对选取的遥感影像做辐射定标、几何校正、影像的拼接裁剪等必要的数据预处理工作，对处理后的影像做波段运算，计算NDWI；然后选取n条空间剖面线，绘制空间剖面线矢量文件；将空间剖面线矢量文件转为与影像空间分辨率相同的栅格文件；结合NDWI影像和空间剖面线栅格文件，通过空间运算，获取剖面线栅格文件的属性表；通过程序运算获取水体提取最佳分割阈值Ｔ；根据阈值，通过阈值分割的方法提取水体边界范围。本发明提取精度可以达到93.7%，可以真实地反映水体范围，在水体范围动态提取的研究中具有较好的应用前景。

Description

一种基于空间剖面线的半自动水体范围提取方法

技术领域

本发明属于空间信息技术在水利工程学领域应用技术领域，涉及一种半自动的水体范围提取方法，具体涉及一种基于空间剖面线的半自动水体范围提取方法。

背景技术

地表水是地球表层生态系统重要的组成部分，具有调节区域气候、维持区域生态平衡和保护生物多样性等生态功能。另外，作为重要的国土资源，地表水具有水力发电、农业灌溉、水产养殖和提供航运等社会经济功能。受人类活动和全球气候的影响，地表水资源如湖泊、湿地等的面积发生剧烈变化，对陆地生态系统造成了严重影响。然而由于地域、时间、经费等问题使得水资源普查与统计工作难以开展，我国地表水面积的动态分析缺乏长时序大范围的监测资料。遥感具有时效性强、覆盖范围广等观测优点，目前已成为研究地表水资源形态变化的有效手段。

遥感技术的应用，使得实时获取水体范围成为可能，可以为汛期的防洪调度、海啸等灾害损失评估以及抗洪抢险等提供数据支持，以便政府等相关部门做出及时准确的决策。此外，遥感提取的水体范围数据可以为水量平衡、库容研究、悬浮泥沙浓度反演、叶绿素浓度反演等提供基础数据。

水体范围提取的方法，国内外学者已经做了大量的研究工作。其中最常用的是指数模型法，最常用的指数为归一化水体指数NDWI(Normalized Difference Water Index)。NDWI的计算公式为：

(ρ_Green-ρ_Red)/(ρ_Green+ρ_Red)；

其中，ρ_Green和ρ_Red分别为绿光和红光波段的反射率。

指数模型法提取水体范围具有简单快捷的优点，但是提取精度的高低取决于阈值的选取。根据阈值是否固定，可将现有方法分为固定阈值法和动态阈值法；根据阈值的确定方式，可将现有研究分为人工目视阈值选取法、半自动阈值选取法和自动阈值选取法。针对水体提取的这一关键性问题，国内外学者做了大量讨论(文献1-8)，有学者通过辐射传输模型研究发现指数模型提取水体范围时，不同影像阈值差距较大，因此固定阈值不适合做长时序多景影像的水体提取。人工目视阈值选取法，通常需要进行多次实验以确定阈值，费时费力，且主观性较强，没有统一客观的标准，不同人甚至同一人在不同的工作状态下选取的阈值都会有差异。自动阈值选取法具有较高的应用价值，然而其缺陷在于，不经过样本选取，缺乏先验知识而使用一种完全基于像元值的统计运算，不能充分考虑到实际应用中地物的复杂性和差异性。半自动化阈值选取法由人和计算机共同干预，常见的做法是，通过人工选取的感兴趣区域的训练样本来确定最佳阈值。

综合以上列举的文献资料，已有不少学者针对水体范围提取的阈值确定方法做了一系列研究。对比不同方法的缺陷和限制，在当前的技术背景下，半自动化阈值选取法具有较好的应用价值。

[文献1]Ji L,Zhang L,Wylie B.2009.Analysis of dynamic thresholds forthe normalized difference water index[J].Photogrammetric Engineering&RemoteSensing,75(11):1307-1317.

[文献2]Liu Y,Song P,PengJ,et al.2012.A physical explanation of thevariation in threshold for delineating terrestrial water surfaces from multi-temporal images:effects of radiometric correction[J].International journal ofremote sensing,33(18):5862-5875.

[文献3]Karsli F,Guneroglu A,Dihkan M.2011.Spatio-temporal shorelinechanges along the southern Black Sea coastal zone[J].Journal of AppliedRemote Sensing,5(1):053545-053545-13.

[文献4]Wang JD,Sheng YW,Tong TSD.2014.Monitoring decadal lakedynamics across the Yangtze Basin downstream of Three Gorges Dam[J].RemoteSensing of Environment,152:251–269.

[文献5]吴川,张玉龙,许秀贞,张全发.2013.基于Landset TM/ETM和HJ-1A/B影像的丹江口水库水域变化监测研究[J].长江流域资源与环境,09:1207-1213.

[文献6]张超,韩琳,陈亮.2010.柴达木盆地湖泊水面变化遥感监测分析[J].人民黄河,12:54-55+58.

[文献7]刘东,李艳.2012.基于遥感技术的鄱阳湖面积库容估算[J].遥感信息,02:57.

[文献8]孙爱民,冯钟葵,葛小青,罗宇,李山山,闫芳,冯旭祥.2015.利用长时间序列Landsat分析博斯腾湖面积变化[J].中国图象图形学报,08:1122-1132.

发明内容

为了解决上述问题，本发明从地物光谱的差异性入手，提出了一种基于空间剖面线的半自动水体范围提取方法；为有效精确实时的获取水体范围提供一种新思路，以便为汛期的防洪调度、海啸等灾害损失的评估以及抗洪抢险等提供数据支持，辅助政府等相关部门做出及时准确的决策。

本发明所采用的技术方案是：一种基于空间剖面线的半自动水体范围提取方法，其特征在于，包括以下步骤：

步骤1：预处理选取的遥感影像；

步骤2：对处理后的遥感影像做波段运算，计算NDWI影像；

步骤3：在经过预处理的遥感影像上，绘制空间剖面线矢量文件；

步骤4：将空间剖面线矢量文件转换为与遥感影像空间分辨率相同的空间剖面线栅格文件；

步骤5：结合步骤2中获得的NDWI影像和步骤4中获得的空间剖面线栅格文件，通过空间运算，获取剖面线栅格文件的属性表；

步骤6：获取水体范围提取最佳分割阈值T；

步骤7：根据最佳分割阈值T，通过阈值分割的方法，提取水体范围。

相对于现有技术，本发明：由人和计算机共同干预，结合先验知识，人工选取样本，充分考虑了实际应用中地物的复杂性和差异性，利用计算机的计算功能快速准确的确定最佳阈值，既保证了结果的客观性和正确性，也提高了工作效率，节省了人力物力财力。

附图说明

图1为本发明实施例的流程图；

图2为本发明实施例的空间剖面线布设示例图；

图3为本发明实施例的空间剖面线NDWI曲线示例图；

图4为本发明实施例的半自动化阈值法与手动阈值法得到的水面积结果的一致性分析示意图；

图5为本发明实施例的半自动化阈值法与手动阈值法得到的水面积结果的同步分析示意图。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合附图及实施例对本发明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不用于限定本发明。

请见图1，本发明提供的一种基于空间剖面线的半自动水体范围提取方法，包括以下步骤：

步骤1，预处理选取的遥感影像；

对选取的遥感影像做辐射定标、几何校正、影像的拼接裁剪等必要的数据预处理工作，对处理后的影像做波段运算，计算NDWI影像；

辐射定标是将传感器记录的电压或数字量化值转换成具有辐射意义的卫星入瞳处亮度值。影像几何校正利用地面控制点和几何校正数学模型来矫正非系统因素产生的误差，由于校正过程中会将坐标系统赋予图像数据，所以此过程包括了地理编码。几何校正有模型：仿射变换(RST)、多项式和局部三角网(Delaunay Triangulation)等。本发明实施例采用2014年GF-1 WFV数据作为实验数据，辐射定标参数来源于中国资源卫星应用中心网站，使用的v1.0版本为2015-10-13更新的定标参数。选取Landsat 8 OLI L1T级几何精纠正数据产品的全色波段(空间分辨率15m)作为参考影像对辐射定标后的GF-1 WFV数据做几何校正。控制点的选取是几何纠正的关键。在本研究中，尽量选择等位置固定不变且易人工准确辨识的点作为控制点。控制点的布设应尽量保证均匀，且不少于20个。配准精度控制在1个像素以内。

指数模型法是水体提取最简单有效的方法，通过计算相应的指数，可以放大影像上水体和非水体的差异。目前使用最广泛的是NDWI。水体的反射率较低，通常低于10％，远低于其他地物。水体在蓝绿波段有较强的反射，在其他可见光波段吸收较强。水体在近红外波段反射率几乎为零，然而植被和裸土在近红外波段有强烈反射。NDWI利用这一差异来增强水体信号。NDWI是最常用的水体提取参数，在水体提取方面更有优势，于是本实例选取NDWI作为模型参数。对上述经过预处理的GF-1号WFV数据做波段运算计算NDWI。

本发明提供的实施例所选区域为鄱阳湖流域的36个大型湖库，所用的遥感影像为2014年的GF-1 WFV数据共计629景，分辨率为16米。对经过辐射校正和几何校正后的影像进行裁剪，将36个湖库依次以矩形裁出，对于鄱阳湖等范围较大，一景影像不足以完整覆盖的湖库，选取拍摄时间相同或相近的多景影像进行拼接。

步骤2，在经过预处理的遥感影像上，绘制空间剖面线矢量文件；

空间剖面线布设示例如附图2所示，在经过预处理的遥感影像上，人工选取空间剖面线，绘制空间剖面线矢量文件。为了获取精确阈值，保证水体范围提取精度，并兼顾计算效率与剖面线的空间布局，针对面积大于1km²，水体长宽比小于10的水体，每幅遥感影像布设不少于8条空间剖面线，布设位置应位于水体常年平均水边线附近，剖面线长度约为水体平均宽度。除此之外，剖面线的规划要遵循以下几个原则：

(1)空间剖面线必须同时覆盖在目标水体和陆地上，不得全部覆盖在单一地物上，可垂直于水体走向布设；

(2)空间剖面线的布设应充分考虑地表类型的复杂性和差异性，尽量均匀布设并覆盖多种地物。

(3)空间剖面线条数不宜太少以获取尽可能小的阈值范围，从而逼近最佳阈值。

按照上述原则，在经过步骤1处理的2014年36个湖库的GF1WFV影像上，针对每幅影像规划空间剖面线。具体实施方法为，根据湖库的实际情况，如考虑每个湖库的空间形态和水体平均宽度，参考空间剖面线的规划原则，在每幅影像上均匀布设8条矢量线，代表该幅影像的8条空间剖面线。

步骤3，数据格式转换，获取空间剖面线栅格文件。

将空间剖面线矢量文件转为与影像空间分辨率相同的栅格文件。

具体实施时将步骤2中以NDWI影像为底图编辑空间剖面线矢量文件通过数据格式转换，批量转换为与GF-1 WFV数据相同空间分辨率(16米)的栅格数据。

步骤4，将空间剖面线矢量文件转换为与遥感影像空间分辨率相同的空间剖面线栅格文件；

结合步骤2中获得的NDWI影像和步骤4中获得的空间剖面线栅格文件，通过空间运算，获取剖面线栅格文件的属性表，即剖面线上各点对应的NDWI像元值，这些像元值可以绘制成NDWI曲线，每一条空间剖面线都对应一条连续的像元值曲线。以其中一条为例，曲线如附图3所示。

具体实施时将步骤3中获取的剖面线栅格数据与NDWI影像叠加，通过空间运算提取剖面线上NDWI的栅格值，每一条空间剖面线对应一条NDWI曲线。由于空间剖面线规划时是同时覆盖水体和陆地覆盖类型，且充分考虑了地表覆盖类的复杂性和差异性，且水体和陆地覆盖类型的NDWI值差异明显(水体为正值，植被为负值)，因此在剖面线从水体穿入陆地或从陆地穿入水体时，其NDWI值会发生跳变，因此获取的每条NDWI曲线都存在一个或者两个跳变点。依据这一原理可以从NDWI曲线上获取对应水体边界的阈值。

步骤5，通过运算获取水体范围提取最佳分割阈值T。

通过运算可以获取曲线上斜率变化最快的两点，记为点A、B，其像素值，分别记作a1、b1(a1＞b1)，则(b1，a1)为通过该条空间剖面线所得到的阈值范围，记作阈值t1∈(b1，a1)。同样的运算方法，对于同一幅影像上的8条空间剖面线，可得到(b1，a1)、(b2，a2)……(b8，a8)共计8个阈值范围。求取a1、a2……a8的最小值及b1、b2……b8的最大值，分别记作a*、b*，因此取值在(b*，a*)之间的阈值t即为通过这8条空间剖面线获取的合理的阈值范围。理论上，空间剖面线条数越多，两个阈值端点序列会不断“夹逼”，最终将得到最小阈值范围(b*，a*)，从而使得t接近于最佳阈值。但是由于人工规划的空间剖面线的条数是有限的，且考虑到剖面线的空间布局和计算量的问题，这里将空间剖面线条数定为8，且将T＝(a*+b*)/2作为最佳分割阈值。

结合NDWI影像和空间剖面线栅格文件，通过空间运算，获取剖面线上各点对应的NDWI像元值，这些像元值可以绘制成NDWI曲线，每一条空间剖面线都对应一条连续的像元值曲线。

步骤6，根据最佳分割阈值T，提取水体边界范围。

按照下述公式，通过阈值分割的方法，提取水体范围。提取结果为目标水位为1，背景为0的二值图。

为了说明本方法的有效性，随机选取了一个水库(社上水库)使用人工目视阈值选取法提取水体范围，与本方法的结果做一致性分析，同步分析和误差分析，定量评估本方法的提取精度。结果分别如附图4、附图5所示。

从一致性分析，同步分析可以看出半自动阈值法与手动阈值法得到的结果具有很好的一致性(R²为0.95)和同步性。相对误差为6.3％，即基于空间剖面线的半自动阈值法水体提取精度可以达到93.7％，具有较好的提取精度。说明基于空间剖面线的半自动化阈值选取方法可以真实地反映水体范围，在水体范围动态提取的研究中具有较好的应用前景。

应当理解的是，本说明书未详细阐述的部分均属于现有技术。

应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。

Claims

1.一种基于空间剖面线的半自动水体范围提取方法，其特征在于，包括以下步骤：

步骤1：预处理选取的遥感影像；

步骤2：对处理后的遥感影像做波段运算，计算NDWI影像；

步骤6：获取水体范围提取最佳分割阈值T；

2.根据权利要求1所述的基于空间剖面线的半自动水体范围提取方法，其特征在于：步骤1中对选取的遥感影像数据做辐射定标、几何校正、影像的拼接裁剪预处理。

3.根据权利要求1所述的基于空间剖面线的半自动水体范围提取方法，其特征在于，步骤2的具体实现过程是：使用经过预处理的遥感影像的绿光波段和红光波段的值，根据NDWI计算公式，通过波段运算，计算NDWI影像。

4.根据权利要求1所述的基于空间剖面线的半自动水体范围提取方法，其特征在于：步骤3中，针对面积大于1km²，水体长宽比小于10的水体；每幅遥感影像布设不少于8条空间剖面线，布设位置位于水体常年平均水边线附近，剖面线长度约为水体平均宽度。

5.根据权利要求1所述的基于空间剖面线的半自动水体范围提取方法，其特征在于，步骤3中，空间剖面线的规划遵循以下几个原则：

原则1：空间剖面线必须同时覆盖在目标水体和陆地上，不得全部覆盖在单一地物上，可垂直于水体走向布设；

原则2：空间剖面线的布设应充分考虑地表类型的复杂性和差异性，尽量均匀布设并覆盖多种地物；

原则3：空间剖面线条数不宜太少以获取尽可能小的阈值范围，从而逼近最佳分割阈值。

6.根据权利要求1所述的基于空间剖面线的半自动水体范围提取方法，其特征在于：步骤5中所述剖面线栅格文件的属性表，即剖面线上各点对应的NDWI像元值，这些像元值能绘制成NDWI曲线，每一条空间剖面线都对应一条连续的像元值曲线。

7.根据权利要求1-6任意一项所述的基于空间剖面线的半自动水体范围提取方法，其特征在于，步骤6的具体实现过程是：通过运算获取NDWI曲线上斜率变化最快的两点A、B，像素值分别记作a1、b1(a1＞b1)，则(b1，a1)为通过该条空间剖面线所得到的阈值范围，记作阈值t1∈(b1，a1)；同理，对于同一幅影像上的n条空间剖面线，可得到(b1，a1)、(b2，a2)……(bn，an)共计n个阈值范围；求取a1、a2……an的最小值及b1、b2……bn的最大值，分别记作a*、b*，因此取值在(b*，a*)之间的阈值即为通过这n条空间剖面线获取的合理的阈值范围，取T＝(a*+b*)/2为最佳分割阈值。

8.根据权利要求1-6任意一项所述的基于空间剖面线的半自动水体范围提取方法，其特征在于：步骤1中所述遥感影像为GF-1WFV的遥感影像。

9.根据权利要求1-6任意一项所述的基于空间剖面线的半自动水体范围提取方法，其特征在于，步骤7中，水体范围提取公式为：

提取结果是目标水位为1，背景为0的二值图。