基于状态空间的惯性测量组合剩余寿命在线预测

引用本文

冯磊, 王宏力, 周志杰, 司小胜, 邹红星. 基于状态空间的惯性测量组合剩余寿命在线预测. 清华大学学报: 自然科学版, 2014, 54(4): 508-514[Lei FENG, Hongli WANG, Zhijie ZHOU, et al. Residual life prediction based on the state space for inertial measurement units[J]. 2014, 54(4): 508-514] 复制到剪切板

Permissions

基于状态空间的惯性测量组合剩余寿命在线预测

冯磊^1,², 王宏力², 周志杰², 司小胜^1,², 邹红星¹

1. 清华大学自动化系, 北京 100084

2. 第二炮兵工程大学, 西安 710025

邹红星, 教授, E-mail:hongxing_zou@mail.tsinghua.edu.cn

作者简介: 冯磊(1983-), 男(汉), 河北, 博士研究生。

基金:国家自然科学基金面上项目(61174030); 国家青年科学基金项目(61004069);

摘要

实时准确的剩余寿命预测能够为惯性测量组合的维护策略安排提供有效的决策支持。由于反映惯性测量组合退化状态的性能指标不能直接监测或直接测量带有噪声,因此需要构建状态空间模型预测惯性测量组合的剩余寿命。考虑到惯性测量组合的性能退化指标随时间呈现非线性特征,首先采用基于非线性漂移的Brown运动(Brownian motion, BM)建模其退化状态,然后基于构建的状态空间模型,利用期望最大化(expectation-maximization, EM)算法和Kalman滤波(Kalman filter)实时估计和更新退化状态和模型未知参数。并且将状态估计的分布函数引入剩余寿命的预测过程,近似得到了剩余寿命分布的解析形式,实现了剩余寿命的实时预测与更新。最后,对惯性测量组合的剩余寿命实时预测问题进行了实验分析,结果表明该方法具有较高的预测精度与较小的预测不确定性。

关键词: 剩余寿命; 预测; 状态空间; 非线性; 期望最大化

中图分类号:TP274⁺.3 文献标志码:A 文章编号:1000-0054(2014)04-0508-07

Residual life prediction based on the state space for inertial measurement units

Lei FENG^1,², Hongli WANG², Zhijie ZHOU², Xiaosheng SI^1,², Hongxing ZOU¹

1. Department of Automation, Tsinghua University, Beijing 100084, China

2. The Second Artillery Engineering University, Xi'an 710025, China

Fund:

Abstract

Real-time and accurate residual life prediction for an inertial measurement unit (IMU) can provide effective decision support for timely and cost-effective maintenance scheduling. The performance index reflecting the degradation of the IMU cannot be observed directly and direct measurements are contaminated by noise. Thus, a state space model was developed to predict the residual life of an IMU. Since the changes in the degradation state of the IMU are nonlinear over time, this analysis was a nonlinear drift-driven Brownian motion (BM) is used to characterize the degradation process, with the expectation maximization (EM) algorithm and the Kalman filter used to jointly estimate and update the state and model parameters. Furthermore, the estimated state distribution is incorporated into the residual life model using an approximate analytical form of the distribution. The approach is validated by comparison with experimental data which indicates that this method gives better prediction accuracies and lower uncertainties.

Keyword: residual life; prediction; state space; nonlinearity; expectation maximization

Show Figures

随着信息化和自动化水平的不断提高,许多工程设备、特别是航空航天设备与军事装备日益复杂。这些设备一旦发生失效将带来灾难性的后果。为了保证设备能够正常运行,避免因为失效引起的经济损失、人员伤亡、任务延误和环境破坏,必须对设备的健康状态即剩余寿命信息进行及时的评估与预测,据此安排必要的维护操作^[1]。捷联惯性导航系统(strapdown inertial navigation system, SINS)在军事装备、航空航天系统等领域有着重要的应用。其通过惯性测量组合(inertial measurement unit, IMU)测量运载体的角速度和线加速度,经过坐标变换与计算,得到运载体的位置、速度与姿态角信息,最终实现自主导航。由于SINS的核心部件IMU直接固定在运载体上,导致陀螺仪和加速度计直接承受运载体的振动、冲击和角运动,工作环境非常恶劣。随着工作时间的积累, IMU的性能会不可避免地发生退化,在这种情况下,及时有效地预测IMU的剩余寿命,适时采取必要的维护手段,对于提高整个系统的可靠性与安全性有着至关重要的作用。

现有的预测方法中,一般假设设备的性能退化变量直接可观测,通过对设备退化数据的监测,建立其退化规律的演化模型,进而预测设备的剩余寿命^[2]。但是,工程实际中,随着设备复杂性的提高,直接测量的退化变量因受到各种噪声与误差的影响,不能精确反映设备的退化状态,这种测量被称为不完整测量^[3,4]。文[4]指出,在测试IMU的陀螺仪漂移系数的过程中,由于受到环境因素、人为因素以及各种随机因素的影响,测量结果不可避免地存在噪声干扰与误差。为解决不完整测量的问题,文[4]利用状态空间模型,对设备实际的退化过程以及退化过程与不完整测量之间的关系建模,通过滤波技术预测设备的剩余寿命。此外,根据状态空间模型的预测方程与更新方程,可以方便地实现在线预测。但是,文[4]仅考虑了退化过程线性变化的情况。现有研究中大部分也都假设设备退化状态的期望即平均退化率随时间线性变化 $\begin{matrix} ^{[58]} \end{matrix}$ ,对于非线性退化过程,由于很难得到其剩余寿命分布的解析形式,研究相对较少,现有方法大多只能通过Monte-Carlo仿真得到寿命分布的数值解^[9,10],计算耗时,不易实现在线预测。如果可以得到剩余寿命分布的解析形式,就可以实现快速的在线预测,从而为后续的维护策略安排提供及时的决策支持。

IMU上的陀螺仪结构复杂,测量得到的漂移系数具有明显的非线性特征^[11],如果利用线性模型表征其漂移系数的变化过程,后续的剩余寿命预测会具有较大的误差与不确定性。本文将IMU上的陀螺仪的漂移系数作为性能退化变量,利用非线性漂移驱动的Brown运动(BM)表征其非线性特征,构建状态空间模型建立不完整测量与实际退化状态的关系,并且利用期望最大化(expectation maximization, EM)算法与Kalman滤波(Kalman filter)联合估计退化状态与未知参数,然后将估计的退化状态分布引入剩余寿命预测过程中,近似得到剩余寿命概率密度函数的参数化解析形式。根据滤波算法的预测方程与更新方程,一旦获取新的观测值,即可对状态与参数更新。相应地,可以更新剩余寿命的分布,最终实现在线剩余寿命预测。实验结果表明了本文所提方法的有效性。

1 基于状态空间模型的剩余寿命预测

1.1 退化模型的构建

假设 X( t)表示时刻 t设备的实际退化状态,对其采用如下非线性漂移驱动的BM建模:

$\begin{matrix} X (t) = X (0) + a t^{b} + σ_{B} B (t) . (1) \end{matrix}$

其中: X(0)表示初始退化状态,假设为0; at^b和 σ_B分别为漂移项和扩散系数, B( t)表示标准BM。可以看出,当 b=1时,式(1)退化为线性漂移驱动的BM, 说明线性情况只是该模型的特例。利用式(1)可以较好地表征退化过程非线性特征^[11,12,13]。

不完整测量广泛存在于工程实际中,为建模直接测量与实际退化之间的关系,采用如下加性噪声模型构建观测方程:

$\begin{matrix} Y (t) = X (t) + ε (t) . (2) \end{matrix}$

其中 ε( t)表示噪声项。假设 ε( t) ~N(0, σ²),且与 B( t)相互独立。

设备的寿命 T可以定义为退化状态首次达到预先设定的失效阈值的时间^[14],即 T可表示为

$\begin{matrix} T = \inf {t : X (t) \geq ω | X (s) < ω, 0 \leq s < t} . (3) \end{matrix}$

其中 ω表示失效阈值。 ω通常根据专家经验或领域标准预先给定。

根据式(3)可知,设备的寿命分布等价于首达时间分布,利用线性漂移驱动的BM建模退化过程,可直接得到首达时间分布的解析形式即逆Gauss分布^[2]。但是,对于非线性漂移驱动的扩散过程而言,其首达时间分布需要求解具有边界约束的FPK(Fokker-Planck-Kolmogorov)程。该方程的求解非常困难,很难得到分布的解析形式。为了简化问题,一些研究人员假设实际的退化过程是严格单调增加的,则寿命分布可直接表示为

$\begin{matrix} F_{T} (t) = \Pr (T \leq t) = \Pr (X (t) \geq ω) . (4) \end{matrix}$

式(4)只适用于退化过程严格单调的情况,对于受各种因素影响具有随机波动的非线性退化过程,利用随机扩散过程描述其特征更加符合工程实际。但是,对于基于BM的模型而言,式(4)的简化计算不适用。或者在这种强假设下,只能得到预测误差较大的近似结果。本文采用文[11]提出的一种弱假设解决上述情况下的寿命分布求解问题。

1.2 剩余寿命的预测

假设退化过程如果在时刻 t到达 ω,则 t之前退化过程到达 ω的概率忽略不计。

上述假设意味着退化过程可以在 t时刻以前到达 ω,也可以在 t之后返回 ω以下,但是出现这种情况的概率很小。因此,该假设不同于单调性假设,具体的解释与实验验证,可参见文[11]。下面简要论述在该假设下,寿命分布解析形式求解的一般思路。

步骤1 问题转化。

首先对退化过程{ X( t), t≥0}进行空间—时间变换^[14] $\begin{matrix} {\dot{X}}_{t} \end{matrix}$ =ψ( t, X_t), $\begin{matrix} \dot{t} \end{matrix}$ =φ( t),将求解非线性退化过程首达 ω的问题转化为标准BM首达依赖于 t的边界函数 S_B( t)的首达时间问题,其中 S_B( t) =( ω-at^b) /σ_B . 变换后, T的概率密度函数 f_T( t)可表示为

$\begin{matrix} f_{T} (t) = p_{B (t)} (S_{B} (t), t) \frac{dφ (t)}{dt} . (5) \end{matrix}$

其中 p_B₍ _t₎( S_B( t), t)为标准BM首达依赖于 t的边界函数 S_B( t)的首达时间概率密度函数。

步骤2 p_B₍ _t₎( S_B( t), t)的求解。

为了得到 p_B₍ _t₎( S_B( t), t)的形式,利用标准BM的Markov性、独立增量性和均方连续性,可以将其进一步分解为两部分,一部分表示 B( t)在边界 S_B( t)上的概率密度,另一部分用来描述退化过程在时刻 t到达边界上的情况下,在时刻 s<t时没有到达边界的概率。通过求解两部分的概率,可以实现对 p_B₍ _t₎( S_B( t), t)的求解。

根据上述求解步骤, f_T( t)可以表示为

$\begin{matrix} f_{T} (t) \approx \frac{1}{\sqrt[]{2 πt}} (\frac{S_{B} (t)}{t} + \frac{1}{σ_{B}} ab t^{b - 1}) \exp (- \frac{S_{B}^{2} (t)}{2 t}) . (6) \end{matrix}$

得到寿命的概率密度函数后,接下来推导如何利用当前时刻 t_i的退化状态 X( t_i)得到当前时刻剩余寿命的概率密度函数。对于 t>t_i,退化过程可表示为 X( t) =X( t_i) +at^b-a $\begin{matrix} t_{i}^{b} \end{matrix}$ +σ_B B( t-t_i) . 则当前时刻,设备的剩余寿命为

$\begin{matrix} L_{i} = \inf {l_{i} : X (t_{i} + l_{i}) \geq ω | X (t_{i}), l_{i} \geq 0} . (7) \end{matrix}$

令 l_i=t-t_i,可以得到一个新的退化过程{ G( l_i), l_i≥0},

$\begin{matrix} G (l_{i}) = G (0) + a (l_{i} + t_{i})^{b} - a t_{i}^{b} + σ_{B} B (l_{i}), l_{i} \geq 0 . (8) \end{matrix}$

其中 G(0) =0 .

根据上述方程可知,时刻 t_i的剩余寿命等价于{ G( l_i), l_i≥0}首达阈值 $\begin{matrix} ω_{t_{i}} \end{matrix}$ =ω-X( t_i)的时间。但是,实际的退化状态是隐含的,难以得到其精确值,因此,只能通过构造的状态空间模型估计时刻 t_i的状态 X( t_i) .为了辨识实际的退化状态,首先将动态系统离散化,得到离散时间点 t_k=kΔ t, k=1,2,…的状态方程和观测方程:

$\begin{matrix} \{\begin{matrix} X_{k} = X_{k - 1} + a t_{k}^{b} - a t_{k - 1}^{b} + σ_{B} \sqrt[]{Δt} ω_{k}, \\ Y_{k} = X_{k} + σ υ_{k} . \end{matrix} (9) \end{matrix}$

其中: Δ t为离散化步长, $\begin{matrix} (ω_{k})_{k \geq 1} \end{matrix}$ 和 $\begin{matrix} (υ_{k})_{k \geq 1} \end{matrix}$ 为独立同分布的噪声序列,且 ω_k~N(0,1), υ_k~N(0,1) .

根据状态空间模型式(9), 可以通过Kalman滤波^[16]估计当前时刻实际的退化状态。首先定义 $\begin{matrix} {\hat{X}}_{i | i} \end{matrix}$ =E( X_i|Y_1: _i)和 P_i|i=Var( X_i|Y_1: _i)分别为利用当前时刻的整个观测序列 Y_1: _i=( Y₁, Y₂,…, Y_i)得到的状态 X_i的条件期望和方差,即 X_i~N( $\begin{matrix} {\hat{X}}_{i | i} \end{matrix}$ , P_i|i)。根据全概率公式,可以得到当前时刻 t_i,剩余寿命概率密度函数为

$\begin{matrix} \begin{matrix} f_{L_{i}} (l_{i} | θ, Y_{1 : i}) \approx \frac{1}{\sqrt[]{2 π l_{i}^{2} (P_{i | i} + σ_{B}^{2} l_{i})}} \\ [ω - aλ (l_{i}) + ab l_{i} (l_{i} + t_{i})^{b - 1} - \\ \frac{P_{i | i} (ω - aλ (l_{i})) + {\hat{X}}_{i | i} σ_{B}^{2} l_{i}}{P_{i | i} + σ_{B}^{2} l_{i}}] \\ \exp (- \frac{(ω - aλ (l_{i}) - {\hat{X}}_{i | i})^{2}}{2 (P_{i | i} + σ_{B}^{2} l_{i})}) . (10) \end{matrix} \end{matrix}$

其中: λ( l_i) = $\begin{matrix} (l_{i} + t_{i})^{b} \end{matrix}$ - $\begin{matrix} t_{i}^{b} \end{matrix}$ , θ=[ a, b, σ_B, σ]表示未知参数向量。

2 参数估计算法

节1推导了时刻 t_i剩余寿命的概率密度函数的解析形式。接下来,采用极大似然估计(maximum likelihood estimation, MLE)策略估计模型的未知参数。

为后续表示方便,首先定义当前时刻的状态序列为 X_1: _i=( X₁, X₂,…, X_i)。由于构建的状态空间模型含有未知状态,不能直接采用MLE估计未知参数,本文采用EM算法解决存在未知状态时的参数极大似然估计问题。通过迭代地计算且最大化包含完整数据集( X_1: _i, Y_1: _i)的对数似然函数的条件期望, EM算法能够产生一个估计序列收敛到参数的极大似然估计值。令 $\begin{matrix} {\hat{θ}}^{(j)} \end{matrix}$ 表示EM算法第 j次迭代所得的估计值,则完整数据对数似然函数的条件期望可表示为

$\begin{matrix} Q (θ, {\hat{θ}}^{(j)}) = E_{{\hat{θ}}^{(j)}} ([lnp (X_{1 : i}, Y_{1 : i} | θ)] | Y_{1 : i}) . (11) \end{matrix}$

其中ln p( X_1: _i, Y_1: _i|θ)为完整数据对数似然函数。总结起来,参数估计程序包含如下两个步骤:

1) E步。

$\begin{matrix} 计算 : Q (θ, {\hat{θ}}^{(j)}) \end{matrix}$

2) M步。

$\begin{matrix} 求解 : {\hat{θ}}^{(j + 1)} = argma x_{θ} Q (θ, {\hat{θ}}^{(j)}) \end{matrix}$

给定初始值 $\begin{matrix} {\hat{θ}}^{(0)} \end{matrix}$ ,上述步骤迭代计算直到满足收敛标准 $\begin{matrix} ^{[1718]} \end{matrix}$ 。为提高参数估计精度,初值 $\begin{matrix} {\hat{θ}}^{(0)} \end{matrix}$ 可根据同类设备的历史失效数据确定。

具体到本文构建的状态空间模型式(9), 完整数据的对数似然函数可以表示为

$\begin{matrix} \begin{matrix} lnp (X_{1 : i}, Y_{1 : i} | θ) = \ln [\overset{i}{\prod_{k = 1}} p (X_{k} | X_{k - 1}; θ) \\ \overset{i}{\prod_{k = 1}} p (Y_{k} | X_{k}; θ)] = - iln (2 π \sqrt[]{Δt}) - iln σ_{B} - \\ ilnσ - \frac{1}{2 σ^{2}} \overset{i}{\sum_{k = 1}} (_{k}^{Y} - X_{k}) 2 - \frac{1}{2 σ_{B}^{2} Δt} \overset{i}{\sum_{k = 1}} [X_{k} - \\ X_{k - 1} - a ({(kΔt)}^{b} - ({(k - 1) Δt)}^{b})]^{2} . (12) \end{matrix} \end{matrix}$

为了计算对数似然函数的条件期望,对于 k=1,2,…, i,首先定义如下变量:

$\begin{matrix} \begin{matrix} {\hat{X}}_{k | i} = E_{{\hat{θ}}^{(j)}} (X_{k} | Y_{1 : i}), \\ P_{k | i} = E_{{\hat{θ}}^{(j)}} (X_{k}^{2} | Y_{1 : i}) - {\hat{X}}_{k | i}^{2}, \\ P_{k, k - 1 | i} = E_{{\hat{θ}}^{(j)}} (X_{k} X_{k - 1} | Y_{1 : i}) - {\hat{X}}_{k | i} {\hat{X}}_{k - 1 | i} . \end{matrix} \end{matrix}$

经过一系列代数运算,可以得到

$\begin{matrix} \begin{matrix} Q (θ, {\hat{θ}}^{(j)}) \propto - \frac{1}{2 σ_{B}^{2} Δt} \overset{i}{\sum_{k = 1}} [A_{k} - 2 a ({(kΔt)}^{b} - \\ ({(k - 1) Δt)}^{b}) B_{k} + a^{2} ({(kΔt)}^{b} - ({(k - 1) Δt)}^{b})^{2}] - \\ iln σ_{B} - ilnσ - \frac{1}{2 σ^{2}} \overset{i}{\sum_{k = 1}} (Y_{k}^{2} + P_{k | i} + {\hat{X}}_{k | i}^{2} - 2 Y_{k} {\hat{X}}_{k | i}) . (13) \end{matrix} \end{matrix}$

其中 A_k=P_k|i+ $\begin{matrix} {\hat{X}}_{k | i}^{2} \end{matrix}$ +P_k-₁ _|i+ $\begin{matrix} {\hat{X}}_{k - 1 | i}^{2} \end{matrix}$ -2 P_k_, _k-₁ _|i-2 $\begin{matrix} {\hat{X}}_{k | i} \end{matrix} \begin{matrix} {\hat{X}}_{k - 1 | i} \end{matrix}$ , B_k= $\begin{matrix} {\hat{X}}_{k | i} \end{matrix}$ - $\begin{matrix} {\hat{X}}_{k - 1 | i} \end{matrix}$ 。

显然,计算完整数据对数似然函数需要得到相对于 $\begin{matrix} {\hat{θ}}^{(j)} \end{matrix}$ 的估计 $\begin{matrix} {\hat{X}}_{k | i} \end{matrix}$ 、 $\begin{matrix} {\hat{X}}_{k - 1 | i} \end{matrix}$ 、 P_k|i、 P_k-₁ _|i和 P_k_, _k-₁ _|i。本文采用Kalman平滑(Kalman smoother)计算这些未知量。Kalman平滑包括两部分: 前向递归即滤波,后向递归即平滑。具体细节可参见文[15]。

根据式(13)可知,目标函数含有4个未知参数,为了降低参数空间的维数,减少计算复杂度,本文采用剖面似然函数技术估计未知参数。首先假设参数 b固定,分别计算式(13)相对于参数 a、 σ_B和 σ的偏导数,并且令这3个偏导数为零,可得

$\begin{matrix} \hat{a} = \frac{\sum_{k}^{=} ({(kΔt)}^{b} - ({(k - 1) Δt)}^{b}) B_{k}}{\sum_{k}^{=} {({(kΔt)}^{b} - ({(k - 1) Δt)}^{b})}^{2}}, (14) \end{matrix}$

$\begin{matrix} {\hat{σ}}_{B}^{2} \end{matrix}$ = $\begin{matrix} \frac{\overset{i}{\sum_{k = 1}} [A_{k} - 2 a ({(kΔt)}^{b} - ({(k - 1) Δt)}^{b}) B_{k} + a^{2} {({(kΔt)}^{b} - ({(k - 1) Δt)}^{b})}^{2}]}{iΔt} \end{matrix}$ ,(15)

$\begin{matrix} {\hat{σ}}^{2} = \frac{\sum_{k}^{=} (Y_{k}^{2} + P_{k | i} + {\hat{X}}_{k | i}^{2} - 2 Y_{k} {\hat{X}}_{k | i})}{i} . (16) \end{matrix}$

其中 $\begin{matrix} \hat{a} \end{matrix}$ 、 $\begin{matrix} {\hat{σ}}_{B} \end{matrix}$ 和 $\begin{matrix} \hat{σ} \end{matrix}$ 分别为 a、 σ_B和 σ的估计值。

将式(14)和(15)代入式(13), 可以得到只含有 b的剖面似然函数。通过一维搜索使此剖面似然函数最大化,即可得到 b的估计值 $\begin{matrix} \hat{b} \end{matrix}$ 。将 $\begin{matrix} \hat{b} \end{matrix}$ 作为 b的值分别代入式(14)和(15), 即可得到 $\begin{matrix} \hat{a} \end{matrix}$ 和 $\begin{matrix} {\hat{σ}}_{B} \end{matrix}$ 。 $\begin{matrix} \hat{σ} \end{matrix}$ 可直接通过计算式(16)得到。

根据上述推导,可以看出,一旦获取新的观测值,首先利用EM算法估计与更新模型的未知参数,然后利用Kalman滤波更新退化状态。参数与状态的更新逐次进行。相应地,根据式(10), 可以更新剩余寿命的概率分布,实现实时剩余寿命预测。

3 惯性测量组合剩余寿命预测实验

3.1 问题描述

IMU是SINS的关键设备,其运行状态对导航精度有着重要的影响。IMU的核心部件为3个陀螺仪和3个加速度计,分别测量运载体的角速度和线加速度。工程实际中,陀螺仪工作在不同环境与负载下,转子高速旋转必然造成转轴的磨损,随着运行时间的增长,磨损不断累积,引起漂移系数的增长,导致系统性能逐渐退化,当漂移系数增大到预先设定的阈值,即认为系统发生失效^[19,20]。因此,陀螺仪漂移是表征IMU性能的一项重要指标,是一种缓变失效。由于陀螺仪敏感轴方向上的漂移系数的变化趋势较为明显,因此本实验选取敏感轴方向的一次项漂移系数作为衡量IMU退化状态的性能指标,记为 K₁。 K₁由每小时的漂移度数表示,单位记为°/h。根据工程实践经验和对惯性测量组合精度的要求,通常假设 K₁的值一旦达到失效阈值0 .38°/h,则表明陀螺仪由于漂移引起的性能退化非常严重,这种情况下,可以认为陀螺仪已经失效。实验中收集的是对某IMU连续工作测试得到的一次项漂移系数。采样间隔为3 h, 共96组数据,如图1所示。

	Figure Option View Download New Window Download As Powerpoint Slide
	图1 一次项漂移系数测试数据

为了便于比较研究,本文采用每一个观测点处预测的剩余寿命分布的均方误差(mean-squared error, MSE)为标准,比较不同模型的预测精度。具体地, MSE定义如下^[21,22,23]:

$\begin{matrix} MS E_{i} = \int_{0}^{\infty} ({\dot{l}}_{i} - l_{i})^{2} f_{L_{i}} (l_{i} | θ, Y_{1 : i}) d l_{i} . (17) \end{matrix}$

其中 $\begin{matrix} {\dot{l}}_{i} \end{matrix}$ 为时刻 t_i设备实际的剩余寿命, $\begin{matrix} f_{L_{i}} \end{matrix}$ ( l_i|θ, Y_1: _i)为每一个状态监测点处根据式(10)得到的剩余寿命概率密度函数。相应地,总均方误差(total mean-squared error, TMSE)可以表示如下:

$\begin{matrix} TMSE = \overset{N}{\sum_{i = 1}} \int_{0}^{\infty} ({\dot{l}}_{i} - l_{i})^{2} f_{L_{i}} (l_{i} | θ, Y_{1 : i}) d l_{i} . (18) \end{matrix}$

其中 N为观测数据的个数。注意到 $\begin{matrix} {\dot{l}}_{i} \end{matrix}$ 为时刻 t_i设备实际的剩余寿命,而本实验所得到的测试数据的第96个数据点为0 .367 7°/h,没有超过失效阈值,因此,IMU在每个监测时间点的真实剩余寿命是未知的。但是,根据测试数据可以推断,在第96次监测之后,漂移系数已经很接近失效阈值,为了便于比较研究,本文假设IMU的寿命近似为294 h。

3.2 实验结果

如前所述,本文利用EM算法和Kalman滤波对参数进行估计与更新。各未知参数的估计与更新结果如图2所示。

	Figure Option View Download New Window Download As Powerpoint Slide
	图2 参数估计与更新结果

由图2可看出,在不同的监测时刻,一旦获取新的观测值,则利用到当前时刻为止的整个观测序列对参数进行估计与更新。随着观测数据的不断增加,各参数的估计值最终收敛到各自的稳定值,收敛速度相对较快。最终 $\begin{matrix} \hat{a} \end{matrix}$ 、 $\begin{matrix} \hat{b} \end{matrix}$ 、 $\begin{matrix} {\hat{σ}}_{B} \end{matrix}$ 和 $\begin{matrix} \hat{σ} \end{matrix}$ 分别为0 .022 7、 0 .601 9、 0 .007 2和0 .004 8。此外,由 $\begin{matrix} \hat{b} \end{matrix}$ 可知,陀螺仪漂移确实具有非线性特征。

为了验证本文方法的有效性,与线性模型和单调性假设下的模型做比较研究。为表示方便,将本文所构建的模型称为模型1, 线性模型为模型2, 单调性假设下的模型为模型3。以第90到96个数据为例,剩余寿命预测结果如图3所示。

	Figure Option View Download New Window Download As Powerpoint Slide
	图3 三种模型的剩余寿命预测结果的比较

由图3可知,模型1由于在预测剩余寿命的过程中考虑了退化状态估计的不确定性,预测得到的剩余寿命的概率密度曲线比另外2种模型的更加紧致,预测结果的不确定性更小。为了进一步说明模型1的预测精度高,分别比较了这3种模型的MSE和TMSE, MSE的比较结果如图4所示。

	Figure Option View Download New Window Download As Powerpoint Slide
	图4 三种模型的剩余寿命预测结果的均方误差

当测试数据较少时,3种模型的MSE都很大,这是因为测试数据少导致参数和实际的退化状态估计不准确,因而预测的剩余寿命也不准确; 并且实际的剩余寿命未知, 294 h只是近似假设。根据式(17)可以看出,预测结果不准确将会导致MSE很大,但是相对于三种模型而言,这是同一标准下的公平比较。由图4可知,从初始时刻开始,模型1的MSE就比模型2的小,说明退化路径存在非线性特征时,如果利用模型2建模这种退化过程,结果有可能具有较大的误差,预测精度不高。此外,大约11个数据点之后,模型1的MSE都比模型3的小,说明如果退化过程非单调,利用单调性假设得到的预测结果比较保守,精度比在本文提出的弱假设下得到的结果要差。模型1、 2和3的预测结果的TMSE分别是3.968 0×10⁵、8.561 1×10⁵和1.855 0×10⁶。

总之,以上实验结果表明本文方法可以很好地解决IMU的剩余寿命预测问题,并且具有较高的预测精度和较小的预测不确定性。

4 结论

为了解决IMU剩余寿命中存在的不完整测量和非线性退化的问题,本文提出了一种基于状态空间的剩余寿命在线预测方法。该方法具有以下特点:

1) 构建了基于状态空间技术的惯性测量组合剩余寿命预测模型,该模型能够较好地建模退化状态以及不完整测量与实际退化状态之间的关系。

2) 利用非线性漂移驱动的BM建模工程实际中广泛存在的非线性退化过程。与单调性假设不同,本文在一种弱假设的条件下,充分考虑状态估计过程中的不确定性,并将其引入后续的剩余寿命预测过程,近似得到了剩余寿命概率密度函数参数化的解析形式。

3) 基于构建的状态空间预测模型,利用EM算法和Kalman滤波,实现了退化状态和参数的联合估计与更新。最终实现了剩余寿命的在线预测。

通过仿真实验分析表明: 本文所提方法无论在可操作性还是准确性、降低不确定性方面均要优于线性模型和单调性假设下的模型,符合工程实际的要求,能够较好地解决IMU剩余寿命预测的相关问题,并且可以为后续的维护策略安排提供准确而实时的预测信息。

The authors have declared that no competing interests exist.

参考文献

View Option

[1]	Zhou Z J, Hu C H, Xu D L, et al. A model for real-time failure prognosis based on hidden Markov model and belief rule base [J]. European Journal of Operational Research, 2010, 207(1): 269-283. [本文引用:1] [JCR: 1.843]
[2]	Si X S, Wang W, Hu C H, et al. Remaining useful life estimation-A review on the statistical data driven approaches[J]. European Journal of Operational Research, 2011, 213(1): 1-14. [本文引用:2] [JCR: 1.843]
[3]	Sun J Z, Zuo H F, Pecht M G. Advances in sequential Monte Carlo methods for joint state and parameter estimation applied to prognostics [C]// Prognostics and System Health Management Conference. Shenzhen, China: IEEE Press, 2011: 1-7. [本文引用:1]
[4]	冯磊, 王宏力, 韩晓霞等. 考虑监测噪声的陀螺仪剩余寿命在线预测[J]. 长春工业大学学报(自然科学版), 2013, 34(2): 155-159. FENG lei, WANG Hongli, HAN Xiaoxia, et al. Real-time prediction on gyro residual life considering the monitoring noise[J]. Journal of Changchun University of Technology (Natural Science Edition), 2013, 34(2): 155-159. (in Chinese) [本文引用:1] [CJCR: 0.331]
[5]	Xu Z G, Ji Y D, Zhou D H. Real-time reliability prediction for a dynamic system based on the hidden degradation process identification[J]. IEEE Transactions on Reliability, 2008, 57(2): 230-242. [本文引用:1] [JCR: 1.657]
[6]	Christer A H, Wang W, Sharp J M. A state space condition monitoring model for furnace erosion prediction and replacement[J]. European Journal of Operational Research, 1997, 101(1): 1-14. [本文引用:1] [JCR: 1.843]
[7]	Batzel T D, Swanson D C. Prognostic health management of aircraft power generators[J]. IEEE Transactions on Aerospace and Electronic Systems, 2009, 45(2): 473-483. [本文引用:1] [JCR: 1.394]
[8]	Gasperin M, Juricic D, Boskoski P. Model-based prognostics of gear health using stochastic dynamical models[J]. Mechanical Systems and Signal Processing, 2010, 25(2): 537-548. [本文引用:1] [JCR: 2.465]
[9]	Cadini F, Zio E, Avram D. Monte Carlo-based filtering for fatigue crack growth estimation[J]. Probabilistic Engineering Mechanics, 2009, 24(3): 367-373. [本文引用:1] [JCR: 1.46]
[10]	Zio E, Peloni G. Particle filtering prognostic estimation of the remaining useful life of nonlinear components[J]. Reliability Engineering and System Safety, 2011, 96(3): 403-409. [本文引用:1] [JCR: 2.048]
[11]	Si X S, Wang W, Hu C H, et al. Remaining useful life estimation based on a nonlinear diffusion degradation process[J]. IEEE Transactions on Reliability, 2012, 61(1): 50-67. [本文引用:2] [JCR: 1.657]
[12]	Zhou Y F, Sun Y, Mathew J, et al. Latent degradation indicators estimation and prediction: A Monte Carlo approach[J]. Mechanical Systems and Signal Processing, 2011, 25(1): 222-236. [本文引用:1] [JCR: 2.465]
[13]	Lei F, Wang H L, Si X S, et al. A state space based prognostic model for hidden and age-dependent nonlinear degradation process[J]. IEEE Transactions on Automation Science and Engineering, 2013, 10(4): 1072-1086. [本文引用:1] [JCR: 2.162]
[14]	Lee M L T, Whitmore G A. Threshold regression for survival analysis: Modeling event times by a stochastic process reaching a boundary[J]. Statistical Science, 2006, 21(4): 501-513. [本文引用:2] [JCR: 1.69]
[15]	Ricciardi L M. On the transformation of diffusion processes into the Wiener process[J]. Journal of Mathematical Analysis and Applications, 1976, 54(1): 185-199. [本文引用:1] [JCR: 1.119]
[16]	王志贤. 最优状态估计与系统辨识 [M]. 西安: 西北工业大学出版社, 2004. WANG Zhixian. Optimal State Estimation and System Identification [M]. Xi'an: Northwestern Polytechnical University Press, 2004. (in chinese) [本文引用:1]
[17]	Dempster A, Laird N, Rubin D. Maximum likelihood from incomplete data via the EM algorithm[J]. Journal of the Royal Statistical Society, Series B, 1977, 39(1): 1-38. [本文引用:1]
[18]	Schon T B. An explanation of the expectation maximization algorithm [Z]. Division of Automatic Control, Linkoping University, Linkoping, Sweden, Tech. Rep. LITH-ISY-R-2915, 2009. [本文引用:1]
[19]	周志杰. 置信规则库在线建模方法与故障预测 [D]. 西安: 第二炮兵工程学院, 2010. ZHOU Zhijie. Online modeling methods of belief rule base and failure prognosis [D]. Xi'an: The Second Artillery Engineering University, 2010. (in chinese) [本文引用:1]
[20]	胡昌华, 司小胜. 基于信度规则库的惯性平台健康状态参数在线估计[J]. 航空学报, 2010, 31(7): 1454-1465. HU Changhua, SI Xiaosheng. Real-time parameters estimation of belief rule base health condition for inertial platform[J]. Acta Aeronautica et Astronautica Sinica, 2010, 31(7): 1454-1465. (in chinese) [本文引用:1] [CJCR: 0.949]
[21]	Carr M J, Wang W. An approximate algorithm for prognostic modelling using condition monitoring information[J]. European Journal of Operational Research, 2011, 211(1): 90-96. [本文引用:1] [JCR: 1.843]
[22]	Si X S, Hu C H, Wang W, et al. An adaptive and nonlinear drift-based Wiener process for remaining useful life estimation [C]//Prognostics and Health Management Conference, Shenzhen, China: IEEE Press, 2011: 1-5. [本文引用:1]
[23]	Carr M J, Wang W. Modeling failure modes for residual life prediction using stochastic filtering theory[J]. IEEE Transactions on Reliability, 2010, 59(2): 90-96. [本文引用:1] [JCR: 1.657]

2010

1.843

0.0

... 为了保证设备能够正常运行,避免因为失效引起的经济损失、人员伤亡、任务延误和环境破坏,必须对设备的健康状态即剩余寿命信息进行及时的评估与预测,据此安排必要的维护操作^[1] ...

2011

1.843

0.0

... 现有的预测方法中,一般假设设备的性能退化变量直接可观测,通过对设备退化数据的监测,建立其退化规律的演化模型,进而预测设备的剩余寿命^[2] ...

... 根据式(3)可知,设备的寿命分布等价于首达时间分布,利用线性漂移驱动的BM建模退化过程,可直接得到首达时间分布的解析形式即逆Gauss分布^[2] ...

2011

0.0

... 但是,工程实际中,随着设备复杂性的提高,直接测量的退化变量因受到各种噪声与误差的影响,不能精确反映设备的退化状态,这种测量被称为不完整测量^[3-4] ...

2013

0.0

0.331

2008

1.657

0.0

1997

1.843

0.0

2009

1.394

0.0

2010

2.465

0.0

2009

1.46

0.0

... 现有研究中大部分也都假设设备退化状态的期望即平均退化率随时间线性变化 58,对于非线性退化过程,由于很难得到其剩余寿命分布的解析形式,研究相对较少,现有方法大多只能通过Monte-Carlo仿真得到寿命分布的数值解^[9-10],计算耗时,不易实现在线预测 ...

2011

2.048

0.0

2012

1.657

0.0

... IMU上的陀螺仪结构复杂,测量得到的漂移系数具有明显的非线性特征^[11],如果利用线性模型表征其漂移系数的变化过程,后续的剩余寿命预测会具有较大的误差与不确定性 ...

... 利用式(1)可以较好地表征退化过程非线性特征^[11-13] ...

2011

2.465

0.0

2013

2.162

0.0

... 利用式(1)可以较好地表征退化过程非线性特征^[11-13] ...

2006

1.69

0.0

... 设备的寿命T可以定义为退化状态首次达到预先设定的失效阈值的时间^[14],即T可表示为 ...

... 首先对退化过程{X(t),t≥0}进行空间—时间变换^[14]X˙t=ψ(t,X_t), t˙=φ(t),将求解非线性退化过程首达ω的问题转化为标准BM首达依赖于t的边界函数S_B(t)的首达时间问题,其中S_B(t)=(ω-at^b)/σ_B ...

1976

1.119

0.0

2004

0.0

... 根据状态空间模型式(9), 可以通过Kalman滤波^[16]估计当前时刻实际的退化状态 ...

1977

0.0

2009

0.0

2010

0.0

... 工程实际中,陀螺仪工作在不同环境与负载下,转子高速旋转必然造成转轴的磨损,随着运行时间的增长,磨损不断累积,引起漂移系数的增长,导致系统性能逐渐退化,当漂移系数增大到预先设定的阈值,即认为系统发生失效^[19-20] ...

2010

0.0

0.949

胡昌华, 司小胜. 基于信度规则库的惯性平台健康状态参数在线估计[J]. 航空学报, 2010, 31(7): 1454-1465.

Changhua

, SI

Xiaosheng.

Real-time parameters estimation of belief rule base health condition for inertial platform[J]. Acta Aeronautica et Astronautica Sinica, 2010, 31(7): 1454-1465. (in chinese)

A real-time and accurate health condition prediction for an inertial platform is essential for cost-effective and timely maintenance planning and scheduling. Due to the fact that the true health condition of the inertial platform cannot be observed directly, it is assumed that the observations of characteristic parameters are available from monitoring, and the characteristic parameters correlate with health condition of the inertial platform. In this article, a health condition prediction system for the inertial platform is established based on belief rule base (BRB), where the characteristic parameters of the inertial platform are used as the inputs of BRB system and the health condition of platform as the output consequence. To overcome the drawbacks of current parameter optimization algorithms for BRB and satisfy real-time prediction, a parameter estimation algorithm is investigated for online updating BRB prediction system based on the expectation maximization (EM) algorithm. When the new input-output information of system operation is available, the model parameter can be updated online. Real-time health condition prediction for the inertial platform system is validated using the established model and the algorithm under investigation. The experimental results show that the proposed method can implement online parameter estimation of health condition prediction for the inertial platform effectively. In addition, compared with offline parameter optimization method, the proposed method can generate better results in terms of prediction accuracy and operating time, and thus has great potential in engineering practice.

302 Unit, Xi’an Institute of Hi-Tech

实时准确的健康状态预测是规划惯性平台系统及时、经济的维修策略的关键技术。由于平台系统的健康状态是不能够直接观测的，假设平台系统的特征参数监测数据是可以获取到的，而且平台系统的健康状态与特征量是相关的。基于信度规则库(BRB)，以平台系统的状态监测特征参数作为BRB系统的输入，以平台的健康状态作为输出结果，组建了惯性平台健康状态预测系统。为了克服现有BRB参数优化方法的不足，实现实时状态预测，基于期望最大化(EM)算法，研究了健康状态预测系统的参数在线估计算法。该算法在获取系统新的输入输出信息后，就对参数进行更新。利用本文提出的方法对惯性平台系统的健康状态实时预测问题进行了实验分析，实验结果表明：该方法可以有效地实现惯性平台系统健康状态预测模型参数实时估计；与参数离线优化方法相比，该方法不仅在预测精度上，而且在运行时间上都具有明显的优势；在工程实际中有良好的应用潜力。

2011

1.843

0.0

... 具体地, MSE定义如下^[21-23]: ...

2011

0.0

2010

1.657

0.0

... 具体地, MSE定义如下^[21-23]: ...