1200字范文 > 【算法竞赛学习】AI助力精准气象和海洋预测

【算法竞赛学习】AI助力精准气象和海洋预测

时间：2024-04-27 07:50:54

赛题简介

赛题背景

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

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

数据描述

数据简介

本次比赛使用的数据包括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=23∗accskill−RMSEScore = \frac{2}{3} * accskill - RMSE Score=32∗accskill−RMSE

其中，

accskill=∑i=124a∗ln(i)∗cori,(i≤,a=1.5;5≤i≤11,a=2;12≤i≤18,a=3;19≤i,a=4)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) accskill=i=1∑24a∗ln(i)∗cori,(i≤,a=1.5;5≤i≤11,a=2;12≤i≤18,a=3;19≤i,a=4)

而：

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

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

线下数据转换

将数据转化为我们所熟悉的形式，每个人的风格不一样，此处可以作为如何将nc文件转化为csv等文件

数据转化

## 工具包导入&数据读取### 工具包导入'''安装工具# !pip install netCDF4 ''' import pandas as pdimport numpy as npimport tensorflow as tffrom tensorflow.keras.optimizers import Adamimport matplotlib.pyplot as pltimport scipy from netCDF4 import Datasetimport netCDF4 as ncimport gc%matplotlib inline

数据读取

SODA_label处理

标签含义

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

将标签转化为我们熟悉的pandas形式

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_indexdf_SODA_label['label']= vsdf_SODA_label.to_csv(label_trans_path + 'df_SODA_label.csv',index = None)

df_SODA_label.head()

转化

SODA_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] = vreturn 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})

%%timedf_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处理

label_path = './data/CMIP_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_CMIP_label= pd.DataFrame({'year_month':year_month_index}) df_CMIP_label['year_month'] = year_month_indexdf_CMIP_label['label']= vsdf_CMIP_label.to_csv(label_trans_path + 'df_CMIP_label.csv',index = None)

df_CMIP_label.head()

CMIP_train处理

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度），所有数据的经纬度范围相同。

CMIP_path = './data/CMIP_train.nc'CMIP_trans_path = './data'nc_CMIP = Dataset(CMIP_path,'r')

nc_CMIP.variables.keys()

dict_keys(['sst', 't300', 'ua', 'va', 'year', 'month', 'lat', 'lon'])

nc_CMIP['t300'][:].shape

(4645, 36, 24, 72)

year_month_index = []years = np.array(nc_CMIP['year'][:])months = np.array(nc_CMIP['month'][:])lats= np.array(nc_CMIP['lat'][:])lons= np.array(nc_CMIP['lon'][:])last_thre_years = 1000for year in years:'''数据的原因，我们'''if year >= 4645 - last_thre_years:for month in months:year_month_index.append('year_{}_month_{}'.format(year,month))df_CMIP_sst = pd.DataFrame({'year_month':year_month_index}) df_CMIP_t300 = pd.DataFrame({'year_month':year_month_index}) df_CMIP_ua = pd.DataFrame({'year_month':year_month_index}) df_CMIP_va = pd.DataFrame({'year_month':year_month_index})

因为内存限制,我们暂时取最后1000个year的数据

def trans_thre_df(df, vals, lats, lons, years, months, last_thre_years = 1000):'''(4645, 36, 24, 72) -- year, month,lat,lon ''' for j,lat_ in (enumerate(lats)):# print(j)for i,lon_ in enumerate(lons):c = 'lat_lon_{}_{}'.format(int(lat_),int(lon_)) v = []for y_,y in enumerate(years):'''数据的原因，我们'''if y >= 4645 - last_thre_years:for m_,m in enumerate(months): v.append(vals[y_,m_,j,i])df[c] = vreturn df

%%timedf_CMIP_sst = trans_thre_df(df = df_CMIP_sst, vals = np.array(nc_CMIP['sst'][:]), lats = lats, lons = lons, years = years, months = months)df_CMIP_sst.to_csv(CMIP_trans_path + 'df_CMIP_sst.csv',index = None)del df_CMIP_sstgc.collect()df_CMIP_t300 = trans_thre_df(df = df_CMIP_t300, vals = np.array(nc_CMIP['t300'][:]), lats = lats, lons = lons, years = years, months = months)df_CMIP_t300.to_csv(CMIP_trans_path + 'df_CMIP_t300.csv',index = None)del df_CMIP_t300gc.collect()df_CMIP_ua = trans_thre_df(df = df_CMIP_ua, vals = np.array(nc_CMIP['ua'][:]), lats = lats, lons = lons, years = years, months = months)df_CMIP_ua.to_csv(CMIP_trans_path + 'df_CMIP_ua.csv',index = None)del df_CMIP_uagc.collect()df_CMIP_va = trans_thre_df(df = df_CMIP_va, vals = np.array(nc_CMIP['va'][:]), lats = lats, lons = lons, years = years, months = months)df_CMIP_va.to_csv(CMIP_trans_path + 'df_CMIP_va.csv',index = None)del df_CMIP_vagc.collect()

(36036, 1729)

数据建模

工具包导入&数据读取

工具包导入

import pandas as pdimport numpy as npimport tensorflow as tffrom tensorflow.keras.optimizers import Adam import joblibfrom netCDF4 import Datasetimport netCDF4 as ncimport gcfrom sklearn.metrics import mean_squared_errorimport numpy as npfrom tensorflow.keras.callbacks import LearningRateScheduler, Callbackimport tensorflow.keras.backend as Kfrom tensorflow.keras.layers import *from tensorflow.keras.models import *from tensorflow.keras.optimizers import *from tensorflow.keras.callbacks import *from tensorflow.keras.layers import Input %matplotlib inline

数据读取

SODA_label处理

标签含义

将标签转化为我们熟悉的pandas形式

df_SODA_label = pd.read_csv('./data/df_SODA_label.csv')df_CMIP_label = pd.read_csv('./data/df_CMIP_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]

模型构建

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) pile(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,424Trainable params: 550,424Non-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)

Metrics

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 = 0rmse_score= 0a = [1.5] * 4 + [2] * 7 + [3] * 7 + [4] * 6y_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 / fenmuaccskill_score += a[i] * np.log(i+1) * cor_irmse_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))

模型预测

模型构建

在上面的部分，我们已经训练好了模型，接下来就是提交模型并在线上进行预测，这块可以分为三步：

导入模型；读取测试数据并且进行预测；生成提交所需的版本；

import tensorflow as tfimport tensorflow.keras.backend as Kfrom tensorflow.keras.layers import *from tensorflow.keras.models import *from tensorflow.keras.optimizers import *from tensorflow.keras.callbacks import *from tensorflow.keras.layers import Input import numpy as npimport osimport zipfiledef RMSE(y_true, y_pred):return tf.sqrt(tf.reduce_mean(tf.square(y_true - y_pred)))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) pile(optimizer=adam, loss=RMSE)return modelfeature_cols = ['month_{}'.format(i+1) for i in range(12)]label_cols = ['month_{}'.format(i+1) for i in range(12, 36)] model = build_model(train_feat=feature_cols, test_feat=label_cols)model.load_weights('./user_data/model_data/model_mlp_baseline.h5')

模型预测

test_path = './tcdata/enso_round1_test_0201/'### 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:continuepathfile = os.path.join(parent, filename)arcname = pathfile[pre_len:].strip(os.path.sep) #相对路径zipf.write(pathfile, arcname)zipf.close()make_zip()