LSTM网络实现洪水预报原理附完整Matlab代码

1 模型原理

LSTM网络，即长短期记忆（Long Short-Term Memory） 网络。其概念最早于上世纪90年代被提出，是循环神经网络（RNN）的一种特殊形式，相较于RNN网络而言，LSTM网络能更好的处理序列中的长短期依赖信息。

1.1 整体结构

LSTM网络引入细胞状态，又称为记忆细胞，同时采用遗忘门、输入门和输出门三个门控单元来控制信息在网络中的流动，决定信息的保留、添加和遗忘等。对于LSTM的某一个细胞而言，输入包括*c**t*-1、*h**t*-1和*x**t*，其中*c**t*-1表示上一时刻的细胞状态，*h**t*-1表示上一时刻的输出值，*x**t*表示此刻的输入值，输出包括两个部分，分别是*c**t*和*h**t**，*其中*c**t*表示此时刻的细胞状态，承载着网络对序列信息的记忆，*h**t*表示此时刻模型的输出值。上一时刻的细胞状态*c**t*-1通过一系列的门控机制及线性运算调控信息的流动，决定信息的遗忘和保留，从而得到该时刻的细胞状态*c**t*输出到下一时刻，从而实现记忆更新。

图1 LSTM网络结构

1.2 遗忘门

遗忘门的作用对象是上一时刻的细胞状态*c**t*-1，是对细胞状态中的信息选择性遗忘，控制上一时刻的细胞状态*c**t*-1到下一时刻的细胞状态*c**t*中信息的保留程度，遗忘门的结构如图2所示。这一部分的输入值为*h**t*-1和*x**t*，输出值为*f**t*。遗忘门的计算见公式。

图2 遗忘门结构

式中，**f****t**表示输出值，σ表示sigmoid函数，sigmoid函数如图3所示，**W****f**表示权重参数，[*h**t*-1，*x**t*]表示由拼接的矩阵，**b****f**表示偏置参数。

图3 sigmoid函数

** ft的维度与c**t-1相同，ft计算结果的取值均在0-1之间，每一个数值都决定*ct*-1对应位置中信息的保留程度，数值越大，保留程度越高，其中0表示完全遗忘，1表示完全保留。

1.3 输入门

输入门同样作用于上一时刻的细胞状态*c**t*-1，是向细胞状态中存储新的信息，输入门的结构如图4所示。输入门可以分为两个部分，第一部分的输入值为*h**t*-1和*x**t*，输出值为*i**t*。这一部分的作用与遗忘门相似，同样是决定信息的保留程度，其计算见公式。

图4 输入门结构

式中，**i****t**表示输出值，σ表示sigmoid函数，**W****i**表示权重参数，[*h**t*-1，*x**t*]表示由拼接的矩阵，**b****i**表示偏置参数。
第二部分的输入值与第一部分的输入值相同，输出值为候选的细胞状态，其计算见公式。

式中，表示输出值，tanh函数如图5所示，**W****c**表示权重参数，[*h**t*-1，*x**t*]表示由拼接的矩阵，**b****c**表示偏置参数。

图5 tanh函数

得到的候选细胞状态能够保留*h**t*-1和*x**t*中的信息，将与**i****t**进行点乘计算，得到的计算结果来决定信息的保留程度，**i****t**的维度也与*c**t*-1相同。同样地，**i****t**中数值越大，保留程度越高，其中0表示完全丢弃，1表示完全保留，最终保留下的信息将添加到新的细胞状态中。

1.4 记忆更新

当前时刻输出的新的细胞状态*c**t*由两部分组成，其中一部分是上一时刻的细胞状态*c**t*-1经过遗忘门后遗留下来的信息，另一部分是输入门中输出的候选细胞状态中新添加的信息，结构图如图6所示。新的细胞状态*c**t*计算见公式。式中各变量的含义与前述相同。

图6 记忆更新

1.5 输出门

输出门以新的细胞状态*c**t*为基础，决定最终的输出值*h**t*，输出门的结构如图7所示。输出门也同样分为两个部分，第一部分的输入值为*h**t*-1和*x**t*，输出值为*o**t*。其计算见公式。

图7 输出门结构

 式中，**o****t**表示输出值，σ表示sigmoid函数，**W****o**表示权重参数，[*h**t*-1，*x**t*]表示由拼接的矩阵，**b****o**表示偏置参数。
第二部分是将新的细胞状态*c**t*和**o****t**进行点乘计算后通过tanh函数进行非线性运算，将输出值变换到-1至1的区间范围内，以此决定新的细胞状态*c**t*中的输出部分，其计算见公式：

式中各变量的含义与前述相同。

2 特征因子的选择

利用神经网络来进行洪水预报时，由于不涉及到物理机制，因此特征因子的选择至关重要。最常采用的因子包括前期降雨量和前期流量，例如在预报t+T（T为预见期）时刻A水文站的流量时，t-a、t-a+1、......t-1、t（a+1为前期时段长度）时刻的降雨量即为前期降雨量，t-a、t-a+1、......t-1、t时刻的流量即为前期流量。其中，降雨量可以是整个流域的平均面雨量，或者对流域分区后计算各个区域的平均面雨量，也可以是多个雨量站的降雨量，根据流域特性自行选择。流量可以采用预报水文站的前期流量，也可以是预报水文站上游的前期流量。除了降雨量和流量以外，蒸发量、气温等也可以成为特征因子。本文将流域平均面雨量和预报水文站流量作为特征因子，以5作为前期时段长度，以此为例介绍模型方法和参考代码。

2.1 模型参数

Matlab中有自带的BP神经网络工具箱，无需自己手写代码，只需要设置好相关的参数即可，包括模型层数、隐藏单元个数、batch大小、损失函数、初始学习率和最大迭代次数等。

2.2 数据预处理

将数据输入进模型前，先要将降雨和流量数据分别进行标准化或归一化处理，便于模型训练。在得到预测结果后，要记得进行反标准化或反归一化（Matlab工具箱有可以直接实现此功能的工具）。

3 注意事项

有了Matlab的工具箱后，代码其实非常简单，但要注意的是输入进去的数据格式。下方代码中XTrain指的是训练集中的前期流量、前期降雨量，YTrain指的是训练集中的预测流量。在这里，训练集中有1000个样本，前期流量和前期降雨量各有5个时段，因此XTrain的size就是10*1000，也就是每一列为1个样本，YTrain的size就是1*1000（如果要同时预测T个时段，那就是T*1000）。XTest也是一样。

XTrain数据格式，P指降雨，Q指流量，滑动窗口为1来构造样本：
P1P2P3......PnP2P3P4......Pn+1..............................
P5
P6P7......Pn+4Q1Q2Q3......Qn..............................Q4Q5Q6......Qn+3Q5Q6Q7......Qn+4
YTrain数据格式，只预测1个时段时：
Q6Q7Q8......Qn+5

4 Matlab代码

clc,clear
%读取降雨量和流量数据
load Ystation.mat
dataTrain = PQ(1:1000,:);%前1000个作为训练数据
dataTest = PQ(1001:1461,:);%后续作为验证数据
window = 5;%例：用5个时段预测
step = 1;%只预测未来1个时段
for i = 1:length(dataTrain)-window-step+1;%间隔一个单位进行滑动
    contem = dataTrain(i:i+window-1,1:2);
%这里因为只有2个特征因子，所以是1:2，如果有n个，就是1:n
    XTrain(:,i) = [contem(:)]%包含前期流量和前期降雨量
    YTrain(:,i) = dataTrain(window+i+step-1,1);%预测时段的真实流量值
end
for i = 1:length(dataTest)-window-step+1
    contem2 = dataTest(i:i+window-1,1:2);
    XTest(:,i) = [contem2(:)];
end
%数据可以提前做好标准化或归一化，这个工具箱也有相关代码可实现这一操作
%确定模型基本参数
numFeatures = (window+1).*1;%XTrian的行数，即特征数
numOutputs = 1;%YTrain的行数，即输出结果的维度
numHiddenUnits = 40;
%搭建模型
layers = [ ...
    sequenceInputLayer(numFeatures)    %输入层
    lstmLayer(numHiddenUnits)  % lstm层，如果是构建多层的LSTM模型，可以修改。
    fullyConnectedLayer(numOutputs)    %为全连接层,是输出的维数。
    regressionLayer];      %其计算回归问题的半均方误差模块
options = trainingOptions('adam', ...
    'MaxEpochs',300, ...
    'GradientThreshold',1, ...
    'InitialLearnRate',0.001, ...      
    'LearnRateSchedule','piecewise', ...%每经过一定数量的epoch，学习率乘一个系数。
    'LearnRateDropPeriod',700, ...      
    'LearnRateDropFactor',0.15, ...     
    'Verbose',0,  ...  
    'Plots','training-progress',...
    'MiniBatchSize', 32);   
net = trainNetwork(XTrain,YTrain,layers,options);
net = predictAndUpdateState(net,XTrain);
YTrainPred = [];
numTimeStepsTrain = numel(XTrain)./numFeatures;
numTimeStepsTest = numel(XTest)./numFeatures;
for i = 2:numTimeStepsTest
    [net,YTestPred(i)] = predictAndUpdateState(net,XTest(:,i),'ExecutionEnvironment','cpu');
end

以上内容仅供参考，如有错误，欢迎批评指正！