航天試驗數據時間序列相似度分析

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

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

0 引言

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

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

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

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

圖1 曲線數據的相似度分析流程

1 符號化表示與DTW距離

1.1 符號化表示

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

標準化曲線數據使其均值變為0、標準差變為1。假設時間序列A=[a1,a2,…,an],標準化后為A′,需要將其轉換為長度為w的字符串S=[s1,s2,…,sw]。PAA算法首先將A′切分為w個片段,每個片段采用取均值的方法求解,即將A′轉換為B=[b1,b2,…,bw]。

(1)

根據正態分布生成α個字符組成的字符集,每個字符對應一個一維區間,每個區間對應的正態分布概率相等。將PAA降維后的序列B中的每一個片段映射到字符區間上,從而轉換成字符串表示。

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

1.2 DTW距離

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

時間序列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算法的任務是在所有彎曲路徑中找一條累積距離最小的路徑,DTW距離是這條最小路徑的累積距離,計算方法如下:

(2)

實現中使用動態規劃對DTW距離進行求解。定義矩陣Md來表示時間序列A和B的最小累積距離,假設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算法主要存在兩個問題:算法時間復雜度較高,處理大規模時間序列的相似度計算速度較慢;算法存在奇異值點,即時間序列A中多個點與時間序列B中的同一個點對齊,這將導致雖然兩條時間序列差異較大,但計算出的DTW距離卻很小。在DTW算法中加入全局約束的方法可以避免奇異值點問題,降低算法時間復雜度。

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

2.1 方法概述

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

圖2 雙重周期性檢測結構圖

2.2 字符串周期性預測方法

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

圖3 字符串不同移位長度對應的距離

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

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

(6)

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

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

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

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

3 非周期曲線分析方法

3.1 雙維度約束的DTW距離

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

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

圖4 S-C全局約束方法

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

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

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

(7)

(8)

3.2 距離分析方法

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

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

(9)

式中:n為更新模板的次數。

4 實驗驗證

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

圖5 某航天試驗曲線數據示例

4.1 周期性的檢測和分析

對于周期曲線,以曲線1為例,前51個數據點的SAX轉換過程如圖6。虛線表示通過等分正態分布劃分值空間,由下往上,每個值區間分別對應字母a至h,前51個數據點轉換生成字符串“adgccffbehadgccff”。

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

對轉換后的字符串進行周期檢測,極小值點判斷閾值設置為0.5。圖7為不同循環移位長度對應的字符串距離,圓圈標注出了距離極小值點,其中實心圓表示周期最小的周期性出現的極小值點,兩個相鄰實心圓的距離即為預估周期。字符串周期檢測結果為10,則預估時間序列的周期為30。

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

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

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

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

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

4.2 非周期曲線相似度分析

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

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

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

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

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

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

5 結論

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