基于比例故障率模型的柴油發動機視情維修決策

王海龍， 張 營， 左付山

(南京林業大學汽車與交通工程學院， 南京 210037)

柴油發動機作為汽車的動力部件，其結構復雜，運動部件多的特點，使得柴油機維修診斷較為困難。目前，對于車輛柴油機的維修主要采用定時維修和事后維修的策略，定時維修具有容易執行等優點，然而并沒有考慮發動機當前運行狀態；事后維修則是在故障發生之后再進行相應的維修，維修成本高，經濟損失大。為了提高經濟效益和社會效益，避免維修不足或者維修過度，相關學者提出了視情維修。視情維修是通過監測設備的狀態參數來反映設備的性能狀態，然后根據相應的狀態在故障發生之前采取一定的維修措施[1]，其最主要的就是建立設備運行狀態參數與設備完好程度之間的對應關系[2]，建立模型的方法包括基于設備可靠性、概率性以及卡爾曼濾波等建模方法。

董思辰等[3]選擇振動信號中的冗余參數作為協變量，利用威布爾比例故障率模型(Weibull proportional hazards model，WPHM)進行齒輪箱的運行可靠性評估，得到可靠度衰減時間與振動信號故障時間相近。劉璐[4]選擇航空機載設備電機軸承的振動信號，通過特征提取后，基于WPHM進行了其可靠性評估和壽命預測。蔣文博等[5]提出了一種基于比例風險模型與機器學習的混合方法，有效的預測了電梯的剩余壽命。通過上述的研究發現，對于機械設備來說，WPHM能綜合考慮設備運行狀態及工作時間，有效反映設備的狀態好壞，從而實現壽命預測與視情維護，因此柴油機作為機械設備，狀態參數復雜并且維修決策較為困難，通過建立WPHM進行視情維修將減少大量維修成本并且提高工作效率。

綜上所述，構建一種基于威布爾比例故障率模型的發動機視情維修策略，首先對監測到的柴油發動機在各個階段的潤滑油金屬含量數據進行了主成分分析(principal component analysis，PCA)數據降維并通過極大似然估計求出了所建立的模型參數；為綜合考慮成本和可用度，采用最小維修成本和最大可用度的綜合維修策略，求得相應的最優預防時間間隔以及失效率閾值，從而根據當前監測到柴油機金屬含量數據判斷是否需要對其進行預防性維修。

1 威布爾比例故障率模型

1.1 比例故障率模型

比例故障率模型前期主要應用于醫學領域[6]，其綜合考慮了設備運行時間和當前狀態信息之間的關系[7]，從而實現了對于設備健康狀況的評估。該模型的性質是不同個體的故障率函數成比例[8]，并且故障率函數不隨時間t的變化而變化。其基本形式為

h(t;Z)=h0(t)exp(rZ)

(1)

式(1)中：h(t;Z)為故障率函數；h0(t)為基本故障率函數，只與時間有關[9]；rZ=r1Z1+r2Z2+…+rpZp，Z為協變量，反映裝備運行狀態，按作用可分為外部協變量與內部協變量[10]，如振動幅度、金屬含量等，r是協變量系數，反映了與之對應的協變量對裝備故障率影響的嚴重程度[11]，ri的值越大，表明對應的協變量對故障率函數影響越大。

1.2 威布爾比例故障率模型的建立

在汽車可靠性研究中，指數分布、威布爾分布、正態分布、對數正態分布等都是較為常用的理論分布類型。由于柴油發動機為機械部件，一般機械部件的故障數據分布規律均服從威布爾分布，同時其對各種類型的數據具有很強的擬合能力，所以應用廣泛，因此上述基本故障率函數選擇較為經典和簡單的兩參數威布爾分布故障率函數作為比例故障率模型的基礎函數[12]，兩參數威布爾故障率函數具體形式為

(2)

式(2)中：β為形狀參數，可改變曲線的陡峭程度；η為尺度參數[13]，可改變曲線中心線位置。

由此，得到了威布爾比例故障率函數[14]為

(3)

故障概率密度函數f(t)、可靠度函數R(t)以及故障率函數h(t)定義如下：f(t)是指所有部件在時間間隔(t，t+Δt)內發生故障的個數占總數的比例；R(t)是指規定時間內完成規定功能的概率，即部件壽命T大于某一時刻t的概率；h(t)是指定義為時間(0，t)內不發生故障的條件下，下一個單位時間(t+Δt)內發生故障的概率。由此，可以得到[15]

(4)

(5)

(6)

式中：N為部件總個數；n(t)為到t時刻發生故障的個數。

對式(6)左右兩邊同時在0～t上積分，得可靠度函數R(t)，即

(7)

(8)

同理可得，相應的威布爾比例故障率可靠度函數R(t;Z)為

(9)

1.3 模型參數求解

為了通過威布爾比例故障率模型進行柴油發動機的可靠性評估，最關鍵的是對于模型中的三個參數進行求解β、η、r。對于三個參數估計的方法主要有最大似然估計法、最小二乘法、矩估計法、區間估計法等。由于極大似然估計法對于設備的不完全壽命數據是處理具有獨特的優勢，使得模型更加精確，因此本文中選擇極大似然估計法對于模型的三個參數進行估計。

極大似然估計的基本原理是通過估計的參數使得樣本出現的概率更大，似然函數的一般形式為

(10)

針對威布爾比例故障率模型，需要通過極大似然法求解，即求參數θ=(β，η，r1，r2，…，rp)何值時，使得似然函數的值最大。根據柴油發動機的失效密度函數f(t;Z)=h(t;Z)R(t;Z)，威布爾比例故障率模型似然函數為

exp(rZ)dt]}

(11)

式(11)中：F為所測數據中故障集；q為故障數據個數；N為數據總集。

為了求得似然函數的最大值，需要先對其進行取對數，得到對數似然函數，然后分別對于β、η、r進行一階求偏導并令其等于0，可以得到對數似然函數方程組。似然函數及其一階方程組[16]為

(12)

(13)

(14)

對于上述對數似然函數方程組求解方法主要有牛頓迭代法、遺傳算法[17]、粒子群法。

牛頓迭代法主要是通過泰勒一階展開，從初始值開始，分別對下一個值進行泰勒展開，不斷循環，直到達到理想精度或最大迭代次數停止計算[18]，從而求得模型參數β、η、r。

原理為

f(x)=f(x0)+(x-x0)f′(x0)

(15)

步驟為

(16)

具體求解：由于上述得到了似然函數一階方程組，所以要求得似然函數方程組的解，還需要對于其進行二次求導[19]，即

(17)

根據牛頓迭代法可得

(18)

式(18)中：β0、η0、r0為初始值。

通過牛頓迭代法求解得到3個待估參數之后，便得到了威布爾比例故障率模型以及相應的可靠度模型，繼而就可以選擇最小成本維修策略或者最大可用度維修策略建立決策模型，求取最優預防時間間隔。

2 綜合維修優化模型

2.1 最小維修成本

針對汽車的柴油發動機，由于經濟性在維修中的是首先被考慮的，其故障后維修費用遠大于預防性維修費用，所以選擇最小維修成本法作為決策目標，其具體模型求解如下所示。

當以單位時間內柴油發動機維修成本最小為目標，制定相應地最優維修策略。單位時間維修成本為E，其計算公式[20-22]為

(19)

式(19)中：CPM為部件的預防性維修成本；CCM為部件故障的時候維修成本；N(t)=(t/η)β為0～t內部件的平均失效數。

為了求得E的最小值，對t進行一階求導，并等于0，得到最優預防時間間隔T為

(20)

2.2 最大可用度法

可用度A(t)表示t時刻系統處于正常狀態的概率，其計算公式[23]為

(21)

(22)

求A(t)的最大值，即轉化為求t取何值時，α可以取到最小值。

2.3 綜合維修優化模型

由于單方面只考慮減少柴油發動機維修成本，有可能會造成可用度無法保障，為了綜合考慮柴油發動機的維修成本最小以及可用度最大，為此提出了以可用度為約束，以單位維修最少為目標，建立柴油發動機綜合維修優化模型[24]，即

(23)

2.4 維修決策模型

根據綜合維修優化模型得到最優預防維修時間間隔T，在故障數據中找出一組與之最接近的失效時間與伴隨時間變量數據，帶入威布爾比例故障率模型，得到失效率閾值下限h。在部件工作正常時，故障率低于h，即

(24)

然后根據監測的周期Δt,繼續代入模型計算，得到失效率模型閾值上限h*，當設備當前狀態及壽命數據代入模型后大于h*,表明設備已經處于維修區，即

(25)

將式(25)經過取對數等變形后，得到了對數風險度rZ，并將維修區域分為正常工作區、加強監測區、維修區[25-26]，即

(26)

基于綜合優化維修決策流程如圖1所示。

3 實例分析

3.1 模型選擇

對于柴油發動機的壽命數據，它有可能服從幾種分布類型，如指數分布、對數正態分布、威布爾分布等[27]。為了判斷壽命數據最符合哪種分布類型，首先需要通過最大似然法估計出各種分布的參數值，然后通過擬合優度檢驗方法進行擬合優度檢驗[28]。然而，有時經過分布的擬合優度檢驗判斷出符合條件的分布類型不止一種，所以引入了赤池信息準則(akaike information criterion，AIC)進行分布模型的最終選擇[29]。AIC信息準則是基于最大熵原理，能夠檢驗出不同模型間差異的顯著性，可通過計算錯判概率進行模型選擇分析，其計算公式為

(27)

為了驗證上文中選擇威布爾模型的有效性，基于柴油發動機壽命數據對于對數分布、對數正態分布以及威布爾分布的AIC值進行了計算與判定。柴油發動機的故障壽命數據如表1所示。

表1 故障壽命數據

首先采用極大似然估計的方法得到了3種分布的參數值，繼而選擇擬合優度檢驗方法中的K-S檢驗[30]對于3種分布的p值和H值進行了求解，求解結果如表2所示。

從表2可以看出，3種分布的K-S檢驗的H值均為0，表示接受原假設，檢驗結果均服從原假設分布。為了進一步確定最佳分布，下面根據AIC公式計算得到了3種分布的AIC值如意表3所示。

表2 參數估計值和K-S檢驗結果

表3 3種分布類型AIC值

根據AIC準則，AIC值最小的分布類型作為最佳的分布類型[31],所以推斷出威布爾分布為故障數據的最佳分布類型。

3.2 模型的建立與求解

在不同階段，柴油發動機的運行狀態與柴油發動機的潤滑油中的各種金屬含量有著直接關系，因此本文選擇其中鐵、鈣、鎂、鉛在潤滑油中的含量作為模型的協變量。在監測的數據中，包括柴油發動機工作時間、換油時間、4種金屬含量等數據。其部分數據如表4所示。

通過表4可以發現，在同一時刻，潤滑油中4種金屬含量數據量級差距較大，為了避免較大數據淹沒其他數，提高模型精度，需要先將原始數據歸一化。同時，為了減少協變量的個數，并得到一組更加全面的有效的特征協變量，可以對于數據進一步的處理，即通過主成分分析得到新的協變量[32-33]，在MATLAB種將主成分貢獻度設置為90%，便得到了三個相互獨立的新的特征協變量Z1、Z2、Z3，回歸變量的相關系數矩陣和相關系數矩陣的特征值與貢獻率分別如表5、表6所示。

表5 回歸變量的相關系數矩陣

表6 相關系數矩陣的特征值與貢獻率

將處理后的監測數據帶入對數似然函數中，并使用牛頓迭代法借助MATLAB對于模型的3個參數進行求解，得到參數β=2.25，η=14 846，r1=-2.883，r2=0.741，r3=-1.633 ,因此威布爾比例故障率模型為

(28)

3.3 維修閾值確定

以綜合維修策略作為決策方法，通過調查及歷史經驗，柴油發動機的預防性維修費用CPM平均為100元，故障后維修費用CCM平均為2 000元，因此CPM/CCM=0.05,將上述結果代入式(19)中，畫出時間t與單位時間維修成本E的曲線如圖2所示，并求出相應的最優預防時間間隔。

圖2 平均時間維修成本曲線

通過圖2以及相應的計算可知，當T=3 486時，單位時間維修成本可以達到最小值，并滿足可用度要求，因此T=3 486為最優預防時間間隔。然后，在歷史監測的壽命與油液金屬含量數據中，找到一組與壽命時間最近的一組數據及其對應的3個新的協變量，將其帶入WPHM模型后，可以得到相應的故障率閾值h=3.32×10-5，再將h代入式(26)中，即可得到決策閾值下限曲線為

0.741×(-0.075 2)-1.636×

(-0.083 3)]=3.32×10-5

(29)

為了進一步確定維修決策閾值上限，設置監測時間間隔時間為1 260 h，即得到時間T=4 746，同樣找到一組與壽命時間最近的一組數據及其對應的3個新的協變量，代入模型后便可得到對應的故障率閾值h*=6.75×10-5，以及決策閾值上限曲線。決策閾值上下限曲線如圖3所示。

圖3 決策閾值上下限曲線

在后續的維修決策中，將監測到的柴油發動機壽命時間和潤滑油金屬含量數據經過數據處理后，帶入決策模型中，如果得到的數值處于決策曲線下方，表明工作正常，繼續進行監測；如果得到的數值處于上限與下限之間，則表明需要加強監測，縮短監測的時間間隔；如果求得的數值處于決策曲線的上方，則需要對于柴油發動機采取相應的維修措施。

3.4 模型驗證

為了驗證上述建立的視情決策維修模型的有效性與可用性，選擇監測到的一組同類型的柴油發動機的潤滑油中金屬含量數據以及對應的運行狀況進行了分析驗證。該發動機從1994年1月開始運行，在運行到1994年12月出現了故障，總共運行了7 523 h。將柴油發動機的金屬含量數據經過PCA后，帶入維修決策模型，得到相應的對數風險值并在圖中畫出，如圖4所示。

圖4 柴油發動機狀態維修曲線圖

從圖4可以看出，該柴油發動機在運行到2 516 h的時候，其最接近決策曲線下限，這時需要加強監測；在運行6 196 h后，其狀態接近決策曲線上限，此時應當采取一定的預防性維修措施。根據該發動機的實際狀態，由于沒有采取相應的維修措施，其在運行到7 523 h發生了故障，該實驗結果表明本文提出的決策方法有效合理，可用于實際的該型發動機視情維修決策。

4 結論

基于收集的柴油發動機壽命數據及其對應的潤滑油金屬含量數據，在對其進行PCA后，選取了累計貢獻率達到90%的3個特征協變量，使用最大似然估計的方法估計得到了威布爾比例危險率模型參數，建立了威布爾比例故障率模型；進而以最小維修成本為目標，以最大可用度作為約束，計算求得了最優預防時間間隔，建立了同類型柴油發動機維修決策模型，為后續的狀態維修提供了有效的決策方法，不僅最大可能地降低了維修成本，又保證了柴油發動機的安全運行和最大可用時間。最后以同類型發動機的狀態及壽命數據為例，驗證了該模型可以有效地通過壽命和金屬含量數據判斷柴油發動機是否需要進行維修，具有較大的應用于視情維修決策的實用價值。