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

(1.水利部松遼水利委員會(huì), 長(zhǎng)春 130021; 2.吉林大學(xué)新能源與環(huán)境學(xué)院，長(zhǎng)春 130021)

由于人類活動(dòng)強(qiáng)度的增大導(dǎo)致地下水污染程度日趨嚴(yán)重，對(duì)飲用水安全和生態(tài)環(huán)境構(gòu)成了威脅。地下水污染具有存在的隱蔽性、發(fā)現(xiàn)的滯后性特點(diǎn)，這給地下水污染修復(fù)方案設(shè)計(jì)、污染風(fēng)險(xiǎn)評(píng)估、污染責(zé)任認(rèn)定都帶來(lái)了很大的困難。地下水污染源識(shí)別是制定合理修復(fù)方案的前提，也是進(jìn)行地下水污染風(fēng)險(xiǎn)評(píng)估及污染責(zé)任認(rèn)定的必要條件[1]。

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

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

1 地下水污染源反演識(shí)別的隨機(jī)統(tǒng)計(jì)方法

運(yùn)用隨機(jī)統(tǒng)計(jì)方法進(jìn)行的地下水污染源反演識(shí)別主要包括二個(gè)方面：第一，運(yùn)用伴隨狀態(tài)方法反演識(shí)別污染源的個(gè)數(shù)及位置；第二，基于上述反演識(shí)別結(jié)果，運(yùn)用貝葉斯方法反演識(shí)別污染源的釋放歷史。

1.1 伴隨狀態(tài)方法

伴隨狀態(tài)方法，又稱伴隨變分方法，是一種借助泛函分析進(jìn)行靈敏度計(jì)算的方法[14]。針對(duì)梯度類優(yōu)化方法中梯度向量計(jì)算量大且難以求解的問(wèn)題，通過(guò)引入伴隨狀態(tài)變量，推導(dǎo)出梯度向量的顯式計(jì)算公式，從而可以一次性地算出靈敏度矩陣H。

在對(duì)流-彌散方程一般形式的基礎(chǔ)上運(yùn)用伴隨變分方法推導(dǎo)出伴隨方程。充分考慮模擬模型初始及邊界條件的基礎(chǔ)上，對(duì)構(gòu)造出的變分等式分部積分，隨之推導(dǎo)出求解伴隨方程所需的伴隨初始和邊界條件。上述過(guò)程中，推導(dǎo)出的伴隨方程和伴隨狀態(tài)下的定解條件共同構(gòu)成了求解地下水污染源識(shí)別反問(wèn)題的數(shù)學(xué)模型，從而可以直接借助數(shù)值計(jì)算方法的求解結(jié)果推斷出污染源的個(gè)數(shù)及位置。

如果定義伴隨狀態(tài)下的時(shí)間變量τ，τ=T-t(τ為逆向時(shí)間，T為水質(zhì)監(jiān)測(cè)末刻時(shí)間)，為了與一般情況下的對(duì)流-彌散方程及其定解條件進(jìn)行區(qū)分，定義基于伴隨狀態(tài)方法變換后的對(duì)流-彌散方程及其定解條件為

(1)

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

說(shuō)明：式(1)中t為正演模型的時(shí)間變量，T為正演模型水質(zhì)監(jiān)測(cè)末刻的時(shí)間(實(shí)為定值，本次研究T=1 440 d，即4年)，τ為反演模型的時(shí)間變量，三個(gè)變量的量綱均為T，但是因?yàn)閼?yīng)用于不同的模型中，因此用3個(gè)變量加以表示。

1.2 貝葉斯方法

地下水污染源反演識(shí)別問(wèn)題，從概率統(tǒng)計(jì)的角度可以看成是一種貝葉斯推理問(wèn)題，即通過(guò)水質(zhì)觀測(cè)數(shù)據(jù)不斷地更新污染源的先驗(yàn)信息，從而得到反演問(wèn)題的解[15]。

(2)

2 地下水污染源反演識(shí)別實(shí)例研究

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

研究區(qū)存在3個(gè)潛在的污染點(diǎn)源，分別記為S1、S2和S3(實(shí)際上S2并未排放污染物)，監(jiān)測(cè)井及污染源的分布如圖1所示。假定S1、S2和S3分時(shí)段排放污染物至含水層中，導(dǎo)致地下水遭受污染?，F(xiàn)定義每個(gè)污染物濃度監(jiān)測(cè)時(shí)段為一個(gè)應(yīng)力期(stress period，SP)，每個(gè)應(yīng)力期的時(shí)長(zhǎng)為4個(gè)月，整個(gè)監(jiān)測(cè)過(guò)程一共持續(xù)了4年，共計(jì)12個(gè)應(yīng)力期(表2)。

圖1 研究區(qū)示意圖

表1 水文地質(zhì)參數(shù)值

2.1 基于伴隨狀態(tài)方法的地下水污染源反演識(shí)別研究

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

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

表2 監(jiān)測(cè)井在各應(yīng)力期的污染物濃度

圖2 伴隨狀態(tài)方法反演得到的污染羽分布圖

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

將伴隨狀態(tài)方法得到的污染源特征與真實(shí)情況比對(duì)，結(jié)論如下：

(1)運(yùn)用伴隨狀態(tài)方法可以確定污染源的個(gè)數(shù)。在得到的污染物初始分布圖中可以明顯識(shí)別出2個(gè)主要的污染源，排除了污染源S2的干擾。實(shí)際上，污染源S2在整個(gè)監(jiān)測(cè)期內(nèi)從未排放過(guò)污染物。

(2)運(yùn)用伴隨狀態(tài)方法還可以確定污染源的位置。根據(jù)各個(gè)網(wǎng)格上的污染物模擬濃度，我們可以推斷出污染源可能位于以(6,2)和(3,2)為中心的區(qū)域內(nèi)。與真實(shí)污染源的位置(6,3)和(3,3)相比，相對(duì)位置誤差為7.7%，可以滿足10%的精度要求。

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

伴隨狀態(tài)方法可以有效識(shí)別出污染源的位置及個(gè)數(shù)。然而隨著問(wèn)題復(fù)雜程度的增加，伴隨狀態(tài)方法已不適用于反演識(shí)別多個(gè)污染源分時(shí)段排放污染物的釋放歷史問(wèn)題?；谪惾~斯推理的污染源識(shí)別方法通過(guò)增加源強(qiáng)的先驗(yàn)信息，縮小了解的空間范圍。同時(shí)通過(guò)構(gòu)造一個(gè)平穩(wěn)分布與目標(biāo)分布相同的馬爾科夫鏈來(lái)得到反演問(wèn)題可能解的樣本，并根據(jù)這些樣本進(jìn)行統(tǒng)計(jì)推斷，進(jìn)而得到反演問(wèn)題的解估計(jì)。

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

2.2.1 模擬模型的初步建立

研究區(qū)地質(zhì)、水文地質(zhì)條件與地下水流場(chǎng)圖同上述情形。研究區(qū)存在3個(gè)污染源，污染源S1、S2、S3的位置如圖3所示。各污染源分時(shí)段排放污染物至含水層中，設(shè)定每個(gè)應(yīng)力期的時(shí)長(zhǎng)為4個(gè)月，污染源活躍于前4個(gè)應(yīng)力期。整個(gè)監(jiān)測(cè)過(guò)程從污染物開(kāi)始排放的時(shí)刻起一共持續(xù)了4年，共計(jì)12個(gè)應(yīng)力期。

圖3 研究區(qū)示意圖

根據(jù)水文地質(zhì)概念模型，建立數(shù)學(xué)模型為

(3)

基于確定的污染源位置和污染場(chǎng)地附加信息初步建立溶質(zhì)運(yùn)移模擬模型。模型的參數(shù)設(shè)置見(jiàn)表1。

2.2.2 數(shù)據(jù)的準(zhǔn)備

數(shù)據(jù)的準(zhǔn)備為替代模型建立所用到的輸入-輸出樣本集的獲取。

假設(shè)各時(shí)段污染源源強(qiáng)M的先驗(yàn)分布為均勻分布U[0,60]，其先驗(yàn)概率密度函數(shù)為

(4)

圖4 水質(zhì)監(jiān)測(cè)數(shù)據(jù)折線圖

污染物濃度數(shù)據(jù)來(lái)源于5口監(jiān)測(cè)井，分別為O1(7,9)、O2(6,4)、O3(5,7)、O4(3,5)、O5(2,7)，如圖4所示。圖4中縱坐標(biāo)表示污染物濃度，折線圖表征O1～O5五口監(jiān)測(cè)井的水質(zhì)監(jiān)測(cè)數(shù)據(jù)；橫坐標(biāo)表示12個(gè)應(yīng)力期(共4年，每個(gè)應(yīng)力期為120 d)，即SP1=120 d，SP2=240 d，…，SP12=1 440 d。

2.2.3 替代模型的建立

替代模型，又稱近似模型，是指能代替原來(lái)的模擬模型近似反映其某種輸入輸出關(guān)系的模型，其本質(zhì)是模擬模型的代替模型。在解決污染源反演識(shí)別問(wèn)題時(shí)用簡(jiǎn)單的替代模型代替復(fù)雜的模擬模型，不僅節(jié)省時(shí)間且減小計(jì)算負(fù)荷，能夠有效地提高反演問(wèn)題求解的計(jì)算效率。

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

圖5 替代模型結(jié)果與模擬模型結(jié)果擬合圖

(2)替代模型的建立?？死锔?Kriging)模型由南非的Danie Krige于1951年提出，它是一種典型的通過(guò)插值方法建立的替代模型。本研究根據(jù)克里格方法的原理基于MATLAB平臺(tái)編寫(xiě)程序分別建立5口監(jiān)測(cè)井的克里格替代模型，并基于自適應(yīng)權(quán)重粒子群優(yōu)化算法進(jìn)行參數(shù)的優(yōu)化。

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

2.2.4 MCMC算法

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

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

圖6 污染源強(qiáng)度的反演迭代曲線

由圖6可以看出，S1在經(jīng)歷了大約3 500步，S2在經(jīng)歷了大約2×104步迭代后(burn-in階段)，各時(shí)段的污染源源強(qiáng)反演結(jié)果開(kāi)始進(jìn)入收斂區(qū)域。表明在此情況下，污染源源強(qiáng)能夠被準(zhǔn)確的識(shí)別。舍棄S1、S2的burn-in階段的結(jié)果，對(duì)其后的鏈進(jìn)行統(tǒng)計(jì)，統(tǒng)計(jì)結(jié)果見(jiàn)表4。

圖7為2個(gè)污染源待反演4個(gè)時(shí)段的源強(qiáng)的后驗(yàn)概率分布直方圖，直觀地表現(xiàn)了在獲得水質(zhì)觀測(cè)數(shù)據(jù)以后，上述待反演源強(qiáng)的分布規(guī)律?？梢钥闯?，后驗(yàn)分布并不服從均勻分布，貝葉斯方法更新了源強(qiáng)的先驗(yàn)分布。

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

表4 污染源源強(qiáng)反演統(tǒng)計(jì)結(jié)果

圖7 污染源源強(qiáng)的后驗(yàn)概率分布直方圖

3 結(jié)論及建議

通過(guò)污染質(zhì)運(yùn)移模擬模型、替代模型、隨機(jī)統(tǒng)計(jì)方法(伴隨狀態(tài)方法及貝葉斯方法)的綜合運(yùn)用，反演識(shí)別地下水污染源的特征。主要得到以下結(jié)論及建議。

(1)伴隨狀態(tài)方法和貝葉斯方法的聯(lián)合運(yùn)用可以反演識(shí)別出污染源的個(gè)數(shù)、位置及釋放歷史。

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

(3)研究將計(jì)算目的層概化為均質(zhì)各向同性的潛水含水層。在今后的研究中，需要針對(duì)實(shí)際的污染場(chǎng)地研究豐富本文的相關(guān)理論，如水質(zhì)監(jiān)測(cè)數(shù)據(jù)的消噪處理、缺失數(shù)據(jù)的內(nèi)插與外延等相關(guān)理論與方法的研究。