基于隨機統計方法的地下水污染源反演識別

張 宇， 呂 軍*， 盧文喜， 劉洪超

(1.水利部松遼水利委員會, 長春 130021; 2.吉林大學新能源與環境學院，長春 130021)

由于人類活動強度的增大導致地下水污染程度日趨嚴重，對飲用水安全和生態環境構成了威脅。地下水污染具有存在的隱蔽性、發現的滯后性特點，這給地下水污染修復方案設計、污染風險評估、污染責任認定都帶來了很大的困難。地下水污染源識別是制定合理修復方案的前提，也是進行地下水污染風險評估及污染責任認定的必要條件[1]。

地下水污染源反演識別是指運用有限且離散的地下水觀測數據，對地下水污染數學模擬模型進行反演求解，識別確定污染源的個數、位置和釋放歷史。地下水污染源反演識別作為一種典型的數理方程的逆問題，其求解的復雜性在于問題的不適定性[2]。目前，地下水污染源反演識別的研究尚處在發展階段，主要包括：地球化學指紋技術和數學模擬方法等。然而，單獨使用地球化學指紋技術不能完全解決確定污染源位置及追溯污染質釋放歷史的問題，且該方法工作量大、經驗性強、實施成本較為昂貴[3-4]。概括而言，地下水污染源反演識別的數學模擬方法可劃分為確定性方法和隨機統計方法兩類[5]。Skaggs等應用吉洪諾夫正則化及準可逆性方法追溯一維均質各向同性含水層介質已知單個污染源位置的污染源釋放歷史[6-7]。然而，該方法采用解析解進行研究，不能滿足實際場地污染源反演識別的需要。Mahar[8]運用模擬-優化模型的求解結果對監測井孔進行優化布設，把水質觀測數據作為反饋信息循環迭代識別污染源的釋放歷史。Bagtzoglou[9]運用隨機統計方法進行污染源反演識別研究。Gzyl等[10]提出了一種基于地質統計學的多步方法識別污染源位置及釋放歷史。與國外的研究相比，中國對于地下水污染源反演識別問題的研究尚處于起步階段。朱嵩[11]運用結合馬爾科夫鏈蒙特卡洛(MCMC)抽樣的貝葉斯方法求解理想算例水動力-水質耦合場條件下的污染源識別反問題；江思珉等[12]基于污染羽形態對比進行地下水污染源識別研究，確定了二維非均質含水層中多個污染源的位置及釋放歷史；肖傳寧等[13]對比分析了基于兩種耦合方法的模擬-優化模型在地下水污染源識別中的計算結果。

目前為止，國內有關地下水污染源識別問題的研究還僅限于運用單一一種確定性方法或隨機統計方法解決污染源反演識別中某一方面的問題。尚未見到將多種方法組合使用以綜合解決污染源反演識別問題方面的研究報道。因此，如何在現場調查的基礎上，通過科學合理的方法反演識別污染源的特征是一個亟待解決且具有重要理論和實際意義的科學問題。鑒于此，以某垃圾填埋場垃圾滲濾液地下泄漏的情形為例，對地下水污染源的反演識別問題進行研究。通過污染質運移模擬模型、替代模型、隨機統計方法(伴隨狀態方法及貝葉斯方法)的綜合運用，反演識別出地下水污染源個數、位置及釋放歷史。

1 地下水污染源反演識別的隨機統計方法

運用隨機統計方法進行的地下水污染源反演識別主要包括二個方面：第一，運用伴隨狀態方法反演識別污染源的個數及位置；第二，基于上述反演識別結果，運用貝葉斯方法反演識別污染源的釋放歷史。

1.1 伴隨狀態方法

伴隨狀態方法，又稱伴隨變分方法，是一種借助泛函分析進行靈敏度計算的方法[14]。針對梯度類優化方法中梯度向量計算量大且難以求解的問題，通過引入伴隨狀態變量，推導出梯度向量的顯式計算公式，從而可以一次性地算出靈敏度矩陣H。

在對流-彌散方程一般形式的基礎上運用伴隨變分方法推導出伴隨方程。充分考慮模擬模型初始及邊界條件的基礎上，對構造出的變分等式分部積分，隨之推導出求解伴隨方程所需的伴隨初始和邊界條件。上述過程中，推導出的伴隨方程和伴隨狀態下的定解條件共同構成了求解地下水污染源識別反問題的數學模型，從而可以直接借助數值計算方法的求解結果推斷出污染源的個數及位置。

如果定義伴隨狀態下的時間變量τ，τ=T-t(τ為逆向時間，T為水質監測末刻時間)，為了與一般情況下的對流-彌散方程及其定解條件進行區分，定義基于伴隨狀態方法變換后的對流-彌散方程及其定解條件為

(1)

式(1)中：Ω為滲流區域；t為時間，s；i、j為直角坐標系坐標x、y、z；ψ*為伴隨狀態變量；c為污染物的溶解濃度，m/L3；ui為滲流速度，L/s；Dij為水動力彌散系數，L2/s；η為含水層介質的有效孔隙度，無量綱；q0、q1分別為單位體積多孔介質污染源項和匯項的體積流量，L3/s；cs為源項污染物濃度，m/L3；c0(x,y,z)為污染物的初始濃度；Γ1為第一類邊界條件；Γ2為第二類邊界條件；Γ3為第三類邊界條件。

說明：式(1)中t為正演模型的時間變量，T為正演模型水質監測末刻的時間(實為定值，本次研究T=1 440 d，即4年)，τ為反演模型的時間變量，三個變量的量綱均為T，但是因為應用于不同的模型中，因此用3個變量加以表示。

1.2 貝葉斯方法

地下水污染源反演識別問題，從概率統計的角度可以看成是一種貝葉斯推理問題，即通過水質觀測數據不斷地更新污染源的先驗信息，從而得到反演問題的解[15]。

(2)

2 地下水污染源反演識別實例研究

以某垃圾填埋場垃圾滲濾液地下泄漏的情形為例，對地下水污染源的反演識別問題進行研究。計算目的層可概化為二維均質各向同性的潛水含水層(1 300 m×800 m)，以100 m×100 m為基本單元，將含水層均勻剖分成104個網格。含水層水流運動為穩定流，含水層厚度為30.5 m。東西邊界為線性變化的給定水頭邊界，南北邊界為隔水邊界，如圖1所示。含水層水文地質參數見表1。

研究區存在3個潛在的污染點源，分別記為S1、S2和S3(實際上S2并未排放污染物)，監測井及污染源的分布如圖1所示。假定S1、S2和S3分時段排放污染物至含水層中，導致地下水遭受污染?，F定義每個污染物濃度監測時段為一個應力期(stress period，SP)，每個應力期的時長為4個月，整個監測過程一共持續了4年，共計12個應力期(表2)。

圖1 研究區示意圖

表1 水文地質參數值

2.1 基于伴隨狀態方法的地下水污染源反演識別研究

根據上述伴隨狀態方法的求解思路和污染物濃度監測數據(表2)，反演追溯各應力期的污染羽分布圖(圖2)。伴隨狀態方法下的污染源追溯初始時間是從正演過程的污染物監測末刻開始，而并非從污染物初始排放時刻算起。所以伴隨狀態方法反演得到的第一個污染羽[圖2(a)]對應的是污染物監測末刻(τ=T-t)，即正演過程的第11個應力期。

圖2 (a)顯示的是在τ=1 440 d(反演末刻)得到的潛在源位置處的污染羽分布，從中可以確定污染源的位置。同樣可以求解出在該時刻各個網格單元內的污染物濃度，表3是104個單元格中污染物濃度較大的8個單元格坐標及污染物濃度。

表2 監測井在各應力期的污染物濃度

圖2 伴隨狀態方法反演得到的污染羽分布圖

表3 反演得到的潛在源位置處的污染物濃度(τ=1 440 d)

將伴隨狀態方法得到的污染源特征與真實情況比對，結論如下：

(1)運用伴隨狀態方法可以確定污染源的個數。在得到的污染物初始分布圖中可以明顯識別出2個主要的污染源，排除了污染源S2的干擾。實際上，污染源S2在整個監測期內從未排放過污染物。

(2)運用伴隨狀態方法還可以確定污染源的位置。根據各個網格上的污染物模擬濃度，我們可以推斷出污染源可能位于以(6,2)和(3,2)為中心的區域內。與真實污染源的位置(6,3)和(3,3)相比，相對位置誤差為7.7%，可以滿足10%的精度要求。

2.2 基于貝葉斯推理的地下水污染源反演識別研究

伴隨狀態方法可以有效識別出污染源的位置及個數。然而隨著問題復雜程度的增加，伴隨狀態方法已不適用于反演識別多個污染源分時段排放污染物的釋放歷史問題?；谪惾~斯推理的污染源識別方法通過增加源強的先驗信息，縮小了解的空間范圍。同時通過構造一個平穩分布與目標分布相同的馬爾科夫鏈來得到反演問題可能解的樣本，并根據這些樣本進行統計推斷，進而得到反演問題的解估計。

通常情況下，基于貝葉斯方法進行地下水污染源識別需要反復多次直接求解模擬模型本身，因此會帶來龐大的計算負荷和冗長的計算時間，嚴重制約地下水污染源識別的計算效率。替代模型作為模擬模型的代替模型，可有效實現模擬模型與反演過程的連接，大幅度減少反演計算的時間。本文對比運用替代模型前后污染源識別的計算時間，檢驗該方法的可行性與有效性。

2.2.1 模擬模型的初步建立

研究區地質、水文地質條件與地下水流場圖同上述情形。研究區存在3個污染源，污染源S1、S2、S3的位置如圖3所示。各污染源分時段排放污染物至含水層中，設定每個應力期的時長為4個月，污染源活躍于前4個應力期。整個監測過程從污染物開始排放的時刻起一共持續了4年，共計12個應力期。

圖3 研究區示意圖

根據水文地質概念模型，建立數學模型為

(3)

基于確定的污染源位置和污染場地附加信息初步建立溶質運移模擬模型。模型的參數設置見表1。

2.2.2 數據的準備

數據的準備為替代模型建立所用到的輸入-輸出樣本集的獲取。

假設各時段污染源源強M的先驗分布為均勻分布U[0,60]，其先驗概率密度函數為

(4)

圖4 水質監測數據折線圖

污染物濃度數據來源于5口監測井，分別為O1(7,9)、O2(6,4)、O3(5,7)、O4(3,5)、O5(2,7)，如圖4所示。圖4中縱坐標表示污染物濃度，折線圖表征O1～O5五口監測井的水質監測數據；橫坐標表示12個應力期(共4年，每個應力期為120 d)，即SP1=120 d，SP2=240 d，…，SP12=1 440 d。

2.2.3 替代模型的建立

替代模型，又稱近似模型，是指能代替原來的模擬模型近似反映其某種輸入輸出關系的模型，其本質是模擬模型的代替模型。在解決污染源反演識別問題時用簡單的替代模型代替復雜的模擬模型，不僅節省時間且減小計算負荷，能夠有效地提高反演問題求解的計算效率。

(1)訓練及檢驗樣本的選取。采用拉丁超立方抽樣方法在[0,60]內抽取120個樣本共8組并進行隨機組合，組成初始輸入樣本集。以產生訓練樣本同樣的方式，生成替代模型的檢驗樣本集，檢驗樣本集的數目為20。將訓練樣本數據輸入到上述初步建立的溶質運移數值模擬模型中，運行模型即得到5口監測井處12個監測時段的污染物濃度，作為相應的輸出響應，即形成輸出數據集。輸入數據集和輸出數據集組合成訓練樣本集，為替代模型的建立準備基礎數據。

圖5 替代模型結果與模擬模型結果擬合圖

(2)替代模型的建立。克里格(Kriging)模型由南非的Danie Krige于1951年提出，它是一種典型的通過插值方法建立的替代模型。本研究根據克里格方法的原理基于MATLAB平臺編寫程序分別建立5口監測井的克里格替代模型，并基于自適應權重粒子群優化算法進行參數的優化。

(3)替代模型的檢驗。將20組檢驗樣本輸入訓練好的克里格替代模型，獲得輸出響應(12個時段的監測井污染物濃度值)，并將其與溶質運移模擬模型的輸出響應對比，對比結果如圖5所示。結果表明：克里格替代模型對模擬模型有較高的逼近精度，其對比結果接近1:1完美線，能夠識別并取代模擬模型的輸入輸出關系。

2.2.4 MCMC算法

貝葉斯定理數學表達式簡單，但其后驗概率密度函數的求解相當困難。為了解決這一問題，本次研究采用Markov Chain Monte Carlo(MCMC)算法中的Metropolis-Hastings搜索策略。

MCMC算法參數設置如下：污染源源強歸一化后的先驗分布為x∈[0,1]上的均勻概率分布。隨機游走的建議分布為q(x*|xt)=U(xt-step,xt+step)，步長step為參數先驗范圍的0.2%，即參數的搜索步長為0.002，迭代進行105次。污染源強度的反演迭代曲線如圖6所示。

圖6 污染源強度的反演迭代曲線

由圖6可以看出，S1在經歷了大約3 500步，S2在經歷了大約2×104步迭代后(burn-in階段)，各時段的污染源源強反演結果開始進入收斂區域。表明在此情況下，污染源源強能夠被準確的識別。舍棄S1、S2的burn-in階段的結果，對其后的鏈進行統計，統計結果見表4。

圖7為2個污染源待反演4個時段的源強的后驗概率分布直方圖，直觀地表現了在獲得水質觀測數據以后，上述待反演源強的分布規律?？梢钥闯觯篁灧植疾⒉环木鶆蚍植?，貝葉斯方法更新了源強的先驗分布。

運用貝葉斯方法反演識別污染源過程的主要計算負荷來自于反復多次調用模擬模型。初步建立的模擬模型在CPU為Intel i7 2.50 GHz，內存為8 GB的計算機上運行一次平均需要6 s。如果用一般的模擬模型來解決此問題，模擬模型需要運行105次。因此，求解此問題需要33.33 h。如果將克里格替代模型引入反演問題的求解中來代替模擬模型，在訓練和檢驗克里格模型的過程中模擬模型需要運行140次，需要14 min。而獲得反演問題的解只需要5.56 h。因此，替代模型的引入大大減小了反演問題求解的計算負荷。

表4 污染源源強反演統計結果

圖7 污染源源強的后驗概率分布直方圖

3 結論及建議

通過污染質運移模擬模型、替代模型、隨機統計方法(伴隨狀態方法及貝葉斯方法)的綜合運用，反演識別地下水污染源的特征。主要得到以下結論及建議。

(1)伴隨狀態方法和貝葉斯方法的聯合運用可以反演識別出污染源的個數、位置及釋放歷史。

(2)克里格模型對模擬模型具有較高的近似精度。貝葉斯方法反演識別污染源過程中，以替代模型作為模擬模型的轉化形式，用替代模型代替模擬模型大幅度地減小了計算負荷，并保持了較高的精度。

(3)研究將計算目的層概化為均質各向同性的潛水含水層。在今后的研究中，需要針對實際的污染場地研究豐富本文的相關理論，如水質監測數據的消噪處理、缺失數據的內插與外延等相關理論與方法的研究。