JP7017363B2

JP7017363B2 - 異常検知装置および異常検知方法

Info

Publication number: JP7017363B2
Application number: JP2017196333A
Authority: JP
Inventors: 久恵渋谷
Original assignee: Hitachi Power Solutions Co Ltd
Current assignee: Hitachi Power Solutions Co Ltd
Priority date: 2017-10-06
Filing date: 2017-10-06
Publication date: 2022-02-08
Anticipated expiration: 2037-10-06
Also published as: JP2019070930A

Description

本発明は、プラントや設備などが出力する複数の時系列センサ信号をもとに異常を早期に検知する異常検知装置および異常検知方法に関する。

電力会社では、ガスタービンの廃熱などを利用して地域暖房用温水を供給したり、工場向けに高圧蒸気や低圧蒸気を供給したりしている。石油化学会社では、ガスタービンなどを電源設備として運転している。このようにガスタービンなどを用いた各種プラントや設備において、設備の不具合またはその兆候を検知する異常検知は、社会へのダメージを最小限に抑えるためにも極めて重要である。

ガスタービンや蒸気タービンのみならず、水力発電所での水車、原子力発電所の原子炉、風力発電所の風車、航空機や重機のエンジン、鉄道車両や軌道、エスカレータ、エレベータ、機器・部品レベルでも、搭載電池の劣化・寿命など、上記のような予防保全を必要とする設備は枚挙に暇がない。

このため、対象設備やプラントでは様々な物理情報を取得する複数のセンサを取り付け、センサ毎の監視基準に従って、対象設備やプラントが正常か異常かを判定される。特許文献１には、過去の正常データの学習に基づいて算出される異常測度をしきい値と比較することによって異常の有無を検知する異常検知方法が開示されている。さらに、検知した異常について、対策・調査などの次の行動の決定を支援する情報を提示することを目的として、センサ信号の二次元の分布に基づいて異常に関連するセンサを特定する方法が開示されている。ここに異常測度とは、複数のセンサによる測定値を１つのベクトル値として表現し、正常状態のベクトル値からの偏移量のことである。

特開２０１４－１４２６９７号公報

特許文献１に記載の二次元の分布に基づいて異常に関連するセンサを特定する方法は、検知された異常に関連するセンサを正しく診断することが可能であるため、作業者はまずその関連を調査し、異常の原因を推定し対策することができる。しかし、予兆の早期検知を狙い検知の感度を上げると、誤報が発生するとともに異常区間が多数出力されることとなるため、どこから調査すべきか判断するのが難しくなる。また、人為的な条件変更や一時操作による変化など対策不要の状態についても異常として検知してしまうため、確認の手間が余計にかかってしまう。

本発明の目的は、上記課題を解決するため、複数の時系列センサ信号に基づく異常検知において、異常区間が多数検知されても異常に関連するセンサの特定を容易にする異常検知装置および異常検知方法を提供することにある。

本発明の異常検知装置は、設備に装着された複数のセンサから出力される複数の時系列のセンサ信号を入力するセンサ信号入力部と、前記センサ信号から時刻毎に特徴ベクトルを抽出する特徴ベクトル抽出部と、指定された学習期間の前記特徴ベクトルを学習データとして各時刻の異常測度を算出する異常測度算出部と、前記異常測度をしきい値と比較することにより各時刻のセンサ信号が正常か異常かを判定する異常検出部と、異常が連続して検知される異常区間を抽出し、前記異常区間毎に異常との関連性の強さを示すパラメータを算出し、前記パラメータに基づき異常検知の原因となった関連センサを特定する関連センサ特定部と、前記異常区間を前記パラメータの類似度に基づいて複数のクラスタに分類し、含まれる異常区間数が多いクラスタほど重要度が高いと推定して、診断支援情報に順位を付けて提示する診断情報提示部と、を備える構成とする。

本発明によれば、重要度が高いと推定される関連センサから優先的に提示するので、ユーザの確認作業の効率を向上させることができる。

本発明に係る異常検知装置の一構成例を示す図である。複数のセンサ信号をリスト化して表形式に表した例を示す図である。異常検知装置の行う全体の処理フローを示す図である。学習時の異常測度算出処理のフローを示す図である。局所部分空間法による異常測度算出処理を説明する図である。学習時の二次元頻度分布算出処理のフローを示す図である。二次元頻度分布画像の例を示す図である。二次元頻度分布画像の例を示す図である。二次元頻度分布画像の例を示す図である。異常検知処理のフローを示す図である。孤立度算出処理のフローを示す図である。孤立度の算出法を具体的に説明する図である。関連センサ特定処理のフローを示す図である。重要度推定処理のフローを示す図である。オフライン解析条件を設定するＧＵＩの例を示す図である。オンライン解析結果の表示対象を指定するＧＵＩの例を示す図である。解析結果全体表示画面の例を示す図である。解析結果拡大表示画面の例を示す図である。クラスタ表示画面の例を示す図である。重要度診断画面の例を示す図である。対策不要と誤報を推定する処理のフローを示す図である。対策不要と推定すべき例を示す図である。複数の設備の異常監視結果を一覧するＧＵＩの例を示す図である。

図１は、本発明に係る異常検知装置の一構成例を示す図である。異常検知装置１００は、検知対象である設備１０１に装着されたセンサから出力されるセンサ信号１０２を、所定時間ごとに（周期的に）取得する。取得したセンサ信号１０２は、一旦センサ信号蓄積部１０３にて蓄積された後、あるいは直接にセンサ信号入力部１０４に入力され、特徴ベクトル抽出部１０５へ送られる。特徴ベクトル抽出部１０５は、センサ信号１０２をもとに特徴ベクトルを抽出し異常測度算出部１０６へ送る。異常測度算出部１０６は、予め指定された学習期間の特徴ベクトルを用いて、所定時間毎（以下、各時刻と表現する場合もある）の特徴ベクトル毎に異常測度を算出する。

しきい値算出部１０７は、異常測度算出部１０６による学習データの異常測度に基づいてしきい値を算出する。二次元分布算出部１０８は学習期間のセンサ信号に基づき、センサ２個の全組合せの二次元頻度分布を算出する。特徴ベクトル抽出部１０５で抽出された学習期間の特徴ベクトル、しきい値算出部１０７で算出されたしきい値、二次元分布算出部１０８で算出された二次元頻度分布ほか、異常検知時に必要となるデータは学習結果として学習結果蓄積部１０９に保存される。異常検出部１１０は、異常測度算出部１０６から送られる各特徴ベクトルの異常測度と、しきい値算出部１０７で算出したしきい値とを比較することで、設備１０１の異常を検出する。関連センサ特定部１１１は、二次元分布算出部１０８で算出された二次元頻度分布を用いて、異常区間毎に異常関連センサを特定する。診断情報提示部１１２は、異常区間をクラスタリングして、高頻度なクラスタほど重要度が高いと推定して重要度順位を決定し、関連センサおよび重要度の情報を含む検知結果１１３を診断支援情報として出力する。

ここで、以下で用いる用語の簡単な説明を行う。特徴ベクトルとは、複数のセンサによる測定値を１つのベクトル値として表現したものである。異常測度とは、注目する特徴ベクトルの、指定された期間の特徴ベクトルからの偏移量のことである。異常区間とは異常が連続して検知される区間のことである。孤立度とは、センサ信号値が学習データからどれだけ乖離しているかを示すパラメータであり、異常との関連性の強さを示すものである。異常関連センサ（関連センサ）とは、異常測度がしきい値を超えて異常として検知された原因となったセンサのことである。クラスタとは、異常区間を孤立度ベクトルの類似度に基づいて分類した集合体で、そこに含まれる異常区間数が多いほど重要度が高いと推定する。

異常検知の対象とする設備１０１は、例えばガスタービンや蒸気タービンなどの設備やプラントである。設備１０１は、その状態を表すセンサ信号１０２を出力する。センサ信号１０２はセンサ信号蓄積部１０３に蓄積される。

図２は、複数のセンサ信号１０２をリスト化して表形式に表した例である。センサ信号１０２は、物理特性の異なる複数の物理情報が所定間隔毎に取得される多次元時系列信号である。図２に示す表の構成は、日時２０１の情報と、複数のセンサのセンサ値２０２を対応させて示している。センサは、数百から数千といった数になる場合もあり、それらの種類によって、例えば、シリンダ、オイル、冷却水などの温度、オイルや冷却水の圧力、軸の回転速度、室温、運転時間などをセンサ値として出力する。センサ値は、設備やプラントなどの出力や状態を表すのみならず、何かの状態をある値（たとえば目標値）に制御するための制御信号の場合もある。

異常検知装置１００の動作には、センサ信号蓄積部１０３に蓄積されたデータを用いて学習データの生成、保存を行う「学習」処理と、入力信号に基づき異常を検知する「異常検知」処理の二つのフェーズがある。基本的に前者はオフラインの処理、後者はオンラインの処理である。ただし、後者をオフラインの処理とすることも可能である。以下の説明では、それらを「学習時」、「異常検知時」という言葉で区別する。

図３は、異常検知装置１００の行う全体の処理フローを示す図である。ここでは処理の概要を記載している。（ａ）は学習時の異常測度算出処理で、学習期間のセンサ信号を入力し（Ｓ３０１）、特徴ベクトルの抽出（Ｓ３０２）と異常測度の算出（Ｓ３０３）を行う。（ｂ）は学習時の学習結果算出処理で、Ｓ３０３で求めた異常測度の値を利用してしきい値の算出（Ｓ３１１）を行い、また、学習期間のセンサ信号を入力し（Ｓ３１２）、センサ２個の全組合せの二次元頻度分布の算出（Ｓ３１３）を行う。（ｃ）は異常検知時の異常判定処理で、検知対象のセンサ信号を入力し（Ｓ３２１）、特徴ベクトルの抽出（Ｓ３２２）と異常測度の算出（Ｓ３２３）を行う。そして、算出した異常測度を、Ｓ３１１で求めたしきい値と比較することにより、設備の正常／異常を判定する（Ｓ３２４）。そして、異常区間を抽出し（Ｓ３２５）、異常区間毎に孤立度ベクトルを算出し（Ｓ３２６）、これをもとに異常関連センサを特定する（Ｓ３２７）。最後に、異常区間のクラスタリングに基づき重要度の推定を行う（Ｓ３２８）。
以下、（ａ）（ｂ）（ｃ）の順に説明するが、それぞれの詳細なフローは、図４Ａ、図５、図７～８、図１０～１１にて説明する。

図４Ａは、学習時の異常測度算出処理のフローを示す図である。最初に、センサ信号入力部１０４において、センサ信号蓄積部１０３に蓄積されたセンサ値のうち指定された期間（学習期間）のセンサ信号を入力する（Ｓ４０１）。学習期間として、設備が正常な状態であった期間を指定するものとする。次に、特徴ベクトル抽出部１０５において、入力されたセンサ信号を正準化する（Ｓ４０２）。センサ信号の正準化は、単位及びスケールの異なる複数のセンサ信号を同様に扱うために行う。具体的には、各センサ信号の、学習期間の平均と標準偏差を用いて、平均が０、分散が１となるように各センサ信号を変換する。異常検知時に同じ変換ができるように、各センサ信号の平均と標準偏差を記憶しておく。または、各センサ信号の、学習期間の最大値と最小値を用いて、最大が１、最小が０となるように各センサ信号を変換する。または、最大値と最小値の代わりに予め設定した上限値と下限値を用いてもよい。この場合は、異常検知時に同様の変換ができるように、各センサ信号の最大値と最小値または上限値と下限値を学習結果蓄積部１０９に記憶しておく。

次に、特徴ベクトル抽出部１０５において、各時刻の特徴ベクトルを抽出する（Ｓ４０３）。特徴ベクトルは、センサ信号を正準化したものをそのまま要素として並べたものである。あるいは、ある時刻に対して±１，±２，・・・のウィンドウを設け、ウィンドウ幅（３，５，・・・）×センサ数の特徴ベクトルとすることで、センサ信号の時間変化を表す特徴を抽出することもできる。また、離散ウェーブレット変換（ＤＷＴ：Discrete Wavelet Transform）を施して、周波数成分に分解してもよい。

次に、異常測度算出部１０６において、学習期間の異常測度を算出する。まず学習期間を複数の区間に分け（Ｓ４０４）、抽出した全特徴ベクトルについて、以下の処理を繰り返す（Ｓ４０５）。複数区間に対応して順次選んだ特徴ベクトルである注目ベクトルと、注目ベクトルと同じ区間を除く学習期間のデータを学習データとする（Ｓ４０６）。注目ベクトルと学習データを用いて異常測度を算出する（Ｓ４０７）。ステップＳ４０４における区間の分割は例えば１日毎とする。あるいは、化学プラントのようなバッチ処理の場合はバッチ毎、加工装置の場合は加工対象個体毎、ＭＲＩのような医療装置の場合は検査対象者毎としてもよい。ステップＳ４０７の異常測度算出処理には、局所部分空間法（ＬＳＣ：Local Sub-space Classifier）や投影距離法（ＰＤＭ：Projection Distance Method）を用いることができる。

図４Ｂは、局所部分空間法による異常測度算出処理を説明する図である。局所部分空間法は、注目ベクトルｑに対するｋ個の近傍ベクトルを選択し、選択したｋ個の近傍ベクトルが張るｋ－１次元のアフィン部分空間へ注目ベクトルｑを投影したときの投影距離を測る方法である。図４Ｂでは、ｋ＝３個の近傍ベクトルｘ１～ｘ３でアフィン部分空間を形成した場合である。そして、注目ベクトルｑに最も近いアフィン部分空間上の点Ｘｂが投影点（基準ベクトル）となり、注目ベクトルｑから基準ベクトルＸｂまでの距離が異常測度である。

具体的な算出法を説明する。評価データｑとそのｋ個の近傍ベクトルｘｉ（ｉ＝１，・・・，ｋ）から、ｑをｋ個並べた行列Ｑとｘｉを並べた行列Ｘを作成し、（１）式から両者の相関行列Ｃを求める。次に、（２）式から近傍ベクトルｘｉの重み付けを表す係数ベクトルｂを計算する。異常測度ｄは、ベクトル（ｑ－Ｘｂ）のノルムまたはその２乗により算出される。

なお、図４Ｂではｋ＝３の場合を説明したが、特徴ベクトルの次元数より十分小さければいくつでもよい。ｋ＝１の場合は、最近傍法と等価の処理になる。

投影距離法は、選択された特徴ベクトルに対し独自の原点をもつ部分空間すなわちアフィン部分空間（分散最大の空間）を作成する方法である。何らかの方法で注目ベクトルに対応する複数の特徴ベクトルを選択し、以下の方法でアフィン部分空間を算出する。

まず、選択された特徴ベクトルの平均μと共分散行列Σを求め、次にΣの固有値問題を解いて、値の大きい方から予め指定したｒ個の固有値に対応する固有ベクトルを並べた行列Ｕをアフィン部分空間の正規直交基底とする。ｒは特徴ベクトルの次元より小さくかつ選択データ数より小さい数とする。またはｒを固定した数とせず、固有値の大きい方から累積した寄与率が予め指定した割合を超えたときの値としてもよい。注目ベクトルから最も近いアフィン部分空間上の点が基準ベクトルとなる。また、注目ベクトルから基準ベクトルを引いたものが残差ベクトルとなり、残差ベクトルのノルムまたはノルムの２乗が異常測度となる。

ここで、複数の特徴ベクトルの選択方法としては、予め指定した数十から数百の数の特徴ベクトルを注目ベクトルから近い順に選択する方法がある。また、学習対象の特徴ベクトルを予めクラスタリングしておき、注目ベクトルに最も近いクラスタに含まれる特徴ベクトルを選択するようにしてもよい。また、注目ベクトルｑのｋ－近傍ベクトルの平均ベクトルまでの距離を異常測度とする局所平均距離法や、ガウシアンプロセスなどを用いてもよい。

次に、図３（ｂ）の学習結果算出処理について説明する。まず、しきい値算出部１０７によるしきい値算出処理（Ｓ３１１）について説明する。このしきい値は、異常検出部１０９に入力する異常測度と比較され、設備の正常／異常を判定するために用いられるものである。しきい値算出部１０７は、正常な学習データを異常と判定しないしきい値を算出する。言い換えれば、正常な学習データから得られる異常測度の最大値をしきい値として算出する。

あるいは、正常な学習データを予め定めた割合より多く正常と判定するしきい値を算出することにしてもよい。この場合は、正常な学習データから得られる異常測度をソートし、異常測度が低い方から前述の予め定めた割合に到達するところの異常測度をしきい値として採用する。

次に、二次元分布算出部１０８による二次元頻度分布算出処理（Ｓ３１３）について説明する。
図５は、学習時の二次元頻度分布算出処理のフローを示す図である。始めに、学習期間のセンサ信号を入力する（Ｓ５０１）。各センサ信号についてステップＳ５０３からＳ５０６までの処理を繰り返す（Ｓ５０２、ループ１）。まず、学習期間のデータの最大値（ＭＡＸ）と最小値（ＭＩＮ）を求める（Ｓ５０３）。次に、最小値から最大値までの範囲を指定された数Ｎで分割する際の刻み幅Ｓを算出する（Ｓ５０４）。なお、Ｓ＝（ＭＡＸ－ＭＩＮ）／Ｎで計算できる。次に、最小値から最大値までの範囲を外側に拡大し、二次元分布算出の処理範囲を算出する（Ｓ５０５）。拡大する範囲は、例えばＭＩＮをＭＩＮ－Ｓ×Ｍ、ＭＡＸをＭＡＸ＋Ｓ×Ｍに変更する。ここでＭは、予め決められた１以上の整数である。

次に学習期間の全データについて、次式によりセンサ信号値（Ｆ）からビン番号（ＢＮＯ）を算出する（Ｓ５０６）。
ＢＮＯ＝ＩＮＴ（Ｎ＊（Ｆ－ＭＩＮ）／（ＭＡＸ－ＭＩＮ））
ただし関数ＩＮＴ（Ｘ）はＸの整数部を表す。ビン番号（ＢＮＯ）を用いることで、各信号値は最小値０～最大値Ｎの（Ｎ＋１）段階の整数値に変換される。

次に、複数のセンサの中から２個のセンサを取り出し、それぞれのセンサ信号の組合せに基づき二次元分布を算出する。これを全てのセンサの組合せについて、ステップＳ５０８からＳ５１０までの処理を繰り返す（Ｓ５０７、ループ２）。ここで２個のセンサの組合せの中には同一センサの組合せを含める。従ってセンサの組合せ数（繰り返し数）は、（センサ数）×（センサ数＋１）／２となる。

まず、二次元分布算出用の二次元配列を確保し、全ての要素に０をセットする（Ｓ５０８）。配列のサイズはＮ＋２Ｍである。学習期間の全データについて、２個のセンサ信号のビン番号ＢＮＯに対応する配列の要素に１を加算する（Ｓ５０９）。すなわち、一方のセンサ信号のビン番号は列の要素に対応させ、他方のセンサ信号のビン番号は行の要素に対応させる。この処理により、センサ２個による信号の二次元の頻度分布（ヒストグラム）が算出される。この頻度分布を画像に変換して保存する（Ｓ５１０）。変換方法については後述する。図示はしていないが、二次元配列のサイズおよびステップＳ５０４およびＳ５０５で算出した各センサ信号の処理対象範囲と刻み幅を、学習結果蓄積部１０９に記録しておく。

ステップＳ５１０における、画像変換方法の例を説明する。始めに配列要素の最大値すなわち最大頻度を求める。画像サイズは配列サイズと同じとし、各要素の値から対応する座標の画素値を例えば、２５５×配列の要素値／最大頻度とする。数値２５５は画素値を８ビットで表す場合の最大値であり、この値を用いれば、そのままビットマップ形式で保存できる。あるいは、画素値を２５５×ＬＯＧ（配列の要素値＋１）／ＬＯＧ（最大頻度＋１）とする。ただし関数ＬＯＧ（Ｘ）はＸの対数を表す。このような変換式を用いれば、最大頻度が大きい場合も、非ゼロの頻度に非ゼロの画素値を対応させることが可能になる。

図６Ａ～図６Ｃは、二次元頻度分布画像のいくつかの例を示す図である。横軸にセンサａの信号値（ビン番号）を、縦軸にセンサｂの信号値（ビン番号）を示す。図５の処理により得られる画像は、二次元の特徴空間上で密度が高いところが高い画素値で表されているため、分布密度画像とも呼ぶことにする。ここでは、画素値の０を白、２５５を黒で表したグレイスケールの画像である。分布密度画像は２つのセンサの相関の強さに応じて、画像のパターンが変化する。図６Ａと図６Ｂは相関が強い場合で、特に図６Ａは時間的な相関が存在する場合である。

頻度分布画像の作り方は、上記方法に限定されない。例えば単純な頻度分布ではなく、１個のデータにガウス分布や他の重みつきフィルタを割り当て、それを重畳するようにしてもよい。あるいは、上記方法で得られた画像に所定サイズの最大値フィルタをかけたり、平均フィルタ、その他の重みつきフィルタをかけたりしてもよい。また、８ビットではなく、１６ビットに変換してもよい。また、必ずしも画像形式で保存する必要はなく二次元配列を変換せずにバイナリあるいはテキスト形式で保存してもよい。

図４Ａおよび図５の学習処理においては、学習結果蓄積部１０９に学習結果を保存しておく。学習結果として保存されるデータには、少なくとも特徴ベクトル抽出のためのパラメータ、異常測度算出のためのパラメータ、センサ正準化のためのパラメータ、抽出した全特徴ベクトルデータ、異常判定しきい値、二次元分布算出のためのパラメータ、二次元頻度分布が含まれる。特徴ベクトル抽出のためのパラメータ及び異常測度算出のためのパラメータは、学習時に指定されたものと共通である。センサ正準化のためのパラメータは、センサ信号入力部１０４がステップＳ４０２の処理で算出した各センサ信号の平均、標準偏差、最大値、最小値などである。二次元分布算出のためのパラメータは、二次元配列のサイズおよび二次元分布算出部１０８がステップＳ５０４およびＳ５０５の処理で算出した各センサ信号の処理対象範囲と刻み幅である。

次に、図３（ｃ）の異常検知時の異常判定処理について、図７から図１１を用いて説明する。
図７は、異常検出部１１０による異常検知処理（Ｓ３２１～Ｓ３２４）のフローを示す図である。ここでは、センサ信号蓄積部１０３に蓄積されたデータのうち指定された期間のデータ、または新たに観測されたデータについて、特徴ベクトルの抽出（特徴ベクトル抽出部１０５）、異常測度の算出（異常測度算出部１０６）を行い、これをしきい値（しきい値算出部１０７）と比較して、異常検出部１１０にて正常か異常かの判定を行う。

異常検出部１１０は、データベースから学習時に保存した学習結果を読み出す（Ｓ７０１）。その際、学習時の異常測度やしきい値に基づいて、ユーザが適切な処理番号を選択し、処理番号に対応付けられた学習結果を用いる。センサ信号入力部１０４は、センサ信号蓄積部１０３または設備１０１からセンサ信号１０２を入力し（Ｓ７０２）、センサ信号毎に正準化する（Ｓ７０３）。このとき、ステップＳ４０２の正準化の処理に用いたパラメータを用いる。次に、特徴ベクトル抽出部１０５は、選択したセンサ信号から、ステップＳ４０３の処理と同じ方法で特徴ベクトルの抽出を行う（Ｓ７０４）。

次に、全特徴ベクトルについてステップＳ７０６およびＳ７０７の処理を行う（Ｓ７０５、ループ）。異常測度算出部１０６は、注目ベクトルと学習データを用いて、異常測度を算出する（Ｓ７０６）。この処理は、図４のステップＳ４０７と同じ方法で行うが、学習データを全て用いることとする。異常検出部１１０は、ステップＳ７０１で読み出したしきい値とステップＳ７０６で算出した異常測度とを比較する。異常測度がしきい値以下であれば設備は「正常」と判定し、異常測度がしきい値より大きければ「異常」と判定する（Ｓ７０７）。

図８は、関連センサ特定部１１１による孤立度算出処理（Ｓ３２５～Ｓ３２７）のフローを示す図である。ここでは、異常が連続して検知された異常区間を抽出し、異常区間毎に学習データから乖離していることを示す孤立度ベクトルを算出し、これをもとに異常関連センサを特定する。

異常検出部１０９にて処理の対象としたセンサ信号を入力し（Ｓ８０１）、全てのデータについて、各センサ信号に対応するビン番号を算出する（Ｓ８０２）。算出の際には、ステップＳ７０１で読み出された二次元分布算出のパラメータ、具体的には二次元配列のサイズとセンサ毎の処理対象範囲と刻み幅を用い、ステップＳ５０６と同様の方法でビン番号を算出する。次に、異常検出部１０９にて算出した異常測度データを入力し（Ｓ８０３）、これをもとに異常が連続して検知されている異常区間を抽出する（Ｓ８０４）。異常区間を求める際には、予め指定された長さ以下の中断は連続しているものとみなす。逆に、日付が変わるなど予め決められたデータの切れ目では、異常検知が続いていても別の異常区間とする。

次に、各異常区間について、ステップＳ８０６からＳ８１３までの孤立度算出および関連センサ特定の処理を繰り返す（Ｓ８０５、ループ１）。まずセンサ２個の全ての組合せについての孤立度を０にリセットし（Ｓ８０６）。全てのセンサｉについて（Ｓ８０７、ループ２）、また全てのセンサｊについて（Ｓ８０８、ループ３）、孤立度を算出する。

各ループでは、着目する異常区間内の全データについて（Ｓ８０９、ループ４）、着目するセンサｉ，ｊの分布密度画像から、ステップＳ８０２で算出したセンサｉ，ｊそれぞれのビン番号に対応する座標の画素値を読み込む（Ｓ８１０）。画素値が０である場合、着目するセンサｉ，ｊの孤立度に１を加算する（Ｓ８１１）。この処理により、２個のセンサの各組合せについて着目する異常区間の孤立度が算出される。孤立度は、二次元分布上で対応する２個のセンサの信号値の組合せが学習データにない場合に高くなる。ここで、ループ４とループ２，３は逆の順番でも構わない。

２個のセンサの各組合せの孤立度をもとに、１個のセンサごとの孤立度を算出する（Ｓ８１２）。例えば、センサｉの孤立度は、センサｉを固定し、全てのセンサｊについてセンサｉ，ｊの孤立度を合計することにより算出する。各センサの孤立度を全センサ分まとめたものを孤立度ベクトルとする。

次に、孤立度をもとに異常関連センサを特定する（Ｓ８１３）。異常が検知されるのは評価対象のデータが学習データから乖離しているためであるから、異常検知された時刻のセンサ信号の孤立度が高くなるセンサを異常に関連するセンサとして抽出する。

上記の処理の中で、孤立度の算出（Ｓ８１０，Ｓ８１１）と異常関連センサの特定（Ｓ８１３）について詳細に説明する。
図９は、孤立度の算出法を具体的に説明する図である。（ａ）は異常検知時に例えば３個のセンサ１～３から入力する信号を示す。各信号の特徴ベクトルから異常測度を算出し、異常測度がしきい値を超えた時点を異常と判定する。（ｂ）は学習データから予め作成した分布密度画像で、２個のセンサの全ての組合せについて学習データの信号分布を画素値で示している。この場合、横軸をセンサｉの信号値、縦軸をセンサｊの信号値とすると、ｉ＝１～３、ｊ＝１～３の９通りの組合せの画像が存在する。

（ａ）で得られた異常検知時の各センサ１～３の信号値（丸印）を、（ｂ）の該当する分布密度画像の該当する座標位置にプロットする。プロット位置において画素値を読み込み、画素値が０であるとき、そのセンサの組合せの孤立度を１とする（×印）。画素値が０以外のときは、孤立度を０とする（丸印）。このようにして、センサごとに組合せの相手を変えて孤立度を合計することで、センサごとの孤立度を算出する。この例では、センサ１とセンサ３の組合せのみにおいて孤立しており、センサ１とセンサ３の孤立度は１、センサ２の孤立度は０となる。

図１０は、孤立度をもとに異常関連センサを特定する処理のフローを示す図である。このフローは、所定回数または条件を満たすまでの繰り返し処理である（Ｓ１００１、ループ１）。最初に、孤立度最大のセンサ、すなわちステップＳ８１２で算出された孤立度ベクトルの最大要素に対応するセンサを探索し、見つかったものをセンサＡとする（Ｓ１００２）。次に、センサＡと組合せて孤立度が最大となる他方のセンサを探索し、見つかったものをセンサＢとする（Ｓ１００３）。ここでの探索対象は、ステップＳ８０７のループ２の処理終了時に算出されているセンサ２個の組合せに対する孤立度である。見つかったセンサＡ，Ｂの孤立度の値をＩＳＯとする。

ＩＳＯが０より大きい場合（Ｓ１００４）、センサＡとＢを関連センサとして抽出する（Ｓ１００５）。また、センサＡとＢの分布密度画像に、異常データをプロットした画像を作成する（Ｓ１００６）。その際、正常データがグレイスケールで表されているのに対し、異常データは彩度の高い色で表す。また、正常データと重なりのない画素と重なりのある画素は異なる色で表す。以下の説明では、正常データのみの分布密度画像と区別するために、「異常プロット画像」と呼ぶこととする。なお、ステップＳ１００４でＩＳＯが０の場合は、ループ１を抜けて、関連センサ特定の処理を終了する。

次に、全てのセンサｉについて繰り返す（Ｓ１００７、ループ２）。センサｉの孤立度からセンサｉ，Ａの孤立度を差し引く（Ｓ１００８）。また、センサＡとＢが異なる場合は（Ｓ１００９）、センサｉの孤立度からセンサｉ，Ｂの孤立度も差し引く（Ｓ１０１０）。さらに、センサＡとＢの孤立度を０とする（Ｓ１０１１）。ステップＳ１００７からＳ１０１１までの処理は、Ｓ１００５で関連センサとして抽出したセンサＡ，Ｂの影響を取り除くために行う。以後、ループ１を繰り返し、残りの孤立度ベクトルのうちで孤立度最大となるセンサＡ，Ｂを新たな関連センサとして抽出する。このようにして、複数の異常関連センサを漏れなく抽出することができる。

次に、図３（ｃ）の重要度推定処理（Ｓ３２８）の詳細について説明する。ここでは診断情報提示部１１２により、異常検知時の最後にクラスタリングに基づく異常区間の重要度推定処理を行う。

図１１は、異常区間の重要度推定処理のフローを示す図である。最初に、全ての異常区間について（Ｓ１１０１、ループ１）、ステップＳ８１２で算出された孤立度ベクトルを要素の合計が１になるよう正規化する（Ｓ１１０２）。ただし、もともとの要素合計が０、すなわち全ての要素が０の場合はそのままとする。次に、全ての異常区間について（Ｓ１１０３、ループ２）、他の異常区間との間で正規化した孤立度ベクトルの類似度を算出する（Ｓ１１０４）。ここで、孤立度ベクトルの要素合計が０の場合は、全て類似度０とする。

類似度の定義としては、以下に示すように距離類似度、コサイン類似度、ヒストグラムインターセクションを利用することができる。以下の式において、Ｓはベクトルａ，ｂ間の類似度を表し、ａ_ｉ，ｂ_ｉはベクトルの要素を表す。

距離類似度は、２個のベクトル間のユークリッド距離を１から減算したものであり、（３）式で定義される。

コサイン類似度は、２個のベクトルのなす角度の余弦であり、（４）式で定義される。

ヒストグラムインターセクションは、頻度の合計が等しいヒストグラムの間の類似度を測る方法であり、２個のヒストグラムを重ねたときの共通部分の面積から算出される。縦軸が頻度でなくても要素数が等しくかつ要素値の合計が等しければ、ベクトル間の類似度算出にも適用可能である。ここではベクトルの要素値を用いて、（５）式で定義される。

いずれの類似度も０から１の値をとり、２個のベクトルが完全に等しいときに１となる。ここまでの処理により、全ての組合せの異常区間間の孤立度ベクトル同士の類似度が算出される。

次に、全ての異常区間を対象として、孤立度ベクトルの類似度に基づき複数のクラスタへの分類（クラスタリング）を行う（Ｓ１１０５）。クラスタ数を予め知ることができず、クラスタに属するデータ数に偏りがあることが予想されるため、ｋ平均法のような分割統治タイプではなく、階層的クラスタリングを適用するとよい。階層的クラスタリングは、個々のデータを１個ずつのクラスタに割り当てるところから開始し、類似したクラスタを再帰的に結合していくものである。結合するクラスタを選択する基準によって、最短距離法、最長距離法、群平均法などの手法がある。それぞれの方法において、クラスタ間の類似度は、クラスタをまたがる区間どうしの類似度の最小値、最大値、平均値で定義される。クラスタ間の類似度が大きいものから順次結合していき、全てのクラスタ間の類似度が予め定めた基準値を下回ったとき、結合を停止する。基準値は０から１の間の実数とするが、０とすると全区間が１個のクラスタに結合され、１とすると全区間が全て異なるクラスタに分かれる。

次に、各クラスタを異常区間の頻度が高い順、すなわちクラスタに含まれる区間数が多い順にソートする（Ｓ１１０６）。これによりクラスタの順位を決定し、順位が高いほど重要度が高いと推定する。これは、異常あるいはその予兆は同じ現象が続くものであり、孤立度ベクトルが類似している区間は同じ現象が発生しているという考えに基づく。

次に、各クラスタについて（Ｓ１１０７、ループ３）、異常測度の最大値または区間内の異常測度累積値が大きい順に、異常区間をソートする（Ｓ１１０８）。そして、１位の異常区間を代表区間とする。これは、クラスタ内では同じ現象が発生していると考えるため、各クラスタで１個の区間の現象を確認すれば十分であり、異常測度が大きいほど確認が容易であると考えられるからである。

以上の処理結果に基づき、重要度が高いと推定した関連センサや異常区間を優先して、診断支援情報を提示する（Ｓ１１０９）。具体的には、提示する重要度の高い関連センサは、ステップＳ１１０６のソートで１位となったクラスタすなわち最も頻度の高いクラスタに含まれ、ステップＳ１１０８のソートで１位となった区間において、図１０に示す関連センサ特定処理にて最初のループのステップＳ１００２で見つかったセンサとする。なお、提示の順序および内容の詳細については、ＧＵＩ（Graphical User Interface）の例とともに後述する。

上記ステップＳ１１０６では、クラスタに含まれる区間数が多いほど重要度が高いと推定したが、クラスタに含まれる区間の長さ（累積時間）が長いほど重要性が高いともいえるので、区間の長さに着目して順位を決定しても良い。

次に、以上の動作を実現するための異常検知装置１００のユーザインタフェース（ＧＵＩ）の例を説明する。
図１２Ａは、オフライン解析実施のための学習期間、及び処理パラメータ含む解析条件を設定するＧＵＩの例である。この画面では、算出された学習結果をレシピとして登録することも可能である。また、過去のセンサ信号１０２は、設備ＩＤ及び時刻と対応付けられてデータベースに保存されているものとする。

オフライン解析条件設定画面１２０１では、対象設備、学習期間、テスト期間、使用センサ、異常測度算出パラメータ、二次元分布算出パラメータ、異常区間抽出パラメータ、クラスタリングパラメータを入力する。設備ＩＤ入力ウィンドウ１２０２には、対象とする設備のＩＤを入力する。設備リスト表示ボタン１２０３の押下により、センサ信号蓄積部１０３に保存されているデータの装置ＩＤのリストが表示されるので、リストから選択入力する。異常検知装置１００につながる設備１０１が１台のみの場合は、設備ＩＤ入力ウィンドウ１２０２は表示されない。

学習期間入力ウィンドウ１２０４には、学習データを抽出したい期間の開始日と終了日を入力する。テスト期間入力ウィンドウ１２０５には、解析対象としたい期間の開始日と終了日を入力する。センサ選択ウィンドウ１２０６には、使用するセンサを入力する。リスト表示ボタン１２０７のクリックによりセンサリスト１２０８が表示されるので、リストから選択入力する。リストから複数選択することも可能である。除外するセンサを指定するようにしてもよい。センサ信号入力時は、ここで選択されたセンサの情報のみが入力される。

異常測度算出パラメータ入力ウィンドウ１２０９には、異常測度算出において使用するパラメータを入力する。図は手法として局所部分空間を採用した場合の例であり、近傍ベクトル数と正則化パラメータを入力する。正則化パラメータは、（２）式において相関行列Ｃの逆行列が求められないことを防ぐため、対角成分に加算する小さい数である。二次元分布算出パラメータ入力ウィンドウ１２１０には、二次元分布算出において使用するパラメータとして、二次元配列のサイズすなわち作成する画像サイズとそのうちの正常範囲に対応するサイズの情報を入力する。異常区間抽出パラメータ入力ウィンドウ１２１１には、異常区間抽出において異常検知が連続しているとみなす中断（正常判定）期間の最大長さを入力する。クラスタリングパラメータ入力ウィンドウ１２１２には、異常区間の階層的クラスタリングにおいて、クラスタの結合を停止する基準となる類似度を入力する。

以上の解析条件の情報が確定したら、実行ボタン１２１４の押下により、オフライン解析を実行する。
まず、学習期間のセンサ信号を用い、図４Ａの処理フロー、続いてステップＳ３１１のしきい値算出処理、続いて図５の処理フローに従って学習を実行する。学習結果として、ステップＳ４０２で算出されたセンサ信号毎の平均と標準偏差、ステップＳ４０３で抽出された学習期間の全特徴ベクトルデータ、ステップＳ３１１で算出されたしきい値、ステップＳ５０４およびＳ５０５で算出した各センサ信号の処理対象範囲と刻み幅、ステップＳ５１０で作成された分布密度画像を保存しておく。

さらに、学習期間およびテスト期間のセンサ信号を用い、図７の処理フローに従って異常測度を算出し、正常か異常かの判定を行い、判定結果を異常測度およびしきい値と併せて保存しておく。ただし、学習期間のデータについては、ステップＳ４０７で算出した異常測度を用いて、正常か異常かの判定を行う。

次に、図８および図１０の処理フローに従って、異常区間を抽出して区間毎に孤立度を算出し、関連センサを特定する。表示のため、各異常区間の開始時刻、終了時刻、孤立度、特定された関連センサ名と異常プロット画像を保存しておく。

さらに、図１１の処理フローに従って、異常区間をクラスタリングして重要度を推定する。処理結果として、各クラスタに含まれる区間番号を保存しておく。保存の際には、ステップＳ１１０８で付けられたクラスタの順位、ステップＳ１１０９で付けられた、クラスタ内の異常区間の順位に従う。

解析終了後、後述する結果表示画面が表示される。ユーザによる確認が終了すると、オフライン解析条件設定画面１２０１に戻ってくる。レシピ名入力ウィンドウ１２１３にレシピ名を入力し、登録ボタン１２１５を押下することにより、設備ＩＤ及びレシピ名と対応付けて学習結果および解析結果を保存し、終了する。ここで、学習結果には、学習の実行により作成保存されたデータのほか、入力ウィンドウ１２０９～１２１２で入力された異常測度算出パラメータ、二次元分布算出パラメータ、異常区間抽出パラメータ、クラスタリングパラメータが含まれる。終了ボタン１２１６が押下された場合は、何もしないで終了する。この場合、学習により作成保存された学習結果および、続く異常検知処理により作成保存された解析結果は、削除されるか次に実行される解析によって上書きされる。

登録された学習結果は、活性か不活性かのラベルをつけて管理され、以降オンラインの解析が実行される。オンライン解析では、新しく入力されたデータに対し、装置ＩＤが一致する活性な学習結果の情報を用いて、図７から図１０に示す処理を行い、結果をレシピ名および処理日時と対応付けて保存しておく。これらの処理は定期的、例えば１日毎に実行する。サンプリング間隔が短い設備やリアルタイム性を求められる設備については、実行の間隔をもっと短くする。続いて図１１に示す処理を行う。この処理も同じタイミングで実行するが、最新の入力データのみを対象とするのではなく、過去に検知された異常区間例えば１か月分を併せてクラスタリングを行い、結果をレシピ名および処理日時と対応付けて保存しておく。

図１２Ｂは、オンライン解析結果の表示対象を指定するためのＧＵＩの例である。ユーザは、表示対象指定画面１２２１から表示対象の設備、レシピ及び期間を指定する。始めに、装置ＩＤ選択ウィンドウ１２２２により設備ＩＤを選択する。次に、レシピ名選択ウィンドウ１２２３により、設備ＩＤ１２２２を対象としたレシピのリストから表示対象のレシピを選択する。データ記録期間表示部１２２４には、入力されたレシピを用いて処理され、記録が残されている期間の開始日と終了日が表示される。結果表示期間指定ウィンドウ１２２５には、結果を表示したい期間の開始日と終了日を入力する。表示ボタン１２２６を押下すると、異常検知処理の結果が表示される。終了ボタン１２２７を押下すると、表示対象を指定する処理を終了する。

図１３Ａ～図１３Ｄは、解析結果をユーザに示すためのＧＵＩの例である。ユーザが各画面の上部に表示されたタブを選択することにより、解析結果全体表示画面１３０１、解析結果拡大表示画面１３０２、クラスタ表示画面１３０３および重要度診断画面１３０４のいずれかに切り換わる。

図１３Ａは、解析結果全体表示画面１３０１の例である。解析結果全体表示画面１３０１には、指定された期間の、異常測度、しきい値、及び判定結果、並びにセンサ信号の時系列グラフが表示される。期間表示ウィンドウ１３０５には、オフライン解析の結果を表示する場合は図１２Ａで指定された学習期間及びテスト期間が表示される。オンライン解析の結果を表示する場合は、図示していないが、図１２Ｂで指定された結果表示期間が表示される。

異常測度表示ウィンドウ１３０６には、指定された学習期間・テスト期間あるいは結果表示期間での異常測度１３０６ａ、しきい値１３０６ｂ（破線）、及び判定結果１３０６ｃが表示される。また、学習に使用した区間に丸印１３０６ｄが表示される。センサ信号表示ウィンドウ１３０７には、指定された学習期間・テスト期間あるいは結果表示期間での指定されたセンサについて、時系列センサ信号１３０７ａが表示される。

センサ選択ウィンドウ１３０８では、ユーザの入力によってセンサを指定する。ただし、ユーザが指定する前は、図１１に示す重要度推定処理で求められたセンサが選択されている。具体的には、ステップＳ１１０６のソートで１位となったクラスタすなわち最も頻度の高いクラスタに含まれ、ステップＳ１１０８のソートで１位となった区間について、図１０に示す関連センサ特定処理において最初のループのステップＳ１００２で見つかったセンサとする。カーソル１３０９は、拡大表示の時の起点を表し、ユーザのマウス操作により移動できる。表示日数指定ウィンドウ１３１０には、解析結果拡大表示画面１３０２での拡大表示の起点から終点までの日数が表示され、この画面で入力することもできる。日付表示ウィンドウ１３１１には、カーソル位置の日付が表示される。終了ボタン１３１２の押下により、解析結果全体表示画面１３０１、解析結果拡大表示画面１３０２、クラスタ表示画面１３０３および重要度診断画面１３０４のいずれもが消去され、解析結果の表示が終了する。

図１３Ｂは、解析結果拡大表示画面１３０２の例である。解析結果拡大表示画面１３０２には、解析結果全体表示画面１３０１においてカーソル１３０９で示された日付を起点とし、表示日数指定ウィンドウ１３１０で指定された日数の期間内の、異常測度、しきい値、判定結果、及びセンサ信号の時系列グラフが表示される。すなわち、異常測度表示ウィンドウ１３０６及びセンサ信号表示ウィンドウ１３０７には、解析結果全体表示画面１３０１と同様の情報が、拡大して表示される。

なお、解析結果拡大表示画面１３０２では、スクロールバー１３１３とスクロールバー領域１３１４を追加表示している。スクロールバー１３１３の長さは表示日数指定ウィンドウ１３１０で指定された日数に、スクロールバー領域１３１４の左端部が拡大表示の起点に対応する。ユーザはスクロールバー１３１３を操作することで、表示の起点を変更することも可能であり、この変更はカーソル１３０９の位置と日付表示ウィンドウ１３１１の表示に反映される。スクロールバー領域１３１４の全体の長さは解析結果全体表示画面１３０１に表示されている期間に相当する。

図１３Ｃは、図１１に示す重要度推定処理のクラスタリングの結果を表示するクラスタ表示画面１３０３の例である。クラスタ表示画面１３０３は、クラスタ時系列情報表示ウィンドウ１３３１、異常測度の時系列グラフ１３３６、孤立度ベクトル表示ウィンドウ１３３７、および終了ボタン１３１２で構成される。終了ボタン１３１２が押下されたときの動作は、他の画面と同様である。

クラスタ時系列情報表示ウィンドウ１３３１には、日付別クラスタ情報が表示され、凡例１３３２、クラスタ順位１３３３、日付別クラスタ有無情報１３３４、診断結果１３３５から構成される。クラスタ順位１３３３は、ステップＳ１１０６のソートの結果に従うものであり、１個の区間のみのクラスタまで、全てを上から順に記載する。凡例１３３２は、クラスタを区別するために使用する枠線の種類（色、太さ、スタイル）である。日付別クラスタ有無情報１３３４は、解析結果全体表示画面１３０１の表示期間に対応しており、各クラスタに含まれる区間がその日にある場合は黒、ない場合は白で表す。有無のみでなく区間数を表示してもよい。

診断結果１３３５は、クラスタ毎に異常か対策不要か誤報かを診断した結果であり、重要度診断画面１３０４でユーザによって入力されるものである。ユーザによる入力が未実施の場合は、空欄にするか、プログラムによる推定結果を網掛けなどで区別して表示する。なお、推定方法については後で説明する。日付別クラスタ有無情報１３３４の日付位置に合わせて、異常測度の時系列グラフ１３３６が表示される。これは、解析結果全体表示画面１３０１の異常測度表示ウィンドウ１３０６に表示される内容と同じものである。

孤立度ベクトル表示ウィンドウ１３３７には、表示期間中の全ての異常区間について、Ｓ８１２で算出された孤立度ベクトルを表す棒グラフが時系列順に表示される。グラフは横軸をセンサ信号名、縦軸を孤立度として描画される。それぞれのグラフは、その区間が含まれるクラスタに応じて、凡例１３３２に表示された種類の枠線が付けられる。また、日付の情報が付加されて表示される。図示していないが、時刻情報、異常測度や累積異常測度などのより詳細な情報を併せて表示してもよい。

図１３Ｄは、図１１に示す診断支援情報提示（Ｓ１１０９）にかかる重要度診断画面１３０４の例である。選択されたクラスタの選択された異常区間について、孤立度ベクトル、異常プロット画像、異常測度・しきい値・判定結果の時系列グラフ、異常関連センサ信号の時系列グラフが表示される。ユーザによる操作がなされていない初期状態では、１位のクラスタの１位の異常区間が選択されている。つまり、重要度が最も高いと推定されるクラスタの代表区間の情報が表示されている。

クラスタ情報表示ウィンドウ１３５１には、クラスタ選択ウィンドウ１３５２と、選択中のクラスタに含まれる区間数と、診断結果入力ウィンドウ１３５３と確認チェックボタン１３５４が表示される。異常区間情報表示ウィンドウ１３５５には、区間選択ウィンドウ１３５６と、選択中の区間の区間番号と日付と時刻が表示される。孤立度ベクトル表示ウィンドウ１３５７には、選択中の区間の孤立度ベクトルを表す棒グラフが表示される。

異常プロット画像表示ウィンドウ１３５８には、選択中の画像が表示される。この画像は、図１０に示す処理フローのステップＳ１００６で作成された、異常プロット画像である。正常データの分布１３６０がグレイスケールで表され、異常データの分布１３６１が彩度の高い色（たとえば赤）で表される。表示する画像は、画像選択ウィンドウ１３５９により選択することが可能である。この番号は、図１０に示す処理フローのループ１において何回目で作成されたのかを表す。初期状態では１回目に作成された画像が選択されている。画像選択ウィンドウ１３５９において数値を入力あるいは矢印ボタンを使って前または後ろの番号を選ぶと、異常プロット画像および関連センサ信号の時系列グラフが更新されるので、関連センサが複数ある場合にも確認することができる。

異常測度表示ウィンドウ１３６２には、選択中の異常区間を含む期間、例えば１日分の異常測度、しきい値、判定結果の時系列グラフが表示される。併せて、異常区間の時刻範囲を表すバー１３６３が表示される。関連センサ（Ａ）信号表示ウィンドウ１３６４と関連センサ（Ｂ）信号表示ウィンドウ１３６５には、選択中の異常プロット画像に対応して、それぞれステップＳ１００２で特定された関連センサＡとステップＳ１００３で特定された関連センサＢのセンサ信号の時系列グラフが表示される。表示期間は異常測度グラフ表示と同じとする。

この重要度診断画面１３０４により、ユーザは検知された異常の重要度診断を行うことができる。同じクラスタに含まれる区間は同じ現象が発生しているという考えから、診断はクラスタ毎に行う。区間選択ウィンドウ１３５６への入力により、異常区間情報表示ウィンドウ１３５５の表示内容とともに、孤立度ベクトル、異常プロット画像、異常測度時系列グラフ、関連センサ信号時系列グラフが全て更新される。この操作により、本当に同じ現象が発生しているかどうかを確認することができるが、通常は代表区間だけ確認すれば十分である。確認ができたら、診断結果入力ウィンドウ１３５３で異常か対策不要か誤報か空欄かをリストから選択し、確認チェックボタン１３５４にチェックを入れる。この操作の結果は、クラスタ時系列情報表示ウィンドウ１３３１の診断結果１３３５に反映される。診断結果入力ウィンドウ１３５３には、初期状態では、診断結果１３３５と同様に空欄またはプログラムによる推定結果が表示されている。クラスタ選択ウィンドウ１３５２の矢印ボタンにより、選択クラスタを上位のものから下位のものに切り替えて、順に診断していくとよい。全て診断結果を入力、あるいは途中であっても、診断確定ボタン１３６６の押下により、診断結果が保存される。

以上の図１３Ａから図１３Ｄに示した例では、クラスタ情報表示ウィンドウ１３５１から、異常測度のグラフ上で目立つ異常のほとんどが同じクラスタに含まれており、その孤立度ベクトル表示ウィンドウ１３３７，１３５７から、異常関連センサは左から９番目のセンサ（Ｓ９）であることが想像できる。次に重要度診断画面１３０４で、異常プロット画像１３５８から、Ｓｅｎｓｏｒ＿Ｙの異常のデータが正常データより高いところに分布していることが分かり、関連センサ信号の波形からもそのことを確認できる。また、解析結果全体表示画面１３０１および解析結果拡大表示画面１３０２のセンサ信号表示ウィンドウ１３０７には、初期状態にはこのセンサ信号が表示される。異常区間で信号値が高くなっていることが、容易に確認できる。このように、本実施例によれば、解析直後の初期状態の画面表示を確認するのみで、主要な異常現象を確認することが可能である。また、表示期間で複数の現象があっても、クラスタリングにより異なるクラスタに分かれることが期待され、別々に診断することが可能である。さらに、低頻度なものは確認作業を省略しても設備を保守する上での危険は小さいと判断し、途中でやめることもできる。逆に、全て確認して誤報や対策不要の発生状況をみることにより、学習期間が適正であるかのチェックに利用することも可能である。

上記実施例において、診断結果入力ウィンドウ１３５３に初期状態で推定結果を表示する場合の、推定方法について説明する。図１１に示す重要度推定処理によれば、重要度が最高のものから最低のものまで順に並べることが可能であるものの、実際に異常（予兆を含む）なのか誤報なのかの推定は行っていないため、このままでは推定結果を表示できない。そこで、重要度順を参考にして、他の特徴も加えて異常か対策不要か誤報かを推定する。

図１４は、対策不要と誤報を推定する処理のフローを示す図である。各クラスタについて以下の推定処理を行う（Ｓ１４０１、ループ）。まず、１位のクラスタで、かつ含まれる区間数が予め指定した基準値（基準１とする）より大きいかどうかを判定する（Ｓ１４０２）。条件を満足すれば（Ｙｅｓ）、対策不要の推定処理を行う。そのクラスタの代表区間の孤立度最大のセンサについて、一定期間毎、例えば１日毎に「統計的特徴」を算出する（Ｓ１４０３）。ここに「統計的特徴」とは、関連センサ信号の平均値、分散値、最大値、最小値、中央値などである。また、頻度分布をもとに最頻値を求めたり、モード分割してモード別の平均、分散を求めたりしてもよい。算出した統計的特徴について、学習期間と１位のクラスタを含む期間を比較し（Ｓ１４０４）、「ステップ状」に変化していれば（Ｓ１４０５）、対策不要と推定する（Ｓ１４０６）。そうでなければ異常と推定する（Ｓ１４０７）。「ステップ状」とは、比較する期間の間では顕著な違いがあり、かつ各期間の中では違いが小さいこと、さらに異常検知した期間の中で上昇も下降もしていないことを意味する。

次に、ステップＳ１４０２の判定で条件にあてはまらない場合（Ｎｏ）は、誤報推定処理を行う。まず、クラスタに含まれる区間数が予め指定した基準値（基準２、ただし基準２＜基準１）より小さいかどうかを判定する（Ｓ１４０８）。基準２より小さい場合（Ｙｅｓ）はステップＳ１４０９へ進み、そうでない場合（Ｎｏ）は異常と推定する（Ｓ１４０７）。これは、同じ現象が何回も起こったのであれば実際に異常である可能性が高いという考えに基づく。次に、代表区間のセンサ信号、異常測度、孤立度から特徴を抽出する（Ｓ１４０９）。ここで抽出する特徴は、異常測度最大値、累積値、区間長さ、孤立度、関連センサ信号の変化速度、変化方向、異常プロット画像上での正常データとの位置関係を表す特徴などである。これら特徴が予め指定された「誤報条件」に適合するかを判定し（Ｓ１４１０）、適合すれば誤報と推定し（Ｓ１４１１）、適合しなければ異常と推定する（Ｓ１４０７）。

ここで「誤報条件」については、異常区間が単発で存在し、異常測度が低く、異常区間長が短い場合は誤報である確率が高いことが経験的に分かっている。従って、異常区間の数、異常区間の長さ、異常測度のいずれか、または全てが基準値以下であれば、異常ではなく誤報であると推定する。さらに、例えば設備の起動やシャットダウンなど状態遷移中は、学習データが少ないため誤報が発生しやすいことも経験的に分かっている。従って、関連センサ信号の変化速度と変化方向から状態遷移中であるかどうかを判断し、遷移中であれば誤報であると推定する。

上記した対策不要の推定処理（Ｓ１４０３～Ｓ１４０６）の具体例について説明する。この処理は、異常ではない人為的な作業、例えばメンテナンスや条件設定変更などによる変化を異常として検知したケースを、対策不要と推定することである。このような場合、作業の時点で急激に変化し、その後は安定した状態となることから、上記の処理により推定が可能である。

図１５は、対策不要と推定すべき例を示す図である。上段は異常測度表示ウィンドウ１３０６で（図１３Ａまたは図１３Ｂ参照）、異常測度１５０１、しきい値１５０２および判定結果１５０３を時系列で示している。（Ａ）点以降の後半には異常区間が多数検出されており、図示していないが、そのほとんどが１位のクラスタに含まれる。下段はセンサ信号表示ウィンドウ１３０７で、１位のクラスタの代表区間の孤立度最大のセンサの時系列信号である。センサ信号１５０４とともに、１日毎の平均値１５０５を重ねて示している。いずれも、（Ａ）点でステップ状に変化しその後は安定していることから、対策不要と推定する。

以下、上記した実施例の変形例を説明する。
上記実施例では１台の設備を対象としたが、複数の設備に適用すれば、複数の設備の異常監視を効率的に行うことが可能となる。その場合は、各設備について別々に図４Ａおよび図５に示す処理により学習を実施し、学習データを蓄積しておく。さらに、各設備のセンサ信号を継続的に入力、蓄積し、定期的に図７、図８、図１０に示す処理を行って異常を検知し、異常区間毎に孤立度を算出し関連センサを特定しておく。同じタイミングで、過去に検出した異常区間と新しく検出した異常区間を合わせて図１１に示す処理により、クラスタリングを行い、クラスタの順位を付ける。そして、各設備の１位のクラスタに含まれる異常区間がある日数を数え、その日数が多い順に設備の順位をつける。

図１６は、複数の設備の異常監視結果を一覧するＧＵＩの例である。設備別クラスタ情報表示ウィンドウ１６０１には、縦軸を設備、横軸を日付とし、１位のクラスタに含まれる異常区間がある日を黒で表している。設備は、前述の順位に従って上から並べている。設備名をクリックすることにより、図１３Ａに示す解析結果全体表示画面が表示されるようにする。あるいは、セルをクリックすることにより、対応する異常区間について、図１３Ｄに示す重要度診断画面が表示されるようにする。従来は異常検出の有無を表示し、連続して発生している場合に注意を喚起することが可能であったが、誤報がある場合や、特定の異常が断続的に発生する場合に重要度の推定を誤る可能性が高かった。これに対し本実施例では、クラスタリングを利用することにより、重要度の高い異常が発生している設備を優先的に調べることができるため、確認作業の効率を向上させることができる。

また、上記実施例では、関連センサ特定部１１１は、異常区間毎に、二次元分布算出部１０８で算出された二次元頻度分布を用いて孤立度ベクトルを算出し、これに基づく異常関連センサを特定した。そして診断情報提示部１１２は、異常区間を孤立度ベクトルの類似度に基づいてクラスタリングした。このように、異常区間や関連センサの重要度を判断するためのパラメータとして孤立度を用いたが、これ以外の方法も可能である。

孤立度以外のパラメータとして、「寄与度ベクトル」を用いる方法を説明する。寄与度ベクトルは、異常測度算出部１０６にて算出する異常測度へどれだけ影響しているかを示す指標であり、異常測度を算出する際に同時に求めることができる。異常測度は、正常を表す基準データから観測データへの距離で定義される。一方寄与度ベクトルは、基準データから観測データへのベクトルであり、異常測度算出に局所部分空間法を採用する場合は、図４Ｂのベクトル（ｑ－Ｘｂ）で定義される。よって、寄与度ベクトルが大きいほど異常測度が大きく、正常なデータから乖離していることになり、重要度が高いものと言える。寄与度ベクトルの類似度算出の際には、要素の合計が１になるように正規化しておく。関連センサ特定においては、寄与度ベクトルの要素の値が大きい順に異常との関連が強いセンサであると見なすことができる。なお、寄与度ベクトルを用いる場合は、二次元頻度分布や孤立度を算出する必要がない。このように、孤立度に代えて寄与度に着目しても、同様の効果を得ることができる。

１００：異常検知装置、１０１：設備、１０２：センサ信号、１０３：センサ信号蓄積部、１０４：センサ信号入力部、１０５：特徴ベクトル抽出部、１０６：異常測度算出部、１０７：しきい値算出部、１０８：二次元分布算出部、１０９：学習データ蓄積部、１１０：異常検出部、１１１：関連センサ特定部、１１２：診断情報提示部、１２０１：オフライン解析条件設定画面、１２２１：表示対象指定画面、１３０１：解析結果全体表示画面、１３０２：解析結果拡大表示画面、１３０３：クラスタ表示画面、１３０４：重要度診断画面。

Claims

設備に装着された複数のセンサから出力される複数の時系列のセンサ信号を入力するセンサ信号入力部と、
前記センサ信号から時刻毎に特徴ベクトルを抽出する特徴ベクトル抽出部と、
指定された学習期間の前記特徴ベクトルを学習データとして用いて、指定された期間のデータの特徴ベクトル、または、新たに観測されたデータの特徴ベクトルから各時刻の異常測度を算出する異常測度算出部と、
前記異常測度をしきい値と比較することにより各時刻のセンサ信号が正常か異常かを判定する異常検出部と、
異常が連続して検知される異常区間を抽出し、前記異常区間毎に前記学習データからの乖離もしくは前記異常測度への影響を示すパラメータを算出し、前記パラメータに基づき異常検知の原因となった関連センサを特定する関連センサ特定部と、
前記異常区間を前記パラメータの類似度に基づいて複数のクラスタに分類し、含まれる異常区間数が多いクラスタほど重要度が高いと推定して、診断支援情報に順位を付けて提示する診断情報提示部と、
を備えることを特徴とする異常検知装置。
請求項１記載の異常検知装置であって、
前記関連センサ特定部が算出する前記学習データからの乖離を示す前記パラメータは、２つのセンサ信号の二次元分布に基づいて算出される孤立度であり、
前記診断情報提示部は前記診断支援情報として、前記クラスタに含まれる異常区間の時刻情報と、孤立度ベクトルを表すグラフと、前記関連センサの二次元分布に前記異常区間のデータをプロットした画像と、前記関連センサの時系列波形と、を表示することを特徴とする異常検知装置。
請求項２記載の異常検知装置であって、
前記診断情報提示部は、前記クラスタに含まれる異常区間の数が予め指定した基準値より大きい場合、前記クラスタに含まれる異常区間の関連センサ信号に基づき一定期間毎に統計的特徴を算出し、前記統計的特徴が学習期間と比較してステップ状に変化している場合は対策不要であると推定する、ことを特徴とする異常検知装置。
請求項２記載の異常検知装置であって、
前記診断情報提示部は、前記クラスタに含まれる異常区間の数、異常区間の長さ、異常測度のいずれか、または全てが予め指定した基準値以下であれば、異常ではなく誤報であると推定する、ことを特徴とする異常検知装置。
請求項１記載の異常検知装置であって、
前記関連センサ特定部が算出する前記異常測度への影響を示す前記パラメータは、前記異常測度算出部にて算出される寄与度であり、正常を表す基準データから観測データへのベクトルとして求められることを特徴とする異常検知装置。
請求項１記載の異常検知装置であって、
前記設備として複数の設備を対象とし、
前記異常検出部は、前記設備ごとに、オンラインで新しく取得したセンサ信号について正常か異常かを判定し、
前記診断情報提示部は、過去に検出した異常区間と新しく検出した異常区間を合わせてクラスタに分類し、含まれる異常区間数が最大のクラスタを求め、前記クラスタに含まれる異常区間がある日数を数え、
前記診断支援情報として、縦軸に設備を前記日数が多い順に並べ、横軸を日付として前記クラスタに含まれる異常区間の有無を表示することを特徴とする異常検知装置。
複数の時系列センサ信号を入力して時刻毎に特徴ベクトルを抽出し、
指定された学習期間の前記特徴ベクトルを学習データとして用いて、指定された期間のデータの特徴ベクトル、または、新たに観測されたデータの特徴ベクトルから各時刻の異常測度を算出し、
前記異常測度に基づいて各時刻のデータが異常か正常かを判定する異常検知方法であって、
異常が連続して検知される異常区間を抽出し、
前記異常区間毎に前記学習データからの乖離もしくは前記異常測度への影響を示すパラメータを算出し、
前記パラメータに基づき異常検知の原因となった関連センサを特定し、
前記異常区間を前記パラメータの類似度に基づいて複数のクラスタに分類し、
含まれる異常区間数が多いクラスタほど重要度が高いと推定して、診断支援情報に順位を付けて提示することを特徴とする異常検知方法。