基于自適應(yīng)多保真度Co-Kriging 代理模型的地下水污染源反演識(shí)別

安永凱,張巖祥,閆雪嫚(1.長(zhǎng)安大學(xué)水利與環(huán)境學(xué)院,陜西 西安 710054；2.長(zhǎng)安大學(xué)旱區(qū)地下水文與生態(tài)效應(yīng)教育部重點(diǎn)實(shí)驗(yàn)室,陜西 西安 710054；.中國(guó)電建集團(tuán)西北勘測(cè)設(shè)計(jì)研究院有限公司,陜西 西安 710065；4.西北大學(xué)城市與環(huán)境學(xué)院,陜西 西安 710127)

地下水污染具有存在的隱蔽性和發(fā)現(xiàn)的滯后性等特點(diǎn),致使人們難以直接了解和掌握地下水污染的相關(guān)狀況(如污染源的空間位置和釋放歷史),這對(duì)地下水污染修復(fù)方案合理制定、風(fēng)險(xiǎn)評(píng)估以及責(zé)任認(rèn)定帶來(lái)極大的困難[1-2].地下水污染源反演識(shí)別能夠根據(jù)污染場(chǎng)地地下水水位和污染質(zhì)濃度監(jiān)測(cè)數(shù)據(jù),結(jié)合現(xiàn)場(chǎng)調(diào)查、動(dòng)態(tài)試驗(yàn)、專(zhuān)業(yè)知識(shí)和專(zhuān)家經(jīng)驗(yàn)等輔助信息,對(duì)描述地下水污染的數(shù)值模擬模型進(jìn)行反演求解,從而識(shí)別確定含水層中污染源個(gè)數(shù)、位置和釋放歷史等地下水污染源參數(shù)[3-4].

現(xiàn)有地下水污染源反演識(shí)別方法主要可分為以下幾類(lèi):解析法、直接法、模擬優(yōu)化法以及隨機(jī)統(tǒng)計(jì)法[5-6].其中,模擬優(yōu)化法是地下水污染源反演識(shí)別研究中應(yīng)用最廣泛的方法[7-8].但模擬優(yōu)化法只能給出污染源參數(shù)的“最優(yōu)解”,無(wú)法充分描述其不確定性[9].與模擬優(yōu)化法不同,應(yīng)用隨機(jī)統(tǒng)計(jì)法開(kāi)展地下水污染源反演識(shí)別研究可獲得未知污染源參數(shù)的后驗(yàn)分布,不僅可充分描述其不確定性,且所得到的信息量更大從而能夠?yàn)榈叵滤廴拘迯?fù)治理、風(fēng)險(xiǎn)評(píng)估等提供更多的參考信息[10-11].然而,應(yīng)用隨機(jī)統(tǒng)計(jì)法進(jìn)行地下水污染源反演識(shí)別需成千上萬(wàn)次地運(yùn)行地下水溶質(zhì)運(yùn)移數(shù)值模擬模型,由此產(chǎn)生了龐大的計(jì)算負(fù)荷及冗長(zhǎng)的計(jì)算時(shí)間[10].建立模擬模型的代理模型是解決該問(wèn)題的有效途徑[12-13].

代理模型可以用非常低的計(jì)算代價(jià)對(duì)復(fù)雜模擬模型的輸入輸出關(guān)系進(jìn)行近似代理[14],根據(jù)其構(gòu)建方式,代理模型可分為兩大類(lèi):第一類(lèi)是響應(yīng)面代理模型,第二類(lèi)是多保真度(MF)代理模型[15].前者是僅利用高保真度(HF)模型(即原始模擬模型)輸入輸出樣本構(gòu)造的數(shù)據(jù)驅(qū)動(dòng)代理模型,具有易于實(shí)施的優(yōu)點(diǎn),在水文地質(zhì)領(lǐng)域得到了廣泛應(yīng)用[2,16].HF 模型精度高但計(jì)算成本也高,導(dǎo)致第一類(lèi)代理模型的構(gòu)建仍需昂貴的計(jì)算成本,因此學(xué)者們提出綜合運(yùn)用HF模型信息和低保真度(LF)模型信息來(lái)構(gòu)建MF代理模型.其中,LF 模型可通過(guò)離散網(wǎng)格粗化、模型降階等簡(jiǎn)化方式得到,相比于HF 模型,LF 模型精度較差但計(jì)算成本也低,可以提供充足的樣本[17].MF 代理模型旨在通過(guò)MF 模擬框架充分挖掘HF 模型和LF 模型之間的相關(guān)關(guān)系,進(jìn)而用少量HF 模型樣本和大量LF 模型樣本構(gòu)建兼顧HF 模型精度與LF 模型效率的代理模型[18].相比于常用的響應(yīng)面代理模型,MF 代理模型能夠進(jìn)一步提高計(jì)算效率且對(duì)模擬模型的逼近精度更高[19-20].但在地下水污染源反演識(shí)別研究中,MF 代理模型卻鮮有報(bào)道.同時(shí),采用自適應(yīng)更新多保真度代理模型策略,即將每一輪地下水污染源反演識(shí)別結(jié)果帶入HF 模型,獲得新的訓(xùn)練樣本,進(jìn)而重新構(gòu)建MF 代理模型,可以有效進(jìn)一步提升MF 代理模型對(duì)模擬模型的逼近精度,以此重新進(jìn)行地下水污染源反演識(shí)別,能夠極大地提高地下水污染源反演識(shí)別精度.

基于此,本文研究應(yīng)用集成差分進(jìn)化算法的Co-Kriging 方法建立地下水?dāng)?shù)值模擬模型的MF 代理模型,在此基礎(chǔ)上,耦合多保真度Co-Kriging 代理模型和MCMC-DREAM(D)算法,并采用自適應(yīng)更新多保真度代理模型策略進(jìn)行地下水污染源反演識(shí)別,旨在為地下水污染修復(fù)治理提供參考.

1 研究方法

技術(shù)路線如圖1 所示.

圖1 技術(shù)路線Fig.1 Technology roadmap

1.1 Co-Kriging 代理模型

Co-Kriging 代理模型是Kriging 代理模型的一種擴(kuò)展性形式,由兩組互不干擾的高、低保真度樣本共同構(gòu)建而成的MF 代理模型[21].相比于只用高保真度樣本建立的Kriging 代理模型,Co-Kriging 代理模型在保持非常高的代理精度的同時(shí),能夠有效降低獲取高保真度樣本所需的時(shí)間和計(jì)算代價(jià)[22].

假設(shè)有兩組高、低保真度樣本,分別為{Xh,Yh}和{Xl,Yl}, 其中Xh和Xl為樣本輸入,且Xh?Xl, Yh和 Yl為樣本輸出.這些樣本集合起來(lái)表示如下:

式中:nh和nl分別為高、低保真度樣本數(shù)量.

使用Zh和Zl分別表示高、低保真度樣本的高斯過(guò)程,兩者之間的關(guān)系可以使用式(3)表示:

式中:ρ為縮放因子,Zd表示Zh與ρZl之間差異的高斯過(guò)程.

參照 Kriging 代理模型中協(xié)方差矩陣的結(jié)構(gòu),Co-Kriging 代理模型的完整協(xié)方差結(jié)構(gòu)如式(4)所示[23]:

式中:ψl和ψd為相關(guān)矩陣,和可以參照Kriging代理模型求解.

協(xié)方差矩陣C中有5個(gè)超參數(shù)需要優(yōu)化,分別為θl、θd、pl、pd和ρ.可以使用極大似然估計(jì)(Maximum Likelihood Estimate, MLE)獲取上述超參數(shù).

通過(guò)使式(5)取得最大值求得θl和pl:

通過(guò)使式(6)取得最大值求得θd、pd和ρ:

在獲得超參數(shù)之后,可以使用式(7)預(yù)測(cè)任意一點(diǎn)在高保真度的高斯過(guò)程中的響應(yīng)值.

1.2 差分進(jìn)化算法

上述Co-Kriging 代理模型的超參數(shù)往往需要配合優(yōu)化算法方可求得,優(yōu)化算法—差分進(jìn)化算法(DE).DE 是一種高效的全局優(yōu)化算法,與遺傳算法的進(jìn)化流程非常相似,都包括變異、交叉和選擇操作,但這些操作具體定義與遺傳算法有所不同[24].如圖2 所示,DE 的進(jìn)化流程如下:(1)確定適應(yīng)度函數(shù)、種群大小、縮放因子、交叉概率和進(jìn)化代數(shù).(2)隨機(jī)產(chǎn)生初始種群.(3)計(jì)算初始種群的適應(yīng)度值.(4)對(duì)初始種群進(jìn)行變異和交叉,得到中間種群.(5)計(jì)算中間種群的適應(yīng)度值.(6)在初始種群和中間種群中選擇個(gè)體,更新初始種群.(7)判斷是否達(dá)到終止條件或最大進(jìn)化代數(shù),若是,則終止,若否,轉(zhuǎn)步驟(4).

圖2 差分進(jìn)化算法的進(jìn)化流程Fig.2 Process of differential evolution algorithm

1.3 貝葉斯推理

貝葉斯推理源于貝葉斯理論.貝葉斯公式是貝葉斯推理的數(shù)學(xué)表達(dá),用于表達(dá)先驗(yàn)分布和后驗(yàn)分布的關(guān)系,目前貝葉斯推理被廣泛應(yīng)用于模型參數(shù)識(shí)別及不確定性分析[25].貝葉斯公式根據(jù)先驗(yàn)知識(shí)設(shè)置參數(shù)的初始分布,通過(guò)分析樣本的概率后,得到參數(shù)的后驗(yàn)分布[26],如式(8)所示:

在實(shí)際應(yīng)用中,除了極個(gè)別特別簡(jiǎn)單的模型可以通過(guò)解析解推斷參數(shù)的后驗(yàn)分布,對(duì)于絕大多數(shù)的復(fù)雜模型,只能通過(guò)采樣的方式獲取參數(shù)的后驗(yàn)分布[27-28].馬爾科夫鏈蒙特卡洛(MCMC)能夠從后驗(yàn)分布采樣并進(jìn)行統(tǒng)計(jì)分析,從而得到后驗(yàn)分布的統(tǒng)計(jì)特征[29-30].目前,MCMC 方法已經(jīng)開(kāi)發(fā)出了眾多采樣算法.其中,DREAM(D)是在DREAM算法的基礎(chǔ)上改進(jìn)而來(lái)的,與原始DREAM 算法不同的是,它在處理連續(xù)型變量的同時(shí)能夠有效處理離散型變量,并且能夠保持細(xì)致的平衡性和遍歷性[31].

然后進(jìn)行馬爾科夫鏈進(jìn)化:

(2)對(duì)于第i(i=1,2, …,N)條馬爾科夫鏈,產(chǎn)生候選種群的公式(9)如下所示:

(3)定義CR(0 ≤CR ≤ 1)為進(jìn)行子空間演化的交叉概率,進(jìn)而決定候選種群是否取代初始種群:

式中:U 為隨機(jī)數(shù),從均勻分布 U( 0,1)采樣得到.

(4)計(jì)算似然度 p( zi),并進(jìn)而計(jì)算接受率α(θi, zi):

根據(jù)接受率 α(θi,zi)的計(jì)算結(jié)果,判斷是否接受zi.

(5)應(yīng)用Gelman-Rubin 收斂診斷方法[32]判斷收斂性,計(jì)算收斂性判斷指標(biāo)R:

式中:g 為DREAM(D)算法中的馬爾科夫鏈長(zhǎng)度;q 為MC 數(shù)目;B 表示q 條馬爾科夫鏈平均值的方差;W 為q 條馬爾科夫鏈內(nèi)方差的平均值.若R≤12,則認(rèn)為采樣過(guò)程已收斂,結(jié)束計(jì)算,統(tǒng)計(jì)分析待求變量θ的后驗(yàn)分布;否則重復(fù)步驟.

2 數(shù)值試驗(yàn)

2.1 數(shù)值算例概述

本文針對(duì)2 個(gè)數(shù)值算例開(kāi)展研究.算例1:考慮一個(gè)二維均質(zhì)各向同性承壓含水層,該含水層?xùn)|西長(zhǎng)為300L,南北寬為240L,東、西兩側(cè)邊界(Γ1和Γ3)均為線性變化的定水頭邊界,南、北邊界(Γ2和Γ4)均為零通量邊界,如圖3(a)所示.算例2:考慮一個(gè)二維非均質(zhì)各向同性承壓水層,該含水層劃分為3 個(gè)參數(shù)分區(qū),僅對(duì)滲透系數(shù)進(jìn)行分區(qū),東西長(zhǎng)為300L,南北寬為200L,西北、東南邊界(Γ1和Γ3)均為定水頭邊界,其余(Γ2和Γ4)均為零通量邊界,如圖3(b)所示.2 個(gè)算例的模擬期均為900T,模擬區(qū)內(nèi)有1 個(gè)地下水污染源,3 個(gè)地下水污染觀測(cè)井.

圖3 算例含水層平面示意Fig.3 Aquifer plan schematic diagram

本研究假設(shè)污染物不會(huì)發(fā)生化學(xué)反應(yīng)和生物降解作用.算例1:假設(shè)污染源初始釋放時(shí)間(τ)和釋放強(qiáng)度(q)未知,并作為待識(shí)別變量,滲透系數(shù)(K)、縱向彌散度(D)、橫向彌散度與縱向彌散度比值(α)和孔隙度(n)等水文地質(zhì)參數(shù)均質(zhì)且已知,待識(shí)別變量先驗(yàn)分布和水文地質(zhì)參數(shù)取值如表1 所示.算例2:假設(shè)污染源初始釋放時(shí)間(τ)、釋放強(qiáng)度(q)、3 個(gè)分區(qū)滲透系數(shù)(K1,K2,K3)和縱向彌散度(D)未知,并作為待識(shí)別變量,橫向彌散度與縱向彌散度比值(α)和孔隙度(n)等水文地質(zhì)參數(shù)已知,待識(shí)別變量先驗(yàn)分布和水文地質(zhì)參數(shù)取值如表2 所示.

表1 算例1 待識(shí)別變量先驗(yàn)分布和水文地質(zhì)參數(shù)值Table 1 Prior distribution of unknown variables and hydrogeological parameter values for case 1

表2 算例2 待識(shí)別變量先驗(yàn)分布和水文地質(zhì)參數(shù)值Table 2 Prior distribution of unknown variables and hydrogeological parameter values for case 2

根據(jù)上述設(shè)定的水文地質(zhì)條件,模擬區(qū)地下水流運(yùn)動(dòng)數(shù)學(xué)模型如式(13)所示:

式中:H 為地下水水頭[L];K 為滲透系數(shù)[L/T];M 為含水層厚度[L];S 為貯水系數(shù)[-];H0和H1為已知函數(shù)[L];Γ1、Γ3為第一類(lèi)邊界條件; Γ2、Γ4為第二類(lèi)邊界條件.為邊界上某點(diǎn)外法線方向上的單位向量.

模擬區(qū)地下水溶質(zhì)運(yùn)移數(shù)學(xué)模型如式(14)所示:

式中:c 為污染質(zhì)濃度[ M / L3];Dxx和Dyy分別為x、y方向上的水動(dòng)力彌散系數(shù)[ L2/ T ];ux和uy分別為x、y 方向上的實(shí)際平均水流速度[ L / T ];f 為源匯項(xiàng)[ M /( L3T )];c0為已知濃度函數(shù)[ M / L3];c1為已知對(duì)流-彌散通量函數(shù)[ M /(L2T )];Γ1、Γ3為第三類(lèi)邊界條件; Γ2、Γ4為第二類(lèi)邊界條件.

分別對(duì)模擬區(qū)進(jìn)行精細(xì)離散和粗化離散,建立HF 模型和LF 模型,然后應(yīng)用MODFLOW 和MT3D分別求解HF 模型和LF 模型.其中,算例1 的HF 模型將模擬區(qū)剖分為 800000 個(gè)網(wǎng)格(長(zhǎng) 1000×寬800),LF 模型將模擬區(qū)剖分為2880 個(gè)網(wǎng)格(長(zhǎng)60×寬48);算例2 的HF 模型將模擬區(qū)剖分為375000 個(gè)網(wǎng)格(長(zhǎng)750×寬500),LF 模型將模擬區(qū)剖分為2400 個(gè)網(wǎng)格(長(zhǎng)60×寬40).

2.2 Co-Kriging 代理模型的建立

本文算例1 待反演識(shí)別變量包括污染源初始釋放時(shí)間(τ)和釋放強(qiáng)度(q),相應(yīng)的代理模型輸入則為上述的2 個(gè)待識(shí)別變量,輸出為3 口監(jiān)測(cè)井在3 個(gè)時(shí)刻(t=780T ,840T ,900T )監(jiān)測(cè)到的污染物濃度值,共9 個(gè)變量;算例2 待反演識(shí)別變量包括污染源初始釋放時(shí)間(τ)、釋放強(qiáng)度(q)、3 個(gè)分區(qū)滲透系數(shù)(K 1, K 2, K 3)和縱向彌散度(D),相應(yīng)的代理模型輸入則為上述的6 個(gè)待求變量,輸出為3 口監(jiān)測(cè)井在3個(gè)時(shí)刻(t=780T ,840T ,900T )監(jiān)測(cè)到的污染物濃度值,共9 個(gè)變量.

Co-Kriging 代理模型的構(gòu)建步驟如下:首先,采用拉丁超立方抽樣方法在待識(shí)別變量的先驗(yàn)分布內(nèi)分別抽取5 個(gè)和200 個(gè)訓(xùn)練樣本,將抽取的5 個(gè)訓(xùn)練樣本代入精細(xì)離散的HF 模型,同時(shí)將抽取的200個(gè)訓(xùn)練樣本代入粗化離散的LF 模型;然后,運(yùn)轉(zhuǎn)HF模型和LF 模型,分別得到對(duì)應(yīng)的模型輸出;接下來(lái),基于5 組HF 模型輸入-輸出樣本構(gòu)建Kriging 代理模型,基于5 組HF 模型和200 組LF 模型輸入-輸出樣本構(gòu)建Co-Kriging 代理模型,即MF 代理模型;最后,重新采用拉丁超立方抽樣方法在待識(shí)別變量的先驗(yàn)分布內(nèi)抽取 5 個(gè)檢驗(yàn)樣本,將其分別代入Kriging 模型和MF 代理模型,運(yùn)轉(zhuǎn)模型,得到相應(yīng)的模型輸出,并將它們與HF 模型的輸出結(jié)果進(jìn)行對(duì)比分析,從而評(píng)估Kriging 代理模型和MF 代理模型的精度.整體的MF 代理模型構(gòu)建流程如圖4 所示.

圖4 Co-Kriging 代理模型的構(gòu)建流程[18]Fig.4 Construction process of Co-Kriging surrogate model[18]

由圖5 可以看出,MF 代理模型平均相對(duì)誤差均小于Kriging 代理模型,這說(shuō)明對(duì)于相同的輸入,MF 代理模型對(duì)HF 模型的逼近精度更高,MF 代理模型比Kriging 代理模型更適合后續(xù)地下水污染源反演識(shí)別.

圖5 Kriging 代理模型與MF 代理模型精度Fig.5 Accuracy of Kriging surrogate model and MF surrogate model

對(duì)于算例1,HF 模型和MF 模型單次運(yùn)行時(shí)間分別為1964.747s 和0.016s;對(duì)于算例2,HF 模型和MF 模型單次運(yùn)行時(shí)間分別為2013.366s 和0.018s.HF 模型單次運(yùn)行時(shí)間較長(zhǎng),在地下水污染源反演識(shí)別過(guò)程中成千上萬(wàn)次調(diào)用HF 模型,需要數(shù)月才能識(shí)別完畢,因此,不合適用于地下水污染源反演識(shí)別;MF 代理模型單次運(yùn)行時(shí)間十分短,在地下水污染源反演識(shí)別過(guò)程中需要成千上萬(wàn)次調(diào)用MF 模型,僅需數(shù)分鐘便能識(shí)別完畢,且對(duì)HF 模型的逼近精度很高,十分適合用于地下水污染源反演識(shí)別.使用MF 代理模型進(jìn)行地下水污染源反演還需要考慮構(gòu)建MF 代理模型的時(shí)間,構(gòu)建MF 代理模型時(shí)間花費(fèi)主要來(lái)自運(yùn)行HF 模型,在本研究中需要運(yùn)行5 次HF 模型,需要數(shù)小時(shí)便能完成.總而言之,在地下水污染源反演識(shí)別過(guò)程中直接調(diào)用MF 代理模型所需時(shí)間遠(yuǎn)遠(yuǎn)小于調(diào)用HF 模型,更適合用于地下水污染源反演識(shí)別.

2.3 污染源反演識(shí)別

在構(gòu)建Co-Kriging 代理模型的基礎(chǔ)上,研究應(yīng)用DREAM(D)算法分別對(duì)兩個(gè)算例的待識(shí)別變量進(jìn)行反演識(shí)別.為保證待識(shí)別變量均收斂于穩(wěn)定的后驗(yàn)分布采樣過(guò)程的收斂性,本文采用Gelman-Rubin收斂診斷方法對(duì)DREAM(D)算法采樣過(guò)程的收斂性進(jìn)行檢驗(yàn).對(duì)于算例1,DREAM(D)算法的采樣過(guò)程如圖6(a)和圖6(b)所示.污染源初始釋放時(shí)間(τ)和釋放強(qiáng)度(q)這2個(gè)待識(shí)別變量均能收斂于穩(wěn)定的后驗(yàn)分布.對(duì)于算例2,DREAM(D)算法的采樣過(guò)程如圖7(a)、(b)、(e)、(f)、(i)、(j)所示.污染源初始釋放時(shí)間(τ)、釋放強(qiáng)度(q)、3 個(gè)分區(qū)滲透系數(shù)(K 1, K 2, K 3)、和縱向彌散度(D)這6 個(gè)待識(shí)別變量均能收斂于穩(wěn)定的后驗(yàn)分布.

圖6 算例1 待識(shí)別變量的采樣迭代過(guò)程及其反演識(shí)別結(jié)果Fig.6 Sampling iteration process and identification results of unknown variables for case 1

在確保采樣過(guò)程收斂性后,算例1 選取后1000組樣本進(jìn)行統(tǒng)計(jì)分析,獲得2 個(gè)待識(shí)別變量的反演識(shí)別結(jié)果;包括其后驗(yàn)分布及統(tǒng)計(jì)特征值.算例2 選取后10000 組樣本進(jìn)行統(tǒng)計(jì)分析,獲得6 個(gè)待求變量的反演識(shí)別結(jié)果,包括其后驗(yàn)分布及統(tǒng)計(jì)特征值.

由圖6 和圖7 可知,待識(shí)別變量后驗(yàn)分布均存在明顯的峰值,這說(shuō)明DREAM(D)算法能夠?qū)崿F(xiàn)對(duì)待求變量的有效反演識(shí)別.此外,可以明顯看出,算例1的2 個(gè)待識(shí)別變量的真實(shí)值均落在其后驗(yàn)分布的高密度區(qū)域;同時(shí),待識(shí)別變量后驗(yàn)分布的最大后驗(yàn)概率(MAP)值與其真實(shí)值很接近.這說(shuō)明,DREAM(D)算法能夠得到較高精度的反演識(shí)別結(jié)果.算例2 的污染源初始釋放時(shí)間(τ)和 3 個(gè)分區(qū)滲透系數(shù)(K 1, K 2, K 3)這4 個(gè)待識(shí)別變量的真實(shí)值均落在其后驗(yàn)分布的高密度區(qū)域;同時(shí),待識(shí)別變量的后驗(yàn)分布的MAP 值與其真實(shí)值很接近,這說(shuō)明,DREAM(D)算法對(duì)這4 個(gè)變量得到良好識(shí)別.污染源釋放強(qiáng)度(q)和縱向彌散度(D)這4 個(gè)待識(shí)別變量的真實(shí)值落在其后驗(yàn)分布的低密度區(qū)域,這2 個(gè)變量未能得到良好識(shí)別.

圖7 算例2 待識(shí)別變量的采樣迭代過(guò)程及其反演識(shí)別結(jié)果Fig.7 Sampling iteration process and identification results of unknown variables for case 2

為了進(jìn)一步提高污染源反演識(shí)別精度,采用自適應(yīng)更新多保真度代理模型策略,即將上一輪污染源反演結(jié)果后驗(yàn)分布的MAP 值代入HF 模型,獲得新的輸入-輸出訓(xùn)練樣本,訓(xùn)練新的多保真度Co-Kriging 模型,進(jìn)而進(jìn)行下一輪污染源反演識(shí)別.分別對(duì)算例1 和算例2 進(jìn)行2 次和3 次自適應(yīng)更新,待識(shí)別變量反演結(jié)果如圖8 和圖9 所示.

圖8 算例1 待識(shí)別變量反演識(shí)別結(jié)果Fig.8 Identification results of unknown variables for case 1

概率密度曲線的MAP 值越靠近真實(shí)值,說(shuō)明反演識(shí)別結(jié)果越好.由圖8 可以看出,經(jīng)過(guò)1 次自適應(yīng)更新后,待識(shí)別變量的MAP 值與其真實(shí)值更加接近,表明自適應(yīng)更新策略能夠顯著提升污染源反演精度.經(jīng)過(guò)2 次自適應(yīng)更新后的污染源反演結(jié)果與經(jīng)過(guò)1 次自適應(yīng)更新后基本一樣,沒(méi)有顯著差別,表明自適應(yīng)更新策略有助于得到高精度且穩(wěn)定的污染源反演結(jié)果.由圖9可以看出,經(jīng)過(guò)3次自適應(yīng)更新后,污染源釋放強(qiáng)度(q)、滲透系數(shù)(K2)和縱向彌散度(D)的MAP值與其真實(shí)值更加接近,表明自適應(yīng)更新策略能夠顯著提升反演精度.滲透系數(shù)(K1,K2)的MAP 值與更新前無(wú)明顯差別,初始釋放時(shí)間(τ)的MAP 值稍微遠(yuǎn)離其真實(shí)值,表明自適應(yīng)更新策略未能夠提升反演精度,原因在于模型輸出值(污染物濃度)對(duì)初始釋放時(shí)間(τ)不敏感和DREAM(D)算法進(jìn)行參數(shù)識(shí)別本身存在一定隨機(jī)性,故而難以提升反演識(shí)別精度.綜合算例1和算例2,采用自適應(yīng)更新多保真度代理模型策略能夠進(jìn)一步提升污染源反演識(shí)別精度.

圖9 算例2 待識(shí)別變量反演識(shí)別結(jié)果Fig.9 Identification results of unknown variables for case 2

因此,基于Co-Kriging 代理模型和DREAM(D)算法的地下水污染源反演識(shí)別是可行有效的.然而,本文構(gòu)建的Co-Kriging 代理模型只適合于當(dāng)前水文地質(zhì)條件和模型時(shí)空離散,當(dāng)水文地質(zhì)條件和模型時(shí)空離散變化時(shí),需要重新構(gòu)建Co-Kriging 代理模型.

3 結(jié)論

3.1 相比僅基于高保真度模型輸入-輸出樣本構(gòu)建的Kriging 代理模型,聯(lián)合運(yùn)用高保真度和低保真度模型輸入-輸出樣本構(gòu)建的Co-Kriging 代理模型對(duì)模擬模型的逼近精度更高.

3.2 耦合多保真度 Co-Kriging 代理模型和MCMC-DREAM(D)算法能夠得到較高精度的污染源反演結(jié)果,且能夠大幅度減小計(jì)算負(fù)荷;同時(shí),采用自適應(yīng)更新多保真度代理模型策略能夠進(jìn)一步提升污染源反演識(shí)別精度.