融合动力学模型的滚动轴承性能退化评估方法

Rolling bearings performance degradation assessment method for integrating dynamic model

摘要

针对在线场景下单一健康指标难以准确反映滚动轴承健康状态的问题， 提出一种融合动力学模型的滚动轴承性能退化评估方法. 在现有轴承动力学模型基础上， 考虑转速变化引起的滑动摩擦力对振动响应的影响， 构建可表征不同转速下滚动轴承健康状态的动力学模型，设计基于卷积神经网络的模型更新方法， 解决动力学模型的在线同步更新问题. 融合在线同步动力学模型， 求解动力学模型接触力特征， 以构建表征滚动轴承健康状态的敏感指标，进而建立基于实测信号与仿真信号相似度动态修正阈值的自适应健康状态评估策略， 实现滚动轴承退化程度的自主划分和健康指标阈值动态更新. 利用西安交通大学、美国国家航空航天局及实验室自建平台滚动轴承振动实验数据集完成方法有效性验证. 实验结果表明：相较于多元多尺度模糊熵、特征相似比等主流方法， 所提出的方法准确率提升了1.67%， 占用内存大小降低了30.4%， 浮点计算量降低了4.55%.

Abstract

Aiming to address the issue that single health indicators are insufficient to accurately reflect the real-time health status of rolling bearings in online scenarios， a performance degradation assessment method for rolling bearings based on a dynamic model is proposed. Built on the existing bearing dynamic model， the influence of sliding friction forces induced by speed variations on vibration response is considered， and a dynamic model capable of characterizing the health state of rolling bearings at different speeds is constructed. A model updating method based on convolutional neural networks is designed to solve the problem of online synchronous updating of dynamic models. By integrating the online synchronous dynamic model， the contact force characteristics of the dynamic model are analyzed to construct a sensitive index for characterizing the health status of the rolling bearing. Furthermore， an adaptive health status assessment strategy is established based on the dynamic correction threshold of the similarity between the measured and simulated signal. This strategy enables autonomous division of the degradation degree of the rolling bearing and dynamic update of the health index threshold. Finally， the validity of the method is verified using rolling bearing vibration experimental datasets from Xi’an Jiaotong University （XJTU-SY）， the National Aeronautics and Space Administration （NASA）， and a laboratory self-built platform. The experimental results demonstrate that， compared with the mainstream methods such as multivariate multi-scale fuzzy entropy and feature similarity ratio， the proposed method improves accuracy by 1.67%， reduces memory usage by 30.4%， and decreases floating-point computations by 4.55%.

Graphical abstract

关键词

滚动轴承 / 卷积神经网络 / 故障诊断 / 动力学模型 / 退化评估

Key words

rolling bearing / convolutional neural network / fault diagnosis / dynamic model / degradation assessment

滚动轴承作为旋转机械的重要部件之一，在设备运转中起着支撑和导向作用，其健康状态影响着设备的运转性能.在轴承设备服役期间，突发性故障或故障恶化会影响设备的整体运行状态，造成巨大的安全隐患和经济损失.因此，滚动轴承健康状态的在线评估有利于保障设备安全运行，对于实现设备的预测性维护具有十分重要的意义^［1-3］.

为实现滚动轴承的健康状态评估， 国内外学者主要围绕动力学建模和健康指标构建来开展研究. 在滚动轴承动力学模型的研究中，WU等^［4］构建了五自由度轴承的动力学模型，研究了振动信号中双冲击特性与轴承健康状态的关系.田晶等^［5］考虑复合故障下的耦合特性，构建了可表征复合故障的中介轴承动力学模型. JIANG等^［6］考虑了故障区域与滚动体之间的三维几何关系对振动响应的影响， 提高了动力学模型的建模精度. 罗茂林等^［7］考虑冲击力对振动响应的影响， 构建了基于双冲击时间间隔表征轴承健康状态的动力学模型. 胡爱军等^［8］考虑时变刚度耦合对振动响应的影响， 研究了不同健康状态下轴承振动响应的多冲击变化特征. 上述动力学模型通过考虑不同的因素来表征滚动轴承的健康状态， 但未考虑变转速复杂工况对动力学建模的影响， 忽略了动力学模型随实时信号同步更新的过程， 模型的准确性和计算速度有待提升.

在动力学模型研究基础上， 国内外学者通过构建健康指标来量化轴承的健康状态. 目前的研究一般直接使用信号时频域特征（峰值、均方根、方差等）作为健康指标来描述轴承的健康状态. 但使用时频域特性直接描述轴承健康状态时， 容易出现虚假波动， 造成健康状态评估不准确. 为进一步提升健康状态评估的准确性， YAO等^［9］提出一种融合熵特征的关键频谱指标， 提升了轴承健康状态评估的准确性. PAN等^［10］提出基于深度置信网络和改进支持向量数据描述的方法来表征轴承的健康状态， 实现了轴承的性能退化评估. 张龙等^［11］通过分析不同时频域指标在相邻样本间的最佳特征相似比来构建健康指标， 并实现了轴承退化阶段的划分. WANG等^［12］通过图谱重构和谱表征的方法来构建轴承健康状态敏感指标， 用于实现轴承健康状态的定量评估. 赵聪聪等^［13］提出了一种基于多元多尺度模糊熵和可拓聚类的方法来表征轴承的健康状态， 进一步提升了健康状态评估的准确性. 上述健康指标的构建方法为轴承的健康状态评估提供了新的思路， 但受时变工况和实时数据流的影响， 在线运行对于健康指标的准确性和计算速度要求较高， 健康指标计算的准确性和计算速度还有待提升.

为提升滚动轴承性能退化评估方法的准确性和计算速度， 本文提出一种融合动力学模型的滚动轴承性能退化评估方法， 主要工作如下：考虑转速变化引起的滑动摩擦力对振动响应的影响， 构建可表征不同转速下滚动轴承健康状态的动力学模型. 设计基于卷积神经网络的模型更新方法， 解决动力学模型在线同步更新问题. 融合在线同步动力学模型， 求解动力学模型接触力特征， 以构建表征滚动轴承健康状态的敏感指标， 进而建立基于实测信号与仿真信号相似度动态修正阈值的自适应健康状态评估策略， 实现滚动轴承退化程度的自主划分和健康指标阈值动态更新. 最后， 利用3个不同实验平台的轴承数据集完成所提方法有效性验证.

1 滚动轴承动力学模型

滚动轴承动力学模型的构建流程包括3部分：首先，构建可用于变转速复杂工况下的动力学模型；其次，基于卷积神经网络实现模型的同步更新；最后，通过与理论值对比来验证模型的准确性.

1.1 动力学模型构建

考虑轴承的内、外圈两个部件在水平和竖直方向上的运动， 建立二维坐标系， 如图1所示.

在此坐标系下，构建非线性动力学方程，如 式（1）所示.

m_{\mathrm{i}}{\ddot{x}}_{\mathrm{i}}+c_{\mathrm{i}}{\dot{x}}_{\mathrm{i}}+k_{\mathrm{i}}{x}_{\mathrm{i}}+f_x=\mathrm{0}

m_{\mathrm{i}}{\ddot{y}}_{\mathrm{i}}+c_{\mathrm{i}}{\dot{y}}_{\mathrm{i}}+k_{\mathrm{i}}{y}_{\mathrm{i}}+f_y+m_{\mathrm{i}}g=\mathrm{0}

m_{\mathrm{o}}{\ddot{x}}_{\mathrm{o}}+c_{\mathrm{o}}{\dot{x}}_{\mathrm{o}}+k_{\mathrm{o}}{x}_{\mathrm{o}}-f_x=\mathrm{0}

m_{\mathrm{o}}{\ddot{y}}_{\mathrm{o}}+c_{\mathrm{o}}{\dot{y}}_{\mathrm{o}}+k_{\mathrm{o}}{y}_{\mathrm{o}}-f_y+m_{\mathrm{o}}g=-F_{\mathrm{r}}（1）

式中： x_i为内圈x方向位移； x_o为外圈x方向位移； y_o为外圈y方向位移； y_i为内圈y方向位移； F_r为外加径向载荷；f_x 为所有滚动体与滚道接触力在水平方向（x方向）的合力分量；f_y 为所有滚动体与滚道接触力在竖直方向（y方向）的合力分量； g为重力加速度. 方程及模型中的其他参数如表1所示.

假设轴承以恒定转速转动， 外圈固定则内圈与滚动体存在公转速度关系， 采用式（2）表示^［4］.

w_{\mathrm{b}}=\frac{w_{\mathrm{i}}}{\mathrm{2}}\left(\mathrm{1}-\frac{D_{\mathrm{b}}\mathrm{c}\mathrm{o}\mathrm{s}\alpha }{D_{\mathrm{p}}}\right)（2）

式中： w_b为滚动体公转角速度； w_i为内圈旋转角速度； D_b为滚动体直径； D_p为节圆直径； α为滚动接触角.

在实际生产中， 轴承的旋转速度较高， 可能存在打滑现象， 式（2）不能反映打滑现象， 不适用于高转速情况， 因此需要对现有的公转速度公式进行修正. 内圈与滚动体间存在动能传递， 内圈保持恒定转速旋转， 通过滚动力作用于滚动体， 使滚动体发生旋转. 当内圈旋转产生的滚动力大于摩擦临界值时， 滚动体与滚道发生相对滑动. 将多个滚动体看成一个整体， 则恒定转速下内圈产生滚动力的计算如式（3）所示. 当滚动体即将发生打滑时， 滚道表面的摩擦由滚动摩擦转变为滑动摩擦， 滑动摩擦力的计算如式（4）所示.

f_{\mathrm{i}}=m_{\mathrm{i}}{w}_{\mathrm{i}}^{\mathrm{2}}（3）

f_{\mathrm{c}}=\mu m_{\mathrm{i}}g\mathrm{c}\mathrm{o}\mathrm{s}\theta +\mu F_{\mathrm{r}}\mathrm{c}\mathrm{o}\mathrm{s}\theta（4）

式中： f_i为内圈旋转产生的滚动力； f_c为滑动摩擦力； μ为滑动摩擦因子； θ为滚动体等效旋转位置角.

基于动能定理可知， 发生相对滑动时动能的变化量为滑动摩擦力做功， 可得如式（5）所示的等量关系.

\frac{\mathrm{1}}{\mathrm{2}}{m}_{\mathrm{i}}{w}_{\mathrm{i}}^{\mathrm{2}}-\frac{\mathrm{1}}{\mathrm{2}}{m}_{\mathrm{b}}{w}_{\mathrm{f}}^{\mathrm{2}}={\int }_{{\theta }_{\mathrm{1}}}^{{\theta }_{\mathrm{2}}}{f}_{\mathrm{c}}\mathrm{d}\theta（5）

式中： {\int }_{{\theta }_{\mathrm{1}}}^{{\theta }_{\mathrm{2}}}{f}_{\mathrm{c}}\mathrm{d}\theta表示单圈内滑动摩擦力的做功， θ₁为0， θ₂为2π. 假设轴承的外部负载分布不变， 则每圈的滑动摩擦相同， 且每圈滑动摩擦所做的功相同， 可推导出相对滑动速度w_f的计算式，如式（6）所示.

w_{\mathrm{f}}=\left\{\begin{array}{l}\sqrt[]{\left(\frac{\mathrm{1}}{\mathrm{2}}{m}_{\mathrm{i}}{w}_{\mathrm{i}}^{\mathrm{2}}-{\int }_{{\theta }_{\mathrm{1}}}^{{\theta }_{\mathrm{2}}}{f}_{\mathrm{c}}\mathrm{d}\theta \right)\frac{\mathrm{2}}{m_{\mathrm{b}}}},f_c\le f_i\\ \mathrm{0},f_c>f_i\end{array}\right.（6）

本文以是否发生相对滑动来界定高速或低速.判定的临界条件为单圈内滑动摩擦力的最小值等于当前转速下内圈提供的滚动力，可推导出临界转速w_c公式，如式（7）所示.若负载不变，转速大于w_c时， 则视为高速运行，转速小于w_c时，则视为低速运行. 滚动体的公转角速度修正如式（8）所示.

w_{\mathrm{c}}=\sqrt[]{\mathrm{m}\mathrm{i}\mathrm{n}(\mu m_{\mathrm{i}}g\mathrm{c}\mathrm{o}\mathrm{s}\theta +\mu F_{\mathrm{r}}\mathrm{c}\mathrm{o}\mathrm{s}\theta )/m_{\mathrm{i}}}（7）

w_{\mathrm{b}}=\frac{w_{\mathrm{i}}}{\mathrm{2}}\left(\mathrm{1}-\frac{D_{\mathrm{b}}\mathrm{c}\mathrm{o}\mathrm{s}\alpha }{D_{\mathrm{p}}}\right)+w_{\mathrm{f}}（8）

基于此， 滚动体的角度位置关系如式（9）所示， t表示当前时刻， w_b表示滚动体公转角速度， θ₀表示自定义的1号滚动体初始位置角， j表示滚动体编号， N_b表示滚动体个数.

{\theta }_j={\theta }_{\mathrm{0}}+w_{\mathrm{b}}t+\frac{\mathrm{2}\mathrm{\pi }(j-\mathrm{1})}{N_{\mathrm{b}}}（9）

基于赫兹接触理论， 可得到滚动体在不同位置的接触变形量［式（10）］和接触力［式（11）］所示^［4］.

{\delta }_j=(x_{\mathrm{o}}-x_i)\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_j+(y_{\mathrm{o}}-y_i)\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_j-C_{\mathrm{r}}（10）

式中：{\delta }_j为第j个滚动体的接触变形量； {\theta }_j为第j个滚动体的角度位置； C_r为径向游隙.

\left\{\begin{array}{c}{f}_x=k_{\mathrm{b}}\sum {\delta }_j^n\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_j\\ f_y=k_{\mathrm{b}}\sum {\delta }_j^n\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_j\end{array}\right.（11）

式中： k_{\mathrm{b}}为滚动体（赫兹接触）等效刚度；n为接触力系数；

当存在外圈故障时， 滚动体经过故障区域时， 会直接影响接触变形量， 在非故障模型的基础上多一个新的变形量H_d， 如式（12）所示.

{\delta }_j=(x_{\mathrm{o}}-x_{\mathrm{i}})\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_j+(y_{\mathrm{o}}-y_{\mathrm{i}})\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_j-C_{\mathrm{r}}-H_{\mathrm{d}}（12）

假设故障为剥落矩形坑， 矩形宽度为L， 矩形深度为H， 则故障对应的圆心角如式（13）所示.

{\varphi }_{\mathrm{s}}=\mathrm{2}\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{i}\mathrm{n}(L/D_{\mathrm{o}})（13）

式中：D_o为外圈滚道直径.

滚动体进出故障区域的过程中所引起的接触变形量不一致， 分为两种情况， 以滚动体是否能触碰到故障底部为判断条件， 使用分段函数来表征.

当故障较小时， 即满足不等式：

变形量的表达式如式（14）所示：

（14）

当故障较大时， 即不满足不等式：

变形量的表达式如式（15）所示：

（15）

滚动体在进入故障区域到完全离开分为多个阶段， φ_j 为滚动体位置角相对于故障起始点的角度， θ_d为滚动体由刚进入切换到已经进入状态的滚动体圆心角的差， 根据几何关系， 两个角度的公式分别如式（16）和式（17）所示.

（16）

{\theta }_{\mathrm{d}}=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{D_{\mathrm{o}}^{\mathrm{2}}/\mathrm{4}+(D_{\mathrm{o}}/\mathrm{2}-D_{\mathrm{b}}{/\mathrm{2}+H)}^{\mathrm{2}}-D_{\mathrm{b}}^{\mathrm{2}}/\mathrm{4}}{D_{\mathrm{o}}(D_{\mathrm{o}}/\mathrm{2}-D_{\mathrm{b}}/\mathrm{2}+H)}\right)（17）

以上是轴承外圈故障的动力学建模， 内圈故障建模与外圈故障类似， 不同的是内圈存在一个角速度， 故障对应的圆心角是时变的. 因此， 用式（18）代替式（13）.

{\varphi }_{\mathrm{s}}=\mathrm{2}\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{i}\mathrm{n}(L/D_{\mathrm{i}})+w_{\mathrm{i}}t（18）

式中： D_i为内圈滚道直径

1.2 模型验证

借鉴文献［14］中的模型验证思路， 引入式（19）的特征频率指标， 对比文献［15］中的模型， 设置的仿真参数如表2所示， 得到的结果如图2所示.

f_{\mathrm{o}}=\frac{w_{\mathrm{i}}}{\mathrm{2}}\left(\mathrm{1}-\frac{D_{\mathrm{b}}}{D_{\mathrm{p}}}\mathrm{c}\mathrm{o}\mathrm{s}\alpha \right)N_{\mathrm{b}}（19）

经计算， 当转速为3 000 r/min时， 理论故障特征频率为154.15 Hz， 对比模型特征频率与理论频率的相对误差为1.849%， 远大于本文模型特征频率与理论频率的相对误差0.097%. 这主要得益于本文模型考虑了转速变化引起的滑动摩擦力对振动响应的影响， 修正了滚动体的公转速度公式， 使得内圈、滚动体、外圈的相对运动关系更接近于实际， 更适用于表征不同转速下的轴承健康状态.

1.3 模型同步更新

为提高在线应用场景下轴承健康状态评估的准确性和计算速度， 需建立可根据实时信号同步更新的模型，其本质上是构建一个振动信号与动力学模型参数之间的实时映射关系， 应基于以下原则进行构建：精准映射、实时计算. 即可通过动力学模型仿真实验分析振动信号的变化规律， 进而辅助映射关系的构建. 动力学模型中的可变参数为故障深度H、故障宽度L， 设置不同的参数如表3所示， 得到如 图3所示的仿真信号.

从图3中可看出，相较于正常情况，故障时仿真信号中会出现更多的周期性冲击信号.为进一步描述不同健康状态下的实时信号与模型参数之间的联系，通过均方根、方差、最大值、波形因子、峰值因子、脉冲因子、裕度因子、均方频率、频率方差等时频域指标来分析信号特征.固定故障深度分别为5.5、2.5、0.5 mm，改变故障宽度（从0 mm扩展到10 mm）进行仿真，计算上述时频域指标，得到如图4所示的变化曲线.固定故障宽度分别为0.5、4.5、8.5 mm，改变故障深度（从0 mm扩展到6 mm）进行仿真，计算上述时频域指标， 得到如图5所示的变化曲线.

通过图4、图5可看出，不同时频域指标表现不一致.图4中， 在故障深度不变时，随故障宽度的增加，均方根、方差、最大值指标存在一个递增的趋势；波形因子、峰值因子、脉冲因子、裕度因子等指标则先递增再递减；而均方频率、频率方差指标则先递增，后在一个范围内波动.图5中，在故障宽度不变的情况下，随故障深度的增加，所有指标均先递增再保持不变.根据曲线的变化规律，可得出定性结论：当故障尺寸较小时，故障深度和宽度都会影响指标趋势，且故障尺寸增加会造成指标递增；当故障尺寸增加到一定程度时，故障深度变化不再影响指标趋势，故障宽度增加会造成方差、均方根、最大值三个指标的递增.相较于其他指标，方差、均方根、最大值三个指标与故障深度、故障宽度的联系更紧密.因此，选择方差、均方根、最大值三个指标作为映射关系构建的基础特征.

为实现映射关系的实时计算， 应选取计算速度快的方法来构建映射关系. 因卷积神经网络（convolutional neural network， CNN）具有可提取局部特征和运行速度快的优势， 适用于实时信号的特征提取与跟踪^［16］. 因此， 采用CNN来实现映射关系的构建. 以实时振动信号的均方根、方差、最大值三个特征作为CNN的输入， 以故障深度H和故障宽度L作为CNN的输出. 其训练过程如下：

Step1：数据准备. 设置不同的故障深度H和故障宽度L， 采用1.1节中的动力学模型进行仿真， 得到对应的仿真信号， 计算仿真信号的均方根、方差、最大值三个特征， 构建一个以振动信号的均方根、方差、最大值作为输入， 以故障深度H和故障宽度L作为输出的训练集. 将训练集按照7∶3∶1的比例划分为训练集、验证集和测试集， 同时对数据进行归一化处理， 以提升模型的计算效率.

Step2：构建模型.按照确定的网络层次结构（包括卷积层、池化层和全连接层）及相关参数进行网络模型构建.其中，全连接层采用ReLU激活函数， 具体网络结构及参数设置如表4所示.

Step3：训练设置.选择均方误差作为损失函数， 优化器为Adam，并设置初始学习率为0.001， batch_size=32，同时使用小的随机值初始化网络权重.使用训练集数据对模型进行训练，在每个epoch中，数据会被送入网络进行前向传播，计算损失函数，然后进行反向传播和权重更新.设置每20个epoch保存一次模型的权重，以便在训练中断时可以从最近的检查点恢复.

Step4：验证调整.在每个epoch后， 使用验证集评估模型的性能.根据验证结果调整网络超参数， 如增加正则化、改变学习率等.模型训练过程中的损失函数曲线如图6所示.

Step5：计算模型每次迭代的准确率和损失函数，当模型满足收敛条件时，完成模型训练并保存模型参数，最终得到训练好的CNN模型.通过图6可发现，随着训练次数的增加，损失函数曲线逐渐下降并收敛，均方误差下降到0.01左右，停止训练模型，保存模型参数.

CNN训练完成后， 使用CNN来实现动力学模型的同步更新， 最后完成接触力的计算， 如图7所示. 在动力学模型使用阶段， 完整的数据流是输入实测信号， 计算均方根、方差、最大值三个特征， 再将特征输入CNN得出故障深度H和故障宽度L， 将故障深度H和故障宽度L输入1.1节所述动力学模型， 最后输出接触力.

2 融合动力学模型的指标构建

为了更准确地反映轴承的健康状态， 本文的健康指标构建包括两部分：

1）计算能表征滚动轴承健康状态的接触力；

2）基于接触力的统计不变特性， 提取接触力的统计特征， 完成健康指标的构建.

2.1 接触力特性分析

通过前文的动力学模型， 可计算得到实时接触力. 从文献［8］可知， 当轴承健康状态发生变化时， 接触力会随之发生改变. 接触力发生改变后， 再与轴承动力学模型中的其他变量发生耦合， 最终造成振动信号发生改变. 因此， 从轴承的振动机理上分析， 接触力与轴承的健康状态相关， 而振动信号包含了与健康状态无关的冗余信息， 相较于直接使用振动信号， 接触力更适用于描述轴承的健康状态.

基于轴承退化数据进行计算， 得到如图8所示的结果. 其中， 图8（a）为轴承全生命周期退化过程的振动信号， 图8（b）为动力学模型计算得到的接触力. 从图8中可知， 振动信号在退化过程中存在虚假波动， 整体变化趋势未随轴承退化程度的增加而单调变化， 而接触力的整体趋势则随着轴承退化程度的增加而递增. 从结果上看， 相较于直接使用振动信号， 接触力更适用于描述轴承的健康状态.

经动力学模型计算得出的接触力的瞬态是一个与所输入振动信号等长的序列， 递增趋势不明显， 借鉴文献［1］中提取信号统计特征的思路， 通过对接触力的瞬态进行统计特征提取，计算接触力的方差、均方根等统计特征， 从而放大递增趋势， 以便更精准地量化滚动轴承的退化趋势.

2.2 接触力统计特征提取与指标构建

通过单调性（Mon）、相关性（Corr）、鲁棒性（Rob）三个方面对接触力的统计特征（方差、均方根、最大值）进行评估和筛选^［1］， 计算公式如式（20）所示， 计算结果如图9所示. 单调性表征的是退化趋势， 鲁棒性表征抗噪声能力， 相关性表征的是接触力与统计特征之间的联系.

\begin{array}{l}\mathrm{M}\mathrm{o}\mathrm{n}(f)=\frac{\mathrm{1}}{N-\mathrm{1}}\left|\sum _k\delta (f_T(t_{k+\mathrm{1}}))-\sum _k\delta (f_T(t_k)-f_T(t_{k+\mathrm{1}}))\right|\\ \mathrm{C}\mathrm{o}\mathrm{r}\mathrm{r}(f)=\frac{\left|N\sum _k{f}_{\mathrm{T}}(t_k)t_k-N\sum _k{f}_{\mathrm{T}}(t_k)\sum _k{t}_k\right|}{\sqrt[]{\left[N\sum _k{f}^{\mathrm{2}}(t_k)-{\left(\sum _k{f}_{\mathrm{T}}(t_k)\right)}^{\mathrm{2}}\right]\left[N\sum _k{t_k}^{\mathrm{2}}-{\left(\sum _k{t}_k\right)}^{\mathrm{2}}\right]}}\\ \mathrm{R}\mathrm{o}\mathrm{b}(f)=\frac{\mathrm{1}}{N}\sum _{k}exp\left(-\left|\frac{f_{\mathrm{R}}(t_k)}{f_{\mathrm{T}}(t_k)}\right|\right)\end{array}（20）

式中： f为待评价的统计特征序列； k为样本序号； N为全生命周期数据的样本总数； t_k 为第k个样本对应的时刻； f（t_k ）为特征f在第k个样本处的取值； f_T（t_k ）为特征f在t_k 处的趋势值； δ（·）为单位阶跃函数.

通过图9可知， 方差在单调性、相关性、鲁棒性上评分最高. 因此， 通过计算接触力的方差来构建健康指标， 如式（21）所示.

H_{\mathrm{I}}=\mathrm{V}\mathrm{a}\mathrm{r}(F_n)（21）

式中： H_I为健康指标； Var（·）为方差运算符； F_n 为由动力学模型计算得到的接触力序列， 由式（11）求得.

3 滚动轴承性能退化评估

轴承性能退化评估框架如图10所示. 实时健康监测是基于动态变化的实时数据流进行故障识别和异常检测， 而性能退化评估是通过健康指标表征滚动轴承完整退化趋势， 完成健康阶段划分^［9］.

3.1 实时健康监测方法

传统的健康监测方法是根据专家经验设置固定阈值， 当健康指标超过固定阈值时触发报警. 这种方法存在的问题是：当信号波动较大时， 即使设备处于正常情况， 其健康指标H_I也可能超出固定阈值引起误报， 或信号波动范围太小， 即使设备处于故障情况H_I也难以达到固定阈值， 造成漏报警. 因此， 在连续的健康监测过程中， 健康指标的阈值需要根据实时信号进行实时更新. 本文的实时健康监测方法如下：

1）确定基准值：确定当前设备的工况， 借鉴文献［17］中对轴承动力学模型中故障尺寸的界定方法， 设置动力学模型故障尺寸为0， 求解动力学模型， 计算H_I值， 记作H_I0.

2）设置区间范围：以H_I0作为基准阈值， 参考文献［18］的思路， 引入拉依达准则设置动态区间， 则故障阈值可表示为[H_{\mathrm{I}\mathrm{0}}-\delta ,H_{\mathrm{I}\mathrm{0}}+\delta ] . 当H_I值大于区间上限时， 认为轴承处于故障状态， 当H_I值小于区间下限时， 轴承处于正常状态. 当H_I值处于区间内时， 将此样本标记为不确定性样本， 通过后续输入的数据继续进行判断.

3）计算H_I值， 更新区间：将实时振动信号输入动力学模型， 计算出仿真信号， 同时计算出H_I值. 参照文献［19］中的信号相似度计算方法， 计算实时信号与仿真信号的余弦相似度， 记作γ. 通过γ来修正故障阈值， 当相似度越高时， 说明此次计算的H_I值可靠性更高， 不确定性的阈值区间范围可以缩小， 反之， 不确定性的阈值区间范围可以扩大. 相似度的计算结果在0~1范围内， 假设当相似度为0.5时， 保持阈值区间不变， 其计算式可表示为[H_{\mathrm{I}\mathrm{0}}-\delta -\mathrm{0.5}+\gamma ,H_{\mathrm{I}\mathrm{0}}+\delta +\mathrm{0.5}-\gamma ].

4）故障识别：当实时数据的H_I大于故障阈值区间的上限时， 触发故障报警， 当H_I小于故障阈值区间的下限时， 为设备正常状态， 若H_I值处于区间[H_{\mathrm{I}\mathrm{0}}-\delta -\mathrm{0.5}+\gamma ,H_{\mathrm{I}\mathrm{0}}+\delta +\mathrm{0.5}-\gamma ]时则进入预警阶段， 通过后续输入的数据进行判断.

3.2 健康阶段自主划分策略

滚动轴承从正常到完全失效的性能退化过程一般分为4个阶段：正常阶段、轻微故障阶段、故障增长阶段、严重故障阶段^［11］. 本文对全生命周期历史数据进行健康阶段划分（health stage， HS）时， 所采用的方法流程如下：

1）计算接触力：将每组振动信号输入动力学模型中进行计算， 得出接触力退化数据.

2）H_I计算：将计算完成的接触力， 按照式（21）计算方差， 得到H_I数据.

3）确定基准阈值：确定当前设备的工况， 借鉴文献［17］中对故障尺寸的界定方法， 分别设置对应不同退化阶段之间临界点的故障尺寸， 求解动力学模型， 计算得到H_I0、H_I1、H_I2， 将H_I0、H_I1、H_I2作为健康阶段划分的基准阈值.

4）设置区间范围：以H_I0、H_I1、H_I2作为基准阈值， 参考文献［18］的思路， 引入拉依达准则， 将故障阈值表示为动态区间［H_{\mathrm{I}i}-\delta ,H_{\mathrm{I}i}+\delta］（i=0，1，2）， 则健康状态的阶段H_S可表示为式（22）.

（22）

5）动态更新区间：将实时振动信号输入至动力学模型， 计算出仿真信号， 计算实时信号与仿真信号的余弦相似度^［19］， 记作γ. 通过γ修正故障阈值， 当相似度越高时， 说明此次计算的H_I值可靠性更高， 不确定性的阈值区间范围可以缩小， 反之， 可以扩大. 相似度的计算结果在0~1范围内， 假设当相似度为0.5时， 保持阈值区间不变， 其计算式可表示为[H_{\mathrm{I}i}-\delta -\mathrm{0.5}+\gamma ,H_{\mathrm{I}i}+\delta +\mathrm{0.5}-\gamma ] （i=0，1，2）. 经修正后， H_S可表示为式（23）所示.

H_{\mathrm{S}}=（23）

6）不确定性样本的处理：部分数据样本计算H_I后， H_I值会处于[H_{\mathrm{I}i}-\delta -\mathrm{0.5}+\gamma ,H_{\mathrm{I}i}+\delta +\mathrm{0.5}-\gamma ]（i=0，1，2）区间内， 导致无法判断该样本处于哪个阶段. 针对此类情况， 可计算当前数据前面一组和后面一组数据的H_S. 若前后两组数据H_S一致， 则用前后两组的H_S代替当前组数据的H_S， 若前后两组数据H_S不一致， 则比较前后两组的γ， 以γ最高的H_S代替当前组数据的H_S.

7）提取临界点：根据上述步骤， 判断H_I曲线的所属范围， 提取出每个阶段的第一组数据作为每个退化阶段的起始点， 最后实现健康阶段划分.

4 实验结果

通过西安交通大学^［20-21］、美国国家航空航天局^［22］及实验室自建平台等得到的滚动轴承振动实验数据集来验证所提方法.

4.1 自建平台数据集实时健康监测实验

为验证本研究实时健康监测方法的有效性， 自主搭建了实验平台， 如图11所示， 相关参数如表5所示. 同时， 为模拟轴承在实际运行中的情况， 设置转速为1 797 r/min， 载荷为1 kN， 在轴承的水平和垂直方向上安装两个PCB 352C33加速度传感器， 获取轴承的振动信号. 采样频率设置为12 kHz， 采样间隔为23 s， 每次采样时长为3 s.

输入每组数据进行实验（后文以第1组故障数据为例进行实验流程说明），如图12所示. 设置动力学模型故障尺寸为0， 求解动力学模型， 得到仿真信号［图12（b）蓝色信号］， 同时得到H_I值为14， 记作H_I0. 以H_I0作为基准阈值， 则故障阈值可表示为［12.6， 15.4］. 计算实时数据的均方根、方差、最大值3个特征， 将结果输入CNN， 得出动力学模型参数. 基于所得参数求解动力学模型， 计算出仿真信号［图12（a）红色信号］， 同时得到H_I值为18， 计算实时信号［图12（a）蓝色信号］， 计算仿真信号的相似度γ为0.9， 通过γ来修正故障阈值， 修正后的故障阈值可表示为［13.0， 15.0］. 由于当前信号的H_I值大于阈值范围， 准确地触发了故障报警.

选取500组数据进行实验， 统计实验结果， 按照式（24）~式（26）分别计算各个评估指标， 得到的结果如表6所示.

\mathrm{A}\mathrm{c}\mathrm{c}\mathrm{u}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{y}=\frac{\mathrm{T}\mathrm{P}+\mathrm{T}\mathrm{N}}{\mathrm{T}\mathrm{P}+\mathrm{T}\mathrm{N}+\mathrm{F}\mathrm{P}+\mathrm{F}\mathrm{N}}（24）

\mathrm{P}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}=\frac{\mathrm{T}\mathrm{P}}{\mathrm{T}\mathrm{P}+\mathrm{F}\mathrm{P}}（25）

\mathrm{R}\mathrm{e}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{l}=\frac{\mathrm{T}\mathrm{P}}{\mathrm{T}\mathrm{P}+\mathrm{F}\mathrm{N}}（26）

式中：TP为真正例数， 即实际为故障且被正确识别为故障的样本数；TN为真负例数， 即实际为正常且被正确识别为正常的样本数；FP为假正例数， 即实际为正常但被误判为故障的样本数；FN为假负例数， 即实际为故障但被误判为正常的样本数；Accuracy、Precision和Recall分别表示准确率、精确率和召回率.

为了验证本文实时健康监测方法的效率，统计了上述实验过程中每组数据计算所用的时间，平均计算时间约为0.14 s，小于数据采集时间0.25 s，可满足在线运行的要求.

在上述实验基础上，对比了均方根等指标在健康监测中的性能，得到如表7所示的结果.本文方法相对于CNN+振动信号方差的方法计算时间成本更高， 但准确率更高，这是由于本文方法进行了接触力计算，增加了模型计算时间，但是融合接触力特征之后，模型能提取到更丰富的滚动轴承退化特征.与BP+接触力方差的方法对比，两种方法都采用了接触力计算的方式，但模型更新机制不同， 本文方法计算时间降低了6.7%，且准确率提高了7%.这主要得益于本文方法采用了CNN进行模型更新，CNN相较于反向传播神经网络（back propagation neural network， BP神经网络）参数辨识和提取能力更强，能更准确地辨识模型参数，从而提高模型更新效率.因此，从实验结果上可看出本文的实时健康监测方法准确率较高，计算速度较快，可满足在线运行的应用场景需求.

4.2 XJTU-SY数据性能退化评估实验

验证性能退化评估方法的数据来自西安交通大学（XJTU-SY）数据集^［20］，实验中交流电机提供动力，运转过程转速处于非平稳状态^［21］，实验中设置采样频率为25.6 kHz，采样间隔为1 min，每次采样时长为1.28 s.以XJTU-SY数据集的Bearing2_5数据进行实验，样本总数339， 实际寿命5 h 39 min.设备如图13所示，设备相关参数如表8所示.

以XJTU-SY数据集的Bearing2_5数据（图14蓝色信号）进行实验， 将完全失效时间（第339组）作为终止故障时间. 根据3.2节中的方法进行计算， 得到全生命周期仿真信号（图14红色信号）， 同时得到第（121，165，184）组数据为健康阶段划分的分界点. 根据分界点进行阶段划分， 得到图15的结果. 如图15所示，可清楚观察到轴承退化的阶段过程：正常状态， 轻微故障， 故障增长， 严重故障.

阶段划分完成后， 通过提取第120组和第121组的故障特征来验证阶段划分的有效性. 由于实验数据采集过程中处于非平稳工况状态， 信号中存在噪声，且初始故障特征比较微弱，会影响故障特征的提取.因此，在提取特征前采用全变分降噪（total variation denoising， TVD）方法进行降噪^［23］，再通过Teager能量算子解调谱提取故障特征，如图16所示. 从图16可知仅在第121组数据中提取到与轴承外圈故障相近的特征，由此可确定初始故障点为121组.

为进一步验证所提阶段划分的合理性， 通过时间序列突变点检测对本文的H_I指标序列进行检测^［24］， 获取突变点， 两种方法的划分结果一致， 佐证了本文划分方法的有效性， 且本文方法计算流程比时间序列突变检测方法简单， 计算速度更快.

4.3 NASA数据性能退化评估实验

为验证性能退化评估方法的效率，采用了美国国家航空航天局（NASA）的轴承数据进行验证分析^［22］.实验中使用NASA数据的长度为20 480个点，采样频率为 20 kHz， 采样点为16 384，驱动电机转速为2 000 r/min，载荷26.6 kN，共984组数据.NASA轴承全寿命周期实验平台如图17所示，设备参数如表9所示，重复上述方法，得到的仿真信号和实际信号如图18所示.

以984组作为完全失效时间，进而计算得到3个转折点（533， 614， 974），3个转折突变点将指标分为呈梯度上升的4个状态.可将滚动轴承的阶段划分为如图19所示的4个阶段：正常状态、轻微故障、故障增长、严重故障.

由于性能退化评估方法的精度和速度受多种因素的影响，因此建立不同的评价指标对模型进行评估，具体包括准确率、占用内存大小、浮点计算量. 准确率表征的是模型计算精度，占用内存大小表征的是模型在占用硬件设备的内存大小，内存占用越小，模型越能部署在轻量化系统中，浮点计算量表征模型的计算速度，计算量越小，模型计算速度越快.由于本文的所有模型和方法都是在Python环境中实现，因此可通过timeit模块直接查看.对比文献［11-13］， 得到表10的结果.

由表10可知，本文方法各项指标均优于对比方法.相较于准确率最高的文献［13］，本文准确率提升1.67%；相较于内存占用最小的文献［11］，本文内存占用降低30.4%；相较于浮点计算量最小的文献［12］，本文计算量降低4.55%. 这主要得益于本文的性能退化评估方法引入了模型同步更新的机制，且引入了相较于振动信号特征更准确的接触力方差，因此更适用于集成到实际的工业健康管理平台中，可满足轻量化运行的要求.

5 结 论

为提升轴承性能退化评估方法的准确性和快速性， 本文提出一种融合动力学模型的滚动轴承性能退化评估方法，主要贡献如下：

1）变转速工况下的高精度动力学建模.本研究在现有动力学模型基础上，创新性地引入了滑动摩擦力修正项，使模型能准确表征不同转速下轴承的动态响应.相较于传统模型，本模型对故障特征频率的相对误差从1.849%降低到0.097%，显著提高了建模精度.

2）基于CNN的轻量级模型同步更新机制.与传统的参数辨识方法相比，所提的CNN映射方法实现了动力学模型的实时更新.相比BP神经网络，所提方法计算时间降低6.7%，准确率提高7%，体现了在计算效率和准确性方面的双重优势.

3）接触力特征驱动的健康评估新范式.不同于直接使用振动信号特征，本文创新性地利用动力学模型计算的接触力特征构建健康指标.与主流方法相比，本方法准确率提升1.67%，内存占用降低30.4%，浮点计算量降低4.55%， 特别适合边缘计算和资源受限场景.

4）相似度驱动的自适应健康评估策略.本研究提出的基于信号相似度动态修正阈值的自适应评估策略，有效解决了固定阈值导致的误报和漏报问题. 该策略实现了91%的监测准确率，为在线健康监测提供了新思路.

综上，本研究为变转速工况下的轴承健康评估领域提供了一个新的解决方案，本文方法不仅在建模精度上有所提升，更在计算效率和实用性方面展现了巨大的优势，为工业设备健康监测技术向智能化、轻量化方向发展提供了一定的理论基础. 在未来的工作中，我们将温度和润滑条件等因素纳入动力学模型， 以构建更全面的多物理场耦合模型， 提高模型的适用性及泛化性.

参考文献

[1] 张小丽， 和飞翔， 梁旺， 等. 基于GADF和CWT并行输入模型的滚动轴承智能诊断研究［J］. 湖南大学学报（自然科学版），2025， 52（2）： 98-108．

[2] ZHANG X L， HE F X， LIANG W，et al ．Research on intelligent diagnosis of rolling bearings based on parallel input model of GADF and CWT［J］. Journal of Hunan University （Natural Sciences）， 2025， 52（2）： 98-108．（in Chinese）

[3] 康守强，周月， 王玉静， 等. 基于改进SAE和双向LSTM的滚动轴承RUL预测方法［J］. 自动化学报，2022，48（9）： 2327-2336．

[4] KANG S Q， ZHOU Y， WANG Y J，et al ．RUL prediction method of a rolling bearing based on improved SAE and Bi-LSTM［J］．Acta Automatica Sinica，2022， 48（9）： 2327-2336．（in Chinese）

[5] 李宝庆， 杨路航， 杨宇， 等 ．基于旋转多面体锥的滚动轴承智能故障诊断［J］. 湖南大学学报（自然科学版）， 2022，49（6）：55-64．

[6] LI B Q， YANG L H， YANG Y，et al ．Rolling bearing intelligent fault diagnosis based on rotated and extended polyhedron cone［J］. Journal of Hunan University （Natural Sciences）， 2022， 49（6）：55-64．（in Chinese）

[7] WU R Q， WANG X F， NI Z X，et al ．Dual-impulse behavior analysis and quantitative diagnosis of the raceway fault of rolling bearing［J］．Mechanical Systems and Signal Processing，2022，169： 108734.

[8] 田晶， 艾辛平， 刘丽丽， 等 ．中介轴承复合故障动力学建模与振动特征分析［J］．振动与冲击， 2022， 41（22）： 144-151．

[9] TIAN J， AI X P， LIU L L，et al ．Dynamic modeling and vibration characteristic analysis of the inter-shaft bearing’s multiple point fault［J］．Journal of Vibration and Shock，2022，41（22）：144-151．（in Chinese）

[10] JIANG Y C， HUANG W T， LUO J N，et al. An improved dynamic model of defective bearings considering the three-dimensional geometric relationship between the rolling element and defect area［J］. Mechanical Systems and Signal Processing，2019，129：694-716.

[11] 罗茂林， 郭瑜， 伍星. 考虑冲击力的球轴承外圈剥落缺陷双冲击现象动力学建模［J］. 振动与冲击， 2019， 38（14）： 48-54．

[12] LUO M L， GUO Y， WU X ．Dynamic modeling of the dual-impulse behavior produced by a spall on the outer race of a ball bearing considering impact forces［J］. Journal of Vibration and Shock，2019，38（14）： 48-54．（in Chinese）

[13] 胡爱军， 吉新星， 向玲， 等. 滚动轴承非线性时变参数动力学模型与故障机理研究［J］. 机械工程学报，2022，58（19）：139-147．

[14] HU A J， JI X X， XIANG L，et al. Nonlinear time-varying parameter dynamic model of rolling bearing and failure mechanism research［J］．Journal of Mechanical Engineering，2022，58（19）：139-147．（in Chinese）

[15] YAO R H， JIANG H K， YANG C X，et al ．An integrated framework via key-spectrum entropy and statistical properties for bearing dynamic health monitoring and performance degradation assessment［J］. Mechanical Systems and Signal Processing，2023，187：10995.

[16] PAN Y N， CHENG D L， WEI T T，et al ．Rolling bearing performance degradation assessment based on deep belief network and improved support vector data description［J］. Mechanical Systems and Signal Processing，2022，181：109458.

[17] 张龙，张号，周建民，等 ．采用显式动力学的轴承性能退化评估指标构建［J］．西安交通大学学报， 2022， 56（8）：11-21．

[18] ZHANG L， ZHANG H， ZHOU J M，et al ．Construction of bearing performance degradation assessment indicator using explicit dynamics［J］．Journal of Xi’an Jiaotong University，2022，56（8）：11-21．（in Chinese）

[19] WANG X， CUI L L， WANG H Q，et al ．A generalized health indicator for performance degradation assessment of rolling element bearings based on graph spectrum reconstruction and spectrum characterization［J］．Measurement，2021，176：109165.

[20] 赵聪聪， 刘玉梅， 赵颖慧， 等. 基于MMFE和可拓k-medoids聚类的轴承性能退化评估［J］. 振动与冲击，2022，41（17）：123-130．

[21] ZHAO C C， LIU Y M， ZHAO Y H，et al ．Evaluation of bearing performance degradation based on MMFE and extensible k-medoids clustering algorithm［J］．Journal of Vibration and Shock，2022，41（17）：123-130．（in Chinese）

[22] 王震，杨正伟，何浩浩， 等 ．非规则轴承故障的动力学建模与仿真［J］．北京航空航天大学学报， 2021， 47（8）： 1580-1593．

[23] WANG Z， YANG Z W， HE H H，et al ．Dynamic modeling and simulation of irregular bearing failure［J］. Journal of Beijing University of Aeronautics and Astronautics，2021，47（8）：1580-1593．（in Chinese）

[24] 李昊泽， 贺雅， 冯坤，等 ．考虑时变激励的滚动轴承局部故障动力学建模［J］．航空学报， 2022， 43（8）： 625176．

[25] LI H Z， HE Y， FENG K，et al. Dynamic modeling of rolling bearing local fault considering time-varying excitation［J］. Acta Aeronautica et Astronautica Sinica， 2022， 43（8）： 625176．（in Chinese）

[26] 孟良，苏元浩，许同乐，等 ．并行卷积神经网络的风电机组故障诊断方法［J］．太阳能学报，2023，44（5）：449-456．

[27] MENG L， SU Y H， XU T L，et al. Wind turbine fault diagnosis method based on parallel convolutional neural network［J］. Acta Energiae Solaris Sinica，2023， 44（5）： 449-456．（in Chinese）

[28] CUI L L， WANG X， WANG H Q，et al ．Remaining useful life prediction of rolling element bearings based on simulated performance degradation dictionary［J］. Mechanism and Machine Theory，2020，153：103967.

[29] XIA J H， ZHANG J B， WANG Y，et al. WC-KNNG-PC：watershed clustering based on k-nearest-neighbor graph and Pauta Criterion［J］. Pattern Recognition，2022，121：108177.

[30] 向家伟 ．数值模型驱动的传动系统故障个性化诊断原理［J］．机械工程学报，2021，57（15）：116-128．

[31] XIANG J W. Numerical model driving personalized diagnosis principle for fault detection in mechanical transmission systems［J］. Journal of Mechanical Engineering，2021，57（15）：116-128．（in Chinese）

[32] WANG B， LEI Y G， LI N P，et al ．A hybrid prognostics approach for estimating remaining useful life of rolling element bearings［J］．IEEE Transactions on Reliability，2020，69（1）：401-412.

[33] 雷亚国，韩天宇，王彪，等 ．XJTU-SY滚动轴承加速寿命试验数据集解读［J］．机械工程学报，2019，55（16）：1-6．

[34] LEI Y G， HAN T Y， WANG B，et al ．XJTU-SY rolling element bearing accelerated life test datasets：a tutorial［J］. Journal of Mechanical Engineering，2019，55（16）：1-6．（in Chinese）

[35] QIU H， LEE J， LIN J，et al. Wavelet filter-based weak signature detection method and its application on rolling element bearing prognostics［J］. Journal of Sound and Vibration，2006，289（4/5）：1066-1090.

[36] 王凡，马军，王晓东， 等 ．一种改进的TVD电机轴承故障特征提取方法［J］．振动与冲击， 2023， 42（10）： 203-214．

[37] WANG F， MA J， WANG X D，et al ．An improved TVD fault feature extraction method for motor bearing［J］. Journal of Vibration and Shock， 2023， 42（10）： 203-214．（in Chinese）

[38] 刘弹，李晓婉，梁霖，等 ．采用时间序列突变点检测的滚动轴承性能退化评价方法［J］．西安交通大学学报，2019，53（12）：10-16．

[39] LIU D， LI X W， LIANG L，et al ．A method for evaluating performance degradation of rolling bearings using detection of time series mutation point［J］. Journal of Xi’an Jiaotong University，2019，53（12）：10-16．（in Chinese）

基金资助

国家自然科学基金资助项目(62163020)

国家自然科学基金资助项目(62173168)

National Natural Science