基于LSTM的多变量多输出温湿度预测

1、模块导入

1. import tensorflow as tf
2. import matplotlib as mpl
3. import matplotlib.pyplot as plt
4. import numpy as np
5. import os
6. import pandas as pd
7. import seaborn as sns
8. from sklearn.model_selection import train_test_split
9. from sklearn.preprocessing import MinMaxScaler
10. from sklearn.metrics import r2_score
11. from tensorflow.keras import utils,losses,layers,Sequential
12. from tensorflow.keras.callbacks import ModelCheckpoint,TensorBoard
15. mpl.rcParams['figure.figsize'] = (10, 8)
16. mpl.rcParams['figure.dpi'] = 150
17. mpl.rcParams['axes.grid'] = False

复制代码

2、加载数据集、预处理

使用 Max Planck Institute for Biogeochemistry 的天气时间序列数据集。 该数据集包含14个不同的特征，例如气温，大气压力和湿度。从2003年开始，每10分钟收集一次。为了提高效率，本文仅使用2009年至2016年之间收集的数据。
下载数据集：

1. zip_path = tf.keras.utils.get_file(
2.     origin='https://storage.googleapis.com/tensorflow/tf-keras-datasets/jena_climate_2009_2016.csv.zip',
3.     fname='jena_climate_2009_2016.csv.zip',
4.     extract=True)
5. csv_path, _ = os.path.splitext(zip_path)

复制代码

加载数据集：

1. df = pd.read_csv(csv_path,parse_dates=['Date Time'],index_col=['Date Time'])
2. df.head()

复制代码

p (mbar)T (degC)Tpot (K)Tdew (degC)rh (%)VPmax (mbar)VPact (mbar)VPdef (mbar)sh (g/kg)H2OC (mmol/mol)rho (g/m**3)wv (m/s)max. wv (m/s)wd (deg)Date Time2009-01-01 00:10:00996.52-8.02265.40-8.9093.33.333.110.221.943.121307.751.031.75152.32009-01-01 00:20:00996.57-8.41265.01-9.2893.43.233.020.211.893.031309.800.721.50136.12009-01-01 00:30:00996.53-8.51264.91-9.3193.93.213.010.201.883.021310.240.190.63171.62009-01-01 00:40:00996.51-8.31265.12-9.0794.23.263.070.191.923.081309.190.340.50198.02009-01-01 00:50:00996.51-8.27265.15-9.0494.13.273.080.191.923.091309.000.320.63214.3如上所示，每10分钟记录一次观测值，一个小时内有6个观测值，一天有144（6x24）个观测值。 给dataset插入新列，列为Data Time列的相应时间值
3、数据可视化

画图看相关性，提取有效特征集
这里选择p (mbar)、Tdew (degC)、max. wv (m/s)作为T (degC)和rh (%)的特征

1. plt.figure(figsize=(16,8))
2. #作图辅助库
3. sns.lineplot(x='p (mbar)',y='T (degC)',data=df[:10000])
4. plt.show()
6. plt.figure(figsize=(16,8))
7. sns.lineplot(x='Tdew (degC)',y='T (degC)',data=df[:10000])
8. plt.show()
10. plt.figure(figsize=(16,8))
11. sns.lineplot(x='max. wv (m/s)',y='T (degC)',data=df[:50000])
12. plt.show()

复制代码

以上是看T (degC)和p (mbar)、Tdew (degC)、max. wv (m/s),的关系
给dataset插入新列，列为Data Time列的相应时间值

1. df['year']=df.index.year
2. df['hour']=df.index.hour
3. df['month']=df.index.month

复制代码

时间与温度的点图

1. plt.figure(figsize=(16,8))
2. sns.pointplot(x='hour',y='T (degC)',data=df[0:50000],hue='month')
3. plt.show()

复制代码

时间与湿度的点图

1. plt.figure(figsize=(16,8))
2. sns.pointplot(x='hour',y='rh (%)',data=df[0:50000],hue='month')
3. plt.show()

复制代码

由于温度与每日的小时变化有关系，而且0-23作为一个循环，所以用三角函数提取周期信息，sin和cos同时使用是因为确保24小时为一个周期

1. df['sin(h)']=[np.sin((x) * (2 * np.pi / 24)) for x in df['hour']]
2. df['cos(h)']=[np.cos((x) * (2 * np.pi / 24)) for x in df['hour']]
3. df

复制代码

4、数据预处理

切分数据集

1. #定义切分函数，x是选取的特征组成的例表，y是标签列（x=dataset[future=] ，y=dataset['T (degC)']）
2. #train_dataset,train_labels=multivariate_data(x_train,y_train,0,100000,3,1,1,True)
3. #上面的一个使用的意思就是：从0开始数到10万，按照3条x数据作为一个元素放入data-》1条y数据作为一个元素存入labels，step=1表示每一条数据就按照上面包装一次，比如data[0]=x[0，1，2]->labels[0]=y[3];data[1]=x[1,2,3]->labels[1]=y[4];
4. #single_step意思是只预测目标的一个未来状态，只预测后1小时，设置为false可以预测未来0到target_size小时内的温度。
5. def multivariate_data(x,y, start_index, end_index, history_size,target_size, step, single_step):
6.     data = []
7.     labels = []
9.     start_index = start_index + history_size
10.    
11.     if end_index is None:
12.         end_index = len(dataset) - target_size
14.     for i in range(start_index, end_index):
15.         indices = range(i-history_size, i, step) # step表示滑动步长
16.         mid_data=x.iloc[indices]
17.         data.append(mid_data)
19.         if single_step:
20.             mid_data=y.iloc[i+target_size]
21.             labels.append(mid_data)
22.         else:
23.             labels.append(y.iloc[i:i+target_size])
25.     return np.array(data), np.array(labels)

复制代码

数据归一化

1. future=['sin(h)','cos(h)','month','max. wv (m/s)','p (mbar)','T (degC)','rh (%)']
2. #数据归一化，由于sin和cos本来就是-1到1，不用归一化
3. for col in future:
4.     scaler=MinMaxScaler()
5.     if(col not in ['sin(h)','cos(h)']):
6.         df[col]=scaler.fit_transform(df[col].values.reshape(-1,1))

复制代码

格式转化与分组和打乱

1. #获取训练特征和训练标签
2. label = ['T (degC)','rh (%)']
3. x=df[future]
4. y=df[label]
6. #通过7-3划分训练集和测试集，70%为训练集 30%为测试集
7. x_train,x_test,y_train,y_test=train_test_split(x,y,test_size=0.3,shuffle=False,random_state=13)
9. #取得训练集，和测试集的格式——》(3,6)->(1,)通过6行历史数据7列目标特征预测1行1列的目标
10. train_dataset,train_labels=multivariate_data(x_train,y_train,0,100000,3,1,1,True)
11. test_dataset,test_labels=multivariate_data(x_test,y_test,0,100000,3,1,1,True)
13. #创建训练组，内部的batch_size，buffer_size，shuffle，batch建议百度
14. #该函数目标是把刚建好的训练集/测试集转化成tensorflow的数据集格式，打乱分组方便训练模型......
15. def create_batch_dataset(x,y,train=True,buffer_size=1000,batch_size=128):
16.     batch_data=tf.data.Dataset.from_tensor_slices((tf.constant(x),tf.constant(y)))
17.     if train:
18.         return batch_data.cache().shuffle(buffer_size).batch(batch_size)
19.     else:
20.         return batch_data.batch(batch_size)
21.         
22. #使用上面函数
23. train_batch_dataset=create_batch_dataset(train_dataset,train_labels)
24. test_batch_dataset=create_batch_dataset(test_dataset,test_labels,train=False)

复制代码

5、模型搭建、编译、训练

1. #建立神经网络模型-3层LSTM和一个输出层
2. model= tf.keras.models.Sequential([
3.    
4.     tf.keras.layers.LSTM(256, input_shape=train_dataset.shape[-2:],return_sequences=True), # input_shape=(20,1) 不包含批处理维度
5.     tf.keras.layers.Dropout(0.4),
6.     tf.keras.layers.LSTM(128, return_sequences=True),
7.     tf.keras.layers.Dropout(0.3),
8.     tf.keras.layers.LSTM(32),
9.     tf.keras.layers.Dense(2)
10. ])
12. #优化器和损失函数设置
13. model.compile(optimizer='adam',loss='mse')
15. #模型保存的相关设置
16. utils.plot_model(model)
17. checkpoint_file='test_model.hdf5'
18. checkpoint_callback=ModelCheckpoint(filepath=checkpoint_file,monitor='loss',moode='min',save_best_only=True,save_weights_only=True)
19. #模型训练
20. history=model.fit(train_batch_dataset,epochs=30,validation_data=test_batch_dataset,callbacks=[checkpoint_callback])

复制代码

通过history获取模型每步训练取得的结果loss和val_loss

1. plt.figure(figsize=(8,8),dpi=200)
2. plt.plot(history.history['loss'])
3. plt.plot(history.history['val_loss'])
4. plt.title('model train vs validation loss')
5. plt.ylabel('loss')
6. plt.xlabel('epoch')
7. plt.legend(['train','validation'], loc='best')
8. plt.show()

复制代码

6、模型验证

1. #通过输入一组数据预测
2. test_preds=model.predict(test_dataset,verbose=1)
3. test_preds[:10]
5. # #将预测后的一组数据转化为1维方便比较
6. TEMP_preds = test_preds[:,0]
7. TEMP_labels = test_labels[:,0]
8. HUM_preds = test_preds[:,1]
9. HUM_labels = test_labels[:,1]

复制代码

温度预测

1. #r2检验
2. score=r2_score(TEMP_labels,TEMP_preds)
3. print(score)

复制代码

0.991836296750627

1. #做出预测结果和实际结果的曲线对比
2. plt.figure(figsize=(16,8))
3. plt.plot(TEMP_labels,label="True value")
4. plt.plot(TEMP_preds,label="Pred value")
5. plt.legend(loc='best')
6. plt.show()

复制代码

湿度预测

1. #r2检验
2. score=r2_score(HUM_labels,HUM_preds)
3. print(score)

复制代码

0.9854786099464197

1. #做出预测结果和实际结果的曲线对比，使用1000次结果对比
2. plt.figure(figsize=(16,8))
3. plt.plot(HUM_labels,label="True value")
4. plt.plot(HUM_preds,label="Pred value")
5. plt.legend(loc='best')
6. plt.show()

复制代码

免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！