Python数据分析案例55——基于LSTM结构自编码器的多变量时间序列异常值监测

案例背景

时间序列的异常值检测是方兴未艾的话题。比如很多单变量的，一条风速，一条用电量这种做时间序列异常值检测，想查看一下哪个时间点的用电量异常。

多变量时间序列由不同变量随时间变化的序列组成，这些时间序列在实际应用中通常
来自不同的传感器或数据源。多变量时间序列异常检测是指在多个相互关联的时间序列中识
别出不符合正常模式的异常点。异常检测的应用场景包括IT 运维监控、交通流量管理、环
境检测、能源管理等。例如在空气质量检测数据中，通过监测PM2.5、二氧化碳等数据的变
化，可以及时发现污染事件，采取相应的环境保护措施。

但是如果变量很多的话，就比如网络上的某些流量他有很多条通道或者是监测的变量不止一个，他们随着时间变化都呈现了不同的趋势，而多变量的异常值检测就也可以用传统的方法来做，比如孤立森林，支持数据描述，异常因子等，使用机器学习的方法无监督进行异常值的标签预测。

但是也可以用深度学习的方法。一般来说都会用自编码器，中间的架构采用多层感知机mlp。但是今天我觉得可以用时间序列里面比较好的lstm循环神经网络去代替多层感知机结构，那么就需要把一个二维数据，也就是行代表不同的时间，列代表不同的变量给它构造为三维的循环神经网络的数据结构，放入神经网络去训练自编码器，然后再去进行解码来预测是否重构误差是否为异常值，这个其中就涉及到数据形状的一个转换的问题。以及在构建滑动窗口的时候会将数据的窗口样本量缺少，而造成前几个点无法预测的问题。今天这些问题都会在这个案例里面一一进行解答和演示。

数据介绍

任务描述：给定一段时期内的服务器性能监测数据，且是随时间记录等间隔的观测数据。
 多变量时间序列异常检测的输入表示为，其中n 表示时间戳的最大长度，k
表示输入特征数量。
 多变量时间序列异常检测任务的输出向量是 ，其中y_{i}表示第i 个时间戳
是否为异常，0 为正常，1 为异常。
数据：本数据集包含 25 个服务器在不同时间点收集的指标数据。每个指标表示服务器的
某种性能度量，例如 CPU 使用率、内存使用率等。数据采样间隔为 5 分钟，每条时间序
列代表一个指标。

train.csv：包含训练数据，数据形状 [时间步数，特征数] ，每个时间步表示一个采样时间
点，每个特征表示一个服务器指标。文件包含 26 列，第 1 列timestamp_(min)表示时间
戳，其余25 列feature_X 表示不同服务器指标数据。训练数据集中不包含异常标签。（因为是无监督）
test.csv：包含测试数据，数据形状 [时间步数，特征数]。其余含义同train.csv 文件说明。

我们需要做的就是预测test里面的全部数据点是否为异常值，给他们打上标签。

当然需要本次案例的数据和全部代码文件的同学还是可以参考：多变量异常值检测。

代码实现

我们采用新手友好的TensorFlow的keras框架，导入包

import os
import pandas as pd
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 pandas.plotting import scatter_matrix
import pickle
import h5py
from scipy import stats
 
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error
import tensorflow as tf
from tensorflow.keras.models import Model, load_model,Sequential
from tensorflow.keras.layers import Input, Dense,LSTM,RepeatVector, BatchNormalization,Dropout
from tensorflow.keras.callbacks import ModelCheckpoint, EarlyStopping
from keras.optimizers import Adam
from tensorflow.keras import regularizers
from tensorflow.keras.utils import plot_model
#sns.set(style='whitegrid', palette='muted', font_scale=1.5)

读取数据

data1=pd.read_csv('train.csv').set_index('timestamp_(min)').iloc[288*300:,:]  # 训练集
data1.head()

由于数据量太大，我这里就只取了288*300之后的数据，288是一天的数据量，我少取了300天的数据，因为还剩下4w多条数据，训练一个小网络是够的。

 描述性统计看看

data1.describe()

可以看到数据基本上都是最大值是一和最小值是零都已经进行了标准化，因此可以直接放入神经网络训练，而不需要进行额外的归一化或者是标准化。

查看前12个特征的时间序列走势图。

data1.iloc[:,:12].plot(figsize=(18,5))

查看后面13个特征的时间序列走势图。

data1.iloc[:,12:].plot(figsize=(18,5))

可以看到数据的范围都是在零跟一之间，并且都是有一定的周期性和季节性的波动。  而有的特征波动性挺大的可能就是异常点。

数据准备

时间序列需要构建滑动窗口特征，因为要构造为三维数据，所以我们会损失几个样本点。但是这样少了几个样本的特征之后，前几个样本的标签就没有办法预测了，因此我们可以用均值进行填充。保证数据的x的长度和y的长度是一模一样的。

window_size=288
def get_traindata(data1,    window_size = window_size):
    # 使用滑动窗口生成时间序列数据
    # 计算均值填充
    mean_values = data1.iloc[:, :].mean()
    padding_data = pd.DataFrame([mean_values] * (window_size - 1))
    
    # 在数据前面添加均值填充的数据
    df_padded = pd.concat([padding_data, data1], ignore_index=True)
    
    # 提取特征数据
    X = df_padded.iloc[:, :].to_numpy()
    
    # 构建滑动窗口
    X_windows = []
    for i in range(len(data1)):
        X_windows.append(X[i:i + window_size])
    X_train = np.array(X_windows) ; X= data1.to_numpy()
    return X_train,X

X_train,X=get_traindata(data1, window_size = 288)
X_train.shape,X.shape

可以看到原来读取的数据是4.6万条，现在还是4.6万条。

测试集也是一样处理

data2=pd.read_csv('test.csv').set_index('timestamp_(min)')#.iloc[:1440,:]
X1_train,X1=get_traindata( data2, window_size = 288)
X1_train.shape,X1.shape

异常值监测

自编码器

我们首先来构造神经网络的架构层

model = Sequential()
model.add(LSTM(16, activation='relu', input_shape=(window_size, X.shape[1]), return_sequences=False))
model.add(BatchNormalization())  # 添加批量归一化层
model.add(Dropout(0.1))  # 添加Dropout层以减少过拟合的风险
model.add(RepeatVector(window_size))
model.add(LSTM(16, activation='relu'))
model.add(BatchNormalization())
model.add(Dropout(0.1))  # 再次添加Dropout层
model.add(Dense(X.shape[1]))
model.compile(optimizer=Adam(learning_rate=0.003), loss='mse')

查看模型信息

model.summary()

因为我安装的tf版本比较高，是不支持GPU的，因此我选择把神经元的数量都改小一点，总共的参数也就5000多个，不怎么多。

可以清楚的看到神经网络的架构主要是数据进入过lstm层，然后再经过层批量归一化，防止梯度消失。然后再经过一个丢包层进行一定的正则化，然后再经过时间的一个向量重复层将它恢复为三维的数据，然后再经过lstm层，再经过成批量归一化，然后再经过一个丢包层，最后用多层感知机进行输出。输出的维度跟x的变量数量是一样的。因为是自编码器嘛，相当于是三维的数据输入，二维的数据输出。相当于是自己用自己进行训练和监督，重构这个误差。

模型训练：

# 训练模型
history =model.fit(X_train, X, epochs=4, batch_size=512, shuffle=False,verbose=1).history

我就训练了四轮，因为每一轮都还有点久的。而且四轮基本上已经收敛了。

产看损失变化图

plt.plot(history['loss'])
#plt.plot(history['val_loss'])
plt.title('model loss')
plt.ylabel('loss')
plt.xlabel('epoch')
plt.legend(['train', 'val'], loc='upper right');

# 计算重构误差

reconstructed_X = model.predict(X_train)
print(reconstructed_X.shape)
mse = np.mean(np.power(X - reconstructed_X, 2), axis=1)
mse

最后得到一个误差列，相当于是每一个样本的误差数值。

我们将其转化为数据框，方便进行存储和后面的筛选

testMSE = mean_squared_error(X.transpose(), reconstructed_X.transpose(), multioutput='raw_values')
error_df = pd.DataFrame({'reconstruction_error': testMSE,'true_value': data1['feature_24']})
print(error_df.shape)
error_df.head()

对误差画个直方图看看。

error_df['reconstruction_error'].plot.hist()

可以看到直方图是一个很明显的右偏分布，它存在很多极大值，这些极大的异常值可能就是所谓的异常点。

我们选择95%作为分位数及95%以下的数据是作为正常值，5%的数据作为异常点，我们去找到这个切割阈值的分位数。

threshold = error_df.reconstruction_error.quantile(q=0.95)
threshold

然后我们按这个阈值标记为正常值和异常值，查看其数量。

mlp_label=np.where( error_df.reconstruction_error.to_numpy()>threshold,1,0)
error_df['pred_class']=mlp_label
error_df['pred_class'].value_counts()

我们对重构误差进行一定的可视化。

groups = error_df.groupby('pred_class')
fig, ax = plt.subplots(figsize=(10,3),dpi=128)
for name, group in groups:
    if name == 1:
        MarkerSize = 3 ; Color = 'orangered' ; Label = 'Fraud' ; Marker = 'd'
    else:
        MarkerSize = 3 ; Color = 'b' ; Label = 'Normal' ; Marker = 'o'
        
    ax.plot(group.index, group.reconstruction_error, linestyle='',color=Color,label=Label,ms=MarkerSize,marker=Marker)
    
ax.hlines(threshold, ax.get_xlim()[0], ax.get_xlim()[1], colors="r", zorder=100, label='Threshold')
ax.legend(loc='upper left', bbox_to_anchor=(0.95, 1))
plt.title("Probabilities of fraud for different classes")
plt.ylabel("Reconstruction error")  ;   plt.xlabel("Data point index")
plt.show()

X轴是代表不同的样本，y轴代表重构的误差，我们把阈值直线也画在了上面，可以清楚的看到重构误差大的一些点，可能就被认为是异常点。

我们把最后一个特征的序列图进行可视化，并在上面标记异常值的点。

def plot_Abnormal(error_df,mode='Autoencoder'):
    plt.figure(figsize=(14, 4),dpi=128)
    plt.plot(error_df.index, error_df["true_value"], label='Power')
    # Plotting the pred_class scatter plot for points where pred_class is 1
    pred_class_1 = error_df[error_df["pred_class"] == 1]
    plt.scatter(pred_class_1.index, pred_class_1["true_value"], color='orange', label='Abnormal value')

    # Adding labels and legend
    plt.xlabel('Time')
    plt.ylabel('Power')
    plt.title(mode)
    plt.legend()
    plt.grid(True)
    plt.show()
plot_Abnormal(error_df,mode='Autoencoder')

可以看到这些较高或者较低的值被标为了黄色的点就是异常值，，由于这是单个特征，所以说不是很明显。只是有的比较异常的位置被标为了异常值，因为异常值的重构误差是所谓的特征一起决定的。但是在图表上只能对于一个特征，一个特征的可视化，这里就只可视化的这一个来观察一下。目前感觉效果是还可以的。

测试集上预测

我们对测试集的数据同样进行重构误差的预测计算，然后将其装入数据框查看。

# 计算重构误差
reconstructed_X1 = model.predict(X1_train)
print(reconstructed_X1.shape)
testMSE1 = mean_squared_error(X1.transpose(), reconstructed_X1.transpose(), multioutput='raw_values')
error_df1 = pd.DataFrame({'reconstruction_error': testMSE1,'true_value': data2['feature_24']})
print(error_df1.shape)
error_df1.head()

error_df1['reconstruction_error'].plot.hist()

可以看到分布跟训练集是一致的，基本上都是一个右偏分布，存在较多的较大的异常值，这些异常值可能就是异常点。

同样我们找到95%的分位的分位数作为阈值。

threshold1 = error_df1.reconstruction_error.quantile(q=0.95)
threshold1

同样我们查看其划分出来的正常值和异常值的数量对比。

mlp_label1=np.where( error_df1.reconstruction_error.to_numpy()>threshold1,1,0)
error_df1['pred_class']=mlp_label1
error_df1['pred_class'].value_counts()

groups1 = error_df1.groupby('pred_class')
fig, ax = plt.subplots(figsize=(10,3),dpi=128)
for name, group in groups1:
    if name == 1:
        MarkerSize = 3 ; Color = 'orangered' ; Label = 'Fraud' ; Marker = 'd'
    else:
        MarkerSize = 3 ; Color = 'b' ; Label = 'Normal' ; Marker = 'o'
        
    ax.plot(group.index, group.reconstruction_error, linestyle='',color=Color,label=Label,ms=MarkerSize,marker=Marker)
    
ax.hlines(threshold1, ax.get_xlim()[0], ax.get_xlim()[1], colors="r", zorder=100, label='Threshold')
ax.legend(loc='upper left', bbox_to_anchor=(0.95, 1))
plt.title("Probabilities of fraud for different classes")
plt.ylabel("Reconstruction error")  ;   plt.xlabel("Data point index")
plt.show()

可以清楚的看到，算出来的很多点误差都还比较高，都是可能潜在的异常值。

同样对最后一个特征的序列图上，把异常点标上进行可视化。

plot_Abnormal(error_df1,mode='Autoencoder')

感觉效果还行，它在某些边缘的极大值，极小值上都标记为了异常点。至于为什么有的地方可能极大值还高一些，它没有标记为异常情况，那可能是因为别的特征一起的作用使得它们的误差并没有那么大。

最后我们将测试集的预测结果储存下来，这样就可以进行提交了。

## 储存
error_df1['pred_class'].to_csv('结果.csv')

总结

以前的自编码器基本上都是使用mlp架构，这里使用了lstm架构完成了一个数据的转换的一个思维。其实也可以进行更高维的转换，比如用二维CNN，或者用同样三维数据结构的一维卷积或者是transformer注意力机制等网络结构。或者是图神经网络等等的方法去，进行不同的这个自编码器架构的构建。

当然由于是无监督的，所以也不好评价哪种方法效果是比较好，只是说这些结构都是可以尝试的。在完成了数据的这种转换之后，神经网络的内部架构都是可以随便更改的，就像搭积木一样是很简单的。也非常容易创新。

创作不易，看官觉得写得还不错的话点个关注和赞吧，本人会持续更新python数据分析领域的代码文章~(需要定制类似的代码可私信)