多元時(shí)間序列缺失數(shù)據(jù)填補(bǔ)方法

李正欣, 張鳳鳴, 王 瑛, 陶 茜, 李 超

(空軍工程大學(xué)裝備管理與無(wú)人機(jī)工程學(xué)院, 陜西 西安 710051)

0 引 言

按時(shí)間順序獲取的一系列觀(guān)測(cè)值xt(j)稱(chēng)為時(shí)間序列,其中,t(t=1,2,…,n)表示第t個(gè)時(shí)刻,j(j=1,2,…,m)表示第j個(gè)變量,xt(j)表示第j個(gè)變量在第t個(gè)時(shí)刻上的記錄值[1]。當(dāng)m=1時(shí),xt(j)為一元時(shí)間序列;當(dāng)m>1時(shí),xt(j)為多元時(shí)間序列。

時(shí)間序列廣泛存在于氣象、水文、醫(yī)學(xué)和經(jīng)濟(jì)等眾多領(lǐng)域[2]。與一元時(shí)間序列相比,多元時(shí)間序列更具一般性,例如,股票交易記錄包含交易量、開(kāi)盤(pán)價(jià)、收盤(pán)價(jià)、最高價(jià)和最低價(jià)[3];“黑匣子”記錄的飛行數(shù)據(jù)包含速度、高度、俯仰角、傾斜角和航向角等。廣義上,任何包含多變量的數(shù)據(jù)存儲(chǔ)都可以視為多元時(shí)間序列[4]。

時(shí)間序列中隱含著深層知識(shí),針對(duì)時(shí)間序列的數(shù)據(jù)挖掘已成為目前數(shù)據(jù)挖掘領(lǐng)域中的研究熱點(diǎn)[5]。然而,在實(shí)際應(yīng)用中,數(shù)據(jù)采集受多種因素的干擾,難免存在缺失數(shù)據(jù),造成序列不連續(xù),從而影響后續(xù)的分析處理。在數(shù)據(jù)挖掘領(lǐng)域UCI公開(kāi)數(shù)據(jù)集中,超過(guò)40%的數(shù)據(jù)庫(kù)都含有缺失數(shù)據(jù)[6]。因此,研究多元時(shí)間序列的缺失數(shù)據(jù)填補(bǔ)方法是分析處理時(shí)間序列的前提、具有重要的應(yīng)用價(jià)值。

1 相關(guān)研究

傳統(tǒng)的數(shù)據(jù)填補(bǔ)方法有均值填補(bǔ)、回歸填補(bǔ)、期望值最大化方法等[7],這些方法計(jì)算簡(jiǎn)單,但當(dāng)缺失數(shù)據(jù)量較多、波動(dòng)較大時(shí),填補(bǔ)精度大幅度下降。

隨后在傳統(tǒng)方法基礎(chǔ)上,又融入了統(tǒng)計(jì)分析、機(jī)器學(xué)習(xí)等方法,在一定程度上提高了填補(bǔ)精度。然而,這些成果主要針對(duì)傳統(tǒng)數(shù)據(jù)類(lèi)型,多元時(shí)間序列在時(shí)間維和變量維上,都具有較強(qiáng)的相關(guān)性,與傳統(tǒng)數(shù)據(jù)具有較大的差異,因此,這些方法并不完全適用于多元時(shí)間序列。

目前,針對(duì)多元時(shí)間序列填補(bǔ)方法的研究成果相對(duì)較少[7],主要有多元回歸、神經(jīng)網(wǎng)絡(luò)、支持向量機(jī)(support vector machine,SVM)和組合填補(bǔ)方法等。

多元自回歸模型[8]是一種經(jīng)典的時(shí)間序列分析方法,它適用于線(xiàn)性模型,但對(duì)非線(xiàn)性的時(shí)間序列效果不佳。

神經(jīng)網(wǎng)絡(luò)[9]可以對(duì)任意非線(xiàn)性函數(shù)進(jìn)行逼近,具有非線(xiàn)性、自學(xué)習(xí)的特點(diǎn),可用于時(shí)間序列缺失數(shù)據(jù)填補(bǔ)。但基于經(jīng)驗(yàn)風(fēng)險(xiǎn)最小化原則,要求較大的樣本規(guī)模,對(duì)于小樣本情況,容易出現(xiàn)過(guò)擬合、泛化能力不足等問(wèn)題。

SVM[10]追求結(jié)構(gòu)風(fēng)險(xiǎn)最小化,專(zhuān)門(mén)針對(duì)有限樣本,得到現(xiàn)有信息下的全局最優(yōu)解,且計(jì)算復(fù)雜度與維度無(wú)關(guān)。然而,大多數(shù)基于SVM的填補(bǔ)方法[11-12],僅關(guān)注缺失數(shù)據(jù)所屬的單變量信息,沒(méi)有針對(duì)多元時(shí)間序列的特點(diǎn)、有效利用其他相關(guān)的多變量信息。文獻(xiàn)[13]提出了一種基于狀態(tài)匹配與SVM的多變量填補(bǔ)方法。首先,判斷存在缺失數(shù)據(jù)的序列所屬的類(lèi);然后,分析該類(lèi)序列中,同缺失數(shù)據(jù)所屬變量相關(guān)的其他變量;最后,在相關(guān)變量中選擇樣本點(diǎn),使用SVM填補(bǔ)缺失數(shù)據(jù)。然而,該方法在判斷序列所屬類(lèi)別及分析相關(guān)變量時(shí),使用人工判讀,對(duì)經(jīng)驗(yàn)依賴(lài)性較強(qiáng)。

文獻(xiàn)[14]提出基于最小二乘支持向量機(jī)(least squares support vector machine,LSSVM)的組合權(quán)重填補(bǔ)方法,同時(shí)利用缺失數(shù)據(jù)所屬的單變量以及與其相關(guān)的多變量信息,對(duì)單變量填補(bǔ)值與多變量填補(bǔ)值加權(quán)求和,作為最終填補(bǔ)結(jié)果。然而,該方法用一個(gè)權(quán)重值處理所有的缺失數(shù)據(jù),難以體現(xiàn)各個(gè)缺失數(shù)據(jù)的差異性;此外,權(quán)重值在真實(shí)數(shù)據(jù)已知的情況下才能得到最優(yōu)結(jié)果,限制了該方法的使用。

時(shí)間序列通常具有非平穩(wěn)、非線(xiàn)性、非結(jié)構(gòu)化等特征,此外,具有典型模式的序列通常樣本規(guī)模較小。針對(duì)以上特征,提出一種基于LSSVM的組合閾值填補(bǔ)方法。首先,以存在缺失數(shù)據(jù)的多元時(shí)間序列為對(duì)象,搜索與其同類(lèi)的相似序列;然后,利用LSSVM分別進(jìn)行多變量填補(bǔ)和單變量填補(bǔ);最后,根據(jù)多變量和單變量填補(bǔ)結(jié)果的差異度,進(jìn)行組合閾值填補(bǔ)。

2 基于LSSVM的組合閾值填補(bǔ)

LSSVM把最小二乘線(xiàn)性系統(tǒng)引入到SVM中,將解二次規(guī)劃問(wèn)題轉(zhuǎn)化為求解線(xiàn)性方程組,能有效提高求解速度與泛化能力[15]。

2.1 LSSVM

訓(xùn)練樣本集D={(xk,yk)}中:xk∈Rn表示輸入、yk∈R表示輸出,原始空間中的函數(shù)估計(jì)轉(zhuǎn)化為求解[14]

s.t.yk=ωTφ(xk)+b+ek,k=1,2,…,N

(1)

式中,γ為正則化參數(shù)；φ(·):Rn→Rnb是核空間映射函數(shù)；權(quán)矢量ω∈Rnb；誤差變量ek∈R；b是偏差量。

可以構(gòu)造拉格朗日方程求解最優(yōu)化問(wèn)題為

L(ω,b,e,α)=J(ω,e)-

(2)

拉格朗日乘子αk∈R。令L對(duì)ω,b,ek,αk的偏導(dǎo)數(shù)分別等于0,消去ek和ω可得

(3)

式中,Ω=K(xk,xl)=φ(xk)Tφ(xl)；y=[y1,y2,…,yN]T;α=[α1,α2,…,αN]T;1v=[1,1,…,1]T。

根據(jù)Mercer條件,存在映射函數(shù)φ和核函數(shù)K(·,·)使得

K(xk,xl)=φ(xk)Tφ(xl)

(4)

則LSSVM的函數(shù)估計(jì)為

(5)

α、b由式(3)求解得到。

2.2 組合閾值填補(bǔ)方法

利用LSSVM填補(bǔ)缺失數(shù)據(jù)有兩種基本思路:單變量填補(bǔ)和多變量填補(bǔ)。單變量填補(bǔ)只利用缺失數(shù)據(jù)所屬的單變量信息,進(jìn)行支持向量回歸,填補(bǔ)缺失數(shù)據(jù)。多變量填補(bǔ)利用與缺失數(shù)據(jù)所屬變量相關(guān)的其他多個(gè)變量信息,對(duì)缺失數(shù)據(jù)進(jìn)行填補(bǔ)。

單變量填補(bǔ)依據(jù)變量自身的內(nèi)在規(guī)律,多變量填補(bǔ)利用變量間的相互關(guān)系。組合閾值填補(bǔ)方法同時(shí)利用單變量和其他多個(gè)相關(guān)變量信息。

不失一般性,設(shè)數(shù)據(jù)集中,序列X1的第i個(gè)變量維度上存在缺失數(shù)據(jù)。組合閾值填補(bǔ)方法的基本流程如圖1所示。

圖1 組合閾值填補(bǔ)方法的基本流程Fig.1 Process of the threshold filling method

(1) 構(gòu)建訓(xùn)練集。在多元時(shí)間序列數(shù)據(jù)集中,搜索與X1相似的其他序列,結(jié)果記為X2,X3,…,Xn,構(gòu)成訓(xùn)練集。

(2) 多變量填補(bǔ)。在X2,X3,…,Xn的每個(gè)時(shí)刻上,以其他(m-1)個(gè)變量維度上的數(shù)據(jù)作為輸入、第i個(gè)變量維度上的數(shù)據(jù)作為輸出,訓(xùn)練LSSVM。然后,針對(duì)X1的缺失數(shù)據(jù)段,利用訓(xùn)練后的LSSVM模型進(jìn)行填補(bǔ),得到多變量填補(bǔ)結(jié)果。

(3) 單變量填補(bǔ)。在訓(xùn)練集的每個(gè)序列中,針對(duì)第i個(gè)變量維度,以k個(gè)變量值作為輸入,第k+1個(gè)變量值作為輸出,訓(xùn)練LSSVM。然后,針對(duì)X1的缺失數(shù)據(jù)段,利用訓(xùn)練后的LSSVM模型進(jìn)行填補(bǔ),得到單變量填補(bǔ)結(jié)果。

(4) 組合閾值填補(bǔ)。以多變量和單變量的填補(bǔ)結(jié)果為基礎(chǔ),按照

(6)

2.3 核函數(shù)與參數(shù)優(yōu)化

LSSVM模型中,核函數(shù)把樣本投影到高維空間,把非線(xiàn)性問(wèn)題轉(zhuǎn)化為線(xiàn)性問(wèn)題。由于高斯徑向基核函數(shù)(radial basis function,RBF)具有良好性能,且只有一個(gè)形狀參數(shù),因此,文中LSSVM采用RBF核函數(shù)。

LSSVM模型涉及2個(gè)參數(shù):正則化參數(shù)γ和RBF核函數(shù)的形狀參數(shù)σ。γ表示擬合誤差權(quán)重,γ越大,擬合函數(shù)越逼近樣本點(diǎn)。σ影響擬合結(jié)果平滑程度,σ越大,擬合結(jié)果越平滑,但會(huì)出現(xiàn)擬合效果不佳的情況;σ過(guò)小,擬合結(jié)果會(huì)出現(xiàn)起伏。為了達(dá)到較好的填補(bǔ)效果,應(yīng)合理選擇參數(shù),保證擬合曲線(xiàn)比較平滑,同時(shí)具有較小的擬合誤差。

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

實(shí)驗(yàn)環(huán)境為:Intel(R) Core(TM) i7-3770 CPU、4G RAM、Windows 7和Matlab R2010a。

用已知分類(lèi)結(jié)果的公開(kāi)數(shù)據(jù)集ASL(Australian sign language)[16]作為實(shí)驗(yàn)數(shù)據(jù)。ASL是一組手語(yǔ)信號(hào),包含22個(gè)變量,左、右手動(dòng)作各用11個(gè)變量描述:5個(gè)變量表示各手指的彎曲程度;6個(gè)變量(6個(gè)自由度)表示手的位置。數(shù)據(jù)集包含95種語(yǔ)意(95個(gè)類(lèi)),每種語(yǔ)意包含27個(gè)序列,一共2 565個(gè)序列。

不失一般性,分別用前兩種語(yǔ)意alive、answer的第27個(gè)序列(記為alive_27、answer_27)作為實(shí)驗(yàn)對(duì)象,序列長(zhǎng)度分別為57、54。刪除第6個(gè)變量維度上的最后10個(gè)數(shù)據(jù),模擬缺失數(shù)據(jù),如圖2所示。比較多變量填補(bǔ)、單變量填補(bǔ)、組合權(quán)重填補(bǔ)[14]和組合閾值填補(bǔ)等4種方法的填補(bǔ)效果,以平均相對(duì)誤差e作為比較依據(jù)。

圖2 存在缺失數(shù)據(jù)的alive_27Fig.2 Sequence containing missing data (alive_27)

(7)

式中,n為缺失數(shù)據(jù)個(gè)數(shù);si,ri分別表示第i個(gè)填補(bǔ)值和真實(shí)值。

3.1 填補(bǔ)效果的比較

為了客觀(guān)地比較填補(bǔ)效果,4種方法均以相似性搜索[17]的結(jié)果作為訓(xùn)練集,本文重點(diǎn)是數(shù)據(jù)填補(bǔ),因此直接用與alive_27、answer_27同類(lèi)的其他26個(gè)序列作為相似性搜索結(jié)果、構(gòu)成訓(xùn)練集。

3.1.1 針對(duì)alive_27的填補(bǔ)

先對(duì)alive_27進(jìn)行填補(bǔ),多變量填補(bǔ)γ=2,σ=5時(shí),e達(dá)到最小值3.6×10-4。針對(duì)不同的步長(zhǎng)k,進(jìn)行單變量填補(bǔ),遍歷γ,σ各種組合,當(dāng)k=10時(shí),e達(dá)到最小值3.3×10-4,對(duì)應(yīng)的γ=38,σ=67。

組合權(quán)重、組合閾值填補(bǔ)均以多變量、單變量填補(bǔ)結(jié)果為基礎(chǔ),涉及參數(shù)分別為權(quán)重ω、閾值t。當(dāng)ω=0.6時(shí),組合權(quán)重填補(bǔ)e=3.16×10-4達(dá)到最小值。用式(6)計(jì)算每個(gè)時(shí)刻,多變量與單變量填補(bǔ)的差異度,如圖3所示。圖3(a)為針對(duì)alive_27的差異度,在10個(gè)填補(bǔ)數(shù)據(jù)中,時(shí)刻50、57上的差異度較大,閾值t=4×10-4時(shí),組合閾值填補(bǔ)用多變量與單變量填補(bǔ)的均值表示兩個(gè)時(shí)刻的填補(bǔ)數(shù)據(jù)。4種方法的填補(bǔ)結(jié)果如圖4所示,誤差如表1所示。

圖3 多變量與單變量填補(bǔ)的差異度Fig.3 Difference between the results of multivariate and univariate filling methods

圖4 針對(duì)alive_27的填補(bǔ)結(jié)果Fig.4 Filling results for alive_27

3.1.2 針對(duì)answer _27的填補(bǔ)

類(lèi)似地,對(duì)answer_27進(jìn)行填補(bǔ),其中,多變量填補(bǔ)γ=38 700、σ=120;單變量填補(bǔ)k=10、γ=130、σ=152;組合權(quán)重填補(bǔ)ω=0.6;根據(jù)圖3(b)針對(duì)answer_27的差異度,組合閾值填補(bǔ)t=1×10-3。填補(bǔ)結(jié)果如圖5所示。

圖5 針對(duì)answer_27的填補(bǔ)結(jié)果Fig.5 Filling results for answer_27

從alive_27、answer_27的填補(bǔ)結(jié)果可以看出,單變量填補(bǔ)體現(xiàn)了缺失數(shù)據(jù)的變化趨勢(shì);多變量填補(bǔ)能夠體現(xiàn)波動(dòng)情況。與組合權(quán)重方法相比,組合閾值方法確定閾值t,更加明確,且填補(bǔ)誤差較小。

3.2 缺失數(shù)據(jù)段長(zhǎng)度對(duì)填補(bǔ)的影響

為了進(jìn)一步比較4種方法的填補(bǔ)效果,增加缺失數(shù)據(jù)段的長(zhǎng)度。分別刪除序列alive_27中第6個(gè)變量維度上的最后20、30、40個(gè)數(shù)據(jù)。填補(bǔ)結(jié)果如圖6所示,填補(bǔ)誤差如表2所示,對(duì)應(yīng)的參數(shù)設(shè)置如表3所示。

隨著缺失數(shù)據(jù)段長(zhǎng)度的增加,4種方法的填補(bǔ)效果同上文類(lèi)似,但多變量填補(bǔ)在部分?jǐn)?shù)據(jù)段上的填補(bǔ)偏差變大。原因在于,單變量填補(bǔ)方法在缺失數(shù)據(jù)所屬的變量維度上,以k個(gè)點(diǎn)為初始值進(jìn)行填補(bǔ),能夠較好地保持序列趨勢(shì);但每次填補(bǔ)都以上次填補(bǔ)的結(jié)果作為輸入、存在累積誤差,難以體現(xiàn)波動(dòng)細(xì)節(jié)。多變量填補(bǔ)方法利用其他多個(gè)相關(guān)變量信息,每次填補(bǔ)都以已知數(shù)據(jù)作為輸入,不存在累積誤差,能夠體現(xiàn)波動(dòng)細(xì)節(jié);但沒(méi)有利用缺失數(shù)據(jù)所屬變量信息,在某些點(diǎn)上存在較大偏差。

圖6 缺失數(shù)據(jù)段長(zhǎng)度對(duì)填補(bǔ)效果的影響Fig.6 Relation between the length of missing data segment and its filling effect

表3 4種方法的參數(shù)設(shè)置

相比之下,組合權(quán)重、組合閾值填補(bǔ)的精度較高,受缺失數(shù)據(jù)段長(zhǎng)度的影響不大。這是因?yàn)閮煞N方法利用單變量和其他相關(guān)多變量信息,較好地體現(xiàn)了波動(dòng)細(xì)節(jié),也有效避免了較大偏差。

然而,組合權(quán)重填補(bǔ)對(duì)單變量和多變量填補(bǔ)結(jié)果加權(quán)平均,用一個(gè)權(quán)重值處理所有的缺失數(shù)據(jù),難以體現(xiàn)各個(gè)缺失數(shù)據(jù)的差異性;此外,權(quán)重值不易確定也限制了其使用。組合閾值填補(bǔ)以多變量填補(bǔ)結(jié)果為基礎(chǔ),當(dāng)多變量與單變量填補(bǔ)的差異度較大時(shí),以?xún)烧叩木底鳛樘钛a(bǔ)結(jié)果,既體現(xiàn)了填補(bǔ)的差異性,也避免了權(quán)重值的選擇。

4 結(jié) 論

在小樣本情況下,針對(duì)多元時(shí)間序列中的缺失數(shù)據(jù),利用缺失數(shù)據(jù)所屬單變量和其他相關(guān)多變量信息,提出了一種基于LSSVM的組合閾值填補(bǔ)方法,并進(jìn)行了實(shí)驗(yàn)驗(yàn)證。結(jié)果表明,它能較好地體現(xiàn)序列的波動(dòng)細(xì)節(jié),且受缺失數(shù)據(jù)段長(zhǎng)度的影響相對(duì)較小。

與組合權(quán)重填補(bǔ)方法相比,組合閾值填補(bǔ)方法的參數(shù)(差異度閾值)較容易確定,且能較好地體現(xiàn)填補(bǔ)的差異性,填補(bǔ)精度較高。

[2] 韓敏, 許美玲, 任偉杰. 多元混沌時(shí)間序列的相關(guān)狀態(tài)機(jī)預(yù)測(cè)模型研究[J]. 自動(dòng)化學(xué)報(bào), 2014, 40(5): 822-829.

HAN M, XU M L, REN W J. Research on multivariate chaotic time series prediction using mRSM model[J]. Acta Automatica Sinica, 2014, 40(5): 822-829.

[3] 吳虎勝, 張鳳鳴, 張超, 等. 多元時(shí)間序列的相似性匹配[J]. 應(yīng)用科學(xué)學(xué)報(bào), 2013, 31(6): 643-649.

WU H S, ZHANG F M, ZHANG C, et al. Matching similar patterns for multivariate time series[J]. Journal of applied sciences, 2013, 31(6): 643-649.

[4] PREE H, HERWIG B, GRUBER T, et al. On general purpose time series similarity measures and their use as kernel functions in support vector machines[J]. Information Sciences, 2014, 281(4): 478-495.

[5] LI H L,GUO C H.Piecewise cloud approximation for time series mining[J]. Knowledge-Based Systems, 2011, 24(4): 492-500.

[6] 馬茜,谷峪,李芳芳,等.順序敏感的多源感知數(shù)據(jù)填補(bǔ)技術(shù)[J].軟件學(xué)報(bào),2016,27(9): 2332-2347.

MA Q, GU Y, LI F F, et al. Order-sensitive missing value imputation technology for multi-source sensory data[J]. Journal of Software, 2016, 27(9): 2332-2347.

[7] 呂政, 趙珺, 劉穎, 等. 基于最大方差權(quán)信息系數(shù)的煤氣數(shù)據(jù)填補(bǔ)[J]. 控制理論與應(yīng)用, 2015, 32(5): 646-654.

Lü Z, ZHAO J, LIU Y, et al. Missing data imputation based on maximal variance weight information coefficient for gas flow in steel industry[J].Control Theory & Applications,2015,32(5):646-654.

[8] 雷春麗, 芮執(zhí)元. 基于多元自回歸模型的電主軸熱誤差建模與預(yù)測(cè)[J]. 機(jī)械科學(xué)與技術(shù), 2012, 31(9): 1526-1529.

LEI C L, RUI Z Y. Thermal error modeling and forecasting based on multivariate autoregressive model for motorized spindle[J].Mechanical Science and Technology for Aerospace Engineering, 2012, 31(9): 1526-1529.

[9] ZOUCAS F A M, BELFIORE P. An empirical analysis of a neural network model for the time series forecasting of different industrial segments[J]. International Journal of Applied Decision Sciences, 2015, 8(3): 261-283.

[10] VAPNIK V. The nature of statistical learning theory[M]. NewYork: Springer-Verlag, 1995.

[11] 許磊, 張鳳鳴. 缺失飛參數(shù)據(jù)填補(bǔ)的組合方法研究[J]. 計(jì)算機(jī)工程與應(yīng)用, 2010, 46(21): 210-212.

XU L, ZHANG F M. Research on combined method in missing flight data complement[J]. Computer Engineering and Applications, 2010, 46(21): 210-212.

[12] CHEN T T, LEE S J. A weighted LS-SVM based learning system for time series forecasting[J]. Information Sciences, 2015, 299(C): 99-116.

[13] 謝川,倪世宏,張宗麟.一種缺失飛行參數(shù)預(yù)處理的新方法[J].計(jì)算機(jī)仿真,2005,22(4): 27-30.

XIE C, NI S H, ZHANG Z L. A new method for preprocessing in absent flight parameter[J].Computer Simulation,2005,22(4): 27-30.

[14] 楊軻, 張曉豐, 趙錄峰, 等. 基于LSSVM的缺失飛行數(shù)據(jù)組合填補(bǔ)方法[J]. 火力與指揮控制, 2013, 38(1): 84-90.

YANG K, ZHANG X F, ZHAO L F, et al. Combined method of complementing missing flight data based on LSSVM[J]. Fire Control & Command Control, 2013, 38(1): 84-90.

[15] LANGONE R, ALZATE C, KETELAERE B D, et al. LS-SVM based spectral clustering and regression for predicting maintenance of industrial machines[J]. Engineering Applications of Artificial Intelligence, 2015, 37(37): 268-278.

[16] MOHAMMED W K. High-quality recordings of australian sign language signs[EB/OL].[2017-01-20].http :∥kdd.ics.uci.edu/databases/ High-quality Australian Sign Language/ High-quality Australian Sign Language. html.

[17] KEOGH E, RATANAMAHATANAC. Exact indexing of dynamic time warping[J]. Knowledge and Information Systems, 2005, 7(3): 358-386.