鏖兵金匮，校战玉堂；
苍鹰赤雀，铁轴牙樯。
沉白马而誓众，负黄龙而渡江，
海潮迎舰，江萍送王。
[toc]
论文介绍
在上一篇文章中，我们介绍了图像显著性检测中的谱残差SR模型，在本文中，我们借用了视觉显著性检测域的光谱残差模型到我们的异常检测应用中。光谱残差（SR）是一种高效的无监督算法，它在视觉显著性检测任务中表现出卓越的性能和鲁棒性。而时间序列异常检测任务本质上也类似于视觉显著性检测问题：显著性是照片或场景中“突出”的东西，我们会专注于最重要的区域，而时间序列曲线中出现的异常就类似于视觉中的突出部分。因此该论文将谱残差SR模型从视觉显著性检测域迁移到时间序列异常检测，创新地将SR和CNN结合在一起，以提高SR模型的性能。整体与上篇论文基本框架类似，只是在最后接了个CNN用来确定异常阈值。
挑战
在线时序异常检测中面临着三大挑战：
-
无标签。
由于KPI序列的标注成本很高，往往都缺乏有效的标注
-
缺乏通用性。
KPI曲线的形式多种多样，目前没有很好的通用解决方法
-
高效率。
因为需要检测的KPI曲线的数量很多，有很多KPI曲线的更新频率都是分钟级别的，所以有实时性的要求
显著图计算
上文讲到，大量的图片做均值之后的log频谱趋近于一条直线，那么一张图片的log频谱减去平均图片的log频谱就是显著部分，左差之后可以通过反傅里叶变换复原图片得到显著区域。对于时间序列（不论是一维还是多维），同样可以使用SR算法计算序列的显著区域： \(\begin{aligned} A(f)&=\text { Amplitude }(\mathfrak{F}(\mathbf{x}))\\ P(f)&=\operatorname{Phrase}(\mathfrak{F}(\mathbf{x}))\\ L(f)&=\log (A(f)) \\ A L(f)&=h_{q}(f) \cdot L(f) \\ R(f)&=L(f)-A L(f) \\ S(\mathbf{x})&=\left\|\mathfrak{F}^{-1}(\exp (R(f)+i P(f)))\right\| \\ \end{aligned}\) 其中$\mathbf{x}$是输入的序列数据，一般是滑动窗口数据。首先通过傅里叶变换之后计算振幅谱$A(f)$，然后计算相位谱$P(f)$（对于傅立叶变换结果中的复数$x+i*y$，相位角$\theta=arctan \frac yx$），然后对振幅谱做Log得到$L(f)$，$AL(f)$是$L(f)$进行均值滤波之后的结果，$R(f)$就是Spectral Residual谱，再进行一个傅里叶反变换就可以得到显著区域。最终的结果$S(\mathbf{x})$称为显著图Saliency map。
其中$h_q(f)$均值滤波器如下： \(h_{q}(f) =\frac{1}{q^{2}}\left[\begin{array}{ccccc} 1 & 1 & 1 & \ldots & 1 \\ 1 & 1 & 1 & \ldots & 1 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & 1 & 1 & \ldots & 1 \end{array}\right] \\\) 下图是论文中给出的一个原始时间序列和经过转换之后得到的显著图的对比：
在上面显著性图中的新奇点（红色显示）比在原始时间序列中要重要得多。
简单阈值规则
根据显著图可以使用简单的阈值规则来识别异常点： \(O\left(x_{i}\right)=\left\{\begin{array}{ll} 1, & \text { if } \frac{S\left(x_{i}\right)-\overline{S\left(x_{i}\right)}} {S\left(x_{i}\right)} > \tau, \\ 0, & \text { otherwise, } \end{array}\right. \\\) 其中$x_i$表示时间序列中的任意点，$S(x_i)$为对应的显著性图，$\overline{S(x_i)}$是$S(x_i)$中前$z$个点的平均，$\tau$为指定的阈值。
序列预测
对于online异常检测一般都是采用滑动窗口的方式，往往我们需要检测的点（也就是数据流中的最新的点）是位于一段序列的末端，而SR算法在当检测的点位于序列中央的时候效果才会比较好，因此在进行SR计算之前需要对序列进行简单的预测进而延长序列，论文中的预测方法也很直接：
\[\begin{aligned} \bar{g}&=\frac{1}{m} \sum_{i=1}^{m} g\left(x_{n}, x_{n-i}\right) \\ x_{n+1} &=x_{n-m+1}+\bar{g} \cdot m\\ &=x_{n-m+1}+\sum_{i=1}^{m} g\left(x_{n}, x_{n-i}\right) \\ &=x_{n-m+1}+ (x_{n}-x_{n-1}+x_{n}-x_{n-2}+\cdots+x_{n}-x_{n-m}) \\ \end{aligned}\]其中$g$为两个点之间的梯度，$m$超参为预测时需要考虑前面的多少个点，在论文中设置为了$5$，也就是拿当前点与前面$m$个点分别求梯度再做平均得到平均梯度，对于需要预测的下个点即拿前面第m个点乘平均梯度得到。论文发现第一个估计点起着决定性的作用，因此，我们只需将$x_{n+1}$复制$k$次作为预测序列，并将预测序列添加到原始序列的尾部即可。
SR-CNN
由于SR方法是通过简单的手动设置阈值进行分类的，因此论文称可以使用CNN这种更加强大的分类器进行分类。但是CNN分类的话需要有明确的标签，论文的解决方法也很粗暴：通过异常注入的方法来制造伪标签。具体来说就是随机选择时间序列中的几个点，计算注入的异常值来替换原始点，并得到其显著性图。异常点的值由以下公式确定： \(x = (\overline{x}+mean)(1+var)\cdot r+x\) 其中$\bar{x}$是前面几个点的平均，$mean$和$var$是当前窗口内所有点的均值和方差，$r$是一个服从标准正态分布$r\sim \mathcal{N}(0,1)$的随机采样值。
其架构如下：
上面是SR-CNN的一个整体结构，整体结构很简单，两层的一维卷积，两层全连接层，最后输出一个概率，使用Cross entropy作为loss函数，SGD作为优化方法。值得注意的是，论文中提到，实验结果的模型使用了6500万个点进行训练，这6500万个点应该是指内部的数据集。
代码实现
根据前面提到的各过程和对应的公式，SR异常检测的代码实现如下：
import numpy as np
from scipy import stats
def series_filter(values, kernel_size=3):
"""
均值滤波
"""
filter_values = np.cumsum(values, dtype=float)
filter_values[kernel_size:] = filter_values[kernel_size:] - \
filter_values[:-kernel_size]
filter_values[kernel_size:] = filter_values[kernel_size:] / kernel_size
for i in range(1, kernel_size):
filter_values[i] /= i + 1
return filter_values
def extrapolate_next(values):
"""
点预测
"""
last_value = values[-1]
slope = [(last_value - v) / i for (i, v) in enumerate(values[::-1])]
slope[0] = 0
next_values = last_value + np.cumsum(slope)
return next_values
def marge_series(values, extend_num=5, forward=5):
"""
序列扩展预测，使要检测的点位于中央
"""
next_value = extrapolate_next(values)[forward]
extension = [next_value] * extend_num
if isinstance(values, list):
marge_values = values + extension
else:
marge_values = np.append(values, extension)
return marge_values
class Silency(object):
def __init__(self, amp_window_size, series_window_size, score_window_size):
self.amp_window_size = amp_window_size # 均值滤波的kernel size
self.series_window_size = series_window_size # 序列预测的长度和参考长度
self.score_window_size = score_window_size #
def transform_silency_map(self, values):
"""
Transform a time-series into spectral residual, which is method in computer vision.
For example, See https://github.com/uoip/SpectralResidualSaliency.
:param values: a list or numpy array of float values.
:return: silency map and spectral residual
"""
freq = np.fft.fft(values)
mag = np.sqrt(freq.real ** 2 + freq.imag ** 2)
spectral_residual = np.exp(
np.log(mag) - series_filter(np.log(mag), self.amp_window_size))
freq.real = freq.real * spectral_residual / mag
freq.imag = freq.imag * spectral_residual / mag
silency_map = np.fft.ifft(freq)
return silency_map
def transform_spectral_residual(self, values):
silency_map = self.transform_silency_map(values)
spectral_residual = np.sqrt(
silency_map.real ** 2 + silency_map.imag ** 2)
return spectral_residual
def generate_anomaly_score(self, values, type="avg"):
"""
未经\tau过滤的异常分数
"""
extended_series = marge_series(
values, self.series_window_size, self.series_window_size)
mag = self.transform_spectral_residual(extended_series)[: len(values)]
if type == "avg":
ave_filter = series_filter(mag, self.score_window_size)
score = (mag - ave_filter) / ave_filter
elif type == "abs":
ave_filter = series_filter(mag, self.score_window_size)
score = np.abs(mag - ave_filter) / ave_filter
elif type == "chisq":
score = stats.chi2.cdf((mag - np.mean(mag))
** 2 / np.var(mag), df=1)
else:
raise ValueError("No type!")
return score
import matplotlib.pyplot as plt
x = np.linspace(0, 10, 2000)
y = np.sin(x)+np.sin(5*x)+np.cos(x)
y[1450:1460] = 2.1
y[600:603] = 2.2
noise = np.random.uniform(0, 0.3, len(y))
y += noise
s = Silency(3, 10, 10)
saliency_map = s.transform_spectral_residual(y)
fig, ax = plt.subplots(2, 1, figsize=(10, 5), sharex=True)
ax[0].plot(x, y, c='#6495ed', label='Time series')
ax[1].plot(x, saliency_map, c='#6495ed', label='Saliency map')
mask = saliency_map > 0.2
ax[0].plot(x[mask], y[mask], marker='x', c='r', linestyle='')
ax[1].plot(x[mask], saliency_map[mask], marker='x', c='r', linestyle='')
for i in ax:
i.legend()
这里可以看到，正如论文中所说：
However, the SR method works better if the target point locates in the center of the sliding window
直接对序列计算剩余谱时，序列两端的计算会比较异常。
因此论文中会先对序列做预测，以使序列在数据的中间：
saliency_map = s.generate_anomaly_score(y)
fig, ax = plt.subplots(2, 1, figsize=(10, 5), sharex=True)
ax[0].plot(x, y, c='#6495ed', label='Time series')
ax[1].plot(x, saliency_map, c='#6495ed', label='Saliency map')
mask = saliency_map > 2
ax[0].plot(x[mask], y[mask], marker='x', c='r', linestyle='')
ax[1].plot(x[mask], saliency_map[mask], marker='x', c='r', linestyle='')
for i in ax:
i.legend()
一维均值滤波器（滑动平均）
论文中提到三处均值滤波：一是对log谱做均值滤波去掉背景即$h_q(f)$，二是对显著图做均值滤波得到异常得分即$\frac{S\left(x_{i}\right)-\overline{S\left(x_{i}\right)}}{S\left(x_{i}\right)}$，三是异常点注入中的$x = (\overline{x}+mean)(1+var)\cdot r+x$。
对于均值滤波实现，下面看三种方法：
列表循环
def average_list(a, win=3):
return [a[:i+1].mean() if i < win else a[i-win+1:i+1].mean() for i, j in enumerate(a)]
pandas rolling
def average_pd(a):
return a.rolling(window=3, min_periods=1).mean()
快速方法
def average_filter(values, n=3):
if n >= len(values):
n = len(values)
res = np.cumsum(values, dtype=float)
res[n:] = res[n:] - res[:-n]
res[n:] = res[n:] / n
for i in range(1, n):
res[i] /= (i + 1)
return res
算法过程如下：
对比
对于上面三种方法的对比如下：
import pandas as pd
data = np.random.random(10_000)
data_series = pd.Series(data)
a = average_list(data)
b = average_list(data_series)
c = average_filter(data)
assert np.allclose(a, b)
assert np.allclose(a, c)
可以看到，三种方法算出来的结果是一样的。
average_filter([1, 2, 3, 4, 5, 6, 7])
array([1. , 1.5, 2. , 3. , 4. , 5. , 6. ])
经过下面的耗时测试，可以发现第三种优化方法计算的速度大大提高。
%timeit average_list(data)
39.4 ms ± 350 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
%timeit average_pd(data_series)
244 µs ± 3.85 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%timeit average_filter(data)
40.6 µs ± 298 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
预测
对于原始序列做预测，可将当前点置于检测序列的中央，具体的预测方法可参考前面的公式，实现如下：
last_value = values[-1]
slope = [(last_value - v) / i for (i, v) in enumerate(values[::-1])]
slope[0] = 0
next_values = last_value + np.cumsum(slope)
其中values为预测参考的当前点及之前的一段序列，values[-1]为当前点。
SR-CNN
使用SR-CNN为有监督模型，需要先准备数据和注入标签，数据准备过程如下：
异常点注入
对于窗口数据data：
data = normalize(data) # 归一化原始数据
num = np.random.randint(1, self.number) # 随机生成要替换的点的个数
ids = np.random.choice(self.win_siz, num, replace=False) # 随机选择要替换哪些点
lbs = np.zeros(self.win_siz, dtype=np.int64) # 标签
mean = np.mean(data) # 样本均值
dataavg = average_filter(data) # 样本点的local average
var = np.var(data) # 样本方差
for id in ids: # 异常点计算
data[id] += (dataavg[id] + mean) * np.random.randn() * min((1 + var), 10)
lbs[id] = 1
tmp.append([data.tolist(), lbs.tolist()])
网络搭建
网络如下：
import torch
from torch import nn
from torchviz import make_dot
class Anomaly(nn.Module):
def __init__(self, window=10):
self.window = window
super(Anomaly, self).__init__()
self.conv1 = nn.Conv1d(
window, window, kernel_size=1, stride=1, padding=0)
self.conv2 = nn.Conv1d(
window, 2 * window, kernel_size=1, stride=1, padding=0)
self.fc1 = nn.Linear(2 * window, 4 * window)
self.fc2 = nn.Linear(4 * window, window)
self.relu = nn.ReLU(inplace=True)
def forward(self, x):
x = x.view(x.size(0), self.window, 1)
x = self.conv1(x)
x = self.relu(x)
x = self.conv2(x)
x = x.view(x.size(0), -1)
x = self.relu(x)
x = self.fc1(x)
x = self.relu(x)
x = self.fc2(x)
return torch.sigmoid(x)
model = Anomaly(window=10)
print(model)
X = torch.normal(0, 1, (1, 10))
y = model(X)
graph = make_dot(y, params=dict(list(model.named_parameters())))
graph.render('net', format='jpg');
Anomaly(
(conv1): Conv1d(10, 10, kernel_size=(1,), stride=(1,))
(conv2): Conv1d(10, 20, kernel_size=(1,), stride=(1,))
(fc1): Linear(in_features=20, out_features=40, bias=True)
(fc2): Linear(in_features=40, out_features=10, bias=True)
(relu): ReLU(inplace=True)
)
