基于二維奇異值分解的多元時間序列相似匹配方法

吳虎勝張鳳鳴鐘 斌

①(空軍工程大學裝備管理與安全工程學院 西安 710051)②(武警工程大學裝備工程學院 西安 710086)

1 引言

多元時間序列類型的數(shù)據(jù)在現(xiàn)實中廣泛存在，其相似性研究的應用領域也十分廣泛，如水文數(shù)據(jù)分析[1]，飛行質(zhì)量評估[2]，圖像匹配[3]，動物行為分析[4]，腦電圖分析[5]等。因此，研究MTS的相似模式匹配具有重要的現(xiàn)實意義和廣泛的應用前景。

當前的相關(guān)研究主要針對一元時間序列，產(chǎn)生了一些比較好的理論和方法并獲得了廣泛應用[6,7]。多元時間序列(Multivariate Time Series, MTS)由于其高維、稀疏、變量間關(guān)系復雜等特點，使得對其進行相似模式匹配的研究有一定困難。現(xiàn)常用的方法有直接歐氏距離法(Eculid Distance, ED)，主成分分析法(Principal Component Analysis, PCA)，基于點分布特征(Point Distribution, PD)的匹配方法，動態(tài)時間彎曲距離法(Dynamic Time Warping,DTW)等。但很多研究已證明 ED方法的魯棒性不好，即對時間序列在垂直和水平方向上的波動的魯棒性都不好，對復雜的 MTS的描述能力有限[8]；PCA方法簡單且無參數(shù)限制，對大規(guī)模的MTS分析具有一定優(yōu)勢[9]。但 PCA 方法的關(guān)鍵在于求取MTS的有效主成分分量，而這通常要求要有足夠的樣本點才行[10]。另外，PCA方法在計算兩主成分的夾角余弦2

cos β時并未考慮兩者的正負方向(例如，夾角60β=°和120β=°時將會得到相同的結(jié)果)，從而導致誤判；PD方法從MTS的形狀特征出發(fā)，通過提取MTS的極大值、極小值等9個局部重要點對其進行模式表示，并借助Eculid距離實現(xiàn)了相似性度量，對于如Robot這樣的小規(guī)模數(shù)據(jù)集的匹配具有一定優(yōu)勢[11]。但該方法未考慮MTS不同變量間量綱和特征差異，且若數(shù)據(jù)集中各MTS樣本的形狀差異較小時其匹配結(jié)果也不一定可靠；DTW 距離經(jīng)擴展后可用于MTS的相似性度量[12]，可支持不同時間跨度的時間序列間的相似性度量，支持時間軸的伸縮和彎曲，但其計算復雜度較高，不能適應大規(guī)模數(shù)據(jù)集的相似匹配[13]。文獻[14]基于DTW距離提出了一種將趨勢距離(Trend Distance, TD)用于MTS的相似匹配，該方法對不同規(guī)模的數(shù)據(jù)集都具有一定匹配能力且相對于DTW距離其計算效率也有所提高，但該方法計算復雜度仍然較大，實際耗時較長。

本文基于2維奇異值分解構(gòu)建了MTS的低階近似特征矩陣作為MTS的模式表示，而后結(jié)合Euclid距離實現(xiàn)了序列的相似模式匹配，并通過多組實驗對比詳細論述了該方法的有效性。

2 基于2DSVD的MTS相似匹配

2.1 多元時間序列

現(xiàn)實中，MTS以各種形式廣泛存在，如腦電圖數(shù)據(jù)(EEG)，飛行數(shù)據(jù)等。這里將變量數(shù) n大于 1的按時間先后順序記錄的值(其中稱為 MTS。 t表示時刻，為第i個變量在t時刻的記錄值。在數(shù)據(jù)挖掘領域，非正則化的數(shù)據(jù)之間比較的意義不大，由于MTS各變量間量綱、特征等方面的差異尤其如此，這里給出同構(gòu)的MTS的定義以限定研究對象。

定義1 同構(gòu)的MTS 同構(gòu)的MTS需要滿足以下條件：(1)對于MTS數(shù)據(jù)集，各序列樣本的變量維數(shù)相同，變量之間一一對應且表示相同的含義；(2)對于某一MTS樣本，各變量數(shù)值的記錄時刻對應；(3)MTS需經(jīng)正則化處理。

這里采用Max-min正則化方法將MTS各變量的值都映射到范圍內(nèi)，變量i的值v轉(zhuǎn)換后得到的值V為

2.2 2維奇異值分解

2維奇異值分解(2D Singular Value Decomposition, 2DSVD)是經(jīng)典奇異值分解的拓展，作為一種很好的低階近似方法，2DSVD可清楚地反應諸如 2維圖片序列，2維衛(wèi)星云圖序列等 2維物體的本質(zhì)[15]。通過 2DSVD 所構(gòu)造的特征是 2維的矩陣而并非 1維的特征向量?？衫迷?MTS樣本直接構(gòu)造其行-行，列-列協(xié)方差矩陣并計算其特征向量實現(xiàn)對 MTS的特征提取。這里首先給出2DSVD的簡單描述：

由上述求解過程可以看出，jM 是從數(shù)據(jù)集整體的角度所得的單個多元時間序列 Sj的特征，反應的是Sj相對于其它MTS的特征，更能體現(xiàn)其區(qū)別于其它MTS的本質(zhì)所在。

2.3 匹配方法及其復雜度分析

基于2DSVD可以得到兩個多元時間序列Sp和Sq之間的距離由如式(6)所示的最佳低階近似表示。

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

表1 基于2DSVD的MTS相似匹配算法流程

對該相似匹配算法進行時間復雜度分析。由文獻[16]可知，對一個aa×的矩陣進行奇異值分解求取其主特征向量的復雜度為由于，其中m為MTS的長度或觀測值數(shù)量，n為MTS的變量數(shù)，則步驟(2)~步驟(4)的時間復雜度為；步驟(5)為提取各MTS的模式表示，其復雜度為，其中c為待匹配的 MTS數(shù)據(jù)集樣本總數(shù)，r為行-行協(xié)方差矩陣 F的主特征向量數(shù)目，s為列-列協(xié)方差矩陣G的主特征向量數(shù)目；步驟(6)，步驟(7)執(zhí)行k-近鄰查詢找出和參考多元時間序列SR最相似的k個MTS的時間復雜度為。因此，算法的復雜度為。一般而言，n,k,r,s相對 m和c而言較小，因此算法的復雜度主要取決于單個MTS的序列長度m和MTS數(shù)據(jù)集樣本總數(shù)c。

3 實驗與分析

3.1實驗方法、數(shù)據(jù)及環(huán)境

為分析各方法性能，采用k-近鄰方法進行實驗，具體描述如下：假定數(shù)據(jù)集中含有n個MTS，任意抽取一個作為輸入樣本X，找出與X最相似的“k個樣本”，一般k取10,5,1。統(tǒng)計找出的“k個樣本”中與X同類別的樣本數(shù)量0n，即可計算準確率。其它序列依次作為輸入樣本，重復以上實驗，可得到n個準確率。如將準確率視為離散隨機變量ε，則準確率的數(shù)學期望*e可按式(8)確定。

選取不同數(shù)據(jù)規(guī)模、已知分類結(jié)果的3組MTS數(shù)據(jù)集作為實驗對象：Robot Execution Failure(記為 REF)數(shù)據(jù)集[17]，Wafer(記為 WA)數(shù)據(jù)集[18]和Electroencephalography(記為 EEG)數(shù)據(jù)集[19]。

實驗環(huán)境：Windows XP系統(tǒng)，CPU 2.00 GHz，內(nèi)存2 G，算法采用Matlab 2008a平臺下的M語言實現(xiàn)。

3.2 相似性匹配方法比較

運用上述實驗方法，采用ED, PCA, PD, TD和2DSVD共5種方法分別對REF, Wafer和EEG 3種不同規(guī)模的數(shù)據(jù)集進行相似性匹配實驗。

現(xiàn)先以Wafer數(shù)據(jù)集為例進行說明，它是一組對半導體加工設備實時監(jiān)控的數(shù)據(jù)集。共采用6個不同的真空傳感器采集數(shù)據(jù)，整個數(shù)據(jù)集分為兩類，normal和abnormal，分別包含1067和127個樣本，隨機抽取其中200個正常和100個不正常樣本共300個樣本進行試驗，樣本的時間跨度為104~198，屬多變量、不等時間跨度的中等規(guī)模的MTS數(shù)據(jù)集。將 300個樣本依次作為輸入樣本進行相似模式匹配。試驗中 PD方法的分割形式為，模式向量則以極小值，5%, 10%, 25%,50%, 75%, 90%, 95%及極大值共9個分位點特征來構(gòu)建[11]；PCA方法提取其前3個主成分作為MTS的模式表示(3個主元的貢獻率Q＆gt;85%)；TD方法中設各變量維度上σ和ω相等，傾斜角差異和時間跨度差異的權(quán)重值分別為0.8ε=,0.2λ=[14];2DSVD方法所涉及到的參數(shù)設置為r=33, s=4。另外，采用ED, 2DSVD和PCA方法計算時，對樣本進行截取使得MTS的時間跨度皆為104，當然也可以采用支持向量機等方法預測使得各 MTS的時間跨度相同。執(zhí)行k近鄰查詢并統(tǒng)計不同準確率上的查詢結(jié)果數(shù)量，結(jié)果如圖1所示。

總體來看，在匹配正確率為小概率事件的情況下(比如e取值為0,0.1,0.2等)，2DSVD和PCA方法，尤其是 2DSVD方法所對應次數(shù)都比其它幾種方法要少；而在正確率為大概率事件時，特別當e=1時，2DSVD方法對應的次數(shù)要比其它幾種方法要多?？梢姡瑢Υ朔N類型的數(shù)據(jù)，2DSVD方法在準確率分布方面具有優(yōu)勢。其次，表2為采用5種方法依據(jù)k近鄰實驗方法分別計算得到的相似匹配準確率的數(shù)學期望和平均計算耗時。對于Wafer數(shù)據(jù)集這種中等規(guī)模的MTS數(shù)據(jù)集，2DSVD, TD和PCA，特別是 2DSVD方法能取得較好的處理效果，其中2DSVD方法的準確率最高，達到了93%左右且計算耗時也相對較短。

圖1 5種模式匹配方法在Wafer數(shù)據(jù)集中的實驗結(jié)果

表2 3種數(shù)據(jù)集的相似匹配準確率的數(shù)學期望e*(%)及平均計算耗時t(ms)

而PD方法的匹配效果卻較差，準確率期望僅在57%左右。為找出原因，隨機選取Wafer數(shù)據(jù)集中第90個樣本(記為WA-90)作為輸入樣本，分別使用5種方法進行模式匹配得到最相似的樣本。如圖2所示，得到各匹配樣本與原樣本在形狀上都具有較好的相似性。筆者又對比了數(shù)據(jù)集中normal類和abnormal類多個樣本，這些樣本都具有較為相似的形狀特征，僅有一些細微差異。PD方法是從MTS的3維空間描述出發(fā)，提取其極大值、極小值等共9個局部形狀重要點來構(gòu)建特征模式向量進而刻畫MTS的特征，是基于MTS的形狀概貌的一種相似匹配方法。而由圖2可見，wafer數(shù)據(jù)集中的各樣本在形狀上很相似，以至于采用PD方法進行相似匹配時出現(xiàn)了較多誤判使得匹配效果較差。

REF數(shù)據(jù)集是對 Robot進行狀態(tài)監(jiān)測的數(shù)據(jù)集，從其5個子數(shù)據(jù)集選取LP1作為實驗數(shù)據(jù)集，LP1共包含 6個變量，分為 nomal, collision, fr_collision和obstruction共4類，共88個樣本，時間跨度均為15，屬于變量維數(shù)低，等時間跨度的小規(guī)模數(shù)據(jù)集；試驗中PD方法的分割形式為:，模式向量構(gòu)建同上，PCA方法提取其前2個主成分作為MTS的模式表示(前2個主元的貢獻率 Q＆gt;85%)；TD方法涉及的參數(shù)同上；2DSVD方法所涉及到的參數(shù)r=3,s=4。如表2所示，對于LP1數(shù)據(jù)集，TD和PCA方法的處理效果不佳，匹配的準確率期望值在55%左右。而2DSVD方法的處理效果雖然優(yōu)于PCA方法，卻不及PD方法。從圖3也可看出，對于此種類型的數(shù)據(jù)，PD方法所得結(jié)果在準確率分布方面具有優(yōu)勢。

為分析原因，隨機選取第 66個樣本(屬obstruction類，記為LP1-66)作為輸入樣本，分別使用5種方法進行模式匹配得到最相似的樣本。如圖4所示，這些MTS樣本都具有不同程度的形狀特征差異且時間序列長度較短，PD方法就是基于MTS的局部形狀重要點的分布特性的，因此對于樣本點較少，突顯局部重要點特征的小規(guī)模MTS數(shù)據(jù)集處理效果最好；TD方法利用擬合線段的傾斜角和時間跨度作為特征，體現(xiàn)的是序列連續(xù)變化趨勢而非狀態(tài)點的具體數(shù)值，因而無法很好地刻畫出時間跨度較小，只體現(xiàn)某些狀態(tài)點的 MTS數(shù)據(jù)的特征，所以匹配效果不佳；PCA方法是基于統(tǒng)計方式的一種相似匹配方法，通常要求足夠的樣本點才能有效求解得到主成分向量。如圖4所示，TD和PCA方法得到的參考樣本與輸入樣本的形狀差異較大，匹配中存在一定的誤判風險，特別對LP1這樣的樣本點數(shù)據(jù)較少的小規(guī)模 MTS數(shù)據(jù)集的相似性模式匹配已不是一種適宜的方法。而 2DSVD方法的特征提取是基于MTS的行-行和列-列協(xié)方差矩陣的主特征向量的，體現(xiàn)的是序列矩陣的2維整體內(nèi)在特征。因此，2DSVD方法對此類數(shù)據(jù)集仍具有一定刻畫能力，所得到的匹配樣本與輸入樣本在形狀上也較為相似。

圖2 Wafer數(shù)據(jù)集的相似模式匹配結(jié)果

圖3 5種模式匹配方法在LP1數(shù)據(jù)集中的實驗結(jié)果

圖4 LP1數(shù)據(jù)集的相似模式匹配結(jié)果

EEG數(shù)據(jù)集是一組腦電圖數(shù)據(jù)。采用 256 Hz的電極同時在64個部位測得，數(shù)據(jù)來源于Alcoholic Subjects和Control Subjects兩種人群，共122個測試者，對每個測試者進行120次測試。不失一般性，實驗分別隨機選取編號為 co2a0000364和co2c0000337兩位實驗者的50次測試作為實驗數(shù)據(jù)集，共 100個樣本，時間跨度均為 256，屬變量維數(shù)高、時間跨度大的較大規(guī)模MTS數(shù)據(jù)集。試驗中PD方法的分割形式為，模式向量構(gòu)建同上；PCA方法提取前 20個主成分作為 MTS的模式表示(20個主元的貢獻率Q＆gt;85%)；TD方法的涉及參數(shù)同上；2DSVD方法所涉及到的參數(shù)r=26, s=5。執(zhí)行k近鄰查詢并統(tǒng)計查詢中不同準確率上的查詢結(jié)果數(shù)量，結(jié)果如圖5所示。總體來看，當匹配正確率為小概率事件的情況下，2DSVD和TD方法，尤其是2DSVD方法所對應次數(shù)都比其它幾種方法要少；而在正確率為大概率事件時，2DSVD方法對應的次數(shù)要比其它幾種方法要多，而ED的準確率分布卻剛好相反?？梢?，對此類型的數(shù)據(jù)，2DSVD方法所得的結(jié)果在準確率分布方面具有優(yōu)勢。

由表2可見，對于EEG數(shù)據(jù)集，從準確率期望來看，除ED和PD外，其它4種方法都能得到很好的相似模式匹配效果。由圖 6也可看出，TD,2DSVD, PCA方法匹配所得樣本與輸入樣本在形狀上比較相似，而ED和PD方法匹配所得樣本與輸入樣本在形狀上差異較大。主要由于ED逐點對齊匹配的方式對于EEG這樣較大規(guī)模的MTS數(shù)據(jù)集已不適用，誤差較大且耗時也相對較長；而PD方法的匹配效果相對另外3種方法又稍差。分析認為EEG的數(shù)據(jù)量較大且其3維形狀特征較為復雜，僅利用極大值點、極小值點等9個局部重要點來描述MTS的特征已顯不足，匹配中不可避免地出現(xiàn)誤判；對于PCA方法，由于樣本點數(shù)量較大，能較為有效地求取MTS的主成分，因而處理效果相對處理LP1數(shù)據(jù)集時要好，準確率期望稍好。TD方法能較好地體現(xiàn)序列的連續(xù)變化趨勢，匹配效果較好，準確率期望達到94%左右。但由于TD方法是基于DTW 的，雖進行了優(yōu)化但計算耗時還是相對其它幾種方法要長很多；2DSVD方法是從MTS數(shù)據(jù)集整體的角度，依據(jù)每個樣本矩陣的本質(zhì)特點計算得到其各自的低階近似，計算簡明，匹配效果最好，對于 EEG數(shù)據(jù)集匹配準確率期望達到 95%，以上且計算耗時也相對較短，較其它幾種方法具有明顯的優(yōu)勢。

綜上可見，2DSVD方法對不同規(guī)模的MTS數(shù)據(jù)集都具有一定的處理能力且計算耗時相對較短，尤其對于如EEG數(shù)據(jù)集這樣的多變量，等時間跨度的較大規(guī)模的MTS數(shù)據(jù)集匹配效果最佳，綜合性能較優(yōu)，優(yōu)于其它4種方法。如表3所示，本文從多個角度對5種方法進行了詳細的對比。

4 結(jié)束語

圖5 5種模式匹配方法在EEG數(shù)據(jù)集中的實驗結(jié)果

圖6 LP1數(shù)據(jù)集的相似模式匹配結(jié)果

表3 5種匹配方法對比

本文基于 2DSVD提出了一種新的，適用于MTS的相似模式匹配方法，通過從MTS數(shù)據(jù)集整體的角度分別計算MTS的行-行和列-列矩陣的協(xié)方差矩陣的主特征向量來構(gòu)建各 MTS的模式表示矩陣，并借助 Eculid距離實現(xiàn)相似性度量。通過與ED, PCA, PD, TD共4種方法對于3種不同規(guī)模的MTS數(shù)據(jù)集的仿真實驗例證了本文方法的有效性，并在試驗的基礎上從模式表示、相似性度量等4個方面總結(jié)了5種方法各自的特點。實驗表明本文方法充分利用了MTS的矩陣2維本質(zhì)特征，不受數(shù)據(jù)集中 MTS樣本的形狀特征限制，對不同規(guī)模的MTS數(shù)據(jù)集均有一定的刻畫能力，且計算耗時相對較短。尤其對如EEG數(shù)據(jù)集這樣的等時間跨度的較大規(guī)模MTS數(shù)據(jù)集的匹配具有明顯優(yōu)勢，準確率期望達到95%左右。但如何依據(jù)本文方法構(gòu)建數(shù)據(jù)索引以利于MTS的快速檢索與查詢，將是我們接下來要著力解決的問題和研究方向。

[1] 李士進, 朱躍龍, 張曉花, 等. 基于BORDA計數(shù)法的多元水文時間序列相似性分析[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] 毛紅保, 張鳳鳴, 馮卉, 等. 多元飛行數(shù)據(jù)相似模式查詢[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] 尹令, 洪添勝, 劉漢興, 等. 結(jié)構(gòu)相似子序列快速聚類算法及其在奶牛發(fā)情檢測中的應用[J]. 農(nóng)業(yè)工程學報, 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 ＆amp; 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] 劉博寧, 張建業(yè), 張鵬, 等. 基于曲率距離的時間序列相似性搜索方法[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 ＆amp; 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 ＆amp;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.