Rolling Element Bearing Fault Diagnosis Using Deep Learning
使用深度学习的滚动轴承故障诊断

Machinery Failure Prevention Technology (MFPT) Challenge Data
机械故障预防技术 （MFPT） 挑战数据

MFPT Challenge data [2] contains 23 data sets collected from machines under various fault conditions. The first 20 data sets are collected from a bearing test rig, with three under good conditions, three with outer race faults under constant load, seven with outer race faults under various loads, and seven with inner race faults under various loads. The remaining three data sets are from real-world machines: an oil pump bearing, an intermediate speed bearing, and a planet bearing. The fault locations are unknown. In this example, you use only the data collected from the test rig with known conditions.
MFPT Challenge 数据 [2] 包含从各种故障条件下的机器收集的 23 个数据集。前 20 个数据集是从轴承试验台收集的，其中 3 个处于良好状态，3 个在恒定载荷下存在外圈故障，7 个在各种载荷下存在外圈故障，7 个在各种载荷下存在内圈故障。其余三个数据集来自真实世界的机器：油泵轴承、中速轴承和行星轴承。故障位置未知。在此示例中，您仅使用从具有已知条件的测试库位收集的数据。

Each data set contains an acceleration signal gs, sampling rate sr, shaft speed rate, load weight load, and four critical frequencies representing different fault locations: ball pass frequency outer race (BPFO), ball pass frequency inner race (BPFI), fundamental train frequency (FTF), and ball spin frequency (BSF). The formulas for BPFO and BPFI are as follows [1].
每个数据集包含一个加速度信号 gs 、采样率 sr 、轴速 rate 、负载重量 load 和四个临界频率，分别代表不同的故障位置：传球频率外圈 （BPFO）、传球频率内圈 （BPFI）、基线序列频率 （FTF） 和球自旋频率 （BSF）。BPFO和BPFI的公式如下[1]。

As shown in the figure, d is the ball diameter and D is the pitch diameter. The variable fr
如图所示，d为球径，D为节距直径。变量 f is the shaft speed, n is the number of rolling elements, and ϕ is the bearing contact angle [1].
是轴转速，n是滚动体的数量，φ是轴承接触角[1]。

Scalogram of Bearing Data
轴承数据的尺度图

To benefit from pretrained CNN deep networks, use the plotBearingSignalAndScalogram helper function to convert 1-D vibration signals in the MFPT dataset to 2-D scalograms. A scalogram is a time-frequency domain representation of the original time-domain signal [3]. The two dimensions in a scalogram image represent time and frequency. To visualize the relationship between a scalogram and its original vibration signal, plot the vibration signal with an inner race fault against its scalogram.
要从预训练的 CNN 深度网络中受益，请使用 plotBearingSignalAndScalogram helper 函数将 MFPT 数据集中的一维振动信号转换为二维标度图。标度图是原始时域信号的时频域表示[3]。标度图图像中的两个维度表示时间和频率。要可视化标度图与其原始振动信号之间的关系，请根据其标度图绘制具有内部竞争故障的振动信号。

% Import data with inner race fault
data_inner = load(fullfile(matlabroot, 'toolbox', 'predmaint', ...
    'predmaintdemos', 'bearingFaultDiagnosis', ...
    'train_data', 'InnerRaceFault_vload_1.mat'));
% Plot bearing signal and scalogram
plotBearingSignalAndScalogram(data_inner)

During the 0.1 seconds shown in the plot, the vibration signal contains 12 impulses because the tested bearing's BPFI is 118.875 Hz. Accordingly, the scalogram shows 12 distinct peaks that align with the impulses in the vibration signal. Next, visualize scalograms for the outer race fault.
在图中所示的 0.1 秒内，振动信号包含 12 个脉冲，因为被测轴承的 BPFI 为 118.875 Hz。 因此，标度图显示了 12 个不同的峰值，这些峰值与振动信号中的脉冲一致。接下来，可视化外竞相断层的标度图。

% Import data with outer race fault
data_outer = load(fullfile(matlabroot, 'toolbox', 'predmaint', ...
    'predmaintdemos', 'bearingFaultDiagnosis', ...
    'test_data', 'OuterRaceFault_3.mat'));
% Plot original signal and its scalogram
plotBearingSignalAndScalogram(data_outer)

The scalogram of the outer race fault shows 8 distinct peaks during the first 0.1 seconds, which is consistent with the ballpass frequencies. Because the impulses in the time-domain signal is not as dominant as in the inner race fault case, the distinct peaks in the scalogram show less contrast with the background. The scalogram of the normal condition does not show dominant distinct peaks.
外圈断层的标度图在前0.1秒内显示8个不同的峰值，这与传球频率一致。由于时域信号中的脉冲不如内部竞争故障情况下的脉冲占主导地位，因此标度图中的不同峰值与背景的对比度较低。正常条件下的标度图不显示显性明显峰。

% Import normal bearing data
data_normal = load(fullfile(matlabroot, 'toolbox', 'predmaint', ...
    'predmaintdemos', 'bearingFaultDiagnosis', ...
    'train_data', 'baseline_1.mat'));
% Plot original signal and its scalogram
plotBearingSignalAndScalogram(data_normal)

The number of distinct peaks is a good feature to differentiate between inner race faults, outer race faults, and normal conditions. Therefore, a scalogram can be a good candidate for classifying bearing faults. In this example, all bearing signal measurements come from tests using the same shaft speed. To apply this example to bearing signals under different shaft speeds, the data needs to be normalized by shaft speed. Otherwise, the number of "pillars" in the scalogram will be wrong.
不同峰的数量是区分内圈断层、外圈断层和正常条件的一个很好的特征。因此，标度图可以很好地对轴承故障进行分类。在此示例中，所有轴承信号测量均来自使用相同轴转速的测试。为了将此示例应用于不同轴转速下的轴承信号，需要按轴速度对数据进行归一化。否则，时标图中的“柱子”数量将是错误的。

Prepare Training Data 准备训练数据

Unzip the downloaded file.
解压缩下载的文件。

if exist('RollingElementBearingFaultDiagnosis-Data-master.zip', 'file')
    unzip('RollingElementBearingFaultDiagnosis-Data-master.zip')
end

The downloaded dataset contains a training dataset with 14 MAT-files (2 normal, 5 inner race fault, 7 outer race fault) and a testing dataset with 6 MAT-files (1 normal, 2 inner race fault, 3 outer race fault).
下载的数据集包含一个包含 14 个 MAT 文件（2 个正常、5 个内部争用故障、7 个外部争用故障）的训练数据集和一个包含 6 个 MAT 文件（1 个正常、2 个内部争用故障、3 个外部争用故障）的测试数据集。

By assigning function handles to ReadFcn, the file ensemble datastore can navigate into the files to retrieve data in the desired format. For example, the MFPT data has a structure bearing that stores the vibration signal gs, sampling rate sr, and so on. Instead of returning the bearing structure itself, the readMFPTBearing function is written so that the file ensemble datastore returns the vibration signal gs inside of the bearing data structure.
通过将 ReadFcn 函数句柄分配给 ，文件集成数据存储可以导航到文件以检索所需格式的数据。例如，MFPT数据具有存储振动信号 gs 、采样率 sr 等的结构 bearing 。编写函数 readMFPTBearing 不是返回轴承结构本身，而是文件集成数据存储返回 bearing 数据结构内部的振动信号 gs 。

fileLocation = fullfile('.', 'RollingElementBearingFaultDiagnosis-Data-master', 'train_data');
fileExtension = '.mat';
ensembleTrain = fileEnsembleDatastore(fileLocation, fileExtension);
ensembleTrain.ReadFcn = @readMFPTBearing;
ensembleTrain.DataVariables = ["gs", "sr", "rate", "load", "BPFO", "BPFI", "FTF", "BSF"];
ensembleTrain.ConditionVariables = ["Label", "FileName"];
ensembleTrain.SelectedVariables = ["gs", "sr", "rate", "load", "BPFO", "BPFI", "FTF", "BSF", "Label", "FileName"]

ensembleTrain = 
  fileEnsembleDatastore with properties:

                 ReadFcn: @readMFPTBearing
        WriteToMemberFcn: []
           DataVariables: [8×1 string]
    IndependentVariables: [0×0 string]
      ConditionVariables: [2×1 string]
       SelectedVariables: [10×1 string]
                ReadSize: 1
              NumMembers: 14
          LastMemberRead: [0×0 string]
                   Files: [14×1 string]

Now, convert the 1-D vibration signals to scalograms and save the images for training. The size of each scalogram is 227-by-227-by-3, which is the same input size required by SqueezeNet. To improve accuracy, the helper function convertSignalToScalogram envelops the raw signal and divides it into multiple segments. After running the following commands, a folder named "train_image" appears in the current folder. All scalogram images of the bearing signals in the "RollingElementBearingFaultDiagnosis-Data-master/train_data" folder are saved in the "train_image" folder.
现在，将一维振动信号转换为标度图并保存图像以进行训练。每个比例尺的大小为 227 x 227 x 3，与 SqueezeNet 所需的输入大小相同。为了提高精度，辅助函数 convertSignalToScalogram 将原始信号封装起来，并将其分成多个段。运行以下命令后，当前文件夹中会出现一个名为“train_image”的文件夹。“RollingElementBearingFaultDiagnosis-Data-master/train_data”文件夹中轴承信号的所有比例图图像都保存在“train_image”文件夹中。

reset(ensembleTrain)
while hasdata(ensembleTrain)
  folderName = 'train_image';
  convertSignalToScalogram(ensembleTrain,folderName);
end

Create an image datastore and split the training data into training and validation data sets, using 80% of the images from the "train_image" folder for training and 20% for validation.
创建图像数据存储并将训练数据拆分为训练和验证数据集，使用“train_image”文件夹中 80% 的图像进行训练，20% 用于验证。

% Create image datastore to store all training images
path = fullfile('.', folderName);
imds = imageDatastore(path, ...
  'IncludeSubfolders',true,'LabelSource','foldernames');
% Use 20% training data as validation set
[imdsTrain,imdsValidation] = splitEachLabel(imds,0.8,'randomize');

Train Network with Transfer Learning
通过迁移学习训练网络

Next, fine-tune the pretrained SqueezeNet convolutional neural network to perform classification on the scalograms. SqueezeNet has been trained on over a million images and has learned rich feature representations. Transfer learning is commonly used in deep learning applications. You can take a pretrained network and use it as a starting point for a new task. Fine-tuning a network with transfer learning is usually much faster and easier than training a network with randomly initialized weights from scratch. You can quickly transfer learned features using a smaller number of training images. Load and view the SqueezeNet network:
接下来，微调预训练的 SqueezeNet 卷积神经网络以对标度图执行分类。SqueezeNet 已经接受了超过 100 万张图像的训练，并学习了丰富的特征表示。迁移学习通常用于深度学习应用。您可以采用预训练网络并将其用作新任务的起点。使用迁移学习微调网络通常比从头开始训练具有随机初始化权重的网络更快、更容易。您可以使用较少数量的训练图像快速传输学习的特征。加载并查看 SqueezeNet 网络：

net = imagePretrainedNetwork("squeezenet")

net = 
  dlnetwork with properties:

         Layers: [68×1 nnet.cnn.layer.Layer]
    Connections: [75×2 table]
     Learnables: [52×3 table]
          State: [0×3 table]
     InputNames: {'data'}
    OutputNames: {'prob_flatten'}
    Initialized: 1

  View summary with summary.

analyzeNetwork(net)

SqueezeNet uses the convolutional layer 'conv10' to extract image features and classify the input image. This layer contains information to combine the features that the network extracts into class probabilities, a loss value, and predicted labels. To retrain SqueezeNet for classifying new images, the convolutional layers 'conv10' need to be replaced with a new layer adapted to the bearing images.
SqueezeNet 使用卷积层 ' conv10 ' 来提取图像特征并对输入图像进行分类。此图层包含用于将网络提取的特征组合为类概率、损失值和预测标注的信息。为了重新训练 SqueezeNet 对新图像进行分类，卷积层 conv10 需要替换为适合方位角图像的新层。

In most networks, the last layer with learnable weights is a fully connected layer. In some networks, such as SqueezeNet, the last learnable layer is a 1-by-1 convolutional layer instead. In this case, replace the convolutional layer with a new convolutional layer with a number of filters equal to the number of classes.
在大多数网络中，具有可学习权重的最后一层是全连接层。在某些网络中，例如 SqueezeNet，最后一个可学习层是 1×1 卷积层。在这种情况下，将卷积层替换为新的卷积层，该卷积层的滤波器数等于类数。

numClasses = numel(categories(imdsTrain.Labels));

newConvLayer = convolution2dLayer([1, 1],numClasses,'WeightLearnRateFactor',10,'BiasLearnRateFactor',10,"Name",'new_conv');
net = replaceLayer(net,'conv10',newConvLayer);

Specify the training options. To slow down learning in the transferred layers, set the initial learning rate to a small value. When you create the convolutional layer, you include larger learning rate factors to speed up learning in the new final layers. This combination of learning rate settings results in fast learning only in the new layers and slower learning in the other layers. When performing transfer learning, you do not need to train for as many epochs. An epoch is a full training cycle on the entire training data set. The software validates the network every ValidationFrequency iterations during training.
指定训练选项。要减慢传输层中的学习速度，请将初始学习速率设置为较小的值。创建卷积层时，需要包括更大的学习率因子，以加快新最终层的学习速度。这种学习速率设置的组合导致仅在新层中快速学习，而在其他层中学习速度较慢。在执行迁移学习时，您不需要训练那么多的时期。epoch 是整个训练数据集上的完整训练周期。该软件在训练期间每次 ValidationFrequency 迭代都会验证网络。

options = trainingOptions('sgdm', ...
  'InitialLearnRate',0.0001, ...
  'MaxEpochs',4, ...
  'Shuffle','every-epoch', ...
  'ValidationData',imdsValidation, ...
  'ValidationFrequency',30, ...
  'Verbose',false, ...
  'MiniBatchSize',20, ...
  'Plots','training-progress');

Train the new network and set the loss function to 'crossentropy' (for classification problems). By default, trainnet uses a GPU if you have Parallel Computing Toolbox™ and a supported GPU device. For information on supported devices, see GPU Computing Requirements (Parallel Computing Toolbox). Otherwise, trainnet uses a CPU. You can also specify the execution environment by using the 'ExecutionEnvironment' name-value argument of trainingOptions.
训练新网络， 'crossentropy' ( 并将分类问题的损失函数设置为）。默认情况下，如果您有 Parallel Computing Toolbox™ 和支持的 GPU 设备， trainnet 则使用 GPU。有关受支持设备的信息，请参阅 GPU 计算要求 （Parallel Computing Toolbox）。否则， trainnet 使用 CPU。还可以使用 的 'ExecutionEnvironment' name-value 参数指定 trainingOptions 执行环境。

net = trainnet(imdsTrain,net,"crossentropy",options);

Validate Using Test Data Sets
使用测试数据集进行验证

Use bearing signals in the "RollingElementBearingFaultDiagnosis-Data-master/test_data" folder to validate the accuracy of the trained network. The test data needs to be processed in the same way as the training data.
使用“RollingElementBearingFaultDiagnosis-Data-master/test_data”文件夹中的方位角信号来验证经过训练的网络的准确性。测试数据的处理方式与训练数据相同。

Create a file ensemble datastore to store the bearing vibration signals in the test folder.
创建文件集成数据存储以将轴承振动信号存储在测试文件夹中。

fileLocation = fullfile('.', 'RollingElementBearingFaultDiagnosis-Data-master', 'test_data');
fileExtension = '.mat';
ensembleTest = fileEnsembleDatastore(fileLocation, fileExtension);
ensembleTest.ReadFcn = @readMFPTBearing;
ensembleTest.DataVariables = ["gs", "sr", "rate", "load", "BPFO", "BPFI", "FTF", "BSF"];
ensembleTest.ConditionVariables = ["Label", "FileName"];
ensembleTest.SelectedVariables = ["gs", "sr", "rate", "load", "BPFO", "BPFI", "FTF", "BSF", "Label", "FileName"];

Convert 1-D signals to 2-D scalograms.
将一维信号转换为二维标度图。

reset(ensembleTest)
while hasdata(ensembleTest)
  folderName = 'test_image';
  convertSignalToScalogram(ensembleTest,folderName);
end

Create an image datastore to store the test images.
创建映像数据存储以存储测试映像。

path = fullfile('.','test_image');
imdsTest = imageDatastore(path, ...
  'IncludeSubfolders',true,'LabelSource','foldernames');

Classify the test image datastore with the trained network.
使用经过训练的网络对测试图像数据存储进行分类。

scores = minibatchpredict(net,imdsTest,'MiniBatchSize',20);
YPred = scores2label(scores,unique(imdsTrain.Labels));

Compute the accuracy of the prediction.
计算预测的准确性。

YTest = imdsTest.Labels;
accuracy = sum(YPred == YTest)/numel(YTest)

accuracy = 0.9807

Plot a confusion matrix.
绘制混淆矩阵。

figure
confusionchart(YTest,YPred)

When you train the network multiple times, you might see some variation in accuracy between trainings, but the average accuracy should be around 98%. Even though the training set is quite small, this example benefits from transfer learning and achieves good accuracy.
当您多次训练网络时，您可能会看到训练之间的准确性存在一些差异，但平均准确率应该在 98% 左右。尽管训练集非常小，但此示例受益于迁移学习并实现了良好的准确性。

Conclusion 结论

This example demonstrates that deep learning can be an effective tool to identify different types of faults in rolling element bearing, even when the data size is relatively small. A deep learning approach reduces the time that conventional approach requires for feature engineering. For comparison, see the example Rolling Element Bearing Fault Diagnosis.
这个例子表明，即使数据量相对较小，深度学习也可以成为识别滚动轴承中不同类型的故障的有效工具。深度学习方法减少了传统方法进行特征工程所需的时间。有关比较，请参阅滚动体轴承故障诊断示例。

References 引用

[1] Randall, Robert B., and Jérôme Antoni. “Rolling Element Bearing Diagnostics—A Tutorial.” Mechanical Systems and Signal Processing 25, no. 2 (February 2011): 485–520. https://doi.org/10.1016/j.ymssp.2010.07.017.
[1] 兰德尔、罗伯特 B. 和杰罗姆·安东尼。“滚动轴承诊断 — 教程。”机械系统与信号处理 25，第 2 期（2011 年 2 月）：485–520。https://doi.org/10.1016/j.ymssp.2010.07.017。

[2] Bechhoefer, Eric. "Condition Based Maintenance Fault Database for Testing Diagnostics and Prognostic Algorithms." 2013. https://www.mfpt.org/fault-data-sets/.
[2] 贝克霍弗，埃里克。“基于状态的维护故障数据库，用于测试诊断和预测算法。”2013. https://www.mfpt.org/fault-data-sets/.

[3] Verstraete, David, Andrés Ferrada, Enrique López Droguett, Viviana Meruane, and Mohammad Modarres. “Deep Learning Enabled Fault Diagnosis Using Time-Frequency Image Analysis of Rolling Element Bearings.” Shock and Vibration 2017 (2017): 1–17. https://doi.org/10.1155/2017/5067651.
[3] Verstraete、David、Andrés Ferrada、Enrique López Droguett、Viviana Meruane 和 Mohammad Modarres。“深度学习支持使用滚动轴承的时频图像分析进行故障诊断。”冲击与振动 2017 （2017）： 1–17.https://doi.org/10.1155/2017/5067651。

Helper Functions 帮助程序函数

function plotBearingSignalAndScalogram(data)
% Convert 1-D bearing signals to scalograms through wavelet transform
fs = data.bearing.sr;
t_total = 0.1; % seconds
n = round(t_total*fs);
bearing = data.bearing.gs(1:n);
[cfs,frq] = cwt(bearing,'amor', fs);

% Plot the original signal and its scalogram
figure
subplot(2,1,1)
plot(0:1/fs:(n-1)/fs,bearing)
xlim([0,0.1])
title('Vibration Signal')
xlabel('Time (s)')
ylabel('Amplitude')
subplot(2,1,2)
surface(0:1/fs:(n-1)/fs,frq,abs(cfs))
shading flat
xlim([0,0.1])
ylim([0,max(frq)])
title('Scalogram')
xlabel('Time (s)')
ylabel('Frequency (Hz)')
end

function convertSignalToScalogram(ensemble,folderName)
% Convert 1-D signals to scalograms and save scalograms as images
data = read(ensemble);
fs = data.sr;
x = data.gs{:};
label = char(data.Label);
fname = char(data.FileName);
ratio = 5000/97656;
interval = ratio*fs;
N = floor(numel(x)/interval);

% Create folder to save images
path = fullfile('.',folderName,label);
if ~exist(path,'dir')
  mkdir(path);
end

for idx = 1:N
  sig = envelope(x(interval*(idx-1)+1:interval*idx));
  cfs = cwt(sig,'amor', seconds(1/fs));
  cfs = abs(cfs);
  img = ind2rgb(round(rescale(flip(cfs),0,255)),jet(320));
  outfname = fullfile('.',path,[fname '-' num2str(idx) '.jpg']);
  imwrite(imresize(img,[227,227]),outfname);
end
end