ES-MDA算法融合ERT數據聯合反演地下水污染源與含水層參數

周念清，張瑞城，江思珉，夏學敏

(1.同濟大學土木工程學院，上海 200092；2.上海理工大學環境與建筑學院，上海 200093)

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

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

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

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

1 研究方法

1.1 地下水污染遷移模型

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

(1)

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

在查明了地下水流動信息的基礎上，再結合MT3DMS程序，對地下水污染物的溶質運移過程進行模擬，其基本控制方程為

(2)

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

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

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

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

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

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

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

(3)

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

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

j=1,2,…,Ne

(4)

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

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

1.3 耦合地下水-地球物理模型的數據同化方法

根據巖石地球物理關系構建聯結地球物理信息與水文地質信息的數值模型，利用Archie公式[24]描述地層電阻率、溶液電阻率與含水飽和度之間關系，其公式基本形式如下：

(5)

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

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

(6)

式中：C為溶液的質量濃度，mg/L；T為溶液溫度，文中假定恒為25 ℃。將式(5)與式(6)聯合，便可構建地層電阻率與溶液質量濃度的數值關系。

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

(7)

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

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

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

2 數值算例

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

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

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

表1 算例設置Tab.1 The setting of case studies

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

(8)

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

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

(9)

3 結果與討論

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

3.1 污染源源強的反演識別

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

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

圖3 污染源源強的反演結果Fig.3 The inversion results of contaminant source strength

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

3.2 滲透系數場的反演識別

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

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

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

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

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

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

4 結論與展望

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

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

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

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

ES-MDA算法僅適用于對符合高斯分布的模型參數的反演求解，而對于非均質性更強、參數維度更高的非高斯場景，可以考慮結合機器學習、深度學習和多點地質統計方法，進一步探求算法的適用性和改進策略。