近年来, 随着航空电子产品可靠性不断提高, 额定应力试验产生的故障数据越来越少甚至无故障, 需采用加速试验的方法, 在保持失效机理不变的条件下, 加速产品失效过程, 以减少试验时间, 缩短研制进度。同时, 由于航空电子产品运行环境复杂, 包括工作应力(如电流、电压等)和环境应力(如温度、湿度、振动和冲击等), 导致传统简单应力试验无法真实复现实际工作环境。因此, 需要研究多应力加速寿命评估方法, 准确获取航空电子产品实际使用环境下的寿命数据。

目前, 许多学者对多应力加速模型进行了研究。魏高乐等^[1-6]提出了一些基于多应力综合加速模型的可靠性试验和评估方法, 进行可靠性分析和寿命估计。Li等^[7-11]结合传统加速模型, 提出适用于特定设备或元器件的双应力和三应力加速寿命模型, 从而预计可靠性指标。但是, 现有多应力加速模型均忽略了应力之间的相互影响, 没有考虑各应力之间的真实相关性, 使得试验预估寿命出现较大偏差。

但是, 被测产品的失效时间是主观记录的, 可能不准确。同时, 当加速寿命数据外推到正常应力时, 模型耦合项选取存在主观因素, 其中的不确定性也会转移并影响评估结果的准确性, 降低评估的可信度。模糊理论是处理不确定性问题的重要方法之一, 目前在可靠性试验领域中被广泛研究。李鹏、Xu等人分别从加速模型、退化模型等不同试验和统计模型出发, 研究模型的模糊化处理方法, 解决模型中的不确定性^[12-17]。因此, 可以在多应力广义耦合加速模型的构建中引入模糊理论进行分析, 以解决模型构建中出现的不确定性问题。

对于在多应力加速试验中获得的试验数据, 需要利用参数估计方法进行数据处理和可靠性指标外推。极大似然估计(maximum likelihood estimation, MLE)提供了一种利用给定观测数据评估模型参数的方法, 被广泛用于估计简单应力下加速模型的未知参数。但多应力模型中参数过多, 多种应力间存在复杂关系, 使得似然函数求解困难, 采用传统计算方法精度不够、效率不高, 甚至无法求解, 因此可以引入智能算法进行计算和参数求解。国内外有许多学者通过智能算法求解参数估计模型, 获取多应力加速模型。Yang等^[18-21]结合神经网络、模拟退火、粒子群优化算法、支持向量回归和遗传算法等进行模型求解, 以预估复杂应力下产品的寿命和测量误差。但很多算法均存在计算效率不高, 容易过早收敛导致结果误差极大等问题。

为解决上述问题, 提出一种考虑广义耦合的多应力加速寿命评估方法, 基于失效机理, 结合模糊映射和模糊相关性明确应力的类型和耦合原则。从试验类别出发, 提出4类MLE方法进行参数估计, 与自适应正弦余弦算法(sine cosine algorithm, SCA)相结合, 求解多应力广义耦合加速模型。以某型低频滤波放大器为例, 设计三应力加速试验, 讨论了考虑实际环境应力耦合和不考虑应力耦合的寿命分布偏差, 得出了考虑实际应力耦合的寿命更短, 符合现实预期。同时, 与传统牛顿迭代法对比, 采用多目标自适应SCA可解决MLE的多参数求解困难问题, 并且求解效率更高, 精度控制在5%以内, 验证了所提方法的有效性。

(1) 单应力加速模型

加速寿命模型的本质是用加大应力而不改变失效机理的方法, 加速产品的失效和退化。根据加速寿命试验的结果, 可以推出正常使用状态或降额使用状态下的产品寿命^[22]。典型简单应力加速模型如下。

1) 温度加速应力模型

Arrhenius模型起源于化学领域, 现今广泛应用于加速试验模型建模^[23], 一般用于描述产品寿命特征量与所施加温度应力之间的关系, 即

(1) $\eta(T)=\alpha_0 \exp \left(\frac{E_a}{k_B T}\right)$

式中: η(·)是产品的寿命特性; α₀为常数, 称产品系数, 与产品自身失效机理和试验方法有关; E_a为化学反应的激活能; k_B是玻尔兹曼常数, 即8.617 eV/℃; T为温度应力, 一般为开尔文温度, 单位为K。

该模型表明寿命与逆反应速率成正比, 为便于进行多应力推广, 温度加速模型可以表示为

(2) $\eta(T)=\alpha_0 \exp \left(\alpha_1 X_1\right)$

式中: X₁为等效的温度应力。

2) 湿度加速应力模型

描述湿度作为加速应力时产品的寿命特征与应力水平之间的关系通常有功率湿度模型、指数湿度模型和倒数指数湿度模型^[24], 分别可表示为

(3) $\left\{\begin{array}{l}\eta(\mathrm{RH})=\lambda(\mathrm{RH})^{-n} \\ \eta(\mathrm{RH})=\exp (\lambda \cdot \mathrm{RH}) \\ \eta(\mathrm{RH})=\exp (\lambda / \mathrm{RH})\end{array}\right.$

式中: RH为所采用的电应力水平; λ和n为与失效模式相关的常数。

3) 电流电压加速应力模型

逆幂律模型通常描述电应力作为加速应力时产品的寿命特征量和应力水平之间的关系^[23], 一般可表达为

(4) $\eta(S)=A S^\gamma$

式中: S为所采用的电应力水平; γ为与产品失效模式相关的常数。

(2) 双应力加速模型

基础的寿命-应力模型只涉及一种应力, 这些模型适用于多数加速试验, 但有些试验包含多种加速应力。因此, 必须推广单应力到双应力乃至多应力, 以更好适应工程需要。在单应力的基础上, 可通过建立物理模型、统计外场及试验数据, 获得一些双应力加速模型, 常用的双应力模型有广义Eyring模型和其相关拓展模型。

广义Eyring模型起源于量子力学和统计力学, 表明了各种外部环境因素对化学反应速率的影响, 在寿命评估中常用于描述含有温度和其他应力与寿命指标^[23]的关系:

(5) $\eta\left(T, X_2\right)=\frac{\alpha_0}{T} \exp \left(\frac{\alpha_1}{k_B T}\right) \cdot \exp \left[X_2\left(\alpha_2+\frac{\alpha_3}{k_B T}\right)\right]$

式中: α₀, α₁, α₂, α₃为待估系数; X₂为非温度应力, 工程应用可省略1/T。同时, 许多学者提出了不同的基于广义Eyring模型的双应力模型, 例如Peck模型、固体断裂Zhurkov模型^[23]等。为便于推广, 将Eyring模型简化为如式(2)的形式:

(6) $\eta\left(X_1, X_2\right)=\alpha_0 \exp \left(\alpha_1 X_1\right) \cdot \exp \left(\alpha_2 X_2\right) \cdot \exp \left(\alpha_3 X_1 X_2\right)$

式中: α₃X₁X₂是X₂与温度应力之间可能的耦合项, 需要根据实际应力情况判定存在与否, 但应力之间的耦合情况不一定只有简单相乘关系, 因此需要推广具有广义耦合关系的多应力加速模型。

(1) 受测产品的寿命分布满足威布尔分布

加速试验一定程度上会使产品的后期失效提前发生, 从而导致恒定故障率的指数分布在加速试验后期偏差过大, 因此失效率随时间递增的威布尔分布更适用于航空电子产品。

在加速寿命试验下, 寿命累积分布函数为

(7) $F(t, \eta, \beta)=1-\exp -\left[\left(\frac{t}{\eta}\right)^\beta\right], t>0$

式中: β是形状参数; η是尺度参数, 一般用来衡量产品的寿命指标。

(2) 试验中多应力水平与寿命满足一定关系

在多应力加速寿命试验中, 要确保加速应力水平不会超过产品可承受的最大应力水平。对于威布尔分布, 特征寿命遵循多应力加速模型:

(8) $\ln \eta=f(\boldsymbol{S})$

式中: f(·)表示应力与对数寿命指标的关系; S为标准化后的多应力向量, S=(S₁, S₂, …, S_N)。由于应力测量单位的不同, 在多应力情况下可能存在很大差异, 无法进行统一度量。因此, 需要变换应力来统一不同测量单位下应力的关系, 消除测量单位的影响^[24]。S_i是第i个应力X_i的标准化形式, 可以表示为

(9) $S_i=S\left(X_i\right)=\frac{X_i-X_{i 0}}{X_{i H}-X_{i 0}}, i \in[1, N]$

式中: X_i0为第i个应力X_i的正常基准应力水平; X_iH为第i个应力X_i的加速基准应力水平。则满足Arrhenius模型、逆幂律模型、指数模型的加速应力标准化过程为

(10) $\left\{\begin{array}{l}X_{i(1)}=\frac{1}{S_i^{\prime}} \\ X_{i(2)}=\ln S_i^{\prime} \\ X_{i(3)}=S_i^{\prime}\end{array}\right.$

式中: S′_i为未标准化应力; 温度应力标准化满足X_i(1), 湿度应力和电流应力标准化满足X_i(2), 某些其他应力标准化满足X_i(3)。

(3) 产品的失效机理在应力水平的各种组合下不发生变化。即在加速模型中, 不同加速应力设置情况下, 同一产品的形状参数β相同。

基于Arrhenius模型, 提出一种考虑广义耦合关系的多应力加速模型, 模型由单应力项和多应力广义耦合项组成。因此, 根据式(9)对应力进行标准化后, 多应力反应速率关系可表示为

(11) $\eta\left(\alpha_0, \alpha_1, \cdots, \alpha_p\right)=\alpha_0 \prod\limits_{k=1}^N \exp \left[\eta_k(S)\right]$

式中: α₀, α₁, …, α_p为未知参数, $p=\sum _{k=1}^N C_N^k$表示所有可能存在的应力单独作用和广义耦合作用的项数; η_k(S)表示k种应力的作用, 如η₁(S)表示N种单应力的单独作用, η₂(S)到η_N(S)表示可能存在的双应力耦合到多应力耦合项。

对等式(11)两边取自然对数, 转换为如式(8)形式的考虑广义耦合的多应力加速模型:

(12) $\ln \eta\left(\alpha_0, \alpha_1, \cdots, \alpha_p\right)=\ln \alpha_0+\sum\limits_{k=1}^N \eta_k(S)$

式中: 应力项表达为

(13) $\left\{ \begin{array}{*{35}{l}} {{\eta }_{1}}(S)=\sum\limits_{r=1}^{N}{{{\alpha }_{r}}}{{S}_{r}} \\ {{\eta }_{2}}(S)=\sum\limits_{\begin{array}{l} r=s=1, r\ne s, r<s, \\ u=N+1, N+2, \cdots , N+C_{N}^{2} \end{array}}^{C_{N}^{2}}{{{\alpha }_{u}}}{{\varphi }_{rs}} \\ {{\eta }_{3}}(S)=\sum\limits_{\begin{array}{l} r=s=v=1, r\ne s\ne v , r <s<v\\ w=N+C_{N}^{2}+1, N+C_{N}^{2}+2, \cdots , N+{{C}_{N}^{3}} \end{array}}^{C_{N}^{3}}{{{\alpha }_{w}}}{{\varphi }_{rsv}} \\ \vdots \\ {{\eta }_{N}}(S)={{\alpha }_{p}}{{\varphi }_{rsv\cdots N}} \\\end{array} \right.$

式中: φ为不同应力间的广义耦合关系; φ_rs为第r和第s个应力间可能存在的广义耦合关系。其他为对应数目应力可能存在的关系, 由于多应力耦合项不一定是应力的单纯相乘关系, 因此引入一种广义表达表示耦合效应。常用函数形式有指数函数和幂函数, 其参数在一定范围变化时, 可以拟合绝大部分统计模型, 广泛用于回归分析和轨迹建模。从而建立广义耦合表达为

(14) $\varphi=\prod\limits_{k_1=1}^{K_1} \exp \left(\omega_{k_1} S_{k_1}\right) \cdot \prod\limits_{k_2=1}^{K_2} S_{k_2}{ }^{\omega_{k_2}} \cdot \prod\limits_{k_3=1}^{K_3}\left(\omega_{k_3}\right)^{S_{k_3}}$

式中: φ为广义耦合表达的通用形式; ω_k1、ω_k2、ω_k3为常数, 表示幂次; k=K₁+K₂+K₃, 为该耦合项中施加应力的数目。耦合项需要结合现实工程情况和产品自身失效机理, 考虑综合应力耦合效应进行判定, 具体形式后续可通过模糊相关性计算确定。

在实际工程应用中, 针对施加应力, 应结合物理失效机理判断应力施加类型和应力耦合项的存在与否。因此, 需分解产品功能结构, 建立可靠性框图^[25], 描述产品中所有要素的相互关系, 从而获取功能-结构和失效-应力的映射关系。

考虑可靠性目标下的层次关系如图 1所示, 其中每层之间的联系相互映射。

通过层次化映射分析可初步定性获取产品的敏感应力, 为进一步确定实际多应力加速试验的应力类型, 需建立定量模型筛选敏感应力: ①调用失效权重向量、各个故障模式的百分比, 确定模糊关系; ②通过模型获取模糊隶属度, 计算模糊映射矩阵, 确定产品对每种应力的敏感度。

(1) 失效-应力模糊关系表示

构造失效-应力模糊映射关系如下: 直积空间X(=应力)×Y(=失效模式), 映射$\widetilde{R}: X \times Y \rightarrow[0, 1]$称应力和失效的模糊关系, 为X×Y上的模糊集, 即$\widetilde{R}$∈F(X×Y)。的隶属函数为$\mu_{\tilde{R}}(x, y)$, 函数值$\mu_{\tilde{R}}(x, y)$代表有序对(x, y)具有关系的程度。如果只取0、1, 则$\widetilde{R}$便是普通关系。如果X=Y, 则称$\widetilde{R}$为X上的二元模糊关系。如果X和Y中的元素是有限的, 即X={x₁, x₂, …, x_m}, Y={y₁, y₂, …, y_n}, 则称X到Y的模糊关系为模糊矩阵, 并表示为

(15) $\widetilde{\boldsymbol{R}}=\left[\begin{array}{cccc}r_{11} & r_{12} & \cdots & r_{1 n} \\ r_{21} & r_{22} & \cdots & r_{2 n} \\ \vdots & \vdots & \ddots & \vdots \\ r_{m 1} & r_{m 2} & \cdots & r_{m n}\end{array}\right]$

式中: $r_{i j}=\mu_{\tilde{R}}\left(x_i, y_j\right)$, i=1, 2, …, m, j=1, 2, …, n。全体m×n模糊矩阵记为M_m×n。

(2) 模糊隶属度和应力敏感度计算

每种应力类型对应的模糊映射定义为

(16) $\widetilde{\boldsymbol{R}}_i=\left[r_{i 1}, r_{i 2}, \cdots, r_{i n}\right]=\left[\mu_{\tilde{R}}\left(x_i, y_j\right)\right]_{1 \times n}$

式中: i=1, 2, …, m, 表示应力类型; j=1, 2, …, n, 表示被检测产品的失效模式；$\widetilde{\boldsymbol{R}}_i$为一个维数为1×n的模糊隶属度矩阵, 表示第i种应力与产品失效模式的模糊映射关系；r_ij=$\mu_{\bar{R}}\left(x_i, y_j\right)$为当施加第i个应力时对第j个失效模式的影响程度, 即模糊隶属度。

根据第1.2.1节, 定义每种故障模式的故障原因分别为疲劳、磨损、结构破坏、腐蚀等K类, 每种故障原因由s(s≤m)个应力导致, 则第j个故障模式的第i个应力的隶属度为

(17) $r_{i j}=\frac{\sum\limits_{k=1}^K c_{j k} l_{s k}}{\sum\limits_{s=1}^m l_{s k}}$

式中: 第j个故障模式的故障原因百分比为c_jk; 第k个故障原因的应力影响程度由影响系数l_sk表示。

定义第s个应力的敏感度为R_s(s≤i), 即

(18) $R_s=\sum\limits_{j=1}^n \lambda_j r_{s j}$

式中: λ_j为第j个失效模式占产品总失效的比例, 即失效模式百分比。应力敏感度可以定量应力对产品的影响程度, 从而确定在多应力加速试验中施加应力的类型。

结合第1.2.1节提出的故障机理和环境应力的映射关系, 可获取产品的失效与应力的联系。当产品的一个失效与两个以上的应力相关时, 同时考虑耦合效应和失效-应力映射关系, 综合判定耦合项的存在与否。参考DO-160G, 给出加速试验中综合应力耦合效应关系示例^[26-27], 如图 2所示。

因此, 不同应力间可能存在一定关系, 考虑到环境效应, 包含3种应力关系^[26-27]: ①应力间存在耦合, 作用于产品失效机理, 表现为加强机械损伤或加速功能退化; ②应力间相互弱化, 影响产品失效速率, 一般不同时施加; ③应力不相关, 同时施加不影响寿命。

由于试验观测等存在随机不确定性, 模型建立等存在认知不确定性, 为消除不确定性对耦合项系数的影响, 采用模糊数来表达施加应力与寿命指标的关系。模糊数包括几种类型, 如区间值数、三角数、梯形数、指数等^[28]。这里采用最常用模糊化方法, 即三角模糊法, 其隶属函数为

(19) $\widetilde{R}_s=u_{A_k}(x)=\left\{\begin{array}{l}\frac{x-a}{b-a}, a \leqslant x<b \\ 1, x=b \\ \frac{c-x}{c-b}, b<x \leqslant c \\ 0, \text { 其他 }\end{array}\right.$

因此, 考虑广义耦合的多应力加速模型的模糊表达为

(20) $\ln \eta\left(\alpha_0, \alpha_1, \cdots, \alpha_p\right)=\ln \alpha_0+\sum\limits_{k=1}^N u_{A_k}(x) \eta_k(S)$

结合试验数据和模糊相关度计算, 可确定具体的应力广义耦合项的表达形式。

相关度衡量样本中两个随机变量是否线性相关, 为确定广义耦合多应力加速模型, 需要计算不同形式耦合项与寿命指标的相关性。传统皮尔森相关系数用来估计一定量实数样本参数之间的相关性, 不适用于模糊数据, 因此提出一种衡量应力与寿命指标模糊相关性的系数计算方法。

对于加速试验获取的寿命数据, 必须考虑模糊区间的长度是否相同, 若将重心和模糊区间长度以相等的权重组合, 则组合可能会超过1或-1, 不满足相关度定义。因此, 为获得合理的模糊相关性, 本文定义所有应力耦合项φ为样本X, 寿命指标ln η为样本Y, 两变量的模糊相关度ε为

(21) $\varepsilon=\beta_1 \varepsilon_1+\beta_2 \varepsilon_2, \beta_1+\beta_2=1$

其中:

(22) $\left\{\begin{array}{l}\varepsilon_1=\frac{\sum\limits_{i=1}^n\left(c x_i-\overline{c x}\right)\left(c y_i-\overline{c y}\right)}{\sqrt{\sum\limits_{i=1}^n\left(c x_i-\overline{c x}\right)^2} \sqrt{\sum\limits_{i=1}^n\left(c y_i-\overline{c y}\right)^2}} \\ \varepsilon_2=\frac{\sum\limits_{i=1}^n\left(\left\|x_i\right\|-\left\|\overline{x_i}\right\|\right)\left(\left\|y_i\right\|-\left\|\overline{y_i}\right\|\right)}{\sqrt{\sum\limits_{i=1}^n\left(\left\|x_i\right\|-\left\|\overline{x_i}\right\|\right)^2} \sqrt{\sum\limits_{i=1}^n\left(\left\|y_i\right\|-\left\|\overline{y_i}\right\|\right)^2}}\end{array}\right.$

式中: 每组样本的模糊数为x_i和y_i, 数据长度以模糊隶属面积衡量, 定义为模糊面积||x_i||和||y_i||; 由此表达模糊重心为(cx_i, cy_i), 面积为||x_i||=area(x_i)、||y_i||=area(y_i)。

同时, (β₁, β₂)是一个模糊概念, 用来衡量模糊重心和面积的权重, 为简单起见可采用(0.9, 0.1)、(0.8, 0.2)、(0.7, 0.3)、(0.6, 0.4)等来表达权重比例, 从而解释两个变量间相关性的变化。本研究认为重心权重更高, 因此应用(β₁, β₂)=(0.7, 0.3)来计算模糊数的相关系数。

在模糊相关度的定义下, ε越接近1, 正关联越强, 两个变量越模糊正相关; 反之ε越接近于-1, 负关联越强, 两个变量越模糊负相关; 当ε接近0时, 两个变量几乎不相关。从而通过模糊相关性选取最相关广义耦合项, 以获取与实际情况最符合的形式。

综合第1.2节和第1.3节可得基于失效机理的多应力耦合加速试验的施加应力类型和模型耦合项确立判定流程。

在实际工程试验中, 一般获取的失效寿命数据有随机截尾数据、全数样本数据、定时截尾数据和定数截尾数据^[22]。因此, 从方法的通用性出发, 提出了基于多应力广义耦合加速寿命多种样本数据的MLE方法。

(1) 完全样本试验

完全样本试验指样本全部失效才停止试验, 这种试验可获得较完整的数据, 统计分析结果也较好, 但所需时间较长。完全样本试验中, 有n个试验样本, 试验到所有器件全部失效t_n时刻停止试验, 依照失效时间先后记录, 试验停止前的失效时间为t₁, t₂, …, t_i, …, t_n, MLE为

(23) $L=\prod\limits_{i=1}^n \frac{\beta}{\eta^\beta} t_i^{\beta-1} \exp \left(\frac{t_i}{\eta}\right)^\beta$

拓展到多应力的情况, 共设计q组多应力试验, 在第h组试验中, 试件有n_h个, 失效时间可表示为t_h1, t_h2, …, t_hi, …, t_hnh(h≤q), 则多应力组合S_h下产品失效的对数似然函数为

(24) $\begin{gathered}\ln L_h=\ln \left\{\prod\limits_{i=1}^{n_h} \frac{\beta}{\eta_h^\beta} t_{h i}^{\beta-1} \exp \left(-\frac{t_{h i}}{\eta_h}\right)^\beta\right\}= \\ n_h \ln \beta-n_h \beta \ln \eta_h+(\beta-1) \sum\limits_{i=1}^{n_h} \ln t_{h i}-\eta_h^{-\beta} \sum\limits_{i=1}^{n_h} t_{h i}^\beta\end{gathered}$

将式(24)代入多应力广义耦合加速模型后, 即可拓展到多应力广义耦合的情况, 获得在完全样本下所有应力组合S_h={S₁, S₂, …, S_q}下的极大似然函数:

(25) $\begin{gathered}\ln L_q\left(\alpha_0, \alpha_1, \cdots, \alpha_p, \beta\right)=\sum\limits_{h=1}^q \ln L_h= \\ \ln \beta \sum\limits_{h=1}^q n_h-\beta \sum\limits_{h=1}^q n_h \ln \eta_h\left(\alpha_0, \alpha_1, \cdots, \alpha_p\right)+ \\ (\beta-1) \sum\limits_{h=1}^q \sum\limits_{i=1}^{n_h} \ln t_{h i}-\sum\limits_{h=1}^q \eta_h\left(\alpha_0, \alpha_1, \cdots, \alpha_p\right)^{-\beta} \sum\limits_{h=1}^q \sum\limits_{i=1}^{n_h} t_{h i}^\beta\end{gathered}$

其中, ln η_h(α₀, α₁, …, α_p)为第h组应力组合下的广义耦合加速模型。

(2) 随机截尾试验

随机截尾试验中, 有n个试验样本, r个样本失效, c个样本未失效, 但中途撤离试验, 称为删失样本, n=r+c, 失效时间为t₁, t₂, …, t_r, 删失时间为τ₁, τ₂, …, τ_c。则似然函数为

(26) $L=\prod\limits_{i=1}^r \frac{\beta}{\eta^\beta} t_i^{\beta-1} \exp \left[-\left(\frac{t_i}{\eta}\right)^\beta\right] \cdot \prod\limits_{j=1}^c \exp \left[-\left(\frac{\tau_j}{\eta}\right)^\beta\right]$

拓展到多应力耦合的情况中, 试件的失效时间和删失时间可分别表示为t_h1, t_h2, …, t_hrh和τ_h1, τ_h2, …, τ_hch, 且n_h=r_h+c_h。

过程同完全样本试验, 式(26)两边取对数, 代入多应力模型, 即可获得在随机截尾试验下多应力广义耦合加速模型的极大似然函数:

(27) $\begin{gathered}\ln L_q\left(\alpha_0, \alpha_1, \cdots, \alpha_p, \beta\right)=\sum\limits_{h=1}^{q} \ln L_h= \\ \ln \beta \sum\limits_{h=1}^q r_h-\beta \sum\limits_{h=1}^q r_h \ln \eta_h\left(\alpha_0, \alpha_1, \cdots, \alpha_p\right)+ \\ (\beta-1) \sum\limits_{h=1}^q \sum\limits_{i=1}^{r_h} \ln t_{h i}-\sum\limits_{h=1}^q \eta_h\left(\alpha_0, \alpha_1, \cdots, \alpha_p\right)^{-\beta} \cdot \\ \sum\limits_{h=1}^q\left(\sum\limits_{i=1}^{r_h} t_{h i}^\beta+\sum\limits_{j=1}^{c_t} t_{h j}^\beta\right)\end{gathered}$

(3) 定时截尾和定数截尾试验

定时截尾试验中, 有n个试验样本, 试验到t₀时刻停止试验。依照失效时间先后记录试验停止前的失效时间为t₁, t₂, …, t_r≤t₀(r < n), 则极大似然函数为

(28) $L=\prod\limits_{i=1}^r \frac{\beta}{\eta^\beta} t_i^{\beta-1} \exp \left\{-\frac{1}{\eta^\beta}\left[\sum\limits_{i=1}^r t_i^\beta+(n-r) t_0^\beta\right]\right\}$

拓展到多应力情况, 试验到t_h0时刻停止试验, 最后一次记录时间为t_hr, 且t_hr≤t_h0。则多应力组合下产品失效的对数似然函数可表示为

(29) $\begin{gathered}\ln L_h=r_h \ln \beta-r_h \beta \ln \eta_h+(\beta-1) \sum\limits_{i=1}^{r_h} \ln t_{h i}- \\ \frac{1}{\eta_h^\beta}\left(\sum\limits_{i=1}^{r_h} t_{h i}^\beta+\left(n_h-r_h\right) t_{h 0}^\beta\right)\end{gathered}$

将式(29)代入多应力耦合加速模型, 得定时截尾试验下, 所有应力组合的极大似然函数模型为

(30) $\begin{gathered}\ln L_q\left(\alpha_0, \alpha_1, \cdots, \alpha_p, \beta\right)=\sum\limits_{h=1}^q \ln L_h= \\ \ln \beta \sum\limits_{h=1}^q r_h-\beta \sum\limits_{h=1}^q r_h \ln \eta_h\left(\alpha_0, \alpha_1, \cdots, \alpha_p\right)+ \\ (\beta-1) \sum\limits_{h=1}^q \sum\limits_{i=1}^{r_h} \ln t_{h i}-\sum\limits_{h=1}^q \eta_h\left(\alpha_0, \alpha_1, \cdots, \alpha_p\right)^{-\beta} \cdot \\ \sum\limits_{h=1}^q\left(\sum\limits_{i=1}^{r_h} t_{h i}^\beta+\left(n_h-r_h\right) t_{h 0}^\beta\right)\end{gathered}$

同理, 在定数截尾试验中, 试验到第r_h个试验件失效后停止试验, 时间为t_hr。则与定时截尾试验类似, 极大似然函数可表示为

(31) $\begin{gathered}\ln L_q\left(\alpha_0, \alpha_1, \cdots, \alpha_p, \beta\right)=\sum\limits_{h=1}^q \ln L_h= \\ \ln \beta \sum\limits_{h=1}^q r_h-\beta \sum\limits_{h=1}^q r_h \ln \eta_h\left(\alpha_0, \alpha_1, \cdots, \alpha_p\right)+ \\ (\beta-1) \sum\limits_{h=1}^q \sum\limits_{i=1}^{r_h} \ln t_{h i}-\sum\limits_{h=1}^q \eta_h\left(\alpha_0, \alpha_1, \cdots, \alpha_p\right)^{-\beta} \cdot \\ \sum\limits_{h=1}^q\left(\sum\limits_{i=1}^{r_h} t_{h i}^\beta+\left(n_h-r_h\right) t_{h r}^\beta\right)\end{gathered}$

综合上述4种数据截尾方式, 可得不同截尾试验下的MLE通用模型为

(32) $\begin{gathered}\ln L_q\left(\alpha_0, \alpha_1, \cdots, \alpha_p, \beta\right)=A \ln \beta-\beta \sum\limits_{h=1}^q B \ln \eta_h\left(\alpha_0, \alpha_1, \cdots, \alpha_p\right)+ \\ C(\beta-1)-D \sum\limits_{h=1}^q\left[\eta_h\left(\alpha_0, \alpha_1, \cdots, \alpha_p\right)\right]^{-\beta}\end{gathered}$

其中, A、B、C、D的形式如表 1所示。

由于模型参数过多, 使用传统方法过程繁琐、计算量大、求解困难, 所得参数估计模型的误差较大, 置信度不足。因此, 引入多目标自适应SCA, 解决多应力广义耦合情况下传统方法难以求解的问题, 同时在一定程度上减小误差。

传统似然函数进行MLE时, 对每个待求参数求偏导, 得到一组非线性方程组, 求解较为困难。因此, 在多参数似然函数的求解过程中, 为提高对模型参数求解效率, 将问题转化为优化算法解决多参数无约束优化问题。

一般优化算法的核心是构造自适应变量目标, 针对非线性参数求解问题, 可以采用令每一个方程的绝对值尽量靠近0的方法来进行多参数优化求解。得到自适应目标后, 可以进行优化迭代, 来最小化目标并找到根, 如下所示：

(33) $\begin{gathered}F=\left[\frac{\partial \ln L_q\left(\alpha_0, \alpha_1, \cdots, \alpha_p, \beta\right)}{\partial \beta}\right]^2+ \\ \sum\limits_{i=0}^p\left[\frac{\partial \ln L_q\left(\alpha_0, \alpha_1, \cdots, \alpha_p, \beta\right)}{\partial \alpha_i}\right]^2\end{gathered}$

算法目的即转化成寻找无约束优化模型的最小值F_min, 建立多参数MLE优化模型, 解决非线性求解困难问题。

SCA主要提供了两种技术, 分别是全局搜索和局部搜索^[29], 两种搜索方法相结合可广泛探索最优解的搜索空间, 通过最小化目标函数进行多参数估计, 避免算法迭代过程中陷入局部最优或者过早收敛, 过程如图 3所示。

(1) 自适应过程控制

为探索和利用之间取得适当的平衡, 以实现全局最优解。SCA进行自适应调整和迭代, 如下所示:

(34) $c_1=b\left(1-\frac{t}{T}\right)$

式中: t是当前迭代的次数; T是最大迭代次数; b是常数值。c₁控制全局搜索和局部探索的转化, 当c₁较大时进行全局搜索, 较小时进行局部探索。

(2) 局部探索控制

正余弦函数在[-2, 2]内的影响用数学模型sin c₂和cos c₂来表示, 该模型描述了正余弦函数的切换范围需要一个解来更新其在不同解中的位置。因此, 定义c₂∈[0, 2π], 以搜索局部空间内外的随机位置, 实现局部探索。

(3) 更新个体位置

为了更新参数估计和迭代的位置, 使用等式更新最优解:

(35) $P_i^{j+1}=\left\{\begin{array}{l}P_i^j+c_1\left|c_3 Q_i^j-P_i^j\right| \sin c_2, c_4<0.5 \\ P_i^j+c_1\left|c_3 Q_i^j-P_i^j\right| \cos c_2, c_4 \geqslant 0.5\end{array}\right.$

式中: P_i^j表示当前解在第j次迭代时在第i维的位置; c₁, c₂, c₃, c₄是4个随机数; Q_i^j是目标位置, 绝对值表示幅度; 参数c₁表示向下一个位置区域的移动方向; 参数c₂在[2, 2π]之间, 表示运动应该朝向或远离目标的程度; 参数c₃表示目标的随机权重分数, 以随机强调(c₃>1)或淡化(c₃ < 1)目标在确定距离方面的影响; 参数c₄是[0, 1]之间的随机数, 表示正弦和余弦函数之间的切换。

(4) 迭代停止

通过一组初始随机解开始优化过程, 然后缓存迄今为止获得的最佳解决方案, 将其指定为目标点, 并在考虑所选最佳解决方案的同时更新其他解决方案。随着迭代次数的增加, 正弦和余弦函数的范围被刷新以保持对搜索空间的利用。当迭代次数达到最大迭代次数时, 停止优化过程。

通过上述过程, 可以获取多应力耦合MLE模型的未知参数, 从而进行寿命评估。

提出基于多应力广义耦合的加速寿命评估方法, 包括三大步骤: ①建立对应含未知参数的多应力广义耦合加速模型, 并确定敏感应力和广义耦合项形式; ②根据不同试验数据类型建立多参数估计模型, 并结合SCA进行参数求解; ③基于加速模型外推正常工作应力下产品的寿命指标。整体思路如图 4所示。

以民用机载电子设备中的射频低噪声放大模块为对象, 进行加速寿命评估和分析。

首先, 进行模糊映射层次分析和敏感应力确定, 确定试验方案, 进行结构-机理-应力的映射分析, 结果如图 5所示。

从可靠性预计手册和外场相似产品数据处获取失效模式和数据, 在工程实际中, 一般采用故障模式、机理及影响分析(failure modes, mechanisms and effects analysis, FMMEA)建立映射关系矩阵。

本文FMMEA结果参考文献[22], 调用该低噪声放大模块各个故障模式的百分比和失效率, 通过结合先验信息和第1.2.1节中的方法, 计算模糊隶属度, 获取失效机理-环境应力模糊映射矩阵:

(36) $\widetilde{\boldsymbol{R}}=\left[\begin{array}{llllll}0.1266 & 0.0951 & 0.1689 & 0.1008 & 0.0422 \\ 0.0646 & 0.0432 & 0.0814 & 0.0517 & 0.0209 \\ 0.0327 & 0.0216 & 0.0272 & 0.0086 & 0.0086 \\ 0.0381 & 0.0252 & 0.0379 & 0.0175 & 0.0109 \\ 0.0571 & 0.0406 & 0.0691 & 0.0450 & 0.0180\end{array}\right]$

最终, 计算得施加的敏感应力结果如表 2所示, 为便于工程试验, 选择施加电流、温度和湿度。

考虑综合应力耦合效应, 对于电子元器件和电路板, 在高温高湿条件下, 高温加速水汽渗透使电路板表面离子的活性增强, 强化化学腐蚀和电化学效应, 高温下电子活性增强, 同时电子元器件对温度敏感, 工作电流存在浮动和变化^[25], 存在耦合关系; 但在潮湿情况下加大电流应力, 一般不会加快损伤或退化进程, 认为不耦合。

因此, 在该三应力试验中, 设置应力耦合关系如表 3所示。

根据加速模型, 考虑实际情况耦合和不考虑耦合的三应力加速模型分别为

(37) $\begin{gathered}\ln L_1\left(\alpha_{10}, \alpha_{11}, \cdots, \alpha_{15}\right)=\ln \alpha_{10}-\alpha_{11} S_1- \\ \alpha_{12} S_2-\alpha_{13} S_3-\alpha_{14} \varphi_{12}-\alpha_{15} \varphi_{13}\end{gathered}$

(38) $\ln L_2\left(\alpha_{20}, \alpha_{21}, \alpha_{22}, \alpha_{23}\right)=\ln \alpha_{20}-\alpha_{21} S_1-\alpha_{22} S_2-\alpha_{23} S_3$

为便于工程计算, 基于试验结果对所有3次方以下的备选模型进行模糊相关性分析, 备选结果如表 4所示。

最后, 确定最优广义耦合项形式为

(39) $\begin{gathered}\ln L_{21}\left(\alpha_{10}, \alpha_{11}, \alpha_{12}, \alpha_{13}, \alpha_{14}, \alpha_{15}\right)=\ln \alpha_{10}-\alpha_{11} S_1-\alpha_{12} S_2- \\ \alpha_{13} S_3-\alpha_{14} S_1^2 \sqrt{S_2}-\alpha_{15} S_1 \exp \left(S_3\right)\end{gathered}$

备选广义耦合项和传统简单耦合项形式分别为

(40) $\begin{gathered}\ln L_{22}\left(\alpha_{20}, \alpha_{21}, \alpha_{22}, \alpha_{23}, \alpha_{24}, \alpha_{25}\right)=\ln \alpha_{20}-\alpha_{21} S_1-\alpha_{22} S_2- \\ \alpha_{23} S_3-\alpha_{24} \exp \left(S_1\right) \sqrt{S_2}-\alpha_{25} \exp \left(S_1\right) \sqrt{S_3}\end{gathered}$

(41) $\begin{gathered}\ln L\left(\alpha_{30}, \alpha_{31}, \alpha_{32}, \alpha_{33}, \alpha_{34}, \alpha_{35}\right)=\ln \alpha_{30}-\alpha_{31} S_1-\alpha_{32} S_2 \\ \alpha_{33} S_3-\alpha_{34} S_1 S_2-\alpha_{35} S_1 S_3\end{gathered}$

给该产品设计三应力加速寿命试验, 3种应力(温度、湿度和电流)的应力设置方法参考文献[4]和文献[30], 共设置6组应力组合。

三应力试验的每个应力组合下的试件数量为20个, 直到所有产品全部失效。该产品的正常工作温度、湿度和电流要求分别为30℃, 50%, 10 mA, 即该组实验的标准法向应力水平。在加速试验中, 施加应力水平如表 5所示。

在三应力试验中投入120件产品, 按评估方案分为6组, 进行定数截尾试验, 直到产品失效截止试验, 统计相关失效数据。

为验证所提出评估方法的有效性, 需要对试验寿命数据进行统计分析, 不同应力水平下的失效散点图如图 6(b)所示。为确保模型正确性, 需要对数据进行拟合优度检验。经赤池信息量准则(Akaike information criterion, AIC)、贝叶斯信息准则(Bayesian information criterion, BIC)和柯尔莫哥洛夫-斯米诺夫(Kolmogorov-Smirnov, K-S)检验综合判断数据是否满足模型前提假设, 分布的概率图、残差及其拟合优度如图 6所示, AIC、BIC检验结果如表 6所示。

同时, 分别对6组试验数据进行K-S检验, 个案数为20, 原假设显著性水平p₀=0.05, 检验结果如表 7所示。

结果表明, 6组数据的显著性水平均大于0.05, 所以接受原假设。通过拟合优度综合判定, 试验数据最满足威布尔分布, 可以进行模型建模、参数估计和寿命评估。

根据式(33)得到MLE模型, 通过SCA和传统牛顿迭代法求解待求参数。

设置算法参数为N=50, a=2, T_max=10 000;搜索空间设置为[-12, 12]。为尽量消除随机数的影响, 进行10次迭代运算, 计算其平均值。图 7给出了三应力最优广义耦合参数求解的所有适应值收敛曲线和平均适应值收敛曲线, 10 000次迭代后平均适应值收敛于1.93×10^-3, 表明该方法求解模型具有良好的迭代收敛性。

经自适应SCA计算, 取10次迭代结果的平均值与传统方法计算结果进行对比, 可得算法计算的结果的均方误差(mean square error, MSE)控制在5%以内, 计算效率、精度较高, 并且收敛稳定, 表 8列出了4种模型的求解结果。经进一步计算, 可得每个应力水平下参数的估计值和标准误差, 如表 9所示。

最后, 可以通过求解出的多应力耦合模型外推在正常应力水平下的寿命指标, 给出了最优广义耦合模型、备选广义耦合模型、传统简单耦合模型的牛顿迭代法和自适应SCA的寿命估计结果, 并与可靠性预计中的历史寿命进行对比, 结果如图 8所示。

如表 10所示, 针对该低噪声放大模块, 进行三应力加速试验, 分别对比3类考虑应力耦合和应力不耦合的模型求解结果, 考虑应力广义耦合作用的结果更贴近于实际预计情况, 且经模糊相关性分析的最优广义耦合模型最满足历史寿命数据。根据先验信息, 该产品预计寿命约为3 200天, 与历史经验寿命数据相比, 不考虑应力间耦合结果存在较大误差, 由此可得出应力间的耦合作用对该产品寿命估计精度有很大影响, 必须考虑实际耦合作用。同时, 基于自适应SCA的多参数估计方法标准误差均在5%以内, 说明了参数预测的准确。

本文提出一种多应力广义耦合加速寿命评估方法, 结合实际航空电子产品复杂运行环境, 首先基于传统加速应力模型, 推导出一种考虑多应力广义耦合的加速模型描述分布参数与应力水平的关系, 并基于模糊映射和模糊相关性, 结合产品和工程应用实际给出广义耦合项的设置判定方法, 为长寿命产品的寿命评估提供理论基础。其次, 提出不同试验数据类型的多应力广义耦合参数MLE方法, 并采用自适应SCA获取未知参数, 解决多应力加速模型多参数估计难题。然后, 以某型低噪声放大模块为例, 对方法的正确性和有效性开展验证, 经验证该方法更贴近于实际工作环境且误差控制在5%以内, 有效解决无法快速准确获取国产航空电子零部件的寿命指标问题。

在小样本和数据缺乏的情况下对样本数据的质量要求较高, 小样本容错率较小; 由于自适应算法的随机性, SCA有可能出现结果偏差较大的问题, 需要进行多次迭代求解。后续有待于对其中存在的不确定性和随机性问题进行进一步的研究和探讨。