检索项 检索词
  土壤学报  2018, Vol. 55 Issue (5): 1108-1119  DOI: 10.11766/trxb201801260022
0

引用本文  

闫一凡, 李晓鹏, 张佳宝, 等. 基于GLUE的土壤溶质运移参数反演及不确定性. 土壤学报, 2018, 55(5): 1108-1119.
YAN Yifan, LI Xiaopeng, ZHANG Jiabao, et al. Parameter Estimation and Uncertainty Evaluation of a Soil Solute Transport Model Using GLUE. Acta Pedologica Sinica, 2018, 55(5): 1108-1119.

基金项目

国家自然科学基金项目(41771265)、国家重点研发计划课题(2016YFD0300601, 2016YFD0200603)和南京土壤研究所“一三五”(领域前沿)项目(ISSASIP1661)资助

通讯作者Corresponding author

刘建立, E-mail: jlliu@issas.ac.cn

作者简介

闫一凡(1989—),女,河南新乡人,博士研究生,主要从事数字农业、土壤水文模型研究。E-mail: yfyan@issas.ac.cn
基于GLUE的土壤溶质运移参数反演及不确定性
闫一凡1,2 , 李晓鹏1 , 张佳宝1 , 刘建立1     
1. 中国科学院南京土壤研究所,南京 210008;
2. 中国科学院大学,北京 100049
摘要:溶质运移模型参数的识别结果常存在较高的不确定性,制约了模型的实际应用。以土壤中Cu2+运移过程为例,采用广义似然不确定性估计(Generalized likelihood uncertainty estimation, GLUE)并引入最大似然值(Maximum Nash-Sutcliffe,MNS)等三种定量指标,探讨了数值反演估计弥散系数等参数的不确定性。结果表明,非线性最小二乘法(Nonlinear least squares,NLLS)得到的唯一“最优”参数组合对Cu2+出流曲线拟合效果很好(R2>0.937),但因“异参同效”,无法刻画预测结果的不确定性。GLUE则可明确溶质运移参数及其响应界面的不确定性,MNS对应的参数组合对Cu2+出流曲线拟合R2>0.937,效果与NLLS的拟合结果高度一致。GLUE计算的95%置信区间覆盖了80%以上的观测点(NLLS为46.3%),其反演参数的取值范围也远大于NLLS的结果。在模型参数及响应界面不确定性分析两方面GLUE方法均优于NLLS方法。
关键词溶质运移模型    参数反演    广义似然不确定性估计(GLUE)方法    不确定性分析    

对流-弥散方程(Convection-dispersion equation,CDE)可用于表征溶质的时空动态特征并揭示溶质在多孔介质中的迁移转化机理,是最常用的描述土壤中溶质运移过程的模型之一[1]。由于溶质运移模型的关键参数(如弥散系数等)存在较强的非线性或尺度依赖性,通常难以通过实验直接准确测定[2]。利用监测数据通过数值模型反演来识别或优化这些参数,是应用CDE解决实际问题的常规做法。传统的参数估计方法主要包括非线性最小二乘法(Nonlinear least squares,NLLS)、高斯-牛顿法(Gauss-Newton)[3]、模拟退火法[4]等。但这些方法容易陷入局部最优,无法应对“异参同效”并进行参数和模拟结果的不确定性分析,在一定程度上制约了模型这一预测和决策工具的潜力发挥[5]

1992年,Beven和Binley[6]提出基于贝叶斯方法的广义似然不确定性估计(Generalized likelihood uncertainty estimation,GLUE)方法。它是一种基于蒙特卡罗模拟,且可识别多种潜在“异参同效”现象的参数优化方法,同时也是不确定性分析的一种重要手段[5]。不同于NLLS这类传统寻优算法,GLUE方法认为存在一个似然度函数,符合一定似然度的参数组合均可接受,且离实测值越近该参数组合的似然度越大,可信度越高。GLUE方法自提出后,已被广泛应用于水文模型、降雨—径流模型以及作物生长等各个领域的模型参数及响应界面的不确定性和敏感性分析中[7-11],在实际应用中表现优良。其似然函数、采样方法及校正标准的选择及其对结果的影响也得到了众多学者的探索[12-15]。与传统寻优算法相比,GLUE方法的优势在于:1)适用于非线性系统;2)无需假设误差的概率分布[12];3)可综合反映所有来源的误差[6];4)可定量分析不确定性。因此,本文选择GLUE方法来探讨数值反演CDE模型参数及模型预测结果的不确定性,期望为数值模拟结果的风险评估提供一定参考。

1 材料与方法 1.1 试验材料

试验所用的土壤采自中国科学院封丘农业生态国家试验站(114.51°E~114.60°E,34.98°N~35.06°N),为黄河冲积成因的典型潮土。自表层向下,土壤质地类型(美国制)分别为砂壤土(0~20 cm)、壤砂土(20~40 cm)和砂黏壤土(40~60 cm),每种质地采集三个重复,共9个土样。实验所用土柱均为高15 cm、直径9 cm的透明有机玻璃柱。采样时,在土柱内壁涂抹少许凡士林润滑,以减少采样阻力、缓减对原状结构的破坏。在采集原状土柱同一地点,取适量扰动土样带回实验室,自然风干后用激光粒度仪(Mastersizer 3000,英国)测定机械组成,烘干法测定容重,电位法测定pH。土壤样本理化性质如表 1

表 1 土壤基本理化性质 Table 1 Physical and chemical properties of soil sample
1.2 静态批量平衡试验及阻滞系数计算

制备5、10、15、20、25 mg·L-1的硝酸铜(Cu(NO3)2)溶液备用。在50 mL的离心管中分别加入2.00 g土壤,加入20 mL制备好的不同浓度梯度的硝酸铜溶液,加塞密封。恒温箱振荡器保持温度23℃±1℃,震荡24 h。然后以4 000 r·min-1的速度离心30 min,过滤分离出上清液,用电感耦合等离子光谱仪(PerkinElmer Optima 8000,美国)测定溶液中的Cu2+浓度。每个浓度处理设置3个重复,并做无土空白实验。利用静态批量平衡实验结果,采用等温吸附模型计算土壤对Cu2+的吸附量,并计算阻滞系数。

采用线性吸附方程,分别对三种质地土壤中Cu2+静态批量平衡试验的数据进行拟合并计算阻滞系数,其公式如下:

$ {\mathit{R}_{\rm{d}}} = 1 + \rho {\mathit{K}_{\rm{d}}}/\theta $ (1)

式中,Rd为阻滞系数;ρ为土壤干容重,g·cm-3Kd为分配系数,L·mg-1θ为土壤体积含水量,cm3·cm-3

1.3 示踪及溶质运移试验

示踪试验在三种不同质地(a砂壤土,b壤砂土,c砂黏壤土)的土柱中进行。试验过程保持室温20℃±1℃。用蠕动泵自下向上缓慢饱和土柱,待形成稳定流场时,输入一个孔隙体积(Pore volume,PV)pH为5.5、浓度为0.05 mol·L-1的溴化钾(KBr)惰性示踪溶液,然后连续洗脱30 h。示踪试验完毕后,以同样的泵入速度向土柱泵入pH为3.5、浓度为0.5 mol·L-1的Cu(NO3)2溶液1 PV,然后连续洗脱30 h。出流溶液用自动部分收集器采集。出流液中Br-浓度由Br-选择电极测定,Cu2+浓度由高效液相色谱-等离子体质谱仪(Agilent 7700x,澳大利亚)测定。

1.4 控制方程

忽略微生物过程,并假定液相中溶质与土壤相间是平衡交换的。考虑线性等温吸附的一维稳态流溶质运移可以用CDE来描述:

$ {\mathit{R}_{\rm{d}}}\frac{{\partial \mathit{C}}}{{\partial \mathit{t}}} = \mathit{D}\frac{{\partial ^2\mathit{C}}}{{\partial {\mathit{x}^{\rm{2}}}}} - \mathit{v}\frac{{\partial \mathit{C}}}{{\partial \mathit{x}}} - \mu \mathit{C} $ (2)

式中,Rd为阻滞系数;C为液相中溶质的浓度,mg·L-1x为距离,cm;t为时间,min;μ为汇项(μ C)的速率系数, min-1,用于描述土壤溶液中基于一阶衰减的不可逆化学持留;v为平均孔隙流速,cm·min-1D为水动力弥散系数,cm2·min-1

1.5 NLLS方法

NLLS算法是以实测值和模型计算值间的误差平方和最小为准则来识别非线性静态模型参数的一种参数估计方法。目前应用最广的基于NLLS的土壤溶质运移参数反演程序是由美国盐土实验室开发的CXTFIT。CXTFIT采用Levenberg- Marquardt算法,通过匹配实测溶质穿透曲线(Breakthrough curve, BTC)来反演弥散系数、孔隙水流速等溶质运移模型参数,并预测溶质浓度随时间和空间的分布规律。其优化目标是寻求一组参数使实测和计算值之间的决定系数R2最大化、残差平方和(SSQ)最小。目标函数为决定系数R2,如下式:

$ {\mathit{R}^{\rm{2}}} = 1 - \frac{{\sum\nolimits_{i = 1}^{\rm{N}} {{{\left( {{\mathit{C}_\mathit{i}} - {\mathit{f}_\mathit{i}}} \right)}^{\rm{2}}}} }}{{\sum\nolimits_{i = 1}^{\rm{N}} {{{\left( {{\mathit{C}_\mathit{i}} - \mathit{\bar C}} \right)}^{\rm{2}}}} }}{\rm{ = }}1{\rm{ - }}\frac{{{\rm{SSQ}}}}{{\sum\nolimits_{i = 1}^{\rm{N}} {{{\left( {{\mathit{C}_\mathit{i}} - \mathit{\bar C}} \right)}^{\rm{2}}}} }} $ (3)

式中,Cifi分别为观测值和模型计算值;C为全部观测值的平均值;n为总观测点数。R2值越接近1,表明拟合效果越好。

1.6 GLUE方法

GLUE方法既考虑到“最优”这一直观事实,又避免了采用单一“最优”参数组合带来的预测风险。考虑到每个参数可能的空间变异和测量误差,将平均孔隙水流速度(v)的初始取值范围定为其测量值ve± 0.5ve cm ·min-1,弥散系数(D)定为0.000 1~0.05 cm2·min-1,阻滞系数(Rd)定为1.0~3.0,汇项速率系数(μ)定为0.001~0.1 min-1。蒙特卡罗采样策略中采用拉丁超立方采样方法(Latin hypercube sampling, LHS)生成随机参数组合。

在溶质运移参数反演时,为避免数值求解困难,一般会利用示踪试验结果先确定vD两个参数,然后将其固定再拟合其他参数。为此,本文中对GLUE方法的应用也分为两个阶段,阶段一:通过对示踪试验数据的拟合确定可接受的参数vD;阶段二:为参数组(vD)随机搭配100组经过LHS采样得到的(Rdμ)重新组合,生成全新的参数组合(vDRdμ)。

为方便与NLLS方法比较,选用Nash–Sutcliffe函数(形同决定系数R2)作为似然函数,定量描述模拟结果与观测值间拟合的优劣程度。

$ \mathit{L}{\left( {{\theta _\mathit{i}}|\mathit{Y}} \right)^{\rm{2}}} = 1 - \frac{{\sum\nolimits_j^n {{{\left( {{\theta _{\mathit{ij}}} - {\theta _{oj}}} \right)}^2}} }}{{\sum\nolimits_j^n {{{\left( {{\theta _{\mathit{oj}}} - {\theta _o}} \right)}^2}} }} = {\mathit{R}^{\rm{2}}} $ (4)

式中,L(θi|Y)为参数组合i的似然函数;θij为参数组合i在时刻j时的模拟值,mg·L-1θoj为时刻j时的观测响应,mg·L-1θo为观测值的平均值,mg·L-1n为总观测点数。对于完美再现观测值的模拟结果,其似然函数值为1。以似然值等于0.9为阈值来区分参数组合为“行为集”(Behavioral set)或“非行为集”(Non-behavioral set)。采用三个指标来定量表征出流液浓度的95%置信区间。分别为平均相对区间长度(Average relative interval length, ARIL)[16],落在95%置信区间内的观测点比例(P95CI)以及最大决定系数值(Maximum determination coefficient),即最大似然值(Maximum Nash-Sutcliffe,MNS)。ARIL表达式如下:

$ {\rm{ARIL = }}\frac{1}{\mathit{n}}\sum {\frac{{{L_{u, t}} - {L_{\mathit{l}, t}}}}{{{R_{o, t}}}}} $ (5)

式中,Lu, tLl, t分别为95%置信区间的出流浓度上下限,mg·L-1n为总观测点数,Ro, t为出流浓度的观测值,mg·L-1

P95CI表达式如下:

$ {{\rm{P}}_{{\rm{95CI}}}}{\rm{ = }}\frac{{{\rm{N}}{{\rm{Q}}_{{\rm{in}}}}}}{n} \times 100\% $ (6)

式中,NQin为落在95%置信区间内的观测点个数,n为总观测点数。

MNS表达式如下:

$ {\rm{MNS = max}}_{i = 1}^{\rm{N}}\left\{ {{\mathit{R}^{\rm{2}}}} \right\}{\rm{ = max}}_{i = 1}^{\rm{N}}{\rm{L}}\left( {{\theta _\mathit{i}}|\mathit{Y}} \right) $ (7)

式中,i为可接受的采样编码,N为可接受的采样数量。模型拟合结果越好,ARIL值越接近于0,P95CI越接近于100%,MNS越接近于1。

2 结果与讨论 2.1 NLLS参数估计及数值模拟

用平衡模型拟合Br-穿透曲线,同时拟合参数vD。由于示踪剂Br-不发生吸附解析、降解沉淀等物理化学反应,其阻滞系数Rd=1,沉淀项速率系数μ=0 min-1。基于NLLS平衡模型对Br-穿透曲线的拟合结果对应的最优参数及决定系数见表 2。由表 2可知,NLLS反演的“最优”参数组合对示踪离子穿透曲线(BTC)拟合效果极佳,R2均大于0.98,均方根误差RMSE<0.046。

表 2 NLLS法对Br-和Cu2+穿透曲线拟合的“最优”参数及拟合效果 Table 2 Optimum parameters of NLLS fitting Br- and Cu2+ BTCs and fitting effect

由于Cu2+运移试验条件均保持与示踪试验时一致,在拟合Cu2+穿透曲线时,仍选用第一阶段识别的vD值,仅对Rdμ进行优化,结果如图 1表 2。观察可知,拟合曲线的峰值略低于观测值,这可能是因为土柱与管壁之间所形成的微弱优势流所致,也可能是实验中存在着一定的测量误差导致的。林青[17]在使用两点非平衡吸附研究重金属运移的过程中也出现类似情况。但总体而言,模拟的结果符合观测值的趋势,R2均大于0.93且RMSE也保持在10-2数量级(表 2)。

图 1 Cu2+穿透曲线(BTCs)观测值与非线性最小二乘法(NLLS)拟合值的对比 Fig. 1 Comparison between measured breakthrough curves (BTCs) of Cu2+ and non-linear least square (NLLS) fitted BTCs
2.2 第一阶段GLUE方法参数估计

以土柱b(壤砂土)的Br-穿透曲线拟合为例,来说明GLUE方法的反演结果。图 2中每一个散点均代表模型运行一次的结果,所以可视作模型响应界面的投影。从指定区域10 000次蒙特卡罗采样中,有818个参数组合的似然值RBr2>0.9,被判定为“行为集”。由图 2表 2可知,最大似然值MNS=0.987 2时的参数(vD)=(0.031 7 cm·min-1,0.003 9 cm2·min-1),与NLLS得到的唯一“最优”解(vD)=(0.032 1 cm·min-1,0.004 2 cm2·min-1),非常接近。且NLLS方法拟合的R2为0.987 4,与最大似然值(MNS)几乎相同。可以预知,蒙特卡罗采样次数越多,模型运行的次数越多,GLUE方法得到的MNS与NLLS寻优方法得到的“最优解”越接近。

图 2 GLUE方法获得的似然值RBr2>0.9的vD散点图 Fig. 2 Dotty plots of v and D of RBr2 > 0.90 acquired by GLUE for bromide BTCs 注:GLUE为广义似然不确定性估计,散点图顶端的线段为NLLS方法的95%置信限。下同 Note:GLUE is the short for generalized likelihood uncertainty estimation and the dashes at the top of the scatter diagram stand for 95% confidence limits determined by NLLS. The same below

与参数的初始取值范围相比,“行为集”的参数取值范围明显缩小。平均孔隙流速v的初始和更新后取值范围分别为[0.015 6,0.046 8] cm·min-1、[0.028 6,0.035 4] cm·min-1;弥散系数D的初始和更新后取值范围分别为[0.000 1,0.050 0] cm2·min-1、[0.000 1,0.022 1] cm2·min-1,区间宽度分别缩小了55.0%和55.8%。但NLLS方法确定的95%置信限(表 2),如图 2中散点图顶端的线段所示,仍较更新后的GLUE的参数取值范围窄很多,说明NLLS方法在Br-运移模型参数反演时存在“异参同效”现象,且其参数vD的置信限远小于实际可被接受的参数范围,导致大量可接受参数被忽略。

2.3 第二阶段GLUE方法参数估计及其不确定性

第二阶段是通过GLUE方法对Cu2+穿透曲线的拟合进行不确定性分析(此处仍以土柱b为例)。为818组RBr2>0.9的“行为集”分别搭配100组通过拉丁超立方采样(LHS)得到的随机参数Rdμ,组成81 800个新的参数组合。逐次运行模型,得到403个参数组合满足Rcu2>0.9,其对应参数的散点图如图 3所示。其他两种质地的土柱也分别得到了448和648组“行为集”。

图 3 GLUE获得的似然值RCu2>0.9的参数的散点图 Fig. 3 Dotty plots of the parameters of RCu2 > 0.90 acquired by GLUE fitting Cu breakthrough

图 3可知,似然值Rcu2最大为0.936 9,对应的参数组合为(v, DRdμ)= (0.034 9 cm·min-1, 0.005 3 cm2·min-1,1.112,0.003 2 min-1)。第二阶段所用的参数vD来自第一阶段对Br-穿透曲线的拟合,因此,这两个参数的可接受范围与图 2一致(无NLLS的95%置信限)。但其在第二阶段的响应界面形状(图 3)与第一阶段不同。这是因为,参数vD之间的交互作用在第二阶段不如第一阶段清晰,第二阶段参数数量增加,彼此之间相互影响,其响应界面也随之改变。参数Rd的可接受范围由[1.0,3.0]缩小为[1.001,1.221],这与根据公式(1),计算得到的静态批量平衡实验的Rd= 1.251也十分接近;参数μ的可接受范围由[0.001, 0.010] min-1缩小为[0.002 6, 0.003 7] min-1,区间宽度分别缩小了89.0%和90.0%。然而,与第一阶段类似,这些可接受的参数范围仍然是通过NLLS方法的95%置信限宽度的5.84倍和5.50倍。这一结果表明,有许多与NLLS寻优得到的“最优解”完全不同的参数组合均可很好地模拟Cu2+的运移过程;换言之,在NLLS方法推算的置信限之外,仍存在众多的参数组合可被接受,并可得到令人满意的模拟结果。而仅通过主观判断,无从确定这些组合中的哪一组参数值才是“真实值”,这也正是模型用于预测和决策时必然存在不确定性的一个重要原因。由此可见,利用NLLS等传统寻优方法得到的“最优解”去预测溶质的运移过程可能存在极大风险,因为其反演结果中蕴含了显著的不确定性。Zhang等[18]用传统寻优法和GLUE方法同时分析了示踪溶液与农药阿特拉津在土壤中的运移,结果也表明在参数估计过程中存在大量的等效参数。

GLUE方法受到争议的原因之一在于选择似然函数和阈值的主观性[19-20]。区分“行为集”和“非行为集”的阈值选择的确存在主观性,但这种选择是建立在对未来模型预测的有效性基础之上的,并非纯粹的主观判断[21]。似然函数的选择亦是如此,本研究中选用的似然函数与NLLS方法的决定系数公式形式相同,目的是便于NLLS和GLUE两种方法的结果对比。近年来,关于如何选定合理的似然函数也得到了更多的关注。如Zhang等[22]为了强调模型在实际应用中的不确定性引入了多种似然函数,并基于概率计算来选定最合理的似然函数。Freni等[15]用多种形式的似然函数估计了城市洪水模型结果的不确定性并对模型进行了校正。这些研究表明,Nash–Sutcliffe函数适用于数据量不大和需要掌握单一参数对模拟结果影响的情况。当观测值不断增多且旨在调整模型适用性时,指数类型的似然函数则更合适。

2.4 模型输出的不确定性

与土柱b类似,土柱a和土柱c示踪和Cu2+运移试验反演结果分别如表 3表 4所示。可以看出,GLUE确定的vDRdμ的置信区间均完全包含NLLS的置信区间。说明在三种土壤中,NLLS的参数置信区间均小于实际可接受的参数范围。若仅使用NLLS算法的置信区间,将导致大量原应被接受的参数组合被舍弃。

表 3 NLLS和GLUE反演得到的Br-运移模型参数及其95%置信区间 Table 3 Parameter of the Br- transport model and 95% confidence intervals obtained using NLLS and GLUE

表 4 NLLS和GLUE反演得到的Cu2+运移模型参数及其95%置信区间 Table 4 Parameter of the Cu2+ transport model and 95% confidence intervals obtained using NLLS and GLUE

GLUE方法对应MNS的参数组合与NLLS“最优”参数组合的模型预测结果如图 4所示。由图可知,在土柱b中,MNS参数组合(vDRdμ)=(0.034 9 cm·min-1,0.005 3 cm2·min-1,1.112,0.003 2 cm-1)与NLLS“最优”参数组合(vDRdμ)=(0.032 1 cm·min-1,0.004 2 cm2·min-1,1.034,0.003 1 cm-1),对观测值的拟合极佳,其R2均为0.937。但MNS与NLLS“最优”的参数组合并不一致,表明不同的参数组合可以达到类似的模拟结果,即“异参同效”现象。在土柱a和土柱c中亦存在相同现象。GLUE方法获取的95%置信区间较NLLS方法获取的置信限宽很多,其覆盖的观测点比例P95CI的均值为84.30%,而NLLS方法则仅为46.05%,这表明NLLS方法获取的置信限无法良好预测Cu2+穿透曲线,特别是峰值区。GLUE方法的95%置信区间对试验观测点的高覆盖率还表明了模型定义及模型结构满足了溶质运移模拟的需求[23]

图 4 NLLS和GLUE两种方法预测Cu2+穿透曲线的不确定性结果对比 Fig. 4 Comparison between NLLS and GLUE used in predicting Cu2+ BTC in uncertainty 注:MNS是最大似然函数,95%CI-NLLS和95%CI-GLUE分别是NLLS和GLUE方法获取的相对浓度的95%置信限 Note:MNS is the short for Maximum Nash-Sutcliffe, which is the largest likelihood value; 95%CI-NLLS and 95%CI-GLUE is separately the 95% confidence interval of relative concentration using NLLS and GLUE

三个定量化指标平均相对区间长度(ARIL)、MNS和置信区间内的观测点比例(P95CI)的统计结果如表 5所示。砂壤土、壤砂土和砂黏壤土中,通过NLLS方法获取的“最优解”的决定系数R2(分别为0.960,0.937,0.954)与GLUE方法获取的MNS的似然值(分别为0.960,0.937,0.953)高度一致,这在图 4也得到直观体现。GLUE与NLLS获得的ARIL均值分别为79.66和135.5,其中壤砂土显著高于其他两种质地土壤。这主要是由于壤砂土内Cu2+穿透曲线的相对浓度在0.000 1以下时,其对应的置信区间和置信限的宽度较其他质地更宽所致。总体而言,GLUE方法对应的ARIL值显著低于NLLS,这说明MNS对应的参数组合对模型的拟合效果可媲美NLLS的“最优”参数组,且GLUE计算的置信区间对观测值的覆盖度更高,平均相对区间宽度也更窄,因此,在参数及模型输出不确定性分析两个方面GLUE方法均优于NLLS方法。

表 5 NLLS和GLUE方法拟合Cu2+穿透曲线的不确定性定量化结果 Table 5 Quantification of the uncertainties in fitting Cu2+ BTCs using NLLS and GLUE

NLLS方法计算得到的“最优”参数组合(表 2)及其对实测值的拟合结果(图 1)均表现优秀。但NLLS方法的残差并非白噪声,这表示该方法的参数估计可能存在偏差[13]。本文的研究结果也证实了该结论。尽管“最优”参数及其置信限的确包含在通过GLUE方法确认的95%置信区间内,但图 2图 3均显示,在NLLS法确定的置信限以外仍有大量参数组合可极好地拟合观测值。这也证实了NLLS寻优的结果——“最优”参数组合并无看似的高“鲁棒性”。再者,NLLS方法对BTC的95%置信限仅覆盖了低于65%的观测点(图 3表 5),该结果再次验证了以上结论。由此可知,试验观测值所能提供的信息可能并不足以识别真实的参数值。在多个可接受的参数向量中,仅通过部分观测值,通常无法判定哪一个参数向量才是参数的真值,对于数据稀缺的野外试验情况更是如此[24]。但依然应该关注和改进传统的寻找单一最优解的算法。其一,对于设计性实验,其理论完备且边界条件等可被严格控制,此时传统寻优算法的单一“最优解”仍有重要意义。其二,GLUE等不确定性分析方法对计算能力要求较高,若计算所有可行性参数组合并运行模型验证超出计算能力,此时就需要单一“最优解”来进行预测。因此,传统寻优算法和GLUE这类承认多解性的算法亦可相辅相成、适当结合。

3 结论

溶质运移模型的参数估计过程和模型预测本身就蕴含不确定性,无论观测数据量多少、数据质量高低均无法完全避免此问题。本文利用传统寻优方法(NLLS)和不确定性分析方法(GLUE)对土壤溶质运移模型的不确定性进行了研究,结果表明,NLLS作为一种参数反演方法简单明了,易于操作,对Br-及Cu2+穿透曲线的拟合效果较好。但“异参同效”现象的存在说明NLLS法的“最优”参数应用于预测时的“鲁棒性”不及预想。GLUE方法清晰直接地表明有多个参数组合能满足拟合要求,其获取的对应于最大似然值MNS的参数向量亦可较好拟合Br-及Cu2+的穿透曲线。由似然值大于0.9的“行为集”计算出各参数的后验取值范围显著大于且包括NLLS方法获取的参数取值范围。NLLS方法较窄的参数置信区间导致许多可被接受的参数组合被摒弃,仅使用单一“最优”解的寻优方法预测溶质运移存在极大的不确定性。

参考文献
[1]
van Genuchten M T, Wagenet R J. 2-site 2-region models for pesticide transport and degradation- theoretical development and analytical solutions . Soil Science Society of America Journal, 1989, 53(5): 1303-1310. DOI:10.2136/sssaj1989.03615995005300050001x (0)
[2]
Elbana T A. Transport and adsorption-desorption of heavy metals in different soils. Alexandria, Egypt: Alexandria University, 2013 (0)
[3]
曾小牛, 刘代志, 牛超, 等. 改进高斯-牛顿法的位场向下延拓. 测绘学报, 2014, 43(1): 37-44.
Zeng X N, Liu D Z, Niu C, et al. Modified Gauss-Newton method for downward continuation of potential field (In Chinese). Acta Geodaetica et Cartographica Sinica, 2014, 43(1): 37-44. (0)
[4]
姚明海, 王娜, 赵连朋. 改进的模拟退火和遗传算法求解TSP问题. 计算机工程与应用, 2013, 49(14): 60-65.
Yao M H, Wang N, Zhao L P. Improved simulated annealing algorithm and genetic algorithm for TSP (In Chinese). Computer Engineering and Applications, 2013, 49(14): 60-65. DOI:10.3778/j.issn.1002-8331.1211-0133 (0)
[5]
Mirzaei M, Huang Y F, EI-Shafie A, et al. Application of the generalized likelihood uncertainty estimation(GLUE)approach for assessing uncertainty in hydrological models:A review . Stochastic Environment Research and Risk Assessment, 2015, 29: 1265-1273. DOI:10.1007/s00477-014-1000-6 (0)
[6]
Beven K, Binley A. The future of distributed models-model calibration and uncertainty prediction . Hydrology Process, 1992, 6(3): 279-298. DOI:10.1002/(ISSN)1099-1085 (0)
[7]
Fraga I, Cea L, Puertas J, et al. Global sensitivity and GLUE-based uncertainty analysis of a 2D-1D dual urban drainage model . Journal of Hydrologic Engineering, 2016, 21(5): 04016004. DOI:10.1061/(ASCE)HE.1943-5584.0001335 (0)
[8]
Dzotsi K A, Basso B, Jones J W. Development, uncertainty and sensitivity analysis of the simple SALUS crop model in DSSAT . Ecological Modelling, 2013, 260: 62-76. DOI:10.1016/j.ecolmodel.2013.03.017 (0)
[9]
Younes A, Mara T A, Fajraoui N, et al. Use of global sensitivity analysis to help assess unsaturated soil hydraulic parameters . Vadose Zone Journal, 2013, 12(1): 1-12. (0)
[10]
Sun M, Zhang X L, Huo Z L, et al. Uncertainty and sensitivity assessments of an agricultural–hydrological model(RZWQM2)using the GLUE method . Journal of Hydrology, 2016, 534: 19-30. DOI:10.1016/j.jhydrol.2015.12.045 (0)
[11]
Sathyamoorthy S, Vogel R M, Chapra S C, et al. Uncertainty and sensitivity analyses using GLUE when modeling inhibition and pharmaceutical cometabolism during nitrification . Environmental Modelling & Software, 2014, 60: 219-227. (0)
[12]
Blasone R S, Vrugt J A, Madsen H, et al. Generalized likelihood uncertainty estimation(GLUE)using adaptive Markov Chain Monte Carlo sampling . Advances in Water Resources, 2008, 31(4): 630-648. DOI:10.1016/j.advwatres.2007.12.003 (0)
[13]
Zhou R R, Li Y, Lu D, et al. An optimization based sampling approach for multiple metrics uncertainty analysis using generalized likelihood uncertainty estimation . Journal of Hydrology, 2016, 540: 274-286. DOI:10.1016/j.jhydrol.2016.06.030 (0)
[14]
Shafii M, Tolson B, Matott L S. Addressing subjective decision-making inherent in GLUE-based multi-criteria rainfall–runoff model calibration . Journal of Hydrology, 2015, 523: 693-705. DOI:10.1016/j.jhydrol.2015.01.051 (0)
[15]
Freni G, Mannina G, Viviani G. Uncertainty in urban stormwater quality modelling:The influence of likelihood measure formulation in the GLUE methodology . Science of the Total Environment, 2009, 408(1): 138-145. DOI:10.1016/j.scitotenv.2009.09.029 (0)
[16]
Jin X, Xu C Y, Zhang Q, et al. Parameter and modeling uncertainty simulated by GLUE and a formal Bayesian method for a conceptual hydrological model . Journal of Hydrology, 2010, 383(3/4): 147-155. (0)
[17]
林青.土壤中重金属运移的数值模拟及不确定性分析.山东青岛: 青岛大学, 2011
Lin Q. Simulation and uncertainty analysis of movement of heavy meatals(In Chinese). Qingdao, Shandong: Qingdao University, 2011 http://www.wanfangdata.com.cn/details/detail.do?_type=degree&id=Y2043387 (0)
[18]
Zhang D, Beven K, Mermoud A. A comparison of non-linear least square and GLUE for model calibration and uncertainty estimation for pesticide transport in soils . Advances in Water Resources, 2006, 29(12): 1924-1933. DOI:10.1016/j.advwatres.2006.02.004 (0)
[19]
Mantovan P, Todini E. Hydrological forecasting uncertainty assessment:Incoherence of the GLUE methodology . Journal of Hydrology, 2006, 330(1/2): 368-381. (0)
[20]
Montanari A. Large sample behaviors of the generalized likelihood uncertainty estimation(GLUE)in assessing the uncertainty of rainfall-runoff simulations . Water Resources Research, 2005, 41(8): 224-236. (0)
[21]
Maulidiani, Rudiyanto, Abas F, et al. Generalized Likelihood Uncertainty Estimation(GLUE)methodology for optimization of extraction in natural products . Food Chemistry, 2018, 250: 37-45. DOI:10.1016/j.foodchem.2018.01.023 (0)
[22]
Zhang Y Q, Liu H H, Houseworth J. Modified generalized likelihood uncertainty estimation(GLUE)methodology for considering the subjectivity of likelihood measure selection . Journal of Hydrologic Engineering, 2011, 16(6): 558-561. DOI:10.1061/(ASCE)HE.1943-5584.0000341 (0)
[23]
Dusek J, Dohnal M, Snehota M, et al. Transport of bromide and pesticides through an undisturbed soil column:A modeling study with global optimization analysis . Journal of Contaminant Hydrology, 2015, 175/176: 1-16. DOI:10.1016/j.jconhyd.2015.02.002 (0)
[24]
Pirot G, Renard P, Huber E, et al. Influence of conceptual model uncertainty on contaminant transport forecasting in braided river aquifers . Journal of Hydrology, 2015, 531: 124-141. DOI:10.1016/j.jhydrol.2015.07.036 (0)
Parameter Estimation and Uncertainty Evaluation of a Soil Solute Transport Model Using GLUE
YAN Yifan1,2 , LI Xiaopeng1 , ZHANG Jiabao1 , LIU Jianli1     
1. Institute of Soil Science, Chinese Academy of Sciences, Nanjing 210008, China;
2. University of Chinese Academy of Sciences, Beijing 100049, China
Abstract: 【Objective】 Computer programs, such as CXTFIT, are commonly used to calibrate soil hydraulic and transport parameters, such as dispersion coefficient and retardation factor. CXTFIT can be used to fit observations quite well, which leave researchers in this aspect such an impression that the "optimum" parameter sets simulated with this program can be used directly for modeling prediction. However, in the process of parameter simulation, inherent uncertainties do exist and are often underestimated. The objectives of this study were to assess and even quantify the uncertainties that may occur in parameter estimation using the convection-dispersion equation (CDE) and in adoption of the parameters in modeling prediction with the non-linear least squares (NLLS) and generalized likelihood uncertainty estimation (GLUE) methods.【Method】 In this study, with the aid of CXTFIT, NLLS and GLUE coupled with the Latin hypercube sampling strategy was used to fit concentrations of bromide and copper nitrate in transport through three oil columns different in texture (i.e. Sandy loam, loamy sand and sandy clay loam), separately. And the parameters were optimized and analyzed to quantify the uncertainties that may occur in these processes by means of three quantitative metrics, that is, MNS (maximum coefficient determination coefficient), P95CI(the percentage of observations included within the 95% confidence intervals) and ARIL (average relative interval length).【Result】 Results show that the only "optimum" parameter set obtained with the NLLS technique fits the curve of solute outflow quite well with determination coefficients (R2) all > 0.98 for fitting Br- transport and > 0.937 for fitting Cu2+ transport, and with root mean square error lingering at the magnitude of 10-2, but it fails to cope with a large number of equivalent parameters. R2 being high in value only indicates the "optimum" parameter set is a proper fit of observation, but it does not mean the "optimum" parameter set is the true characterization of solute transport. The parameter set corresponding to the MNS of solute outflow fitted with can be used to simulate the observation as well as NLLS (R2 > 0.937). But the value range of acceptable parameters determined by GLUE are much wider than that of NLLS (the length of 95% confidence intervals of GLUE is about more than 5 times as high as that of NLLS), which means that a large number of parameter sets that are high in likelihood value fall outside of the 95% confidence intervals determined by NLLS. The 95% confidence intervals of outflow concentration determined by NLLS covered 28.13%, 64.00%, and 46.01% of the data observed separately in the soil columns different in texture, leaving almost half uncovered on average, whereas those determined by GLUE did 87.62%, 80.93%, and 84.3%, separately, which indicates that it is not a good choice to use NLLS to optimize parameters and uncertainties of the model output.【Conclusion】 To put all into a nutshell, GLUE performs better than NLLS in both parameter and response surface uncertainty analysis, for NLLS underestimates significantly the uncertainties in estimation of major transport parameters. GLUE has a much wider acceptable parameter valuation range and 95% confidence intervals for outflow concentration, which indicates that the "optimum" solution acquired by NLLS does not show any robustness as the solution acquired by CXTFIT does. So the usage of only "optimum" parameter sets to predict solute transport has to face high risk and high uncertainty.
Key words: Solute transport model    Parameter estimation    Generalized likelihood uncertainty estimation (GLUE)    Uncertainty analysis