使用机器学习预测莱茵河水位

Long Short Term Memory (LSTM) models are a powerful type of neural network ideally suited to predict time-dependent data. Rhine water levels fit right into this category: they vary over time, depending on a range of variables such as rain, temperatures and snow cover in the Alps.

长短期记忆(LSTM)模型是一种功能强大的神经网络，非常适合预测与时间有关的数据。 莱茵河的水位正好属于这一类：它们随时间而变化，取决于阿尔卑斯山的降雨，温度和积雪等一系列变量。

The Rhine is Europe’s lifeblood. For centuries it has been used as a major artery for shipping goods into Germany, France, Switzerland and Central Europe. However, with climate change, water levels on the river are likely to become more variable. Forecasting the river’s level accurately is therefore a primary concern for a whole range of actors, from shipping companies to commodity traders and industrial conglomerates.

莱茵河是欧洲的命脉。 几个世纪以来，它一直被用作将货物运往德国，法国，瑞士和中欧的主要动力。 然而，随着气候变化，河流上的水位可能变得更加可变。 因此，准确地预测河流的水位是从航运公司到商品贸易商和工业集团等所有参与者的首要考虑。

A barge laden with coal navigating the Rhine (source: https://commons.wikimedia.org/wiki/File:Coal_barge_Chilandia_on_Rhine_-_looking_south.jpg) 一艘载有煤炭的驳船驶向莱茵河(来源： https : //commons.wikimedia.org/wiki/File : Coal_barge_Chilandia_on_Rhine _-_ looking_south.jpg )

Unlike classical regression-based models, LSTMs are able to capture non-linear relationships between different variables; more precisely, the sequence dependence among these variables. This blog focuses on the problem of Rhine river forecasting using LSTMs, rather than the theory behind these models.

与经典的基于回归的模型不同，LSTM能够捕获不同变量之间的非线性关系。 更确切地说，这些变量之间的序列依赖性。 该博客主要关注使用LSTM预测莱茵河的问题，而不是这些模型背后的理论。

眼前的问题 (Problem at hand)

The problem we are looking to solve here is the following: we would like to forecast next-day water levels at Kaub, a key chokepoint in western Germany, with the highest possible accuracy.

我们要在此处解决的问题如下：我们希望以尽可能最高的精度预测第二天德国西部的主要阻塞点Kaub的水位。

We have historical daily data from 2 January 2000 to 27 July 2020, equivalent to 7513 observations. The dataset includes 15 different categories, displayed as columns:

我们拥有2000年1月2日至2020年7月27日的每日历史数据，相当于7513次观测。 数据集包括15个不同的类别，显示为列：

‘date’: the date of the observation

“日期”：观察日期 ‘Kaub’: the day-on-day difference in Kaub water level, in centimetres — this is the ‘y’ value we are trying to forecast (source: WSV)

“ Kaub”：Kaub水位的每日差异，以厘米为单位-这是我们要预测的“ y”值(来源：WSV)‘Rheinfelden’: the absolute value of the water flow at Rheinfelden, in Switzerland, in cubic metres per second (source: BAFU)

“莱茵费尔登”：瑞士莱茵费尔登的水的绝对值，以立方米/秒为单位(来源：BAFU) ‘Domat’: the absolute value of the water flow at Domat, near the source of the Rhine, in cubic metres per second (source: BAFU)

“ Domat”：靠近莱茵河源头的Domat水流的绝对值，以立方米/秒为单位(来源：BAFU) ‘precip_middle’: the average daily amount of rain recorded at 20 weather stations along the Rhine, in millimetres (source: DWD)

'precip_middle'：沿莱茵河的20个气象站记录的日平均降雨量，以毫米为单位(来源：DWD) ‘avgtemp_middle’: the average temperature recorded at the same stations, in degrees Celsius

'avgtemp_middle'：同一站记录的平均温度，以摄氏度为单位 ‘maxtemp_middle’: the maximum temperature recorded at the same stations

'maxtemp_middle'：同一站记录的最高温度‘mintemp_middle’: the minimum temperature recorded at the same stations

'mintemp_middle'：同一站记录的最低温度‘precip_main’: the average daily amount of rain recorded at 8 weather stations along the Main, a major tributary of the Rhine, in millimetres (source: DWD)

'precip_main'：沿美因河(莱茵河的主要支流)的8个气象站记录的日平均降雨量，以毫米为单位(来源：DWD)‘avgtemp_main’: the average temperature recorded at the same stations, in degrees Celsius

'avgtemp_main'：同一站点记录的平均温度，以摄氏度为单位 ‘maxtemp_main’: the maximum temperature recorded at the same stations

'maxtemp_main'：在相同站点上记录的最高温度‘mintemp_main’: the minimum temperature recorded at the same stations

'mintemp_main'：同一站记录的最低温度‘precip_neckar’: the average daily amount of rain recorded at 7 weather stations along the Neckar, also a major tributary of the Rhine, in millimetres (source: DWD)

“ precip_neckar”：内卡河(也是莱茵河的主要支流)内的7个气象站记录的日平均降雨量，以毫米为单位(来源：DWD)‘avgtemp_neckar’: the average temperature recorded at the same stations, in degrees Celsius

'avgtemp_neckar'：同一站记录的平均温度，以摄氏度为单位 ‘maxtemp_neckar’: the maximum temperature recorded at the same stations

'maxtemp_neckar'：在同一站点记录的最高温度‘mintemp_neckar’: the minimum temperature recorded at the same stations

'mintemp_neckar'：在同一站点记录的最低温度

Note that the choice of variables is entirely mine and is based on my experience dealing with Rhine analysis. Selecting the right inputs is one of the most important steps in time series analysis, whether you are using classical regression models or neural networks. If you select too few variables, the model may not capture the full complexity of the data (this is called underfitting). By contrast, if you choose too many inputs, the model is likely to overfit the training set. This is bad as well as it could mean the model struggles to generalise to a new dataset, which is essential for predictions.

请注意，变量的选择完全是我的，并且基于我对莱茵河分析的经验。 无论您使用的是经典回归模型还是神经网络，选择正确的输入都是时间序列分析中最重要的步骤之一。 如果选择的变量太少，则该模型可能无法捕获数据的全部复杂性(这称为欠拟合)。 相反，如果您选择太多的输入，则该模型可能会过度拟合训练集。 这是不好的，也可能意味着模型难以将其推广到一个新的数据集，这对于预测至关重要。

First, let’s load all the libraries we will need for this exercice:

首先，让我们加载此练习所需的所有库：

import datetimeimport matplotlib as mplimport matplotlib.pyplot as pltimport numpy as npimport pandas as pdimport seaborn as snsimport tensorflow as tffrom sklearn.preprocessing import LabelEncoderfrom sklearn.preprocessing import StandardScalerfrom sklearn.metrics import mean_squared_errorimport joblib

Below is a sample of the first few lines of the dataset. We can load it easily with the Pandas library:

以下是数据集前几行的示例。 我们可以使用Pandas库轻松加载它：

# first, we import data from excel using the read_excel functiondf = pd.read_excel('RhineLSTM.xlsx')# then, we set the date of the observation as the indexdf.set_index('date', inplace=True)df.head()

Once loaded, we can plot the dataset using the Matplotlib library:

加载后，我们可以使用Matplotlib库绘制数据集：

# specify columns to plotcolumns = [0, 1, 2, 3, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]i = 1values = df.values# define figure object and sizeplt.figure(figsize=(9,40))# plot each column with a for loopfor variable in columns: plt.subplot(len(columns), 1, i) plt.plot(values[:, variable]) plt.title(df.columns[variable], y=0.5, loc='right') i += 1plt.show()

It’s also generally a good idea to plot histograms of the variables:

通常，绘制变量的直方图也是一个好主意：

# histograms of the variablesdf.hist(figsize=(9,18))plt.show()

Using the Seaborn library, you can create a violin plot to understand the distribution of each variable:

使用Seaborn库，您可以创建一个小提琴图来了解每个变量的分布：

# calculate dataset mean and standard deviationmean = df.mean()std = df.std()# normalise dataset with previously calculated valuesdf_std = (df - mean) / std# create violin plotdf_std = df_std.melt(var_name='Column', value_name='Normalised')plt.figure(figsize=(12, 6))ax = sns.violinplot(x='Column', y='Normalised', data=df_std)_ = ax.set_xticklabels(df.keys(), rotation=90)

评估不同的模型 (Evaluating different models)

I have trained different types of models on the Rhine dataset to establish which one fits best:

我已经在Rhine数据集上训练了不同类型的模型，以确定哪种模型最合适：

The baseline model, also known as persistence model, returns the current rate of change in Kaub water level as the prediction (essentially predicting “no change”). This is a reasonable baseline as Rhine water levels typically change over a number of days due to wider weather phenomena (eg, slow melt in the Alps progressively making its way downstream).

基线模型(也称为持久性模型)将Kaub水位的当前变化率作为预测值(实质上是预测“无变化”)。 这是一个合理的基准，因为莱茵河的水位通常会由于天气现象的变化而在几天内发生变化(例如，阿尔卑斯山的缓慢融化逐渐向下游移动)。 The simplest model you can train assumes a linear relationship between the input variables and the predicted output. Its main advantage over more complex models is that it is easy to interpret, however it performs only marginally better than the baseline network.

您可以训练的最简单模型假设输入变量和预测输出之间存在线性关系。 与更复杂的模型相比，它的主要优点是易于解释，但其性能仅比基线网络好一点。 A dense network is more powerful, but cannot see how the input variables are changing over time. This shortcoming is addressed by multi-step dense and convolution neural networks, which take multiple time steps as input for each prediction.

密集的网络功能更强大，但无法看到输入变量随时间的变化。 多步密集和卷积神经网络解决了该缺点，该网络将多个时间步长作为每个预测的输入。 The LSTM model comes out on top, with a lower mean absolute error rate on the validation and test sets.

LSTM模型名列前茅，在验证和测试集上具有较低的平均绝对错误率。

The Python code required to create this chart is too long for this blog, but you can access it here, applied to a different dataset.

创建此图表所需的Python代码对于此博客而言太长，但是您可以在此处访问它，并将其应用于其他数据集。

LSTM回归模型 (LSTM model for regression)

LSTM networks are a form of recurrent neural network that can learn long sequences of data. Instead of neurons, they are made of memory blocks connected to each other via layers. A memory block contains gates (input, forget, output) that manage its state and output, and enable it to be smarter than a typical neuron.

LSTM网络是递归神经网络的一种形式，可以学习较长的数据序列。 它们不是神经元，而是由通过层相互连接的存储块组成。 一个存储块包含管理其状态和输出的门(输入，忘记，输出)，并使之比典型的神经元更智能。

They have been used extensively in academic circles to forecast river height and have been proven to outperform classical hydrological models in certain situations.

它们已在学术界广泛用于预测河流高度，并已证明在某些情况下优于经典水文模型。

Let’s start all over again with the code from the beginning of this article to load the Rhine database:

让我们从本文开头的代码重新开始，以加载Rhine数据库：

import datetimeimport matplotlib as mplimport matplotlib.pyplot as pltimport numpy as npimport pandas as pdimport seaborn as snsimport tensorflow as tffrom sklearn.preprocessing import LabelEncoderfrom sklearn.preprocessing import StandardScalerfrom sklearn.metrics import mean_squared_errorimport joblib# first, we import data from excel using the read_excel functiondf = pd.read_excel('RhineLSTM.xlsx', sheet_name='Detailed4_MAIN’)# then, we set the date of the observation as the indexdf.set_index('date', inplace=True)

Data preparation

资料准备

To build a functioning LSTM network, the first (and most difficult) step is to prepare the data.

要构建运行良好的LSTM网络，第一步(也是最困难的一步)是准备数据。

We will frame the problem as predicting today’s rate of change in Kaub water level (t) given the weather and Swiss upstream flows of today and the previous 6 days (backward_steps = 7).

考虑到今天和前6天的天气和瑞士的上游流量，我们将以此框架来预测今天的考布水位(t)的变化率。

The dataset is standardised using the StandardScaler() function in the Scikit-Learn library. For each column of the dataframe, each value in the column has the mean value subtracted, and then divided by the standard deviation of the whole column. This is a pretty ordinary step for most machine learning models and allows the whole network to learn faster (more on this below).

使用Scikit-Learn库中的StandardScaler()函数对数据集进行标准化。 对于数据框的每一列，该列中的每个值均减去平均值，然后除以整个列的标准差。 对于大多数机器学习模型来说，这是一个非常普通的步骤，它使整个网络学习得更快(下面有更多内容)。

Then, the dataframe is passed through a transformation function. For each column, we create a copy of each of the previous 7 days’ values (15 * 7 = 120 columns). The shape of the resulting dataframe is 7506 rows x 120 columns.

然后，数据帧通过转换函数传递。 对于每列，我们为前7天的每个值创建一个副本(15 * 7 = 120列)。 所得数据框的形状为7506行x 120列。

A lot of the code was inspired by this fantastic blog post.

许多代码是从这篇精彩的博客文章中获得启发的。

# load datasetvalues = df.values# ensure all data is floatvalues = values.astype('float32')# normalise each feature variable using Scikit-Learnscaler = StandardScaler()scaled = scaler.fit_transform(values)# save scaler for later usejoblib.dump(scaler, 'scaler.gz')# specify the number of lagged steps and featuresbackward_steps = 7n_features = df.shape[1]# convert series to supervised learningdef series_to_supervised(data, n_in=1, n_out=1, dropnan=True): n_vars = 1 if type(data) is list else data.shape[1] df = pd.DataFrame(data) cols, names = list(), list() # input sequence (t-n, ... t-1) for i in range(n_in, 0, -1): cols.append(df.shift(i)) names += [('var%d(t-%d)' % (j+1, i)) for j in range(n_vars)] # forecast sequence (t, t+1, ... t+n) for i in range(0, n_out): cols.append(df.shift(-i)) if i == 0: names += [('var%d(t)' % (j+1)) for j in range(n_vars)] else: names += [('var%d(t+%d)' % (j+1, i)) for j in range(n_vars)] # put it all together agg = pd.concat(cols, axis=1) agg.columns = names # drop rows with NaN values if dropnan: agg.dropna(inplace=True) return agg# frame as supervised learningreframed = series_to_supervised(scaled, backward_steps, 1)

Define training and test datasets

定义训练和测试数据集

We must split the prepared dataframe into training and test datasets to allow a fair evaluation of our results. The training dataset represents 80% of our values and we will use the remaining 20% for evaluation. As we are dealing with data ordered through time, it is a very bad idea to shuffle the dataset, so we keep it as is. Next, we reshape our training and test datasets into three dimensions for later use.

我们必须将准备好的数据框分为训练和测试数据集，以便对我们的结果进行公正的评估。 训练数据集代表我们价值的80％，我们将使用剩余的20％进行评估。 由于我们要处理按时间排序的数据，因此改组数据集是一个非常糟糕的主意，因此我们将其保持原样。 接下来，我们将训练和测试数据集重塑为三个维度，以供以后使用。

# split into train and test setsvalues = reframed.valuesthreshold = int(0.8 * len(reframed))train = values[:threshold, :]test = values[threshold:, :]# split into input and outputsn_obs = backward_steps * n_featurestrain_X, train_y = train[:, :n_obs], train[:, -n_features]test_X, test_y = test[:, :n_obs], test[:, -n_features]print(train_X.shape, len(train_X), train_y.shape)# reshape input to be 3D [samples, timesteps, features]train_X = train_X.reshape((train_X.shape[0], backward_steps, n_features))test_X = test_X.reshape((test_X.shape[0], backward_steps, n_features))print(train_X.shape, train_y.shape, test_X.shape, test_y.shape)

Fit model

拟合模型

Finally, we are able to fit our LSTM network. Thanks to the TensorFlow/Keras library, this only requires a few lines of code. I have chosen to fit 64 memory blocks in batch sizes of 72. I use the Adam optimisation algorithm, which is more efficient than the classical gradient descent procedure.

最后，我们能够适应我们的LSTM网络。 感谢TensorFlow / Keras库，这仅需要几行代码。 我选择将批量大小为72的64个存储块装入。我使用Adam优化算法，该算法比经典的梯度下降过程更有效。

# design networkmodel = tf.keras.models.Sequential()model.add(tf.keras.layers.LSTM(64, input_shape=(train_X.shape[1], train_X.shape[2])))model.add(tf.keras.layers.Dense(1))model.compile(loss='mae', optimizer='adam')# define early stopping parametercallback = tf.keras.callbacks.EarlyStopping(monitor='loss', patience=3)# fit networkhistory = model.fit(train_X, train_y, epochs=25, callbacks=[callback], batch_size=72, validation_data=(test_X, test_y), verbose=2, shuffle=False)# plot historyplt.figure(figsize=(12, 6))plt.plot(history.history['loss'], label='train')plt.plot(history.history['val_loss'], label='test')plt.ylabel('mean absolute error [Kaub, normalised]')plt.legend()plt.show()

模型结果 (Model results)

After the model is fit, we can launch the forecast and invert the scaling to obtain our final results. We can then calculate an error score for the model. Here, the model achieved a root mean squared error (RMSE) of 6.2 centimetres, which is good but can probably be improved on.

模型拟合后，我们可以启动预测并反转比例以获取最终结果。 然后，我们可以计算模型的错误评分。 在此，该模型实现了6.2厘米的均方根误差(RMSE)，这很好，但可能会得到改善。

# make a predictionyhat = model.predict(test_X)test_X = test_X.reshape((test_X.shape[0], backward_steps*n_features))# invert scaling for forecastinv_yhat = np.concatenate((yhat, test_X[:, -(n_features - 1):]), axis=1)inv_yhat = scaler.inverse_transform(inv_yhat)inv_yhat = inv_yhat[:,0]# invert scaling for actualtest_y = test_y.reshape((len(test_y), 1))inv_y = np.concatenate((test_y, test_X[:, -(n_features - 1):]), axis=1)inv_y = scaler.inverse_transform(inv_y)inv_y = inv_y[:,0]# calculate RMSErmse = np.sqrt(mean_squared_error(inv_y, inv_yhat))print('Test RMSE: %.3f' % rmse)

I have also checked the model’s performance during extreme weather events, such as the extensive floods recorded in Europe in May-June 2016 which caused more than Eur1 billion of damage in Bavaria alone. The model was generally able to track the rise in water levels, however during two particular peaks (one in mid-April and the other in early June) it gave low figures.

我还检查了该模型在极端天气事件下的性能，例如2016年5月至6月在欧洲记录的大范围洪水，仅在巴伐利亚州就造成了10亿欧元以上的损失。 该模型通常能够跟踪水位的上升，但是在两个特定的高峰期(一个在4月中旬，另一个在6月初)，它给出的数值很低。

The trained LSTM network also performed well during storm Axel, in May 2019, which caused a very rapid rise in water height of more than 1 metre on two consecutive days. Once again, however, it gave slightly lower estimates than actual figures when floods peaked.

经过训练的LSTM网络在2019年5月的Axel风暴中也表现良好，导致连续两天水位Swift上升超过1米。 但是，当洪灾达到顶峰时，它给出的估算值略低于实际数字。

In both cases, I suspect the low figures are due to one major Rhine tributary, the Moselle, which was left out of the analysis for lack of reliable meteorological data. The Moselle joins the Rhine at Koblenz, just north of Kaub. This is a possible area of improvement in future.

在这两种情况下，我都怀疑数字偏低是由于莱茵河的一个主要支流摩泽尔河，由于缺乏可靠的气象数据而被排除在分析之外。 摩泽尔河(Moselle)在考布(Kaub)以北的科布伦茨(Koblenz)加入莱茵河。 这是将来可能需要改进的领域。

模型表现 (Model performance)

Finding the best hyperparameters is one of the most important (and time-consuming) tasks in machine learning. I ran many different versions of the network in order to find the best possible setup:

在机器学习中，找到最佳的超参数是最重要(且耗时)的任务之一。 我运行了许多不同版本的网络，以便找到可能的最佳设置：

Pre-processing: For those familiar with Scikit-Learn, I found the StandardScaler() to be more suitable for the Rhine dataset than the MinMaxScaler(), which normalises all values to between 0 and 1. Normalising values with the MinMaxScaler() gave erratic performance. A quick look at the histograms plotted at the beginning of this blog post indicates that most of the Rhine input variables follow a Gaussian distribution, so this is logical.

预处理：对于熟悉Scikit-Learn的人，我发现StandardScaler()比MinMaxScaler()更适合莱茵数据集，后者将所有值归一化为0到1。使用MinMaxScaler()归一化值性能不稳定。 快速浏览一下此博客文章开始处绘制的直方图，表明大多数莱茵河输入变量都遵循高斯分布，因此这是合乎逻辑的。

Backward steps: I found 7 days to be a reasonable time frame for weather analysis and given the time required for water to flow down the Alps. There is no meaningful improvement in the model’s performance by increasing the steps beyond this figure.

后退步骤：我发现7天是进行天气分析的合理时间范围，并且考虑到水从阿尔卑斯山流下所需的时间。 通过增加超出此数字的步骤，模型的性能没有任何有意义的改善。

Neurons: Small networks of 16 memory cells perform less well, but larger networks of 32 neurons and beyond are fairly similar. I picked 64 memory cells for my model.

神经元：由16个记忆细胞组成的小型网络性能较差，但由32个及以上的神经元组成的较大网络则非常相似。 我为模型选择了64个存储单元。

Epochs: I have chosen to stop training once the algorithm has gone through the dataset 25 times. Beyond this threshold, performance begins to plateau and the model starts to overfit the training set.

纪元：算法选择经过25次数据集后，我便选择停止训练。 超过此阈值，性能开始趋于平稳，模型开始过度拟合训练集。

I have deployed a version of this model online. Its performance can be monitored here.

我已经在线部署了该模型的一个版本。 它的性能可以在这里监控。

I’ll try to improve this model in future by adding/removing variables and tweaking hyperparameters further. If you have any thoughts, don’t hesitate to get in touch with me.

将来，我将尝试通过添加/删除变量并进一步调整超参数来改进此模型。 如果您有任何想法，请随时与我联系。

翻译自: https://towardsdatascience.com/using-machine-learning-to-predict-rhine-water-levels-44afce697074