基于Bootstrap的小樣本可靠性評(píng)估方法

張 震, 劉儉輝, 趙 成, 剡昌鋒

(1. 蘭州理工大學(xué) 機(jī)電工程學(xué)院, 甘肅 蘭州 730050; 2. 林德液壓(中國(guó))有限公司, 山東 濰坊 261061)

目前，各種設(shè)備朝著精密化、高速化、低能耗的方向發(fā)展，對(duì)各種機(jī)械零部件的可靠性要求也越來(lái)越高，導(dǎo)致可靠性評(píng)估成為設(shè)備研發(fā)過(guò)程中不可或缺的重要組成部分.然而，由于研發(fā)成本、實(shí)驗(yàn)周期、實(shí)驗(yàn)條件等各方面因素的限制，有些產(chǎn)品無(wú)法進(jìn)行大量實(shí)驗(yàn)，投試數(shù)量非常小，只能針對(duì)小樣本數(shù)據(jù)進(jìn)行可靠性評(píng)估.對(duì)小樣本數(shù)據(jù)進(jìn)行分析得到設(shè)備壽命分布是可靠性評(píng)估的關(guān)鍵，因此提出一種優(yōu)良的小樣本數(shù)據(jù)參數(shù)求解方法是解決可靠性評(píng)估的前提.

針對(duì)小樣本問(wèn)題，Bayes法[1]、Bootstrap法[2]、Bayes Bootstrap法[3]、Monte Carlo法[4]以及Grey Model法[5]是實(shí)際工程中的常用方法.國(guó)內(nèi)外學(xué)者針對(duì)小樣本可靠性評(píng)估問(wèn)題進(jìn)行了大量的研究[6-13].Yang等[7]根據(jù)Bayes理論建立了小樣本的制動(dòng)系統(tǒng)防滑閥、重型數(shù)控鏜銑床、航天器可靠性評(píng)估模型.Liu等[9]基于Bootstrap的思想求解了航空電動(dòng)滑環(huán)壽命分布參數(shù)，結(jié)合傳統(tǒng)的可靠性預(yù)測(cè)方法，得到了可靠性指標(biāo).Chen等[10]將Bayes方法融入Bootstrap的思想中，解決了小樣本可靠性模型在實(shí)際應(yīng)用中可能存在偏差的問(wèn)題.Zhou等[12]將灰色理論運(yùn)用到小樣本故障數(shù)據(jù)的預(yù)測(cè)中，分別從模糊理論、白化方程本身、Bootstrap法對(duì)灰色模型進(jìn)行改進(jìn)，提高了故障預(yù)測(cè)的可靠性.除此以外，有些學(xué)者還提出了其他小樣本數(shù)據(jù)處理方法.Zou[14]提出了一種基于最小二乘支持向量機(jī)(LSSVM)的小數(shù)據(jù)樣本壽命預(yù)測(cè)方法，處理指數(shù)分布的小數(shù)據(jù)樣本時(shí)，有很好的推廣性和壽命預(yù)測(cè)的準(zhǔn)確性.任博等[15]針對(duì)故障數(shù)據(jù)稀少難以準(zhǔn)確確定數(shù)據(jù)分布類型的問(wèn)題，提出了Kriging方法建立部件故障不確定性的精確描述模型.譚巍等[16]通過(guò)融合性能退化可靠性理論和隨機(jī)過(guò)程，分析了航空發(fā)動(dòng)機(jī)的性能退化過(guò)程，并運(yùn)用Weibull分布提出了基于隨機(jī)過(guò)程的航空發(fā)動(dòng)機(jī)性能退化過(guò)程中的可靠性評(píng)估模型.

小樣本可靠性研究常采用定時(shí)截尾實(shí)驗(yàn)方法[17-19].馮自立等[17]對(duì)Weibull分布下的小樣本定時(shí)截尾實(shí)驗(yàn)數(shù)據(jù)的處理方法進(jìn)行研究，發(fā)現(xiàn)形狀參數(shù)的矩估計(jì)量中值與其真值的比值僅取決于樣本結(jié)構(gòu)，不受分布母體的影響，據(jù)此提出了中值無(wú)偏矩估計(jì)法并借助Bootstrap法對(duì)該算法進(jìn)行了改進(jìn).

中值無(wú)偏矩估計(jì)法對(duì)于定時(shí)截尾數(shù)據(jù)的處理適用性很好，但不能完整處理實(shí)驗(yàn)數(shù)據(jù).因此，本文在極大似然估計(jì)法的基礎(chǔ)上引入Bootstrap的思想，對(duì)樣本進(jìn)行重抽樣，使小樣本問(wèn)題轉(zhuǎn)化為大樣本問(wèn)題，并采用核密度估計(jì)法(Kernel density estimation)求解得到參數(shù)的概率密度函數(shù)，進(jìn)而得到一種適用于小樣本定時(shí)截尾數(shù)據(jù)和完整實(shí)驗(yàn)數(shù)據(jù)的可靠性評(píng)估方法，并運(yùn)用Monte Carlo法進(jìn)行了驗(yàn)證.

1 基本理論

1.1 Weibull分布

疲勞壽命可靠性研究最多的是正態(tài)分布和Weibull分布，其中，工程中大多數(shù)機(jī)械產(chǎn)品的壽命都滿足Weibull分布，其分布概率密度函數(shù)的表達(dá)式為

(1)

式中：β為形狀參數(shù)，它的大小決定了分布曲線的離散程度，β值越大函數(shù)分布越集中，在可靠性評(píng)估中，產(chǎn)品的可靠壽命和特征壽命隨著形狀參數(shù)的增大而提高，β值一般與材料有關(guān)；η為尺度參數(shù)，表示產(chǎn)品可靠度R=0.367 9時(shí)對(duì)應(yīng)的壽命值，即t0.367 9=η，能夠放大或縮小坐標(biāo)尺度，但不能影響Weibull分布的形狀；γ為位置參數(shù)，即最小壽命，其決定了分布曲線的起始位置；t為產(chǎn)品的服役時(shí)間.

當(dāng)γ=0時(shí)，式(1)就退化為兩參數(shù)Weibull分布，兩參數(shù)Weibull分布的累計(jì)分布函數(shù)、可靠度函數(shù)、失效率函數(shù)表達(dá)式分別為

1.2 兩參數(shù)Weibull分布極大似然估計(jì)

形狀參數(shù)β與尺度參數(shù)η的威布爾分布概率密度函數(shù)為

(5)

對(duì)n個(gè)產(chǎn)品進(jìn)行定時(shí)截尾實(shí)驗(yàn)，其中有r個(gè)產(chǎn)品失效，失效數(shù)據(jù)為(t1,t2,…,tr)，截尾樣本的個(gè)數(shù)q=n-r，截尾時(shí)間為tc，則截尾數(shù)據(jù)為(tr+1,tr+2,…,tn)=(tc,tc,…,tc)，完整失效數(shù)據(jù)是一種特殊情況，數(shù)據(jù)的似然函數(shù)為

(6)

似然函數(shù)取對(duì)數(shù)為

(7)

1.3 Bootstrap抽樣

(10)

(11)

重復(fù)以上步驟，能夠得到多個(gè)θ(F)，就可以統(tǒng)計(jì)求出未知參數(shù)的θ分布和特征值.

2 基于Bootstrap的極大似然估計(jì)改進(jìn)

2.1 方法原理

針對(duì)極大似然估計(jì)在對(duì)失效數(shù)據(jù)進(jìn)行Weibull參數(shù)估計(jì)時(shí)，樣本量不足造成的誤差過(guò)大且不能得到參數(shù)置信區(qū)間的問(wèn)題，提出運(yùn)用Bootstrap抽樣方法對(duì)失效數(shù)據(jù)進(jìn)行再生，擴(kuò)大樣本信息，基于此再運(yùn)用極大似然估計(jì)求解Weibull分布的參數(shù)β和η，得到多組θ的估計(jì)值(為方便敘述，用θ統(tǒng)一表示β和η)，這是求解概率密度函數(shù)的基礎(chǔ).為避免概率密度求解主觀假設(shè)造成的誤差，采用核密度估計(jì)來(lái)求解θ的概率密度函數(shù)，直接從抽樣結(jié)果的樣本數(shù)據(jù)出發(fā)擬合.用f(θ)統(tǒng)一表示f(β)和f(η)，即β和η的概率密度函數(shù)為

(12)

式中：n為樣本個(gè)數(shù)；h為帶寬；K(·)為核函數(shù)，需滿足：

(13)

h的取值會(huì)直接影響概率密度函數(shù)的光滑程度.取值較大，則有更多的樣本點(diǎn)對(duì)概率密度估計(jì)值產(chǎn)生影響，生成曲線更光滑，但也丟失了部分?jǐn)?shù)據(jù)信息；取值較小，則生成曲線為不光滑折線，但能反映每個(gè)數(shù)據(jù)所包含的信息.因此，選擇合適的帶寬是核密度估計(jì)的保證.定義h的誤差函數(shù)為

(14)

(15)

式中：σ為樣本數(shù)據(jù)的標(biāo)準(zhǔn)差.

將hMISE代入式(12)可得到各樣本點(diǎn)的概率密度估計(jì)值，但卻不能得到總體的概率密度函數(shù).為了得到概率密度函數(shù)，采用混合高斯模型對(duì)概率密度估計(jì)值進(jìn)行擬合，混合高斯模型的表達(dá)式為

(16)

式中：N為所用高斯分布的個(gè)數(shù);ai、bi、ci為第i個(gè)高斯分布的參數(shù).

當(dāng)給定顯著性水平α?xí)r，參數(shù)θ估計(jì)區(qū)間為

[θL,θU]=[θα/2,θ1-α/2]

(17)

式中：θL和θU為參數(shù)的下限和上限.

則：

定義可靠度R(t)置信區(qū)間為[RL(t),RU(t)]，下限RL(t)為

min{exp[-(t/ηL)βL],exp[-(t/ηL)βU],…,

exp[-(t/ηU)βL],exp[-(t/ηU)βU]}

(20)

上限RU(t)為

max{exp[-(t/ηL)βL],exp[-(t/ηL)βU],…,

exp[-(t/ηU)βL],exp[-(t/ηU)βU]}

(21)

本文基于Bootstrap法和極大似然估計(jì)法，通過(guò)Bootstrap重抽樣擴(kuò)大了原樣本，然后采用極大似然估計(jì)求解樣本參數(shù)，得到參數(shù)的多組估計(jì)值，采用核密度估計(jì)求解參數(shù)的概率密度函數(shù)，在此基礎(chǔ)上進(jìn)一步得到參數(shù)的置信區(qū)間以及可靠度的置信區(qū)間.

2.2 Monte Carlo驗(yàn)證

令β=2.5，η=200，生成一組符合兩參數(shù)Weibull分布的數(shù)據(jù)作為失效時(shí)間，從中選取一定量的數(shù)據(jù)作為樣本，樣本數(shù)為10 000組，得到2種方法下參數(shù)的計(jì)算結(jié)果，如表1所列.為闡述方便，稱本文提出的方法為B-mle法.

將表1中的計(jì)算結(jié)果轉(zhuǎn)化為圖1.

表1 B-mle法和Mle法參數(shù)結(jié)果對(duì)比

圖1 不同樣本量的計(jì)算結(jié)果對(duì)比Fig.1 Comparison of calculation results of different sample sizes

由表1和圖1可以看出，隨著樣本量的增加，B-mle法與Mle法所求參數(shù)的誤差均減小.樣本量越大，計(jì)算結(jié)果與真值越接近.在相同樣本量的情況下，B-mle法比Mle法求得參數(shù)的計(jì)算精度更高.在樣本量為8時(shí)，相較于Mle法求得的β值誤差，B-mle法減小了22.51%，η值誤差減小了2.61%.在90%置信水平下，參數(shù)的置信區(qū)間如表2所列.

表2 90%置信水平參數(shù)置信區(qū)間

由表2可以看出，B-mle方法不僅可以得到參數(shù)的點(diǎn)估計(jì)，還可以得到給定置信水平下的置信區(qū)間，而極大似然估計(jì)法只能得到參數(shù)的點(diǎn)估計(jì)，證明B-mle法優(yōu)于Mle法.

2.3 2種方法評(píng)估效果對(duì)比

為了對(duì)比2種方法的評(píng)估效果，在時(shí)間[0,400]的范圍內(nèi)每間隔1取作一個(gè)樣本點(diǎn)，計(jì)算2種方法的可靠度估計(jì)值與真值的均方根誤差RMSE，計(jì)算公式為

(22)

由式(22)可以看出，RMSE的值越小，估計(jì)值與真值的偏差越小，評(píng)估效果越好.計(jì)算不同樣本下的RMSE，結(jié)果如表3所列.

表3 2種方法可靠度估計(jì)值的RMSE對(duì)比

由表3可以看出，樣本量相同的情況下B-mle法得到的可靠度估計(jì)誤差總小于Mle法的，再次證明了B-mle法的正確性與優(yōu)越性.

3 案例分析

3.1 柱塞泵實(shí)驗(yàn)數(shù)據(jù)

以某型液壓系統(tǒng)柱塞泵為研究對(duì)象，對(duì)其在現(xiàn)場(chǎng)使用情況進(jìn)行記錄，投入使用的該型號(hào)柱塞泵總體為N=14，柱塞泵失效數(shù)據(jù)X=[15,25,29,40,60,35,50,90,110,130,150,160,170,180]，以上數(shù)據(jù)來(lái)自于文獻(xiàn)[15]，對(duì)該數(shù)據(jù)從小到大進(jìn)行排序，如表4所列.

表4 柱塞泵故障數(shù)據(jù)

3.2 經(jīng)驗(yàn)可靠度評(píng)估

對(duì)故障數(shù)據(jù)進(jìn)行經(jīng)驗(yàn)可靠性評(píng)估，根據(jù)Johnson方法[20]，將可靠性經(jīng)驗(yàn)值表示為中位秩，即

(23)

式中：r(ti)為經(jīng)驗(yàn)可靠度函數(shù);ti為X中第i個(gè)壽命數(shù)據(jù);i為壽命數(shù)據(jù)ti的序號(hào);n為X的壽命數(shù)據(jù)個(gè)數(shù).

由式(23)得到故障數(shù)據(jù)的經(jīng)驗(yàn)分布函數(shù)，并用兩參數(shù)Weibull分布去描述故障數(shù)據(jù)，如圖2所示.

圖2 柱塞泵故障數(shù)據(jù)經(jīng)驗(yàn)分布函數(shù)Fig.2 Empirical distribution function of piston pump failure data

由圖2可以看出，兩參數(shù)Weibull分布能夠?qū)?jīng)驗(yàn)分布函數(shù)進(jìn)行很好的描述，對(duì)表4中的故障數(shù)據(jù)用極大似然估計(jì)求解形狀參數(shù)η與尺度參數(shù)β，結(jié)果分別是99.029與1.556.

3.3 數(shù)據(jù)的K-S檢驗(yàn)

Kolmogorov-Smirnov檢驗(yàn)法用來(lái)檢驗(yàn)樣本是否服從指定的分布.在進(jìn)行K-S檢驗(yàn)時(shí)，首先將n個(gè)數(shù)據(jù)由小到大排列，根據(jù)假設(shè)的分布計(jì)算每個(gè)數(shù)據(jù)對(duì)應(yīng)的分布函數(shù)G(ti).提出假設(shè),H0表示數(shù)據(jù)符合假設(shè)的分布,H1表示數(shù)據(jù)不符合假設(shè)的分布.接著通過(guò)樣本的經(jīng)驗(yàn)分布函數(shù)r(ti)和假設(shè)的分布函數(shù)G(t)構(gòu)造統(tǒng)計(jì)量D，即

D=max|r(t)-G(t)|

(24)

通過(guò)查閱K-S檢定表，由樣本大小n和顯著性水平α，找出臨界值D(n,α)，當(dāng)D>D(n,α)時(shí)，拒絕H0.

借助MATLAB中kstest函數(shù)來(lái)進(jìn)行K-S檢驗(yàn)，檢驗(yàn)值和臨界值如表5所列.

表5 K-S檢驗(yàn)

由表5可以看出，檢驗(yàn)值0.141小于臨界值0.349，接受原假設(shè)H0，數(shù)據(jù)符合形狀參數(shù)η與尺度參數(shù)β為99.029與1.556的Weibull分布.但此時(shí)參數(shù)是基于小樣本數(shù)據(jù)計(jì)算的，有可能存在誤差過(guò)大的問(wèn)題.

本文應(yīng)用Bootstrap、極大似然估計(jì)法對(duì)所求尺度參數(shù)η進(jìn)行統(tǒng)計(jì)，直方圖如圖3所示.通過(guò)核密度估計(jì)與混合高斯模型擬合得到尺度參數(shù)η的概率密度函數(shù),如圖4所示.

圖3 尺度參數(shù)η頻率直方圖Fig.3 Frequency histogram of scale parameters η

圖4 尺度參數(shù)η混合高斯分布擬合Fig.4 Fitting of mixed Gaussian distribution of scale parameters η

由圖4可以看出，4個(gè)混合高斯分布的擬合效果最好，故選用Gauss-4作為密度函數(shù)進(jìn)行后續(xù)計(jì)算.

形狀參數(shù)β的概率密度函數(shù)求解過(guò)程與尺度參數(shù)η類似,頻率直方圖和混合高斯分布擬合如圖5、圖6所示.

圖5 形狀參數(shù)β頻率直方圖

圖6 形狀參數(shù)β混合高斯分布擬合Fig.6 Fitting of mixed Gaussian distribution of shape parameters β

由圖6可以看出，Gauss-4作為密度函數(shù)擬合效果較好.在90%置信度下，尺度參數(shù)η的置信區(qū)間為[71.722, 127.386]，形狀參數(shù)β的置信區(qū)間為[1.418, 2.729].尺度參數(shù)η的期望值為99.328，形狀參數(shù)β的期望值為1.723.當(dāng)給定時(shí)間時(shí)，將參數(shù)值代入式(16)和式(17)即可得到可靠度的置信區(qū)間.

4 結(jié)論

本文通過(guò)Bootstrap抽樣擴(kuò)充小數(shù)據(jù)樣本，在此基礎(chǔ)上采用極大似然估計(jì)、核密度估計(jì)等方法對(duì)數(shù)據(jù)進(jìn)行處理，提出了適用于小樣本完整數(shù)據(jù)和非完整數(shù)據(jù)的B-mle方法，并得出以下結(jié)論：

1) 通過(guò)Monte Carlo模擬仿真，B-mle法求解的參數(shù)相較于Mle法，誤差減小了22.51%和2.61%，證明了B-mle法的正確性與優(yōu)越性.

2) 樣本量為8時(shí)，B-mle法與Mle法的均方根誤差RMSE分別為0.100與0.067；樣本量為10時(shí)，B-mle法與Mle法的均方根誤差RMSE分別為0.019與0.016.再次證明了B-mle法比Mle法在處理小樣本數(shù)據(jù)時(shí)具有更高的精度.

致謝：本文得到蘭州理工大學(xué)紅柳優(yōu)青人才資助項(xiàng)目(2020062001)的資助，在此表示感謝.