CN104915550B - 基于多时间尺度建模的轴承寿命预测方法 - Google Patents

Abstract

本发明公开了基于多时间尺度建模的轴承寿命预测方法，属于轴承寿命预测技术领域，其将能够描述系统快变时间尺度行为的数据驱动模型和能够描述系统慢变时间尺度行为的物理模型有机地结合起来，提出了基于改进相空间弯曲算法和帕里斯裂纹扩展模型的轴承剩余寿命预测算法；利用多维AR模型来构建提升参考模型，改进了传统相空间弯曲算法；同时，为了有效结合相空间弯曲算法和帕里斯模型，提出了基于时间分段策略的改进帕里斯裂纹扩展模型；这种基于多时间尺度建模的预测方法对轴承的剩余使用寿命取得了良好的预测效果，为机械旋转部件的剩余使用寿命预测研究提供了一种新的研究思路。

Description

基于多时间尺度建模的轴承寿命预测方法

技术领域

本发明属于轴承寿命预测技术领域，具体涉及基于多时间尺度建模的轴承寿命预测方法。

背景技术

众所周知，由轴承故障导致的旋转机械不能正常工作的现象屡见不鲜(比如由通用电器和IEEE工业应用学会联合主持的电机可靠性学习中发现由轴承引起的问题占了整个机器故障的40％以上)。因而轴承的故障诊断和寿命预测研究得到了工业界和学术界的普遍关注。传统的轴承故障诊断是针对已出现故障的轴承信号进行识别与分类，而轴承的寿命预测则是在故障发生之前对机械部件进行连续跟踪监测，并基于提取的跟踪特征建立模型对轴承的剩余寿命提前进行预测。因此，轴承的寿命预测研究对于降低停机时间、保证产品质量、提高生产效率有着极为重要的意义。

在现有轴承寿命预测技术当中，基于数据驱动模型的预测技术和基于物理模型的预测技术是两个重要的组成部分。其中，数据驱动模型是直接通过传感器实测信号构建得到，描述的是快变时间尺度下获取的信息，反映的是一个快变时间尺度系统的动态特征。物理模型可以通过动力学分析轴承物理状态的变化以及实际损伤的演化情况，而这种故障状态的演化通常是一个缓慢的过程，可以被描述为一个慢变时间尺度系统，物理模型反映的就是这样一个慢变时间尺度系统的动态特征。实际上，一个现实中的机械动态系统往往可以被看做一个复合系统，由上述的快变时间尺度系统和慢变时间尺度系统相互联系耦合而成。单独使用数据驱动模型或者物理模型往往只能反映整个动态系统的一个方面的特征表现，无法全面地描述快变时间尺度系统和慢变时间尺度系统的动态特征。

发明内容

发明目的：本发明提供了基于多时间尺度建模的轴承寿命预测方法，能够将数据驱动模型以及物理模型有机结合，实现两种预测模型优势互补的多时间尺度模型。

技术方案：本发明的基于多时间尺度建模的轴承寿命预测方法，具体包括以下几个步骤：

1)在初始时间t＝0时刻，通过数据采集设备采集N个轴承振动数据点{x_R(1),x_R(2),...,x_R(N)}，N为数据采集设备的采样频率，利用互信息法(选取延迟时间参数τ，利用虚假邻近点法(该方法由Cao在Practicalmethodfordeterminingtheminimumembeddingdimensionofascalartimeseries一文中提出)选取嵌入维数参数d，采用公式(1)重构N个轴承振动数据点的相空间，并将该相空间作为参考相空间：

$\left\{\begin{array}{c}{y}_R\left(1\right)={\{x_R\left(1\right),x_R(1+\mathrm{\τ}),...,x_R(1+(d-1)\mathrm{\τ})\}}^{\text{'}}\\ .\\ .\\ .\\ y_R\left(n\right)={\{x_R\left(n\right),x_R(n+\mathrm{\τ}),...,x_R(n+(d-1)\mathrm{\τ})\}}^{\text{'}}\\ .\\ .\\ .\\ y_R(N-(d-1)\mathrm{\τ})={\{x_R(N-(d-1)\mathrm{\τ}),x_R(N-(d-2)\mathrm{\τ}),...,x_R\left(N\right)\}}^{\text{'}}\end{array}\right.,n=1,...,N-(d-1)\mathrm{\τ}---\left(1\right);$

2)t＝t+1，通过数据采集设备采集t时刻N个轴承振动数据点{x_t(1),x_t(2),...,x_t(N)}，利用公式(2)，采用与参考相空间相同的延迟时间τ和嵌入维数d重构t时刻N个轴承振动数据点的相空间：

$\left\{\begin{array}{c}{y}_t\left(1\right)={\{x_t\left(1\right),x_t(1+\mathrm{\τ}),...,x_t(1+(d-1)\mathrm{\τ})\}}^{\text{'}}\\ .\\ .\\ .\\ y_t\left(n\right)={\{x_t\left(n\right),x_t(n+\mathrm{\τ}),...,x_t(n+(d-1)\mathrm{\τ})\}}^{\text{'}}\\ .\\ .\\ .\\ y_t(N-(d-1)\mathrm{\τ})={\{x_t(N-(d-1)\mathrm{\τ}),x_t(N-(d-2)\mathrm{\τ}),...,x_t\left(N\right)\}}^{\text{'}}\end{array}\right.,n=1,...,N-(d-1)\mathrm{\τ}---\left(2\right);$

3)对于t时刻相空间中的某一个向量y_t(n)，在参考相空间中寻找p(p>2(τ×d))个与其距离最近的向量y_R(k),y_R(k+1),...,y_R(k+p-1)，其中k＝1,...,N-(d-1)τ，令

$L_k=\left(\begin{array}{c}{y}_R(k+\mathrm{\τ})\\ .\\ .\\ .\\ y_R(k+p+\mathrm{\τ}-1)\end{array}\right),B_k=\left(\begin{array}{c}{y}_R^{\text{'}}(k+\mathrm{\τ}-1),...,y_R^{\text{'}}\left(k\right)\\ .\\ .\\ .\\ y_R^{\text{'}}(k+p+\mathrm{\τ}-2),...,y_R^{\text{'}}(k+p-1)\end{array}\right)---\left(3\right)$

利用公式(4)求出多维自回归模型参数：

${\widehat{\mathrm{\Φ}}}_k={\left(B_k^{\text{'}}{B}_k\right)}^{-1}{B}_k^{\text{'}}{L}_k---\left(4\right)$

利用公式(5)对y_t(n)向量构建跟踪函数：

$e_t\left(n\right)=y_t(n+\mathrm{\τ})-{\left((y_t^{\text{'}}(n+\mathrm{\τ}-1),...,y_t^{\text{'}}\left(n\right)){\widehat{\mathrm{\Φ}}}_k\right)}^{\text{'}}---\left(5\right)$

另外，计算y_t(n)与向量y_R(k),y_R(k+1),...,y_R(k+p-1)的距离，记其中最远的距离为r_n；

4)n＝n+1，若n<N-(d-1)τ，返回步骤3)；否则进入步骤5)；

5)利用t时刻相空间中所有向量对应的所有跟踪函数，使用公式(6)计算t时刻相空间的跟踪指标：

$e_t=\sqrt{\frac{\underset{n=1}{\overset{M}{\mathrm{\&Sigma;}}}q\left(n\right){\left|\right|e_t\left(n\right)\left|\right|}^2}{\underset{n=1}{\overset{M}{\mathrm{\&Sigma;}}}q\left(n\right)}}---\left(6\right)$

其中，q(n)为权值函数，m为t时刻相空间的关联维数；若t<400，返回步骤2)；若t＝400，进入步骤6)；若t>400，进入步骤7)；

6)计算e₁、e₂…e₄₀₀这400个跟踪指标的平均值μ_h及标准差σ_h，并设定μ_h+5σ_h为健康阈值th₁；

7)如果e_t<th₁，返回步骤2)；如果e_t＝th₁，设定t_initial为轴承初始故障出现时间，返回步骤2)；否则，进一步判断e_t是否小于1.5th₁，若是，则此时轴承进入初始裂纹扩展阶段，利用公式(7)计算此时轴承剩余寿命N_t：

$N_t=D_1\ln \frac{3{\mathrm{th}}_1}{e_t}-(t-t_{\mathrm{initial}})(1-\frac{e_t-{\mathrm{th}}_1}{{\mathrm{th}}_1})---\left(7\right)$

其中，D₁为利用轴承初始裂纹扩展阶段的历史经验数据确定的常数；若e_t大于或等于1.5th₁，则此时轴承进入快速裂纹扩展阶段，利用公式(8)计算此时轴承剩余寿命：

$N_t=D_2\ln \frac{3{\mathrm{th}}_1}{e_t}---\left(8\right)$

其中，D₂为利用轴承快速裂纹扩展阶段的历史经验数据确定的常数；

8)若e_t>3th₁，认为此时轴承已损坏，预测结束；否则返回步骤2)。

有益效果：与现有技术相比，本发明具有以下优点：

1、本发明中提出的相空间弯曲算法可以通过实测数据的相空间轨迹变化情况构建非线性跟踪指标，并且将该跟踪指标与实际轴承的物理损伤变量建立线性关系，因此相空间弯曲算法得到的跟踪指标就可以进一步作为帕里斯模型的输入进行轴承剩余寿命预测的研究；这种基于相空间弯曲和帕里斯模型的多时间尺度建模方法，可以有效地将数据驱动模型和物理模型结合起来，对轴承的剩余使用寿命取得了良好的预测效果，为轴承剩余使用寿命预测研究提供了一种新的研究思路；

2、本发明中提出的基于多维AR模型的提升参考模型构建方法使得参考模型从简单的一元模型提升为多元模型，增强了描述相空间弯曲程度的能力，有效改进了传统相空间弯曲算法，提高了其对于轴承退化跟踪的准确性；

3、本发明中提出的基于时间分段的帕里斯模型改进策略根据数据驱动模型得到的跟踪指标的特性改进了传统帕里斯模型的寿命预测公式，使得数据驱动模型和物理模型能够有效结合，显著提高了帕里斯模型预测轴承剩余寿命的准确性。

附图说明

图1为基于多时间尺度建模的轴承寿命预测算法流程图；

图2为改进相空间弯曲算法流程图；

图3为应用多时间尺度建模方法进行轴承寿命预测结果图。

具体实施方式

下面结合附图及具体实例对本发明作进一步说明。

如图1所示，基于多时间尺度建模的轴承寿命预测方法，具体包括以下几个步骤：

1)在初始时间t＝0时刻，通过数据采集设备采集N个轴承振动数据点{x_R(1),x_R(2),...,x_R(N)}，N为数据采集设备的采样频率，利用互信息法选取延迟时间参数τ，利用虚假邻近点法选取嵌入维数参数d，采用公式(1)重构N个轴承振动数据点的相空间，并将该相空间作为参考相空间；

$\left\{\begin{array}{c}{y}_R\left(1\right)={\{x_R\left(1\right),x_R(1+\mathrm{\&tau;}),...,x_R(1+(d-1)\mathrm{\&tau;})\}}^{\text{'}}\\ .\\ .\\ .\\ y_R\left(n\right)={\{x_R\left(n\right),x_R(n+\mathrm{\&tau;}),...,x_R(n+(d-1)\mathrm{\&tau;})\}}^{\text{'}}\\ .\\ .\\ .\\ y_R(N-(d-1)\mathrm{\&tau;})={\{x_R(N-(d-1)\mathrm{\&tau;}),x_R(N-(d-2)\mathrm{\&tau;}),...,x_R\left(N\right)\}}^{\text{'}}\end{array}\right.,n=1,...,N-(d-1)\mathrm{\&tau;}---\left(1\right)$

2)采用改进相空间弯曲算法求跟踪指标的方法参见图2：t＝t+1，通过数据采集设备采集t时刻N个轴承振动数据点{x_t(1),x_t(2),...,x_t(N)}，利用公式(2)，采用与参考相空间相同的延迟时间τ和嵌入维数d重构t时刻N个轴承振动数据点的相空间：

$\left\{\begin{array}{c}{y}_t\left(1\right)={\{x_t\left(1\right),x_t(1+\mathrm{\&tau;}),...,x_t(1+(d-1)\mathrm{\&tau;})\}}^{\text{'}}\\ .\\ .\\ .\\ y_t\left(n\right)={\{x_t\left(n\right),x_t(n+\mathrm{\&tau;}),...,x_t(n+(d-1)\mathrm{\&tau;})\}}^{\text{'}}\\ .\\ .\\ .\\ y_t(N-(d-1)\mathrm{\&tau;})={\{x_t(N-(d-1)\mathrm{\&tau;}),x_t(N-(d-2)\mathrm{\&tau;}),...,x_t\left(N\right)\}}^{\text{'}}\end{array}\right.,n=1,...,N-(d-1)\mathrm{\&tau;}---\left(2\right)$

3)对于t时刻相空间中的某一个向量y_t(n)，在参考相空间中寻找p(p>2(τ×d))个与其距离最近的向量y_R(k),y_R(k+1),...,y_R(k+p-1)，其中k＝1,...,N-(d-1)τ，令

$L_k=\left(\begin{array}{c}{y}_R(k+\mathrm{\&tau;})\\ .\\ .\\ .\\ y_R(k+p+\mathrm{\&tau;}-1)\end{array}\right),B_k=\left(\begin{array}{c}{y}_R^{\text{'}}(k+\mathrm{\&tau;}-1),...,y_R^{\text{'}}\left(k\right)\\ .\\ .\\ .\\ y_R^{\text{'}}(k+p+\mathrm{\&tau;}-2),...,y_R^{\text{'}}(k+p-1)\end{array}\right)---\left(3\right)$

利用公式(4)求出多维自回归模型参数：

${\widehat{\mathrm{\&Phi;}}}_k={\left(B_k^{\text{'}}{B}_k\right)}^{-1}{B}_k^{\text{'}}{L}_k---\left(4\right)$

利用公式(5)对y_t(n)向量构建跟踪函数：

$e_t\left(n\right)=y_t(n+\mathrm{\&tau;})-{\left((y_t^{\text{'}}(n+\mathrm{\&tau;}-1),...,y_t^{\text{'}}\left(n\right)){\widehat{\mathrm{\&Phi;}}}_k\right)}^{\text{'}}---\left(5\right)$

另外，计算y_t(n)与向量y_R(k),y_R(k+1),...,y_R(k+p-1)的距离，记其中最远的距离为r_n；

4)n＝n+1，若n<N-(d-1)τ，返回步骤3)；否则进入步骤5)；

5)利用t时刻相空间中所有向量对应的所有跟踪函数，使用公式(6)计算t时刻相空间的跟踪指标：

$e_t=\sqrt{\frac{\underset{n=1}{\overset{M}{\mathrm{\&Sigma;}}}q\left(n\right){\left|\right|e_t\left(n\right)\left|\right|}^2}{\underset{n=1}{\overset{M}{\mathrm{\&Sigma;}}}q\left(n\right)}}---\left(6\right)$

其中，q(n)为权值函数，m为t时刻相空间的关联维数；若t<400，返回步骤2)；若t＝400，进入步骤6)；若t>400，进入步骤7)；

6)计算e₁、e₂…e₄₀₀这400个跟踪指标的平均值μ_h及标准差σ_h，并设定μ_h+5σ_h为健康阈值th₁；

7)如果e_t<th₁，返回步骤2)；如果e_t＝th₁，设定t_initial为轴承初始故障出现时间，返回步骤2)；否则，开始采用改进帕里斯裂纹扩展模型进行轴承剩余寿命预测，进一步判断e_t是否小于1.5th₁，若是，则此时轴承进入初始裂纹扩展阶段，利用公式(7)计算此时轴承剩余寿命N_t：

$N_t=D_1\ln \frac{3{\mathrm{th}}_1}{e_t}-(t-t_{\mathrm{initial}})(1-\frac{e_t-{\mathrm{th}}_1}{{\mathrm{th}}_1})---\left(7\right)$

其中，D₁为利用轴承初始裂纹扩展阶段的历史经验数据确定的常数；若e_t大于或等于1.5th₁，则此时轴承进入快速裂纹扩展阶段，利用公式(8)计算此时轴承剩余寿命：

$N_t=D_2\ln \frac{3{\mathrm{th}}_1}{e_t}---\left(8\right)$

其中，D₂为利用轴承快速裂纹扩展阶段的历史经验数据确定的常数；

8)若e_t>3th₁，认为此时轴承已损坏，预测结束；否则返回步骤2)。

参见图3，在ZA2115轴承的疲劳实验中，从时刻586开始对轴承进行寿命预测，图中横轴代表每一次寿命预测的起始点，纵轴代表对应预测起始点的剩余使用寿命，单位均为10分钟/点。真实剩余寿命是轴承在预测起始点的实际剩余寿命，预测剩余寿命为利用多时间尺度建模方法预测得到的轴承剩余寿命。其中，公式(7)中D₁通过ZA2115轴承的历史疲劳实验得到的初始故障出现时间对应的轴承剩余寿命N_th求得：

D₁＝N_th/ln3

相似的，公式(8)中D₂通过ZA2115轴承的历史疲劳实验得到的进入快速裂纹扩展阶段时对应的轴承剩余寿命N_th2求得：

D₂＝N_th2/ln2

通过对轴承在382个时刻点进行寿命预测，得到平均预测误差为4.61小时。

应理解上述实施例仅用于说明本发明技术方案的具体实施方式，而不用于限制本发明的范围。在阅读了本发明之后，本领域技术人员对本发明的各种等同形式的修改和替换均落于本申请权利要求所限定的保护范围。

Claims (1)

1.基于多时间尺度建模的轴承寿命预测方法，其特征在于，该方法具体包括以下几个步骤：

1)在初始时间t＝0时刻，通过数据采集设备采集N个轴承振动数据点{x_R(1),x_R(2),...,x_R(N)}，N为数据采集设备的采样频率，利用互信息法选取延迟时间参数τ，利用虚假邻近点法选取嵌入维数参数d，采用公式(1)重构N个轴承振动数据点的相空间，并将该相空间作为参考相空间：

$\left\{\begin{array}{c}{y}_R\left(1\right)={\{x_R\left(1\right),x_R(1+\mathrm{\&tau;}),...,x_R(1+(d-1\left)\mathrm{\&tau;}\right)\}}^{\&prime;}\\ \begin{array}{c}.\\ .\\ .\end{array}\\ y_R\left(n\right)={\{x_R\left(n\right),x_R(n+\mathrm{\&tau;}),...,x_R(n+(d-1\left)\mathrm{\&tau;}\right)\}}^{\&prime;}\\ \begin{array}{c}.\\ .\\ .\end{array}\\ y_R(N-(d-1\left)\mathrm{\&tau;}\right)={\{x_R(N-(d-1\left)\mathrm{\&tau;}\right),x_R(N-(d-2\left)\mathrm{\&tau;}\right),...,x_R(N\left)\right\}}^{\&prime;}\end{array}\right.,n=1,...,N-(d-1)\mathrm{\&tau;}---\left(1\right)$

进入步骤2)；

2)t＝t+1，通过数据采集设备采集t时刻N个轴承振动数据点{x_t(1),x_t(2),...,x_t(N)}，利用公式(2)，采用与参考相空间相同的延迟时间τ和嵌入维数d重构t时刻N个轴承振动数据点的相空间：

$\left\{\begin{array}{c}{y}_t\left(1\right)={\{x_t\left(1\right),x_t(1+\mathrm{\&tau;}),...,x_t(1+(d-1\left)\mathrm{\&tau;}\right)\}}^{\&prime;}\\ \begin{array}{c}.\\ .\\ .\end{array}\\ y_t\left(n\right)={\{x_t\left(n\right),x_t(n+\mathrm{\&tau;}),...,x_t(n+(d-1\left)\mathrm{\&tau;}\right)\}}^{\&prime;}\\ \begin{array}{c}.\\ .\\ .\end{array}\\ y_t(N-(d-1\left)\mathrm{\&tau;}\right)={\{x_t(N-(d-1\left)\mathrm{\&tau;}\right),x_t(N-(d-2\left)\mathrm{\&tau;}\right),...,x_t(N\left)\right\}}^{\&prime;}\end{array}\right.,n=1,...,N-(d-1)\mathrm{\&tau;}---\left(2\right)$

进入步骤3)；

3)对于t时刻相空间中的某一个向量y_t(n)，在参考相空间中寻找p个与其距离最近的向量y_R(k),y_R(k+1),...,y_R(k+p-1)，其中p>2(τ×d),k＝1,...,N-(d-1)τ，令

$L_k=\left(\begin{array}{c}{y}_R(k+\mathrm{\&tau;})\\ \begin{array}{c}.\\ .\\ .\end{array}\\ y_R(k+p+\mathrm{\&tau;}-1)\end{array}\right),B_k=\left(\begin{array}{c}{y}_R^{\&prime;}(k+\mathrm{\&tau;}-1),...,y_R^{\&prime;}\left(k\right)\\ \begin{array}{c}.\\ .\\ .\end{array}\\ y_R^{\&prime;}(k+p+\mathrm{\&tau;}-2),...,y_R^{\&prime;}(k+p-1)\end{array}\right)---\left(3\right)$

利用公式(4)求出多维自回归模型参数：

${\widehat{\mathrm{\&Phi;}}}_k={\left(B_k^{\&prime;}{B}_k\right)}^{-1}{B}_k^{\&prime;}{L}_k---\left(4\right)$

利用公式(5)对y_t(n)向量构建跟踪函数：

$e_t\left(n\right)=y_t(n+\mathrm{\&tau;})-{\left(\right(y_t^{\&prime;}\left(n+\mathrm{\&tau;}-1\right),...,y_t^{\&prime;}\left(n\right)\left){\widehat{\mathrm{\&Phi;}}}_k\right)}^{\&prime;}---\left(5\right)$

另外，计算y_t(n)与向量y_R(k),y_R(k+1),...,y_R(k+p-1)的距离，记其中最远的距离为r_n，进入步骤4)；

4)n＝n+1，若n<N-(d-1)τ，返回步骤3)；否则进入步骤5)；

5)利用t时刻相空间中所有向量对应的所有跟踪函数，使用公式(6)计算t时刻相空间的跟踪指标：

$e_t=\sqrt{\frac{\underset{n=1}{\overset{M}{\&Sigma;}}q\left(n\right)\left|\right|e_t\left(n\right)|{|}^2}{\underset{n=1}{\overset{M}{\&Sigma;}}q\left(n\right)}}---\left(6\right)$

其中，q(n)为权值函数，m为t时刻相空间的关联维数；若t<400，返回步骤2)；若t＝400，进入步骤6)；若t>400，进入步骤7)；

6)计算e₁、e₂…e₄₀₀这400个跟踪指标的平均值μ_h及标准差σ_h，并设定μ_h+5σ_h为健康阈值th₁，进入步骤7)；

7)步骤6)中所述的th_1，如果e_t<th₁，返回步骤2)；如果e_t＝th₁，设定t_initial为轴承初始故障出现时间，返回步骤2)；否则，进一步判断e_t是否小于1.5th₁，若是，则此时轴承进入初始裂纹扩展阶段，利用公式(7)计算此时轴承剩余寿命N_t：

$N_t=D_1\ln \frac{3{\mathrm{th}}_1}{e_t}-(t-t_{initial})(1-\frac{e_t-{\mathrm{th}}_1}{{\mathrm{th}}_1})---\left(7\right)$

其中，D₁为利用轴承初始裂纹扩展阶段的历史经验数据确定的常数；若e_t大于或等于1.5th₁，则此时轴承进入快速裂纹扩展阶段，利用公式(8)计算此时轴承剩余寿命：

$N_t=D_{2}ln\frac{3{\mathrm{th}}_1}{e_t}---\left(8\right)$

其中，D₂为利用轴承快速裂纹扩展阶段的历史经验数据确定的常数，进入步骤8)；

8)若e_t>3th₁，认为此时轴承已损坏，预测结束；否则返回步骤2)。