基于DEM的资料短缺地区设计洪水推求

张 婷，冯 平

(天津大学水利工程仿真与安全国家重点实验室，天津 300072)

1 前言

资料短缺地区的水文研究从20世纪90年代开始逐步被全球的水文工作者所重视［1］。国际水文科学协会(International Association of Hydrological Sciences，IAHS)于2003年启动了无资料(资料不足)流域(地区)的水文预报，即一个简称为PUB(prediction in ungauged basins)的国际水文计划，意在未来10年全面开展资料短缺流域的水文计算研究［2］。

目前国内外学者已经对PUB展开了积极的工作。李红霞等基于目前常用的几种区域化方法，提出了“综合相似法”来选择参证流域［3］。柴晓玲等研究了IHACRES模型在资料短缺地区径流模拟中的应用［4］。陈志明指出了地貌瞬时单位线采用地貌特征信息模拟资料短缺地区径流时存在的一些问题［5］。在国外，典型的研究主要有用分解法进行资料短缺流域的水文预报，它实际上是我国地区综合法的一个特例。

对于资料短缺地区的水文计算，可通过数字高程模型(digital elevation model，DEM)技术与水文模型耦合的方法，将详细的地理信息(如高程、地理条件、土壤类型、植被、土地利用等)用于模型的参数估算［6］。这些地貌特征信息的输入提高了资料短缺流域水文模拟模型的确定性。

2 流域概况

福州市位于福建省的东部、闽江中下游，总面积为12184 km2，海拔多为600～l 000 m，属于典型河口盆地地貌。全区地貌类型以山地、丘陵为主，占全区总面积的72．7%。福州地区年平均气温为19．1 ～20．1 ℃，年平均降水量为1200 ～1700 mm［7］。

文章的研究流域为福州市北部山区，东起磨洋河，西至过溪水库，面积约为44 km2，环绕福州江北中心城区。北部山区洪水通过过溪、新店溪、马沙溪等11条河道排入江北中心城区水系，山区洪水与城区洪水相遇后一起排入闽江。因此对北部山区进行设计洪水的估算对江北城区的防洪排涝规划管理至关重要。

根据水系分布特点将其划分为11个子流域，利用DEM图提取各子流域(见图1)，其中过溪、新店溪、登云溪上游有过溪、八一、登云3座水库对其水量进行调蓄。该流域有赤桥(1966—2008年)、叶洋(1970—2008年)、岭头(1970—2008年)等雨量站的降雨观测资料(地理位置见图1)和赤桥水文站的各时段最大降水量(1955—2008年)以及文山里水文站的蒸发资料(1979—2008年)可供使用，而只有赤桥水文站(1966—2008年)能提供流量资料，其位置在新店溪子流域北部，其他子流域均无流量资料。可见，福州市北部山区的设计洪水计算问题属于资料短缺地区的设计洪水计算。

图1 研究流域DEMFig．1 DEM of the research basin

3 方案设计

由于研究流域内径流资料短缺，雨量资料相对充足，故采用通过暴雨资料推求设计洪水的方法。本流域地处南方湿润地区，植被良好，雨量充沛，地下水丰富，符合蓄满产流机制条件，故可采用新安江三水源模型计算产流。具体研究方案如下:

3．1 设计暴雨计算

根据短历时最大降水量统计资料，采用Pearson-Ⅲ型曲线进行暴雨频率分析［8］。给出不同标准下(频率为2%、5%、10%和20%)的各时段暴雨设计值。然后选取典型暴雨，根据设计暴雨值和典型暴雨的时程分配推求不同频率下的设计暴雨过程。

3．2 产流计算方案

采用新安江三水源蓄满产流模型。由于赤桥流域与研究的各子流域地理位置相近，而且地形、地貌、植被、土壤、水文、地质条件等均非常相似，因此，模型参数可通过赤桥水文站的流量资料进行率定，再移植到其他子流域。

3．3 汇流计算方案

根据净雨在流域上汇流的途径不同，将汇流过程分成地面汇流、壤中流汇流和地下汇流3个部分。壤中流和地下径流分别采用不同的线性调蓄水库模拟其汇流过程。由于缺少径流资料，文章通过DEM推求等流时线，然后采用等流时线法模拟流域地面汇流过程。

4 设计暴雨推求

采用Pearson-Ⅲ型曲线，由北部山区短历时最大降水量统计资料，分别绘制其最大1 h、3 h、6 h、12 h和24 h降雨量频率分布曲线，如图2是最大24 h降雨量的频率分布曲线。然后，给出各时段相应频率为2%、5%、10%和20%的设计暴雨值(见表1)。

图2 北部山区最大24 h降雨量频率分布曲线Fig．2 The frequency distribution curve of the biggest 24 h rainfall in northern mountain area

选取1966年9月3日10:00至9月4日10:00的24 h暴雨为典型暴雨。选择原因为:该场降雨雨量大，强度也大，暴雨核心部分出现在后期，形成的洪水主峰出现较迟。根据设计暴雨值和典型暴雨的时程分配推求出不同频率下的设计暴雨过程。

5 模型计算

5．1 产流计算

新安江三水源产流模型主要由三部分组成:蓄满产流计算、流域水源划分和流域蒸散发计算［9］。

5．1．1 蓄满产流

设R是净雨量，降雨量为P，蒸发量为E，蓄水容量为Wm，降雨开始时的实际蓄水量为W，对于流域中某点而言，蓄满前，R=0;蓄满后:

表1 不同频率下的设计暴雨值Table 1 Design rainstorm in different frequencies

5．1．2 划分三水源

设RS、RSS、RG分别为地面径流、壤中流、地下径流，则有:

5．1．3 蒸散发计算

设E为流域蒸发量，EU、EL、ED分别为上层、下层、深层蒸发量，则有:

5．2 汇流计算

5．2．1 地面径流汇流

采用等流时线法。随着流域的地形地貌逐渐可用DEM来表示［10］，于是便可以通过流域的DEM而不是地图来自动生成等流时线分布图［11］。具体提取的方法如下［12］。

1)首先对原始DEM数据进行洼地填充，得到无洼地的 DEM［13］。

2)通过Arcgis软件或者Swat软件提取各个子流域。

3)以栅格为单位，通过Arcgis软件的Flow Direction工具提取水流方向，用Slope工具提取坡度。

4)通过经验公式来计算流经每个栅格的速度vi，如:

式(4)中，i为栅格序号;Si为某一水流出流方向上的坡度;a为参数，具有速度的量纲;b为幂指数，反映坡度大小对流速的影响。在利用式(4)来计算流经栅格的速度时，必须先合理地确定公式中的参数a和b。方法有两种，一种是通过水力学方法来直接测定，另外一种是通过水文学中的单位线理论来间接率定［14］。文章采用第二种方法，利用赤桥水文站资料率定得到参数 a=1．55、b=0．3，然后移用到其他子流域。

5)根据式(4)计算出速度图层后，利用栅格计算器计算该图层的倒数，即1/vi。各个栅格的汇流时间可以通过提取水流长度(Flow Length工具)，在“输入权重光栅”(Input weight raster)一栏中选择1/vi图层作为权重得到。各个子流域的汇流时间为15～40 min，按照5 min的间距reclass进行分类，通过查询属性表中每个时段的栅格个数即可得到各个子流域的等流时面积。图3为第5子流域的等流时面积分布。

图3 第5子流域等流时面积分布Fig．3 The isochrones-divided area of the fifth sub basin

等流时线汇流具体计算方法:假设流域被划分为N块等流时面积，记为fi(i=1，2，…，N)，一场空间分布均匀的净雨有M个时段，记为hj(j=1，2，…，M)，计算时段及等流时线时均距为Δt，此处Δt=5 min，出口断面有流量的时段数 T=M+N－1，每个时段的流域出口断面流量(地面径流)记为QS(k)(k=1，2，…，T)，计算公式为:

5．2．2 壤中流和地下径流汇流

壤中流和地下径流分别采用线性水库调蓄模型和线性水库蓄泄模型模拟其汇流过程，设QSS(k)、RSS(k)分别为第k个时段的壤中流流量和平均壤中流净雨深，设Qg(k)、Rg(k)分别为第k个时段的地下径流流量和平均地下净雨深，KKSS、KKG分别为壤中流的消退系数和地下径流的消退系数，f是子流域面积，则有:

5．2．3 总径流量

第k个时段总径流量Q(k)的计算式为:

5．3 模型参数率定

选取赤桥水文站1979—2008年中较典型的15场洪水进行参数的率定和验证(前10场用于率定，后5场用于验证)，得到的模型参数见表2。图4为其中900731次洪水过程的模拟结果。图5为第5子流域不同标准下的设计洪水。表3是赤桥站模拟与实测峰量对比情况。

表2 模型参数率定值Table 2 Adjusted parameters of the model

续表

图4 900731次洪水过程模拟Fig．4 Simulation of the 900731 flood process

图5 第5子流域不同标准下的设计洪水Fig．5 The design flood in different frequencies of the fifth sub basin

表3 实测与模拟峰量对比表Table 3 The peak flow and quantity of floods obtained by measurement and simulation

表3结果表明径流的合格率较为理想，其中10场洪水峰值相对误差在±10%以内，12场洪水在±20%以内，平均相对误差为10%。另外，其中有7场洪水洪量相对误差在±10%以内。多数场次洪水的确定性系数在0．65以上。总体来说，模拟效果比较好，其中洪峰的模拟比洪量要好，这更有利于小流域设计洪水的计算。

6 设计洪水推求

根据计算得到的不同频率下的设计暴雨过程，利用产流参数在各子流域上进行产流计算，然后利用笔者建议的基于DEM方法进行汇流计算，就可以给出各子流域的设计洪水过程，其洪峰和洪量结果如表4所示。

表4 各子流域设计洪峰洪量统计表Table 4 The design peak flow and quantity of floods in the sub basins

7 对比验证

为了验证表4给出的设计洪水的合理性，也采用水文比拟法计算了各子流域的设计洪峰值。根据赤桥水文站资料分析，流域洪峰流量随面积的递减指数n=0．835。这样，选择赤桥站为参证站，利用式(9)就可以推求不同频率的设计洪峰流量(见表5)。

表5 水文比拟法与DEM方法计算洪峰流量结果对比(单位:m3/s)Table 5 The result of flood peak flow calculated by hydrological analogy approach and DEM approach(unit:m3/s)

对比二者计算结果:当设计标准较低时，基于DEM方法的结果偏大;当设计标准较高时，水文比拟法结果偏大，但是二者基本相近。水文比拟法因未考虑地形、地质、人类活动等的影响，存在一定缺陷。因此，对于资料短缺地区的设计洪水计算，基于DEM的方法可能更合理一些。

8 结语

以缺少径流资料的福州市北部山区为研究流域，根据水系特点将其划分为11个子流域，研究了资料短缺地区的设计洪水推求方法。根据流域资料情况，建议通过暴雨资料推求设计洪水，主要做了以下工作:

1)结合研究流域的实际情况，通过DEM技术与水文模型相结合的方法，探讨资料短缺地区设计洪水的计算方法。

2)通过参证站水文模型参数的移用，解决流域内其他子流域的产汇流参数的估算问题。特别是在汇流计算中，将近年来发展非常迅速的地理信息系统工具运用到研究当中，通过DEM推求等流时线，避免了在地图上人工勾绘等流时线的繁琐工作。随着全国各地数字高程模型的建立，将其运用到水文计算当中有重要的现实意义。

文章采用的这种推求设计洪水的方法为其他资料短缺地区的设计洪水推求提供了一个范例。另外，通过DEM推求等流时线的方法也可与其他水文模型相结合，用于其他流域的径流模拟和设计洪水的推求。

［1］谈 戈，夏 军，李 新．资料短缺地区水文预报研究的方法与出路［J］．冰川冻土，2004，26(2):192－196．

［2］刘苏峡，夏 军，莫兴国．资料短缺流域水文预报(PUB计划)研究进展［J］．水利水电技术，2005，36(2):9 －12．

［3］李红霞，张永强，敖天其，等．无资料地区径流预报方法比较与改进［J］．长江科学院院报，2010，27(2):11－15．

［4］柴晓玲，郭生练，彭定志，等．IHACRES模型在无资料地区径流模拟中的应用研究［J］．水文，2006，26(2):30－33．

［5］陈志明．流域地貌瞬时单位线法剖析［J］．水电能源科学，1993，11(6):105 －111．

［6］李红霞．无径流资料流域的水文预报研究［D］．大连:大连理工大学，2009．

［7］戴志忠．福州地区水文特性分析［J］．水利科技，2004(3):9－11．

［8］詹道江，叶守泽．工程水文学［M］．北京:中国水利水电出版社，2007．

［9］赵人俊．流域水文模拟［M］．北京:水利电力出版社，1984．

［10］任立良，刘新仁．数字高程模型在流域水系拓扑结构计算中的应用［J］．水科学进展，1999，10(2):129 －134．

［11］Maidment D R．Developing a spatially distributed unit hydrograph by using GIS［C］//Kovar K，Nachtnebel H P．Hydro GIS 93，Application of GIS in Hydrology and Water Resources:Proceedings of the Vienna Conference．IAHS，1993(211):181 －192．

［12］熊立华，彭定志．基于数字高程模型的等流时线推求与应用［J］．武汉大学学报(工学版)，2003，36(3):1 －3．

［13］汤国安，杨 昕．Arcgis地理信息系统空间分析实验教程［M］．北京:科学出版社，2006．

［14］熊立华，郭生练．分布式流域水文模型［M］．北京:中国水利水电出版社，2004．