基于GIS和DEM数据的流域水系提取与优化分析

2021-03-05王成文黄小琴徐兆祥

中国农村水利水电 2021年2期

王成文，李英,2,3，黄小琴，张勃，徐兆祥，李阳

(1.宁夏回族自治区水文环境地质勘察院，银川 750011；2.中国地质大学(北京)水资源与环境学院，北京 100083；3.宁夏回族自治区地质局，银川 750021)

近年来，随着GIS技术的不断发展，区域地形地表水系的提取研究逐步由传统的野外测量和地学统计方法转向高精度的DEM(数字高程模型)数据数字化和自动化提取工作中来。在流域水系数字化提取研究中，如何利用DEM数据更精确真实地反映实际地表水系特征，是当前水文模型研究热点和难点，同时水系的数值化提取研究也是数据信息反演解译实际的重要组成部分[1]。目前，对于流域水系的自动提取研究方法丰富多样，总体可分为两种思路，一种是基于已有遥感图像进行自动提取，另一种是利用DEM数据进行水系自动提取。由于DEM数据分辨率不断更新提升，且能够表达丰富的地貌信息而被广泛研究和应用[2]。

大量研究实践发现[3-8]，在实际地形中凹陷洼地和平缓地带的是普遍存在的，DEM数据中对于栅格单元局部凹陷和平坦地带细节处理不够完善，因此在提取水系之前必须要对DEM进行预处理。此外，基于GIS技术提取水系过程中，如何合理的设定网格数和阈值范围是关乎河流数字化成败的关键因素[9-11]，而以往研究中大多是进行网格和阈值的大范围值域的模拟、筛选，甄选出最佳值。本文在参阅国内外文献和反复实践的基础上，利用ArcGIS10.4软件中水文分析模型，进行基于DEM数据的流域水系提取与优化处理分析，探讨流域水系提取的优化方法，特别是对网格数和阈值的选择进行数学拟合分析，选择合理的网格和阈值，为流域水系特征的数字化提取研究提供一定参考。

1 研究区概况

本文拟以宁夏海原县境内范围为研究区，进行山区流域水系提取。海原县位于宁夏回族自治区中南部，六盘山西北麓，行政区划属中卫市管辖，地理位置介于东经105°09′～106°10′，北纬36°06′～37°04′之间，东西宽约100 km，南北长约95 km，面积4 990 km2，海拔在1 400～2 800 m之间，最高点位于南、西华山山区，最高海拔2 954 m，见图1。

图1 海原县等高线示意图Fig.1 Contour lines diagram in Haiyuan county

2 研究方法与过程

2.1 数据及预处理

本文中的DEM数据是源于“地理空间数据云”信息中心的高程数据，它是集结了ASTER GDEM第一版本(V1)的数据连接而成，数字高程分辨率为30 m。由于ASTER GDEM V1数据中有个别数据异常现象，需结合美国NASA的SRTM高程数据，分辨率为90 m。研究中以30 m分辨率的DEM数据为主，90 m辅助使用，影像图见图2。

图2 研究区DEM影像图Fig.2 DEM image of study area

由于DEM数据是由等高线生成，使得栅格数据中带有一个或一组凹陷点，表现为四周高中间低的洼地，在进行水系提取前必须填充，对DEM数据进行预处理，称为“填洼”处理[12]。利用ArcGIS中水文分析的Fill Sink工具进行填洼，先对每一个栅格单元进行扫描、收集，然后运用Zlimit选项确定合适的填充阈值，该值即为被填补深度的临界值(地形标高)，只有当洼地低于该临界值时才被识别、填充，而超过该临界值的地形标高保留不变以达修复目的，修复结果见图3。

图3 DEM的预处理影像Fig.3 DEM preprocessing images

2.2 流域水系的提取

基于处理后的DEM数据，进行水系提取工作，采用目前广泛应用D8单流向算法，该方法是通过DEM数据判别流域水系水流方向，计算流域汇流累积阈[13]。水系提取过程：①地形分析，②水流方向判别，③栅格数的确定和汇流累计量算，④设定合理集流阈值，提取水系信息，⑤栅格水系矢量化。水系提取技术路线见图4。

图4 水系提取实现技术路线Fig.4 Technical route for basin river extraction

2.2.1 研究区地形分析

坡度、坡向的地形特征是流域水系流向分析的基础。坡度表示地表倾斜程度，坡向是表征某点高程值变幅的最大变化方向[14]。一般运用拟合曲线面法求解，并结合ArcGIS中3D analyst模块对坡度、坡向进行分析、提取，公式如下：

(1)

A=Sx/Sy

(2)

式中：S为坡度；A为坡向；Sx为x方向坡度；Sy为y方向坡度。

坡度、坡向提取结果见图5，研究区内坡度在0.01°～36.72°之间，由西南往北东坡度呈减小趋势，到东北一带坡度<20°，坡向取值在0～360°之间。从坡度坡向可以看出，研究区坡度较大、坡向多变，地形极其复杂，中部、南部等大部分地区为纵深很长的山地，坡向清晰可辨；东部坡度较为平缓，辨析度降低。

图5 研究区坡向坡度影像分析Fig.5 Image analysis of slope direction and gradient in study area

2.2.2 水流方向确定

运用D8算法[11,15]来确定水流方向，首先计算DEM数据每个栅格单元与四周的坡度和坡向关系，然后选择最陡坡度，设定最陡坡度为该单元的水流流向，每个栅格单元中的水流向四周有8个方向流出可能，并用1、2、4、8、16、32、64、128这8个流向编码分别表示各个流出方向，通过对比坡度，选择某一方向，确定水流方向，水流流向编码见图6 ，其最大优势是在基于ArcGIS 等的环境下能够较为快速、准确地运算。

图6 水流流向编码Fig.6 Flow direction code

2.2.3 汇流累积量的计算

利用ArcGIS 水文分析模块下的Flow Accumulation函数确定流域水系的汇流量，设置栅格临界值，识别有效栅格数，计算出有效栅格上累积的汇流栅格数，确定汇水面积，汇流累积量为栅格单元数目与栅格单元面积乘积之和[16,17]。计算识别的有效栅格数的多少可以表征流域水系的汇流能力，即某流域单元内上的汇流栅格数越大，形成的汇流面积越大，形成地表的径流越明显。

2.2.4 流域水系的提取

科学确定符合研究区地形条件的阈值是提取流域水系的关键,以往研究一般设定阈值主要是通过反复对比不同阈值提取的流网，主观判断最佳阈值[18]。本文水系的阈值设定是根据需要对不同级别的河流设定不同阈值，要综合考虑研究区的基本状况、流域地形地貌等方面，分别将汇流量阈值范围为2 500～25 000，提取河网水系，并运用多种数学函数拟合确定合理阈值提取水系。

3 结果分析与优化

3.1 结果验证与误差分析

本文以野外实测数据为准值，通过野外实地勘察测量，室内数字化得到河流的中心线，并以该中心线的长度为真实河长。由实测结果可知，海原县境内共发育清水河和祖历河两大水系，均属黄河支流，按其支流集水区域的不同，自北向南又可分为西河、苋麻河和杨明河等流域区。

本文以海原县境内河流为研究对象，结合ArcGIS中的水文分析模块绘制其流域信息，其流域长度见表1。在基于30 m DEM数据基础上，利用D8算法提取了4级水系，其中一级水系为11条，二级水系为15条，三级水系为8条，四级水系均为11条。提取的水系总长度为350.05 km。为了定量比较提取的结果，表1给出了基于30 m DEM数据提取值和野外实测值的数据对比，研究区流域水系提取效果见图7。对比可得，基于GIS和DEM数据的流域水系提取计算结果与野外实测量算成果较为接近，尤其大中河道长度偏差仅在7.5%以内，中小河流河道长度偏差在9.2%以内，而小河流河道长度的偏差值可达20%以上，可见提取的低级别水系数目的误差较大。总之，从提取的水系长度及数量看，基于GIS和DEM数据提取的流域水系有一定参考意义。

表1 水系提取与高精度DEM测值对比Tab.1 Comparison of basin river extraction and high-precision DEM measurements

图7 研究区流域水系提取影像图Fig.7 Image of basin river extraction in the study area

另外，运用SPSS统计软件分析，对比研究区提取的36条河流，发现DEM数据在流域水系提取应用中，在地形起伏较大的山地、山谷区域的提取识别较好，误差也小，一般水系提取的误差小于5.32%；而在起伏变化不大的丘陵、平原地带提取识别效果一般，且误差变大，提取的误差值一般介于5.5%～12.8%，个别河流误差大于20%，见图8。

图8 不同地区DEM数据提取误差分析Fig.8 Error analysis of DEM data extraction in different regions

3.2 网格数的确定与改进

合理识别有效DEM数据的栅格单元数是影响汇流网格数的主要因素，而网格数的确定是提取水系流域面积的决定性因素，真实流域的水系与基于GIS提取水系的距离误差是随着网格数的变化而变化，通过拟合网格数与距离误差的函数，得到距离误差最小值所对应的网格数，并通过此网格数提取流域水系。在2 000～40 000网格数范围内，选取若干网格数分别生成水系图。在Arcgis10.4下加载DEM数据和提取的水系图，逐一量测提取水系河源与实际河源坐标的距离误差，从而得到距离误差与网格数的相关关系，如图9所示。

图9 距离误差随网格数的变化曲线Fig.9 Range error curve along with grid numbers

由图9可知，当网格数在10 000～15 000范围时，实际河源与提取水系河源之间的距离误差变化呈最小趋势，而当网格数从15 000左右继续增加时，距离误差又呈增大的变化趋势。为了进一步确定最小误差及其对应的网格数，通过Origin9.0中的函数拟合，运用多种函数及偏导数的计算，最终选取采用多项式拟合，函数最次幂为3次，相关系数R=0.976。同时，求取函数式的一阶偏导数，并令其一阶偏导数等于0，再根据网格数的变化区间，确定误差最小时的网格数为12 245，即网格数为12 245时所对应的水系即为该区的水系(见表2)。

表2 距离误差拟合曲线参数Tab.2 Range error fitting curve parameters

3.3 集流阈值格数的合理取值

利用Hydrology水文处理工具，分别设定2 500，5 000，8 000，12 000，18 000，25 000 共6个汇流累积量阈值，生成栅格河网，见图10，并结合不同阈值的流域水系特征，对比实际水系流量，计算其河流长度、汇流面积与河源数等水文信息，运用数学函数拟合分析，综合确定合理阈值。

图10 不同阈值条件下的流域水系提取Fig.10 Basin river extraction under different threshold value conditions

以研究区南部水系的杨明河为例，从累积汇流量图层中提取各阈值下的流域水系，计算不同阈值下的河流长度、汇流面积与河源支流个数等水文信息，详见表3。随着阈值的不断增加，河流长度、汇流面积与支流个数均呈减小的趋势，且不同水文参数的变化趋势各有差异。结合各阈值流域特征拟合曲线分析(图11)，当阈值在8 000之内时，河流长度随阈值的增加呈骤降趋势，由464.36 km急降到88.72 km；当阈值大于8 000时，河流长度随阈值的增加呈缓慢减少趋势，表明河网密度变化较平缓，总体变幅较小。当阈值在12 000内时，流域面积随阈值的增加呈骤降趋势，由3 545.32 km2骤降到552.52 km2；当阈值大于12 000 时，流域面积随阈值的增加呈缓慢减少趋势，表明河流面积变化较平缓，变幅较小。同样，河流支流个数同样呈相识的规律变化，在阈值12 000为河流支流个数的变化拐点，呈“骤降-缓减”的变化趋势。可见阈值在8 000和12 000为该研究区的关节点，尤其是阈值12 000，对整个河流的影响呈显著影响。

表3 不同阈值的流域水系特征(杨明河流域)Tab.3 Basin river characteristic with different threshold values (Yangming River)

图11 不同流域特征因数随阈值变化的拟合曲线Fig.11 Fitting curve of different basin characteristic factors with threshold value

另外，运用函数拟合手段，分别采用幂函数、指数函数、多项式等多种函数对河流长度、流域面积及支流个数进行趋势拟合分析，以确定提取水系流网的最佳阈值。通过多次试验对比，发现河流长度和支流个数采用幂函数(y=AxB)拟合效果较好，相关系数分别为0.981和0.973，曲线的阈值变化拐点范围均为8 000～12 000；而流域面积则更适合采用指数函数拟合，其相关系数为0.998，拟合效果最好，曲线的阈值变化拐点范围同样是在8 000～12 000，详见表4。

表4 不同流域特征因数随阈值变化的拟合参数Tab.4 Fitting parameters of different basin characteristic factors with threshold value

结合研究区地形气候特征和实际测量分析，海原地区为西北干旱山区，常年地表蒸发远远大于降水补给，易选取阈值较大值，即选择水系流网较稀疏，水流面积较小值，更符合干旱山区的实际水系特征。由此可以得出，选取的阈值设为12 000时，模拟的海原地区水系河网与实际水系吻合效果最佳，同时更科学地指导基于DEM数据源的水系和实测水系之间的差异分析。