基于Gamma退化過程的剩余壽命預測及維修決策優化模型研究

谷玉波，賈云獻，張英波

(軍械工程學院 裝備指揮與管理系，石家莊 050003)

機械設備在使用過程中，某些部件的運行狀態會逐漸退化。狀態的退化可以表現為機械部件的磨損、裂紋的增長、腐蝕程度的加深等，這些都是一系列理化效應的結果。產品狀態退化的最終結果會導致功能故障的發生，而且在產品的退化過程中能夠直接表征退化狀態的狀態指標(磨損、裂紋、腐蝕等)經常是不可測量或難于直接進行測量的。為獲得設備狀態退化程度的數據，可以運用CBM狀態監控技術對設備的狀態數據進行提取分析，間接地反映出設備的健康狀態，實現對設備剩余壽命的預測，最終實現維修決策的優化，使維修效果達到最佳。在此，運用Gamma退化過程狀態空間模型來描述設備狀態性能的退化，通過該模型實現了剩余壽命的預測，并運用剩余壽命預測結果對維修決策模型進行優化。

1 Gamma退化過程的狀態空間模型

在常用的分布中，Gamma分布可以用來描述連續累積磨損過程，并且具有退化過程建模所需的所有屬性，即非負、平穩增長和獨立增量過程[1]。因此，運用Gamma分布建立狀態空間模型，可以描述設備的退化過程，并實現剩余壽命的預測。

1.1 狀態空間模型

狀態空間模型(State Space Model，SSM)是一種用來描述系統演化過程的模型，可以用來建立產品的狀態退化模型。狀態空間模型假設系統狀態隨時間的演化可以由一個不可觀測或難以觀測的狀態序列確定，與狀態序列伴隨的是一個可以觀測的測量序列，兩者之間的關系可以由狀態空間模型來描述，其中應用較為普遍的狀態空間模型是典型相關(Canonical Correlation)方法[2]。通過建立狀態方程和觀測方程，狀態空間模型為描述系統的動態演化特性建立了一個有效的模型框架，一些復雜的問題就可以用比較簡單的形式表示。狀態空間模型的基本方程為

xi+1=Fi(xi,ui,ζ,θ1)，

(1)

yi=Hi(xi,ui,ε,θ2)，

(2)

式中：Fi為狀態轉移函數，Hi為測量函數，一般而言Fi和Hi假設是已知的；ui為系統的輸入變量或控制變量；θ1,θ2分別為狀態方程和觀測方程的靜態參數；ζ和ε分別為系統噪聲和觀測噪聲，一般假設ζ和ε相互獨立，且服從某一特定分布。

1.2 Gamma過程

Gamma過程是一個具有獨立、非負增量的隨機過程。運用Gamma過程對一些退化過程建立模型，能夠使建模過程中的數學計算更加明確易懂[3]，適合描述在時域上[4]單調累積微小增量的退化過程，比如磨損、腐蝕、裂紋增長等。

根據Gamma過程，定義非負隨機過程{X(t);t≥0}，其概率密度函數為

(3)

2 Gamma-SSM退化模型

2.1 模型的建立

假設設備狀態x的退化符合Gamma過程，y表示與x相對應的觀測變量，xi表示設備ti時刻所對應狀態，則建立的狀態空間模型的狀態方程和觀測方程為

xi-xi-1～Gamma(α(ti)-α(ti-1),λ)，

(4)

yi=Hi(xi)+ε。

(5)

假設當系統狀態x達到故障閾值Xf時系統發生故障，如圖1所示。則系統狀態從0首次達到故障閾值所需的時間為

圖1 系統退化過程

Tf=inf{t:x=Xf,t>0}。

(6)

為簡化問題，假設觀測噪聲是可以加的。假設設備狀態退化過程為平穩Gamma過程，形狀參數α(t)為線性函數，即α(t)=a·t，并且假設觀測量y與狀態量x之間的關系為y=c·x+ε，則在獲得觀測序列y0:n的條件下，退化模型可以表示為

xi+1-xi～Gamma(α·(ti+1-ti),λ)，

(7)

yi=c·xi+ε。

(8)

由假設可知，觀測噪聲ε服從均值為0的正態分布。所以，y～N(c·x,σ)。因此只要獲得參數a,λ,c,σ的值，就可以確定該退化模型的形式。

2.2 模型參數估計

在對上述模型中參數進行求解的過程中，會出現新的需要估計的參數，導致運算量增大。為解決此問題，運用文獻[5]提到的經驗最大化(Experience Maximization，EM)算法與粒子濾波(Particle Filtering，PF)算法相結合的方法求解模型中出現的參數。

EM算法是進行極大似然估計的一種有效方法，簡單來說，這一算法的每一次迭代主要包括2個步驟：(1)求期望(Expectation Step)，稱為E步；(2)求極大值(Maximization)，稱為M步。EM算法通過假設潛在變量的存在，來簡化似然函數極大估計[3]。

PF算法是一種基于Monte Carlo方法和遞推Bayes估計的統計濾波方法[6]，建立在序貫重要性采樣(Sequence Importance Sampling，SIS)和Bayes理論的基礎上，對于解決非線性、非Gauss問題有很好的效果[7]。PF算法的基本思想是：首先根據系統狀態向量的經驗分布在狀態空間產生一組隨機樣本集合(即粒子)，然后根據觀測數據不斷調整粒子的權重和位置，通過調整后的粒子信息修正最初的經驗分布。其實質是用由粒子及權重組成的離散隨機測度近似相關的概率分布，并根據算法遞推更新離散隨機測度。當粒子容量足夠大時，就近似于狀態變量真實的后驗概率分布。

由于EM-PF算法兼有兩種算法各自的優點，可以更加有效地對模型中的參數進行求解。模型參數的求解過程如圖2所示。

圖2 模型參數求解過程

3 剩余壽命預測

剩余壽命是指產品從被檢測的某一時刻起到該產品發生故障的時間長度[8]，設備剩余壽命的預測是維修管理中的一個重要環節，由于剩余壽命預測對產品使用過程中的安全性、經濟性和任務性具有重要的影響，因此它是制定維修決策的重要依據。

通過確定狀態空間模型參數，就能得到Gamma狀態空間模型的具體形式，進而得到產品的剩余壽命分布函數。令Ns為粒子個數；Tf為設備故障時間；wi為第i個粒子的權重；Xf為故障閾值；y0:c為直至當前時刻所獲得的觀測序列；y0:i為當前獲得的觀測值序列；ti為第i次狀態檢測時間；xi為ti時刻的狀態值；x1:n為n次檢測所獲得的狀態序列，即x1:n={x1,x1,x2,…,xn} ；yi為ti時刻的退化量觀測值；y0:n為n次檢測所獲得的觀測序列，即y0:n={y0,y1,y2,…,yn},則ti時刻的剩余壽命τi分布函數為

(9)

p(τi+ti≥Tf|xi)=p(x(τi+ti)≥

(10)

(11)

另外，預測tk時刻的狀態概率密度函數可表示為

(12)

(13)

(14)

(15)

4 維修決策優化模型

構建剩余壽命概率密度函數的最終目的是根據剩余壽命與實際需求做出維修決策，確定合理的維修或更換時間。平時對設備進行維修管理的過程中，人們總是希望設備的維修管理費用盡可能達到最低。以費用最小為目標的維修決策模型可以根據設備當前的運行狀態確定在下一間隔期內是否進行維修以及何時進行維修，可最大程度地避免故障的發生。

在維修策略中，產品有2種維修方式：(1)修復性維修，即產品在達到預定維修時刻之前發生故障而進行的維修或更換活動，也稱故障維修或故障更換；(2)預防性維修，即產品在達到預定的維修時刻進行的維修或更換活動。因此按照維修決策的目標，并根據剩余壽命概率密度函數，建立以最小費用為目標的維修決策模型，其決策過程如圖3所示。

圖3 以費用最小為目標的維修決策示意圖

4.1 假設和參數

(1)每隔Δt進行一次狀態檢測并獲取當前時刻的狀態數據。

(2)如果在檢測間隔期內發生故障，則立刻進行修復性維修。

(3)設備的預防性維修費用小于故障后的修復性維修費用，即Tp

(4)預防性維修和修復性維修可以使設備恢復初始狀態，即修復如新。

(5)參數：E(C)為更新周期內的期望總費用；E(T)為期望更新周期長度；Δt為狀態檢測間隔期，Δt=ti-ti-1；f(τi|y0:i)為ti時刻狀態信息為y0:i時產品的剩余壽命概率密度分布函數；Pp為部件進行預防性維修的概率；Pf為部件進行修復性維修的概率；cp為預防性維修費用；cf為修復性維修費用；TR為最佳維修更換時間；C(TR)為部件維修更換時間為TR時，長期使用下的單位時間費用；ti為當前時刻的狀態監控點，t0=0；Tp為預防性維修所需的平均時間；Tc為修復性維修所需的平均時間；ci為每次實施狀態檢測的費用。

4.2 以費用為目標的決策模型

設備部件在長期使用下的單位時間費用可表示為

(16)

E(C)=cfPf+cpPp+nci=cff(τi

ti|y0:i)+cp[1-f(τi

cp+(cf-cp)f(τi

(17)

(18)

(19)

將(17)～(18)式代入(16)式中，則費用模型為

(20)

其中，nci為更新周期內監控所產生的費用，由于狀態檢測中檢測費用大多為設備投入費用，nci常常可以忽略不計，因此在實際決策中，近似取nci=0。如果Tp和Tc相對于TR來說很小，在費用模型中也可以將其忽略。

利用數值方法最小化(20)式，得出設備進行基于狀態的維修時可達到的最小單位費用及所對應的維修時刻TR。當獲取新的狀態信息值y0:i+1時，需要將數據代入(20)式重新計算更新結果。

4.3 最佳維修時間確定方法

通過費用決策模型可以得到各檢測時刻對應的維修時間TR，使得單位時間維修費用最小。當TR-ti>Δt時，不對部件進行維修，直到下次正常檢測時刻；當首次出現TR-ti≤Δt時，則到達TR時對部件進行維修或更換，此時的TR即為基于狀態維修的最佳維修時間。

5 實例分析

試驗采用3套Rexnord ZA-2155雙列滾子軸承，對2套軸承進行全壽命試驗，用于模型參數估計和模型驗證，對第3套軸承進行截尾壽命試驗，用于維修決策，分別記為軸承1、軸承2和軸承3。試驗得到2組全壽命試驗數據和1組截尾壽命數據。

全壽命試驗過程中，軸承的轉速和載荷保持在2 000 r/min和10 000 LB(約44 500 N)，并采用PCB加速度傳感器對軸承的振動信號進行采集，每1 h采集1次信號，采樣頻率為20 kHz，長度為1 s。

試驗過程中，軸承1和軸承2分別在966 h,982 h出現強烈的噪聲，此時認為軸承故障。用振動信號的能量均值作為間接狀態觀測數據y，磨損量為狀態數據x，2組全壽命數據的振動信號能量均值在整個壽命周期上的變化如圖4所示。

圖4 試驗軸承平均能量的變化過程

應用軸承1的間接狀態觀測數據對模型參數進行估計，并采用EM-PF算法對剩余壽命模型中的參數進行計算，結果見表1。

表1 模型參數計算結果

確定模型參數之后就能得到基于Gamma退化過程的狀態空間模型，通過計算求得軸承剩余壽命概率密度函數f(τi|y0:i)。對于軸承2的觀測數據，運用Matlab 7.0對剩余壽命預測方法進行實現，不同狀態檢測時刻(500～980 h)的剩余壽命概率密度函數如圖5所示，其中*表示剩余壽命的估計值，?表示剩余壽命真實值。剩余壽命預測值與實際值見表2。

圖5 不同檢測時刻剩余壽命概率密度

表2 剩余壽命預測值與實際值比較(軸承2) h

由上述結果可以看出，該模型對于此類試驗背景下的軸承剩余壽命預測具有一定的實用性，下面利用該模型對軸承3進行維修決策，確定使得單位時間費用最小的最佳維修時間。假設cp=600元，cf=1 200元，在得到軸承的剩余壽命概率密度函數之后，將f(τi|y0:i)代入到(2)式中，計算各維修時刻的單位時間費用，得到每次檢測時費用最低的更換時間TR，如圖6所示。圖中可以看出，在某一檢測時刻，隨著更換時間的增加，單位時間費用先減少后增加，在極值點處可以取得該次檢測的最佳維修時間。圖中?表示各檢測時刻使單位時間費用最低的維修時間。

圖6 不同檢測時刻所確定的維修更換時間

結合最佳維修時間確定的方法，通過Matlab計算，可以得出首次滿足TR-ti≤Δt條件時，TR=701 h。所以在701 h對軸承進行更換，可以使單位時間費用C(TR)達到最低,約為289.6元。

6 結束語

為解決維修決策的實際問題，重點研究了剩余壽命預測方法和維修決策優化模型。首先根據產品狀態的退化過程的特征建立了基于Gamma 退化過程的狀態空間模型，并對模型參數的估計方法進行了說明，進而通過該模型實現剩余壽命的預測，得到產品的剩余壽命概率密度函數。最后以滾子軸承壽命試驗得到的數據為例，驗證了剩余壽命預測模型，并以費用最小為目標建立了維修決策模型，對另一軸承進行了維修決策的優化，實現了剩余壽命預測與維修決策的有機結合，實例證明了所建立模型的有效性和可行性。