【短临预报系列第一篇】Convolutional LSTM Network: A Machine LearningApproach for Precipitation Nowcasting_convolutional lstm network: a machine learning app

前言：

这是我研究短临预报看的第一篇文章，施行健博士所写，他提出的ConvLSTM模型从2015年到现在还有人使用，引用量达到了5000+，可以说是开山之作。

文章地址：https://arxiv.org/pdf/1506.04214.pdf

摘要：

降水预报的目标是预测一个地区在较短时间内的降雨强度。以前很少有研究从机器学习的角度来研究这一关键和具有挑战性的天气预报问题。在本文中，我们将降水临近预报定义为一个时空序列预测问题，其中输入和预测目标都是时空序列。通过扩展全连通LSTM(FC-LSTM)，使其在输入到状态和状态到状态的转换中都具有卷积结构，我们提出了卷积LSTM(ConvLSTM)，并利用它建立了降水临近预报问题的端到端可训练模型。实验表明，我们的ConvLSTM网络能够更好地捕捉时空相关性，且始终优于FC-LSTM和ROVER算法。

介绍：

首先，介绍了降水预报的意义，对于某个地区的精确预测可以有效避免因降雨引发的泥石流，洪水等自然灾害，同时还可以做好必要的防范措施，降水临近预报问题相当具有挑战性，已成为气象界的热门研究课题。

现代的降水预报方法大致分为两类，数值天气预报（NWP）和雷达外推方法。数值天气预报（NWP）是指根据大气实况，在一定的初值和边界条件下，通过大型计算机作数值计算，求解描写天气演变过程的流体力学和热力学的方程组，预测未来一定时间段的大气运动状态和天气现象的方法。雷达回波外推主通过模型训练来预测云层的移动，可以充分的利用雷达回波图历史数据，有效的捕获时空相关性。光流法以平移为主，无法模拟云团的变化。

这些技术问题可以通过从机器学习的角度来看这个问题来解决。本质上，降水临近预报是一个以输入过去雷达地图的序列，输出未来雷达地图的固定数（通常大于1）的序列的时空序列预测问题。然而，这种学习问题，无论它们的确切应用如何，首先都不重要，因为时空序列的高维性，特别是当必须进行多步预测时，除非预测模型能很好地捕捉到数据的时空结构。此外，由于大气的混沌特性，建立一个有效的雷达回波数据预测模型更具有挑战性。

在本文中，我们提出了一种新的卷积LSTM(ConvLSTM)网络用于降水临近预报。我们将降水临近预报表示为一个时空序列预测问题，可以在[23]中提出的一般序列到序列学习框架下得到解决。为了很好地模拟时空关系，我们将FC-LSTM的思想扩展到ConvLSTM，它在输入到状态和状态到状态的转换中都具有卷积结构。通过叠加多个ConvLSTM层，形成编码-预测结构，可以建立降水临近预报的端到端可训练模型。为了评估，我们创建了一个新的真实雷达回波数据集，可以促进进一步的研究，特别是设计问题的机器学习算法的研究。当在一个合成的Moving-MNIST数据集[21]和雷达回波数据集上进行评估时，我们的ConvLSTM模型始终优于FC-LSTM和ROVER算法。

Preliminaries：

降水临近预报的目标是利用先前观测到的雷达回波序列来预测当地地区（如香港、纽约或东京）的未来雷达地图的固定长度。在实际应用中，雷达地图通常每6-10分钟从天气雷达上获取一次，并在接下来的1-6小时内进行临近预报，即预测前方的6-60帧。降水预报实际上是时空预测问题。

 

 

 每个正方形小格点代表11公里可以看做一个地区。p值为降雨量，也可以为雷达反射率，它和降雨量之间有关系可以转换。这里肯定有朋友要问为什么要有个P，这个P的维度是干嘛的？对于雷达回波图来说每一个图片的对应像素点通过一定的变换其实是可以直接得到这个点的一些物理量的，所以这个P实际作者想说的就是这个物理量的维度。

 对于降水临近预报，每个时间戳的观测是二维雷达回波图。如果我们将地图划分为平铺的不重叠的斑块，并将斑块内的像素视为其测量值，临近预测问题自然成为一个时空序列预测问题。我们注意到，我们的时空序列预测问题不同于一步时间序列预测问题，因为我们的问题的预测目标是一个同时包含空间和时间结构的序列。

LSTM Model

LSTM模型，是循环神经网络的一种变体，可以很有效的解决简单循环神经网络的梯度爆炸或消失问题。

 

遗忘门ft，控制上一时刻的内部状态ct-1需要遗忘多少信息。

输入门it，控制当前时刻的候选状态有多少信息需要保存。

（ct 非线性函数得到的候选状态）

输出门ot，控制当前时刻的内部状态ct有多少信息需要输出给外部状态ht。

新的内部状态LSTM网络引入一个新的内部状态Ct专门进行线性的循环信息传递，同时（非线性）输出信息给隐藏层的外部状态ht。LSTM网络中的门是一种软门，取值在（0，1）之间，表示以一定的比例运行信息通过。

Convolutional LSTM

FC-LSTM在处理时空数据方面的主要缺点是它在没有进行空间信息编码的输入到状态和状态到状态的转换中使用了完全连接,缺乏空间信息的提取。

为了克服这个问题，我们设计的一个显著特点是，所有的输入都是，单元输

隐藏状态，和ConvLSTM的门都是3D张量，最后两个维度是行和列，ConvLSTM通过其本地邻居的输入和过去的状态来确定网格中某个单元格的未来状态。这可以很容易地通过在状态到状态和输入到状态的转换中使用卷积算子来实现

 

 从图中可以看出输入和隐藏状态都是卷积操作，有效提取空间信息。

为了确保状态具有与输入相同的行数和相同数量的列，在应用卷积操作之前需要padding。在这里，隐藏状态的填充在可以将边界点视为使用外部世界的状态进行计算。通常，在第一个输入到来之前，我们将LSTM的所有状态初始化为零，这对应于未来的“完全无知”。类似地，如果我们对隐藏状态执行零填充(本文中使用的)，我们实际上是将外部世界的状态设置为零，并且假定没有关于outsid的先验知识。通过在状态上填充，我们可以对边界点进行不同的处理，这在许多情况下是有帮助的。例如，想象一下我们正在观察的系统是一个被墙包围的移动球。虽然我们看不到这些墙，但我们可以通过一次又一次地发现球而推断它们的存在，当边界点与内点具有相同的状态转移动力学时，很难做到这一点。

预测-编码结构

与FC-LSTM一样，ConvLSTM也可以作为更复杂结构的构建块。对于我们的时空序列预测问题，我们使用图3所示的结构。两个网络、一个编码网络和一个预测网络。类似于中，预测网络的初始状态和单元输出是从编码网络的最后一个状态复制出来的。这两个网络都是通过堆叠几个ConvLSTM层来形成的。由于我们的预测目标具有与输入相同的维数，所以我们将预测网络中的所有状态连接起来，并将它们送入 1 × 1 卷积层，以生成最终预测。

 

我们可以使用与[23]类似的观点来解释此结构。编码LSTM将整个输入序列压缩成一个隐藏状态张量，预测LSTM将这个隐藏状态展开，得到最终的预测结果:

这种结构也类似于[21]中的LSTM未来预测器模型，但我们的输入和输出单元都是保持所有空间信息的三维张量。由于网络具有多个堆叠的ConvLSTM层，它具有很强的代表性，适用于我们研究的降水预报等复杂动力系统的预报。

 实验

我们首先将我们的ConvLSTM网络与一个合成MovingMNIST数据集中的FC-LSTM网络进行比较，以获得对我们模型的行为的一些基本理解。我们用不同的号码运行我们的模型根据[21]中的层和内核大小，并研究一些“域外域”情况。为了验证我们的模型在更具有挑战性的降水预报问题上的有效性，我们建立了一个新的雷达回波数据集，并将我们的模型与基于几种常用降水预报指标的最先进的ROVER者算法进行了比较。在这两个数据集上进行的实验结果得出以下结论：

• ConvLSTM在处理时空相关性方面优于FC-LSTM.
• 使得状态到状态卷积核的大小大于1对于捕获时空运动模式是必要的。
• 更深层次的模型可以较少的参数产生更好的结果。
• ConvLSTM在降水预报方面的性能优于ROVER。

我们的模型实现是在Theano [3，1]的帮助下用Python实现的。我们用一台NVIDIA K20 GPU在一台计算机上运行所有的实验。此外，附录中还包括更具说明性的“gif” 例子

Moving-MNIST数据集

对于这个合成数据集，我们使用类似于[21]中描述的生成过程。数据集中的所有数据实例长20帧(输入10帧，预测10帧)。并包含在64×64补丁内弹跳的两个手写数字。移动数字是从MNIST数据集中500位的子集中随机选择的。在[3，5]中，开始位置和速度方向都是随机选择的，速度振幅是随机选择的。该生成过程重复了15000次，得到了包含10000条训练序列、2000条验证序列和3000条测试序列的数据集。我们通过使用反向传播时间(BPTT)[2]最小化交叉熵损失来训练所有LSTM模型，并且RMSProp[24]，学习速率为10×3，衰减率为0.9。此外，我们在验证集上执行早期停止。

 

尽管生成过程很简单，但结果在所得到的数据集中存在强的非线性，因为移动的数字可能表现出复杂的外观，并且在它们的运动期间将阻塞和反弹。在不了解系统内部动力学的情况下，模型很难对测试集进行准确的预测。

对于FC-LSTM网络，我们使用了与[21]中具有两个2048节点的LSTM层的无条件未来预测器模型相同的结构。对于我们的ConvLSTM网络，我们将补丁大小（patch size）设置为4×4，这样每个64×64帧由16×16×16张量表示。我们用不同的层数测试我们的模型的三种变体。1层网络包含一个包含256个隐藏状态的ConvLSTM层。2层网络有两个ConvLSTM层，每个层有128个隐藏状态，而3层网络在三个ConvLSTM层中分别有128、64和64隐藏状态。所有输入到状态和状态到状态的内核大小为5×5。我们的实验表明，ConvLSTM 网络的性能始终优于 FC-LSTM 网络. 此外，此外，更深层次的模型可以提供更好的结果，尽管两层和三层网络之间的改进并不显著。此外，我们还尝试了其他网络配置，将2层和3层网络的状态到状态和输入到状态的内核分别改变为1×1和9×9。虽然新的2层网络的参数与原网络接近，但由于仅仅使用1×1的state-to-state的转变难以捕捉到时空运动模式，结果更加糟糕。同时，新的三层网络比新的两层网络性能更好，因为较高层可以看到更宽的输入范围。然而，它的性能不如具有较大的状态到状态内核大小的网络。这提供了更大的状态核更适合捕捉时空相关性的证据。事实上，对于1×1内核，状态的接受场不会随着时间的推移而增长。但是对于较大的核，后期的状态有更大的接受场，并且与更大的范围有关。所有的信息。 测试集上各算法的平均交叉熵损耗(每个序列的交叉熵损耗) 如表1 所示。 我们需要指出的是，我们的实验设置与[21]不同，在[21]中，假定有无限多的训练数据可用。选择当前的脱机设置是为了了解在数据不足的情况下，不同的模型是如何执行的。3层ConvLSTM和FC-LSTM在线环境下的比较见附录。

接下来，我们在一些“领域外”的输入上测试我们的模型。我们再生成3，000个移动数字序列，使用从500 MNIST的不同子集随机抽取的与训练集不重叠的数字。由于模型从来没有见过任何三位数的系统，这样的“域外”运行是对模型泛化能力的很好的测试[21]。平均该数据集上三层模型的交叉熵误差为6379.42.通过对部分预测结果的观察，发现该模型能较好地分离重叠位，预测t值。他的整体运动，虽然预测的数字是相当模糊的。一个“领域外”的预测例子如图10所示。

雷达回波数据集

本文所使用的雷达回波数据集是2011年至2013年在香港收集的3年天气雷达强度的子集。由于不是每天都下雨，我们的预报目标是降水，我们选择前97个雨天形成我们的数据集。对于预处理，我们首先通过设置 将强度值Z 转换为灰度像素P ，并在中心 330 × 330 区域裁剪雷达图。之后，我们将带有半径10的磁盘过滤器应用到100x100，并调整雷达图的大小。为了减少测量仪器带来的噪声，我们进一步去除了一些噪声区域的像素值，这些区域是通过将K-均值聚类方法应用于每月像素平均值来确定的。天气雷达数据每6分钟记录一次，因此每天有240帧。为了获得不相交的训练、测试和验证子集，我们将每个日序列划分为40个不重叠的帧块，并随机分配4个块进行训练，1块用于测试并且1块用于验证。数据实例是使用一个20帧宽的滑动窗口从这些块中分割出来的。因此，我们的雷达回波数据集包含8148个训练序列，2037个测试序列和2037个验证序列，所有序列都有20帧长(输入5帧，预测15帧)。虽然从同一天分割的训练和测试实例可能有一些依赖关系，但这种拆分策略仍然是合理的，因为在现实生活中，我们可以访问所有以前的数据，包括同一天的数据，让我们可以在线微调模型。这种数据分割可能被看作是该应用程序的实际“微调”设置的近似值。

我们将补丁大小(patch size)设为2，并训练一个包含64个隐藏状态和3×3个内核的2层ConvLSTM网络。对于ROVER算法，我们在验证集上调整光流估计器的参数，并使用最佳参数(如附录所示)报告测试结果。此外，我们还尝试了三种不同的Rover初始化方案：ROVER 1计算最后两个观测帧的光流，然后进行半拉格朗日平流；ROVER 2以最后两个流场的平均值初始化速度；ROVER 3给出最后三个流场的加权平均值(权重分别为0.7、0.2和0.1)初始化。此外，我们还训练了一个具有两个 2000 节点 LS TM层的 FC-LSTM 网络。无论是 ConvLSTM 网络还是 FC-LSTM 网络都优化了 15 个预测的交叉熵误差。

我们使用几种常用的降水预报指标，即雨量均方误差(雨量均方误差)、关键成功指数(CSI)、虚警率(Far)、检测概率(POD)以及相关性来评价这些方法。降雨-MSE度量被定义为预测降雨与地面真相之间的平均平方误差。由于我们的预测是在像素级完成的，所以我们将它们投射回雷达回波强度，并使用z-r关系计算网格的每个单元的降雨量[15]：Z = 10 log a + 10b log R，其中Z是以分贝为单位的雷达回波强度，R是以毫米/小时为单位的雨量，a，b是两个常数，a=118.239，b=1.5241。CSI、FAR和POD是类似于机器学习研究人员常用的精确和回忆的技能分数。我们使用0.5mm/h降雨率(表示是否下雨)将预测和地面真相转换为0/1矩阵，并计算命中(预测=1，真相=1)，失败(预测)ON=0，true=1)和假警报(预测=1，true=0)。三个技能分数定义为：，，。预测框架P与地面真框架T的相关性定义为：其中，。

关于为什么选择97个雨天，如果无雨样本太多训练出来效果会出现全黑的图片，我们至少要保证每张图大于百分20的降雨量。和这张图差不多，不然预测效果会不理想。

 

 

 

ROVER2与ConvLSTM的对比预测图。
最后作者总结了一下这个实验室的结论：

所有结果如表2和图5所示。我们可以发现，FC-LSTM网络的性能并不好地完成这个任务，这主要是由于雷达地图上的强空间相关性，即云的运动在一个局部区域是高度一致的。全连接的结构有太多的冗余连接，这使得优化不太可能捕获这些局部一致性。此外，还可以看出，ConvLSTM的性能优于基于光流的漫游者算法，这主要有两个原因。首先，ConvLSTM能够很好地处理边界条件。在现实生活中，有很多情况下，边界处突然聚集的云，这表明有些云是从外部来自的。如果ConvLSTM网络在训练过程中看到了类似的模式，它可以发现编码网络中这种类型的突然变化，并在预测网络中给出合理的预测。然而，这一点很难用基于光流和基于半拉格朗日平流的方法来实现。另一个原因是，ConvLSTM被端到端训练为该任务和数据中的一些复杂的时空模式

结论与展望

在本文中，我们已经成功地将机器学习方法，特别是深度学习，应用到具有挑战性的降水临近预报问题上，而迄今为止还没有受益于复杂的机器学习技术。我们将降水临近预报定义为一个时空序列预测问题，并提出了一个新的LSTM扩展，称为ConvLSTM来解决这一问题。ConvLSTM层不仅保留了FC-LSTM的优势，而且由于其固有的卷积结构，也适用于时空数据。通过将ConvLSTM纳入编码预测结构，我们建立了一个端到端可训练的预报模型。对于未来的工作，我们将研究如何将ConvLSTM应用于基于视频的动作识别。一种想法是在卷积神经网络生成的空间特征图上添加ConvLSTM，并使用ConvLSTM的隐藏状态进行最终的分类。

个人认为本文最大的创新是提出了ConvLSTM结构并用于降水预测。

代码

import torch
import torch.nn as nn
from torch.autograd import Variable
 
 
class ConvLSTMCell(nn.Module):
    def __init__(self, input_channels, hidden_channels, kernel_size):
        super(ConvLSTMCell, self).__init__()
 
        assert hidden_channels % 2 == 0
 
        self.input_channels = input_channels
        self.hidden_channels = hidden_channels
        self.kernel_size = kernel_size
        self.num_features = 4
 
        self.padding = int((kernel_size - 1) / 2)
 
        self.Wxi = nn.Conv2d(self.input_channels, self.hidden_channels, self.kernel_size, 1, self.padding, bias=True)
        self.Whi = nn.Conv2d(self.hidden_channels, self.hidden_channels, self.kernel_size, 1, self.padding, bias=False)
        self.Wxf = nn.Conv2d(self.input_channels, self.hidden_channels, self.kernel_size, 1, self.padding, bias=True)
        self.Whf = nn.Conv2d(self.hidden_channels, self.hidden_channels, self.kernel_size, 1, self.padding, bias=False)
        self.Wxc = nn.Conv2d(self.input_channels, self.hidden_channels, self.kernel_size, 1, self.padding, bias=True)
        self.Whc = nn.Conv2d(self.hidden_channels, self.hidden_channels, self.kernel_size, 1, self.padding, bias=False)
        self.Wxo = nn.Conv2d(self.input_channels, self.hidden_channels, self.kernel_size, 1, self.padding, bias=True)
        self.Who = nn.Conv2d(self.hidden_channels, self.hidden_channels, self.kernel_size, 1, self.padding, bias=False)
 
        self.Wci = None
        self.Wcf = None
        self.Wco = None
 
    def forward(self, x, h, c):
        ci = torch.sigmoid(self.Wxi(x) + self.Whi(h) + c * self.Wci)
        cf = torch.sigmoid(self.Wxf(x) + self.Whf(h) + c * self.Wcf)
        cc = cf * c + ci * torch.tanh(self.Wxc(x) + self.Whc(h))
        co = torch.sigmoid(self.Wxo(x) + self.Who(h) + cc * self.Wco)
        ch = co * torch.tanh(cc)
        return ch, cc
 
    def init_hidden(self, batch_size, hidden, shape):
        if self.Wci is None:
            self.Wci = nn.Parameter(torch.zeros(1, hidden, shape[0], shape[1])).cuda()
            self.Wcf = nn.Parameter(torch.zeros(1, hidden, shape[0], shape[1])).cuda()
            self.Wco = nn.Parameter(torch.zeros(1, hidden, shape[0], shape[1])).cuda()
        else:
            assert shape[0] == self.Wci.size()[2], 'Input Height Mismatched!'
            assert shape[1] == self.Wci.size()[3], 'Input Width Mismatched!'
        return (Variable(torch.zeros(batch_size, hidden, shape[0], shape[1])).cuda(),
                Variable(torch.zeros(batch_size, hidden, shape[0], shape[1])).cuda())
 
 
class ConvLSTM(nn.Module):
    # input_channels corresponds to the first input feature map
    # hidden state is a list of succeeding lstm layers.
    def __init__(self, input_channels, hidden_channels, kernel_size, step=1, effective_step=[1]):
        super(ConvLSTM, self).__init__()
        self.input_channels = [input_channels] + hidden_channels
        self.hidden_channels = hidden_channels
        self.kernel_size = kernel_size
        self.num_layers = len(hidden_channels)
        self.step = step
        self.effective_step = effective_step
        self._all_layers = []
        for i in range(self.num_layers):
            name = 'cell{}'.format(i)
            cell = ConvLSTMCell(self.input_channels[i], self.hidden_channels[i], self.kernel_size)
            setattr(self, name, cell)
            self._all_layers.append(cell)
 
    def forward(self, input):
        internal_state = []
        outputs = []
        for step in range(self.step):
            x = input
            for i in range(self.num_layers):
                # all cells are initialized in the first step
                name = 'cell{}'.format(i)
                if step == 0:
                    bsize, _, height, width = x.size()
                    (h, c) = getattr(self, name).init_hidden(batch_size=bsize, hidden=self.hidden_channels[i],
                                                             shape=(height, width))
                    internal_state.append((h, c))
 
                # do forward
                (h, c) = internal_state[i]
                x, new_c = getattr(self, name)(x, h, c)
                internal_state[i] = (x, new_c)
            # only record effective steps
            if step in self.effective_step:
                outputs.append(x)
 
        return outputs, (x, new_c)
 
 
if __name__ == '__main__':
    # gradient check
    convlstm = ConvLSTM(input_channels=512, hidden_channels=[128, 64, 64, 32, 32], kernel_size=3, step=5,
                        effective_step=[4]).cuda()
    loss_fn = torch.nn.MSELoss()
 
    input = Variable(torch.randn(1, 512, 64, 32)).cuda()
    target = Variable(torch.randn(1, 32, 64, 32)).double().cuda()
 
    output = convlstm(input)
    output = output[0][0].double()
    res = torch.autograd.gradcheck(loss_fn, (output, target), eps=1e-6, raise_exception=True)
    print(res)

convlstm代码：GitHub - automan000/Convolutional_LSTM_PyTorch: Multi-layer convolutional LSTM with Pytorch