基于改進(jìn)稀疏網(wǎng)格替代模擬的地下水DNAPLs運(yùn)移不確定性分析

高鑫宇，曾獻(xiàn)奎，吳吉春

(南京大學(xué)地球科學(xué)與工程學(xué)院/表生地球化學(xué)教育部重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210023)

近年來(lái)，隨著石油化工產(chǎn)品和農(nóng)藥的廣泛使用，大量有機(jī)污染物進(jìn)入地下水中，造成地下水資源的污染。非水相流體(Non-Aqueous Phase Liquids, NAPLs)有機(jī)污染物是當(dāng)前地下水污染研究的重要對(duì)象[1]，根據(jù)其密度相對(duì)水的密度大小，分為輕非水相流體(LNAPLs)和重非水相流體(DNAPLs)。DNAPLs難溶于水，且密度比水大，可能附著于含水層而成為地下水的長(zhǎng)期污染源[2-3]。DNAPLs在地下水中存在水相、油相和氣相，它們之間相互轉(zhuǎn)化構(gòu)成多相流體系統(tǒng)。由于DNAPLs運(yùn)移過(guò)程復(fù)雜、處理費(fèi)用高昂，通過(guò)數(shù)值模擬技術(shù)定量描述DNAPLs的分布和遷移轉(zhuǎn)化，成為當(dāng)前DNAPLs場(chǎng)地污染修復(fù)與風(fēng)險(xiǎn)評(píng)價(jià)的重要手段[4-6]。

由于地質(zhì)環(huán)境的復(fù)雜性(如含水介質(zhì)的非均質(zhì)性、地層結(jié)構(gòu)的分布等)、有限的觀測(cè)數(shù)據(jù)以及DNAPLs的復(fù)雜運(yùn)移規(guī)律，導(dǎo)致地下水DNAPLs運(yùn)移模型的參數(shù)、結(jié)構(gòu)難以進(jìn)行有效識(shí)別，從而不能精確刻畫其遷移轉(zhuǎn)化過(guò)程，模擬結(jié)果具有不確定性。地下水模型的不確定性分析成為DNAPLs運(yùn)移定量刻畫的重要前提，基于貝葉斯理論的MCMC (Markov Chain Monte Carlo)模擬是當(dāng)前主要的不確定性定量分析方法。Liu[7]應(yīng)用MCMC對(duì)DNAPLs溶解及運(yùn)移的半解析解模型進(jìn)行了參數(shù)估計(jì)，并用于DNAPLs污染修復(fù)模型的污染源識(shí)別[8]。不確定性分析通常需要重復(fù)調(diào)用DNAPLs運(yùn)移模型幾千至幾萬(wàn)次，模型運(yùn)行時(shí)間成為當(dāng)前DNAPLs運(yùn)移模擬不確定性分析的重要問(wèn)題。替代模型是當(dāng)前克服不確定性分析計(jì)算耗時(shí)問(wèn)題的有效途徑，即通過(guò)建立與原始模型具有相同或接近模擬精度的模型代替原始模型。DNAPLs模擬中常見(jiàn)的替代模型技術(shù)有支持向量機(jī)[9-10]、克里金插值[11-12]、徑向基函數(shù)人工神經(jīng)網(wǎng)絡(luò)[9]和核極限學(xué)習(xí)機(jī)[13]等。

通過(guò)MCMC進(jìn)行不確定性分析時(shí)，需要不斷運(yùn)行模型計(jì)算似然函數(shù)，MCMC結(jié)果對(duì)似然函數(shù)高度敏感。因此，MCMC對(duì)替代模型的精度要求很高，當(dāng)前研究中尚未見(jiàn)基于似然函數(shù)替代模擬的DNAPLs運(yùn)移不確定性分析。本次研究擬采用稀疏網(wǎng)格(sparse grid, SG)替代模型技術(shù)，該方法適用于構(gòu)建具有復(fù)雜響應(yīng)關(guān)系的替代模型。SG替代模型構(gòu)建過(guò)程中，為提高構(gòu)建效率、減少參數(shù)空間網(wǎng)格冗余，一般有兩種自適應(yīng)方式，局部自適應(yīng)(local adaptive, LA)和維數(shù)自適應(yīng)(dimension adaptive, DA)。如Zhang等[14]通過(guò)一種高階LA-SG技術(shù)建立了鈾素反應(yīng)性溶質(zhì)運(yùn)移替代模型，并通過(guò)MCMC識(shí)別得到模型參數(shù)后驗(yàn)分布。Zeng等[15]基于DA-SG技術(shù)建立了溶質(zhì)運(yùn)移替代模型，并通過(guò)MCMC進(jìn)行了地下水污染源識(shí)別。替代模型的構(gòu)建成本(即需運(yùn)行一定數(shù)量的原始模型)與替代對(duì)象的復(fù)雜度密切相關(guān)，且隨著模型參數(shù)維數(shù)的增加而顯著增加。Zeng[15]等利用DA-SG建立5維溶質(zhì)運(yùn)移替代模型需運(yùn)行原模型1 545次。Zhang等[14]通過(guò)高階LA-SG建立6維反應(yīng)性溶質(zhì)運(yùn)移替代模型的成本是運(yùn)行原模型3 765次。Zeng等[16]使用LA-SG及約1萬(wàn)次原模型運(yùn)行建立了10維非均質(zhì)地下水流替代模型。因此，高效率、高精度的替代模型技術(shù)對(duì)于復(fù)雜、高維DNAPLs運(yùn)移模擬的不確定性分析具有重要的理論和實(shí)際意義。

本次研究擬基于SG技術(shù)，將替代模型參數(shù)空間的兩種網(wǎng)格自適應(yīng)生成方式(LA、DA)相結(jié)合，進(jìn)一步提高替代模型的構(gòu)建效率和替代精度。此外，通過(guò)兩個(gè)解析案例(5維rastrign函數(shù)、20維banana函數(shù))和一個(gè)DNAPLs運(yùn)移室內(nèi)模擬模型，驗(yàn)證所提出的維數(shù)、局部自適應(yīng)稀疏網(wǎng)格替代模型(DA-LA-SG)的效率和精度。

1 研究方法

1.1 稀疏網(wǎng)格替代模型

1.1.1基本原理

稀疏網(wǎng)格(SG)是一種基于Smolyak[17]規(guī)則的疊層拉格朗日插值方法，具有方法原理簡(jiǎn)單、插值精度高等優(yōu)點(diǎn)，在水資源與水環(huán)境模擬中應(yīng)用廣泛。對(duì)于一維非線性函數(shù)f(x)，其疊層拉格朗日插值項(xiàng)U定義為：

(1)

m0=1,m1=2,

mi=2i-1ifi≥2

(2)

SG插值點(diǎn)的生成規(guī)則有Clenshaw-Curits規(guī)則、Gauss-Patterson規(guī)則、Uniform規(guī)則等[18]，本研究中使用Uniform規(guī)則，在[0,1]間的坐標(biāo)表達(dá)式為：

(3)

(4)

式中：L——最大插值級(jí)數(shù)；

mi——對(duì)應(yīng)級(jí)數(shù)i的稀疏網(wǎng)格插值點(diǎn)的數(shù)量；

φ(x)——插值基函數(shù)，本次研究使用三次多項(xiàng)式基函數(shù)[18]；

c——插值系數(shù)，也叫基函數(shù)的疊層盈余，表示插值公式與目標(biāo)函數(shù)在x處的差值。

當(dāng)函數(shù)f(x)相對(duì)平滑時(shí)，隨著插值級(jí)數(shù)i增加，疊層盈余將接近于零，所以c可以作為局部自適應(yīng)的誤差標(biāo)準(zhǔn)。

基于一維疊層插值方法，多維稀疏網(wǎng)格插值為：

i=(i1，…，iD)

(5)

式中：D——模型參數(shù)x=(x1,…,xD) 的維度；

i——ΔVi,D[f(x)]的多維插值級(jí)數(shù)向量；

ik(k=1,…,D)——第k維插值級(jí)數(shù)。

ΔVi,D[f(x)]定義為：

ΔVi,D[f(x)]=ΔUi1?…?ΔUiD[f(x)]=

(6)

式中：j——向量(jl=1,…,mil,l=1,…,D)；

1.1.2自適應(yīng)稀疏網(wǎng)格(Adaptive SG)

普通SG方法的插值級(jí)數(shù)向量的集合為：β={i∈Nd:λ(i)≤L}。相對(duì)于全張量網(wǎng)格節(jié)點(diǎn)，緩解了“維數(shù)詛咒”效應(yīng)，但對(duì)于高維參數(shù)空間，替代模型的成本仍不可忽視。根據(jù)目標(biāo)函數(shù)參數(shù)和輸出之間的響應(yīng)關(guān)系，通過(guò)網(wǎng)格自適應(yīng)的生成方式，減少模型參數(shù)空間的冗余網(wǎng)格節(jié)點(diǎn)，可以進(jìn)一步降低SG替代模型的成本。

(1)維數(shù)自適應(yīng)稀疏網(wǎng)格(DA-SG)

DA-SG的基本思路[19]可以總結(jié)為：① 從級(jí)數(shù)向量(0, 0, …, 0)開(kāi)始，在當(dāng)前級(jí)數(shù)向量i的前鄰域{i+ej: 1 ≤j≤D}(ej表示第j維為1、其他為0的單位向量)中搜索候選級(jí)數(shù)向量A，候選級(jí)數(shù)向量條件為i-ej∈Ofor1≤j≤D,ij>1，O為已有級(jí)數(shù)向量集合；②在A中找到使測(cè)試點(diǎn)的誤差項(xiàng)趨于最大的級(jí)數(shù)向量k，并將其加入已有級(jí)數(shù)向量集合O中；③將k作為當(dāng)前已有級(jí)數(shù)向量，重復(fù)步驟①和②，直到誤差低于預(yù)定義的誤差標(biāo)準(zhǔn)。這種方法結(jié)合了貪婪算法的思想，可在構(gòu)建稀疏網(wǎng)格的同時(shí)確定對(duì)目標(biāo)函數(shù)變化有最大影響的級(jí)數(shù)向量，提高了構(gòu)建替代模型的效率。

圖1所示為DA-SG插值節(jié)點(diǎn)的自適應(yīng)生成過(guò)程。圖1(a)～(c)所示為當(dāng)前插值級(jí)數(shù)向量，圖1(d)～(f)表示相應(yīng)級(jí)數(shù)向量對(duì)應(yīng)的網(wǎng)格點(diǎn)分布圖。圖1(a)中候選級(jí)數(shù)向量A包含(4，0)、(1，1)和(0，3)。首先，將A中誤差項(xiàng)最大的級(jí)數(shù)向量i即(1，1)加入已有級(jí)數(shù)向量集合O中。然后，對(duì)該級(jí)數(shù)向量i的所有前鄰域級(jí)數(shù)向量(即(2，1)、(1，2))進(jìn)行候選條件測(cè)試，結(jié)果發(fā)現(xiàn)均滿足條件，將其加入集合A中，此時(shí)A包含4個(gè)級(jí)數(shù)向量。最后，重復(fù)上述過(guò)程，尋找誤差最大項(xiàng)，一直循環(huán)計(jì)算直至誤差達(dá)到標(biāo)準(zhǔn)。

圖1 二維維數(shù)自適應(yīng)稀疏網(wǎng)格(DA-SG)示意圖Fig.1 Diagrammatic sketch of two dimensional DA-SG

(2)局部自適應(yīng)稀疏網(wǎng)格(LA-SG)

圖2 二維局部自適應(yīng)稀疏網(wǎng)格(LA-SG)示意圖[18]Fig.2 Diagrammatic sketch of two dimensional LA-SG

(3)維數(shù)-局部自適應(yīng)稀疏網(wǎng)格(DA-LA-SG)

DA-SG方法通過(guò)識(shí)別對(duì)目標(biāo)函數(shù)變化有顯著影響的級(jí)數(shù)向量，提高了SG構(gòu)建效率，但當(dāng)目標(biāo)函數(shù)變化區(qū)域集中在參數(shù)空間較小區(qū)域時(shí)，DA-SG方法的效率較低。與DA-SG相比，LA-SG方法通過(guò)在非線性程度低的區(qū)域減少插值節(jié)點(diǎn)，克服了DA-SG的不足，但LA-SG在所有級(jí)數(shù)向量空間均構(gòu)造插值節(jié)點(diǎn)。因此，提出將DA、LA結(jié)合的改進(jìn)稀疏網(wǎng)格替代模型方法DA-LA-SG，能夠進(jìn)一步提到SG替代模型的構(gòu)建效率。如圖3所示，首先對(duì)目標(biāo)函數(shù)進(jìn)行維數(shù)自適應(yīng)，尋找誤差最大項(xiàng)；其次，在計(jì)算插值系數(shù)時(shí)進(jìn)行局部自適應(yīng)，判斷局部自適應(yīng)的標(biāo)準(zhǔn)，舍棄已滿足誤差標(biāo)準(zhǔn)節(jié)點(diǎn)的下級(jí)節(jié)點(diǎn)；最后，重復(fù)進(jìn)行上述過(guò)程，直至維數(shù)自適應(yīng)誤差達(dá)到預(yù)設(shè)的標(biāo)準(zhǔn)，結(jié)束替代模型運(yùn)行。

圖3 DA-LA-SG算法流程圖Fig.3 Algorithm flow chart of DA-LA-SG

1.2 貝葉斯不確定性分析

貝葉斯不確定性分析方法的基礎(chǔ)是貝葉斯公式，可表示[22]為：

(7)

式中：p(θ|y)——參數(shù)θ的后驗(yàn)概率分布；

p(θ)——θ的先驗(yàn)概率分布；

p(y|θ)——θ的似然函數(shù)；

常用的似然函數(shù)為高斯函數(shù)，其自然對(duì)數(shù)表達(dá)式為：

(8)

式中：Σ——觀測(cè)誤差的協(xié)方差矩陣；

|Σ|——行列式。

由于p(θ|y)難以直接計(jì)算，一般采用統(tǒng)計(jì)抽樣的方式間接獲取，如MCMC模擬。MCMC方法利用構(gòu)建馬爾科夫鏈來(lái)反演模型參數(shù)θ的后驗(yàn)分布，通過(guò)在θ的概率空間內(nèi)持續(xù)搜索，逐漸收斂于一個(gè)穩(wěn)定的后驗(yàn)分布π(θ)。本次研究中采用最新的抽樣算法DREAMzs[23]進(jìn)行參數(shù)不確定性研究，基本步驟如下：①基于先驗(yàn)信息，設(shè)置模型參數(shù)的先驗(yàn)分布和馬爾科夫鏈的起點(diǎn)[θi,i=1,…,N]，N表示平行馬爾科夫鏈的數(shù)量；②將式(8)作為似然函數(shù)，利用觀測(cè)數(shù)據(jù)y，通過(guò)DREAMzs抽樣得到參數(shù)后驗(yàn)樣本；③對(duì)采樣得到的結(jié)果(模型參數(shù)和輸出的后驗(yàn)分布)進(jìn)行統(tǒng)計(jì)分析，定量分析不確定性大小。

1.3 替代模型的精度評(píng)估

采用目標(biāo)函數(shù)和替代模型間的相關(guān)系數(shù)R2和相對(duì)均方根誤差NRMSE[24]作為替代模型精度的評(píng)價(jià)指標(biāo)。

(1)相關(guān)系數(shù)(R2)

(9)

式中：n——樣本個(gè)數(shù)；

yi——目標(biāo)函數(shù)的輸出值；

R2越接近于1，替代模型效果越好。

(2)相對(duì)均方根誤差(NRMSE)

(10)

式中：|f|max——目標(biāo)函數(shù)輸出最大絕對(duì)值。

2 算例分析

2.1 解析案例

2.1.15維rastrign函數(shù)

選擇5維rastrign函數(shù)作為本次替代模型的目標(biāo)函數(shù)，用于測(cè)試DA-SG、LA-SG和DA-LA-SG等3個(gè)替代模型的效率和精度。Rastrign函數(shù)的解析表達(dá)式為：

其中，A=10,-5.12 ≤xi≤ 5.12,n=5。

圖4 (a)Rastrign函數(shù)在2維參數(shù)空間上的響應(yīng)面; (b)各替代模型NRMSE隨樣本數(shù)量的變化Fig.4 (a)Response surface of the Rastrign function in 2-D space; (b)The NRMSE of surrogates vary with the number of samples

圖4(a)為目標(biāo)函數(shù)在二維空間的響應(yīng)面示意圖。圖4(b)顯示DA-SG、LA-SG和DA-LA-SG等3個(gè)替代模型隨著插值節(jié)點(diǎn)數(shù)量的增多，替代誤差(NRMSE)均逐漸降低。從圖中可以看出，在替代模型建立的初期(插值節(jié)點(diǎn)少于100)，LA-SG的替代效率(達(dá)到一定替代精度所需的構(gòu)建成本)高于DA-SG和DA-LA-SG，但是插值節(jié)點(diǎn)數(shù)量的增多，DA-SG、DA-LA-SG的替代效率逐漸高于LA-SG。這是由于在替代模型的構(gòu)建初期，插值節(jié)點(diǎn)分布較少，LA-SG能夠更有效識(shí)別目標(biāo)函數(shù)的局部特征，而隨著插值節(jié)點(diǎn)數(shù)量增多，DA-SG方法能夠有效識(shí)別對(duì)目標(biāo)函數(shù)變化有顯著影響的級(jí)數(shù)向量，進(jìn)行插值空間的降維處理。當(dāng)插值節(jié)點(diǎn)數(shù)量低于2 600時(shí)，DA-SG與DA-LA-SG替代模型表現(xiàn)基本相同，這是由于替代模型前期插值系數(shù)(誤差)較大，達(dá)不到局部自適應(yīng)誤差標(biāo)準(zhǔn)。當(dāng)插值節(jié)點(diǎn)數(shù)量超過(guò)2 600，DA-LA-SG表現(xiàn)優(yōu)于DA-SG，DA-LA-SG中的部分插值節(jié)點(diǎn)因局部自適應(yīng)而被舍棄。因此，在替代成本相同的條件下，DA-LA-SG的精度高于DA-SG。

2.1.220維banana函數(shù)

為了進(jìn)一步測(cè)試高維情況下DA-SG、LA-SG和DA-LA-SG等3個(gè)替代模型的效率和精度，選擇20維banana函數(shù)作為本次替代模型的目標(biāo)函數(shù)，banana函數(shù)的解析表達(dá)式為：

-10≤xi≤10,1≤i≤n

圖5(a)為目標(biāo)函數(shù)在二維參數(shù)空間的響應(yīng)面示意圖。圖5(b)所示為DA-SG、LA-SG和DA-LA-SG等3個(gè)替代模型隨著插值節(jié)點(diǎn)數(shù)量的增多，替代誤差(NRMSE)逐漸降低的過(guò)程。從圖中可以看出，在替代模型的建立初期(插值節(jié)點(diǎn)少于1 000)，LA-SG的替代效率比DA-SG和DA-LA-SG高，但隨著插值節(jié)點(diǎn)數(shù)量的增多，DA-SG、DA-LA-SG的替代效率逐漸比LA-SG高。這是由于在替代模型的構(gòu)建初期，插值節(jié)點(diǎn)較少，LA-SG通過(guò)識(shí)別目標(biāo)函數(shù)的局部特征進(jìn)行插值節(jié)點(diǎn)的自適應(yīng)，而隨著插值節(jié)點(diǎn)的增多，DA-SG通過(guò)識(shí)別重要的級(jí)數(shù)向量進(jìn)行降維處理，其替代效率逐漸升高。插值節(jié)點(diǎn)數(shù)量低于2 000時(shí)，DA-SG與DA-LA-SG表現(xiàn)基本相同，這是由于替代模型前期的插值誤差較大，局部自適應(yīng)過(guò)程很少發(fā)生。當(dāng)插值節(jié)點(diǎn)數(shù)量超過(guò)2 000，DA-LA-SG表現(xiàn)優(yōu)于DA-SG，這是由于局部自適應(yīng)DA-LA-SG舍棄了部分冗余插值節(jié)點(diǎn)，從而提高了替代效率。

圖5 (a) Banana函數(shù)在2維空間上的響應(yīng)面；(b)各替代模型NRMSE隨樣本數(shù)量的變化Fig.5 (a) Response surface of the banana function in 2-D space; (b) The NRMSE of surrogates vary with the number of samples

2.2 二維砂箱DNAPLs運(yùn)移室內(nèi)實(shí)驗(yàn)

本次實(shí)驗(yàn)選用四氯乙烯(PCE)作為典型DNAPLs污染物，通過(guò)向砂箱中注入PCE的方式模擬污染物的泄露。圖6(a)為二維砂箱多孔介質(zhì)中PCE的運(yùn)移模擬室內(nèi)實(shí)驗(yàn)示意圖。砂箱的長(zhǎng)、寬、厚分別為0.6 m×0.45 m×0.016 m。砂箱左側(cè)為入流邊界，流速保持7.2×10-3m3/d，由蠕動(dòng)泵控制。砂箱背景介質(zhì)為均質(zhì)多孔介質(zhì)，內(nèi)部分布5個(gè)低滲透介質(zhì)作為透鏡體模擬介質(zhì)非均質(zhì)性對(duì)PCE運(yùn)移的影響。砂箱實(shí)驗(yàn)中，PCE注入點(diǎn)位于x=0.3 m，y=0.4 m，z=0.008 m處。通過(guò)光透法[25](LTM)實(shí)時(shí)監(jiān)測(cè)砂箱中的PCE飽和度，且用作本次實(shí)驗(yàn)的觀測(cè)數(shù)據(jù)。此外，本次案例研究采用T2VOC程序[26]對(duì)污染源注入過(guò)程的PCE運(yùn)移進(jìn)行數(shù)值模擬，PCE注入時(shí)間共持續(xù)80 min。圖6(b)顯示PCE停止注入時(shí)刻(即t=80 min)，砂箱中PCE的飽和度分布。

圖6 (a)二維砂箱PCE運(yùn)移實(shí)驗(yàn)示意圖；(b)停止注入時(shí)刻PCE飽和度分布圖Fig.6 (a)Diagrammatic sketch of the 2-D sandbox PCE experiment; (b)Saturation distribution of PCE after stopping injection

在實(shí)際情況下，由于有限的觀測(cè)信息，一些水文地質(zhì)參數(shù)難以準(zhǔn)確給定，例如滲透率、孔隙度等。對(duì)于PCE運(yùn)移模擬模型，PCE的濃度分布對(duì)含水層的滲透性十分敏感[27]。因此，本次研究中選擇兩個(gè)區(qū)域(背景介質(zhì)和低滲透體)的滲透率(k1,k2)和孔隙度(μ1,μ2) 作為未知參數(shù)，同時(shí)將T2VOC程序中相對(duì)滲透率函數(shù)中的擬合系數(shù)(n1,n2)作為未知參數(shù)。通過(guò)觀測(cè)數(shù)據(jù)(PCE的飽和度分布)，對(duì)本次PCE運(yùn)移模擬的6個(gè)未知參數(shù)進(jìn)行識(shí)別。

將本次砂箱實(shí)驗(yàn)處理為二維PCE運(yùn)移問(wèn)題，模擬區(qū)域在x,y,z方向上離散為35×38×1個(gè)網(wǎng)格單元。PCE運(yùn)移數(shù)值模型中，采用Parker本構(gòu)模型。此外，對(duì)PCE運(yùn)移模型進(jìn)行參數(shù)不確定性分析時(shí)，未知模型參數(shù)的先驗(yàn)分布均假定為均勻分布，范圍見(jiàn)表1。

表1 PCE運(yùn)移模型中未知參數(shù)的取值范圍

2.2.1建立似然函數(shù)的替代模型

本次案例分析中，分別采用LA-SG、DA-SG和DA-LA-SG等3種方法建立lgp(y|θ)的替代模型，并對(duì)這三種替代方法的效率和精度進(jìn)行對(duì)比分析。其中，替代效率表示建立替代模型的成本，替代精度為替代模型與原始T2VOC結(jié)果的接近程度。

圖7 PCE運(yùn)移模擬對(duì)數(shù)似然函數(shù)替代模型的NRMSEFig.7 NRMSE of surrogates for the logarithmic likelihood function

圖7為砂箱PCE運(yùn)移模擬對(duì)數(shù)似然函數(shù)的替代模型(DA-SG、LA-SG和DA-LA-SG)隨著插值節(jié)點(diǎn)數(shù)量的增多，替代精度逐漸升高的過(guò)程。從圖中可以看出，在替代模型的構(gòu)建初期(插值節(jié)點(diǎn)數(shù)量小于等于100)，3個(gè)替代模型具有相接近的替代效率，替代模型的NRMSE隨插值點(diǎn)數(shù)量的衰減速度相似。隨著樣本數(shù)量的增多，替代模型的替代精度均逐漸升高，但NRMSE的衰減速度與替代模型初期相比，均有所降低。DA-LA-SG具有最高的替代效率，即在相同的替代成本條件下，DA-LA-SG具有最高的替代精度，其次為DA-SG和LA-SG。在替代模型的構(gòu)建初期，各參數(shù)維度上的插值節(jié)點(diǎn)數(shù)量較少，進(jìn)行維數(shù)自適應(yīng)時(shí)，所舍棄的插值節(jié)點(diǎn)并不明顯。隨著插值水平(Level)的增加，各參數(shù)維度所需的插值節(jié)點(diǎn)數(shù)量迅速增加，進(jìn)行維數(shù)自適應(yīng)能夠顯著減少插值節(jié)點(diǎn)的數(shù)量，提高替代模型的效率。

圖8顯示了本次PCE運(yùn)移模型的對(duì)數(shù)似然函數(shù)lgp(y|θ)的響應(yīng)面及其替代模型的效果。圖8(a)～(c)分別表示lgp(y|θ)在不同2維參數(shù)空間上的響應(yīng)面，(d)表示采用DA-LA-SG方法建立的lgp(y|θ)替代模型的效果。結(jié)果顯示，建立該替代模型共需要運(yùn)行原始PCE運(yùn)移模型745次，原始模型與替代模型間的相關(guān)系數(shù)R2為0.999 7，替代模型的RMSE為3.63，表明通過(guò)DA-LA-SG方法建立的lgp(y|θ)替代模型具有較高的替代效率和精度。本次室內(nèi)砂箱實(shí)驗(yàn)PCE運(yùn)移模擬參數(shù)不確定性分析中，采用DA-LA-SG建立的似然函數(shù)替代模型，而無(wú)需調(diào)用原始的PCE運(yùn)移模型。

圖8 (a)～(c)對(duì)數(shù)似然函數(shù)在二維參數(shù)空間上的響應(yīng)面；(d) DA-LA-SG替代模型效果圖Fig.8 (a)～(c) Response surfaces of the 2-D logarithmic likelihood function; (d) Performance of DA-LA-SG

2.2.2基于DA-LA-SG替代模型的PCE運(yùn)移模擬參數(shù)不確定性分析

圖9 PCE運(yùn)移模型參數(shù)的后驗(yàn)概率分布Fig.9 Posterior probability density of the PCE migration parameters

采用基于DREAMzs算法的MCMC模擬識(shí)別6個(gè)未知模型參數(shù)的后驗(yàn)分布，參數(shù)先驗(yàn)分布見(jiàn)表1。為了保證搜索效率，使用3條平行馬爾科夫鏈，每條鏈長(zhǎng)度為30 000，預(yù)熱期長(zhǎng)度為10 000，剩下的樣本用于生成參數(shù)的后驗(yàn)分布。在MCMC運(yùn)行過(guò)程中，調(diào)用基于DA-LA-SG建立的似然函數(shù)替代模型不需運(yùn)行原始模型。圖9為識(shí)別得到的6個(gè)模型參數(shù)的后驗(yàn)概率分布，坐標(biāo)橫軸區(qū)間表示相應(yīng)參數(shù)的先驗(yàn)分布范圍。經(jīng)過(guò)MCMC參數(shù)識(shí)別，k1、k2和n1的后驗(yàn)概率密度呈正態(tài)分布，相對(duì)于先驗(yàn)分布，參數(shù)的后驗(yàn)分布范圍明顯縮小，說(shuō)明這3個(gè)參數(shù)對(duì)觀測(cè)數(shù)據(jù)(PCE的飽和度)比較敏感。相比之下，模型參數(shù)μ1、μ2和n2后驗(yàn)概率分布大致呈均勻分布，且與先驗(yàn)分布范圍相同，因此這3個(gè)參數(shù)對(duì)觀測(cè)數(shù)據(jù)不夠敏感，參數(shù)識(shí)別效果較差。本次DNAPL運(yùn)移模擬不確定性分析結(jié)果表明，砂箱實(shí)驗(yàn)中DNAPL的運(yùn)移主要受介質(zhì)的滲透性影響[28-29]，DNAPL在運(yùn)移過(guò)程中受到弱透鏡體的阻礙，由于無(wú)法克服毛細(xì)壓力而在其上方聚積并側(cè)向運(yùn)移，孔隙度對(duì)DNAPL運(yùn)移的影響較小。因此，對(duì)于實(shí)際DNAPL污染場(chǎng)地的評(píng)估與修復(fù)，場(chǎng)地含水層的滲透性質(zhì)特別是弱透水介質(zhì)分布的精確刻畫，是建立準(zhǔn)確、可靠DNAPLS運(yùn)移模擬模型的重要條件。

3 結(jié)論

為克服DNAPLs運(yùn)移模擬不確定性分析中的計(jì)算耗時(shí)問(wèn)題，本次研究在傳統(tǒng)稀疏網(wǎng)格(SG)替代模型的基礎(chǔ)上，提出了一種將局部自適應(yīng)(LA)和維數(shù)自適應(yīng)(DA)兩種網(wǎng)格生成方式耦合的改進(jìn)替代模型DA-LA-SG。通過(guò)兩個(gè)解析案例和一個(gè)DNAPL運(yùn)移室內(nèi)實(shí)驗(yàn)，對(duì)比分析了DA-SG、LA-SG和DA-LA-SG三種替代模型技術(shù)的表現(xiàn)，并將DA-LA-SG替代模型用于DNAPL運(yùn)移數(shù)值模擬不確定性分析中，得到以下結(jié)論：

(1)在替代模型的構(gòu)建初期，LA-SG的替代效率(一定替代精度條件下的替代成本)優(yōu)于或接近于DA-SG和DA-LA-SG，隨著插值節(jié)點(diǎn)數(shù)量的增多，DA-SG和DA-LA-SG的替代效率逐漸優(yōu)于LA-SG，且DA-LA-SG具有最高的替代效率。

(2)使用改進(jìn)的維數(shù)、局部自適應(yīng)稀疏網(wǎng)格技術(shù)(DA-LA-SG)能夠高效、高精度的建立PCE運(yùn)移模型的似然函數(shù)替代模型。構(gòu)建替代模型需運(yùn)行原始模型745次，原始模型與替代模型間的相關(guān)系數(shù)R2達(dá)到0.999 7。

(3)將建立的似然函數(shù)替代模型用于PCE運(yùn)移模擬參數(shù)不確定性分析，克服了因重復(fù)調(diào)用PCE運(yùn)移模型導(dǎo)致的計(jì)算耗時(shí)問(wèn)題，通過(guò)MCMC模擬獲得了模型參數(shù)的后驗(yàn)分布。結(jié)果發(fā)現(xiàn)，背景介質(zhì)滲透率k1、擬合系數(shù)n1和透鏡體滲透率k2的可識(shí)別性較強(qiáng)，背景介質(zhì)孔隙度μ1、透鏡體孔隙度μ2和擬合系數(shù)n2的識(shí)別效果較差，這些參數(shù)對(duì)飽和度觀測(cè)數(shù)據(jù)不敏感。