羅成明,盧文喜*,王梓博,常振波

地下水污染模擬模型的不確定性分析

羅成明1,2,盧文喜1,2*,王梓博1,2,常振波1,2

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

為同時(shí)分析源匯項(xiàng)和水文地質(zhì)參數(shù)不確定性對(duì)地下水污染數(shù)值模擬模型輸出結(jié)果的影響,以撫順市某煤矸石堆放場(chǎng)為研究實(shí)例展開研究.首先以硫酸根離子作為模擬因子,建立該場(chǎng)地地下水污染數(shù)值模擬模型.然后,分別采用局部靈敏度分析和全局靈敏度分析兩種方法對(duì)模擬模型參數(shù)進(jìn)行靈敏度分析并對(duì)比二者的結(jié)果,最終篩選出對(duì)模型輸出結(jié)果影響較大的兩個(gè)參數(shù)作為模型的隨機(jī)參數(shù).為減少反復(fù)調(diào)用模擬模型產(chǎn)生的計(jì)算負(fù)荷,分別對(duì)3口觀測(cè)井應(yīng)用克里格、核極限學(xué)習(xí)機(jī)、支持向量機(jī)和BP神經(jīng)網(wǎng)絡(luò)4種方法建立模擬模型的替代模型,根據(jù)這4種替代模型在不同井的擬合效果,為每口井優(yōu)選出一個(gè)表現(xiàn)最好的替代模型,并利用優(yōu)選出的替代模型完成蒙特卡洛隨機(jī)模擬.最后,對(duì)隨機(jī)模擬的結(jié)果進(jìn)行統(tǒng)計(jì)分析與風(fēng)險(xiǎn)評(píng)估.結(jié)果表明,在置信度為80%時(shí),1,2,3號(hào)三口井濃度的置信區(qū)間分別為:211.48～845.04mg/L,0～406.98mg/L,231.42～958.37mg/L.此外,依據(jù)《地下水質(zhì)量標(biāo)準(zhǔn)》和各井的污染物濃度分布函數(shù)曲線得出:1號(hào)井和3號(hào)井的水質(zhì)達(dá)標(biāo)Ⅴ類水的概率分別為90.1%和93.1%,2號(hào)井達(dá)標(biāo)Ⅲ類水的概率為80.7%,藉此為地下水資源管理和污染防治提供合理依據(jù).

地下水污染隨機(jī)模擬；靈敏度分析；替代模型；不確定性分析；風(fēng)險(xiǎn)評(píng)估

地下水模擬模型可以用來描述地下水水流運(yùn)動(dòng)和溶質(zhì)運(yùn)移規(guī)律,是地下水資源管理和污染防治的重要技術(shù)手段.然而,在人為因素以及自然因素的雙重影響下,地下水?dāng)?shù)值模擬模型存在不確定性,導(dǎo)致模型輸出結(jié)果存在不確定性[1].

近年來地下水溶質(zhì)運(yùn)移的不確定分析發(fā)展迅速.Li等[2]將不確定性分析應(yīng)用到風(fēng)險(xiǎn)評(píng)估方面,分析了天津市南水北調(diào)工程取水區(qū)淺層和深層地下水的風(fēng)險(xiǎn)值.王梓博等[1]、常振波等[3]通過蒙特卡羅方法分析了水文地質(zhì)參數(shù)不確定性;李久輝等[4]考慮了邊界條件不確定性對(duì)污染質(zhì)運(yùn)移結(jié)果的影響;袁乾等[5]同時(shí)考慮了水文地質(zhì)參數(shù)及源匯項(xiàng)中“源”的不確定性.在實(shí)際情況中,水文地質(zhì)參數(shù)及源匯項(xiàng)的不確定性對(duì)地下水模擬模型輸出結(jié)果的不確定性影響顯著.然而,在目前地下水溶質(zhì)運(yùn)移不確定性的研究中,考慮單因素不確定性的研究較多,而同時(shí)考慮多因素不確定性的研究較少,本文同時(shí)考慮水文地質(zhì)參數(shù)及源匯項(xiàng)中的“匯”進(jìn)行不確定性分析.

在不確定性分析的過程中將會(huì)多次調(diào)用模擬模型,這會(huì)帶來極大的計(jì)算負(fù)荷,而替代模型不僅可以大幅減少計(jì)算負(fù)荷還可以保證較高的逼近精度.此前,高鑫宇等[6]將改進(jìn)稀疏網(wǎng)格替代模型應(yīng)用于地下水DNAPLs運(yùn)移數(shù)值模擬的不確定性中,降低了計(jì)算耗時(shí).侯澤宇等[7]同時(shí)應(yīng)用克里格、核極限學(xué)習(xí)機(jī)和支持向量機(jī)3種方法建立集對(duì)權(quán)組合替代模型,大幅減少模擬-優(yōu)化過程的計(jì)算負(fù)荷.葛淵博等[8]將機(jī)器學(xué)習(xí)方法引用于不確定分析中,利用BP神經(jīng)網(wǎng)絡(luò)方法建立替代模型,保證了較高的精度.在前人的研究中,對(duì)于場(chǎng)地多口觀測(cè)井替代模型的建立采用的都是相同的方法.然而,在研究問題較為復(fù)雜且存在多口監(jiān)測(cè)井的情況下,某種單一方法所建立的替代模型并不具備普適性,很難同時(shí)滿足所有井的精度要求.針對(duì)該問題,本文分別對(duì)研究區(qū)內(nèi)3口井利用克里格、核極限學(xué)習(xí)機(jī)、支持向量機(jī)和BP神經(jīng)網(wǎng)絡(luò)4種方法建立替代模型,經(jīng)過對(duì)比分析,為每口井分別優(yōu)選出表現(xiàn)最高的一種替代模型,以此提高模擬模型的逼近精度.

綜上,本文以撫順市某煤矸石堆放場(chǎng)為研究實(shí)例,同時(shí)分析源匯項(xiàng)和水文地質(zhì)參數(shù)不確定性對(duì)地下水污染數(shù)值模擬模型輸出結(jié)果的影響.根據(jù)場(chǎng)地的實(shí)際情況及污染物特征,本文選取硫酸根離子為模擬因子,建立該場(chǎng)地的地下水污染數(shù)值模擬模型.同時(shí)應(yīng)用局部和全局靈敏度分析兩種方法對(duì)參數(shù)進(jìn)行靈敏度分析并對(duì)比二者的分析結(jié)果,篩選出靈敏度較大的兩個(gè)參數(shù)作為模擬模型的隨機(jī)參數(shù).為減少多次調(diào)用模擬模型產(chǎn)生的巨大計(jì)算負(fù)荷,分別應(yīng)用克里格、核極限學(xué)習(xí)機(jī)、支持向量機(jī)和BP神經(jīng)網(wǎng)絡(luò)4種方法建立模擬模型的替代模型,根據(jù)這四種替代模型在不同井的擬合效果,分析對(duì)比篩選出場(chǎng)地各井精度最高的替代模型,并利用優(yōu)選出的替代模型完成蒙特卡洛隨機(jī)模擬.最后對(duì)該研究區(qū)隨機(jī)模擬的結(jié)果進(jìn)行統(tǒng)計(jì)分析與區(qū)間估計(jì),進(jìn)行地下水污染風(fēng)險(xiǎn)評(píng)估,為地下水污染防治提供科學(xué)依據(jù).

1 研究方法

1.1 蒙特卡羅法

蒙特卡羅方法是以概率統(tǒng)計(jì)理論為基礎(chǔ)的一種方法.在具體的問題中,通常采用成百數(shù)千次的模擬,統(tǒng)計(jì)得到各個(gè)事件發(fā)生的概率或者各事件觀測(cè)值的算術(shù)平均值,借此得到問題的解.利用蒙特卡羅方法進(jìn)行不確定性研究可將地下水?dāng)?shù)值模擬各因素(參數(shù),源匯項(xiàng)等)的不確定性轉(zhuǎn)化成輸出結(jié)果的不確性,利于問題的研究.

1.2 靈敏度分析法

利用靈敏度分析法來篩選出對(duì)模擬模型結(jié)果影響最大的參數(shù),以該參數(shù)為隨機(jī)變量不僅會(huì)減少替代模型的維度,還可以減少計(jì)算負(fù)載[9].在目前的靈敏度分析方法中,包括局部靈敏度分析和全局靈敏度分析[10].本文將采用這兩種方法,分別求得地下水模擬模型參數(shù)的靈敏度,并將二者的結(jié)果進(jìn)行對(duì)比分析,局部靈敏度分析方法作為試驗(yàn)方法,全局靈敏度方法作為驗(yàn)證方法,更準(zhǔn)確的選出靈敏度較大的參數(shù).

(1)局部靈敏度分析主要目的是選擇影響最大的因素,其公式如下:

式中:X為靈敏度系數(shù);為了方便計(jì)算以及不同參數(shù)靈敏度的比較,通常采用無量綱的形式:

式中:X表示參數(shù)α變化時(shí)對(duì)輸出結(jié)果的影響程度,X越大說明參數(shù)的變化對(duì)輸出結(jié)果影響越大.

(2)全局靈敏度分析是基于局部靈敏度分析的不足進(jìn)行改進(jìn).相比于局部靈敏度分析,全局靈敏度能夠考慮不同參數(shù)對(duì)研究對(duì)象的共同影響,這也更加的契合于實(shí)際問題.目前全局靈敏度的分析方法較多,如基于差分分析的sobol法[11],多元回歸法[12]以及Morris方法[13]等.本文采用Morris方法進(jìn)行模型參數(shù)的全局靈敏度分析.Morris方法著重在于參數(shù)的取值以及依據(jù)參數(shù)的重要性進(jìn)行排序[14],該方法取值抽樣時(shí),一次只改變一個(gè)參數(shù),假設(shè)地下水模型需要分析的參數(shù)個(gè)數(shù)為個(gè),則會(huì)輪流計(jì)算個(gè)參數(shù)的目標(biāo)函數(shù)值[15].根據(jù)以下公式可以計(jì)算每個(gè)參數(shù)的基效應(yīng):

EE為第個(gè)參數(shù)變化引起的變化效應(yīng).為了方便不同參數(shù)的對(duì)模型的影響度進(jìn)行比較,對(duì)影響度進(jìn)行修正:

(4)

每個(gè)參數(shù)可有多種取值,最終取均值:

式中:D為某一參數(shù)的變化;為單個(gè)參數(shù)單次取樣的影響值;為單個(gè)參數(shù)多次取樣后的影響均值.

1.3 替代模型

替代模型在多學(xué)科領(lǐng)域內(nèi)應(yīng)用廣泛,其基本原理是通過已知的輸入-輸出樣本集構(gòu)造同模型相擬合的近似函數(shù),借此得到模型的輸出結(jié)果.在分析地下水模擬不確定性的過程中,為了減少反復(fù)調(diào)用模擬模型產(chǎn)生的大量計(jì)算負(fù)荷選擇建立模擬模型的替代模型.本文分別采用BP神經(jīng)網(wǎng)絡(luò)(BPNN)方法、核極限學(xué)習(xí)機(jī)(KELM)方法、支持向量機(jī)(SVR)方法以及克里格(KRG)方法建立模擬模型的替代模型,并分別對(duì)比4種方法,針對(duì)各個(gè)觀測(cè)井的精度好壞,選出每個(gè)觀測(cè)井精度最好的替代模型,進(jìn)行后續(xù)的隨機(jī)模擬研究.

1.3.1 BP神經(jīng)網(wǎng)絡(luò)方法 BP神經(jīng)網(wǎng)絡(luò)的核心是根據(jù)每次訓(xùn)練得到的結(jié)果與預(yù)想結(jié)果進(jìn)行誤差分析,之后修改閾值和權(quán)重,不斷反復(fù)多次的訓(xùn)練,進(jìn)而得到同模擬模型輸出值相吻合的結(jié)果.淺顯的理解就是根據(jù)反饋不斷更新迭代,以達(dá)到優(yōu)化的目的.BP神經(jīng)網(wǎng)絡(luò)通常由三層組成:輸入層、隱含層和輸出層;計(jì)算過程分為兩個(gè)階段:信息正向傳播與誤差的反向傳播過程[16].其簡易結(jié)構(gòu)如圖1所示.

(1)信息正向傳播過程:讓信息通過輸入層進(jìn)入網(wǎng)絡(luò),每一層網(wǎng)絡(luò)都將對(duì)其進(jìn)行計(jì)算,然后得到最終輸出結(jié)果的過程BP神經(jīng)網(wǎng)絡(luò)常采用的傳遞函數(shù)Sigmoid函數(shù),此函數(shù)是非線性變換函數(shù),一般采用:

正向傳輸?shù)妮敵鰧佑?jì)算公式如下:

式中:O代表神經(jīng)元的輸出;W為神經(jīng)元與神經(jīng)元連接的權(quán)重.

(2)逆向傳播過程:逆向反饋從最后一層即輸出層開始,第一次的正向反饋的權(quán)重隨機(jī)賦值,為了達(dá)到目標(biāo)優(yōu)化的效果,此時(shí)需要調(diào)整整個(gè)網(wǎng)絡(luò)的參數(shù)即權(quán)重值,對(duì)網(wǎng)絡(luò)調(diào)整的依據(jù)是輸出層輸出值與類別之間的差異.

對(duì)于輸出層:

式中:E表示第個(gè)節(jié)點(diǎn)的誤差值;O表示第個(gè)節(jié)點(diǎn)的輸出值;T記錄輸出值.

中間隱層通過下層節(jié)點(diǎn)總誤差按權(quán)重累加:

式中:E為下一層節(jié)點(diǎn)的誤差率;W為當(dāng)層節(jié)點(diǎn)到下層節(jié)點(diǎn)的權(quán)重值.

計(jì)算誤差后可利用誤差率對(duì)權(quán)重進(jìn)行更新:

式中:h為學(xué)習(xí)速率,取值范圍在(0.1),一般小于0.5[17].

1.3.2 核極限學(xué)習(xí)機(jī)方法 核極限學(xué)習(xí)機(jī)(KELM)基于極限學(xué)習(xí)機(jī)(ELM),將核函數(shù)引入后者的訓(xùn)練中,將傳統(tǒng)的隨機(jī)映射替換成映射,大幅提高了后者的擬合及分類的能力[18].KELM的基本原理如下:

對(duì)于訓(xùn)練樣本,ELM的輸出y可表示為:

經(jīng)過一系列的優(yōu)化及化簡,KELM輸出函數(shù)的表達(dá)公式如下:

1.3.3 克里格方法 克里格方法的核心是在范圍有限的研究區(qū)內(nèi),利用變量的相互契合性和隨機(jī)變化性搭建成該方法的基本結(jié)構(gòu),然后對(duì)區(qū)域中變量的取值進(jìn)行誤差最小的優(yōu)化估計(jì).該方法將輸入與輸出變量之間的關(guān)系以回歸方程的形式表達(dá)出來,其公式如下:

目前常用的模型為高斯模型,模型如下:

1.3.4 支持向量機(jī)方法 支持向量機(jī)所采用的并非為傳統(tǒng)統(tǒng)計(jì)學(xué)上的經(jīng)驗(yàn)風(fēng)險(xiǎn)最小化原理,而是一種全新的通用機(jī)器學(xué)習(xí)方法:結(jié)構(gòu)風(fēng)險(xiǎn)最小化原理,該種方法也更加的貼合目前階段的科學(xué)研究.支持向量機(jī)的核心在于回歸的思想,而回歸的核心是借以非線性映射函數(shù)把輸入數(shù)據(jù)映射到高維空間進(jìn)行線性回歸.方程可表示為:

最后得到:

2 應(yīng)用研究

2.1 概況

本文以撫順市某煤矸石堆至渾河圈定的地區(qū)作為本次地下水?dāng)?shù)值模擬的研究對(duì)象.研究區(qū)東西長約9700m,南北寬約6300m,總面積約52.35km2,其地勢(shì)東高西低,地下水流向大體上自東向西.將研究區(qū)概化為非均質(zhì)各向同性含水層,主要研究對(duì)象為一層潛水含水層,厚度約8m,水流為二維穩(wěn)定流,根據(jù)研究區(qū)含水層的滲透性可分為兩個(gè)區(qū)域.研究區(qū)示意圖詳見圖2.

圖2 研究區(qū)概況

研究區(qū)北部邊界AB為渾河,概化為已知水頭邊界;東北部邊界BC和南部邊界DE平行于區(qū)內(nèi)地下水流線,概化為零通量邊界;將研究區(qū)東部邊界CD概化為側(cè)向徑流補(bǔ)給邊界;西南部邊界EA概化為側(cè)向徑流排泄邊界,側(cè)向徑流邊界概化為已知流量邊界.在溶質(zhì)運(yùn)移模型中,研究區(qū)北部邊界AB和西南部邊界EA概化為已知對(duì)流-彌散通量邊界;東北部邊界BC和南部邊界DE由于沒有對(duì)流作用,概化為零通量的水動(dòng)力彌散通量邊界;東部邊界CD概化為已知濃度邊界.將歷史水位資料的均值作為水流模型的初始條件;將硫酸根離子的環(huán)境基底值(零濃度)作為溶質(zhì)運(yùn)移模型的初始條件.

此次的研究考慮源匯項(xiàng)的影響,根據(jù)實(shí)際的踏勘情況以及當(dāng)?shù)鼐用竦挠盟闆r設(shè)置三口觀測(cè)井兼抽水井(W1～W3).考慮實(shí)際場(chǎng)地污染物所占的比重以及物理化學(xué)性質(zhì),選定硫酸根離子作為此次污染物研究對(duì)象.對(duì)實(shí)際場(chǎng)地的地下水流模型以及溶質(zhì)運(yùn)移模型的校正得表1參數(shù)取值,污染源硫酸根離子的釋放強(qiáng)度為7000mg/d.抽水井用水量取值如表2所示.

表1 模型水文地質(zhì)參數(shù)校正后取值

表2 抽水井年用水量

2.2 模型建立

2.2.1 數(shù)值模擬 充分考慮研究區(qū)的水文地質(zhì)條件,并在地下水概念模型的基礎(chǔ)上,建立地下水水流模型和溶質(zhì)運(yùn)移模型,后者是在前者的基礎(chǔ)上建立的.

(1)水流模型:

(2)溶質(zhì)運(yùn)移模型:

本文的地下水水流模型和溶質(zhì)運(yùn)移模型是運(yùn)用GMS中MODEFLOW板塊和MT3D板塊求解.給模型輸入校正后的水文地質(zhì)參數(shù),抽水井W1的抽水量為18250m3/a, 抽水井W2的抽水量為25550m3/a, 抽水井W1的抽水量為36500m3/a,預(yù)報(bào)20a后的地下水污染狀況(圖3).

圖3 污染物運(yùn)移分布

2.2.2 靈敏度分析 篩選隨機(jī)變量及方法對(duì)比本文將同時(shí)采用局部靈敏度分析和全局靈敏度分析,篩選出對(duì)模型影響較大的參數(shù).

(1)局部靈敏度分析:對(duì)于污染場(chǎng)地的三口監(jiān)測(cè)井,將孔隙度、滲透系數(shù)、降雨入滲系數(shù)以及縱向彌散度L這4個(gè)參數(shù)在波動(dòng)范圍內(nèi)取均值輸入到模擬模型中,得到一組標(biāo)準(zhǔn)輸出值,后分別將各參數(shù)增加和減少10%,20%,且保證每次只變化一個(gè)參數(shù),控制其他的參數(shù)不變,輸入模型中得到污染物濃度的輸出值,利用式(2)計(jì)算得到三口井各個(gè)參數(shù)的靈敏度,然后將三口井的參數(shù)靈敏度取均值,篩選出靈敏度最大的兩個(gè)參數(shù),結(jié)果如表3所示.

表3 靈敏度計(jì)算結(jié)果

根據(jù)3口監(jiān)測(cè)井平均靈敏度分析示意圖(圖4),局部靈敏度分析法所篩選出的參數(shù)分別是孔隙度和滲透系數(shù).

圖4 三口監(jiān)測(cè)井平均靈敏度分析示意

(2)全局靈敏度分析:同局部靈敏度分析選取同樣的水文地質(zhì)參數(shù),進(jìn)行全局靈敏度分析,每個(gè)參數(shù)分別有3種取值,即參數(shù)原值以及增加和減少10%的值.將4個(gè)參數(shù)進(jìn)行組合,共有15種組合方式,如表4所示.

將表4的15組組合的參數(shù)變化后輸入到模擬模型中,得到污染物的濃度值作為輸出,為了方便不同參數(shù)靈敏度的比較,采用公式(3)、(4)、(5)計(jì)算得到各組在不同取值下的靈敏度,單組取多次平均值,即可得15組參數(shù)組合的全局靈敏度.

表4 參數(shù)組合方式

圖5 參數(shù)全局靈敏度分析

如圖5所示,靈敏度最高的兩個(gè)參數(shù)分別是孔隙度和滲透系數(shù)該結(jié)果與局部靈敏度分析結(jié)果一致.但考慮不同參數(shù)的相互作用時(shí),其結(jié)果存在差異,如圖5中的組合6,在同時(shí)考慮孔隙度和滲透系數(shù)變化時(shí),整體的靈敏度比二者單獨(dú)的靈敏度要低,再如第15組,同時(shí)考慮四個(gè)參數(shù)同時(shí)變化所帶來的影響,其結(jié)果并非是單參數(shù)影響的疊加,整體而言全局靈敏度分析更具實(shí)際意義.利用兩種靈敏度分析方法(前者為試驗(yàn)方法,后者為檢驗(yàn)方法)進(jìn)行參數(shù)篩選其結(jié)果更具可靠性.

利用局部靈敏度分析和全局靈敏度分析共同篩選出對(duì)模型影響最大的兩個(gè)參數(shù):孔隙度和滲透系數(shù)將這兩個(gè)參數(shù)作為隨機(jī)變量,利用拉丁超立方抽樣法進(jìn)行抽樣,其余的參數(shù)作為確定性變量.

2.2.3 拉丁超立方抽樣 本文采用拉丁超立方抽樣法分別在孔隙度、滲透系數(shù)以及三口抽水井抽水量的合理變化范圍內(nèi)進(jìn)行抽樣,其中考慮場(chǎng)地分區(qū)的差異性,對(duì)場(chǎng)地孔隙度及滲透系數(shù)進(jìn)行2次抽樣,分配給東西兩個(gè)區(qū)域,后將所有抽樣值隨機(jī)組合,抽樣的區(qū)間及服從的分布狀態(tài)如表5所示.該抽樣過程通過MATLAB軟件實(shí)現(xiàn),對(duì)上述參數(shù)分別抽取110組,隨機(jī)組合得到輸入樣本集.

表5 參數(shù)概率分布及取值情況

2.2.4 建立替代模型 將上述2個(gè)分區(qū)的2個(gè)水文地質(zhì)參數(shù)以及3口抽水井的抽水量共7維數(shù)據(jù)輸入到地下水模擬模型中,運(yùn)用GMS求解得到110組污染物的濃度即輸出值.其中前100組作為訓(xùn)練樣本建立替代模型,后10組作為檢驗(yàn)樣本檢驗(yàn)替代模型的精度.

本文3口監(jiān)測(cè)井分別采用BP神經(jīng)網(wǎng)絡(luò)(BPNN)方法、核極限學(xué)習(xí)機(jī)(KELM)方法、支持向量機(jī)(SVR)方法以及克里格(KRG)方法建立替代模型,選用確定性系數(shù)2以及平均相對(duì)百分比誤差(MAPE)檢驗(yàn)替代模型的精度,分別挑選出各井最適合的替代模型.

如下為評(píng)定精度指標(biāo)的計(jì)算公式:

(1)確定性系數(shù)2:2的值越接近1則代表替代模型越精確.

(2)平均相對(duì)百分比誤差(MAPE):表示替代模型的輸出結(jié)果同模擬模型輸出結(jié)果的平均偏差.

利用2和MAPE兩個(gè)指標(biāo)分別評(píng)定各個(gè)抽水井不同替代模型的精度,其結(jié)果如表6所示.

采用2進(jìn)行精度評(píng)價(jià)時(shí):通過計(jì)算可知抽水井W1最適用的替代模型是克里格替代模型,2達(dá)0.8214;抽水井W2最適用的替代模型是BP神經(jīng)網(wǎng)絡(luò)替代模型,2達(dá)0.9965;抽水井W3最適用的替代模型是克里格替代模型,2達(dá)0.9972.以2為評(píng)價(jià)指標(biāo)為研究區(qū)3口井選定的替代模型依次為:克里格、BP神經(jīng)網(wǎng)絡(luò)及克里格替代模型.采用MAPE進(jìn)行精度評(píng)價(jià)時(shí):計(jì)算可知抽水井W1最適用的替代模型是核極限學(xué)習(xí)機(jī)替代模型,MAPE達(dá)10.49%;抽水井W2最適用的替代模型是BP神經(jīng)網(wǎng)絡(luò)替代模型,MAPE達(dá)4.12%;抽水井W3最適用的替代模型是克里格替代模型,MAPE達(dá)1.16%.以MAPE為評(píng)價(jià)指標(biāo)為研究區(qū)3口井選定的替代模型依次為:核極限學(xué)習(xí)機(jī)、BP神經(jīng)網(wǎng)絡(luò)及克里格替代模型.

表6 三口井不同替代模型精度

2.2.5 替代模型精度檢驗(yàn) 為驗(yàn)證3口井單獨(dú)建立替代模型后的綜合精度,本文將通過如下步驟進(jìn)行分析:(1)場(chǎng)地的3口井同時(shí)應(yīng)用一種方法建立替代模型,共4種方法,因此建立4組替代模型;(2)井1,2,3分別建立克里格、BP神經(jīng)網(wǎng)絡(luò)及克里格替代模型;(3)井1,2,3分別建立核極限學(xué)習(xí)機(jī)、BP神經(jīng)網(wǎng)絡(luò)及克里格替代模型;(4)以確定性系數(shù)2及平均相對(duì)百分比誤差(MAPE)為評(píng)定指標(biāo)評(píng)價(jià)上述各方法建立的替代模型的綜合精度,結(jié)果如表7所示.

表7 不同類型替代模型精度

通過對(duì)比分析發(fā)現(xiàn),單獨(dú)建立各井的最優(yōu)替代模型,其模擬精度有顯著提升.單替代模型中MAPE最高的是BP神經(jīng)網(wǎng)絡(luò)替代模型,為6.68%,2最高的是克里格替代模型、支持向量機(jī)及核極限學(xué)習(xí)機(jī)替代模型,為0.96.依據(jù)單井最優(yōu)建立的KRG-BPNN- KRG替代模型和KELM-BPNN-KRG替代模型的2接近,為0.97;前者M(jìn)APE為5.45%,后者M(jìn)APE為5.26%,整體的精度提升顯著.但KELM-BPNN-KRG替代模型的MAPE值更低,因此選用KELM- BPNN-KRG替代模型進(jìn)行后續(xù)的隨機(jī)模擬.

2.2.6 蒙特卡羅隨機(jī)模擬 本次研究采用蒙特卡羅法對(duì)研究場(chǎng)地的地下水污染進(jìn)行隨機(jī)模擬.第一步利用拉丁超立方抽樣方法對(duì)前文敘述的7維數(shù)據(jù)進(jìn)行抽樣,得到1000組參數(shù),然后將參數(shù)依次輸入到KELM-BPNN-KRG替代模型中得到3口監(jiān)測(cè)井污染物濃度輸出值,最后對(duì)各井所對(duì)應(yīng)的輸出結(jié)果進(jìn)行統(tǒng)計(jì)分析.

3 結(jié)果與討論

3.1 統(tǒng)計(jì)與分析

利用3口監(jiān)測(cè)井的1000組數(shù)據(jù)繪制污染物濃度累計(jì)頻次直方圖并對(duì)其進(jìn)行統(tǒng)計(jì)分析.如圖6所示,1號(hào)井污染物濃度值集中在350～650mg/L,2號(hào)井污染物濃度值集中在50～270mg/L,3號(hào)井污染物濃度值集中在400～800mg/L.3口監(jiān)測(cè)井的污染物濃度均值分別為528.26mg/L、167.83mg/L、594.90mg/L.其中3號(hào)井的污染最為嚴(yán)重,這是因?yàn)?號(hào)距離污染源較近且該井的抽水量最大,在一定程度上影響了污染物的運(yùn)移;1號(hào)井同污染源的距離與3號(hào)井同污染源的距離接近,但其抽水量較少,其污染物的濃度比3號(hào)井污染物的濃度略低;2號(hào)井的污染物濃度最低,首先是因?yàn)槠渚辔廴驹摧^遠(yuǎn),其次3號(hào)井位于其東側(cè),3號(hào)井抽水量較大,使得后續(xù)污染物向西運(yùn)移的濃度降低.

表8 三口監(jiān)測(cè)井污染物輸出值統(tǒng)計(jì)指標(biāo)

3.2 區(qū)間估計(jì)

區(qū)間估計(jì)的主要目的是對(duì)未知參數(shù)真值范圍在一定置信水平下進(jìn)行估計(jì),相比確定的結(jié)果,該種方法更具有實(shí)際意義.本文應(yīng)用切比雪夫不等式,對(duì)污染場(chǎng)地3口監(jiān)測(cè)井污染物濃度在一定置信水平下進(jìn)行區(qū)間估計(jì)[5].

式中:()是各時(shí)段監(jiān)測(cè)井污染物濃度均值;()為污染物濃度值的方差;是置信水平;為任意給定的正數(shù).利用公式(20)計(jì)算得到三口監(jiān)測(cè)井不同置信水平下污染物濃度的區(qū)間估計(jì)值,如表9所示,分析發(fā)現(xiàn),置信水平越高,區(qū)間范圍越大,根據(jù)實(shí)際的研究可以采取不同置信度下的污染物濃度范圍.

表9 各井不同置信水平污染物濃度值區(qū)間估計(jì)

3.3 風(fēng)險(xiǎn)評(píng)估

表10 污染物達(dá)標(biāo)概率預(yù)測(cè)

對(duì)1000組監(jiān)測(cè)數(shù)據(jù)利用MATLAB軟件繪制污染物濃度分布函數(shù)曲線,進(jìn)行結(jié)果分析,其中縱坐標(biāo)代表分布的概率,橫坐標(biāo)為濃度值(圖7).考慮到污染場(chǎng)地的地理位置且有較多居民居住,參照《地下水質(zhì)量標(biāo)準(zhǔn)》(GB/T 14848-2017)[19],利用概率分布曲線分析評(píng)定3口監(jiān)測(cè)井中硫酸根離子達(dá)各級(jí)標(biāo)準(zhǔn)的概率,如表10所示.分析可知,硫酸根離子污染20a后,3口井達(dá)到Ⅰ類標(biāo)準(zhǔn)的概率依次為:0%、11%、0%;達(dá)到Ⅱ類標(biāo)準(zhǔn)的概率依次為:0%、40%、0%;達(dá)到Ⅲ類標(biāo)準(zhǔn)的概率依次為:2.2%、29.7%、1.6%;達(dá)到Ⅳ類標(biāo)準(zhǔn)的概率依次為:7.7%、13.1%、5.3%;達(dá)到Ⅴ類標(biāo)準(zhǔn)的概率依次為:90.1%、6.2%、93.1%.通過如上分析可知1號(hào)井和3號(hào)井污染嚴(yán)重,不易作為飲用水,2號(hào)井經(jīng)適當(dāng)處理可作為生活用水.

4 結(jié)論

4.1 利用局部靈敏度分析方法和全局靈敏度分析方法同時(shí)篩選出研究區(qū)靈敏度較大的兩個(gè)參數(shù):孔隙度和滲透系數(shù),其中局部靈敏度分析方法作為試驗(yàn)方法,全局靈敏度方法作為驗(yàn)證方法,這使得篩選出的參數(shù)更加可靠且具有實(shí)際意義.

4.2 單獨(dú)建立各井的最優(yōu)替代模型,其模擬精度有顯著提升.經(jīng)過優(yōu)選1,2,3號(hào)井的替代模型依次為:核極限學(xué)習(xí)機(jī)、BP神經(jīng)網(wǎng)絡(luò)及克里格替代模型,通過此方法既減少了反復(fù)調(diào)用模擬模型產(chǎn)生的巨大的計(jì)算負(fù)荷,又保持了較好的精度.

4.3 本文同時(shí)考慮了水文地質(zhì)參數(shù)和源匯項(xiàng)的不確定性,對(duì)這兩個(gè)因素應(yīng)用蒙特卡羅方法對(duì)地下水?dāng)?shù)值模擬模型進(jìn)行了不確定性分析,并估算了20年后各井水質(zhì)達(dá)標(biāo)情況.其中1號(hào)井和3號(hào)井達(dá)到Ⅴ類水的概率分別為90.1%和93.1%,污染嚴(yán)重.2號(hào)井達(dá)標(biāo)Ⅲ類水的概率為80.7%,適當(dāng)處理可作為生活用水.借此分析為地下水資源的利用以及污染防治提供合理依據(jù).

5 展望

5.1 由于硫酸根離子在地下水環(huán)境中吸附和降解的作用較微弱,故本文只考慮了硫酸根離子保守遷移的情況,將在未來的研究中全面考慮該種情況.

5.2 一般來說,對(duì)非均質(zhì)性越強(qiáng)的含水層建立替代模型,其精度會(huì)有一定程度的下降.本文根據(jù)水文地質(zhì)條件實(shí)際情況,將研究區(qū)劃分為兩個(gè)滲透系數(shù)分區(qū),非均質(zhì)性較弱,在日后的研究中將探索高非均質(zhì)性含水層建模任務(wù).

[1] 王梓博,盧文喜,王 涵,等.某鉬礦尾礦庫地下水污染的隨機(jī)模擬[J]. 中國環(huán)境科學(xué), 2020,40(5):2124-2131.

Wang Z B, Lu W X, Wang H, et al.Stochastic simulation of groundwater pollution in a molybdenum tailings reservoir [J]. China Environmental Science, 2020,40(5):2124-2131.

[2] Li F, Zhu J, Deng X, et al. Assessment and Uncertainty Analysis of Groundwater Risk Lifawem [J]. Environmental Research, 2018,160: 140-151.

[3] 常振波,盧文喜,辛 欣,等.基于靈敏度分析和替代模型的地下水污染風(fēng)險(xiǎn)評(píng)價(jià)方法[J]. 中國環(huán)境科學(xué), 2017,37(1):167-173.

Chang Z B, Lu W X, Xin X, et al. Groundwater pollution risk assessment method based on sensitivity analysis and surrogate model [J]. China Environmental Science, 2017,37(1):167-173.

[4] 李久輝,盧文喜,辛 欣,等.考慮邊界條件不確定性的地下水污染風(fēng)險(xiǎn)分析[J]. 中國環(huán)境科學(xué), 2018,38(6):2167-2174.

Li J H, Lu W X, Xin X, et al. Groundwater pollution risk analysis considering uncertainty of boundary conditions [J]. China Environmental Science, 2018,38(6):2167-2174.

[5] 袁 乾,盧文喜,范 越,等.基于替代模型的煤矸石堆地下水污染隨機(jī)模擬[J]. 中國環(huán)境科學(xué), 2019,39(6):2444-2451.

Yuan Q, Lu W X, Fan Y, et al. Stochastic simulation of groundwater pollution in gangue dump based on surrogate model [J]. China Environmental Science, 2019,39(6):2444-2451.

[6] 高鑫宇,曾獻(xiàn)奎,吳吉春.基于改進(jìn)稀疏網(wǎng)格替代模擬的地下水DNAPLs運(yùn)移不確定性分析[J]. 水文地質(zhì)工程地質(zhì), 2020,47(1):1-10.

Gao X Y, Zeng X K, Wu J C. Uncertainty analysis of groundwater DNAPLs migration based on improved sparse grid instead of simulation [J]. Hydrogeology Engineering Geology, 2020,47(1):1-10.

[7] 侯澤宇,王 宇,盧文喜.地下水DNAPLs污染修復(fù)多相流模擬的替代模型[J]. 中國環(huán)境科學(xué), 2019,39(7):2913-2920.

Hou Z Y, Wang Y, Lu W X. An surrogate model for simulation of multiphase flow for groundwater DNA PLs contamination remediation [J]. China Environmental Science, 2019,39(7):2913-2920.

[8] 葛淵博,盧文喜,王梓博,等.基于BP神經(jīng)網(wǎng)絡(luò)替代模型的地下水污染隨機(jī)模擬[J]. 中國農(nóng)村水利水電, 2022,(3):107-113,119.

Ge Y B, Lu W X, Wang Z B, et al. Stochastic simulation of groundwater pollution based on BP neural network surrogate model [J]. Rural Water Resources and Hydropower in China, 2022,(3):107-113,119.

[9] 王 涵,盧文喜,李久輝,等.地下水DNAPLs污染多相流的隨機(jī)模擬及其不確定性分析[J]. 中國環(huán)境科學(xué), 2018,38(7):2572-2579.

Wang H, Lu W X, Li J H, et al. Stochasticsimulation of multiphase flow contaminated by groundwater DNAPLs and analysis of uncertainty [J]. China Environmental Science, 2018,38(7):2572-2579.

[10] Luo J, Lu W. Sobol'sensitivity analysis of NAPL-contaminated aquifer remediation process based on multiple surrogates [J]. Computers & Geosciences, 2014,67:110-116.

[11] SOBOL'I M. Sensitivity estimates for nonlinear mathematicalmodels [J]. Mathematical Modelling Computer Experiments, 1993,1(4):407-414.

[12] 王浩昌,杜鵬飛,趙冬泉,等.城市降雨徑流模型參數(shù)全局靈敏度分析[J]. 中國環(huán)境科學(xué), 2008,(8):725-729.

Wang H C, Du P F, Zhao D Q, et al. Global sensitivity analysis of urban rainfall-runoff model parameters [J]. China Environmental Science, 2008,(8):725-729.

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

Shu L C, Wang M M, Liu R G, et al. Sensitivity analysis of parameters in groundwater numerical simulation [J]. Journal of Hohai University (Natural Science Edition), 2007,(5):491-495.

[14] 劉 波,金愛兵,高永濤,等.基于Morris法的單裂隙巖體溫度場(chǎng)參數(shù)靈敏度分析[J]. 采礦與安全工程學(xué)報(bào), 2016,33(1):152-157.

Liu B, Jin A B, Gao Y T, et al. Sensitivity analysis of temperature field parameters of single fractured rock mass based on Morris method [J]. Journal of Mining and Safety Engineering, 2016,33(1):152-157.

[15] 賈順卿,盧文喜,李久輝,等.基于U-D分解卡爾曼濾波地下水污染源溯源辨識(shí)[J]. 中國環(huán)境科學(xué), 2021,41(2):713-719.

Jia S Q, Lu W X, Li J H, et al. Inverseidentification of groundwater pollution source based on U-D decomposition Kalman filter [J]. China Environmental Science , 2021,41(2):713-719.

[16] 李 萍,曾令可,稅安澤,等.基于MATLAB的BP神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)系統(tǒng)的設(shè)計(jì)[J]. 計(jì)算機(jī)應(yīng)用與軟件, 2008,(4):149-150,184.

Li P, Zeng L K, Shui A Z, et al. Design of BP neural network prediction system based on MATLAB [J]. Computer Applications and Software, 2008,(4):149-150,184.

[17] 潘紫東,盧文喜,范 越,等.基于模擬-優(yōu)化方法的地下水污染源溯源辨識(shí)[J]. 中國環(huán)境科學(xué), 2020,40(4):1698-1705.

Pan Z D, Lu W X, Fan Y, et al. Inverseidentification of groundwater pollution sources based on simulation-optimization method [J]. China Environmental Science, 2020,40(4):1698-1705.

[18] Wu S, Wang Y, Cheng S . Extreme learning machine based wind speed estimation and sensorless control for wind turbine power generation system [J]. Neurocomputing, 2013,102(FEB.15):163-175.

[19] GB/T 14848-2017 地下水質(zhì)量標(biāo)準(zhǔn)[S].

GB/T 14848-2017 Groundwater quality standard [S].

Uncertainty analysis of groundwater pollution simulation model.

LUO Cheng-ming1,2, LU Wen-xi1,2*, WANG Zi-bo1,2, CHANG Zhen-bo1,2

(1.Key Laboratory of Groundwater Resources and Environmental, Ministry of Education, Jilin University, Changchun 130012, China；2.College of New Energy and Environment, Jilin University, Changchun 130012, China)., 2022,42(7)：3224～3233

In order to simultaneously analyze the effect of source and sink items and hydrogeological parameters uncertainty on the output of groundwater pollution numerical simulation model, a gangue dump site in Fushun City was taken as an example to study. Firstly, the groundwater pollution numerical simulation model of the site was established with sulfate ion as the simulation factor. Then, local sensitivity analysis and global sensitivity analysis were used to analyze the sensitivity of the simulation model parameters, and the results of thetwowere compared. Finally, two parameters that have a great impact on the model output were selected as the random parameters of the model.In order to reduce the calculation load caused by repeatedly calling the simulation model, four methods of Kriging(KRG), kernel extremelearning machine(KELM), support vector machine(SVM) and BP neural network (BPNN) were used to establish the surrogate model of the simulation model for three observation wells respectively. According to the fitting effect of the four surrogate models in different wells, a surrogate model with the best performance was selected for each well, Monte Carlo stochastic simulation was completed by using the optimized surrogate model. Finally, statistical analysis and risk assessment were carried out on the results of random simulation. The results show that when the confidence was 80%, the confidence intervals of the concentration values of the well 1,2,3 were 211.48～845.04mg/L, 0～406.98mg/L, 231.42～958.37mg/L. In addition, according to the groundwater quality standard and the pollutant concentration distribution function curve of each well, the probability that the water quality of well 1 and well 3 met the class V water standard was 90.1% and 93.1% respectively, and the probability that well 2 met the class III water standard was 80.7%, so as to provide a reasonable basis for groundwater resource management and pollution prevention and control.

stochastics simulation of groundwater pollution；sensitivity analysis；surrogate model；uncertainty analysis；risk assessment

X523 文獻(xiàn)識(shí)別碼：A

1000-6923(2022)07-3224-10

羅成明(2000-),男,四川南充人,吉林大學(xué)碩士研究生,主要從事地下水?dāng)?shù)值模擬方面研究.

2021-11-29

國家自然科學(xué)基金資助項(xiàng)目(41972252);國家重點(diǎn)研發(fā)計(jì)劃(2018YFC1800405)

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