梁洁琳，丁 康，何国林，2，林慧斌，蒋 飞

（1.华南理工大学机械与汽车工程学院，广东 广州 510640；2.琶洲实验室，广东 广州 510335）

引言

滚动轴承作为旋转机械的关键部件，其退化状态影响机械设备的安全可靠性能。准确的轴承寿命预测有利于预防性维护系统，避免轴承突然失效而导致设备发生致命故障［1-4］。其中，运行状态划分和剩余使用寿命预测是轴承寿命预测框架中的两个重要问题。

目前常用的轴承正常运行和损伤运行阶段划分的方法为通过合理的阈值准则选取退化初始时间（DIT）［5］。Niu 等［6］应用国际规范的轴承阈值基准线和报警系数获取融合指标的DIT。申彦斌等［7］提出一种基于卷积自编码器的轴承健康指标（Health Indicator，HI）构建方法，根据所提HI 的变化规律人为划分三个运行阶段；Hong 等［8］利用自组织映射网络构成全寿命退化评估曲线并手动设定正常和早期退化的阈值界限。然而，以上方法均基于大量统计特性或数据经验主观选取，没有进一步验证DIT 结果的合理性，对于单个系统而言，所选DIT 可能会滞后或提前，导致检修延迟或误预警。

基于数据驱动的轴承寿命预测方法主要包括机器学习和统计模型两种方法。前者结合大量可用数据样本预测退化趋势，适用于处理复杂机械系统的RUL 预测问题［9-10］。但受限于成本和时间耗费，实际应用中难以获取高质量数据样本，增加了RUL 预测难度。而统计模型预测方法的优势在于同时结合先验数学［11-12］或物理［13-14］模型以及最新数据样本以更新未知参数，相对简单易行。其中，基于指数模型的传统粒子滤波算法广泛应用于锂电池、飞行器等［15-16］不同对象的RUL 预测领域。滚动轴承的全寿命退化过程也与指数模型相似，但由于受损伤演变特性［17］或随机误差的影响，表征轴承退化过程的观测指标会局部偏离全局指数模型，呈上下非线性波动，导致传统粒子滤波的预测曲线偏离实际退化趋势，降低寿命预测精度。

综上，本文提出了联合频域相关分析和改进粒子滤波算法驱动的滚动轴承寿命预测方法。基于滚动轴承在退化过程中振动响应信号频域特征具有短期相似和长期差异这一特性，构建平均相关系数曲线，识别合适的早期退化时间，并结合小波降噪和希尔伯特包络解调方法验证所选DIT 的准确性。由于观测值的局部波动，传统粒子滤波的粒子后验分布也会相应偏离真实分布，导致RUL 预测精度降低。因此构建了一个同时考虑全局指数式退化趋势与局部波动双重因素的改进粒子滤波算法。所提方法通过新构造的粒子约束器校正粒子分布，引导粒子的后验概率密度朝真实退化值更新，实现更高精度的寿命预测。

1 ACC-DFPF 滚动轴承寿命预测方法

1.1 基于频域相关分析的ACC 区间划分

均方根（Root Mean Square，RMS）具有较好的单调性，且能有效反映振动能量，该时域特征常作为轴承的健康指标。然而，RMS 在早期故障时间前后的幅值变化不明显，难以直接准确识别DIT。实际上，故障特征频率作为轴承发生故障的机理特性之一，频域幅值和频率倍频会在轴承不同的退化状态存在差异。基于此，通过构建平均相关系数（ACC）曲线和所设阈值识别RMS 的DIT，具体步骤如图1所示。

图1 滚动轴承状态区间划分流程Fig.1 Process of rolling bearing state division

a）对滚动轴承的全寿命振动响应时域信号做滑窗采样，并相应进行快速傅里叶变换得到幅值谱Pk(k=1，2，…，T)，T为幅值谱总个数；

b）运用Pearson 相关系数评价任意两幅值谱的相关程度，如下式：

构建相关系数矩阵M：

式中ρst为第s和第t个幅值谱之间的相关系数，s，t=1，2，…，T；Ps，n，Pt，n分别对应第s和第t个幅值谱的第n条谱线频率值和分别为 每个幅 值谱的 谱线数，第s和第t个幅值谱的谱线频率均值。

c）为直观地反映起始健康信号与后续每段信号的相关程度变化，计算一维ACC 曲线，计算式为：

结合3σ 阈值识别DIT，最后用小波降噪和包络解调谱验证结果的准确性。

1.2 传统粒子滤波(PF)算法框架

粒子滤波通过蒙特卡罗方法（Monte Carlo，MC）避免贝叶斯递推过程中的高维积分计算，得到粒子后验概率分布和样本均值，从而估计非线性系统的状态。滚动轴承的退化趋势常用指数模型描述，其时变系统的状态方程和观测方程分别为：

式中uk和νk分别对应k时刻未知模型参数ak和bk的过程误差；yk和xk分别为k时刻的系统观测指标和模型状态值，两者存在观测误差χk；uk，νk，χk分别服从标 准差为σv，σq，σr的正态分布

传统粒子滤波算法基于式（4）进行预测和更新。

a）预测过程。已知第1 至k-1 时刻的系统观测值y1∶k-1和参数向量Θk-1，采用一组粒子数为Np的粒子描述k时刻参数的先验概率分布P(Θk|Θk-1，y1∶k-1）。基于式（4），第j个粒子的参数预测向量和模型状态值分别为：

b）更新过程。第k时刻获取观测值yk时，采用下式更新和归一化每个预测粒子的权重

此时输出k时刻的指标估计值

1.3 改进粒子滤波驱动的轴承RUL 预测

由式（6）可知，PF 算法的粒子权重更新主要依赖最新观测值。当输入的最后观测值偏离全局退化值时，PF 算法后续预测的轴承退化轨迹会偏离真实趋势，如图2 所示。A1和分别为实际观测值和全局模型值，当A1＞时，PF 预测曲线比真实趋势更早到达失效阈值，高估了实际故障严重程度；反之预测值比真实趋势更平缓，以上两种局部波动观测值均降低了传统PF 算法的轴承RUL 预测精度。

图2 传统粒子滤波的预测趋势Fig.2 Trends predicted by conventional particle filter

因此，在原有PF 方法的更新过程中，新构造一个同时考虑全局模型参数和局部观测值的约束器，形成改进粒子滤波。所提约束器引导粒子分布朝向全局退化值更新，抑制局部波动观测值对轴承HI 预测结果的干扰。具体步骤为：

1）构建k时刻的粒子损失函数：

2）k时刻的权重系数λk∈[0，1]表示为：

式中c和δ均为常 数因子。令x0，k=a0exp(b0k)为k时刻的全局模型估计值。

3）第j个粒子的参数预测向量校正为：

利用Θ0和yk计算损失函数，目的是通过平滑的全局模型趋势和输入观测值同时约束粒子分布，减弱输入样本对后续预测曲线的影响。此外，依据yk和x0，k的差异程度，权重系 数λk可自适 应调整式（10）的两项权重。

综上，所提ACC-DFPF 滚动轴承寿命预测方法流程如图3 所示，具体说明如下：

图3 ACC-DFPF 滚动轴承寿命预测方法流程Fig.3 Flowchart of ACC-DFPF method for rolling bearing life prediction

a）采集轴承数据，对各段频谱进行相关分析，构建ACC 曲线，结合所设阈值识别DIT，完成HI 指标的正常和损伤状态划分，并依据小波降噪包络解调谱的故障频率成分验证DIT 的准确性。

b）采用最小二乘法拟合轴承损伤区间的HI 指标，计算全局模型经验参数Θ0=[a0b0]Τ，Θ0也作为DFPF 算法的初始化参数输入DFPF 算法。当k时刻获取yk时，采用式（10）～（12）所提约束器校正粒子预测过程的先验分布，随后进行粒子权重更新和重采样，直至再无观测值输入时，粒子校正停止。

c）继续采用式（5）和式（8）计算，输出轴承从最后观测值时刻kh至预测失效时刻的HI 指标估计值和估计 寿命值

2 试验验证

2.1 ACC 划分区间方法的分析对比

选择XJTU-SY 滚动轴承全寿命试验数据集［19］验证所提方法的有效性。测试平台的采样频率为2.56 kHz，每间隔1 min 采样1.28 s，在3 类工况下各重复5 次试验，共计15 组轴承全寿命时间序列样本。其中，bearing3-1 发生外圈磨损。

以轴承bearing3-1 为例，所得ACC 退化曲线和DIT 点如图4 所示。初始运行时系统不稳定，对应ACC 值产生幅度波动。健康状态阶段的轴承平均相关系数相对平稳，随后在初始故障位置开始减小，并在损伤阶段随退化程度加深而整体减小。图4中，ACC 值在500～1500 min 内相对稳定，由此得到3σ 阈值下限为0.6106，在 第2371 min 首次超过阈值，表明此时该轴承发生初始故障。

图4 ACC 识别起始故障点结果Fig.4 The DIT identification result with ACC model

进一步采用小波降噪和包络解调验证所选DIT的准确性。考虑早期轴承故障特征微弱，先用小波降噪减弱信号样本的噪声干扰。由图5 降噪包络谱可知，第2369 min 未出现故障特征频率成分，而第2371 和第2373 min 均出现了与外圈故障特征频率理论值123.32 Hz 接近的频谱成分，且第2530 min轴承故障程度严重，故障频率与倍频的幅值都增大，验证了轴承在第2371 min 存在早期故障。

图5 不同时刻的轴承包络谱Fig.5 Envelope spectrums of bearing at different moments

此外，将所提方法与时域特征RMS、峭度kurtosis 和频域特征频谱重心（Frequency Center，FC）识别DIT 方法的结果进行对比。bearing3-1 的RMS 在第2340 min 后均超过上限阈值，如图6 所示，说明RMS 比ACC 方法提早31 min 推断轴承发生早期故障。但图7 中的第2340 和2350 min 降噪解调谱未有明显的故障特征频率谱峰。此外，RMS整体比ACC 曲线振荡，且在第2340 min 前有多处异常RMS 超过阈值造成DIT 的错误预警。

图6 RMS 识别起始故障点结果Fig.6 The DIT identification result with RMS model

图7 第2340 和2350 min 的降噪包络谱Fig.7 Denoising envelope spectrum at 2340 and 2350 min

峭度方法识别DIT 的结果如图8 所示，在第2523 min 后均超过上限阈值，推断该时刻发生了轴承初始故障，比ACC 方法延迟152 min，对应检验了包络谱有突出的故障频率谱峰。此外，在第2371 min 前有多处超过所设阈值的剧烈波动峭度，和RMS 方法存在类似的误预警现象。

图8 Kurtosis 识别起始故障点结果Fig.8 The DIT identification result with kurtosis model

频谱重心识别DIT 如图9 所示。分别求取［0，2000］ min 和［2400，2454］ min 上升和 下降区间的FC 值斜率，交叉点位于第2346 min，初步推断该时刻发生轴承初始故障，相比ACC 方法提前25 min识别了DIT。此外，该轴承的FC 指标相较于ACC和以上时域特征整体单调性较差。

图9 FC 识别起始故障点结果Fig.9 The DIT identification result with FC model

综上，ACC 方法可以比时域RMS，kurtosis 和频域FC 方法更准确地识别出轴承初始故障点，且ACC 曲线整体更平滑，避免错误预警早期故障点的情况，为下一步进行轴承剩余寿命预测提供起始点依据。

2.2 DFPF 方法预测轴承RUL 的分析对比

选取包括同一尺寸3 个工况的bearing1-3，bearing2-1 和bearing3-1 作为算法验证数据集。结合3 类运行工况共6 组历史轴承的退化规律选取统一的失效阈值。对RMS 健康指标进行最大、最小值归一化，避免了不同工况下因轴承HI 数值差异较大导致失效阈值γ难选择的问题，如图10 所示。图10（b）的RMS 值在0.6 以上均变化陡峭，表明该阈值后轴承加速劣化至完全损坏，因此设定失效阈值γ=0.6。

图10 6 组轴承数据的HI 曲线Fig.10 Six sets of bearing health indicator curves

基于ACC 方法和γ确定表1 验证数据集的DIT和实际失效时刻kT，采用最小二乘法拟合对应HI 指标值得到全局经验模型参数a0，b0以及标准差σq，σr和σv。其中，初始权重系数取λ0=1 表示没有观测值输入时，粒子的校正完全依赖所得全局经验模型参数。DFPF 算法的主要结构参数如表2 所示。

表1 验证数据集的DIT 和失效时刻Tab.1 DIT and actual degradation time of the testing data

表2 DFPF 算法主要结构参数Tab.2 The main structural parameters of DFPF algorithm

分别输入特定观测值比例，采用所提DFPF 算法得到一系列指标预测值，与PF 算法结果相比，如图11 所示。

图11 DFPF 和PF 算法所得性能退化预测曲线Fig.11 Performance degradation prediction curves obtained by DFPF and PF algorithm

对于轴承bearing2-1，50%的样本输入比例对应最后输入观测值的时刻约为第469 min，该局部观测数据可能受测量过程的噪声随机影响或是退化演变波动特性，偏离了全局经验模型。由于PF算法输出的预测值主要依赖于最后输入的观测值，该预测失效时刻比实际失效时刻提早12 min，而所提DFPF 方法更为有效，仅提早了5 min。原因在于所提DFPF方法有全局模型参数和局部观测值这两项共同约束粒子分布，降低了后续预测结果对局部波动数据的敏感性，使得轴承未来状态的预测值比传统PF 算法更趋近真实观测曲线。同理对于轴承bearing1-3 和bearing3-1，DFPF 预测曲线的95%置信区间间距更窄，RUL 相对误差更小，具体如表3 所示。由以上分析可知，所提DFPF 算法预测轴承剩余寿命精度优于PF 算法。

表3 特定时间段下轴承RUL 预测误差结果Tab.3 Results of bearing RUL prediction error under a specific time period

为了更直观地说明不同时间段下所提DFPF 算法的预测优势，按5%，10%，15%，…，95%共19 组不同的观测值输入比例预测轴承剩余使用寿命，如图12 所示。

图12 不同起始预测时刻的轴承RUL 曲线Fig.12 RUL curves of bearings at different initial prediction moments

此外，分别用均方根误差Ermse、平均绝对百分比误差Emape和平均绝对误差Emae指标评价两种算法的RUL 预测精度，具体计算公式为：

以上指标越小代表预测精度越高，具体结果对比如表4 所示。

表4 不同起始预测时刻的轴承RUL 误差Tab.4 RUL errors for bearings at different initial prediction moments

当输入观测值比例为较少的5%时，轴承RUL预测结果对应输出图12（a）～（c）的左起第一个点，可直观看出其接近实际值。其原因可能是预测结果受初始模型参数影响，而以上轴承的初始模型参数均由对应的全局指标数据拟合得出，准确度较高，因此即使观测样本较少，RUL 预测结果仍会接近真实寿命。此外，对于bearing2-1 和bearing1-3，占比为10%至30%对应的预测RUL 也接近真实值，进一步分析原因可能是该比例下的实际观测值波动较小，整体较接近全局经验模型。以上分析可知在输入观测值个数不足的条件下，观测值和初始模型参数会影响轴承RUL 预测。

3 结论

（1）构造了表征轴承退化性能的ACC 曲线，准确地定位了轴承起始故障点，并采用小波降噪和希尔伯特包络解调方法辅助验证所选点的早期故障情况。由XJTU-SY 轴承数据集的验证结果对比可知，所提ACC 方法相较于均方根和峭度方法能够更准确找到早期故障点，实现轴承正常与损伤两个状态阶段的有效划分。

（2）所提DFPF 算法构造了含有全局模型参数和局部实际观测值的约束器用于校正粒子分布，克服传统粒子滤波算法受观测值波动干扰的缺陷。XJTU-SY 轴承数据集的试验结果表明，所提DFPF算法后续预测的健康指标值比PF 算法更接近全局运行状态，置信区间宽度和相对误差Erh均更小，实现了更准确的轴承寿命估计。

（3）输入不同数目观测数据的试验结果表明，DFPF 算法可以有效减少轴承局部异常波动观测数据对轴承RUL 预测结果的影响，其RUL 的Ermse，Emape和Emae误差均小于PF 算法，实现了更可靠的寿命预测。此外，受限于试验样本的数量，全局经验模型参数的调整缺少更多全局退化趋势相似的轴承离线观测数据，而目前在实际应用中难以获取。今后如何根据有限的测量数据获取该参数信息，是寿命预测研究中需要重点考虑的问题之一。