4.4 多维聚合
单维度异常发现算法能够得到单个时间序列的报警强度,通过某种方法对多个维度的报警强度进行“聚合”,就能得到多维时间序列的报警强度。我们仍然介绍简单朴素的方法来完成“聚合”,那就是对各维度加权平均。
AWi=sum(Wri**Hri)
其中AWi是第i个时刻的的聚合报警强度,Wri是第i个时刻m维时间序列的单维算法得到的报警强度,Hri是第i个时刻m维时间序列对应的权重。
权重从何而来?
(一)工艺知识
工业生产中存在大量的工艺知识,相互关联的时间序列中,哪个更重要,哪个不重要,在工艺知识中会有所涉及,查阅这些资料或者询问工艺专家,是可以得到相对准确的权重的。
(二)算法计算
当然,工艺知识不可能特别全面,从历史数据中学习权重是更可靠的办法。不过在介绍权重计算方法前先来看一个故事。
上图是二战时盟军返航飞机弹孔分布图,从图中可以看到,这些弹孔分布并不均匀,翅膀上比较多,引擎上比较少。当时军方普遍认为,应该减少装甲总量,然后在受攻击最多的部位增加装甲,这样飞机可以轻一点,但是防护作用不会减弱,因为防御的效率提高了。但是,这些部位需要增加多少装甲,他们并不清楚,于是找到瓦尔德(哥伦比亚大学的统计研究学家),希望得到答案。但是,瓦尔德彻底否定了他们的想法,给出了相反的答案。
瓦尔德认为,需要加装装甲的地方不应该是留有弹孔的地方,反而是没有弹孔的地方,即飞机的引擎。
瓦尔德讲,飞机各部位被击中的概率应该是均等的,但是引擎上的弹孔却比其余部位少,这说明那些被击中引擎的飞机根本没有机会返航。我们看到的数据,都来自成功返航的飞机,这说明即便翅膀被打得千疮百孔,仍能安全返航。
军方马上按照瓦尔德的建议改进了飞机,取得了良好的效果。
这就是“幸存者偏差”。我们不能只考虑看到的数据(返航的飞机),更应该考虑看不到的数据(未返航的飞机)。
为了避免幸存者偏差,多维时间序列权重的分配方法应该遵循这样的原则:历史上经常发生报警且报警强度大的维度(相当于机翼)权重小,不常发生报警且报警强度小的维度(相当于引擎)权重大。
多维时间序列X
X中每一行Xri是某个时刻各个维度时间序列的值,每一列Xcj是一个时间序列。
利用单维度异常发现算法E(…)计算得到各个维度的报警强度矩阵W,我们已经知道W中的元素可能大于1,而且理论上可以是无穷大,为了避免某个很大的数值影响到整体的计算,再考虑到报警强度为1时就已经表示强度很大了,所以将W中大于1的元素强制等于1。
wij=if(wij>1,1,wij)
每个时刻都有一个权重序列,即权重也是个矩阵H:
当前时刻的权重Hri是用之前一个区间的报警强度Wr[-k]i计算得到的,即:
Hri=F(Wr[-k]i)
其中F(…)是根据前面原则计算权重的方法。
现在我们来研究F(…) 的具体计算方法,即如何计算第i个时刻的权重序列Hri。
令Z =W[-k]i
其中zij是Z第i个时刻第j个维度的报警强度。
1. 计算总报警强度SZ:
SZ=[sz1,sz2,…,szm]
szj=sum(Zcj),j∈[1,m]
szj是SZ的第j个维度的总报警强度。
2. 计算第i个时刻的权重Hri:
每个时刻都能得到一个权重序列SZ。
SZ=[sz1,sz2,…,szm]
假设Hri=Y
Y=[y1,y2,…,ym]
其中yj是第j维时间序列的权重。
为了避免“幸存者偏差”,szj越大,yj越小。对SZ中的元素求倒数再做线性变换就能得到权重Y。
SZ中可能所有元素都是0,即所有维度都没有报警,强制令各维度权重相等,SZ中还可能部分元素是0,即szj可能是0,没有倒数,可以强制令它的倒数等于SZ中最大非0元素倒数的2倍,考虑到这些,yj需要调整一下:
其中1/SZ是SZ中非0元素的倒数序列。
有了权重Hri,聚合报警强度AWi就很容易计算了:
AWi=sum(Wri**Hri)
SPL例程
A |
B |
C |
|
1 |
=file(C1).import@tci().(if(~>1,1,~)) |
/第一维报警强度Xc1 |
|
2 |
=file(C2).import@tci().(if(~>1,1,~)) |
/第二维报警强度Xc2 |
|
3 |
1000 |
/学习区间长度k |
|
4 |
[0,0.01,0.1,0.5,1] |
||
5 |
=to(A3) |
||
6 |
=A1(A5) |
/第一维度学习数据 |
|
7 |
=A2(A5) |
/第二维度学习数据 |
|
8 |
=A6.sum() |
/SZ1 |
|
9 |
=A7.sum() |
/SZ2 |
|
10 |
=[A8,A9] |
||
11 |
if A10.sum()==0 |
=Hri=A10.(1/A10.len()) |
/SZ元素全为0的权重Hri |
12 |
else |
=A10.select(~!=0) |
|
13 |
=B12.max(1/~) |
||
14 |
=A10.(if(~==0,B13*2,1/~)) |
/非0元素处理 |
|
15 |
=s=B14.sum(),Hri=B14.(~/s) |
/权重Hri |
|
16 |
=[A1(A3+1),A2(A3+1)] |
/Wri |
|
17 |
=sum(A16**Hri) |
/AWi |
计算结果示例:
两个维度的报警强度图:
图中横轴都是索引,纵轴是报警强度。
图(a)(b)分别是第一维度和第二维度的报警强度走势图。从图中可以看出学习区间([1,1000])内,第二维度相较于第一维度的报警强度次数多且强度大,权重应该小。例程中计算出的权重也符合预期,分别是0.88和0.12。
最后一个点即当前时刻的聚合报警强度是0.065。