田 騰,石茂林,宋學官,馬 躍,馮翔宇
(大連理工大學機械工程學院,遼寧大連 116024)
隨著物聯(lián)網(wǎng)和檢測技術(shù)的發(fā)展,工程系統(tǒng)中的大量時間序列數(shù)據(jù)被記錄下來,挖掘這類數(shù)據(jù)對于提升工程系統(tǒng)的設計、分析、運行和維護水平具有重要的意義[1-2]。然而,受制于檢測系統(tǒng)的精度和工程系統(tǒng)嘈雜的作業(yè)環(huán)境,所獲得的工程數(shù)據(jù)往往存在大量噪聲或異常數(shù)據(jù)[3-5]。因此,在進行數(shù)據(jù)挖掘與分析之前,就有必要進行數(shù)據(jù)的異常識別與歸類。
目前,時間序列的異常值檢測方法主要包括基于統(tǒng)計的方法[6-7]、基于時序模型的方法[8]、基于滑動窗口[9]的方法。相比于其他方法,基于滑動窗口的方法能夠?qū)r間序列進行分割,實現(xiàn)時間序列數(shù)據(jù)的降維處理,降低了計算復雜度,同時該方法便于理解規(guī)則又計算簡單[10],因此得到了廣泛應用。基于滑動窗口的檢測方法常用的子序列特征指標有:數(shù)據(jù)的空間距離、數(shù)據(jù)的統(tǒng)計分布和數(shù)據(jù)的斜率。然而基于數(shù)據(jù)空間距離、數(shù)據(jù)的統(tǒng)計分布的方式存在數(shù)據(jù)中心難以有效提取[11-12]、難以描述子序列數(shù)據(jù)結(jié)構(gòu)特征[13-14]的問題,基于數(shù)據(jù)斜率的方式能夠反映數(shù)據(jù)的結(jié)構(gòu)變化,進而得到了研究人員更多的關注。文獻[15]計算子序列內(nèi)相鄰兩點的斜率與設定的斜率范圍進行比較,以此判斷異常值,實現(xiàn)了電流信號的異常檢測。文獻[16]分別計算目標點與前、后數(shù)據(jù)的斜率進行異常判斷,實現(xiàn)了電網(wǎng)異常數(shù)據(jù)的檢測。文獻[17]計算子序列內(nèi)的最后一點與第一點的斜率變化作為子序列的特征信息,借助K-Mean實現(xiàn)異常檢測。文獻[18]采用子序列的斜率、均值、極值差等信息對子序列進行分段線性表示,使用層次聚類算法實現(xiàn)了水文異常數(shù)據(jù)的檢測。文獻[19]采用最小二乘法擬合序列斜率,根據(jù)斜率的相對變化值判定異常數(shù)據(jù),實現(xiàn)了大型立式磨床的故障預警。
工程數(shù)據(jù)中異常數(shù)據(jù)分布隨機性大、波動性大,而基于斜率的傳統(tǒng)檢測方法存在數(shù)據(jù)形態(tài)信息利用率低、數(shù)據(jù)結(jié)構(gòu)特征描述不全面、特征提取不準確等弊端,導致其難以適用于工程數(shù)據(jù)的異常檢測。為此,本文采用基于子序列斜率的置信區(qū)間的方式進行特征提取,并提出基于滑動窗口的時間序列異常檢測方法。仿真數(shù)據(jù)檢測結(jié)果表明,本文的方法提高了檢測精度。工程數(shù)據(jù)的異常檢測實驗表明,提出的算法能夠準確檢測異常值信息,能夠勝任工程數(shù)據(jù)的檢測工作。
該方法首先將時序數(shù)據(jù)分割為若干個子序列,隨后提取子序列的數(shù)據(jù)特征并用于識別異常子序列,最后采用Gath-Geva聚類算法識別異常子序列中的異常值。具體步驟介紹如下。
一條時間長度為n的時間序列表示為
Y=(y(t1),y(t2),…,y(tn))
式中y(ti)為ti時刻記錄的數(shù)據(jù),采集時間ti是嚴格增加的。
在時間序列的子序列分割中,采用長度為l(l< 圖1 滑動窗口工作原理圖 置信區(qū)間是數(shù)據(jù)統(tǒng)計分析中常用的技術(shù)手段,能夠以一定的可靠程度估計總體數(shù)據(jù)的所在區(qū)間。斜率能夠表述子序列的結(jié)構(gòu)特征,但是斜率信息分布隨機,為了準確描述子序列斜率信息分布情況,采用斜率的置信區(qū)間距離半徑作為子序列的特征表示,特征提取步驟如下。 假設其中某一個子序列為Yj(1≤j≤(n-l)/r+1),根據(jù)式(1)計算該子序列中任意相鄰2個數(shù)據(jù)點之間的斜率: (1) 則在該子序列Yj中,共得到l-1個斜率值,可得斜率的置信區(qū)間距離半徑為 (2) 置信上限為 (3) 置信下限為 (4) 基于提取的子序列特征,給出疑似異常子序列的判斷模型:假設含有n個變量的時間序列Y=(y(t1),y(t2),…,y(tn)),存在第j個子序列yj(1≤j≤(n-l)/r+1)斜率的置信區(qū)間距離半徑dj>γ,其中γ為閾值,則認為時間序列Y中的第j個子序列yj為異常子序列,其中含有異常值,待進一步實現(xiàn)異常值與正常值的歸類。 針對每個識別出的異常子序列,采用GG聚類算法對數(shù)據(jù)進行劃分,將異常值與正常值進行歸類。 設聚類樣本集合X={x1,x2,…,xN},現(xiàn)將數(shù)據(jù)X劃分為C類(2≤C≤N),設其聚類中心為V={v1,v2,v3,…,vC};設隸屬度矩陣為U=[uik]C×N,其中元素uik∈[0,1]代表第k個樣本數(shù)據(jù)屬于第i個簇的概率(1≤k≤N,1≤i≤C)。通過迭代調(diào)整(U,V)使目標函數(shù)取得最小值: (5) 具體步驟如下。首先計算聚類中心V={v1,v2,v3,…,vC}: (6) 式中h為迭代次數(shù)。 然后更新隸屬度矩陣: (7) 直到‖U(h)-U(h-1)‖<ε為止,其中,ε為允許誤差;否則繼續(xù)令h=h+1,重復上述步驟,直至滿足條件。 輸入:長度為n的時間序列Y=(y(t1),y(t2),…,y(tn)),異常子序列判斷閾值γ,滑動窗口長度l。 輸出:時間序列Y的異常值信息。 步驟1:使用長度為l的滑動窗口對時間序列Y進行等長度的劃分,分為若干子序列yj(1≤j≤(n-l)/r+1); 步驟2:采用式(1)計算每個子序列中相鄰兩點的斜率變化ki,采用式(2)計算斜率的置信區(qū)間距離半徑dj作為子序列特征; 步驟3:將距離半徑dj與閾值γ進行比較,初步確定含有異常值的異常子序列; 步驟4:通過GG聚類算法對異常子序列進行正常值與異常值的歸類與劃分,輸出時間序列Y中的異常值信息。 本節(jié)通過仿真數(shù)據(jù)和真實數(shù)據(jù)驗證所提出的基于滑動窗口的時間序列異常檢測方法的有效性,并通過如下指標評價實驗結(jié)果[17]: (1)查全率: (8) (2)查準率: (9) 仿真數(shù)據(jù)集為包含4個變量的長度為800的時間序列,表示如下: Y={y1,y2,y3,y4} y1=-sin(1+0.5x)+e,x∈[0.05,10]; y2=cos(1+0.5x)+e,x∈[0.05,10]; y3=log2(1+0.5x)+e,x∈[0.05,10]; y4=(-0.2x+1)2+e,x∈[0.05,10]。 式中:e為服從均值為0、方差為1的正態(tài)分布的隨機數(shù)。 在數(shù)據(jù)集中隨機插入20個異常值點Z~N(0,25)。異常點位置依次為3、9、29、60、72、164、235、244、358、475、518、540、549、565、606、614、624、652、678、704。 由于數(shù)據(jù)集中存在噪聲e,插入的異常值波動較大,為防止“大數(shù)吃小數(shù)”的現(xiàn)象發(fā)生,在實驗前對數(shù)據(jù)集進行歸一化處理,如下所示: (10) 首先研究窗口長度l對提出算法的影響。在實驗中,將步長r設定為1,即每完成一個子序列特征提取,就向后移動一個數(shù)據(jù)點。圖2為不同窗口長度條件下的實驗結(jié)果,從圖2中可以看出,所得結(jié)果的查全率和查準率均隨著窗口長度的增加而呈現(xiàn)先增長后下降的變化趨勢;在窗口長度l∈[7,11]時,實驗結(jié)果趨于平穩(wěn)。解釋如下:當窗口長度較小時,子序列中的相鄰數(shù)據(jù)點數(shù)量較少,導致難以通過數(shù)據(jù)斜率變化信息獲取子序列精確的內(nèi)部特征,從而導致子序列特征提取不準確,進而造成檢測精度較低;當窗口長度l∈[9,11]時,子序列的數(shù)據(jù)量充足,能夠以斜率信息準確表示子序列的數(shù)據(jù)形態(tài),通過GG聚類能夠精準識別異常值,因此查準率最高;當窗口長度l>11時,查全率R以及查準率P呈現(xiàn)下降趨勢,主要是因為檢測出的異常子序列中包含的數(shù)據(jù)較多,在異常子序列中不僅包含異常值,而且還包含了若干相鄰的正常值,當通過GG聚類算法進行數(shù)據(jù)劃分時,容易將部分正常值劃分為異常值,也容易將少量異常值劃分為正常值,因此查全率、查準率降低。綜合考慮以上實驗結(jié)果,最終滑動窗口長度選取為7,異常子序列判斷結(jié)果如圖3所示。圖4為刪除異常值以后的時間序列,可以看出,刪除異常值后的時間序列數(shù)據(jù)變得更加平穩(wěn),數(shù)據(jù)流中并未出現(xiàn)較大的峰值。 圖2 滑動窗口長度l對異常檢測的影響 圖3 異常子序列判斷結(jié)果圖 圖4 刪除異常值后的時間序列 表1中顯示的不同檢測算法的檢測結(jié)果,對比方法分別為文獻[13]提出的基于方差的時間序列異常檢測方法(AD-Variance)、文獻[17]提出的基于斜率的數(shù)據(jù)異常檢測方法(AD-Slope),滑動窗口長度設定為7。由于子序列特征提取方式不一致,造成子序列特征數(shù)據(jù)分布存在一定的差異,為準確識別異常子序列,分別設定最優(yōu)閾值:本方法閾值為γ=0.1,AD-Variance方法閾值為γ=0.05,AD-Slope方法閾值為γ=0.03,10次重復實驗結(jié)果如表1所示。 表1 不同檢測算法結(jié)果 從表1中的R、P值的方差信息來看,本文方法以及AD-Slope方法都出現(xiàn)了小幅波動,而AD-Variance方法方差為0。AD-Variance方法僅從閾值判斷異常值,所以當閾值設置完畢以后,就不會出現(xiàn)任何的波動情況,所以其方差為0。而本文方法和AD-Slope方法均采用聚類算法進行異常值的劃分,而聚類算法結(jié)果受隨機生成的初始矩陣影響,導致檢測結(jié)果存在一定的波動。本文方法的方差數(shù)值小于AD-Slope方法結(jié)果,表明GG聚類算法與K-Mean算法相比運行結(jié)果更穩(wěn)定。AD-Variance方法采用了置信區(qū)間的閾值判斷,方法較為簡單,所以其運行成本優(yōu)于本文方法。同時AD-Slope方法采用子序列的斜率信息進行閾值判斷,并未計算置信區(qū)間進行特征提取,所以其運行成本低于本文的運行成本。 從表1中R、P結(jié)果來看,本文與AD-Variance方法相比,查全率R以及查準率P分別提高了6.9%和4.3%,由于本文以子序列斜率的置信區(qū)間距離半徑作為子序列的結(jié)構(gòu)特征,更好地表示了子序列的特征變化,而AD-Variance方法采取子序列均值與方差作為子序列的特征,不能很好反映子序列的結(jié)構(gòu)變化情況,因此,在R、P評價指標數(shù)值上低于本文的方法。與AD-Slope方法相比,本文方法的查全率R以及查準率P分別提高了46.3%和67.4%,本文采用子序列中的全部斜率信息進行特征提取,而AD-Slope方法僅采用子序列的第一個以及最后一個數(shù)據(jù)求解斜率,因而本文方法能夠更好地保留了子序列中全部數(shù)據(jù)的形態(tài)信息,提高了子序列特征提取的準確性。從以上實驗結(jié)果可以看出,本文方法的查準率R以及查全率P優(yōu)于其余2種方法,但計算成本有所提升,考慮到本文提出方法具有較好的異常識別率,所提出的方法依舊具有較好的實用性。 在本小節(jié)中,將所提出的檢測方法應用于隧道掘進機(tunnel boring machine,TBM)實測數(shù)據(jù)的異常檢測,以驗證提出方法的工程可用性。該數(shù)據(jù)來源于國內(nèi)某城市地鐵隧道標段,隧道長度約2 km,直徑6.4 m,采用推進速度作為檢測對象,數(shù)據(jù)長度為1 800。 不同時間序列的子序列的特征數(shù)據(jù)分布具有差異性,因此本小節(jié)的最優(yōu)閾值設為γ=0.04;滑動窗口長度設為l=7,所得實驗結(jié)果如圖5所示??梢钥闯?,子序列特征曲線具有6個明顯的異常峰值,即異常子序列。針對檢測出來的子序列標號,形成如圖6所示的異常子序列。從圖6中可以看出,每一個異常子序列都包含了數(shù)據(jù)的異常值。 圖5 子序列閾值判斷 圖6 檢測出的異常子序列 針對這些異常子序列通過GG聚類算法實現(xiàn)異常值與正常值的歸類與劃分,表2中R、P為進行10次實驗的平均值;D1,D2分別為查全率R以及查準率P的10次實驗的方差,數(shù)值表示算法的波動情況。從表2中可以看出R、P都處于84%以上,且波動較小。針對隧道掘進機實測推進速度的異常檢測結(jié)果,本文提出的算法具有較高的檢測精度以及穩(wěn)定性,表明該方法具有良好的工程可用性。 表2 推進速度異常值檢測結(jié)果 本文提出了一種基于滑動窗口的時間序列異常檢測方法,通過斜率置信區(qū)間的方式解決了子序列特征難以提取的問題,采用滑動窗口方法和Gath-Geva聚類算法實現(xiàn)了異常值檢測。仿真數(shù)據(jù)的實驗結(jié)果表明,與以子序列方差信息和傳統(tǒng)的斜率信息提取子序列特征的檢測方法相比,本文提出的方法提高了檢測精度。隧道掘進機實測數(shù)據(jù)的異常檢測實驗表明,提出的算法能夠準確檢測異常值信息,能夠勝任工程數(shù)據(jù)的檢測工作。1.2 子序列特征提取與識別
1.3 Gath-Geva聚類算法識別異常值
1.4 異常檢測算法步驟
2 實驗仿真
2.1 仿真數(shù)據(jù)驗證
2.2 工程數(shù)據(jù)驗證
3 結(jié)論