基于EEMD的固定分段數(shù)分段線性表示方法

劉學(xué)彬,梁智飛,朱衛(wèi)平,祝 凱*

(1.青島理工大學(xué) 信息與控制工程學(xué)院,山東 青島 266000;2.中石油煤層氣有限責(zé)任公司,北京 102200)

0 引 言

由于時(shí)間序列是高維且存在大量噪音的,直接在原始序列上進(jìn)行預(yù)測、模式發(fā)現(xiàn)和分類等挖掘任務(wù)的效率較低,同時(shí)也會(huì)影響挖掘結(jié)果的精度和可信度。因此,使用特征表示方法將時(shí)間序列從高維度轉(zhuǎn)換到低維度,這種方法可以在降低時(shí)間序列復(fù)雜度的同時(shí),保留時(shí)間序列的主要信息,為進(jìn)一步深入研究時(shí)間序列奠定基礎(chǔ)[1]。

目前國內(nèi)外有不少學(xué)者致力于時(shí)間序列特征表示方法的研究,時(shí)間序列特征表示方法的主要代表有:基于域變換的表示方法(離散傅里葉變換[2]和離散小波變換[3];符號(hào)化表示方法,其中應(yīng)用最廣泛的是符號(hào)聚合近似方法[4-5];分段累計(jì)近似方法[6]和分段線性表示(Piecewise Linear Representation,PLR)[7]。其中PLR具有簡單、直觀的特點(diǎn),能夠有效保留原序列的形態(tài)信息以減少擬合誤差,是一種應(yīng)用廣泛的時(shí)間序列特征表示方法。因此,該文著眼于分段線性表示方法的研究和改進(jìn)。

目前,PLR的研究主要集中于解決分段數(shù)和分段點(diǎn)的選擇問題上。為了解決這些問題,時(shí)序的分段表示方法可以分為以下幾種:(1)限制分段數(shù):主要代表是分段累計(jì)近似方法,但該方法沒有考慮實(shí)際序列形態(tài),不能很好地保留原始序列特征;(2)限制分段誤差:主要代表性算法有自頂向下[8]、自底向上[9]、滑動(dòng)窗口[10]。限制分段誤差方法對一些狀態(tài)變化的拐點(diǎn)不敏感,不能保證每一分段只具有一種基本趨勢。針對上述問題,近年來不少學(xué)者提出了一些改進(jìn)方法。例如,尚福華[11]和廖俊[12]提出基于趨勢轉(zhuǎn)折點(diǎn)的分段線性表示方法;陳帥飛[13]提出基于關(guān)鍵點(diǎn)的分段線性表示方法;劉意楊[14]提出基于轉(zhuǎn)折點(diǎn)和趨勢段的分段線性表示方法等。但是,這些方法使用單一的啟發(fā)式規(guī)則,難以適用于數(shù)據(jù)分布復(fù)雜的時(shí)間序列,進(jìn)而導(dǎo)致算法出現(xiàn)局部最優(yōu)化問題,而且不能靈活控制壓縮率,不能適應(yīng)后期要求分段數(shù)一定的應(yīng)用[15]。

針對上述方法存在局部最優(yōu)化和不能預(yù)計(jì)分段數(shù)的問題,提出了基于EEMD的固定分段數(shù)分段線性表示方法。首先,通過模態(tài)重構(gòu)思想過濾掉細(xì)節(jié)信息,提取到全局性分段點(diǎn);然后,根據(jù)各初始分段子序列的波動(dòng)程度,確定子序列段內(nèi)分段點(diǎn)數(shù)量分布;最后,采用基于分段數(shù)閾值的自底向上方法將子序列合并到要求的分段數(shù)。

1 分段線性表示相關(guān)概念及問題描述

1.1 分段線性表示相關(guān)概念

(1)

(2)

{xpq-1≤xpq}∩{xpq+1

或

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

(3)

此外,規(guī)定一個(gè)有限長度的時(shí)間序列起點(diǎn)和終點(diǎn)為重要點(diǎn)。由式(3)得到m個(gè)重要點(diǎn),則重要點(diǎn)序列表示為:

(4)

1.2 問題描述

傳統(tǒng)的算法采用單一的啟發(fā)式規(guī)則提取局部特征點(diǎn),當(dāng)原始時(shí)間序列波動(dòng)頻率較為劇烈且集中時(shí),容易出現(xiàn)多個(gè)點(diǎn)的斜率變化近似。時(shí)間序列如圖1所示。

圖1 斜率波動(dòng)頻繁劇烈的情況

圖1中序列點(diǎn)a,b,c,d,e,f點(diǎn)斜率變化近似,當(dāng)通過調(diào)節(jié)斜率變化閾值d使得達(dá)到要求的壓縮率時(shí),會(huì)出現(xiàn)臨界閾值,如下:

其中,下標(biāo)L表示左,R表示右。由上述公式和圖1知,c,d兩點(diǎn)作為反映序列整體趨勢的特征點(diǎn)因斜率變化小而“漏提取”,即分段方法的結(jié)果遺漏掉能夠反映整體特征的數(shù)據(jù)點(diǎn);由此可認(rèn)為,b,e兩點(diǎn)為“過提取”,即分段方法的結(jié)果提取到不能反映整體特征的數(shù)據(jù)點(diǎn),導(dǎo)致算法陷入局部最優(yōu)化。

2 集合經(jīng)驗(yàn)?zāi)B(tài)分解和改進(jìn)的自底向上分段

2.1 集合經(jīng)驗(yàn)?zāi)B(tài)分解

Huang[17]提出了經(jīng)驗(yàn)?zāi)B(tài)分解(Empirical Mode Decomposition,EMD)。該方法的核心思想是將復(fù)雜的信號(hào)分解為有限個(gè)頻率從高到低的本征模態(tài)函數(shù)(Intrinsic Mode Functions,IMF),對于某時(shí)間序列{x(t)}經(jīng)驗(yàn)?zāi)B(tài)分解的具體步驟如下:

(1)求出{x(t)}中所有的極值。

(2)采用3次樣條函數(shù)進(jìn)行插值擬合上包絡(luò)線bmax(t)和下包絡(luò)線bmin(t)。

(3)計(jì)算上下包絡(luò)線平均值m(t):

m(t)=[bmax(t)+bmin(t)]

(5)

(4)從時(shí)間序列中提取均值并將x(t)和m(t)的差定義為:

d(t)=x(t)-m(t)

(6)

(5)檢查d(t)的屬性:如果滿足IMF分量條件,則將d(t)表示為第k個(gè)IMF,并將x(t)替換為殘差r1(t)=x(t)-d(t)。第k個(gè)IMF分量通常表示為ck(t);如果不滿足,則將x(t)替換為d(t)。

(6)重復(fù)步驟(1)～(5)直到殘差為單調(diào)函數(shù)為止。

原始時(shí)間序列可以表示為若干個(gè)IMF和一個(gè)殘差的線性組合:

(7)

其中,x(t)表示1維信號(hào);ck(t)表示第k個(gè)IMF分量;r(t)表示殘余。

當(dāng)時(shí)間序列的時(shí)間尺度呈現(xiàn)跳躍性時(shí),采用EMD對其進(jìn)行分解,將會(huì)產(chǎn)生一個(gè)IMF分量包含不同時(shí)間尺度特征成分的情況,這種現(xiàn)象被稱為模態(tài)混疊[18],它使得EMD得到的分解結(jié)果的可靠性和可解性受到影響。Wu[18]提出了集合經(jīng)驗(yàn)?zāi)B(tài)分解(Ensemble Empirical Mode Decomposition,EEMD)解決這一問題?；舅枷胧菍⒉煌自肼暥啻渭尤朐紩r(shí)間序列以消除模態(tài)混疊現(xiàn)象。

如圖2(a)所示,對1組示例時(shí)間序列進(jìn)行EEMD分解,得到了6個(gè)IMF分量和1個(gè)RES殘余,如圖2(b)所示。

2.2 IMF重構(gòu)

Zhang等人[19]采用EEMD技術(shù)來分析石油價(jià)格

(a)示例時(shí)間序列

(b)EEMD分解結(jié)果圖2 示例時(shí)間序列集合經(jīng)驗(yàn)?zāi)B(tài)分解

變化。他們發(fā)現(xiàn),經(jīng)本征模態(tài)函數(shù)重構(gòu)后的序列可以很好地反映序列的關(guān)鍵轉(zhuǎn)折點(diǎn)和整體趨勢變化?；谶@項(xiàng)研究,該文使用EEMD技術(shù)對時(shí)間序列進(jìn)行分解,并將分解得到的IMF分為高頻部分、低頻部分和殘余。前兩個(gè)成分能夠揭示時(shí)間序列所蘊(yùn)含的物理意義,并發(fā)現(xiàn)時(shí)序的一些新特征。對EEMD分解得到的N個(gè)IMF,求出每個(gè)IMF的平均值,得到用于分解高頻和低頻分量的K函數(shù)。以圖2(a)的時(shí)間序列為例,構(gòu)建的K函數(shù)及高、低頻和殘余分量如圖3所示。

(a)分解高頻和低頻分量的K函數(shù)

(b)原始時(shí)序和3個(gè)分量圖3 K函數(shù)及對應(yīng)的3個(gè)分量

由圖3(a)知,在IMF5處,平均值開始偏離零點(diǎn),因此使用IMF1～I(xiàn)MF4的部分重構(gòu)表示高頻分量,使用IMF5和IMF6的部分重構(gòu)表示低頻分量,殘余單獨(dú)處理。圖3(b)顯示了原始時(shí)間序列和3個(gè)分量。殘余反映時(shí)間序列長期緩慢變化;低頻分量的每次急劇上升或下降可能對應(yīng)1個(gè)物理事件或是某種程度上的噪聲表征;而高頻分量通過去除大量的小幅波動(dòng)使得可以反映時(shí)間序列的整體變化趨勢。下面給出模態(tài)重構(gòu)序列的定義。

定義4(模態(tài)重構(gòu)序列):對于某時(shí)間序列X,對X進(jìn)行EEMD分解得到N個(gè)IMF,定義參與重構(gòu)的起始IMF索引為s,終止索引為e,重構(gòu)序列XR表示為:

(8)

在高頻分量基礎(chǔ)上,提取全局特征點(diǎn),實(shí)現(xiàn)時(shí)間序列的初始分段。

(9)

根據(jù)上式,對圖3(b)中的高頻分量提取全局特征點(diǎn),實(shí)現(xiàn)時(shí)間序列的初始分段。

由圖4知,原始時(shí)間序列被全局特征點(diǎn)分割為12段子序列,每段子序列都保持整體上升、下降、保持三種基本趨勢,有效去除大量小幅波動(dòng),反映時(shí)間序列整體變化趨勢。

圖4 全局特征點(diǎn)初始分段

2.3 時(shí)間序列符號(hào)化

假設(shè)在序列中需要查找N個(gè)分段點(diǎn),上節(jié)已提取了M個(gè)全局特征點(diǎn),并將原時(shí)間序列分成了M+1個(gè)初始段。接下來,采用廖俊[12]提到的時(shí)間序列點(diǎn)間的模式變化提取剩下的N-M個(gè)分段點(diǎn),如圖5所示。

為了反映時(shí)間序列內(nèi)的模式變化,將所有時(shí)間序列數(shù)據(jù)點(diǎn)符號(hào)化[20]。在時(shí)間序列X中,給定某一序列點(diǎn)xj,然后分別用前一點(diǎn)xi和后一點(diǎn)xk與該點(diǎn)做差分,即xk-xj=Q和xj-xi=P。具體步驟如下[1]:

(3)將不符合上述條件的點(diǎn)用“0”表示。

(4)遍歷整個(gè)序列,得到符號(hào)化序列。

圖5 時(shí)間序列3點(diǎn)之間的模式變化

其中δ為自定義閾值,將所有符號(hào)化的子序列分別求和,存入Hi中,得到長度為M+1的序列:H={H1,H2,…,HM+1},通過以下公式:

(10)

得到M+1個(gè)子序列內(nèi)分段點(diǎn)的分布數(shù)量:

C={C1,C2,…,CM+1}

(11)

2.4 固定分段數(shù)的自底向上分段

經(jīng)典的自底向上方法由Keogh等人[7]提出,該方法的基本思想是通過循環(huán)地合并誤差最小的相鄰分段,直到所有的擬合誤差均不小于分段閾值為止。該算法存在偶數(shù)限制的不足,為了解決該問題,孫煥良在其[21]研究中提出了優(yōu)化的PLR_BU算法,但是仍無法準(zhǔn)確地預(yù)測時(shí)間序列的分段數(shù)。

針對這一缺陷,該文在優(yōu)化的PLR_BU算法基礎(chǔ)上進(jìn)行了改進(jìn),提出了固定的PLR_BU算法。該算法的基本思想:首先將長度為n的時(shí)間序列依次相連前后兩點(diǎn),然后給定分段數(shù)閾值,循環(huán)地執(zhí)行下述過程:(1)計(jì)算相鄰的分段合并后的擬合誤差;(2)查找擬合誤差最小的相鄰分段{xi,xj,xk},移除此相鄰分段的中心點(diǎn)xj,序列長度減1;(3)計(jì)算新生成的段與前后分段的擬合誤差。重復(fù)上述過程,直到合并到滿足設(shè)定的分段數(shù)為止。固定分段數(shù)的PLR_BU算法偽代碼如表1所示。

2.5 時(shí)間序列分段線性表示方法

該文提出基于EEMD的固定分段數(shù)分段線性表示方法,具體算法步驟如下:

給定時(shí)間序列X={x1,x2,…,xn},斜率變化閾值δ,分段數(shù)N。

(3)時(shí)間序列符號(hào)化和確定初始分段段內(nèi)分段點(diǎn)分布。根據(jù)斜率變化將時(shí)間序列轉(zhuǎn)換成由“0”和“1”組成的符號(hào)化序列。計(jì)算符號(hào)化后各初始分段內(nèi)數(shù)據(jù)和,得到H={H1,H2,…,HM+1},根據(jù)式(10)和(11),得到最終子序列段內(nèi)分段點(diǎn)分布數(shù)量序列C={C1,C2,…,CM+1}。

(4)固定分段數(shù)的PLR_BU算法確定最終分段點(diǎn)。根據(jù)改進(jìn)的PLR_BU算法對子序列繼續(xù)分段,直到分段數(shù)為N,最終分段點(diǎn)序列:X={x1,xi1,…,xiN-2,xn}。

表1 固定分段數(shù)的PLR_BU算法

參數(shù)說明:在文中方法中,斜率變化閾值δ是主要的參數(shù)。設(shè)置閾值δ的目的是按照斜率變化過濾數(shù)據(jù)點(diǎn),δ值過小時(shí),會(huì)將斜率變化相對較小的數(shù)據(jù)點(diǎn)也轉(zhuǎn)換為“1”,導(dǎo)致相對平緩的序列段分段點(diǎn)的分布數(shù)量也較多;δ值過大時(shí),會(huì)將斜率變化相對較大的數(shù)據(jù)點(diǎn)轉(zhuǎn)換為“1”,導(dǎo)致只有斜率變化相對較大的序列段分段點(diǎn)的分布數(shù)量才會(huì)較多。由于P和Q表示相鄰點(diǎn)的差分,所以δ取值范圍為0<δ

3 實(shí)驗(yàn)與分析

3.1 實(shí)驗(yàn)對比方法

該文選擇以下3種時(shí)間序列分段線性表示方法作為比較對象。

(1)PAA[6]。

(2)PLR_ITTP[12]。

(3)PLR_TRIP[14]。

3.2 仿真數(shù)據(jù)驗(yàn)證

X=

(12)

式中,t為整數(shù),共600個(gè)數(shù)據(jù)。選擇該仿真序列作實(shí)驗(yàn),是由于這個(gè)序列的重要點(diǎn)比較清晰,且沒有噪聲的干擾,因此,重要點(diǎn)更容易被發(fā)現(xiàn)。對序列X加上均值(μ=0)、方差(σ設(shè)為0.5～3)、步長為0.5的隨機(jī)誤差,對比在不同噪聲的情況下,各個(gè)分段方法的抗噪聲能力,檢驗(yàn)是否可以準(zhǔn)確提取全局特征點(diǎn)的有效性,對于對比分析時(shí)間序列消除局部最優(yōu)化問題具有參考意義。

根據(jù)仿真數(shù)據(jù)的實(shí)際趨勢,將序列分為13段,將擬合誤差作為評(píng)價(jià)標(biāo)準(zhǔn),顯然,提取的分段點(diǎn)越接近原始序列趨勢分段點(diǎn),擬合誤差越小。實(shí)驗(yàn)結(jié)果如圖6和表2所示。

表2 不同噪聲下不同PLR的擬合誤差

圖6是不同噪聲情況下,不同分段方法的分段擬合結(jié)果。

圖6 不同噪聲情況下的不同PLR結(jié)果

由圖6可知,隨著σ的逐漸增大,分段算法對能夠反映整體趨勢的分段點(diǎn)的識(shí)別越來越困難,而文中方法相比于其他3種方法,對該序列中能夠反映整體趨勢的分段點(diǎn)的識(shí)別較為準(zhǔn)確,尤其在高噪聲情況下更為明顯,而其他3種方法均提取了錯(cuò)誤的分段點(diǎn)。

由表2知,文中方法雖受噪聲干擾,但總的來說,抗噪聲干擾的能力比其他3種方法有所加強(qiáng),可以非常準(zhǔn)確地提取反映整體趨勢的分段點(diǎn),而其他方法則極易受到噪聲的干擾,導(dǎo)致陷入局部最優(yōu)狀態(tài)。

由圖6和表2可知,盡管噪聲的增加對PAA的擬合誤差影響較小,但其整體的擬合誤差相對較大,這是由于PAA采用等長分段,分段點(diǎn)的選取不會(huì)受到噪聲的影響。而PLR_ITTP、PLR_TRIP都通過某種抗噪機(jī)制削弱了噪聲的干擾,使得擬合誤差相對較小。PLR_ITTP在噪聲較低時(shí),擬合誤差較小,但在噪聲較高時(shí),擬合誤差隨之變大,這是因?yàn)镻LR_ITTP只重視時(shí)間序列的局部特征,而忽略了全局意義下的時(shí)間轉(zhuǎn)折點(diǎn),這將導(dǎo)致在高噪聲的時(shí)候,算法會(huì)錯(cuò)誤地提取關(guān)鍵點(diǎn)。PLR_TRIP在不同噪聲下擬合誤差都較大,并且出現(xiàn)了震蕩上升的情況,由于PLR_TRIP的角度閾值和趨勢段閾值的組合選取是復(fù)雜的,不同的閾值組合會(huì)造成擬合誤差的大幅變化,同時(shí)PLR_TRIP提出趨勢段的概念來削弱噪聲的干擾也是只關(guān)注局部信息。而文中方法先通過模態(tài)重構(gòu)方法得到全局分段點(diǎn),使得文中方法有效克服上述方法存在的局部最優(yōu)化缺陷,之后使用自底向上方法進(jìn)行融合,保留了基于分段誤差的分段算法擬合誤差較小的優(yōu)點(diǎn)。

3.3 工業(yè)實(shí)例應(yīng)用

壓裂施工過程中,通過記錄不同時(shí)間段的施工壓力、泵注排量和加砂體積分?jǐn)?shù)獲得壓裂施工曲線,有效利用壓裂施工曲線并提取有效信息,不僅能夠?qū)?chǔ)層及裂縫參數(shù)再認(rèn)識(shí),而且對于指導(dǎo)壓裂施工以及調(diào)整壓裂設(shè)計(jì)方案、提高壓裂技術(shù)水平和施工效果有重要的借鑒作用。因此,對壓裂施工曲線的挖掘和分析有著重要的工程意義和應(yīng)用價(jià)值[22]。但壓裂施工曲線是一種高維數(shù)據(jù),為了方便后續(xù)存儲(chǔ)和挖掘需要對其進(jìn)行壓縮表示。該文選用某區(qū)塊的壓裂施工數(shù)據(jù),時(shí)間間隔為1 s,共2 500個(gè)數(shù)據(jù)值。

為了比較文中算法與其他算法的優(yōu)劣,并考慮到實(shí)際壓裂施工曲線的高維性特點(diǎn),將考察所有方法在壓縮率分別為90%、92%、94%、96%、98%時(shí)的擬合誤差,實(shí)驗(yàn)結(jié)果如圖7所示。

圖7 壓裂施工曲線不同壓縮率下的擬合誤差

由圖7可知,文中方法在不同壓縮率下的擬合誤差都是最小的,優(yōu)于其他3種方法。

圖8 96%的壓縮率下不同方法分段擬合效果對比

由圖8可知,對壓裂施工時(shí)間序列這種高噪聲且分布復(fù)雜的序列進(jìn)行趨勢提取時(shí),PLR_ITTP、PLR_TRIP因存在局部最優(yōu)化的問題而導(dǎo)致“漏提取”“過提取”現(xiàn)象較為嚴(yán)重,丟失了反映整體趨勢的重要點(diǎn);而文中算法能夠有效地去除時(shí)間序列中的噪聲,準(zhǔn)確提取反映時(shí)間序列整體趨勢的分段點(diǎn)。

4 結(jié)束語

針對現(xiàn)有方法的不足,該文提出一種基于EEMD的固定分段數(shù)表示方法。仿真實(shí)驗(yàn)結(jié)果表明:該方法的擬合誤差分別比PAA、PLR_ITTP、PLR_TRIP平均減小了80%、59%、78%,其有效地解決了現(xiàn)有分段方法存在的問題,極大地削弱了噪聲干擾,從而能夠準(zhǔn)確地找到反映整體趨勢的分段點(diǎn)。最后將該方法應(yīng)用于壓裂施工數(shù)據(jù),其擬合誤差分別比PAA、PLR_ITTP、PLR_TRIP減小了86%、84%、89%,再次證明了所提方法對趨勢提取的有效性和準(zhǔn)確性。