基于GLUE方法的饱和多孔介质中溶质运移模型参数不确定性分析

摘要：模型参数的不确定性分析是模型不确定性研究的重要内容之一。本文以示踪剂Br和反应性溶质Cu在石英砂中的运移为例，采用GLUE方法探讨了多孔介质中溶质运移模型参数的不确定性。研究结果表明，仅对水力学参数θs和λ进行识别时，θs和λ的可识别性较强。对耦合Freundlich等温吸附的模型参数进行识别时，由于参数间的相互作用，θs和λ的可识别性降低；吸附特性参数kF 的后验分布基本呈均匀分布，可识别性较差，吸附特性参数β、ω、f 的取值区间则相对收敛，可识别性较强。K-S 检验结果表明，参数区域敏感度由高到底的排序为f、ω、β、kF、λ、θs，主要是因为石英砂对Cu的吸附以动力学反应为主，而f 和ω是与动力学吸附反应相关的两个参数。上述结论有助于加深对溶质运移模型参数的理解和提高模型预测的可靠性。关键词：GLUE方法；溶质运移模型；参数；不确定性；区域敏感度

中图分类号：O17；X53 文献标识码：A 文章编号：0559-9350（2012）09-1017-08

1 研究背景

随着工农业的迅速发展，土壤污染日益严重，为了更好地管理这些已进入地下环境中的化学物质，使其对地下水的潜在危害降到最低限度，必须准确地描述污染物在地下环境中的物理化学反应和运移过程，由此出现了各式各样的描述多孔介质中溶质运移的数学模型。参数识别是数学模型应用的前提，一个能够反映客观实际的数学模型，只有被赋予合理的参数值，才具有较强的预测能力，因而模型参数识别长期以来都是一个研究的热点问题。

优化方法、区域灵敏度分析（Regionalized Sensitivity Analysis，RSA）方法和普适似然不确定估计（General Likelihood Uncertainty Estimation，GLUE）方法是目前环境数学模型应用较多的参数识别方法[1]。优化方法认为观测数据是系统真值的最好估计，在特定模型结构下只有唯一一组最优参数与之对应。但是由于模型结构复杂性与输入数据的不确定性，仅仅依赖优化方法得到的最优参数进行预测无法保证模型应用的精度和预测结果的可靠性。20世纪70年代末，Hornberger和Spear将过于强硬的优化条件弱化，即将其转化为一些可以用定量或定性的语言描述的条件来决定参数的取舍，从而在一定程度上克服了采用优化方法进行参数识别带来的不确定性问题，这就是RSA 方法[2]。在RSA 方法的基础上，Beven于1992年提出了GLUE方法[3]，它吸收了RSA 方法和模糊数学方法的优点，认为存在一个似然度函数，符合一定似然度的结果都是可接受的，且离实测值越近，似然度越大，可信度越高。GLUE方法既考虑到最优这一直观事实，也避免了采用单一的最优值进行预测而带来的风险。在环境系统分析中，该方法曾被用于城市降雨径流模型、洪水频率模型、地下水补给模型、水质模型及土壤地球化学模型等多个环境系统的模拟预测。Rankinen 等运用GLUE 方法分析了半分布式INCA-N（Integrated Nutrients in Catchments—Nitrogen）模型结构及参数的不确定性[4]。Mertens等在模拟土壤水流动时运用GLUE方法揭示了参数先验分布对模型预测不确定性的影响[5]。卫晓靖等探讨了融合马尔科夫链-蒙特卡洛算法的改进GLUE 方法在SMAR（The Soil Moisture Accounting and Routing Model）模型中的应用，指出该方法能真实合理的反映水文模型的不确定性[6]。近年来，GLUE方法多被应用于流域水文模型的不确定性研究中，但其在土壤溶质运移模拟中的应用还少见报道。已有研究表明土壤中溶质运移的数学模拟也存在着不确定性，故本文试图采用GLUE方法对多孔介质中溶质运移模型参数的不确定性进行分析，以拓展GLUE方法的应用范围，并加深对模型参数的理解和提高模型预测的可靠性。

2 材料与方法

实验以石英砂为介质，其主要矿物成分为SiO2（约占99.5%），另外，还含有氧化铁（<0.02%）、黏、云母和有机杂质等。将石英砂过20目筛，用1∶2的硝酸溶液浸泡24h，然后水洗，最后用去离子水洗至pH 值不再升高，在110℃下烘干。然后把石英砂装入高20cm，内径为5cm的有机玻璃柱中，砂柱分5次装填，每次装填4cm，用塑料压实器压实，以保证砂柱内的全部石英砂颗粒分布均匀。在砂柱两端的接口处各加入一薄层脱脂棉，防止石英砂颗粒堵塞进出水孔。

用蠕动泵由下至上缓慢地注入去离子水饱和砂柱，待水流在砂柱中形成稳定的流场后，输入一个孔隙体积（pv）配制好的溶液，然后用去离子水连续洗脱，整个过程中用自动部分收集器定时采集出流溶液。实验中用到的溶液及研究内容如下：（1）以pH 为5.5，浓度为0.05mol/L的惰性KBr溶液作为示踪剂，以脉冲方式输入1pv的示踪溶液，调节蠕动泵的转速控制流速约为27cm/h，用Br选择电极测量出流液浓度；（2）配制pH为5.5，浓度为5mg/L的Cu（NO3）2 溶液，以脉冲方式输入1pv的溶液。调节蠕动泵的转速控制流速约为27cm/h，用原子吸收分光光度计分析出流液中Cu的浓度。以上实验均设两个重复。文中pv=v*t/l；其中v 为平均孔隙水流速，cm/min；t 为时间，min；l 为土柱长度，cm。

3 溶质运移模型

目前，常用的描述溶质在土壤中运移的模型主要有两大类，即确定性模型和随机模型。确定性模型，即对流-弥散方程（Convection Dispersion Equation，CDE），该模型可以较好地揭示溶质在多孔介质中的运移机理及溶质在时间、空间上的变化特征，具有坚实的物理基础[7-9]。随机模型将溶质在土壤中的运移过程作为一个随机过程来处理，其一是在确定性模型中引入随机参数[10-11]；其二是完全随机模型，也就是说当制约溶质运移的机理不清楚时，在某一时刻从给定土壤容积输出的溶质质量速率可看作是前一时刻输入的随机函数[12]。由于本文研究的是饱和稳定流条件下，溶质在均匀石英砂介质中的运移，故采用CDE模型对实验结果进行模拟分析。

3.1 控制方程

运移方程中包括了水动力弥散和化学非平衡吸附过程，方程如下：

式中：c 为液相中溶质的浓度，mg/cm3；se、sk分别为平衡与非平衡吸附时单位质量石英砂所吸附的溶质质量，mg/g；D为水动力弥散系数，cm2 /min；D=λv ；λ为弥散度（cm）；v 为平均孔隙水流速（cm/min）；ρ为容重，g/cm3；θs为饱和含水量，cm3/cm3；x 为距离，cm；t 为时间，min。

根据实验结果，采用Freundlich等温吸附方程来描述溶质的吸附，则：

式中：f 为吸附平衡时发生瞬时吸附的交换点位所占的分数（-）；ω为一阶动力学速率系数，/min；kF（cm3/g）、β（-）为吸附平衡的经验系数。

3.2 定解条件

实验开始时以脉冲的方式加入一个pv示踪剂和重金属溶液，此时，溶质运移的上边界条件为溶质通量边界，为一常数；随后是一个用去离子水连续洗脱的过程，此时，上边界为零通量边界。在整个模拟过程中，模拟区域的下边界定在溶质的出流端，边界条件为零浓度梯度。上述定解条件的数学表达式为：

式中：q0（0，t）为加入石英砂柱中的水通量，cm/min；对于脉冲阶段，c0（0，t）分别为示踪剂（mol/L）和重金属（mg/cm3）的浓度，对于洗脱阶段，c0（0，t）取值为0。

以上方程通过HYDRUS-1D 求解[13]，获取溶质出流浓度随时间的变化情况，即穿透曲线（BTCs）。

4 GLUE方法

GLUE方法认为模型模拟结果的好坏并不是由模型中的某个参数所决定，而是由一组模型参数来决定。在预先设定的参数分布空间内，通过随机采样的方法获取模型的参数组合并运行模型，为了提高采样效率，本文采用LHS（Latin Hypercube Sampling）方法[14]。选定似然度函数，计算模拟结果与观测值的似然程度。在所有的似然值中，设定一个临界值，低于该临界值的参数组的似然值被赋值为0，对高于该临界值的所有似然值重新归一化，按照归一化处理后似然值的大小，得到参数的后验分布，分析模型参数的不确定性，求出在某置信度下模型预报的不确定性范围。

4.1 似然度函数

似然度函数主要用来表征模拟结果与实测结果的拟合程度，它的选择应该遵循似然函数值随模拟结果相似性增加而单调上升的原则。似然度函数的选取具有一定的主观性，最常用的似然度函数是确定性系数（R2）：

式中：θi 为第i 组参数；Y 为对应参数组θi 的取值；L（θi|Y）为第i 组参数的似然值，此处为确定性系数R2；Qij为第i 组参数在j 时刻的模拟值；Qoj为j 时刻的实测值；Qo 为实测值的平均值；n 为时间序列的数目。R2越接近于1，模拟结果越好，参数的似然值越高。

4.2 参数概率分布

模型参数的不确定性分析，需要确定参数的取值范围和先验分布形式，通常情况下，模型参数的先验分布难以确定。在模型参数的相关研究中，通常采用均匀分布来描述参数的先验分布，本文亦采用均匀分布来描述溶质运移模型参数的先验分布。

模型模拟后参数的后验分布基于贝叶斯理论计算得到。贝叶斯公式的基本形式如下：

式中：P（B_i|A）为参数Bi的后验分布；P（Bi）为参数的先验分布；P（A|Bi）通常为似然度值。

4.3 参数区域敏感度分析

GLUE 方法可利用累积似然度进行区域敏感度分析。如果参数对目标函数没有显著影响那么参数累积概率密度分布应接近于原始分布，即均匀分布；如果参数取值对似然度的影响较大，则参数累积概率密度分布与原始分布相差较大。本文采用Kolmogorov-Smirnov（K-S）统计量衡量似然度累积分布与原始分布的差别，其通过比较2个分布的最大垂直距离来计算2个分布的差异。2个分布之间的K-S距离越大，该参数的区域敏感度也就越高。

D = max|F | n (x) - Fm (x) （7）

式中：D 为2个分布的距离；m、n 为样本的个数；Fn、Fm 为样本分布函数。

5 结果与讨论

5.1 水力学参数的不确定性

选取石英砂中Br运移的实验结果，对模型中的饱和含水量θs 和弥散度λ进行不确定性分析。因为Br为非反应性溶质，故模型中的吸附特性参数kF、β和ω都为0。根据实验条件及模型参数的物理意义，确定参数的取值范围：θs 为0.2~0.8cm3/cm3，λ为0~1cm。通过LHS 方法，得到10 000 组在参数取值范围内均匀分布的参数组，运行模型，共获得1 035 组可接受的参数（R2>0.9），参数及R2的分布如图1所示。

图1 参数-似然值分布散点

经过Br 出流浓度率定后，参数的取值范围缩小，θs 为0.421~0.510cm3/cm3，λ为0.003~0.845cm，显示出明显的峰值区域。另外，在高的似然值区，存在很多参数组，即出现“异参同效（equifinality）”现象，说明存在很多等效的参数组。

基于贝叶斯理论，对于可接受的参数组的确定性系数R2重新归一化，根据式（6）计算各参数的概率密度分布和累积概率密度分布，如图2所示。从中可以看出，饱和含水量θs 和弥散度λ的概率密度均呈类正态分布，参数的可识别性较高，θs 在0.44~0.50 区间的概率密度最大，λ在0.05~0.35 区间的概率密度最大。K-S检验结果表明，θs （0.429）的区域敏感度大于λ（0.199），对饱和多孔介质中非反应性污染物运移的影响较大。

图2 参数后验分布

把可接受的参数组的似然值（R2>0.9）重新归一化，按似然值的大小排序，求出置信度90%下模型预测的不确定性范围。图3给出了模型模拟90%置信度的不确定性范围，其中包括实测值和不确定性的上限和下限。从图中可以看出，除峰值浓度大于1的实测值（由测量误差引起）落在90%的置信区间外，其余的点均落在90%的置信区间内。

图3 Br穿透曲线的实测值和模拟值

5.2 反应性参数的不确定性

在水力学参数θs、λ识别的基础上，选取石英砂中Cu运移的实验结果，对耦合Freundlich吸附的溶质运移模型中的吸附特性参数kF、β、ω、f 加以识别。将经GLUE方法率定后的θs 和λ两个参数的后验取值范围作为新的参数范围，根据实验条件及模型的物理意义，kF、β、ω、f 的取值范围分别为：0~1cm3/g，0.3~1.5，0~0.1/min，0~1。由于一般选取的总的有效模拟次数在一百次到几百次之间[3]，故通过多次模拟运算，决定采用LHS方法产生20 000组在上述参数取值范围内均匀分布的参数组，并以R2>0.7作为临界值，共获得377组可接受的参数组，参数及R2的分布如图4所示。

图4 参数-似然值分布散点

从图中可以看出，θs、λ、kF 的取值均存在很大的不确定性，即这些参数在较大的区间内任意取值，都可能获得较好的模拟结果。θs和λ的取值范围虽然没有发生变化，但是散点图的形状与Br时的明显不同，说明了参数间的相互作用降低了θs和λ的可识别性，故在暂不考虑溶质的吸附作用时，该模型能够比较准确地模拟污染物的运移。β、ω、f 的取值区间则相对收敛，显示出似然值的峰值区域或集中分布区域，说明参数取值的微小变化会引起模型输出的较大差异，因此在模型参数率定阶段应加强对这3个参数的识别。另外，在高的似然值区，存在很多参数组，出现“异参同效”现象，说明存在很多等效的参数组。参数等效性的产生是由于模型参数在高维空间具有的复杂相关性，是模型参数不确定性的集中体现。

基于贝叶斯理论，对可接受的参数组的确定性系数R2重新归一化，根据式（6）计算获得的各参数的概率密度分布和累积概率密度分布如图5所示。θs、λ、kF 的概率密度基本呈均匀分布，可识别性较差；β、ω、f 的概率密度均呈现较明显的类正态分布，可识别型较强。其中β在1.08~1.26区间的概率密度最大，ω在0.005~0.020 区间的概率密度最大，f 在0~0.10 区间的概率密度最大，其中f 在0.6~1.0区间的概率密度较低模型参数识别时可以不予考虑。

对参数进行K-S 检验，检验结果如表1 所示。参数的区域敏感度由高到低排序为f、ω、β、kF、λ、θs，与以往所研究的参数局部敏感度分析结果有较大不同[15]。由于参数局部敏感度分析没有考虑参数间的相互作用[16]，而GLUE 方法中的参数是在封闭的区间随机取值，参数在整个区间的取值不太可能被其它参数同等地抵消，在一定程度上排除了参数相关性的影响[17]，因此区域敏感度分析结果更具有可靠性。θs和λ的区域敏感度比单独识别时有较大幅度的降低，说明参数间的相互作用导致θs 和λ的不确定性增加。同时，θs 和λ的敏感度明显低于吸附特性参数kF、β、ω、f，说明吸附作用对多孔介质中反应性溶质运移的影响较大。f 和ω的敏感度较高可能是因为石英砂对Cu的吸附以动力学反应为主，而f 和ω是与动力学吸附作用相关的两个重要参数，在参数率定过程中应仔细识别。

图5 参数后验分布

表1 可接受参数的K-S统计量

注：θ_s*、λ*为单独识别时的K-S统计量。

图6 Cu穿透曲线的实测值和模拟值

把可接受的参数组的似然值（R2>0.7）重新归一化，按似然值的大小排序，求出在置信度90%下模型预测的不确定性范围，如图6所示，其中包括实测值和不确定性的上限和下限。

从中可以看出，大部分实测值落在90%的置信区间内，只有部分拖尾处的点落在置信区间外，高估了Cu的解吸浓度。造成这种现象的原因可能是：（1）模型本身的问题；（2）简单地认为先验参数为均匀分布；（3）有限的模拟次数。

6 结论

采用GLUE 方法对溶质运移模型参数的不确定性进行了分析，选用LHS 采样方法对参数区间进行均匀采样，以确定性系数R2作为似然度函数，并基于贝叶斯理论得到参数的后验分布，方法直观清晰，表明GLUE方法可以有效地分析多孔介质中溶质运移模型参数的不确定性。

（1）通过对CDE 模型中水力学参数的不确定性分析，发现饱和含水量θs 和弥散度λ可识别性较强，似然值显示出明显的峰值区域，参数的后验分布呈类正态分布。

（2）耦合Freundlich 等温吸附的模型参数的不确定性分析结果表明，由于参数间的相互作用，θs、λ的可识别性降低。吸附特性参数kF 的后验分布基本呈均匀分布，可识别性较差；而β、ω、f 的取值区间相对收敛，具有较强的确定性，对模型输出结果影响较大，在参数率定过程中应仔细识别。K-S检验结果表明，参数区域敏感度由高到底的排序为f、ω、β、kF、λ、θs，这主要是因为石英砂对Cu的吸附以动力学反应为主，而f 和ω是与动力学吸附反应相关的两个重要参数。

（3）计算得到置信度90%的模型预测的上、下限，大部分实测值落在90%的置信区间内，只有部分点落在置信区间外。

参考文献：

[1]邓义祥，王琦，赖斯芸，等. 优化、RSA 和GLUE 方法在非线性环境模型参数识别中的比较[J]. 环境科学，2003，24（6）：9-15 .

[2]Hornberger G M，Spear R C . An approach to the preliminary of environmental systems[J]. Journal of Environmental Management，1981，12：7-18 .

[3]Beven K J，Binley A . The future of distributed models：model calibration and uncertainty prediction[J]. Hydrological Processes，1992，6：279-298 .

[4]Rankinen K，Karvonen T，Butterfield D . An application of the GLUE methodology for estimating the parameters of the INCA-N mode[l J]. Science of the Total Environment，2006，365：123-139 .

[5]Mertens J，Madsen H，Feyen L，Jacques D，Feyen J . Including prior information in the estimation of effectivesoil parameters in unsaturated zone modeling[J]. Journal of Hydrology，2004，294（4）：251-269 .

[6]卫晓靖，熊立华，万民，等. 融合马尔科夫链-蒙特卡洛算法的改进通用似然不确定性估计方法在流域水文模型中的应用[J]. 水利学报，2009，40（4）：464-473 .

[7]Van Genuchten M Th，Wagenet R J . Two-site/two-region models for pesticides transport and degradation：theoretical development and analytical solutions[J]. Soil Science Society of America Journal，1989，53（5）：1303-1310 .

[8]冯绍元，齐志明，王亚平. 排水条件下饱和土壤中镉运移实验及其数值模拟[J]. 水利学报，2004（10）：89-94 .

[9]李志新，许迪，李益农，等. 畦灌施肥地表水流与非饱和土壤水流-溶质运移集成模拟：Ⅰ模型[J]. 水利学报，2009，40（6）：673-678 .

[10]Amoozegar-Fard A，Nielsen D R，Warrick A W . Soil solute concentration distributions for spatially varying pore water velocities and apparent diffusion coefficients[J]. Soil Science Society of America Journal，1982，46：3-8 .

[11]杨金忠，蔡树英，黄冠华，等. 多孔介质中水分及溶质运移的随机理论[M]. 北京：科学出版社，2000 .

[12]Jury W A . Simulation of solute transport using a transfer function model[J]. Water Resources Research，1982，18（2）：363-368 .

[13]?imu。nek J，van Genuchten M Th，?ejna M . The HYDRUS-1D software package for simulating the one-dimensional movement of water，heat，and multiple solutes in variably-saturated media[R]. U . S . Salinity Laboratory，U . S . Department of Agriculture，California，2008 .

[14]McKay M D ，Beckman R J ，Conover W J . A comparison of three methods for selecting values of input variables in the analysis of output from a computer code[J]. Technometrics，1979，21：239-245 .

[15]林青，徐绍辉. 饱和多孔介质中重金属运移参数的敏感度分析[J]. 环境科学学报，2011，31（1）：136-143 .

[16]郑春苗，Bennett G D . 地下水污染物迁移模拟[M]. 北京：高等教育出版社，2009 .

[17]Montanari A . Large sample behaviors of the generalized likelihood uncertainty estimation（GLUE）in assessing the uncertainty of rainfall-runoff simulations[J]. Water Resources Research，2005，41（8）：W08406，doi：10.1029/2004WR003826 .

作者简介：林青（1981-），女，山东人，博士，主要从事地下环境中水流与溶质运移的数值模拟研究。