具体实施方式

下面结合具体实施例对本发明进行详细说明。以下实施例将有助于本领域的技术人员进一步理解本发明，但不以任何形式限制本发明。应当指出的是，对本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变化和改进。这些都属于本发明的保护范围。

如图1-图2所示，本发明提供了一种产品可靠性综合评估方法，从各阶段各类型数据的特点入手，对实际工程中存在的可靠性数据进行分类，给出了一种多阶段多类型数据的分类和预处理方法，结合Bayes及其相关模型，提出了一套融合多阶段多类型数据的产品可靠性综合评估流程。包括如下步骤：

步骤S1：将多阶段多类型数据进行分类与预处理，具体包括以下步骤：

步骤S1.1：多阶段多类型数据进行分类，数据类型如下：

产品现场类数据：产品实际使用或服役环境条件下的寿命型、成败型和退化型数据；

产品异环境类数据：产品在寿命周期各阶段开展的不同试验条件下不同于现场使用环境条件下的寿命型、成败型和退化型数据；各阶段包括可靠性试验、寿命试验、加速寿命试验、功能试验、环境试验、定期试验等阶段。

相似产品数据：又称为异总体数据，包括相似产品现场类数据和相似产品异环境类数据，相似产品现场类数据是指相似产品实际使用或服役环境条件下的寿命型、成败型和退化型数据，相似产品异环境类数据是指相似产品在寿命周期各阶段开展的不同试验条件下不同于现场使用环境条件下的寿命型、成败型和退化型数据；

部件级数据：产品各组成部件的现场类数据和异环境类数据，部件级现场类数据是指产品各组成部件实际使用或服役环境条件下的寿命型、成败型和退化型数据，部件级异环境类数据是指产品各组成部件在寿命周期各阶段开展的不同试验条件下不同于现场使用环境条件下的寿命型、成败型和退化型数据。

专家型数据：利用专家对产品可靠度的经验估计值来表示的数据。

步骤S1.2：对多阶段多类型数据进行预处理，为了便于数据的表达和预处理，对上述数据进一步归类，多阶段多类型数据可分为四类基本数据，即寿命型数据、成败型数据、退化型数据和专家型数据。四类基本数据类型的统一格式分别见表1至4：

表1寿命型数据的统一格式

表2成败型数据的统一格式

表3退化型数据的统一格式

表4专家型数据的统一格式

表1中“寿命类型”由寿命分布类型检验预处理模块确定，“故障时间”和“截尾时间”用于记录不同截尾方式(包括定时、定数、定总、随机截尾)对应的样品故障时间和截尾时间。表2中的“试验次数”和“失败次数”用于记录样品成败型数据的总试验次数和失败次数。表3中的“观测时刻”和“性能退化量”用于记录样品退化型数据的观测时间点和性能退化量。表4中的可靠度的“点估计”和“置信下限估计”(给定置信度)为专家对可靠度给出的主观估计值。

表1、表2和表3中的“相似因子”用于反映被评估产品与样品之间的相似程度，由相似产品数据预处理模块确定。当相似因子＝1时，表示样品技术状态与被评估产品完全相同；当0<相似因子<1时，表示样品的可靠性较被评估产品低；当相似因子>1时，表示样品的可靠性较被评估产品高。

表1、表2和表3中的“环境因子”是异环境数据与现场数据的环境折算因子或加速系数。当环境因子＝1时，表示该数据为现场数据；当环境编号>0且≠1时，则表示该数据为异环境数据。

进一步明确分类之后，分别针对寿命型数据、异环境类数据、部件级数据、退化型数据、相似产品数据进行预处理。

针对寿命型数据的预处理：利用寿命分布类型检验和参数评估方法进行预处理，得到寿命类型及其参数估计；寿命类型包括指数寿命型、威布尔寿命型、正态寿命型、对数正态型等。

针对异环境类数据的预处理：利用环境因子估算方法进行预处理，给出异环境类数据相对于现场数据的“环境因子”；异环境类数据包括产品异环境类数据、部件级异环境类数据、相似产品异环境类数据，环境因子估算方法是对不同分布类型的环境因子进行估算。

针对部件级数据的预处理：利用系统可靠性综合方法进行预处理，将部件级数据折合成产品数据；部件级数据包括部件级现场类数据和异环境数据，系统可靠性综合方法包括MML、SR、L-M、CMSR、Bayes法等。

针对退化型数据的预处理：利用退化建模方法进行预处理，将退化型数据转换成伪寿命型数据；退化型数据包括产品现场类退化型数据、部件级现场类退化型数据、部件级异环境类退化型数据，退化建模方法包括基于退化量分布、基于Wiener、Gamma等随机过程进行退化量分布。

针对相似产品数据的预处理：利用相似因子估算方法进行预处理，给出相似产品数据的“相似因子”；相似产品数据包括相似产品寿命型数据、相似产品成败型数据、相似产品异环境类寿命型数据和相似产品异环境类成败型数据，相似因子估算方法包括综合评判法、专家打分法等。

各类数据的处理流程如图3至图9所示。

步骤S2：将分类与预处理之后的多阶段多类型数据进行一致性检验；由于现场数据能够真实地反映产品的可靠性水平，因此在一致性检验时，以产品现场数据为一致性检验基准，将其他多阶段多类型数据与现场数据进行寿命分布类型和失效机理一致性检验。只有通过一致性检验的数据才能用于产品的可靠性综合评估。基本方法如下：

步骤S2.1：对于寿命型数据，根据寿命类型，判断多阶段多类型寿命数据(包括异环境寿命数据、相似产品寿命数据、伪寿命数据)与产品现场数据是否服从同一种寿命分布。

步骤S2.2：对失效机理进行一致性检验，对于异环境寿命数据(环境因子≠1)，还需要进行失效机理的一致性检验。根据寿命分布的类型及其失效机理不变条件，见表5，判断异环境寿命数据是否满足失效机理不变条件，即环境因子不变原则。只有满足失效机理不变条件，才能将异环境数据和现场数据进行融合。

表5常见寿命分布的环境因子与失效机理不变的条件

步骤S3：一致性检验后对多阶段多类型数据进行折算和综合。具体包括如下步骤：

步骤S3.1：对异环境数据的折算和综合：对于环境因子K_i0≠1的异环境数据，根据数据类型将异环境数据折算成现场数据，再将折算后的现场数据进行综合。基本方法如下：

①指数寿命型数据

保持故障数z不变，试验总时间τ乘以加速因子K_i0或其置信下限即

将所有异环境折算后的数据进行综合，其办法是：

(z′₀,τ′₀)＝(z₁,τ₁·K₁₀)+(z₂,τ₂·K₂₀)+…+(z_q,τ_q·K_q0) (2)

②其他寿命型数据(包括威布尔、正态、对数正态寿命型)

所有故障时间数据和截尾数据乘以加速因子K_i0或其置信下限即

其中：t_i为环境S_i下的试验数据，包括每个样本的故障时间或截尾时间。

对所有异环境折算后的数据进行综合，其办法是：

t′₀＝(t′₁₀,t′₂₀,…,t′_q0) (4)

③成败型数据

保持失败数f_i不变，试验次数n_i乘以加速因子K_i0

对所有异环境折算后的数据进行综合，其办法是：

(f′₀,n′₀)＝(f′₁,n′₁)+(f′₂,n′₂)+…+(f′_q,n′_q) (6)

步骤S3.2：对相似产品数据的折算和综合：对于相似因子F_i0≠1的相似产品数据，根据数据类型将相似产品数据折算成被评估产品数据，再将折算后的被评估产品数据进行综合。基本方法与异环境数据的折算和综合方法相似，区别在于将环境因子换成相似因子。

步骤S3.3：根据数据类型，将折算和综合后的异环境数据和相似产品数据再进行综合，得到寿命型数据和成败型数据的验前信息。

步骤S4：多阶段多类型数据进行折算和综合之后根据数据类型确定先验分布和超参数。具体包括如下步骤：

步骤S4.1：根据总体分布中未知参数的取值范围选取先验分布π(θ)，基本方法如下：

对于[0,1]取值的参数，选用贝塔分布B(a,b)作为先验分布，超参数为(a,b)；对于(-∞,+∞)取值的参数，选用正态分布N(μ,σ²)作为先验分布，超参数为(μ,σ²)；对于(0,+∞)取值的参数，选用伽马分布Γ(z,τ)或倒伽马分布IΓ(z,τ)作为先验分布，超参数为(z,τ)。

各类验前分布的密度函数、期望和方差如表6所示。

表6常见验前分布的概率密度函数、核、期望和方差

步骤S4.2：根据验前数据，确定验前分布的超参数，基本方法为利用验前信息，产生自助样本，得到各样本对应分布的参数的极大似然估计值样本，计算出参数的均值和方差，该均值和方差应等于对应参数的验前分布的均值和方差，联立方程组，得到寿命分布参数的验前分布的超参数。

具体地，①对指数寿命型数据选取先验分布和超参数：

指数分布的参数失效率λ,λ∈(0,+∞)的共轭验前分布为伽马分布，即π(λ)～Γ(z,τ)，超参数(z,τ)分别表示指数型验前试验的故障数和试验总时间，其值由验前信息确定。

确定超参数(z,τ)的方法如下：

利用先验信息和Bootstrap法获得失效率的若干估计值，记为计算先验均值和先验方差并令它们等于伽马分布Γ(z,τ)的期望和方差，联立方程组解之，得到超参数的估计值：

②对正态或对数正态型数据选取先验分布和超参数：

正态或对数正态分布的参数μ和σ²未知且相关时，其参数(μ,σ²)的联合共轭分布是正态—逆伽马分布，即v₀,μ₀,k₀为超参数，由验前数据确定。

确定超参数的方法如下：

对于均值，先利用验前信息和Bootstrap法得到均值μ的验前均值和验前方差分别令其等于μ的验前分布的期望和方差

同理，对于方差，利用验前信息和Bootstrap法得到方差σ²的验前均值和验前方差分别令其等于σ²的验前分布的期望和方差，即

联立上面两个方程组，解得

③对威布尔寿命型数据选取先验分布和超参数：

威布尔分布的参数有(m,η)。根据m的取值范围，分两种情况讨论：

情况一：若已知失效率是递减的，则0<m<1，选取贝塔分布为m的先验分布，即取m～β(a,b)，λ＝(1/η)^m～Γ(z,τ)，则在m和λ是相互独立的假设下，它们的先验分布为：

π(λ,m)∝λ^z-1e^-τλm^a-1(1-m)^b-1 (10)

超参数(a,b)和(z,τ)的确定方法如下：

对于参数m，利用验前信息和Bootstrap法得到m的验前均值和方差分别令其等于m的验前分布m～β(a,b)的期望和方差，即

求上述方程组，得到超参数(a,b)：

同样，对于参数λ，利用验前信息和Bootstrap法的λ的验前均值和方差分别令其等于λ的验前分布λ＝(1/η)^m～Γ(z,τ)的期望和方差，即

求上述方程组，得到超参数(z,τ)：

情况二：已知失效率是递增的，则m>1，则可选取如下先验分布，使m′＝m-1服从伽马分布。由于m>1，则令m′＝m-1的先验分布为伽马分布，即取m′＝m-1～Γ(a,b)，λ＝(1/η)^m～Γ(d,τ)，则在m和λ是相互独立的假设下，它们的先验分布为

π(λ,m)∝λ^z-1e^-τλ(m-1)^a-1e^-b(m-1) (14)

超参数(a,b)和(z,τ)的确定方法与情况一相同。

④对成败型数据选取先验分布和超参数：

二项分布的参数为成功率或可靠度R：共轭验前分布为贝塔分布，即π(R)～B(a₀,b₀)，超参数a₀,b₀分别表示验前试验的成功次数和失败次数，其值由验前信息确定。验前分布密度函数：

确定超参数(a₀,b₀)的方法如下：

利用先验信息和Bootstrap法获得成功概率(可靠度)的若干估计值，记为计算先验均值和先验方差并令它们等于贝塔分布B(α,β)的期望和方差，利用如下方程组得到超参数(a₀,b₀)：

⑤当验前信息为专家信息时，选取先验分布方法如下：

当现场数据为成败型时，选取贝塔分布为先验分布；

当现场数据为寿命型时，选取负对数伽马分布为先验分布。

确定超参数的方法：

已知专家对可靠度的点估计和置信下限R_0L(置信度1-β)，当R的先验分布为贝塔分布时，由下列方程组确定超参数(a,b)

当R的先验分布为负对数伽马分布时，令任务时间t₀，τ为试验总时间，则等效任务数为η＝τ/t₀。由下列方程组确定超参数(z,η)

步骤S5：确定先验分布π(θ)和超参数之后，根据Bayes理论及相关模型，进行数据融合和可靠性评估。具体包括如下步骤：

步骤S5.1：在已知总体分布类型和现场数据的情况下，确定各分布的似然函数p(x|θ)。

各类总体分布的似然函数如表7所示。

表7常见总体分布的似然函数

步骤S5.2：根据Bayes理论，确定验后分布的核或分布密度函数。

验后分布密度函数为：

其中，c＝1/∫_Θp(x|θ)π(θ)dθ为一个与θ无关的正则常数，p(x|θ)π(θ)为验后分布的核，即

π(θ|x)∝p(x|θ)π(θ)

步骤S5.3：根据验后分布的核，判断验后分布的类型，得到参数的后验分布密度函数。

具体地，①对于二项分布，试验数据：(n,s)，n为试验次数，s为成功次数。

参数：可靠度R,R∈[0,1]，其共轭验前分布为贝塔分布，即π(R)～Be(s₀,f₀)，超参数s₀,f₀分别表示验前试验的成功次数和失败次数，其值由验前信息确定。

验前分布密度函数

似然函数为

验后分布的核

π(R|n,s)∝p(n,s|R)π(R)

R的验后分布：也为贝塔分布，即π(R|n,s)～Be(s₀+s,f₀+f)，其pdf为

②指数分布

试验数据：(z,τ)，z＝r为失效数，τ＝T_r为试验总时间

参数：失效率λ,λ∈(0,+∞)，其共轭验前分布为伽马分布，即π(λ)～Ga(z₀,τ₀)，超参数z₀,τ₀分别表示指数型验前试验的失效数和试验总时间，其值由验前信息确定。

似然函数为：

p(z,τ|λ)＝

λ验后分布：也为伽马分布，即，其pdf为

③正态分布和对数正态分布

试验数据：

正态x₍₁₎≤x₍₂₎≤…≤x_(n)

对数正态lnx₍₁₎≤lnx₍₂₎≤…≤lnx_(n)

若为寿命数据，则均应大于0。这里只考虑完全样本的情况。

参数：存在以下三种情况：

情况(一)：μ未知(σ²已知)

μ的验前共轭分布为正态分布，即超参数μ₀,分别表示验前数据的均值和方差，由验前数据估计得到。

验前分布的核为：

似然函数为：

μ的共轭验后分布仍为正态分布，即其中：

情况(二)：σ²未知(μ已知)

σ²的验前共轭分布为逆伽马分布，即π(σ²)～IGa(a,b)，超参数a,b由验前数据估计得到。

验前分布密度函数：

似然函数：

σ²的共轭验后分布仍为逆伽马分布

情况(三)：μ和σ²未知且相关

参数(μ,σ²)的联合共轭分布是正态—逆伽马分布，即v₀,μ₀,k₀为超参数，由验前数据确定。

其中：即

则联合验前分布

似然函数为

则联合后验分布也为正态－逆伽马分布

其中：

μ的后验分布为自由度为v_n的学生t分布，即

σ²的后验分布为逆伽马分布，即

④威布尔分布

试验数据：

t₁≤t₂≤…≤t_r，即n个样本的前r个故障时间

两参数威布尔分布的分布函数：

参数：分以下三种情况讨论：

情况(一)：形状参数m已知，η未知

试验数据可以表示成则这些数据为来自指数分布Exp(λ)，λ＝(1/η)^m的前r个次序统计量，则参数λ的验后分布可以按照指数分布给出。

情况(二)：失效率是递减的，即0＜m＜1，选取贝塔分布为m的先验分布。

取m～β(a,b)，λ＝(1/η)^m～Γ(d,τ)，则在m和λ是相互独立的假设下，先验分布为

π(λ,m)∝λ^d-1e^-τλm^a-1(1-m)^b-1

给定试验样本数据条件下，似然函数为

其中：表示以为失效时间的总试验时间。

则m和λ的联合后验分布为

即

其中：C₁为与m和λ无关的常数，其值为

由此可得，m的后验分布为

λ的后验分布为

情况(三)：失效率是递增的，即m＞1，则可选取如下先验分布，使m′＝m-1服从伽马分布，即m′＝m-1～Γ(a,b)，选取λ＝(1/η)^m～Γ(d,τ)，则在m和λ是相互独立的假设下，先验分布为

π(λ,m)∝λ^d-1e^-τλ(m-1)^a-1e^-b(m-1)

给定试验样本数据条件下，似然函数为

其中：表示以为失效时间的总试验时间。

则m和λ的联合后验分布为

即

其中：C₂为与m和λ无关的常数，其值为

由此可得，m的后验分布为

λ的后验分布为

步骤S5.4：在已知参数的后验分布密度函数的情况下，根据总体分布，计算各可靠性评估指标的后验期望估计值，作为其贝叶斯估计值。

具体地，①对于二项分布，R的贝叶斯估计(后验期望)：

②指数分布

λ的贝叶斯估计(后验期望)：

平均寿命的贝叶斯估计(后验期望)：

任务时间t₀时的可靠度的贝叶斯估计(后验期望)：

③正态分布和对数正态分布

情况(一)

μ的贝叶斯估计：

任务时间t₀时的可靠度的点估计：

情况(二)

σ²的贝叶斯估计：

任务时间t₀时的可靠度的贝叶斯估计(后验期望)：

其中：c₁为与s²有关的归一化常数。

情况(三)

μ的期望后验估计：

σ²的期望后验估计：

R的Bayes估计：

其中：c₂为一个与s²有关的常数。

④威布尔分布

情况(一)

同指数分布。

情况(二)

m的贝叶斯估计：

λ的贝叶斯估计：

可靠度的贝叶斯估计：

情况(三)

m的贝叶斯估计：

λ的贝叶斯估计：

可靠度的贝叶斯估计：

本发明还提供了一种产品可靠性综合评估系统，包括：

模块M1：将多阶段多类型数据进行分类与预处理；

模块M2：将分类与预处理之后的多阶段多类型数据进行一致性检验；

模块M3：一致性检验后对多阶段多类型数据进行折算和综合；

模块M4：多阶段多类型数据进行折算和综合之后根据数据类型确定先验分布和超参数；

模块M5：确定先验分布和超参数之后，根据Bayes理论及相关模型，进行数据融合和可靠性评估。

本发明又提供了一种存储有计算机程序的计算机可读存储介质，计算机程序被处理器执行时实现上述方法的步骤。

本发明继续提供的一种产品可靠性综合评估设备，包括上述的产品可靠性综合评估系统或者上述的存储有计算机程序的计算机可读存储介质。

本领域技术人员知道，除了以纯计算机可读程序代码方式实现本发明提供的系统及其各个装置、模块、单元以外，完全可以通过将方法步骤进行逻辑编程来使得本发明提供的系统及其各个装置、模块、单元以逻辑门、开关、专用集成电路、可编程逻辑控制器以及嵌入式微控制器等的形式来实现相同功能。所以，本发明提供的系统及其各项装置、模块、单元可以被认为是一种硬件部件，而对其内包括的用于实现各种功能的装置、模块、单元也可以视为硬件部件内的结构；也可以将用于实现各种功能的装置、模块、单元视为既可以是实现方法的软件模块又可以是硬件部件内的结构。

以上对本发明的具体实施例进行了描述。需要理解的是，本发明并不局限于上述特定实施方式，本领域技术人员可以在权利要求的范围内做出各种变化或修改，这并不影响本发明的实质内容。在不冲突的情况下，本申请的实施例和实施例中的特征可以任意相互组合。