CN109145502B - 一种威布尔型单元寿命分布参数估计方法 - Google Patents

Abstract

本发明涉及一种威布尔型单元寿命分布参数估计方法，该方法首先根据威布尔型单元的使用寿命数据生成n组候选的分布参数，然后初始化似然度，再根据k次检查结果获取的完好单元的数量、故障单元的数量以及检查时刻，依次更新似然度；最后在更新后的似然度中找到最大似然度，则最大似然度对应的威布尔分布尺度参数和形状参数即为估计结果。本方法的参数估计结果在总体上能“跟随”理论成熟方法的参数估计结果，其估计精度能满足工程要求。

Description

一种威布尔型单元寿命分布参数估计方法

技术领域

本发明涉及产品质量检测技术领域，具体涉及一种威布尔型单元寿命分布参数估计方法。

背景技术

产品可靠性是一种描述产品质量的核心属性，常用产品寿命的分布类型和参数来表达产品的可靠性。准确获知产品的可靠性，是开展产品的可靠性增长、维修性/保障性设计等工作的前提。在专门的可靠性试验中，一般能实时、在线监测产品的完好状态：一旦产品发生故障能立刻被发现，因此能获得产品寿命X的准确数值。在获得足够数量的寿命数据后即可分析出产品的寿命分布类型及参数。但在工作场景下，并不见得能针对产品配备在线监测设备，因而不能实时监测产品的完好状态。工作场景下更常见的做法是定期或不定期对产品做完好性检查。假定产品投入使用时刻为零时刻，如果在检查时刻Tc，产品状态为完好，这意味着该产品的寿命X大于Tc；如果在检查时刻Tc，产品状态为故障，这意味着该产品的寿命X小于Tc。与有寿命准确数值的X相比，[检查时刻Tc 状态(完好或故障)]是删失了部分信息的寿命数据。目前，在理论上还没有利用这种删失型数据准确估计寿命分布参数的方法。

相对于标准的可靠性试验场景，当工作环境、使用方式等发生变化时，产品的实际寿命往往可能随之发生变化，因此，即便掌握了产品在可靠性试验场景下的寿命分布规律，也仍有必要去了解产品在工作场景下的实际寿命分布规律。

发明内容

本发明针对现有技术中存在的技术问题，提出了一种利用删失型数据估计寿命分布参数的近似方法，其估计精度能满足工程要求。

产品由各种单元组成。机电件寿命一般服从威布尔型分布，如：滚珠轴承、继电器、蓄电池、液压泵、齿轮、材料疲劳件等，该分布适于描述老化导致的故障。威布尔型单元指寿命服从威布尔分布的单元，寿命X分布记作X～W(α,b)，其中尺度参数α＞0，在工程上形状参数b≥1，X的密度函数为

假定：单元投入使用时刻为零时刻，同批次的单元同时投入使用，且各批次单元的工作场景相似。在第i次检查时，检查时刻记为Tc_i，该批次产品中，完好单元的数量记为Nr_i，故障单元的数量记为Nf_i。共完成了k次检查。

基于该假定，本发明解决上述技术问题的技术方案如下：一种威布尔型单元寿命分布参数估计方法。

本方法包含以下步骤：

步骤1，根据威布尔型单元的使用寿命数据生成n组候选的分布参数(α_j,b_j),1≤j≤n，其中，α_j表示威布尔分布的尺度参数，b_j表示威布尔分布的形状参数，n为正整数；

步骤2，初始化似然度P_j，令

步骤3，根据k次检查结果获取的完好单元的数量Nr_i、故障单元的数量Nf_i以及检查时刻Tc_i，依次更新似然度P_j；

步骤4，在更新后的似然度P_j(1≤j≤n)中找到最大似然度，记为P_M，则似然度P_M对应的α_M、b_M分别为威布尔分布尺度参数和形状参数的估计结果。

具体的，所述步骤1包括：

步骤1.1，确定威布尔分布的尺度参数α1_j1＝α_min+(j1-1)d1,1≤j1≤n1，其中，

α_max表示威布尔型单元寿命分布的尺度参数上限，α_min表示威布尔型单元寿命分布的尺度参数下限，n1为正整数，且n1≥2；

步骤1.2，确定威布尔分布的形状参数b1_j2＝b_min+(j2-1)d2,1≤j2≤n2，其中，

b_max表示威布尔型单元寿命分布的形状参数上限，b_min表示威布尔型单元寿命分布的形状参数下限，n2为正整数，且n2≥2；

步骤1.3，取n＝n1×n2，由α1_j1和b1_j2进行遍历组合获得n组候选的分布参数(α_j,b_j),1≤j≤n。

优选地，步骤1.3中所述的遍历采用以下方式实现：

令j＝1；

第一层循环中遍历j1＝1:n1，第二层循环中遍历j2＝1:n2，

令

α_j＝α1_j1；b_j＝b1_j2；j＝j+1；

其中，α_max≥α1_j1≥α_min，b_max≥b1_j2≥b_min。

进一步，所述步骤3具体包括：

步骤3.1，令i＝1，i表示检查次数；

步骤3.2，遍历计算W_j，1≤j≤n，令

其中

Tc_i表示第i次检查时的检查时刻，Nr_i为第i次检查时完好单元的数量，Nf_i为第i次检查时故障单元的数量；

步骤3.3，遍历更新似然度P_j，令

步骤3.4，令i＝i+1，若i≤k则转步骤3.2，否则转步骤4，k为总检查次数。

本发明的有益效果是：本方法的参数估计结果在总体上能“跟随”理论成熟方法的参数估计结果，其估计精度能满足工程要求。

附图说明

图1为本发明方法流程图；

图2为采用理论成熟方法和本发明方法的仿真结果示意图。

具体实施方式

以下结合实施例对本发明的原理和特征进行描述，所举实例只用于解释本发明，并非用于限定本发明的范围。

产品由各种单元组成。机电件寿命一般服从威布尔型分布，如：滚珠轴承、继电器、蓄电池、液压泵、齿轮、材料疲劳件等，该分布适于描述老化导致的故障。威布尔型单元指寿命服从威布尔分布的单元，寿命X分布记作X～W(α,b)，其中尺度参数α＞0，在工程上形状参数b≥1，X的密度函数为

假定：单元投入使用时刻为零时刻，同批次的单元同时投入使用，且各批次单元的工作场景相似。在第i次检查时，检查时刻记为Tc_i，该批次产品中，完好单元的数量记为Nr_i，故障单元的数量记为Nf_i。共完成了k次检查。

基于该假定：

实施例1

本方法，如图1所示，步骤如下：

1、确定候选的寿命分布参数

根据威布尔型单元的使用寿命数据生成n组候选的分布参数(α_j,b_j),1≤j≤n，其中，α_j表示威布尔分布的尺度参数，b_j表示威布尔分布的形状参数，n为正整数；

生成候选分布参数的具体实现方法为：

步骤1.1，确定威布尔分布的尺度参数α1_j1＝α_min+(j1-1)d1,1≤j1≤n1，其中，

α_max表示威布尔型单元寿命分布的尺度参数上限，α_min表示威布尔型单元寿命分布的尺度参数下限，n1为正整数，且n1≥2；

步骤1.2，确定威布尔分布的形状参数b1_j2＝b_min+(j2-1)d2,1≤j2≤n2，其中，

b_max表示威布尔型单元寿命分布的形状参数上限，b_min表示威布尔型单元寿命分布的形状参数下限，n2为正整数，且n2≥2；

步骤1.3，取n＝n1×n2，由α1_j1和b1_j2进行遍历组合获得n组候选的分布参数(α_j,b_j),1≤j≤n。其中遍历的方式可通过以下方式实现：

令j＝1；

第一层循环中遍历j1＝1:n1，第二层循环中遍历j2＝1:n2，

令

α_j＝α1_j1；b_j＝b1_j2；j＝j+1；

其中，α_max≥α1_j1≥α_min，b_max≥b1_j2≥b_min。

2、初始化似然度

初始化似然度P_j,1≤j≤n，令

3、遍历调整似然度

步骤3.1，令i＝1，i表示检查次数；

步骤3.2，遍历计算W_j，1≤j≤n，令

其中

Tc_i表示第i次检查时的检查时刻，Nr_i为第i次检查时完好单元的数量，Nf_i为第i次检查时故障单元的数量；

步骤3.3，遍历更新似然度P_j，令

步骤3.4，令i＝i+1，若i≤k则转步骤3.2，否则转步骤4，k为总检查次数。

4、输出寿命分布参数估计结果

在所有似然度P_j(1≤j≤n)中找到最大似然度，记为P_M，则最大似然度P_M对应的α_M、b_M分别为威布尔分布尺度参数和形状参数的估计结果。

实施例2

某威布尔型单元的10次状态检查结果如表1所示，试估计其寿命分布的形状参数和尺度参数。

表1

检查序号	检查时刻h	故障单元的数量	完好单元的数量
				1	630	3	2
2	250	1	4
				3	880	3	2
4	500	0	5
				5	130	0	5
6	1250	4	1
				7	1000	3	2
8	380	0	5
				9	1130	3	2
10	750	1	4

计算过程如下：

1、确定候选的寿命分布参数

按照以往经验，估计该单元的尺度参数在500～3000范围内，以500为步长；估计该单元的形状参数在1.1～2.6范围内，以0.5为步长；共生成24个候选的分布参数(α_j,b_j),1≤j≤24。

2、初始化似然度

初始化似然度P_j,1≤j≤24，令

3、遍历调整似然度

3.1令i＝1

3.2遍历计算W_j,1≤j≤24，令

其中

3.3遍历更新似然度P_j,1≤j≤24，令

3.4更新i，令i＝i+1，若i≤10则转3.2，否则转4。表2列出了第i次更新后的似然度。

4、输出寿命分布参数估计结果

在所有似然度P_j(1≤j≤24)中最大似然度为P₇，则α₇＝1000,b₇＝2.1分别为威布尔分布尺度参数和形状参数的估计结果。

表2第i次更新后的似然度

实施例3

可建立以下仿真模型来模拟对单元的检查过程。

假定单元的实际寿命服从威布尔分布W(α,b)，共进行k次检查，记第i次检查时刻为Tc_i，同批次的单元同时投入使用，第i批次的单元数量为N_i。

1)令i＝1

2)随机产生N_i个随机数simT_ij,1≤j≤N_i，这些随机数服从威布尔分布W(α,b)。

3)在simT_ij(1≤j≤N_i)中，找到大于Tc_i的随机数，其数量为完好单元的数量记为Nr_i，故障单元数量Nf_i为N_i-Nr_i。

4)更新i，令i＝i+1。若i≤k则转2)，否则本次模拟k次检查结束。

对于以上仿真模型得到的Tc_i、Nr_i、Nf_i，可采用本发明方法用于估计单元寿命的分布参数。

对于以上仿真模型得到的simT_ij，可采用理论上成熟的方法用于估计单元寿命的分布参数。

以单元的实际寿命服从威布尔分布W(1000,2.1)为例，共进行20次检查，同批次的单元同时投入使用，第i批次的单元数量为6，利用上述仿真模型开展多次仿真，得到大量模拟检查结果并进行分布参数估计，并对多次估计结果进行统计分析。

针对模拟的删失型检查结果数据Tc_i、Nr_i、Nf_i，采用本文方法估计的威布尔分布尺度参数α的均值为1013.0、根方差为66.1，形状参数b的均值为2.20、根方差为0.39。

针对模拟的寿命数据simT_ij，采用理论方法估计的威布尔分布尺度参数α的均值为1000.5、根方差为45.2，形状参数b的均值为2.12、根方差为0.16。

图2显示了利用10次仿真结果，分别针对寿命型数据和对应的删失型数据，各自采用理论成熟方法和本文方法得到的寿命分布参数估计结果。从图2来看，本文方法的参数估计结果在总体上能“跟随”理论成熟方法的参数估计结果。

大量仿真验证结果表明，本文方法有较好的估计精度，满足工程应用要求。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (3)

1.一种威布尔型单元寿命分布参数估计方法，其特征在于，包括以下步骤：

步骤1，根据威布尔型单元的使用寿命数据生成n组候选的分布参数(α_j,b_j),1≤j≤n，其中，α_j表示威布尔分布的尺度参数，b_j表示威布尔分布的形状参数，n为正整数；

步骤2，初始化似然度P_j，令

步骤3，根据k次检查结果获取的完好单元的数量Nr_i、故障单元的数量Nf_i以及检查时刻Tc_i，依次更新似然度P_j；

步骤4，在更新后的似然度P_j中找到最大似然度，记为P_M，则似然度P_M对应的α_M、b_M分别为威布尔分布尺度参数和形状参数的估计结果；

所述步骤3具体包括：

步骤3.1，令i＝1，i表示检查次数；

步骤3.2，遍历计算W_j，1≤j≤n，令

其中

Tc_i表示第i次检查时的检查时刻，Nr_i为第i次检查时完好单元的数量，Nf_i为第i次检查时故障单元的数量；

步骤3.3，遍历更新似然度P_j，令

步骤3.4，令i＝i+1，若i≤k则转步骤3.2，否则转步骤4，k为总检查次数。

2.根据权利要求1所述一种威布尔型单元寿命分布参数估计方法，其特征在于，所述步骤1具体包括：

步骤1.1，确定威布尔分布的尺度参数α1_j1＝α_min+(j1-1)d1,1≤j1≤n1，其中，

α_max表示威布尔型单元寿命分布的尺度参数上限，α_min表示威布尔型单元寿命分布的尺度参数下限，n1为正整数，且n1≥2；

步骤1.2，确定威布尔分布的形状参数b1_j2＝b_min+(j2-1)d2,1≤j2≤n2，其中，

b_max表示威布尔型单元寿命分布的形状参数上限，b_min表示威布尔型单元寿命分布的形状参数下限，n2为正整数，且n2≥2；

步骤1.3，取n＝n1×n2，由α1_j1和b1_j2进行遍历组合获得n组候选的分布参数(α_j,b_j),1≤j≤n。

3.根据权利要求2所述一种威布尔型单元寿命分布参数估计方法，其特征在于，步骤1.3中所述的遍历采用以下方式实现：

令j＝1；

第一层循环中遍历j1＝1:n1，第二层循环中遍历j2＝1:n2，

令

α_j＝α1_j1；b_j＝b1_j2；j＝j+1；

其中，α_max≥α1_j1≥α_min，b_max≥b1_j2≥b_min。

Images