本发明涉及产品可靠性预测技术领域，具体涉及一种基于备件保障数据的伽玛型单元寿命分布参数估计方法。
背景技术：
产品由各种单元组成。通常以产品的寿命分布函数来定量描述产品的可靠性。伽玛分布常用来描述类似“冲击”引起的故障，假若单元能经受若干次外界冲击，但当单元受冲击次数累积到一定次数时就产生故障。例如电网中存在着电涌现象，一些电子器件当承受的电涌冲击次数超过一定数量时会发生故障。伽玛型单元指寿命服从伽玛分布的单元，寿命x服从伽玛分布记作x～ga(α,b)，其中α＞0为形状参数，b＞0为尺度参数，x的密度函数为
寿命分布函数由函数类别和分布参数决定，本质上是一种描述统计规律的函数。因此，只有当寿命样本值达到一定数量规模时，该函数才能较为准确的描述产品寿命的分布规律。在实际工作中，工作环境往往对产品可靠性有较大影响，因此需要关注工作环境下产品的寿命分布规律。在实际工作中，产品工况数据的记录质量与专门的可靠性试验相比往往较差，常常由于记录不及时或记录不全，导致故障发生时间、产品累积工作时间等直接反映产品寿命值的相关数据丢失，只留有保障任务开始前的备件数量、任务成功与否、任务计划时间等数据。这使得基于寿命值数据的分布参数常规估计方法难以从这种备件保障数据中分析产品的可靠性。
技术实现要素：
本发明针对现有技术中存在的技术问题，提供一种基于备件保障数据的伽玛型单元寿命分布参数估计方法，该方法约定：共有k次保障任务的相关数据，第i次保障任务时该单元的计划工作时间记为twi(以下简称为任务时间)，任务开始前的备件数量记为nbji，保障任务结果记为fi，1≤i≤k。在保障任务期间，以备件更换故障件的方式来排除故障，若发生的故障次数不大于备件数量，则本次保障任务成功(fi记为1)，否则视为失败(fi记为0)。即：以形如[twinbjifi]的备件保障数据来描述一次保障任务。在积累了k组备件保障数据后，估计寿命分布参数。
本发明解决上述技术问题的技术方案如下：
一种基于备件保障数据的伽玛型单元寿命分布参数估计方法，包括以下步骤：
步骤1，根据k次保障任务的相关数据生成n组候选分布参数(αj,bj),1≤j≤n，其中，αj表示伽玛分布的形状参数，bj表示伽玛分布的尺度参数，k、n为正整数；
步骤2，初始化似然值lj，令lj＝0，1≤j≤n；
步骤3，遍历候选参数，依次计算每次保障任务的概率p1；并利用所述k次保障任务的概率计算更新似然值lj；
步骤4，在更新后的似然值lj(1≤j≤n)中找到最大似然值，记为lm，则似然值lm对应的αm、bm分别为该单元寿命分布参数的估计值。
进一步的，所述步骤1包括：
步骤1.1，确定伽玛分布的形状参数α1j1＝αmin+(j1-1)d1,1≤j1≤n1，其中，
步骤1.2，确定伽玛分布的尺度参数b1j2＝bmin+(j2-1)d2,1≤j2≤n2，其中，
步骤1.3，取n＝n1×n2，由α1j1和b1j2进行遍历组合获得n组候选的分布参数(αj,bj),1≤j≤n。
进一步的，步骤1.3中所述的遍历通过以下方式实现：
令j＝1；
第一层循环中遍历j1＝1:n1，第二层循环中遍历j2＝1:n2，
令αj＝α1j1；bj＝b1j2；j＝j+1；
其中，αmax≥α1j1≥αmin，bmax≥b1j2≥bmin。
进一步的，所述步骤3包括：
步骤3.1，令候选参数序号j＝1；
步骤3.2，令任务序号i＝1，1≤i≤k，k为保障任务最大次数；
步骤3.3，对于第i次保障任务，利用下式计算保障任务概率p1i；
其中，twi为第i次保障任务时该单元的计划工作时间，nbji为任务开始前的备件数量；fi表示保障任务结果，保障任务成功则fi＝1，保障任务失败则fi＝0；
步骤3.4，更新i＝i+1，若i≤k则跳转至步骤3.3，否则跳转至步骤3.5；
步骤3.5，令
步骤3.6，令j＝j+1，若j≤n则跳转至步骤3.2，否则跳转至步骤4。
本发明的有益效果是：本发明提出的寿命分布参数估计方法，能在因缺少专门的可靠性试验导致缺少寿命值数据时，有效利用实际工作产生的备件保障数据，较为准确地估计寿命分布参数。
具体实施方式
以下结合实施例对本发明的原理和特征进行描述，所举实例只用于解释本发明，并非用于限定本发明的范围。
产品由各种单元组成。伽玛分布常用来描述类似“冲击”引起的故障，假若单元能经受若干次外界冲击，但当单元受冲击次数累积到一定次数时就产生故障。例如电网中存在着电涌现象，一些电子器件当承受的电涌冲击次数超过一定数量时会发生故障。伽玛型单元指寿命服从伽玛分布的单元，寿命x服从伽玛分布记作x～ga(α,b)，其中α＞0为形状参数，b＞0为尺度参数，x的密度函数为
约定：共有k次保障任务的相关数据，第i次保障任务时该单元的计划工作时间记为twi(以下简称为任务时间)，任务开始前的备件数量记为nbji，保障任务结果记为fi，1≤i≤k。在保障任务期间，以备件更换故障件的方式来排除故障，若发生的故障次数不大于备件数量，则本次保障任务成功(fi记为1)，否则视为失败(fi记为0)。即：以形如[twinbjifi]的备件保障数据来描述一次保障任务。在积累了k组备件保障数据后，具体估计寿命分布参数的方法如下：
1确定候选的寿命分布参数
生成n组候选分布参数(αj,bj),1≤j≤n，其中，αj表示伽玛分布的形状参数，bj表示伽玛分布的尺度参数，n为正整数；
生成候选分布参数的具体实现方法为：
1)确定伽玛分布的形状参数α1j1＝αmin+(j1-1)d1,1≤j1≤n1，其中，
2)确定伽玛分布的尺度参数b1j2＝bmin+(j2-1)d2,1≤j2≤n2，其中，
3)取n＝n1×n2，由α1j1和b1j2进行遍历组合获得n组候选的分布参数(αj,bj),1≤j≤n。其中，遍历的方式可以通过以下方式实现：
令j＝1；第一层循环中遍历j1＝1:n1，第二层循环中遍历j2＝1:n2，令αj＝α1j1；bj＝b1j2；j＝j+1。其中，αmax≥α1j1≥αmin，bmax≥b1j2≥bmin。
2初始化似然值
初始化似然值lj,1≤j≤n，令lj＝0。
3、遍历计算似然值
3.1令候选参数序号j＝1；
3.2令任务序号i＝1；
3.3对于第i次保障任务，计算保障任务概率p1i：
3.4更新i＝i+1，若i≤k则转3.3，否则转3.5；
3.5令
3.6令j＝j+1，若j≤n则转3.2，否则转4。
4、输出结果
在lj(1≤j≤n)中找到最大值，记其为lm，则αm、bm为该单元寿命分布参数的估计值。
实施例
某单元的[twinbjifi]型备件保障数据如表1所示，从工程经验可知该单元的寿命服从伽玛分布ga(α,b)，试估计其寿命分布参数。
表1某单元的备件保障数据
计算过程如下：
1确定候选的寿命分布参数
从以往经验得知，该单元的形状参数在1.1～2.7范围内，以0.4为步长；尺度参数在100～1700范围内，以400为步长，共生成25个候选的分布参数αj、bj,1≤j≤25。
2初始化似然值
初始化似然值lj,1≤j≤25，令lj＝0。
3、遍历计算似然值
3.1令候选参数序号j＝1
3.2令任务序号i＝1
3.3对于第i次保障任务，计算保障任务概率p1i：
3.4更新i＝i+1，若i≤15则转3.3，否则转3.5。
表2是j＝1～25时，15次保障任务概率结果。
表2保障任务概率p1i
3.5令
表3似然值
3.6令j＝j+1，若j≤25则转3.2，否则转4。
4、输出结果
从表3可知，α13＝1.9、b13＝900对应的似然值最大，因此将其作为该单元的寿命分布参数估计值。
可建立以下仿真模型模拟备件保障过程：
1)产生1+nbj个随机数simti,1≤i≤1+nbj用于模拟装备该单元和nbj个备件的寿命，simti服从该类单元的寿命分布规律；
2)令累积时间
3)若simtw≥tw，则本次保障任务成功，令f＝1；否则，本次保障任务失败，令f＝0。
运行上述仿真模型，可以得到大量形如[twinbjifi]的备件保障数据，对本发明的方法进行仿真验证。以α＝1.9,b＝900,k＝50为例进行大量仿真。如果基于备件保障数据，用本文方法得到的该单元的寿命分布形状参数统计结果的均值为2.10、根方差为0.68，尺度参数统计结果的均值为926.4、根方差为363.6。如果基于寿命数据simti，采用理论方法计算得到的形状参数统计结果的均值为1.93、根方差为0.19，尺度参数统计结果的均值为895.8、根方差为99.6。二者的差异在工程允许范围以内。
以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。