魯法明,王 琳,包云霞,李 昂,曾慶田
(1.山東科技大學 計算機科學與工程學院,山東 青島 266590;2.山東科技大學 數(shù)學與系統(tǒng)科學學院,山東 青島 266590)
時間序列是對某個物理量進行等時間間隔觀測所得到的數(shù)值序列,時間序列分析廣泛應用于工業(yè)制造、金融、氣候監(jiān)測和醫(yī)療等領域。例如,在工業(yè)制造領域,采集并分析生產(chǎn)過程中的設備狀態(tài)時序數(shù)據(jù),可以對生產(chǎn)線異常做出預警,提高產(chǎn)品質量和生產(chǎn)效率[1-3];在經(jīng)濟金融領域,股票的價格走勢、利率等信息也構成時間序列,對其進行分析可以在一定程度上預測股票收益或檢測股票操縱異常[4-6];在氣候環(huán)境領域,對降雨量、河流水位、溫度等時間序列進行分析可以為氣候應對等問題提供依據(jù)[7-9];在醫(yī)療領域,持續(xù)監(jiān)測病人的心電(Electrocardiogram, ECG)和腦電(Electroencephalogram, EEG)活動情況便得到相應的時間序列,對其進行分析可以及時發(fā)現(xiàn)病情的異常情況[10-11];在公共管理領域,對公交系統(tǒng)客流量和車流量時序數(shù)據(jù)進行分析可以為公交管理和路線調整提供依據(jù)[12]。
時間序列數(shù)據(jù)通常具有維度高和數(shù)據(jù)量大等特點,直接分析和挖掘原始數(shù)據(jù)的運算量較大,因此時間序列的降維表示成為業(yè)界研究熱點之一。目前,比較成熟的時間序列降維表示方法有離散傅里葉變換(Discrete Fourier Transform, DFT)[13]、離散小波變換(Discrete Wavelet Transformation, DWT)[14]、分段線性近似(Piecewise Linear Approximation, PLA)[15]、分段聚合近似(Piecewise Aggregate Approximation, PAA)[16]和符號聚合近似(Symbolic Aggregate approximation, SAX)[17]等。其中,SAX廣泛應用于時間序列數(shù)據(jù)的挖掘與分析,其將時間序列轉換為字符串序列,在達到降維效果的同時較好地保持了時間序列的原始信息[16]。
SAX[16-17]將原始時間序列分為若干時間段,每個時間段內的觀測值用一個符號聚合表示,從而快速有效地對時間序列降維。降維后的符號序列仍然可以反映時間序列的總體變化態(tài)勢,且具有快速高效等優(yōu)點。然而,SAX在時間序列的符號化過程中丟失了原始時序數(shù)據(jù)的許多統(tǒng)計和形態(tài)信息。統(tǒng)計信息主要指極值、方差等信息,形態(tài)信息主要包括變化趨勢和波動信息等。
為減少統(tǒng)計信息的丟失量,LKHAGVA等[18]提出一種基于極大值、極小值和均值的擴展符號聚合近似(Extension of Symbolic Aggregate approximation, ESAX)方法,該方法同樣將每個子序列的最大值和最小值加入符號序列,使原來每個子序列對應的符號變成3個。ESAX彌補了部分丟失的統(tǒng)計信息,強調極值信息的重要性,適用于部分對數(shù)據(jù)極值敏感的場景。
同樣是針對統(tǒng)計信息丟失的問題,鐘清流等[19]在SAX的基礎上加入方差來描述時間序列的發(fā)散程度。該方法將時間序列的符號特征看做矢量,將各個子序列的方差和均值特征分別看做組成該矢量的兩個分量,在一定程度上彌補了時序數(shù)據(jù)在方差方面丟失統(tǒng)計信息的問題,然而當某些應用場景中存在形態(tài)特征不同、均值和方差相同的時序數(shù)據(jù)時,該方法需要進一步優(yōu)化。
針對SAX在符號化過程中丟失形態(tài)特征信息的問題,李海林等[20]提出SF_SAX。該方法在SAX的基礎上,采用最小二乘法對子序列進行直線擬合,用直線的斜率近似值表示該子序列的形態(tài)特征,然后按文獻[19]的方法將均值和斜率表示為符號矢量。SF_SAX能較好地識別時間序列在形態(tài)方面的差異,然而當壓縮比較大且序列波動較大時,擬合效果需要進一步改進。
為捕捉時間序列中數(shù)據(jù)值的變化趨勢信息,SUN等[21]提出一種基于趨勢距離的時間序列符號化(SAX_TD)方法,該方法將每個子序列的起點值和終點值與其均值的差稱為趨勢因子,通過計算趨勢因子之間的距離衡量時間序列趨勢的變化,能較好地識別時間序列的變化趨勢;與文獻[21]類似,同樣使用子序列的起始點來構建趨勢,季海娟等[22]提出基于始末距離的時間序列符號聚合近似表示方法(SAX_SM),該方法在計算子序列均值時去掉起始點,而將起始點加入子序列映射的字符序列中;李海林等[23]用數(shù)值導數(shù)描述子序列的趨勢特征,導數(shù)為正數(shù)表示子序列呈上升趨勢,為負數(shù)表示子序列呈下降趨勢,為0則表示無變化。這3種方法在一定程度上彌補了符號化過程中趨勢信息丟失的問題。
前述各種改進的SAX在對時序數(shù)據(jù)進行降維的同時,分別從不同角度捕捉傳統(tǒng)SAX丟失的統(tǒng)計信息和形態(tài)信息,各有其適用的場合。本文著眼于彌補時間序列符號化過程中丟失的數(shù)值波動信息。具體而言,本文提出一種融合波動信息的時間序列符號聚合近似方法,定義了一種新的波動率指標來同時刻畫時間序列的波動幅度和變化趨勢等形態(tài)特征,用融合了波動率的符號矢量近似刻畫子序列,在此基礎上給出一種新的時間序列距離度量方法。然后,以該度量方法為基礎,提出時間序列的相似性計算和分類方法,并在公開的加利福尼亞大學河濱分校(University of California, Riverside,UCR)時序數(shù)據(jù)集和凱斯西儲大學(Case Westem Reserve University, CWRU)的軸承故障診斷數(shù)據(jù)集上對該分類學習的效果進行實驗,對所提方法的適用場景進行分析。
SAX是LIN等[16-17]提出的基于PAA的符號化表示方法。其在時間序列數(shù)據(jù)近似服從正態(tài)分布的前提下,使用符號序列聚合表示原始時間序列,在降維的同時保留時間序列的總體變化態(tài)勢,并保證符號空間相似模式之間的距離滿足真實距離的下界要求。
簡而言之,SAX將時間序列轉換為字符序列主要經(jīng)過以下3步:
步驟1將原始時間序列規(guī)格化,轉換成均值為0、標準差為1的序列。
(1)
表1 字符集大小α=3,…,10時各分割點的取值情況
假設某時間序列數(shù)據(jù)經(jīng)過規(guī)格化處理后的結果如圖1所示,圖中在128個觀測點取值。按照SAX設置分段數(shù)量為8,字符集大小為3,首先采用PAA算法得到各段子序列的均值,然后通過符號聚合近似后得到符號序列abcccbaa,該符號序列對應的圖形表示如圖2所示。對比圖1和圖2可見,采用SAX進行降維后,數(shù)據(jù)維度從原始的128個觀測值降到只有8個觀測值,而時間序列的總體變化態(tài)勢在降維前后基本保持一致。
(2)
(3)
式(3)中β取值如表1所示。
由上述步驟可見,SAX采用PAA計算出的子序列均值代替這一段子序列,容易丟失原序列部分信息(如方差、變化趨勢和波動情況等),而且數(shù)據(jù)的壓縮比越大,信息丟失得越多。另外,該方法只能保留原時間序列的總體變化趨勢,無法描述各段的局部形態(tài)信息。針對SAX的不足,下面給出融合波動率的時間序列符號聚合近似方法。波動率從一定程度上同時刻畫時間序列的局部波動幅度、變化趨勢等特征,可以彌補SAX部分丟失的信息。
時間序列的波動情況同時隱含時間序列的數(shù)值和形態(tài)方面的特征,這些特征往往與時間序列的狀態(tài)或模式變化有關,捕捉時序數(shù)據(jù)的波動信息對時間序列分析和異常檢測等應用具有重要作用。觀察圖1所示的時間序列,可以直觀地發(fā)現(xiàn)各個子序列局部的狀態(tài)或模式變化,但SAX處理后僅保留各個子序列的均值,導致丟失了時間序列信息,從圖形上看,圖2相比圖1丟失了很多信息。為更好地刻畫時間序列的波動信息,本文在SAX的基礎上加入波動率來量化時間序列的波動特征,其定義和計算公式如下:
定義1給定時間間隔0,記xi和xi+1分別為第i時刻和第i+1時刻的觀測值,則稱式(4)的計算結果為本時間段內的波動率。
(4)
由式(4)可見,當時間序列各個觀測點的取值整體呈上升或下降趨勢時,波動率的取值較大;當沒有明顯的上升和下降趨勢,但局部波動幅度較大時,波動率的取值也較大。反之,當某個時間段內各個觀測點的取值均穩(wěn)定在某個恒定值附近時,波動率的取值接近零。顯然,該式客觀地刻畫了時間序列的波動情況,后文實驗結果也表明這一方法在捕捉時間序列形態(tài)信息方面的有效性。
結合上述波動率計算公式,融合波動率之后的時間序列符號化方法稱為融合波動率的時間序列符號聚合近似方法(time series Symbol Aggregation approximation method for fusion VOLAtility,SAX_VOLA),其具體過程如下:
給定長度為n的時間序列C={c1,c2,c3,…,cn},SAX_VOLA對時間序列進行符號表示時主要有以下4步:
步驟1與SAX相同,首先對原時間序列進行規(guī)格化,將其轉換成均值為0、標準差為1的序列。
例如,對于圖1中的時間序列,仍然設置分段數(shù)量為8,字符集大小為3,先采用PAA算法得到各段子序列的均值,再進行符號聚合近似,得到的符號序列仍然為abcccbaa。SAX_VOLA針對每個子序列計算其波動率,并將波動率取值與對應子序列的符號組合為一個二元矢量,得到的二元符號矢量序列為
c·i+0.334·j,c·i+0.375·j,
c·i+0.365·j,b·i+0.426·j,
a·i+0.279·j,a·i+0.246·j}。
(5)
(6)
上述時間序列的距離度量方法在傳統(tǒng)符號距離的基礎上融入波動率距離,而波動率在一定程度上同時刻畫了時間序列的波動幅度和變化趨勢,因此相比傳統(tǒng)的符號化方法,該方法通常能更加準確地度量時間序列的距離。為驗證這種距離度量的準確性,后文將從某公開的數(shù)據(jù)源中尋找含有分類標簽的多個數(shù)據(jù)集,基于該相似性度量指標對時間序列進行分類學習,通過分類的準確性評估本文所提符號聚合方法和時間序列距離度量方法的有效性。
UCR時間序列檔案庫[24]是使用最為廣泛的時間序列數(shù)據(jù)源之一,它含有128個時間序列數(shù)據(jù)集,具體包括ECG數(shù)據(jù)集、Trace故障數(shù)據(jù)集、Wafer傳感器數(shù)據(jù)集等。其中的Trace數(shù)據(jù)集[25]是流程工業(yè)領域核電站儀器故障的綜合數(shù)據(jù)集(完整的數(shù)據(jù)集一共有16類,此處UCR時間序列檔案庫收錄了Trace數(shù)據(jù)集中的4類[26]);Wafer數(shù)據(jù)集[27]是在半導體制造領域的硅晶片加工期間,通過各種傳感器收集的線上測量值所構成的數(shù)據(jù)集,分為正常和異常兩類。文獻[21-22]以UCR為數(shù)據(jù)源,對多種時間序列符號聚合近似表示方法在時序數(shù)據(jù)分類方面的準確性進行對比。為驗證本文所提SAX_VOLA的有效性,同時便于不同方法間進行比較,采用文獻[21-22]采納的20個時序數(shù)據(jù)集進行實驗評估。具體選擇的UCR中的數(shù)據(jù)集如表2所示,每個數(shù)據(jù)集包括訓練集和測試集兩部分,不同數(shù)據(jù)集類別標簽的數(shù)目從2~50不等,時間序列長度從60~637不等,這些數(shù)據(jù)集在類別數(shù)量和序列長度上的多樣性可以更客觀地對各種符號化方法的適用性做出評估。
表2 選擇UCR的20組數(shù)據(jù)集信息表
續(xù)表2
本文借助時間序列數(shù)據(jù)分類的準確性對比不同符號聚合近似表示方法。在分類器的選擇方面,鑒于K最鄰近(K-Nearest Neighbor, KNN)分類器[28]原理簡單,差錯率較低,本文用KNN分類器作時間序列數(shù)據(jù)的分類算法??紤]到1-NN分類器直接按照訓練集中與其距離最近的時間序列類別進行分類,采用KNN分類器中的1-NN分類算法對時間序列進行分類。
就時間序列距離度量而言,分類過程中分別選取原始時間序列的歐式距離,以及SAX,SAX_TD[21],SAX_SM[22]和本文所提SAX_VOLA方法中給出的距離計算公式作為距離度量的依據(jù)。其中,SAX_TD用每段子序列的起止點與均值的差構建趨勢距離,進而計算出兩條時間序列符號化后的距離;SAX_SM直接用子序列之間起止點的差構建趨勢距離,進而求得兩條時間序列符號化之后的距離。
除了前述不同的時間序列距離計算方法,符號化表示時字符集的大小α和時間序列的分段數(shù)w也會對分類的準確率產(chǎn)生影響。為減小這些參數(shù)對評估結果的干擾,對各種不同的距離計算方法進行多次實驗,設α=3,…,10,w=2,…,n/2(n為時間序列的長度,w每次取值為前一次的2倍),選擇分類準確率最高的結果作為相應符號聚合近似表示的最終實驗結果(如果不同參數(shù)值取得了相同的分類準確率,則選w較小的參數(shù))。表3所示為不同方法在不同數(shù)據(jù)集上得到的分類準確率,其中Eucild表示歐式距離,加粗表示各數(shù)據(jù)集上所取得的最好分類準確率。
表3 各種方法在不同數(shù)據(jù)集上的分類準確率
續(xù)表3
由表3可見,本文所提SAX_VOLA及其對應的時間序列距離度量方法,在19組數(shù)據(jù)集上比直接用原始時序數(shù)據(jù)的歐式距離進行分類的效果好,在制造業(yè)相關的Trace數(shù)據(jù)集上提升效果最好,達到17%,在Yoga數(shù)據(jù)集上提升最低,僅為0.1%;在19組數(shù)據(jù)集上優(yōu)于SAX,在OliveOil數(shù)據(jù)集上提升較大,高達73.3%,在Wafer數(shù)據(jù)集上提升較小,只有0.1%;在17組數(shù)據(jù)集上優(yōu)于SAX_TD,1組數(shù)據(jù)集上持平,在Trace數(shù)據(jù)集上提升最高,為18%,在Wafer數(shù)據(jù)集上提升最小,為2%;在18組數(shù)據(jù)集上略優(yōu)于SAX_SM方法,在Trace數(shù)據(jù)集上提升最高,為22.7%,在Wafer數(shù)據(jù)集上提升最小,僅1%;在17組數(shù)據(jù)集上優(yōu)于所有方法。具體對比如圖3所示,可見基于SAX_VOLA的時間序列分類效果在本文所列多數(shù)數(shù)據(jù)集上優(yōu)于傳統(tǒng)的歐氏距離、SAX、SAX_TD和SAX_SM方法。
由圖3和表3可見,本文所提方法在Coffee數(shù)據(jù)集上比其他方法的分類準確率略低。該數(shù)據(jù)集中時間序列的圖像如圖4所示(圖中橫坐標表示時間序列的長度,縱坐標表示時間序列上每個點的對應值),可見該時間序列局部波動不顯著。作為對比,圖5和圖6給出了SAX_VOLA表現(xiàn)較好的Wafer和Trace數(shù)據(jù)集中的時間序列圖像,其最明顯的特點就是局部波動明顯,或者有明顯的上升或下降趨勢,與前面提及的波動率能夠有效捕捉子序列的上升/下降趨勢信息及波動幅度信息的事實相符。一般而言,當原始時間序列有明顯的上升、下降趨勢或局部波動幅度較大時,本文所提SAX_VOLA的效果更好。這種現(xiàn)象在制造業(yè)領域生產(chǎn)線的監(jiān)測數(shù)據(jù)上較常見,因此本文方法在該類情況下會取得更準確的分析結果。
下面分析本算法的時間性能。實驗所用計算機的配置為CPU i5-4200M、8 G內存、Windows操作系統(tǒng)。圖7所示為前述4種符號化方法在ECG200,GunPoint,OliveOil,Trace數(shù)據(jù)集上的時間成本(包括符號化降維和相似度計算)對比圖,可見本文所提SAX_VOLA的時間成本略高于其他3種方法,原因是本文提出的波動率指標計算復雜度高于始末距離、趨勢距離和SAX無附加指標的計算復雜度,這也是為了抽取更多序列信息而付出的代價。然而,隨著分段數(shù)的增加,由于段內采樣點的減少會降低波動率計算的代價,SAX_VOLA與SAX_TD的時間成本差距逐漸縮小。當然,無論本文所提SAX_VOLA還是傳統(tǒng)的符號化方法,由于取值點大大減少,在進行時間序列相似度計算和分類等任務時,時間效率均明顯優(yōu)于未符號化時對原始時間序列的處理效率。
CWRU滾動軸承數(shù)據(jù)中心的軸承故障診斷數(shù)據(jù)集[29]是世界公認的軸承診斷標準數(shù)據(jù)集之一,為進一步驗證SAX_VOLA在工業(yè)領域的應用,本次實驗選取CWRU軸承數(shù)據(jù)中的驅動端數(shù)據(jù)。被診斷軸承型號為深溝球軸承SKF6205,系統(tǒng)采樣頻率為12 kHz,負載為1 HP,選取兩個周期。被診斷軸承的缺陷位置有滾動體損傷、外圈損傷和內圈損傷3種,損傷直徑分別為0.007 inch,0.014 inch,0.02 inch,共9種損傷狀態(tài),加上正常狀態(tài)共計10種狀態(tài),每種狀態(tài)選取100個樣本,然后隨機挑選30%的數(shù)據(jù)作為測試集,字符集α和w的設置與各種方法在UCR時間序列檔案庫中數(shù)據(jù)集上的設置相同。表4所示為各種方法在軸承故障數(shù)據(jù)集上的最高分類準確率與取得最高分類準確率所耗費的時間,可見SAX_VOLA取得的分類準確率略優(yōu)于其他方法,時間成本比略低于SAX_TD。導致這一現(xiàn)象的原因與3.1節(jié)相同,為了保證時間序列分類的準確性,需要在符號化過程中盡量抽取原始時間序列更多的信息,因此將耗費更多的計算時間。
表4 各種方法在軸承故障數(shù)據(jù)集上的分類準確率和耗時
針對制造業(yè)等領域的時間序列數(shù)據(jù)降維問題,本文提出一種新的時間序列符號聚合近似方法,通過引入波動率指標同時量化時間序列的局部波動幅度和變化趨勢等信息,彌補了傳統(tǒng)SAX在符號化過程中丟失的波動信息。實驗結果表明,當時間序列有明顯的上升和下降趨勢,或者局部存在頻繁的波動時,本文方法在時間序列分類問題上的準確率上通常優(yōu)于傳統(tǒng)方法。
然而,在時間序列的降維過程中,本文采用等長分割的辦法,可能導致識別出的序列模式不準確,后續(xù)工作可以嘗試根據(jù)波動率對時間序列進行不等長分割,進而更加精確地捕捉時間序列模式和狀態(tài)的變化信息。另外,時間序列分析已經(jīng)應用于智能制造[30-31]、業(yè)務過程管理[32-34]等諸多領域,如何在這些領域推廣和應用本文方法也是下一步研究的重點。