구름 요소 베어링 결함 진단

라이브 스크립트 열기

이 예제에서는 특히 다른 기계 컴포넌트에서 강한 마스킹 신호가 발생하는 경우에 가속도 신호를 기반으로 구름 요소(rolling element) 베어링의 결함 진단을 수행하는 방법을 보여줍니다. 이 예제에서는 포락선 스펙트럼 분석 및 스펙트럼 첨도를 적용하여 베어링 결함을 진단하는 방법과 빅데이터 응용으로까지 확장하는 방법을 보여줍니다.

문제 개요

구름 요소(rolling element) 베어링은 외륜, 내륜, 케이지 또는 구름 요소에서 국소적인 결함이 발생할 수 있습니다. 구름 요소가 외륜이나 내륜의 국소 결함 부위에 부딪힐 때 또는 구름 요소의 결함 부위가 외륜이나 내륜에 부딪힐 때 베어링과 응답 변환기 사이에 고주파수 공명이 발생합니다 [1]. 다음 그림은 구름 요소가 내륜의 국소 결함 부위에 부딪히는 모습을 보여줍니다. 여기서 문제는 다양한 유형의 결함을 어떻게 검출하고 식별할 것인가입니다.

MFPT(Machinery Failure Prevention Technology) Challenge 데이터

MFPT Challenge 데이터 [4]는 다양한 결함 상태를 보이는 기계들로부터 수집한 23개의 데이터 세트를 포함합니다. 처음 20개의 데이터 세트는 베어링 시험 장치에서 수집했습니다. 3개는 양호한 상태, 3개는 일정 하중이 가해진 결함 외륜, 7개는 변동 하중이 가해진 결함 외륜, 7개는 변동 하중이 가해진 결함 내륜에서 수집한 것입니다. 나머지 3개의 데이터 세트는 실제 기계(오일 펌프 베어링, 중간 속도 베어링, 유성기어 베어링)에서 수집되었습니다. 결함 위치는 알려지지 않았습니다. 이 예제에서는 알려진 상태를 갖는 시험 장치에서 수집된 데이터만 사용됩니다.

각 데이터 세트는 가속도 신호 "gs", 샘플링 레이트 "sr", 축 속도 "rate", 하중 무게 "load" 그리고 서로 다른 결함 위치를 나타내는 4개의 임계 주파수 BPFO(외륜 볼 통과 주파수), BPFI(내륜 볼 통과 주파수), FTF(기본 열 주파수), BSF(볼 자전 주파수)를 포함합니다. 다음은 네 개의 임계 주파수의 수식입니다 [1].

외륜 볼 통과 주파수(BPFO)

$B P F O = \frac{n f_{r}}{2} (1 - \frac{d}{D} \cos ϕ)$

내륜 볼 통과 주파수(BPFI)

$B P F I = \frac{n f_{r}}{2} (1 + \frac{d}{D} \cos ϕ)$

기본 열 주파수(FTF) - 케이지 속도라고도 함

$F T F = \frac{f_{r}}{2} (1 - \frac{d}{D} \cos ϕ)$

볼(롤러) 자전 주파수

$B S F = \frac{D f_{r}}{2 d} [1 - {(\frac{d}{D} \cos ϕ)}^{2}]$

그림에서 볼 수 있듯이 $d$ 는 볼의 직경이고 $D$ 는 피치의 직경입니다. 변수 $f_{r}$ 은 축 속도, $n$ 은 구름 요소의 개수, $ϕ$ 는 베어링 접촉각입니다 [1].

베어링 진단을 위한 포락선 스펙트럼 분석

MFPT 데이터 세트에서 축 속도는 일정하므로 축 속도 변동의 영향을 제거하기 위한 전처리 단계로서 차수 추적을 수행할 필요가 없습니다.

구름 요소가 외륜이나 내륜의 국소 결함 부위에 부딪히거나 구름 요소의 결함 부위가 외륜이나 내륜에 부딪힐 때 그 충격이 BPFO, BPFI, FTF, BSF 같은 임계 주파수를 변조시킵니다. 따라서 진폭 복조를 거쳐 생성된 포락선 신호는 원시 신호의 스펙트럼 분석에서 제공할 수 없는 진단 정보를 더 많이 전달합니다. MFPT 데이터셋의 내륜 결함 신호를 예로 들어 보겠습니다.

dataInner = load(fullfile(matlabroot, 'toolbox', 'predmaint', ...
    'predmaintdemos', 'bearingFaultDiagnosis', ...
    'train_data', 'InnerRaceFault_vload_1.mat'));

내륜 결함 원시 데이터를 시간 영역에서 시각화합니다.

xInner = dataInner.bearing.gs;
fsInner = dataInner.bearing.sr;
tInner = (0:length(xInner)-1)/fsInner;
figure
plot(tInner, xInner)
xlabel('Time, (s)')
ylabel('Acceleration (g)')
title('Raw Signal: Inner Race Fault')
xlim([0 0.1])

원시 데이터를 주파수 영역에서 시각화합니다.

figure
[pInner, fpInner] = pspectrum(xInner, fsInner);
pInner = 10*log10(pInner);
plot(fpInner, pInner)
xlabel('Frequency (Hz)')
ylabel('Power Spectrum (dB)')
title('Raw Signal: Inner Race Fault')
legend('Power Spectrum')

이번에는 원시 신호의 파워 스펙트럼에서 저주파수 범위로 확대하여 BPFI와 그 처음 몇 개의 고조파에서 주파수 응답을 자세히 살펴봅니다.

figure
plot(fpInner, pInner)
ncomb = 10;
helperPlotCombs(ncomb, dataInner.BPFI)
xlabel('Frequency (Hz)')
ylabel('Power Spectrum (dB)')
title('Raw Signal: Inner Race Fault')
legend('Power Spectrum', 'BPFI Harmonics')
xlim([0 1000])

BPFI와 그 고조파에서 명확한 패턴을 볼 수가 없습니다. 원시 신호의 주파수 분석은 유용한 진단 정보를 제공하지 않습니다.

시간 영역 데이터를 보면 원시 신호의 진폭이 특정 주파수에서 변조되며 변조의 기본 주파수가 대략 1/0.009Hz $\approx$ 111Hz임을 알 수 있습니다. 구름 요소가 내륜의 국소 결함 부위에 부딪칠 때의 주파수, 다시 말해서 BPFI는 118.875Hz인 것으로 알려져 있습니다. 이는 베어링에 내륜 결함이 있을 수 있음을 나타냅니다.

figure
subplot(2, 1, 1)
plot(tInner, xInner)
xlim([0.04 0.06])
title('Raw Signal: Inner Race Fault')
ylabel('Acceleration (g)')
annotation('doublearrow', [0.37 0.71], [0.8 0.8])
text(0.047, 20, ['0.009 sec \approx 1/BPFI, BPFI = ' num2str(dataInner.BPFI)])

변조된 진폭을 추출하기 위해 원시 신호의 포락선을 계산하고 아래쪽 서브플롯에 시각화합니다.

subplot(2, 1, 2)
[pEnvInner, fEnvInner, xEnvInner, tEnvInner] = envspectrum(xInner, fsInner);
plot(tEnvInner, xEnvInner)
xlim([0.04 0.06])
xlabel('Time (s)')
ylabel('Acceleration (g)')
title('Envelope signal')

이번에는 포락선 신호의 파워 스펙트럼을 계산하고 BPFI와 그 고조파에서 주파수 응답을 살펴봅니다.

figure
plot(fEnvInner, pEnvInner)
xlim([0 1000])
ncomb = 10;
helperPlotCombs(ncomb, dataInner.BPFI)
xlabel('Frequency (Hz)')
ylabel('Peak Amplitude')
title('Envelope Spectrum: Inner Race Fault')
legend('Envelope Spectrum', 'BPFI Harmonics')

대부분의 에너지가 BPFI와 그 고조파에 집중된 것을 볼 수 있습니다. 이는 베어링의 내륜에 결함이 있음을 나타내며, 데이터의 결함 유형과 일치합니다.

다른 결함 유형에 포락선 스펙트럼 분석 적용하기

이번에는 정상 데이터와 외륜 결함 데이터에 대해 동일한 포락선 스펙트럼 분석을 반복합니다.

dataNormal = load(fullfile(matlabroot, 'toolbox', 'predmaint', ...
    'predmaintdemos', 'bearingFaultDiagnosis', ...
    'train_data', 'baseline_1.mat'));
xNormal = dataNormal.bearing.gs;
fsNormal = dataNormal.bearing.sr;
tNormal = (0:length(xNormal)-1)/fsNormal;
[pEnvNormal, fEnvNormal] = envspectrum(xNormal, fsNormal);

figure
plot(fEnvNormal, pEnvNormal)
ncomb = 10;
helperPlotCombs(ncomb, [dataNormal.BPFO dataNormal.BPFI])
xlim([0 1000])
xlabel('Frequency (Hz)')
ylabel('Peak Amplitude')
title('Envelope Spectrum: Normal')
legend('Envelope Spectrum', 'BPFO Harmonics', 'BPFI Harmonics')

예상대로, 정상 베어링 신호의 포락선 스펙트럼은 BPFO나 BPFI에서 유의미한 피크를 보이지 않습니다.

dataOuter = load(fullfile(matlabroot, 'toolbox', 'predmaint', ...
    'predmaintdemos', 'bearingFaultDiagnosis', ...
    'train_data', 'OuterRaceFault_2.mat'));
xOuter = dataOuter.bearing.gs;
fsOuter = dataOuter.bearing.sr;
tOuter = (0:length(xOuter)-1)/fsOuter;
[pEnvOuter, fEnvOuter, xEnvOuter, tEnvOuter] = envspectrum(xOuter, fsOuter);

figure
plot(fEnvOuter, pEnvOuter)
ncomb = 10;
helperPlotCombs(ncomb, dataOuter.BPFO)
xlim([0 1000])
xlabel('Frequency (Hz)')
ylabel('Peak Amplitude')
title('Envelope Spectrum: Outer Race Fault')
legend('Envelope Spectrum', 'BPFO Harmonics')

외륜 결함 신호의 경우에도 BPFO 고조파에서 뚜렷한 피크가 보이지 않습니다. 포락선 스펙트럼 분석은 외륜 결함이 있는 베어링과 정상 베어링을 구분하지 못하는 것일까요? 잠시 멈추고 시간 영역의 신호를 여러 다른 상태에서 다시 살펴보겠습니다.

먼저 신호를 다시 시간 영역에서 시각화하고 첨도를 계산합니다. 첨도는 확률 변수의 4차 표준 모멘트입니다. 이것으로 신호의 임펄스 정도 또는 확률 변수의 꼬리 부분이 얼마나 두터운지를 특성화합니다.

figure
subplot(3, 1, 1)
kurtInner = kurtosis(xInner);
plot(tInner, xInner)
ylabel('Acceleration (g)')
title(['Inner Race Fault, kurtosis = ' num2str(kurtInner)])
xlim([0 0.1])

subplot(3, 1, 2)
kurtNormal = kurtosis(xNormal);
plot(tNormal, xNormal)
ylabel('Acceleration (g)')
title(['Normal, kurtosis = ' num2str(kurtNormal)])
xlim([0 0.1])

subplot(3, 1, 3)
kurtOuter = kurtosis(xOuter);
plot(tOuter, xOuter)
xlabel('Time (s)')
ylabel('Acceleration (g)')
title(['Outer Race Fault, kurtosis = ' num2str(kurtOuter)])
xlim([0 0.1])

내륜 결함 신호의 임펄스 정도가 유의미하게 큰 것을 알 수 있으며, 이로 인해 포락선 스펙트럼 분석은 BPFI에서 결함 시그니처를 효과적으로 캡처할 수 있습니다. 외륜 결함 신호의 경우에도 BPFO에서 진폭 변조를 어느 정도 식별할 수 있으나, 이는 강한 잡음에 의해 가려져 있습니다. 정상 신호는 어떠한 진폭 변조도 보이지 않습니다. 포락선 스펙트럼 분석에서는 BPFO에서 진폭 변조가 있는 임펄스 신호를 추출하는 작업 (또는 신호 대 잡음비를 향상하는 작업)이 주요 전처리 단계입니다. 다음 섹션에서는 가장 높은 첨도를 갖는 신호를 추출하는 커토그램 및 스펙트럼 첨도를 소개하고, 필터링된 신호에 포락선 스펙트럼 분석을 수행합니다.

대역 선택을 위한 커토그램 및 스펙트럼 첨도

커토그램 및 스펙트럼 첨도는 주파수 대역 내에서 국소적으로 첨도를 계산합니다. 이들은 가장 높은 첨도를 갖는(또는 가장 높은 신호 대 잡음비를 갖는) 주파수 대역을 찾는 데 사용할 수 있는 강력한 툴입니다 [2]. 가장 높은 첨도를 갖는 주파수 대역을 찾은 후에는 원시 신호에 대역통과 필터를 적용하여 포락선 스펙트럼 분석을 위해 더욱 임펄스 정도가 큰 신호를 얻을 수 있습니다.

level = 9;
figure
kurtogram(xOuter, fsOuter, level)

이 커토그램은 중심 주파수가 2.67kHz이고 대역폭이 0.763kHz인 주파수 대역에서 가장 높은 첨도가 2.71임을 나타냅니다.

이번에는 커토그램에서 제안하는 최적의 윈도우 길이를 사용하여 스펙트럼 첨도를 계산합니다.

figure
wc = 128;
pkurtosis(xOuter, fsOuter, wc)

스펙트로그램에 주파수 대역을 시각화하기 위해 스펙트로그램을 계산하고 그 옆에 스펙트럼 첨도를 표시합니다. 스펙트럼 첨도를 또 다른 방식으로 해석하자면, 스펙트럼 첨도 값이 높다는 것은 해당 주파수에서 전력 변동이 크다는 것을 나타내기 때문에 스펙트럼 첨도는 신호에서 비정상(Nonstationary) 성분을 찾을 때 유용한 툴입니다[3].

helperSpectrogramAndSpectralKurtosis(xOuter, fsOuter, level)

제안된 중심 주파수와 대역폭을 사용하여 신호에 대역통과 필터를 적용하면 첨도가 개선될 수 있으며 외륜 결함의 변조된 진폭을 찾을 수 있습니다.

[~, ~, ~, fc, ~, BW] = kurtogram(xOuter, fsOuter, level);

bpf = designfilt('bandpassfir', 'FilterOrder', 200, 'CutoffFrequency1', fc-BW/2, ...
    'CutoffFrequency2', fc+BW/2, 'SampleRate', fsOuter);
xOuterBpf = filter(bpf, xOuter);
[pEnvOuterBpf, fEnvOuterBpf, xEnvOuterBpf, tEnvBpfOuter] = envspectrum(xOuter, fsOuter, ...
    'FilterOrder', 200, 'Band', [fc-BW/2 fc+BW/2]);

figure
subplot(2, 1, 1)
plot(tOuter, xOuter, tEnvOuter, xEnvOuter)
ylabel('Acceleration (g)')
title(['Raw Signal: Outer Race Fault, kurtosis = ', num2str(kurtOuter)])
xlim([0 0.1])
legend('Signal', 'Envelope')

subplot(2, 1, 2)
kurtOuterBpf = kurtosis(xOuterBpf);
plot(tOuter, xOuterBpf, tEnvBpfOuter, xEnvOuterBpf)
ylabel('Acceleration (g)')
xlim([0 0.1])
xlabel('Time (s)')
title(['Bandpass Filtered Signal: Outer Race Fault, kurtosis = ', num2str(kurtOuterBpf)])
legend('Signal', 'Envelope')

대역통과 필터링 후에 첨도 값이 증가한 것을 볼 수 있습니다. 이번에는 주파수 영역에서 포락선 신호를 시각화합니다.

figure
plot(fEnvOuterBpf, pEnvOuterBpf);
ncomb = 10;
helperPlotCombs(ncomb, dataOuter.BPFO)
xlim([0 1000])
xlabel('Frequency (Hz)')
ylabel('Peak Amplitude')
title('Envelope Spectrum of Bandpass Filtered Signal: Outer Race Fault ')
legend('Envelope Spectrum', 'BPFO Harmonics')

커토그램과 스펙트럼 첨도에서 제안하는 주파수 대역으로 원시 신호에 대역통과 필터를 적용한 결과 포락선 스펙트럼 분석이 BPFO와 그 고조파에서 결함 시그니처를 드러낼 수 있음을 볼 수 있습니다.

일괄 처리

이번에는 파일 앙상블 데이터저장소를 사용하여 배치 단위의 훈련 데이터에 알고리즘을 적용해 보겠습니다.

이 툴박스에서는 데이터셋의 일부를 제공합니다. 데이터셋을 현재 폴더에 복사하고 쓰기 권한을 활성화합니다.

copyfile(...
    fullfile(matlabroot, 'toolbox', 'predmaint', 'predmaintdemos', ...
    'bearingFaultDiagnosis'), ...
    'RollingElementBearingFaultDiagnosis-Data-master')
fileattrib(fullfile('RollingElementBearingFaultDiagnosis-Data-master', 'train_data', '*.mat'), '+w')
fileattrib(fullfile('RollingElementBearingFaultDiagnosis-Data-master', 'test_data', '*.mat'), '+w')

전체 데이터셋을 사용하려면 https://github.com/mathworks/RollingElementBearingFaultDiagnosis-Data에서 전체 리포지토리를 zip 파일로 다운로드하고 라이브 스크립트와 동일한 디렉터리에 저장하십시오. 다음 명령을 사용하여 파일의 압축을 풉니다.

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

이 예제의 결과는 전체 데이터셋에서 생성한 것입니다. 전체 데이터셋은 14개의 mat 파일(정상 2개, 내륜 결함 4개, 외륜 결함 7개)을 갖는 훈련 데이터셋과 6개의 mat 파일(정상 1개, 내륜 결함 2개, 외륜 결함 3개)을 갖는 테스트 데이터셋을 포함합니다.

ReadFcn 과 WriteToMemberFcn에 함수 핸들을 할당하면 파일 앙상블 데이터저장소가 파일을 탐색하여 원하는 형식으로 데이터를 가져올 수 있습니다. 예를 들어, 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.WriteToMemberFcn = @writeMFPTBearing;
ensembleTrain.SelectedVariables = ["gs", "sr", "rate", "load", "BPFO", "BPFI", "FTF", "BSF", "Label", "FileName"]

ensembleTrain = 
  fileEnsembleDatastore with properties:

                 ReadFcn: @readMFPTBearing
        WriteToMemberFcn: @writeMFPTBearing
           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]

ensembleTrainTable = tall(ensembleTrain)

Starting parallel pool (parpool) using the 'local' profile ...
connected to 6 workers.

ensembleTrainTable =

  M×10 tall table

           gs             sr      rate    load     BPFO      BPFI      FTF       BSF           Label                   FileName        
    _________________    _____    ____    ____    ______    ______    ______    _____    __________________    ________________________

    [146484×1 double]    48828     25       0     81.125    118.88    14.838    63.91    "Inner Race Fault"    "InnerRaceFault_vload_1"
    [146484×1 double]    48828     25      50     81.125    118.88    14.838    63.91    "Inner Race Fault"    "InnerRaceFault_vload_2"
    [146484×1 double]    48828     25     100     81.125    118.88    14.838    63.91    "Inner Race Fault"    "InnerRaceFault_vload_3"
    [146484×1 double]    48828     25     150     81.125    118.88    14.838    63.91    "Inner Race Fault"    "InnerRaceFault_vload_4"
    [146484×1 double]    48828     25     200     81.125    118.88    14.838    63.91    "Inner Race Fault"    "InnerRaceFault_vload_5"
    [585936×1 double]    97656     25     270     81.125    118.88    14.838    63.91    "Outer Race Fault"    "OuterRaceFault_1"      
    [585936×1 double]    97656     25     270     81.125    118.88    14.838    63.91    "Outer Race Fault"    "OuterRaceFault_2"      
    [146484×1 double]    48828     25      25     81.125    118.88    14.838    63.91    "Outer Race Fault"    "OuterRaceFault_vload_1"
            :              :       :       :        :         :         :         :              :                        :
            :              :       :       :        :         :         :         :              :                        :

분석의 마지막 섹션에서, BPFO 및 BPFI에서의 대역통과 필터링된 포락선 스펙트럼 진폭이 베어링 결함 진단을 위한 두 개의 상태 지표임을 알 수 있습니다. 따라서 다음 단계는 모든 훈련 데이터에서 이 두 상태 지표를 추출하는 것입니다. 알고리즘을 더 견고하게 만들려면 BPFO와 BPFI 주변에 협소한 대역(대역폭 = $10 Δ f$ 로, 여기서 $Δ f$ 는 파워 스펙트럼의 주파수 분해능)을 설정하고 이 협소한 대역 내부에서 최대 진폭을 찾으십시오. 알고리즘은 아래에 나열된 bearingFeatureExtraction 함수에 포함되어 있습니다. 이 예제의 나머지 부분에서는 BPFI와 BPFO 주변의 포락선 스펙트럼 진폭을 "BPFIAmplitude"와 "BPFOAmplitude"라고 칭합니다.

% To process the data in parallel, use the following code
% ppool = gcp;
% n = numpartitions(ensembleTrain, ppool);
% parfor ct = 1:n
%     subEnsembleTrain = partition(ensembleTrain, n, ct);
%     reset(subEnsembleTrain);
%     while hasdata(subEnsembleTrain)
%         bearingFeatureExtraction(subEnsembleTrain);
%     end
% end
ensembleTrain.DataVariables = [ensembleTrain.DataVariables; "BPFIAmplitude"; "BPFOAmplitude"];
reset(ensembleTrain)
while hasdata(ensembleTrain)
    bearingFeatureExtraction(ensembleTrain)
end

파일 앙상블 데이터저장소에 새로운 상태 지표를 추가한 후에는 SelectedVariables가 파일 앙상블 데이터저장소에서 관련 데이터를 읽어 들이도록 지정하고 추출된 상태 지표를 포함하는 특징 테이블을 만듭니다.

ensembleTrain.SelectedVariables = ["BPFIAmplitude", "BPFOAmplitude", "Label"];
featureTableTrain = tall(ensembleTrain);
featureTableTrain = gather(featureTableTrain);

Evaluating tall expression using the Parallel Pool 'local':
- Pass 1 of 1: 0% complete
Evaluation 0% complete

- Pass 1 of 1: Completed in 3 sec
Evaluation completed in 3 sec

featureTableTrain

featureTableTrain=14×3 table
    BPFIAmplitude    BPFOAmplitude          Label       
    _____________    _____________    __________________

        0.33918         0.082296      "Inner Race Fault"
        0.31488         0.026599      "Inner Race Fault"
        0.52356         0.036609      "Inner Race Fault"
        0.52899         0.028381      "Inner Race Fault"
        0.13515         0.012337      "Inner Race Fault"
       0.004024          0.03574      "Outer Race Fault"
      0.0044918           0.1835      "Outer Race Fault"
      0.0074993          0.30166      "Outer Race Fault"
       0.013662          0.12468      "Outer Race Fault"
      0.0070963          0.28215      "Outer Race Fault"
      0.0060772          0.35241      "Outer Race Fault"
       0.011244          0.17975      "Outer Race Fault"
      0.0036798        0.0050208      "Normal"          
        0.00359        0.0069449      "Normal"

이렇게 만든 특징 테이블을 시각화합니다.

figure
gscatter(featureTableTrain.BPFIAmplitude, featureTableTrain.BPFOAmplitude, featureTableTrain.Label, [], 'ox+')
xlabel('BPFI Amplitude')
ylabel('BPFO Amplitude')

BPFI 진폭과 BPFO 진폭의 상대 값은 서로 다른 유형의 결함을 효과적으로 나타낼 수 있습니다. 두 개의 기존 특징에 로그 비율을 적용한 특징을 새로 생성하고, 각 결함 유형별로 그룹화한 히스토그램으로 시각화합니다.

featureTableTrain.IOLogRatio = log(featureTableTrain.BPFIAmplitude./featureTableTrain.BPFOAmplitude);
figure
hold on
histogram(featureTableTrain.IOLogRatio(featureTableTrain.Label=="Inner Race Fault"), 'BinWidth', 0.5)
histogram(featureTableTrain.IOLogRatio(featureTableTrain.Label=="Outer Race Fault"), 'BinWidth', 0.5)
histogram(featureTableTrain.IOLogRatio(featureTableTrain.Label=="Normal"), 'BinWidth', 0.5)
plot([-1.5 -1.5 NaN 0.5 0.5], [0 3 NaN 0 3], 'k--')
hold off
ylabel('Count')
xlabel('log(BPFIAmplitude/BPFOAmplitude)')
legend('Inner Race Fault', 'Outer Race Fault', 'Normal', 'Classification Boundary')

히스토그램은 세 개의 서로 다른 베어링 상태가 확실히 구분된 것을 보여줍니다. BPFI 진폭과 BPFO 진폭의 로그 비율은 베어링 결함을 분류하는 데 유효한 특징입니다. 예제를 단순화하기 위해 매우 간단한 분류기를 도출합니다. 즉, $\log (\frac{BPFIAmp litude}{BPFOAmp litude}) \leq - 1.5$ 이면 베어링에 외륜 결함이 있고, $- 1.5 < \log (\frac{BPFIAmp litude}{BPFOAmp litude}) \leq 0.5$ 이면 베어링이 정상이고, $\log (\frac{BPFIAmp litude}{BPFOAmp litude}) > 0.5$ 이면 베어링에 내륜 결함이 있습니다.

테스트 데이터 세트를 사용한 검증

이제 테스트 데이터 세트에 워크플로를 적용하여, 직전 섹션에서 얻은 분류기를 검증해 보겠습니다. 테스트 데이터는 1개의 정상 데이터 세트, 2개의 내륜 결함 데이터 세트, 3개의 외륜 결함 데이터 세트를 포함합니다.

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.WriteToMemberFcn = @writeMFPTBearing;
ensembleTest.SelectedVariables = ["gs", "sr", "rate", "load", "BPFO", "BPFI", "FTF", "BSF", "Label", "FileName"]

ensembleTest = 
  fileEnsembleDatastore with properties:

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

ensembleTest.DataVariables = [ensembleTest.DataVariables; "BPFIAmplitude"; "BPFOAmplitude"];
reset(ensembleTest)
while hasdata(ensembleTest)
    bearingFeatureExtraction(ensembleTest)
end

ensembleTest.SelectedVariables = ["BPFIAmplitude", "BPFOAmplitude", "Label"];
featureTableTest = tall(ensembleTest);
featureTableTest = gather(featureTableTest);

Evaluating tall expression using the Parallel Pool 'local':
- Pass 1 of 1: Completed in 1 sec
Evaluation completed in 1 sec

featureTableTest.IOLogRatio = log(featureTableTest.BPFIAmplitude./featureTableTest.BPFOAmplitude);

figure
hold on
histogram(featureTableTrain.IOLogRatio(featureTableTrain.Label=="Inner Race Fault"), 'BinWidth', 0.5)
histogram(featureTableTrain.IOLogRatio(featureTableTrain.Label=="Outer Race Fault"), 'BinWidth', 0.5)
histogram(featureTableTrain.IOLogRatio(featureTableTrain.Label=="Normal"), 'BinWidth', 0.5)

histogram(featureTableTest.IOLogRatio(featureTableTest.Label=="Inner Race Fault"), 'BinWidth', 0.1)
histogram(featureTableTest.IOLogRatio(featureTableTest.Label=="Outer Race Fault"), 'BinWidth', 0.1)
histogram(featureTableTest.IOLogRatio(featureTableTest.Label=="Normal"), 'BinWidth', 0.1)
plot([-1.5 -1.5 NaN 0.5 0.5], [0 3 NaN 0 3], 'k--')
hold off
ylabel('Count')
xlabel('log(BPFIAmplitude/BPFOAmplitude)')
legend('Inner Race Fault - Train', 'Outer Race Fault - Train', 'Normal - Train', ...
    'Inner Race Fault - Test', 'Outer Race Fault - Test', 'Normal - Test', ...
    'Classification Boundary')

테스트 데이터 세트의 BPFI 진폭과 BPFO 진폭의 로그 비율은 훈련 데이터 세트의 로그 비율과 일치하는 분포를 보입니다. 직전 섹션에서 얻은 나이브 분류기는 테스트 데이터 세트에서 완벽한 정확도를 달성했습니다.

잘 일반화되는 분류기를 얻으려면 대개 하나의 특징만으로는 충분하지 않다는 사실에 유의하십시오. 데이터를 여러 개의 조각으로 나누고(데이터 점을 더 많이 만들기 위해), 진단 관련 특징을 여러 개 추출하고, 특징의 중요도 순위에 따라 특징의 서브셋을 선택한 다음 Statistics & Machine Learning Toolbox의 분류 학습기 앱을 사용하여 다양한 분류기를 훈련시키면 보다 정교한 분류기를 얻을 수 있습니다. 이 워크플로의 자세한 내용을 보려면 "Using Simulink to generate fault data" 예제를 참조하십시오.

요약

이 예제에서는 커토그램, 스펙트럼 첨도 및 포락선 스펙트럼을 사용하여 구름 요소 베어링에 있는 여러 유형의 결함을 식별하는 방법을 보여주었습니다. 알고리즘을 디스크에 있는 한 배치의 데이터 세트에 적용한 결과 BPFI와 BPFO의 대역통과 필터링된 포락선 스펙트럼의 진폭이 베어링 진단을 위한 두 개의 중요한 상태 지표임을 알 수 있었습니다.

참고 문헌

[1] Randall, Robert B., and Jerome Antoni. "Rolling element bearing diagnostics—a tutorial." Mechanical Systems and Signal Processing. Vol. 25, Number 2, 2011, pp. 485–520.

[2] Antoni, Jérôme. "Fast computation of the kurtogram for the detection of transient faults." Mechanical Systems and Signal Processing. Vol. 21, Number 1, 2007, pp. 108–124.

[3] Antoni, Jérôme. "The spectral kurtosis: a useful tool for characterising non-stationary signals." Mechanical Systems and Signal Processing. Vol. 20, Number 2, 2006, pp. 282–307.

[4] Bechhoefer, Eric. "Condition Based Maintenance Fault Database for Testing Diagnostics and Prognostic Algorithms." 2013. https://www.mfpt.org/fault-data-sets/

헬퍼 함수

function bearingFeatureExtraction(ensemble)
% Extract condition indicators from bearing data
data = read(ensemble);
x = data.gs{1};
fs = data.sr;

% Critical Frequencies
BPFO = data.BPFO;
BPFI = data.BPFI;

level = 9;
[~, ~, ~, fc, ~, BW] = kurtogram(x, fs, level);

% Bandpass filtered Envelope Spectrum
[pEnvpBpf, fEnvBpf] = envspectrum(x, fs, 'FilterOrder', 200, 'Band', [max([fc-BW/2 0]) min([fc+BW/2 0.999*fs/2])]);
deltaf = fEnvBpf(2) - fEnvBpf(1);

BPFIAmplitude = max(pEnvpBpf((fEnvBpf > (BPFI-5*deltaf)) & (fEnvBpf < (BPFI+5*deltaf))));
BPFOAmplitude = max(pEnvpBpf((fEnvBpf > (BPFO-5*deltaf)) & (fEnvBpf < (BPFO+5*deltaf))));

writeToLastMemberRead(ensemble, table(BPFIAmplitude, BPFOAmplitude, 'VariableNames', {'BPFIAmplitude', 'BPFOAmplitude'}));
end

참고 항목

fileEnsembleDatastore