李久輝,盧文喜*,常振波,王 涵,范 越

考慮參數(shù)不確定性的地下水污染源識(shí)別

李久輝1,2,3,盧文喜1,2,3*,常振波1,2,3,王 涵1,2,3,范 越1,2,3

(1.吉林大學(xué)地下水資源與環(huán)境教育部重點(diǎn)實(shí)驗(yàn)室,吉林 長(zhǎng)春 130012；2.吉林大學(xué)吉林省水資源與水環(huán)境重點(diǎn)實(shí)驗(yàn)室,吉林 長(zhǎng)春 130012；3.吉林大學(xué)新能源與環(huán)境學(xué)院,吉林 長(zhǎng)春 130012)

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

地下水污染源；模擬-優(yōu)化；不確定性；替代模型；推算模型

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

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

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

1 研究方法

1.1 局部靈敏度分析方法

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

(2)

1.2 克里格方法

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

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

式中:為個(gè)樣本點(diǎn)之間的相關(guān)矩陣,其表達(dá)式如下:

式5中的可以通過(guò)求解無(wú)約束優(yōu)化規(guī)劃計(jì)算得到,無(wú)約束優(yōu)化規(guī)劃目標(biāo)函數(shù)見(jiàn)式9.在得到后,和可以相繼被計(jì)算出來(lái).應(yīng)用克里格方法建立的替代模型,在點(diǎn)處的響應(yīng)值()的估計(jì)值為:

1.3 蒙特卡羅方法

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

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

2 案例研究

2.1 研究區(qū)概況

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

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

本文借鑒文獻(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ù)見(jiàn)表1.地下水中污染物質(zhì)的初始濃度為100mg/L.污染物質(zhì)遷移的總模擬時(shí)間為10a,共有20個(gè)模擬期(每6個(gè)月為1個(gè)模擬期,每月30d).假設(shè)污染物質(zhì)是不經(jīng)過(guò)生物轉(zhuǎn)化或化學(xué)變化的保守污染物.污染源的位置已知,在第1個(gè)模擬期初至第4個(gè)模擬期末,污染物被持續(xù)釋放到地下水中,在后來(lái)的模擬期不再釋放污染物質(zhì)(表2).區(qū)域有7口觀測(cè)井.含水層中5個(gè)參數(shù)區(qū)的分布,觀測(cè)井和污染源的位置見(jiàn)圖1.

表1 含水層參數(shù)

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

基于以上的研究案例,本研究的技術(shù)路線見(jiàn)圖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代表縱向方向(長(zhǎng)度單位為m).

圖3 水流空間分布

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

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

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

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

與實(shí)際問(wèn)題不同,假想例子沒(méi)有實(shí)際監(jiān)測(cè)數(shù)據(jù).所以將真實(shí)的污染源特征,輸入污染物質(zhì)運(yùn)移模擬模型,正向運(yùn)行模擬模型得到觀測(cè)井的污染物質(zhì)濃度,作為反向識(shí)別污染源特征時(shí)的實(shí)測(cè)數(shù)據(jù).水流和污染物質(zhì)在第1800和3600d的空間分布情況分別見(jiàn)圖3和圖4.觀測(cè)井的污染物質(zhì)濃度見(jiàn)圖5.

圖5 觀測(cè)井的污染物質(zhì)濃度實(shí)測(cè)數(shù)據(jù)

2.3 靈敏度分析

參數(shù)不確定性會(huì)影響污染物質(zhì)運(yùn)移模擬模型的輸出結(jié)果.考慮模型所有參數(shù)對(duì)輸出結(jié)果的影響,又會(huì)導(dǎo)致模型過(guò)于復(fù)雜.因此,應(yīng)用局部靈敏度分析方法[21],篩選出對(duì)模擬模型輸出結(jié)果影響最大的參數(shù),作為輸入模擬模型的隨機(jī)變量,其它參數(shù)作為模擬模型的確定性變量.

篩選隨機(jī)變量時(shí),首先,將輸入模擬模型的參數(shù)取均值(表3)輸入模型,運(yùn)行模型,得到觀測(cè)井的污染物質(zhì)濃度.然后,將輸入模型的某個(gè)參數(shù)加減10%,20%,其它參數(shù)保持不變.再次運(yùn)行模型,得到某一參數(shù)變化后,對(duì)應(yīng)的觀測(cè)井污染物質(zhì)濃度.利用式(2)計(jì)算出各參數(shù)的靈敏度系數(shù),選擇靈敏度系數(shù)最大的參數(shù),作為輸入模型的隨機(jī)變量.靈敏度分析結(jié)果見(jiàn)圖6.

表3 參數(shù)服從的概率分布及取值范圍

由圖6可以看出,靈敏度系數(shù)最大的參數(shù)是滲透系數(shù).因此,將滲透系數(shù)作為輸入模型的隨機(jī)變量,其它參數(shù)作為確定性變量,取定值輸入模型.根據(jù)前人經(jīng)驗(yàn)[24],滲透系數(shù)服從的概率分布見(jiàn)表3.利用拉丁超立方抽樣方法[22-23]對(duì)滲透系數(shù)在給定范圍內(nèi)抽樣600組,每組參數(shù)樣本具有5個(gè)值,分別對(duì)應(yīng)研究區(qū)的5個(gè)參數(shù)分區(qū).同時(shí),利用該方法對(duì)污染物質(zhì)釋放強(qiáng)度抽樣600組,每組樣本具8個(gè)值,分別是2個(gè)污染源在4個(gè)釋放時(shí)段的污染物質(zhì)釋放強(qiáng)度.使?jié)B透系數(shù)和釋放強(qiáng)度隨機(jī)組合,得到600組由滲透系數(shù)和釋放強(qiáng)度組成的輸入樣本.

圖6 參數(shù)靈敏度分析

應(yīng)用模擬-優(yōu)化方法識(shí)別地下水污染源時(shí),污染物質(zhì)運(yùn)移模擬模型會(huì)作為等式約束連接到優(yōu)化模型中,求解優(yōu)化模型時(shí),成百上千次的調(diào)用模擬模型,會(huì)產(chǎn)生大量的計(jì)算負(fù)荷,浪費(fèi)大量的計(jì)算時(shí)間.建立模擬模型的替代模型,能有效解決這一問(wèn)題.

2.4 建立替代模型

應(yīng)用克里格方法建立污染物質(zhì)運(yùn)移模擬模型的替代模型.參考克里格方法的原理,利用MATLAB編寫克里格替代模型的訓(xùn)練程序.然后,建立污染物質(zhì)運(yùn)移模擬模型的替代模型.建立替代模型的步驟如下:

1.將600組輸入樣本依次輸入污染物質(zhì)運(yùn)移模擬模型,運(yùn)行模擬模型,可以得到與600組輸入樣本一一對(duì)應(yīng)的600組觀測(cè)井的污染物質(zhì)濃度.

2.利用540組輸入樣本作為替代模型的輸入,污染物質(zhì)的濃度作為替代模型的輸出,訓(xùn)練替代模型.

3.將余下的60組輸入樣本,依次輸入替代模型,得到替代模型輸出的60組觀測(cè)井的污染物質(zhì)濃度.利用60組替代模型的輸出和模擬模型的輸出,檢驗(yàn)替代模型的精度.

應(yīng)用確定性系數(shù)和平均相對(duì)誤差檢驗(yàn)替代模型的精度,它們的計(jì)算方法見(jiàn)式14和15:

2.5 建立優(yōu)化模型

建立替代模型后,需要建立識(shí)別地下水污染源特征的非線性優(yōu)化模型.非線性優(yōu)化模型包括目標(biāo)函數(shù)、決策變量和約束條件3部分.將觀測(cè)井的污染物質(zhì)實(shí)測(cè)濃度與模擬計(jì)算濃度擬合(替代模型的輸出濃度)最小平方誤差和,作為目標(biāo)函數(shù);將污染物質(zhì)的釋放強(qiáng)度,作為優(yōu)化模型的決策變量;建立識(shí)別地下水污染源特征的非線性優(yōu)化模型.優(yōu)化模型的約束條件包括污染源釋放強(qiáng)度不等式約束、5個(gè)分區(qū)滲透系數(shù)和地下水溶質(zhì)運(yùn)移規(guī)律等式約束;優(yōu)化模型表示如下:

利用遺傳算法[25]求解某滲透系數(shù)組約束下的優(yōu)化模型,只能得到與該組滲透系數(shù)影響下的污染源識(shí)別結(jié)果.而研究滲透系數(shù)不確定性對(duì)污染源識(shí)別結(jié)果的影響,需要統(tǒng)計(jì)多組滲透系數(shù)樣本影響下的污染源識(shí)別結(jié)果.因此,應(yīng)用蒙特卡羅方法進(jìn)行考慮滲透系數(shù)不確定性的污染源識(shí)別.

2.6 建立推算模型

應(yīng)用蒙特卡羅的思想,將優(yōu)化模型作為“試驗(yàn)”發(fā)生器,將某滲透系數(shù)樣本影響下的污染源識(shí)別結(jié)果,作為某次隨機(jī)“試驗(yàn)”的結(jié)果.通過(guò)對(duì)成千上萬(wàn)組隨機(jī)“試驗(yàn)”結(jié)果的統(tǒng)計(jì)與分析,考慮滲透系數(shù)不確定性對(duì)污染源識(shí)別結(jié)果的影響.

建立優(yōu)化模型后,再次應(yīng)用拉丁超立方方法,對(duì)滲透系數(shù)抽樣90組.將第1組滲透系數(shù),賦值到優(yōu)化模型中,然后求解優(yōu)化模型(整個(gè)求解過(guò)程中,滲透系數(shù)一直保持不變),得到與第1組滲透系數(shù)影響下的污染物質(zhì)釋放強(qiáng)度(包括S1的4個(gè)污染物質(zhì)釋放強(qiáng)度值和S2的4個(gè)污染物質(zhì)釋放強(qiáng)度值).對(duì)第2組滲透系數(shù)做與第1組滲透系數(shù)相同處理,直到第90組滲透系數(shù).這樣就依次求解了90個(gè)優(yōu)化模型,獲得了90組滲透系數(shù)影響下的污染物質(zhì)釋放強(qiáng)度.

但是,僅僅考慮90組滲透系數(shù)對(duì)污染源識(shí)別結(jié)果的影響,很難具有代表性.而求解成千上萬(wàn)個(gè)優(yōu)化模型,又會(huì)產(chǎn)生海量的計(jì)算負(fù)荷,花費(fèi)大量的計(jì)算時(shí)間.因此提出了應(yīng)用克里格方法,建立能夠描述滲透系數(shù)與污染物質(zhì)釋放強(qiáng)度之間關(guān)系的推算模型,來(lái)推算不同滲透系數(shù)取值影響下的污染物質(zhì)釋放強(qiáng)度.將70組滲透系數(shù)和與之影響下的污染物質(zhì)釋放強(qiáng)度,分別作為推算模型的輸入與輸出,訓(xùn)練推算模型.利用余下的20組滲透系數(shù)影響下的污染物質(zhì)釋放強(qiáng)度作為檢驗(yàn)數(shù)據(jù)(每組包括S1的4個(gè)污染物質(zhì)釋放強(qiáng)度值和S2的4個(gè)污染物質(zhì)釋放強(qiáng)度值),檢驗(yàn)推算模型的精度.推算模型與替代模型的精度評(píng)價(jià)指標(biāo)相同.

3 結(jié)果與討論

3.1 替代模型的精度評(píng)價(jià)

表4 替代模型的R2和MRE

圖7 7口觀測(cè)井的污染物質(zhì)濃度擬合

(a)~(g)分別表示第1到第7口觀測(cè)井

對(duì)基于克里格方法替代模型的精度進(jìn)行了評(píng)價(jià).通常情況下,如果替代模型的2高于0.9,MRE在10%以內(nèi),則認(rèn)為替代模型精度滿足研究需要.結(jié)合表4和圖7,圖8可以看出,替代模型與模擬模型輸出濃度的平均相對(duì)誤差均小于2%,確定性系數(shù)都大于0.99,接近1,替代模型與模擬模型的輸出結(jié)果擬合程度非常好,替代模型的精度很高,可以代替模擬模型連接到優(yōu)化模型中.

由圖8可知,利用7口觀測(cè)井濃度數(shù)據(jù)進(jìn)行替代模型的精度檢驗(yàn)時(shí),出現(xiàn)了個(gè)別的異常值,但是異常值的相對(duì)誤差均未超過(guò)4%,在10%以內(nèi),可見(jiàn)即使對(duì)于異常值,替代模型對(duì)模擬模型的擬合程度仍然很高.出現(xiàn)異常值的可能原因是,建立替代模型的訓(xùn)練數(shù)據(jù),在個(gè)別的輸入樣本取值及取值領(lǐng)域內(nèi)覆蓋的不夠全面,導(dǎo)致替代模型在該輸入樣本取值處泛化性能欠佳.因此,在獲取訓(xùn)練數(shù)據(jù)時(shí),盡量去覆蓋輸入樣本的取值范圍是很有必要的.

圖8 7口觀測(cè)井的污染物質(zhì)濃度相對(duì)誤差箱線圖

3.2 識(shí)別結(jié)果的不確定性分析

圖9 推算模型擬合曲線

應(yīng)用克里格方法,建立了描述滲透系數(shù)與污染物質(zhì)釋放強(qiáng)度之間關(guān)系的推算模型.由圖9可知, 推算模型的確定性系數(shù)達(dá)到了0.9895,與1非常接近;平均相對(duì)誤差為4.51%,在10%以內(nèi);推算模型的精度很高,可以用來(lái)推算任意滲透系數(shù)影響下的污染物質(zhì)釋放強(qiáng)度.

考慮滲透系數(shù)不確定性的地下水污染源識(shí)別,需要統(tǒng)計(jì)成百上千組滲透系數(shù)影響下的污染源識(shí)別結(jié)果.具體做法是:再次對(duì)滲透系數(shù)抽樣8000組,將其全部輸入推算模型,可以得到各組滲透系數(shù)影響下的污染物質(zhì)釋放強(qiáng)度.對(duì)8000組污染物質(zhì)釋放強(qiáng)度進(jìn)行統(tǒng)計(jì)與分析.分析因?yàn)闈B透系數(shù)不確定性,導(dǎo)致的污染源識(shí)別結(jié)果的不確定性.

由圖10可以看出,受滲透系數(shù)不確定性影響,污染源的識(shí)別結(jié)果也具有不確定性.應(yīng)用切比雪夫不等式[22]估計(jì)不同置信水平對(duì)應(yīng)的污染源識(shí)別結(jié)果置信區(qū)間.由表5可知,不同置信程度對(duì)應(yīng)的置信區(qū)間,均包含2個(gè)污染源的實(shí)際釋放強(qiáng)度.隨著置信程度增大,置信區(qū)間變大,反之則反.決策者可以根據(jù)不同的管理目標(biāo),選擇相信不同置信程度對(duì)應(yīng)的置信區(qū)間,做為污染源特征的識(shí)別結(jié)果.也可以選擇概率密度最大的釋放強(qiáng)度作為污染源識(shí)別結(jié)果(識(shí)別結(jié)果見(jiàn)表6).

(a)~(d)是S1各個(gè)釋放時(shí)段的污染物質(zhì)釋放強(qiáng)度;(e)~(h)是S2各個(gè)釋放時(shí)段的污染物質(zhì)釋放強(qiáng)度

表5 不同置信水平對(duì)應(yīng)的污染源識(shí)別結(jié)果置信區(qū)間

由表6可知,識(shí)別得到的S1和S2的釋放強(qiáng)度與真實(shí)的污染源特征雖然有一定的誤差,但是整體上對(duì)污染源的真實(shí)特征逼近程度較好.其中,S1第4時(shí)段的釋放強(qiáng)度,與其真實(shí)釋放強(qiáng)度差距最大,導(dǎo)致這樣結(jié)果的原因可能有以下2方面:(1)替代模型和推算模型均具有誤差,會(huì)在一定程度上影響識(shí)別結(jié)果的準(zhǔn)確性.(2)抽樣的過(guò)程具有隨機(jī)性,抽樣的樣本不一定會(huì)完美覆蓋到真實(shí)的含水層滲透系數(shù).含水層滲透系數(shù)偏離真值,會(huì)導(dǎo)致污染源的識(shí)別結(jié)果偏離污染源特征的真實(shí)結(jié)果.

求解一個(gè)優(yōu)化模型,進(jìn)行1000次迭代計(jì)算,需要的時(shí)間約為0.45h,求解8000個(gè)優(yōu)化模型需要3600h.而建立推算模型需要約40.5h,應(yīng)用推算模型推算8000組滲透系數(shù)影響下的污染源特征,僅需要不到1s(瞬間輸出).應(yīng)用推算模型可以節(jié)省約99%的計(jì)算時(shí)間和計(jì)算負(fù)荷.

表6 概率密度最大的污染源識(shí)別結(jié)果

4 結(jié)論

4.1 基于克里格方法建立的替代模型具有較高的精度,確定性系數(shù)高于0.99,平均相對(duì)誤差小于1.2%.將克里格替代模型做為等式約束條件連接到優(yōu)化模型中,供求解優(yōu)化模型調(diào)用,在保證一定精度的前提下,可以減少大量的負(fù)荷和時(shí)間.

4.2 應(yīng)用克里格方法,建立了描述滲透系數(shù)與污染物質(zhì)釋放強(qiáng)度之間關(guān)系的推算模型.建立的推算模型具有較高的精度,確定性系數(shù)和平均相對(duì)誤差分別為0.9895和4.51%.運(yùn)用克里格推算模型計(jì)算了8000滲透系數(shù)影響下的污染源特征,節(jié)省了約99%的計(jì)算負(fù)荷和時(shí)間.

4.3 對(duì)8000組污染源特征進(jìn)行定量的統(tǒng)計(jì)與分析,統(tǒng)計(jì)得到了置信水平為80%、60%、40%和20%對(duì)應(yīng)的污染源識(shí)別結(jié)果置信區(qū)間;分析出S1在4個(gè)釋放時(shí)段概率密度最大的釋放強(qiáng)度分別是74.22×105、53.67×105、40.35×105和16.36 ×105mg/d,S2在4個(gè)釋放時(shí)段概率密度最大的釋放強(qiáng)度分別是55.08× 105、50.50×105、39.41×105和21.01×105mg/d,將概率密度最大的釋放強(qiáng)度做為污染源識(shí)別結(jié)果,對(duì)污染源的真實(shí)特征有較好的逼近程度.

4.4 將模擬-優(yōu)化方法、靈敏度分析方法、蒙特卡羅方法和克里格方法結(jié)合,進(jìn)行考慮滲透系數(shù)不確定性的污染源識(shí)別研究,更加符合實(shí)際情況,而且能夠得到更加豐富的污染源識(shí)別結(jié)果,為決策者提供更多的參考依據(jù).

[1] 王景瑞,胡立堂.地下水污染源識(shí)別的數(shù)學(xué)方法研究進(jìn)展[J]. 水科學(xué)進(jìn)展, 2017,28(6):943-952. Wang J R, Hu L T. Advances in mathematical methods of groundwater pollution source identification [J]. Advances in Water Science, 2017, 28(6):943-952.

[2] Snodgrass, M F, Kitanidis P K. A geostatistical approach to contaminant source identification [J]. Water Resources Research, 1997,33(4):537-546.

[3] Lapworth D J, Baran N, Stuart M E, et al. Emerging organic contaminants in groundwater: A review of sources, fate and occurrence [J]. Environmental Pollution, 2012,163(4):287-303

[4] Xing Z X, Qu R, Zhao Y, et al. Identifying the release history of a groundwater contaminant source based on an ensemble surrogate model [J]. Journal of Hydrology, 2019,572:501-516.

[5] Gorelick S M, Evans B, Remson I. Identifying sources of groundwater pollution: An optimization approach [J]. Water Resources Research, 1983,19(3):779-790.

[6] Mirghani B Y, Mahinthakumar K G, Tryby M E. A parallel evolutionary strategy based simulation–optimization approach for solving groundwater source identification problems [J]. Advances in Water Resources, 2009,32(9):1373–1385.

[7] Datta B, Chakrabarty D, Dhar A. Identification of unknown groundwater pollution sources using classical optimization with linked simulation [J]. Journal of Hydro-environment Research, 2011,5(1): 25–36.

[8] 江思珉,張亞力,蔡 奕,等.單純形模擬退火算法反演地下水污染源強(qiáng)度[J]. 同濟(jì)大學(xué)學(xué)報(bào)(自然科學(xué)版),2013,41(2):253-257. Jiang S M, Zhang Y L, Cai Y, et al. Groundwater contamination identification by hybrid simplex method of simulated annealing [J]. Journal of Tongji university (Natural science), 2013,41(2):253-257.

[9] 江思珉,張亞力,周念清,等.基于污染羽形態(tài)對(duì)比的地下水污染源識(shí)別研究[J]. 水利學(xué)報(bào), 2014,45(6):735-74. Jiang S M, Zhang Y L, Zhou N Q, et al.Groundwater contaminant source identification based on the morphological comparison of pollution plume [J] Journal of Hydraulic Engineering, 2014,45(6): 735-74.

[10] 肖傳寧,盧文喜,趙 瑩,等.基于徑向基函數(shù)模型的優(yōu)化方法在地下水污染源識(shí)別中的應(yīng)用[J]. 中國(guó)環(huán)境科學(xué), 2016,36(7):2067-2072.Xiao C N, Lu W X, Zhao Y, et al. Optimization method of identification of groundwater pollution sources based on radial basis function model [J]. China Environmental Science, 2016,36(7):2067- 2072.

[11] 侯澤宇,盧文喜,王 宇.基于替代模型的地下水DNAPLs污染源反演識(shí)別[J]. 中國(guó)環(huán)境科學(xué), 2019,39(1):188-195. Hou Z Y, Lu W X, Wang Y. Surrogate-based source identification of DNAPLs-contaminated groundwater [J]. China Environmental Science, 2019,39(1):188-195.

[12] Zhao Y, Lu W X, An Y K. Surrogate model-based simulation- optimization approach for groundwater source identification problems [J]. Environmental Forensics, 2015,16(3):296-303

[13] 肖傳寧,盧文喜,安永凱,等.基于兩種耦合方法的模擬-優(yōu)化模型在地下水污染源識(shí)別中的對(duì)比[J]. 中國(guó)環(huán)境科學(xué), 2015,35(8): 2393-2399.Xiao C N, Lu W X, An Y K, et al. Simulation-optimization model of identification of groundwater pollution sources based on two coupling method [J]. China Environmental Science,2015,35(8):2393-2399.

[14] 束龍倉(cāng),陶玉飛,劉佩貴.考慮水文地質(zhì)參數(shù)不確定性的地下水補(bǔ)給量可靠度計(jì)算[J]. 水利學(xué)報(bào), 2008,(3):346-350. Shu L C, Tao Y F, Liu P G. Reliability calculation method for groundwater recharge in consideration of uncertainty of hydrogeological parameters [J]. Journal of Hydraulic Engineering, 2008,(3):346-350.

[15] 吳雯倩,靳孟貴.淮北市地下水流數(shù)值模擬及水文地質(zhì)參數(shù)不確定性分析[J]. 水文地質(zhì)工程地質(zhì), 2014,41(3):21-28.Wu W Q, Jin M G. An analysis of period features of groundwater micro-dynamics in the North China Plain [J].Hydrogeology & Engineering Geology,2014,41(3):21-28.

[16] 吳義忠,陳立平.多領(lǐng)域物理系統(tǒng)的仿真優(yōu)化方法[M]. 科學(xué)出版社, 2011. Wu Y Z, Cheng L P.Simulation optimization method of multi domain physical system [M].Science Press, 2011.

[17] 李久輝,盧文喜,常振波,等.基于不確定性分析的地下水污染超標(biāo)風(fēng)險(xiǎn)預(yù)警[J]. 中國(guó)環(huán)境科學(xué), 2017,37(6):2270-2277. Li J H, Lu W X, Chang Z B, et al. Risk prediction of groundwater pollution based on uncertainty analysis [J]. China Environmental Science, 2017,37(6):2270-2277.

[18] 尹增謙,管景峰,張曉宏,等.蒙特卡羅方法及應(yīng)用[J]. 物理與工程, 2002,(3):45-49. Yin Z Q, Guan J F, Zhang X H, et al. The Monte Carlo method and its application [J].Physics and Engineering, 2002,(3):45-49.

[19] 朱 輝,劉義保,游運(yùn).蒙特卡羅方法與擬蒙特卡羅方法的歷史、現(xiàn)狀及展望[J]. 東華理工大學(xué)學(xué)報(bào)(自然科學(xué)版), 2010,33(4):357- 362. Zhu H, Liu Y B, You Y. Monte Carlo method and quasi Monte Carlo method [J]. Journal of East China University of Technology (Natural Science), 2010,33(4):357-362.

[20] Ayvaz M T. A linked simulation–optimization model for solving the unknown groundwater pollution source identification problems [J]. Journal of Contaminant Hydrology, 2010,117(1–4):46–59.

[21] 常振波,盧文喜,辛 欣,等.基于靈敏度分析和替代模型的地下水污染風(fēng)險(xiǎn)評(píng)價(jià)方法[J]. 中國(guó)環(huán)境科學(xué), 2017,37(1):167-173. Chang Z B, Lu W X, Xin X, et al. Risk prediction of groundwater pollution based on uncertainty analysis [J]. China Environmental Science, 2017,37(1):167-173.

[22] Parnianifard A, Azfanizam A S, Ariffin M K A, et al. Comparative study of metamodeling and sampling design for expensive and semi- expensive simulation models under uncertainty [J]. Simulation- Transactions of the Society for Modeling and Simulation Interational, 2020,96(1):89-110.

[23] 施小清,吳吉春,姜蓓蕾,等.基于LHS方法的地下水流模型不確定性分析[J]. 水文地質(zhì)工程地質(zhì), 2009,36(2):1-6. Shi X Q, Wu J C, Jiang B L, et al. Uncertainty analysis of groundwater models based on the Latin Hypercube sampling technique [J].Hydrogeology & Engineering Geology,2009,36(2):1-6.

[24] 苗添升,盧文喜,歐陽(yáng)琦,等.地下水?dāng)?shù)值模擬的不確定性分析在水質(zhì)預(yù)測(cè)中的應(yīng)用[J]. 水電能源科學(xué), 2016,34(8):20-23,44.Miao T S, Lu W X, Ouyang Q, et al. Application of uncertainty analysis of groundwater numerical simulation in water quality prediction [J].Water Resources and Power,2016,34(8):20-23,44.

[25] Guo J Y, Lu W X, Yang Q C, et al. The application of 0~1mixed integer nonlinear programming optimization model based on a surrogate model to identify the groundwater pollution source [J]. Journal of Contaminant Hydrology, 2018,220:18-25.

Identification of groundwater contamination sources considering parameter uncertainty.

LI Jiu-hui1,2,3, LU Wen-xi1,2,3*, CHANG Zhen-bo1,2,3, WANG Han1,2,3, FAN Yue1,2,3

(1.Key Laboratory of Groundwater Resources and Environment, Jilin University, Ministry of Education, Changchun 130021, China；2.Jilin Provincial Key Laboratory of Water Resources and Environment, Jilin University, Changchun 130021, China；3.College of New Energy and Environment, Jilin University, Changchun 130021, China)., 2021,41(4)：1711~1722

To analyze the influence of parameter uncertainty on the identification of groundwater contamination sources, a calculation model describing the relationship between hydraulic conductivity and contaminant release intensity was established through the comprehensive application of the simulation-based optimization method, sensitivity analysis method, Monte Carlo method and Kriging method. The results showed that the calculation model had high accuracy, the certainty coefficient and average relative error were 0.9895 and 4.51%, respectively; the contamination sources identification results under the influence of 8000 groups of hydraulic conductivity were calculated with the new model, which brought about 99% saving in the calculation load and time. Through quantitative statistical analysis of 8000 groups of contamination sources identification results, the contamination sources with the highest probability density were identified, along with the confidence intervals of contamination sources identification results corresponding to 80%, 60%, 40% and 20%. This study improved the applicability of simulation-based optimization method to identify groundwater contamination sources when facing the disadvantage caused by the parameter uncertainty, thus, can provide more reference for decision makers.

groundwater contamination sources；simulation-optimization；uncertainty；surrogate model；calculation model

X523

A

1000-6923(2021)04-1711-12

李久輝(1992-),女,內(nèi)蒙古赤峰人, 吉林大學(xué)博士研究生,主要從事地下水?dāng)?shù)值模擬方面的研究.發(fā)表論文18余篇.

2020-09-08

國(guó)家自然科學(xué)基金資助項(xiàng)目(41672232);吉林省科技發(fā)展項(xiàng)目(20170101066JC)

* 責(zé)任作者, 教授, luwenxi@jlu.edu.cn