基于趨勢符號聚合近似的衛星時序數據分類方法

阮輝，劉雷，胡曉光，*

（1.北京航空航天大學 自動化科學與電氣工程學院，北京100083； 2.北京機電工程研究所，北京100074）

時間序列是按照時間排序的一組隨機變量，其通常是在相等間隔的時間段內依照給定的采樣率對某種潛在過程進行觀測的結果［1］。在衛星的測控管理過程中，會產生大量的遙測數據，它們以時間序列的形式存儲在數據庫中。而運行狀態監測系統傳感器產生的監測數據通過遙測系統傳輸至地面控制中心，此類數據是地面判斷在軌衛星運行和健康狀態的唯一依據［2］。同時，這些海量的時間序列蘊含可用于衛星故障診斷的規律和知識。通過數據挖掘，可以提取衛星各器件的信息，發現異常、關聯、模式、趨勢等知識。有效掌握和利用這些信息規律對衛星異常檢測和關聯分析、故障診斷、監測預警，對衛星測管管理與決策活動，如改進衛星設計、提高測試及監測自動化等工作具有特別重要的意義［3-7］。

此外，在 金 融［8］、天 氣 觀 察［9］、生 物 醫 學 測量［10］領域同樣大量產生此類型的數據。因此，在過去幾十年中，時間序列數據挖掘的研究在開發新算法和改進可用算法以滿足當前需求方面一直非?；钴S。

由于原始時間序列具有高維、高音量和大量噪聲的特征，計算機要對原始時間序列進行分類面臨許多挑戰。

現有的分類算法可以分為2類，即基于形狀的分類算法和基于結構的分類算法。

基于形狀的分類算法以最近鄰（One-Nearest-Neighbor，1-NN）分類器作為基礎，時序數據的表示方法則為歐氏距離（Euclidean Distance，ED）［11］和動態時間規整（Dynamic Time Warping，DTW）［12］。

在短時間序列數據集上，1-NN DTW 已被證明是具有高度代表性和競爭力的分類器。基于形狀的技術的不足之處在于：當使用噪聲或包含特征子結構的長數據對數據進行分類時，它們顯示的性能很差。

基于結構的技術有2個步驟。首先，這些方法包括離散傅里葉變換（Discrete Fourier Transform，DFT）［13］，可索引分段線性逼近［14］用于提取特征向量；然后，采用經典數據挖掘算法（如決策樹），支持向量機被用來做分類任務。

目前，符號化方法也被廣泛使用。因為它具有簡單、可讀性、高效率之外，還可以使用其他領域的算法，如信息檢索或文本處理等。其中，Lin等［15］提出的符號聚合近似（Symbolic Aggregate Approximation，SAX）方法是時下流行的符號化表示方法。SAX基于分段總體逼近（Piecewise Aggregate Approximation，PAA）方法［16］，并假定PAA值遵循高斯分布。SAX根據高斯曲線下的等大小區域離散化PAA值，從而產生斷點，后采用計算開銷小的低邊界計算準確的距離。

模式 袋 模 型（Bag-of-Patterns，BOP）［17］提 取子結構作為時間序列的高級特征，將這些子結構轉化為SAX，并采用ED作為度量距離進行相似性度量的相關應用，如分類、聚類等。BOP方法不考慮子結構之間的順序，只是將時間序列看成是一些SAX字符串出現概率的集合，每個SAX字符串是相互獨立的，類似于直方圖的統計表示。

但SAX表示的質量取決于PAA系數，即時間序列分割的段數，用于量化的符號數量（字母大小數）和高斯假設。同時，SAX中的符號是根據每個分段的均值繪制的，不能表示分段的趨勢。

為此，一些文獻嘗試解決SAX的這些問題。Pham等［18］通過引入對時間序列截的長度和字母大小起作用的自適應斷點矢量來緩解高斯假設。但是，通過使用聚類方法引入預處理階段，就失去了SAX的簡單性。為描述SAX的趨勢信息，文獻［19］中提出的擴展符號聚合近似（Extended SAX，ESAX）將PAA段的符號最小值和最大值與相關的SAX符號以及它們的出現順序相關聯。這定義了一個抽象形狀，可以在時間序列檢索中提供更好的結果。但是，ESAX 表示的大小是SAX表示的大小的3倍。從效率的角度來看，筆者沒有將文獻［19］的方法與相同大小的SAX表示進行比較。文獻［20］提出的基于趨勢的符號聚合逼近（Trend-based Symbolic Aggregate approXimation，TSAX），在每個分段中添加了2個趨勢指標，TSAX可以區分2個不同趨勢之間的差異，但仍不能反映同一趨勢中的差異。另外，TSAX表示的大小與ESAX相同，是SAX的3倍。

提高分類算法準確率、減少運算時間是數據符號化表示的主要目標。本文提出一種改進的符號表示——趨勢符號聚合近似（Trend Symbolic Aggregate Approximation，TrSAX），集成SAX與最小二乘法，用以描述時間序列的均值和斜率，在不增加運算時間的情況下，進一步提高了分類的準確率。

1 衛星時序數據

1.1 衛星遙測數據來源

根據任務的不同，一顆衛星所含分系統略有差別，但通常包括熱控分系統、姿軌控分系統、星務管理分系統、電源分系統、推進分系統、結構分系統、測控分系統、總體電路分系統等。

為監測各分系統的工作狀態，設計時在各分系統中設置了大量的測試傳感器，采集測試信息傳遞到地面，形成衛星遙測數據。

1.2 衛星遙測數據特點

衛星的測試數據一般分為數字量和模擬量。數字量可細分為獨立數字量、關聯數字量和狀態數字量。模擬量可細分為恒定模擬量、區間模擬量和趨勢變化模擬量。其中，數字量主要為指令、計數和狀態等，模擬量主要為電流、電壓、角度、溫度、壓力等。此外，由于衛星是按照軌道周期運行，衛星遙測數據存在周期性。

綜上，不同衛星不同分系統的具體遙測數據需要在進行充分分析、數據預處理之后，進行分類、聚類、異常檢測等相關研究。

1.3 真實衛星遙測數據分析

文獻［21］中對“風云三號”衛星遙測數據進行了分析，其中含133個測試參量，時間跨度為500天。并對3種具有周期特性的測試序列（電流、角度、轉速）為對象進行了詳細研究。

由于衛星遙測數據具有明顯的周期性，通過對遙測數據的每個周期進行分析，可知衛星在該周期之內的運行狀態是否正常。

首先，需要對遙測數據進行分段，具體根據衛星工作方式進行分段，根據幅角測試參量的測試值為0～360不斷循環，具有周期性，且周期分隔點明確。同時，幅角測試值能夠反映出衛星軌道的運行位置，因此，將測試值由360跳變到0的點作為幅角突變點進行分段。對于因數據缺失而不完整的序列，在實驗中進行剔除，確保各分段子序列完整。

如圖1～圖3所示，以角度測試時間序列的幅角突變點為標識進行分段，將每段序列進行疊加繪圖。樣本數據樣本選用的是3種遙測數據測試序列的前4 000個有效樣本點，每個周期的數據點為86，樣本數據包含47個周期。圖中顯示各個分段子序列之間的耦合度高，分段合理。

圖1 衛星某角度測試時間序列分段疊加結果Fig.1 Satellite angle test time series segmentation merging results

圖2 衛星某轉速測試時間序列分段疊加結果Fig.2 Satellite rotation speed test time series segmentation merging results

圖3 衛星某電流測試時間序列分段疊加結果Fig.3 Satellite electric current test time series segmentation merging results

圖中分段子序列具有以下特點：

1）各子序列整體變化趨勢相同，局部有細小差別，波形不完全重合。

2）測試子序列波動趨勢不同，有的較為平滑如圖2所示，有的較為激烈如圖3所示。

2 相關定義

2.1 時間序列

時間序列數據［15］是指將使用同一統計指標的數值按照時間先后順序排列而成的數列：

式中：時間序列T的長度為n。

2.2 截

使用滑動窗口函數［17］，將時間序列T=（t1，t2，…，tn）分為固定長度為ω的多個截Si；ω=（ti，ti+1，…，ti+ω-1）。

式中：時間序列T分成n-ω+1個截。

2.3 段

截Si；ω=（ti，ti+1，…，ti+ω-1）被分成m個等長的段，表 示 為 集 合｛sj，j=1，2，…，m｝，sj=｛t（j-1）×L+l，l=1，2，…，L｝，L=floor（ω/m），floor函數為向下取整函數。

圖4為顯示了時間序列劃分的示意圖。長度為512的時間序列樣本被分為449個截，每個截被等分為4段。

圖4 時間序列劃分示意圖Fig.4 Schematic diagram of time series division

3 符號化表示

3.1 段平均值

段sj=｛t（j-1）×L+l，l=1，2，…，L｝的平均值ˉsj計算如下：

表1列出了字母數為3～9的斷點，斷點字母符號的范圍為3～5。使用φ-1個斷點將分布空間劃分為φ個等概率區域。N（0，1）高斯曲線下的B=β1，β2，…，βi的斷點間面積等于1/φ。

表1 字母數為3～9的斷點查找表Table 1 Look up table from breakpoints with alphabet sizes from 3 to 5

圖5 φ=3和φ=4時，各平均值和小寫字母之間的對應關系Fig.5 Corresponding relationship between average values and lower-case letters for φ=3 and φ=4

3.2 段斜率值

SAX忽略了時間序列段中對于分析時間序列分類和相似性至關重要的趨勢信息這一重要特征。為了描述趨勢信息，本文使用最小二乘法來計算每個時間序列的斜率值。

式中：L=floor（ω/m），為一個截中段的個數；j表示第j個段。

采用0°和±45°三個角度的斜率作為角度間隔值將角度空間（-90°，90°）劃分為5個非重疊間隔，如圖6所示，將計算出的斜率用大寫字母表示，該字母對應于它們的駐留區域。通過這種方法，將角度空間（-90°，90°）轉換為5個大寫字母，代表5種情況：速降（A）、緩降（B）、水平（C）、緩增（D）和速增（E）。

每段用2個符號表示，大寫字母表示斜率值，小寫字母表示平均值。將計算出的大寫字母和小寫字母組合起來代表每個時間序列截。這些為每個時間截組成的大小字母串，稱為趨勢符號聚合近似（TrSAX）表示。圖7為長度64的時間序列截的TrSAX字母串，對于φ=4，截的TrSAX字為AcAaDaEa。

圖6 角度間隔和字母的對映關系Fig.6 Angle interval and corresponding letters

圖7 趨勢符號聚合近似（TrSAX）表示示意圖Fig.7 Schematic diagram of trend symbolic aggregate approximation（TrSAX）representation

4 TrSAX詞袋

與BOP方法類似，本文提取時間序列中的截作為子序列，并將其轉換為TrSAX詞，則一條時間序列可轉為一個TrSAX詞袋（Bag-of-TrSAX，BOTS）。

4.1 BOTS參數設定

BOTS有4個參數，分別為截長度ω和3個用于表示截的TrSAX參數，即字母符號數φ、該截中段的數量m和大寫字母數γ。

4.2 BOTS工作流程

1）使用4.1節的TrSAX參數，為數據集建立一個大小為s=（γφ）m的TrSAX“詞匯表”。

2）使用滑動窗口函數，將長度為n的原始時間序列T轉換為n-ω+1個固定大小的截。

3）將這些時間序列截進行（0，1）標準化。

4）將時間截轉換為TrSAX字。

5）時間序列T轉換成為一個TrSAX字袋。

如果時間序列沒有噪聲或尖峰，是平滑的，則許多連續的截將轉換為相同的TrSAX字。若對所有這些出現的TrSAX單詞進行計數，則計數過多。

為此，為了防止這種情況，本文使用數據規約［17］，僅統計重復出現的TrSAX詞中第一個，并忽略所有重復項，直到新詞出現。數據規約如下所示，首次出現的TrSAX詞標粗。

假設獲得如下TrSAX字符串S：

數據規約之后，字符串S迭代為S′：

然后，構造直方圖H以計算出現在字符串S′中的TrSAX單詞的數量。

在每個時間序列獲得直方圖H之后，建立TrSAX字矩陣M，其中每一行對應于時間序列數據，每一列對應于TrSAX字。表2給出了一個時間序列的BOTS表示形式的虛擬示例。

對于獲得的矩陣M，時間序列相似性度量為兩行中列的累加平方值的平方根。給定矩陣Mm×n，Ri=［mi1，mi2，…，min］和 Rj=［mj1，mi2，…，mjn］為矩陣的第i行和第j行。

表2 時間序列的BOTS表示形式的虛擬示例Table 2 Visual example of BOTS representation for time series

兩者從時間序列Ti和Tj轉換而來，由BOTS表示，BOTS距離定義為

式中：

5 實驗分析

5.1 實驗數據集

目前的衛星遙測數據各周期內的正常模式、異常模式、故障模式均沒有明確資料可參考，因而選用具有衛星遙測參數數據類似特性的公開數據集進行驗證。從UCR公開數據集中篩選出與衛星遙測數據中的角度序列類似的SonyAIBORobot-Surface數據集，與衛星遙測數據中的轉速序列類似的Fish數據集，與衛星遙測數據中的電流序列類似的FaceFour數據集進行實驗驗證，3個時間序列數據集的信息如表3所示。

表3 時間序列數據集詳細信息Table 3 Details on time series datasets

5.2 參數訓練

每個數據集包含一個訓練集和一個測試集。采用訓練集訓練參數。如4.1節所述，BOTS分類器需要4個參數。

首先，使用網格搜索方法結合交叉驗證從訓練集中訓練了這些參數。使用固定的大寫字母數γ=5和字母符號數φ=4。對于段的數量m，本文選擇｛3，4，5｝。對于截長度ω，都從16開始，每次加倍，直到小于或等于n/2。訓練算法如下所示。

訓練階段的目的是找到一組參數值，使得訓練的分類準確率最大。如果2組參數具有相同的分類準確率，則選擇較小的參數。在訓練階段的最后，獲得了訓練數據集的最佳參數值集合｛γ，φ，m，ω｝，可直接用在測試集上進行分類。

算法1BOTS分類算法。

From sk.learn.model＿selection import cross＿val＿score

colLengh=len（trainInput［0］）

end=int（math.log（colLengh/2，2））

windowList=list（map（int，np.logspace（4，end，end-3，base=2））） ＃截長度從16開始，每次加倍，直到小于或等于n/2

wordList=［3，4，5］ ＃段的數量

tmpErrorRate，optimalWindow，optimalWord=1，None，None

for window in windowList：

for word in wordList：

testErrorRates=

cross＿val＿score（trainAndPredict（window，word，train Input，trainLabels，testInput，testLabels），X＿trainval，Y＿trainval，cv=5） ＃5折交叉驗證

testErrorRate = testErrorRates.mean（） ＃取平均數

logRecord.append（［window，word，testErrorRate］）

tmpErrorRate，optimalWindow，optimalWord=testErrorRate，window，word ＃得到最優的截長度、段的數量

5.3 TrSAX各項參數的影響

所有的參數值集合在測試集上運行實驗，分類錯誤率結果及運行時間如圖8和圖9所示。

如圖8所示，段的數量m從3增加到4，分類錯誤率隨著截長的變化，出現非線性的變化。但是當m值為5時，分類錯誤率明顯增加。m的值對分類錯誤率有影響，但沒有一定規律。通常，對于本文中的大多數數據集的實驗，值為3或4效果很好。

圖8 三個數據集的分類錯誤率結果Fig.8 Classification Error rate results for three datasets

圖9 三個數據集的運行時間結果Fig.9 Running time results for three datasets

如圖9所示，段的數量m對運行時間的影響有2個規則。首先，m從3增加到5，運行時間明顯增加。原因是TrSAX詞匯量為s=20m，隨m的增加呈指數增長，并且與運行時間成正相關。另外，在m相同的情況下，運行時間隨截長的增加而減少，并且當字典大小較大時，這種情況更加明顯。這是因為詞匯量大且稀疏，并且截長ω的增加減少了映射單詞的數量，同時減小了所獲得的矩陣M 的大小，從而減少了運行時間。

5.4 對比算法及實驗結果

由于本文方法旨在通過增加趨勢信息來提高SAX的分類準確率。因此，將與1-NN SAX（或BOP）進行比較。另外，選擇了ESAX和TSAX作為比較算法。如引言所述，基于ED和DTW 建立的1-NN分類器在短時間序列上具有很高的代表性和競爭力。因此，將本文提出的BOTS與這5種算法進行比較，如表4所示。

通過對3個數據集進行預定義分類的實驗，驗證了BOTS的有效性。實驗結果顯示，BOTS分類算法在3個數據集中分類錯誤率都是最低，表現最佳。在平均排名方面，基本順序是BOTS＞1-NN TSAX＞1-NN SAX＞1-NN DTW ＞1-NN ED＞1-NN ESAX。

表4 不同表示算法的分類結果Table 4 Classification results of different representation algorithms

在所選擇的3個數據集中，1-NN DTW 和1-NN ED表現較好，略低于1-NN SAX，優于1-NN ESAX。而BOTS和1-NN TSAX 表現優于1-NN SAX，而1-NN ESAX表現排名最差，為第6。說明趨勢信息對降低分類錯誤率影響明顯，對SAX增加合適的趨勢信息可以明顯降低分類錯誤率，而不合理的趨勢信息則會明顯增加分類錯誤率。

與SAX相比，TrSAX的優點在于：它可以通過每個段的趨勢特征描述具有相同均值但趨勢不同或趨勢方向相同，但斜率有差異的某些模式（如緩增或速增）來提高分類的準確性。另外，對于某種時間序列，趨勢信息對于專家而言非常重要，并且許多決策是通過趨勢分析確定的。

TrSAX也有一些缺點：它僅是通過0°和±45°3個角度間隔值簡單地將趨勢特征區分為速增、緩增、水平、緩降和速降5個狀態，這對于一些數據集有用。對于一些特殊的數據集，不同的角度間隔會使得分類準確率不同。如何合理確定角度間隔值的問題仍未解決。

6 結 論

1）本文提出了一種基于趨勢符號聚合近似的時間序列表示方法TrSAX，可以將原始時間序列轉換為矩陣結構，該矩陣結構可以使用最小二乘法結合SAX表示來描述每個時間序列段的平均值和斜率值。

2）對文獻［21］中的“風云三號”衛星遙測數據進行了分析，得出角度序列、轉速序列、電流序列3種遙測數據具有明顯的周期性，以角度測試時間序列的幅角突變點為標識進行分段，進行疊加繪圖，可以得到耦合度較高的周期曲線。在UCR公共數據集中篩選出與角度序列、轉速序列、電流序列類似的3個數據集進行了實驗驗證。

3）通過實驗分析段的數量m，截長度ω的各種值對BOTS分類性能的影響。期望較大m和較小的ω來描述原始時間序列的盡可能多的細節，并為時間序列分類帶來更好的結果。但是，實驗結果證明這一期望沒有實現。隨著m的增加或ω的減少，分類錯誤率不會持續下降。為了確定m和ω的最佳組合，需要在提供的訓練數據集上進行訓練。此外，m對分類錯誤率更敏感，并且m比ω具有更大的影響。對于大多數數據集，m等于3或4可以產生令人滿意的結果。

4）對比實驗結果顯示，與衛星遙測參數中角度序列、轉速序列、電流序列類似的3個公開數據集中，BOTS分類錯誤率明顯低于其他5種分類方法，對未來衛星模擬量遙測參數的分類提供了一種新的思路。

在未來的研究中，將研究參數m 和ω的設置，以加速和簡化參數的確定，同時需要對角度間隔值進行進一步的優化。此外，還將探索基于TrSAX的其他任務，如聚類以及相關性分析等。