本发明涉及可靠性测试技术领域，特别涉及一种估计寿命服从伽玛型分布单元可靠性参数的方法。
背景技术：
可靠性是描述产品质量的核心属性，通常用寿命的分布规律(分布类型和参数)来定量描述可靠性。理论上，针对产品开展大量的可靠性试验，可以获得足够数量的产品寿命数据，然后可采用成熟的数理统计方法来估计出产品寿命的分布类型和参数。但在实际工作中，针对产品开展大量可靠性试验，往往意味着高昂的经济成本和漫长的试验耗时，因此更常见的做法是利用“少量的可靠性试验数据+在产品研制、生产、使用等阶段产生的大量数据”，估算产品寿命的分布规律。在产品的可靠性试验中，一般配备有专门的在线检测设备，用于实时监测产品的完好性状态，及时记录产品的故障时刻，因此可以获得产品寿命的数值。但在产品研制、生产、使用等这些非可靠性试验场景下，不一定配备有专门的在线检测设备，只能定期或不定期地对产品进行完好性检查，因而不能准确获知产品的故障时刻，也就无法获得寿命的数值信息。
技术实现要素：
为了克服背景技术中存在的缺陷，本发明提供一种估计伽玛型单元可靠性参数的方法。
为了实现上述目的，本发明采用的技术方案为：一种估计伽玛型单元可靠性参数的方法，包括以下步骤：
步骤一：确定候选的寿命分布参数，根据已有伽玛型单元可靠性数据分布规律，初步确定待估伽玛型单元寿命伽玛分布的形状参数上限αmax和形状参数下限αmin及待估伽玛型单元寿命伽玛分布的尺度参数上限bmax和尺度参数下限bmin；
在确定的伽玛型单元伽玛形状参数区间内，等间隔生成n1个候选形状参数，候选参数中每个相邻参数之间的步长相等为d1，根据伽玛分布形状参数的步长依次计算出n1个伽玛分布的形状参数α1j1，其中1≤j1≤n1；
在确定的伽玛型单元伽玛尺度参数区间内，等间隔生成n2个候选尺度参数，候选参数中每个相邻参数之间的步长相等为d2，根据伽玛分布尺度参数的步长依次计算出n2个伽玛分布的尺度参数b1j2，其中1≤j2≤n2；
取n＝n1×n2，由α1j1和b1j2进行遍历组合获得n组候选的分布参数(αj,bj),1≤j≤n；
步骤二：遍历寿命形状参数为αj和尺度参数为bj的组合计算似然值，对于每个候选参数组合，针对一组包含m个伽玛型单元检测信息的数据组，根据第i个检测信息所包含的在ti时刻的单元状态信息fi确定其对应的计算系数wi，依据m个检测数据不断迭代更新似然值lj，迭代之初，设定每个候选参数组合对应的似然值初始值为0，在每个候选参数组合所对应的迭代完毕后的似然值中，寻找最大值即lm，则该最大值对应的αm为伽玛型单元形状参数的估计值，bm为尺度参数的估计值。
在上述方案中，所述骤一中候选的寿命分布参数的计算过程如下：
1)确定伽玛分布的形状参数α1j1和步长d1的计算过程如下：
其中，αmax表示伽玛分布的形状参数上限，αmin表示伽玛分布的形状参数下限，n1为正整数，且n1≥2；
2)确定伽玛分布的尺度参数b1j2和步长d2的计算过程如下：
其中，bmax为伽玛分布的尺度参数上限，bmin为伽玛分布的尺度参数下限，n2为正整数，且n2≥2。
3)α1j1和b1j2遍历组合的计算方式为：
令j＝1，第一层循环中遍历j1＝1:n1，第二层循环中遍历j2＝1:n2；令αj＝α1j1，bj＝b1j2，j＝j+1；其中αmax≥α1j1≥αmin，bmax≥b1j2≥bmin。
在上述方案中，所述步骤二中计算系数wi和似然值lj的计算公式如下：
其中，wi为计算系数，log(*)为自然对数函数，αj为伽玛分布的形状参数，bj为伽玛分布的尺度参数，γ(α)为伽玛函数，且
在上述方案中，所述步骤二中似然值lj遍历计算过程如下：
1)令j＝1；
2)令i＝1，lj＝0；
3)计算系数
其中，log(*)为自然对数函数，αj为伽玛分布的形状参数，bj为伽玛分布的尺度参数，γ(α)为伽玛函数，且
4)更新i＝i+1，若i≤m则转3)，否则转5)；
5)更新j＝j+1，若j≤n则转2)，否则6)；
6)在lj(1≤j≤n)中找到最大值，记其为lm，则αm为伽玛型单元形状参数的估计值，bm为尺度参数的估计值。
与现有技术相比，本发明的有益效果是：利用“少量的可靠性试验数据+在产品研制、生产、使用等阶段产生的大量数据”，估算产品寿命的分布规律，避免针对产品开展大量的可靠性试验造成的人力、物力及财力的消耗。
具体实施方式
以下结合某伽玛型单元案例，对本发明作进一步的详细描述。
本发明一种估计伽玛型单元可靠性参数的方法，包括以下步骤：
步骤一：确定候选的寿命分布参数，根据已有伽玛型单元可靠性数据分布规律，初步确定待估伽玛型单元伽玛分布的形状参数上限αmax和形状参数下限αmin及待估伽玛型单元伽玛分布的尺度参数上限bmax和尺度参数下限bmin；
在确定的伽玛型单元伽玛形状参数区间内，等间隔生成n1个候选形状参数，候选参数中每个相邻参数之间的步长相等为d1，根据伽玛分布形状参数的步长依次计算出n1个伽玛分布的形状参数α1j1，其中1≤j1≤n1；
在确定的伽玛型单元伽玛尺度参数区间内，等间隔生成n2个候选尺度参数，候选参数中每个相邻参数之间的步长相等位d2，根据伽玛分布尺度参数的步长依次计算出n2个伽玛分布的尺度参数b1j2，其中1≤j2≤n2；
取n＝n1×n2，由α1j1和b1j2进行遍历组合获得n组候选的分布参数(αj,bj),1≤j≤n；
步骤二：遍历形状参数为αj和尺度参数为bj的组合计算似然值，对于每个候选参数组合，针对一组包含m个伽玛型单元检测信息的数据组，根据第i个检测信息所包含的在ti时刻的单元状态信息fi确定其对应的计算系数wi，依据m个检测数据不断迭代更新似然值lj，迭代之初，设定每个候选参数组合对应的似然值初始值为0，在每个候选参数组合所对应的迭代完毕后的似然值中，寻找最大值记为lm，则该最大值对应的αm为伽玛型单元形状参数的估计值，bm为尺度参数的估计值。
在上述方案中，所述骤一中候选的寿命分布参数的计算过程如下：
1)确定伽玛分布的形状参数α1j1和步长d1的计算过程如下：
其中，αmax表示伽玛分布的形状参数上限，αmin表示伽玛分布的形状参数下限，n1为正整数，且n1≥2。
2)确定伽玛分布的尺度参数b1j2和步长d2的计算过程如下：
其中，bmax为伽玛分布的尺度参数上限，bmin为伽玛分布的尺度参数下限，n2为正整数，且n2≥2。
3)α1j1和b1j2遍历组合的计算方式为：
令j＝1，第一层循环中遍历j1＝1:n1，第二层循环中遍历j2＝1:n2；令αj＝α1j1，bj＝b1j2，j＝j+1；其中αmax≥α1j1≥αmin，bmax≥b1j2≥bmin。
在上述方案中，所述步骤二中计算系数wi和似然值lj的计算公式如下：
其中，wi为计算系数，log(*)为自然对数函数，log(*)为自然对数函数，αj为伽玛分布的形状参数，bj为伽玛分布的尺度参数，γ(α)为伽玛函数，且
在上述方案中，所述步骤二中似然值lj遍历计算过程如下：
1)令j＝1；
2)令i＝1，lj＝0；
3)计算系数
其中，log(*)为自然对数函数，αj为伽玛分布的形状参数，bj为伽玛分布的尺度参数，γ(α)为伽玛函数，且
4)更新i＝i+1，若i≤m则转3)，否则转5)；
5)更新j＝j+1，若j≤n则转2)，否则6)；
6)在lj(1≤j≤n)中找到最大值，记其为lm，则αm为伽玛型单元形状参数的估计值，bm为尺度参数的估计值。
实施例，某伽玛型单元的[ft]型可靠性数据如下表，试估计其形状参数和尺度参数。
从以往经验得知，该伽玛型单元的形状参数在1.1～3.9范围内，以0.4为步长；尺度参数在500～2900范围内，以400为步长，共生成56个候选的分布参数(αj,bj),1≤j≤56，计算结果如下表所示：
从上表中可以看出，在lj(1≤j≤56)中的最大值为l25，则α25＝2.3，b25＝1700为该伽玛型单元寿命分布参数的估计值。
为进一步验证本发明方法的可行性，建立了如下仿真模型。
假定某伽玛型单元的寿命服从伽玛分布ga(α,b)。
(1)产生k1个随机数simti(1≤i≤k1)，simti服从伽玛分布ga(α,b)，用于模拟单元的寿命值。令fi＝0,ti＝simti得到k1组[fiti],1≤i≤k1。
(2)产生k2个随机数simti(k1+1≤i≤k1+k2)，simti服从伽玛分布ga(α,b)，用于模拟单元的寿命值。
(3)产生k2个均匀随机数simtci(k1+1≤i≤k1+k2)，用于模拟检查时刻。
(4)在k1+1≤i≤k1+k2范围内，令
利用以上仿真得到的k1+k2组寿命数据[fiti]后，再应用本文方法可得到分布参数的估计值。以α＝2.3,b＝1700,k1＝5,k2＝15为例，大量仿真后用本文方法得到的该伽玛型单元形状参数统计结果的均值为2.68、根方差为0.94，尺度参数统计结果的均值为1709.0、根方差为718.7。如果使用k1+k2组寿命数据simti的话，采用理论方法计算得到的形状参数统计结果的均值为2.77、根方差为0.97，尺度参数统计结果的均值为1611.2、根方差为538.8，二者的差异在工程允许范围以内。
以上仅为本发明的实施方式，并非因此限制本发明的专利范围，凡是利用本发明说明书内容所作的等效结构或等效流程变换，或直接或间接运用在其他的相关技术领域，均同理包括在本发明的专利保护范围内。