(57)【要約】【課題】  多次元時系列センサ信号に基づく異常検知において、対策、調査など次のアクションを決めるために、どのセンサが異常に関連しているのかを特定する必要があるが、基準値と観測値の差が大きいものを関連センサとする従来の方法では、誤る場合があった。【解決手段】  センサ信号から時刻毎に多次元特徴ベクトルを抽出し、予め指定された学習期間の特徴ベクトルの集合と各時刻の特徴ベクトルに基づいて各時刻の基準特徴ベクトルを算出し、各時刻の特徴ベクトルと基準特徴ベクトルの差に基づき異常測度を算出し、異常測度と所定のしきい値との比較により異常を検出し、異常が検出された時刻について、特徴値の２次元の分布密度に基づいて異常に関連するセンサを特定するようにした。【選択図】図１
2014142697.tif 000002
【特許請求の範囲】
【請求項１】
設備または装置に装着された複数のセンサから出力される複数のセンサ信号に基づいて前記設備または装置の異常を診断する方法であって、予め指定された学習期間における前記複数のセンサ信号の特徴値から多次元の特徴ベクトルを求め、前記診断時に前記複数のセンサ信号の特徴値から時刻毎に多次元特徴ベクトルを抽出し、前記求めた予め指定された学習期間の前記多次元の特徴ベクトルの集合と前記診断時に抽出した時刻毎の前記多次元の特徴ベクトルに基づいて前記時刻毎の基準特徴ベクトルを算出し、前記抽出した時刻毎の多次元特徴ベクトルと前記算出した時刻毎の基準特徴ベクトルとの差に基づき異常測度を算出し、前記算出した異常測度と所定のしきい値とを比較することにより異常を検出し、前記異常が検出された時刻における前記複数のセンサ信号の特徴値の２次元の分布密度に基づいて前記検出した異常に関連するセンサを前記複数のセンサの中から特定することを特徴とする異常診断方法。
【請求項２】
上記異常診断方法において、前記所定のしきい値を予め指定された学習期間の異常測度に基づき算出することを特徴とする請求項１に記載の異常診断方法。
【請求項３】
設備または装置に装着された複数のセンサから出力される複数のセンサ信号をもとに学習データを作成して蓄積する工程と、前記設備または装置に装着された複数のセンサから新たに出力された複数のセンサ信号の異常診断を行う工程とを含む設備の異常を診断する方法であって、
前記学習データを作成して蓄積する工程は、前記複数のセンサ信号から特徴ベクトルを抽出する工程と、前記抽出した特徴ベクトルを学習データとして蓄積する工程と、前記学習データとして蓄積した特徴ベクトルの各々について前記特徴ベクトルに応じて前記蓄積した学習データの中から所定数の特徴ベクトルを選択する工程と、前記選択された所定数の特徴ベクトルを用いて学習用の基準ベクトルを作成する工程と、前記作成した学習用の基準ベクトルに基づいて前記学習データとして蓄積した特徴ベクトルの各々について異常測度を算出する工程と、前記算出した前記異常測度に基づきしきい値を算出する工程と、前記特徴ベクトルから全組合せの２次元の特徴分布密度を算出する工程とを有し、
前記センサ信号の異常診断を行う工程は、前記設備または装置に装着された前記複数のセンサから出力される複数のセンサ信号から特徴ベクトルを抽出して観測ベクトルとする工程と、前記抽出した観測ベクトルに応じて前記蓄積した学習データの中から所定数の特徴ベクトルを選択する工程と、前記学習データの中から選択された所定数の特徴ベクトルを用いて異常診断用の基準ベクトルを作成する工程と、前記作成した異常診断用の基準ベクトルに基づいて前記観測ベクトルの異常測度を算出する工程と、前記算出した異常測度と前記学習データを作成して蓄積する工程において算出したしきい値に基づき前記観測ベクトルが異常か正常かを判定する工程と、前記判定する工程において異常と判定された観測ベクトルに対応するセンサ信号が検出された時刻における前記観測ベクトルと前記２次元の特徴分布密度に基づき異常に関連するセンサを特定する工程とを有する
ことを特徴とする異常診断方法。
【請求項４】
前記学習用の基準ベクトルを作成する工程と、前記異常診断用の基準ベクトルを作成する工程とは、何れも局所部分空間法を利用すること特徴とする請求項３に記載の異常診断方法。
【請求項５】
前記学習データを作成して蓄積する工程において算出される前記２次元の特徴分布密度は、前記特徴ベクトルの集合の２次元の特徴の頻度分布に基づき算出され、濃淡画像に変換して保存されることを特徴とする請求項３に記載の異常診断方法。
【請求項６】
前記センサ信号の異常診断を行う工程における前記異常に関連するセンサを特定する工程は、前記観測ベクトルと前記基準ベクトルに基づき観測値と基準値の差を算出する工程と、前記算出された差が最大となる特徴を構成するセンサを第一の関連センサとする工程と、前記第一の関連センサにかかわる前記２次元の特徴分布密度に基づき、前記観測ベクトルに対応する値の密度が最小となる前記第一の関連センサのペアとなる特徴を構成するセンサを第二の関連センサとする工程からなる請求項３に記載の異常診断方法。
【請求項７】
前記センサ信号の異常診断を行う工程における前記異常に関連するセンサを特定する工程は、全ての前記２次元の特徴分布密度に基づき、前記観測ベクトルに対応する値の密度を算出する工程と、前記特徴毎に前記算出した密度を合計した値が最小となる特徴を構成するセンサを第一の関連センサとする工程と、前記第一の関連センサにかかわる前記２次元の特徴分布密度に基づき、前記観測ベクトルに対応する値の密度が最小となる前記第一の関連センサのペアとなる特徴を構成するセンサを第二の関連センサとする工程からなる請求項３に記載の異常診断方法。
【請求項８】
前記センサ信号の異常診断を行う工程における前記異常に関連するセンサを特定する工程は、前記学習データのうち前記第一の関連センサのセンサ信号が異常と判定された観測ベクトルと近いデータを用いて前記第一の関連センサを除く全組合せの２次元の部分的特徴分布密度を算出する工程と、全ての前記２次元の部分的特徴分布密度に基づき、前記観測ベクトルに対応する値の密度を算出する工程と、特徴毎に前記密度を合計した値が最小となる特徴を構成するセンサを第二の関連センサとする工程と、前記第二の関連センサにかかわる前記２次元の部分的特徴分布密度に基づき、前記観測ベクトルに対応する値の密度が最小となる前記第二の関連センサのペアとなる特徴を構成するセンサを第三の関連センサとする工程とをさらに含む請求項６または請求項７に記載の異常診断方法。
【請求項９】
前記センサ信号の異常診断を行う工程における前記異常に関連するセンサを特定する工程は、連続して異常検出された区間を特定する工程と、前記区間に含まれる全ての観測ベクトルについて、前記観測ベクトルと前記基準ベクトルに基づき観測値と基準値の差の絶対値の和を算出する工程と、前記算出された和が最大となる特徴を構成するセンサを第一の関連センサとする工程と、前記区間に含まれる全ての観測ベクトルについて、前記第一の関連センサにかかわる前記２次元の特徴分布密度に基づき、前記観測ベクトルに対応する値の密度の和および密度が非０となる回数を算出する工程と、前記算出された回数が最小かつそのうち前記算出された和が最小となる前記第一の関連センサのペアとなる特徴を構成するセンサを第二の関連センサとする工程からなる請求項３に記載の異常診断方法。
【請求項１０】
前記異常測度および前記しきい値および前記判定結果の時系列グラフと、前記特定された異常に関連するセンサから出力されたセンサ信号の時系列グラフと、前記特定された２次元の特徴の学習データの分布密度を表す画像に、前記異常と判定された観測ベクトルに対応する点を重ねてプロットした画像を表示する工程をさらに含むことを特徴とする請求項３に記載の異常診断方法。
【請求項１１】
設備または装置に装着された複数のセンサから出力される複数のセンサ信号をもとに学習データを作成して蓄積する工程と、前記設備または装置に装着された複数のセンサから新たに出力された複数のセンサ信号の異常診断を行う工程とを含む設備の異常を診断する方法であって、
前記学習データを作成して蓄積する工程は、設備または装置から出力されるイベント信号に基づいて稼動状態別のモード分割を行う工程と、前記複数のセンサ信号から特徴ベクトルを抽出する工程と、前記抽出した特徴ベクトルを学習データとして蓄積する工程と、前記学習データとして蓄積した特徴ベクトルのおのおのについて前記特徴ベクトルに応じて前記学習データの中から所定数の特徴ベクトルを選択する工程と、前記選択された所定数の特徴ベクトルを用いて学習用の基準ベクトルを作成する工程と、前記作成した学習用の基準ベクトルに基づき前記学習データとして蓄積した特徴ベクトルの異常測度を算出する工程と、前記算出した異常測度に基づき前記モード毎にしきい値を算出する工程と、前記学習データとして蓄積した特徴ベクトルから全組合せの２次元の特徴分布密度を前記モード毎に算出する工程を有し、
前記センサ信号の異常診断を行う工程は、前記設備または装置に装着された前記複数のセンサから出力される複数のセンサ信号から特徴ベクトルを抽出して観測ベクトルとする工程と、前記抽出した観測ベクトルに応じて前記蓄積した学習データの中から所定数の特徴ベクトルを選択する工程と、前記学習データの中から選択された所定数の特徴ベクトルを用いて異常診断用の基準ベクトルを作成する工程と、前記作成した異常診断用の基準ベクトルに基づいて前記観測ベクトルの異常測度を算出する工程と、前記算出した異常測度と前記分割したモードと前記分割したモード別に算出した前記しきい値に基づき前記観測ベクトルが異常か正常かを判定する工程と、前記判定により異常と判定された観測ベクトルに対応するセンサ信号が検出された時刻における前記観測ベクトルと前記モードと前記モード別に算出した前記２次元の特徴分布密度に基づき異常に関連するセンサを特定する工程とを有する
ことを特徴とする異常診断方法。
【請求項１２】
設備または装置に装着されたセンサから出力されるセンサ信号に基づいて前記設備または装置の異常を診断する装置であって、前記センサ信号を蓄積する生データ蓄積手段と、前記センサ信号の特徴値から特徴ベクトルを抽出する特徴ベクトル抽出手段と、予め指定された学習期間において前記特徴値から抽出した前記特徴ベクトルの集合と各時刻において前記特徴値から抽出した前記特徴ベクトルに基づいて前記各時刻における基準特徴ベクトルを算出する基準ベクトル算出手段と、前記各時刻の前記特徴ベクトルと前記基準特徴ベクトルの差に基づき異常測度を算出する異常測度算出手段と、前記予め指定された学習期間における前記異常測度に基づいてしきい値を算出するしきい値算出手段と、前記異常測度と前記しきい値との比較により異常を検出する異常検出手段と、前記予め指定された学習期間における前記特徴値の２次元の分布密度を算出する分布密度算出手段と、前記異常検出手段により検出された異常に対応するセンサ信号が検出された時刻における前記特徴値の２次元の分布密度に基づいて異常に関連するセンサを特定する関連センサ特定手段とを備えたことを特徴とする異常診断装置。
【請求項１３】
設備または装置に装着されたセンサから出力されるセンサ信号およびイベント信号に基づいて前記設備または装置の異常を診断する装置であって、前記センサ信号を蓄積する生データ蓄積手段と、前記イベント信号に基づいて稼動状態別のモード分割を行うモード分割手段と、前記センサ信号の特徴値から特徴ベクトルを抽出する特徴ベクトル抽出手段と、予め指定された学習期間において前記特徴値から抽出した前記特徴ベクトルの集合と各時刻において前記特徴値から抽出した前記特徴ベクトルに基づいて前記各時刻における基準特徴ベクトルを算出する基準ベクトル算出手段と、前記各時刻の前記特徴ベクトルと前記基準特徴ベクトルの差に基づき異常測度を算出する異常測度算出手段と、前記予め指定された学習期間における前記異常測度に基づいて前記モード分割手段で分割したモード毎にしきい値を算出するしきい値算出手段と、前記異常測度と前記算出したモード毎のしきい値との比較により異常を検出する異常検出手段と、前記予め指定された学習期間における前記特徴値の２次元の分布密度を前記モード毎に算出する分布密度算出手段と、前記異常検出手段により検出された異常に対応するセンサ信号が検出された時刻における前記モードに対応する前記特徴値の２次元の分布密度に基づいて異常に関連するセンサを特定する関連センサ特定手段とを備えたことを特徴とする異常診断装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、プラントや設備などの出力する多次元時系列データをもとに異常を早期に検知し、検知された異常に関連するセンサを特定する異常診断方法およびその装置に関する。
【背景技術】
【０００２】
電力会社では、ガスタービンの廃熱などを利用して地域暖房用温水を供給したり、工場
向けに高圧蒸気や低圧蒸気を供給したりしている。石油化学会社では、ガスタービンなどを電源設備として運転している。このようにガスタービンなどを用いた各種プラントや設備において、設備の不具合あるいはその兆候を検知する異常検知は、社会へのダメージを最小限に抑えるためにも極めて重要である。
【０００３】
ガスタービンや蒸気タービンのみならず、水力発電所での水車、原子力発電所の原子炉、風力発電所の風車、航空機や重機のエンジン、鉄道車両や軌道、エスカレータ、エレベータ、機器・部品レベルでも、搭載電池の劣化・寿命など、上記のような予防保全を必要とする設備は枚挙に暇がない。
【０００４】
このため、対象設備やプラントに複数のセンサを取り付け、センサ毎の監視基準に従って正常か異常かを判断することが行われている。米国特許第６，９５２，６６２号（特許文献１）には、観測データと過去の学習データとの類似度を独自の方法で計算し、類似度に応じて重み付けされたデータの線形結合により推定値を算出して、推定値と観測データのはずれ度合をもとに異常を検出する方法が開示されている。また、特開２０１１−７０６３５号公報（特許文献２）には、過去の正常データから作成されたモデルとの比較によって算出される異常測度に基づいて異常の有無を検知する異常検知方法において、正常モデルを局所部分空間法によって作成することが開示されている。
【０００５】
しかし、異常を検知するのみでは対策や調査などの次のアクションが決められないため、どのセンサが異常に関連しているのかを診断したいというニーズがある。
【０００６】
このようなニーズに対応する技術として、特開２０１２−１３８０４４号公報（特許文献３）には、複数の測定変数から算出されるＴ２統計量およびＱ統計量に基づいて異常を検出し、異常検出後にこれらの統計量に対する各測定変数の寄与量に基づいて異常要因候補を列挙するプロセス状態監視装置が開示されている。
【先行技術文献】
【特許文献】
【０００７】
【特許文献１】米国特許第６，９５２，６６２号明細書
【特許文献２】特開２０１１−７０６３５号公報
【特許文献３】特開２０１２−１３８０４４号公報
【発明の概要】
【発明が解決しようとする課題】
【０００８】
センサ毎の監視基準に従って正常か異常かを判断する方法では、監視基準を超えたセンサがすなわち異常に関連するセンサであるが、この方法では２種以上のセンサ信号の組合せとしての異常はそもそも検知できないという問題がある。
【０００９】
特許文献１や特許文献２に記載の方法によれば、学習データとして正常時のデータを与えることにより、学習にない観測データが観察されると、これらを異常として検出することができる。複数のセンサ信号をまとめて演算するため、センサ信号の組合せとしての異常も検出できる。しかし、異常に関連するセンサを正しく診断する方法については記載がない。
【００１０】
特許文献３に記載の方法によれば、複数の測定項目のＴ２統計量やＱ統計量を異常の指標とすることにより、２種以上のセンサ信号の組合せとしての異常も検知でき、統計量に対する各測定項目の寄与量に基づいて異常に関連するセンサを診断することが可能である。しかし、Ｔ２統計量やＱ統計量を指標とした異常検出は、運転と停止の切り替りや運転条件の切り替りが頻繁かつ不規則にある設備やプラントに適用すると、誤報が多くなってしまい、そのままでは適用困難である。特許文献２に記載の局所部分空間法は、学習データを網羅的に集めることにより、上記のようなケースでも誤報を抑制することが可能である。
【００１１】
しかし、特許文献１に記載の異常検知手法に、特許文献２に記載の寄与量に基づいて異常に関連するセンサを診断する方法を適用すると、誤る場合がある。この理由について簡単に説明する。局所部分空間法は、図２に示す通り、観測データ(q)に近い学習データ(x1，x2，x3)を用いて作られるアフィン部分空間への投影距離を異常の指標としている。つまり、観測データ(q)からアフィン部分空間への垂線の足(b)を基準として、bからqのベクトルの各要素の大きさが寄与量となる。あるセンサ信号に異常がある場合、これに基準(b)が追随して異常とは関係ない他のセンサの寄与量が大きくなることが起こりうる。
【００１２】
単純化した例を図３Ａ及び図３Ｂに示す。図３Ａにおいて、abcの３次元のセンサ信号において正常なデータが図示するように分布しており、観測データがq→q’→qと変化した場合、基準データはb→b’→bと変化する。時系列のグラフで表すと図３Ｂに示す通りとなるため、aとbの寄与量が大きくなるが、実際にはcが異常に関連していると考えられるため、誤って診断されることになる。このことは、局所部分空間法のみでなく、ＡＲＴ２のようなクラスタに基づく方法や特許文献１に記載の類似度ベースモデリングと呼ばれる手法など、観測データによって基準が変化する事例ベースの手法すべてにあてはまる。
【００１３】
そこで、本発明の目的は、上記課題を解決し、多次元時系列センサ信号に基づいて高感度に異常を検出し、検出した異常に関連するセンサを正しく診断する異常診断方法および装置を提供することにある。
【課題を解決するための手段】
【００１４】
上記課題を解決するために、本発明では、設備または装置に装着された複数のセンサから出力される複数のセンサ信号に基づいて設備または装置の異常を診断する方法及びその装置において、予め指定された学習期間における複数のセンサ信号の特徴値から多次元の特徴ベクトルを求め、診断時に複数のセンサ信号の特徴値から時刻毎に多次元特徴ベクトルを抽出し、求めた予め指定された学習期間の多次元の特徴ベクトルの集合と診断時に抽出した時刻毎の多次元の特徴ベクトルに基づいて時刻毎の基準特徴ベクトルを算出し、抽出した時刻毎の多次元特徴ベクトルと算出した時刻毎の基準特徴ベクトルとの差に基づき異常測度を算出し、算出した異常測度と所定のしきい値とを比較することにより異常を検出し、異常が検出された時刻における複数のセンサ信号の特徴値の２次元の分布密度に基づいて検出した異常に関連するセンサを複数のセンサの中から特定するようにした。
【００１５】
又、上記課題を解決するために、本発明では、設備または装置に装着された複数のセンサから出力される複数のセンサ信号をもとに学習データを作成して蓄積する工程と、設備または装置に装着された複数のセンサから新たに出力された複数のセンサ信号の異常診断を行う工程とを含む設備の異常を診断する方法及びその装置において、学習データを作成して蓄積する工程は、複数のセンサ信号から特徴ベクトルを抽出する工程と、抽出した特徴ベクトルを学習データとして蓄積する工程と、学習データとして蓄積した特徴ベクトルの各々について特徴ベクトルに応じて蓄積した学習データの中から所定数の特徴ベクトルを選択する工程と、選択された所定数の特徴ベクトルを用いて学習用の基準ベクトルを作成する工程と、作成した学習用の基準ベクトルに基づいて学習データとして蓄積した特徴ベクトルの各々について異常測度を算出する工程と、算出した異常測度に基づきしきい値を算出する工程と、特徴ベクトルから全組合せの２次元の特徴分布密度を算出する工程とを有し、センサ信号の異常診断を行う工程は、設備または装置に装着された複数のセンサから出力される複数のセンサ信号から特徴ベクトルを抽出して観測ベクトルとする工程と、抽出した観測ベクトルに応じて蓄積した学習データの中から所定数の特徴ベクトルを選択する工程と、学習データの中から選択された所定数の特徴ベクトルを用いて異常診断用の基準ベクトルを作成する工程と、作成した異常診断用の基準ベクトルに基づいて観測ベクトルの異常測度を算出する工程と、算出した異常測度と学習データを作成して蓄積する工程において算出したしきい値に基づき観測ベクトルが異常か正常かを判定する工程と、判定する工程において異常と判定された観測ベクトルに対応するセンサ信号が検出された時刻における観測ベクトルと２次元の特徴分布密度に基づき異常に関連するセンサを特定する工程とを有して構成した。
【００１６】
さらにまた、上記課題を解決するために、本発明では、設備または装置に装着された複数のセンサから出力される複数のセンサ信号をもとに学習データを作成して蓄積する工程と、前記設備または装置に装着された複数のセンサから新たに出力された複数のセンサ信号の異常診断を行う工程とを含む設備の異常を診断する方法及びその装置において、学習データを作成して蓄積する工程は、設備または装置から出力されるイベント信号に基づいて稼動状態別のモード分割を行う工程と、複数のセンサ信号から特徴ベクトルを抽出する工程と、抽出した特徴ベクトルを学習データとして蓄積する工程と、学習データとして蓄積した特徴ベクトルのおのおのについて特徴ベクトルに応じて学習データの中から所定数の特徴ベクトルを選択する工程と、選択された所定数の特徴ベクトルを用いて学習用の基準ベクトルを作成する工程と、作成した学習用の基準ベクトルに基づき学習データとして蓄積した特徴ベクトルの異常測度を算出する工程と、算出した異常測度に基づきモード毎にしきい値を算出する工程と、学習データとして蓄積した特徴ベクトルから全組合せの２次元の特徴分布密度をモード毎に算出する工程を有し、センサ信号の異常診断を行う工程は、設備または装置に装着された複数のセンサから出力される複数のセンサ信号から特徴ベクトルを抽出して観測ベクトルとする工程と、抽出した観測ベクトルに応じて蓄積した学習データの中から所定数の特徴ベクトルを選択する工程と、学習データの中から選択された所定数の特徴ベクトルを用いて異常診断用の基準ベクトルを作成する工程と、作成した異常診断用の基準ベクトルに基づいて観測ベクトルの異常測度を算出する工程と、算出した異常測度と分割したモードと分割したモード別に算出したしきい値に基づき観測ベクトルが異常か正常かを判定する工程と、判定により異常と判定された観測ベクトルに対応するセンサ信号が検出された時刻における観測ベクトルとモードとモード別に算出した２次元の特徴分布密度に基づき異常に関連するセンサを特定する工程とを有して構成した。
【発明の効果】
【００１７】
本発明によれば、各時刻のデータの異常検出を、学習データとその時刻のデータから算出される基準データからの距離に基づいて行うため、高感度な異常検出が可能である。異常が検出された時刻の関連センサ特定において、寄与量が大きいセンサをそのまま関連センサとするのではなく、センサ信号の２次元の分布密度に基づいて異常に関連するセンサを特定するため、正しく診断することができる。
【００１８】
以上の手法を適用したシステムにより、ガスタービンや蒸気タービンなどの設備のみならず、水力発電所での水車、原子力発電所の原子炉、風力発電所の風車、航空機や重機のエンジン、鉄道車両や軌道、エスカレータ、エレベータ、そして機器・部品レベルでは、搭載電池の劣化・寿命など、種々の設備・部品において、対象の異常を早期に検出するのみならず、関連センサに応じた調査あるいは対策を迅速に行うことが可能となる。
【図面の簡単な説明】
【００１９】
【図１】本発明の実施例１に関する設備状態監視システムの概略の構成を示すブロック図である。
【図２】局所部分空間法を説明する図である。
【図３Ａ】寄与量に基づくセンサ信号特定における誤りの例を説明する図で、正常データの分布を示す３次元のグラフである。
【図３Ｂ】寄与量に基づくセンサ信号特定における誤りの例を説明する図で、３つのセンサ信号の基準値に対する観測値の変化を示すセンサ信号のグラフである。
【図４】センサ信号の例を示す信号リストの表である。
【図５】センサ信号入力からしきい値設定までの処理の流れを示すフロー図である。
【図６】２次元の特徴の分布密度を算出する処理の流れを示すフロー図である。
【図７】センサ信号入力から異常判定までの処理の流れを示すフロー図である。
【図８】異常に関連するセンサを特定する処理の流れを示すフロー図の例である。
【図９】異常に関連するセンサを特定する処理の流れを示すフロー図の別の例である。
【図１０】異常に関連する３個のセンサを特定する処理の流れを示すフロー図である。
【図１１】異常が連続する区間毎に関連するセンサを特定する処理の流れを示すフロー図である。
【図１２Ａ】分布密度画像に異常区間の対応点をプロットした例を表す図である。
【図１２Ｂ】分布密度画像に異常区間の対応点をプロットした例を表す図である。
【図１２Ｃ】分布密度画像に異常区間の対応点をプロットした例を表す図である。
【図１２Ｄ】分布密度画像に異常区間の対応点をプロットした例を表す図である。
【図１３】レシピ設定のためのＧＵＩの１例を表す表示画面の正面図である。
【図１４Ａ】結果表示画面の例で複数の時系列データを表示した画面の正面図である。
【図１４Ｂ】結果表示画面の例で複数の時系列データを拡大表示した画面の正面図である。
【図１４Ｃ】結果表示画面の例で異常関連センサの分布密度と複数の時系列データを表示した画面の正面図である。
【図１５】レシピ設定のためのテスト結果の一覧を表示した画面の正面図である。
【図１６】結果表示対象指定のためのＧＵＩの１例を表す表示画面の正面図である。
【図１７】結果表示画面に含まれる期間表示ウィンドウの表示例である。
【図１８】本発明の実施例２に関する設備状態監視システムの概略の構成を示すブロック図である。
【図１９Ａ】本発明の実施例２に関するイベント信号の例を示す信号リストの表である。
【図１９Ｂ】本発明の実施例２に関するイベント信号に基づくモード分割処理の流れを示すフロー図である。
【図１９Ｃ】本発明の実施例２に関する設備の可動状態を分割して４種のモードの何れかに分類した状態を示すイベント信号の模式図である。
【発明を実施するための形態】
【００２０】
本発明は、設備または装置に装着された複数のセンサから出力される複数のセンサ信号に基づいて設備または装置の異常を診断する方法及びその装置において、センサ信号から時刻毎に多次元の特徴ベクトルを抽出し、予め指定された学習期間の特徴ベクトルの集合と各時刻の特徴ベクトルに基づいて各時刻の基準ベクトルを算出し、各時刻の特徴ベクトルと基準特徴ベクトルとの差に基づいて異常測度を算出し、異常測度と所定のしきい値との比較により異常を検出し、異常が検出された時刻における特徴値の２次元の分布密度に基づいて異常に関連するセンサを特定できるようにしたものである。
以下に、本発明の実施例を図を用いて説明する。
【実施例１】
【００２１】
以下、図面を用いて本発明の内容を詳細に説明する。
図１に、本発明の異常診断方法を実現するシステムの一構成例を示す。
本システムは、設備１０１から出力されるセンサ信号１０２を蓄積するセンサ信号蓄積部１０３、センサ信号１０２をもとに特徴ベクトルを抽出する特徴ベクトル抽出部１０４、予め指定された学習期間の特徴ベクトルの集合と各時刻の特徴ベクトルに基づいて各時刻の基準特徴ベクトルを算出する基準ベクトル算出部１０５、各時刻の特徴ベクトルと基準特徴ベクトルの差に基づき異常測度を算出する異常測度算出部１０６、予め指定された学習期間の異常測度に基づきしきい値を算出するしきい値算出部１０７、異常測度と算出されたしきい値との比較により異常を検出する異常検出部１０８、学習期間のセンサ信号の２次元の分布密度を算出する分布密度算出部１０９、異常が検出された時刻について、異常に関連するセンサを特定する関連センサ特定部１１０を備えて構成される。
【００２２】
本システムの動作には、蓄積されたデータを用いて学習データの生成、保存を行う「学習」と入力信号に基づき異常を検出し、関連センサの特定を行う「異常診断」の二つのフェーズがある。基本的に前者はオフラインの処理、後者はオンラインの処理である。ただし、後者をオフラインの処理とすることも可能である。以下の説明では、それらを学習時、異常診断時という言葉で区別する。
【００２３】
状態監視の対象とする設備１０１は、ガスタービンや蒸気タービンなどの設備やプラントである。設備１０１は、その状態を表すセンサ信号１０２を出力する。センサ信号１０２はセンサ信号蓄積部１０３に蓄積されている。センサ信号１０２をリスト化して表形式に表した例を図４に示す。センサ信号１０２は一定間隔毎に取得される多次元時系列信号であり、それをリスト化した表は、図４に示すように、日時の欄４０１と設備１０１に設けられた複数のセンサ値のデータの欄４０２からなる。センサの種類は、数百から数千といった数になる場合もあり、例えば、シリンダ、オイル、冷却水などの温度、オイルや冷却水の圧力、軸の回転速度、室温、運転時間などである。出力や状態を表すのみならず、何かをある値に制御するための制御信号の場合もある。
【００２４】
学習時の処理の流れを図５を用いて説明する。学習時はセンサ信号蓄積部１０３に蓄積されたデータのうち指定された期間のデータを用いて、特徴ベクトルを抽出する。ここで抽出されたデータを学習データと呼ぶ。交差検証によって学習データの異常測度を算出し、全学習データの異常測度に基づき異常判定のしきい値を算出する。また、分布密度算出部１０９において特徴値の２次元の分布密度を算出する。
【００２５】
まず、特徴ベクトル抽出部１０４、基準ベクトル算出部１０５、異常測度算出部１０６、しきい値算出部１０７における処理の流れを図５を用いて説明する。はじめに、特徴ベクトル抽出部１０４において、センサ信号蓄積部１０３から学習期間として指定された期間のセンサ信号１０２を入力し（Ｓ５０１）、センサ信号毎に正準化した後（Ｓ５０２）、特徴ベクトルの抽出を行う（Ｓ５０３）。次に、基準ベクトル算出部１０５において、抽出された特徴ベクトルに交差検証用のグループ番号を付与する（Ｓ５０４）。グループは、１日分を１グループとするなど、期間を指定して決めるか、予め指定されたグループ数に等分に分けて決める。
【００２６】
抽出された特徴ベクトルから、１個目の特徴ベクトルに注目し（Ｓ５０５）、注目ベクトルと異なるグループの特徴ベクトルの中から、注目ベクトルに近い順に予め指定された数の特徴ベクトルを選択し（Ｓ５０６）、選択した特徴ベクトルを用いて基準ベクトルを作成する（Ｓ５０７）。異常測度算出部１０６において、注目した特徴ベクトルの基準ベクトルまでの距離に基づいて異常測度を算出する（Ｓ５０８）。全ベクトルの異常測度算出が終了していれば（Ｓ５０９）、しきい値算出部１０７において、全ベクトルの異常測度に基づいてしきい値を設定する（Ｓ５１１）。ステップＳ５０９において全ての異常測度算出が終了していなければ、次の特徴ベクトルに注目し（Ｓ５１０）、ステップＳ５０６からＳ５０９の処理を繰り返す。
【００２７】
次に、各ステップについて詳細に説明する。  
ステップＳ５０２においては、各センサ信号の正準化を行う。例えば、指定された期間の平均と標準偏差を用いて、平均を０、分散を１となるように変換する。異常診断時に同じ変換ができるよう、各センサ信号の平均と標準偏差を記憶しておく。あるいは、各センサ信号の指定された期間の最大値と最小値を用いて最大が１、最小が０となるように変換する。あるいは、最大値と最小値の代わりに予め設定した上限値と下限値を用いてもよい。異常診断時に同じ変換ができるよう、各センサ信号の最大値と最小値または上限値と下限値を記憶しておく。センサ信号の正準化は、単位およびスケールの異なるセンサ信号を同時に扱うためのものである。
【００２８】
ステップＳ５０３においては、時刻毎に特徴ベクトル抽出を行う。センサ信号を正準化したものをそのまま並べることが考えられるが、ある時刻に対して±1，±2，…のウィンドウを設け，[ウィンドウ幅（3，5，…）×センサ数]の特徴ベクトルにより、データの時間変化を表す特徴を抽出することもできる。また、離散ウェーブレット変換(DWT: Discrete Wavelet Transform)を施して、周波数成分に分解してもよい。さらに、ステップＳ５０３において、特徴選択を行う。最低限の処理として、分散が非常に小さいセンサ信号および単調増加するセンサ信号を除く必要がある。
【００２９】
また、相関解析による無効信号を削除することも考えられる。これは、多次元時系列信号に対して相関解析を行い、相関値が１に近い複数の信号があるなど、極めて類似性が高い場合に、これらは冗長だとして、この複数の信号から重複する信号を削除し、重複しないものを残す方法である。このほか、ユーザが指定するようにしてもよい。また、長期変動が大きい特徴を除くことも考えられる。長期変動が大きい特徴を用いることは正常状態の状態数を多くすることにつながり、学習データの不足を引き起こすためである。例えば、１周期期間毎の平均と分散を算出し、それらのばらつきによって長期変動の大きさを推定できる。
【００３０】
基準ベクトル作成手法としては、局所部分空間法(LSC: Local Sub-space Classifier)や投影距離法(PDM: Projection Distance Method)が考えられる。
【００３１】
局所部分空間法は、注目ベクトルqのk−近傍ベクトルを用いてｋ−１次元のアフィン部分空間を作成する方法である。図２にｋ＝３の場合の例を示す。k-近傍ベクトルはすなわちステップＳ５０６において選択された特徴ベクトルである。指定する数は、ｋの値である。
【００３２】
図２に示すように、異常測度は図に示す投影距離で表されるため、注目ベクトルqに最も近いアフィン部分空間上の点bを求めればよい。評価データqとそのk−近傍ベクトルxi( i＝1，…，k )からbを算出するには、qをk個並べた行列Qとxiを並べた行列Xから  
【００３３】
【数１】
2014142697.tif 000003
により相関行列Cを求め、
【００３４】
【数２】
2014142697.tif 000004
によりbを計算する。ここまでが、ステップＳ５０７における基準ベクトル作成にあたる。
【００３５】
異常測度ｄはｑとｂの間の距離であるから次式で表される。  
【数３】
2014142697.tif 000005
【００３６】
なお、図２ではｋ＝３の場合を説明したが、特徴ベクトルの次元数より十分小さければいくつでもよい。ｋ＝１の場合は、最近傍法と等価の処理になる。
【００３７】
投影距離法は、選択された特徴ベクトルに対し独自の原点をもつ部分空間すなわちアフィン部分空間（分散最大の空間）を作成する方法である。ステップＳ５０６において指定する数はいくつでも良いが、大きすぎるとベクトルの選択、部分空間の算出ともに時間がかかってしまうので数十から数百とするとよい。
【００３８】
アフィン部分空間の算出方法について説明する。まず、選択された特徴ベクトルの平均μと共分散行列Σ を求め、次にΣの固有値問題を解いて値の大きい方から予め指定したr個の固有値に対応する固有ベクトルを並べた行列Ｕをアフィン部分空間の正規直交基底とする。ｒは特徴ベクトルの次元より小さくかつ選択データ数より小さい数とする。あるいはｒを固定した数とせず、固有値の大きい方から累積した寄与率が予め指定した割合を超えたときの値としてもよい。異常測度は、注目ベクトルのアフィン部分空間への投影距離とする。
【００３９】
この他、注目ベクトルqのk-近傍ベクトルの平均ベクトルまでの距離を異常測度とする局所平均距離法や、ガウシアンプロセスなどを用いてもよい。
【００４０】
ステップＳ５１１におけるしきい値を設定する方法について説明する。学習期間の全特徴ベクトルの異常測度を昇順にソートし、予め指定した１に近い比率に到達する値を求める。この値を基準としてオフセットを加える、定数倍するなどの処理によりしきい値を算出する。オフセット０、倍率を１とすれば、この値がしきい値となる。算出したしきい値は、図示していないが、学習データと対応付けて記録しておく。
【００４１】
上記に示したように、注目ベクトルの近傍ベクトルを用いて基準ベクトルを作成することにより、状態が複雑に変化する設備に対しても適切な基準により精度の高い異常測度を算出することが可能である。学習データの交差検証により算出された異常測度に基づいてしきい値を算出するため、誤報を抑制することができる。
【００４２】
次に、分布密度算出部１０９における処理の流れを図６を用いて説明する。  
始めに、学習期間の特徴ベクトルを入力する（Ｓ６０１）。最初の特徴に注目し（Ｓ６０２）、学習期間のデータの最大値（ＭＡＸ）と最小値（ＭＩＮ）を求める（Ｓ６０３）。次に、最小値から最大値を指定された数Ｎで分割する際の刻み幅Ｓを算出する（Ｓ６０４）。Ｓ＝（ＭＡＸ−ＭＩＮ）／Ｎ  で計算できる。次に、最小値と最大値から外側に範囲を広げて分布密度算出の処理範囲を算出する（Ｓ６０５）。広げる範囲は例えばＭＩＮをＭＩＮ−Ｓ×Ｍ、ＭＡＸをＭＡＸ＋Ｓ×Ｍに変更する。ここでＭは予め決められた１以上の整数とする。
【００４３】
次に学習期間の全データについて、特徴値（Ｆ）からビン番号（ＢＮＯ）を次式で算出する（Ｓ６０６）。  
ＢＮＯ＝ＩＮＴ（（Ｆ−ＭＩＮ）・Ｎ／（ＭＡＸ−ＭＩＮ））
ただしＩＮＴ（Ｘ）はＸの整数部を表す。
【００４４】
以上の処理について未処理の特徴が残っている場合は（Ｓ６０７）、次の特徴に注目し（Ｓ６０８）、ステップＳ６０３からＳ６０６の処理を行う。ステップＳ６０７において全特徴について処理していれば、ステップＳ６０９に進む。
【００４５】
ステップＳ６０９においては、特徴２個の最初の組合せに注目する。２個の特徴は同じものであってもよい。次に分布密度算出用の２次元配列を確保し、すべての要素に０をセットする（Ｓ６１０）。配列のサイズはＮ＋２Ｍである。学習期間の全データについて、２個の特徴のビン番号に対応する配列の要素に１を加算する（Ｓ６１１）。この処理により、特徴２個による２次元の頻度分布（ヒストグラム）が算出される。この頻度分布を画像に変換して保存する（Ｓ６１２）。変換方法については後述する。上記処理が特徴２個の全ての組合せについて実施済みであれば（Ｓ６１３）、処理を終了し（Ｓ６１５）、そうでなければ次の組合せに注目し（Ｓ６１４）、ステップＳ６１０からＳ６１２までの処理を行う。図示はしていないが、２次元配列のサイズおよびステップＳ６０５で算出した各特徴の最小値と最大値を記録しておく。
【００４６】
ステップＳ６１２における、画像変換方法の例を説明する。始めに配列要素の最大値すなわち最大頻度を求める。画像サイズは配列サイズと同じとし、各要素の値から対応する座標の画素値を例えば、[２５５×配列の要素値／最大頻度]とする。２５５は画素値を８ビットで表す場合の最大値であり、この値を用いれば、そのままビットマップ形式で保存できる。あるいは、画素値を[２５５×ＬＯＧ（配列の要素値＋１）／ＬＯＧ（最大頻度＋１）]とする。ただしＬＯＧ（Ｘ）はＸの対数を表す。このような変換式を用いれば、最大頻度が大きい場合も非ゼロの頻度に非ゼロの画素値を対応させることが可能になる。
【００４７】
上記処理により得られた画像は、２次元の特徴空間上で密度が高いところが高い画素値で表されているため、分布密度画像と呼ぶこととする。画像の作り方は、上記方法に限定されない。例えば単純な頻度分布ではなく、１個のデータにガウス分布や他の重みつきフィルタを割り当て、それを重畳するようにしてもよい。あるいは、上記方法で得られた画像に所定サイズの最大値フィルタをかけたり、平均フィルタ、その他の重みつきフィルタをかけたりしてもよい。また必ずしも画像形式で保存する必要はなく２次元配列をテキスト形式で保存してもよい。
【００４８】
異常診断時の処理の流れを図７ないし図９を用いて説明する。異常診断時はセンサ信号蓄積部１０３に蓄積されたデータのうち指定された期間のデータあるいは新たに観測されたデータの異常測度を算出し、正常か異常かの判定を行い、異常と判定された時刻について関連するセンサを特定する。
【００４９】
図７は特徴ベクトル抽出部１０４、基準ベクトル算出部１０５、異常測度算出部１０６、異常検出部１０８における異常診断時の処理の流れを説明する図である。始めに、特徴ベクトル抽出部１０４において、センサ信号蓄積部１０３または設備１０１からセンサ信号を入力し（Ｓ７０１）、センサ信号毎に正準化した後（Ｓ７０２）、特徴ベクトルの抽出を行う（Ｓ７０３）。センサ信号の正準化には、図５のステップＳ５０２に示す処理において、学習データの正準化に用いたものと同じパラメータを用いる。特徴ベクトルの抽出は、ステップＳ５０３の処理と同じ方法で行う。したがって、ステップＳ５０３において特徴選択を実行した場合は同じ特徴を選択する。ここで抽出された特徴ベクトルを、学習データと区別するために観測ベクトルと呼ぶこととする。
【００５０】
次に、基準ベクトル作成部１０５において、学習データの特徴ベクトルの中から、観測ベクトルに近い順に予め指定された数の特徴ベクトルを選択し（Ｓ７０４）、それらの特徴ベクトルを用いて基準ベクトルを作成する（Ｓ７０５）。異常測度算出部１０６において、観測ベクトルの基準ベクトルまでの距離に基づいて異常測度を算出する（Ｓ７０６）。異常検出部１０８において、学習時に算出したしきい値と異常測度を比較して、異常測度がしきい値より大きければ異常、そうでなければ正常と判定する（Ｓ７０７）。
【００５１】
図８は、関連センサ特定部１１０における処理の流れの一実施例を説明する図である。この処理は、ステップＳ７０７において、異常と判定されたデータについて実施する。始めに、観測ベクトルと基準ベクトルから、各特徴の観測値と基準値の差を算出し（Ｓ８０１）、これが最大となる特徴を探す（Ｓ８０２）。これを特徴Ａとしておく。次に観測ベクトルから、図６で示した分布密度算出処理で使用した２次元配列のサイズおよびステップＳ６０５で算出した各特徴の最小値と最大値を用いて、各特徴のビン番号を算出する（Ｓ８０３）。一方の特徴が特徴Ａである分布密度画像を読み込み、各特徴と特徴Ａのビン番号に対応する座標の画素値を読み（Ｓ８０４）、これが最小となる特徴を探す（Ｓ８０５）。これを特徴Ｂとし、特徴Ａと特徴Ｂの分布密度画像に観測ベクトルの対応点をプロットする（Ｓ８０６）。例えば、カラー画像に変換して、特徴Ａと特徴Ｂのビン番号に対応する座標に使用していない色をつける。Ｓ８０４における画素値が０であれば、あるいは予め決められた値より小さければ、特徴Ａと特徴Ｂを関連センサとする（Ｓ８０７）。
【００５２】
図９は、関連センサ特定部１１０における処理の流れの別の一実施例を説明する図である。図８の処理と同様、ステップＳ７０７において、異常と判定されたデータについて実施する。始めに、観測ベクトルから、図６で示した分布密度算出処理で使用した２次元配列のサイズおよびステップＳ６０５で算出した各特徴の最小値と最大値を用いて、各特徴のビン番号を算出する（Ｓ９０１）。
【００５３】
次に全ての分布密度画像を読み込み、それぞれの分布密度画像について２個のビン番号に対応する座標の画素値を読み（Ｓ９０２）、特徴毎にその値の合計siを算出する（Ｓ９０３）。つまり、特徴iと特徴jの分布密度画像から読み取った値をI(i，j)とすると、siは次式で計算される。  
【００５４】
【数４】
2014142697.tif 000006
【００５５】
このsiが最小となる特徴を探し（Ｓ９０４）、これを特徴Ｃとしておく。次に一方の特徴が特徴Ｃである分布密度画像について、各特徴と特徴Ｃのビン番号に対応する座標の画素値を読み（Ｓ９０５）、これが最小となる特徴を探す（Ｓ９０６）。これを特徴Ｄとし、ステップＳ８０６と同様の方法で特徴Ｃと特徴Ｄの分布密度画像に観測ベクトルの対応点をプロットする（Ｓ９０７）。Ｓ９０５における画素値が０であれば、あるいは予め決められた値より小さければ、特徴Ｃと特徴Ｄを関連センサとする（Ｓ９０８）。
【００５６】
また別の実施例では、ステップＳ９０１で各特徴のビン番号を算出し、ステップＳ９０２で全ての分布密度画像について観測ベクトルに対応する座標の画素値を読んだ後、これが最小となる特徴の組合せを探し、これらを特徴Ｅおよび特徴Ｆとする。ステップＳ８０６と同様の方法で特徴Ｅと特徴Ｆの分布密度画像に観測ベクトルの対応点をプロットし、特徴Ｅと特徴Ｆの分布密度画像から読んだ画素値が０であれば、あるいは予め決められた値より小さければ、特徴Ｅと特徴Ｆを関連センサとする。
【００５７】
以上、説明した関連センサ特定の方法はいずれかひとつを実施してもよいし、関連センサが特定できるまでシーケンシャルに実施してもよいし、全て実施して観測ベクトルをプロットした分布密度画像を表示して、ユーザに決めさせてもよい。
【００５８】
さらに、上記方法で関連センサが特定できない場合に詳細に解析する方法について、図１０を用いて説明する。始めに、学習期間の特徴ベクトルを入力する（Ｓ１００１）。次に学習期間の全データについて、図６で示した分布密度算出処理で使用した２次元配列のサイズおよびステップＳ６０５で算出した各特徴の最小値と最大値を用いて、各特徴のビン番号を算出する（Ｓ１００２）。同様に、観測ベクトルについて各特徴のビン番号を算出する（Ｓ１００３）。
【００５９】
次に、学習期間の特徴ベクトルから前述の特徴Ａのビン番号が観測ベクトルと同じデータのみ抽出し（Ｓ１００４）、これらを用いてステップＳ６１０からＳ６１２に従い、特徴Ａを除く組合せの部分的分布密度画像を作成する（Ｓ１００５）。全ての部分的分布密度画像について２個のビン番号に対応する座標の画素値を読み（Ｓ１００６）、特徴毎にその値の合計を算出する（Ｓ１００７）。これが最小となる特徴を探し（Ｓ１００８）、これを特徴Ｇとしておく。
【００６０】
一方の特徴が特徴Ｇである部分的分布密度画像について、各特徴と特徴Ｇのビン番号に対応する座標の画素値を読み（Ｓ１００９）、これが最小となる特徴を探し特徴Ｈとする（Ｓ１０１０）。特徴Ｇと特徴Ｈの部分的分布密度画像に観測ベクトルの対応点をプロットする（Ｓ１０１１）。Ｓ１００９における画素値が０であれば、あるいは予め決められた値より小さければ、特徴Ａと特徴Ｇと特徴Ｈを関連センサとする（Ｓ１０１２）。なお、これらの処理では特徴Ａのかわりに前述の特徴Ｃを用いてもよい。
【００６１】
上記方法は、異常検出した時刻毎に関連センサを特定する方法であるが、連続して異常検出された場合に適用すると冗長な処理になってしまう。そこで、異常が連続して検出されている間はバッファに情報をためておき、まとめて処理を行うことが考えられる。この場合、異常の継続時間や累積異常度でセンサ特定を行う対象を絞り込んでもよい。
【００６２】
連続して検出された一連の異常区間について関連するセンサを特定する方法を、図１１を用いて説明する。まず、異常検出部１０８においてステップＳ７０７で異常判定された結果をもとに、連続して異常検出された区間を特定する（Ｓ１１０１）。その間、観測ベクトルと基準ベクトルはバッファに蓄積しておく。なお「連続」というのは必ずしも厳密な意味ではなく、予め決められた長さの中断を含んでいてもよい。次に、上記区間に含まれる全ての観測ベクトルについて、特徴毎に観測値と基準値の差の絶対値の和を算出し（Ｓ１１０２）、これが最大となる特徴を探す（Ｓ１１０３）。これを特徴Ａとしておく。
【００６３】
次に、上記区間に含まれる全ての観測ベクトルについて、図６で示した分布密度算出処理で使用した２次元配列のサイズおよびステップＳ６０５で算出した各特徴の最小値と最大値を用いて、各特徴のビン番号を算出する（Ｓ１１０４）。一方の特徴が特徴Ａである分布密度画像を読み込み、上記区間に含まれる全ての観測ベクトルについて、各特徴と特徴Ａのビン番号に対応する座標の画素値を読む（Ｓ１１０５）。上記区間に含まれる全ての観測ベクトルについて特徴毎に非０の画素をカウントする（Ｓ１１０６）。特徴iについてカウントした値をcnt(i)とする。また、上記区間に含まれる全ての観測ベクトルについて特徴毎に画素値の合計を算出する（Ｓ１１０７）。特徴iについて算出した値をsum(i)とする。cnt(i)が最小となる特徴を探し、複数ある場合はそのうちsum(i)が最小となる特徴を探す（Ｓ１１０８）。これを特徴Ｂとし、特徴Ａと特徴Ｂの分布密度画像に異常区間に含まれる観測ベクトルの対応点をプロットする（Ｓ１１０９）。Ｓ１１０６における非０の画素数が区間に含まれる観測ベクトル数より小さければ、つまり０の画素数が１以上あれば、特徴Ａと特徴Ｂを関連センサとする（Ｓ１１１０）。
【００６４】
あるいは、cnt(i)が最小となる特徴が複数ある場合は、学習データの相関の高い特徴を探して特徴Ｂとしてもよい。相関の高さは例えば特徴Ａともう一方の特徴の相関値に基づいて決める。または、分布密度画像の非０の画素の数が少ないものを相関が高いと考える。
【００６５】
なお、上記に説明した図１１の処理フローは、図８に示す処理を複数の観測ベクトルに拡張した処理であるが、図９または図１０に示す処理を複数の観測ベクトルに拡張して用いてもよいし、それらを組み合わせて用いてもよい。
【００６６】
ステップＳ１１０９にて作成される分布密度画像の例を図１２Ａ〜Ｃに示す。横軸と縦軸はそれぞれの特徴の相対的な値を表し、色の濃い部分が学習データの密度が高い部分である。星印は異常区間に含まれる観測データを表す。実際には色をつけて表示すればよい。
【００６７】
図１２Ａは、特徴ａ、ｂの分布密度を表す。学習データの分布からは、特徴ａ、ｂとも低い状態と特徴ａ、ｂとも高い状態の定常状態があり状態遷移時は特徴ａに遅れて特徴ｂが上昇していることがわかる。時間ずれの関係は安定している。特徴ａ上昇のタイミングに対して特徴ｂの上昇が早く起こっていることが異常として検出されている。あるいは下降タイミングのずれを検出しているかもしれないが、どちらであるかはセンサ信号の時系列グラフから確認できる。
【００６８】
図１２Ｂは、特徴ｃ、ｄの分布密度を表す。特徴ｃ、ｄとも低い状態、特徴ｃが高く特徴ｄが低い状態、特徴ｃ、ｄとも高い状態とそれらの中間の状態はあるが特徴ｄが高く特徴ｃが低い状態にはならない。そのような状態が異常として検出されている。
【００６９】
図１２Ｃは、横軸も縦軸も特徴eの値を表したものである。実質的には２次元ではないが、同じ形式の画像として表すことにより、表示やセンサ特定の処理を簡便にすることができる。学習データは直線上に分布する。特徴eの値が学習データのどの値よりも高いことが異常として検出されている。
【００７０】
上記の通り、予め学習データを用いて全ての組合せの２次元特徴分布密度を算出しておき、異常が検出されたときの観測ベクトルの対応点の分布密度に基づいて異常に関連するセンサを特定するため、図３Ａ及び図３Ｂに示すような状況である場合でも正しいセンサを特定することが可能である。また、２次元特徴分布密度を画像で表し、これをユーザに提示することにより、それらのセンサが異常に関連することが理解されやすくなる。また、学習データの分布の傾向を知ることができる。例えば、図１２Ｄに示すように、学習データの密度の低い部分が広い範囲に存在する場合に、星印１２０１で示す異常を検出したとしてもそれが本当に異常かどうかは疑わなくてはならない。その場合、学習期間を増やして密度が高くなるようにするか、この組合せでは一緒に処理しないことにする、具体的には特徴ｇを処理対象からはずすようにするとよい。
【００７１】
以上の方法を実現するシステムのＧＵＩの実施例を説明する。  
学習期間および処理パラメータ設定のためのＧＵＩの例を、図１３に示す。以下の説明ではこの設定のことを単にレシピ設定と呼ぶことにする。また、過去のセンサ信号１０２は設備ＩＤおよび時刻と対応付けられてデータベースに保存されているものとする。レシピ設定画面１３０１では、対象装置、学習期間、使用センサ、基準算出パラメータ、しきい値設定パラメータを入力する。設備ＩＤ入力ウィンドウ１３０２には、対象とする設備のＩＤを入力する。設備リスト表示ボタン１３０３押下により図示はしていないがデータベースに保存されているデータの装置ＩＤのリストが表示されるので、リストから選択入力する。
【００７２】
学習期間入力ウィンドウ１３０４には、学習データを抽出したい期間の開始日と終了日を入力する。センサ選択ウィンドウ１３０５には、使用するセンサを入力する。リスト表示ボタン１３０６のクリックによりセンサリスト１３０７が表示されるので、リストから選択入力する。リストから複数選択することも可能である。基準算出パラメータ入力ウィンドウ１３０８には、正常モデル作成において指定するパラメータを入力する。図は正常モデルとして局所部分空間を採用した場合の例であり、モデル作成に使う近傍ベクトル数と正則化パラメータを入力する。正則化パラメータは、数２において相関行列Cの逆行列が求められないことを防ぐため、対角成分に加算する小さい数である。
【００７３】
しきい値設定パラメータ入力ウィンドウ１３０９には、図５に示す処理における交差検証のグループをどのように決めるかをラジオボタンで選択する。１日を１グループとするか、指定された数のグループに分割するかを選び、後者の場合はグループ数を入力する。また、ステップＳ５１１のしきい値設定処理において累積ヒストグラムに適用する比率を値入力する。レシピ名入力ウィンドウ１３１０には、入力された情報に対応付けるユニークな名前を入力する。全ての情報を入力したらテスト期間入力ウィンドウ１３１１にテスト対象期間を入力し、テストボタン１３１２の押下により、レシピのテストを行う。
【００７４】
この操作により、同じレシピ名で実行したテストの通し番号が採番される。次に、特徴ベクトル抽出部１０４において、指定した学習期間のセンサ信号１０２から特徴ベクトルの抽出を行う。図５で説明したステップＳ５０２の正準化においては、指定した学習期間の全データを用いて平均と標準偏差を求める。この平均と標準偏差の値もセンサ毎にレシピ名およびテスト番号に対応付けて保存しておく。
【００７５】
図５で説明したステップＳ５０７に対応して、交差検証により基準ベクトル算出部１０５において各時刻の基準ベクトルを作成し、ステップＳ５０８に対応して、異常測度算出部１０６において異常測度を算出する。更に、ステップＳ５１１に対応して、しきい値算出部１０７において異常判定しきい値を算出し、レシピ名およびテスト番号と対応付けて保存する。分布密度算出部１０９において、図６に示した処理により２次元の特徴の分布密度画像を作成し、保存する。併せて、ステップＳ６０５で算出した処理範囲を全ての特徴について保存しておく。その他、装置ＩＤ情報、使用センサ情報、学習期間、特徴ベクトル抽出に用いるパラメータ、基準作成パラメータをレシピ名およびテスト番号と対応付けて保存する。
【００７６】
次に、指定したテスト期間のセンサ信号１０２を入力として、図７に示した異常検知の処理を行う。検出された異常について、図８乃至図１１の何れかに示した手順で関連センサ特定の処理を行う。検出された異常は区間毎に通し番号が付加される。区間毎に、通し番号、区間開始時刻、区間終了時刻、複数の関連センサ名と観測ベクトルの対応点がプロットされた分布密度画像名が記録される。
【００７７】
テストの結果をユーザに示すためのＧＵＩの例を図１４Ａおよび図１４Ｂおよび図１４Ｃに示す。各画面の上部に表示されたタブを選択することにより、結果全体表示画面１４０１と結果拡大表示画面１４０２と異常診断結果表示画面１４０３を切り換えることができる。
【００７８】
図１４Ａには、結果全体表示画面１４０１を示す。結果全体表示画面１４０１には、指定された全期間の異常測度、しきい値、判定結果とセンサ信号の時系列グラフを表示する。期間表示ウィンドウ１４０４には、指定された学習期間およびテスト期間が表示される。異常測度表示ウィンドウ１４０５には、指定された学習期間およびテスト期間の異常測度としきい値と判定結果が表示される。センサ信号表示ウィンドウ１４０６には、指定された期間の指定されたセンサの出力値が表示される。センサの指定は、センサ名選択ウィンドウ１４０７への入力によって行う。ただし、ユーザが指定する前は、先頭のセンサが選択されている。カーソル１４０８は、拡大表示の時の起点を表し、マウス操作により移動できる。
【００７９】
表示日数指定ウィンドウ１４０９には、この画面では使用しないが、結果拡大表示画面１４０２での、拡大表示の起点から終点までの日数が表示される。この画面で入力することもできる。日付表示ウィンドウ１４１０には、カーソル位置の日付が、表示される。終了ボタン１４１１押下により結果全体表示画面１４０１、結果拡大表示画面１４０２、異常診断結果表示画面１４０３とも消去し終了する。
【００８０】
図１４Ｂには、結果拡大表示画面１４０２を示す。結果拡大表示画面１４０２には、結果全体表示画面１４０１において、カーソル１４０８で示された日付を起点として、指定された日数の異常測度、しきい値、判定結果とセンサ信号の時系列グラフを表示する。期間表示ウィンドウ１４０４には、結果全体表示画面１４０１と同じ情報が表示される。異常測度表示ウィンドウ１４０５およびセンサ信号表示ウィンドウ１４０６には、結果全体表示画面１４０１と同様の情報が、拡大表示される。表示日数指定ウィンドウ１４０９で、拡大表示の起点から終点までの日数を指定する。
【００８１】
日付表示ウィンドウ１４１０には、拡大表示の起点の日付が表示されている。スクロールバー１４１２で表示の起点を変更することも可能であり、この変更はカーソル１４０８の位置と日付表示ウィンドウ１４１０の表示に反映される。スクロールバー表示領域１４１３の全体の長さは結果全体表示画面１４０１に表示されている全期間に相当する。また、スクロールバー１４１２の長さは表示日数指定ウィンドウ１４０９で指定された日数に相当し、スクロールバー１４１２の左端部が拡大表示の起点に対応する。異常が検出された時刻には、異常測度グラフの対応する位置に異常区間番号１４１４がバルーン表示される。終了ボタン１４１１押下により終了する。
【００８２】
図１４Ｃには、異常診断結果表示画面１４０３の例を示す。異常に関連するセンサの２次元の特徴分布密度と時系列グラフを併せて表示する。異常診断結果表示画面１４０３に切り替ったときには、結果拡大表示画面１４０２に表示されている異常のうち最も番号の若い異常を表示対象とする。日付表示ウィンドウ１４１５には、表示対象の異常の日付が表示される。異常区間番号表示ウィンドウ１４１６には、表示対象の異常の区間番号が表示される。ウィンドウ１４１６では数値入力により、異常区間番号を指定することができる。時刻表示ウィンドウ１４１７には、表示対象の異常の区間開始時刻と区間終了時刻が表示される。
【００８３】
分布密度表示ウィンドウ１４１８には、例えば図１１のステップＳ１１０９において保存された画像を表示する。この画像は、ステップＳ１１１０において特定された２個のセンサの２次元頻度分布を表したものであり、横軸は”SensorX”の規格化された値、縦軸は”SensorY” の規格化された値を表す。１４１９は、学習期間の分布を表しており、濃く表示されているほど高頻度であることを意味する。１４２０は、表示対象の異常区間の観測ベクトルに対応するデータである。異常表示指示ウィンドウ１４２２で、チェックボックスのチェックをはずすことにより、表示しないようにすることもでき、その特徴値の学習期間のデータの分布密度を確認できる。処理範囲表示ウィンドウ１４２１には２個のセンサのステップＳ６０５において算出された処理範囲が表示される。これらの値は、それぞれ分布密度画像の左端に対応する”SensorX”の値、右端に対応する”SensorX”の値、下端に対応する”SensorY”の値、上端に対応する”SensorY”の値である。
【００８４】
異常測度表示ウィンドウ１４２３には、表示対象の異常区間を含む期間の異常測度としきい値と判定結果が表示される。表示期間は信号の変化が観測可能なよう十分に拡大される長さを予め決めておく。第一関連センサ信号表示ウィンドウ１４２４には、ステップＳ１１１０において特定された２個のセンサのうち一方の、異常測度表示ウィンドウ１４２３と同じ期間の出力値がグラフ表示される。第二関連センサ信号表示ウィンドウ１４２５には、もう一方のセンサ信号出力値がグラフ表示される。基準表示指示ウィンドウ１４２６のチェックボックスにチェックを入れることにより、図示していないが第一関連センサ信号表示ウィンドウ１４２４および第二関連センサ信号表示ウィンドウ１４２５に基準値を異なる色で重ねて表示することができる。
【００８５】
ボタン１４２７及び１４２８は、表示する異常診断結果を切り替えるためのボタンで、ボタン１４２７「次」をクリックすると、現在表示されている異常のデータよりも異常診断結果表示画面１４０３において次に若い番号の異常のデータが表示される。また、ぼたん１４２８「前」をクリックすると、現在表示されている異常のデータよりも異常診断結果表示画面１４０３において１つ若い番号の異常のデータが表示される。
【００８６】
又、図示はしていないが、分布密度表示ウィンドウ１４１８には、ステップＳ１１０９において保存された表数の画像を並べて表示するようにしてもよい。
【００８７】
本図では２個のセンサの時系列グラフを別々に表示したが、１個のウィンドウに異なる色で重ねて表示してもよい。終了ボタン１４１１押下により終了する。
【００８８】
図１４Ａ〜Ｃに示すいずれかの画面で、終了ボタン１４１１押下により異常検出結果および診断結果の確認が終了したら、図１３に示すレシピ設定画面１３０１の表示に戻る。テスト番号表示ウィンドウ１３１３には、上記のテストで採番された番号が表示されている。確認した内容に問題があれば、学習期間や選択センサ、パラメータなどを変更し、テストボタン１３１２の押下により、再度テストを行う。あるいは、一度行ったテストの結果を再度確認することもできる。テスト番号表示ウィンドウ１３１３からテスト番号を選択入力し、表示ボタン１３１４を押下する。この操作により、レシピ名とテスト番号に対応付けて保存された情報をロードし結果全体表示画面１４０１を表示する。タブの切り替えにより結果拡大表示画面１４０２または異常診断結果表示画面１４０３を表示させることもできる。確認が済んだら終了ボタン１４１１押下により、レシピ設定画面１３０１の表示に戻る。
【００８９】
登録ボタン１３１５の押下により、上記レシピ名とテスト番号表示ウィンドウ１３１３に表示中のテスト番号に対応付けて保存されている情報をレシピ名と対応付けて登録し、終了する。キャンセルボタン１３１６が押下された場合は、何も保存しないで終了する。
【００９０】
また、テスト結果一覧ボタン１３１７が押下された場合は、図１５に示す、テスト結果一覧表示画面１５０１を表示する。テスト結果リスト１５０２には、全てのテストの学習期間、テスト期間、選択センサ番号、基準作成パラメータ、しきい値設定パラメータなどのレシピ情報と、しきい値、異常区間数などのテスト結果情報を表示する。リストの左端に選択チェックボタンがあり、いずれか一つのみ選択することができる。詳細表示ボタン１５０３押下により、レシピ名とテスト番号に対応付けて保存された情報をロードし、結果全体表示画面１４０１を表示する。タブの切り替えにより結果拡大表示画面１４０２または異常診断結果表示画面１４０３を表示させることもできる。確認が済んだら終了ボタン１４１１押下により、テスト結果一覧表示画面１５０１の表示に戻る。登録ボタン１５０４の押下により、選択中のテスト番号に対応付けて保存されている情報をレシピ名と対応付けて登録し、テスト結果一覧表示画面１５０１の表示およびレシピ設定画面１３０１の表示を終了する。戻るボタン１５０５が押下された場合は、レシピの登録は行わずにレシピ設定画面１３０１の表示に戻る。
【００９１】
登録されたレシピは、活性か不活性かのラベルをつけて管理され、新しく観測されたデータに対しては、装置ＩＤが一致する活性なレシピの情報を用いて図７ないし図１１を用いて説明した何れかの特徴ベクトル抽出から異常検出、関連センサ特定までの処理を行い、結果をレシピ名と対応付けて保存しておく。
【００９２】
以上の異常検診断処理の結果をユーザに示すためのＧＵＩの例を、図１６に示す。  
図１６は、表示対象を指定するＧＵＩの例である。表示対象指定画面１６０１から表示対象の設備、レシピおよび期間を指定する。初めに、装置ＩＤ選択ウィンドウ１６０２により装置ＩＤを選択する。次に、レシピ名選択ウィンドウ１６０３により、装置ＩＤを対象としたレシピのリストから表示対象のレシピを選択する。データ記録期間表示部１６０４には、入力されたレシピを用いて処理され、記録が残されている期間の開始日と終了日が表示される。結果表示期間指定ウィンドウ１６０５には、結果を表示したい期間の開始日と終了日を入力する。表示センサ指定ウィンドウ１６０６には、表示したいセンサの名を入力する。表示ボタン１６０７押下により図１４Ａに示す結果全体表示画面１４０１を表示する。終了ボタン１６０８押下により終了する。
【００９３】
結果表示にかかわるＧＵＩの画面および操作は、図１４Ａ〜Ｃに示すテスト結果表示にかかわるＧＵＩとほぼ同じであるため、異なる部分のみ説明する。結果全体表示画面１４０１、結果拡大表示画面１４０２における期間表示ウィンドウ１４０４に対応する表示ウィンドウとしては、図１７に示すように結果表示期間指定ウィンドウ１６０５で指定された表示期間が表示される。結果全体表示画面１４０１に対応する画面において、異常測度表示ウィンドウ１４０５に対応するウィンドウには、指定された表示期間の異常測度としきい値と判定結果が表示される。センサ信号表示ウィンドウ１４０６に対応するウィンドウには、指定された期間の、表示センサ指定ウィンドウ１６０６により指定されたセンサの出力値が表示される。表示対象センサは、センサ名選択ウィンドウ１４０８に対応するウィンドウへの入力によって変更することも可能である。
【００９４】
上記実施例は学習データ設定をオフライン、異常診断処理をリアルタイム、結果表示をオフラインでそれぞれ処理するものであるが、結果表示もリアルタイムに行うことが可能である。その場合、表示期間の長さ、表示対象とするレシピ、表示対象とする情報を予め定めておき、一定時間毎に最新の情報を表示するよう構成すればよい。
【００９５】
逆に、任意の期間を設定し、レシピを選択して、オフラインで異常診断処理を行う機能を付加したものも本発明の範囲に含まれる。
【００９６】
本実施例によれば、状態が複雑に偏けする設備に対しても精度の高い異常測度を算出することができる。又、誤報の発生を抑制することができる。
【００９７】
更に、本実施例によれば、正常なセンサの信号から得られるデータが異常なセンサの信号のデータの影響を受けるような場合であっても、検出した異常に関連するセンサを正しく診断することが可能になった。
【００９８】
また、検査の結果得られる２次元特徴分布密度を画像で表示するので、異常に関連するセンサの見分けが容易になる。さらに、学習データの分布の傾向を画面上で確認することができるようになった。
【実施例２】
【００９９】
以上、設備から出力されるセンサ信号に基づき異常診断する方法の実施例を説明したが、別の実施例として、さらに、設備から出力されるイベント信号も利用して異常診断する方法を説明する。
【０１００】
図１８に、本発明の異常診断方法を実現するシステムの第二の実施例を示す。本実施例における構成は、実施例１で説明した図１に示した構成に、モード分割部１８０２を加えたものとなっている。実施例１と同じ構成については同じ番号を付してある。モード分割部１８０２は、設備からイベント信号１８０１を入力し、これに基づき設備の稼動状態を表すモードに分割する。モード分割の結果は、しきい値算出部１０７´に入力され、図５で説明した処理フローが分割されたモードごとに実行されてステップＳ５１１のしきい値設定がモード別に実施される。
【０１０１】
また、図７で説明した処理フローも分割されたモードごとに実行されて、異常検出部１０８´におけるステップＳ７０７の異常判定を、対応するモードのしきい値を用いて行う。また、分布密度算出部１０９´でもモード分割の結果を入力し、分布密度画像作成をモード毎に行い、関連センサ特定部１１０´におけるセンサ特定の処理をモード毎に行う。
【０１０２】
上記以外は、全て前述の実施例１で説明した方法と同様であるため、実施例１と異なるイベント信号に基づくモード分割方法の実施例を図１９Ａ〜Ｃを用いて説明する。イベント信号の例を図１９Ａに示す。図１９Ａに示したリスト１９２０は、不定期に出力される設備の操作・故障・警告を表す信号であり、時刻と操作・故障・警告を表す文字列またはコードからなる。
【０１０３】
図１９Ｂに示すように、このイベント信号をモード分割部１８０２に入力し（Ｓ１９０１）、所定の文字列またはコードの検索により起動シーケンスと停止シーケンスの切り出しを行う（Ｓ１９０２）。その結果をもとに、停止シーケンスの終了時刻から起動シーケンスの開始時刻までの「定常OFF」モード１９１１、起動シーケンス中の「起動」モード１９１２、起動シーケンスの終了時刻から停止シーケンスの開始時刻までの「定常ON」モード１９１３、停止シーケンス中の「停止」モード１９１４の４つの稼動状態に分割する（Ｓ１９０３）。図１９Ｃに例を示す。
【０１０４】
シーケンス切り出しのためには、予めシーケンスの開始イベントおよび終了イベントを指定しておき、イベント信号の先頭から最後まで以下の要領でスキャンしながら切り出していく。
(1)シーケンスの途中でない場合は、開始イベントを探索する。見つかったらシーケンスの開始とする。
(2)シーケンスの途中の場合は、終了イベントを探索する。見つかったらシーケンスの終了とする。ここで終了イベントとは、指定の終了イベントのほか、故障、警告、指定の開始イベントとする。
【０１０５】
以上のように、イベント信号を利用することにより、多様な稼動状態を正確に分けることができ、モード別にしきい値を設定することにより、「起動」モード１９１２および「停止」モード１９１４の過渡期において学習データ不足により感度を落とす必要がある場合でも、「定常OFF」モード１９１１および「定常ON」モード１９１３では高感度な異常検知が可能になる。
【０１０６】
即ち、本実施例によれば、設備の多様な稼働状態に応じてしきい値を設定することができ、高感度な異常検出を可能にし、更に検出した異常に関連するセンサを正しく診断することが可能になった。また、学習データの分布密度をモード毎に算出することにより、異なる状態のデータが含まれなくなるため、診断が容易になり、分布密度画像の理解も容易になる。
【符号の説明】
【０１０７】
１０１・・・設備  １０２・・・センサ信号  １０３・・・センサ信号蓄積部  １０４・・・特徴ベクトル抽出部  １０５・・・基準ベクトル作成部  １０６・・・異常測度算出部  １０７・・・しきい値算出部  １０８・・・異常検出部  １０９・・・分布密度算出部  １１０・・・関連センサ特定部  １３０１・・・レシピ設定画面    １４０１・・・結果全体表示画面  １４０２・・・結果拡大表示画面  １４０３・・・異常診断結果表示画面    １５０１・・・テスト結果一覧表示画面    １６０１・・・表示対象指定画面    １８０１・・・イベント信号  １８０２・・・モード分割部。
【図１】
2014142697.tif 000007
【図２】
2014142697.tif 000008
【図３Ａ】
2014142697.tif 000009
【図３Ｂ】
2014142697.tif 000010
【図４】
2014142697.tif 000011
【図５】
2014142697.tif 000012
【図６】
2014142697.tif 000013
【図７】
2014142697.tif 000014
【図８】
2014142697.tif 000015
【図９】
2014142697.tif 000016
【図１０】
2014142697.tif 000017
【図１１】
2014142697.tif 000018
【図１２Ａ】
2014142697.tif 000019
【図１２Ｂ】
2014142697.tif 000020
【図１２Ｃ】
2014142697.tif 000021
【図１２Ｄ】
2014142697.tif 000022
【図１３】
2014142697.tif 000023
【図１４Ａ】
2014142697.tif 000024
【図１４Ｂ】
2014142697.tif 000025
【図１４Ｃ】
2014142697.tif 000026
【図１５】
2014142697.tif 000027
【図１６】
2014142697.tif 000028
【図１７】
2014142697.tif 000029
【図１８】
2014142697.tif 000030
【図１９Ａ】
2014142697.tif 000031
【図１９Ｂ】
2014142697.tif 000032
【図１９Ｃ】
2014142697.tif 000033