作為RNN的一種變體，LSTM廣泛用于時間序列的預(yù)測涣狗。本文結(jié)合EMD（empirical mode decomposition）算法及LSTM提出了EMD-LSTM算法用于空氣質(zhì)量預(yù)測袜匿。結(jié)果表明更啄，僅使用LSTM算法時，預(yù)測結(jié)果具有滯后性居灯，與LSTM相比祭务，EMD-LSTM不僅能夠大幅提高RMSE（root-mean-square-error），且能夠改善預(yù)測值滯后于真實值現(xiàn)象怪嫌。

LSTM簡介

LSTM為RNN的一種變體义锥，能夠解決RNN在訓(xùn)練過程中出現(xiàn)的梯度消失問題，可以記憶長距離的依賴信息岩灭。詳情請點擊該鏈接（https://zybuluo.com/hanbingtao/note/581764）

LSTM預(yù)測

使用LSTM進(jìn)行預(yù)測需將時間序列數(shù)據(jù)轉(zhuǎn)換成監(jiān)督學(xué)習(xí)數(shù)據(jù)拌倍，具體過程可參考鏈接（https://machinelearningmastery.com/multivariate-time-series-forecasting-lstms-keras/），代碼如下：

from keras import Sequential
from keras.layers import LSTM, Dense, Dropout
from keras.regularizers import l2
from Air_Pollution_Forcast_Beijing.model.data_tranform import scaler, test_x, train_X, test_X, train_y, test_y, n_hours, \
    n_features
import matplotlib.pyplot as plt
from numpy import concatenate  # 數(shù)組拼接
from math import sqrt
from sklearn.metrics import mean_squared_error
from scipy.interpolate import spline
import numpy as np

model = Sequential()
model.add(LSTM(20, input_shape=(train_X.shape[1], train_X.shape[2]), return_sequences=True, kernel_regularizer=l2(0.005),
               recurrent_regularizer=l2(0.005)))
model.add(LSTM(20, kernel_regularizer=l2(0.005), recurrent_regularizer=l2(0.005)))
model.add(Dense(1))
model.compile(loss='mae', optimizer='adam')
history = model.fit(train_X, train_y, epochs=100, batch_size=2**8, validation_data=(test_X, test_y))

'''
    對數(shù)據(jù)繪圖
'''
plt.plot(history.history['loss'], label='train')
plt.plot(history.history['val_loss'], label='test')
plt.legend()
plt.show()

# make the prediction,為了在原始數(shù)據(jù)的維度上計算損失噪径，需要將數(shù)據(jù)轉(zhuǎn)化為原來的范圍再計算損失
yHat = model.predict(test_X)
y = model.predict(train_X)
test_X = test_X.reshape((test_X.shape[0], n_hours*n_features))

'''
    這里注意的是保持拼接后的數(shù)組  列數(shù)  需要與之前的保持一致
'''
inv_yHat = concatenate((yHat, test_X[:, -7:]), axis=1)   # 數(shù)組拼接
inv_yHat = scaler.inverse_transform(inv_yHat)
inv_yHat = inv_yHat[:, 0]

test_y = test_y.reshape((len(test_y), 1))
inv_y = concatenate((test_y, test_X[:, -7:]), axis=1)
inv_y = scaler.inverse_transform(inv_y)    # 將標(biāo)準(zhǔn)化的數(shù)據(jù)轉(zhuǎn)化為原來的范圍
inv_y = inv_y[:, 0]


model.summary()
rmse = sqrt(mean_squared_error(inv_yHat, inv_y))
print('Test RMSE: %.3f' % rmse)
raw = inv_y.size
inv_y = inv_y[-24*3:]
inv_yHat = inv_yHat[-24*3:]
plt.plot(inv_yHat, label='forecast')
plt.plot(inv_y, label='observation')
plt.ylabel('pm2.5')
plt.legend()
plt.show()

LSTM預(yù)測結(jié)果

預(yù)測結(jié)果RMSE為22.9柱恤，從圖中可以看到，預(yù)測值滯后于真實值且當(dāng)前時刻的預(yù)測值幾乎等于上一時刻的真實值找爱。這種現(xiàn)象可能是由于時間序列的非平穩(wěn)性導(dǎo)致的梗顺，需要對時間序列進(jìn)行平穩(wěn)性處理。對時間序列進(jìn)行平穩(wěn)性處理的方法包括對數(shù)變換车摄、平滑法寺谤、差分及分解（emd、小波變換等）等方法吮播，本文采用emd分解算法對時間序列進(jìn)行分解得到序列的imf（Intrinsic Mode Function变屁，本征模函數(shù)）及殘差，再對各分量數(shù)據(jù)分別進(jìn)行LSTM預(yù)測薄料，最后將各分量預(yù)測結(jié)果進(jìn)行疊加敞贡，得到最終預(yù)測結(jié)果。

EMD-LSTM預(yù)測

時間序列信號經(jīng)emd分解得到的各分量數(shù)據(jù)如下所示：

emd分解

對各分量數(shù)據(jù)分別進(jìn)行LSTM預(yù)測摄职，代碼如下：

import pandas as pd
import numpy as np
from Air_Pollution_Forcast_Beijing.util import PROCESS_LEVEL1
from sklearn.preprocessing import LabelEncoder
from sklearn.preprocessing import MinMaxScaler
from Air_Pollution_Forcast_Beijing.model.emd import emd
from pyhht.visualization import plot_imfs
from Air_Pollution_Forcast_Beijing.model.series_to_supervised_learning import series_to_supervised
pd.options.display.expand_frame_repr = False
from keras import Sequential
from keras.layers import LSTM, Dense, Dropout
from keras.regularizers import l2
from numpy import concatenate
from sklearn.metrics import mean_squared_error
from math import sqrt
import matplotlib.pyplot as plt

dataset = pd.read_csv(PROCESS_LEVEL1, header=0, index_col=0)
dataset_columns = dataset.columns
values = dataset.values
# values[:, 0] = pd.Series(values[:, 0]).diff(1)
# values = pd.DataFrame(values).dropna().values
'''
計算時間序列的自相關(guān)圖及偏相關(guān)圖
'''
# fig = plt.figure(figsize=(12,8))
# ax1=fig.add_subplot(211)
# fig = plot_acf(values[:, 0], lags=10, ax=ax1)
# ax2 = fig.add_subplot(212)
# fig = plot_pacf(values[:, 0], lags=10, ax=ax2)


# 對第四列（風(fēng)向）數(shù)據(jù)進(jìn)行編碼誊役，也可進(jìn)行 啞編碼處理
encoder = LabelEncoder()
values[:, 4] = encoder.fit_transform(values[:, 4])
values = values.astype('float32')

# 對數(shù)據(jù)進(jìn)行歸一化處理, valeus.shape=(, 8),inversed_transform時也需要8列
scaler = MinMaxScaler(feature_range=(0, 1))
scaled = scaler.fit_transform(values)
# scaleds = []

'''
進(jìn)行emd分解
'''
pollution = values[:, 0]
imfs = emd(pollution)
plot_imfs(pollution, np.array(imfs))
imfsValues = []
for imf in imfs:
    values[:, 0] = imf
    imfsValues.append(values.copy())
inv_yHats = []
inv_ys  = []
for imf in imfsValues:
    scaler = MinMaxScaler(feature_range=(0, 1))
    scaled = scaler.fit_transform(imf)
    # scaleds.append(scaled)
    n_hours = 4
    n_features = 8
    reframed = series_to_supervised(scaled, n_hours, 1)
    values = reframed.values
    n_train_hours = 365 * 24 * 4
    train = values[:n_train_hours, :]
    test = values[n_train_hours:, :]

    # 監(jiān)督學(xué)習(xí)結(jié)果劃分,test_x.shape = (, 8)
    n_obs = n_hours * n_features
    train_x, train_y = train[:, :n_obs], train[:, -n_features]
    test_x, test_y = test[:, :n_obs], test[:, -n_features]

    # 為了在LSTM中應(yīng)用該數(shù)據(jù)，需要將其格式轉(zhuǎn)化為3D format谷市，即[Samples, timesteps, features]
    train_X = train_x.reshape((train_x.shape[0], n_hours, n_features))
    test_X = test_x.reshape((test_x.shape[0], n_hours, n_features))

    model = Sequential()
    model.add(LSTM(20, input_shape=(train_X.shape[1], train_X.shape[2]), return_sequences=True, kernel_regularizer=l2(0.005),
                   recurrent_regularizer=l2(0.005)))
    model.add(LSTM(20, kernel_regularizer=l2(0.005), recurrent_regularizer=l2(0.005)))
    model.add(Dense(1))
    model.compile(loss='mae', optimizer='adam')
    history = model.fit(train_X, train_y, epochs=500, batch_size=2 ** 10, validation_data=(test_X, test_y))

    # make the prediction,為了在原始數(shù)據(jù)的維度上計算損失蛔垢，需要將數(shù)據(jù)轉(zhuǎn)化為原來的范圍再計算損失
    yHat = model.predict(test_X)
    y = model.predict(train_X)
    test_X = test_X.reshape((test_X.shape[0], n_hours * n_features))

    '''
        這里注意的是保持拼接后的數(shù)組  列數(shù)  需要與之前的保持一致
    '''
    inv_yHat = concatenate((yHat, test_X[:, -7:]), axis=1)  # 數(shù)組拼接
    inv_yHat = scaler.inverse_transform(inv_yHat)
    inv_yHat = inv_yHat[:, 0]
    inv_yHats.append(inv_yHat)

    test_y = test_y.reshape((len(test_y), 1))
    inv_y = concatenate((test_y, test_X[:, -7:]), axis=1)
    inv_y = scaler.inverse_transform(inv_y)  # 將標(biāo)準(zhǔn)化的數(shù)據(jù)轉(zhuǎn)化為原來的范圍
    inv_y = inv_y[:, 0]
    inv_ys.append(inv_y)

inv_yHats = np.array(inv_yHats)
inv_yHats = np.sum(inv_yHats, axis=0)
inv_ys = np.array(inv_ys)
inv_ys = np.sum(inv_ys, axis=0)
rmse = sqrt(mean_squared_error(inv_yHats, inv_ys))
print('Test RMSE: %.3f' % rmse)
inv_y = inv_ys[-24*3:]
inv_yHat = inv_yHats[-24*3:]
plt.plot(inv_yHat, label='forecast')
plt.plot(inv_y, label='observation')
plt.ylabel('pm2.5')
plt.legend()
plt.show()

預(yù)測結(jié)果rmse為18.7，高于僅使用LSTM進(jìn)行預(yù)測的情形且滯后現(xiàn)象得到較大程度的改善迫悠。

emd_lstm預(yù)測結(jié)果