貝葉斯多模型分析方法的對比研究

薛　亮,夏　強

(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