基于Arcgis-python断面汇水面积批量提取方法研究

洪明海，汪仕伟，曾文治，黄介生，杨荣芳，苏海鹏，李析男

(1. 贵州省水利水电勘测设计研究院，贵阳 550002；2. 武汉大学水资源与水电工程科学国家重点实验室，武汉 430072)

0 引 言

汇水面积是由分水线包围的集水面积，在闭合的汇水面积内，构成一定形状的水路网系统[1]，是设计桥梁涵洞的孔径大小、水库大坝的坝高以及水库来水量等的重要参数。因此，准确地勾绘出设计断面的汇水面积对于水库工程设计、生态环境治理以及水文统计分析提供数据尤为重要[2, 3]。

在地理信息系统软件出现之前，汇水边界主要以地形图上的分水线进行划分[4]，在地形图上依据等高线分布绘制流域的分界线，然后量算出该分界线所包围的面积，效率和精度依赖于人员操作熟练度；同时，对于地形复杂，山脊线和分水岭等分界线难以区分，人工勾绘极易出错进而严重影响汇水面积的成果，卢洋暘等[5]的研究也表明人工勾绘汇水面积的方法在处理复杂山地地形和范围较大地形时存在较大的局限性；此外，人工勾绘的方法对于河道断面位置无相应的地形图或者已有地形图的精度难以满足工程计算的需求时显得力不从心，李晓等[6]也认为需要一种新的方法应用于无地形图地区的工程位置处流域面积的提取和计算。而现在，计算机程序可以从DEM中提取地形要素特征来生成汇水边界，应用计算机技术，只需要传统方法的一小部分时间就可以得到汇水区域的边界[7]。目前，Arcgis、SWAT、MIKE、Mapgis、Civil 3D等众多软件都可以实现基于DEM的水文研究工作，比如汇水面积的勾绘和计算、提取河网信息以及流域特征分析等[8-12]。基于DEM分析计算的方法不依赖地形图，不仅节约时间成本，而且受人为因素影响较小，且在地形复杂地区和大范围地区均具有较好的适用性，但目前该方法多是针对单个断面的一系列复杂操作，自动化程度还不高，对于短时间内需要同时处里大量的断面汇水面积并获得数据结果仍是一个不小的挑战。因此，有必要对该方法做进一步的拓展和研究以满足工程实践的需要。

本文拟采用ESRI公司的Arcgis10.4软件来处理芙蓉江流域内河道断面的汇水面积，同时通过python调用Arcigs中相应的水文工具形成汇水面积处理模型，批量获取河道断面的汇水面积，将计算的结果和已有的数据进行比较分析，同时将该方法应用到洪渡河流域、桐梓河流域、綦江流域、赤水河流域、大渡口以上区域、大渡口以下干流区域以及六盘水区域等7个研究区进行验证，为河道断面汇水面积的勾绘提供一种新的借鉴方法。

1 材料与方法

1.1 研究区概况和数据来源

1.1.1 研究区概况

本次计算的断面均位于贵州省遵义市的芙蓉江流域，根据《中国河湖大典·长江卷》[13]和《贵州河湖》记载，芙蓉江(107°02′E～107°52′E，28°08′N～29°20′N)发源于贵州省遵义市绥阳县枧坝镇，由西南向东北流经贵州省遵义市正安县、道真县和重庆市武隆区、彭水县，于武隆区江口镇汇入乌江，全长234 km，流域面积7 406 km2，西邻綦江水系的桐梓河、松坎河，北靠乌江，南邻湄江，东邻洪渡河；南北长约160 km，东西宽约195 km，流域地势起伏大，东低西高，上游地形较开阔平坦，沿河发育有河谷阶地和河漫滩，中下游河谷深切，河道狭窄，滩多流急，岩溶发育，多洞穴、井泉和伏流暗河段，属中山峡谷及岩溶地貌特征，平均高程约850 m(黄海基面，下同)，总落差1 124 m，河流平均比降5.76‰，流域面积大于1 000 km2的支流有清溪河、三江和梅江，流域面积大于100 km2的支流有新民河、关坝等15条河流。流域属亚热带湿润季风气候区，年平均气温16.4 ℃，年无霜期281.8 d，年日照时数为1 057 h，多年平均降水量1 106 mm，降水多集中在4-10月，占年降水量的73%，年水面蒸发量为603.2 mm。河口多年平均流量为169 m3/s，多年平均悬移质输沙量为277.3 万t。流域多年平均水资源总量53.3 亿m3，干流水力资源理论蕴藏量47.3 万kW。流域内无大型工矿企业，工业污水排放量小，污染物主要来自农田的农药、化肥、养殖业的污染负荷和生活污水。

1.1.2 数据来源

本文使用的dem数据为Aster传感器采集的30 m分辨率的栅格数据，WGS84椭球，数据下载于地理空间数据云(http:∥www.gscloud.cn/search)，由中国科学院计算机网络信息中心提供。河流矢量数据来源于各水务局搜集资料整理，已经项目三级校审制度审核。

根据已有的资料整理，本文分析选取的河道断面均为贵州省遵义市境内的芙蓉江流域上50 km2以上支流的断面，共有38个，分三类，其中，一类是支流汇口断面，如梅江汇口断面、罗家河汇口断面，共有29个，占比76.32%；第二类是已建成水库断面，以坝址为控制断面，如沙坝水库断面、庆垭电站水库断面，共有6个，占总数的15.79%；第三类是河流上的水文站控制断面，有3个，占总数的7.89%，它们是旺草水文站断面、长坝水文站断面和五家院子水文站断面，从图1可以看到本次分析选取的三种类型共38个河道断面的分布情况，其中水文站类型的3个断面有在芙蓉江干流上游的(旺草水文站断面)，有在芙蓉江干流下游的(长坝水文站断面)，还有位于芙蓉江一级支流梅江上的(五家院子水文站断面)；6个水库型断面则分布在贵州省境内的芙蓉江最大的三条一级支流清溪河、三江河梅江上；支流汇口型断面则广泛分布于芙蓉江流域的支流上，因此，本次计算分析选取的河道断面具有较强的代表性。

研究区域的河流水系情况和本次分析选取的38个河道断面分布情况如图 1所示。

38个断面的汇水面积数据均来源于水利普查资料，如表 2中“已有资料面积”栏所示，其汇水面积最大值为长坝水文站，为5 454 km2，最小值为铁钉岩水库，汇水面积为7 km2。本文根据中华人民共和国国家标准 GB/T 4883-2008(数据的统计处理和解释·正态样本离群值的判断和处理)[14]，采用标准推荐的在未知总体标准差的情形下的离群值处理方法—格拉布斯(Grubbs)检验法对已有资料的流域面积进行检验，使用双侧检验，其计算公式如下所示：

图1 研究区域河流水系及河道断面分布情况Fig.1 Distribution of river systems and river sections of study area

(1)

(2)

(3)

(4)

表1 格拉布斯双侧检验计算结果表Tab.1 Table of calculation results of grubbs bilateral test

此外，本文提出的方法还在洪渡河流域、桐梓河流域、綦江流域、赤水河流域、大渡口以上区域、大渡口以下干流区域以及六盘水区域等7个研究区进行验证，该7个研究区的汇水面积数据也包含水文站、水库和汇口等3种类型的断面数据，且所有7个研究区的数据均采用格拉布斯双侧检验法进行检验，保证对比资料在统计学上真实有效。

1.2 计算原理

基于DEM流域划分和河流水系的提取一般都会对DEM栅格数据进行填洼处理以及计算流向，此外，当获取的外部矢量数据(比如水文站等断面数据)和DEM不是来自同一份数据的话，矢量数据和DEM之间就有产生一定的误差，这时就需要使用捕捉倾泻点功能或者编辑功能对其适当修正。

洼地(或突起)是由于数据的分辨率或者将高程四舍五入到中心像元最近的整数值而产生的常见错误，应对洼地进行填充以确保盆地和河流得以正确划界，如果未对洼地进行填充，则生成的河流水系网络可能会呈现不连续性。填洼工具使用与焦点流、流向、洼地、分水岭和区域填充等工具等效的功能来定位和填充洼地，该工具的执行过程会进行迭代，直到指定 z限制内的所有洼地均填充完毕，在填充洼地的同时，可能会在填充区域的边界处创建其他洼地，这些洼地将在下个迭代中移除[15]，其原理示意图如图2。

图2 填洼原理示意图Fig.2 Schematic diagram of fill principle

图3 流向原理示意图Fig.3 Schematic diagram of flow direction

1.3 研究技术路线及步骤

本文的研究技术路线如图 4所示，其主要操作步骤如下：

(1)对研究区DEM进行裁剪投影处理。

(2)根据DEM提取研究区的河流水系，其中填洼处理和计算流向以及设置流量阈值生成河流是比较关键的步骤，填洼影响流向的形成，计算流向决定了流量的汇集方向，而流量阈值的设置直接决定着河流水系的密度。

(3)对37个河道断面数据进行处理分析，一般导入Arcgis的矢量河道断面数据和所使用的DEM数据或多或少都会有一定的偏差，这时需要使用捕捉倾泻点工具或者编辑工具并参照已有的矢量河流来适当调整以满足软件的计算需要。

(4)以上3步都处理好之后就可以使用分水岭工具开始计算断面的汇水区域了，注意断面汇水区域之间的重叠情况，这里使用python批量分割出每个断面单独计算之后再合并到一起，避免了上游断面的汇水区域未包含在下游的汇水区域中的特殊情况，此外计算得到的断面汇水区域是栅格，需要将其转为矢量格式。

(5)计算矢量汇水区域的各种属性值，这里主要计算38个断面的汇水面积值，需要在投影坐标系下方可计算。

图4 本文研究技术路线图Fig.4 Technology roadmap

1.4 分析指标

评价指标选择绝对误差Δ、相对误差m、决定系数R2、均方根误差RMSE、相对均方根误差RRMSE、相对分析误差RPD、偏差Bias：

Δ=Xgis-Xdata

(5)

(6)

(7)

R2=R×R

(8)

(9)

(10)

(11)

(12)

2 结果分析

2.1 汇水面积分析

表2给出的是python调用Arcgis工具批量计算的芙蓉江流域内37个河道断面汇水面积(下称计算面积)和已有资料的汇水面积(下称已有面积)结果，其中，按照公式(5)和公式(6)分别计算了计算面积和已有面积的绝对误差和相对误差。表2中，旺草水文站(8号)位于芙蓉江干流上，但其位于上游，从图1中看到其上只有罗家河一条支流汇入，其汇水面积小得多(计算面积为401.83 km2，已有面积为413 km2)(表 2)，相对误差为2.7%；此外，我们在芙蓉江支流梅江上也选取了一个水文站控制断面—五家院子水文站断面(图 1)，其计算面积为1 212.14 km2，已有面积为1 230 km2，相对误差为1.45%(表 2)；2个水文站控制断面中，绝对误差最大的是五家院子水文站控制断面，为17.86 km2，相对误差最大的则是旺草水文站控制断面，为2.7%(表2)。

而6个水库型断面均位于芙蓉江的支流上，贵州省境内芙蓉江最大的三条一级支流清溪河(清溪水库断面)、三江(庆垭电站水库断面)和梅江(洋渡电站水库断面和大河坝水库断面)均有分布(图 1)，其中铁钉岩水库断面位于巴渔河上(伏流型河流)。6个水库型断面中，汇水面积最大的是清溪水库断面，计算面积为1 237.44 km2，已有面积为1 261 km2(表 2)；绝对误差最小的巴渔河上的铁钉岩水库断面，为0.54 km2(表2)；相对误差最小的为洋渡电站水库断面，为0.79%，最大的为铁钉岩水库断面，为7.66%。

本次计算分析中最多的是支流型汇口断面，共有29个，占比76.32%，广泛分布在芙蓉江的各支流上。其中汇水面积最小的是芙蓉江一级支流茨竹河汇口断面(图 1)，计算面积为48.64 km2，已有面积为50 km2(表 2)；而最大的则是清溪河汇口断面，计算面积为1 493.35 km2，已有面积为1 504 km2(表 2)。29个支流型汇口断面中，绝对误差最小的是池武溪汇口断面，为0.08 km2，相对误差最小的是梅江汇口断面，为0.09%(表2)。王峥等[18]人采用DEM对泾河流域的1 000 km2以上部分支流的汇水面积进行提取并做了比较分析，绝对误差也达到17 km2。其余河道断面对应的汇水面积以及误差等数据请参见表 2。

2.2 计算方法的验证和误差分析

37个河道断面的汇水面积数据的相对误差中，最小的为梅江汇口断面，为0.09%，最大的则是铁钉岩水库断面，为7.66%。图5给出的是已有资料的汇水面积和python批量计算的汇水面积相关性曲线图，从图中可以看到基本上所有的数据点均落在1∶1线上，计算表明其决定系数R2为0.99，而决定系数是用来衡量理论值与计算值的离散程度，其值范围0～1，越接近于1，效果越好，因此通过决定系数来看通过python调用Arcgis工具批量处理的方法精度很高。此外，本文也计算了37个数据点的均方根误差RRMSE、相对分析误差RPD以及偏差Bias，相对均方根误差是用来衡量残差的相对大小，该值越小，说明计算结果越好[19]，而数据点计算得到的相对均方根误差仅为2.74%，表明计算结果很好；相对分析误差大于2，说明计算结果很好，相对分析误差小于1.4，说明计算较差，相对分析误差介于1.4～2，说明计算较可靠，还可以进一步调整[20]，37个数据点的相对分析误差为49.01，远大于2，表明计算结果很好；37个河道断面的汇水面积数据的偏差为2.88，偏差为理论值与计算值的偏差程度，大的偏差值表明计算值没有很好地逼近理论值，计算效果差。

表2 已有资料和python批量计算的断面汇水面积误差表Tab.2 Table of error of available data and python-calculation

注：类型栏1代表水文站，2代表水库，3代表汇口。

图5 已有资料和python批量计算的汇水面积对比Fig.5 Comparison of catchment area of available data and python-calculation

同时，本文还将该方法应用其他研究区域以进一步验证通过python调用Arcgis批量计算汇水面积的可靠性。除本文建立方法使用的芙蓉江流域以外，还在洪渡河流域、桐梓河流域、綦江流域、赤水河流域、大渡口以上区域、大渡口以下干流区域以及六盘水区域等7个研究区广泛选取河道控制断面来验证该方法，选取的断面类型包含汇口、水库和水文站等3类。同样，本文也采用格拉布斯双侧检验法对新选取的研究区已有的汇水面积资料进行检验，排除统计离群值的存在，保证用于比较分析的汇水面积在统计学上真实有效。之后，根据选取的断面数据采用本文提出的方法批量计算了汇水面积数据，然后计算了平均相对误差(m)、决定系数(R2)、相对均方根误差(RRMSE)、对分析误差(RPD)以及偏差(Bias)列于表3。从表3中可以看出，相对误差平均值均较小，最大的仅为7.76%，而决定系数均在0.95以上，表明计算得到的汇水面积数据和已有的汇水面积数据具有很高的相关性，同时相对均方根误差也很小，最小的为桐梓河流域2.1%，7个验证研究区的相对分析误差均大于2，说明该方法计算结果好，计算方法合理，可用于对断面汇水面积的计算。其余统计数据请参见表3。

尽管本次选取的37个河道断面的汇水面积和已有资料的汇水面积数据对比结果很好，相对误差小，决定系数高，但在模型搭建和运行处理过程中也发现了许多需要进一步完善和改进的地方，误差产生来源主要有以下4点，列于表4：第一，所使用的DEM产生的误差。本文采用的DEM分辨率为30 m×30 m，和真实的现实地形相比，在某个像元处可能其地形会被概化为一个值，会产生不小的误差，如果条件允许，可以考虑10 m×10 m、1 m×1 m甚至更高分辨率的DEM数据，但这样也会增加模型的运行时间，同时还会对计算机等硬件性能有更高的需求。第二，对比值产生的误差。本文所分析的两组河道断面的汇水面积数据，一组来源于我们通过python调用Arcgis计算得出，而另一组对应河道断面的汇水面积数据则来源于水利普查等相关资料，其源头是人工在地形图上勾绘计算得到，这种就会存在人为误差，在小的汇水区域可能误差较小，但是在汇水面积很大，比如超过1 000 km2的断面，河流水系错综复杂，地形图上的山脊线和山谷线不易分辨，这样就可能会产生较大的误差。第三，Arcgis本身计算原理的误差甚至可能是错误。在提取河道断面的汇水面积过程中，其中有两个步骤是计算流向和填洼，关于其计算原理在前面有所介绍。填洼是如果某中心像元值均低于其周围的像元值，那么软件会认为这个中心像元值是一个错误值，从而会对该像元值进行填洼，使其像元值等于周围值的最小值，但这有时候并不符合现实情况，比如说自然或者人工形成的洼地，这是客观存在的，但是软件会将其“补平”；另外一个步骤是计算流向，Arcgis虽然提供了三种不同的流向计算方法，但是默认的(或者当您打开计算流向工具时该方法即为D8，不可更改)算法是D8算法，该算法认为任意一个像元值仅存在一个流向，从而会计算出某一点的唯一流向，但是现实中，某点的水流可能会有两个流向，只是每个流向分配的流量不一样而已，因此此种算法在某些情况下可能会造成很大的错误。第四，工具参数引起的误差。在使用Arcgis自带的工具“捕捉倾泻点”时，对于搜索距离的控制也是相当的关键，必要的时候还需要打开相应的矢量河流数据人工校对软件捕捉的结果是否准确，尤其是在河流水系复杂的支流汇口断面，张潇戈[21]等的研究也表明用DEM提取的河网汇口和真实河网汇口存在一定的偏差。

表3 计算方法在其他研究区的验证Tab.3 Validation of python-calculation in other study area

表4 研究方法误差来源分析Tab.4 Analysis of error sources of python-calculation

3 结 语

本文通过python调用Arcgis相应的工具来批量计算芙蓉江流域50 km2以上河流的河道断面的汇水面积，河道断面类型包含水库、支流汇口以及水文站等3种控制断面，分布广泛，代表性强，其主要结论如下。

(1)计算结果准确性高。计算的河道断面的汇水面积绝对误差Δ中最小的仅为0.08 km2，相对误差m最小的仅为0.09%，决定系数R2高达0.99，相对均方根误差RRMSE仅有2.74%，相对分析误差RPD为49.01，远大于2，偏差Bias为2.88，表明本方法不管汇水面积大小都能获得很好的计算结果。

(2)计算速度快，可极大地节省生产时间。本文计算的河道断面均位于贵州省遵义市境内的芙蓉江流域，根据已有的数据资料，为了方便计算结果的验证和比较，本文选取的河道断面共有37个，通过计算机计算汇水面积每分钟可以计算7个，大约需要6 min，而通过人工在地形图上勾绘，一个断面勾绘汇水面积大约需要25 min，38个断面共需要950 min，可见，通过python批量处理的时间仅为人工的0.63%，这在生产实践中极大地节约工作时间，提高工作效率，赢得时间成本。

(3)为河流控制断面的汇水面积勾绘提供一种新的计算方法。采用python调用Arcgis批量计算汇水面积，不仅可以极大地提高工作效率，还可以方便快速地查看和计算其他特征参数，为控制断面流量的计算和坝址库容曲线的计算分析奠定基础。