好久没更新了最近上班也接触了不少新的模型和方法现在先把手上的有价值的案例再做一个总结太简单的案例就不写了。后面再把新东西做个总结。背景本文使用时空卷积神经网络模型进行多个坐标轴上空间站点的风速进行预测主要进行高维特征数据集的预处理整体的建模训练和预测流程然后进行实验对比分析。由于存在多个经纬度坐标轴上的风速数据传统的循环神经网络系列的模型(RNN,LSTM,GRU)无法考虑地理位置空间上的信息。所以需要引入新的模型方法。时空卷积神经网络ConvLSTM是一种结合卷积神经网络CNN和长短期记忆网络LSTM的混合深度学习模型专门用于处理具有时空依赖性的序列数据。其核心思想是通过卷积操作提取空间特征同时利用LSTM的门控机制建模时间动态从而实现对时空数据的端到端学习。本章的多个空间站点风速预测任务中输入数据表示为五维张量样本量时间步长纬度经度特征数其中时间步长代表历史风速序列的长度纬度和经度构成空间网格特征数可包含风速、风向等变量。其结构如图所示数据简介风电数据集选取本文选取一处风力发电场数据。该地区位于****有多个不同经纬度上的坐标轴上的不同的风速观测情况的变量。在时间跨度上本章选取数据集涵盖了从2023年1月1日00:00到2023年3月31日23:00的时间段时间分辨率保持在1小时。2880个样本中包括年月日时信息、10米高度的东向风速分量和北向风速分量、2米高度的气温、总降水量5个特征变量目标变量为风速。家用电脑跑不了多少数据量只能弄个2000多的小样本本文的全部数据和代码获取参考时空卷积原始数据坐标为维度40个和经度40个总共1600个点的风速数据但是原始数据过于庞大本文只选取最靠前的2个维度和2个经度总共4个坐标上相邻的风电场的数据情况因为相邻的风电场可以更好的使得模型中学习到地理空间上的关联信息。经度取值为-0.25和-0.5纬度取值为53.75和54.0。交叉组合为如下4个地理位置[-0.25, 53.75]、[-0.25, 54.0]、[-0.5,53.75]和[-0.5, 54.0]。同样数据集通常划分为训练集、验证集和测试集本章依旧采用70%的数据作为训练集、20%的数据作为验证集用于模型调参、10%的数据测试。数据预处理导入包本文深度学习用的是keras框架。import os,time import datetime import random as rn import numpy as np import pandas as pd import matplotlib.pyplot as plt import seaborn as sns plt.rcParams [font.sans-serif] SimHei #显示中文 plt.rcParams [axes.unicode_minus]False #显示负号 from sklearn.model_selection import train_test_split from sklearn.preprocessing import MinMaxScaler,StandardScaler from sklearn.metrics import mean_absolute_error from sklearn.metrics import mean_squared_error,r2_score #import keras #import keras.backend as K import tensorflow as tf from keras.layers import Layer from keras.models import Model, Sequential from keras.layers import GRU, Dense,Conv1D,Dropout,Flatten,SimpleRNN,LSTM, ConvLSTM2D ,Reshape, TimeDistributed #MaxPooling1D,GlobalMaxPooling1D,Embedding #from keras.callbacks import EarlyStopping #from tensorflow.keras import regularizers #from keras.utils.np_utils import to_categorical from tensorflow.keras import optimizers读取数据dfpd.read_csv(output__2023_part1.csv,parse_dates[valid_time]).set_index(valid_time).iloc[:2880,:] df.shape展示前5行df.head()这个数据6400列就是不同地理坐标上的各种特征月日时信息、10米高度的东向风速分量和北向风速分量、2米高度的气温、总降水量5个特征变量很混乱所以这个数据需要做预处理整理为我们需要的样式。## 1. 变量名称分开 split_cols df.columns.str.extract(r\(([^,]),([^)])\)_([^_])).set_axis([lat, lon, var],axis1) # 步骤2创建多层列索引 multi_cols pd.MultiIndex.from_arrays([ split_cols[lat], split_cols[lon], split_cols[var]], names[Latitude, Longitude, Variable]) df.columns multi_cols# 步骤3堆叠地理坐标到行索引df df.stack([Latitude, Longitude]).reset_index() df.head(2)# 步骤4设置最终索引结构df df.set_index([valid_time, Latitude, Longitude]).sort_index() df.index.names [Time, Latitude, Longitude]## 5.计算风速df[wind_speed]np.sqrt(df[u10]**2 df[v10]**2)# 查看原始坐标轴 df.index.get_level_values(1).unique() , df.index.get_level_values(2).unique()我们只选取4个点位也就是横2个点。纵坐标2个点2*2组合4个地理位置。太多了计算会很慢## 只采样2 个点 lat_use_listdf.index.get_level_values(1).unique().to_numpy()[::20] lon_use_listdf.index.get_level_values(2).unique().to_numpy()[::20] lat_dimlen(lat_use_list) ; lon_dimlen(lon_use_list) lat_use_list,lon_use_list## 过滤掉其他坐标 df1df.reset_index() df1df1[df1[Latitude].isin(lat_use_list )] df1df1[df1[Longitude].isin(lon_use_list )] df1.shape### 预处理格式完成df1df1.set_index([Time, Latitude, Longitude]) df1前面的时间经纬度设置为了索引后面5列是我们的变量。做一下索引交换经度在前维度在后df_winddf1.unstack([Latitude, Longitude]).swaplevel(axis1).sort_index(axis1, level[1, 2], ascending[True, True])[wind_speed] df_wind.columns划分训练集验证集测试集train_ratio 0.7 val_ratio 0.2 n_samples df_wind.shape[0] train_end int(n_samples * train_ratio) val_end train_end int(n_samples * val_ratio)画图展示# 定义颜色和标签 colors {train: gold, val: lightblue, test: lightpink} labels {train: Train, val: Validation, test: Test} # 创建2x2子图 plt.figure(figsize(16, 10),dpi128) for i, ((lon, lat), col) in enumerate(df_wind.items()): plt.subplot(2, 2, i1) # 分割数据 train df_wind.iloc[:train_end, i] val df_wind.iloc[train_end:val_end, i] test df_wind.iloc[val_end:, i] # 绘制三条线 plt.plot(train.index, train, colorcolors[train], labellabels[train]) plt.plot(val.index, val, colorcolors[val], labellabels[val]) plt.plot(test.index, test, colorcolors[test], labellabels[test]) # 设置标题和图例 plt.title(fLon: {lon}°, Lat: {lat}°, fontsize16) #if i 0: # 只在第一个子图显示图例 plt.legend() plt.tight_layout() plt.show()由于构建的是时空卷积神经网络的模型数据的输入形状要构建为五维张量其维度的含义分别是样本量时间步长纬度经度特征数。本章节数据依旧采用滑动窗口法进行构建数据集。选用滑动窗口24过去1天的风速数据加上其他特征一起预测下一时刻的风速。整体而言对于单个样本Xxt-24,xt-23,xt-22,…,xt-1,xt去预测Yxt1的数据。其中X为四维数组单个xt为三维数组表示为该t时刻的在不同经度和不同维度上的特征向量单个Yxt1为向量表示不同坐标轴上的t1时刻的风速按照如上的方法进行窗口滑动构建样本最终形成的特征变量X是五维张量维度分别是样本量时间步长时间步长纬度经度特征数。Y是二维矩阵。数据标准化datadf1.to_numpy() scaler MinMaxScaler() scaler scaler.fit(data[:,:-1]) Xscaler.transform(data[:,:-1]) y_scaler MinMaxScaler() y_scaler y_scaler.fit(data[:,-1].reshape(-1,1)) yy_scaler.transform(data[:,-1].reshape(-1,1)) X.shape,y.shapex和y标准化好了合并df_datapd.DataFrame(np.c_[X,y], indexdf1.index, columnsdf1.columns, )我们要把这个表格数据处理为5维的张量## 转为5纬数据 def create_convlstm_dataset(df, window_size6, target_colwind_speed): # 步骤1: 获取空间坐标的唯一值按原始顺序 #unique_coords df.index.get_level_values([Latitude, Longitude]).unique() unique_lats df.index.get_level_values(Latitude).unique().sort_values() unique_lons df.index.get_level_values(Longitude).unique().sort_values() # 步骤2: 将数据重塑为4D张量 (时间步, 纬度, 经度, 特征) # 先 unstack 空间坐标到列 tensor_4d ( df.unstack([Latitude, Longitude]) .swaplevel(axis1).sort_index(axis1, level[1, 2], ascending[True, True]) .values.reshape( len(df.index.unique(Time)), len(unique_lats), len(unique_lons), len(df.columns) ) ) # 步骤3: 创建滑动窗口序列 samples [] ; targets [] for i in range(len(tensor_4d) - window_size): # 输入窗口: [i, iwindow_size) samples.append(tensor_4d[i:iwindow_size]) # 目标值: 窗口后一时刻的全空间风速 target_time_idx i window_size target tensor_4d[target_time_idx, :, :, df.columns.get_loc(target_col)] targets.append(target.flatten()) # 展平为 (lat * lon,) # 转换为 numpy 数组 X np.array(samples) # (n_samples, window_size, lat, lon, features) y np.array(targets) # (n_samples, lat * lon) return X, y window_size 24 # 时间窗口大小 X, y create_convlstm_dataset(df_data, window_sizewindow_size) # 打印维度 print(f输入维度: {X.shape} → (样本数, 时间步长, 纬度数, 经度数, 特征数)) print(f输出维度: {y.shape} → (样本数, 纬度*经度))时间步长24是我们设定的2,2是横纵两个维度也就是空间上的经纬度。5就是特征数量y自己也是可以拿来当特征的用自己的前一秒的数值预测下一秒的数值并不算穿越泄露的问题。y就是4个点也就是下一个时刻4个坐标上的每个点位的风速。# 按时间顺序划分 训练集验证集测试集# 按时间顺序划分 训练集验证集测试集 train_ratio 0.7 val_ratio 0.2 n_samples X.shape[0] train_end int(n_samples * train_ratio) val_end train_end int(n_samples * val_ratio) X_train, y_train X[:train_end], y[:train_end] X_val, y_val X[train_end:val_end], y[train_end:val_end] X_test, y_test X[val_end:], y[val_end:] # 转换为 float32 节省内存 X_train X_train.astype(float32) y_train y_train.astype(float32) print(训练集形状:, X_train.shape, y_train.shape) print(验证集形状:, X_val.shape, y_val.shape) print(测试集形状:, X_test.shape, y_test.shape)定义随机数种子评估函数。def set_my_seed(): os.environ[PYTHONHASHSEED] 0 np.random.seed(1) rn.seed(12345) tf.random.set_seed(123) def evaluation(y_test, y_predict): mae mean_absolute_error(y_test, y_predict) mse mean_squared_error(y_test, y_predict) rmse np.sqrt(mse) mape(abs(y_predict -y_test)/ y_test).mean() r_2r2_score(y_test, y_predict) return mae,rmse, mape ,r_2构建模型就是四个模型MLP,LSTM,GRU,ConvLSTM2Ddef build_model(X_train, modeLSTM, hidden_dim[32, 16], lat_dim2, lon_dim2): # 自动计算输入形状 timesteps, height, width, channels X_train.shape[1:] input_shape (timesteps, height, width, channels) model Sequential() if mode MLP: # 直接展平所有维度 model.add(Flatten(input_shapeinput_shape)) model.add(Dense(hidden_dim[0], activationrelu)) model.add(Dense(hidden_dim[1], activationrelu)) model.add(Dense(lat_dim * lon_dim, activationlinear)) elif mode in [LSTM, GRU]: # 先用 MLP 编码每个时间步的空间特征 model.add(TimeDistributed(Flatten(), input_shapeinput_shape)) # (timesteps, height*width*channels) model.add(TimeDistributed(Dense(hidden_dim[0], activationrelu))) # 空间特征编码 # 再输入时序层 if mode LSTM: model.add(LSTM(hidden_dim[1], return_sequencesFalse)) else: model.add(GRU(hidden_dim[1], return_sequencesFalse)) model.add(Dense(lat_dim * lon_dim, activationlinear)) elif mode ConvLSTM2D: # 原始 ConvLSTM2D 处理 model.add(ConvLSTM2D( filtershidden_dim[0], kernel_size(3, 3), paddingsame, activationtanh, input_shapeinput_shape, return_sequencesFalse)) model.add(Flatten()) model.add(Dense(lat_dim * lon_dim, activationlinear)) else: raise ValueError(无效模型名称. 只能选择从 MLP, LSTM, GRU, ConvLSTM2D) model.compile(optimizeradam, lossmse,metrics[tf.keras.metrics.RootMeanSquaredError(),mape,mae] ) return model画图函数def plot_loss(hist, imfname): plt.figure(figsize(16,2),dpi100) keys [k for k in hist.history.keys() if not k.startswith(val_)] for i, key in enumerate(keys): plt.subplot(1, len(keys), i1) plt.plot(hist.history[key], k-, labelfTraining {key}) if fval_{key} in hist.history: plt.plot(hist.history[fval_{key}], r--, labelfValidation {key}) plt.title(f{imfname} {key.capitalize()}) plt.xlabel(Epochs) plt.ylabel(key) plt.legend() plt.tight_layout() plt.show() def plot_fit(y_test, y_pred, col): plt.figure(figsize(6,2)) plt.plot(y_test, colorred, labelactual) plt.plot(y_pred, colorblue, labelpredict) plt.title(f{col}坐标上的 拟合值和真实值对比,fontsize10) plt.xlabel(Time) plt.ylabel(wind) plt.xticks(fontsize8, colork ) plt.legend() plt.show()评估结果的时候我们需要把标准化的数据逆转回去def inverse_transform_and_to_dataframe(y_pred, y_scaler): 对预测结果进行逆标准化并转换为DataFrame 参数: y_pred (np.ndarray) y_scaler (StandardScaler): 已拟合的标准化器 返回: pd.DataFrame: 包含逆变换结果的DataFrame y_pred np.asarray(y_pred) # 确保输入是numpy数组 y_inversed np.zeros_like(y_pred) # 对每一列进行逆标准化变换 for col in range(y_pred.shape[1]): # 需要reshape为2D数组 (n_samples, 1) y_inversed[:, col] y_scaler.inverse_transform( y_pred[:, col].reshape(-1, 1) ).flatten() # 转换为DataFrame df_y_predpd.DataFrame(y_inversed ,columnsdf_sample.columns ,indexdf_sample.index) cols df_y_pred.columns.map(lambda x: ,.join( if Unnamed in i else i for i in x)) cols[f({c}) for c in cols ] df_y_pred.columnscols return df_y_pred训练评估函数这个函数包括训练评估可视化一体的model_list [MLP,LSTM,GRU, ConvLSTM2D] location [f({lon},{lat}) for lon in lon_use_list for lat in lat_use_list] index1 pd.MultiIndex.from_product([model_list, location], names[Model, location]) df_preds_all pd.DataFrame(columnsindex1) df_eval_allpd.DataFrame(columns[MAE,RMSE,MAPE,R2],indexindex1) df_sampledf1.unstack([Latitude, Longitude]).swaplevel(axis1).sort_index(axis1, level[1, 2], ascending[True, True])[wind_speed].iloc[-len(y_test):,:] ### 训练函数 def train_fun(modeLSTM,batch_size32,epochs50,hidden_dim[32,16],lat_dim2,lon_dim2,verbose0,show_lossTrue,show_fitTrue): #构建 和训练模型 s time.time() set_my_seed() modelbuild_model(X_train, modemode, hidden_dimhidden_dim, lat_dimlat_dim, lon_dimlon_dim) #earlystop EarlyStopping(monitorloss, min_delta0, patience5) hist model.fit( X_train, y_train, epochsepochs, batch_sizebatch_size, validation_data(X_val, y_val), verboseverbose) if show_loss: plot_loss(hist) #预测 y_predmodel.predict(X_test) df_y_pred inverse_transform_and_to_dataframe(y_pred, y_scaler) df_y_test inverse_transform_and_to_dataframe(y_test, y_scaler) #数据转化 df_y_testdf_y_test[df_y_pred.columns] print(数据表格的列名称是否相同为True则可以继续不然预测值和真实值可能不是一个坐标上的) print( locationdf_y_test.columns ) #print(f真实y的形状{df_y_test.shape},预测y的形状{df_y_pred.shape}) etime.time() print(f运行时间为{round(e-s,3)}) # 查看预测效果 if show_fit: for col in location: plot_fit(df_y_test[col], df_y_pred[col],col) #储存预测结果 和评价指标 df_preds_all.loc[:,(mode,location)]np.array(df_y_pred) for col in location: scorelist (evaluation(df_y_test[col], df_y_pred[col])) df_eval_all.loc[(mode,col),:]score s[round(i,3) for i in score] print(f{mode}在{col}坐标上的预测效果为MAE:{s[0]},RMSE:{s[1]},MAPE:{s[2]},R2:{s[3]}) print(运行结束) return s[-1]参数设置window_size24 batch_size64 epochs20 hidden_dim[32,16] verbose0 show_fitTrue show_lossTrue modeLSTM #MLP,GRU先训练MLPtrain_fun(modeMLP,batch_sizebatch_size,epochsepochs,hidden_dimhidden_dim, lat_dimlat_dim,lon_dimlon_dim,verbose0,show_lossTrue,show_fitTrue)可以展示所有的坐标点上的拟合效果还有整体的评估指标。下面训练lstm优于上面的训练函数都定义好了改个参数就行。train_fun(modeLSTM,batch_sizebatch_size,epochsepochs,hidden_dimhidden_dim, lat_dimlat_dim,lon_dimlon_dim,verbose0,show_lossTrue,show_fitTrue)GRU模型train_fun(modeGRU,batch_size90,epochsepochs,hidden_dimhidden_dim, lat_dimlat_dim,lon_dimlon_dim,verbose0,show_lossTrue,show_fitTrue)二维卷积模型train_fun(modeConvLSTM2D,batch_sizebatch_size,epochsepochs,hidden_dimhidden_dim, lat_dimlat_dim,lon_dimlon_dim,verbose1,show_lossTrue,show_fitTrue)运行时间是LSTM,GRU的30多倍。。。虽然效果更好一点。调整一下预测结果的数据框的索引下面进行评估。df_preds_all.indexdf_sample.index结果对比可视化评价指标的表df_eval_all可以清楚的看到每个模型在每个坐标点的风速预测的误差指标。预测结果的表df_preds_all.head(5)这些表都是多层索引的。列上面的第一层是模型第二层是不同坐标轴。获取真实值的表df_y_testdf_sample.copy() cols df_y_test.columns.map(lambda x: ,.join( if Unnamed in i else i for i in x)) cols[f({c}) for c in cols ] df_y_test.columnscols df_y_test预测表交换一下索引坐标作为第一层模型作为第二层。df_multistep_predsdf_preds_all.T.swaplevel(0,1).sort_index().T df_multistep_preds.head(3)方便可视化下面就可以把不同坐标上的不同模型预测效果进行可视化Forecasting_horizon[ f坐标{i} for i in location] timeslocation colors3 [ tomato, green, darkviolet, darkorange, royalblue, gold, deepskyblue, crimson, lime, orchid, sienna, cyan, magenta, yellowgreen, slateblue, chocolate, teal, firebrick, dodgerblue, olive ] plt.subplots(2,2,figsize(12,7),dpi256) for i1,ti in enumerate(times): nint(str(22)str(i11)) plt.subplot(n) data_modeldf_multistep_preds[ti] for i2,col in enumerate(data_model.columns): plt.plot(data_model.index,data_model[col],labelcol,colorcolors3[i2]) plt.plot(data_model.index,df_y_test[ti],labelActual,colork,linestyle:,lw2) plt.ylabel(Wind,fontsize12) plt.xlabel(Time,fontsize12) plt.title(f{Forecasting_horizon[i1]},fontsize16) plt.legend(fontsize9,locupper right) plt.tight_layout() plt.show()评价指标表做一下预处理columns pd.MultiIndex.from_tuples(df_eval_all.index, names[Model, location]) df_plot pd.DataFrame(df_eval_all.values, columns[MAE, RMSE, MAPE, R2], indexcolumns).T # 使用链式方法转换数据框 df_plotdf_plot.stack(level1).swaplevel(0,1).sort_index(axis1)df_plot.swaplevel(0,1).sort_index().loc[R2,:].loc[location,:]可视化bar_width 0.09 ; marks [**,xxx,,o,-] Forecasting_horizon[ f坐标{i} for i in location] timeslocation m np.arange(len(times)) plt.subplots(4,1,figsize(5,5),dpi256) for i1,ti in enumerate([MAE,RMSE,MAPE,R2]): nint(str(41)str(i11)) plt.subplot(n) df_one_modeldf_plot.swaplevel(0,1).sort_index().loc[ti,:].loc[location,:] for i,col in enumerate(df_one_model.columns): plt.bar(xmbar_width*i, heightdf_one_model[col], labelcol, widthbar_width,colorcolors3[i])#hatchmarks[i], names[f {index} for index in Forecasting_horizon] plt.xticks(range(0, len(names)),names,fontsize7) plt.ylabel(ti,fontsize10) if i10: plt.legend(fontsize8,bbox_to_anchor(0.9,1.35),ncol4) #,ncollen(df_one_model.columns) plt.tight_layout() plt.show()从上面的预测对比图和评价指标的图可以看到二维时空卷积可以做到比LSTM和GRU更好的效果MLP效果最差。这种二维时空卷积不仅可以放在这种风速预测场景用还可以放在手机设备反欺诈场景五个维度例如样本量时间步长手机x轴手机y轴特征[触摸压力触摸半径]等等