简介面向工业物联网与设备健康管理场景的轴承寿命预测MATLAB代码包适合从事故障诊断、预测性维护研究的工程师及高校学生。压缩包仅1KB包含3个m文件代码结构精炼覆盖数据读取、预处理、时域特征提取与寿命预测建模等关键环节可帮助读者快速上手轴承剩余寿命分析流程。内容结合轴承振动信号的均方根、峭度、幅值等时域变换指标演示如何识别早期故障模式并预测失效风险对理解设备健康管理与工业4.0智能决策具有直接参考价值。资源已有1134人学习轻量实用适合作为算法验证与教学示例。1. 轴承寿命预测中的时域变换先压缩信号再追踪退化轴承寿命预测中最容易走的弯路是直接拿原始振动波形去拟合退化趋势。原始时域信号每秒有数千甚至数万个采样点对完整生命周期而言既有退化趋势也夹杂大量与故障无关的噪声和工况干扰。时域变换要解决的就是在保留退化信息的前提下压缩维度把振动信号变换成物理意义明确的统计量或包络成分再构建成健康指标曲线最终映射到剩余使用寿命。这篇以“时域变换”为主线从基础统计量提取到包络解调与 STFT再到退化模型拟合和预测验证覆盖轴承寿命预测从特征到结论的完整链路。适合正在建设设备健康管理系统PHM的工程师以及把剩余寿命预测落到产线数据上的算法开发人员。2. 时域变换第一步把振动波形转换为可追踪的健康指标2.1 六个基础时域统计量的退化敏感度对照滚动轴承在退化过程中振动信号的能量、分布形态和冲击成分都会发生规律性变化。时域变换的第一层是从原始波形中提取能反映这些变化的统计量。工程上经常成套使用均值、均方根、峭度、峰值因子、脉冲因子和裕度因子六个特征它们的退化敏感度和适用阶段差异很大。统计量计算公式退化过程中的变化趋势敏感阶段均值$\bar{x}\frac{1}{N}\sum_{i1}^{N}x_i$基本稳定严重磨损时出现偏移无明显指示均方根RMS$x_{rms}\sqrt{\frac{1}{N}\sum_{i1}^{N}x_i^2}$随退化能量增大单调上升中后期峭度$K\frac{1}{N}\sum_{i1}^{N}\left(\frac{x_i-\mu}{\sigma}\right)^4$早期点蚀冲击使峭度跳升严重退化后回落早期预警峰值因子$C\frac{x_{peak}}{x_{rms}}$早期出现冲击尖峰时升高早期脉冲因子$I\frac{x_{peak}}{\frac{1}{N}\sum_{i1}^{N}\lvert x_i\rvert}$对单个冲击更敏感早期裕度因子$L\frac{x_{peak}}{\left(\frac{1}{N}\sum_{i1}^{N}\sqrt{\lvert x_i\rvert}\right)^2}$对幅值变化放大程度最高极早期表中六个统计量可以分成两类来理解。RMS 是能量型指标反映退化过程中整体冲击能量的提升趋势稳定但到中后期才明显加速。峭度是分布形态指标对偶发冲击极其敏感外圈点蚀或滚动体剥落刚出现时峭度值能从正常的 3 左右迅速升到 6 甚至 10 以上这比 RMS 能更早捕捉到故障苗头。两者互补实际项目里很少单独使用其中一个。这个表格对部署还有一个提示不要盲目把全部统计量塞进模型。退化早期多个统计量之间的相关性很高全部放进去会产生共线性反而降低模型稳定性。我通常先做相关性分析保留彼此相关系数低于 0.7 的 2 到 3 个特征进入下游的寿命预测模型。2.2 RMS 与峭度的组合判据为什么单指标总在关键节点失效如果只用 RMS 一个指标早期微弱剥落产生的能量增量会被正常振动淹没等 RMS 出现明显抬升时轴承往往已经进入了加速退化阶段。反过来如果只用峭度轴承进入严重磨损后冲击变成连续状态信号分布趋向平稳峭度反而会从高位回落出现“退化了指标反而变好”的假象。两个单一指标都在退化链路的关键节点上失效只是失效的位置不同。常用的做法是把两者组合成双阈值判据峭度作为早期检测通道RMS 作为趋势跟踪通道。当峭度持续超过 3.5 且出现向上趋势时判定进入早期故障阶段此时不急于估算剩余寿命而是缩短采样周期、持续记录特征当 RMS 超过初始健康值的 2 到 3 倍时判定进入加速退化阶段再启动寿命预测模型。这种分级策略能规避单一指标在退化中途失效的问题也让寿命预测模型只在数据形态合适的区间内工作。2.3 用 Python 滑动窗口批量提取时域特征与参数设定实际振动采集系统给出的是一整段时域波形需要按时间窗把波形切块后逐个统计。下面是一个可以直接用于项目的滑动窗口特征提取函数import numpy as np import pandas as pd def extract_time_features(signal, fs, window_s1.0, overlap0.5): 从原始时域振动信号中滑动提取统计特征 Parameters ---------- signal : np.ndarray 原始加速度振动信号 fs : int 采样率单位 Hz window_s : float 单窗口时长单位秒 overlap : float 相邻窗口重叠率取值 0~1 Returns ------- pd.DataFrame 每个窗口的起始时间及对应的时域特征 win_len int(fs * window_s) # 一个窗口对应的采样点数 step int(win_len * (1.0 - overlap)) # 窗口滑动步长 n_windows (len(signal) - win_len) // step 1 # 能切出的窗口总数 records [] for i in range(n_windows): seg signal[i * step : i * step win_len] std np.std(seg, ddof1) rms np.sqrt(np.mean(seg ** 2)) mean np.mean(seg) kurt np.mean(((seg - mean) / std) ** 4) if std 1e-8 else 0.0 peak np.max(np.abs(seg)) records.append({ time: i * step / fs, # 窗口起始的绝对时间 rms: rms, kurtosis: kurt, crest_factor: peak / rms if rms 0 else 0.0, }) return pd.DataFrame(records) # 用法示例采样率 25600 Hz加速度信号为 sig # features extract_time_features(sig, fs25600, window_s1.0, overlap0.75)代码里值得注意的参数设置。window_s1.0意味着在 25600 Hz 采样率下单个窗口有 25600 个采样点统计量不会因为样本数太少而抖动。overlap0.75表示窗口每次只滑动 0.25 秒特征曲线的时间粒度更细代价是相邻特征点强相关需要降低计算量时改成 0.5 即可。峭度计算前先判断std是否接近零否则轴承静止或传感器未挂载时会出现除零或无穷大这是时域特征提取最常见的崩溃点。另一个容易忽略的细节是np.std默认的自由度是 0除以 N短窗口内会低估标准差窗口样本数不足 5000 时建议用ddof1。这个函数输出按时间排列的 RMS 和峭度序列正是第 4 章寿命预测模型的输入进入模型前通常还要对 RMS 做中值滤波去毛刺对峭度做去趋势处理。3. 时域信号的高阶变换包络解调与 STFT 时频分析3.1 希尔伯特变换解包络从时域冲击序列中提取故障特征频率基础统计量能告诉工程师“出问题了”却不能直接回答“故障发生在外圈、内圈还是滚动体”。这个问题需要靠故障特征频率来回答而特征频率的提取依赖对时域信号做包络变换。轴承故障会产生周期性冲击冲击激起的结构共振与轴承固有频率混叠在原始波形里。直接对原始信号做频谱分析共振峰会掩盖故障冲击形成的中低频周期成分。工程上的标准流程是先用带通滤波截取共振频带再用希尔伯特变换构造解析信号并取模得到包络波形最后对包络做 FFT。整个流程从时域冲击序列出发得到包络时域波形再得到包络频谱故障冲击的周期在包络波形里清晰可见。给定转频 $f_r$、滚动体个数 $Z$、滚动体直径 $d$、节圆直径 $D$ 和接触角 $\alpha$四类典型故障的通过频率如下故障位置通过频率计算公式外圈$f_{BPFO} \frac{Z}{2} f_r (1 - \frac{d}{D} \cos\alpha)$内圈$f_{BPFI} \frac{Z}{2} f_r (1 \frac{d}{D} \cos\alpha)$滚动体$f_{BSF} \frac{D}{2d} f_r (1 - (\frac{d}{D}\cos\alpha)^2)$保持架$f_{FTF} \frac{f_r}{2} (1 - \frac{d}{D} \cos\alpha)$在寿命预测场景里包络谱的价值不只是识别故障类型更重要的是确认故障是否处于发展期。当包络谱中出现清晰的 BPFO 谱线且幅值随运行时间单调升高说明退化在稳定推进此时启动寿命预测模型是有意义的如果谱线频繁漂移、幅值波动大说明工况不稳定预测模型输出的置信度必须下调。实际处理时还要注意许多代码直接用scipy.signal.hilbert对原始波形取包络跳过带通滤波这一步这在噪声较大的现场会失效。如果不先做带通滤波希尔伯特变换得到的是全频段包络故障冲击信息会被低频振动和高频噪声同时干扰。带通滤波的中心频率一般选在敲击试验或锤击试验得到的共振峰附近带宽取 1~2 kHz。3.2 STFT 时频图退化过程中频率成分的时域演化轨迹如果希望同时看到“什么时候出现”和“哪个频率在变化”短时傅里叶变换STFT是最直观的选择。STFT 的本质是在时域信号上滑动有限长窗对每个窗内信号做 FFT得到频率随时间变化的二维矩阵from scipy.signal import stft import numpy as np fs 25600 # 采样率 25.6 kHz nperseg 1024 # 每窗 1024 点对应 0.04 秒 noverlap 512 # 相邻窗重叠 512 点 # signal 为一段足够长的时域波形Zxx 为复数时频矩阵 f, t, Zxx stft(signal, fsfs, windowhann, npersegnperseg, noverlapnoverlap) spec np.abs(Zxx) # 幅值谱shape (len(f), len(t)) # 对频率轴求和得到每个时间点的总能量即退化趋势 energy_trend np.sum(spec, axis0) # 只取 3k~4kHz 共振频带观察其能量随时间变化 band_mask (f 3000) (f 4000) band_energy np.sum(spec[band_mask, :], axis0)STFT 有一个无法回避的权衡频率分辨率 $\Delta f fs / nperseg$时间分辨率等于窗的时长。nperseg1024在 25600 Hz 采样率下给出 25 Hz 的频率分辨率和 0.04 秒的时间分辨率分离轴承故障特征频率通常足够了因为相邻谱线间距一般大于 25 Hz。如果把nperseg提高到 2560频率分辨率变成 10 Hz但时间窗变成 0.1 秒对早期瞬态冲击的捕捉会钝化。band_energy这类频带聚合指标在实际系统中非常实用。很多系统不使用整张时频图而是把共振频带能量单独提取出来作为与 RMS 并列的退化特征。它排除了无关频段的干扰退化响应通常比全频带 RMS 更敏感。注意windowhann对振动信号是一个稳妥选择旁瓣抑制好但代价是主瓣变宽等效频率分辨率比矩形窗略差选窗时不要只看旁瓣高度。3.3 包络谱与 STFT 在不同退化阶段的使用边界包络谱和 STFT 不是替代关系而是在寿命预测链路里各管一段分析手段主要输出适合阶段典型局限包络谱故障特征频率谱线早期故障识别与定位多故障混合时谱线混淆STFT 时频图频率−时间二维能量图中期退化过程观测频率分辨率与时间分辨率矛盾包络谱序列特征频率幅值随时间的变化退化趋势定量跟踪依赖带通滤波参数误设会衰减特征在实际项目里我一般把包络谱当作“定性”工具先用它确认故障模式和特征频带再把 STFT 的频带能量当作“定量”工具把共振频带能量序列送进寿命预测模型。两个环节叠加使用比任何单一方法都更早锁定退化起点。4. 将时域特征转化为剩余寿命 RUL模型选择与参数设定4.1 从健康指标曲线到 RUL 的三种建模路线及其适用条件时域变换得到的 RMS、峭度或频带能量是随时间变化的曲线接下来的核心任务是把这条曲线映射为剩余使用寿命。根据数据量和工况复杂度有三类常用路线。第一类是经验退化模型最典型的是指数模型。轴承磨损阶段的振动能量增长近似指数形态拟合出指数参数后可以直接外推寿命。这类模型只需要几十个历史数据点就能启动部署初期数据积累不足时是最可靠的选择。第二类是纯数据驱动的时序模型比如 LSTM它需要同工况下多台轴承的全生命周期数据数据量不够时很容易过拟合优点是无需假设退化形态。第三类是混合路线先拟合指数退化模型得到退化速率参数再把这些参数作为特征输入回归模型或做查表映射物理可解释性保留得比较好。对于大多数工业现场推荐从第一类路线切入。第二类路线需要足够的历史故障样本只有少数头部企业具备这样的数据积累。下面用指数模型展开说明参数具体怎么定。4.2 指数退化模型的参数拟合与滚动更新一段可直接套用的代码常用的退化方程是 $HI(t)a e^{bt}c$其中 a 为幅值漂移b 为退化速率c 为背景噪声水平。最小二乘拟合用 SciPy 的 curve_fit 就能完成不需要引入深度学习框架import numpy as np from scipy.optimize import curve_fit def exp_model(t, a, b, c): 指数退化模型HI(t) a * exp(b * t) c return a * np.exp(b * t) c def fit_rul_model(t_hist, hi_hist, threshold, p0None): 拟合指数退化模型并计算剩余寿命 Parameters ---------- t_hist : np.ndarray 历史运行时间序列单位小时 hi_hist : np.ndarray 由时域变换得到的健康指标序列 threshold : float 失效阈值健康指标超过该值视为寿命终止 p0 : list 拟合初值默认根据数据自动估计 Returns ------- rul : float 剩余寿命预测值 perr : np.ndarray 参数估计标准差用于评估不确定性 if p0 is None: p0 [np.max(hi_hist) - np.min(hi_hist), 0.005, np.min(hi_hist)] popt, pcov curve_fit(exp_model, t_hist, hi_hist, p0p0, maxfev20000) a_fit, b_fit, c_fit popt # 求解 HI(t) threshold 对应的失效时间 if a_fit 0: t_fail np.log((threshold - c_fit) / a_fit) / b_fit else: t_fail np.inf rul t_fail - t_hist[-1] perr np.sqrt(np.diag(pcov)) return max(rul, 0), perr # 每次采集到新数据段用累计数据重新拟合 # rul, perr fit_rul_model(t_hist, hi_hist, threshold0.45)这段代码里有两个关键点需要展开。第一是threshold的取值直接影响 RUL常见做法是取健康阶段初始 RMS 的 2.5 到 4 倍或者取历史故障样本在失效前最后一个窗口的 RMS 值。阈值取低会让预测寿命偏保守取高则可能错过故障点。第二是p0初值初始化不好时curve_fit会收敛到局部极小值尤其当 a 的初值远大于实际值时。处理方式是先用前几个采样点的均值估计 c 的上界再估计 a 为健康指标的实际变化幅度最后才给 b 设定初值。perr是参数标准差可以用来近似计算 RUL 的置信区间。把 a、b、c 的误差通过失效时间公式做一次误差传播就能得到 RUL 的上界和下界。工程上通常以 1.96 倍标准差作为 95% 置信区间预测结果在区间内的概率越高运维决策越有依据。提示curve_fit默认假设残差服从高斯分布当退化曲线存在较长平台期时残差会明显偏离正态。此时改用scipy.optimize.minimize配合 Huber 损失进行鲁棒拟合比强行换模型函数更有效。4.3 滚动更新策略与关键工程参数表实际部署时不会等数据全部到位才拟合一次而是每采集到一段新数据就滚动更新模型。推荐设定一个最小拟合长度有效退化数据少于健康指标总数的 30% 时不输出 RUL只做状态监测超过 30% 后才启动模型并在每次更新后记录 b 值的稳定性。如果 b 值在连续三次更新中的浮动超过 20%说明数据噪声太大或工况不稳定此时模型给出的 RUL 不可信需要降级为趋势预警输出而不是具体数值。工程参数推荐取值说明最小拟合点数退化段窗口数的 30%不足时不输出 RUL退化速率稳定性阈值连续 3 次更新偏差 20%超限时标记低置信度更新间隔每新增 1 段窗口即更新一次指数模型拟合在毫秒级无计算压力失效阈值初始 RMS 的 2.5~4 倍或取历史失效样本的 90 分位数5. 时域变换轴承寿命预测的验证方法与三个易踩的坑5.1 用不对称评分函数检验时域变换效果验证时域变换和预测模型的效果不建议只盯 RMSE。同样大小的预测误差提前预测和滞后预测的工程代价完全不同滞后预测意味着设备已经失效而模型仍报告“还能运行”可能直接导致停机事故。工业界普遍采用不对称评分函数当预测剩余寿命大于实际值时给较小惩罚当预测值小于实际值时惩罚按指数放大。实测中先把每个样本的误差记为 d RUL_pred − RUL_true再按分段原则累积总分。同一个模型用 RMSE 看表现尚可、用不对称评分却很差时说明误差偏向滞后侧需要检查退化速率 b 是否被拟合得过于激进。5.2 数据泄漏、窗口长度与趋势漂移三个必须检查的位置第一是数据泄漏。这个坑在时域特征提取阶段出现得最多用全段信号计算全局均值和方差做标准化然后才划分训练集和测试集等于把测试集的统计信息带进了训练过程。正确做法是只在健康阶段信号上计算均值和方差保存后应用到后续全部数据。第二是窗口长度。窗口太长会平滑掉退化曲线中的加速段延后模型感知到退化的时间窗口太短又让 RMS 和峭度剧烈抖动拟合出的指数参数不稳定。判断方式直接看特征曲线的毛刺程度毛刺过多就把窗口拉长到 2~5 秒并配合中值滤波。第三是趋势漂移。传感器受温度、安装松动影响健康阶段的 RMS 基线可能整体抬升这与退化无关。处理方式是对 RMS 做差分或扣除滑动平均基线但要避免对早期退化信号过度处理否则会把有效退化趋势一起滤掉。5.3 从预测残差反推问题出在特征还是模型模型效果差时先看残差形态再决定改哪里能省去大量盲目调参时间。如果残差随时间单调增大说明退化模型函数形式偏差指数模型不适用时改试幂函数或其他曲线族如果残差出现周期性波动说明工况中存在周期性负载变化处理方向是工况归一化而不是改模型结构如果残差随预测提前量的增加而扩大说明特征提取环节噪声过大此时应回到时域变换链路检查包络谱的带通参数或 STFT 的频带选择是否偏离了实际共振区间。按这个顺序逐级排查通常能在不更换模型的前提下解决大部分预测偏差。本文还有配套的精品资源点击获取