基于比例关系加速退化建模的设备剩余寿命在线预测

2021-01-26王泽洲陈云翔蔡忠义项华春王莉莉

系统工程与电子技术 2021年2期

王泽洲, 陈云翔, 蔡忠义, 项华春, 王莉莉

(空军工程大学装备管理与无人机工程学院, 陕西西安 710051)

0 引言

加速退化试验是获取设备退化信息高效且经济的手段。利用加速退化数据预测设备剩余寿命,并以此为依据制定维修保障策略已成为当前预测与健康管理领域的重要研究内容[1-5]。

受环境和自身因素的影响,设备的退化过程呈现出显著的随机性,而加速试验更进一步加剧了退化的不确定性,因此采用具备时变不确定特征的Wiener过程刻画加速退化过程具备合理性,并受到研究人员的广泛关注,现已在加速度计、激光器、发光二极管(light emitting diode, LED)等设备的加速退化建模中得到应用和验证[6-10]。唐圣金等[11]率先采用Wiener过程建立了设备的加速退化模型,进而实现了对LED步进加速退化试验的最优设计,然而该方法仅采用了基本线性Wiener过程模型,未考虑到非线性情况与个体差异等因素,导致该方法的适用性不强。Tang等[12]考虑到不同设备间退化过程具有差异性,在文献[11]退化模型的基础上引入随机漂移系数,并基于激光器的加速退化数据预测了其在常应力条件下的剩余寿命,提升了剩余寿命预测的准确性,但该方法仍无法适用于非线性退化设备。进一步,Liao等[13]构建了具备非线性和个体差异性的复合Wiener退化模型,并将其拓展应用于加速退化场景,取得了良好效果,但该模型未能分析测量误差对设备剩余寿命预测的影响。针对文献[13]的不足,Hao等[14]综合分析了非线性、个体差异性以及测量不确定性对设备退化过程的影响,并建立了考虑上述3类因素的加速退化模型,从而有效提升了模型的准确性,拓展了方法的应用范围。

然而,上述研究均未能采用目标设备的现场检测数据来同步更新设备的退化状态,因而限制了剩余寿命预测准确性的提升。蔡忠义等[15]利用贝叶斯原理在线更新设备的漂移系数,实现了对加速退化设备剩余寿命的在线预测。但该方法将扩散系数视为固定常数,忽略了应力越大退化不确定性越强的客观规律,导致该方法的预测准确性不高。进一步,Wang等[16]考虑到加速应力对漂移系数和扩散系数的双重影响,并基于贝叶斯原理实现了对漂移系数和扩散系数的同步更新,在一定程度上提升了剩余寿命预测的准确性。但该方法存在一个潜在假设,即漂移系数与扩散系数必须满足特定共轭先验分布,导致其在不满足先验分布假设时预测结果较差。在上述研究的基础上,文献[17]基于加速因子不变原则建立了考虑漂移系数与扩散系数比例关系的退化模型,并对漂移系数与扩散系数进行了同步更新,且不必满足特定共轭先验部分假设,从而拓展了方法的应用范围,进一步提升了剩余寿命预测的准确性。然而,该方法仅考虑了常应力条件下的退化建模,无法适用于加速退化条件下设备的剩余寿命预测。

基于上述分析可知,目前针对基于加速退化建模的设备剩余寿命在线预测方法研究尚不完善,退化状态更新机制亟待改进。为了实现对常应力条件下设备剩余寿命的准确预测,本文将漂移系数与扩散系数的比例关系引入加速退化建模,并基于卡尔曼滤波原理对其进行同步更新,进一步提升了方法的科学性。针对退化模型参数与滤波先验参数,本文提出基于两步极大似然原理的参数估计方法,确保了参数估计的合理性。此外,制定了基于加速因子不变原则的退化数据折算规则,实现了对常应力条件下设备剩余寿命的准确预测。最后,利用加速度计加速退化试验实测数据,对本文方法进行了验证。

1 带比例关系的设备退化建模

传统的Wiener过程模型多用于描述线性退化设备的性能退化轨迹,其基本表达式为

X(t)=X(0)+λt+σBB(t)

(1)

式中,X(t)为t时刻设备的性能退化量;X(0)为设备的初始性能状态量,且常令X(0)=0;λ为漂移系数,用以表征设备的退化速率;σB为扩散系数,主要体现设备退化轨迹的波动程度;B(t)为标准布朗运动,用以表示退化过程的时变不确定性,且B(t)～N(0,t)。

考虑到测量手段的不确定性,众多研究人员将测量误差引入退化建模,建立了带测量误差的Wiener退化模型:

Y(t)=X(t)+ε

(2)

随着设备的功能结构日益复杂,退化过程呈现出显著的非线性特征,导致传统线性Wiener模型精度逐步下降,进而促使了非线性Wiener模型的产生,其具体表达式为

Y(t)=X(0)+λΛ(t|b)+σBB(Λ(t|b))+ε

(3)

基于加速退化试验过程中设备退化机理不变准则,文献[18-19]提出了加速因子不变原则。在此基础上,文献[17]证明了式(3)中漂移系数与扩散系数的比例关系,具体可表示为

(4)

式中,k为未知常数,且与应力大小无关。

(5)

2 加速退化建模

2.1 加速模型

加速模型是用以描述设备退化特征量与敏感应力之间关系的数学模型。根据敏感应力类型的不同,常采用的加速模型有:

(1) Arrhenius模型(敏感应力:温度)

(6)

(2) Eyring模型(敏感应力:温度)

(7)

(3) Inverse Power模型(敏感应力:电压)

λ(S)=αSβ

(8)

(4) Exponential模型(敏感应力:电压)

λ(S)=αexp(βS)

(9)

不失一般性,加速模型的一般形式可表示为

λ(S)=αρ(S|β)

(10)

式中,ρ(S|β)为应力S的函数。

2.2 加速退化模型

加速退化试验主要包含恒定应力加速退化试验、步进应力加速退化试验和序进应力加速退化试验3类[15,20-21]。其中，恒定应力和步进应力试验是应用最为广泛的加速退化试验方式,其理论方法成熟,评估效果较好;序进应力试验目前仍处于初步研究阶段,该方法对试验环境要求严苛,统计模型过于复杂,试验成本高,应用较少。因此,本文主要针对恒定应力和步进应力加速退化试验进行分析,并基于本文提出的带比例关系的设备退化模型构建设备加速退化模型。

(1) 恒定应力加速退化模型

恒定应力加速退化试验中,全体样本被分为若干测试组,并置于不同应力条件下进行试验,具体如图1所示。

整个试验过程中,不同测试组对应的应力条件恒定不变。在满足失效机理不变的前提下,假设恒定应力加速退化试验共有M个应力,且S1

(11)

图1 恒定应力加速退化模型Fig.1 Constant stress accelerated degradation model

(2) 步进应力加速退化模型

步进应力加速退化试验中,全体样本被置于相同的应力条件下进行试验,而试验应力随时间逐步升高或降低,且在两次变化之间应力条件保持恒定,具体如图2所示。

图2 步进应力加速退化模型Fig.2 Step stress accelerated degradation model

在满足失效机理不变的前提下,假设步进应力加速退化试验共有M个应力,且S1

(12)

3 参数估计

假设共有N个测试样本,则Ym,i, j=Y(tm,i, j)表示第m个应力条件下,第i个测试样本在第j个监测时刻对应的性能退化数据。其中,i=1,2,…,nm;而j=1,2,…,lm,i,lm,i表示第m个应力条件下第i个样本的总监测次数。令Ym,i表示第m个应力条件下第i个样本的全部性能退化数据,则Ym=[Ym,1,Ym,2,…,Ym,nm]表示第m个应力条件下试验样本的全部退化数据。则Y=[Y1,Y2,…,YM]可表示所有样本总体退化数据。令ΔYm,i, j=Y(tm,i, j)-Y(tm,i, j-1),则ΔYm,i=[ΔYm,i,1,ΔYm,i,2,…,ΔYm,i,lm,i]T。

基于上述分析,若设备退化过程满足如式(5)所示带比例关系的非线性Wiener过程,则基于Wiener过程的独立增量特性,易知ΔYm,i服从多元正态分布N(μm,i,Σm,i),且其对应的期望和协方差矩阵分别为

μm,i=λm,iΔTm,i

(13)

(14)

式中,

ΔTm,i=[ΔTm,i,1,ΔTm,i,2,…,ΔTm,i,lm,i]T

(15)

(16)

(17)

(18)

由此,可得设备退化数据Y对应的轮廓对数似然函数:

(19)

(20)

(21)

(22)

将式(21)和式(22)代入式(20)可得

(23)

若令

(24)

(25)

(26)

(27)

将式(26)代入式(24)可得

(28)

4 退化数据折算

加速退化试验的最终目的是为了得到正常应力条件下设备的寿命、剩余寿命、退化模型参数等特征参量。因此,在预测常应力条件下设备的剩余寿命时,首先需要将加速应力条件下设备的退化数据折算为常应力条件下的退化数据。

针对式(3)所示非线性退化模型,基于加速因子不变原则,可得不同加速应力下加速因子[18-19]的定义式为

(29)

为将加速应力条件下设备的退化数据折算为常应力条件下的退化数据,给出引理1。

引理 1若应力Sp对应力Sq的加速因子为Ap,q,则应力Sp下的退化数据可折算为应力Sq条件下的退化数据,具体折算公式如下:

(30)

(31)

引理1的证明过程详见文献[16],在此不加赘述。

5 退化状态在线更新

基于目标设备的现场监测数据实时更新退化模型参数,可以突出个体退化过程的差异性,能够有效提升剩余寿命预测的准确性。本文基于卡尔曼滤波原理对设备退化状态进行在线更新。

(32)

(33)

(34)

(35)

L=[1,0]

(36)

进而可将式(32)转换为

(37)

为对式(37)进行卡尔曼滤波处理,首先需定义参数如下:

(1) 真实退化状态滤波均值

(38)

(2) 真实退化状态滤波方差

(39)

(3) 真实退化状态一步预测均值

(40)

(4) 真实退化状态一步预测方差

(41)

(5)Wj-1的协方差

(42)

其中，E(·)表示求期望值；D(·)表示求方差；cov(·)表示求协方差。

基于上述分析,利用卡尔曼滤波原理,即可实现对设备退化状态的在线更新,具体步骤如下:

步骤 1状态预测

(43)

步骤 2协方差预测

(44)

步骤 3滤波增益

(45)

步骤 4状态更新

(46)

步骤 5协方差更新

Pj|j=Pj|j-1-KjLPj|j-1

(47)

则利用式(33)～式(47),即可实现常应力条件下目标设备退化状态的在线更新。

6 剩余寿命自适应预测

设备的寿命被通常定义为从初始状态运行直至首次失效所经历的总时间,也称首达时间[22],其定义式为

T=inf{t:X(t)≥D|D>0}

(48)

式中,D为预先设置的设备失效阈值,通常基于工程经验确定。

基于寿命的定义,易得任意运行时刻tj对应的剩余寿命lj,定义式为

L=inf{lj:X(tj+lj)≥D|D>0}

(49)

文献[16]给出了非线性Wiener退化模型的剩余寿命概率密度函数:

(50)

(51)

基于卡尔曼滤波更新机制,可知参数α与设备真实退化量折算值X*满足二维正态分布,进而可得对应的条件分布为

(52)

(53)

利用全概率公式,将式(52)和式(53)代入式(51),即可得到常应力条件下设备剩余寿命的概率密度函数:

(54)

通常情况下,可采用剩余寿命预测的期望作为剩余寿命的预测值,具体计算过程为

7 实例分析

加速度计是重要的机载设备,是飞机惯导系统的核心组成部分,其性能状态的好坏将对飞机的战斗力产生直接影响。通过监测加速度计的性能状态,合理预测其剩余寿命,有助于科学制定维修策略,确保飞行安全与任务遂行。机载加速度计属于高可靠长寿命产品,常应力条件下退化极为缓慢,导致传统退化试验难以适用,因而需采用加速退化试验来高效经济地掌握其退化规律。本文采用文献[10]给出的某型加速度计恒定应力加速退化试验数据进行分析,具体试验条件如下:

(1) 加速度计退化指标选取一次项标度因数K1,且失效阈值为0.002 5g/V。

(2) 加速度计其退化过程满足非线性Wiener过程,且Λ(t|b)=tb;

(3) 加速度计对温度应力敏感,且其加速模型满足Arrhenius模型,即ρ(Sm|β)=exp(-β/S);

(4) 加速试验共有S1=65 ℃,S2=75 ℃,S3=85 ℃这3组应力水平,并且每个应力条件下设置6个试验样本。

具体退化实验数据如图3所示。

图3 加速退化数据Fig.3 Accelerated degradation data

由图3可知,仅有75 ℃条件下的第2个加速度计样本的性能退化量超过了失效阈值,表明该设备发生失效,而其他试验样本则均未发生失效。因此,本文选取该设备作为目标设备进行分析,以验证本文所提方法的准确性,其余加速度计退化试验数据则用于进行模型参数估计。

进一步,基于本文提出的两步极大似然估计法,即可对退化模型参数进行估计,具体参数估计结果如表1所示。由于fminsearch函数初值的选取对求解似然函数最大值的影响极大,因此本文采用随机生成初值,再带入目标函数进行对比分析的方法来确定迭代初值。

表1 退化模型参数估计值

为证明漂移系数与扩散系数比例关系的真实存在,且与加速应力大小无关,本文利用极大似然估计法对不同应力条件下退化模型的参数进行估计,具体结果如表2所示。

表2 不同应力条件下的漂移/扩散系数估计

通常情况下,加速度计正常工作的环境温度约为20 ℃。则由式(29)可得加速应力S2(75 ℃)对常应力S0(20 ℃)的加速因子为

(56)

进一步,基于本文给出的加速退化数据折算方法,可得目标设备在常应力条件下的退化数据折算值,具体如表3所示。

表3 退化数据折算

基于上述分析,利用本文提出的基于卡尔曼滤波的退化状态在线更新方法,即可在线更新常应力条件下设备的退化状态,如图4所示,进而实现对设备剩余寿命的在线预测。

图4 退化模型参数在线更新Fig.4 Online updating of degradation model parameter

为便于对比分析,将本文所提方法记为M1；将文献[16]提出的考虑漂移系数与扩散系数同步更新,但需满足特定共轭先验分布假设的方法记为M2；将文献[15]提出的仅更新漂移系数而不更新扩散系数的方法记为M3。不同方法对应剩余寿命预测结果如表4所示。

表4 剩余寿命预测结果

由表4可知,M1对应的剩余寿命预测值与真实剩余寿命最为接近,且其对应的95%置信区间可以实现在寿命周期全阶段对目标设备真实剩余寿命的全覆盖,而M2与M3则无法实现对目标设备剩余寿命的全覆盖,表明本文所提方法的准确性更高。进一步分析可以发现,M2需要满足特定共轭先验分布假设,而当假设条件不能成立时,将产生较大的预测误差,可能造成对剩余寿命的乐观估计,进而引起维修和更换的延后,将会对飞行安全与任务完成产生严重威胁。而M3未能同步更新漂移系数与扩散系数,将造成寿命预测值偏小,可能导致设备的提前维修和更换,增加了维修成本,降低了保障效能。此外,由表4可知,M3通过固定扩散系数获得了更低的预测不确定性,从而使得剩余寿命预测置信区间较窄,有助于提升预测精度,但该方法可能导致置信区间过窄以至于无法包含真实剩余寿命,导致预测准确性难以满足要求。而M1与M2均将漂移系数与扩散系数视为动态随机变量,从而增大了预测的不确定性,致使M1与M2的置信区间较M3更宽,但M2需要满足特定共轭先验分布假设才能成立,当该假设条件不满足时,剩余寿命预测的准确性将难以保证。

为了进一步比较不同方法的差异性,本文以均方误差(mean square error, MSE)值作为衡量剩余寿命预测准确性的标准。MSE值可以直观体现剩余寿命预测值偏离真实剩余寿命的情况,且MSE值越小,表明剩余寿命预测的准确性越高,反之则准确性越低。MSE的定义为