賽題簡介

賽題背景

發生在熱帶太平洋上的厄爾尼諾-南方濤動(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=\frac{2}{3}?accskill?RMSE$$

其中，

$$\begin{gathered} accskill=\sum _{i=1}^{24}a?ln(i)?co{r}_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) \end{gathered}$$

而：

$$cor=\frac{\sum (X?\overline{(}X))\sum (Y?\overline{(}Y)}{\sqrt{\sum (X?\bar{X}{)}^2)\sum (Y?\bar{Y}{)}^2)}}$$

$$RMSE=\sum _{i=1}^{24}rms{e}_i,$$

線下數據轉換

• 將數據轉化為我們所熟悉的形式，每個人的風格不一樣，此處可以作為如何將nc文件轉化為csv等文件

數據轉化

## 工具包導入&數據讀取 ### 工具包導入''' 安裝工具 # !pip install netCDF4 ''' import pandas as pd import numpy as np import tensorflow as tf from tensorflow.keras.optimizers import Adam import matplotlib.pyplot as plt import scipy from netCDF4 import Dataset import netCDF4 as nc import 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_index df_SODA_label['label'] = vsdf_SODA_label.to_csv(label_trans_path + 'df_SODA_label.csv',index = None) df_SODA_label.head() year_monthlabel01234

year_1_month_1	-0.40720701217651367
year_1_month_2	-0.20244435966014862
year_1_month_3	-0.10386104136705399
year_1_month_4	-0.02910841442644596
year_1_month_5	-0.13252995908260345

轉化

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}) %%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處理

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_index df_CMIP_label['label'] = vsdf_CMIP_label.to_csv(label_trans_path + 'df_CMIP_label.csv',index = None) df_CMIP_label.head() year_monthlabel01234

year_1_month_1	-0.26102548837661743
year_1_month_2	-0.1332537680864334
year_1_month_3	-0.014831557869911194
year_1_month_4	0.10506672412157059
year_1_month_5	0.24070978164672852

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 = 1000 for 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 %%time df_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_sst gc.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_t300 gc.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_ua gc.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_va gc.collect() (36036, 1729)

數據建模

工具包導入&數據讀取

工具包導入

import pandas as pd import numpy as np import tensorflow as tf from tensorflow.keras.optimizers import Adam import joblib from netCDF4 import Dataset import netCDF4 as nc import gc from sklearn.metrics import mean_squared_error import numpy as np from tensorflow.keras.callbacks import LearningRateScheduler, Callback import tensorflow.keras.backend as K from 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處理

標簽含義

標簽數據為Nino3.4 SST異常指數，數據維度為（year,month）。 CMIP(SODA)_train.nc對應的標簽數據當前時刻Nino3.4 SST異常指數的三個月滑動平均值，因此數據維度與維度介紹同訓練數據一致 注：三個月滑動平均值為當前月與未來兩個月的平均值。

將標簽轉化為我們熟悉的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) 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)

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 tf import tensorflow.keras.backend as K from 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 np import os import 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) model.compile(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_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:continuepathfile = os.path.join(parent, filename)arcname = pathfile[pre_len:].strip(os.path.sep) #相對路徑zipf.write(pathfile, arcname)zipf.close() make_zip()

提升方向

• 模型角度：我們只使用了簡單的MLP模型進行建模，可以考慮使用其它的更加fancy的模型進行嘗試；
• 數據層面：構建一些特征或者對數據進行一些數據變換等；
• 針對損失函數設計各種trick的提升技巧；

創作挑戰賽新人創作獎勵來咯，堅持創作打卡瓜分現金大獎

總結

以上是生活随笔為你收集整理的【算法竞赛学习】AI助力精准气象和海洋预测的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。