CN109388860B - 一种伽玛型单元寿命分布参数估计方法 - Google Patents

Abstract

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

Description

一种伽玛型单元寿命分布参数估计方法

技术领域

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

背景技术

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

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

发明内容

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

伽玛分布常用来描述类似“冲击”引起的故障，假若单元能经受若干次外界冲击，但当单元受冲击次数累积到一定次数时就产生故障。例如电网中存在着电涌现象，一些电子器件当承受的电涌冲击次数超过一定数量时会发生故障。伽玛型单元指寿命服从伽玛分布的单元，寿命X服从伽玛分布记作X～Ga(α,b)，其中α＞0为形状参数，b＞0为尺度参数，X的密度函数为

其中Γ(α)为伽玛函数，且

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

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

本方法包含以下步骤：

步骤1，根据伽玛型单元使用寿命数据生成n组候选的分布参数(α2_j,b2_j),1≤j≤n，其中，α2_j表示伽玛分布的形状参数，b2_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对应的α2_M、b2_M分别为伽玛分布形状参数和尺度参数的估计结果。

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

步骤1.1，确定伽玛分布的形状参数α_j1＝α_min+(j1-1)d1,1≤j1≤n1，其中，

α_max表示伽玛型单元寿命分布的形状参数上限，α_min表示伽玛型单元寿命分布的形状参数下限，n1为正整数，且n1≥2；

步骤1.2，确定伽玛分布的尺度参数b_j2＝b_min+(j2-1)d2,1≤j2≤n2，其中，

b_max表示伽玛型单元寿命分布的尺度参数上限，b_min表示伽玛型单元寿命分布的尺度参数下限，n2为正整数，且n2≥2；

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

进一步，步骤1.3中所述的遍历采用以下方式实现：

令j＝1；

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

令

α2_j＝α_j1；b2_j＝b_j2；j＝j+1；

其中，α_max≥α_j1≥α_min，b_max≥b_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～Ga(α,b)，其中α＞0为形状参数，b＞0为尺度参数，X的密度函数为

其中Γ(α)为伽玛函数，且

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

基于该假定：

实施例1

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

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

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

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

1.1，确定伽玛分布的形状参数α_j1＝α_min+(j1-1)d1,1≤j1≤n1，其中，

α_max表示伽玛型单元寿命分布的形状参数上限，α_min表示伽玛型单元寿命分布的形状参数下限，n1为正整数，且n1≥2；

1.2，确定伽玛分布的尺度参数b_j2＝b_min+(j2-1)d2,1≤j2≤n2，其中，

b_max表示伽玛型单元寿命分布的尺度参数上限，b_min表示伽玛型单元寿命分布的尺度参数下限，n2为正整数，且n2≥2；

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

其中遍历的方式采用以下方式实现：

令j＝1；

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

令

α2_j＝α_j1；b2_j＝b_j2；j＝j+1；

其中，α_max≥α_j1≥α_min，b_max≥b_j2≥b_min。

2、初始化似然度

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

3、遍历调整似然度

3.1令i＝1

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

其中

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

3.4更新i，令i＝i+1，若i≤k则转步骤3.2，否则转步骤4。

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

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

实施例2

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

表1

检查序号	检查时刻h	故障单元的数量	完好单元的数量
				1	2330	2	3
2	780	1	4
				3	2590	5	0
4	1560	4	1
				5	2070	3	2
6	520	1	4
				7	260	0	5
8	1040	2	3
				9	1820	3	2
10	1300	1	4

计算过程如下：

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

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

2、初始化似然度

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

3、遍历调整似然度

3.1令i＝1

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

其中

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

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

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

在所有似然度P_j(1≤j≤28)中最大似然度为P₁₀，则α2₁₀＝1.7,b2₁₀＝1000分别为伽玛分布形状参数和尺度参数的估计结果。

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

实施例3

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

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

1)令i＝1

2)随机产生N_i个随机数simT_ij,1≤j≤N_i，这些随机数服从伽玛分布Ga(α,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，可采用理论上成熟的方法用于估计单元寿命的分布参数。

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

针对模拟的删失型检查结果数据Tc_i、Nr_i、Nf_i，采用本文方法估计的伽玛分布形状参数α的均值为1.68、根方差为0.40，尺度参数b的均值为1080.0、根方差为272.7。

针对模拟的寿命数据simT_ij，采用理论方法估计的伽玛分布形状参数α的均值为1.72、根方差为0.21，尺度参数b的均值为996.0、根方差为131.1。

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

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

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

Claims (3)

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

步骤1，根据伽玛型单元使用寿命数据生成n组候选的分布参数(α2_j,b2_j),1≤j≤n，其中，α2_j表示伽玛分布的形状参数，b2_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对应的α2_M、b2_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，确定伽玛分布的形状参数α_j1＝α_min+(j1-1)d1,1≤j1≤n1，其中，

α_max表示伽玛型单元寿命分布的形状参数上限，α_min表示伽玛型单元寿命分布的形状参数下限，n1为正整数，且n1≥2；

步骤1.2，确定伽玛分布的尺度参数b_j2＝b_min+(j2-1)d2,1≤j2≤n2，其中，

b_max表示伽玛型单元寿命分布的尺度参数上限，b_min表示伽玛型单元寿命分布的尺度参数下限，n2为正整数，且n2≥2；

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

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

令j＝1；

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

令

α2_j＝α_j1；b2_j＝b_j2；j＝j+1；

其中，α_max≥α_j1≥α_min，b_max≥b_j2≥b_min。

Images