基于PEST的HSPF模型径流模拟优化
杨博1
,
陈莹1,2,3
,
陈兴伟1,2,3,
刘梅冰1,2,3,
高路1,2,3
1. 福建师范大学地理科学学院, 350007, 福州;
2. 湿润亚热带山地生态国家重点实验室培育基地, 350007, 福州;
3. 福建省陆地灾害监测评估工程技术研究中心, 350007, 福州
收稿日期：2017-06-15; 修回日期：2018-03-05
项目名称：国家自然科学基金"土地利用变化对流域洪水风险的影响及其不确定性分析"（41301031）；福建省科技厅省属公益类科研专项"晋江流域非一致性洪水风险演变与归因"（2018R1101006-3）
摘要：参数率定是水文模型构建的关键步骤之一，为了提高模型的率定效率和模拟效果，针对HSPF模型，以晋江山美水库流域为例，将PEST应用于日径流模型的率定，分别采用单目标函数和多目标函数率定参数，并对日流量偏差这一目标函数进行改进，划分为丰水期和枯水期的日流量偏差，对比分析了单目标与多目标率定结果的优劣及多目标函数改进前与改进后率定结果的差异。结果表明：1）相比单目标率定，多目标率定的NSE、R2、RMSE和PBIAS等结果均得到不同程度的提高。2）2种多目标率定的结果，在对日流量偏差改进后，枯水期日径流的模拟结果得到显著提高，率定期和验证期的NSE分别由-0.04和0.43提升至0.25和0.54，R2分别由0.37和0.59提升至0.47和0.64，PBIAS分别由-33.57%和-22.72%提升至-24.52%和-16.79%，RMSE分别由9.54和8.89 m3/s降低至8.06和7.96 m3/s。因此，相比单目标率定，多目标率定能更好的把握径流变化趋势、径流总量和流量保证率的特征；并且将日流量偏差划分丰、枯水期之后，研究区的枯水期径流模拟效果得到显著提高。
关键词：HSPF PEST 径流 率定
HSPF runoff simulation and optimization based on PEST automatic calibration
1. College of Geographical Science, Fujian Normal University, 350007, Fuzhou, China;
2. Cultivation Base of State Key Laboratory of Humid Subtropical Mountain Ecology, 350007, Fuzhou, China;
3. Fujian Provincial Engineering Research Center for Monitoring and Assessing Terrestrial Disasters, 350007, Fuzhou, China
地表径流是造成土壤侵蚀的主要原因之一[1]。近年来，由于气候变化和人类活动的影响，流域水循环和水资源问题趋于复杂化。水文模型是对自然界中复杂水文现象的近似模拟[2]，可以实现流域径流定量化研究，是水文科学研究的一种手段和方法。HSPF模型是半分布式水文模型的优秀代表，在水文过程模拟、气候与土地利用变化的水文响应研究等方面得到了广泛应用[3]。参数率定是水文模拟中不可避免的重要环节[4]。目前HSPF模型主要通过GA (Genetic Algorithm)[5]、SCE(Shuffled Complex Evolution Algorithm)[6]、GLUE(Generalized Likelihood Uncertainty Estimation)[7]以及人工率定等方法，且主要基于单目标函数实现参数的率定。PEST(Parameter Estimation)自动率定程序是由Doherty博士基于GML(Gauss-Marquardt-Levenberg)算法开发的，其具有逆海森方法和最速下降法的优点，可以通过较少的模型运行次数，得到最优的参数结果[8]。此外，与其他自动率定方法相比，PEST对建模者编程能力的要求较低，不需要通过编写代码即可实现水文模型和优化算法间的耦合，且基于WINDOWS平台运行的特性让操作简单易行，已在地下水模型中广泛应用[9-10]。目前，基于PEST实现HSPF模型多目标函数率定的相关研究较少[11]，且在率定中多侧重于洪水流量，对枯水的关注不足；因此，本研究以山美水库流域为例，将PEST应用于HSPF模型中，开展单目标与多目标率定的对比研究，并根据丰水期和枯水期对多目标率定的目标函数进行优化，分析其对径流模拟的影响，为流域水文模型的参数率定提供参考。
1 研究区概况
山美水库流域位于福建省晋江流域的东北部(E 118°3′~180°30′，N 25°9′~25°32′)，海拔在60~1360m之间，流域面积1023km2，占晋江流域总面积的18.2%。流域属于亚热带季风气候，年均温约为20.9，年降水量约为1600mm，干湿季分明，降水空间分布上自西北向东南逐渐减少。在土地利用方面，研究区主要以林地为主(64.4%)，其次为园地(14.3%)、耕地(12.1%)、建设用地(5.1%)、水域(2.2%)、未利用地(1.1%)和草地(0.8%)。
2 研究方法
2.1 数据来源
HSPF模型的构建需要涉及气象、水文、DEM和土地利用等数据。笔者使用的气象数据包括2001—2010年永春、德化2个气象站日最高气温、最低气温、相对湿度和风速数据，以及15个雨量站日降雨量数据；水文数据包括2001—2010年山美水库日入库、出库流量和龙门滩月调水数据；DEM数据分辨率为30m×30m，从“中国科学院国际科学数据服务平台”获得；土地利用数据来自山美水库2006年TM遥感影像的解译，共划分水域、草地、未利用地、建设用地、耕地、园地和林地7种类型。
2.2 HSPF模型简介
HSPF是由美国国家环境保护局于1981年开发的，采用Fortran语言编写的半分布式综合流域模型。模型由透水地段水文水质模拟模块(PERLND)、不透水地段水文水质模拟模块(IMPLND)以及地表水体模拟模块(RCHRES)3部分组成[12]。
2.3 PEST简介与模型率定
PEST自动率定是模型文件输入、输出与用户交互的过程，通过时间序列处理程序TSPROC提取HSPF模型的输入与输出文件，使用非线性评价方法GML进行迭代运算，得到参数之后代入HSPF模型进行模拟，直至满足迭代终止条件，获得最优参数结果[13]。
PEST的目标函数主要涉及了不同时间尺度的径流总量、各种流量组分(基流、洪峰)和不同频率流量保证率的模拟效果[14]，笔者选择日流量偏差、月流量偏差和1%、5%、10%、25%、50%、75%、90%、95%超流量天数偏差作为目标函数，公式如下：
|
$
{f_1}\left( \theta \right) = \sum\limits_{i = 1}^n {{{\{ [{Q_i}-Q{\prime _i}\left( \theta \right)]{w_{1, i}}\} }^2};}
$
|
(1) |
|
$
{f_2}\left( \theta \right) = \sum\limits_{k = 1}^N {{{\left\{ {\left[{\sum\limits_{i = 1}^n {{Q_i}-\sum\limits_{i = 1}^n {Q{\prime _i}\left( \theta \right)} } } \right]{w_{2, y}}} \right\}}^2}} ;{\rm{ }}
$
|
(2) |
|
$
{f_3}\left( \theta \right) = \sum\limits_{i = 1}^{{E_{\rm{c}}}} {{{\{ [{E_i}-E{\prime _i}\left( \theta \right)]{w_{3, i}}\} }^2}} 。$
|
(3) |
式中：f1(θ)、f2(θ)、f3(θ)分别为日流量、月流量、超流量天数偏差的目标函数；Qi和Q′i为日流量的实测值和模拟值；Ei和E′i为日流量实测值和模拟值中超过某一流量值的时间，d；n为模拟时间，d；Ny为模拟月数，Ny=12；n为每月的时间，d；Ec为模拟的超流量天数类型数量；θ为率定参数集；w为目标函数权重。
针对率定中对枯水流量考虑的不足，本研究对日流量偏差f1(θ)这一目标函数进行改进，将时间序列划分为丰水期(4—6月)和枯水期(11月—翌年3月)，f4(θ)和f5(θ)分别计算丰水期和枯水期的日流量偏差，公式如下：
|
$
{f_4}\left( \theta \right) = \sum\limits_{i = 1}^{{n_{\rm{a}}}} {{{\{ [{Q_i}-Q{\prime _i}\left( \theta \right)]{w_{4, i}}\} }^2}} ;
$
|
(4) |
|
$
{f_5}\left( \theta \right) = \sum\limits_{i = 1}^{{n_{\rm{b}}}} {{{\{ [{Q_i}-Q{\prime _i}\left( \theta \right)]{w_{5, i}}\} }^2}} 。$
|
(5) |
式中：na为丰水期模拟时间，d；nb为枯水期模拟时间，d。
笔者选择日流量偏差，即将f1(θ)作为单目标率定的目标函数；F1(θ)和F2(θ)作为2个多目标率定的目标函数，具体公式如下：
|
$
{F_1}\left( \theta \right) = {f_1}\left( \theta \right) + {f_2}\left( \theta \right) + {f_3}\left( \theta \right);
$
|
(6) |
|
$
{F_2}\left( \theta \right) = {f_2}\left( \theta \right) + {f_3}\left( \theta \right) + {f_4}\left( \theta \right) + {f_5}\left( \theta \right)。$
|
(7) |
多目标函数设置不同的权重会直接影响率定结果，应保证各目标函数的偏差在同一个数量级内，本研究将权重设置为各目标函数初始偏差的倒数[15]。为下文方便，本文将基于F1(θ)和F2(θ)多目标函数的分析结果分别称为多目标和多目标(丰枯)的率定结果。
笔者将2001—2005年作为模型的率定期，2006—2010年为验证期，并选取纳什效率系数NSE、决定系数R2、均方根误差RMSE和比例偏差PBIAS对模型日、月径流的模拟结果进行评估。其中NSE和R2表示径流模拟的拟合程度，RMSE和PBIAS均反映模拟值与实测值的偏差。一般认为，若0.5 < NSE≤0.75，模拟结果为良好，NSE>0.75为优秀，R2与NSE相同，PBIAS < 10%为优秀，10%≤PBIAS < 15%为良好，若PBIAS为负值则相反[16]。
3 结果与分析
3.1 参数敏感性分析与参数率定
参考相关研究确定研究区相关模型参数的率定范围[17]，采用Morris筛选法[18]对参数进行敏感性分析，并选取NSE作为评价指标，筛选出9个敏感参数(表 1)。
表 1(Table 1)
表 1 单目标与多目标函数的参数率定结果
Table 1 Estimated parameter values by calibration based on single objective or multi-objective functions
参数
Parameter |
率定范围
Calibration range |
单目标
Single objective |
多目标
Multi-objective |
多目标(丰枯)
Multi-objective (Wet and dry) |
| LZSN/in |
2~14 |
4.2 |
5.39 |
2 |
| INFILT/(in·hr-1) |
0.05~0.6 |
0.340~0.6* |
0.05~0.405* |
0.05~0.6* |
| KVARY/in-1 |
0~1 |
0.00191 |
0.00238 |
0.08282 |
| AGWRC/d-1 |
0.5~0.99 |
0.83 |
0.94 |
0.98 |
| DEEPFR |
0.01~0.5 |
0.393 |
0.01 |
0.011 |
| BASETP |
0.01~0.2 |
0.01 |
0.021 |
0.111 |
| UZSN/in |
0.04~0.6 |
0.049~0.6* |
0.04~0.6* |
0.04~0.137* |
| INTFW |
1~7 |
4.3588 |
7 |
7 |
| IRC/d-1 |
0.2~0.4 |
0.286~0.4* |
0.4 |
0.4 |
| 注：参数使用英制单位。LZSN为额定的下土壤层蓄积；INFILT为土壤渗透率；KVARY为地下水退水系数影响因子；AGWRC为地下水退水系数；DEEPFR为深层地下水入渗系数；BASETP为基流蒸散发系数；UZSN为额定的上土壤层蓄积；INTFW为壤中流入渗系数；IRC为壤中流退水系数。*随土地利用类型变化取值不同。Notes: Parameters use imperial units. LZSN is lower zone nominal storage. INFILT is soil infiltration capacity index. KVARY is a parameter of affecting the behavior of groundwater recession flow. AGWRC is groundwater recession coefficient. DEEPFR is fraction of groundwater inflow lost to deep groundwater. BASETP is fraction of potential ET that can be satisfied from base flow. UZSN is upper zone nominal storage. INTFW is interflow inflow parameter. IRC is interflow recession parameter. * indicates that value varies with the change of land use type. |
|
表 1 单目标与多目标函数的参数率定结果
Table 1 Estimated parameter values by calibration based on single objective or multi-objective functions
|
3.2 模拟结果比较
3.2.1 基于单目标函数与多目标函数参数率定的模拟结果比较
从日径流模拟结果(表 2)来看，在率定期，就R2而言，单目标和多目标率定的结果相差不大；就NSE、RMSE和PBIAS而言，基于多目标率定的结果更优。单目标率定的NSE为0.70，2个多目标率定的NSE均达到0.79，表明2个多目标率定的结果对日径流过程的拟合程度更优。单目标率定的RMSE为32.85m3/s，2个多目标率定分别为27.56和27.42m3/s，单目标率定的PBIAS为-28.4%，2个多目标率定分别为-6.92%和-6.65%，显然，多目标率定可以有效降低日径流模拟值与实测值的偏差。此外，验证期结论与率定期相同，多目标率定的日径流模拟结果更优。
表 2(Table 2)
表 2 单目标与多目标率定的日、月径流过程拟合优度结果
Table 2 Goodness-of-fit of single objective and multi-objective calibration for daily and monthly flows
目标函数
Objective function |
参数名称
Name |
率定期Calibration(2001—2005) |
|
验证期Validation(2006—2010) |
| 日流量Daily flow |
月流量Monthly flow |
日流量Daily flow |
月流量Monthly flow |
| 单目标Single objective |
NSE |
0.70 |
0.81 |
|
0.80 |
0.82 |
| RMSE/(m3·s-1) |
32.85 |
13.17 |
29.99 |
17.52 |
| R2 |
0.78 |
0.96 |
0.88 |
0.98 |
| PBIAS/% |
-28.40 |
-28.42 |
-29.52 |
-29.52 |
| 多目标Multi-objective |
NSE |
0.79 |
0.95 |
|
0.89 |
0.97 |
| RMSE/(m3·s-1) |
27.56 |
6.52 |
21.85 |
7.01 |
| R2 |
0.80 |
0.97 |
0.90 |
0.98 |
| PBIAS/% |
-6.92 |
-6.97 |
-8.71 |
-8.73 |
多目标(丰枯)Multi-objective
(Wet and dry) |
NSE |
0.79 |
0.96 |
|
0.90 |
0.97 |
| RMSE/(m3·s-1) |
27.42 |
5.92 |
21.75 |
7.02 |
| R2 |
0.80 |
0.97 |
0.90 |
0.98 |
| PBIAS/% |
-6.65 |
-6.72 |
-8.10 |
-8.15 |
| 注：NSE为纳什效率系数；RMSE为均方根误差；PBIAS为比例偏差。Notes: NSE is Nash-Sutcliffe efficiency coefficient. RMSE is root mean square error. PBIAS is percent bias. |
|
表 2 单目标与多目标率定的日、月径流过程拟合优度结果
Table 2 Goodness-of-fit of single objective and multi-objective calibration for daily and monthly flows
|
从月径流模拟结果来看(表 2)，NSE、RMSE、R2和PBIAS等4个评价指标的结果与日径流的结论一致，相比单目标率定，2个多目标率定的NSE和R2数值提高，RMSE和PBIAS数值降低，说明2个多目标率定的月径流拟合程度更优，模拟值与实测值的偏差更小。从月径流过程(图 1)也可见，相对于多目标率定而言，单目标率定的月径流对实测值存在更大程度的低估，对径流总量和变化趋势的模拟较差。此外，基于多目标与多目标(丰枯)率定的月径流模拟结果存在差异，并且差异主要出现在枯水期。
从流量历时曲线(图 2)来看，较单目标率定而言，2个多目标率定的FDC曲线整体拟合程度更优，且在低流量部分模拟结果提高的更为显著。为进一步比较各目标函数率定的模拟结果，以Q10、Q50和Q90等流量特征值为代表，统计模拟值与实测值的PBIAS(表 3)。单目标率定的Q10、Q50和Q90的PBIAS在率定期分别为-21.40 %、-25.25%和-54.51%，验证期分别为-30.81%、-23.21%和-43.39%，与实测值相差较大；而多目标率定由于考虑了不同频率的超流量天数偏差，其Q10、Q50和Q90的PBIAS均显著降低，更接近实测值，但Q90仍存在较大的偏差。此外，2个多目标率定的模拟结果也存在较大差异。较多目标率定而言，虽然多目标(丰枯)率定的Q10偏差有所提高，但是Q90的偏差却显著降低，率定期和验证期的PBIAS结果分别从-45.39%和-32.86%提高到-27.03%和-13.20%。显然，多目标(丰枯)率定更能捕捉到低流量的流量特征，可以显著减小枯水流量模拟值与实测值的偏差。
表 3(Table 3)
表 3 单目标与多目标率定的日流量不同频率的PBIAS结果
Table 3 PBIAS of daily flow using single objective and multi-objective calibration at different frequencies
时段
Period |
频率
Frequency/% |
实测值
Observed/(m3·s-1) |
单目标
Single objective |
|
多目标
Multi-objective |
|
多目标(丰枯)
Multi-objective(Wet and dry) |
模拟值Simulated/
(m3·s-1) |
PBIAS/% |
模拟值Simulated/
(m3·s-1) |
PBIAS/% |
模拟值Simulated/
(m3·s-1) |
PBIAS/% |
| 率定期Calibration |
10 |
69.10 |
54.31 |
-21.40 |
|
67.46 |
-2.38 |
|
62.80 |
-9.11 |
| 50 |
22.98 |
17.18 |
-25.25 |
21.83 |
-5.02 |
22.61 |
-1.63 |
| 90 |
9.75 |
4.44 |
-54.51 |
5.32 |
-45.39 |
7.11 |
-27.03 |
| 验证期Validation |
10 |
93.32 |
64.57 |
-30.81 |
|
81.59 |
-12.57 |
|
77.77 |
-16.66 |
| 50 |
27.58 |
21.18 |
-23.21 |
28.24 |
2.38 |
27.50 |
-0.28 |
| 90 |
9.18 |
5.20 |
-43.39 |
6.16 |
-32.86 |
7.97 |
-13.20 |
|
表 3 单目标与多目标率定的日流量不同频率的PBIAS结果
Table 3 PBIAS of daily flow using single objective and multi-objective calibration at different frequencies
|
3.2.2 基于2个多目标函数参数率定的模拟结果比较
针对多目标与多目标(丰枯)率定流量特征值的差异，进一步采用区分丰、枯水期的方式，分析比较2个多目标函数的率定结果。
从丰水期的日、月径流模拟结果(表 4)来看，相比多目标(丰枯)率定，多目标率定的RMSE和PBIAS取得了更好的模拟效果，但是两者在数值上差别较小。同时，从丰水期的月径流过程(图 1)也可以看出，多目标与多目标(丰枯)率定的月径流过程差别不显著。
表 4(Table 4)
表 4 2种多目标率定的丰水期与枯水期日、月径流过程拟合优度结果
Table 4 Goodness-of-fit of two multi-objective calibration approaches for daily and monthly flows during wet and dry seasons
时期
Period |
目标函数
Objective function |
参数名称
Name |
率定期Calibration
(2001—2005) |
|
验证期Validation
(2006—2010) |
| 日流量Daily flow |
月流量Monthly flow |
日流量Daily flow |
月流量Monthly flow |
| 丰水期Wet season |
多目标
Multi-objective |
NSE |
0.78 |
0.96 |
|
0.89 |
0.97 |
| RMSE/(m3·s-1) |
35.09 |
6.15 |
27.55 |
7.89 |
| R2 |
0.79 |
0.96 |
0.89 |
0.98 |
| PBIAS/% |
-0.84 |
-0.78 |
-6.11 |
-6.13 |
多目标(丰枯)
Multi-objective(Wet and dry) |
NSE |
0.78 |
0.96 |
|
0.89 |
0.96 |
| RMSE/(m3·s-1) |
35.18 |
6.28 |
27.61 |
8.65 |
| R2 |
0.78 |
0.96 |
0.89 |
0.98 |
| PBIAS/% |
-2.57 |
-2.55 |
-6.49 |
-6.53 |
| 枯水期Dry season |
多目标
Multi-objective |
NSE |
-0.04 |
-1.03 |
|
0.43 |
0.34 |
| RMSE/(m3·s-1) |
9.54 |
7.00 |
8.89 |
5.54 |
| R2 |
0.37 |
0.49 |
0.59 |
0.78 |
| PBIAS/% |
-33.57 |
-33.73 |
-22.72 |
-22.68 |
多目标(丰枯)
Multi-objective(Wet and dry) |
NSE |
0.25 |
-0.20 |
|
0.54 |
0.71 |
| RMSE/(m3·s-1) |
8.06 |
5.38 |
7.96 |
3.67 |
| R2 |
0.47 |
0.60 |
0.64 |
0.93 |
| PBIAS/% |
-24.52 |
-24.77 |
-16.79 |
-16.81 |
|
表 4 2种多目标率定的丰水期与枯水期日、月径流过程拟合优度结果
Table 4 Goodness-of-fit of two multi-objective calibration approaches for daily and monthly flows during wet and dry seasons
|
从枯水期日径流的模拟结果(表 4)来看，在率定期，相比多目标率定，多目标(丰枯)率定的NSE和R2分别从-0.04提高到0.25、从0.37提高到0.47，说明多目标(丰枯)率定的日径流拟合效果更优；就RMSE和PBIAS而言，相比多目标率定，多目标(丰枯)率定RMSE从9.54降低到8.06m3/s，PBIAS从-33.57%提高到-24.52%，日径流模拟值与实测值的偏差降低。验证期的结论与率定期相同，相比多目标率定，多目标(丰枯)率定的NSE和R2提高，而RMSE和PBIAS降低。以率定期的枯水期日径流过程为例(图 3)，多目标率定明显低估了枯水流量，对枯水期径流的模拟较差；而多目标(丰枯)率定可以捕捉到枯水的流量特征，缩小了模拟值与实测值的差距。月径流NSE、RMSE、R2、PBIAS的结果(表 4)与日径流相同，即划分了丰水期和枯水期之后，HSPF模型月径流拟合程度和流量偏差均得到优化。从枯水期的月径流过程(图 1)也可以看出，相比于多目标率定，多目标(丰枯)率定对于枯水期月径流也有一个更好的把握，拟合效果更优。
4 结论
1) 在水文模型自动率定中，目标函数对模拟结果具有很大影响。相比仅采用日流量偏差的单目标率定，使用日流量、月流量和超流量天数偏差作为目标函数的多目标率定在日、月径流变化趋势、径流总量和不同频率的流量保证率偏差等方面均具有更好的模拟效果；因此，多目标率定综合了不同目标函数的特点，对不同的流量特征均具有较好的模拟效果。
2) 以日流量、月流量和超流量天数偏差作为目标函数的多目标率定，虽然整体模拟效果已经较好，但是对枯水的模拟较差，率定期和验证期的PBIAS分别达到了-33.57%和-22.72%。而将日流量偏差划分丰、枯水期之后，枯水期日流量偏差作为1个单独的目标函数独立出来，并通过权重将4个目标函数设置在同一个数量级内，进而弥补了传统多目标率定忽视枯水流量的不足，显著提高了枯水期日、月径流的模拟效果。
同时，水文模型的参数率定是一个复杂的过程，今后仍需要对目标函数选择及权重设置深入研究，以提高模型的模拟精度及准确性。
2018, Vol. 16 
