随着材料科学与生产制造工艺的发展进步，军工领域出现了越来越多的高可靠性长寿命电子产品。此类产品在额定工作应力下性能退化缓慢，仅仅依靠性能检测数据已经难以快速有效地进行可靠性评估。因而加速退化试验成为可靠性评估的有效途径。通常认为温度是导致电子部件退化失效的最主要环境应力类型，温度对多种失效机理具有较好的加速效果，如电子元器件的二次俘获、表面电荷扩散、电迁移等，因而常将温度作为加速应力。

通常认为，由于生产工艺、贮存条件、人员操作等多方面主客观因素，产品的个体差异性是客观存在的。目前，考虑个体差异的随机过程模型在可靠性领域得到了较多研究应用，如最常用的3类随机过程模型：Wiener过程^[1-2]、Gamma过程^[3]、Inverse Gaussian过程模型^[4-5]等，但以上研究都是在额定应力下的产品可靠性评估。如若考虑温度加速应力作用下的电子部件二元退化建模，考虑个体差异的随机过程可靠性评估则需要解决以下2个重要问题：一是如何确定此类模型参数与加速应力的约束关系；二是如何进行相关性建模。

确定模型参数与加速应力的约束关系是基于加速退化数据可靠性评估的前提，也是一个研究难点。目前，采用的模型通常都是随机过程基本模型，即针对的都是不考虑个体差异的随机过程模型参数与加速应力的约束关系。此种情况下，解决方法大多是根据经验假定模型参数与加速应力的关系，比如Wiener过程模型，就存在着2个不同假定：假定漂移参数与应力有关而扩散参数与加速应力无关^[6-9]；假定漂移参数和扩散参数都与加速应力相关^[10-11]等。显然，这易导致评估结果可信度不高。另一种解决方法就是通过加速因子不变原则推导模型参数与加速应力的约束关系，加速因子不变原则是指为保证加速因子具有工程应用性，则加速因子应当是一个不随时间变化的常数，比如文献[12-13]推导得到了Wiener过程模型漂移参数和扩散参数与加速应力的约束关系。以上，都是加速应力作用下，针对不考虑个体差异的随机过程模型的，而文献[14]采用加速因子不变原则推导得到了加速应力下，考虑个体差异的Wiener过程模型参数与应力的约束关系。

二元加速退化建模中，相关性建模通常采用基于Copula函数的建模方法^[15]，也有采用基于马氏距离的二元降维建模方法^[16]的。在基于Copula函数的建模方法中，基本都将Copula函数参数视为固定值，即忽略了产品个体之间存在的随机相关性。文献[17]提出了随机相关的建模思想，但只应用在额定应力下，没有进一步考虑加速应力作用下Copula函数参数与加速应力的关系。

综上，为解决电子部件二元加速退化可靠性评估难题，建立基于贝叶斯方法的二元随机相关性能退化可靠性评估框架。采用Wiener RDV(Random Drift-Volatility)模型进行边缘分布建模，在此基础上，构建基于Copula函数的多个随机相关模型，采用两阶段贝叶斯参数估计方法进行参数估计(通过OpenBUGS软件实现基于MCMC(Markov Chain Monte Carlo)的参数估计)，进行相关性模型选择，最后基于蒙特卡罗方法进行可靠度计算，并应用实例验证了所提方法有效性。

1 二元随机相关可靠性分析框架

基于Wiener过程的二元加速性能退化可靠性分析，以边缘退化建模为基础，利用Copula函数实现二元相关关系建模。主要包括3部分内容：①边缘分布建模和相关性建模；②两阶段贝叶斯参数估计；③基于蒙特卡罗方法的可靠性评估。首先，利用样本性能退化数据，采用RDV模型建立边缘分布模型，并考虑产品间具有随机相关性，建立基于Copula函数的相关性模型；然后，采用两阶段贝叶斯参数估计方法对边缘分布模型进行参数估计，在此基础上得到边缘分布退化增量的累积概率分布值，将其作为Copula函数的输入值，并通过Copula函数模型选择，确定相对最优的相关模型；最后，在参数估计和模型选择基础上，采用蒙特卡罗方法建立可靠度函数曲线，完成产品的可靠性评估。基本框架如图 1所示。

本文采用如下3个假定。

假定1  同型产品，不仅在性能退化阶段具有个体差异，而且在产品之间的相关关系上也可能具有个体差异性。

假定2  同型产品相关关系的差异性，是指Copula函数参数的差异性，即同型产品具有相同类型的Copula函数，但Copula函数的参数值不同。

假定3  加速应力可能会影响相关关系，但是应力不改变Copula函数类型，只改变Copula函数参数值。

2 Wiener RDV模型

文献[14]给出了一种考虑个体差异的Wiener RDV模型，该模型将性能退化过程的均值和方差都确定为随机变量，该模型描述为

式中：X(t)为退化量; B(·)为标准Brown运动函数; μ为漂移参数; σ为扩散参数; Λ(t)=t^c为时间函数；a、b、e、h为超参数；N(h, )表示均值为h，方差为的正态分布。模型的累积分布函数为

其中：D为阈值；F_2a(·)为自由度为2a的t分布函数。文献[14]采用加速因子不变原则推导得到了该模型参数与加速应力的约束关系：

式中：a_k、b_k、c_k、e_k和h_k为应力s_k的模型参数；a_l、b_l、c_l、e_l和h_l为应力s_l的模型参数；AF为应力s_k相对于应力s_l的加速因子。由式(3)可采用Arrhenius模型将温度应力s下参数表示为

式中：γ_RDV(1)、γ_RDV(2)、γ_RDV(3)和γ_RDV(4)为待定参数。

3 二元随机相关加速退化数据建模

假定产品具有2个寿命表征参数X₁和X₂，边缘性能退化过程都服从Wiener RDV模型。t_ijk为第i个样本第k个参数第j次测量时刻，X_ijk为t_ijk时刻测量值，k=1, 2；i=1, 2, …, N；j=1, 2, …, M；N为样本数量，M为测量次数。对所有样本进行了相同次数的测量，不同性能参数在相同时刻进行测量。则易得性能退化数据结构为

式中：X_k为第k个性能参数的所有样本退化数据。时间函数Λ_k(t)=t^c_k，Λ_k(0)=0，ΔΛ_k(t)=(t+Δt)_k^c_k－t_k^c_k，ΔX_k为描述第k个性能退化过程在[t, t+Δt]增量区间的退化增量。边缘分布函数分别为F₁(ΔX_ij1)和F₂(ΔX_ij2)，联合分布函数为H(ΔX_ij1, ΔX_ij2)，则由SKlar定理，存在一个Copula函数C(u, v)使得H(ΔX_ij1, ΔX_ij2)=C(F₁(ΔX_ij1), F₂(ΔX_ij2))成立。

考虑随机相关性，可建立基于Copula函数的6个相关模型。模型所用符号说明如下：s_i为作用在第i个样本的应力。模型参数分为性能退化模型参数α和Copula相关性模型参数β，并定义α_k^F、α_k^R和α_k^H分别为第k个性能参数退化模型的固定参数、随机参数和超参数；β^F、β^R和β^H分别为Copula相关性模型参数的固定参数、随机参数和超参数。

1) 模型A

不考虑相关关系个体差异性，Copula函数参数θ是一个固定参数，这是目前通常采用的模型，即

性能退化过程参数α_k^F={α_μk^F, α_σk^F, c_k}，α_μk^F={γ_RDV(2), γ_RDV(3), γ_RDV(4)}，α_σk^F={γ_RDV(1), γ_RDV(2)}，α_k^R={μ_k, σ_k}，α_μk^H={h_{k, si}, e_{k, si}}，α_σk^H={a_σk, b_{σk, si}}，h_{k, si}、e_{k, si}和b_{σk, si}为第i个样本相应的参数，下标s_i表示参数直接与应力相关，i=1, 2, …，N。相关性模型中β_A只含固定参数θ。

模型B~模型F的性能退化模型部分与模型A相同，因而以下内容只描述相关性模型部分。

2) 模型B

不考虑相关关系个体差异性，但是Copula函数参数θ(s_i; β_B^F)与加速应力相关，即

Copula函数模型参数只有固定参数集β_B^F={θ_si, γ_B(1), γ_B(2)}。θ_si为第i个样本在s_i应力下的模型参数，i=1, 2, …, N。

3) 模型C

考虑相关关系存在个体差异性，Copula函数参数是服从正态分布的随机变量，即

式中：θ为随机变量，θ~N(a_θ, b_θ²)；β_C={β_C^R, β_C^H}，β_C^R={θ_i}，β_C^H={a_θ, b_θ}，θ_i为第i个样本参数，i=1, 2, …, N。

4) 模型D

考虑相关关系存在个体差异性，Copula函数参数是服从正态分布的随机变量，同时考虑Copula函数参数的超参数a_θ(s_i; β_D^F)与加速应力相关，即

式中：θ~N(a_θ(s_i; β_D^F), b_θ²)，a_θ(s_i; β_D^F)=exp(γ_D(1)－γ_D(2)/s_i)；β_D={β_D^F, β_D^R, β_D^H}，β_D^R={θ_i}，β_D^H={a_{θ, si}, b_θ}，β_D^F={γ_D(1), γ_D(2)}，a_{θ, si}为第i个样本应力s_i下的a_θ参数，i=1, 2, …, N。

5) 模型E

考虑相关关系存在个体差异性，Copula函数参数是服从正态分布的随机变量，同时考虑Copula函数参数的超参数b_θ(s_i; β_E^F)与加速应力相关，即

式中：θ~N(a_θ, b_θ²(s_i; β_E^F))，b_θ(s_i; β_E^F)=exp(γ_E(1)－γ_E(2)/s_i)；β_E={β_E^F, β_E^R, β_E^H}，β_E^R={θ_i}，β_E^H={a_θ, b_{θ, si}}，β_E^F={γ_E(1), γ_E(2)}，b_{θ, si}为第i个样本应力s_i下的b_θ参数，i=1, 2, …, N。

6) 模型F

考虑相关关系存在个体差异性，Copula函数参数是服从正态分布的随机变量，同时考虑Copula函数参数的超参数a_θ(s_i; β_F^F)和b_θ(s_i; β_F^F)都与加速应力相关，即

式中：θ~N(a_θ(s_i; β_aθ^F), b_θ²(s_i; β_bθ^F))，a_θ(s_i; β_aθ^F)=exp(γ_F(1)－γ_F(2)/s_i)，b_θ(s_i; β_bθ^F)=exp(γ_F(3)－γ_F(4)/s_i)；β_F={β_aθ^F, β_bθ^F, β_F^R, β_F^H}，β_F^H={a_{θ, si}, b_{θ, si}}，β_F^R={θ_i}，β_aθ^F={γ_F(1), γ_F(2)}，β_bθ^F={γ_F(3), γ_F(4)}、a_{θ, si}、b_{θ, si}为第i个样本的a_θ、b_θ参数，i=1, 2, …, N。

4 参数估计及模型选择

采用两阶段贝叶斯参数估计方法分别对边缘性能退化过程以及相关性模型进行参数估计。首先利用样本退化数据估计边缘分布参数，当得到边缘分布参数估计值、后，进而可计算得到退化增量的累积分布函数值F₁(ΔX_ij1|s_i, )、F₂(ΔX_ij2|s_i, )。第2阶段的参数估计就是以F₁(ΔX_ij1|s_i, )、F₂(ΔX_ij2|s_i, )为样本估计Copula函数参数。所有模型后验分布的求解，则采用基于MCMC的参数后验分布抽样方法，具体可以通过OpenBUGS软件编程实现。各参数的先验分布都采用无信息先验(合理区间的均匀分布)。

本文采用常用的Gaussian、Frank、Gumbel和Clayton Copula函数进行模型选择，其分布函数及参数如表 1所示。这4类Copula函数描述了不同的相关关系结构，Gaussian Copula和Frank Copula描述的是对称相关结构，Gumbel Copula函数描述较强上尾相关特征的变量间的相关关系，Clayton Copula函数则描述具有较强下尾相关特征的变量相关关系。

Copula函数模型选择方法，可以采用散点图，偏差信息准则(DIC)值以及Kendall τ的非参数估计值等方法。

当获取参数估计值后，边缘退化过程的累积分布函数值(F₁(ΔX_ij1|s_i, ), F₂(ΔX_ij2|s_i, ))散点图很好地体现了不同Copula函数的相关关系结构。因而根据散点图，可对相关关系进行判断。样本的Kendall τ非参数估计值可通过式(12)计算:

式中：(F_1i, F_2i)表示n个Copula函数的样本。当样本量足够大，无限接近真实值。

DIC值是常用的用于模型选择的一种定量方法，定义为

式中：f(t|)为似然函数；p_D为待估参数数目。即DIC是模型期望偏差与表示模型复杂度的估计参数数目之和。DIC值越小，表明模型拟合越好。

5 基于蒙特卡罗的可靠度计算

根据文献[17]思路，为避免可靠度计算对模型方差值的依赖性，采用蒙特卡罗方法进行可靠度计算。仿真生成足够数量的样本退化数据，在任一时刻，定义该时刻未失效样本数V与总样本数N的比值为该时刻可靠度值。V是指t时刻所有边缘退化过程的退化量都小于阈值的样本数量。

经过两阶段贝叶斯参数估计和模型选择后，边缘分布函数分别为F₁(t|)、F₂(t|)，Copula函数模型为C(F₁(ΔX_ij1), F₂(ΔX_ij2)|)，则基于蒙特卡罗的可靠度算法步骤如下：

步骤1  生成N个二元性能退化过程仿真样本。

① 根据相关性模型C(F₁(ΔX_ij1), F₂(ΔX_ij2)|)，生成M组(F₁(ΔX_ij1), F₂(ΔX_ij2))仿真数据。

② 根据边缘分布模型，由F₁^－1(ΔX_ij1)、F₂^－1(ΔX_ij2)分别计算得到ΔX_ij1、ΔX_ij2。

③ 计算X_ijk=X_i(j－1)k+ΔX_ijk。

④ 重复步骤①~步骤③N次，得到N个二元性能退化过程。

步骤2  计算t_j时刻可靠度值。

① 比较第i个样本退化量X_ij1与D₁，X_ij2与D₂大小。若X_ij1 < D₁, X_ij2 < D₂同时成立，则V=V+1；若X_ij1≥D₁或X_ij2≥D₂有一个满足，则该样本失效，并将阈值作为该样本后续时刻退化量值。

② 对N个样本逐次进行步骤①，得到V最终的累加值。

③ 计算t_j时刻可靠度值为V/N。

④ 重复步骤①~步骤③M次，得到所有测量时刻的可靠度值。

步骤3  重复步骤1和步骤2 L次，得到L个性能退化过程的可靠度曲线。

步骤4  由L个可靠度曲线，进一步可求取可靠度曲线的100(1－α)%置信区间。

① 将L个t_j时刻可靠度值R₁(t_j), R₂(t_j), …R_L(t_j)按升序排列，表示为R₍₁₎(t_j), R₍₂₎(t_j), …R_(L)(t_j)，均值为R(t_j)。

② 计算t_j时刻100(1－α)%置信区间，置信区间下限R_(l)(t_j)，l=L·Φ(2Φ^－1(p)+Φ^－1(α/2))，上限R_(u)(t_j)，u=L·Φ(2Φ^－1(p)+Φ^－1(1－α/2))，p为R₍₁₎(t_j), R₍₂₎(t_j), …R_(L)(t_j)小于R(t_j)的比例。

③ 重复步骤①和步骤②，得到M个测量时刻点可靠度均值及其置信区间R(t₁)[R_(l)(t₁), R_(u)(t₁)], R(t₂)[R_(l)(t₂), R_(u)(t₂)], …, R(t_M)·[R_(l)(t_M), R_(u)(t_M)]。

6 实例应用

电子部件通常随导弹、鱼水雷等装备长期贮存，在长贮期间一个显著特点就是产品个体经历了不同条件的贮存环境、不等次数的检测监测、不同时长的战备值班等，这些差异使得同批次电子部件在性能退化和多元相关关系等方面都呈现出个体差异性。此外，尽管装备型号较多、结构组成差异较大，但是较多的电子器件在不同电子部件的共用程度也较高，本文基于数据拟合的可靠性评估方法有较广的适用性，因而下文将以某型装备滚控电路板为例验证所提方法。

滚控电路板属于板级电子部件，具有多个关键性能参数，是较为典型的多元性能退化产品。文献[18]通过对滚控电路板失效模式分析指出负脉宽与正脉宽之差ΔT和周期T是其2个关键性能参数，并且由预试验分析结果得到加速应力下电路板失效机理不改变的最高应力为130℃。因而文献[18]选择9个样本进行恒定加速退化试验，确定3个加速应力为S₁=90℃，S₂=103℃，S₃=120℃，每个加速应力下3个样本。限于篇幅，详细数据见文献[18]。

6.1 边缘分布参数估计

采用OpenBUGS软件实现贝叶斯参数估计，所得结果如表 2所示。利用边缘退化过程随机参数的后验分布样本值，绘制箱线图，如图 2所示，2个边缘退化过程参数都具有很明显的个体差异性，表明采用RDV模型是合理的。

6.2 随机相关模型参数估计与模型选择

由6.1节得到边缘分布参数估计值后，就可得到边缘退化增量累积分布函数(CDFs)取值散点(F₁(ΔX_ij1|), F₂(ΔX_ij2|))。图 3为3个不同加速应力下样本的散点图，图 4为所有样本的散点图。从图 4(a)和图 4(b)中的直方图表明散点具有较明显的对称性和正相关性，而图 4(c)散点相对分散，对称性和正相关性不易判断。图 4呈现出了较为明显的对称性和正相关性。因而，由假定3“加速应力不改变Copula函数类型”，可确定在待选的4类Copula函数中，能够描述对称正相关结构的Gaussian Copula或Frank Copula模型可能是相对合适的模型。

进一步地，分别采用Gaussian模型A，Frank模型A，Clayton模型A和Gumbel模型A进行相关性建模和参数估计，参数估计值如表 3所示，结果表明Frank Copula函数的DIC值最小。采用式(12)估计得到Kendall τ的非参数估计值为0.272 2，而表 3同时给出了各模型τ值，可见Frank Copula函数τ值最为接近τ的非参数估计值。因而，不论是DIC值还是τ的非参数估计值，Frank Copula函数都是相对最优的。

确定采用Frank Copula函数后，进而考虑随机相关模型。采用模型A~模型F进行建模、参数估计，估计值如表 4所示。由DIC值可以看出，模型C的DIC值最小，但是6个模型的DIC值相差都在5以内，不宜做出判断。采用箱线图进一步判断。图 5(a)表明，无论是模型C、模型D、模型E或者模型F，Copula函数参数都具有明显的个体差异性；而图 5(b)~图 5(d)表明模型B、模型D、模型E和模型F，应力变化对模型参数值影响很小。综上，模型C是6个模型中相对最优的模型。

6.3 可靠性评估

额定应力为25℃。将表 2中RDV模型参数代入式(3)和式(4)，并由式(2)得到边缘退化过程的可靠度函数R₁(t)=1－F_Wrdv1(t)，R₂(t)=1－F_Wrdv2(t)。将R₁(t)、R₂(t)及表 4模型C中参数估计值代入式(8)，可得可靠度函数的解析结果R_a(t)为

式中：

按照第5节基于蒙特卡罗的算法计算可靠度值：采用模型C的参数估计值，并在区间[0, 2×10⁷]上设置测量时刻数M=4 000，单次仿真样本数量N=200，循环次数L=200，得到仿真结果如图 6所示。图 6(a)中R_s(t)是仿真所得L条可靠度曲线，R(t)则是所有仿真结果R_s(t)各个时刻的可靠度均值曲线，R_u(t)和R_l(t)分别是仿真结果R_s(t)各个时刻的可靠度值的95%置信区间的上下限。R_a(t)是解析结果的可靠度曲线。尽管各个测量时刻的置信区间上下限值连线并不严格降低，但整体呈明显下降趋势。

图 6(a)表明，虽然解析结果在仿真结果置信区间范围内，但相比于仿真结果均值曲线，明显低估可靠度值。对于板级电路板等高可靠性长寿命产品的基础组成部件，对其可靠度往往要求保持在0.99，甚至更高，因而计算得到了R(t)和R_a(t)可靠度值降至0.99时的寿命值分别约为194 214 h(约22年)，69 972 h(约8年), 如图 6(b)所示。文献[18]以测试值通常都应服从正态分布为依据，采用了二元联合正态分布进行可靠性评估。可见文献[18]对加速退化数据模型分析过于主观，其计算得到该产品可靠度为0.95时的可靠寿命约为22.69年，相比于本文结果，相对保守。

7 结论

本文研究了加速应力下电子部件二元随机相关性能退化可靠性评估，结论如下：

1) 采用考虑随机影响的Wiener过程模型进行边缘分布建模，能描述性能退化过程的个体差异性，提高性能退化建模精度。Wiener过程模型随机参数的箱线图，能较好地反映性能退化过程的个体差异性，是判断产品是否具有个体差异性的有效方法。

2) 产品个体差异性不仅仅体现在性能退化过程中，而且体现在产品间相关关系的差异。随机相关模型能较好地描述相关关系的差异性，是对考虑随机影响的随机过程可靠性建模地进一步完善。

3) 采用散点图、DIC值和Kendall τ的非参数估计值能较好地实现基于Copula函数的随机相关模型选择。

4) 基于蒙特卡罗的可靠度计算方法，能避免可靠度计算对模型方差值的依赖，具有更好的适用性和准确性。