基于DRASTIC- LU參數(shù)和數(shù)據(jù)驅動模型的地下水硝酸鹽脆弱性分區(qū)

Seyed Ahmad Eslaminezhad, Mobin Eftekhari, Mohammad Akbari

(1.德黑蘭大學工程學院測繪與地理信息工程系， 伊朗 德黑蘭 1417466191；2.伊斯蘭阿扎德大學馬什哈德分校土木工程系， 伊朗 馬什哈德 9187147587；3.比爾詹德大學土木工程系， 伊朗 比爾詹德 9717434765)

1 研究背景

現(xiàn)代社會對自然資源的過度使用和污染物的大量排放已經威脅到地下水資源并造成嚴重污染. 自然和人為污染過程已導致全球的水質迅速下降，在發(fā)展中國家尤為嚴重[1]. 當前，地表水和地下水污染有關問題的水質調查廣泛展開[2]. 由于含水層中污染物容量高、持續(xù)時間長和自然凈化弱，一旦含水層受到污染，污染將持續(xù)存在[3-4]. 此外，改善地下水質量的含水層凈化是一個復雜的過程，需要時間和資金. 因此，控制或減少地下水污染是地下水資源管理的重要原則之一[5]. 1960年法國學者首次提出脆弱性的概念. 脆弱性是一種相對屬性，無法察覺和測量，取決于地質和水文地質環(huán)境的含水層特性[6]. 為了評估脆弱性，已經提出多種方法，包括基于過程的評價法、統(tǒng)計方法和重疊評價法[7]. 其中，基于過程的評價法利用仿真模型來估計污染物的運移；統(tǒng)計方法分析地下水空間變量與污染物量之間的相關性；重疊評價法綜合考慮污染物從地表運移到地下飽水區(qū)的控制參數(shù)，并計算一個區(qū)域不同部分的脆弱性指數(shù)[7]. 最著名的重疊指數(shù)方法之一是Aller等[8]在1987年為美國環(huán)境保護署開發(fā)的DRASTIC模型，該模型將通過不同參數(shù)獲得的信息利用地理信息系統(tǒng)(geographic information system,GIS)進行組合分析. DRASTIC模型基于水文地質概念，將特定區(qū)域的物理和水文地質特征相結合[9]. 本文DRASTIC- LU模型的主要優(yōu)點之一是使用多個信息層進行脆弱性評估，減少了錯誤信息或難以預料問題對最終結果的影響[10]. 此外，本文的模型比基本的DRASTIC模型具有更多的土地利用層[10]. DRASTIC- LU模型包括含水層埋深、凈補給量、含水層介質、土壤介質、地形、滲流帶的影響、水力傳導率和土地利用等參數(shù). 根據(jù)前人的研究，地下水脆弱性圖繪制分析可以分為知識驅動和數(shù)據(jù)驅動兩大類方法[1,11]. 數(shù)據(jù)驅動的方法對于研究充分區(qū)域或統(tǒng)計數(shù)據(jù)數(shù)量(參考數(shù)據(jù))足夠的區(qū)域非常有效，該方法目的是確定更詳細的工作位置. 知識驅動的方法對于初步研究地區(qū)、目標或數(shù)據(jù)很少的區(qū)域有效. 權重估計和分類估計是基于專家的判斷，不需要已知的數(shù)據(jù). 基于知識驅動和數(shù)據(jù)驅動集成方法在地下水脆弱性圖中已有大量研究，僅舉幾例：

Sener等[1]建立了基于Egirdir湖流域的GIS的DRASTIC模型. 在本案例研究中，使用DRASTIC和DRASTIC- AHP模型做出地下水脆弱性圖，并用于線性回歸來確定地下水脆弱性(硝酸鹽質量濃度)與含水層脆弱性區(qū)域之間的關系. 研究結果表明，與DRASTIC模型相比，DRASTIC- AHP模型與該區(qū)域的硝酸鹽質量濃度具有更大的相關性. Pisciotta等[2]使用SINTACS模型和DRASTIC模型評估地下水中硝酸鹽污染的內在脆弱性，并比較意大利地中海地區(qū)西西里島農業(yè)過度發(fā)展帶來的環(huán)境影響. 通過皮爾遜相關方法將獲得的硝酸鹽質量濃度值與SINTACS和DRASTIC值進行比較，結果表明SINTACS模型與硝酸鹽質量濃度值具有更大的相關性. Khan等[6]使用DRASTIC模型、改進的DRASTIC- LU模型、DRASTIC- AHP模型和改進的DRASTIC- LU- AHP模型來評估Raipur市的地下水脆弱性. 他們使用這4個模型將地下水脆弱性分為非常低、低、中等、高和非常高5個類別. 為了確定DRASTIC模型的準確性，總共使用了50個季風前和季風后季節(jié)的硝酸鹽質量濃度的地下水樣，并發(fā)現(xiàn)改進的DRASTIC- LU- AHP模型最準確并且最適合該研究區(qū)域. Chamanehpour等[12]使用DRASTIC模型、SINTACTS模型和FUZZY- DRASTIC模型評估Nehbandan市地下水固有脆弱性. 研究結果表明，評估該地區(qū)地下水脆弱性的最佳方法是FUZZY- DRASTIC模型，該模型與該地區(qū)硝酸鹽質量濃度的相關系數(shù)高于DRASTIC和SINTACS模型. Koh等[3]使用普通最小二乘(ordinary least squares,OLS)回歸和地理加權回歸(geographically weighted regression,GWR)方法來確定整個島嶼硝酸鹽質量濃度與各種參數(shù)(地形、水文和土地利用)之間的關系. 比較OLS模型和GWR模型發(fā)現(xiàn)， GWR模型比OLS模型能更好地反映地下水脆弱性. 這項研究的結果證明，GWR模型是研究地下水質量與環(huán)境因素之間不同空間關系的有效工具.

迄今為止，地下水脆弱性圖的制作受到廣泛關注，但是在對其進行的研究中，有一些問題沒有被充分關注. 首先，這些研究都沒有為地下水脆弱性分區(qū)提供足夠的參數(shù)組合. 其次，尚未使用適當?shù)姆治鰜泶_定有效DRASTIC- LU參數(shù)的最佳組合，并基于有效DRASTIC- LU參數(shù)做出地下水脆弱性分區(qū)圖. 本研究中，通過改進的DRASTIC模型(將土地利用圖添加到模型中)，除了確定Birjand平原含水層對污染的脆弱性外，還考慮地下水資源管理中影響含水層脆弱性的最有效參數(shù). 因此，為了達到本研究的目的，根據(jù)該地區(qū)硝酸鹽質量濃度(參考數(shù)據(jù))的適合性，結合空間和非空間數(shù)據(jù)驅動方法，包括人工神經網絡(artificial neural network,ANN)和GWR結合二元粒子群算法(binary particle swarm optimization,BPSO)，在確定有效DRASTIC- LU參數(shù)的最佳組合的基礎上，制作地下水脆弱性圖.

2 方法

2.1 研究區(qū)域

研究區(qū)域為Birjand，位于Bagheran高地的北部，經度為59°13′，緯度為32°53′. 該平原近似長方形狀，北至Moli高地和Markooh城堡，南至 Bagheran高地和Rach山，東至Bandar山和Amin Abad村莊，西至Korang和Chang高地，并在其中心區(qū)域形成沖積含水層. Birjand平原含水層補給區(qū)的總面積約為3 425 km2，其中980 km2是平原，其余為高地. 年平均降水量為170 mm，Birjand平原的最高和最低海拔分別為1 491 m和1 295 m. 該地區(qū)的土地包括旱地、荒地、農田和居民區(qū). Birjand含水層的形狀幾乎為矩形，平均長度為55 km，寬度為5 km[13]. 圖1為本研究中數(shù)學方法的流程圖.

圖1 本研究的數(shù)學方法流程

2.2 DRASTIC- LU模型

DRASTIC模型用于利用水文地質參數(shù)評估地下水相對污染脆弱性的能力. 該模型由7個可測量的水文地質參數(shù)組合而成，這些參數(shù)影響污染物運移到地下水中，這些參數(shù)包括含水層埋深、凈補給量、滲流層影響、土壤介質、含水層介質、水力傳導率和地形. 在本研究中，增加了土地利用層以提高模型的準確性.

地表和地下水位之間的距離稱為含水層埋深，表示污染物到達地下水必須經過的深度，m[14]. 凈補給量指從地表滲入地下水位的水量，mm/a[15]，此外，通過增加某個地區(qū)的凈補給量，該地區(qū)的地下水污染潛力將會增加. 滲流層指地下水位和非飽和土壤環(huán)境之間的區(qū)域[16]，該區(qū)域是不飽和的或不連續(xù)飽和的，并控制污染物的通過和稀釋. 土壤介質是非飽和區(qū)之上的區(qū)域，并且接近植物根部滲透或存在微生物活性的區(qū)域[15]. 含水層介質與飽和帶的組分特征有關，例如孔隙率、顆粒類型和顆粒大小[17]. 水力傳導率反映含水層輸水的水文地質特性和能力[16]，直接影響污染物傳輸. 該參數(shù)控制從飽和帶注入點開始的污染物的傳輸和分布[15]. 水力傳導性越大，污染物在含水層中流動的可能性就越大，因此脆弱性就越大. 地形參數(shù)是地表的坡度. 地表坡度影響污染物的運動和滲透[17-18]. 土地坡度越小，污染物與地表的接觸就越多，隨后，污染物滲入土地的可能性就越大. 因此，通過降低地表坡度，可以增加含水層脆弱性[16]. 人口增長、農業(yè)生產、工業(yè)生產和城市活動導致土地利用變化和地下水資源使用增加，不僅使這些寶貴資源的量減少了，水質也變差了[19]. 在本研究中，作者將土地利用作為DRASTIC模型的擴展參數(shù). 實際上，由于研究案例中農業(yè)和工業(yè)生產的增加，土地使用對水資源污染也具有重要影響. 這項參數(shù)已合并到DRASTIC模型中.

2.3 人工神經網絡

人工神經網絡是受人腦神經系統(tǒng)啟發(fā)的計算方法之一. 這種網絡的顯著特征之一是具有學習能力和自學習能力，并且由于該特征，使學習理解模式成為可能[20]. 相比用于預測和模擬模式的回歸方法，人工神經網絡最重要的優(yōu)勢在于，在鏈接輸入和輸出數(shù)據(jù)時不需要初始模型[21]. 根據(jù)數(shù)據(jù)之間的固有關系，在自變量和因變量之間建立線性或非線性模型.

在本研究中，多層感知神經網絡用于模擬地下水脆弱性. 這種類型的神經網絡由一組連續(xù)排列在不同層中的神經元組成. 多層感知學習定律稱為逼近訓練規(guī)則，用于估計未知網絡參數(shù). 多層感知器的工作方式是將模式提供給網絡并計算輸出. 實際輸出和所需輸出會導致網絡系數(shù)發(fā)生變化，這樣可以在以后的階段獲得更準確的輸出. 為了成功進行網絡訓練，必須逐漸使其輸出接近所需的輸出并且降低誤差. 本研究中使用的人工神經網絡設計如圖2所示.

圖2 本研究中人工神經網絡架構應用[21]

2.4 地理加權回歸

根據(jù)空間數(shù)據(jù)的空間自相關性和空間非平穩(wěn)性，使用普通最小二乘法[22]等基礎全局回歸可能性較小. 在這種方法中，事件之間的空間相關性被視為權重矩陣，由于環(huán)境因素的差異性和局部變化的存在，用于觀測的GWR模型的回歸系數(shù)需要局部測量[23-24]. GWR模型的計算公式[24]為

(1)

式中：yi是因變量(硝酸鹽質量濃度)；xj是自變量(DRASTIC- LU參數(shù))；n是總隨機點數(shù)；εi是GWR模型殘差；(ui，vi)表示第i個點在坐標中的坐標；βj(ui，vi)是第i個點坐標的第j個回歸系數(shù).為了計算空間權重矩陣，有必要指定所需的核函數(shù).根據(jù)先前的研究，本研究使用指數(shù)和雙平方這2個核函數(shù)[24-25]，分別是

(2)

(3)

(4)

式中：dij是2個觀測值i和j之間的歐幾里得距離值，并且是2個觀測值i和j沿水平和垂直軸之間的距離；b是帶寬值.每個位置的回歸系數(shù)都不同，因此在GWR模型中，可以根據(jù)等式通過標準偏差函數(shù)公式[24]獲得回歸系數(shù)的局部變化

(5)

式中：βij是觀測值i中因子j的回歸系數(shù)；βj是總觀測值中j因子的平均回歸系數(shù)；n是總觀測值.

為了評估ANN和GWR模型的輸出，通常使用確定系數(shù)(R2)來衡量擬合度，ERMS值則用于衡量觀測值的殘差分布[24]

(6)

(7)

式中：n是總觀測值；yi是觀測值i的值；i是觀測值i的預測值；是總觀測值的平均值.

2.5 二進制粒子群優(yōu)化算法

微粒群優(yōu)化算法(particle swarm optimization,PSO)算法是一種優(yōu)化算法，它不太可能在局部最小值處被捕獲，并且可以根據(jù)概率規(guī)則搜索不確定和復雜的區(qū)域[26]. 同樣在該算法中，提出的路徑的解不依賴于初始種群，并且從搜索空間的每個點開始，收斂到最優(yōu)解[26]. 隨后Kennedy等[27]引入二進制PSO算法，與PSO的連續(xù)版本不同，該算法僅限于0和1(二進制)變量，并且速度值可以將粒子從0更改為1. 根據(jù)本研究的目的，選用BPSO算法.

在此算法中，用于更新每個粒子的速度和位置的公式[27]為

Vi(t+1)=wVi(t)+c1r1(pb-Xi(t))+
c2r2(gb-Xi(t))

(8)

Xi(t+1)=Xi(t)+Vi(t+1)

(9)

(10)

式中：Vi(t)是質點i的速度；Xi(t)是質點i的位置；Vi(t+1)是質點i在下一位置的速度；Xi(t+1)為粒子i在下一個位置中的位置；pb是粒子i歷程的最佳位置；gb是所有粒子中歷程的最佳位置；c1是個體學習系數(shù)；c2是集體學習系數(shù)；w是慣性權重；r1、r2和ρ是[0，1]范圍內的隨機數(shù).在這項研究中，BPSO算法的步驟(結合ANN和GWR模型)如圖3所示.

圖3 本文模型的計算步驟

步驟1隨機賦予粒子群的位置和速度的初始值.

步驟2訓練ANN和GWR模型，并計算總體中每個粒子的適應度函數(shù)(R2).

步驟3如果達到停止標準(100次迭代)，則停止BPSO算法，否則轉到步驟4. 如果算法達到停止條件，則選定的DRASTIC- LU參數(shù)與利用地下水脆弱性圖的有效參數(shù)編制的有效參數(shù)相同.

步驟4確定粒子的pb和gb.

步驟5計算每個粒子的速度，并根據(jù)相互關系移動到下一個位置(轉到步驟2).

3 結果與討論

3.1 定義因變量和自變量

根據(jù)表1，在本研究中，DRASTIC- LU參數(shù)為自變量，這些參數(shù)按表1中順序構成BPSO算法的粒子維數(shù).

表1 本研究中的自變量

圖4說明了該案例研究的DRASTIC- LU參數(shù). 通過安裝在含水層中的15個壓力計獲得2016年10月至2017年9月的水深數(shù)據(jù). 凈補給量借助觀測井中的水位數(shù)據(jù)確定含水層地下水量的分區(qū)方法. 含水層介質和滲流層影響數(shù)據(jù)通過South Khorasan區(qū)域供水組織在該地區(qū)的19口觀測井和生產井記錄確定. 區(qū)域地形使用具有15 m空間分辨率的數(shù)字高程模型圖. 水力傳導率是根據(jù)South Khorasan區(qū)域供水組織進行的抽水試驗數(shù)據(jù)計算得出. 土壤類型使用South Khorasan農業(yè)組織的遙感和地圖數(shù)據(jù)得出. 土地利用是基于2017年Google Earth Engine中的Landsat 8衛(wèi)星圖像. DRASTIC- LU地圖繪制采用克里金插值法[13]，因為該方法具有最小的誤差.

圖4 DRASTIC- LU參數(shù)

在該試驗中，采用Birjand自來水和廢水處理部門的水質實驗室提供的Birjand平原含水層中的21口井的硝酸鹽質量濃度的數(shù)據(jù). 應用的數(shù)據(jù)是2014年至2019年5 a間硝酸鹽質量濃度的平均值. 圖5(a)顯示了采用克里金插值法生成的整個Birjand含水層中硝酸鹽質量濃度.

圖5 基于觀察井和隨機點的因變量

根據(jù)圖5(b)，將DRASTIC- LU參數(shù)作為自變量、硝酸鹽質量濃度作為因變量計算后，通過隨機采樣方法在研究區(qū)域中創(chuàng)建了1 000個點并用于GWR- BPSO模型[28]. 然后計算這些點的所有可用信息參數(shù)(自變量和因變量)的值(以標準化的方式). 利用公式

(11)

1～8含義參見表1. 圖6 DRASTIC- LU參數(shù)之間的相關矩陣

3.2 實施數(shù)據(jù)驅動的模型

為了實現(xiàn)ANN和GWR模型，總數(shù)據(jù)的70%用于訓練，總數(shù)據(jù)的30%用于測試，并且在進入算法之前將所有數(shù)據(jù)進行標準化. 根據(jù)先前的研究和反復試驗的方法，選擇70∶30的比例[30]. 在本文中，該比率可提供最佳性能結果. 由于評估GWR模型的最重要參數(shù)之一(模型與數(shù)據(jù)的兼容性)是確定系數(shù)參數(shù)(R2)，因此選擇BPSO算法適應度函數(shù)以最小化1-R2的值[24]. BPSO算法的初始參數(shù)的最佳值是根據(jù)從不同迭代獲得的試驗以及根據(jù)表2的反復試驗來選擇的. 為簡化實現(xiàn)過程而停止的條件是具體執(zhí)行的次數(shù).

表2 BPSO算法中的給定參數(shù)

圖7顯示本研究中BPSO算法的群結構，表1中提到的DRASTIC- LU參數(shù)構成BPSO算法的粒子維.

圖7 BPSO算法的Swarm結構

由于BPSO算法的隨機性，在前人研究的基礎上，將該算法在期望的迭代次數(shù)下重復執(zhí)行10次，并將這10次執(zhí)行中最好的1次作為最終輸出[31]. 根據(jù)圖8(a)，通過將ANN模型與BPSO算法相結合，得到適應度函數(shù)(1-R2)的最佳值為0.106 0(10次執(zhí)行中的最佳值). 根據(jù)圖8(b)的ANN模型，硝酸鹽質量濃度的有效DRASTIC- LU參數(shù)是通過含水層介質、土地利用、土壤介質、地形和水力傳導率5個參數(shù)確定. 實際上，圖8(b)為在BPSO算法的第100次迭代中，所有粒子(30個粒子)的適應度函數(shù)(R2)最優(yōu)粒子.

圖8 基于ANN模型和BPSO算法計算的相關適應度函數(shù)值

然后，將具有指數(shù)核和雙平方核的GWR模型與BPSO算法相結合，以確定預測硝酸鹽質量濃度的有效DRASTIC- LU參數(shù). 為了實現(xiàn)GWR模型，將隨機點的坐標用作權重矩陣的輸入. 根據(jù)圖9，得到將具有指數(shù)核和雙平方核的GWR模型和BPSO算法組合的適應度函數(shù)(1-R2)的最佳值分別為0.074 5和0.006 5(10次執(zhí)行中的最佳值). 同樣，根據(jù)圖10，對于具有指數(shù)核的GWR模型，確定了含水層介質、含水層埋深、滲流層影響、土地利用、土壤介質和水力傳導率等6個有效參數(shù)；對于雙平方核，確定了含水層介質、含水層埋深、滲流層影響、土地利用、壤介質，地形和水力傳導率等7個有效參數(shù). 這些含水層參數(shù)確定為預測硝酸鹽質量濃度的有效DRASTIC- LU參數(shù).

圖9 通過結合GWR模型和BPSO算法獲得適應度函數(shù)的最佳值

圖10 通過結合GWR模型和BPSO算法預測硝酸鹽質量濃度的有效DRASTIC- LU參數(shù)

圖11顯示了ANN和GWR模型的R2和RMSE值. 因此， DRASTIC- LU參數(shù)基礎上雙平方核在預測硝酸鹽質量濃度方面具有更高的準確性.

圖11 依據(jù)R2和RMSE進行ANN和GWR模型比較

根據(jù)圖12，結合ANN和GWR模型以及BPSO算法預測的硝酸鹽質量濃度(基于DRASTIC- LU參數(shù))的圖譜顯示在[0,1]范圍內. 根據(jù)等間隔分類方法，將估計的硝酸鹽質量濃度分為5個輸出類別，從極低的硝酸鹽質量濃度到非常高的硝酸鹽質量濃度定性顯示. 根據(jù)從R2值(擬合好的)和RMSE值(觀測值的殘差分布以及模型的準確性)獲得的結果，將GWR模型與雙平方核和BPSO算法結合使用具有更高的估算硝酸鹽質量濃度的能力(地下水脆弱性)，見圖12(c).

圖12 估算的硝酸鹽質量濃度

如上所述，由于GWR模型中每個位置的回歸系數(shù)都不同，因此可以通過標準偏差函數(shù)獲得回歸系數(shù)的局部變化和空間非平穩(wěn)性. 作者采用前述雙核的GWR方法，研究了DRASTIC- LU參數(shù)與地下水脆弱性隨位移變化的最大值和最小值. 圖13顯示了用于計算局部變化率和空間非平穩(wěn)性的回歸系數(shù)GWR模型(具有指數(shù)核和雙平方核)的標準偏差.

圖13 具有指數(shù)核和雙平方核的回歸系數(shù)GWR模型的標準偏差

根據(jù)圖13，對于具有指數(shù)核的GWR模型，滲流層影響參數(shù)影響與地下水脆弱性(硝酸鹽質量濃度)之間的關系最大，而土地利用參數(shù)與地下水脆弱性(硝酸鹽質量濃度)之間的關系則最小. 同樣，在具有雙平方核的GWR模型中，滲流層影響參數(shù)與地下水脆弱性(硝酸鹽質量濃度)的影響之間的關系最大，而土壤介質參數(shù)與地下水脆弱性(硝酸鹽質量濃度)之間的關系最小. 最后，使用全局Moran’s指數(shù)確定GWR模型殘差的空間自相關，該空間自相關公式為[32]

(12)

表3 具有2個指數(shù)和雙平方核的GWR模型殘差的Moran指數(shù)值

4 結論

1) 采用基于DRASTIC- LU參數(shù)的空間和非空間數(shù)據(jù)驅動的方法，包括GWR和ANN模型，來預測硝酸鹽質量濃度. 結果表明，基于空間自相關和空間非平穩(wěn)性的特點，所使用的基于DRASTIC- LU參數(shù)的GWR模型預測硝酸鹽質量濃度方面具有更高的準確性.

2) 將二進制粒子群優(yōu)化算法與ANN和GWR模型結合使用，表明DRASTIC- LU參數(shù)對研究區(qū)域的硝酸鹽質量濃度預測具有顯著影響. 重要的是該方法不限于本研究，而且可以用于預測不同類型區(qū)域中的硝酸鹽質量濃度.

3) 對于未來的工作，建議延長統(tǒng)計周期，如果可能的話，使用更多的取樣井可以幫助更好地建立和評估模型. 此外，建議采用遺傳算法、蟻群算法等其他進化方法. 由于對帶有土地利用參數(shù)的改進后的DRASTIC模型在研究半干旱地區(qū)取得了較好的結果，因此建議結合其他指標對該模型進一步進行改進.

1 Introduction

Nowadays, overuse of natural resources and abundant production of waste in modernized societies have threatened groundwater resources and caused extreme contamination. Natural and anthropogenic pollution processes have caused a rapid decline in water quality globally and mostly in developing countries[1]. Today, water quality surveys have become widespread and include issues related to surface water and groundwater pollution[2]. Once an aquifer is contaminated, the contamination becomes persistent as a result of high capacity, longer persistence time, and lack of physical accessibility[3-4]. On the other hand, decontamination which struggles to improve the groundwater quality is a complicated process that needs time and financial resources. Hence, groundwater contamination control or reduction is one of the important principles in groundwater resource management[5]. The concept of vulnerability was first introduced in France in 1960. Vulnerability is a type of relative attribute, without smell and non-measurable, and depends on the aquifer properties of the geological and hydrogeological environment[6]. Various methods have been presented in order to evaluate the vulnerability including the processing methods, statistical methods, and overlapping methods[7]. The processing methods use simulation models to estimate contaminant movement. The statistical methods use the correlation between spatial variables and the quantity of the contaminants in groundwater. Also, overlapping methods integrate the controller parameters of contaminant movement from the ground surface into the saturation zone and determine an index called the vulnerability index in different parts of a territory[7]. One of the most well-known overlap index methods is the DRASTIC model developed by Aller in 1987 for the United States Environmental Protection Agency[8]. This model is an overlapping method in which information obtained from different parameters is analyzed in a combined way and then processed by GIS. This model is also based on the concept of hydrogeological status and can combine the physical and hydrogeological characteristics of a specific area[9]. One of the main advantages of the DRASTIC-LU model is the vulnerability evaluation using several numbers of informational layers, which limits the effect of errors or unpredictable problems on the final output[10]. Moreover, the suggested model has a layer of land use more than the basic DRASTIC one. The DRASTIC-LU model includes the parameters of depth to water, net recharge, aquifer media, soil media, topography, impact of vadose zone, hydraulic conductivity, and land use. According to previous research, groundwater vulnerability map production can be classified into two general categories of knowledge-driven and data-driven methods[1,11]. Data-driven methods are highly effective in known areas or areas where the number of known evidence is statistically sufficient (references data). In these methods, the goal is to identify new locations for more detailed work, while in knowledge-driven methods, they are effective in environments that are less known or where there are few targets in the area. Weight estimates and class estimates are based on expert judgment and do not require evidence of an answer. Numerous studies have been conducted on groundwater vulnerability maps with approaches based on knowledge-driven and data-driven integration methods, to name a few:

Sener et al.[1]used the DRASTIC model based on a geographic information system (GIS) in Egirdir Lake basin. Therefore, DRASTIC and DRASTIC-AHP models were used to prepare groundwater vulnerability maps in the case study. Finally, they are used as linear regression to examine the relationship between groundwater vulnerability (nitrate concentration) and the aquifer vulnerability areas. The results showed that the DRASTIC-AHP model has more correlation with the nitrate concentration in the region compared to the DRASTIC model. Pisciotta et al.[2]used two models of the SINTACS and DRASTIC to evaluate the intrinsic vulnerability against nitrate contamination vulnerability in groundwater and to compare the environmental effect of excessive agricultural development in Sicily, the Mediterranean region of Italy. The obtained nitrate concentration values were compared with the values of the SINTACS and DRASTIC by Pearson Correlation Method and the results demonstrated that the SINTACS model has a further correlation with nitrate values. Khan et al.[6]used the DRASTIC, modified DRASTIC-LU, DRASTIC-AHP, and modified DRASTIC-LU-AHP models to assess the groundwater vulnerable map in Raipur city. Using these four models, they have divided the groundwater vulnerability maps into 5 classes, namely, very low, low vulnerability, moderate, high and very high. To determine the accuracy of DRASTIC models, a total of 50 groundwater samples of nitrate concentration for pre-monsoon and post-monsoon seasons were used and it was observed that the modified DRASTIC-LU-AHP model is the most accurate and appropriate study area. Chamanehpour et al.[12]in their study evaluated the imperative and undeniable vulnerability of groundwater in Nehbandan city using three models DRASTIC, SINTACTS and FUZZY-DRASTIC. The results showed that the best method for estimating the groundwater vulnerability of the region is the FUZZY-DRASTIC model, which has a higher correlation coefficient with the nitrate concentration of the region than the DRASTIC and SINTACS models. Koh et al.[3]used ordinary least squares (OLS) regression and geographically weighted regression (GWR) methods to determine the relationships between nitrate concentrations and various parameters (topography, hydrology, and land use) across the island. Comparison between OLS and GWR models showed that the GWR model has better performance in preparing the groundwater vulnerability map than the OLS model. The results of this study verified that the GWR model is a useful tool to investigate the different spatial relationships between groundwater quality and environmental factors.

Groundwater vulnerability map production is a topic that has received a lot of attention so far; but among the studies conducted, there are points that have received less attention. First, none of these studies provide an adequate combination of parameters for groundwater vulnerability zoning. Second, proper analysis has not been used to determine the optimal combination of effective DRASTIC-LU parameters and to prepare a groundwater vulnerability zoning map based on the effective DRASTIC-LU parameters. In this study, by modifying the DRASTIC model (adding land use maps to model maps), in addition to recognizing the vulnerability of Birjand plain aquifer to pollution, the most effective parameters affecting the aquifer vulnerability to manage groundwater resources are also considered. Therefore, to achieve the purpose of this study, due to the availability of nitrate concentration in the region (reference data), the combination of spatial and non-spatial data-driven methods including artificial neural network (ANN) and geographically weighted regression (GWR) with binary particle swarm optimization (BPSO) were used to prepare a groundwater vulnerability map based on determining the optimal combination of effective DRASTIC-LU parameters.

2 Methods

2.1 Study Area

Birjand is the region where the experiment was conducted, and it is located in the northernpart of the Bagheran highlands at the longitude of 59°13′ and the latitude of 32°53′. This plain is rectangular and is limited to Moli highlands and Markooh in the north, Bagheran highlands and Rach Mountain in the south, Bandar Mountain and Amin Abad in the east, and Korang and Chang highlands in the west, as well as the alluvial aquifer which has been formed in its central zone. The total area of Birjand plain aquifer basin is about 3 425 km2of which 980 km2is plain and the rest is occupied by highlands. The average precipitation has been 170 millimeter per year, and the maximum and minimum of Birjand plain altitude is 1 491 m and 1 295 m respectively. The region land use includes the dryland, abandoned lands, farmlands, and residential areas. In total, Birjand aquifer is almost rectangular in shape, which its average length is 55 km and its width is 5 km[13]. Also, Fig.1 shows the flowchart of the proposed implementation method in this study.

Fig.1 Flowchart showing the methodology used in this study

2.2 DRASTIC-LU model

DRASTIC model has been used to evaluate the capability of groundwater relative contamination vulnerability, using hydrogeological parameters. This model has been constituted from the combination of seven measurable hydrogeological parameters, which are effective on the transfer of contaminants into groundwater, and comprising the depth to water, net recharge, impact of vadose zone, soil media, aquifer media, hydraulic conductivity and topography. At the present study, the land use layer has been added to enhance the accuracy of the model as mentioned earlier.

The distance between the land surface and the water table is called depth to water. In other words, this parameter indicates the depth that a contaminant must pass through to reach groundwater, m[14]. Net recharge is the amount of water that penetrates from the ground surface and reaches into a water table, mm/a[15]. Furthermore, through the increase of net recharge value in a region, the contamination potential of groundwater will rise in that region. The impact of vadose zone comprises the zone between the water table and the unsaturated soil environment[16]. This zone is unsaturated or saturated discontinuously and controls the contaminant passing and its dilution. The soil media is the area above the non-saturated zone and approaches the zone where plant roots penetrate in, or the activity of microorganisms exists[15]. The aquifer media is associated with the component characteristics of the saturated zone, such as the porosity rate, particle types and size[17]. The capability of the hydrogeological characteristics of the aquifer in water transferring and associated contaminants is called hydraulic conductivity[16]. This parameter controls the transport and distribution of contaminants from the injection point inside the saturated zone[15]. The further the hydraulic conductivity is the more will be the possibility of contaminants flowing in the aquifer, so the vulnerability will be more. The topography parameter is the slope of the land surface. The land surface slope affects the movement and penetration of contaminants[17-18]. The less the land slope is, the more will be the contact of a contaminant with the land surface, and subsequently, the more likely will be its penetration into the land. Hence through the decrease of land surface slope, the possibility of aquifer vulnerability increases[16]. Population growth, agricultural, industrial and urban activities and thereupon increasing the land-use change and the usage of groundwater resources, not only has diminished but also has degraded the quality of these valuable resources[19]. In this study, land use as an extension parameter to DRASTIC model has been considered. Actually, due to the increased amount of agricultural and industrial activities in the research case study, it seems that land use has an essential impact on water resource pollutions. Then, this parameter has been integrated into the DRASTIC model.

2.3 Artificial neural network

Artificial neural networks are one of the computational methods inspired by the neural system of the human brain. One of the remarkable characteristics of this type of network is their ability to learn and the ability to generalize this learning, and because of this feature, they make it possible to learn to understand patterns[20]. The most important advantage of artificial neural networks over regression methods for predicting and modeling a pattern is that there is no need for an initial model in linking input and output data[21]. Based on the intrinsic relationships between data, a linear or nonlinear model is established between independent and dependent variables.

In this study, a multilayer perceptron neural network has been used to model groundwater vulnerability. This type of neural network consists of a set of neurons arranged in different layers in a row. The law of multilayer perceptron learning is called the error propagation rule, which is used to estimate unknown network parameters. The multilayer perceptron works in such a way that a pattern is supplied to the network and its output is calculated. Actual output values and desired output cause the network coefficients to change; in such a way that a more accurate output is obtained in later stages. To succeed in network training, its output must be gradually brought closer to the desired output and the error rate must be reduced. The design ANN used in this study is illustrated in Fig.2.

Fig.2 ANN architecture applied in this study[21]

2.4 Geographically Weighted Regression

According to spatial autocorrelation and spatial non-stationarity properties for spatial data, it is less possible to use basic global regressions such as ordinary least square[22]. In this method, the spatial dependencies between the events are considered as weight matrices, and due to the heterogeneity of the environmental factors and the existence of local variation, regression coefficients of the GWR model for observation are measured locally[23-24]. The GWR model is calculated as equation (1)[24]:

(1)

Whereyiis the dependent variable (nitrate concentration);xjis the independent variables (DRASTIC-LU parameters);nis the total random points;εiis the residual GWR model; (ui,vi) denotes the coordinates of the ith point in space andβj(ui,vi) is the regression coefficient for coordinates of the ith point. To calculate the spatial weight matrix, it is necessary to specify the desired kernel function. According to previous research, this study used two kernels including exponential and bi-square these two kernels which are calculated as equation (2) and equation (3), respectively[24-25]:

(2)

(3)

(4)

(5)

Whereβijis the regression coefficient for the factorjin the observationi;βjis the mean regression coefficient of factorjin the total observations andnis the total observations.

To evaluate the ANN and GWR models output the coefficient of determination (R2) is usually used to measure the goodness of fit and theERMSvalue measure the residuals distribution of the observation, which are obtained based on equation (6) and equation (7)[24]

(6)

(7)

Wherenis the total observations;yiis the value for observationi;iis the predicted value for observationiandis the mean value for total observations.

2.5 Binary particle swarm optimization

The particle swarm optimization (PSO) algorithm is an optimization algorithm that makes it less likely to be captured at a local minimum and can search uncertain and complex areas based on probabilistic rules[26]. Also in this algorithm, the solution of the proposed path is not dependent on the initial population and starting from each point in the search space, the solution converges to the optimal solution[26]. After a while, Kennedy and Eberhart[27]introduced the binary PSO algorithm, which, unlike the continuous version of it, is limited to having zero and one (binary) variables and the velocity value can change a particle from zero to one. According to the purpose of this study, the BPSO algorithm has been used. In this algorithm, equation (8)-(10) are used to update the velocity and position of each particle[27]:

Vi(t+1)=wVi(t)+c1r1(pb-Xi(t))+
c2r2(gb-Xi(t))

(8)

Xi(t+1)=Xi(t)+Vi(t+1)

(9)

(10)

WhereVi(t) is the velocity of the particlei;Xi(t) is the position of the particlei;Vi(t+1) is the velocity of the particleiin the next position;Xi(t+1) is the position of the particleiin the next position;pbis the best position of the experience for the particlei;gbis the best position experienced in all particles;c1is the personal learning coefficient;c2is the collective learning coefficient;wis the inertia weight andr1,r1andρare random numbers in the range [0.1]. In this study the steps of the BPSO algorithm (in combination with the ANN and GWR models) are as follows which showed in Fig.3:

Fig.3 Calculation steps of the recommended models

Step1Give the initial value to a population of particles with random positions and velocities.

Step2Training ANN and GWR models and calculating the fitness function (R2) of each particle in this population.

Step3Stop the BPSO algorithm if it reaches the stop criterion (100 iterations), otherwise go to step 4. If the algorithm reaches the condition of stopping, then the selected DRASTIC-LU parameters are the same effective parameters that are prepared with the help of effective parameters of the groundwater vulnerability map.

Step4Determine the pbest and gbest for particles.

Step5Calculate the velocity of each particle and move to the next position based on the relations (go to step 2).

3 Results and Discussion

3.1 Define dependent and independent variables

According to Table 1, in this study, the DRASTIC-LU parameters are considered as independent variables. These parameters, in the order mentioned, form the particle dimension of the BPSO algorithm.

Table 1 Independent variables in this study

Fig.4 illustrates the DRASTIC-LU parameters of the case study. Depth of water has been obtained through 15 piezometers installed in the aquifer in a one-year period from October 2016 to September 2017. In this study, in order to obtain net recharge, the zoning method of changes in aquifer groundwater volume with the help of water level data in the observation wells has been used. Aquifer media and impact of vadose zone were achieved through a log of 19 observation and operation wells in the region through the Regional Water Organization of Southern Khorasan Province. The topography of area has been obtained using a digital elevation model map with a spatial resolution of 15 meters. Hydraulic conductivity has been calculated from the pumping experiment data performed by the Regional Water Organization of South Khorasan Province. For soil type, remote sensing and map data prepared by Agriculture Organization of Southern Khorasan Province were used. Finally, land use was based on Landsat 8 satellite imagery in the Google Earth Engine environment captured in 2017. For the production of DRASTIC-LU maps, a kriging interpolation has been used since this method has the minimum error[13].

Fig.4 Maps of DRASTIC-LU parameters

In this experiment, the nitrate contents of 21 wells located in Birjand plain aquifer based on the data from the water quality lab of the Water and wastewater department of Birjand city have been used. The applied data is an average of nitrate concentration during five years from 2014-2019. To show the nitrate concentrations all over the Birjand aquifer (10-6, unit); a kriging interpolation has been used (Fig.5(a)).

Fig.5 Map of dependent variable based on observation wells and random points

According to Fig.5(b), after calculating the DRASTIC-LU parameters as independent variables and nitrate concentration as dependent variable, 1 000 points created in the study area by a random sampling method to extract values and use in the GWR-BPSO model[28]. Then the values of all available information parameters (independent and dependent variables) for these points were calculated (in a normalized way). The correlation between the DRASTIC-LU parameters from equation (11) was examined[29]:

(11)

1-8 shown as Table 1. Fig.6 Correlation matrix between DRASTIC-LU parameters

3.2 Implement data-driven models

For the implementation of the ANN and GWR models, 70% of the total data was used for training and 30% of the total data was used for testing, and all data were normalized before entering the algorithms. Based on previous research and trial and error method, a ratio of 70∶30 was selected[30]. In article this ratio gives the best performance result. Due to the fact that one of the most important parameters for evaluating the GWR model (model compatibility with data) is the coefficient of determination parameter (R2), therefore, the BPSO algorithm fitness function has been selected to minimize the value of 1-R2[24]. The optimal values of the initial parameters of the BPSO algorithm were selected based on the experiments obtained from different iterations and through trial and error according to Table 2. The condition for stopping to simplify the implementation process is the number of specific executions.

Table 2 Set parameters in the BPSO algorithm

Fig.7 shows the swarm structure of the BPSO algorithm in this study, which the DRASTIC-LU parameters mentioned in Table 1 form the particle dimension of the BPSO algorithm.

Fig.7 Swarm structure of the BPSO algorithm

Due to the random nature of the BPSO algorithm and based on previous research, this algorithm with the desired number of iterations was repeated 10 executions and the best of these 10 executions was considered as the final output[31]. According to Fig.8(a), by performing the combination of the ANN model and the BPSO algorithm, the best value of fitness function (1-R2) was obtained 0.106 0 (the best of 10 executions). Also, according to Fig.8(b) for the ANN model, 5 parameters of aquifer media, land use, soil media, topography and hydraulic conductivity were determined as effective DRASTIC-LU parameters in predicting nitrate concentration. In fact, Fig.8(b) shows the best particle in terms of fitness function (R2) among all particles (30 particles) in the 100th iteration of the BPSO algorithm.

Fig.8 Results of ANN+BPSO model

Then, the GWR model with two exponential and bi-square kernels and BPSO algorithm was combined to determine the effective DRASTIC-LU parameters in predicting nitrate concentration. To implement the GWR model, the coordinates of the random points were used as inputs in its weight matrix. According to Fig.9, the best value of fitness function (1-R2) for combination the GWR model with two exponential and bi-square kernels and the BPSO algorithm, was obtained 0.074 5 and 0.006 5 (the best of 10 executions), respectively. Also, according to Fig.10 for the GWR model with the exponential kernel, 6 parameters of aquifer media, depth of water, impact of vadose zone, land use, soil media and hydraulic conductivity and for the bi-square kernel, 7 parameters of aquifer media, depth of water, impact of vadose zone, land use, soil media, topography and hydraulic conductivity were determined as an effective DRASTIC-LU parameters in predicting nitrate concentration.

Fig.9 The best value of fitness function by combining of GWR model and BPSO algorithm

Fig.10 Effective DRASTIC-LU parameters in predicting nitrate concentration by combining of GWR model and BPSO algorithm

In Fig.11, the values ofR2and RMSE for the ANN and GWR models are shown. Accordingly, the bi-square kernel has higher accuracy in predicting nitrate concentration based on the DRASTIC-LU parameters.

Fig.11 Comparison of ANN and GWR models in terms of R2 and RMSE

According to Fig.12, the maps of predicted nitrate concentration (based on the DRASTIC-LU parameters) by combining ANN and GWR models and BPSO algorithm showed in the range [0,1]. The estimated nitrate concentration is classified into five output classes according to the equal interval classification method which is shown qualitatively from very low nitrate concentration to very high nitrate concentration. According to the results obtained fromR2value (goodness of fit) and the RMSE value (residuals distribution of the observation and accuracy of model), the combination of GWR model with bi-square kernel and BPSO algorithm has a higher ability to estimate nitrate concentration (groundwater vulnerability), which showed in Fig.12(c).

Fig.12 Map of estimated nitrate concentration

As mentioned, since the regression coefficients are different for each location in the GWR model,local variation and spatial non-stationarity of the regression coefficients can be obtained by the standard deviation function. Therefore, for the GWR method with the two kernels mentioned, the most and least variation between the DRASTIC-LU parameters and groundwater vulnerability with displacement have been investigated. Fig.13 shows the standard deviation of regression coefficients GWR model (with two exponential and bi-square kernels) for calculating the rate of local variation and spatial non-stationarity.

Fig.13 Standard deviation of regression coefficients GWR model with exponential and bi-square kernels

According to Fig.13, for the GWR model with the exponential kernel, the relationship between impact of vadose zone and groundwater vulnerability (nitrate concentration) with displacement has the most variation and the relationship between land use parameter and groundwater vulnerability (nitrate concentration) has the least variation. Also, in the GWR model with the bi-square kernel, the relationship between impact of vadose zone parameter and groundwater vulnerability (nitrate concentration) with displacement has the most variation and the relationship between soil mediaparameter and groundwater vulnerability (nitrate concentration) has the least variation. Finally, global Moran’s index was used to determine the spatial autocorrelation of GWR model residuals, which is calculated from equation (12)[32]:

(12)

Table 3 Values of Moran’s index for GWR model residuals with two exponential and bi-square kernels

4 Conclusions

1) To achieve the main purpose of this study, the spatial and non-spatial data-driven methods including GWR and ANN model were used to predict the nitrate concentration based on the DRASTIC-LU parameters. The results showed that the GWR model used, taking into account the characteristics of spatial autocorrelation and spatial non-stationarity, has higher accuracy in predicting nitrate concentration based on the DRASTIC-LU parameters.

2) In this study, the binary particle swarm optimization algorithm was used in combination with the ANN and GWR models, which showed that the DRASTIC-LU parameters have a significant effect on predicting nitrate concentration in the study area. The important point is that the mentioned method is not limited to this case study and can be used to predict the nitrate concentration in various types of regions.

3) For future works, we propose a more extended statistical period and, if possible, more sampling wells can help building and evaluating the model better. Also we suggest using other evolutionary methods such as genetic algorithm, ant bee colony algorithm, etc. Since the improvement of the DRASTIC model with land use parameter has led to better results in the semi-arid region studied, further improvements of the model in combination with other indices are suggested.