基于ASTER GDEM V2黑龙江省牡丹江流域水系提取与分级

2019-09-27李志飞徐金忠杨庆楠

水利科学与寒区工程 2019年5期

李志飞，徐金忠，杨庆楠

(黑龙江省水利科学研究院，黑龙江哈尔滨 150080)

随着“河长制”工作在黑龙江省的全面推行，河流生态系统保护和修复的现实需求更加迫切，河流生态系统健康评价与生态修复成为当前的研究热点[1-4]。河流是一个完整的连续体，上下游、左右岸构成一个完整的体系，需要运用系统思维综合治理。水土保持是我国生态文明建设的重要组成部分，是江河治理的根本，河流两岸水土保持生态体系对河流生态保育与“河长制”实施的科学支撑也成为研究趋势与重点之一。

牡丹江为松花江干流右岸最大支流，位于黑龙江省东南部山地区，流域内水土流失较为严重，牡丹江的生态安全对流域内防洪安全、饮水安全、生态环境等具有重要作用。本文基于ASTER GDEM V2数据，利用ArcGIS的Hydrology水文处理工具对黑龙江省牡丹江流域内水系分布与分级进行了研究，为进一步研究不同尺度河流生态结构与功能、河流生态保护和修复、水土流失防治等提供基础依据。

1 研究区概况

牡丹江为松花江干流右岸最大支流，发源于长白山脉白头山之北的牡丹岭，流经吉林省敦化市和黑龙江省宁安市、牡丹江市、海林市、林口市、依兰县等县(市)，在依兰县城西注入松花江。主要支流有珠尔多河、海浪河、五林河、乌斯浑河、三道河等。流域总面积37 023 km2，黑龙江省境内流域面积28 543 km2。按照全国水土保持三级区划，属于“长白山山地水源涵养减灾区”；按照黑龙江省水土保持三级区划，属于“东南部山地水源涵养减灾区”。

研究区位于黑龙江省东南部，属于中温带大陆性季风气候，春季短，回暖快，风大易旱；夏季温热多雨；秋季短，降温快；冬季漫长寒冷。年均气温4.3 ℃，年均降雨量540 mm左右，年均径流深150 mm左右。海拔介于0～1710 m，地形以山地和丘陵为主，呈中山、低山、丘陵、河谷盆地 4 种形态。土壤类型主要有暗棕壤、白浆土、草甸土等。

2 数据下载

本文数据采用ASTER GDEM数据，数据来源于中国科学院计算机网络信息中心地理空间数据云平台(http://www.gscloud.cn)。ASTER GDEM数据产品是基于“先进星载热发射和反辐射计(ASTER)”数据计算生成，是目前唯一覆盖全球陆地表面的高分辨率高程影像数据。ASTER GDEM V2全球数字高程数据于2015年1月6日正式发布。数据时期：2009年；数据类型：TIFF；投影：经纬度；空间分辨率：30 m；共享方式：免费。

通过查询流域范围，共下载ARTGTM_N43E128 、ARTGTM_N43E129 、ARTGTM_N43E130、 ARTGTM_N43E131、 ARTGTM_N44E128 、ARTGTM_N44E129 、ARTGTM_N44E130 、ARTGTM_N44E131 、ARTGTM_N45E128 、ARTGTM_N45E129 、ARTGTM_N45E130 、ARTGTM_N45E131 、ARTGTM_N46E129 、ARTGTM_N46E130等14幅图。把上述图件拼接，用黑龙江省牡丹江流域界线裁剪，得到原始DEM数据。

3 水系提取

水系提取的基本过程分为以下4步：(1)原始DEM洼地填充。(2)计算水流方向。(3)计算汇流累积量。(4)水系的生成。

3.1 填洼(Fill)

填洼最主要的目的是计算栅格流向，保证水系的完整与准确。由于DEM误差以及一些特殊地形的影响，使得DEM表面存在一些凹陷的区域，导致在水流方向计算时，得到不合理的甚至错误的水流方向。因此，在水流方向计算之前，应该对原始DEM数据进行洼地填充，得到无洼地的DEM[8]。

(3) 单边供电模式下2列AW3车同时起动电流上升率(0.38 A/ms)，电流上升斜率较低，且电流上升持续时间为秒级数据，远远低于保护整定数据，未达到电流增量保护的整定值(25 A/ms)。

对DEM数据进行填洼，输出无洼地的demfill栅格数据。

操作步骤：ArcMap→ArcToolbox→Spatial Analyst Tools→Hydrology→Fill。

3.2 水流方向的提取(Flow Direction)

水流方向通过计算中心格网与邻域格网的最大距离权落差来确定。ArcGIS中水流方向是利用D8算法(最大坡降法)，用邻域栅格的坡降来确定水流方向，坡降最大的单元格即为水流出的方向。每一格网的水流方向指水流离开此网格的指向。流向的生成是个自动的过程，运算的时间跟电脑性能和DEM图的精度大小有关。

用填洼后的demfill栅格数据提取水流方向，输出Fdirection栅格数据。

操作步骤：ArcMap→ArcToolbox→Spatial Analyst Tools→Hydrology→Flow Direction。

3.3 汇流累积量计算(Flow Accumulation)

汇流累积量是基于水流方向计算得出的。汇流累积量的基本思想是以规格网表示的DEM每点处有一个单位水量，按照水流由高处流向低处的自然规律，根据区域流向栅格计算每点流过的水量，从而得到该区域的汇流累积量。当汇流量累积到一定程度的时候，地表就有水流产生，所有汇流量超过临界值的栅格就是潜在的水流路径，这些水流路径构成水系。目前，常用的水系提取方法是地表径流漫流模型[8]。

用Fdirection栅格数据计算汇流累积量，输出Faccumulation栅格数据。

操作步骤：ArcMap→ArcToolbox→Spatial Analyst Tools→Hydrology→Flow Accumulation。

3.4 使用栅格计算器(Raster Calculator)，获得水系

集水阈值的确定影响着水系提取的精度，是自动提取的水系是否与实际河道相符的关键。阈值越大，水系越稀疏，会导致部分常流河道缺失；阈值越小，水系越密集，容易生成伪河道。“试误法”是确定集水阈值常用的方法，其缺点是受人为影响较大，取值较为随意，直接影响水系提取效果。采用水系密度变化率等于集水阈值变化率的数学方法确定最佳集水阈值，即将水系密度与集水阈值的幂函数同切线方程联立，求解唯一值作为水系提取的阈值，这种方法避免了人为干扰，相对客观，效果最佳[8]。

选取0.005 km2、0.01 km2、0.02 km2、0.05 km2、0.1 km2、0.2 km2、0.4 km2、0.6 km2、0.8 km2、1.0 km2、1.2 km2、1.4 km2、1.6 km2、1.8 km2、2.0 km2共15个集水阈值，分别计算出水系密度，集水阈值与水系密度拟合成幂函数见图1。

图1 集水阈值与水系密度的关系

采用公式(1)为切线方程：

y=-x+a

(1)

式中：y为水系密度；x为集水阈值；a为截距。

集水阈值与水系密度的幂函数同切线方程联立：

当x有唯一解时，所取的值即为最佳集水阈值，经计算，当a=0.188 05时，x有唯一解0.062 77。因此，选取最佳集水面积阈值为62 770 m2。

使用栅格计算器，利用Con(“Faccumulation” >= 62770,1)，输出stream栅格数据。

操作步骤：ArcMap→ArcToolbox→Spatial Analyst Tools→Map Algebra→Raster Calculator。

4 水系分级

根据生成的水系栅格数据和水流方向栅格数据，进行水系分级。最小的、长期有水流动的河段是1级河流；两条1级河流汇流成2级河流；两条2级河流汇流成4级河流；下一级河流的等级为上一级河流等级相加，即Shreve分级法。

4.1 水系分级(Stream Order)

水系分级是对一个线性的河流网络进行分级别的数字标识。在 ArcGIS 的水文分析中，提供两种常用的水系分级方法: Strahler 分级和 Shreve 分级。两者的不同之处在于 Strahler 方法中当两条支流汇合时如级别不同则下一级河段的级别不增加，而保留上一级的最大级别，只有当级别相同时才进一级。Shreve 方法中支流汇合后，下一级河流的等级为上一级支流等级相加的结果。本文采用Shreve分级法。

利用Stream和Fdirection，输出streamorder栅格数据，分级方法选Shreve。

操作步骤：ArcMap→ArcToolbox→Spatial Analyst Tools→Hydrology→Stream Order。