航空发动机作为飞机的“心脏”，不仅为飞机提供动力，同时也为航空事业的发展注入强大动力。然而，航空发动机由于长期运行在复杂多变的恶劣环境下，其性能会随之产生衰减，并在累积到某种量级时，导致发动机失效^[1]。目前，航空发动机在航空领域发挥着重要的作用，一旦发生由失效引起的事故，所造成的后果将难以估计。因此，通过监测的航空发动机性能衰减数据，建立相应的性能衰减模型，进而实现航空发动机的剩余寿命(Remanining Useful Life, RUL)预测是近年来热门的一项关键技术^[2]。然而，航空发动机的寿命有多种定义，如首翻寿命、大修寿命、总寿命等，同时又受多种因素约束，如性能衰减、关键件寿命等。鉴于此，本文主要研究的是在首翻寿命的定义下，基于性能衰减建模的航空发动机剩余寿命预测。

目前，基于性能衰减建模的剩余寿命预测方法已经广泛地应用于航空发动机^[3-5]、导弹^[6-8]、电机^[9]等。例如，黄亮等^[3]基于多阶段Wiener过程建立了航空发动机的性能衰减模型，通过融合发动机的历史性能衰减数据和实时监测数据，实现对发动机的剩余寿命预测。王浩伟等^[6]针对导弹多个部件间的竞争退化，综合考虑了退化失效和突发失效2种竞争模式，提高了对导弹剩余寿命预测的精度。袁庆洋等^[9]针对电机的退化特性，建立了多段Wiener过程模型，实现了电机的剩余寿命预测。

然而，对于航空发动机这类复杂的机械系统，在性能衰减过程中普遍存在着三源不确定性，即时变不确定性^[10]、个体差异不确定性^[11]、测量不确定性^[12]。时变不确定性是指性能衰减过程随着时间变化的固有不确定性；个体差异不确定性是指同类型设备之间的性能衰减差异，即同类型设备的不同个体在性能衰减过程中虽然具有相似的衰减路径，但实际的衰减率是不相同的；测量不确定性是指测量的性能衰减数据与实际的性能衰减状态之间的差异。在实际中，由于受噪声、干扰、仪器精度等影响，想要对性能衰减状态的精确测量几乎是不可能的，因此测量中难免会引入测量误差，造成测量的不确定性。

目前，有关同时考虑三源不确定性的研究较少。为此，文献[13]基于线性的Wiener过程建立了同时考虑三源不确定性的性能衰减模型，并使得预测的剩余寿命分布可以随着性能衰减数据在线更新。然而，对于航空发动机而言，在性能衰减过程中不仅具有三源不确定性，还具有强非线性。鉴于此，文献[14]在非线性下考虑了三源不确定性，并且证明了同时考虑三源不确定性可以有效地提高剩余寿命预测的精度。文献[15]在此基础上，利用粒子滤波实现了对锂电池的剩余寿命预测。但是，以上有关考虑三源不确定性的工作中存在一个共性问题，即假定有多组同类型设备的历史性能衰减数据来离线地辨识模型初始参数，并且这些参数一旦确定，就不再随着新的性能衰减数据进行更新。因此，剩余寿命预测的精度受限于数据量，并且预测的结果在很大程度上依赖于选择的模型初始参数。然而，针对新研发的航空发动机，往往是缺乏此类相关的历史性能衰减数据和先验信息，因此，以上方法具有一定的局限性，且目前相关的研究较少。

对于目前航空发动机剩余寿命预测存在的问题，本文建立了同时考虑三源不确定性和非线性的性能衰减模型，并在首达时间下推导出剩余寿命的分布。然后，针对新研发的航空发动机，提出了一种基于Kalman滤波和条件期望最大化(Expectation Conditional Maximization, ECM)算法的参数估计方法，以克服缺乏历史性能衰减数据和先验信息的问题，并当新的监测数据可用时，能够自适应地估计模型参数，进而更新剩余寿命的分布，实现了对新研发航空发动机剩余寿命的在线预测。最后，通过实例验证了本文方法可以有效地提高预测的精度并降低预测的不确定性，且更加适用于新研发的航空发动机。

1 模型描述

针对航空发动机在性能衰减过程中存在的非线性和不确定性，本文采用一类一般性的非线性Wiener过程来建模航空发动机的性能衰减过程。具体地，性能衰减过程{X(t), t≥0}由标准Brown运动B(t)驱动，可以表示为

$ X\left( t \right) = {x_0} + \lambda \cdot \int_0^t {\rho \left( {\tau ;\theta } \right){\rm{d}}\tau } + {\sigma _{\rm{B}}}B\left( t \right) $	（1）

式中：λ为漂移系数；σ_B为扩散系数；ρ(t; θ)为时间t的非线性函数，用以描述航空发动机性能衰减的非线性特性；θ为未知参数，与σ_B共同表征同类型航空发动机之间的共性特征；B(t)为标准Brown运动，且有σ_BB(t)~N(0, σ_B²t)，用于刻画性能衰减过程固有的时变不确定性；x₀为初始性能衰减状态，不失一般性，假设x₀=0。进一步，ρ(t; θ)的不同函数形式可以描述不同的非线性性能衰减过程。显然，如果ρ(t; θ)=ρ，则式(1)转化为线性的性能衰减模型。因此，本文建立的性能衰减模型更具有一般性。

从工程实际出发，另外2种不确定性的主要来源是：①即使是同类型的航空发动机，受外部环境和内部因素的相互影响，不同个体之间的性能衰减过程往往也存在着差异。因此，有必要在性能衰减建模中考虑个体之间的差异不确定性。根据文献[13-15]，本文将漂移系数λ作为随机参数，以描述个体差异不确定性。进一步，假设λ~N(μ_λ, σ_λ²)，且与标准Brown运动B(t)独立。②由于受测量仪器精度、外部环境噪声等影响，理想的精确测量是很难实现的。因此，监测的性能衰减数据将不可避免地存在测量误差。在这种情况下，为了表征这种测量不确定性，令测量过程{Y(t), t>0}为

$ Y\left( t \right) = X\left( t \right) + \varepsilon $	（2）

式中：ε为随机测量误差，且ε~N(0, σ_ε²)。进一步假设ε、λ和B(t)相互独立。以上所有假设被广泛应用于性能衰减建模和剩余寿命预测领域^[11-16]。

在以上框架下，本文的主要目标是针对缺乏历史数据和先验信息的新研发航空发动机，利用实时监测的性能衰减数据来自适应地预测航空发动机的剩余寿命，并且在获得新的性能衰减数据后实现对剩余寿命的在线更新。具体地，首先假定获得的性能衰减数据是在离散时间点0=t₀ < t₁ < … < t_n监测的。为了方便起见，令y_n=Y(t_n)表示在t_n时刻监测的性能衰减数据，那么到t_n时刻所有性能衰减数据的集合可以表示为Y_0:n={y₀, y₁, …, y_n}。同样地，潜在性能衰减状态的集合可以表示为X_0:n={x₀, x₁, …, x_n}，其中，x_n=X(t_n)。因此，式(2)可以进一步写成y_n=x_n+ε_n，其中测量误差ε_n为ε在t_n时刻的具体实现。

为了实现三源不确定性下剩余寿命的自适应预测，建立随机参数λ的更新机制为λ_n=λ_n-1+γ，其中γ~N(0, σ_γ²)，初始随机参数λ₀~N(μ_λ, σ_λ²)。因此，考虑三源不确定性的性能衰减模型在状态空间模型的框架下可以构造为

$ \left\{ \begin{array}{l} {x_n} = {x_{n - 1}} + {\lambda _n}\mathit{\Lambda }\left( {{t_n};\theta } \right) - {\lambda _{n - 1}}\mathit{\Lambda }\left( {{t_{n - 1}};\theta } \right) + {\alpha _n}\\ {\lambda _n} = {\lambda _{n - 1}} + \gamma \\ {y_n} = {x_n} + {\varepsilon _n} \end{array} \right. $	（3）

式中：α_n=σ_B[B(t_n)-B(t_n-1)]，ε_n~N(0, σ_ε²)，Λ(t_n; θ)=$\int_{0}^{t_{n}} \rho(\tau; \theta) \mathrm{d} \tau$。根据标准Brown运动的性质，进一步有α_n~N[0, σ_B²(t_n-t_n-1)]。通过调查相关文献发现^[13-17]，大多数工作在进行性能衰减建模时存在一个潜在假设，即在t_n时刻估计的随机参数λ_n与t_n-1时刻随机参数的后验估计λ_n-1完全相等。然而，这导致了与随机参数的更新机制λ_n=λ_n-1+γ相矛盾。因此，为了解决这种潜在假设，本文建立了一个新的性能衰减模型，即式(3)。到此，三源不确定性已经被融入到非线性的性能衰减模型中。

2 剩余寿命自适应预测

考虑到表征设备健康水平的潜在性能衰减状态{X(t), t≥0}一旦首次达到设定的失效阈值ω(一般由工业标准设定)时，则认为设备失效。因此，基于首达时间的概念^[1-5]，设备在t_n时刻的剩余寿命L_n定义为

$ {L_n} = \inf \left\{ {{l_n} > 0:X\left( {{l_n} + {t_n}} \right) \ge \omega } \right\} $	（4）

相应地，剩余寿命L_n的概率密度函数(PDF)为f_Ln|Y_0:n(l_n|Y_0:n)。因此，预测航空发动机剩余寿命的关键是基于监测的性能衰减数据Y_0:n={y₀, y₁, …, y_n}，求解f_Ln|Y_0:n(l_n|Y_0:n)，并在获得新的性能衰减数据后，实现剩余寿命分布的在线更新。

为了实现剩余寿命的自适应预测，首先将式(3)进一步改写为

$ \left\{ {\begin{array}{*{20}{l}} {{x_n} = {x_{n - 1}} + {\lambda _{n - 1}}\int_{{t_{n - 1}}}^{{t_n}} \rho \left( {\tau ;\theta } \right){\rm{d}}\tau + {\beta _n}}\\ {{\lambda _n} = {\lambda _{n - 1}} + \gamma }\\ {{y_n} = {x_n} + {\varepsilon _n}} \end{array}} \right. $	（5）

式中：β_n=γ·Λ(t_n; θ+α_n)，根据正态分布的性质，可以进一步得到β_n~N(0, Q_n)，其中Q_n=σ_γ²[Λ(t_n; θ)²+σ_B²(t_n-t_n-1)。基于此，将潜在性能衰减状态x_n和随机参数λ_n视作“隐含”状态，进一步整理式(5)可以得到

$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{z}}_n} = {\mathit{\boldsymbol{A}}_n}{\mathit{\boldsymbol{z}}_{n - 1}} + {\mathit{\boldsymbol{\delta }}_n}}\\ {{y_n} = \mathit{\boldsymbol{C}}{\mathit{\boldsymbol{z}}_n} + {\varepsilon _n}} \end{array}} \right. $	（6）

式中：δ_n~N(0, φ_n)，且有

$ \begin{array}{l} {\mathit{\boldsymbol{z}}_n} = \left[ {\begin{array}{*{20}{c}} {{x_n}}\\ {{\lambda _n}} \end{array}} \right],{\mathit{\boldsymbol{\delta }}_n} = \left[ {\begin{array}{*{20}{c}} {{\beta _n}}\\ \gamma \end{array}} \right],{\mathit{\boldsymbol{A}}_n} = \left[ {\begin{array}{*{20}{c}} 1&{\int_{{t_{n - 1}}}^{{t_n}} \rho \left( {\tau ;\theta } \right){\rm{d}}\tau }\\ 0&1 \end{array}} \right],\\ \mathit{\boldsymbol{C}} = {\left[ {\begin{array}{*{20}{l}} 1\\ 0 \end{array}} \right]^{\rm{T}}},{\mathit{\boldsymbol{\varphi }}_n} = \left[ {\begin{array}{*{20}{l}} {{Q_n}}&{{S_n}}\\ {{S_n}}&{\sigma _\gamma ^2} \end{array}} \right] \end{array} $

其中:S_n为β_n和γ的协方差，根据公式可得S_n=σ_γ²Λ(t_n; θ)。

定义基于监测数据Y_0:n估计的z_n的期望和协方差分别为

$ \left\{ \begin{gathered} {{\mathit{\boldsymbol{\hat z}}}_{n\left| n \right.}} = \left[ {\begin{array}{*{20}{l}} {{{\hat x}_{n\left| n \right.}}} \\ {{{\hat \lambda }_{n\left| n \right.}}} \end{array}} \right] = \mathbb{E}\left( {{\mathit{\boldsymbol{z}}_n}\left| {{Y_{0;n}}} \right.} \right) \hfill \\ {\mathit{\boldsymbol{P}}_{n\left| n \right.}} = \left[ {\begin{array}{*{20}{c}} {\kappa _{x,n}^2}&{\kappa _{x\lambda ,n}^2} \\ {\kappa _{x\lambda ,n}^2,}&{\kappa _{\lambda ,n}^2} \end{array}} \right] = cov\left( {{\mathit{\boldsymbol{z}}_n}\left| {{Y_{0;n}}} \right.} \right) \hfill \\ \end{gathered} \right. $

（7）

式中：$\hat{x}_{n | n}=\mathbb{E}\left(x_{n} | Y_{0; n}\right)$; $\hat{\lambda}_{n | n}=\mathbb{E}\left(\lambda_{n} | Y_{0; n}\right)$; $\kappa_{x, n}^{2}=\operatorname{var}\left(x_{n} | Y_{0; n}\right)$; $\kappa_{\lambda, n}^{2}=\operatorname{var}\left(\lambda_{n} | Y_{0; n}\right)$; $\kappa_{x \lambda, n}^{2}=\operatorname{cov}\left(x_{n} \lambda_{n} | Y_{0, n}\right)$。

同样地，定义一步前向的期望和协方差为

$ \left\{ \begin{gathered} {{\mathit{\boldsymbol{\hat z}}}_{n\left| {n - 1} \right.}} = \left[ {\begin{array}{*{20}{l}} {{{\hat x}_{n\left| {n - 1} \right.}}} \\ {{{\hat \lambda }_{n\left| {n - 1} \right.}}} \end{array}} \right] = \mathbb{E}\left( {{\mathit{\boldsymbol{z}}_n}\left| {{Y_{0;n - 1}}} \right.} \right) \hfill \\ {\mathit{\boldsymbol{P}}_{n\left| {n - 1} \right.}} = \left[ {\begin{array}{*{20}{c}} {\kappa _{x,n\left| {n - 1} \right.}^2}&{\kappa _{x\lambda ,n\left| {n - 1} \right.}^2} \\ {\kappa _{x\lambda ,n\left| {n - 1} \right.}^2}&{\kappa _{\lambda ,n\left| {n - 1} \right.}^2} \end{array}} \right] = {\rm{cov}}\left( {{\mathit{\boldsymbol{z}}_n}\left| {{Y_{0;n - 1}}} \right.} \right) \hfill \\ \end{gathered} \right. $

（8）

根据以上设定，假定初始的潜在性能衰减状态x₀~N(0, σ_x²)，则由潜在性能衰减状态和随机参数组成的z_n可以通过Kalman滤波进行联合估计，实现过程如下：

状态估计

$ \left\{ \begin{array}{l} {{\mathit{\boldsymbol{\hat z}}}_{n\left| {n - 1} \right.}} = {\mathit{\boldsymbol{A}}_n}{{\mathit{\boldsymbol{\hat z}}}_{n - 1\left| {n - 1} \right.}}\\ {{\mathit{\boldsymbol{\hat z}}}_{n\left| n \right.}} = {{\mathit{\boldsymbol{\hat z}}}_{n\left| {n - 1} \right.}} + {\mathit{\boldsymbol{K}}_n}\left( {{y_n} - \mathit{\boldsymbol{C}}{{\mathit{\boldsymbol{\hat z}}}_{n\left| {n - 1} \right.}}} \right)\\ {\mathit{\boldsymbol{K}}_n} = {\mathit{\boldsymbol{P}}_{n\left| {n - 1} \right.}}{\mathit{\boldsymbol{C}}^{\rm{T}}}{\left( {\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{P}}_{n\left| {n - 1} \right.}}{\mathit{\boldsymbol{C}}^{\rm{T}}} + \sigma _\varepsilon ^2} \right)^{ - 1}}\\ {\mathit{\boldsymbol{P}}_{n\left| {n - 1} \right.}} = {\mathit{\boldsymbol{A}}_n}{\mathit{\boldsymbol{P}}_{n - 1\left| {n - 1} \right.}}\mathit{\boldsymbol{A}}_n^{\rm{T}} + {\mathit{\boldsymbol{\varphi }}_n} \end{array} \right. $

（9）

协方差更新

$ {\mathit{\boldsymbol{P}}_{n\left| n \right.}} = {\mathit{\boldsymbol{P}}_{n\left| {n - 1} \right.}} - {\mathit{\boldsymbol{K}}_n}\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{P}}_{n\left| {n - 1} \right.}} $	（10）

滤波初值

$ {{\mathit{\boldsymbol{\hat z}}}_{0\left| 0 \right.}} = \left[ {\begin{array}{*{20}{c}} 0\\ {{\mu _\lambda }} \end{array}} \right],{\mathit{\boldsymbol{P}}_{0\left| 0 \right.}} = \left[ {\begin{array}{*{20}{c}} {\sigma _x^2}&0\\ 0&{\sigma _\lambda ^2} \end{array}} \right] $	（11）

根据Kalman滤波的高斯性质，z_n的后验PDF是双变量高斯分布的，即z_n~$N\left(\hat{\boldsymbol{z}}_{n | n}, \boldsymbol{P}_{n | n}\right)$。因此，z₀~N(μ₀, P₀)，其中μ₀=$\hat{\boldsymbol{z}}_{0|0}$，P₀=P_0|0。

根据建立的状态空间模型式(6)，基于监测的性能衰减数据Y_0:n，计算剩余寿命的PDF为

$ \begin{array}{l} {f_{{L_n}\left| {{Y_{0,n}}} \right.}}\left( {{l_n}|{Y_{0,n}}} \right) \cong \\ \frac{{{{\cal B}_n}\left( {{\cal D}_n^2\kappa _{\lambda ,n}^2 + {{\cal F}_n}} \right) - {{\cal A}_n}{{\cal C}_n}{\cal D}_n^2\kappa _{\lambda ,n}^2 - {{\cal C}_n}{{\cal F}_n}{{\hat \lambda }_{n\left| n \right.}}}}{{{{\cal F}_n}\sqrt {2{\rm{ \mathsf{ π} }}{{\left( {{\cal D}_n^2\kappa _{\lambda ,n}^2 + {{\cal F}_n}} \right)}^3}} }} \times \\ \exp \left[ { - \frac{{{{\left( {\omega - {{\hat x}_{n\left| n \right.}} - {{\hat \lambda }_{n\left| n \right.}}\int_{{t_n}}^{{l_n} + {t_n}} {\rho \left( {\tau ,\theta } \right){\rm{d}}\tau } } \right)}^2}}}{{2\left( {{\cal D}_n^2\kappa _{\lambda ,n}^2 + {{\cal F}_n}} \right)}}} \right] \end{array} $

（12）

式中：$\mathcal{A}_{n}=\omega-\hat{x}_{n | n}+\varphi \hat{\lambda}_{n | n}$，$\mathcal{B}_{n}=\left(\omega-\hat{x}_{n | n}+\varphi \hat{\lambda}_{n | n}\right) \sigma_{\mathrm{B}}^{2}$，$\mathcal{C}_{n}\left[\varphi+\int_{t_{n}}^{l_{n}+t_{n}} \rho(\tau; \theta) \mathrm{d} \tau+l_{n}\left(\int_{t_{n}}^{l_{n}+t_{n}} \rho(\tau; \theta) \mathrm{d} \tau\right)^{\prime}\right] \sigma_{\mathrm{B}}^{2}$-$\left(\int_{t_{n}}^{l_{n}+t_{n}} \rho(\tau; \theta) \mathrm{d} \tau\right)^{\prime}\left(\kappa_{x, n}^{2}-\varphi \kappa_{x \lambda, n}^{2}\right)$，$\varphi=\frac{\kappa_{x \lambda, n}^{2}}{\kappa_{\lambda, n}^{2}}$，$\mathcal{D}_{n}=$$\int_{t_{n}}^{l_{n}+t_{n}} \rho(\tau; \theta) \mathrm{d} \tau+\varphi, \mathcal{F}_{n}=\kappa_{x, n}^{2}-\varphi \kappa_{x \lambda, n}^{2}+\sigma_{\mathrm{B}}^{2} l_{n}$。

这里，式(12)的证明过程可以参照文献[14]，故在此省略。基于以上结果，在获得新的性能衰减数据后，可以通过式(12)对剩余寿命进行实时预测。然而，其中的模型参数μ_λ、σ_λ²、σ_x²、σ_B²、σ_γ²、σ_ε²、θ均是未知的。通常，许多相关文献利用多组同类型设备的历史性能衰减数据来确定模型参数，并且这些参数一旦确定，就不再随着新的性能衰减数据进行在线更新^{[2-6, 10-12, 14-16]}。然而，在实际中经常会出现缺乏此类历史性能衰减数据的情况，尤其是本文所研究的新研发航空发动机。因此，需要利用实时监测的性能衰减数据对模型参数进行自适应估计，使得预测的剩余寿命更加准确，有利于掌握设备当前的性能衰减程度。

3 基于ECM算法的模型参数估计

为了实现模型参数自适应估计，首先令Θ=[μ₀^T, vec(P₀)^T, σ_B², σ_γ², σ_ε², θ]^T表示模型参数向量，其中vec(·)是矩阵拉直运算符，并且

$ {\mathit{\boldsymbol{\mu }}_0} = \left[ {\begin{array}{*{20}{c}} 0\\ {{\mu _\lambda }} \end{array}} \right],{\mathit{\boldsymbol{P}}_0} = \left[ {\begin{array}{*{20}{c}} {\sigma _x^2}&0\\ 0&{\sigma _\lambda ^2} \end{array}} \right] $	（13）

根据极大似然估计(MLE)理论，性能衰减数据Y_0:n关于模型参数向量Θ的对数似然函数为

$ {{\cal L}_n}\left( \mathit{\boldsymbol{ \boldsymbol{\varTheta} }} \right) = \ln p\left( {{Y_{0:n}}\left| \mathit{\boldsymbol{ \boldsymbol{\varTheta} }} \right.} \right) $	（14）

式中：p(Y_0:n|Θ)为联合PDF。进一步，模型参数向量Θ的MLE值${\mathit{\boldsymbol{ \boldsymbol{\hat \varTheta} }}_n}$可以通过最大化式(14)得到，具体为

$ {{\mathit{\boldsymbol{ \boldsymbol{\hat \varTheta} }}}_n} = \arg \mathop {\max }\limits_\mathit{\boldsymbol{ \boldsymbol{\varTheta} }} {{\cal L}_n}\left( \mathit{\boldsymbol{ \boldsymbol{\varTheta} }} \right) $	（15）

然而，由于将潜在性能衰减状态x_n和随机参数λ_n视作“隐含”状态，因此无法直接最大化式(15)。为解决这一问题，自然会采用期望最大化(EM)算法进行参数估计，通过最大化联合似然函数p(z_n, Y_0:n|Θ)去逼近模型参数的MLE。EM算法的估计过程可以分为以下2步。

E步：

$ \ell \left( {\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}\left| {\mathit{\boldsymbol{ \boldsymbol{\hat \varTheta} }}_n^{\left( m \right)}} \right.} \right) = {\mathbb{E}_{{z_n}\left| {{Y_{0:n}}} \right.,\mathit{\boldsymbol{ \boldsymbol{\hat \varTheta} }}_n^{\left( m \right)}}}\left\{ {\ln p\left( {{\mathit{\boldsymbol{z}}_n},{Y_{0:n}}\left| \mathit{\boldsymbol{ \boldsymbol{\varTheta} }} \right.} \right)} \right\} $

（16）

式中：${\mathit{\boldsymbol{ \boldsymbol{\hat \varTheta} }}_n}^{\left(m \right)}$为基于性能衰减数据Y_0:n在第m次迭代时得到的结果。

M步：

$ \mathit{\boldsymbol{ \boldsymbol{\hat \varTheta} }}_n^{\left( {m + 1} \right)} = \arg \mathop {\max }\limits_\mathit{\boldsymbol{ \boldsymbol{\varTheta} }} \left\{ {\ell \left( {\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}\left| {\mathit{\boldsymbol{ \boldsymbol{\hat \varTheta} }}_n^{\left( m \right)}} \right.} \right)} \right\} $	（17）

根据EM算法理论^[18]，假定从第m次迭代得到的${\mathit{\boldsymbol{ \boldsymbol{\hat \varTheta} }}_n}^{\left(m \right)}$出发，通过计算E步和M步，可以在第m+1次迭代时得到一个更好的MLE值${\mathit{\boldsymbol{ \boldsymbol{\hat \varTheta} }}_n}^{\left(m+1 \right)}$。在实际应用中，仅执行一次迭代很难得到令人满意的估计值。因此，EM算法一般需要执行多次迭代。

首先，为了计算对数联合似然函数，可以通过条件概率公式得到

$ \begin{array}{l} 2\ln p\left( {{\mathit{\boldsymbol{z}}_n},{Y_{0:n}}\left| \mathit{\boldsymbol{ \boldsymbol{\varTheta} }} \right.} \right) = \\ \;\;\;2\left[ {\ln p\left( {{Y_{0:n}}|{\mathit{\boldsymbol{z}}_n},\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}} \right) + \ln p\left( {{\mathit{\boldsymbol{z}}_n}|\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}} \right)} \right] = \\ \;\;\; - \ln \left| {{\mathit{\boldsymbol{P}}_0}} \right| - {\left( {{\mathit{\boldsymbol{z}}_0} - {\mathit{\boldsymbol{\mu }}_0}} \right)^{\rm{T}}}\mathit{\boldsymbol{P}}_0^{ - 1}\left( {{\mathit{\boldsymbol{z}}_0} - {\mathit{\boldsymbol{\mu }}_0}} \right) - \\ \;\;\;\sum\limits_{i = 1}^n {\left[ {\ln \left| {{\mathit{\boldsymbol{\varphi }}_i}} \right| + {{\left( {{\mathit{\boldsymbol{z}}_i} - {\mathit{\boldsymbol{A}}_i}{\mathit{\boldsymbol{z}}_{i - 1}}} \right)}^{\rm{T}}}\mathit{\boldsymbol{\varphi }}_i^{ - 1}\left( {{\mathit{\boldsymbol{z}}_i} - {\mathit{\boldsymbol{A}}_i}{\mathit{\boldsymbol{z}}_{i - 1}}} \right)} \right]} - \\ \;\;\;n\ln \sigma _\varepsilon ^2 - \sum\limits_{i = 1}^n {\frac{{{{\left( {{y_i} - \mathit{\boldsymbol{C}}{\mathit{\boldsymbol{z}}_i}} \right)}^2}}}{{\sigma _\varepsilon ^2}}} \end{array} $

（18）

然后，通过E步对式(18)求条件期望，有

$ \begin{gathered} \ell \left( {\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}\left| {\mathit{\boldsymbol{ \boldsymbol{\hat \varTheta} }}_n^{\left( m \right)}} \right.} \right) = \hfill \\ \;\;{\mathbb{E}_{{z_n}\left| {{Y_{0:n}}} \right.,\mathit{\boldsymbol{ \boldsymbol{\hat \varTheta} }}_n^{\left( m \right)}}}\left\{ {2\ln p\left( {{\mathit{\boldsymbol{z}}_n},{Y_{0:n}}\left| \mathit{\boldsymbol{ \boldsymbol{\varTheta} }} \right.} \right)} \right\} = \hfill \\ \;\;\; - Tr\left\{ {\mathit{\boldsymbol{P}}_0^{ - 1}{\mathbb{E}_{{z_n}\left| {{Y_{0:n}}} \right.,\mathit{\boldsymbol{ \boldsymbol{\hat \varTheta} }}_n^{\left( m \right)}}}\left[ {\left( {{\mathit{\boldsymbol{z}}_0} - {\mathit{\boldsymbol{\mu }}_0}} \right){{\left( {{\mathit{\boldsymbol{z}}_0} - {\mathit{\boldsymbol{\mu }}_0}} \right)}^{\text{T}}}} \right]} \right\} \hfill \\ \;\;\; - \sum\limits_{i = 1}^n {\left\{ {\ln \left| {{\mathit{\boldsymbol{\varphi }}_i}} \right| + Tr \left\{ {\mathit{\boldsymbol{\varphi }}_i^{ - 1}\left[ {\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} - \mathit{\boldsymbol{ \boldsymbol{\varXi} A}}_i^{\text{T}} - {\mathit{\boldsymbol{A}}_i}{\mathit{\boldsymbol{ \boldsymbol{\varXi} }}^{\text{T}}} + } \right.} \right.} \right.} \hfill \\ \;\;\;\left. {\left. {\left. {{\mathit{\boldsymbol{A}}_i}\mathit{\boldsymbol{ \boldsymbol{\varPhi} A}}_i^{\text{T}}} \right]} \right\}} \right\} - Tr\left\{ {\sigma _\varepsilon ^{ - 2}\sum\limits_{i = 1}^n {\left[ {{{\left( {{y_i} - \mathit{\boldsymbol{C}}{\mathit{\boldsymbol{z}}_{i\left| n \right.}}} \right)}^2} + } \right.} } \right. \hfill \\ \;\;\;\left. {\left. {\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{P}}_{i\left| n \right.}}{\mathit{\boldsymbol{C}}^{\text{T}}}} \right]} \right\} - \ln \left| {{\mathit{\boldsymbol{P}}_0}} \right| - n\ln \sigma _\varepsilon ^2 \hfill \\ \end{gathered} $

（19）

式中：

$ \left\{ \begin{gathered} \mathit{\boldsymbol{ \boldsymbol{\varOmega} }} = {\mathbb{E}_{{z_n}\left| {{Y_{0:n}}} \right.,\mathit{\boldsymbol{ \boldsymbol{\hat \varTheta} }}_n^{\left( m \right)}}}\left( {{\mathit{\boldsymbol{z}}_i}\mathit{\boldsymbol{z}}_i^{\text{T}}} \right) \hfill \\ \mathit{\boldsymbol{ \boldsymbol{\varXi} }} = {\mathbb{E}_{{z_n}\left| {{Y_{0:n}}} \right.,\mathit{\boldsymbol{ \boldsymbol{\hat \varTheta} }}_n^{\left( m \right)}}}\left( {{\mathit{\boldsymbol{z}}_i}\mathit{\boldsymbol{z}}_{i - 1}^{\text{T}}} \right) \hfill \\ \mathit{\boldsymbol{ \boldsymbol{\varPhi} }} = {\mathbb{E}_{{z_n}\left| {{Y_{0:n}}} \right.,\mathit{\boldsymbol{ \boldsymbol{\hat \varTheta} }}_n^{\left( m \right)}}}\left( {{\mathit{\boldsymbol{z}}_{i - 1}}\mathit{\boldsymbol{z}}_{i - 1}^{\text{T}}} \right) \hfill \\ \end{gathered} \right. $

（20）

显然，为了推导出E步结果，需要先计算式(20)。这里，本文采用Rauch-Tung-Striebel(RTS)平滑滤波算法^[19]计算Ω、Ξ、Φ：

平滑滤波

$ \left\{ {\begin{array}{*{20}{l}} {{\Im _i} = {\mathit{\boldsymbol{P}}_{i|i}}\mathit{\boldsymbol{A}}_{i + 1}^{\rm{T}}\mathit{\boldsymbol{P}}_{i + 1|i}^{ - 1}}\\ {{{\mathit{\boldsymbol{\hat z}}}_{i|n}} = {{\mathit{\boldsymbol{\hat z}}}_{i|i}} + {\Im _i}\left( {{{\mathit{\boldsymbol{\hat z}}}_{i + 1|n}} - {{\mathit{\boldsymbol{\hat z}}}_{i + 1|i}}} \right)}\\ {{\mathit{\boldsymbol{P}}_{i|n}} = {\mathit{\boldsymbol{P}}_{i|i}} + {\Im _i}\left( {{\mathit{\boldsymbol{P}}_{i + 1|n}} - {\mathit{\boldsymbol{P}}_{i + 1|i}}} \right)\Im _i^{\rm{T}}} \end{array}} \right. $

（21）

协方差初值

$ {\mathit{\boldsymbol{M}}_{n|n}} = \left( {{\mathit{\boldsymbol{I}}_{2 \times 2}} - {\mathit{\boldsymbol{K}}_n}\mathit{\boldsymbol{C}}} \right){\mathit{\boldsymbol{A}}_n}{\mathit{\boldsymbol{P}}_{n - 1|n - 1}} $	（22）

后向迭代

$ {\mathit{\boldsymbol{M}}_{i|n}} = {\mathit{\boldsymbol{P}}_{i|i}}\Im _{i - 1}^{\rm{T}} + {\Im _i}\left( {{\mathit{\boldsymbol{M}}_{i + 1|n}} - {\mathit{\boldsymbol{A}}_{i + 1}}{\mathit{\boldsymbol{P}}_{i|i}}} \right)\Im _{i - 1}^{\rm{T}} $	（23）

式中：M_i|n=cov(z_i, z_i-1|Y_0:n)。基于RTS算法，条件期望Ω、Ξ、Φ可以计算为

$ \left\{ \begin{array}{l} \mathit{\boldsymbol{ \boldsymbol{\varOmega} }} = {{\mathit{\boldsymbol{\hat z}}}_{i|n}}\mathit{\boldsymbol{\hat z}}_{i|n}^{\rm{T}} + {\mathit{\boldsymbol{P}}_{i|n}}\\ \mathit{\boldsymbol{ \boldsymbol{\varXi} }} = {{\mathit{\boldsymbol{\hat z}}}_{i|n}}\mathit{\boldsymbol{\hat z}}_{i - 1|n}^{\rm{T}} + {\mathit{\boldsymbol{M}}_{i|n}}\\ \mathit{\boldsymbol{ \boldsymbol{\varPhi} }} = {{\mathit{\boldsymbol{\hat z}}}_{i - 1|n}}\mathit{\boldsymbol{\hat z}}_{i - 1|n}^{\rm{T}} + {\mathit{\boldsymbol{P}}_{i - 1|n}} \end{array} \right. $

（24）

基于以上结果，可以得到$\ell \left({\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}|\mathit{\boldsymbol{ \boldsymbol{\hat \varTheta} }}_n^{(m)}} \right)$的解析表达式，因此，接下来利用M步计算模型参数向量Θ的MLE值。为了最大化$\ell \left({\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}|\mathit{\boldsymbol{ \boldsymbol{\hat \varTheta} }}_n^{(m)}} \right)$，直接的方法是令$\partial \ell \left({\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}|\mathit{\boldsymbol{ \boldsymbol{\hat \varTheta} }}_n^{(m)}} \right)$/$\partial \mathit{\boldsymbol{ \boldsymbol{\varTheta} }}$=0，那么在第m+1次迭代时，可以得到以下结果:

$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{\mu }}_{0,n}^{\left( {m + 1} \right)} = {{\mathit{\boldsymbol{\hat z}}}_{0|n}}}\\ {\mathit{\boldsymbol{P}}_{0,n}^{\left( {m + 1} \right)} = {\mathit{\boldsymbol{P}}_{0|n}}}\\ {\sigma _{\varepsilon ,n}^{{2^{(m + 1)}}} = \frac{1}{n}\sum\limits_{i = 1}^n {\left[ {{{\left( {{y_i} - \mathit{\boldsymbol{C}}{\mathit{\boldsymbol{z}}_{i|n}}} \right)}^2} + \mathit{\boldsymbol{C}}{\mathit{\boldsymbol{P}}_{i|n}}{\mathit{\boldsymbol{C}}^{\rm{T}}}} \right]} } \end{array}} \right. $

（25）

式中：μ_{0, n}^(m+1)、P_{0, n}^(m+1)、σ_{ε, n}^2(m+1)为$\arg \mathop {\max }\limits_\mathit{\boldsymbol{ \boldsymbol{\varTheta} }} \left\{ {\ell \left({\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}|\mathit{\boldsymbol{ \boldsymbol{\hat \varTheta} }}_n^{(m)}} \right)} \right\}$的全局最优解。

但需要注意的是，这里仅有μ₀、P₀和σ_ε²可以得到解析表示，这是由于性能衰减模型式(6)中的矩阵A_n和φ_n受参数θ的影响，其他参数σ_B²、σ_γ²和θ无法通过求解$\partial \ell \left({\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}|\mathit{\boldsymbol{ \boldsymbol{\hat \varTheta} }}_n^{(m)}} \right)/\partial \mathit{\boldsymbol{ \boldsymbol{\varTheta} }} = 0$得到。如果直接采用启发式的优化算法进行求解，不仅无法保证算法的收敛性，而且可能会导致较差的在线能力。鉴于此，本文采用ECM算法来克服传统EM算法的缺陷^[20-21]。根据ECM算法将参数估计过程分为以下2个CM步。

CM1：首先固定参数θ，利用上一步迭代值θ_n^(m)进行替代，即可利用式(19)对模型参数向量Θ求导。这里，注意式(19)中的矩阵A_i和φ_i是与监测时刻相关的，因此先将式(19)中的第2项改写为

$ \begin{array}{l} - \sum\limits_{i = 1}^n {\left\{ {\ln \left| {{\mathit{\boldsymbol{\varphi }}_i}} \right| + Tr\left\{ {\mathit{\boldsymbol{\varphi }}_i^{ - 1}\left[ {\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} - \mathit{\boldsymbol{ \boldsymbol{\varXi} A}}_i^{\rm{T}} - {\mathit{\boldsymbol{A}}_i}{\mathit{\boldsymbol{ \boldsymbol{\varXi} }}^{\rm{T}}} + } \right.} \right.} \right.} \\ \;\;\;\;\left. {\left. {\left. {{\mathit{\boldsymbol{A}}_i}\mathit{\boldsymbol{ \boldsymbol{\varPhi} A}}_i^{\rm{T}}} \right]} \right\}} \right\} = \\ - \sum\limits_{i = 1}^n {\left\{ {\ln \left| {{\mathit{\boldsymbol{\varphi }}_i}} \right| + Tr\left\{ {\mathit{\boldsymbol{\varphi }}_i^{ - 1}{\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}_i}\left[ {\theta _n^{(m)}} \right]} \right\}} \right\}} = \\ \;\;\;\; - \sum\limits_{i = 1}^n {\left\{ {\ln \left[ {\sigma _{\rm{B}}^2\sigma _\gamma ^2\left( {{t_i} - {t_{i - 1}}} \right)} \right]} \right. + } \\ \;\;\;\;\frac{{{d_i}\left[ {\theta _n^{(m)}} \right]{{\left[ {\mathit{\Lambda }\left( {{t_i};\theta _n^{(m)}} \right)} \right]}^2}}}{{\sigma _{\rm{B}}^2\left( {{t_i} - {t_{i - 1}}} \right)}} + \\ \;\;\;\;\frac{{{a_i}\left[ {\theta _n^{(m)}} \right] - \left( {{b_i}\left[ {\theta _n^{(m)}} \right] + } \right.}}{{\sigma _{\rm{B}}^2\left( {{t_i} - {t_{i - 1}}} \right)}}\\ \;\;\;\;\frac{{\left. {{c_i}\left[ {\theta _n^{(m)}} \right]} \right)\mathit{\Lambda }\left( {{t_i};\theta _n^{(m)}} \right)}}{{\sigma _{\rm{B}}^2\left( {{t_i} - {t_{i - 1}}} \right)}} + \\ \;\;\;\;\left. {\frac{{{d_i}\left[ {\theta _n^{(m)}} \right]}}{{\sigma _\gamma ^2}}} \right\} \end{array} $

（26）

式中：Λ(t_i; θ_n^(m))=$\int_{0}^{t_{i}} \rho\left(\tau; \theta_{n}^{(m)}\right) \mathrm{d} \tau$，且Ψ_i[θ_n^(m)]为参数θ_n^(m)的函数矩阵，表示为

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}_i}\left[ {\theta _n^{(m)}} \right] = \mathit{\boldsymbol{ \boldsymbol{\varOmega} }} - \mathit{\boldsymbol{ \boldsymbol{\varXi} A}}_i^{\rm{T}} - {\mathit{\boldsymbol{A}}_i}{\mathit{\boldsymbol{ \boldsymbol{\varXi} }}^{\rm{T}}} + {\mathit{\boldsymbol{A}}_i}\mathit{\boldsymbol{ \boldsymbol{\varPhi} A}}_i^{\rm{T}} = }\\ {\left[ {\begin{array}{*{20}{c}} {{a_i}\left[ {\theta _n^{\left( m \right)}} \right]}&{{b_i}\left[ {\theta _n^{\left( m \right)}} \right]}\\ {{c_i}\left[ {\theta _n^{\left( m \right)}} \right]}&{{d_i}\left[ {\theta _n^{\left( m \right)}} \right]} \end{array}} \right]} \end{array} $

（27）

其中：a_i[θ_n^(m)]、b_i[θ_n^(m)]、c_i[θ_n^(m)]、d_i[θ_n^(m)]为矩阵Ψ_i[θ_n^(m)]的元素值。

基于式(19)和式(26)，在参数θ固定的条件下，σ_B²和σ_γ²可以通过令$\partial \ell \left({\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}|\mathit{\boldsymbol{ \boldsymbol{\hat \varTheta} }}_n^{(m)}} \right)/\partial \mathit{\boldsymbol{ \boldsymbol{\varTheta} }} = 0$得到解析表达式，即

$ \begin{array}{l} \sigma _{{\rm{B}},n}^{{2^{\left( {m + 1} \right)}}} = \frac{1}{n}\sum\limits_{i = 1}^n {\left\{ {\frac{{{a_i}\left[ {\theta _n^{\left( m \right)}} \right] - {b_i}\left[ {\theta _n^{\left( m \right)}} \right]\mathit{\Lambda }\left( {{t_i};\theta _n^{(m)}} \right)}}{{{t_i} - {t_{i - 1}}}} - \frac{{{c_i}\left[ {\theta _n^{(m)}} \right]\mathit{\Lambda }\left( {{t_i};\theta _n^{(m)}} \right) - {d_i}\left[ {\theta _n^{\left( m \right)}} \right]{{\left[ {\mathit{\Lambda }\left( {{t_i};\theta _n^{(m)}} \right)} \right]}^2}}}{{{t_i} - {t_{i - 1}}}}} \right\}} \\ \sigma _{\gamma ,n}^{{2^{\left( {m + 1} \right)}}} = \frac{1}{n}\sum\limits_{i = 1}^n {\left\{ {{d_i}\left[ {\theta _n^{\left( m \right)}} \right]} \right\}} \end{array} $

（28）

CM2：为了基于性能衰减数据Y_0:n估计参数θ，首先将式(25)和式(28)一起代入到式(19)中，便可以通过最大化$\ell$(θ|μ_{0, n}^(m+1), P_{0, n}^(m+1), σ_{B, n}^2(m+1), σ_{γ, n}^2(m+1), σ_{ε, n}^2(m+1))来估计参数θ_n^(m+1)，即

$ \begin{array}{l} \theta _n^{\left( {m + 1} \right)} = \arg \mathop {\max }\limits_\theta \left\{ {\ln \left| {{\mathit{\boldsymbol{P}}_{0|n}}} \right| - 2 - n - n\ln \sigma _{\varepsilon ,n}^{{2^{(m + 1)}}} - \sum\limits_{i = 1}^n {\left\{ {\ln \left[ {\sigma _{{\rm{B}},n}^{{2^{(m + 1)}}}\sigma _{\gamma ,n}^{{2^{(m + 1)}}}\left( {{t_i} - {t_{i - 1}}} \right)} \right.} \right.} } \right\} - \\ \;\;\;\;\;\;\;\left. {\sum\limits_{i = 1}^n {\left\{ {\frac{{{a_i}(\theta ) - \left[ {{b_i}(\theta ) + {c_i}(\theta )} \right]\mathit{\Lambda }\left( {{t_i};\theta } \right) + {d_i}(\theta ){{\left[ {\mathit{\Lambda }\left( {{t_i};\theta } \right)} \right]}^2}}}{{\sigma _{{\rm{B}},n}^{{2^{(m + 1)}}}\left( {{t_i} - {t_{i - 1}}} \right)}} + \frac{{{d_i}(\theta )}}{{\sigma _{\gamma ,n}^{{2^{(m + 1)}}}}}} \right\}} } \right\} \end{array} $

（29）

由此可见，与传统的EM算法相比，本文的ECM算法通过CM1和CM2两步，可以将原本需要优化的3个参数简化为一个参数θ，降低了计算复杂度，提高了在线性能。此外，本文提出的参数估计方法可以在监测到一个新的性能衰减数据时，实现对模型参数的自适应估计和在线更新，克服了新研发航空发动机历史数据和先验信息不足的限制，这对于工程实际具有重要意义。

4 实例验证

根据相关文献可知，航空发动机的性能衰减和故障通常表现在排气温度裕度(EGTM)的衰减^[22-23]。因此，EGTM是表征航空发动机健康状态的一个重要指标。EGTM的定义为航空发动机出厂前测定的排气温度(EGT)限制值与全推力(或全功率)起飞时实际的EGT之间的差值。随着航空发动机工作时间的增加，EGT会逐渐升高，反映到EGTM就会逐渐衰退，因此准确预测EGTM对提高航空发动机的可靠性至关重要。

实际中，由于外部环境和内部因素的随机性，航空发动机在某一时刻的性能衰减量也是随机的，而且同类型航空发动机EGTM的衰减轨迹虽然具有相似性，但实际上不同个体之间仍然存在着差异性。因此，监测的EGTM数据不仅会受到时变不确定性的影响，还会受到个体差异不确定性的影响。此外，EGT数据来源于航空发动机上安装的传感器，经过空地数据链、报文解析程序等步骤进入到性能监测系统，在这整个过程中不可避免地会引入测量误差，造成测量的不确定性。根据EGTM的定义可知，EGTM是基于EGT计算的，因此EGTM数据中必然存在测量误差。通过上述分析发现，监测的EGTM数据中存在着三源不确定性共存的现象，因此同时考虑三源不确定性是必要的，也是更加符合实际需求的。下面以某型号航空发动机在6年内监测的EGTM数据来验证本文方法的有效性，收集的性能衰减数据如图 1所示。根据工程经验，设定EGTM的失效阈值ω=0 ℃，并且考虑以下幂函数形式的性能衰减模型：

$ X(t) = {x_0} + \lambda \cdot {t^\theta } + {\sigma _{\rm{B}}}B(t) $	（30）

这里，根据EGTM数据的特点，选择式(1)中非线性函数的形式为ρ(t; θ)=θ·t^θ-1。到目前为止，此类模型已经在航空发动机、轴承、电池等复杂系统的性能衰减建模中得到了广泛的应用^{[3, 12, 17]}。

如图 1所示，监测间隔为100次循环，在监测到第4 000次循环时，EGTM首次达到失效阈值ω，期间共监测了41次。显然，随着循环次数的增加，航空发动机的性能衰减路径整体呈下降趋势。然而，本文所提出的剩余寿命预测方法是基于性能衰减路径呈上升趋势的情况下推导出来的，因此需要对原始性能衰减数据进行转换。这里，本文利用EGTM的初始值减去其他所有的性能衰减数据，相应地，失效阈值转换为初始值减去ω。下面，基于转换后的性能衰减数据来验证本文方法的有效性和优越性。随着性能衰减数据的积累，采用所提出的参数估计方法可以对模型参数Θ在每一个监测点进行自适应地估计和更新，结果如图 2所示。

图 2的结果表明，模型参数可以随着性能衰减数据的积累很快地收敛，并且在获得一个新的性能衰减数据时，模型参数和剩余寿命的PDF都可以进行自适应地估计和更新。因此，克服了历史数据和先验信息不足的限制，能够适用于新研发航空发动机的剩余寿命预测。

为了进一步证明本文方法的优越性，选择以下2种方法进行对比研究：

1) M₁^[14]：同时考虑了三源不确定性，但在参数估计中，需要假定存在多组同类型设备的历史数据来离线地估计模型的初始参数，并且这些参数一旦确定，就不再随着性能衰减数据的积累进行更新。

2) M₂^[17]：在参数估计中，模型参数可以实现在线更新。但在性能衰减建模中，存在与随机参数的更新机制相矛盾的潜在假设，并且仅考虑了时变不确定性和个体差异不确定性，忽略了实际中普遍存在的测量不确定性。

图 3给出了3种方法在最后6个监测点预测的剩余寿命的PDF曲线，从图中可以看出，3种方法得到的PDF曲线均包含了实际的剩余寿命，且随着循环次数的增加逐渐变窄，这意味着预测的剩余寿命的不确定性越来越低。此外，由于本文方法在性能衰减建模中同时考虑了三源不确定性，并且在每一个监测点可以对模型参数进行自适应地估计和在线更新，因此本文方法得到的PDF曲线更加窄而尖锐，说明本文方法在剩余寿命预测方面具有更高的精度和更好的性能。

为了进一步量化比较结果，本文采用可靠性领域中常用的均方误差(MSE)指标，该指标同时考虑了剩余寿命预测的精度和不确定性，在t_n时刻的定义式为

$ {v_{{\rm{MSE}},n}} = \int_0^\infty {{{\left( {{l_n} - {{\tilde l}_n}} \right)}^2}} {f_{{L_n}|{Y_{0,n}}}}\left( {{l_n}|{Y_{0:n}}} \right){\rm{d}}{l_n} $	（31）

式中：$\tilde{l}_{n}$为t_n时刻航空发动机实际的剩余寿命，f_Ln|Y0:n(l_n|Y_0:n)为剩余寿命的PDF。MSE越小，说明预测的结果越准确。

图 4给出了3种方法在所有监测点的MSE比较结果。由于本文方法和M₂方法均不依赖于历史数据和先验信息，且模型的初始参数是随机给定的，因此在监测初期性能衰减数据较少时，剩余寿命预测的MSE高于M₁，但随着性能衰减数据的积累，本文方法的MSE明显低于M₁和M₂。为进一步验证，图 5给出了4个不同监测点的箱形图。其中，箱子的大小表示预测的剩余寿命的不确定性，若箱子越大，则预测的不确定性越大，精度越低；点划线表示当前监测点的实际剩余寿命；箱子中的实线表示预测的剩余寿命的期望，实线越接近点划线，说明预测的精度越高，结果越准确。从图 5(a)中可以看出，此时由于性能衰减数据较少，M₁优于本文方法和M₂，但是随着性能衰减数据的积累，本文方法的箱子更小且更加接近实际的剩余寿命，如图 5(b)~图 5(d)所示。此外，引入相对误差(RE)指标来进一步量化剩余寿命预测的精度，在寿命的20%、45%、70%和95%分位点上分别给出了3种方法RE的比较结果，如表 1所示。

通过表 1的结果可以看出，本文方法可以有效减少剩余寿命预测的RE，进而提高了剩余寿命预测的精度，尤其在寿命的95%分位点处，剩余寿命预测的RE仅为1.68%，明显小于其他2种方法。以上的实验结果表明，本文方法具有更好的预测能力，剩余寿命预测的精度更高，且不依赖于历史数据和先验信息，更加适用于新研发的航空发动机。

5 结论

1) 本文在三源不确定性下建立了非线性的性能衰减模型，不仅克服了现有方法的潜在假设，而且可以较好地建模航空发动机的性能衰减轨迹，并能够将三源不确定性的影响纳入到预测的剩余寿命分布中。

2) 利用Kalman滤波技术可以对潜在的性能衰减状态进行实时估计，并基于首达时间的概念，得到了同时考虑三源不确定性的剩余寿命分布。

3) 在获得一个新的性能衰减数据时，基于RTS平滑和ECM算法对模型参数进行自适应地估计和在线更新，克服了新研发航空发动机缺乏历史数据和先验信息的问题，降低了参数估计的不确定性，显著提高了剩余寿命预测的精度，尤其在寿命的95%分位点处，剩余寿命预测的RE仅为1.68%。