論文
隠れマルコフモデルとエリアセンシング手法を用いた未知の異常データ検出
2020 年 106 巻 2 号 p. 91-99
詳細
2020 年 106 巻 2 号 p. 91-99
In this paper, we propose an anomaly detection method for equipment abnormalities using data measured by AreaSensing techniques. The data are gathered from the vibration of equipment installed in a wide area or large scales such as conveyor equipment and bridges. It is difficult to know which data of displacement, velocity, and acceleration is appropriate in advance because each frequency component is different. So, we apply the Hidden Markov Model, which estimates the latent state for each decomposition level by continuously frequency-resolving time series data using Wavelet method. The analysis results show that the normal and abnormally states are estimated. However, as the problem of this method, it is not possible to compare and apply the information criteria such as AIC to the Wavelet decomposed data as it is to appropriately decide which data and model parameters should be used. To overcome the defects, we propose a new evaluation function and developed a method to find a model that can stably estimate the normal and abnormal state transition, even for data separated into different frequencies. Besides, when the current measurement data contains no abnormal state, there was a problem of extracting multiple latent states that are normal but different. We focus on the difference between the state transition probabilities of the normal and unknown model. As the experimental result, the effectiveness of the proposed method has been confirmed. By using the method, it is possible to continuously diagnose abnormalities using vibration measurement data measured by AreaSensing techniques.
コンベア設備や橋梁など，広範囲あるいは大規模に敷設されている設備の点検保守作業には，多大な人手がかかっている。また，高温や振動を伴う設備では，異常を検出するためのセンサーの設置が困難な場合がある。このような設備に対して，カメラでの高速撮影画像1)やモアレ検出機能2)を利用して遠隔で振動計測を行うエリアセンシング技術が開発されてきた。高速撮影が可能なアクティブ・ビジョン・システムは，高速カメラによって撮影された輝度値の時間変動が，振動現象に起因する情報を含むものと想定し，カメラによって振動計測を行うシステムである。構造物の複数の点で動的に変化する変位と振動を測定できるマルチスレッド・アクティブ・ビジョン・センシング技術によって，ミリ秒レベルでの超高速視点切り替え技術を用いて，同時に複数の振動測定が可能であり，広範囲にわたって分布する小さな振動を同時に測定することができる。サンプリングモアレカメラは，モアレ現象を利用して，測定対象物と離れた場所から高速変位測定ができる機能を持ち，構造物の振動波形データを高精度・リアルタイムに監視することができる。本稿では，これらのエリアセンシング技術によって測定された振動データから，コンベア設備の未知の異常状態を発見する手法の提案とその実験結果について報告する。
本稿では，2章で関連研究について述べた後，3章で提案するモデル選択手法について説明する。そして，4章では，コンベアデータからの異常検出について説明し，5章で提案手法の課題と改善結果を示し，最後に6章で結論をまとめる。
生産設備に限らず経済事象や社会現象など，様々な状態を観測した時系列データから異常検知を行う計算手法に関する取り組みは多く行われてきた。ガウス混合分布などを用いて外れ値を検出する方法や，異常変数の選択手法に対応したマハラノビス=タグチ法，時系列上の急激な変化を観測する変化点検出などが一般的に適用されてきた3–8)。また，観測空間では認識できない潜在的な異常検知の手法として，潜在変数を利用した事例も報告されている9,10)。
鉄鋼業などの設備異常に関しても，状態空間モデルや主成分分析を適用して溶鉱炉のモニタリングを行う研究11–13)や，ウェーブレットによる時間・周波数解析を適用した研究14–17)，サポートベクターマシンによる故障診断手法の研究18,19)，ニューラルネットワークを使った連続鋳造工程の異常発見手法の研究20)などが報告されている。また，橋梁振動の異常検出にベイズ推論を用いた研究21)，道路や橋梁などの社会基盤施設の異常発見に統計的変化点検出や特徴パラメータ診断の手法を用いた研究22–25)，k-近傍データ距離を利用した局所部分空間法を用いた異常検知技術26)，定常状態監視データの発生確率を利用した異常検知システム27)，エネルギー監視データを利用した設備異常検知28)など，様々な手法による異常検知技術が提案されてきた。
また，隠れマルコフモデルを利用した異常検知技術も発展してきており，ベアリング異常に適用した事例29)や，周波数帯に分けて隠れマルコフモデルを木構造にして異常判定する研究30)，銀行株変化率の時系列データを対象として隠れマルコフモデルを適用し，その変化の背後にある状態として，2値の離散値として銀行株の強気・弱気相場を推定する手法などが提案されている31)。これらの手法は，正常・異常状態の背後にある真の状態を推定することができる点で優れているが，周波数帯域や分布モデルなどのハイパーパラメータの最適化に課題があり，適切なモデル選択が難しい。そこで，本研究では，これらの課題を解決することを目的として，新たなモデル評価基準の提案と実データでの検証を行う。
2・2 時系列異常データ推定モデル時系列データの特徴パターンをモデル化する手法の一つに，隠れマルコフモデル(Hidden Markov Model)がある32–34)。このモデルの基盤となっているマルコフ連鎖モデルは，確率過程として離散的な状態を扱い，未来の状態が現在の状態の値だけで決定される独立なものと考えるモデルであるが，これは，現在の状態が直接観測可能であることを前提としている。通常の混合モデルでは，パラメータθが与えられた後の各データX={x1,...,xN}の分布に対して，条件付き独立性が仮定されていた。しかし，異常診断などの場合では，時系列データを発する背後に，正常か異常かの状態があり，多くの場合それらを直接観測することができない。これに対して，隠れマルコフモデルは，状態は直接観測される必要はなく，確率事象だけが観測されることをモデル化している。この隠れた状態を表す確率過程を潜在変数と呼び，それぞれは独立ではなく，マルコフ過程を通じて連携する混合分布モデルとなっている。各潜在変数S={s1,...,sN}の発生過程に依存性を持たせることで，隠れマルコフモデルを表現できる(式(1))。
| (1) |
ここで，潜在変数Sはマルコフ過程の状態数Kの状態系列であり，状態遷移確率行列Aを持つ。πは，先頭の状態S1を決める初期確率パラメータである。
時系列データを対象とした隠れマルコフモデル(HMM)の一般的な構成は，Fig.1のようになっている。この図は，left-to-right HMMと呼ばれるモデルで，初期状態確率p(S1)=1，p(Sj)=0(j≠1)となっており，全ての列は状態j=1から始まるように制限されている。
Hidden Markov Model (left-to-right).
一方で，1つ隣の状態にのみ依存する1次マルコフ連鎖の場合は，状態遷移確率AはK×Kの行列で表現され，状態iから状態jへの遷移確率をaijで表すことができる。本モデルでは，異常状態から正常状態へ戻るような逆方向の遷移も考えられる。
ところで，コンベア設備の振動データなどの場合，モーターやベアリング，支柱などの振動の発生源や，固有振動数によって，観測される周波数が異なってくることが予想される。しかし，事前にどの周波数領域が怪しいのかを特定することは困難であり，また，どの時間帯で異常状態に遷移するのかも予測ができない。通常のFFT解析では，時系列のある部分を対象に周波数解析をする必要があり，ウィンドウを動かして時系列の全領域を探索する場合でも，ウィンドウ幅の粒度での分解となる。そこで，本報告では，長時間継続する時系列データを連続的に周波数解析ができるウェーブレット(Wavelet)解析手法を採用する。この手法は，マザー・ウェーブレットと呼ばれる小さな波形ψ(x)を平行移動と伸縮をして原波形に掛け合わせることで時間周波数解析を行う。これを計算が容易になるように離散問題に適用した離散ウェーブレット変換の一つに，多重解像度解析(Multiresolution analysis)がある。これは，基底関数となるマザー・ウェーブレットを用いて，2倍毎の解像度波形を連続的に得る手法である。Haarのスケーリング関数φHを式(2)に定義する。
| (2) |
レベルjの数列{ck(j)}が知られていて，この数列からレベルj-1の数列を次のように求めれば，
| (3) |
| (4) |
これを利用して，関数fjは式(5)のように分解できる。
| (5) |
ここでfj(x)とgj(x)は式(6)，(7)で与えられる。
| (6) |
| (7) |
それぞれのkについて計算するこれらの式は，分解アルゴリズムと呼ばれる。
Fig.2は，Wavelet変換手法の一つである多重解像度解析(離散ウェーブレット変換とも呼ばれる)を行なったコンベア振動データの解析結果である。最下層が加速度センサーを使って測定されたコンベア振動データで，上位層に行くに従って，スケーリング関数を使って信号列を半分の解像度に分解した時系列波形となっている。生波形とは異なる波形の特徴が，分解された波形では現れている。そこで，Wavelet変換を行ったそれぞれの解像度波形に対して隠れマルコフモデルを適用することで，未知の異常データの検出を行う。これは，観測された周波数別時系列データの背後にある，時間経過と共にマルコフ過程で遷移する未知の定常状態の推定を行うことで，設備振動データからの異常状態を発見しようとするものである。
Multiresolution analysis of conveyor oscillation data.
しかしながら，本報告で用いる手法にはモデル選択において以下の課題がある。一般的なモデル選択手法の場合，AICやBICといった情報量基準を用いて，推定精度とモデルの複雑さを考慮した適切なモデルを選択することができた。そのためには，分析対象とするデータが各モデルで同一であることが前提となっていた。ところが，本提案手法で用いる多重解像度解析は，信号列としてのデータそのものを解像度ごとに分解してしまうため，モデルに投入するデータ自体が変わってしまい，情報量基準をモデル選択の比較に用いることができない35)。この課題を解決するために，次節で新たな手法を提案する。
本節では，隠れマルコフモデル分析における異なるデータを対象としたモデル選択の手法について述べる。前説で指摘したように，情報量基準のみを使ったモデル選択は，離散ウェーブレット変換を行なったデータに対して単純には使用できない。また，多重解像度解析の解像度レベルや隠れマルコフモデルに用いる分布の形状と分布の数など，多くのハイパーパラメータが存在し，その全組み合わせを試行錯誤する必要がある。この労力を削減し，かつ異常発見に適したモデルパラメータを選択することを目的に，変数選択・ハイパーパラメータ選択のための評価関数について検討する。
3・1 評価関数異常診断の観点から，あまり頻繁な状態遷移は好ましくないと考えられる。そこで，ある状態iから別の状態jへ遷移する確率が，状態iから同じ状態iに遷移する確率よりも十分に小さいHMMモデルが望ましいと考え，以下の式(8)を評価関数として用いることを提案する。
| (8) |
ここで，Sは状態の集合，p(i, j)は状態iからjへの遷移確率であり，evalは値が小さいほど優れている。
3・2 数値実験前節で提案した評価関数の妥当性を検証するため，橋梁モデル振動計測データとコンベア振動計測データを対象に，HMMモデルのハイパーパラメータおよび，HMMモデルの学習に用いる変数を網羅的に変更して，評価関数の妥当性を検証する。
分析では，広島大学石井研究室で開発されたアクティブ・ビジョン・システムによる橋梁模型からの高速視覚振動データを使用した。マーカーが取り付けられた4 m長のトラス構造ブリッジモデルを振動させ，33.3 fpsで動作する15台の仮想カメラが動作するアクティブ・ビジョン・システムが，ブリッジ構造の変形を観測することで，橋梁の模擬振動データを取得した。測定されたデータの中から，正常な状態の振動データと，橋梁の中央部分の構造の一部の継ぎ目を外し，構造物のクラックの発生を模擬した2種類のデータを使用した。振動発生はハンマー殴打によるインパルス応答データとして測定されているため，応答波形部分を切り出し，クラックあり・なしのそれぞれ5回の繰り返し振動データを作成し，隠れマルコフモルを用いた潜在状態推定を行った。Fig.3に分析に用いた波形を示す。前半1/3がクラックあり，後半2/3がクラックなしの5回繰り返しデータである。
Oscillatory waveform of a bridge.
前節で示した橋梁モデル振動計測データを対象に，確率分布としてdist={正規分布(NOR)，混合正規分布(MIX)}，状態数としてnStates={2, 3, 4}，確率分布が混合正規分布の場合の正規分布の数としてnMixt={2, 3, 4}，確率分布が正規分布の場合の初期化方法としてinit={RANDOM，KMEANS}の全組合せについて，evalおよびBICを計算する。
Table 1に，評価関数evalの観点から優れたもの上位7個のモデルのevalおよびBIC，Viterbiアルゴリズムを適用した際のグラフを示す。1，2位のモデルは，前半のクラックありの状態と後半のクラックなしの状態を正しく分離できている。この結果は，評価関数の妥当性を示唆していると考えられる。
Results of applying the proposed evaluation function to the data of the bridge model.
コンベア振動計測データを対象に，ウェーブレット変換および隠れマルコフモデルを用いた異常検出を試みる。データは，エリアセンシング技術の一つとして福井大学藤垣研究室で開発されたサンプリングモアレカメラ2)によるベルトコンベア振動計測データを用いた。これは，ベルトコンベアの支柱に貼られた格子状のターゲットを，5 mの距離からサンプリングモアレカメラで計測し，変位を観測したもので，サンプリング周期は80 Hzとなっている。どちらもコンベアは動作中，データの前半は在荷，後半は空荷となっている。前半部分と後半部分のそれぞれ30秒間である。Fig.4に測定した振動データの変位，速度，加速度の時系列グラフを示す。また，加速度データに対してWavelet分解を行なった結果をFig.5に示す。
Vibration data (load/empty), displacement (above), velocity (middle), acceleration (lower).
Multiresolution analysis of conveyor oscillation data (load/empty).
上記のデータに対して，4・1節で定義した評価関数evalを用いて，変位，速度，加速度の各計測データ，それぞれのWavelet分解レベル，隠れマルコフモデルのハイパーパラメータ(推定分布種類，分布数，分布の初期値)の組み合わせの中から，最もeval値が小さいものを選択することを行なった。前述のように，eval関数を用いると，ある状態iから別の状態jへ遷移する確率が，状態iから同じ状態iに遷移する確率よりも小さくなるモデルを選択することができる。このことにより，正常状態と異常状態を時系列のある時点で明確に分離することが可能となる。実験では，測定したコンベア振動計測データに対して，Wavelet関数で5レベルまで周波数分解した。そして，隠れマルコフモデルの推定分布として，正規分布と混合正規分布を用いてeval関数で，全ての組み合わせで評価を行った。その結果，最良のeval値になったのは，変位データの場合はWaveletレベルが0(分解なし)，速度および加速度データの場合はWaveletレベルが3，そして，いずれの場合も，隠れマルコフモデルは正規分布および潜在状態数が2となった。加速度データにおける隠れマルコフモデルの状態遷移図をFig.6，それぞれの最小eval値とその時の分析モデルのハイパーパラメータをTable 2に示す。
HMM State transition model.
| Wavelet level | HMM distribution | Hidden state | Eval | BIC | |
|---|---|---|---|---|---|
| displacement | 0 | Normal | 2 | 0.0004 | –17308 |
| velocity | 3 | Normal | 2 | 0.026 | –3598 |
| acceleration | 3 | Normal | 2 | 0.012 | –3852 |
Table 3に状態遷移確率を，Fig.7には，変位，速度，加速度データでの潜在状態の推定結果を示す。速度データの場合は，前半の空荷の状態でも，状態推定が安定していない。変位と加速度データの場合は，前半と後半でほぼ正確に推定できていることがわかる。設備診断の場合，振動周波数によって変位振幅，速度振幅，加速度振幅の応答が違うため，軸受の傷による振動や，アンバランスや緩み，回転軸のふれなど，どの振動が増加するのかを検出したい異常に応じて選択する必要がある。今回のコンベア振動計測の場合は，変位と加速度に変換した振動データが，積荷の有無を正確に推定することができた。
| state | S1 | S2 | |
|---|---|---|---|
| displacement | S1 | 0.946 | 0.054 |
| S2 | 0.578 | 0.942 | |
| velocity | S1 | 0.999 | 0.001 |
| S2 | 0.003 | 0.997 | |
| acceleration | S1 | 0.974 | 0.026 |
| S2 | 0.063 | 0.937 |
The result of HMM analysis. displacement (above), velocity (middle), acceleration (lower). (Online version in color.)
情報量基準を使った手法では，分析対象のデータが周波数分解などで変形された場合，モデルの適切さを比較することが不可能であった。本手法では，eval関数の導入によって，Waveletによって分解されたコンベア振動データの潜在的な2状態を正しく分解することが可能となった。しかし，これは正常な状態(積荷)と異常な状態(空荷)の2状態が存在するという情報を，事前に得た上でモデル化した分析である。実際には，正常か異常かがわからない状態で日々の監視を行わなければならない。測定したデータの全てが正常状態であるかもしれないし，全てが異常状態かもしれない。あるいは，正常と異常の両方が含まれているデータなのかもしれない。そこで，これらの課題に対処するために次のような分析手法を検討する。
5・1 未知の振動データの異常推定法隠れマルコフモデルの状態遷移確率および各状態からの系列データ出力確率をデータとして利用して，既知の正常なデータから作成した時系列モデルと，正常か異常かがわからないデータから作成した時系列モデルの比較を行い，正常状態と異常状態を判断する。Fig.8にその考え方を示す。
Decision model of normal and abnormal state. Mk: Model of known normal data, Mu: Model of unknown data. (Online version in color.)
正常であることがわかっている時系列データだけから作成したモデルをMki，i=1...mとする。このモデルは，データとして変位，速度，加速度と，Waveletの各レベル，隠れマルコフモデルの分布などのハイパーパラメータの全組み合わせ(全m個)に対して分析を行い，それぞれの状態遷移確率を保存しておく。次に，正常か異常かが未知の時系列データからevalが最小となるモデルMuを作成する。そして，そのモデルと同じハイパーパラメータを持つMkiと比較する。もし，未知の時系列が正常ならば2つのモデルの状態遷移確率は近似し，異常ならば2つのモデルは異なることが予想される。そこで，以下の式(9)のように，状態遷移確率の平均自乗誤差MSEを定義する。
| (9) |
ここで，skij，suijは既知モデルと未知モデルの状態遷移確率，kは潜在状態数である。このMSEを使って，正常な状態から測定した二つの時系列データ間のMSEと，正常状態と未知状態から測定した二つの時系列間のMSEを比較する。
5・2 未知の振動データの異常推定実験結果Fig.9にこの手法を用いて分析した結果を示す。図の左から，変位データにおける正常モデルと未知モデルのMSE(dis:abnormal)，変位データにおける正常モデル間のMSE(dis:normal)，速度データにおける正常モデルと未知モデルのMSE(velo:abnormal)，速度データにおける正常モデル間のMSE(velo:normal)，加速度データにおける正常モデルと未知モデル間のMSE(acc:abnormal)，加速度データにおける正常モデル間のMSE(acc:normal)である。どれもt検定の結果はp<0.01であった。これらは，測定データに未知の異常データが含まれた場合は，有意な差が出ていることを示しており，提案モデルのよって解析した状態遷移確率から，未知の正常ではない状態を発見することができている。
Comparison of state transition probability MSE. dis.abnormal: Displacement abnormal data, dis.normal: Displacement normal data, velo.abnormal: Velocity abnormal data, velo.normal: Velocity normal data, acc.abnormal: Acceleration abnormal data, acc.normal: Acceleration normal data. (Online version in color.)
本研究では，コンベア設備や橋梁など，広範囲あるいは大規模に敷設されている設備の振動を，エリアセンシング技術のアクティブ・ビジョン・システムやサンプリングモアレカメラで計測したデータを用いて，設備異常発見を行う新しい手法を提案した。振動データは，回転軸受けの傷やアンバランス，ビビリ現象や回転のふれなどで，変位，速度，加速度のどのデータが適切か事前に知ることができない。これは，それぞれの周波数成分が異なることによる。そこで，本提案では，Wavelet手法を用いて，時系列データを連続して周波数分解し，それぞれの分解レベルに対して，潜在状態を推定できる隠れマルコフモデルを適用することを試みた。分析の結果，正常状態と異常状態を適切に推定できることを示した。しかし，この手法の課題として，Wavelet分解をしたデータに対しては，AICやBICといった情報量基準をそのまま適用して比較することができず，どのデータとモデルパラメータを用いるべきかを適切に決めることができない。そこで，新たにeval関数を提案し，異なる周波数に分解されたデータに対しても，正常・異常の状態遷移を安定して推定することができるモデルを発見する手法を開発した。また，現在の測定データに異常状態が含まれていない場合，どれも正常だが性質の異なる複数の潜在状態を抽出してしまう課題に対して，正常モデルと未知モデルの状態遷移確率の違いに着目した手法を提案し，その有効性を確認した。
本提案手法を用いることで，エリアセンシング技術で計測した振動計測データを用いて，連続的に異常診断を行うことが可能となる。今後の課題として，多種多様な設備でのコンベア振動データに対して，その有効性を確認することと，モデルのハイパーパラメータがより複雑になるモデルに対して，連続変数と質的変数を含むブラックボックス関数最適化のための進化計算手法を適用して，より高速にモデル推定ができる手法を確立することを検討している。
本研究は，一般社団法人日本鉄鋼協会の支援により実施された。ここに謝意を表する。
