基于DEM的空间流速场构建方法及应用

谢志刚，石 朋,2，纪小敏，赵兰兰，瞿思敏，崔彦萍

(1.河海大学水文水资源学院，南京 210098；2.河海大学 水文水资源与水利工程科学国家重点实验室，南京 210098；3.江苏省水文水资源勘测局，南京 210029；4.水利部信息中心，北京 100053)

已有研究表明[1,2]，流域汇流可被视为净雨在地貌扩散和水动力扩散综合作用下的汇集过程。Redriguez-Iturbe和Gupta等[3,4]提出假定速度项表现为指数分布，则其线性反应函数为时间指数分布，由此流域地貌瞬时单位线可以解释为一滴雨抵达流域上任何处，其汇流时间的概率密度函数。芮孝芳等[5]研究发现，确定等待时间概率密度函数的实质就是确定水质点流速的概率密度函数，故把推求地貌瞬时单位线问题转化为确定流速的概率密度函数问题。

计算流速的方法有很多，主要的流速计算公式包括曼宁公式、梅德门特公式、谢才公式等[6]。1993年，梅德门特教授[7]基于栅格假定了一个时不变的空间流速场；1995年，Muzik[8]在梅德门特教授的基础上，结合曼宁公式与连续方程进行流速计算，并利用GIS确定流速场；2003年，Assefa M.Melesse等[9]人总结出结合曼宁公式与运动波近似方程来估算坡面流速，联合连续方程与曼宁公式求解河道流速，以此构建流速场；2006年，石朋等[10,11]采用曼宁公式计算栅格单元流速并以此建立了一个在空间上变化而时间上不变的空间流速场，以沿渡河流域为研究对象，取得了较好的计算效果。孔凡哲等[12]区分坡地单元与河道网格单元，以SCS公式计算坡地流速，以梅德门特公式计算河道流速，由此计算出整个流域的流速分布。考虑到在一次降雨径流过程中，雨强对流速的影响是不容忽视的[13]，本文拟通过将雨强纳入考虑范围并改进原有的曼宁公式，计算流域各网格内流速，并应用于定安河流域，以提高洪水模拟过程的精度。

1 研究思路

1.1 流速计算方法

流域下垫面条件存在空间差异性，不同地形地貌条件下流域各处的流速大小与方向各不相同。流速大小的计算是本文的核心所在。本文主要选用了两种流速计算方法，分别是：传统曼宁公式及考虑雨强的曼宁公式。

1.1.1 曼宁公式

传统的曼宁公式如下：

(1)

式中：n为地表粗糙系数，可以通过土地利用情况获得；S为地表坡度，可以通过DEM获得；R为水力半径，可以近似以水深代替。

由已有研究可知[14]，水深的计算可以建立水深与汇水面积的函数关系：

H=φpAψp

(2)

式中：H为具有超越概率p的栅格平均水深；A为栅格上游流域汇水面积；φp为网格常数，与洪水频率有关；ψp为几何常数，与洪水频率有关。

1.1.2 考虑雨强的曼宁公式

实际洪水过程中存在降雨过程，不同的降雨强度对汇流速度也存在影响，因此采用考虑雨强的曼宁公式进行流速计算，可提高计算精度。计算公式如下：

(3)

式中：i为雨强大小的无量纲因子；b为经验系数。

雨强作为影响流速大小的影响因子纳入曼宁公式中，如公式(3)所示。在GIS中利用克里金插值法得出流域雨强栅格图，即可知每个栅格内的雨强大小，再由已知的糙率系数、坡度和水力半径通过公式(3)可求出受雨强影响的栅格流速。由于整个洪水过程中存在降雨空间分布不均匀的情况[15]，不同地区雨强大小不同，并且在同一时段可能只有部分地区降雨，因此要对雨强栅格图进行0处理，即利用栅格计算器把所有0栅格替换成1，在无降雨地区i取值1，等价于直接用曼宁公式计算。

1.2 汇流时间计算

通过上述计算所得栅格流速，依据栅格大小可以计算出水流在栅格内的滞留时间。计算公式如下：

Δτ=L/v

(4)

(5)

式中：Δτ为栅格内滞留时间；L为栅格的边长。

公式(4)表示的是水流方向同栅格边线平行的情况，公式(5)表示的是水流方向同栅格对角线平行的情况。在GIS中通过D8算法可以确定水流方向，即确定每个栅格内水流方向，通过公式(4)和公式(5)计算各栅格内滞留时间，在GIS中沿水流方向累积各栅格到达流域出口的时间即可确定各栅格汇流累积时间，计算公式如下：

(6)

式中：τ为某栅格到达流域出口的累积汇流时间；m为某栅格径流路径上的栅格数。

1.3 单位线的确定

通过以上两种流速计算方法在GIS中可利用栅格计算器计算出每个栅格的流速大小，以此可得出两种不同的空间流速场，利用空间流速场以公式(4)、(5)可以计算得到每个栅格的滞留时间，通过公式(6)沿水流方向得出各栅格汇流累积时间，在GIS中利用重分类模块统计时段内通过流域出口的栅格数，即得到时段内出流面积。以汇流时间为横坐标，以时段内出流面积与总面积比值为纵坐标，可得到流域汇流时间的概率密度分布。根据径流过程形成的“粒子学说”[2-4]，即水质点在弱相互作用情况下的汇流时间概率密度分布函数等价于地貌瞬时单位线，可知以上所得汇流时间概率密度分布即为地貌瞬时单位线，将此应用于实例当中，对比径流模拟计算结果，分析两种流速计算方法存在的问题与优点。

2 实例应用

2.1 流域资料处理

本文所选流域是万泉河水系内的定安河流域，流域面积1 333 km2，内设长征、大墩、合罗、加报、罗担、木色、琼中、石古、思河、乌坡和乌石11个雨量站，具有完善的降雨资料。流域地处热带季风气候区，常年雨量充沛，多年平均降雨量约1 639 mm。流域是万泉河水系的一级支流，受水利工程影响较小，具有多年连续的降雨径流资料，适合本文模拟研究计算。

本文基于流域DEM(90×90)进行空间流速场的构建，推导出相应的地貌瞬时单位线，并应用于降雨径流模拟分析。流域坡度分布如图1所示。

图1 流域坡度分布Fig.1 Slope distribution of River Basin

根据流域土地利用情况确定了糙率系数n的栅格数据图；利用GIS中水文分析模块可计算各栅格上游汇水面积，以此根据公式(2)可求出流域水深栅格数据图；通过克里金插值法得出流域时段雨强栅格数据图。利用GIS中栅格计算器，以上述得出的栅格数据通过公式(1)～(3)计算每个栅格的流速，分别构建出以曼宁公式计算出的空间流速场和以考虑雨强的曼宁公式计算出的空间流速场，雨强曼宁公式以20010825次洪为例。图2(a)为曼宁公式构建出空间上变化而时不变流速场，图2(b)为2001年8月25日次洪降雨第2个时段雨强的空间流速场栅格图。

图2 计算所得空间流速场Fig.2 Calculated spatial velocity field

根据汇流时间计算公式(4)～(6)可以计算出相应的汇流累积时间栅格分布图，图3(a)为曼宁公式计算出的累积时间栅格图，图3(b)为2001年8月25日次洪降雨第2个时段雨强的计算出的累积时间栅格图。

图3为沿水流方向到达流域出口的汇流累积时间分布图，距离汇流出口越远所需汇流时间越长，与实际情况相符合，据此统计汇流时间的概率密度分布函数[16]，得出地貌瞬时单位线。雨强曼宁公式由于各时段雨强不同所得出的单位线也不同，故采用时段出流方式进行汇流计算。

图3 汇流累积时间分布Fig.3 Accumulation time distribution of confluence

2.2 降雨径流模拟计算

利用上述所得地貌瞬时单位线转换成1 h时段单位线，应用于定安河流域进行降雨径流模拟分析，选取该流域2000-2013年中的7场洪水进行模拟计算，流域每年洪涝灾害多发生于5-11月之间，强降雨主要集中在7-10月，考虑雨强对流速的影响，故选取此期间7场洪水进行模拟计算，所得计算结果如表1所示。经过对比分析，由曼宁公式计算值相比较于考虑雨强的曼宁公式计算值，7场洪水平均径流相对误差由11.52%下降到5.93%，平均洪峰相对误差由6.32%下降到3.85%，平均确定性系数由0.761上升到0.840。由此可知考虑雨强大小影响的曼宁公式计算所得结果更符合实际情况，其径流深、洪峰流量、确定性系数模拟效果更优[17]。图4为20010825和20071011两场洪水由曼宁公式计算的洪水过程线和考虑雨强的曼宁公式计算的洪水过程性与实测洪水过程线的对比图。

表1 定安河径流计算结果Tab.1 Runoff calculation results of the Ding An River

注：此平均值为绝对值平均值。

图4 洪水模拟计算对比Fig.4 Comparison of flood simulation

3 结 语

本文核心是在流域流速计算方法上，而传统地貌瞬时单位线的确定是建立在流速空间分布均匀的情况下，实际流域中流速大小受地形地貌影响以及水力条件影响较大，故本文提出了考虑雨强影响的曼宁公式计算流域流速，建立了整个流域的空间流速场。曼宁公式充分考虑了土地利用情况、坡度以及水深对流速的影响，由此可以确定一个空间上变化而时不变的空间流速场，但是实际洪水过程中存在降雨空间分布不均匀性，不同区域雨强不同亦对流速产生重大影响，故此加入雨强这一变量因素，以此改进曼宁公式求得流域空间流速场，应用于定安河流域，取得了比较好的模拟效果，故此可推广应用于缺资料地区的汇流计算过程中。

本文所取雨强为1 h时段雨强，未考虑不同时段大小雨强对流速的影响，后期重点研究不同时段大小雨强对流速计算精度的影响，提高径流模拟效果。本文所取栅格大小是90 m×90 m，未考虑不同分辨率大小情况下坡度值精度不同，由此产生流速计算误差，故在后续研究过程中会进一步考虑分辨率对流速计算产生的影响。