夏心蕾,劉俊陽(yáng),彭漢章,韓翔宇

北京航天自動(dòng)控制研究所,北京 100854

0 引言

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

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

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

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

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

1 符號(hào)化表示與DTW距離

1.1 符號(hào)化表示

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

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

(1)

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

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

1.2 DTW距離

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

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

(2)

實(shí)現(xiàn)中使用動(dòng)態(tài)規(guī)劃對(duì)DTW距離進(jìn)行求解。定義矩陣Md來(lái)表示時(shí)間序列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算法主要存在兩個(gè)問(wèn)題:算法時(shí)間復(fù)雜度較高,處理大規(guī)模時(shí)間序列的相似度計(jì)算速度較慢;算法存在奇異值點(diǎn),即時(shí)間序列A中多個(gè)點(diǎn)與時(shí)間序列B中的同一個(gè)點(diǎn)對(duì)齊,這將導(dǎo)致雖然兩條時(shí)間序列差異較大,但計(jì)算出的DTW距離卻很小。在DTW算法中加入全局約束的方法可以避免奇異值點(diǎn)問(wèn)題,降低算法時(shí)間復(fù)雜度。

2 周期性曲線綜合檢測(cè)和分析

2.1 方法概述

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

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

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

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

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

2.3 時(shí)間序列的周期相似度分析

通過(guò)DTW距離測(cè)量某一時(shí)間序列每?jī)蓚€(gè)相鄰周期的相似度,通過(guò)周期間的相似度進(jìn)一步確認(rèn)該時(shí)間序列是否具有周期性。將字符串周期預(yù)測(cè)得到的字符串周期與每個(gè)字符代表的時(shí)間序列長(zhǎng)度相乘,得到時(shí)間序列預(yù)測(cè)周期。設(shè)時(shí)間序列預(yù)測(cè)周期為T(mén)ts,則

(6)

將時(shí)間序列向前移動(dòng)一個(gè)預(yù)測(cè)周期的長(zhǎng)度,則向前移動(dòng)后的時(shí)間序列長(zhǎng)度為lts-Tts。將向前移動(dòng)后的時(shí)間序列和原序列前l(fā)ts-Tts長(zhǎng)度的序列進(jìn)行比較,使用DTW距離度量?jī)烧叩南嗨贫?。只有?dāng)該距離小于設(shè)定閾值,認(rèn)為該時(shí)間序列為周期序列。

2.4 周期曲線的頻域特征值計(jì)算

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

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

3 非周期曲線分析方法

3.1 雙維度約束的DTW距離

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

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

圖4 S-C全局約束方法

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

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

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

(7)

(8)

3.2 距離分析方法

3.3 非周期時(shí)間序列模板生成

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

(9)

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

4 實(shí)驗(yàn)驗(yàn)證

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

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

4.1 周期性的檢測(cè)和分析

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

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

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

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

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

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

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

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

4.2 非周期曲線相似度分析

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

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

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

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

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

表3 不同約束DTW度量方法耗時(shí)

5 結(jié)論

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