基于長短時記憶神經(jīng)網(wǎng)絡的鄱陽湖水位預測*

郭 燕,賴錫軍

(1:中國科學院南京地理與湖泊研究所，中國科學院流域地理學重點實驗室, 南京 210008) (2:中國科學院大學資源與環(huán)境學院，北京 100049)

鄱陽湖是中國第一大淡水湖，水量主要來自流域“五河”. 汛期，長江水可能倒灌鄱陽湖. 鄱陽湖是典型的季節(jié)性吞吐型湖泊，水位年內年際變化劇烈，具有“洪水一片，枯水一線”的特征. 每年4-6月為鄱陽湖主汛期，湖水位隨“五河”來水的增加而升高；7-9月是長江主汛期，受到長江高水位的頂托作用，湖水位進一步升高；每年10月-次年3月鄱陽湖屬于枯水期，水位較低. 水位變化會導致濕地植被分布以及候鳥棲息等也相應發(fā)生改變[1-3]. 近十幾年來，在氣候變化和高強度的人類活動影響下，鄱陽湖水情發(fā)生了劇烈的變化，特別是枯水問題突出，出現(xiàn)了枯水發(fā)生時間提前、持續(xù)時間延長的現(xiàn)象. 湖區(qū)生態(tài)環(huán)境因水量平衡的變化面臨威脅[4-5].

鄱陽湖水量平衡近年成為了學界關注熱點. 為闡明鄱陽湖的水量變化，開展了大量卓有成效的研究，建立了諸多水文水動力模型和數(shù)理統(tǒng)計模型，分析了氣候變化和不同類型的人類活動對鄱陽湖水量的影響過程以及江湖交互作用下湖泊水位的變化特征[6-9]. 水動力模型可實現(xiàn)過程的精細模擬，但是水動力系統(tǒng)龐大、需要非常完善的基礎資料、計算復雜、耗時較長，對水位高效預測仍存挑戰(zhàn)[10-11]. 數(shù)理統(tǒng)計方法(如：多元線性回歸(MLR)、累積距平曲線(CPC)、概率分布函數(shù)(PDF)、泰爾森方法(TSA))簡單易行，主要用于長時間尺度的變化檢驗分析[12-16]. 鄱陽湖水位受流域“五河”和長江來水多重控制，其變化與“五河”來水和長江來水有著非常復雜的非線性關系[17]，難以采用簡單的統(tǒng)計模型模擬預測水位和來水的響應關系. 人工神經(jīng)網(wǎng)絡方法作為一種數(shù)據(jù)驅動的自適應性方法，沒有先驗假設，可較好地模擬預測復雜的非線性作用關系，也在鄱陽湖得到了應用[18]. 近年來，基于人工神經(jīng)網(wǎng)絡的深度學習算法得到了飛速發(fā)展. 特別是長短時記憶神經(jīng)網(wǎng)絡方法(long short-term memory,LSTM)，它通過循環(huán)反饋結構存儲歷史信息，具有較強的時間序列問題求解能力，在模擬預測水文時間序列問題中受到關注. Zhang等[19]基于雨量器和水位傳感器的在線數(shù)據(jù)，比較了不同神經(jīng)網(wǎng)絡結構在模擬和預測挪威Drammen聯(lián)合下水道結構水位方面的性能，證實了LSTM比沒有顯式細胞記憶的傳統(tǒng)架構更適合于多步預測. Zhang等[20]使用LSTM預測農(nóng)業(yè)地區(qū)的地下水位，并將基于LSTM方法的預測結果與傳統(tǒng)神經(jīng)網(wǎng)絡的預測結果進行了比較，發(fā)現(xiàn)前者的性能優(yōu)于后者. Lee等[21]利用基于物理的水文模型SWAT和數(shù)據(jù)驅動的深度學習算法LSTM，在湄公河下游Kratie站進行了徑流模擬. Liu等[22]提出了改進的CA-LSTM上下文感知神經(jīng)網(wǎng)絡，基于收集到的洪水因子序列數(shù)據(jù)對洪水進行了預測.

鄱陽湖水位受江湖共同制約，其變化與流域“五河”來水量及長江來水量具有非常復雜的非線性關系. 本文選取鄱陽湖為研究對象，采用LSTM神經(jīng)網(wǎng)絡方法，建立鄱陽湖水位日尺度的預測模型，探討其用于預測江湖交互作用下鄱陽湖水位變化過程的潛力，為鄱陽湖水量平衡研究提供一個快速有效的預測方法.

1 數(shù)據(jù)和方法

1.1 研究區(qū)概況

鄱陽湖(28°24′～29°46′N，115°49′～116°46′E)位于長江中下游南岸，北部較平坦，東西南三面環(huán)山，全流域總體地勢是南高北低，南部地形起伏變化較大，形成了一個以鄱陽湖為底，向北開口連接長江的巨大盆地. 鄱陽湖流域由贛江、撫河、信江、饒河和修水5個相互獨立的子流域組成，與湖區(qū)共同形成了一個完整的流域系統(tǒng)，流域總面積為16.2×104km2[16, 23]. 多年平均徑流量為1491億m3，水位年內變幅為9.59～15.36 m.

1.2 研究數(shù)據(jù)

選用“五河”流域外洲站、李家渡站、梅港站、渡峰坑站、虎山站、萬家埠站和長江中游漢口站7站逐日平均流量數(shù)據(jù)以及湖口站、星子站、都昌站、吳城站和康山站5站逐日平均水位數(shù)據(jù)建模，時間序列長度均為1956-2000年,監(jiān)測站點分布見圖1.

圖1 鄱陽湖監(jiān)測站點的空間分布Fig.1 Location of sampling sites in the Lake Poyang

1.3 LSTM基本原理

1.3.1 LSTM前向計算 LSTM是1997年Hochreiter和Schmidhuber[24]為了解決RNN模型的梯度消失或爆炸缺陷而開發(fā)的一種復雜的遞歸模型. Gers等[25]正式提出LSTM網(wǎng)絡層是由遺忘門f、輸入門i、記憶單元C和輸出門o組成(圖2). LSTM的關鍵是記憶單元，它像傳送帶一樣，將信息從上一個單元傳遞到下一個單元.

圖2 LSTM內部結構示意圖Fig.2 Schematic diagram of the internal structure of LSTM

首先，決定從記憶單元中增減多少信息. 遺忘門本質上是一個σ神經(jīng)網(wǎng)絡層，根據(jù)前一時刻的輸出ht-1和當前的輸入Xt，產(chǎn)生一個介于0和1之間的ft值：

ft=σ(Wxf·Xt+Whf·ht-1+bf)

(1)

接下來，確定將在記憶單元中存儲哪些新信息. 這一部分由一個σ神經(jīng)網(wǎng)絡層和一個tanh神經(jīng)網(wǎng)絡層兩部分組成. 輸入門根據(jù)前一時刻的輸出ht-1和當前的輸入Xt，產(chǎn)生一個介于0和1之間的it值：

it=σ(Wxi·Xt+Whi·ht-1+bi)

(2)

(3)

從而更新記憶單元：

(4)

最后，我們決定輸出什么信息. 首先，運行一個σ神經(jīng)網(wǎng)絡層來決定輸出記憶單元的哪些部分. 然后，通過tanh神經(jīng)網(wǎng)絡層將Ct值調整為-1～1，并將其乘以σ神經(jīng)網(wǎng)絡層的輸出，最終得到目標信息.

ot=σ(Wxo·Xt+Who·ht-1+bo)

(5)

ht=ot·tanh(Ct)

(6)

式(1～6)中，f、i、o、C和h分別是遺忘門、輸入門、輸出門、記憶單元和輸出信息，W是相應的權重矩陣，b是偏差矩陣，σ和tanh是激活函數(shù).

1.3.2 LSTM反向訓練 神經(jīng)網(wǎng)絡訓練的過程是去尋找最優(yōu)參數(shù)，使得模型收斂，即損失函數(shù)達到極小甚至最小的過程. 網(wǎng)絡通過反向傳播損失函數(shù)，利用梯度下降法迭代更新網(wǎng)絡的權重[26-27]. 我們選取最常用的均方誤差計算損失函數(shù)：

(7)

1.4 模型評估標準

LSTM是通過已有的訓練樣本(即已知數(shù)據(jù)及對應的輸出)去學習得到一個最優(yōu)模型，再利用該模型將所有的輸入映射為相應的輸出. 我們采用均方根誤差(RMSE)、納什效率系數(shù)(NSE)和決定系數(shù)(R2)3個指標對模型模擬的效果進行評價[28].

(8)

RMSE值為0表示觀測值與預測值完全吻合.

(9)

NSE是一個參數(shù)，它決定了剩余方差(噪聲)相對于測量數(shù)據(jù)(信息)中的方差的相對重要性，取值范圍為負無窮至1. 接近1，表示模擬效果好，模型可信度高；接近0，表示模擬結果接近觀測值的平均值水平，即總體結果可信，但過程模擬誤差大.

(10)

2 鄱陽湖水位預報模型構建

2.1 模型結構

采用單層LSTM循環(huán)神經(jīng)網(wǎng)絡建模，輸入長江和“五河”來水量，其中長江水量以漢口站流量為代表，“五河”水量以贛江外洲站、撫河李家渡站、信江梅港站、饒河渡峰坑站和虎山站、修水萬家埠站流量為代表，輸出湖口站、星子站、都昌站、吳城站和康山站5個代表站的水位，實現(xiàn)根據(jù)“五河”來水量和長江來水量組合預測鄱陽湖湖區(qū)未來1 d水位的空間分布，神經(jīng)網(wǎng)絡結構見圖3，表達式為：

(11)

圖3 單層長短時記憶網(wǎng)絡LSTM模型結構Fig.3 Structure of a single layer LSTM model

數(shù)據(jù)之間的差異性會對模型的學習能力產(chǎn)生負面影響. 因此，為了保證構建的模型中參數(shù)能夠穩(wěn)定收斂，在神經(jīng)網(wǎng)絡訓練之前需對其進行預處理，以確保所有變量保持在相同的量表上[29-30]. 我們對所有數(shù)據(jù)進行[0,1]歸一化處理：

(12)

式中，X是歸一化后的數(shù)據(jù)，Xori是原始數(shù)據(jù)，Xmax和Xmin分別是原始數(shù)據(jù)的最大值和最小值.

模型損失函數(shù)選用均方根誤差，優(yōu)化選取基于梯度下降的ADAM算法，采用Dropout正則化方法防止模型過擬合[31-33]. 訓練集為1956-1980年的水文時間序列數(shù)據(jù)，驗證集數(shù)據(jù)長度為1981-2000年. 在LSTM神經(jīng)網(wǎng)絡中，輸入時間窗、隱藏神經(jīng)元數(shù)目、初始學習率大小等重要參數(shù)會直接影響到模型預測效果. 所以，接下來我們將對模型中的這些參數(shù)進行優(yōu)選.

2.2 模型參數(shù)優(yōu)選

2.2.1 輸入時間窗長短對模型模擬效果的影響 鄱陽湖作為大型通江洪泛湖泊，其入湖流量與湖泊水位關系呈現(xiàn)明顯的非線性特征，加之長江與湖泊之間相互作用的復雜性，湖泊流量與水位之間呈現(xiàn)多種對應關系[34]. “五河”水量經(jīng)一段時間傳輸入湖，鄱陽湖接受“五河”及長江來水進行調蓄作用，因而水位和流量并不總是同步變化，存在明顯的相位滯后效應，輸入時間窗的優(yōu)選計算就是確定流量-水位之間滯后時長的過程. 我們計算并比較了不同長短的輸入時間窗(1～10 d，即分別將當前時刻的流量、前2 d內的流量、……、前10 d內的流量作為輸入變量)下，根據(jù)“五河”來和長江來水量，組合預測鄱陽湖湖區(qū)未來1 d的水位空間分布. 同時，為了確保模型較高的準確率，隱藏神經(jīng)元數(shù)目應盡量多取，設置了100個；初始學習率大小應定義為較小的值，設置為0.001. 不同輸入時間窗模式下模型擬合評價指標計算結果如圖4所示.

圖4 不同輸入時間窗下模型訓練和驗證過程的模擬結果Fig.4 Simulation results of model training and verification under different input time windows

雖然神經(jīng)網(wǎng)絡模型的外延效果一般不是很理想，但可通過反復對模型調試計算、選取最適宜的輸入時間窗來提高模型預測效果. 圖4結果表明：不同輸入時間窗下，湖口站水位預測的效果均最好，其次分別是星子站、都昌站和吳城站，康山站水位預測效果為五站中最差. 利用當前時刻的七站流量來預測未來1 d的鄱陽湖五站水位時，湖口站訓練和驗證過程的RMSE均較大，分別為0.91 m和0.79 m.NSE和R2相對較小，訓練階段NSE為0.94，R2為0.95；驗證階段NSE和R2均為0.96. 隨著輸入時間窗逐步延長，模型訓練和驗證階段的模擬效果都穩(wěn)步提高. 當輸入時間窗為7 d時，模型各項性能指標已達到很好的效果，繼續(xù)增大輸入時間窗，模型預測效果沒有明顯提高. 因此，我們選擇利用“五河”六站及長江干流漢口站這七站前7 d內的流量來預測鄱陽湖湖口、星子、都昌、吳城和康山站未來1 d的水位. 其中，湖口站訓練和驗證階段的RMSE分別為0.62和0.55 m，NSE分別為0.97和0.98，R2均為0.98.

2.2.2 隱藏神經(jīng)元數(shù)目對模型模擬效果的影響 LSTM網(wǎng)絡層中隱藏神經(jīng)元數(shù)目是影響模型預測結果準確率的重要參數(shù)之一. 若數(shù)量太少，網(wǎng)絡不能很好地學習，需要訓練的次數(shù)也多，訓練精度也不高；若數(shù)量太多，訓練時間又較長，甚至可能導致模型不收斂[35-36]. 因此，我們進一步計算了隱藏神經(jīng)元數(shù)分別為1、5、10、20、30、40、50、100、200、500共10種模式下，根據(jù)前7 d“五河”六站及漢口站流量組合預測鄱陽湖5站未來 1 d 的水位空間分布，相應的訓練和驗證階段模型擬合評價指標的計算結果如圖5所示.

圖5 隱藏神經(jīng)元數(shù)目對模型訓練和驗證過程的模擬結果的影響Fig.5 The effect of the number of hidden nodes on the simulation results of model training and verification

計算結果表明，一定范圍內隨著隱藏神經(jīng)元數(shù)目的增加，模型在訓練和驗證階段的水位預測效果均穩(wěn)步提高，但當隱藏神經(jīng)元數(shù)目增加到一定值之后，繼續(xù)增加其數(shù)量模型預測效果變化不大. 隱藏神經(jīng)元數(shù)為1個時，各站水位預測效果最差. 湖口站水位訓練和驗證階段的RMSE均高達1.0 m以上，分別為1.02 和1.13 m，NSE和R2均相對較低. 隱藏神經(jīng)元數(shù)目為50個時，各站水位預測效果的評估指標均達到平穩(wěn). 其中，作為5個水位站中預測效果最差的康山站，RMSE僅為0.45 m，NSE為0.95，R2為0.96. 考慮到模型預測精度以及計算時間，我們設置了50個隱藏神經(jīng)元節(jié)點. 模型訓練階段的RMSE為0.61 m，NSE為0.97，R2為0.98；驗證階段的RMSE為0.53 m，NSE為0.98，R2也為0.98.

2.2.3 初始學習率大小對模型模擬效果的影響 學習率表示每次迭代后權重的更新量，學習率太小，模型更新速度慢；學習率過大，模型可能錯過最優(yōu)解. 為了找到模型的全局最優(yōu)解，我們借助基于梯度下降的具有相當魯棒性的ADAM自適應學習率優(yōu)化算法作為優(yōu)化器[37]，算法描述為：

mt=β1·mt-1+(1-β1)gt

(13)

(14)

(15)

(16)

(17)

圖6 初始學習率大小對模型訓練和驗證過程的模擬結果的影響Fig.6 The effect of the initial learning rate on the simulation results of model training and verification

計算了4種初始學習率(分別為0.1、0.01、0.001和0.0001)下模型訓練和驗證階段損失函數(shù)的大小(圖6). 我們發(fā)現(xiàn)，初始學習率為0.1時，損失函數(shù)曲線隨著迭代次數(shù)的增加發(fā)生不同幅度的震蕩，此時學習率選擇過大；當初始學習率為0.01時，訓練階段的損失函數(shù)具有較好的學習過程，但驗證階段的損失函數(shù)曲線迅速下降，模型學得太快，基于對模型精度和速度的綜合考慮，認為初始學習率過大；當初始學習率為 0.0001 時，則經(jīng)過較長時間模型才得以收斂. 因此，我們選取訓練和驗證階段均有很好的學習過程曲線的 0.001 作為最適初始學習率，且模型快速地收斂于0.03. 此時，五站中水位預測效果最好的湖口站，訓練和驗證階段的RMSE分別為0.58和0.50 m，NSE均為0.98，R2也均為0.98. 預測效果相對最差的康山站，驗證階段的RMSE為0.42 m，NSE和R2分別可達0.95和0.96. 康山站的RMSE值反而比湖口站的低，是因為各站年內水位變幅存在差異，康山站多年年內變幅均值為6.03 m，而湖口站可達12.00 m.

3 水位預測

3.1 預測結果分析

綜上所述，最終我們建立了含50個隱藏神經(jīng)元的單層LSTM神經(jīng)網(wǎng)絡，采用基于梯度下降的自適應學習率的ADAM優(yōu)化算法，初始學習率為0.001，利用“五河”六站和長江干流共7個流量站前7 d內的日均流量，組合預測鄱陽湖從北到南5個水位站未來1 d的日均水位. 將1956-1980年水文序列數(shù)據(jù)作為訓練集，1981-2000年數(shù)據(jù)對模型進行驗證. 構建的LSTM模型對湖口、星子、都昌、康山和吳城站的水位預測精度較高，模型擬合的各項評價指標如表1所示，訓練階段和驗證階段RMSE都很小，范圍分別為0.48～0.60 m和0.41～0.50 m. 同時，訓練和驗證階段的NSE和R2都很大，幾乎接近于1. 其中，訓練階段和驗證階段的NSE范圍分別為0.92～0.98和0.96～0.98；兩個階段的R2范圍則分別為0.93～0.98和0.96～0.98. 湖口站和星子站1981-2000年水位預測的驗證結果如圖7所示，以兩站2000年的水位過程線為例，湖口站年均水位實測值和預測值分別為12.98和13.26 m. 最高日均水位實測值為18.10 m，預測值較之大0.62 m，但與年均水位相比，相差較小，具有較強的可靠性. 日均最低水位的實測值與預測值分別為7.73和7.61 m，較年均和年最高水位值的預測結果更精確. 構建的LSTM模型對湖泊最下游湖口站的水位預測精度最高，從最上游康山站至湖口站，模型對五站水位預測精度有逐步增強的趨勢，這可能和湖區(qū)各站受長江的頂托關系逐漸趨弱有關.

表1 最佳模型評價指標計算結果

圖7 1981-2000年湖口站(a)和星子站(b)水位預測結果Fig.7 Water level prediction results of Hukou station (a) and Xingzi station (b) from 1981 to 2000

3.2 模型能力比較

為了充分驗證本文構建的LSTM模型的水位模擬精度，突出模型的優(yōu)勢以及后續(xù)模型應用，我們也用一般的循環(huán)神經(jīng)網(wǎng)絡BP神經(jīng)網(wǎng)絡對鄱陽湖水位進行模擬，并與LSTM模型的模擬精度進行對比. BP神經(jīng)網(wǎng)絡同樣是將1956-1980年水文序列數(shù)據(jù)作為訓練集，1981-2000年數(shù)據(jù)對模型進行驗證，利用“五河”六站和長江干流共7個流量站前7 d內的日均流量，組合預測鄱陽湖從上游到下游5個水位站未來1 d的日均水位. 比較了兩種方法的各項評價指標(表2).

表2 兩種神經(jīng)網(wǎng)絡方法驗證階段評價指標計算結果對比

圖8 三種訓練數(shù)據(jù)集下模型驗證結果比較Fig.8 Comparison of model verification results under three training data sets

對五站水位的預測，LSTM模型得到的最佳RMSE值均低于0.50 m，而BPNN模型得到的最佳RMSE值除最上游康山站為0.77 m，其他四站的RMSE值均高于0.80 m；LSTM模型得到的NSE，預測精度最低的康山站水位也達到0.96，而BPNN得到的NSE，預測精度最高的湖口站水位也僅有0.95. 不難看出，我們構建的LSTM模型模擬精度明顯強于BPNN模型. 除了模擬精度外，計算速度也是衡量模型性能的一個重要指標. 計算了兩種模型的計算速度，結果表明，LSTM模型的計算速度(21.66 s)明顯快于BPNN的計算速度(138.82 s). 綜上，構建的LSTM模型在精度和速度上較BPNN模型均具有明顯的優(yōu)勢.

3.3 訓練數(shù)據(jù)集的影響

實際情況中，應盡可能選取足夠多的數(shù)據(jù)來建模，但由于各種不可抗因素，往往無法獲取完整的水文時間序列數(shù)據(jù). 為此，我們需要考慮選用不同的訓練數(shù)據(jù)集，分析模型的效果. 我們基于水文時間序列數(shù)據(jù)自身的特點，利用上述參數(shù)設計好的LSTM模型，設計3種方案分別對五站水位進行預測，分別為：①訓練集數(shù)據(jù)為長時間序列(1956-1980年)；②訓練集數(shù)據(jù)較短時間序列(1956-1960年)；③訓練集數(shù)據(jù)時間序列為5個典型年：1954年、1973年、1974年、1977年和1978年. 根據(jù)《水文情報預報規(guī)范》(GB/T 22482-2008)中的距平百分率劃分徑流豐平枯的標準，1954年和1973年是典型的豐水年，1977年是平水年，1974年和1978年為典型的枯水年. 利用1981-2000年的數(shù)據(jù)進行驗證，比較3種訓練集數(shù)據(jù)模式下模型擬合的效果.

圖9 三種訓練集訓練模式下都昌站1981-2000年水位預測結果Fig.9 Water level prediction results of Duchang station from 1981 to 2000 under three training sets

五站水位預測效果的各項評價指標計算結果如圖8所示，同種預測方案中，湖口站水位預測效果均最好，康山站均最差. 方案①中，湖口站驗證階段RMSE、NSE和R2分別為0.51 m、0.98和0.98；康山站相應的指標計算結果分別為0.43 m、0.95和0.95. 方案②中，五站水位預測的效果明顯降低. 湖口站RMSE、NSE和R2分別為0.93 m、0.94和0.97；康山站相應的指標計算結果分別為0.88 m、0.80和0.91. 方案①和方案②中訓練集的樣本量分別為9126和1827. 而將樣本量為1820個的5個典型年的數(shù)據(jù)作為訓練數(shù)據(jù)進行模型訓練時，五站水位的預測效果與方案①相比較差，但明顯優(yōu)于方案②的預測效果. 湖口站RMSE、NSE和R2分別為0.76 m、0.96和0.96；康山站相應的指標計算結果分別為0.56 m、0.92和0.93. 所以，我們在利用LSTM神經(jīng)網(wǎng)絡預測湖泊水位時，應盡可能選取足夠長時間序列的數(shù)據(jù). 若因為不可抗因素無法獲取完整的數(shù)據(jù)序列，選取涵蓋各種代表性數(shù)據(jù)的典型年數(shù)據(jù)進行訓練，也可以獲得較好的模型預測效果.

都昌站代表鄱陽湖主湖區(qū)水位，3種方案下其水位預測的結果如圖9所示(其他水位站預測結果與其相似，故圖省略)，當1956-1980年數(shù)據(jù)作為訓練時間序列時，五站多年水位的預測值與真實值之間具有非常高的對應關系. 1956-1960年是鄱陽湖典型的枯水期，因而模型通過訓練能夠較好地反映低水位特征. 5個典型年涵蓋了幾個重要的豐、枯水年信息，總體預測效果介于前兩者之間.