數(shù)據(jù)驅(qū)動的設(shè)備時變可靠度評價及故障預(yù)測

林子昕，王少琦，安維中，別海燕

（中國海洋大學(xué)化學(xué)化工學(xué)院，山東青島266100）

化工行業(yè)復(fù)雜的生產(chǎn)工藝條件對設(shè)備的要求很高，設(shè)備是工業(yè)生產(chǎn)運(yùn)行的基礎(chǔ)，設(shè)備運(yùn)行的可靠性直接影響裝置的連續(xù)穩(wěn)定運(yùn)行。由于設(shè)備的缺陷直接造成事故或者間接造成的事故會給企業(yè)造成巨大的損失，同時由于設(shè)備缺陷造成裝置停工，影響企業(yè)正常生產(chǎn)，帶來的經(jīng)濟(jì)損失也是巨大的。因此，保證化工生產(chǎn)裝置的長周期穩(wěn)定運(yùn)行越來越重要，設(shè)備的可靠性必須得到重視。

可靠性評估是指根據(jù)設(shè)備在整個壽命周期內(nèi)的可靠性結(jié)構(gòu)、分布模型以及有關(guān)的可靠性信息，利用數(shù)理統(tǒng)計(jì)的方法對所評定產(chǎn)品的可靠性指標(biāo)給出估計(jì)的過程[1]?？煽啃栽u估分為單元可靠性評估和系統(tǒng)可靠性評估，一般來說，“單元”是指作為單獨(dú)研究或單獨(dú)試驗(yàn)對象的任何元器件、零件，甚至一臺完整的設(shè)備；而“系統(tǒng)”則包括了組成系統(tǒng)的各單元以及各單元間的可靠性結(jié)構(gòu)[2]。近年來，可靠性評估技術(shù)迅速發(fā)展，現(xiàn)有的單元可靠性評價方法可大致分為基于物理失效模型、數(shù)據(jù)驅(qū)動的失效統(tǒng)計(jì)模型和基于經(jīng)驗(yàn)?zāi)Ｐ偷目煽啃栽u估方法。在實(shí)際應(yīng)用中，由于難以獲取復(fù)雜、高可靠性設(shè)備失效機(jī)理的物理模型，數(shù)據(jù)驅(qū)動的可靠性評估與故障預(yù)測方法成為近年來研究的主流，通過分析設(shè)備失效數(shù)據(jù)發(fā)現(xiàn)其失效分布類型，從而進(jìn)行統(tǒng)計(jì)推斷。

相比于電子和微電子產(chǎn)品的可靠性分布模型近似服從于指數(shù)分布[3]，一些學(xué)者[4-5]通過研究指出，機(jī)械設(shè)備產(chǎn)品的可靠性分布模型大多數(shù)服從威布爾分布模型。文獻(xiàn)[6]采用威布爾分布與極大似然估計(jì)相結(jié)合的方法建立數(shù)學(xué)模型來間接評估尾軸承的可靠壽命；文獻(xiàn)[7]通過改進(jìn)經(jīng)驗(yàn)分布函數(shù)的建立途徑，優(yōu)化威布爾分布參數(shù)的求解過程，進(jìn)而開發(fā)一種液壓支架可靠性評估的改進(jìn)威布爾分布法。威布爾分布能夠很好地擬合不同類型失效率的產(chǎn)品壽命數(shù)據(jù)，在設(shè)備可靠性建模中應(yīng)用最廣泛。

在工程領(lǐng)域中，各種主客觀因素的影響，結(jié)構(gòu)的抗力性能在服役期間逐漸減弱，常把這種考慮結(jié)構(gòu)抗力隨時間變化的可靠度稱為時變可靠度[8]。近年來，時變可靠度的計(jì)算日益受到重視。文獻(xiàn)[9]考慮了一般大氣環(huán)境下銹蝕因素對鋼筋混凝土梁可靠度的影響，根據(jù)極限狀態(tài)方程建立時變可靠度模型。文獻(xiàn)[10]研究了飛機(jī)鎖緊機(jī)構(gòu)轉(zhuǎn)動關(guān)節(jié)的磨損問題，結(jié)果表明關(guān)節(jié)磨損隨時間的變化關(guān)系對鎖緊機(jī)構(gòu)的可靠性有重要影響。

由于化工設(shè)備的運(yùn)行工況、載荷、應(yīng)力、強(qiáng)度等參數(shù)隨時間或載荷作用次數(shù)等壽命指標(biāo)在不斷變化，化工設(shè)備的可靠性是一個典型的時變過程，傳統(tǒng)的可靠性分析理論與方法不能較好地體現(xiàn)設(shè)備的動態(tài)特性，難以準(zhǔn)確地反映設(shè)備的可靠性隨壽命指標(biāo)的變化規(guī)律。本文在傳統(tǒng)威布爾（Weibull）分布可靠性評估模型基礎(chǔ)上引入隨時間變化的設(shè)備可靠性影響因素，建立多變量故障速率模型，對設(shè)備進(jìn)行時變可靠度評價以及故障預(yù)測。

1 基于威布爾分布可靠性的評估模型

用于評估機(jī)械設(shè)備或零件可靠性的數(shù)學(xué)模型包括指數(shù)分布、正態(tài)分布、對數(shù)正態(tài)分布和威布爾分布，對比其他數(shù)學(xué)模型，威布爾分布兼容性好，對各種類型的數(shù)據(jù)擬合能力強(qiáng)，可以較全面地描述設(shè)備不同失效期的失效過程與特征。

1.1 威布爾分布數(shù)學(xué)模型

二參數(shù)威布爾分布模型的失效率函數(shù)按式(1)計(jì)算。

式中，t為故障時間，d（或h、s 等）；β為威布爾分布形狀參數(shù)；η為威布爾分布尺度參數(shù)。

在威布爾分布中最重要的參數(shù)是形狀參數(shù)，它可以決定概率密度曲線的基本形狀，尺度參數(shù)起放大或縮小曲線的作用，但不影響分布曲線的形狀。通過改變形狀參數(shù)可以表示不同階段設(shè)備的失效情況，同時，改變形狀參數(shù)的值威布爾分布也可以作為其他分布的近似，如可將形狀參數(shù)設(shè)為合適的值近似正態(tài)、對數(shù)正態(tài)、指數(shù)等分布。

1.2 參數(shù)估計(jì)

計(jì)算威布爾分布的模型參數(shù)有多種方法，目前常用的方法有極大似然估計(jì)法、矩估計(jì)法和最小二乘法估計(jì)法等。文獻(xiàn)[11]比較了以上各種計(jì)算方法的特點(diǎn)，結(jié)果表明極大似然估計(jì)法計(jì)算得出的數(shù)據(jù)更能滿足實(shí)際要求，因此本文用極大似然估計(jì)法來求解基于威布爾分布可靠性評估模型的形狀參數(shù)β和位置參數(shù)η。

設(shè)總體X～Weibull(β,η),X1,X2,…,Xn為來自總體的n個獨(dú)立的樣本，則其似然函數(shù)見式(3)。

將式(3)取對數(shù)，可得對數(shù)似然函數(shù)，見式(4)。

使用Newton-Raphson 迭代法計(jì)算上述非線性方程組，即可求出極大似然法估計(jì)的模型形狀參數(shù)β和尺度參數(shù)η，最終得到基于威布爾分布的可靠性評價模型。

2 時變可靠度評價模型

在化工行業(yè)中，即使一批結(jié)構(gòu)形式、尺寸、材料和加工方法上都相同的設(shè)備，它們的使用壽命卻各有長短，不盡相同，其原因是多方面的，工作環(huán)境如溫度和壓力的變化會對設(shè)備的壽命產(chǎn)生很大的影響。因此，在傳統(tǒng)威布爾分布的可靠度評價模型的基礎(chǔ)上進(jìn)行改進(jìn)，引入設(shè)備可靠性的影響因素，采用數(shù)據(jù)驅(qū)動的方法，建立一種基于威布爾分布的設(shè)備時變可靠度評價模型。

2.1 灰色關(guān)聯(lián)分析

灰色關(guān)聯(lián)分析是對一個系統(tǒng)發(fā)展變化態(tài)勢的定量描述和比較的方法，研究系統(tǒng)中特征因素與影響因素之間的關(guān)系，計(jì)算各影響因素序列與特征序列的相關(guān)程度，依照關(guān)聯(lián)度大小得出結(jié)論?；ぴO(shè)備生產(chǎn)運(yùn)行中影響因素較多，本文通過灰色關(guān)聯(lián)分析計(jì)算關(guān)聯(lián)度來確定影響設(shè)備可靠性的主要因素。關(guān)聯(lián)系數(shù)計(jì)算公式見式(6)。

式(7)為數(shù)列xi對特征數(shù)據(jù)序列x0的關(guān)聯(lián)度。即把各個關(guān)聯(lián)系數(shù)集中為一個平均值，以此來反映整體關(guān)聯(lián)程度[12]。

2.2 響應(yīng)面法

近年來，元模型被廣泛用于通過使用代理模型預(yù)測系統(tǒng)響應(yīng)來降低可靠性分析的計(jì)算成本[13]，響應(yīng)面法等基于回歸的元模型代理模型得到了廣泛應(yīng)用。響應(yīng)面法（Response Surface Method）是利用數(shù)學(xué)和統(tǒng)計(jì)學(xué)技術(shù)解決復(fù)雜系統(tǒng)的輸入（隨機(jī)變量）與輸出（系統(tǒng)響應(yīng)）關(guān)系的方法。

響應(yīng)面法分為多項(xiàng)式響應(yīng)面法和神經(jīng)網(wǎng)絡(luò)響應(yīng)面法，神經(jīng)網(wǎng)絡(luò)響應(yīng)面法計(jì)算過程十分復(fù)雜，工程上的計(jì)算多采用多項(xiàng)式響應(yīng)面法。多項(xiàng)式響應(yīng)模型見式(8)。

式中，ai,k，bi,j,k,n(i=1,2,...,n)為表達(dá)式的待定系數(shù)。

從響應(yīng)面函數(shù)表達(dá)式可看出，若隨機(jī)變量個數(shù)取N，共有(N+1)(N+2)/2 個待定系數(shù)。綜合響應(yīng)式函數(shù)的計(jì)算量、計(jì)算精度和計(jì)算穩(wěn)定性等方面因素，響應(yīng)面函數(shù)最常用的是二次多項(xiàng)式形式，式(8)中K= 2即為二次多項(xiàng)式響應(yīng)模型。

通過可決系數(shù)R2對響應(yīng)關(guān)系式進(jìn)行檢驗(yàn)，見式(9)。

一些時變可靠性分析方法采用替代模型來逼近復(fù)雜系統(tǒng)的隱式狀態(tài)函數(shù)[14]，本文采用響應(yīng)面法建立失效率模型中形狀參數(shù)β和尺度參數(shù)η與設(shè)備運(yùn)行影響因素的關(guān)系，將影響因素引入設(shè)備失效率計(jì)算模型，建立多變量失效速率模型。

2.3 時變可靠度模型

在傳統(tǒng)威布爾分布的可靠性評價模型中，用極大似然法可計(jì)算得到設(shè)備的形狀參數(shù)β和尺度參數(shù)η，根據(jù)響應(yīng)面法可計(jì)算得到模型參數(shù)與影響因素的響應(yīng)關(guān)系式，分別見式(10)和式(11)。

將式(10)和式(11)代入威布爾分布的失效率模型可得出故障速率模型，見式(12)。

計(jì)算時變可靠度需考慮設(shè)備在不同時間段內(nèi)的運(yùn)行工況，運(yùn)行條件不同，故障速率模型也不同，設(shè)備時變可靠度計(jì)算模型如下。

若t∈[0,t1],可直接計(jì)算可靠度，見式(13)。

3 應(yīng)用實(shí)例

在化工生產(chǎn)過程中，機(jī)泵設(shè)備始終都發(fā)揮著至關(guān)重要的作用，以某公司催化裂化裝置中的12 類機(jī)泵運(yùn)行數(shù)據(jù)為例，對設(shè)備時變可靠度進(jìn)行評價。根據(jù)機(jī)泵的歷史故障時間數(shù)據(jù)，采用極大似然法計(jì)算失效率函數(shù)中的形狀參數(shù)β和尺度參數(shù)η以及平均故障間隔時間，如表1所示。

機(jī)泵運(yùn)行過程中影響因素較多，如泵的揚(yáng)程，泵中流體的流量，流體的溫度、密度，泵的進(jìn)壓、出壓以及泵的功率等，機(jī)泵運(yùn)行狀況見表2。采用灰色關(guān)聯(lián)分析法比較影響因素與模型參數(shù)β和η之間關(guān)聯(lián)度的大小，選擇關(guān)聯(lián)度較大的因素為模型的主要影響因素。

以機(jī)泵的形狀參數(shù)β和尺度參數(shù)η為特征參考數(shù)列，分別計(jì)算其與泵的進(jìn)壓、泵的出壓、泵的功率、泵的揚(yáng)程、進(jìn)出泵的流體壓差、流體溫度和流體流量等影響因素數(shù)列的關(guān)聯(lián)度。

圖1 為形狀參數(shù)β與各影響因素的關(guān)聯(lián)系數(shù)圖。圖2 為尺度參數(shù)η與各影響因素的關(guān)聯(lián)系數(shù)圖。

從圖中可以看出，除了泵的進(jìn)壓與模型參數(shù)的關(guān)聯(lián)度較低，其他影響因素與模型參數(shù)的關(guān)聯(lián)度均較高，通過計(jì)算平均關(guān)聯(lián)度，得到影響因素與模型參數(shù)的關(guān)聯(lián)度。模型參數(shù)與各影響因素關(guān)聯(lián)度大小順序?yàn)槊芏圈眩境鰤簆2＞揚(yáng)程H＞功率W＞流量V＞溫度T＞進(jìn)壓p1。

考慮到12 類機(jī)泵數(shù)據(jù)有限，選擇泵中流體的密度、泵的出壓p2和揚(yáng)程H，作為機(jī)泵運(yùn)行中的主要影響因素，采用響應(yīng)面法建立選擇的影響因素數(shù)據(jù)與設(shè)備形狀參數(shù)和尺度參數(shù)之間的響應(yīng)關(guān)系式。響應(yīng)關(guān)系式中變量個數(shù)為3，需要求解的二次響應(yīng)式的響應(yīng)系數(shù)個數(shù)是10，在響應(yīng)關(guān)系式中，aij為變量一次項(xiàng)和平方項(xiàng)的系數(shù)；bij為變量交叉項(xiàng)的系數(shù)，c為常數(shù)項(xiàng)。在Matlab中，計(jì)算得到響應(yīng)系數(shù)，通過Matlab擬合的響應(yīng)關(guān)系式計(jì)算的數(shù)據(jù)與設(shè)備形狀參數(shù)和尺度參數(shù)數(shù)據(jù)的對比如圖3和圖4所示。

圖3 響應(yīng)關(guān)系式擬合的形狀參數(shù)與原始數(shù)據(jù)對比圖

圖4 響應(yīng)關(guān)系式擬合的尺度參數(shù)與原始數(shù)據(jù)對比圖

由圖可知，擬合數(shù)據(jù)曲線與原始數(shù)據(jù)曲線基本一致，分別計(jì)算曲線的可決系數(shù)R2，分別為R2（β響應(yīng)式）=0.9862、R2（η響應(yīng)式）=0.9037?？蓻Q系數(shù)均大于0.90，響應(yīng)關(guān)系式擬合效果較好。

根據(jù)響應(yīng)系數(shù)，得出威布爾分布形狀參數(shù)β和尺度參數(shù)η與影響因素（出壓、揚(yáng)程、壓差）的關(guān)系式，見式(19)和式(20)。

在失效率函數(shù)中，引入設(shè)備運(yùn)行的影響因素，將所求得的關(guān)系式代入設(shè)備失效率模型中，可得出故障速率模型和時變可靠度模型，分別見式(21)和式(22)。

若機(jī)泵的運(yùn)行條件隨時間發(fā)生變化，可根據(jù)上述兩個模型計(jì)算設(shè)備時變可靠度。表3為機(jī)泵1的運(yùn)行條件，根據(jù)所給條件計(jì)算機(jī)泵1的時變可靠度和MTBF。假設(shè)機(jī)泵在0～100d運(yùn)行參數(shù)不變，其中密度為1000g/cm3，出壓為6.2MPa，揚(yáng)程為600m，在100～200d 和200～300d 時，由于物性參數(shù)和操作參數(shù)發(fā)生變化，3 種影響因素均有一定程度的增大，計(jì)算其模型參數(shù)，時變可靠度和平均剩余壽命如表4所示。

表3 機(jī)泵1運(yùn)行條件

表4 機(jī)泵的可靠性計(jì)算

從表4 中可以看出，機(jī)泵運(yùn)行200～300d 后，時變可靠度為0.3373，此時機(jī)泵的可靠度較低，在200～300d 時剩余壽命為負(fù)值，說明此時機(jī)泵可能已達(dá)到故障狀態(tài)。

若不考慮機(jī)泵運(yùn)行參數(shù)的變化，在泵中流體密度1000g/cm3、泵的出壓6.2MPa、揚(yáng)程600m操作條件下，運(yùn)行300d后，計(jì)算機(jī)泵的可靠度為0.4757，在運(yùn)行工況近似的情況下，時變可靠度計(jì)算結(jié)果為0.3373，證明了時變可靠度計(jì)算方法的可行性。

4 結(jié)論

本文在傳統(tǒng)威布爾分布的基礎(chǔ)上，改進(jìn)模型參數(shù)的計(jì)算方式，建立了數(shù)據(jù)驅(qū)動的設(shè)備時變可靠度評價模型，并對設(shè)備的故障進(jìn)行預(yù)測?？捎糜谟?jì)算新設(shè)備可靠度和MTBF的計(jì)算以及運(yùn)行工況隨時間發(fā)生變化的設(shè)備時變可靠度的計(jì)算和故障預(yù)測。在運(yùn)行工況發(fā)生改變的情況下，設(shè)備的時變可靠度計(jì)算模型可對設(shè)備故障時間進(jìn)行預(yù)測。

符號說明

H—— 揚(yáng)程，m

p1—— 泵的進(jìn)壓，MPa

p2—— 泵的出壓，MPa

Δp—— 泵的壓差，MPa

R—— 可靠度

T—— 流體溫度，°C

t—— 設(shè)備故障時間，d

V—— 體積流量，L/h

W—— 泵的功率，W

λ—— 失效率，%

β—— 威布爾分布形狀參數(shù)

η—— 威布爾分布尺度參數(shù)，d