基于異方差高斯過程的時間序列數據離群點檢測

嚴 宏，楊 波，楊紅雨

(1．中國民用航空飛行學院計算機學院，四川廣漢618307; 2．四川大學國家空管自動化系統技術重點實驗室，成都610064)(*通信作者電子郵箱yanhonggh@163．com)

0 引言

在數據挖掘和機器學習領域，異常檢測是指監測并判定數據集中與預期行為或模式不相符的單點數據、數據集合或數據序列［1］。異常數據通常與數據集中其他大部分數據在某種相似性度量上，比如歐氏距離(Euclidean distance)、馬氏距離(Mahalanobis distance)、皮爾遜相關系數(Pearson correlation)等，存在預定義或顯著的差異，往往預示著風險、安全事件或事故的發生，有助于采取預防或修正措施避免或減少由異常導致的損失，因此被廣泛應用于醫療、交通、金融、電力和氣象等行業。

時序數據是按時間先后順序測量或記錄的序列數據，蘊含著事物發展和變化的運行模式和內在規律。時間序列數據的異常檢測具有廣泛的應用場景，如網絡流量監測［2］、醫療數據分析［3］、水文時序數據分析［4］以及瓦斯濃度監測［5］等，這些場景中根據應用需求的差異所需要檢測的異常種類有所不同。根據異常特征以及表現形式的不同，時間序列數據異常可以分為數據點異常、前后關聯異常和子序列異常［1］。本文關注的是數據點異常，對其的檢測通常稱為離群點檢測。

目前，針對時間序列數據的離群點檢測方法多種多樣，其中主要的檢測方法有3種:1)基于分類的方法［6－8］，這類方法使用已標注是否異常的數據訓練模型，能夠有效利用訓練樣本的特性進行離群點判定，但是受限于訓練樣本的數據分布，當某些時間點上的樣本缺失或不足時，將影響模型的異常檢測性能。另外，該類方法需要對訓練數據進行異常標注，大量數據下的標注工作將會增加該類方法的投入成本，降低可行性。2)基于聚類的方法［9－11］自動將數據集分為若干簇類，不屬于任何簇類或簇類中數據數目較少的情況判別為離群點，該類方法屬于無監督學習，避免了模型訓練時對標注數據的依賴，但是依賴于對聚類模型參數的合理設置，不適當的參數設置將會嚴重地影響異常檢測效果。3)基于預估模型的方法［2，12－13］通過構建模型計算正常模式下檢測時間點的數據預測值，然后與實測值比較完成離群點判定。該類方法具有較好的直觀性和解釋性，但是難點在于需要使模型的預估數據與正常情況下的實際數據相符，或者是偏差小于合理的閾值。求取檢測時間點的數據預測值屬于回歸問題，采用線性回歸、多項式回歸以及支持向量回歸等模型只能求取檢測時間點的單一預測值。在考慮容許一定范圍偏差的實際應用情境中，如采用固定大小的閾值判定偏差有異常，那么又將面臨閾值設置合理性以及閾值固定不變帶來的問題。針對上述問題，本文提出一種基于高斯過程對時間序列數據正常模式進行建模的方法，利用高斯過程模型的特性求取正常情況下標準值和偏差數據的概率分布，構建的預估模型能夠輸出具有容差區間的預測值，并使用公開的真實網絡流量數據進行實際應用場景下的有效性和性能驗證。

1 問題及方法描述

1．1 問題描述

為了便于描述，采用如下符號對研究問題中涉及的數據進行統一的形式化描述。使用符號t通用地表示時間序列中任一時間點，本文將t由協調世界時(Coordinated Universal Time，UTC)之類的標準時間轉換為以某個時間點為基準的相對時間，所測量的數據值由y表示。在構建預估模型之前，需要獲取用于訓練模型的數據集 D={(t1，y1)，(t2，y2)，…，(tN，yN)}，該訓練集共有N組數據，每組數據中兩個數據項y和t一一對應。如果進一步對數據進行人工標注，將每組數據(ti，yi)(i=1，2，…，N)，使用對應符號δ表明是否異常，δi為1表示數據異常，δi為0則表示數據正常，那么將會得到標注訓練集 D={(t1，y1，δ1)，(t2，y2，δ2)，…，(tN，yN，δN)}。使用訓練集對模型進行訓練后，對于需要檢測的時間點t*，模型會輸出表示時間序列數據正常行為或模式下t*對應的y*;但是不同于常見的回歸問題求解，由于通常情況下存在各種不可避免的干擾因素，測量的數據容許出現一定范圍的偏差，因此，模型最終需要輸出的并非是單一數值，而是一個容差區間，表示正常情況下y*的取值范圍。

1．2 方法概述

本文借助于高斯過程對時間序列數據的正常行為或模式進行建模，高斯過程近些年取得了豐富的研究成果，不僅在機器學習領域近些年的重要著作［14－15］中得以著重論述，Williams等［16］還為機器學習中的高斯過程撰寫了專著。高斯過程作為一類隨機過程，常用于處理非線性等復雜回歸問題，通常按如下形式定義:

其中m(t)和k(t，t')分別為高斯過程的均值函數和協方差函數，定義如下:

通常情況下由于信息有限，無法預先獲知相應的均值函數，常見方法是使用常值函數0作為均值函數，這樣整個高斯過程交由協方差函數確定，這種做法同樣能夠較好地描述整個過程中變量分布的變化情況，并且還在一定程度上降低了計算的復雜度［16］。由于高斯過程使用貝葉斯推斷理論不僅可以計算得出預測值的最大后驗估計，還能夠得到一個估計的預測值概率分布。本文正是利用高斯過程該項特點計算得出正常情況下需要檢測的時間點上數據值的預計分布情況，然后確定容差區間進行離群點檢測。

2 離群點檢測模型

2．1 時間序列數據分解

由于傳感器誤差、事物內在可變性、外界環境影響等不可避免的干擾因素，時間序列數據中測量的數據值往往存在合理范圍內的偏差，為此建模過程中將實際測量的數據y作如下分解:

其中:si表示的是時間點ti上測量值在無偏差時應當得到的理想標準值，ri即為yi中的偏差項。通過此方法完成分解后，時間序列中測量數據的形成過程更加明確，具有可解釋性，然后分別對兩個構成部分建模求取正常情況下的數據分布情況。

2．2 標準值建模

模型構建過程中使用高斯過程從部分已知標準值的基礎上求解需要檢測的時間點t*上標準值s*的后驗估計，具體方法如下:

假定時間點集t={ti|i=1，2，…，N}對應的數據標準值s={si|i=1，2，…，N}已知并構建對應向量，如上所述使用均值函數為常值0的高斯過程對標準值建模，協方差函數選用常見的squared exponential核函數［16］，那么由高斯過程的概率分布定義可知，s和s*滿足如下的多元高斯分布:

其中KN+1為此多元高斯分布的協方差矩陣，具有如下的形式:

其中

由此可以應用高斯過程的推導結果［14－16］計算出需要檢測的時間點t*上標準值s*的后驗分布如下:

2．3 偏差項建模

從建模過程中時間序列數據的分解可知，式(1)中的偏差項r類似于大多數高斯過程模型中引入的噪聲項［16］，但是通常情況下都假設各個時間點的噪聲項符合高斯分布，均值為0并且其方差恒定不變。這種假設主要是為了簡化后續推導，使得后驗估計能夠直接求得解析解，然而在現實中存在不同時間點上受干擾程度不一樣導致合理的數據偏差程度發生變化的場景，為此本文假設偏差項r滿足均值為0的高斯分布，但是其方差不再為固定數值。為了求得偏差項r與時間點t的函數關系，本文使用另外一個高斯過程對偏差項建模，更確切地講是對偏差項方差的對數建模，主要是為了保證偏差項方差的非負特性，如下所示:

其中kr(t，t')是偏差項方差對應高斯過程定義所使用的協方差函數，同樣選用squared exponential核函數，而常值μr用于表示其偏差項方差的平均水平。

2．4 后驗分布求解

通過上述對標準值和偏差項分別建模后，根據線性高斯模型相關理論［15］可知，式(1) 表明向量 y=(y1，y2，…，yN)T和向量s=(s1，s2，…，sN)T之間滿足線性關系，已知的測量數據y1，y2，…，yN和檢測時間點t*對應的標準值s*仍然滿足多元高斯分布，其協方差矩陣為:

其中 KD是關于時間點 t1，t2，…，tN上偏差項的對角矩陣diag(r)，對角元素r=(r1，r2，…，rN)T。至此，可以再次利用高斯過程的推導結論［14－16］，可以計算出關于 s*的后驗分布p(s*|t*，y，t，r，r*)，附加上偏差項 r*就可以得到時間點 t*上正常情況下預估數據y*的后驗分布:

其中

由于式(4)和(5)中KD涉及的r以及式(5)中的r*未完成求解，因此為了完成對y*的后驗分布求解，需要進一步對其中的r和r*進行積分，然后得到:

僅從目前的模型假設以及已知時間序列數據中無法求解p(r，r*|t*，y，t)，盡管可以采用文獻［17］中提到的采樣方法求取近似解，但是存在計算量大、耗時長的缺點。本文采用文獻［18］中的思路利用變分推斷求取近似解，使用變分推斷中常用的平均場(mean field)方法［15］，對邊緣概率p(y)的對數作如下分解:

其中KL(·‖·)表示KL散度(Kullback-Leibler divergence)，而L(q(s)，q(r))是ln(p(y))的下界，當尋求分布q(s)和q(r)最大化該下界時，也就會使式(6)中的KL散度最小化，從而優化p(s，r|y)在分解形式q(s)q(r)下的近似。采用文獻［18］中下界L(q(s)，q(r))關于q(s)最大化的推導結論，可以進一步得到一個僅依賴于q(r)的近似下界:

其中

對于q(r)的分布采用變分推斷中常用的多元高斯分布，其均值向量和協方差矩陣分別用μq和Σq表示，結合之前對s多元高斯分布的設定以及式(2)和(3)，可以進一步推導式(7)中的下界為如下形式:

其中:tr(·)表示矩陣的跡，R為一對角矩陣，其對角元素為［R］ii=exp(［μq］i－ ［Σq］ii)，而 Kr是使用協方差函數 kr(t，t')計算的協方差矩陣。根據式(8)的極值點與偏導數的關系，可以求得:

其中Λ表示半正定對角矩陣。至此，模型的訓練可以歸結為對最大化式(8)表示的下界，其中需要優化的參數包括標準值s和偏差項r對應的兩個高斯過程協方差函數中的參數、矩陣Λ中的對角線元素以及用于控制偏差項方差平均水平的μr，本文選用模型訓練中常見的共軛梯度法［19］進行參數優化。使用訓練集完成參數優化后，利用變分推斷結果可以得到s*的后驗分布:

其中

根據式(2)和(3)中對偏差項方差的假設，代入其在時間點t*上的期望值也是最大概率值作為近似求解y*的偏差項估計，得到如下結果:

其中y*的均值與式(9)等同，而方差為式(10)中的標準值方差加上偏差項引入的方差:

其中 kr*

2．5 離群點判定

通過上述推導過程得到檢測時間點t*上y*的后驗分布后，可以選取其中特定范圍的數據分布區間作為正常情況時間點t*上測量數據的容差區間，如果需要檢測的實際數據y不位于該區間，那么就將其判定為離群點數據，相應的判別函數可用如下形式描述:

其中關于a*和d*來自于式(9)和式(11)的推導結果，而zα由模型中設置容差區間占比的超參數α確定，超參數α并不從訓練數據中學習，而是根據實際應用需求進行調整設置。例如，使用高斯分布中常用的95%置信區間作為容差區間時，即將超參數α設置為95%，那么zα可近似取值為1．96。

2．6 異常訓練數據過濾

在訓練集D中數據沒有異常標注項δ的情況下，包含在其中的異常數據會引起模型訓練中常見的噪聲數據干擾問題。由于高斯過程本身對噪聲數據具有一定的抗干擾能力［16］，本文利用該項特性設計了一種迭代訓練模型的方法，基本思路是通過當前模型輸出反復過濾當前訓練集然后重新訓練直至過濾的異常數據小于既定的容差比例或者達到最大訓練次數為止，具體操作如下描述:

輸入 訓練集 D={(ti，yi)|i=1，2，…，N}，容差區間占比α，最大迭代訓練次數m。

輸出 離群點檢測模型參數θ。

過程 離群點檢測模型訓練。

1) 設置集合Ds=D，Da=

2) 隨機初始化模型參數θ

3) t=0

4) repeat

5) 設置當前訓練集Ds=DsDa

6) 使用Ds訓練模型，更新參數θ

7) 當前異常數據集Da=

8) for(tj，yj) ∈ Dsdo

9) 計算時間 tj的容差區間［aj－ zαdj，aj+zαdj］

10) if yj［aj－ zαdj，aj+zαdj］then

11) 添加異常數據Da=Da∪(tj，yj)

12) end if

13) end for

14) t=t+1

15)until(|Da|/|Ds|)≤(1－α)∨(t＞m)

上述訓練方法能夠使模型輸出的容差區間反映出無標注數據集中大部分數據的分布情況，這將在下一章中得以驗證。

3 實驗結果及分析

3．1 效果驗證

為了驗證對時間序列數據中離群點數據的檢測效果，使用雅虎公司Webscope項目［20］提供的公開時序數據進行測試，該數據來源于雅虎公司真實場景中網絡節點的流量統計，用于研究網絡流量時序數據中的異常檢測，內容包括全天24小時以整點開始每個小時內的流量統計，并且包含每個流量數據是否異常的人工標注。實驗中使用Matlab編寫模型代碼，版本為 Matlab R2013a，操作系統使用 Windows 7 Professional，硬件環境為 Intel Core i7-6700K 4．0 GHz 以及16 GB DDR4-2133。由于篇幅有限，選取Webscope項目中三個場景作為案例演示，首先基于人工標注的正常流量數據，分別采用高斯過程建模中常見的偏差項方差恒定不變和異方差偏差項兩種方法對正常流量數據的容差區間進行建模，然后為了驗證針對無人工標注是否異常的模型訓練方法效果，將人工標注的異常數據和正常數據全部用于容差區間的模型訓練，分析對比了三種方法在各個場景下容差區間建模和離群點檢測效果。

1)場景一。

圖1分別顯示場景一中的網絡流量數據分布和使用上述三種方法完成的容差區間建模結果，由于Webscope項目原始數據中沒有流量單位說明故而縱坐標流量單位未標注，數據對應橫坐標為流量統計的整點時刻。

圖1 場景一的實驗效果Fig．1 Experimental results of scenario 1

從圖中數據分布可以看出該場景下正常流量數據在各個時間點上邊界值以及分布區間長度方面較為相似，使用三種方法計算得出的容差區間在直觀效果上幾乎相同，對于遠離正常流量數據分布區間的異常數據都能實現成功檢測，對于孤立的邊界正常數據也都出現誤報的現象，但這也從一定程度上體現了這些孤立邊界數據和集中分布正常數據之間的差異性。實驗結果表明在各個時間點上正常數據分布相似的情況下，采用異方差高斯過程建模能夠取得與常見的偏差項方差恒定假設類似的建模效果，計算得出的容差區間長度在各個時間點上近似相同，與正常數據在各個時間點的分布情況相符。

2)場景二。

圖2將分別顯示了場景二中的網絡流量數據分布和使用上述三種方法完成的容差區間建模結果。

不同于場景一，該場景中正常流量數據在各個時間點上差異較大，在使用標注正常數據訓練的情況下，若采用高斯過程建模中常見的偏差項方差恒定的假設，必然會使計算得出的容差區間長度在各個時間點上保持不變，與實際正常數據分布情況不符，這在圖2(a)中得到較為直觀的驗證，導致時間點5、6、7和8時，容差區間過小而出現較多正常數據誤報。使用異方差高斯過程建模后，盡管還是出現了部分孤立邊界數據誤報的情況，但是在時間點5、6、7和8時，圖2(a)中部分誤報的數據在圖2(b)中處于容差區間內，誤報率得以下降，并且從圖2(b)中可以看出，計算得到的容差區間在各個時間點上的變化與該場景下正常流量數據在各個時間點上的合理偏差趨勢更為相符。此外，對于全部數據無標注是否異常的情況，如圖2(c)所示，本文訓練方法計算的容差區間相對于圖2(b)在時間點7、8、16、17和18上的容差區間更小，出現了更多的誤報數據，但是計算得出的容差區間其長度在各個時間點上的變化還是與正常數據的合理偏差趨勢基本相符，并且完成了所有標注異常數據的檢測。

圖2 場景二的實驗效果Fig．2 Experimental results of scenario 2

3)場景三。

圖3將分別顯示了場景三中的網絡流量數據分布和使用上述三種方法完成的容差區間建模結果。

與場景二類似，該場景中實驗結果再次表明:采用異方差高斯過程構建預估模型的離群點檢測方法，可以如實地描述正常流量數據在各個時間點上的分布變化情況;而對于無標注的流量數據，本文模型訓練方法存在對訓練數據的異常判定和迭代過濾，相對于直接使用標注正常數據訓練模型，會使得較為稀疏的正常邊界數據中更多數據點被判定為異常，這也正是圖3(c)所示在時間點4～17上的容差區間相對于圖3(b)更小，出現較多誤報數據的原因，但還是與正常數據偏差范圍的變化趨勢基本相符。

圖3 場景三的實驗效果Fig．3 Experimental results of scenario 3

與場景二不同的是，該場景下高斯過程建模中常見的偏差項方差恒定假設不僅導致時間點10～20上誤報率過高，而且在時間點2和3上的標注異常數據還被判定為正常，出現了異常漏報的情況，進一步體現了假設偏差項方差恒定的方法無法描述正常數據偏差范圍變化趨勢的缺陷。

3．2 性能對比

為了進一步與目前時序數據離群點檢測的常用方法進行性能比較，選取相關研究中基于一類支持向量機(One-class SVM)［10－11］、自 回 歸 積 分 滑 動 平 均 模 型 (Autoregressive Integrated Moving Average Model，ARIMA)［2，12］以及基于密度并伴隨噪聲的空間聚類算法(Density-based Spatial Clustering of Application with Noise，DBSCAN)［9］三種方法作為對比，選用異常檢測中常用的誤報率、召回率和F1-score作為性能對比的度量指標。實驗環境中使用的硬件和軟件與3．1節中的描述相同。實驗中本文模型中容差區間占比α設置為95%，鑒于此基于One-class SVM的模型其訓練樣本異常檢出比率上界設置為5%，核函數采用常用的徑向基函數(Radial Basis Function，RBF)，使用通常的網格搜索進行參數調優，經過預先的參數設置對誤報率、召回率和F1-score的影響測試后，將RBF核函數參數 γ 搜索區間設定為［0．000 1，0．01］，搜索步長為0．0001，最終選取三次參數設置作為示例說明，包括使得F1-score在其中取值最高的參數設置。基于ARIMA的模型根據ADF檢驗(Augmented Dickey-Fuller test)確定差分階數d，通過貝葉斯信息準則(Bayesian Information Criterion，BIC)確定自回歸和移動平均階數p和q，同樣鑒于本文模型的容差區間占比設置將預測值的95%置信區間作為正常數據區間，不位于該置信區間的數據判定為離群點。DBSCAN算法雖然不像k-means等聚類算法需要預先確定聚類簇數，但其效果還是與鄰域半徑ε和核心對象鄰域內最小對象個數MinPts兩個參數密切相關，為了尋求能使反映誤報率和召回率兩方面的綜合性能指標F1-score取值最大的參數設置，同樣借鑒網格搜索的思路進行多組參數的比較，鄰域半徑ε搜索區間設置［0．01，0．1］，搜索步長為 0．01，鄰域內最小對象個數MinPts搜索區間設置［5，15］，搜索步長為1，最后除了選取F1-score取值最大的參數設置進行模型對比，還選取了其他4組能使誤報率或召回率取得較好結果的參數設置作為示例說明。實驗中對雅虎公司Webscope項目中選用的50個時序數據集進行數據歸一化預處理，避免各個時序數據集數值范圍不同造成的影響，各個對比模型、參數說明和性能指標數據如表1所示。

表1 性能指標對比Tab．1 Comparison of performance indicators

基于One-class SVM的模型在網格搜索參數過程中，當RBF核函數參數γ設置為0．001時F1-score性能指標數值最高，在所有模型對比中異常檢測效果也相對較好，但該模型會受到核函數參數設置的影響。從表1中實驗數據可以看出當參數γ增加到0．01時會提高召回率，但也導致誤報率過高，另外當參數γ降低到0．0005時出現了參數設置不當造成召回率降低并且誤報率反而增加的情況。基于DBSCAN的模型的性能指標同樣受到參數設置的影響，使用(ε=0．05，MinPts=12)一組參數設置時才取得誤報率和召回率之間較好的權衡考量，F1-score性能指標在該模型所有參數設置中最好，使用(ε =0．06，MinPts=12)和(ε =0．05，MinPts=10)兩組參數設置時能進一步降低誤報率，但也導致召回率下降，使得 F1-score性能指標降低，使用(ε =0．04，MinPts=12)和(ε=0．05，MinPts=15)兩組參數設置時召回率取得較高數值，但卻出現了過高的誤報率，難以投入實際應用。基于異方差高斯過程模型的離群點檢測方法在使用標注正常數據訓練模型的情況下，與高斯過程建模中常見的偏差項方差恒定不變的方法相比，能夠計算得出更加合理的容差區間，取得顯著的性能提升。盡管該模型在誤報率和召回率兩個單項指標上沒有取得所有實驗結果中的最高數值，但在綜合指標F1-score上相對于其他模型都取得了一定程度的提升。在使用全部無標注數據訓練的情況下，該模型在召回率和F1-score性能指標上也取得較為滿意的結果。此外，該模型另一優勢在于高斯過程模型相關參數通過訓練集優化確定，避免了其他模型中出現的因參數設置不當造成誤報率過高或召回率過低的情況。

4 結語

本文基于預估模型檢測時序數據離群點檢測方法并沒有直接針對監測數據進行數學建模，而是首先將監測數據分解為標準值和偏差項兩個部分，這種做法與常見的高斯過程建模加入噪聲項的方法類似，但是區別在于并沒有假定偏差項獨立同分布以致于方差恒定不變，而是再次使用高斯過程對各個時間點的偏差項建模，從而能夠基于異方差高斯過程對不同時間點上正常數據合理偏差范圍變化實現有效的數學描述。通過實驗數據表明，本文的離群點檢測方法能夠取得誤報率和召回率兩個方面較好的權衡，并且無需考慮關鍵模型參數的人工設置，避免參數設置不當對性能指標的嚴重影響。在之后的研究工作中，還需要進一步考慮不同應用場景下高斯過程協方差函數的選取以及容差區間占比設置對于離群點檢測性能的影響以及改進。