夏心蕾,劉俊陽,彭漢章,韓翔宇

北京航天自動控制研究所,北京 100854

0 引言

在航天、電信、醫(yī)療等多個領(lǐng)域中,曲線類型數(shù)據(jù)都有著廣泛應(yīng)用。如何挖掘曲線類型數(shù)據(jù)是一項重要研究內(nèi)容。在航天領(lǐng)域,飛行器的各種數(shù)據(jù)是判斷飛行器狀態(tài)是否正常的一項重要依據(jù)。曲線數(shù)據(jù)的正確性判斷通常需要專業(yè)人員快速給出結(jié)論,因而不僅對專業(yè)水平要求較高,也存在一定的人為誤判、漏判風(fēng)險,所以自動判讀曲線數(shù)據(jù)正確性對提高航天試驗的判讀效率具有重要意義。

時間序列判讀方法主要分為基于統(tǒng)計特征的判讀方法[1-2]、基于相似度度量的判讀方法[3-4]和人工智能判讀方法[5-6]。基于統(tǒng)計特征的判讀方法通過構(gòu)建數(shù)據(jù)的統(tǒng)計分布模型,檢測時間序列異常。該方法難點在于數(shù)據(jù)分布特性的識別與描述[1]。基于相似度的判讀方法通過計算序列間的距離檢測異常。常用的相似度度量方法有歐氏距離、弗雷歇距離、DTW(Dynamic Time Warping)距離等。歐式距離是最常用的計算距離的方法,計算方法簡單、時間成本較低,但無法處理長度不同的時間序列。弗雷歇距離和DTW距離均可對長度不同的時間序列進行處理,弗雷歇距離求解的是兩個時間序列的最大距離而DTW距離求解的是累積距離。文獻[3]對馬氏距離、DTW距離、夾角距離等多種相似性度量方法在衛(wèi)星遙測數(shù)據(jù)異常檢測中的效果展開研究,驗證了基于DTW距離的異常檢測方法效果較好,能夠識別差異較小的異常序列。文獻[4]提出先使用SAX降低序列長度,再使用Fast-DTW進行衛(wèi)星異常檢測。但這些研究中沒有將周期數(shù)據(jù)和非周期數(shù)據(jù)進行區(qū)分處理,當(dāng)檢測周期數(shù)據(jù)時,無法識別出周期信號的幅度這個重要特征的正確性,尚不能完全滿足周期信號的判讀要求。人工智能判讀方法的典型算法包括最近鄰、支持向量機、神經(jīng)網(wǎng)絡(luò)等。這類算法需要對每種類型曲線進行訓(xùn)練,對于繁多的曲線類型,工作量巨大,并且無法完全保證被測數(shù)據(jù)集分析結(jié)果的準(zhǔn)確性,因此在安全關(guān)鍵領(lǐng)域應(yīng)用難度較大[5],應(yīng)用廣度受到限制。

由于某航天試驗數(shù)據(jù)中周期和非周期曲線的判讀方法不同,僅采用同一種相似度度量方法,會帶來較大判讀偏差,因此需要首先對周期曲線和非周期曲線進行分類。文獻[4]采用基于極值的分段方法對偽周期時間序列進行分割。文獻[7-8]提出基于趨勢的符號化表示(TSX)分割方法,從而支持相似搜索。但上述這些研究是針對偽周期時間序列的周期估計,并不能很好地識別非周期曲線。

本文首先說明曲線數(shù)據(jù)的周期性檢測方法,再針對周期和非周期兩類曲線分別構(gòu)建相似度分析方法,通過分析與模板序列的相似度,預(yù)估曲線的正確性,其流程如圖1所示。本文使用符號近似方法SAX和DTW距離計算相融合的方法用于周期性檢測,對周期和非周期數(shù)據(jù)進行分類;使用帶有雙維度約束的DTW距離計算方法計算非周期信號相似度,使用與無約束DTW距離對比的方法判別曲線是否相似;最后提出了曲線模板更新方法,自動生成標(biāo)準(zhǔn)曲線。

圖1 曲線數(shù)據(jù)的相似度分析流程

1 符號化表示與DTW距離

1.1 符號化表示

符號集合近似(SAX)是一種將時間序列離散化為字符串的方法。該方法首先假設(shè)時間序列服從正態(tài)分布,通過等分正態(tài)分布(用字母表示)劃分時間序列值空間,將落在對應(yīng)分布空間的均值映射為對應(yīng)的字母。算法思路分為2步:1)標(biāo)準(zhǔn)化轉(zhuǎn)換并利用分段聚合近似(PAA)進行數(shù)據(jù)降維;2)將PAA的結(jié)果用字符串表示。

標(biāo)準(zhǔn)化曲線數(shù)據(jù)使其均值變?yōu)?、標(biāo)準(zhǔn)差變?yōu)?。假設(shè)時間序列A=[a1,a2,…,an],標(biāo)準(zhǔn)化后為A′,需要將其轉(zhuǎn)換為長度為w的字符串S=[s1,s2,…,sw]。PAA算法首先將A′切分為w個片段,每個片段采用取均值的方法求解,即將A′轉(zhuǎn)換為B=[b1,b2,…,bw]。

(1)

根據(jù)正態(tài)分布生成α個字符組成的字符集,每個字符對應(yīng)一個一維區(qū)間,每個區(qū)間對應(yīng)的正態(tài)分布概率相等。將PAA降維后的序列B中的每一個片段映射到字符區(qū)間上,從而轉(zhuǎn)換成字符串表示。

符號空間距離度量方法認為相同字符的距離為0,不同字符的距離等于所對應(yīng)區(qū)間較大的字符對應(yīng)的區(qū)間下界減去所對應(yīng)區(qū)間較小的字符對應(yīng)的區(qū)間上界。

1.2 DTW距離

DTW計算距離的主要思路為彎曲時間軸,通過尋找兩條時間序列的最佳對應(yīng)關(guān)系,獲得最小累積距離。

時間序列A=[a1,a2,…,an],時間序列B=[b1,b2,…,bm],令矩陣Mb記錄基距離,即Mb(i,j)表示ai和bj的基距離,度量方式可采用曼哈頓距離。若Mb中路徑W=[w1,w2,…,wK]滿足:1)max(m,n)≤K≤m+n-1;2)W起始于w1=(0,0),終止于wK=(n,m);3)對于wi,1≤i≤K-1,令wi=(xi,yi),wi+1=(xi+1,yi+1),若0≤xi+1-xi≤1,0≤yi+1-yi≤1,則稱路徑W為彎曲路徑。DTW算法的任務(wù)是在所有彎曲路徑中找一條累積距離最小的路徑,DTW距離是這條最小路徑的累積距離,計算方法如下:

(2)

實現(xiàn)中使用動態(tài)規(guī)劃對DTW距離進行求解。定義矩陣Md來表示時間序列A和B的最小累積距離,假設(shè)Md的i行j列元素的值為γi,j:

(3)

γi,j=Mb(i,j)+min{γi-1,j,γi,j-1,γi-1,j-1}

(4)

dDTW(A,B)=γn,m

(5)

其中:i和j的取值范圍分別為0

DTW算法主要存在兩個問題:算法時間復(fù)雜度較高,處理大規(guī)模時間序列的相似度計算速度較慢;算法存在奇異值點,即時間序列A中多個點與時間序列B中的同一個點對齊,這將導(dǎo)致雖然兩條時間序列差異較大,但計算出的DTW距離卻很小。在DTW算法中加入全局約束的方法可以避免奇異值點問題,降低算法時間復(fù)雜度。

2 周期性曲線綜合檢測和分析

2.1 方法概述

周期性檢測可以通過SAX方法,將時間序列轉(zhuǎn)換為字符串,通過字符串的周期性間接判斷時間序列的周期性,但由于SAX中采用了PAA方法對數(shù)據(jù)降維,損失了時間序列的精度,因此這種判斷方法精度較低,存在誤將非周期時間序列判斷為周期時間序列的可能。周期性檢測還可以通過DTW方法判斷時間序列各段周期是否相似,從而判斷曲線周期性。但由于時間序列的周期未知,需要逐個值代入測試,而周期可能的取值范圍大,導(dǎo)致該方法時間復(fù)雜度高。綜合考慮上述情況,本文提出雙重周期性檢測方法,結(jié)構(gòu)圖如圖2所示。第1步周期性檢測采用SAX方法,對字符串的周期性進行檢測,從而初步排除非周期曲線并預(yù)估出周期; 第2步檢測實際上是在第1步檢測結(jié)果的基礎(chǔ)上再進一步進行周期性確認,即使用第1步檢測預(yù)估的周期,采用DTW方法測量每兩個相鄰周期時間序列的相似性。只有兩步檢測結(jié)果都為有周期性的時間序列才被確認為周期信號。

圖2 雙重周期性檢測結(jié)構(gòu)圖

2.2 字符串周期性預(yù)測方法

字符串周期預(yù)測計算的主要思路為先將原字符串進行不同長度lshift的循環(huán)移位,分別計算移位后的字符串與原字符串的距離dstr(無量綱);然后根據(jù)所記錄的不同移位長度對應(yīng)的距離,判斷是否存在T,使得距離在移位長度為c×T,(c=1,2,…)時循環(huán)達到極小值,曲線周期性震蕩。若存在,則預(yù)測T為字符串周期。周期曲線滿足上述規(guī)律,如圖3(a);非周期曲線無此規(guī)律,如圖3(b)。極小值判斷設(shè)定閾值θxtrm,當(dāng)該點小于其前后點的值,并且差值大于θxtrm時,認為該點為極小值點。

圖3 字符串不同移位長度對應(yīng)的距離

2.3 時間序列的周期相似度分析

通過DTW距離測量某一時間序列每兩個相鄰周期的相似度,通過周期間的相似度進一步確認該時間序列是否具有周期性。將字符串周期預(yù)測得到的字符串周期與每個字符代表的時間序列長度相乘,得到時間序列預(yù)測周期。設(shè)時間序列預(yù)測周期為Tts,則

(6)

將時間序列向前移動一個預(yù)測周期的長度,則向前移動后的時間序列長度為lts-Tts。將向前移動后的時間序列和原序列前l(fā)ts-Tts長度的序列進行比較,使用DTW距離度量兩者的相似度。只有當(dāng)該距離小于設(shè)定閾值,認為該時間序列為周期序列。

2.4 周期曲線的頻域特征值計算

在上述周期性分析中對序列進行了降維,因此得出的預(yù)測周期,往往不是原時間序列的真實周期,可能為真實周期倍數(shù)。周期曲線的真實周期值需要通過離散傅立葉變換(DFT)對其頻域特征值進行分析獲得。此外,根據(jù)某航天試驗數(shù)據(jù)要求,周期曲線判讀標(biāo)準(zhǔn)為待測數(shù)據(jù)曲線的幅度與理論幅度相差不大于Δ。本文通過對判讀結(jié)果為正確的數(shù)據(jù)序列集的幅度取均值,動態(tài)生成幅度模板,將待測曲線幅頻特征與模板幅頻特征進行比較,通過兩者的一致性來初步判斷周期曲線的正確性。在對曲線相頻特性有分析要求時,則需要進一步分析判斷。

除了通過頻域特征分析周期曲線外,還可以增加對曲線形狀的分析,通過記錄曲線單個周期數(shù)據(jù),生成單周期模板,比較模板與待測曲線的形狀相似性,進一步精確判斷周期信號的正確性。

3 非周期曲線分析方法

3.1 雙維度約束的DTW距離

對于非周期的航天試驗數(shù)據(jù)曲線,允許同一工況下不同的試驗過程之間的測試值在時間序列上伸縮對齊,也就是允許一定程度上時間序列回卷后相似,這種特性的數(shù)據(jù)適合使用DTW距離進行分析。但由于某航天試驗數(shù)據(jù)判讀對測試值在時間以及幅值上的偏移距離均有限制,因此需要在DTW距離算法的基礎(chǔ)上分別加入時間和幅值偏移約束,形成帶雙維度約束的DTW距離計算方法,來滿足本領(lǐng)域曲線的相似性評價標(biāo)準(zhǔn)。

時間約束使用S-C全局約束方法。設(shè)置時間約束值rt,rt為1的S-C全局約束方法如圖4。陰影部分表示彎曲窗口,限制得到的最小累積距離彎曲路徑在這個彎曲窗口內(nèi)。

圖4 S-C全局約束方法

在S-C約束下,加入對幅值的約束。假設(shè)A為模板序列,B為待測序列。設(shè)置數(shù)據(jù)點bi的上下界,根據(jù)S-C全局約束方法,序列A中只有數(shù)據(jù)點集{ai-rt,ai-rt+1,…,ai+rt}可與序列B的第i個數(shù)據(jù)點bi對齊,以數(shù)據(jù)點集的最大值作為上界,最小值作為下界。同時,加入幅值約束松緊度θtyt,動態(tài)調(diào)整幅值約束的松緊度。若滿足式(7),則認為bi滿足幅值約束。

在用動態(tài)規(guī)劃算法計算雙維度約束的DTW距離時,定義矩陣Mdual表示雙維度約束下的最小累積距離矩陣,如式(8)。

通過全局約束限制最小累積距離彎曲路徑的偏移量,可以避免奇異值問題,時間復(fù)雜度可以降為O(n)。

(7)

(8)

3.2 距離分析方法

3.3 非周期時間序列模板生成

由于非周期時間序列分析需要模板曲線,本文根據(jù)計算DTW距離時得出的最小彎曲路徑,動態(tài)更新標(biāo)準(zhǔn)曲線。當(dāng)該待測曲線相似度分析結(jié)果為相似時更新模板。假設(shè)舊模板中第i個點為oi,在最小彎曲路徑中檢測曲線A與oi對齊的所有點的均值為aavg,則新模板中第i個點ui的表達式為:

(9)

式中:n為更新模板的次數(shù)。

4 實驗驗證

為驗證本文所提方法的正確性,使用某運載項目航天試驗數(shù)據(jù)進行了實驗分析,待測曲線示例如圖5所示。下面對3種典型曲線進行分析。

圖5 某航天試驗曲線數(shù)據(jù)示例

4.1 周期性的檢測和分析

對于周期曲線,以曲線1為例,前51個數(shù)據(jù)點的SAX轉(zhuǎn)換過程如圖6。虛線表示通過等分正態(tài)分布劃分值空間,由下往上,每個值區(qū)間分別對應(yīng)字母a至h,前51個數(shù)據(jù)點轉(zhuǎn)換生成字符串“adgccffbehadgccff”。

圖6 SAX轉(zhuǎn)換示意圖(曲線1)

對轉(zhuǎn)換后的字符串進行周期檢測,極小值點判斷閾值設(shè)置為0.5。圖7為不同循環(huán)移位長度對應(yīng)的字符串距離,圓圈標(biāo)注出了距離極小值點,其中實心圓表示周期最小的周期性出現(xiàn)的極小值點,兩個相鄰實心圓的距離即為預(yù)估周期。字符串周期檢測結(jié)果為10,則預(yù)估時間序列的周期為30。

圖7 字符串周期預(yù)測示意圖(曲線1)

周期時間序列DTW相似度計算設(shè)定閾值為0.1,曲線1的DTW距離為0.0408,小于設(shè)定閾值,周期性檢測結(jié)果為周期曲線。通過離散傅里葉變化計算得出時間序列周期為10,與預(yù)期相符。

對于非周期曲線,以曲線2為例,截取前51個數(shù)據(jù)點的SAX轉(zhuǎn)換過程如圖8。轉(zhuǎn)換生成的字符串為“dddddhhhhhhhhhfdd”。對轉(zhuǎn)換后的字符串進行周期檢測,不同移位長度對應(yīng)的字符串距離如圖9,由于沒有滿足周期性出現(xiàn)特征條件的極小值點,檢測結(jié)果為非周期序列,符合預(yù)期。

圖8 SAX轉(zhuǎn)換示意圖(曲線2)

圖9 字符串周期預(yù)測示意圖(曲線2)

4.2 非周期曲線相似度分析

以曲線2為例,如圖5(b),將曲線2的第400點加入幅值為1的脈沖擾動,構(gòu)造得到加擾曲線,不同約束DTW度量方法測量出的被測曲線和模板曲線的距離如表1,其中S-C約束與雙維度約束的時間約束值均為1,雙維度約束的幅值約束松緊度設(shè)置為0.2(θtyt=0.2)。加擾曲線的無約束DTW距離和帶有S-C約束的DTW距離都可能小于正常曲線累積距離閾值,即存在部分曲線相似度分析結(jié)果與實際不符的可能。而雙維度約束DTW距離算法的計算結(jié)果因為有幅值約束,因此結(jié)果為無窮大,從而識別出曲線異常。

表1 曲線2不同約束DTW度量方法結(jié)果

以曲線3為例,如圖5(c),將曲線左移3個單位構(gòu)造得到異常曲線數(shù)據(jù)。采用不同約束DTW度量方法測量出的上述兩條曲線與模板曲線的距離如表2。本試驗原曲線的雙維度約束和無約束DTW距離的比值為5.19,小于θcmlt(本實驗中設(shè)定為10);異常曲線的雙維度約束和無約束DTW距離的比值為13.02,大于θcmlt,偏差顯著。因此,采用本文雙維度約束和無約束DTW距離比較的距離分析方法,可以在一定條件下較為有效地識別由于時間彎曲過大帶來累積距離的異常。

表2 曲線3不同約束DTW度量方法結(jié)果

下面對3種約束DTW度量方法進行耗時統(tǒng)計,如表3,雙維度約束方法計算時長明顯小于無約束方法。對于正常曲線,由于雙維度約束中加入了幅值約束的判斷,因此其計算時長略大于S-C約束方法的計算時長。對于異常曲線,雙維度約束計算時長小于S-C約束方法的計算時長。顯然,雙維度約束算法與S-C約束相當(dāng),比無約束快約一百倍。

表3 不同約束DTW度量方法耗時

5 結(jié)論

本文針對某航天數(shù)據(jù)的時間序列判讀的難點問題,提出了SAX和DTW融合的雙重周期性檢測方法和帶有雙維度約束的DTW非周期相似度計算方法。在某航天試驗數(shù)據(jù)集上進行分析實驗驗證,結(jié)果表明本文提出的分析方法正確有效。