基于模態(tài)重構(gòu)與多維評價的時間序列趨勢提取

杜加礎(chǔ)，車文剛，程文輝

(昆明理工大學(xué) 信息工程與自動化學(xué)院，昆明 650500)

0 引 言

時間序列是按照時間順序排列、隨時間變化且相互關(guān)聯(lián)的有序數(shù)據(jù)集合。隨著信息技術(shù)的不斷發(fā)展和進步，時間序列的數(shù)據(jù)量和復(fù)雜度與日俱增，并存在于社會經(jīng)濟生活的各個方面。時間序列的主要研究方向有數(shù)據(jù)預(yù)處理、數(shù)據(jù)分割、相似度度量、特征提取、維度規(guī)約、分類、聚類和規(guī)則識別等，對于時間序列的研究成果廣泛運用于電信、網(wǎng)絡(luò)安全、醫(yī)療衛(wèi)生、電子商務(wù)和工業(yè)工程等領(lǐng)域。

趨勢特征是時間序列的一個重要信息，它直觀反映了時間序列的變化趨勢和模式特征，對于時間序列的特征提取一直是時間序列研究的一個重要方向。為提取時間序列的趨勢特征，人們提出了時間序列特征表示方法，能夠?qū)崿F(xiàn)對時間序列的降維，并保留其主要的趨勢特征。目前已經(jīng)提出的表示方法主要分為4類[1]：①基于頻域變換的表示方法，主要有離散傅里葉變換[2-3]和離散小波變換[4-6]；②基于模型的表示方法，主要有自回歸模型[7]、隱馬爾科夫模型[8]和核模型[9]；③基于符號的表示方法，其中使用最廣泛的是符號聚合近似方法[10-11]；④基于分段的特征表示方法，主要有分段聚合近似(piecewise aggregate approximation, PAA)[12]和分段線性表示(piecewise linear representation, PLR)[13-14]。其中PLR方法因其具有簡潔直觀、處理效率較高、產(chǎn)生的表示結(jié)果更符合人類視覺體驗等優(yōu)點，得到了廣泛的應(yīng)用。傳統(tǒng)的PLR方法主要有自頂向下(top-down, TD)的PLR-TD方法[15]、自底向上(bottom-up, BU)的PLR-BU方法[16]、基于滑動窗口(sliding window, SW)的PLR-SW方法[17]和基于重要點(important points segment, IPS)的PLR-IPS方法[18]。為了解決傳統(tǒng)PLR方法分段擬合難以達到理想效果的問題，近年來出現(xiàn)了不少改進的PLR方法，例如文獻[19]提出的基于轉(zhuǎn)折點和趨勢段(turning point and trend segment, TPTS)的分段線性表示方法(PLR-TPTS)、文獻[20]提出的基于重要點雙重評價(double Evaluation, DE)的分段線性表示方法(PLR-DE)等，但是，這些方法仍然存在錯誤提取分段點、擬合誤差較大、難以反映整體趨勢等問題。

文獻[21]提出了經(jīng)驗?zāi)B(tài)分解(empirical mode decomposition, EMD)，該方法主要用于處理非線性、非穩(wěn)態(tài)的數(shù)據(jù)。近年來，EMD逐漸從信號分析處理領(lǐng)域被引入到對時間序列的分析中，但是EMD存在模態(tài)混疊的問題。文獻[22]提出了集合經(jīng)驗?zāi)B(tài)分解(ensemble empirical mode decomposition, EEMD)，通過對原始信號多次加入不同的白噪聲來消除EMD出現(xiàn)的模態(tài)混疊現(xiàn)象。文獻[23]提出了自適應(yīng)噪聲的完備經(jīng)驗?zāi)B(tài)分解(complete ensemble empirical mode decomposition with adaptive noise，CEEMDAN)方法，通過加入自適應(yīng)噪聲，解決了EEMD重構(gòu)誤差不為零、分解效率低等問題，且進一步減小了模態(tài)混疊效應(yīng)。

本文在PLR-IPS的基礎(chǔ)上，將模態(tài)分解方法和模態(tài)重構(gòu)思想引入時間序列的分段線性表示研究中，提出一種基于CEEMDAN分解重構(gòu)與多維評價的重要點評價策略，并結(jié)合智能尋優(yōu)算法提出一種新的時間序列趨勢提取算法。新方法從單點、局部和整體3個維度綜合評價重要點，通過模態(tài)重構(gòu)思想選取全局性分段點，克服了局部分析的缺陷；通過對重要點的分類討論，使得對重要點的評價更為精準(zhǔn)。同時，本文方法使用智能尋優(yōu)算法，避免了基于經(jīng)驗選擇參數(shù)而導(dǎo)致的擬合誤差上下波動問題。

1 基于PLR-IPS的趨勢提取方法

1997年，文獻[14]將分段線性表示算法引入時間序列數(shù)據(jù)挖掘領(lǐng)域；2002年，文獻[18]提出了基于重要點的PLR-IPS方法，它能夠直觀反映時間序列的變化趨勢，且具有線性的時間復(fù)雜度和較高的精確度。

時間序列S是由n項數(shù)據(jù)值按照時間上的先后順序組成的有限集合，可以表示為

S={(x1,t1),…,(xi,ti),…,(xn,tn)}

(1)

(1)式中，(xi,ti)為時間序列中的第i個元素，表示ti時刻的數(shù)據(jù)值為xi，且i

(2)

{xpq-1≤xpq}∩{xpq

{xpq-1

{xpq-1≥xpq}∩{xpq>xpq+1}∪

{xpq-1>xpq}∩{xpq≥xpq+1}

(3)

此外，規(guī)定時間序列的兩個端點為重要點，由(3)式得到m個重要點，則重要點序列可以表示為

(4)

由上述的定義可知，本文的重要點序列由時間序列的端點和內(nèi)部的離散局部極值點組成。

(5)

(5)式中，d為分段點數(shù)量，且一般情況下有d

給定原始時間序列S[n]及其對應(yīng)的PLR序列SP[d]，則時間序列分段線性表示后的壓縮率C可以表示為

(6)

(7)

基于PLR-IPS的方法對時間序列趨勢提取的目標(biāo)就是在指定壓縮率的情況下，使得分段表示后的擬合誤差達到最低，減小擬合誤差的關(guān)鍵在于采用合理的評價策略從重要點序列中選取分段點。

2 模態(tài)重構(gòu)及多維評價

2.1 CEEMDAN分解與IMF重構(gòu)

第1步使用EMD對S[n]+ε0ωi[n]進行I次分解得到I個模態(tài)分量，則CEEMDAN分解產(chǎn)生的第1個模態(tài)分量(k=1)可以表示為

(8)

第2步計算得到第1個殘差為

(9)

第3步使用EMD對r1[n]+ε1E1(ωi[n])進行I次分解得到I個模態(tài)分量，則CEEMDAN分解產(chǎn)生的第2個模態(tài)分量(k=2)可以表示為

(10)

第4步同理，對于k=2,…,K，其對應(yīng)的第k個殘差可以表示為

(11)

第5步使用EMD對rk[n]+εkEk(ωi[n])進行I次分解得到I個模態(tài)分量，則CEEMDAN分解產(chǎn)生的第k+1個模態(tài)分量可以表示為

(12)

第6步重復(fù)執(zhí)行第4步和第5步，直到得到的殘差不能再被分解，則最終的殘差可以表示為

(13)

在得到CEEMDAN方法的所有K個模態(tài)分量后，給定的原始時間序列可以表示為

(14)

算法實現(xiàn)中，每次實驗中的系數(shù)εk能夠允許在每一個模態(tài)分解階段選擇合適的信噪比。

實驗表明，CEEMDAN能精確重構(gòu)原始信號，在EEMD和EMD基礎(chǔ)上重構(gòu)誤差顯著減小，故本文選用CEEMDAN對時間序列進行分解和重構(gòu)。

設(shè)置信噪比ε=0.2，噪聲添加次數(shù)I=500，對含有400個數(shù)據(jù)的示例時間序列S[n]進行CEEMDAN分解后得到8個模態(tài)分量IMF和最終殘差RES，CEEMDAN分解的結(jié)果如圖1所示。

2007年，文獻[24]在使用EEMD研究原油價格變化時，發(fā)現(xiàn)IMF重構(gòu)后的序列能夠反應(yīng)原油價格的關(guān)鍵轉(zhuǎn)折點和整體變化趨勢?；谶@個研究成果，本文使用CEEMDAN方法對時間序列進行分解，并探討重構(gòu)IMF數(shù)量對重構(gòu)序列的影響。從低頻到高頻每次多加入一個IMF進行多次重構(gòu)，可以得到多個重構(gòu)的時間序列。IMF多次重構(gòu)的結(jié)果如圖2所示。

圖1 CEEMDAN分解結(jié)果Fig.1 Results of CEEMDAN decomposition

圖2表明，從低頻IMF開始，每增加一個更高頻的IMF參與重構(gòu)，重構(gòu)序列就能保留更多原始時間序列的細節(jié)；每減少一個高頻的IMF，重構(gòu)序列就能去除更多小幅波動而能夠更加反映原始時間序列的整體變化趨勢。將參與重構(gòu)的IMF數(shù)量定義為重構(gòu)參數(shù)θ，下面給出模態(tài)重構(gòu)序列的定義。

(15)

2.2 基于模態(tài)重構(gòu)序列的整體評價因子

(16)

圖2 IMF多次重構(gòu)的結(jié)果Fig.2 Results of multiple IMFs reconstructions

但是CEEMDAN分解重構(gòu)后會造成一些重要點的漂移，故需要進行漂移點修正。對比重構(gòu)重要點序列SRZ[mθ]和原始時間序列重要點序列SZ[m]，若重構(gòu)重要點序列中的重要點不存在于重要點序列中，則進行漂移修正；若重構(gòu)重要點為極大值重要點，則將此重構(gòu)重要點修正為與其距離最近的SZ[m]中的極大值重要點；若為極小值重要點則修正為與其距離最近的SZ[m]中的極小值重要點。在對重構(gòu)重要點序列進行漂移點修正后，所有重構(gòu)重要點均能與原始時間序列重要點對應(yīng)，所以有

(17)

(18)

(18)式中，當(dāng)θ

根據(jù)定義4可知，全局因子只有K個值，只能將m個重要點評價為?K/2」個重要等級，無法精確評價各個重要點的重要程度，本文繼續(xù)給出在單點和局部維度上對重要點進行評價的因子。

2.3 單點評價因子與局部評價因子

根據(jù)(3)式，重要點可以分為極大值重要點和極小值重要點。圖3為極大值重要點的AL與AR示意圖。圖3中，xpq-1、xpq+1為極大值重要點xpq的相鄰極小值重要點，xpq-2、xpq+2為相鄰極大值重要點，AL與AR為重要點與其相鄰2個重要點構(gòu)成的三角形的面積。重要點對時間序列趨勢的影響越大，其與相鄰重要點的差異程度越大。即重要點xpq與左邊的相鄰重要點xpq-1差異越大，則AL越大；與右邊的相鄰重要點xpq+1差異越大，則AR越大，據(jù)此給出特征因子的定義。

圖3 極大值重要點的AL與AR示意圖Fig.3 Illustration for ALand ARof maximum important point

1)當(dāng)q=1,m時有

FC(q)=max(FC)

(19a)

2)當(dāng)q=2時有

(19b)

3)當(dāng)q=m-1時有

(19c)

4)當(dāng)q=3,4,…,m-2時有

(19d)

特征因子是聚焦于重要點本身特征的因子，極大值重要點特征因子的4種類型如圖4所示。

1)當(dāng)q=1,m時有

FC(q)=max(FC)

(20a)

2)當(dāng)q=2時有

(20b)

3)當(dāng)q=m-1時有

(20c)

4)當(dāng)q=3,4,…,m-2時有

FC(q)=

(20d)

本文特征因子參考了文獻[25]中的三角形中心距離，但是上述方法均只采用了相鄰兩點來度量重要點的差異程度。文獻[20]中的“距離因子”只考慮了相鄰的2個重要點，本文的特征因子將其擴展到相鄰2個極大值重要點和2個極小值重要點并進行分類討論，對重要點的評價結(jié)果更為可靠。下面給出邊界因子評價重要點對局部趨勢的影響范圍。

圖4 極大值重要點特征因子的4種類型Fig.4 Four types of characteristic factor of maximum important point

(21)

(xpq-l

(xpq+r

則有BL=pq-l，BR=pq+r，邊界因子著重于重要點的局部影響力范圍，可以理解為該重要點的“勢力范圍”。

圖5為極大值重要點xpq的邊界因子示意圖，左邊界點滿足條件xpq-l+1

特征因子和邊界因子分別反映了重要點在自身維度和在局部維度上的重要程度，下面結(jié)合全局因子給出綜合評價因子的定義。

圖5 極大值重要點的邊界因子示意圖Fig.5 Illustration for boundary factor of maximum important point

(22)

(22)式中：η為去噪閾值，表示去噪強度；λ為加權(quán)參數(shù)，表示邊界因子與全局因子相加的權(quán)值；Mean(FC)為所有重要點特征因子的平均值。

2.4 時間序列趨勢提取算法

1)獲取重要點序列。根據(jù)定義1提取時間序列S[n]的重要點序列SZ[m]作為分段點的備選集。

3)計算3個評價因子。根據(jù)定義4，定義5和定義6計算得到時間序列所有m個重要點的全局因子FG[m]，特征因子FC[m]和邊界因子FB[m]。

5)提取時間序列趨勢。將SP[d]中的所有分段點用直線段連接起來，得到時間序列的線性表示，完成對時間序列的特征趨勢提取。

本文采用遺傳算法結(jié)合十折交叉驗證算法來確定綜合評價模型的去噪閾值η和加權(quán)參數(shù)λ，選取使得擬合誤差E最小的一組參數(shù)值作為最終模型參數(shù)，簡要步驟如下。

1)進行數(shù)據(jù)預(yù)處理，將時間序列S[n]按照時間順序平均分為10等份，輪流將其中的9份作為訓(xùn)練集數(shù)據(jù)，1份作為測試集數(shù)據(jù)。

2)按照本文時間序列趨勢提取算法，分別計算得到各組訓(xùn)練集的全局因子FG、特征因子FC和邊界因子FB，并根據(jù)(22)式將去噪閾值η和加權(quán)參數(shù)λ設(shè)為決策變量，以獲取最優(yōu)PLR序列使得擬合誤差E最小化為最終目標(biāo)構(gòu)建適應(yīng)度函數(shù)。

3)初始化遺傳算法，設(shè)置種群大小為50，精英數(shù)目和交叉后代比例分別為30和0.75，最大進化代數(shù)、停止代數(shù)和適應(yīng)度函數(shù)偏差分別為60、60和1e-6，使得算法能夠在進化60代后停止。

本文趨勢提取算法的流程如圖6所示。

3 實驗與分析

3.1 實驗對比方法

本文將選擇上文介紹的4種時間序列的趨勢提取方法作為比較對象。

1)PAA方法[12]。根據(jù)壓縮率將時間序列分成等長的線性段，并獲取每段內(nèi)所有數(shù)據(jù)的均值。

2)PLR-BU方法[16]。先對時間序列上的所有數(shù)據(jù)點進行線性分割，依次對具有最小合并代價的相鄰線性段合并直至滿足要求。

3)PLR-TPTS方法[19]。先根據(jù)數(shù)據(jù)的角度變化值選擇轉(zhuǎn)折點，然后用轉(zhuǎn)折點表征趨勢段，通過尋找極值趨勢段進而選擇全局性的分段點。

4)PLR-DE方法[20]。先找到時間序列上所有重要點，并給出2個評價指標(biāo)對重要點的重要程度進行評價，最后根據(jù)評價結(jié)果來選取分段點。

圖6 本文趨勢提取算法流程圖Fig.6 Flow chart of the trend extraction algorithm

3.2 數(shù)值仿真實驗

仿真實驗序列為

(23)

(23)式中，t為時間序列整數(shù)，共600個數(shù)據(jù)。該序列包含了平穩(wěn)、不同趨勢的上升和下降等模式，可以比較全面準(zhǔn)確地測試算法性能。原始仿真數(shù)據(jù)示意圖如圖7所示。

圖7 原始仿真數(shù)據(jù)示意圖Fig.7 Illustration for raw simulation data

實驗中對序列S[600]加上均值μ=0，方差σ為0.5、1.0、1.5、2.0、2.5、3.0的噪聲，設(shè)定壓縮率C要求分別為75%、80%、85%、95%，使用遺傳尋優(yōu)算法得到去噪閾值η和加權(quán)參數(shù)λ，考察不同方法趨勢提取結(jié)果的擬合誤差E。顯然，分段表示序列越接近原始時間序列則其擬合誤差越小，對應(yīng)的算法性能也越好。不同壓縮率下的擬合誤差和如表1所示。由表1可知，隨著壓縮率的增加，各算法的擬合誤差逐漸增大，在考慮所有壓縮率的情況下，本文方法的擬合誤差總和分別比PAA、PLR-BU、PLR-TPTS和PLR-DE減小了66%、9%、64%和54%。本文方法在不同的壓縮率下都保持了最小的擬合誤差。

表1 不同壓縮率下的擬合誤差和Tab.1 Sum of fitting errors under different compression rates

圖8為各算法在不同壓縮率下添加不同噪聲的具體擬合誤差。可以看出，在一定的壓縮率情況下，各算法擬合誤差會隨著噪聲的增大而增大，本文方法雖受到噪聲的影響，但是擬合誤差在大多數(shù)噪聲下比其他方法小，去除噪聲的能力更強。

PAA方法擬合誤差雖然受到噪聲的影響很小，但是整體的擬合誤差較大，這是因為PAA方法的分段點是固定的，噪聲的增加并不影響分段點的選擇。PLR-TPTS方法由于其角度閾值是按照經(jīng)驗選取的，不同的角度閾值會造成擬合誤差的大幅變化。PLR-DE方法由于重要點評價因子沒有根據(jù)具體情況分類討論，在面對復(fù)雜的時間序列時擬合誤差會出現(xiàn)震蕩上升的情況。這使得上述3種方法的擬合誤差在壓縮率和噪聲的變化下出現(xiàn)了交叉。PLR-BU方法在低壓縮率和低噪聲時擬合誤差較小，但是在高壓縮率和高噪聲情況下會出現(xiàn)較大擬合誤差，這是因為PLR-BU方法只關(guān)注時間序列局部的信息而忽視了整體的趨勢，在高壓縮率高噪聲情況下會錯誤提取重要點，出現(xiàn)與本文方法的擬合誤差曲線出現(xiàn)交叉的情況。本文方法先考慮重要點在整體上的重要程度，并在考慮局部和單點情況時根據(jù)實際情況進行分類討論，讓評價因子更為精確，同時使用智能尋優(yōu)算法確定最優(yōu)參數(shù)，使得本文方法有效克服了上述方法存在的缺陷，在不同壓縮率的情況下都能保持較小擬合誤差，且在相同壓縮率下隨著噪聲的增大擬合誤差也比其他方法小。

圖8 不同壓縮率下的擬合誤差Fig.8 Fitting errors under different compression rates

3.3 實際應(yīng)用實驗

心電圖(ECG)是臨床上最常用的檢查之一，也是時間序列挖掘和模式識別研究中的重要領(lǐng)域。但是心電時間序列是一種高維度的數(shù)據(jù)，為了方便存儲、查詢和挖掘，需要對其進行壓縮表示。本文選用來自MIT-BIH正常竇性節(jié)律心電圖的數(shù)據(jù)庫，截取第1條通道的前10 s心電圖數(shù)據(jù)作為本文實際應(yīng)用驗證的時間序列，此序列包含16個心跳周期，共1 280個數(shù)據(jù)值。

實驗仍然采用上述4種方法作為本文方法的比較對象，考慮實際應(yīng)用中高壓縮率的需求，考察不同方法在壓縮率C分別為90%、91%、92%、93%、94%、95%時的擬合誤差E。

圖9為不同方法在不同壓縮率情況下對心電圖序列進行趨勢提取結(jié)果的擬合誤差?？梢钥闯?，本文方法在不同壓縮率情況下均保持了較小的擬合誤差。PAA方法的擬合誤差仍然對于壓縮率的變化不敏感，PLR-DE方法由于其對于重要點的評價不夠精準(zhǔn)，在高噪聲情況下隨著壓縮率增加，擬合誤差有大幅度的增加并超過了PLR-TPTS。PLR-BU方法在高噪聲情況下由于優(yōu)先考慮局部趨勢的缺陷，其擬合誤差在壓縮率增加下逐漸增加并在93%之后超過了PAA方法。擬合誤差的具體數(shù)值如表2所示。

表2 不同壓縮率下的擬合誤差具體數(shù)值Tab.2 Specific value of fitting errors under different compression rates

圖9 心電圖序列趨勢提取的擬合誤差Fig.9 Fitting errors of ECG trend extraction

圖10為在壓縮率為90%的情況下，不同方法的擬合曲線對比。由圖10可見，在對心電圖這種高噪聲的時間序列進行趨勢提取時，其他方法都出現(xiàn)了錯誤的提取，而本文方法能夠有效濾除干擾，準(zhǔn)確提取了所有最關(guān)鍵的重要點，使得擬合曲線準(zhǔn)確描述了原始時間序列的趨勢特征，保持了較小的擬合誤差。

圖10 90%的壓縮率下不同方法的擬合曲線對比Fig.10 Comparison of fitting curves of algorithms under 90% compression rate

4 結(jié) 論

本文在基于重要點的分段線性表示方法基礎(chǔ)上，提出一種基于模態(tài)分量重構(gòu)及多維評價的時間序列趨勢提取算法。首先，使用CEEMDAN對原始時間序列進行分解，并進行多次模態(tài)重構(gòu)；其次，基于模態(tài)重構(gòu)序列組給出重要點的全局因子，在整體維度上對時間序列進行刻畫；然后，使用特征因子結(jié)合去噪閾值濾除噪聲的干擾，使用邊界因子結(jié)合加權(quán)參數(shù)與全局因子進行運算得到綜合評價因子；最后，根據(jù)壓縮率的要求綜合評價重要點序列確定分段點實現(xiàn)時間序列的趨勢提取。仿真實驗表明，本文方法的擬合誤差分別比PAA、PLR-BU、PLR-TPTS和PLR-DE減小了66%、9%、64%和54%，其克服了現(xiàn)有趨勢提取方法的缺陷，有效實現(xiàn)了去噪，能準(zhǔn)確地提取時間序列的趨勢。將本文方法應(yīng)用于心電圖序列趨勢提取，其擬合誤差分別比PAA、PLR-BU、PLR-TPTS和PLR-DE減小了61%、55%、71%和73%，這進一步表明，在高噪聲情況下，本文所提方法具有良好的去噪能力，能夠全面準(zhǔn)確地提取時間序列的趨勢特征。