1. 引言
流域水文模型是进行水文科学研究的重要手段之一,是探究水文物理变化规律、解决实际水文生产实践问题的重要工具 [1]。水文模型主要分为集总式水文模型和分布式水文模型。集总式模型将流域视为一个整体,并将降雨、蒸发等水文气象条件均匀分布于整个流域,主要考虑水的垂向运动,计算相对简便。而分布式模型则将流域划分为一个个网格,充分考虑了水文气象信息及土壤下垫面特征在空间上分布的不均匀性,将水的垂向运动与水平向运动均考虑在内,具有更强的物理机制。
随着3S技术(遥感技术(Remote sensing, RS)、地理信息系统(Geography information systems, GIS)、全球定位系统(Global positioning systems, GPS))的成熟和广泛应用,水文气象资料在空间上的分布信息越来越容易获取,分布式水文模型也获得了更广泛的发展与应用 [2]。熊立华等 [1] 于2004年提出了一个基于DEM的分布式降雨径流模型(DEM-based distributed rainfall-runoff model, DDRM)。该模型曾成功应用于英国的SlaptonWood流域 [3]、中国的清江流域 [4]、渠江流域 [5]、赣江流域 [6]、旬河流域 [7]、东江流域 [8] 和北江流域 [9],其改良版还成功应用于典型喀斯特流域——西江流域 [10]。本文将DDRM应用于武隆站上游的乌江流域进行径流模拟,并和三水源新安江模型 [11] (下文简称XAJ模型)的径流模拟效果进行对比分析。
2. 分布式水文模型DDRM
2.1. 模型介绍
DDRM是一个基于栅格单元产流计算和基于栅格流向进行逐栅格分级汇流演算的分布式水文模型 [8]。构建DDRM时,首先利用GIS平台将研究的目标流域划分为相同大小的栅格单元,然后根据已知的水文站点资料将流域划分为数个子流域(图1),从而更加充分的考虑流域内不同区域下垫面特征的差异并使模型运算过程并行化以提高模型运算效率。在模型产流部分,DDRM假定流域的产流机制为蓄满产流,当发生降雨时,雨水直接渗入土壤,地下土壤蓄水量在完成降雨、蒸散发、地下水入流及出流的过程后发生变化,超过栅格土壤蓄水能力的部分将涌出地面形成浅层地表水并在重力作用下产生坡面流汇入栅格河道。在模型汇流部分,DDRM首先利用熊立华等于2007年提出的分级确定法 [12] 来确定各个子流域内的栅格汇流演算顺序,从上游栅格开始,各栅格的产流量都根据栅格汇流演算顺序不断向下游进行演算至各个子流域的出口栅格(节点d和e)处;然后再根据子流域之间的河网拓扑关系对各个子流域出口栅格进行河网汇流演算至全流域出口栅格(节点f)。
2.2. 栅格产流
DDRM的栅格产流计算基于蓄满产流机制进行,即扣除实际蒸散发后的净雨量直接汇入土壤,在完成地下水的水平向交换之后,超出土壤蓄水能力的土壤水涌出地面形成浅层地表水,并在重力作用下产生坡面流汇入河道(图2)。
![](//html.hanspub.org/file/6-2410954x10_hanspub.png)
Figure 1. Discretization of catchment into sub-catchments and grids in DDRM
图1. DDRM栅格及子流域划分示意图 [13]
![](//html.hanspub.org/file/6-2410954x11_hanspub.png)
Figure 2. Hydrological processes involved in a grid in DDRM
图2. DDRM栅格单元产流计算示意图 [13]
为了定义每个栅格处的土壤蓄水能力从而充分考虑流域下垫面条件空间分布的不均匀性,DDRM利用地形指数来计算栅格i处的土壤蓄水能力如下(后文公式中参数的下标i表示栅格i处的参数,下标t表示t时刻时的参数):
(1)
(2)
式中:
为栅格地形指数;
为栅格集水面积;
为地表坡度;
为土壤蓄水能力;S0为全流域最小土壤蓄水能力;
、
分别为全流域最小和最大地形指数;SM为全流域土壤蓄水能力变化幅度。其中S0、n、SM三个参数需通过算法进行优选。
基于蓄满产流机制,记栅格i处t时刻的地下土壤水蓄量的计算值为
,当其大于
时,超出的部分冒出地表产流形成浅层地表水
。
浅层地表水在重力作用下形成坡面流
并汇入栅格河道,根据线性水库法计算如下:
(3)
式中:TP为反映坡面流形成的一个时间参数,需通过算法进行优选。
栅格i处的地下水出流量
与该处的地下土壤水蓄量和地下水力坡度有关,DDRM用地表坡度来近似表示地下水力坡度,
计算如下:
(4)
式中:
为地下水出流阈值,只有当地下土壤水蓄量超出该阈值时才会产生地下水出流,
为反映地下水出流特性的一个参数;TS为反映地下水出流的一个时间参数;
为该处的栅格平均坡度;b为反映坡度对地下水出流影响程度的一个参数。其中
、TS、b需通过算法进行优选。
栅格i处的实际蒸散发量
计算如下:
(5)
式中:
为3.3.2中利用Blaney-Criddle (B-C)公式及反距离加权插值法计算所得的栅格i处t时刻的潜在蒸散发量。
栅格i处的地下水入流量
为其相邻的所有上游栅格的地下水出流量之和。基于水量平衡的原理,根据时刻t时的地下土壤水蓄量,在考虑降雨、实际蒸散发、地下水出流和入流的变化后即可计算出栅格i在下一时刻
时的地下土壤水蓄量。
2.3. 栅格汇流
栅格i内的河道入流量
为其相邻上游各栅格河道出流量
之和。栅格河道的汇流演算采用马斯京根法,在考虑坡面流的情况下,栅格河道出流量计算如下:
(6)
式中:
、
为马斯京根汇流参数,需通过算法进行优选。
2.4. 河网汇流
模型在每个子流域上分别进行栅格产汇流计算得到子流域的出口流量,然后根据河网拓扑关系利用马斯京根法从河段上游节点的入流过程
演算至下游节点的出流过程
从而实现河网汇流演算,其计算公式为:
(7)
式中:
、
为河网马斯京根汇流参数,需通过算法进行优选。
如图1所示,有些栅格可能同时是多条河流的汇流出口,此时可采用线性叠加法进行计算,即认为该栅格处的流量为上游多条河流独立向下游进行演算所得的流量之和,图1中节点f的总出流量
计算如下:
(8)
(9)
(10)
式中:
、
、
分别为节点f、d、e的总出流量;
、
分别为节点d、e演算到节点f的流量;
、
分别为节点d至节点f的河网马斯京根汇流参数,需通过算法优选;
、
分别为节点e至节点f的河网马斯京根汇流参数,需通过算法优选;
为节点f所对应的区间出流量。
若某个节点对应的子流域是头流域(head watershed)时(即图1中的节点d和e),该节点的总出流量即等于其对应的区间出流量,即:
(11)
(12)
式中:
、
分别为节点d、e对应的区间出流量。
3. DDRM在乌江流域的构建
3.1. 研究区域概况
乌江是长江右岸上游最大的一条支流,发源于贵州省境内威宁县香炉山花鱼洞,至涪陵汇入长江。乌江流域上、中游居于云贵高原东部,下游局部位于四川盆地、武陵山区和鄂西南山地,整个流域位于东经104˚18'~109˚22'、北纬26˚7'~30˚22'之间。乌江全长1037 km,其流域面积约为87,920 km2。乌江河段上支流众多,仅在贵州省境内,流域面积大于100 km2的支流就有75条,在整个乌江流域内,有15条一级支流的流域面积超过了1000 km2 [14]。乌江流域大部分地区都属于亚热带季风气候,其年平均气温为10℃~19℃,总体呈现西冷东热的特征,年均降雨量为1100~1900 mm,其汛期(5~10月)的降雨量占全年降雨量的75%以上。
本文选取乌江武隆站以上区域作为研究对象,其面积为83,035 km2。图3给出了研究区域内武隆站及25个气象观测站点的位置分布。研究采用的水文气象资料包括:武隆站1976~1990年的逐日径流序列和25个气象观测站点(气象观测站点根据泰森多边形法选取) 1976~1990年的逐日降雨量和日平均气温序列。并将武隆站的径流序列作为本文的模拟对象。
![](//html.hanspub.org/file/6-2410954x65_hanspub.png)
Figure 3. Topography of Wujiang basin and the locations of a hydrological station and meteorological stations
图3. 乌江流域及水文气象站点图
3.2. 数字流域信息提取
本文采用从http://srtm.csi.cgiar.org/网站上下载的分辨率为3'' (约为88 m)的DEM数据,利用ArcGIS软件将其重采样至1 km的分辨率后进行DEM的预处理,包括填洼、计算栅格流向、计算栅格集水面积、提取河网水系和确定流域边界等过程,然后根据栅格的坡度和集水面积来计算栅格地形指数。由于仅有武隆站一个水文站的日径流序列,因此本文不对乌江流域进行子流域的划分。
3.3. 水文气象输入数据
3.3.1. 降雨
分布式水文模型相对于集总式水文模型而言最主要的优势就是能够充分考虑水文气象资料在空间上的不均匀分布。根据研究表明,降雨及时空分配是水文模型最为关键的输入数据 [15],10%的降雨误差可以导致35%的径流误差 [16],因此在空间上构建合理的降雨数据对提升模型的模拟效果至关重要。由于站点降雨数据是离散的,为了能让流域内的降雨数据在空间上连续分布,本文采用反距离加权插值法 [17] 对站点的降雨数据进行插值计算。反距离加权插值法根据相近相似的原则,每个站点的数据对流域内的每个点都有一定的权重,这个权重随着距离的增加而减小,即站点数据对距离越近的点贡献越大,对距离越远的点贡献越小。本文以处理后的DEM栅格为基准进行降雨数据的空间插值。
3.3.2. 潜在蒸散发
由于缺少实测蒸散发数据,本文采用以温度为基础的B-C公式来计算潜在蒸散发量。t时刻的潜在蒸散发量
计算如下:
(13)
式中:k为一经验系数,取0.85;
代表t时刻的白昼时间占全年白昼时间的比例;
为t时刻的日平均温度。
在计算得到所有气象站点的逐日潜在蒸散发数据之后,用反距离加权插值法对站点潜在蒸散发数据进行插值计算从而获得在空间上连续分布的潜在蒸散发数据。
3.4. 模型参数率定
选用由美国亚利桑那大学的Duan等 [18] 于1992年提出的SCE-UA算法,SCE-UA算法基于复合型搜索和自然界生物进化的原理提出,是一种全局优化算法,能快速、高效地搜索到参数的全局最优解 [19]。
由于缺少水库调度的实际数据,为了尽量避免修建水库等其他人为因素的影响,本文在已有资料的基础上尽量选择早期数据进行研究,因此本文选取1976年1月1日至1979年12月31日作为模型率定期,1980年1月1日至1981年12月31日作为模型验证期,将模型预热期设置为100d,并采用模拟日径流的纳什效率系数(Nash-Sutcliffe Efficiency, NSE)作为模型率定的目标函数,NSE的值越接近1则代表模型的模拟效果越好,其计算式如下:
(14)
除了以NSE作为模型模拟效果的评定指标外,还选取径流深相对误差(Relative Error, RE)作为评定指标,其计算式如下:
(15)
式中:T表示总时段数目;
为时刻t的模拟径流量;
为时刻t的实测径流量;
为全时段实测径流量的平均值。
DDRM各参数的物理意义及率定结果见表1,由于本文未对乌江流域划分子流域,无需进行河网汇流过程,故不考虑
和
两个参数。
此外,本文还在乌江流域构建XAJ模型进行日径流模拟,其原理参考文献 [11],并与DDRM的模拟结果进行对比,XAJ模型的参数率定同样采用SCE-UA算法,其参数的物理意义及率定结果见表2。
![](Images/Table_Tmp.jpg)
Table 1. Descriptions of the DDRM parameters
表1. DDRM参数表
![](Images/Table_Tmp.jpg)
Table 2. Descriptions of the XAJ model parameters
表2. XAJ模型参数表
4. 结果分析
DDRM及XAJ模型在率定期和验证期的径流模拟效果见表3。
![](Images/Table_Tmp.jpg)
Table 3. Evaluation of the model-simulated runoff from DDRM and Xinanjiang model
表3. DDRM和XAJ模型径流效果评价表
分析表3中数据:从NSE的角度来说,在率定期4年数据中,DDRM和XAJ模型的NSE均大于0.7,且DDRM的NSE有2年大于0.8,其4年的模拟效果整体较XAJ模型更优,在验证期2年数据中,DDRM和XAJ模型在1980年的模拟效果较好,但1981年的模拟情况较差,其原因将在后文进行猜测分析,总体来说DDRM和XAJ模型对于乌江流域这一湿润流域模拟效果较好,且DDRM模拟效果略优于XAJ模型;从RE的角度来说,在全部6年的数据中,DDRM和XAJ模型互有优劣,总体来说模拟效果基本相当。
在6年的模拟中,1981年的模拟情况非常差,该年出现了实测径流深异常小的情况,在其年降雨量及平均气温与其他年份大致相当的情况下,该年出现极端气候现象的可能性较小,因此推测该异常现象由人类活动导致。根据文献 [14] 中的资料显示:乌江干流梯级开发的第六个电站,也是第一座大型水电站,乌江渡水电站,在1981年完成了第2台机组的建成投产。推测可能是由于水电站的建设过程中修建引水隧洞、大江截流等施工过程以及投产后的调度过程等大规模人类活动从一定程度上破坏了短期内水文资料的一致性从而导致1981年模拟效果不佳。
![](//html.hanspub.org/file/6-2410954x82_hanspub.png)
Figure 4. Time series of observed runoff and simulated runoff over the flood season in 1977
图4. 1977年洪水季节实测与模拟日径流量过程对比图
![](//html.hanspub.org/file/6-2410954x83_hanspub.png)
Figure 5. Time series of observed runoff and simulated runoff over the flood season in 1978
图5. 1978年洪水季节实测与模拟日径流量过程对比图
在DDRM和XAJ模型模拟效果都较好的年份中选取径流深最大和最小的年份(1977年和1978年)中洪水季节的模拟日径流序列进行绘图展示(图4、图5)。从图4和图5中也可以发现,DDRM与XAJ模型在峰现时间的模拟上与实测径流基本吻合,但在洪峰量级的模拟上,DDRM略优于XAJ模型。
除了能够模拟流域出口的日径流序列,DDRM还具有模拟流域内土壤含水率
和径流量Q在时空间上分布的能力。本文选取19780525场次洪水进行研究,DDRM对本次洪水模拟的NSE为0.72,洪峰流量相对误差为−9.06%,模拟峰现时间与实测峰现时间仅相差1 d,模拟效果较好。本次洪水的流域出口在整个洪水期间的径流模拟情况如图6所示。
![](//html.hanspub.org/file/6-2410954x85_hanspub.png)
Figure 6. Time series of observed runoff and simulated runoff from DDRM during “19780525” flood
图6. 19780525号场次洪水实测流量与DDRM模拟流量过程对比图
图7给出了该次洪水在乌江流域中不同时刻的日土壤含水率的空间分布情况。在涨水时刻(19780528),流域大部分区域土壤含水率在0.8左右。在洪峰时刻(19860531),由于连续的强降雨,此时流域中几乎大部分地区都趋于饱和。在退水时刻(19780604),整个乌江流域的土壤含水率明显降低。
![](//html.hanspub.org/file/6-2410954x86_hanspub.png)
Figure 7. The spatial distribution of simulated soil moisture contents from DDRM on the four days during “19780525” flood
图7. 乌江流域在19780525场次洪水中不同时刻土壤含水率的空间分布图
图8则给出了该次洪水在乌江流域中不同时刻的日径流量空间分布情况,根据图中信息可以发现,径流量在每一段河网上都呈现从上游往下游逐渐增大的趋势,且属于河网水系的栅格处的径流量明显大于非河网水系处栅格的径流量。从整个时段来看,洪峰时刻(19780531)的径流量略大于涨水时刻(19780528),明显大于退水时刻(19780604),在退水阶段,随着时间的推移,全流域的径流量逐渐减小,径流量较高的区域基本属于乌江河网水系部分。
![](//html.hanspub.org/file/6-2410954x87_hanspub.png)
Figure 8. The spatial distribution of simulated runoff values from DDRM on the four days during “19780525” flood
图8. 乌江流域在19780525场次洪水中不同时刻日径流量的空间分布图
5. 结语
本文以武隆站上游的乌江流域为研究对象构建DDRM进行1976~1981年的日径流模拟,以纳什效率系数NSE和径流总量相对误差RE这两个指标进行精度评定,并选用三水源新安江模型对比分析。得出结论如下:在保证水文资料一致性不被破坏的前提下,DDRM径流模拟的NSE较新安江模型略优,径流总量相对误差RE的模拟效果与新安江模型基本相当。总体而言,DDRM在乌江流域的模拟效果略优于新安江模型。
DDRM除了能模拟流域出口处全时段的径流数据,也能够模拟流域内任一时刻土壤含水率及径流量的空间分布,本文选取19780525场次洪水对DDRM的这一能力进行了展示。
由于缺少水库调度的数据,本文未能进一步模拟水库投产运行后的日径流序列,这一不足之处可以作为以后的研究内容。
总的来说,DDRM结构简单、参数较少且物理过程明确,在已经成功应用于多个湿润地区大中型流域的基础上,本文将其成功应用于乌江流域,进一步表明了DDRM具有较大的推广价值,值得在其他流域进一步探究其应用前景。
基金项目
国家重点研发计划项目(2017YFC0405901)和国家自然科学基金项目(51525902)。