基于DEM建立小型水库流域R-V地貌瞬时单位线之初探

刘 芝

(四川省资阳市雁江区水务局，四川 资阳，641300)

1 概述

在水库洪水预报中，传统的预报方法是绘制水库洪水预报综合查算图[1]，但该查算图需要根据本流域大量的实测降雨流量资料进行配线，可移用性较差，必须做到一库一图，而对于实测资料缺乏的小型水库，只能借用邻近流域特性参数或地区综合经验参数绘制查算图，预报结果的可信度不高。随着社会经济的发展和社会管理水平的提高，尤其在“智慧地球”、“智慧城市”、“智慧水利”的发展框架下，国家在逐渐筑牢大江大河防洪体系的同时，对水库洪水预报能力的要求也在不断提高。研究和揭示水库洪水与流域地形地貌的关系，建立揭示流域水文响应与地形地貌关系的地貌瞬时单位线(GIUH)及具有物理基础的分布式水文模型，是破解小型水库洪水预报难题的必然出路。随着计算机技术和地理信息技术的发展，尤其是数字高程模型(DEM)的广泛应用，提取流域河网及地形地貌信息成为可便捷实现的过程。地貌瞬时单位线(GIUH)和分布式水文模型[2]在小型水库洪水预报中的推广和应用，成为未来的必然趋势。

2 R-V GIUH基本理论

Rodriguez-Iturbe和Valdes于1979年首先提出了地貌瞬时单位线(GIUH)的概念，并将瞬时均匀降落在流域上的净雨看作由无数的水质点组成，用概率论的方法来研究流域上水质点的运动，当组成净雨的水质点间呈弱相互作用时，流域瞬时单位线等价于水质点持留时间的概率密度函数，这就是经典的R-V GIUH理论。

2.1 河流分级

R-V GIUH理论是基于Horton-Strahler河流分级法提出的。Horton-Strahler河流分级的具体方法为：将所有河网弧段中没有支流的河网弧段定为第1级；两个1级河网弧段汇成的河网弧段定为第2级；两个2级河网弧段汇成的河网弧段定为第3级；如此下去，当且仅当同级别的两条河网弧段汇成的河网弧段级别增加一级，直到河网出水口。

2.2 Horton地貌参数

R-V GIUH理论将瞬时单位线参数同流域地貌参数结合，提供了流域响应函数的地貌学解释。自然界的水系在发育过程中遵循一定的规律。由于受地质构造和自然环境的影响，河网形态在平面上表现为有规律的排列并在不同尺度上具有一定的自相似性。Horton通过大量研究发现[3]，任何一个流域内，各级河道的数目、平均长度、平均集水面积、平均纵比降与河道的Horton-Strahler河流分级级别之间存在几何级数关系，并据此提出了Horton地貌定律，即河数定律、河长定律、面积定律、比降定律。而Horton地貌定律中的河数率RB、河长率RL、面积率RA是推求R-V地貌瞬时单位线时重要的地貌参数。

Horton地貌参数(河数率RB、河长率RL、面积率RA)的计算公式如下：

RB=Ni-1/Nii=2，3，…，Ω

一般来说，不同级别的河流的河数率、河长率、面积率并不相同，应用时可取各级河流的平均值。且大量研究成果显示：天然河网的Horton地貌参数存在一个较窄的变化范围，河数率范围为3～5，河长率范围为1.5～3.5，面积率范围为3～6。

2.3 R-V GIUH基本公式推导

基于Horton-Strahler河流分级方法，把水质点到达的某一级河流或坡面称作状态xi，把水质点从注入点到出口断面顺序经过的各状态称作路径si，则水质点沿路径si流至出口断面的汇流时间为Tsi=Tx1+Tx2+…+TxΩ。根据流域瞬时单位线等价于水质点持留时间的概率密度函数，可导出流域地貌瞬时单位线基本表达式[4]：

μ(t)=fB(t)=∑siSfx1(t)×fx2(t)×…×fxΩ(t)·p(si)

(1)

式中：fx1(t)，fx2(t)，…，fxΩ(t)表示水质点在当前状态下的等待时间概率密度函数；p(si)表示水质点取路径si的概率；S表示所有路径si的集合，有S={s1，s2，……}；Ω表示河流最高级别；符号×表示卷积相乘。

式(1)中的路径概率p(si)有：

p(si)=θx1px1x2px2x3…pxΩ-1xΩ

(2)

式中：θxi为水质点处于初始状态的概率，简称初始概率；pxixj(i

(3)

(4)

式(1)中等待时间概率密度函数fxi(t)的确定，实质上就是确定水质点流速的概率密度函数。常用下列单参数指数型函数来进行拟合：

(5)

式中：Ki为水质点在状态xi的平均持留时间，等于状态的空间尺度与水质点在该状态的平均流速的比值，可用地形地貌学方法确定。

将R-V GIUH基本表达式导出为数学计算公式的过程比较复杂，且当河流最高级别Ω不同时，导出的计算公式形式亦不相同，使用时需要分河网级别予以导出。R-V GIUH理论原作者给出了三级河网的R-V GIUH计算公式，由文康[6]等于1988年给出了四级河网和五级河网的R-V GIUH计算公式。

2.4 流域平均流速的推求

(6)

式中：v为流域平均流速；L为主河道长度；τ为流域汇流时间，可用下列公式推求：

(7)

(8)

式中：ψ为洪峰径流系数；F为流域集水面积；L为主河道长度；J为河道平均比降；S为暴雨雨力；n为暴雨公式指数；m为汇流参数。以上参数均能利用《四川省水文手册》进行查算。

3 基于DEM提取流域河网和Horton地貌参数

基于DEM提取流域河网和地貌参数可在ArcGIS软件中实现。ArcGIS软件是美国环境系统研究所ESRI开发的一款无缝扩展GIS软件，具有强大的地理信息分析处理能力。ArcGIS软件中的Hydrology工具集能够方便地提取河网结构等大量水文信息。

3.1 DEM的获取

数字高程模型(DEM)是表示地面高程随空间变化的数学模型，是对地形地貌的一种离散的数字表达。DEM的获取方式主要有三种：一是通过互联网上既有的DEM数据库直接下载，常用的DEM产品有ASTERGDEM、SRTM、GLS2005、TANDEM等，分辨率多为90m和30m；二是通过现有的地形等高线图生成DEM，即将现有地形等高线图数字矢量化后，进行TIN建模，进而生成DEM；三是根据实测数据源生成DEM，DEM数据源的获取方式包括地面测量、数字摄影、激光雷达、立体遥感、GPS等。普通的研究和应用采用实测数据源生成DEM的可能性不高，通常采用第一、二种DEM获取方式，且以互联网下载方式最为便捷。

3.2 DEM流域位置定位

互联网下载获得的DEM是分区块的，在区块DEM数据上定位到具体研究区域的位置，需要借助于行政边界矢量数据。在ArcGIS软件中打开DEM数据，并将下载到的行政边界矢量数据加载到DEM上(需注意坐标系统的一致性)，然后使用ArcToolbox中的Extract by Mask工具进行剪切，即得到行政边界区域内的DEM数据。将对应行政区域的水系图通过空间校正叠加到DEM数据上，可定位出具体流域的位置，再沿集水区域边界剪切，便可得到具体流域范围的DEM数据。

3.3 DEM洼地填充处理

受DEM分辨率的限制，以及DEM生产过程中的系统误差，原始DEM表面存在着一些凹陷或尖峰的区域，这些是DEM数据的缺陷，会导致不合理的甚至是错误的水流方向计算结果。但并非所有的洼地都是数据误差造成的，也有可能是真实存在的地形。故需要进行洼地计算，计算出洼地的位置、区域、深度，以此判断哪些洼地是数据误差，哪些洼地是真实地形。具体计算步骤如下：①用Hydrology工具集下的FlowDirection工具计算出水流方向；②依据水流方向，用Sink工具计算出洼地位置；③用Watershed工具计算出洼地贡献区域；④用Zonal工具集下的ZonalStatistic工具计算出洼地贡献区域最低高程；⑤用Zonal Fill工具计算出洼地贡献区域出口最低高程；⑥用MapAlgebra工具集下的Raster Calculator工具计算出洼地深度，洼地深度=洼地贡献区域出口最低高程-洼地贡献区域最低高程。依据计算出的洼地深度设置合理的洼地填充阈值，用Hydrology工具集下的Fill工具对洼地进行填平处理，得到接近于真实地形的无洼地DEM。

3.4 河网提取

河网是水在地形地貌上的产物，同时水的汇集与流动又受到河网的约束。河网特征反映一个流域的综合水文特征。常用的河网提取方法为连续河网生成法，即采用D8法按最陡坡度确定DEM中每个网格单元的水流方向，根据水流方向计算出每个网格单元的上游集水面积(即汇流累积量)，设置一个集水面积阈值，以该阈值作为河道认定的门槛，由汇流累积量大于该阈值的网格单元组成河网。河网提取具体步骤如下：①用Hydrology工具集下的Flow Direction工具计算出无洼地DEM的水流方向；②依据无洼地DEM水流方向，用Flow Accumulation工具计算出汇流累积量；③设置适当的集水面积阈值，用Map Algebra工具集下的RasterCalculator工具生成栅格河网；④用Hydrology工具集下的Stream to Feature工具将栅格河网转换为矢量河网。

3.5 流域划分

流域又称集水区域，是指由分水线所包围的河流集水区，流经其中的水流从一个公共的出水口流出。流域的边界可依据DEM的数字地形分析确定。流域划分具体步骤如下：①依据提取出的流域河网，用Hydrology工具集下的Stream Link工具计算出河网结点；②依据无洼地DEM水流方向，用Basin工具划分出流域盆地；③依据河网结点计算结果，用Snap Pour Point工具计算出汇水区出水口；④用Watershed工具生成集水流域；⑤用Conversion Tools工具集下的Raster to Polygon工具对栅格集水流域进行矢量化，得到矢量集水流域。

3.6 Horton地貌参数的提取

Horton地貌定律是以Horton-Strahler河流分级为基础的。河网的分级包含重要的水文信息，可用以研究水流的运行和汇流模式等。河网分级具体步骤为：①根据提取的栅格河网，用Hydrology工具集下的Stream Order工具，选择Strahler分级法，完成河网的分级；②用Stream to Feature工具对栅格分级河网进行矢量化，得到矢量分级河网。

Horton地貌参数的提取步骤具体如下：①打开矢量分级河网的属性列表，添加列，用Calculate Geometry工具计算出河段长度，导出含河段长度的属性列表，便可获得每一级别河段的条数和长度等信息；②打开矢量集水流域的属性列表，添加列，用Calculate Geometry工具计算出集水流域面积，导出含集水流域面积的属性列表，便可获得每一级河流的集水面积等信息；③依据Horton地貌参数计算公式，计算得到河数率RB、河长率RL、面积率RA等地貌参数值。

4 应用实例

4.1 研究区域概况

雁江区位于四川盆地中部，全区共有小型水库78座，水库流域面积介于0.1km2～450km2之间，普遍没有降雨和流量观测站，受预报方法和管理水平的限制，其水库洪水预报工作基本上是空白的。本次研究以高板桥水库流域作为研究对象，进行流域河网和Horton地貌参数的提取，并建立水库流域的R-V地貌瞬时单位线。高板桥水库位于雁江区东峰镇境内，属沱江水系，地理坐标为E104°53′18.3″、N30°3′50.3″，总库容402万m3，坝址以上集雨面积33.77km2，主河道长13.355km，河道平均坡比降1.317‰，是一座以防洪为主、兼顾农业灌溉等综合效益的小(1)型水库。

4.2 流域河网和地貌参数的提取

研究使用的DEM为从地理空间数据云免费下载的ASTER GDEM V2 30m分辨率DEM，雁江区的行政边界矢量数据为从全国地理信息资源目录服务系统免费下载获得，流域定位所用的雁江区水系图比例尺为1∶85000。通过图层叠加和剪切，可分别获得雁江区和高板桥水库流域的DEM，见图1。

(a)雁江区DEM (b)高板桥水库流域DEM

按照洼地计算方法对原始DEM进行洼地计算，得到洼地的位置、区域和深度。经计算，高板桥水库流域共有洼地21处，洼地深度最大值25m，最小值1m。由于本研究使用的原始DEM分辨率只有30m，故结合实际地形考虑，认为以上洼地均是由于数据误差产生的伪洼地，故将所有洼地予以填平处理，得到研究区域的无洼地DEM。

设置适当的集水面积阈值，利用ArcGIS软件中的Hydrology工具集进行流域河网提取，得到高板桥水库流域的流域河网，见图2。

(a)雁江区流域河网 (b)高板桥水库流域河网

对提取出的高板桥水库流域河网进行Strahler分级运算，获得具有河流分级属性的流域河网数据，并以此计算出各级河流的条数、长度、集水面积。根据Horton地貌参数计算公式，可计算得到高板桥水库流域的河数率、河长率、面积率等地貌参数，见表1。

表1 高板桥水库流域Horton地貌参数

可见，基于DEM提取出的高板桥水库流域河网共分成4级，总集水面积42.0783km2，与从万分之一地形图上量算到的水库坝址以上集雨面积33.77km2存在一定差异。高板桥水库流域的河数率为4.3667、河长率为2.3821、面积率为4.9665，均在自然河流Horton地貌参数的取值范围内，可用于进行流域R-V地貌瞬时单位线的推求。

4.3 R-V地貌瞬时单位线推求

根据《四川省水文手册》用推理公式求得高板桥水库流域的平均汇流时间约为8h，进而求得高板桥水库流域平均流速约1.67m/s。

将提取的高板桥水库流域地貌参数和计算的流域平均流速带入四级河网流域地貌瞬时单位线公式进行计算，发现初始概率θx4<0，这是不合理的，此处采用文康[6]等推荐的处理方式，即当不满足RB/RA<0.8这个约束条件，初始概率出现负值时，直接以不包括低一级河道集水面积在内的各级河网集水面积与流域总面积之比表示初始概率。

最后，可得出高板桥水库流域的地貌瞬时单位线为：

μ(t)=-0.1671e-9.2492t-20.5336e-4.5271t+6.5354e-0.9442t+17.2233e-3.058t

将上述地貌瞬时单位线与适当的产流模型结合，就能实现高板桥水库流域的水库洪水预报。

5 结论和讨论

(1)本文以雁江区高板桥水库流域为实例，基于DEM提取出流域河网和Horton地貌参数，并成功建立了水库流域的R-V地貌瞬时单位线，为进一步探索小型水库流域分布式水文预报模型奠定了基础，对推动小型水库洪水预报方法研究具有积极意义。

(2)本研究采用中小流域暴雨洪水推理公式计算流域平均流速，计算公式中不仅包含下垫面特性参数，还充分考虑了降雨强度对流域平均流速的影响，且其参数值均可通过《四川省水文手册》查算，具有一定的实用和推广价值。

(3)本研究所使用的DEM为从互联网免费下载获得，分辨率仅为30m，且经过源数据插值处理，数据精度欠佳，以致提取出的流域参数与实际流域有一定出入。在实际的水库洪水预报应用中，为更准确的获取流域地形地貌信息，尚需采用精度更高的DEM产品。

(4)研究中发现，栅格河网生成时，设置不同的集水面积阈值，可得到不同密度的河网，主观随意性比较大。可尝试选取不同地区不同尺度的多个流域做专题研究，探寻河网生成阈值的固有特性，寻求不同区域不同尺度下的河网生成阈值常值。

(5)R-V GIUH理论未考虑河网的水动力特性空间变异化，即认为水质点在流域中的传播速度是固定不变的，然而实践证明水流在河网中的传播是一个非线性过程，且水流在从上游向下游传播的过程中，其流速缓慢增加。同时，R-V GIUH理论几乎完全忽略了坡面汇流过程，将坡面汇流时间忽略不计，这在坡面汇流滞时比重较大的小流域上，是无法真实模拟流域汇流过程的。故此，R-V地貌瞬时单位线用于小型水库流域汇流计算，是存在一定理论缺陷的，尚需探索相关方面的改进。