(57)【要約】      （修正有）【課題】的確な基準を以って高感度の異常検知を可能とし周期運動体の劣化状況を的確に判定することができる状態監視方法を提供する。【解決手段】時系列データをＡＲＭＡモデルに当てはめ、パラメータの推定値と正常値との乖離度をシステムパラメータ間距離として簡易計算手法により求め、周期運動体の状態を監視する。【選択図】図４
2008020430.tif 000002
【特許請求の範囲】
【請求項１】
周期運動体の振動の変位を示すN個の信号からなる信号列を取得し、前記信号列中に正
常信号列よりも絶対値がS倍大きい大信号がm個の信号ごとに含まれるときに、ＡＲＭＡ（自己回帰移動平均）モデルを用いて、前記信号列のシステムパラメータ間距離J（但し，
【数１】
2008020430.tif 000181
正常時は、
【数２】
2008020430.tif 000182
とし、
系の異常が進展していった時のパラメータ推定値を
【数３】
2008020430.tif 000183
とし、ＡＲパラメータａ1，ａ2は，
【数４】
2008020430.tif 000184
一方、ＭＡパラメータｂ1，ｂ2については、
【数５】
2008020430.tif 000185
【数６】
2008020430.tif 000186
ARモデルの１次遅れ、２次遅れ、３次遅れ、４次遅れの自己相関係数ρ１、ρ２、ρ３、ρ４は、
【数７】
2008020430.tif 000187
【数８】
2008020430.tif 000188
【数９】
2008020430.tif 000189
【数１０】
2008020430.tif 000190
によって求められる。）を求め，
求めた距離Ｊの値に基づいて周期運動体の状態を監視することを特徴とする、周期運動体の状態監視方法。
【請求項２】
周期運動体は、転がり軸受け又は歯車からなる請求項１に記載の方法。
【請求項３】
周期運動体の振動の変位を示すN個の信号からなる信号列を取得し、前記信号列中に正
常信号列よりも絶対値がS倍大きい大信号がm個の信号ごとに含まれるときに、ＡＲＭＡ（自己回帰移動平均）モデルを用いて、前記信号列のシステムパラメータ間距離J（但し，
【数１１】
2008020430.tif 000191
正常時は、
【数１２】
2008020430.tif 000192
とし、
系の異常が進展していった時のパラメータ推定値を
【数１３】
2008020430.tif 000193
とし、ＡＲパラメータａ1，ａ2は，
【数１４】
2008020430.tif 000194
一方、ＭＡパラメータｂ1，ｂ2については、
【数１５】
2008020430.tif 000195
【数１６】
2008020430.tif 000196
ARモデルの１次遅れ、２次遅れ、３次遅れ、４次遅れの自己相関係数ρ１、ρ２、ρ３、ρ４は、
【数１７】
2008020430.tif 000197
【数１８】
2008020430.tif 000198
【数１９】
2008020430.tif 000199
【数２０】
2008020430.tif 000200
によって求められる。）を求め，
求めた距離Ｊの値に基づいて周期運動体の状態を監視する状態監視装置を備えた周期運動体の状態監視システム。
【請求項４】
状態監視装置は、システムパラメータ間距離Ｊが所定の基準値を超えているかどうかの判定を行う請求項３記載のシステム。
【請求項５】
表示部をさらに備え、前記状態監視装置は、パラメータ間距離Ｊ又は前記判定の結果を表示部に出力する請求項４記載のシステム。
【請求項６】
周期運動体の振動の変位を示すN個の信号からなる信号列を取得し、前記信号列中に正
常信号列よりも絶対値がS倍大きい大信号がm個の信号ごとに含まれるときに、ＡＲＭＡ（自己回帰移動平均）モデルを用いて、前記信号列のシステムパラメータ間距離J（但し，
【数２１】
2008020430.tif 000201
正常時は、
【数２２】
2008020430.tif 000202
とし、
系の異常が進展していった時のパラメータ推定値を
【数２３】
2008020430.tif 000203
とし、ＡＲパラメータａ1，ａ2は，
【数２４】
2008020430.tif 000204
一方、ＭＡパラメータｂ1，ｂ2については、
【数２５】
2008020430.tif 000205
【数２６】
2008020430.tif 000206
ARモデルの１次遅れ、２次遅れ、３次遅れ、４次遅れの自己相関係数ρ１、ρ２、ρ３、ρ４は、
【数２７】
2008020430.tif 000207
【数２８】
2008020430.tif 000208
【数２９】
2008020430.tif 000209
【数３０】
2008020430.tif 000210
によって求められる。）を求め，
求めた距離Ｊの値に基づいて周期運動体の状態を監視する工程をコンピュータに実行させることを特徴とする、周期運動体の状態監視プログラム。
【請求項７】
周期運動体の振動の変位を示すN個の信号からなる信号列を取得し、前記信号列中に正
常信号列よりも絶対値がS倍大きい大信号がm個の信号ごとに含まれるときに、ＡＲＭＡ（自己回帰移動平均）モデルを用いて、前記信号列のシステムパラメータ間距離J（但し，
【数３１】
2008020430.tif 000211
正常時は、
【数３２】
2008020430.tif 000212
とし、
系の異常が進展していった時のパラメータ推定値を
【数３３】
2008020430.tif 000213
とし、ＡＲパラメータａ1，ａ2，ａ3は，
【数３４】
2008020430.tif 000214
一方、ＭＡパラメータｂ1，ｂ2，ｂ3については、
【数３５】
2008020430.tif 000215
（ただし、
【数３６】
2008020430.tif 000216
ARモデルの１次遅れ、２次遅れ、３次遅れ、４次遅れ５次遅れ、６次遅れの自己相関係数ρ１、ρ２、ρ３、ρ４、ρ５、ρ６は、
【数３７】
2008020430.tif 000217
によって求められ、
Ｓ＝１が正常値として、
【数３８】
2008020430.tif 000218
また、
【数３９】
2008020430.tif 000219
となる。）を求め，
求めた距離Ｊの値に基づいて周期運動体の状態を監視することを特徴とする、周期運動体の状態監視方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、周期運動体の状態監視方法、状態監視システム及び状態監視プログラムに関する。本発明が対象とする周期運動体には、例えば、工場内の設備に用いられる転がり軸受けや歯車等があり、本発明はこのような周期運動体の状態が正常であるかどうか簡易に判断するための方怯、システム、プログラムを提供するものである。
【背景技術】
【０００２】
大型設備を有する鉄鋼業界等では，突発的に発生する設備故障でラインが停止すると，設備稼働率の低下，次工程への材料供給不足，納期が切迫している受注物件の納期遅れ等多大な損害が発生する。
【０００３】
そこで，これらの事態を防止するための設備異常検知は，重要な役割を果たす。従来は時間基準保全（Time Based Maintenance : TBM ）が主流をなしていたが、近年は設備監
視のハードウェア・ソフトウェアの性能アップも相まって、状態基準保全（ Condition Base d Maintenance : CBM ）が主流となってきた。こちらの方が、部品コスト低減、保全コスト低減、故障率低減につながるからである。保全を行なうと、保全後の初期故障を生じる確率が高くなるため、保全しなくても良いものも定期保全で保全したために初期故障を生じたりすることがある。部品の状態がよければできるだけ保全しない方がよいといえる。 CBM の方が求められる所以である。
【０００４】
CBM に移行してくると，異常の兆候をできるだけ速やかに捉えることが求められる。そのための手法として，現在様々なものが検討されているが，業態や分野によってその指標も異なってくる。
【０００５】
CBM においては常に設備を監視して，センサ類で収集したデーク群を処理することになる。そこではオンラインか擬似オンラインでの処理が必要である。異常検知手法としては従来感度の良い指標として Kurtosis 値、 Bicoherence 値(特許文献１、特許文献２)な
どのほか、時系列をAR（ Auto Regressive ；自己回帰）モデルに当てはめ，システムパ
ラメータを推定し，そのシステムパラメータの正常値からの季離度などを計算する方法（特許文献３、非特許文献１）やARMA( Auto Regressive Moving Average Model ；自己回帰移動平均）モデル（特許文献4、特許文献5、特許文献６）が用いられてきた。
【特許文献１】竹安数博：周期運動体の監視方法，特公昭 62- 60011 , ( 1987 )
【特許文献２】竹安数博：周期運動体の監視方法，特公昭 64-4611 , ( 1989 )
【特許文献３】竹安数博：周期運動体の監視方法，特公平 1-24246、(1989)
【特許文献４】竹安数博：パラメータ推定方法、状態監視方法、パラメータ推定装置、状態監視装置、状態監視システム、コンピュータプログラム、及び記録媒体、特開2004-295639、(2004)
【特許文献５】竹安数博：パラメータ推定方法、データ予測方法、パラメータ推定装置、データ予測装置、コンピュータプログラム、及び記録媒体、特開2004-302690、(2004)
【特許文献６】竹安数博：パラメータ計算方法、状態監視方法、パラメータ計算装置、状態監視装置、状態監視システム及びコンピュータプログラム、特開2004-295639、(2004)
【非特許文献１】K.Takeyasu, Y.Ishii "System Parameter Distance for Machine Diagnosis", IFORS, 2005
【発明の開示】
【発明が解決しようとする課題】
【０００６】
本発明は、時系列解析においては普遍的とも言えるARMAモデルによるシステムパラメータ推定時に簡易計算方法を導入し、現場でも簡便に設備診断ができる方法を導出・提示する。
ARMAモデルにおいては、そのシステム同定はARモデルと異なりYule-Walker方程式を解
いて解析的に求めることはできない。ブートストラップ法等を用いた繰り返し計算アルゴリズムによることになる。しかし、その方法の収束性は必ずしも良くなかった。
一方キュムラントを用いれば、3次白色の性質を利用して3次相関のYule-Walker方程式
で解析的に解くことができる。しかしながらMA部分は繰り返し計算となる。
本発明は、キュミュラントが4次白色であることに着目し、4変数の関係を用いれば、キュムラントも自己相関関数の組み合わせで解析的に計算できることを明らかにする。また数値計算で用いた２次のモデルではMA部分も解析的に導出されることを示す。
【課題を解決するための手段】
【０００７】
本発明の第１の解決手段は、離散時系列における（p、q）次のＡＲＭＡ（自己回帰移動平均）モデル
【数１】
2008020430.tif 000003
正常時は、
【数２】
2008020430.tif 000004
とし、
系の異常が進展していった時のパラメータ推定値を
【数３】
2008020430.tif 000005
とおき、
正常時のシステムパラメータと異常推定時のシステムパラメータ間距離Jをもって、監
視するシステムにおいて、周期運動体の振動の変位を示すN個の信号からなる信号列を取
得し、前記信号列中に正常信号列よりも絶対値がS倍大きい大信号がm個の信号ごとに含まれるときに、前記信号列のシステムパラメータ間距離Jを、４変数の関数用いて、
【数４】
2008020430.tif 000006
で求めて、周期運動体の状態を監視することを特徴とする、周期運動体の状態監視方法とするものである。
【０００８】
前項において、ARモデルの１次遅れ２次遅れ、３次遅れ、４次遅れの自己相関係数ρ１、ρ２、ρ３、ρ４は、
【数５】
2008020430.tif 000007
【数６】
2008020430.tif 000008
【数７】
2008020430.tif 000009
【数８】
2008020430.tif 000010
【数９】
2008020430.tif 000011
一方、ＭＡパラメータｂ１，ｂ２については、
【数１０】
2008020430.tif 000012
【数１１】
2008020430.tif 000013
【数１２】
2008020430.tif 000014
から得られた、Ｊ値を用いて
周期運動体の状態を監視することを特徴とする、周期運動体の状態監視方法に関するものである。
【０００９】
第２の解決手段は、これらの周期運動体の監視方法を転がり軸受け又は歯車に適応することとする。
【００１０】
第３の解決手段は、第１、第２の解決手段を周期運動体の状態監視装置を備えたシステムにすることである。
【００１１】
第４の解決手段は、第１、第２の解決手段を周期運動体の状態監視装置及び該状態の表示部を備え、監視する工程を備えたプログラムとすることにある。
更に、ｐ＝ｑ＝３の場合にでも、適用できるように、１次から６次の自己相関関数を用いたＪ値でもって、周期運動体の状態監視方法を提供する。
【発明の効果】
【００１２】
本発明では、現場でも効果的に活用できる設備診断の簡易手法として、時系列データを
ARMA モデルに当てはめ、パラメータの推定値と正常値との乖離度をシステムパラメータ間距離として、簡易計算手法により導き出した。
ARMA モデルにおいては、そのシステム同定は AR モデルと異なりYule-Walker方程式を
解いて解析的に求めることはできない。ブートストラップ法等を用いた繰り返し計算アルゴリズムによることになる。しかし、その方法の収束性は必ずしも良くなかった。
一方キュムラントを用いれば、 3 次白色の性質を利用して 3 次相関のYule-Walker方程式で解析的に解くことができる。しかしながら MA 部分は繰り返し計算となる。
本発明では、キュミュラントが 4 次白色であることに着目し、 4 変数の関係を用いれば、キュムラントも自己相関関数の組み合わせで解析的に計算できることを見出した。また数値計算で用いた 2 次のモデルでは MA 部分も解析的に導出されることを示した。そし
て結果の分析により手法の有効性を確認した。このため、周期運動体の劣化状況を的確に判定することができ、高感度の異常検知が可能となり、今後の設備診断の信頼性向上に寄与することができるものと考えられる。本手法は電卓を用いてもできる簡便なものであり、現場において効果的に活用することができる。また、マイコンチップ等に組み込むことによりリアルタイムに簡易に設備異常検知を行うことができる。
【発明を実施するための最良の形態】
【００１３】
時系列データを離散時間線形モデルへあてはめて分析するときは、 AR モデルや ARMA モデルを活用することが多い。本論文では時系列解析においては普遍的とも言える ARMA モデルを用いて検討する。
(ｐ、ｑ)次の ARMA（Ahtoregressive Moving Average Model；自己回帰移動平均）モデ
ルは一般に下記ように表記される。
【数１３】
2008020430.tif 000015
また、(1) 式における MA 過程は可逆条件を満たすものとする。
次に、システムパラメータ間距離を計算する。
現代制御理論において、系の状態空間表現を行い、 ARMA 過程表現に落とし込むと p≦q となることは周知である。  以下、簡単のため p = q とする。
【数１４】
2008020430.tif 000016
【数１５】
2008020430.tif 000017
【数１６】
2008020430.tif 000018
上記の J 値の推移をみることによって正常時から所定の乖離度となったとき系を異常と
判断する（図 1 ）。このようにして異常検知を行なうことができる。
【００１４】
基礎的な準備として、関係式を示す。
軸受、歯車等の回転体においては劣化が進行するに従って振動が大きくなる。また、据付等が不適切な場合にも振動が大きくなることは一般的によく知られている。計測対象から
2008020430.tif 000019
【数１７】
2008020430.tif 000020
分散σ２は
【数１８】
2008020430.tif 000021
で表される。
3 次モーメント MT ( 3)、 4 次モーメントMT(4)はそれぞれ次のように表される。
【数１９】
2008020430.tif 000022
【数２０】
2008020430.tif 000023
3 次モーメント、 4 次モーメントを正規化したものは、それぞれ歪度 Skewness (SK)、
尖度 Kurtosis(KT)として知られている。
【数２１】
2008020430.tif 000024
【数２２】
2008020430.tif 000025
これらを離散時間系で記述すると次のようになる。
サンプリング時間間隔をΔtとすると、離散時間データは
【数２３】
2008020430.tif 000026
【数２４】
2008020430.tif 000027
【数２５】
2008020430.tif 000028
【数２６】
2008020430.tif 000029
【数２７】
2008020430.tif 000030
【数２８】
2008020430.tif 000031
【数２９】
2008020430.tif 000032
となる。
【００１５】
一方、キュムラントについて、特性関数Ψ( u ）は次のように定義される。
【数３０】
2008020430.tif 000033
特性関数を Taylor展開し、各項、係数を c ( n ) とすると
【数３１】
2008020430.tif 000034
【数３２】
2008020430.tif 000035
となる。この c (n）をキュムラントと呼ぶ。キュムラントとモーメントとの間には次の
ような関係が存在する。
【数３３】
2008020430.tif 000036
【数３４】
2008020430.tif 000037
【数３５】
2008020430.tif 000038
【数３６】
2008020430.tif 000039
確率密度関数が正規分布のとき、 3 次以上の高次キュムラントは 0 になることが周知である。
【００１６】
２次の自己相関関数について、４つの変数に関する関数式
【数３７】
2008020430.tif 000040
の関係を用いて、x ( t ) , x ( t ＋τ1) , x ( t ＋τ2) , x ( t ＋τ3) について考
えると、
【数３８】
2008020430.tif 000041
【数３９】
2008020430.tif 000042
である。
ここで、τ2＝τ1、τ3＝０とおくと、
【数４０】
2008020430.tif 000043
とおくと、
【数４１】
2008020430.tif 000044
となる。これは相関関数の相関関数、つまり 2 次の相関関数ともみることができる。
【００１７】
４次キュムラントと自己相関関数を活用したARMAモデルのパラメータ推定について、説明すると、
まず、4次キュムラントに対するYule-Walker方程式を用いたARパラメータが推定について述べ、ついで自己相関関数を用いたARMAパラメータ推定アルゴリズムを示す。
ARパラメータ推定について
確率変数が正規分布に従うとき4次キュムラントは0となる。
2008020430.tif 000045
【数４２】
2008020430.tif 000046
（42）式を書き直すと
【数４３】
2008020430.tif 000047
を得る。
ここで
【数４４】
2008020430.tif 000048
の4次キュムラントである。
2008020430.tif 000049
【数４５】
2008020430.tif 000050
を得る。
これを
【数４６】
2008020430.tif 000051
【数４７】
2008020430.tif 000052
【数４８】
2008020430.tif 000053
【数４９】
2008020430.tif 000054
（1）2次の場合
【数５０】
2008020430.tif 000055
【００１８】
自己相関関数の簡易計算方法について、
振動信号を離散時間系で記述すると、信号データをサンプリングしたものとして、
2008020430.tif 000056
【数５１】
2008020430.tif 000057
【数５２】
2008020430.tif 000058
である。
自己相関関数は次のように定義される。
【数５３】
2008020430.tif 000059
そのため式(５３)は
【数５４】
2008020430.tif 000060
と書き変えることができる。Ｎが十分大きいとき0次の自己相関関数は分散にnearly equalである。自己相関係数は
【数５５】
2008020430.tif 000061
と表される。
回転体に傷がついた場合等には、回転周期ごとにピーク波形が生じる。サンプリングしたデータの回ごとに通常のＳ倍のピークをもつ信号が現れるものと仮定する。なお、サンプリング間隔の定め方についてはサンプリング定理にもとづく決定方法が周知であるが、ここでは、明確化するために単純化している。
ｍ回ごとに通常のＳ倍の信号が生ずると仮定する。また、特別なピーク（Ｓ倍の信号）時以外の平均、分散は通常時の平均、分散と同じであると仮定する。
2008020430.tif 000062
【数５６】
2008020430.tif 000063
【数５７】
2008020430.tif 000064
を得る。
【００１９】
2008020430.tif 000065
【数５８】
2008020430.tif 000066
【数５９】
2008020430.tif 000067
【数６０】
2008020430.tif 000068
【数６１】
2008020430.tif 000069
【数６２】
2008020430.tif 000070
【数６３】
2008020430.tif 000071
となる。<Case 1>と同様にして
【数６４】
2008020430.tif 000072
【表１】
2008020430.tif 000073
【００２０】
2008020430.tif 000074
【数６５】
2008020430.tif 000075
【数６６】
2008020430.tif 000076
【表２】
2008020430.tif 000077
【数６７】
2008020430.tif 000078
【数６８】
2008020430.tif 000079
【表３】
2008020430.tif 000080
【００２１】
2008020430.tif 000081
【数６９】
2008020430.tif 000082
【数７０】
2008020430.tif 000083
【表４】
2008020430.tif 000084
【００２２】
一般的にARMAモデルのMA部分は下記のように非線形方程式となる。
(1)式において
【数７１】
2008020430.tif 000085
とおくと
【数７２】
2008020430.tif 000086
【数７３】
2008020430.tif 000087
【数７４】
2008020430.tif 000088
【数７５】
2008020430.tif 000089
【数７６】
2008020430.tif 000090
【数７７】
2008020430.tif 000091
【数７８】
2008020430.tif 000092
【数７９】
2008020430.tif 000093
【数８０】
2008020430.tif 000094
【数８１】
2008020430.tif 000095
したがって
【数８２】
2008020430.tif 000096
【数８３】
2008020430.tif 000097
式（71）より
【数８４】
2008020430.tif 000098
【数８５】
2008020430.tif 000099
【数８６】
2008020430.tif 000100
【数８７】
2008020430.tif 000101
が得られる。
以上のように (2,2)次のARMAモデルを解析的に解くことができる。ARMAモデルにおいて１次から４次の自己相関関数を計算するだけで求めることができる。
【００２３】
ＡＲパラメータの推定について、
2008020430.tif 000102
【表５】
2008020430.tif 000103
【００２４】
MAパラメータの推定について、
2008020430.tif 000104
【数８８】
2008020430.tif 000105
したがって、式（55）より
【数８９】
2008020430.tif 000106
したがって
【数９０】
2008020430.tif 000107
【表６】
2008020430.tif 000108
【００２５】
いま、パラメータの正常値ai、biを、推定値ai∧、bi∧をとし、評価関数として前記Ｊを用いる。
【数９１】
2008020430.tif 000109
【数９２】
2008020430.tif 000110
したがって
【数９３】
2008020430.tif 000111
【数９４】
2008020430.tif 000112
の関係を利用して
【数９５】
2008020430.tif 000113
【数９６】
2008020430.tif 000114
【数９７】
2008020430.tif 000115
【数９８】
2008020430.tif 000116
となる。
【００２６】
2008020430.tif 000117
【表７】
2008020430.tif 000118
【００２７】
（２）3次の場合
式（４４）より、i＝0とすると
【数９９】
2008020430.tif 000119
（９９）
したがって
【数１００】
2008020430.tif 000120
（１００）
以上の関係を用いるとパラメータaの値は、下記のように表すことができる。
【数１０１】
2008020430.tif 000121
（１０１）
【００２８】
次に、5次遅れの自己相関係数について、記載する。
2008020430.tif 000122
【数１０２】
2008020430.tif 000123
（１０２）
2008020430.tif 000124
【数１０３】
2008020430.tif 000125
（１０３）
S＝2,4,6の場合の簡易数値計算を同様に行い、その値の変遷結果を表８に示す。
いま、前節の仮定の下で、m＝12，N＝100とすると結果は次表のようになる。
【表８】
2008020430.tif 000126
【００２９】
次に、6次遅れの自己相関係数について、記載する。
2008020430.tif 000127
【数１０４】
2008020430.tif 000128
（１０４）
2008020430.tif 000129
【数１０５】
2008020430.tif 000130
（１０５）
S＝2,4,6の場合の簡易数値計算を同様に行い、その値の変遷結果を表９に示す。
いま、前節の仮定の下で、m＝12，N＝100とすると結果は次表のようになる。
【表９】
2008020430.tif 000131
【００３０】
次に、ＡＲＭＡモデルのＭＡ部のパラメータ推定について、記載する。
一般的にARMAモデルのMA部分は下記のように非線形方程式となる。
（１） 式において
【数１０６】
2008020430.tif 000132
（１０６）
とおくと
【数１０７】
2008020430.tif 000133
（１０７）
2008020430.tif 000134
【数１０８】
2008020430.tif 000135
（１０８)
【数１０９】
2008020430.tif 000136
（１０９）
となることは周知である。
【００３１】
これは一般的には繰り返し演算法で解くなどのアルゴリズムが開発されている。
今回必要となるのはパラメータがb1b2のみであるので、下記のように解くことができる。
q＝3の場合
【数１１０】
2008020430.tif 000137
（１１０）
【数１１１】
2008020430.tif 000138
（１１１）
【数１１２】
2008020430.tif 000139
（１１２）
【数１１３】
2008020430.tif 000140
（１１３）
【数１１４】
2008020430.tif 000141
（１１４）
2008020430.tif 000142
【数１１５】
2008020430.tif 000143
（１１５）
【数１１６】
2008020430.tif 000144
（１１６）
【数１１７】
2008020430.tif 000145
（１１７）
【数１１８】
2008020430.tif 000146
（１１８）
【００３２】
したがって
【数１１９】
2008020430.tif 000147
（１１９）
【数１２０】
2008020430.tif 000148
（１２０）
【数１２１】
2008020430.tif 000149
（１２１）
ここで、式（１０６）より
【数１２２】
2008020430.tif 000150
（１２２）
【数１２３】
2008020430.tif 000151
（１２３）
【数１２４】
2008020430.tif 000152
（１２４）
【００３３】
これらの値を式（１１９）（１２０）（１２１）に代入して、b1,b2,b3の値を求めることができる。
以上のように (3,3)次のARMAモデルを解析的に解くことができる。ARMAモデルにおいて１次から6次の自己相関関数を計算するだけで求めることができる。
【００３４】
2008020430.tif 000153
【表１０】
2008020430.tif 000154
【００３５】
2008020430.tif 000155
定常エルゴード的正規過程より、平均０、分散１と仮定しても一般性を失わないため、
【数１２５】
2008020430.tif 000156
（１２５）
したがって、式（５５）より
【数１２６】
2008020430.tif 000157
（１２６）
したがって
【数１２７】
2008020430.tif 000158
（１２７）
2008020430.tif 000159
なお、正常時に解を持つ後述の式（139）（140）の場合のみ記載。
【表１１】
2008020430.tif 000160
【００３６】
次に、正常値との乖離度の計算について、記載する。
2008020430.tif 000161
p＝3であるから
【数１２８】
2008020430.tif 000162
（１２８）
上記のJ値の推移をみることによって正常時から所定の乖離度となったとき系を異常と判断する（図1）。このようにして異常検知を行なう。
【００３７】
（1）ARパラメータ
2008020430.tif 000163
【数１２９】
2008020430.tif 000164
(129)
【数１３０】
2008020430.tif 000165
(130)
【数１３１】
2008020430.tif 000166
(131)
S＝1のとき、式（64）、（66）、（68）より
【数１３２】
2008020430.tif 000167
（１３２）
のため
【数１３３】
2008020430.tif 000168
（１３３）
のとき
【数１３４】
2008020430.tif 000169
（１３４）
の関係を利用して
【数１３５】
2008020430.tif 000170
（１３５）
【数１３６】
2008020430.tif 000171
（１３６）
【数１３７】
2008020430.tif 000172
（１３７）
となる。
【００３８】
（2）ＭＡパラメータ
S＝1のとき正常値であると考えられる。
ここで、式（126）、（127）ならびに式（64）、（66）、（68）、（70）、(103)、(105)より
【数１３８】
2008020430.tif 000173
（１３８）
したがって、式（122）、（123）、(124)に式（135）、（136）、(137)のa1,a2,a3の値を代入し、上記と同様にして、式（119）（120）（121）からパラメータb1,b2,b3の正常値を求めることができる。
【００３９】
式（119）において、次式の場合分子、分母とも０となり、式（133）及び（134）の関係式を用いて、
【数１３９】
2008020430.tif 000174
（１３９）
の場合に解を持ち、そのときは、次式となる。
【数１４０】
2008020430.tif 000175
（１４０）
以上の式（139）（140）に式（133）（134）の関係式を用いて、
2008020430.tif 000176
【数１４１】
2008020430.tif 000177
（１４１）
【数１４２】
2008020430.tif 000178
（１４２）
【数１４３】
2008020430.tif 000179
（１４３）
となる。
【００４０】
評価関数の計算としては、表１０、表１１の値を用いて評価関数Ｊ値の計算を、式（128）を用いて行うと、次表１２のようになる。
【表１２】
2008020430.tif 000180
表１２および図３の評価関数Ｊの推移をみると、Ｓの値が大きくなるにしたがい、評価関数であるシステムパラメータ間距離Ｊは順次大きくなっていき、したがって異常検知の指標として用いることができる。
【００４１】
以下に、本件プログラムの実施の一例を示す。
図４は、本発明の状態監視システムの構成を示すブロック図である。本発明の状態監視システムは、転がり軸受け又は歯車等の周期運動体である工場内の設備の状態を監視し、設備に異常が生じた場合は警報を発するべく運用される。設備には、振動などの設備の状態を計測するセンサ 31 が設けられている。センサ 31 は、データ取得装置 32 に接続されており、計測データをデータ取得装置 32 へ入力すべく構成されている。データ取得装置 32 は、センサ 31 から入力された計測データを所定の周期でサンプリングし、複数の信号からなる信号列を作成し、作成した信号列から各種のデータを取得する機能を有している。データ取得装置 32 は、工場内に備えられた通信ネットワークＮＷに接続され、通信ネットワークＮＷは本発明の状態監視装置 1 に接続されており、データ取得装置 32 は、取得したデータを通信ネットワークＮＷを介して状態監視装置 1 へ送信する。
【００４２】
状態監視装置１は、コンピュータを用いて構成され、演算を行う CPU （演算部） 11 と、演算に伴って発生する一時的な情報を記憶する RAM （記憶部） 12 と、 CD−ROM ドライブ等の外部記憶装置 13 と、ハードディスク等の内部記憶装置 14 とを備えており、 CPU１１ は、 CD−ROM 等の本発明の記録媒体 2 から本発明のコンピュータプログラム 20 を外部記憶装置 13 にて読み取り、読み取ったコンピュータプログラム 20 を内部記憶装置 14 に記憶し、 RAM12 にコンピュータプログラム 20 をロードし、ロードしたコンピュータプログラム 20 に基づいて状態監視装置 1 に必要な処理を実行する。また、状態監視装置 1 は、工場内の通信ネットワーク NW に接続された入力部 15（受付部）を備えており、 CPU１１は、通信ネットワーク NW を介してデータ取得装置 32 から送信されたデータを入力部 15 にて受信する。更に、状態監視装置 1 は、情報を外部へ出力する出力部 16 を備えており、出力部 16 は、警報装置 4 に接続され、 CPU  11 は、設備の異常を示す情報を出力部 16 から警報装置 4 へ送信する。警報装置 4 は、ブザー、ランプ、又は警報の内容を表示する表示部などを備え、状態監視装置 1 から受信した情報に従って設備の異常を報知する。
【００４３】
なお、状態監視装置 1 は、通信ネットワーク NW に接続されている図示しない外部のサーバー装置から本発明に係るコンピュータプログラム 20 をダウンロードし、 CPU  11 にて処理を実行する形態であってもよい。
【００４４】
図５は、１実施形態に係る本発明の状態監視システムが行う動作を示すフローチャートである。センサ 31 は、設備の稼働に伴った図６に示す如き振動などのデータを計測し、データ取得装置 32 は、センサ 31 から入力された計測データをサンプリングし（ S1 01 ）、複数の信号からなる信号列を取得する。データ取得装置 32 は、サンプリングの結果取得した信号列に前記信号が所定の数 N 個蓄積されたか否かを判定し（ S102 ）、信号が N 個蓄積されていない場合は（ S102 : NO ）、ステップ S101 へ処理を戻してサンプリングを継続し、信号が N 個蓄積されている場合は（ S102 : YES ）、取得した信号列における信号の絶対値の平均の所定倍などの所定値よりも大きい絶対値を有する大信号が前記信号列に含まれているか否かを判定する（ S103 ）。大信号が信号列に含まれていない場合には（ S103 : NO ）、データ取得装置 32 は、設備は正常であるとして処理を終了し、大信号が信号列に含まれていた場合には（ S103 : YES ）、データ取得装置 32 は、取得した信号列にて、他の信号の絶対値に対する大信号の絶対値の倍率S 、及び一の大信号当たりに信号列に含まれる信号数 m を計測し（ S104 ）、取得した信号列に含まれる信号数 N 、倍率S 及び一の大信号当たりの信号数 m を、通信ネットワーク NW を介して状態監視装置 1 へ送信する（ S105 ）。
【００４５】
状態監視装置１の CPU  11 は、 N 、Ｓ 及び m を入力部 15 にて受信し（ Ｓ106 ）、コンピュータプログラム 20 を RAM12 へロードし、ロードしたコンピュータプログラム 20 に従って、受信した N 、Ｓ及び m を RAM12 に読み込み、監視対象の設備から得られた N 個の信号からなる信号列の 1次及び 2 次遅れの自己相関係数ρ１及びρ２を、数式 ( 61 ）及び（ 66 ）を用いて計算する（ Ｓ107 )。
【００４６】
CPU１１は、次に、 RAM12 にロードしたコンピュータプログラム 20 に従って、数式（50）を用いてパラメータ a 1 , a２ を計算し（ Ｓ108 ）、更に、数式(86),(87)を用いてb1、b2を計算し、数式(91)を用いて、システムパラメータ間距離 J を計算する（ Ｓ109 ) 。 CPU  １１ は、 RAM12 にロードしたコンピュータプログラム 20 に従って、計算したシステムパラメータ間距離 J を予め定めてある所定値と比較し、 J が所定値よりも大きいか否かを判定する（S110）。
【００４７】
ステップ S110 にて J が所定値よりも大きくない場合には（ S110 : NO ）、 CPU 11 は、 RAM12 にロードしたコンピュータプログラム 20 に従って、監視対象の設備は正常であると判定して、処理を終了する。 J が所定値より大きい場合には（ S110 : YES ）、 CPU11は、 RAM12 にロードしたコンピュータプログラム 20 に従って、監視対象の設備が異常であると判定して、 J の値に応じた異常の度合いを示す異常情報を、出力部 16 から警報装置 4 へ送信し（ S111 ）、処理を終了する。警報装置 4 は、状態監視装置 1 から受信した異常情報に従って、ブザーを鳴らす、ランプを点灯させる、又は表示部に異常情報の内容を表示する等、設備の異常を報知し、異常であると判定された設備が手動で停止されるなどの処置が行われる。
【００４８】
なお、前述の処理では、一の大信号当たりに信号列に含まれる信号数 m をデータ取得装置 32 にて計測する処理を含んでいるが、初期異常時に発生する一の大信号当たりの信号数 m は、転がり軸受けの転動体の数または歯車の数などの既知の値に対応しているため、 m の値を予め状態監視装置 1 の内部記憶装置 14 に記憶しておき、記憶している m の値を用いてシステムパラメータ間距離 J を計算する処理を用いてもよい。また、 CPU  １１は、 S110 を行わずに、 S109 で計算したシステムパラメータ間距離 J を表示部に出力し、表示部がこの値を表示するようにしてもよい。この場合、オペレータがこの値を確認して設備の状態を判断することができる。
【００４９】
以上詳述した如く、本発明においては、取得した信号列に大信号が含まれている場合に、信号列が含む信号数 N 、大信号の絶対値の倍率S及び一の大信号当たりに信号列に含まれる信号数 m を用いて、簡易的に 1 次及び 2 次遅れの自己相関係数ρ１及びρ２を計算し、計算したρ１及びρ２を用いてシステムパラメータ間距離 J を計算し、計算した J の値に設備の状態を判定する。この J の値は、周期運動体が正常な状態から離れるに従って大きくなるものであり、 J の値を参照することによって周期運動体の状態を適切に把握することができ、初期の段階で異常を発見することが可能になる。
【００５０】
本実施形態においては、一のセンサ 31 が計測したデータに基づいて設備の状態を監視する形態を示しているが、これに限るものではなく、一の設備、又は複数の設備の夫々に複数のセンサ 31 を備え、複数のセンサ 31 の夫々をデータ取得装置 32 及び通信ネットワーク NW を介して状態監視装置 1 に接続させ、夫々のセンサ 31 が計測した夫々 のデータに基づいて一又は複数の設備の状態を監視する形態としてもよい。また、本実施形態においては、データ取得装置 32 が取得した信号を蓄積する形態を示しているが、データ取得装置 32 は取得した信号を蓄積せずに状態監視装置 1 へ送信し、状態監視装置 1 は、受信した複数の信号を蓄積し、所定数の信号が蓄積されたときに、 N , S , m の値を RAM 12 に読み込んでシステムパラメータ間距離 J を計算する処理を行う形態としてもよい。 また、本実施形態においては、状態監視装置 1 の出力部 16 に警報装置 4 が接続されている形態を示しているが、状態監視装置 1 により設備が異常であると判定された場合に設備を停止させる等の設備の制御を行う制御装置を接続させる形態としてもよい。
【００５１】
また、本実施形態においては、監視対象の設備が転がり軸受けであるとして主に説明を行っているが、歯車またはチェーン等、周期的な運動を行う他の設備に対しても本発明は適用できる。
【００５２】
図５は、本発明の実施形態における状態監視方法を示す概念図である。本実施形態においては、電卓などの簡易的な計算装置を用いて設備の状態を把握する。設備に設けられたセンサ 31 には、オシロスコープ等のデータ表示装置 51 が接続されており、データ表示装置 51 は、センサ 31 が計測した計測データを表示する。データ表示装置 ５１には、図 ７ ( a ）、（ b ）に示す如きデータが表示され、設備に異常が発生した場合には、図７ ( b ）に示す如き大信号が含まれるデータが表示される。図６に示すように、設備の作業員などのオペレータは、データ表示装置 51 に表示されたデータから大信号の絶対値の倍率Sを目測し、信号数 N 、倍率S及び一の大信号当たりにデータに含まれる信号数ｍを電卓などの計算装置 52 へ入力し、数式（61）,（66）,(68)及び(70)を用いて簡易的に １次、２次、３次及び4次遅れの自己相関係数を計算し、数式（50）を用いて 2 次の AR モデルのパラメータであるa1、ａ２を計算し、また、数式(86),(87)を用いてb1、b2を計算し、数式(91)を用いて、システムパラメータ間距離 J を計算する。計算した J の値を以って設備の状態を判定する。
【００５３】
数式(61),(66),(68),(70)の簡易的な式を用いることにより、簡単に自己相関係数ρ１、ρ２ , ρ3及びρ4を計算することが可能となり、このため、簡単にシステムパラメータ間距離 J を計算することが可能なる。従って、センサ 31 が計測した計測データを直接に観測したオペレータが電卓などの簡易的な計算装置を用いてシステムパラメータ間距離 J を計算して、設備が稼働している現場において簡単に設備の状態を把握することができる。
【図面の簡単な説明】
【００５４】
【図１】本発明における評価関数による異常検知の方法の概念図である。
【図２】本発明による一実施例での評価指標値Ｊの推移である。
【図３】本発明他実施例での評価指標値Ｊの推移である。
【図４】本発明の状態監視システムの構成を示すプロック図である。
【図５】本発明一実施例に係る状態監視システムが行う動作を示すフローチャートである。
【図６】本発明の一実施形態における状態監視方法を示す概念図である。
【図７】周期振動体から得られるデータの例を示す模式的特性図である。
【符号の説明】
【００５５】
１     状態監視装置
１１   CPU （演算部）
１２   RAM （記憶部）
１５   入力部（受付部）
２     記録媒体
２０   コンピュータプログラム
３１   センサ
３２   データ取得装置
４     警報装置
５１   データ表示装置
５２   計算装置
【図１】
2008020430.tif 000220
【図２】
2008020430.tif 000221
【図３】
2008020430.tif 000222
【図４】
2008020430.tif 000223
【図５】
2008020430.tif 000224
【図６】
2008020430.tif 000225
【図７】
2008020430.tif 000226