圖神經(jīng)網(wǎng)絡(luò)驅(qū)動的流域洪水預(yù)報技術(shù)

馬森標

（福建中銳網(wǎng)絡(luò)股份有限公司研發(fā)中心，福建 福州，350108）

0 引言

流域洪水預(yù)報是防汛搶險和防洪系統(tǒng)調(diào)度運用的決策依據(jù)，傳統(tǒng)的洪水預(yù)報采用的機理模型[1]屬于半經(jīng)驗半物理過程的預(yù)報模型，有一定的物理基礎(chǔ)，但模型參數(shù)的率定十分復雜和困難[2]，基于數(shù)據(jù)驅(qū)動的模型可根據(jù)海量的歷史水文數(shù)據(jù)挖掘其潛在變化規(guī)律，具有比較高的性能和價值。

由于數(shù)據(jù)驅(qū)動模型的優(yōu)點，目前國內(nèi)已有很多專家學者將數(shù)據(jù)驅(qū)動的人工智能模型引入到洪水預(yù)報中。劉冬英等[3]提出了改進的BP 網(wǎng)絡(luò)模型在洪水預(yù)報中的應(yīng)用，取得較高精度的預(yù)報結(jié)果；金保明等[4]利用RBF 神經(jīng)網(wǎng)絡(luò)預(yù)報山區(qū)流域洪水的方法，取得了不錯的預(yù)報精度效果；在閩江十里庵流量預(yù)報中應(yīng)用BP 神經(jīng)網(wǎng)絡(luò)，預(yù)報精度符合要求[5]；李曉麗等[6]提出了不確定支持向量機在洪水預(yù)報模型中的應(yīng)用，獲得了較高的洪水預(yù)報精度；馬森標等[7]利用長短時記憶神經(jīng)網(wǎng)絡(luò)（LSTM）附加Attention 機制進行水庫水位預(yù)報，取得了較好的性能。雖然以上算法在一定程度上提高了洪水的預(yù)報精確度，但是對于輸入的數(shù)據(jù)特征，尤其是關(guān)鍵性的降雨、水位數(shù)據(jù)特征挖掘明顯不足，關(guān)注了時序特征，卻忽略了其本身具有的空間分布特征。同時，降雨量對水位的影響受到流域的地形、土壤、植被等自然因素的作用，表現(xiàn)出了顯著的滯后性、持續(xù)性等特點。同樣，受地理位置和流域地形的作用，泄洪量對水位的影響，及上游水位對下游水位的影響同樣存在滯后性、持續(xù)性等特點。

圖卷積神經(jīng)網(wǎng)絡(luò)（GCN）[8]是一種在圖結(jié)構(gòu)數(shù)據(jù)上運行的基于深度學習的擴展方法，可以挖掘符合非歐幾里得空間域的流域數(shù)據(jù)結(jié)構(gòu)，它能夠從底層數(shù)據(jù)中提取出網(wǎng)絡(luò)節(jié)點復雜的關(guān)聯(lián)性[9]。

基于上述出發(fā)點，本研究提出了一個圖指導的時空關(guān)聯(lián)預(yù)報模型（GSCPM，graph-guided spatiotemporal correlation prediction model），用以預(yù)報流域水位。具體而言，該模型首先通過多個長短期記憶網(wǎng)絡(luò)（LSTM）編碼每個監(jiān)測點歷史屬性的時間關(guān)聯(lián)特征，隨后利用挖掘監(jiān)測點間的地理空間依賴，最后通過一個全連接神經(jīng)網(wǎng)絡(luò)進行回歸預(yù)報，可解決流域洪水預(yù)報中的時空關(guān)系建模和滯后影響問題。

1 模型概述

在本節(jié)中，我們將詳細介紹所提出的模型的實現(xiàn)細節(jié)，該模型的概述如圖1 所示。

圖1 GSCPM 模型

1.1 問題定義

流域水位預(yù)報問題的本質(zhì)為多變量時間序列預(yù)報問題。設(shè)流域中有K 個監(jiān)測點，第i 個監(jiān)測點的可觀測數(shù)據(jù)被定義為，其中t 為當前時間，d為過去時間天數(shù)，x∈Rn為該觀測點包含的n 個變量，包括水位、降雨量和水庫泄洪量。給定所有監(jiān)測點的可觀測數(shù)據(jù)，我們的目標是預(yù)報未來1 天的監(jiān)測點水位

1.2 時間依賴編碼

在獲取K 個監(jiān)測點的歷史數(shù)據(jù)后，我們先通過LSTM 編碼其時間依賴關(guān)系。由于每個輸入點的特征維度可能不盡相同，并且只使用一個LSTM 可能無法區(qū)分不同監(jiān)測點數(shù)據(jù)的特異性，我們?yōu)槊總€監(jiān)測點設(shè)置了一個LSTM 網(wǎng)絡(luò)，從而得到包含時間依賴的監(jiān)測點嵌入：

其中，每個LSTM 的具體執(zhí)行過程如下：

式中：Wf、Wi、WO、WC為可訓練權(quán)重矩陣；bf、bi、bO、bC為可訓練偏置矩陣；σ（·）和tanh（·）為非線性激活函數(shù)；ht-1為前一時刻的LSTM 單元的輸出；[ ]∵ 為向量連接操作。

1.3 空間依賴編碼

為編碼獲得時空融合特征，我們在時間依賴編碼后，使用GCN 對所得嵌入進行基于地理空間關(guān)系的特征聚合。

一個地理信息圖可由G=（V，E，C）描述。其中，V 表示圖的節(jié)點集合，其中節(jié)點vi∈V 表示第i 個監(jiān)測點；E表示邊的集合，其中邊eij∈E 表示節(jié)點vi和vj中存在關(guān)系；C 為節(jié)點的特征集合，其中ci∈C 為第i 個節(jié)點的屬性值，我們將其設(shè)置為第i 個監(jiān)測點的嵌入表示，即ci=ei。為使地理信息圖能表示實際監(jiān)測點間的地理關(guān)系，我們提出了一個簡單有效的地理信息圖構(gòu)造方法。具體而言，遍歷每個監(jiān)測點，將其和流域河道中相鄰的監(jiān)測點都進行連接，如果該條邊已經(jīng)存在，則遍歷下一個監(jiān)測點。圖2 展示了構(gòu)圖的一個示例。

圖2 拓撲結(jié)構(gòu)圖生成示例

經(jīng)過上述構(gòu)圖操作，流域監(jiān)測點的地理信息被融入圖的拓撲結(jié)構(gòu)之中。隨后，使用一個GCN 聚合節(jié)點屬性，從而將圖的拓撲信息融入節(jié)點嵌入中。一個圖卷積層的執(zhí)行過程如下：

式中：D和A分別為圖G 的度矩陣和鄰接矩陣；W和b 分別是可訓練的權(quán)重和偏置矩陣；l 和L 分別是GCN 執(zhí)行的當前層數(shù)和最大層數(shù)。Zl為第l 層輸出節(jié)點嵌入矩陣，同時作為下一層的輸入。第一層GCN 的輸入Z0為節(jié)點特征集C 中每個元素拼接而成的矩陣，GCN 的最終輸出ZL被按行拆分K 個向量，代表每個監(jiān)測點包含時空信息的嵌入表達。

1.4 回歸預(yù)測

為進行最終的水位預(yù)報，K 個監(jiān)測點的嵌入向量被按列拼接為新的向量Z。隨后，采用一個全連接網(wǎng)絡(luò)進行回歸預(yù)測：

式中：WM和bM為可訓練的權(quán)重和偏置矩陣；y 為最終預(yù)報的水位值集合，包含未來1 天的K 個監(jiān)測點的水位值。最后，采用均方誤差作為損失函數(shù)，用以反向傳播訓練：

式中：Y 為實際的水位值集合。

2 溪源流域水位預(yù)報研究

2.1 數(shù)據(jù)基本情況

本文選取福州溪源流域為研究對象，該流域共包含溪源水庫及下游的侯官和融僑3 個監(jiān)測點，其中溪源水庫監(jiān)測站點有泄洪量、降雨量數(shù)據(jù)，侯官和融僑站點有水位數(shù)據(jù)。通過章節(jié)1.3 所提出的構(gòu)圖方法，對3 個監(jiān)測點進行構(gòu)圖，所得流域的時空示意圖[9]如圖3 所示。

圖3 局部流域時空示意圖

2.2 特征工程

通過機理和數(shù)學方法創(chuàng)造的特征，可以考慮到降雨量、泄洪量對水位的影響及上游水位對下游水位的影響都受到流域的地形、土壤、植被諸多自然因素的作用，構(gòu)造的特征將融合滯后性、持續(xù)性特點，根據(jù)這些特征訓練的模型具有更好的性能，這也是數(shù)據(jù)驅(qū)動模型優(yōu)于機理模型的重要原因之一。

2.2.1 滯后特征構(gòu)造

目前，在機器學習領(lǐng)域，降雨量、泄洪量和上游水位對流域水位的滯后效應(yīng)，很難通過計算得到，通過構(gòu)造相應(yīng)的滯后特征，來模擬這些滯后影響，簡化計算方法，提高建模效率。

一階和二階滯后特征的構(gòu)造過程如圖4 所示。假設(shè)原序列為｛xt-d…，xt-1，xt｝，其m 階段滯后則為其向左平移m 個時間點的序列，表示為｛xt-d-m…，xt-1-m，xt-m｝。在此基礎(chǔ)上，設(shè)Rt(i)為t 時刻滯后（i-1）個時間單位的降雨量序列；Ft(i)為t 時刻滯后（i-1）個時間單位的泄洪量序列；Wut(i)為t 時刻滯后（i-1）個時間單位的流域上游水位序列；Wdt(i)為t 時刻滯后（i-1）個時間單位的流域下游水位序列，其中i≥2。

圖4 滯后因素構(gòu)造過程示例

2.2.2 降雨量滯后特征與水位特征的相關(guān)性分析

流域監(jiān)測點在某一時刻的水位跟之前的降雨量密切相關(guān)，相關(guān)程度可以通過相關(guān)性的強弱來間接表示（見圖5 圖6）。

圖5 雨量曲線圖

圖6 m 階降雨量滯后特征與上、下游水位相關(guān)系數(shù)曲線圖

從表1 的相關(guān)系數(shù)可以看出，降雨量對水位的影響具有滯后效應(yīng)，降雨量對水位的影響是逐漸加強的，隨后逐漸減退，與蓄滿產(chǎn)流的自然規(guī)律具有一致性。

表1 m 階降雨量滯后特征與上、下游水位相關(guān)系數(shù)表

2.2.3 泄洪量滯后特征與水位特征的相關(guān)性分析

泄洪量對水位的影響巨大，在大壩某一時刻段的泄洪量，將持續(xù)影響流域監(jiān)測點的水位，這個影響可以通過泄洪量和水位的相關(guān)性進行衡量，影響程度可以通過相關(guān)性的強弱來間接表示（見圖7 圖8）。

圖7 泄洪量曲線圖

圖8 m 階泄洪量滯后特征與上、下游水位相關(guān)系數(shù)曲線圖

從表2 的相關(guān)系數(shù)可以看出，泄洪量對上、下游水位的影響具有持續(xù)性，泄洪量對上、下游水位的持續(xù)影響是漸漸變?nèi)醯模季哂袕娤嚓P(guān)性，符合水庫泄洪河道水位陡漲的自然客觀規(guī)律。

表2 m 階泄洪量滯后特征與上、下游水位相關(guān)系數(shù)表

2.2.4 上游水位滯后特征與下游水位的相關(guān)性分析

在某一時刻流域下游監(jiān)測點的水位跟該時刻前的上游水位密切相關(guān)，相關(guān)程度可以通過相關(guān)性的強弱來間接表示（見圖9 圖10）。

圖9 上、下游水位曲線圖

圖10 m 階上游水位滯后特征與下游水位相關(guān)系數(shù)曲線圖

從表3 的相關(guān)系數(shù)可以看出，下游水位跟上游水位具有較強的相關(guān)性，隨著時間的推移，相關(guān)性逐漸衰減，符合某一時刻的洪水波在傳播過程中發(fā)生衰減演進的客觀規(guī)律。

表3 m 階上游水位滯后特征與下游水位相關(guān)系數(shù)表

2.3 數(shù)據(jù)集構(gòu)建及數(shù)據(jù)預(yù)處理

基于章節(jié)2.2 中所展示的結(jié)論，我們添加相關(guān)度最高的滯后特征作為模型的輸入。添加特征后，在溪源水庫監(jiān)測點包含泄洪量、泄洪量的1 階滯后特征、降雨量和降雨量的6 階滯后特征，共四個特征值；在上游流域監(jiān)測點，共包含流域水位和流域水位的1 階滯后特征兩個特征值；在下游流域監(jiān)測點，包含流域水位一個特征值。

由于水庫水文數(shù)據(jù)包含水位、泄洪量、降雨量等數(shù)據(jù)，數(shù)值差異大，需要進行量綱的統(tǒng)一，以此提高模型的收斂速度和計算精度，因此，本文對序列的每一個特征值進行最大最小值歸一化處理，即：

式中：x 為預(yù)處理序列的取值；xmin為序列的最小值；xmax為序列的最大值。

2.4 評價指標

本研究所采用的評價指標包括平均絕對誤差（MAE，Mean Absolute Error）、均方誤差（MSE，Mean Square Error）、平均絕對百分比誤差（MAPE，Mean Absolute Percentage Error）三種，其定義如下所示：

2.5 實現(xiàn)細節(jié)

本研究的所有實驗均用Pytorch（V1.7.0）[10]進行實現(xiàn)。為找到最優(yōu)參數(shù)，我們使用網(wǎng)格搜索進行尋優(yōu)。具體而言，學習率在{0.01，0.005，0.0025}區(qū)間進行搜索；正則化系數(shù)在{1e-5，1e-7，1e-9}區(qū)間搜索；歷史時間長度在{2，4，16}中搜索；嵌入維度和LSTM 單元數(shù)在{8，16，32，64，128}區(qū)間搜索。此外，每次訓練執(zhí)行300 個epoch，并使用Adam[11]優(yōu)化器優(yōu)化模型參數(shù)。

2.6 實驗結(jié)果

2.6.1 對比分析

為驗證所提出的GSCPM 的有效性，我們將其與RNN、LSTM 和GRU 網(wǎng)絡(luò)進行對比，表4 展示了這4 種方法的對比結(jié)果?？梢钥闯?，GSCPM 在三個指標上均明顯優(yōu)于其他三個方法。在所有方法中，RNN 表現(xiàn)得最差，LSTM 在解決RNN 梯度爆炸問題后得到了性能提升是明顯的，GRU 在簡化LSTM 后也帶來一定的性能提升。然而，這些方法都無法考慮不同監(jiān)測點間的空間相關(guān)性，因而限制了性能。將GSCPM 與LSTM 對比，LSTM 在添加GCN 后帶來的性能提升是明顯的，說明了時空關(guān)聯(lián)預(yù)報對流域水位預(yù)報的重要性。

表4 GSCPM 及其對比方法的實驗結(jié)果

2.6.2 消融分析

為驗證空間依賴性和所構(gòu)建的滯后特征的有效性，本研究設(shè)置了4 個子方法來與進行對比：

1）GSCPM（-rain）：將降雨量滯后特征刪除后的GSCPM。

2）GSCPM（-flood）：將泄洪量滯后特征刪除后的GSCPM。

3）GSCPM（-up）：將上游水位滯后特征刪除后的GSCPM。

4）GSCPM（-all）：刪除所有滯后特征后的GSCPM。

表5 展示了所有子方法的實驗結(jié)果。從表中可看出，GSCPM 展示了最佳的性能和添加滯后特征的有效性。將GSCPM（-rain）、GSCPM（-flood）、GSCPM（-up）和GSCPM（-all）進行對比，它們?nèi)叨嫉玫搅瞬煌潭鹊男阅芴嵘f明所提議的三種滯后特征都能單獨提升模型性能。此外，將GSCPM 與GSCPM（-rain）、GSCPM（-flood）、GSCPM（-up）進行對比，GSCPM 仍然具有性能提升，說明滯后特征的組合使用會帶來更大的性能增益。

表5 GSCPM 及其子方法的實驗結(jié)果

2.6.3 可視化分析

為進一步理解所提出的GSCPM 模型，我們將該模型及對比模型在測試集上的預(yù)報結(jié)果進行了可視化，其呈現(xiàn)如圖11～圖18 所示。可以看出，該模型在向前一步預(yù)報的情況下，基本擬合出了未來的變化趨勢，只在少數(shù)突變情況下造成較大誤差。值得注意的是，下游水位的真實值波動劇烈，然而GSCPM 仍能捕捉到大部分波動，形成較為精準的預(yù)報。相比于其他三個方法的可視化結(jié)果，GSCPM 能更好地擬合真實水位曲線，也能更好地預(yù)報出突變情況，這也表明為什么我們的方法在眾多評價指標上表現(xiàn)得最好。

圖11 GSCPM 上游預(yù)報結(jié)果可視化

圖12 GSCPM 下游預(yù)報結(jié)果可視化

圖13 RNN 上游預(yù)報結(jié)果可視化

圖14 RNN 下游預(yù)報結(jié)果可視化

圖15 LSTM 上游預(yù)報結(jié)果可視化

圖16 LSTM 下游預(yù)報結(jié)果可視化

圖17 GRU 上游預(yù)報結(jié)果可視化

圖18 GRU 下游預(yù)報結(jié)果可視化

3 結(jié)論

本研究提出了一種圖指導的時空關(guān)聯(lián)預(yù)報網(wǎng)絡(luò)模型，用以針對性地解決流域洪水預(yù)報中的時空關(guān)系建模問題和滯后影響問題。實踐證明，該模型可以有效降低對流域數(shù)據(jù)維度的要求，并挖掘流域洪水跟降雨量、泄洪量、上游水位和下游水位之間復雜的時空關(guān)聯(lián)性。本文提供的是流域2 個節(jié)點的洪水預(yù)報模型，接下來將研究通過圖神經(jīng)網(wǎng)絡(luò)節(jié)點的擴展，通過一個圖神經(jīng)網(wǎng)絡(luò)模型，對整個流域的洪水進行預(yù)報的可行性問題，以促進流域洪水預(yù)報的廣泛實踐。