吳虎勝張鳳鳴鐘 斌
①(空軍工程大學裝備管理與安全工程學院 西安 710051)②(武警工程大學裝備工程學院 西安 710086)
多元時間序列類型的數據在現實中廣泛存在,其相似性研究的應用領域也十分廣泛,如水文數據分析[1],飛行質量評估[2],圖像匹配[3],動物行為分析[4],腦電圖分析[5]等。因此,研究MTS的相似模式匹配具有重要的現實意義和廣泛的應用前景。
當前的相關研究主要針對一元時間序列,產生了一些比較好的理論和方法并獲得了廣泛應用[6,7]。多元時間序列(Multivariate Time Series, MTS)由于其高維、稀疏、變量間關系復雜等特點,使得對其進行相似模式匹配的研究有一定困難。現常用的方法有直接歐氏距離法(Eculid Distance, ED),主成分分析法(Principal Component Analysis, PCA),基于點分布特征(Point Distribution, PD)的匹配方法,動態時間彎曲距離法(Dynamic Time Warping,DTW)等。但很多研究已證明 ED方法的魯棒性不好,即對時間序列在垂直和水平方向上的波動的魯棒性都不好,對復雜的 MTS的描述能力有限[8];PCA方法簡單且無參數限制,對大規模的MTS分析具有一定優勢[9]。但 PCA 方法的關鍵在于求取MTS的有效主成分分量,而這通常要求要有足夠的樣本點才行[10]。另外,PCA方法在計算兩主成分的夾角余弦2
cos β時并未考慮兩者的正負方向(例如,夾角60β=°和120β=°時將會得到相同的結果),從而導致誤判;PD方法從MTS的形狀特征出發,通過提取MTS的極大值、極小值等9個局部重要點對其進行模式表示,并借助Eculid距離實現了相似性度量,對于如Robot這樣的小規模數據集的匹配具有一定優勢[11]。但該方法未考慮MTS不同變量間量綱和特征差異,且若數據集中各MTS樣本的形狀差異較小時其匹配結果也不一定可靠;DTW 距離經擴展后可用于MTS的相似性度量[12],可支持不同時間跨度的時間序列間的相似性度量,支持時間軸的伸縮和彎曲,但其計算復雜度較高,不能適應大規模數據集的相似匹配[13]。文獻[14]基于DTW距離提出了一種將趨勢距離(Trend Distance, TD)用于MTS的相似匹配,該方法對不同規模的數據集都具有一定匹配能力且相對于DTW距離其計算效率也有所提高,但該方法計算復雜度仍然較大,實際耗時較長。
本文基于2維奇異值分解構建了MTS的低階近似特征矩陣作為MTS的模式表示,而后結合Euclid距離實現了序列的相似模式匹配,并通過多組實驗對比詳細論述了該方法的有效性。
現實中,MTS以各種形式廣泛存在,如腦電圖數據(EEG),飛行數據等。這里將變量數 n大于 1的按時間先后順序記錄的值(其中稱為 MTS。 t表示時刻,為第i個變量在t時刻的記錄值。在數據挖掘領域,非正則化的數據之間比較的意義不大,由于MTS各變量間量綱、特征等方面的差異尤其如此,這里給出同構的MTS的定義以限定研究對象。
定義1 同構的MTS 同構的MTS需要滿足以下條件:(1)對于MTS數據集,各序列樣本的變量維數相同,變量之間一一對應且表示相同的含義;(2)對于某一MTS樣本,各變量數值的記錄時刻對應;(3)MTS需經正則化處理。
這里采用Max-min正則化方法將MTS各變量的值都映射到范圍內,變量i的值v轉換后得到的值V為

2維奇異值分解(2D Singular Value Decomposition, 2DSVD)是經典奇異值分解的拓展,作為一種很好的低階近似方法,2DSVD可清楚地反應諸如 2維圖片序列,2維衛星云圖序列等 2維物體的本質[15]。通過 2DSVD 所構造的特征是 2維的矩陣而并非 1維的特征向量。可利用原 MTS樣本直接構造其行-行,列-列協方差矩陣并計算其特征向量實現對 MTS的特征提取。這里首先給出2DSVD的簡單描述:



由上述求解過程可以看出,jM 是從數據集整體的角度所得的單個多元時間序列 Sj的特征,反應的是Sj相對于其它MTS的特征,更能體現其區別于其它MTS的本質所在。
基于2DSVD可以得到兩個多元時間序列Sp和Sq之間的距離由如式(6)所示的最佳低階近似表示。


相似性的判別主要是基于距離的思想,即依據一定的距離公式(如歐氏距離,DTW距離,最長公共子串距離,編輯距離等)計算所比較兩序列之間的距離,認為距離越小則兩時間序列越相似,反之亦然。這里采用計算簡便、符合三角不等式的歐氏距離,如式(7)中表示歐氏范數。通過計算距離就可以度量兩序列Sp和Sq之間的距離,再借助于k-近鄰方法就可實現對MTS的相似匹配,具體算法流程如表1所示。

表1 基于2DSVD的MTS相似匹配算法流程
對該相似匹配算法進行時間復雜度分析。由文獻[16]可知,對一個aa×的矩陣進行奇異值分解求取其主特征向量的復雜度為由于,其中m為MTS的長度或觀測值數量,n為MTS的變量數,則步驟(2)~步驟(4)的時間復雜度為;步驟(5)為提取各MTS的模式表示,其復雜度為,其中c為待匹配的 MTS數據集樣本總數,r為行-行協方差矩陣 F的主特征向量數目,s為列-列協方差矩陣G的主特征向量數目;步驟(6),步驟(7)執行k-近鄰查詢找出和參考多元時間序列SR最相似的k個MTS的時間復雜度為。因此,算法的復雜度為。一般而言,n,k,r,s相對 m和c而言較小,因此算法的復雜度主要取決于單個MTS的序列長度m和MTS數據集樣本總數c。
為分析各方法性能,采用k-近鄰方法進行實驗,具體描述如下:假定數據集中含有n個MTS,任意抽取一個作為輸入樣本X,找出與X最相似的“k個樣本”,一般k取10,5,1。統計找出的“k個樣本”中與X同類別的樣本數量0n,即可計算準確率。其它序列依次作為輸入樣本,重復以上實驗,可得到n個準確率。如將準確率視為離散隨機變量ε,則準確率的數學期望*e可按式(8)確定。

選取不同數據規模、已知分類結果的3組MTS數據集作為實驗對象:Robot Execution Failure(記為 REF)數據集[17],Wafer(記為 WA)數據集[18]和Electroencephalography(記為 EEG)數據集[19]。
實驗環境:Windows XP系統,CPU 2.00 GHz,內存2 G,算法采用Matlab 2008a平臺下的M語言實現。
運用上述實驗方法,采用ED, PCA, PD, TD和2DSVD共5種方法分別對REF, Wafer和EEG 3種不同規模的數據集進行相似性匹配實驗。
現先以Wafer數據集為例進行說明,它是一組對半導體加工設備實時監控的數據集。共采用6個不同的真空傳感器采集數據,整個數據集分為兩類,normal和abnormal,分別包含1067和127個樣本,隨機抽取其中200個正常和100個不正常樣本共300個樣本進行試驗,樣本的時間跨度為104~198,屬多變量、不等時間跨度的中等規模的MTS數據集。將 300個樣本依次作為輸入樣本進行相似模式匹配。試驗中 PD方法的分割形式為,模式向量則以極小值,5%, 10%, 25%,50%, 75%, 90%, 95%及極大值共9個分位點特征來構建[11];PCA方法提取其前3個主成分作為MTS的模式表示(3個主元的貢獻率Q>85%);TD方法中設各變量維度上σ和ω相等,傾斜角差異和時間跨度差異的權重值分別為0.8ε=,0.2λ=[14];2DSVD方法所涉及到的參數設置為r=33, s=4。另外,采用ED, 2DSVD和PCA方法計算時,對樣本進行截取使得MTS的時間跨度皆為104,當然也可以采用支持向量機等方法預測使得各 MTS的時間跨度相同。執行k近鄰查詢并統計不同準確率上的查詢結果數量,結果如圖1所示。
總體來看,在匹配正確率為小概率事件的情況下(比如e取值為0,0.1,0.2等),2DSVD和PCA方法,尤其是 2DSVD方法所對應次數都比其它幾種方法要少;而在正確率為大概率事件時,特別當e=1時,2DSVD方法對應的次數要比其它幾種方法要多。可見,對此種類型的數據,2DSVD方法在準確率分布方面具有優勢。其次,表2為采用5種方法依據k近鄰實驗方法分別計算得到的相似匹配準確率的數學期望和平均計算耗時。對于Wafer數據集這種中等規模的MTS數據集,2DSVD, TD和PCA,特別是 2DSVD方法能取得較好的處理效果,其中2DSVD方法的準確率最高,達到了93%左右且計算耗時也相對較短。

圖1 5種模式匹配方法在Wafer數據集中的實驗結果

表2 3種數據集的相似匹配準確率的數學期望e*(%)及平均計算耗時t(ms)
而PD方法的匹配效果卻較差,準確率期望僅在57%左右。為找出原因,隨機選取Wafer數據集中第90個樣本(記為WA-90)作為輸入樣本,分別使用5種方法進行模式匹配得到最相似的樣本。如圖2所示,得到各匹配樣本與原樣本在形狀上都具有較好的相似性。筆者又對比了數據集中normal類和abnormal類多個樣本,這些樣本都具有較為相似的形狀特征,僅有一些細微差異。PD方法是從MTS的3維空間描述出發,提取其極大值、極小值等共9個局部形狀重要點來構建特征模式向量進而刻畫MTS的特征,是基于MTS的形狀概貌的一種相似匹配方法。而由圖2可見,wafer數據集中的各樣本在形狀上很相似,以至于采用PD方法進行相似匹配時出現了較多誤判使得匹配效果較差。
REF數據集是對 Robot進行狀態監測的數據集,從其5個子數據集選取LP1作為實驗數據集,LP1共包含 6個變量,分為 nomal, collision, fr_collision和obstruction共4類,共88個樣本,時間跨度均為15,屬于變量維數低,等時間跨度的小規模數據集;試驗中PD方法的分割形式為:,模式向量構建同上,PCA方法提取其前2個主成分作為MTS的模式表示(前2個主元的貢獻率 Q>85%);TD方法涉及的參數同上;2DSVD方法所涉及到的參數r=3,s=4。如表2所示,對于LP1數據集,TD和PCA方法的處理效果不佳,匹配的準確率期望值在55%左右。而2DSVD方法的處理效果雖然優于PCA方法,卻不及PD方法。從圖3也可看出,對于此種類型的數據,PD方法所得結果在準確率分布方面具有優勢。
為分析原因,隨機選取第 66個樣本(屬obstruction類,記為LP1-66)作為輸入樣本,分別使用5種方法進行模式匹配得到最相似的樣本。如圖4所示,這些MTS樣本都具有不同程度的形狀特征差異且時間序列長度較短,PD方法就是基于MTS的局部形狀重要點的分布特性的,因此對于樣本點較少,突顯局部重要點特征的小規模MTS數據集處理效果最好;TD方法利用擬合線段的傾斜角和時間跨度作為特征,體現的是序列連續變化趨勢而非狀態點的具體數值,因而無法很好地刻畫出時間跨度較小,只體現某些狀態點的 MTS數據的特征,所以匹配效果不佳;PCA方法是基于統計方式的一種相似匹配方法,通常要求足夠的樣本點才能有效求解得到主成分向量。如圖4所示,TD和PCA方法得到的參考樣本與輸入樣本的形狀差異較大,匹配中存在一定的誤判風險,特別對LP1這樣的樣本點數據較少的小規模 MTS數據集的相似性模式匹配已不是一種適宜的方法。而 2DSVD方法的特征提取是基于MTS的行-行和列-列協方差矩陣的主特征向量的,體現的是序列矩陣的2維整體內在特征。因此,2DSVD方法對此類數據集仍具有一定刻畫能力,所得到的匹配樣本與輸入樣本在形狀上也較為相似。

圖2 Wafer數據集的相似模式匹配結果

圖3 5種模式匹配方法在LP1數據集中的實驗結果

圖4 LP1數據集的相似模式匹配結果
EEG數據集是一組腦電圖數據。采用 256 Hz的電極同時在64個部位測得,數據來源于Alcoholic Subjects和Control Subjects兩種人群,共122個測試者,對每個測試者進行120次測試。不失一般性,實驗分別隨機選取編號為 co2a0000364和co2c0000337兩位實驗者的50次測試作為實驗數據集,共 100個樣本,時間跨度均為 256,屬變量維數高、時間跨度大的較大規模MTS數據集。試驗中PD方法的分割形式為,模式向量構建同上;PCA方法提取前 20個主成分作為 MTS的模式表示(20個主元的貢獻率Q>85%);TD方法的涉及參數同上;2DSVD方法所涉及到的參數r=26, s=5。執行k近鄰查詢并統計查詢中不同準確率上的查詢結果數量,結果如圖5所示。總體來看,當匹配正確率為小概率事件的情況下,2DSVD和TD方法,尤其是2DSVD方法所對應次數都比其它幾種方法要少;而在正確率為大概率事件時,2DSVD方法對應的次數要比其它幾種方法要多,而ED的準確率分布卻剛好相反。可見,對此類型的數據,2DSVD方法所得的結果在準確率分布方面具有優勢。
由表2可見,對于EEG數據集,從準確率期望來看,除ED和PD外,其它4種方法都能得到很好的相似模式匹配效果。由圖 6也可看出,TD,2DSVD, PCA方法匹配所得樣本與輸入樣本在形狀上比較相似,而ED和PD方法匹配所得樣本與輸入樣本在形狀上差異較大。主要由于ED逐點對齊匹配的方式對于EEG這樣較大規模的MTS數據集已不適用,誤差較大且耗時也相對較長;而PD方法的匹配效果相對另外3種方法又稍差。分析認為EEG的數據量較大且其3維形狀特征較為復雜,僅利用極大值點、極小值點等9個局部重要點來描述MTS的特征已顯不足,匹配中不可避免地出現誤判;對于PCA方法,由于樣本點數量較大,能較為有效地求取MTS的主成分,因而處理效果相對處理LP1數據集時要好,準確率期望稍好。TD方法能較好地體現序列的連續變化趨勢,匹配效果較好,準確率期望達到94%左右。但由于TD方法是基于DTW 的,雖進行了優化但計算耗時還是相對其它幾種方法要長很多;2DSVD方法是從MTS數據集整體的角度,依據每個樣本矩陣的本質特點計算得到其各自的低階近似,計算簡明,匹配效果最好,對于 EEG數據集匹配準確率期望達到 95%,以上且計算耗時也相對較短,較其它幾種方法具有明顯的優勢。
綜上可見,2DSVD方法對不同規模的MTS數據集都具有一定的處理能力且計算耗時相對較短,尤其對于如EEG數據集這樣的多變量,等時間跨度的較大規模的MTS數據集匹配效果最佳,綜合性能較優,優于其它4種方法。如表3所示,本文從多個角度對5種方法進行了詳細的對比。

圖5 5種模式匹配方法在EEG數據集中的實驗結果

圖6 LP1數據集的相似模式匹配結果

表3 5種匹配方法對比
本文基于 2DSVD提出了一種新的,適用于MTS的相似模式匹配方法,通過從MTS數據集整體的角度分別計算MTS的行-行和列-列矩陣的協方差矩陣的主特征向量來構建各 MTS的模式表示矩陣,并借助 Eculid距離實現相似性度量。通過與ED, PCA, PD, TD共4種方法對于3種不同規模的MTS數據集的仿真實驗例證了本文方法的有效性,并在試驗的基礎上從模式表示、相似性度量等4個方面總結了5種方法各自的特點。實驗表明本文方法充分利用了MTS的矩陣2維本質特征,不受數據集中 MTS樣本的形狀特征限制,對不同規模的MTS數據集均有一定的刻畫能力,且計算耗時相對較短。尤其對如EEG數據集這樣的等時間跨度的較大規模MTS數據集的匹配具有明顯優勢,準確率期望達到95%左右。但如何依據本文方法構建數據索引以利于MTS的快速檢索與查詢,將是我們接下來要著力解決的問題和研究方向。
[1] 李士進, 朱躍龍, 張曉花, 等. 基于BORDA計數法的多元水文時間序列相似性分析[J]. 水利學報, 2009, 40(3): 378-384.Li Shi-jin, Zhu Yue-long, Zhang Xiao-hua, et al.. BORDA counting method based similarity analysis of multivariate hydrological time series[J]. Journal of SHUILI, 2009, 40(3):378-384.
[2] 毛紅保, 張鳳鳴, 馮卉, 等. 多元飛行數據相似模式查詢[J].計算機工程與應用, 2011, 47(16): 151-155.Mao Hong-bao, Zhang Feng-ming, Feng Hui, et al..Similarity-based pattern querying in multivariate flight data[J]. Computer Engineering and Application, 2011, 47(16):151-155.
[3] 周瑜, 劉俊濤, 白翔. 形狀匹配方法研究與展望[J]. 自動化學報, 2012, 38(6): 889-910.Zhou Yu, Liu Jun-tao, and Bai Xiang. Research and perspective on shape matching[J]. Acta Automatica Sinica,2012, 38(6): 889-910.
[4] 尹令, 洪添勝, 劉漢興, 等. 結構相似子序列快速聚類算法及其在奶牛發情檢測中的應用[J]. 農業工程學報, 2012, 28(15):107-112.Yin Ling, Hong Tian-sheng, Liu Han-xing, et al..Subsequence clustering algorithm based on structural similarity and its application in cow estrus detection[J].Transaction of the Chinese Society of Agricultural Engineering, 2012, 28(15): 107-122.
[5] Hanna G, Marek D, Leszek K, et al.. Detection of similar sequences in EEG maps series using correlation coefficients matrix[J]. Machine Graphics & Vision International Journal,2011, 20(1): 73-92.
[6] Keogh E, Li W, Xi X P, et al.. Supporting exact indexing of arbitrarily rotated shapes and periodic time series under Euclidean and warping distance measures[J]. The VLDB Journal, 2009, 18(3): 611-630.
[7] 劉博寧, 張建業, 張鵬, 等. 基于曲率距離的時間序列相似性搜索方法[J]. 電子與信息學報, 2012, 34(9): 2200-2207.Liu Bo-ning, Zhang Jian-ye, Zhang Peng, et al.. Similarity search method in time series based on curvature distance[J].Journal of Electronics & Information Technology, 2012, 34(9):2200-2207.
[8] Li Hai-lin, Guo Chong-hui, and Qiu Wang-ren. Similarity measure based on piecewise linear approximation and derivative dynamic time warping for time series mining[J].Expert Systems with Applications, 2011, 38(12): 14732-14743.
[9] Yang K and Shahabi C. A pca-based similarity measure for multivariate time series[C]. Proceedings of the Second ACM International workshop on Multimedia Databases,Washington DC, USA, 2004: 65-74.
[10] Singhal A and Seborg D E. Pattern matching in multivariate time series databases using a moving window approach[J].Industrial and Engineering Chemistry Research, 2002, 41(16):3822-3828.
[11] 管河山, 姜青山, 王聲瑞. 基于點分布特征的多元時間序列模式匹配方法[J]. 軟件學報, 2009, 20(1): 67-79.Guan He-shan, Jiang Qing-shan, and Wang Sheng-rui.Pattern matching method based on point distribution for multivariate time series[J]. Journal of Software, 2009, 20(1):67-79.
[12] Zoltan B and Janos A. Correlation based dynamic time warping of multivariate time series[J]. Expert Systems with Applications, 2012, 39(17): 12814-12823.
[13] Stephan S, Brijnesh J, and Ernesto W D. Pattern recognition in multivariate time series[C]. Proceedings of the fourth workshop on workshop for Ph. D. students in information &knowledge management, Glasgow, UK, 2011: 27-34.
[14] 李正欣, 張鳳鳴, 李克武. 多元時間序列模式匹配方法研究[J].控制與決策, 2011, 26(4): 565-570.Li Zheng-xin, Zhang Feng-ming, and Li Ke-wu. Research on pattern matching method for multivariate time series[J].Control and Decision, 2011, 26(4): 565-570.
[15] Chris D and Ye J. Two-dimensional singular decomposition(2DSVD) for 2D maps and images[C]. Proceedings of the fifth IEEE International Conference on Data Mining. Houston,USA, 2005: 32-43.
[16] Brand M. Incremental singular value decomposition of uncertain data with missing values[C]. Proceedings of the 2002 European Conference on Computer Vision, Copenhagen,Denmark, 2002, 3: 1-12.
[17] Lopes L S and Camarinha L M. Robot execution failures[OL].http://kdd.ics.uci.edu/databases/robotfailure/robotfailure.html. 2012.
[18] Bobski. Wafer database [OL]. http://www.cs.cmu.edu/bobski.html. 2013.
[19] Begleiter H.EEG database [OL]. http://kdd.ics.uci.edu/databases/eeg/eeg.html. 2012.