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

薛　亮,夏　強(qiáng)

(1.中國石油大學(xué)石油工程學(xué)院油氣田開發(fā)工程系,北京102249；2.成都理工大學(xué)地質(zhì)災(zāi)害防治與地質(zhì)環(huán)境保護(hù)國家重點(diǎn)實(shí)驗(yàn)室,四川成都610059；3.油氣資源與探測國家重點(diǎn)實(shí)驗(yàn)室,北京102249)

?

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

薛亮1,3,夏強(qiáng)2

(1.中國石油大學(xué)石油工程學(xué)院油氣田開發(fā)工程系,北京102249；2.成都理工大學(xué)地質(zhì)災(zāi)害防治與地質(zhì)環(huán)境保護(hù)國家重點(diǎn)實(shí)驗(yàn)室,四川成都610059；3.油氣資源與探測國家重點(diǎn)實(shí)驗(yàn)室,北京102249)

摘要：多模型分析能夠考慮模型本身存在的不確定性,在決策分析和風(fēng)險(xiǎn)評估中具有越來越重要的作用。對具有嚴(yán)謹(jǐn)統(tǒng)計(jì)分析理論基礎(chǔ)的貝葉斯模型平均法和極大似然貝葉斯多模型平均法做了詳細(xì)介紹,并改進(jìn)了傳統(tǒng)極大似然貝葉斯多模型平均法不能考慮參數(shù)不確定性的不足,使極大似然貝葉斯多模型平均法對貝葉斯模型平均法近似得更為準(zhǔn)確。地質(zhì)統(tǒng)計(jì)多模型對地層滲透系數(shù)的預(yù)測分析結(jié)果表明,2種多模型分析方法在參數(shù)空間確定、模型后驗(yàn)權(quán)重和滲透系數(shù)預(yù)測方面都具有很好的一致性。極大似然貝葉斯多模型平均法能和水文學(xué)現(xiàn)有的參數(shù)估計(jì)方法很好結(jié)合,且計(jì)算量小。

關(guān)鍵詞：不確定性分析；貝葉斯模型平均法；極大似然估計(jì)；地質(zhì)統(tǒng)計(jì)

0引言

地下水資源管理及污染防治一直以來都是國內(nèi)外學(xué)者關(guān)心的熱點(diǎn)問題。合理的地下水資源管理規(guī)劃和有效的污染防治方案制定都需要建立相應(yīng)的模型,以對地下水的流動(dòng)動(dòng)態(tài)和污染物運(yùn)移的趨勢做出相關(guān)預(yù)測,并依據(jù)預(yù)測結(jié)果制定相應(yīng)的決策。然而,地下水系統(tǒng)是一個(gè)復(fù)雜的自然系統(tǒng),使用單一概念的模型對地下水系統(tǒng)進(jìn)行描述,在解決決策分析問題中存在著諸多不足。Neuman在多模型分析研究中提出,使用單一概念模型存在2種類型的錯(cuò)誤：一是由于依賴不準(zhǔn)確模型而造成的統(tǒng)計(jì)偏差；二是由于模型空間考慮不充分而造成決策風(fēng)險(xiǎn)低估[1]。多模型分析研究已經(jīng)成為國內(nèi)外水文領(lǐng)域研究的前沿?zé)狳c(diǎn)問題[2- 6]。

多模型分析方法主要有3種：廣義似然不確定性估計(jì)法(GeneralizedLikelihoodUncertaintyEstimation,GLUE)、貝葉斯模型平均法(BayesianModelAveraging,BMA)和極大似然貝葉斯模型平均法(MaximumLikelihoodBayesianModelAveraging,MLBMA)。GLUE最早在水文學(xué)中提出,現(xiàn)已成為一種常用的多模型分析方法[7]。然而,GLUE通常被認(rèn)為缺乏嚴(yán)格的統(tǒng)計(jì)理論依據(jù),似然方程和模型保留準(zhǔn)則需主觀確定,參數(shù)與預(yù)測的統(tǒng)計(jì)性質(zhì)不統(tǒng)一,沒有在預(yù)測概率分布中明確考慮模型誤差等[8-9]。BMA由統(tǒng)計(jì)學(xué)家在貝葉斯分析框架下提出,理論方法嚴(yán)謹(jǐn),能考慮模型、參數(shù)、觀測等諸多因素導(dǎo)致的不確定性,應(yīng)用于諸多領(lǐng)域。Hoeting等人對BMA做了詳細(xì)的研究和探討[10]。為使BMA更易于與水文地質(zhì)領(lǐng)域完善的極大似然模型反演方法相結(jié)合,Neuman建立了MLBMA[1],并證明了多模型方法相較于單一模型在解決統(tǒng)計(jì)決策分析問題上的優(yōu)勢[2],但該方法忽略了參數(shù)的不確定性。

本文對傳統(tǒng)的MLBMA進(jìn)行了改進(jìn),考慮參數(shù)的不確定性,使MLBMA在考慮不確定性因素方面獲得和BMA相對等的功能。對BMA和MLBMA理論方法分別進(jìn)行介紹,對比2種方法的實(shí)現(xiàn)特點(diǎn)和執(zhí)行流程,并通過地質(zhì)統(tǒng)計(jì)多模型分析，對2種多模型方法的分析結(jié)果進(jìn)行比較。

1貝葉斯多模型分析方法

1.1貝葉斯模型平均法(BMA)

預(yù)測量的多模型后驗(yàn)概率分布是組成模型集的單個(gè)模型給出的預(yù)測量后驗(yàn)概率分布的加權(quán)平均，進(jìn)行加權(quán)平均時(shí)的權(quán)重為該給定單個(gè)模型的后驗(yàn)概率。假設(shè)預(yù)測量為Δ,K個(gè)模型都可以對預(yù)測量進(jìn)行估計(jì)預(yù)測,這些模型組成模型集合M={M1,…,MK},模型集中的每個(gè)模型Mk都由參數(shù)向量θk進(jìn)行表征,觀測數(shù)據(jù)為D,則條件于該觀測數(shù)據(jù)的預(yù)測量多模型后驗(yàn)概率分布可由下式定義

(1)

考慮參數(shù)存在不確定性時(shí),給定模型的預(yù)測量后驗(yàn)概率為

(2)

模型后驗(yàn)概率的確定是多模型分析的關(guān)鍵點(diǎn),可以通過貝葉斯定理進(jìn)行計(jì)算，公式為

(3)

邊際似然函數(shù)可以通過下式進(jìn)行計(jì)算

(4)

得到模型的后驗(yàn)概率后,預(yù)測量的多模型均值為

(5)

(6)

(7)

(8)

BMA多模型分析方法的執(zhí)行步驟為：①定義每個(gè)模型中參數(shù)的先驗(yàn)概率分布,并隨機(jī)取樣產(chǎn)生1組參數(shù)實(shí)現(xiàn)；②對所有模型中每組參數(shù)的每個(gè)實(shí)現(xiàn)對應(yīng)的似然函數(shù)進(jìn)行計(jì)算；③通過式(4)計(jì)算邊際似然函數(shù),并計(jì)算考慮參數(shù)不確定性的預(yù)測均值和方差；④通過式(5)、(6)計(jì)算考慮模型不確定性的均值和方差。

1.2極大似然貝葉斯模型平均法(MLBMA)

對于每個(gè)給定的模型而言,在定義模型參數(shù)的概率密度函數(shù)并據(jù)此生成模型參數(shù)實(shí)現(xiàn)的過程中,并沒有考慮觀測數(shù)據(jù)對參數(shù)分布的影響,參數(shù)的分布是依據(jù)經(jīng)驗(yàn)主觀定義的,這會導(dǎo)致由于有效計(jì)算量的減少造成總體計(jì)算量的增加。如果參數(shù)分布定義不準(zhǔn)確,還可能導(dǎo)致最終的統(tǒng)計(jì)分析結(jié)果存在偏差。此外,式(3)、(4)中積分的準(zhǔn)確計(jì)算,都需要通過模特卡羅方法對每個(gè)模型的參數(shù)進(jìn)行大量取樣,這也會造成計(jì)算量的增加。

(9)

式(2)可由按照此概率分布產(chǎn)生的參數(shù)實(shí)現(xiàn),通過模特卡羅積分進(jìn)行計(jì)算。

在極大似然估計(jì)理論框架下,模型的后驗(yàn)概率或者模型權(quán)重可以通過下式進(jìn)行估計(jì)

(10)

KIC=-2ln[

(11)

MLBMA多模型分析方法的執(zhí)行步驟為：①通過極大似然估計(jì)方法得到每個(gè)模型的參數(shù)估計(jì)值及其估計(jì)協(xié)方差,并得到相應(yīng)的似然函數(shù)值；②依據(jù)參數(shù)的多元正態(tài)分布式(9)隨機(jī)取樣產(chǎn)生1組參數(shù)實(shí)現(xiàn),并計(jì)算考慮參數(shù)不確定性的預(yù)測均值和方差；③通過式(10)計(jì)算后驗(yàn)?zāi)Ｐ透怕剩虎芡ㄟ^式(5)、(6)計(jì)算考慮模型不確定性的均值和方差。

2算例研究

為對比BMA和MLBMA這2種多模型方法,本文使用地質(zhì)統(tǒng)計(jì)模型對地層滲透系數(shù)的空間非均質(zhì)性進(jìn)行分析。二維50×50單位長度的參考或真實(shí)滲透系數(shù)場由高斯型截?cái)鄡绾瘮?shù)地質(zhì)統(tǒng)計(jì)模型(TruncatedPowerVariogramwithGaussianmode,TPVG)產(chǎn)生,可表示為

(12)

生成的滲透系數(shù)參考場見圖1。生成參考場模型的參數(shù)值為θ=(A,H,λu)T=(0.1,0.25,5)T。100個(gè)隨機(jī)選擇數(shù)據(jù)點(diǎn)上的滲透系數(shù)值被用作觀測數(shù)據(jù)進(jìn)行多模型分析,其空間位置見圖1中的“D”字符所示。

圖1　生成的滲透系數(shù)參考場

(13)

球狀模型：γ(s)=

式中,σ2為滲透系數(shù)場的方差；R為滲透系數(shù)場的變程。

此次研究選取了3個(gè)參數(shù)的截?cái)鄡绾瘮?shù)地質(zhì)統(tǒng)計(jì)模型作為生成參考場的真實(shí)模型,比包含2個(gè)參數(shù)的指數(shù)、高斯和球狀模型更為復(fù)雜,并能更靈活地刻畫滲透系數(shù)場的空間變異性。另外,為使構(gòu)建的算例能反映真實(shí)的模擬需要,生成參考場的截?cái)鄡绾瘮?shù)模型并沒有放到進(jìn)行多模型分析的模型集合內(nèi)。因?yàn)橥ǔＴ趯?shí)際情況下,真實(shí)模型一般都比選擇的模型復(fù)雜,且真實(shí)模型都是未知的。

在執(zhí)行BMA分析時(shí),假設(shè)這3種選擇模型的參數(shù)服從多元均勻分布,方差的取值范圍為0～3,變程為0～30,通過拉丁超立方方法產(chǎn)生10 000組參數(shù)集合,并對每組參數(shù)依據(jù)觀測數(shù)據(jù)計(jì)算相應(yīng)的似然函數(shù)值和參數(shù)權(quán)重。在執(zhí)行MLBMA時(shí),3種選擇模型都通過100個(gè)觀測數(shù)據(jù)進(jìn)行參數(shù)校正,得到參數(shù)估計(jì)值和估計(jì)協(xié)方差,假設(shè)校正后的參數(shù)服從多元正態(tài)分布,并通過式(9)產(chǎn)生2 000組參數(shù)實(shí)現(xiàn)。2種方法的數(shù)值計(jì)算結(jié)果都是通過MATLAB自行編制代碼,結(jié)合模型參數(shù)校正軟件PEST計(jì)算獲得。3種模型在BMA和MLBMA中的參數(shù)對比結(jié)果見圖2。圖2中,黑點(diǎn)為在MLBMA中產(chǎn)生的2 000組方差和變程組成的參數(shù),背景等值線為在BMA中10 000組方差和變程參數(shù)對應(yīng)的插值參數(shù)權(quán)重。

從圖2可以看出,MLBMA產(chǎn)生的參數(shù)集與BMA具有很強(qiáng)的一致性,主要表現(xiàn)在MLBMA產(chǎn)生的參數(shù)區(qū)域能夠很好地覆蓋BMA具有較高權(quán)重的參數(shù)區(qū)域。BMA首先通過先驗(yàn)的猜測選定參數(shù)取值區(qū)域,然后通過觀測數(shù)據(jù)計(jì)算似然函數(shù)值,并得到該參數(shù)的權(quán)重。只有當(dāng)所有的生成參數(shù)都計(jì)算完畢后,才能得到在加權(quán)平均中的起重要作用的有效參數(shù)。在此例中,有效參數(shù)主要存在于等值線圖的黑色區(qū)域,而灰色區(qū)域的參數(shù)由于權(quán)重很低,在得到滲透系數(shù)預(yù)測統(tǒng)計(jì)結(jié)果的加權(quán)平均過程中起不到重要作用。因此,BMA要求的參數(shù)取樣量較大(此處為10 000組),這也造成了計(jì)算量相對較高。MLBMA通過參數(shù)校正過程,利用觀測數(shù)據(jù)首先確定參數(shù)產(chǎn)生的區(qū)域。由于多元正態(tài)分布的假設(shè),其在極大似然估計(jì)參數(shù)值附近區(qū)域會產(chǎn)生較多的取樣點(diǎn),隨著向區(qū)域外部的擴(kuò)展,取樣點(diǎn)數(shù)隨之較少,每個(gè)參數(shù)取樣點(diǎn)處的滲透系數(shù)預(yù)測計(jì)算結(jié)果都是等權(quán)重平均。因此,此處2 000組參數(shù)在計(jì)算最終滲透系數(shù)預(yù)測統(tǒng)計(jì)結(jié)果時(shí)起的作用都是一樣的,避免了無效的預(yù)測計(jì)算,節(jié)省了計(jì)算量。

圖2　3種模型在BMA和MLBMA中的參數(shù)對比

后驗(yàn)?zāi)Ｐ透怕驶蛘吣Ｐ蜋?quán)重是多模型分析的重要基礎(chǔ)。BMA和MLBMA的模型權(quán)重對比見圖3。從圖3可以看出,2個(gè)多模型分析方法計(jì)算的模型權(quán)重結(jié)果具有很好的一致性。2種方法判定指數(shù)模型都具有最高的模型權(quán)重,而高斯模型的模型權(quán)重最低。通過2種方法計(jì)算的模型權(quán)重值略有差別,原因主要有：第一,在理論上MLBMA計(jì)算模型權(quán)重的方法是BMA的近似,近似過程中會產(chǎn)生相應(yīng)的截?cái)嗾`差；第二,在此算例中由于球狀模型分段函數(shù)的特點(diǎn),其參數(shù)權(quán)重分布具有一定雙峰的特征,這點(diǎn)從圖2c中的BMA參數(shù)權(quán)重等值線圖可以較為明顯的觀察到。MLBMA產(chǎn)生的參數(shù)點(diǎn)雖然可以較好地覆蓋優(yōu)勢參數(shù)區(qū)域,但還是一定程度上漏失掉了次優(yōu)勢區(qū)域。因此,使用MLBMA雖然有計(jì)算量小,無需對參數(shù)的分布進(jìn)行先驗(yàn)主觀猜測等優(yōu)勢,但如果參數(shù)權(quán)重的分布呈現(xiàn)明顯的雙峰或者多峰的特點(diǎn),則MLBMA的計(jì)算準(zhǔn)確性比BMA會有所降低。

圖3　BMA和MLBMA的模型權(quán)重對比

圖4　BMA和MLBMA的滲透系數(shù)多模型預(yù)測均值對比

BMA和MLBMA的滲透系數(shù)多模型預(yù)測均值、方差對比分別見圖4、5。從對比結(jié)果看,BMA和MLBMA對預(yù)測統(tǒng)計(jì)信息的刻畫具有很好的共同性,主要變現(xiàn)在圖4a、4b表示的滲透系數(shù)均值和圖5a、5b表示的滲透系數(shù)方差分布都非常接近。另外,圖4預(yù)測均值與圖1中的真實(shí)滲透系數(shù)場比較可知,相較于真實(shí)滲透系數(shù)場,雖然由于平均的作用使預(yù)測結(jié)果相對較為平滑,但是2種方法都可以很好地預(yù)測滲透系數(shù)的空間分布規(guī)律。對預(yù)測不確定性,即預(yù)測方差的確定也表明在有觀測數(shù)據(jù)的地方預(yù)測誤差為0,隨著距觀測點(diǎn)的距離越來越大,預(yù)測的方差或不確定性也隨之變大。

圖5　BMA和MLBMA的滲透系數(shù)多模型預(yù)測方差對比

3結(jié)語

(1)MLBMA是BMA在極大似然估計(jì)理論基礎(chǔ)上的近似,通過改進(jìn)傳統(tǒng)的MLBMA不能考慮參數(shù)不確定性的局限,使MLBMA對BMA近似得更為準(zhǔn)確。

(2)MLBMA由于在確定參數(shù)空間時(shí)已經(jīng)考慮了觀測數(shù)據(jù)的影響,并確定了優(yōu)勢有效參數(shù)的產(chǎn)生區(qū)域,而BMA則需要對參數(shù)空間進(jìn)行主觀的先驗(yàn)猜測,并在對先驗(yàn)參數(shù)空間進(jìn)行充分取樣計(jì)算后,才能依據(jù)后驗(yàn)參數(shù)權(quán)重得到優(yōu)勢參數(shù)分布區(qū)域。因此,MLBMA在能夠保證計(jì)算精度的同時(shí),避免了大量的無效或者低效的模型預(yù)測計(jì)算,提高了計(jì)算效率。

(3)多模型地質(zhì)統(tǒng)計(jì)模型對地層滲透系數(shù)空間分布預(yù)測的算例表明,BMA和MLBMA在后驗(yàn)?zāi)Ｐ蜋?quán)重及滲透系數(shù)預(yù)測統(tǒng)計(jì)規(guī)律的計(jì)算結(jié)果都具有很好的一致性。但是,如果模型參數(shù)權(quán)重的分布呈現(xiàn)出明顯的多峰趨勢,則MLBMA的計(jì)算精度有可能降低。

參考文獻(xiàn)：

[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]王慧亮, 李敘勇, 解瑩. 多模型方法在非點(diǎn)源污染負(fù)荷中的應(yīng)用展望[J]. 水科學(xué)進(jìn)展, 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.

(責(zé)任編輯楊健)

收稿日期：2015- 10- 07

基金項(xiàng)目：國家自然科學(xué)基金資助項(xiàng)目(41402199,41502237)；中國石油大學(xué)(北京)引進(jìn)人才科研啟動(dòng)基金項(xiàng)目(2462014YJRC038)；油氣資源與探測國家重點(diǎn)實(shí)驗(yàn)室青年人才培育課題(PRP/indep- 4- 1409)

作者簡介：薛亮(1983—),男,河北邢臺人,副教授,博士,主要從事孔隙及裂隙介質(zhì)中滲流規(guī)律、不確定性分析和數(shù)據(jù)同化等方面的研究．

中圖分類號：P641

文獻(xiàn)標(biāo)識碼：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