Main Content

転動体ベアリングの故障診断

この例では、特に他のマシン コンポーネントからの強いマスキング信号がある場合に、加速度信号に基づいて転動体ベアリングの故障診断を行う方法を説明します。ここでは包絡線スペクトル解析とスペクトル尖度を適用してベアリング故障を診断する方法を示します。この方法はビッグ データ アプリケーションにスケール アップできます。

問題の概要

転動体ベアリングにおける局所的な故障は、外輪、内輪、ケージ、または転動体に発生する可能性があります。転動体が外輪または内輪上の局所的な故障部位に当たったり、転動体上の故障が外輪や内輪に当たると、ベアリングと応答変換器間に高周波共振が励起されます [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=Dfr2d[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 は、2.67 kHz を中心として 0.763 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 つ) からなるテスト データセットが含まれています。

ReadFcn WriteToMemberFcn に関数ハンドルを割り当てることで、ファイル アンサンブルのデータストアがファイル内に移動してデータを必要な形式で取得できるようにします。たとえば、MFPT データには構造体 bearing があり、これは振動信号 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 and Machine Learning Toolbox の分類学習器アプリを使用してさまざまな分類器に学習させることで、より洗練された分類器を得ることができます。このワークフローの詳細については、「Simulink を使用した故障データの生成」の例を参照してください。

まとめ

この例では、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://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

参考

関連するトピック