劉 云，王梓宇

(昆明理工大學(xué)信息工程與自動化學(xué)院，云南 昆明 650500)

1 引言

在許多數(shù)據(jù)挖掘任務(wù)[1 - 3]所用的原始數(shù)據(jù)中，或多或少都會存在一些數(shù)據(jù)異常的部分，這些異常數(shù)據(jù)可能是由噪聲或觀測誤差等因素造成的無價值干擾，也可能蘊含著由罕見事件或復(fù)雜過程引起的有價值信息[4 - 6]，檢測并提取這些異常部分的方法在多個領(lǐng)域有著廣闊的應(yīng)用前景。然而大多數(shù)現(xiàn)有的逐點分析異常檢測技術(shù)不能準(zhǔn)確地反映數(shù)據(jù)的異常特征，部分異常區(qū)間檢測技術(shù)，要么不適用于多變量時空數(shù)據(jù)，要么不能較好地檢測大小變化的異常區(qū)間[7,8]。因此，有效檢測多變量時空序列的異常區(qū)間是一個有價值的研究方向。

Keogh等[9]提出了一種使用符號聚合近似的啟發(fā)式有序時間序列HOT SAX(Heuristically Ordered Time series using Symbolic Aggregate approXimation)算法。HOT SAX算法通過將長度為d的子序列表示為一個d維向量，計算該向量到其d維空間中最近鄰的歐氏距離作為異常得分來提取時間序列中的異常子序列。該算法將檢測序列的樣本數(shù)(以下統(tǒng)稱序列長度)事先指定為一個固定參數(shù)，由于時空序列中異常區(qū)間的長度通常不為定值，該固定參數(shù)的設(shè)定將影響算法的檢測精度。Kim等[10]提出了一種魯棒核密度估計RKDE(Robust Kernel Density Estimation)算法。將基于徑向半正定核的核密度估計器KDE(Kernel Density Estimation)解釋為相關(guān)再生核希爾伯特空間中的樣本均值，用經(jīng)典M估計來估計這些均值，可以得到加權(quán)核密度估計器(RKDE)，通過收斂于RKDE的核化迭代加權(quán)最小二乘IRWLS(Iteratively Re-Weighted Least Square)算法計算權(quán)值,從而使RKDE能穩(wěn)健地估計異常值。雖然KDE是一個靈活的模型，但其不考慮數(shù)據(jù)屬性間的相關(guān)性且不能很好地擴展到長時間序列。

為了準(zhǔn)確地找到多變量時空數(shù)據(jù)中的異常區(qū)間，本文提出了一種無偏KL散度UKLD(Unbiased KL Divergence)算法。定義時空數(shù)據(jù)中的異常區(qū)間后，算法首先為時間序列嵌入時間延遲，考慮樣本之間的相關(guān)性；隨后將檢測區(qū)間和剩余區(qū)間估計為高斯分布，并通過累計和來加快參數(shù)估計過程；最后使用無偏KL散度計算分布之間的差異水平從而避免由區(qū)間長度產(chǎn)生的計算偏差。仿真分析結(jié)果表明，與HOT SAX算法和RKDE算法相比，UKLD算法在精度指標(biāo)方面有優(yōu)化提升。

2 模型

2.1 定義

[l,r)表示離散區(qū)間{t∈T|l≤t

r≤n+1∧a≤r-l≤b}

(1)

于是，時間序列X中所有符合給定約束大小A=(at,ax,ay,az),B=(bt,bx,by,bz)的子塊集合定義為：

(2)

其中，It、Ix、Iy和Iz分別表示時間點t和空間位置坐標(biāo)x、y、z的符合約束大小(at,bt)和(ax,bx)、(ay,by)、(az,bz)的區(qū)間。

以下將JA,B省略索引簡記為J。

給定任何區(qū)間I∈J，除了該特定范圍外的時間序列的剩余部分定義為：

(3)

以下簡記為Ω。

定義的時間序列模型如圖1所示。道格拉斯·霍金斯(Douglas Hawkins)把異常定義為“一個與其他觀察結(jié)果偏離巨大的觀察，以至于懷疑它由一個不同的機制產(chǎn)生”[11]。受該定義啟發(fā)，相較于該時間序列的剩余部分Ω，時間序列X的一個子塊I∈J由一個不同機制產(chǎn)生。于是，子塊I稱為時空時間序列X的一個異常區(qū)間，Ω為剩余區(qū)間。I與Ω的差異由差異水平度量D(pI,pΩ)定義。pI和pΩ分別為異常區(qū)間I和剩余區(qū)間Ω的概率分布。

Figure 1 An illustrative example of notations on time-series data圖1 時間序列數(shù)據(jù)符號的說明性示例

2.2 KL散度

散度KL(Kullback-Leibler)[12,13]是衡量2個概率分布之間差異水平的度量指標(biāo)。2個分布的差異越大，KL散度值越大。區(qū)間I與Ω的分布pI、pΩ的KL散度定義如下：

(4)

其中，pI(Xi)、pΩ(Xi)分別表示時間序列x中區(qū)間I與Ω的樣本Xi中區(qū)間I與Ω的分布。

3 無偏KL散度(UKLD)算法

3.1 預(yù)處理

(5)

從而將先前時間步長中的上下文合并到時間序列的每個樣本中。其中嵌入維度κ是堆疊的樣本數(shù)量，時間延遲τ是被列為上下文的2個連續(xù)時間步長間的間隔。

為了靈活有效地估計嵌入時延的時間序列中異常區(qū)間I和剩余區(qū)間Ω的分布，在考慮數(shù)據(jù)屬性間相關(guān)性的同時能擴展到長時間序列，假設(shè)區(qū)間I和Ω中的數(shù)據(jù)都服從多元正態(tài)高斯分布N(μI,SI)和N(μΩ,SΩ)[15]。

高斯分布模型需要對相應(yīng)區(qū)間中所有樣本求和來估計分布參數(shù)[15]。由于其中一些區(qū)間相互重疊，對多個區(qū)間執(zhí)行求和是多余的，并且對大規(guī)模數(shù)據(jù)而言由于計算量過于巨大，這顯然是不可行的。但是，使用累計和[16]后可以顯著加快計算過程。

(6)

3.2 無偏KL散度

本文把無偏KL散度作為分布pI和pΩ的差異水平的度量，這是無偏KL散度算法UKLD的核心。首先分析KL散度DKL(pI,pΩ)的非無偏性，隨后介紹無偏KL散度DU-KL(pI,pΩ)的理論推導(dǎo)。

(7)

由于時間序列的所有樣本都服從正態(tài)分布，其均值也服從正態(tài)分布：

(8)

(9)

所以，平均向量μ的所有維度都是獨立相同的正態(tài)分布變量，其差異也服從正態(tài)分布：

(10)

DKL(pI,pΩ)=

(11)

其中，(μΩ-μI)i表示取向量μΩ-μI的第i維。

理論分析表明，KL散度不是無偏的，其系統(tǒng)地偏向于較小長度的區(qū)間。由于這種偏差與數(shù)據(jù)本身無關(guān)，與區(qū)間長度相關(guān)，較短區(qū)間將自動比較長區(qū)間獲得更高的得分，所以異常值將被劃分為多個連續(xù)的小區(qū)間，這明顯降低了區(qū)間檢測的質(zhì)量。

當(dāng)時間序列的長度n非常大時，卡方分布的尺度接近于：

(12)

(13)

漸進(jìn)服從自由度為d+(d(d+1))/2的卡方分布。對于KL散度，樣本集是異常區(qū)間I中的數(shù)據(jù)，用于檢驗該數(shù)據(jù)分布的參數(shù)由剩余區(qū)間Ω中的數(shù)據(jù)估計而來。于是檢驗統(tǒng)計量可用來度量檢驗區(qū)間I中的數(shù)據(jù)與基于時間序列剩余部分?jǐn)?shù)據(jù)而建立的模型的擬合程度。

由式(11)定義的換算系數(shù)對KL散度歸一化后，檢驗統(tǒng)計量λ與KL散度的關(guān)系為：

λ=2m·DKL(pI,pΩ)

(14)

因此，得到無偏KL散度為：

(15)

該無偏散度由屬性的數(shù)量d影響而不再是區(qū)間長度m，這能顯著提高算法檢測長時間序列中的長度變化異常區(qū)間的準(zhǔn)確性。

得到所有區(qū)間的異常得分并按降序排序后，使用非極大值抑制方法以僅得到非重疊區(qū)間。對于樣本超過10 000的長時間序列，使用近似非極大值抑制方法，通過維持一個固定大小的當(dāng)前最高得分的非重疊區(qū)間列表，從而避免存儲所有區(qū)間的得分。算法最后返回一個區(qū)間排名，輸出前k個異常區(qū)間。

3.3 算法流程

UKLD算法的流程如算法1中所示。首先為時間序列嵌入時間延遲，使用高斯分布估計檢測區(qū)間I和剩余區(qū)間Ω的分布并通過累計和加快高斯分布的參數(shù)估計過程(步驟1～步驟3)；隨后使用無偏KL散度計算2個區(qū)間分布之間的差異水平，將該差異水平作為檢測區(qū)間的異常得分(步驟4～步驟5)。最后使用非極大值抑制方法得到非重疊的異常區(qū)間(步驟6～步驟7)。

算法1UKLD算法

輸出：非重疊異常區(qū)間。

步驟3pI(Xi)=N(μI,SI)，pΩ(Xi)=N(μΩ,SΩ);

步驟4概率分布=?;

步驟5對于?i∈{1,2,…,n}:

步驟6用非極大值抑制方法確定非重疊異常區(qū)間I的集合；

步驟7返回異常區(qū)間I的集合。

3.4 算法的復(fù)雜度分析

4 仿真分析

4.1 數(shù)據(jù)集

與其他類似算法所采用的標(biāo)準(zhǔn)數(shù)據(jù)集一致，本文依據(jù)2017年1月1日至2017年9月30日間的丹麥歐登塞市[18]交通流計數(shù)來評估算法的性能。數(shù)據(jù)中的時間序列從高斯過程GP(h,K)中采樣，h為零均值函數(shù)h(x) =0，K為平方指數(shù)協(xié)方差函數(shù)：

K(xt,xt′)=(2πl(wèi)2)-1/2·

(16)

GP的長度尺度為：l2=0.01，噪聲參數(shù)為：σ2=0.001，δ(t,t′)為克羅內(nèi)克δ函數(shù)。時間序列中異常區(qū)間的長度在時間序列長度的5%～20%變化，異常區(qū)間中的值是從高斯過程中采樣的另一個函數(shù)值。測試數(shù)據(jù)集共有100個時間序列和190個異常。

仿真分析方法包括示例分析和量化分析，其中示例分析在截取一段代表性時間序列上比較3種算法提取的異常區(qū)間，通過對比該段時間序列上3種算法提取的同一異常區(qū)間的特點作判斷。用方框表示異常區(qū)間，提取的異常區(qū)間越完整越精煉說明算法的精度越高。隨后在數(shù)據(jù)集上做量化分析，在整個數(shù)據(jù)集上定量比較3種算法的精度，由于時空數(shù)據(jù)中異常區(qū)間檢測是檢測任務(wù)而非分類任務(wù)，把由交并比[19]量化定義的平均精度AP(Average Precision)作為精度指標(biāo)。

4.2 固定序列長度的精度分析

從圖2中可以看出，檢測序列長度為125時，HOT SAX、RKDE和UKLD 3種算法都能找到該段時間序列中的3個異常區(qū)間。圖2中，f(x)表示由高斯過程采樣后時間序列的函數(shù)值，方框表示算法檢測到的異常區(qū)間。圖2a和圖2c說明HOT SAX算法和UKLD算法的精度相差不大，UKLD算法略優(yōu)于HOT SAX算法。圖2b說明雖然RKDE算法找到的第1個和第3個異常區(qū)間與HOT SAX算法和UKLD算法的無明顯差異，但第2個異常區(qū)間被RKDE算法分為2個連續(xù)區(qū)間，這會導(dǎo)致以后進(jìn)行異常區(qū)間分析時結(jié)果不準(zhǔn)確，所以RKDE算法的精度不如HOT SAX算法和UKLD算法的。

Figure 2 Anomalous intervals when length of the sequence is fixed圖2 序列長度為定值時的異常區(qū)間

因此，當(dāng)檢測序列長度為125時，UKLD算法的精度最優(yōu)，HOT SAX算法次之，RKDE算法較差。

從圖3中可以看出，當(dāng)檢測序列長度設(shè)為25,75,125,150,200,250(歸一化為0.1,0.3,0.5,0.6,0.8,1)時，隨著檢測序列長度的逐漸變大，3種算法的AP值逐漸升高，在歸一化值為0.5(序列長度為125)附近到達(dá)最高,隨后逐漸降低，算法的精度先升高后降低。對于檢測序列的固定長度的不同取值(25～250)，UKLD算法的AP值都要高于HOT SAX算法和RKDE算法的。在AP值最高點處，UKLD算法的精度在0.7附近，高于HOT SAX算法和RKDE算法的0.6附近的水平。即使在AP值最低(歸一化值為1，序列長度為250)時，UKLD算法的精度仍在0.5附近，高于HOT SAX算法和RKDE算法的0.4附近的水平。

Figure 3 Accuracy comparison when length of the sequence is fixed圖3 序列長度為定值時的量化精度

可見，對于檢測序列的固定長度的不同取值，UKLD算法的精度高于HOT SAX算法和RKDE算法的。

4.3 區(qū)間序列長度的精度分析

從圖4中可以看出，當(dāng)HOT SAX算法設(shè)定檢測序列長度為150，RKDE和UKLD算法設(shè)定序列長度取值為[25,175)時，3種算法都能找到該段時間序列中的3個異常區(qū)間。圖4a和圖4c說明UKLD算法較HOT SAX算法找到的異常區(qū)間的無用邊界部分更少，算法精度更高。圖4b說明雖然RKDE算法找到的第1個和第3個異常區(qū)間的精度介于UKLD算法與HOT SAX算法之間，但第2個異常區(qū)間被RKDE算法分為2個連續(xù)區(qū)間，這會使以后進(jìn)行異常區(qū)間分析時結(jié)果不準(zhǔn)確。

Figure 4 Anomalous intervals when length of the sequence is variable圖4 序列長度為取值區(qū)間時的異常區(qū)間

可見，當(dāng)檢測序列長度為取值區(qū)間時，UKLD算法的精度較HOT SAX算法和RKDE算法有明顯提升。

從圖5中可以看出，當(dāng)HOT SAX算法設(shè)定檢測序列長度為50,100,125,150,200，250，RKDE算法和UKLD算法設(shè)定序列長度取值為[25,75),[25,125),[25,150),[25,175),[25,225)和[25,275)時(以上都?xì)w一化為0.2,0.4,0.5,0.6,0.8，1)，隨著檢測序列長度取值區(qū)間逐漸變寬，UKLD算法和RKDE算法的AP值逐漸升高并趨于平穩(wěn)，算法的精度逐漸變高并趨于穩(wěn)定。由于HOT SAX算法只能把檢測序列的長度設(shè)為定值，故算法的AP值折線圖同圖3無異。因此,對于檢測序列長度為取值區(qū)間的情況，隨著取值區(qū)間逐漸變寬，序列長度固定的HOT SAX算法的精度明顯不如RKDE算法和UKLD算法的。

Figure 5 Accuracy comparison when length of the sequence is variable圖5 序列長度為取值區(qū)間時的量化精度

對于檢測序列長度的任意取值區(qū)間，UKLD算法的AP值都明顯高于RKDE算法的。隨著AP值趨于平穩(wěn)后，UKLD算法的精度穩(wěn)定于0.85附近，仍明顯高于RKDE算法的0.7附近的水平。

可見，對于檢測序列長度為取值區(qū)間的情況，UKLD算法和RKDE算法的精度指標(biāo)較序列長度固定的HOT SAX算法有明顯提升。對于檢測序列長度的任意取值區(qū)間，UKLD算法的精度明顯高于RKDE算法的。

5 結(jié)束語

為了從大量的多變量時空數(shù)據(jù)中發(fā)現(xiàn)數(shù)據(jù)異常的部分，本文提出了一種無偏KL散度算法(UKLD)。通過定義時空時間序列中的差異區(qū)間，在嵌入時間延遲后估計區(qū)間的分布，用累計和來加快參數(shù)估計過程，最后通過推導(dǎo)定義的無偏KL散度計算區(qū)間之間的差異水平來降低KL散度的計算偏差，將該差異水平作為檢測區(qū)間的異常得分，從而得到時空數(shù)據(jù)中的異常區(qū)間。仿真結(jié)果表明，UKLD算法在精度指標(biāo)上有明顯優(yōu)化。未來工作的方向主要是對時間間隔本身進(jìn)行深入分析來降低算法的計算成本，同時深入探索對于高維數(shù)據(jù)的有效概率密度估計來進(jìn)一步提升算法的準(zhǔn)確性。