带测量误差的设备非线性退化建模与剩余寿命在线预测

Nonlinear Degradation Modeling and Online Prediction of Remaining Life for Equipment with Measurement Errors

摘要

现有的剩余寿命在线预测方法通常基于贝叶斯理论更新随机退化模型的漂移参数，但未更新扩散参数，为此，提出一种同时更新漂移与扩散参数的新方法。建立了考虑多种退化模式的随机退化模型，并依据首达时间原理推导出寿命及剩余寿命概率密度函数。先采用极大似然法离线估计模型的初始参数，再结合贝叶斯原理与期望最大化算法在线更新漂移参数与扩散参数。电容退化数据、陀螺仪漂移数据、铝合金构件裂纹增长数据验证了所提方法的有效性。

Abstract

The existing online prediction methods for remaining life typically updated the drift parameters of stochastic degradation models based on Bayesian theory， while did not update the diffusion parameters. So a new method was proposed to simultaneously update both drift parameters and diffusion parameters. A stochastic degradation model was established considering multiple degradation modes， and the probability density functions of lifetime and remaining life were derived based on the first-passage-time principle. The initial parameters of the model were estimated offline by maximum likelihood method. Subsequently， the drift parameters and diffusion parameters were updated online by integrating Bayesian theory and expectation maximization algorithm. The effectiveness of the proposed method was validated by capacitor degradation data， gyroscope drift data， and aluminum alloy components crack growth data.

Graphical abstract

关键词

退化建模 / 参数估计 / 参数更新 / 剩余寿命预测

Key words

degradation modeling / parameter estimation / parameter updating / remaining life prediction

0 引言

航空、航天、军工等对可靠性要求较高的领域需要大量结构复杂、寿命长的可靠设备。设备复杂度的提高加大了运行过程中的不确定性，对设备的安全稳定运行构成巨大挑战。对成本高、数量少、失效模式复杂的设备而言，剩余寿命预测是重要研究方向^［1-4］。传统的寿命预测方法高度依赖失效数据，而高可靠性、长寿命设备在常规环境下或加速寿命试验中难以获取充足的失效样本，严重制约了寿命预测方法的适用性。设备寿命与性能退化过程强关联，因此，基于性能退化数据的寿命预测方法为高可靠性、长寿命设备的寿命估计提供了新的技术路径^［5-6］。退化数据可全面表征设备全生命周期内的退化特征、个体差异、环境因素对设备寿命的影响规律，因此基于性能退化数据的寿命预测方法受到国内外学者的广泛关注。

实际工程中，设备的性能退化过程具有显著的随机性，而随机过程能更真实地反映设备性能退化的轨迹特征，因此当前研究多采用随机过程描述设备的动态演变规律。实际工程中，由于设备载荷的差异及内外环境的变化，退化数据往往表现出非单调性特征。Wiener过程作为一种随机过程，能有效描述具有递增或递减趋势的非单调退化行为，因而在退化建模中得到广泛应用。

针对具有线性特征的设备性能退化轨迹，KAHLE等^［7］提出基于线性Wiener过程的退化模型，并给出参数估计方法与设备寿命的定义。为进一步刻画同一批次设备因工艺、设计、功能差异导致的个体性能波动，厉海涛等^［8］建立考虑随机漂移参数的线性Wiener过程退化模型来预测动量轮寿命。SI等^［9］将测量误差作为模型不确定性来源，并基于贝叶斯推理与期望最大化（expectation-maximization，EM）算法推导出剩余寿命的概率密度函数。PENG等^［10］提出漂移因子服从斜正态分布的线性Wiener退化模型。GUAN等^［11］针对线性分段退化设备，提出了考虑测量误差的两阶段线性Wiener过程退化模型，并将其用于轴承寿命预测。线性Wiener模型仅适用于具有线性退化趋势的情形，而实际设备的退化多具有非线性特征，因此提出一系列非线性Wiener退化模型。SI等^［12］构建了指数型与幂函数型的非线性模型，并推导出寿命与剩余寿命的解析表达式。彭才华等^［13］融合幂函数型模型与线性模型，提出一种非线性退化模型来预测激光发生器的寿命。YU等^［14］建立了考虑三源不确定性的非线性Wiener过程退化模型，推导了剩余寿命的概率密度函数，并通过涡扇发动机数据验证了所提方法的有效性。LIN等^［15］针对设备不同时段的退化速率不一致的问题，提出两阶段非线性Wiener过程建模方法，提高了寿命预测的准确性。

随着模型结构的不断丰富与完善，研究者将目光投向如何利用实时监测数据动态更新模型参数，以提高剩余寿命预测的准确性。在此背景下，贝叶斯理论凭借其能自然融合先验知识与实时观测数据的核心优势而广泛用于模型参数的动态更新，获得更准确的剩余寿命估计值。杨升^［16］深入研究了基于贝叶斯理论的数控车床刀具剩余寿命在线预测。任宏宇等^［17］构建了基于长短期记忆神经网络的绝缘栅双极型晶体管（insulated gate bipolar transistor，IGBT）寿命预测模型，并引入贝叶斯理论优化网络参数。XU等^［18］针对温度变化对锂电池放电容量和老化的影响，提出融合Arrhenius 模型与Wiener过程的新型退化模型，并基于贝叶斯方法实现了锂电池剩余寿命的在线预测。FAN等^［19］针对存在竞争失效的设备，构建了基于序贯贝叶斯理论的剩余寿命预测模型。这些研究共同印证了贝叶斯理论在动态更新与寿命预测中的有效性与灵活性。

然而，现有基于退化数据的剩余寿命在线预测方法大多侧重于利用贝叶斯理论更新随机退化模型的漂移系数，往往忽视了对扩散系数的实时修正，更缺乏在同时考虑测量误差的情况下，系统研究漂移与扩散参数的协同更新。实际上，受设备个体差异与环境噪声等不确定性因素的影响，基于历史数据估计的扩散参数可能与在线监测数据的估算值存在偏差。为实现更可靠的参数估计，有必要利用在线监测信息动态更新扩散参数。为此，本文提出一种能同时考虑测量误差，并实现在线联合更新漂移参数与扩散参数的剩余寿命预测方法，以进一步提高预测的准确性与可靠性。

1 退化建模

大多数设备的性能退化数据具有非线性特征，因此本文采用非线性随机过程描述设备的退化过程，选用扩散过程刻画设备性能的退化现象。扩散过程是由标准布朗运动驱动的时间与状态都连续的随机过程，但扩散过程的漂移系数与扩散系数都不局限于常数。基于上述特性，本文建立的扩散模型为^［12］

\mathrm{d}X(t)=\mu (t,\theta )\mathrm{d}t+{\delta }_{\mathrm{B}}\mathrm{d}B(t)（1）

式中：X（t）表示退化过程；μ（t，θ）为漂移系数；θ为未知参数；δ_B为扩散系数；B（t）表示标准布朗运动。

设备的退化数据一般由带误差的仪器测得，而随机变化的外部环境影响测量值，导致真实的退化数据难以获得，但含有误差的测量数据依然能反映设备的退化趋势。因此本文建立带测量误差的退化模型：

\left.\begin{array}{l}X(t)=a\phi (t,b)+{\delta }_{\mathrm{B}}B(t)\\ Y(t)=X(t)+\epsilon \end{array}\right\}（2）

其中，aφ（t，b）为退化模型的漂移项，若φ（t，b）=t，则式（2）为线性退化模型；a为刻画同类设备个体间差异的漂移系数，且服从均值{\mu }_a、方差{\delta }_a^{\mathrm{2}}的正态分布；b为固定参数，描述同类设备的共性；ε为误差项，且服从均值0、方差{\delta }_{\epsilon }^{\mathrm{2}}的正态分布。为简化后续的寿命与剩余寿命函数的推导及参数估计，本文假设a、b、B（t）、ε相互独立。为充分验证所提方法的有效性与普适性，本文选用模型M₁即非线性退化模型φ₁（t，b）=t ^b +exp（bt）、模型M₂即指数型退化模型φ₂（t，b）=exp（bt）、模型M₃即线性退化模型φ₃（t，b）=t分别进行测试分析。

2 寿命预测

本文将设备的性能退化过程｛X（t），t≥0｝中首次穿越失效阈值ω的时间定义为设备的寿命^{［11，13-14］}，因此，寿命T可描述为

T=inf｛t：X（t）≥ω| X（0）<ω｝（3）

为不失一般性，令初始退化值X（t）=0，失效阈值ω>0。为求解寿命T的概率密度函数f_T （t）的解析表达式，先给出以下假设^［20］：

假设1 依据首达时间定义的设备寿命，如果设备在t时刻正常运行，则认为t时刻以前未发生失效。

假设2 如果退化过程｛X（t）\left|t\ge \mathrm{0}\}\right.在t时刻穿越失效阈值即X（t）=ω，那么｛X（t）\left|t\ge \mathrm{0}\}\right.在t时刻以前穿越失效阈值ω的概率可忽略不计。

为求解f_T （t）的表达式，先考虑退化模型中不含随机数即a为常数的情况。基于假设1、假设2、式（1）可得定理1：

定理1^［21］ 如果μ（t，θ）是t≥0上的连续函数，在给定参数a时，｛X（t）\left|t\ge \mathrm{0}\}\right.穿越失效阈值ω的首达时间的概率密度函数为

f_{T\left|a\right.}(t\left|a\right.)\cong \frac{S(t)+\mu (t,\theta )}{t\sqrt[]{\mathrm{2}\mathrm{\pi }t}}\mathrm{e}\mathrm{x}\mathrm{p}(-\frac{S^{\mathrm{2}}(t)}{\mathrm{2}t})（4）

S(t)=\frac{\mathrm{1}}{{\delta }_{\mathrm{B}}}(\omega -{\int }_{\mathrm{0}}^t\mu (\tau ,\theta )\mathrm{d}\tau )

a为随机参数时，可通过全概率公式：

f_T(t)={\int }_a^{}{f}_{T\left|a\right.}(t\left|a\right.)p(a)\mathrm{d}a=E_a[f_{T\left|a\right.}(t\left|a\right.)]（5）

式中：p（a）为参数a的概率密度函数；E_a ［*］表示对随机变量a求期望。

得到f_T （t）的表达式。当a为随机变量时，为简化求解寿命概率密度函数的过程，先给出引理1：

引理1^［12］ 若a\sim N(\mu ,{\delta }_{}^{\mathrm{2}})，A,B,C,\omega \in \mathbf{R}，则下式成立：

E_a[(\omega -Aa)\mathrm{e}\mathrm{x}\mathrm{p}(-\frac{{(\omega -Ba)}^{\mathrm{2}}}{\mathrm{2}C})]=\sqrt[]{\frac{C}{B^{\mathrm{2}}{\delta }^{\mathrm{2}}+C}}\cdot

(\omega -A\frac{B{\delta }^{\mathrm{2}}\omega +\mu C}{B^{\mathrm{2}}{\delta }^{\mathrm{2}}+C})\mathrm{e}\mathrm{x}\mathrm{p}(-\frac{{(\omega -B\mu )}^{\mathrm{2}}}{\mathrm{2}(B^{\mathrm{2}}{\delta }^{\mathrm{2}}+C)})（6）

通过定理1、全概率公式（式（5））、引理1可推出模型M₁的寿命概率密度函数f_T （t），如定理2所示。

定理2 对于退化模型M₁，a\sim N({\mu }_a,{\delta }_a^{\mathrm{2}})时，｛X（t）\left|t\ge \mathrm{0}\}\right.穿越失效阈值ω的首达时间寿命概率密度函数为

f_T\left(t\right)\cong \frac{\mathrm{1}}{t\sqrt[]{\mathrm{2}\mathrm{\pi }({\delta }_a^{\mathrm{2}}{\beta }^{\mathrm{2}}(t)+{\delta }_{\mathrm{B}}^{\mathrm{2}}t)}}\mathrm{e}\mathrm{x}\mathrm{p}(-\frac{(\omega -{\mu }_a\beta {(t))}^{\mathrm{2}}}{\mathrm{2}({\delta }_a^{\mathrm{2}}{\beta }^{\mathrm{2}}(t)+{\delta }_{\mathrm{B}}^{\mathrm{2}}t)})\cdot[\omega -(\beta (t)-b{t}^b-bt\mathrm{e}\mathrm{x}\mathrm{p}(bt))\frac{{\delta }_a^{\mathrm{2}}\omega \beta (t)+{\mu }_a{\delta }_{\mathrm{B}}^{\mathrm{2}}t}{{\delta }_a^{\mathrm{2}}{\beta }^{\mathrm{2}}(t)+{\delta }_{\mathrm{B}}^{\mathrm{2}}t})]（7）

\beta \left(t\right)=t^b+\mathrm{e}\mathrm{x}\mathrm{p}(bt)-\mathrm{1}

在t_i 时刻，根据首达时间的定义得到设备的剩余寿命

L_{t_i}=\mathrm{i}\mathrm{n}\mathrm{f}\{l_{t_i}>\mathrm{0}:X(l_{t_i}+t_i)\ge \omega \}（8）

同样通过定理1、全概率公式（式（5））、引理1可推出模型M₁的剩余寿命概率密度函数：

\begin{array}{c}f(l_{t_i})\cong ({\omega }_{t_i}-\eta (l_{t_i})\frac{{\delta }_a^{\mathrm{2}}{\omega }_{t_i}\gamma (l_{t_i})+{\mu }_a{l}_{t_i}{\delta }_{\mathrm{B}}^{\mathrm{2}}}{{\delta }_a^{\mathrm{2}}{\gamma }^{\mathrm{2}}(l_{t_i})+l_{t_i}{\delta }_{\mathrm{B}}^{\mathrm{2}}})\cdot \\ \mathrm{e}\mathrm{x}\mathrm{p}(-\frac{({\omega }_{t_i}-{\mu }_a\gamma (l_{t_i}{))}^{\mathrm{2}}}{\mathrm{2}({\delta }_a^{\mathrm{2}}{\gamma }^{\mathrm{2}}(l_{t_i})+l_{t_i}{\delta }_{\mathrm{B}}^{\mathrm{2}})})(\mathrm{2}\mathrm{\pi }{l}_{t_i}^{\mathrm{2}}\{[(l_{t_i}+t_i{)}^b+\end{array}

\mathrm{e}\mathrm{x}\mathrm{p}(b(l_{t_i}+t_i))-t_i^b-\mathrm{e}\mathrm{x}\mathrm{p}{(bt)]}^{\mathrm{2}}+l_{t_i}{\delta }_B^{\mathrm{2}}{\})}^{-\mathrm{1}/\mathrm{2}}（9）

\gamma (l_{t_i})=(l_{t_i}+t_i{)}^b+\mathrm{e}\mathrm{x}\mathrm{p}(b(l_{t_i}+t_i))-(t_i^b+\mathrm{e}\mathrm{x}\mathrm{p}(b{l}_{t_i}))

\begin{array}{l}\eta (l_{t_i})=(t_i+l_{t_i}{)}^b+(\mathrm{1}-b{l}_{t_i})\mathrm{e}\mathrm{x}\mathrm{p}(b(t_i+l_{t_i}))-\\ (t_i^b+\mathrm{e}\mathrm{x}\mathrm{p}(b{t}_i))-b{l}_{t_i}[b(t_i+l_{t_i}{)]}^{b-\mathrm{1}}\end{array}

{\omega }_{t_i}=\omega -X\left(t_i\right)

式中：X（t_i ）为t_i 时刻的退化值。

模型M₂、M₃的剩余寿命概率密度函数可参考文献［8］、［12］，不再赘述。

3 参数估计

本文采用极大似然法估计退化模型中的未知参数^［13］。假设Y_h （t_j ）为设备h在第j时刻测量得到的性能退化值，则退化轨迹可表示为

Y_h(t_j)=a_h\phi (t_{h,j},b)+{\delta }_{\mathrm{B}}B(t_{h,j})+{\epsilon }_{h,j}（10）

h=1，2，…，n j=1，2，…，m_h

式中：n为退化设备的总数；m_h 为第h个设备的测量数据量；a_h 独立同分布且满足a_h\sim N({\mu }_a,{\delta }_a^{\mathrm{2}})。

令{\mathit{Y}}_h=(Y_h(t_{\mathrm{1}}),Y_h(t_{\mathrm{2}}),\cdots ,Y_h(t_m{))}^{\mathrm{T}}表示n个设备退化数据，{\mathit{T}}_h=(T_{h,\mathrm{1}},T_{h,\mathrm{2}},\cdots ,T_{h.m_h}{)}^{\mathrm{T}}为时间向量，其中，T_{h.j}=\phi \left(t_{h,j},b\right)。根据式 （10） 和独立性假设可知Y_h 服从多变量正态分布，其均值 μ₀和协方差 Σ_h 分别为

\left.\begin{array}{l}{\mathit{\mu }}_{\mathrm{0}}={\mu }_a{\mathit{T}}_h\\ {\mathit{\Sigma }}_h={\delta }_a^{\mathrm{2}}{\mathit{T}}_h{\mathit{T}}_h^{\mathrm{T}}+{\mathit{\Omega }}_h\end{array}\right\}（11）

{\mathit{\Omega }}_h={\delta }_{\mathrm{B}}^{\mathrm{2}}{\mathit{Q}}_h+{\delta }_{\epsilon }^{\mathrm{2}}\mathit{F}（12）

\left.\begin{array}{l}{\mathit{Q}}_h=\left[\begin{array}{ccccc}{t}_{h,\mathrm{1}}& t_{h,\mathrm{1}}& t_{h,\mathrm{1}}& \cdots & t_{h,\mathrm{1}}\\ t_{h,\mathrm{1}}& t_{h,\mathrm{2}}& t_{h,\mathrm{2}}& \cdots & t_{h,\mathrm{2}}\\ t_{h,\mathrm{1}}& t_{h,\mathrm{2}}& t_{h,\mathrm{3}}& \cdots & t_{h,\mathrm{3}}\\ ⋮& ⋮& ⋮& & ⋮\\ t_{h,\mathrm{1}}& t_{h,\mathrm{2}}& t_{h,\mathrm{3}}& \cdots & t_{h,m}\end{array}\right]\\ \mathit{F}=[F_{ij}{]}_{m_h\times m_h}\end{array}\right\}（13）

F_{ij}=\left\{\begin{array}{cc}\mathrm{1}& i=j=\mathrm{1}\\ \mathrm{2}& i=j=\mathrm{2,3},\cdots ,m_h\\ -\mathrm{1}& i=j+\mathrm{1}\mathrm{或}j=i+\mathrm{1}\end{array}\right.

采用极大似然法估计退化模型中的未知参数μ_a 、δ_a 、δ_B、δ_ε 、b。由上述分析可得Y对应的对数似然函数为

l\left(Y\right)=-\mathrm{l}\mathrm{n}\mathrm{2}\mathrm{\pi }{\displaystyle \sum _{h=\mathrm{1}}^n}{m}_h-\frac{\mathrm{1}}{\mathrm{2}}{\displaystyle \sum _{h=\mathrm{1}}^n}[\mathrm{l}\mathrm{n}\left|{\Sigma }_h\right|+({\mathit{Y}}_h-

{\mu }_a{\mathit{T}}_h{)}^{\mathrm{T}}{\mathit{\Sigma }}^{-\mathrm{1}}({\mathit{Y}}_h-{\mu }_a{\mathit{T}}_h)]（14）

\left|{\mathit{\Sigma }}_h\right|=\left|{\mathit{\Omega }}_h\right|(\mathrm{1}+{\delta }_a^{\mathrm{2}}{\mathit{T}}_h^{\mathrm{T}}{\mathit{\Omega }}_h^{-\mathrm{1}}{\mathit{T}}_h)（15）

{\mathit{\Sigma }}_h^{-\mathrm{1}}={\mathit{\Omega }}_h^{-\mathrm{1}}-\frac{{\delta }_a^{\mathrm{2}}}{\mathrm{1}+{\delta }_a^{\mathrm{2}}{\mathit{T}}_h^{\mathrm{T}}{\mathit{\Omega }}_h^{-\mathrm{1}}{\mathit{T}}_h}{\mathit{\Omega }}_h^{-\mathrm{1}}{\mathit{T}}_h{\mathit{T}}_h^{\mathrm{T}}{\mathit{\Omega }}_h^{-\mathrm{1}}（16）

对式（14）求μ_a 的一阶偏导，并令其等于零，可得关于μ_a 的极大似然估计：

{\widehat{\mu }}_a={\displaystyle \sum _{h=\mathrm{1}}^n}{\mathit{T}}_h^{\mathrm{T}}{\mathit{\Sigma }}_h^{-\mathrm{1}}{\mathit{Y}}_h/({\displaystyle \sum _{h=\mathrm{1}}^n}{\mathit{T}}_h^{\mathrm{T}}{\mathit{\Sigma }}_h^{-\mathrm{1}}{\mathit{T}}_h)（17）

将式（17）代入式（14）可得关于参数δ_a 、δ_B、δ_ε 、b的剖面对数似然函数：

\begin{array}{c}l({\delta }_a,{\delta }_{\mathrm{B}},b,{\delta }_{\epsilon }\left|Y,{\widehat{\mu }}_a\right.)=-\mathrm{l}\mathrm{n}(\mathrm{2}\mathrm{\pi }){\displaystyle \sum _{h=\mathrm{1}}^n}{m}_h-\frac{\mathrm{1}}{\mathrm{2}}{\displaystyle \sum _{h=\mathrm{1}}^n}ln\left|{\mathit{\Sigma }}_h\right|-\\ \frac{\mathrm{1}}{\mathrm{2}}[{\displaystyle \sum _{h=\mathrm{1}}^n}{\mathit{Y}}_h^{\mathrm{T}}{\mathit{\Omega }}_h^{-\mathrm{1}}{\mathit{Y}}_h-{({\displaystyle \sum _{h=\mathrm{1}}^n}{\mathit{T}}_h^{\mathrm{T}}{\mathit{\Sigma }}_h^{-\mathrm{1}}{\mathit{Y}}_h)}^{\mathrm{2}}(\mathrm{2}{\mathit{T}}_h^{\mathrm{T}}{\mathit{\Sigma }}_h^{-\mathrm{1}}{\mathit{T}}_h-\end{array}{\displaystyle \sum _{h=\mathrm{1}}^n}{\mathit{T}}_h^{\mathrm{T}}{\mathit{\Omega }}_h^{-\mathrm{1}}{\mathit{T}}_h)/({\mathit{T}}_h^{\mathrm{T}}{\mathit{\Sigma }}_h^{-\mathrm{1}}{\mathit{T}}_h{)}^{\mathrm{2}}]（18）

通过极大化式（18）便可得出δ_a 、δ_B、δ_ε 、b的估计值，再将δ_a 、δ_B、δ_ε 、b的估计值代入式 （17）即可得出μ_a 的估计值。

4 参数在线更新

令 Y_1：k =（y₁，y₂，…，y_k ）为 t₁到t_k 时刻设备的历史状态退化数据，则t_k 时刻的漂移系数a的后验概率为

p(a\left|Y_{\mathrm{1}:k}\right.)=\frac{p(Y_{\mathrm{1}:k}\left|a\right.)p(a)}{p(Y_{\mathrm{1}:k}\left|a\right.)}\propto p(Y_{\mathrm{1}:k}\left|a\right.)（19）

式中：p（a）为漂移系数a的先验概率密度函数；p（ Y_1：k |a）为在给定参数a后的联合概率密度函数。

前文已假设参数a服从正态分布，因此有

p(a)=\frac{\mathrm{1}}{\sqrt[]{\mathrm{2}\mathrm{\pi }{\delta }_a^{\mathrm{2}}}}\mathrm{e}\mathrm{x}\mathrm{p}(-\frac{a-{\mu }_a}{\mathrm{2}{\delta }_a^{\mathrm{2}}})（20）

由Wiener过程的性质可知，模型M₁的联合概率密度函数为

p({\mathit{Y}}_{\mathrm{1}:k}\left|a\right.)={\displaystyle \prod _{i=\mathrm{1}}^k}\frac{\mathrm{1}}{\sqrt[]{\mathrm{2}\mathrm{\pi }(t_i-t_{i-\mathrm{1}})}}\mathrm{e}\mathrm{x}\mathrm{p}(-\frac{Y(t_i)-Y(t_{i-\mathrm{1}})}{\mathrm{2}{\delta }_B^{\mathrm{2}}(t_i-t_{i-\mathrm{1}})}-

\frac{a(t_i^b+\mathrm{e}\mathrm{x}\mathrm{p}(b{t}_i)-t_{i-\mathrm{1}}^b-\mathrm{e}\mathrm{x}\mathrm{p}(b{t}_{i-\mathrm{1}}{))}^{\mathrm{2}}}{\mathrm{2}{\delta }_{\mathrm{B}}^{\mathrm{2}}(t_i-t_{i-\mathrm{1}})})（21）

因此，可得

p(a\left|Y_{\mathrm{1}:k}\right.)\propto p(Y_{\mathrm{1}:k}\left|a\right.)p(a)\propto \mathrm{e}\mathrm{x}\mathrm{p}(-\frac{(a-{\mu }_a{)}^{\mathrm{2}}}{\mathrm{2}{\delta }_a^{\mathrm{2}}}-{\displaystyle \sum _{i=\mathrm{1}}^k}[-\frac{a(t_i^b+\mathrm{e}\mathrm{x}\mathrm{p}(b{t}_i)-t_{i-\mathrm{1}}^b-\mathrm{e}\mathrm{x}\mathrm{p}(b{t}_{i-\mathrm{1}}{))}^{\mathrm{2}}}{\mathrm{2}{\delta }_{\mathrm{B}}^{\mathrm{2}}(t_i-t_{i-\mathrm{1}})}+\frac{Y(t_i)-Y(t_{i-\mathrm{1}})}{\mathrm{2}{\delta }_{\mathrm{B}}^{\mathrm{2}}(t_i-t_{i-\mathrm{1}})})]（22）

根据共轭分布性质可知，如果参数a的先验分布是正态分布，则参数a的后验分布也是正态分布，即a\left|Y_{\mathrm{1}:k}\right.~N({\mu }_{a,k},{\delta }_{a,k}^{\mathrm{2}})。整理式（22）可得

{\mu }_{a,k}=\frac{{\delta }_{\mathrm{B}}^{\mathrm{2}}{\mu }_{a,\mathrm{0}}+\mathrm{\Delta }Y{\displaystyle \sum _{i=\mathrm{1}}^k}\frac{t_i^b+\mathrm{e}\mathrm{x}\mathrm{p}(b{t}_i)-t_{i-\mathrm{1}}^b-\mathrm{e}\mathrm{x}\mathrm{p}(b{t}_{i-\mathrm{1}})}{t_i-t_{i-\mathrm{1}}}}{{\delta }_{\mathrm{B}}^{\mathrm{2}}+{\delta }_{a,\mathrm{0}}^{\mathrm{2}}{\displaystyle \sum _{i=\mathrm{1}}^k}\frac{(t_i^b+\mathrm{e}\mathrm{x}\mathrm{p}(b{t}_i)-t_{i-\mathrm{1}}^b-\mathrm{e}\mathrm{x}\mathrm{p}(b{t}_{i-\mathrm{1}}{))}^{\mathrm{2}}}{t_i-t_{i-\mathrm{1}}}}（23）

\mathrm{\Delta }Y=Y\left(t_i\right)-Y\left(t_{i-\mathrm{1}}\right)

{\delta }_{a,k}^{\mathrm{2}}=\frac{{\delta }_{\mathrm{B}}^{\mathrm{2}}{\delta }_{a,\mathrm{0}}^{\mathrm{2}}}{{\delta }_{\mathrm{B}}^{\mathrm{2}}+{\delta }_{a,\mathrm{0}}^{\mathrm{2}}{\displaystyle \sum _{i=\mathrm{1}}^k}\frac{t_i^b+\mathrm{e}\mathrm{x}\mathrm{p}(b{t}_i)-t_{i-\mathrm{1}}^b-\mathrm{e}\mathrm{x}\mathrm{p}(b{t}_{i-\mathrm{1}}{)}^{\mathrm{2}}}{t_i-t_{i-\mathrm{1}}}}（24）

模型M₂与M₃的漂移参数更新方法与此类似，故不再赘述。式（18）的扩散参数δ_B通过历史退化数据估算得出。设备个体差异和环境噪声等不确定性因素导致扩散参数δ_B随时间发生变化。为准确估算设备的剩余寿命，有必要更新扩散参数δ_B。本文构建的退化模型中，设备个体差异与环境噪声等不确定性因素在数学上可表征为参数估计过程中的隐变量，故扩散参数δ_B的在线更新在本质上是含有隐变量的极大似然估计问题。EM算法通过“期望步（E步）-最大化步（M步）”的迭代框架高效求解此类问题，因此，本文采用EM算法实现δ_B的递推更新，以确保参数估计结果始终精准贴合设备当前的实际运行状态。为便于表达，将t_k 时刻{\delta }_{\mathrm{B}}^{\mathrm{2}}的估计值记为{\delta }_{\mathrm{B},k}^{\mathrm{2}}。对于模型M₁，已知 Y_1：k 的对数似然函数

\mathrm{l}\mathrm{n}p({\mathit{Y}}_{\mathrm{1}:k},a|{\delta }_{\mathrm{B}}^{\mathrm{2}})=\mathrm{l}\mathrm{n}p({\mathit{Y}}_{\mathrm{1}:k}|{\delta }_{\mathrm{B}}^{\mathrm{2}},a)+\mathrm{l}\mathrm{n}p(a|{\delta }_{\mathrm{B}}^{\mathrm{2}})=

-\frac{k+\mathrm{1}}{\mathrm{2}}\mathrm{l}\mathrm{n}\mathrm{2}\mathrm{\pi }-\frac{\mathrm{1}}{\mathrm{2}}{\displaystyle \sum _{i=\mathrm{1}}^k}ln(t_i-t_{i-\mathrm{1}})-\frac{k+\mathrm{1}}{\mathrm{2}}\mathrm{l}\mathrm{n}{\delta }_{\mathrm{B}}{\delta }_a-\frac{(a-{\mu }_a{)}^{\mathrm{2}}}{\mathrm{2}{\delta }_a^{\mathrm{2}}}-{\displaystyle \sum _{i=\mathrm{1}}^k}[\frac{(Y(t_i)-Y(t_{i-\mathrm{1}})}{\mathrm{2}{\delta }_{\mathrm{B}}^{\mathrm{2}}(t_i-t_{i-\mathrm{1}})}-

\frac{a(t_i^b+\mathrm{e}\mathrm{x}\mathrm{p}(b{t}_i)-t_{i-\mathrm{1}}^b-\mathrm{e}\mathrm{x}\mathrm{p}(b{t}_{i-\mathrm{1}}{)))}^{\mathrm{2}}}{\mathrm{2}{\delta }_{\mathrm{B}}^{\mathrm{2}}(t_i-t_{i-\mathrm{1}})}]（25）

再求\mathrm{l}\mathrm{n}p({\mathit{Y}}_{\mathrm{1}:k},a|{\delta }_{\mathrm{B}}^{\mathrm{2}})的期望，可得

\begin{array}{l}{E}_{a|Y_{\mathrm{1}:k},{\delta }_{\mathrm{B}}^{\mathrm{2}}}[\mathrm{l}\mathrm{n}p(Y_{\mathrm{1}:k},a|{\delta }_{\mathrm{B}}^{\mathrm{2}})]=-\frac{k+\mathrm{1}}{\mathrm{2}}\mathrm{l}\mathrm{n}\mathrm{2}\mathrm{\pi }-k\mathrm{l}\mathrm{n}{\delta }_{\mathrm{B}}-\\ \mathrm{l}\mathrm{n}{\delta }_a-\frac{({\mu }_{a,k}^{}-{\mu }_a^{}{)}^{\mathrm{2}}+{\delta }_{a,k}^{\mathrm{2}}}{\mathrm{2}{\delta }_a^{\mathrm{2}}}-\frac{\mathrm{1}}{\mathrm{2}}{\displaystyle \sum _{i=\mathrm{1}}^k}ln(t_i-t_{i-\mathrm{1}})-\\ {\displaystyle \sum _{i=\mathrm{1}}^k}\frac{(Y(t_i)-Y(t_{i-\mathrm{1}}{))}^{\mathrm{2}}}{\mathrm{2}{\delta }_{\mathrm{B}}^{\mathrm{2}}(t_i-t_{i-\mathrm{1}})}+{\mu }_{a,k}{\displaystyle \sum _{i=\mathrm{1}}^k}[\frac{Y(t_i)-Y(t_{i-\mathrm{1}})}{{\delta }_{\mathrm{B}}^{\mathrm{2}}(t_i-t_{i-\mathrm{1}})}·\\ (t_i^b-t_{i-\mathrm{1}}^b+\mathrm{e}\mathrm{x}\mathrm{p}(b{t}_i)-\mathrm{e}\mathrm{x}\mathrm{p}(b{t}_{i-\mathrm{1}}))]+\end{array}

({\mu }_{a,k}^{\mathrm{2}}+{\delta }_{a,k}^{\mathrm{2}}){\displaystyle \sum _{i=\mathrm{1}}^k}\frac{(t_i^b-t_{i-\mathrm{1}}^b+\mathrm{e}\mathrm{x}\mathrm{p}(b{t}_i)-\mathrm{e}\mathrm{x}\mathrm{p}(b{t}_{i-\mathrm{1}}{))}^{\mathrm{2}}}{\mathrm{2}{\delta }_{\mathrm{B}}^{\mathrm{2}}(t_i-t_{i-\mathrm{1}})}（26）

对式（26）关于δ{}_{\mathrm{B}}^{\mathrm{2}}求偏导并令偏导数为零可得δ{}_{\mathrm{B}}^{\mathrm{2}}的更新表达式：

\begin{array}{c}{\delta }_{\mathrm{B},k}^{\mathrm{2}}=\frac{{\mu }_{a,k}^{\mathrm{2}}+{\delta }_{a,k}^{\mathrm{2}}}{k}{\displaystyle \sum _{i=\mathrm{1}}^k}\frac{(t_i^b-t_{i-\mathrm{1}}^b+\mathrm{e}\mathrm{x}\mathrm{p}(b{t}_i)-\mathrm{e}\mathrm{x}\mathrm{p}(b{t}_{i-\mathrm{1}}{))}^{\mathrm{2}}}{t_i-t_{i-\mathrm{1}}}+\\ \frac{\mathrm{1}}{k}\{{\displaystyle \sum _{i=\mathrm{1}}^k}\frac{(Y(t_i)-Y(t_{i-\mathrm{1}}{))}^{\mathrm{2}}}{t_i-t_{i-\mathrm{1}}}-\mathrm{2}{\mu }_{a,k}{\displaystyle \sum _{i=\mathrm{1}}^k}[\frac{Y(t_i)-Y(t_{i-\mathrm{1}})}{t_i-t_{i-\mathrm{1}}}·\end{array}(t_i^b-t_{i-\mathrm{1}}^b+\mathrm{e}\mathrm{x}\mathrm{p}(b{t}_i)-\mathrm{e}\mathrm{x}\mathrm{p}(b{t}_{i-\mathrm{1}}))]\}（27）

模型M₂与M₃的扩散参数更新方法与此类似，不再赘述。

5 实验研究

本节通过测量得到的设备性能退化数据验证提出的方法。为直观展示所提方法的预测效果，计算出所提方法预测的总体均方误差E_TMSE、均方根误差E_RMSE、平均绝对误差E_MAE^{［15，22］}。E_TMSE可通过均方误差E_MSE计算得出。t_k 时刻的均方误差为

E_{\mathrm{M}\mathrm{S}\mathrm{E},k}={\int }_{\mathrm{0}}^{\infty }(l_k-{\tilde{l}}_k{)}^{\mathrm{2}}{f}_{L_k\left|X_{\mathrm{1}:k}\right.}(L_k\left|X_{\mathrm{1}:k}\right.)d{l}_k（28）

式中：l_k、{\tilde{l}}_k分别为t_k 时刻设备剩余寿命的实际值与预测值；f_{L_k\left|X_{\mathrm{1}:k}\right.}(L_k\left|X_{\mathrm{1}:k}\right.)为设备剩余寿命的概率密度函数。

由此，可计算出

E_{\mathrm{T}\mathrm{M}\mathrm{S}\mathrm{E}}=\frac{\mathrm{1}}{n}{\displaystyle \sum _{k=\mathrm{1}}^n}{E}_{\mathrm{M}\mathrm{S}\mathrm{E},k}（29）

E_{\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}=\sqrt[]{\frac{\mathrm{1}}{n}{\displaystyle \sum _{k=\mathrm{1}}^n}(l_k-{\tilde{l}}_k{)}^{\mathrm{2}}}（30）

E_{\mathrm{M}\mathrm{A}\mathrm{E}}=\frac{\mathrm{1}}{n}{\displaystyle \sum _{k=\mathrm{1}}^n}|l_k-{\tilde{l}}_k|（31）

式中：n为设备剩余寿命预测的采样点数。

E_MAE越小，预测值越接近真实值，模型的预测精度越高。为验证所提方法的有效性与优势，设置退化模型无测量误差且无参数更新（方法Ⅰ）、退化模型含测量误差但无参数更新（方法Ⅱ）、退化模型含测量误差且仅更新漂移参数（方法Ⅲ）、退化模型含测量误差且同时更新漂移参数与扩散参数（方法Ⅳ）进行仿真分析。

5.1 高压脉冲电容退化数据

将某型号的5个高压脉冲电容避光贮存在有干燥剂的干燥皿内，跟踪监测其电容量的相对变化量，如图1所示^［23］，可以看出，电容量的相对变化量呈现明显的非线性特征，适用于非线性退化建模。

将电容的退化数据代入模型M₁，采用极大似然法估计退化模型参数的初始值，再通过贝叶斯原理与EM算法更新退化模型的漂移参数和扩散参数，将更新后的参数代入剩余寿命概率密度函数，得到电容的剩余寿命分布，其概率密度曲线如图2所示，剩余寿命预测结果如图3所示。

由图2可以看出，4种方法都能预测电容的剩余寿命，但预测结果不同。方法Ⅰ的预测值与实际值相差较大，方法Ⅱ～Ⅳ都能较准确地预测电容的剩余寿命。图2e所示为4种预测方法在最后一个（第14个）测量点的剩余寿命概率密度函数曲线，方法Ⅳ预测的剩余寿命概率密度函数曲线比方法Ⅰ、Ⅱ及Ⅲ的曲线峰值更大、形态更窄，这表明方法Ⅳ的预测结果更集中，其预测不确定性更小，即本文方法（方法Ⅳ）的预测结果更精确。减小剩余寿命预测的不确定性可降低设备的维护成本、提高决策结果的置信度^［21］，对设备的健康管理非常重要。由图2可知，方法Ⅱ对电容剩余寿命的预测精度高于方法Ⅰ。这说明测量误差是预测设备剩余寿命时必须考虑的因素。由图3可知，本文所提方法（方法Ⅳ）预测的电容剩余寿命更接近于实际的剩余寿命。为直观说明4种方法预测的差异，绘制4种方法预测的均方误差曲线（图4），计算4种方法预测的E_TMSE、E_RMSE、E_MAE（表1）。由图4可看出，方法Ⅰ的曲线最高、均方误差E_MSE最大，方法Ⅳ的曲线最低、均方误差最小。对比方法Ⅰ与方法Ⅱ的曲线可知，考虑测量误差后，预测的均方误差减小。由表1可知，方法Ⅰ的E_TMSE、E_RMSE、E_MAE最大，方法Ⅳ的E_TMSE、E_RMSE、E_MAE最小，说明提出的更新扩散参数方法更优。

将电容退化数据代入模型M₂中，采用极大似然法估计出模型参数后，通过贝叶斯原理与EM算法更新退化模型中的漂移参数和扩散参数，将更新后的参数代入剩余寿命概率密度函数，并绘制出其曲线，如图5所示。随着测量点的增多，剩余寿命的概率密度函数曲线的峰值增大、分布宽度变小，预测不确定性减小。方法Ⅰ预测的剩余寿命曲线与实际剩余寿命曲线相距较远，因此预测的电容剩余寿命误差较大。方法Ⅱ～Ⅳ预测的电容剩余寿命误差较小，预测的不确定性更小。由图5e可知，方法Ⅳ预测的电容剩余寿命概率密度函数曲线峰值最大、分布宽度最小，预测的不确定性最小。对比图5a、图5b可知，考虑测量误差既减小了电容剩余寿命的预测误差，又减小了剩余寿命预测的不确定性，说明考虑测量误差在预测设备剩余寿命的重要性。

由图6可知，方法Ⅳ预测的剩余寿命更接近实际剩余寿命。图7给出了4种方法的预测均方误差曲线，表2给出了4种方法预测的E_TMSE、E_RMSE、E_MAE。由图7可以看出，方法Ⅰ的均方误差最大，方法Ⅳ的均方误差最小。方法Ⅱ的均方误差小于方法Ⅰ，说明考虑测量误差后，预测的均方误差减小。由表2可知，方法Ⅱ的E_TMSE、E_RMSE、E_MAE小于方法Ⅰ，再次说明考虑测量误差的重要性。方法Ⅲ、Ⅳ的E_TMSE、E_RMSE、E_MAE小于方法Ⅱ，说明利用监测得到的退化数据实时更新退化模型参数提高了参数的估算精度，减小了剩余寿命的预测误差。方法Ⅳ的E_TMSE、E_RMSE、E_MAE最小，误差最小，预测结果最准确，再次说明更新退化模型扩散系数方法的优越性。对比表1、表2可知，模型M₁预测的电容剩余寿命误差小于模型M₂，说明模型M₁具有实用性。

将电容退化数据代入模型M₃，沿用前文所述的参数估计与更新方法，得到剩余寿命概率密度函数曲线（图8）与剩余寿命曲线（图9）。

对比图2、图5、图8可知，模型M₃预测的电容剩余寿命概率密度函数曲线的峰值最小，曲线形态最为平坦。该特征在图2e、图5e、图8e中表现得尤为明显，这说明模型M₃预测的电容剩余寿命不确定性更大。出现上述现象的原因是电容退化过程具有较强的非线性动态特征，而线性模型难以充分捕捉此类特征，导致预测结果分布较为发散。因此，进行电容剩余寿命预测时，优先选用非线性模型。由图8e可知，基于模型M₃的4种预测方法中，方法Ⅰ预测的剩余寿命概率密度函数曲线峰值小、分布宽，表明其预测的不确定性大；方法Ⅳ预测的剩余寿命概率密度函数曲线峰值高、分布窄，说明其预测的剩余寿命不确定性小。图9所示结果表明，方法Ⅳ所预测的剩余寿命最接近设备的真实剩余寿命。4种方法预测的均方误差曲线如图10所示。4种方法的E_TMSE、E_RMSE、E_MAE如表3所示。

由图10可以看出，方法Ⅲ与方法Ⅳ的均方误差随测量点的增多逐渐减小，方法Ⅳ的均方误差小于方法Ⅲ。由表3可知，方法Ⅰ的E_TMSE、E_RMSE、E_MAE最大，方法Ⅳ的E_TMSE、E_RMSE、E_MAE最小，这说明提出的参数更新方法也适用线性退化模型。对比表1、表2与表3可知，模型M₃的E_TMSE、E_RMSE、E_MAE均大于模型M₁与M₂，说明M₃预测的误差更大，这表明M₃与实际退化数据特征之间存在明显失配。电容退化过程是非线性的，线性模型难以准确刻画其动态演变，从而导致预测不确定性增大、精度下降，因此，应优先选用非线性模型预测电容的剩余寿命。

5.2 陀螺仪退化数据

受工作环境、电机转子磨损等影响，陀螺仪的漂移系数在使用过程中逐渐增大，表现为性能的渐进式退化。退化积累至失效阈值时，陀螺仪将无法维持正常工作，即发生失效。本文基于某陀螺仪的退化试验数据开展剩余寿命预测的研究。试验样本为同型号、同批次的6台陀螺仪。在相同的试验环境下，每隔100h采集陀螺仪的漂移系数，所得数据如图11所示^［24］。

将陀螺仪退化数据代入非线性退化模型M₁，依然先采用极大似然法估计模型参数的初始值，再采用贝叶斯方法与EM算法更新模型的漂移参数和扩散参数，最后将更新的参数代入剩余寿命概率密度函数，画出剩余寿命概率密度函数曲线，如图12所示；同时，其剩余寿命预测结果如图13所示。由图12可看出，4种方法预测的剩余寿命概率密度函数曲线差异较小，都能比较准确地预测陀螺仪的剩余寿命。由图13可看出，方法Ⅳ预测的陀螺仪剩余寿命更接近实际值。4种方法的均方误差曲线见图14，E_TMSE、E_RMSE、E_MAE如表4所示。由图14可以看出，方法Ⅰ、方法Ⅱ的均方误差曲线几乎重合，说明陀螺仪退化数据包含的测量误差较小；方法Ⅳ的均方误差曲线最低，误差最小。由表4可知，方法Ⅰ的E_TMSE、E_RMSE、E_MAE最大，方法Ⅳ的E_TMSE、E_RMSE、E_MAE最小，这说明方法Ⅰ的预测误差最大，方法Ⅳ的预测误差最小、效果最好，再次说明所提方法的有效性。

将陀螺仪退化数据代入模型M₂，同样进行参数的估计与更新，再将更新后的参数代入剩余寿命概率密度函数并画出函数曲线，如图15所示。由图15可看出，4种方法都能预测陀螺仪的剩余寿命。4种方法预测的剩余寿命概率密度函数曲线差异较小，因此通过分析均方误差曲线（图17）以及E_TMSE、E_RMSE、E_MAE（表5）来评估不同方法的预测性能。由图17可看出，方法Ⅰ与方法Ⅱ的均方误差曲线几乎重合，再次说明陀螺仪退化数据包含的测量误差较小。采用方法Ⅲ与方法Ⅳ更新模型的参数后，均方误差逐渐减小，但方法Ⅳ的均方误差最小。由表5可知，方法Ⅰ的E_TMSE、E_RMSE、E_MAE最大，说明方法Ⅰ预测的误差最大；方法Ⅱ的E_TMSE、E_RMSE、E_MAE小于方法Ⅰ，说明退化数据虽包含测量误差，但测量误差较小；方法Ⅳ的E_TMSE、E_RMSE、E_MAE最小，说明方法Ⅳ预测的误差最小、预测效果最好。

将陀螺仪退化数据代入线性退化模型M₃，采用上述参数的估计与更新方法画出对应的剩余寿命概率密度函数曲线与剩余寿命曲线，如图18、图19所示。由图18可看出，4种方法都能较准确地预测陀螺仪的剩余寿命，但方法Ⅳ预测的剩余寿命不确定性更小。对比分析图12、图15、图18可知，图18中的剩余寿命概率密度函数曲线峰值更大、分布宽度更小，说明M₃预测的剩余寿命不确定性更小。这是因为陀螺仪的退化数据具有线性特征，更适合线性退化模型。

由图19可看出，基于模型M₃的4种方法预测的剩余寿命均与实际值接近，但方法Ⅳ的预测值最接近实际值，即方法Ⅳ具有最高的预测精度。4种方法的均方误差曲线见图20，4种方法的E_TMSE、E_RMSE、E_MAE如表6所示。由图20可知，方法Ⅰ与方法Ⅱ的均方误差曲线非常接近，说明退化数据包含的测量误差较小。方法Ⅳ的均方误差最小，表现最优。对比分析表4～表6可知，基于M₃预测的E_TMSE、E_RMSE、E_MA均小于基于M₁、M₂的预测结果。该结论与前文关于剩余寿命预测不确定性的分析结果相互印证。说明该退化数据更适合线性退化模型。由表6可知，方法Ⅰ的E_TMSE、E_RMSE、E_MA最大，误差最大；方法Ⅳ的E_TMSE、E_RMSE、E_MA最小，预测的误差最小、预测效果最好。

5.3 铝合金构件疲劳裂纹的增长数据

采用试验获取的铝合金构件疲劳裂纹增长数据预测构件的剩余寿命^［12］。所用数据为区间截尾数据，覆盖范围为1.5×10⁵～2.4×10⁵个循环周期（Cycle）。每组（共4组）数据包含10个连续的裂纹观测值，4组裂纹增长数据如图21所示。由图21可以看出，铝合金构件疲劳裂纹增长数据呈现明显的非线性特征，更适合非线性退化建模。

将该数据代入模型M₁，采用极大似然法、贝叶斯原理与EM算法分别估计与更新模型的漂移参数和扩散参数后，将参数代入剩余寿命概率密度函数，画出剩余寿命概率密度函数曲线（图22）、剩余寿命预测结果（图23）。由图22可看出，随着监测点数的增加，剩余寿命概率密度函数曲线的峰值逐渐增大、分布变窄，剩余寿命预测的不确定性减小。4种预测方法中，方法Ⅰ的剩余寿命概率密度曲线峰值最小、分布最宽，如图22e所示，这说明方法Ⅰ预测的不确定性最大。方法Ⅱ、方法Ⅲ的剩余寿命概率密度函数曲线基本一致，但高于方法Ⅰ的剩余寿命概率密度函数曲线。方法Ⅳ的剩余寿命概率密度函数曲线峰值最大、分布最窄，预测的不确定性最小。由图23可知，4种方法预测的剩余寿命均与实际寿命较接近，但方法Ⅳ的预测值最接近实际值，说明方法Ⅳ预测的精度最高。4种方法的均方误差曲线见图24，预测结果的E_TMSE、E_RMSE、E_MAE如表7所示。由图24可看出，方法Ⅰ与方法Ⅱ的均方误差曲线比较接近，说明铝合金构件疲劳裂纹增长数据包含的测量误差较小；相较于方法Ⅰ与方法Ⅱ，方法Ⅲ、方法Ⅳ的均方误差更小，说明更新模型参数的重要性；方法Ⅳ的均方误差最小，说明更新退化模型漂移系数的重要性。由表7可知，方法Ⅰ的E_TMSE、E_RMSE、E_MAE最大，方法Ⅳ的E_TMSE、E_RMSE、E_MAE最小；方法Ⅲ的E_TMSE、E_RMSE、E_MAE小于方法Ⅱ，说明动态更新退化模型中的参数有效提高了估计的准确性，减小了预测误差。本文提出的更新退化模型扩散参数方法更能减小剩余寿命的预测误差，提高剩余寿命预测精度。

将铝合金构件疲劳裂纹增长数据代入模型M₂，采用上述的参数估计与更新方法后，画出对应的剩余寿命概率密度函数曲线（图25）与剩余寿命曲线（图26）。由图25可看出，方法Ⅳ的剩余寿命概率密度函数曲线峰值最大、分布最窄，预测的不确定性最小；方法Ⅱ预测的不确定性小于方法Ⅰ，说明考虑测量误差减小了剩余寿命预测的不确定性；方法Ⅲ预测的不确定性小于方法Ⅱ，说明实时更新退化模型参数的重要性。由图26可知，方法Ⅳ预测的剩余寿命最接近于真实值。4种方法的均方误差曲线见图27，4种方法的E_TMSE、E_RMSE、E_MAE如表8所示。由图27可以看出，方法Ⅰ与方法Ⅱ的均方误差曲线几乎重合，说明铝合金构件疲劳裂纹增长数据包含的测量误差较小；方法Ⅳ的均方误差值更小，说明更新退化模型漂移系数可减小预测的均方误差。由表8可知，方法Ⅳ的E_TMSE、E_RMSE、E_MAE最小，说明提出的更新随机退化模型扩散参数方法减小了设备剩余寿命预测误差。

将铝合金构件疲劳裂纹增长数据代入模型M₃，画出线性退化模型对应的剩余寿命概率密度函数曲线（图28）与剩余寿命曲线（图29）。由图28可看出，4种方法预测的铝合金构件疲劳剩余寿命误差均较大，这是因为裂纹增长数据具有明显的非线性特征，更适用合非线性退化模型。方法Ⅰ预测的剩余寿命概率密度函数曲线峰值最小、分布最广，预测的不确定性最大。方法Ⅱ预测的剩余寿命不确定性小于方法Ⅰ，方法Ⅲ、方法Ⅳ预测的剩余寿命不确定性更小。对比图 16、图 18、图20的预测结果可知，模型M₃预测的剩余寿命概率密度函数曲线峰值最小、分布最广，说明模型M₃预测的不确定性大于模型M₁与M₂。因此，当退化数据具有非线性特征时，建立非线性退化模型才能准确估计设备的剩余寿命。由图29可知，方法Ⅳ预测的剩余寿命更接近于真实值。4种方法的均方误差曲线见图30，E_TMSE、E_RMSE、E_MAE如表9所示。

由图30可知，铝合金构件疲劳裂纹增长数据包含的测量误差较小，方法Ⅳ的均方误差最小。由表9可知，方法Ⅳ的E_TMSE、E_RMSE、E_MAE最小，预测误差最小，再次说明更新退化模型的漂移系数可减小预测误差，证明了本文所提方法的有效性。

6 结论

为充分利用设备的性能退化数据、准确估计退化模型的参数、减小设备的剩余寿命预测误差、提高剩余寿命的预测精度，本文提出一种同时在线更新随机退化模型的漂移系数和扩散系数的新方法。首先建立带测量误差的新型非线性退化模型、指数非线性退化模型与线性退化模型，根据首达时间原理、Wiener过程的性质推导出退化模型对应的寿命与剩余寿命概率密度函数的解析表达式。先采用极大似然法估计出退化模型的初始参数，再采用贝叶斯理论与EM算法同时更新退化模型的漂移参数和扩散参数。以高压脉冲电容退化数据、惯性导航平台陀螺漂移数据、铝合金构件疲劳裂纹增长数据为基础，分别对不含测量误差且不更新参数、含测量误差但不更新参数、含测量误差而仅更新漂移参数、含测量误差且同时更新漂移参数与扩散参数（所提方法）的退化模型进行仿真分析。仿真结果表明，所提方法能有效减小剩余寿命的预测误差，提高预测准确性。

参考文献

[1] 刘小峰， 亢莹莹， 柏林. 轴承自驱式独立退化轨迹构建与剩余寿命灰色预测［J］. 中国机械工程， 2024， 35（9）： 1613-1621.

[2] LIU Xiaofeng， KANG Yingying， BO Lin. Self-driven Independent Degradation Trajectory Construction and Remaining Life Gray Prediction for Bearings［J］. China Mechanical Engineering， 2024， 35（9）： 1613-1621.

[3] 王宇， 刘秋发， 彭一真. 非均匀监测条件下滚动轴承剩余寿命预测方法［J］. 机械工程学报， 2023， 59（23）： 96-104.

[4] WANG Yu， LIU Qiufa， PENG Yizhen. A Remaining Useful Life Prediction Approach with Nonuniform Monitoring Conditions for Rolling Bearings［J］. Journal of Mechanical Engineering， 2023， 59（23）： 96-104.

[5] WANG Chen， ZHANG Liming， CHEN Ling， et al. Remaining Useful Life Prediction of Nuclear Reactor Control Rod Drive Mechanism Based on Dynamic Temporal Convolutional Network［J］. Reliability Engineering & System Safety， 2025， 253： 110580.

[6] ZHOU Jianghong， QIN Yi. A Continuous Remaining Useful Life Prediction Method with Multistage Attention Convolutional Neural Network and Knowledge Weight Constraint［J］. IEEE Transactions on Neural Networks and Learning Systems， 2025， 36（7）： 11847-11860.

[7] MEEKER W Q， HAMADA M. Statistical Tools for the Rapid Development and Evaluation of High-reliability Products［J］. IEEE Transactions on Reliability， 1995， 44（2）： 187-198.

[8] LU C， MEEKER W， ESCOBAR L. A Comparison of Degradation and Failure-time Analysis Methods for Estimating a Time-to-failure Distribution［J］. Statistica Sinica， 1996， 6（3）： 531-546.

[9] KAHLE W， LEHMANN A. The Wiener Process as a Degradation Model： Modeling and Parameter Estimation［M］.Boston：Birkhäuser， 2010.

[10] 厉海涛， 金光， 周经伦， 等. 动量轮维纳过程退化建模与寿命预测［J］. 航空动力学报， 2011， 26（3）： 622-627.

[11] LI Haitao， JIN Guang， ZHOU Jinglun， et al. Momentum Wheel Wiener Process Degradation Modeling and Life Prediction［J］. Journal of Aerospace Power， 2011， 26（3）： 622-627.

[12] SI Xiaosheng， CHEN Maoyin， WANG Wenbin， et al. Specifying Measurement Errors for Required Lifetime Estimation Performance［J］. European Journal of Operational Research， 2013， 231（3）： 631-644.

[13] PENG C Y， TSENG S T. Statistical Lifetime Inference with Skew-Wiener Linear Degradation Models［J］. IEEE Transactions on Reliability， 2013， 62（2）： 338-350.

[14] GUAN Qingluan， WEI Xiukun， BAI Wenfei， et al. Two-stage Degradation Modeling for Remaining Useful Life Prediction Based on the Wiener Process with Measurement Errors［J］. Quality and Reliability Engineering International， 2022， 38（7）： 3485-3512.

[15] SI Xiaosheng， WANG Wenbin， HU Changhua， et al. Remaining Useful Life Estimation Based on a Nonlinear Diffusion Degradation Process［J］. IEEE Transactions on Reliability， 2012， 61（1）： 50-67.

[16] 彭才华， 钱富才， 杜许龙. 新型非线性退化建模与剩余寿命预测［J］. 计算机集成制造系统， 2019， 25（7）： 1647-1654.

[17] PENG Caihua， QIAN Fucai， DU Xulong. New Nonlinear Degradation Modeling and Residual Life Prediction［J］. Computer Integrated Manufacturing Systems， 2019， 25（7）： 1647-1654.

[18] YU Wennian， TU Wenbing， KIM I Y， et al. A Nonlinear-drift-driven Wiener Process Model for Remaining Useful Life Estimation Considering Three Sources of Variability［J］. Reliability Engineering & System Safety， 2021， 212： 107631.

[19] LIN Jingdong， LIAO Guobo， CHEN Min， et al. Two-phase Degradation Modeling and Remaining Useful Life Prediction Using Nonlinear Wiener Process［J］. Computers & Industrial Engineering， 2021， 160： 107533.

[20] 杨升. 基于贝叶斯更新的数控车床刀具剩余寿命预测研究［D］. 长春： 吉林大学， 2023.

[21] YANG Sheng. Research on cutting tool remaining useful life prediction of CNC lathe based on Bayesian updating［D］. Changchun： Jilin University， 2023.

[22] 任宏宇， 余瑶怡， 杜雄， 等. 基于优化长短期记忆神经网络的IGBT寿命预测模型［J］. 电工技术学报， 2024， 39（4）： 1074-1086.

[23] REN Hongyu， YU Yaoyi， DU Xiong， et al. IGBT Lifetime Prediction Model Based on Optimized Long Short-term Memory Neural Network［J］. Transactions of China Electrotechnical Society， 2024， 39（4）： 1074-1086.

[24] XU Xiaodong， TANG Shengjin， YU Chuanqiang， et al. Remaining Useful Life Prediction of Lithium-ion Batteries Based on Wiener Process under Time-varying Temperature Condition［J］. Reliability Engineering & System Safety， 2021， 214： 107675.

[25] FAN Mengfei， ZENG Zhiguo， ZIO E， et al. A Sequential Bayesian Approach for Remaining Useful Life Prediction of Dependent Competing Failure Processes［J］. IEEE Transactions on Reliability， 2019， 68（1）： 317-329.

[26] YE Zhisheng， WANG Yu， TSUI K L， et al. Degradation Data Analysis Using Wiener Processes with Measurement Errors［J］. IEEE Transactions on Reliability， 2013， 62（4）： 772-780.

[27] 司小胜， 胡昌华， 周东华. 带测量误差的非线性退化过程建模与剩余寿命估计［J］. 自动化学报， 2013， 39（5）： 530-541.

[28] SI Xiaosheng， HU Changhua， ZHOU Donghua. Nonlinear Degradation Process Modeling and Remaining Useful Life Estimation Subject to Measurement Error［J］. Acta Automatica Sinica， 2013， 39（5）： 530-541.

[29] PARK C， PADGETT W J. Stochastic Degradation Models with Several Accelerating Variables［J］. IEEE Transactions on Reliability， 2006， 55（2）： 379-390.

[30] 叶海福， 姚永和， 于成龙. 高压脉冲电容器电容量漂移的研究［J］. 高电压技术， 2007， 33（6）： 37-41.

[31] YE Haifu， YAO Yonghe， YU Chenglong. Capacitance Variability of High Voltage Pulse Capacitors［J］. High Voltage Engineering， 2007， 33（6）： 37-41.

[32] 王泽洲， 陈云翔， 蔡忠义， 等. 考虑随机失效阈值的设备剩余寿命在线预测［J］. 系统工程与电子技术， 2019， 41（5）： 1162-1168.

[33] WANG Zezhou， CHEN Yunxiang， CAI Zhongyi， et al. Real-time Prediction of Remaining Useful Lifetime for Equipment with Random Failure Threshold［J］. Systems Engineering and Electronics， 2019， 41（5）： 1162-1168.

基金资助

国家重点研发计划(2019YFB1707303)

国家重点研发计划(2020YFB17-13603)