CN106125713A - 一种区间删失情况下可靠性增长的评估与预测方法 - Google Patents

Abstract

一种区间删失情况下可靠性增长的评估与预测方法，构造一个服从幂律分布的次序统计量模型；根据设备的先验知识确定可靠性增长过程参数的初值，基于服从幂律分布的次序统计模型，随机地插补出缺失的故障时间；在缺失的故障时间填补完全的条件下，采用期望最大方法，迭代地更新可靠性增长过程的参数；评估当前时刻的平均无故障时间、预测未来故障的发生时间。该方法能够融合已经观测到的数据迭代地估计出失效时间分布并据此随机插补出缺失的故障时间，然后在完整的数据集上实现了系统可靠性增长规律的确定，实现了当前故障间隔时间的评估和未来故障发生时间的预测，弥补了统计分析上的不便，并能在非高斯的假设下有效地减小数据缺失造成的偏差。

Description

一种区间删失情况下可靠性增长的评估与预测方法

技术领域

本发明属于可靠性工程领域，具体涉及一种区间删失情况下可靠性增长的评估与预测方法。

背景技术

一个复杂的设备或系统，在投入正式使用的初期通常会不可避免地不断暴露出各种缺陷和故障。这些缺陷通常与该产品的设计方案、制造工艺、操作方法、维护技术以及管理水平等息息相关。针对这些软缺陷，采用“运行–暴露–改进–运行”的策略及时采取改进措施，能不断提高设备的可靠性，这就是设备的可靠性增长过程。在这一过程中，必须采用“可靠性增长评估方法”判断设备的可靠性是否增长并定量地评估及预测其增长的速度，以适时调整各类改进措施，以期能有效地控制可靠性增长过程的时间和费用。

工程实际中，已经逐渐发展出各种映射模型以解决复杂系统可靠性增长的评估问题。在现有的大量方法中，以Crow-AMSAA(army materiel systems analysis activity)模型的应用最为广泛。该模型本质上是将可靠性增长曲线用失效率为非时齐泊松过程模型来描述，是基于时间函数可靠性增长评估模型。在利用这些模型对系统可靠性进行动态地评定时，必须要有系统历史故障时间的精确记录。但由于测试仪器故障或人为疏漏等多种原因，故障时间记录不完整，尤其是区间数据删失的情况在可靠性增长试验阶段时常发生。这些缺失数据对可靠性的评估影响很大，处理不当极可能使各类可靠性指标的评估受到不同程度的歪曲。

通常，解决上述问题主要采用全信息极大似然的方式，即是直接构造数据缺失情况下模型中未知参数或者可靠性指标的条件概率分布，通过极大化对数似然函数直接对各类可靠性指标做出统计推断。然而这种方法存在的局限在于：在数据缺失的条件下，诸多统计推断方法会变得非常繁琐，且需在正态假设下下才可能得到无偏的估计。这就对基于现场数据的分析造成了困难，同时也为数据缺失条件下的可靠性增长评估提出挑战。

发明内容

针对上述困难和挑战，本发明的目的在于提供一种区间删失情况下可靠性增长的评估与预测方法，该方法能够融合已经观测到的数据迭代地估计出失效时间分布并据此随机插补出缺失的故障时间，然后在完整的数据集上实现了系统可靠性增长规律的确定，并运用广义Gamma分布实现了当前故障间隔时间的评估和未来故障发生时间的预测。这种方法一方面弥补了统计分析上的不便，并能在非高斯的假设下有效地减小数据缺失造成的偏差。

为达到以上目的，本发明的技术方案为：

一种区间删失情况下可靠性增长的评估与预测方法，包括以下步骤：

1)首先，假设可靠性增长过程服从一个非时齐的泊松过程，采用条件概率定理融合先验的可靠性增长趋势和已观测到的区间累计失效数，推断出缺失数据集的联合条件概率密度函数；

2)然后，基于该联合概率密度，构造一个服从幂律分布的次序统计量模型；

根据设备的先验知识确定可靠性增长过程参数的初值，基于服从幂律分布的次序统计模型，随机地插补出缺失的故障时间；

3)最后，引入蒙特卡罗期望最大化的推断架构；

在缺失的故障时间填补完全的条件下，采用期望最大方法，迭代地更新可靠性增长过程的参数；

运用广义Gamma分布评估当前时刻的平均无故障时间、预测未来故障的发生时间。

本发明进一步的改进在于，步骤2中具体步骤如下：

首先，设随机变量X服从0到1上的均匀分布，记为X～U[0,1]，运用逆变换方法，得到一个服从幂律分布的随机变量t:

$t=(t_{i+1}^b-t_i^b)x+t_i^b,x~U\[0,1\]$

其中，t_i+1和t_i是已知区间的两个端点，t_i的取值范围为任意非负实数且应满足t_i＜t_i+1；

则随机变量t服从幂律分布，且具体分布律p(t)如下：

$p\left(t\right)=\frac{{\mathrm{abt}}^{b-1}}{({\mathrm{at}}_{i+1}^{b-1}-{\mathrm{at}}_i^{b-1})}$

完成服从幂律分布的次序统计量模型的构造；

其次，从这一个服从幂律分布的总体中独立抽取n_i个样本，按升序排列，得到一组次序统计量；该次序统计量能够随机模拟某一区间缺失的故障时间；

然后，对每个区间重复以上过程，得到多组次序统计量，依据区间顺序对多组次序统计量进行升序排列，得到一个新的次序统计量，该次序统计量能够插补整个测试过程中的缺失故障时间。

本发明进一步的改进在于，步骤3中在缺失的故障时间填补完全的条件下，采用期望最大方法，迭代地更新可靠性增长过程的参数的具体过程如下：

当缺失的故障时间填补完全时，对参数a,b的对数极大似然函数在缺失故障时间概率分布上求期望，如下式：

$\begin{array}{c}Q(a,b|a^{\left(i\right)},b^{\left(i\right)},t_{\left(1\right)},t_{\left(2\right)},\mathrm{...},t_{\left(n\right)})\\ =\underset{t_{\left(1\right)},\mathrm{...},t_{\left(n\right)}}{\∫}\log \[p(a,b|a^{\left(i\right)},b^{\left(i\right)},t_{\left(1\right)},t_{\left(2\right)},\mathrm{...},t_{\left(n\right)})\]p\left(Z\right|a^{\left(i\right)},b^{\left(i\right)},Y;N\left(t_{i+1}\right)-N\left(t_i\right)=n_i)d{t}_{\left(1\right)}\mathrm{...}d{t}_{\left(n\right)}\end{array}$

运用蒙特卡罗法求上述期望，上式简化如下：

$Q(a,b|a^{\left(i\right)},b^{\left(i\right)},t_{\left(1\right)},t_{\left(2\right)},\mathrm{...},t_{\left(n\right)})=\frac{1}{m}\underset{j=1}{\overset{m}{\Σ}}\left(n\log \right(a^{\left(i\right)}{b}^{\left(i\right)})+(-a^{\left(i\right)}{t}_n^{b^{\left(i\right)}})+\underset{i=1}{\overset{n}{\Σ}}(\left(b^{\left(i\right)}-1\right)\log t_i\left)\right)$

然后，上式分别相对a⁽ⁱ⁾,b⁽ⁱ⁾求偏导并令其偏导等于0，得到a⁽ⁱ⁾,b⁽ⁱ⁾的新值，将该值作为次序统计量中参数的初值，反复迭代直至满足停止条件。

本发明进一步的改进在于，a⁽ⁱ⁾,b⁽ⁱ⁾的新值计算公式如下：

$\frac{\∂Q}{\∂a}=\frac{n}{a}-\frac{1}{m}\underset{j=1}{\overset{m}{\Σ}}{t}_{nj}^b=0$

$\frac{\∂Q}{\∂b}=\frac{n}{b}-{\mathrm{at}}_n^b\log t_n+\frac{1}{m}\underset{j=1}{\overset{m}{\Σ}}\underset{i=1}{\overset{n}{\Σ}}\log t_{ij}=0.$

本发明进一步的改进在于，步骤3中运用广义Gamma分布评估当前时刻的平均无故障时间(Mean Time betweenFailure,MTBF)、预测未来故障的发生时间的具体过程如下：

基于上述确定的可靠性增长过程关键参数，得到失效时间的联合概率密度函数，对其进行依次积分，得故障率为λ(t)＝abt^b-1的非时齐泊松过程的第n次失效时间服从广义Gamma分布，记为：

$t_n~G\mathrm{\Γ}\left(t\right|n,a,b)=\left\{\begin{array}{c}\frac{{\mathrm{ba}}^n}{\mathrm{\Γ}\left(n\right)}{t}^{nb-1}\exp (-a{t}^b),t>0\\ 0,t\≤0\end{array}\right.$

由此，评估系统平均故障间隔时间M[t_n]和预测未来第k次故障发生时间如下：

$M\[t_n\]={\left(t_n\right)}^{1-\stackrel{\‾}{b}}/\left(\stackrel{\‾}{a}\stackrel{\‾}{b}\right)$

${\stackrel{\‾}{t}}_{n+k}=E\left(t_{n+k}\right)=a^{(-b^{-1})}\frac{\mathrm{\Γ}(n+k+b^{-1})}{\mathrm{\Γ}(n+k)}.$

与现有技术比较，本发明的优势为：

1)本发明采用非时齐泊松过程对可靠性增长过程进行建模，基于条件概率定理融合区间故障计数信息，能够得到缺失故障时间的分布规律，由此所进行的数据插补能有效地减小数据缺失造成的统计偏差；

2)本发明采用蒙特卡罗期望最大方法，借助一个服从幂律分布的次序统计量模型，改进了缺失故障时间的插补过程，该方法能在增长模型参数未知的条件下，简便地实现缺失数据的插补；

3)本发明整个分析与评估过程实现了整体信息和局部信息的融合，迭代地推断出可靠性增长过程的参数，为数据缺失条件下可靠性评估提供一种有效的新技术。

4)本发明的可靠性评估方法简单可靠，便于在工程实践中使用。

5)该方法能够融合已经观测到的数据迭代地估计出失效时间分布并据此随机插补出缺失的故障时间，然后在完整的数据集上实现了系统可靠性增长规律的确定，并运用广义Gamma分布实现了当前故障间隔时间的评估和未来故障发生时间的预测。这种方法一方面弥补了统计分析上的不便；另一方面，而且更为重要的是，它能在非高斯的假设下有效地减小数据缺失造成的偏差。

附图说明

图1为管道压缩机辅助系统故障停机次数的统计直方图；

图2为管道压缩机辅助系统累计失效数-累计试验时间图；

图3为本发明针对的故障时间缺失情况示意图；

图4为发明提出可靠性评估预测方法的流程图；

图5为本发明第50次故障到第91次故障时间预测结果图。

图6为本发明第69次故障到第91次故障时间预测结果图。

图7为本发明第73次故障到第91次故障时间预测结果图。

具体实施方式

下面结合附图对本发明做详细描述。

本发明提供了一种区间删失情况下可靠性增长的评估与预测方法，包括以下步骤：

1)收集整个可靠性增长测试过程中的数据，将缺失的设备故障时间序列构成一组取值分布未知的随机向量。

首先，假设可靠性增长过程服从一个以λ(t)＝abt^b-1为强度函数非时齐的泊松过程，采用条件概率定理融合先验的可靠性增长趋势和已观测到的区间累计失效数，推断出缺失数据集的联合条件概率密度函数；

推断出缺失数据集的联合条件概率密度函数的具体过程如下：

(1)假设任一区间上的失效累计数服从一个以λ(t)＝abt^b-1为强度函数的非时齐的泊松过程，其中t为非负实数且a,b是强度函数中的参数。即区间[t_i,t_i+1)上发生n次失效的概率P{N(t_i+1)-N(t_i)＝n}如下：

$P\{N\left(t_{i+1}\right)-N\left(t_i\right)=n\}=\frac{\exp (-\underset{t_i}{\overset{t_{i+1}}{\∫}}{\mathrm{abs}}^{b-1}ds){\[\underset{t_i}{\overset{t_{i+1}}{\∫}}{\mathrm{abs}}^{b-1}ds\]}^n}{n!}$

其中，t_i+1和t_i是已知区间的两个端点，t_i的取值范围为任意非负实数且应满足t_i＜t_i+1；

(2)若在区间[t_i,t_i+1)上，假设缺失的故障时间分别为[t_(i,1),t_(i,2),…,t_(i,n)],基于非时齐的泊松过程的性质，计算故障时间的联合概率密度如下：

$p_{t_{(i,1)},...,t_{(i,n)}}(t_{(i,1)},...,t_{(i,n)})=\underset{k=1}{\overset{n}{\Π}}\mathrm{\λ}\left(t_{\left(k\right)}\right)\exp (-({\mathrm{at}}_{i+1}^b-{\mathrm{at}}_i^b\left)\right)$

(3)若在区间[t_i,t_i+1)上，基于区间中累计失效数为N(t_i+1)-N(t_i)＝n_i，重新修正缺失故障时间的联合概率密度如下：

$p_{t_{(i,1)},...,t_{(i,n_i)}|N(t_i+1)-N\left(t_i\right)=n_i}(t_{(i,1)},...,t_{(i,n)}|n(t_i+1)-n\left(t_i\right)=n_i)=n_i!\underset{k=1}{\overset{n_i}{\Π}}(\frac{{\mathrm{abt}}_{(i,k)}^{b-1}}{({\mathrm{at}}_{i+1}^b-{\mathrm{at}}_i^b)})$

该式即为融合了非时齐泊松过程和区间累计故障数的缺失数据集的联合条件概率密度。

2)然后，基于该联合概率密度，构造一个服从幂律分布的次序统计量模型，服从幂律分布的次序统计量模型能够将多维且相关的缺失数据集转化为多个一维且相互独立的数据，由此可以简化复杂缺失数据集的插补过程，实现高效地插补。具体过程如下：

(1)设随机变量X服从0到1上的均匀分布，记为X～U[0,1]，运用逆变换方法，经过如下变换可得到一个服从幂律分布的随机变量t:

$t=(t_{i+1}^b-t_i^b)x+t_i^b,x~U\[0,1\]$

其中，t_i+1和t_i是已知区间的两个端点。

其中，服从幂律分布的随机变量t服从幂律分布，且具体分布律p(t)如下：

$p\left(t\right)=\frac{{\mathrm{abt}}^{b-1}}{({\mathrm{at}}_{i+1}^{b-1}-{\mathrm{at}}_i^{b-1})}$

完成服从幂律分布的次序统计量模型的构造。

(2)从这一个服从幂律分布的总体中独立抽取n_i个样本，按升序排列，得到一组次序统计量。

次序统计量的分布定理：设T₁,...,T_n来自某总体的一个样本，总体的分布函数为F(t)，概率密度函数为p(t)。n个次序统计量t₍₁₎,...,t_(n)的联合概率密度为：

g(t₍₁₎,...t_(n))＝n！p(t₍₁₎)p(t₍₂₎)...p(t_(n))，x₁＜x₂＜...＜x_n

该定理可以保证，上述顺序统计量的概率密度与缺失故障时间的联合概率密度相同。

该次序统计量即可以实现某一区间缺失的故障时间的插补。

(3)对每个区间重复以上(1)、(2)过程，依据区间顺序对次序统计量进行排列，得到一个新的次序统计量，该含序统计量即可实现整个测试过程中的故障时间的插补。

3)最后，引入蒙特卡罗期望最大化的推断架构。

根据设备的先验知识确定可靠性增长过程参数的初值，基于步骤2)中的次序统计模型，随机地插补出缺失的故障时间；在缺失故障时间填补完全的条件下，采用期望最大方法，迭代地更新可靠性过程的参数；具体过程如下：

(1)将未知的故障时间作为有添条件，据此计算可靠性增长过程中关于未知参数a和b的对数极大似然函数如下：

$log\[p(a,b|Y,Z)\]=nlog\left(ab\right)+(-{\mathrm{at}}_n^b)+\underset{i=1}{\overset{n}{\Σ}}\left(\right(b-1\left)\log t_i\right)$

其中为Y已知条件，表示各区间的累计失效数，即：

Y:{N(t₁)-N(t₀)＝n₁,...,N(t_i)-N(t_i-1)＝n_i}

其中Z为缺失的数据，具体表示各个未知的具体失效时间，即：

Z:{t₍₁₎,t₍₂₎,...,t_(n)}

(2)将上述对数极大似然函数相对于缺失数据的概率密度函数求积分，求出关于缺失数据的条件期望的表达式，该表达式实际上是关于参量的函数，如下：

$\begin{array}{c}Q\left(a,b|a^{\left(i\right)},b^{\left(i\right)}\right)\\ =\underset{Z}{\∫}\log \[p(a,b|a^{\left(i\right)},b^{\left(i\right)},Z)\]p\left(Z\right|a^{\left(i\right)},b^{\left(i\right)},Y)dZ\end{array}$

(3)利用缺失故障时间的模拟值，基于蒙特卡罗法对上式求期望，运用蒙特卡罗法求上述期望，上式可以简化如下：

$Q(a,b|a^{\left(i\right)},b^{\left(i\right)},t_{\left(1\right)},t_{\left(2\right)},\mathrm{...},t_{\left(n\right)})=\frac{1}{m}\underset{j=1}{\overset{m}{\Σ}}\left(n\log \right(a^{\left(i\right)}{b}^{\left(i\right)})+(-a^{\left(i\right)}{t}_n^{b^{\left(i\right)}})+\underset{i=1}{\overset{n}{\Σ}}(\left(b^{\left(i\right)}-1\right)\log t_i\left)\right)$

然后，上式分别相对a⁽ⁱ⁾,b⁽ⁱ⁾求偏导并令其偏导等于0，得到a⁽ⁱ⁾,b⁽ⁱ⁾的新值，若新值与旧值之差大于设定阈值时，返回至次序统计量赋初值阶段并将新值设定为步骤2)所需初值，若新值与旧值之差小于设定阈值时，结束循环并将该值输出进行后续的可靠性的评估与预测。

新值计算公式如下：

$\frac{\∂Q}{\∂a}=\frac{n}{a}-\frac{1}{m}\underset{j=1}{\overset{m}{\Σ}}{t}_{nj}^b=0$

$\frac{\∂Q}{\∂b}=\frac{n}{b}-{\mathrm{at}}_n^b\log t_n+\frac{1}{m}\underset{j=1}{\overset{m}{\Σ}}\underset{i=1}{\overset{n}{\Σ}}\log t_{ij}=0$

由此可得到参数a和b的更新值和缺失故障时间的随机模拟值。将这些更新值作为初值再次代入步骤2)和3)中进行迭代，直至满足停止条件。

运用广义Gamma分布评估当前时刻的平均无故障时间(Mean TimebetweenFailure,MTBF)、预测未来故障的发生时间。具体过程如下：

基于上述确定的可靠性增长过程关键参数，可得到失效时间的联合概率密度函数，对其进行依次积分，可得故障率为λ(t)＝abt^b-1的非时齐泊松过程的第n次失效时间服从广义Gamma分布，记为：

$t_n~G\mathrm{\Γ}\left(t\right|n,a,b)=\left\{\begin{array}{c}\frac{{\mathrm{ba}}^n}{\mathrm{\Γ}\left(n\right)}{t}^{nb-1}\exp (-a{t}^b),t>0\\ 0,t\≤0\end{array}\right.$

由此,可以评估系统平均故障间隔时间M[t_n]和预测未来第k次故障发生时间如下：

$M\[t_n\]={\left(t_n\right)}^{1-\stackrel{\‾}{b}}/\left(\stackrel{\‾}{a}\stackrel{\‾}{b}\right)$

${\stackrel{\‾}{t}}_{n+k}=E\left(t_{n+k}\right)=a^{(-b^{-1})}\frac{\mathrm{\Γ}(n+k+b^{-1})}{\mathrm{\Γ}(n+k)}$

下面结合附图对本发明的内容作进一步详细说明：

实施例：

该实施案例结合某管道压缩机辅助系统的停机数据验证了该发明的有效性。

某管道压缩机辅助系统从2012年至2014年间共有91次故障停机记录，具体数据如下表所示：

表1 管道压缩机辅助系统2012年至2014年故障停机情况

图1是管道压缩机辅助系统故障停机次数的统计直方图，由此可发现压缩机辅助系统故障次数具有波动降低的趋势。

图2是管道压缩机辅助系统累计失效数-累计试验时间图，该曲线向上凸，此时相邻故障间隔时间增大，产品可靠性增长。

图3是故障时间数据的缺失模式示意图。它也是现场数据采集中经常发生的数据缺失模式。该模式下，整个测试过程被分隔为若干个区间，区间内的累计运行时间和相应的累计故障停机数都被精确的记录下，但每次故障发生的具体时间缺失。分析发现，表1中的数据缺失模式即属于此类。

图4是本发明提出的一种区间删失情况下可靠性增长的评估与预测方法的处理流程图，下面依照该流程处理管道压缩机辅助系统故障停机。

某管道压缩机辅助系统具有明显的可靠性增长趋势，但由于缺失具体的失效时间数据，可靠性增长的精确建模比较困难。鉴于此，可靠性增长评估与预测的主要框架包括：建立缺失故障时间的概率分布，模型的参数推断，可靠性评估与预测等三个部分。

(1)首先设可靠性增长过程服从故障率为λ(t)＝abt^b-1的非时齐泊松过程，其中a,b是未知参数。整体的增长规律融合11个区间上的故障数，分别建立每个区间上缺失故障时间的概率密度函数如下，该式即为缺失故障时间的概率分布：

$p(t_{(i,1)},\mathrm{...},t_{(i,n_i)}|n\left(t_{i+1}\right)-n\left(t_i\right)=n_i)=\frac{\left(\underset{m=1}{\overset{n_i}{\Π}}\mathrm{\λ}\right(t_{(i,m)}\left)\exp \right(-\left({\mathrm{at}}_{i+1}^b-{\mathrm{at}}_i^b\right)\left)\right)}{\frac{\exp (-\underset{t_i}{\overset{t_{i+1}}{\∫}}\mathrm{\λ}(s\left)dx\right){\[\underset{t_i}{\overset{t_{i+1}}{\∫}}\mathrm{\λ}\left(s\right)ds\]}^n}{n_i!}}=n_i!\underset{m=1}{\overset{n_i}{\Π}}\left(\frac{{\mathrm{abt}}_{(i,m)}^{b-1}}{({\mathrm{at}}_{i+1}^b-{\mathrm{at}}_i^b)}\right)$

其中，i取值为0到10，分别代表11个区间；n_i代表第i个区间上的故障累计数；t_(i,m)代表第i个区间上第m次故障的发生时间。

(2)然后，对第i区间，从分布律为的幂律分布总体中独立抽取n_i个独立同分布的随机变量，并依升序排列成一组次序统计量。对第1个到第11个区间进行相同的操作，再将这些次序统计量依区间顺序排列，即可得到一组新的次序统计量该组次序统计量即可模拟缺失故障时间。

(3)再次，对次序统计量中的未知参数给定一组初值a₀,b₀，一般而言，可靠性发生增长时b₀小于1，这里将b₀设置为0.5。根据参数已确定的次序统计量模型产生多组缺失数据的插补值在此基础上利用蒙特卡洛方法对参数a,b的对数极大似然函数在缺失故障时间概率分布上求期望，即是处理以下含参数的求和式：

$Q(a,b|a^{\left(i\right)},b^{\left(i\right)},t_{\left(1\right)},t_{\left(2\right)},\mathrm{...},t_{\left(n\right)})=\frac{1}{m}\underset{j=1}{\overset{m}{\Σ}}\left(n\log \right(a^{\left(i\right)}{b}^{\left(i\right)})+(-a^{\left(i\right)}{t}_n^{b^{\left(i\right)}})+\underset{i=1}{\overset{n}{\Σ}}(\left(b^{\left(i\right)}-1\right)\log t_i\left)\right)$

对上式分别相对a⁽ⁱ⁾,b⁽ⁱ⁾求偏导并令其等于0，可得到a⁽ⁱ⁾,b⁽ⁱ⁾的新值，将该值作为次序统计量中的参数初值，反复迭代直至满足停止条件。新值计算公式如下：

$\frac{\∂Q}{\∂a}=\frac{n}{a}-\frac{1}{m}\underset{j=1}{\overset{m}{\Σ}}{t}_{nj}^b=0$

$\frac{\∂Q}{\∂b}=\frac{n}{b}-{\mathrm{at}}_n^b\log t_n+\frac{1}{m}\underset{j=1}{\overset{m}{\Σ}}\underset{i=1}{\overset{n}{\Σ}}\log t_{ij}=0$

$\⇒\frac{\∂Q}{\∂b}=\frac{n}{b}-\frac{n}{\underset{j=1}{\overset{m}{\Σ}}{t}_{n_j}^b}\×\underset{i=1}{\overset{n}{\Σ}}\underset{j=1}{\overset{m}{\Σ}}\left(t_{ij}^b\ln t_{ij}\right)+\frac{1}{m}\underset{i=1}{\overset{n}{\Σ}}\underset{j=1}{\overset{m}{\Σ}}\left(\ln t_{ij}\right)=0$

由此，计算出a＝0.70；b＝0.38。

(4)最后，由此计算得到的a，b的值,可以估计当前系统平均故障间隔时间和未来第k次故障发生时间如下：

$M\[t_n\]={\left(t_n\right)}^{1-\stackrel{\‾}{b}}/\left(\stackrel{\‾}{a}\stackrel{\‾}{b}\right)$

${\stackrel{\‾}{t}}_{n+k}=E\left(t_{n+k}\right)=a^{(-b^{-1})}\frac{\mathrm{\Γ}(n+k+b^{-1})}{\mathrm{\Γ}(n+k)}$

图5记录了从第49次故障发生开始预测其后42次故障发生的时间；图6记录了从第69次故障发生开始预测其后22次故障发生的时间；图7记录了从第73次故障发生开始预测其后18次故障发生时间。结果表明，从第73次失效开始，本发明提出方法已经能够很好的预测了未来18次失效的发生时间，同时对比其他方法本发明提出方法的可靠性评估结果更准确，故障时间预测精度更高。

本发明提出的数据不完备情况下可靠性增长的评估方法，能够方便地融合全局可靠性增长的趋势以及局部可观察数据，通过随机迭代地方式实现了缺失时间的分布规律的估计和缺失时间数据的随机模拟并最终实现系统的可靠性准确评估和预测。

针对仪器故障、人员疏漏和不完善的维修策略而导致的区间数据删失问题，首先利用可靠性增长整体趋势服从非时齐泊松过程这一假设，融合局部可观察数据推导出有序缺失数据集的联合概率分布。其次，基于上述联合概率分布，提出一个新的次序统计量模型，可将多维且相关的缺失数据集等价地转化为多个一维且相互独立数据，以实现高效插补。最后，结合该次序统计量模型和蒙塔卡罗-期望最大方法，迭代地进行多重插补和未知参数推断，实现系统可靠性增长规律的确定和未来故障发生时间的预测。该方法计算简便、效果显著，为区间数据删失情况下可靠性增长的评估与预测提供了一种有效的方法。

Claims (5)

1.一种区间删失情况下可靠性增长的评估与预测方法，其特征在于，包括以下步骤：

1)首先，假设可靠性增长过程服从一个非时齐的泊松过程，采用条件概率定理融合先验的可靠性增长趋势和已观测到的区间累计失效数，推断出缺失数据集的联合条件概率密度函数；

2)然后，基于该联合概率密度，构造一个服从幂律分布的次序统计量模型；

根据设备的先验知识确定可靠性增长过程参数的初值，基于服从幂律分布的次序统计模型，随机地插补出缺失的故障时间；

3)最后，引入蒙特卡罗期望最大化的推断架构；

在缺失的故障时间填补完全的条件下，采用期望最大方法，迭代地更新可靠性增长过程的参数；

运用广义Gamma分布评估当前时刻的平均无故障时间、预测未来故障的发生时间。

2.根据权利要求1所述的一种区间删失情况下可靠性增长的评估与预测方法，其特征在于，步骤2中具体步骤如下：

首先，设随机变量X服从0到1上的均匀分布，记为X～U[0,1]，运用逆变换方法，得到一个服从幂律分布的随机变量t:

$t=(t_{i+1}^b-t_i^b)x+t_i^b,x~U\[0,1\]$

其中，t_i+1和t_i是已知区间的两个端点；t_i的取值范围为任意非负实数且应满足t_i<t_i+1；

其中，服从幂律分布的随机变量t服从幂律分布，且具体分布律p(t)如下：

$p\left(t\right)=\frac{{\mathrm{abt}}^{b-1}}{({\mathrm{at}}_{i+1}^{b-1}-{\mathrm{at}}_i^{b-1})}$

完成服从幂律分布的次序统计量模型的构造；

其次，从这一个服从幂律分布的总体中独立抽取n_i个样本，按升序排列，得到一组次序统计量；该次序统计量能够随机模拟某一区间缺失的故障时间；

然后，对每个区间重复以上过程，得到多组次序统计量，依据区间顺序对多组次序统计量进行升序排列，得到一个新的次序统计量，该次序统计量能够插补整个测试过程中的缺失故障时间。

3.根据权利要求1所述的一种区间删失情况下可靠性增长的评估与预测方法，其特征在于，步骤3中在缺失的故障时间填补完全的条件下，采用期望最大方法，迭代地更新可靠性增长过程的参数的具体过程如下：

当缺失的故障时间填补完全时，对参数a,b的对数极大似然函数在缺失故障时间概率分布上求期望，如下式：

$\begin{array}{c}Q(a,b|a^{\left(i\right)},b^{\left(i\right)},t_{\left(1\right)},t_{\left(2\right)},...,t_{\left(n\right)})\\ =\underset{t_{\left(1\right)},...,t_{\left(n\right)}}{\&Integral;}log\&lsqb;p(a,b|a^{\left(i\right)},b^{\left(i\right)},t_{\left(1\right)},t_{\left(2\right)},...,t_{\left(n\right)})\&rsqb;p\left(Z\right|a^{\left(i\right)},b^{\left(i\right)},Y;N\left(t_{i+1}\right)-N\left(t_i\right)=n_i){\mathrm{dt}}_{\left(1\right)}...{\mathrm{dt}}_{\left(n\right)}\end{array}$

运用蒙特卡罗法求上述期望，上式简化如下：

$Q(a,b|a^{\left(i\right)},b^{\left(i\right)},t_{\left(1\right)},t_{\left(2\right)},...,t_{\left(n\right)})=\frac{1}{m}\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}\left(n\log \right(a^{\left(i\right)}{b}^{\left(i\right)})+(-a^{\left(i\right)}{t}_n^{b^{\left(i\right)}})+\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}(\left(b^{\left(i\right)}-1\right)\log t_i\left)\right)$

然后，上式分别相对a⁽ⁱ⁾,b⁽ⁱ⁾求偏导并令其偏导等于0，得到a⁽ⁱ⁾,b⁽ⁱ⁾的新值，将该值作为次序统计量中参数的初值，反复迭代直至满足停止条件。

4.根据权利要求3所述的一种区间删失情况下可靠性增长的评估与预测方法，其特征在于，a⁽ⁱ⁾,b⁽ⁱ⁾的新值计算公式如下：

$\frac{\&part;Q}{\&part;a}=\frac{n}{a}-\frac{1}{m}\underset{j=1}{\overset{m}{\&Sigma;}}{t}_{nj}^b=0$

$\frac{\&part;Q}{\&part;b}=\frac{n}{b}-{\mathrm{at}}_n^b\log t_n+\frac{1}{m}\underset{j=1}{\overset{m}{\&Sigma;}}\underset{i=1}{\overset{n}{\&Sigma;}}log{t}_{ij}=0.$

5.根据权利要求1所述的一种区间删失情况下可靠性增长的评估与预测方法，其特征在于，步骤3中运用广义Gamma分布评估当前时刻的平均无故障时间(Mean TimebetweenFailure,MTBF)、预测未来故障的发生时间的具体过程如下：

基于上述确定的可靠性增长过程关键参数，得到失效时间的联合概率密度函数，对其进行依次积分，得故障率为λ(t)＝abt^b-1的非时齐泊松过程的第n次失效时间服从广义Gamma分布，记为：

$t_n~G\mathrm{\&Gamma;}\left(t\right|n,a,b)=\left\{\begin{array}{c}\frac{{\mathrm{ba}}^n}{\mathrm{\&Gamma;}\left(n\right)}{t}^{nb-1}\exp (-a{t}^b),t>0\\ 0,t\&le;0\end{array}\right.$

由此，评估系统平均故障间隔时间M[t_n]和预测未来第k次故障发生时间如下：

$M\&lsqb;t_n\&rsqb;={\left(t_n\right)}^{1-\stackrel{\&OverBar;}{b}}/\left(\stackrel{\&OverBar;}{a}\stackrel{\&OverBar;}{b}\right)$

${\stackrel{\&OverBar;}{t}}_{n+k}=E\left(t_{n+k}\right)=a^{(-b^{-1})}\frac{\mathrm{\&Gamma;}(n+k+b^{-1})}{\mathrm{\&Gamma;}(n+k)}.$