摘要：為分析參數不確定性對地下水污染源識別的影響,本文通過模擬-優化方法、靈敏度分析方法、蒙特卡羅方法和克里格方法的綜合運用,建立了描述滲透系數與污染物質釋放強度之間關系的推算模型,進行了考慮參數不確定性的地下水污染源識別研究.研究結果表明,推算模型具有較高的精度,確定性系數和平均相對誤差分別為0.9895和4.51%;運用推算模型推算了8000組滲透系數影響下的污染源識別結果,節省了約99%的計算負荷和時間;對8000組污染源識別結果進行了定量的統計與分析,得到了概率密度最大的污染源識別結果和置信水平分別為80%、60%、40%和20%對應的污染源識別結果置信區間.本研究改善了應用模擬-優化方法進行地下水污染源識別時,難以考慮參數不確定性的缺點,可以為決策者提供更多的參考依據.

關鍵詞：地下水污染源；模擬-優化；不確定性；替代模型；推算模型

地下水污染具有存在的隱蔽性和發現的滯后性特點,致使人們對于地下水污染源的特征缺乏了解和掌握.這給地下水污染修復方案的合理設計、污染責任認定和污染風險評估都帶來了很大的困難[1-3].因此,關于地下水污染源識別的研究就顯得格外重要.

地下水污染源識別興起于20世紀80年代,發展到今天應用于地下水污染源識別的方法包括直接方法、概率和地統計模擬方法、模擬-優化方法和地球物理探測法等[4].其中,模擬-優化方法,近些年來被廣泛應用于地下水污染源識別[5-7].江思珉等[8]運用單純形模擬退火混合算法求解優化模型,識別地下水污染源強度.隨后,又將模擬-優化方法與卡爾曼濾波方法結合,進行基于污染羽形態對比的地下水污染源識別研究[9].肖傳寧等[10]應用基于徑向基函數替代模型的模擬-優化方法,識別地下水污染源的污染物質泄漏量.侯澤宇等[11]應用基于核極限學習機替代模型的模擬-優化方法,對地下水重非水相流體(DNAPLs)污染源及含水層參數的進行同步識別.

盡管應用模擬-優化方法進行地下水污染源識別,取得了豐碩的研究成果.但是,應用模擬-優化方法進行地下水污染源識別研究,只會得到唯一的污染源識別結果(某一組參數取值影響下的污染源特征)[12-13].因此,基于模擬-優化方法,考慮參數不確定性的污染源識別研究,實施起來尤為困難.而參數的不確定性客觀存在[14-15],會影響污染源的識別結果.

本研究提出將模擬-優化方法、靈敏度分析方法、蒙特卡羅方法和克里格方法結合,進行考慮參數不確定性的地下水污染源識別研究.首先根據研究區的具體條件,建立地下水污染物質運移數值模擬模型.運用靈敏度分析方法篩選出對模擬模型輸出結果影響最大的參數.為減少調用模擬模型耗費的計算時間和負荷,基于克里格方法建立了模擬模型的替代模型.然后,確定決策變量、目標函數和約束條件,建立識別污染源的優化模型.將替代模型做為等式約束連接到優化模型中.對篩選出的參數抽樣90組,把所有參數都依次做為約束條件賦值到優化模型中并求解優化模型,得到所有參數影響下的污染源識別結果.最后,基于多組參數及其影響下的污染源特征,利用克里格方法建立描述參數與污染源特征之間關系的推算模型,運用推算模型計算成千上萬組參數影響下的污染源識別結果(蒙特卡羅模擬),并對識別結果進行定量的統計與分析.

1 研究方法

1.1 局部靈敏度分析方法

局部靈敏度分析方法,是用來篩選模型參數對模型輸出結果影響大小的方法.它的原理是利用模型輸出結果對輸入模型的參數求偏導數,利用偏導數大小,來判斷輸入模型的參數對模型輸出結果影響程度的大小(式1),為了方便計算,式1可以轉換為式2:

式中:Sk是靈敏度系數;xk是輸入模型的第k個參數;xk是輸入模型的第k個參數的變化量;是參數變化時,模型輸出結果;是參數為xk時模型輸出結果;n是區域觀測井的數量.靈敏度系數越大,說明該參數對模型輸出結果影響越大.式2中符號單位根據具體分析參數而定.

1.2 克里格方法

克里格方法是地統計學的一種插值方法.近些年來,克里格方法被延伸為一種建立替代模型的方法,被應用于多個工程領域[16-17].克里格方法的原理如下:

式中:qk為待定參數;和分別是第i和第j個樣本的k維取值.

b基函數的待定參數,可以通過最優線性無偏估計可以求得:

1.3 蒙特卡羅方法

蒙特卡羅方法,又稱隨機抽樣或統計試驗方法.蒙特卡羅方法的基本思想是通過“試驗”的方法,統計某個隨機事件發生的頻率,或是某個隨機變量的一些數字特征,以這種事件出現的頻率估計這一隨機事件的概率,或者將某些數字特征,作為隨機事件的解.

蒙特卡羅方法的實現一般包括以下幾步:(1)基于需要解決的問題,構造易于實現的概率統計模型或者隨機過程.(2)確定輸入模型的隨機變量和隨機變量抽樣方法.(3)基于概率統計模型或隨機過程,進行多次模擬試驗,對試驗的結果進行統計與分析,將某一結果的發生頻率或者某些數字特征(包括均值和方差、標準差等)作為問題的解.蒙特卡羅方法的基本原理,詳見尹增謙等[18]和朱輝等[19]文章.

2 案例研究

2.1 研究區概況

圖1 研究區域概況

Fig.1 Overview of the study area

S1和S2分別代表第1和第2個污染源;O1和O7分別代表第1到第7口觀測井

本文借鑒文獻[4,20]的案例進行研究.研究區是具有5個參數分區,邊界不規則的二維非均質各項同性承壓含水層,地下水流為非穩定流,水流方向由AB邊界指向CD邊界(圖1);AB和CD邊界是已知水頭邊界,水頭分別為100 和80m;AC和BD邊界為隔水邊界;研究區在垂直方向上接受均勻的水量補給,補給率為0.0000864m/d.含水層的水文地質參數見表1.地下水中污染物質的初始濃度為100mg/L.污染物質遷移的總模擬時間為10a,共有20個模擬期(每6個月為1個模擬期,每月30d).假設污染物質是不經過生物轉化或化學變化的保守污染物.污染源的位置已知,在第1個模擬期初至第4個模擬期末,污染物被持續釋放到地下水中,在后來的模擬期不再釋放污染物質(表2).區域有7口觀測井.含水層中5個參數區的分布,觀測井和污染源的位置見圖1.

表1 含水層參數

Table 1 Aquifer parameters

表2 污染物質在各釋放時段的釋放強度

Table 2 Contaminant release intensity in each release periods

基于以上的研究案例,本研究的技術路線見圖2.

2.2 數值模擬模型

根據研究區的具體條件,建立水文地質概念模型.在概念模型的基礎上,建立水流和溶質運移數值模擬模型.描述二維承壓含水層系統中非穩態流的地下水水流的控制偏微分方程如下:

式中:K是滲透系數,m/d;H是水頭,m;W是水流模擬模型的源匯項,m/d;ms是貯水率或釋水率,m-1;t是時間,d;x代表橫向方向,y代表縱向方向(長度單位為m).

圖2 技術路線圖

Fig.2 Technology Roadmap

圖3 水流空間分布

Fig.3 Spatial distribution map of water flow

描述溶質運移的控制偏微分方程如下:

(12)

式中:q是孔隙度,無量綱;C是污染物質濃度,mg/L;ux是橫向水流實際平均流速(縱向水流實際平均流速與橫向計算方法相同),m/d;D是彌散度,m2/d;R是溶質運移模擬模型的源匯項,mg/d.

圖4 污染物質空間分布

Fig.4 Spatial distribution map of contaminant

達西定律可用于計算式12中的ui,如式13所示:

(13)

建立水流和溶質運移模擬模型后,使用GMS軟件中MODFLOW和MT3DMS工具箱來模擬水流和污染物質運移的過程.