贝叶斯多模型分析方法的对比研究

薛　亮,夏　强

(1.中国石油大学石油工程学院油气田开发工程系,北京102249；2.成都理工大学地质灾害防治与地质环境保护国家重点实验室,四川成都610059；3.油气资源与探测国家重点实验室,北京102249)



贝叶斯多模型分析方法的对比研究

薛亮1,3,夏强2

(1.中国石油大学石油工程学院油气田开发工程系,北京102249；2.成都理工大学地质灾害防治与地质环境保护国家重点实验室,四川成都610059；3.油气资源与探测国家重点实验室,北京102249)

摘要：多模型分析能够考虑模型本身存在的不确定性,在决策分析和风险评估中具有越来越重要的作用。对具有严谨统计分析理论基础的贝叶斯模型平均法和极大似然贝叶斯多模型平均法做了详细介绍,并改进了传统极大似然贝叶斯多模型平均法不能考虑参数不确定性的不足,使极大似然贝叶斯多模型平均法对贝叶斯模型平均法近似得更为准确。地质统计多模型对地层渗透系数的预测分析结果表明,2种多模型分析方法在参数空间确定、模型后验权重和渗透系数预测方面都具有很好的一致性。极大似然贝叶斯多模型平均法能和水文学现有的参数估计方法很好结合,且计算量小。

关键词：不确定性分析；贝叶斯模型平均法；极大似然估计；地质统计

0引言

地下水资源管理及污染防治一直以来都是国内外学者关心的热点问题。合理的地下水资源管理规划和有效的污染防治方案制定都需要建立相应的模型,以对地下水的流动动态和污染物运移的趋势做出相关预测,并依据预测结果制定相应的决策。然而,地下水系统是一个复杂的自然系统,使用单一概念的模型对地下水系统进行描述,在解决决策分析问题中存在着诸多不足。Neuman在多模型分析研究中提出,使用单一概念模型存在2种类型的错误：一是由于依赖不准确模型而造成的统计偏差；二是由于模型空间考虑不充分而造成决策风险低估[1]。多模型分析研究已经成为国内外水文领域研究的前沿热点问题[2- 6]。

多模型分析方法主要有3种：广义似然不确定性估计法(GeneralizedLikelihoodUncertaintyEstimation,GLUE)、贝叶斯模型平均法(BayesianModelAveraging,BMA)和极大似然贝叶斯模型平均法(MaximumLikelihoodBayesianModelAveraging,MLBMA)。GLUE最早在水文学中提出,现已成为一种常用的多模型分析方法[7]。然而,GLUE通常被认为缺乏严格的统计理论依据,似然方程和模型保留准则需主观确定,参数与预测的统计性质不统一,没有在预测概率分布中明确考虑模型误差等[8-9]。BMA由统计学家在贝叶斯分析框架下提出,理论方法严谨,能考虑模型、参数、观测等诸多因素导致的不确定性,应用于诸多领域。Hoeting等人对BMA做了详细的研究和探讨[10]。为使BMA更易于与水文地质领域完善的极大似然模型反演方法相结合,Neuman建立了MLBMA[1],并证明了多模型方法相较于单一模型在解决统计决策分析问题上的优势[2],但该方法忽略了参数的不确定性。

本文对传统的MLBMA进行了改进,考虑参数的不确定性,使MLBMA在考虑不确定性因素方面获得和BMA相对等的功能。对BMA和MLBMA理论方法分别进行介绍,对比2种方法的实现特点和执行流程,并通过地质统计多模型分析，对2种多模型方法的分析结果进行比较。

1贝叶斯多模型分析方法

1.1贝叶斯模型平均法(BMA)

预测量的多模型后验概率分布是组成模型集的单个模型给出的预测量后验概率分布的加权平均，进行加权平均时的权重为该给定单个模型的后验概率。假设预测量为Δ,K个模型都可以对预测量进行估计预测,这些模型组成模型集合M={M1,…,MK},模型集中的每个模型Mk都由参数向量θk进行表征,观测数据为D,则条件于该观测数据的预测量多模型后验概率分布可由下式定义

(1)

考虑参数存在不确定性时,给定模型的预测量后验概率为

(2)

模型后验概率的确定是多模型分析的关键点,可以通过贝叶斯定理进行计算，公式为

(3)

边际似然函数可以通过下式进行计算

(4)

得到模型的后验概率后,预测量的多模型均值为

(5)

(6)

(7)

(8)

BMA多模型分析方法的执行步骤为：①定义每个模型中参数的先验概率分布,并随机取样产生1组参数实现；②对所有模型中每组参数的每个实现对应的似然函数进行计算；③通过式(4)计算边际似然函数,并计算考虑参数不确定性的预测均值和方差；④通过式(5)、(6)计算考虑模型不确定性的均值和方差。

1.2极大似然贝叶斯模型平均法(MLBMA)

对于每个给定的模型而言,在定义模型参数的概率密度函数并据此生成模型参数实现的过程中,并没有考虑观测数据对参数分布的影响,参数的分布是依据经验主观定义的,这会导致由于有效计算量的减少造成总体计算量的增加。如果参数分布定义不准确,还可能导致最终的统计分析结果存在偏差。此外,式(3)、(4)中积分的准确计算,都需要通过模特卡罗方法对每个模型的参数进行大量取样,这也会造成计算量的增加。

(9)

式(2)可由按照此概率分布产生的参数实现,通过模特卡罗积分进行计算。

在极大似然估计理论框架下,模型的后验概率或者模型权重可以通过下式进行估计

(10)

KIC=-2ln[

(11)

MLBMA多模型分析方法的执行步骤为：①通过极大似然估计方法得到每个模型的参数估计值及其估计协方差,并得到相应的似然函数值；②依据参数的多元正态分布式(9)随机取样产生1组参数实现,并计算考虑参数不确定性的预测均值和方差；③通过式(10)计算后验模型概率；④通过式(5)、(6)计算考虑模型不确定性的均值和方差。

2算例研究

为对比BMA和MLBMA这2种多模型方法,本文使用地质统计模型对地层渗透系数的空间非均质性进行分析。二维50×50单位长度的参考或真实渗透系数场由高斯型截断幂函数地质统计模型(TruncatedPowerVariogramwithGaussianmode,TPVG)产生,可表示为

(12)

生成的渗透系数参考场见图1。生成参考场模型的参数值为θ=(A,H,λu)T=(0.1,0.25,5)T。100个随机选择数据点上的渗透系数值被用作观测数据进行多模型分析,其空间位置见图1中的“D”字符所示。

图1　生成的渗透系数参考场

(13)

球状模型：γ(s)=

式中,σ2为渗透系数场的方差；R为渗透系数场的变程。

此次研究选取了3个参数的截断幂函数地质统计模型作为生成参考场的真实模型,比包含2个参数的指数、高斯和球状模型更为复杂,并能更灵活地刻画渗透系数场的空间变异性。另外,为使构建的算例能反映真实的模拟需要,生成参考场的截断幂函数模型并没有放到进行多模型分析的模型集合内。因为通常在实际情况下,真实模型一般都比选择的模型复杂,且真实模型都是未知的。

在执行BMA分析时,假设这3种选择模型的参数服从多元均匀分布,方差的取值范围为0～3,变程为0～30,通过拉丁超立方方法产生10 000组参数集合,并对每组参数依据观测数据计算相应的似然函数值和参数权重。在执行MLBMA时,3种选择模型都通过100个观测数据进行参数校正,得到参数估计值和估计协方差,假设校正后的参数服从多元正态分布,并通过式(9)产生2 000组参数实现。2种方法的数值计算结果都是通过MATLAB自行编制代码,结合模型参数校正软件PEST计算获得。3种模型在BMA和MLBMA中的参数对比结果见图2。图2中,黑点为在MLBMA中产生的2 000组方差和变程组成的参数,背景等值线为在BMA中10 000组方差和变程参数对应的插值参数权重。

从图2可以看出,MLBMA产生的参数集与BMA具有很强的一致性,主要表现在MLBMA产生的参数区域能够很好地覆盖BMA具有较高权重的参数区域。BMA首先通过先验的猜测选定参数取值区域,然后通过观测数据计算似然函数值,并得到该参数的权重。只有当所有的生成参数都计算完毕后,才能得到在加权平均中的起重要作用的有效参数。在此例中,有效参数主要存在于等值线图的黑色区域,而灰色区域的参数由于权重很低,在得到渗透系数预测统计结果的加权平均过程中起不到重要作用。因此,BMA要求的参数取样量较大(此处为10 000组),这也造成了计算量相对较高。MLBMA通过参数校正过程,利用观测数据首先确定参数产生的区域。由于多元正态分布的假设,其在极大似然估计参数值附近区域会产生较多的取样点,随着向区域外部的扩展,取样点数随之较少,每个参数取样点处的渗透系数预测计算结果都是等权重平均。因此,此处2 000组参数在计算最终渗透系数预测统计结果时起的作用都是一样的,避免了无效的预测计算,节省了计算量。

图2　3种模型在BMA和MLBMA中的参数对比

后验模型概率或者模型权重是多模型分析的重要基础。BMA和MLBMA的模型权重对比见图3。从图3可以看出,2个多模型分析方法计算的模型权重结果具有很好的一致性。2种方法判定指数模型都具有最高的模型权重,而高斯模型的模型权重最低。通过2种方法计算的模型权重值略有差别,原因主要有：第一,在理论上MLBMA计算模型权重的方法是BMA的近似,近似过程中会产生相应的截断误差；第二,在此算例中由于球状模型分段函数的特点,其参数权重分布具有一定双峰的特征,这点从图2c中的BMA参数权重等值线图可以较为明显的观察到。MLBMA产生的参数点虽然可以较好地覆盖优势参数区域,但还是一定程度上漏失掉了次优势区域。因此,使用MLBMA虽然有计算量小,无需对参数的分布进行先验主观猜测等优势,但如果参数权重的分布呈现明显的双峰或者多峰的特点,则MLBMA的计算准确性比BMA会有所降低。

图3　BMA和MLBMA的模型权重对比

图4　BMA和MLBMA的渗透系数多模型预测均值对比

BMA和MLBMA的渗透系数多模型预测均值、方差对比分别见图4、5。从对比结果看,BMA和MLBMA对预测统计信息的刻画具有很好的共同性,主要变现在图4a、4b表示的渗透系数均值和图5a、5b表示的渗透系数方差分布都非常接近。另外,图4预测均值与图1中的真实渗透系数场比较可知,相较于真实渗透系数场,虽然由于平均的作用使预测结果相对较为平滑,但是2种方法都可以很好地预测渗透系数的空间分布规律。对预测不确定性,即预测方差的确定也表明在有观测数据的地方预测误差为0,随着距观测点的距离越来越大,预测的方差或不确定性也随之变大。

图5　BMA和MLBMA的渗透系数多模型预测方差对比

3结语

(1)MLBMA是BMA在极大似然估计理论基础上的近似,通过改进传统的MLBMA不能考虑参数不确定性的局限,使MLBMA对BMA近似得更为准确。

(2)MLBMA由于在确定参数空间时已经考虑了观测数据的影响,并确定了优势有效参数的产生区域,而BMA则需要对参数空间进行主观的先验猜测,并在对先验参数空间进行充分取样计算后,才能依据后验参数权重得到优势参数分布区域。因此,MLBMA在能够保证计算精度的同时,避免了大量的无效或者低效的模型预测计算,提高了计算效率。

(3)多模型地质统计模型对地层渗透系数空间分布预测的算例表明,BMA和MLBMA在后验模型权重及渗透系数预测统计规律的计算结果都具有很好的一致性。但是,如果模型参数权重的分布呈现出明显的多峰趋势,则MLBMA的计算精度有可能降低。

参考文献：

[1]NEUMANS.MaximumlikelihoodBayesianaveragingofuncertainmodelpredictions[J].StochasticEnvironmentalResearchandRiskAssessment, 2003, 17(5)： 291- 305.

[2]YEMing,NEUMANSP,MEYERPD.MaximumlikelihoodBayesianaveragingofspatialvariabilitymodelsinunsaturatedfracturedtuff[J].WaterResourcesResearch, 2004, 40(5)： 1- 17.

[3]NOWAKW,DeBARROSFPJ,RUBINY.Bayesiangeostatisticaldesign：Task-drivenoptimalsiteinvestigationwhenthegeostatisticalmodelisuncertain[J].WaterResourcesResearch, 2010, 46(3)： 3535- 3552.

[4]XUELiang,ZHANGDongxiao.AmultimodeldataassimilationframeworkviatheensembleKalmanfilter[J].WaterResourcesResearch, 2014, 50(5)： 4197- 4219.

[5]XUELiang,ZHANGDongxiao,GUADAGNINIA.MultimodelBayesiananalysisofgroundwaterdataworth[J].WaterResourcesResearch, 2014, 50(11)： 8481- 8496.

[6]王慧亮, 李叙勇, 解莹. 多模型方法在非点源污染负荷中的应用展望[J]. 水科学进展, 2011, 22(5)： 727- 732.

[7]BEVENK,BINLEYA.Thefutureofdistributedmodels:modelcalibrationanduncertaintyprediction[J].Hydrologicalprocesses, 1992, 6(3)： 279- 298.

[8]MONTANARIA.Largesamplebehaviorsofthegeneralizedlikelihooduncertaintyestimation(GLUE)inassessingtheuncertaintyofrainfall-runoffsimulations[J].WaterResourcesResearch, 2005, 41(8)： 8406- 8419.

[9]STEDINGERJR,VOGELRM,LEESU,etal.Appraisalofthegeneralizedlikelihooduncertaintyestimation(GLUE)method[J].Waterresourcesresearch, 2008, 44(12)： 6- 23.

[10]HOETINGJA,MADIGAND,RAFTERYAE,etal.Bayesianmodelaveraging：atutorial[J].Statisticalscience, 1999, 14(4)： 382- 401.

(责任编辑杨健)

收稿日期：2015- 10- 07

基金项目：国家自然科学基金资助项目(41402199,41502237)；中国石油大学(北京)引进人才科研启动基金项目(2462014YJRC038)；油气资源与探测国家重点实验室青年人才培育课题(PRP/indep- 4- 1409)

作者简介：薛亮(1983—),男,河北邢台人,副教授,博士,主要从事孔隙及裂隙介质中渗流规律、不确定性分析和数据同化等方面的研究．

中图分类号：P641

文献标识码：A

文章编号：0559- 9342(2016)04- 0031- 05

ComparativeStudyonBayesianMultimodelAnalysisMethods

XUELiang1,3,XIAQiang2

(1.DepartmentofOil-gasFieldDevelopment,CollegeofPetroleumEngineering,ChinaUniversityofPetroleum,Beijing102249,China; 2.StateKeyLaboratoryofGeohazardPreventionandGeoenvironmentProtection,ChengduUniversityofTechnology,Chengdu610059,Sichuan,China; 3.StateKeyLaboratoryofPetroleumResourcesandProspecting,ChinaUniversityofPetroleum,Beijing102249,China)

Abstract:Multimodel analysis plays more important role in decision-making and risk assessment in recent years due to its capability to take conceptual model uncertainty into account. Two commonly used multimodel analysis methods, Bayesian model averaging method and its maximum likelihood version, are introduced. The maximum likelihood Bayesian averaging method has been improved to take parameter uncertainty into account. These methods are applied to analyze the spatial distribution of log hydraulic conductivity. The results show that these two methods are consistent with each other in terms of identifying the parameter space, determining the posterior model weights and predicting the log hydraulic conductivity distribution. The maximum likelihood Bayesian averaging method can be incorporated with the well-developed inverse modeling methods in hydro-geological researches．

Key Words:uncertainty analysis; Bayesian model averaging; maximum likelihood estimation; geostatistics