MATLAB环境基于健康指标（Health indicator）的滚动轴承故障诊断

轴承的剩余使用寿命RUL预测过程一般包括以下三个步骤：（1）数据采集，（2）健康指标HI构建，（3）RUL预测。在预测过程中，RUL并不能直接依靠观测得到，其主要原因是轴承的工作环境复杂，噪声、温度、湿度等其他因素往往对振动信号产生一定的干扰。同时由于轴承的退化过程通常是非线性的，其退化模式难以通过人为方式捕捉，因此HI的构建显得非常重要。HI被用来作为反映退化过程的数据指标，和RUL存在一一对应的关系。性能优越的健康指标在表征轴承退化趋势的过程中发挥着重要的作用，可以降低神经网络学习故障特征的困难度和复杂度，缩短模型训练时间，减小网络误差，提高RUL预测精度。健康指标作为评估轴承健康状态以及运行状况的重要参数，应能够表征轴承的退化过程或退化程度，这就要求从机械设备的原始监测数据中提取轴承退化特征、挖掘时序信息，构建能够充分、全面、准确地表征轴承健康状态和运行状况的健康指标曲线，从而有效地评估轴承健康状态和预测轴承的RUL。

HI的构建对剩余使用寿命的预测发挥着至关重要的作用，尽管随着机器学习技术和深度学习技术的发展和应用，现有的健康指标构建方法可以有效地构建出反映轴承退化过程的HI，在一定程度上提高了健康指标的单调性和相关性，为RUL预测奠定了良好的基础，但仍有以下几个问题需要重点关注：

（1）轴承特征参数的数量较多，在前期进行特征筛选时多依赖于专家经验和先验知识，缺乏一定的评价标准，为此，需要构建一种较为完善和准确的参数筛选和健康指标评价标准，实现原始特征参数的筛选和健康指标性能的评估。

（2）基于机器学习和深度学习技术构建HI的方法，主要是通过特征提取和特征融合，建立健康指标和剩余使用寿命之间的非线性关系，然而在特征提取和融合时存在着信息丢失的现象，导致HI不能充分反映轴承的退化过程，影响预测结果。

（3）由于轴承工作环境复杂，工作状况易受温度、湿度、噪声、负载和其他因素的干扰，在不同运行阶段表现出不同的退化特征，基于有监督的方法构建健康指标时，需合理划分健康阶段，为网络输入数据设置科学合理的标签。

（4）轴承数据量较大，模型训练时容易出现训练速度较低、过拟合的现象，如何提高网络运算速度、避免过拟合现象的发生是亟需解决的问题。

鉴于此，提出一种基于健康指标（Health indicator）的滚动轴承故障诊断方法，运行环境为MATLAB R2018A，轴承特征参数及HI构建方法如下：

 % features extraction
    shaft_speed = 1 ; % in the cycle domain the shaft speed is 1.
    [ftf, bsf, bpfo, bpfi] = calc_bearing_tones(shaft_speed, bearing_specifications.num_balls, ...
    bearing_specifications.ball_diameter, bearing_specifications.pitch_diameter, bearing_specifications.bearing_contact_angle) ;
    bearing_tones = [ftf, bsf, bpfo, bpfi] ;
    
    order_sensativity = 0.02 ; % sensitivity of the order value for searching the bearing tone pick
    first_ftf_value = extract_specific_bearing_tone(dorder, sig_env_order, ftf, order_sensativity) ;
    first_bsf_value = extract_specific_bearing_tone(dorder, sig_env_order, bsf, order_sensativity) ;
    first_bpfo_value = extract_specific_bearing_tone(dorder, sig_env_order, bpfo, order_sensativity) ;
    first_bpfi_value = extract_specific_bearing_tone(dorder, sig_env_order, bpfi, order_sensativity) ;

    % HI calculation
    hi_ftf = mean(abs(first_ftf_value - healthy_ftf_average)./healthy_ftf_std) ;
    hi_ftf_faulty_vctr(sig_num) = hi_ftf ;
    hi_bsf = mean(abs(first_bsf_value - healthy_bsf_average)./healthy_bsf_std) ;
    hi_bsf_faulty_vctr(sig_num) = hi_bsf ;
    hi_bpfo = mean(abs(first_bpfo_value - healthy_bpfo_average)./healthy_bpfo_std) ;
    hi_bpfo_faulty_vctr(sig_num) = hi_bpfo ;
    hi_bpfi = mean(abs(first_bpfi_value - healthy_bpfi_average)./healthy_bpfi_std) ;
    hi_bpfi_faulty_vctr(sig_num) = hi_bpfi ;

所使用数据如下：

主代码如下：

% load data
load('bearing.mat')

% processing of the healthy signals
num_healthy_sigs = size(sigs_healthy_t, 2) ; % number of healthy signals
% pre-allocation
healthy_ftf_values = zeros(num_healthy_sigs, 1) ; hi_ftf_healthy_vctr = zeros(num_healthy_sigs, 1) ;
healthy_bsf_values = zeros(num_healthy_sigs, 1) ; hi_bsf_healthy_vctr = zeros(num_healthy_sigs, 1) ;
healthy_bpfo_values = zeros(num_healthy_sigs, 1) ; hi_bpfo_healthy_vctr = zeros(num_healthy_sigs, 1) ;
healthy_bpfi_values = zeros(num_healthy_sigs, 1) ; hi_bpfi_healthy_vctr = zeros(num_healthy_sigs, 1) ;
for sig_num = 1 : num_healthy_sigs
    
    sig_healthy_t = sigs_healthy_t(:, sig_num) ; % signal in the time domain
    speed_healthy_gear_t = speed_healthy_gear(:, sig_num) ; % speed vector of the gear
    speed_healthy_bearing_t = speed_healthy_bearing(:, sig_num) ;  % speed vector of the bearing
    t = [0 : dt : (length(sig_healthy_t)-1)*dt].' ; % time vector
    
    % dephase
    num_of_fllwing_sgmnts_2_average = 20 ;
    sig_after_dephase = dephase(t, speed_healthy_gear_t, sig_healthy_t, num_of_fllwing_sgmnts_2_average) ;

    % angular resampling
    [sig_cyc, cyc_fs] = angular_resampling(t, speed_healthy_bearing_t, sig_after_dephase) ;
    dcyc = 1 / cyc_fs ;
    
    % envelope calculation
    sig_env_cyc = abs(hilbert(sig_cyc).^2) ;

    % convert envelope from the cycle domain to the order domain
    sig_env_order = fft(sig_env_cyc) ;
    dorder = 1 / (length(sig_env_cyc)*dcyc) ;

    % features extraction
    shaft_speed = 1 ; % in the cycle domain the shaft speed is 1.
    [ftf, bsf, bpfo, bpfi] = calc_bearing_tones(shaft_speed, bearing_specifications.num_balls, ...
    bearing_specifications.ball_diameter, bearing_specifications.pitch_diameter, bearing_specifications.bearing_contact_angle) ;
    bearing_tones = [ftf, bsf, bpfo, bpfi] ;
    
    order_sensativity = 0.02 ; % sensitivity of the order value for searching the bearing tone pick
    first_ftf_value = extract_specific_bearing_tone(dorder, sig_env_order, ftf, order_sensativity) ;
    first_bsf_value = extract_specific_bearing_tone(dorder, sig_env_order, bsf, order_sensativity) ;
    first_bpfo_value = extract_specific_bearing_tone(dorder, sig_env_order, bpfo, order_sensativity) ;
    first_bpfi_value = extract_specific_bearing_tone(dorder, sig_env_order, bpfi, order_sensativity) ;
    
    healthy_ftf_values(sig_num) = first_ftf_value ;
    healthy_bsf_values(sig_num) = first_bsf_value ;
    healthy_bpfo_values(sig_num) = first_bpfo_value ;
    healthy_bpfi_values(sig_num) = first_bpfi_value ;

end % of for

healthy_ftf_average = mean(healthy_ftf_values) ; healthy_ftf_std = std(healthy_ftf_values) ;
healthy_bsf_average = mean(healthy_bsf_values) ; healthy_bsf_std = std(healthy_bsf_values) ;
healthy_bpfo_average = mean(healthy_bpfo_values) ; healthy_bpfo_std = std(healthy_bpfo_values) ;
healthy_bpfi_average = mean(healthy_bpfi_values) ; healthy_bpfi_std = std(healthy_bpfi_values) ;

for sig_num = 1 : num_healthy_sigs
    
    % HI calculation
    hi_ftf = mean(abs(healthy_ftf_values(sig_num) - healthy_ftf_average)./healthy_ftf_std) ;
    hi_ftf_healthy_vctr(sig_num) = hi_ftf ;
    hi_bsf = mean(abs(healthy_bsf_values(sig_num) - healthy_bsf_average)./healthy_bsf_std) ;
    hi_bsf_healthy_vctr(sig_num) = hi_bsf ;
    hi_bpfo = mean(abs(healthy_bpfo_values(sig_num) - healthy_bpfo_average)./healthy_bpfo_std) ;
    hi_bpfo_healthy_vctr(sig_num) = hi_bpfo ;
    hi_bpfi = mean(abs(healthy_bpfi_values(sig_num) - healthy_bpfi_average)./healthy_bpfi_std) ;
    hi_bpfi_healthy_vctr(sig_num) = hi_bpfi ;

end % of for

% processing of the faulty signals
num_faulty_sigs = size(sigs_faulty_t, 2) ;
hi_ftf_faulty_vctr = zeros(num_faulty_sigs, 1) ;
hi_bsf_faulty_vctr = zeros(num_faulty_sigs, 1) ;
hi_bpfo_faulty_vctr = zeros(num_faulty_sigs, 1) ;
hi_bpfi_faulty_vctr = zeros(num_faulty_sigs, 1) ;
for sig_num = 1 : num_faulty_sigs
    
    sig_faulty_t = sigs_faulty_t(:, sig_num) ; % signal in the time domain
    speed_faulty_gear_t = speed_faulty_gear(:, sig_num) ; % speed vector of the gear
    speed_faulty_bearing_t = speed_faulty_bearing(:, sig_num) ; % speed vector of the bearing
    t = [0 : dt : (length(sig_faulty_t)-1)*dt].' ; % time vector
    
    % dephase
    num_of_fllwing_sgmnts_2_average = 20 ;
    sig_after_dephase = dephase(t, speed_faulty_gear_t, sig_faulty_t, num_of_fllwing_sgmnts_2_average) ;

    % angular resampling
    [sig_cyc, cyc_fs] = angular_resampling(t, speed_faulty_bearing_t, sig_after_dephase) ;
    dcyc = 1 / cyc_fs ;
    
    % envelope calculation
    sig_env_cyc = abs(hilbert(sig_cyc).^2) ;

    % convert envelope from the cycle domain to the order domain
    sig_env_order = fft(sig_env_cyc) ;
    dorder = 1 / (length(sig_env_cyc)*dcyc) ;

    % features extraction
    shaft_speed = 1 ; % in the cycle domain the shaft speed is 1.
    [ftf, bsf, bpfo, bpfi] = calc_bearing_tones(shaft_speed, bearing_specifications.num_balls, ...
    bearing_specifications.ball_diameter, bearing_specifications.pitch_diameter, bearing_specifications.bearing_contact_angle) ;
    bearing_tones = [ftf, bsf, bpfo, bpfi] ;
    
    order_sensativity = 0.02 ; % sensitivity of the order value for searching the bearing tone pick
    first_ftf_value = extract_specific_bearing_tone(dorder, sig_env_order, ftf, order_sensativity) ;
    first_bsf_value = extract_specific_bearing_tone(dorder, sig_env_order, bsf, order_sensativity) ;
    first_bpfo_value = extract_specific_bearing_tone(dorder, sig_env_order, bpfo, order_sensativity) ;
    first_bpfi_value = extract_specific_bearing_tone(dorder, sig_env_order, bpfi, order_sensativity) ;

    % HI calculation
    hi_ftf = mean(abs(first_ftf_value - healthy_ftf_average)./healthy_ftf_std) ;
    hi_ftf_faulty_vctr(sig_num) = hi_ftf ;
    hi_bsf = mean(abs(first_bsf_value - healthy_bsf_average)./healthy_bsf_std) ;
    hi_bsf_faulty_vctr(sig_num) = hi_bsf ;
    hi_bpfo = mean(abs(first_bpfo_value - healthy_bpfo_average)./healthy_bpfo_std) ;
    hi_bpfo_faulty_vctr(sig_num) = hi_bpfo ;
    hi_bpfi = mean(abs(first_bpfi_value - healthy_bpfi_average)./healthy_bpfi_std) ;
    hi_bpfi_faulty_vctr(sig_num) = hi_bpfi ;

end % of for


% ----------------------------------------------------------------------- %
% Part for figures

axis_font_size = 15 ;
title_font_size = 30 ;
axis_name_font_size = 25 ;
lgnd_font_size = 20 ;

figure
plot([1:num_healthy_sigs], hi_ftf_healthy_vctr, 'LineWidth', 3, 'Color', 'g') ;
hold on
plot([num_healthy_sigs+1:num_healthy_sigs+num_faulty_sigs], hi_ftf_faulty_vctr, 'LineWidth', 3, 'Color', [0.7 0.7 0.7]) ;
hold off
ax = gca;
ax.FontSize = axis_font_size;
title('Health indicator of the cage', 'FontName', 'Times New Roman', 'FontSize', title_font_size)
xlabel('Record number', 'FontName', 'Times New Roman', 'FontSize', axis_name_font_size)
ylabel('HI value', 'FontName', 'Times New Roman', 'FontSize', axis_name_font_size)
legend('Healthy records', 'Faulty records');

figure
plot([1:num_healthy_sigs], hi_bsf_healthy_vctr, 'LineWidth', 3, 'Color', 'g') ;
hold on
plot([num_healthy_sigs+1:num_healthy_sigs+num_faulty_sigs], hi_bsf_faulty_vctr, 'LineWidth', 3, 'Color', [0.7 0.7 0.7]) ;
hold off
ax = gca;
ax.FontSize = axis_font_size;
title('Health indicator of the rolling elements', 'FontName', 'Times New Roman', 'FontSize', title_font_size)
xlabel('Record number', 'FontName', 'Times New Roman', 'FontSize', axis_name_font_size)
ylabel('HI value', 'FontName', 'Times New Roman', 'FontSize', axis_name_font_size)
legend('Healthy records', 'Faulty records');

figure
plot([1:num_healthy_sigs], hi_bpfo_healthy_vctr, 'LineWidth', 3, 'Color', 'g') ;
hold on
plot([num_healthy_sigs+1:num_healthy_sigs+num_faulty_sigs], hi_bpfo_faulty_vctr, 'LineWidth', 3, 'Color', 'r') ;
hold off
ax = gca;
ax.FontSize = axis_font_size;
title('Health indicator of the outer race', 'FontName', 'Times New Roman', 'FontSize', title_font_size)
xlabel('Record number', 'FontName', 'Times New Roman', 'FontSize', axis_name_font_size)
ylabel('HI value', 'FontName', 'Times New Roman', 'FontSize', axis_name_font_size)
legend('Healthy records', 'Faulty records');

figure
plot([1:num_healthy_sigs],hi_bpfi_healthy_vctr, 'LineWidth', 3, 'Color', 'g') ;
hold on
plot([num_healthy_sigs+1:num_healthy_sigs+num_faulty_sigs], hi_bpfi_faulty_vctr, 'LineWidth', 3, 'Color', [0.7 0.7 0.7]) ;
hold off
ax = gca;
ax.FontSize = axis_font_size;
title('Health indicator of the inner race', 'FontName', 'Times New Roman', 'FontSize', title_font_size)
xlabel('Record number', 'FontName', 'Times New Roman', 'FontSize', axis_name_font_size)
ylabel('HI value', 'FontName', 'Times New Roman', 'FontSize', axis_name_font_size)
legend('Healthy records', 'Faulty records');

出图如下：