紀(jì)國(guó)良，周 曼，劉 濤，胡騰騰，丁 勇

(中國(guó)長(zhǎng)江三峽集團(tuán)有限公司， 湖北 宜昌 443100)

1 研究背景

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

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

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

2 傳統(tǒng)水動(dòng)力學(xué)方法

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

(1)

(2)

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

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

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

作為比較，本文亦應(yīng)用傳統(tǒng)水動(dòng)力學(xué)法計(jì)算重要站點(diǎn)水位，假設(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é)點(diǎn)的隱藏狀態(tài)，x是輸入向量，下標(biāo)t表示時(shí)刻。它的計(jì)算順序遵從時(shí)間的先后順序，從左至右逐步合成特征向量的信息，完成特征的提取過(guò)程。其中每個(gè)矩形代表一個(gè)計(jì)算節(jié)點(diǎn)，常用的計(jì)算節(jié)點(diǎn)有傳統(tǒng)節(jié)點(diǎn)、LSTM(Long Short-Term Memory)[7]節(jié)點(diǎn)和GRU(Gated Recurrent Unit)[8]節(jié)點(diǎn)，具體計(jì)算方式如下。

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

傳統(tǒng)計(jì)算節(jié)點(diǎn)公式為

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

(3)

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

LSTM計(jì)算節(jié)點(diǎn)公式為：

(4)

GRU計(jì)算節(jié)點(diǎn)公式為

(5)

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

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

(6)

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

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

對(duì)于復(fù)雜型數(shù)據(jù)，單層循環(huán)神經(jīng)網(wǎng)絡(luò)有時(shí)候難以學(xué)習(xí)數(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ò)是多個(gè)單層循環(huán)神經(jīng)網(wǎng)絡(luò)的疊加，下一層的輸出是上一層的輸入，先計(jì)算第一層，計(jì)算完畢后，第一層所有節(jié)點(diǎn)的隱藏狀態(tài)作為第二層的輸入，然后再計(jì)算第二層，依此類推，直到最后一層計(jì)算完畢。最后一層最后一個(gè)節(jié)點(diǎn)的輸出是整個(gè)網(wǎng)絡(luò)的輸出，記為hm,n，表示第m層的第n個(gè)節(jié)點(diǎn)的輸出。與式(6)相同，目標(biāo)函數(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)的建模語(yǔ)言TensorFlow[10]實(shí)現(xiàn)。TensorFlow是基于python[11]開(kāi)發(fā)的軟件包，是一個(gè)采用數(shù)據(jù)流圖架構(gòu)設(shè)計(jì)、用于數(shù)值計(jì)算的開(kāi)源軟件庫(kù)，能夠自動(dòng)地根據(jù)變量之間的依賴關(guān)系運(yùn)用鏈?zhǔn)椒▌t求導(dǎo)梯度，方便模型建立和求解。

5 三峽水庫(kù)長(zhǎng)壽站水位計(jì)算

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

5.1 樣本構(gòu)造

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

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

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

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

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

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

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

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

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

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

為了使數(shù)據(jù)樣本涵蓋不同種類的水雨情，本文收集2009年1月1日0時(shí)至2019年11月18 日23時(shí)的歷史運(yùn)行數(shù)據(jù)，共95 376條小時(shí)數(shù)據(jù)，其中2009—2018年的數(shù)據(jù)用于訓(xùn)練和驗(yàn)證模型，2019年的數(shù)據(jù)用于測(cè)試模型。由于不同時(shí)期的壩前水位和入庫(kù)流量差別較大，因此需要根據(jù)不同的壩前水位和入庫(kù)流量進(jìn)行分類，針對(duì)每個(gè)類別單獨(dú)訓(xùn)練模型。三峽水庫(kù)2009—2019年入庫(kù)流量的取值范圍是[3 320, 71 200] m3/s，壩前水位的范圍是[145, 175]m，綜合考慮不同水情和模型對(duì)數(shù)據(jù)量的要求，共分為7類，其中樣本根據(jù)5.1節(jié)的描述生成，訓(xùn)練集和驗(yàn)證集是對(duì)樣本(由2009—2018年數(shù)據(jù)生成)按照90%和10%隨機(jī)劃分生成；測(cè)試集是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個(gè)類別，需要對(duì)每個(gè)類別設(shè)置網(wǎng)絡(luò)的層數(shù)、每層隱藏節(jié)點(diǎn)個(gè)數(shù)、優(yōu)化算法學(xué)習(xí)率α、正則項(xiàng)系數(shù)λ、批量訓(xùn)練樣本數(shù)等超參數(shù)。超參數(shù)的選擇首先根據(jù)經(jīng)驗(yàn)確定范圍，然后在驗(yàn)證集上測(cè)試其效果，選擇驗(yàn)證集上效果最好的超參數(shù)。表4給出了最佳的超參數(shù)設(shè)置。

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

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

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

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

表5 不同類別的訓(xùn)練集、驗(yàn)證集和測(cè)試集誤差范圍Table 5 Error ranges of training, validation,and test data sets

圖3 2009—2019年長(zhǎng)壽站水位計(jì)算值與實(shí)測(cè)值對(duì)比和計(jì)算值與實(shí)測(cè)值誤差分布Fig.3 Comparison between calculated and measured water levels of Changshou station from 2009 to 2019

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

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

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

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

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

圖4 不同年份水動(dòng)力學(xué)法與循環(huán)神經(jīng)網(wǎng)絡(luò)法對(duì)比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ò)的輸入實(shí)際上采用的是水動(dòng)力學(xué)方法的邊界條件，它的有效性可以歸納為以下兩點(diǎn)：

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

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

6 結(jié) 論

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

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

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