依赖于时间和状态的长程相关退化过程剩余寿命预测方法

摘要

本发明公开了一种依赖于时间和状态的长程相关退化过程剩余寿命预测方法，属于健康管理技术领域，包括以下步骤：首先，收集等间隔采样的传感器数据；根据退化数据特征构造基于分形布朗运动的退化模型；利用基于二次变分的方法估计模型中的Hurst指数；最大化利用Radon‑Nikodym导数构造的似然比函数来估计漂移项未知参数；利用极大似然方法估计扩散项未知参数；然后通过弱收敛理论将原来的退化过程近似为一个具有时变扩散项系数的基于布朗运动的随机过程；通过一组变换将退化过程进一步简化；最后得到解析的剩余寿命分布。本发明可较为准确地预测剩余寿命分布。

说明书

技术领域

本发明属于健康管理的技术领域，具体涉及一种依赖于时间和状态的长程相关退化过程剩余寿命预测方法。

背景技术

工业设备的正常运行是保障生产安全和企业效益的根本前提。工业设备一旦失效，不仅会给生产带来巨大的经济损失，还有可能引发严重的事故。因此，健康管理对工业设备具有重要的意义。剩余寿命预测是健康管理的重要组成部分。准确的剩余寿命预测结果能够有效指导设备的维护策略和备件供应，从而避免过度维修带来的浪费以及故障发生带来的损失。

为了获得准确的寿命预测结果，需要建立一个能够尽可能完全描述设备退化的模型。目前，大多数退化模型都是仅依赖于时间的非长程相关过程。然而，在实际设备的退化过程中，退化过程不仅依赖与时间，并且依赖于设备当前的退化状态。另一方面，环境影响、人为操作以及退化机制本身都有可能使得退化过程具有长程相关性，即不同时刻的退化之间存在相关性。虽然，目前有少量的工作关注于同时依赖于时间和状态的退化过程建模以及长程相关退化过程建模，但是并没有一种能够同时考虑依赖于时间和状态的长程相关退化过程建模及剩余寿命预测方法。

依赖于时间和状态的长程相关退化过程的剩余寿命预测主要有以下难点：第一，由于退化状态本身包含在退化模型中，因此无法写出解析的似然函数，难以使用传统估计方法进行参数估计。第二，由于退化过程不是Markov过程也不是半鞅，因此难以求得解析的首达时间分布。

发明内容

针对现有技术中存在的上述技术问题，本发明提出了一种依赖于时间和状态的长程相关退化过程剩余寿命预测方法，设计合理，克服了现有技术的不足，具有良好的推广效果。

为了实现上述目的，本发明采用如下技术方案：

依赖于时间和状态的长程相关退化过程剩余寿命预测方法，按照如下步骤依次进行：

步骤1：分别在等间隔采样时刻t₀,t₁,t₂,...,t_k，收集设备传感器数据x₀,x₁,x₂,...,x_k，其中，k为采样个数；

步骤2：根据收集的设备传感器数据特征建立基于分形布朗运动的退化模型，如公式(1)所示：

dX(t)＝μ[X(t),t；θ]dt+σ_HdB_H(t)>

其中，X(t)为退化过程，μ[X(t),t；θ]为漂移项系数，θ为漂移项系数中包含的未知参数构成的向量，σ_H为扩散项系数，B_H(t)为标准分形布朗运动；

步骤3：利用基于二次变分的方法，根据公式(2)估计退化模型中的Hurst指数：

其中，γ₁,γ₂,...,γ_p为基于Symlets小波函数的小波分解高通滤波器系数，p为该小波函数消失矩阶数；

步骤4：构造分形布朗运动B_H(t)的基本鞅M^H(t)，如公式(3)所示：

其中，Γ为Gamma函数；

步骤5：根据公式(4)计算分形布朗运动B_H(t)基本鞅M^H(t)的二次变分v^H(t)：

步骤6：构造退化过程的Radon-Nikodym导数R(θ)，如公式(5)所示：

其中，τ为采样间隔；

步骤7：最大化似然比函数ln[R(θ)]，即求得ln[R(θ)]关于漂移项各未知参数的偏导数并令偏导数等于零，解得漂移项各未知参数的估计值，如公式(6)所示：

步骤8：利用极大似然方法，根据公式(7)计算漂移系数估计值：

其中，δ_k＝[δ₁>2 …>k]，Q_k是一个k×k的

矩阵，其第i行第j列元素为其中，τ为采样间隔；

步骤9：将退化模型(1)近似为一个具有时变扩散项系数的基于标准布朗运动的随机过程，如公式(8)所示：

其中，

步骤10：在时刻t_k，重新构造随机过程Y(l_k)：

dY(l_k)＝η(Y(l_k),l_k；θ)dl_k+σ(l_k)dB(l_k)>

其中，l_k＝t-t_k；

步骤11：根据已知失效阈值ω，得到退化过程首达时间的近似解析分布，如公式(12)所示：

其中，ω_k＝ω-x_k。

本发明所带来的有益技术效果：

本发明提出的一种依赖于时间和状态的长程相关退化模型比现有的退化模型更加一般化，更加贴近实际设备的退化；本发明利用Radon-Nikodym导数作为似然比函数估计退化模型未知参数，解决了同时依赖于时间和状态的随机过程的解析似然函数无法得到的问题，估计结果比传统的近似极大似然方法更加准确；为了求解退化过程首达时间，本发明首先利用分形布朗运动的弱收敛理论将退化过程近似为一个具有时变扩散系数的Wiener过程，然后通过一个时间-空间变换，将退化过程首达时间转换为标准布朗运动首达时变阈值的时间，进而得到了退化过程剩余寿命近似的解析分布。

附图说明

图1是本发明退化过程的剩余寿命预测方法的流程图。

图2是本发明示例中高炉炉壁温度退化曲线示意图。

图3是本发明示例中高炉炉壁剩余寿命预测结果示意图。

具体实施方式

下面结合附图以及具体实施方式对本发明作进一步详细说明：

本发明提出的一种依赖于时间和状态的长程相关退化过程的剩余寿命预测方法，流程如图1所示：

步骤1：分别在等间隔采样时刻t₀,t₁,t₂,...,t_k，收集设备传感器数据x₀,x₁,x₂,...,x_k，其中，k为采样个数；

步骤2：根据收集的设备传感器数据特征建立基于分形布朗运动的退化模型，如公式(1)所示：

dX(t)＝μ[X(t),t；θ]dt+σ_HdB_H(t)>

其中，X(t)为退化过程，μ[X(t),t；θ]为漂移项系数，θ为漂移项系数中包含的未知参数构成的向量，σ_H为扩散项系数，B_H(t)为标准分形布朗运动；

步骤3：利用基于二次变分的方法，根据公式(2)估计退化模型中的Hurst指数：

其中，γ₁,γ₂,...,γ_p为基于Symlets小波函数的小波分解高通滤波器系数，p为该小波函数消失矩阶数；

步骤4：构造分形布朗运动B_H(t)的基本鞅M^H(t)，如公式(3)所示：

其中，Γ为Gamma函数；

步骤5：根据公式(4)计算分形布朗运动B_H(t)基本鞅M^H(t)的二次变分v^H(t)：

步骤6：构造退化过程的Radon-Nikodym导数R(θ)，如公式(5)所示：

其中，τ为采样间隔；

步骤7：最大化似然比函数ln[R(θ)]，即求得ln[R(θ)]关于漂移项各未知参数的偏导数并令偏导数等于零，解得漂移项各未知参数的估计值，如公式(6)所示：

步骤8：利用极大似然方法，根据公式(7)计算漂移系数估计值：

其中，δ_k＝[δ₁>2 …>k]，Q_k是一个k×k的

矩阵，其第i行第j列元素为其中，τ为采样间隔；

步骤9：将退化模型(1)近似为一个具有时变扩散项系数的基于标准布朗运动的随机过程，如公式(8)所示：

其中，

步骤10：在时刻t_k，重新构造随机过程Y(l_k)：

dY(l_k)＝η(Y(l_k),l_k；θ)dl_k+σ(l_k)dB(l_k)>

其中，l_k＝t-t_k；

步骤11：根据已知失效阈值ω，得到退化过程首达时间的近似解析分布，如公式(12)所示：

其中，ω_k＝ω-x_k。

为了帮助理解本发明以及展示其剩余寿命预测的效果，下面对一示例进行详细说明。本示例基于MATLAB工具，利用酒泉钢铁2号高炉的数据对本发明进行说明，结合附图展示本发明的效果。

步骤1：收集高炉炉壁内部温度传感器数据，选择位于高8.134米深0.4米的热电偶，连续采样200天的温度数据(采样间隔为8小时，即天，对应的数据长度为600)，记为x₀,x₁,x₂,...,x₅₉₉，退化数据如图2所示；

步骤2：将退化过程建模为如下随机过程：

dX(t)＝[a+btX(t)]dt+σ_HdB_H(t)

其中，a和b为漂移项系数中包含的未知参数，σ_H为扩散项系数，B_H(t)为标准分形布朗运动；

步骤3：在第k个采样时刻进行剩余寿命预测，利用x₀,x₁,x₂,...,x_k，根据式(2)估计退化模型中的Hurst指数，其估计值

其中，γ₁,γ₂,...,γ_p为基于Symlets小波函数的小波分解高通滤波器系数，p为该小波函数消失矩阶数；

步骤4：构造分形布朗运动B_H(t)的基本鞅M^H(t)，如公式(3)所示：

其中，Γ为Gamma函数；

步骤5：根据公式(4)及计算分形布朗运动B_H(t)基本鞅M^H(t)的二次变分v^H(t)：

步骤6：构造退化过程的Radon-Nikodym导数R(θ)，如公式(5)所示：

其中，τ为采样间隔；

步骤7：分别求取R(θ)关于a和b的偏导数，并令偏导数等于0，解得估计值：

步骤8：利用极大似然方法，根据公式(7)计算漂移系数估计值：

其中，δ_k＝[δ₁>2 …>k]，Q_k是一个k×k的

矩阵，其第i行第j列元素为其中，τ为采样间隔；

步骤9：根据式(9)得到离散形式的函数σ(t)，并将退化过程近似为如式(8)所示的随机过程。

其中，

步骤10：在时刻t_k，根据式(11)，构造如下随机过程：

其中，l_k＝t-t_k；

步骤11：假设失效阈值ω为260℃，根据式(12)得到第k个采样时刻剩余寿命分布的预测结果：

其中，η(y,s；θ)＝a+b(t_k+s)(X(t_k)+y)，ω_k＝ω-x_k。

根据高炉炉壁的温度数据可知，首次到达阈值260℃的时间为第588个采样时刻(即第4704小时)。分别在第450、460、470、480、490、500、510、520、530、540、550、560个采样时刻(即第3600、3680、3760、3840、3920、4000、4080、4160、4240、4320、4400、4480小时)进行剩余寿命预测，得到的结果如图3所示。

当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。