JP7014080B2 - 情報処理装置、情報処理方法及びプログラム - Google Patents
情報処理装置、情報処理方法及びプログラム Download PDFInfo
- Publication number
- JP7014080B2 JP7014080B2 JP2018141513A JP2018141513A JP7014080B2 JP 7014080 B2 JP7014080 B2 JP 7014080B2 JP 2018141513 A JP2018141513 A JP 2018141513A JP 2018141513 A JP2018141513 A JP 2018141513A JP 7014080 B2 JP7014080 B2 JP 7014080B2
- Authority
- JP
- Japan
- Prior art keywords
- matrix
- eigenvalues
- transition data
- autocorrelation
- data
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Landscapes
- Testing Of Devices, Machine Parts, Or Other Structures Thereof (AREA)
- Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
- Complex Calculations (AREA)
Description
本発明は、情報処理装置、情報処理方法及びプログラムに関する。
周期的な運動を行う物体の異常診断を行う技術として、運動を行っている物体から得られる信号を測定し、その信号における振幅の時間的な推移に基づいて、その物体が正常であるか異常であるかの判定を行う方法がある。
特許文献１には、検出した信号から診断に必要な周波数帯域の信号を取り出し、取り出された信号のエンベロープ信号を求め、得られたエンベロープ信号の周波数帯域を低帯域化し、低帯域化した後のエンベロープ信号を間引きし、間引きした後のエンベロープ信号を周波数解析し、解析結果に基づいて異常を診断するシステムが開示されている。
特許文献１には、検出した信号から診断に必要な周波数帯域の信号を取り出し、取り出された信号のエンベロープ信号を求め、得られたエンベロープ信号の周波数帯域を低帯域化し、低帯域化した後のエンベロープ信号を間引きし、間引きした後のエンベロープ信号を周波数解析し、解析結果に基づいて異常を診断するシステムが開示されている。
特許文献１に開示された方法では、ローパスフィルタにかけることで低帯域化されたエンベロープ信号に対して周波数解析を行うことで、診断を行う。ローパスフィルタにおける遮断周波数は、ユーザによって定性的に決定される。しかし、ローパスフィルタを用いる場合、遮断周波数によっては、重要な信号が減衰し、その信号を検知できなくなる可能性がある。そのため、特許文献１では、診断の精度には限界があった。
本発明は、以上のような問題点に鑑みてなされたものであり、周期的な運動を行う物体の異常診断をより精度よく行うことができるようにすることを目的とする。
本発明は、以上のような問題点に鑑みてなされたものであり、周期的な運動を行う物体の異常診断をより精度よく行うことができるようにすることを目的とする。
本発明の情報処理装置は、周期的な運動を行う物体の前記周期的な運動に係る計測データに基づいて、前記計測データの振幅の時間的な推移の様子を示す推移データを取得する取得手段と、前記取得手段により取得された前記推移データに基づいて、修正された自己回帰モデルにおける係数を決定する決定手段と、前記決定手段により決定された前記係数に基づいて、前記物体の異常を診断する診断手段と、を有し、前記修正された自己回帰モデルは、前記推移データの実績値と、前記実績値に対する前記係数と、を用いて、前記推移データの予測値を表す式であり、前記決定手段は、前記推移データから導出される自己相関行列を特異値分解することで導出される前記自己相関行列の固有値を対角成分とする対角行列と、前記自己相関行列の固有ベクトルを列成分とする直交行列と、から導出される第１の行列を係数行列とし、前記推移データから導出される自己相関ベクトルを定数ベクトルとする方程式を用いて、前記係数を決定し、前記自己相関ベクトルは、時差が１から前記修正された自己回帰モデルで用いられる前記実績値の数であるｍまでの前記推移データの自己相関を成分とするベクトルであり、前記自己相関行列は、時差が０からｍ－１までの前記推移データの自己相関を成分とする行列であり、前記第１の行列は、前記自己相関行列の固有値のうち支配的な固有値として決定されたｓ個の固有値と前記対角行列とから導出される第２の行列Σsと、前記ｓ個の固有値と前記直交行列とから導出される第３の行列Ｕsと、から導出される行列ＵsΣsＵs
Tであり、前記第２の行列は、前記対角行列の部分行列であって、前記ｓ個の固有値を対角成分とする行列であり、前記第３の行列は、前記直交行列の部分行列であって、前記ｓ個の固有値に対応する固有ベクトルを列成分ベクトルとする行列である。
本発明によれば、周期的な運動を行う物体の異常診断をより精度よく行うことができる技術を提供することができる。
＜実施形態＞
以下、本発明の一実施形態について図面に基づいて説明する。
以下、本発明の一実施形態について図面に基づいて説明する。
（本実施形態の処理の概要）
本実施形態では、診断システムが鉄道台車に利用される軸受から測定された振動の信号に基づいて、その軸受の異常を診断する処理について説明する。
診断システムは、軸受の振動に応じた信号を測定し、測定された信号に基づいて、その信号における振幅の時間的な推移の様子を示す信号である推移データを取得し、取得した推移データを自己回帰モデルで表現し、自己回帰モデルの係数に関する条件式における行列を特異値分解し、その行列の固有値を求め、求めた固有値から支配的な固有値を決定し、決定した支配的な固有値を用いて、自己回帰モデルの係数を決定する。そして、診断システムは、決定した係数で特定される修正された自己回帰モデルにより推移データの予測値を求め、この推移データの予測値からなるデータのパワースペクトルに基づいて、軸受の異常を診断する。
本実施形態では、診断システムが鉄道台車に利用される軸受から測定された振動の信号に基づいて、その軸受の異常を診断する処理について説明する。
診断システムは、軸受の振動に応じた信号を測定し、測定された信号に基づいて、その信号における振幅の時間的な推移の様子を示す信号である推移データを取得し、取得した推移データを自己回帰モデルで表現し、自己回帰モデルの係数に関する条件式における行列を特異値分解し、その行列の固有値を求め、求めた固有値から支配的な固有値を決定し、決定した支配的な固有値を用いて、自己回帰モデルの係数を決定する。そして、診断システムは、決定した係数で特定される修正された自己回帰モデルにより推移データの予測値を求め、この推移データの予測値からなるデータのパワースペクトルに基づいて、軸受の異常を診断する。
（軸受異常の診断実験）
鉄道台車に利用されている軸受の異常診断精度の評価のために行った実験について説明する。
図１（ａ）は、本実験の状況を説明する図である。図１（ａ）の状況は、実験室に鉄道台車における車輪の駆動機構が設置されている状況である。駆動モータ１０１は、駆動モータ軸１０２、歯車継手１０３、小歯車軸１００を介して、小歯車１０５を回転させる。小歯車１０５は、回転することで、大歯車１０７を回転させる。それに合わせて、車軸１０９が回転し、車軸１０９に接続される車輪が回転することになる。本実験では、車軸１０９には、車輪の代わりにダイナモ１１０が接続されている。ダイナモ１１０は、疑似的な走行負荷を与えるために接続されている発電機である。本実験では、ダイナモ１１０を回転させるための負荷を、走行の際の負荷として想定する。
小歯車１０５が嵌められた小歯車軸１００と大歯車１０７が嵌められた車軸１０９とは、歯車箱１０４に固定されている。そして、歯車箱１０４には、小歯車軸との間に軸受１０６が装着されており、車軸１０９との間に軸受１０８が装着されている。図１（ｂ）は、歯車箱１０４に取り付けられている軸受を側面から見た様子を示す図である。
鉄道台車に利用されている軸受の異常診断精度の評価のために行った実験について説明する。
図１（ａ）は、本実験の状況を説明する図である。図１（ａ）の状況は、実験室に鉄道台車における車輪の駆動機構が設置されている状況である。駆動モータ１０１は、駆動モータ軸１０２、歯車継手１０３、小歯車軸１００を介して、小歯車１０５を回転させる。小歯車１０５は、回転することで、大歯車１０７を回転させる。それに合わせて、車軸１０９が回転し、車軸１０９に接続される車輪が回転することになる。本実験では、車軸１０９には、車輪の代わりにダイナモ１１０が接続されている。ダイナモ１１０は、疑似的な走行負荷を与えるために接続されている発電機である。本実験では、ダイナモ１１０を回転させるための負荷を、走行の際の負荷として想定する。
小歯車１０５が嵌められた小歯車軸１００と大歯車１０７が嵌められた車軸１０９とは、歯車箱１０４に固定されている。そして、歯車箱１０４には、小歯車軸との間に軸受１０６が装着されており、車軸１０９との間に軸受１０８が装着されている。図１（ｂ）は、歯車箱１０４に取り付けられている軸受を側面から見た様子を示す図である。
歯車箱１０４上の位置１１１は、振動計測装置が設置される位置である。振動計測装置は、加速度センサ、レーザ変位計等のセンサを含み、センサを介して物体の振動を検知し、検知した振動に応じた信号を出力する。位置１１１に設置された振動計測装置が位置１１１に伝達される振動を計測し、計測した振動を示す信号を外部の情報処理装置等に有線又は無線による通信を介して出力する。本実験では、振動計測装置は、位置１１１に設置されるとするが、軸受１０６に起因する振動が伝達される位置であるならば、任意の位置に設置することとしてもよい。
ここで、軸受１０６について、説明する。図２は、軸受１０６の一例の構成を説明する図である。軸受１０６は、外輪、内輪、転動体、保持器の４つの部分から構成される。図２には、軸受１０６の外輪、内輪、転動体、保持器の概要が示されている。軸受１０６は、外輪と内輪との間に保持器により保持された転動体が挟まれる構成となっている。
本実験では、小歯車１０５に装着されている軸受１０６として、疵のない正常な軸受、外輪に疵のある軸受のそれぞれを利用した。外輪に疵のある軸受については、他の種類の疵がないことも確認した。それぞれの軸受１０６毎に、駆動モータ１０１を駆動させ、位置１１１に設置された振動計測装置が振動を計測することで、軸受１０６の異常診断を行うための計測データを取得した。本実験では、駆動モータ１０１は、軸受１０６の内輪の回転数が３０００ｒｐｍ（ｒｏｕｎｄ ｐｅｒ ｍｉｎｕｔｅ）となるように駆動した。また、軸受１０６の外輪の疵の回転周期に対応する基本周波数は、約２６０Ｈｚである。また、軸受１０６の外輪につけた疵に係る信号は、この基本周波数の信号と、この基本周波数の信号の整数倍の信号（高調波）と、の重ね合された信号（例えば、矩形波等）となることが知られている。そのため、約２６０Ｈｚの信号と、約２６０Ｈｚの整数倍の信号と、が検知されると、軸受１０６の外輪に疵が存在すると診断できることとなる。
以上のように、軸受１０６が正常な場合、軸受１０６の外輪に疵がある場合について、位置１１１における振動の計測データを取得した。
ここで、軸受１０６について、説明する。図２は、軸受１０６の一例の構成を説明する図である。軸受１０６は、外輪、内輪、転動体、保持器の４つの部分から構成される。図２には、軸受１０６の外輪、内輪、転動体、保持器の概要が示されている。軸受１０６は、外輪と内輪との間に保持器により保持された転動体が挟まれる構成となっている。
本実験では、小歯車１０５に装着されている軸受１０６として、疵のない正常な軸受、外輪に疵のある軸受のそれぞれを利用した。外輪に疵のある軸受については、他の種類の疵がないことも確認した。それぞれの軸受１０６毎に、駆動モータ１０１を駆動させ、位置１１１に設置された振動計測装置が振動を計測することで、軸受１０６の異常診断を行うための計測データを取得した。本実験では、駆動モータ１０１は、軸受１０６の内輪の回転数が３０００ｒｐｍ（ｒｏｕｎｄ ｐｅｒ ｍｉｎｕｔｅ）となるように駆動した。また、軸受１０６の外輪の疵の回転周期に対応する基本周波数は、約２６０Ｈｚである。また、軸受１０６の外輪につけた疵に係る信号は、この基本周波数の信号と、この基本周波数の信号の整数倍の信号（高調波）と、の重ね合された信号（例えば、矩形波等）となることが知られている。そのため、約２６０Ｈｚの信号と、約２６０Ｈｚの整数倍の信号と、が検知されると、軸受１０６の外輪に疵が存在すると診断できることとなる。
以上のように、軸受１０６が正常な場合、軸受１０６の外輪に疵がある場合について、位置１１１における振動の計測データを取得した。
次に、振動計測装置が計測した計測データを用いた軸受１０６の異常診断の方法を説明する。以下では、振動計測装置が計測した計測データを、計測データｖとする。この計測データｖから、計測データｖにおける振幅の時間的な推移を示す推移データｙを求めることとした。本実験では、推移データｙとして、計測データｖのエンベロープを求めた。
ここで、計測信号ｖに基づいて、計測信号ｖのエンベロープを求める方法について説明する。
エンベロープとは、計測信号が振幅と角周波数が時刻によって変化する実信号（時刻ｔの実関数）であるとしたときに、この計測信号の解析信号の瞬時振幅として定義されるものである。計測信号ｖの解析信号ｖaとは、以下の式１で定義される複素信号（時刻ｔの複素関数）である。
ここで、計測信号ｖに基づいて、計測信号ｖのエンベロープを求める方法について説明する。
エンベロープとは、計測信号が振幅と角周波数が時刻によって変化する実信号（時刻ｔの実関数）であるとしたときに、この計測信号の解析信号の瞬時振幅として定義されるものである。計測信号ｖの解析信号ｖaとは、以下の式１で定義される複素信号（時刻ｔの複素関数）である。
ここで、ｉは、虚数単位であり、ｖ＾（式において＾はｖの上に付けて表記）は、以下の式２で定義される実信号であって、ｖのヒルベルト変換と呼ばれるものである。
ここで、πは、円周率である。そして、計測信号ｖのエンベロープｖmとは、以下の式３で定義される実信号である。
解析信号ｖaは、必ずしも式２の計算を行わなくともよく、フーリエ変換（ＦＦＴ）を用いて、以下の手順で比較的簡単に算出することができる。先ず、計測信号ｖのフーリエ変換Ｖを求める。次に、以下の式４を用いて、解析信号ｖaのフーリエ変換Ｖaを求める。最後に、Ｖaに逆フーリエ変換を施して解析信号ｖaを求める。
式３に示されるｖm（ｔ）は、ｖのエンベロープと呼ばれる信号であり、計測信号ｖにおける振幅の時間的な推移の様子を示す信号の一例である。
本実験では、計測データｖから上述のフーリエ変換を用いる方法で解析信号ｖaを算出し、式３を用いて、ｖのエンベロープ（ｖm（ｔ））を、推移データｙとして求めた。推移データには、エンベロープ以外にも以下のような信号がある。
エンベロープ（ｖm（ｔ））に正の実数を乗じたデータは、推移データの一例である。推移データは、時間的な値の推移の様子が示せればよいため、各時刻における値自体は重要でないためである。また、エンベロープ（ｖm（ｔ））の各時刻の値に、予め定められた値を足したデータについても、推移データの一例である。
また、例えば、ｖ（ｔ）における値がマイナスの部分について、－１を乗じて正の値にした関数ｖ’（ｔ）＝｜ｖ（ｔ）｜を想定し、このｖ’（ｔ）について、予め定められた閾値以上の周波数の信号をカットするようなローパスフィルタをかけることで求まる信号も、ｖの振幅の時間的な推移の様子を示す推移データの一例である。また、関数ｖ’（ｔ）における山の部分を特定し、特定した部分を補間（例えば、多項式補間）することで、求まる関数が示す信号も、ｖの振幅の時間的な推移の様子を示す推移データの一例である。
本実験では、計測データｖから上述のフーリエ変換を用いる方法で解析信号ｖaを算出し、式３を用いて、ｖのエンベロープ（ｖm（ｔ））を、推移データｙとして求めた。推移データには、エンベロープ以外にも以下のような信号がある。
エンベロープ（ｖm（ｔ））に正の実数を乗じたデータは、推移データの一例である。推移データは、時間的な値の推移の様子が示せればよいため、各時刻における値自体は重要でないためである。また、エンベロープ（ｖm（ｔ））の各時刻の値に、予め定められた値を足したデータについても、推移データの一例である。
また、例えば、ｖ（ｔ）における値がマイナスの部分について、－１を乗じて正の値にした関数ｖ’（ｔ）＝｜ｖ（ｔ）｜を想定し、このｖ’（ｔ）について、予め定められた閾値以上の周波数の信号をカットするようなローパスフィルタをかけることで求まる信号も、ｖの振幅の時間的な推移の様子を示す推移データの一例である。また、関数ｖ’（ｔ）における山の部分を特定し、特定した部分を補間（例えば、多項式補間）することで、求まる関数が示す信号も、ｖの振幅の時間的な推移の様子を示す推移データの一例である。
推移データｙは、計測データｖにおける比較的高い周波数の成分が除かれ、比較的低い周波数の成分の特徴が残っているデータとみなすことができる。そのため、軸受１０６に異常が生じた場合、推移データｙでは、例えば、軸受１０６に生じた異常に対応する計測データｖにおける比較的低い周波数の信号（例えば、疵の回転周期に係る信号等）が強調されることとなる。そこで、例えば、軸受１０６の異常により生じた信号が比較的高い周波数に現れず、疵の固有振動数と比較し、比較的低い周波数に現れる場合、推移データｙを用いることで、軸受１０６に生じた異常を診断できることとなる。
続いて、本実験では、推移データｙから、修正された自己回帰モデルにより推移データの予測値を求め、この推移データの予測値からなるデータから特徴量を算出することとした。
ある時刻ｋ（１≦ｋ≦Ｍ）における推移データｙの値をｙkとする。Ｍは、推移データｙがどの時刻までのデータを含むかを示す数であり、予め設定されている。ｙkを近似する自己回帰モデルは、例えば、以下の式５のようになる。式５に示すように、自己回帰モデルとは、時系列データにおけるある時刻ｋ（ｍ＋１≦ｋ≦Ｍ）のデータの予測値ｙ＾k（式において＾はｙの上に付けて表記）を、時系列データにおけるその時刻よりも前の時刻ｋ－ｌ（１≦ｌ≦ｍ）のデータの実績値ｙk-lを用いて表す式である。
続いて、本実験では、推移データｙから、修正された自己回帰モデルにより推移データの予測値を求め、この推移データの予測値からなるデータから特徴量を算出することとした。
ある時刻ｋ（１≦ｋ≦Ｍ）における推移データｙの値をｙkとする。Ｍは、推移データｙがどの時刻までのデータを含むかを示す数であり、予め設定されている。ｙkを近似する自己回帰モデルは、例えば、以下の式５のようになる。式５に示すように、自己回帰モデルとは、時系列データにおけるある時刻ｋ（ｍ＋１≦ｋ≦Ｍ）のデータの予測値ｙ＾k（式において＾はｙの上に付けて表記）を、時系列データにおけるその時刻よりも前の時刻ｋ－ｌ（１≦ｌ≦ｍ）のデータの実績値ｙk-lを用いて表す式である。
式５におけるαは、自己回帰モデルの係数である。また、ｍは、自己回帰モデルにおいてある時刻ｋに対応する推移データｙの値であるｙkを、その時刻よりも前の過去幾つのデータを用いて近似するかを示すＭ未満の整数であり、本実験では、５００とする。
続いて、最小二乗法を用いて、自己回帰モデルによる予測値ｙ＾kが実測値であるｙkに近似するための条件式を求める。自己回帰モデルによる予測値ｙ＾kが実測値であるｙkに近似するための条件として、式５で与えられるｙ＾kとｙkとの二乗誤差を最小化することを考える。即ち、本実験では、自己回帰モデルによる予測値ｙ＾kをｙkに近似するために最小二乗法を用いる。以下の式６は、推移データと自己回帰モデルによる予測値との二乗誤差を最小にするための条件式である。
続いて、最小二乗法を用いて、自己回帰モデルによる予測値ｙ＾kが実測値であるｙkに近似するための条件式を求める。自己回帰モデルによる予測値ｙ＾kが実測値であるｙkに近似するための条件として、式５で与えられるｙ＾kとｙkとの二乗誤差を最小化することを考える。即ち、本実験では、自己回帰モデルによる予測値ｙ＾kをｙkに近似するために最小二乗法を用いる。以下の式６は、推移データと自己回帰モデルによる予測値との二乗誤差を最小にするための条件式である。
式６より、以下の式７の関係を満たす。
また、式７を変形（行列表記）することで、以下の式８のようになる。
式８におけるＲjlは推移データｙの自己相関と呼ばれるもので、以下の式９で定義される値である。このときの｜ｊ－ｌ｜を時差という。
式８を基に、以下の自己回帰モデルの係数に関する関係式である式１０を考える。式１０は、自己回帰モデルによる推移データの予測値と、その予測値に対応する時刻における推移データと、の誤差を最小化する条件から導出される方程式で、ユール・ウォーカー（Ｙｕｌｅ－Ｗａｌｋｅｒ）方程式と呼ばれるものである。また、式１０は自己回帰モデルの係数から成るベクトルを変数ベクトルとする線形方程式で、式１０における左辺の定数ベクトルは、時差が１からｍまでの推移データの自己相関を成分とするベクトルで、以下では、自己相関ベクトルとする。また、式１０における右辺の係数行列は、時差が０からｍ－１までの推移データの自己相関を成分とする行列であり、以下では、自己相関行列とする。
また、式１０における右辺の自己相関行列（Ｒjlで構成されるｍ×ｍの行列）を、以下の式１１のように、自己相関行列Ｒと表記する。
一般に、自己回帰モデルの係数を求める際には、式１０を係数αについて解くという方法が用いられる。式１０では、自己回帰モデルで導出されるある時刻ｋにおける推移データの予測値ｙ＾kが、その時刻ｋにおける推移データの実績値ｙkにできるだけ近づくように係数αを導出する。よって、自己回帰モデルの周波数特性に、各時刻における推移データの実績値ｙkに含まれる多数の周波数成分が含まれる。したがって、例えば、推移データｙに含まれるノイズが多い場合には、軸受１０６の振動に係る信号を抽出できなかったり、軸受１０６の故障の態様に応じた特徴を抽出することができなかったりするという問題が生じる。
そこで、本発明者らは、自己回帰モデルの係数αに乗算される自己相関行列Ｒに着目し、鋭意検討した結果、自己相関行列Ｒの固有値の一部を用いて、推移データｙに含まれるノイズの影響が低減され、軸受に生じた異常（例えば、疵、汚れ等）に係る信号成分が強調される（ＳＮ比を高める）ように自己相関行列Ｒを書き換えればよいという着想に至った。
以下で、このことの具体例を説明する。
自己相関行列Ｒを特異値分解する。自己相関行列Ｒの成分は、対称であるので、自己相関行列Ｒを特異値分解すると以下の式１２のように、直交行列Ｕと、対角行列Σと、直交行列Ｕの転置行列との積となる。
そこで、本発明者らは、自己回帰モデルの係数αに乗算される自己相関行列Ｒに着目し、鋭意検討した結果、自己相関行列Ｒの固有値の一部を用いて、推移データｙに含まれるノイズの影響が低減され、軸受に生じた異常（例えば、疵、汚れ等）に係る信号成分が強調される（ＳＮ比を高める）ように自己相関行列Ｒを書き換えればよいという着想に至った。
以下で、このことの具体例を説明する。
自己相関行列Ｒを特異値分解する。自己相関行列Ｒの成分は、対称であるので、自己相関行列Ｒを特異値分解すると以下の式１２のように、直交行列Ｕと、対角行列Σと、直交行列Ｕの転置行列との積となる。
式１２の行列Σは、式１３に示すように、対角成分が自己相関行列Ｒの固有値となる対角行列である。対角行列Σの対角成分を、σ11、σ22、・・・、σmmとする。また、行列Ｕは、各列成分ベクトルが自己相関行列Ｒの固有ベクトルとなる直交行列である。直交行列Ｕの列成分ベクトルを、ｕ1、ｕ2、・・・、ｕmとする。自己相関行列Ｒの固有ベクトルｕjに対する固有値がσjjという対応関係が有る。自己相関行列Ｒの固有値は、自己回帰モデルによる推移データの予測値の時間波形に含まれる各周波数の成分の強度に反映する変数である。
図３は、軸受１０６に外輪に疵がある軸受を用いた場合に計測された計測データから得られた推移データについて、上記手順で求められた自己相関行列Ｒの固有値の分布を示す図である。
図３のグラフは、自己相関行列Ｒを特異値分解して得られた固有値σ11～σmmを昇順に並べ替え、プロットしたグラフであり、横軸が固有値のインデックス、縦軸が固有値の値を示す。
図３を見ると、他よりも顕著に高い値をもつ固有値が１つあるのが分かる。また、その他の固有値の値を見てみると２つの固有値が、他の固有値よりも比較的高い値を持つことが分かる。
図３のグラフは、自己相関行列Ｒを特異値分解して得られた固有値σ11～σmmを昇順に並べ替え、プロットしたグラフであり、横軸が固有値のインデックス、縦軸が固有値の値を示す。
図３を見ると、他よりも顕著に高い値をもつ固有値が１つあるのが分かる。また、その他の固有値の値を見てみると２つの固有値が、他の固有値よりも比較的高い値を持つことが分かる。
また、軸受１０６に疵がない正常な軸受を用いた場合に計測された計測データから得られた推移データについて、上記手順で求められた自己相関行列Ｒの固有値の分布についても、図３と同様に、他よりも顕著に高い値をもつ固有値が１つあり、その他の固有値のうち、幾つかの固有値が他よりも比較的高い値を持つような分布となった。
何れの場合でも、自己相関行列Ｒの固有値のうち値が他よりも顕著に高い固有値の数は、固有値全体の数よりも顕著に少ないことが分かる。そのため、固有値全体の数よりも顕著に少ない個数の、他よりも値が顕著に高い固有値が、自己相関行列Ｒの支配的な成分に対応する固有値であると判断できる。行列の支配的な成分とは、その行列における大部分を占める成分である。そこで、発明者らは、全ての固有値のうち、支配的な固有値だけを用いて、他の固有値を用いずに、自己相関行列Ｒを修正（行列Ｒ’を生成）することで、支配的な成分を残したまま、余計なノイズ等の成分を除外することができるという着想を得た。行列の支配的な固有値とは、その行列の支配的な成分に対応する固有値である。
何れの場合でも、自己相関行列Ｒの固有値のうち値が他よりも顕著に高い固有値の数は、固有値全体の数よりも顕著に少ないことが分かる。そのため、固有値全体の数よりも顕著に少ない個数の、他よりも値が顕著に高い固有値が、自己相関行列Ｒの支配的な成分に対応する固有値であると判断できる。行列の支配的な成分とは、その行列における大部分を占める成分である。そこで、発明者らは、全ての固有値のうち、支配的な固有値だけを用いて、他の固有値を用いずに、自己相関行列Ｒを修正（行列Ｒ’を生成）することで、支配的な成分を残したまま、余計なノイズ等の成分を除外することができるという着想を得た。行列の支配的な固有値とは、その行列の支配的な成分に対応する固有値である。
そこで、行列Ｒの固有値σ11、σ22、・・・、σmmを求めた後、σ11、σ22、・・・、σmmのうち、支配的な固有値を決定した。
本実験では、σ11、σ22、・・・、σmmの全ての平均を求めて、求めた平均よりも大きい固有値の全てを、支配的な固有値として決定した。本実験では、疵のない正常な軸受１０６を用いた場合、外輪に疵のある軸受１０６を用いた場合の双方において、σ11、σ22、σ33の３個の固有値を、支配的な固有値として決定した。即ち、使用する固有値の数である使用固有値数ｓの値は、３となる。
また、支配的な固有値を、以下のようにして決定することもできる。例えば、選択した固有値の合計のσ11、σ22、・・・、σmmの総和に対する割合が予め定められた閾値（例えば、９５パーセント等）以上となるように、大きいものから選択した固有値を、支配的な固有値として決定できる。また、選択した固有値の合計のσ11、σ22、・・・、σmmの総和に対する割合が予め定められた閾値（例えば、７５パーセント等）以上となるように、選択した予め定められた閾値（例えば、１個、３個、５個、固有値全体の個数の１パーセント等）以下の個数の固有値を、支配的な固有値として決定することもできる。
自己相関行列Ｒの特異値分解の結果得られる対角行列Σの対角成分であるσ11、σ22、・・・、σmmの値は、数式の表記を簡略にするために降順とする。これらの自己相関行列Ｒの固有値のうち、支配的な固有値として決定したｓ個の固有値を用いて、以下の式１４のように、行列Ｒ’を定義する。行列Ｒ’は、自己相関行列Ｒの固有値のうち使用固有値数ｓ個の固有値用いて自己相関行列Ｒを近似した行列である。
本実験では、σ11、σ22、・・・、σmmの全ての平均を求めて、求めた平均よりも大きい固有値の全てを、支配的な固有値として決定した。本実験では、疵のない正常な軸受１０６を用いた場合、外輪に疵のある軸受１０６を用いた場合の双方において、σ11、σ22、σ33の３個の固有値を、支配的な固有値として決定した。即ち、使用する固有値の数である使用固有値数ｓの値は、３となる。
また、支配的な固有値を、以下のようにして決定することもできる。例えば、選択した固有値の合計のσ11、σ22、・・・、σmmの総和に対する割合が予め定められた閾値（例えば、９５パーセント等）以上となるように、大きいものから選択した固有値を、支配的な固有値として決定できる。また、選択した固有値の合計のσ11、σ22、・・・、σmmの総和に対する割合が予め定められた閾値（例えば、７５パーセント等）以上となるように、選択した予め定められた閾値（例えば、１個、３個、５個、固有値全体の個数の１パーセント等）以下の個数の固有値を、支配的な固有値として決定することもできる。
自己相関行列Ｒの特異値分解の結果得られる対角行列Σの対角成分であるσ11、σ22、・・・、σmmの値は、数式の表記を簡略にするために降順とする。これらの自己相関行列Ｒの固有値のうち、支配的な固有値として決定したｓ個の固有値を用いて、以下の式１４のように、行列Ｒ’を定義する。行列Ｒ’は、自己相関行列Ｒの固有値のうち使用固有値数ｓ個の固有値用いて自己相関行列Ｒを近似した行列である。
式１４における行列Ｕsは、式１２の直交行列Ｕの左からｓ個の列成分ベクトル（使用される固有値に対応する固有ベクトル）により構成されるｍ×ｓ行列である。つまり、行列Ｕsは、直交行列Ｕから左のｍ×ｓの成分を切り出して構成される部分行列である。また、Ｕs
TはＵsの転置行列であり、式１２の行列ＵTの上からｓ個の行成分ベクトルにより構成されるｓ×ｍ行列である。式１４における行列Σsは、式１２の対角行列Σの左からｓ個の列と上からｓ個の行により構成されるｓ×ｓ行列である。つまり、行列Σsは、対角行列Σから左上のｓ×ｓの成分を切り出して構成される部分行列である。
行列Σs及び行列Ｕsを行列成分表現すれば、以下の式１５のようになる。
行列Σs及び行列Ｕsを行列成分表現すれば、以下の式１５のようになる。
自己相関行列Ｒの代わりに行列Ｒ’を用いることで、式１０の関係式を、以下の式１６のように書き換える。
式１６を変形することで、係数αを求める式１７が得られる。式１７によって求められた係数αを用いて、式５により予測値ｙ＾kを算出するモデルを「修正された自己回帰モデル」とする。
これまで対角行列Σの対角成分であるσ11、σ22、・・・、σmmの値を降順として説明したが、係数αの算出過程において対角行列Σの対角成分は降順である必要はなく、その場合には、行列Ｕsは直交行列Ｕから左のｍ×ｓの成分を切り出すのではなく、使用される固有値に対応する列成分ベクトル（固有ベクトル）を切り出して構成される部分行列であり、行列Σsは対角行列Σから左上のｓ×ｓの成分を切り出すのではなく、使用される固有値を対角成分とするように切り出される部分行列である。
式１７は、修正された自己回帰モデルの係数の決定に利用される方程式である。式１７の行列Ｕsは、自己相関行列Ｒの特異値分解により得られる直交行列の部分行列であって、利用される固有値に対応する固有ベクトルを列成分ベクトルとする行列である第３の行列の一例である。また、式１７の行列Σsは、自己相関行列Ｒの特異値分解により得られる対角行列の部分行列であって、利用される固有値を対角成分とする行列である第２の行列の一例である。そして、式１７の行列ＵsΣsＵs Tは、行列Σsと行列Ｕsとから導出される行列である第１の行列の一例である。
これまで対角行列Σの対角成分であるσ11、σ22、・・・、σmmの値を降順として説明したが、係数αの算出過程において対角行列Σの対角成分は降順である必要はなく、その場合には、行列Ｕsは直交行列Ｕから左のｍ×ｓの成分を切り出すのではなく、使用される固有値に対応する列成分ベクトル（固有ベクトル）を切り出して構成される部分行列であり、行列Σsは対角行列Σから左上のｓ×ｓの成分を切り出すのではなく、使用される固有値を対角成分とするように切り出される部分行列である。
式１７は、修正された自己回帰モデルの係数の決定に利用される方程式である。式１７の行列Ｕsは、自己相関行列Ｒの特異値分解により得られる直交行列の部分行列であって、利用される固有値に対応する固有ベクトルを列成分ベクトルとする行列である第３の行列の一例である。また、式１７の行列Σsは、自己相関行列Ｒの特異値分解により得られる対角行列の部分行列であって、利用される固有値を対角成分とする行列である第２の行列の一例である。そして、式１７の行列ＵsΣsＵs Tは、行列Σsと行列Ｕsとから導出される行列である第１の行列の一例である。
式１７の右辺を計算することにより、修正された自己回帰モデルの係数αが求まる。以上、修正された自己回帰モデルの係数の導出方法の一例について説明してきたが、その基となる自己回帰モデルの係数の導出については直感的に分かり易いように予測値に対して最小二乗法を用いる方法で説明した。しかしながら、一般的には確率過程という概念を用いて自己回帰モデルを定義し、その係数を導出する方法が知られている。その場合に、自己相関は確率過程（母集団）の自己相関で表現されており、この確率過程の自己相関は時差の関数として表されるものである。従って、本実施形態における計測データの自己相関は、確率過程の自己相関を近似するものであれば他の計算式で算出した値に代えても良く、例えば、Ｒ22～Ｒmmは時差が０の自己相関であるが、これらをＲ11に置き換えてもよい。
次に、求めた係数αで特定される修正された自己回帰モデルにより推移データの予測値を求め、この推移データの予測値からなるデータについて、フーリエ変換を施すことで、この推移データの予測値からなるデータのパワースペクトルを求める。
次に、求めた係数αで特定される修正された自己回帰モデルにより推移データの予測値を求め、この推移データの予測値からなるデータについて、フーリエ変換を施すことで、この推移データの予測値からなるデータのパワースペクトルを求める。
本実験では、得られた推移データ毎に、式９と式１１を用いて自己相関行列Ｒを求め、式１２で表される特異値分解を行うことによって自己相関行列Ｒの固有値を求め、自己相関行列Ｒの固有値の分布を求めた。また、得られた推移データ毎に、推移データの波形を求めた。また、得られた推移データ毎に、上記手順で求められた自己相関行列Ｒの特異値分解から式１５を用いて行列Σsと行列Ｕsとを求め、式９と式１７を用いて修正された自己回帰モデルの係数αを求め、求めた係数αで特定される修正された自己回帰モデルの波形を求めた。そして、得られた推移データ毎に、求めた修正された自己回帰モデルの波形について、フーリエ変換を施すことで、この波形のパワースペクトルを求めた。
また、比較の対象として、疵のない正常な軸受１０６を用いた場合、外輪に疵のある軸受１０６を用いた場合のそれぞれについて、計測データｖの解析信号ｖaから式３を用いて求められた推移データｙ（エンベロープ）それぞれについて、フーリエ変換を施すことで、推移データｙのパワースペクトルを求めた。
また、比較の対象として、疵のない正常な軸受１０６を用いた場合、外輪に疵のある軸受１０６を用いた場合のそれぞれについて、計測データｖの解析信号ｖaから式３を用いて求められた推移データｙ（エンベロープ）に対して、ローパスフィルタをかけて求まる信号について、フーリエ変換を施すことで、このローパスフィルタをかけて求まる信号のパワースペクトルを求めた。ここで、推移データｙに、１５００Ｈｚを遮断周波数としたローパスフィルタをかけた場合、８００Ｈｚを遮断周波数としたローパスフィルタをかけた場合、６００Ｈｚを遮断周波数としたローパスフィルタをかけた場合それぞれについて、パワースペクトルを求めた。
また、比較の対象として、疵のない正常な軸受１０６を用いた場合、外輪に疵のある軸受１０６を用いた場合のそれぞれについて、計測データｖの解析信号ｖaから式３を用いて求められた推移データｙ（エンベロープ）それぞれについて、フーリエ変換を施すことで、推移データｙのパワースペクトルを求めた。
また、比較の対象として、疵のない正常な軸受１０６を用いた場合、外輪に疵のある軸受１０６を用いた場合のそれぞれについて、計測データｖの解析信号ｖaから式３を用いて求められた推移データｙ（エンベロープ）に対して、ローパスフィルタをかけて求まる信号について、フーリエ変換を施すことで、このローパスフィルタをかけて求まる信号のパワースペクトルを求めた。ここで、推移データｙに、１５００Ｈｚを遮断周波数としたローパスフィルタをかけた場合、８００Ｈｚを遮断周波数としたローパスフィルタをかけた場合、６００Ｈｚを遮断周波数としたローパスフィルタをかけた場合それぞれについて、パワースペクトルを求めた。
求められた実験結果を、図４～９を用いて説明する。
図４は、疵のない正常な軸受１０６を用いた場合、外輪に疵のある軸受１０６を用いた場合のそれぞれについて計測された計測データｖを示す図である。
図４（ａ）のグラフは、疵のない正常な軸受１０６を用いた場合に計測された計測データｖの波形を示す。図４（ａ）のグラフの横軸は、時間を示し、縦軸は、振幅を示す。図４（ｂ）のグラフは、外輪に疵のある軸受１０６を用いた場合に計測された計測データｖの波形を示す。図４（ｂ）のグラフの横軸は、時間を示し、縦軸は、振幅を示す。
図４（ａ）のグラフと図４（ｂ）のグラフとを見ると、似通った波形をしている。そのため、これらの波形から、軸受１０６が正常であるか、異常であるかを診断することが困難である。
図４は、疵のない正常な軸受１０６を用いた場合、外輪に疵のある軸受１０６を用いた場合のそれぞれについて計測された計測データｖを示す図である。
図４（ａ）のグラフは、疵のない正常な軸受１０６を用いた場合に計測された計測データｖの波形を示す。図４（ａ）のグラフの横軸は、時間を示し、縦軸は、振幅を示す。図４（ｂ）のグラフは、外輪に疵のある軸受１０６を用いた場合に計測された計測データｖの波形を示す。図４（ｂ）のグラフの横軸は、時間を示し、縦軸は、振幅を示す。
図４（ａ）のグラフと図４（ｂ）のグラフとを見ると、似通った波形をしている。そのため、これらの波形から、軸受１０６が正常であるか、異常であるかを診断することが困難である。
図５は、疵のない正常な軸受１０６を用いた場合に、計測データｖの解析信号ｖaから式３を用いて求められた推移データｙについてのパワースペクトルを示す図である。図５のパワースペクトルは、０Ｈｚ～１６００Ｈｚまでの周波数の信号のパワーを示すパワースペクトルである。パワースペクトルの横軸は、周波数を示し、縦軸は、対応する周波数の信号の強度を示す。
図６は、疵のない正常な軸受１０６を用いた場合に求められた推移データｙから、上記手順で求められた係数αで特定される修正された自己回帰モデルにより求められる推移データの予測値からなるデータについてのパワースペクトルを示す図である。図６のパワースペクトルは、０Ｈｚ～１６００Ｈｚまでの周波数の信号のパワーを示すパワースペクトルである。
図５と図６とのパワースペクトルを見てみると、図６のパワースペクトルでは、図５のパワースペクトルと比べて約８００Ｈｚ以上の周波数の信号が顕著に減衰していることが分かる。図６のパワースペクトルでは、推移データｙについての自己相関行列の固有値のうち、支配的なものだけを用いて求められた修正された自己回帰モデルについてのパワースペクトルである。そのため、推移データｙについて、本質的な成分は、０Ｈｚ～約８００Ｈｚの信号であると判断できる。そこで、軸受１０６に何らかの異常が生じて、推移データｙがその異常の影響を受ける場合、その影響は、推移データｙにおける本質的な成分である０Ｈｚ～約８００Ｈｚの信号（修正された自己回帰モデルにおける本質的な成分）に表れると判断できる。
図６は、疵のない正常な軸受１０６を用いた場合に求められた推移データｙから、上記手順で求められた係数αで特定される修正された自己回帰モデルにより求められる推移データの予測値からなるデータについてのパワースペクトルを示す図である。図６のパワースペクトルは、０Ｈｚ～１６００Ｈｚまでの周波数の信号のパワーを示すパワースペクトルである。
図５と図６とのパワースペクトルを見てみると、図６のパワースペクトルでは、図５のパワースペクトルと比べて約８００Ｈｚ以上の周波数の信号が顕著に減衰していることが分かる。図６のパワースペクトルでは、推移データｙについての自己相関行列の固有値のうち、支配的なものだけを用いて求められた修正された自己回帰モデルについてのパワースペクトルである。そのため、推移データｙについて、本質的な成分は、０Ｈｚ～約８００Ｈｚの信号であると判断できる。そこで、軸受１０６に何らかの異常が生じて、推移データｙがその異常の影響を受ける場合、その影響は、推移データｙにおける本質的な成分である０Ｈｚ～約８００Ｈｚの信号（修正された自己回帰モデルにおける本質的な成分）に表れると判断できる。
図７は、外輪に疵のある軸受１０６を用いた場合に、計測データｖの解析信号ｖaから式３を用いて求められた推移データｙについてのパワースペクトルを示す図である。図７（ａ）のパワースペクトルは、０Ｈｚ～４００００Ｈｚまでの周波数の信号のパワーを示すパワースペクトルである。図７（ｂ）のパワースペクトルは、０Ｈｚ～１８０００Ｈｚまでの周波数の信号のパワーを示すパワースペクトルであり、図７（ａ）のパワースペクトルの一部を拡大したパワースペクトルである。図７（ｃ）のパワースペクトルは、０Ｈｚ～１６００Ｈｚまでの周波数の信号のパワーを示すパワースペクトルであり、図７（ａ）のパワースペクトルの一部を拡大したパワースペクトルである。
図８は、外輪に疵のある軸受１０６を用いた場合に求められた推移データｙから、上記手順で求められた係数αで特定される修正された自己回帰モデルにより求められる推移データの予測値からなるデータについてのパワースペクトルを示す図である。図８（ａ）のパワースペクトルは、０Ｈｚ～４００００Ｈｚまでの周波数の信号のパワーを示すパワースペクトルである。図８（ｂ）のパワースペクトルは、０Ｈｚ～１６００Ｈｚまでの周波数の信号のパワーを示すパワースペクトルであり、図８（ａ）のパワースペクトルの一部を拡大したパワースペクトルである。
図８は、外輪に疵のある軸受１０６を用いた場合に求められた推移データｙから、上記手順で求められた係数αで特定される修正された自己回帰モデルにより求められる推移データの予測値からなるデータについてのパワースペクトルを示す図である。図８（ａ）のパワースペクトルは、０Ｈｚ～４００００Ｈｚまでの周波数の信号のパワーを示すパワースペクトルである。図８（ｂ）のパワースペクトルは、０Ｈｚ～１６００Ｈｚまでの周波数の信号のパワーを示すパワースペクトルであり、図８（ａ）のパワースペクトルの一部を拡大したパワースペクトルである。
図７（ｃ）のパワースペクトルを見ると、２６０Ｈｚ付近と、その倍の周波数である５２０Ｈｚ付近と、その３倍の周波数である７８０Ｈｚ付近と、にピークが見られ、軸受１０６の外輪の疵に対応する信号が出ているのが分かる。しかし、図７（ｂ）のパワースペクトルを見ると、推移データｙの本質的な成分と判断される０Ｈｚ～約８００Ｈｚの範囲以上の周波数において、０Ｈｚ～約８００Ｈｚの範囲で見られるピークと同程度の強度のピークが多数存在していることが分かる。また、図７（ａ）のパワースペクトルを見ると、図７（ｂ）に示される周波数（１８０００Ｈｚ）以上の周波数においても、ピークが多数存在していることが分かる。疵の種類によっては１８０００Ｈｚ以上の範囲においてピークが生じる場合もあることから、推移データｙのパワースペクトルに基づいて軸受１０６の異常を診断する場合には、異常が過検知されてしまう可能性がある。
図８（ｂ）のパワースペクトルと図６のパワースペクトルとを比べてみると、図８（ｂ）のパワースペクトルには、軸受１０６の外輪の疵の回転周期に対応する周波数である２６０Ｈｚ付近と、その倍の５２０Ｈｚ付近と、にピークが立っていることが分かる。即ち、図８（ｂ）のパワースペクトルに、軸受１０６の外輪の疵に係る回転周波数に対応するピークが生じていることとなる。
また、図８（ａ）のパワースペクトルと、図７（ａ）及び（ｂ）のパワースペクトルと、を見ると、図８（ａ）のパワースペクトルにおける８００Ｈｚ以上の周波数の成分が、図７（ａ）及び（ｂ）のパワースペクトルに比べて、顕著に減衰していることが分かる。そのため、発明者らは、図８のパワースペクトルに基づいて軸受１０６の異常を診断することで、図７のパワースペクトルに基づいて軸受１０６の異常を診断する場合に比べて、異常を過検知してしまうような事態を防止できるという知見を得た。
図８（ｂ）のパワースペクトルと図６のパワースペクトルとを比べてみると、図８（ｂ）のパワースペクトルには、軸受１０６の外輪の疵の回転周期に対応する周波数である２６０Ｈｚ付近と、その倍の５２０Ｈｚ付近と、にピークが立っていることが分かる。即ち、図８（ｂ）のパワースペクトルに、軸受１０６の外輪の疵に係る回転周波数に対応するピークが生じていることとなる。
また、図８（ａ）のパワースペクトルと、図７（ａ）及び（ｂ）のパワースペクトルと、を見ると、図８（ａ）のパワースペクトルにおける８００Ｈｚ以上の周波数の成分が、図７（ａ）及び（ｂ）のパワースペクトルに比べて、顕著に減衰していることが分かる。そのため、発明者らは、図８のパワースペクトルに基づいて軸受１０６の異常を診断することで、図７のパワースペクトルに基づいて軸受１０６の異常を診断する場合に比べて、異常を過検知してしまうような事態を防止できるという知見を得た。
図９は、外輪に疵のある軸受１０６を用いた場合に求められた推移データｙにローパスフィルタをかけたデータについてのパワースペクトルを示す図である。
図９（ａ）のパワースペクトルは、推移データｙに対して、１５００Ｈｚ以上の周波数をカットするよう設計されたローパスフィルタをかけたデータについてのパワースペクトルである。図９（ｂ）のパワースペクトルは、推移データｙに対して、８００Ｈｚ以上の周波数をカットするよう設計されたローパスフィルタをかけたデータについてのパワースペクトルである。図９（ｃ）のパワースペクトルは、推移データｙに対して、６００Ｈｚ以上の周波数をカットするよう設計されたローパスフィルタをかけたデータについてのパワースペクトルである。
図９（ａ）のパワースペクトルを見ると、２６０Ｈｚ付近と、５２０Ｈｚ付近と、において、他のピークに比べて顕著に高いピークが存在し、軸受１０６の外輪の疵に対応する信号が出ているのが分かる。
図９（ａ）のパワースペクトルは、推移データｙに対して、１５００Ｈｚ以上の周波数をカットするよう設計されたローパスフィルタをかけたデータについてのパワースペクトルである。図９（ｂ）のパワースペクトルは、推移データｙに対して、８００Ｈｚ以上の周波数をカットするよう設計されたローパスフィルタをかけたデータについてのパワースペクトルである。図９（ｃ）のパワースペクトルは、推移データｙに対して、６００Ｈｚ以上の周波数をカットするよう設計されたローパスフィルタをかけたデータについてのパワースペクトルである。
図９（ａ）のパワースペクトルを見ると、２６０Ｈｚ付近と、５２０Ｈｚ付近と、において、他のピークに比べて顕著に高いピークが存在し、軸受１０６の外輪の疵に対応する信号が出ているのが分かる。
図９（ｂ）のパワースペクトルを見ると、２６０Ｈｚ付近において、他のピークに比べて顕著に高いピークが存在することが分かる。しかし、５２０Ｈｚ付近のピークは、０Ｈｚ～２００Ｈｚ付近で見られるピークと同程度のピークとなっており、０Ｈｚ～２００Ｈｚ付近で見られるピークとの違いが明確でないため、５２０Ｈｚ付近のピークが未検知となってしまう可能性がある。
図９（ｃ）のパワースペクトルを見ると、２６０Ｈｚ付近において、他のピークに比べて顕著に高いピークが存在することが分かる。しかし、５２０Ｈｚ付近のピークは、０Ｈｚ～２００Ｈｚ付近で見られるピークよりも低いピークとなっており、０Ｈｚ～２００Ｈｚ付近で見られるピークとよりも小さい重要でないピークとして、５２０Ｈｚ付近のピークが未検知となってしまう可能性がある。
図９（ａ）のパワースペクトルには、軸受１０６の外輪の疵に対応する信号のピークが見られるため、図９（ａ）のパワースペクトルに基づいて、軸受１０６の診断を適切に行うことができる。しかし、図９（ｂ）、（ｃ）のパワースペクトルでは、軸受１０６の外輪の疵に対応する信号のピークのうち５２０Ｈｚ付近のピークが未検知となる可能性があるため、軸受１０６の診断を適切に行うことができない可能性がある。
図９（ｃ）のパワースペクトルを見ると、２６０Ｈｚ付近において、他のピークに比べて顕著に高いピークが存在することが分かる。しかし、５２０Ｈｚ付近のピークは、０Ｈｚ～２００Ｈｚ付近で見られるピークよりも低いピークとなっており、０Ｈｚ～２００Ｈｚ付近で見られるピークとよりも小さい重要でないピークとして、５２０Ｈｚ付近のピークが未検知となってしまう可能性がある。
図９（ａ）のパワースペクトルには、軸受１０６の外輪の疵に対応する信号のピークが見られるため、図９（ａ）のパワースペクトルに基づいて、軸受１０６の診断を適切に行うことができる。しかし、図９（ｂ）、（ｃ）のパワースペクトルでは、軸受１０６の外輪の疵に対応する信号のピークのうち５２０Ｈｚ付近のピークが未検知となる可能性があるため、軸受１０６の診断を適切に行うことができない可能性がある。
このように、ローパスフィルタを用いる場合、カットする周波数によっては、診断に重要な成分の信号が減衰することで、その信号が検知できないようになる可能性があることが確認された。また、ローパスフィルタをどの周波数をカットするように設計するかについては、ユーザの経験による判断によって決定せざるを得ない。そのため、ローパスフィルタを適切な周波数をカットするように設計することは困難である。
対して、本実験で行った方法では、計測データｖから推移データｙを取得し、取得した推移データｙの自己相関行列Ｒの固有値から定量的に選択された支配的な固有値を用いて、推移データｙの修正された自己回帰モデルの係数αを求め、求めた係数αに基づいて診断を行った。このように、本実験で行った方法では、ユーザによる定性的な判断によるパラメータの決定を要しないため、ユーザの判断によって診断の質が低下するような事態を防止できる。また、図８（ｂ）のパワースペクトルを見てみると、自己相関行列Ｒの固有値から定量的に選択された支配的な固有値を用いて決定された係数αで特定される自己回帰モデルにより、診断に重要な信号（２６０Ｈｚ付近と５２０Ｈｚ付近とのピークに対応する信号）が維持されていることが分かる。
対して、本実験で行った方法では、計測データｖから推移データｙを取得し、取得した推移データｙの自己相関行列Ｒの固有値から定量的に選択された支配的な固有値を用いて、推移データｙの修正された自己回帰モデルの係数αを求め、求めた係数αに基づいて診断を行った。このように、本実験で行った方法では、ユーザによる定性的な判断によるパラメータの決定を要しないため、ユーザの判断によって診断の質が低下するような事態を防止できる。また、図８（ｂ）のパワースペクトルを見てみると、自己相関行列Ｒの固有値から定量的に選択された支配的な固有値を用いて決定された係数αで特定される自己回帰モデルにより、診断に重要な信号（２６０Ｈｚ付近と５２０Ｈｚ付近とのピークに対応する信号）が維持されていることが分かる。
このように本発明者らは、推移データｙの自己相関行列Ｒの固有値のうち、支配的な固有値のみを用いて、推移データｙの修正された自己回帰モデルの係数αを求めることで、推移データｙに含まれる成分のうち、物体の異常診断に有用な成分を抽出したデータを求めることができるとの知見を得ることができた。
本実施形態の処理は、以下のような処理である。即ち、本実験で得られた知見を基づいて、物体の振動の計測データの振幅の時間的な推移の様子を示す推移データから求められる自己相関行列Ｒを特異値分解し、得られる固有値のうちの支配的な固有値を用いて、推移データを近似する修正された自己回帰モデルの係数を求める。次に、求めた係数で特定される修正された自己回帰モデルから、式５の右辺の計算を行って推移データの予測値を求め、この推移データの予測値からなるデータについて、フーリエ変換を施すことで、この推移データの予測値からなるデータのパワースペクトルを取得する。そして、取得したパワースペクトルに基づいて、物体の異常を診断する処理である。
本実施形態の処理は、以下のような処理である。即ち、本実験で得られた知見を基づいて、物体の振動の計測データの振幅の時間的な推移の様子を示す推移データから求められる自己相関行列Ｒを特異値分解し、得られる固有値のうちの支配的な固有値を用いて、推移データを近似する修正された自己回帰モデルの係数を求める。次に、求めた係数で特定される修正された自己回帰モデルから、式５の右辺の計算を行って推移データの予測値を求め、この推移データの予測値からなるデータについて、フーリエ変換を施すことで、この推移データの予測値からなるデータのパワースペクトルを取得する。そして、取得したパワースペクトルに基づいて、物体の異常を診断する処理である。
（システム構成）
図１０は、本実施形態の診断システムのシステム構成の一例を示す図である。診断システムは、周期的な運動を行う物体についての異常診断を行うシステムである。診断システムは、情報処理装置１０００、振動計測装置１００１を含む。本実施形態では、診断システムは、鉄道台車に利用されている軸受の異常診断を行う。
情報処理装置１０００は、振動計測装置１００１により計測された信号に基づいて、軸受の異常診断を行うパーソナルコンピュータ（ＰＣ）、サーバ装置、タブレット装置等の情報処理装置である。また、情報処理装置１０００は、電車に組み込まれたコンピュータ等であってもよい。
振動計測装置１００１は、加速度センサ等のセンサを含み、センサを介して物体の振動を検知し、検知した振動に応じた信号を有線又は無線で情報処理装置１０００等の外部の装置に出力する計測装置である。
図１０は、本実施形態の診断システムのシステム構成の一例を示す図である。診断システムは、周期的な運動を行う物体についての異常診断を行うシステムである。診断システムは、情報処理装置１０００、振動計測装置１００１を含む。本実施形態では、診断システムは、鉄道台車に利用されている軸受の異常診断を行う。
情報処理装置１０００は、振動計測装置１００１により計測された信号に基づいて、軸受の異常診断を行うパーソナルコンピュータ（ＰＣ）、サーバ装置、タブレット装置等の情報処理装置である。また、情報処理装置１０００は、電車に組み込まれたコンピュータ等であってもよい。
振動計測装置１００１は、加速度センサ等のセンサを含み、センサを介して物体の振動を検知し、検知した振動に応じた信号を有線又は無線で情報処理装置１０００等の外部の装置に出力する計測装置である。
（情報処理装置の機能構成）
図１０を参照しながら、情報処理装置１０００が有する機能の一例を説明する。
情報処理装置１０００は、取得部１０１０、決定部１０２０、診断部１０３０、出力部１０４０を含む。
≪取得部１０１０≫
取得部１０１０は、周期的な運動を行う物体について、この周期的な運動に係る計測データｖを取得し、取得した計測データｖに基づいて、計測データｖにおける振幅の時間的な推移の様子を示す推移データｙを取得する。
本実施形態では、取得部１０１０は、小歯車軸１００を回転させ、位置１１１に設置された振動計測装置１００１により計測された振動の計測データｖを取得し、取得した計測データｖに基づいて、計測データにおける振幅の時間的な推移の様子を示す推移データｙを取得する。
図１０を参照しながら、情報処理装置１０００が有する機能の一例を説明する。
情報処理装置１０００は、取得部１０１０、決定部１０２０、診断部１０３０、出力部１０４０を含む。
≪取得部１０１０≫
取得部１０１０は、周期的な運動を行う物体について、この周期的な運動に係る計測データｖを取得し、取得した計測データｖに基づいて、計測データｖにおける振幅の時間的な推移の様子を示す推移データｙを取得する。
本実施形態では、取得部１０１０は、小歯車軸１００を回転させ、位置１１１に設置された振動計測装置１００１により計測された振動の計測データｖを取得し、取得した計測データｖに基づいて、計測データにおける振幅の時間的な推移の様子を示す推移データｙを取得する。
≪決定部１０２０≫
決定部１０２０は、取得部１０１０により取得された推移データｙに基づいて、修正された自己回帰モデルにおける係数αを決定する。修正された自己回帰モデルは、推移データｙの実績値ｙk-1～ｙk-mと、この実績値に対する係数α（＝α1～αm）と、を用いて、推移データｙの予測値ｙ＾kを表す式５である。決定部１０２０は、推移データｙから求まる式１１で表される自己相関行列Ｒの固有値から支配的な固有値を決定して、決定した支配的な固有値に基づいて、自己相関行列Ｒから本質的な成分を抽出した式１４に示される行列Ｒ’を生成する。行列Ｒ’は、第１の行列の一例である。決定部１０２０は、修正された自己回帰モデルにおける係数αを、一般的に知られている自己回帰モデルにおける係数の決定に用いる式１０（ユール・ウォーカー方程式）において、式１１で表される自己相関行列Ｒの代わりに、生成した行列Ｒ’（式１４）を用いる式である式１６を用いて決定する。式１０は、時差が０からｍ－１までの推移データｙの自己相関を成分とする自己相関行列Ｒを係数行列とし、時差が１からｍまでの推移データｙの自己相関を成分とする自己相関ベクトルを定数ベクトルとする方程式である。
式１０は、式５で算出される推移データｙの予測値ｙ＾kと、推移データｙの予測値ｙ＾kに対応する時刻ｋにおける推移データｙの実測値ｙkとの二乗誤差を最小化する条件を示す条件式として導出することができる。行列Ｒ’は、自己相関行列Ｒを特異値分解（式１２）することで導出される自己相関行列Ｒの固有値を対角成分とする対角行列Σと、自己相関行列Ｒの固有ベクトルを列成分とする直交行列Ｕと、から導出される。行列Ｒ’は、自己相関行列Ｒの固有値のうち支配的な固有値として決定されたｓ個の固有値を用いて、対角行列Σの部分行列であって、ｓ個の固有値を対角成分とする行列Σsと、直交行列Ｕの部分行列であって、ｓ個の固有値に対応する固有ベクトルを列成分ベクトルとする行列Ｕsと、から導出される行列ＵsΣsＵs Tである。
決定部１０２０は、取得部１０１０により取得された推移データｙに基づいて、修正された自己回帰モデルにおける係数αを決定する。修正された自己回帰モデルは、推移データｙの実績値ｙk-1～ｙk-mと、この実績値に対する係数α（＝α1～αm）と、を用いて、推移データｙの予測値ｙ＾kを表す式５である。決定部１０２０は、推移データｙから求まる式１１で表される自己相関行列Ｒの固有値から支配的な固有値を決定して、決定した支配的な固有値に基づいて、自己相関行列Ｒから本質的な成分を抽出した式１４に示される行列Ｒ’を生成する。行列Ｒ’は、第１の行列の一例である。決定部１０２０は、修正された自己回帰モデルにおける係数αを、一般的に知られている自己回帰モデルにおける係数の決定に用いる式１０（ユール・ウォーカー方程式）において、式１１で表される自己相関行列Ｒの代わりに、生成した行列Ｒ’（式１４）を用いる式である式１６を用いて決定する。式１０は、時差が０からｍ－１までの推移データｙの自己相関を成分とする自己相関行列Ｒを係数行列とし、時差が１からｍまでの推移データｙの自己相関を成分とする自己相関ベクトルを定数ベクトルとする方程式である。
式１０は、式５で算出される推移データｙの予測値ｙ＾kと、推移データｙの予測値ｙ＾kに対応する時刻ｋにおける推移データｙの実測値ｙkとの二乗誤差を最小化する条件を示す条件式として導出することができる。行列Ｒ’は、自己相関行列Ｒを特異値分解（式１２）することで導出される自己相関行列Ｒの固有値を対角成分とする対角行列Σと、自己相関行列Ｒの固有ベクトルを列成分とする直交行列Ｕと、から導出される。行列Ｒ’は、自己相関行列Ｒの固有値のうち支配的な固有値として決定されたｓ個の固有値を用いて、対角行列Σの部分行列であって、ｓ個の固有値を対角成分とする行列Σsと、直交行列Ｕの部分行列であって、ｓ個の固有値に対応する固有ベクトルを列成分ベクトルとする行列Ｕsと、から導出される行列ＵsΣsＵs Tである。
≪診断部１０３０≫
診断部１０３０は、決定部１０２０により決定された係数αに基づいて、物体の異常を診断する。診断部１０３０は、決定部１０２０により決定された係数αで特定される修正された自己回帰モデルから、推移データの予測値を求め、この推移データの予測値からなるデータについて、フーリエ変換を施すことで、この推移データの予測値からなるデータのパワースペクトルを求め、このパワースペクトルに基づいて、物体の異常を診断してもよい。更に、診断部１０３０は、パワースペクトルにおけるピークを示す周波数に基づいて、物体における異常の有無を診断の結果としてもよいし、物体における異常を有する部位を診断の結果としてもよい。
本実施形態では、軸受１０６における疵の有無を診断の結果とする。
≪出力部１０４０≫
出力部１０４０は、診断部１０３０による診断の結果に係る情報を出力する。
診断部１０３０は、決定部１０２０により決定された係数αに基づいて、物体の異常を診断する。診断部１０３０は、決定部１０２０により決定された係数αで特定される修正された自己回帰モデルから、推移データの予測値を求め、この推移データの予測値からなるデータについて、フーリエ変換を施すことで、この推移データの予測値からなるデータのパワースペクトルを求め、このパワースペクトルに基づいて、物体の異常を診断してもよい。更に、診断部１０３０は、パワースペクトルにおけるピークを示す周波数に基づいて、物体における異常の有無を診断の結果としてもよいし、物体における異常を有する部位を診断の結果としてもよい。
本実施形態では、軸受１０６における疵の有無を診断の結果とする。
≪出力部１０４０≫
出力部１０４０は、診断部１０３０による診断の結果に係る情報を出力する。
（情報処理装置のハードウェア構成）
図１１は、情報処理装置１０００のハードウェア構成の一例を示す図である。
情報処理装置１０００は、ＣＰＵ１１００、主記憶装置１１０１、補助記憶装置１１０２、入出力Ｉ／Ｆ１１０３を含む。各構成要素は、システムバス１１０４を介して、相互に通信可能に接続されている。
ＣＰＵ１１００は、情報処理装置１０００を制御する中央演算装置である。主記憶装置１１０１は、ＣＰＵ１１００のワークエリアやデータの一時的な保管場所として機能するＲＡＭ（Ｒａｎｄｏｍ Ａｃｃｅｓｓ Ｍｅｍｏｒｙ）等の記憶装置である。補助記憶装置１１０２は、各種のプログラム、設定データ、振動計測装置１００１から出力される計測データ、診断情報等を記憶するＲＯＭ（Ｒｅａｄ Ｏｎｌｙ Ｍｅｍｏｒｙ）、ＨＤＤ（Ｈａｒｄ Ｄｉｓｋ Ｄｒｉｖｅ）、ＳＳＤ（Ｓｏｌｉｄ Ｓｔａｔｅ Ｄｒｉｖｅ）等の記憶装置である。入出力Ｉ／Ｆ１１０３は、振動計測装置１００１等の外部の装置との間での情報のやり取りに利用されるインターフェースである。
ＣＰＵ１１００が、補助記憶装置１１０２等に記憶されたプログラムに基づき処理を実行することによって、図１０で説明した情報処理装置１０００の機能及び図１２で後述するフローチャートの処理等が実現される。
図１１は、情報処理装置１０００のハードウェア構成の一例を示す図である。
情報処理装置１０００は、ＣＰＵ１１００、主記憶装置１１０１、補助記憶装置１１０２、入出力Ｉ／Ｆ１１０３を含む。各構成要素は、システムバス１１０４を介して、相互に通信可能に接続されている。
ＣＰＵ１１００は、情報処理装置１０００を制御する中央演算装置である。主記憶装置１１０１は、ＣＰＵ１１００のワークエリアやデータの一時的な保管場所として機能するＲＡＭ（Ｒａｎｄｏｍ Ａｃｃｅｓｓ Ｍｅｍｏｒｙ）等の記憶装置である。補助記憶装置１１０２は、各種のプログラム、設定データ、振動計測装置１００１から出力される計測データ、診断情報等を記憶するＲＯＭ（Ｒｅａｄ Ｏｎｌｙ Ｍｅｍｏｒｙ）、ＨＤＤ（Ｈａｒｄ Ｄｉｓｋ Ｄｒｉｖｅ）、ＳＳＤ（Ｓｏｌｉｄ Ｓｔａｔｅ Ｄｒｉｖｅ）等の記憶装置である。入出力Ｉ／Ｆ１１０３は、振動計測装置１００１等の外部の装置との間での情報のやり取りに利用されるインターフェースである。
ＣＰＵ１１００が、補助記憶装置１１０２等に記憶されたプログラムに基づき処理を実行することによって、図１０で説明した情報処理装置１０００の機能及び図１２で後述するフローチャートの処理等が実現される。
（異常診断処理）
本実施形態では、診断システムは、図１と同様の状況で位置１１１に設置された振動計測装置１００１により計測された振動の計測データに基づいて計測データの振幅の時間的な推移の様子を示す推移デ―タを取得し、取得した推移データに基づいて、軸受１０６の異常を診断することとする。本実施形態では、診断システムは、軸受１０６に疵が存在するか否かを特定することで、軸受１０６の異常を診断することとする。本実施形態では、診断システムは、予め振動計測装置１００１により測定された測定データに基づいて、軸受１０６の異常を診断することとする。
図１２は、情報処理装置１０００の処理の一例を示すフローチャートである。
Ｓ１２０１において、取得部１０１０は、振動計測装置１００１により計測された振動の計測データｖとして、時刻１～Ｍまでの計測データであるｖ1～ｖMを取得する。そして、取得部１０１０は、式４を用いて、取得した計測データｖからフーリエ変換を用いる方法で解析信号ｖaを取得する。取得部１０１０は解析信号ｖaに基づいて、式３を用いて、計測データｖのエンベロープを、推移データｙとして取得する。取得部１０１０は、推移データｙとして、時刻１～Ｍまでの推移データであるｙ1～ｙMを取得する。
また、取得部１０１０は、例えば、ｖ（ｔ）における値がマイナスの部分について、－１を乗じて正の値にした関数ｖ’（ｔ）＝｜ｖ（ｔ）｜を求め、求めたｖ’（ｔ）について、予め定められた閾値以上の周波数の信号をカットするようなローパスフィルタをかけることで求まる信号を、推移データｙとして取得することとしてもよい。また、取得部１０１０は、関数ｖ’（ｔ）における山の部分を特定し、特定した部分を補間（例えば、多項式補間）することで求まる関数が示す信号を、推移データｙとして取得することとしてもよい。
本実施形態では、診断システムは、図１と同様の状況で位置１１１に設置された振動計測装置１００１により計測された振動の計測データに基づいて計測データの振幅の時間的な推移の様子を示す推移デ―タを取得し、取得した推移データに基づいて、軸受１０６の異常を診断することとする。本実施形態では、診断システムは、軸受１０６に疵が存在するか否かを特定することで、軸受１０６の異常を診断することとする。本実施形態では、診断システムは、予め振動計測装置１００１により測定された測定データに基づいて、軸受１０６の異常を診断することとする。
図１２は、情報処理装置１０００の処理の一例を示すフローチャートである。
Ｓ１２０１において、取得部１０１０は、振動計測装置１００１により計測された振動の計測データｖとして、時刻１～Ｍまでの計測データであるｖ1～ｖMを取得する。そして、取得部１０１０は、式４を用いて、取得した計測データｖからフーリエ変換を用いる方法で解析信号ｖaを取得する。取得部１０１０は解析信号ｖaに基づいて、式３を用いて、計測データｖのエンベロープを、推移データｙとして取得する。取得部１０１０は、推移データｙとして、時刻１～Ｍまでの推移データであるｙ1～ｙMを取得する。
また、取得部１０１０は、例えば、ｖ（ｔ）における値がマイナスの部分について、－１を乗じて正の値にした関数ｖ’（ｔ）＝｜ｖ（ｔ）｜を求め、求めたｖ’（ｔ）について、予め定められた閾値以上の周波数の信号をカットするようなローパスフィルタをかけることで求まる信号を、推移データｙとして取得することとしてもよい。また、取得部１０１０は、関数ｖ’（ｔ）における山の部分を特定し、特定した部分を補間（例えば、多項式補間）することで求まる関数が示す信号を、推移データｙとして取得することとしてもよい。
Ｓ１２０２において、決定部１０２０は、Ｓ１２０１で取得した推移データｙと、予め設定された定数Ｍと、修正された自己回帰モデルにおいて、ある時刻のデータを過去の幾つのデータで近似するかを示す数ｍと、に基づいて、式９と式１１とを用いて自己相関行列Ｒを生成する。決定部１０２０は、予め記憶されているＭ、ｍの情報を読み出すことで、Ｍとｍとを取得する。本実施形態では、ｍの値は、５００である。また、Ｍは、ｍよりも大きい整数である。
Ｓ１２０３において、決定部１０２０は、Ｓ１２０２で生成した自己相関行列Ｒを特異値分解することで、式１２の直交行列Ｕと対角行列Σを取得し、対角行列Σから自己相関行列Ｒの固有値σ11～σmmを取得する。
Ｓ１２０３において、決定部１０２０は、Ｓ１２０２で生成した自己相関行列Ｒを特異値分解することで、式１２の直交行列Ｕと対角行列Σを取得し、対角行列Σから自己相関行列Ｒの固有値σ11～σmmを取得する。
Ｓ１２０４において、決定部１０２０は、自己相関行列Ｒの固有値のうちの支配的な固有値を決定する。決定部１０２０は、例えば、行列Ｒの全ての固有値の平均値を求め、行列Ｒの固有値のうち、求めた平均値よりも大きい固有値を全て特定し、特定した固有値を支配的な固有値として決定する。そして、決定部１０２０は、使用固有値数ｓを、支配的な固有値として選択した固有値の数に決定する。
また、決定部１０２０は、例えば、選択される固有値の合計のσ11、σ22、・・・、σmmの総和に対する割合が予め定められた閾値（例えば、９５パーセント等）以上となるように、固有値を選択（例えば、値の大きいものから選択）し、選択した固有値を支配的な固有値として決定してもよい。また、決定部１０２０は、選択された固有値の合計のσ11、σ22、・・・、σmmの総和に対する割合が予め定められた閾値（例えば、７５パーセント等）以上となるように、予め定められた閾値（例えば、１個、３個、５個、固有値全体の個数の１パーセントの個数等）以下の個数の固有値を選択し、選択した固有値を支配的な固有値として決定してもよい。
また、決定部１０２０は、例えば、選択される固有値の合計のσ11、σ22、・・・、σmmの総和に対する割合が予め定められた閾値（例えば、９５パーセント等）以上となるように、固有値を選択（例えば、値の大きいものから選択）し、選択した固有値を支配的な固有値として決定してもよい。また、決定部１０２０は、選択された固有値の合計のσ11、σ22、・・・、σmmの総和に対する割合が予め定められた閾値（例えば、７５パーセント等）以上となるように、予め定められた閾値（例えば、１個、３個、５個、固有値全体の個数の１パーセントの個数等）以下の個数の固有値を選択し、選択した固有値を支配的な固有値として決定してもよい。
Ｓ１２０５において、決定部１０２０は、計測データｙと、Ｓ１２０４で決定した支配的な固有値σ11～σssと、自己相関行列Ｒの特異値分解により得られた直交行列Ｕと、に基づいて、式１７を用いて、修正された自己回帰モデルの係数αを決定する。
Ｓ１２０６において、診断部１０３０は、Ｓ１２０５で決定した係数αにより特定される修正された自己回帰モデルから、推移データの予測値を求め、この推移データの予測値からなるデータについて、フーリエ変換を施すことで、この推移データの予測値からなるデータのパワースペクトルを取得する。
Ｓ１２０６において、診断部１０３０は、Ｓ１２０５で決定した係数αにより特定される修正された自己回帰モデルから、推移データの予測値を求め、この推移データの予測値からなるデータについて、フーリエ変換を施すことで、この推移データの予測値からなるデータのパワースペクトルを取得する。
Ｓ１２０７において、診断部１０３０は、Ｓ１２０６で取得したパワースペクトルに基づいて、軸受１０６の異常を診断する。
本実施形態では、診断部１０３０は、例えば、Ｓ１２０６で取得したパワースペクトルにおいて、予め定められた閾値以上の強度を示す信号に対応する周波数を特定する。そして、診断部１０３０は、特定した周波数の中に、他の周波数の整数倍となっている周波数が存在するか否かを判定する。診断部１０３０は、例えば、特定した周波数の中から順に１つずつ選択して、選択した周波数ごとに、選択した周波数を整数倍（例えば、２倍）した値と、特定した周波数に含まれる他の周波数それぞれの値と、の差を取得する。そして、診断部１０３０は、取得した差の中に予め定められた閾値（例えば、１、５等）以下のものがあったら、特定した周波数の中に、他の周波数の整数倍となっている周波数が存在すると判定し、軸受１０６に疵が生じていると診断する。
また、診断部１０３０は、取得した差の中に予め定められた閾値以下のものがなかったら、特定した周波数の中に、他の周波数の整数倍となっている周波数が存在しないと判定し、軸受１０６が未知の状態にあると診断する。
また、診断部１０３０は、閾値以上の強度を示す信号に対応する周波数が存在しない場合、軸受１０６に異常が生じておらず、正常であると診断する。
Ｓ１２０８において、出力部１０４０は、軸受１０６の異常の診断の結果を示す情報を出力する。出力部１０４０は、例えば、表示装置に、軸受１０６の異常の診断の結果を示す情報を表示することで出力する。また、出力部１０４０は、例えば、軸受１０６の異常の診断の結果を示す情報を、補助記憶装置１１０２に記憶することで出力することとしてもよい。また、出力部１０４０は、例えば、外部のサーバ装置等の設定された送信先に、軸受１０６の異常の診断の結果を示す情報を送信することで出力することとしてもよい。
本実施形態では、診断部１０３０は、例えば、Ｓ１２０６で取得したパワースペクトルにおいて、予め定められた閾値以上の強度を示す信号に対応する周波数を特定する。そして、診断部１０３０は、特定した周波数の中に、他の周波数の整数倍となっている周波数が存在するか否かを判定する。診断部１０３０は、例えば、特定した周波数の中から順に１つずつ選択して、選択した周波数ごとに、選択した周波数を整数倍（例えば、２倍）した値と、特定した周波数に含まれる他の周波数それぞれの値と、の差を取得する。そして、診断部１０３０は、取得した差の中に予め定められた閾値（例えば、１、５等）以下のものがあったら、特定した周波数の中に、他の周波数の整数倍となっている周波数が存在すると判定し、軸受１０６に疵が生じていると診断する。
また、診断部１０３０は、取得した差の中に予め定められた閾値以下のものがなかったら、特定した周波数の中に、他の周波数の整数倍となっている周波数が存在しないと判定し、軸受１０６が未知の状態にあると診断する。
また、診断部１０３０は、閾値以上の強度を示す信号に対応する周波数が存在しない場合、軸受１０６に異常が生じておらず、正常であると診断する。
Ｓ１２０８において、出力部１０４０は、軸受１０６の異常の診断の結果を示す情報を出力する。出力部１０４０は、例えば、表示装置に、軸受１０６の異常の診断の結果を示す情報を表示することで出力する。また、出力部１０４０は、例えば、軸受１０６の異常の診断の結果を示す情報を、補助記憶装置１１０２に記憶することで出力することとしてもよい。また、出力部１０４０は、例えば、外部のサーバ装置等の設定された送信先に、軸受１０６の異常の診断の結果を示す情報を送信することで出力することとしてもよい。
（まとめ）
以上、本実施形態では、診断システムは、物体の周期的な運動から計測された計測データに基づいて、計測データの振幅の時間的な推移の様子を示す推移データを取得し、取得した推移データから、自己相関行列Ｒを生成し、自己相関行列Ｒを特異値分解し、得られた固有値から定量的に支配的な固有値を決定した。そして、診断システムは、決定した支配的な固有値を用いて、推移データを近似する修正された自己回帰モデルの係数αを決定した。そして、診断システムは、決定した係数αで特定される修正された自己回帰モデルにより推移データの予測値を求め、この推移データの予測値からなるデータのパワースペクトルを求め、求めたパワースペクトルに基づいて、軸受１０６の異常を診断した。診断システムは、自己相関行列Ｒの固有値のうち一部の固有値を用いることで、軸受１０６の異常診断に有用な成分が残り、有用でない成分が残らないように、推移データを近似する修正された自己回帰モデルの係数αを決定できる。これにより、診断システムは、推移データの軸受１０６の異常診断に有用な成分に基づいて、より精度よく異常診断を行うことができる。
また、本実施形態の処理により、診断システムは、ローパスフィルタを用いる場合のようにカットする周波数を定性的に決定する必要もないため、ローパスフィルタを用いた場合のように診断に重要な信号が検知できない程減衰する事態を防止できる。
また、本実施形態の処理により、診断システムは、予め異常診断に用いる信号がどのような周波数の信号であるか等の想定を要することなく、推移データから異常診断に有用な成分を、抽出できる。
以上、本実施形態では、診断システムは、物体の周期的な運動から計測された計測データに基づいて、計測データの振幅の時間的な推移の様子を示す推移データを取得し、取得した推移データから、自己相関行列Ｒを生成し、自己相関行列Ｒを特異値分解し、得られた固有値から定量的に支配的な固有値を決定した。そして、診断システムは、決定した支配的な固有値を用いて、推移データを近似する修正された自己回帰モデルの係数αを決定した。そして、診断システムは、決定した係数αで特定される修正された自己回帰モデルにより推移データの予測値を求め、この推移データの予測値からなるデータのパワースペクトルを求め、求めたパワースペクトルに基づいて、軸受１０６の異常を診断した。診断システムは、自己相関行列Ｒの固有値のうち一部の固有値を用いることで、軸受１０６の異常診断に有用な成分が残り、有用でない成分が残らないように、推移データを近似する修正された自己回帰モデルの係数αを決定できる。これにより、診断システムは、推移データの軸受１０６の異常診断に有用な成分に基づいて、より精度よく異常診断を行うことができる。
また、本実施形態の処理により、診断システムは、ローパスフィルタを用いる場合のようにカットする周波数を定性的に決定する必要もないため、ローパスフィルタを用いた場合のように診断に重要な信号が検知できない程減衰する事態を防止できる。
また、本実施形態の処理により、診断システムは、予め異常診断に用いる信号がどのような周波数の信号であるか等の想定を要することなく、推移データから異常診断に有用な成分を、抽出できる。
本実施形態では、診断システムは、鉄道台車の軸受の異常診断を行うこととした。しかし、診断システムは、歯車等の回転体、振動子等の振動体、バネ等の伸縮体、生物の心臓等の周期的な運動を行う他の物体についての異常診断を行うこととしてもよい。
本実施形態では、ｍは、予め設定された５００という数であるとした。しかし、診断対象に応じて実験を行う等して適切に計測データを近似できる過去のデータの数を特定し、特定した数をｍの値としてもよい。例えば、情報処理装置１０００は、診断対象から計測された計測データから求められた推移データについて、ある時刻の推移データを、その時刻よりも過去の推移データを複数用いて自己回帰モデル等を利用して近似して、近似した値とその時刻の推移データとの差分が設定された閾値未満となる場合の近似に用いられた過去の推移データの数を、ｍとして決定してもよい。
本実施形態では、診断システムは、予め振動計測装置１００１により測定された測定データに基づいて、軸受１０６の異常を診断することとした。しかし、診断システムは、リアルタイムで振動計測装置１００１により測定されて出力され続けている測定データを用いて、稼働中の軸受１０６の異常を診断することとしてもよい。即ち、情報処理装置１０００は振動計測装置１００１から新たな計測データｖが取得される度に、取得された計測データｖから推移データｙ1～ｙMを更新して、図１２の処理を行うようにしてもよい。
本実施形態では、情報処理装置１０００は、補助記憶装置１１０２に記憶されたプログラムを実行することで、図１２の処理を実現することとした。しかし、情報処理装置１０００は、外付けの記憶媒体や外部の記憶サーバ等に記憶されたプログラムを実行することで、図１２の処理を実現することとしてもよい。
また、例えば、上述した情報処理装置１０００の機能の一部又は全てをハードウェアとして情報処理装置１０００に実装してもよい。
また、以上説明した本発明の実施形態は、何れも本発明を実施するにあたっての具体化の例を示したものに過ぎず、これらによって本発明の技術的範囲が限定的に解釈されてはならないものである。即ち、本発明はその技術思想、またはその主要な特徴から逸脱することなく、様々な形で実施することができる。
本実施形態では、ｍは、予め設定された５００という数であるとした。しかし、診断対象に応じて実験を行う等して適切に計測データを近似できる過去のデータの数を特定し、特定した数をｍの値としてもよい。例えば、情報処理装置１０００は、診断対象から計測された計測データから求められた推移データについて、ある時刻の推移データを、その時刻よりも過去の推移データを複数用いて自己回帰モデル等を利用して近似して、近似した値とその時刻の推移データとの差分が設定された閾値未満となる場合の近似に用いられた過去の推移データの数を、ｍとして決定してもよい。
本実施形態では、診断システムは、予め振動計測装置１００１により測定された測定データに基づいて、軸受１０６の異常を診断することとした。しかし、診断システムは、リアルタイムで振動計測装置１００１により測定されて出力され続けている測定データを用いて、稼働中の軸受１０６の異常を診断することとしてもよい。即ち、情報処理装置１０００は振動計測装置１００１から新たな計測データｖが取得される度に、取得された計測データｖから推移データｙ1～ｙMを更新して、図１２の処理を行うようにしてもよい。
本実施形態では、情報処理装置１０００は、補助記憶装置１１０２に記憶されたプログラムを実行することで、図１２の処理を実現することとした。しかし、情報処理装置１０００は、外付けの記憶媒体や外部の記憶サーバ等に記憶されたプログラムを実行することで、図１２の処理を実現することとしてもよい。
また、例えば、上述した情報処理装置１０００の機能の一部又は全てをハードウェアとして情報処理装置１０００に実装してもよい。
また、以上説明した本発明の実施形態は、何れも本発明を実施するにあたっての具体化の例を示したものに過ぎず、これらによって本発明の技術的範囲が限定的に解釈されてはならないものである。即ち、本発明はその技術思想、またはその主要な特徴から逸脱することなく、様々な形で実施することができる。
１０００ 情報処理装置
１０１０ 取得部
１０２０ 決定部
１０３０ 診断部
１０４０ 出力部
１００１ 振動計測装置
１１００ ＣＰＵ
１０１０ 取得部
１０２０ 決定部
１０３０ 診断部
１０４０ 出力部
１００１ 振動計測装置
１１００ ＣＰＵ
Claims (11)
- 周期的な運動を行う物体の前記周期的な運動に係る計測データに基づいて、前記計測データの振幅の時間的な推移の様子を示す推移データを取得する取得手段と、
前記取得手段により取得された前記推移データに基づいて、修正された自己回帰モデルにおける係数を決定する決定手段と、
前記決定手段により決定された前記係数に基づいて、前記物体の異常を診断する診断手段と、
を有し、
前記修正された自己回帰モデルは、前記推移データの実績値と、前記実績値に対する前記係数と、を用いて、前記推移データの予測値を表す式であり、
前記決定手段は、前記推移データから導出される自己相関行列を特異値分解することで導出される前記自己相関行列の固有値を対角成分とする対角行列と、前記自己相関行列の固有ベクトルを列成分とする直交行列と、から導出される第１の行列を係数行列とし、前記推移データから導出される自己相関ベクトルを定数ベクトルとする方程式を用いて、前記係数を決定し、
前記自己相関ベクトルは、時差が１から前記修正された自己回帰モデルで用いられる前記実績値の数であるｍまでの前記推移データの自己相関を成分とするベクトルであり、
前記自己相関行列は、時差が０からｍ－１までの前記推移データの自己相関を成分とする行列であり、
前記第１の行列は、前記自己相関行列の固有値のうち支配的な固有値として決定されたｓ個の固有値と前記対角行列とから導出される第２の行列Σsと、前記ｓ個の固有値と前記直交行列とから導出される第３の行列Ｕsと、から導出される行列ＵsΣsＵs Tであり、
前記第２の行列は、前記対角行列の部分行列であって、前記ｓ個の固有値を対角成分とする行列であり、
前記第３の行列は、前記直交行列の部分行列であって、前記ｓ個の固有値に対応する固有ベクトルを列成分ベクトルとする行列である情報処理装置。 - 前記取得手段は、前記計測データに基づいて、前記計測データのエンベロープのデータを、前記推移データとして取得する請求項１記載の情報処理装置。
- 前記ｓ個の固有値は、前記自己相関行列のうち、前記自己相関行列の全ての固有値の平均値よりも大きい固有値である請求項１又は２記載の情報処理装置。
- 前記ｓ個の固有値は、予め定められた閾値以下の個数の固有値であり、
前記ｓ個の固有値の合計の前記自己相関行列の固有値の全体の合計に対する割合は、予め定められた閾値以下となる請求項１又は２記載の情報処理装置。 - 前記自己相関行列は、前記推移データの予測値と前記推移データの予測値に対応する時刻における前記推移データの実測値との二乗誤差を最小化する条件を示す条件式に含まれる行列である請求項１乃至４何れか１項記載の情報処理装置。
- 前記診断手段は、前記修正された自己回帰モデルにより求められる前記推移データの予測値からなるデータのパワースペクトルを導出し、前記パワースペクトルに基づいて、前記物体の異常を診断する請求項１乃至５何れか１項記載の情報処理装置。
- 前記計測データは、前記物体の振動に係る計測データである請求項１乃至６何れか１項記載の情報処理装置。
- 前記診断手段による診断の結果に係る情報を出力する出力手段を更に有する請求項１乃至７何れか１項記載の情報処理装置。
- 前記物体は、鉄道台車に利用される軸受であり、
前記診断手段は、前記決定手段により決定された前記係数に基づいて、前記物体における疵の存在を含む異常を診断する請求項１乃至８何れか１項記載の情報処理装置。 - 情報処理装置が実行する情報処理方法であって、
周期的な運動を行う物体の前記周期的な運動に係る計測データに基づいて、前記計測データの振幅の時間的な推移の様子を示す推移データを取得する取得ステップと、
前記取得ステップで取得された前記推移データに基づいて、修正された自己回帰モデルにおける係数を決定する決定ステップと、
前記決定ステップで決定された前記係数に基づいて、前記物体の異常を診断する診断ステップと、
を含み、
前記修正された自己回帰モデルは、前記推移データの実績値と、前記実績値に対する前記係数と、を用いて、前記推移データの予測値を表す式であり、
前記決定ステップでは、前記推移データから導出される自己相関行列を特異値分解することで導出される前記自己相関行列の固有値を対角成分とする対角行列と、前記自己相関行列の固有ベクトルを列成分とする直交行列と、から導出される第１の行列を係数行列とし、前記推移データから導出される自己相関ベクトルを定数ベクトルとする方程式を用いて、前記係数を決定し、
前記自己相関ベクトルは、時差が１から前記修正された自己回帰モデルで用いられる前記実績値の数であるｍまでの前記推移データの自己相関を成分とするベクトルであり、
前記自己相関行列は、時差が０からｍ－１までの前記推移データの自己相関を成分とする行列であり、
前記第１の行列は、前記自己相関行列の固有値のうち支配的な固有値として決定されたｓ個の固有値と前記対角行列とから導出される第２の行列Σsと、前記ｓ個の固有値と前記直交行列とから導出される第３の行列Ｕsと、から導出される行列ＵsΣsＵs Tであり、
前記第２の行列は、前記対角行列の部分行列であって、前記ｓ個の固有値を対角成分とする行列であり、
前記第３の行列は、前記直交行列の部分行列であって、前記ｓ個の固有値に対応する固有ベクトルを列成分ベクトルとする行列である情報処理方法。 - コンピュータを、請求項１乃至９何れか１項記載の情報処理装置の各手段として機能させるためのプログラム。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2018141513A JP7014080B2 (ja) | 2018-07-27 | 2018-07-27 | 情報処理装置、情報処理方法及びプログラム |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2018141513A JP7014080B2 (ja) | 2018-07-27 | 2018-07-27 | 情報処理装置、情報処理方法及びプログラム |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JP2020016611A JP2020016611A (ja) | 2020-01-30 |
| JP7014080B2 true JP7014080B2 (ja) | 2022-02-01 |
Family
ID=69580336
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2018141513A Active JP7014080B2 (ja) | 2018-07-27 | 2018-07-27 | 情報処理装置、情報処理方法及びプログラム |
Country Status (1)
| Country | Link |
|---|---|
| JP (1) | JP7014080B2 (ja) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP7484476B2 (ja) | 2020-06-18 | 2024-05-16 | 株式会社ジェイテクト | 測定装置、転がり軸受の荷重測定方法、及び歯車機構 |
| JP7383367B1 (ja) | 2023-03-07 | 2023-11-20 | 新川センサテクノロジ株式会社 | 回転機器の振動データ分析方法及び分析システム |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2002090266A (ja) | 2000-09-20 | 2002-03-27 | Mitsui Eng & Shipbuild Co Ltd | 余寿命予測装置 |
| JP2009257862A (ja) | 2008-04-15 | 2009-11-05 | Original Engineering Consultants Co Ltd | 回転機械等の設備の音信号による健全性診断方法 |
| WO2016117358A1 (ja) | 2015-01-21 | 2016-07-28 | 三菱電機株式会社 | 検査データ処理装置および検査データ処理方法 |
| JP2017181203A (ja) | 2016-03-29 | 2017-10-05 | 巴バルブ株式会社 | バルブ診断方法及びバルブ診断装置 |
-
2018
- 2018-07-27 JP JP2018141513A patent/JP7014080B2/ja active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2002090266A (ja) | 2000-09-20 | 2002-03-27 | Mitsui Eng & Shipbuild Co Ltd | 余寿命予測装置 |
| JP2009257862A (ja) | 2008-04-15 | 2009-11-05 | Original Engineering Consultants Co Ltd | 回転機械等の設備の音信号による健全性診断方法 |
| WO2016117358A1 (ja) | 2015-01-21 | 2016-07-28 | 三菱電機株式会社 | 検査データ処理装置および検査データ処理方法 |
| JP2017181203A (ja) | 2016-03-29 | 2017-10-05 | 巴バルブ株式会社 | バルブ診断方法及びバルブ診断装置 |
Also Published As
| Publication number | Publication date |
|---|---|
| JP2020016611A (ja) | 2020-01-30 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| JP6919397B2 (ja) | 情報処理装置、情報処理方法及びプログラム | |
| JP6848813B2 (ja) | 情報処理装置、情報処理方法及びプログラム | |
| US10520397B2 (en) | Methods and apparatuses for defect diagnosis in a mechanical system | |
| EP2904368B1 (en) | Turbine blade fatigue life analysis using non-contact measurement and dynamical response reconstruction techniques | |
| JP6196093B2 (ja) | 軸受装置の振動解析方法、軸受装置の振動解析装置、および転がり軸受の状態監視装置 | |
| US7430483B2 (en) | Method and apparatus for diagnosing a mechanism | |
| KR20100094452A (ko) | 롤링 베어링 대미지 검측 및 자동 식별 방법 | |
| JPH10258974A (ja) | 非定常信号解析装置及び非定常信号解析プログラムを記録した媒体 | |
| JP7470849B2 (ja) | 異常検出装置、異常検出方法、及びプログラム | |
| JP7014080B2 (ja) | 情報処理装置、情報処理方法及びプログラム | |
| US11209487B2 (en) | Rotor diagnostic apparatus, rotor diagnostic method, and rotor diagnostic program | |
| CN112955839A (zh) | 异常检测装置、异常检测方法和程序 | |
| EP3480455B1 (en) | Wind turbine monitoring device, wind turbine monitoring method, wind turbine monitoring program, and storage medium | |
| JP6511573B1 (ja) | 転がり軸受の異常診断方法及び異常診断装置、異常診断プログラム | |
| KR101543146B1 (ko) | 진동 장치의 상태 판단 방법 | |
| JP7188143B2 (ja) | 異常予兆検出システム、異常予兆検出方法 | |
| Pawlik | The use of the acoustic signal to diagnose machines operated under variable load | |
| JP2003344202A (ja) | データ処理装置及びデータ処理方法 | |
| WO2021215169A1 (ja) | データ処理装置、データ処理方法及びプログラム | |
| KR19990067540A (ko) | 비정상신호 해석장치 및 비정상신호 해석프로그램을 기록한매체 | |
| JP7350667B2 (ja) | 異常検出装置、回転機械、異常検出方法、及びプログラム | |
| CN114061947A (zh) | 基于稀疏时频分析的齿轮箱变转速故障诊断方法及*** | |
| JP6192413B2 (ja) | 軸受装置の振動解析方法、軸受装置の振動解析装置、および転がり軸受の状態監視装置 | |
| JP5476413B2 (ja) | 回転機械の健全性診断方法 | |
| Nacib et al. | A comparative study of various methods of gear faults diagnosis |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20210303 |
|
| TRDD | Decision of grant or rejection written | ||
| A977 | Report on retrieval |
Free format text: JAPANESE INTERMEDIATE CODE: A971007 Effective date: 20211220 |
|
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20211221 |
|
| A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20220103 |