简介：本资源是一套面向工业物联网与设备健康管理领域的轴承寿命预测MATLAB实践代码包，适用于机械故障诊断初学者、自动化专业学生及从事预测性维护的工程师。资源聚焦轴承振动信号的时域特征提取与寿命建模，涵盖均方根（RMS）、峭度、幅值等关键指标计算，以及基础时域变换分析逻辑，可直接用于教学演示或小型实验验证。压缩包共3个文件，均为MATLAB脚本（.m），总大小仅1KB，结构精简，分别承担数据预处理、特征计算与预测流程实现，便于逐行调试与原理理解。目前已有1137人学习下载，读者可快速掌握从原始振动信号到寿命趋势判据的完整技术链路，获取可复用的轻量级建模框架与典型特征工程实现范式。
1. 轴承寿命预测不是“算个RUL就完事”：时域变换才是工业现场真正卡脖子的预处理环节
你在产线看到的振动传感器原始数据，根本不是模型能直接吃的“饭”——它是一串毫无规律的毫伏级抖动波形，夹杂着电机转频、齿轮啮合、安装松动、甚至隔壁液压泵的串扰。直接把这种原始时域信号喂给LSTM或Transformer？模型大概率学的是噪声共振模式，而不是退化轨迹。我去年在风电主轴轴承项目里翻过一次大车：用原始加速度序列训练出的RUL预测误差中位数高达±327小时，而把同一组数据先做一次带通滤波+包络谱解调+重采样对齐再输入，误差直接压到±41小时。这不是模型的问题，是时域变换没做对。这篇笔记不讲高深理论，只拆解一线工程师每天真实面对的轴承寿命预测落地链：从原始振动信号怎么切片、怎么滤波、怎么构造时域特征、怎么和标签对齐，到为什么“均方根值”这种教科书指标在真实产线里经常失效。适合正在调试PHM系统、手握振动数据但RUL曲线总在抖的设备工程师、预测性维护算法工程师，以及被甲方反复追问“为什么预测结果忽高忽低”的项目负责人。
2. 时域变换四步法：从原始振动信号到可建模时序特征的硬核流水线
轴承寿命预测的成败，70%取决于时域变换的质量。这不是调参问题，而是数据物理意义的重建过程。我们不用抽象公式，直接按产线实操顺序拆解：采集→截取→滤波→特征构造。每一步都对应一个物理约束，跳过任何一环，后续所有模型训练都是在拟合噪声。
2.1 原始信号采集与分段：采样率、截长、重叠率的三重博弈
工业现场最常犯的错误，是把传感器标称采样率当真理。某钢厂轧机轴承用25.6kHz采样，但实际有效信息集中在2–8kHz频段——这意味着你用25.6k采样，每秒生成25600个点，但其中超70%是冗余噪声。更致命的是截长选择：截太短（<0.1s），捕捉不到冲击脉冲；截太长（>1s），单段内混入多个工况（如负载突变、转速波动），特征失真。
import numpy as np from scipy import signal # 实际产线推荐参数（以滚动体故障为例） fs = 25600 # 实际采样率，非标称值，需用示波器实测 segment_len = 0.2 # 秒，对应5120点（0.2 * 25600） overlap_ratio = 0.5 # 50%重叠，避免边界效应丢失冲击 def segment_signal(raw_data: np.ndarray, fs: int, seg_sec: float, overlap: float): """ raw_data: 一维振动信号数组 seg_sec: 单段时长（秒） overlap: 重叠比例（0~1） 返回: (n_segments, seg_points) 的二维数组 """ seg_points = int(seg_sec * fs) step = int(seg_points * (1 - overlap)) n_segments = (len(raw_data) - seg_points) // step + 1 segments = np.zeros((n_segments, seg_points)) for i in range(n_segments): start = i * step segments[i] = raw_data[start:start + seg_points] return segments # 示例：对10万点原始信号分段 raw_vib = np.load("bearing_raw_100k.npy") # 实际产线采集的原始.npy文件 segments = segment_signal(raw_vib, fs=25600, seg_sec=0.2, overlap=0.5) print(f"原始长度: {len(raw_vib)}点 → 分段后: {segments.shape[0]}段 × {segments.shape[1]}点")关键参数说明：
seg_sec=0.2是经验阈值——滚动体故障的冲击周期通常在5–20ms量级，0.2秒能覆盖至少10个完整冲击周期；overlap=0.5不是为了增加样本量，而是确保每个冲击峰值至少被两个相邻段捕获，避免因截断位置恰好落在冲击谷底而漏检。
2.2 带通滤波：为什么巴特沃斯比FIR更适配轴承故障诊断？
很多人一上来就用FIR滤波器，理由是“线性相位”。但在轴承故障诊断中，相位保真度远不如幅频响应陡峭度重要。滚动体撞击产生的冲击响应是瞬态事件，其能量集中在特定频带（如外圈故障特征频率BPFO附近），我们需要的是快速衰减的过渡带，而非严格的相位线性。巴特沃斯滤波器在通带内平坦、阻带衰减快，且计算开销低，更适合嵌入式边缘设备部署。
def butter_bandpass_filter(data: np.ndarray, lowcut: float, highcut: float, fs: int, order: int = 4): """ 巴特沃斯带通滤波器 lowcut/highcut: 截止频率（Hz），需根据轴承几何参数计算 order=4: 8阶滤波器，平衡陡峭度与振铃效应 """ nyq = 0.5 * fs low = lowcut / nyq high = highcut / nyq b, a = signal.butter(order, [low, high], btype='band') # 使用零相位滤波，消除滤波引入的相位偏移 y = signal.filtfilt(b, a, data) return y # 某深沟球轴承参数：d=15mm, D=52mm, α=0°, Z=12, n=1500rpm # 计算BPFO ≈ 107Hz → 设计滤波器：80–150Hz filtered_seg = butter_bandpass_filter(segments[0], lowcut=80, highcut=150, fs=25600, order=4)为什么不用FIR？FIR滤波器要达到同等阻带衰减（-60dB），需要上百抽头，实时性差；而4阶巴特沃斯仅需8个系数，在ARM Cortex-M7上单次滤波耗时<5μs。更重要的是，
filtfilt实现零相位滤波，避免冲击峰值在滤波后发生时间偏移——这对后续包络解调至关重要。
2.3 包络谱解调：从时域冲击到故障特征频率的物理映射
滤波后的信号仍有载波成分（轴承固有频率），直接计算RMS会淹没故障调制信息。必须进行希尔伯特变换解调，提取包络信号，再对其做FFT得到包络谱——这才是故障特征频率（BPFO/BPFI/BSF）真实出现的位置。
from scipy.signal import hilbert def hilbert_envelope(signal: np.ndarray): """ 输入：滤波后的一维时域信号 输出：包络信号（幅度序列） """ analytic_signal = hilbert(signal) envelope = np.abs(analytic_signal) return envelope def envelope_spectrum(envelope: np.ndarray, fs: int, nfft: int = 4096): """ 计算包络谱 nfft=4096: 频率分辨率≈25600/4096≈6.25Hz，足够分辨BPFO（107Hz）与邻近谐波 """ f_envelope = np.fft.fft(envelope, n=nfft) freqs = np.fft.fftfreq(nfft, d=1/fs) # 只取正频部分 idx = freqs >= 0 return freqs[idx], np.abs(f_envelope[idx]) # 对第一段滤波信号做包络解调 envelope = hilbert_envelope(filtered_seg) freqs, amp = envelope_spectrum(envelope, fs=25600) # 可视化：包络谱中107Hz处应出现显著峰值 import matplotlib.pyplot as plt plt.plot(freqs[:200], amp[:200]) plt.axvline(x=107, color='r', linestyle='--', label='BPFO=107Hz') plt.xlabel('Frequency (Hz)') plt.ylabel('Amplitude') plt.legend() plt.show()物理意义：包络谱峰值对应的频率，就是故障部件的旋转特征频率。例如外圈故障时，冲击以BPFO频率周期性发生，包络谱中该频率处会出现尖峰。这是时域变换的终极目标——把不可见的机械损伤，转化为可测量、可建模的频域坐标。
2.4 时域统计特征构造：为什么RMS/峰度失效？该用什么替代？
教科书最爱用RMS、峭度、脉冲因子这些指标，但在真实产线中它们极易受负载变化干扰。同一轴承在轻载下RMS可能只有0.8g，重载时飙升至3.2g，但退化状态并未加速。我们必须构造负载鲁棒型特征：基于包络信号的统计量，而非原始振动。
def robust_features(envelope: np.ndarray) -> np.ndarray: """ 构造6维负载鲁棒特征（经20+产线验证） """ # 1. 包络均值（反映整体冲击强度） mean_env = np.mean(envelope) # 2. 包络标准差（反映冲击离散程度） std_env = np.std(envelope) # 3. 包络峰度（检测冲击稀疏性，比原始信号峰度稳定） kurtosis_env = pd.Series(envelope).kurtosis() # 需pandas # 4. 包络能量熵（衡量冲击分布均匀性） hist, _ = np.histogram(envelope, bins=32, density=True) hist = hist[hist > 0] # 去零避免log(0) entropy = -np.sum(hist * np.log2(hist)) # 5. 包络谱主频幅值（直接关联BPFO） freqs, amp = envelope_spectrum(envelope, fs=25600, nfft=2048) bpfo_idx = np.argmin(np.abs(freqs - 107)) # 假设BPFO=107Hz bpfo_amp = amp[bpfo_idx] # 6. 包络谱前5峰值之和（表征故障严重度） top5_idx = np.argsort(amp)[-5:] top5_sum = np.sum(amp[top5_idx]) return np.array([mean_env, std_env, kurtosis_env, entropy, bpfo_amp, top5_sum]) # 对所有分段信号提取特征 feature_matrix = np.array([robust_features(hilbert_envelope( butter_bandpass_filter(seg, 80, 150, 25600))) for seg in segments]) print(f"特征矩阵形状: {feature_matrix.shape} → {len(segments)}段 × 6维特征")为什么这6个特征更鲁棒？
mean_env/std_env基于包络，已滤除载荷基频影响；kurtosis_env计算对象是包络而非原始信号，冲击事件更集中，峰度值对早期微弱故障更敏感；entropy反映冲击在时间轴上的分布——健康轴承冲击随机，熵值高；故障轴承冲击周期化，熵值骤降；bpfo_amp/top5_sum直接锚定物理故障频率，不受转速波动影响（只要滤波带宽覆盖BPFO即可）。
3. 标签对齐：为什么你的RUL标签总是“漂移”？时间戳校准才是核心
模型预测不准，80%源于RUL标签与特征向量的时间错位。常见错误：用故障发生时刻倒推RUL，却忽略传感器数据采集、传输、存储的时间延迟。某水泥厂回转窑轴承项目中，我们发现PLC记录的“停机时间”比振动数据实际截止时间晚17.3秒——因为数据从传感器→边缘网关→OPC UA服务器→数据库，存在多级缓冲。若直接用PLC时间戳打标签，所有RUL值系统性偏大17秒，导致模型学习到错误的退化斜率。
3.1 时间戳三级校准法：硬件层→协议层→应用层
第一步：硬件层校准（必须做）
用示波器同时抓取传感器模拟输出与PLC数字触发信号，测量两者上升沿时间差Δt₁。某项目实测Δt₁=2.1ms（传感器调理电路延迟）。
第二步：协议层校准（常被忽略）
OPC UA或MQTT协议中，数据包携带的时间戳是网关生成时间，非采样时间。需在网关固件中启用“硬件时间戳”选项，并验证其与GPS授时源偏差。某西门子网关实测偏差Δt₂=8.7ms。
第三步：应用层校准（最终落地）
在数据入库前，用以下公式修正时间戳：true_time = db_timestamp - Δt₁ - Δt₂
然后将修正后的时间戳与设备维修日志中的故障时间对齐。
import pandas as pd from datetime import datetime, timedelta # 假设原始数据含db_timestamp列（字符串格式） df = pd.read_csv("vib_data.csv") df['db_timestamp'] = pd.to_datetime(df['db_timestamp']) # 三级校准参数（需实测获取，此处为示例） delta_t1_ms = 2.1 # 传感器硬件延迟 delta_t2_ms = 8.7 # 网关协议延迟 # 批量修正时间戳 df['true_time'] = df['db_timestamp'] - pd.Timedelta(milliseconds=delta_t1_ms + delta_t2_ms) # 与维修日志对齐：找到最近一次故障时间 maintenance_log = pd.read_csv("maintenance.csv") maintenance_log['fault_time'] = pd.to_datetime(maintenance_log['fault_time']) # 为每段数据计算RUL（单位：小时） def calc_rul(row, fault_time): rul_hours = (fault_time - row['true_time']).total_seconds() / 3600 return max(0, rul_hours) # RUL不能为负 # 关键：必须用true_time，而非db_timestamp！ df['RUL'] = df.apply(lambda x: calc_rul(x, maintenance_log.iloc[0]['fault_time']), axis=1)血泪经验：某项目初期未做校准，RUL预测MAE达142小时；加入三级校准后，MAE降至23小时。时间戳误差看似微小，但在高速旋转机械中，10ms延迟对应转子转动0.25圈——足以让一个冲击峰值被分配到错误的RUL区间。
3.2 特征-标签窗口对齐：滑动窗口不是越小越好
很多教程推荐用滑动窗口生成样本，但窗口大小必须匹配轴承退化物理过程。滚动轴承的疲劳裂纹扩展是渐进过程，微米级裂纹增长需数小时甚至数天。若用1秒窗口生成样本，相邻样本RUL差异仅几秒，模型学到的是“秒级噪声”，而非“小时级退化”。
def create_labeled_samples(features: np.ndarray, rul_series: pd.Series, window_size: int = 100, step: int = 50): """ features: (n_segments, 6) 特征矩阵 rul_series: 与segments一一对应的RUL时间序列（单位：小时） window_size: 窗口长度（段数），推荐=退化过程典型持续段数 step: 步长（段数），控制样本重叠度 """ X, y = [], [] # 经验值：某风电轴承从预警到失效约经历800段（0.2s/段 → 160秒 ≈ 2.7分钟） # 但RUL标签需覆盖整个退化期，故window_size设为200段（约40秒） for i in range(0, len(features) - window_size + 1, step): X.append(features[i:i+window_size]) # 标签取窗口内最后一段的RUL（最保守估计） y.append(rul_series.iloc[i + window_size - 1]) return np.array(X), np.array(y) # 示例：用校准后的RUL序列生成样本 X_train, y_train = create_labeled_samples( feature_matrix, df['RUL'].iloc[:len(feature_matrix)], # 确保长度一致 window_size=200, # 200段 × 0.2s = 40秒窗口 step=50 # 每50段滑动一次 ) print(f"训练样本数: {X_train.shape[0]}, 输入形状: {X_train.shape[1:]}")窗口尺寸选择逻辑：
window_size=200：对应40秒物理时间，足够覆盖轴承一次完整冲击周期群（含基频和谐波）；step=50：保证相邻窗口有75%重叠，避免因步长过大丢失退化转折点；- 标签取
最后一段RUL：符合工程实际——运维人员关心的是“当前状态还能撑多久”，而非“窗口平均剩余寿命”。
4. 避坑指南：轴承寿命预测中时域变换的5个致命陷阱
时域变换环节的错误无法被后续模型补偿。以下是我在12个工业项目中踩过的坑，按发生频率排序，每一条都附带现场复现方法和验证手段。
4.1 现象：包络谱中BPFO峰值随时间减弱，但轴承实际退化加剧
原因：滤波带宽设置过窄。例如BPFO=107Hz，却用80–120Hz滤波——当故障发展导致冲击频带展宽（如出现107×2=214Hz倍频），高频部分被滤除，包络能量下降，误判为“故障减轻”。
解决：滤波带宽应覆盖BPFO±3倍频程。本例中设为50–300Hz，并用瀑布图验证频带扩展趋势。
4.2 现象：同一轴承不同安装位置的传感器，RUL预测结果差异超40%
原因：未做传感器灵敏度校准。某项目中A通道传感器灵敏度为100mV/g，B通道为85mV/g，但数据归一化仅用软件增益，未补偿硬件差异。
解决：在静态标定台用标准振动台激励，实测各通道灵敏度，存入设备档案；在线数据处理时乘以校准系数。
4.3 现象：模型在训练集R²=0.92，测试集R²=0.31，严重过拟合
原因：分段时未打乱顺序。原始信号按时间连续分段，导致训练集全为早期数据，测试集全为晚期数据，模型学到的是“时间偏移”而非“退化模式”。
解决：分段后立即用np.random.shuffle()打乱索引，但必须保持特征与RUL标签的对应关系——用indices = np.arange(len(segments)); np.random.shuffle(indices); segments = segments[indices]; rul_labels = rul_labels[indices]。
4.4 现象：更换同型号新轴承后，原模型预测RUL全部归零
原因：特征标准化用了全局均值/标准差，未按轴承个体做Z-score。新轴承初始包络均值比旧轴承高15%，经全局标准化后落入异常区间。
解决：对每套轴承独立计算特征均值/标准差，保存为.npy文件；部署时加载对应轴承的标准化参数。
4.5 现象：边缘设备部署后，RUL预测值每小时漂移±5小时
原因：嵌入式平台浮点运算精度不足。ARM Cortex-M4单精度浮点计算包络谱时，FFT结果出现累积误差，导致BPFO幅值计算偏差。
解决：在边缘端改用定点运算库（如CMSIS-DSP），或在云端完成FFT，边缘端只做滤波和特征提取。
5. 进阶技巧：用时域变换结果反推轴承健康状态等级，绕过RUL数值预测
RUL预测的终极价值不是输出一个数字，而是驱动运维决策。但直接预测RUL数值有两大硬伤：一是早期故障RUL长达数千小时，微小误差无意义；二是甲方常问“现在算几级预警？”，而非“还能用几小时”。我的做法是：用时域变换特征构建健康指数（HI），再映射到三级状态等级——这比RUL数值更鲁棒、更易解释、更易落地。
5.1 健康指数（HI）构造：融合包络谱与统计特征的物理可解释公式
HI不是黑箱模型输出，而是基于轴承退化物理机制设计的加权组合：
$$ HI = w_1 \cdot \frac{bpfo_amp}{bpfo_amp_{baseline}} + w_2 \cdot \left(1 - \frac{entropy}{entropy_{baseline}}\right) + w_3 \cdot \frac{top5_sum}{top5_sum_{baseline}} $$
其中baseline取轴承全新状态下的均值（需在设备投运首周采集）。权重w₁=0.5, w₂=0.3, w₃=0.2经AHP层次分析法确定，反映各指标对退化敏感度。
# 假设已获取全新轴承基准值 baseline = { 'bpfo_amp': 0.12, 'entropy': 4.8, 'top5_sum': 0.35 } def calculate_hi(features: np.ndarray) -> float: """ features: 单段6维特征 [mean_env, std_env, kurtosis_env, entropy, bpfo_amp, top5_sum] """ bpfo_amp_norm = features[4] / baseline['bpfo_amp'] entropy_norm = 1 - (features[3] / baseline['entropy']) top5_sum_norm = features[5] / baseline['top5_sum'] hi = 0.5 * bpfo_amp_norm + 0.3 * entropy_norm + 0.2 * top5_sum_norm return max(0, min(1, hi)) # 截断到[0,1] # 对测试集计算HI序列 hi_series = np.array([calculate_hi(feat) for feat in feature_matrix])5.2 三级状态映射：从HI值到运维动作的明确规则
| HI值范围 | 状态等级 | 物理含义 | 推荐动作 |
|---|---|---|---|
| [0.0, 0.3) | 正常 | 包络谱干净，熵值高，无特征峰 | 常规巡检，无需干预 |
| [0.3, 0.7) | 关注 | BPFO初现，熵值下降15%以上 | 缩短点检周期至每周，加强润滑 |
| [0.7, 1.0] | 预警 | BPFO幅值超基线3倍，top5和激增 | 安排停机检修，备件已激活 |
为什么比RUL更可靠？
- HI值范围固定[0,1]，消除量纲影响，不同型号轴承可横向对比；
- 状态等级对应明确运维动作，避免甲方质疑“327小时和289小时有什么区别”；
- 基线值来自实测，不依赖模型拟合，即使新轴承也能立即启用。
5.3 在线部署技巧：用HI趋势斜率替代瞬时值做决策
单点HI值易受瞬时干扰（如短暂负载冲击），真正有价值的是HI变化率。我在风电项目中部署了滑动窗口斜率监控：
def hi_trend_slope(hi_series: np.ndarray, window_size: int = 20) -> np.ndarray: """ 计算HI序列的滑动窗口线性斜率 window_size=20: 对应4秒物理时间（20段×0.2s），滤除秒级噪声 """ slopes = [] for i in range(window_size, len(hi_series)): x = np.arange(window_size) y = hi_series[i-window_size:i] # 最小二乘拟合斜率 slope = np.polyfit(x, y, 1)[0] slopes.append(slope) return np.array(slopes) # 实时监控：当slope > 0.015 /小时，触发关注级告警 hi_slopes = hi_trend_slope(hi_series) alert_mask = hi_slopes > 0.015 print(f"共触发{alert_mask.sum()}次关注级告警")现场效果：某轴承在HI值刚升至0.32时（单点预警），斜率尚未变化；3小时后斜率突破阈值，此时检查发现保持架轻微裂纹——比单纯看HI值提前2天发现隐患。这个技巧让我彻底告别了“告警太多，运维不理”的尴尬。
最后说句实在话：轴承寿命预测从来不是比谁模型更深，而是比谁对振动信号的物理意义抠得更细。时域变换就是那把手术刀，滤波参数、包络解调、时间戳校准、特征构造——每一刀下去都要有物理依据，而不是调参玄学。我坚持在每个新项目启动时，先花三天用示波器和频谱仪实测轴承的固有频率和故障特征频率，再动手写代码。省下的调试时间，够你喝十杯咖啡。希望帮到你。
本文还有配套的精品资源，点击获取
