基于LSTM网络的实时入库洪水预报方法

Real-Time Reservoir Flood Forecasting Method Based on LSTM Network

摘要: 随着流域水文气象资料观测手段的进步和数据的不断积累，采用数据驱动的方式来研究水库流域入库径流问题，已成为一种可行的思路。本文建立了一种基于LSTM (Long Short-term Memory)网络深度学习算法的入库洪水预报模型，该模型以水库前期入库流量序列、流域前期降雨序列和上游站前期出库流量序列等多个维度作为模型输入，水库预报入库流量作为输出，进行实时洪水预报。安康水库洪水预报的实例验证表明，LSTM网络模型可以很好地模拟入库洪水过程，并在实时洪水预报中也有较好的表现。本研究可为水库防洪及兴利调度提供技术支撑。

Abstract: With the progress of hydrological and meteorological data observation and the accumulation of data, it has become a feasible way to study the reservoir inflow forecasting by data driven method. This paper established a flood forecasting model based on LSTM (long short term memory) network and depth learning algorithm. The model takes the reservoir early stage inflow, rainfall and upstream sequence as the inputs, and the reservoir inflow as output. The application of Ankang reservoir shows that LSTM network model can simulate the flood hydrograph very well, and has good performance in real-time flood forecasting. This study can provide technique support for reservoir flood control and beneficial operation.

1. 引言

水库的洪水预报系统是整个水库防洪和发电调度决策系统的核心组成部分，可靠的洪水预报是做好水库防洪及兴利调度的前提。入库洪水预报通常基于流域降雨、蒸发、下垫面条件等建立的洪水预报模型来进行 [1]。目前，洪水预报模型主要分为3种类型：一是通过对流域产汇流机理与过程进行概化描述而建立的概念性或集总式水文模型，二是基于数学物理方程、以空间离散方式详细描述径流形成与转化过程的分布式水文模型，三是通过建立相关前期影响因素与当前径流要素之间的“输入–输出”关系进行预测的数据驱动模型 [2]。

概念性水文模型和分布式水文模型都是依据历史数据率定模型参数然后用于预报，由于模型结构及参数所限，往往难以完全反映流域水文规律，且所需的下垫面条件等数据不易实时获取。同时，流域产汇流过程是高度复杂和不确定的，复杂的水文模型并不一定能显著提高预报精度，还会遇到“异参同效”等问题 [3]。

随着信息技术的进步和推广应用，流域水文气象资料的观测手段也已取得了很大提升，一些学者尝试采用数据驱动的方式来研究流域产汇流问题。数据驱动方法不着重于描述具体的水文物理机制，而是依赖于历史观测资料，特别适合已积累了较多实测数据、但流域产流机制复杂难以理清的水库入库洪水预报。

早在20世纪90年代初期，人工神经网络(Artificial Neural Network, ANN)就被应用于降雨径流模拟当中 [4] [5]。近几年来，由于递归神经网络(Recurrent Neural Network, RNN)和长短期记忆网络(Long-Short Term Memory, LSTM)结构特别适用于时间序列数据的关系模拟，也被引入流域降雨产流模拟的研究中。但是，在水文预报领域，基于LSTM网络和深度学习的洪水预报研究才刚刚兴起，不同的研究者在不同的流域、从不同的角度做了许多探索性的尝试 [6] [7] [8] [9]。

可以预期，随着大数据分析技术和深度学习智能算法的不断完善，基于大数据、传统水文模型以及深度学习算法相结合的洪水预报理论及方法将会逐步发展起来，并在未来有着广阔的应用前景 [10] [11] [12]。

2. 研究方法

长短期记忆网络是一种在循环神经网络基础上构建长时记忆结构的特殊神经网络。RNN是一种具有重复神经网络模块的链式网络，但其在处理长序列数据时会存在梯度消失、梯度爆炸和时序相关性较差等问题，而LSTM在RNN的基础上增加了输入与输出的控制转换规则来处理记忆单元状态，将先前的序列样本特征与当前时间步长的输入信息建立相关性，从而在一定程度上弥补了循环神经网络的缺陷。

LSTM网络的详细结构如图1所示，网络内部结构通常由输入层、多个LSTM单元以及输出层组成。其中多个LSTM单元共同组成了网络中的隐藏层，每个单元内部拥有三个已定义的门控结构，用于提取和保存时序特征，并通过逻辑运算调整各单元的相关状态 [8]。

Figure 1. Internal structure diagram of LSTM network

图1. LSTM网络内部结构示意图

设输入的时间序列样本为 $X=\left\{x_0,x_1,\cdots ,x_{n+1}\right\}$，定义时间步长为n + 1，则对于任意时间步n，其输入 $x_n$ 将由输入层保存到LSTM计算单元内。另外，网络会将前一个时间步n − 1得到的单元状态(Cell State) $C_{n-1}$ 与输出值(HiddenState) $h_{n-1}$ 同时作为输入应用到本次的计算流程中。计算单元内的遗忘门负责从输出单元状态 $h_{n-1}$ 与当前输入 $x_n$ 中决定删除或保留哪些信息，其计算公式如下：

$f_n=\sigma \left(W_{fh}\cdot h_{n-1}+W_{fx}\cdot x_n+b_f\right)$ (1)

其中： $W_{fh}$ 和 $W_{fx}$ 分别为遗忘权重矩阵和遗忘隐藏权重矩阵， $b_f$ 是遗忘门的偏置，而 $\sigma$ 为sigmoid函数且 $\sigma \left(x\right)={\left(1+e^{-x}\right)}^{-1}$，该激活函数会将输出矩阵中的各元素取值范围将被约束在区间[0, 1]之间，通过与 $C_{n-1}$ 按位乘操作后，从而决定了长期记忆流的保留程度。

输入门负责将哪些信息添加到当前时间步的细胞状态 $C_n$ 中，具体计算公式如下：

$i_n=\sigma \left(W_{ix}\cdot x_n+W_{ih}\cdot h_{n-1}+b_i\right)$ (2)

${\tilde{C}}_n=\tanh \left(W_{cx}\cdot x_n+W_{ch}\cdot h_{n-1}+b_c\right)$ (3)

其中： $W_{ih}$ 和 $W_{ix}$ 分别为输入门中不共享参数的权重矩阵， $b_c$ 是输入门的偏置。通过将 $x_n$ 与 $h_{n-1}$ 带入式(2)可计算求得输入门矩阵或向量 $i_n$，用于决定候选记忆单元 $C_n$ 的更新程度。 $C_n$ 由式(3)求得，其计算原理与门限控制单元相同，并在最后使用双曲正切激活函数代替sigmoid用于对状态和输出进行限定，使矩阵中各元素的取值范围约束在[−1, 1]区间内，从而使候选记忆单元 $C_n$ 中各元素的值便于后续处理。

接下来，LSTM单元状态会通过公式(4)进行参数更新。首先上一时刻的记忆单元 $C_{n-1}$ 与遗忘门计算得到的 $f_n$ 临时变量之间进行逐元素乘法运算，从而可将无用的长时信息遗忘且维持有用信息，同理 $i_n\cdot {\tilde{C}}_n$ 可以保留当前时间步中输入的有效信息，最后对齐进行加法运算，使得历史信息与当前信息共同保留在新的矩阵 $C_n$ 中，此时更新LSTM中时间步长为n的记忆单元状态为 $C_n$。

$C_n=f_n\cdot C_{n-1}+i_n\cdot {\tilde{C}}_n$ (4)

算法最后需要根据记忆单元状态求出输出值 $h_n$。该过程同样分为两个步骤，具体如下所示：

第一步：输出门会读取当前输入 $x_n$ 和上一时间步长的输出值 $h_{n-1}$，经过sigmoid激活函数得到一个值为0~1之间的激活向量或矩阵，以此来决定输入值以及上一时间步长的输出值中哪些部分可以输出；

第二步：利用tanh函数得到−1~1范围之间的激活向量或矩阵来处理LSTM单元当前的记忆状态单元 $C_n$，并与第一步输出的结果进行矢量逐元素相乘，最终得到输出 $h_n$，该值中包含有n个时间步长内的所有隐层信息。

$o_n=\sigma \left(W_{ox}\cdot x_n+W_{oh}\cdot h_{n-1}+b_o\right)$ (5)

$h_n=o_n\cdot \tanh \left(C_n\right)$ (6)

因此LSTM模型在洪水预报中的应用可以看作是较为复杂的时间窗函数，依赖给定的历史洪水数据，并以有监督的训练方式来迭代更新计算单元内的权重参数，最终得到最拟合训练数据的一种表征算法，并进一步将其应用于新数据的预测。

3. LSTM预报模型的建立

3.1. 模型结构

LSTM网络同其他循环神经网络一样，需要根据样本量大小和数据稀疏程度来合理搭建网络结构，再传入相应的数据进行批处理计算。本研究选用TensorFlow2.0深度框架来搭建LSTM网络模型，可以直接调用其高级的API接口Keras。整个网络模型搭建过程可分为以下几个步骤：

1) 数据预处理

研究中使用的数据主要是流量和降雨数据，由于不同类型的数据取值范围差异较大，为使算法更快且更易正确地收敛到最优解，使模型具有更好的模拟性能，需要对数据进行归一化处理，归一化处理一般有0均值归一化和最大最小值归一化两种方法，其计算公式见下：

$最大最小值归一化：X_{norm}=\frac{X-X_{\min }}{X_{\max }-X_{\min }}$ (7)

$0均值归一化：Z=\frac{X-\mu }{\sigma }$ (8)

经过多次试验，发现最大最小值归一化的建模能力更强，且使用该类型线性的变换后对原始数据中洪峰和水量平衡的影响更小。因此本文最终采用的是最大最小值归一化公式(7)。

2) 输入层参数的选取

神经网络的输入层往往决定了网络的预测能力，通过选用不同的输入层参数组合，上游水库的出库径流、流域平均降雨、流域平均蒸发、之前时刻的入库径流等，不同的组合方法具有不同的效果。

3) LSTM网络结构的搭建

网络结构主要是指网络层数以及不同层数的大小，先后尝试过一层网络、二层网络和三层网络，经过多次试验后，发现两层LSTM层加一层全连接层的网络模拟效果最好，激活函数选用ReLU。

4) 损失函数的选取

深度学习中可供选择的损失函数有多种，本次研究最终选择了较为常用的两种，均方误差损失(Mean Squared Error, MSE)和平均绝对百分误差(Mean Absolute Percentage Error, MAPE)。其误差计算公式见下：

$\text{MSE}=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(y-y_{per}\right)}^2}$ (9)

$\text{MAPE}=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}\left|x_i-m\left(x\right)\right|}$ (10)

同样的，通过多次试验，最终选择了表现较好的均方误差损失(MSE)作为模型的损失函数。

5) 优化器的选择

合适的优化器可以减少参数更新带来的变化，从而带来更加稳定的收敛。本次研究选择了较为常见的两种，分别是RMSProp算法和Adam算法。经过多次试验后发现差别并不大，最终选择了相对稳定的Adam算法。

3.2. 模型参数

在LSTM网络的计算中往往还需要通过堆积LSTM隐藏层来模拟较为复杂的水文时间序列。本研究在使用的双层LSTM计算结构如图2所示，输入前n个时段的特征值向量，计算第n + 1个时段的目标值。

Figure 2. Computing structure diagram of LSTM network

图2. LSTM网络计算结构图

经过比选，研究中采用安康水库前期入库流量序列、前期面雨量序列和上游站(石泉水库)前期出库流量序列3个维度作为模型输入，以下一时段安康预报入库流量作为输出。

为使得模型能够正常运行，还需要初始化其他的一些参数。单批训练样本个数为50，为保证充分收敛，迭代代数设置为300。

3.3. 评价指标

预报计算过程与实测过程之间的吻合程度可用确定性系数DC来评定。DC是一个标准化统计指标，其值介于(−∞, 1]之间，等于1表示模型结果完美拟合实测值。计算公式如下：

$DC=1-\frac{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left[y_c\left(i\right)-y_0\left(i\right)\right]}^2}}{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left[y_0\left(i\right)-{\bar{y}}_0\right]}^2}}$ (11)

式中：DC为确定性系数，又称纳西效率系数(NSE)； $y_0\left(i\right)$ 为实测值(m³/s)； $y_c\left(i\right)$ 为计算值(m³/s)； ${\bar{y}}_0$ 为实测值的均值(m³/s)；n为资料序列长度。

另外，还选择洪峰相对误差和洪量相对误差两个指标来评定预报精度。根据《水文情报预报规范(GB/T22482-2008)》 [13] 及安康预报工作实际状况，研究中以确定性系数 ≥ 0.7、洪峰相对误差 ≤ 20%、洪量相对误差 ≤ 20%作为合格要求。

4. 研究区域及数据

4.1. 研究区概况

安康水库位于汉江上游陕西省境内，是一座以发电为主，兼顾防洪、航运等综合利用的大型水利枢纽。坝址位于安康市城西18 km处，控制流域集水面积35,700 km²。地理上属南北过度的亚热带地区，夏季受太平洋副高控制，多暴雨。汉江安康水库坝址以上洪水特别频繁、陡涨陡落，同流域的地理条件以及形成洪水的暴雨移动方向密切相关。由于流域内水系条件较为特殊，常遇降雨往往形成稀遇洪水。

汉江上游支流众多，水系分布不对称，安康水库上游共有河流221条(不含汉江干流)，其中一级支流62 条，二级以下支流159条。流域面积200 km²以上的一级支流共26条，其中褒河、湑水河、子午河、牧马河、任河、岚河6条较大一级支流流域面积均在2000 km²以上，年径流量都超过10亿m³。任河是汉江石泉–安康区间的最大支流，流域面积4870 km²，年均径流量41.3亿m³，在重庆市城口县境内建有巴山水电站工程。岚河流域面积2130 km²，年均径流量13.50亿m³，干流上建有控制性工程蔺河水库。石泉–安康流域主要水系情况，可概化如图3所示，其中括号内的数字为多年平均流量(m³/s)。

Figure 3. General map of main rivers between Shiquanand-Ankang section

图3. 石泉–安康区间主要水系概化图

安康水库的入库流量组成大体可划分为2个部分，即：石泉以上流量和石泉–安康区间流量。由于石泉以上部分的洪水受石泉水库控制，且喜河水库又无调节能力，故石泉出库流量即可代表石泉上游的产流。而石泉–安康区间支流众多，少数支流上虽设有水文站，但在做实时预报时无法及时获取其数据，故研究中暂未将其纳入模型计算。因此，整个区间按一个产流单元考虑。

4.2. 数据资料

研究中收集到安康水库2010~2019年间共21场的洪水资料，包括石泉–安康区间流域内各遥测站的降雨量、石泉水库出库流量及安康水库入库流量等。各测站雨量按算术平均法转换为区间面雨量，流量均按1 h整理。为了率定参数和检验效果，将前15场洪水用于模型训练，后面6场洪水用于检验。

4.3. 入库流量平滑处理

在数据整理过程中，发现安康水库逐小时入库流量资料锯齿状严重，有的时段甚至出现负值，与实际物理过程不符，通过下面方式进行平滑修正处理：

1) 首先进行3 h滑动平均处理，如果处理之后还不满足要求，就继续做4 h平滑处理，直到流量曲线相对光滑，一般最长不会超过6 h；

2) 然后将滑动之后的流量进行修正，保证滑动之后的洪峰流量与原洪峰流量发生时间是同一时刻；

3) 最后，缩放处理洪峰流量使之与原洪峰流量保持一致，并保证整体水量平衡，即该场次洪水的洪量修正前后也保持一致。

5. 精度评价及分析

在模型训练时，所有时段的安康入库流量、区间面雨量和石泉出库流量均为已知值，但在实时预报时，后期的安康入库是未知的。为便于区别两种计算模式，本文将前者称为模拟计算，后者称为预报计算。模拟计算的任务主要是根据建立的预报模型，输入实际发生的降雨量、前期入库流量等数据，通过深度学习后，得到LSTM洪水预报模型神经元参数并保存为记忆文件，以此参数作为预报计算的基础。预报计算时，输入区间面雨量和石泉出库流量，而安康入库则滚动使用前期预报计算得到的流量。因此，由于两种计算模式的输入不同，所以结果的精度会有差异。

5.1. 模拟计算精度评价

5.1.1. 率定期结果

对率定期的15场洪水模拟计算，确定性系数DC最高为0.993，最低为0.919，平均为0.967。洪峰误差最大值为−15.1%，最小值为0.2%，误差绝对值平均为4.8%，误差小于20%的有15个场次，占100%。洪量误差最大值为3.6%，最小值为0.03%，误差绝对值平均为1.4%，误差小于20%的有15个场次，占100.0%。整体模拟精度较好，其中20150629号洪水模拟结果见图4。

Figure 4. Simulation results of No. 20150629 flood

图4. 20150629号洪水模拟计算结果图

5.1.2. 验证期结果

对验证期的6场洪水模拟计算，确定性系数DC最高为0.981，最低为0.783，平均为0.931。洪峰误差最大值为−34.4%，最小值为−1.8%，误差绝对值平均为11.6%，误差小于20%的有5个场次，占83.3%。洪量误差最大值为−8.3%，最小值为1.0%，误差绝对值平均为3.8%，误差小于20%的有6个场次，占100%。整体模拟精度较率定期有所下降，其中20190915号洪水模拟结果见图5。

Figure 5. Simulation results of No. 20190915 flood

图5. 20190915号洪水模拟计算结果图

从率定期和验证期的模拟计算结果来看，建立的LSTM网络洪水预报模型能很好地模拟历史洪水过程，接下来将其应用于实时洪水预报计算。

5.2. 预报计算精度对比分析

实时预报计算时，除初始几个时段外，安康入库改为输入前期预报计算得到的流量。将全部21场洪水均用来进行预报计算，确定性系数DC最高为0.991，最低为0.371，平均为0.863。洪峰误差最大值为−46.9%，最小值为−0.2%，误差绝对值平均为11.5%，误差小于20%的有16个场次，占76.2%。洪量误差最大值为28.6%，最小值为0.7%，误差绝对值平均为8.6%，误差小于20%的有19个场次，占90.5%。精度较模拟计算有所下降，但整体上也可满足安康水库实时预报精度的要求。模拟计算与预报计算的误差比较情况如表1所示。

由于模拟计算的区间面雨量、石泉出库、前期安康入库等输入均为已发生值，故LSTM网络模型能很好地模拟历史洪水过程。而预报计算时，安康入库的输入改用前期预报得到的流量，与实际值本就存在一定误差，因此预报精度有所降低，特别是洪峰相对误差偏大。

预见期方面，在模拟计算时，用前n个时段的多维输入预报第n + 1个时段的入库流量，故预见期只有一个时段。而预报计算时，利用预报流量不断向后滚动计算，可以明显延长预见期。另一个优势是，在实际预报时LSTM模型可以直接利用前期多个时段的已发生数据，不断对后期预报进行修正，进行实时滚动预报。

Table 1. Comparison table of effects of flood simulation and forecast calculation

表1. 场次洪水模拟与预报计算效果对比表

6. 结论

本文基于LSTM网络建立了入库洪水预报模型，并应用于安康水库实际洪水预测。确定了模型输入、输出、隐层结构及主要参数，比较了有无实际前期入库输入影响的两种预报计算模式，得到主要结论如下：

1) LSTM网络模型可以很好地模拟逐小时入库洪水过程，模拟计算精度较高。

2) 实时预报时将前期计算得到的入库流量作为后期预报的输入，预报计算精度会有所下降。

3) 基于LSTM网络的实时入库洪水预报方法具有一定实用性。

当然，受限于实时预报时获取资料的时效问题，本研究依然存在不足。如模型可进一步增加输入维度，水库区间可考虑划分为多个子流域等，这些都有待下一步深入开展研究。

参考文献