使用LSTM模型进行时序数据预测

1. LSTM模型简介

        LSTM，全程长短期记忆网络（Long Short-Term Memory），是在循环神经网络基础上发展起来的，主要是为了解决长序列训练过程中的梯度消失和梯度爆炸问题。LSTM神经网络能够处理长时间序列数据，具有长时记忆能力，因此成为解决时间序列数据预测问题的有效方法。其他关于LSTM的发展过程和相关原理就不多赘述了，本文重点是结合实例介绍如何使用LSTM进行预测。

2. 研究目标

        使用临近若干时刻的气象观测数据，预测未来3个时刻的能见度。

3. 准备工作

        3.1 环境准备

        python 3.x环境，必要库：pytorch、numpy、pandas、scikit-learn

        3.2 数据准备

        全国气象站逐小时观测数据，要素包括：温度、露点温度、相对湿度、饱和水汽压差、能见度。每个时刻的观测数据为一个csv文件，每个文件中的行代表不同的站点观测值，列代表不同的要素属性，如下所示：

        数值为999999表示缺测，这是我们要清洗掉的数据。

4. 数据预处理

        4.1 读取数据

        首先从文件中读取数据：

import pandas as pd
import os
import numpy as np
 
indir = r'./filepath'    #保存所有数据文件的文件夹
file_list = os.listdir(indir)    #获得所有数据文件的文件名
feature_vars = ['TEM', 'DPT', 'RHU', 'VAP', 'VIS']
label_vars = 'VIS'
 
dataset = []
for f in range(len(file_list)):
    df = pd.read_csv(os.path.join(indir, f), na_values=[999999, 999107], usecols=['Station_Id_C', 'Lat', 'Lon', 'TEM', 'DPT', 'RHU', 'VAP', 'VIS'])
    df = df.dropna(axis=0, how='any').query('Lat>17 & Lat<29 & Lon>108 & Lon<123')    #剔除所有存在空值的行，并提取目标空间范围内的数据
    dataset.append(df.loc[:, feature_vars].values)
dataset = np.array(dataset, dtype=np.float32)    #(times, stations, features)
 

        上步将所有数据都读进了dataset三维数组中，三个维度分别是时间、站点、要素。dropna函数用于剔除所有存在空值的行，query函数用于筛选位于指定范围内的站点。接下来需要根据要预测的时刻（pred_step）和用于预测的时刻（time_step）将数据构造为样本集。

        4.2 创建样本集

        （1）什么是样本？

        准确理解样本的含义对于接下来的工作非常重要！

        样本就是一组数据，它包括两部分内容：特征数据和标签数据。标签数据很简单，就是你的预测目标，而特征数据则是用来计算标签数据的数据。举个例子：

        假设我想要用当前时刻的温度和气压去预测风速，那么温度和气压组成的数据就是特征数据，风速则是标签数据，特征数据和标签数据又组成一个样本。如果我有100个时刻的温度、气压、风速数据，那么就能形成100个样本，这便是样本集。如果只有单站点数据，那么样本数量把往往和时间长度有关；如果是多站点数据，那么时间长度与站点数量都会影响最终的样本总数。

        （2）构造样本集

        上面对于样本的例子是最简单的一种，现在来个复杂点的：我们需要用最近5个时刻的温度和气压数据，去预测未来3个时刻的风速数据。

        显然，现在我们有了三个“时间”，分别是用于预测的时间（time_step）、预测的目标时间（pred_step）、样本数据的时间（times）。为了便于理解，不妨假设原始观测资料共有times=100个连续时刻：

如上图所示，在时间序列上从左向右滑动，每次滑动都能生成一个样本，显然时序总长为100，time_step=5，pred_step=3时，总共能够生成93个样本，每个样本包含：

特征数据：time_step*features，即5个时刻*2个要素

标签数据：pred_step*labels，即3个时刻*1个要素

经过这一步处理后，原来的时间长度缩短了（times - time_step - pred_step + 1)，新时间维的长度就是这个站点所能够生成的样本数，如果总共有N个站点，那么样本总数就是N*(times - time_step - pred_step + 1)。代码如下：

time_step = 5
pred_step = 3
samples = []
X = []
Y = []
 
for i in range(dataset.shape[0] - time_step - pred_step + 1):
    X.append(dataset[i:i+time_step, :, :-1])    #(time_step, stations, features)
    Y.append(dataset[i+time_step:i+time_step+pred_step, :, -1])    #(pred_step, stations)
X = np.array(X).transpose(0, 2, 1, 3)    #(new_times, stations, time_step, features)
Y = np.array(Y).transpose(0, 2, 1)    #(new_times, stations, pred_step)

        注意，这里由于标签要素的值恰好在dataset的最后一列，所以直接切片，不然的话得根据位置调整。

        到了这一步，就完全可以把新的时间维度与站点维度合并成样本维度了：

samples_x = X.reshape(-1, time_step, features)    #形状为(samples, time_step, features)
samples_y = Y.reshape(-1, pred_step)    #形状为(samples, pred_step)

        样本集创建完成，样本总数为（times - time_step - pred_step + 1) * stations。

        4.3 拆分样本集

        将样本集拆分成训练集、验证集、测试集三部分：

from sklearn.model_selection import train_test_split
 
train_x, tmp_x, train_y, tmp_y = train_test_split(data_x, data_y, train_size=0.6, random_state=42, shuffle=True)
val_x, test_x, val_y, test_y = train_test_split(tmp_x, tmp_y, train_size=0.5, random_state=42, shuffle=True)

        4.4 样本归一化

        不同要素的数值范围可能又很大差别，如果直接使用会导致它们对模型预测的影响度不同，因此需要先进行归一化。

from sklearn.preprocessing import MinMaxScaler
 
scale_x = MinMaxScaler()
train_x = scale_x.fit_transform(train_x.reshape(-1, len(feature_vars)).reshape(-1, time_step, len(feature_vars))
val_x = scale_x.transform(val_x.reshape(-1, len(feature_vars))).reshape(-1, time_step, len(feature_vars))
test_x = scale_x.transform(test_x.reshape(-1, len(feature_vars))).reshape(-1, time_step, len(feature_vars))
 
scale_y = MinMaxScaler()
train_y = scale_y.fit_transform(train_y)
val_y = scale_y.transform(val_y)
test_y = scale_y.transform(test_y)

        之所以不先归一化再拆分样本集，是为了避免归一化的过程中训练样本获取到了测试集的数据信息，影响测试效果。这里记得保留scale_y这个归一化器，后面在模型预测完成后可以用它对预测值逆归一化，就能使预测数值重新变得可解释。

5. 模型训练

        5.1 定义LSTM模型

        首先需要根据我们的数据和任务需求定义一个lstm模型，pytorch已经提供了很便捷的方法，我们直接调用即可：

import torch
import torch.nn as nn
import torch.optim as optim
 
class LSTMModel(nn.Module):
    '''
    long-short term model
    '''
    def __init__(self, input_size, hidden_size, output_size, num_layers=1, dropout_rate=0.2):
        super(LSTMModel, self).__init__()
        self.input_size = input_size
        self.hidden_size = hidden_size
        self.output_size = output_size
        self.num_layers = num_layers
        self.dropout_rate = dropout_rate
        self.lstm = nn.LSTM(input_size, hidden_size, num_layers, batch_first=True)
        self.dropout = nn.Dropout(dropout_rate)
        self.fc1 = nn.Linear(hidden_size, output_size)
        self.relu = nn.ReLU()
 
    def forward(self, x):
        h0 = torch.zeros(self.num_layers, x.size(0), self.hidden_size).to(x.device)
        c0 = torch.zeros(self.num_layers, x.size(0), self.hidden_size).to(x.device)
        out, _ = self.lstm(x, (h0, c0)) # LSTM层
        out = self.dropout(out[:, -1, :])    #非必需
        out = self.fc1(out) # 全连接层
        out = self.relu(out)    #非必需
        return out

上面的模型包含两个函数块，__init__函数为LSTMModel类的初始化函数，对模型的各项参数进行了定义。forward函数定义了数据在模型中的流动、传递和操作，就像流水线加工一样，先经过第一道函数处理，再传给第二个，以此类推......

数据在模型中的流动过程为：

（1）首先特征数据x被输入到模型中，先通过lstm层，这一层接受形如[batch_size, seq_len, input_size]的三维张量。其中batch_size为每个批次输入的样本数（后面会再解释相关概念）；seq_len表示样本的序列长度，对应于time_step；input_size表示样本的特征数量，对应于len(feature_vars）。

lstm默认的输入维度为[seq_len, batch_size, input_size]，但我们可以通过令batch_first=True将前两维交换位置。

值得注意的是，如果数据的time_step=1（即每个样本只有1个时刻的特征要素数据），或者len(feature_vars)=0（即每个样本只有1个特征要素），那也必须保持三维，只是将对应的维度长度设置为1，即输入的张量x变成[batch_size, 1, input_size]或[batch_size, seq_len, 1]。

lstm层输出为out=[batch_size, seq_len, hidden_size]，显然x的最后一个维度被变成了hidden_size，其他维度长度没变化；但是再seq_len维度上，out储存了模型每个时间步的预测值，其中只有最后一个时间步，即out[:, -1, :]才是我们需要的最终预测值。

（2）lstm的输出被传入dropout层，这一步主要是为了随机丢弃一些神经元，避免模型过拟合，非必需。

（3）上步输出结果传入fc1全连接层，得到[batch_size, output_size]，这里的output_size根据任务需求由自己定义。比如我需要预测1个要素的未来3个时刻数值，则定义output_size=3。

（4）上步结果通过relu函数，去除负值。这一步也非必需。

        5.2 模型训练

import torch.optim as optim
from torch.utils.data import DataLoader, TensorDataset
import matplotlib.pyplot as plt
 
def set_loader(x, y, batch_size):
    '''
    输入要素数据和标签数据，构造数据加载器
    参数:
    -----
    x: 要素数据。数组
    y: 标签数据。数组
    '''
    
    tensor_x = torch.from_numpy(x)
    tensor_y = torch.from_numpy(y)
    loader = DataLoader(TensorDataset(tensor_x, tensor_y), batch_size=batch_size, shuffle=True)
    
    return loader
 
def lstm_train(model, epochs, train_loader, val_loader, learning_rate=0.01, plot_loss=False):
    '''
    LSTM模型训练
    '''
    # 定义损失函数和优化器
    loss_function = nn.MSELoss()
    optimizer = optim.Adam(model.parameters(), lr=learning_rate)
 
    train_losses = []
    val_losses = []
 
    for epoch in range(epochs):
        model.train()
        train_loss = 0.0
        for x_batch, y_batch in train_loader:
            optimizer.zero_grad()   #每一次循环前，梯度清零
            outputs = model(x_batch)
            loss = loss_function(outputs, y_batch)
            loss.backward()         #反向传播
            optimizer.step()        #梯度下降
            train_loss += loss.item() * x_batch.size(0)
        train_loss /= len(train_loader.dataset)
        train_losses.append(train_loss)
 
        model.eval()
        val_loss = 0.0
        with torch.no_grad():
            for x, y in val_loader:
                outputs = model(x)
                loss = loss_function(outputs, y)
                val_loss += loss.item() * x.size(0)
        val_loss /= len(val_loader.dataset)
        val_losses.append(val_loss)
 
        print(f'Epoch [{epoch+1}/{epochs}], Train Loss: {train_loss:.4f}, Val Loss: {val_loss:.4f}')
 
    if plot_loss:
        #查看损失函数曲线
        fig = plt.figure(figsize=[8,6])
        ax = fig.add_subplot(111)
        ax.plot(train_losses, 'b', label='train_losses')
        ax.plot(val_losses, 'r', label='val_losses')
        ax.legend()
        ax.set_title(f'Epochs:{epochs}  learning_rate:{learning_rate}')
        plt.show()
    
    return model
 
hidden_size = 16
input_size = len(feature_vars)
output_size = pred_time
batch_size = 16
epochs = 500
 
train_loader = set_loader(train_x, train_y, batch_size)
val_loader = set_loader(val_x, val_y, batch_size)
 
model = LSTMModel(input_size, hidden_size, output_size, num_layers=1, dropout_rate=0.3)
 
model = lstm_train(model, epochs, train_loader, val_loader, learning_rate=0.001, plot_loss=False)

set_loader函数用于根据样本数据集创建数据加载器，主要做了两件事：（1）把数组转成torch.tensor格式，这是必要的；（2）根据设置好的batch_size，把样本分批次。

batch_size的意义在于每次向模型投入不同数量的样本用于训练，例如样本总数100，batch_size=10，则每次训练都要分10次投放样本。一般来说，batch_size习惯设置为2的幂次数，当然非必须，最小可以设置为1，即每个只投放1个样本。batch_size越大，模型训练越快，但是收敛也会比较慢。

lstm_train函数用于训练模型，需要输入的参数包括：

model：实例化的LSTMModel。

epochs：模型训练次数。

train_loader：训练样本集的数据加载器。

val_loader：验证样本集的数据加载器。

learning_rate：学习率，数值范围为0~1。又是一个需要反复尝试的参数。学习率过大，模型会无法收敛；学习率太小，模型则容易过拟合。

plot_loss：bool型。若为True，则模型训练完成后会弹出损失函数曲线图，用于辅助调参。

再具体看lstm_train函数内部：

loss_function = nn.MSELoss()为模型损失函数，用于衡量预测值和真实值的差距。

optimizer = optim.Adam(model.parameters(), lr=learning_rate)定义了一个优化器。

这两个都可以根据实际需求修改。

for epoch in range(epochs)开始进行多次模型训练：

for x_batch, y_batch in train_loader:

            optimizer.zero_grad()   #每一次循环前，梯度清零

            outputs = model(x_batch)

            loss = loss_function(outputs, y_batch)

            loss.backward()         #反向传播

            optimizer.step()        #梯度下降

            train_loss += loss.item() * x_batch.size(0)

梯度清零、反向传播、梯度下降是3个固定步骤，train_loss则是用于保存模型再训练集上预测的损失函数值，便于诊断模型性能。

后面在验证集上对模型执行了同样的操作，并记录了验证集损失函数值。

两个损失函数的变化可以反应模型的学习效果与泛化能力：

理想情况下，两条曲线应该是先快速下降，然后慢慢稳定下来。

如果曲线一直在下降，那么可能是训练次数不够，需要调高epochs。

如果train_loss下降/稳定，而val_loss却稳定/上升，那么说明模型过拟合了，可以简化一下网络结构或者增加dropout层。

        5.3 模型预测

        上步，我们进行了模型训练和验证，并完成了调优，接下来就可以在测试集上进行预测了。

model.eval()
with torch.no_grad():
    out = model(torch.from_numpy(test_x))    #模型预测
pred = scale_y.inverse_transform(out)    #结果逆归一化

        原数据的变化范围太大，预测结果在极大极小的表现较差，但是基本能够预测出变化趋势。