关于机器学习:人工智能创新挑战赛助力精准气象和海洋预测Baseline2数据探索性分析温度风场可视化CNNLSTM模型建模

“AI Earth”人工智能翻新挑战赛：助力精准气象和陆地预测 Baseline[2]：数据探索性剖析（温度风场可视化）、CNN+LSTM 模型建模

1. 气象陆地预测 - 数据分析

数据分析是解决一个数据挖掘工作的重要一环，通过数据分析，咱们能够理解标签的散布、数据中存在的缺失值和异样值、特色与标签之间的相关性、特色之间的相关性等，并依据数据分析的后果，领导咱们后续的特色工程以及模型的抉择和设计。

在本次工作中，将摸索赛题中给出的两份训练数据，可视化剖析四个气象特色的散布状况，思考如何进行特色工程以及如何抉择或设计模型来实现咱们的预测工作。

• 学习指标

  1. 学习如何摸索并可视化剖析气象数据。

  2. 依据数据分析后果思考以下两个问题：

     • 是否结构新的特色？
     • 抉择或设计什么样的模型进行预测？

本赛题应用的训练数据包含 CMIP5 中 17 个模式提供的 140 年的历史模仿数据、CMIP6 中 15 个模式提供的 151 年的历史模仿数据和美国 SODA 模式重建的 100 年的历史观测异化数据，采纳 nc 格局保留，其中 CMIP5 和 CMIP6 别离是世界气象钻研打算 (WCRP) 的第 5 次和第 6 次耦合模式比拟打算，这二者都提供了多种不同的气象模式对于多种气象变量的模仿数据。这些数据蕴含四种气象变量：海表温度异样(SST)、热含量异样(T300)、纬向风异样(Ua)、经向风异样(Va)，数据维度为(year, month, lat, lon)，对于训练数据提供对应月份的 Nino3.4 指数标签数据。简而言之，提供的训练数据中的每个样本为某年、某月、某个维度、某个经度的 SST、T300、Ua、Va 数值，标签为对应年、对应月的 Nino3.4 指数。

须要留神的是，样本的第二维度 month 的长度不是 12 个月，而是 36 个月，对应从以后 year 开始间断三年的数据，例如 SODA 训练数据中 year 为 0 时蕴含的是从第 1 – 第 3 年逐月的历史观测数据，year 为 1 时蕴含的是从第 2 年 – 第 4 年逐月的历史观测数据，也就是说，样本在工夫上是有穿插的。

另外一点须要留神的是，Nino3.4 指数是 Nino3.4 区域从以后月开始间断三个月的 SST 平均值，也就是说，咱们也能够不间接预测 Nino3.4 指数，而是以 SST 为预测指标，间接求得 Nino3.4 指数。测试数据为国内多个陆地材料异化后果提供的随机抽取的 $N$ 段长度为 12 个月的工夫序列，数据采纳 npy 格局保留，维度为(12, lat, lon, 4)，第一维度为间断的 12 个月份，第四维度为 4 个气象变量，按 SST、T300、Ua、Va 的程序寄存。测试集文件序列的命名如 test_00001_01_12.npy 中 00001 示意编号，01 示意起始月份，12 示意终止月份。

在本次工作中咱们须要用到 global_land_mask 这个 Python 库，这个库能够依据输出的纬度和经度判断该点是否在海洋上。

我的项目链接以及码源见文末

# !unzip -d input /home/aistudio/data/data221707/enso_round1_train_20210201.zip

Archive:  /home/aistudio/data/data221707/enso_round1_train_20210201.zip
   creating: input/enso_round1_train_20210201/
  inflating: input/enso_round1_train_20210201/CMIP_train.nc  
  inflating: input/enso_round1_train_20210201/.DS_Store  
  inflating: input/enso_round1_train_20210201/readme.txt  
  inflating: input/enso_round1_train_20210201/SODA_train.nc  
  inflating: input/enso_round1_train_20210201/SODA_label.nc  
  inflating: input/enso_round1_train_20210201/CMIP_label.nc  

# 装置 global_land_mask
!pip install global_land_mask

import numpy as np
import matplotlib.pyplot as plt
import random
import netCDF4
import seaborn as sns
from global_land_mask import globe
from scipy import interpolate
plt.rcParams['font.sans-serif'] = ['SimHei'] #中文反对
%matplotlib inline

1.1 以 SODA 数据集为例，摸索标签的散布状况

# 读取 SODA 数据

# 存放数据的门路
path = '/home/aistudio/'
data = netCDF4.Dataset(path + 'SODA_train.nc')
label = netCDF4.Dataset(path + 'SODA_label.nc')
label = np.array(label.variables['nino'])

print(data.variables['sst'].shape)

(100, 36, 24, 72)

label.shape

(100, 36)

能够看到，数据集中的每个样本是以某一年为起始的接下来 36 个月的气象数据，同样的，标签也是以这一年为起始的接下来 36 个月中每个月的 Nino3.4 指数。然而一年只有 12 个月，怎么会有 36 个月的数据呢？

咱们无妨来查看一下，将每个样本中的 12 个月进行拼接时，Nino3.4 指数的变动曲线。

# 别离将样本中的 0 -12 月、12-24 月、24-36 月进行拼接，绘制 Nino3.4 指数的变动曲线
plt.figure(figsize=(8, 10))
plt.subplots_adjust(hspace=0.4, wspace=0.4)
for i in range(3):
    label_all = [label[j, 12*i:12*(i+1)] for j in range(label.shape[0])]
    label_all = np.concatenate(label_all, axis=0)
    plt.subplot(3, 1, i+1)
    plt.plot(label_all, 'k', linewidth=1)
    plt.xlabel('Time / month')
    plt.ylabel('Nino')
    plt.title('nino: month{}-{}'.format(12*i, 12*(i+1)))
plt.show()

/opt/conda/envs/python35-paddle120-env/lib/python3.7/site-packages/matplotlib/font_manager.py:1331: UserWarning: findfont: Font family ['sans-serif'] not found. Falling back to DejaVu Sans
  (prop.get_family(), self.defaultFamily[fontext]))

从图中能够看出，三条曲线的变化趋势是完全相同的，只是在工夫上有 12 个月的位移。这阐明，重叠局部的标签是雷同的，给出这样的样本的目标是能够以前 12 个月作为模型的输出 X、后 24 个月为预测指标 Y 构建训练样本。

进一步思考，将每个样本构建一个训练样本，那么咱们所能失去的训练数据量就只有 4645（CMIP 数据）+100（SODA 数据）=4745 条，这样小的数据量用于模型训练必然是不够的，那么如何减少数据量呢？这个问题留待大家思考。

同样的，咱们能够随机抽取三个样本，查看每个样本中的 Nino3.4 指数变动曲线。

# 随机抽取三个样本，绘制 Nino3.4 指数变动曲线
plt.figure(figsize=(8, 10))
plt.subplots_adjust(hspace=0.4, wspace=0.4)
for i in range(3):
    y = random.randint(0, label.shape[0])
    plt.subplot(3, 1, i+1)
    plt.plot(label[y], 'k', linewidth=1)
    plt.xlabel('Time / month')
    plt.ylabel('Nino')
    plt.title('nino: year{}'.format(y))
plt.show()

1.2 以 SST 特色为例，进行海陆掩膜和插值剖析

在给定数据中，经度和纬度坐标都是离散的，每隔 5 度有一个坐标点，在这样的经纬度坐标下的 SST 值也是离散的，因而咱们以样本 0 第 0 月的 SST 数据为例，用插值函数来拟合经纬度坐标与 SST 值之间的函数关系，失去平滑的 SST 散布。

lon = np.array(data.variables['lon'])
lat = np.array(data.variables['lat'])

x = lon
y = lat
# 以纬度和经度生成网格点坐标矩阵
xx, yy = np.meshgrid(x, y)
# 取样本 0 第 0 月的 SST 值
z = data.variables['sst'][0, 0]
# 采纳三次多项式插值，失去 z = f(x, y)的函数 f
f = interpolate.interp2d(x, y, z, kind = 'cubic')

x

array([  0.,   5.,  10.,  15.,  20.,  25.,  30.,  35.,  40.,  45.,  50.,
        55.,  60.,  65.,  70.,  75.,  80.,  85.,  90.,  95., 100., 105.,
       110., 115., 120., 125., 130., 135., 140., 145., 150., 155., 160.,
       165., 170., 175., 180., 185., 190., 195., 200., 205., 210., 215.,
       220., 225., 230., 235., 240., 245., 250., 255., 260., 265., 270.,
       275., 280., 285., 290., 295., 300., 305., 310., 315., 320., 325.,
       330., 335., 340., 345., 350., 355.])

经度的理论取值是从 -180°到 0°到 180°，而给定的数据中经度的取值是 0°到 360 度每距离 5°取一个坐标值，因而在之后判断是否为海洋时须要转换为理论的经度。

y

array([-55., -50., -45., -40., -35., -30., -25., -20., -15., -10.,  -5.,
         0.,   5.,  10.,  15.,  20.,  25.,  30.,  35.,  40.,  45.,  50.,
        55.,  60.])

数据中纬度的取值也是每距离 5°取一个坐标值。判断每个经纬度坐标点是否在海洋上，用空值遮蔽海洋局部，便于察看陆地上 SST 的散布。

# 将经度 x 转换为理论经度从新生成网格点坐标矩阵
lon_grid, lat_grid = np.meshgrid(x-180, y)
# 判断坐标矩阵上的网格点是否为海洋
is_on_land = globe.is_land(lat_grid, lon_grid)
is_on_land = np.concatenate([is_on_land[:, x >= 180], is_on_land[:, x < 180]], axis=1)
# 进行海洋掩膜，将海洋的 SST 设为空值
z[is_on_land] = np.nan
# 可视化掩膜后果，黄色为海洋，紫色为陆地
plt.imshow(is_on_land[::-1, :])

<matplotlib.image.AxesImage at 0x7f5d009a5950>

z.shape

(24, 72)

# 可视化陆地上的 SST 散布
plt.imshow(z[::-1, :], cmap=plt.cm.RdBu_r)

<matplotlib.image.AxesImage at 0x7f5d0090c390>

由上图能够看到，SST 的散布是离散的，咱们用之前失去的插值函数来平滑 SST 值，可视化平滑后的 SST 散布。

# 设置距离为 1°的经纬度坐标网格，用插值函数失去该坐标网格点的 SST 值
xnew = np.arange(0, 356, 1)
ynew = np.arange(-65, 66, 1)
znew = f(xnew, ynew)
lon_grid, lat_grid = np.meshgrid(xnew-180, ynew)
is_on_land = globe.is_land(lat_grid, lon_grid)
is_on_land = np.concatenate([is_on_land[:, xnew >= 180], is_on_land[:, xnew < 180]], axis=1)
# 同样进行海洋掩膜
znew[is_on_land] = np.nan
# 绘制平滑后的 SST 分布图
plt.imshow(znew[::-1, :], cmap=plt.cm.RdBu_r)

<matplotlib.image.AxesImage at 0x7f5d008efd50>

咱们用同样的办法绘制样本 0 中每个月的 SST 分布图，察看 SST 散布的变动。

# 绘制 0 年 36 个月的海陆掩膜
for i in range(1):
    plt.figure(figsize=(15, 18))
    for j in range(36):
        x = lon
        y = lat
        xx, yy = np.meshgrid(x, y)
        z = data.variables['sst'][i, j]
        f = interpolate.interp2d(x, y, z, kind='cubic')
        
        xnew = np.arange(0, 356, 1)
        ynew = np.arange(-65, 66, 1)
        znew = f(xnew, ynew)
        
        lon_grid, lat_grid = np.meshgrid(xnew-180, ynew)
        is_on_land = globe.is_land(lat_grid, lon_grid)
        is_on_land = np.concatenate([is_on_land[:, xnew >= 180], is_on_land[:, xnew < 180]], axis=1)
        znew[is_on_land] = np.nan
        plt.subplot(9, 4, j+1)
        plt.imshow(znew[::-1, :], cmap=plt.cm.RdBu_r)
        plt.title('sst - year:{}, month:{}'.format(i+1, j+1))

能够看到，SST 在每 12 个月中的前 4 个月和后 4 个月都较低，在两头 4 个月时较高，这阐明，海表温度在秋季和夏季较低，在冬季和秋季出现逐步升高到最高点而后逐步升高的变化趋势，这与咱们的认知常识是相符的。

大家也能够用同样的办法察看剖析其它三个气象特色的变化趋势。

1.3 以 CMIP 数据集为例，进行缺失值剖析

# 读取 CMIP 数据

# 存放数据的门路
path = '/home/aistudio/input/enso_round1_train_20210201/'
data = netCDF4.Dataset(path + 'CMIP_train.nc')
label = netCDF4.Dataset(path + 'CMIP_label.nc')
label = np.array(label.variables['nino'])

# 取得海洋的掩膜
lon_grid, lat_grid = np.meshgrid(x-180, y)
is_on_land = globe.is_land(lat_grid, lon_grid)
is_on_land = np.concatenate([is_on_land[:, x >= 180], is_on_land[:, x < 180]], axis=1)
mask = np.zeros(data.variables['sst'].shape, dtype=int)
mask[:, :, :, :] = is_on_land[np.newaxis, np.newaxis, :, :]

# 查看 SST 特色的缺失值数量
name = 'sst'
data_ = np.array(data.variables[name])
before_nan = np.sum(np.isnan(data_))
print('before:', before_nan)

before: 0

# 查看 T300 特色的缺失值数量
name = 't300'
data_ = np.array(data.variables[name])
before_nan = np.sum(np.isnan(data_))
print('before:', before_nan)

before: 3055032

# 查看 Va 特色的缺失值数量
name = 'va'
data_ = np.array(data.variables[name])
before_nan = np.sum(np.isnan(data_))
print('before:', before_nan)

before: 13921123

# 查看 Ua 特色的缺失值数量
name = 'ua'
data_ = np.array(data.variables[name])
before_nan = np.sum(np.isnan(data_))
print('before:', before_nan)

before: 13921123

四个气象特色中，SST 特色不存在缺失值，Va 和 Ua 特色中的缺失值数量最多。

接下来以 Ua 特色为例，可视化剖析缺失值的状况。

# 统计每年每月中 Ua 特色的缺失值数量
m = np.zeros(data_.shape[0:2])
for i in range(data_.shape[0]):
    for j in range(data_.shape[1]):
        if np.sum(np.isnan(data_[i][j])) != 0:
            m[i, j] = np.sum(np.isnan(data_[i, j]))

# 计算每一年的缺失值
before = np.sum(m, axis=1)
# 可视化每一年的缺失值数量
plt.plot(before, 'k')
plt.ylabel('nan count')
plt.xlabel('year')
plt.show()

能够看到在某些年份中存在较多缺失值。

# 查看 Ua 特色中存在缺失值的年数
len(np.where(before!=0)[0])

755

咱们取样本 1900 来察看 Ua 特色的散布。

# 可视化样本 1900 中 Ua 特色的散布
plt.imshow(data_[1900, 0][::-1, :])

上图中红色局部即为缺失值，能够看到，缺失值少数散布在海洋上，咱们将海洋局部进行填充，察看填充后 Ua 的散布。

# 将海洋地位填 0
data_[mask==1] = 0

# 可视化填充后样本 1900 中 Ua 特色的散布
plt.imshow(data_[1900, 0][::-1, :])

<matplotlib.image.AxesImage at 0x7f5c18dc4d10>

对海洋局部进行填充后缺失值数量大大减少。

# 统计填充后缺失值的数量
after_nan = np.sum(np.isnan(data_))

print('before: %d \nafter: %d \npercentage: %f'%(before_nan, after_nan, 1 - float(after_nan) / before_nan))

before: 13921123 
after: 2440742 
percentage: 0.824673

海洋局部填充解决了 82% 的缺失值。

# 可视化填充后每一年的缺失值数量
m = np.zeros(data_.shape[0: 2])
for i in range(data_.shape[0]):
    for j in range(data_.shape[1]):
        if np.sum(np.isnan(data_[i, j])) != 0:
            m[i, j] = np.sum(np.isnan(data_[i, j]))
after = np.sum(m, axis=1)
plt.plot(after, 'k')
plt.ylabel('nan count')
plt.xlabel('year')
plt.show()

# 比照填充前后每一年缺失值的数量
plt.plot(before, 'k')
plt.plot(after, 'r')
plt.legend(['before', 'after'])
plt.title(name)
plt.ylabel('nan count')
plt.xlabel('year')
plt.show()

1.4 温度场和风场可视化

在气象问题中，温度与风向往往是密切相关的。当温度越高时，气压越低，空气向上流动，温度越低时，气压越高，空气向下流动，于是温度高的中央上方的空气就会向温度低的中央流动，造成风。因而在剖析气象问题时，咱们往往会把温度和风向放在一起进行可视化。

如何把风向可视化呢？这里咱们要用到 plt.quiver()这个函数。

plt.quiver()用于绘制二维的向量场，次要输出参数有：

• X：向量起始点的 X 轴坐标
• Y：向量起始点的 Y 轴坐标
• U：向量的 X 轴重量
• V：向量的 Y 轴重量

具体的用法能够参考官网：https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.q…

# 对温度场 SST 进行插值，失去插值函数
x = lon
y = lat
xx, yy = np.meshgrid(x, y)
z = data.variables['sst'][0, 0]
f = interpolate.interp2d(x, y, z, kind='cubic')

# 取得海洋掩膜
lon_grid, lat_grid = np.meshgrid(x-180, y)
is_on_land = globe.is_land(lat_grid, lon_grid)
is_on_land = np.concatenate([is_on_land[:,x>=180], is_on_land[:,x<180]], axis=1)

# 对 Ua 和 Va 进行海洋掩膜
ua = data.variables['ua'][0, 0]
ua[is_on_land] = np.nan
va = data.variables['va'][0, 0]
va[is_on_land] = np.nan

# 插值后生成平滑的 SST 散布
xnew = np.arange(0, 356, 1)
ynew = np.arange(-65, 66, 1)
znew = f(xnew, ynew)

# 对平滑后的 SST 进行海洋掩膜
lon_grid, lat_grid = np.meshgrid(xnew-180, ynew)
is_on_land = globe.is_land(lat_grid, lon_grid)
is_on_land = np.concatenate([is_on_land[:, xnew >= 180], is_on_land[:, xnew < 180]], axis=1)
znew[is_on_land] = np.nan

# 绘制温度场
plt.figure(figsize=(15, 10))
plt.imshow(znew[::-1, :], cmap=plt.cm.RdBu_r)
plt.colorbar(orientation='horizontal') # 显示程度色彩条
# 绘制风向场，其实这里精确来说绘制的是风向异样的向量，而非理论的风向
plt.quiver(lon, lat+65, ua[::-1, :], va[::-1, :], alpha=0.8) # 在坐标 (lon, lat) 处绘制与 sqrt(ua^2, va^2)成正比长度的箭头
plt.title('year0-month0: SST/UA-VA')
plt.show()

从上图中能够看出，温度异样 SST 在 0 值左近时没有显著的风向异样，而在其余区域风向异样通常由 SST 值大的中央指向 SST 值小的中央。

ENSO 景象是指在温度场上赤道东太平洋温度继续异样增暖，在风向场上寒带东太平洋与寒带西太平洋气压变动（体现为风向）相同的景象。在上图这个样本中没有呈现 ENSO 景象，大家能够用同样的办法绘制并察看存在 ENSO 景象（Nino3.4 指数间断 5 个月超过 0.5℃）的样本的温度和风场。

1.5 小结

在以上的数据分析中，不难看出咱们在剖析气象问题时，采纳的依然是通用的数据分析思路：剖析标签 -> 剖析特色（包含特色散布、特色与特色的关系、特色与标签的关系）-> 剖析数据的根本状况（包含缺失值、异样值、反复值等）。这个思路大家能够灵便利用到各种问题的剖析中去，不至于拿到数据后无从下手。

通过以上的数据分析，咱们能够失去以下论断：

1. 重叠局部的标签是雷同的，为了减少数据量，咱们能够从每条数据中取固定的 12 个月拼接起来，用滑窗构建训练数据集。
2. SST 特色中没有缺失值，在其余特色中，缺失值基本上在海洋局部，将海洋局部用 0 填充能够解决绝大部分的缺失值。

当初咱们回到开篇的学习指标中提出的第一个问题：是否结构新的特色？目前的答案是不能。从各位 TOP 选手的计划以及相干的 ENSO 预测的论文来看，大家的眼光都聚焦在如何构建模型上，而鲜少有人会去结构新的特色。这与其说是不能，更不如说是不必要，因为个别咱们结构新的特色是为了从给出的特色中失去与预测指标更相干的特色，由此进步模型的学习效果，然而就本赛题而言，结构统计特色或者其余新的特色收效不高，咱们更心愿通过模型来开掘给定数据之间在工夫和空间上的依赖关系。

于是，解题的重点就落到了第二个问题上：抉择或设计什么样的模型进行预测？在接下来的三个工作中，咱们就要来学习 TOP 选手们的建模计划。

参考文献

1. 优良选手的 EDA 分享：https://gitee.com/Little_Six/aiweather-ocean-forecasts/blob/m…
2. Climate Data Guide：https://climatedataguide.ucar.edu/climate-data/nino-sst-indic…
3. 中国气象局国家气象核心：https://cmdp.ncc-cma.net/pred/cn_enso.php?product=cn_enso_nin…
4. WORLD-CLASS RESEARCH IN EARTH SYSTEM SCIENCE(世界一流的地球零碎科学研究): https://ncar.ucar.edu/

2. 气象陆地预测 - 模型建设之 CNN+LSTM

本次工作咱们将学习来自 TOP 选手“学习 AI 的打工人”的建模计划，该计划中采纳的模型是 CNN+LSTM。

在本赛题中，咱们结构的模型须要实现两个工作，开掘空间信息以及开掘工夫信息。那么，说到开掘空间信息的模型，咱们会很天然的想到 CNN，同样的，开掘工夫信息的模型咱们会很容易想到 LSTM，咱们本次学习的这个 TOP 计划正是结构了 CNN+LSTM 的串行构造。

• 学习指标

  1. 学习 TOP 计划的数据处理办法。
  2. 学习 TOP 计划的模型构建办法。

2.1 数据处理

该 TOP 计划的数据处理次要包含四局部：

1. 减少月特色。将序列数据的起始月份作为新的特色。
2. 数据扁平化。将序列数据按月拼接起来通过滑窗减少数据量。
3. 空值填充。
4. 结构数据集。随机采样结构数据集。

import netCDF4 as nc
import random
import os
from tqdm import tqdm
import math
import pandas as pd
import numpy as np

import seaborn as sns
color = sns.color_palette()
sns.set_style('darkgrid')
import matplotlib.pyplot as plt
%matplotlib inline

import torch
from torch import nn, optim
import torch.nn.functional as F
from torch.utils.data import Dataset, DataLoader

from sklearn.metrics import mean_squared_error

Cannot run import torch because of system compatibility. AI Studio prepared an entire environment based on PaddlePaddle already. Please use PaddlePaddle to build your own model or application.

# 固定随机种子
SEED = 22

def seed_everything(seed=42):
    random.seed(seed)
    os.environ['PYTHONHASHSEED'] = str(seed)
    np.random.seed(seed)
    torch.manual_seed(seed)
    torch.cuda.manual_seed(seed)
    torch.backends.cudnn.deterministic = True
    
seed_everything(SEED)

# 查看 CUDA 是否可用
train_on_gpu = torch.cuda.is_available()

if not train_on_gpu:
    print('CUDA is not available.  Training on CPU ...')
else:
    print('CUDA is available!  Training on GPU ...')

CUDA is available!  Training on GPU ...

# 读取数据
# 存放数据的门路
path = '/kaggle/input/ninoprediction/'
soda_train = nc.Dataset(path + 'SODA_train.nc')
soda_label = nc.Dataset(path + 'SODA_label.nc')
cmip_train = nc.Dataset(path + 'CMIP_train.nc')
cmip_label = nc.Dataset(path + 'CMIP_label.nc')

2.1.2 减少月特色

本赛题的线上测试集是任意选取某个月为起始的长度为 12 的序列，因而该计划中减少了起始月份作为新的特色。然而应用整数 1~12 不能反映 12 月与 1 月相邻这一特点，因而须要借助三角函数的周期性，同时思考到独自应用 sin 函数或 cos 函数会存在某些月份的函数值雷同的景象，因而同时应用 sin 函数和 cos 函数作为两个新增月份特色，保障每个起始月份的这两个特色组合都是举世无双的，并且又可能很好地体现出月份的周期性特色。

咱们能够通过可视化直观地感触下每个月份所结构的月份特色组合。

months = range(0, 12)
month_labels = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sept', 'Oct', 'Nov', 'Dec']
# sin 月份特色
months_sin = map(lambda x: math.sin(2 * math.pi * x / len(months)), months)
# cos 月份特色
months_cos = map(lambda x: math.cos(2 * math.pi * x / len(months)), months)

# 绘制每个月的月份特色组合
plt.figure(figsize=(20, 5))
x_axis = np.arange(-1, 13, 1e-2)
sns.lineplot(x=x_axis, y=np.sin(2 * math.pi * x_axis / len(months)))
sns.lineplot(x=x_axis, y=np.cos(2 * math.pi * x_axis / len(months)))
sns.scatterplot(x=months, y=months_sin, s=200)
sns.scatterplot(x=months, y=months_cos, s=200)
plt.xticks(ticks=months, labels=month_labels)
plt.show()

结构 SODA 数据的 sin 月份特色。

# 结构一个维度为 100*36*24*72 的矩阵，矩阵中的每个值为所在月份的 sin 函数值
soda_month_sin = np.zeros((100, 36, 24, 72))
for y in range(100):
    for m in range(36):
        for lat in range(24):
            for lon in range(72):
                soda_month_sin[y, m, lat, lon] = math.sin(2 * math.pi * (m % 12) / 12)
                
soda_month_sin.shape

(100, 36, 24, 72)

结构 SODA 数据的 cos 月份特色。

# 结构一个维度为 100*36*24*72 的矩阵，矩阵中的每个值为所在月份的 cos 函数值
soda_month_cos = np.zeros((100, 36, 24, 72))
for y in range(100):
    for m in range(36):
        for lat in range(24):
            for lon in range(72):
                soda_month_cos[y, m, lat, lon] = math.cos(2 * math.pi * (m % 12) / 12)
                
soda_month_cos.shape

(100, 36, 24, 72)

结构 CMIP 数据的 sin 月份特色。

# 结构一个维度为 4645*36*24*72 的矩阵，矩阵中的每个值为所在月份的 sin 函数值
cmip_month_sin = np.zeros((4645, 36, 24, 72))
for y in range(4645):
    for m in range(36):
        for lat in range(24):
            for lon in range(72):
                cmip_month_sin[y, m, lat, lon] = math.sin(2 * math.pi * (m % 12) / 12)
                
cmip_month_sin.shape

(4645, 36, 24, 72)

结构 CMIP 数据的 cos 月份特色。

# 结构一个维度为 4645*36*24*72 的矩阵，矩阵中的每个值为所在月份的 cos 函数值
cmip_month_cos = np.zeros((4645, 36, 24, 72))
for y in range(4645):
    for m in range(36):
        for lat in range(24):
            for lon in range(72):
                cmip_month_cos[y, m, lat, lon] = math.cos(2 * math.pi * (m % 12) / 12)
                
cmip_month_cos.shape

(4645, 36, 24, 72)

2.1.2 数据扁平化

在 Task2 中咱们发现，赛题中给出的数据量非常少，如何减少数据量呢？对于时序数据，一种罕用的做法就是滑窗。

因为每条数据在工夫上有重叠，咱们取数据的前 12 个月拼接起来，就失去了长度为（数据条数×12 个月）的序列数据，如图 1 所示：

而后咱们以每个月为起始月，接下来的 12 个月作为模型输出 X，后 24 个月的 Nino3.4 指数作为预测指标 Y 构建训练样本，如图 2 所示：

须要留神的是，CMIP 数据提供了不同的拟合模式，只有在同种模式下各个年份的数据在工夫上是间断的，因而同种模式的数据能力在工夫上拼接起来，除去最初 11 个月不能形成训练样本外，滑窗最终能取得的训练样本数量能够按以下形式计算失去：

• SODA：1 种模式×（100 年×12-11）=1189 条样本
• CMIP6：15 种模式×（151 年×12-11）=27015 条样本
• CMIP5：17 种模式×（140 年×12-11）=28373 条样本

在上面的代码中，咱们只将各个模式的数据拼接起来而没有采纳滑窗，这是因为思考到采纳滑窗失去的训练样本维度是（数据条数×12×24×72），须要占用大量的内存资源。咱们在之后构建数据集时，随机抽取了局部样本，大家在理论问题中，如果资源足够的话，能够采纳滑窗构建的全副的数据，不过须要留神数据量大的状况下能够思考构建更深的模型来开掘更多信息。

# 数据扁平化
def make_flatted(train_ds, label_ds, month_sin, month_cos, info, start_idx=0):
    keys = ['sst', 't300', 'ua', 'va']
    label_key = 'nino'
    # 年数
    years = info[1]
    # 模式数
    models = info[2]
    
    train_list = []
    label_list = []
    
    # 将同种模式下的数据拼接起来
    for model_i in range(models):
        blocks = []
        
        # 对每个特色，取每条数据的前 12 个月进行拼接
        for key in keys:
            block = train_ds[key][start_idx + model_i * years: start_idx + (model_i + 1) * years, :12].reshape(-1, 24, 72, 1).data
            blocks.append(block)
        # 减少 sin 月份特色
        block_sin = month_sin[start_idx + model_i * years: start_idx + (model_i + 1) * years, :12].reshape(-1, 24, 72, 1)
        blocks.append(block_sin)
        # 减少 cos 月份特色
        block_cos = month_cos[start_idx + model_i * years: start_idx + (model_i + 1) * years, :12].reshape(-1, 24, 72, 1)
        blocks.append(block_cos)
        
        # 将所有特色在最初一个维度上拼接起来
        train_flatted = np.concatenate(blocks, axis=-1)
        
        # 取 12-23 月的标签进行拼接，留神加上最初一年的最初 12 个月的标签（与最初一年 12-23 月的标签独特形成最初一年前 12 个月的预测指标）label_flatted = np.concatenate([label_ds[label_key][start_idx + model_i * years: start_idx + (model_i + 1) * years, 12: 24].reshape(-1).data,
            label_ds[label_key][start_idx + (model_i + 1) * years - 1, 24: 36].reshape(-1).data
        ], axis=0)
        
        train_list.append(train_flatted)
        label_list.append(label_flatted)
        
    return train_list, label_list

soda_info = ('soda', 100, 1)
cmip6_info = ('cmip6', 151, 15)
cmip5_info = ('cmip5', 140, 17)

soda_trains, soda_labels = make_flatted(soda_train, soda_label, soda_month_sin, soda_month_cos, soda_info)
cmip6_trains, cmip6_labels = make_flatted(cmip_train, cmip_label, cmip_month_sin, cmip_month_cos, cmip6_info)
cmip5_trains, cmip5_labels = make_flatted(cmip_train, cmip_label, cmip_month_sin, cmip_month_cos, cmip5_info, cmip6_info[1]*cmip6_info[2])

# 失去扁平化后的数据维度为（模式数×序列长度×纬度×经度×特色数），其中序列长度 = 年数×12
np.shape(soda_trains), np.shape(cmip6_trains), np.shape(cmip5_trains)

((1, 1200, 24, 72, 6), (15, 1812, 24, 72, 6), (17, 1680, 24, 72, 6))

del soda_month_sin, soda_month_cos
del cmip_month_sin, cmip_month_cos
del soda_train, soda_label
del cmip_train, cmip_label

2.1.3 空值填充

在 Task2 中咱们发现，除 SST 外，其它特色中都存在空值，这些空值根本都在海洋上，因而咱们间接将空值填充为 0。

# 填充 SODA 数据中的空值
soda_trains = np.array(soda_trains)
soda_trains_nan = np.isnan(soda_trains)
soda_trains[soda_trains_nan] = 0
print('Number of null in soda_trains after fillna:', np.sum(np.isnan(soda_trains)))

Number of null in soda_trains after fillna: 0

# 填充 CMIP6 数据中的空值
cmip6_trains = np.array(cmip6_trains)
cmip6_trains_nan = np.isnan(cmip6_trains)
cmip6_trains[cmip6_trains_nan] = 0
print('Number of null in cmip6_trains after fillna:', np.sum(np.isnan(cmip6_trains)))

Number of null in cmip6_trains after fillna: 0

# 填充 CMIP5 数据中的空值
cmip5_trains = np.array(cmip5_trains)
cmip5_trains_nan = np.isnan(cmip5_trains)
cmip5_trains[cmip5_trains_nan] = 0
print('Number of null in cmip6_trains after fillna:', np.sum(np.isnan(cmip5_trains)))

Number of null in cmip6_trains after fillna: 0

2.2 结构数据集

在划分训练 / 验证集时，一个须要思考的问题是训练集、验证集、测试集三者的散布是否是统一的。在本赛题中咱们拿到的是两份数据，其中 CMIP 数据是 CMIP5/ 6 模式模仿的历史数据，SODA 数据是由 SODA 模式重建的的历史观测异化数据，线上测试集则是来自国内多个陆地材料的异化数据，由此看来，SODA 数据和线上测试集的散布是较为统一的，CMIP 数据的散布则与测试集不同。在三者不统一的状况下，咱们通常会尽可能使验证集与测试集的散布统一，这样当模型在验证集上有较好的体现时，在测试集上也会有较好的体现。

因而，咱们从 CMIP 数据的每个模式中各抽取 100 条数据作为训练集（这里抽取的样本数只是作为一个示例，理论模型训练的时候应用多少样本须要综合思考可用的资源条件和构建的模型深度），从 SODA 模式中抽取 100 条数据作为验证集。有的同学可能会纳闷，既然这里只用了 100 条 SODA 数据，那么为什么还要对 SODA 数据扁平化后再抽样而不间接用原始数据呢，因为间接取原始数据的前 12 个月作为输出，后 24 个月作为标签所失去的验证集每一条都是从 0 月开始的，而线上的测试集起始月份是随机抽取的，因而这里依然要尽可能保障验证集与测试集的数据分布统一，使构建的验证集的起始月份也是随机的。

咱们这里没有结构测试集，因为线上的测试集曾经公开了，能够间接应用，在较量时，线上的测试集是窃密的，须要结构线下的测试集来评估模型成果，同时须要留神线下的评估后果和线上的提交后果是否差距不大或者变化趋势是统一的，如果不是就须要调整线下的测试集，保障它和线上测试集的散布尽可能统一，可能较为精确地批示模型的调整方向。

# 结构训练集

X_train = []
y_train = []
# 从 CMIP5 的 17 种模式中各抽取 100 条数据
for model_i in range(17):
    samples = np.random.choice(cmip5_trains.shape[1]-12, size=100)
    for ind in samples:
        X_train.append(cmip5_trains[model_i, ind: ind+12])
        y_train.append(cmip5_labels[model_i][ind: ind+24])
# 从 CMIP6 的 15 种模式种各抽取 100 条数据
for model_i in range(15):
    samples = np.random.choice(cmip6_trains.shape[1]-12, size=100)
    for ind in samples:
        X_train.append(cmip6_trains[model_i, ind: ind+12])
        y_train.append(cmip6_labels[model_i][ind: ind+24])
X_train = np.array(X_train)
y_train = np.array(y_train)

# 结构测试集

X_valid = []
y_valid = []
samples = np.random.choice(soda_trains.shape[1]-12, size=100)
for ind in samples:
    X_valid.append(soda_trains[0, ind: ind+12])
    y_valid.append(soda_labels[0][ind: ind+24])
X_valid = np.array(X_valid)
y_valid = np.array(y_valid)

# 查看数据集维度
X_train.shape, y_train.shape, X_valid.shape, y_valid.shape

((3200, 12, 24, 72, 6), (3200, 24), (100, 12, 24, 72, 6), (100, 24))

del cmip5_trains, cmip5_labels
del cmip6_trains, cmip6_labels
del soda_trains, soda_labels

# 保留数据集
np.save('X_train_sample.npy', X_train)
np.save('y_train_sample.npy', y_train)
np.save('X_valid_sample.npy', X_valid)
np.save('y_valid_sample.npy', y_valid)

2.3 模型构建

在模型构建局部的通用流程是：结构评估函数 -> 构建并训练模型 -> 模型评估，后两步是循环的，能够依据评估后果从新调整并训练模型，再从新进行评估。

2.3.1 结构评估函数

模型的评估函数通常就是官网给出的评估指标，不过在较量中常常会呈现线下的评估后果和提交后的线上评估后果不统一的状况，这通常是线下测试集和线上测试集散布不统一造成的。

# 读取数据集
X_train = np.load('../input/ai-earth-task03-samples/X_train_sample.npy')
y_train = np.load('../input/ai-earth-task03-samples/y_train_sample.npy')
X_valid = np.load('../input/ai-earth-task03-samples/X_valid_sample.npy')
y_valid = np.load('../input/ai-earth-task03-samples/y_valid_sample.npy')

X_train.shape, y_train.shape, X_valid.shape, y_valid.shape

((3200, 12, 24, 72, 6), (3200, 24), (100, 12, 24, 72, 6), (100, 24))

# 结构数据管道
class AIEarthDataset(Dataset):
    def __init__(self, data, label):
        self.data = torch.tensor(data, dtype=torch.float32)
        self.label = torch.tensor(label, dtype=torch.float32)

    def __len__(self):
        return len(self.label)
    
    def __getitem__(self, idx):
        return self.data[idx], self.label[idx]

batch_size = 32

trainset = AIEarthDataset(X_train, y_train)
trainloader = DataLoader(trainset, batch_size=batch_size, shuffle=True)

validset = AIEarthDataset(X_valid, y_valid)
validloader = DataLoader(validset, batch_size=batch_size, shuffle=True)

2.3.2 模型结构与训练

这部分是赛题的重点，该 TOP 计划采纳的是 CNN+LSTM 的串行构造，其中 CNN 用来提取空间信息，LSTM 用来提取工夫信息，模型构造如下图所示。

• CNN 局部

CNN 罕用于解决图像信息，它在解决空间信息上也有很好的体现。CNN 的输出尺寸是（N,C,H,W），其中 N 是批量梯度降落中一个批次的样本数量，H 和 W 别离是输出图像的高和宽，C 是输出图像的通道数，对于本题中的空间数据，H 和 W 就对应数据的纬度和经度，C 对应特色数。咱们的训练样本中还多了一个工夫维度，因而须要用将输出数据的格局 (N,T,H,W,C) 转换为(N×T,C,H,W)。

BatchNormalization（前面简称 BN）是批标准化层，通常放在卷积层后用于标准化数据的散布，可能缩小各层不同数据分布之间的相互影响和依赖，具备放慢模型训练速度、防止梯度爆炸、在肯定水平上能加强模型泛化能力等长处，是神经网络问题中罕用的“大杀器”。不过目前对于 BN 层和 ReLU 激活函数的搁置程序孰先孰后的问题七嘴八舌，具体还是看模型的成果。对于这个问题的探讨能够参考 https://www.zhihu.com/question/283715823

总体来看 CNN 这一部分采纳的是比拟通用的构造，第一层采纳比拟大的卷积核（7×7），前面接多层的小卷积核（3×3），并用 BN 晋升模型成果，用池化层缩小模型参数、扩充感触野，池化层罕用的有 MaxPooling 和 AveragePooling，通常 MaxPooling 成果更好，不过具体看模型成果。模型的次要难点就在于调参，目前模型调参没有规范的答案，更多地是参考前人的教训以及一直地尝试。

• LSTM 局部

CNN 局部通过 Flatten 层将除工夫维度以外的维度压平（即除工夫步长 12 外的其它维度大小相乘，例如 CNN 局部最初的池化层输入维度是（N,T,C,H,W），则压平后的维度是（N,T,C×H×W）），输出 LSTM 层。LSTM 层承受的输出维度为（Time_steps，Input_size），其中 Time_steps 就是工夫步长 12，Input_size 是压平后的维度大小。Pytorch 中 LSTM 的主要参数是 input_size、hidden_size（隐层节点数）、batch_first（一个批次的样本数量 N 是否在第 1 维度），batch_first 为 True 时输出和输入的数据格式为(N,T,input_size/hidden_size)，为数据格式为(T,N,input_size/hidden_size)，须要留神的一点是 LSTM 的输入模式是 tensor 格局的 output 和 tuple 格局的(h_n,c_n)，其中 output 是所有工夫步的输入(N,T,hidden_size)，h_n 是隐层的输入（即最初一个工夫步的输入，格局为(1,N,hidden_size)），c_n 是记忆细胞 cell 的输入。因为咱们通过多层 LSTM 要取得的并非一个工夫序列，而是要抽取出一个对于输出序列的特色表白，因而最初咱们应用最初一个 LSTM 层的隐层输入 h_n 作为全连贯层的输出。LSTM 的应用办法能够参考 https://pytorch.org/docs/stable/generated/torch.nn.LSTM.html?…

因为 LSTM 有四个门，因而 LSTM 的参数量是 4 倍的 Input_size×hidden_size，参数量过多就容易过拟合，同时因为数据量也较少，因而该计划中只重叠了两个 LSTM 层。

# 结构模型
class Model(nn.Module):
    def __init__(self):
        super().__init__()
        self.conv1 = nn.Conv2d(6, 32, kernel_size=7, stride=2, padding=3)
        self.conv2 = nn.Conv2d(32, 32, kernel_size=3, stride=1, padding=1)
        self.bn = nn.BatchNorm2d(32)
        self.avgpool = nn.AvgPool2d(kernel_size=2, stride=2)
        self.flatten = nn.Flatten()
        self.lstm1 = nn.LSTM(3456, 2048, batch_first=True)
        self.lstm2 = nn.LSTM(2048, 1024, batch_first=True)
        self.fc = nn.Linear(1024, 24)
        
    def forward(self, x):
        # 转换输出形态
        N, T, H, W, C = x.shape
        x = x.permute(0, 1, 4, 2, 3).contiguous()
        x = x.view(N*T, C, H, W)
        
        # CNN 局部
        x = self.conv1(x)
        x = F.relu(self.bn(x))
        x = self.conv2(x)
        x = F.relu(self.bn(x))
        x = self.avgpool(x)
        x = self.flatten(x)

        # 留神 Flatten 层后输入为(N×T,C_new)，须要转换成(N,T,C_new)
        _, C_new = x.shape
        x = x.view(N, T, C_new)
        
        # LSTM 局部
        x, h = self.lstm1(x)
        x, h = self.lstm2(x)
        # 留神这里只应用隐层的输入
        x, _ = h
        
        x = self.fc(x.squeeze())
        
        return x

def rmse(y_true, y_preds):
    return np.sqrt(mean_squared_error(y_pred = y_preds, y_true = y_true))

# 评估函数
def score(y_true, y_preds):
    # 相关性技巧评分
    accskill_score = 0
    # RMSE
    rmse_scores = 0
    a = [1.5] * 4 + [2] * 7 + [3] * 7 + [4] * 6
    y_true_mean = np.mean(y_true, axis=0)
    y_pred_mean = np.mean(y_preds, axis=0)
    for i in range(24):
        fenzi = np.sum((y_true[:, i] - y_true_mean[i]) * (y_preds[:, i] - y_pred_mean[i]))
        fenmu = np.sqrt(np.sum((y_true[:, i] - y_true_mean[i])**2) * np.sum((y_preds[:, i] - y_pred_mean[i])**2))
        cor_i = fenzi / fenmu
        accskill_score += a[i] * np.log(i+1) * cor_i
        rmse_score = rmse(y_true[:, i], y_preds[:, i])
        rmse_scores += rmse_score
    return 2/3.0 * accskill_score - rmse_scores

model = Model()
print(model)

Model((conv1): Conv2d(6, 32, kernel_size=(7, 7), stride=(2, 2), padding=(3, 3))
  (conv2): Conv2d(32, 32, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1))
  (bn): BatchNorm2d(32, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True)
  (avgpool): AvgPool2d(kernel_size=2, stride=2, padding=0)
  (flatten): Flatten(start_dim=1, end_dim=-1)
  (lstm1): LSTM(3456, 2048, batch_first=True)
  (lstm2): LSTM(2048, 1024, batch_first=True)
  (fc): Linear(in_features=1024, out_features=24, bias=True)
)

思考到本次工作的评估指标 score=2/3×accskill-RMSE，其中 RMSE 是 24 个月的 rmse 的累计值，咱们这里能够自定义评估指标中的 RMSE 作为损失函数。

# 采纳 RMSE 作为损失函数
def RMSELoss(y_pred,y_true):
    loss = torch.sqrt(torch.mean((y_pred-y_true)**2, dim=0)).sum()
    return loss

model_weights = './task03_model_weights.pth'
device = 'cuda' if torch.cuda.is_available() else 'cpu'
model = Model().to(device)
criterion = RMSELoss
optimizer = optim.Adam(model.parameters(), lr=1e-3, weight_decay=0.001)  # weight_decay 是 L2 正则化参数
epochs = 10
train_losses, valid_losses = [], []
scores = []
best_score = float('-inf')
preds = np.zeros((len(y_valid),24))

for epoch in range(epochs):
    print('Epoch: {}/{}'.format(epoch+1, epochs))
    
    # 模型训练
    model.train()
    losses = 0
    for data, labels in tqdm(trainloader):
        data = data.to(device)
        labels = labels.to(device)
        optimizer.zero_grad()
        pred = model(data)
        loss = criterion(pred, labels)
        losses += loss.cpu().detach().numpy()
        loss.backward()
        optimizer.step()
    train_loss = losses / len(trainloader)
    train_losses.append(train_loss)
    print('Training Loss: {:.3f}'.format(train_loss))
    
    # 模型验证
    model.eval()
    losses = 0
    with torch.no_grad():
        for i, data in tqdm(enumerate(validloader)):
            data, labels = data
            data = data.to(device)
            labels = labels.to(device)
            pred = model(data)
            loss = criterion(pred, labels)
            losses += loss.cpu().detach().numpy()
            preds[i*batch_size:(i+1)*batch_size] = pred.detach().cpu().numpy()
    valid_loss = losses / len(validloader)
    valid_losses.append(valid_loss)
    print('Validation Loss: {:.3f}'.format(valid_loss))
    s = score(y_valid, preds)
    scores.append(s)
    print('Score: {:.3f}'.format(s))
    
    # 保留最佳模型权重
    if s > best_score:
        best_score = s
        checkpoint = {'best_score': s,
                      'state_dict': model.state_dict()}
        torch.save(checkpoint, model_weights)

# 绘制训练 / 验证曲线
def training_vis(train_losses, valid_losses):
    # 绘制损失函数曲线
    fig = plt.figure(figsize=(8,4))
    # subplot loss
    ax1 = fig.add_subplot(121)
    ax1.plot(train_losses, label='train_loss')
    ax1.plot(valid_losses,label='val_loss')
    ax1.set_xlabel('Epochs')
    ax1.set_ylabel('Loss')
    ax1.set_title('Loss on Training and Validation Data')
    ax1.legend()
    plt.tight_layout()

training_vis(train_losses, valid_losses)

咱们通常会绘制训练 / 验证曲线来察看模型的拟合状况，上图中咱们别离绘制了训练过程中训练集和验证集损失函数变动曲线。能够看到，训练集的损失函数降落很快，然而验证集的损失函数是震荡的，没有显著的降落，这阐明模型的学习效果较差，并存在过拟合问题，须要调整相干的参数。

2.3.3 模型评估

最初，咱们在测试集上评估模型的训练后果。

# 加载最佳模型权重
checkpoint = torch.load('../input/ai-earth-model-weights/task03_model_weights.pth')
model = Model()
model.load_state_dict(checkpoint['state_dict'])

<All keys matched successfully>

# 测试集门路
test_path = '../input/ai-earth-tests/'
# 测试集标签门路
test_label_path = '../input/ai-earth-tests-labels/'

import os

# 读取测试数据和测试数据的标签，并记录每个测试样本的起始月份用于之后结构月份特色
files = os.listdir(test_path)
X_test = []
y_test = []
first_months = []  # 样本起始月份
for file in files:
    X_test.append(np.load(test_path + file))
    y_test.append(np.load(test_label_path + file))
    first_months.append(int(file.split('_')[2]))

X_test = np.array(X_test)
y_test = np.array(y_test)
X_test.shape, y_test.shape

((103, 12, 24, 72, 4), (103, 24))

# 结构一个维度为 103*12*24*72 的矩阵，矩阵中的每个值为所在月份的 sin 函数值
test_month_sin = np.zeros((103, 12, 24, 72, 1))
for y in range(103):
    for m in range(12):
        for lat in range(24):
            for lon in range(72):
                test_month_sin[y, m, lat, lon] = math.sin(2 * math.pi * ((m + first_months[y]-1) % 12) / 12)
                
test_month_sin.shape

(103, 12, 24, 72, 1)

# 结构一个维度为 103*12*24*72 的矩阵，矩阵中的每个值为所在月份的 cos 函数值
test_month_cos = np.zeros((103, 12, 24, 72, 1))
for y in range(103):
    for m in range(12):
        for lat in range(24):
            for lon in range(72):
                test_month_cos[y, m, lat, lon] = math.cos(2 * math.pi * ((m + first_months[y]-1) % 12) / 12)
                
test_month_cos.shape

(103, 12, 24, 72, 1)

# 结构测试集
X_test = np.concatenate([X_test, test_month_sin, test_month_cos], axis=-1)
X_test.shape

(103, 12, 24, 72, 6)

testset = AIEarthDataset(X_test, y_test)
testloader = DataLoader(testset, batch_size=batch_size, shuffle=False)

# 在测试集上评估模型成果
model.eval()
model.to(device)
preds = np.zeros((len(y_test),24))
for i, data in tqdm(enumerate(testloader)):
    data, labels = data
    data = data.to(device)
    labels = labels.to(device)
    pred = model(data)
    preds[i*batch_size:(i+1)*batch_size] = pred.detach().cpu().numpy()
s = score(y_test, preds)
print('Score: {:.3f}'.format(s))

4it [00:00, 65.03it/s]

Score: 14.946

2.4 总结

• 该计划在数据处理局部采纳了滑窗来结构数据集，这是序列预测问题中罕用的减少数据量的办法。另外，该计划中减少了一组月份特色，集体认为在时序场景中减少的这组特色收益不高，更多的是通过模型开掘序列中的依赖关系，并且因为维度减少会使得训练数据占用的资源大大增加，对模型的成果晋升不显著。不过在其余场景中这种特色构造方法依然是值得借鉴的。
• 该计划没有抉择时空序列预测畛域的现有模型，而是抉择本人设计模型，计划中的这种结构模型的思路非常适合初学者学习，灵便地将不同模型串行或并行组合可能联合模型各自的劣势，这种模型构造方法须要留神的是一个模型的输入维度与另一个模型承受的输出维度要互相匹配。
• 尝试在不减少月份特色的状况下应用 CNN+LSTM 模型进行预测，同时尝试批改模型参数或层数比拟模型在测试集上的评分。
• 未能解决的问题

正如其余参赛者一样，咱们碰到最大的问题就是问题稳定问题。即便什么都不变，只是随机数的变动都会使得线上问题有十分大的变动(10 分高低都是可能的)。咱们尝试了许多伎俩，可能最容易想到的就是 SODA 和其余两个模式应该不在一个域上，所以须要通过迁徙学习或者域自适应来解决这个问题。然而尝试均以失败告终。

• 其余尝试

  • 减少模型层数(对于该问题，层数越多貌似成果越不好)
  • 批改模型中的卷积层，池化层等
  • 更换损失函数和优化器
  • 将 SST,T300,Ua,Va 分成 4 通道别离输出模型
  • 只应用 SST 的特色作为模型的输出
  • 应用域自适应做成一个双输入的模型：1 个失常的输入，另一个输入用来辨别 SODA，CMIP6 还是 CMIP5。
  • 尝试自注意力机制
  • 收集几个线上分数比拟好的模型，而后交融

我的项目链接以及码源见文末

链接跳转到文末，欢送订阅

更多文章请关注公重号：汀丶人工智能