周念清，張瑞城，江思珉，夏學(xué)敏

(1.同濟(jì)大學(xué)土木工程學(xué)院，上海 200092；2.上海理工大學(xué)環(huán)境與建筑學(xué)院，上海 200093)

地下水資源由于具有分布廣泛、水質(zhì)良好、不易受污染等特點，在人類日常生活和工農(nóng)業(yè)生產(chǎn)中被廣泛利用。然而，隨著城市化進(jìn)程的加快和現(xiàn)代工業(yè)的迅速發(fā)展，越來越多的污染物因泄漏或排放不當(dāng)而進(jìn)入地下水體中，導(dǎo)致地下水污染問題日趨嚴(yán)重。為了有效地治理和修復(fù)地下水污染，首先需要對污染源信息(包括污染源強(qiáng)度和位置等)進(jìn)行準(zhǔn)確識別。傳統(tǒng)的地下水監(jiān)測采樣方法由于成本較高且獲取的觀測數(shù)據(jù)較少，往往難以直接獲取污染源信息以及場地的水文地質(zhì)參數(shù)[1-2]。現(xiàn)階段的研究多通過構(gòu)建地下水?dāng)?shù)值模擬模型并用求解逆問題的方式對污染源進(jìn)行參數(shù)反演[3-4]。

參數(shù)反演法是通過分析多源觀測信息，并與已知的參照值進(jìn)行比較，不斷迭代更新模型參數(shù)，使正演模型輸出的觀測值與參照值不斷趨近，從而得到未知模型參數(shù)的近似估計，將其應(yīng)用于地下水污染溯源問題中，可實現(xiàn)對未知污染源參數(shù)的模擬預(yù)測[5]。近年來，貝葉斯推斷已被廣泛地應(yīng)用于地下水參數(shù)反演問題之中[6-8]。其中，基于馬爾科夫鏈蒙特卡羅(Markov chain Monte Carlo，MCMC)方法和基于集合的數(shù)據(jù)同化方法(ensemble-based data assimilation)是貝葉斯推斷的兩種常見應(yīng)用[9-11]。MCMC方法通過構(gòu)造合適的馬爾科夫鏈進(jìn)行抽樣并使用蒙特卡洛方法進(jìn)行積分計算，最終收斂到平穩(wěn)分布的參數(shù)后驗分布，但使用MCMC方法通常需要大量調(diào)用系統(tǒng)模型才能達(dá)到收斂，并且模型的調(diào)用次數(shù)將隨著參數(shù)維度的增加而大幅增加，給計算帶來沉重的負(fù)擔(dān)[12]?；诩系臄?shù)據(jù)同化算法在求解參數(shù)反演問題時，需要調(diào)用系統(tǒng)模型的次數(shù)相對較少，因而比MCMC算法的計算效率更高，其中集合卡爾曼濾波[13](ensemble Kalman filter，EnKF)和集合平滑器[14](ensemble smoother，ES)是其典型代表，二者的區(qū)別在于EnKF需要同時更新輸入?yún)?shù)和狀態(tài)變量，而ES僅進(jìn)行參數(shù)更新，因此可以有效避免EnKF中參數(shù)和狀態(tài)的不一致問題[15]。為了提高ES方法的計算效率和對高維問題的適用性，Emerick等[16]提出了一種多重數(shù)據(jù)同化集合平滑器(ensemble smoother with multiple data assimilation,ES-MDA)，利用添加擾動后的觀測誤差協(xié)方差矩陣對觀測值進(jìn)行多次數(shù)據(jù)同化，提高了算法在求解高維參數(shù)反演問題時的性能。

在利用數(shù)值模型對地下水污染物遷移轉(zhuǎn)化過程進(jìn)行模擬預(yù)測的過程中，模型參數(shù)的準(zhǔn)確性是影響預(yù)測精度的關(guān)鍵。然而，由于模型參數(shù)存在空間變異性特點，對于非均質(zhì)性很大的污染場地，難以利用稀疏的觀測數(shù)據(jù)來準(zhǔn)確反演估計污染場地的模型參數(shù)[17]。為了解決傳統(tǒng)的勘測方法成本較高且獲取的觀測數(shù)據(jù)較為稀疏的問題，近年來有學(xué)者[18-19]將地球物理方法(如高密度電阻率法、地質(zhì)雷達(dá)法等)引入到污染場地調(diào)查中，能夠以較低的成本獲取大量的連續(xù)觀測數(shù)據(jù)，有效提高了參數(shù)反演方法在應(yīng)對高維參數(shù)空間時的適用性。

本文提出一種基于ES-MDA算法求解地下水污染溯源問題的數(shù)據(jù)同化方法，并融合由高密度電阻率(electrical resistivity tomography,ERT)法采集的ERT觀測數(shù)據(jù)進(jìn)行參數(shù)反演，然后將該反演結(jié)果與傳統(tǒng)的以濃度和水頭為觀測值方法而非采用地球物理方法得到的結(jié)果進(jìn)行比較，驗證該方法在求解高維參數(shù)反演問題時的有效性。

1 研究方法

1.1 地下水污染遷移模型

利用MODFLOW程序[20]與MT3DMS程序[21]構(gòu)建地下水污染遷移的數(shù)值模擬模型，其中，MODFLOW程序可以求解飽和多孔介質(zhì)中的地下水流問題，用基本微分方程表示為

(1)

式中：Kx、Ky、Kz分別為滲透系數(shù)在X、Y、Z方向上的分量，m/d；H為水頭，m；W為單位時間從單位體積含水層中流入或流出的水量，d-1；Ss為貯水率，m-1；t為時間，d。

在查明了地下水流動信息的基礎(chǔ)上，再結(jié)合MT3DMS程序，對地下水污染物的溶質(zhì)運(yùn)移過程進(jìn)行模擬，其基本控制方程為

(2)

式中：C為質(zhì)量濃度，mg/L；t為時間，d；Dij為水動力彌散系數(shù)，m2/d；vi為含水介質(zhì)中的實際水流速度，m/d；qs為單位體積含水層的源匯項流量，d-1；Cs為源匯項中溶質(zhì)質(zhì)量濃度，mg/L；∑Rn為化學(xué)反應(yīng)項總和，mg/(L·d)。

1.2 多重數(shù)據(jù)同化集合平滑器(ES-MDA)

ES-MDA算法是地下水污染溯源問題中一種常用的數(shù)據(jù)同化方法，其基本原理與ES算法類似，區(qū)別在于ES-MDA算法在運(yùn)算過程中利用觀測誤差的協(xié)方差矩陣對觀測值進(jìn)行多次數(shù)據(jù)同化，實現(xiàn)了對參數(shù)樣本的迭代更新，以更好地應(yīng)對高維參數(shù)空間反演問題[22-23]。ES-MDA算法的基本步驟如下：

第一步確定數(shù)據(jù)同化的迭代次數(shù)(Na)和每次迭代對應(yīng)的膨脹系數(shù)ai(其中，i=1,2,…，Na)。

第二步從模型參數(shù)的先驗分布中選取Ne個樣本，組成初始樣本集合，M1=[m1,1,m2,1,…，mNe,1]。

第三步從i=1開始執(zhí)行多次迭代計算，在第i次迭代中，對樣本集合Mi=[m1，i,m2,i,…,mNe,i]中的每一個樣本運(yùn)行正演模型以獲取該樣本對應(yīng)的預(yù)測值dj,i，正演模型的表達(dá)式如下：

dj,i=f(mj,i),j=1,2,…，Ne

(3)

第四步根據(jù)公式(4)對樣本集合進(jìn)行更新：

mj，i+1=mj，i+CMD(CDD+aiCD)-1(dobsj-dj,i),

j=1,2,…,Ne

(4)

式中：CMD為模型參數(shù)mi與預(yù)測值di之間的協(xié)方差矩陣;CDD為預(yù)測值的自協(xié)方差矩陣;CD為觀測誤差的協(xié)方差矩陣;dobsj為根據(jù)aiCD添加擾動后的觀測值。

第五步重復(fù)執(zhí)行第3步和第4步直到完成第Na次迭代，得到最終的參數(shù)樣本集合MNa=[m1,Na,m2,Na,…,mNe,Na]，并從中獲得對參數(shù)后驗的反演估計。

1.3 耦合地下水-地球物理模型的數(shù)據(jù)同化方法

根據(jù)巖石地球物理關(guān)系構(gòu)建聯(lián)結(jié)地球物理信息與水文地質(zhì)信息的數(shù)值模型，利用Archie公式[24]描述地層電阻率、溶液電阻率與含水飽和度之間關(guān)系，其公式基本形式如下：

(5)

式中：σt為地層總電阻率，Ω·m；σw為溶液電阻率，Ω·m；φ為有效孔隙度；Sw為含水飽和度；m為黏結(jié)指數(shù)；n為飽和度指數(shù)。

溶液電阻率與溫度和溶液質(zhì)量濃度等因素有關(guān)，Sen[17]綜合考慮了這2種因素對溶液電阻率的影響，提出了如下模型公式：

(6)

式中：C為溶液的質(zhì)量濃度，mg/L；T為溶液溫度，文中假定恒為25 ℃。將式(5)與式(6)聯(lián)合，便可構(gòu)建地層電阻率與溶液質(zhì)量濃度的數(shù)值關(guān)系。

高密度電阻率法(ERT法)作為一種常用的地球物理方法，可以通過供電電極(A、B)向地下空間輸入穩(wěn)定的電流I，然后根據(jù)測量電極(M、N)處的電位值V得到視電阻率的觀測數(shù)據(jù)ρa(bǔ)，再據(jù)此推算出電阻率的空間分布。利用開源的Pygimli程序[25]構(gòu)建地球物理模型并求解ERT正演問題。ERT法的正演模型可表示為

(7)

式中：ρ為地層電阻率的空間分布，Ω·m；V為電位值，V；I為電流值，A；δ(r)為狄拉克δ函數(shù)。

以上述方法為基礎(chǔ)，提出基于ES-MDA算法融合地球物理數(shù)據(jù)的數(shù)據(jù)同化方法，通過耦合地下水-地球物理模型，將質(zhì)量濃度數(shù)據(jù)轉(zhuǎn)換為ERT觀測值，然后利用ES-MDA算法融合ERT觀測數(shù)據(jù)對未知污染源參數(shù)進(jìn)行反演估計。該方法的流程圖見圖1。

圖1 基于ES-MDA算法的數(shù)據(jù)同化方法流程Fig.1 The flowchart of the data assimilation method based on ES-MDA algorithm

2 數(shù)值算例

圖2 污染場地示意圖Fig.2 The sketch map of contaminated site

如圖2所示，污染場地存在2個點源污染，污染物成分相同，根據(jù)場地調(diào)查情況確定位置為S1和S2，假定在模擬時間內(nèi)其質(zhì)量濃度保持不變，點源S1處的污染物質(zhì)量濃度為60 mg/L，點源S2處的污染物質(zhì)量濃度為30 mg/L。場地內(nèi)共布設(shè)8個水文監(jiān)測點(W1至W8)，用于測量質(zhì)量濃度與水頭數(shù)據(jù)，并在模擬期t為40、60、80、100、120、140、160 d時進(jìn)行數(shù)據(jù)采集，共獲得8個固定水頭數(shù)據(jù)與56個質(zhì)量濃度數(shù)據(jù)。同時，在地表布設(shè)1條水平測線(51個電極，間距2 m)與3條垂直方向的測線(位于x=24、48，72 m處)，每條測線包含20個間距為1 m的電極，共布設(shè)111個電極，并采用dipole-dipole的方式在相同的時間點測量電位數(shù)據(jù)(共有7次電流注入，電流為1 A)，共獲得763個ERT觀測數(shù)據(jù)?？紤]到觀測誤差對參數(shù)反演結(jié)果的影響，設(shè)置質(zhì)量濃度與水頭觀測值相對誤差為2%，ERT觀測值相對誤差為5%。

針對不同類型的觀測數(shù)據(jù)，設(shè)計了3組算例來對比其參數(shù)反演效果。3組算例均采用ES-MDA算法作為數(shù)據(jù)同化方法：算例1使用質(zhì)量濃度與水頭數(shù)據(jù)作為觀測數(shù)據(jù)；算例2僅使用ERT數(shù)據(jù)作為觀測數(shù)據(jù)；算例3使用質(zhì)量濃度、水頭和ERT數(shù)據(jù)作為觀測數(shù)據(jù)。3組算例觀測數(shù)據(jù)的采集時間相同，具體設(shè)置見表1。

表1 算例設(shè)置Tab.1 The setting of case studies

考慮到含水層的空間各向異性，在對滲透系數(shù)場進(jìn)行參數(shù)反演時計算代價往往較高。本算例中的滲透系數(shù)場根據(jù)模型的網(wǎng)格剖分維度為1 000維，屬于維度較高的反演問題，如果直接進(jìn)行求解效率會很低。因此，采用Karhunen-Loeve展開[6](KL展開)的方法對滲透系數(shù)場進(jìn)行降維。

(8)

式中：ζi為獨(dú)立標(biāo)準(zhǔn)高斯分布的隨機(jī)數(shù)；λi和si(x,y)為特征值和特征向量；NKL為KL展開保留的項數(shù)，本算例中取NKL=60，可以保留對數(shù)滲透系數(shù)場95%以上的變異性。

利用ES-MDA算法進(jìn)行參數(shù)反演的估計精度可以使用均方根誤差(ERMS)來量化，該值反映了模型參數(shù)的估計值與真實值之間的差距，ERMS的值越小(趨近于0)，參數(shù)反演的精度越高。ERMS的計算公式為

(9)

3 結(jié)果與討論

在3組算例中，潛在污染點源S1和S2的質(zhì)量濃度在模擬期內(nèi)保持不變，共有2個未知質(zhì)量濃度參數(shù)需要識別。滲透系數(shù)場采用KL展開的方式對其降維，其中KL展開項數(shù)NKL=60，因此對高維參數(shù)場(1 000維)的識別降維成對60個高斯分布的KL展開項的反演。由此得出，對于3組算例，總共需要進(jìn)行數(shù)據(jù)同化的未知模型參數(shù)數(shù)量均為62個。根據(jù)以往設(shè)計數(shù)值算例的經(jīng)驗，將ES-MDA算法的基礎(chǔ)參數(shù)設(shè)置為：樣本集合數(shù)Ne=1 000;迭代次數(shù)Na=7;膨脹系數(shù)ai=14、12、10、8、6、4、2(其中，i=1,2,…,Na)。

3.1 污染源源強(qiáng)的反演識別

圖3(a)至3(c)分別為3組算例對污染源源強(qiáng)的反演結(jié)果。圖3(a)是將質(zhì)量濃度與水頭數(shù)據(jù)作為觀測值得到的源強(qiáng)反演結(jié)果，經(jīng)過ES-MDA算法7次迭代后，污染源S1處的源強(qiáng)已基本收斂，污染源S2處的源強(qiáng)與真值仍有較大偏差，其主要原因是觀測井?dāng)?shù)量有限，導(dǎo)致采集到的觀測數(shù)據(jù)較少。圖3(b)是將ERT數(shù)據(jù)作為觀測值得到的源強(qiáng)反演結(jié)果，S1與S2處的源強(qiáng)參數(shù)均已基本收斂，且與算例1相比，反演結(jié)果更加趨近于真值，說明利用地球物理方法獲取的大量ERT觀測數(shù)據(jù)能有效改善傳統(tǒng)觀測方法獲取數(shù)據(jù)稀疏的不足。圖3(c)是將質(zhì)量濃度與水頭數(shù)據(jù)和ERT數(shù)據(jù)結(jié)合，共同作為觀測值得到的源強(qiáng)反演結(jié)果，經(jīng)過7次迭代后S1與S2處的源強(qiáng)參數(shù)明顯收斂，且其與真值的擬合程度均優(yōu)于算例1和算例2。

由此可得，融合ERT數(shù)據(jù)的ES-MDA算法可以更加精確地反演污染源源強(qiáng)，并且在此基礎(chǔ)上添加質(zhì)量濃度與水頭觀測數(shù)據(jù)，將使反演結(jié)果得到進(jìn)一步優(yōu)化，說明觀測信息的數(shù)量會直接影響參數(shù)反演的效果。利用地球物理方法獲取的ERT數(shù)據(jù)雖然測量精度不如質(zhì)量濃度與水頭數(shù)據(jù)，但可以便捷地獲取大量觀測信息，本案例中采集到的ERT觀測數(shù)據(jù)量(763個)遠(yuǎn)多于算例1中采集到的質(zhì)量濃度與水頭數(shù)據(jù)量(64個)。并且，在對污染源源強(qiáng)的反演識別上，通過融合多源觀測數(shù)據(jù)能夠有效提升參數(shù)反演的精度。然后，通過量化的方式進(jìn)一步比較使用不同類型的觀測數(shù)據(jù)反演污染源源強(qiáng)的效果，將3組算例污染源源強(qiáng)反演結(jié)果的ERMS值列于表2中。由表2可知，3組算例是針對任一個污染源的ERMS值均呈現(xiàn)逐漸縮小的趨勢，證明了融合多源觀測數(shù)據(jù)的算例3對污染源強(qiáng)度的反演精度要優(yōu)于算例1和算例2。

圖3 污染源源強(qiáng)的反演結(jié)果Fig.3 The inversion results of contaminant source strength

表2 污染源參數(shù)反演結(jié)果的均方根誤差Tab.2 Root-mean-square error of Case 1,Case 2 and Case 3

3.2 滲透系數(shù)場的反演識別

圖4、圖5和圖6分別為3組算例對滲透系數(shù)場的反演結(jié)果。作為參照場的lnK真實分布如圖4(a)、圖5(a)和圖6(a)；在反演開始階段生成的初始隨機(jī)場如圖4(b)、圖5(b)和圖6(b)，為了對比3組算例的反演效果，采用相同的初始隨機(jī)場。對滲透系數(shù)場進(jìn)行參數(shù)反演得到的后驗均值場如圖4(c)、圖5(c)和圖6(c)；反映集合離散程度的方差場如圖4(d)、圖5(d)和圖6(d)。

圖4 算例1滲透系數(shù)場的反演結(jié)果Fig.4 The inversion results of hydraulic conductivity in Case 1

如圖4所示，在使用質(zhì)量濃度與水頭數(shù)據(jù)作為觀測值的算例1中，其后驗均值場基本反映出lnK場的高、低值分布情況，但是對lnK場主體形態(tài)的描述還不夠精細(xì)，其反演結(jié)果與參照場仍存在不小的差距。在使用ERT數(shù)據(jù)作為觀測值的算例2(圖5)中，其后驗均值場不僅能反映出lnK場的高、低值分布區(qū)域，還能較好地刻畫滲透系數(shù)場的主體形態(tài)，因此，算例2與參照場的擬合程度要優(yōu)于算例1，但是在對滲透系數(shù)場分布細(xì)節(jié)的刻畫上仍然不夠精細(xì)。由圖5(d)可以得出算例2的估計方差已趨近于0，說明繼續(xù)進(jìn)行迭代反演對結(jié)果的提升空間有限。在算例3中，結(jié)合ERT數(shù)據(jù)和質(zhì)量濃度與水頭數(shù)據(jù)作為觀測值，得出的結(jié)果見圖6。從圖6可以明顯地看出后驗均值場在主體形態(tài)與細(xì)節(jié)上均較為精細(xì)地刻畫出了lnK場的空間分布情況，其與參照場的擬合程度要明顯優(yōu)于算例1和算例2。

圖5 算例2滲透系數(shù)場的反演結(jié)果Fig.5 The inversion results of hydraulic conductivity in Case 2

通過上述對比分析，在對高維滲透系數(shù)場反演識別上，融合ERT數(shù)據(jù)的ES-MDA算法可以得到更加精確的結(jié)果，并且在此基礎(chǔ)上融合質(zhì)量濃度與水頭數(shù)據(jù)作為觀測值，可以有效地提升滲透系數(shù)場的反演精度，說明利用多源觀測數(shù)據(jù)在求解參數(shù)反演問題上具有明顯的優(yōu)勢。

圖6 算例3滲透系數(shù)場的反演結(jié)果Fig.6 The inversion results of hydraulic conductivity in Case 3

4 結(jié)論與展望

提出了基于ES-MDA算法融合ERT數(shù)據(jù)的地下水參數(shù)反演方法，并通過數(shù)值算例研究，得出以下結(jié)論：

利用ES-MDA算法融合質(zhì)量濃度與水頭數(shù)據(jù)或融合ERT數(shù)據(jù)，均可以實現(xiàn)對污染源強(qiáng)度和滲透系數(shù)場的聯(lián)合反演，并且隨著迭代次數(shù)的增加，反演結(jié)果將在一定程度上更加趨近于真實值，證明了該數(shù)據(jù)同化算法在求解地下水參數(shù)反演問題上的有效性。

通過3組數(shù)值算例來比較融合不同類型觀測數(shù)據(jù)的ES-MDA算法的參數(shù)反演效果。算例1至算例3在反演污染源源強(qiáng)時，對任一污染源的ERMS值逐漸遞減，說明反演精度從優(yōu)到劣依次為算例3、算例2、算例1。算例1至算例3在反演滲透系數(shù)場時融合了多源觀測數(shù)據(jù)(質(zhì)量濃度、水頭與ERT數(shù)據(jù))的算例3對于lnK場主體形態(tài)和局部細(xì)節(jié)的表征要明顯優(yōu)于算例1和算例2；算例2對lnK場主體形態(tài)的表征優(yōu)于算例1，但對lnK場細(xì)節(jié)的刻畫程度與算例3仍存在差距：證明了融合多源觀測數(shù)據(jù)的ES-MDA算法在求解參數(shù)反演問題上的優(yōu)越性，能夠有效提升參數(shù)反演的精度。

算例2中的反演結(jié)果要明顯地優(yōu)于算例1，其主要原因是算例2中的ERT觀測數(shù)據(jù)的數(shù)量遠(yuǎn)多于算例1中的質(zhì)量濃度與水頭數(shù)據(jù)，即使精度略低也能夠為參數(shù)反演提供更多的信息。高密度電阻率法能夠快速而便捷地獲取大量具有空間連續(xù)性的ERT觀測數(shù)據(jù)，而傳統(tǒng)觀測方法只能在有限的觀測井處獲取少量不連續(xù)的觀測數(shù)據(jù)。因此，在實際場地調(diào)查中ERT方法將更具優(yōu)勢。后續(xù)研究考慮將數(shù)據(jù)同化方法與其他地球物理方法相結(jié)合，探索更為精確高效的算法。

ES-MDA算法僅適用于對符合高斯分布的模型參數(shù)的反演求解，而對于非均質(zhì)性更強(qiáng)、參數(shù)維度更高的非高斯場景，可以考慮結(jié)合機(jī)器學(xué)習(xí)、深度學(xué)習(xí)和多點地質(zhì)統(tǒng)計方法，進(jìn)一步探求算法的適用性和改進(jìn)策略。