随着机械产品或系统的可靠性要求不断提升，传统的基于大样本失效数据的统计研究方法相较于小样本特征出现了更多的不适用性。在航空航天、军工、核电等领域，小样本情形下系统失效存在很大的随机性，失效数据会越来越趋向于小样本化，这给可靠性研究带来更多的困难。可靠性研究的关键问题在于确定反映产品或系统失效特征的失效分布，在实际研究中，通常将失效数据拟合成具体的失效分布模型，并估计其相应的分布参数，然后进行可靠性评估。威布尔分布作为一种研究较为深入的失效分布模型，在机械可靠性工程领域有着广泛的应用。威布尔分布的概率密度函数为

f ( t ) = β η ( t − γ η ) β − 1 exp [ − ( t − γ η ) β ] , t ≥ γ 。

(1)

式中， β 为形状参数，其决定了分布曲线的形状，在可靠性评估中，产品的可靠寿命和特征寿命随着形状参数的增大而提高^[1]；η 为尺度参数，能够放大或缩小坐标尺度，但不能影响威布尔分布的形状； γ为位置参数，即最小寿命，其决定了分布曲线的起始位置。以上为三参数威布尔分布模型。当 γ=0时，式(1)就简化为两参数威布尔分布。目前关于威布尔分布的研究主要有3个方向：威布尔分布的参数估计研究^[2-3]、运用威布尔分布对无失效数据分析研究^[4-5]和小样本特征下威布尔分布的可靠性研究^[6-7]。在机械可靠性评估研究中，两参数威布尔分布因其不能充分描述其失效特征，小样本失效数据下对参数估计也存在较大误差，因而本文采用三参数威布尔分布模型进行机械可靠性评估研究。相较于两参数威布尔分布，三参数威布尔分布在机械可靠性评估中有如下优点：首先，机械产品或系统失效都表现为具有损耗特征的渐变性失效，磨损、疲劳和老化都是慢慢积累而成以致最后失效，失效通常都需要一定的时间过程，例如滚动轴承在工作中产生的点蚀会逐渐渐变为凹坑，当凹坑面积达到相应的限定值时便会引起失效，因此威布尔分布中位置参数的存在是有其实际意义的；其次，在于机械可靠性评估中，研究中除了关注其可靠寿命、平均寿命等特征量外，还可能关注其最小寿命，而位置参数在某种意义上就代表了最小寿命，这与机械设备的保养等直接相关；最后三参数威布尔分布模型适应性更强，相较于两参数威布尔分布具有较小的误差，对于机械产品或系统的失效模型拟合程度更好。

对于小样本失效数据下三参数威布尔分布的可靠性评估问题，目前研究要么集中于小样本问题的解决方法^[8]，要么集中于对于三参数威布尔分布的参数估计研究^[9]，而对于小样本失效数据下的三参数威布尔分布的机械可靠性评估问题，目前鲜有集中全面的研究。本文以贝叶斯方法为解决小样本问题的主要方法。现有研究中关于贝叶斯方法在可靠性研究中的运用路径一般都是先求出关键的先验分布，然后结合小样本数据计算其后验分布，最后得到可靠性参数估计值。但将贝叶斯方法运用于威布尔分布特别是3个参数的复杂情形时，计算其解析解存在较大的困难^{[6, 10]}。灰色估计法^[11]因其可以一次性估计出威布尔分布的3个参数，且在小样本情形下相对于其他方法也有较高的精度，所以本文采用灰色估计法对小样本失效数据进行初步拟合。但又由于小样本数据的数据量较少，直接将灰色估计后的初步结果应用于机械可靠性评估中时，必然使得计算结果准确度较低；而将威布尔分布直接应用于贝叶斯分析又缺少相应的共轭先验分布^[12]。因此，为确定小样本失效数据的三参数威布尔先验分布，本文采用蒙特卡洛方法对灰色估计得到的威布尔分布函数进行随机抽样，以获得多组再生样本作为Bootstrap子样，运用贝叶斯公式计算抽样结果的离散概率分布，最后通过计算全概率得到经贝叶斯方法修正后的可靠性指标。该方法极大地简化了三参数威布尔分布直接运用贝叶斯公式时存在解析解难以计算的问题，为小样本下三参数威布尔分布的可靠性评估提供了新的思路。

1 贝叶斯-蒙特卡洛方法及先验分布的确定

贝叶斯-蒙特卡洛方法(Bayesian Monte-Carlo方法，简称BMC方法)^[13]是一种通过将蒙特卡洛仿真结果纳入先验信息来执行贝叶斯公式的有效方法。BMC方法在保留了蒙特卡洛方法特征的基础上提升了计算效率，减少了方法的执行费用，对于小样本问题的机械可靠性评估具有一定针对性。本节在描述贝叶斯-蒙特卡洛方法的基本步骤的基础上，用灰色估计法和蒙特卡洛抽样确定三参数威布尔分布参数的贝叶斯先验分布。相较于文献[13]介绍的BMC方法框架，本文将灰色估计法与蒙特卡洛抽样得到Bootstrap子样相结合来确定先验分布，在参数估计更准确的同时，也得到了更为可靠的先验信息。

1.1 BMC方法简介和贝叶斯公式

BMC方法在三参数威布尔分布的机械可靠性评估中的运用过程为：1)选取威布尔分布为机械可靠性研究对象，通过两参数与三参数对比确定三参数威布尔分布为研究对象；2)利用灰色估计法对小样本数据进行参数估计并检验，对估计结果进行蒙特卡洛抽样，将得到Bootstrap子样同样利用灰色估计法进行三参数的估计，以此作为先验信息来确定威布尔分布参数的贝叶斯先验分布；3)结合先验分布和小样本数据，利用贝叶斯公式计算后验分布；4)利用威布尔分布参数的后验分布计算可靠度、失效率、平均无故障工作时间等可靠性指标，以此对产品或系统进行可靠性评估。

贝叶斯理论的核心在于贝叶斯公式，贝叶斯公式综合了先验信息、总体信息和样本信息，其表达式为

π ( θ | x ) = π ( θ ) p ( x | θ ) m ( x ) 。

(2)

其中， θ 为总体分布中的未知参数，即所要估计的威布尔分布中的形状参数 β 、尺度参数 η 和位置参数 γ ； π ( θ ) 为参数的先验分布； x = ( x 1 , x 2 , ⋯ , x n ) 为来自总体 p ( x | θ ) 的一个随机样本。

m ( x ) = ∫ Θ π ( θ ) p ( x | θ ) d θ 。

(3)

m ( x ) 为与 θ 无关的 x 的边际概率密度函数，其不包括 θ 的任何信息。 π ( θ | x ) 为估计参数的后验概率分布。

p ( x | θ ) = ∏ i = 1 n p ( x i | θ ) 。

(4)

1.2 先验分布的计算

执行贝叶斯方法的关键在于先验分布的选取，不同的先验分布确定方法最后的可靠性评估结果可能存在较大差异。目前较为普遍的先验分布确定方法包括贝叶斯假设法、杰弗莱原则^[10]以及共轭分布法等。三参数威布尔分布作为一种形式上相对复杂的失效分布，并不存在相应的共轭分布，这使得其先验分布的确定较为困难。本文结合小样本失效数据特征和威布尔分布缺少共轭先验分布的实际，采用灰色估计法初步拟合三参数威布尔分布并进行分布拟合检验，对拟合后的分布函数采用蒙特卡洛抽样得到的Bootstrap子样来作为威布尔分布参数的贝叶斯先验分布。该先验分布确定方法克服了小样本数据量不足的缺陷，同时也避免了因威布尔分布参数多而导致的计算复杂的问题。

由式(1)可知三参数威布尔分布的分布函数为

F ( x | β , η , γ ) = 1 − e x p [ − ( x − γ η ) β ] , x ≥ γ 。

(5)

式(5)对应的可靠度函数为

R ( x | β , η , γ ) = exp [ − ( x − γ η ) β ] , x ≥ γ 。

(6)

由文献[14]可得 ( x i , t i ) 的一般时间序列为

x i = η exp ( t i β ) + γ 。

(7)

式中， t i = ln [ ln ( 1 / R ( x i ) ) ] 。灰色模型GM(1,1)^[15]的基本形式为

x = c exp ( a t i ) + b 。

(8)

具体建模步骤^[16]如下。1) 将试验得到的失效数据按从小到大排序得 x 1 ≤ x 2 ≤ ⋯ ≤ x n ；2)利用中位秩法计算可靠度 R ( x i ) = 1 − ( i − 0.3 ) / ( n + 0.4 ) ；3)令 t i = ln ln ( 1 / R ( x i ) ) ，将 ( x i , t i ) 视为一般时间序列，可得到一阶灰色模型的时间响应模型为 x ^ ( t ) = c e − a t + b ，最后比对三参数威布尔分布可得 β = − 1 / a , η = c , γ = b 。灰色估计方法经过一次建模过程即可得到威布尔分布的3个参数，且在小样本情形下仍具有较高精度，是一种较好的三参数威布尔分布的参数估计方法。为使小样本数据拟合得到的三参数威布尔分布在一定程度上具有稳健性，需要对其进行拟合优度检验。常用的拟合优度检验方法包括 χ 2 检验法、正态分布检验法和K-S检验法。K-S检验法能较好地适应小样本下概率分布的检验，因而本文采用K-S检验法对三参数威布尔分布进行拟合优度检验^[17]。

为计算分布参数的贝叶斯先验分布，本文以小样本失效数据作为全样本，通过灰色估计法得到威布尔分布的参数估计值 θ ^ ( x ) ，将估计值 θ ^ ( x ) 代入三参数威布尔分布函数中，利用蒙特卡洛抽样得到 m 组Bootstrap子样，然后对这 m 组子样继续使用灰色估计法进行威布尔分布拟合，得到的参数估计值记为 θ ^ 1 ( x * ) , θ ^ 2 ( x * ) , ⋯ , θ ^ m ( x * ) ，将 θ ^ 1 ( x * ) , θ ^ 2 ( x * ) , ⋯ , θ ^ m ( x * ) 作为参数的贝叶斯先验分布。该方法解决了三参数威布尔分布贝叶斯分析中不存在共轭分布的问题，也充分利用了先验信息。蒙特卡洛方法抽样得到的Bootstrap子样也降低了对小样本失效数据的依赖性，相较于直接的Bootstrap法提高了参数估计的稳健性。

假设小样本失效数据 t 1 , t 2 , ⋯ , t r , ⋯ , t n 为全样本无截尾数据，即 n = r ， t 1 ≤ t 2 ≤ ⋯ t n 。一般认为样本个数 n ＜ 30 为小样本数据^[18]。将前述由灰色估计法计算出的威布尔分布参数值 β 、 η 、 γ 代入抽样公式进行抽样，产生 m 组每组样本容量为 n 的Bootstrap子样。抽样先产生在[0,1]区间上均匀分布的伪随机数数组 u [ n ] ，再由反函数抽样法抽得服从三参数威布尔分布的随机变量，该过程可借助Matlab软件来实现。抽样公式为

t [ n ] = η ( − ln u [ n ] ) 1 / β + γ 。

(9)

将抽到的 m 组子样数据 t [ n ] 按从小到大排列，对每组数据分别利用灰色估计法进行三参数威布尔分布的参数估计，假设得到的各组参数值记为 β i 、 η i 、 γ i ，将其估计值分别保存于数组 β [ i ] 、 η [ i ] 、 γ [ i ] 中。在 β i 、 η i 、 γ i 相互独立的情况下，可得 β 、 η 、 γ 的离散先验联合分布列为 p i j k = P ( β = β i , η = η j , γ = γ k ) ，由概率论相关理论可知 β i 、 η i 、 γ i 的联合概率密度为常数^[17]。

π ( β = β i , η = η j , γ = γ k ) = 1 / m 3 , ​ ​ i , j , k = 1 , 2 , ⋯ , m 。

(10)

式(10)即为所要求解的三参数威布尔分布的贝叶斯先验分布。

2 后验分布及可靠性指标的计算

由式(1)可得三参数威布尔分布的失效概率密度函数的参数条件密度函数为

f ( t i | β = β i , η = η j , γ = γ k ) = β i η j ( t [ i ] − γ k η j ) β i − 1 exp [ − ( t [ i ] − γ k η j ) β i ] 。

(11)

似然函数为

f ( x | β = β i , η = η j , γ = γ k ) = ∏ i = 1 n f ( x i | β = β i , η = η j , γ = γ k ) = β i n η j − n β i U β i − 1 exp ( − V / η j β i ) 。

(12)

式中， U = ∏ i = 1 n ( t [ i ] − γ k ) ， V = ∑ i = 1 n ( t [ i ] − γ k ) β i 。

三参数威布尔分布的参数先验分布为

π ( β = β i , η = η j , γ = γ k ) = 1 / m 3 。

(13)

由贝叶斯公式可得参数 β 、 η 、 γ 的后验分布为

π ( β = β i , η = η j , γ = γ k | x ) = f ( x | β ​ = β i , η ​ = ​ η j , γ = ​ γ k ) π ( β = β i , η = η j , γ = γ k ) ∫ f ( x | β ​ = ​ β i , η = η j , γ = γ k ) π ( β = β i , η = η j , γ = γ k ) d ​ x = β i n η j − n β i U β i − 1 exp ( − V / η j β i ) ∑ i = 1 m ∑ j = 1 m ∑ k = 1 m β i n η j − n β i U β i − 1 exp ( − V / η j β i )

(14)

式中， i , j , k = 1 , 2 , ⋯ , m ，计算结果即为贝叶斯公式修正后的三参数威布尔分布参数的后验概率分布。由式(17)和全概率公式即可计算可靠度函数 R ( t ) 、失效概率密度函数 f ( t ) 、失效率函数 λ ( t ) 和平均无故障工作时间(MTBF)的点估计值。

可靠度函数的点估计为

R ( t ) = ∑ i = 1 m ∑ j = 1 m ∑ k = 1 m R i j k ( t ) π [ i ] [ j ] [ k ] = ∑ i = 1 m ∑ j = 1 m ∑ k = 1 m π [ i ] [ j ] [ k ] exp [ − ( t − γ k η j ) β i ] 。

(15)

失效概率密度函数的点估计为

f ( t ) = ∑ i = 1 m ∑ j = 1 m ∑ k = 1 m f i j k ( t ) π [ i ] [ j ] [ k ] = ∑ i = 1 m ∑ j = 1 m ​ ∑ k = 1 m ​ π [ i ] [ j ] [ k ] β i η j ( ​ t ​ − γ k η j ) β i ​ − 1 exp [ − ( t ​ − ​ γ k η j ) β i ] 。

(16)

失效率函数的点估计为

λ ( t ) = f ( t ) / R ( t ) 。

(17)

平均无故障工作时间为

M T B F = ∫ 0 ∞ R ( t ) d t 。

(18)

3 数值算例

文献[19]给出了某计算机数控系统的失效数据，本文以该失效数据为样本进行BMC方法的计算。表1给出了该计算机数控系统半年使用时间的故障数据。

3.1 与两参数威布尔分布的比较

为了体现三参数威布尔分布在机械产品失效数据拟合效果上优于两参数威布尔分布，对表1数据分别采用两参数威布尔分布模型和三参数威布尔分布模型进行拟合，两参数威布尔分布模型采用极大似然法进行估计，三参数威布尔分布模型则采用前述的灰色估计法进行估计，结果如表2所示。

结合表2的参数估计结果和文献[20]的结论可知，三参数估计的形状参数 β 总比两参数的估计值要小一些，尺度参数 γ 表明了三参数威布尔分布中的机械最小寿命的存在，使其在小样本失效数据下更能体现实际的失效特征。

3.2 三参数威布尔分布机械可靠性评估

将得到参数初步估计值代入三参数威布尔分布模型中并进行K-S分布拟合检验。K-S检验值0.183 35< D n , α =0.296 8,说明原假设服从该参数的威布尔分布模型成立，所以三参数威布尔分布的初步估计模型为

F ( t ) = 1 − exp [ − ( t − 54.46 595.13 ) 1.37 ] 。

(19)

从而根据式(9)可得蒙特卡洛抽样公式为

t [ n ] = 595.13 ( − ln u [ n ] ) 1 / 1.37 + 54.46 。

(20)

由式(20)抽得 m 组每组容量为21的Bootstrap子样，抽样组数 m 通常越大越能使得结果稳定，但抽样费用也相应增大。本文参考文献[18]的做法，认为 m 取决于可靠度曲线的斜率，这里取 m =30，即用蒙特卡洛方法抽样产生30组每组样本容量为21的Bootstrap子样，对30组子样继续使用灰色估计方法进行参数估计，将得到的各子样组的 β 、 η 、 γ 参数估计值分别保存于数组 β [ i ] 、 η [ i ] 、 γ [ i ] 中，具体计算结果如表3所示。

利用表3的数据结合式(14)可计算出 β 、 η 、 γ 的后验概率分布，然后根据公式(15)~(18)可分别计算出典型运行时间点的可靠度、失效概率密度、失效率和平均无故障工作时间等可靠性指标值，具体如表4所示。

该数控系统的可靠度、失效概率密度以及失效率随时间变化如图1~图3所示。

从表4和图1~图3可以看出，使用本文方法进行可靠性评估的结果与实际期望结果较为一致。系统的可靠度随时间的增加而逐渐降低，与实际情况较为吻合，失效概率密度和失效率在投入使用后较短的一段时间内并不为零，也反映了该系统可能存在早期投入使用时设备不稳定产生的早期失效情形。

为展现经贝叶斯公式修正后的三参数威布尔分布参数的后验分布相较于单纯的灰色估计法更加精确，选取MTBF用作比较的可靠性指标。由原始小样本失效数据经灰色估计后得到三参数估计值 β = 1.37 , η = 595.13 , γ = 54.46 ，其MTBF大约为500 h，单独的灰色估计得到参数后计算的MTBF为^[20]

M T B F ​ = ​ 54.46 + 595.13 Γ ( 1 + 1 / 1.37 ) ​ = ​ 597.58 h 。

(21)

由式(18)计算得到MTBF=465.34 h，可见经贝叶斯定理修正后的可靠性指标值更接近于实际值，相较于单纯的灰色估计法进行参数估计后的计算更为准确。蒙特卡洛抽样组数 m 、初始小样本失效数据数 n 等都会对可靠性评估结果产生影响，抽样产生的Bootstrap子样作为贝叶斯先验分布，经修正后可在一定程度上用于可靠性评估。

4 结语

本文结合当前复杂产品或系统“高可靠，少失效”的小样本失效数据的实际，运用贝叶斯-蒙特卡洛方法对小样本下三参数威布尔分布的可靠性评估进行了分析。相较于大样本失效数据可靠性评估方法，该方法综合运用了蒙特卡洛抽样产生自助子样来扩大小样本和贝叶斯方法结合子样信息作为先验信息，不需要很多失效数据即可完成对小样本下复杂系统可靠性评估，对实际的可靠性分析工作有一定的参考价值。而三参数威布尔分布相较于两参数虽然更加复杂但更能反映实际的失效特征和描绘失效分布，用BMC方法对其进行分析评估能产生较好的效果。本文在对三参数威布尔分布进行初步拟合时利用了灰色估计方法，因其能一次性估计出威布尔分布的3个参数，本文认为该方法相对于其他参数估计方法有一定的优势，但估计精度等其他问题还需要进一步的研究。