聚類HMM模型在QAR數據分析中的應用研究

楊 慧 毛好好 霍緯綱

(中國民航大學計算機科學與技術學院 天津 300300)

0 引 言

QAR是飛機標準配置組件，真實地記錄了飛行過程中的各種參數。隨著民航運輸業整體規模的擴大, 航空安全更加引人關注。近年來隨著技術的不斷進步，以QAR為代表的飛行數據記錄設備不斷更新升級，數據采集數量與質量都有很大提高，日益引起各航空運營商的高度重視[1]。針對QAR數據的研究，目前常用的方法是序列相似性匹配、小波變換、決策樹以及關聯規則等。如之前所做的QAR數據多維子序列的相似性搜索，基于 FP-Tree 的 QAR 數據故障檢測研究等。但是這些方法主要關注在異常數據的檢測與診斷，而忽略了故障或異常發生前后的狀態。

隱馬爾可夫模型(HMM)作為動態時間序列統計模型，具有嚴謹的數據結構和可靠的計算性能，早期成為語音識別的主流技術[2]。近幾十年來， 隱馬爾科夫模型被廣泛地應用于各領域中，比如語音識別、行為識別、 生物學、 控制和經濟領域等。黃曉彬等[3]對中國股市進行信息探測， 使用貝葉斯和馬爾科夫鏈蒙特卡洛的方法，驗證了該模型對市場信息的識別能力較強。隨著研究的深入，國內外已開始把 HMM 方法引入到狀態監測和故障診斷領域[4]。于天劍等[5]的HMM在電機軸承上的故障診斷，利用HMM方法與ANFIS方法相結合的情況來對故障問題進行建模，在頻域內提取特征值，用ACC算法進行故障分類、最后用隱馬爾可夫模型方法進行了故障預測。柳姣姣等的基于隱馬爾科夫模型的時空序列預測方法[6]，提出了基于時空密度聚類的隱馬爾科夫模型。本文提出一種基于聚類的HMM模型，針對QAR數據特點，分析發生故障或異常時QAR數據中不同屬性的變化特點，提取主要影響屬性進行分析。通過對其聚類進行數據離散化，并將其分為多個狀態，對故障或異常發生的過程進行HMM建模。以狀態序列的形式描述故障或異常發生的過程，以飛機空中顛簸故障為例，建立空中顛簸故障HMM模型，并檢驗了該模型的有效性。

1 相關工作

1.1 QAR數據

本文的目標是分析QAR數據，由于QAR數據是從傳感器獲取的流數據，其采樣頻率一般為每秒采樣一次，數據量非常巨大。

QAR數據存在著以下幾個特點：(1) 數據量巨大；(2) 時間性很強；(3) 參數種類眾多；(4) 干擾多，隨機性強，具有很多的不確定因素；(5) 需要很強的經驗知識。

1.2 隱馬爾可夫模型(HMM)

一個HMM即為五元組λ={N,M,A,B,π}，N為模型狀態數，M為觀察狀態數。為了簡便，常記為三元組λ={A,B,π}來表示一個隱馬爾可夫模型。

隱馬爾可夫模型(HMM)通常解決三個基本問題，針對三個問題并提出相對應的解決算法：

問題一：評估問題，即某個模型下觀察值序列的概率計算問題。

問題二：解碼問題，即隱馬爾可夫模型的最優狀態序列問題。

問題三：學習問題，即隱馬爾可夫模型各種參數的訓練問題。

HMM能夠對非平穩信號變化的規律進行統計并建立參數化模型。QAR為時序數據，在飛機從正常到發生故障的過程中, 每個階段的QAR數據在不斷地發生變化。若對QAR的這些觀測數據進行處理，通過學習訓練出λ={A,B,π}，用 HMM來統計性地描述這些QAR數據的改變,以狀態序列的形式描述故障的演化過程，并通過解碼，得到最優的狀態趨勢序列，可以更好地識別飛機的狀態趨勢。

以空中顛簸故障為例，空中顛簸是威脅飛行安全的故障之一，而且空中顛簸的發生具有隨機性，每次顛簸的強度也具有隨機性。空中顛簸影響QAR數據不斷發生變化，若把空中顛簸數據的狀態趨勢分成三個等級：嚴重顛簸、輕微顛簸和正常。由于空中顛簸故障發生在飛機飛行途中，地面人員無法直接得知顛簸程度，故可將空中顛簸故障發生時的狀態趨勢作為隱馬爾可夫模型的隱狀態?？罩蓄嶔す收习l生時，QAR數據的多個屬性值受到影響而發生異常，QAR數據是可以直接觀測到的。故可以作為隱馬爾可夫模型的觀測值，通過分析QAR數據，也就是分析觀測值，訓練得出隱馬爾可夫模型參數，即隱馬爾可夫模型的學習問題。狀態轉換如圖1所示。

圖1 狀態轉換圖

訓練得到空中顛簸故障的隱馬爾可夫模型參數后，給定一個新的QAR時間序列，即可推測序列中數據所對應的最優的空中顛簸狀態趨勢，進行狀態趨勢分析，即隱馬爾可夫模型的解碼問題。但是由上面對QAR數據特點可知，QAR數據是高維數據，數據量大，且不同屬性的數據類型也不相同，同一屬性的觀測值繁多。而隱馬爾可夫模型的觀測值的狀態是有限的，故現在面臨的問題是如何對QAR進行處理，使其符合隱馬爾可夫模型中對觀測值的要求。

2 基于聚類的HMM模型介紹

通過以上分析可知，由于QAR數據具有高維度、數據量龐大和數據復雜等特點，想要建立空中顛簸故障的隱馬爾可夫模型，需對QAR數據進行處理，使其離散化為有限個狀態值，從而符合隱馬爾可夫模型對觀測值的要求。本文選擇用聚類算法對其QAR數據進行處理，但是直接對原始數據聚類，效率太低，在聚類之前先對QAR數據進行分段線性表示。然后利用改進的廣度優先鄰居搜索聚類算法對QAR數據進行聚類，根據聚類結果對數據進行離散化。最后利用訓練數據得到隱馬爾科夫模型參數。

2.1 基于固定分段數的自底向上分段算法

分段線性表示有數據壓縮和過濾的作用。時間序列的線性分段表示方法很多，各類時間序列的特征不同，因此，分段表示及分析的方法也就不同。其次，我們對各類時間序列分析的目的也不相同，針對QAR數據的分析，因為正常情況的時間序列是平穩波動、有規律的，但發生異常或故障時就會出現突然的、零星的尖峰，這些異常則是我們關注的重點。針對眾多時間序列的分段線性算法，選擇合適QAR數據的分段算法，才能更好表示數據。選擇合適的時間序列分段表示方法，不僅可以壓縮飛行參數的數據量，同時還可以有效剔除原始參數序列中的噪聲干擾。

廖俊等[7]提出的基于趨勢轉折點的分段線性表示，是只通過計算某個數據與前后兩個數據的差值或距離來判斷趨勢點或特征點，都是局部算法，其精度和壓縮比都不能很好地平衡。

如圖2所示，在第400數據點附近，當出現一段數據連續上升或連續下降，且每次上升和下降的幅度較小時，存在的轉折點或特征點無法識別，故不適合QAR數據。

圖2 基于特征點的分段線性表示

自底向上分段算法，設置擬合誤差閾值α，首先將將N個待分段的時間序列數據點兩兩連接 ,劃分成不重合的N/2個初始分段即劃分為相鄰點的短序列，此時的擬合誤差為0。然后將相鄰兩個短序列連接起來成為一個較長序列 ,再次計算中間所有點的擬合誤差，找到誤差最小且小于閾值α的分段即為第二個分段。依此方式循環，從中選擇擬合代價最小的 ,如果該最小值小于分段閾值α, 則合并對應的兩個相鄰段, 直到所有分段的擬合誤差都小于閾值α，分段完成。

但是自底向上算法進行分段時，設定誤差閾值，雖然可以保證分段的精確度，但是所得的分段數無法確定，不利于進一步的分析。故本文借鑒固定分段數的思想，提出基于固定分段數的自底向上(BU)線性分段算法。

算法的主要思想:設置最大分段數N，擬合誤差閾值e，進行自底向上的線性分段時，依次遍歷，分別計算兩點中間的所有點與該兩點之間的最小二乘擬合線段的誤差平方和，找到誤差平方和最小且小于誤差閾值e的線段，即為新的一個分段，用M記錄已經得到的分段數，當M大于N時，停止分段。這樣既保證了分段精度，也避免了由于精度值設置太大出現太多分段的現象，而且還可以根據實際情況的需要來確定分段數目。

擬合誤差的平方和公式為：

Err=∑(y-y′)2

(1)

式中：y和y′分別表示兩個分段點之間在同一時間點tc上的原始時間序列值和擬合值。其中擬合值是在兩個分段點連接起來的線段上tc時刻對應的值，計算公式如下(其中k代表線段的斜率)：

k=(y1-y2)/(t1-t2)

(2)

y′=y1-k×t1+k×tc

(3)

在QAR數據中，參數眾多，下面主要以06AccelVert屬性的數據作為實驗數據，具體說明實驗結果。對數據進行預處理，提取所需要的巡航階段QAR數據。采用上文提到的基于固定分段數的自底向上線性分段算法對1 000條飛機巡航階段數據進行線性分段。規定分段的擬合誤差閾值為0.5，分段數為20，分段數達到20時，停止分段。

圖3中x軸表示時間，y軸表示當前時間的飛機06AccelVert屬性的取值。

圖3 基于固定分段數的BU分段結果圖(長度1 000)

時間序列分段后，需要選擇合適的模式來表示數據，不同的表示方式對后面的聚類也會有所影響。在大多數的線性分段算法方法中，常用線段平均值描述線段，也有其他的描述方法。章登義等[9]用斜率和長度描述線段?；谥匾c的飛機發動機數據異常子序列檢測[9]中則是用中值描述線段。針對QAR數據，分別用對兩種表示方法進行分析。這兩種方法描述QAR數據都有一定的局限。

第一，QAR數據正常時波動不大，波動較大時也不會超過0.4，因此會出現很多斜率很接近，但長度相差較大的線段所表示的狀態是相同的。而針對線段斜率和長度進行聚類時，由于長度差異較大，這些狀態相同的線段無法聚到一個類中。

第二，QAR數據異常時，數據波動較大，線段端點值較正常值差別很大，當用線段中點值來表示數據時，端點值被抵消，嚴重降低了異常數據與正常數據的差別。

根據以上的分析，使用線段的斜率和長度，或者斜率與中值都無法很好地表示線段。故本文針對QAR數據的特點對線段的表示方法進行改進，采用線段的斜率和線段端點值差的絕對值來表示。

圖4中x軸是分段后所的線段的序列，一共20條線段，y軸是線段的值，其中虛線表示線段的斜率，實線表示線段的兩端點的差的絕對值。可以看出，斜率變化較大的線段兩端差的絕對值也較大，斜率變化小的線段兩端差的絕對值小，該方法能很好地表示數據，更利于線段的聚類。

圖4 所得分段的斜率和端點差的絕對值

2.2 二次廣度優先鄰居搜索聚類

廣度優先鄰居搜索聚類[10]即廣度優先搜索某對象的直接鄰居和間接鄰居,將找到符合條件的所有鄰居合并聚成一類。接著重復該步驟完成所有對象的聚類。

算法的基本思路是從待聚類的對象中的任意對象開始,基于廣度優先和距離參數r,依次搜索該對象的直接鄰居和間接鄰居,對找到的每一個對象進行評估。如果符合設定的參數θ標準 ,則將它們聚為一類。接著重復該步驟直到所有對象都已聚類完畢。

廣度優先鄰居搜索聚類存在一定的局限性，因為采用相似矩陣表示對象之間的相似性，當數據量太大時，算法效率降低。QAR數據是時序數據，數據量巨大，針對這一局限性，對廣度優先鄰居搜索聚類進行改進以更好地適用于QAR數據。

QAR數據為流數據，即時序數據。關于時序數據的聚類，倪巍偉等[11]曾提出一種基于分區的思想，即首先按照時間順序對數據進行分區，然后對分區中的數據進行聚類。本文對廣度優先鄰居搜索聚類的算法進行改進時，借鑒了分區的思想。即將QAR數據劃分為數據塊，每個數據塊包含相同數目的數據點，分別對每個區塊內的QAR數據進行廣度優先鄰居搜索聚類。得到每個區塊內的聚類結果后，計算聚類中心，再對所得的聚類結果進行廣度優先鄰居搜索聚類，合并相似度高的類，得到最終的聚類結果。即二次廣度優先鄰居搜索聚類。

定義1設QAR時序數據X(X1,X2,…,XN)長度為N，數據劃分長度為S，則數據被劃分成M塊。M=N/S。

用上文中提出的基于固定分段數的自底向上分段算法對每個區塊內的數據進行分段。

定義2設某區塊的分段數為c，用線段斜率和端點差的絕對值表示數據，數據表示為((k1,d1),(k2,d2),…,(kc,dc))。

飛機發生故障時表現為QAR數據波動異常，異常點往往是孤立點或者離群點，與同類算法相比, 改進的廣度優先搜索聚類算法具有實現簡單、復雜度低和容易設定最佳參數等優點, 并且能夠很好地識別孤立點,發現任意形狀的聚類。就QAR數據特點而言，該算法合適對QAR數據進行聚類。

QAR數據是一種時序序列，針對時間序列的數據挖掘，通常要將數據點密集的長時間的序列轉換成相對稀疏的具有特定意義的符號序列，即轉化為相對簡單且易于分析處理的形式[12]。采用二次廣度優先鄰居搜索聚類算法對序列進行聚類，分析聚類結果，給每類賦予相應的類別標簽，并用類標簽代替原始數值序列，實現數據離散化。為隱馬爾科夫模型的建立做準備。

2.3 模型建立

QAR數據的隱馬爾可夫模型的模型建立大致分為以下三個階段：第一階段，線性分段，對QAR數據進行壓縮；第二階段：聚類，得到聚類結果，分析聚類結果，實現數據符號化，得到觀測狀態；第三階段：使用前向后向算法訓練數據建立HMM模型。

具體實現步驟如下：

步驟1進行QAR數據預處理，得到最初的原始數據，并將其分為訓練數據和測試數據。

步驟2對所有數據用基于固定分段數的自底向上分段算法進行線性分段表示。

步驟3對分段后所得線段用二次廣度優先鄰居聚類算法進行聚類，得到結果，分析聚類結果，實現數據符號化，得到觀測狀態，將數據序列都轉化為有限個觀測態的序列。

步驟4用MATLAB中estimate方法初始化模型參數λ0。

步驟5訓練數據的序列作為離散隱馬爾科夫的觀測序列，用EM算法訓練隱馬爾科夫模型參數，即隱馬爾可夫模型的學習問題，獲得最優模型λ。

3 實 驗

實驗工具選用MATLAB R2010a，實驗平臺配置為2.13 GHz的CPU、6 GB內存的PC，操作系統為Windows 7，實驗數據某B737型飛機記錄的數據。

3.1 QAR數據預處理

本文的QAR數據，來自于已經過QAR譯碼的CSV文件。CSV文件的QAR數據是已經過譯碼后的工程值文件。 通過對CSV文件的研究可知，飛機發生的許多故障通常和一部分QAR數據屬性相關，因此可以從發生的故障類型和發生故障的機型著手，去查看QAR數據CSV文件中的相關屬性列，并選取某些屬性和某些故障作針對性的研究?？罩蓄嶔す收舷嚓P的屬性主要有EPR(發動機壓比)、EGT(排氣溫度)、ALV(垂直加速度)、VIB(振動)、N1(低壓轉子轉速)和N2(高壓轉子轉速)。按照層次分析法，專家對空中顛簸故障發生相關屬性重要程度進行了統計和調查。根據專家經驗可得飛機在空中是否有空中顛簸的故障發生主要與飛機的垂直加速度的屬性值有關，故本文只提取垂直加速度ALV數據來分析空中顛簸故障。數據預處理后得到飛機巡航階段的垂直加速度(ALV)數據，用其中 57 576條數據作為訓練數據。

3.2 QAR數據分析

運用改進的廣度優先鄰居搜索聚類算法對數據聚類。首先對訓練數據進行分區，這里取每個區長度分別為S=1 000，對每個數據區進行基于固定分段數的自底向上線性分段，并用每個線段的斜率和端點值差的絕對值表示線段，即此時每條線段可用二維上的點來表示。下面以某個區為代表進行說明，圖5為57區數據，圖6為第57區的聚類結果。

圖5 第57區的線段

圖6 第57區線段聚類結果

每個分區聚類結束后都要計算并記錄下每一類的聚類中心，聚類中心即為每類中線段的平均值，存儲在新的數據文件中。用分區聚類結束后，所有類的聚類中心在同一個數據文件中。針對這些聚類中心進行二次聚類，即并合并相似度高的類，聚類結果如圖7所示。

圖7 對區塊類的聚類中心的聚類結果

對聚類結果進行分析，圖7中斜率接近0的線段，線段的端點值差的絕對值也小，這樣的線段被聚為一類，這類數據幾乎沒有變化，即數據趨勢為平穩；圖7中斜率小于0，且端點差的絕對值在0.2附近的線段被聚為一類，這些數據趨勢為輕微下降；同樣圖7中斜率大于0，且端點差的絕對值在0.2附近的哪一類線段的數據趨勢為輕微上升；斜率小于0，且端點差絕對值變化較大的數據，其數據趨勢為突然下降；同理，斜率大于0，且端點差絕對值變化較大的數據，其數據趨勢為突然上升。

對所得聚類結果進行統計，根據以上分析、匯總可得到訓練數據的趨勢情況，匯成表格結果見表1。

表1 線段趨勢情況

通過以上對聚類結果的分析與統計，采用效果最合適的聚類對其數據進行符號化，用5種趨勢的標簽代替原始數據的值作為離散隱馬爾科夫的觀測序列。故此時隱馬爾科夫的觀測狀態共有5種，即平穩、輕微上升、輕微下降、突然上升和突然下降。分別用0、1、-1、2、-2表示。圖8為訓練數據的線段的趨勢。

圖8 線段趨勢圖

由分析可知，飛機發生空中顛簸故障時，主要影響到的是飛機垂直加速度(ALV)的屬性值。故將飛機的運行狀態作為隱馬爾科夫的隱狀態，即設定隱馬爾可夫模型的隱狀態有三種，分別是正常、輕微顛簸和嚴重顛簸，用數值1、2、3表示。根據聚類中心值可初步確定觀測序列所對應的初始的狀態序列。在已知觀測序列和預設狀態序列的情況下，可用MATLAB中的hmmestimate方法估計狀態轉移概率矩陣和發射概率矩陣。將此作為EM算法的初始值。

根據得出的觀測序列及與其對應的狀態序列估計其狀態轉移概率矩陣和發射概率矩陣，并采用前向后向算法得出訓練后的狀態轉移概率矩陣和發射概率矩陣。分別如下：

訓練前：

訓練后：

用訓練所得HMM模型預測測試數據中飛機的運行狀態。數據同樣來自某 B737型飛機記錄的并經過數據預處理后得到的飛機巡航階段的垂直加速度(ALV)數據，從中提取四組數據，分別是2 000條、2 000條、5 000條、5 000條。首先對四組測試數據結合之前的聚類方法，運用k-近鄰算法得出其分類標簽，完成數據的符號化，即用線段的趨勢來表示符號化數據。然后采用Viterbi算法預測其狀態。

下面為測試數據中不同的兩段數據的預測情況對比。圖9為數據1中某段原始數據，圖10為該段數據的預測結果。

圖9 數據1中某段原始數據

圖10 數據1中某段預測結果

圖9為原始序列圖，橫坐標為數據序列號，縱坐標為垂直加速度屬性值。圖10是預測結果圖，其中橫坐標還是數據序列號，縱坐標為飛機運行狀態。從中可以看出，在飛機出現空中顛簸時，垂直加速度數據異常，該模型能夠預測出異常數據，并判斷其顛簸的程度。圖10識別出了兩段嚴重顛簸狀態。

四組數據的具體預測結果如表2所示。

表2 預測情況統計

通過以上大量的測試實驗，該模型能很好地預測原數據的狀態趨勢。但是當需要識別的數據長度增加時，其預測精度有所降低，但大部分的狀態能夠預測。可知通過對QAR的數據進行趨勢分析，以平穩、輕微上升、輕微下降、突然上升和突然下降這5種趨勢來表示QAR數據，將QAR數據轉化為有限個觀測值，作為隱馬爾可夫模型的觀測序列。經過學習訓練建立的隱馬爾可夫模型能很好地預測空中顛簸數據的狀態趨勢。雖然該模型目前能夠預測出顛簸狀態，但是飛機出現空中顛簸是不確定性的，因此需要大量的監控歷史數據來不斷訓練學習，才能保證更加精確地預測飛機的運行狀態。

4 結 語

本文針對傳統的QAR數據分析方法只關注異常點而忽略了QAR的狀態趨勢特征的問題。通過分析QAR數據與飛機空中顛簸故障，針對大量QAR數據，提出一種基于聚類的HMM模型，分析發生故障或異常時QAR數據中不同屬性的變化特點。通過對聚類進行數據離散化得出數據的狀態趨勢，并用類標簽代替原始數值序列，進行數據離散化，對飛機巡航階段空顛簸故障數據進行了有效處理，建立了有效的HMM模型。用HMM來統計性地描述故障發生時QAR數據的改變,以狀態序列的形式描述故障的演化過程，并能夠根據觀測數據有效識別飛機所處運行狀態趨勢，成功驗證了該方法的有效性。

[1] 曹惠玲,周百政.QAR數據在航空發動機監控中的應用研究[J].中國民航大學學報,2010,28(3):15-19.

[2] Jiang H,Li X,Liu C.Large margin hidden Markov models for speech recognition[J].IEEE Transactions on Audio Speech & Language Processing,2005,14(5):1584-1595.

[3] 黃曉彬,王春峰,房振明,等.基于隱馬爾科夫模型的中國股票信息探測[J].系統工程理論與實踐,2012,32(4):713-720.

[4] 張璇,周峰.隱馬爾科夫模型應用領域、熱點及趨勢分析——基于CiteSpaceII[J].現代商貿工業,2015,36(15):63-65.

[5] 于天劍,陳特放,陳雅婷,等.HMM在電機軸承上的故障診斷[J].哈爾濱工業大學學報,2016,48(2):184-188.

[6] 柳姣姣,禹素萍,吳波,等.基于隱馬爾科夫模型的時空序列預測方法[J].微型機與應用,2016,35(1):74-76.

[7] 廖俊,于雷,羅寰,等.基于趨勢轉折點的時間序列分段線性表示[J].計算機工程與應用,2010,46(30):50-53.

[8] 章登義,歐陽黜霏,吳文李.針對時間序列多步預測的聚類隱馬爾科夫模型[J].電子學報,2014(12):2359-2364.

[9] 楊慧,王光霞.基于重要點的飛機發動機異常子序列檢測[J].計算機工程與設計,2016,37(9):2543-2547.

[10] 錢江波,董逸生.一種基于廣度優先搜索鄰居的聚類算法[J].東南大學學報(自然科學版),2004,34(1):109-112.

[11] 倪巍偉,陸介平,陳耿,等.基于k均值分區的數據流離群點檢測算法[J].計算機研究與發展,2006,43(9):1639-1643.

[12] 任江濤,何武,印鑒,等.一種時間序列快速分段及符號化方法[J].計算機科學,2005,32(9):166-169.