《1 前言》

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 流域概况》

2 流域概况

福州市位于福建省的东部、闽江中下游，总面积为 12 184 km2 ，海拔多为600 ~1 000 m，属于典型河口盆地地貌。 全区地貌类型以山地、丘陵为主，占全区总面积的 72.7 %。福州地区年平均气温为19.1 ~20.1 ℃，年平均降水量为1 200 ~1 700 mm[7]。

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

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

《图1》

图1　研究流域 DEM

Fig.1　DEM of the research basin

《3 方案设计》

3 方案设计

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

《3.1 设计暴雨计算》

3.1 设计暴雨计算

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

《3.2 产流计算方案》

3.2 产流计算方案

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

《3.3 汇流计算方案》

3.3 汇流计算方案

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

《4 设计暴雨推求》

4 设计暴雨推求

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

《图2》

图2　北部山区最大24 h 降雨量频率分布曲线

Fig.2　The frequency distribution curve of the biggest 24 h rainfall in northern mountain area

《表1》

表1　不同频率下的设计暴雨值

Table 1　Design rainstorm in different frequencies

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

《5 模型计算》

5 模型计算

《5.1 产流计算》

5.1 产流计算

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

5.1.1 蓄满产流

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

5.1.2 划分三水源

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

5.1.3 蒸散发计算

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

《5.2 汇流计算》

5.2 汇流计算

5.2.1 地面径流汇流

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

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

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

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

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

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

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

《图3》

图3　第5 子流域等流时面积分布

Fig.3　The isochrones-divided area of the fifth sub basin

等流时线汇流具体计算方法：假设流域被划分为 N 块等流时面积，记为 （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 模型参数率定》

5.3 模型参数率定

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

《表2》

表2　模型参数率定值

Table 2　Adjusted parameters of the model

《图4》

图4　900731 次洪水过程模拟

Fig.4　Simulation of the 900731 flood process

《图5》

图5　第 5 子流域不同标准下的设计洪水

Fig.5　The design flood in different frequencies of the fifth sub basin

《表3》

表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 设计洪水推求》

6 设计洪水推求

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

《表4》

表4　各子流域设计洪峰洪量统计表

Table 4　The design peak flow and quantity of floods in the sub basins

《7 对比验证》

7 对比验证

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

式（9）中，分别为设计断面和参证站洪峰流量均值 , m3/s；F设 、F参 分别为设计断面和参证站汇水面积，km2 。

《表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 结语》

8 结语

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

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

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

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