地下水污染具有存在的隱蔽性和發(fā)現(xiàn)的滯后性特點,致使人們對于地下水污染源的特征缺乏了解和掌握.這給地下水污染修復(fù)方案的合理設(shè)計、污染責(zé)任認定和污染風(fēng)險評估都帶來了很大的困難[1-3].因此,關(guān)于地下水污染源識別的研究就顯得格外重要.
地下水污染源識別興起于20世紀80年代,發(fā)展到今天應(yīng)用于地下水污染源識別的方法包括直接方法、概率和地統(tǒng)計模擬方法、模擬-優(yōu)化方法和地球物理探測法等[4].其中,模擬-優(yōu)化方法,近些年來被廣泛應(yīng)用于地下水污染源識別[5-7].江思珉等[8]運用單純形模擬退火混合算法求解優(yōu)化模型,識別地下水污染源強度.隨后,又將模擬-優(yōu)化方法與卡爾曼濾波方法結(jié)合,進行基于污染羽形態(tài)對比的地下水污染源識別研究[9].肖傳寧等[10]應(yīng)用基于徑向基函數(shù)替代模型的模擬-優(yōu)化方法,識別地下水污染源的污染物質(zhì)泄漏量.侯澤宇等[11]應(yīng)用基于核極限學(xué)習(xí)機替代模型的模擬-優(yōu)化方法,對地下水重非水相流體(DNAPLs)污染源及含水層參數(shù)的進行同步識別.
盡管應(yīng)用模擬-優(yōu)化方法進行地下水污染源識別,取得了豐碩的研究成果.但是,應(yīng)用模擬-優(yōu)化方法進行地下水污染源識別研究,只會得到唯一的污染源識別結(jié)果(某一組參數(shù)取值影響下的污染源特征)[12-13].因此,基于模擬-優(yōu)化方法,考慮參數(shù)不確定性的污染源識別研究,實施起來尤為困難.而參數(shù)的不確定性客觀存在[14-15],會影響污染源的識別結(jié)果.
本研究提出將模擬-優(yōu)化方法、靈敏度分析方法、蒙特卡羅方法和克里格方法結(jié)合,進行考慮參數(shù)不確定性的地下水污染源識別研究.首先根據(jù)研究區(qū)的具體條件,建立地下水污染物質(zhì)運移數(shù)值模擬模型.運用靈敏度分析方法篩選出對模擬模型輸出結(jié)果影響最大的參數(shù).為減少調(diào)用模擬模型耗費的計算時間和負荷,基于克里格方法建立了模擬模型的替代模型.然后,確定決策變量、目標函數(shù)和約束條件,建立識別污染源的優(yōu)化模型.將替代模型做為等式約束連接到優(yōu)化模型中.對篩選出的參數(shù)抽樣90組,把所有參數(shù)都依次做為約束條件賦值到優(yōu)化模型中并求解優(yōu)化模型,得到所有參數(shù)影響下的污染源識別結(jié)果.最后,基于多組參數(shù)及其影響下的污染源特征,利用克里格方法建立描述參數(shù)與污染源特征之間關(guān)系的推算模型,運用推算模型計算成千上萬組參數(shù)影響下的污染源識別結(jié)果(蒙特卡羅模擬),并對識別結(jié)果進行定量的統(tǒng)計與分析.
1 研究方法
1.1 局部靈敏度分析方法
局部靈敏度分析方法,是用來篩選模型參數(shù)對模型輸出結(jié)果影響大小的方法.它的原理是利用模型輸出結(jié)果對輸入模型的參數(shù)求偏導(dǎo)數(shù),利用偏導(dǎo)數(shù)大小,來判斷輸入模型的參數(shù)對模型輸出結(jié)果影響程度的大小(式1),為了方便計算,式1可以轉(zhuǎn)換為式2:
式中:Sk是靈敏度系數(shù);xk是輸入模型的第k個參數(shù);xk是輸入模型的第k個參數(shù)的變化量;是參數(shù)變化時,模型輸出結(jié)果;是參數(shù)為xk時模型輸出結(jié)果;n是區(qū)域觀測井的數(shù)量.靈敏度系數(shù)越大,說明該參數(shù)對模型輸出結(jié)果影響越大.式2中符號單位根據(jù)具體分析參數(shù)而定.
1.2 克里格方法
克里格方法是地統(tǒng)計學(xué)的一種插值方法.近些年來,克里格方法被延伸為一種建立替代模型的方法,被應(yīng)用于多個工程領(lǐng)域[16-17].克里格方法的原理如下:
式中:qk為待定參數(shù);和分別是第i和第j個樣本的k維取值.
b基函數(shù)的待定參數(shù),可以通過最優(yōu)線性無偏估計可以求得:
1.3 蒙特卡羅方法
蒙特卡羅方法,又稱隨機抽樣或統(tǒng)計試驗方法.蒙特卡羅方法的基本思想是通過“試驗”的方法,統(tǒng)計某個隨機事件發(fā)生的頻率,或是某個隨機變量的一些數(shù)字特征,以這種事件出現(xiàn)的頻率估計這一隨機事件的概率,或者將某些數(shù)字特征,作為隨機事件的解.
蒙特卡羅方法的實現(xiàn)一般包括以下幾步:(1)基于需要解決的問題,構(gòu)造易于實現(xiàn)的概率統(tǒng)計模型或者隨機過程.(2)確定輸入模型的隨機變量和隨機變量抽樣方法.(3)基于概率統(tǒng)計模型或隨機過程,進行多次模擬試驗,對試驗的結(jié)果進行統(tǒng)計與分析,將某一結(jié)果的發(fā)生頻率或者某些數(shù)字特征(包括均值和方差、標準差等)作為問題的解.蒙特卡羅方法的基本原理,詳見尹增謙等[18]和朱輝等[19]文章.
2 案例研究
2.1 研究區(qū)概況
圖1 研究區(qū)域概況
Fig.1 Overview of the study area
S1和S2分別代表第1和第2個污染源;O1和O7分別代表第1到第7口觀測井
本文借鑒文獻[4,20]的案例進行研究.研究區(qū)是具有5個參數(shù)分區(qū),邊界不規(guī)則的二維非均質(zhì)各項同性承壓含水層,地下水流為非穩(wěn)定流,水流方向由AB邊界指向CD邊界(圖1);AB和CD邊界是已知水頭邊界,水頭分別為100 和80m;AC和BD邊界為隔水邊界;研究區(qū)在垂直方向上接受均勻的水量補給,補給率為0.0000864m/d.含水層的水文地質(zhì)參數(shù)見表1.地下水中污染物質(zhì)的初始濃度為100mg/L.污染物質(zhì)遷移的總模擬時間為10a,共有20個模擬期(每6個月為1個模擬期,每月30d).假設(shè)污染物質(zhì)是不經(jīng)過生物轉(zhuǎn)化或化學(xué)變化的保守污染物.污染源的位置已知,在第1個模擬期初至第4個模擬期末,污染物被持續(xù)釋放到地下水中,在后來的模擬期不再釋放污染物質(zhì)(表2).區(qū)域有7口觀測井.含水層中5個參數(shù)區(qū)的分布,觀測井和污染源的位置見圖1.
表1 含水層參數(shù)
Table 1 Aquifer parameters
表2 污染物質(zhì)在各釋放時段的釋放強度
Table 2 Contaminant release intensity in each release periods
基于以上的研究案例,本研究的技術(shù)路線見圖2.
2.2 數(shù)值模擬模型
根據(jù)研究區(qū)的具體條件,建立水文地質(zhì)概念模型.在概念模型的基礎(chǔ)上,建立水流和溶質(zhì)運移數(shù)值模擬模型.描述二維承壓含水層系統(tǒng)中非穩(wěn)態(tài)流的地下水水流的控制偏微分方程如下:
式中:K是滲透系數(shù),m/d;H是水頭,m;W是水流模擬模型的源匯項,m/d;ms是貯水率或釋水率,m-1;t是時間,d;x代表橫向方向,y代表縱向方向(長度單位為m).
圖2 技術(shù)路線圖
Fig.2 Technology Roadmap
圖3 水流空間分布
Fig.3 Spatial distribution map of water flow
描述溶質(zhì)運移的控制偏微分方程如下:
(12)
式中:q是孔隙度,無量綱;C是污染物質(zhì)濃度,mg/L;ux是橫向水流實際平均流速(縱向水流實際平均流速與橫向計算方法相同),m/d;D是彌散度,m2/d;R是溶質(zhì)運移模擬模型的源匯項,mg/d.
圖4 污染物質(zhì)空間分布
Fig.4 Spatial distribution map of contaminant
達西定律可用于計算式12中的ui,如式13所示:
(13)
建立水流和溶質(zhì)運移模擬模型后,使用GMS軟件中MODFLOW和MT3DMS工具箱來模擬水流和污染物質(zhì)運移的過程.
與實際問題不同,假想例子沒有實際監(jiān)測數(shù)據(jù).所以將真實的污染源特征,輸入污染物質(zhì)運移模擬模型,正向運行模擬模型得到觀測井的污染物質(zhì)濃度,作為反向識別污染源特征時的實測數(shù)據(jù).水流和污染物質(zhì)在第1800和3600d的空間分布情況分別見圖3和圖4.觀測井的污染物質(zhì)濃度見圖5.
圖5 觀測井的污染物質(zhì)濃度實測數(shù)據(jù)
Fig.5 Measured data of contaminant concentration in observation wells
2.3 靈敏度分析
參數(shù)不確定性會影響污染物質(zhì)運移模擬模型的輸出結(jié)果.考慮模型所有參數(shù)對輸出結(jié)果的影響,又會導(dǎo)致模型過于復(fù)雜.因此,應(yīng)用局部靈敏度分析方法[21],篩選出對模擬模型輸出結(jié)果影響最大的參數(shù),作為輸入模擬模型的隨機變量,其它參數(shù)作為模擬模型的確定性變量.
篩選隨機變量時,首先,將輸入模擬模型的參數(shù)取均值(表3)輸入模型,運行模型,得到觀測井的污染物質(zhì)濃度.然后,將輸入模型的某個參數(shù)加減10%,20%,其它參數(shù)保持不變.再次運行模型,得到某一參數(shù)變化后,對應(yīng)的觀測井污染物質(zhì)濃度.利用式(2)計算出各參數(shù)的靈敏度系數(shù),選擇靈敏度系數(shù)最大的參數(shù),作為輸入模型的隨機變量.靈敏度分析結(jié)果見圖6.
表3 參數(shù)服從的概率分布及取值范圍
Table 3 Probability distribution and value range of parameters
由圖6可以看出,靈敏度系數(shù)最大的參數(shù)是滲透系數(shù).因此,將滲透系數(shù)作為輸入模型的隨機變量,其它參數(shù)作為確定性變量,取定值輸入模型.根據(jù)前人經(jīng)驗[24],滲透系數(shù)服從的概率分布見表3.利用拉丁超立方抽樣方法[22-23]對滲透系數(shù)在給定范圍內(nèi)抽樣600組,每組參數(shù)樣本具有5個值,分別對應(yīng)研究區(qū)的5個參數(shù)分區(qū).同時,利用該方法對污染物質(zhì)釋放強度抽樣600組,每組樣本具8個值,分別是2個污染源在4個釋放時段的污染物質(zhì)釋放強度.使?jié)B透系數(shù)和釋放強度隨機組合,得到600組由滲透系數(shù)和釋放強度組成的輸入樣本.
圖6 參數(shù)靈敏度分析
Fig.6 Parameter sensitivity analysis graph
應(yīng)用模擬-優(yōu)化方法識別地下水污染源時,污染物質(zhì)運移模擬模型會作為等式約束連接到優(yōu)化模型中,求解優(yōu)化模型時,成百上千次的調(diào)用模擬模型,會產(chǎn)生大量的計算負荷,浪費大量的計算時間.建立模擬模型的替代模型,能有效解決這一問題.
2.4 建立替代模型
應(yīng)用克里格方法建立污染物質(zhì)運移模擬模型的替代模型.參考克里格方法的原理,利用MATLAB編寫克里格替代模型的訓(xùn)練程序.然后,建立污染物質(zhì)運移模擬模型的替代模型.建立替代模型的步驟如下:
1.將600組輸入樣本依次輸入污染物質(zhì)運移模擬模型,運行模擬模型,可以得到與600組輸入樣本一一對應(yīng)的600組觀測井的污染物質(zhì)濃度.
2.利用540組輸入樣本作為替代模型的輸入,污染物質(zhì)的濃度作為替代模型的輸出,訓(xùn)練替代模型.
3.將余下的60組輸入樣本,依次輸入替代模型,得到替代模型輸出的60組觀測井的污染物質(zhì)濃度.利用60組替代模型的輸出和模擬模型的輸出,檢驗替代模型的精度.
應(yīng)用確定性系數(shù)和平均相對誤差檢驗替代模型的精度,它們的計算方法見式14和15:
式中:yi是模擬模型輸出的第i個時段末的污染物質(zhì)濃度,mg/L;是替代模型輸出的第i個時段末的污染物質(zhì)濃度,mg/L;是模擬模型輸出的各時段末的污染物質(zhì)濃度均值,mg/L;n是觀測樣本的總數(shù)量.
2.5 建立優(yōu)化模型
建立替代模型后,需要建立識別地下水污染源特征的非線性優(yōu)化模型.非線性優(yōu)化模型包括目標函數(shù)、決策變量和約束條件3部分.將觀測井的污染物質(zhì)實測濃度與模擬計算濃度擬合(替代模型的輸出濃度)最小平方誤差和,作為目標函數(shù);將污染物質(zhì)的釋放強度,作為優(yōu)化模型的決策變量;建立識別地下水污染源特征的非線性優(yōu)化模型.優(yōu)化模型的約束條件包括污染源釋放強度不等式約束、5個分區(qū)滲透系數(shù)和地下水溶質(zhì)運移規(guī)律等式約束;優(yōu)化模型表示如下:
式中:Ki是某個參數(shù)分區(qū)的滲透系數(shù),m/d;qm是污染物質(zhì)釋放強度,mg/d;是在觀測井處污染物質(zhì)的模擬計算濃度,mg/L;是在觀測井處的污染物質(zhì)實測濃度,mg/L.是替代模型.
利用遺傳算法[25]求解某滲透系數(shù)組約束下的優(yōu)化模型,只能得到與該組滲透系數(shù)影響下的污染源識別結(jié)果.而研究滲透系數(shù)不確定性對污染源識別結(jié)果的影響,需要統(tǒng)計多組滲透系數(shù)樣本影響下的污染源識別結(jié)果.因此,應(yīng)用蒙特卡羅方法進行考慮滲透系數(shù)不確定性的污染源識別.
2.6 建立推算模型
應(yīng)用蒙特卡羅的思想,將優(yōu)化模型作為“試驗”發(fā)生器,將某滲透系數(shù)樣本影響下的污染源識別結(jié)果,作為某次隨機“試驗”的結(jié)果.通過對成千上萬組隨機“試驗”結(jié)果的統(tǒng)計與分析,考慮滲透系數(shù)不確定性對污染源識別結(jié)果的影響.
建立優(yōu)化模型后,再次應(yīng)用拉丁超立方方法,對滲透系數(shù)抽樣90組.將第1組滲透系數(shù),賦值到優(yōu)化模型中,然后求解優(yōu)化模型(整個求解過程中,滲透系數(shù)一直保持不變),得到與第1組滲透系數(shù)影響下的污染物質(zhì)釋放強度(包括S1的4個污染物質(zhì)釋放強度值和S2的4個污染物質(zhì)釋放強度值).對第2組滲透系數(shù)做與第1組滲透系數(shù)相同處理,直到第90組滲透系數(shù).這樣就依次求解了90個優(yōu)化模型,獲得了90組滲透系數(shù)影響下的污染物質(zhì)釋放強度.
但是,僅僅考慮90組滲透系數(shù)對污染源識別結(jié)果的影響,很難具有代表性.而求解成千上萬個優(yōu)化模型,又會產(chǎn)生海量的計算負荷,花費大量的計算時間.因此提出了應(yīng)用克里格方法,建立能夠描述滲透系數(shù)與污染物質(zhì)釋放強度之間關(guān)系的推算模型,來推算不同滲透系數(shù)取值影響下的污染物質(zhì)釋放強度.將70組滲透系數(shù)和與之影響下的污染物質(zhì)釋放強度,分別作為推算模型的輸入與輸出,訓(xùn)練推算模型.利用余下的20組滲透系數(shù)影響下的污染物質(zhì)釋放強度作為檢驗數(shù)據(jù)(每組包括S1的4個污染物質(zhì)釋放強度值和S2的4個污染物質(zhì)釋放強度值),檢驗推算模型的精度.推算模型與替代模型的精度評價指標相同.
3 結(jié)果與討論
3.1 替代模型的精度評價
表4 替代模型的R2和MRE
Table 4R2and MRE of surrogate model
圖7 7口觀測井的污染物質(zhì)濃度擬合
Fig.7 Fitting graph of contaminant concentration of seven observation wells
(a)~(g)分別表示第1到第7口觀測井
對基于克里格方法替代模型的精度進行了評價.通常情況下,如果替代模型的R2高于0.9,MRE在10%以內(nèi),則認為替代模型精度滿足研究需要.結(jié)合表4和圖7,圖8可以看出,替代模型與模擬模型輸出濃度的平均相對誤差均小于2%,確定性系數(shù)都大于0.99,接近1,替代模型與模擬模型的輸出結(jié)果擬合程度非常好,替代模型的精度很高,可以代替模擬模型連接到優(yōu)化模型中.
由圖8可知,利用7口觀測井濃度數(shù)據(jù)進行替代模型的精度檢驗時,出現(xiàn)了個別的異常值,但是異常值的相對誤差均未超過4%,在10%以內(nèi),可見即使對于異常值,替代模型對模擬模型的擬合程度仍然很高.出現(xiàn)異常值的可能原因是,建立替代模型的訓(xùn)練數(shù)據(jù),在個別的輸入樣本取值及取值領(lǐng)域內(nèi)覆蓋的不夠全面,導(dǎo)致替代模型在該輸入樣本取值處泛化性能欠佳.因此,在獲取訓(xùn)練數(shù)據(jù)時,盡量去覆蓋輸入樣本的取值范圍是很有必要的.
圖8 7口觀測井的污染物質(zhì)濃度相對誤差箱線圖
Fig.8 The relative error box plot of the contaminant concentration of seven observation wells
3.2 識別結(jié)果的不確定性分析
圖9 推算模型擬合曲線
Fig.9 Fitting graph of the inference model
應(yīng)用克里格方法,建立了描述滲透系數(shù)與污染物質(zhì)釋放強度之間關(guān)系的推算模型.由圖9可知, 推算模型的確定性系數(shù)達到了0.9895,與1非常接近;平均相對誤差為4.51%,在10%以內(nèi);推算模型的精度很高,可以用來推算任意滲透系數(shù)影響下的污染物質(zhì)釋放強度.
考慮滲透系數(shù)不確定性的地下水污染源識別,需要統(tǒng)計成百上千組滲透系數(shù)影響下的污染源識別結(jié)果.具體做法是:再次對滲透系數(shù)抽樣8000組,將其全部輸入推算模型,可以得到各組滲透系數(shù)影響下的污染物質(zhì)釋放強度.對8000組污染物質(zhì)釋放強度進行統(tǒng)計與分析.分析因為滲透系數(shù)不確定性,導(dǎo)致的污染源識別結(jié)果的不確定性.
由圖10可以看出,受滲透系數(shù)不確定性影響,污染源的識別結(jié)果也具有不確定性.應(yīng)用切比雪夫不等式[22]估計不同置信水平對應(yīng)的污染源識別結(jié)果置信區(qū)間.由表5可知,不同置信程度對應(yīng)的置信區(qū)間,均包含2個污染源的實際釋放強度.隨著置信程度增大,置信區(qū)間變大,反之則反.決策者可以根據(jù)不同的管理目標,選擇相信不同置信程度對應(yīng)的置信區(qū)間,做為污染源特征的識別結(jié)果.也可以選擇概率密度最大的釋放強度作為污染源識別結(jié)果(識別結(jié)果見表6).
圖10 不確定性分析
Fig.10 Uncertainty analysis diagram
(a)~(d)是S1各個釋放時段的污染物質(zhì)釋放強度;(e)~(h)是S2各個釋放時段的污染物質(zhì)釋放強度
表5 不同置信水平對應(yīng)的污染源識別結(jié)果置信區(qū)間
Table 5 Confidence interval of contamination source identification results corresponding to different confidence levels
由表6可知,識別得到的S1和S2的釋放強度與真實的污染源特征雖然有一定的誤差,但是整體上對污染源的真實特征逼近程度較好.其中,S1第4時段的釋放強度,與其真實釋放強度差距最大,導(dǎo)致這樣結(jié)果的原因可能有以下2方面:(1)替代模型和推算模型均具有誤差,會在一定程度上影響識別結(jié)果的準確性.(2)抽樣的過程具有隨機性,抽樣的樣本不一定會完美覆蓋到真實的含水層滲透系數(shù).含水層滲透系數(shù)偏離真值,會導(dǎo)致污染源的識別結(jié)果偏離污染源特征的真實結(jié)果.
求解一個優(yōu)化模型,進行1000次迭代計算,需要的時間約為0.45h,求解8000個優(yōu)化模型需要3600h.而建立推算模型需要約40.5h,應(yīng)用推算模型推算8000組滲透系數(shù)影響下的污染源特征,僅需要不到1s(瞬間輸出).應(yīng)用推算模型可以節(jié)省約99%的計算時間和計算負荷.
表6 概率密度最大的污染源識別結(jié)果
Table 6 Identification results of contamination sources with the highest probability density
4 結(jié)論
4.1 基于克里格方法建立的替代模型具有較高的精度,確定性系數(shù)高于0.99,平均相對誤差小于1.2%.將克里格替代模型做為等式約束條件連接到優(yōu)化模型中,供求解優(yōu)化模型調(diào)用,在保證一定精度的前提下,可以減少大量的負荷和時間.
4.2 應(yīng)用克里格方法,建立了描述滲透系數(shù)與污染物質(zhì)釋放強度之間關(guān)系的推算模型.建立的推算模型具有較高的精度,確定性系數(shù)和平均相對誤差分別為0.9895和4.51%.運用克里格推算模型計算了8000滲透系數(shù)影響下的污染源特征,節(jié)省了約99%的計算負荷和時間.
4.3 對8000組污染源特征進行定量的統(tǒng)計與分析,統(tǒng)計得到了置信水平為80%、60%、40%和20%對應(yīng)的污染源識別結(jié)果置信區(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,將概率密度最大的釋放強度做為污染源識別結(jié)果,對污染源的真實特征有較好的逼近程度.
4.4 將模擬-優(yōu)化方法、靈敏度分析方法、蒙特卡羅方法和克里格方法結(jié)合,進行考慮滲透系數(shù)不確定性的污染源識別研究,更加符合實際情況,而且能夠得到更加豐富的污染源識別結(jié)果,為決策者提供更多的參考依據(jù).
聲明:本文來源:土行者 作者:李久輝 盧文喜 常振波 王涵 范越
本網(wǎng)站出于傳遞更多信息而非盈利之目的,內(nèi)容僅供參考。版權(quán)歸原作者所有,若有侵權(quán),請聯(lián)系我們刪除。
山東銘翔環(huán)??萍加邢薰?/strong>
市場部電話:
0531-86989660
傳真:
0531-86982068
郵箱:
mingxiang6688@163.com
地址:
濟南市歷下區(qū)工業(yè)南路106號?