近年来, 受剧烈的区域土地利用变化和全球气候变异的综合影响, 地表水环境系统受到人为活动的影响越来越深刻^{[1, 2]}.由于自然和社会经济等诸多因素的作用, 水资源分布及变化特征的准确量化成为水文水资源领域研究的热点和难点之一.因此, 为实现区域水资源可持续发展, 合理评价水资源的时空特征对于提高区域水资源利用效率具有重要而深刻的现实意义和价值.

数值模拟是进行水文分析、洪水预报和水资源评价与管理的最重要和有效的手段之一^[3].水文模型能够定量描述地表复杂水文循环过程, 是模拟水文现象特征与变化过程的最有效的工具.常见的水文模型有SWAT(soil water assessment tool)模型^[4]、HSPF(hydrological simulation program-Fortran)模型^[5]、TOPMODEL模型^[6]、SWMM(storm water management model)模型^[7]、SLURP(semi-land use runoff program)模型^[8]、MIKE-SHE模型^[9]和新安江模型^[10]等.其中, HSPF模型以Stanford水文模型为基础, 能够综合模拟地表径流、土壤流失、污染物传输、河道水力等过程, 并广泛应用于下垫面差异及气候变化对流域地表水环境影响研究^[11], 诸如土地利用变化和气候变异对区域水文过程的影响^[12], 工业点源污染评价和农业非点源污染模拟等^{[13, 14]}.但是, 受限于模型结构简化、输入观测数据和模型参数等因素的不确定性影响, 流域模型模拟结果存在较大的不确定性^{[15, 16]}.国内外学者对HSPF模型的不确定性问题进行了深入研究.例如, Xie等^[17]采用多准则和GLUE算法研究了SWAT和HSPF模型的不确定性, 结果表明参数校准过程对HSPF模型模拟精度的影响更明显, 而SWAT模型则更适合于数据缺乏区域; Fonseca等^[18]采用GLUE算法研究了HSPF模型参数敏感性和不确定性, 结果表明径流模拟结果对土壤、土地利用和气象因素变化较为敏感, 敏感参数INFILT和LZSN与土壤和土地利用密切相关; 程晓光等^[19]以永定河水系的妫水河流域为例, 采用GLUE算法探究了半干旱地区HSPF模型的敏感性和不确定性, 结果表明HSPF模型参数间存在复杂的“异参同效”和相关性现象, 另外降雨量是影响月尺度模型不确定性的最重要的因素.但是, 上述研究主要从单一时间尺度(例如月尺度)的视角研究模型的不确定性, 关于多时间尺度变化对HSPF模型水文模拟的不确定性的影响研究却鲜见报道, 使得选择何种时间尺度能够有效地降低水环境管理决策过程中存在的不确定性, 是一个学术界和管理部门不可回避的科学问题.基于此, 本研究通过构建潮河流域HSPF模型, 采用GLUE算法研究了模型主要参数的敏感性、不确定性以及参数间的相关性, 并重点分析了年尺度、季节尺度和月尺度下地表径流预测的不确定性.

密云水库是北京市最重要的地表饮用水源地.自1999年以来, 受气候变化和人为活动的综合作用, 密云水库流域地表径流量减少趋势极为显著, 因而水资源缺乏已经成为影响本地区经济发展的主要限制条件之一^[20].本文以国内外广泛使用的流域水文模型HSPF模型为研究对象, 探究模型参数的敏感性和复杂的“异参同效”现象, 以及不同时间尺度下地表径流预测的不确定性, 以期为流域水资源综合调控和生态建设提供技术支持.

潮河流域位于北京市东北部, 地理位置介于E 116°27′~117°34′, N 40°20′~41°27′, 河流长度239.50 km, 流域面积4 888 km², 河网密度0.05 km·km^-2 (图 1), 辖3县24乡镇, 其中北京市密云区北部共计3个乡镇.气候为典型的中温带向暖温带过渡, 半干旱向半湿润过渡的大陆性季风气候, 年平均降水量494 mm, 主要集中在6~10月.棕壤和褐土广布, 约占全流域面积的80%;植被类型多样, 以针阔混交林为主; 流域总人口37万人, 人口密度约为75.70人·hm^-2, 以农业人口为主, 约占77%;产业结构以畜牧业和种植业为主, 主要畜牧类型有大牲畜、猪、羊和家禽; 主要农作物类型为玉米和小麦, 农田管理方式和化肥管理粗放, 平均化肥施肥量225~300 kg·hm^-2.近50年来, 受流域生态恢复、水利设施建设和跨流域调水等的影响, 潮河流域地表径流量总体呈现减少趋势, 且自1999年以后减少趋势愈来愈显著^[21].

图 1

Fig. 1

图 1 研究区位置示意 Fig. 1 Location of study area

本研究所需的数据包括数字高程模型(digital elevation model, DEM)、土地利用数据、气象数据、水文观测数据等, 详细数据信息参见表 1.其中空间数据均采用Albers等积圆锥投影对数据进行预处理, 以统一各类空间数据的投影和地理坐标, 所有操作均由ArcGIS 10.1完成.

表 1 (Table 1)

表 1 研究所需数据及来源 Table 1 Content of data and sources 年份 类型 内容 尺度 格式 来源 1980 数字高程模型 高程、坡度和坡向 30 m 栅格 中科院资源环境科学数据中心 2012 土地利用 土地利用及空间分布 1:10万 矢量 中科院资源环境科学数据中心和21世纪空间技术有限公司 2005 土壤类型 土种类型及理化性质 1:100万 矢量 中科院资源环境科学数据中心 1998~2010 降雨量 降雨量 逐日 txt 中国气象科学数据共享服务网 1998~2010 气象数据 气温、风速、太阳辐射和露点温度 逐日 txt 中国气象科学数据共享服务网 1998~2010 水文数据 月均径流总量 逐月 xls 密云水库管理处下会水文站

表 1 研究所需数据及来源 Table 1 Content of data and sources

HSPF模型是由美国环保署于20世纪80年代开发的一款优秀的半分布式流域水文水质模型, 能够适用于不同气候情景和土地利用影响下的透水地面、不透水地面、混合型河湖水库等水体水文水质过程模拟^[22].目前, HSPF模型已经集成于BASINS系统中, 可以对地形数据、土地利用和土壤类型数据进行预处理和实现子流域单元的刻画; 利用WDMUtil工具将建模所需的气象和水文数据整合生成WDM文件, 其中不同子流域气象数据的分配方法采用泰森多边形法, 再导入重分类的土地利用数据, 生成潮河流域HSPF模型所需数据库, 进而构建WinHSPF工程.本研究采用BASINS 4.0系统自带的流域划分模块, 将潮河流域分为21个子流域.基于此, 对潮河流域多年径流过程的不确定性进行量化研究.

GLUE算法(generalized likelihood uncertainty estimation), 即普适似然不确定分析方法, 该算法由Beven等于1992年首次提出^[23], 并广泛应用于模型参数敏感性和不确定性分析研究^{[24, 25]}.其基本思想是决定模型模拟效果优劣的并不是模型单个参数, 而是模型的参数组合. GLUE算法的主要步骤：①选择特定的似然函数及似然值阈值, 以此判断参数组的优劣; ②随机生成若干参数组合代入模型中, 对于低于该似然值阈值的似然值赋值为0, 高于似然值阈值的似然值作归一化处理; ③将符合条件的参数模拟结果按照从大到小顺序排列, 计算得到不同置信水平下(诸如90%或95%水平)模型预报的不确定性范围.本研究以纳什系数为似然函数, 选择0.2作为判别模型参数有效性的阈值标准, 采用拉丁超立方技术取样方法随机选择3万组参数, 代入HSPF模型进行参数敏感性和不确定性分析.

本研究采用纳什系数(Nash-Sutcliffe efficiency, NSE)和相对误差(relative error, RE)作为模型校准和不确定分析的评价指标, 以表征模型模拟值与实际观测值间的匹配程度.

(1)

(2)

式中, P_i为模拟值, O_i为观测值, O_avg为观测值的均值, n为模拟时段长度.其中NSE∈(-∞, 1], 其值越接近于1表示模拟结果越好; RE值越趋近于0, 表明模拟值与监测值越接近, 误差越小.一般地, 纳什系数NSE大于0.5, 相对误差RE小于15%, 则认为模型模拟结果较好.

由于模型结构复杂参数众多, 为提高HSPF模型率定的时效性, 本文采用参数敏感性和估计软件PEST筛选出HSPF模型中影响径流的主要敏感性参数(见表 2).以潮河流域下会水文站径流量为基础, 对HSPF模型进行参数校准与验证.其中, 1998~1999年为模型的预热期, 2000~2005年为模型的校准期, 2006~2010年为验证期.参考相关研究结果, 通过手动调整影响径流量的主要参数, 使径流的模拟值与实测值相接近.如图 2可知, 校准期和验证期的纳什系数(NSE)分别为0.84和0.53, 相对误差(RE)分别为4.40%和2.80%, 参考相关模型评价标准, 表明率定后的HSPF模型能够较好地反映潮河流域径流量变化规律^[26].

表 2 (Table 2)

表 2 HSPF模型敏感参数及取值范围 Table 2 Sensitive parameters and value range for the HSPF model 参数名称 中文释义 单位 率定值 取值范围 物理意义 LZSN 下层土壤含水量 In 5.25 2~15 LZSN参数值越大, 土壤蓄水能力越大, 潜在蒸发量越大, 引起产流下降 UZSN 上层土壤含水量 In 1.14 0.05~2 UZSN参数值增加, 上层蓄水容量和蒸发能力提高, 坡面流减少 INFILT 入渗能力 In·h -1 0.11 0.001~0.5 INFILT参数控制着降雨在地表径流、地下径流及土壤含水量之间的分配比例, 参数值增加, 地表径流减小, 基流增加, 反之, 壤中流和坡面流增加 IRC 壤中流消退系数 d-1 0.52 0.30~0.85 IRC参数是当前日壤中流出流与24 h前壤中流出流的比率, 影响峰值到基流的落水阶段, 影响单次暴雨时间 BASETP 基流蒸散发参数 — 0.12 0.001~0.2 BASETP参数增加, 径流减少 AGWETP 地下水蒸散发参数 — 0.13 0.001~0.2 AGWETP参数增加, 径流减少 AGWRC 地下水退水率 d-1 0.89 0.85~0.999 AGWRC参数控制着地下水退水过程, 是地下水出流速率与1天前速率之比, 该值越高, 产流越小 CEPSC 截留量 In 0.25 0.01~0.4 — INTFW 壤中流入流参数 — 5.00 1~10 INTFW参数是决定截留水进入地下成为壤中流的系数, 该参数增加, 壤中流增加, 坡面流减少, 影响暴雨峰值, 不影响基流 DEEPF R 地下水出流进入深层地下水的比例 — 0.085 0.001~0.5 DEEPFR参数值越大, 产流越小

表 2 HSPF模型敏感参数及取值范围 Table 2 Sensitive parameters and value range for the HSPF model

图 2

Fig. 2

图 2 HSPF模型径流率定与验证过程线 Fig. 2 Runoff hydrograph for calibration and validation of the HSPF model

以纳什系数为似然目标函数, 采用拉丁超立方抽样方法随机抽取3万组样本, 对HSPF模型进行参数敏感性和不确定性分析.本研究选取似然函数临界值为0.2, 低于该值的参数似然值赋值为0, 高于该值的似然值重新归一化并以似然值大小进行排序.后续的参数间相关性和模型预测的不确定性也是基于该临界值展开分析.为直观展示参数取值与似然函数关系, 二者散点图中似然函数临界值设定为0(图 3).根据模型参数取值区间内似然函数值散点分布差异性, 本文将模型参数分为全局敏感性参数、局部敏感性参数和不敏感参数.

图 3

Fig. 3

图 3 模型参数与似然值散点图 Fig. 3 Scatter plot of model parameters and likelihood value

全局敏感性参数：由图 3可知, 参数LZSN、INFILT、IRC和AGWRC的似然值在整个区间内变化幅度较大, 说明径流模拟结果对上述4个参数变化的响应较为敏感, 其中LZSN在2~15间随着参数值增大似然值呈减小趋势; INFILT在0.001~0.07间似然值随着参数的增大而增大, 从0.07~0.5间似然值随着参数的增大而减小; AGWRC在0.995~0.999间似然值随着参数的增大而迅速增大, 在0.98~0.995间散点较为稀疏, 究其原因是参数值在这个区间内的似然值以负值为主.

局部敏感性参数：参数UZSN的敏感性具有一定的区间特征, 即在0.05~1.2间似然值随着参数的增大而增大, 在1.2~2间似然值相对稳定.

不敏感性参数：参数DEEPFR、BASETP、AGWEPT、INTFW和CEPSC无明显的变化规律, 表明上述参数对潮河流域径流模拟影响并不显著, 对HSPF模型径流模拟预测的不确定性影响较小.因而, 可以将上述参数视为不敏感参数.综上可知, LZSN、INFILT、IRC和AGWRC这4个参数是对潮河流域水文过程和水平衡影响最大的参数, 这些参数与其他学者研究结果相类似^[27].但是, DEEPFR、BASETP、AGWEPT、INTFW和CEPSC等参数对该流域水文过程影响并不显著, 究其原因可能是流域植被覆盖度是影响参数敏感性的重要因子之一, 覆盖度通过控制植被截留量和基流蒸发量, 进而影响流域地表径流过程.另外, 流域地表特征及地形起伏也会导致影响参数敏感性差异(如UZSN)^[28].

目前, 模型不确定性分析研究大多假定参数间是相互独立的, 研究表明参数间复杂的相关性可能导致参数冗余, 是参数不确定性的来源之一^[19].基于此, 本研究对所识别出的5个敏感参数进行相关性分析(表 3).结果表明, LZSN与INFILT呈极显著负相关; 与UZSN呈极显著正相关; INFILT与UZSN呈极显著负相关, 与AGWRC呈显著正相关关系; UZSN和AGWRC呈极显著负相关关系.由此可见, 忽略不同参数间相关性差异也可能会对模型不确定性产生较大影响.

表 3 (Table 3)

表 3 HSPF模型敏感性参数相关系数矩阵1) Table 3 Correlation matrix of sensitive parameters in the HSPF model LZSN INFILT UZSN IRC AGWRC LZSN 1.000 INFILT -0.458** 1.000 UZSN 0.290** -0.696** 1.000 IRC 0.041 0.055 -0.007 1.000 AGWRC -0.071 0.155** -0.109** -0.061 1.000 1)**代表在0.01水平显著相关, 其他代表不相关

表 3 HSPF模型敏感性参数相关系数矩阵1) Table 3 Correlation matrix of sensitive parameters in the HSPF model

由于模型结构的不完备性和监测资料误差, 因此通过参数调整得到的最优参数并不能代表模型的“真实”参数值, 它只是在特定的目标函数下的似然估计的优化参数值.有研究表明, 无论是人工调试还是计算机自动优选都不可避免“异参同效”现象, 即不同组参数组合得到相同模拟效果^[29].由于模型“异参同效”参数组众多, 在综合考虑参数空间长度和跨度的基础上, 本文中等效参数组的选取尽量兼顾每个参数的最大取值空间, 同时限于篇幅原因, 在此给出了潮河流域5组具有代表性的等效参数组(表 4).这5组参数组合均可使日径流纳什系数达到0.55, 说明率定产生的“最优”参数往往不是唯一的, 即具有很强的不确定性.其中, 空间长度是参数最初取值范围长度, 跨度是参数似然值为0.55时参数取值的最大值与最小值之差, 百分比是跨度与空间长度的比值.以参数LZSN为例, 其参数空间为^{[2, 15]}, 即空间长为13;当似然值为0.55时, 最大值为3.73, 最小值为2.11, 其跨度为最大值与最小值之差(1.62);百分比就是跨度与空间长度的百分比(12.46%).百分比小说明散点图中高似然值区间参数取值较为集中, 参数敏感性程度较高. LZSN、INFILT、IRC和AGWRC的百分比均低于15%, 而参数UZSN的百分比相对偏大, 但是其值小于25%.其他不敏感参数的百分比均大于25%.结合4.2结果来看, 参数LZSN、INFILT、IRC和AGWRC属于全局敏感性参数, 参数UZSN是局部敏感性参数, 其余为不敏感参数, 二者结论基本一致.

表 4 (Table 4)

表 4 异参同效参数组及统计信息 Table 4 Statistical information for the equifinality parameter sets 组别 LZSN INFILT DEEPFR BASETP AGWETP UZSN INTFW CEPSC IRC AGWRC NSE 1 2.78 0.07 0.49 0.09 0.17 1.71 6.21 0.20 0.57 0.99 0.55 2 2.23 0.06 0.30 0.07 0.14 1.27 2.75 0.20 0.56 0.99 0.55 3 2.11 0.10 0.41 0.12 0.12 1.49 8.46 0.11 0.60 0.99 0.55 4 3.36 0.03 0.49 0.04 0.13 1.75 7.27 0.27 0.59 0.99 0.55 5 3.73 0.06 0.32 0.17 0.12 1.52 3.58 0.04 0.55 0.99 0.55 空间长 13.0 0.49 0.49 0.19 0.15 1.95 9.00 0.39 0.55 0.15 — 跨度 1.62 0.07 0.19 0.12 0.04 0.48 5.70 0.22 0.04 0.0013 — 百分比/% 12.46 14.12 39.09 61.38 26.50 24.95 63.37 58.13 8.77 0.87 —

表 4 异参同效参数组及统计信息 Table 4 Statistical information for the equifinality parameter sets

采用GLUE方法研究模型模拟的不确定性, 需要确定似然函数(纳什系数)的临界值, 对高于该值的似然值作归一化处理, 将模拟值从大到下排序, 进而估算模型模拟的不确定性时间序列.本文采用90%置信水平的模拟值代表其预测的不确定性的界限, 得到模型在90%置信水平下的年径流、季径流和月径流的不确定性预测结果(图 4).由图 4可知, 不同时间尺度下径流模拟的不确定性具有一定的相似性, 即不确定性与降雨量和径流量大小密切相关.从年尺度来看, 2001、2005、2008及2010年降雨充沛径流量大, 模拟的不确定范围也较大; 从季节尺度来看, 夏季降雨集中径流量偏高, 导致模拟的不确定性范围大于其它季节; 从月尺度来看, 每年的6~8月是潮河流域的雨季, 降雨量径流量均较大, 使得不确定性预测的范围也较大; 由此可见, 模型的不确定性预测范围随降雨量变化而变化, 与降雨量呈明显的正相关关系, 即丰水期径流量偏高, 相应的模型预报的不确定性偏大; 反之, 枯水期径流量偏低, 模型预报的不确定性偏小.其它水文模型研究结果也证实了降雨量变异是影响水文模拟不确定性的重要因素之一^{[30, 31]}, 与本研究结果基本一致.

图 4

Fig. 4

图 4 不同时间尺度径流预报不确定性范围 Fig. 4 Uncertainty range of runoff simulation at various temporal scales

从模型不确定性可接受范围来看, 年尺度下不确定性范围包含了81.80%的观测值, 季节尺度下不确定性范围包含了78.70%的观测值, 月尺度不确定范围包含了80.56%的观测值, 属于可接受范围.另外, 尚有部分径流实测值未落于不确定性范围内, 说明模型模拟并不能完全还原流域水文过程, 究其原因可能是一方面HSPF模型结构与假设条件与复杂的流域水文地质特征存在一定差异, 另一方面似然值的临界值选取具有较大的主观性, 不同的临界值会产生不同的不确定性范围^[19].

(1) 对潮河流域而言, 所构建的HSPF模型径流模拟取得了较好的效果, 率定期和验证期的纳什系数分别为0.84和0.55, 相对误差分别为4.40%和2.80%, 表明HSPF模型能够真实再现潮河流域历史水文过程.

(2) 影响HSPF模型的主要参数可以分为3类, 即全局敏感性参数(LZSN、INFILT、IRC和AGWRC)、局部敏感性参数(UZSN)和不敏感参数(DEEPFR、BASETP、AGWEPT、INTFW和CEPSC).

(3) 不同敏感性参数存在复杂的相关关系, 参数组(LZSN与INFILT)、(INFILT与UZSN)和(UZSN与AGWRC)间均呈极显著负相关关系; (LZSN与UZSN)和(UZSN与AGWRC)呈极显著正相关.

(4) HSPF模型参数间存在大量复杂的“异参同效”现象, 进一步证实了影响模拟结果优劣的是各参数组合而非某一单个参数值.

(5) 模型预测的不确定性范围与降雨量紧密相关, 降雨量越大不确定性就越大, 反之亦然, 由此说明降雨季节差异是引起模型预报不确定性偏大的主要因素.

(6) 不同时间尺度下HSPF模型的径流模拟效果总体较好, 但也存在一定差异(年尺度、季节尺度和月尺度下不确定性范围分别包含了81.80%、78.70%和80.56%的观测值), 即年尺度效果略优于月尺度和季节尺度.