地下水模擬不確定性問題的多模型分析

宋 凱, 劉 丹, 劉 建

(西南交通大學地球科學與環境工程學院, 四川 成都 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.