JP2019045484A - 状態監視方法および状態監視装置 - Google Patents

Abstract

【課題】周波数帯域の最適化を予め行なう必要がなく、かつ被試験対象物の異常の有無の判定精度の高い状態監視方法を提供する。【解決手段】状態監視方法は、測定データに対して、互いに異なる中心周波数を有する複数のバンドパスフィルタの各々を用いたフィルタ処理を行なうことにより、複数の第１波形を生成する第１工程と、複数の第１波形の各々に対して第１特徴量を算出する第２工程と、中心周波数を横軸とし、第１特徴量を縦軸とするグラフに、複数のバンドパスフィルタの各々について、当該バンドパスフィルタの中心周波数と、当該バンドパスフィルタを用いて生成された時間波形の第１特徴量とをプロットすることにより得られる第２波形の特徴を示す第２特徴量を算出する第３工程と、テストデータに対して第１〜第３工程を行なうことにより算出された第２特徴量に基づいて、被試験対象物が異常か否かを判定する第４工程とを備える。【選択図】図６

Description

本発明は、回転機械等の被試験装置の状態監視方法および状態監視装置に関する。

従来、回転機械等の被試験対象物では、各種センサが測定した物理量の測定データに基づいて、被試験対象物の状態を監視している。具体的には、測定データによって示される波形の特徴を示す特徴量が算出され、算出された特徴量を用いて被試験対象物の異常の有無が判定される。

異常の有無の判定精度は、特徴量の異常に対する感度に依存する。そのため、様々な特徴量を用いた状態監視方法が開発されている。たとえば、特開平１０−２７４５５８号公報（特許文献１）には、波形データのスペクトルの時間変化から求めたピークが生じる周波数と、当該ピークが生じる時間間隔との組を用いて異常の有無を判定する異常診断方法が開示されている。特開２００６−３００８９５号公報（特許文献２）には、クラスタリングマップの各ニューロンと設備の運転時の信号から抽出した周波数成分に対応するニューロンとの最小距離を用いて設備の異常の有無を判定する設備監視方法が開示されている。

特開平１０−２７４５５８号公報 特開２００６−３００８９５号公報

しかしながら、被試験対象物の異常の有無の判定の精度をさらに向上させることが望まれている。

この発明は、上記の課題を解決するためになされたものであって、その目的は、被試験対象物の異常の有無の判定精度の高い状態監視方法および状態監視装置を提供することである。

この発明は、被試験対象物に設置したセンサから取得した測定データを用いて被試験対象物の状態を監視する状態監視方法である。状態監視方法は、第１〜第４工程を備える。第１工程では、測定データに対して、互いに異なる中心周波数を有する複数のバンドパスフィルタの各々を用いたフィルタ処理を行なうことにより、複数のバンドパスフィルタにそれぞれ対応する複数の第１波形を生成する。第２工程では、複数の第１波形の各々に対して第１特徴量を算出する。第３工程では、中心周波数を横軸とし、第１特徴量を縦軸とするグラフに、複数のバンドパスフィルタの各々について、当該バンドパスフィルタの中心周波数と、当該バンドパスフィルタを用いて生成された時間波形に対して算出された第１特徴量とをプロットすることにより得られる第２波形の特徴を示す第２特徴量を算出する。第４工程では、測定データのうち被試験対象物の診断時にセンサから得られたテストデータに対して第１〜第３工程を行なうことにより算出された第２特徴量に基づいて、被試験対象物が異常か否かを判定する。

好ましくは、複数の第１波形は時間波形である。第１特徴量は、複数の第１波形の中の当該第１特徴量に対応する第１波形と、対応する第１波形をフーリエ変換して得られる第３波形と、第３波形をフーリエ変換して得られる第４波形とのいずれか１つの波形の特徴を示す。

好ましくは、複数の第１波形は時間波形である。第２工程は、複数の第１波形の中の第１特徴量に対応する第１波形と、対応する第１波形をフーリエ変換して得られる第３波形と、第３波形をフーリエ変換して得られる第４波形との少なくとも２つの波形の各々に対して、当該波形の特徴を示す第３特徴量を算出する工程と、少なくとも２つの波形の各々に対して算出した第３特徴量から第１特徴量を算出する工程とを含む。

好ましくは、状態監視方法は、測定データのうち被試験対象物の正常時にセンサから得られた正常データに対して第１〜第３工程を行なうことにより得られた第２特徴量を学習データとして機械学習を行なうことにより、被試験対象物が異常か否かを判定するためのアルゴリズムを決定する第５工程をさらに備える。第４工程では、アルゴリズムに従って前記被試験対象物が異常か否かを判定する。

好ましくは、第２特徴量は、第２波形の実効値、最大値、波高率、尖度、歪度、変動係数および寄与率のいずれかである。

この発明は、他の局面においては、上記のいずれかの方法を用いて、被試験対象物を診断する、状態監視装置である。

本発明の状態監視方法は、被試験装置の異常の有無の判定精度を高めることができる。

本実施の形態に係る状態監視装置の構成を示すブロック図である。 データ演算部の詳細を示すブロック図である。 図１のデータ取得部が行なう処理を説明するためのフローチャートである。 図１のデータ取得部が取得した測定データの一例を示す図である。 図１のデータ取得部が取得した測定データの別の例を示す図である。 図２のフィルタ処理部および特徴量算出部が行なう処理を説明するためのフローチャートである。 図４の測定データから生成された３つの時間波形を示す。 図５の測定データから生成された３つの時間波形を示す。 中心周波数波形の例を示す図である。 図２の学習部が行なう処理を説明するためのフローチャートである。 図２の判定部が行なう処理を説明するためのフローチャートである。 実施例と比較例との異常度の変化を示す図である。 検証実験２における、実施例のデータ演算方法により算出された第２特徴量のｔ検定結果を示す図である。 比較例のデータ演算方法により算出された比較用特徴量のｔ検討結果を示す図である。 変形例の第１特徴量の一例を示す図である。 検証実験３における第２特徴量のｔ検定結果を示す図である。 検証実験４における第２特徴量のｔ検定結果を示す図である。

以下、本発明の実施の形態について図面を参照しつつ説明する。なお、以下の図面において同一または相当する部分には同一の参照番号を付し、その説明は繰返さない。

［状態監視装置の基本構成］
図１は、本実施の形態に係る状態監視装置の構成を示すブロック図である。図１を参照して、状態監視装置１００は、被試験装置１０に設置された振動センサ２０から振動加速度信号を受けて、被試験装置１０の状態を監視し、異常の有無を判定する。被試験装置１０は、例えば工場や発電所などに設置された回転機械を含む設備であり、振動センサ２０は、回転時に生じる異常振動を検出することができる。なお、本実施の形態では、振動加速度信号を出力する振動センサ２０を例示するが、設備の運転状況を確認できる出力信号であれば振動センサ以外の検出信号であっても良い。たとえば、速度、変位、音、ＡＥ（Acoustic Emission）、温度、負荷トルク、モータ電力等の物理量を示す信号を出力するセンサを振動センサ２０に代えて使用してもよい。

状態監視装置１００は、Ａ／Ｄコンバータ１１０と、データ取得部１２０と、記憶装置１３０と、データ演算部１４０と、表示部１５０とを含む。

Ａ／Ｄコンバータ１１０は、振動センサ２０の振動加速度信号を受け、デジタル形式の信号に変換して出力する。データ取得部１２０は、Ａ／Ｄコンバータ１１０からデジタル形式の振動加速度信号を受け、ノイズ除去処理を行なって、所定時間（たとえば２０秒）の振動加速度の変化を示す測定データを生成する。データ取得部１２０は、生成した測定データを記憶装置１３０に記録する。データ演算部１４０は、記憶装置１３０から正常時に測定しておいた測定データ（正常データ）を読み出して、異常の有無を判定するためのアルゴリズムを決定したり、当該アルゴリズムを用いてテスト時に測定した測定データ（テストデータ）から被試験装置１０の異常の有無を判定したりする。データ演算部１４０は、異常の有無を判定した場合、表示部１５０に判定結果を表示させる。

図２は、データ演算部の詳細を示すブロック図である。データ演算部１４０は、フィルタ処理部１４１と、特徴量算出部１４２と、学習部１４３と、しきい値記憶部１４４と、判定部１４６とを含む。

フィルタ処理部１４１は、記憶装置１３０に格納された測定データに対して、互いに異なる中心周波数を有する複数のバンドパスフィルタ（ＢＰＦ（Band-pass Filter））の各々を用いたフィルタ処理を行なうことにより、複数のバンドパスフィルタにそれぞれ対応する複数の時間波形（第１波形）を生成する。

特徴量算出部１４２は、フィルタ処理部１４１により生成された複数の時間波形の各々に対して第１特徴量を算出する。さらに、特徴量算出部１４２は、測定データごとに、バンドパスフィルタの中心周波数の低い順に第１特徴量を並べることによって得られる波形（以下、「中心周波数波形」という）（第２波形）を生成する。特徴量算出部１４２は、測定データごとに、当該測定データに対応する中心周波数波形の特徴を示す第２特徴量を算出する。

学習部１４３は、被試験装置１０が正常時（初期状態など）に取得された複数の測定データ（正常データ）に対してそれぞれ算出された複数の第２特徴量を学習データとして機械学習を行なうことにより、被試験装置１０が異常か否かを判定するためのアルゴリズムを決定する。学習部１４３は、決定したアルゴリズムをしきい値記憶部１４４に記憶させる。

判定部１４６は、しきい値記憶部１４４が記憶するアルゴリズムに従って、被試験装置１０の診断時に得られた測定データ（テストデータ）に対して算出された第２特徴量を用いて被試験装置１０の異常の有無を判定する。

［データ取得部の処理］
図３は、図１のデータ取得部が行なう処理を説明するためのフローチャートである。データ取得部１２０は、ステップＳ１において、デジタル形式の振動加速度信号を受信し、所定時間（たとえば２０秒）の振動加速度の変化を示すデータを生成する。ステップＳ２において、データ取得部１２０は、生成したデータに対して、ローパスフィルタ、バンドパスフィルタ、ハイパスフィルタ等のうち適切なフィルタ処理を施して、基本的なノイズを除去した測定データを生成する。ステップＳ３において、データ取得部１２０は、生成した測定データを記憶装置１３０に記憶する。

なお、データ取得部１２０は、被試験装置１０の初期状態、修理完了時などの正常動作することが分かっている時に取得した測定データを正常データとして記憶装置１３０に記憶する。データ取得部１２０は、被試験装置１０の使用中に診断を行ないたい時に取得した測定データをテストデータとして記憶装置１３０に記憶する。データ取得部１２０は、タイマーなどで指定された時間（たとえば２時間ごと）に自動的にテストデータを取得する。

図４は、データ取得部１２０が取得した測定データの一例を示す図である。図５は、データ取得部１２０が取得した測定データの別の例を示す図である。図４および図５には、被試験装置１０が軸受であるときの振動加速度の時間変化が示される。図４には、軸受の軌道面に、直径１．３５ｍｍの円筒状の損傷を人工的に設けたときの測定データが示され、図５には、軸受の軌道面にさらに大きい損傷を人工的に設けたときの測定データが示される。

［フィルタ処理部および特徴量算出部の処理］
図６は、図２のフィルタ処理部１４１および特徴量算出部１４２が行なう処理を説明するためのフローチャートである。

フィルタ処理部１４１は、ステップＳ１１において、測定データに対して、互いに中心周波数の異なる複数のバンドパスフィルタの各々を用いたフィルタ処理を行なう。複数のバンドパスフィルタのバンド幅は特に限定されるものではなく、一定（たとえば１ｋＨｚ）であってもよいし、互いに異なっていてもよい。複数のバンドパスフィルタの中心周波数は、一定間隔（たとえば１ｋＨｚ）であってもよいし、変則的な間隔であってもよい。さらに、複数のバンドパスフィルタの中の１つのバンドパスフィルタの帯域は、残りのいずれのバンドパスフィルタの帯域とも重なっていなくてもよいし、残りのいずれかのバンドパスフィルタの帯域と重なっていてもよい。フィルタ処理部１４１は、複数のバンドパスフィルタにそれぞれ対応する複数の時間波形を生成する。

図７は、図４の測定データから生成された３つの時間波形を示す。図８は、図５の測定データから生成された３つの時間波形を示す。図７および図８の（ａ）は、０〜１ｋＨｚを通過帯域とするバンドパスフィルタを用いて生成された時間波形を示す。図７および図８の（ｂ）は、１〜２ｋＨｚを通過帯域とするバンドパスフィルタを用いて生成された時間波形を示す。図７および図８の（ｃ）は、９〜１０ｋＨｚを通過帯域とするバンドパスフィルタを用いて生成された時間波形を示す。図７および図８に示されるように、周波数帯域によって時間波形が異なる。

図６に戻って、次にステップＳ１２において、特徴量算出部１４２は、複数の時間波形の各々に対して第１特徴量を算出する。具体的には、第１特徴量は、対応する時間波形の特徴を示す時間領域特徴量と、当該時間波形をフーリエ変換して得られる周波数波形（第３波形）の特徴を示す周波数領域特徴量と、当該周波数波形を時間波形と見なしてフーリエ変換して得られるケフレンシ波形（第４波形）の特徴を示すケフレンシ領域特徴量とのいずれかである。第１特徴量は、たとえば、波形の実効値、最大値、波高率、尖度、歪度、ピークトゥピーク値、変動係数、標準偏差、分散、相関係数、寄与率、エンベロープ処理後の波形の実効値である。もしくは、第１特徴量は、これらパラメータのいずれかの値の常用対数であってもよい。相関係数とは、２つの確率変数（ここでは、時間、周波数またはケフレンシと加速度）の分布法則の関係（たとえば線形関係）の程度を示す係数である。寄与率は、相関係数の二乗値である。

特徴量算出部１４２は、１つの時間波形に対して、１種類の第１特徴量を算出してもよいし、複数種類の第１特徴量を算出してもよい。たとえば、特徴量算出部１４２は、時間波形、周波数波形およびケフレンシ波形の３つの波形の各々に対して、実効値および尖度を算出してもよい。この場合、特徴量算出部１４２は、１つの時間波形に対して、３×２＝６種類の第１特徴量を算出することになる。

次にステップＳ１３において、特徴量算出部１４２は、バンドパスフィルタの中心周波数を変数とする中心周波数波形を生成する。つまり、特徴量算出部１４２は、バンドパスフィルタの中心周波数を横軸とし、第１特徴量を縦軸とするグラフに、複数のバンドパスフィルタの各々について、当該バンドパスフィルタの中心周波数と、当該バンドパスフィルタに対応する時間波形に対して算出された第１特徴量とを示す点をプロットすることにより中心周波数波形を生成する。

図９は、中心周波数波形の例を示す図である。図９には、第１特徴量が周波数波形の最大値であるときの中心周波数波形が示される。図９において、「損傷サイズ１．３５ｍｍ」と記載された中心周波数波形は、図４の測定データから生成された波形を示し、「損傷なし」と記載された中心周波数波形は、被試験装置１０である軸受に損傷が見られないときの測定データから生成された波形を示す。図９に示されるように、損傷がないときの中心周波数波形の形状と、損傷があるときの中心周波数波形の形状とが異なる。このことは、中心周波数波形の形状が損傷の有無に影響されることを意味している。

特徴量算出部１４２は、複数種類の第１特徴量を算出した場合、当該複数種類の第１特徴量のそれぞれに対応する複数の中心周波数波形を生成する。

図６に戻って、次にステップＳ１４において、特徴量算出部１４２は、中心周波数波形の特徴を示す第２特徴量を算出する。第２特徴量は、たとえば、中心周波数波形の実効値、最大値、波高率、尖度、歪度、ピークトゥピーク値、変動係数、標準偏差、分散、相関係数および寄与率、ならびにエンベロープ処理後の中心周波数波形の実効値のいずれかである。

相関係数とは、２つの確率変数（ここでは、中心周波数と第１特徴量）の分布法則の関係（たとえば線形関係）の程度を示す係数である。寄与率は、相関係数の二乗値である。相関係数として標本相関係数を用いることができる。標本相関係数ρは、バンドパスフィルタの中心周波数の値をｘ_ｉ（ｉ＝１、２、・・・、ｎ）、当該バンドパスフィルタを用いて生成された時間波形に対応する第１特徴量の値をｙ_ｉ（ｉ＝１、２、・・・、ｎ）とするとき、以下の式１で算出される。式１において、ｘ_ａｖｅは、ｘ_ｉ（ｉ＝１、２、・・・、ｎ）の平均値であり、ｙ_ａｖｅは、ｙ_ｉ（ｉ＝１、２、・・・、ｎ）の平均値である。

なお、第１特徴量が中心周波数の特定範囲において大きく変化しやすことが予めわかっている場合には、中心周波数が当該特定範囲内の第１特徴量の値を用いて、相関係数が算出されてもよい。

特徴量算出部１４２は、１つの中心周波数波形に対して、１種類の第２特徴量を算出してもよいし、複数種類の第２特徴量を算出してもよい。たとえば、特徴量算出部１４２は、中心周波数波形に対して、実効値、最大値、波高率、尖度、歪度、変動係数および寄与率の７種類の第２特徴量を算出してもよい。

［学習部の処理］
図１０は、図２の学習部１４３が行なう処理を説明するためのフローチャートである。学習部１４３は、まずステップＳ２１において、被試験装置１０が正常である第１正常期間に測定された複数の正常データからそれぞれ算出された第２特徴量を取得する。

次にステップＳ２２において、学習部１４３は、ステップＳ２１で取得した第２特徴量を学習データとして機械学習を行ない、異常の有無を判定するための異常度を算出するための演算式を決定する。具体的には、学習部１４３は、既知の機械学習手法（One Class Support Vector Machine：ＯＣ−ＳＶＭ）を用いて、学習データから異常と正常との分類境界を決定する。そして、学習部１４３は、決定した分類境界からの距離を異常度とし、当該異常度を算出するための演算式を決定する。

上述したように、特徴量算出部１４２は、１つの測定データに対してｍ種類（ｍは１以上の整数）の第１特徴量を算出可能であり、かつ、１種類の第１特徴量に対してｎ種類（ｎは１以上の整数）の第２特徴量を算出可能である。そのため、１つの測定データに対して特徴量算出部１４２から出力される第２特徴量の種類は、ｍ×ｎ種類となる。学習部１４３は、１つの正常データに対して算出されるｍ×ｎ種類の第２特徴量の各々の成分とする特徴量ベクトルを生成する。学習部１４３は、ＯＣ−ＳＶＭを用いて、正常データの個数と同じ個数の特徴量ベクトルから分類境界を決定し、決定した分類境界からの距離である異常度を算出するための演算式を決定する。

次にステップＳ２３において、学習部１４３は、被試験装置１０が正常である第２正常期間に測定された複数の正常データからそれぞれ算出された第２特徴量を取得する。第２正常期間は、第１正常期間の後の期間である。ステップＳ２４において、学習部１４３は、ステップＳ２２で決定された演算式に従って、正常データごとに、ステップＳ２３で取得した第２特徴量を用いて異常度を算出する。ステップＳ２５において、学習部１４３は、正常データごとに算出された異常度のばらつきに基づいて、異常判定しきい値を決定する。たとえば、学習部１４３は、異常度の平均値に異常度の標準偏差の定数倍（たとえば３倍）の値を加算した値を異常判定しきい値として決定する。ステップＳ２６において、学習部１４３は、ステップＳ２２で決定した異常度の演算式と、ステップＳ２５で決定した異常判定しきい値とを、被試験装置１０が異常か否かを判定するためのアルゴリズムとしてしきい値記憶部１４４に記憶させる。

なお、ステップＳ２３を省略し、学習部１４３は、ステップＳ２４において、ステップＳ２２で決定された演算式に従って、正常データごとに、ステップＳ２１で取得した第２特徴量を用いて異常度を算出してもよい。

［判定部の処理］
図１１は、図２の判定部１４６が行なう処理を説明するためのフローチャートである。判定部１４６は、まずステップＳ３１において、被試験装置１０の診断時に測定された測定データ（テストデータ）から算出された第２特徴量を取得する。ステップＳ３２において、判定部１４６は、しきい値記憶部１４４が記憶する演算式に従って、ステップＳ３１で取得した第２特徴量を用いて異常度を算出する。ステップＳ３３において、判定部１４６は、ステップＳ３２で算出した異常度と、しきい値記憶部１４４が記憶する異常判定しきい値とを比較し、比較結果に基づいて、被試験装置１０の異常の有無を判定する。

［本実施の形態の効果の検証実験１］
上記の状態監視方法によって算出される異常度が公知の方法によって算出される異常度に比べて優れていることの検証実験を行なった。

被試験装置として損傷のない新品の軸受を用い、異常が生じて使用不可能となるまで当該軸受を連続運転させた。軸受に設置した振動センサによって測定された振動加速度の測定データを２時間ごとに取得した。運転条件および測定条件は以下の通りである。

＜運転条件＞
軸受：アンギュラ玉軸受（型番7216：内径80mm、外径140mm、幅26mm）
ラジアル負荷：1.3kN
アキシアル負荷：1.3kN
回転速度：1500回転／分
潤滑方式：グリース
＜測定条件＞
測定データ：振動加速度
データ長さ：20秒
サンプリング速度：50kHz。

このような運転条件および測定条件により取得された測定データ群に対して、実施例のデータ演算方法に従って算出された異常度の変化と、比較例のデータ演算方法に従って算出された異常度の変化とを比較した。

実施例におけるデータ演算方法を以下のようにした。
・フィルタ処理部１４１による処理方法：各バンドパスフィルタのバンド幅１ｋＨｚ、複数のバンドパスフィルタの中心周波数の間隔１ｋＨｚ、周波数レンジ０〜２０ｋＨｚとし、１つの測定データから２０個の時間波形を生成する。
・特徴量算出部１４２が１つの測定データに対して算出する第１特徴量：６種類（時間波形、周波数波形およびケフレンシ波形の各々の実効値および尖度）。
・特徴量算出部１４２が１つの中心周波数波形に対して算出する第２特徴量：７種類（実効値、最大値、波高率、尖度、歪度、変動係数、寄与率）。寄与率は、中心周波数が１．５〜９．５ｋＨｚの範囲における、中心周波数と第１特徴量との相関係数の二乗である。
・学習部１４３の機械学習の手法：ＯＳ−ＳＶＭ。

実施例のデータ演算方法では、上記の測定データ群のうち１番目から所定番目までの初期の測定データ（正常データ）の各々から算出された６×７＝４２種類の第２特徴量を成分とする特徴量ベクトルを学習データとして機械学習を行なうことにより、異常度の演算式を決定した。

一方、比較例のデータ演算方法を以下のようにした。
・使用する特徴量：測定データで示される時間波形と、当該時間波形をフーリエ変換して得られる周波数波形と、当該周波数波形を時間波形と見なしてフーリエ変換して得られるケフレンシ波形との各々の実効値、最大値、波高率、尖度、歪度の計１５種類。
・機械学習の手法：ＯＳ−ＳＶＭ。

比較例のデータ演算方法は、上記の測定データ群のうち１番目から所定番目までの初期の測定データ（正常データ）の各々から算出された上記１５種類の特徴量を成分とする特徴量ベクトルを学習データとして機械学習を行なうことにより、異常度の演算式を決定した。

図１２は、実施例と比較例との異常度の変化を示す図である。図１２に示されるように、同一の測定データ群から算出した異常度であるにもかかわらず、実施例の異常度が上昇し始める時点は、比較例の異常度が上昇し始める時点に比べて約３ケ月早い。つまり、実施例によって算出された異常度は、軸受のわずかな異常を反映した値であることがわかる。このように、実施例によって算出された異常度を用いることにより、軸受の異常の有無を精度良く判定できることが確認された。

［本実施の形態の効果の検証実験２］
次に、被試験装置の異常に対する上記の状態監視方法で用いた第２特徴量の感度が優れていることの検証実験を行なった。

被試験装置として、損傷のない軸受と、人工の損傷を設けた軸受とを準備した。人工の損傷は、軸受の外輪軌道に放電加工で設けた微細な円筒穴である。これらの軸受を、ラジアル負荷およびアキシアル負荷のかかる状態で、一定の回転速度で運転した時の振動加速度を測定した。回転速度は３条件とした。放電穴直径（以下、損傷サイズ）は以下に示す５種類（損傷サイズ0.00mmは損傷のない軸受を示す）である。各損傷サイズで３６回振動加速度を測定した。運転条件および測定条件は次の通りである。

＜運転条件＞
軸受：アンギュラ玉軸受（型番7216：内径80mm、外径140mm、幅26mm）
ラジアル負荷：1.3kN
アキシアル負荷：1.3kN
回転速度：1000回転／分、1500回転／分、2000回転／分
潤滑方式：循環給油
損傷サイズ：0.00mm（損傷なし）、φ0.34mm、φ0.68mm、φ1.02mm、φ1.35mm
＜測定条件＞
測定データ：振動加速度
データ長さ：20秒
サンプリング速度：50kHz
測定回数：36回。

実施例におけるデータ演算方法を以下のようにした。
・フィルタ処理部１４１による処理方法：各バンドパスフィルタのバンド幅１ｋＨｚ、複数のバンドパスフィルタの中心周波数の間隔１ｋＨｚ、周波数レンジ０〜２０ｋＨｚとし、１つの測定データから２０個の時間波形を生成する。
・特徴量算出部１４２が１つの測定データに対して算出する第１特徴量：６種類（時間波形、周波数波形およびケフレンシ波形の各々の実効値および尖度）。
・特徴量算出部１４２が１つの中心周波数波形に対して算出する第２特徴量：７種類（実効値、最大値、波高率、尖度、歪度、変動係数、寄与率）。寄与率は、中心周波数が１．５〜９．５ｋＨｚの範囲における、中心周波数と第１特徴量との相関係数の二乗とした。

一方、比較例のデータ演算方法を以下のようにした。
・使用する特徴量（以下、「比較用特徴量」という）：測定データで示される時間波形と、当該時間波形をフーリエ変換して得られる周波数波形と、当該周波数波形を時間波形と見なしてフーリエ変換して得られるケフレンシ波形との各々の実効値、最大値、波高率、尖度、歪度の計１５種類。

損傷なしの軸受に対応する測定データから算出した特徴量（実施例では第２特徴量、比較例では比較用特徴量）と、損傷ありの軸受に対応する測定データから算出した特徴量（実施例では第２特徴量、比較例では比較用特徴量）との差をｔ検定を用いて評価した。

図１３は、実施例のデータ演算方法により算出された第２特徴量のｔ検定結果を示す図である。図１４は、比較例のデータ演算方法により算出された比較用特徴量のｔ検討結果を示す図である。図１３および図１４において、検定結果は常用対数化した値の絶対値である。ここでは、１０以上の数値である場合に、損傷なしの軸受に対応する測定データから算出した特徴量と、損傷ありの軸受に対応する測定データから算出した特徴量とに明らかな差があると言える。すなわち、１０以上の数値を示す特徴量は、軸受の異常に対して良好な感度を有する。そのため、１０以上の数値を示す特徴量を用いることにより、軸受の異常の有無の判定精度が向上する。

図１３および図１４に示されるように、実施例のｔ検定結果は、全般的に、比較例のｔ検定結果よりも優れている。

たとえば、回転速度1000回／分、損傷サイズφ1.35ｍｍの条件において、ｔ検定結果が３０以上である特徴量は、比較例のデータ演算方法では１種類であるのに対し、実施のデータ演算方法では７種類であった。具体的には、周波数波形の尖度を第１特徴量とする中心周波数波形の「波高率」、「尖度」、「歪度」および「変動係数」と、時間波形の実効値を第１特徴量とする中心周波数波形の「寄与率」と、周波数波形の実効値を第１特徴量とする中心周波数波形の「寄与率」と、ケフレンシ波形の実効値を第１特徴量とする中心周波数波形の「寄与率」とである。この中でも、周波数波形の尖度を第１特徴量とする中心周波数波形の「尖度」および「変動係数」と、周波数波形の実効値を第１特徴量とする中心周波数波形の「寄与率」と、ケフレンシ波形の実効値を第１特徴量とする中心周波数波形の「寄与率」とは、４０以上のｔ検討結果を示し、軸受の異常に対して優れた感度を有した。

また、比較例のデータ演算方法では、損傷サイズが最大のφ1.35mmの軸受であっても、３つの回転速度の全てについてｔ検討結果１０以上を示す比較用特徴量はなかった。これに対し、実施例のデータ演算方法では、損傷サイズが最大のφ1.35mmの軸受に対して、６種類の第２特徴量が、３つの回転速度の全てについてｔ検討結果１０以上を示した。具体的には、周波数波形の尖度を第１特徴量とする中心周波数波形の「実効値」、「最大値」および「変動係数」と、時間波形の実効値を第１特徴量とする中心周波数波形の「寄与率」と、周波数波形の実効値を第１特徴量とする中心周波数波形の「寄与率」と、ケフレンシ波形の実効値を第１特徴量とする中心周波数波形の「寄与率」とである。これらの第２特徴量は、回転速度にかかわらず、軸受の異常に対して優れた感度を有する。この中でも、周波数波形の実効値を第１特徴量とする中心周波数波形の「寄与率」は、損傷サイズがφ0.6mm以上の軸受に対しても、３つの回転速度の全てについてｔ検討結果１０以上を示した。

このように、実施例のデータ演算方法により算出される第２特徴量の中には、軸受の異常に対して優れた感度を示す特徴量が含まれることが確認された。

なお、４２種類の第２特徴量の中には、損傷サイズおよび回転速度がいずれであっても、１０未満のｔ検討結果を示すものが含まれる。しかしながら、上述したように、従来特徴量よりも優れたｔ検討結果を示す第２特徴量が複数存在している。そのため、４２種類の第２特徴量を成分とする特徴量ベクトルを用いて異常度を演算することにより、軸受の異常の有無を精度良く判定することが可能となる。

［本実施の形態の効果の検証実験３］
以下の（ａ）に記す点のみを変更し、上記の検証実験２と同様の検証実験３を行なった。
（ａ）特徴量算出部１４２が１つの中心周波数波形に対して算出する第２特徴量：１種類（寄与率）。寄与率は、バンドパスフィルタの中心周波数（中心周波数の範囲は限定しない）と第１特徴量との相関係数の二乗とした。

図１６は、検証実験３における第２特徴量のｔ検定結果を示す図である。図１６に示す検定結果は、図１３と同様に、常用対数化した値の絶対値である。図１６に示されるように、検証実験３のｔ検定結果も、全般的に、比較例のｔ検定結果（図１４参照）よりも優れている。すなわち、寄与率を求める際にバンドパスフィルタの中心周波数の範囲を限定しなくても、寄与率は、軸受の異常に対して優れた感度を有した。

［本実施の形態の効果の検証実験４］
以下の（ｂ）（ｃ）に記す点のみを変更し、上記の検証実験２と同様の検証実験４を行なった。
（ｂ）特徴量算出部１４２が１つの測定データに対して算出する第１特徴量：６種類（時間波形、周波数波形およびケフレンシ波形の各々の実効値および尖度の常用対数）。なお、時間波形の尖度については負の値をとる場合があるため、全データの最小値を各データから差し引く処理を行なうことにより得られる値の常用対数を用いた。
（ｃ）特徴量算出部１４２が１つの中心周波数波形に対して算出する第２特徴量：１種類（寄与率）。寄与率は、バンドパスフィルタの中心周波数（中心周波数の範囲は限定しない）と第１特徴量との相関係数の二乗とした。

図１７は、検証実験４における第２特徴量のｔ検定結果を示す図である。図１７に示す検定結果は、図１３と同様に、常用対数化した値の絶対値である。図１７に示されるように、検証実験４のｔ検定結果も、全般的に、比較例のｔ検定結果（図１４参照）よりも優れている。すなわち、第１特徴量としてパラメータの値の常用対数値を用いても、寄与率は、軸受の異常に対して優れた感度を有した。

［変形例］
上記の説明では、特徴量算出部１４２は、バンドパスフィルタを用いて生成された時間波形の特徴を示す時間領域特徴量と、当該時間波形をフーリエ変換して得られる周波数波形の特徴を示す周波数領域特徴量と、当該周波数波形をフーリエ変換して得られるケフレンシ波形の特徴を示すケフレンシ領域特徴量とのいずれかを第１特徴量として算出した。しかしながら、ステップＳ１２（図６参照）において、特徴量算出部１４２は、上記の時間領域特徴量、周波数領域特徴量およびケフレンシ領域特徴量の少なくとも２つの特徴量を算出する工程と、当該少なくとも２つの特徴量から第１特徴量を算出する工程とを行なってもよい。たとえば、特徴量算出部１４２は、上記の時間領域特徴量（第３特徴量）、周波数領域特徴量（第３特徴量）およびケフレンシ領域特徴量（第３特徴量）の少なくとも２つの相関関係を示す値（和、差、比、平均）などを第１特徴量として算出する。

図１５は、変形例の第１特徴量の一例を示す図である。図１５において、（ａ）は、周波数領域特徴量の１つである、周波数波形の最大値ｐを示し、（ｂ）は、ケフレンシ領域特徴量の１つである、ケフレンシ波形の最大値ｑを示し、（ｃ）は、第１特徴量を示す。図１５に示されるように、特徴量算出部１４２は、複数のバンドパスフィルタの各々について、当該バンドパスフィルタを用いて生成された時間波形から得られた周波数波形の最大値ｐと、当該周波数波形から得られたケフレンシ波形の最大値ｑとの比（ｑ／ｐ）を第１特徴量として算出する。

上記の説明では、学習部１４３は、ＯＣ−ＶＳＭを用いて機械学習を行なった。しかしながら、学習部１４３は、ＳＶＭ以外にも、ランダムフォレスト、ロジスティック回帰、決定木、ニューラルネットワークを用いて機械学習を行なってもよい。

図１３に示されるように、ｔ検討結果が優れた特定の第２特徴量（図１３では、第１特徴量を「周波数波形の実効値」とするときの寄与率）が予め分かっている場合には、特徴量算出部１４２は、当該特定の第２特徴量のみを算出してもよい。この場合、学習部１４３は、異常判定しきい値のみを決定する。判定部１４６は、テストデータに対して算出された第２特徴量と異常判定しきい値との比較結果に基づいて、被試験装置１０の異常の有無を判定すればよい。

［作用・効果］
本実施の形態に係る状態監視方法は、被試験装置（被試験対象物）１０に設置した振動センサ２０から取得した測定データを用いて被試験装置１０の状態を監視する方法であり、第１〜第４工程を備える。第１工程（Ｓ１１）では、測定データに対して、互いに異なる中心周波数を有する複数のバンドパスフィルタの各々を用いたフィルタ処理を行なうことにより、複数のバンドパスフィルタにそれぞれ対応する複数の時間波形（第１波形）を生成する。第２工程（Ｓ１２）では、複数の時間波形の各々に対して第１特徴量を算出する。第３工程（Ｓ１３，Ｓ１４）では、中心周波数を横軸とし、第１特徴量を縦軸とするグラフに、複数のバンドパスフィルタの各々について、当該バンドパスフィルタの中心周波数と、当該バンドパスフィルタを用いて生成された時間波形に対して算出された第１特徴量とをプロットすることにより得られる中心周波数波形（第２波形）の特徴を示す第２特徴量を算出する。第４工程（Ｓ３２、Ｓ３３）では、測定データのうち被試験装置１０の診断時にセンサから得られたテストデータに対して第１〜第３工程を行なうことにより算出された第２特徴量に基づいて、被試験装置１０が異常か否かを判定する。

上記の構成によれば、被試験装置１０の異常に対する感度が良好な第２特徴量に基づいて被試験装置１０が異常か否かを判定する。これにより、被試験装置１０の異常の有無の判定精度を高めることができる。その結果、被試験装置１０の異常を早期に発見することができる。

被試験装置１０に異常が生じた場合、測定データの特定の周波数帯域に変化が見られることが多い。そのため、当該特定の周波数帯域のみを抽出することにより、異常の有無を判定しやすくなる。そこで、従来の状態監視方法では、異常の有無の判定の精度を上げるために、異常の有無の判定に有効な周波数帯域を予め最適化していた。ただし、異常の内容によって、変化が見られる周波数帯域が異なるため、被試験装置１０ごとに周波数帯域の最適化を行なう必要があった。

上記の構成によれば、第２特徴量は、中心周波数波形の特徴を示す。中心周波数波形は、中心周波数を横軸とし、第１特徴量を縦軸とするグラフに、複数のバンドパスフィルタの各々の中心周波数と対応する第１特徴量とをプロットすることにより得られる。被試験装置１０の異常によって変化が生じる特定の周波数帯域が不明であったとしても、中心周波数波形における当該特定の周波数帯域に含まれる中心周波数の箇所に、異常の程度に応じた変化が見られる。そのため、当該箇所の変化が反映された第２特徴量を確認することにより、被試験装置１０の異常の有無を容易に判定することができる。

状態監視方法は、測定データのうち被試験装置１０の正常時に振動センサ２０から得られた正常データに対してステップＳ１１〜Ｓ１４を行なうことにより得られた第２特徴量を学習データとして機械学習を行なうことにより、被試験装置１０が異常か否かを判定するためのアルゴリズムを決定する第５工程（Ｓ２２，Ｓ２４，Ｓ２５）をさらに備える。ステップＳ３２、Ｓ３３では、判定部１４６は、アルゴリズムに従って被試験装置１０が異常か否かを判定する。

これにより、第２特徴量の種類が複数であったとしても、これらを総合的に考慮したアルゴリズムに従って異常か否かを判定することができる。

第１特徴量は、たとえば、対応する時間波形と、当該時間波形をフーリエ変換して得られる周波数波形（第３波形）と、周波数波形を時間波形と見なしてフーリエ変換して得られるケフレンシ波形（第４波形）とのいずれか１つの波形の特徴を示す。

あるいは、ステップＳ１２は、複数の時間波形の中の第１特徴量に対応する時間波形と、当該対応する時間波形をフーリエ変換して得られる周波数波形と、周波数波形を時間波形と見なしてフーリエ変換して得られるケフレンシ波形との少なくとも２つの波形の各々に対して、当該波形の特徴を示す特徴量（第３特徴量）を算出する工程と、当該少なくとも２つの波形の各々に対して算出した特徴量から第１特徴量を算出する工程とを含んでもよい。

第２特徴量は、たとえば、前記第２波形の実効値、最大値、波高率、尖度、歪度、変動係数および寄与率のいずれかである。

今回開示された実施の形態は、すべての点で例示であって制限的なものではないと考えられるべきである。本発明の範囲は、上記した実施の形態の説明でなくて特許請求の範囲によって示され、特許請求の範囲と均等の意味および範囲内でのすべての変更が含まれることが意図される。

１０ 被試験装置、２０ 振動センサ、１００ 状態監視装置、１１０ Ａ／Ｄコンバータ、１２０ データ取得部、１３０ 記憶装置、１４０ データ演算部、１４１ フィルタ処理部、１４２ 特徴量算出部、１４３ 学習部、１４４ しきい値記憶部、１４６ 判定部、１５０ 表示部。

Claims (6)

1. 被試験対象物に設置したセンサから取得した測定データを用いて前記被試験対象物の状態を監視する状態監視方法であって、
   前記測定データに対して、互いに異なる中心周波数を有する複数のバンドパスフィルタの各々を用いたフィルタ処理を行なうことにより、前記複数のバンドパスフィルタにそれぞれ対応する複数の第１波形を生成する第１工程と、
   前記複数の第１波形の各々に対して第１特徴量を算出する第２工程と、
   中心周波数を横軸とし、第１特徴量を縦軸とするグラフに、前記複数のバンドパスフィルタの各々について、当該バンドパスフィルタの前記中心周波数と、当該バンドパスフィルタを用いて生成された時間波形に対して算出された前記第１特徴量とをプロットすることにより得られる第２波形の特徴を示す第２特徴量を算出する第３工程と、
   前記測定データのうち前記被試験対象物の診断時に前記センサから得られたテストデータに対して前記第１〜第３工程を行なうことにより算出された前記第２特徴量に基づいて、前記被試験対象物が異常か否かを判定する第４工程とを備える状態監視方法。
2. 前記複数の第１波形は時間波形であり、
   前記第１特徴量は、前記複数の第１波形の中の当該第１特徴量に対応する第１波形と、前記対応する第１波形をフーリエ変換して得られる第３波形と、前記第３波形をフーリエ変換して得られる第４波形とのいずれか１つの波形の特徴を示す、請求項１に記載の状態監視方法。
3. 前記複数の第１波形は時間波形であり、
   前記第２工程は、
   前記複数の第１波形の中の前記第１特徴量に対応する第１波形と、前記対応する第１波形をフーリエ変換して得られる第３波形と、前記第３波形をフーリエ変換して得られる第４波形との少なくとも２つの波形の各々に対して、当該波形の特徴を示す第３特徴量を算出する工程と、
   前記少なくとも２つの波形の各々に対して算出した前記第３特徴量から前記第１特徴量を算出する工程とを含む、請求項１に記載の状態監視方法。
4. 前記測定データのうち前記被試験対象物の正常時に前記センサから得られた正常データに対して前記第１〜第３工程を行なうことにより得られた前記第２特徴量を学習データとして機械学習を行なうことにより、前記被試験対象物が異常か否かを判定するためのアルゴリズムを決定する第５工程をさらに備え、
   前記第４工程では、前記アルゴリズムに従って前記被試験対象物が異常か否かを判定する、請求項１から３のいずれか１項に記載の状態監視方法。
5. 前記第２特徴量は、前記第２波形の実効値、最大値、波高率、尖度、歪度、変動係数および寄与率のいずれかである、請求項１から４のいずれか１項に記載の状態監視方法。
6. 請求項１から５のいずれか１項に記載の方法を用いて、前記被試験対象物を診断する、状態監視装置。