基于靈敏度分析和替代模型的地下水污染風(fēng)險(xiǎn)評價(jià)方法

常振波,盧文喜*,辛 欣,顧文龍,崔尚進(jìn)(1.吉林大學(xué),地下水與資源環(huán)境教育部重點(diǎn)實(shí)驗(yàn)室,吉林長春 130012；2.吉林大學(xué)環(huán)境與資源學(xué)院,吉林 長春 130012)

基于靈敏度分析和替代模型的地下水污染風(fēng)險(xiǎn)評價(jià)方法

常振波1,2,盧文喜1,2*,辛 欣1,2,顧文龍1,2,崔尚進(jìn)1,2(1.吉林大學(xué),地下水與資源環(huán)境教育部重點(diǎn)實(shí)驗(yàn)室,吉林長春 130012；2.吉林大學(xué)環(huán)境與資源學(xué)院,吉林 長春 130012)

采用蒙特卡洛方法,借助隨機(jī)模型進(jìn)行地下水污染風(fēng)險(xiǎn)評價(jià),模型中隨機(jī)變量利用靈敏度分析的方法確定,使地下水風(fēng)險(xiǎn)評價(jià)結(jié)果更為可靠,并借助一個(gè)假想例子來說明評價(jià)過程.結(jié)果表明,模擬輸出數(shù)據(jù)符合正態(tài)分布規(guī)律,對正態(tài)分布概率密度函數(shù)積分可以得到污染風(fēng)險(xiǎn),井1、井2和井3的污染風(fēng)險(xiǎn)分別為0%、78.52%和100%;根據(jù)整個(gè)模擬區(qū)的污染風(fēng)險(xiǎn)分布圖可以劃分出具有不同污染風(fēng)險(xiǎn)程度的子區(qū)域,藉此能夠定量評價(jià)模擬區(qū)不同子區(qū)域的污染風(fēng)險(xiǎn)程度.

地下水污染風(fēng)險(xiǎn)評價(jià)；靈敏度分析；蒙特卡洛；替代模型；概率統(tǒng)計(jì)

地下水人為污染事件時(shí)有發(fā)生[1],相比地表水污染,地下水污染修復(fù)所需費(fèi)用更高.地下水污染具有發(fā)生的隱蔽性和發(fā)現(xiàn)的滯后性[2-3],污染一經(jīng)發(fā)現(xiàn)難以修復(fù).為了緩解地下水質(zhì)惡化的現(xiàn)狀,可靠的地下水污染風(fēng)險(xiǎn)評價(jià)工作是十分必要的.決策者可以根據(jù)風(fēng)險(xiǎn)評價(jià)結(jié)果,通過人為規(guī)劃和調(diào)控降低污染發(fā)生的可能性[4-5].

風(fēng)險(xiǎn)評價(jià)實(shí)質(zhì)就是不確定性分析,沒有不確定性,就沒有風(fēng)險(xiǎn)[6].因此,對地下水溶質(zhì)運(yùn)移中不確定性因素進(jìn)行定量分析,并在結(jié)果中體現(xiàn)風(fēng)險(xiǎn)程度,可以使地下水污染風(fēng)險(xiǎn)評價(jià)結(jié)果更加科學(xué).對此國內(nèi)外展開了大量的研究,Goodrich等[7]應(yīng)用蒙特卡洛模擬,研究由于參數(shù)不確定性導(dǎo)致地下水流和溶質(zhì)運(yùn)移過程的不確定性;Bennett等[8]將蒙特卡洛方法得到的一系列污染物遷移結(jié)果用于暴露評價(jià);梁婕[9]利用貝葉斯方法推斷出滲透系數(shù)隨機(jī)場的相關(guān)參數(shù)的后驗(yàn)分布,基于參數(shù)的后驗(yàn)分布研究地下水溶質(zhì)運(yùn)移規(guī)律;申升[10]進(jìn)一步考慮了縱向彌散度對地下水溶質(zhì)運(yùn)移的影響.

前人研究中隨機(jī)變量多為事先確定,然而對于不同的模型參數(shù)對結(jié)果的影響是不同的,引入靈敏度分析方法,篩選對結(jié)果影響較大的參數(shù)作為隨機(jī)變量,利用蒙特卡洛方法對地下水溶質(zhì)運(yùn)移進(jìn)行不確定性分析,并以濃度值超過某一限定值的概率作為污染風(fēng)險(xiǎn),對整個(gè)模擬區(qū)的污染風(fēng)險(xiǎn)進(jìn)行定量評價(jià).為了減少大量重復(fù)調(diào)用模擬模型而產(chǎn)生的計(jì)算負(fù)荷,本文借助替代模型完成蒙特卡洛過程.

1 研究方法

1.1 蒙特卡洛法

蒙特卡洛法作為最常用的不確定性分析方法之一被廣泛地應(yīng)用于地下水?dāng)?shù)值模擬中[11-13].蒙特卡洛方法就是通過設(shè)置一系列平行的確定性實(shí)驗(yàn)來模擬隨機(jī)問題.在設(shè)置平行實(shí)驗(yàn)時(shí)選擇待求因素作為隨機(jī)變量,研究隨機(jī)變量變化對結(jié)果的影響.蒙特卡洛法將參數(shù)的不確定性轉(zhuǎn)化為溶質(zhì)運(yùn)移的不確定性,能夠很好地表征污染風(fēng)險(xiǎn)[14-15].本文以蒙特卡洛方法作為框架,具體研究思路如圖1所示.

1.2 靈敏度分析

靈敏度分析可以定量地評價(jià)參數(shù)不確定性對結(jié)果造成的影響,本文利用靈敏度分析的方法分析參數(shù)對輸出濃度的影響,挑選兩個(gè)靈敏度系數(shù)較大的參數(shù)作為隨機(jī)變量,進(jìn)行地下水污染風(fēng)險(xiǎn)評價(jià),這樣會使得風(fēng)險(xiǎn)評價(jià)更為可靠.分析方法可分為全局靈敏度分析和局部靈敏度分析,其中局部靈敏度分析能得到單一參數(shù)變化對結(jié)果的影響.靈敏度利用結(jié)果對參數(shù)求偏導(dǎo)數(shù),反映參數(shù)變化對結(jié)果的影響程度[16-17],即:

式中:Sk表示當(dāng)參數(shù)αk變化時(shí)對輸出結(jié)果y的影響程度,也就是靈敏度系數(shù).

對于某一特定參數(shù)靈敏度系數(shù)的求解,可通過以下方法近似,使該參數(shù)的值由αk變化為αk+Δαk,結(jié)果則由 yi( αk)變化為 yi( αk+Δαk),可用以下公式求得近似值,即:

1.3 替代模型

少量的數(shù)據(jù)無法支撐地下水污染風(fēng)險(xiǎn)評價(jià),而大量數(shù)據(jù)的獲取需要大量的時(shí)間,其中時(shí)間主要花費(fèi)在模擬模型的求解上.為了減少模擬模型的使用次數(shù),引入替代模型.替代模型能夠在某些行為上逼近模擬模型,且調(diào)用簡便,運(yùn)行時(shí)間短.蒙特卡羅模擬與替代模型結(jié)合能有效減小計(jì)算負(fù)荷[18-19].

克里格方法利用協(xié)方差的變化來表達(dá)空間的變化,現(xiàn)被延伸為一種替代模型的建立方法,廣泛地應(yīng)用于多個(gè)行業(yè),是一種黑箱模型[20].

克里格模型的方程可以表示成以下形式:

上述對 y( x)的估計(jì)值y?( x)可以分為兩部分,其中前一部分為線性回歸部分,后一部分為隨機(jī)部分.其中 f( x) =[f1( x) , f2(x) ,… , fk(x )]為已知回歸模型的基函數(shù);待定參數(shù) β = [β1, β2,… ,βk],可以通過訓(xùn)練數(shù)據(jù)求得;z( x)為隨機(jī)部分,其方差為σ2,均值為0,協(xié)方差:

R( xi, xj)為采樣點(diǎn)xi和點(diǎn)xj的關(guān)聯(lián)函數(shù),有多種類型可供選擇,本文采用較為常用的高斯模型:

式中:θk為待定參數(shù),通過求解優(yōu)化模型求得;xki為第i個(gè)樣本的k維坐標(biāo).

根據(jù)克里格模型,在預(yù)測點(diǎn) x處的響應(yīng)值y( x)的預(yù)測估計(jì)值為:

式中:f為基函數(shù),為方便起見,本文選擇常數(shù)型,f為一常數(shù)列向量;r( x)為點(diǎn)x與n個(gè)訓(xùn)練樣本采樣點(diǎn) (x1, x2,… ,xn)之間的相關(guān)向量,;y為與n個(gè)采樣點(diǎn)對應(yīng)的響應(yīng)值,為 n×1的向量;β為線性回歸部分的待定參數(shù),可以通過最優(yōu)線性無偏估計(jì)求得:

R為n個(gè)采樣點(diǎn)的相關(guān)矩陣:

方差σ2的估計(jì)值為:

待定參數(shù) θk可以通過一個(gè)無約束優(yōu)化問題求得:

通過MATLAB軟件實(shí)現(xiàn)上述替代模型的建立過程.

2 假想算例

2.1 概況

圖1 技術(shù)路線流程Fig.1 The flow chart of technical route

模擬區(qū)內(nèi)存在二維均質(zhì)各向同性潛水含水層,含水層信息如圖1所示,含水層厚為12m,模擬區(qū)初始水位為8m,水力梯度約為0.007,年平均降水量730mm,降雨入滲系數(shù)為0.2,選擇滲透系數(shù)、給水度、縱向彌散度和孔隙度4個(gè)參數(shù)進(jìn)行靈敏度分析,根據(jù)經(jīng)驗(yàn)給出其取值范圍和概率分布,如表1所示,其中因橫向彌散度與縱向彌散度的比值為 0.1,故只考慮縱向彌散度.區(qū)內(nèi)共設(shè)三口抽水井,定流量抽水,各井抽水量均為300m3/d,持續(xù)時(shí)間 1a,同時(shí)將抽水井作為水質(zhì)觀測井.現(xiàn)假設(shè)在模擬區(qū)內(nèi)建造一座工廠,工廠每天向含水層排放400m3的污水,污染物濃度500mg/L,污染質(zhì)為保守物質(zhì),遷移過程中不發(fā)生物理化學(xué)反應(yīng).含水層污染物的初始濃度為零,定水頭邊界可以視為零濃度邊界,隔水邊界可以視為零通量邊界.

表1 參數(shù)的概率分布及取值情況Table 1 The probability distribution and values information of parameters

2.2 隨機(jī)模型建立

2.2.1 數(shù)學(xué)模型 地下水溶質(zhì)運(yùn)移模型是在地下水水流模型的基礎(chǔ)上建立的,研究區(qū)地下水為二維均質(zhì)各向同性潛水非穩(wěn)定流,可以建立數(shù)學(xué)模型:

地下水溶質(zhì)運(yùn)移數(shù)學(xué)模型:

地下水水流模型和溶質(zhì)運(yùn)移模型通過達(dá)西定律聯(lián)系起來:

式中:K為含水層的滲透系數(shù),m/d;H為潛水水位, m;B為潛水含水層底板高程,m;P 為人工開采強(qiáng)度,m/d;R為降水入滲補(bǔ)給量,m/d;μ為含水層給水度,無量綱;G為模擬區(qū)范圍;S1、 S2為含水層邊界,在水流模型中為定水頭邊界,在溶質(zhì)運(yùn)移模型中為定濃度邊界;S3、 S4為含水層邊界,在水流模型中為零流量邊界,在溶質(zhì)運(yùn)移模型中為零彌散通量邊界;Kn為邊界法向量上的滲透系數(shù),m;n為含水層介質(zhì)的孔隙度,無量綱;c為污染質(zhì)濃度,mg/L;Dx、Dy為水動力彌散系數(shù)在x、y方向的分量,m2/d;vx、 vy為滲透流速v在x、y方向上的分量,m/d.

圖2 研究算例平面示意Fig.2 Plan view of the aquifer configuration

數(shù)學(xué)模型是借助GMS軟件的MODFLOW模塊和MT3DMS模塊進(jìn)行求解的.

2.2.2 隨機(jī)變量確定 本文利用靈敏度分析的方法確定隨機(jī)變量,參與分析的參數(shù)有滲透系數(shù)K、給水度μ、縱向彌散度 αL和孔隙度n.由圖3可知,4個(gè)參數(shù)的靈敏度系數(shù)由大到小依次為SK、 Sα、Sn、 Sμ,滲透系數(shù)和縱向彌散度的靈敏度系數(shù)較大,選擇這兩個(gè)參數(shù)作為隨機(jī)變量,按照其經(jīng)驗(yàn)概率分布進(jìn)行取值,剩下的兩個(gè)參數(shù)作為確定值輸入模型.

2.3 替代模型的建立

利用拉丁超立方抽樣的方法對隨機(jī)變量滲透系數(shù)和縱向彌散度按照其概率分布進(jìn)行抽樣,得到50組輸入數(shù)據(jù)集,將50組參數(shù)組合分別輸入模擬模型,利用GMS求解,得到50組輸出結(jié)果.將參數(shù)組合和輸出結(jié)果中三口觀測井的末時(shí)刻濃度值分別作為輸入輸出,利用40組輸入和輸出建立克里格替代模型,并利用 10組數(shù)據(jù)(每一組數(shù)據(jù)均包括3口觀測井的濃度值)驗(yàn)證替代模型的可靠性.圖4為模擬模型輸出結(jié)果和替代模型輸出結(jié)果對比圖,可以看出二者的相對誤差在2%以下,可認(rèn)為替代模型滿足誤差要求,并能夠替代模擬模型進(jìn)行調(diào)用.

圖3 參數(shù)靈敏度曲線Fig.3 Sensitivity curves of parameters

2.4 蒙特卡羅模擬

利用蒙特卡羅方法模擬500次試驗(yàn),對結(jié)果進(jìn)行統(tǒng)計(jì).即首先根據(jù)滲透系數(shù)和縱向彌散度的概率分布,按照拉丁超立方抽樣的方法,得到500組輸入數(shù)據(jù)集;然后將500組參數(shù)組合輸入替代模型,得到 500組輸出數(shù)據(jù);對輸出數(shù)據(jù)進(jìn)行分析.

3 結(jié)果分析

3.1 單井污染風(fēng)險(xiǎn)分析

利用SPSS軟件中的單樣本K-S檢驗(yàn)分析500組輸出,假設(shè)三口井的濃度值為正態(tài)分布,可以得到結(jié)果如表 2所示,其中當(dāng) P值大于 0.05時(shí)認(rèn)為假設(shè)成立,得到三口井的分布均為正態(tài)分布.由此畫出三口井濃度值的累積頻率直方圖和其對應(yīng)的正態(tài)分布概率密度函數(shù),如圖5(a)～(c)所示.

圖4 模型預(yù)測結(jié)果對比Fig.4 Comparison of different models’ result

表2 輸出數(shù)據(jù)單樣本K-S分析結(jié)果Table 2 The results of K-S single sample analysis on output data

圖5 觀測井污染物濃度值直方圖及正態(tài)分布概率密度函數(shù)Fig.5 Observation wells pollutant concentration histogram and normal probability density function

分析三口供水井的污染風(fēng)險(xiǎn),以污染物濃度值超過100mg/L的地下水為是被污染的,將超過100mg/L概率作為污染風(fēng)險(xiǎn)大小,利用如下公式可計(jì)算三口井污染風(fēng)險(xiǎn):

式中:P為單井污染風(fēng)險(xiǎn)大小;()f x為單井濃度值的概率密度函數(shù);σ為正態(tài)分布的標(biāo)準(zhǔn)差;μ為概率分布的均值.井1、井2和井3的污染風(fēng)險(xiǎn)分別為0%、78.52%和100%.由此可以得到,井2和井3水質(zhì)均有較大的污染風(fēng)險(xiǎn),為此可能要對工廠排放污水的濃度做出一定的限制或者尋找新的供水井,下面通過對區(qū)域的污染風(fēng)險(xiǎn)進(jìn)行評價(jià)來為尋找新的供水井提出一些建議.

3.2 區(qū)域污染風(fēng)險(xiǎn)分析

圖6 模擬區(qū)污染風(fēng)險(xiǎn)分布Fig.6 the pollution risk distribution map of the whole simulation area

為更加全面地分析整個(gè)研究區(qū)的污染風(fēng)險(xiǎn),計(jì)算出每個(gè)離散網(wǎng)格所對應(yīng)位置的污染風(fēng)險(xiǎn),并繪制等值線圖見圖 6.可以看出污染風(fēng)險(xiǎn)在不同空間位置上具有差異性,研究區(qū)上方區(qū)域污染風(fēng)險(xiǎn)較小且范圍較小,下方區(qū)域污染風(fēng)險(xiǎn)較大且范圍較大.井1位于白色區(qū)域即污染低風(fēng)險(xiǎn)區(qū)域,在此區(qū)域的地下水受到污染的風(fēng)險(xiǎn)較小,低于20%;井2位于淺色區(qū)域即污染中等風(fēng)險(xiǎn)區(qū)域,此區(qū)域地下水受到污染的風(fēng)險(xiǎn)在20%到80%;井3位于深色區(qū)域即污染高風(fēng)險(xiǎn)區(qū)域,此區(qū)域地下水受到污染的風(fēng)險(xiǎn)在80%以上.井2和井3均有較高的污染風(fēng)險(xiǎn),為了使地下水水質(zhì)能夠維持使用標(biāo)準(zhǔn),需要尋找新的供水井.圖中白色區(qū)域風(fēng)險(xiǎn)小,且靠近定水頭邊界水量充足,可以在此區(qū)域內(nèi)選擇新供水井.

4 結(jié)論

4.1 在文中假設(shè)情境下靈敏度分析方法篩選出滲透系數(shù)和縱向彌散度作為隨機(jī)參數(shù),利用這兩個(gè)參數(shù)作為隨機(jī)變量得到的地下水污染風(fēng)險(xiǎn)評價(jià)結(jié)果更加科學(xué)和可靠.

4.2 利用本文的方法進(jìn)行地下水風(fēng)險(xiǎn)評價(jià),能夠得到整個(gè)模擬區(qū)污染風(fēng)險(xiǎn)分布圖,根據(jù)此圖可以確定污染風(fēng)險(xiǎn)小的新供水井位置.

4.3 本文利用蒙特卡羅的方法與克里格替代模型,大大減少了計(jì)算負(fù)荷,而且克里格替代模型預(yù)測結(jié)果與模擬模型模擬結(jié)果相比,相對誤差均小于 2%,說明替代模型可以較好地模仿模擬模型的行為.

[1] 董文福,傅德黔.近年來我國環(huán)境污染事故綜述 [J]. 環(huán)境科學(xué)與技術(shù), 2009,32(7):75-77.

[2] 王 昊.我國地下水污染成因及防治措施研究 [J]. 資源節(jié)約與保護(hù), 2013(5):88-97.

[3] 王焰新.地下水污染與防治 [M]. 北京:高等教育出版社, 2007.

[4] 張麗君.地下水脆弱性和風(fēng)險(xiǎn)性評價(jià)研究進(jìn)展綜述 [J]. 水文地質(zhì)工程地質(zhì), 2006,33(6):113-119.

[5] 劉增超,董 軍,何連生,等.基于過程模擬的地下水污染風(fēng)險(xiǎn)評價(jià)方法研究 [J]. 中國環(huán)境科學(xué), 2013,33(6):1120-1126.

[6] 馬祿義,許學(xué)工,徐麗芬.中國綜合生態(tài)風(fēng)險(xiǎn)評價(jià)的不確定性分析 [J]. 北京大學(xué)學(xué)報(bào):自然科學(xué)版, 2011(5):893-900.

[7] Goodrich M T, Mccord J T. Quantification of uncertainty in exposure assessments at hazardous waste sites [J]. Ground Water, 1995,33(5):727-732.

[8] Bennett D H, James A L, Mckone T E, et al. On uncertainty in remediation analysis: variance propagation from subsurfacetransport to exposure modeling [J]. Reliability Engineering and System Safety, 1998,62(1/2):117-129.

[9] 梁 婕.基于不確定理論的地下水溶質(zhì)運(yùn)移及污染風(fēng)險(xiǎn)研究[D]. 長沙:湖南大學(xué), 2009.

[10] 申 升.基于貝葉斯理論的地下水溶質(zhì)運(yùn)移的不確定性研究[D]. 長沙:湖南大學(xué), 2012.

[11] 王 宇,盧文喜,卞建民,等.基于小波神經(jīng)網(wǎng)絡(luò)的地下水流數(shù)值模擬模型的替代模型研究 [J]. 中國環(huán)境科學(xué), 2015,35(1):139-146.

[12] 歐陽琦,盧文喜,侯澤宇,等.基于替代模型的地下水溶質(zhì)運(yùn)移不確定性分析 [J]. 中國環(huán)境科學(xué), 2016,36(4):1119-1124.

[13] 曾獻(xiàn)奎,王 棟,吳吉春.地下水流概念模型的不確定性分析 [J].南京大學(xué)學(xué)報(bào)(自然科學(xué)), 2012,48(6):746-752.

[14] 姚磊華.白楊水源地潛水水質(zhì)的 Monte—Carlo隨機(jī)模擬 [J].長春科技大學(xué)學(xué)報(bào), 1998(3):296-302.

[15] Hamed M M, Conte J P, Bedient P B. Probabilistic screening tool for ground-water contamination assessment [J]. Journal of Environmental Engineering, 1995,121(11):767-775.

[16] 束龍倉,劉佩貴,劉 波,等.傍河水源地?cái)?shù)學(xué)模型的參數(shù)靈敏度分析—以遼寧省北票市某傍河水源地為例 [J]. 工程勘察, 2006,(8):29-31.

[17] 束龍倉,王茂枚,劉瑞國,等.地下水?dāng)?shù)值模擬中的參數(shù)靈敏度分析 [J]. 河海大學(xué)學(xué)報(bào)(自然科學(xué)版), 2007,35(5):491-495.

[18] 安永凱,盧文喜,董海彪,等.基于克里格法的地下水流數(shù)值模擬模型的替代模型研究 [J]. 中國環(huán)境科學(xué), 2014,34(4):1073-1079.

[19] 顧文龍,盧文喜,馬洪云,等.地下水?dāng)?shù)值模擬分析中降水入滲補(bǔ)給強(qiáng)度及滲透系數(shù)不確定性評價(jià) [J]. 水電能源科學(xué), 2015(11): 45-48.

[20] 吳義忠.多領(lǐng)域物理系統(tǒng)的仿真優(yōu)化方法 [M]. 北京:科學(xué)出版社, 2011.

Groundwater contamination risk assessment method based on sensitivity analysis and surrogate model.

CHANG Zhen-bo1,2, LU Wen-xi1,2*, XIN Xin1,2, GU Wen-long1,2, CUI Shang-jin1,2
(1.Key Laboratory of Groundwater Resources and Environment Ministry of Education, Jilin University, Changchun 130012, China；2.College of Environment and Resources, Jilin University, Changchun 130012, China). China Environmental Science, 2017,37(1)：167～173

Stochastic model with Monte Carlo method was used to assess the risk of groundwater pollution. The random variables in the model were determined by the sensitivity analysis, which made the result of the groundwater risk assessment more reliable, and the evaluation process was illustrated by a hypothetical example. The results showed that the simulated output data was in accord with the normal distribution law, and normal probability density function of the integrator was the risk of contamination. The pollution risks of well 1, 2 and 3 were 0%, 78.52% and 100%, respectively. Sub regions with different risk of pollution were divided according to the pollution risk distribution map of the whole simulation area, so as to quantitatively evaluate the pollution risk of different sub area in the simulation area.

groundwater contamination risk assessment；sensitivity analysis；Monte Carlo；surrogate model；probability statistics

X703

A

1000-6923(2017)01-0167-07

常振波(1993-),男,河南開封人,吉林大學(xué)碩士研究生,主要從事地下水溶質(zhì)運(yùn)移的不確定性研究.

2016-05-15

國家自然科學(xué)基金項(xiàng)目(41372237);吉林大學(xué)研究生創(chuàng)新基金資助項(xiàng)目(2016207,2016100)

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