时间序列预测(2)：AI助力精准气象和海洋预测

引言 本文主要是记载学习阿里天池大数据时间序列预测竞赛的baseline模型的构建过程，以及对baseline进行改进的一些思路，后续将会根据这些思路对模型进行调优。

一、背景

发生在热带太平洋上的厄尔尼诺-南方涛动(ENSO)现象是地球上最强、最显著的年际气候信号。通过大气或海洋遥相关过程，经常会引发洪涝、干旱、高温、雪灾等极端事件，对全球的天气、气候以及粮食产量具有重要的影响。准确预测ENSO，是提高东亚和全球气候预测水平和防灾减灾的关键。

本次赛题是一个时间序列预测问题。基于历史气候观测和模式模拟数据，利用T时刻过去12个月(包含T时刻)的时空序列（气象因子），构建预测ENSO的深度学习模型，预测未来1-24个月的Nino3.4指数，如下图所示：

二、数据预处理

本次的数据是由阿里天池大数据官网提供，共包含以下几个文件：

readme.txt
CMIP_label.nc  
CMIP_train.nc
SODA_train.nc
SODA_label.nc

.nc（network Common Data Format）文件是气象上常用的数据格式，python上读取.nc使用较多的库为netCDF4这个库，为了便于建模使用，首先需要对这些数据进行处理，将其转换为我们常用的.csv格式。

数据简介

本次比赛使用的数据包括CMIP5/6模式的历史模拟数据和美国SODA模式重建的近100多年历史观测同化数据。每个样本包含以下气象及时空变量：海表温度异常(SST)，热含量异常(T300)，纬向风异常（Ua），经向风异常（Va），数据维度为（year,month,lat,lon）。对于训练数据提供对应月份的Nino3.4 index标签数据。

每个数据样本第一维度（year）表征数据所对应起始年份，对于CMIP数据共4645年，其中1-2265为CMIP6中15个模式提供的151年的历史模拟数据（总共：151年 *15 个模式=2265）；2266-4645为CMIP5中17个模式提供的140年的历史模拟数据（总共：140年 *17 个模式=2380）。对于历史观测同化数据为美国提供的SODA数据。

其中每个样本第二维度（mouth）表征数据对应的月份，对于训练数据均为36，对应的从当前年份开始连续三年数据（从1月开始，共36月），比如：

• SODA_train.nc中[0,0:36,:,:]为第1-第3年逐月的历史观测数据；SODA_train.nc中[1,0:36,:,:]为第2-第4年逐月的历史观测数据；…,SODA_train.nc中[99,0:36,:,:]为第100-102年逐月的历史观测数据。

• CMIP_train.nc中[0,0:36,:,:]为CMIP6第一个模式提供的第1-第3年逐月的历史模拟数据；
  …,
  CMIP_train.nc中[150,0:36,:,:]为CMIP6第一个模式提供的第151-第153年逐月的历史模拟数据；
  CMIP_train.nc中[151,0:36,:,:]为CMIP6第二个模式提供的第1-第3年逐月的历史模拟数据；
  …,
  CMIP_train.nc中[2265,0:36,:,:]为CMIP5第一个模式提供的第1-第3年逐月的历史模拟数据；
  …,
  CMIP_train.nc中[2405,0:36,:,:]为CMIP5第二个模式提供的第1-第3年逐月的历史模拟数据；
  …,
  CMIP_train.nc中[4644,0:36,:,:]为CMIP5第17个模式提供的第140-第142年逐月的历史模拟数据。

其中每个样本第三、第四维度分别代表经纬度（南纬55度北纬60度，东经0360度），所有数据的经纬度范围相同。

训练数据标签说明

标签数据为Nino3.4 SST异常指数，数据维度为（year,month）。

CMIP(SODA)_train.nc对应的标签数据当前时刻Nino3.4 SST异常指数的三个月滑动平均值，因此数据维度与维度介绍同训练数据一致。

注：三个月滑动平均值为当前月与未来两个月的平均值。

测试数据说明

测试用的初始场（输入）数据为国际多个海洋资料同化结果提供的随机抽取的n段12个时间序列，数据格式采用NPY格式保存，维度为（12，lat，lon, 4）,12为t时刻及过去11个时刻，4为预测因子，并按照SST,T300,Ua,Va的顺序存放。

测试集文件序列的命名规则：test_编号_起始月份_终止月份.npy，如test_00001_01_12_.npy。

模型评估指标

评分细则说明： 根据所提供的n个测试数据，对模型进行测试，得到n组未来1-24个月的序列选取对应预测时效的n个数据与标签值进行计算相关系数和均方根误差，如下图所示。并计算得分。计算公式为:

Score = \frac{2}{3} * accskill - RMSE

其中，

accskill = \sum_{i=1}^{24} a * ln(i) * cor_i, \\
(i \le,a = 1.5; 5 \le i \le 11, a= 2; 12 \le i \le 18,a=3;19 \le i, a = 4)

而：

cor = \frac{\sum(X-\bar(X))\sum(Y-\bar(Y)}{\sqrt{\sum(X-\bar{X})^2)\sum(Y-\bar{Y})^2)}}

RMSE = \sum_{i=1}^{24} rmse_i, 

数据的转换

首先处理SODA_label数据吗，将其转换为csv格式数据：

label_path       = './data/SODA_label.nc'
label_trans_path = './data/' 
nc_label         = Dataset(label_path,'r')
 
years            = np.array(nc_label['year'][:])
months           = np.array(nc_label['month'][:])

year_month_index = []
vs               = []
for i,year in enumerate(years):
    for j,month in enumerate(months):
        year_month_index.append('year_{}_month_{}'.format(year,month))
        vs.append(np.array(nc_label['nino'][i,j]))

df_SODA_label               = pd.DataFrame({'year_month':year_month_index}) 
df_SODA_label['year_month'] = year_month_index
df_SODA_label['label']      = vs

df_SODA_label.to_csv(label_trans_path + 'df_SODA_label.csv',index = None)

其次，处理SODA_train数据，此处需要注意，train数据的结构有所不同。

SODA_train.nc中[0,0:36,:,:]为第1-第3年逐月的历史观测数据；

SODA_train.nc中[1,0:36,:,:]为第2-第4年逐月的历史观测数据；

…,

SODA_train.nc中[99,0:36,:,:]为第100-102年逐月的历史观测数据。

SODA_path        = './data/SODA_train.nc'
nc_SODA          = Dataset(SODA_path,'r') 

• 自定义抽取对应数据转化为df的形式；

index为年月; columns为lat和lon的组合

def trans_df(df, vals, lats, lons, years, months):
    '''
        (100, 36, 24, 72) -- year, month,lat,lon 
    ''' 
    for j,lat_ in enumerate(lats):
        for i,lon_ in enumerate(lons):
            c = 'lat_lon_{}_{}'.format(int(lat_),int(lon_))  
            v = []
            for y in range(len(years)):
                for m in range(len(months)): 
                    v.append(vals[y,m,j,i])
            df[c] = v
    return df

year_month_index = []

years              = np.array(nc_SODA['year'][:])
months             = np.array(nc_SODA['month'][:])
lats             = np.array(nc_SODA['lat'][:])
lons             = np.array(nc_SODA['lon'][:])


for year in years:
    for month in months:
        year_month_index.append('year_{}_month_{}'.format(year,month))

df_sst  = pd.DataFrame({'year_month':year_month_index}) 
df_t300 = pd.DataFrame({'year_month':year_month_index}) 
df_ua   = pd.DataFrame({'year_month':year_month_index}) 
df_va   = pd.DataFrame({'year_month':year_month_index})

分别对SODA_train数据的不同维度（sst, t300, ua, uv）应用转换函数，并将转换后的结果保存为csv文件。

%%time
df_sst = trans_df(df = df_sst, vals = np.array(nc_SODA['sst'][:]), lats = lats, lons = lons, years = years, months = months)
df_t300 = trans_df(df = df_t300, vals = np.array(nc_SODA['t300'][:]), lats = lats, lons = lons, years = years, months = months)
df_ua   = trans_df(df = df_ua, vals = np.array(nc_SODA['ua'][:]), lats = lats, lons = lons, years = years, months = months)
df_va   = trans_df(df = df_va, vals = np.array(nc_SODA['va'][:]), lats = lats, lons = lons, years = years, months = months)

label_trans_path = './data/'
df_sst.to_csv(label_trans_path  + 'df_sst_SODA.csv',index = None)
df_t300.to_csv(label_trans_path + 'df_t300_SODA.csv',index = None)
df_ua.to_csv(label_trans_path   + 'df_ua_SODA.csv',index = None)
df_va.to_csv(label_trans_path   + 'df_va_SODA.csv',index = None)

针对CMIP_label，CMIP_train数据的处理与上述处理方法类似，分别得到转换后的数据集：df_CMIP_sst, df_CMIP_t300, df_CMIP_ua, df_CMIP_va。

三、模型的构建和预测

构建训练集

根据df_SODA_label数据构建训练集：

df_SODA_label = pd.read_csv('./data/df_SODA_label.csv')

df_SODA_label['year']  = df_SODA_label['year_month'].apply(lambda x: x[:x.find('m') - 1])
df_SODA_label['month'] = df_SODA_label['year_month'].apply(lambda x: x[x.find('m') :])

df_train = pd.pivot_table(data = df_SODA_label, values = 'label',index = 'year', columns = 'month')
year_new_index    = ['year_{}'.format(i+1)  for i in range(df_train.shape[0])]
month_new_columns = ['month_{}'.format(i+1) for i in range(df_train.shape[1])]
df_train = df_train[month_new_columns].loc[year_new_index]

训练模型

再利用tensorflow构建一个多层感知机模型(MLP)：

def RMSE(y_true, y_pred):
    return tf.sqrt(tf.reduce_mean(tf.square(y_true - y_pred)))

def RMSE_fn(y_true, y_pred):
    return np.sqrt(np.mean(np.power(np.array(y_true, float).reshape(-1, 1) - np.array(y_pred, float).reshape(-1, 1), 2)))

def build_model(train_feat, test_feat): #allfeatures, 
    inp    = Input(shape=(len(train_feat)))  
    
    x = Dense(1024, activation='relu')(inp)  
    x = Dropout(0.25)(x) 
    x = Dense(512, activation='relu')(x)   
    x = Dropout(0.25)(x)  
    output = Dense(len(test_feat), activation='linear')(x)   
    model  = Model(inputs=inp, outputs=output)

    adam = tf.optimizers.Adam(lr=1e-3,beta_1=0.99,beta_2 = 0.99) 
    model.compile(optimizer=adam, loss=RMSE)

    return model

feature_cols = ['month_{}'.format(i+1) for i in range(12)]
label_cols   = ['month_{}'.format(i+1) for i in range(12, df_train.shape[1])] 

model_mlp = build_model(feature_cols, label_cols)
model_mlp.summary()

Model: "model"
_________________________________________________________________
Layer (type)                 Output Shape              Param #   
=================================================================
input_1 (InputLayer)         [(None, 12)]              0         
_________________________________________________________________
dense (Dense)                (None, 1024)              13312     
_________________________________________________________________
dropout (Dropout)            (None, 1024)              0         
_________________________________________________________________
dense_1 (Dense)              (None, 512)               524800    
_________________________________________________________________
dropout_1 (Dropout)          (None, 512)               0         
_________________________________________________________________
dense_2 (Dense)              (None, 24)                12312     
=================================================================
Total params: 550,424
Trainable params: 550,424
Non-trainable params: 0
_________________________________________________________________

将训练好的模型参数保存到文件中：

tr_len = int(df_train.shape[0] * 0.8)
tr_fea     = df_train[feature_cols].iloc[:tr_len,:].copy()
tr_label   = df_train[label_cols].iloc[:tr_len,:].copy()
 
val_fea     = df_train[feature_cols].iloc[tr_len:,:].copy()
val_label   = df_train[label_cols].iloc[tr_len:,:].copy() 


model_weights = './user_data/model_data/model_mlp_baseline.h5'

checkpoint = ModelCheckpoint(model_weights, monitor='val_loss', verbose=0, save_best_only=True, mode='min',
                             save_weights_only=True)

plateau        = ReduceLROnPlateau(monitor='val_loss', factor=0.5, patience=5, verbose=1, min_delta=1e-4, mode='min')
early_stopping = EarlyStopping(monitor="val_loss", patience=20)
history        = model_mlp.fit(tr_fea.values, tr_label.values,
                    validation_data=(val_fea.values, val_label.values),
                    batch_size=4096, epochs=200,
                    callbacks=[plateau, checkpoint, early_stopping],
                    verbose=2) 

模型评估

根据官方给出的评价规则对模型进行评价：

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_score     = 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_true,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]) 
    return  2 / 3.0 * accskill_score - rmse_score

y_val_preds = model_mlp.predict(val_fea.values, batch_size=1024)
print('score', score(y_true = val_label.values, y_preds = y_val_preds))

模型预测

根据官方给出的测试数据集，用前面训练出来的模型进行预测：

test_path = './tcdata/enso_round1_test_20210201/'

### 0. 模拟线上的测试集合
# for i in range(10):
#     x = np.random.random(12) 
#     np.save(test_path + "{}.npy".format(i+1),x)

### 1. 测试数据读取
files = os.listdir(test_path)
test_feas_dict = {}
for file in files:
    test_feas_dict[file] = np.load(test_path + file)
    
### 2. 结果预测
test_predicts_dict = {}
for file_name,val in test_feas_dict.items():
    test_predicts_dict[file_name] = model.predict(val.reshape([-1,12]))
#     test_predicts_dict[file_name] = model.predict(val.reshape([-1,12])[0,:])

### 3.存储预测结果
for file_name,val in test_predicts_dict.items(): 
    np.save('./result/' + file_name,val)

需要将预测结果打包到run.sh目录下方：

#打包目录为zip文件
def make_zip(source_dir='./result/', output_filename = 'result.zip'):
    zipf = zipfile.ZipFile(output_filename, 'w')
    pre_len = len(os.path.dirname(source_dir))
    source_dirs = os.walk(source_dir)
    print(source_dirs)
    for parent, dirnames, filenames in source_dirs:
        print(parent, dirnames)
        for filename in filenames:
            if '.npy' not in filename:
                continue
            pathfile = os.path.join(parent, filename)
            arcname = pathfile[pre_len:].strip(os.path.sep)   #相对路径
            zipf.write(pathfile, arcname)
    zipf.close()
make_zip() 

四、模型的改进方向

• 模型角度：baseline中只使用了简单的MLP模型进行建模，可以考虑使用其它的更加fancy的模型进行尝试；比如LSTM和transformer等。
• 数据层面：构建一些特征或者对数据进行一些数据变换等，本次构建模型由于本地算力的原因只使用了其中的部分数据，且还可以从时间的角度构造更多的特征；
• 针对损失函数设计各种trick的提升技巧；