ドキュメンテーション

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

回転要素のベアリングの故障診断

この例は、特に他のマシン コンポーネントから強いマスキング信号が発生している場合に加速度信号に基づく回転要素のベアリングの故障診断を実行する方法を説明します。この例では、包絡線スペクトル解析とスペクトル尖度を適用してベアリングの故障を診断する方法と、それがビッグ データ事例に拡張可能であることを示します。

問題概要

回転要素のベアリングの局所化した故障は、外輪、内輪、ケージ、または回転要素で発生する可能性があります。ベアリングと応答変換器間の高周波数の応答は、回転要素が外輪と内輪上で局所的な故障を発生させたとき、または回転要素での故障が外輪または内輪に発生したときに励起されます [1]。次の図は、内輪の局所的故障を発生させている回転要素を示します。問題は、さまざまなタイプの故障をどのように検出して特定するかです。

Machinery Failure Prevention Technology (MFPT) Challenge データ

MFPT Challenge データ [4] には、さまざまな故障条件の下でマシンから収集された 23 個のデータセットが含まれます。最初の 20 個のデータセットはベアリング テスト装置から収集され、良好な状態が 3 個、定数負荷での外輪の故障が 3 個、さまざまな負荷での外輪の故障が 7 個、さまざまな負荷での内輪の故障が 7 個あります。残りの 3 個のデータセット (オイル ポンプのベアリング、中間速度のベアリング、およびプラネット ベアリング) は、実際のマシンから収集されます。故障の場所は不明です。この例では、既知の条件でテスト装置から収集されたデータのみが使用されます。

各データセットには、加速度信号 "gs"、サンプリング レート "sr"、シャフト速度 "rate"、負荷の重み "load"、およびさまざまな故障場所を表す 4 つの臨界周波数 (外輪転動体通過周波数 (BPFO)、内輪転動体通過周波数 (BPFI)、保持器周波数 (FTF)、転動体回転周波数 (BSF)) が含まれます。ここに、臨界周波数の式を示します [1]。

  • 転動体通過周波数、外輪 (BPFO)

BPFO=nfr2(1-dDcosϕ)

  • 転動体通過周波数、内輪 (BPFI)

BPFI=nfr2(1+dDcosϕ)

  • 保持器周波数 (FTF)、ケージ速度とも呼ばれる

FTF=fr2(1-dDcosϕ)

  • 転動体 (ローラー) 回転周波数

BSF=D2d[1-(dDcosϕ)2]

図に示すように、d はボールの直径、D はピッチの直径です。変数 fr はシャフト速度、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.009 Hz 111 Hz あたりであると観測されます。回転要素が内輪で局所的故障を発生する周波数、つまり BPFI が 118.875 Hz であることがわかっています。これは、ベアリングに潜在的に内輪の故障があることを示しています。

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 で振幅変調を使用したインパルス信号の抽出 (S/N 比の強調) は、包絡線スペクトル解析前に行う主要な前処理手順です。次のセクションでは、kurtogram とスペクトル尖度を導入し、最高尖度を含む信号を抽出し、フィルター処理された信号の包絡線スペクトル解析を実行します。

帯域選択の Kurtogram とスペクトル尖度

Kurtogram とスペクトル尖度は周波数帯域内で局所的に尖度を計算します。最高尖度 (あるいは最高の S/N 比) を持つ周波数帯域を検出するための強力なツールがあります [2]。最高尖度を持つ周波数帯域を正確に特定した後、バンドパス フィルターは生の信号に適用され、包絡線スペクトル解析のインパルス信号をさらに取得できます。

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

kurtogram は、0.763 kHz の帯域幅で、中心が 2.67 kHz の周波数帯域に 2.71 の最高尖度があることを示します。

kurtogram によって提案される最適なウィンドウの長さを使用して、スペクトル尖度を計算します。

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

スペクトログラムに周波数帯域を視覚化するために、スペクトル尖度を中心から離れた所に置きます。別の方法でスペクトル尖度を解釈するには、高いスペクトル尖度の値は対応する周波数でパワーの高い分散を示します。これにより、スペクトル尖度は信号の非定常成分を検出するための有用なツールになります [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')

kurtogram とスペクトル尖度によって提案される周波数帯域を使用して生の信号をフィルター処理するバンドパスによって、包絡線スペクトル解析が 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 個) のテスト データセットが含まれます。

ReadFcnWriteToMemberFcn に関数ハンドルを割り当てることによって、ファイル アンサンブル データストアは、ファイルに移動して、目的の形式にデータを取得できます。たとえば、MFPT データには、バイブレーション信号 gs、サンプリング レート sr などを保存する構造体 bearing があります。ベアリング構造自体を返す代わりに、ファイル アンサンブル データストアが bearing データ構造内部のバイブレーション信号 gs を返すために readMFPTBearing 関数が書き込まれます。

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 で包絡線スペクトルの振幅をフィルター処理したバンドパスがベアリング故障診断の 2 つの条件インジケーターであることがわかります。そのため、次の手順は、すべての学習データからの 2 つの条件インジケーターの抽出です。アルゴリズムをさらにロバストにするには、BPFO と BPFI 周辺の狭い帯域 (帯域幅 = 10Δf。ここで、Δf はパワー スペクトルの周波数分解能です) を設定した後、この狭い帯域幅内部の最大振幅を検出します。アルゴリズムは、下に示す関数 bearingFeatureExtraction に含まれます。BPFI と BPFO の周辺の包絡線スペクトルの振幅は、残りの例では "BPFI 振幅" と "BPFO 振幅" と呼ばれます。

% 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 振幅の相対値は、異なる故障のタイプの有効なインジケーターとなる可能性があります。ここで、新しい特徴が作成されます。これは、2 つの既存の特徴の対数比で、異なる故障のタイプによってグループ化されるヒストグラムで可視化されます。

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')

ヒストグラムは、3 つの異なるベアリング条件間で明確な間隔を示します。BPFI 振幅と BPFO 振幅間の対数比は、ベアリング故障を分類するための有効な特徴です。例を単純化するには、非常に簡単な分類器が次のようには得られます。log(BPFIAmplitudeBPFOAmplitude)-1.5 の場合、ベアリングに外輪の故障があります。-1.5<log(BPFIAmplitudeBPFOAmplitude)0.5 の場合、ベアリングは正常です。log(BPFIAmplitudeBPFOAmplitude)>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」の例を参照してください。

まとめ

この例は、kurtogram、スペクトル尖度および包絡線スペクトルを使用して、回転要素のベアリングで故障の異なるタイプを特定する方法を示します。アルゴリズムは、ディスクのデータセットのバッチに適用され、BPFI と BPFO で包絡線スペクトルをフィルター処理したバンドパスの振幅がベアリング診断の 2 つの重要な条件インジケーターであることを示すのに役立ちます。

参照

[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://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