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

關(guān)鍵詞：地下水污染源；模擬-優(yōu)化；不確定性；替代模型；推算模型

地下水污染具有存在的隱蔽性和發(fā)現(xiàn)的滯后性特點(diǎn),致使人們對于地下水污染源的特征缺乏了解和掌握.這給地下水污染修復(fù)方案的合理設(shè)計(jì)、污染責(zé)任認(rèn)定和污染風(fēng)險(xiǎn)評估都帶來了很大的困難[1-3].因此,關(guān)于地下水污染源識(shí)別的研究就顯得格外重要.

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

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

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

1 研究方法

1.1 局部靈敏度分析方法

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

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

1.2 克里格方法

克里格方法是地統(tǒng)計(jì)學(xué)的一種插值方法.近些年來,克里格方法被延伸為一種建立替代模型的方法,被應(yīng)用于多個(gè)工程領(lǐng)域[16-17].克里格方法的原理如下:

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

b基函數(shù)的待定參數(shù),可以通過最優(yōu)線性無偏估計(jì)可以求得:

1.3 蒙特卡羅方法

蒙特卡羅方法,又稱隨機(jī)抽樣或統(tǒng)計(jì)試驗(yàn)方法.蒙特卡羅方法的基本思想是通過“試驗(yàn)”的方法,統(tǒng)計(jì)某個(gè)隨機(jī)事件發(fā)生的頻率,或是某個(gè)隨機(jī)變量的一些數(shù)字特征,以這種事件出現(xiàn)的頻率估計(jì)這一隨機(jī)事件的概率,或者將某些數(shù)字特征,作為隨機(jī)事件的解.

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

2 案例研究

2.1 研究區(qū)概況

圖1 研究區(qū)域概況

Fig.1 Overview of the study area

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

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

表1 含水層參數(shù)

Table 1 Aquifer parameters

表2 污染物質(zhì)在各釋放時(shí)段的釋放強(qiáng)度

Table 2 Contaminant release intensity in each release periods

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

2.2 數(shù)值模擬模型

根據(jù)研究區(qū)的具體條件,建立水文地質(zhì)概念模型.在概念模型的基礎(chǔ)上,建立水流和溶質(zhì)運(yùn)移數(shù)值模擬模型.描述二維承壓含水層系統(tǒng)中非穩(wěn)態(tài)流的地下水水流的控制偏微分方程如下:

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

圖2 技術(shù)路線圖

Fig.2 Technology Roadmap

圖3 水流空間分布

Fig.3 Spatial distribution map of water flow

描述溶質(zhì)運(yùn)移的控制偏微分方程如下:

(12)

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

圖4 污染物質(zhì)空間分布

Fig.4 Spatial distribution map of contaminant

達(dá)西定律可用于計(jì)算式12中的ui,如式13所示:

(13)

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