特許 6658250 異常診断方法、異常診断装置、及び異常診断プログラム

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】6658250

(24)【登録日】2020年2月10日

(45)【発行日】2020年3月4日

(54)【発明の名称】異常診断方法、異常診断装置、及び異常診断プログラム

(51)【国際特許分類】

   G01M 99/00 20110101AFI20200220BHJP

   G05B 23/02 20060101ALI20200220BHJP

【ＦＩ】

   G01M99/00 Z

   G05B23/02 302S

【請求項の数】5

【全頁数】19

(21)【出願番号】特願2016-84511(P2016-84511)

(22)【出願日】2016年4月20日

(65)【公開番号】特開2017-194341(P2017-194341A)

(43)【公開日】2017年10月26日

【審査請求日】2019年2月21日

(73)【特許権者】

【識別番号】000000099

【氏名又は名称】株式会社ＩＨＩ

(74)【代理人】

【識別番号】100088155

【弁理士】

【氏名又は名称】長谷川 芳樹

(74)【代理人】

【識別番号】100113435

【弁理士】

【氏名又は名称】黒木 義樹

(74)【代理人】

【識別番号】100176245

【弁理士】

【氏名又は名称】安田 亮輔

(74)【代理人】

【識別番号】100171583

【弁理士】

【氏名又は名称】梅景 篤

(72)【発明者】

【氏名】鈴木 由宇

(72)【発明者】

【氏名】河野 幸弘

【審査官】 川瀬 正巳

(56)【参考文献】

【文献】 特開２００３−３１５２１１（ＪＰ，Ａ）

【文献】 国際公開第２０１３／１８７２９５（ＷＯ，Ａ１）

【文献】 特開２０１６−０２３４８９（ＪＰ，Ａ）

【文献】 米国特許出願公開第２０１５／０１１２８９１（ＵＳ，Ａ１）

(58)【調査した分野】（Int.Cl.，ＤＢ名）

Ｇ０１Ｍ ９９／００

Ｇ０５Ｂ ２３／０２

(57)【特許請求の範囲】

【請求項1】

診断対象システムに関する複数の計測項目のデータに基づいて前記診断対象システムの異常を診断する異常診断方法であって、
前記診断対象システムが正常に稼働している時間帯における複数の第１計測時点において取得された前記計測項目のデータである正常稼働データに基づいて、目的変数が前記計測項目のうちの一の診断項目であり、説明変数が他の前記計測項目である回帰式の回帰係数の確率分布を表す値を算出する第１ステップと、
前記回帰式、前記第１ステップで算出した前記回帰係数の確率分布を表す値、及び診断の対象となる時間帯における第２計測時点において取得された前記計測項目のデータである診断対象データに基づいて、前記第２計測時点における前記診断項目の予測値の確率分布を表す値を算出する第２ステップと、
前記診断対象データ及び前記第２ステップで算出した前記診断項目の予測値の確率分布を表す値に基づいて前記第２計測時点における前記診断項目の異常度である第１異常度を算出し、前記第１異常度に基づいて前記診断項目の異常判定を行う第３ステップと、を含み、
前記第３ステップでは、前記第２計測時点について、前記診断項目の計測値と前記診断項目の予測値の確率分布における平均値との差分を当該確率分布の標準偏差で除した値を前記第１異常度として算出する、異常診断方法。

【請求項2】

前記第１ステップでは、前記回帰係数の確率分布が正規分布に従うと仮定した場合の前記回帰係数の確率分布を表す値を算出し、
前記第２ステップでは、前記回帰係数の確率分布が正規分布に従うと仮定した場合の前記第２計測時点における前記診断項目の予測値の確率分布を表す値を算出する、請求項１記載の異常診断方法。

【請求項3】

前記正常稼働データ、前記回帰式、及び前記第１ステップで算出した前記回帰係数の確率分布を表す値に基づいて、各前記第１計測時点における前記診断項目の予測値の確率分布を表す値を算出する第４ステップと、
前記正常稼働データ及び前記第４ステップで算出した前記診断項目の予測値の確率分布を表す値に基づいて、各前記第１計測時点における前記診断項目の異常度である第２異常度を算出する第５ステップと、を更に含み、
前記第３ステップでは、前記第１異常度及び前記第２異常度に基づいて前記診断項目の異常判定を行う、請求項１又は２記載の異常診断方法。

【請求項4】

診断対象システムに関する複数の計測項目のデータに基づいて前記診断対象システムの異常を診断する異常診断装置であって、
前記診断対象システムが正常に稼働している時間帯における複数の第１計測時点において取得された前記計測項目のデータである正常稼働データに基づいて、目的変数が前記計測項目のうちの一の診断項目であり、説明変数が他の前記計測項目である回帰式の回帰係数の確率分布を表す値を算出する第１算出部と、
前記回帰式、前記第１算出部により算出された前記回帰係数の確率分布を表す値、及び診断の対象となる時間帯における第２計測時点において取得された前記計測項目のデータである診断対象データに基づいて、前記第２計測時点における前記診断項目の予測値の確率分布を表す値を算出する第２算出部と、
前記診断対象データ及び前記第２算出部により算出された前記診断項目の予測値の確率分布を表す値に基づいて前記第２計測時点における前記診断項目の異常度である第１異常度を算出し、前記第１異常度に基づいて前記診断項目の異常判定を行う異常判定部と、を含み、
前記異常判定部は、前記第２計測時点について、前記診断項目の計測値と前記診断項目の予測値の確率分布における平均値との差分を当該確率分布の標準偏差で除した値を前記第１異常度として算出する、異常診断装置。

【請求項5】

診断対象システムに関する複数の計測項目のデータに基づいて前記診断対象システムの異常を診断する異常診断方法をコンピュータに実行させる異常診断プログラムであって、
前記診断対象システムが正常に稼働している時間帯における複数の第１計測時点において取得された前記計測項目のデータである正常稼働データに基づいて、目的変数が前記計測項目のうちの一の診断項目であり、説明変数が他の前記計測項目である回帰式の回帰係数の確率分布を表す値を算出する第１ステップと、
前記回帰式、前記第１ステップで算出した前記回帰係数の確率分布を表す値、及び診断の対象となる時間帯における複数の第２計測時点において取得された前記計測項目のデータである診断対象データに基づいて、前記第２計測時点における前記診断項目の予測値の確率分布を表す値を算出する第２ステップと、
前記診断対象データ及び前記第２ステップで算出した前記診断項目の予測値の確率分布を表す値に基づいて前記第２計測時点における前記診断項目の異常度である第１異常度を算出し、前記第１異常度に基づいて前記診断項目の異常判定を行う第３ステップと、を前記コンピュータに実行させ、
前記第３ステップでは、前記第２計測時点について、前記診断項目の計測値と前記診断項目の予測値の確率分布における平均値との差分を当該確率分布の標準偏差で除した値を前記第１異常度として算出する、異常診断プログラム。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、異常診断方法、異常診断装置、及び異常診断プログラムに関する。

【背景技術】

【0002】

異常診断に関する技術として、例えば特許文献１には、プラントに関する複数の計測項目のデータに基づいてプラントの運転状態を監視する方法が記載されている。特許文献１に記載された方法では、プラントが正常に稼働している時間帯に取得された正常稼働データに基づいて算出されるマハラノビス距離（ＭＤ値）が異常度として用いられている。ＭＤ値は、正常範囲からの乖離の大きさを表す値である。また、特許文献２には、重回帰分析により正常稼働データに基づいて計測項目の予測値を算出し、算出した予測値に基づいてガスタービンの性能を診断する方法が記載されている。

【先行技術文献】

【特許文献】

【0003】

【特許文献1】特開２０１０−１９１６４１号公報

【特許文献2】特許第３９２２４２６号公報

【発明の概要】

【発明が解決しようとする課題】

【0004】

上記特許文献１に記載された方法において異常度として用いられるＭＤ値は、全ての計測項目のデータが混合されて得られる値である。そのため、異常に寄与した計測項目の特定、すなわち異常箇所の特定が困難である。また、上記特許文献２に記載された方法を異常診断方法に適用する場合、計測値と算出した予測値との間の乖離の大きさに基づいて異常判定を行うことが考えられる。しかし、乖離がどの程度の大きさになれば異常が発生したと判定すればよいかを決定することが困難となるおそれがある。

【0005】

本発明は、異常判定の容易化及び異常箇所の特定の容易化が図られた異常診断方法、異常診断装置、及び異常診断プログラムを提供することを目的とする。

【課題を解決するための手段】

【0006】

本発明の一側面に係る異常診断方法は、診断対象システムに関する複数の計測項目のデータに基づいて診断対象システムの異常を診断する異常診断方法であって、診断対象システムが正常に稼働している時間帯における複数の第１計測時点において取得された計測項目のデータである正常稼働データに基づいて、目的変数が計測項目のうちの一の診断項目であり、説明変数が他の計測項目である回帰式の回帰係数の確率分布を表す値を算出する第１ステップと、回帰式、第１ステップで算出した回帰係数の確率分布を表す値、及び診断の対象となる時間帯における第２計測時点において取得された計測項目のデータである診断対象データに基づいて、第２計測時点における診断項目の予測値の確率分布を表す値を算出する第２ステップと、診断対象データ及び第２ステップで算出した診断項目の予測値の確率分布を表す値に基づいて第２計測時点における診断項目の異常度である第１異常度を算出し、第１異常度に基づいて診断項目の異常判定を行う第３ステップと、を含む。

【0007】

本発明の一側面に係る異常診断装置は、診断対象システムに関する複数の計測項目のデータに基づいて診断対象システムの異常を診断する異常診断装置であって、診断対象システムが正常に稼働している時間帯における複数の第１計測時点において取得された計測項目のデータである正常稼働データに基づいて、目的変数が計測項目のうちの一の診断項目であり、説明変数が他の計測項目である回帰式の回帰係数の確率分布を表す値を算出する第１算出部と、回帰式、第１算出部により算出された回帰係数の確率分布を表す値、及び診断の対象となる時間帯における第２計測時点において取得された計測項目のデータである診断対象データに基づいて、第２計測時点における診断項目の予測値の確率分布を表す値を算出する第２算出部と、診断対象データ及び第２算出部により算出された診断項目の予測値の確率分布を表す値に基づいて第２計測時点における診断項目の異常度である第１異常度を算出し、第１異常度に基づいて診断項目の異常判定を行う異常判定部と、を含む。

【0008】

本発明の一側面に係る異常診断プログラムは、診断対象システムに関する複数の計測項目のデータに基づいて診断対象システムの異常を診断する異常診断方法をコンピュータに実行させる異常診断プログラムであって、診断対象システムが正常に稼働している時間帯における複数の第１計測時点において取得された計測項目のデータである正常稼働データに基づいて、目的変数が計測項目のうちの一の診断項目であり、説明変数が他の計測項目である回帰式の回帰係数の確率分布を表す値を算出する第１ステップと、回帰式、第１ステップで算出した回帰係数の確率分布を表す値、及び診断の対象となる時間帯における第２計測時点において取得された計測項目のデータである診断対象データに基づいて、第２計測時点における診断項目の予測値の確率分布を表す値を算出する第２ステップと、診断対象データ及び第２ステップで算出した診断項目の予測値の確率分布を表す値に基づいて第２計測時点における診断項目の異常度である第１異常度を算出し、第１異常度に基づいて診断項目の異常判定を行う第３ステップと、をコンピュータに実行させる。

【0009】

これらの異常診断方法、異常診断装置、及び異常診断プログラムでは、診断の対象となる時間帯の第２計測時点における診断項目の異常度である第１異常度が算出される。これにより、診断の対象となる時間帯に診断項目に異常が生じていたかどうかを判定することができる。また、第１異常度は、第２計測時点における診断項目の予測値の確率分布を表す値に基づいて算出される。これにより、計測値と予測値の確率分布との比較に基づいて異常判定を行うことができるので、異常判定を容易に行うことが可能となる。また、診断項目ごとに第１異常度が算出される。これにより、例えば複数の診断項目が設定されている場合でも、異常箇所の特定を容易に行うことが可能となる。このように、上記異常診断方法、異常診断装置、及び異常診断プログラムでは、異常判定の容易化及び異常箇所の特定の容易化が図られている。

【0010】

第１ステップでは、回帰係数の確率分布が正規分布に従うと仮定した場合の回帰係数の確率分布を表す値を算出し、第２ステップでは、回帰係数の確率分布が正規分布に従うと仮定した場合の第２計測時点における診断項目の予測値の確率分布を表す値を算出してもよい。これにより、第２計測時点における診断項目の予測値の確率分布を表す値の算出に要する時間を低減することができる。

【0011】

第３ステップでは、第２計測時点について、診断項目の計測値と診断項目の予測値の確率分布における平均値との差分を当該確率分布の標準偏差で除した値を第１異常度として算出してもよい。これにより、診断項目の予測値の確率分布における標準偏差が小さい場合に、計測値と平均値との差分の大きさが第１異常度に鋭敏に反映される。その結果、診断項目の予測に関するばらつきの大きさが第１異常度に反映されるので、診断項目の異常判定における精度を向上することができる。

【0012】

本発明の一側面に係る異常診断方法は、正常稼働データ、回帰式、及び第１ステップで算出した回帰係数の確率分布を表す値に基づいて、各第１計測時点における診断項目の予測値の確率分布を表す値を算出する第４ステップと、正常稼働データ及び第４ステップで算出した診断項目の予測値の確率分布を表す値に基づいて、各第１計測時点における診断項目の異常度である第２異常度を算出する第５ステップと、を更に含み、第３ステップでは、第１異常度及び第２異常度に基づいて診断項目の異常判定を行ってもよい。これにより、第１異常度だけでなく、診断対象システムが正常に稼働している時間帯の各第１計測時点における診断項目の異常度である第２異常度にも基づいて診断項目の異常判定が行われる。そのため、診断項目の異常判定における精度を一層向上することができる。

【発明の効果】

【0013】

本発明によれば、異常判定の容易化及び異常箇所の特定の容易化が図られる。

【図面の簡単な説明】

【0014】

【図1】本発明の一実施形態に係る異常診断装置の構成を示すブロック図である。

【図2】（ａ）及び（ｂ）は、回帰係数の確率分布を表す値の算出を説明するための図である。

【図3】（ａ）及び（ｂ）は、回帰係数の確率分布を表す値の算出を説明するための図である。

【図4】第１異常度を説明するための図である。

【図5】図１の異常診断装置により実行される処理を示すフローチャートである。

【図6】第１算出処理の詳細を示すフローチャートである。

【図7】第２算出処理の詳細を示すフローチャートである。

【図8】異常判定処理の詳細を示すフローチャートである。

【図9】異常診断ブログラムを説明するための図である。

【発明を実施するための形態】

【0015】

以下、本発明の一実施形態について、図面を参照しつつ詳細に説明する。なお、以下の説明において、同一又は相当要素には同一符号を用い、重複する説明を省略する。

【0016】

図１に示される異常診断装置１は、診断対象となるシステムである診断対象システムの異常を診断する装置である。診断対象システムとしては、例えば、ガスタービン、航空エンジン、又は真空炉等の機械システムが挙げられる。このような診断対象システムには、システムの動作状況を確認するために、システムに関する複数の計測項目を計測するための複数のセンサが取り付けられている。計測項目には、例えば、温度、圧力、又は回転数等が含まれる。異常診断装置１は、これらの計測項目のデータに基づいて診断対象システムの異常を診断する。

【0017】

異常診断装置１は、入力装置２及び出力装置３と通信可能に接続される。入力装置２は、例えば、キーボード、マウス、操作ボタン、又はタッチパネル等により構成され、異常診断装置１のユーザが異常診断装置１に対して操作を加える場合等に用いられる。入力装置２には、例えば、後述する正常稼働データＤ１、診断対象データＤ２、及び設定パラメータＤ３等が外部から入力される。入力装置２は、入力されたデータを異常診断装置１に送信する。出力装置３は、例えばディスプレイ等により構成され、異常診断装置１による診断結果を表示（出力）する。なお、異常診断装置１が入力装置２及び出力装置３を含んで構成されてもよい。

【0018】

異常診断装置１は、物理的には、ＣＰＵ（Central Processing Unit）１０と、記憶装置２０と、を備えている。記憶装置２０は、例えば、ＲＡＭ（Random Access Memory）、フラッシュメモリ、又はハードディスク装置等のデータの読み書きが可能な記録媒体によって構成される。記憶装置２０は、例えば、正常稼働データＤ１、診断対象データＤ２、設定パラメータＤ３、算出データＤ４、診断結果データＤ５、及び異常診断プログラムＰ等を記憶している。

【0019】

正常稼働データＤ１は、診断対象システムが正常に稼働している第１時間帯における複数の第１計測時点において取得された計測項目のデータである。診断対象データＤ２は、診断の対象となる第２時間帯における複数の第２計測時点において取得された計測項目のデータである。第１時間帯の長さと第２時間帯の長さとは、互いに同一であってもよく、異なっていてもよい。第１計測時点の間隔及び第２計測時点の間隔は、任意に設定されてよい。設定パラメータＤ３は、後述の処理に用いられる複数のパラメータを含む。設定パラメータＤ３としては、例えば、いずれも後述される基底関数、超パラメータの最大値、擬似逆行列の算出に関する閾値、収束判定の閾値、収束計算の実行回数の上限値、異常判定の閾値、及び全体異常度の重み係数等が挙げられる。算出データＤ４は、後述の処理において算出される複数の変数の値を含む。診断結果データＤ５は、後述の処理により得られる診断結果を含む。異常診断プログラムＰは、後述の処理を実行するためのプログラムである。

【0020】

ＣＰＵ１０は、機能的には、入力部１１と、第１算出部１２と、第２算出部１３と、異常判定部１４と、出力部１５と、を備えている。ＣＰＵ１０が記憶装置２０に記憶されている異常診断プログラムＰを読み出して実行することにより、入力部１１、第１算出部１２、第２算出部１３、異常判定部１４、及び出力部１５の各機能が実現される。異常診断プログラムＰの具体的な構成については後述する。

【0021】

入力部１１は、入力装置２から送信された正常稼働データＤ１、診断対象データＤ２、及び設定パラメータＤ３を受信し、記憶装置２０に記憶する。

【0022】

第１算出部１２は、正常稼働データＤ１に基づいて、目的変数が計測項目のうちの一の診断項目であり、説明変数が他の計測項目である回帰式の回帰係数の確率分布を表す値を算出する。確率分布を表す値とは、例えば、確率分布の平均若しくは分散等の統計量、又は各確率変数についての確率の値等である。本実施形態では、第１算出部１２は、計測項目のそれぞれが診断項目とされた場合について、説明変数が残りの計測項目の全てである回帰式の回帰係数の確率分布を表す値を算出する。

【0023】

第１算出部１２による回帰係数の確率分布を表す値の算出について更に説明する。正常稼働データＤ１は、

【数1】

と表される。すなわち、正常稼働データＤ１には、Ｎ個の正常稼働データセット{ｘ}が含まれる。正常稼働データセット｛ｘ｝は、同一の第１計測時点において取得された各計測項目であるＤ個の要素｛ｘ_１，ｘ_２，…，ｘ_ｉ，…，ｘ_Ｄ｝を有する列ベクトルである。ここで、要素番号ｉは、１〜Ｄの整数値を取り得る変数である。

【0024】

下記式（２）に示されるように、診断項目ｘ_ｉが目的変数であり、他の計測項目｛ｚ_ｉ｝が説明変数である回帰式が設定される。計測項目｛ｚ_ｉ｝は、(Ｄ−１)個の要素｛ｘ_１，…，ｘ_ｉ−１，ｘ_ｉ＋１，…，ｘ_Ｄ｝を有する列ベクトルである。基底関数｛Φ（ｚ_ｉ）｝は、回帰式の形状を表す関数であり、Ｍ個の要素｛Φ_１（ｚ_ｉ），…，Φ_Ｍ（ｚ_ｉ）｝を有する列ベクトルである。基底関数｛Φ（ｚ_ｉ）｝は、設定パラメータＤ３として記憶装置２０に予め記憶されている。基底関数｛Φ（ｚ_ｉ）｝、すなわち回帰式の形状は、診断項目ｘ_ｉごとに設定されてよく、各診断項目ｘ_ｉに共通であってもよい。回帰式は、任意の形状であってよい。回帰式は、例えば、目的変数が説明変数と回帰係数の線形結合により得られる線形式であってよいし、二次式等を含む非線形な式であってもよい。回帰係数｛ｗ｝は、回帰式の回帰係数であるＭ個の要素｛ｗ_１，…，ｗ_Ｍ｝を有する列ベクトルである。回帰係数｛ｗ｝は、診断項目ｘ_ｉごとに設定される。精度βは、回帰係数｛ｗ｝の確率分布を表す値の算出に用いられる変数である。精度βは、診断項目ｘ_ｉごとに設定される。式（２）では、誤差εが、平均が０で標準偏差が精度βの逆数である正規分布に従うと仮定されている。なお、誤差εは、正規分布に従うと仮定されなくてもよい。

【数2】

【0025】

診断項目ｘ_ｉの予測値の確率分布は、下記式（３）により表される。

【数3】

正常稼働データＤ１

【数4】

について、

【数5】

【数6】

【数7】

とすると、正常稼働データＤ１に関する尤度関数は、下記式（８）により表される。

【数8】

【0026】

異常診断装置１では、下記参考文献１に記載されたスパースベイジアン学習（ＳＢＬ：SparseBayesian Learning）により、回帰係数｛ｗ｝の確率分布を表す値が算出される。ＳＢＬは、回帰係数の確率分布が正規分布に従うと仮定した場合の回帰係数の確率分布を表す値を算出する手法である。ＳＢＬでは、回帰係数｛ｗ｝に対して事前分布

【数9】

が設定される。ここで、行列｛Ａ｝は、下記式（１０）により表される。

【数10】

上記式（９）に示されるように、回帰係数｛ｗ｝の確率分布は、正規分布の形式で表わされる。超パラメータ｛α｝は、回帰係数｛ｗ｝の確率分布に関して設定されたパラメータであるＭ個の要素｛α_１，…，α_Ｍ｝を有する列ベクトルである。超パラメータ｛α｝は、回帰係数｛ｗ｝ごとに設定される。
（参考文献１）Michael E. Tipping, “Sparse Bayesian Learning and the Relevance Vector Machine,” Journal of Machine Learning Research 1 (2001) 211-244

【0027】

回帰係数｛ｗ｝の事後分布は、ベイズの定理により、

【数11】

となる。回帰係数｛ｗ｝の事後分布の平均｛ｍ｝は、下記式（１２）により表され、分散｛Σ｝は、下記式（１３）により表される。平均｛ｍ｝は、Ｍ個の要素を有する列ベクトル、分散｛Σ｝は、Ｍ行Ｍ列の行列である。回帰係数｛ｗ｝の事後分布は、平均｛ｍ｝及び分散｛Σ｝により表される。すなわち、この例では、回帰係数｛ｗ｝の確率分布を表す値は、平均｛ｍ｝及び分散｛Σ｝である。

【数12】

【数13】

【0028】

第１算出部１２は、正常稼働データＤ１に基づいて、回帰係数｛ｗ｝に関する畳み込み積分

【数14】

を最大化する超パラメータ｛α｝及び精度βを算出する。第１算出部１２は、算出した超パラメータ｛α｝及び精度βを用いて平均｛ｍ｝及び分散｛Σ｝を算出する。第１算出部１２は、算出した超パラメータ｛α｝、精度β、平均｛ｍ｝、及び分散｛Σ｝を算出データＤ４として記憶装置２０に記憶させる。

【0029】

上記式（１２）により表される畳み込み積分を最大化する超パラメータ｛α｝及び精度βを算出するアルゴリズムを説明する。このアルゴリズムのインプットは、正常稼働データＤ１、基底関数｛Φ（ｚ_ｉ）｝及び超パラメータの最大値α_ｍａｘである。この超パラメータの最大値α_ｍａｘは、設定パラメータＤ３として記憶装置２０に予め記憶されている。最大値α_ｍａｘは、例えば１０^８程度の値に設定されてよい。アウトプットは、各診断項目ｘ_ｉについての超パラメータ｛α｝、精度β、平均｛ｍ｝、及び分散｛Σ｝である。以下では、１つの診断項目ｘ_ｉについての超パラメータ｛α｝、精度β、平均｛ｍ｝、及び分散｛Σ｝を算出する処理を説明するが、他の診断項目についても同様に算出される。

【0030】

まず、第１算出部１２は、行列｛Ａ｝及び精度βを初期化する（ステップＳ１）。続いて、第１算出部１２は、基底関数｛Φ（ｚ_ｉ）｝を用いて、

【数15】

を算出する（ステップＳ２）。続いて、第１算出部１２は、平均｛ｍ｝及び分散｛Σ｝を上記式（１２）及び式（１３）により算出する（ステップＳ３）。ここで、分散Σの算出に際して逆行列計算が必要になる。この逆行列計算は、擬似逆行列で代用されてよい。擬似逆行列の算出に関する閾値は、設定パラメータＤ３として記憶装置２０に予め記憶されている。続いて、第１算出部１２は、超パラメータ｛α｝及び精度βを所定の更新式により算出（更新）する（ステップＳ４）。超パラメータ｛α｝及び精度βの更新式（算出式）としては、例えば下記参考文献１に記載されたものが用いられてよい。
（参考文献１）Michael E. Tipping, “SparseBayesian Learning and the Relevance Vector Machine,” Journal ofMachine Learning Research 1 (2001) 211-244

【0031】

続いて、第１算出部１２は、超パラメータ｛α｝及び精度βの収束判定を行う（ステップＳ５）。具体的には、第１算出部１２は、更新前の超パラメータ｛α｝の値と更新後の超パラメータ｛α｝の値との差の絶対値が所定の閾値よりも小さく、且つ更新前の精度βと更新後の精度βとの差の絶対値が所定の閾値よりも小さいか否かを判定する。これらの閾値は、設定パラメータＤ３として記憶装置２０に予め記憶されている。第１算出部１２は、判定の結果、超パラメータ｛α｝及び精度βのいずれかが収束していないと判定した場合（ステップＳ５でＮＯ）、ステップＳ２に戻る。一方、第１算出部１２は、超パラメータ｛α｝及び精度βの両方が収束したと判定した場合（ステップＳ５でＹＥＳ）、現在の超パラメータ｛α｝及び精度βにより平均｛ｍ｝及び分散｛Σ｝を算出する（ステップＳ６）。以上の処理により、診断項目ｘ_ｉについての超パラメータ｛α｝、精度β、平均｛ｍ｝、及び分散｛Σ｝が算出される。なお、上記収束計算の繰り返し回数が実行回数の上限値に達した場合には、収束計算が停止されてよい。この上限値は、設定パラメータＤ３として記憶装置２０に予め記憶されている。

【0032】

図２及び図３を参照して、回帰係数｛ｗ｝の確率分布を表す値の算出について更に説明する。図２に示されるように、超パラメータα_ｉが大きいほど、当該超パラメータα_ｉに対応する回帰係数ｗ_ｉの事後分布が値０の近傍により集中する。回帰係数ｗ_ｉが値０に近づくほど、当該回帰係数ｗ_ｉに対応する基底関数Φ_ｉ（ｚ_ｉ）が回帰式に及ぼす影響が小さくなる。図３（ｂ）では、計測項目ｗ_ｉ^（２）に対応する基底関数Φ_ｉ（ｚ_ｉ^（２））が回帰式に及ぼす影響が、図３（ａ）の場合の当該影響よりも小さくなっている。上述したようなアルゴリズムでは、最適な超パラメータ｛α｝の値は、上記式（１４）の畳み込み積分を最大化する超パラメータ｛α｝である。上記式（１４）の畳み込み積分を超パラメータα_ｉの関数とみなしたとき、上記式（１４）は入力データである上記式（４）の値に応じて、単調増加関数又は極大値を1つだけ持つ関数となる。上記式（１４）が単調増加関数である場合、最適な超パラメータα_ｉは無限大になるので、超パラメータ｛α｝の算出過程において、超パラメータα_ｉが大きくなり易い。これにより、後述する診断項目ｘ_ｉの予測値の確率分布を表す値の算出の際に、当該超パラメータα_ｉに対応する基底関数Φ_ｉ（ｚ_ｉ）についての計算が簡易化される。その結果、診断項目ｘ_ｉの予測値の確率分布を表す値の算出に要する時間が低減される。

【0033】

第２算出部１３は、上記回帰式、第１算出部１２により算出された回帰係数｛ｗ｝の確率分布を表す値、及び診断対象データＤ２に基づいて、各第２計測時点における各診断項目ｘ_ｉの予測値の確率分布を表す値を算出する。具体的には、第２算出部１３は、診断対象データＤ２

【数16】

の各第２計測時点ｔにおける各診断項目ｘ_ｉの予測値の確率分布を表す値を、下記式（１７）により算出する。第２算出部１３は、算出した予測値の確率分布を表す値を算出データＤ４として記憶装置２０に記憶させる。

【数17】

診断項目ｘ_ｉの予測値の確率分布の予測平均μは、下記式（１８）により表され、予測分散σ^２は、下記式（１９）により表される。各診断項目ｘ_ｉの予測値の確率分布は、予測平均μ及び予測分散σ^２により表される。すなわち、この例では、回帰係数｛ｗ｝の確率分布を表す値は、予測平均μ及び予測分散σ^２である。

【数18】

【数19】

【0034】

異常判定部１４は、診断対象データＤ２及び第２算出部１３により算出された各診断項目ｘ_ｉの予測値の確率分布を表す値に基づいて、各第２計測時点ｔにおける診断項目ｘ_ｉの異常度である第１異常度ｄ_ｉを算出する。具体的には、異常判定部１４は、各第２計測時点ｔにおける各診断項目ｘ_ｉの第１異常度ｄ_ｉを、下記式（２０）により算出する。すなわち、本実施形態では、異常判定部１４は、各第２計測時点ｔについて、診断項目ｘ_ｉの計測値と予測平均μとの差分を当該確率分布の予測標準偏差σで除した値を第１異常度ｄ_ｉとして算出する。異常判定部１４は、算出した第１異常度ｄ_ｉを算出データＤ４として記憶装置２０に記憶させる。

【数20】

【0035】

異常判定部１４は、第１異常度ｄ_ｉに基づいて各第２計測時点ｔにおける診断項目ｘ_ｉの異常判定を行う。異常判定部１４は、例えば、各第２計測時点ｔおける第１異常度ｄ_ｉが所定の閾値（例えば、値３）よりも大きいか否かを判定する。一の第２計測時点ｔにおいて第１異常度ｄ_ｉが当該閾値よりも大きいと判定した場合、当該第２計測時点ｔにおいて診断項目ｘ_ｉに異常が発生していたと判定する。一方、第１異常度ｄ_ｉが当該閾値以下であると判定した場合、当該第２計測時点ｔにおいて診断項目ｘ_ｉに異常が発生していなかったと判定する。この異常判定に用いられる閾値は、設定パラメータＤ３として記憶装置２０に予め記憶されている。異常判定部１４は、異常判定の結果を診断結果データＤ５として記憶装置２０に記憶させる。

【0036】

異常判定部１４は、第１異常度ｄ_ｉに基づいて各第２計測時点ｔにおける全体異常度を算出する。一の第２計測時点ｔにおける全体異常度は、例えば、当該第２計測時点ｔにおける各診断項目ｘ_ｉの重み付き和である。この算出に用いられる重み係数は、設定パラメータＤ３として記憶装置２０に予め記憶されている。重み係数は、例えば使用者によって設定される。異常判定部１４は、全体異常度を診断結果データＤ５として記憶装置２０に記憶させる。

【0037】

異常判定部１４は、第１異常度ｄ_ｉに基づいて、異常箇所を特定するための情報を算出（生成）する。この異常箇所の特定情報は、例えば、各第２計測時点ｔにおいて第１異常度ｄ_ｉが最も大きかった診断項目ｘ_ｉと、当該診断項目ｘ_ｉの第１異常度ｄ_ｉとを示す情報である。異常箇所の特定情報は、各第２計測時点ｔにおいて第１異常度ｄ_ｉが大きかった順に並べられた複数の診断項目ｘ_ｉと、当該診断項目ｘ_ｉの第１異常度ｄ_ｉとを示す情報であってもよい。

【0038】

図４を参照して、第１異常度ｄ_ｉについて更に説明する。図４に示されるように、第１異常度ｄ_ｉは、診断項目ｘ_ｉの計測値が予測平均μから予測標準偏差σの何倍離れているかを表す。図４では、一例として、予測平均μを示す線分、及び予測平均μから予測標準偏差σの３倍離れた線分が示されている。異常診断装置１では、診断項目ｘ_ｉの予測値の確率分布における予測標準偏差σが小さい場合に、診断項目ｘ_ｉの計測値と予測平均μとの差分の大きさが第１異常度ｄ_ｉに鋭敏に反映される。その結果、診断項目ｘ_ｉの予測に関するばらつきの大きさが第１異常度ｄ_ｉに反映されるので、診断項目ｘ_ｉの異常判定における精度を向上することができる。

【0039】

出力部１５は、算出データＤ４及び診断結果データＤ５を出力する出力手段として機能する。出力部１５は、記憶装置２０から算出データＤ４及び診断結果データＤ５を読み出し、出力装置３に送信する。出力装置３は、出力部１５から算出データＤ４及び診断結果データＤ５を受信すると、受信した算出データＤ４及び診断結果データＤ５を診断結果として出力する。

【0040】

次に、図５〜図８を参照して、異常診断装置１が行う異常診断方法の一例を説明する。図５に示されるように、この異常診断方法では、まず、第１算出部１２が第１算出処理を実行する（ステップＳ１１、第１ステップ）。続いて、第２算出部１３が第２算出処理を実行する（ステップＳ１２、第２ステップ）。続いて、異常判定部１４が異常判定処理を実行する（ステップＳ１３、第３ステップ）。続いて、出力部１５が出力処理を実行する（ステップＳ１４）。これらステップＳ１１〜Ｓ１４の処理は、任意のタイミングで実行されてよく、診断対象システムの停止中及び稼働中のいずれに実施されてもよい。

【0041】

図６を参照して第１算出処理（学習フェーズ）を説明する。まず、第１算出部１２は、記憶装置２０から正常稼働データＤ１を読み込む（ステップＳ２１）。続いて、第１算出部１２は、第１算出部１２は、正常稼働データＤ１の前処理を行う（ステップＳ２２）。具体的には、第１算出部１２は、外れ値除去処理及び正規化処理を行う。外れ値除去処理では、第１算出部１２は、例えば、予め設定された範囲から計測値が外れている計測項目のデータが除去する。この範囲は、設定パラメータＤ３として記憶装置２０に予め記憶されていてよい。正規化処理では、第１算出部１２は、例えば、各第１計測時点における各計測項目の計測値を、平均が値０で分散が値１となるように正規化する。このとき、第１算出部１２は、正常稼働データＤ１の平均及び分散を設定パラメータＤ３として記憶装置２０に記憶させてよい。続いて、第１算出部１２は、各診断項目ｘ_ｉについて、回帰式の回帰係数｛ｗ｝の確率分布を表す値を算出する（ステップＳ２３）。続いて、第１算出部１２は、算出した回帰係数｛ｗ｝の確率分布を表す値を算出データＤ４として記憶装置２０に記憶させる（ステップＳ２４）。

【0042】

図７を参照して第２算出処理（診断フェーズ）を説明する。まず、第２算出部１３は、診断対象データＤ２及び第１算出部１２により算出された回帰係数｛ｗ｝の確率分布を表す値を記憶装置２０から読み込む（ステップＳ３１）。続いて、第２算出部１３は、診断対象データＤ２の前処理を行う（ステップＳ３１）。具体的には、第２算出部１３は、外れ値除去処理及び正規化処理を行う。外れ値除去処理では、第２算出部１３は、予め設定された範囲から計測値が外れている計測項目のデータを除去する。この範囲は、例えば、上述したステップＳ２２の範囲と同一である。正規化処理では、上記ステップＳ２２において記憶装置２０に記憶された平均及び分散が用いられる。続いて、第２算出部１３は、各診断項目ｘ_ｉについて、各第２計測時点ｔにおける各診断項目ｘ_ｉの予測値の確率分布を表す値を算出する（ステップＳ３３）。続いて、第２算出部１３は、算出した予測値の確率分布を表す値を算出データＤ４として記憶装置２０に記憶させる（ステップＳ３４）。

【0043】

図８を参照して異常判定処理（異常判定フェーズ）を説明する。まず、異常判定部１４は、診断対象データＤ２及び第２算出部１３により算出された各診断項目ｘ_ｉの予測値の確率分布を表す値を記憶装置２０から読み込む（ステップＳ４１）。続いて、異常判定部１４は、各第２計測時点ｔにおける各診断項目ｘ_ｉの第１異常度ｄ_ｉを算出する（ステップＳ４２）。続いて、異常判定部１４は、第１異常度ｄ_ｉに基づいて各第２計測時点ｔにおける各診断項目ｘ_ｉの異常判定を行う（ステップＳ４３）。続いて、異常判定部１４は、全体異常度を算出する（ステップＳ４４）。続いて、異常判定部１４は、異常箇所の特定情報を算出する（ステップＳ４５）。続いて、異常判定部１４は、第１異常度ｄ_ｉ、異常判定の結果、全体異常度、及び異常箇所の特定情報を診断結果データＤ５として記憶装置２０に記憶させる。

【0044】

出力処理では、出力部１５は、記憶装置２０から診断結果データＤ５を読み出し、出力装置３に送信する（ステップＳ１４）。

【0045】

次に、図９を参照しながら、コンピュータを異常診断装置１として機能させるための異常診断プログラムＰを説明する。

【0046】

異常診断プログラムＰは、メインモジュールＰ１０、入力モジュールＰ１１、第１算出モジュールＰ１２、第２算出モジュールＰ１３、異常判定モジュールＰ１４、及び出力モジュールＰ１５を備える。メインモジュールＰ１０は、異常診断に係る処理を統括的に制御する部分である。入力モジュールＰ１１、第１算出モジュールＰ１２、第２算出モジュールＰ１３、異常判定モジュールＰ１４、及び出力モジュールＰ１５を実行することにより実現される機能はそれぞれ、上記各実施形態における入力部１１、第１算出部１２、第２算出部１３、異常判定部１４、及び出力部１５の機能と同様である。すなわち、異常診断プログラムＰは、第１算出処理、第２算出処理、異常判定処理、及び出力処理をコンピュータに実行させる。

【0047】

異常診断プログラムＰは、ＣＤ−ＲＯＭ（Compact Disk Read Only Memory）、ＤＶＤ−ＲＯＭ（DigitalVersatile Disk Read Only Memory）、又は半導体メモリ等の有形の記録媒体に固定的に記録された状態で提供されてよい。あるいは、異常診断プログラムＰは、搬送波に重畳されたデータ信号として通信ネットワークを介して提供されてもよい。

【0048】

以上説明した本実施形態の異常診断方法、異常診断装置１、及び異常診断プログラムＰでは、診断の対象となる第２時間帯の各第２計測時点ｔにおける診断項目ｘ_ｉの異常度である第１異常度ｄ_ｉが算出される。これにより、診断の対象となる第２時間帯に各診断項目ｘ_ｉに異常が生じていたかどうかを判定することができる。また、第１異常度ｄ_ｉは、各第２計測時点ｔにおける各診断項目ｘ_ｉの予測値の確率分布を表す値に基づいて算出される。これにより、計測値と予測値の確率分布との比較に基づいて異常判定を行うことができるので、異常判定を容易に行うことが可能となる。また、診断項目ｘ_ｉごとに第１異常度ｄ_ｉが算出される。これにより、上記実施形態のように複数の診断項目ｘ_ｉが設定されている場合でも、異常箇所の特定を容易に行うことが可能となる。このように、本実施形態の異常診断方法、異常診断装置１、及び異常診断プログラムＰでは、異常判定の容易化及び異常箇所の特定の容易化が図られている。

【0049】

上述したように、上記特許文献１に記載された方法において異常度として用いられるＭＤ値は、全ての計測項目のデータが混合されて得られる値である。そのため、異常度の意味が使用者にとって直感的に分かり難い。これに対して、本実施形態の異常診断方法、異常診断装置１、及び異常診断プログラムＰでは、第１異常度ｄ_ｉは計測値と予測値の確率分布との比較に基づくので、第１異常度ｄ_ｉが使用者に呈示されることにより、使用者は、計測値と予測値との間の乖離の大きさを直感的に把握することが可能となる。例えば、上述したように、本実施形態では、第１異常度ｄ_ｉは、診断項目ｘ_ｉの計測値が予測平均μから予測標準偏差σの何倍離れているかが表す。そのため、第１異常度ｄ_ｉが使用者に呈示されることにより、使用者は、計測値と予測値との間の乖離の大きさを直感的に把握することができる。なお、例えば、診断項目ｘ_ｉの予測値の確率分布と計測値とを使用者に呈示することによっても、使用者は、計測値と予測値との間の乖離の大きさを直感的に把握することができる。

【0050】

第１算出処理では、回帰係数｛ｗ｝の確率分布が正規分布に従うと仮定した場合の回帰係数｛ｗ｝の確率分布を表す値が算出され、第２算出処理では、回帰係数｛ｗ｝の確率分布が正規分布に従うと仮定した場合の各第２計測時点ｔにおける各診断項目ｘ_ｉの予測値の確率分布を表す値が算出される。これにより、各第２計測時点ｔにおける各診断項目ｘ_ｉの予測値の確率分布を表す値の算出に要する時間を低減することができる。更に、診断項目ｘ_ｉの予測値の確率分布を表す値の算出過程において当該確率分布の予測平均μ及び予測標準偏差σが算出されるので、第１異常度ｄ_ｉの算出が容易化される。

【0051】

異常判定処理では、各第２計測時点ｔについて、各診断項目ｘ_ｉの計測値と各診断項目ｘ_ｉの予測値の確率分布における予測平均μとの差分を当該確率分布の予測標準偏差σで除した値が第１異常度ｄ_ｉとして算出される。これにより、診断項目ｘ_ｉの予測値の確率分布における予測標準偏差σが小さい場合に、計測値と予測平均μとの差分の大きさが第１異常度ｄ_ｉに鋭敏に反映される。その結果、予測標準偏差σの大きさ、すなわち診断項目ｘ_ｉの予測に関するばらつきの大きさが第１異常度ｄ_ｉに反映されるので、診断項目ｘ_ｉの異常判定における精度を向上することができる。

【0052】

本実施形態の異常診断方法、異常診断装置１、及び異常診断プログラムＰでは、第１異常度ｄ_ｉが使用者に呈示されることにより、使用者は、各第２計測時点ｔにおける各診断項目ｘ_ｉの異常の程度を把握することができる。異常判定の結果が使用者に呈示されることにより、使用者は、各第２計測時点ｔにおける各診断項目ｘ_ｉの異常の有無を把握することができる。全体異常度が使用者に呈示されることにより、使用者は、各第２計測時点ｔにおける診断対象システム全体の異常の程度を把握することができる。また、全体異常度が呈示されることにより、使用者が監視すべき異常度の数を減らすことができる。また、使用者は、重み係数を調整することにより、重点的に監視すべき診断項目ｘ_ｉと監視の優先度が低い診断項目ｘ_ｉとを選択することができる。異常箇所の特定情報が使用者に呈示されることにより、使用者は、診断対象システムのどの箇所において第１異常度ｄ_ｉが高くなっていたかを把握することができる。また、異常箇所の特定情報が使用者に呈示されることにより、使用者は、異常が発生している場合に、異常箇所を特定することができる。

【0053】

以上、本発明の一実施形態について説明したが、本発明は、上記実施形態に限られない。上記実施形態において、第１算出部１２は、正常稼働データＤ１、回帰式、及びステップＳ２３で算出した回帰係数｛ｗ｝の確率分布を表す値に基づいて、各第１計測時点における各診断項目ｘ_ｉの予測値の確率分布を表す値を算出してもよい（第４ステップ）。この算出は、ステップＳ３３の算出と同様に行うことができる。第４ステップは、例えばステップＳ２３の後に実行される。第１算出部１２は、正常稼働データＤ１、及び第４ステップで算出した各診断項目ｘ_ｉの予測値の確率分布を表す値に基づいて、各第１計測時点における診断項目ｘ_ｉの異常度である第２異常度ｅ_ｉを算出してもよい（第５ステップ）。この算出は、ステップＳ４２の算出と同様に行うことができる。第５ステップは、第４ステップの後に実行される。ステップＳ２４では、第１算出部１２は、上述した処理に加えて、算出した予測値の確率分布を表す値及び第２異常度ｅ_ｉを算出データＤ４として記憶装置２０に記憶させてよい。ステップＳ４３では、異常判定部１４は、第１異常度ｄ_ｉ及び第２異常度ｅ_ｉに基づいて各第２計測時点ｔにおける各診断項目ｘ_ｉの異常判定を行ってもよい。異常判定部１４は、例えば、各第１計測時点における診断項目ｘ_ｉの異常度である第２異常度ｅ_ｉの平均及び分散に基づいて、第１異常度ｄ_ｉの異常判定に用いる閾値を決定してもよい。これにより、第１異常度ｄ_ｉだけでなく、診断対象システムが正常に稼働している第１時間帯の各第１計測時点における診断項目ｘ_ｉの異常度である第２異常度ｅ_ｉにも基づいて各診断項目ｘ_ｉの異常判定が行われる。その結果、診断項目ｘ_ｉの異常判定における精度を一層向上することができる。

【0054】

回帰式は、目的変数が計測項目のうちの一の診断項目であり、説明変数が他の計測項目の一部であるものであってもよい。例えば、相関が無いことが既知である計測項目については、説明変数から除外されてもよい。すなわち、説明変数は、他の計測項目の少なくとも１つであればよい。上記実施形態では、計測項目の全てが診断項目とされたが、計測項目の一部が診断項目とされてもよい。つまり、診断項目とされない計測項目が存在してもよい。診断項目とされる計測項目は、例えば使用者によって選択されてよい。診断項目は、１つであってもよい。第２計測時点は、１つの時点であってもよい。

【0055】

回帰係数｛ｗ｝の確率分布の算出にＳＢＬが用いられなくてもよい。つまり、回帰係数の確率分布が正規分布に従うと仮定した場合の回帰係数の確率分布を表す値及び予測値の確率分布を表す値が算出されなくてもよい。この場合、例えば、マルコフ連鎖モンテカルロ法等により回帰係数の確率分布を表す値及び予測値の確率分布を表す値が数値的に算出されてもよい。この場合、回帰係数の確率分布を表す値及び予測値の確率分布を表す値は、各確率変数についての確率の値であってよい。また、回帰係数の確率分布を表す値及び予測値の確率分布を表す値は、平均、分散、及び各確率変数についての確率の値に限られず、確率分布を表現する他の統計量等であってもよい。また、第１異常度ｄ_ｉは、例えば、予測値の確率分布の代表値（例えば、平均値、最頻値、又は中央値等）と計測値との差分を定数倍したものであってよい。第１異常度ｄ_ｉは、診断項目ｘ_ｉの計測値と予測平均μとの差分を定数倍したものであってもよい。第１異常度ｄ_ｉに基づいて診断項目ｘ_ｉの異常判定が行われればよく、第１異常度ｄ_ｉは使用者に呈示されなくてもよい。また、全体異常度は算出されなくてもよく、異常箇所の特定情報は算出されなくてもよい。

【0056】

基底関数は、カーネル法を用いてカーネル関数として設定されてもよい。この場合、カーネル関数は、線形カーネル、多項式カーネル、又はガウシアンカーネル等であってよい。ＳＢＬにおいて基底関数としてカーネル関数が設定される場合の学習法は、関連ベクトルマシン（Relevance Vector Machines：ＲＶＭｓ）と呼ばれる場合がある。

【符号の説明】

【0057】

１ 異常診断装置
１２ 第１算出部
１３ 第２算出部
１４ 異常判定部
Ｄ１ 正常稼働データ
Ｄ２ 診断対象データ
ｄ 第１異常度
ｅ 第２異常度
Ｐ 異常診断プログラム
ｘ 診断項目
ｗ 回帰係数
μ 予測平均
σ 予測標準偏差
σ 予測標準偏差
Φ 基底関数

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】