葛淵博，盧文喜，王梓博，王 涵，常振波

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

0 引 言

自2009年起，我國水利工程行業(yè)邁入了快速發(fā)展階段，水利工程對我國社會、經(jīng)濟(jì)的發(fā)展及環(huán)境保護(hù)等方面大有裨益［1］。但是，水利工程的建設(shè)同樣會對地下水水位及水質(zhì)造成極大的影響［2］。地下水污染一旦發(fā)生，其危害巨大且深遠(yuǎn)，地下水的污染防治迫在眉睫［3］。

針對地下水污染的預(yù)測防治，一般采用地下水?dāng)?shù)值模擬方法分析污染物在地下水中的遷移轉(zhuǎn)化規(guī)律以及分布情況。然而，由于自然因素和人為因素的雙重影響，數(shù)值模擬模型中的水文地質(zhì)參數(shù)具有高度的不確定性，從而導(dǎo)致模型輸出結(jié)果的不確定性［4］。

近年來，不確定性分析方法在地下水溶質(zhì)運(yùn)移方面發(fā)展十分迅速。蒙特卡羅模擬方法是一種國內(nèi)外常見的分析地下水?dāng)?shù)值模擬不確定性的方法，該方法適用于討論多變量模型中的不確定性問題，并且可以將數(shù)值模擬模型的參數(shù)不確定性直接轉(zhuǎn)化為模擬結(jié)果的不確定性，具有很廣泛的適用性和簡便性［5］。Hassan等［6］應(yīng)用蒙特卡羅方法評價了Alaska試驗場地下水模擬模型參數(shù)的不確定性；Wu等［7］提出將地下水?dāng)?shù)值模擬不確定性分析結(jié)果與風(fēng)險評估相結(jié)合并且取得了很好的結(jié)果；王涵［8］將不確定性分析方法與地下水DNAPLs污染多相流的隨機(jī)模擬相結(jié)合，并利用污染物濃度分布函數(shù)估算了單井遭受污染的風(fēng)險；常振波等［9］分析了參數(shù)不確定性對地下水溶質(zhì)運(yùn)移數(shù)值模擬模型輸出結(jié)果的影響，并從風(fēng)險評估的角度對不確定性分析結(jié)果進(jìn)行了解釋。

但是，進(jìn)行蒙特卡羅模擬時需要反復(fù)調(diào)用數(shù)值模擬模型，將會產(chǎn)生巨大的計算負(fù)荷。而替代模型既可以減少計算負(fù)荷，又可以保持較高的逼近精度。辛欣等［10］運(yùn)用多元回歸分析方法建立了多相流模擬模型的替代模型——雙響應(yīng)面模型，且替代模型精度較高；王梓博等［11］分別運(yùn)用克里格方法和支持向量機(jī)法建立了模擬模型的替代模型，并選擇精度較高的克里格替代模型完成了隨機(jī)模擬。但是前人很少有人嘗試將機(jī)器學(xué)習(xí)與隨機(jī)模擬相結(jié)合應(yīng)用于地下水溶質(zhì)運(yùn)移過程中。BP 神經(jīng)網(wǎng)絡(luò)作為應(yīng)用最廣泛的機(jī)器學(xué)習(xí)算法，開始在地下水領(lǐng)域內(nèi)嶄露頭角。

胡偉等［12］建立BP 神經(jīng)網(wǎng)絡(luò)對平谷平原地下水水質(zhì)進(jìn)行快速評價，且基于BP 神經(jīng)網(wǎng)絡(luò)法的評價模型擬合度較好；潘紫東等［13］采用BP 神經(jīng)網(wǎng)絡(luò)方法建立了模擬模型的替代模型，并且使用模擬退火算法求解優(yōu)化模型得到了反演識別結(jié)果。但是前人針對BP神經(jīng)網(wǎng)絡(luò)中隱含層節(jié)點(diǎn)數(shù)很少做具體研究。

本文根據(jù)一假想算例建立地下水污染數(shù)值模擬模型，運(yùn)用靈敏度分析法篩選出對模擬模型輸出結(jié)果影響較大的參數(shù)作為隨機(jī)變量，為減小反復(fù)調(diào)用模擬模型產(chǎn)生的巨大計算負(fù)荷，再運(yùn)用BP 神經(jīng)網(wǎng)絡(luò)以及克里格方法建立替代模型并比較二者精度。在運(yùn)用BP 神經(jīng)網(wǎng)絡(luò)方法時，采用三分法算法快速地確定了使替代模型誤差最小的隱含層節(jié)點(diǎn)數(shù)。然后選取精度比較高的替代模型進(jìn)行輸入輸出從而完成蒙特卡羅隨機(jī)模擬。最后，對隨機(jī)模擬的結(jié)果進(jìn)行統(tǒng)計分析與區(qū)間估計，完成地下水污染風(fēng)險評價，從而對地下水污染防治提供科學(xué)依據(jù)。

1 研究方法

1.1 蒙特卡羅法

蒙特卡羅方法是一種基于隨機(jī)數(shù)的計算方法，其主要思想是為了得到某個事件出現(xiàn)的概率或者某一個變量的期望值，可以通過做多次試驗的方法，統(tǒng)計事件出現(xiàn)的頻率或者變量的平均值，用該數(shù)字特征作為研究問題的解。在研究過程中無論隨機(jī)事件的結(jié)果滿足何種分布，在模擬次數(shù)足夠多的情況下，都可以得到一個比較精確的概率分布。

1.2 靈敏度分析法

靈敏度分析法用于分析模型中某些系數(shù)或者參數(shù)的變化對最優(yōu)解的影響。

在地下水模擬模型中，一個參數(shù)的靈敏度越大，則證明該參數(shù)相比其他參數(shù)更能影響模擬模型輸出結(jié)果的準(zhǔn)確性。

目前一般使用全局靈敏度分析和局部靈敏度分析兩種分析方法，其中，局部靈敏度分析的計算公式如下［14］：

式中：Xk為靈敏度系數(shù)。計算Xk也可使用如下公式：

也可采取標(biāo)準(zhǔn)化無量綱形式：

本文將會采用局部靈敏度分析法評價單一參數(shù)的變化對模型輸出結(jié)果造成的影響。

1.3 替代模型

進(jìn)行蒙特卡羅模擬時需要大量的統(tǒng)計試驗來完成模擬，如果直接調(diào)用數(shù)值模擬模型進(jìn)行統(tǒng)計實驗，會帶來巨大的計算負(fù)荷。而替代模型可以用來代替數(shù)值模擬模型的輸入輸出響應(yīng)關(guān)系，能夠以更小的計算量得到與模擬模型相近的結(jié)果［15］。

1.3.1 基于BP算法的三層神經(jīng)網(wǎng)絡(luò)

1968年由Rumehart 和Mcclelland 為首的科學(xué)家們首次提出了BP 神經(jīng)網(wǎng)絡(luò)這一科學(xué)概念［16］，這是一種基于誤差反向傳播算法進(jìn)行訓(xùn)練的多層前饋神經(jīng)網(wǎng)絡(luò)。單隱層的BP 神經(jīng)網(wǎng)絡(luò)具有較高的訓(xùn)練效率，可以逼近任何一個閉區(qū)間內(nèi)的連續(xù)函數(shù)，完成任意大小的n維到L維的映射。

經(jīng)典的BP 神經(jīng)網(wǎng)絡(luò)通常由三層組成：輸入層，隱含層與輸出層。計算過程分為兩個階段：前饋過程（Feed-Forward）與反向傳播過程（BackPropagation）［17］。

（1）前饋過程：每個隱含層和輸出層神經(jīng)元輸出與輸入的函數(shù)關(guān)系為：

圖1 BP神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)Fig.1 BP neural network structure

式中：Wij表示神經(jīng)元i與神經(jīng)元j之間連接的權(quán)重；Oj代表神經(jīng)元j的輸出；sigmod 函數(shù)稱為神經(jīng)元的激勵函數(shù)（activation function），除了sigmod函數(shù)外，常用還有tanh和ReLU函數(shù)。

將n個輸入向量依次送入輸入神經(jīng)元，隱含層神經(jīng)元獲得輸入層的輸出并計算出輸出值，最后輸出層的神經(jīng)元根據(jù)隱含層輸出計算出回歸值。

（2）反向傳播過程：根據(jù)誤差調(diào)整連接權(quán)重Wij的過程稱之為反向傳播過程。首先隨機(jī)初始化連接權(quán)重Wij，對某一訓(xùn)練樣本進(jìn)行一次前饋過程得到各神經(jīng)元的輸出。再計算輸出層的誤差：

式中：Ej代表神經(jīng)元j的誤差；Oj表示神經(jīng)元j的輸出；Tj表示當(dāng)前訓(xùn)練樣本的參考輸出；sigmod'(x)是上文sigmod 函數(shù)的一階導(dǎo)數(shù)。

1.3.2 克里格（Kriging）方法

克里格方法將輸入與輸出變量之間的關(guān)系以回歸方程的形式表達(dá)出來，是一種常用的地質(zhì)學(xué)統(tǒng)計方法，并且是一種黑箱模型。

其公式如下［18］：

式中：R(xi，xj)代表xi和xj的關(guān)聯(lián)函數(shù)。

高斯模型的公式如下所示：

式中：θk可以通過算法求出。

1.4 三分法算法

在BP 神經(jīng)網(wǎng)絡(luò)中，隱含層節(jié)點(diǎn)數(shù)是一個非常重要的參數(shù)，若隱含層的節(jié)點(diǎn)數(shù)太少，網(wǎng)絡(luò)將無法正確學(xué)習(xí)；若節(jié)點(diǎn)數(shù)太多，盡管在一定程度上可以減小系統(tǒng)誤差，但網(wǎng)絡(luò)學(xué)習(xí)的時間也會隨之加長，而且會使訓(xùn)練掉入局部極小點(diǎn)，產(chǎn)生過擬合現(xiàn)象［19］。所以說隱含層節(jié)點(diǎn)數(shù)的合理設(shè)計對BP 神經(jīng)網(wǎng)絡(luò)的性能有很大影響。在本文中根據(jù)前人經(jīng)驗使用三分法算法［20］確定隱含層節(jié)點(diǎn)數(shù)，最終得到使BP 神經(jīng)網(wǎng)絡(luò)擬合效果最好的、精度最高的隱含層節(jié)點(diǎn)數(shù)，進(jìn)而使用BP神經(jīng)網(wǎng)絡(luò)建立替代模型進(jìn)行求解。

“三分法”算法流程如下所示：

式中：k為樣本量；n1為隱層節(jié)點(diǎn)數(shù)；n為輸入層節(jié)點(diǎn)數(shù)。當(dāng)i＞n1時= 0。

式中：n為輸入層節(jié)點(diǎn)數(shù)；m 為輸出層節(jié)點(diǎn)數(shù)；a為1～10 之間的常數(shù)。

根據(jù)以上經(jīng)驗公式可得隱含層節(jié)點(diǎn)數(shù)的取值為［a，b］，且根據(jù)二次函數(shù)的圖像性質(zhì)，當(dāng)最低點(diǎn)位于區(qū)間內(nèi)時，那么該區(qū)間的斜率絕對值最小，并且左邊區(qū)間斜率為負(fù)，右邊區(qū)間斜率為正。

以下為三分法算法流程圖見圖2。

圖2 三分法算法流程圖Fig.2 Flowchart of the trinity algorithm

2 案例應(yīng)用

2.1 問題概述

本文針對假想算例展開研究：研究區(qū)長約544 m，南北寬約555 m，面積約為0.302 km2。研究的主要對象為一層潛水含水層，將研究區(qū)概化為非均質(zhì)各向同性含水層，水流為二維非穩(wěn)定流。水流方向由北至南，區(qū)域內(nèi)地表平坦，北側(cè)Γ1及南側(cè)Γ2均概化為定水頭邊界，東側(cè)Γ3及西側(cè)Γ4均概化為零流量邊界。研究區(qū)內(nèi)地下水補(bǔ)給項主要為降水入滲補(bǔ)給，年平均降水量約為800 mm，降水入滲系數(shù)為0.4。地下水排泄方式主要是人工開采，區(qū)內(nèi)設(shè)置三口抽水井，分別是井1、2、3，且三口抽水井的抽水量均為500 m3/d，并將其作為觀測井。區(qū)內(nèi)有一固體廢棄物堆放，經(jīng)過降水淋濾，向地下水排放污染物量為5 kg/d，假設(shè)污染物不會發(fā)生化學(xué)變化以及生物轉(zhuǎn)化與遷移。潛水含水層中污染物初始濃度為0。將研究區(qū)按照潛水含水層滲透系數(shù)、縱向彌散度以及給水度劃分為南北兩個區(qū)域（圖3 中東西向?qū)嵕€為南北兩地區(qū)的分界線）。研究區(qū)示意圖及各水文地質(zhì)參數(shù)表如圖3所示。

圖3 研究區(qū)示意圖Fig.3 Schematic diagram of the study area

表1 研究區(qū)水文地質(zhì)參數(shù)數(shù)值Tab.1 Values of hydrogeological parameters in the study area

2.2 模型建立

2.2.1 數(shù)值模擬

（1）地下水水流模型。

式中：K為含水層的滲透系數(shù)，m/d；H為地下水水位，m；B為隔水底板高程，m；w為源匯項；μ為含水層給水度；Γ1，Γ2為第一類邊界，Γ3，Γ4為第二類邊界為邊界上某點(diǎn)（x，y）處外法線方向上的單位向量；φ(x，y，t)為已知水位函數(shù)。

（2）地下水溶質(zhì)運(yùn)移模型。

式中：t為時間；c為溶質(zhì)的濃度，mg/L；Dx、Dy為水動力彌散系數(shù)，L2/T；ux為實際平均流速向量u在x方向上的分量；uy為實際平均流速向量u在y方向上的向量；n為含水層孔隙度；b為含水層厚度；I為源匯項；Γ1為已知濃度邊界；Γ2為已知對流-彌散通量邊界；Γ3，Γ4為零通量的水動力彌散通量邊界；c0(x，y)，c1(x，y，t)，c2(x，y，t) 是已知濃度函數(shù)。

在已知研究區(qū)概況以及污染質(zhì)運(yùn)移數(shù)值模型基礎(chǔ)之上，利用GMS 軟件將研究區(qū)剖分為2886 個有限差分網(wǎng)格，然后利用MODFLOW 工具箱以及MT3DMS 工具箱模擬10年后研究區(qū)污染質(zhì)運(yùn)移情況。由于本文算例中固體廢棄物的滲濾液以固定速率排出，所以僅需一個應(yīng)力期并劃分了10 個時間步長。10年后污染物濃度分布如圖4所示。

圖4 10年后污染物濃度分布圖Fig.4 Pollutant concentration distribution after 10 years

2.2.2 靈敏度分析法篩選隨機(jī)變量

靈敏度分析法用來篩選出對數(shù)值模擬模型輸出結(jié)果影響較大的參數(shù)，并將其作為模型中的隨機(jī)變量。

本文選取滲透系數(shù)K1、K2，給水度μ1、μ2，縱向彌散度α1、α2以及孔隙度n七個水文地質(zhì)參數(shù)進(jìn)行靈敏度分析：

首先將各參數(shù)在各自的數(shù)值參考范圍內(nèi)取均值輸入模擬模型，得到第一組輸出，再將所有水文地質(zhì)參數(shù)的均值分別增加、減少10%和20%，在保證只改變其中一個參數(shù)的情況下輸入到模擬模型中，得到第二組輸出結(jié)果，運(yùn)用公式計算三口觀測井各參數(shù)的靈敏度系數(shù)值，見表2。

表2 靈敏度近似計算解Tab.2 Sensitivity approximate calculation solution

由圖5可知，兩個分區(qū)的滲透系數(shù)K1、K2以及兩個分區(qū)的縱向彌散度α1、α2的靈敏度系數(shù)相比其他3 個水文地質(zhì)參數(shù)的靈敏度系數(shù)較大。接下來將滲透系數(shù)K1、K2，縱向彌散度α1、α2作為隨機(jī)變量，其他它參數(shù)作為確定性變量。

圖5 靈敏度分析結(jié)果圖Fig.5 Results of sensitivity analysis

2.2.3 拉丁超立方抽樣

本文選擇拉丁超立方抽樣方法對三口觀測井內(nèi)的4個隨機(jī)變量分別抽樣取值并隨機(jī)組合。根據(jù)前人經(jīng)驗總結(jié)，各水文地質(zhì)參數(shù)的取值及服從分布如表3所示。

表3 參數(shù)概率分布及取值情況Tab.3 Probability distribution and value of parameters

運(yùn)用MATLAB軟件對上述參數(shù)分別抽樣60組，并將其隨機(jī)組合，再運(yùn)用數(shù)值模擬模型求解，得到井1、2、3 的60 組污染質(zhì)濃度輸出結(jié)果。

2.2.4 建立替代模型

將上述四個水文地質(zhì)參數(shù)輸入地下水模擬模型中，計算得到相應(yīng)的60 組輸出。再將前50 組輸入輸出作為訓(xùn)練樣本，后10 組輸入輸出作為檢驗樣本分別利用BP 神經(jīng)網(wǎng)絡(luò)方法及克里格方法建立替代模型并檢驗替代模型的精度。

本文將基于以下3個評價指標(biāo)比較兩個替代模型的精度：

（1）均方誤差（MSE）。均方誤差表示預(yù)測值與真實值之差平方的期望值，MSE可以用來評價數(shù)據(jù)的變化程度，計算公式如下所示：

MSE的值越小，說明替代模型擁有更好的精確度。

（2）確定性系數(shù)R2（R-Square）。計算公式如下所示：

R2的取值可以用來判斷替代模型的好壞，其取值范圍為［0，1］，取值越接近1，則擬合效果越優(yōu)。

（3）平均相對誤差（MAPE）。MAPE是一個百分比值，表示預(yù)測結(jié)果較真實結(jié)果的平均偏離，計算公式如下所示：

MAPE的值越小，說明替代模型擁有更好的精確度。

上式中，m為樣本個數(shù)，yi為模擬模型的輸出值為替代模型的輸出值為模擬模型輸出值的平均值。

首先基于BP 神經(jīng)網(wǎng)絡(luò)建立數(shù)值模擬模型的替代模型，將研究區(qū)模型中滲透系數(shù)K1、K2，縱向彌散度α1、α2作為輸入向量輸入，將3口監(jiān)測井得到的污染物濃度作為輸出向量，建立了一個三層BP神經(jīng)網(wǎng)絡(luò)。

三分法算法程序運(yùn)行過程中a，b，c，d（m，n均設(shè)置為0）的變化是：

（1）a=3，c=6，d=9，b=12；

（2）a=6，b=9。

根據(jù)三分法算法最終確定隱含層節(jié)點(diǎn)數(shù)為6，并利用matlab 軟件建立了如圖6 所示的BP 神經(jīng)網(wǎng)絡(luò)替代模型，其中Hidden Layer 以及Output Layer 中的W與b分別代表BP 神經(jīng)網(wǎng)絡(luò)中的參數(shù)權(quán)重與神經(jīng)元偏置。

圖6 BP神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)圖Fig.6 BP neural network structure diagram

再利用克里格方法建立替代模型，根據(jù)檢驗組的輸入輸出，最終得到兩種替代模型各自的評估指標(biāo)，如表4所示。

表4 替代模型精度分析Tab.4 Precision analysis of substitution model

可以明顯看出，基于BP 神經(jīng)網(wǎng)絡(luò)建立的替代模型的精度高于基于克里格方法建立的替代模型。所以本文使用BP 神經(jīng)網(wǎng)絡(luò)替代模型進(jìn)行隨機(jī)模擬。

2.3 蒙特卡羅隨機(jī)模擬

本文利用蒙特卡羅方法對地下水污染進(jìn)行隨機(jī)模擬。首先采用拉丁超立方抽樣方法抽取上述1 000 組靈敏度較高的水文地質(zhì)參數(shù)，然后將其輸入到BP 神經(jīng)網(wǎng)絡(luò)替代模型中，得到相應(yīng)的1、2、3號井的1 000組濃度輸出，最后對這三口井進(jìn)行統(tǒng)計分析。

3 結(jié)果與討論

3.1 統(tǒng)計分析

本文利用SPSS 軟件中的K-S 檢驗方法分別對三口觀測井的污染物濃度進(jìn)行分布檢驗，包括正態(tài)分布、指數(shù)分布、均勻分布以及泊松分布4種，其中正態(tài)分布檢驗結(jié)果如表5所示。

由表5 可以看出三口井污染物濃度的檢驗統(tǒng)計分別為0.019、0.023、0.018，對應(yīng)的顯著性水平均大于0.05，即K-S 檢驗接受零假設(shè)，服從指定分布，即三口井的污染物濃度均符合正態(tài)分布。

表5 單樣本柯爾莫戈洛夫-斯米諾夫檢驗結(jié)果表 mg/LTab.5 Results of Kolmogorov-Sminov test for a single sample

分別繪制三口井的污染物濃度頻數(shù)分布直方圖見圖7。

由圖7 可知，一號井污染物濃度在370～410 mg/L 時的概率最大，為72%，二號井污染物濃度在205～225 mg/L 時的概率最大，為68%，三號井污染物濃度在120～130 mg/L 時的概率最大，為65%。

圖7 三口觀測井濃度頻數(shù)分布直方圖Fig.7 Histogram of concentration frequency distribution of three observation Wells

利用SPSS 軟件得到如表6 所示三口觀測井的各項統(tǒng)計指標(biāo)。

表6 三口觀測井污染物濃度統(tǒng)計Tab.6 Pollutant concentration statistics of the three observation Wells

由表6可以看出，三口觀測井的變異系數(shù)均小于15%，可認(rèn)為三口觀測井污染物濃度輸出的結(jié)果離散性較小，不確定性也很小。

3.2 區(qū)間估計

區(qū)間估計是參數(shù)估計的一種形式，指通過從總體中抽取樣本，根據(jù)一定的正確度要求，構(gòu)造出適當(dāng)?shù)膮^(qū)間以作為總體分布參數(shù)的真值所在區(qū)間的估計。本文將其引申，以三口觀測井的濃度輸出值作為要觀測的參數(shù)，計算其在不同置信水平下的區(qū)間范圍。

已知三口觀測井的污染物濃度輸出近似服從正態(tài)分布，故利用z分布進(jìn)行區(qū)間估計。由表7 可以看出，置信水平越高，區(qū)間范圍越大；置信水平越低，區(qū)間范圍越小，且越集中在平均值附近。

表7 三口觀測井污染物濃度區(qū)間估計Tab.7 Estimation of pollutant concentration interval for three observation Wells

3.3 風(fēng)險評估

根據(jù)替代模型的輸出結(jié)果，本文繪制出分布函數(shù)曲線圖對三口井的污染物濃度輸出結(jié)果進(jìn)行分析。

根據(jù)三口井的不同用途，假設(shè)當(dāng)三口觀測井的濃度輸出結(jié)果同時滿足一號井污染物濃度大于400 mg/L，二號井污染物濃度大于220 mg/L，三號井污染物濃度大于127 mg/L 時，證明研究區(qū)地下水已受到嚴(yán)重污染。根據(jù)累積分布函數(shù)曲線圖可以看出，一號井污染物濃度大于400 mg/L 的風(fēng)險是0.66，二號井污染物濃度大于220 mg/L 的風(fēng)險是0.60，三號井污染物濃度大于127 mg/L的風(fēng)險是0.58。

圖8 三口觀測井污染物濃度分布函數(shù)圖Fig.8 Distribution function of pollutant concentration in three

4 結(jié) 論

（1）本文運(yùn)用靈敏度分析法在模擬模型的眾多參數(shù)中篩選出了對模擬模型輸出結(jié)果影響較大的兩個參數(shù)，極大程度上降低了替代模型的維數(shù)并減少了計算負(fù)荷。

（2）本文運(yùn)用三分法算法很快地確定出了使BP 神經(jīng)網(wǎng)絡(luò)精度最高的隱含層節(jié)點(diǎn)數(shù)。

（3）本文分別采用克里格方法以及BP 神經(jīng)網(wǎng)絡(luò)方法建立了替代模型，且確定了最優(yōu)隱含層節(jié)點(diǎn)數(shù)的BP 神經(jīng)網(wǎng)絡(luò)替代模型具有更高的逼近精度。

（4）本文利用蒙特卡羅方法對地下水?dāng)?shù)值模擬模型中的參數(shù)進(jìn)行不確定性分析，并根據(jù)隨機(jī)模擬的結(jié)果完成了地下水污染風(fēng)險評價，為地下水污染的防治提供了更加合理的依據(jù)。 □