基于多層自助最大熵法的可靠性評估

夏新濤，葉亮，李云飛，常振

（河南科技大學機電工程學院，河南洛陽471003）

基于多層自助最大熵法的可靠性評估

夏新濤，葉亮，李云飛，常振

（河南科技大學機電工程學院，河南洛陽471003）

在小樣本無任何先驗信息的條件下，提出多層自助最大熵評估模型分析機械產(chǎn)品壽命可靠性。用自助法對當前無失效數(shù)據(jù)樣本進行再抽樣，獲得足夠的樣本數(shù)據(jù)?；谧畲箪胤?，改變抽樣個數(shù)得到不同的拉格朗日乘子。再次運用自助法對拉格朗日乘子的小樣本數(shù)據(jù)進行再抽樣，基于最大熵法獲得拉格朗日乘子的區(qū)間估計。對每個拉格朗日乘子的上下限進行排列組合，得到多個概率密度函數(shù)和可靠性函數(shù)，運用最小不確定性原理得到可靠性函數(shù)的區(qū)間估計。試驗研究表明，多層自助最大熵評估模型可以有效地解決概率分布已知或未知的小樣本無失效數(shù)據(jù)的可靠性評估問題。

系統(tǒng)評估與可行性分析；可靠性評估；多層自助最大熵法；乏信息；無失效數(shù)據(jù)；拉格朗日乘子

0 引言

目前，無失效數(shù)據(jù)的可靠性評估方法主要有經(jīng)典統(tǒng)計方法和貝葉斯統(tǒng)計方法。經(jīng)典統(tǒng)計理論認為，通過實驗獲得的失效數(shù)據(jù)越多，對產(chǎn)品進行可靠性評估的結(jié)果越準確。但是，對于許多高可靠性或危險的試驗來說，很難獲得失效數(shù)據(jù)，例如高端武器和航天航空飛行器的試驗等。因此，無失效數(shù)據(jù)的可靠性評估問題日益引起學術界的關注。在機械產(chǎn)品可靠性試驗中，獲得的數(shù)據(jù)通常有3類，即失效數(shù)據(jù)、不完全失效數(shù)據(jù)、無失效數(shù)據(jù)。目前，對失效數(shù)據(jù)的研究已經(jīng)達到比較高的水平。然而對無失效數(shù)據(jù)的研究還不夠深入，因此這也是長期以來國內(nèi)外研究的熱點與難點。

關于無失效數(shù)據(jù)的可靠性分析，很多研究工作已經(jīng)取得了相應的成果。如傅惠民等［1］提出一種Weibull分布定時無失效數(shù)據(jù)可靠性分析方法，在形狀參數(shù)下限已知的情況下，給出了可靠度和使用壽命的單側(cè)置信下限，從而能夠根據(jù)定時無失效數(shù)據(jù)對產(chǎn)品進行高置信水平的可靠性評定；王憑慧等［2］在高置信水平下提出了一種衛(wèi)星推力器無失效數(shù)據(jù)的可靠性分析方法；李寶盛等［3］求出了工作壽命可靠度置信區(qū)間的解析公式，計算和分析了不同試驗樣本工作壽命可靠度的置信限，得出了對可靠性評估有指導意義的結(jié)論；蒲星［4］運用Bayes方法對無失效數(shù)據(jù)問題進行了有效的可靠性研究，分析了無失效數(shù)據(jù)產(chǎn)生的原因，并提出了運用黃金分割系數(shù)求失效概率的方法。但召江等［5］在形狀參數(shù)先驗分布分別為均勻分布與擬合出的概率分布時，利用無失效試驗數(shù)據(jù)得出失效率和形狀參數(shù)的Bayes估計，進而計算出Weibull分布的特征壽命；Jiang等［6］提出了運用修整極大似然方法估計威布爾分布的參數(shù)。在國外，Nam等［7］以深溝球軸承為例，研究了軸承在高溫條件下的壽命可靠性；Bailey［8］基于無失效數(shù)據(jù)，建立了二項概率分布的預測模型；Kwon［9-10］基于Bayes方法，分別針對壽命服從Weibull分布和對數(shù)正態(tài)分布的產(chǎn)品，提出了無失效數(shù)據(jù)樣本的可靠性驗證方法。

對于未知概率分布的小樣本數(shù)據(jù)的可靠性評估問題，現(xiàn)有的可靠性評估方法難以解決，到目前為止還是一個重要的科學技術難題。鑒于此，夏新濤等［11］以灰色系統(tǒng)的基本原理為依據(jù)，構(gòu)建出特殊的灰相似空間，允許概率分布未知和數(shù)據(jù)個數(shù)很少（n≥4）。賈波等［12］在仿真實現(xiàn)過程中，發(fā)現(xiàn)最大熵方法存在溢出的問題，通過變量變換法成功解決了此問題。夏新濤等［13］將灰自助原理融入泊松過程，提出灰自助泊松方法，以預測滾動軸承振動性能可靠性的變異過程?；谧畲箪胤?，胡松偉［14］提出了火工品類關鍵單元可靠性驗證的LQB型方案。在國外，Yeom等［15］提出改進的最大熵抽樣方法，從而使大多數(shù)抽樣點落在一個用復雜的非線性函數(shù)表示的可行域內(nèi)；Young等［16］基于最大熵原理，求解第一結(jié)尾時間序列的可靠性邊界。

根據(jù)測量理論和統(tǒng)計理論，任何參數(shù)估計都伴隨不確定性，可靠性評估中的不確定性可以用可靠性估計區(qū)間來表示。現(xiàn)有可靠性理論把拉格朗日乘子當作常數(shù)，而Bayes理論認為統(tǒng)計學上的參數(shù)可當作變量，可以對參數(shù)進行區(qū)間估計，因而也能對可靠性函數(shù)做出區(qū)間估計。因此本文對于只有無失效數(shù)據(jù)而沒有概率分布任何先驗信息條件下機械產(chǎn)品的可靠性評估問題，提出了多層自助最大熵評估模型來分析機械產(chǎn)品壽命可靠性。具體過程如下:首先，運用自助法對小數(shù)據(jù)樣本進行等概率可放回再抽樣，獲得大量樣本數(shù)據(jù)。改變每次抽樣個數(shù)，基于最大熵法，可得到多個不同的拉格朗日乘子。然后，再次利用自助法對所有的拉格朗日乘子進行再抽樣，多次使用最大熵法獲得每個拉格朗日乘子的區(qū)間估計。最后，對各個拉格朗日乘子的上下限進行排列組合，得到多個概率密度函數(shù)和可靠性函數(shù)，根據(jù)最小不確定性原理獲得可靠性函數(shù)的區(qū)間估計。通過試驗證明運用該方法能夠?qū)煽啃院瘮?shù)做出合理的區(qū)間估計，分析結(jié)果真實可信。

1 多層自助最大熵法原理

1.1 自助法原理

假設通過試驗獲得一組無失效數(shù)據(jù)序列X，用向量表示為

式中:xl為第l個無失效數(shù)據(jù)；n為無失效數(shù)據(jù)的個數(shù)。

從無失效數(shù)據(jù)樣本X中等概率可放回地抽樣，每次抽取t個數(shù)據(jù)，連續(xù)重復抽取B次，可以得到自助樣本Xt為

式中:t為每次抽取原始樣本數(shù)據(jù)的個數(shù)；k為生成自助樣本的數(shù)據(jù)序號；xt（k）為每次抽取t個數(shù)據(jù)生成自助樣本的第k個數(shù)據(jù)。

1.2 基于最大熵方法求解概率密度函數(shù)

最大熵方法能夠?qū)ξ粗母怕史植甲龀鲋饔^偏見為最小的最佳估計。在求解過程中，引入拉格朗日乘子，從而把概率分布問題轉(zhuǎn)化為拉格朗日乘子的求解問題。

為了敘述方便，將離散的無失效數(shù)據(jù)序列X連續(xù)化，定義最大熵的表達式為

式中:f（x）為連續(xù)化后的數(shù)據(jù)序列的概率密度函數(shù)；Inf（x）為概率密度函數(shù)的對數(shù)；S為積分區(qū)間。

通過調(diào)整f（x）可以使熵達到最大值，同時采取拉格朗日乘子法求解。設為拉格朗日函數(shù)，可得

式中:v為原點矩階數(shù)，常用v=5；mλ為第λ階原點矩；xλ為求解第λ階原點矩時f（x）的系數(shù)；cλ為第λ+1個拉格朗日乘子，λ=0，1，…，v.

其他m個拉格朗日乘子應滿足

1.3 積分區(qū)間的映射

為了使求解收斂，將無失效數(shù)據(jù)按遞增的順序排列并分成j組，畫出直方圖，同時可得到組中值zq和頻數(shù)pq，q=2，3，…，j+1.然后將直方圖擴展成j+ 2組，并令p1=pj+2.將原始數(shù)據(jù)區(qū)間S映射到區(qū)間［-e，e］中，具體如下:

令

式中:a、b為映射參數(shù)；w為所要變換的自變量；x∈［-e，e］，

由dx=dw/a可得

因此，最大熵概率分布密度函數(shù)由（6）式變換為

1.4 求解概率密度函數(shù)的上下限

現(xiàn)有可靠性理論把拉格朗日乘子當做常數(shù)，不能獲得概率密度函數(shù)的上下限。因此，本文根據(jù)Bayes原理，把拉格朗日乘子當做變量，運用多層自助最大熵方法求出概率密度函數(shù)的上下限。

改變抽樣個數(shù)k，基于自助最大熵法可得到不同的拉格朗日乘子Cλ和映射參數(shù)a、b.

式中:cλ（t）為每次抽取t個數(shù)據(jù)時的第λ+1個拉格朗日乘子系數(shù)；a（t）、b（t）為每次抽取t個數(shù)據(jù)時的區(qū)間映射參數(shù)。

1.4.1 各個拉格朗日乘子的區(qū)間估計

將各個拉格朗日乘子作為源信息樣本，再進行等概率可放回地抽樣?；谧畲箪胤椒?，可求解得各個拉格朗日乘子的區(qū)間估計。

假設顯著性水平為α，α∈［0，1］，則置信水平P為

各個拉格朗日乘子cλ在置信水平P下的下邊界值設為cλL，則有

式中:cλ0為積分變量初始值。

因此，參數(shù)的估計區(qū)間為

1.4.2 映射參數(shù)的點估計

1.4.3 概率密度的上下界函數(shù)求解

分別將拉格朗日乘子的上限或下限值，映射參數(shù)的點估計值代入（13）式中，由排列組合原理，可得2v+1條概率密度函數(shù)曲線的表達式為

式中:i為概率密度函數(shù)的序號；c0i為第i個概率密度函數(shù)的第1個拉格朗日系數(shù)，c0i=c0L或c0i=c0U；cλi為第i個概率密度函數(shù)的第λ+1個拉格朗日系數(shù)，cλi=cλL或cλi=cλU，cλU、cλL分別為拉格朗日乘子cλ的上下限。

定義概率密度估計真值函數(shù)f0（x）為

式中:f0（x）為每次抽取原始樣本數(shù)據(jù)個數(shù)n時的概率密度函數(shù)；c00、cλ0分別為每次抽取n個數(shù)據(jù)時的第1個和第λ+1個拉格朗日系數(shù)；a0、b0為每次抽取n個數(shù)據(jù)時的區(qū)間映射參數(shù)。

根據(jù)最小不確定性原理，可以從概率密度函數(shù)曲線中得到距離概率密度估計真值函數(shù)f0（x）最近的上下兩條曲線，即最大熵概率密度分布的上下界函數(shù)fU（x），fL（x）.

1.5 可靠性的上下界函數(shù)求解

1.5.1 可靠性估計真值函數(shù)

將每次抽取原始無失效數(shù)據(jù)樣本個數(shù)n時計算所得的可靠性函數(shù)作為其估計真值函數(shù)R0（x），具體計算如下:對最大熵概率密度估計真值函數(shù)f0（x）積分，得到最大熵概率分布估計真值函數(shù)F0（x）為

式中:S0為積分區(qū)間下限值。

設可靠性系數(shù)為rc，由個體無失效數(shù)據(jù)個數(shù)n可得總體的失效概率估計真值函數(shù)P0（x）為

總體可靠性估計真值函數(shù)也即最大熵可靠性估計真值函數(shù)R0（x）為

1.5.2 可靠性函數(shù)的區(qū)間估計

將最大熵概率密度分布的上下界函數(shù)fU（x）、fL（x）分別代入式（26）、（27）、（28）中，可求解出可靠性函數(shù)的上下界曲線RU（x）、RL（x）為

給定置信水平P，根據(jù)可靠性函數(shù)的真值估計曲線R0（x），得到估計真值x0.將其代入到可靠性函數(shù)的上下限中，計算出可靠性函數(shù)的區(qū)間估計為

2 試驗研究

案例1 根據(jù)文獻［17］，認為軸承的壽命概率分布為威布爾分布。對20套軸承進行壽命試驗，獲得的無失效數(shù)據(jù)如下:

X=（422，422，539，539，539，539，602，602，

770，770，770，770，847，847，847，

847，924，924，924，924）.

用自助法每次抽取20個數(shù)據(jù)，共抽取30 000次，所得數(shù)據(jù)如圖1所示。

圖1 自助法獲取軸承樣本數(shù)據(jù)Fig.1 Sample data of bearings obtained by bootstrap method

基于最大熵法計算可得:映射參數(shù)a0=0.016 8，b0=-11.974 0；拉格朗日乘子（c00，c10，c20，c30，c40，c50）=（-4.601 9，0.351 5，-1.135 0，-0.022 6，-0.025 2，-0.003 8）.

由（25）式可計算出概率密度估計真值函數(shù)f0（x），將f0（x）代入（26）式、（27）式、（28）式中，可得可靠性估計真值函數(shù)R0（x），如圖2所示。由圖2可得:R0（677）=0.994 5，與文獻［17］中運用Bayes方法計算出的R（677）=0.9762相差0.0183.這表明自助最大熵方法的可靠性估計真值結(jié)果與用Bayes方法的可靠性估計真值結(jié)果相差很小?？梢娮灾畲箪胤椒ǐ@得的可靠性估計真值函數(shù)的擬合效果跟Bayes方法的估計結(jié)果基本一致。但是自助最大熵方法對概率分布沒有要求，而且能解決小樣本數(shù)據(jù)的可靠性評估問題。因此，用自助最大熵方法獲得的可靠性估計真值函數(shù)評估機械產(chǎn)品的可靠性是有效且可行的。

改變每次抽樣個數(shù)t，共抽取30 000次，可得拉格朗日乘子如表1所示。

圖2 可靠性估計真值函數(shù)曲線Fig.2 Estimated true value function curve of reliability

表1 改變抽樣個數(shù)所得各個拉格朗日乘子Tab.1 Lagrange multipliers obtained by changing the number of samples

設置信水平P=90%，即α=0.10，將所得的17個c0值作為樣本進行自助再抽樣可得:置信區(qū)間的下邊界值c0L=-5.156 4；置信區(qū)間的上邊界值c0U=-4.869 1.

因此，參數(shù)c0的估計區(qū)間:［c0L，c0U］=［-5.156 4，-4.8691］.

同樣算得參數(shù)c1的估計區(qū)間:［c1L，c1U］=［0.4367，0.634 0］；參數(shù)c2的估計區(qū)間:［c2L，c2U］=［-1.004 8，-0.842 8］；參數(shù)c3的估計區(qū)間:［c3L，c3U］=［-0.015 7，0.009 9］；參數(shù)c4的估計區(qū)間:［c4L，c4U］=［-0.044 1，-0.028 4］；參數(shù)c5的估計區(qū)間:［c5L，c5U］=［-0.007 6，-0.002 3］.

改變每次抽樣個數(shù)，共抽取30 000次，可得映射參數(shù)如表2所示。

經(jīng)過計算，參數(shù)a、b的估計值分別為

分別將拉格朗日乘子的上限或下限值，映射參數(shù)的點估計值代入（24）式中，由排列組合原理和最大熵方法，可得64個概率密度函數(shù)。

表2 改變抽樣個數(shù)所得映射參數(shù)Tab.2 Mapping parameters obtained by changing the number of samples

經(jīng)過計算，由最小不確定性原理可得:當c0= c0L=-5.156 4，c1=c1U=0.634 0，c2=c2U=-0.842 8，c3=c3U=0.009 9，c4=c4U=-0.028 4，c5=c5U= -0.002 3時，可得到可靠性函數(shù)的上限曲線RU（x）；當c0=c0U=-4.869 1，c1=c1L=0.436 7，c2=c2L= -1.004 8，c3=c3L=-0.015 7，c4=c4L=-0.044 1，c5=c5L=-0.007 6時，可得到可靠性函數(shù)的下限曲線RL（x）.如圖3所示。

圖3 可靠性函數(shù)曲線Fig.3 Curves of reliability functions

由圖3可知可靠性函數(shù)取值隨著自變量時間的增大而逐漸減小，這是符合工程實際的。而且在x= 700 h之前，可靠性函數(shù)上下限曲線圖均與可靠性估計真值函數(shù)曲線圖基本完全重合。在x=700 h之后，隨著時間的增加，可靠性函數(shù)上下限曲線均越來越偏離可靠性估計真值函數(shù)曲線。這是因為在當前時間段內(nèi)，對可靠性函數(shù)估計的不確定度較小。隨著時間變量的增大，可靠性函數(shù)的估計難度逐漸增大，對可靠性函數(shù)估計的不確定度也逐漸增大，即上下限值越來越偏離估計真值。

設可靠性系數(shù)rc=0.1，置信水平P=90%，即顯著性水平α=0.1，計算出估計真值x0:x0=826 h.可靠性函數(shù)的區(qū)間估計為

案例2 根據(jù)文獻［18］，認為導彈壽命的概率分布為指數(shù)分布。對19個導彈進行壽命試驗，獲得的無失效數(shù)據(jù)如下:

用自助法每次抽取19個數(shù)據(jù)，共抽取30 000次，所得數(shù)據(jù)如圖4所示。

圖4 自助法獲取導彈樣本數(shù)據(jù)Fig.4 Sample data of missiles obtained by bootstrap method

基于最大熵法計算可得:映射參數(shù)a0=3.575 6，b0=-6.774 8.拉格朗日乘子（c00，c10，c20，c30，c40，c50）=（0.730 4，-0.150 8，-0.939 5，0.009 6，-0.103 9，-0.013 1）.

由（25）式可計算出概率密度估計真值函數(shù)f0（x），將f0（x）代入（26）式、（27）式、（28）式中，可得可靠性估計真值函數(shù)R0（x），如圖5所示。

文獻［18］中運用Bayes方法計算出R（3）= 0.916 4，運用多層自助最大熵方法計算出的R0（3）=0.898 0，二者相差0.018 4.這表明多層自助最大熵方法的可靠性估計真值結(jié)果與用Bayes方法的可靠性估計真值結(jié)果誤差很小?？梢娮灾畲箪胤椒ǐ@得的可靠性估計真值函數(shù)的擬合效果跟Bayes方法的估計結(jié)果基本一致。但是自助最大熵方法對概率分布沒有要求，而且能解決小樣本數(shù)據(jù)（n=19）的可靠性評估問題。因此，用多層自助最大熵方法評估機械產(chǎn)品的可靠性是有效且可行的。

改變每次抽樣個數(shù)t，共抽取30 000次，所得拉格朗日乘子如表3所示。

設置信水平P=90%，即α=0.10，將所得的16個c0值作為樣本進行自助再抽樣可得:置信區(qū)間的下邊界值c0L=0.359 4，置信區(qū)間的上邊界值c0U=0.592 4.

圖5 可靠性估計真值函數(shù)曲線Fig.5 Estimated true value function curve of reliability

表3 改變抽樣個數(shù)所得各個拉格朗日乘子Tab.3 Lagrange multipliers obtained by changing the number of samples

因此，參數(shù)c0的估計區(qū)間:［c0L，c0U］=［0.359 4，0.592 4］.

同樣算得參數(shù)c1的估計區(qū)間:［c1L，c1U］=［0.1934，0.363 0］；參數(shù)c2的估計區(qū)間:［c2L，c2U］=［-1.213 1，-0.998 8］；參數(shù)c3的估計區(qū)間:［c3L，c3U］=［-0.016 3，0.010 3］；參數(shù)c4的估計區(qū)間:［c4L，c4U］=［-0.043 2，0.011 1］；參數(shù)c5的估計區(qū)間:［c5L，c5U］=［-0.019 3，0.007 6］.

改變每次抽樣個數(shù)，共抽取30 000次，可得映射參數(shù)如表4所示。

經(jīng)過計算，參數(shù)a、b的估計值分別為

分別將拉格朗日乘子的上限或下限值，映射參數(shù)的點估計值代入（24）式中，由排列組合原理和最大熵方法，可得64個概率密度函數(shù)。

經(jīng)過計算，由最小不確定性原理可得:當c0= c0L=0.359 4，c1=c1U=0.363 1，c2=c2U=-0.998 8，c3=c3U=0.010 3，c4=c4U=0.011 1，c5=c5L= -0.019 3時，可得到可靠性函數(shù)的上限曲線RU（x）；

當c0=c0U=0.592 4，c1=c1L=0.193 4，c2=c2L= -1.2131，c3=c3L=-0.016 3，c4=c4L=-0.0432，c5=c5L=-0.019 3時，可得到可靠性函數(shù)的下限曲線RL（x）.如圖6所示。

表4 改變抽樣個數(shù)所得映射參數(shù)Tab.4 Mapping parameters obtained by changing the number of samples

圖6 可靠性函數(shù)曲線Fig.6 Curves of reliability functions

由圖6可知:可靠性函數(shù)取值隨著自變量時間的增大而逐漸減小，這是符合工程實際的。而且在x=1.9 a之前，可靠性函數(shù)上下限曲線圖均與可靠性估計真值函數(shù)曲線圖基本完全重合。在x=1.9 a之后，隨著時間的增加，可靠性函數(shù)上下限曲線均越來越偏離可靠性估計真值函數(shù)曲線。這是因為在當前時間段內(nèi)，對可靠性函數(shù)估計的不確定度較小。隨著時間變量的增大，可靠性函數(shù)的估計難度逐漸增大，對可靠性函數(shù)估計的不確定度也逐漸增大，即上下限越來越偏離估計真值。

設可靠性系數(shù)rc=0.1，置信水平P=90%，即顯著性水平為α=0.1，計算出估計真值x0:x0= 2.375 a.可靠性函數(shù)的區(qū)間估計為

案例3 某電子產(chǎn)品壽命的概率分布未知，對其進行模擬實驗，得到的無失效數(shù)據(jù)如下:

用自助法每次抽取10個數(shù)據(jù)，共抽取30 000次，所得數(shù)據(jù)如圖7所示。

圖7 自助法獲取樣本數(shù)據(jù)Fig.7 Sample data obtained by bootstrap method

基于最大熵法可得:映射參數(shù)a0=0.150 4，b0=-5.412 8.拉格朗日乘子（c00，c10，c20，c30，c40，c50）=（-2.473 8，-0.321 1，-0.994 8，-0.023 4，-0.031 2，0.004 7）.

由（25）式可計算出概率密度估計真值函數(shù)f0（x），將f0（x）代入（26）式、（27）式、（28）式中，可得可靠性估計真值函數(shù)R0（x），如圖8所示。

改變每次抽樣個數(shù)t，共抽取30 000次，可得拉格朗日乘子如表5所示。

設置信水平P=90%，即α=0.1，將所得的7個c0值作為樣本進行自助再抽樣可得:置信區(qū)間的下邊界值c0L=-2.772 4，置信區(qū)間的上邊界值c0U= -2.561 2.

因此，參數(shù)c0的估計區(qū)間:［c0L，c0U］=［-2.772 4，-2.561 2］.

同樣算得參數(shù)c1的估計區(qū)間:［c1L，c1U］=［-0.224 7，-0.095 6］；參數(shù)c2的估計區(qū)間:［c2L，c2U］=［-0.985 3，-0.785 4］；參數(shù)c3的估計區(qū)間:［c3L，c3U］=［0.002 0，0.032 5］；參數(shù)c4的估計區(qū)間:［c4L，c4U］=［-0.058 1，0.022 5］；參數(shù)c5的估計區(qū)間:［c5L，c5U］=［-0.002 6，0.003 1］.

表5 改變抽樣個數(shù)所得各個拉格朗日乘子Tab.5 Lagrange multipliers obtained by changing the number of samples

圖8 可靠性估計真值函數(shù)曲線Fig.8 Estimated true value function curve of reliability

改變每次抽樣個數(shù)，共抽取30 000次，可得映射參數(shù)如表6所示。

表6 改變抽樣個數(shù)所得映射參數(shù)Tab.6 Mapping parameters obtained by changing the number of samples

經(jīng)過計算，參數(shù)a、b的估計值分別為

分別將拉格朗日乘子的上限或下限值，映射參數(shù)的點估計值代入（24）式中，由排列組合原理和最大熵方法，可得64個概率密度函數(shù)。

經(jīng)過計算，由最小不確定性原理可得:當c0=c0L= -2.772 4，c1=c1U=-0.095 6，c2=c2U=-0.785 4，c3=c3U=0.032 5，c4=c4U=0.022 5，c5=c5L= -0.002 6時，可得到可靠性函數(shù)的上限曲線RU（x）；當c0=c0U=-2.561 2，c1=c1L=-0.224 7，c2=c2L=-0.985 3，c3=c3U=0.002 0，c4=c4L= -0.058 1，c5=c5U=0.003 1時，可得到可靠性函數(shù)的下限曲線RL（x）.如圖9所示。

圖9 可靠性函數(shù)曲線Fig.9 Curves of reliability functions

由圖9可知:可靠性函數(shù)取值隨著自變量時間的增大而逐漸減小，這是符合工程實際的。而且在x=33 h之前，可靠性函數(shù)上下限曲線均與可靠性估計真值函數(shù)曲線基本完全重合。在x=33 h之后，隨著時間的增加，可靠性函數(shù)上下限曲線均越來越偏離可靠性估計真值函數(shù)曲線。這是因為在當前時間段內(nèi)，對可靠性函數(shù)估計的不確定度較小。隨著時間變量的增大，可靠性函數(shù)的估計難度逐漸增大，對可靠性函數(shù)估計的不確定度也逐漸增大，即上下限值越來越偏離估計真值。

取可靠性系數(shù)rc=0.1，置信水平P=90%，即顯著性水平為α=0.1，計算出估計真值x0:x0= 47.5 h.

可靠性函數(shù)的區(qū)間估計:

案例4 根據(jù)文獻［19］，機電設備壽命T服從對數(shù)正態(tài)分布LN（μ，σ2）時。4臺電機分別工作到3 782.2 h、4 212.6 h、4 219.2 h和4 476.2 h均未失效。由于試驗數(shù)據(jù)收集的困難性，希望利用上述數(shù)據(jù)對該型號電機的可靠性進行評估。

用自助法每次抽取4個數(shù)據(jù)，共抽取10 000次，所得數(shù)據(jù)如圖10所示。

圖10 自助法獲取樣本數(shù)據(jù)Fig.10 Sample data of bearings obtained by bootstrap method

基于最大熵法計算可得:映射參數(shù)a0=0.003 0，b0=-12.476 7.拉格朗日乘子（c00，c10，c20，c30，c40，c50）=（-7.272 5，1.176 6，0.038 2，-0.637 4，-0.053 8，0.071 0）.

由（25）式可計算出概率密度估計真值函數(shù)f0（x），將f0（x）代入（26）式、（27）式、（28）式中，可得可靠性估計真值函數(shù)R0（x），如圖11所示。

圖11 可靠性估計真值函數(shù)曲線Fig.11 Estimated true value function curve of reliability

運用多層自助最大熵方法計算出的R0（4 000）= 0.955 8，R0（4 500）=0.925 0.文獻［19］中運用Bayes方法計算出當λ=0.2時，μ的估計值為10.721 8，σ的估計值為1.767 3，從而可得R（4 000）=0.915 2；R（4 500）=0.904 4.可見兩種方法計算結(jié)果分別相差0.040 6和0.020 6，多層自助最大熵方法的可靠性估計真值結(jié)果與用Bayes方法的可靠性估計真值結(jié)果誤差很小。因此自助最大熵方法獲得的可靠性估計真值函數(shù)的擬合效果跟Bayes方法的估計結(jié)果基本一致。自助最大熵方法可以自動識別樣本數(shù)據(jù)內(nèi)部規(guī)律，從而計算出小樣本數(shù)據(jù)（n=4）的概率分布函數(shù)。在數(shù)據(jù)處理過程中，并未利用已知的概率分布先驗信息。因此，用多層自助最大熵方法評估機械產(chǎn)品的可靠性是有效且可行的。

值得注意的是，由于樣本數(shù)據(jù)較少，無法通過改變抽樣個數(shù)獲得拉格朗日乘子的樣本數(shù)據(jù)，因此無法進行可靠性區(qū)間估計。進行可靠性區(qū)間估計需要原樣本數(shù)據(jù)個數(shù)n≥7.

案例5 針對案例4，進行模擬試驗，得到的無失效數(shù)據(jù)如下:

用自助法每次抽取7個數(shù)據(jù)，共抽取10 000次，所得數(shù)據(jù)如圖12所示。

圖12 自助法獲取樣本數(shù)據(jù)Fig.12 Sample data obtained by bootstrap method

基于最大熵法可得:映射參數(shù)a0=0.009 2，b0=-38.354 0.拉格朗日乘子（c00，c10，c20，c30，c40，c50）=（-5.653 8，1.294 6，-1.345 5，0.412 5，0.045 0，0.092 4）.

由（25）式可計算出概率密度估計真值函數(shù)f0（x），將f0（x）代入（26）式、（27）式、（28）式中，可得可靠性估計真值函數(shù)R0（x），如圖13所示。

改變每次抽樣個數(shù)t，共抽取30 000次，可得拉格朗日乘子如表7所示。

設置信水平P=90%，即α=0.1，將所得的7個c0值作為樣本進行自助再抽樣可得:置信區(qū)間的下邊界值c0L=-5.804 9，置信區(qū)間的上邊界值c0U= -5.703 2.

因此，參數(shù)c0的估計區(qū)間:［c0L，c0U］=［-5.804 9，-5.703 2］.

同樣算得參數(shù)c1的估計區(qū)間:［c1L，c1U］=［1.248 6，1.624 0］；參數(shù)c2的估計區(qū)間:［c2L，c2U］=［-1.582 4，-0.704 4］；參數(shù)c3的估計區(qū)間:［c3L，c3U］=［0.002 0，0.032 5］；參數(shù)c4的估計區(qū)間:［c4L，c4U］=［-0.098 7，0.082 6］；參數(shù)c5的估計區(qū)間:［c5L，c5U］=［-0.065 2，0.091 5］.

改變每次抽樣個數(shù)，共抽取30 000次，可得映射參數(shù)如表8所示。

圖13 可靠性估計真值函數(shù)曲線Fig.13 Estimated true value function curve of reliability

表7 改變抽樣個數(shù)所得各個拉格朗日乘子Tab.7 Lagrange multipliers obtained by changing the number of samples

表8 改變抽樣個數(shù)所得映射參數(shù)Tab.8 Mapping parameters obtained by changing thenumber of samples

經(jīng)過計算，參數(shù)a、b的估計值分別為

分別將拉格朗日乘子的上限或下限值，映射參數(shù)的點估計值代入（24）式中，由排列組合原理和最大熵方法，可得64個概率密度函數(shù)。

經(jīng)過計算，由最小不確定性原理可得:當c0= c0L=-5.804 9，c1=c1U=1.624 0，c2=c2U=-0.704 4，c3=c3U=0.032 5，c4=c4U=0.082 6，c5=c5U=0.091 5時，可得到可靠性函數(shù)的上限曲線RU（x）；當c0= c0U=-5.703 2，c1=c1L=1.248 6，c2=c2U=-0.704 4，c3=c3L=0.002 0，c4=c4L=-0.098 7，c5=c5L= -0.065 2時，可得到可靠性函數(shù)的下限曲線RL（x）.如圖14所示。

圖14 可靠性函數(shù)曲線Fig.14 Curves of reliability functions

由圖14可知:可靠性函數(shù)取值隨著自變量時間的增大而逐漸減小，這是符合工程實際的。而且在x=4 230 h之前，可靠性函數(shù)上下限曲線均與可靠性估計真值函數(shù)曲線基本完全重合。在x=4 230 h之后，隨著時間的增加，可靠性函數(shù)上下限曲線均越來越偏離可靠性估計真值函數(shù)曲線。這是因為在當前時間段內(nèi)，對可靠性函數(shù)估計的不確定度較小。隨著時間變量的增大，可靠性函數(shù)的估計難度逐漸增大，對可靠性函數(shù)估計的不確定度也逐漸增大，即上下限值越來越偏離估計真值。

取可靠性系數(shù)rc=0.1，置信水平P=90%，即顯著性水平為α=0.1，計算出估計真值x0:x0= 4 425.4 h.可靠性函數(shù)的區(qū)間估計:

3 討論

1）對于概率分布已知的產(chǎn)品壽命可靠性評估問題，自助最大熵方法的評估誤差很小。案例1用多層自助最大熵法計算出的R0（577）=0.995 2與文獻［17］中運用Bayes方法計算出的R（577）= 0.984 5相差0.010 7；R0（677）=0.994 5與文獻［17］中運用Bayes方法計算出的R（677）=0.976 2相差0.018 3.文獻［18］中運用Bayes方法計算出R（3）=0.916 4，案例2運用多層自助最大熵法計算出的R0（3）=0.898 03，二者相差0.018 37.案例4運用多層自助最大熵方法計算出的R0（4 000）= 0.955 8，R0（4 500）=0.925 0.文獻［19］中運用Bayes方法計算出R（4 000）=0.915 2，R（4 500）= 0.904 4，兩種方法計算結(jié)果分別相差0.040 6和0.020 6.試驗案例1、案例2、案例4表明運用自助最大熵方法得到的可靠性評估結(jié)果與運用Bayes方法得到的評估結(jié)果基本相同。這是因為最大熵方法能夠?qū)ξ粗母怕史植甲龀鲋饔^偏見為最小的最佳估計，可以自動識別樣本數(shù)據(jù)內(nèi)部規(guī)律，從而計算出樣本數(shù)據(jù)個數(shù)極少（n≥4）條件下的概率分布函數(shù)。多層自助最大熵法在數(shù)據(jù)處理過程中，并未利用已知概率分布（Weibull分布、指數(shù)分布、對數(shù)正態(tài)分布）這一信息。

2）對于概率分布未知的產(chǎn)品壽命可靠性評估問題，古典統(tǒng)計理論無法解決。而案例3運用多層自助最大熵方法，計算得到R0（56.97）=0.898 9；案例5運用多層自助最大熵方法計算得到R0（4 500）= 0.897 5.試驗案例表明多層自助最大熵方法不但適用于概率分布已知的產(chǎn)品壽命可靠性評估問題，也適用于概率分布未知的產(chǎn)品壽命可靠性評估問題。

3）對于可靠性函數(shù)的區(qū)間估計問題，現(xiàn)有方法無法解決。本文中的多層自助最大熵方法把拉格朗日乘子當作變量，對可靠性函數(shù)進行區(qū)間估計，得出結(jié)論:可靠性函數(shù)取值隨著自變量時間的增大而逐漸減小，在一定范圍內(nèi)，可靠性函數(shù)上下限曲線圖與可靠性估計真值函數(shù)曲線圖基本完全重合。超出該范圍，可靠性函數(shù)上下限曲線均越來越偏離可靠性估計真值函數(shù)曲線。

4）實際應用中，一般在達到產(chǎn)品無失效數(shù)據(jù)樣本最大值的一半之前，可靠性函數(shù)上下限曲線圖與可靠性估計真值函數(shù)曲線圖基本完全重合。如案例1中x=700 h＞1/2×924=462 h，案例2中x= 1.9 a＞1/2×3=1.5 a，案例3中x=33 h＞1/2× 56.97=28.485 h.

4 結(jié)論

將拉格朗日乘子當做變量，提出多層自助最大熵方法，以解決機械產(chǎn)品的乏信息可靠性評估問題。該方法對概率分布沒有要求，可以有效地解決乏信息條件下無失效數(shù)據(jù)的可靠性評估問題與可靠性函數(shù)的區(qū)間估計問題，試驗證明該方法是對現(xiàn)有可靠性評估方法的有益補充。

（References）

［1］ 傅惠民，張勇波.Weibull分布定時無失效數(shù)據(jù)可靠性分析方法［J］.航空動力學報，2010，25（12）:2807-2810. FU Hui-min，ZHANG Yong-bo.Method of reliability analysis for time truncated zero-failure data based on Weibull distribution［J］. Journal of Aerospace Power，2010，25（12）:2807-2810.（in Chinese）

［2］ 王憑慧，范本堯，傅惠民.衛(wèi)星推力器可靠性評估和壽命預測［J］.航空動力學報，2004，19（6）:745-748. WANG Ping-hui，F(xiàn)AN Ben-yao，F(xiàn)U Hui-min.Reliability assessment and life prediction for satellite engine［J］.Journal of Aerospace Power，2004，19（6）:745-748.（in Chinese）

［3］ 李寶盛，何洪慶.“小子樣、零失效”情況下壽命可靠度的置信分析方法［J］.兵工學報，2001，22（2）:234-237. LI Bao-sheng，HE Hong-qing.Method of confidence analysis for life reliability on small and zero failure samples［J］.Acta Armamentarii，2001，22（2）:234-237.（in Chinese）

［4］ 蒲星.一種無失效數(shù)據(jù)可靠性分析的方法［J］.電子質(zhì)量，2014（7）:29-31. PU Xing.A method for reliability analysis of zero-failure data［J］. Journal of Electronic Quality，2014（7）:29-31.（in Chinese）

［5］ 但召江，樓洪梁，李興林，等.無失效數(shù)據(jù)下滾動軸承的可靠性估計［J］.軸承，2013（9）:22-24. DAN Shao-jiang，LOU Hong-liang，LI Xing-lin，et al.Reliability estimation for rolling bearings based on zero-failure data［J］.Journal of Bearing，2013（9）:22-24.（in Chinese）

［6］ Jiang P，Lim J H，Zuo M J，et al.Reliability estimation in a Weibull lifetime distribution with zero-failure field data［J］. Quality and Reliability Engineering International，2010，26（7）: 691-701.

［7］ Nam J S，Kim H E，Kim K U.A new accelerated zero-failure test model for rolling bearings under elevated temperature conditions［J］.Journal of Mechanical Science and Technology，2013，27（6）:1801-1807.

［8］ Bailey R T.Estimation from zero-failure data［J］.Risk Analysis，1997，17（3）:375-380.

［9］ Kwon Y.Design of Bayesian zero-failure reliability demonstration test for products with Weibull lifetime distribution［J］.Journal of Applied Reliability，2014，14（4）:220-224.

［10］ Kwon Y.The effect of scale parameter in designing reliability demonstration test for lognormal lifetime distribution［J］.Journal of Applied Reliability，2014，14（1）:53-57.

［11］ 夏新濤，陳曉陽，張永振，等.兩個數(shù)據(jù)序列的灰關系及其工程應用的兩個例子［J］.兵工學報，2008，29（1）:73-77. XIA Xin-tao，CHEN Xiao-yang，ZHANG Yong-zhen，et al.E-valuation for stability of manufacturing process based on grey relation［J］.Acta Armamentarii，2008，29（1）:73-77.（in Chinese）

［12］ 賈波，劉福，雷正偉.基于最大熵方法對測量數(shù)據(jù)估計的改進方法研究［J］.價值工程，2010（28）:228-230. JIA Bo，LIU Fu，LEI Zheng-wei.Study on evaluation of measurement result based on improved maximum entropy method［J］. Journal of Value Engineering，2010（28）:228-230.（in Chinese）

［13］ 夏新濤，孟艷艷，邱明.用灰自助泊松方法預測滾動軸承振動性能可靠性的變異過程［J］.機械工程學報，2015，51（9）: 97-103. XIA Xin-tao，MENG Yan-yan，QIU Ming.Forecasting for variation process of reliability of rolling bearing vibration performance using grey bootstrap Poisson method［J］.Journal of Mechanical Engineering，2015，51（9）:97-103.（in Chinese）

［14］ 胡松偉.最大熵法在發(fā)射裝置可靠性評估中的應用［J］.四川兵工學報，2014，35（5）:31-33. HU Song-wei.Maximum entropy method application on missile launching system reliability evaluation［J］.Journal of Sichuan Ordnance，2014，35（5）:31-33.（in Chinese）

［15］ Yeom D H，Bae S G.A constraint-based maximum entropy sampling method for Kriging models in fuel cell applications［J］. Quality and Reliability Engineering International，2013，29（6）: 819-828.

［16］ Young H S，Byung M K.Reliability bound based on the maximum entropy principle with respect to the first truncated moment［J］.Journal of Mechanical Science and Technology，2010，24（9）:1891-1900.

［17］ 茆詩松，夏劍鋒，管文琪.軸承壽命試驗中無失效數(shù)據(jù)的處理［J］.應用概率統(tǒng)計，1993，9（3）:326-331. MOU Shi-song，XIA Jian-feng，GUAN Wen-qi.Processing of zero-failure data in bearing life test［J］.Chinese Journal of Applied Probability and Statistics，1993，9（3）:326-331.（in Chinese）

［18］ 孫亮，徐廷學，王冬梅.某型導彈無失效數(shù)據(jù)的處理方法［J］.戰(zhàn)術導彈技術，2004（3）:29-32. SUN Liang，XU Ting-xue，WANG Dong-mei.Process method for zero-failure data of certain type missile［J］.Tactical Missile Technology，2004（3）:29-32.（in Chinese）

［19］ 孫建中.機電設備無失效數(shù)據(jù)的可靠性評估［J］.艦船科學技術，2006，28（1）:50-53. SUN Jian-zhong.Reliability estimation of zero-failure data for mechanical and electrical products［J］.Ship Science and Technology，2006，28（1）:50-53.（in Chinese）

Reliability Evaluation Based on Hierarchical Bootstrap Maximum Entropy Method

XIA Xin-tao，YE Liang，LI Yun-fei，CHANG Zhen
（School of Mechatronical Engineering，Henan University of Science and Technology，Luoyang 471003，Henan，China）

A hierarchical bootstrap maximum entropy evaluation model is proposed to analyze the life reliability of mechanical products under the condition of small samples without any prior information.Adequate sample data is obtained by using bootstrap method to re-sample the current zero-failure data samples.Based on maximum entropy method，the different Lagrange multipliers can be obtained by changing the number of samples.In order to get the interval estimation values of Lagrange multipliers，the bootstrap method is used again to re-sample the small sample data of Lagrange multipliers.The probability density functions and reliability functions are achieved by carrying on permutation and combination for the upper and lower limit values of each Lagrange multiplier，so the interval estimation values of reliability functions can be gained using minimum uncertainty principle.Experimental investigation shows that the hierarchical bootstrap maximum entropy evaluation model can effectively solve the reliability evaluation problem for zero-failure data of small samples with known or unknown probability distributions.

system assessment and feasibility analysis；reliability evaluation；hierarchical bootstrap maximum entropy method；poor information；zero-failure data；Lagrange multipliers

TB114.3

A

1000-1093（2016）07-1317-13

10.3969/j.issn.1000-1093.2016.07.022

2016-01-05

國家自然科學基金項目（51475144、51075123）

夏新濤（1957—），男，教授，博士生導師，博士。E-mail:xiaxt1957@163.com；

葉亮（1990—），男，碩士研究生。E-mail:172682823@qq.com