考慮突變狀態(tài)檢測的齒輪實時剩余壽命預(yù)測

石 慧， 曾建潮,2

(1.太原科技大學(xué) 工業(yè)與系統(tǒng)工程研究所，太原 030024； 2.中北大學(xué) 計算機與控制工程學(xué)院，太原 030051)

考慮突變狀態(tài)檢測的齒輪實時剩余壽命預(yù)測

石 慧1， 曾建潮1,2

(1.太原科技大學(xué) 工業(yè)與系統(tǒng)工程研究所，太原 030024； 2.中北大學(xué) 計算機與控制工程學(xué)院，太原 030051)

為解決齒輪疲勞退化過程中狀態(tài)突變后剩余壽命難以準確預(yù)測問題，提出一種考慮退化突變點檢測與剩余壽命預(yù)測相關(guān)聯(lián)的齒輪疲勞實時剩余壽命預(yù)測新方法。針對齒輪磨損退化過程建立狀態(tài)空間預(yù)測模型，利用接收到的齒輪實時監(jiān)測振動信息實時更新模型參數(shù)，同時對退化過程中的突變狀態(tài)點進行檢測，并根據(jù)突變點所提供的壽命信息采用卡爾曼前向濾波及平滑算法結(jié)合期望最大化參數(shù)估計算法在濾波的同時不斷對狀態(tài)空間模型參數(shù)進行修正，改變退化突變后的濾波效果，進行實時狀態(tài)預(yù)測與壽命估計。運用齒輪疲勞壽命試驗臺的實時監(jiān)測數(shù)據(jù)對預(yù)測模型進行驗證，結(jié)果表明利用突變點信息對預(yù)測模型進行修正后可以更快的對系統(tǒng)的動態(tài)變化進行跟蹤，提高預(yù)測齒輪退化狀態(tài)及實時剩余壽命的準確度。

剩余壽命預(yù)測；狀態(tài)空間建模；卡爾曼濾波；突變狀態(tài)檢測；模型修正

齒輪是如汽輪機、離心壓縮機、風(fēng)機等工業(yè)部門中應(yīng)用最為廣泛的旋轉(zhuǎn)機械設(shè)備傳動系統(tǒng)中的關(guān)鍵部件。齒輪嚙合受力情況較為復(fù)雜，在變工況、變載荷的情況下運行，更容易發(fā)生故障。當齒輪發(fā)生斷齒、齒面疲勞、膠合等故障時，常常會引起整個機械設(shè)備災(zāi)難性的破壞[1]。同時隨著信息傳感設(shè)備的發(fā)展，可對齒輪的運行狀態(tài)進行實時監(jiān)測，傳統(tǒng)的故障診斷方法雖然可以實現(xiàn)齒輪故障的有效診斷，但存在一定的局限性。因此對齒輪這類易磨損部件的退化狀態(tài)預(yù)測與健康管理已成為業(yè)界日益關(guān)注的問題。利用接收到的大量實時監(jiān)測信息更為準確預(yù)測系統(tǒng)的退化狀態(tài)及其剩余壽命，可以提供有關(guān)健康狀態(tài)的關(guān)鍵信息，進而識別和管理故障的發(fā)生、規(guī)劃維修活動，為更合理地制定基于狀態(tài)的維護維修策略提供依據(jù)[2-4]。

齒輪的磨損退化過程是一個動態(tài)隨機過程。近年來隨著信息傳感技術(shù)的發(fā)展，可對齒輪進行實時監(jiān)測，并利用接收到大量的反映設(shè)備退化性能的監(jiān)測信息進行更準確的剩余壽命預(yù)測，因此基于統(tǒng)計學(xué)理論和人工智能理論進行剩余壽命預(yù)測建模的數(shù)據(jù)驅(qū)動剩余壽命預(yù)測方法被廣泛應(yīng)用于齒輪的預(yù)測與健康管理領(lǐng)域。它包括狀態(tài)空間模型、隨機濾波模型、自回歸滑動平均模型、支持向量機、人工神經(jīng)網(wǎng)絡(luò)及其衍生模型和灰色模型等剩余壽命方法。但在系統(tǒng)運行過程中，可監(jiān)測到的數(shù)據(jù)往往只是與狀態(tài)有關(guān)的一些間接觀測數(shù)據(jù)，例如振動信號、溫度和壓力等信號。而且觀測數(shù)據(jù)往往因受到隨機噪聲的干擾而出現(xiàn)觀測誤差。根據(jù)觀測數(shù)據(jù)構(gòu)造狀態(tài)空間模型對系統(tǒng)的狀態(tài)進行估計與預(yù)測可以很好地解決這一問題。卡爾曼濾波是一種線性最優(yōu)無偏估計器，隨著監(jiān)測數(shù)據(jù)的增加，濾波精度越來越高，可對突變狀態(tài)進行跟蹤。Ga?perin等[5]針對齒輪的磨損退化建立非線性狀態(tài)方程，進行線性化并采用卡爾曼濾波狀態(tài)估計。Carr等[6]采用擴展卡爾曼濾波的方式進行設(shè)備退化狀態(tài)及剩余壽命預(yù)測。并假設(shè)剩余壽命與退化狀態(tài)之間滿足對數(shù)關(guān)系進而估計剩余壽命分布，并制定預(yù)防性替換的維修策略。Sun等[7]采用卡爾曼濾波的方式進行設(shè)備剩余壽命預(yù)測。提出健康指標(Health Index，HI)的概念，設(shè)備新時HI=1，發(fā)生故障HI=0。采用時間序列蒙特卡洛仿真(SMC)進行剩余壽命預(yù)測。Bressel等[8]提出擴展卡爾曼濾波方式進行質(zhì)子交換膜燃料電池退化狀態(tài)預(yù)測，但其降低模型不確定度，不能很好跟蹤系統(tǒng)的動態(tài)變化。

還有一些文獻提出了兩階段和多階段的退化狀態(tài)預(yù)測模型[9]。Lall等[10]考慮針對正常和潛在故障兩個退化狀態(tài)，其中正常狀態(tài)不進行預(yù)測，潛在故障階段建立卡爾曼預(yù)測模型，進行剩余壽命預(yù)測。但文中退化閾值是給定的，沒有針對退化的突變狀態(tài)進行識別檢測及動態(tài)跟蹤。Wang等[11]考慮設(shè)備連續(xù)運行過程，假設(shè)已通過各種統(tǒng)計方法檢測到退化異常點，分為正常退化和潛在故障兩個階段，認為在潛在故障階段退化速度增大，僅在此階段預(yù)測退化狀態(tài)求其剩余壽命。而齒輪的退化過程隨著設(shè)備負載及工況的變化而變化，開始齒輪嚙合階段與后期磨損加速階段是復(fù)雜多變的多階段退化過程，沒有恒定不變的階段，難以用一個退化模型來描述。胡昌華等[12]提出了一種多階段的預(yù)測模型，但其僅是根據(jù)歷史監(jiān)測數(shù)據(jù)分段擬合退化模型建立含可變參數(shù)的狀態(tài)預(yù)測模型，但在擬合分布過程中需要大量破壞性實驗數(shù)據(jù)，很難收集或成本太高。同時系統(tǒng)運行中其退化狀態(tài)很難直接測量，李鑫等[13]提出多種退化模式下動態(tài)轉(zhuǎn)移的健康狀態(tài)自適應(yīng)預(yù)測，由健康狀態(tài)評估及選定的退化閾值決定改變退化模式進而采用不同的剩余壽命預(yù)測方式。但其考慮退化階段時假設(shè)已通過統(tǒng)計方法檢測出異常突變點，直接基于新的階段進行壽命預(yù)測及維護維修策略的制定，沒有考慮檢測出異常點后如何根據(jù)突變點提供的信息相應(yīng)改變預(yù)測方法以更好地跟蹤剩余壽命的動態(tài)變化。

在系統(tǒng)實際運行過程中，往往由于其退化磨損過程由漸變、量變發(fā)展為突變、質(zhì)變，導(dǎo)致系統(tǒng)的材料、結(jié)構(gòu)等發(fā)生變化，出現(xiàn)退化突變點。突變點后的退化速率增大進而使剩余壽命發(fā)生較大變化。突變點是退化狀態(tài)發(fā)生變化的起始點，它所包含的大量剩余壽命信息能夠為剩余壽命預(yù)測及預(yù)防性維護維修策略的制定提供有用信息[14-15]。本文針對齒輪疲勞壽命試驗的實時監(jiān)測數(shù)據(jù)建立可變參數(shù)的系統(tǒng)退化狀態(tài)空間模型，同時考慮異常點檢測與剩余壽命預(yù)測的關(guān)聯(lián)性，對齒輪退化過程中的突變狀態(tài)點進行檢測，突變點后采用卡爾曼前向濾波及平滑算法結(jié)合期望最大化參數(shù)估計算法在濾波的同時不斷對未知的系統(tǒng)模型參數(shù)進行修正，增加突變點后觀測數(shù)據(jù)的權(quán)重，減小突變點之前的監(jiān)測信息產(chǎn)生的誤差。同時利用檢測到的突變點信息修正狀態(tài)初始值及參數(shù)初始值的選擇，在突變點修正濾波增益，通過大的濾波增益及濾波估計誤差改變?yōu)V波效果，重新進行狀態(tài)估計及壽命預(yù)測。將基于退化狀態(tài)突變點檢測和剩余壽命的預(yù)測進行聯(lián)合研究，可以利用突變點信息更快的對系統(tǒng)的動態(tài)變化進行跟蹤，提高齒輪退化狀態(tài)預(yù)測精度。

1 齒輪疲勞試驗

1.1疲勞試驗臺架

齒輪疲勞壽命試驗采用了如圖1所示的封閉試驗臺架。試驗臺的中心距為150 mm，電機轉(zhuǎn)速為1 200 r/min。試驗過程對箱體振動、油溫和噪音等進行監(jiān)測。

圖1 齒輪試驗臺架

疲勞試驗臺原理圖如圖2(a)所示。本試驗共布置11個傳感器，主試箱傳感器位置圖如圖2(b)所示。加速度傳感器安裝在主試箱的軸承座位置，這樣可以使得接收到的振動信號的衰減最小[16]，在齒輪箱內(nèi)安裝溫度傳感器，在主試箱的正上方安裝噪聲傳感器。1#～8#為加速度傳感器(1#～4#布置在主試箱軸承座的徑向，7#和8#布置在主試箱的軸向，5#和6#布置在陪試箱軸承座的徑向)；9#和10#為聲傳感器，分別布置在主試箱和陪試箱正上方約40 cm處；11#表示溫度傳感器，布置在主試箱體內(nèi)，試驗中測試潤滑油溫度。采樣頻率25.6 kHz，采樣時間60 s，采樣間隔9 min。本試驗中采用常規(guī)成組法即恒載的方式進行，轉(zhuǎn)矩為822.7 N·M。當試驗齒輪發(fā)生斷齒時即判定該齒輪失效，如圖3所示。

(a) 齒輪試驗臺架

(b) 主試箱傳感器

Fig.2 The schematic of test bench and sensors’ placement scheme of the main gearbox

圖3 主試箱齒輪發(fā)生故障斷齒

1.2試驗齒輪的安裝

試驗中采用材料為合金鋼，齒面硬度為58-61HRC的硬齒面齒輪，表面處理為滲碳淬火。采用正反面交錯搭接嚙合方式(如圖4所示)，電機轉(zhuǎn)速為1 200 r/min。主試箱齒輪模數(shù)m=3，齒數(shù)為z1=z2=50，壓力角α=20°，齒寬29 mm，實際工作齒寬13～14 mm；陪試箱齒輪齒數(shù)為z3=z4=24，壓力角α=20°，工作齒寬20 mm。因而齒輪嚙合頻率有2個，分別為1 000 Hz(主試箱齒輪)和480 Hz(陪試箱齒輪)。本次實驗潤滑油采用L-CKC320工業(yè)閉式齒輪油。

圖4 試驗齒輪正反交錯搭接安裝圖

2 特征提取

對齒輪動態(tài)監(jiān)測信息采用有效的方法進行疲勞狀態(tài)的特征提取是實現(xiàn)狀態(tài)預(yù)測和維護維修的關(guān)鍵[17]。實際應(yīng)用中很難采集到表示系統(tǒng)特性的直接狀態(tài)數(shù)據(jù)，隨著傳感技術(shù)的發(fā)展，可接收到準確的振動信號，振動信號分析已成為齒輪疲勞狀態(tài)特征提取最有效的方法之一。

目前，較常用的振動信號分析方法包括原始時域、包絡(luò)時域、原始頻譜、解調(diào)頻譜以及細化頻譜分析法，好的評價指標不僅能真實地反映齒輪運行過程中狀態(tài)性能的變化，而且易于計算。本文利用均方幅值(Root Mean Square, RMS)對齒輪磨損退化性能進行衰退評估。均方幅值作為有量綱的統(tǒng)計特征值，不受齒輪個體差異的影響，隨疲勞狀態(tài)改變呈現(xiàn)遞增趨勢，能較好地反映各采樣時刻振動能量動態(tài)的變化，而且易于計算。對于每次采樣時間Δt長度內(nèi)，離散隨機信號的時間序列均方幅值可表示為：

(1)

式中：Δt為采樣時間；Fs為采樣頻率；n為采樣點數(shù)，n=Fs×Δt。

3 基于狀態(tài)空間模型的退化狀態(tài)建模

在壽命預(yù)測中，如何建立一個適當?shù)哪Ｐ褪鞘滓獑栴}。可以通過非線性狀態(tài)空間建模進行狀態(tài)估計。模型的一般形式為：

Xt=f(Xt-1)+wt-1

(2)

Yt=h(Xt)+vt

(3)

式中：Xt為t時刻(t≥1)的系統(tǒng)狀態(tài)向量；f(·)為系統(tǒng)狀態(tài)轉(zhuǎn)移模型；Yt為t時刻的系統(tǒng)觀測向量；h(·)為系統(tǒng)觀測模型。式(2)稱為狀態(tài)轉(zhuǎn)移方程，式(3)稱為觀測方程。系統(tǒng)噪聲wt和觀測噪聲vt服從零均值的高斯分布(其協(xié)方差陣為：Var(wt)=Qt，Var(vt)=Rt)。

非線性模型可以很好的描述動態(tài)過程及一系列廣泛的動態(tài)行為。然而，在實際應(yīng)用中對非線性系統(tǒng)的分析和辨識相對較復(fù)雜。式和式在許多情況下，通過采用線性化進行簡化，線性化后的狀態(tài)空間模型為：

Xt=Ftt-1Xt-1+wt-1

(4)

Yt=HtXt+vt

(5)

4 突變狀態(tài)檢測

當系統(tǒng)運行一段時間，退化狀態(tài)可能發(fā)生突變，退化速率增大，退化狀態(tài)向上偏移，如圖5所示。利用突變點信息對預(yù)測系統(tǒng)剩余壽命的狀態(tài)空間模型進行修正時首先需要進行退化狀態(tài)突變點的檢測。突變狀態(tài)的檢測定義為在某些未知時刻檢測到隨機序列分布發(fā)生變化。其廣泛應(yīng)用于計算機網(wǎng)絡(luò)入侵檢測和安全系統(tǒng)以及各類設(shè)備的故障檢測中。

圖5 退化狀態(tài)偏移的突變狀態(tài)點檢測示意圖

累積和控制圖(Cumulative Sum, CUSUM)是基于似然比的一種接近最優(yōu)的非參數(shù)化統(tǒng)計方法，對事件的突變狀態(tài)進行累積，將過程中的小偏移量累加起來，求其累積進而判斷是否發(fā)生突變。該算法要求的假定條件較少，能有效反映過程變化的靈敏性，非常適合用在系統(tǒng)退化過程中的突變檢測[18]。

max(yn-μ0-k)

(6)

(7)

c=-(Δμ)(h/σ+1.166)=

-2k(h/σ+1.166)

(8)

5 退化狀態(tài)預(yù)測模型修正

如圖6所示，假設(shè)退化狀態(tài)不發(fā)生突變，A點之后系統(tǒng)退化狀態(tài)為Ⅰ曲線。但狀態(tài)突變點A之后退化速率增加，實際退化狀態(tài)可能變?yōu)镮II曲線。系統(tǒng)接收到大量監(jiān)測信息，利用累積求和算法(CUSUM)檢測到狀態(tài)突變點A，盡管隨著時間的推移，實時監(jiān)測數(shù)據(jù)的增多，濾波估計的精度越來越高?？柭鼮V波值將以曲線II的趨勢逐漸穩(wěn)定并收斂于真實值。但在收斂的過程中對突變點之后的退化狀態(tài)無法進行快速跟蹤，會出現(xiàn)狀態(tài)估計預(yù)測誤差Δx，影響實時剩余壽命的預(yù)測和系統(tǒng)維護維修決策的制定。

圖6 退化狀態(tài)預(yù)測模型修正示意圖

由以上分析可知采用卡爾曼濾波對系統(tǒng)進行狀態(tài)估計及預(yù)測時，當系統(tǒng)出現(xiàn)狀態(tài)突變點后，仍然存在以下問題：

(2)在正常的卡爾曼濾波遞推過程中，由于初始值的準確值不能確切知道，通常選取都是不準確的，只能假定一初值。當濾波時間充分長后，它的卡爾曼最優(yōu)濾波值將漸近穩(wěn)定而不依賴于濾波初值的選擇，但初始值選擇誤差大，其收斂于真值的速度慢。系統(tǒng)退化狀態(tài)出現(xiàn)突變點后，初始值選擇不準確會增加狀態(tài)估計與預(yù)測過程中的誤差。

(3)在卡爾曼濾波遞推過程中，開始的幾步遞推計算中得到的狀態(tài)估計值十分粗略，這一階段稱為粗估階段。在卡爾曼濾波粗估階段，Kt的值和濾波估計誤差Pt很大，可以對狀態(tài)預(yù)測進行較好的修正；之后進入精估階段，Kt和Pt的值逐漸減小，使得觀測值對狀態(tài)估計值的修正作用變小。

因此，采用狀態(tài)空間模型進行狀態(tài)估計與預(yù)測時，在突變狀態(tài)之后對狀態(tài)空間模型進行修正，期望觀測信息對狀態(tài)估計值產(chǎn)生較強的修正作用，使得狀態(tài)估計較快地收斂于部件退化的真實狀態(tài)。可以更好的對突變狀態(tài)進行跟蹤，減少預(yù)測誤差，提高收斂速度和預(yù)測精度。需對模型進行以下修正：

(1)狀態(tài)估計時增加突變點后的監(jiān)測信息，修正系統(tǒng)模型參數(shù)，進行自適應(yīng)的狀態(tài)估計。

(2)利用突變點提供的信息在突變點后進行系統(tǒng)模型參數(shù)修正的同時，重新選擇濾波初始值。使得狀態(tài)預(yù)測值較快的接近被估計的真實狀態(tài)。

(3)在突變點修正濾波增益Kt，通過大的濾波增益Kt可以對狀態(tài)估計進行較強的修正。

5.1期望最大化參數(shù)估計

當系統(tǒng)在運行過程中，隨著外界環(huán)境的變化，模型參數(shù)可能發(fā)生變化，因此可以利用監(jiān)測數(shù)據(jù)對模型參數(shù)進行估計。設(shè)式(4)和式(5)中的未知參數(shù)可用一個參數(shù)向量θ表示，即θ={Ftt-1,Ht,Qt,Rt}，在可實時監(jiān)測條件下，需進行狀態(tài)空間模型參數(shù)的估計及退化狀態(tài)的預(yù)測。

期望最大化(Expectation-Maximization，EM)算法主要用來計算基于不完全數(shù)據(jù)的極大似然估計，尤其適用于估計存在隱藏變量下的狀態(tài)空間模型參數(shù)估計。因此首先采用期望最大化算法從初始時刻以迭代的方式實時估計狀態(tài)空間模型參數(shù)。該算法包括兩個步驟：E步驟和M步驟。

(1)E步驟

θ0={F0,H0,Q0,R0}

(9)

為使退化狀態(tài)初值選擇更為合理，本文采用卡爾曼前向濾波及固定區(qū)間反向濾波相結(jié)合的方法在參數(shù)估計的同時來選取最優(yōu)的退化狀態(tài)初值。

已知在ti+n時刻，接收到i+n個監(jiān)測點的監(jiān)測信息，求得其對數(shù)最大似然函數(shù)為：

(10)

(11)

采用期望最大化算法估計參數(shù)，需首先計算式的期望值，記為EXY,θk[lnL(θ)]。

(12)

卡爾曼光滑算法即為求期望值的過程[21]，首先求得：

(13)

(14)

(15)

(16)

(17)

(18)

將式(13)～式(18)代入式(12)，可求得：

(19)

(2)M步驟

通過求解最大似然估計算法求得各參數(shù)的估計值：

(20)

EM算法不斷進行E步驟和M步驟的迭代，直到似然函數(shù)值穩(wěn)定不再變化，即求解到ti+n時刻的參數(shù)估計值θti+n。

5.2突變點后初始值的修正

θti+n={Fti+nti+n-1,Hti+n,Qti+n,Rti+n}

(21)

如圖6所示，當在ti+n時刻檢測到ti時刻存在突變點A，需修正突變點A之后的卡爾曼濾波模型初始值及參數(shù)。

利用ti時刻及之前接收到的所有監(jiān)測信息采用EM參數(shù)估計算法求得的參數(shù)估計值記為：

θti={Ftiti-1,Hti,Qti,Rti}

(22)

(23)

式中：β與γ為參數(shù)修正因子，且β+γ=1。選擇γ<β，通過參數(shù)修正因子來加大突變點之后的數(shù)據(jù)的權(quán)重，相應(yīng)減少突變點之前的數(shù)據(jù)的權(quán)重，重新選擇突變點之后參數(shù)估計初值。

(24)

5.3濾波增益的修正

檢測到突變點后，通過自適應(yīng)增益函數(shù)增加突變點后的增益矩陣Kt的模值，使得狀態(tài)估計比較快的收斂于部件退化的真實狀態(tài)。假設(shè)檢測到ti處退化狀態(tài)發(fā)生突變，構(gòu)造一自適應(yīng)增益函數(shù)gi(ti)為：

(25)

式中：a與b為自適應(yīng)增益函數(shù)的參數(shù)，增大濾波增益的方法是以犧牲濾波的最優(yōu)性為代價的，因此并非濾波增益收斂到1的速度越大越好，因此自適應(yīng)增益函數(shù)的參數(shù)a與b的選擇應(yīng)同時考慮其收斂速度與濾波最優(yōu)性的問題。

(26)

5.4退化狀態(tài)突變后自適應(yīng)濾波

(27)

(28)

(29)

(30)

(31)

(32)

式中ti時刻需修正卡爾曼濾波增益，修正后的卡爾曼濾波增益記為：

(33)

6 試驗結(jié)果

6.1振動信號的特征提取

齒輪箱中的齒輪發(fā)生故障時，安裝在軸承位置的傳感器經(jīng)過軸和軸承接收到振動信號，信號衰減最小，可以很好地反映齒輪箱的振動特性。因此在進行壽命預(yù)測時選擇齒輪疲勞試驗中離斷齒位置最近，且在軸承座位置布置的4#傳感器輸出的463組振動信號進行特征提取。試驗中采樣頻率25.6 kHz，采樣時間60 s，采樣間隔9 min，將采樣點次數(shù)折算為監(jiān)測時間，就可得到如圖7所示的4#傳感器輸出的特征值隨監(jiān)測時間變化曲線圖。該曲線包含了從開始磨合階段到完成疲勞試驗共77 h的振動信號均方幅值。

整個RMS曲線在去除個別奇異點后整體變化趨勢可反映試驗齒輪各監(jiān)測點磨損情況與振動能量對應(yīng)的關(guān)系。整個RMS曲線在開始嚙和階段逐漸趨于穩(wěn)定，在60 h后齒輪磨損增大，幅值顯著上升直到發(fā)生斷齒。

圖7 4#傳感器特征值提取RMS曲線

(1) 當監(jiān)測時間t∈[0,10]時，齒輪處于嚙合階段，RMS呈下降趨勢，表明齒輪在開始嚙合階段均方幅值比正常疲勞磨損要大。

(2)t∈[10,68]時，RMS值在逐漸增加，表明齒輪嚙合后齒面進入正常磨損階段。由于電壓值不穩(wěn)定，導(dǎo)致轉(zhuǎn)速發(fā)生變化在大約30 h處和55 h處RMS突然增大，然后RMS又趨于穩(wěn)定。

(3)t∈[68,77]時，RMS急劇增大，從最后齒輪疲勞失效后結(jié)果得知，此時的齒輪磨損開始加劇，直到在77.2 h處發(fā)生斷齒，此時退化狀態(tài)故障閾值為xf=77.375 mm/s2。

6.2狀態(tài)空間建模及突變狀態(tài)檢測

6.2.1 狀態(tài)空間建模及初始參數(shù)的選擇

根據(jù)式(13)和(14)建立齒輪疲勞試驗中4#傳感器輸出的觀測值與系統(tǒng)狀態(tài)之間狀態(tài)空間模型為：

x(t)=Ftt-1x(t-1)+w(t-1)

(34)

y(t)=Htx(t)+v(t)

(35)

6.2.2 退化狀態(tài)突變點的檢測

系統(tǒng)開始運行經(jīng)過齒輪嚙合階段進入齒輪正常退化階段后進行突變點檢測，選擇誤警率η=0.25%，平均運行鏈長ARL=400，取最小的預(yù)警退化量為Δμ=20，退化均值μ0=51.952，方差σ2=3.044，突變預(yù)警閾值由式和式可求解突變檢測預(yù)警閾值h=59.702 3，得到累積和統(tǒng)計量如圖8所示。

圖8 報警時CUSUM統(tǒng)計量變化圖

圖8中A點為突變點，B點為預(yù)警點。檢測到退化狀態(tài)的突變點發(fā)生在在66.3 h處，在69 h處累積和統(tǒng)計量超過預(yù)警閾值h進行報警，由于齒輪磨損退化狀態(tài)發(fā)生變化。

6.3預(yù)測模型的修正

根據(jù)檢測到的突變點對壽命預(yù)測的模型進行修正可以更好的跟蹤系統(tǒng)的動態(tài)變化，為剩余壽命的實時預(yù)測及維護維修策略的制定提供依據(jù)。因此利用檢測出的66.3 h處突變點提供的退化信息對預(yù)測模型進行修正，在69 h處開始預(yù)測其退化值與實際接收到的齒輪退化值進行比較，計算其預(yù)測誤差。

6.3.1 突變點后初始值的修正

圖9 預(yù)測模型修正前后所估計的參數(shù)F的取值

Fig.9 Values of the estimated parameters for modified prediction model comparison with original

由圖9可知突變點前，狀態(tài)空間模型修正前后的模型參數(shù)F取值變化較接近，在突變點處，根據(jù)突變點信息對模型進行修正后的參數(shù)變化較大，而模型修正前的參數(shù)變化相對滯后。

6.3.2 濾波增益參數(shù)的選擇

由圖10可知隨著ti的增大，b為定值時，增大參數(shù)a的值，或在a為定值時，減小b的取值，gi(ti)收斂到1的速度在增大，(增大濾波增益的方法是以犧牲濾波的最優(yōu)性為代價的，因此并非gi(ti)收斂到1的速度越大越好，因此自適應(yīng)增益函數(shù)的參數(shù)a與b的選擇應(yīng)同時考慮其收斂速度與濾波最優(yōu)性的問題。本文中濾波增益參數(shù)取值為a=8×10-2，b=2時，則自適應(yīng)增益函數(shù)gi(ti)，如圖11所示。

(a) kti=1時，改變參數(shù)a的自適應(yīng)增益函數(shù)取值

(b) kti=1時，改變參數(shù)b的自適應(yīng)增益函數(shù)取值

圖11 參數(shù)為a=0.08，b=2的自適應(yīng)增益函數(shù)

6.4剩余壽命預(yù)測

利用修正后的狀態(tài)空間方程可對系統(tǒng)的監(jiān)測點之后的退化狀態(tài)進行預(yù)測，通過預(yù)測到的退化狀態(tài)值及已知的退化狀態(tài)故障閾值可求解首次到達故障閾值的時間(First Passage Time, FPT)[9]。進而以蒙特卡洛的分析方法求得在tp時刻所預(yù)測的系統(tǒng)的平均剩余壽命uN(tp)為：

(36)

式中：Tfm為第m次預(yù)測退化狀態(tài)到達利用齒輪疲勞試驗所得系統(tǒng)故障閾值xf=77.375 mm/s2的時間，N為蒙特卡洛仿真次數(shù)，用最小方差控制N的取值，即N滿足：

(37)

試驗中ε的取值為：ε=e-10。

圖12給出在69 h處狀態(tài)空間模型修正前后通過預(yù)測得到齒輪磨損退化狀態(tài)均值。

(a) 狀態(tài)空間模型未修正時預(yù)測退化狀態(tài)

(b) 狀態(tài)空間模型修正后預(yù)測退化狀態(tài)

如圖13所示可知在不同時刻修正前后的狀態(tài)空間模型所預(yù)測的退化狀態(tài)曲線均不同，69 h處預(yù)測狀態(tài)值偏離實際值較多，誤差較大。隨著監(jiān)測信息增多，濾波時間充分長后，修正前后的狀態(tài)空間模型所預(yù)測的退化狀態(tài)均接近真實值，誤差減小。在相同的預(yù)測時刻，比較修正前后的狀態(tài)空間模型預(yù)測結(jié)果，修正后的狀態(tài)空間模型所預(yù)測的齒輪退化狀態(tài)更接近真實值，能很好地預(yù)測系統(tǒng)退化狀態(tài)發(fā)生突變后其退化狀態(tài)的變化以及故障發(fā)生時間。未修正的狀態(tài)空間模型所預(yù)測到的特征不能很好地突變后的退化狀態(tài)變化，預(yù)測剩余壽命誤差較大。由圖13(e)可以看到在剛檢測到突變狀態(tài)點時，由于修正后的狀態(tài)空間模型剛剛修正初值，開始的遞推計算所得到的狀態(tài)估計值比較粗略，估計誤差較大，但經(jīng)過幾步遞推運算之后狀態(tài)估計值比較快的靠近被估計的真實狀態(tài)，能夠很好地跟蹤其動態(tài)變化。

(a) 69 h處模型修正前后預(yù)測

(c) 73 h處模型修正前后預(yù)測

(d) 75 h處模型修正前后預(yù)測

(e) 76 h狀態(tài)空間模型修正前后處

(f) 73 h處模型修正前后預(yù)測

由蒙特卡羅法求得狀態(tài)空間模型修正前后不同時刻所預(yù)測的剩余壽命概率密度函數(shù)如圖14所示?？芍谕活A(yù)測時刻，狀態(tài)空間模型修正后所預(yù)測的平均剩余壽命更接近真實的剩余壽命，且其概率密度函數(shù)與修正前的狀態(tài)空間模型所預(yù)測的剩余壽命概率密度函數(shù)相比預(yù)測方差減小。隨著預(yù)測時間的推移，接收到監(jiān)測信息的增加，在76 h利用狀態(tài)空間模型所預(yù)測到的剩余壽命概率密度函數(shù)，模型修正前后剩余壽命概率密度函數(shù)峰值均逐漸接近真實值。

為進一步比較狀態(tài)空間模型修正前后預(yù)測剩余壽命的準確性，引入預(yù)測準確度的概念[23]。預(yù)測準確度PA為：

(38)

式中：‖表示絕對值；uN(tp)表示tp時刻預(yù)測求得的平均剩余壽命；Ta為齒輪實際剩余壽命；T*=77.2為齒輪實際故障時間，則Ta=T*-tp。由此可計算狀態(tài)空間模型修正前后其預(yù)測準確度的比較如表1所示。

如表1所示，隨著監(jiān)測信息的增多，修正前后的狀態(tài)空間模型所預(yù)測剩余壽命預(yù)測誤差均是逐漸減小，預(yù)測準確度逐漸提高。但在不同的預(yù)測點69 h、70 h以及73 h和75 h處相比，狀態(tài)空間模型修正后所預(yù)測的剩余壽命準確度均高于利用修正前的狀態(tài)空間模型所預(yù)測求得的預(yù)測準確度。在76 h處未修正預(yù)測模型剩余壽命逐漸接近模型修正后所預(yù)測到的剩余壽命，為進一步說明考慮突變狀態(tài)檢測的齒輪實時剩余壽命預(yù)測有效性，建立基于反向傳播(Back Propagation，BP)神經(jīng)網(wǎng)絡(luò)齒輪壽命預(yù)測模型并與本文所提預(yù)測方法進行比較。首先將齒輪疲勞試驗中突變點之后的退化狀態(tài)特征值作為網(wǎng)絡(luò)訓(xùn)練輸入向量；然后建立三層BP神經(jīng)網(wǎng)絡(luò)并對網(wǎng)絡(luò)進行訓(xùn)練；最后利用訓(xùn)練好的網(wǎng)絡(luò)得到當前時刻退化狀態(tài)的預(yù)測曲線。

仍以齒輪疲勞壽命試驗4#傳感器輸出的463組振動信號進行特征提取，通過仿真算例將預(yù)測值和實際值相比較(如圖15所示)可知：在不同的預(yù)測點70 h、 73 h、75 h和76 h處相比，隨著訓(xùn)練數(shù)據(jù)增加，神經(jīng)網(wǎng)絡(luò)預(yù)測退化狀態(tài)所預(yù)測剩余壽命預(yù)測誤差逐漸減小。如圖15(d)所示在76 h處通過神經(jīng)網(wǎng)絡(luò)預(yù)測模型所預(yù)測的齒輪故障時間為80.1 h，基于狀態(tài)空間模型突變點修正后所預(yù)測的齒輪故障時間為77 h，齒輪實際故障時間T*=77.2 h?？芍谙嗤念A(yù)測點，基于狀態(tài)空間模型突變點修正后所預(yù)測的剩余壽命準確度高于利用BP神經(jīng)網(wǎng)絡(luò)所預(yù)測求得的預(yù)測準確度。

(a) 70 h修正前剩余壽命概率密度函數(shù)

(b) 70 h修正后剩余壽命概率密度函數(shù)

(c) 73 h修正前剩余壽命概率密度函數(shù)

(d) 73 h修正后剩余壽命概率密度函數(shù)

(e)75 h修正前剩余壽命概率密度函數(shù)

(f)75 h修正后剩余壽命概率密度函數(shù)

(g) 76 h修正前剩余壽命概率密度函數(shù)

(h) 76 h修正后剩余壽命概率密度函數(shù)

預(yù)測時刻tp/h修正前uN(tp)/h修正后uN(tp)/h實際剩余壽命Ta/h修正前預(yù)測準確度PA/%修正后預(yù)測準確度PA′/%6914．15．78．228．9569．517012．35．07．229．1769．44736．54．04．245．2495．24751．32．32．259．0995．45761．01．01．295．2495．24

由此可知當接收到的監(jiān)測信息足夠多時，利用數(shù)據(jù)驅(qū)動的剩余壽命預(yù)測方法可以較準確地預(yù)測剩余壽命。但是基于狀態(tài)的預(yù)防性維護維修更關(guān)心系統(tǒng)整個運行過程中的實時剩余壽命預(yù)測的準確度，尤其系統(tǒng)退化性能發(fā)生突變后準確的剩余壽命預(yù)測可以為及時調(diào)整預(yù)防性維護維修策略提供必要的依據(jù)。因此本文所提基于狀態(tài)空間的齒輪實時剩余壽命預(yù)測方法，在系統(tǒng)運行過程中發(fā)生突變后，利用突變點所提供的信息對模型進行修正可以更好地跟蹤系統(tǒng)的退化狀態(tài)的變化，提高系統(tǒng)運行過程中實時剩余壽命的預(yù)測準確度。

7 結(jié) 論

本文根據(jù)接收到的齒輪實時監(jiān)測信息建立系統(tǒng)狀態(tài)空間預(yù)測模型。同時考慮異常點檢測與剩余壽命預(yù)測的關(guān)聯(lián)性，對齒輪磨損退化過程中的突變狀態(tài)點進行檢測后，利用退化狀態(tài)突變點所提供的壽命信息采用卡爾曼前向濾波及平滑算法結(jié)合EM參數(shù)估計算法在濾波的同時修正狀態(tài)空間模型，進行狀態(tài)估計及壽命預(yù)測。將退化狀態(tài)突變點檢測和剩余壽命的預(yù)測進行聯(lián)合研究，通過齒輪疲勞壽命試驗驗證其修正后的預(yù)測模型可以更快的對系統(tǒng)的動態(tài)變化進行跟蹤，提高實時剩余壽命預(yù)測準確度。

(a) 70 h處神經(jīng)網(wǎng)絡(luò)預(yù)測退化狀態(tài)

(b) 73 h處神經(jīng)網(wǎng)絡(luò)預(yù)測退化狀態(tài)

(c) 75 h處神經(jīng)網(wǎng)絡(luò)預(yù)測退化狀態(tài)

(d) 76 h處神經(jīng)網(wǎng)絡(luò)預(yù)測退化狀態(tài)

作者的下一步工作是針對齒輪利用多個傳感器接收到的數(shù)據(jù)，考慮數(shù)據(jù)融合并對預(yù)測模型修正，進行更為準確的實時壽命預(yù)測。

[1] 馮志鵬, 秦嗣峰. 基于 Hilbert 振動分解和高階能量算子的行星齒輪箱故障診斷研究[J]. 振動與沖擊, 2016, 35(5): 47-54.

FENG Zhipeng,QIN Sifeng. Planetary gearbox fault diagnosis based on Hilbert vibration decomposition and higher order differential energy operator[J]. Journal of Vibration and Shock, 2016, 35(5): 47-54.

[2] SI X, WANG W, CHEN M, et al. A degradation path-dependent approach for remaining useful life estimation with an exact and closed-form solution[J]. European Journal of Operational Research, 2013, 226(1): 53-66.

[4] 孫強, 岳繼光. 基于不確定性的故障預(yù)測方法綜述[J]. 控制與決策, 2014, 29(5): 769-778.

SUN Qiang, YUE Jiguang. Review on fault prognostic methods based on uncertainty[J]. Control and Decision, 2014, 29(5): 769-778.

[6] CARR M J, WANG W. An approximate algorithm for prognostic modelling using condition monitoring information[J]. European Journal of Operational Research, 2011, 211(1): 90-96.

[7] SUN J, ZUO H, WANG W, et al. Application of a state space modeling technique to system prognostics based on a health index for condition-based maintenance[J]. Mechanical Systems and Signal Processing, 2012, 28: 585-596.

[8] BRESSEL M, HILAIRET M, HISSEL D. Extended kalman filter for prognostic of proton exchange membrane fuel cell[J]. Applied Energy, 2016, 164: 220-227.

[9] TANG D, MAKIS V, JAFARI L, et al. Optimal maintenance policy and residual life estimation for a slowly degrading system subject to condition monitoring[J]. Reliability Engineering & System Safety, 2015, 134: 198-207.

[10] LALL P, LOWE R, GOEBEL K. Extended kalman filter models and resistance spectroscopy for prognostication and health monitoring of leadfree electronics under vibration[J]. Reliability, IEEE Transactions on, 2012, 61(4): 858-871.

[11] WANG W. A model to predict the residual life of rolling element bearings given monitored condition information to date[J]. IMA Journal of Management Mathematics, 2002, 13(1): 3-16.

[12] 胡昌華, 王志遠, 周志杰. 基于隨機濾波理論的剩余壽命預(yù)測方法研究[J]. 系統(tǒng)仿真技術(shù), 2011, 7(2): 83-88.

HU Changhua, WANG Zhiyuan, ZHOU Zhijie. A study on residual life prediction method based on stochastic filtering theory[J]. System Simulation Technology, 2011, 7(2):83-88.

[13] 李鑫, 呂琛, 王自力, 等. 考慮退化模式動態(tài)轉(zhuǎn)移的健康狀態(tài)自適應(yīng)預(yù)測[J]. 自動化學(xué)報, 2014, 40(9): 1889-1895.

LI Xin, Lü Chen, WANG Zili, et al. Self-adaptive health condition prediction considering dynamic transfer of degradation mode[J]. Acta Automatica Sinica, 2014, 40(9):1889-1895.

[14] 周東華, 魏慕恒, 司小勝. 工業(yè)過程異常檢測, 壽命預(yù)測與維修決策的研究進展[J]. 自動化學(xué)報, 2013, 39(6): 711-722.

ZHOU Donghua, WEI Muheng, SI Xiaosheng. A survey on anomaly detection, life prediction and maintenance decision for industrial processes[J]. Acta Automatica Sinica, 2013, 39(6): 711-722.

[15] SON J, ZHANG Y, SANKAVARAM C, et al. RUL prediction for individual units based on condition monitoring signals with a change point[J]. Reliability, IEEE Transactions on, 2015, 64(1): 182-196.

[16] 劉亞瓊, 王鐵, 武志斐, 等. QC商用車變速器疲勞試驗振動分析[J]. 機械傳動, 2014, 38(7): 104-105.

LIU Yaqiong, WANG Tie, WU Zhifei, et al. Vibration analysis of QC commercial vehicle transmission fatigue test[J]. Mechanical Drive, 2014, 7(38): 104-105.

[17] HINES J A, MARK M D. Bending-fatigue damage-detection on notched-tooth spiral-bevel gears using the average-log-ratio, ALR, algorithm[J]. Mechanical Systems and Signal Processing, 2014, 43(1): 44-56.

[18] UNNIKRISHNAN J A, VEERAVALLI V V, MEYN S P. Minimax robust quickest change detection[J]. Information Theory, IEEE Transactions on, 2011, 57(3): 1604-1614.

[19] 伊廷華, 郭慶, 李宏男. 基于控制圖的 GPS 異常監(jiān)測數(shù)據(jù)檢驗方法研究[J]. 工程力學(xué), 2013, 30(8): 133-141.

YI Tinghua, GUO Qing, LI Hongnan. The research on detection methods of gps abnormal monitoring data based on control chart[J].Engineering Mechanics, 2013,30(8):133-141.

[20] SIEGMUND D. Error probabilities and average sample number of the sequential probability ratio test[J]. Journal of the Royal Statistical Society. Series B (Methodological), 1975: 394-401.

[21] WANG W, CARR M, XU W, et al. A model for residual life prediction based on Brownian motion with an adaptive drift[J]. Microelectronics Reliability, 2011, 51(2): 285-293.

[22] KHAN M E, DUTT D N. An expectation-maximization algorithm based kalman smoother approach for event-related desynchronization (ERD) estimation from EEG[J]. Biomedical Engineering, IEEE Transactions on, 2007, 54(7): 1191-1198.

[23] LU C, TAO L, FAN H. An intelligent approach to machine component health prognostics by utilizing only truncated histories[J]. Mechanical Systems and Signal Processing, 2014, 42(1): 300-313.

Modelfortherealtimeremainingusefullifepredictionofgearsbasedontheabruptchangedetection

SHIHui1,ZENGJianchao1,2

(1.Division of Industrial and System Engineering, Taiyuan University of Science and Technology, Taiyuan 030024，China; 2. School of Computer Science and Control Engineering, North University of China, Taiyuan 030051, China)

In order to accurately predict the gear remaining useful life in the degradation process, a new method for the real-time prediction of gear contact fatigue remaining useful life was put forward, which is a method integrating the abrupt change detection and remaining life prediction. A state-space model for predicting degradation states of gear wear was established by using the real time monitoring vibration information to update the model parameters. The Kalman forward filtering and smoothing algorithm combined with the parameter estimation by the expectation-maximization algorithm was made in use and the prediction model was modified incessantly to change the filtering effect according to the life information from the abrupt change detection. The real-time monitoring data of the contact fatigue life collected on a gear test rig were used to verify the model proposed. The results show that a revised prediction model using the abrupt point information can achieve the prediction faster than a dynamic tracking system, and can improve the accuracy of gear degradation state and real-time remaining useful life prediction.

remaining useful life prediction; state-space models; Kalman filtering; abrupt change detection; model correction

TH163.5；TP277

A

10.13465/j.cnki.jvs.2017.21.026

2016-01-11 修改稿收到日期：2016-08-05

石慧 女，博士，副教授，1979年生

曾建潮 男，教授，博士生導(dǎo)師，1963年生