考慮突變狀態檢測的齒輪實時剩余壽命預測

石 慧， 曾建潮,2

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

考慮突變狀態檢測的齒輪實時剩余壽命預測

石 慧1， 曾建潮1,2

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

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

剩余壽命預測；狀態空間建模；卡爾曼濾波；突變狀態檢測；模型修正

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

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

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

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

1 齒輪疲勞試驗

1.1疲勞試驗臺架

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

圖1 齒輪試驗臺架

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

(a) 齒輪試驗臺架

(b) 主試箱傳感器

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

圖3 主試箱齒輪發生故障斷齒

1.2試驗齒輪的安裝

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

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

2 特征提取

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

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

(1)

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

3 基于狀態空間模型的退化狀態建模

在壽命預測中，如何建立一個適當的模型是首要問題。可以通過非線性狀態空間建模進行狀態估計。模型的一般形式為：

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

(2)

Yt=h(Xt)+vt

(3)

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

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

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

(4)

Yt=HtXt+vt

(5)

4 突變狀態檢測

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

圖5 退化狀態偏移的突變狀態點檢測示意圖

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

max(yn-μ0-k)

(6)

(7)

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

-2k(h/σ+1.166)

(8)

5 退化狀態預測模型修正

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

圖6 退化狀態預測模型修正示意圖

由以上分析可知采用卡爾曼濾波對系統進行狀態估計及預測時，當系統出現狀態突變點后，仍然存在以下問題：

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

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

因此，采用狀態空間模型進行狀態估計與預測時，在突變狀態之后對狀態空間模型進行修正，期望觀測信息對狀態估計值產生較強的修正作用，使得狀態估計較快地收斂于部件退化的真實狀態?？梢愿玫膶ν蛔儬顟B進行跟蹤，減少預測誤差，提高收斂速度和預測精度。需對模型進行以下修正：

(1)狀態估計時增加突變點后的監測信息，修正系統模型參數，進行自適應的狀態估計。

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

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

5.1期望最大化參數估計

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

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

(1)E步驟

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

(9)

為使退化狀態初值選擇更為合理，本文采用卡爾曼前向濾波及固定區間反向濾波相結合的方法在參數估計的同時來選取最優的退化狀態初值。

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

(10)

(11)

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

(12)

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

(13)

(14)

(15)

(16)

(17)

(18)

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

(19)

(2)M步驟

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

(20)

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

5.2突變點后初始值的修正

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

(21)

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

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

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

(22)

(23)

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

(24)

5.3濾波增益的修正

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

(25)

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

(26)

5.4退化狀態突變后自適應濾波

(27)

(28)

(29)

(30)

(31)

(32)

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

(33)

6 試驗結果

6.1振動信號的特征提取

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

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

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

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

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

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

6.2狀態空間建模及突變狀態檢測

6.2.1 狀態空間建模及初始參數的選擇

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

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

(34)

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

(35)

6.2.2 退化狀態突變點的檢測

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

圖8 報警時CUSUM統計量變化圖

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

6.3預測模型的修正

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

6.3.1 突變點后初始值的修正

圖9 預測模型修正前后所估計的參數F的取值

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

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

6.3.2 濾波增益參數的選擇

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

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

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

圖11 參數為a=0.08，b=2的自適應增益函數

6.4剩余壽命預測

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

(36)

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

(37)

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

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

(a) 狀態空間模型未修正時預測退化狀態

(b) 狀態空間模型修正后預測退化狀態

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

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

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

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

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

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

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

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

(38)

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

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

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

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

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

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

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

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

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

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

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

預測時刻tp/h修正前uN(tp)/h修正后uN(tp)/h實際剩余壽命Ta/h修正前預測準確度PA/%修正后預測準確度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

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

7 結 論

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

(a) 70 h處神經網絡預測退化狀態

(b) 73 h處神經網絡預測退化狀態

(c) 75 h處神經網絡預測退化狀態

(d) 76 h處神經網絡預測退化狀態

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

[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] 孫強, 岳繼光. 基于不確定性的故障預測方法綜述[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] 胡昌華, 王志遠, 周志杰. 基于隨機濾波理論的剩余壽命預測方法研究[J]. 系統仿真技術, 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] 李鑫, 呂琛, 王自力, 等. 考慮退化模式動態轉移的健康狀態自適應預測[J]. 自動化學報, 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] 周東華, 魏慕恒, 司小勝. 工業過程異常檢測, 壽命預測與維修決策的研究進展[J]. 自動化學報, 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 異常監測數據檢驗方法研究[J]. 工程力學, 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年生

曾建潮 男，教授，博士生導師，1963年生