秦犖晟， 陳曉陽， 沈雪瑾

(上海大學 機械自動化工程系,上海 200072)

小樣本下基于競爭失效的軸承可靠性評估

秦犖晟， 陳曉陽， 沈雪瑾

(上海大學 機械自動化工程系,上海 200072)

軸承是機械產品中重要的組成部分，其性能和壽命與機械設備的運行壽命密切相關。為獲得準確的軸承可靠性評估結果，應當綜合考慮不同失效模式對軸承可靠性的影響。針對軸承試驗中完全失效的數據，利用Bootstrap法構造軸承壽命分布參數的先驗分布，依據Bayes法估計出相應的后驗分布，并對后驗期望修偏即可獲得軸承壽命分布參數。通過進一步分析軸承的振動性能退化數據，獲得軸承局部失效的壽命分布。利用Copula函數對軸承的完全失效的壽命分布模型和局部失效壽命模型進行綜合分析，借助試驗數據點的經驗Kendall相關秩估計Copula函數的相關參數，求得軸承競爭失效下的可靠性分析結果，獲得的分析結果有助于尋找軸承設計過程中存在的缺陷，提高軸承的可靠度。

軸承；Bootstrap方法；Bayes法；性能退化；Copula函數；可靠性

近年來，軸承產品的使用性能以及軸承的可靠性分析得到了越來越多的關注，如何提高滾動軸承的使用性能、使用壽命，以及如何提高軸承可靠性評估精度成為熱點。Dong等[1]基于支持向量機和馬爾科夫模型對軸承退化過程做了預測，并通過實例證明了所提出方法的有效性；Edwin等[2]介紹了三種可以通過軸承振動退化數據估計球軸承剩余壽命的方法；Zhang等[3]介紹了如何選取合適的軸承運行參數確定軸承的剩余使用壽命。

基于多種失效模式競爭失效的可靠性分析可以提高軸承可靠性評估的準確性。王華偉等[4]通過構造混合Weibull分布模型對航空發(fā)動機的競爭失效進行了分析；唐家銀等[5]通過使用Copula函數對多故障模式相關情況下的長壽命產品進行了可靠性分析；常春波等[6]基于系統自然退化過程和沖擊過程的分布函數推導出基于競爭失效的系統可靠度函數。Bocchetti等[7]通過假設退化型失效和突發(fā)型失效相互獨立的情況下，給出了產品在競爭失效下的可靠性。

產品的可靠性試驗常常存在樣本量少的問題，這會造成對產品可靠性模型參數的估計產生較大誤差?？煽啃灶I域中的小樣本問題得到了越來越多的重視，蔣喜等[8]應用虛擬樣本增廣法解決了極小子樣下電主軸性能退化的可靠性分析；Xiong 等[9]針對小樣本的測試數據提出了一種計算方法，求解產品結構安全的可靠性；Picheny 等[10]通過采用Bootstrap法對有限樣本下的試驗數據進行分析求得了產品的可靠度。近年來，Bayes方法也被廣泛應用于對Weibull分布參數的估計[11-12]。金光[13]提出了一種綜合性能與壽命數據的Bayes-Bootstrap法，通過將Bayes方法與Bootstrap方法相結合，計算小樣本問題的可靠度。

本文將結合Bootstrap方法以及Bayes分析方法，分析處理軸承試驗中的小樣本數據，實現小樣本試驗數據下壽命分布參數的估計；運用Copula函數建立滾動軸承在不同失效模式下競爭失效的可靠性評估模型，并結合Kendall相關秩給出一種Copula函數參數估計的簡便方法，實現對競爭失效下滾動軸承的可靠性評估，提高滾動軸承可靠性評估的準確度。

1 機械產品失效模式的分類

Blanche等[14]提出的失效模式分類，指出產品喪失既定功能的失效被稱為擴展失效，擴展失效可分為完全失效(Complete Failures)和局部失效(Partial Failures)。完全失效即產品完全喪失既定功能的失效，局部失效即產品喪失部分既定功能的失效。

在機械產品實際使用過程中，機械產品發(fā)生的失效可以認為是完全失效和局部失效的競爭失效。

2 多失效模式的軸承可靠性評估模型

2.1 軸承壽命分布模型

由于Weibull分布可以很好地模擬多種分布，因此假設滾動軸承的試驗壽命t服從參數為(β,η)的二參數Weibull分布，其中β為Weibull分布的形狀參數，η為Weibull分布的尺度參數，Weibull分布的概率密度函數如式(1)所示。

(1)

則產品試驗壽命t的分布函數如式(2)所示。

(2)

2.2 完全失效的滾動軸承可靠性評估模型

在小樣本情況下，應用Bayes-Bootstrap法來估計Weibull分布參數，求出滾動軸承完全失效下的壽命分布。

2.2.1 Bayes-Bootstrap方法的實現

Bootstrap方法是一種再抽樣理論，肖支才等提出了一種非參數的Bootstrap抽樣方法，這種抽樣方法獲得的再生樣本點不僅包含原始樣本點，還適當地在原始樣本點的基礎上進行拓展。抽樣方法的具體步驟如下：

(1)利用計算機產生[0,1]區(qū)間均勻分布的隨機數δ;

(2)令θ=(n-1)δ，i=[θ]+1，[·]為向下取正；

(3)令t*=ti+(θ-i+1)(ti+1-ti)；

其中ti是原始樣本t1,t2,…,tn按從小到大排序后得到的第i個樣本數據，則t*即為所產生的再生樣本數據。

重復上述過程B次，就可以獲得B組Weibull分布參數β和η的估計值，從而構造β和η的先驗概率密度π1(β)和π2(η)。本文假設參數(β,η)是相互獨立的，且β和η的先驗概率分布均為正態(tài)分布，則β和η的聯合先驗概率密度為

π(β,η)=π1(β)π2(η)

(3)

假設某一產品試驗壽命值為t1,t2,…,tn，該產品的試驗壽命值服從Weibull分布，則該樣品的壽命t1,t2,…,tn在給定參數β和η條件下的極大似然估計可以表示為

(4)

(5)

在β和η的聯合先驗概率密度為π(β,η)=π1(β)π2(η)的情況下，根據Bayes公式，有

f(β,η|t1,t2,…tn)=

(6)

式(6)中，Θ表示(β,η)的取值空間。由式(3)、式(6)可知，Weibull分布參數β和η的后驗概率密度分別為

(7)

(8)

式(7)、(8)中，Θ1表示β的取值空間，Θ2表示η的取值空間。由于式(7)、(8)中的積分不容易求解，因此參數β和η的后驗分布很難得到具體的形式，可以用β和η的后驗期望作為β和η的Bayes點估計值。由(7)、(8)可以得到β和η的Bayes點估計值:

(9)

(10)

在可靠性試驗的評估結果中，β和η的值分別存在一定有界的范圍[β1,β2]和[η1,η2]，使得式(11)和式(12)成立。

(11)

(12)

這樣式(9)、(10)可以被進一步改寫成式(13)、(14)的形式。

(13)

(14)

對式(13)、(14)作有限范圍內的數值積分，求出β和η后驗分布的期望值作為Weibull分布參數的估計值。

(15)

(16)

2.2.2 模型的驗證

利用計算機隨機生成5個服從分布參數已知的Weibull分布隨機數，運用2.2.1提出的方法分別估計5個隨機數的Weibull分布參數。針對不同的分布參數分別進行試驗，獲得的結果如表1所示。

表1Weibull分布參數估計結果對比

Tab.1ThecomparisonoftheWeibulldistributionparametersestimatedbydifferentmethods

真實值(β,η)極大似然估計值(β,η)2.2.1方法計算結果(β,η)(1.5,1500)(1.08,1196.08)(1.25,1343.6)(1.8,1800)(1.49,1248.83)(1.75,1467.73)(2,1600)(1.42,1083.27)(1.6,1282.31)(2.5,2000)(1.79,1745.78)(1.94,1923.36)(2.5,2200)(2.07,1690.88)(2.31,1937.76)(3,2000)(2.49,1606.1)(3.00,1939.39)

由表1的計算結果可知在樣本量較小的情況下，2.2.1的方法求得的結果更準確。

2.3 局部失效的滾動軸承可靠性評估模型

軸承的局部失效壽命數據可以通過追蹤其性能參數獲得。基于性能退化的產品可靠性評估目前已有多種較為成熟的計算模型[17-18]，對于滾動軸承產品，也可以選取其性能指標，考慮采用擬合退化軌跡的方法估計出滾動軸承達到失效閾值所需要的時間，如圖1所示。

圖1 性能退化軌跡圖Fig.1 The figure of the performance degradation trace

通過退化軌跡所求得的軸承達到失效閾值的時間即軸承的失效時間，由此運用2.2提出的計算方法即可計算軸承的性能退化壽命分布。

2.4 競爭失效的滾動軸承可靠性評估模型

針對多失效模式競爭可靠性分析，運用Copula函數分析求解基于多失效模式滾動軸承可靠度評估模型。

由Sklar定理[19]可知，聯合概率密度f(x,y)可以表示為式(17)的形式：

f(x,y)=c(F(x),G(y))f(x)g(y)

(17)

式中，c(F(x)和G(y))表示相應的Copula密度函數，f(x)和g(y)分別為F(x)、G(y)的概率密度函數。

在常見的Copula函數形式中，Gumbel Copula函數結構與Lee[20]提出的一類二元Weibull分布的生存函數結構相似，因此選用Gumbel Copula對軸承進行競爭失效分析。令F(x)=u，G(y)=v，則Gumbel Copula的函數結構如式(18)所示。

θ∈(0,1]

(18)

對Copula函數的參數估計有多種估計方法，但計算過程復雜，不便于工程應用。然而利用試驗數據點的經驗Kendall相關秩估計Copula函數的相關參數，可以簡化計算，便于工程應用。

Kendall相關秩被描述為任意兩對數據和諧的概率與不和諧的概率之差[19]。如果用τ表示Kendall相關秩，假設數據組{(x1,y1),(x2,y2),…(xn,yn)}是連續(xù)隨機變量向量(X,Y)的n個隨機觀測值，如果(xi-xj)(yi-yj)>0，則稱數據對(xi,yi)與數據對(xj,yj)是和諧的；如果(xi-xj)(yi-yj)<0,則稱(xi,yi)與(xj,yj)是不和諧的，用P[(xi-xj)(yi-yj)>0]表示兩個數據對和諧的概率，P[(xi-xj)(yi-yj)<0]表示兩個數據對不和諧的概率，則Kendall相關秩τ可以表示為

(19)

由文獻[21]可知，Gumbel Copula函數的參數θ與其相應的Kendall相關秩τ存在如下關系：

τ=1-θ

(20)

因此，可以將試驗數據點的經驗Kendall相關秩看作Gumbel Copula函數的Kendall相關秩τ，這樣就可以快速而又簡便的估計出Gumbel Copula函數參數θ的值。

3 實例分析

針對某型號滾動軸承進行4 000 h的定時截尾壽命試驗。試驗采用CGS3-8型高溫高速微型軸承試驗機，試驗轉速58 000 r/min，軸向載荷5 N。在試驗的過程中測量軸承的振動加速度值，根據客戶的要求，試驗軸承的振動加速度值超過1.5 g時，就認定軸承已失效，并停止試驗。試驗中投入了7套軸承，共失效5套，軸承的壽命數據如表2所示，其中樣品1和樣品2均因為振動加速度值超過閾限停止試驗，樣品6、樣品7均未發(fā)生失效，所注壽命時間為截尾時間。

表2 軸承試驗數據Tab.2 Test data of the ball bearings

為預測樣品3、樣品4、樣品5、樣品6、樣品7的局部失效發(fā)生時間，截取這5套軸承失效前600 h的振動數據，如圖2、圖3、圖4所示；截取樣品6和樣品7截尾時間前600 h的振動數據，如圖5、圖6所示。

圖2 樣品3振動退化數據Fig.2 The vibration data of Sample 3

圖3 樣品4振動退化數據Fig.3 The vibration data of Sample 4

圖4 樣品5振動退化數據Fig.4 The vibration data of Sample 5

圖5 樣品6振動退化數據Fig.5 The vibration data of Sample 6

圖6 樣品7振動退化數據Fig.6 The vibration data of Sample 7

依據文中局部失效和完全失效的定義可判定樣品1和樣品2失效類型為局部失效，而樣品3、樣品4、樣品5在發(fā)生失效時振動未超過閾值，因此可判定樣品3、樣品4、樣品5的失效類型是完全失效。

3.1 軸承完全失效的壽命分布計算

由試驗數據可知樣品3、樣品4、樣品5為完全失效，失效時間分別為2 472 h、2 506 h、3 382 h，而樣品1在1 313 h未發(fā)生完全失效，樣品2在2 288 h未發(fā)生完全失效，樣品6和樣品7在4 000 h均未發(fā)生失效。結合數據類型，可將滾動軸承完全失效的壽命數據看作隨機截尾試驗，由參考文獻[22]可知，壽命分布參數β和η可通過求解方程(21)獲得。

(21)

式中:β和η分別是Weibull分布的形狀參數和尺度參數;ti是失效數據;則λi=1；如果ti是截尾數據，則λi=0。

結合參考文獻[23]對小樣本下定時截尾試驗壽命數據的處理方法，應用Bootstrap方法重復抽樣5 000次。再利用2.2中的壽命分布參數估計方法對Weibull分布參數進行估計，得到軸承完全失效的壽命曲線如圖7所示。

圖7 軸承完全失效壽命分布曲線

Fig.7 Life distribution curve of bearings under complete failure

3.2 軸承局部失效的壽命分布計算

這里選用指數形式的退化軌跡y=a·exp(bt)來擬合軸承的振動數據。計算所得軸承振動值達到失效閾值的時間即軸承的局部失效壽命時間，如表3所示，并采用2.2.1提出的方法進行計算，得到軸承局部失效的壽命曲線如圖8所示。

表3 軸承局部失效的壽命數據Tab.3 Lifetime data of the ball bearings under partial failure

圖8 軸承局部失效壽命曲線Fig.8 Life distribution curve of bearings under partial failure

3.3 競爭失效下軸承可靠性評估

4 結果與討論

4.1 基于競爭失效評估結果的比較

文中2.3提出的模型充分考慮了不同失效模式之間的相關性，不同失效模式競爭下的軸承可靠度計算結果如表4和圖10所示。

圖9 競爭失效的軸承可靠性曲線

Fig.9 Reliability curve of bearings under competing failure modes

從圖10中可以看出不考慮失效模式之間的相關性會使得可靠性的評估結果偏低，結合可靠性試驗，該型號7套軸承在試驗前1 300 h沒有發(fā)生失效，由表4的計算結果可知利用Copula函數求得的可靠性評估結果更符合試驗情況。

4.2 軸承可靠性評估結果的分析

圖11截取了前1 500 h完全失效與局部失效的概率進行對比，發(fā)現在軸承運轉的過程中完全失效與局部失效都有存在。在1 100 h之前，軸承發(fā)生局部失效的可能性相對較大，而在1 100 h之后，軸承發(fā)生完全失效的可能性相對較大。由此可知，在軸承運行的早期，軸承更容易發(fā)生局部失效，需要進一步尋找引起失效的原因，提高軸承的早期的可靠度；在1 100 h后，軸承的失效率主要由完全失效引起的，需要進一步尋找軸承完全失效的具體原因，以提高軸承整體的可靠度。

表4 可靠度計算結果對比Tab.4 The comparison of the reliability assessment

圖10 不同計算方法的可靠性結果Fig.10 Reliability assessment result calculated by different ways

圖11 完全失效與局部失效失效率對比

Fig.11 The comparison of the failure rate under complete failure and partial failure

5 結 論

(1) 軸承的失效是多種失效模式競爭失效的結果，在對軸承的可靠性評估中需要考慮多種失效模式對軸承可靠度的影響以及不同失效模式之間的相關性。

(2) 針對小樣本試驗數據，運用Bootstrap方法構造試驗樣本的先驗信息，采用Bayes的方法估計出相應的后驗結果，并與原始試驗數據的計算結果融合，獲得小樣本下分布參數的估計結果。

(3) 應用Gumbel Copula函數模型，實現對滾動軸承競爭失效可靠度的評估。通過試驗數據點的經驗Kendall相關秩估計Copula函數的相關參數，計算過程簡便，易于工程應用。

(4) 通過對滾動軸承實例的分析，發(fā)現如果不考慮失效模式之間的相關性可能會低估軸承運轉的可靠性。

(5) 對滾動軸承競爭失效下的可靠性分析有助于提高滾動軸承可靠性評估的準確度，為制定提高軸承可靠性的相關決策提供依據。

[1] DONG S J, YIN S R, TANG B P, et al. Bearing degradation process prediction based on the support vector machine and Markov model[J]. Shock and Vibration, 2014(1): 1-15.

[2] SUTRISNO E, OH H, VASAN A S S, et al. Estimation of remaining useful life of ball bearings using data driven methodologies[C]. IEEE Conference on Prognostics & Health Management, 2012: 1-7.

[3] ZHANG B, ZHANG L J, XU J W. Degradation feature selection for remaining useful life prediction of rolling element bearings[J]. Quality and Reliability Engineering International, 2016, 32(3): 547-554.

[4] 王華偉，高軍，吳海橋. 基于競爭失效的航空發(fā)動機剩余壽命預測[J].機械工程學報，2014，50(6): 197-205.

WANG Huawei, GAO Jun, WU Haiqiao. Residual remaining life prediction based on competing failures for aircraft engines[J]. Journal of Mechanical Engineering, 2014, 50(6): 197-205.

[5] 唐家銀，何平，梁紅琴，等. 多故障模式高長壽命產品相關性失效的綜合可靠性評估[J].機械工程學報，2013，49(12): 176-182.

TANG Jiayin, HE Ping, LIANG Hongqin, et al. Comprehensive reliability assessment of long-life products with correlated multiple failure modes[J]. Journal of Mechanical Engineering, 2013, 49(12): 176-182.

[6] 常春波，曾建潮.δ沖擊條件下相關性競爭失效過程的系統可靠性建模[J]. 振動與沖擊，2015，34(8): 203-208.

CHANG Chunbo, ZENG Jianchao. Reliability modeling for dependent competing failure processes underδ-shock[J]. Journal of Vibration and Shock, 2015, 34(8): 203-208.

[7] BOCCHETTI D, GIORGIO M, GUIDA M, et al. A competing risk model for the reliability of cylinder liners in marine Diesel engines[J]. Reliability Engineering and System Safety, 2009, 94(8): 1299-1307.

[8] 蔣喜，劉宏昭，劉麗蘭，等. 基于偽壽命分布的電主軸極小子樣可靠性研究[J]. 振動與沖擊，2013，32(19): 80-85.

JIANG Xi, LIU Hongzhao, LIU Lilan, et al. Extremely small-scale sample’s reliability of an electric spindle based on distribution of false lifetime[J]. Journal of Vibration and Shock, 2013, 32(19): 80-85.

[9] XIONG J, SHENOI R A, GAO Z. Small sample theory for reliability design[J]. The Journal of Strain Analysis for Engineering Design, 2002, 37(1): 87-92.

[10] PICHENY V, KIM N H, HAFTKA R T. Application of the bootstrap method in conservative estimation of reliability with limited samples[J]. Structural & Multidisciplinary Optimization, 2010, 41(2): 205-217.

[11] KURZ D, LEWITSCHNIG H, PILZ J. Advanced Bayesian estimation of Weibull early life failure distributions[J]. Quality and Reliability Engineering International, 2014, 30(3): 363-373.

[12] ERTO P. New practical bayes estimators for the 2-parameter Weibull distribution[J]. IEEE Transactions on Reliability, 1982, 31(2): 194-197.

[13] 金光. 一種綜合性能與壽命試驗數據的Bayes-Bootstrap方法[J]. 宇航學報，2007，28(3): 731-734.

JIN Guang. A Bayes-Bootstrap method synthesizing performance and life data[J]. Journal of Astronautics, 2007, 28(3): 731-734.

[14] RAUSAND M. 系統可靠性理論：模型、統計方法及應用[M].郭強，王秋芳，劉樹林譯. 2版.北京：國防工業(yè)出版社，2011.

[15] 肖支才，高華明，丁通，等. 自助法在小樣本數據均值估計中的應用[J].海軍航空工程學院學報，2009，24(5): 563-572.

XIAO Zhicai, GAO Huaming, DING Tong, et al. Research on data anlysis of small sample based on bootstrap method[J]. Journal of Naval Aeronautical and Astronautical University, 2009, 24(5): 563-572.

[16] EFRON B, TIBSHIRANI R. An introduction to the Bootstrap[M]. New York: Chapman hall, 1993.

[17] LU C J, MEEKER W Q. Using degradation measures to estimation time-to-failure distribution[J]. Technometrics, 1993, 35(2): 161-174.

[18] 趙建印，劉芳，孫權，等. 基于性能退化數據的金屬化膜電容器可靠性評估[J]. 電子學報，2005，33(2): 378-381.

ZHAO Jianyin, LIU Fang, SUN Quan, et al. Reliability estimate of metallized-film pulse capacitor from degradation data[J]. Acta Electronica Sinca, 2005, 33(2): 378-381.

[19] NELSEN R B. An introduction to Copulas. (Second Edition)[M]. Springer, 2005.

[20] LEE L. Multivariate distributions having Weibull properties[J].Journal of Multivariate Analysis, 1979, 9(2):267-277.

[21] 唐家銀，何平，陳崇雙. 相關性失效機械系統的可靠性分析方法[M].北京：國防工業(yè)出版社，2014.

[22] 陳家鼎. 隨機截尾情形下Weibull分布參數的最大似然估計的相合性[J]. 應用概率統計，1989，5(3): 226-223.

CHEN Jiading. On the consistency of maximum likelihood estimators based on randomly censored data in the Weibull distribution case[J]. Chinese Journal of Applied Probability and Statistics, 1989, 5(3): 226-223.

[23] 夏新濤，章寶明，徐永智. 滾動軸承性能與可靠性乏信息變異過程評估[M].北京：科學出版社，2013.

Reliabilityassessmentofbearingsbasedoncompetingfailureundersmallsampledata

QIN Luosheng, CHEN Xiaoyang, SHEN Xuejin

(Department of Mechanical Automation Engineering, Shanghai University, Shanghai 200072, China)

Bearings are important parts in machinery products. Their performance and life are closely related to operating lives of mechanical systems. It is necessary to consider the effects of different failure modes on rolling bearings’ reliability. Here, aiming at the full failure data in bearings’ reliability tests, the prior distributions of bearing life distribution parameters were established with Bootstrap method. Then, the corresponding posterior distributions were estimated using Bayes method. The bearing life distribution parameters were obtained through the posterior expectation reduction. The life distribution model of bearing local failure was gained through further analyzing bearings’ vibration performance degradation data. Copula function was used to analyze comprehensively the life distribution model of bearing full failure and that of bearing local failure. The relative parameters of Copula function were estimated with the experience Kendall relative ranks of test data. Finally, the reliability assessment results of bearings under competing failure were achieved. The results were helpful to find defects in bearing design and improve the bearing reliability.

bearings; Bootstrap method;Bayes method; performance degradation; Copula function; reliability

國家“十二五”規(guī)劃項目(D.50-0109-15-001,D.71-0109-14-197)

2016-05-30 修改稿收到日期：2016-09-12

秦犖晟 男，碩士生,1992年生

沈雪瑾 女,博士,教授,1963年生

TH133.33;TB114.3

A

10.13465/j.cnki.jvs.2017.23.036