基于人工神經網絡的晉祠泉水位模擬研究

邢立文，崔寧博

（1.山西省水利水電科學研究院，山西 太原030000；2.四川大學水利水電學院，四川成都610065）

我國北方地區(qū)分布有大片碳酸鹽巖巖溶地層，這些地區(qū)的巖溶地下水資源非常豐富。與南方巖溶地區(qū)的巖溶地下水不同，北方的巖溶大泉具有匯水面積廣，徑、補、蓄、排獨立，泉水排泄集中，流量穩(wěn)定且集中，水質優(yōu)良，地下調節(jié)庫容巨大的特點，多數具備建立大型供水水源地的有利條件，在供水方面有很大優(yōu)越性，是北方很多城市及工農業(yè)供水水源地。山西、陜西、山東、河南4省巖溶地下水資源開采量已達24.39億m3/a[1]。改革開放以來，我國北方地區(qū)城鎮(zhèn)化進程發(fā)展迅速、工農業(yè)及第三產業(yè)高速發(fā)展，巖溶地下水資源過量開發(fā)和煤炭產業(yè)粗放式擴張，對巖溶地下水循環(huán)條件造成了嚴重破壞和影響，很多巖溶大泉出現斷流、水位下降、污染等問題，甚至引發(fā)環(huán)境地質災害。我國北方巖溶地下水資源迫切需要從單純的開發(fā)利用向開發(fā)利用與保護管理并重的可持續(xù)發(fā)展轉變。

人工神經網絡（ANN）是由大量神經元按一定方式連接而成用來模仿人腦神經網絡結構和功能建立的靈活而復雜的數學網絡。神經網絡算法具備線性和非線性學習算法的能力，根據已經收集的信息來傳輸，收集和學習信息。人工神經網絡根據已知的輸入和輸出數據進行訓練，可以預測未來的輸入對應的輸出。

1 研究區(qū)概況

山西是我國北方巖溶泉水最多的省份，巖溶泉水占地下水資源總量的33.52%[2]，為山西省能源基地的建設提供了有力保障。據統計，流量大于0.10 m3/s的巖溶泉水有83處，其中流量大于1.00 m3/s的有這18處，這18處泉水多年平均總流量為98.84 m3/s。巖溶地下水在山西省生產、生活和文化建設中具有重要支撐作用，但20世紀80年代后，受自然環(huán)境變化、人類活動等因素影響，巖溶地下水日趨嚴重的環(huán)境問題已成為山西省國民經濟可持續(xù)發(fā)展的制約因素。晉祠泉屬汾河流域，是太原市排名第二的巖溶大泉，20世紀50年代平均流量為1.99 m3/s，年際泉水流量變化幅度為0.20 m3/s左右，屬穩(wěn)定型泉水。晉祠泉具有保障太原市經濟建設供水的重任[3]，但是大量開采地下巖溶水資源導致晉祠泉水位、流量大幅度衰減，年均泉水流量由20世紀60年代的1.69 m3/s下降至70年代的1.13 m3/s，再由80年代的0.46 m3/s下降至90年代初的0.15 m3/s，最終于1994年4月30日斷流，嚴重影響當地工農業(yè)生產，同時使旅游景觀大為遜色，并危及千年文物古跡，引起了有關部門的高度重視。

2 研究方法

2.1 泉水位影響因素及其計算方法

晉祠巖溶水系統補給來源為巖溶裸露區(qū)降水入滲補給和汾河在巖溶河段的河道滲漏補給。晉祠巖溶水系統排泄主要分為人類活動排泄和自然排泄兩類，其中自然因素包括晉祠泉天然排泄、以潛流形式向平泉和太原盆地孔隙水排泄，人類活動排泄包括晉祠泉域內巖溶水開采和礦坑排水。

（1）降水入滲補給。降水入滲是晉祠泉域巖溶地下水系統的主要補給來源，晉祠泉域面積為2 030 km2，其中裸露可溶巖面積391 km2，泉域內其他地區(qū)入滲條件較差，所以晉祠泉域入滲補給只考慮裸露區(qū)降水入滲補給量。

式中：Q為降水入滲補給量，萬m3/a；α為降水入滲補給系數；P為裸露區(qū)降水量，mm；F為裸露可溶巖面積，km2。

（2）滲漏補給。滲漏補給量采用河道退水系數法[4]計算：

式中：Q滲為區(qū)間河道滲漏量；m為河道滲漏系數；Q上為上游斷面實測流量；η為河道退水系數，η=Q上/Q下，其中Q下為下游斷面實測流量。

（3）巖溶水資源開發(fā)利用。晉祠泉域2012年有149口巖溶水井，取水量為1 605萬m3[5]。根據晉祠泉域巖溶水開發(fā)利用情況，1994—2012年巖溶地下水開采量見表1。

表1 晉祠泉域1994—2012年地下水位與各均衡要素匯總

（4）煤礦礦井排水量。1994年泉域排泄區(qū)巖溶水的礦井排水量為0.14 m3/s，隨著煤層開采深度的增加排水量逐漸增大，2006年達到最大值[2]。2008年山西省煤礦兼并重組后，礦井排水量逐年減小，見表1。

（5）側向排泄量分析。側向排泄量是泉域巖溶地下水循環(huán)的一個重要因素。根據晉祠泉域的水文地質條件，除了以泉水形式排泄的巖溶水外，還有部分巖溶水橫向排泄變?yōu)榕璧乜紫端推饺?。晉祠泉域側向排泄量與水力梯度和斷面滲透系數等有關，因此采用Darcy公式計算側向排泄量，結果見表1。

為了晉祠泉水復流和改善晉祠泉生態(tài)環(huán)境，從2013年開始，山西省委、省政府及有關部門積極采取了多種有效治理措施，諸如關停小煤礦和實施煤礦兼并重組，減少采礦排水；水源置換和關井壓采；汾河清水復流等。晉祠泉水復流工程實施前后地下水要素見表2。

表2 晉祠泉2013—2017年地下水位與各均衡要素匯總

2.2 泉水水位預測模型

基于前饋神經網絡Feed-forward-net、動態(tài)遞歸神經網絡Elman-net、時延神經網絡Time-delay-net、非線性動態(tài)自回歸神經網絡Narx-net、級聯神經網絡Cascade-net共5種人工神經網絡構建晉祠泉水位預測模型。

（1）前饋神經網絡。前饋神經網絡是一種方便、直觀、有效的基于誤差反向傳播算法的比較成熟的人工神經網絡模型。在這種網絡中，學習過程由正向和反向傳播組成。在正向傳播過程中，輸入信號從輸入層經隱含層單元傳向輸出層，每一層神經元的狀態(tài)只影響下一層神經元的狀態(tài)。若在輸出層不能得到期望輸出時，則轉入反向傳播，將輸出信號的誤差沿原來的連接通路返回，通過修改各層神經元的權值和閾值，使得網絡全局誤差最小[6]。常見的3層前饋神經網絡模型結構見圖1（wij為輸入層第i個神經元與隱含層第j個神經元之間的連接權值，vjt為隱含層第j個神經元與輸出層第t個神經元之間的連接權值）。

圖1 前饋神經網絡結構

設隱含層第j個神經元的閾值為θj，輸出層第t個神經元的閾值為γt，則隱含層第j個神經元的輸入sj、輸出bj為

式中：xi為輸入層第i個神經元的輸入；n為輸入層神經元個數；f（）為模擬神經元非線性特性的激勵函數。

輸出層第t個神經元的輸入Lt、輸出Ct為

式中：p為隱含層神經元個數。

（2）動態(tài)遞歸神經網絡。動態(tài)遞歸神經網絡作為一階時滯算子，達到了記憶的目的[7]，因此具有適應時變特性的能力，增強了網絡和網絡本身的全局穩(wěn)定性。其處理動態(tài)信息的能力強于前饋神經網絡，可以用來解決快速優(yōu)化的問題。動態(tài)遞歸神經網絡結構見圖2。

圖2 動態(tài)遞歸神經網絡結構

動態(tài)遞歸神經網絡的數學表達式：

式中：y（t）為 m 維輸出節(jié)點向量；x（t）為 n維中間層節(jié)點單元向量；u（t）為 r維輸入向量；xζ（t）為 n 維反饋狀態(tài)向量；W3為隱含層到輸出層的連接權值；W2為輸入層到隱含層的連接權值；W1為承接層到隱含層的連接權值；g[]為輸出神經元的傳遞函數，是隱含層輸出的線性組合；f｛｝為隱含層神經元的傳遞函數。

（3）時延神經網絡。時延神經網絡在結構上與普通神經網絡相似，但是每個神經元可以處理一個輸入子集，并且對于具有不同時滯的輸入數據具有不同的權重[8]。時延神經元TDN是時延神經網絡的基本單元。 每一個 TDN 有 M 個輸入 I1（t）、I2（t）、…、IM（t）和一個輸出O（t）。每個輸入Ii（t）包含N步時延，用以保存當前時刻之前的N個時刻的輸入Ii（t-d）（d=1，

式中：bi為第i個輸入的偏置；f｛｝為激勵函數。

時延神經網絡結構見圖3，其中隱含層有J個神經元，輸出層有R個神經元。輸出層各個神經元輸出為

式中：N1為輸出層節(jié)點的時延步數；r為輸出層節(jié)點編號；為第r個輸出層節(jié)點與隱含層節(jié)點間的連接權重；為第r個輸出層節(jié)點的偏置；（t）為第 j個隱含層的輸出；N2為隱含層節(jié)點的時延步數；j為隱含層節(jié)點編號；Ii（t-d）為第 i個輸入層節(jié)點；為第j個隱含層節(jié)點與輸入層節(jié)點的連接權重；為第j個隱含層節(jié)點的偏置。

圖3 時延神經網絡結構

（4）非線性動態(tài)自回歸神經網絡。Narx神經網絡將前饋神經網絡的輸出向量引入到輸入向量中，在保持時延后通過外部反饋。本文采用Narx神經網絡代替三層前饋神經網絡進行長期預測，給出了一種非線性動態(tài)自回歸神經網絡，Narx神經網絡的結構見圖4。其中：R=nu+ny，表示兩層Narx神經網絡的輸入個數，nu為非線性系統的輸入階數，ny為非線性系統的輸出階數，nu、ny的值可以通過模型階次辨識方法估計得出；j（圖4 中 i=1）為 Narx神經網絡隱含層中第 i個神經元和輸入向量中第j個元素之間的連接權值；S1為Narx神經網絡隱含層中神經元的個數；為Narx神經網絡隱含層中第i個神經元的偏置；為Narx神經網絡隱含層中第i個神經元的凈輸入；f1為Narx神經網絡隱含層神經元的傳輸函數；為Narx神經網絡隱含層中第i個神經元的輸出值；為 Narx 神經網絡的輸出神經元和隱含層中第i個神經元之間的連接權值；為Narx神經網絡輸出神經元的偏置；為Narx神經網絡輸出神經元的凈輸入；f2為Narx神經網絡輸出神經元的傳輸函數，在此使用純線性函數；既表示Narx神經網絡輸出神經元的輸出值，又表示Narx神經網絡在k+1時刻的預測輸出值ym（k+1）。

圖4 非線性動態(tài)自回歸神經網絡結構

（5）級聯神經網絡。在分類和預測方法的基礎上，選擇BP神經網絡作為級聯神經網絡的前級網絡，RBF徑向基網絡作為后級網絡。充分利用BP神經網絡學習時間短、RBF網絡訓練速度快的優(yōu)點，建立了基于BP神經網絡和RBF網絡的級聯神經網絡預測模型。將自然因素作為BP神經網絡的輸入變量，利用自然因素預測BP神經網絡的輸出。以BP神經網絡的輸出和人類因素為輸入，通過RBF網絡輸出自然和人為聯合作用下的泉水位預測值。該方法簡化了BP神經網絡和RBF網絡的拓撲結構，同時可以取得滿意的預測精度和計算效率。因篇幅有限，不再贅述RBF網絡的原理和結構，級聯神經網絡結構見圖5。

圖5 級聯神經網絡結構

2.3 降水量預測模型

LSTM深度學習模型是Hochreiter＆Schmidhuber為了避免長依賴問題結合適當的梯度學習算法設計的一種特殊的循環(huán)神經網絡（RNN），該模型具有刪除或添加信息到元胞（cell）狀態(tài)的門（Gate），可以學習長期依賴的歷史信息。LSTM深度學習模型與RNN模型最主要的區(qū)別在于其算法中加入了一個判斷信息有用與否的元胞結構，而LSTM的關鍵正是該元胞狀態(tài)——元胞的狀態(tài)類似于輸送帶，元胞的狀態(tài)在整個鏈上運行，只有一些小的線性操作作用其上，信息很容易保持不變的流過整個鏈，即激活函數的結構均為線性的，LSTM神經網絡只是對其進行傳遞，沒有做復雜的變換，不會更改信息。LSTM深度學習模型的每個判斷元胞當中放置了4扇閥門，分別為輸入閥門、遺忘閥門、更新閥門和輸出閥門，這些閥門可以打開或關閉，用于判斷模型網絡的記憶態(tài)在該層輸出的結果是否達到閾值從而加入到當前該層的計算中。當一個輸入信息進入LSTM的判斷元胞當中，該元胞即可根據既定規(guī)則來判斷該值是否有用。只有符合算法認證的信息才會留下并被記憶，不符合的信息則通過遺忘閥門被遺忘，從而有效地解決循環(huán)神經網絡（RNN）長期依賴的問題。LSTM深度學習模型結構見圖6。

圖6 LSTM深度學習模型元胞結構示意

（1）遺忘閥門Ft。遺忘閥門通過sigmoid激活函數讀取舊元胞的輸出值Yt-1和當前元胞的輸入值Xt，并結合偏置矢量和相應的權重來判斷是否保留當前元胞的狀態(tài) St，即“1”表示“保留“，“0”表示“刪除”。 該閥門的計算原理：

式中：σ為sigmoid激活函數；Wf為遺忘閥門的權重；Yt-1為上一個元胞的輸出值；Xt為當前元胞的輸入值；bf為遺忘閥門的偏置矢量。

（2）輸入閥門It。輸入閥門用來判斷輸入值Xt是否應該用于更新當前元胞的狀態(tài)St。輸入閥門一方面由sigmoid激活函數負責決定當前元胞的輸出值Yt，另一方面由tanh激活函數負責創(chuàng)建一個當前元胞的候選矢量Ut，該矢量由更新閥門判斷是否應將其加入到當前元胞的狀態(tài)St中。該閥門的計算原理：

式中：Wt為輸入閥門的權重；Wc為候選矢量的權重；bt為輸入閥門的偏置矢量；bc為候選矢量的偏置矢量；tanh表示雙曲正切激活函數。

（3）更新閥門Ut。該閥門通過將遺忘閥門和輸入閥門相乘來判斷當前元胞是否應該將舊元胞的狀態(tài)由St-1更新到St。該閥門的原理：

式中：St-1為舊元胞的狀態(tài)。

（4）輸出閥門Ot。輸出閥門負責計算當前元胞的輸出值Yt。輸出閥門通過運行sigmoid激活函數，判斷舊元胞的輸出值和當前元胞的輸入值來確定是否輸出當前的元胞狀態(tài)。接著把當前元胞狀態(tài)St通過tanh激活函數進行處理，并將它和sigmoid激活函數的輸出Ot相乘，最終確定當前元胞的輸出值Yt。輸出閥門的原理：

式中：Wo為輸出閥門的權重；bo為輸出閥門的偏置矢量。

2.4 模型訓練算法

晉祠泉水位預測模型訓練算法選擇了traingd（梯度下降法）、traingdm（有動量的梯度下降法）、traingda（自適應lr梯度下降法）、traingdx（自適應lr動量梯度下降法）、trainrp（彈性梯度下降法）、trainbr（貝葉斯正則化算法）、traincgf（Fletcher-Reeves共軛梯度法）、traincgp（Ploak-Ribiere共軛梯度法）、traincgb（Powell-Beale共軛梯度法）、trainscg（量化共軛梯度法）、trainbfg（擬牛頓算法）、trainosss（一步正割算法）、trainlm（Levenberg-Marquardt萊文貝格-馬夸特算法）、trainr（學習功能的隨機順序增量算法）共14種算法對5種人工神經網絡進行訓練。

2.5 模型精度評價

晉祠泉水位模型精度采用均方根誤差（RMSE）、效率系數（NSE）、平均絕對誤差（MAE）、確定系數（R2）、整體評價指標（Gpi）[9]進行評價，其計算方法：

式中：Xi、Yi分別為第i個時段泉水位的實測值和估算值；X′、Y′分別為泉水位實測值和估算值的平均值；n為時段數；Ti為 RMSE、MAE、R2、NSE 歸一化后的值；T′i為對應參數的中位數；當Tj為RMSE和 MAE時，αj取-1，其他情況取1。

3 結果與討論

3.1 訓練算法效率

通過5種神經網絡類型和14種訓練算法組合，建立了70種不同的網絡模型來預測晉祠泉水位。通過分析發(fā)現：在晉祠泉水位估算中，traincgb、trainrp、traincgf、traincgp等訓練算法效果比較理想；trainbfg、traingd、trainlm、trainr等算法無法模擬得到晉祠泉水位；其他算法大多能有效估算晉祠泉水位，但存在訓練次數多、隱含層節(jié)點數多、占用CPU大等問題，見表3。

表3 神經網絡類型、訓練算法運行效率比較

3.2 模型精度

基于5種人工神經網絡構建了晉祠泉水位預測模型，采用1994—2012年實測降水入滲補給量、河道滲漏補給量、巖溶水開采量、煤礦礦井排水量、側向排泄量和巖溶水位數據為訓練集對神經網絡模型進行訓練，以2013—2017年為神經網絡模型測試集驗證模型的預測精度，預測精度指標包括效率系數NSE（該值越接近1，表明模型可信度越高）、均方根誤差RMSE（該值越小，表明模型偏差越小）、平均絕對誤差MAE（該值越小，表明模型精度越好）、確定系數R2（該值越接近1，表明模型吻合度越高）。為了統一模型精度評價體系，引入整體評價指標（Gpi）評價模型適用性，結果見表4、表5。

表4 人工神經網絡泉水位模擬結果預測模型適應性評價

表5 人工神經網絡泉水位預測模型適應性評價

由表5可以看出，本文選用的5種人工神經網絡模型的各項適應性評價指標差異明顯：NSE、RMSE值動態(tài)遞歸神經網絡模型最優(yōu)，前饋神經網絡模型最差；MAE值級聯神經網絡模型最優(yōu)，前饋神經網絡模型最差；R2值前饋神經網絡模型最優(yōu)，非線性動態(tài)自回歸神經網絡模型最差。模型的適應性排名并不一致，無法單靠某一指標判斷模型的精度和適應性，從而選擇最合適的泉水位預測模型。因此引入了整體評價指標（Gpi）來綜合評價模型的適用性，5種人工神經網絡模型的Gpi值動態(tài)遞歸神經網絡模型＞級聯神經網絡模型＞時延神經網絡模型＞非線性動態(tài)自回歸神經網絡模型＞前饋神經網絡模型，就此可以判斷動態(tài)遞歸神經網絡模型為預測晉祠泉水位變化的最合適模型。

由圖7可以看出，晉祠泉在1994—2012年的近20 a間水位呈現谷狀曲線，在2008年左右晉祠泉水位達到最低值774.94 m后，2009—2012年呈現逐漸上升趨勢，最高值達到788.21 m。降水入滲量在1994—2008年雖然變化幅度不大但仍呈現整體上升的態(tài)勢，這與晉祠泉水位下降的趨勢相反，可見降水入滲并非影響晉祠泉水位變化的關鍵因素。同時，晉祠泉開采量呈現1994—2008年逐年增加繼而到2012年逐年減少的態(tài)勢，這與晉祠泉水位變化一致，可以判斷出晉祠泉開采是晉祠泉水位變化的關鍵影響因素。

圖7 晉祠泉域水位人工神經網絡模型訓練、驗證及預測

2013年后，為晉祠泉水復流和改善晉祠泉生態(tài)環(huán)境，山西省政府采取了多種積極有效的措施，即使2013—2017年降水量下降明顯，晉祠泉水位仍然保持了逐年上升的良好態(tài)勢，這也印證了晉祠泉開采等人類因素才是影響泉水位變化的關鍵因素。

3.3 未來10 a晉祠泉水位預測

為了預測2018—2027年晉祠泉域地下水位變化趨勢，采用動態(tài)遞歸神經網絡模型預測晉祠泉域2018—2027年地下水位的變化趨勢并測定了晉祠泉的復流時間（見圖7）。其中基于晉祠泉域1956—2017年共62 a的降水系列資料，采用LSTM人工神經網絡模型預報未來2018—2027年（未來10 a）晉祠泉域的降水量（見圖7），并根據所預測的降水量采用式（1）計算降水入滲補給量，采用式（2）計算河道滲漏補給量，采用Darcy公式計算側向排泄量，巖溶水資源開發(fā)利用量和煤礦礦井排水量則采用《太原市晉祠泉域水資源保護條例》中的限定值。

從圖7可以看出，由于人類積極影響因素的介入，晉祠泉水位保持平穩(wěn)上升的良好態(tài)勢，泉水位在2019年可上升到802.66 m，超過了晉祠泉復流的最低水位802.59 m[10]，泉水位在未來10 a最高可上升至803.74 m，表明多種積極有效措施的實施可以取得良好的泉水復流效果。

4 結 論

人類因素是影響晉祠泉水位的關鍵因素，人類消極干預嚴重是造成晉祠泉水位下降、斷流的主要原因。在人類積極干預介入下，泉水位完全具備復流的條件，晉祠泉在2018年左右可以復流。各種訓練算法并非均能有效適用于晉祠泉水位預測，其中 trainbfg、traingd、trainlm、trainr等算法無法模擬晉祠泉水位，traincgb、trainrp、traincgf、traincgp 等訓練算法效果較好，其他算法大多存在訓練次數多、占用CPU大等問題。動態(tài)遞歸神經網絡模型為預測晉祠泉水位變化的最合適模型。