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

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

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

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

(1.吉林大學地下水資源與環(huán)境教育部重點實驗室,吉林 長春 130012；2.吉林大學吉林省水資源與水環(huán)境重點實驗室,吉林 長春 130012；3.吉林大學新能源與環(huán)境學院,吉林 長春 130012)

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

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

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

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

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

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

1 研究方法

1.1 局部靈敏度分析方法

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

(2)

1.2 克里格方法

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

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

式中:為個樣本點之間的相關矩陣,其表達式如下:

式5中的可以通過求解無約束優(yōu)化規(guī)劃計算得到,無約束優(yōu)化規(guī)劃目標函數(shù)見式9.在得到后,和可以相繼被計算出來.應用克里格方法建立的替代模型,在點處的響應值()的估計值為:

1.3 蒙特卡羅方法

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

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

2 案例研究

2.1 研究區(qū)概況

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

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

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

表1 含水層參數(shù)

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

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

2.2 數(shù)值模擬模型

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

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

圖3 水流空間分布

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

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

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

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

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

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

2.3 靈敏度分析

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

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

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

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

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

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

2.4 建立替代模型

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

1.將600組輸入樣本依次輸入污染物質運移模擬模型,運行模擬模型,可以得到與600組輸入樣本一一對應的600組觀測井的污染物質濃度.

2.利用540組輸入樣本作為替代模型的輸入,污染物質的濃度作為替代模型的輸出,訓練替代模型.

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

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

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

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

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

2.6 建立推算模型

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

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

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

3 結果與討論

3.1 替代模型的精度評價

表4 替代模型的R2和MRE

圖7 7口觀測井的污染物質濃度擬合

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

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

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

圖8 7口觀測井的污染物質濃度相對誤差箱線圖

3.2 識別結果的不確定性分析

圖9 推算模型擬合曲線

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

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

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

(a)~(d)是S1各個釋放時段的污染物質釋放強度;(e)~(h)是S2各個釋放時段的污染物質釋放強度

表5 不同置信水平對應的污染源識別結果置信區(qū)間

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

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

表6 概率密度最大的污染源識別結果

4 結論

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

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

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

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

[1] 王景瑞,胡立堂.地下水污染源識別的數(shù)學方法研究進展[J]. 水科學進展, 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] 江思珉,張亞力,蔡 奕,等.單純形模擬退火算法反演地下水污染源強度[J]. 同濟大學學報(自然科學版),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)對比的地下水污染源識別研究[J]. 水利學報, 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)化方法在地下水污染源識別中的應用[J]. 中國環(huán)境科學, 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污染源反演識別[J]. 中國環(huán)境科學, 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)化模型在地下水污染源識別中的對比[J]. 中國環(huán)境科學, 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] 束龍倉,陶玉飛,劉佩貴.考慮水文地質參數(shù)不確定性的地下水補給量可靠度計算[J]. 水利學報, 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ù)值模擬及水文地質參數(shù)不確定性分析[J]. 水文地質工程地質, 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] 吳義忠,陳立平.多領域物理系統(tǒng)的仿真優(yōu)化方法[M]. 科學出版社, 2011. Wu Y Z, Cheng L P.Simulation optimization method of multi domain physical system [M].Science Press, 2011.

[17] 李久輝,盧文喜,常振波,等.基于不確定性分析的地下水污染超標風險預警[J]. 中國環(huán)境科學, 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] 尹增謙,管景峰,張曉宏,等.蒙特卡羅方法及應用[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] 朱 輝,劉義保,游運.蒙特卡羅方法與擬蒙特卡羅方法的歷史、現(xiàn)狀及展望[J]. 東華理工大學學報(自然科學版), 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] 常振波,盧文喜,辛 欣,等.基于靈敏度分析和替代模型的地下水污染風險評價方法[J]. 中國環(huán)境科學, 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]. 水文地質工程地質, 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] 苗添升,盧文喜,歐陽琦,等.地下水數(shù)值模擬的不確定性分析在水質預測中的應用[J]. 水電能源科學, 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-),女,內蒙古赤峰人, 吉林大學博士研究生,主要從事地下水數(shù)值模擬方面的研究.發(fā)表論文18余篇.

2020-09-08

國家自然科學基金資助項目(41672232);吉林省科技發(fā)展項目(20170101066JC)

* 責任作者, 教授, luwenxi@jlu.edu.cn