基于DEM的西天山中部河流提取研究

2020-10-19杜青松李国玉李金明

中国农村水利水电 2020年10期

杜青松，李国玉，李金明，周宇

(1.中国科学院西北生态环境资源研究院冻土工程国家重点实验室，兰州 730000；2.中国科学院西北生态环境资源研究院大兴安岭冻土工程与环境观测试验研究站，黑龙江加格达奇 165000；3.中国科学院大学，北京 100049)

0 引言

河流作为一种降水的自然特征，与降雨量、地形等因素息息相关。数字高程模型(Digital elevation model，DEM)是一种地表地形精确有效的数字表达，包含丰富的地形、地貌和水文信息[1-5]，无洼地DEM数据是提取自然界河流的基础数据[3,6,7]。基于ArcGIS的水文分析扩展模块，可以快速提取自然界的河流信息[3, 5]，这一方法是目前河网提取最行之有效的方法之一。自O′callagehan和Mark提出河网识别法[8-10]以来，学者们通过坡面流模拟，设定合适阈值快速提取自然界河网这一方法在国内外已开展大量工作[11-13]并取得大量成果[1,14,15]。通过DEM数据提取河网信息不仅适用于自然山区[14]，对喀斯特地貌的小流域提取[15]也同样适用。但在平坦地形地区，该方法提取的河网不理想，存在伪河流现象。

基于DEM数据提取河网的难点在于集水(阈值)面积的确定，对于同一流域利用不同集水面积阈值得到的河网不同[16]，合适的集水面积对准确河流地提取至关重要。集水阈值的确定方法大致可分为试错法[17, 18]、数理统计法[3,14,19-21]、地形分维数法[22, 23]等。各方法在不同区域均得到应用，但都存在一定的人为主观因素影响。其中，由于试错法以定性的方式去描述提取的河网，需要大量的经验，主观性较强，产生的误差最大。

本文在综合前人研究的基础上，基于数理统计方法，以试错的方式提取不同集水面积下的河网，拟合集水面积阈值-河网密度曲线，采用当河网密度变化与集水面积变化相等[24]时求解方程确定集水阈值，用该阈值提取研究区的河网信息，将提取的河网与实际河流对比作精度验证，计算河流的霍顿参数并对该研究区域内矿区通过河流对周围环境造成的影响做简要分析。

1 研究区概况

1.1 研究区边界确定

在传统方法中，对于研究区边界的确定一般是以行政单位为边界，或直接画一个包含研究目标物的矩形框作为研究区的边界。河流是一种自然要素，为了保证其本身的自然属性，通过上诉的两种方法均不能得到自然界完整的流域，因此为了保证研究区河网的完整性，本文对研究区的边界获取通过“盆域”工具获取。“盆域”工具通过DEM数据和流向数据，可以提取河流完整的流域边界，将通过“盆域”分析后得到的栅格矢量化可得到研究区的边界。

1.2 研究区自然概况

研究区位于西天山中部，行政区域隶属于新疆巴音郭楞蒙古自治州和静县。地理位置(见图1)42°59′59.96″～43°34′14.13″N，84°51′22.21″～85°28′33.32″E。面积1 280.84 km2,海拔2 462～5 179 m，平均海拔3 406.79 m，属于高寒高海拔地区。现场观测资料表明该区年最高温度16 ℃，最低气温-30 ℃左右，日最大降水量为 146 mm，年降雨量超过1 000 mm，蒸发量425 mm，最大风速 12 m/s，平均湿度约 43%，风向以NEE为主。地表覆盖以高山草甸为主，以草本为主，近流域范围内存在少量灌木，无高大乔木。海拔3 700 m向上为冰原带，几乎无植被生长。草原是下游牧民放牧的根基，草生长的好坏直接决定牧民的收入水平，草甸生长需要的水、畜牧用水主要来源于高海拔地区的冰雪融水。

图1 研究区DEM分布图Fig.1 The DEM data of study area

研究区被一条分水岭分为两个流域，分水岭以北的玛纳斯河流域和以南的众多细小河流流域，最终汇入伊犁河。分水岭处有一个大型露天铁矿(新疆诺尔湖铁矿)，铁矿由位于河流上游的矿区、选矿区和位于河流下游的生活区组成。在矿区的东北方向有一个自然湖泊诺尔湖，湖泊里的水由西向东注入玛纳斯河。河水主要来源于冰雪融水和位于地表下的石冰川和冻土融化，河流汛期为每年8月，10月后温度降低，冻土冻结，活动层变浅，矿区以固态降水为主，河流进入枯水期。

2 数据与方法

2.1 数据源及数据预处理

利用SRTM v4 DEM数据(http:∥srtm.csi.cgiar.org/SELECTION/inputCoord.asp)，对下载的数据进行数据镶嵌、坐标转换、数据裁剪等预处理，得到与研究区大小一致的DEM原始影像。SRTM v4 DEM作为最常用的地表数字高程之一，被广泛用于河网提取[25, 26]研究。

由于受空间分辨率和系统误差的影响，DEM数据在生成过程中，会产生一些的错误的“洼地”，这会造成在河流流向提取时，DEM中洼地水流方向的不正确，进而影响河网提取的准确性。因此需要对原始DEM进行“填洼”处理以消除洼地的影响。但并不是所有的洼地都是由数据误差造成的，有些洼地是真实地形的反应，如图1中诺尔湖所在的位置，基于此理论，在对原始DEM作“填洼”处理时要设置合理的填洼阈值。

填洼阈值可通过实际计算得到，具体步骤分为：①对原始DEM应用“流向”工具计算流向数据；②通过“汇”工具流向数据为输入数据计算DEM中的所有洼地；③利用“分水岭”工具提取每个洼地的贡献区域；④使用“分区统计”工具计算出每个洼地贡献区域的最低高程和最低高程；⑤运行“栅格计算器”工具将最、最低高程值作差得到洼地填充阈值；⑥最后操纵“填洼”工具对原始DEM进行 “填洼”迭代运算，得到无洼地DEM。

2.2 研究方法

本文首先设置不同的集水面积阈值，并提取每一集水面积对应的河网信息，计算出其河网密度、河流长度、流域面积等，然后对每组数据作相关性分析，拟合合适的模型。最后借鉴前人研究，采用河网密度变化率等于集水值变化率[24]的数值方法确定该区域正确的集水阈值，该方法最大优点是可以有效地避免人为主观性的影响，求解的集水面积为唯一确定值。其原理是：随着集水面积的增大，位于坡面的河道起点会向流域内相邻的地势平坦区域靠近，位于坡面上的河道数减少，伪河道也会被删除[1]。将河网密度与集水阈值的幂函数同切线方程联立，求解方程组的唯一解作为河网提取的阈值，此时的阈值为最佳河流积水阈值[14, 24]。河流汇流面积通过公式计算：

S=NPix

(1)

式中：S为汇流面积，单位由栅格分辨率决定，小流域单位为m2，大流域单位为km2；N为流量栅格数；Pix为栅格分辨率。

通过ArcGIS10.6软件的水文工具箱分别提取不同集水面积下的河网，在等面积投影坐标系下利用几何计算器统计矢量化后的河网长度和流域面积等，再利用统计工具OriginLab计算河网密度，构建集水面积-河网密度曲线。将拟合的函数与切线方程联立，求解方程组的唯一解作为河网的最佳集水面积。将这一集水面积下提取的河网叠加在Google Earth影像上，与实际河流对比以验证其提取的准确性。方法流程图(图2)，具体操作步骤如下：

图2 方法流程图Fig 2 Method flow

(1)利用水文工具箱的“填洼”工具对原始DEM进行填洼处理，得到无洼地DEM；

(2)通过“流向”工具，采用D8算法提取河流方向；

(3)应用“流量”工具提取河流累计栅格河网栅格；

(4)在“栅格计算器”工具设置不同的汇流栅格数(条件函数等)，提取对应的栅格河网，同时计算该栅格数对应的集水面积；

(5)将(4)得到的栅格河流利用“河流链接”和“河流分级”工具进行河网分级处理，分级采用斯特拉勒(STRAHLER)分级法；

(6)输入分级的河网数据，通过“河网矢量化”工具将栅格河网转为矢量数据；

(7)使用“分水岭”工具提取流域面积，同时将得到的结果矢量化；

(8)计算河网密度，拟合集水面积-河网密度曲线，联立方程，求解集水阈值；

(9)提取最优集水面积对应的河网，验证其精度。

3 结果与分析

3.1 不同集水阈值的河网提取

ArcGIS水文分析中的流量定义与水文学中的流量含义不一样，它不代表实际流量，该流量仅是流量分析结果的一个栅格累计计数。流量栅格中每一个像元记录的是流向该点的栅格数量总和，流域集水面积与流量栅格数的关系为式1，通过该式可实现流量栅格数与集水面积的相互转换。

以500为间隔，分别设置流量栅格数为500～8 500 共计17组数据，对应集水面积为0.45～7.65 km2之间的以0.45 m2为间隔的17个数据，同时提取各集水面积对应的河网，统计分析各集水阈值与河流数、河流长度、河流流域面积和河网密度之间的关系。

集水面积与河流数的关系如图3所示，集水面积与总河流数和一级河流数相关性较好，与其他级别河流数相关性不明显。随着集水阈值的增加，各级河流数和总河流数均减少，一开始降幅大，在集水阈值为5.4 km2时，总河流数和一级河流条数趋于稳定；在集水面积为1.35 km2时，最高河流等级由5级减少为4级，总河道数由2 039 减少到628条，减少69.2%；当集水面积超过5.85 km2后，最高级的4级河流条数稳定为3条，不再随集水面积的增大而变化。

图3 集水阈值(面积)与河流数关系Fig.3 Relationship between catchment threshold (area) and number of rivers

图4描述的是集水面积与河流长度的关系，分析得出，集水面积与河流长度的关系和集水面积与河流条数的规律类似，总河流长度和一级河流长度受集水面积影响大，集水面积对其他各级河流长度影响小，随着集水面积的变化，1～4级河流长度变化不明显，趋于一个稳定值。总的来说，河流长度随着集水阈值的增大总体呈减小趋势，但减小幅度慢慢变缓；当集水面积超过5.85 km2时，4级河流长度保持3.79 km不变，不再随集水面积的增大而变化，此时与河道数保持不变相对应。

图4 集水阈值(面积)与河长数的关系Fig.4 Relationship between catchment threshold (area) and river length

从图5可以看出，集水面积与流域面积的相关性不高，变化规律不明显，除一级河流流域面积总体上随着集水面积的增大而减小外，其他各级流域对应的流域面积变化无规则，相关性极差，各级河流流域面积变化情况也不一致，有增有减。

图5 集水阈值(面积)与流域面积的关系Fig.5 Relationship between catchment threshold (area) and watershed area

河网密度为单位面积上的河流长度，是河流的一个基本参数，与河流阈值变化息息相关。图6显示集水面积-河网密度曲线，可以得出，集水面积与河网密度的相关性极高(R2=0.999 78)，河网密度随着集水面积的增大而减小，最终趋于稳定。集水面积-河网密度呈现为幂函数图像(图6)，函数表达式为：

y=0.828 24x-0.488 93(R2=0.999 8)

(2)

式中：y为河网密度，km-1；x为集水面积，km2。

图6 集水阈值(面积)-河网曲线Fig.6 Catchment threshold (area)-river network curve

3.2 集水面积阈值的确定

采用河网密度变化与集水面积变化相等时的数值方法确定最佳集水面积阈值。将式子y=-x+a与式(2)联立，求解方程组，得出的唯一解x为最佳阈值。

(3)

求解方程组(3)，有a=1.659 4时，方程组的未知数x存在唯一解0.544 9，由此确定研究区最佳集水面积为0.544 9 km2，计算得出对应的流量累积栅格数为605。同时提取该集水面积阈值下的河网，结果经矢量化后如图7所示。

图7 最佳集水面积下提取的河网Fig.7 River network extracted under optimal catchment area

以最优集水面积0.544 9 km2提取的河网，河流由5种级别组成，1级河流为最基础河流，5级为最高等级。提取的河网存在大量的基础河流，这与研究区的高海拔、起伏大的地形地貌特征相符，即使是在相对平坦的区域也几乎不存在平直的伪河流，这一结果与实际调研情况也相符合。

通过河流与矿区、生活区等的空间拓扑关系及河流的流向可知：矿区开采使用大量的炸药、运输车辆存在漏油、燃油加铅等现象，采矿会对环境造成污染，是一大的污染源。经过河流的作用会将这一点源污染扩展为面源污染。其中采矿区对诺尔湖沿线的水体造成污染严重并通过河流将污染物带入玛纳斯河流域，生活区和下游草地主要受选矿区影响。河流流向对环境污染研究采样点的确定有指导意义，通过对河流沿线和污染源地区采样对比研究，可大致了解研究区域的污染特征。

河流霍顿参数[22]是表征河流与地貌关系的基本参数，在水文模型模拟研究中具有重大意义，是提取是水文模型构建的基础[19, 22]，研究区的河流基础数据和地貌特征参数如表1所示。

研究区总河流共计1 661 条，基础河流发育，占50%；一级河流面积占63%；河网密度为1.11 km-1，河源密度0.65 条/km2；河流的河数率、河长率、面积率分别为1.66、1.83和2.12。该区域由于海拔高、地形起伏大，有利于河流的形成；同时，存在大量冰雪、冻土和石冰川，受气温影响大，故存在大量由固态水融化而形成的细小河流。

3.3 提取河网的精度验证

研究区范围较小，缺乏能与之相匹配的地表基础数据，如河流数据等。故很难采用传统方法将提取的河网与“蓝线河网”对比定量研究[21, 27, 28]得出河流提取的误差。为验证河网提取精度，将提取的河流与Google影像叠加(图8)，对比其与Google影像中实际河流之间的误差。通过对比发现提取的河网与实际河流基本一致，几乎不存在伪河网或河网提取不足的情况。