深度学习神经网络实践：基于LSTM的电离层预报

1 引言

在本科的毕业设计中，我的题目就是基于深度学习的电离层预报，在初期的测试了中，试了挺多的网络，也试了挺多的预报方法和输出的方式，多步预测和单步预测也都试了，最后选择了用LSTM网络进行预报，从0基础到最后写完毕设的论文，也不算容易。所以这里写篇文章记录一下，中间涉及到一些文件数据的读取处理及转换这边就不说了，毕竟主要想讲一讲我怎么跑通LSTM网络。
代码参考的就是Keras官方提供的文档，还是比较通俗易懂的。

2 环境配置

Python 3.8 + TensorFlow 2.3 + Keras 2.4，keras还是比较适合小白的，上手简单，网络只需要暴力叠加。
对于tensorflow，我当时分别配了两个环境，一个是windows10我自己的电脑上的，另一个是在Linux超算服务器上配置的环境。
当时因为我做的是全球的电离层预报，里面涉及到的数据比较多，经常几千个程序一起运行，数据量也很大，所以我拿自己的电脑进行调代码，做测试，代码写好后就拿服务器来进行结果的运算。配环境的过程我都记录了一下，参考下面：

• windows10 配置 TensorFlow 环境：
  https://blog.csdn.net/qq_40185784/article/details/111165906?spm=1001.2014.3001.5501

• 超算平台/服务器配置 TensorFlow 环境：
  https://blog.csdn.net/qq_40185784/article/details/119113802?spm=1001.2014.3001.5501

3 数据和方法

我选了两年多的数据，2018.1.1-2020.3.1。训练和验证时候用的2018和2019两年的数据，而2020年两个月的数据用来最后结果的评估分析。具体我在做的时候还用了几种不同方法来对全球的电离层进行建模，来对比同一个深度学习网络预报不同的东西得到最后全球的电离层结果的差异，这里就不展开说了。
评估的指标选的是比较常用的几个公式：（d_i 为网络输出的预报值，^d_i 为真值，n 代表统计天数）
$$\{\begin{array}{l}\begin{array}{l}\begin{array}{l}relativeerror=\frac{1}{n}{\sum }_{i=1}^n\frac{|d_i-{\hat{d}}_i|}{{\hat{d}}_i}\times 100\%\\ MAD=\frac{1}{n}{\sum }_{i=1}^n|d_i-{\hat{d}}_i|\\ bias=\frac{1}{n}{\sum }_{i=1}^n\left(d_i-{\hat{d}}_i\right)\end{array}\\ STD=\sqrt{\frac{{\sum }_{i=1}^n{\left(d_i-{\hat{d}}_i-bias\right)}^2}{n}}\end{array}\\ RMS=\sqrt{\frac{{\sum }_{i=1}^n{\left(d_i-{\hat{d}}_i\right)}^2}{n}}\end{array}$$

4 网络结构和一些公式

介绍一些LSTM的内部结构和我搭建的网络

LSTM内部结构

LSTM神经网络是循环神经网络(RNN)的改进网络，是为了解决RNN模型梯度弥散的问题而提出的。该网络可基于输入的时间序列、通过内在遗忘的方式来学习历史数据并输出预测数据，其内部结构如图4-1所示。

LSTM内部中共有三个控制门：输入门、输出门、遗忘门。细胞内通过这三个控制门选择性的记忆反馈实现时间上的记忆功能，并运用于进行预测。输入门确定需要在细胞状态中要保存的信息和待更新的信息：
$$\begin{gathered} {\Gamma }_t^i=\sigma \left(W_i\cdot \left[h_{t-1},X_t\right]+b_i\right) \\ {\tilde{C}}_t=\tan \text{h}\left(W_C\cdot \left[h_{t-1},X_t\right]+b_C\right) \end{gathered}$$
式中输入的是上一个细胞元的计算结果(上一次的隐向量) 和新一次的训练数据 。第一部分得到 确定哪些信息需要更新，第二部分通过一个tanh单元构建备选向量 ，用来添加至细胞状态中。为了得到细胞状态保存信息的 矩阵，引入sigmoid单元。

sigmoid单元的作用是将结果输出归算至0至1之间，表达对数据的更新信息。其中σ函数是sigmoid单元的形式，首先计算输入的线性组合，其次通过σ函数作为阈值单元来输出结果。σ函数的形式如下：
$$\sigma \left(y\right)=\frac{1}{1+e^{-y}}$$
而对于tanh单元，使用tanh函数可以将输入映射到[-1,1]范围内。其次，遗忘门要决定哪些信息要丢失：
$${\Gamma }_t^f=\sigma \left(W_f\cdot \left[h_{t-1},X_t\right]+b_f\right)$$
接收参数同样是*h_t-1*和 X_t，并且使用sigmoid单元对细胞状态 中每一个参数的输出值归算至0至1之间，表达对数据的“遗忘程度”，1表示“完全接受”，0表示“完全忽略”。通过遗忘门后生成的向量 与上一个神经细胞的细胞状态向量 相乘以决定丢失选择忘记的数据。然后添加上 与新的备选向量 的乘积对之前的状态进行更新，得到更新后的细胞状态 ：
$$C_t={\Gamma }_t^f\otimes C_{t-1}+{\Gamma }_t^i\otimes {\tilde{C}}_t$$
最后通过输出门来确定细胞输出的隐向量 ，输出依赖于本次的细胞状态：
$$\begin{gathered} {\Gamma }_t^O=\sigma \left(W_O\left[h_{t-1},X_t\right]+b_O\right) \\ {h}_t=O_t\times \tan \text{h}\left(C_t\right) \end{gathered}$$
确定细胞状态能够输出的部分 ${\Gamma }_t^O$ 后，将更新后的细胞状态 $C_t$ 输入 函数使得数据调整到[-1,1]，之后与输出门输出结果相乘，得到隐向量 。至此LSTM一个细胞的数据更新完毕。对于整个训练过程而言，输入数据 与上一个细胞元的隐向量分别于权重矩阵 $W_i,W_C,W_f,W_O$ 进行矩阵点乘操作，分别得到结果矩阵 ${\Gamma }_t^i,{\Gamma }_t^f,{\Gamma }_t^O,{\tilde{C}}_t$ ，其中 $W_i,W_C,W_f,W_O$ 矩阵内的参数就是该LSTM细胞的核心权重参数，训练的目的也就是通过结果来训练调整这些参数的值使得输出的结果符合预期与要求。

用于预报的网络

输入层我选了一些和电离层关系比较密切的参数，经过处理后输入网络。
隐藏层叠加了两层的LSTM，神经元用的都是32，最后叠一个dense是为了把输出的结果展开为一个纬度为24的矩阵，因为我输出的数据是连续24个的时间序列，如果是只预报一个时刻，这里需要改成dense=1。

输出层就是把输出的数据还原，因为前面有数据预处理的过程，这里要返回真实数据。

4 训练代码

首先是读取数据的代码，csv文件存储的是我处理好后的数据

import pandas as pd
folder = r'C:\Users\GLASSYM\Desktop\LSTM\\'

csvname = 'TEC_map.csv'
input_dex = 'input_dex.csv'

dex_data = pd.read_csv(folder+input_dex)
ionex_data = pd.read_csv(folder+csvname)

print(ionex_data.shape[0],ionex_data.shape[1])

18720 5187

之后是进行一个简单的数据查看

import numpy as np
pre_ion = ionex_data.iloc[:,2580].values
pre_ion = pre_ion[:,np.newaxis];
pre_dex = dex_data.iloc[:,3:7].values
pre_data = np.hstack((pre_ion,pre_dex))
print(pre_data, len(pre_data))

[[ 12.5 -5. 66.8 0. 1. ]
[ 10. -7. 66.8 0.25881905 0.96592583]
…
[ 24.2 -21. 69.4 -0.5 0.8660254 ]
[ 23. -22. 69.4 -0.25881905 0.96592583]] 18720

进行数据预处理，归一化

#两年数据作为训练集
data_mean = pre_data.mean(axis=0)
data_std = pre_data.std(axis=0)

deal_data = (pre_data - data_mean) / data_std
print(data_mean, data_std)
print(deal_data)

[ 1.45496314e+01 -5.39118590e+00 6.97937179e+01 4.60220656e-18 -3.70074342e-18]
[ 8.91233958 11.478995 2.85311328 0.70710678 0.70710678]
[[-2.29976808e-01 3.40784100e-02 -1.04928114e+00 -6.50850293e-18 1.41421356e+00]
…
[ 1.08280979e+00 -1.35977184e+00 -1.37995905e-01 -7.07106781e-01 1.22474487e+00]
[ 9.48165014e-01 -1.44688748e+00 -1.37995905e-01 -3.66025404e-01 1.36602540e+00]]

准备数据生成器，生成对应的训练集验证集和测试集，返回的是你用于训练的数据和对应的标签
这里要注意一下，你需要生成标签的矩阵维度就是你叠加在网络里面dense层的值，如果你预报一个值，那么标签里面一定有一个维度（代码里面是delay）为1，dense=1。如果连续预报12小时，那么标签也有个连续12小时的数据，那么有一个维度为12，dense=12。

# data ：浮点数据原始数组
# lookback ：输入数据包含的时间步
# delay ：目标在多少个时间步之后
# min_index，max_index ：data数组的索引，用于界定需要抽取哪些时间步。保存一部分用于验证，一部分用于测试
# shuffle ：决定是否打乱样本；TRUE：打乱 ；FLASE：按顺序处理样本
# batch_size ：每个批量的样本数
def generator(data, lookback, delay, min_index, max_index,
              shuffle=False, batch_size=128):
    if max_index is None:
        max_index = len(data) - delay - 1
    i = min_index + lookback
    while 1:
        if shuffle:
            rows = np.random.randint(
                min_index + lookback, max_index, size=batch_size)
        else:
            if i + batch_size >= max_index:
                i = min_index + lookback
            rows = np.arange(i, min(i + batch_size, max_index))
            i += len(rows)

        samples = np.zeros((len(rows),
                           lookback,
                           data.shape[-1]))
        targets = np.zeros((len(rows),delay))
        
        for j, row in enumerate(rows):
                indices_sam = range(rows[j] - lookback, rows[j])
                indices_tar = range(rows[j] + 1, rows[j] + 1 + delay)
                samples[j] = data[indices_sam]
                targets[j] = data[indices_tar][:,0]

        yield samples, targets

生成对应的训练集验证集和测试集，lookback，delay，batch_size这些数据后续我都是在脚本外传入程序的。

lookback = 72
delay = 24
batch_size = 128

train_gen = generator(deal_data,
                      lookback=lookback,
                      delay=delay,
                      min_index=0,
                      max_index=17016,
                      shuffle=True,
                      batch_size=batch_size)
val_gen = generator(deal_data,
                    lookback=lookback,
                    delay=delay,
                    min_index=17017,
                    max_index=18720,
                    batch_size=batch_size)
test_gen = generator(deal_data,
                     lookback=lookback,
                     delay=delay,
                     min_index=17449,
                     max_index=None,
                     batch_size=batch_size)

# This is how many steps to draw from `val_gen`
# in order to see the whole validation set:
val_steps = (18720 - 17017 - lookback) // batch_size

# This is how many steps to draw from `test_gen`
# in order to see the whole test set:
test_steps = (len(deal_data) - 17449 - lookback) // batch_size

构建训练网络，两层LSTM层，mode.fit里面的的verbose参数是控制日志输出的格式。
激活函数我用的是默认的tanh，按理说relu函数应该效果会好些，但是我测试出来relu和tanh区别不大，而用relu的时候cudnn无法进行加速，训练的速度就很慢（这个应该是keras函数内部的一些原因），但是默认tanh的话用cudnn加速训练速度提升了90%，我只能说无敌。
给的参数是 steps_per_epoch=200, epochs=20 ，测试的时候训练还是挺快的。

from keras.models import Sequential
from keras import layers
from keras.optimizers import RMSprop
from keras.layers import LeakyReLU
import numpy as np

model = Sequential()
model.add(layers.LSTM(32,
            return_sequences = True,
            input_shape = (None,deal_data.shape[-1])))
model.add(layers.LSTM(32))
model.add(layers.Dense(24))

model.compile(optimizer=RMSprop(), loss='mae')
history = model.fit(train_gen,
                    steps_per_epoch=200,
                    epochs=20,
                    validation_data=val_gen,
                    validation_steps=val_steps,
                    verbose = 1)

画出损失值的结果散点图

import matplotlib.pyplot as plt

loss = history.history['loss']
val_loss = history.history['val_loss']

epochs = range(len(loss))

plt.figure()
plt.subplots(figsize=(20,10))
plt.plot(epochs, loss, 'ro', label='Training loss')
plt.plot(epochs, val_loss, 'b', label='Validation loss')
plt.title('Training and validation loss')
plt.legend(fontsize=24)

plt.show()

对训练完的模型输入测试集test进行预测并输出损失值

results=model.predict(test_gen, 
                      steps=test_steps,
                      verbose=1)
loss=model.evaluate(test_gen,
                    steps=test_steps,)
print(loss)

还原数据，画结果对比图，注意这里是隔24个数据读取一次，即连续的预报24小时

real_data = pre_data[17449+72:,0] 
pre_test = []
len_result = len(results)
for i in range(int(len_result/24)):
    pre_one_day = results[i*24] * data_std[0] + data_mean[0]
    pre_test = np.append(pre_test,pre_one_day)
#pre_test = results[:,0] * data_std[0] + data_mean[0]

print(real_data,len(real_data))
print(pre_test,len(pre_test))

[13.5 9.3 7.4 … 23.7 24.2 23. ] 1199
[ 9.86116982 8.84894657 8.43008041 … 20.85349083 17.75438309 14.54913998] 1152

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(20,10),)
plt.plot(real_data[:300], 'o', color='r', label='- really -', markersize=10)
plt.plot(pre_test[:300], 'o', color='b', label='- predict -', markersize=10)
plt.legend(fontsize=20)

这里图片画了300个散点，可以看出结果还是比价好的，如果要放论文的画必须得调整一下大小和图例文字分辨率等，这边就随便展示一下。

保存模型

# save model
LSTM_model = 'C:\Users\GLASSYM\Desktop\IONEX_file\model_lstm.h5'
model.save(LSTM_model)

这样基本上一个简单的LSTM训练和预报就完成了，现在回想起来，其实整个代码的难度不高，都很简单，只是刚开始看的时候容易漏掉很多细节，老出bug，其实注意一下细节就好。

5 小结

• 使用的数据是最近几年的数据，而近几年太阳均处在平静期，太阳活动不强。对于太阳活动活跃期，网络的结构是否需要相应的变化以及网络的预报效果仍有待检验与测试。
• 实验均为短期(1天)的预报，能否进行一个长期连续的预测预报，以及能否根据对历史数据的训练来对历史缺失的数据进行修复仍有待研究。

上面的东西都是之前写的和做的，暑假最近在跑一个新的方法，看看能不能在空间上对预报效果进行提升。也许有用也许没用，反正我觉的深度学习这东西就是个玄学，练练丹。也不知道未来三年的研究生能研究出个啥，也许会结合定位做点东西。Wish me good luck!