基于Bayes理论的田间层状土壤水分运动参数识别及不确定性分析

林 青，徐绍辉

（青岛大学 环境科学与工程学院，山东 青岛 266071）

1 研究背景

目前，常用的描述非饱和土壤中水分运动的基本方程是Richards方程，该方程的应用需要获取土壤水分特征曲线及导水率与土壤水势关系的相关信息，由于两者关系的高度非线性，在实验室以及野外直接测量这些属性参数通常是耗时和昂贵的。另外，由于大孔隙、空间异质性、测量尺度等问题的存在，针对点尺度样本的室内稳态试验所得到的土壤水分运动参数往往不能直接反映天然条件下的田间尺度土壤水分运动特征［1］，因此根据原位观测的土壤含水率和土壤水势等数据反演非饱和带土壤水分运动参数的研究一直备受关注。基于最优化理论的反演方法可依据特定的初、边值条件给出最能匹配观测数据的模型参数，为参数的估计提供了一种简单便捷的间接估算方法，很多研究者采用该方法估计了土壤水力特征［2-3］。但是，采用该方法进行反演计算时要求模型满足适定性（可识别性、唯一性和稳定性），其取决于测量变量的个数、边界条件、测量误差及土壤类型等因素［4-5］。虽然已有研究表明，额外变量的测量（如基质势）可以在一定程度上减少不确定性问题［6-7］，但额外变量的测定在实验室内容易实现，而在野外是较难进行的。另外，实际上模型可能有多组参数解，即“异参同效”现象，由于反演问题的不确定性获得的最优解只是这多组参数中的任意一个。因此，最优化反演方法不能很好地处理观测数据和参数的不确定性，存在非唯一解问题；而如果再考虑普遍存在的土壤空间异质性，如层状土壤，将使参数反演问题维数升高，不确定性更为复杂。目前，关于参数不确定性的定量描述及其对模型预测不确定性的影响评价，在国内外已成为研究热点，这有助于我们更为深刻地理解和认识现实世界和模型系统特征之间的本质区别。

近十几年来，贝叶斯方法在模型参数估计及预报中得到广泛应用。但是，应用贝叶斯方法的最大障碍在于后验统计推断中分母积分的计算，尤其是对于具有高维参数向量的复杂总体分布，即便是用数值积分方法，后验分布计算式中分母的计算也比较困难，而迭代逼近算法使这一难题得到了有效解决，其中的马尔科夫链蒙特卡罗（MCMC）方法就是较好的方法。已有研究表明，MCMC法对参数后验分布的搜索，无论是搜索性能还是搜索效率，均表现出了独特的优越性，其在贝叶斯推断的参数识别过程中得到广泛应用［8］。

贝叶斯方法可以最大限度地利用现有资料和已知信息，通过似然函数用观测资料来修正输入的参数分布，使得最终的参数分布与实际相吻合。因此，本文采用Hydrus_1d模拟自然条件下包气带土壤水分运动，基于贝叶斯推断，选取基于DREAMZS（Differential Evolution Adaptive Metropolis）抽样的MCMC方法，对大沽河流域农田层状分布土壤的水力特性参数进行识别，获得参数后验分布的抽样样本，并对其模拟预测的不确定性进行分析，探讨造成模拟结果不确定性的因素，有助于提高模型应用的精度和预测结果的可靠性。

2 研究方法

2.1 试验数据 试验区域位于青岛即墨移风店镇上泊村（36°33'N，120°12'E），地处大沽河中游河谷冲积平原，面积约8 hm2，东西长275.0 m，南北长290.0 m。试验区农作物主要为小麦-玉米轮作，土壤类型为棕壤，土壤结构分层明显，结合地质工程勘察结果，将土壤剖面划分为5层，对应的深度分别为0～30、30～60、60～100、100～200和200～450 cm，其中0～100 cm为耕土层，100～200 cm为粉质黏土层，200～450 cm为中细砂。对耕层土壤用圆盘渗透仪（Disc permeameter）进行饱和导水率Ks的原位测定，人工开挖剖面，用环刀取原状土壤进行水分特征曲线和容重的测定，并取扰动土样进行土壤粒径分析。采用土钻法对土壤含水量进行长期监测，监测点分别设在20、40、80和150 cm处，监测时间从2013年4月1日至2014年10月5日，2013年每7 d监测1次，生育期、雨后或灌溉后加密监测1次，2014年每10 d监测一次。为了降低土壤含水量观测数据的误差，共设置8个观测点（南北各3个，中部2个），取其平均值进行参数反演。地下水位每5 d观测一次。2013年为模型的率定期，利用这一阶段的监测数据进行参数识别和不确定性分析，2014年为验证期，用于预测计算和模型验证。

2.2 作物生长条件下土壤水分运动模型 本文采用经典的Richards方程模拟变饱和土壤中垂向一维水流运动，如式（1）所示，利用Hydrus_1d模型对方程进行求解。

式中：θ为土壤体积含水量，cm3/cm3；h为压力水头，cm；K（h）为非饱和土壤导水率，cm/d；S（z，t）为t时刻z深度处根系吸水速率，表示单位时间单位土体的根系吸水量，cm3/(cm3·d)，采用van Genu⁃chten宏观根系吸水模型计算［9］。其中，t为时间，d；z为土壤深度，坐标向下为正，cm。土壤水分胁迫函数选取S-Shape方程，方程中只有两个参数p0和p50，对冬小麦和夏玉米而言，p0≈3；p50为作物潜在蒸腾率减少50%时相应的土壤基质势，cm。p50与作物生理特性有关，其绝对值越大，作物耐旱吸水能力越强，这里取为-1600 cm。

土壤的持水和导水特征由VGM（van Genuchten-Mualem）模型来计算，

式中：θr、θs分别为残余含水量和饱和含水量，cm3/cm3；α（1/cm）、n和m为曲线形状参数，m=1-1/n；为有效饱和度。在模拟过程中l取经验值0.5。因此，在模型中仅对各层土壤的θr，θs，α，n和Ks进行率定。

模拟期内试验区地下水位变化情况如图1所示，地形标高为17.8 m，地下水埋深为2.7～4.2 m，故取模拟土层厚度为4.5 m，剖分为不等间隔的100个单元。

图1 2013—2014年试验区地下水位

模型初始条件为：

式中：h0（z）为初始压力值，设为-100 cm，目的是为了避免由于层状非均质土壤含水量不连续造成模型不收敛，而且为了降低模拟结果对初始条件的依赖性，在模型校准期设置30 d预热期（spin-up）；L为土体深度，cm。

边界条件为：

式中：hL（t）为下边界处（L=4.5m）的压力水头，模拟期间地下水最大埋深为4.2 m，模拟土层的底部始终位于地下水面以下，故以变水头作为模型的下边界条件；上边界采用通量已知的第二类边界条件，在模拟时段内逐日输入通过上界面的变量值，包括降水量、灌溉量、作物潜在蒸散发量，E（t）表示土壤水分最大蒸发或最大入渗强度，cm/d；ha为地表允许的最小压力水头，设定为-16 000 cm。降雨量数据来自移风镇雨量站，冬小麦和夏玉米的潜在蒸散发率参照FAO推荐的Hargreaves公式来求得，如图2所示，其中作物系数、叶面积指数、根系深度等参照廖凯华在大沽河流域土壤水资源评价研究中的数据［10］。

2.3 基于Bayes理论的参数反演

式中：y为所有观测点的土壤含水量数据；p（y|θ）为给定参数θ时数据y被观测到的概率，即似然函数；p（y）为观测数据的产生概率，该值的计算比较困难，通常的做法是把其从方程中去掉，则：

图2 2013—2014年试验区冬小麦-夏玉米逐日潜在蒸散发率及降雨量

假设各观测点的数据相互独立，则所有数据对应的似然函数是各观测点数据对应的似然函数的乘积，则似然函数可表示为：

式中：n为土壤剖面观测点的个数。为简化计算和达到数值的稳定性，文中似然函数取其对数值。

2.3.2 基于DREAMZS抽样的MCMC模拟 MCMC方法是在贝叶斯理论框架下，建立一个平稳分布为所求后验分布p（θ|y）的随机样本，应用中若产生转移轨迹的马尔科夫链的建议分布选取不适当，会造成该方法收敛速度较慢。DREAM算法作为一种自适应抽样方法，结合了差分进化算法和自适应Me⁃tropolis算法的优点，可同时运行多条马尔科夫链进行全局搜索，并可以自动调整建议分布的范围和方向，在复杂高维、非线性和多峰目标分布问题中表现出极高的效率［13-14］。抽样过程中，假设第i条链的当前状态用d维（待估参数个数）向量θi（i=1，…，N）表示，则其候选样本可根据下式计算：

式中：δ为用于产生候选样本的平行链对数；其中d′为更新的参数维数；r1（j）和r2（h）∈{1，…，N}，且r1（j）≠r2（h）≠i；e和ε为随机产生的很小的数，且e服从均匀分布Ud（-b，b），ε服从于正态分布Nd（0，b*），b、b*为定义的极小值。

2.3.3 模型参数先验分布的确定 本课题组曾对试验区所在大沽河流域上100个样点的土壤饱和导水率Ks进行了测定，统计结果显示其服从对数正态分布。已有研究表明，VGM模型中的形状特性参数α和n分别服从对数正态分布和正态分布，而其中的两个水力特性参数（θs和θr）都遵从正态分布［17-19］，因此，在这里我们设定参数α服从对数正态分布，n、θr和θs服从正态分布，如表1所示。对于0～200 cm深度内的土壤水分运动参数，其先验分布由试验区8个样点32个土样的颗粒组成和容重通过土壤转换函数（Rosetta）获得，由于200～450 cm缺少先验信息，参数设为均匀分布，取值范围根据经验并参考相关文献给出［9，20］。

表1 不同深度处土壤水分运动参数的先验分布和极限值

2.4 不确定性预测区间的评价指标 通常评价预测置信区间优良性指标为预测区间包含的实测数据点比例（PCI）和平均相对区间宽度（ARIL）。PCI值越大，预测区间覆盖实测数据点的比例越高，模型可靠性越大。ARIL值越小，预测的平均相对区间宽度越窄，精度越高。研究表明，随着PCI的增高，ARIL也在增大，即两者往往是相互矛盾的。因此，在前两个指标的基础上，本文还采用单位平均相对宽度所包含的实测点数据比例（PUCI）这个综合指标［21］，来评价模型的可信性，PUCI值越大，置信区间的性能越优良。ARIL和PUCI表达式如下：

式中：LimitUpper，t、LimitLower，t分别为t时刻模拟值95%置信区间的上限和下限；K为观测值的个数；Qobs，t为t时刻的观测值。

3 结果与讨论

3.1 参数的后验分布和统计特征 采用DREAMZS抽样方法，基于计算时间（约5 h）与样本收敛性（R＜1.2）的考虑，设定3条马尔科夫链，每条链最大样本数为50 000个，按1∶10进行稀释。图3为收敛准则R随进化代数的变化曲线，不同颜色的曲线代表了不同参数采样结果的收敛性变化过程。可以看出，当平行马尔科夫链进化代数达到45 000代时，各参数的R值均已小于1.2的基准值，说明各参数都收敛于稳定的后验分布。本文取采样结果最后25%的样本（3750个）作为参数后验分布样本进行统计分析。

图3 DREAMzs算法采样过程的收敛性

表2 参数后验分布的均值、标准差、变异系数、最大后验概率参数（MAP）及95%的置信区间

表2列出了水分运动参数的后验分布特征值及95%的置信区间。从后验均值来看，各层土壤的参数之间差别较明显，尤其是Ks，不同土层之间相差接近10倍，反映出土壤剖面持水与导水特征的变异性。从置信区间和变异系数来看，Ks最不敏感，置信区间较宽，变异系数较大，θs最为敏感，可识别性较高，置信区间较窄，变异系数较小。研究区8个点位耕层土壤Ks的原位分层测定结果分别为5.69～15.12、8.47～23.71和5.91～19.74 cm/d，实测值全部位于95%置信区间内，说明参数Ks后验分布的可靠性较高；但是无论是实测值还是后验分布的标准差都较大，说明Ks的不确定性较大。对比分析室内水分特征曲线反演获得的水分特征参数（如表2试验值所示），发现反演得到的θs与田间获得的参数相差不大，认为可把室内获得的θs用于田间水分的模拟，这与Kumar等的研究结果相似［22］；其它参数则差别较大，这主要是由于田间土壤的大孔隙、空间变异、尺度效应等因素的影响，因此，在使用这些参数对田间土壤水分进行模拟预测时要慎重考虑其适用性。

3.2 最大后验概率参数（MAP）模拟分析 采用后验分布中的最大概率参数组对2013—2014年（2013年为率定期和2014年为验证期）土壤水分动态变化进行模拟分析，不同深度处土壤含水量的模拟值和实测值如图4所示。从图中可以看出，模拟结果基本能反映出各深度处土壤水分的动态变化趋势，相较于率定期而言，验证期的水分变化相对平缓，但是模型仍能较好的模拟预测验证期水分的变化。另外，虽然率定期缺少220 cm深度处的实测含水量数据，但仍能获得该深度处土壤的水力学参数，且验证期采用该参数的模拟值与实测值吻合程度较好。然而，由于土壤含水量受到降水、蒸发等外界因素干扰影响较大，而且土壤水分状况实测值还受到一些不确定因素的影响，加上一些仪器和人为的误差，模拟值和实测值存在一定的偏差，尤其是降雨后土壤水分的实测值和模拟值差别较大。就土壤含水量随土层变化而言，受降雨、灌水及蒸散发的影响，上层土壤含水量变化剧烈，而深层土壤含水量变化相对平缓且接近饱和。

图4 率定期（2013）和验证期（2014）不同深度处土壤含水量模拟值和实测值的对比

采用纳什系数（NSE）、均方根误差（RMSE）、平均绝对误差（MAE）和偏差百分比（PBIAS）4个评价指标来综合衡量模拟的效果。NSE越接近于1，MAE、RMSE和PBIAS越接近于0，说明模拟效果越好，评价指标值如表3所示。

就纳什系数NSE而言，NSE＞0.5时为可接受模拟值。通过计算，不同深度处土壤水分模拟值和实测值的NSE均在可接受的范围内（由于率定期缺乏220 cm深度处的实测数据，故验证期该处的NSE较低为0.498，但是其模拟结果仍能反应土壤水分的变化趋势）。另外，计算得到的RMSE、MAE和PBIAS指标值都较小，其中RMSE小于0.02 cm3/cm3，MAE小于0.015 cm3/cm3，PBIAS处于-1.57%～4.32%之间，接近于0，皆说明模拟值和实测值吻合良好。所以，基于DREAMZS方法的Bayes估计得到的参数可实现对田间层状分布土壤水分动态的良好模拟，尤其是在实测数据不足（缺乏底层土壤含水量实测数据）的情况下，其仍能实现参数的识别及土壤水分的模拟预测。从不同土层来看，20、40、80 cm剖面处模拟值和实测值的NSE相差不大（0.613～0.661），而150 cm深度处的NSE则较小，主要是因为该深度正处于粉质黏土层，黏土水分运动具有明显的滞后性［23］，而模型没有考虑这种滞后性。RMSE、MAE和|PBIAS|3个评价指标都表现出随着土层深度增大而减少的趋势，所以，总体而言模型对深层土壤含水量的模拟精度高于表层，主要原因在于深层土壤含水量受外界条件影响较小，变化相对表层比较平缓。

表3 不同深度处模拟值与实测值吻合度的评价指标值

3.3 模拟结果的不确定性分析 不确定性理论认为，对存在不确定性的事件进行数值模拟时，一条最优的拟合曲线并不能真实反映客观事件本身的规律，科学的模拟应该是确定一组合理的取值区间。因此，针对率定期（2013）模型预测结果的不确定性，用采样的后25%的参数集合运行土壤水分运动模型，根据模拟的结果分析模型预测的参数不确定性及总不确定性。这里参数不确定性指由参数非唯一性所造成的模拟结果的不确定性，主要源自于模型的复杂性、参数间的相互作用等；总不确定性是指由参数的不确定性、模型结构的不确定性及输入数据的不确定性等所引起的模拟预测的不确定性，其由模拟结果加上一个正态分布的随机误差项获得，随机误差的标准差为最大后验概率参数（MAP）经Hydrus_1d模拟得到的均方根误差。图5为95%置信区间下模型预测值的上下限及土壤含水量的实测值，黑色区间为参数不确定性所引起的模型预测结果的不确定性，灰色区间为由参数、模型结构、边界条件等所引起的模型预测的总不确定性，点为实测值。

图5 不同深度处土壤含水量模拟值的95%的置信区间和实测值

从图5可以看出，除了220 cm处的预测区间较宽外，其它深度处的预测区间相对较窄，参数不确定性区域（黑色区间带）没有包含所有观测值，总的不确定性区域（灰色区间带）包含了所有的观测值，但是部分时间段内预测的不确定性范围较大，说明在目前进行土壤水分预测时模型结构、边界条件及输入数据等造成的不确定性是值得考虑的。对比降雨数据（图1），可以看出，模型低估了由于降雨所引起的土壤水分变化预测的不确定性，而高估了无雨期间土壤水分变化的不确定性，表明在进行田间土壤水分模拟的不确定性分析时，采用标准可加误差模型可能是不适当的，这在Jiang等的研究中也有提到［21］。

模型预测置信区间的评价指标值如表4所示。从表中可以看出，由模型参数的不确定性所造成的模拟结果的不确定性较小，预测带平均相对宽度为0.026～0.135，预测带对实测值的覆盖度为0.545～0.773。综合评价指标值PUCI为4.407～26.385，且随着土层深度的增加而增大，表明随着土壤含水量预测深度的增加，模型预测的性能（可靠性和精度的综合表现）越高。模型预测结果的总不确定性明显高于参数预测结果的不确定性，ARIL较大为0.086～0.472，PUCI较小为2.013～11.047，说明模型结构等不确定性造成了很大的模型预测结果的不确定性。相对于模型参数和观测数据而言，基于现有科学认知体系构建的模型结构往往是模拟过程中不确定性的根本来源［24］，所以，未来要提高模型对土壤水分动态模拟的精度，需要对模型结构做进一步的改进。

表4 土壤含水量模拟值95%置信区间评价指标值

4 结论

本文基于贝叶斯理论，采用基于自适应差分演化抽样的MCMC方法，对大沽河流域农田层状分布土壤的水力特性参数进行识别，并分析了土壤水分模拟预测的不确定性，得到了以下几点结论：

（1）采用Bayes方法对模型参数进行识别，不仅得到了确定的模型输出（最大后验概率参数组的模拟结果），而且还获得由参数不确定性和总不确定性所造成的模型预测的区间范围，相比最优化方法获得的一组结果，Bayes方法的模拟结果给出了更多的信息，提高了预测的可靠性。

（2）土壤水分运动参数识别结果显示，Ks最不敏感，置信区间较宽，θs最为敏感，置信区间较窄，可识别性较高。通过与室内试验反演获得的参数相比，发现室内反演得到的θs可用于田间水分的模拟，而其他参数差别较大，在进行田间水分模拟时要谨慎使用。

（3）最大后验概率参数模拟结果基本能反应出各深度处土壤水分的动态变化趋势，即使是在缺失率定期（220 cm深度处）实测数据的情况下，其验证期的模拟值和实测值的吻合程度也较好。所以，基于DREAMZS方法的Bayes估计得到的参数可实现对田间层状分布土壤水分动态的良好模拟。

（4）随着土壤含水量预测深度的增加，PUCI值越大，表明模型预测的性能（可靠性和精度的综合表现）越高，未来要提高模型预测的精度，需要对模型结构做适当的修改完善。

参 考 文 献：

［1］WÖHLING T，VRUGT J A.Combining multiobjective optimization and Bayesian model averaging to calibrate forecast ensembles of soil hydraulic models［J］.Water Resources Research，2008，44（12）：181-198.

［2］查元源，周发超，杨金忠.一种由土壤剖面含水率估算土壤水力参数的方法［J］.水利学报，2011，42（8）：883-891.

［3］李春友，任理，李保国.利用优化方法求算Van Genuchten方程参数［J］.水科学进展，2001，12（4）：473-478.

［4］CARRERA J，NEUMAN S P.Estimation of aquifer parameters under transient and steady state conditions：1.maximum likelihood method incorporating prior information［J］.Water Resources Research，1986，22（2）：199-210.

［5］RUSSO D，BRESLER E，SHANI U，et al.Analysis of infiltration events in relation to determining soil hydraulic properties by inverse problems methodology［J］.Water Resources Research，1991，27（6）：1361-1373.

［6］KOOL J B，PARKER J C.Analysis of the inverse problem for transient unsaturated flow［J］.Water Resources Re⁃search，1988，24：817-830.

［7］PARKER J C，KOOL J B，van GENUCHTEN M T.Determining soil hydraulic properties from one-step outflow experiments by parameter estimation.Part 2：experimental studies［J］.Soil Science Society of America Journal，1985，49：1354-1359.

［8］朱嵩，毛根海，刘国华，等.改进的MCMC方法及其应用［J］.水利学报，2009，40（8）：1019-1023.

［9］ŠIMŮNEK J，ŠEJNA M，SAITO H，et al.The HYDRUS-1D Software Package for Simulating the One-Dimen⁃sional Movement of Water，Heat，and Multiple Solutes in Variably-Saturated Media：version 4.17［EB/OL］.［2017-09-25］.https：//www.pc-progress.com/Downloads/Pgm_hydrus1D/HYDRUS1D-4.17.pdf.

［10］廖凯华.大沽河流域土壤水资源评价及农业节水灌溉模式研究［D］.青岛：青岛大学，2009.

［11］WÖHLING T，VRUGT J A.Multiresponse multilayer vadose zone model calibration using Markov chain Monte Carlo simulation and field water retention data［J］.Water Resources Research，2011，47（4）：W04510，doi：10.1029/2010WR009265.

［12］KAVETSKI D，KUCZERA G，FRANKS S W.Bayesian analysis of input uncertainty in hydrological modeling：1.Theory［J］.Water Resources Research，2006，42（3）：W03407，doi：10.1029/2005WR004368

［13］SCHARNAGL B，VRUGT J A，VEREECKEN H，et al.Information content of incubation experiments for inverse estimation of pools in the Rothamsted carbon model：a Bayesian perspective［J］.Biogeosciences Discussions，2009，6（5）：1736-1751.

［14］SCHARNAGL B，VRUGT J A，VEREECKEN H，et al.Inverse modelling of in situ soil water dynamics：investi⁃gating the effect of different prior distributions of the soil hydraulic parameters［J］.Hydrology&Earth System Sci⁃ences，2011，15（10）：3043-3059.

［15］LALOY E，VRUGT J A.High-dimensional posterior exploration of hydrologic models using multiple-try DREAM（ZS）and high performance computing［J］.Water Resources Research，2012，48（1）：182-205.

［16］GELMAN A，RUBIN D B.Inference from iterative simulation using multiple sequences［J］.Statistical Science，1992，7（4）：457-472.

［17］MOREIRA P H S，van GENUCHTEN M T，ORLANDE H R B，et al.Bayesian estimation of the hydraulic and solute transport properties of a small-scale unsaturated soil column［J］.Journal of Hydrology&Hydromechan⁃ics，2016，64（1）：30-44.

［18］RUSSO D，BOUTON M.Statistical analysis of spatial variability in unsaturated flow parameters［J］.Water Re⁃sources Research，1992，28（7）：1911-1925.

［19］杨金忠.多孔介质中水分及溶质运移的随机理论［M］.北京：科学出版社，2000.

［20］山东省环境水文地质总站.山东省青岛市大沽河地下水库勘察报告［R］.1987.

［21］JIANG S，JOMAA S，BÜTTNER O，et al.Multi-site identification of a distributed hydrological nitrogen model using Bayesian uncertainty analysis［J］.Journal of Hydrology，2015，529（20）：940-950.

［22］KUMAR S，SEKHAR M，REDDY D V，et al.Estimation of soil hydraulic properties and their uncertainty：com⁃parison between laboratory and field experiment［J］.Hydrological Processes，2010，24（23）：3426-3435.

［23］王金生，杨志峰，陈家军，等.包气带土壤水分滞留特征研究［J］.水利学报，2000（2）：1-7.

［24］BECK M B.Water quality modeling：A review of the analysis of uncertainty［J］.Water Resources Research，1987，23（8）：1393-1442.