汪 云，楊海博，徐 建，鄭夢琪，韓智昕，趙 耘，趙 耀

(中國冶金地質(zhì)總局山東正元地質(zhì)勘查院，山東 泰安 271000)

0 引 言

地下水是地球水資源的重要組成部分，是農(nóng)業(yè)灌溉、工礦業(yè)和城市的重要水源之一。但在一定的條件下，比如降雨、蒸發(fā)、徑流和人工開采等，都會造成地下水位的變化，從而有可能產(chǎn)生地下水降落漏斗、地面沉降等環(huán)境地質(zhì)問題，對生態(tài)和人類生產(chǎn)生活造成重大影響。所以，對地下水的動態(tài)變化趨勢進行研究具有相當(dāng)重要的意義。

地下水位的變化是一個復(fù)雜的自然過程，通過預(yù)測地下水位，在一定程度上可以反映出未來地下水的動態(tài)變化規(guī)律。目前，主要采用數(shù)學(xué)模型來進行地下水位的預(yù)測，應(yīng)用于地下水位預(yù)測的數(shù)學(xué)模型可分為確定性模型[1]和隨機性模型[2]。確定性模型通常應(yīng)用有限的物理學(xué)規(guī)律來描述水文過程，是指不包含任何隨機成分的模型，只要設(shè)定了輸入和各個輸入之間的關(guān)系，其輸出也是確定的，該方法需要長期的水文氣象資料，對資料數(shù)據(jù)的數(shù)量和精度要求較高，過程較為復(fù)雜，計算量較大，模型參數(shù)的優(yōu)選和識別有一定難度。隨機性模型應(yīng)用隨機過程描述水文環(huán)節(jié)，包括回歸分析、時間序列分析等，還有模糊識別法、灰色模型、小波神經(jīng)網(wǎng)絡(luò)分析等新興的方法。張展羽等[3]針對地下水位在時間序列上表現(xiàn)出高度的隨機性和滯后性，建立了基于主成分分析與多變量時間序列模型耦合的地下水位預(yù)報型，將該模應(yīng)用于濟南市陡溝灌區(qū)地下水位預(yù)測。陸海濤和梁富山[4]將支持向量機應(yīng)用于地下水位預(yù)測中，并與GM(1，1)模型對比。徐強等[5]考慮到BP算法存在極易收斂于局部極小點與過擬合等缺點，構(gòu)建小波神經(jīng)網(wǎng)絡(luò)基礎(chǔ)上并引入遺傳算法加以優(yōu)化，以解決上述不足，并與BP和WNN對比預(yù)測了天津市深層承壓水水位。李慧等[6]針對不恰當(dāng)?shù)剡x取RBF神經(jīng)網(wǎng)絡(luò)的網(wǎng)絡(luò)結(jié)構(gòu)和參數(shù)會使網(wǎng)絡(luò)收斂慢的問題，采用粒子群優(yōu)化算法對RBF神經(jīng)網(wǎng)絡(luò)參數(shù)進行優(yōu)化，建立了基于粒子群優(yōu)化算法的RBF神經(jīng)網(wǎng)絡(luò)模型(POS-RBF模型)，對涇惠渠灌區(qū)地下水位埋深進行了模擬和預(yù)測。徐源蔚等[7]針對地下水位樣本之間的相似性及影響因子與地下水位之間具有的確定、不確定性特征，提出將集對分析與相似預(yù)測相結(jié)合，從同、異、反三方面定量刻畫地下水位的當(dāng)前樣本與歷史樣本之間的相似性，建立了基于集對分析的地下水位相似預(yù)測模型。周振民等[8]對GM(1，1)模型的結(jié)構(gòu)進行了研究和改進，選用灤河下游地區(qū)節(jié)水工程改造后地下水位資料，對灤河下游地區(qū)地下水動態(tài)進行分析和數(shù)值模擬。這些預(yù)測模型各具優(yōu)點，無論是傳統(tǒng)方式還是新模型，在實際應(yīng)用中均取得了一定的成果。但實踐表明，單一的預(yù)測方式在模擬地下水位復(fù)雜的動態(tài)特征方面仍舊不可避免地存在著一定的局限性。

神經(jīng)網(wǎng)絡(luò)是20世紀80年代以來人工智能領(lǐng)域興起的研究熱點，相比于傳統(tǒng)的統(tǒng)計分析模型，具有更好的實時預(yù)測性，能夠解決多個自變量和因變量的問題，目前廣泛應(yīng)用于地下水位預(yù)測中。常見的有BP神經(jīng)網(wǎng)絡(luò)模型[9-11]，RBF神經(jīng)網(wǎng)絡(luò)模型等[12]。近些年，隨著深度學(xué)習(xí)技術(shù)的高速發(fā)展，深度學(xué)習(xí)模型逐漸被應(yīng)用于時間序列數(shù)據(jù)的研究中。其中，循環(huán)神經(jīng)網(wǎng)絡(luò)將時間序列引入神經(jīng)網(wǎng)絡(luò)中，使模型在時序數(shù)據(jù)分析上表現(xiàn)出更強的適應(yīng)性。但是也存在梯度消失或者梯度爆炸的問題。長短期記憶神經(jīng)網(wǎng)絡(luò)模型通過遺忘門和更新門來調(diào)控梯度，將反向傳播累乘作用變成累加作用解決了梯度消失或爆炸的問題。本文將長短期記憶神經(jīng)網(wǎng)絡(luò)模型應(yīng)用于地下水水位預(yù)測研究，以期獲得一種更好的地下水位預(yù)測方法。

1 研究方法

1.1 RNN循環(huán)神經(jīng)網(wǎng)絡(luò)

傳統(tǒng)神經(jīng)網(wǎng)絡(luò)是通過輸入層，然后存在隱含層，再到輸出層，實現(xiàn)了層與層之間的鏈接，但是每一層的節(jié)點卻是相互獨立，沒有鏈接的，因此傳統(tǒng)神經(jīng)網(wǎng)絡(luò)在處理序列數(shù)據(jù)時，結(jié)果往往有很大的偏差。不同于傳統(tǒng)神經(jīng)網(wǎng)絡(luò)，循環(huán)神經(jīng)網(wǎng)絡(luò)會對之前的信息進行記憶并應(yīng)用到當(dāng)前的輸出計算中，即隱含層之間的節(jié)點為連接的，并且隱含層的輸入包括輸入層的輸出與上一時刻隱含層的輸出。RNN前向輸出如圖1所示。

圖1 RNN前向輸出流程注：x表示輸入量；h表示隱含層單元；o表示輸出量；L表示損失函數(shù)；y表示訓(xùn)練集標簽；V、W、U表示共享權(quán)重系數(shù)矩陣。

對于時刻t：

h(t)=f(Ux(t)+Wh(t-1)+b)

(1)

輸出量：

o(t)=Vh(t)+c

(2)

模型預(yù)測輸出：

(3)

式中：soft max為激活函數(shù)。

基于處理時序數(shù)據(jù)的特點，RNN常見訓(xùn)練方法為BPTT算法，訓(xùn)練尋優(yōu)的參數(shù)為V、W、U。但在尋優(yōu)過程中，可能存在“梯度消失”或“梯度爆炸”現(xiàn)象[13]。

1.2 LSTM神經(jīng)網(wǎng)絡(luò)

為解決RNN上述問題，Hochreiter和Schmidhuber于1997年提出了LSTM[14]，通過門控制將短期與長期記憶結(jié)合起來，從而在一定程度上解決梯度消失的問題，LSTM結(jié)構(gòu)如圖2所示。

圖2 LSTM輸出結(jié)構(gòu)

LSTM通過讀取h(t-1)與x(t)，輸出0或1的數(shù)值給c(t-1)中的數(shù)字，其中“1”表示保留，“0”表示丟棄，步驟如下：

ft=σ(W(f)·[h(t-1)，x(t)]+b(f))

(4)

式中：h(t-1)表示上一時刻的單元輸出；b表示對應(yīng)的偏置項。

輸入門決定著讓多少信息加入到單元中，通過sigmoid層的信息決定與tanh層的信息生成，聯(lián)合對單元狀態(tài)進行更新。輸出門決定著當(dāng)前單元狀態(tài)哪一部分信息被輸出，依然通過sigmoid層與tanh層完成。輸入門步驟如式(5)、(6)、(7)所示，輸出門步驟如式(8)、(9)所示[11]：

i(t)=σ(W(i)·[h(t-1),x(t)]+b(i))

(5)

(6)

(7)

o(t)=σ(W(o)·[h(t-1),x(t)+b(o)])

(8)

h(t)=o(t)·tanh(c(t))

(9)

式中：f、i、c、o分別表示遺忘門、輸入門、單元狀態(tài)、輸出門；h(t)表示當(dāng)前時刻的輸出；σ、tanh分別表示激活函數(shù)sigmoid與雙曲正切函數(shù)。

2 研究區(qū)概況與數(shù)據(jù)來源

2.1 研究區(qū)概況

泰安市地處山東省中部，地理坐標：東經(jīng)116°04′～118°00′，北緯35°38′～36°29′，總面積7 762 km2，耕地面積316 907 hm2。屬于華北暖溫帶半濕潤大陸性季風(fēng)氣候區(qū)，多年平均降水量為681.6 mm，多集中在6-9月份，占全年降水量的75%，多年平均蒸發(fā)量1 813 mm。境內(nèi)水系主要屬于大汶河水系和淮河水系，地形東高西低，為侵蝕或剝蝕構(gòu)造地形，海拔37.5～1 532 m；塊狀地貌是研究區(qū)地貌主要特征，斷塊山與斷陷盆地發(fā)育，按地貌成因，可劃分為侵蝕構(gòu)造中度切割中山，侵蝕、剝蝕、溶蝕低山，侵蝕、剝蝕、溶蝕丘陵，剝蝕堆積山前傾斜平原，山間盆地沖積平原和黃河沖積平原6種類型。

研究區(qū)地層由老到新分屬太古界侵入巖及變質(zhì)巖系(Ar、)，古生界寒武系(∈)、奧陶系(O)，上古生界石炭系(C)、二疊系(P)，中生界侏羅系(J)、白堊系(K)，新生界古近系(E)、第四系(Q)。區(qū)內(nèi)地質(zhì)構(gòu)造較為發(fā)育，有泰山和徂徠山復(fù)式背斜，泰山斷裂、新泰-蒙陰斷裂、蒙山斷裂、南留弧型斷裂、夏張斷裂、肥城弧型斷裂。研究區(qū)地下水類型可劃分為松散巖類孔隙水、碎屑巖類孔隙裂隙水、碳酸鹽巖類裂隙巖溶水和基巖裂隙水。研究的泰安市岱岳區(qū)滿莊鎮(zhèn)姜家園村046J地下水屬于松散巖類孔隙水，水位變化主要受氣象因素影響。

2.2 研究數(shù)據(jù)來源

地下水位數(shù)據(jù)為地處泰安市岱岳區(qū)滿莊鎮(zhèn)姜家園村046J地下水監(jiān)測井。選取2001-2016年該井逐月平均地下水位觀測資料和該地氣候數(shù)據(jù)作為研究數(shù)據(jù)，地下水位變化如圖3所示。

圖3 鎮(zhèn)姜家園村046J井歷年地下水位變化

3 地下水位預(yù)測模型構(gòu)建

3.1 模型框架

根據(jù)地下水水位的時空變化特點以及循環(huán)神經(jīng)網(wǎng)絡(luò)從簡設(shè)計的原則，本研究構(gòu)建的預(yù)測模型框架如圖4所示，該框架大致可分為輸入層、隱含層、模型訓(xùn)練、輸出層4個部分，各模塊功能如下。

圖4 LSTM預(yù)測模型結(jié)構(gòu)圖

3.2 數(shù)據(jù)預(yù)處理

為消除量綱及數(shù)量級的影響，提高模型預(yù)測精度，加快模型訓(xùn)練收斂速度，需要對原始變量時序集進行標準化處理，公式可表示為：

x*=(x-xmin)/(xmax-xmin)

(10)

式中：x*表示標準化后的值；x表示標準化前的值；xmax表示原始數(shù)據(jù)集中最大值；xmin表示原始數(shù)據(jù)集中最小值。

3.3 模型訓(xùn)練

根據(jù)LSTM模型結(jié)構(gòu)，網(wǎng)絡(luò)訓(xùn)練步驟如下：

(1)在輸入層中，定義原始變量數(shù)據(jù)集為：

Fn=(x1,x2,…,xt,…,xn)

(11)

式中：xt=(at,bt,ct,dt,et,ft,gt),at,bt,ct,dt,et,ft,gt分別表示t時刻對應(yīng)的地下水位、蒸發(fā)量、降雨量、氣溫、氣壓、相對濕度、日照時長。

(2)劃分訓(xùn)練集與測試集，并對原始變量數(shù)據(jù)集按照式(10)進行標準化處理，可表示為：

(12)

(13)

X={x1′,x2′,…,xL′}

(14)

X(t)={xt′,x(t+1)′,…,x(m-L+t-1)′} 1≤t≤L,t,L∈N

(15)

(4)隱含層理論輸出為：

Y={y1,y2,…,yL}

(16)

Y(t)={y(t+1)′,y(t+2)′,…,y(m-L+t)′}

(17)

實際輸出為：

P={p(1),p(2),…,p(L)}

(18)

P(t)=LSTM(x(t),c(t-1),h(t-1))

(19)

4 實例驗證

將上述訓(xùn)練好的模型進行地下水位預(yù)測，選取2001-2014年地下水位、蒸發(fā)量、降雨量、氣溫、氣壓、相對濕度、日照等數(shù)據(jù)作為訓(xùn)練集，2015-2016年數(shù)據(jù)作為預(yù)測集，模型預(yù)測精度以均方根誤差作為評價指標。

首先固定非關(guān)鍵參數(shù)取值，選擇批次大小為168，根據(jù)地下水水位變化特點，選擇時間步長為12，確定重復(fù)訓(xùn)練次數(shù)為1 000，其次確定重要參數(shù)的范圍以簡化調(diào)參過程，隱含層節(jié)點數(shù)經(jīng)驗公式如下式所示，學(xué)習(xí)率一般按3的倍數(shù)來調(diào)整。

(20)

(21)

式中：m為隱含層節(jié)點數(shù)；n為輸入層節(jié)點數(shù)；l為輸出層節(jié)點數(shù)，a∈[1,10]，a∈N。因此，隱含層節(jié)點數(shù)m∈[3,13]，學(xué)習(xí)率lr∈{0.001,0.003,0.01,0.03,0.1,0.3}。相關(guān)組合均方根誤差如表1所示。

表1 相關(guān)參數(shù)組合對應(yīng)誤差

由表1可得，當(dāng)lr=0.03，m=13時，均方根誤差最小，為1.436 5，選取該參數(shù)對模型訓(xùn)練，整個訓(xùn)練、預(yù)測過程在python環(huán)境中進行。

為驗證多變量LSTM神經(jīng)網(wǎng)絡(luò)預(yù)測地下水水位的優(yōu)點，本章節(jié)使用單變量LSTM神經(jīng)網(wǎng)絡(luò)、BP神經(jīng)網(wǎng)絡(luò)對地下水水位進行預(yù)測，并作對比分析，已有大量學(xué)者對BP神經(jīng)網(wǎng)絡(luò)做研究，本章不再贅述，其中，單變量LSTM神經(jīng)網(wǎng)絡(luò)僅以地下水水位作為變量輸入，三者訓(xùn)練期、預(yù)測期擬合曲線如圖5、6、7所示。

圖5 多變量LSTM預(yù)測結(jié)果

圖6 單變量LSTM預(yù)測結(jié)果

圖7 BP神經(jīng)網(wǎng)絡(luò)預(yù)測結(jié)果

由預(yù)測結(jié)果可知，多變量LSTM模型的預(yù)測精度高于單變量LSTM模型、BP神經(jīng)網(wǎng)絡(luò)模型，且擬合程度也較為理想，得益于多變量LSTM對時間序列訓(xùn)練的同時，引入了的相關(guān)影響變量。單變量LSTM模型僅能對時間序列訓(xùn)練，當(dāng)外界條件有較大變動時，預(yù)測精度將受到較大影響，因此單變量LSTM模型的擬合、預(yù)測曲線呈現(xiàn)出較為平穩(wěn)的周期性變化。BP神經(jīng)網(wǎng)絡(luò)預(yù)測結(jié)果往往取決于經(jīng)驗值，因此在訓(xùn)練期擬合曲線趨勢大致與歷史相近，但對于未來水位的預(yù)測誤差很大，究其原因，主要是歷史訓(xùn)練數(shù)據(jù)過少且無法學(xué)習(xí)歷史時序規(guī)律?；诙嘧兞康腖STM預(yù)測模型擬合曲線并不完全完美，在某些峰值處誤差較大，特別是在地下水水位變化規(guī)律出現(xiàn)波動時，但總體效果較為穩(wěn)定。分析地下水水位變化曲線可知，研究區(qū)的地下水水位變化規(guī)律性較差，剖析其原因，最大影響因素為降雨量的變化波動較大，其他影響變量波動則較為平緩，地下水水位與降雨量的動態(tài)變化如圖8所示。分析兩者動態(tài)變化規(guī)律可知，研究區(qū)地下水水位的變化往往滯后于降雨量的變化，且降雨量變化規(guī)律與水位變化規(guī)律基本一致，充分說明了降雨量是影響該地區(qū)地下水水位的重要因素。曲線中2010-2016年降雨量略小于過去幾年，地下水水位卻高于往年，分析其原因可能是由于人工開采量的減少導(dǎo)致水位的上升，在沒有開采量數(shù)據(jù)的情況下，預(yù)測精度依然較高則是因為LSTM預(yù)測時，考慮周邊影響變量的同時，也考慮了歷史數(shù)據(jù)的時空變化規(guī)律，也說明了研究區(qū)近年地下水開采量變化較為穩(wěn)定，沒有發(fā)生突變。

圖8 地下水水位與降雨量動態(tài)變化

5 結(jié) 語

通過構(gòu)建基于長短期記憶神經(jīng)網(wǎng)絡(luò)(LSTM)的地下水水位預(yù)測模型并與BP神經(jīng)網(wǎng)絡(luò)預(yù)測模型與單變量LSTM預(yù)測模型對比，得出以下結(jié)論。

(1)基于多變量的LSTM神經(jīng)網(wǎng)絡(luò)模型可以作為預(yù)測地下水水位動態(tài)變化簡單有效的工具，特別是在一些水文地質(zhì)資料比較匱乏的地區(qū)。

(2)在少量歷史數(shù)據(jù)的情況下，基于多變量的LSTM預(yù)測模型亦能做到有效的預(yù)測，而且對歷史數(shù)據(jù)重模擬的擬合效果較為理想。單變量LSTM預(yù)測模型僅僅呈現(xiàn)出周期性變化，同時存在較為嚴重的滯后性，無法為一些外界條件變化較大的地區(qū)提供水資源管理參考。BP神經(jīng)網(wǎng)絡(luò)對歷史數(shù)據(jù)學(xué)習(xí)性比較強大，但對于未來地下水水位的預(yù)測則需要大量的歷史數(shù)據(jù)訓(xùn)練作為支撐。

(3)地下水水位動態(tài)變化與氣候、下墊面等自然因素相關(guān)的同時，很大程度上還受人工開采量的影響，若在缺乏開采量資料的地區(qū)進行預(yù)測，且研究區(qū)開采量波動較大時，就會影響模型原有的預(yù)測模式，使預(yù)測結(jié)果產(chǎn)生較大偏差。