ρ=0.998 の壁 — 百円台の 8-bit マイコンに 2D マハラノビス距離を乗せて、ベアリング劣化を聴く

「TinyML が来る」と言われて、もう数年たちます。STM32 H7 とか ESP32 とか、Cortex-M 系のけっこう立派なマイコンに、CNN や決定木を載せる事例は、わりと見るようになりました。

ただ、現場にはもっと安いマイコンがたくさん転がっています。製造ラインの古いコンベアの片隅、ビル設備の天井裏、農業ハウスのファンの根元。量産ベースで 1 個百円かそこら、現役の 8-bit マイコンが、何十年も働き続けている。これらに「ちょっとした目」を足したい、というニーズは、たぶん相当ある。

この記事はその素振りです。振動の 2 軸データから、ベアリングの劣化を「分布から外れた」と判定する処理を、PIC16F1503（14-pin DIP、SRAM 128 バイト、フラッシュ 2K words、ハードウェア乗算器なし）の上に乗せられるかを、NASA 公開のベアリング run-to-failure データで試してみました。

結論を先に言うと、乗ります。1 評価あたりおよそ 2,800 サイクル、フラッシュ消費 11%、SRAM 消費 21%。ただしそこに至るには、健常時の共分散行列が ρ_xy ≈ 0.998 でほとんど特異という壁に、一回正面から向き合う必要がありました。今日はその壁の話を書きます。

何のデータで遊んでいるか — NASA IMS Bearing Data

NASA の Prognostics Center of Excellence は、機械の故障予知研究用に、何種類かのデータセットを公開しています。今回使う IMS Bearing Data（University of Cincinnati 提供）は、ベアリングをわざと寿命が尽きるまで回し続けた実験データで、データセット界隈ではかなり古典の部類です。

Set 1（今回使うやつ）の中身は:

• 4 個のベアリング × 各 2 軸（x, y）の加速度 = 8 ch
• サンプリング 20 kHz、約 1 秒（20,480 点）のスナップショットをおおむね 10 分ごと（最初の 43 ファイルだけ 5 分間隔）
• 2003-10-22 から 2003-11-25 まで、約 35 日連続、合計 2,156 ファイル
• 期末に Bearing 3 で内輪レース欠陥（inner race defect）、Bearing 4 でローラー欠陥

今回は Bearing 3（Ch 5, 6）の 2 軸を主役にします。なぜかと言えば「inner race defect」はベアリング異常検知の教科書例で、構造振動として比較的綺麗に出てくる種類だから。読者の頭の中で「これは古典の罠案件だ」と思ってもらえるネタが、最初の一本としては都合がいい。

何を計算するか — 2D マハラノビス距離

「異常」の定義は突き詰めるとややこしいので、ここではいちばん素朴な統計距離を採ります。マハラノビス距離 $d^2$ は、観測点 $\mathbf{x}$ がベースライン分布の中心 $\boldsymbol{\mu}$ からどれだけ「分布の形を考慮して」離れているか、を測る量です。

2 次元（x, y 軸の RMS）なら、$\Sigma$ は対称 2×2 行列で、3 つの独立成分しかありません。式を展開すると

で済みます。割り算は 1 個だけ、しかも学習時に係数を事前計算してしまえば実行時にはゼロにできます。平方根も要りません（しきい値も $d^2$ で持てばいい）。これは 8-bit マイコンに乗せるにはありがたい性質です。

入力は、生波形そのものではなくウィンドウごとの特徴量 (RMS_x, RMS_y) に置きました。20 kHz の生サンプルに $d^2$ を毎回当てるのは数学的にも実装的にも無駄が多く、ベアリング監視の文脈では RMS は事実上業界の共通言語に近い。ここを変に凝っても自己満足になるので、定石どおりです。

16-bit に乗せようとして詰まった話 — ρ=0.998 の壁

健常 215 窓（実験開始からの最初のひとかたまり）から $\boldsymbol{\mu}$ と $\Sigma$ を学習すると、こうなりました。

mu  = (0.139, 0.142) g
Σ   = [[1.53e-4, 1.53e-4],
       [1.53e-4, 1.53e-4]]
ρ_xy = 0.998
det Σ = 9.98e-11

$\rho_{xy} = 0.998$。ベアリングの x 軸と y 軸の 窓ごとの RMS が、健常時にはほぼ同じ向きに動いていることを意味しています。一つの軸回転の励振を 2 方向から拾っているのだから、強さの増減はそろって出る、というのは物理的にも自然な見え方です（ただし瞬時の波形まで完全に同期、という強い意味ではありません）。

ところがここで $\Sigma^{-1}$ を計算すると、各要素が 1.5×10⁶ オーダーに跳ね上がります。det Σ が 10⁻¹⁰ レベルで小さく、行列がほぼ特異だからです。多重共線性で逆行列がノイズを過剰に増幅する、典型的な数値悪条件です。Q14 で扱いたい係数を 16-bit signed に詰めようとすると、

k_a = round(s_xx / det) ≈ 1,535,541   → 21 bit。16-bit に入らない。

詰めません。

ここで普通は「じゃ 32-bit でいいか」と進みがちですが、もう一手考えてみる価値があります。対角正則化（Tikhonov / shrinkage）: $\Sigma$ を $\Sigma + \varepsilon I$ に置き換える、よく知られた数値安定化です。逆行列が過敏になっている方向にわずかな対角成分を足すことで、過剰なノイズ増幅を抑えにいく操作で、別に変なハックではありません。

$\varepsilon$ を少しずつ大きくしながら係数が 16-bit signed に収まる最小値を探すと、$\varepsilon = 2\times10^{-5}$ で

k_a = 26188 (0x664C)
k_b = 23062 (0x5A16)
k_c = 26081 (0x65E1)

きれいに 16-bit signed に収まりました。代償は $\det\Sigma$ が約 66 倍に持ち上がり、$d^2$ のスケールが少し圧縮されることです。これが異常検知の感度をどれだけ削ったかは、後で実データで確かめます。

ここで決めた数値モデル全体:

mu_x_q = 2276    (= 0.139 g × 2^14)
mu_y_q = 2320
k_a    = 26188   (s_xx / det × 1)
k_b    = 23062
k_c    = 26081
T_q    = 11898   (健常 99.9 %ile を 32-bit acc 空間に換算したしきい値)

これだけ覚えてもらえれば、あとは整数演算だけで $d^2$ が転がせます。

PIC16F1503 に詰める — ハードウェア乗算器なし 8-bit の現実

PIC16F1503 のスペック:

• Enhanced mid-range 8-bit core, 14-bit 命令長、49 命令
• フラッシュ 2K words、SRAM 128 バイト（うち 80 バイトが Bank 0 GPR）
• ハードウェア乗算器なし（PIC18F 系には MULWF があるが、PIC16F1xxx は持っていない）
• 内部 HFINTOSC 8 MHz × 4 PLL = 32 MHz、命令サイクル = 4 クロック = 8 MIPS

つまり 8×8 → 16 bit の乗算もソフトウェアで書きます。これが面倒くさそうに見えて、案外そうでもない、というのが今回の発見でした。

中核は素朴な右シフト累積です（pseudo-asm）:

mul8x8u:        ; W = A,  mul_b = B,  pp_hi:pp_lo = A*B
        movwf   tmp_a
        clrf    pp_lo
        clrf    pp_hi
        movlw   .8
        movwf   cnt
m8_loop:
        rrf     tmp_a, f        ; A >>= 1, C = old bit 0
        btfsc   STATUS, C
        goto    m8_add
        bcf     STATUS, C
        goto    m8_shift
m8_add:
        movf    mul_b, w
        addwf   pp_hi, f        ; partial += B
m8_shift:
        rrf     pp_hi, f        ; result >>= 1, carry in
        rrf     pp_lo, f
        decfsz  cnt, f
        goto    m8_loop
        return

ループ 8 回で 16-bit 積が出ます。これを 4 回呼び合わせれば 16×16 → 32-bit 符号無し、符号管理を被せれば符号付きまで届きます。

$d^2$ 評価本体は、上で決めた係数を流して

acc = 0
prod = dx*dx (signed 32-bit)
prod = k_c * (prod >> 16)        ← 32-bit 中の上位 16-bit を取って再乗算
acc += prod
prod = dy*dy
prod = k_a * (prod >> 16)
acc += prod
prod = 2 * dx * dy
prod = k_b * (prod >> 16)
acc -= prod
anomaly := (acc >= T_q)

の 3 項を順に積みます。中間で >> 16 しているのが今回の妥協点で、32-bit 累算器に収めるための割り切りです（48-bit にすれば全精度で持てる、これは後段の話）。

gpasm でアセンブルしてサイズを測ると:

Program Memory Words Used:   230 / 2048   (11.2%)
RAM bytes used:                ~27 / 128  (21%)

これだけ。残りは ADC ISR と通信に回せます。

gpsim で 1 サイクルずつ追って、最悪値で測る

gpsim を使えば、組んだ HEX を 1 サイクル単位で踏みながら、レジスタとサイクルカウンタを追えます（実機 PIC を持っていなくても PC 内で完結します）。

worst-case として「データセット末期の Bearing 3 が壊れかけている瞬間」の入力（rms_x = 0.5936 g, rms_y = 0.4960 g）を与え、mah_evaluate ルーチンの入口と出口で cycles を読みました。

cycles at mah_evaluate entry  =   12
cycles at mah_evaluate return = 2805
→ 1 評価 = 2793 サイクル

32 MHz / 4 = 8 MIPS で割ると 約 349 μs。1 秒のセンサ窓 1 個に対して 0.035 % CPU しか使いません。残りはほぼ「何もしていない」時間で、超低消費電力モードで眠らせられます。百円台のチップでこれが回るのは、正直ちょっとびっくりしました。

結果のアキュムレータ値と異常フラグも、上で書いた Python 等価実装とビットレベルで一致することを確認しています（acc = 0x004DF697 = 5,109,399、result = 1 (anomaly)）。これは「コードを書いたつもりだったが実は別物だった」を防ぐ、最低限の儀式です。

35 日 × 2,156 窓を流して、判定の地に足を付ける

1 点だけ合ったから OK、ではちょっと心許ない。学習に使った最初の 215 窓を除いた残り全部、2,156 窓に上の整数演算を走らせて、正則化したフロート参照値（PIC が「理想的に」出したいはずの値）と異常判定を突き合わせます。

全 2156 窓:
  acc range          = [-16,917, +5,083,318]
  異常判定一致率      = 77.46 %
  miss (B あり, asm なし) = 486
  false positive      = 0
  最初に異常と判定した窓:
      参照 (B 系列)   = window 214
      asm (PIC 同等)  = window 246  (差 32 窓 = 約 5.3 時間)

ここで起きていることを言葉にすると、

• 空振り（false positive）はゼロ。これは産業現場では地味だが重要な性質で、誤報で人を動かさない。
• 見逃し（miss）は 486 窓あるが、これは「分布の中心からまだそれほど離れていない、しきい値ぎりぎり」のところで失っている。>> 16 で下位 16 ビットを刈った代償です。
• 異常の最初の検出は参照より 5.3 時間遅れる。run 全体が約 800 時間あって、Bearing 3 はその末尾 30 時間あたりで急速に死んでいくので、「最初の異常を 5 時間遅く検出する」は、catastrophic failure に至るまでの保守可能ウィンドウ内ではしばしば許容範囲、というのが私の読みです（が、ここは現場の SLA・許容停止時間で要相談）。

実際の検出の様子はこんな形になりました。上が x/y 軸の RMS、下が PIC 等価の $d^2$（log スケール）で、しきい値とそこを最初に超えた点を縦線で重ねています。

上のパネルを見るかぎり、RMS だけで「異常が始まった」と判断するのは、最初の検出時点（10/31）では正直けっこう難しい。0.13 → 0.17 g の変化は、現場の振動計の表示誤差レベルだと思います。一方、$d^2$ は分布の形（健常時の x/y 共動）を学習しているので、RMS には載らないわずかな崩れを拾います。これがマハラノビス距離を 2 軸に対して効かせている本体です。

下のパネルで $d^2$ が早期に何度か跳ね上がっているのは、上で書いた >> 16 の量子化ノイズが閾値付近で揺れているためです。実機で運用するなら、ここに「N 回連続で閾値超え」のヒステリシスを足すか、ローパスをかけて静かにします（asm 本体に対する追加コストは数十サイクルで済むはず）。

ちなみに asm の出力をそのまま物理スケールに戻すと、参照値とのずれは 平均 4.1、最大 10.7（参照値の値域 0.2〜1248）。大きい異常では十分小さい、小さい異常では効きが出る、という典型的な固定小数の挙動でした。

32-bit と 48-bit のあいだの一線 — ここが本当の設計判断

>> 16 を入れていなければ、上のずれはほぼゼロになります。$d^2$ の係数 $k$ と 32-bit の中間積を素直に掛け合わせると 48-bit の結果になるので、それを 6 バイト累算器にそのまま積めば、参照との一致率は

  異常判定一致率   = 99.72 %
  miss             = 5
  false positive   = 1
  最初の異常検出   = 参照と同じ window 214 (差 0)

ここまで上がります。代償は、累算器が 2 バイト増え、k × p の乗算が 16×16 から 16×32 に変わるので、1 評価あたりのサイクル数が 約 4,300（550 μs）に増えること。RAM は +2 byte、コードはおそらく +50 word 前後でしょう。

この「32-bit acc で 77% 一致 / FP=0 と、48-bit acc で 99.7% 一致 / FP=1」の選択は、用途で割れます。

• 電池駆動・ロギング寄り — 32-bit で十分。FP=0 なので運用が静か。
• 早期警報が命の用途（軸受の自動切替を握っているなど） — 48-bit に拡張する。RAM と命令空間の 1.5% を払って 5 時間早い検出を買う。

どちらも 百円台のチップに乗ります。これが今日いちばん書きたかったことです。

振り返って、何が見えたか

20 年くらい現場の電気・組込を見てきて、感じることがあります。「足りない計算量」より、「足りない統計知識」のほうが、設計の自由度を狭めている気がする。

今回の場合、計算量はそもそも問題ではありませんでした（2,800 サイクル）。詰まったのは $\rho=0.998$ という、データの素性のほう。これは現場でセンサを 2 個並べると当たり前に出てくる構造で、「素朴に Σ⁻¹ を計算したら、係数が爆発した」という現象は、 たぶん多くの組込屋がどこかで一度踏みます。

そのときに、

• 「8-bit には乗らないやめよう」と引き返すのは早い。
• 「32-bit に上げて押し切る」もたぶん早い。
• 正則化を 1 行入れて、削った感度を実データで定量する、までセットでやって初めて、削った/削ってないが議論になる。

異常検知に限らず、信号処理を MCU に詰めるときの設計判断は、たいていここに似たかたちで還ってくるように思います。

百円かそこらの 8-bit マイコンでも、35 日間のベアリング run-to-failure に対して、最終的な故障に到達するよりだいぶ前のタイミングで、健常分布から外れたことを示すフラグを上げられそうだ — 単発の早期警報を保証するというよりは、ここに運用ロジック（持続性、ヒステリシス、複数センサ間の票決）を一段被せれば、現場の判断材料には足る — という見通しは、私のような中年エンジニアにとっては、ちょっとした希望です。古い回転機、古いポンプ、古い空調機。あれらの上に、こういう小さな目を一個ずつ足していく仕事は、たぶんこれから 10 年くらい残っているのではないでしょうか。

コード: 今回の PIC asm（mahal.asm）、numpy リファレンス、Q-format 量子化と ε 探索、全 2,156 窓に対する sweep スクリプト、図の生成までを一通り github.com/logicia32/pic-bearing-anomaly に置いてあります。MIT。

データ: NASA Prognostics Center of Excellence の IMS Bearing Data。公式の入口は PCoE Data Set Repository で、配布されている zip の中身は 7z（さらに RAR）の入れ子になっています。README に解凍手順を書いておきました。