彭沛，趙永平，*，王雨瑋

1.南京航空航天大學(xué) 能源與動(dòng)力學(xué)院，南京 210016

2.北京航空工程技術(shù)研究中心，北京 100076

在大數(shù)據(jù)時(shí)代，航空發(fā)動(dòng)機(jī)采集的傳感器信息的應(yīng)用由最初的事故原因調(diào)查擴(kuò)展到發(fā)動(dòng)機(jī)設(shè)計(jì)性能評(píng)估、飛行訓(xùn)練考核輔助評(píng)價(jià)、協(xié)助地勤人員進(jìn)行狀態(tài)監(jiān)控及視情維護(hù)等方面。伴隨著飛機(jī)整體性能的提升和發(fā)動(dòng)機(jī)結(jié)構(gòu)的復(fù)雜化，人們對(duì)飛行安全及發(fā)動(dòng)機(jī)可靠性方面將會(huì)提出更高的要求，因而充分挖掘并利用飛行中發(fā)動(dòng)機(jī)采集到的數(shù)據(jù)信息將會(huì)成為后續(xù)研究的必然趨勢(shì)。傳感器的測(cè)量過程由于發(fā)動(dòng)機(jī)工作的復(fù)雜環(huán)境而易受隨機(jī)噪聲的影響，這就導(dǎo)致地面研究人員對(duì)原始獲取的數(shù)據(jù)難以直觀地分析研究。因此，如何將外場(chǎng)獲取的數(shù)據(jù)最大程度還原出發(fā)動(dòng)機(jī)在飛機(jī)飛行時(shí)真實(shí)的工作狀態(tài)，這對(duì)后續(xù)基線模型的建立、實(shí)施狀態(tài)監(jiān)控及維護(hù)決策有重要的意義［1-4］。但研究人員處理發(fā)動(dòng)機(jī)獲取的傳感器數(shù)據(jù)目前仍采用人工判讀的方法，尤其在批量分析時(shí)會(huì)因判讀數(shù)據(jù)量大而耗時(shí)耗力。另外，發(fā)動(dòng)機(jī)環(huán)境參數(shù)的不同和個(gè)體制造的差異會(huì)導(dǎo)致判讀規(guī)律的變化，這也增加了研究人員對(duì)工作模式判讀的難度。

國(guó)內(nèi)外為解決上述問題做出很多嘗試。王奕惟等［5］將快速存儲(chǔ)記錄器（Quick Access Re‐coder，QAR）的數(shù)據(jù)圖像化并采用卷積神經(jīng)網(wǎng)絡(luò)來判別狀態(tài)。該方法本質(zhì)上是對(duì)時(shí)間序列人為采樣的單點(diǎn)識(shí)別而未考慮到傳感器數(shù)據(jù)的動(dòng)力學(xué)性質(zhì)，這割裂了時(shí)序上的連續(xù)性，并且在識(shí)別過程中該方法也未考慮過渡態(tài)。CWT-CNN 方法［6］結(jié)合了小波變換和卷積神經(jīng)網(wǎng)絡(luò)，但其識(shí)別的單維時(shí)序樣本并非含多種模式的信號(hào)，而是人為分割后單一模式的信號(hào)。此外，Cynthia 等［7］提出了將剪枝精確線性時(shí)間法［8］和自組織神經(jīng)網(wǎng)絡(luò)結(jié)合并給出過渡態(tài)的判別方法。該模型成功地識(shí)別了飛機(jī)的多個(gè)飛行狀態(tài)，但僅涉及單維時(shí)序信號(hào)如油門桿角度，這樣就忽視了真實(shí)數(shù)據(jù)中的時(shí)滯效應(yīng)，即按油門桿角度切換模式時(shí)，其他監(jiān)視參數(shù)如轉(zhuǎn)速、瞬時(shí)耗量等不會(huì)同步地發(fā)生改變。

多維時(shí)序數(shù)據(jù)的挖掘在機(jī)械零件的傳感器分析、金融財(cái)政等早有廣泛的應(yīng)用。此類問題可細(xì)分為時(shí)序分割、模式發(fā)現(xiàn)、模式匹配等任務(wù)。例如基于動(dòng)態(tài)規(guī)劃的算法［9］、剪枝精確線性時(shí)間法、基于模式的隱馬爾可夫模型［10］、基于高斯假設(shè)的貪婪分割模型［11］。這幾類方法僅能由不同模式來分割序列，并不能直接識(shí)別出模式，而且這些方法對(duì)參數(shù)調(diào)節(jié)很敏感，例如片段數(shù)和相關(guān)閾值等。在航空發(fā)動(dòng)機(jī)收集到的大批量多維的時(shí)序數(shù)據(jù)下，不同航段采集的時(shí)序數(shù)據(jù)包含的工作模式數(shù)量和類別均不同，不同時(shí)序數(shù)據(jù)的要調(diào)節(jié)的參數(shù)有很大差別且隨著數(shù)據(jù)量的增加調(diào)節(jié)參數(shù)往往會(huì)耗費(fèi)大量的時(shí)間。在這種場(chǎng)景下理想的方法應(yīng)該能尋找任意的模式并不受參數(shù)的影響，即無參數(shù)化。AutoPlait［12］是基于多級(jí)鏈隱馬爾可夫模型，它利用無損壓縮原則實(shí)現(xiàn)模型無參數(shù)化。Auto-Cyclone［13］也采用相似的思路實(shí)現(xiàn)模式季節(jié)性挖掘。但上述2 種方法并不適合于發(fā)動(dòng)機(jī)的工作狀態(tài)的挖掘，它們更傾向于挖掘有季節(jié)性的模式而模糊了發(fā)動(dòng)機(jī)的穩(wěn)態(tài)模式和過渡態(tài)模式的界限，并且復(fù)雜的過渡態(tài)模式也無法給出很好的處理方法。事實(shí)上，上述2 種無參數(shù)方法更適合于飛行動(dòng)作識(shí)別，但模型訓(xùn)練也需要豐富的飛行姿態(tài)數(shù)據(jù)。

考慮到現(xiàn)有方法在實(shí)際工程中遇到的諸多問題，本文針對(duì)航空發(fā)動(dòng)機(jī)工作時(shí)產(chǎn)生的時(shí)序數(shù)據(jù)，設(shè)計(jì)了一個(gè)快速自動(dòng)挖掘工作模式的方法，即AutoMiner。該方法自然地集成了時(shí)序分割、模式發(fā)現(xiàn)和模式識(shí)別任務(wù)，它以最原始的全航段多維時(shí)序數(shù)據(jù)為研究對(duì)象，采用最小編碼代價(jià)的原則進(jìn)行模型無參數(shù)化，以應(yīng)對(duì)多變的時(shí)序數(shù)據(jù)。并且該方法還借鑒了標(biāo)記傳播的思想［14］來解決過渡態(tài)的識(shí)別問題。在模型訓(xùn)練方面，該方法支持并行化并利用儲(chǔ)存單元省去多次重復(fù)計(jì)算，這大大縮短了訓(xùn)練時(shí)間。本文模型的設(shè)計(jì)初衷為航空發(fā)動(dòng)機(jī)工作模式的挖掘，但由于模型的可遷移性，還能被應(yīng)用于更高級(jí)模式的識(shí)別，如飛行動(dòng)作。最后，在準(zhǔn)確度、迭代效率、可遷移性等方面，本文開展一系列實(shí)驗(yàn)，通過與多種方法對(duì)比可視化地展示了AutoMiner 方法的優(yōu)越性。

1 AutoMiner 方法優(yōu)化目標(biāo)

該節(jié)介紹工作模式的線性模型，由最小化編碼損失的方式提出優(yōu)化目標(biāo)，并給出利用標(biāo)記傳播法識(shí)別過渡態(tài)在內(nèi)的發(fā)動(dòng)機(jī)工作模式的方法。

1.1 工作模式的線性模型

從航空發(fā)動(dòng)機(jī)上獲取指定的傳感器數(shù)據(jù)，可組成多維時(shí)間序列X=[x1，x2，…，xn]，xi為第i時(shí)刻的d維列向量。假設(shè)模型根據(jù)X的本征模式能將其分割成K段，即S={s1，s2，…，sK}，第i段si可視為X中由斷點(diǎn)τi和τi+1分割的部分，則斷點(diǎn)序號(hào)集組成τ={τ1，τ2，…，τK+1}，其中τ1=1，τK+1=n。S中每段按照工作模式可以賦予相應(yīng)標(biāo)記并組成F={f1，f2，…，fK}，fi∈{1，2，…，c}表示第i段si屬于工作模式fi，c為工作模式總數(shù)。

當(dāng)航空發(fā)動(dòng)機(jī)處于相同工作模式時(shí)，內(nèi)在物理規(guī)律應(yīng)是相似的，工作模式的轉(zhuǎn)變會(huì)引起這種內(nèi)在物理規(guī)律的變化，而這種規(guī)律體現(xiàn)到可測(cè)傳感器數(shù)據(jù)上就是多維時(shí)序數(shù)據(jù)的變化。由此思路，對(duì)S中每個(gè)片段建立數(shù)學(xué)模型，自然地認(rèn)為相同工作模式的數(shù)學(xué)模型是相似的，而不同工作模式的模型會(huì)有顯著區(qū)別。對(duì)第i段si可建模為

1.2 總編碼代價(jià)的計(jì)算

在計(jì)算過程中，編碼描述數(shù)據(jù)X及儲(chǔ)存算法模型需要付出一定的代價(jià)，而最小編碼代價(jià)對(duì)應(yīng)的模型是最理想的，這個(gè)思想與經(jīng)典的奧卡姆剃刀定律相似。于是，可建立優(yōu)化目標(biāo)為

式中：CT表示編碼的總代價(jià)，即本文所需最小化的優(yōu)化目標(biāo)，分別由描述算法模型的編碼代價(jià)和算法模型應(yīng)用到數(shù)據(jù)X后產(chǎn)生的編碼代價(jià)組成。

假設(shè)本文討論的編碼過程均是無損的，定義符號(hào)log*(a) 為對(duì)整數(shù)a的編碼長(zhǎng)度［15］，即log*(a)=log2(a)+log2(log2(a))+…，按此規(guī)律求和且只有正數(shù)項(xiàng)包含在求和項(xiàng)中。AutoMiner算法模型簡(jiǎn)式為：{Θ1，Θ2，…，ΘK，S，F(xiàn)，K，r}，r代表S中片段的總類別數(shù)，區(qū)別于1.1 節(jié)中提及的工作模式數(shù)c。例如本文涉及到8 個(gè)穩(wěn)態(tài)模式，則包含過渡態(tài)后總的模式數(shù)c為9，在優(yōu)化目標(biāo)中總類別數(shù)r的計(jì)算方法為

式中：IF(?)為判斷函數(shù)，例如對(duì)IF(l∈F)而言，若語句l∈F為真，則取值1，反之為0。式（3）中r將前8 個(gè)穩(wěn)態(tài)按是否出現(xiàn)計(jì)算，第9 個(gè)表示的過渡態(tài)按出現(xiàn)次數(shù)計(jì)算。

式中：cF為浮點(diǎn)系數(shù)［16］，其默認(rèn)取值為4×8，即32 bit；q為線性模型的階數(shù)。

由①～⑤便可獲得算法模型編碼所需的代價(jià)，即式（2）第1 項(xiàng)，而式（2）求和中的第2 項(xiàng)，體現(xiàn)模型在給定數(shù)據(jù)X完成任務(wù)的能力，即序列分割及模式挖掘的準(zhǔn)確度。本文采用哈夫曼編碼，將編碼代價(jià)轉(zhuǎn)化為線性模型的似然概率，得到

式中：Θ={Θ1，Θ2，…，ΘK}；P(X|Θ) 為似然函數(shù)值。

考慮到同一工作模式中線性模型的相關(guān)參數(shù)應(yīng)盡可能相似，利用式（1）中偏差項(xiàng)服從正態(tài)分布的性質(zhì)，構(gòu)造對(duì)數(shù)似然函數(shù)為

式中：det 表示對(duì)矩陣取行列式；l(i)表示由第i個(gè)片段計(jì)算得到的對(duì)數(shù)似然函數(shù)。

在序列分割過程中，當(dāng)片段長(zhǎng)度|si|小于維數(shù)d時(shí)會(huì)導(dǎo)致矩陣求逆出現(xiàn)病態(tài)，本文引入正則化的方法將式（6）進(jìn)行改造，如式（7）所示，這種正則化方法常用于解決高維時(shí)間序列的問題［17］。

式中：λ為正則化系數(shù)，常取較小值，本文取10?5；tr 表示矩陣的求跡運(yùn)算。

通過極大化對(duì)數(shù)似然函數(shù)能求解出參數(shù)Θ，令式（7）對(duì)Θ的偏導(dǎo)數(shù)分別為0，得

式中：Ιd表示d階單位矩陣。當(dāng)不考慮正則化時(shí)即λ=0，Σi與相等。式（7）求偏導(dǎo)數(shù)涉及到諸多矩陣運(yùn)算，對(duì)式（8）更詳細(xì)的推導(dǎo)可見附錄A。

于是，將式（8）結(jié)果重新代入式（7），可得到最大似然值的結(jié)果，并結(jié)合式（5）和式（6）可得

綜合式（5）、式（9）可得式（2）第2 項(xiàng)的表達(dá)式為

式中：為便于后續(xù)對(duì)編碼代價(jià)迭代求解斷點(diǎn)，將結(jié)果化為?(τi，τi+1)的求和，且該項(xiàng)只與斷點(diǎn)有關(guān)。

最后可以通過①～④及式（4）、式（10）獲得編碼的總代價(jià)CT，即模型的優(yōu)化目標(biāo)

式中：優(yōu)化目標(biāo)僅為斷點(diǎn)集τ的函數(shù)，隨后便可尋優(yōu)出合適的斷點(diǎn)組合，并使總編碼代價(jià)最小。

1.3 基于標(biāo)記傳播的含過渡態(tài)時(shí)序挖掘方法

在計(jì)算式（11）中編碼的總代價(jià)CT前，還需在尋優(yōu)迭代過程中對(duì)片段集S賦予合適的標(biāo)記集F，即找到每段隸屬的工作模式。在計(jì)算出F后便可由式（3）計(jì)算出總類別數(shù)r。

根據(jù)本文含過渡態(tài)模式的工程場(chǎng)景可做出如下描述：在F標(biāo)記集中，第i段si可以劃分為模式fi∈{1，2，…，9}，其中前8 個(gè)模式為穩(wěn)定態(tài)，第9 個(gè)模式為過渡態(tài)。假設(shè)已獲得待識(shí)別片段序列為s0，對(duì)該段序列建立線性模型可獲得參數(shù)Θ0并組成樣本點(diǎn)。參照行業(yè)工程手冊(cè)以及專家建議可按同樣方式選取特征生成帶標(biāo)記的參照樣本集表示采集第i個(gè)參照片段數(shù)據(jù)后維的特征向量，實(shí)驗(yàn)中為Θ0包含的參數(shù)總數(shù)。記參照樣本集的對(duì)應(yīng)標(biāo)記矩陣為YL=[y1，y2，…，ym]，9 維的標(biāo)記向量yi中除第l個(gè)元素yil為1，其他項(xiàng)均為0，即表示參照數(shù)據(jù)屬于第l類模式。由于可獲取的參照樣本僅為前8類的穩(wěn)態(tài)模式，則參考樣本的標(biāo)記中yi9均為0。

式中：迭代過程中的初始標(biāo)記值Y(0)被賦值為先驗(yàn)標(biāo)記矩陣Y，其值被預(yù)先人為設(shè)置。式（13）中第2 項(xiàng)表示先驗(yàn)信息在傳播時(shí)的參與程度。為讓傳播過程中僅改變標(biāo)記y0，比例矩陣Dα設(shè)置為分塊矩陣，α越大即待標(biāo)記數(shù)據(jù)會(huì)接收更多其他點(diǎn)的信息，本文默認(rèn)取值為0.99［18］。

不難證明式（13）的迭代必收斂于式（14）。

由式（14）獲得軟標(biāo)記y0，其每個(gè)元素表示樣本隸屬對(duì)應(yīng)模式的概率，選取最大概率對(duì)應(yīng)的序號(hào)l即為所屬的模式。本質(zhì)上，參照樣本集雖只含穩(wěn)態(tài)模式，但該方法利用標(biāo)記傳播更新生成了軟標(biāo)記，使得過渡態(tài)模式對(duì)應(yīng)的軟標(biāo)記值在這里充當(dāng)著剩余概率的角色，即樣本不屬于前8 種穩(wěn)態(tài)模式的概率。

時(shí)序模式識(shí)別步驟如算法1 所示。為盡可能降低計(jì)算復(fù)雜度，在采用式（13）計(jì)算W時(shí)對(duì)已知的距離矩陣無需重復(fù)計(jì)算。于是由算法1 步驟2 可將算法1 的計(jì)算復(fù)雜度由O(m2)降至O(m)。

算法1 工作模式時(shí)序挖掘算法Algorithm 1 Working mode time-series mining algorithm

2 斷點(diǎn)尋優(yōu)的啟發(fā)式算法

式（12）得到總編碼代價(jià)僅與斷點(diǎn)τ有關(guān)，問題轉(zhuǎn)化為尋找最優(yōu)的斷點(diǎn)組合使總編碼代價(jià)最小，在該斷點(diǎn)組合下對(duì)應(yīng)的片段集S以及標(biāo)記集F即是AutoMiner 方法挖掘時(shí)序數(shù)據(jù)X獲得的最優(yōu)方案。在AutoMiner 方法中提出的斷點(diǎn)優(yōu)化算法是種啟發(fā)式算法，它無需全局遍歷所有的解而是在迭代范圍內(nèi)找尋局部最優(yōu)的組合。每一次迭代時(shí)添加一個(gè)斷點(diǎn)，并局部調(diào)整斷點(diǎn)組合以保證總編碼代價(jià)總是維持著最小值的狀態(tài)，直到找到總編碼代價(jià)不再減小的斷點(diǎn)組合。這個(gè)思路類似標(biāo)準(zhǔn)的“自上而下”方法［19］，但該方法容易陷入局部最優(yōu)并且承擔(dān)著較高的計(jì)算成本。本文將提出的啟發(fā)式算法細(xì)分為圖1 所示的3 個(gè)步驟來分別介紹。

圖1 斷點(diǎn)尋優(yōu)示意圖Fig.1 Schematic diagram of breakpoint optimization

2.1 時(shí)序片段內(nèi)最優(yōu)斷點(diǎn)選擇策略

啟發(fā)式方法如圖1（a）所示。從式（12）中可以看出，每次迭代添加的斷點(diǎn)僅對(duì)斷點(diǎn)所在的片段產(chǎn)生影響。于是做出下述約定：對(duì)于片段si而言，在其內(nèi)添加的一個(gè)斷點(diǎn)t并將片段si分別分割為左右2 個(gè)子片段，其總編碼代價(jià)的變化量為

式中：g是K0、r0、K1和r1的函數(shù)，它由式（4）和式（11）計(jì)算獲得。式（10）中展示了與斷點(diǎn)有關(guān)的函數(shù)?的表達(dá)式。K0和r0為添加斷點(diǎn)t前的總斷點(diǎn)數(shù)和總類別數(shù)，同理添加斷點(diǎn)后可表示為K1和r1。易知K1=K0+1，而r0和r1可由式（3）獲得。

在si內(nèi)尋找最優(yōu)斷點(diǎn)tc的思路如算法2 所示，其中算法中總斷點(diǎn)數(shù)K、添加斷點(diǎn)前的標(biāo)記集F等均默認(rèn)已知。算法2 步驟1 表示從τi+T到τi+1?T以尋優(yōu)步長(zhǎng)T來迭代選取t，設(shè)置較低步長(zhǎng)有利于精確的搜索但會(huì)增加計(jì)算時(shí)間。

此外，算法2 還采取了其他減少計(jì)算開銷的方法。由于對(duì)于算法2 步驟1～步驟7 的循環(huán)之間并未存在聯(lián)系，即可支持并行化計(jì)算來節(jié)省大量計(jì)算時(shí)間。引入儲(chǔ)存單元Pi是為了將此次尋優(yōu)結(jié)果保存下來，防止后續(xù)迭代過程循環(huán)調(diào)用算法2。而選擇總代價(jià)的變化量作為判斷指標(biāo)，也從一定程度上，避免了直接比較總編碼代價(jià)帶來的計(jì)算高額量。若考慮并行化，算法2 對(duì)于si的計(jì)算復(fù)雜度可由降為O(md2)。

算法2 時(shí)序片段內(nèi)部最優(yōu)斷點(diǎn)選擇策略Algorithm 2 Internal optimal breakpoint selection strategy for time-series segments

2.2 斷點(diǎn)的動(dòng)態(tài)更新策略

由于尋優(yōu)算法的策略是啟發(fā)式的，為了最大程度避免陷入局部最優(yōu)解，在每次迭代添加斷點(diǎn)后需要?jiǎng)討B(tài)回顧全局，再次更新斷點(diǎn)組合以保證求解的斷點(diǎn)組合的總編碼代價(jià)保持最低值。參考圖1（b），更新策略如算法3 所示。當(dāng)K不大時(shí)更新次數(shù)將是適度的，而隨著K的增大，每次更新只影響局部斷點(diǎn)且多次迭代后斷點(diǎn)集趨于最優(yōu)組合，所以在容錯(cuò)閾值范圍ξ內(nèi)更新一定會(huì)在適度次數(shù)后收斂，實(shí)驗(yàn)中ξ=5。在算法3 步驟2中，棧儲(chǔ)存添加點(diǎn)或變化點(diǎn)，隨后每次僅對(duì)棧彈出點(diǎn)的鄰域最近兩斷點(diǎn)進(jìn)行更新。在斷點(diǎn)改變后也會(huì)同步更新原來的片段，這時(shí)與片段對(duì)應(yīng)儲(chǔ)存單元中的內(nèi)容將會(huì)清除并等待下次調(diào)用算法2 時(shí)再次賦值。

算法3 添加斷點(diǎn)后的動(dòng)態(tài)更新策略Algorithm 3 Dynamic update strategy after addition of breakpoints

2.3 斷點(diǎn)的尋優(yōu)算法

依照?qǐng)D1（c），算法4 中展示了尋優(yōu)的總體思路。算法4 中步驟3～步驟12 的循環(huán)是獨(dú)立運(yùn)行的并支持并行化運(yùn)算。在算法4 步驟2～步驟27中，若片段si未因斷點(diǎn)更新改變且該時(shí)序片段的最優(yōu)斷點(diǎn)未被此次迭代選中，下一次迭代開始時(shí)便無需對(duì)si重復(fù)尋優(yōu)，這也是設(shè)置存儲(chǔ)單元Pi的目的，因此每個(gè)存儲(chǔ)單元總是唯一地指向片段si。若斷點(diǎn)平均更新次數(shù)為L(zhǎng)，算法3 的計(jì)算復(fù)雜度為，K為求解獲得的最優(yōu)斷點(diǎn)數(shù)?？紤]并行化后，計(jì)算復(fù)雜度為。

算法4 斷點(diǎn)尋優(yōu)的啟發(fā)式算法Algorithm 4 Breakpoint-seeking heuristic algorithm

在某航段上截取部分簡(jiǎn)單的三維時(shí)間序列并可視化地展示AutoMiner 方法斷點(diǎn)尋優(yōu)的工作流程圖，如圖2 所示，其中M2 為過渡態(tài)模式，而M1、M3、M4對(duì)應(yīng)不同的穩(wěn)態(tài)模式。算法初始目標(biāo)是添加單個(gè)合適的斷點(diǎn)使得總代價(jià)下降，這僅需對(duì)時(shí)序數(shù)據(jù)掃描一次便可得圖2（a）。隨后，圖2（b）～圖2（d）為逐步添加斷點(diǎn)的過程。由于算法4步驟22 加入了終止條件，算法將在第6 次迭代時(shí)終止并將第5 次的迭代結(jié)果作為最優(yōu)斷點(diǎn)組合輸出。如圖2（e）所示，若不加以約束添加斷點(diǎn)總代價(jià)反而會(huì)隨之增加。

圖2 AutoMiner 方法尋優(yōu)流程Fig.2 Optimization process of AutoMiner method

3 實(shí)驗(yàn)分析

圖3 是將AutoMiner 方法運(yùn)用于實(shí)際工程場(chǎng)景中的示意圖?？梢钥闯鲈摲椒▽耐鈭?chǎng)采集發(fā)動(dòng)機(jī)上原始的傳感器數(shù)據(jù)進(jìn)行自動(dòng)處理以獲得大批量數(shù)據(jù)能進(jìn)行后續(xù)的分析。參考樣本集按研究的時(shí)序?qū)ο蠖兓？赏ㄟ^專家研究設(shè)計(jì)手冊(cè)或借鑒歷史已標(biāo)注工作模式的數(shù)據(jù)來獲得。此處歷史數(shù)據(jù)中包含了不同退化狀態(tài)下的工作模式數(shù)據(jù)信息，但實(shí)際中退化量與出廠的理論值通常偏差不大，選擇手冊(cè)的信息作為參照也是合理的。

圖3 實(shí)際工程場(chǎng)景應(yīng)用的示意圖Fig.3 Application diagram of engineering scenario

3.1 實(shí)驗(yàn)背景

本文以某型機(jī)載航空發(fā)動(dòng)機(jī)在飛行過程中采集的多組傳感器數(shù)據(jù)為研究對(duì)象，實(shí)驗(yàn)中涉及3 組飛機(jī)全航段工作數(shù)據(jù)，此處一個(gè)航段記錄為飛機(jī)的一次長(zhǎng)時(shí)間飛行任務(wù)。這些數(shù)據(jù)囊括了實(shí)驗(yàn)考慮的全部工作模式，在飛機(jī)飛行時(shí)由4 臺(tái)同型號(hào)發(fā)動(dòng)機(jī)的傳感器進(jìn)行同步測(cè)量。

對(duì)于實(shí)驗(yàn)參數(shù)的篩選需要參考實(shí)驗(yàn)手冊(cè)，選取隨工作模式的變化量最大的多個(gè)物理指標(biāo)，比如絕對(duì)氣壓高度體現(xiàn)了飛機(jī)起飛降落及巡航狀態(tài)，從而也反映發(fā)動(dòng)機(jī)工作模式的變化。在工作模式需要改變時(shí)，飛行員會(huì)通過外部干預(yù)調(diào)節(jié)油門桿角度，這會(huì)引起發(fā)動(dòng)機(jī)轉(zhuǎn)子轉(zhuǎn)速、瞬時(shí)耗量等物理量的變化。為解決參數(shù)的時(shí)滯效應(yīng)，本文選擇上述4 類物理指標(biāo)作為不同工作模式的分辨依據(jù)，其中油門桿角度、轉(zhuǎn)速和瞬時(shí)耗量均由2 組以上傳感器參與測(cè)量，收集以上數(shù)據(jù)可共同組成描述發(fā)動(dòng)機(jī)工作狀態(tài)的多維時(shí)間序列。后續(xù)實(shí)驗(yàn)將證明，考慮發(fā)動(dòng)機(jī)傳感器的多維時(shí)序數(shù)據(jù)比常規(guī)只考慮單組油門桿角度作為辨別依據(jù)的方案更能適應(yīng)復(fù)雜的數(shù)據(jù)情況。本文涉及的發(fā)動(dòng)機(jī)工作模式的人為編號(hào)如表1所示。

表1 某型航空發(fā)動(dòng)機(jī)工作模式說明Table 1 Description of operating mode of aero-engine

3.2 數(shù)據(jù)預(yù)處理

飛機(jī)常需要應(yīng)對(duì)復(fù)雜的飛行條件甚至完成某些高難度的動(dòng)作，這使得發(fā)動(dòng)機(jī)上測(cè)量的參數(shù)不可避免地會(huì)產(chǎn)生波動(dòng)，并且熱力及電磁環(huán)境干擾和不穩(wěn)定的傳感器測(cè)量狀態(tài)也會(huì)引入噪聲。為了減少噪聲以提高挖掘時(shí)序數(shù)據(jù)信息的準(zhǔn)確度，便需要對(duì)采集到的數(shù)據(jù)進(jìn)行平滑降噪。選取某航段上一臺(tái)發(fā)動(dòng)機(jī)的油門桿角度進(jìn)行比較，實(shí)驗(yàn)方法包括滑動(dòng)中值濾波［20］、滑動(dòng)均值濾波以及小波軟閾值濾波［21］。對(duì)比實(shí)驗(yàn)結(jié)果如圖4 所示。數(shù)據(jù)點(diǎn)之間采樣間隔為1 s，考慮到濾波時(shí)盡可能保留過渡態(tài)模式信息，選擇滑動(dòng)窗口的長(zhǎng)度為60 s，同時(shí)采用8 層小波分解來計(jì)算近似系數(shù)。

圖4 3 種濾波方法的對(duì)比結(jié)果Fig.4 Comparative results of three filtering methods

由于中位數(shù)面對(duì)異常數(shù)據(jù)以及高頻噪聲更加魯棒，從圖4 中可以看出相比均值濾波而言，中值濾波能得到更適合識(shí)別工作模式的時(shí)序數(shù)據(jù)。在小波軟閾值濾波中，不同分解的層數(shù)及合適的小波基會(huì)隨著時(shí)序信號(hào)的改變而變化，因此不適用于大批量自動(dòng)處理數(shù)據(jù)的場(chǎng)景。

3.3 具體工程場(chǎng)景下的算法實(shí)驗(yàn)

為展示AutoMiner 方法在航空發(fā)動(dòng)機(jī)工作模式挖掘場(chǎng)景下的優(yōu)越性，本節(jié)提取部分航段中傳感器采集的時(shí)序數(shù)據(jù)與多種時(shí)序分割算法對(duì)比，包括基于高斯核的KCP 方法［22］和基于高斯模型的GSS 方法［11］。為在模式挖掘和無參數(shù)化算法方面進(jìn)一步比較，對(duì)比算法中還包括成功運(yùn)用于飛行模式挖掘的PELT-SOM 算法［7］，以及基于多級(jí)鏈隱馬爾可夫模型的AutoPlait無參數(shù)化算法［9］。

3.3.1 時(shí)序分割能力的對(duì)比評(píng)估

在時(shí)序分割方面，合適的評(píng)價(jià)指標(biāo)至今仍是廣泛討論的內(nèi)容。例如，專家評(píng)判分割斷點(diǎn)在第1 000 s，模型求解出斷點(diǎn)為第1 005 s 是否有效。若容錯(cuò)閾值為10，模型得到的斷點(diǎn)為第1 011 s是否應(yīng)該被懲罰。嚴(yán)格的懲罰制度得到準(zhǔn)確度、混淆矩陣等指標(biāo)并不適合時(shí)序分割問題［23］。上述問題歸結(jié)于對(duì)評(píng)判邊界并未做出合適的處理，本文借鑒文獻(xiàn)［24］的思路選擇評(píng)價(jià)指標(biāo)為R分?jǐn)?shù)，其表達(dá)式為

式中：τ(1)為模型求解的斷點(diǎn)集。數(shù)據(jù)真實(shí)的斷點(diǎn)集τ(2)中距離斷點(diǎn)最近的斷點(diǎn)表示為，此處討論不包含時(shí)序數(shù)據(jù)的起點(diǎn)和終點(diǎn)。指標(biāo)R的取值范圍為0～1，其值越大說明該模型分割效果越好。

采集部分飛行航段上發(fā)動(dòng)機(jī)傳感器數(shù)據(jù)進(jìn)行了2 組實(shí)驗(yàn)，實(shí)驗(yàn)中各物理量?jī)H展示1 個(gè)傳感器的測(cè)量結(jié)果并通過一定比例地?zé)o量綱化操作繪制出如圖5、圖6 所示的曲線，其中曲線圖中的分割線為人工結(jié)合工程經(jīng)驗(yàn)給出，即可視為真實(shí)的分割斷點(diǎn)位置。另外，圖5（b）、圖5（c）、圖5（f）和圖6（b）、圖6（c）、圖6（f）涉及算法均可支持工作模式的識(shí)別，而圖5（d）、圖5（e）和圖6（d）、圖6（e）中僅支持對(duì)時(shí)序數(shù)據(jù)的分割，其中對(duì)比算法中參數(shù)β為懲罰系數(shù)，λ為正則化系數(shù)，M 表示某些對(duì)比算法對(duì)該工程場(chǎng)景失效而僅識(shí)別出的未知模式。

圖5 在航段1 中工作模式挖掘的可視化圖Fig.5 Visualization of pattern mining in Flight segment 1

圖6 在航段2 中工作模式挖掘的可視化圖Fig.6 Visualization of pattern mining in Flight segment 2

將時(shí)序分割結(jié)果代入式（16）中可得到表2。理想的算法應(yīng)該具有較高的R分?jǐn)?shù)且生成斷點(diǎn)數(shù)接近于真實(shí)斷點(diǎn)數(shù)。由于圖5、圖6 中可以看出航段2 相比航段1 包含更豐富的工作模式，在復(fù)雜的數(shù)據(jù)組成表2 中航段2 的R分?jǐn)?shù)也會(huì)相較航段1 偏低。對(duì)于PELT-SOM 和KCP 算法而言，懲罰系數(shù)越小，意味著允許生成更多的斷點(diǎn)，一方面會(huì)使R分?jǐn)?shù)在一定程度上增大，但過多的斷點(diǎn)會(huì)增加計(jì)算時(shí)長(zhǎng)，同時(shí)使模型結(jié)果過于復(fù)雜而不利于后續(xù)分析。從圖5（c）和圖6（c）中可以看出，若懲罰系數(shù)過小，算法會(huì)忽視某些必要的過渡模式。就GSS 算法而言，相比上述2 種算法表現(xiàn)相對(duì)更好，但其不支持模式挖掘功能并且需要提前調(diào)節(jié)超參數(shù)，它并不是此場(chǎng)景下最好的選擇。由于AutoPlait 算法設(shè)計(jì)初衷為捕捉更高級(jí)復(fù)雜的模式，且并不能適應(yīng)過渡態(tài)和數(shù)據(jù)的無規(guī)律變化，這導(dǎo)致該算法在2 個(gè)航段上的實(shí)驗(yàn)并不理想。綜合表2 和圖5、圖6 不難看出，AutoMiner 方法在航段1 和航段2 上較其他幾類方法表現(xiàn)更優(yōu)越。

表2 在航段1 和航段2 上分割評(píng)價(jià)指標(biāo)R 分?jǐn)?shù)的對(duì)比Table 2 Comparison of split evaluation indicators R score on Flight segments 1 and 2

3.3.2 算法可遷移性的實(shí)驗(yàn)解釋

為研究AutoMiner 算法在更高級(jí)模式的挖掘能力，隨機(jī)選擇有規(guī)律性的時(shí)序數(shù)據(jù)，改變模型階數(shù)可獲得圖7 中的結(jié)果。由航段3 的實(shí)驗(yàn)可以看出，復(fù)雜的模型可以挖掘出由多個(gè)工作模式組成的更高級(jí)的模式。圖7（c）、圖7（d）中斜線塊所示為高級(jí)的模式，由于未提供相應(yīng)高級(jí)模式的參照數(shù)據(jù)集，其自然地被視為過渡態(tài)模式。模型階次越高意味著算法采用了更復(fù)雜的模型，這可以挖掘出更復(fù)雜的模型狀態(tài)。如圖7（e）的算法中選擇多級(jí)鏈馬爾可夫模型便可以挖掘出圖中M1～M5 所示的更高級(jí)的工作模式。

圖7 在航段3 中算法遷移性實(shí)驗(yàn)Fig.7 Algorithm transportability experiments in Flight segment 3

3.3.3 實(shí)驗(yàn)迭代過程的討論

在未采用并行化措施時(shí)，AutoMiner 算法計(jì)算復(fù)雜度為。從航段2 的實(shí)驗(yàn)可以看出，隨著迭代次數(shù)的增加，每步迭代需要的時(shí)間將快速減少，這是因?yàn)殡S著斷點(diǎn)的增多對(duì)應(yīng)時(shí)序片段長(zhǎng)度較少，加上引入存儲(chǔ)單元也大大縮短了后續(xù)迭代所需的計(jì)算時(shí)間。圖8 中T為采樣步長(zhǎng)，圖例中括號(hào)為與之對(duì)應(yīng)的評(píng)價(jià)指標(biāo)R分?jǐn)?shù)。從圖8 中可以看出適當(dāng)增大采樣步長(zhǎng)會(huì)使算法總耗時(shí)急劇減少，但在一定程度內(nèi)分割指標(biāo)的減少量并不大，所以可以根據(jù)實(shí)際情況適當(dāng)增加采樣步長(zhǎng)犧牲較少的精度以換取較低計(jì)算成本。

圖8 航段2 中算法總耗時(shí)隨迭代次數(shù)的變化Fig.8 Variation of total time consumed by algorithm with number of iterations in Flight segment 2

圖9 為AutoMiner 算法在航段1～航段3 上求解斷點(diǎn)組合時(shí)總代價(jià)隨迭代次數(shù)的變化圖，其中選用的模型默認(rèn)為一階線性模式。為獲得完整的總代價(jià)趨勢(shì)變化圖，將算法4 步驟22 的中止功能暫時(shí)關(guān)閉。在圖9 中變化曲線存在明顯的拐點(diǎn)，在迭代達(dá)到曲線極小值點(diǎn)后，進(jìn)一步增加斷點(diǎn)數(shù)量反而會(huì)使總代價(jià)增大。另外，通過式（10）的分析不難發(fā)現(xiàn)，當(dāng)模型和時(shí)序數(shù)據(jù)間的模型誤差逐漸減小，即足夠小時(shí)，編碼代價(jià)會(huì)在數(shù)值上變?yōu)樨?fù)數(shù)，但這與總代價(jià)最小化的尋優(yōu)策略并不矛盾。

圖9 AutoMiner 算法在迭代過程中總代價(jià)的變化Fig.9 Variation of total cost of AutoMiner algorithm during iteration process

3.3.4 模式挖掘效果的對(duì)比評(píng)估

現(xiàn)今模式挖掘效果的評(píng)估方法常采用準(zhǔn)確度、F1分?jǐn)?shù)，其中F1分?jǐn)?shù)是由混淆矩陣得到的精確率和召回率的調(diào)和平均數(shù)。在航段1～航段3來源更完整的3 組航段上繼續(xù)開展實(shí)驗(yàn)，將Au‐toMiner 算法挖掘的工作模式與PELT-SOM 算法的結(jié)果進(jìn)行比較可以得到表3 所示結(jié)果。從中可以看出過渡態(tài)模式上2 種方法的各項(xiàng)指標(biāo)均弱于穩(wěn)態(tài)模式，而相較之下AutoMiner 方法會(huì)獲得更好的結(jié)果。若時(shí)序數(shù)據(jù)位于2 種工作模式之間，例如圖5 中前1 000 s 對(duì)應(yīng)發(fā)動(dòng)機(jī)停轉(zhuǎn)，PELT-SOM 會(huì)將其與最相似的穩(wěn)態(tài)模式即慢車模式優(yōu)先匹配從而導(dǎo)致判斷失誤。AutoMiner 方法通過標(biāo)記傳播獲得軟標(biāo)記，比較屬于各模式的概率并把不屬于穩(wěn)態(tài)模式的概率視為過渡態(tài)從而解決了此類問題。

表3 模式挖掘評(píng)價(jià)指標(biāo)的對(duì)比Table 3 Comparison of evaluation metrics for pattern mining

為了對(duì)模型獲取時(shí)序片段的工作模式開展可視化分析，選用t-SNE 方法［25］將從原始時(shí)間序列中提取的多維特征降至二維的特征圖上并繪制如圖10 所示散點(diǎn)圖，圖中按模式總量排布并且括號(hào)中表示通過不同模型挖掘出的對(duì)應(yīng)工作模式的時(shí)序片段數(shù)量。圖10（b）中軟標(biāo)記最大值代表對(duì)應(yīng)工作模式的概率值，它在圖10（b）中表示為散點(diǎn)的形狀大小。由圖10 可以看出AutoMiner 模型獲得的工作模式更加接近于真實(shí)的工作模式。

圖10 所有時(shí)序片段的工作模式的可視化分布Fig.10 Visualization of pattern distribution of whole time-series segmentation

4 結(jié)論

1）AutoMiner 方法通過對(duì)模型進(jìn)行編碼，利用編碼代價(jià)最小的思想來獲取最優(yōu)的時(shí)序片段組合，并在分段時(shí)利用標(biāo)記傳播來解決過渡態(tài)模式識(shí)別的問題。相較單變量的方法，本文采用的多變量方法能更好地解決狀態(tài)參數(shù)的時(shí)滯問題。

2）通過在某型發(fā)動(dòng)機(jī)全航段數(shù)據(jù)上的實(shí)驗(yàn)表明，在與諸多算法對(duì)比后，AutoMiner 方法在分割指標(biāo)以及模式挖掘指標(biāo)上均表現(xiàn)優(yōu)異。

3）AutoMiner 方法支持并行化計(jì)算并采取諸多措施加速迭代過程，實(shí)驗(yàn)還證明了該方法在更復(fù)雜工作模式方面具有很好的可遷移性。

附錄A：

1）式（7）對(duì)θi求偏導(dǎo)數(shù)并令之為0，求解θi。

求解對(duì)θi的偏導(dǎo)數(shù)時(shí)，將Σi視為已知常數(shù)矩陣，并對(duì)式（7）的跡取微分可得：

依照微分和偏導(dǎo)數(shù)轉(zhuǎn)化公式并結(jié)合跡的性質(zhì)可將式（A1）化為

對(duì)式（A2）化簡(jiǎn)后得到

令式（A2）為0，即可化簡(jiǎn)出式（8）中θi的結(jié)果。

2）式（7）對(duì)Σi求偏導(dǎo)數(shù)并令之為0，求解Σi。

同樣，求解對(duì)Σi的偏導(dǎo)數(shù)時(shí)，將θi視為已知常數(shù)矩陣，并對(duì)式（7）的跡取微分可得：

式（A4）可化為含A矩陣的形式，可以直觀看出A矩陣為對(duì)稱矩陣，便得到所求偏導(dǎo)數(shù)如式（A5）所示，A矩陣在此情形下可以視為雅可比矩陣。

式中：Id為d階單位矩陣；符號(hào)“°”表示矩陣的哈達(dá)瑪積，即矩陣點(diǎn)乘。若矩陣A由元素aij組成，令式（A5）結(jié)果為0，可得

式（A6）反映出A=0，化簡(jiǎn)可得式（8）中Σi的結(jié)果。