基于多阶段试验数据折合的航天产品可靠性评估方法与流程

本发明属于可靠性评估领域，具体涉及一种基于折合因子法对多阶段试验中出现无失效数据的航天产品进行可靠性评估的方法。

背景技术：

火箭发动机涡轮泵、轴承等的航天产品具有高可靠性、小子样的特点，同时又受到时间和预算的限制，因此在可靠性试验过程中可能出现极少数失效甚至无失效发生。此外，工程实践中通常在完成一个阶段的试验之后对产品进行设计改进，弥补产品缺陷，提高产品可靠性。因此，不同阶段之间的产品属于不同的寿命分布母体，其寿命分布的参数是变化的。

自二十世纪中期提出weibull分布以来，其应用遍布在工程制造、可靠性、失效分析等领域。假设航天产品寿命服从weibull分布，其累积分布函数、概率密度函数以及可靠度函数如下所示：

其中，t表示工作时间，m是与产品材料相关的形状参数，η是尺度参数，也被称为特征寿命。

对于多阶段产品可靠性增长和评估，目前已有的方法可以分为两大类：经典参数可靠性模型和可靠性增长bayes模型。经典参数可靠性模型中较为常见的包括duane模型，amsaa模型，fries模型和非齐次几何模型。这些经典方法不能在小样本情况下描述产品可靠性增长，也不能充分利用之前阶段数据提高可靠性评估的精度。现有的可靠性增长bayes模型主要有smith模型，barlow-scheuer模型，logistic模型，折合因子法等。但是现有研究没有充分考虑将之前阶段的数据信息融合到验前分布中以提高评估精度，仅仅进行当前阶段的可靠性评估，也不能预测产品的最终可靠度。

根据是否需要验前信息，无失效数据可靠性评估方法分为两大类：经典分析计算方法和bayes理论方法。经典分析计算方法，例如修正极大似然函数法、样本空间排序法等，无需验前信息就可进行可靠性评估，故而其评估精度可能差强人意。bayes理论方法包括bayes估计、多层bayes估计以及e-bayes估计方法等，该类方法融合多源的验前信息，提高了评估精度，也因此，先验分布的选择成了决定性的也是最为关键的环节。

技术实现要素：

针对现有技术存在的缺陷，本发明提供一种基于多阶段试验数据折合的航天产品可靠性评估方法。本发明基于weibull分布的阶段间数据折合方法，联合折合后的数据和当前阶段数据形成的新数据得到bayes联合验后分布，进而采用蒙特卡洛链抽样法(metropol-hastings(m-h)算法)计算可靠性参数的估计值，并进一步得到可靠度的区间估计，完成可靠性评估，提高火箭发动机涡轮泵、轴承等的航天产品可靠性评估精度。

为实现上述技术目的，本发明采用的技术方案是：

一种基于多阶段试验数据折合的航天产品可靠性评估方法，包括：

步骤一，进行航天产品的两阶段可靠性试验，获取航天产品第一阶段以及第二阶段的产品试验寿命数据；

步骤二，对第一阶段的产品试验寿命数据进行折合；

步骤三，确定第二阶段weibull寿命分布参数的验前分布，并得到联合分布数据下的bayes验后分布核函数；

步骤四，采用基于bayes的m-h抽样方法进行可靠度的区间估计。

进一步地，本发明步骤一中，航天产品第一阶段可靠性试验中，有n1个样本投入试验，在截尾时间ta时出现r个失效，其中失效数据为t1，t2，...，tr，无失效数据为tr+1＝tr+2＝…＝tn1＝ta，此阶段的weibull寿命分布参数为(m1，η1)。

航天产品第二阶段可靠性试验中，有n2个样本投入试验，在截尾时间t2时无失效数据产生，此阶段的weibull寿命分布参数为(m2，η2)。

进一步地，本发明步骤二，包括：

(2.1)对于第一阶段的产品试验寿命数据，已知其失效数据为t1，t2，...，tr以及形状参数m1的值，使用极大似然估计方法得到尺度参数η1的估计值；

两变量的第一阶段的产品试验寿命数据的极大似然函数为：

式中，d1表示第一阶段的产品试验寿命数据；l(d1|m1，η1)表示第一阶段的产品试验寿命数据的极大似然函数，其中包括两个变量m1和η1。

根据极大似然估计方法的原理，得到η1的估计值；

完成对产品第一阶段可靠性试验的可靠性评估，其寿命分布函数为已知。

(2.2)计算第一阶段产品的失效概率；

其中：p1j表示在时间t＝tj(j＝1，...，n1)时第一阶段产品的失效概率。

(2.3)计算第二阶段产品的失效概率；

其中p2表示第二阶段产品的失效概率。

(2.4)进行数据折合，得到联合数据；

tj→t1j＝tjcj

其中cj(j＝1，...，n1)表示t＝tj时第一阶段数据向第二阶段数据折合的折合因子；第一阶段折合后的数据为t11，t12，...，t1r，t1a，...，t1a，与第二阶段数据合并通称为联合数据d2。

进一步地，本发明步骤三的实现方法包括：

(3.1)对于第二阶段形状参数m2，给出取值区间范围[ma，mb]，选择均匀分布作为其先验分布，则m2的先验分布函数为：

(3.2)对于给定的任务工作时间τ，依据可靠性增长，认为第二阶段的可靠度大于第一阶段的可靠度，第二阶段产品在任务工作时间τ时的可靠度rτ的下限rl确定为：

选择均匀分布作为任务时间可靠度的先验分布，rτ的先验分布函数表达式如下所示：

(3.3)第二阶段产品的两参数weibull分布极大似然函数为：

其中，l(d2|m2，η2)表示联合数据d2在m2和η2为未知变量下的极大似然函数；对于weibull寿命分布，f(ti)表示在t＝ti时的概率密度，r(t1a)和r(tb)分别表示在t＝t1a和t＝tb时的产品可靠度；

在(3.2)中已经确定任务时间可靠度的先验分布，故而需要将下式代入极大似然函数中，得到联合数据d2在m2和rτ为未知变量下的极大似然函数l(d2|m2，rτ)；

根据bayes定理得到联合数据条件下的两变量联合验后分布函数。

进一步地，本发明步骤四的实现方法包括：

(4.1)验后分布的核函数如下：

(4.2)采用m-h采样的分块式算法计算得到m2、rτ的估计值。

(4.3)进行可靠度的区间估计；

定义γ为置信水平，则通过解下式可以得到任务时间可靠度的置信下限rll；

在本发明步骤(4.2)中采用m-h采样的分块式算法计算得到m2、rτ的估计值，步骤如下：

首先，定义一个函数bivexp(m2，rτ)：

初始化t＝50000；

θ＝1；

定义两个数组m2和rτ，其中下标范围为0～(t-1)；

m2[0]为均匀分布区间[1，4]上的随机数；(其中[0]表示该数组中的第一个数据，索引下标为0)

rτ[0]为均匀分布区间[0.9，1]上的随机数；

设置t＝0，同时每次迭代t一步递增，直至t＝t，进行以下迭代过程：

m2_star为均匀分布区间[1，4]上的随机数；

rτ_star为均匀分布区间[0.9，1]上的随机数；

定义：

定义：u为均匀分布区间(0，1)上的随机数；

如果u≤α：

m2[t]＝m2_star；rτ[t]＝r_sta；

否则：

m2[t]＝m2[t-1]；rτ[t]＝rτ[t-1]；

分别得到两个1*50000的数组m2和rτ；对50000个数据求和取平均值即为所求变量(m2，rτ)的估计值。

另一方面，本发明提供一种基于多阶段试验数据折合的航天产品可靠性评估系统，包括：

两阶段可靠性试验寿命数据采集模块，用于获取航天产品第一阶段以及第二阶段的产品试验寿命数据，并予以保存。

折合模块，用于对第一阶段的产品试验寿命数据进行折合。

确定模块，用于确定第二阶段weibull寿命分布参数的验前分布，并得到联合分布数据下的bayes验后分布核函数。

可靠度估计模块，采用基于bayes的m-h抽样方法进行可靠度的区间估计。

本发明还提供一种计算机设备，包括存储器和处理器，所述存储器存储有计算机程序，所述处理器执行所述计算机程序时实现上述任一种基于多阶段试验数据折合的航天产品可靠性评估方法的步骤。

本发明还提供一种计算机可读存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现上述任一种基于多阶段试验数据折合的航天产品可靠性评估方法的步骤。

与现有技术相比，本发明提出了基于weibull寿命分布的阶段间折合因子计算方法，给出了完整的bayes可靠性评估过程，并给出任务时间可靠度的区间估计。该方法不受样本量与失效数据类型的限制，将当前阶段无失效数据联合先前阶段的折合数据进行可靠性评估，提高了评估的精确度。本发明实施例方法思考成熟、可行性与精确度良好，为多阶段试验航天产品可靠性评估的研究打下基础。

附图说明

图1是本发明一实施例中的流程图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚明白，下面将以附图及详细叙述清楚说明本发明所揭示内容的精神，任何所属技术领域技术人员在了解本发明内容的实施例后，当可由本发明内容所教示的技术，加以改变及修饰，其并不脱离本发明内容的精神与范围。本发明的示意性实施例及其说明用于解释本发明，但并不作为对本发明的限定。

实施例1：

参照图1，提供一种基于多阶段试验数据折合的航天产品可靠性评估方法，包括以下步骤：

步骤一，进行航天产品的两阶段可靠性试验，获取航天产品第一阶段以及第二阶段的产品试验寿命数据：

航天产品第一阶段可靠性试验中：有n1个样本投入试验，在截尾时间ta时出现r个失效，其中失效数据为t1，t2，...，tr，无失效数据为tr+1＝tr+2＝…＝tn1＝ta，此阶段的weibull寿命分布参数为(m1，η1)。

航天产品第二阶段可靠性试验中：有n2个样本投入试验，在截尾时间tb时无失效数据产生，此阶段的weibull寿命分布参数为(m2，η2)。

步骤二，对第一阶段的产品试验寿命数据进行折合。

(2.1)对于第一阶段的产品试验寿命数据，已知其失效数据为t1，t2，...，tr1，并给定形状参数m1的值，故而可以使用极大似然估计(下文简称mle)方法得到尺度参数(特征寿命)η1的估计值。

两变量的mle极大似然函数为：

式中，d1表示第一阶段的产品试验寿命数据；l(d1|m1，η1)表示第一阶段的产品试验寿命数据的极大似然函数，其中包括两个变量m1和η1。

根据mle方法的原理，可以进一步得到η1的估计值。

完成对产品第一阶段可靠性试验的可靠性评估，其寿命分布函数为已知。

(2.2)计算第一阶段产品的失效概率；

其中：p1j表示在时间t＝tj(j＝1，...，n1)时第一阶段产品的失效概率。

(2.3)计算第二阶段产品的失效概率；

其中p2表示第二阶段产品的失效概率。

(2.4)进行数据折合，得到联合数据；

t1→t1j＝tjcj

其中cj(j＝1，...，n1)表示t＝tj时第一阶段数据向第二阶段数据折合的折合因子；第一阶段折合后的数据为t11，t12，...，t1r，t1a，...，t1a，与第二阶段数据合并通称为联合数据d2。

步骤三，确定第二阶段分布参数的验前分布，并得到联合分布数据下的bayes验后分布核函数。

经过步骤二，第一阶段数据已被折合为与第二阶段数据同寿命分布母体的数据。第二阶段产品的可靠性评估使用bayes方法进行。

(3.1)对于第二阶段形状参数m2，给出取值区间范围[ma，mb]，选择均匀分布作为其先验分布，则m2的先验分布函数为：

(3.2)对于给定的任务工作时间τ，依据可靠性增长，认为第二阶段的可靠度大于第一阶段的可靠度，第二阶段产品在任务工作时间τ时的可靠度rτ的下限rl确定为：

同样的，没有更多关于可靠度rτ的先验信息。因此，选择均匀分布作为任务时间可靠度的先验分布，rτ的先验分布函数表达式如下所示：

(3.3)第二阶段产品的两参数weibull分布极大似然函数为：

其中，l(d2|m2，η2)表示联合数据d2在m2和η2为未知变量下的极大似然函数；对于weibull寿命分布，f(ti)表示在t＝ti时的概率密度，r(t1a)和r(tb)分别表示在t＝t1a和t＝tb时的产品可靠度；

在(3.2)中已经确定任务时间可靠度的先验分布，故而需要将下式代入极大似然函数中，得到联合数据d2在m2和rτ为未知变量下的极大似然函数l(d2|m2，rτ)；

根据bayes定理得到联合数据条件下的两变量联合验后分布函数。

步骤四，采用基于bayes的m-h抽样方法进行参数估计。

由于联合数据条件下的两变量联合验后分布函数中包含二重复杂积分以及多元求解，传统的计算手段已经不再适用。这里采用抽样方法进行多元求解，并利用积分限法进行可靠度区间估计。

(4.1)验后分布的核函数如下：

(4.2)采用m-h算法计算得到m2、rτ的估计值。

对于m-h采样的分块式(block-wise)算法，在选择已知分布时，需要选择与目标分布具有相同维度的分布，即在上一步中获得的核函数。具体的，步骤如下：

首先，定义一个函数bivexp(m2，rτ)：

初始化t＝50000；

θ＝1；

定义两个数组m2和rτ，其中下标范围为0～(t-1)；

m2[0]为均匀分布区间[1，4]上的随机数，其中[0]表示该数组中的第一个数据，索引下标为0；

rτ[0]为均匀分布区间[0.9，1]上的随机数；

设置t＝0，同时每次迭代t一步递增，直至t＝t，进行以下迭代过程：

m2_star为均匀分布区间[1，4]上的随机数；

rτ_star为均匀分布区间[0.9，1]上的随机数；

定义：

定义：u为均匀分布区间(0，1)上的随机数；

如果u≤α：

m2[t]＝m2_star；rτ[t]＝rτ_star；

否则：

m2[t]＝m2[t-1]；rτ[t]＝rτ[t-1]；

分别得到两个1*50000的数组m2和rτ；对50000个数据求和取平均值即为所求变量(m2，rτ)的估计值。

(4.3)进行可靠度的区间估计。

定义γ为置信水平，则通过解下式可以得到任务时间可靠度的置信下限rll。

实施例2

在本实施例中根据某单位工程设定火箭发动机轴承在第一阶段中的weibull分布参数为m1＝1.5，η0＝2000，并有10个产品投入试验(n1＝n2＝10)，试验截尾时间ta＝tb＝500s，额定任务工作时间τ＝500s。依据此参数进行仿真得到试验过程中出现一个失效数据t1＝316s的试验数据。

采用本发明提供的基于多阶段试验数据折合的航天产品可靠性评估方法，对其进行可靠度估计，结果如下：

第一阶段可靠性评估结果为：m1＝1.5，η1＝2243.1。

折合因子计算结果为：c1＝0.8431，c2＝…＝cn＝1.6357。

折合后试验数据为：t11＝266.39，t12＝…＝t1n＝t1a＝817.84。

第二阶段先验分布结果为：

m-h抽样方法得到参数估计值为：m2＝1.99，rτ＝0.96。

不同置信水平γ得到可靠度的置信下限结果如下：

当γ＝0.75时，rll＝0.9380；

当γ＝0.90时，rll＝0.9190；

当γ＝0.95时，rll＝0.9108。

以上所述实施例仅表达了本申请的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本申请构思的前提下，还可以做出若干变形和改进，这些都属于本申请的保护范围。因此，本申请专利的保护范围应以所附权利要求为准。