轴承故障诊断：振动信号分析与Fisher判别实战

1. 轴承故障诊断技术概述

在工业设备维护领域，轴承故障诊断一直是个既基础又关键的课题。作为一名长期从事设备状态监测的工程师，我见证过太多因为轴承故障未及时发现而导致的生产事故。记得去年某化工厂的主风机轴承突然失效，导致整条生产线停机48小时，直接经济损失超过200万元。这类惨痛教训让我深刻意识到：轴承虽小，却是旋转机械的"关节"，其健康状态直接关系到整个生产系统的安危。

传统的人工巡检方式存在明显局限：

• 依赖经验：老师傅的"听音辨病"难以量化
• 时效性差：定期检查可能错过早期故障征兆
• 主观性强：不同人员判断标准不一致

随着传感器技术和数据分析方法的发展，基于振动信号的智能诊断已成为行业主流。其中，时域特征提取+Fisher判别分析的组合方案，因其实现简单、效果稳定，特别适合中小型企业的技术改造。这个方案的核心优势在于：

1. 硬件成本低：只需常规振动传感器
2. 计算效率高：适合在线监测
3. 解释性强：特征具有明确的物理意义

2. 技术实现原理详解

2.1 振动信号采集规范

优质的原始数据是分析的基础。在实际项目中，我们采用以下采集标准：

传感器安装要点：

• 位置选择：优先测量径向振动，传感器应安装在轴承座刚性最好的部位
• 方向确定：主振动方向（通常为垂直方向）必须与传感器轴线对齐
• 耦合方式：磁座适用于临时检测，永久监测应采用螺纹固定

采样参数设置：

matlab复制% 典型参数设置示例
fs = 12000;  % 采样频率(Hz)，根据轴承特征频率确定
N = 8192;    % 采样点数，保证包含至少5个完整的轴承旋转周期

关键经验：采样频率应至少为轴承外圈故障特征频率的10倍，对于常见工业轴承，12kHz通常足够。

2.2 时域特征工程实践

时域特征直接从振动波形中提取，计算效率高且物理意义明确。以下是我们在多个项目中验证有效的特征集：

基本统计特征：

• 峰值（Peak）：反映瞬时冲击强度
• 均方根值（RMS）：表征整体振动能量
• 峭度（Kurtosis）：对冲击敏感，早期故障指示性好

无量纲特征：

matlab复制function [CF, IF, SF] = dimensionlessFeatures(x)
    peak = max(abs(x));
    rms = sqrt(mean(x.^2));
    CF = peak/rms;              % 峰值因子
    IF = peak/mean(abs(x));     % 脉冲因子
    SF = rms/mean(abs(x));      % 波形因子
end

实际应用心得：

1. 不同故障类型对特征的敏感度不同：
   • 剥落故障：重点关注峰值和脉冲因子
   • 磨损故障：观察RMS值的趋势变化
2. 特征值需做归一化处理，消除量纲影响：

   matlab复制X_normalized = (X - min(X)) / (max(X) - min(X));
   

3. 建议构建包含8-12个特征的特征池，后续通过Fisher判别进行筛选

2.3 Fisher判别分析实战

Fisher判别的核心思想是最大化类间离散度与最小化类内离散度的比值。其实现步骤包括：

1. 计算散度矩阵：

matlab复制% 输入参数说明：
% X: n×d特征矩阵（n样本数，d特征维数）
% y: n×1标签向量

classes = unique(y);
Sw = zeros(d); Sb = zeros(d);  % 初始化类内、类间散度矩阵
mu = mean(X);                  % 全局均值

for c = classes'
    Xc = X(y==c,:);
    n_c = size(Xc,1);
    mu_c = mean(Xc);
    Sw = Sw + (Xc - mu_c)'*(Xc - mu_c);
    Sb = Sb + n_c*(mu_c - mu)'*(mu_c - mu);
end

2. 求解投影向量：

matlab复制[V, D] = eig(Sb, Sw);          % 广义特征值分解
[~, idx] = sort(diag(D), 'descend');
W = V(:, idx(1:k));            % 取前k个最大特征值对应向量

工程应用技巧：

• 维度选择：通常保留2-3个判别分量即可达到90%以上的分类准确率
• 可视化分析：将数据投影到二维平面观察聚类效果
• 正则化处理：当Sw接近奇异时，加入小扰动项Sw = Sw + λ*eye(d)

3. MATLAB实现全流程

3.1 数据准备与预处理

我们使用美国凯斯西储大学（CWRU）轴承数据集进行演示，这是行业公认的基准数据：

数据加载与分割：

matlab复制load('bearing_data.mat');  % 假设已整理好的数据文件
[trainData, testData] = splitData(data, 0.7);  % 70%训练，30%测试

% 数据标准化（重要！）
[trainData.X, mu, sigma] = zscore(trainData.X);
testData.X = (testData.X - mu) ./ sigma;

标签处理技巧：

• 将健康状态编码为数字标签：正常=1，内圈故障=2，外圈故障=3，滚动体故障=4
• 对于不平衡数据集，可采用SMOTE过采样或调整类别权重

3.2 特征提取模块实现

完整特征提取函数示例：

matlab复制function features = extractTimeFeatures(signal)
    features = zeros(1, 10);  % 预分配空间
    
    % 时域统计特征
    features(1) = max(abs(signal));                  % 峰值
    features(2) = rms(signal);                      % RMS
    features(3) = kurtosis(signal);                 % 峭度
    features(4) = skewness(signal);                 % 偏度
    
    % 无量纲指标
    peak = features(1);
    rms_val = features(2);
    mean_abs = mean(abs(signal));
    
    features(5) = peak / rms_val;                   % 峰值因子
    features(6) = peak / mean_abs;                  % 脉冲因子
    features(7) = rms_val / mean_abs;               % 波形因子
    features(8) = peak / (mean(signal.^2)^0.5);     % 裕度因子
    
    % 其他特征
    features(9) = entropy(signal);                  % 近似熵
    features(10) = sum(abs(diff(signal)));          % 绝对平均斜率
end

3.3 分类模型构建与评估

完整训练流程：

matlab复制% 特征提取
trainFeatures = cell2mat(arrayfun(@(i) extractTimeFeatures(trainData.signal{i}),...
    1:numel(trainData.signal), 'UniformOutput', false));

% Fisher判别分析
[W, lambda] = fisherDiscriminant(trainFeatures, trainData.label);

% 降维投影
trainScore = trainFeatures * W(:,1:2);  % 取前两个判别向量

% SVM分类器训练
mdl = fitcecoc(trainScore, trainData.label,...
    'Learners', templateSVM('KernelFunction','gaussian'));

模型评估方法：

matlab复制% 测试集处理
testFeatures = cell2mat(arrayfun(@(i) extractTimeFeatures(testData.signal{i}),...
    1:numel(testData.signal), 'UniformOutput', false));
testScore = testFeatures * W(:,1:2);

% 预测与评估
pred = predict(mdl, testScore);
accuracy = sum(pred == testData.label) / numel(testData.label);
confmat = confusionmat(testData.label, pred);

4. 工程应用中的挑战与对策

4.1 噪声干扰处理方案

现场环境噪声是影响诊断精度的主要因素。我们总结的降噪策略包括：

硬件层面：

• 选用IEPE型加速度传感器，自带信号调理
• 安装抗电磁干扰的屏蔽电缆
• 在传感器与安装面之间使用绝缘垫片

软件处理：

matlab复制% 小波阈值降噪示例
[cleanSignal, ~, ~] = wden(signal, 'rigrsure', 's', 'sln', 5, 'db4');

频带选择技巧：

1. 计算原始信号的包络谱
2. 识别轴承特征频率所在的频带
3. 设计带通滤波器提取目标频段

4.2 特征选择优化方法

当特征维度较高时，可采用以下方法优化：

序列前向选择(SFS)流程：

1. 初始化空特征集F=∅
2. 每次从未选特征中加入一个使分类准确率提升最大的特征
3. 直到准确率不再显著提高

基于Fisher得分的筛选：

matlab复制% 计算每个特征的Fisher得分
score = (mean(feature1) - mean(feature2))^2 / (var(feature1) + var(feature2));

4.3 在线监测系统设计要点

将算法部署到实际监测系统中需考虑：

实时性保障：

• 采用滑动窗口处理：窗口长度1024点，重叠率50%
• 特征计算优化：预先计算好统计量，避免重复运算

状态判定逻辑：

mermaid复制graph TD
    A[原始振动信号] --> B[特征提取]
    B --> C{特征值超过阈值?}
    C -->|是| D[触发报警]
    C -->|否| E[继续监测]
    D --> F[保存故障数据]

系统集成建议：

• 数据采集：NI CompactDAQ或类似设备
• 通信协议：采用OPC UA实现与PLC的数据交互
• 界面显示：用趋势图展示关键特征值变化

5. 效果验证与案例分享

在某水泥厂风机轴承监测项目中，我们实施了这套方案，取得了显著成效：

实施效果对比：

典型故障识别案例：

1. 早期内圈剥落：

   • 第3天：峭度值开始缓慢上升（从3.2→4.5）
   • 第7天：脉冲因子超过阈值1.5倍
   • 第10天：肉眼可见磨损痕迹

2. 润滑不良诊断：

   • 特征表现：RMS值渐进升高，同时峰值因子下降
   • 处理建议：检查润滑系统，补充专用油脂

这套方法在近两年的实践中不断完善，现在已经形成标准化实施流程。对于想入门设备智能运维的工程师，我的建议是：

1. 先从CWRU标准数据集练手
2. 重点理解每个特征的物理意义
3. 逐步尝试添加自己的新特征
4. 最后再考虑引入深度学习等复杂方法