随着加速试验技术的快速发展,利用加速退化建模技术获取产品性能退化数据、外推产品寿命信息成为当前可靠性工程领域的研究热点。类似常规应力退化建模,加速应力退化建模技术主要有基于退化轨迹、退化量分布和随机过程的建模^{[1, 2, 3, 4, 5]}。

基于退化轨迹和退化量分布2种建模技术使用最早且较为成熟,但因其将产品性能退化描述为确定性过程,未考虑产品个体随时间发生退化的随机性,使其工程应用范围受限。以Wiener过程为代表的随机过程模型,充分考虑了产品在正常工作或试验环境中所发生性能退化的随机性和动态性,可以较好地描述产品的真实退化过程,成为当前加速退化建模的重要技术手段。

Wiener过程是一种具有平稳、便于计算与分析特性的随机过程模型,具有广泛的应用前景。当产品的性能退化过程总体均匀变化,而个体差异随时间增加而逐步增大时,都可用Wiener过程来描述^{[6, 7]}。目前,基于Wiener过程的加速退化建模的重点及难点主要在于以下2点:

1) 如何处理非线性数据。主要解决方法是进行坐标变换,将非线性数据组变换为线性数据组。Whitmore和Schenkelberg^[8]最早提出2类时间尺度变换模型来解决Wiener过程在非线性退化数据建模中的问题,实例表明了模型具有较好的实用性。在此基础上,文献[9, 10, 11]采用时间尺度变换模型,进一步开展了基于非线性Wiener过程的加速退化数据建模与应用研究。

2) 如何考虑个体之间的性能退化差异。由于产品在制造、材料及环境中受随机因素的影响表现出个体之间退化速率不一,存在一定差异性,因此需要研究个体差异对产品寿命评估的影响。主要方法是对Wiener过程的系数进行随机化处理。文献[12, 13]对Wiener过程的2个系数都进行随机化处理,假定扩散系数服从逆伽马分布的前提下漂移系数服从正态分布,但假定的分布类型未经分布假设检验,通过最大期望(Expectation Maximization,EM)算法估计参数,计算复杂,难以推广;文献[14, 15]仅将Wiener过程的漂移系数看作服从正态分布的随机变量,给出了考虑个体差异的产品寿命模型,模型拟合性好,便于计算。

基于此,本文在步进应力加速退化试验场合下,采用时间尺度变换模型将非线性退化数据转换为线性数据,将Wiener过程的漂移系数随机化处理,建立考虑个体差异的加速退化数据模型,采用两步极大似然估计法确定模型中的未知参数,较EM算法更适合处理加速退化场合下的退化数据。

假设产品的性能退化过程X(t)可以用Wiener过程进行描述:

由于X(t)符合齐次马尔可夫过程,X(t)的均值和方差都是时间的线性函数,因此,Wiener过程一般用于描述产品线性退化过程,即退化量X(t)与时间t呈线性关系。

设产品退化失效的阈值为l。当X(t)首次达到l时,认定产品失效,则产品寿命T满足:

经推导可知,首次达到l时的寿命服从逆高斯分布,则产品的可靠度和概率密度函数分别为

1) 试样在一组加速应力(即温度)作用下性能发生退化且退化过程不可逆。

2) 试样在各加速应力下的退化失效机理、模式都保持不变。

3) 试样在各加速应力下的性能退化参数仅有一个且测量间隔和次数相同。

4) 某时刻试样的剩余寿命仅与该时刻已累计退化失效水平和当时应力有关。

根据试验方案,预先确定一组高于正常应力S₀的加速应力S₁ < S₂ < … < S_n,试验初期,将若干试样放置在初始应力S₁下进行试验并记录相关性能退化参数,达到时刻t₁后将应力提高至S₂下继续试验,达到时刻t₂后将应力提高至下一应力,如此进行,直到将应力提高至S_n下完成试验,其中,t_i-1(i=1,2,…,n)为应力S_i-1到S_i的转换时刻。试验应力施加过程如图 1所示。

图 1 步进应力加速退化试验载荷历程Fig. 1 Step-stress accelerated degradation test loading process

图选项

已知现有m个样本在n个加速应力下进行步进应力加速退化试验,每个应力水平下各测量K次性能退化参数,则在应力S_i下第j个试样第k次测量时刻为t_i,k^j,性能退化量为X(t_i,k^j)。

经步进应力加速退化过程分析可知(见图 2),第j个试样在应力S_i下的性能退化量初始值是前一步应力下性能退化量的末尾值,即X(t_i,0^j)=X(t_i-1,k^j),则基于Wiener过程的步进应力加速退化数据建模如下:

图 2 步进应力加速退化过程Fig. 2 Step-stress accelerated degradation process

图选项

在应力S_i下第j个试样第k次测量的性能退化量较第k-1次测量的性能退化量增加了ΔX(t_i,k^j)=X(t_i,k^j)-X(t_i,k-1^j),时间增量为Δt_i,k^j=t_i,k^j-t_i,k-1^j,则根据Wiener过程性质可知:

对于非线性退化数据,可采用时间尺度变换模型,将其转换为线性数据。Whitmore和Schenkelberg^[8]首次提出了时间尺度变换模型应为时间t的非负单调递增函数,常见函数为

本文采用时间尺度变换模型,将非线性数据组[t,X(t)]转换为线性数据组[τ,Y(τ)],则可将式(1)改写为

令τ=Λ(t),Y(τ)=X(t),则产品的寿命可靠度和概率密度函数分别为

为了体现出试样个体之间的差异化,将漂移系数λ随机化处理,即假定λ服从正态分布N(μ_λ,σ_λ²),可将式(8)改写为

对于步进应力加速退化试验而言,一般认为Wiener过程的漂移系数λ与应力S有关,而扩散系数σ与应力S无关,则加速退化模型一般可表示为^[16]

则温度应力下的步进应力加速退化模型为

分析式(16)可知,同一温度应力作用下的试样退化速度并无差异。而考虑到试样个体之间的差异客观存在,应采用基于随机变量的Arrhenius模型来描述试样个体之间的退化差异性,则在应力S_i下第j个试样的加速退化模型可表示为

令a~N(μ_a,σ_a²),则应力S_i下的考虑了个体差异的漂移系数λ_i可表示为

若不考虑个体差异,令σ_a²=0,则式(17)变为传统的Arrhenius模型。

针对经时间尺度变换模型变换后的、考虑个体差异的线性退化数据[τ,Y(τ)],由式(6)可知：

根据式(19)、式(20),建立极大似然估计函数:

采用时间尺度变换模型Ⅰ,将τ=t^c、ΔY(τ_i,k^j)=ΔX(t_i,k^j)与式(17)代入式(21),可得

分别令式(22)关于a_j、σ_B²的1阶偏导数为零,即

理论上,根据收集到的数据组[t_i,k^j,X(t_i,k^j)]及所建的模型,可以解出式(24)、式(25),求出a_j、σ_B²。但进一步分析可知,式(24)、式(25)的求解其实依赖于b、c的取值。因此,采用两步极大似然估计法来求解未知参数Θ。

第1步 估计b、c。采用MATLAB软件中的fminsearch函数可解决上述问题。以b、c为变量，lnL(Θ)为优化函数,分别先给b、c赋初值b₀、c₀,然后进行二维遍历搜索,直到函数lnL(Θ)取得最大值时停止搜索,此时所返回的b、c值即为所求的b、c。

第2步 估计a_j(j=1,2,…,m)、σ_B²。将第1步中所确定的b、c代入式(24)、式(25)中,即可求出a_j、σ_B²。同时,μ_a、σ_a²估计公式如下:

将文献[17]进行步进应力加速退化试验激光器的数据信息记为真实值,经仿真得到的性能退化数据见表 1,性能退化轨迹见图 3。已知激光器以电流为其性能特征参数,在工作过程中受温度应力影响较为敏感,其性能退化速度与温度之间的关系符合Arrhenius模型。激光器的正常工作温度为25℃,现有12个试样进行步进应力加速退化试验,步进温度应力分别为25、50和75℃,每个温度应力下分别测量5次数据,每次测量时间间隔为150h。当激光器工作电流增加量达到规定阈值(l=10mA)时,可认定其退化失效。

图 3 步进应力加速退化轨迹Fig. 3 Step-stress accelerated degradation trajectories

图选项

1) 判断标准

将本文提出的考虑个体差异且可处理非线性数据的方法记为M1;将文献[14]提出的考虑个体差异且只适用于线性数据的方法记为M2;将文献[18]提出的未考虑个体差异且只适用于线性数据的方法记为M3。

依据赤池信息量准则(Akaike Information Criterion,AIC)和总体均方误差(TMSE)判断各方法的优劣^[14]。当AIC或TMSE的值越小时,说明方法的拟合程度越好,计算公式如下:

2) 两步极大似然估计

令Arrhenius模型中参数b的初值为2600,时间尺度变换模型Ⅰ中参数c的初值为1,采用MATLAB软件中的fminsearch函数进行遍历搜索,再将所返回的b、c值代入式(24)~式(26),求得未知参数Θ,见表 2。

由表 2可知,当c取1.035时,可将非线性退化数据转换为线性退化数据。同时,采用CvM检验法^[19],对转换前后的性能退化增量数据,进行显著水平为0.05的正态分布假设检验,结果表明,仅有转换后的数据接受假设,进一步验证了对非线性退化数据进行转换的必要性。

3) 评估结果

根据本文所建模型,计算出Wiener过程的漂移系数、扩散系数、产品寿命的点估计值(MTTF)和置信度为95%的区间估计(CI)，见表 3。

4) 模型优劣判断

分别计算出不种方法的AIC值(见表 4),发现M1的AIC最小,说明M1的模型拟合性更好。在时间0~10000h内分别取10、100、1000个测量点,计算不种方法的TMSE的值(见表 4),发现M1的TMSE值始终最小,说明M1的寿命估计值更接近真实值。

另外,从可靠度曲线(见图 4)可以看出,M1可靠度曲线早期几乎与真实值曲线重合且低于M2和M3可靠度曲线;中期低于真实值、M2和M3可靠度曲线;末期略低于真实值和M2曲线,略高于M3可靠度曲线,说明M1的模型拟合性较好且可靠性评估略保守,与工程实践做法相符。

图 4 不同方法的可靠度曲线Fig. 4 Reliability curves of different methods

图选项

从概率密度曲线(见图 5)可以看出,M1概率密度曲线较M2概率密度曲线更接近于真实值曲线,而M3概率密度曲线较真实值曲线偏差较大,与3种方法所计算出的TMSE值相符,说明3种方法的模型拟合优劣排序为:M1,M2,M3,M1最优。

图 5 不同方法的概率密度曲线Fig. 5 Probatility density curves of different methods

图选项

1) 考虑了个体退化差异,将Wiener过程中的漂移系数随机化处理并假定为正态分布,这种处理方法较传统未考虑个体差异和未进行漂移系数随机化处理的方法更优。

2) 采用两步极大似然估计法,对模型中未知参数进行估计,可以较好地克服传统极大似然估计法的局限性,求解未知参数的最优值。

3) 优先采用时间尺度变换模型对非线性退化数据进行线性化处理,使得处理后的数据更符合Wiener过程特性。

4) 实例分析中,从3种不同方法的AIC值、TMSE值、可靠度曲线以及概率密度曲线进行对比分析,验证了本文所提出的方法的优势,具有一定的工程应用价值。