紀國良，周 曼，劉 濤，胡騰騰，丁 勇

(中國長江三峽集團有限公司， 湖北 宜昌 443100)

1 研究背景

大型水庫重要站點的水位預測是水庫防洪中的重要問題，目前主要采用水動力學方法計算。水動力學方法基于嚴格的圣維南方程組原理[1]，對輸入邊界條件的準確性有較高要求，當邊界條件清晰時，其計算精度很高；否則，容易造成較明顯的計算誤差。在大型水庫的防洪調(diào)度中，干流和較大支流的入庫流量通常有較為可靠的預報，而水庫區(qū)間來流則是根據(jù)流域降雨情況進行估算，其空間的分布通常難以準確預測，因此區(qū)間流量成為最難準確界定的邊界條件，使得水動力學方法在計算某些站點的水位精度時難以滿足實時應(yīng)用需求。本文針對這一問題開展研究， 使用基于數(shù)據(jù)驅(qū)動的循環(huán)神經(jīng)網(wǎng)絡(luò)算法，在區(qū)間來流空間分布未知的情況下，提高水庫重要站點的水位計算精度。

目前已有許多工作致力于提高水動力學模型的計算精度。長江科學院黃仁勇等[2]在計算三峽水庫洪水傳播過程時，一方面增加了水庫斷面測量點數(shù)量，降低了斷面概化誤差；另一方面不斷提升入庫流量的預報水平并將更多支流加入到模型中。在糙率率定方面，中山大學楊世孝等[3]提出了一種使用最優(yōu)化方法反求糙率的方式；陳素紅等[4]建立了基于多親遺傳算法的河道糙率率定模型；陳一帆等[5]以糙率和水力狀態(tài)量作為河網(wǎng)非線性動態(tài)系統(tǒng)狀態(tài)變量，采用擴展卡爾曼濾波構(gòu)建了結(jié)合糙率動態(tài)校正的河網(wǎng)水情數(shù)據(jù)同化模型。這些研究已經(jīng)取得了很大的進展，但在區(qū)間來流空間分布難以確定的情況下，其計算精度仍有進一步改進的空間。

本文采用循環(huán)神經(jīng)網(wǎng)絡(luò)模型[6]，研究干、支流入庫流量、壩前水位等因素到重要站點水位的直接映射關(guān)系。在應(yīng)用部分，以三峽水庫長壽站水位計算為目標，使用2009—2019年的數(shù)據(jù)訓練、驗證和測試模型，應(yīng)用結(jié)果顯示在不同的水情分類下，循環(huán)神經(jīng)網(wǎng)絡(luò)的預測效果優(yōu)于傳統(tǒng)水動力學方法。

2 傳統(tǒng)水動力學方法

傳統(tǒng)水動力學方法主要是通過求解圣維南方程組[1]來獲取斷面的水位和流量，一維圣維南方程描述如下。水流連續(xù)方程見式(1)，水流運動方程見式(2)。

(1)

(2)

式中：Z表示水位；t為時間；B為河寬；Q為流量；q、v分別為匯入口處的支流流量、流速；A為過水斷面面積；g為重力加速度；K為斷面的流量模數(shù)；0≤x≤l，l表示河長。對于單河道計算，還需要如下初邊值條件:

①初始條件。在t=0時，河道每個斷面的水位和流量Z|t=0=Z(x)，Q|t=0=Q(x)。

②邊值條件。河道的上游入口和下游出口分別給出在任意時刻t的流量過程線Q=Q(t)和水位過程線Z=Z(t)，或者兩個邊值條件均為流量與水位之間的關(guān)系Q=Z(t)；如果是河網(wǎng)計算，還需要給出每個支流的入流過程線Q支=Q支(t)。在水庫計算中，一般是應(yīng)用上游流量過程和下游水位過程作為邊值條件。

作為比較，本文亦應(yīng)用傳統(tǒng)水動力學法計算重要站點水位，假設(shè)區(qū)間流量沿程均勻分布。

3 循環(huán)神經(jīng)網(wǎng)絡(luò)方法

3.1 單層循環(huán)神經(jīng)網(wǎng)絡(luò)

圖1展示了單層循環(huán)神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)，其中W和U是網(wǎng)絡(luò)參數(shù)，h是相應(yīng)節(jié)點的隱藏狀態(tài)，x是輸入向量，下標t表示時刻。它的計算順序遵從時間的先后順序，從左至右逐步合成特征向量的信息，完成特征的提取過程。其中每個矩形代表一個計算節(jié)點，常用的計算節(jié)點有傳統(tǒng)節(jié)點、LSTM(Long Short-Term Memory)[7]節(jié)點和GRU(Gated Recurrent Unit)[8]節(jié)點，具體計算方式如下。

圖1 單層循環(huán)神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)Fig.1 Structure of recurrent neural network of a single layer

傳統(tǒng)計算節(jié)點公式為

ht=σ(Wxt+Uht-1+b) 。

(3)

式中:b是偏置項參數(shù);σ(x)=1/(1+e-x)是sigmoid函數(shù)[9]。

LSTM計算節(jié)點公式為：

(4)

GRU計算節(jié)點公式為

(5)

與傳統(tǒng)計算節(jié)點相比，LSTM和GRU節(jié)點可以有效避免梯度爆炸或者梯度消失現(xiàn)象[7-8]，能夠處理更長序列的數(shù)據(jù)；LSTM與GRU相比，前者的計算量比較大，對數(shù)據(jù)量的要求較高，但是二者在許多任務(wù)上的效果近似?？紤]到本文研究問題的數(shù)據(jù)規(guī)模，選擇GRU節(jié)點進行實驗。

在預測任務(wù)中，將網(wǎng)絡(luò)中最后一個計算節(jié)點的隱藏向量hn看作模型的預測值，假設(shè)共有N個樣本，hn,i表示第i個樣本的預測值，yi表示第i個樣本的實測值，則定義目標函數(shù)為

(6)

式中：λ是系數(shù)，第一項是最小二乘項，使得預測值盡可能接近實測值；第二項是正則項，能有效防止模型過擬合。

3.2 深層循環(huán)神經(jīng)網(wǎng)絡(luò)

對于復雜型數(shù)據(jù)，單層循環(huán)神經(jīng)網(wǎng)絡(luò)有時候難以學習數(shù)據(jù)之間的映射關(guān)系，需要加深網(wǎng)絡(luò)層次。圖2是深層循環(huán)神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)圖，深層循環(huán)神經(jīng)網(wǎng)絡(luò)是多個單層循環(huán)神經(jīng)網(wǎng)絡(luò)的疊加，下一層的輸出是上一層的輸入，先計算第一層，計算完畢后，第一層所有節(jié)點的隱藏狀態(tài)作為第二層的輸入，然后再計算第二層，依此類推，直到最后一層計算完畢。最后一層最后一個節(jié)點的輸出是整個網(wǎng)絡(luò)的輸出，記為hm,n，表示第m層的第n個節(jié)點的輸出。與式(6)相同，目標函數(shù)定義為

圖2 深層循環(huán)神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)Fig.2 Structure of recurrent neural network of multilayer

(7)

4 模型的搭建和求解方法

本文在構(gòu)建模型和求解方面使用基于圖結(jié)構(gòu)的建模語言TensorFlow[10]實現(xiàn)。TensorFlow是基于python[11]開發(fā)的軟件包，是一個采用數(shù)據(jù)流圖架構(gòu)設(shè)計、用于數(shù)值計算的開源軟件庫，能夠自動地根據(jù)變量之間的依賴關(guān)系運用鏈式法則求導梯度，方便模型建立和求解。

5 三峽水庫長壽站水位計算

三峽水庫蓄水運用以來，根據(jù)以往的研究成果和實際運行管理經(jīng)驗，長壽區(qū)域是回水敏感區(qū)域，該區(qū)域長壽站是庫尾重要的水位控制性站點，是研究的重點，因此本文選擇長壽站作為研究對象。

5.1 樣本構(gòu)造

三峽水庫地域廣，涉及到的支流眾多，結(jié)合圣維南方程組的計算原理和實際經(jīng)驗，選取影響長壽站水位的主要因素為三峽入庫流量、寸灘流量、武隆流量、出庫流量和壩前水位。其中寸灘處于長壽的上游，是上游來水的代表站；“三峽入庫流量”指的是三峽水庫庫總的來流量，包含區(qū)間流量，這里不考慮區(qū)間來流的空間分布；武隆站是長壽下游烏江的控制站點，其來水對長壽站水位具有較顯著的頂托作用；出庫流量和壩前水位反映的是水庫出口的控制條件。

在主要因素選取之后，需要考慮影響當前長壽站水位的時段數(shù)?？紤]到洪水調(diào)度實際情況，水文預報是間隔6 h一個值，因此相鄰時段的間隔取6 h；由于寸灘流量對長壽站水位的影響最大，且水流從寸灘到達長壽站平均時間約6 h，因此長壽站水位預報計算的輸入取2個時段數(shù)據(jù)為宜，即當前和6 h前兩個時刻。表1和表2給出了一個樣本的示例，表1是模型的輸入，表2是模型的輸出，即若要預測2018年1月1日 10時長壽站的水位，則需要給出2018年1月1 日4時和2018年1月1日10時三峽入庫流量、寸灘流量、武隆流量、三峽出庫流量和鳳凰山水位(壩前水位)。在實際使用時，均使用水文預報數(shù)據(jù)，因此本方法的預見期同水文預報預見期。

表1 樣本輸入示例Table 1 Feature vectors of a sample

表2 樣本輸出示例Table 2 Outputs of a sample

根據(jù)第3節(jié)描述的數(shù)學形式，預測2018年1月1 日10時長壽站水位時模型的輸入為

同時，模型輸出為173.75 m?？紤]到流量和水位的單位不同，以及神經(jīng)網(wǎng)絡(luò)的輸入值不宜過大(否則其非線性項梯度接近0，模型無法優(yōu)化)，需要將向量進行歸一化處理，同時消除單位和數(shù)值大小的影響，流量和水位的歸一化處理如下：

歸一化值=(當前值-最小值)/(最大值-最小值)。

在收集的數(shù)據(jù)中，流量最大值為71 200 m3/s，最小值為995 m3/s；最高水位為177.09 m，最低水位為145 m，因此歸一化后得到模型的輸入向量為

對應(yīng)的模型輸出為0.898 4。經(jīng)過上述處理，水位、流量的數(shù)值均無量綱，且壓縮在區(qū)間[0,1]內(nèi)，在預測應(yīng)用時按照逆過程還原即可獲得水位值。

5.2 數(shù)據(jù)收集和水情分類

為了使數(shù)據(jù)樣本涵蓋不同種類的水雨情，本文收集2009年1月1日0時至2019年11月18 日23時的歷史運行數(shù)據(jù)，共95 376條小時數(shù)據(jù)，其中2009—2018年的數(shù)據(jù)用于訓練和驗證模型，2019年的數(shù)據(jù)用于測試模型。由于不同時期的壩前水位和入庫流量差別較大，因此需要根據(jù)不同的壩前水位和入庫流量進行分類，針對每個類別單獨訓練模型。三峽水庫2009—2019年入庫流量的取值范圍是[3 320, 71 200] m3/s，壩前水位的范圍是[145, 175]m，綜合考慮不同水情和模型對數(shù)據(jù)量的要求，共分為7類，其中樣本根據(jù)5.1節(jié)的描述生成，訓練集和驗證集是對樣本(由2009—2018年數(shù)據(jù)生成)按照90%和10%隨機劃分生成；測試集是2019年數(shù)據(jù)生成的樣本。生成樣本的數(shù)據(jù)分類情況如表3所示。

表3 數(shù)據(jù)分類情況Table 3 Data of different classes

5.3 模型超參數(shù)設(shè)置

在5.2節(jié)中，根據(jù)不同的水情將數(shù)據(jù)分成了7個類別，需要對每個類別設(shè)置網(wǎng)絡(luò)的層數(shù)、每層隱藏節(jié)點個數(shù)、優(yōu)化算法學習率α、正則項系數(shù)λ、批量訓練樣本數(shù)等超參數(shù)。超參數(shù)的選擇首先根據(jù)經(jīng)驗確定范圍，然后在驗證集上測試其效果，選擇驗證集上效果最好的超參數(shù)。表4給出了最佳的超參數(shù)設(shè)置。

表4 超參數(shù)設(shè)置Table 4 Setting of hyper-parameter

5.4 試驗結(jié)果和誤差分析

5.4.1 試驗結(jié)果的整體比較

圖3(a)是2009—2019年所有樣本(包含所有類別的訓練集、驗證集和測試集 )的計算水位和實測水位對比圖，具體時間是從2009年1月1日 6時至2019年11月18日 17時。圖3(b)是對應(yīng)的誤差分布，表5是7個類別的訓練集、驗證集和測試集的誤差上下限。

表5 不同類別的訓練集、驗證集和測試集誤差范圍Table 5 Error ranges of training, validation,and test data sets

圖3 2009—2019年長壽站水位計算值與實測值對比和計算值與實測值誤差分布Fig.3 Comparison between calculated and measured water levels of Changshou station from 2009 to 2019

由圖3可知，預測水位曲線與實測水位曲線距離很近，說明模型對數(shù)據(jù)的擬合和預測效果較好。

在水位較高時，預測線和實測線幾乎重合，計算精度較高，滿足實際應(yīng)用的需求。從表5中可以看出，較大誤差主要發(fā)生在Ⅰ、Ⅱ和Ⅶ類中，且Ⅰ和Ⅶ類均出現(xiàn)了超過1 m以上的誤差。Ⅰ、Ⅱ類誤差較大的原因是壩前水位較低，河槽較窄，流量略微增大后水位就會有顯著變化，因而計算難度較大；Ⅶ類是三峽入庫流量在25 000 m3/s以上，壩前水位在170 m以下的條件，此時入庫流量較大，水位變化幅度大，預測難度相應(yīng)增大，因此最大誤差有-1.27 m。以上2種情況都發(fā)生在長壽站水位較低時，因此不影響實際中的應(yīng)用。

5.4.2 與傳統(tǒng)水動力學法比較

2010年、2012年、2014年、2017年是長江流域典型洪水年，其中2010年、2012年、2014年洪水發(fā)生在汛期，2017年發(fā)生在蓄水期。

圖4分別是2010年、2012年、2014年、2017年長壽站實測水位、水動力學法計算水位和循環(huán)神經(jīng)網(wǎng)絡(luò)計算水位的對比。從圖4可以看出，較大誤差多發(fā)生在水位低的情況下；在高水位時2種計算方法的誤差都顯著縮小；比較2種計算方法的誤差，循環(huán)神經(jīng)網(wǎng)絡(luò)的誤差在絕大多數(shù)情況下小于水動力學方法，說明了本文方法的有效性。

圖4 不同年份水動力學法與循環(huán)神經(jīng)網(wǎng)絡(luò)法對比Fig.4 Comparison of the prediction result between hydrodynamic method and recurrent neural network under different water regimens in different years

5.4.3 循環(huán)神經(jīng)網(wǎng)絡(luò)方法有效性的原因分析

循環(huán)神經(jīng)網(wǎng)絡(luò)的輸入實際上采用的是水動力學方法的邊界條件，它的有效性可以歸納為以下兩點：

(1)循環(huán)神經(jīng)網(wǎng)絡(luò)能有效避開區(qū)間流量空間分布未知的問題。在水動力學方法中，圣維南方程組經(jīng)過離散和線性化以后，求解的過程是依靠壩前水位逐步向上游計算各斷面的水位，上游斷面的水位對下游斷面水位具有很強的依賴關(guān)系，當區(qū)間來水空間分布未知時，水庫中間的部分斷面水位計算就會產(chǎn)生誤差，由于遞推關(guān)系，這些誤差會向上游傳播，最終導致長壽站水位計算誤差較大。循環(huán)神經(jīng)網(wǎng)絡(luò)方法只使用了區(qū)間來流的總流量(包含在入庫流量中)，忽略了其空間分布，反而模擬效果更好，主要原因是在大多數(shù)情況下，由于距離較遠，區(qū)間來流對長壽站水位的頂托作用較小，盡管區(qū)間來流會影響庫區(qū)中間的水位，但是模型的輸入只依賴壩前水位，并不依賴中間水位，所以能有效避開區(qū)間來流的空間分布問題。誠然，當區(qū)間流量空間分布未知時，上述2種方法都只適用于區(qū)間來流占比較小的情況，占比較大時都會產(chǎn)生較大的計算誤差。

(2)循環(huán)神經(jīng)網(wǎng)絡(luò)沒有使用地形資料和糙率參數(shù)，避免了數(shù)據(jù)測量誤差和率定誤差，因此能提高預測精度。

6 結(jié) 論

(1)在大型水庫重要站點水位計算中，循環(huán)神經(jīng)網(wǎng)絡(luò)方法是直接挖掘輸入數(shù)據(jù)到輸出數(shù)據(jù)的關(guān)系，降低了邊界條件準確性的要求，只需要使用庫周入庫流量，有效避開了區(qū)間來流空間分布未知的問題；同時也不需要使用地形資料和糙率參數(shù)等，沒有測量和率定誤差，因此它的預測精度較高，并且精度會隨著運行數(shù)據(jù)的積累而不斷提高。

(2)循環(huán)神經(jīng)網(wǎng)絡(luò)可應(yīng)用于其他問題的研究。由于水文領(lǐng)域的數(shù)據(jù)均具有時間屬性，非常適合此種網(wǎng)絡(luò)結(jié)構(gòu)，因此它可在流量、水文、泥沙等方面進行推廣。

(3)考慮到三峽庫區(qū)一直處于淤積過程中，局部集中淤積會對庫水位產(chǎn)生較大影響，神經(jīng)網(wǎng)絡(luò)模型訓練使用的數(shù)據(jù)包含了淤積前和淤積后，數(shù)據(jù)前后的一致性略有差別，本文下一步工作是研究泥沙淤積對水位計算的影響。