MATLAB环境下基于电机电流特征分析的齿轮故障检测

做多了基于机器学习和深度学习的机械故障诊断，感觉没意思了，换个口味，写一下基于现代信号处理的齿轮故障检测。本例来源于Mathworks公司的某一个大佬，主要讲解如何利用电机的电流信号识别齿轮系中的故障。

运行环境：MATLAB R2021B，低于此版本可能不能运行成功

完整代码面包多下载：

https://mianbaoduo.com/o/bread/Y5icmp5r

使用传统振动传感器进行齿轮的故障监测具有一定的挑战性，尤其是当振动传感器不易安装时。由于电流信号很容易获取，因此本例说明如何应用电机电流特征分析提取频谱特征，以检测电机驱动齿轮中的故障。

硬件

本例中的电流数据是Futaba S3003伺服器中采集的，该伺服器经过修改可以实现连续旋转。伺服系统将内部直流电机的高转速转换为输出的高扭矩。伺服系统由直流电机、一组尼龙金属齿轮和控制电路组成。本例移除了控制电路，允许直接监控直流电机的电流信号，使用红外光断续器采集伺服器输出的转速信号。伺服器和红外光断续器由3D 打印支架固定。

直流电机以恒定5伏电压驱动，四对齿轮提供 278:1 的减速比，花键处的轴转速约为 50 rpm。 由于电流测量值的变化太小而无法检测，因此使用单电源传感器接口放大器放大了电流信号。 然后使用抗混叠五阶椭圆低通滤波器对放大的电流信号进行滤波，使其平滑并消除噪声，然后通过模数转换器 (ADC) 将其发送到 Arduino Uno。

如上面流程图所示，首先分别使用放大器和抗混叠低通滤波器对电流信号进行放大和滤波。 Arduino Uno 通过 ADC 以 1.5 kHz 采样电流信号，并将其与转速计脉冲一起作为串行数据以 115200 的波特率传输到计算机。MATLAB 脚本从 Arduino Uno 获取串行数据并将其写入逗号分隔值 (CSV) 文件。 然后使用 MATLAB 读取和处理CSV 文件以提取频谱指标。

伺服齿轮系

Futaba S3003 由四对尼龙齿轮组成，如上图所示。 直流电机轴上的小齿轮 P1 与阶梯齿轮 G1 啮合，小齿轮P2与阶梯齿轮G2啮合，小齿轮P3与阶梯齿轮G3啮合，小齿轮 P4 与连接到输出花键的齿轮 G4 啮合。 阶梯式齿轮组 G1 和 P2、G2 和 P3 以及 G3 和 P4 是自由旋转齿轮，即它们不连接到各自的轴上。 当电机以 5 伏电压驱动时，该组驱动齿轮提供了 278:1 的减速比，从 13901 rpm 的电机速度到输出花键的50 rpm转速。 下表概述了每个齿轮啮合处的输出速度、齿轮啮合频率和累积齿轮减速的齿数和理论值。

在将故障引入阶梯齿轮 G2 和 G3 之前，共收集了 10 个健康数据集。 由于齿轮是用尼龙模制而成的，因此在齿空间中切割槽，在两个齿轮中引入了模拟裂纹。 齿距是沿正齿轮的节圆测量的两个相邻齿之间的间隙。 槽深度约为齿轮半径的 70%。 引入齿轮G2和G3的故障后，共记录了10个故障数据集。

数据可视化

考虑一组健康和故障数据，可视化当前信号和转速脉冲

healthyData = servoData.healthyData{1,1};
faultyData = servoData.faultyData{1,1};
healthyCurrent = healthyData.Signal;
faultyCurrent = faultyData.Signal;
healthyTacho = healthyData.Pulse;
faultyTacho = faultyData.Pulse;
tHealthy = healthyData.Time;
tFaulty = faultyData.Time;

figure
ax1 = subplot(221);
plot(tHealthy,healthyCurrent)
title('Current Signal - Healthy Gears')
ylabel('Current (mA)')
ax2 = subplot(222);
plot(tFaulty,faultyCurrent)
title('Current Signal - Faulty Gears')
ylabel('Current (mA)')
ax3 = subplot(223);
plot(tHealthy,healthyTacho)
title('Tachometer Pulse - Healthy Gears')
ylabel('Pulses, 8/rev')
ax4 = subplot(224);
plot(tFaulty,faultyTacho)
title('Tachometer Pulse - Faulty Gears')
ylabel('Pulses, 8/rev')
linkaxes([ax1,ax2,ax3,ax4],'x');
ax1.XLim = seconds([0 2]);

输出花键每转八个脉冲，对应于伺服轮上的八个孔。

计算额定转速

计算额定速度，以检测齿轮系统中感兴趣的频率，并将其与功率谱上的频率正确匹配，感兴趣的频率是以Hz为单位的实际输出速度值。使用1500Hz的采样频率值，可视化转速计信号，并使用tachorpm函数计算额定转速。

fs = 1500;
figure
tachorpm(healthyTacho,fs,'PulsesPerRev',8)
figure
tachorpm(faultyTacho,fs,'PulsesPerRev',8)
rpmHealthy = mean(tachorpm(healthyTacho,fs,'PulsesPerRev',8))
rpmFaulty = mean(tachorpm(faultyTacho,fs,'PulsesPerRev',8))

rpmHealthy = 49.8550

rpmFaulty = 49.5267

在正常数据集和故障数据集之间，输出轴速度有非常小的差异。实际额定rpm值也接近理论值50 rpm。

构造频带

构造频带是计算频谱指标的重要前提。利用齿轮系中驱动齿轮的齿数和名义转数，首先计算感兴趣的频率，感兴趣的频率是以Hz为单位的实际输出速度值，其值接近下表中列出的理论值。

%齿数
G4 = 41; G3 = 35; G2 = 50; G1 = 62;
P4 = 16; P3 = 10; P2 = 10; P1 = 10;
rpm = rpmHealthy;
%感兴趣的频率
FS5 = mean(rpm)/60;
FS4 = G4/P4*FS5;
FS3 = G3/P3*FS4;
FS2 =  G2/P2*FS3;
FS1 =  G1/P1*FS2;
FS = [FS1,FS2,FS3,FS4,FS5]

接下来，使用faultBands函数构造所有输出速度的频带，其中包括以下感兴趣的频率。

FS1 at 231 Hz, its first two harmonics, and 0:1 sidebands of FS2

FS2 at 37.3 Hz, its first two harmonics, and 0:1 sidebands of FS3

FS3 at 7.5 Hz and its first two harmonics

FS4 at 2.1 Hz and its first two harmonics

[FB1,info1] = faultBands(FS1,1:2,FS2,0:1);
[FB2,info2] = faultBands(FS2,1:2,FS3,0:1);
[FB3,info3] = faultBands(FS3,1:2);
[FB4,info4] = faultBands(FS4,1:2);
FB = [FB1;FB2;FB3;FB4]

输出FB是一个16 × 2的数组，包含以上感兴趣的频率。

提取功率谱密度(PSD)数据

计算并可视化运行正常和故障数据的功率谱，利用结构信息中的信息在频谱图上绘制感兴趣的频率，并放大部分频率

蓝色表示正常频谱，红色表示故障数据频谱。从图中可以观察到故障频率振幅的增加:

1FS1 at 231 Hz, its second harmonic 2FS1 at 462 Hz, and their respective sidebands

放大图，观察以下频率的振幅增加情况:

1FS2 at 37.2 Hz and its sidebands

1FS3 at 7.5 Hz and its second harmonic 2FS3 at 15 Hz

1FS4 at 2.1 Hz and its second harmonic 2FS4 at 4.2 Hz

利用pspectrum函数分别计算健康信号和故障信号的PSD

[psdHealthy,fHealthy] = pspectrum(healthyCurrent,fs,'FrequencyResolution',0.5);
[psdFaulty,fFaulty] = pspectrum(faultyCurrent,fs,'FrequencyResolution',0.5);

计算频谱指标

使用频带数据和PSD数据，使用faultBandMetrics函数计算健康数据和故障数据的频谱指标

spectralMetrics = faultBandMetrics({[psdHealthy,fHealthy],[psdFaulty,fFaulty]},FB)

spectralMetrics=2×49 table
    PeakAmplitude1    PeakFrequency1    BandPower1    PeakAmplitude2    PeakFrequency2    BandPower2    PeakAmplitude3    PeakFrequency3    BandPower3    PeakAmplitude4    PeakFrequency4    BandPower4    PeakAmplitude5    PeakFrequency5    BandPower5    PeakAmplitude6    PeakFrequency6    BandPower6    PeakAmplitude7    PeakFrequency7    BandPower7    PeakAmplitude8    PeakFrequency8    BandPower8    PeakAmplitude9    PeakFrequency9    BandPower9    PeakAmplitude10    PeakFrequency10    BandPower10    PeakAmplitude11    PeakFrequency11    BandPower11    PeakAmplitude12    PeakFrequency12    BandPower12    PeakAmplitude13    PeakFrequency13    BandPower13    PeakAmplitude14    PeakFrequency14    BandPower14    PeakAmplitude15    PeakFrequency15    BandPower15    PeakAmplitude16    PeakFrequency16    BandPower16    TotalBandPower
    ______________    ______________    __________    ______________    ______________    __________    ______________    ______________    __________    ______________    ______________    __________    ______________    ______________    __________    ______________    ______________    __________    ______________    ______________    __________    ______________    ______________    __________    ______________    ______________    __________    _______________    _______________    ___________    _______________    _______________    ___________    _______________    _______________    ___________    _______________    _______________    ___________    _______________    _______________    ___________    _______________    _______________    ___________    _______________    _______________    ___________    ______________

      0.0020712           193.75         0.010881        0.50813            231.06         0.46652        0.0019579            272.5         0.011841       0.0020486           424.06         0.011225        0.54982               462         0.89544        0.0024293           493.69        0.0091045        0.002966           29.812        0.0035485        0.015582           37.25          0.011132       0.0028865           44.688        0.0041317        0.011896            67.062          0.0072753        0.059126             74.5           0.033568         0.013218                82          0.0079886         5.7904             7.4375           2.3115          0.068452            14.938          0.027653          0.79006             2.125           0.14382          0.09849             4.25            0.01806          3.9737    
       0.009804           192.44         0.017916         3.6921            229.44          2.9975        0.0035204           266.44         0.015639       0.0057056           421.75         0.019293         1.2974            459.69          3.2185        0.0053261           495.88         0.016296       0.0031674           28.938        0.0044271        0.023983              37          0.014447          0.0136           44.438        0.0089119        0.011419            66.625          0.0077035          0.0684               74           0.037016         0.012244            81.438          0.0075805         7.7931              7.375           3.0193           0.15692            14.812          0.058255           2.4211             2.125            0.4407          0.55167             4.25            0.10029          9.9838    

输出是FB中频率范围的一个2 × 49的指标表。第一行包含健康数据的指标，第二行包含故障数据的指标。注意，以下指标对于故障数据的值要比健康数据的值高得多:

PeakAmplitude2 for the frequency at 231 Hz with a difference of 3.1842 units

TotalBandPower with a difference of 6.01 units

因此，可以使用这两个指标创建散点图，分别对故障数据和健康数据进行分组。

创建散点图

创建散点图，使用表servoData中的所有20个数据集的两个谱指标PeakAmplitude2和TotalBandPower对健康数据和故障数据进行分类。

plotData = zeros(10,4);
for n = 1:max(size(servoData))
    hC = servoData.healthyData{n,1}.Signal;
    fC = servoData.faultyData{n,1}.Signal;
    
    [psdH,freqH] = pspectrum(hC,fs,'FrequencyResolution',0.5);
    [psdF,freqF] = pspectrum(fC,fs,'FrequencyResolution',0.5);
    
    sMetrics = faultBandMetrics({[psdH,freqH],[psdF,freqF]},FB);
    plotData(n,:) = [sMetrics{:,4}',sMetrics{:,49}'];
end
figure
scatter(plotData(:,1),plotData(:,3),[],'blue')
hold on;
scatter(plotData(:,2),plotData(:,4),[],'red')
legend('Healthy Data','Faulty Data','Location','best')
xlabel('Peak Amplitude 2')
ylabel('Total Band Power')
hold off

观察健康数据集和故障数据集被分组在散点图的不同区域。因此，通过分析伺服电机的电流特征，可以对故障和健康数据进行分类。