JP2019079352A - 情報処理装置、情報処理方法及びプログラム - Google Patents
情報処理装置、情報処理方法及びプログラム Download PDFInfo
- Publication number
- JP2019079352A JP2019079352A JP2017206597A JP2017206597A JP2019079352A JP 2019079352 A JP2019079352 A JP 2019079352A JP 2017206597 A JP2017206597 A JP 2017206597A JP 2017206597 A JP2017206597 A JP 2017206597A JP 2019079352 A JP2019079352 A JP 2019079352A
- Authority
- JP
- Japan
- Prior art keywords
- matrix
- measurement data
- bearing
- coefficient
- autocorrelation
- 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.)
- Granted
Links
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01H—MEASUREMENT OF MECHANICAL VIBRATIONS OR ULTRASONIC, SONIC OR INFRASONIC WAVES
- G01H13/00—Measuring resonant frequency
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01H—MEASUREMENT OF MECHANICAL VIBRATIONS OR ULTRASONIC, SONIC OR INFRASONIC WAVES
- G01H17/00—Measuring mechanical vibrations or ultrasonic, sonic or infrasonic waves, not provided for in the preceding groups
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01M—TESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
- G01M13/00—Testing of machine parts
- G01M13/04—Bearings
- G01M13/045—Acoustic or vibration analysis
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Acoustics & Sound (AREA)
- Testing Of Devices, Machine Parts, Or Other Structures Thereof (AREA)
- Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
- Complex Calculations (AREA)
Abstract
【課題】周期的な運動を行う物体の異常診断をより精度よく行うことを目的とする。【解決手段】周期的な運動を行う物体の前記周期的な運動に係る計測データを取得する取得手段と、前記取得手段により取得された前記計測データに基づいて、修正された自己回帰モデルにおける係数である修正係数を決定する決定手段と、前記決定手段により決定された前記修正係数のパターンである計測パターンに基づいて、前記物体の異常を診断する診断手段と、を有し、前記修正された自己回帰モデルは、前記物体の前記周期的な運動に係る計測データと、前記計測データに対する前記修正係数と、に基づいて、前記計測データの予測値を導出する式であり、前記決定手段は、前記計測データから導出される自己相関行列を特異値分解することで導出される自己相関行列の固有値に基づいて、前記修正係数を決定する。【選択図】図１５
Description
本発明は、情報処理装置、情報処理方法及びプログラムに関する。
周期的な運動を行う物体の異常診断を行う技術として、運動を行っている物体から得られる信号を測定し、その信号の特徴量を算出し、その特徴量から正常であるか異常であるかの判定を行う方法がある。信号の周波数特性を自己回帰モデルで同定することで、より高精度なスペクトル特性を取得する技術がある（特許文献１）。
特許文献１に開示された方法では、様々な周波数域の信号成分を含んだ信号に対して、周波数特性を同定するため、物体の異常に起因する周波数等の特定したい周波数を正確に特定することが困難であり、物体の異常を検知出来ない可能性があるという問題がある。
本発明は、以上のような問題点に鑑みてなされたものであり、周期的な運動を行う物体の異常診断をより精度よく行うことができるようにすることを目的とする。
本発明は、以上のような問題点に鑑みてなされたものであり、周期的な運動を行う物体の異常診断をより精度よく行うことができるようにすることを目的とする。
本発明の情報処理装置は、周期的な運動を行う物体の前記周期的な運動に係る計測データを取得する取得手段と、前記取得手段により取得された前記計測データに基づいて、修正された自己回帰モデルにおける係数である修正係数を決定する決定手段と、前記決定手段により決定された前記修正係数のパターンである計測パターンに基づいて、前記物体の異常を診断する診断手段と、を有し、前記修正された自己回帰モデルは、前記計測データの実績値と、前記実績値に対する前記修正係数と、を用いて、前記計測データの予測値を表す式であり、前記決定手段は、前記計測データから導出される自己相関行列を特異値分解することで導出される前記自己相関行列の固有値を対角成分とする対角行列と、前記自己相関行列の固有ベクトルを列成分とする直交行列と、から導出される第１の行列を係数行列とし、前記計測データから導出される自己相関ベクトルを定数ベクトルとする方程式を用いて、前記修正係数を決定し、前記自己相関ベクトルは、時差が１から前記修正された自己回帰モデルで用いられる前記実績値の数であるｍまでの前記計測データの自己相関を成分とするベクトルであり、前記自己相関行列は、時差が０からｍ−１までの前記計測データの自己相関を成分とする行列であり、前記第１の行列は、１以上且つｍ未満の設定された数であるｓに対して、前記自己相関行列のｓ個の固有値と前記対角行列とから導出される第２の行列Σsと、前記ｓ個の固有値と前記直交行列とから導出される第３の行列Ｕsと、から導出される行列ＵsΣsＵs Tであり、前記第２の行列は、前記対角行列の部分行列であって、前記ｓ個の固有値を対角成分とする行列であり、前記第３の行列は、前記直交行列の部分行列であって、前記ｓ個の固有値に対応する固有ベクトルを列成分ベクトルとする行列である。
本発明によれば、周期的な運動を行う物体の異常診断をより精度よく行うことができる技術を提供することができる。
＜実施形態１＞
以下、本発明の一実施形態について図面に基づいて説明する。
以下、本発明の一実施形態について図面に基づいて説明する。
（本実施形態の処理の概要）
本実施形態では、診断システムが鉄道台車に利用される軸受から測定された振動の信号に基づいて、その軸受の異常を診断する。
診断システムは、軸受の振動に応じた信号を測定し、測定された信号を自己回帰モデルで表現し、自己回帰モデルの係数に関する条件式における行列を特異値分解し、その行列の固有値を求め、求めた固有値のうち最大のものから設定された数のみを用いて、自己回帰モデルの係数を決定する。そして、診断システムは、決定した係数のパターンに基づいて、軸受の異常を診断する。
本実施形態では、診断システムが鉄道台車に利用される軸受から測定された振動の信号に基づいて、その軸受の異常を診断する。
診断システムは、軸受の振動に応じた信号を測定し、測定された信号を自己回帰モデルで表現し、自己回帰モデルの係数に関する条件式における行列を特異値分解し、その行列の固有値を求め、求めた固有値のうち最大のものから設定された数のみを用いて、自己回帰モデルの係数を決定する。そして、診断システムは、決定した係数のパターンに基づいて、軸受の異常を診断する。
（軸受異常の診断実験）
鉄道台車に利用されている軸受の異常診断精度の評価のために行った実験について説明する。
図１（ａ）は、本実験の状況を説明する図である。図１（ａ）の状況は、実験室に鉄道台車における車輪の駆動機構が設置されている状況である。駆動モータ１０１は、駆動モータ軸１０２、歯車継手１０３、小歯車軸１００を介して、小歯車１０５を回転させる。小歯車１０５は、回転することで、大歯車１０７を回転させる。それに合わせて、車軸１０９が回転し、車軸１０９に接続される車輪が回転することになる。本実験では、車軸１０９には、車輪の代わりにダイナモ１１０が接続されている。ダイナモ１１０は、疑似的な走行負荷を与えるために接続されている発電機である。本実験では、ダイナモ１１０を回転させるための負荷を、走行の際の負荷として想定する。
小歯車１０５が嵌められた小歯車軸１００と大歯車１０７が嵌められた車軸１０９とは、歯車箱１０４に固定されている。そして、歯車箱１０４には、小歯車軸との間に軸受１０６が装着されており、車軸１０９との間に軸受１０８が装着されている。図１（ｂ）は、歯車箱１０４に取り付けられている軸受を側面から見た様子を示す図である。
鉄道台車に利用されている軸受の異常診断精度の評価のために行った実験について説明する。
図１（ａ）は、本実験の状況を説明する図である。図１（ａ）の状況は、実験室に鉄道台車における車輪の駆動機構が設置されている状況である。駆動モータ１０１は、駆動モータ軸１０２、歯車継手１０３、小歯車軸１００を介して、小歯車１０５を回転させる。小歯車１０５は、回転することで、大歯車１０７を回転させる。それに合わせて、車軸１０９が回転し、車軸１０９に接続される車輪が回転することになる。本実験では、車軸１０９には、車輪の代わりにダイナモ１１０が接続されている。ダイナモ１１０は、疑似的な走行負荷を与えるために接続されている発電機である。本実験では、ダイナモ１１０を回転させるための負荷を、走行の際の負荷として想定する。
小歯車１０５が嵌められた小歯車軸１００と大歯車１０７が嵌められた車軸１０９とは、歯車箱１０４に固定されている。そして、歯車箱１０４には、小歯車軸との間に軸受１０６が装着されており、車軸１０９との間に軸受１０８が装着されている。図１（ｂ）は、歯車箱１０４に取り付けられている軸受を側面から見た様子を示す図である。
歯車箱１０４上の位置１１１は、振動計測装置が設置される位置である。振動計測装置は、加速度センサ、レーザ変位計等のセンサを含み、センサを介して物体の振動を検知し、検知した振動に応じた信号を出力する。位置１１１に設置された振動計測装置が位置１１１に伝達される振動を計測し、計測した振動を示す信号を外部の情報処理装置等に有線又は無線による通信を介して出力する。本実験では、振動計測装置は、位置１１１に設置されるとするが、軸受１０６に起因する振動が伝達される位置であるならば、任意の位置に設置することとしてもよい。
ここで、軸受１０６について、説明する。図２は、軸受１０６の一例の構成を説明する図である。軸受１０６は、外輪、内輪、転動体、保持器の４つの部分から構成される。図２には、軸受１０６の外輪、内輪、転動体、保持器の概要が示されている。軸受１０６は、外輪と内輪との間に保持器により保持された転動体が挟まれる構成となっている。
本実験では、小歯車１０５に装着されている軸受１０６として、疵のない正常な軸受、内輪に疵のある軸受、保持器に疵のある軸受のそれぞれを利用した。それぞれの軸受１０６毎に、駆動モータ１０１を駆動させ、位置１１１に設置された振動計測装置が振動を計測することで、軸受１０６の異常診断を行うための計測データを取得した。本実験では、駆動モータ１０１は、軸受１０６の内輪の回転数が２９５４．７ｒｐｍ（ｒｏｕｎｄ ｐｅｒ ｍｉｎｕｔｅ）となるように駆動した。
以上のように、軸受１０６が正常な場合、軸受１０６の内輪に疵がある場合、軸受１０６の保持器に疵がある場合、軸受１０６の転動体に疵がある場合、軸受１０６の外輪に疵がある場合について、位置１１１における振動の計測データを取得した。
ここで、軸受１０６について、説明する。図２は、軸受１０６の一例の構成を説明する図である。軸受１０６は、外輪、内輪、転動体、保持器の４つの部分から構成される。図２には、軸受１０６の外輪、内輪、転動体、保持器の概要が示されている。軸受１０６は、外輪と内輪との間に保持器により保持された転動体が挟まれる構成となっている。
本実験では、小歯車１０５に装着されている軸受１０６として、疵のない正常な軸受、内輪に疵のある軸受、保持器に疵のある軸受のそれぞれを利用した。それぞれの軸受１０６毎に、駆動モータ１０１を駆動させ、位置１１１に設置された振動計測装置が振動を計測することで、軸受１０６の異常診断を行うための計測データを取得した。本実験では、駆動モータ１０１は、軸受１０６の内輪の回転数が２９５４．７ｒｐｍ（ｒｏｕｎｄ ｐｅｒ ｍｉｎｕｔｅ）となるように駆動した。
以上のように、軸受１０６が正常な場合、軸受１０６の内輪に疵がある場合、軸受１０６の保持器に疵がある場合、軸受１０６の転動体に疵がある場合、軸受１０６の外輪に疵がある場合について、位置１１１における振動の計測データを取得した。
次に、振動計測装置が計測した計測データを用いた軸受１０６の異常診断の方法を説明する。以下では、振動計測装置が計測した計測データを、計測データｙとする。
本実験では、計測データｙを、修正された自己回帰モデルで近似し、この修正された自己回帰モデルの係数から特徴量を算出することとした。
ある時刻ｋ（１≦ｋ≦Ｍ）における計測データｙの値をｙ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、σ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に置き換えてもよい。
本実験では、軸受１０６が正常な場合、軸受１０６の内輪に疵がある場合、軸受１０６の保持器に疵がある場合、軸受１０６の転動体に疵がある場合、軸受１０６の外輪に疵がある場合それぞれについて、得られた計測データ毎に、式５と式７を用いて自己相関行列Ｒを求め、式８で表される特異値分解を行うことによって自己相関行列Ｒの固有値を求め、自己相関行列Ｒの固有値の分布を求めた。
次に、求めた係数αで特定される修正された自己回帰モデルの周波数特性を求める方法を説明する。
式１で与えられる計測データの予測値ｙ＾kの予測誤差ｘkが白色雑音になるという性質を使うと、修正された自己回帰モデルは白色雑音である予測誤差ｘkを入力とし、計測データの実測値ｙkを出力とする線形時不変システムとみなせるので、以下の手順で周波数特性の算出式を導出することができる。以下では、修正された自己回帰モデルを線形時不変モデルとみなしたものを、単にシステムと称することにする。予測誤差ｘkは、以下の式１４で表される。
次に、求めた係数αで特定される修正された自己回帰モデルの周波数特性を求める方法を説明する。
式１で与えられる計測データの予測値ｙ＾kの予測誤差ｘkが白色雑音になるという性質を使うと、修正された自己回帰モデルは白色雑音である予測誤差ｘkを入力とし、計測データの実測値ｙkを出力とする線形時不変システムとみなせるので、以下の手順で周波数特性の算出式を導出することができる。以下では、修正された自己回帰モデルを線形時不変モデルとみなしたものを、単にシステムと称することにする。予測誤差ｘkは、以下の式１４で表される。
式１４の両辺にｚ変換を施すと、以下の式１５が得られる。
式１５から、システムのインパルス応答のｚ変換である伝達関数Ｈ（ｚ）は、以下の式１６のように求まる。
システムの周波数特性は、正弦波入力に対する出力の振幅と位相との変化として現れ、インパルス応答のフーリエ変換で求められる。言い換えると、ｚが複素平面の単位円上を回転した時の伝達関数Ｈ（ｚ）が周波数特性となる。ここで、式１６におけるｚを、以下の式１７のように置くことを考える。
その場合、Ｈ（ｚ）の振幅特性（システムの周波数特性）は、以下の式１８のように表すことができる。式１８におけるωＴは、０〜２πの範囲で変化するものとする。
本実験では、得られた計測データ毎に、式５と式７を用いて自己相関行列Ｒを求め、式８で表される特異値分解を行うことによって自己相関行列Ｒの固有値を求め、自己相関行列Ｒの固有値の分布を求めた。また、得られた計測データ毎に、計測データの波形を求めた。また、得られた計測データ毎に、上記手順で求められた自己相関行列Ｒの特異値分解から式１１を用いて行列Σsと行列Ｕsとを求め、式５と式１３を用いて修正された自己回帰モデルの係数αを求め、求めた係数αで特定される修正された自己回帰モデルの波形を求めた。そして、得られた計測データ毎に、式１８を用いて、修正された自己回帰モデルの周波数特性等を求めた。
求められた実験結果を、図３〜８を用いて説明する。
求められた実験結果を、図３〜８を用いて説明する。
図３は、計測データそれぞれについて、上記手順で求められた自己相関行列Ｒの固有値の分布を示す図である。
図３（ａ）のグラフは、軸受１０６に疵がない正常な軸受を用いた場合に計測された計測データを利用した際の自己相関行列Ｒの固有値の分布を示すグラフである。図３（ｂ）のグラフは、軸受１０６に内輪に疵がある軸受を用いた場合に計測された計測データを利用した際の自己相関行列Ｒの固有値の分布を示すグラフである。図３（ｃ）のグラフは、軸受１０６に保持器に疵がある軸受を用いた場合に計測された計測データを利用した際の自己相関行列Ｒの固有値の分布を示すグラフである。図３（ｄ）のグラフは、軸受１０６に転動体に疵がある軸受を用いた場合に計測された計測データを利用した際の自己相関行列Ｒの固有値の分布を示すグラフである。図３（ｅ）のグラフは、軸受１０６に外輪に疵がある軸受を用いた場合に計測された計測データを利用した際の自己相関行列Ｒの固有値の分布を示すグラフである。
図３の各グラフは、自己相関行列Ｒを特異値分解して得られた固有値σ11〜σmmを昇順に並べ替え、プロットしたグラフであり、横軸が固有値のインデックス、縦軸が固有値の値を示す。
図３（ａ）のグラフは、軸受１０６に疵がない正常な軸受を用いた場合に計測された計測データを利用した際の自己相関行列Ｒの固有値の分布を示すグラフである。図３（ｂ）のグラフは、軸受１０６に内輪に疵がある軸受を用いた場合に計測された計測データを利用した際の自己相関行列Ｒの固有値の分布を示すグラフである。図３（ｃ）のグラフは、軸受１０６に保持器に疵がある軸受を用いた場合に計測された計測データを利用した際の自己相関行列Ｒの固有値の分布を示すグラフである。図３（ｄ）のグラフは、軸受１０６に転動体に疵がある軸受を用いた場合に計測された計測データを利用した際の自己相関行列Ｒの固有値の分布を示すグラフである。図３（ｅ）のグラフは、軸受１０６に外輪に疵がある軸受を用いた場合に計測された計測データを利用した際の自己相関行列Ｒの固有値の分布を示すグラフである。
図３の各グラフは、自己相関行列Ｒを特異値分解して得られた固有値σ11〜σmmを昇順に並べ替え、プロットしたグラフであり、横軸が固有値のインデックス、縦軸が固有値の値を示す。
図３（ａ）を見ると、他よりも顕著に高い値をもつ固有値が５つあるのが分かる。その５つの固有値の中でも特に２つの固有値が他の３つよりも高い値となっていることが分かる。
何れの場合にも、自己相関行列Ｒの固有値のうち値が他よりも顕著に高い固有値の数は、固有値全体の数よりも顕著に少ないことが分かる。
何れの場合にも、自己相関行列Ｒの固有値のうち値が他よりも顕著に高い固有値の数は、固有値全体の数よりも顕著に少ないことが分かる。
また、本実験では、得られた計測データ毎に、計測データの波形を求めた。
図４は、計測データそれぞれの波形を示す図である。
図４（ａ）のグラフは、軸受１０６に疵がない正常な軸受を用いた場合に計測された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの実績値を示す。図４（ｂ）のグラフは、軸受１０６に内輪に疵がある軸受を用いた場合に計測された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの実績値を示す。図４（ｃ）のグラフは、軸受１０６に保持器に疵がある軸受を用いた場合に計測された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの実績値を示す。図４（ｄ）のグラフは、軸受１０６に転動体に疵がある軸受を用いた場合に計測された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの実績値を示す。図４（ｅ）のグラフは、軸受１０６に外輪に疵がある軸受を用いた場合に計測された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの実績値を示す。 図４の各グラフを見ると、例えば、図４（ａ）及び図４（ｃ）が似通った波形をしているのが分かる。また、例えば、図４（ｄ）及び図４（ｅ）が似通った波形をしているのも分かる。そのため、図４のグラフから軸受１０６に生じた異常を診断することは困難である。
図４は、計測データそれぞれの波形を示す図である。
図４（ａ）のグラフは、軸受１０６に疵がない正常な軸受を用いた場合に計測された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの実績値を示す。図４（ｂ）のグラフは、軸受１０６に内輪に疵がある軸受を用いた場合に計測された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの実績値を示す。図４（ｃ）のグラフは、軸受１０６に保持器に疵がある軸受を用いた場合に計測された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの実績値を示す。図４（ｄ）のグラフは、軸受１０６に転動体に疵がある軸受を用いた場合に計測された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの実績値を示す。図４（ｅ）のグラフは、軸受１０６に外輪に疵がある軸受を用いた場合に計測された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの実績値を示す。 図４の各グラフを見ると、例えば、図４（ａ）及び図４（ｃ）が似通った波形をしているのが分かる。また、例えば、図４（ｄ）及び図４（ｅ）が似通った波形をしているのも分かる。そのため、図４のグラフから軸受１０６に生じた異常を診断することは困難である。
図５は、使用固有値数ｓを１として求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示す図である。
図５（ａ）のグラフは、使用固有値数ｓを１として、軸受１０６に疵がない正常な軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図５（ｂ）のグラフは、使用固有値数ｓを１として、軸受１０６に内輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図５（ｃ）のグラフは、使用固有値数ｓを１として、軸受１０６に保持器に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図５（ｄ）のグラフは、使用固有値数ｓを１として、軸受１０６に転動体に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図５（ｅ）のグラフは、使用固有値数ｓを１として、軸受１０６に外輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。
図５の各グラフは、図４の各グラフと比べると、ノイズ成分が低減されていることが分かる。
図５（ａ）のグラフは、使用固有値数ｓを１として、軸受１０６に疵がない正常な軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図５（ｂ）のグラフは、使用固有値数ｓを１として、軸受１０６に内輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図５（ｃ）のグラフは、使用固有値数ｓを１として、軸受１０６に保持器に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図５（ｄ）のグラフは、使用固有値数ｓを１として、軸受１０６に転動体に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図５（ｅ）のグラフは、使用固有値数ｓを１として、軸受１０６に外輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。
図５の各グラフは、図４の各グラフと比べると、ノイズ成分が低減されていることが分かる。
図６は、使用固有値数ｓを１として求められた係数αで特定される修正された自己回帰モデルの周波数特性を示す図である。図６の各グラフは、図５の各波形についての周波数特性を示すグラフである。
図６（ａ）のグラフは、使用固有値数ｓを１として、軸受１０６に疵がない正常な軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフであり、横軸が周波数、縦軸が信号強度を示す。図６（ｂ）のグラフは、使用固有値数ｓを１として、軸受１０６に内輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフであり、横軸が周波数、縦軸が信号強度を示す。図６（ｃ）のグラフは、使用固有値数ｓを１として、軸受１０６に保持器に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフであり、横軸が周波数、縦軸が信号強度を示す。図６（ｄ）のグラフは、使用固有値数ｓを１として、軸受１０６に転動体に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフであり、横軸が周波数、縦軸が信号強度を示す。図６（ｅ）のグラフは、使用固有値数ｓを１として、軸受１０６に外輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフであり、横軸が周波数、縦軸が信号強度を示す。
図６（ａ）のグラフは、使用固有値数ｓを１として、軸受１０６に疵がない正常な軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフであり、横軸が周波数、縦軸が信号強度を示す。図６（ｂ）のグラフは、使用固有値数ｓを１として、軸受１０６に内輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフであり、横軸が周波数、縦軸が信号強度を示す。図６（ｃ）のグラフは、使用固有値数ｓを１として、軸受１０６に保持器に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフであり、横軸が周波数、縦軸が信号強度を示す。図６（ｄ）のグラフは、使用固有値数ｓを１として、軸受１０６に転動体に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフであり、横軸が周波数、縦軸が信号強度を示す。図６（ｅ）のグラフは、使用固有値数ｓを１として、軸受１０６に外輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフであり、横軸が周波数、縦軸が信号強度を示す。
図６（ａ）のグラフを見ると、周波数８２４Ｈｚの部分にピークが立っているのが分かる。また、図６（ｂ）のグラフを見ると、周波数１６４９Ｈｚの部分にピークが立っているのが分かる。また、図６（ｃ）のグラフを見ると、周波数１６４６Ｈｚの部分にピークが立っているのが分かる。
図６（ｄ）のグラフを見ると、周波数１９８５３Ｈｚの部分にピークが立っているのが分かる。また、図６（ｅ）のグラフを見ると、周波数２３００７Ｈｚの部分にピークが立っているのが分かる。
図６の各グラフにおけるピークが立っている周波数を見ると、軸受１０６に疵がない場合は、軸受１０６に疵がある場合に比べて小さい周波数でピークが立っていることが分かる。即ち、使用固有値数ｓを１として、求められた係数αで特定される修正された自己回帰モデルの周波数特性から、軸受１０６における疵の有無を判断することが可能であることが分かる。また、図６（ｄ）よりも、図６（ｅ）の方が高い周波数であることが分かる。また、図６（ｄ）、（ｅ）の各グラフにおけるピークが立っている周波数は、図６（ａ）〜（ｃ）の各グラフにおけるピークが立っている周波数よりも顕著に高いことが分かる。そのため、使用固有値数ｓを１として、求められた係数αで特定される修正された自己回帰モデルの周波数特性から、軸受１０６の転動体又は外輪の疵の有無を判断することが可能であることが分かる。
以上より、自己相関行列Ｒの１個の固有値に基づいて求められた係数αで特定される修正された自己回帰モデルの周波数特性に基づいて、軸受１０６の異常診断を行うことができるという知見を得ることができた。
図６（ｄ）のグラフを見ると、周波数１９８５３Ｈｚの部分にピークが立っているのが分かる。また、図６（ｅ）のグラフを見ると、周波数２３００７Ｈｚの部分にピークが立っているのが分かる。
図６の各グラフにおけるピークが立っている周波数を見ると、軸受１０６に疵がない場合は、軸受１０６に疵がある場合に比べて小さい周波数でピークが立っていることが分かる。即ち、使用固有値数ｓを１として、求められた係数αで特定される修正された自己回帰モデルの周波数特性から、軸受１０６における疵の有無を判断することが可能であることが分かる。また、図６（ｄ）よりも、図６（ｅ）の方が高い周波数であることが分かる。また、図６（ｄ）、（ｅ）の各グラフにおけるピークが立っている周波数は、図６（ａ）〜（ｃ）の各グラフにおけるピークが立っている周波数よりも顕著に高いことが分かる。そのため、使用固有値数ｓを１として、求められた係数αで特定される修正された自己回帰モデルの周波数特性から、軸受１０６の転動体又は外輪の疵の有無を判断することが可能であることが分かる。
以上より、自己相関行列Ｒの１個の固有値に基づいて求められた係数αで特定される修正された自己回帰モデルの周波数特性に基づいて、軸受１０６の異常診断を行うことができるという知見を得ることができた。
しかし、図６（ｂ）と図６（ｃ）とのグラフに示されるように、軸受１０６の内輪に疵がある場合と、軸受１０６の保持器に疵がある場合と、では、似通った値の周波数において周波数特性のピークが立っている。そのため、使用固有値数ｓを１として、求められた係数αで特定される修正された自己回帰モデルの周波数特性から、軸受１０６に内輪の又は保持器の疵の有無を判断することは可能であるが、内輪の疵であるか、保持器の疵であるか、両方の疵であるか、を判断することは困難である。
そこで、本実験では、更に、使用固有値数ｓの値を５に増やして、自己相関行列Ｒの固有値のうち、最大のものから５つを用いて係数αを求め、求めた係数αで特定される修正された自己回帰モデルの波形、周波数特性を求めた。
そこで、本実験では、更に、使用固有値数ｓの値を５に増やして、自己相関行列Ｒの固有値のうち、最大のものから５つを用いて係数αを求め、求めた係数αで特定される修正された自己回帰モデルの波形、周波数特性を求めた。
図７は、使用固有値数ｓを５として求められた係数αで特定される修正された自己回帰モデルの波形を示す図である。
図７（ａ）のグラフは、使用固有値数ｓを５として、軸受１０６に疵がない正常な軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図７（ｂ）のグラフは、使用固有値数ｓを５として、軸受１０６に内輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図７（ｃ）のグラフは、使用固有値数ｓを５として、軸受１０６に保持器に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図７（ｄ）のグラフは、使用固有値数ｓを５として、軸受１０６に転動体に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図７（ｅ）のグラフは、使用固有値数ｓを５として、軸受１０６に外輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。
図７（ａ）のグラフは、使用固有値数ｓを５として、軸受１０６に疵がない正常な軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図７（ｂ）のグラフは、使用固有値数ｓを５として、軸受１０６に内輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図７（ｃ）のグラフは、使用固有値数ｓを５として、軸受１０６に保持器に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図７（ｄ）のグラフは、使用固有値数ｓを５として、軸受１０６に転動体に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図７（ｅ）のグラフは、使用固有値数ｓを５として、軸受１０６に外輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。
図８は、使用固有値数ｓを５として求められた係数αで特定される修正された自己回帰モデルの周波数特性を示す図である。図８の各グラフは、図７の各波形についての周波数特性を示すグラフである。
図８（ａ）のグラフは、使用固有値数ｓを５として、軸受１０６に疵がない正常な軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフであり、横軸が周波数、縦軸が信号強度を示す。図８（ｂ）のグラフは、使用固有値数ｓを５として、軸受１０６に内輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフであり、横軸が周波数、縦軸が信号強度を示す。図８（ｃ）のグラフは、使用固有値数ｓを５として、軸受１０６に保持器に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフであり、横軸が周波数、縦軸が信号強度を示す。図８（ｄ）のグラフは、使用固有値数ｓを５として、軸受１０６に転動体に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフであり、横軸が周波数、縦軸が信号強度を示す。図８（ｅ）のグラフは、使用固有値数ｓを５として、軸受１０６に外輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフであり、横軸が周波数、縦軸が信号強度を示す。
図８（ａ）のグラフは、使用固有値数ｓを５として、軸受１０６に疵がない正常な軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフであり、横軸が周波数、縦軸が信号強度を示す。図８（ｂ）のグラフは、使用固有値数ｓを５として、軸受１０６に内輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフであり、横軸が周波数、縦軸が信号強度を示す。図８（ｃ）のグラフは、使用固有値数ｓを５として、軸受１０６に保持器に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフであり、横軸が周波数、縦軸が信号強度を示す。図８（ｄ）のグラフは、使用固有値数ｓを５として、軸受１０６に転動体に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフであり、横軸が周波数、縦軸が信号強度を示す。図８（ｅ）のグラフは、使用固有値数ｓを５として、軸受１０６に外輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフであり、横軸が周波数、縦軸が信号強度を示す。
図８（ａ）のグラフを見ると、周波数８２１Ｈｚと１０４７Ｈｚの部分にピークが立っているのが分かる。また、図８（ｂ）のグラフを見ると、周波数１６４８Ｈｚ、２０２１Ｈｚ、６４７４Ｈｚの部分にピークが立っているのが分かる。また、図８（ｃ）のグラフを見ると、周波数１６４６Ｈｚ、３２９１Ｈｚ、１７４６Ｈｚの部分にピークが立っているのが分かる。図８（ｄ）のグラフを見ると、周波数１９８５３Ｈｚの部分にピークが立っているのが分かる。また、図８（ｅ）のグラフを見ると、周波数２２７８５Ｈｚ、２３０２０Ｈｚの部分にピークが立っているのが分かる。
図８の各グラフにおけるピークと、図６の各グラフにおけるピークと、を見比べてみると、図６のグラフでピークが立っている周波数近辺においては、図８のグラフでもピークが立っているのが分かる。図８のグラフでは、更に、図６のグラフで見られなかったピークが見られるようになっている。
このように、係数αを求めるのに利用される自己相関行列Ｒの固有値を増やすことで、修正された自己回帰モデルによる予測値において、図６でピークが立っていた周波数以外の周波数の成分の信号強度が増加したことが分かる。即ち、係数αを求めるのに利用される自己相関行列Ｒの固有値は、求められた係数αで特定される修正された自己回帰モデルによる計測データの予測値の時間波形に含まれる各周波数の成分の強度と相関があることが見出された。
また、図８（ｂ）、図８（ｃ）のグラフを見てみると、図６（ｂ）、図６（ｃ）のグラフと同様に、１６４０Ｈｚ付近にピークが立っているのが分かる。しかし、図８（ｂ）では、更に、２０２１Ｈｚ、６４７４Ｈｚにおいてもピークが立っており、図８（ｃ）では、更に、３２９１Ｈｚ、１７４６Ｈｚにおいてもピークが立っていることが分かる。このように、図８では、内輪に疵のある軸受１０６と保持器に疵のある軸受１０６とで、ピークが立っている周波数に顕著な違いが見て取れる。即ち、使用固有値数ｓを５にすることで、軸受の疵が、内輪の疵であるか、保持器の疵であるか、両方の疵であるか、を見分けることができる。
図８の各グラフにおけるピークと、図６の各グラフにおけるピークと、を見比べてみると、図６のグラフでピークが立っている周波数近辺においては、図８のグラフでもピークが立っているのが分かる。図８のグラフでは、更に、図６のグラフで見られなかったピークが見られるようになっている。
このように、係数αを求めるのに利用される自己相関行列Ｒの固有値を増やすことで、修正された自己回帰モデルによる予測値において、図６でピークが立っていた周波数以外の周波数の成分の信号強度が増加したことが分かる。即ち、係数αを求めるのに利用される自己相関行列Ｒの固有値は、求められた係数αで特定される修正された自己回帰モデルによる計測データの予測値の時間波形に含まれる各周波数の成分の強度と相関があることが見出された。
また、図８（ｂ）、図８（ｃ）のグラフを見てみると、図６（ｂ）、図６（ｃ）のグラフと同様に、１６４０Ｈｚ付近にピークが立っているのが分かる。しかし、図８（ｂ）では、更に、２０２１Ｈｚ、６４７４Ｈｚにおいてもピークが立っており、図８（ｃ）では、更に、３２９１Ｈｚ、１７４６Ｈｚにおいてもピークが立っていることが分かる。このように、図８では、内輪に疵のある軸受１０６と保持器に疵のある軸受１０６とで、ピークが立っている周波数に顕著な違いが見て取れる。即ち、使用固有値数ｓを５にすることで、軸受の疵が、内輪の疵であるか、保持器の疵であるか、両方の疵であるか、を見分けることができる。
このように本発明者らは、使用固有値数ｓが１の場合には見分けることが困難であった疵の存在する部位を、使用固有値数ｓを増やすことで、見分けることができるようになるという知見を得た。一方で、本発明者らは、使用固有値数ｓを増やし過ぎると修正された自己回帰モデルは、実際の計測データに近づいていき、軸受１０６の異常診断に有用でないノイズ成分等についても信号強度が増加されてしまい、式１８で得られる周波数特性にもノイズ成分に起因するピークが生じるという知見を得た。即ち、自己相関行列Ｒの固有値のうち、相対的に値の高い一部の固有値を用いて係数αを求めることで、計測データに含まれる成分のうち、物体の異常診断に有用な成分を近似した修正された自己回帰モデルを求めることができると考えられる。
以上より、係数αを求めるのに利用される自己相関行列Ｒの固有値の数である使用固有値ｓを調整することで、より精度よく、軸受１０６の異常診断を行うことができるという知見を得ることができた。
以上より、係数αを求めるのに利用される自己相関行列Ｒの固有値の数である使用固有値ｓを調整することで、より精度よく、軸受１０６の異常診断を行うことができるという知見を得ることができた。
また、本実験では、軸受１０６が疵のない正常な状態である場合、内輪に疵がある状態である場合、保持器に疵がある状態である場合、転動体に疵がある状態である場合、外輪に疵がある状態である場合それぞれについて、以下の作業をｑ回行った。本実験では、ｑは５０としたが、３０や１００等の他の値としてもよい。即ち、軸受１０６について、連続したＭ個の計測データを取得し、取得した連続したＭ個の計測データ毎に、式１１を用いて行列Σsと行列Ｕsとを求め、式５と式１３を用いて修正された自己回帰モデルの係数αを求めた。使用固有値数ｓは、３とした。
これにより、軸受１０６が正常な状態である場合、内輪に疵がある状態である場合、保持器に疵がある状態である場合、転動体に疵がある状態である場合、外輪に疵がある状態である場合の５つの場合それぞれについて、ｑ個の修正された自己回帰モデルの係数αを取得した。
軸受１０６が取り得る状態は、疵のない正常な状態、内輪に疵がある状態、保持器に疵がある状態、転動体に疵がある状態、外輪に疵がある状態の５つの状態であるとする。以下では、軸受１０６が取り得る状態の数をｐとする。本実験では、ｐは５となる。
これにより、軸受１０６が正常な状態である場合、内輪に疵がある状態である場合、保持器に疵がある状態である場合、転動体に疵がある状態である場合、外輪に疵がある状態である場合の５つの場合それぞれについて、ｑ個の修正された自己回帰モデルの係数αを取得した。
軸受１０６が取り得る状態は、疵のない正常な状態、内輪に疵がある状態、保持器に疵がある状態、転動体に疵がある状態、外輪に疵がある状態の５つの状態であるとする。以下では、軸受１０６が取り得る状態の数をｐとする。本実験では、ｐは５となる。
図９は、軸受１０６が正常な状態である場合、内輪に疵がある状態である場合、保持器に疵がある状態である場合、転動体に疵がある状態である場合、軸受１０６の外輪に疵がある場合それぞれについての、上記手順で求められた係数αのパターンを示す図である。係数αのパターンとは、係数αの成分を並べた際の形状を示す情報である。本実験では、ｍは５００であるため、係数αは、α1〜α500の５００個の成分から成る。図９に示されるパターンそれぞれは、αi（１＜＝ｉ＜＝５００）をインデックスｉ順に並べたパターンである。このインデックスｉは、式１において係数の成分がかかる計測データｙの計測時点を示すインデックスとみなすことができる。即ち、図９に示されるパターンそれぞれは、係数αの各成分を、式１において各成分がかかる計測データｙの計測時点に基づいて並べたパターンの一例である。
係数αのパターンには、例えば、係数αの成分をインデックス順に並べたベクトル、係数αの成分をインデックス順に予め定められた数（例えば、１等）ずつとばしに並べたベクトル、係数αの成分をインデックスと逆順に並べたベクトル等がある。また、係数αのパターンには、例えば、各成分に同じ値が加算された係数αの成分をインデックス順に並べたベクトル、各成分に同じ値が加算された係数αの成分をインデックス順に予め定められた数（例えば、１等）ずつとばしに並べたベクトル、各成分に同じ値が加算された係数αの成分をインデックスと逆順に並べたベクトル等がある。また、係数αのパターンには、例えば、各成分に同じ値が掛けられた係数αの成分をインデックス順に並べたベクトル、各成分に同じ値が掛けられた係数αの成分をインデックス順に予め定められた数（例えば、１等）ずつとばしに並べたベクトル、各成分に同じ値が掛けられた係数αの成分をインデックスと逆順に並べたベクトル等がある。また、係数αの成分を並べた際の形状を示す情報である係数αのパターンは、ベクトルでなく、係数αの成分を並べたグラフの画像であるとしてもよい。
係数αのパターンには、例えば、係数αの成分をインデックス順に並べたベクトル、係数αの成分をインデックス順に予め定められた数（例えば、１等）ずつとばしに並べたベクトル、係数αの成分をインデックスと逆順に並べたベクトル等がある。また、係数αのパターンには、例えば、各成分に同じ値が加算された係数αの成分をインデックス順に並べたベクトル、各成分に同じ値が加算された係数αの成分をインデックス順に予め定められた数（例えば、１等）ずつとばしに並べたベクトル、各成分に同じ値が加算された係数αの成分をインデックスと逆順に並べたベクトル等がある。また、係数αのパターンには、例えば、各成分に同じ値が掛けられた係数αの成分をインデックス順に並べたベクトル、各成分に同じ値が掛けられた係数αの成分をインデックス順に予め定められた数（例えば、１等）ずつとばしに並べたベクトル、各成分に同じ値が掛けられた係数αの成分をインデックスと逆順に並べたベクトル等がある。また、係数αの成分を並べた際の形状を示す情報である係数αのパターンは、ベクトルでなく、係数αの成分を並べたグラフの画像であるとしてもよい。
図９に示されるパターンを見てみると、以下のようなことが見てとれる。軸受１０６の外輪に疵がある場合、ｑ個の係数αのパターンそれぞれは、上下対称な２つの波線で囲まれた領域を示すようなパターンとなっている。また、軸受１０６の転動体に疵がある場合、ｑ個の係数αのパターンそれぞれは、複数の波線が重なり合っているようなパターンとなっている。また、軸受１０６の内輪に疵がある場合、ｑ個の係数αのパターンそれぞれは、４つのピークのある山形のパターンが連続しているパターンとなっている。また、軸受１０６に疵がない場合、ｑ個の係数αのパターンそれぞれは、徐々に減衰していく波型のパターンとなっている。また、軸受１０６の保持器に疵がある場合、ｑ個の係数αのパターンそれぞれは、２つのピークのある山形のパターンが連続しているパターンとなっている。
このように、軸受１０６が正常な場合、軸受１０６の内輪に疵がある場合、軸受１０６の保持器に疵がある場合、軸受１０６の転動体に疵がある場合、軸受１０６の外輪に疵がある場合それぞれについて、修正された自己回帰モデルの係数αが、特徴的なパターンを示すことが見て取れる。
そこで、発明者らは、修正された自己回帰モデルの係数αのパターンが、どのようなパターンを示すのかを特定することで、軸受１０６の診断を行うことができるという知見を得た。
このように、軸受１０６が正常な場合、軸受１０６の内輪に疵がある場合、軸受１０６の保持器に疵がある場合、軸受１０６の転動体に疵がある場合、軸受１０６の外輪に疵がある場合それぞれについて、修正された自己回帰モデルの係数αが、特徴的なパターンを示すことが見て取れる。
そこで、発明者らは、修正された自己回帰モデルの係数αのパターンが、どのようなパターンを示すのかを特定することで、軸受１０６の診断を行うことができるという知見を得た。
また、本実験では、以下のようにして、軸受１０６が正常な場合、軸受１０６の内輪に疵がある場合、軸受１０６の保持器に疵がある場合、軸受１０６の転動体に疵がある場合、軸受１０６の外輪に疵がある場合それぞれについて、修正された自己回帰モデルの係数αの代表的なパターンを学習により特定した。
まず、各列が図９に示される係数αそれぞれを示すｍ×ｄのサイズの行列Ｙを、以下の式１９のように生成する。ｄは、ｐ×ｑであり、本実験では２５０となる。
まず、各列が図９に示される係数αそれぞれを示すｍ×ｄのサイズの行列Ｙを、以下の式１９のように生成する。ｄは、ｐ×ｑであり、本実験では２５０となる。
式１９のαk、(i,j)は、ｉに対応する状態の軸受１０６から計測されたｑ個の計測データのうちｊ番目の計測データから決定された係数αの第ｋ成分αkを示す。本実験では、外輪に疵がある状態の軸受１０６から計測された計測データから決定された係数αの成分は、式１９におけるαk、(0,j)（１＜＝ｋ＜＝ｍ、１＜＝ｊ＜＝ｑ）である。また、転動体に疵がある状態の軸受１０６から計測された計測データから決定された係数αの成分は、式１９におけるαk、(1,j)（１＜＝ｋ＜＝ｍ、１＜＝ｊ＜＝ｑ）である。また、内輪に疵がある状態の軸受１０６から計測された計測データから決定された係数αの成分は、式１９におけるαk、(2,j)（１＜＝ｋ＜＝ｍ、１＜＝ｊ＜＝ｑ）である。また、疵がない正常な状態の軸受１０６から計測された計測データから決定された係数αの成分は、式１９におけるαk、(3,j)（１＜＝ｋ＜＝ｍ、１＜＝ｊ＜＝ｑ）である。また、保持器に疵がある状態の軸受１０６から計測された計測データから決定された係数αの成分は、式１９におけるαk、(4,j)（１＜＝ｋ＜＝ｍ、１＜＝ｊ＜＝ｑ）である。
本実験では、この行列Ｙを非負値行列因子分解することで、軸受１０６の状態それぞれについて、修正された自己回帰モデルの係数αの代表的なパターンを特定する。行列Ｙを非負値行列因子分解するためには、行列Ｙの各列の係数αのパターンを変えずに、行列Ｙ自体の各成分を、非負値にするよう調整する必要がある。そこで、例えば、行列Ｙの０未満の成分のうち、最小のものの絶対値を、行列Ｙの各成分に加算することで、行列Ｙの各成分が非負値となるようにする。以下では、各成分が非負値である行列を、非負値行列とする。
以下の式２０に示すように、この行列Ｙを、非負値行列である行列Ａと行列Ｐとに非負値行列因子分解する。行列Ａは、ｍ×ｐのサイズの行列であり、行列Ｐは、ｐ×ｄのサイズの行列である。行列Ａの列数は、非負値行列因子分解における基底の数であり、本実験では、軸受１０６が取り得る状態の数としてｐ（５）とする。これにより、行列Ａの各列には、行列Ｙから導出されたｐ（５）個の係数αのパターン（基底ともいう）が格納されることとなる。
以下の式２０に示すように、この行列Ｙを、非負値行列である行列Ａと行列Ｐとに非負値行列因子分解する。行列Ａは、ｍ×ｐのサイズの行列であり、行列Ｐは、ｐ×ｄのサイズの行列である。行列Ａの列数は、非負値行列因子分解における基底の数であり、本実験では、軸受１０６が取り得る状態の数としてｐ（５）とする。これにより、行列Ａの各列には、行列Ｙから導出されたｐ（５）個の係数αのパターン（基底ともいう）が格納されることとなる。
以下の式２１のようにコスト関数Ｊを定義する。
式２１の右辺第１項は、行列Ａと行列Ｐとの積と、行列Ｙと、が最小二乗的意味で一致するための条件を示す。式２１の右辺第１項の右下のＦ記号は、フロベニウスノルムであることを示す記号である。即ち、式２１の右辺第１項内の｜｜ＡＰ−Ｙ｜｜は、フロベニウスノルムであり、行列（ＡＰ−Ｙ）の各成分の絶対値の２乗の合計の平方根となる。式２１の右辺第２項は、行列Ｐのスパース性を要求する条件を示す。この条件を付加することによって、非負値行列因子分解が過学習となることを防ぐことができる。
式２１の右辺第２項におけるインデックスｉは、行列の行を示すインデックスである。また、式２１の右辺第２項におけるインデックスｋは、行列の列を示すインデックスである。式２１の右辺第３項は、行列Ａの異なる列に格納されるパターン同士の直交性を近似的に要求する条件を示し、行列Ａのｐ（５）個の列に格納されるパターンが可能な限り線形従属となることを防ぐ効果がある。行列Ａの列に格納されるパターンが線形従属となれば、あるパターンを行列Ａの列に格納されるパターンの線形結合で近似するときの係数（重みともいう）が一通りではなくなるため、パターンの同定が難しくなる。また、式２１の右辺第３項における行列ＡTは、行列Ａの転置行列であり、行列Ｉは、単位行列である。また、式２１の右辺第３項におけるインデックスｋ、ｋ-（式中では、上にバーが引かれたｋとして表現）は、それぞれ行列の行と、列と、を示すインデックスである。
ここで、コスト関数Ｊとして、式２１の右辺第１項のみを用いても、非負値行列因子分解は機能するので、コスト関数Ｊとして、式２１の右辺第１項のみを用いてもよい。また、上述の過学習や線形従属といった問題が生じる場合に、コスト関数Ｊとして、式２１の右辺第１項に式２１の右辺第２項、第３項を追加した式を用いることとしてもよい。そして、右辺第２項、第３項を追加する場合、軸受１０６の各状態と、図１４等で後述する方法で決定される学習パターンと、を良好に対応付けることが出来るように、例えば、ケーススタディにより係数βと係数γとを予め決定しておく。本実施形態では、βを０．０３５、γを０とした。
このコスト関数Ｊを最小化するように、行列Ｙを、行列Ａと行列Ｐと分解する。即ち、以下の式２２で定義されるコスト関数Ｊの最小化問題を解くこととなる。
式２１の右辺第２項におけるインデックスｉは、行列の行を示すインデックスである。また、式２１の右辺第２項におけるインデックスｋは、行列の列を示すインデックスである。式２１の右辺第３項は、行列Ａの異なる列に格納されるパターン同士の直交性を近似的に要求する条件を示し、行列Ａのｐ（５）個の列に格納されるパターンが可能な限り線形従属となることを防ぐ効果がある。行列Ａの列に格納されるパターンが線形従属となれば、あるパターンを行列Ａの列に格納されるパターンの線形結合で近似するときの係数（重みともいう）が一通りではなくなるため、パターンの同定が難しくなる。また、式２１の右辺第３項における行列ＡTは、行列Ａの転置行列であり、行列Ｉは、単位行列である。また、式２１の右辺第３項におけるインデックスｋ、ｋ-（式中では、上にバーが引かれたｋとして表現）は、それぞれ行列の行と、列と、を示すインデックスである。
ここで、コスト関数Ｊとして、式２１の右辺第１項のみを用いても、非負値行列因子分解は機能するので、コスト関数Ｊとして、式２１の右辺第１項のみを用いてもよい。また、上述の過学習や線形従属といった問題が生じる場合に、コスト関数Ｊとして、式２１の右辺第１項に式２１の右辺第２項、第３項を追加した式を用いることとしてもよい。そして、右辺第２項、第３項を追加する場合、軸受１０６の各状態と、図１４等で後述する方法で決定される学習パターンと、を良好に対応付けることが出来るように、例えば、ケーススタディにより係数βと係数γとを予め決定しておく。本実施形態では、βを０．０３５、γを０とした。
このコスト関数Ｊを最小化するように、行列Ｙを、行列Ａと行列Ｐと分解する。即ち、以下の式２２で定義されるコスト関数Ｊの最小化問題を解くこととなる。
式２２における第１の制約条件は、行列Ａと行列Ｐとが非負値行列であることを示す条件である。式２２における第２の制約条件は、行列Ａを基準化（ｎｏｒｍａｌｉｚｅ）するための条件である。また、式２２の第２の制約条件におけるインデックスｋは、行列の行と、列と、を示すインデックスである。
本実験では、コスト関数Ｊを最小化するために、以下のような処理を行う。即ち、以下の式２３と式２４とを用いて行列Ａと行列Ｐとの更新を繰り返すことで、コスト関数Ｊを最小化する。また、コスト関数Ｊを最小化するために、以下に記載する方法以外の方法（例えば、最急降下法、ニュートン法、確率的勾配降下法等の最適化手法）を用いてもよい。
本実験では、コスト関数Ｊを最小化するために、以下のような処理を行う。即ち、以下の式２３と式２４とを用いて行列Ａと行列Ｐとの更新を繰り返すことで、コスト関数Ｊを最小化する。また、コスト関数Ｊを最小化するために、以下に記載する方法以外の方法（例えば、最急降下法、ニュートン法、確率的勾配降下法等の最適化手法）を用いてもよい。
式２３のηPは、行列Ｐの更新に用いられる緩和係数である。式２４のηAは、行列Ａの更新に用いられる緩和係数である。
式２３におけるコスト関数Ｊの行列Ｐによる偏微分は、以下の式２５の通りとなる。
式２３におけるコスト関数Ｊの行列Ｐによる偏微分は、以下の式２５の通りとなる。
式２５の行列Ｆは、サイズがｐ×ｄの全ての成分が１の行列である。
式２４におけるコスト関数Ｊの行列Ａによる偏微分は、以下の式２６の通りとなる。
式２４におけるコスト関数Ｊの行列Ａによる偏微分は、以下の式２６の通りとなる。
式２６の行列Ｅは、サイズがｐ×ｐの全ての成分が１の行列である。
式２５を用いて、式２３を書き換えると以下の式２７が求まる。
式２５を用いて、式２３を書き換えると以下の式２７が求まる。
式２７のインデックスｉは、行列の行を示すインデックスである。式２７のインデックスｋは、行列の列を示すインデックスである。
行列Ｐは、非負値であるため、式２７の値は、非負値である必要がある。式２７の右辺の第３項は、必ず非負値となる項である。即ち、式２７の右辺から第３項を除いた式の値が０であるならば、必ず式２７は、非負値となる。そのため、以下の式２８に示すように、式２７の右辺から第３項を除いた式が０であるという条件を満たすように、行列Ｐの更新を行うこととする。
行列Ｐは、非負値であるため、式２７の値は、非負値である必要がある。式２７の右辺の第３項は、必ず非負値となる項である。即ち、式２７の右辺から第３項を除いた式の値が０であるならば、必ず式２７は、非負値となる。そのため、以下の式２８に示すように、式２７の右辺から第３項を除いた式が０であるという条件を満たすように、行列Ｐの更新を行うこととする。
式２８のインデックスｉは、行列の行を示すインデックスである。式２８のインデックスｋは、行列の列を示すインデックスである。
式２８から、式２３における緩和係数ηPの値が、以下の式２９のように求まる。
式２８から、式２３における緩和係数ηPの値が、以下の式２９のように求まる。
式２９のインデックスｉは、行列の行を示すインデックスである。式２９のインデックスｋは、行列の列を示すインデックスである。
式２９を用いて、式２３を書き換えると、以下の式３０のようになる。
式２９を用いて、式２３を書き換えると、以下の式３０のようになる。
式３０のインデックスｉは、行列の行を示すインデックスである。式３０のインデックスｋは、行列の列を示すインデックスである。この式３０を用いて、行列Ｐの各成分を更新することで、行列Ｐを更新していくこととなる。
また、式２６を用いて、式２４を書き換えると以下の式３１が求まる。
また、式２６を用いて、式２４を書き換えると以下の式３１が求まる。
式３１のインデックスｋは、行列の行を示すインデックスである。式２７のインデックスｊは、行列の列を示すインデックスである。
行列Ａは、非負値であるため、式３１の値は、非負値である必要がある。式３１の右辺の第３項は、必ず非負値となる項である。即ち、式３１の右辺から第３項を除いた式の値が０であるならば、必ず式２７は、非負値となる。そのため、以下の式３２に示すように、式３１の右辺から第３項を除いた式が０であるという条件を満たすように、行列Ａの更新を行うこととする。
行列Ａは、非負値であるため、式３１の値は、非負値である必要がある。式３１の右辺の第３項は、必ず非負値となる項である。即ち、式３１の右辺から第３項を除いた式の値が０であるならば、必ず式２７は、非負値となる。そのため、以下の式３２に示すように、式３１の右辺から第３項を除いた式が０であるという条件を満たすように、行列Ａの更新を行うこととする。
式３２のインデックスｋは、行列の行を示すインデックスである。式３２のインデックスｊは、行列の列を示すインデックスである。
式３２から、式２４における緩和係数ηAの値が、以下の式３３のように求まる。
式３２から、式２４における緩和係数ηAの値が、以下の式３３のように求まる。
式３３のインデックスｋは、行列の行を示すインデックスである。式３３のインデックスｊは、行列の列を示すインデックスである。
式３３を用いて、式２４を書き換えると、以下の式３４のようになる。
式３３を用いて、式２４を書き換えると、以下の式３４のようになる。
式３４のインデックスｋは、行列の行を示すインデックスである。式３４のインデックスｊは、行列の列を示すインデックスである。この式３４を用いて、行列Ａの各成分を更新することで、行列Ａを更新していくこととなる。
また、式２２における第２の制約条件を満たすようにするため、以下の式３５に示すように、行列Ａを更新する。
また、式２２における第２の制約条件を満たすようにするため、以下の式３５に示すように、行列Ａを更新する。
式３５のｋは、行列の行を示すインデックスである。式３５のｊは、行列の行又は列を示すインデックスである。
コスト関数Ｊの値が、収束するまで、式３０を用いて行列Ｐを更新し、式３４、式３５を用いて行列Ａを更新することを繰り返す。例えば、行列Ｐと行列Ａとの更新前後で、コスト関数Ｊの値の変動が予め定められた閾値以下となることが、予め定められた回数連続したら、コスト関数Ｊの値が収束したとする。そして、その際の行列Ａと行列Ｐとをコスト関数Ｊを最小化させる行列とする。
コスト関数Ｊの値が、収束するまで、式３０を用いて行列Ｐを更新し、式３４、式３５を用いて行列Ａを更新することを繰り返す。例えば、行列Ｐと行列Ａとの更新前後で、コスト関数Ｊの値の変動が予め定められた閾値以下となることが、予め定められた回数連続したら、コスト関数Ｊの値が収束したとする。そして、その際の行列Ａと行列Ｐとをコスト関数Ｊを最小化させる行列とする。
コスト関数Ｊが収束した際の行列Ａの各列が示すパターンを確認した。図１０に、コスト関数Ｊが収束した際の行列Ａの各列が示すパターンの一例を示す。行列Ａ列には、列毎に、ｍ（５００）個の値が含まれ、係数αのパターンを示す。
図１０の一番上のグラフは、行列Ａの１列目のｍ個の値（Ａi、1（１＜＝ｉ＜＝ｍ））を、行を示すインデックス順に並べたパターンを示す。図１０の上から２つ目のグラフは、行列Ａの２列目のｍ個の値（Ａi、2（１＜＝ｉ＜＝ｍ））を、行を示すインデックス（ｉ）順に並べたパターンを示す。図１０の上から３つ目のグラフは、行列Ａの３列目のｍ個の値（Ａi、3（１＜＝ｉ＜＝ｍ））を、行を示すインデックス順に並べたパターンを示す。図１０の下から２つ目のグラフは、行列Ａの４列目のｍ個の値（Ａi、4（１＜＝ｉ＜＝ｍ））を、行を示すインデックス順に並べたパターンを示す。図１０の一番下のグラフは、行列Ａの５列目（ｐ列目）のｍ個の値（Ａi、5（１＜＝ｉ＜＝ｍ））を、行を示すインデックス順に並べたパターンを示す。
図１０の一番上のグラフは、行列Ａの１列目のｍ個の値（Ａi、1（１＜＝ｉ＜＝ｍ））を、行を示すインデックス順に並べたパターンを示す。図１０の上から２つ目のグラフは、行列Ａの２列目のｍ個の値（Ａi、2（１＜＝ｉ＜＝ｍ））を、行を示すインデックス（ｉ）順に並べたパターンを示す。図１０の上から３つ目のグラフは、行列Ａの３列目のｍ個の値（Ａi、3（１＜＝ｉ＜＝ｍ））を、行を示すインデックス順に並べたパターンを示す。図１０の下から２つ目のグラフは、行列Ａの４列目のｍ個の値（Ａi、4（１＜＝ｉ＜＝ｍ））を、行を示すインデックス順に並べたパターンを示す。図１０の一番下のグラフは、行列Ａの５列目（ｐ列目）のｍ個の値（Ａi、5（１＜＝ｉ＜＝ｍ））を、行を示すインデックス順に並べたパターンを示す。
また、コスト関数Ｊが収束した際の行列Ｐの各成分の値を確認した。式２０から得られる以下の式３６に示されるように、行列Ｐの各列（第ｊ列Ｐk,j（１＜＝ｋ＜＝ｐ））は、行列Ｙの各列（第ｊ列Ｙi,j（１＜＝ｉ＜＝ｍ））に示されるパターンを行列Ａの各列に示されるパターン（基底）の線形結合で近似するときの係数（重み）を示している。
図１１に、コスト関数Ｊが収束した際の行列Ｐの各行の値を示すグラフを示す。図１１のグラフの横軸は、行列Ｐの列のインデックスを示す。図１１のグラフの縦軸は、行列Ｐの各成分の値を示す。
図１１のライン１１０１は、行列Ｐの１行目のｄ（２５０）個の値を結んだラインである。ライン１１０１は、横軸が０以上、５０未満の範囲（１列目からｑ（５０）列目までの範囲）では、約１２となり、それ以外の範囲では、１２と比べて非常に微小な値（約０）となっている。
図１１のライン１１０２は、行列Ｐの２行目のｄ個の値を結んだラインである。ライン１１０２は、横軸が５０以上、１００未満の範囲（ｑ＋１列目から２ｑ列目までの範囲）では、約１２となり、それ以外の範囲では、１２と比べて非常に微小な値（約０）となっている。
図１１のライン１１０１は、行列Ｐの１行目のｄ（２５０）個の値を結んだラインである。ライン１１０１は、横軸が０以上、５０未満の範囲（１列目からｑ（５０）列目までの範囲）では、約１２となり、それ以外の範囲では、１２と比べて非常に微小な値（約０）となっている。
図１１のライン１１０２は、行列Ｐの２行目のｄ個の値を結んだラインである。ライン１１０２は、横軸が５０以上、１００未満の範囲（ｑ＋１列目から２ｑ列目までの範囲）では、約１２となり、それ以外の範囲では、１２と比べて非常に微小な値（約０）となっている。
図１１のライン１１０３は、行列Ｐの４行目のｄ個の値を結んだラインである。ライン１１０３は、横軸が１００以上、１５０未満の範囲（２ｑ＋１列目から３ｑ列目までの範囲）では、約１２となり、それ以外の範囲では、１２と比べて非常に微小な値（約０）となっている。
図１１のライン１１０４は、行列Ｐの５行目のｄ個の値を結んだラインである。ライン１１０４は、横軸が１５０以上、２００未満の範囲（３ｑ＋１列目から４ｑ列目までの範囲）では、値が約１２となり、それ以外の範囲では、１２と比べて非常に微小な値（約０）となっている。
図１１のライン１１０５は、行列Ｐの３行目のｄ個の値を結んだラインである。ライン１１０５は、横軸が２００以上、２５０未満の範囲（４ｑ＋１列目から５ｑ列目までの範囲）では、約１２となり、それ以外の範囲では、１２と比べて非常に微小な値（約０）となっている。
図１１のライン１１０４は、行列Ｐの５行目のｄ個の値を結んだラインである。ライン１１０４は、横軸が１５０以上、２００未満の範囲（３ｑ＋１列目から４ｑ列目までの範囲）では、値が約１２となり、それ以外の範囲では、１２と比べて非常に微小な値（約０）となっている。
図１１のライン１１０５は、行列Ｐの３行目のｄ個の値を結んだラインである。ライン１１０５は、横軸が２００以上、２５０未満の範囲（４ｑ＋１列目から５ｑ列目までの範囲）では、約１２となり、それ以外の範囲では、１２と比べて非常に微小な値（約０）となっている。
即ち、行列Ｐの各列の値は、１列目からｑ（５０）列目までの範囲では、１行目の値が約１２となり、２〜５行目の値が約０となる。また、行列Ｐの各列の値は、ｑ＋１列目から２ｑ列目までの範囲では、２行目の値が約１２となり、１行目、３〜５行目の値が約０となる。また、行列Ｐの各列の値は、２ｑ＋１列目から３ｑ列目までの範囲では、４行目の値が約１２となり、１〜３行目、５行目の値が約０となる。また、行列Ｐの各列の値は、３ｑ＋１列目から４ｑ列目までの範囲では、５行目の値が約１２となり、１〜４行目の値が約０となる。また、行列Ｐの各列の値は、４ｑ＋１列目から５ｑ列目までの範囲では、３行目の値が約１２となり、１、２、４、５行目の値が約０となる。
そのため、行列ＡＰの１列目からｑ列目までの各列のデータが示すパターンは、図１０の一番上のグラフに示されるパターンに約１２の重み、図１０のそれ以外のグラフに示されるパターンに約０の重み、で図１０の各パターンが重ね合わされたパターンであり、図１０の一番上のグラフに示されるパターンに類似するパターンとなる。また、行列ＡＰのｑ＋１列目から２ｑ列目までの各列のデータが示すパターンは、図１０の上から２つ目のグラフに示されるパターンに約１２の重み、図１０のそれ以外のグラフに示されるパターンに約０の重み、で図１０の各パターンが重ね合わされたパターンであり、図１０の上から２つ目のグラフに示されるパターンに類似するパターンとなる。また、行列ＡＰの２ｑ＋１列目から３ｑ列目までの各列のデータが示すパターンは、図１０の下から２つ目のグラフに示されるパターンに約１２の重み、図１０のそれ以外のグラフに示されるパターンに約０の重み、で図１０の各パターンが重ね合わされたパターンであり、図１０の下から２つ目のグラフに示されるパターンに類似するパターンとなる。また、行列ＡＰの３ｑ＋１列目から４ｑ列目までの各列のデータが示すパターンは、図１０の一番下のグラフに示されるパターンに約１２の重み、図１０のそれ以外のグラフに示されるパターンに約０の重み、で図１０の各パターンが重ね合わされたパターンであり、図１０の一番下のグラフに示されるパターンに類似するパターンとなる。また、行列ＡＰの４ｑ＋１列目から５ｑ列目までの各列のデータが示すパターンは、図１０の上から３番目のグラフに示されるパターンに約１２の重み、図１０のそれ以外のグラフに示されるパターンに約０の重み、で図１０の各パターンが重ね合わされたパターンであり、図１０の上から３つ目のグラフに示されるパターンに類似するパターンとなる。
そのため、行列ＡＰの１列目からｑ列目までの各列のデータが示すパターンは、図１０の一番上のグラフに示されるパターンに約１２の重み、図１０のそれ以外のグラフに示されるパターンに約０の重み、で図１０の各パターンが重ね合わされたパターンであり、図１０の一番上のグラフに示されるパターンに類似するパターンとなる。また、行列ＡＰのｑ＋１列目から２ｑ列目までの各列のデータが示すパターンは、図１０の上から２つ目のグラフに示されるパターンに約１２の重み、図１０のそれ以外のグラフに示されるパターンに約０の重み、で図１０の各パターンが重ね合わされたパターンであり、図１０の上から２つ目のグラフに示されるパターンに類似するパターンとなる。また、行列ＡＰの２ｑ＋１列目から３ｑ列目までの各列のデータが示すパターンは、図１０の下から２つ目のグラフに示されるパターンに約１２の重み、図１０のそれ以外のグラフに示されるパターンに約０の重み、で図１０の各パターンが重ね合わされたパターンであり、図１０の下から２つ目のグラフに示されるパターンに類似するパターンとなる。また、行列ＡＰの３ｑ＋１列目から４ｑ列目までの各列のデータが示すパターンは、図１０の一番下のグラフに示されるパターンに約１２の重み、図１０のそれ以外のグラフに示されるパターンに約０の重み、で図１０の各パターンが重ね合わされたパターンであり、図１０の一番下のグラフに示されるパターンに類似するパターンとなる。また、行列ＡＰの４ｑ＋１列目から５ｑ列目までの各列のデータが示すパターンは、図１０の上から３番目のグラフに示されるパターンに約１２の重み、図１０のそれ以外のグラフに示されるパターンに約０の重み、で図１０の各パターンが重ね合わされたパターンであり、図１０の上から３つ目のグラフに示されるパターンに類似するパターンとなる。
行列Ｙの１列目からｑ列目までの各列のデータが示すパターンは、図９に示す外輪に疵がある軸受１０６に対応するパターンである。また、行列Ｙのｑ＋１列目から２ｑ列目までの各列のデータが示すパターンは、図９に示す転動体に疵がある軸受１０６に対応するパターンである。また、行列Ｙの２ｑ＋１列目から３ｑ列目までの各列のデータが示すパターンは、図９に示す内輪に疵がある軸受１０６に対応するパターンである。また、行列Ｙの３ｑ＋１列目から４ｑ列目までの各列のデータが示すパターンは、図９に示す疵のない正常な軸受１０６に対応するパターンである。また、行列Ｙの４ｑ＋１列目から５ｑ列目までの各列のデータが示すパターンは、図９に示す保持器に疵がある軸受１０６に対応するパターンである。
行列Ｙと行列Ａとを見比べてみると、双方の各列のデータが示すパターンは、類似するパターンとなっていることが分かる。即ち、行列Ａの各列のデータが示すパターンは、それぞれ、各状態の軸受１０６に対応する係数αのパターンとなっていることが分かる。行列Ａの１列目のデータは、外輪に疵がある状態の軸受１０６に対応する係数αの代表的なパターンを示す。行列Ａの２列目のデータは、転動体に疵がある状態の軸受１０６に対応する係数αの代表的なパターンを示す。行列Ａの３列目のデータは、保持器に疵がある状態の軸受１０６に対応する係数αの代表的なパターンを示す。行列Ａの４列目のデータは、内輪に疵のある状態の軸受１０６に対応する係数αの代表的なパターンを示す。行列Ａの５列目のデータは、疵のない正常な状態の軸受１０６に対応する係数αの代表的なパターンを示す。
以上より、行列Ｙを生成し、生成した行列Ｙを、式２２を満たすように、非負値行列因子分解により学習することで、行列Ａの各列（基底）として、各状態の軸受１０６に対応する係数αのパターンが特定されたことが確認された。
行列Ｙと行列Ａとを見比べてみると、双方の各列のデータが示すパターンは、類似するパターンとなっていることが分かる。即ち、行列Ａの各列のデータが示すパターンは、それぞれ、各状態の軸受１０６に対応する係数αのパターンとなっていることが分かる。行列Ａの１列目のデータは、外輪に疵がある状態の軸受１０６に対応する係数αの代表的なパターンを示す。行列Ａの２列目のデータは、転動体に疵がある状態の軸受１０６に対応する係数αの代表的なパターンを示す。行列Ａの３列目のデータは、保持器に疵がある状態の軸受１０６に対応する係数αの代表的なパターンを示す。行列Ａの４列目のデータは、内輪に疵のある状態の軸受１０６に対応する係数αの代表的なパターンを示す。行列Ａの５列目のデータは、疵のない正常な状態の軸受１０６に対応する係数αの代表的なパターンを示す。
以上より、行列Ｙを生成し、生成した行列Ｙを、式２２を満たすように、非負値行列因子分解により学習することで、行列Ａの各列（基底）として、各状態の軸受１０６に対応する係数αのパターンが特定されたことが確認された。
本実施形態の処理は、以下のような処理である。即ち、予め、外輪に疵がある状態、転動体に疵がある状態、内輪に疵がある状態、疵のない正常な状態、保持器に疵がある状態のそれぞれの状態である軸受１０６に対して、軸受１０６の振動の計測データを取得し、取得した計測データそれぞれについて、計測データから求められる自己相関行列Ｒを特異値分解し、得られる固有値の一部を用いて、計測データを近似する修正された自己回帰モデルの係数を求める。そして、各列のデータが求めた係数のそれぞれである行列Ｙを生成し、生成した行列Ｙを、非負値行列因子分解により学習することで、各状態の軸受１０６に対応する修正されタ自己回帰モデルにおける係数のパターンを特定する。
そして、異常診断の対象となる軸受１０６の振動の計測データを取得し、取得した計測データから求められる自己相関行列Ｒを特異値分解し、得られる固有値の一部を用いて、計測データを近似する修正された自己回帰モデルの係数を求める。求めた係数のパターンと、学習により特定した係数のパターンと、に基づいて、軸受１０６の異常を診断する処理である。
そして、異常診断の対象となる軸受１０６の振動の計測データを取得し、取得した計測データから求められる自己相関行列Ｒを特異値分解し、得られる固有値の一部を用いて、計測データを近似する修正された自己回帰モデルの係数を求める。求めた係数のパターンと、学習により特定した係数のパターンと、に基づいて、軸受１０６の異常を診断する処理である。
（システム構成）
図１２は、本実施形態の診断システムのシステム構成の一例を示す図である。診断システムは、周期的な運動を行う物体についての異常診断を行うシステムである。診断システムは、情報処理装置１２００、振動計測装置１２０１を含む。本実施形態では、診断システムは、鉄道台車に利用されている軸受の異常診断を行う。
情報処理装置１２００は、振動計測装置１２０１により計測された信号に基づいて、軸受の異常診断を行うパーソナルコンピュータ（ＰＣ）、サーバ装置、タブレット装置等の情報処理装置である。また、情報処理装置１２００は、電車に組み込まれたコンピュータ等であってもよい。
振動計測装置１２０１は、加速度センサ等のセンサを含み、センサを介して物体の振動を検知し、検知した振動に応じた信号を有線又は無線で情報処理装置１２００等の外部の装置に出力する計測装置である。
図１２は、本実施形態の診断システムのシステム構成の一例を示す図である。診断システムは、周期的な運動を行う物体についての異常診断を行うシステムである。診断システムは、情報処理装置１２００、振動計測装置１２０１を含む。本実施形態では、診断システムは、鉄道台車に利用されている軸受の異常診断を行う。
情報処理装置１２００は、振動計測装置１２０１により計測された信号に基づいて、軸受の異常診断を行うパーソナルコンピュータ（ＰＣ）、サーバ装置、タブレット装置等の情報処理装置である。また、情報処理装置１２００は、電車に組み込まれたコンピュータ等であってもよい。
振動計測装置１２０１は、加速度センサ等のセンサを含み、センサを介して物体の振動を検知し、検知した振動に応じた信号を有線又は無線で情報処理装置１２００等の外部の装置に出力する計測装置である。
（情報処理装置の機能構成）
図１２を参照しながら、情報処理装置１２００が有する機能の一例を説明する。
情報処理装置１２００は、取得部１２１０、決定部１２２０、学習部１２３０、診断部１２４０、出力部１２５０を含む。
≪取得部１２１０≫
取得部１２１０は、周期的な運動を行う物体について、この周期的な運動に係る計測データｙを取得する。
本実施形態では、取得部１２１０は、軸受１０６の内輪を回転させ、位置１１１に設置された振動計測装置１２０１により計測された振動の計測データｙを取得する。
図１２を参照しながら、情報処理装置１２００が有する機能の一例を説明する。
情報処理装置１２００は、取得部１２１０、決定部１２２０、学習部１２３０、診断部１２４０、出力部１２５０を含む。
≪取得部１２１０≫
取得部１２１０は、周期的な運動を行う物体について、この周期的な運動に係る計測データｙを取得する。
本実施形態では、取得部１２１０は、軸受１０６の内輪を回転させ、位置１１１に設置された振動計測装置１２０１により計測された振動の計測データｙを取得する。
≪決定部１２２０≫
決定部１２２０は、取得部１２１０により取得された計測データｙに基づいて、修正された自己回帰モデルにおける係数αを決定する。修正された自己回帰モデルは、計測データｙの実績値ｙ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を取得する。そして、決定部１２２０は、Ｓ１５０１で取得した計測データｙと、予め設定された定数Ｍと、修正された自己回帰モデルにおいて、ある時刻のデータを過去の幾つのデータで近似するかを示す数ｍと、に基づいて、式５と式７とを用いて自己相関行列Ｒを生成する。決定部１２２０は、予め記憶されているＭ、ｍの情報を読み出すことで、Ｍとｍとを取得する。本実施形態では、ｍの値は、５００である。また、Ｍは、ｍよりも大きい整数である。
本実施形態では、診断システムは、図１と同様の状況で位置１１１に設置された振動計測装置１２０１により計測された振動の計測データに基づいて、軸受１０６の異常を診断することとする。
図１４は、学習処理の一例を示すフローチャートである。図１４を用いて、学習パターンを学習により特定する処理を説明する。
本実施形態では、情報処理装置１２００は、軸受１０６が外輪に疵のある状態である場合、軸受１０６が転動体に疵のある状態である場合、軸受１０６が内輪に疵のある状態である場合、軸受１０６に疵がなく正常な状態である場合、軸受１０６が外輪に疵のある状態である場合のそれぞれについて、予め、ｑ（５０）回、以下の処理を繰り返す。
即ち、取得部１２１０は、振動計測装置１２０１により計測された振動の計測データｙを取得する。取得部１２１０は、計測データｙとして、ランダムに設定した時刻を時刻１とした場合の時刻１〜Ｍまでの計測データであるｙ1〜ｙMを取得する。そして、決定部１２２０は、Ｓ１５０１で取得した計測データｙと、予め設定された定数Ｍと、修正された自己回帰モデルにおいて、ある時刻のデータを過去の幾つのデータで近似するかを示す数ｍと、に基づいて、式５と式７とを用いて自己相関行列Ｒを生成する。決定部１２２０は、予め記憶されているＭ、ｍの情報を読み出すことで、Ｍとｍとを取得する。本実施形態では、ｍの値は、５００である。また、Ｍは、ｍよりも大きい整数である。
そして、決定部１２２０は、生成した自己相関行列Ｒを特異値分解することで、式８の直交行列Ｕと対角行列Σを取得し、対角行列Σから自己相関行列Ｒの固有値σ11〜σmmを取得する。決定部１２２０は、自己相関行列Ｒの複数の固有値であるσ11〜σmmのうち、最大のものから使用固有値数ｓだけの固有値σ11〜σssを、修正された自己回帰モデルの係数αを求めるのに利用する自己相関行列Ｒの固有値として選択する。本実施形態では、使用固有値数ｓは、３とするが、２以下としてもよいし、４以上としてもよい。そして、決定部１２２０は、計測データｙと、固有値σ11〜σssと、自己相関行列Ｒの特異値分解により得られた直交行列Ｕと、に基づいて、式１３を用いて、修正された自己回帰モデルの係数αを決定する。そして、決定部１２２０は、決定した係数αを、非負値化して（例えば、係数αのそれぞれの成分に、係数αに含まれる０未満の最低の成分の値に−１をかけた値を足す等する）、学習データとして補助記憶装置１３０２に記憶する。
これにより、情報処理装置１２００は、軸受１０６が外輪に疵のある状態である場合、軸受１０６が転動体に疵のある状態である場合、軸受１０６が内輪に疵のある状態である場合、軸受１０６に疵がなく正常な状態である場合、軸受１０６が外輪に疵のある状態である場合のそれぞれについて、ｑ（５０）個ずつ修正された自己回帰モデルにおける係数αのデータを、学習データとして取得し、記憶することができる。
これにより、情報処理装置１２００は、軸受１０６が外輪に疵のある状態である場合、軸受１０６が転動体に疵のある状態である場合、軸受１０６が内輪に疵のある状態である場合、軸受１０６に疵がなく正常な状態である場合、軸受１０６が外輪に疵のある状態である場合のそれぞれについて、ｑ（５０）個ずつ修正された自己回帰モデルにおける係数αのデータを、学習データとして取得し、記憶することができる。
本実施形態では、軸受１０６が取り得る状態は、疵のない正常な状態、内輪に疵がある状態、保持器に疵がある状態、転動体に疵がある状態、外輪に疵がある状態の５つであるとする。即ち、軸受１０６が取り得る状態の数ｐは、５となる。
Ｓ１４０１において、学習部１２３０は、補助記憶装置１３０２に記憶されている学習データを取得する。
Ｓ１４０２において、学習部１２３０は、Ｓ１４０１で取得した学習データに基づいて、式２３の行列Ｙを生成する。学習部１２３０は、ｍ（５００）×ｄ（２５０）のサイズの行列を生成し、各列のデータを、学習データに含まれる各係数αの値とすることで、行列Ｙを生成する。また、ｍ×ｐのサイズの行列Ａを生成し、行列Ａの各成分に、予め定められた初期値を設定する。また、ｐ×ｄのサイズ行列Ｐを生成し、行列Ｐの各成分に、予め定められた初期値を設定する。
Ｓ１４０３において、学習部１２３０は、Ｓ１４０２で生成した行列Ｙと、行列Ａと、行列Ｐと、予め定められた係数βと、予め定められた係数γと、に基づいて、式２１を用いて、コスト関数Ｊの値を、コストとして求める。学習部１２３０は、求めたコストを、更新前コストとして、主記憶装置１３０１等に記憶する。
Ｓ１４０４において、学習部１２３０は、行列Ｙと、行列Ａと、行列Ｐと、係数γと、サイズがｐ×ｐの全ての成分が１の行列である行列Ｅと、に基づいて、式３４を用いて、行列Ａを更新する。また、学習部１２３０は、更新した行列Ａを、式３５を用いて、更に更新する。また、学習部１２３０は、行列Ｙと、行列Ａと、行列Ｐと、係数βと、サイズがｍ×ｐの全ての成分が１の行列である行列Ｆと、に基づいて、式３０を用いて、行列Ｐを更新する。
Ｓ１４０１において、学習部１２３０は、補助記憶装置１３０２に記憶されている学習データを取得する。
Ｓ１４０２において、学習部１２３０は、Ｓ１４０１で取得した学習データに基づいて、式２３の行列Ｙを生成する。学習部１２３０は、ｍ（５００）×ｄ（２５０）のサイズの行列を生成し、各列のデータを、学習データに含まれる各係数αの値とすることで、行列Ｙを生成する。また、ｍ×ｐのサイズの行列Ａを生成し、行列Ａの各成分に、予め定められた初期値を設定する。また、ｐ×ｄのサイズ行列Ｐを生成し、行列Ｐの各成分に、予め定められた初期値を設定する。
Ｓ１４０３において、学習部１２３０は、Ｓ１４０２で生成した行列Ｙと、行列Ａと、行列Ｐと、予め定められた係数βと、予め定められた係数γと、に基づいて、式２１を用いて、コスト関数Ｊの値を、コストとして求める。学習部１２３０は、求めたコストを、更新前コストとして、主記憶装置１３０１等に記憶する。
Ｓ１４０４において、学習部１２３０は、行列Ｙと、行列Ａと、行列Ｐと、係数γと、サイズがｐ×ｐの全ての成分が１の行列である行列Ｅと、に基づいて、式３４を用いて、行列Ａを更新する。また、学習部１２３０は、更新した行列Ａを、式３５を用いて、更に更新する。また、学習部１２３０は、行列Ｙと、行列Ａと、行列Ｐと、係数βと、サイズがｍ×ｐの全ての成分が１の行列である行列Ｆと、に基づいて、式３０を用いて、行列Ｐを更新する。
Ｓ１４０５において、学習部１２３０は、行列Ｙと、Ｓ１４０３で更新した行列Ａと、Ｓ１４０３で更新した行列Ｐと、係数βと、係数γと、に基づいて、式２１を用いて、コスト関数Ｊの値を、コストとして求める。学習部１２３０は、求めたコストを、更新後コストとして、主記憶装置１３０１等に記憶する。
Ｓ１４０６において、学習部１２３０は、更新後コストと、更新前コストと、の差分の絶対値が、予め定められた閾値以下であるか否かを判定する。学習部１２３０は、更新後コストと、更新前コストと、の差分が、予め定められた閾値以下であると判定した場合、Ｓ１４０７の処理に進む。また、学習部１２３０は、更新後コストと、更新前コストと、の差分の絶対値が、予め定められた閾値よりも大きいと判定した場合、Ｓ１４０８の処理に進む。学習部１２３０は、Ｓ１４０６の処理の後、更新後コストを、新たな更新前コストとして、主記憶装置１３０１に記憶する。
Ｓ１４０７において、学習部１２３０は、主記憶装置１３０１等に記憶されたカウンタの値に１を加える。学習部１２３０は、このカウンタの値を、図１４の処理の開始前に０に初期化している。
Ｓ１４０８において、学習部１２３０は、主記憶装置１３０１等に記憶されたカウンタの値を、０に初期化して、Ｓ１４０４の処理に進む。
Ｓ１４０９において、学習部１２３０は、主記憶装置１３０１等に記憶されたカウンタの値が、予め定められた閾値以上であるか否かを判定する。学習部１２３０は、主記憶装置１３０１等に記憶されたカウンタの値が、予め定められた閾値以上であると判定した場合、コスト関数Ｊの値が極小値に収束したとして、Ｓ１４１０の処理に進む。学習部１２３０は、主記憶装置１３０１等に記憶されたカウンタの値が、予め定められた閾値未満であると判定した場合、Ｓ１４０４の処理に進む。
Ｓ１４０６において、学習部１２３０は、更新後コストと、更新前コストと、の差分の絶対値が、予め定められた閾値以下であるか否かを判定する。学習部１２３０は、更新後コストと、更新前コストと、の差分が、予め定められた閾値以下であると判定した場合、Ｓ１４０７の処理に進む。また、学習部１２３０は、更新後コストと、更新前コストと、の差分の絶対値が、予め定められた閾値よりも大きいと判定した場合、Ｓ１４０８の処理に進む。学習部１２３０は、Ｓ１４０６の処理の後、更新後コストを、新たな更新前コストとして、主記憶装置１３０１に記憶する。
Ｓ１４０７において、学習部１２３０は、主記憶装置１３０１等に記憶されたカウンタの値に１を加える。学習部１２３０は、このカウンタの値を、図１４の処理の開始前に０に初期化している。
Ｓ１４０８において、学習部１２３０は、主記憶装置１３０１等に記憶されたカウンタの値を、０に初期化して、Ｓ１４０４の処理に進む。
Ｓ１４０９において、学習部１２３０は、主記憶装置１３０１等に記憶されたカウンタの値が、予め定められた閾値以上であるか否かを判定する。学習部１２３０は、主記憶装置１３０１等に記憶されたカウンタの値が、予め定められた閾値以上であると判定した場合、コスト関数Ｊの値が極小値に収束したとして、Ｓ１４１０の処理に進む。学習部１２３０は、主記憶装置１３０１等に記憶されたカウンタの値が、予め定められた閾値未満であると判定した場合、Ｓ１４０４の処理に進む。
Ｓ１４１０において、学習部１２３０は、行列Ｐに基づいて、行列Ａの各列のデータが示すパターンが、どの状態である軸受１０６に対応する修正された自己回帰モデルにおける係数αのパターンであるかを特定する。
行列Ｙのある列のデータは、行列Ａと、行列Ｐにおけるその列と、をかけたデータとして求まる。即ち、行列Ｙのある列のデータは、行列Ａの各列のデータを、行列Ｐのその列における各行の値を重みとして、重ね合せたデータとなる。そこで、学習部１２３０は、Ｓ１４１０で以下のようにして、行列Ａの各列のデータが示すパターンが、どの状態である軸受１０６に対応する修正された自己回帰モデルにおける係数αのパターンであるかを特定する。
学習部１２３０は、行列Ｙの各列のうち、ある１つの状態である軸受１０６に対応する係数αを示すデータを格納する列を特定する。例えば、学習部１２３０は、行列Ｙの各列のうち、外輪に疵がある状態の軸受１０６に対応する係数αを示すデータを格納する列である１列目〜ｑ列目を特定する。そして、学習部１２３０は、行列Ｐにおける特定した列（１列目〜ｑ列目）について、行毎に値を集計する。そして、学習部１２３０は、集計した値のうち、最も大きいものに対応する行を特定する。学習部１２３０は、特定した行列Ｐの行とかけ合わされる行列Ａの列を特定する。例えば、学習部１２３０は、行列Ｐの１行目を特定した場合、行列Ｐの１行目とかけ合わされる行列Ａの列は、１列目なので、１列目を特定する。そして、学習部１２３０は、行列Ａの特定した列が示すデータを、その状態である軸受１０６に対応する係数αの代表的なパターンを示すデータとして特定する。
学習部１２３０は、軸受１０６が取り得る各状態について、以上の処理を行うことで、行列Ａの各列のデータが示すパターンが、どの状態である軸受１０６に対応する修正された自己回帰モデルにおける係数αのパターンであるかを特定する。
学習部１２３０は、特定した各状態の軸受１０６に対応する修正された自己回帰モデルにおける係数αの代表的なパターンを示すデータを、学習パターンとして、補助記憶装置１３０２等に記憶する。
行列Ｙのある列のデータは、行列Ａと、行列Ｐにおけるその列と、をかけたデータとして求まる。即ち、行列Ｙのある列のデータは、行列Ａの各列のデータを、行列Ｐのその列における各行の値を重みとして、重ね合せたデータとなる。そこで、学習部１２３０は、Ｓ１４１０で以下のようにして、行列Ａの各列のデータが示すパターンが、どの状態である軸受１０６に対応する修正された自己回帰モデルにおける係数αのパターンであるかを特定する。
学習部１２３０は、行列Ｙの各列のうち、ある１つの状態である軸受１０６に対応する係数αを示すデータを格納する列を特定する。例えば、学習部１２３０は、行列Ｙの各列のうち、外輪に疵がある状態の軸受１０６に対応する係数αを示すデータを格納する列である１列目〜ｑ列目を特定する。そして、学習部１２３０は、行列Ｐにおける特定した列（１列目〜ｑ列目）について、行毎に値を集計する。そして、学習部１２３０は、集計した値のうち、最も大きいものに対応する行を特定する。学習部１２３０は、特定した行列Ｐの行とかけ合わされる行列Ａの列を特定する。例えば、学習部１２３０は、行列Ｐの１行目を特定した場合、行列Ｐの１行目とかけ合わされる行列Ａの列は、１列目なので、１列目を特定する。そして、学習部１２３０は、行列Ａの特定した列が示すデータを、その状態である軸受１０６に対応する係数αの代表的なパターンを示すデータとして特定する。
学習部１２３０は、軸受１０６が取り得る各状態について、以上の処理を行うことで、行列Ａの各列のデータが示すパターンが、どの状態である軸受１０６に対応する修正された自己回帰モデルにおける係数αのパターンであるかを特定する。
学習部１２３０は、特定した各状態の軸受１０６に対応する修正された自己回帰モデルにおける係数αの代表的なパターンを示すデータを、学習パターンとして、補助記憶装置１３０２等に記憶する。
（異常診断処理）
本実施形態では、診断システムは、軸受１０６のどこに疵が存在するかを特定することで、軸受１０６の異常を診断することとするが、軸受１０６に疵が存在するか否かを特定することで、軸受１０６の異常を診断することとしてもよい。診断システムは、予め振動計測装置１２０１により測定された測定データに基づいて、軸受１０６の異常を診断することとしてもよいし、リアルタイムで振動計測装置１２０１により測定されて出力され続けている測定データを用いて、稼働中の軸受１０６の異常を診断することとしてもよい。
図１５は、異常診断処理の一例を示すフローチャートである。
Ｓ１５０１において、取得部１２１０は、振動計測装置１２０１により計測された振動の計測データｙを取得する。取得部１２１０は、計測データｙとして、時刻１〜Ｍまでの計測データであるｙ1〜ｙMを取得する。
Ｓ１５０２において、決定部１２２０は、Ｓ１５０１で取得した計測データｙと、予め設定された定数Ｍと、修正された自己回帰モデルにおいて、ある時刻のデータを過去の幾つのデータで近似するかを示す数ｍと、に基づいて、式５と式７とを用いて自己相関行列Ｒを生成する。決定部１２２０は、予め記憶されているＭ、ｍの情報を読み出すことで、Ｍとｍとを取得する。本実施形態では、ｍの値は、５００である。また、Ｍは、ｍよりも大きい整数である。
本実施形態では、診断システムは、軸受１０６のどこに疵が存在するかを特定することで、軸受１０６の異常を診断することとするが、軸受１０６に疵が存在するか否かを特定することで、軸受１０６の異常を診断することとしてもよい。診断システムは、予め振動計測装置１２０１により測定された測定データに基づいて、軸受１０６の異常を診断することとしてもよいし、リアルタイムで振動計測装置１２０１により測定されて出力され続けている測定データを用いて、稼働中の軸受１０６の異常を診断することとしてもよい。
図１５は、異常診断処理の一例を示すフローチャートである。
Ｓ１５０１において、取得部１２１０は、振動計測装置１２０１により計測された振動の計測データｙを取得する。取得部１２１０は、計測データｙとして、時刻１〜Ｍまでの計測データであるｙ1〜ｙMを取得する。
Ｓ１５０２において、決定部１２２０は、Ｓ１５０１で取得した計測データｙと、予め設定された定数Ｍと、修正された自己回帰モデルにおいて、ある時刻のデータを過去の幾つのデータで近似するかを示す数ｍと、に基づいて、式５と式７とを用いて自己相関行列Ｒを生成する。決定部１２２０は、予め記憶されているＭ、ｍの情報を読み出すことで、Ｍとｍとを取得する。本実施形態では、ｍの値は、５００である。また、Ｍは、ｍよりも大きい整数である。
Ｓ１５０３において、決定部１２２０は、Ｓ１５０２で生成した自己相関行列Ｒを特異値分解することで、式８の直交行列Ｕと対角行列Σを取得し、対角行列Σから自己相関行列Ｒの固有値σ11〜σmmを取得する。
Ｓ１５０４において、決定部１２２０は、自己相関行列Ｒの複数の固有値であるσ11〜σmmのうち、最大のものから使用固有値数ｓだけの固有値σ11〜σssを、修正された自己回帰モデルの係数αを求めるのに利用する自己相関行列Ｒの固有値として選択する。本実施形態では、使用固有値数ｓは、３とするが、２以下としてもよいし、４以上としてもよい。そして、決定部１２２０は、計測データｙと、固有値σ11〜σssと、自己相関行列Ｒの特異値分解により得られた直交行列Ｕと、に基づいて、式１３を用いて、修正された自己回帰モデルの係数αを決定する。Ｓ１５０４で決定された係数αは、計測データから決定された修正された自己回帰モデルにおける係数である修正係数の一例である。Ｓ１５０４で決定された係数αのパターンは、計測パターンの一例である。
Ｓ１５０５において、診断部１２４０は、補助記憶装置１３０２から、Ｓ１４０４で最終的に特定された学習パターンを取得する。
Ｓ１５０４において、決定部１２２０は、自己相関行列Ｒの複数の固有値であるσ11〜σmmのうち、最大のものから使用固有値数ｓだけの固有値σ11〜σssを、修正された自己回帰モデルの係数αを求めるのに利用する自己相関行列Ｒの固有値として選択する。本実施形態では、使用固有値数ｓは、３とするが、２以下としてもよいし、４以上としてもよい。そして、決定部１２２０は、計測データｙと、固有値σ11〜σssと、自己相関行列Ｒの特異値分解により得られた直交行列Ｕと、に基づいて、式１３を用いて、修正された自己回帰モデルの係数αを決定する。Ｓ１５０４で決定された係数αは、計測データから決定された修正された自己回帰モデルにおける係数である修正係数の一例である。Ｓ１５０４で決定された係数αのパターンは、計測パターンの一例である。
Ｓ１５０５において、診断部１２４０は、補助記憶装置１３０２から、Ｓ１４０４で最終的に特定された学習パターンを取得する。
Ｓ１５０６において、診断部１２４０は、ｄを１として、Ｓ１５０４で決定された係数αを用いて、式１９のように、ｍ×ｄのサイズの行列Ｙを生成する。即ち、行列Ｙの１列が、Ｓ１５０４で決定された係数αのパターンを示すこととなる。また、診断部１２４０は、行列Ｙの各成分を、非負値化（例えば、０未満の行列Ｙの成分のうち最低のものに−１をかけた値を、行列Ｙの各成分に足す等する）する。
Ｓ１５０７において、診断部１２４０は、ｍ×ｐのサイズの行列Ａを生成する。そして、診断部１２４０は、行列Ａの各列に、学習パターンのデータを設定する。また、診断部１２４０は、ｐ×ｄ（１）のサイズの行列Ｐを生成し、各成分に予め定められた初期値を設定する。以下では、診断部１２４０は、Ｓ１５０７で生成した行列Ａの値を、確定した値として、行列Ｐを更新することで、行列Ｙを、非負値行列因子分解する。
Ｓ１５０８において、診断部１２４０は、行列Ｙと、行列Ａと、行列Ｐと、予め定められた係数βと、予め定められた係数γと、に基づいて、式２１を用いて、コスト関数Ｊの値を、コストとして求める。診断部１２４０は、求めたコストを、更新前コストとして、主記憶装置１３０１等に記憶する。
Ｓ１５０９において、診断部１２４０は、行列Ｙと、行列Ａと、行列Ｐと、係数βと、サイズがｍ×ｐの全ての成分が１の行列である行列Ｆと、に基づいて、式３０を用いて、行列Ｐを更新する。
Ｓ１５０７において、診断部１２４０は、ｍ×ｐのサイズの行列Ａを生成する。そして、診断部１２４０は、行列Ａの各列に、学習パターンのデータを設定する。また、診断部１２４０は、ｐ×ｄ（１）のサイズの行列Ｐを生成し、各成分に予め定められた初期値を設定する。以下では、診断部１２４０は、Ｓ１５０７で生成した行列Ａの値を、確定した値として、行列Ｐを更新することで、行列Ｙを、非負値行列因子分解する。
Ｓ１５０８において、診断部１２４０は、行列Ｙと、行列Ａと、行列Ｐと、予め定められた係数βと、予め定められた係数γと、に基づいて、式２１を用いて、コスト関数Ｊの値を、コストとして求める。診断部１２４０は、求めたコストを、更新前コストとして、主記憶装置１３０１等に記憶する。
Ｓ１５０９において、診断部１２４０は、行列Ｙと、行列Ａと、行列Ｐと、係数βと、サイズがｍ×ｐの全ての成分が１の行列である行列Ｆと、に基づいて、式３０を用いて、行列Ｐを更新する。
Ｓ１５１０において、診断部１２４０は、行列Ｙと、行列Ａと、Ｓ１５０９で更新した行列Ｐと、係数βと、係数γと、に基づいて、式２１を用いて、コスト関数Ｊの値を、コストとして求める。診断部１２４０は、求めたコストを、更新後コストとして、主記憶装置１３０１等に記憶する。
Ｓ１５１１において、診断部１２４０は、更新後コストと、更新前コストと、の差分の絶対値が、予め定められた閾値以下であるか否かを判定する。診断部１２４０は、更新後コストと、更新前コストと、の差分の絶対値が、予め定められた閾値以下であると判定した場合、Ｓ１５１２の処理に進む。また、診断部１２４０は、更新後コストと、更新前コストと、の差分の絶対値が、予め定められた閾値よりも大きいと判定した場合、Ｓ１５１３の処理に進む。診断部１２４０は、Ｓ１５１１の処理の後、更新後コストを、新たな更新前コストとして、主記憶装置１３０１に記憶する。
Ｓ１５１２において、診断部１２４０は、主記憶装置１３０１等に記憶されたカウンタの値に１を加える。診断部１２４０は、このカウンタの値を、図１５の処理の開始前に０に初期化している。
Ｓ１５１３において、診断部１２４０は、主記憶装置１３０１等に記憶されたカウンタの値を、０に初期化して、Ｓ１５０９の処理に進む。
Ｓ１５１４において、診断部１２４０は、主記憶装置１３０１等に記憶されたカウンタの値が、予め定められた閾値以上であるか否かを判定する。診断部１２４０は、主記憶装置１３０１等に記憶されたカウンタの値が、予め定められた閾値以上であると判定した場合、コスト関数Ｊの値が極小値に収束したとして、Ｓ１５１５の処理に進む。診断部１２４０は、主記憶装置１３０１等に記憶されたカウンタの値が、予め定められた閾値未満であると判定した場合、Ｓ１５０９の処理に進む。
Ｓ１５１１において、診断部１２４０は、更新後コストと、更新前コストと、の差分の絶対値が、予め定められた閾値以下であるか否かを判定する。診断部１２４０は、更新後コストと、更新前コストと、の差分の絶対値が、予め定められた閾値以下であると判定した場合、Ｓ１５１２の処理に進む。また、診断部１２４０は、更新後コストと、更新前コストと、の差分の絶対値が、予め定められた閾値よりも大きいと判定した場合、Ｓ１５１３の処理に進む。診断部１２４０は、Ｓ１５１１の処理の後、更新後コストを、新たな更新前コストとして、主記憶装置１３０１に記憶する。
Ｓ１５１２において、診断部１２４０は、主記憶装置１３０１等に記憶されたカウンタの値に１を加える。診断部１２４０は、このカウンタの値を、図１５の処理の開始前に０に初期化している。
Ｓ１５１３において、診断部１２４０は、主記憶装置１３０１等に記憶されたカウンタの値を、０に初期化して、Ｓ１５０９の処理に進む。
Ｓ１５１４において、診断部１２４０は、主記憶装置１３０１等に記憶されたカウンタの値が、予め定められた閾値以上であるか否かを判定する。診断部１２４０は、主記憶装置１３０１等に記憶されたカウンタの値が、予め定められた閾値以上であると判定した場合、コスト関数Ｊの値が極小値に収束したとして、Ｓ１５１５の処理に進む。診断部１２４０は、主記憶装置１３０１等に記憶されたカウンタの値が、予め定められた閾値未満であると判定した場合、Ｓ１５０９の処理に進む。
Ｓ１５１５において、診断部１２４０は、行列Ｐに基づいて、軸受１０６の状態を診断する。診断部１２４０は、例えば、行列Ｐの各行の値のうち、最大のものに対応する行を特定する。診断部１２４０は、特定した行列Ｐの行とかけあわされる行列Ａの列を特定する。そして、診断部１２４０は、軸受１０６の状態を、特定した行列Ａの列が示す学習パターンに対応する状態となっていると診断する。
また、診断部１２４０は、行列Ｐの各行の値のうち、予め定められた閾値以上のものに対応する行を特定し、特定した行列Ｐの行とかけあわされる行列Ａの列を特定し、軸受１０６の状態を、特定した行列Ａの列が示す学習パターンに対応する状態となっていると診断してもよい。また、診断部１２４０は、行列Ｐの各行の値のうち、予め定められた閾値以上のものに対応する行として、複数の行を特定した場合、特定した複数の行とかけあわされる行列Ａの複数の列を特定し、軸受１０６の状態を、特定した行列Ａの複数の列それぞれが示す学習パターンに対応する複数の状態の複合状態となっていると診断してもよい。また、診断部１２４０は、例えば、行列Ｐの各行の値に、予め定められた閾値以上となる値がない場合、軸受１０６に未知の状態（例えば、未知の異常が発生している状態）となっていると診断してもよい。
Ｓ１５１６において、出力部１２５０は、Ｓ１５１５での軸受１０６の異常の診断の結果を示す情報を出力する。出力部１２５０は、例えば、表示装置に、軸受１０６の異常の診断の結果を示す情報を表示することで出力する。また、出力部１２５０は、例えば、軸受１０６の異常の診断の結果を示す情報を、補助記憶装置１３０２に記憶することで出力することとしてもよい。また、出力部１２５０は、例えば、外部のサーバ装置等の設定された送信先に、軸受１０６の異常の診断の結果を示す情報を送信することで出力することとしてもよい。
また、診断部１２４０は、行列Ｐの各行の値のうち、予め定められた閾値以上のものに対応する行を特定し、特定した行列Ｐの行とかけあわされる行列Ａの列を特定し、軸受１０６の状態を、特定した行列Ａの列が示す学習パターンに対応する状態となっていると診断してもよい。また、診断部１２４０は、行列Ｐの各行の値のうち、予め定められた閾値以上のものに対応する行として、複数の行を特定した場合、特定した複数の行とかけあわされる行列Ａの複数の列を特定し、軸受１０６の状態を、特定した行列Ａの複数の列それぞれが示す学習パターンに対応する複数の状態の複合状態となっていると診断してもよい。また、診断部１２４０は、例えば、行列Ｐの各行の値に、予め定められた閾値以上となる値がない場合、軸受１０６に未知の状態（例えば、未知の異常が発生している状態）となっていると診断してもよい。
Ｓ１５１６において、出力部１２５０は、Ｓ１５１５での軸受１０６の異常の診断の結果を示す情報を出力する。出力部１２５０は、例えば、表示装置に、軸受１０６の異常の診断の結果を示す情報を表示することで出力する。また、出力部１２５０は、例えば、軸受１０６の異常の診断の結果を示す情報を、補助記憶装置１３０２に記憶することで出力することとしてもよい。また、出力部１２５０は、例えば、外部のサーバ装置等の設定された送信先に、軸受１０６の異常の診断の結果を示す情報を送信することで出力することとしてもよい。
（まとめ）
以上、本実施形態では、診断システムは、予め、複数の状態の軸受１０６の周期的な運動から計測された計測データから、自己相関行列Ｒを生成し、自己相関行列Ｒを特異値分解し、得られた固有値のうち、最大のものから設定された数を用いて、計測データを近似する修正された自己回帰モデルの係数αを、学習データとして決定した。そして、診断システムは、学習データを用いて、各状態の軸受１０６に対応する代表的な係数αのパターンを、学習パターンとして学習により特定した。そして、診断システムは、診断対象となる軸受１０６の周期的な運動から計測された計測データから修正された自己回帰モデルの係数αを決定し、決定した係数αのパターンと、学習により特定した学習パターンと、に基づいて、軸受１０６の異常を診断した。診断システムは、自己相関行列Ｒの固有値のうち一部の固有値を用いることで、軸受１０６の異常診断に有用な成分が残り、有用でない成分が残らないように、計測データを近似する修正された自己回帰モデルの係数αを決定できる。この係数αのパターンは、軸受１０６の状態毎に、それぞれ特徴的なパターンを示すことが実験により見出され、本実施形態では、診断システムは、この係数αのパターンに基づいて、軸受１０６の異常を診断した。これにより、診断システムは、計測データの軸受１０６の異常診断に有用な成分に基づいて、より精度よく異常診断を行うことができる。
以上、本実施形態では、診断システムは、予め、複数の状態の軸受１０６の周期的な運動から計測された計測データから、自己相関行列Ｒを生成し、自己相関行列Ｒを特異値分解し、得られた固有値のうち、最大のものから設定された数を用いて、計測データを近似する修正された自己回帰モデルの係数αを、学習データとして決定した。そして、診断システムは、学習データを用いて、各状態の軸受１０６に対応する代表的な係数αのパターンを、学習パターンとして学習により特定した。そして、診断システムは、診断対象となる軸受１０６の周期的な運動から計測された計測データから修正された自己回帰モデルの係数αを決定し、決定した係数αのパターンと、学習により特定した学習パターンと、に基づいて、軸受１０６の異常を診断した。診断システムは、自己相関行列Ｒの固有値のうち一部の固有値を用いることで、軸受１０６の異常診断に有用な成分が残り、有用でない成分が残らないように、計測データを近似する修正された自己回帰モデルの係数αを決定できる。この係数αのパターンは、軸受１０６の状態毎に、それぞれ特徴的なパターンを示すことが実験により見出され、本実施形態では、診断システムは、この係数αのパターンに基づいて、軸受１０６の異常を診断した。これにより、診断システムは、計測データの軸受１０６の異常診断に有用な成分に基づいて、より精度よく異常診断を行うことができる。
＜変形例＞
本実施形態では、診断システムは、鉄道台車の軸受１０６の異常診断を行うこととした。しかし、診断システムは、歯車等の回転体、振動子等の振動体、バネ等の伸縮体、生物の心臓等の周期的な運動を行う他の物体についての異常診断を行うこととしてもよい。
本実施形態では、ｍは、予め設定された５００という数であるとした。しかし、診断対象に応じて実験を行う等して適切に計測データを近似できる過去のデータの数を特定し、特定した数をｍの値としてもよい。例えば、情報処理装置１２００は、診断対象から計測された計測データについて、ある時刻の計測データを、その時刻よりも過去の計測データを複数用いて自己回帰モデル等を利用して近似して、近似した値とその時刻の計測データとの差分が設定された閾値未満となる場合の近似に用いられた過去の計測データの数を、ｍとして決定してもよい。
本実施形態では、情報処理装置１２００は、使用固有値数ｓを３とした。しかし、情報処理装置１２００は、以下のようにして、使用固有値数ｓの値を決定してもよい。即ち、情報処理装置１２００は、まず、自己相関行列Ｒの固有値全ての合計値を求める。そして、情報処理装置１２００は、自己相関行列Ｒの固有値のうち、最大のものからａ個の合計が、求めた合計値の設定された割合（例えば９０％）以上となるａのうち、最小のものを使用固有値数ｓとして決定してもよい。
本実施形態では、情報処理装置１２００は、自己相関行列Ｒの固有値のうち、最大のものから、使用固有値数ｓ個の固有値を用いて、修正された自己回帰モデルにおける係数αを決定することとした。しかし、情報処理装置１２００は、自己相関行列Ｒの固有値のうち、任意のｓ個の固有値を用いて、修正された自己回帰モデルにおける係数αを決定することとしてもよい。例えば、情報処理装置１２００は自己相関行列Ｒの固有値のうち、最大のものと、３つ目に大きいものと、の２つを利用して、修正された自己回帰モデルにおける係数αを決定することとしてもよい。
本実施形態では、診断システムは、鉄道台車の軸受１０６の異常診断を行うこととした。しかし、診断システムは、歯車等の回転体、振動子等の振動体、バネ等の伸縮体、生物の心臓等の周期的な運動を行う他の物体についての異常診断を行うこととしてもよい。
本実施形態では、ｍは、予め設定された５００という数であるとした。しかし、診断対象に応じて実験を行う等して適切に計測データを近似できる過去のデータの数を特定し、特定した数をｍの値としてもよい。例えば、情報処理装置１２００は、診断対象から計測された計測データについて、ある時刻の計測データを、その時刻よりも過去の計測データを複数用いて自己回帰モデル等を利用して近似して、近似した値とその時刻の計測データとの差分が設定された閾値未満となる場合の近似に用いられた過去の計測データの数を、ｍとして決定してもよい。
本実施形態では、情報処理装置１２００は、使用固有値数ｓを３とした。しかし、情報処理装置１２００は、以下のようにして、使用固有値数ｓの値を決定してもよい。即ち、情報処理装置１２００は、まず、自己相関行列Ｒの固有値全ての合計値を求める。そして、情報処理装置１２００は、自己相関行列Ｒの固有値のうち、最大のものからａ個の合計が、求めた合計値の設定された割合（例えば９０％）以上となるａのうち、最小のものを使用固有値数ｓとして決定してもよい。
本実施形態では、情報処理装置１２００は、自己相関行列Ｒの固有値のうち、最大のものから、使用固有値数ｓ個の固有値を用いて、修正された自己回帰モデルにおける係数αを決定することとした。しかし、情報処理装置１２００は、自己相関行列Ｒの固有値のうち、任意のｓ個の固有値を用いて、修正された自己回帰モデルにおける係数αを決定することとしてもよい。例えば、情報処理装置１２００は自己相関行列Ｒの固有値のうち、最大のものと、３つ目に大きいものと、の２つを利用して、修正された自己回帰モデルにおける係数αを決定することとしてもよい。
本実施形態では、学習データから非負値行列因子分解により得られる学習パターンと各状態との直接的な対応関係が存在した。即ち、各状態の軸受１０６に対応する係数αのパターンを学習パターンの線形結合で近似するとき、各状態に対応する学習パターンに重みが集中した。しかしながら、このような学習パターンにかかる重みの集中が見られなくても、各学習パターンにかかる重みの組み合わせと、各状態の軸受１０６に対応する係数αのパターンとの対応関係が存在すれば、この重みを用いた学習パターンの線形結合を各状態に対応する学習パターンとして新たに取り直せばよい。
本実施形態では、診断システムは、軸受１０６が５つの状態それぞれである場合について、修正された自己回帰モデルにおける係数αを代表するパターンを５つ学習により特定することとした。しかし、診断システムは、軸受１０６が取り得る状態が６つ以上ある場合、診断システムは、軸受１０６の状態それぞれについて、修正された自己回帰モデルにおける係数αを代表するパターンを６つ以上学習により特定することとしてもよい。
この場合、学習部１２３０は、例えば、ｐの値を６以上の値として、式１９を用いて行列Ｙを生成し、コスト関数Ｊを最小化するように、ｍ×ｐのサイズの行列Ａと、ｐ×ｄのサイズの行列Ｐと、を更新することで、行列Ａの各列として、各状態の軸受１０６に対応する係数αのパターンを学習により特定する。
軸受１０６が、例えば、外輪に擦り疵がある状態、外輪に凹凸がある状態、内輪に擦り疵がある状態、内輪に凹凸がある状態、転動体に擦り疵がある状態、転動体に凹凸がある状態、保持器に擦り疵がある状態、保持器に凹凸がある状態、疵のない正常な状態の９個の状態を取り得る場合、学習部１２３０は、以下のような処理を行うこととしてもよい。即ち、学習部１２３０は、各状態の軸受１０６から取得された計測データから決定された修正された自己回帰モデルにおける係数αを学習データとして、各状態の軸受１０６に対応する９個の係数αのパターンを学習により特定してもよい。
本実施形態では、診断システムは、軸受１０６が５つの状態それぞれである場合について、修正された自己回帰モデルにおける係数αを代表するパターンを５つ学習により特定することとした。しかし、診断システムは、軸受１０６が取り得る状態が６つ以上ある場合、診断システムは、軸受１０６の状態それぞれについて、修正された自己回帰モデルにおける係数αを代表するパターンを６つ以上学習により特定することとしてもよい。
この場合、学習部１２３０は、例えば、ｐの値を６以上の値として、式１９を用いて行列Ｙを生成し、コスト関数Ｊを最小化するように、ｍ×ｐのサイズの行列Ａと、ｐ×ｄのサイズの行列Ｐと、を更新することで、行列Ａの各列として、各状態の軸受１０６に対応する係数αのパターンを学習により特定する。
軸受１０６が、例えば、外輪に擦り疵がある状態、外輪に凹凸がある状態、内輪に擦り疵がある状態、内輪に凹凸がある状態、転動体に擦り疵がある状態、転動体に凹凸がある状態、保持器に擦り疵がある状態、保持器に凹凸がある状態、疵のない正常な状態の９個の状態を取り得る場合、学習部１２３０は、以下のような処理を行うこととしてもよい。即ち、学習部１２３０は、各状態の軸受１０６から取得された計測データから決定された修正された自己回帰モデルにおける係数αを学習データとして、各状態の軸受１０６に対応する９個の係数αのパターンを学習により特定してもよい。
また、診断システムは、軸受１０６が取り得る状態が４つ以下しかない場合、診断システムは、軸受１０６の状態それぞれについて、修正された自己回帰モデルにおける係数αを代表するパターンを４つ以下学習により特定することとしてもよい。
この場合、学習部１２３０は、例えば、ｐの値を４以下の値として、式１９を用いて行列Ｙを生成し、コスト関数Ｊを最小化するように、ｍ×ｐのサイズの行列Ａと、ｐ×ｄのサイズの行列Ｐと、を更新することで、行列Ａの各列として、各状態の軸受１０６に対応する係数αのパターンを学習により特定する。
軸受１０６が、例えば、疵のない正常な状態、何れかの部位に疵のある状態の２個の状態を取り得る場合、学習部１２３０は、それぞれの軸受１０６から取得された計測データから決定された修正された自己回帰モデルにおける係数αを学習データとして、各状態の軸受１０６に対応する２個の係数αのパターンを学習により特定してもよい。
この場合、学習部１２３０は、例えば、ｐの値を４以下の値として、式１９を用いて行列Ｙを生成し、コスト関数Ｊを最小化するように、ｍ×ｐのサイズの行列Ａと、ｐ×ｄのサイズの行列Ｐと、を更新することで、行列Ａの各列として、各状態の軸受１０６に対応する係数αのパターンを学習により特定する。
軸受１０６が、例えば、疵のない正常な状態、何れかの部位に疵のある状態の２個の状態を取り得る場合、学習部１２３０は、それぞれの軸受１０６から取得された計測データから決定された修正された自己回帰モデルにおける係数αを学習データとして、各状態の軸受１０６に対応する２個の係数αのパターンを学習により特定してもよい。
本実施形態では、診断システムは、軸受１０６が予め定められた５つの状態（疵のない正常な状態、外輪に疵がある状態、内輪に疵がある状態、転動体に疵がある状態、保持器に疵がある状態）のうちのどれかを特定することで、異常を診断した。しかし、軸受１０６がどのような状態を取り得るかを正確に把握できないが、疵のない正常な状態か否かを診断したいというような場合がある。このような場合、診断システムは、学習パターンのうち、正常な状態の軸受１０６に対応するもののみを用いて、以下のような処理を行うことで、軸受１０６が正常であるか否かを診断してもよい。
即ち、診断システムは、図１５の処理において、Ｓ１５０９で、行列Ｐを更新する際に、行列Ａの成分のうち、疵のない正常な状態の軸受１０６に対応する学習パターンを格納する列以外の成分を、式３４を用いて更新する。そして、診断システムは、Ｓ１５１５で、行列Ａの成分のうち、疵のない正常な状態の軸受１０６に対応する学習パターンを格納する列と、かけあわされる行列Ｐの行の値が、行列Ｐの行の中で最大であるなら、軸受１０６が正常な状態であると診断することとしてもよい。
このように、診断システムは、学習により特定した学習パターンのうちの一部を用いて、軸受１０６の異常を診断してもよい。
即ち、診断システムは、図１５の処理において、Ｓ１５０９で、行列Ｐを更新する際に、行列Ａの成分のうち、疵のない正常な状態の軸受１０６に対応する学習パターンを格納する列以外の成分を、式３４を用いて更新する。そして、診断システムは、Ｓ１５１５で、行列Ａの成分のうち、疵のない正常な状態の軸受１０６に対応する学習パターンを格納する列と、かけあわされる行列Ｐの行の値が、行列Ｐの行の中で最大であるなら、軸受１０６が正常な状態であると診断することとしてもよい。
このように、診断システムは、学習により特定した学習パターンのうちの一部を用いて、軸受１０６の異常を診断してもよい。
本実施形態では、式１９に示す行列Ｙの各列には、複数の状態の軸受１０６に対応する計測データから決定された修正された自己回帰モデルにおける係数αが、状態毎に、ｑ個ずつ格納されることとした。しかし、行列Ｙの各列には、複数の状態の軸受１０６に対応する計測データから決定された修正された自己回帰モデルにおける係数αが、状態毎に、異なる数ずつ格納されることとしてもよい。
本実施形態では、診断システムは、行列Ｙを、同じ状態の軸受１０６に対応する係数αを格納する列が連続するように、生成した。しかし、診断システムは、同じ状態の軸受１０６に対応する係数αを格納する列が連続しないように、行列Ｙを生成してもよい。
本実施形態では、診断システムは、各列に各状態の軸受１０６に対応する係数αが格納された行列Ｙを生成し、行列Ｙを非負値行列因子分解により学習することで、各状態の軸受１０６に対応する係数αの代表的なパターンを特定した。しかし、診断システムは、例えば、学習データを取得する際の処理と同様の処理で、ある状態の軸受１０６に対応する係数αを複数決定し、決定した複数の係数αの平均を求めることにより学習することで、求めた平均の係数が示すパターンを、その状態の軸受１０６に対応する係数αの代表的なパターンとして特定してもよい。
本実施形態では、診断システムは、行列Ｙを、同じ状態の軸受１０６に対応する係数αを格納する列が連続するように、生成した。しかし、診断システムは、同じ状態の軸受１０６に対応する係数αを格納する列が連続しないように、行列Ｙを生成してもよい。
本実施形態では、診断システムは、各列に各状態の軸受１０６に対応する係数αが格納された行列Ｙを生成し、行列Ｙを非負値行列因子分解により学習することで、各状態の軸受１０６に対応する係数αの代表的なパターンを特定した。しかし、診断システムは、例えば、学習データを取得する際の処理と同様の処理で、ある状態の軸受１０６に対応する係数αを複数決定し、決定した複数の係数αの平均を求めることにより学習することで、求めた平均の係数が示すパターンを、その状態の軸受１０６に対応する係数αの代表的なパターンとして特定してもよい。
本実施形態では、診断システムは、軸受１０６の測定データから決定された修正された自己回帰モデルにおける係数αのパターンに基づいて、軸受１０６の異常を診断した。しかし、診断システムは、Ｓ１５１５で、軸受１０６の状態を特定できなかった場合、例えば、式１８を用いて、修正された自己回帰モデルの周波数特性を求め、求めた周波数特性に基づいて、軸受１０６の異常を診断してもよい。また、診断システムは、例えば、式１８を用いて、修正された自己回帰モデルの周波数特性を求め、求めた周波数特性に基づいて、軸受１０６の異常を診断し、診断に失敗した場合、図１５の処理を行うこととしてもよい。
このように、診断システムは、軸受１０６の測定データから決定された修正された自己回帰モデルにおける係数αのパターンと、軸受１０６の測定データから決定された修正された自己回帰モデルの周波数特性と、に基づいて、異常を診断することで、より精度よく軸受１０６の異常を診断できる。
このように、診断システムは、軸受１０６の測定データから決定された修正された自己回帰モデルにおける係数αのパターンと、軸受１０６の測定データから決定された修正された自己回帰モデルの周波数特性と、に基づいて、異常を診断することで、より精度よく軸受１０６の異常を診断できる。
本実施形態では、情報処理装置１２００は、補助記憶装置１３０２に記憶されたプログラムを実行することで、図１４、１５の処理を実現することとした。しかし、情報処理装置１２００は、外付けの記憶媒体や外部の記憶サーバ等に記憶されたプログラムを実行することで、図１４、１５の処理を実現することとしてもよい。
また、例えば、上述した情報処理装置１２００の機能の一部又は全てをハードウェアとして情報処理装置１２００に実装してもよい。
また、以上説明した本発明の実施形態は、何れも本発明を実施するにあたっての具体化の例を示したものに過ぎず、これらによって本発明の技術的範囲が限定的に解釈されてはならないものである。即ち、本発明はその技術思想、またはその主要な特徴から逸脱することなく、様々な形で実施することができる。
また、例えば、上述した情報処理装置１２００の機能の一部又は全てをハードウェアとして情報処理装置１２００に実装してもよい。
また、以上説明した本発明の実施形態は、何れも本発明を実施するにあたっての具体化の例を示したものに過ぎず、これらによって本発明の技術的範囲が限定的に解釈されてはならないものである。即ち、本発明はその技術思想、またはその主要な特徴から逸脱することなく、様々な形で実施することができる。
１２００ 情報処理装置
１２１０ 取得部
１２２０ 決定部
１２３０ 学習部
１２４０ 診断部
１２５０ 出力部
１２０１ 振動計測装置
１３００ ＣＰＵ
１２１０ 取得部
１２２０ 決定部
１２３０ 学習部
１２４０ 診断部
１２５０ 出力部
１２０１ 振動計測装置
１３００ ＣＰＵ
Claims (11)
- 周期的な運動を行う物体の前記周期的な運動に係る計測データを取得する取得手段と、
前記取得手段により取得された前記計測データに基づいて、修正された自己回帰モデルにおける係数である修正係数を決定する決定手段と、
前記決定手段により決定された前記修正係数のパターンである計測パターンに基づいて、前記物体の異常を診断する診断手段と、
を有し、
前記修正された自己回帰モデルは、前記計測データの実績値と、前記実績値に対する前記修正係数と、を用いて、前記計測データの予測値を表す式であり、
前記決定手段は、前記計測データから導出される自己相関行列を特異値分解することで導出される前記自己相関行列の固有値を対角成分とする対角行列と、前記自己相関行列の固有ベクトルを列成分とする直交行列と、から導出される第１の行列を係数行列とし、前記計測データから導出される自己相関ベクトルを定数ベクトルとする方程式を用いて、前記修正係数を決定し、
前記自己相関ベクトルは、時差が１から前記修正された自己回帰モデルで用いられる前記実績値の数であるｍまでの前記計測データの自己相関を成分とするベクトルであり、
前記自己相関行列は、時差が０からｍ−１までの前記計測データの自己相関を成分とする行列であり、
前記第１の行列は、１以上且つｍ未満の設定された数であるｓに対して、前記自己相関行列のｓ個の固有値と前記対角行列とから導出される第２の行列Σsと、前記ｓ個の固有値と前記直交行列とから導出される第３の行列Ｕsと、から導出される行列ＵsΣsＵs Tであり、
前記第２の行列は、前記対角行列の部分行列であって、前記ｓ個の固有値を対角成分とする行列であり、
前記第３の行列は、前記直交行列の部分行列であって、前記ｓ個の固有値に対応する固有ベクトルを列成分ベクトルとする行列である情報処理装置。 - 前記ｓ個の固有値は、前記自己相関行列の固有値のうち、値が最大の固有値を含む請求項１記載の情報処理装置。
- 前記ｓ個の固有値は、前記自己相関行列の固有値の中から、値が大きいものから順に選ばれた固有値である請求項２記載の情報処理装置。
- 前記自己相関行列は、前記計測データの予測値と当該計測データの予測値に対応する時刻における前記計測データの実測値との二乗誤差を最小化する条件を示す条件式に含まれる行列である請求項１乃至３何れか１項記載の情報処理装置。
- 予め定められた複数の状態それぞれにある前記物体から決定された一つ以上の前記修正された自己回帰モデルにおける係数である学習データに基づいて、前記複数の状態それぞれに対応する前記修正された自己回帰モデルにおける係数のパターンを、学習パターンとして、学習により特定する学習手段を更に有し、
前記診断手段は、前記決定手段により決定された前記計測パターンと、前記学習手段により学習により特定された前記学習パターンと、に基づいて、前記物体の異常を診断する請求項１乃至４何れか１項記載の情報処理装置。 - 前記学習手段は、各列に前記学習データが格納された行列を、前記複数の状態の数を基底の数として、非負値行列因子分解することで、前記学習パターンを学習により特定する請求項５記載の情報処理装置。
- 前記計測データは、前記物体の振動に係る計測データである請求項１乃至６何れか１項記載の情報処理装置。
- 前記診断手段による診断の結果に係る情報を出力する出力手段を更に有する請求項１乃至７何れか１項記載の情報処理装置。
- 前記物体は、鉄道台車に利用される軸受であり、
前記診断手段は、前記決定手段により決定された前記計測パターンに基づいて、前記物体における疵の存在を含む異常を診断する請求項１乃至８何れか１項記載の情報処理装置。 - 情報処理装置が実行する情報処理方法であって、
周期的な運動を行う物体の前記周期的な運動に係る計測データを取得する取得ステップと、
前記取得ステップで取得された前記計測データに基づいて、修正された自己回帰モデルにおける係数である修正係数を決定する決定ステップと、
前記決定ステップで決定された前記修正係数のパターンである計測パターンに基づいて、前記物体の異常を診断する診断ステップと、
を含み、
前記修正された自己回帰モデルは、前記計測データの実績値と、前記実績値に対する前記修正係数と、を用いて、前記計測データの予測値を表す式であり、
前記決定ステップでは、前記計測データから導出される自己相関行列を特異値分解することで導出される前記自己相関行列の固有値を対角成分とする対角行列と、前記自己相関行列の固有ベクトルを列成分とする直交行列と、から導出される第１の行列を係数行列とし、前記計測データから導出される自己相関ベクトルを定数ベクトルとする方程式を用いて、前記修正係数を決定し、
前記自己相関ベクトルは、時差が１から前記修正された自己回帰モデルで用いられる前記実績値の数であるｍまでの前記計測データの自己相関を成分とするベクトルであり、
前記自己相関行列は、時差が０からｍ−１までの前記計測データの自己相関を成分とする行列であり、
前記第１の行列は、１以上且つｍ未満の設定された数であるｓに対して、前記自己相関行列のｓ個の固有値と前記対角行列とから導出される第２の行列Σsと、前記ｓ個の固有値と前記直交行列とから導出される第３の行列Ｕsと、から導出される行列ＵsΣsＵs Tであり、
前記第２の行列は、前記対角行列の部分行列であって、前記ｓ個の固有値を対角成分とする行列であり、
前記第３の行列は、前記直交行列の部分行列であって、前記ｓ個の固有値に対応する固有ベクトルを列成分ベクトルとする行列である情報処理方法。 - コンピュータを、請求項１乃至９何れか１項記載の情報処理装置の各手段として機能させるためのプログラム。
Priority Applications (4)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2017206597A JP6848813B2 (ja) | 2017-10-25 | 2017-10-25 | 情報処理装置、情報処理方法及びプログラム |
| PCT/JP2018/028942 WO2019026980A1 (ja) | 2017-08-04 | 2018-08-01 | 情報処理装置、情報処理方法及びプログラム |
| EP18840775.3A EP3663741A4 (en) | 2017-08-04 | 2018-08-01 | INFORMATION PROCESSING DEVICE, INFORMATION PROCESSING METHOD AND PROGRAM |
| CN201880052067.6A CN110998272A (zh) | 2017-08-04 | 2018-08-01 | 信息处理装置、信息处理方法及程序 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2017206597A JP6848813B2 (ja) | 2017-10-25 | 2017-10-25 | 情報処理装置、情報処理方法及びプログラム |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JP2019079352A true JP2019079352A (ja) | 2019-05-23 |
| JP6848813B2 JP6848813B2 (ja) | 2021-03-24 |
Family
ID=65232868
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2017206597A Active JP6848813B2 (ja) | 2017-08-04 | 2017-10-25 | 情報処理装置、情報処理方法及びプログラム |
Country Status (4)
| Country | Link |
|---|---|
| EP (1) | EP3663741A4 (ja) |
| JP (1) | JP6848813B2 (ja) |
| CN (1) | CN110998272A (ja) |
| WO (1) | WO2019026980A1 (ja) |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2021135683A (ja) * | 2020-02-26 | 2021-09-13 | 株式会社東芝 | 学習装置、推論装置、学習方法及び推論方法 |
| CN113899550A (zh) * | 2021-09-28 | 2022-01-07 | 新黎明科技股份有限公司 | 迭代高阶能量算子融合防爆电机轴承故障诊断方法及系统 |
| WO2022044609A1 (ja) * | 2020-08-24 | 2022-03-03 | 古野電気株式会社 | 船舶航行支援装置、船舶航行支援方法、および、船舶航行支援プログラム |
| WO2022044608A1 (ja) * | 2020-08-24 | 2022-03-03 | 古野電気株式会社 | 船舶航行支援装置、船舶航行支援方法、および、船舶航行支援プログラム |
Families Citing this family (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2021117752A1 (ja) * | 2019-12-11 | 2021-06-17 | Ntn株式会社 | 転がり軸受の状態監視方法及び転がり軸受の状態監視装置 |
| JPWO2022059719A1 (ja) * | 2020-09-15 | 2022-03-24 | ||
| WO2022059720A1 (ja) * | 2020-09-15 | 2022-03-24 | 国立大学法人京都大学 | 構造物診断システム、構造物診断方法、および構造物診断プログラム |
| CN116558828B (zh) * | 2023-07-10 | 2023-09-15 | 昆明理工大学 | 基于自相关系数稀疏度特征的滚动轴承健康状态评估方法 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JPH0124246B2 (ja) * | 1981-01-30 | 1989-05-10 | Sumitomo Metal Ind | |
| JP2003050158A (ja) * | 2001-08-06 | 2003-02-21 | Mitsui Eng & Shipbuild Co Ltd | 周期的運動体の異常診断方法および装置 |
| WO2004017038A1 (en) * | 2002-08-14 | 2004-02-26 | Katsumi Hiramatsu | Detector of defects for rotating machinery |
| JP2012159298A (ja) * | 2011-01-28 | 2012-08-23 | Mitsubishi Heavy Ind Ltd | 健全性評価装置及びその方法並びにプログラム |
Family Cites Families (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP3108405B2 (ja) * | 1998-07-23 | 2000-11-13 | 核燃料サイクル開発機構 | 機器の診断方法 |
| JP6196093B2 (ja) * | 2012-12-25 | 2017-09-13 | Ntn株式会社 | 軸受装置の振動解析方法、軸受装置の振動解析装置、および転がり軸受の状態監視装置 |
| EP3029449A4 (en) * | 2013-08-01 | 2017-03-15 | NTN Corporation | Bearing-device vibration analysis method, bearing-device vibration analysis device, and rolling-bearing status-monitoring device |
| CN104655425B (zh) * | 2015-03-06 | 2017-05-03 | 重庆大学 | 基于稀疏表示和大间隔分布学习的轴承故障分类诊断方法 |
| CN105021399B (zh) * | 2015-06-26 | 2017-12-05 | 长安大学 | 一种基于单通道信号盲分离滚动轴承的特征提取方法 |
| CN105300692B (zh) * | 2015-08-07 | 2017-09-05 | 浙江工业大学 | 一种基于扩展卡尔曼滤波算法的轴承故障诊断及预测方法 |
| CN105606363B (zh) * | 2016-01-29 | 2017-11-07 | 济南大学 | 一种基于域自适应的轴承故障诊断方法 |
| CN105806604B (zh) * | 2016-03-18 | 2018-10-19 | 北京唐智科技发展有限公司 | 一种机车车辆走行部轴承保持架故障预报警方法 |
| CN106895975B (zh) * | 2017-01-17 | 2019-03-15 | 苏州大学 | 基于Stacked SAE深度神经网络的轴承故障诊断方法 |
| CN106885697B (zh) * | 2017-03-17 | 2019-09-20 | 华东交通大学 | 基于fcm-hmm的滚动轴承的性能退化评估方法 |
-
2017
- 2017-10-25 JP JP2017206597A patent/JP6848813B2/ja active Active
-
2018
- 2018-08-01 CN CN201880052067.6A patent/CN110998272A/zh active Pending
- 2018-08-01 WO PCT/JP2018/028942 patent/WO2019026980A1/ja unknown
- 2018-08-01 EP EP18840775.3A patent/EP3663741A4/en active Pending
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JPH0124246B2 (ja) * | 1981-01-30 | 1989-05-10 | Sumitomo Metal Ind | |
| JP2003050158A (ja) * | 2001-08-06 | 2003-02-21 | Mitsui Eng & Shipbuild Co Ltd | 周期的運動体の異常診断方法および装置 |
| WO2004017038A1 (en) * | 2002-08-14 | 2004-02-26 | Katsumi Hiramatsu | Detector of defects for rotating machinery |
| JP2012159298A (ja) * | 2011-01-28 | 2012-08-23 | Mitsubishi Heavy Ind Ltd | 健全性評価装置及びその方法並びにプログラム |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2021135683A (ja) * | 2020-02-26 | 2021-09-13 | 株式会社東芝 | 学習装置、推論装置、学習方法及び推論方法 |
| WO2022044609A1 (ja) * | 2020-08-24 | 2022-03-03 | 古野電気株式会社 | 船舶航行支援装置、船舶航行支援方法、および、船舶航行支援プログラム |
| WO2022044608A1 (ja) * | 2020-08-24 | 2022-03-03 | 古野電気株式会社 | 船舶航行支援装置、船舶航行支援方法、および、船舶航行支援プログラム |
| CN113899550A (zh) * | 2021-09-28 | 2022-01-07 | 新黎明科技股份有限公司 | 迭代高阶能量算子融合防爆电机轴承故障诊断方法及系统 |
| CN113899550B (zh) * | 2021-09-28 | 2022-08-02 | 新黎明科技股份有限公司 | 迭代高阶能量算子融合防爆电机轴承故障诊断方法及系统 |
Also Published As
| Publication number | Publication date |
|---|---|
| WO2019026980A1 (ja) | 2019-02-07 |
| EP3663741A4 (en) | 2021-05-19 |
| EP3663741A1 (en) | 2020-06-10 |
| CN110998272A (zh) | 2020-04-10 |
| JP6848813B2 (ja) | 2021-03-24 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| JP2019079352A (ja) | 情報処理装置、情報処理方法及びプログラム | |
| Cofre-Martel et al. | Deep convolutional neural network-based structural damage localization and quantification using transmissibility data | |
| JP6919397B2 (ja) | 情報処理装置、情報処理方法及びプログラム | |
| Arora | Comparative study of finite element model updating methods | |
| JP6237774B2 (ja) | 情報処理システム、情報処理方法及びプログラム | |
| US7532988B2 (en) | Virtual load monitoring system and method | |
| JP6679166B2 (ja) | 安全性診断装置、安全性診断方法及び安全性診断プログラム | |
| US9690749B2 (en) | Smart data sampling and data reconstruction | |
| US20220067586A1 (en) | Imaging systems with hybrid learning | |
| WO2019235614A1 (ja) | 関係性分析装置、関係性分析方法および記録媒体 | |
| US20200363287A1 (en) | Damage diagnosing device, damage diagnosing method, and recording medium having damage diagnosing program stored thereon | |
| Karg et al. | Clinical gait analysis: comparing explicit state duration HMMs using a reference-based index | |
| Nestorović et al. | Identification of modal parameters for complex structures by experimental modal analysis approach | |
| WO2019235603A1 (ja) | 関係性分析装置、関係性分析方法および記録媒体 | |
| JP7036233B2 (ja) | 情報処理装置、情報処理方法及びプログラム | |
| JP7014080B2 (ja) | 情報処理装置、情報処理方法及びプログラム | |
| CN110100164A (zh) | 疲劳限度应力确定系统、疲劳限度应力确定装置以及疲劳限度应力确定方法 | |
| Parajuli et al. | Sparsity and biomechanics inspired integration of shape and speckle tracking for cardiac deformation analysis | |
| RU2559401C1 (ru) | Способ диагностирования технического состояния служебных систем летательного аппарата | |
| Voitikova et al. | Linear regression in hemodynamics | |
| JP7363889B2 (ja) | 学習装置、学習方法、コンピュータプログラム及び記録媒体 | |
| JP7099296B2 (ja) | 評価方法、システム構築方法、及び評価システム | |
| Tsinias et al. | Efficient experimental identification of three-dimensional tyre structural properties | |
| US10867369B2 (en) | Image data restoration apparatus and image data restoration method | |
| US10692256B2 (en) | Visualization method, visualization device, and recording medium |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20200603 |
|
| TRDD | Decision of grant or rejection written | ||
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20210202 |
|
| A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20210215 |
|
| R151 | Written notification of patent or utility model registration |
Ref document number: 6848813 Country of ref document: JP Free format text: JAPANESE INTERMEDIATE CODE: R151 |