地下水模拟不确定性问题的多模型分析

宋 凯, 刘 丹, 刘 建

(西南交通大学地球科学与环境工程学院, 四川 成都 610031)

近年来,地下水模拟已成为研究地下水环境各类问题的主要方法,地下水模拟的高度概化、水文地质条件的不协调与研究问题本身的复杂性,致使地下水模拟不确定性问题的出现,直观表征为模拟预报结果与实际情况的偏差.现广泛应用的确定性数值模型仅能获得唯一解,未考虑模拟不确定性对模型预测结果的影响,依据此预测结果进行决策存在风险.因此,十分必要对模型进行不确定分析获取优化模型,提高模型精度.

地下水模拟不确定性根据其来源可分为:参数不确定性、模型不确定性和资料不确定性[1].主要通过参数识别方法研究参数不确定性问题,如单纯形法、最速下降法、共轭梯度法、高斯-牛顿法、遗传算法、模拟退火算法、蒙特卡罗法及贝叶斯方法等[2];一般通过多模型方法探讨概念模型不确定性问题;可通过连续长期观测资料,不依赖于某时刻或单一空间数据来克服资料不确定性.地下水模型参数的不确定性问题己经获得了广泛的关注,如Beven和Binley提出的GLUE (generalized likelihood uncertainty estimation)方法对水文模型的参数不确定性进行估计[3-4];马尔科夫链-蒙特卡洛(Markov Chain Monte Carlo, MCMC)也是一种重要的不确定性分析方法,它在蒙特卡洛模拟框架内不断演化马尔科夫链,使采集的样本参数收敛于模型参数的后验概率分布[5-6],能够有效地探索参数分布空间的高概率密度区域,并反映出参数后验概率的分布特征[7].Hassan等[5]使用MCMC方法对Alaska、Amchitka Island的Milrow试验场地下水模型参数进行了不确定性评价.Kuczera和Parent[8]、Rojas[9]、陆乐[2]和刑贞相[10]等依托不同的试验场地及水文模型研究各类参数不确定性对模型的影响,并应用于模型参数的识别及地下水环境的风险分析.有关地下水概念模型的不确定性分析起步较晚,部分研究针对模型结构的不确定性进行了分析,如Rojas等[11]提出了GLUE与贝叶斯模型平均(Bayesian model averaging)结合的方法,分别对模型参数和概念模型的不确定性进行了统计.此外,Neuman[12]、Ye[13]、吴吉春[1]、曾献奎[14]等亦对概念模型不确定性进行了相关研究.这些研究多侧重于对多个概念模型模拟结果的综合分析,以获取输出变量的分布特征;或假设以理想模型的不同边界条件构建多模型进行影响分析研究.然而,各类不确定因素对模型影响的敏感度及概念模型的可靠性等方面缺乏系统的分析研究.

本文应用调整接受条件后的自适应采样(adaptive metropolis, A-M)算法,以大量水文地质试验数据为参数分布的先验信息构建多模型,根据模型输出数据进行参数识别,并结合基于AICc准则的多模型分析方法,研究参数不确定性及不同结构概念模型对模拟输出的影响及其敏感性.

1 地下水模拟的多模型分析方法

1.1 多模型的构建思路

传统的多模型分析遵循的主要步骤:(1) 考虑构建模拟区多个可能的模型;(2) 在相同观测数据条件下校正模型;(3) 使用某种准则对模型进行排序;(4) 去掉可能性小的模型;(5) 对余下模型得到的预测值与统计量进行权重分析[15-16].

通过调整A-M接受样本条件,对多模型分析步骤进行调整:(1) 将不同概念模型结合服从某分布的随机抽样参数样本,构建研究区多个可能的计算模型;(2) 在相同观测数据条件下,通过调整后A-M采样的样本接受条件,直接剔除预测值与观测值偏差较大模型;(3) 对余下模型得到的预测值与统计量进行权重分析.

1.2 AM-MCMC的优化

MCMC是一种重要的不确定性分析方法,该方法的效率很大程度上取决于其采样的算法.常用的算法有:metropolis-hastings(M-H)算法、吉布斯(Gibbs)采样[17]、A-M算法[18]及single component adaptive metropolis(SCAM)算法[19]等.相比传统的M-H与Gibbs采样,A-M不再需要确定变量的推荐分布,而是决定于初始抽样的协方差,将先验的推荐分布定义为空间的多维正态分布形式,其初始协方差可根据先验信息确定,因此大量先验数据成为A-M采样算法效率及准确性的基础.A-M及SCAM采样原理相近,但若参数组中包含较多维向量,需要分析全局最优解为多维向量参数组时,A-M算法较SCAM更适用.

A-M算法是将参数组看成多维的向量,第i步参数样本推荐服从第i-1次采样所得的向量θi-1为均值,协方差矩阵为Ci的多元正态分布.在初始i0次抽样中,协方差矩阵Ci取固定值C0,C0的确定可依据先验信息,之后自适应更新.协方差矩阵计算如式(1).

(1)

式中:

C0为初始协方差矩阵;

COV(θ0， …，θi-1)为已有的所有样本向量的协方差矩阵;

ε为较小的正数,本次研究取值10-5,为确保Ci不成为奇异矩阵;

sd为比例因子,依赖于参数空间维度d,以确保接受率在一个合适的范围内,sd=2.42/d;Id为d维的单位矩阵[20].

Ci+1为参数i+1次采样的协方差矩阵,由式(1)推出协方差公式为

(2)

式中:

A-M算法采样具体步骤如下:

步骤1按先验分布随机产生初始样本向量θ0;

步骤2利用公式计算Ci;

步骤3产生的参数样本θ*～N(θi,Ci);

步骤4计算接受概率α,

(3)

若接受产生样本θ*,令θi+1=θ*;传统的参数样本接受与否,通过模型的计算值与实际观测资料计算得来的接受概率判定,如式(3).A-M算法不再依赖于参数的推荐分布,假设参数样本先验及后验分布均服从于多元正态分布,以大量实测参数数据为先验信息,随机采集样本的“失真”及后验分布的不收敛基本可忽略.后期依据AICc信息量准则统计模型的预测值进行权重分析,获取最优参数区间.

为提高采样效率,尝试将原接受条件调整为模型的各项计算值与对应观测值的残差均值在K范围内,对模型进行初步筛选,即将步骤4接受样本条件由式(3)调整为式(4).

(4)

式中:

y1、y2分别为实际观测及模型计算值;

n为观测数据个数.

步骤5重复步骤(2)～(4),直到取得足够多的样本.

1.3 基于AICc准则的多模型分析

基于AICc准则的参数识别是利用模型预测结果来计算模型平均预测值及模型残差,并通过排列模型,计算模型概率或权重来分析模型的最优取值区间.单个模型的权重通常是由信息量准则来确定的.信息量准则是信息理论和似然理论结合的成果,日本统计学家赤池弘次首先提出了赤池信息量准则AIC (Akaike’s information criteron),将信息理论中的Kullback-Leibler距离和Fisher极大似然函数联系起来[21].继AIC之后提出修正的AIC信息量准则(AICc).AICc信息量的计算如式(5)～(7).

(5)

(6)

(7)

式中:

Ai为AICc信息量的计算值;

k为待估参数个数.

AICc与AIC的不同仅在于它多了式(5)右边第3项.这项是标识因观测数据较少产生的二阶偏差,当n/k<40时,尤其需要考虑到二阶偏差对多模型分析的影响.在对地下水系统进行建模分析时,n/k<40的情况很普遍[22],因此推荐使用AICc准则.得到AICc值之后,用模型的AICc值减去所有备选模型中的AICc最小值,计算每个模型的Delta值Δi、Δj(i≠j).最后根据Delta值计算模型的后验概率ωi,R是参加多模型分析得到的备选模型总数[23].

(8)

2 实例计算研究

2.1 研究区概况及水文地质条件

研究区位于我国西南某平原区,区内主要分布第四系松散沉积砂砾卵石孔隙潜水含水层.模拟区位于该平原区某河流右岸的一级阶地,地下水类型为第四系全新统冲洪积层砂卵砾石孔隙潜水,根据钻探及模拟区原位水文地质试验成果,含水层厚度约为25 m,渗透系数介于31.34～138.96 m/d,东侧为当地最低侵蚀基准面,地下水由西北向东南径流.区内1995年起运营生活垃圾填埋场,2010年停止堆填并采取封场措施,至今原始地形地貌已然发生改变,根据收集的原始地形资料及现有堆填区实测地形数据分析,堆填区经过削坡、整形等封场措施后仍高于原始地形约9 m.地表高程及坡降等因素的改变将对地下水补给条件产生影响,本文将依据上述差异构建2组不同概念的模型.

2.2 多模型构建

根据原始地形资料及实测的地形数据,构建2组地形存在差异,其余水文地质条件相同的概念模型,1#模型考虑堆填体对原始地形的改变,2#模型为原始地形(图1).依靠先验信息即原位水文地质试验获取的渗透系数数据,生成渗透系数对数的初始平均值及初始方差,依次按AM-MCMC采样方法生成每组样本.每组样本中包含对应网格数的渗透系数对数值,将每组对数值进行转化并输入1#模型进行模拟计算.根据接受条件筛选100组(每组含2 250个数据)参数并获取相应输出数据;将筛选好的100组参数输入2#模型,对应获取相应输出数据;最后依据AICc准则进行多模型分析.

2.2.1模型概化及边界条件设置

根据水文地质条件及钻孔信息构建的2组概念模型范围均为X(1 000 m)×Y(900 m)(X为南北向,Y为东西向),网格为10 m×10 m,含水层厚度约为25 m,东侧边界为当地最低侵蚀基准面,设置为河流边界,西侧及北侧设置为定水头边界,同时,模拟区内设置14个水位观测孔.以原始地形数据及堆填后实测地形数据分别输入模型来刻画2组模型中地形地貌差异(图1).

(a) 1#模型

(b) 2#模型图1 根据不同地形地貌条件概化的2组模型Fig.1 Groups hydrogeological conceptual model

2.2.2AM-MCMC采样

渗透系数为待估参数,视为随机变量.按前述接受条件调整后的A-M算法进行采样.先验信息是参数随机采样的基础,其来源的可靠性将决定采样的合理性.

研究区主要地下水类型为松散堆积孔隙潜水,主要由:山前扇状冲洪积砂砾卵石层孔隙潜水,河道漫滩、一级阶地冲洪积砂卵砾石层孔隙潜水,河间二级阶地冰-水堆积泥质砂砾卵石层孔隙潜水构成,3类孔隙潜水分布于平原坝区,相互叠置,介质类型相似,其间无明显的隔水层,地下水有着密切的水力联系,构成了研究区上部潜水含水层组.

因此,本次先验信息由模拟区周边上述3类含水介质的原位水文地质试验数据构成[24-26],其中包括129个钻孔的174组抽水试验数据(图2)及模拟区5个钻孔的14组数据.经统计,先验信息参数对数取值范围为1.369～5.583,初始参数样本服从均值3.407,协方差C0为0.713的正态分布.依次通过A-M算法采集参数样本,耦合地下水数值模拟软件Modflow输出模拟结果进行不确定性分析.

2.3 地下水流场模拟不确定性分析

文中A-M算法中的接受条件由似然函数求解的后验概率修正为更直接的残差平方的接受值域范围.基于此接受条件的修改,首先需检验条件改变后参数取值的遍历性及收敛性.

注:Q4al+pl为第四系全新统河道漫滩、一级阶地冲洪积砂卵砾石层孔隙潜水;Q4alp为第四系全新统山前扇状冲洪积砂卵砾石层孔隙潜水;Q3fgl+al为第四系上更新统河间二级阶地冰-水堆积泥质砂卵砾石层孔隙潜水;Q1+2fgl+al为第四系中、下更新统泥卵砾石孔隙潜水.图2 研究区同类含水介质水文地质试验孔分布Fig.2 Distribution of boreholes at similar aquifers

2.3.1参数取值遍历性分析

采用A-M算法对参数样本进行采样,根据接受条件共筛选100组参数,每组含2 250个样本值.图3(a)、(b)为参数在采样过程中均值与方差的迭代迹线.当取样20组(样本个数达到45 000个以上)时,参数的均值和方差趋于平稳.图3(c)、(d)分别为45 000个参数样本的采样过程中样本值遍历参数的可能取值范围,通过自适应更新,样本值取样波动逐步减弱,采样过程基本稳定.综合考虑均值、方差迭代迹线和样本采样过程,调整接受条件后的A-M采样方法并没有对参数后验分布的遍历性及收敛性产生影响,更直接的接受条件可修正样本取样空间,提高样本采集效率.

2.3.2多模型分析

多模型的构建旨在分析参数不确定性与模型不确定性对模拟结果的影响.在剔除不符合接受条件的模型后,2组不同的概念模型分别在其中选取100个计算模型.并依据AICc准则多模型分析方法计算得到各模型的AICc值、Delta值及模型的后验概率ωi.

运用AICc准则分析参数不确定性对模拟预测结果的影响时,Burnham和Anderson建议如果模型的后验概率超过0.9时,可视其为最佳模型用于预测.而在地下水模拟不确定性研究中,是不易输出如此高概率模型从而获取单个的最优解,输出的是一系列拟合程度相近的模型,即分析所得的是参数样本的最优取值区间而非唯一解.多模型分析方法可通过对预测平均值或预测值置信区间的分析来反映预测值的范围及其与参数的不确定性关系.运用AICc分析认为Delta值小于2的模型是较好的模型,Delta值介于4～7的模型为经验推荐模型,而Delta值大于10的模型可以舍去[17].甄选1#模型中Delta值小于10的计算模型,对输出观测孔地下水位序列值进行统计分析,可计算各观测孔的众数和置信区间.观测孔水位众数、95%置信区间和实际观测值的序列如图4表示,再进一步剔除Delta值大于10的模型后,剩余计算模型的水位众数与观测值几乎重合.

(a) 参数对数均值(b) 方差(c) 随机样本(d) 相对频率图3 采样过程Fig.3 Sampling process

根据1#模型与2#模型的输出结果排序,排名前三的参数样本相同,排名前十的计算模型中有6组相同,表明在本文设定的不确定性条件下,较之概念模型的不确定性,参数的不确定性的敏感性更高,对模型输出结果的精度更具控制性.虽地下水流场模拟过程中不仅会出现“异参同效”,甚至存在“异参、异构同效”,但仍能从高精度模型出现频次及最优解区间区间范围等方面分析,考虑地形实际变化的1#模型优于2#,反映模型的不确定性仍对模拟结果有着明显的影响.根据输出结果统计分析:(1) 高精度模型出现频次的不同,1#、2#概念模型分别有65%、46%方差值介于1～2之间;(2) 最优解区间取值范围及概率的不同,剔除Delta值大于10的计算模型后,1#模型中仅保留前10个模型,累计后验概率为0.996;2#模型保留前21个模型,而其前10个模型的累计后验概率仅为 0.884(图5).

图4 1#模型计算模型地下水位众数、95%置信区间(阴影区域)与观测值Fig.4 95% confidence intervals (shadow areas), observations and mean values

图5 不同精度模型比例Fig.5 Comparative precision of the different models

3 结 论

地下水模拟不确定性与模型的输入参数、模型结构等因素有关.研究表明:

(1) 由于影响因素间的相互补偿致使模型的输出存在“异参同效”甚至“异参、异构同效”,因此,综合考虑参数和模型结构因素而获取的取值区间应是更合理的.

(2) 在充实的先验数据及参数分布特征既定的条件下,将A-M采样算法中的接受条件调整为模型输出值与实际值方差的接受值域,不会对参数样本的遍历性及收敛性产生影响.

(3) 文中构建的多模型分析方法可识别不同影响因素的敏感性,经过参数样本接受条件及基于AICc准则的多模型分析的双重筛选,能较为高效及准确地获得参数最优区间,同时,亦可完成较优概念模型的甄选识别.

[1] 吴吉春,陆乐. 地下水模拟不确定性分析[J]. 南京大学学报:自然科学,2011,47(3): 227-234.

WU Jichun, LU Le. Uncertainty analysis for groundwater modeling[J]. Journal of Nanjing University: Natural Sciences, 2011, 47(3): 227-234.

[2] 陆乐,吴吉春,陈景雅. 基于贝叶斯方法的水文地质参数识别[J]. 水文地质工程地质,2008(5): 58-63.

LU Le, WU Jichun, CHEN Jingya. Identification of hydrogeological parameters based on the Bayesian method[J]. Hydrogeology and Engineering Geology, 2008(5): 58-63.

[3] BEVEN K, BINLEY A. The future of distributed models-model calibration and uncertainty prediction[J]. Hydrological Processes, 1992, 6(3): 279-98.

[4] BEVEN K, FREER J. Equifinality, data assimilation, and uncertainty estimation in mechanistic modelling of complex environmental systems using the GLUE methodology[J]. Journal of Hydrology, 2001, 249(1): 11-29.

[5] HASSAN A E, BEKHIT H M, CHAPMAN J B. Using Markov Chain Monte Carlo to quantify parameter uncertainty and its effect on predictions of a groundwater flow model[J]. Environmental Moddelling & Software, 2009, 24(6): 749-63.

[6] ROJAS R, KAHUNDE S, PETERS L, et al. Application of a multimodel approach to account for conceptual model and scenario uncertainties in groundwater modelling[J]. Journal of Hydrology, 2010, 394(3): 416-35.

[7] BLASONE R S, VRUGT J A, MADSEN H, et al. Generalized likelihood uncertainty estimation(GLUE) using adaptive Markov Chain Monte Carlo sampling[J]. Advances in Water Resources, 2008, 31(4): 630-48.

[8] KUCZERA G, PARENT E. Monte Carlo assessment of parameter uncertainty in conceptual catchment models: the metropolis algorithm[J]. Journal of Hydrology, 1998, 211(1): 69-85.

[9] ROJAS R, FEYEN L, BATCLAAN O, et al. On the value of conditioning data to reduce conceptual model uncertainty in groundwater modeling[J]. Water Resources Research, 2010, 46: W08520-1-W08520-75.

[10] 刑贞相,芮孝芳,崔海燕,等. 基于AM-MCMC算法的贝叶斯概率洪水预报模型[J]. 水利学报,2007,38(12): 1500-1506.

XING Zhenxiang, RUI Xiaofang, CUI Haiyan, et al. Bayesian probabilistic flood forecasting model based on adaptive metropolis-MCMC algorithm[J]. Journal of Hydraulic Engineering, 2007, 38(12): 1500-1506.

[11] ROJAS R, FEYEN L, DASSARGUES A. Conceptual model uncertainty in groundwater modeling: Combining generalized likelihood uncertainty estimation and Bayesian model averaging[J]. Water Resources Research, 2008, 44: 12418.

[12] NEUMAN S P. Maximum likelihood Bayesian averaging of uncertain model predictions[J]. Stochastic Environmental Research and Risk Assessment, 2003, 17(5): 291-305.

[13] YE M, NEUMAN S P, MEYER P D. Maximum likelihood Bayesian averaging of spatial variability models in unsaturated fractured tuff[J]. Water Resources Research, 2004, 40: W05113-1-W05113-21.

[14] 曾献奎,王栋,吴吉春. 地下水流概念模型的不确定性分析[J]. 南京大学学报:自然科学,2012,48(6): 746-752.

ZENG Xiankui, WANG Dong, WU Jichun. Uncertainty analysis of groundwater flow conceptual model[J]. Journal of Nanjing University: Natural Sciences, 2012, 48(6): 746-753.

[15] NEUMAN S P. Maximum likelihood Bayesian averaging of alternative conceptual mathematical models[J]. Stochastic Environmental Research and Risk Assessment, 2003, 17(5): 291-305.

[16] REFSGAARD J C, SLUIJS J P V D , BROWN J, et al. A framework for dealing with uncertainty due to model structure error[J]. Advances in Water Resources, 2006, 29: 1586-1597.

[17] GILKS W R, RICHARDSON S, SPIEGELHALTER D J. Markov chain monte carlo in practice[M]. London: Chapman & Hall, 1996: 112-119.

[18] HAARIO H, SAKSMAN E, TAMMINEN J. An adaptive metropolis algorithm[J]. Bernoulli, 2001, 7(2): 223-242.

[19] HAARIO H, SAKSMAN E, TANMIINEN J. Componentwise adaptation for high dimensional MCMC[J]. Computational Statistics, 2005, 20(2): 265-273.

[20] GEHNAN A, CARLIN J B, STREN H.S, et al. Bayesian data analysis[M]. London: Chapmann and Hall, 1995: 142-151.

[21] BURNHAM K P, ANDERSON D R. Model selection and multi-model inference: a practical information-theoretic approach[M]. New York: Springer-Verlag, 2002: 163-177.

[22] POETER E P, ANDERSON D. Multi-model ranking and inference in groundwater modeling[J]. Ground Water, 2005, 43(4): 597-605.

[23] 夏强. 地下水不确定性问题的多模型分析方法及应用[D]. 北京:中国地质大学,2011.

[24] 四川省地质局. 成都幅水文地质报告[R]. 成都:四川省地质局,1977.

[25] 四川省地质局. 都江堰幅水文地质报告[R]. 成都:四川省地质局,1977.

[26] 四川省地质矿产局. 成都平原水文地质工程地质综合勘察评价报告[R]. 成都:四川省地质矿产局,1985.