栅格流向算法在山丘区水系提取中的应用评价

2019-10-23崔寅鹤吴鹏飞汪银奎刘金涛

水力发电 2019年7期

崔寅鹤，吴鹏飞，汪银奎，刘金涛，3

(1.河海大学水文水资源学院，江苏南京210098； 2.西藏自治区水文水资源勘测局，西藏拉萨850000； 3.河海大学水文水资源与水利工程科学国家重点实验室，江苏南京210098)

0 引言

地表流径描述了水、土壤、污染物等受重力影响在陆地表面的运移路线[1]，在水文建模[2]、土壤侵蚀[3]等研究领域应用广泛。除少数基于不规则三角网[4]、等高线[5]等数据的算法外，目前主流的流径提取算法普遍以栅格数字高程模型(Digital elevation model，DEM)作为数据基础[6- 7]。使用DEM提取流径的算法普遍包括填洼、流向计算、汇流累积3个步骤[8]。在山丘区，流径提取的结果主要受流向算法影响。流向算法分为单流向算法和多流向算法两类[9]。以3×3网格为基础，单流向算法将本栅格产生的水流和上游栅格汇流至本栅格的水流作为一个整体，认为其只能流向8个与本栅格相邻栅格中的1个；多流向算法则将水流按比例分配给多个相邻栅格。目前，应用较多的单流向算法包括D8算法[10]、D8-LTD算法[11]、GRASS软件rwflood模块使用的快速算法[12]等，多流向算法包括MFD算法[13]、Dinf算法[14]、MDinf算法[15]等。其中，Dinf算法虽然仅计算得到一个流向，但在进行汇流时需将流量按比例分配给两个相邻单元，故此处将其定义为多流向算法。

表2 流向算法原理

目前对各种流向算法的对比分析大多聚焦于D8、MFD、Dinf等经典算法[16- 17]，对较新的rwflood、MDinf等算法的功能评价较少。本文选取3处典型高山、丘陵地形，使用刘学军等[17]提出的定量化对比方法分析了6种流向算法(D8、D8-LTD、rwflood、MFD、Dinf、MDinf)的差异性，并对比了这些算法对真实河道的还原效果以及计算速率，以期为相关研究的流向算法选取提供参考。

1 材料与方法

1.1 研究材料

为了更好地研究流向算法在不同山丘区的应用效果，本文分别在青藏高原、东南丘陵、黄土高原选取了姜湾研究区、夺底沟研究区和杜家湾研究区(见图1)。其中，夺底沟研究区坡度较大，平均坡度25.58°；姜湾研究区山地坡度居中，平均坡度20.23°；杜家湾研究区为黄土丘陵，山地坡度最小，平均坡度13.23°。

图1 研究区位置及DEM数据

本文所选用的DEM数据源为SRTM V4.1 30 m分辨率数据。SRTM数据通过航天飞机搭载雷达进行测绘， V4.1版本采用了新型的插值方法，让数据的准确度有了较大的提升。本文使用的全部数据下载自美国地质调查局网站(http://earthexplorer. usgs.gov/)。

1.2 流向算法

流向算法旨在确定水流流动、土壤和污染物输移路径，体现在DEM数据中就是要解决各栅格向下游栅格分配流量的问题。基于这个问题可以将流向计算分为单流向算法和多流向算法两大类。

1.3 研究方法

首先对各研究区DEM数据进行填洼(rwflood算法不需要填洼)，利用填洼后数据进行流向计算和回流累积；再根据汇流累积结果设定最大流向的1/100为阈值提取水系。本文使用的填洼算法是LIU等[18]提出的算法，对平坦地区赋予1×10-5m的梯度，汇流算法为Su等提出的基于树状图的新算法[19]。这两种算法均为目前计算速度较快的算法。

本文采用了刘学军等提出的相似度量化分析方法，将汇水面积(Total Catchment Area， TCA)选为分析对象；并且将TCA进行了对数变换以避免由TCA分布不均匀造成误差的现象[17]。此外，本文还对各算法在实际应用方面的效果进行了评估。

1.3.1累计频率分布

在分析和表示地形属性分布特征时，频率和频率分布图是一种统计分析方法。TCA的数量表现和算法的不同通常会导致使用不同算法的TCA累计频率出现差异。这种差异可以直观地反映为频率和频率分布曲线形状的不同。从差异特征中便可分析得出算法的机理以及TCA的空间分布。

1.3.2相对差系数

应用不同机理的径流算法所得到的DEM在相同单元格内TCA值存在差异。为了便于比较这种差异，故将其以相对差系数α来表示[17]。α的定义为

(1)

表3 各研究区TCA分布

当α≤0时，则两种算法之间没有可比性；当0<α<1时，两种算法存在差异且α值越小差异越大；当α=1时，两种算法之间没有差异。由式(1)可看出，相对差系数是一种可将不同算法进行两两相比的数值指标，而在计算过程中需要假定其中一种算法作为基准算法；因此，α的大小取决于基准算法的选择。由于基准算法的不同，每一对进行比较的算法会产生2个α。

1.3.3Google Earth实际效果对比

使用不同算法得到的研究区河网图像在保证河道连续，没有异常汇流现象的同时还要考虑其结果与实际地形的吻合程度。Google Earth作为一款虚拟地球软件，可以清晰地看到全球各地的卫星照片。通过ArcGIS将栅格中的最大值除以100作为阈值，提取大于阈值的栅格；再将值为1的点过滤，结果即为河道；之后再将河道数据转换为用Google Earth可以打开的格式，便可将河网置入卫星照片中以比较其与实际地形的匹配程度。

1.3.4运行时间对比

不同算法的实际应用价值，其运行时间也是评价的一个重要方面。因此，本文对每种算法在填洼(r.watershed算法中不包含填洼步骤)、流向计算、汇流计算3个步骤的用时也做了统计和比较，以便综合实际效果等方面对每种算法做一个应用评价。为了减小误差，采用运行3次求取平均值的方法。

2 结果与分析

2.1 运行成果对比

D8、D8-LTD、rwflood、MFD、Dinf和MDinf算法在3个研究区的DEM上计算得到的TCA图(见表3)，TCA值越大，灰度越暗。由对比分析可看出，同一研究区，单流向算法较多流向算法水流路径更加细密，D8和D8-LTD两种算法最为细致且差异较小，多流向算法的水流路径则比较宽。之所以会有这样的差异是因为两类算法在向下游分配水流时的方式不同，多流向算法向下游多个甚至所有的栅格进行分配，使区域内多数栅格获得的流量较单流向算法多，集水面积的差别减小。对于多流向算法，由于Dinf和MDinf算法分配给下游最陡方向的流量比MFD多，故灰度图颜色相对较淡。另外，不同算法之间提取主要水系差异性较小，但在各栅格流量数值上存在较大差异，可见在应用过程中，流向算法对于注重整体研究区特征的研究影响不大，对于注重各栅格TCA数值的研究则存在较大影响。

2.2 累积频率分布分析

图2为3个研究区的TCA累积频率分布图。在ln(A)≤10.5时，6条曲线逐渐趋于重合，当ln(A)≥10.5时，6条曲线几乎完全重合。

图2 3个研究区ln(A)累积频率分布

在图2中，单流向算法(D8、D8-LTD、rwflood)与多流向算法(MFD、Dinf、MDinf)的图线存在比较明显的差别，易于区分。单流向算法中汇水面积较小的栅格所占比例远大于多流向算法。这与2.1节相同，同样是由单流向算法和多流向算法流量分配策略的差异导致的。

MFD在相同频率下的汇水面积最大。这与其对流量的分配方式有关。所有的下游栅格都可按比例分配到一定的上游水量；从而使得汇水面积较其他算法更大。Dinf和MDinf在图线的走向上基本相同，数值介于单流向算法与MFD算法之间。另外，从图2中也可看出夺底沟研究区6种算法之间的差异略大于另外两研究区，这是由于在实际地形中，夺底沟研究区下游存在部分平坦地区，多流向算法使水流发散的特征被显著增强。

2.3 相对差系数分析

为了分析不同地形单元[17]上6种算法所计算的TCA之间存在的差异，将地形分为坡面区域(ln(A)≤10.5)和汇流区域(ln(A)>10.5)两部分，分别对全研究区、坡面区域、汇流区域计算相对差系数(见表4)。表3中3个数据的含义分别为全研究区计算的相对差系数，山坡区域相对差系数，汇流区域相对差系数。

从表4中可以看出，除rwflood算法外，两种单流向算法或者两种多流向算法进行比较时计算得到的α值较大，而单流向算法与多流向算法相互比较时得到的α值较小。由于Dinf和MDinf算法在流量分配过程中分给最陡坡度的流量极大，因此与单流向算法相似度较高。D8-LTD算法是对D8算法基于全局偏差进行校正，大多数栅格流向与D8相同，所以相似度高。另外，相互比较的两种算法在基准算法不同的情况下得到的山坡区域α值差异不大。

Dinf算法由于实际仅计算一个流向，具备单流向算法的部分优势，在与单流向算法相比时α值较其他多流向算法高。rwflood算法由于流向判断基于高程而非物理意义更合理的坡度，得到的所有相对差系数均较小，与其他算法差异大；因而，在追求精度时的流向计算中应该尽量避免使用rwflood算法。此外，MFD算法在多流向算法中α值整体较小，即与其他算法差异较大，这是由于MFD算法计算中更加强调水流弥散的影响。以上分析结论验证了刘学军等[17]的成果。

2.4 Google Earth实际效果对比

表5展示了将应用6种算法提取出的全研究区的河道导入进Google Earth软件后的主要河道分布图，以及将筛选数据的阈值适当调小后得到的局部山区的水系分布图。单流向算法获取的河道利用ArcGIS转成了线，多流向算法由于可以反应河道宽度故保留了原有的栅格形式。由于不同算法提取出的河道有所重叠，故图线不能全部显示。

从表5中全研究区图中可以看出，除rwflood算法存在提取的主河道过于平行、顺直，与实际河道偏离较多外，其他5种算法在河谷中的差异较小，提取出的主河道基本重合且与实际河道吻合程度较高，3种多流向算法反映的河道宽度较准确。在夺底沟研究区存在平坦地区，各种算法的结果存在比较大的差异，且准确性不高，这是由于平坦地区受填洼算法选取、人类活动等影响，河道的还原难度大。对于局部山坡，单流向算法得到的水系存在较多的过直河道，rwflood算法的过直河道最多，与其他算法差异最大，甚至出现了多处河道翻越山脊的错误情形。

表4 各研究区6种算法相对差系数

表5 Google Earth全研究区及山区局部导出图

整体而言，rwflood算法对河道的还原效果最差，多流向算法反映真实河道宽度、走向的效果优于单流向算法，D8-LTD算法较D8算法更接近于真实河道。相同阈值下，单流向算法得到的河长大于多流向算法，MFD算法因水流发散程度最严重得到的河长最短。在TCA值勉强达到阈值的源头栅格，一旦出现水流发散的情况，存在该栅格的所有下游栅格TCA值均低于阈值的可能，因此多流向算法提取的可视化河道在水系上游存在断流现象。