ベアリング診断の「精度 99%」はどこから来るのか ―― 定番データセットで確かめる

機械に後付けして、振動と電流から故障の兆候を知らせる小さな装置を作りたいと考えています。業種を問わず、現場で動いている機械にはほぼ例外なくモータとベアリングが入っていて、ベアリングが壊れればラインが止まります。

ただ、手元にまだ実機がありません。壊れかけのベアリングも、実験に使える機械もない。そこで先に、公開されている実測データで、振動からどこまで故障が見えるのかを確かめておくことにしました。

以前、百円マイコンにマハラノビス距離を載せて、ベアリングの劣化を「いつもと違う」と検出する記事を書きました。前回は異常の検出まで。今回はその次の段階、どこが壊れたかの特定を試します。記事の後半は、この分野の論文や解説でよく見かける「診断精度 99%」という数字の中身を、自分の手で確かめる話になります。

定番データセット CWRU とそのクセ

使ったのは、ベアリング故障診断のベンチマークとしていちばん名の通った CWRU のデータです。2 馬力のモータのベアリングに放電加工で人工キズを入れたもので、場所は内輪・転動体（玉）・外輪の 3 種、キズの直径は数段階。これをモータ負荷 0〜3 HP で回して振動を収録してあります。今回はドライブエンド 12 kHz 収録のサブセットから、キズ径 0.007 / 0.014 / 0.021 インチの 3 段階、正常 4 本と故障 36 本の計 40 レコードを使いました。

ダウンロードは .mat ファイルの直リンクです。ただ、20 年以上使われてきたデータセットだけあって、いくつかクセがあります。

• あるファイルには前の計測の変数が混入しています（99.mat の中に X098_DE_time が同居している）。「_DE_time で終わる最初の変数を拾う」ような素朴なパーサだと、隣のレコードを黙って読みます。
• 回転数の変数が入っていないファイルがあります。実測回転数も公称値と数 rpm ずれている（公称 1730 rpm に対し実測 1721 rpm など）ので、周波数の照合には実測値を使います。
• 正常データのサンプリング周波数が公式ページに明記されていません。データから推定したのですが、ここで一度判定を誤りました。結論を先に書くと正常は 48 kHz で、この記事では故障側と揃えるため 12 kHz に間引いています。

三つ目は失敗談として書いておきます。最初、12 kHz と仮定してスペクトルを描いたら、実測回転数どおりの位置（軸 1 次 29.9 Hz）にピークが立ったので、正常も 12 kHz だと判定しました。ところがこのモータは 4 極で、同期 1800 rpm のすぐ下で回っています。48 kHz のデータを 12 kHz として読むと周波数軸が 4 分の 1 に縮むため、実 120 Hz、つまり電源周波数の 2 倍に出る電磁加振の線が、30.0 Hz の位置に現れます。軸 1 次の 29.9 Hz とほぼ同じ位置で、分解能では区別できません。

判定の決め手になったのは、回転が下がる負荷 3 HP のレコード（軸 28.75 Hz）でした。軸回転のピークだと思っていた線は、このレコードでも 30.0 Hz から動きません。回転に追従しない、電気由来の線だったわけです。本物の 12 kHz である故障レコードなら、同じ負荷でピークはそれぞれの実測回転数（28.7〜28.8 Hz）に追従します。

決定的に見えた根拠が偶然の一致だったという構図は、後半の評価の話でもう一度出てきます。

キズの場所ごとに、欠陥周波数は計算で決まる

ベアリング診断で使われるのは、キズの場所ごとに振動へ出る周波数が幾何から計算できるという性質です。玉が外輪の一点のキズを叩く頻度、内輪のキズが玉の列に当たる頻度、保持器の公転、玉の自転。どれも玉の数と径、ピッチ円径、接触角、回転数から決まります。

このモータのベアリングについて、公開されている欠陥周波数（軸回転を 1 とした次数）は次のとおりです。

鵜呑みにせず、幾何の式から検算してみました。BPFO の式を玉数について逆に解くと 9.000 個と整数で出て、その玉数から 4 つの周波数がすべて小数 4 桁まで再現できました。ここでひとつ注意点があります。公式表の転動体の値 4.7135 次は、玉の自転周波数（2.3567 次）のちょうど 2 倍です。玉のキズは 1 自転のあいだに内輪と外輪の両方に 1 回ずつ当たるため、照合対象は自転の 2 倍になるわけです（実際にはこれに保持器 FTF の変調が側帯波として付くともされます）。資料によって BSF がどちらを指すか揺れる、昔からの紛らわしいポイントです。

エンベロープ解析の手順

キズの衝撃はごく短いパルスです。その繰り返し周波数は生のスペクトルでは見えにくく、代わりに数 kHz の構造共振を励振します。そこで定石どおり、共振帯だけを帯域通過で残し、ヒルベルト変換で衝撃の繰り返し（包絡）を取り出し、その包絡のスペクトルを描きます。

def envelope_spectrum(x, fs=12_000):
    b, a = butter(4, (2500, 5500), btype="bandpass", fs=fs)  # 共振帯だけ残す
    env = np.abs(hilbert(filtfilt(b, a, x)))                 # 衝撃の「繰り返し」を包絡で
    env -= env.mean()
    return welch(env, fs=fs, nperseg=1 << 15)                # 繰り返しの周波数を見る

判定は機械学習なしです。計算で求めた位置（各欠陥周波数の 1〜3 次）にピークが立っているかを、ピークが近傍のノイズの床からどれだけ dB で浮いているかで点数化し、最も高い場所を答えとします。

内輪キズは赤の破線（BPFI）の位置に、外輪キズは青の破線（BPFO）の位置に、高い線が立っています。回転数と寸法から事前に計算した位置と、ぴったり重なりました。正常にも山はありますが（このレコードで最も高いのは軸回転のほぼ 4 倍の位置でした）、こちらはどの欠陥周波数とも合いません。

一方、転動体のキズだけは、緑の破線（2×BSF）に何も立っていません。

40 レコードの判定結果

故障 36 レコードをこの方法で判定した結果です。スコアは負荷 4 条件の平均で、正常 4 本は 3 部位のスコアがどれも一桁 dB に収まるため、しきい値を設ければ除外できます。

内輪   : 12/12  （全サイズで BPFI がくっきり）
外輪   :  8/12  （0.007 と 0.021 は完璧、0.014 だけ 0/4）
転動体 :  2/12  （当たった 2 本もスコアはノイズ床レベル）

転動体の判定がほぼ全部外れたのは帯域の選び方の問題かと思い、共振帯の仮定を 6 通りに変えて試しましたが、どの帯域でも 2×BSF の線は立ちませんでした。これはこの分野では知られた結果で、CWRU データを系統的に検証した Smith と Randall のベンチマーク研究（2015）も、転動体キズの多くと外輪 0.014 インチは古典的なエンベロープ解析では診断が難しいと報告しています。手元の結果も同じ傾向でした。物理的な理由としては、玉のキズは自転しながら保持器で運ばれてレースに当たる向きが変わり続けるため、衝撃の間隔が一定になりにくいからだと言われています。

この時点では、幾何から計算した周波数だけで外輪と内輪はかなり特定できる一方、転動体はこのデータをこの方法で見る限りほぼ特定できない、という結果です。

「診断精度 99%」はどこから来るのか

この分野の論文やブログでは、CWRU で 99% 台の診断精度をよく見かけます。定番のタスク設定は、正常と 3 部位 × 3 サイズを合わせた 10 クラス分類です。信号を 1 秒の窓に切り（重なりなし、各レコード最大 10 枚で計 395 窓）、窓単位でランダムに訓練とテストへ分けます。

再現してみました。特徴量はスペクトルを 16 帯域に割った log パワー、分類器はロジスティック回帰で、深層学習は使いません。

結果は 99.7% でした（395 窓。窓を 5 組に分けて順にテストへ回す、いわゆる 5-fold です）。

エンベロープ解析で「転動体は見えない」と確認したばかりの、同じデータです。要因は評価の切り方にあります。窓をランダムに分けると、同じ収録の隣どうしの窓が訓練とテストの両側に混ざります。こうなると分類器は、故障の物理をとらえなくても、収録ごとの固有のパターンを暗記するだけで高い精度を出せてしまいます。

さらに根本的なのは、CWRU の故障クラスがどれも物理的に 1 個のベアリング、つまり 1 つの加工キズ個体だという点です。負荷 0〜3 HP といっても、同じキズ個体を条件だけ変えて収録したものにすぎません。そのため「訓練は負荷 0〜2、テストは負荷 3」という一見問題のない分け方でも、テスト側にいるのは訓練済みの同じ個体です。試しに RMS と尖度など時間統計 4 個だけで評価したところ、この負荷またぎでもあっさり 100% が出ました。

別の個体でテストするには、このデータではキズサイズをまたぐしかありません。0.007 と 0.021 で訓練して 0.014 でテストする、という分け方です。サイズが違えば別の個体、別の加工です。ただしサイズをまたぐ場合、10 クラスのままでは評価が組めません。たとえば 0.014 をテストに出すと、「内輪 0.014」のようなクラスを訓練側が一度も見ていないことになるからです。そこでラベルを部位 3 クラス（内輪・転動体・外輪）に落とし、ここからは故障レコードだけの 3 クラス分類として、同じ特徴量・同じ分類器で、ランダム分割とサイズまたぎを比較します。

収録のパターンを拾いやすい帯域エネルギー特徴は、100% から 68.1% に下がります。内訳を見ると、テストサイズが 0.007 のときは 100%、0.014 のときは 33.3% で、後者は 3 クラスのあてずっぽうと同じです。精度が高く出るか低く出るかすら、未知の個体が来るまで分かりません。時間統計は全体で 41.7% でした。前半のサンプリング周波数の誤判定と同じ構図で、強く見える結果ほど、別の説明と切り分けてから信じる必要があります。

一方、欠陥周波数のスコアを特徴にした物理特徴は 68.6% でした。数字としては高くありません。ただしテストサイズを変えても 57.5〜76.7% の範囲に収まり、誤りは転動体と外輪 0.014 に集中します。エンベロープ解析で特定できないと分かった場所と同じです。

物理特徴の価値は数字の高さではなく崩れ方に出る

両者を比べると差がはっきりします。収録のパターンを拾う特徴は、リークのある評価では高い数字が出ますが、未知の個体では下がり、どこで下がるかも事前に分かりません。物理に根ざした特徴は、数字は控えめでも、外す場所がほぼ決まっていて（転動体と外輪 0.014）、その理由に物理の説明がつきます。

現場に置く装置にとってどちらが重要かと言えば、私は後者だと考えています。必要なのは見かけの高い精度ではなく、何が検出できて何ができないかの線引きです。外輪と内輪の異常はこの方法で拾える、転動体は振動だけでは難しいので別の情報源を足すか、検出（いつもと違う、の判定）までで割り切る。そうした設計判断は、リークのない評価からしか得られません。

省いていること

• CWRU のキズは放電加工の人工キズです。実際の疲労剥離は形も進行も違います（実損傷を含むデータセットもあり、次回扱います）。
• 回転数ほぼ一定・単一故障・センサはハウジング直付きという、診断には恵まれた条件です。現場では回転数が変わり、故障は重なり、センサは付けたい場所に付けられるとは限りません。
• 転動体レコードの包絡には、2×BSF の代わりに保持器（FTF）系の線が見えるという観察はありました。ただ手元の転動体は各サイズ 4 レコードしかなく、特徴に加えると 4 本への過適合になりかねないため、今回は使っていません。
• 正常データのサンプリング周波数の判定は、回転が同期速度に近い負荷 0〜1 のレコードでは原理的にできません（軸 1 次と電源 2 倍の線がほぼ同じビンに落ちるため）。今回の判定は、回転が十分下がる負荷 3 HP のレコードと、対照の故障レコードに依っています。
• 評価は 1 秒窓・重なりなし・計 400 窓弱の小規模なものです。数字は傾向を見るための目安で、細かい差を論じられる規模ではありません。

おわりに

幾何から計算した周波数だけで、外輪と内輪のキズは特定できました。機械学習は使っていません。一方の「99%」は、少なくともこのデータでは、故障の物理ではなく収録ごとのパターンを学習した結果でした。同じデータ、同じ分類器でも、評価の切り方で数字は別物になります。作る装置に載せるべきものも見えてきました。見かけの精度ではなく、性能の落ち方が緩やかな物理特徴と、リークのない評価です。

次回は電流を足します。振動と電流を同期収録した公開データセット（Paderborn 大学）で、センサを増やすと検出の範囲がどう広がるかを確かめる予定です。

ダウンロードスクリプト（出典・ファイル番号の対応表つき）と、検証・診断・評価のコードは、そのまま動く形で公開しています。

• リポジトリ: logicia32/bearing-diag-lab
• データ出典: Case Western Reserve University Bearing Data Center
• 参考文献: W. A. Smith and R. B. Randall, "Rolling element bearing diagnostics using the Case Western Reserve University data: a benchmark study," Mechanical Systems and Signal Processing, 2015