基于PEST的MODFLOW-HYDRUS耦合模型参数优化

胡 丹，伍靖伟，黄介生

(武汉大学水资源与水电工程科学国家重点实验室，武汉 430072)

0 引 言

土壤水和地下水运动是水文循环的重要环节，能够准确模拟土壤水、地下水运动对水资源高效利用具有重要意义[1,2]。地下水与土壤水之间有着直接的水力联系和水通量交换，为了更好地模拟区域尺度土壤水、地下水的运动，近年来发展了一些饱和-非饱和耦合模型，这些模型的主要区别在于对土壤水运动的模拟，其主要方式有3种：第一种是水均衡模型，如MODFLOW中的REC和ET模块[3]，以土壤水质量平衡为原理，独立考虑蒸发、根系吸水等水文过程，由于水均衡模型不涉及水头，无法与地下水的计算直接联系，过分简化了土壤水运动，故会带来较大误差，但求解相对简单；第二种是基于动力波方程，如MODFLOW-UZF1[4]，模拟土壤水运动时只考虑重力作用，而忽略毛管力作用；第三种是基于理查兹方程，如三维非饱和-三维饱和的耦合模型MODFLOW-VSF[5]、一维非饱和-三维饱和的耦合模型MODFLOW-HYDRUS[6]，理查兹方程是由达西定律和质量守恒定律推导而得，对土壤水运动的描述更加准确，但由于具有高度非线性，计算成本也更大。饱和-非饱和耦合模型能否做到模拟精度和计算成本的平衡直接关系到模型的运用推广。Twarakavi等[7]对MODFLOW-UZF1、MODFLOW-VSF和MODFLOW-HYDRUS三种耦合模型进行了比较，发现进行区域模拟时，MODFLOW-HYDRUS和MODFLOW-VSF的模拟精度比MODFLOW-UZF1更高，在计算要求方面，MODFLOW-HYDRUS和MODFLOW-UZF1的计算量比MODFLOW-VSF更小，尤其对于浅地下水、蒸发强烈的地区，MODFLOW-HYDRUS进行土壤水-地下水区域模拟时，相比MODFLOW-UZF1和MODFLOW-VSF，在模拟精度和计算量方面能做到更好的平衡。

饱和-非饱和耦合模型进行区域水分运动模拟时，往往需要输入较多的参数，且水文地质参数和土壤水力参数空间变异性较大，导致模拟结果存在较大不确定性。模型参数的获得目前主要有实测法和参数率定法，其中实验法测得的数据由于尺度效应的存在，很难直接运用到模型中。参数率定法包括试错法、参数自动优化方法以及试错法与参数自动优化法联合优化。传统的试错法依赖使用者对模型的了解及经验，主观性太强，花费时间较长，且难以达到参数最优化。参数自动优化方法能够较快得到优化解，但是对初值的依赖较大，不同的初值可能导致差别较大的率定结果[8]。试错法与参数自动优化法联合优化既能发挥模型使用者的知识和经验，又能利用先进的优化技术[9]。

PEST[10]是独立于计算模型的参数优化方法，已经分别在地下水、土壤水运动模型中得到广泛运用[11-13]，但是在耦合模型中运用效果和效率如何尚未有研究。

本文以内蒙古河套灌区永联试验区地下水观测井数据和土壤含水率实测数据为基础，采用试错法和PEST相结合的方法，对MODFLOW-HYDRUS模型的水文地质参数和土壤水力参数进行优化，以提升区域尺度上饱和-非饱和带中水分运动的模拟效果。

1 材料与方法

1.1 研究区域及数据

研究区域为内蒙古河套灌区五原县永联村境内的永联试验区(见图1)，东经108°00′～108°01′，北纬41°04′～41°05′，南北长约13 km，东西宽约5 km，总面积约2 875 hm2，其中灌溉面积约1 453 hm2，土地利用类型主要是耕地、裸地(包括荒地、渠道、排水沟等)和村庄。试验区地势自南向北稍有倾斜，中间较高，东西两侧略低，整体较为平坦，地面高程1 028.6～1 025.6 m。试验区边界条件较为清楚，北边是六排干排水边界，南边是皂火渠，东西两边分别是永什分干沟和乃永分干沟(见图1)。

试验区内布置有11口地下水观测井(见图1)，每5 d观测1次地下水位，每个地下水观测井附近布置有剖面土壤含水率监测点，每10 d进行1次取样测定，取样深度为0～100 cm，分为0～10，10～30，30～50，50～70，70～100 cm 5个层次。地下水位观测采用测钟法，土壤含水率采用烘干法测定。

图1 研究区示意图Fig.1 Sketch map of study area

1.2 饱和-非饱和耦合模型

非饱和带水流运动考虑为一维垂向运动，用一维Richards方程来描述：

(1)

式中：θ为土壤的体积含水量，cm3/cm3；t为时间，d；h为压力水头，m；z为空间坐标，向上为正，m；S为汇项，m3/(m3·d)；K(h)为非饱和水力传导度，m/d。

饱和带水流运动为三维运动，其控制偏微分方程为：

(2)

式中：h为水头，m；x，y，z为空间坐标，m；Kxx，Kyy，Kzz为渗透系数在x，y，z方向上的分量，这里假定渗透系数主轴方向与坐标轴方向一致，m/d；W为汇项，m3/(m3·d)；t为时间，d；SS为给水度，无量纲。

非饱和带与饱和带的耦合过程参考Seo等[5]的介绍。

本文采用2007年5月1日至11月30日(共214 d)期间实测的地下水位和土壤含水率数据对模型参数进行率定。

(1)网格划分与土壤剖面设置。将研究区域划分为25×50×2网格，其中有效网格为742个，水平方向上网格大小为160 m×264 m，垂直方向上，第一、二层底部高程分别为1 020.29 m和973.75 m。根据地面高程、地下水位、土地利用类型等情况，采用k均值聚类分析方法，在研究区内设置25个土壤剖面(见图2)，其中1～7号土壤剖面为耕地，8～14号为排水沟，15～16号为人民支渠，17号为皂火渠，18～21号为荒地，22～25号为村庄。土壤剖面的高度为地面到1 020.29 m高程处。

图2 网格与土壤剖面Fig.2 Grids and soil profiles

(2)初始输入数据。以2007年5月1日实测的地下水位数据插值得到整个区域的地下水位作为初始地下水位[见图3(b)]；土壤剖面的初始压力水头分布按照Seo等[5]的方法计算得到。

图3 地面高程和初始地下水位Fig. 3 Ground surface elevation and initial groundwater table elevation

(3)边界条件。上边界条件为变化的大气边界和灌溉补给。从五原县气象站获得该研究区域气象资料，包括降雨量、气温、湿度、风速、日照时数等，采用FAO-56[14]介绍的方法计算每日的潜在蒸发量和蒸腾量，将人民支渠口部的净引水量作为区域内的灌溉水量，而灌溉渗漏量概化为人民支渠上的线性补给。耕地(1～7号土壤剖面)的上边界条件包括降雨、蒸发、蒸腾和灌溉[见图4(a)]，荒地(18～21号土壤剖面)、排水沟(8～14号土壤剖面)和皂火渠(17号土壤剖面)的上边界条件包括降雨和蒸发[见图4(b)]，人民支渠(15～16号土壤剖面)的上边界条件包括降雨、蒸发和灌溉渗漏量[见图4(c)]，村庄(22～25号土壤剖面)的上边界条件只有降雨[见图4(d)]。第二含水层下是约1 m的黏土层，故将下边界条件视为零通量边界。将六排干、永什分干和乃永分干三条排水沟设置为排水边界，皂火渠灌溉渠道设置为河流边界。

图4 不同土壤剖面的上边界条件Fig.4 Upper boundary conditions of different soil profiles

(4)模型参数。研究区域两个含水层分别是沙壤土和砂土，水文地质参数的初始值根据五原县永联村钻孔水文地质报告取值(见表1)，van Genuchten[15]模型中土壤水力参数的初始值根据Carsel and Parrish[16]取值(见表1)。然后采用试错法对模型各参数进行调整(见表1)，最后进行参数自动优化。

表1 土壤水力参数和水文地质参数取值Tab.1 Values of soil hydraulic parameters and hydrogeological parameters

1.3 PEST参数优化算法

PEST算法采用高斯牛顿法与梯度下降法结合快速收敛性，又具备梯度下降法的全局搜索性，它能够使目标函数快速有效地收敛，进行参数优化时模型运行的次数比其他算法更少[17]。

PEST算法的核心是对目标函数求解最小值，目标函数是实验观测值与计算值的带权重差异函数，根据参数优化的需求和实验观测值的种类，可将观测值分组并赋予相应的权重，以达到模拟值与实测值的最佳匹配。本研究包含地下水位、土壤含水率2组目标函数，其公式如下：

(3)

式中：Φ为目标函数；Oi和Si分别表示第i个实验观测值和第i个模拟值；p表示实验观测值的个数；ω表示权重系数；下标wt和sm分别对应地下水位和土壤含水率。

1.4 评价方法

采用决定系数(R2)、纳什系数(NSE)、均方根误差(RMSE)、标准化均方根误差(NRMSE)4种指标来评价模型的模拟效果：

(7)

决定系数(R2)用来衡量模拟值和实测值的相关密切程度，R2越接近1表示相关性越好；纳什系数(NSE)可用来评价模型的模拟质量，NSE越接近1说明模拟质量越好；均方根误差(RMSE)用来衡量模拟值与实测值之间的偏差，RMSE越接近0，表示模拟值与观测值的差异越小，模拟效果越好；标准化均方根误差(NRMSE)也是用来衡量模拟误差的大小，NRMSE值越小，说明模拟效果越好。

2 结果与分析

参数自动优化共选取19个参数作为优化对象，包括15个土壤水力参数和4个水文地质参数，参数优化初值为试错法调整后的结果，各个参数的取值范围和优化结果见表1。

将试错法和PEST参数优化后的地下水位模拟值和观测井实测值进行比较(见图5)，同时将模拟开始后第46、86、128、168天土壤剖面含水率的模拟值和实测值进行比较(见图6)。将所有的地下水位观测值与对应时刻的模拟值绘制在一张图中(见图7)，将所有的土壤含水率实测值与对应时刻的模拟值绘制在一张图中(见图7)，并计算评价模型效果的4个指标(见表2)。

采用2008年5月1日至11月31日期间的地下水位和土壤含水率实测值对模型进行验证，将地下水位观测值与对应时刻的模拟值绘制在一张图中(见图8)，将土壤含水率实测值与对应时刻的模拟值绘制在一张图中(见图8)，并计算评价模型效果的4个指标(见表2)。

图5 地下水位模拟值与实测值对比Fig.5 Computed vs. observed groundwater table

图6 不同时刻土壤剖面含水率模拟值与实测值对比Fig.6 Computed vs. observed soil moisture of different time

模拟期观测变量试错法R2NSERMSENRMSEPESTR2NSERMSENRMSE率定期地下水位0．7650．7310．5210．1040．8080．7860．4640．093土壤含水率0．4160．4080．0640．1370．5290．5280．0570．122验证期地下水位0．8040．6030．5800．1290．8240．6410．5520．123土壤含水率0．3250．0320．0660．1740．4380．3950．0530．140

图7 率定期地下水位与土壤含水率散点图Fig. 7 Scatter plots of groundwater table and soil moisture for calibration period

图8 验证期地下水位与土壤含水率散点图Fig. 8 Scatter plots of groundwater table and soil moisture for validation period

从图5可以看出，地下水位模拟值与观测值的变化趋势基本一致；从图6可以看出，模拟开始后不同时刻土壤剖面含水率模拟值与实测值较为接近；从图7、图8可以发现，地下水位和土壤含水率散点都较好地集中在1:1直线附近。

表2的模型效果评价指标显示，PEST算法进行参数优化，提高了模型的模拟效果。

在率定期，地下水位的决定系数由0.765提高到0.808，提高了6%；纳什系数由0.731提高到0.786，提高了8%；均方根误差由0.521减小到0.464，减小了11%；标准化均方根误差由0.104减小到0.093，减小了11%。

在率定期，土壤含水率的决定系数由0.416提高到0.529，提高了27%；纳什系数由0.408提高到0.528，提高了29%；均方根误差由0.064减小到0.057，减小了11%；标准化均方根误差由0.137减小到0.122，减小了11%。

在验证期，地下水位的决定系数由0.804提高到0.824，提高了2%；纳什系数由0.603提高到0.641，提高了6%；均方根误差由0.580减小到0.552，减小了5%；标准化均方根误差由0.129减小到0.123，减小了5%。

在验证期，土壤含水率的决定系数由0.325提高到0.438，提高了35%；纳什系数由0.032提高到0.395，提高了11.3倍；均方根误差由0.066减小到0.053，减小了20%；标准化均方根误差由0.174减小到0.140，减小了20%。

同时从表2可以看出，耦合模型对地下水位的模拟效果要优于土壤含水率。从试错法的结果看，地下水位在率定期和验证期决定系数分别为0.765和0.804，纳什系数分别为0.731和0.603；而土壤含水率在率定期和验证期决定系数分别为0.416和0.325，纳什系数分别为0.408和0.032。通过PEST参数优化后，地下水位在率定期和验证期决定系数分别为0.808和0.824，纳什系数分别为0.786和0.641；而土壤含水率在率定期和验证期决定系数分别为0.529和0.438，纳什系数分别为0.528和0.395。土壤含水率的各项评价指标都低于地下水位的相应指标。对土壤水运动模拟效果较差，可能的原因是本次模拟将整个土壤剖面视为均质土壤，与实际土壤剖面为非均质土壤的情况不符。

3 结 语

(1)本文采用试错法和PEST优化算法相结合的方法对MODFLOW-HYDRUS耦合模型的水文地质参数和土壤水力参数进行优化，提高了模拟精度。与试错法的结果相比，PEST参数优化后决定系数提高了2%～35%，纳什系数提高了6～11.3倍，均方根误差减小了5%～20%，标准化均方根误差减小了5%～20%，其中验证期土壤含水率的各项指标提高或减小的幅度均为最大，决定系数、纳什系数、均方根误差和标准化均方根误差变化幅度分别为35%、11.3倍、20%和20%。

(2)MODFLOW-HYDRUS耦合模型模拟区域尺度上地下水、土壤水运动时，在地下水运动模拟上的表现要优于土壤水运动的模拟。从试错法与PEST参数优化后结果来看，地下水位的决定系数在0.765～0.824之间，纳什系数在0.603～0.786之间；而土壤含水率的决定系数在0.325～0.529之间，纳什系数在0.032～0.528之间，明显低于地下水位相对应的指标。导致土壤水运动模拟效果不佳的原因可能是将整个土壤剖面视为均质土壤，与实际土壤剖面非均质土壤的情况不符，在今后运用模型时可以考虑按照实际的非均质土壤剖面来设置。

□

[1] 杨欣坤, 王 宇, 赵兰坡, 等. 土壤水动力学参数及其影响因素研究进展[J]. 中国农学通报, 2014,30(3):38-43.

[2] 薛禹群. 中国地下水数值模拟的现状与展望[J]. 高校地质学报, 2010,16(1):1-6.

[3] Harbaugh A W, E R Banta, M C Hill, et al. MODFLOW-2000, the U.S. Geological Survey modular ground-water model user guide to modularization concepts and the ground-water fl ow process.USGS, Denver, CO,2000.

[4] Niswonger R G, D E Prudic, R S Regan. Documentation of the Unsaturated-Zone Flow (UZF1) package for modeling unsaturated flow between the land surface and the water table with MODFLOW-2005. Techniques and Methods 6-A19. Available at http:∥pubs.usgs.gov/tm/2006/tm6a19/ (verified 14 Apr. 2008). USGS, Denver, CO,2006.

[5] Thoms R B, R L Johnson, R W Healy. User’s guide to the Variably Saturated Flow (VSF) process for MODFLOW. Techniques and Methods 6-A18. Available at http://pubs.usgs.gov/tm/2006/ tm6a18/ (verified 15 Apr. 2008). USGS, Reston, VA, 2006.

[8] 舒晓娟, 陈洋波, 黄锋华, 等. PEST在WetSpa分布式水文模型参数率定中的应用[J]. 水文, 2009,29(5):45-49.

[9] 章四龙, 刘九夫. 通用模型参数率定技术研究[J]. 水文, 2005,25(1):9-12.

[10] Doherty J. PEST Model-Independent Parameter Estimation User Manual[M]. 5th Edition. Brisbane, Australia: Watermark Numerical Computing, 2004.

[11] Zyvoloski G, E Kwicklis, A A. Eddebbarh, et al. The site-scale saturated zone flow model for Yucca Mountain: calibration of different conceptual models and their impact on flow paths[J]. Journal of Contaminant Hydrology, 2003,62-63(2):731-750.

[12] 董艳辉, 李国敏, 郭永海, 等. 应用并行PEST 算法优化地下水模型参数[J]. 工程地质学报, 2010,18(1):140-145.

[13] Fang Q X, Green T R, Ma L, et al. Optimizing soil hydraulic parameters in RZWQM2 under fallow conditions[J]. Soil Science Society of America Journal, 2010,74(6):1 897-1 913.

[14] Allen R G, Pereira L S, Raes D, et al. Crop Evapotranspiration: Guidelines for Computing Crop Water Requirements[D]. Rome, Italy: United Nations Food and Agriculture Organization, 1998.

[15] van Genuchten, M.Th. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[J]. Soil Science Society of America Journal, 1980,44(44):892-898.

[16] Carsel R F, Parrish R S. Developing joint probability distributions of soil water retention characteristics[J]. Water Resources Research, 1988,24(5):755-769.

[17] Bahremand A, De Smedt F. Distributed hydrological modeling and sensitivity analysis in Torysa Watershed, Slovakia[J]. Water Resources Management, 2008,22(3):393-408.