劉柏兵， 宋 東， 李春曉

(西北工業(yè)大學航空工程系 西安，710000)

引 言

隨著現(xiàn)代復(fù)雜系統(tǒng)的功能要求、性能要求和使用環(huán)境要求的不斷提升，對系統(tǒng)的故障預(yù)測和健康管理的要求也在不斷提升。在PHM技術(shù)中，一項基本而重要的任務(wù)就是了解系統(tǒng)的健康狀況及其發(fā)生性能衰退/健康狀態(tài)降級的相關(guān)規(guī)律，積累系統(tǒng)的性能衰退/健康降級模式，進而建立其全面的系統(tǒng)衰退模型[1]。建立衰退模型的一種方法，就是對狀態(tài)監(jiān)測數(shù)據(jù)進行時間序列數(shù)據(jù)挖掘[2]。

國外較早開展該研究的機構(gòu)是IBM公司的Pazzani和Agrawal 研究小組。在時間序列特征表示方面，目前已有不少方法，其中分段線性表示就是一種簡單直觀的特征表示方法，它使用線性模型對時間序列進行分割， 根據(jù)不同的分割方法可以使用不同的分割策略來實現(xiàn)[3]， 如滑動窗口、自底向上和自頂向下?；瑒哟翱诤妥缘紫蛏戏椒ǖ臅r間復(fù)雜度為序列長度的平方階， 而自頂向下的時間復(fù)雜度為線性階。其中針對滑動窗口方法對時間序列擬合效果不佳這一問題，Shatkay作了對應(yīng)的解釋和分析。Park等[4]提出通過采用單調(diào)變換來改進該方法，雖然能有效地分割平滑的時間序列，但對時間序列含有噪聲則無能為力。李愛國等[5]提出了一種對時間序列進行在線分割的遞推算法，能夠?qū)崟r發(fā)現(xiàn)和預(yù)測時態(tài)模式，包括了分段聚合近似方法、符號化表示方法、基于域變換的表示方法、奇異值分解和基于模型的表示方法等。

針對基于狀態(tài)監(jiān)測數(shù)據(jù)的衰退模式挖掘中存在的噪聲、序列時域形變或漂移不識別等問題，筆者提出了一種P-D-H聚類方法,其中：P是指通過分段聚合近似(PAA)方法對由狀態(tài)監(jiān)測數(shù)據(jù)形成的退化軌跡時間序列進行模式表示，這樣能夠有效地對時間序列進行降維，在降低了時間消耗的基礎(chǔ)上充分反映了時間序列的信息，提高算法的效率；D是指采用動態(tài)時間彎曲距離(DTW)代替歐式距離作為模式序列的相似性度量，克服了歐式距離當序列在時間軸上發(fā)生輕微偏移時變化很大的缺點，提高了算法的準確性；H是指采用層次聚類方法實現(xiàn)衰退模式聚類，可以表征聚類的全過程，并且詳盡展示了聚類簇間的距離關(guān)系。該方法充分利用現(xiàn)有數(shù)據(jù)資源，對系統(tǒng)衰退模式進行挖掘，為系統(tǒng)的健康預(yù)測提供依據(jù)，具有一定的理論價值和應(yīng)用價值。

1 基于P-D-H的系統(tǒng)衰退模式挖掘方法

1.1 退化軌跡的建立

產(chǎn)品在存儲或使用中發(fā)生的性能衰退表現(xiàn)為產(chǎn)品的某些性能特征參數(shù)不斷退化并超過其規(guī)定范圍，這種用于描述產(chǎn)品衰退過程的性能特征參數(shù)變化曲線/序列，稱為觀測產(chǎn)品的“退化軌跡”[6]。基于狀態(tài)數(shù)據(jù)的退化軌跡P-D-H聚類方法的衰退模式挖掘，首先就是要從狀態(tài)監(jiān)測數(shù)據(jù)中提取系統(tǒng)的健康特征參量，特征參量形成的時間序列趨勢成分，可以作為退化軌跡表征系統(tǒng)的衰退過程[7]。處理過程如圖1所示(其中MA表示采用移動平均的方法)：a.對每組信號進行特征提取，得到不同特征的變化序列；b.提取特征序列，采用平滑法進行趨勢分析；c.確定健康特征，獲得退化軌跡。

圖1 退化軌跡建立流程圖Fig.1 Deletion trajectory build process

此處給出本研究所采用特征時間序列趨勢成分提取方法以及斯皮爾曼等級相關(guān)系數(shù)衰退特征選取方法。

1.1.1移動平滑法

假設(shè)時間序列

F(t)={Fi(t)|i=1,2,…,m;t=1,2,…,n}

是狀態(tài)監(jiān)測參數(shù)在t個時間點下，m個特征參數(shù)的時間序列集，它可以表示為Fi(t)=Xi(t)+μi(t)，其中：Xi(t)為特征序列Fi(t)采用移動平均平滑法進行處理得到的趨勢成分序列；μi(t)為去除趨勢成分后的Fi(t)殘差序列。則平滑后的趨勢序列Xi(t)滿足

(1)

1.1.2斯皮爾曼等級相關(guān)系數(shù)衰退特征選取方法

斯皮爾曼等級相關(guān)系數(shù)適合于定序變量或不滿足正態(tài)分布假設(shè)的等間隔數(shù)據(jù)，它利用兩變量的秩次大小作線性相關(guān)分析，對原始變量的分布不作要求，屬于非參數(shù)統(tǒng)計方法，適用范圍比較廣。 采用斯皮爾曼等級相關(guān)系數(shù)ρ作為衰退特征選取指標的思路是，計算狀態(tài)監(jiān)測參數(shù)X(t)={Xi(t),i=1,2,…,m}中每個Xi(t)與時間序列t=1,2,…,n之間的斯皮爾曼等級相關(guān)系數(shù)ρi，選取ρi值最大的若干Xi(t)作為系統(tǒng)的退化軌跡，記為Di(t)。實際應(yīng)用中，變量接的連結(jié)是無關(guān)緊要的，于是可以通過簡單的步驟計算ρ。對于樣本容量為n的樣本，n個原始數(shù)據(jù)Xi和Yi被轉(zhuǎn)換成等級數(shù)據(jù)xi和yi[8],得到被觀測的兩個等級差值di=xi-yi，則ρ為

(2)

1.2 退化軌跡集聚類

通過上述退化軌跡建立方法對狀態(tài)監(jiān)測數(shù)據(jù)建立退化軌跡后，得到樣本個數(shù)為N的監(jiān)測數(shù)據(jù)退化軌跡時間序列集D(t)={Di(t),i=1,2,…,N}，采用PAA的方法對退化軌跡D(t)進行近似表征，然后通過DTW度量近似序列之間的距離，以此作為退化軌跡序列之間的相似度度量準則，最后采用自底向上的層次聚類算法對退化軌跡集合中的樣本進行聚類。

1.2.1退化軌跡的分段聚合近似

PAA是一種時間序列的基于劃分的近似表示，其意義簡單明確，較好地保持了時間序列在時域上的變化趨勢特征，降低了時間序列中包含的噪聲信息。由于本研究中的衰退模式挖掘目的在于發(fā)現(xiàn)時間序列變化趨勢的相似性，采用分段線性表示的方法對時間序列降維和特征表示是一種可行的方法[9]。將PAA近似表征方法與DTW距離度量相結(jié)合，可以克服范數(shù)距離度量對序列時域形變或漂移不識別的問題，其方法原理如下。

(3)

為了使數(shù)據(jù)的維數(shù)由n維向量降為N維向量,首先需要將數(shù)據(jù)分割成等長的N小段,之后需要計算出每一段的均值,再用這些數(shù)值重新組成向量作為時間序列降維后的表示[10]。

1.2.2DTW距離

在進行時間序列樣本相似性度量的時候，最基本距離函數(shù)主要為普通范數(shù)距離和動態(tài)時間彎曲距離(DTW)。普通范數(shù)距離Lp在計算時要求兩條時間序列等長，且計算方法嚴格按照點與點的順序一一對應(yīng)，對不等長的時間序列無法計算其普通范數(shù)距離。與普通的范數(shù)距離不同，DTW在聚類時不要求時間序列之間的點與點匹配，當兩條序列長度不同時，允許時間序列點自我復(fù)制后再進行對齊匹配。DTW對時間序列的突變或異常點不敏感，更加適用于此類時間序列的距離度量，且可以實現(xiàn)異步相似性比較，普通范數(shù)距離只能實現(xiàn)同步比較。因此，筆者采用DTW距離作為相似性度量的準則,其原理如下。

對于兩條退化軌跡序列Q=q1,q2,…,qn和C=c1,c2,…,cm，定義一個n行m列矩陣

(4)

將兩個時間序列分別置于二維坐標的兩軸，如圖2所示，進而可以定義彎曲路徑：在兩個不同時間序列間的距離矩陣中，定義時間序列間相異性關(guān)系的K個連續(xù)的矩陣元素的集合，稱為彎曲路徑。

W=w1,w2,…,wi,…,wk

(5)

圖2 彎曲路徑示意圖Fig.2 Schematic diagram of curved path

將退化軌跡序列Q，C之間的動態(tài)時間彎曲距離DTW定義為兩序列距離矩陣中彎曲路徑最小的路徑長度，稱為最優(yōu)彎曲路徑，即

(6)

其中：分母K是為了在比較不同長度的路徑時有統(tǒng)一的標準。

DTW距離的求解采用累積距離矩陣的動態(tài)規(guī)劃方法來計算。累積距離矩陣實際上是一種遞推關(guān)系，對于上文中定義的兩個時間序列Q和C，累積距離定義為

(i=2,3,…,n;j=2,3,…,m)

(7)

DTW距離為兩條時間序列對齊匹配關(guān)系的最佳值給出了定義，同時為時間序列的時間彎曲和不同長度的時間序列之間的相似性度量提供了支持[11]。

1.2.3P-D-H層次聚類算法流程

本研究采用的層次聚類不需要預(yù)先設(shè)定聚類數(shù)目，允許從任意維度對數(shù)據(jù)集進行劃分，它的實現(xiàn)建立于對數(shù)據(jù)集中的兩兩樣本進行相似性(距離)計算的基礎(chǔ)上，通過不斷合并相似性最高(距離最小)的樣本來構(gòu)成新的集群，將合并得到的集群看成一個新的樣本，再次進行相似性計算和樣本的合并，如此反復(fù)進行運算，直至所有樣本合并為一類。層次聚類方法的優(yōu)點在于可以表征聚類的全過程，并且詳盡地展示了聚類簇間的距離關(guān)系[12]。其具體流程如下：

1) 對輸入樣本集D={x1,x2,…,xm},定義樣本間距離函數(shù)dist(xi,xj)和類間距離函數(shù)d(Ci,Cj);

2) 將樣本集中的每一個樣本化為一簇，即Cj={xj}；

3) 對于聚類集合Cj={xj}進行兩兩的距離計算，構(gòu)建距離矩陣M，M(i,j)=M(j,i)=d(Ci,Cj)；

4) 找出距離最近的兩個聚類簇Ci*和Cj*，形成新的聚類簇及編號；

5) 刪除距離矩陣M的相應(yīng)行列，計算新的聚類簇與其他簇的距離，更新距離矩陣；

6) 循環(huán)進行第3步～第5步，直到全部樣本聚為一類；

7) 確定劃分的聚類數(shù)目k。

這里樣本間距離為dist(xi,xj)，采用DTW距離計算方法，記為distDTW；類間距離d(Ci,Cj)采用關(guān)于集合的某種距離即可，可選用的類間距離有最小距離、最大距離和平均距離。

由上述可知，最小距離由兩個類別間的最近樣本決定，最大距離由兩個類別間的最遠樣本決定，而平均距離則由兩個類別中的所有樣本共同決定。最小距離和最大距離的計算比較簡單，但對簇的整體特性刻畫比較片面，平均距離能更好地刻畫簇的特征，但計算量稍大一些。

1.2.4數(shù)據(jù)標準化

1) Z-score標準化

(8)

其中：μ和σ分別為X的均值和標準差。

2) Min-max標準化

(9)

3) 規(guī)范化

(10)

2 試驗分析

通過對軸承磨損衰退模式挖掘來驗證基于P-D-H的系統(tǒng)衰退模式挖掘方法的有效性，本節(jié)的數(shù)據(jù)來源于IEEE可靠性學會與FEMTO-ST Institute共同組織的2012IEEE PHM大會，目的是用于軸承剩余壽命預(yù)測比賽。數(shù)據(jù)采集于FEMTO-ST Institute的PRONOSTIA試驗臺[13]進行的17個軸承加速退化試驗[14]。

試驗在數(shù)小時內(nèi)完成了17個軸承的加速退化，失效判據(jù)統(tǒng)一為監(jiān)測量的振動加速度，超過20g時停止試驗。工況與試驗件的信息如下。

工況1：1800r/min，4000N；試驗軸承編號1-1～1-7，共7個；

工況2：1650r/min，4200N；試驗軸承編號2-1～2-7，共7個；

工況3：1500r/min，5000N；試驗軸承編號3-1～3-3，共3個。

采樣方式和參數(shù)單位如表1所示，即振動加速度信號每間隔10s進行1次采集，每次記錄2560個數(shù)據(jù)點；溫度信號每秒記錄1個數(shù)據(jù)點。

表1 采集參數(shù)詳情表

2.1 退化軌跡的建立

2.1.1 特征提取

基于對振動加速度信號的觀測和一般了解，采取如下的特征提取方法：首先，對每組信號直接提取2個有量綱統(tǒng)計特征和2個為無量綱統(tǒng)計特征，包括均方根(root mean square,簡稱RMS)、峰值(peak)、峰值因子(crest factor，簡稱CF)、譜峭度(spectrum kurtosis，簡稱SK)；其次，對振動信號進行帶寬為500 Hz的數(shù)字帶通濾波，得到頻段0～12 kHz的24個濾波信號，對每個濾波信號同樣提取4個統(tǒng)計量。峰值特征在計算的時候，將每次觀測數(shù)據(jù)中的3個最大值進行平均來得到峰值特征，以消除異常值的影響[15]。各個統(tǒng)計特征量計算公式如表2所示。

表2 統(tǒng)計特征量計算方法

表2中：xi為觀測值；N為觀測點數(shù)，此處為2 560；i∈[1,N]，為樣本均值；Xord(ti)為降序排列的觀測序列。

2.1.2 平滑趨勢

對特征序列進行平滑處理，以得到其平滑趨勢。以軸承1-1為例，觀測記錄434組數(shù)據(jù)，對未濾波的單次2 560個采樣值提取以上4個統(tǒng)計量，得到的特征序列如圖3所示。此處選取窗寬n=20，平滑后序列計算方法參見式(1)，平滑后的趨勢序列如圖4所示。

圖3 軸承1-1特征時間序列Fig.3 Bearing 1-1 characteristic time series

圖4 軸承1-1特征平滑時間序列Fig.4 Bearing 1-1 features smooth time series

2.1.3 退化軌跡的建立

通過式(2)及相關(guān)定義，對平滑趨勢后的特征序列與其對應(yīng)的記錄時間進行斯皮爾曼等級相關(guān)性系數(shù)計算。綜合所有數(shù)字特征和相關(guān)系數(shù)計算結(jié)果，發(fā)現(xiàn)均方根和峰值因子的與時間序列相關(guān)度系數(shù)均在0.5以下，因此不作為特征參數(shù)參與建立退化軌跡。在所有軸承信號分析中，峰值信號與時間的相關(guān)度系數(shù)值大多保持在0.8以上，表明了其良好的單調(diào)性，在軸承1～3中高達0.893 2，且各個軸承峰值時間序列數(shù)字范圍均在0～20之內(nèi)，具備較好的統(tǒng)一性。譜峭度與時間的相關(guān)性低于峰值信號，保持在0.6以上，數(shù)值范圍起于3止于4。因此選取未濾波信號峰值(以下簡稱峰值P)和5.5～6 kHz頻段下濾波信號的譜峭度(以下簡稱譜峭度K)作為兩個不同的健康特征量。這也符合工程實際中的一般認識：軸承在健康狀態(tài)下的振動信號是服從高斯分布的，其譜峭度值為3；隨著試驗的進行，軸承的磨損逐漸加重，其振動信號的譜峭度值越來越偏離3，表現(xiàn)出漸變的退化趨勢。圖5和圖6分別為軸承1-1水平振動加速度信號峰值特征序列和譜峭度特征序列。

圖5 峰值特征平滑序列Fig.5 Peak feature smoothing sequence

圖6 譜峭度特征平滑序列Fig.6 Spectrum kurtosis feature smoothing sequence

2.2 聚類過程

選取了兩種特征量作為健康表征參數(shù)，故對兩種特征參數(shù)下的樣本退化軌跡分別進行聚類，聚類的總樣本個數(shù)為34，樣本時間序列長度不完全相等。分段聚合數(shù)目N=100聚類過程按照2.2節(jié)所述的聚類過程，這里直接給出得到兩種特征下的樣本聚類譜系分別如圖7和圖8所示。

圖7 基于譜峭度的聚類譜系圖Fig.7 Based on the clustering spectrum of spectral kurtosis

圖8 基于峰值的聚類譜系圖Fig.8 Based on the clustering spectrum of peak

由圖7、圖8可觀察到，在譜峭度和峰值特征表征下，大多數(shù)的樣本退化軌跡比較相似，基于譜峭度特征的樣本在更小的距離上聚合，離群樣本與其他樣本差異性較大；而基于峰值特征的聚類對樣本區(qū)分度更加明確。綜合以上兩圖，將最終聚類數(shù)目確定為4，得到聚類劃分結(jié)果如表3和表4所示。

表3 基于譜峭度的模式聚類結(jié)果表

表4 基于峰值的模式聚類結(jié)果表

由于兩種不同的特征參數(shù)與不同故障模式相關(guān)系數(shù)有所不同，且試驗中可能存在噪聲的干擾，因此，兩種不同特征參數(shù)下的模式3和模式4聚類的結(jié)果可能存在差異。

圖9選用樣本1,6,20,23來分別表示4種衰退模式之間的差異，其中4條不同顏色的曲線分別代表了不同的樣本，圖中橫軸為觀測時間點，縱軸為SK值。同理，圖10選用樣本1,6,19,20來分別表示4種衰退模式之間的差異，圖中橫軸為觀測時間點，縱軸為峰值。

圖9 基于譜峭度的衰退模式樣本退化軌跡Fig.9 Sample degradation trajectory of degradation pattern based on spectrum kurtosis

圖10 基于峰值的不同模式樣本退化軌跡Fig.10 Sample degradation trajectory of different patterns based on peaks

由聚類結(jié)果可以看到，樣本6和樣本20在兩種特征下的退化軌跡都表現(xiàn)出與大多數(shù)樣本不一致的特性。在兩種健康特征下，聚類獲得的4種退化模式下的退化曲線均有明顯的形狀差別。樣本6基于譜峭度的退化軌跡與基于峰值的退化軌跡區(qū)別明顯，得到的退化軌跡和聚類結(jié)果可用于建立系統(tǒng)的衰退模型，具有較好的效果。

3 結(jié)束語

筆者 針對PHM技術(shù)中系統(tǒng)衰退模型的建立進行了深入的研究，對比了多種現(xiàn)有的基于狀態(tài)監(jiān)測數(shù)據(jù)的退化軌跡建立和時間序列聚類的模式挖掘方法，總結(jié)了其中的缺點和不足。在此基礎(chǔ)上提出了一種P-D-H聚類方法，利用此方法對IEEE可靠性學會與FEMTO-ST Institute共同組織的軸承剩余壽命預(yù)測比賽數(shù)據(jù)進行了衰退模式挖掘，獲得了4種退化模式下的退化曲線，驗證了方法的有效性，為狀態(tài)監(jiān)測數(shù)據(jù)的系統(tǒng)健康衰退模式挖掘提供了一種有效的方法。