一维信号频域特征提取在轴承故障诊断与趋势预测中的应用

轴承故障诊断和趋势预测是工业设备健康管理的核心内容，频域特征提取在这方面发挥着至关重要的作用。

1. 频域分析的基本原理

轴承振动信号的频域分析基于傅里叶变换，将时域信号转换为频域表示，从而揭示信号的频率组成特征。轴承故障会产生特定的频率成分，这些成分可以作为故障诊断的依据。

1.1 轴承故障特征频率计算

不同类型的轴承故障会产生不同的特征频率：

• 内圈故障频率 (BPFI): BPFI = (n/2) * f_r * (1 + (d/D) * cosφ)
• 外圈故障频率 (BPFO): BPFO = (n/2) * f_r * (1 - (d/D) * cosφ)
• 滚动体故障频率 (BSF): BSF = (D/(2*d)) * f_r * (1 - (d/D)^2 * cos^2φ)
• 保持架故障频率 (FTF): FTF = (1/2) * f_r * (1 - (d/D) * cosφ)

其中：

• n = 滚动体数量
• f_r = 轴旋转频率 (Hz)
• d = 滚动体直径
• D = 轴承节径
• φ = 接触角

2. MATLAB实现频域特征提取

2.1 基本频域特征提取函数

function features = extractFrequencyFeatures(signal, fs)
    % 提取信号的频域特征
    % 输入:
    %   signal - 输入信号
    %   fs - 采样频率
    % 输出:
    %   features - 包含各种频域特征的结构体
    
    % 计算FFT
    N = length(signal);
    fft_data = fft(signal);
    fft_magnitude = abs(fft_data(1:floor(N/2)+1);
    fft_magnitude = fft_magnitude / max(fft_magnitude); % 归一化
    freq_axis = (0:floor(N/2)) * fs / N;
    
    % 计算功率谱密度 (PSD)
    [psd, freq_psd] = pwelch(signal, hanning(256), 128, 256, fs);
    
    % 初始化特征结构体
    features = struct();
    
    % 1. 重心频率 (Spectral Centroid)
    features.centroid = sum(freq_axis .* fft_magnitude) / sum(fft_magnitude);
    
    % 2. 均方频率 (Mean Square Frequency)
    features.mean_square_freq = sum(freq_axis.^2 .* fft_magnitude) / sum(fft_magnitude);
    
    % 3. 频率方差 (Frequency Variance)
    features.freq_variance = sum((freq_axis - features.centroid).^2 .* fft_magnitude) / sum(fft_magnitude);
    
    % 4. 频率标准差 (Frequency Standard Deviation)
    features.freq_std = sqrt(features.freq_variance);
    
    % 5. 频谱偏度 (Spectral Skewness)
    features.skewness = sum((freq_axis - features.centroid).^3 .* fft_magnitude) / ...
                       (sum(fft_magnitude) * features.freq_std^3);
    
    % 6. 频谱峭度 (Spectral Kurtosis)
    features.kurtosis = sum((freq_axis - features.centroid).^4 .* fft_magnitude) / ...
                       (sum(fft_magnitude) * features.freq_variance^2);
    
    % 7. 频带能量特征
    bands = [0, 0.1*fs/2; 0.1*fs/2, 0.2*fs/2; 0.2*fs/2, 0.3*fs/2; 
             0.3*fs/2, 0.4*fs/2; 0.4*fs/2, 0.5*fs/2];
    
    for i = 1:size(bands, 1)
        band_mask = freq_axis >= bands(i,1) & freq_axis <= bands(i,2);
        features.(sprintf('band_energy_%d', i)) = sum(fft_magnitude(band_mask));
    end
    
    % 8. 峰值频率及其幅度
    [pks, locs] = findpeaks(fft_magnitude, 'SortStr', 'descend', 'NPeaks', 5);
    if length(pks) >= 5
        for i = 1:5
            features.(sprintf('peak_freq_%d', i)) = freq_axis(locs(i));
            features.(sprintf('peak_amp_%d', i)) = pks(i);
        end
    end
    
    % 9. 频谱熵 (Spectral Entropy)
    pdf = fft_magnitude / sum(fft_magnitude); % 概率密度函数
    features.spectral_entropy = -sum(pdf .* log2(pdf + eps));
    
    % 10. 包络谱特征 (Envelope Spectrum Features)
    [env_features, env_freq] = extractEnvelopeFeatures(signal, fs);
    features.envelope_peak_ratio = env_features.peak_ratio;
    features.envelope_centroid = env_features.centroid;
end

function [env_features, freq_axis] = extractEnvelopeFeatures(signal, fs)
    % 提取包络谱特征
    % 希尔伯特变换获取包络
    analytic_signal = hilbert(signal);
    envelope = abs(analytic_signal);
    
    % 计算包络谱
    N = length(envelope);
    env_spectrum = abs(fft(envelope - mean(envelope)));
    env_spectrum = env_spectrum(1:floor(N/2)+1);
    freq_axis = (0:floor(N/2)) * fs / N;
    
    % 归一化
    env_spectrum = env_spectrum / max(env_spectrum);
    
    % 提取包络谱特征
    env_features.centroid = sum(freq_axis .* env_spectrum) / sum(env_spectrum);
    
    % 计算峰值比例 (故障特征)
    [pks, locs] = findpeaks(env_spectrum, 'MinPeakHeight', 0.1);
    if length(pks) > 1
        sorted_pks = sort(pks, 'descend');
        env_features.peak_ratio = sorted_pks(2) / sorted_pks(1);
    else
        env_features.peak_ratio = 0;
    end
end

2.2 轴承故障诊断系统

function bearingFaultDiagnosis()
    % 轴承故障诊断主函数
    
    % 加载数据（这里使用模拟数据，实际应用中应替换为真实数据）
    [healthy_data, faulty_data, fs] = loadBearingData();
    
    % 提取健康轴承特征
    healthy_features = [];
    for i = 1:size(healthy_data, 2)
        features = extractFrequencyFeatures(healthy_data(:, i), fs);
        healthy_features = [healthy_features; struct2array(features)];
    end
    
    % 提取故障轴承特征
    faulty_features = [];
    for i = 1:size(faulty_data, 2)
        features = extractFrequencyFeatures(faulty_data(:, i), fs);
        faulty_features = [faulty_features; struct2array(features)];
    end
    
    % 创建标签
    labels = [zeros(size(healthy_features, 1), 1); ones(size(faulty_features, 1), 1)];
    all_features = [healthy_features; faulty_features];
    
    % 特征标准化
    all_features = zscore(all_features);
    
    % 训练分类器（使用SVM）
    rng(1); % 设置随机种子以确保可重复性
    cv = cvpartition(labels, 'HoldOut', 0.3);
    idx_train = training(cv);
    idx_test = test(cv);
    
    % 使用PCA进行特征降维
    [coeff, score, ~, ~, explained] = pca(all_features(idx_train, :));
    num_components = find(cumsum(explained) >= 95, 1); % 保留95%方差
    
    train_features = score(:, 1:num_components);
    test_features = (all_features(idx_test, :) - mean(all_features(idx_train, :))) * coeff(:, 1:num_components);
    
    % 训练SVM分类器
    svm_model = fitcsvm(train_features, labels(idx_train), ...
        'KernelFunction', 'rbf', 'Standardize', true, ...
        'KernelScale', 'auto', 'BoxConstraint', 1);
    
    % 预测
    predictions = predict(svm_model, test_features);
    
    % 评估性能
    accuracy = sum(predictions == labels(idx_test)) / length(labels(idx_test));
    confusion_mat = confusionmat(labels(idx_test), predictions);
    
    fprintf('诊断准确率: %.2f%%\n', accuracy * 100);
    disp('混淆矩阵:');
    disp(confusion_mat);
    
    % 可视化特征空间
    visualizeFeatures(train_features, labels(idx_train), svm_model);
end

function [healthy_data, faulty_data, fs] = loadBearingData()
    % 加载轴承数据（这里使用模拟数据）
    % 在实际应用中，应替换为真实数据加载代码
    fs = 12000; % 采样频率 12kHz
    t = 0:1/fs:1-1/fs;
    
    % 生成健康轴承数据（主要是随机噪声）
    healthy_data = zeros(length(t), 10);
    for i = 1:10
        healthy_data(:, i) = 0.1 * randn(size(t)) + 0.05 * sin(2*pi*50*t);
    end
    
    % 生成故障轴承数据（包含故障特征频率）
    faulty_data = zeros(length(t), 10);
    fault_freq = 100; % 假设故障特征频率为100Hz
    
    for i = 1:10
        % 基础噪声
        base_signal = 0.1 * randn(size(t)) + 0.05 * sin(2*pi*50*t);
        
        % 添加故障冲击（周期性冲击）
        impulse_train = zeros(size(t));
        impulse_period = round(fs / fault_freq);
        impulse_train(1:impulse_period:end) = 1;
        
        % 冲击响应（衰减正弦波）
        impulse_response = exp(-500*t) .* sin(2*pi*2000*t);
        fault_signal = conv(impulse_train, impulse_response, 'same');
        
        faulty_data(:, i) = base_signal + 0.3 * fault_signal / max(abs(fault_signal));
    end
end

function visualizeFeatures(features, labels, svm_model)
    % 可视化特征空间和分类边界
    
    figure;
    
    % 如果特征维度大于2，使用前两个主成分
    if size(features, 2) > 2
        [~, score] = pca(features);
        features_2d = score(:, 1:2);
    else
        features_2d = features;
    end
    
    % 创建网格用于绘制决策边界
    h = 0.02;
    x1_min = min(features_2d(:, 1)) - 1;
    x1_max = max(features_2d(:, 1)) + 1;
    x2_min = min(features_2d(:, 2)) - 1;
    x2_max = max(features_2d(:, 2)) + 1;
    [xx1, xx2] = meshgrid(x1_min:h:x1_max, x2_min:h:x2_max);
    
    % 预测网格点的类别
    Z = predict(svm_model, [xx1(:), xx2(:)]);
    Z = reshape(Z, size(xx1));
    
    % 绘制决策区域
    contourf(xx1, xx2, Z, 'AlphaData', 0.3);
    colormap([0.8 0.9 0.9; 0.9 0.8 0.9]);
    hold on;
    
    % 绘制数据点
    gscatter(features_2d(:, 1), features_2d(:, 2), labels, 'br', 'xo', 10);
    
    xlabel('主成分 1');
    ylabel('主成分 2');
    title('轴承状态分类结果');
    legend('健康', '故障', 'Location', 'best');
    grid on;
    hold off;
end

2.3 趋势预测与健康指标构建

function healthIndicatorTrendPrediction()
    % 轴承健康指标构建与趋势预测
    
    % 生成模拟时间序列数据（实际应用中应使用真实数据）
    [time_series, time] = generateTimeSeriesData();
    
    % 提取每个时间点的频域特征
    features_over_time = [];
    for i = 1:size(time_series, 2)
        feature_vector = extractFrequencyFeatures(time_series(:, i), 12000);
        features_over_time = [features_over_time; struct2array(feature_vector)];
    end
    
    % 构建健康指标（使用多个特征的加权组合）
    % 这里使用第一主成分作为健康指标
    [coeff, score, ~, ~, explained] = pca(features_over_time);
    health_indicator = score(:, 1);
    
    % 健康指标归一化（0-1范围，1表示健康，0表示故障）
    health_indicator = 1 - (health_indicator - min(health_indicator)) / ...
        (max(health_indicator) - min(health_indicator));
    
    % 趋势预测（使用ARIMA模型）
    num_train = floor(0.7 * length(health_indicator));
    train_data = health_indicator(1:num_train);
    test_data = health_indicator(num_train+1:end);
    
    % 训练ARIMA模型
    arima_model = arima(5, 1, 2); % ARIMA(5,1,2)
    estimated_model = estimate(arima_model, train_data);
    
    % 预测
    [yF, yMSE] = forecast(estimated_model, length(test_data), 'Y0', train_data);
    lower = yF - 1.96 * sqrt(yMSE);
    upper = yF + 1.96 * sqrt(yMSE);
    
    % 计算预测误差
    mse = mean((test_data - yF).^2);
    rmse = sqrt(mse);
    mae = mean(abs(test_data - yF));
    
    fprintf('趋势预测性能:\n');
    fprintf('RMSE: %.4f\n', rmse);
    fprintf('MAE: %.4f\n', mae);
    
    % 可视化结果
    plotTrendResults(time, health_indicator, train_data, test_data, yF, lower, upper);
    
    % 预测剩余使用寿命 (RUL)
    predictRUL(health_indicator, time);
end

function [time_series, time] = generateTimeSeriesData()
    % 生成模拟时间序列数据（轴承退化过程）
    fs = 12000;
    measurement_interval = 3600; % 每小时的测量间隔（秒）
    total_time = 30 * 24 * 3600; % 30天的总时间（秒）
    
    time = 0:measurement_interval:total_time;
    time_series = zeros(floor(fs * 0.5), length(time)); % 每次测量0.5秒数据
    
    % 生成退化过程数据
    for i = 1:length(time)
        % 基础信号（随着时间衰减）
        t_signal = 0:1/fs:0.5-1/fs;
        base_signal = 0.1 * randn(size(t_signal)) + 0.05 * sin(2*pi*50*t_signal);
        
        % 故障发展（随时间增强）
        degradation_level = 1 - exp(-0.0001 * time(i)); % 退化水平
        fault_freq = 100 + 50 * degradation_level; % 故障频率随时间变化
        
        % 故障冲击
        impulse_train = zeros(size(t_signal));
        impulse_period = round(fs / fault_freq);
        impulse_train(1:impulse_period:end) = 1;
        
        impulse_response = exp(-500*t_signal) .* sin(2*pi*2000*t_signal);
        fault_signal = conv(impulse_train, impulse_response, 'same');
        
        % 组合信号
        time_series(:, i) = base_signal + degradation_level * 0.5 * fault_signal / max(abs(fault_signal));
    end
end

function plotTrendResults(time, health_indicator, train_data, test_data, yF, lower, upper)
    % 绘制趋势预测结果
    
    figure;
    
    % 转换时间为天
    time_days = time / (24 * 3600);
    num_train = length(train_data);
    
    % 绘制实际健康指标
    plot(time_days, health_indicator, 'b-', 'LineWidth', 1.5, 'DisplayName', '实际健康指标');
    hold on;
    
    % 绘制预测结果
    plot(time_days(num_train+1:end), yF, 'r--', 'LineWidth', 1.5, 'DisplayName', '预测值');
    fill([time_days(num_train+1:end), fliplr(time_days(num_train+1:end))], ...
         [lower', fliplr(upper')], 'r', 'FaceAlpha', 0.2, 'EdgeColor', 'none', ...
         'DisplayName', '95%置信区间');
    
    % 添加分界线
    plot([time_days(num_train), time_days(num_train)], [0, 1], 'k--', 'LineWidth', 1, 'DisplayName', '训练/测试分界');
    
    xlabel('时间 (天)');
    ylabel('健康指标');
    title('轴承健康趋势预测');
    legend('show');
    grid on;
    hold off;
end

function predictRUL(health_indicator, time)
    % 预测剩余使用寿命 (RUL)
    
    % 设置故障阈值
    failure_threshold = 0.2;
    
    % 找到健康指标首次低于阈值的时间点
    failure_idx = find(health_indicator < failure_threshold, 1);
    
    if isempty(failure_idx)
        fprintf('在当前监测期内未检测到故障发展至阈值\n');
        return;
    end
    
    % 使用线性回归预测RUL
    last_n_points = 20; % 使用最后20个点进行预测
    if length(health_indicator) < last_n_points
        last_n_points = length(health_indicator);
    end
    
    % 选择用于拟合的数据点
    fit_indices = (length(health_indicator)-last_n_points+1):length(health_indicator);
    x = time(fit_indices)';
    y = health_indicator(fit_indices)';
    
    % 线性拟合
    p = polyfit(x, y, 1);
    y_fit = polyval(p, x);
    
    % 求解达到故障阈值的时间
    % p(1)*t + p(2) = failure_threshold
    failure_time = (failure_threshold - p(2)) / p(1);
    
    % 计算RUL
    current_time = time(end);
    rul = failure_time - current_time;
    
    fprintf('剩余使用寿命预测:\n');
    fprintf('当前时间: %.2f 天\n', current_time/(24*3600));
    fprintf('预测故障时间: %.2f 天\n', failure_time/(24*3600));
    fprintf('剩余使用寿命: %.2f 天\n', rul/(24*3600));
    
    % 绘制RUL预测图
    figure;
    plot(time/(24*3600), health_indicator, 'b-', 'LineWidth', 1.5, 'DisplayName', '健康指标');
    hold on;
    plot(x/(24*3600), y_fit, 'r--', 'LineWidth', 1.5, 'DisplayName', '趋势拟合');
    plot([min(x), failure_time]/(24*3600), [failure_threshold, failure_threshold], 'k:', 'LineWidth', 1, 'DisplayName', '故障阈值');
    plot([failure_time, failure_time]/(24*3600), [0, failure_threshold], 'k:', 'LineWidth', 1);
    
    xlabel('时间 (天)');
    ylabel('健康指标');
    title(sprintf('剩余使用寿命预测: %.2f 天', rul/(24*3600)));
    legend('show');
    grid on;
    hold off;
end

参考代码 一维信号频域特征提取可用于轴承故障诊断和趋势预测 www.3dddown.com/cnb/50931.html

3. 实际应用建议

3.1 数据预处理注意事项

1. 降噪处理：在实际工业环境中，振动信号通常包含大量噪声，建议使用小波降噪或自适应滤波等方法进行预处理。

2. 转速变化处理：如果设备转速变化较大，应考虑使用阶比分析(order analysis)代替传统的FFT分析。

3. 数据标准化：不同传感器、不同时间采集的数据可能存在尺度差异，应进行适当的标准化处理。

3.2 特征选择与优化

1. 相关性分析：计算各个特征与故障状态的相关性，选择相关性高的特征。

2. 递归特征消除：使用递归特征消除(RFE)方法选择最优特征子集。

3. 领域知识引导：结合轴承故障机理知识，优先选择与故障物理过程直接相关的特征。

3.3 系统部署考虑

1. 实时性要求：根据实际应用的实时性要求，优化算法复杂度。

2. 自适应阈值：设置自适应故障阈值，避免因工况变化导致的误报。

3. 模型更新机制：建立模型定期更新机制，适应设备老化等长期变化。

4. 总结

频域特征提取在轴承故障诊断和趋势预测中具有重要作用。通过提取频谱重心、频带能量、包络谱特征等多种频域特征，并结合机器学习方法，可以构建高效的故障诊断和预测系统。