邢貞相，曲睿卓，趙 瑩，紀(jì) 毅，張 涵

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

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

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

1 研究方法

1.1 克里格（Kriging）

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

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

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

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

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

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

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

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

1.3 最小二乘支持向量機(jī)（LSSVM）

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

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

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

根據(jù)（Karush-Kuhn-Tucker）KKT條件，LSSVM方程兩邊偏導(dǎo)數(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ù)具有廣泛收斂域和強(qiáng)泛化能力，為理想回歸核函數(shù)?？杀硎救缦拢?/p>

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

1.4 評(píng)價(jià)指標(biāo)

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

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

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

②相對(duì)誤差（Relativeerror,RE）

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

2.1 替代模型建立

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

2.2 替代模型驗(yàn)證

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

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

2.3 地下水污染源識(shí)別模型建立

完成Kriging、LSSVM和GRNN替代模型建立后，將其替代模擬模型與優(yōu)化模型結(jié)合，反演地下水污染源釋放過(guò)程。運(yùn)用Matlab中GA優(yōu)化工具箱求解未知參數(shù)，種群數(shù)為20，迭代200次達(dá)到收斂。當(dāng)觀測(cè)數(shù)據(jù)存在噪聲時(shí)，優(yōu)化過(guò)程增加至迭代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í)驗(yàn)室研究均可能存在測(cè)量誤差，觀察濃度具有不確定性。因此，實(shí)際濃度如下式：

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

3 案例分析

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

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

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

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

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

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

圖2 不同替代模型對(duì)案例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ī)則形狀，含水層不均勻，平面圖見(jiàn)圖3。該含水層有5個(gè)不同水力傳導(dǎo)帶，每個(gè)區(qū)域水力傳導(dǎo)率值恒定。水力傳導(dǎo)率值為：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個(gè)污染源和7個(gè)取樣點(diǎn)。模擬時(shí)間為10年（120個(gè)月），且分為20個(gè)相等壓力期，每個(gè)應(yīng)力期持續(xù)6個(gè)月。假設(shè)源在前4個(gè)應(yīng)力期釋放污染物。

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

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

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

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

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

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

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

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

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

圖6 不同替代模型對(duì)案例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é)果誤差增加。在無(wú)噪聲情況下，運(yùn)用GRNN替代模型反演精度較低，反演結(jié)果較合理。當(dāng)觀測(cè)值存在噪聲情況下，即使在低噪聲水平下，反演結(jié)果RE高達(dá)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，說(shuō)明在中等噪聲水平下，運(yùn)用Kriging替代模型反演精度良好。此外，由圖7可知，即使在高水平噪聲下，Kriging替代模反演精度被接受。

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

4 結(jié)論

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

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

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

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