不同替代模型在地下水污染源釋放歷史反演中適用性研究

邢貞相，曲睿卓，趙 瑩，紀 毅，張 涵

（1.東北農(nóng)業(yè)大學水利與土木工程學院，哈爾濱 150030；2.農(nóng)業(yè)部農(nóng)業(yè)水資源高效利用重點實驗室，哈爾濱 150030）

隨著城鎮(zhèn)化建設和經(jīng)濟快速發(fā)展，地下水污染日益嚴重。地下水系統(tǒng)復雜，污染源信息不易獲取[1]，難以在源頭治理地下水污染，治理成本高、難度大[2]，治理修復措施無法完全發(fā)揮作用[3]。因此，降低地下水污染風險，識別地下水污染源特征成為當前研究重點。Bagtzoglou和Atmadja運用直接水文反演方法，重建保守污染羽空間分布[4]。Alapati和Kabala使用非線性最小二乘法識別地下水污染物羽流釋放歷史[5]。優(yōu)化算法是解決地下水反演問題最普遍且高效的方法[6-7]。優(yōu)化算法分為局部優(yōu)化算法[8]和啟發(fā)式算法[9]。相對于局部優(yōu)化算法，啟發(fā)式算法不會陷入局部極小值，可解決非凸模型優(yōu)化問題，廣泛應用于地下水污染源識別[10]。遺傳算法（Genetic algorithm，GA）是最常見啟發(fā)式算法，應用于解決二維含水層中污染源反演問題，可提高模型計算效率[11]。利用GMS軟件中MODFLOW和MT3D兩個模塊可建立地下水污染運移過程模擬模型，將優(yōu)化模型與模擬模型結(jié)合，通過反演求得污染源釋放過程。遺傳算法可提高計算效率，但反演過程反復調(diào)用模擬模型，耗時較長，適用性較弱。替代模型計算時間短，功能與模擬模型相似。替代模型研究中，曾主要以人工神經(jīng)網(wǎng)絡（Artificial neural network,ANN）為主[12]。Zhao等研究發(fā)現(xiàn)克里格（Kriging）模型在相同地質(zhì)條件下，反演速度和準確性均優(yōu)于ANN[13]。其他替代模型包括雙響應面法[14]、BP神經(jīng)網(wǎng)絡模型[15]、徑向基函數(shù)神經(jīng)網(wǎng)絡（RBF）[16]、支持向量機法（SVM）[17]等。

替代模型是解決地下水污染源反演有效方法，但不同替代模型對水文地質(zhì)條件區(qū)域適用性不同。因此，本文選用Kriging、廣義回歸神經(jīng)網(wǎng)絡（GRNN）和最小二乘支持向量機法（LSSVM）3種應用較廣泛替代模型，對比分析不同復雜水文地質(zhì)條件下其對地下水污染源釋放歷史反演適用性。目前Kriging和RBF替代模型較廣泛。但RBF替代模型存在難以確定隱層節(jié)點等問題且適用性一般，故選用GRNN替代模型。GRNN替代模型在學習速率與逼近能力方面明顯優(yōu)于RBF替代模型[18]。選擇LSSVM作為替代模型速度更快、精確度更高[19]。為探究3種替代模型適用性，將選擇相同輸入輸出訓練模型，從模型學習能力和反演時間、準確性和穩(wěn)定性等方面評價。

1 研究方法

1.1 克里格（Kriging）

Kriging模型是依據(jù)協(xié)方差函數(shù)對隨機過程空間建模和預測回歸的算法[20]。其基本形式為：

其中x代表輸入數(shù)據(jù)，y（x）代表模型輸出，g（x）代表原始函數(shù)全局模型，Z（x）代表高斯隨機函數(shù)。Z（x）協(xié)方差可計算如下：

其中R代表相關函數(shù)，r代表對稱相關矩陣，n代表樣本數(shù)。

1.2 廣義回歸神經(jīng)網(wǎng)絡（GRNN）

廣義回歸神經(jīng)網(wǎng)絡具有很強非線性擬合能力和柔性網(wǎng)絡結(jié)構(gòu)[21]，GRNN模型由輸入層、模式層、求和層和輸出層組成。輸入層神經(jīng)元數(shù)目等于學習樣本中輸入向量維數(shù)。模式層為隱含層，其神經(jīng)元數(shù)目等于樣本數(shù)n，神經(jīng)元i傳遞函數(shù)為：

式中，xi為各單元核函數(shù)中心向量。

在求和層中，計算模式層各單元輸出之和與模式層各單元輸出加權，其傳遞函數(shù)分別為：

式中，yi為各訓練樣本加權和權重值；最后由輸出層計算輸出，輸出層神經(jīng)元數(shù)目等于樣本中輸出向量維數(shù)[22]，公式為：

1.3 最小二乘支持向量機（LSSVM）

LSSVM用于解決模式分類和函數(shù)估計問題，通過計算樣本最小平方誤差擬合。其功能描述如下[23]：

訓練樣本集是{( xi,yi)|i=1,2,…,n} ，xi是輸入變量值，yi是對應輸出變量。其中x是輸入向量，ω是權重向量，d是偏差，φ是高維特征空間映射。根據(jù)結(jié)構(gòu)風險最小化原則，尋求ω最小化方程為：

其中γ是懲罰因子，ei是誤差控制函數(shù)。將約束條件在方程式（8）中給出：

根據(jù)（Karush-Kuhn-Tucker）KKT條件，LSSVM方程兩邊偏導數(shù)為：

式中，y=[y1,y2, …,yn]T，a=[a1,a2, …,an]T，lv=[1,1，...，1]T，I為單位矩陣，g為相鄰矩陣。g（xi,xj）=φ（xi）Tφ（xj）,i,j=1,…,n,為滿足Mercer理論核矩陣[24]。徑向基函數(shù)具有廣泛收斂域和強泛化能力，為理想回歸核函數(shù)?？杀硎救缦拢?/p>

式中，σ-核寬度參數(shù)，影響LSSVM在特征空間中獲得最優(yōu)分類超平面泛化能力。LSSVM預測模型可表示為：

1.4 評價指標

為驗證替代模型精度，采用確定性系數(shù)（R2）、相對誤差（RE）和均方根誤差（RMSE）3種精度評估指標。R2越接近1，表示模型擬合精度越高。RE和RMSE越接近0，表示模型精度擬合精度越高。其公式如下：

①確定性系數(shù)（Coefficient of determination,R2）

式中，yi為模擬模型輸出值，即真實觀測濃度值；y?i為替代模型輸出值；yˉi為模擬模型輸出值平均值。

②相對誤差（Relativeerror,RE）

2 地下水流數(shù)值模擬模型建立

2.1 替代模型建立

利用拉丁超立方抽樣方法對模型輸入抽樣，即為地下水污染源釋放濃度。本文研究3種不同水文地質(zhì)條件案例分別有500組輸入數(shù)據(jù)。分別輸入GMS模擬模型中，得到500組輸出值，即得到觀測井觀測濃度值。將其中相同450組輸入和輸出作為3個替代模型訓練樣本，利用剩余50組輸入和輸出數(shù)據(jù)作為檢驗組驗證替代模型精度。

2.2 替代模型驗證

3種替代模型不同情況下模擬表現(xiàn)存在差別，不同案例情況下3種替代模型得到訓練，驗證模型精度，分析模型學習能力。結(jié)果見表1和2。

確定性系數(shù)反映模型預測值與實際觀測值擬合度，其可檢驗替代模型模擬能力，但無法檢驗預測極值能力。而相對誤差對模擬值中極大或極小值誤差反映敏感。通過模擬3個案例觀測值發(fā)現(xiàn)，3種替代模型R2均在0.9以上，RE均在20%以內(nèi)。說明建立這3種替代模型有效，可運用于后續(xù)反演。對比發(fā)現(xiàn)，GRNN替代模型學習能力較弱，但表現(xiàn)較穩(wěn)定，在案例1和案例2條件下R2和RE較接近。在案例3情況下，模擬精度下降。Kriging替代模型學習能力較強，R2均在0.99以上，但對于極小值和位置較遠觀測井觀測值模擬能力較弱。在案例3中，第6、7觀測井位置RE高于其他口井位置。LSSVM替代模型在水文地質(zhì)較為簡單情況下，預測能力較強，隨水文地質(zhì)條件變化，R2降低，RE增高。通過案例1和案例3對比可知，LSSVM在復雜水文地質(zhì)條件適用性較差。

2.3 地下水污染源識別模型建立

完成Kriging、LSSVM和GRNN替代模型建立后，將其替代模擬模型與優(yōu)化模型結(jié)合，反演地下水污染源釋放過程。運用Matlab中GA優(yōu)化工具箱求解未知參數(shù)，種群數(shù)為20，迭代200次達到收斂。當觀測數(shù)據(jù)存在噪聲時，優(yōu)化過程增加至迭代500次。

表1 替代模型在案例1中模擬結(jié)果Table1 Simulation resultsof surrogatemodelsin case 1

表2 替代模型在案例2和案例3中模擬結(jié)果Table 2 Simulation resultsof surrogate modelsin case 2 and case3

2.4 測量誤差

由于野外和實驗室研究均可能存在測量誤差，觀察濃度具有不確定性。因此，實際濃度如下式：

其中C為觀測污染物濃度。ε是誤差矩陣，表示數(shù)據(jù)中存在噪聲水平。如果a＜0.10，則噪聲水平低。如果0.10≤a≤0.15，則噪聲水平適中。如果a>0.15，則噪聲水平高。

3 案例分析

3.1 二維規(guī)則邊界穩(wěn)流均質(zhì)含水層（案例1）

本研究通過3種不同水文地質(zhì)條件案例，測試反演模型性能。第1種水文地質(zhì)研究區(qū)域平面圖見圖1。

圖1 案例1含水層結(jié)構(gòu)平面Fig.1 Aquifer structureplane of case 1

含水層模型邊界條件是東側(cè)和西側(cè)特定水頭線性變化，北側(cè)和南側(cè)無流動。該含水層有3個污染源和4個采樣點，每個采樣點采集污染濃度。模擬時間為5年（60個月），分為20個相等壓力期，每個壓力期持續(xù)3個月。假設源在前4個壓力期間釋放污染物。表3列出實際污染源釋放濃度。根據(jù)上述建立模擬-優(yōu)化模型對案例1作地下水污染源釋放過程反演，得到結(jié)果見圖2，誤差分析結(jié)果見表4。

由圖2可知，3種替代模型反演結(jié)果與真實地下水污染源釋放濃度較接近。R2均在0.99以上，RE均在5%以內(nèi)，RMSE均為1.5。說明案例1水文地質(zhì)情況下，3種替代模型反演精度較高，適用性較好。

表3 案例1地下水污染源流量實際值Table 3 Actual value of groundwater pollution sourceflux in Case 1

圖2 不同替代模型對案例1污染源釋放歷史反演結(jié)果Fig.2 Retrieval results of different surrogatemodels on the release history of case 1 pollution sources

表4 案例1反演結(jié)果誤差Table 4 Retrieval result error in case 1

3.2 二維非規(guī)則邊界非均質(zhì)穩(wěn)流含水層（案例2）

案例2研究區(qū)域是幾何不規(guī)則形狀，含水層不均勻，平面圖見圖3。該含水層有5個不同水力傳導帶，每個區(qū)域水力傳導率值恒定。水力傳導率值為：K1=0.0004 m·s-1，K2=0.0002 m·s-1，K3=0.0001 m·s-1，K4=0.0003 m·s-1，K5=0.0007 m·s-1s。

該含水層有2個污染源和7個取樣點。模擬時間為10年（120個月），且分為20個相等壓力期，每個應力期持續(xù)6個月。假設源在前4個應力期釋放污染物。

表5列出實際污染源釋放濃度。反演結(jié)果見圖4，誤差分析結(jié)果見表6。

圖3 案例2含水層結(jié)構(gòu)平面Fig.3 Aquifer structure plane of case 2

表5 案例2地下水污染源流量實際值Table 5 Actual value of groundwater pollution source flow in case 2

圖4 不同替代模型對案例2污染源釋放歷史反演結(jié)果Fig.4 Retrieval results of different surrogate models on the release history of case 2 pollution sources

表6 案例2反演結(jié)果誤差Table6 Retrieval result error in case 2

由表6可知，在案例2較復雜含水層條件下，Kriging和LSSVM替代模型表現(xiàn)較穩(wěn)定。與案例1相比，Kriging和LSSVM替代模型R2下降0.02%，RMSE升至2.20，RE均小于10%，表明兩種模型對案例2含水層條件適用性較好。而GRNN替代模型RE和RMSE均較高，且R2下降0.07%。雖然反演結(jié)果在合理范圍內(nèi)，但相較于其他兩種模型，反演結(jié)果誤差較大。

3.3 二維非規(guī)則邊界非穩(wěn)流非均質(zhì)含水層（案例3）

由圖5可知，含水層其他特征與第2種情況相同。與案例2不同的是，1個抽水井存在于地質(zhì)層中心，模擬地下水水流瞬態(tài)流動和瞬態(tài)傳輸條件。所需反演實際值與案例2相同。反演結(jié)果如圖6所示，誤差分析結(jié)果如表7所示。

由圖6和表7可知，Kriging替代模型反演結(jié)果較其他2種替代模型更優(yōu)，R2接近0.99，RE在10%以內(nèi)且RMSE為2.8。運用LSSVM替代模型反演結(jié)果不如Kriging替代模型，但相較于案例2，表現(xiàn)較穩(wěn)定。其RE在10%以內(nèi)，RMSE在5以內(nèi)，說明在案例3情況下，LSSVM替代模型適用性較好。運用GRNN替代模型反演結(jié)果精度較低，R2小于0.9，RMSE為6.62，RE高達20%以上，誤差較大。說明GRNN替代模型對于案例3水文地質(zhì)情況適用性較差。

為進一步檢驗3種替代模型結(jié)合優(yōu)化模型對地下水污染源釋放歷史反演適用性和穩(wěn)定性，反演結(jié)果和誤差分析見圖7。

圖5 案例3含水層結(jié)構(gòu)平面Fig.5 Aquifer structureplaneof case 3

圖6 不同替代模型對案例3污染源釋放歷史反演結(jié)果Fig.6 Retrieval results of different surrogate models on the release history of case 3 pollution sources

表7 案例3反演結(jié)果誤差Table7 Retrieval result error in case 3

圖7 在不同噪聲水平下替代模型反演結(jié)果誤差分析Fig.7 Error analysis of the retrieval results of surrogate models at different noiselevels

隨噪聲水平增加，反演結(jié)果誤差增加。在無噪聲情況下，運用GRNN替代模型反演精度較低，反演結(jié)果較合理。當觀測值存在噪聲情況下，即使在低噪聲水平下，反演結(jié)果RE高達26.88。在高水平噪聲條件下，R2降至0.7，反演結(jié)果較差。而LSSVM替代模型，隨噪聲水平增加，反演結(jié)果誤差變大，但該模型表現(xiàn)穩(wěn)定，誤差值上升平穩(wěn)。在中等噪聲水平或以下，反演結(jié)果和Kriging替代模型反演結(jié)果接近。與前兩種替代模型相比，Kriging替代模型在噪聲水平a=10%情況下，R2在0.96以上，RE為10%，RMSE小于5，說明在中等噪聲水平下，運用Kriging替代模型反演精度良好。此外，由圖7可知，即使在高水平噪聲下，Kriging替代模反演精度被接受。

綜合以上3個案例反演結(jié)果發(fā)現(xiàn)，Kriging替代模型適用于所有案例，觀測值存在較高水平噪聲時可獲得精度較高反演結(jié)果。LSSVM替代模型適用于所有案例，但在觀測值噪聲水平較高情況下適用性次于Kriging替代模型。而GRNN替代模型僅在案例1和案例2水文地質(zhì)條件下適用性優(yōu)良，在案例3中反演精度較低。尤其在觀測值存在噪聲情況下，反演結(jié)果可能不被接受。提出替代模型主要目的是提高計算效率。在3個替代模型反演過程中，分別記錄反演所需時間。對比可知，運用Kriging替代模型反演過程所需計算時間是LSSVM和GRNN反演所需時間2倍。

4 結(jié)論

本文分別運用Kriging、LSSVM和GRNN建立替代模型，結(jié)合遺傳算法反演3種不同含水層條件和多個噪聲水平下地下水污染源釋放過程。重點討論相同情況下3種替代模型反演精度和各單一替代模型在不同復雜情況適用性，研究成果可為解決案例相似條件下地下水污染源釋放歷史反演現(xiàn)實問題中替代模型有效選擇提供試驗支撐。本文在研究過程中設定3種案例中監(jiān)測污染物穩(wěn)定不可降解，但監(jiān)測污染物可被降解和吸附，后續(xù)將研究不穩(wěn)定污染物污染源反演。

a.GRNN替代模型僅在簡單水文地質(zhì)情況下適用性較好，其他較復雜情況下反演結(jié)果較差，尤其在非穩(wěn)定流地下水含水層，但反演時間短。

b.LSSVM替代模型雖在復雜水文地質(zhì)情況適用性優(yōu)良且反演時間短，但當觀測值存在較大觀測誤差時，反演精度下降，適用性較差。

c.在3種替代模型中Kriging替代模型最為穩(wěn)定，在復雜水文地質(zhì)情況適用性更強，更適用于觀測值存在噪聲情況。但Kriging替代模型反演所需時間成本較高，現(xiàn)實應用應根據(jù)實際情況選擇替代模型。