一种结合光谱特征和汇流累积量的流域水体提取方法与流程

本发明涉及遥感影像水体信息解译技术领域，特别是一种结合光谱特征和汇流累积量的流域水体提取方法。

背景技术：

地表水体是地球水循环的重要组成部分，获取流域的地表水体分布对水文分析至关重要。由于实地标注的巨大工作量，流域地表水体提取一般通过数字方法进行提取。目前光学遥感影像已成为地表水体监测的常规手段，其采集的数据能够提供宏观的、实时的、动态的水体分布信息且具有很高的成本效益。

mcfeeters在1996年提出了归一化差异型水体指数(ndwi,normalizeddifferencewaterindex)，该方法利用水体在不同波段反射率的不同，通过可见光波段和近红外波段的反差构成ndwi。ndwi可以快速地提取影像中的水体信息，但是此种方法受到密集建筑物的干扰，限制了其实际应用。徐涵秋在2006年提出了修正的归一化差异水体指数(mndwi,modifiednormalizeddifferencewaterindex)，使用短波红外波段来替代近红外波段，能更好地抑制植被信息，同时还可以减小水体与建筑物之间的光谱混淆，更有利于城市水体的提取，但其无法有效区分阴影、道路和其他深色物体。feyisa在2014年提出了自动水体提取指数(awei,automatedwaterextractionindex)，该指数能根据不同背景条件进行形式变换，从而在城市区和山区的精度都优于mndwi。

尽管水体指数是一种直观、简洁的水体提取方法，但由于水体指数不可能完全将水体与其他地物分离，直接应用某一阈值进行水体提取往往会产生大量的错误。针对此问题，isikdogan在2017年提出一种自动水系提取的引擎(rivamap)，该模型基于canny边缘检测理论，使用多尺度奇异指数对水体指数图进行强化，有效减少了裸岩和植被对水体提取的影响。然而rivamap应用于10m分辨率的sentinel-2影像时会将大量陆地单元识别为水体单元，产生了很大的误差。

流域数字水系提取也可通过数字高程模型(dem,digitalelevationmodel)完成，常见的方法是对dem进行处理获得其汇流累积量，然后设定一个阈值进行水系提取。此种方法简单易行，可以直接生成连通的河网，但是提取水系的最佳阈值往往难以确定，且受限于dem精度，此方法提取的水系常常不能很好地反应实际水系分布的情况。现有基于光谱特征的水体提取方法会产生大量假阳性误差(fp,falsepositive)的问题。

技术实现要素：

本发明所要解决的技术问题是克服现有技术的不足而提供一种结合光谱特征和汇流累积量的流域水体提取方法，本发明提高了流域水体提取的准确度和科学性。

本发明为解决上述技术问题采用以下技术方案：

根据本发明提出的一种结合光谱特征和汇流累积量的流域水体提取方法，包括以下步骤：

步骤1、获取哨兵2号多光谱影像的第3波段、第8波段、第11波段和第12波段，对上述波段进行预处理并生成自动水体提取指数图aweinsh；其中，第3波段为绿色波段，第8波段是近红外波段，第11波段是短波红外1波段，第12波段是短波红外2波段；

步骤2、对aweinsh分析并提取流域的水体分布，输出流域的水体分布中心线图worigin；

步骤3、利用dem数据计算流域内各个栅格的坡度及汇流累积量并建立汇流累积量a与平均坡度的相关关系曲线；

步骤4、以相关关系曲线第一个转折点处的汇流累积量作为形成水体的最低限度进行水体分布中心线图worigin的校正，输出最终的水体分布中心线图wrevised。

作为本发明所述的一种结合光谱特征和汇流累积量的流域水体提取方法进一步优化方案，所述步骤1中进行预处理的过程，具体如下：

使用双线性插值法将短波红外1波段和短波红外2波段重采样至与绿色波段和近红外波段相同的空间分辨率。

作为本发明所述的一种结合光谱特征和汇流累积量的流域水体提取方法进一步优化方案，所述步骤2中对aweinsh分析并提取流域的水体分布，输出流域的水体分布中心线图worigin，具体如下：

将aweinsh输入至基于python的rivamap工具包进行分析并提取流域的水体分布，输出流域的水体分布中心线图worigin；

作为本发明所述的一种结合光谱特征和汇流累积量的流域水体提取方法进一步优化方案，所述步骤3利用dem数据计算流域内各个栅格的坡度及汇流累积量并建立汇流累积量a与平均坡度的相关关系曲线，具体如下：

步骤3.1、计算流域各个栅格单元的坡度s和汇流累积量a；

a＝{a¹，a²，...，aⁿ}

式中n为汇流累积量不同的栅格单元总个数，a^j为第j个汇流累积量，j＝1，2...n；

步骤3.2、将各个栅格单元的汇流累积量按数值的从小到大进行顺序，求汇流累积量相同单元的平均坡度

式中i为某一汇流累积量对应的n个栅格单元中的栅格编号，si为该汇流累积量下的第i个栅格单元对应的坡度，j为对应n个汇流累积量的编号，为对应第j个汇流累积量的平均坡度；

步骤3.3、采用双对数坐标绘制汇流累积量a-平均坡度相关关系曲线。

作为本发明所述的一种结合光谱特征和汇流累积量的流域水体提取方法进一步优化方案，所述步骤4具体如下：

步骤4.1、观察汇流累积量-平均坡度相关关系曲线，以相关关系曲线的第一个转折点处的汇流累积量作为形成水体的最小汇流累积量acritical；

步骤4.2、以acritical作为汇流累积量的阈值，生成临界水体图wcritical，其中wcritical表示为

步骤4.3、利用双线性插值将wcritical重采样至与worigin相同的空间分辨率；

步骤4.4、输出最终水体分布中心线图wrevised，计算方法如下

wrevised＝worigin×wcritical。

本发明采用以上技术方案与现有技术相比，具有以下技术效果：

(1)本发明通过数字高程模型确定流域生成水体的最低汇流累积量要求，以此校正基于光谱特征提取的水体，这样能有效减少水体提取的假阳性，得到更加精确的流域水体图；

(2)本发明精度高且适用性良好；利用10m分辨率的哨兵2号多光谱影像可以提取出更精密的水体分布，同时本发明针对不同类型、不同大小的流域均能准确有效地提取其水体分布；

(3)本发明所应用的数据公开可靠并易于获取；哨兵2号多光谱影像及数字高程模型数据均全球覆盖，同时方法涉及的流程简洁、关系明确，保证了计算效率的同时也符合了客观规律。

附图说明

图1为本发明的整体流程图。

图2为本发明计算生成的aweinsh图。

图3为本发明计算生成的初始水体中心线图。

图4为本发明计算生成的坡度图。

图5为本发明计算生成的汇流累积量图。

图6为本发明建立的平均坡度-汇流累积量相关关系曲线图。

图7为本发明计算生成的临界水体图。

图8为本发明计算生成的修正后的水体中心线图。

图9为本发明在夺底沟流域的计算结果与哨兵2号rgb影像对比的示意图。

具体实施方式

为了使本发明的目的、技术方案和优点更加清楚，下面将结合附图及具体实施例对本发明进行详细描述。

针对背景技术中的问题，结合“坡度-汇流累积量相关关系图的第一个转折点象征着发散地形向收敛地形的转变，也表明了水体主导生成过程的变化(kirkby，1986；tarbotonetal.，1992；georgiou，1993)”这一论述，可认为坡度-汇流累积量相关关系图的第一个转折点处的汇流累积量代表着生成水体的最低汇流累积量限制，以此对基于光谱特征生成的水体图进行校正，就可以减少水体图中的识别错误，增加水体识别的准确度和科学性。

如图1所示，本发明提出的一种结合光谱特征和汇流累积量的流域水体提取方法，包括以下步骤：

步骤s1、通过谷歌地球引擎(gee，googleearthengine)获取哨兵2号多光谱影像的第3波段、第8波段、第11波段和第12波段，对上述波段进行预处理并生成自动水体提取指数图aweinsh，其中，第3波段为绿色波段，第8波段是近红外波段，第11波段是短波红外1波段，第12波段是短波红外2波段，包括以下步骤：

①使用双线性插值法将短波红外1波段和短波红外2波段重采样至与绿色波段和近红外波段相同的空间分辨率

②计算水体指数aweinsh，对第3波段-绿色波段、第8波段-近红外波段、第11波段-短波红外1波段、第12波段-短波红外2波段进行组合，突出水体并抑制噪声，计算公式如下

步骤s2、对aweinsh分析并提取流域的水体分布，输出流域的水体分布中心线图worigin(如图3)，包括以下步骤：

将aweinsh输入至基于python的rivamap工具包进行分析并提取流域的水体分布，输出流域的水体分布中心线图worigin；

步骤s3、利用dem数据计算流域内各个栅格的坡度(如图4)及汇流累积量(如图5)并建立平均坡度与汇流累积量a的相关关系曲线(如图6)，包括以下步骤：

①计算流域各个栅格单元的坡度s和汇流累积量a。

a＝{a¹，a²，...，aⁿ}(2)

式中n为汇流累积量不同的栅格单元总个数，a^j为第j个汇流累积量，j＝1,2...n。

②将各个栅格单元的汇流累积量按数值的从小到大进行顺序，求汇流累积量相同单元的坡度平均值

式中i为某一汇流累积量对应的n个栅格单元中的栅格编号，si为该汇流累积量下的第i个栅格单元对应的坡度，j为对应n个汇流累积量的编号，为对应第j个汇流累积量的平均坡度。

③采用双对数坐标绘制平均坡度-汇流累积量a相关关系曲线。

步骤s4、以相关关系曲线第一个转折点处的汇流累积量作为形成水体的最低限度(如图7)进行水体图的校正，输出最终的水体分布中心线图wrevised(如图8)，包括以下步骤：

①观察平均坡度-汇流累积量相关关系曲线，以相关关系曲线的第一个转折点处的汇流累积量作为形成水体的最小汇流累积量acritical。

②以acritical作为汇流累积量的阈值，生成临界水体图wcritical，其中wcritical可表示为

③利用双线性插值将wcritical重采样至与worigin相同的空间分辨率。

④输出最终水体分布中心线图wrevised，计算方法如下

wrevised＝worigin×wcritical(6)

以拉萨夺底沟流域为例，流域面积68.1km²。本例的研究区中运用到sentinel-2msi10m/20m分辨率光学影像与srtm30m分辨率数字高程数据，数据均来源于googleearthengine(https://earthengine.google.com/)。

步骤一、通过谷歌地球引擎(gee，googleearthengine)获取哨兵2号多光谱影像的第3波段、第8波段、第11波段和第12波段，对上述波段进行预处理并生成自动水体提取指数图aweinsh，其中，第3波段为绿色波段，第8波段是近红外波段，第11波段是短波红外1波段，第12波段是短波红外2波段，包括以下步骤：

①将夺底沟流域边界引入gee并获取其对应区域的sentinel-2msilevel-2a的绿色波段、近红外波段、短波红外1波段、短波红外2波段和srtmdem30m数据，对短波红外1波段和短波红外2波段进行双线性插值使其分辨率与绿色波段一致(20m→10m)。

②根据公式(1)计算aweinsh(如图2)。

步骤二、对aweinsh分析并提取流域的水体分布，输出流域的水体分布中心线图worigin(图3)，包括以下步骤：

将aweinsh输入rivamap代替默认的mndwi，默认高斯卷积核大小和处理尺度，设置输出选项，生成流域水体分布中心线图worigin。

步骤三、利用dem数据计算流域内各个栅格的坡度(如图4)及汇流累积量(如图5)并建立平均坡度-汇流累积量a的相关关系曲线(如图6)，包括以下步骤：

①使用arcgis导入dem数据，进行填洼后，采用d8流向算法生成水流流向图并以此计算汇流累积量，之后使用坡度工具计算栅格单元的坡度。

②导出流域各栅格单元的汇流累积量和坡度值，并按照公式(2)(3)(4)所示，按照汇流累积量进行排序，计算相同汇流累积量对应坡度的平均值。

③根据平均坡度和其对应的汇流累积量生成散点图，采用双对数坐标。

步骤四、以相关关系曲线第一个转折点处的汇流累积量作为形成水体的最低限度(如图7)进行水体图的校正，输出最终的水体分布图wrevised(如图8)，包括以下步骤：

①分析平均坡度-汇流累积量相关关系图的第一个转折点，一般出现在汇流累积量较小的位置，代表了控制水流生成的主导过程发生了变化，以该位置处的汇流累积量作为生成水体的最低阈值acritical。

②以acritical为阈值生成wcritical，计算公式如(5)所示，计算结果如图7所示。

③将wcritical与rivamap生成的水体分布中心线图相乘，生成最终的水体分布图wrevised(如图8)，计算公式如(6)。

④将wrevised与对应流域的哨兵2号rgb影像进行对比(如图9)，可以发现本发明有效降低了原水体分布中心线的假阳性并保持了水体分布的准确性。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围内。