一．技術分野
本発明は設備診断、医療診断などにおいて、対象物の状態変化の有無を判定するための状態判定装置及びオンライン状態監視・診断システムに関するものである。
二．背景技術
（１）従来、状態監視の対象物から測定した波形データの確率密度関数は、正常状態時に「正規分布」に従うと仮定して状態判定を行う（［１］特許公開２０００−１７１２９１、［２］豊田利夫、二保知也：「正常時の振動波形のみを用いた回転機械の異常診断」、日本設備管理学会誌、Ｖｏｌ．１１，Ｎｏ．１，１９９９年、ｐ．４−１１．［３］劉 信芳、豊田利夫、陳 鵬、馮 芳、二保知也：「Ｉｎｆｏｒｍａｔｉｏｎ Ｄｉｖｅｒｇｅｎｃｅによる回転機械の異常診断」、精密工学会誌、Ｖｏｌ．６６、Ｎｏ．１、２０００年、ｐ．１５７−１６２．）
（２）対象物の状態変化の有無を判定するために用いる特徴パラメータを正規分布の確率変数と仮定して統計検定を行う（［４］河部佳樹、豊田利夫、江口透、福田和久：「無次元兆候パラメータを用いた回転機械の劣化傾向管理（ＩＩＩ）」、平成７年度日本設備管理学会秋季研究発表大会論文集、Ｖｏｌ．２、１９９５年、ｐ．３２．）
状態監視の対象物から測定した信号が脈動信号の場合は、特徴パラメータや確率密度分布関数などにより状態判定が困難である。従来、脈動信号の中に含まれる特異成分（あるいは、異常成分）を検出する方法として、
（３）時間−周波数解析手法（ウェーブレット解析、ウィグナ分布解析、短時間ＦＦＴなど）を用いた局所異常の検出法がある（［５］章 忠、中堀智之、川畑洋昭：「高速ウェーブレット変換およびその脳波解析への応用」、日本機械学会論文集（Ｃ編）、Ｖｏｌ．６５，Ｎｏ．６３２，１９９９年、ｐ．１９１５−１９２１．）
（４）脈動信号のピーク値などに合わせた所定の間隔でサンプリングした値は所定の値以下、または、前回サンプリングした値との差が所定の値以上の場合に異常と判定する。（［６］特開平６−３１７２１５）
（５）脈動信号のピーク値などに合わせた所定の間隔でサンプリング毎の代表値のばらつき（標準偏差）が所定のしきい値以上の場合に異常と判定する。（［７］特開平７−２９３３１１）
三．発明が解決しようとする課題
しかしながら、前記従来の諸手法には次のような問題があった。
第１と第２の方法では、波形データと特徴パラメータの値が正規分布に従うと仮定しているが、実際に正常状態でも必ずしも正規分布に従うとは限らないから、波形データと特徴パラメータの値が正規分布に従うと仮定して対象物の状態変化を統計検定により判定すると、誤判定をもたらす。
第３の方法では、処理に必要な時間が長いため、アルタイムの特異点の検出が困難である。
第４と第５の方法では、波形データを所定の間隔でサンプリングする場合、回転数の変化や負荷の変化などによるピーク位置が変化した時、サンプリング値は正常でも大きく変化するから、誤判定をもたらす。また、サンプリング毎の代表値のばらつき（標準偏差）を計算すると、リアルタイムの特異検出が困難である。
四．課題を解決するための手段
上記に述べたような問題点を解決するために、本発明においては、測定した波形データを既知確率分布（例えば、正規分布）の波形データに変換した後、あるいは、波形データから算出された特徴パラメータを既知確率分布（例えば、正規分布）の特徴パラメータに変換した後、統計検定や可能性理論や情報理論などにより対象物の状態判定を行う。
また、測定した対象物の信号が脈動信号である場合は、特徴パラメータやスペクトルなどによる特異成分の抽出及び状態判定が困難であるので、雑音を除去した信号の包絡線波形データを求め、正規化された包絡線波形データおよび周期パルスにより信号中の特異成分を検出して状態判定を行う。または、上記に雑音を除去した脈動信号のピーク波形データを求め、ピーク波形データの正規化処理を行い、正規化されたピーク波形データおよび周期パルス波形データにより脈動信号中の特異成分を検出して状態判定を行う。
五．発明の効果
対象物の状態監視のために測定された波形データや特徴パラメータは必ずしも正規分布に従うと限らないため、これらの波形データあるいは特徴パラメータを正規分布に従うと仮定して対象物の状態変化の有無を統計検定により判定し、或いは、状態予測を行うと、大きな誤差をもたらす。本発明においては、測定した波形データを既知確率分布（例えば、正規分布）の波形データに変換した後、あるいは、波形データから算出された特徴パラメータを既知確率分布（例えば、正規分布）の特徴パラメータに変換した後、統計検定や可能性理論や情報理論などにより対象物の状態判定を行う。従って、本発明の手法は状態判定、或いは、状態予測の精度が従来正規分布に従うと仮定した場合より高い。
脈動信号の場合、測定した脈動信号の波形データには一部のピークが欠落し、または波形データの局所に微小な特異成分が存在する時、従来の判定法（特徴パラメータやスペクトルや確率密度関数など）による検出が困難である。本発明においては、雑音を除去した脈動信号の包絡線波形データを求め、正規化された包絡線波形データおよび周期パルスにより信号中の特異成分を検出して状態判定を行う。または、上記に雑音を除去した脈動信号のピーク波形データを求め、ピーク波形データの正規化処理を行い、正規化されたピーク波形データおよび周期パルス波形データにより脈動信号中の特異成分を検出して状態判定を行う。本発明の方法では、主な信号処理はハードウェアでも実現でき、数値処理装置（計算機）にかかる負担が小さいから、リアルタイムに異常の検出が実現できる。また、正常状態の波形データを基準としないため、脈動周期が変化しても検出と判定の結果に与える影響が少ない。
六．発明を実施するための最良の形態
１．特徴パラメータについて
状態判定のために用いる特徴パラメータは時間領域と周波数領域の特徴パラメータがある。周波数領域の特徴パラメータは（参考文献１）で定義されているものがある。（［８］陳 鵬，豊田利夫：遺伝的プログラミングによる周波数領域の特徴パラメータの自己再組織化，日本機械論文集（Ｃ編），Ｖｏｌ．６５ Ｎｏ．６３３，ｐｐ．１９４６−１９５３，１９９８．）ここで時間領域の特徴パラメータについて詳述する。
対象物の状態変化の有無を判定するために用いる時間領域の特徴パラメータは次の通りである。
１）無次元特徴パラメータ
測定した時系列波形データから、フィルタを用いて、低、中、高周波数領域の波形データを抽出する。抽出した波形データｘ（ｔ）を次式で正規化する。
ここで、ｘ’ｉはＡ／Ｄ変換後のｘ（ｔ）の離散値であり、μとＳはそれぞれｘ’ｉの平均値と標準偏差である。
従来用いられる無次元特徴パラメータは式（２）〜式（１３）に示す（［９］Ｐｅｎｇ ＣＨＥＮ，Ｔｏｓｈｉｏ ＴＯＹＯＴＡ，Ｙｕｅｔｏｎ ＬＩＮ，Ｆｅｉｙｕｅ ＷＡＮＧ：ＦＡＩＬＵＲＥ ＤＩＡＧＮＯＳＩＳ ＯＦ ＭＡＣＨＩＮＥＲＹ ＢＹ ＳＥＬＦ−ＲＥＯＲＧＡＮＩＺＡＴＩＯＮ ＯＦ ＳＹＭＰＴＯＭ ＰＡＲＡＭＥＴＥＲＳ ＩＮ ＴＩＭＥ ＤＯＭＡＩＮ ＵＳＩＮＧ ＧＥＮＥＴＩＣ ＡＬＧＯＲＩＴＨＭＳ，Ｉｎｔｅｒｎａｔｉｏｎａｌ Ｊｏｕｒｎａｌ ｏｆ Ｉｎｔｅｌｌｉｇｅｎｔ Ｃｏｎｔｒｏｌ ａｎｄ Ｓｙｓｔｅｍ，Ｖｏｌ．３，Ｎｏ．４，ｐｐ．５７１−５８５，１９９９．）
ここで、
は絶対平均値で，Ｎはデータの総数である。
は標準偏差である。
ここで、μｐは波形の極大値（ピーク値）の平均値である。
ここで、｜μｍａｘ｜は波形の１０個の最大値の平均値である。
ここで、σｐは極大値の標準偏差値である。
ここで、μＬとσＬはそれぞれ極小値（谷値）の平均値と標準偏差値である。
式（２）〜式（１３）は従来の特徴パラメータであるが、数値計算で容易に高速計算するために、「区間特徴パラメータ」は式（１４）〜式（２１）のように新たに提案する。
但し、ｘｉ≧ｋσ、ｋは任意に設定できるが、例えば、ｋ＝０．５、１、２。μｋ１はｘｉの平均値である。ｔは任意に設定できるが、例えば、ｔ＝０．５、１、２、３、４。
但し、ｘｈ≦−ｈσ、ｈは任意に設定できるが、例えば、ｋ＝０．５、１、２。μｈ１はｘｈの平均値である。ｔは任意に設定できるが、例えば、ｔ＝０．５、１、２、３、４。
但し、ｈ０は時系列波形が単位時間あたり０レベルをクロースする頻度、ｈｐは単位時間あたり時系列波形のピークの数である。
但し、ｈｎσは波形が単位時間あたりｎσレベルをクロースする頻度である。ｎは任意に設定できるが、例えば、ｎ＝０．５、１、２。
但し、ｈ−ｎσは波形が単位時間あたり−ｎσレベルをクロースする頻度。である。ｎは任意に設定できるが、例えば、ｎ＝０．５、１、２。
２）有次元特徴パラメータ
有次元特徴パラメータを計算するとき、測定した波形データに対して、式（１）のような正規化をしない。
ここで、｜ｘｉ｜ｐは波形データの絶対値のピーク値（極大値）、Ｎｐはピーク値の総数である。
なお、上記の特徴パラメータ以外にも多くの特徴パラメータが定義できるが、本方法を応用するとき、まず上記の特徴パラメータで試し、もし状態識別の効果が良くなければ、更に他の特徴パラメータを追加定義すればよい。
２．波形データと特徴パラメータを指定の確率分布への変換
測定した波形データをｘ＊ ｉで、波形データから算出した特徴パラメータをｐ＊ ｉで表す。ｘ＊ ｉあるいはｐ＊ ｉを用いて、統計理論により状態判定や状態予測を行う場合、ｘ＊ ｉあるいはｐ＊ ｉがどのような確率分布に従うか事前に知る必要がある。しかし、ｘ＊ ｉあるいはｐ＊ ｉがどのような確率分布従うかは事前に知らない場合が多い。そこで、指定の既知確率分布関数をΣとすると、次式を用いてｘ＊ ｉあるいはｐ＊ ｉをΣに従う確率変数ｘｉあるいはｐｉに変換することができる。
ここで、Ｆｘｉ（ｘ＊ ｉ）とＦｐｉ（ｐ＊ ｉ）はそれぞれｘ＊ ｉとｐ＊ ｉとの累積確率分布（或いは、累積度数分布）であり、Σ−１はΣの逆関数である。例えば、Σは正規分布、ワイブル分布、指数分布、ガンマ分布などである。ｘ＊ ｉあるいはｐ＊ ｉをｘｉあるいはｐｉに変換した後、統計検定などを用いて状態判定や状態予測を行う。
なお、元の波形データｘ＊ ｉは図１のように４種類に分けられる。つまり、平均値より大きいデータｘｉ＋、平均値より小さいデータｘｉ−、式（１）による正規化後の絶対値データ｜ｘｉ｜及び全体波形データｘ＊ ｉＡであり、それぞれの累積確率分布（或いは、累積度数分布）はＦｘｉ＋（ｘｉ＋）、Ｆｘｉ−（ｘｉ−）、Ｆ｜ｘｉ｜（｜ｘｉ｜）及びＦｘｉＡ（ｘ＊ ｉＡ）で表すが、以下、特に指定しなければ、統一にＦｘｉ（ｘ＊ ｉ）で表す。また、特徴パラメータｐ＊ ｉは正規化後のｘｉ＋とｘｉ−、｜ｘｉ｜及びｘ＊ ｉＡのいずれかで計算される。
ここで、例として、Σを正規分布に指定した場合、ｘ＊ ｉあるいはｐ＊ ｉを正規分布に従うように変換する方法について詳述する。
正規分布の確率密度関数ｆ（ｔ）は次式で表す。
ここで、μは確率変数ｔの平均値で、σは標準偏差である。
（１）基準状態を参照する場合
測定対象の基準状態、例えば、１回目測定した時の状態、を決めて、この時の波形データｘ＊ ｉ０あるいは特徴パラメータｐ＊ ｉ０の確率密度関数と累積確率分布関数をそれぞれｆｘ０（ｘ＊ ｉ０）とＦｘ０（ｘ＊ ｉ０）あるいはｆｐ０（ｐ＊ ｉ０）とＦｐ０（ｐ＊ ｉ０）とする。平均値μが０で標準偏差σが１である標準正規分布の確率密度関数をφ（ｘｉ）、標準正規分布の確率分布関数をΦ（ｘｉ）とする。なお、離散データであるｘ＊ ｉ０あるいはｐ＊ ｉ０の確率密度関数の替わりに「頻度分布関数」や「ヒスとグラム」を用いてもよいが、以下には「確率密度関数」を用いて説明する。
１）正規分布に従う平均値に変換する方法
状態ｋでのｘ＊ ｉｋとｐ＊ ｉｋを次式でそれぞれ正規分布に従う平均値μｘｉｋｏとμｐｉｋ０に変換される。なお、状態ｋは任意の状態であり、基準状態も含む。
ここで、Φ−１はΦの逆関数で、σｘｉ０とσｐｉ０はそれぞれ正規分布に変換されたｘ＊ ｉｋとｐ＊ ｉｋとの標準偏差であり、次式で求められる。
平均値μｘｉｋ０とμｐｉｋ０を用いて状態判定や状態予測を行う。
２）直接変換法
状態ｋでのｘ＊ ｉｋとｐ＊ ｉｋを次式で正規分布の確率変数に変換する。
ここで、ＳｘｋとＳｐｋはそれぞれｘ＊ ｉｋとｐ＊ ｉｋの標準偏差で、μｘｋとμｐｋはそれぞれｘ＊ ｉｋとｐ＊ ｉｋの平均値である。ｘ’ｉｋ０あるいはｐ’ｉｋ０を用いて状態判定や状態予測を行う。
（２）基準状態を参照しない場合
まず、状態ｋでの波形データｘ＊ ｉｋと特徴パラメータｐ＊ ｉｋとの確率密度関数（あるいは、頻度分布）をそれぞれｆｘｋ（ｘ＊ ｉｋ）とｆｐｋ（ｐ＊ ｉｋ）とし、確率分布関数（あるいは、累積度数分布）をＦｘｋ（ｘ＊ ｉｋ）とＦｐｋ（ｐ＊ ｉｋ）とする。
１）直接変数変換法
波形データｘ＊ ｉｋと特徴パラメータｐ＊ ｉｋを次式で正規分布の確率変数に変換する。
ここで、ＳｘｋとＳｐｋはそれぞれｘ＊ ｉｋとｐ＊ ｉｋとの標準偏差で、μｘｋとμｐｋはｘ＊ ｉｋとｐ＊ ｉｋとの平均値である。ｘ’ｉｋあるいはｐ’ｉｋを用いて状態判定や状態予測を行う。
２）正規分布に従う平均値に変換する方法
正規分布に従う平均値は次式で求められる。
ここで、σｘｉｋとσｐｉｋは次式で求められる。
μｘｉｋとμｐｉｋは正規分布に従うので、μｘｉｋあるいはμｐｉｋを用いて状態判定や状態予測を行う。
３）間接変換法
Ｍ個の特徴パラメータｐ＊ ｉｋの最小値と最大値をそれぞれ（ｐ＊ ｉｋ）ｍｉｎと（ｐ＊ ｉｋ）ｍａｘとする。（ｐ＊ ｉｋ）ｍｉｎから（ｐ＊ ｉｋ）ｍａｘまでＮ個の等間隔区間に分割する。各区間の代表値（例えば、中央値）をｐ＊ ｉｋｊとする。ここで、ｊ＝１〜Ｎ。ｐ＊ ｉｋｊをｐ＊ ｉｋの替わりに式（３６）または式（３８）に代入すると、Ｎ個のｐ’’ｉｋまたはμ’ｐｉｋが得られる。ｐ’’ｉｋまたはμ’ｐｉｋを用いて状態判定や状態予測を行う。
同様に、Ｍ個の波形データｘ＊ ｉｋの最小値と最大値をそれぞれ（ｘ＊ ｉｋ）ｍｉｎと（ｘ＊ ｉｋ）ｍａｘとする。（ｘ＊ ｉｋ）ｍｉｎから（ｘ＊ ｉｋ）ｍａｘまでＮ個の等間隔区間に分割する。各区間の代表値（例えば、中央値）をｘ＊ ｉｋｊとする。ここで、ｊ＝１〜Ｋ。ｘ＊ ｉｋｊをｘ＊ ｉｋの替わりに式（３５）または式（３４）に代入すると、Ｋ個のｘ’’ｉｋまたはμ’ｘｉｋが得られる。ｘ’’ｉｋまたはμ’ｘｉｋを用いて状態判定や状態予測を行う。
４）特徴パラメータｐ＊ ｉｋの区間平均値を求める方法
Ｎ個特徴パラメータｐ＊ ｉｋを求めた後、Ｍ組に分割して第ｊ組の平均値は次のように求める。
ここで、Ｎｊは第ｊ組にあるｐ＊ ｉｋの数である。μｋ （ｊ）は近似的に正規分布に従うので、μｋ （ｊ）を用いて状態判定や状態予測を行う。
波形データｘ＊ ｉｋを式（２９）、式（３３）、式（３５）、式（３７）で正規確率分布の波形データに変換したμｘｉｋ０、ｘ’ｉｋ０、ｘ’ｉｋ、μｘｉｋ、ｘ’’ｉｋ、μ’ｘｉｋは「正規分布の波形データ」と呼ぶ。また、特徴パラメータｐ＊ ｉを式（３０）、式（３４）、式（３６）、式（３８）で正規確率変数に変換したμｐｉｋ０、ｐ’ｉｋ０、ｐ’ｉｋ、μｐｉｋ、ｐ’’ｉｋ、μ’ｐｉｋ、μｋ （ｊ）は「正規分布の特徴パラメータ」と呼ぶ。
３．正規分布の特徴パラメータによる状態判別方法
ここで、正規分布の特徴パラメータを用いて、対象物の状態を判別する方法を述べる。
（１）統計理論による判別
１）正規分布の特徴パラメータの平均値の検定
状態ｋと状態ｙで求めた正規分布の特徴パラメータをそれぞれｐｉｋとｐｉｙとする。ここで、ｉ＝１〜Ｍ、Ｍは使用する正規分布の特徴パラメータの総数を表す。ｐｉｋとｐｉｙとの平均値をそれぞれμｉｋとμｉｙ、ｐｉｋとｐｉｙとの標準偏差をそれぞれＳｉｋとＳｉｙとする。一般にＪ個のｐｊの平均値μと標準偏差Ｓは次の式で計算する。
μｉｋとμｉｙが等しいか否かの検定は次のように行う（参考文献３参照）
（参考文献３）Ｋ．Ａ．Ｂｒｏｗｎｌｅｅ．Ｓｔａｔｉｓｔｉｃａｌ Ｔｈｅｏｒｙ ａｎｄ Ｍｅｔｈｏｄｏｌｏｇｙ ｉｎ Ｓｃｉｅｎｃｅ ａｎｄ Ｅｎｇｉｎｅｅｒｉｎｇ，Ｓｅｃｏｎｄ Ｅｄｉｔｉｏｎ，Ｔｈｅ Ｕｎｉｖｅｒｓｉｔｙ Ｃｈｉｃａｇｏ，１９６５
が成立すれば、有意水準αで「μｉｋとμｉｙとが等しくない」と判定する。ここで、ｔα／２（Ｊ−１）は自由度Ｊ−１のｔ分布の確率密度関数が下側確率α／２に対するパーセント点である。
２）正規分布の特徴パラメータの分散の検定
ＳｉｋとＳｉｙが同じか否かの検定は次のように行う（参考文献４参照）。
（参考文献４）Ｋ．Ａ．Ｂｒｏｗｎｌｅｅ．Ｓｔａｔｉｓｔｉｃａｌ Ｔｈｅｏｒｙ ａｎｄ Ｍｅｔｈｏｄｏｌｏｇｙ ｉｎ Ｓｃｉｅｎｃｅ ａｎｄ Ｅｎｇｉｎｅｅｒｉｎｇ，Ｓｅｃｏｎｄ Ｅｄｉｔｉｏｎ，Ｔｈｅ Ｕｎｉｖｅｒｓｉｔｙ Ｃｈｉｃａｇｏ，１９６５
が成立すれば、有意水準αで「ＳｉｋとＳｉｙとが等しくない」と判定する。ここで、Ｆα／２（Ｊ−１，Ｊ−１）は自由度Ｊ−１のＦ分布の確率密度関数が下側確率α／２に対するパーセント点である。
有意水準αを変えたとき、式（４４）または式（４５）を満足するか否かを確認することにより、状態ｙが状態ｋに対する状態変化の程度を決める。状態変化の程度を有意水準αにより決定する例は表１に示す。なお、設備診断の場合、状態ｋを正常状態とすれば、状態ｙは正常状態か、注意状態か、危険状態かの判別については、表１のように「正常」（α１）、「注意」（α２）、及び「危険」（α３）のように設定し検定することができる。つまり、式（４４）または式（４５）はα１の時に成立しなければ「正常」と判定する。また、式（４４）または式（４５）はα２の時に成立すれば「注意」と、α３の時に成立すれば「危険」と判定する。なお、表１中のαの数値範囲は例であり、設備の重要度などによって決められる。
数個の特徴パラメータを用いて状態判定を行う場合、判定の結果は最も状態変化が大きいと示した特徴パラメータの判定結果に従われる。例えば、３つの特徴パラメータｐ１、ｐ２、ｐ３を用いて判定するとき、ｐ１の判定結果が「注意」、ｐ２の判定結果が「正常」、ｐ３の判定結果が「危険」である場合、最終の判定結果は「危険」とする。
３）信頼区間による判定
基準時点で測定した波形データから求めた正規分布の特徴パラメータの平均値をμｉ０とし、他の時点で測定した波形データから求めた正規分布の特徴パラメータの平均値をμｉｋとすると、μｉ０の信頼区間は次式で与えられる。
ここで、ｔα／２（Ｊ−１）は自由度Ｊ−１のｔ分布の確率密度関数が下側確率α／２に対するパーセント点である。Ｓｉ０は波形データから求めたｐｉ０の標準偏差である。μｉｋは式（４６）に示す区間内にあれば、１−αの確率でμｉ０との差がないという。μｉ０の９９％信頼区間は、Ｊ＞１０のとき、近似的に次のようになる。
従って、μｉｋは式（４７）の範囲を超えれば、９９％の確率でμｉ０と違うと判定する。更に、測定した波形データで求めたμｉｋの信頼区間は、次式で求められる。
ここで、Ｓｉｋは波形データから求めたｐｉｋの標準偏差である。
表１のαを式（４６）に代入すれば、次のような信頼区間が得られる。
μｉｋがこれらの区間内にあるか否かによって状態の判定を行う。
数個の特徴パラメータを用いて状態判定を行う場合、最終の判定結果は式（４９），（５０），（５１）に記述される内容と同じである。
（２）可能性理論による判別
１）可能性分布関数の作成
状態ｋにおける波形データで正規分布の特徴パラメータｐｉの値を算出した後、ｐｉの確率密度関数ｆｋ（ｐｉ）から可能性分布関数Ｐｋ（ｐｉ）を式（５２）で求める。可能性理論によれば、ｐｉがどのような確率分布に従っても、その可能性分布関数が求められる。ｐｉが正規分布に従う場合、Ｎ段の可能性分布関数ｐｋ（ｐｉ）は次のように求める（参考文献５参照）。
ここで、
但し、上式中、ｐｉｘ＝ｍｉｎ｛ｐｉ｝＋ｘ×（ｍａｘ｛ｐｉ｝−ｍｉｎ｛ｐｉ｝）／Ｎ、ｘ＝１〜Ｎ，Ｓｉはｐｉの標準偏差，μｉはのｐｉ平均値である。
（参考文献５）Ｌ．Ｄａｖｉｓ：ＨＡＮＤＢＯＯＫ ＯＦ ＧＥＮＥＴＩＣ ＡＬＧＯＲＩＴＨＭＳ，Ｖａｎ Ｎｏｓｔｒａｎｄ Ｒｅｉｎｈｏｌｄ，Ａ Ｄｉｖｉｓｉｏｎ ｏｆ Ｗａｄｓｗｏｒｔｈ，Ｉｎｃ（１９９０）
２）可能性の求め方
図２のように、状態ｋと状態ｙで求めた正規分布の特徴パラメータｐｉの可能性分布関数をＰｋ（ｐｉ）とＰｙ（ｐｉ）とし、状態ｙで求めた正規分布の特徴パラメータの値をｐ’ｉとすると、「状態ｙが状態ｋと同じである」という可能性ｗは次のように求められる。
ａ）ｐ’ｉの平均値ｐ’ｉｍｅａｎとＰｋ（ｐｉ）とのマッチングによるｗの決定、
ｂ）Ｐｙ（ｐｉ）とＰｋ（ｐｉ）とのマッチングによるｗの決定。
なお、Ｐｙ（ｐｉ）とＰｋ（ｐｉ）とのマッチングによるｗを求める式は次に示す。
３）状態変化の判別
状態ｋで求めた正規分布の特徴パラメータｐｉの可能性分布関数ｐｋ（ｐｉ）が得られた後、左右両側の「状態変化が小さい」の可能性分布関数（ｐｃ１（ｐｉ）とｐｃ２（ｐｉ））、および「状態変化が大きい」の可能性分布関数（ｐｄ１（ｐｉ）とｐｄ２（ｐｉ））は図２に示すように決定する。境界値の
におけるｉ，ｊはユーザ入力により決定するが、標準値としてｉ＝３、ｊ＝６とする。
設備診断の場合、正常状態の可能性分布関数をｐｋ（ｐｉ）、注意状態の可能性分布関数をｐｃ１（ｐｉ）とｐｃ２（ｐｉ）、危険状態の可能性分布関数をｐｄ１（ｐｉ）とｐｄ２（ｐｉ）とする。実際の識別時に得られた「正常」と「注意」と「危険」との可能性は図２のように表示する。また、「危険」と判定した場合、警報を出すことも可能である。
（３）情報理論による状態判別法
対象物の基準状態における正規分布の特徴パラメータの確率密度関数をｆｐ０（ｐｉ）し、基準状態以外における正規分布の特徴パラメータの確率密度関数をｆｐｋ（ｐｉ）とする。なお、基準状態以外の時点での状態を「テスト状態」と呼ぶ。テスト状態が基準状態と同じ状態であるか否かは次の「Ｋｕｌｌｂａｃｋ−Ｌｅｉｂｌｅｒ Ｉｎｆｏｒｍａｔｉｏｎ（カルバック情報量、ＫＩ）」と「Ｉｎｆｏｒｍａｔｉｏｎ Ｄｉｖｅｒｇｅｎｃｅ（情報ダイバージェンス、ＩＤ）」で判定できる。
ＫＩｐとＩＤｐによる状態判別法は［１０］に詳細に記述しているので、ここで省略する。（［１０］劉 信芳、豊田利夫、陳 鵬、馮 芳、二保知也：「Ｉｎｆｏｒｍａｔｉｏｎ Ｄｉｖｅｒｇｅｎｃｅによる回転機械の異常診断」、精密工学会誌、Ｖｏｌ．６６、Ｎｏ．１、２０００年、ｐ．１５７−１６２．）
（４）複数の特徴パラメータを統合することによる状態判別法
数個の特徴パラメータを統合して、状態変化の有無や状態予測を行うこともできる。特徴パラメータの統合法は、主成分分析法やＫＬ展開法などがあるが、特徴パラメータの統合法により求めた新たな特徴パラメータは「統合特徴パラメータ」とよぶ。ここで、主成分分析法の例を示す。（［１１］大津，栗田，関田著：パターン認識，朝倉書店、１９９６．）
例えば、設備診断の場合、設備の正常状態下で求めた無次元特徴パラメータ（ｐ１，ｐ２，・・・，ｐｍ）に対して、ｍ個の主成分は次のように表されます。
各主成分ｚ１〜ｚｍは「統合特徴パラメータ」とも呼ぶ。
相関行列は次のように求められます。
ここで、正常状態で収録したｎ組のデータを｛ｐ１ｋ，ｐ２ｋ，・・・，ｐｍｋ｝、ｋ＝１，２，・・・，ｎとすると、
相関行列Ｒの固有ベクトル、λ＝（λ１，λ２，・・・，λｍ）、λ１，≧λ２≧・・・≧λｍとすると、
固有値λｉに対応する式（５９）の係数が求まり、第ｉ番目の主成分は次のように求められます。
主成分を用いて次式のように状態判定を行う。
ここで、Ｋは使用する主成分の数、αは有意水準、χ２（Ｋ，α）は自由度Ｋのカイ２乗分布の確率密度関数が上側確率αに対するパーセント点である。αの決定法は表１同じである。なお、Ｋ＝３に対してα＝０．０５のとき、χ２（３，０．０５）＝７．８１５である。
なお、式（６４）に示す判定方法以外に、特徴パラメータの統合法により求めた統合特徴パラメータ、例えば、式（６３）に示す主成分、を正規分布の確率変数に変換した後、統計検定や可能性理論により状態判定を行うこともできる。
ここで、実例を示す。図３はある回転機械の正常状態時（図３（ａ））と回転軸ミスアライメント状態時（図３（ｂ））に測定した振動の加速度波形データである。図４（ａ）及び図５（ａ）は式（２）〜式（１３）及び式（１８）〜式（２１）に示す１４個の無次元特徴パラメータを６０回求め、正規分布変換を行う前に、Ｋ＝３のときに求めた式（６４）右側の値（６０個）を示す。図４（ｂ）及び図５（ｂ）は式（３１）に基づいて正規分布変換を行った後、Ｋ＝３のときに求めた式（６４）右側の値（９個）を示す。これらの図により、正規分布変換後の結果は正規分布変換前の結果より良好であることがわかる。なお、「良好である」とは、式（６４）左側の値は正常状態の時に小さく、異常状態の時に大きいということである。
４．正規分布の波形データによる状態判別方法
（１）特徴パラメータによる状態判別法
対象物から計測された波形データｘ＊ ｉ（例えば、図６（ａ）、（ｃ））を式（３３）により正規分布の波形データｘ’ｉｋ０（例えば、図６（ｂ）、（ｄ））に変換した後、正規分布の波形データｘ’ｉｋ０を用いて前述の正規分布の特徴パラメータによる状態判別方法で状態判定を行うことができる。なお、この方法は正規分布の波形データμｘｉｋ０とｘ’ｉｋ０に適用できる。
（２）情報理論による状態判別法
基準状態の正規分布の波形データの確率密度関数をｆｘ０（ｘｉ）とし、基準状態以外の時点での正規分布の波形データの確率密度関数をｆｘｋ（ｘｉ）とする。なお、基準状態以外の時点での状態を「テスト状態」と呼ぶ。テスト状態が基準状態と同じ状態であるか否かは次の「Ｋｕｌｌｂａｃｋ−Ｌｅｉｂｌｅｒ Ｉｎｆｏｒｍａｔｉｏｎ（カルバック情報量、ＫＩ）」と「Ｉｎｆｏｒｍａｔｉｏｎ Ｄｉｖｅｒｇｅｎｃｅ（情報ダイバージェンス、ＩＤ）」で判定できる。
ＫＩｘとＩＤｘによる状態判別法はに詳細に記述しているので、ここで省略する。（［１２］劉 信芳、豊田利夫、陳 鵬、馮 芳、二保知也：「Ｉｎｆｏｒｍａｔｉｏｎ Ｄｉｖｅｒｇｅｎｃｅによる回転機械の異常診断」、精密工学会誌、Ｖｏｌ．６６、Ｎｏ．１、２０００年、ｐ．１５７−１６２．）
ここで、実例を示す。図６（ａ）はある回転機械の正常状態時に測定した振動加速度の波形データである。図６（ｃ）は同じ回転機械のアンバランス時に測定した波形データである。
平均値以上の波形データｘｉ＋、平均値以下の波形データｘｉ−、及び絶対値データ｜ｘｉ｜を正規分布の波形データに変換した後、それぞれｘ’ｉ＋、ｘ’ｉ−、及び｜ｘｉ｜’で表す。ｘ’ｉ＋、ｘ’ｉ−、及び｜ｘｉ｜’の平均値と標準偏差は次の通りである。
ｘ’ｉ＋：正常時の平均値＝０．６９、アンバランス時の平均値＝１．９６
正常時の標準偏差＝１．０３、アンバランス時の標準偏差＝２．３７
ｘ’ｉ−：正常時の平均値＝−０．３７、アンバランス時の平均値＝−０．６６
正常時の標準偏差＝０．４７、アンバランス時の標準偏差＝０．６９
｜ｘｉ｜’：正常時の平均値＝１．４５、アンバランス時の平均値＝３．２９
正常時の標準偏差＝１．５２、アンバランス時の標準偏差＝３．５３
（参考文献８）によると、過検出率α＝０．１５、見逃率β＝０．１５とすると、ＫＩ＞１．２１或いはＩＤ＞２．４３ならば、８５％ｎの確率で「テスト状態は基準状態と違う」と判定できる。上記のｘ’ｉ＋、ｘ’ｉ−、及び｜ｘｉ｜’のＫＩとＩＤは次の通りである。
ｘ’ｉ＋：ＫＩ＝０．５１９， ＩＤ＝２．２４
ｘ’ｉ−：ＫＩ＝０．５３， ＩＤ＝０．４０
｜ｘｉ｜’：ＫＩ＝０．５７， ＩＤ＝２．６７
｜ｘｉ｜’のＩＤ＞２．４３から、８５％ｎの確率で「図６（ａ）の状態は図６（ｃ）の状態と違う」と判定できる。
また、式（４４）と式（４５）に示す平均値と分散の検定により状態判定を行ってもよい。
更に、上記のｘ’ｉ＋、ｘ’ｉ−、及び｜ｘｉ｜’を用いて、式（２）〜式（２５）で特徴パラメータを求めて、統計検定や可能性理論及び特徴パラメータの統合法により状態判定を行うこともできる。
５．状態予測
正規分布の特徴パラメータｐｉが求められた後、従来の状態予測手法を用いて、測定対象の状態を予測することが行える（［１３］石川真澄、武藤博道：予測手法、計測と制御、１９８２．３．［１４］Ｏｇａｗａ，Ｍ．：Ｔｉｍｅ ｓｅｒｉｅｓ ａｎａｌｙｓｉｓ ａｎｄ ｓｔｏｃｈａｓｔｉｃ ｐｒｅｄｉｃｔｉｏｎ，Ｂｕｌｌ．Ｍａｔｈ．Ｓｔａｔ．，８，８−７２，１９５８．［１５］Ｂ．ピカニオル著、小村賢二、柴山宮恵子訳：予測のための統計学、晃洋書房、１９８７．）
図７は状態予測の方法を示す。各測定時点（ｘ１〜ｘ７）で測定した特徴パラメータ値、或いは、主成分値を正規確率分布に従うように変換した後、回帰解析により予測曲線とその信頼区間を求め、「寿命限界」との交点で「最短寿命」、「平均寿命」及び「最長寿命」を求める。
図８はある回転機械が正常状態からアンバランス状態へ変化していく間に８回測定した波形例である。測定１は正常状態であり、測定８は最も程度の重いアンバランス状態である。図９は式（２）〜式（１３）で求めた、各測定における特徴パラメータの平均値を示す。この例でわかるように、特徴パラメータは異常状態の程度が重くなっていくにつれて値が単調増加の特徴パラメータ（例えば、ｐ４，ｐ５）があれば、値が単調減少の特徴パラメータ（例えば、ｐ２，ｐ８，ｐ９，ｐ１０）もある。また、異常状態の程度が重くなっても値がほぼ不変な特徴パラメータ（例えば、ｐ１，ｐ７）もある。従って、値が単調増加（或いは、単調減少）の特徴パラメータを選んで、寿命予測を行うべきである。また、主成分により寿命予測を行う時、値が単調増加の特徴パラメータのみ（或いは、単調減少の特徴パラメータのみ）を選んで、主成分を求めて寿命予測を行うべきである。
上記に述べた波形データと特徴パラメータを正規分布の確率変数に変換する手法、状態判定手法及び状態予測手法を実現するための計測と処理の流れは図１６（ａ）に示す。また、図１６（ａ）を実現するための波形データ計測及び状態判定装置の回路は図１７に示す。
６．脈動信号から特異成分の検出法及び状態判定法
ここで、ガスエンジン失火診断の例を用いて説明する。
図１０（ａ）、図１１（ａ）、図１２（ａ）、図１３（ａ）はガスエンジン失火（図中のピーク値の異常個所）が発生したときのシリンダ圧力波形データである。なお、図１２（ａ）、図１３（ａ）は多数の異常ピークが連続して発生する例を示す。このような脈動信号に発生する局所異常（特異成分）に対して、オンライン特異点の検出及び状態判定の流れは図１５に示す。以下、その手順を説明する。
１）．診断ための準備（図１５（ａ））
（１）診断対象の脈動信号と回転パルス信号を同時に計測する。なお、回転パルス信号は周期パルス信号とも言い、脈動信号の各ピーク値のタイミングを決めるときに用いられる。
（２）雑音を除去するために，ローパスフィルタリングを行う。ローパスフィルタのカットオフ周波数ｆＬは次のように決定する。
ここで、ｎは軸の回転数（ｒｐｍ）、ｚは１回転毎のピーク数（ピーク／１回転）、ｆ０は余裕周波数（＞ｎ／６０、フィルタリング後のノイズ除去効果を観察することにより決定する）である。
（３）ローパスフィルタリング後の包絡線波形データ、または、ピーク波形データ、または、移動平均波形データを求める。
（４）上記の包絡線波形データ、または、ピーク波形データ、または、移動平均波形データの正規化を次式のように行う。
ここで、ｘ（ｔ）は正規化後の波形データ、ｘ’（ｔ）は元の波形データ、μ’（ｔ）はｘ’（ｔ）の平均値、ｓはｘ’（ｔ）の標準偏差である。なお、μ’（ｔ）とｓを求めるとき、正常状態の波形データが望ましいが、正常状態の波形データでなくてもよい。また、μ’（ｔ）とｓを求めるとき、運転条件（回転数や負荷）をできるだけ判定のときと同じように設定する。但し、運転条件が多少変化しても判定結果に大きく影響しない。
（５）回転パルス信号のピークと脈動信号のピークとの相対位置関係を調べて置き、異常を識別するための閾値を決定する。閾値の決定は次のように考える。
のとき、異常が発生したと判定する。ここで、ｋ（＝２〜４）は診断対象物によって決定する。
｜ｘ（ｔ）｜がｋσより大きい間、回転パルス信号に対応する脈動信号のピークの所に異常があると判定する。
２）．オンライン診断（図１５（ｂ））
（１）判定対象物の脈動信号と回転パルス信号を同時に計測する。
（２）ローパスフィルタリングを行う。ローパスフィルタのカットオフ周波数ｆＬは式（６６）のように決定する。
（３）ローパスフィルタリング後の包絡線波形データ、または、ピーク波形データを求める。
（４）上記の包絡線波形データ、または、ピーク波形データの正規化は式（６８）により求める。
（５）正規化された包絡線波形データ、または、ピーク波形データの絶対値は閾値（ｋσ）より大きいか否かを監視する。大きくなければ、異常がないと判定する。大きければ、異常が発生したと判定し、｜ｘ（ｔ）｜がｋσより大きい間、前もって測って置いた回転パルス信号と脈動信号のピークとの関係により、異常の箇所を判定する。
上記の手順に基づく信号処理の例を図１０、図１１、図１２、図１３に示す。
図１０（ａ）（ｂ）、図１１（ａ）（ｂ）、図１２（ａ）（ｂ）、図１３（ａ）（ｂ）ガスエンジン失火（図中の異常箇所）が発生したときのシリンダ圧力の時系列波形データとスペクトルである。
図１０（ｃ）（ｄ）、図１１（ｃ）（ｄ）、図１２（ｃ）（ｄ）、図１３（ｃ）（ｄ）はローパスフィルタにより雑音を除去した時系列波形データとスペクトルである。
図１０（ｅ）、図１１（ｅ）、図１２（ｅ）、図１３（ｅ）は回転パルス波形データである。
図１０（ｆ）、図１２（ｆ）は包絡線波形データである。
図１１（ｆ）、図１３（ｆ）はピーク波形データである。
図１０，１１，１２，１３の（ｆ）で分かるように、異常（ピークの欠落、または、小さくなる）が発生した場合、その絶対値｜ｘ（ｔ）｜が２σより大きくなることが判明できる。
また、絶対値｜ｘ（ｔ）｜が２σより大きくなる箇所と回転パルス波形データとの対応関係により異常箇所を判明できる。
上記に包絡線波形データとピーク波形データの例を示したが、図１４に示すように、移動平均波形データ（ｕｉ）を求めたら、ｘ（ｔ）の替わりにｕ（ｔ）を用いて上記の手順と同様に脈動信号の特異点検出ができる。なお、カットオフ周波数（ｆｃ）は式（６７）により求めたら、移動平均の点数（Ｍ）が求められる。
上記に述べた脈動信号の特異成分の検出と判定方法を実現するための計測と処理の流れは図１６（ｂ）に示す。また図１６（ｂ）を実現するための脈動信号計測及び状態判定装置の回路は図１８に示す。
７．状態判定装置、または、オンライン状態判定システムの処理流れ
状態判定装置装置、または、オンライン状態判定システムの処理流れは図１６に示す。図１６（ａ）に示すように、測定波形データからノイズを除去し、特徴パラメータを求め、その特徴パラメータを正規分布の確率変数に変換した後、統計検定や可能性理論や特徴パラメータの統合により状態変化の有無判定を行う。これらの処理は計算機または専用装置で実現できる。また、図１６（ｂ）に示すように、脈動信号からの特異成分検出及び状態判定は、ローパスフィルタ、包絡線（または、ピーク値）はハードウェアで実現できる。数値演算装置（または、計算機）は正規化したｘ（ｔ）、｜ｘ（ｔ）＞ｋσ｜の判定、特異箇所の検出と判定および判定結果の表示を行えばよいから、リアルタイムに行える。one. Technical field
The present invention relates to a state determination device and an online state monitoring / diagnosis system for determining whether or not a state of an object has changed in equipment diagnosis, medical diagnosis, and the like.
two. Background art
(1) Conventionally, state determination is performed on the assumption that the probability density function of waveform data measured from an object to be monitored follows a “normal distribution” in a normal state ([1] Patent Publications 2000-171291, [2]. Toshio Toyoda, Tomoya Niho: “Abnormal diagnosis of rotating machinery using only normal vibration waveform”, Journal of Japan Institute of Equipment Management, Vol.11, No.1, 1999, p.4-11. [3] Nobuyoshi Liu, Toshio Toyoda, Satoshi Chen, Yoshiaki Tsuji, Tomoya Niho: “Abnormal Diagnosis of Rotating Machinery by Information Divergence”, Journal of Precision Engineering, Vol. 66, No. 1, 2000, p.157-162.)
(2) Statistical test is performed on the assumption that the characteristic parameter used to determine whether there is a change in the state of the object is a normally distributed random variable ([4] Yoshiki Kawabe, Toshio Toyoda, Toru Eguchi, Kazuhisa Fukuda: (Management of deterioration tendency of rotating machinery using dimensionless sign parameters (III) ", Proceedings of the 1995 Annual Conference of the Japan Institute of Equipment Management, Vol. 2, 1995, p. 32.)
When the signal measured from the state monitoring object is a pulsation signal, it is difficult to determine the state by using a characteristic parameter, a probability density distribution function, or the like. Conventionally, as a method of detecting a specific component (or abnormal component) included in a pulsation signal,
(3) Local anomaly detection methods using time-frequency analysis methods (wavelet analysis, Wigna distribution analysis, short-time FFT, etc.) ([5] Akira Tadashi, Tomoyuki Nakabori, Hiroaki Kawabata: "High-speed wavelet transform and its Application to electroencephalogram analysis ”, Transactions of the Japan Society of Mechanical Engineers (C), Vol. 65, No. 632, 1999, p. 1915-1921.)
(4) A value sampled at a predetermined interval according to the peak value of the pulsation signal or the like is determined to be abnormal when the difference from the previously sampled value is equal to or less than a predetermined value. ([6] JP-A-6-317215)
(5) When the variation (standard deviation) of the representative value for each sampling is greater than or equal to a predetermined threshold at a predetermined interval according to the peak value of the pulsation signal, it is determined that there is an abnormality. ([7] Japanese Patent Laid-Open No. 7-293311)
three. Problems to be solved by the invention
However, the conventional methods have the following problems.
In the first and second methods, it is assumed that the waveform data and the characteristic parameter values follow a normal distribution. However, the waveform data and the characteristic parameter values do not always follow the normal distribution even in a normal state. If it is assumed that the normal distribution is followed and the state change of the object is determined by the statistical test, an erroneous determination is caused.
In the third method, since the time required for the processing is long, it is difficult to detect the singular point of the real time.
In the fourth and fifth methods, when waveform data is sampled at a predetermined interval, when the peak position changes due to a change in the rotation speed or a change in load, the sampling value changes greatly even if it is normal. Bring. Further, if the variation (standard deviation) of the representative value for each sampling is calculated, real-time singular detection is difficult.
Four. Means for solving the problem
In order to solve the problems as described above, in the present invention, after the measured waveform data is converted into waveform data of a known probability distribution (for example, normal distribution), or a feature calculated from the waveform data After the parameters are converted into characteristic parameters of a known probability distribution (for example, normal distribution), the state of the object is determined by statistical testing, possibility theory, information theory, or the like.
Also, if the signal of the measured object is a pulsation signal, it is difficult to extract the singular component and determine the state using feature parameters, spectrum, etc. The state is determined by detecting a singular component in the signal from the envelope waveform data and the periodic pulse. Alternatively, the peak waveform data of the pulsation signal from which the noise has been removed is obtained, the peak waveform data is normalized, and a specific component in the pulsation signal is detected from the normalized peak waveform data and periodic pulse waveform data. Determine the state.
Five. The invention's effect
Since the waveform data and feature parameters measured for monitoring the condition of the object do not necessarily follow a normal distribution, it is assumed that these waveform data or feature parameters follow a normal distribution, and whether there is a change in the state of the object If it is determined by the test or the state is predicted, a large error is caused. In the present invention, after the measured waveform data is converted into waveform data of a known probability distribution (for example, normal distribution), or the characteristic parameter calculated from the waveform data is converted to the characteristic parameter of the known probability distribution (for example, normal distribution). After conversion to, the state of the object is determined by statistical tests, possibility theory, information theory, and the like. Therefore, the method of the present invention has a higher accuracy of state determination or state prediction than when it is assumed that the conventional normal distribution is followed.
In the case of a pulsation signal, when the measured waveform data of the pulsation signal is missing some peaks, or when a small singular component exists locally in the waveform data, the conventional judgment method (feature parameter, spectrum, probability density Etc.) are difficult to detect. In the present invention, the envelope waveform data of the pulsation signal from which noise has been removed is obtained, and the state is determined by detecting the singular component in the signal from the normalized envelope waveform data and the periodic pulse. Alternatively, the peak waveform data of the pulsation signal from which the noise has been removed is obtained, the peak waveform data is normalized, and a specific component in the pulsation signal is detected from the normalized peak waveform data and periodic pulse waveform data. Determine the state. In the method of the present invention, main signal processing can be realized by hardware, and the burden on the numerical processing device (computer) is small, so that it is possible to detect abnormality in real time. Further, since the waveform data in the normal state is not used as a reference, even if the pulsation period changes, the influence on the detection and determination results is small.
Six. BEST MODE FOR CARRYING OUT THE INVENTION
1. About feature parameters
The feature parameters used for state determination include time domain and frequency domain feature parameters. Some frequency domain characteristic parameters are defined in (Reference 1). ([8] Chen, Toshio Toyoda: Self-reorganization of frequency domain feature parameters by genetic programming, Japan Machine Papers (C), Vol. 65 No. 633, pp. 1946-1953, 1998.) Here, characteristic parameters in the time domain will be described in detail.
The time domain feature parameters used to determine whether or not the state of the object has changed are as follows.
1) Dimensionless feature parameters
From the measured time-series waveform data, waveform data in the low, medium, and high frequency regions are extracted using a filter. The extracted waveform data x (t) is normalized by the following equation.
The dimensionless feature parameters used in the past are shown in Equations (2) to (13) ([9] Peng CHEN, Toshio TOYOTA, Yueton LIN, Feiyu WANG: FAIURE DIAGNOSIS OF MACHINEYN IN MION IN MISTION OF MISON GENETIC ALGORITHMS, International Journal of Intelligent Control and System, Vol. 3, No. 4, pp. 571-585, 1999.)
here,
Is an absolute average value, and N is the total number of data.
Is the standard deviation.
Where μpIs the average value of the maximum value (peak value) of the waveform.
Where | μmax| Is an average value of 10 maximum values of the waveform.
Where σpIs the standard deviation of the maximum value.
Where μLAnd σLAre the average value and standard deviation value of the local minimum values (valley values).
Equations (2) to (13) are conventional feature parameters, but “section feature parameters” are newly proposed as in Equations (14) to (21) for easy high-speed calculation by numerical calculation. To do.
Where xi≧ kσ and k can be set arbitrarily. For example, k = 0.5, 1, and 2. μk1Is xiIs the average value. t can be set arbitrarily. For example, t = 0.5, 1, 2, 3, 4.
Where xhAlthough ≦ −hσ and h can be set arbitrarily, for example, k = 0.5, 1, and 2. μh1Is xhIs the average value. t can be set arbitrarily. For example, t = 0.5, 1, 2, 3, 4.
However, h0Is the frequency at which the time series waveform closes to 0 level per unit time, hpIs the number of peaks of the time series waveform per unit time.
However, hnσIs the frequency at which the waveform closes the nσ level per unit time. Although n can be set arbitrarily, for example, n = 0.5, 1, 2.
However, h-NσIs the frequency at which the waveform closes the -nσ level per unit time. It is. Although n can be set arbitrarily, for example, n = 0.5, 1, 2.
2) Dimensional feature parameters
When the dimensional feature parameter is calculated, the measured waveform data is not normalized as in Expression (1).
Where | xi｜pIs the peak value (maximum value) of the absolute value of the waveform data, NpIs the total number of peak values.
Many feature parameters can be defined in addition to the above feature parameters. However, when applying this method, try the above feature parameters first, and if the effect of state identification is not good, add other feature parameters. Define it.
2. Convert waveform data and feature parameters to specified probability distribution
X measured waveform data* iThe characteristic parameter calculated from the waveform data is p* iRepresented by x* iOr p* iWhen performing state determination and state prediction by statistical theory using* iOr p* iIt is necessary to know in advance what probability distribution follows. But x* iOr p* iIt is often not known in advance what kind of probability distribution follows. Therefore, if the specified known probability distribution function is Σ, x* iOr p* iA random variable x following ΣiOr piCan be converted to
Where Fxi(X* i) And Fpi(P* i) Are each x* iAnd p* iCumulative probability distribution (or cumulative frequency distribution) and Σ-1Is the inverse function of Σ. For example, Σ is a normal distribution, Weibull distribution, exponential distribution, gamma distribution, or the like. x* iOr p* iXiOr piAfter conversion to, state determination and state prediction are performed using statistical tests and the like.
The original waveform data x* iAre divided into four types as shown in FIG. That is, data x greater than the average valuei +, Data smaller than average xi-, Absolute value data after normalization by equation (1) | xi| And whole waveform data x* iAEach cumulative probability distribution (or cumulative frequency distribution) is Fxi +(Xi +), Fxi-(Xi-), F| Xi |(| Xi|) And FxiA(X* iA), But below, unless otherwise specified,xi(X* i). The characteristic parameter p* iIs the normalized xi +And xi-, | Xi| And x* iACalculated with either
Here, as an example, if Σ is specified as a normal distribution, x* iOr p* iWill be described in detail.
The probability density function f (t) of the normal distribution is expressed by the following equation.
Here, μ is an average value of the random variable t, and σ is a standard deviation.
(1) When referring to the standard condition
Determine the reference state of the measurement target, for example, the state at the first measurement, and the waveform data x at this time* i0Or feature parameter p* i0The probability density function and cumulative probability distribution function ofx0(X* i0) And Fx0(X* i0) Or fp0(P* i0) And Fp0(P* i0). The probability density function of a standard normal distribution with an average value μ of 0 and a standard deviation σ of 1 is represented by φ (xi), The standard normal distribution probability distribution function is Φ (xi). Note that x is discrete data.* i0Or p* i0“Frequency distribution function” and “His and Gram” may be used in place of the probability density function, but will be described below using “probability density function”.
1) Method to convert to mean value according to normal distribution
X in state k* ikAnd p* ikIs the mean value μ according to the normal distributionxikoAnd μpik0Is converted to Note that the state k is an arbitrary state and includes a reference state.
Where Φ-1Is the inverse function of Φ and σxi0And σpi0Are each converted to a normal distribution x* ikAnd p* ikThe standard deviation is obtained from the following equation.
Average value μxik0And μpik0The state judgment and state prediction are performed using.
2) Direct conversion method
X in state k* ikAnd p* ikIs converted to a normally distributed random variable by the following equation.
Where SxkAnd SpkAre each x* ikAnd p* ikStandard deviation of μxkAnd μpkAre each x* ikAnd p* ikIs the average value. x ’ik0Or p ’ik0The state judgment and state prediction are performed using.
(2) When reference conditions are not referenced
First, waveform data x in state k* ikAnd feature parameter p* ikAnd the probability density function (or frequency distribution) ofxk(X* ik) And fpk(P* ik) And the probability distribution function (or cumulative frequency distribution) as Fxk(X* ik) And Fpk(P* ik).
1) Direct variable conversion method
Waveform data x* ikAnd feature parameter p* ikIs converted to a normally distributed random variable by the following equation.
Where SxkAnd SpkAre each x* ikAnd p* ikWith the standard deviation ofxkAnd μpkIs x* ikAnd p* ikAnd the average value. x ’ikOr p ’ikThe state judgment and state prediction are performed using.
2) Method to convert to average value according to normal distribution
The average value according to the normal distribution is obtained by the following equation.
Where σxikAnd σpikIs obtained by the following equation.
μxikAnd μpikFollows a normal distribution, so μxikOr μpikThe state judgment and state prediction are performed using.
3) Indirect conversion method
M feature parameters p* ikThe minimum and maximum values of (p* ik)minAnd (p* ik)maxAnd (P* ik)minTo (p* ik)maxIs divided into N equally spaced intervals. The representative value (for example, median) of each section is p* ikjAnd Here, j = 1 to N. p* ikjP* ikSubstituting into equation (36) or equation (38) instead of N p ″ikOr μ ’pikIs obtained. p ″ikOr μ ’pikThe state judgment and state prediction are performed using.
Similarly, M pieces of waveform data x* ikThe minimum and maximum values of (x* ik)minAnd (x* ik)maxAnd (X* ik)minTo (x* ik)maxIs divided into N equally spaced intervals. The representative value (eg, median) of each section is x* ikjAnd Here, j = 1 to K. x* ikjX* ikSubstituting into equation (35) or equation (34) instead ofikOr μ ’xikIs obtained. x ″ikOr μ ’xikThe state judgment and state prediction are performed using.
4) Feature parameter p* ikTo find the interval average of
N feature parameters p* ikThen, the average value of the j-th set is obtained as follows by dividing into M sets.
Where NjIs p in the jth set* ikIs the number of μk (J)Since approximately follows a normal distribution, μk (J)The state judgment and state prediction are performed using.
Waveform data x* ikIs converted into waveform data of a normal probability distribution by Equation (29), Equation (33), Equation (35), and Equation (37).xik0, X ’ik0, X ’ik, Μxik, X ’ik, Μ ’xikIs called “normally distributed waveform data”. The characteristic parameter p* iIs converted to a normal random variable by Equation (30), Equation (34), Equation (36), and Equation (38).pik0, P ’ik0, P ’ik, Μpik, P "ik, Μ ’pik, Μk (J)Is called “normal distribution feature parameter”.
3. State discrimination method using feature parameters of normal distribution
Here, a method for discriminating the state of the object using the characteristic parameter of the normal distribution will be described.
(1) Discrimination based on statistical theory
1) Test of mean value of feature parameter of normal distribution
The characteristic parameters of the normal distribution obtained in state k and state y are pikAnd piyAnd Here, i = 1 to M and M represent the total number of feature parameters of the normal distribution to be used. pikAnd piyAnd the average value ofikAnd μiy, PikAnd piyAnd the standard deviation of eachikAnd SiyAnd In general, J pjThe average value μ and the standard deviation S are calculated by the following equations.
μikAnd μiyThe equality test is performed as follows (see Reference 3)
(Reference 3) A. Brownlee. Statistical Theory and Methodology in Science and Engineering, Second Edition, The University Chicago, 1965
Is established, the significance level αikAnd μiyAre not equal. " Where tα / 2(J-1) is a percentage point with respect to the lower probability α / 2 of the probability density function of the t distribution with the degree of freedom J-1.
2) Test of variance of characteristic parameter of normal distribution
SikAnd SiyAre tested as follows (see Reference 4).
(Reference 4) A. Brownlee. Statistical Theory and Methodology in Science and Engineering, Second Edition, The University Chicago, 1965
Is established, the significance level α is “SikAnd SiyAre not equal. " Where Fα / 2(J-1, J-1) is a percentage point with respect to the lower probability α / 2 of the probability density function of the F distribution with the degree of freedom J-1.
By checking whether the expression (44) or the expression (45) is satisfied when the significance level α is changed, the degree of the state change of the state y with respect to the state k is determined. Table 1 shows an example in which the degree of state change is determined by the significance level α. In the case of equipment diagnosis, if the state k is a normal state, the determination of whether the state y is a normal state, a caution state, or a dangerous state is “normal” (α1), “Caution” (α2) And “danger” (α3) Can be set and verified. That is, Formula (44) or Formula (45) is α1If it does not hold at this time, it is determined as “normal”. Moreover, Formula (44) or Formula (45) is α2If it is established at3If it is established at this time, it is determined as “dangerous”. In addition, the numerical range of α in Table 1 is an example, and is determined by the importance of the facility.
When the state determination is performed using several feature parameters, the determination result follows the determination result of the feature parameter indicated that the state change is the largest. For example, three feature parameters p1, P2, P3When judging using1The judgment result of is “Caution”, p2The judgment result is “normal”, p3If the determination result is “danger”, the final determination result is “danger”.
3) Judgment by confidence interval
The average value of the normal distribution feature parameters obtained from the waveform data measured at the reference time is μi0The average value of the normal distribution feature parameters obtained from the waveform data measured at other times is μikThen μi0The confidence interval is given by
Where tα / 2(J-1) is a percentage point with respect to the lower probability α / 2 of the probability density function of the t distribution with the degree of freedom J-1. Si0Is obtained from waveform datai0Is the standard deviation. μikIf within the interval shown in Equation (46), μ has a probability of 1−α.i0There is no difference. μi0The 99% confidence interval is approximately as follows when J> 10.
Therefore, μikExceeds the range of equation (47), with a probability of 99% μi0Judge that it is different. Furthermore, the μ obtained from the measured waveform dataikIs obtained by the following equation.
Where SikIs obtained from waveform dataikIs the standard deviation.
By substituting α in Table 1 into equation (46), the following confidence interval is obtained.
μikThe state is determined depending on whether or not is in these sections.
When state determination is performed using several feature parameters, the final determination result is the same as the contents described in the equations (49), (50), and (51).
(2) Discrimination based on possibility theory
1) Creation of possibility distribution function
Normal distribution feature parameter p with waveform data in state kiAfter calculating the value of piProbability density function fk(Pi) Possibility distribution function Pk(Pi) Is obtained by equation (52). According to the possibility theory, piThe probability distribution function is obtained for any probability distribution. piIs a normal distribution, the N-stage probability distribution function pk(Pi) Is determined as follows (see Reference 5).
here,
However, in the above formula, pix= Min {pi} + Xx (max {pi} -Min {pi}) / N, x = 1 to N, SiIs piStandard deviation of μiHano piAverage value.
(Reference 5) Davis: HANDBOOK OF GENETIC ALGORITHMS, Van Northland Reinhold, A Division of Wadsworth, Inc (1990)
2) How to find the possibility
As shown in FIG. 2, the characteristic parameter p of the normal distribution obtained in the state k and the state yiThe probability distribution function of Pk(Pi) And Py(Pi) And the characteristic parameter value of the normal distribution obtained in the state y is p ′iThen, the possibility w that “the state y is the same as the state k” is obtained as follows.
a) p ’iAverage value of p ′imean and Pk(Pi) To determine w by matching with
b) Py(Pi) And Pk(Pi) To determine w.
Py(Pi) And Pk(PiThe formula for obtaining w by matching with) is as follows.
3) State change discrimination
Feature parameter p of normal distribution obtained in state kiPossibility distribution function pk (pi), The probability distribution function (pc1(Pi) And pc2(Pi)), And the probability distribution function of “large state change” (pd1(Pi) And pd2(Pi)) Is determined as shown in FIG. Boundary value
I and j are determined by user input, but i = 3 and j = 6 are standard values.
In the case of equipment diagnosis, the probability distribution function in the normal state is expressed as pk (pi), The probability distribution function of attention state pc1(Pi) And pc2(Pi), P.d1(Pi) And pd2(Pi). The possibility of “normal”, “caution”, and “danger” obtained at the time of actual identification is displayed as shown in FIG. If it is determined as “dangerous”, an alarm can be issued.
(3) State discrimination method by information theory
The probability density function of the characteristic parameter of the normal distribution in the reference state of the object is fp0(PiAnd the probability density function of the characteristic parameter of the normal distribution in a state other than the reference state is expressed as fpk(Pi). A state at a time other than the reference state is referred to as a “test state”. Whether or not the test state is the same as the reference state can be determined by the following “Kullback-Leibler Information (Kalback information amount, KI)” and “Information Divergence (information divergence, ID)”.
KIpAnd IDpSince the state determination method based on is described in detail in [10], it is omitted here. ([10] Liu Nobuyoshi, Toyoda Toshio, Chen Hao, Zhao Fang, Niho Tomoya: “Abnormal Diagnosis of Rotating Machinery by Information Divergence”, Journal of Precision Engineering, Vol. 66, No. 1, 2000, p.157- 162.)
(4) State determination method by integrating a plurality of feature parameters
By integrating several feature parameters, it is also possible to perform the presence / absence of state change and state prediction. The feature parameter integration method includes a principal component analysis method and a KL expansion method, and a new feature parameter obtained by the feature parameter integration method is referred to as an “integrated feature parameter”. Here, an example of the principal component analysis method is shown. ([11] Otsu, Kurita, Sekida: Pattern recognition, Asakura Shoten, 1996.)
For example, in the case of equipment diagnosis, the dimensionless feature parameter (p1, P2, ..., pm), M principal components are expressed as follows.
Each main component z1~ ZmIs also called “integrated feature parameter”.
The correlation matrix is obtained as follows:
Here, n sets of data recorded in a normal state are represented by {p1k, P2k, ..., pmk}, K = 1, 2,..., N,
Eigenvector of correlation matrix R, λ = (λ1, Λ2, ..., λm), Λ1, ≧ λ2≧ ・ ・ ・ ≧ λmThen,
Eigenvalue λiThe coefficient of equation (59) corresponding to is obtained, and the i-th principal component is obtained as follows.
The state is determined using the principal component as in the following equation.
Where K is the number of principal components used, α is the significance level, χ2(K, α) is a percentage point of the probability density function of the chi-square distribution with K degrees of freedom relative to the upper probability α. The method for determining α is the same as in Table 1. When α = 0.05 with respect to K = 3, χ2It is (3,0.05) = 7.815.
In addition to the determination method shown in Formula (64), after converting the integrated feature parameters obtained by the feature parameter integration method, for example, the principal component shown in Formula (63), into random variables of normal distribution, State determination can also be made by possibility theory.
Here is an example. FIG. 3 shows acceleration waveform data of vibrations measured in a normal state (FIG. 3A) and a rotational axis misalignment state (FIG. 3B) of a certain rotating machine. 4 (a) and 5 (a) obtain the 14 dimensionless feature parameters shown in equations (2) to (13) and (18) to (21) 60 times, and perform normal distribution conversion. Before, the value (60 pieces) on the right side of Expression (64) obtained when K = 3 is shown. FIG. 4B and FIG. 5B show values (9) on the right side of Expression (64) obtained when K = 3 after performing normal distribution conversion based on Expression (31). From these figures, it can be seen that the result after the normal distribution conversion is better than the result before the normal distribution conversion. “Good” means that the value on the left side of the equation (64) is small in the normal state and large in the abnormal state.
4). State discrimination method using normal distribution waveform data
(1) State discrimination method using feature parameters
Waveform data measured from the object x* i(For example, FIGS. 6A and 6C are converted into the waveform data x ′ of the normal distribution by the equation (33).ik0After conversion into (for example, FIGS. 6B and 6D), waveform data x ′ having a normal distributionik0The state determination can be performed by the above-described state determination method using the characteristic parameter of the normal distribution. This method uses normal distribution waveform data μxik0And x ’ik0Applicable to.
(2) State discrimination method by information theory
The probability density function of the normal distribution waveform data in the reference state is expressed as fx0(Xi) And the probability density function of the waveform data of the normal distribution at a time other than the reference state is fxk(Xi). A state at a time other than the reference state is referred to as a “test state”. Whether or not the test state is the same as the reference state can be determined by the following “Kullback-Leibler Information (Kalback information amount, KI)” and “Information Divergence (information divergence, ID)”.
KIxAnd IDxSince the state discrimination method based on is described in detail in, it is omitted here. ([12] Liu Nobuyoshi, Toyoda Toshio, Chen Sung, Zhao Fang, Futoshi Tomoya: “Abnormal Diagnosis of Rotating Machinery by Information Divergence”, Journal of Precision Engineering, Vol. 66, No. 1, 2000, p.157- 162.)
Here is an example. FIG. 6A shows vibration acceleration waveform data measured when a certain rotating machine is in a normal state. FIG. 6C shows waveform data measured when the same rotating machine is unbalanced.
Waveform data above average value xi +Waveform data x below average valuei-And absolute value data | xi| Is converted into normal distribution waveform data, and then x ′i +, X ’i-And | xi| '. x ’i +, X ’i-And | xiThe average value and standard deviation of | 'are as follows.
x ’i +: Average value at normal time = 0.69, Average value at unbalance = 1.96
Standard deviation at normal time = 1.03, standard deviation at unbalance = 2.37
x ’i-: Average value at normal time = -0.37, average value at unbalance = -0.66
Standard deviation at normal time = 0.47, standard deviation at unbalance = 0.69
| Xi| ′: Average value at normal time = 1.45, Average value at unbalance = 3.29
Standard deviation at normal time = 1.52, standard deviation at unbalance = 3.53
According to (Reference 8), if the over-detection rate α = 0.15 and the oversight rate β = 0.15, if KI> 1.21 or ID> 2.43, the probability of 85% n It can be determined that the state is different from the reference state. X ’abovei +, X ’i-And | xiThe KI and ID of | 'are as follows.
x ’i +: KI = 0.519, ID = 2.24
x ’i-: KI = 0.53, ID = 0.40
| Xi| ′: KI = 0.57, ID = 2.67
| XiFrom the ID of | ′> 2.43, it can be determined that “the state of FIG. 6A is different from the state of FIG. 6C” with a probability of 85% n.
Further, the state determination may be performed by testing the average value and variance shown in Expression (44) and Expression (45).
Furthermore, the above x 'i +, X ’i-And | xi| 'May be used to obtain the characteristic parameters by the equations (2) to (25), and the state determination may be performed by statistical test, possibility theory, and feature parameter integration method.
5). Condition prediction
Normal distribution feature parameter piCan be predicted using a conventional state prediction method ([13] Masumi Ishikawa, Hiromichi Muto: Prediction method, measurement and control, 1982.3. [14] Ogawa) , M .: Time series analysis and stochastic prediction, Bull. Math. Stat., 8, 8-72, 1958. [15] B. Picaniol, Kenji Komura, Keiko Shibayama Translation: Statistics for Prediction, Yoyo Shobo, 1987.)
FIG. 7 shows a state prediction method. Each measurement time point (x1~ X7) Measured characteristic parameter values or principal component values so as to follow a normal probability distribution, then a prediction curve and its confidence interval are obtained by regression analysis, and the “shortest life”, “ Obtain the “average life” and “longest life”.
FIG. 8 is an example of a waveform measured eight times while a certain rotating machine changes from a normal state to an unbalanced state. Measurement 1 is normal and measurement 8 is the most severe unbalance. FIG. 9 shows the average value of the characteristic parameter in each measurement, which is obtained by the equations (2) to (13). As can be seen from this example, the characteristic parameter has a monotonically increasing characteristic parameter (for example, p) as the degree of abnormal state increases.4, P5) Is a feature parameter whose value is monotonically decreasing (eg, p2, P8, P9, P10There is also. In addition, a characteristic parameter whose value is almost unchanged even if the degree of abnormal state becomes heavy (for example, p1, P7There is also. Therefore, the lifetime prediction should be performed by selecting a characteristic parameter whose value is monotonously increasing (or monotonically decreasing). In addition, when performing life prediction based on the principal component, only the characteristic parameter whose value is monotonously increasing (or only the characteristic parameter being monotonically decreasing) should be selected to obtain the principal component and perform life prediction.
FIG. 16A shows the flow of measurement and processing for realizing the method for converting the waveform data and the characteristic parameters described above into a random variable of normal distribution, the state determination method, and the state prediction method. Moreover, the circuit of the waveform data measurement and state determination apparatus for implement | achieving Fig.16 (a) is shown in FIG.
6). Detection method of singular component from pulsation signal and state determination method
Here, it demonstrates using the example of a gas engine misfire diagnosis.
FIGS. 10 (a), 11 (a), 12 (a), and 13 (a) are cylinder pressure waveform data when a gas engine misfire (abnormal portion of the peak value in the figure) occurs. FIGS. 12A and 13A show an example in which a large number of abnormal peaks occur continuously. FIG. 15 shows a flow of online singularity detection and state determination for such a local abnormality (single component) occurring in the pulsation signal. The procedure will be described below.
1). Preparation for diagnosis (FIG. 15A)
(1) The pulsation signal and rotation pulse signal to be diagnosed are simultaneously measured. The rotation pulse signal is also called a periodic pulse signal, and is used to determine the timing of each peak value of the pulsation signal.
(2) Perform low-pass filtering to remove noise. The cut-off frequency fL of the low-pass filter is determined as follows.
Here, n is the number of rotations of the shaft (rpm), z is the number of peaks per rotation (peak / 1 rotation), f0 is determined by observing the margin frequency (> n / 60, noise removal effect after filtering). ).
(3) Obtain envelope waveform data, peak waveform data, or moving average waveform data after low-pass filtering.
(4) Normalization of the envelope waveform data, peak waveform data, or moving average waveform data is performed as follows.
Here, x (t) is the normalized waveform data, x ′ (t) is the original waveform data, μ ′ (t) is the average value of x ′ (t), and s is the standard of x ′ (t). Deviation. Note that when μ ′ (t) and s are obtained, the waveform data in the normal state is desirable, but it may not be the waveform data in the normal state. Further, when obtaining μ ′ (t) and s, the operating conditions (the number of revolutions and the load) are set as much as possible in the determination. However, even if the operating conditions change slightly, the determination result is not greatly affected.
(5) The relative positional relationship between the peak of the rotation pulse signal and the peak of the pulsation signal is examined and determined, and a threshold value for identifying an abnormality is determined. The determination of the threshold is considered as follows.
At this time, it is determined that an abnormality has occurred. Here, k (= 2 to 4) is determined depending on the diagnostic object.
While | x (t) | is larger than kσ, it is determined that there is an abnormality at the peak of the pulsation signal corresponding to the rotation pulse signal.
2). Online diagnosis (Figure 15 (b))
(1) The pulsation signal and rotation pulse signal of the determination target are measured simultaneously.
(2) Perform low-pass filtering. The cut-off frequency fL of the low-pass filter is determined as shown in Expression (66).
(3) Obtain envelope waveform data or peak waveform data after low-pass filtering.
(4) Normalization of the above envelope waveform data or peak waveform data is obtained by the equation (68).
(5) Monitor whether the absolute value of the normalized envelope waveform data or peak waveform data is greater than a threshold value (kσ). If not, it is determined that there is no abnormality. If it is larger, it is determined that an abnormality has occurred. While | x (t) | is greater than kσ, the location of the abnormality is determined based on the relationship between the rotation pulse signal measured in advance and the peak of the pulsation signal.
Examples of signal processing based on the above procedure are shown in FIG. 10, FIG. 11, FIG. 12, and FIG.
10 (a) (b), FIG. 11 (a) (b), FIG. 12 (a) (b), FIG. 13 (a) (b) when a gas engine misfire (abnormal part in the figure) occurs It is time series waveform data and a spectrum of cylinder pressure.
10 (c) (d), FIG. 11 (c) (d), FIG. 12 (c) (d), and FIG. 13 (c) (d) are time-series waveform data and spectrum from which noise is removed by a low-pass filter. is there.
FIG. 10E, FIG. 11E, FIG. 12E, and FIG. 13E show rotation pulse waveform data.
FIG. 10F and FIG. 12F are envelope waveform data.
FIG. 11 (f) and FIG. 13 (f) show peak waveform data.
As can be seen from (f) of FIGS. 10, 11, 12 and 13, when an abnormality (missing peak or smaller) occurs, the absolute value | x (t) | is found to be larger than 2σ. it can.
Further, the abnormal location can be determined from the correspondence between the location where the absolute value | x (t) | is greater than 2σ and the rotation pulse waveform data.
Although the example of envelope waveform data and peak waveform data was shown above, as shown in FIG. 14, moving average waveform data (ui), The singular point of the pulsation signal can be detected using u (t) instead of x (t) as in the above procedure. The cutoff frequency (fc) Is obtained by equation (67), the moving average score (M) is obtained.
FIG. 16B shows the flow of measurement and processing for realizing the method for detecting and determining the singular component of the pulsation signal described above. FIG. 18 shows a circuit of a pulsation signal measurement and state determination apparatus for realizing FIG.
7). Process flow of status determination device or online status determination system
The processing flow of the state determination device or the online state determination system is shown in FIG. As shown in FIG. 16A, noise is removed from the measured waveform data, a feature parameter is obtained, the feature parameter is converted into a random variable of normal distribution, and then statistical testing, possibility theory, and feature parameter integration are performed. The presence / absence of state change is determined. These processes can be realized by a computer or a dedicated device. As shown in FIG. 16B, the singular component detection and state determination from the pulsation signal can be realized by a low-pass filter and the envelope (or peak value) by hardware. Since the numerical arithmetic unit (or the computer) has only to determine normalized x (t), | x (t)> kσ |, detect and determine a singular part, and display the determination result, it can be performed in real time.
図１は、測定した波形データを４種類、つまり、平均値より大きいデータｘｉ＋、平均値より小さいデータｘｉ−、式（１）による正規化後の絶対値データ｜ｘｉ｜及び全体波形データｘ＊ ｉＡに分ける例を示すグラフである。
図２は可能性分布関数の例を示すグラフである。
図３は正常状態と回転軸ミスアライメント状態の振動波形データの例を示すグラフである。
図４は正常状態の振動波形データに対する主成分分析結果を示すグラフである。
図５は回転軸ミスアライメント状態の振動波形データに対する主成分分析結果を示すグラフである。
図６は正常状態とアンバランス状態の振動波形データの例を示すグラフである。
図７は状態予測方法を示すグラフである。
図８は正常状態からアンバランス状態へ変化していく時の振動波形データの例を示すグラフである。
図９は各測定の波形の特徴パラメータ値を示すグラフである。
図１０は包絡線波形データによる１個のピークが欠落した場合の信号処理例を示すグラフである。
図１１はピーク波形データによる１個のピークが欠落した場合の信号処理例を示すグラフである。
図１２は包絡線波形データによる多数のピーク異常の場合の信号処理例を示すグラフである。
図１３はピーク波形データによる多数のピーク異常の場合の信号処理例を示すグラフである。
図１４は移動平均波形データの求め方を示すグラフである。
図１５は脈動信号の特異成分検出及び異常診断処理の流れを示すフローチャートである。
図１６は状態判定装置、または、オンライン状態判定システムの処理流れを示すフローチャートである。
図１７は信号計測及び状態判定装置の回路の一例を示す回路図であり、図中の符号は次の通りである。
１ センサ、２ アンプ、３ フィルタ、４ 信号処理・計算装置、５ 表示出力装置、６ データ用ＲＡＭ、７ ＡＤ変換器、８ ＤＣポート、９ ＳＣＩ（シリアル・コミュニケーション・インタフェース）、１０ ＣＰＵ、１１ フラッシュＲＯＭ、１２ 外部コンピュータ．
図１８は脈動信号計測及び状態判定装置の回路の一例を示す回路図であり、図中の符号は次の通りである。
１ センサ、２ アンプ、３ フィルタ、４ 包絡線処理部、５ 周期パルスセンサ、６ 表示出力装置、７ ＡＤ変換器、８ ＤＣポート、９ ＳＣＩ（シリアル・コミュニケーション・インタフェース）、１０ ＣＰＵ、１１ フラッシュＲＯＭ、１２ 外部コンピュータ、１３ データ用ＲＡＭ、１４ 信号処理・計算装置．FIG. 1 shows four types of measured waveform data, that is, data x i + larger than the average value, data x i− smaller than the average value, absolute value data | x i | after normalization by the equation (1), and the entire waveform It is a graph which shows the example divided into data x * iA .
FIG. 2 is a graph showing an example of the possibility distribution function.
FIG. 3 is a graph showing an example of vibration waveform data in a normal state and a rotational axis misalignment state.
FIG. 4 is a graph showing the principal component analysis results for vibration waveform data in a normal state.
FIG. 5 is a graph showing the principal component analysis results for the vibration waveform data in the rotational axis misalignment state.
FIG. 6 is a graph showing an example of vibration waveform data in a normal state and an unbalanced state.
FIG. 7 is a graph showing a state prediction method.
FIG. 8 is a graph showing an example of vibration waveform data when changing from a normal state to an unbalanced state.
FIG. 9 is a graph showing the characteristic parameter values of the waveform of each measurement.
FIG. 10 is a graph showing an example of signal processing when one peak is missing from the envelope waveform data.
FIG. 11 is a graph showing an example of signal processing when one peak is missing from the peak waveform data.
FIG. 12 is a graph showing an example of signal processing in the case of a large number of peak abnormalities based on envelope waveform data.
FIG. 13 is a graph showing an example of signal processing in the case of a number of peak abnormalities based on peak waveform data.
FIG. 14 is a graph showing how to obtain moving average waveform data.
FIG. 15 is a flowchart showing the flow of singular component detection and abnormality diagnosis processing of a pulsation signal.
FIG. 16 is a flowchart showing a processing flow of the state determination device or the online state determination system.
FIG. 17 is a circuit diagram showing an example of a circuit of the signal measurement and state determination apparatus, and the reference numerals in the figure are as follows.
1 sensor, 2 amplifier, 3 filter, 4 signal processing / calculation device, 5 display output device, 6 data RAM, 7 AD converter, 8 DC port, 9 SCI (serial communication interface), 10 CPU, 11 flash ROM, 12 External computer.
FIG. 18 is a circuit diagram showing an example of the circuit of the pulsation signal measurement and state determination apparatus, and the reference numerals in the figure are as follows.
1 sensor, 2 amplifier, 3 filter, 4 envelope processing unit, 5 period pulse sensor, 6 display output device, 7 AD converter, 8 DC port, 9 SCI (serial communication interface), 10 CPU, 11 flash ROM , 12 External computer, 13 Data RAM, 14 Signal processing / calculation device.
