(57)【要約】【課題】周期的な運動を行う物体の異常診断をより精度よく行うことを目的とする。【解決手段】周期的な運動を行う物体の前記周期的な運動に係る計測データを取得する取得手段と、前記取得手段により取得された前記計測データに基づいて、修正された自己回帰モデルにおける係数を決定する決定手段と、前記決定手段により決定された前記係数に基づいて、前記物体の異常を診断する診断手段と、を有し、前記修正された自己回帰モデルは、前記物体の前記周期的な運動に係る計測データと、前記計測データに対する前記係数と、に基づいて、前記計測データの予測値を導出する式であり、前記決定手段は、前記計測データから導出される自己相関行列を特異値分解することで導出される自己相関行列の固有値に基づいて、前記係数を決定する。【選択図】図１１
2018163135.tif 000002
【特許請求の範囲】
【請求項１】
周期的な運動を行う物体の前記周期的な運動に係る計測データを取得する取得手段と、
前記取得手段により取得された前記計測データに基づいて、修正された自己回帰モデルにおける係数を決定する決定手段と、
前記決定手段により決定された前記係数に基づいて、前記物体の異常を診断する診断手段と、
を有し、
前記修正された自己回帰モデルは、前記計測データの実績値と、前記実績値に対する前記係数と、を用いて、前記計測データの予測値を表す式であり、
前記決定手段は、前記計測データから導出される自己相関行列を特異値分解することで導出される前記自己相関行列の固有値を対角成分とする対角行列と、前記自己相関行列の固有ベクトルを列成分とする直交行列と、から導出される第１の行列を係数行列とし、前記計測データから導出される自己相関ベクトルを定数ベクトルとする方程式を用いて、前記係数を決定し、
前記自己相関ベクトルは、時差が１から前記修正された自己回帰モデルで用いられる前記実績値の数であるｍまでの前記計測データの自己相関を成分とするベクトルであり、
前記自己相関行列は、時差が０からｍ−１までの前記計測データの自己相関を成分とする行列であり、
前記第１の行列は、１以上且つｍ未満の設定された数であるｓに対して、前記自己相関行列のｓ個の固有値と前記対角行列とから導出される第２の行列Σsと、前記ｓ個の固有値と前記直交行列とから導出される第３の行列Ｕsと、から導出される行列ＵsΣsＵsTであり、
前記第２の行列は、前記対角行列の部分行列であって、前記ｓ個の固有値を対角成分とする行列であり、
前記第３の行列は、前記直交行列の部分行列であって、前記ｓ個の固有値に対応する固有ベクトルを列成分ベクトルとする行列である情報処理装置。
【請求項２】
前記ｓ個の固有値は、前記自己相関行列の固有値のうち、値が最大の固有値を含む請求項１記載の情報処理装置。
【請求項３】
前記ｓ個の固有値は、前記自己相関行列の固有値の中から、値が大きいものから順に選ばれた固有値である請求項２記載の情報処理装置。
【請求項４】
前記診断手段による前記物体の異常の診断が成功したか否かを判定する判定手段を更に有し、
前記決定手段は、前記判定手段により、前記診断手段による前記物体の異常の診断が成功していないと判定されると、ｓを設定し直して、前記係数を決定し直し、
前記診断手段は、更に、前記判定手段により、前記診断手段による前記物体の異常の診断が成功していないと判定されると、前記決定手段により決定し直された係数に基づいて、前記周期的な運動を行う物体の異常を診断する請求項１乃至３何れか１項記載の情報処理装置。
【請求項５】
前記自己相関行列は、前記計測データの予測値と当該計測データの予測値に対応する時刻における前記計測データの実測値との二乗誤差を最小化する条件を示す条件式に含まれる行列である請求項１乃至４何れか１項記載の情報処理装置。
【請求項６】
前記診断手段は、前記修正された自己回帰モデルの周波数の分布を示す周波数特性を、前記決定手段により決定された係数を用いて導出し、前記周波数特性に基づいて、前記物体の異常を診断する請求項１乃至５何れか１項記載の情報処理装置。
【請求項７】
前記診断手段は、前記周波数特性におけるピークを示す周波数に基づいて、前記物体における異常の有無を特定し、特定した結果を診断結果とすることで、前記物体の異常を診断する請求項６記載の情報処理装置。
【請求項８】
前記診断手段は、前記周波数特性においてピークを示す周波数に基づいて、前記物体における異常を有する部位を特定し、特定した結果を診断結果とすることで、前記物体の異常を診断する請求項６又は７記載の情報処理装置。
【請求項９】
前記計測データは、前記物体の振動に係る計測データである請求項１乃至８何れか１項記載の情報処理装置。
【請求項１０】
前記診断手段による診断の結果に係る情報を出力する出力手段を更に有する請求項１乃至９何れか１項記載の情報処理装置。
【請求項１１】
前記物体は、鉄道台車に利用される軸受であり、
前記診断手段は、前記決定手段により決定された前記係数に基づいて、前記物体における疵の存在を含む異常を診断する請求項１乃至１０何れか１項記載の情報処理装置。
【請求項１２】
前記周期的な運動は、前記軸受における内輪の回転運動であり、
前記計測データは、前記軸受の内輪の回転数（ｒｐｍ）の１５．６倍以上のサンプリング周波数（Ｈｚ）で計測された前記軸受の振動に係る計測データである請求項１１記載の情報処理装置。
【請求項１３】
情報処理装置が実行する情報処理方法であって、
周期的な運動を行う物体の前記周期的な運動に係る計測データを取得する取得ステップと、
前記取得ステップで取得された前記計測データに基づいて、修正された自己回帰モデルにおける係数を決定する決定ステップと、
前記決定ステップで決定された前記係数に基づいて、前記物体の異常を診断する診断ステップと、
を含み、
前記修正された自己回帰モデルは、前記計測データの実績値と、前記実績値に対する前記係数と、を用いて、前記計測データの予測値を表す式であり、
前記決定ステップでは、前記計測データから導出される自己相関行列を特異値分解することで導出される前記自己相関行列の固有値を対角成分とする対角行列と、前記自己相関行列の固有ベクトルを列成分とする直交行列と、から導出される第１の行列を係数行列とし、前記計測データから導出される自己相関ベクトルを定数ベクトルとする方程式を用いて、前記係数を決定し、
前記自己相関ベクトルは、時差が１から前記修正された自己回帰モデルで用いられる前記実績値の数であるｍまでの前記計測データの自己相関を成分とするベクトルであり、
前記自己相関行列は、時差が０からｍ−１までの前記計測データの自己相関を成分とする行列であり、
前記第１の行列は、１以上且つｍ未満の設定された数であるｓに対して、前記自己相関行列のｓ個の固有値と前記対角行列とから導出される第２の行列Σsと、前記ｓ個の固有値と前記直交行列とから導出される第３の行列Ｕsと、から導出される行列ＵsΣsＵsTであり、
前記第２の行列は、前記対角行列の部分行列であって、前記ｓ個の固有値を対角成分とする行列であり、
前記第３の行列は、前記直交行列の部分行列であって、前記ｓ個の固有値に対応する固有ベクトルを列成分ベクトルとする行列である情報処理方法。
【請求項１４】
コンピュータを、請求項１乃至１２何れか１項記載の情報処理装置の各手段として機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、情報処理装置、情報処理方法及びプログラムに関する。
【背景技術】
【０００２】
周期的な運動を行う物体の異常診断を行う技術として、運動を行っている物体から得られる信号を測定し、その信号の特徴量を算出し、その特徴量から正常であるか異常であるかの判定を行う方法がある。このとき、正常時と異常時との特徴量の差異を明確にすることによって異常診断を精度よく行うために、その物体の運動の基本周波数の整数倍に対応する狭帯域周波数領域に属する成分を測定された信号から抽出し、抽出した成分の信号から得られる自己回帰モデルの係数を算出し、この自己回帰モデルの係数から特徴量を算出する技術がある（特許文献１）。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】特開２００３−５０１５８号公報
【発明の概要】
【発明が解決しようとする課題】
【０００４】
特許文献１に開示された方法では、測定された信号に対して、予めフィルタ処理等を施すことで、特徴量を算出するための信号を抽出する。しかしながら、特許文献１に開示された方法では、通過させるフィルタの周波数領域を予め設定しておく必要があるため、物体の異常を反映している信号成分の周波数が想定外で、予め設定した周波数領域から外れた場合は、正常時と異常時とで特徴量に差異が現れないため、異常を検知出来ない可能性があるという問題がある。
本発明は、以上のような問題点に鑑みてなされたものであり、周期的な運動を行う物体の異常診断をより精度よく行うことができるようにすることを目的とする。
【課題を解決するための手段】
【０００５】
本発明の情報処理装置は、周期的な運動を行う物体の前記周期的な運動に係る計測データを取得する取得手段と、前記取得手段により取得された前記計測データに基づいて、修正された自己回帰モデルにおける係数を決定する決定手段と、前記決定手段により決定された前記係数に基づいて、前記物体の異常を診断する診断手段と、を有し、前記修正された自己回帰モデルは、前記計測データの実績値と、前記実績値に対する前記係数と、を用いて、前記計測データの予測値を表す式であり、前記決定手段は、前記計測データから導出される自己相関行列を特異値分解することで導出される前記自己相関行列の固有値を対角成分とする対角行列と、前記自己相関行列の固有ベクトルを列成分とする直交行列と、から導出される第１の行列を係数行列とし、前記計測データから導出される自己相関ベクトルを定数ベクトルとする方程式を用いて、前記係数を決定し、前記自己相関ベクトルは、時差が１から前記修正された自己回帰モデルで用いられる前記実績値の数であるｍまでの前記計測データの自己相関を成分とするベクトルであり、前記自己相関行列は、時差が０からｍ−１までの前記計測データの自己相関を成分とする行列であり、前記第１の行列は、１以上且つｍ未満の設定された数であるｓに対して、前記自己相関行列のｓ個の固有値と前記対角行列とから導出される第２の行列Σsと、前記ｓ個の固有値と前記直交行列とから導出される第３の行列Ｕsと、から導出される行列ＵsΣsＵsTであり、前記第２の行列は、前記対角行列の部分行列であって、前記ｓ個の固有値を対角成分とする行列であり、前記第３の行列は、前記直交行列の部分行列であって、前記ｓ個の固有値に対応する固有ベクトルを列成分ベクトルとする行列である。
【発明の効果】
【０００６】
本発明によれば、周期的な運動を行う物体の異常診断をより精度よく行うことができる技術を提供することができる。
【図面の簡単な説明】
【０００７】
【図１】図１は、実験の状況を示す図である。
【図２】図２は、軸受の詳細の一例を説明する図である。
【図３】図３は、固有値の分布の一例を示す図である。
【図４】図４は、計測データの波形の一例を示す図である。
【図５】図５は、修正された自己回帰モデルの波形の一例を示す図である。
【図６】図６は、修正された自己回帰モデルの周波数特性の一例を示す図である。
【図７】図７は、修正された自己回帰モデルの波形の一例を示す図である。
【図８】図８は、修正された自己回帰モデルの周波数特性の一例を示す図である。
【図９】図９は、診断システムのシステム構成の一例を示す図である。
【図１０】図１０は、情報処理装置のハードウェア構成の一例を示す図である。
【図１１】図１１は、情報処理装置の処理の一例を示すフローチャートである。
【図１２】図１２は、固有値の分布の一例を示す図である。
【図１３】図１３は、計測データの波形の一例を示す図である。
【図１４】図１４は、修正された自己回帰モデルの波形の一例を示す図である。
【図１５】図１５は、修正された自己回帰モデルの周波数特性の一例を示す図である。
【図１６】図１６は、修正された自己回帰モデルの波形の一例を示す図である。
【図１７】図１７は、修正された自己回帰モデルの周波数特性の一例を示す図である。
【発明を実施するための形態】
【０００８】
＜実施形態１＞
以下、本発明の一実施形態について図面に基づいて説明する。
【０００９】
（本実施形態の処理の概要）
本実施形態では、診断システムが鉄道台車に利用される軸受から測定された振動の信号に基づいて、その軸受の異常を診断する。
診断システムは、軸受の振動に応じた信号を測定し、測定された信号を自己回帰モデルで表現し、自己回帰モデルの係数に関する条件式における行列を特異値分解し、その行列の固有値を求め、求めた固有値のうち最大のものから設定された数のみを用いて、自己回帰モデルの係数を決定する。そして、診断システムは、決定した係数から、自己回帰モデルの周波数特性を求め、求めた周波数特性に基づいて、軸受の異常を診断する。
【００１０】
（軸受異常の診断実験）
鉄道台車に利用されている軸受の異常診断精度の評価のために行った実験について説明する。
図１（ａ）は、本実験の状況を説明する図である。図１（ａ）の状況は、実験室に鉄道台車における車輪の駆動機構が設置されている状況である。駆動モータ１０１は、駆動モータ軸１０２、歯車継手１０３、小歯車軸１００を介して、小歯車１０５を回転させる。小歯車１０５は、回転することで、大歯車１０７を回転させる。それに合わせて、車軸１０９が回転し、車軸１０９に接続される車輪が回転することになる。本実験では、車軸１０９には、車輪の代わりにダイナモ１１０が接続されている。ダイナモ１１０は、疑似的な走行負荷を与えるために接続されている発電機である。本実験では、ダイナモ１１０を回転させるための負荷を、走行の際の負荷として想定する。
小歯車１０５が嵌められた小歯車軸１００と大歯車１０７が嵌められた車軸１０９とは、歯車箱１０４に固定されている。そして、歯車箱１０４には、小歯車軸との間に軸受１０６が装着されており、車軸１０９との間に軸受１０８が装着されている。図１（ｂ）は、歯車箱１０４に取り付けられている軸受を側面から見た様子を示す図である。
【００１１】
歯車箱１０４上の位置１１１は、振動計測装置が設置される位置である。振動計測装置は、加速度センサ、レーザ変位計等のセンサを含み、センサを介して物体の振動を検知し、検知した振動に応じた信号を出力する。位置１１１に設置された振動計測装置が位置１１１に伝達される振動を計測し、計測した振動を示す信号を外部の情報処理装置等に有線又は無線による通信を介して出力する。本実験では、振動計測装置は、位置１１１に設置されるとするが、軸受１０６に起因する振動が伝達される位置であるならば、任意の位置に設置することとしてもよい。
ここで、軸受１０６について、説明する。図２は、軸受１０６の一例の構成を説明する図である。軸受１０６は、外輪、内輪、転動体、保持器の４つの部分から構成される。図２には、軸受１０６の外輪、内輪、転動体、保持器の概要が示されている。軸受１０６は、外輪と内輪との間に保持器により保持された転動体が挟まれる構成となっている。
本実験では、小歯車１０５に装着されている軸受１０６として、疵のない正常な軸受、内輪に疵のある軸受、保持器に疵のある軸受のそれぞれを利用した。それぞれの軸受１０６毎に、駆動モータ１０１を駆動させ、位置１１１に設置された振動計測装置が振動を計測することで、軸受１０６の異常診断を行うための計測データを取得した。本実験では、駆動モータ１０１は、軸受１０６の内輪の回転数が２９５４．７ｒｐｍ（ｒｏｕｎｄ  ｐｅｒ  ｍｉｎｕｔｅ）となるように駆動した。
以上のように、軸受１０６が正常な場合、軸受１０６の内輪に疵がある場合、軸受１０６の保持器に疵がある場合について、位置１１１における振動の計測データを取得した。
【００１２】
次に、振動計測装置が計測した計測データを用いた軸受１０６の異常診断の方法を説明する。以下では、振動計測装置が計測した計測データを、計測データｙとする。
本実験では、計測データｙを、修正された自己回帰モデルで近似し、この修正された自己回帰モデルの係数から特徴量を算出することとした。
ある時刻ｋ（１≦ｋ≦Ｍ）における計測データｙの値をｙkとする。Ｍは、計測データｙがどの時刻までのデータを含むかを示す数であり、予め設定されている。ｙkを近似する自己回帰モデルは、例えば、以下の式１のようになる。式１に示すように、自己回帰モデルとは、時系列データにおけるある時刻ｋ（ｍ＋１≦ｋ≦Ｍ）のデータの予測値ｙ＾k（式において＾はｙの上に付けて表記）を、時系列データにおけるその時刻よりも前の時刻ｋ−ｌ（１≦ｌ≦ｍ）のデータの実績値ｙk-lを用いて表す式である。
【００１３】
【数１】
2018163135.tif 000003
【００１４】
式１におけるαは、自己回帰モデルの係数である。また、ｍは、自己回帰モデルにおいてある時刻ｋにおける計測データｙの値であるｙkを、その時刻よりも前の過去幾つのデータを用いて近似するかを示すＭ未満の整数であり、本実験では、１５００とする。
続いて、最小二乗法を用いて、自己回帰モデルによる予測値ｙ＾kが実測値であるｙkに近似するための条件式を求める。自己回帰モデルによる予測値ｙ＾kが実測値であるｙkに近似するための条件として、式１で与えられるｙ＾kとｙkとの二乗誤差を最小化することを考える。即ち、本実験では、自己回帰モデルによる予測値ｙ＾kをｙkに近似するために最小二乗法を用いる。以下の式２は、計測データと自己回帰モデルによる予測値との二乗誤差を最小にするための条件式である。
【００１５】
【数２】
2018163135.tif 000004
【００１６】
式２より、以下の式３の関係を満たす。
【００１７】
【数３】
2018163135.tif 000005
【００１８】
また、式３を変形（行列表記）することで、以下の式４のようになる。
【００１９】
【数４】
2018163135.tif 000006
【００２０】
式４におけるＲjlは計測データｙの自己相関と呼ばれるもので、以下の式５で定義される値である。このときの｜ｊ−ｌ｜を時差という。
【００２１】
【数５】
2018163135.tif 000007
【００２２】
式４を基に、以下の自己回帰モデルの係数に関する関係式である式６を考える。式６は、自己回帰モデルによる計測データの予測値と、その予測値に対応する時刻における計測データと、の誤差を最小化する条件から導出される方程式で、ユール・ウォーカー（Ｙｕｌｅ−Ｗａｌｋｅｒ）方程式と呼ばれるものである。また、式６は自己回帰モデルの係数から成るベクトルを変数ベクトルとする線形方程式で、式６における左辺の定数ベクトルは、時差が１からｍまでの計測データの自己相関を成分とするベクトルで、以下では、自己相関ベクトルとする。また、式６における右辺の係数行列は、時差が０からｍ−１までの計測データの自己相関を成分とする行列であり、以下では、自己相関行列とする。
【００２３】
【数６】
2018163135.tif 000008
【００２４】
また、式６における右辺の自己相関行列（Ｒjlで構成されるｍ×ｍの行列）を、以下の式７のように、自己相関行列Ｒと表記する。
【００２５】
【数７】
2018163135.tif 000009
【００２６】
一般に、自己回帰モデルの係数を求める際には、式６を係数αについて解くという方法が用いられる。式６では、自己回帰モデルで導出されるある時刻ｋにおける計測データの予測値ｙ＾kが、その時刻ｋにおける計測データの実績値ｙkにできるだけ近づくように係数αを導出する。よって、自己回帰モデルの周波数特性に、各時刻における計測データの実績値ｙkに含まれる多数の周波数成分が含まれる。したがって、例えば、計測データｙに含まれるノイズが多い場合には、軸受１０６の振動に係る信号を抽出できなかったり、軸受１０６の故障の態様に応じた特徴を抽出することができなかったりするという問題が生じる。
そこで、本発明者らは、自己回帰モデルの係数αに乗算される自己相関行列Ｒに着目し、鋭意検討した結果、自己相関行列Ｒの固有値の一部を用いて、計測データｙに含まれるノイズの影響が低減され、軸受の振動に係る信号成分が強調される（ＳＮ比を高める）ように自己相関行列Ｒを書き換えればよいという着想に至った。
以下で、このことの具体例を説明する。
自己相関行列Ｒを特異値分解する。自己相関行列Ｒの要素は、対称であるので、自己相関行列Ｒを特異値分解すると以下の式８のように、直交行列Ｕと、対角行列Σと、直交行列Ｕの転置行列との積となる。
【００２７】
【数８】
2018163135.tif 000010
【００２８】
式８の行列Σは、式９に示すように、対角成分が自己相関行列Ｒの固有値となる対角行列である。対角行列Σの対角成分を、σ11、σ22、・・・、σmmとする。また、行列Ｕは、各列成分ベクトルが自己相関行列Ｒの固有ベクトルとなる直交行列である。直交行列Ｕの列成分ベクトルを、ｕ1、ｕ2、・・・、ｕmとする。自己相関行列Ｒの固有ベクトルｕjに対する固有値がσjjという対応関係が有る。自己相関行列Ｒの固有値は、自己回帰モデルによる計測データの予測値の時間波形に含まれる各周波数の成分の強度に反映する変数である。
【００２９】
【数９】
2018163135.tif 000011
【００３０】
自己相関行列Ｒの特異値分解の結果得られる対角行列Σの対角成分であるσ11、σ22、・・・、σmmの値は、数式の表記を簡略にするために降順とする。これらの自己相関行列Ｒの固有値のうち、最大のものから１以上且つｍ未満の設定された数である使用固有値数ｓ個の固有値を用いて、以下の式１０のように、行列Ｒ’を定義する。行列Ｒ’は、自己相関行列Ｒの固有値のうち使用固有値数ｓ個の固有値用いて自己相関行列Ｒを近似した行列である。
【００３１】
【数１０】
2018163135.tif 000012
【００３２】
式１０における行列Ｕsは、式８の直交行列Ｕの左からｓ個の列成分ベクトル（使用される固有値に対応する固有ベクトル）により構成されるｍ×ｓ行列である。つまり、行列Ｕsは、直交行列Ｕから左のｍ×ｓの要素を切り出して構成される部分行列である。また、ＵsTはＵsの転置行列であり、式８の行列ＵTの上からｓ個の行成分ベクトルにより構成されるｓ×ｍ行列である。式１０における行列Σsは、式８の対角行列Σの左からｓ個の列と上からｓ個の行により構成されるｓ×ｓ行列である。つまり、行列Σsは、対角行列Σから左上のｓ×ｓの要素を切り出して構成される部分行列である。
行列Σs及び行列Ｕsを行列要素表現すれば、以下の式１１のようになる。
【００３３】
【数１１】
2018163135.tif 000013
【００３４】
自己相関行列Ｒの代わりに行列Ｒ’を用いることで、式６の関係式を、以下の式１２のように書き換える。
【００３５】
【数１２】
2018163135.tif 000014
【００３６】
式１２を変形することで、係数αを求める式１３が得られる。式１３によって求められた係数αを用いて、式１により予測値ｙ＾kを算出するモデルを「修正された自己回帰モデル」とする。
これまで対角行列Σの対角成分であるσ11、σ22、・・・、σmmの値を降順として説明したが、係数αの算出過程において対角行列Σの対角成分は降順である必要はなく、その場合には、行列Ｕsは直交行列Ｕから左のｍ×ｓの要素を切り出すのではなく、使用される固有値に対応する列成分ベクトル（固有ベクトル）を切り出して構成される部分行列であり、行列Σsは対角行列Σから左上のｓ×ｓの要素を切り出すのではなく、使用される固有値を対角成分とするように切り出される部分行列である。
式１３は、修正された自己回帰モデルの係数の決定に利用される方程式である。式１３の行列Ｕsは、自己相関行列Ｒの特異値分解により得られる直交行列の部分行列であって、利用される固有値に対応する固有ベクトルを列成分ベクトルとする行列である第３の行列である。また、式１３の行列Σsは、自己相関行列Ｒの特異値分解により得られる対角行列の部分行列であって、利用される固有値を対角成分とする行列である第２の行列である。そして、式１３の行列ＵsΣsＵsTは、行列Σsと行列Ｕsとから導出される行列である第１の行列である。
【００３７】
【数１３】
2018163135.tif 000015
【００３８】
式１３の右辺を計算することにより、修正された自己回帰モデルの係数αが求まる。以上に、修正された自己回帰モデルの係数の導出方法の一例について説明してきたが、その基となる自己回帰モデルの係数の導出については直感的に分かり易いように予測値に対して最小二乗法を用いる方法で説明した。しかしながら、一般的には確率過程という概念を用いて自己回帰モデルを定義し、その係数を導出する方法が知られている。その場合に、自己相関は確率過程（母集団）の自己相関で表現されており、この確率過程の自己相関は時差の関数として表されるものである。従って、本実施形態における計測データの自己相関は、確率過程の自己相関を近似するものであれば他の計算式で算出した値に代えても良く、例えば、Ｒ22〜Ｒmmは時差が０の自己相関であるが、これらをＲ11に置き換えても良い。
次に、求めた係数αで特定される修正された自己回帰モデルの周波数特性を求める方法を説明する。
式１で与えられる計測データの予測値ｙ＾kの予測誤差ｘkが白色雑音になるという性質を使うと、修正された自己回帰モデルは白色雑音である予測誤差ｘkを入力とし、計測データの実測値ｙkを出力とする線形時不変システムとみなせるので、以下の手順で周波数特性の算出式を導出することができる。以下では、修正された自己回帰モデルを線形時不変モデルとみなしたものを、単にシステムと称することにする。予測誤差ｘkは、以下の式１４で表される。
【００３９】
【数１４】
2018163135.tif 000016
【００４０】
式１４の両辺にｚ変換を施すと、以下の式１５が得られる。
【００４１】
【数１５】
2018163135.tif 000017
【００４２】
式１５から、システムのインパルス応答のｚ変換である伝達関数Ｈ（ｚ）は、以下の式１６のように求まる。
【００４３】
【数１６】
2018163135.tif 000018
【００４４】
システムの周波数特性は、正弦波入力に対する出力の振幅と位相との変化として現れ、インパルス応答のフーリエ変換で求められる。言い換えると、ｚが複素平面の単位円上を回転した時の伝達関数Ｈ（ｚ）が周波数特性となる。ここで、式１６におけるｚを、以下の式１７のように置くことを考える。
【００４５】
【数１７】
2018163135.tif 000019
ここで、ｊは虚数単位、ωは角周波数、Ｔは標本間隔である。
【００４６】
その場合、Ｈ（ｚ）の振幅特性（システムの周波数特性）は、以下の式１８のように表すことができる。式１８におけるωＴは、０〜２πの範囲で変化するものとする。
【００４７】
【数１８】
2018163135.tif 000020
【００４８】
本実験では、得られた計測データ毎に、式５と式７を用いて自己相関行列Ｒを求め、式８で表される特異値分解を行うことによって自己相関行列Ｒの固有値を求め、自己相関行列Ｒの固有値の分布を求めた。また、得られた計測データ毎に、計測データの波形を求めた。また、得られた計測データ毎に、上記手順で求められた自己相関行列Ｒの特異値分解から式１１を用いて行列Σsと行列Ｕsとを求め、式５と式１３を用いて修正された自己回帰モデルの係数αを求め、求めた係数αで特定される修正された自己回帰モデルの波形を求めた。そして、得られた計測データ毎に、式１８を用いて、修正された自己回帰モデルの周波数特性等を求めた。
求められた実験結果を、図３〜９を用いて説明する。
【００４９】
図３は、計測データそれぞれについて、上記手順で求められた自己相関行列Ｒの固有値の分布を示す図である。
図３（ａ）のグラフは、軸受１０６に疵がない正常な軸受を用いた場合に計測された計測データを利用した際の自己相関行列Ｒの固有値の分布を示すグラフである。図３（ｂ）のグラフは、軸受１０６に内輪に疵がある軸受を用いた場合に計測された計測データを利用した際の自己相関行列Ｒの固有値の分布を示すグラフである。図３（ｃ）のグラフは、軸受１０６に保持器に疵がある軸受を用いた場合に計測された計測データを利用した際の自己相関行列Ｒの固有値の分布を示すグラフである。
図３の各グラフは、自己相関行列Ｒを特異値分解して得られた固有値σ11〜σmmを昇順に並べ替え、プロットしたグラフであり、横軸が固有値のインデックス、縦軸が固有値の値を示す。
【００５０】
図３（ａ）を見ると、他よりも顕著に高い値をもつ固有値が５つあるのが分かる。その５つの固有値の中でも特に２つの固有値が他の３つよりも高い値となっていることが分かる。
何れの場合にも、自己相関行列Ｒの固有値のうち値が他よりも顕著に高い固有値の数は、固有値全体の数よりも顕著に少ないことが分かる。
【００５１】
図４は、計測データそれぞれの波形を示す図である。
図４（ａ）のグラフは、軸受１０６に疵がない正常な軸受を用いた場合に計測された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの実績値を示す。図４（ｂ）のグラフは、軸受１０６に内輪に疵がある軸受を用いた場合に計測された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの実績値を示す。図４（ｃ）のグラフは、軸受１０６に保持器に疵がある軸受を用いた場合に計測された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの実績値を示す。
図４の各グラフを見ると、例えば、図４（ａ）及び図４（ｃ）が似通った波形をしているのが分かる。そのため、図４のグラフから軸受１０６に生じた異常を診断することは困難である。
【００５２】
図５は、使用固有値数ｓを１として求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示す図である。
図５（ａ）のグラフは、使用固有値数ｓを１として、軸受１０６に疵がない正常な軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図５（ｂ）のグラフは、使用固有値数ｓを１として、軸受１０６に内輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図５（ｃ）のグラフは、使用固有値数ｓを１として、軸受１０６に保持器に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。
図５の各グラフは、図４の各グラフと比べると、ノイズ成分が低減されていることが分かる。
【００５３】
図６は、使用固有値数ｓを１として求められた係数αで特定される修正された自己回帰モデルの周波数特性を示す図である。図６の各グラフは、図５の各波形についての周波数特性を示すグラフである。
図６（ａ）のグラフは、使用固有値数ｓを１として、軸受１０６に疵がない正常な軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフであり、横軸が周波数、縦軸が信号強度を示す。図６（ｂ）のグラフは、使用固有値数ｓを１として、軸受１０６に内輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフであり、横軸が周波数、縦軸が信号強度を示す。図６（ｃ）のグラフは、使用固有値数ｓを１として、軸受１０６に保持器に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフであり、横軸が周波数、縦軸が信号強度を示す。
【００５４】
図６（ａ）のグラフを見ると、周波数８２４Ｈｚの部分にピークが立っているのが分かる。また、図６（ｂ）のグラフを見ると、周波数１６４９Ｈｚの部分にピークが立っているのが分かる。また、図６（ｃ）のグラフを見ると、周波数１６４６Ｈｚの部分にピークが立っているのが分かる。
図６の各グラフにおけるピークが立っている周波数を見ると、軸受１０６に疵がない場合は、軸受１０６に疵がある場合に比べて小さい周波数でピークが立っていることが分かる。即ち、使用固有値数ｓを１として、求められた係数αで特定される修正された自己回帰モデルの周波数特性から、軸受１０６における疵の有無を判断することが可能であることが分かる。
以上より、自己相関行列Ｒの１個の固有値に基づいて求められた係数αで特定される修正された自己回帰モデルの周波数特性に基づいて、軸受１０６の異常診断を行うことができるという知見を得ることができた。
【００５５】
しかし、図６（ｂ）と図６（ｃ）とのグラフに示されるように、軸受１０６の内輪に疵がある場合と、軸受１０６の保持器に疵がある場合と、では、似通った値の周波数において周波数特性のピークが立っている。そのため、使用固有値数ｓを１として、求められた係数αで特定される修正された自己回帰モデルの周波数特性から、軸受１０６に内輪の又は保持器の疵の有無を判断することは可能であるが、内輪の疵であるか、保持器の疵であるか、両方の疵であるか、を判断することは困難である。
そこで、本実験では、更に、使用固有値数ｓの値を５に増やして、自己相関行列Ｒの固有値のうち、最大のものから５つを用いて係数αを求め、求めた係数αで特定される修正された自己回帰モデルの波形、周波数特性を求めた。
【００５６】
図７は、使用固有値数ｓを５として求められた係数αで特定される修正された自己回帰モデルの波形を示す図である。
図７（ａ）のグラフは、使用固有値数ｓを５として、軸受１０６に疵がない正常な軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図７（ｂ）のグラフは、使用固有値数ｓを５として、軸受１０６に内輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図７（ｃ）のグラフは、使用固有値数ｓを５として、軸受１０６に保持器に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。
【００５７】
図８は、使用固有値数ｓを５として求められた係数αで特定される修正された自己回帰モデルの周波数特性を示す図である。図８の各グラフは、図７の各波形についての周波数特性を示すグラフである。
図８（ａ）のグラフは、使用固有値数ｓを５として、軸受１０６に疵がない正常な軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフであり、横軸が周波数、縦軸が信号強度を示す。図８（ｂ）のグラフは、使用固有値数ｓを５として、軸受１０６に内輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフであり、横軸が周波数、縦軸が信号強度を示す。図８（３）のグラフは、使用固有値数ｓを５として、軸受１０６に保持器に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフであり、横軸が周波数、縦軸が信号強度を示す。
【００５８】
図８（ａ）のグラフを見ると、周波数８２１Ｈｚと１０４７Ｈｚの部分にピークが立っているのが分かる。また、図８（ｂ）のグラフを見ると、周波数１６４８Ｈｚ、２０２１Ｈｚ、６４７４Ｈｚの部分にピークが立っているのが分かる。また、図８（ｃ）のグラフを見ると、周波数１６４６Ｈｚ、３２９１Ｈｚ、１７４６Ｈｚの部分にピークが立っているのが分かる。
図８の各グラフにおけるピークと、図６の各グラフにおけるピークと、を見比べてみると、図６のグラフでピークが立っている周波数近辺においては、図８のグラフでもピークが立っているのが分かる。図８のグラフでは、更に、図６のグラフで見られなかったピークが見られるようになっている。
このように、係数αを求めるのに利用される自己相関行列Ｒの固有値を増やすことで、修正された自己回帰モデルによる予測値において、図６でピークが立っていた周波数以外の周波数の成分の信号強度が増加したことが分かる。即ち、係数αを求めるのに利用される自己相関行列Ｒの固有値は、求められた係数αで特定される修正された自己回帰モデルによる計測データの予測値の時間波形に含まれる各周波数の成分の強度と相関があることが見出された。
また、図８（ｂ）、図８（ｃ）のグラフを見てみると、図６（ｂ）、図６（ｃ）のグラフと同様に、１６４０Ｈｚ付近にピークが立っているのが分かる。しかし、図８（ｂ）では、更に、２０２１Ｈｚ、６４７４Ｈｚにおいてもピークが立っており、図８（ｃ）では、更に、３２９１Ｈｚ、１７４６Ｈｚにおいてもピークが立っていることが分かる。このように、図８では、内輪に疵のある軸受１０６と保持器に疵のある軸受１０６とで、ピークが立っている周波数に顕著な違いが見て取れる。即ち、使用固有値数ｓを５にすることで、軸受の疵が、内輪の疵であるか、保持器の疵であるか、両方の疵であるか、を見分けることができる。
【００５９】
このように本発明者らは、使用固有値数ｓが１の場合には見分けることが困難であった疵の存在する部位を、使用固有値数ｓを増やすことで、見分けることができるようになるという知見を得た。一方で、本発明者らは、使用固有値数ｓを増やし過ぎると修正された自己回帰モデルは、実際の計測データに近づいていき、軸受１０６の異常診断に有用でないノイズ成分等についても信号強度が増加されてしまい、式１８で得られる周波数特性にもノイズ成分に起因するピークが生じるという知見を得た。即ち、自己相関行列Ｒの固有値のうち、相対的に値の高い一部の固有値を用いて係数αを求めることで、計測データに含まれる成分のうち、物体の異常診断に有用な成分を近似した修正された自己回帰モデルを求めることができると考えられる。
以上より、係数αを求めるのに利用される自己相関行列Ｒの固有値の数である使用固有値ｓを調整することで、より精度よく、軸受１０６の異常診断を行うことができるという知見を得ることができた。
本実施形態の処理は、以下のような処理である。即ち、本実験で得られた知見を基づいて、物体の振動の計測データから求められる自己相関行列Ｒを特異値分解し、得られる固有値の一部を用いて、計測データを近似する修正された自己回帰モデルの係数を求める。次に、求めた係数で特定される修正された自己回帰モデルから、式１の右辺の計算を行って計測データの予測値を求める訳ではなく、求めた係数を用いて、修正された自己回帰モデルの周波数特性を取得する。そして、取得した周波数特性に基づいて、物体の異常を診断する処理である。ただし、本実施形態の処理に、計測データの予測値（式１の左辺の値）を求める処理を加えてもよい。
【００６０】
（システム構成）
図９は、本実施形態の診断システムのシステム構成の一例を示す図である。診断システムは、周期的な運動を行う物体についての異常診断を行うシステムである。診断システムは、情報処理装置９００、振動計測装置９０１を含む。本実施形態では、診断システムは、鉄道台車に利用されている軸受の異常診断を行う。
情報処理装置９００は、振動計測装置９０１により計測された信号に基づいて、軸受の異常診断を行うパーソナルコンピュータ（ＰＣ）、サーバ装置、タブレット装置等の情報処理装置である。また、情報処理装置９００は、電車に組み込まれたコンピュータ等であってもよい。
振動計測装置９０１は、加速度センサ等のセンサを含み、センサを介して物体の振動を検知し、検知した振動に応じた信号を有線又は無線で情報処理装置９００等の外部の装置に出力する計測装置である。
【００６１】
（情報処理装置の機能構成）
図９を参照しながら、情報処理装置９００が有する機能の一例を説明する。
情報処理装置９００は、取得部９１０、決定部９２０、診断部９３０、判定部９４０、出力部９５０を含む。
≪取得部９１０≫
取得部９１０は、周期的な運動を行う物体について、この周期的な運動に係る計測データｙを取得する。
本実施形態では、取得部９１０は、軸受１０６を回転させ、位置１１１に設置された振動計測装置９０１により計測された振動の計測データｙを取得する。
【００６２】
≪決定部９２０≫
決定部９２０は、取得部９１０により取得された計測データｙに基づいて、修正された自己回帰モデルにおける係数αを決定する。修正された自己回帰モデルは、計測データｙの実績値ｙk-1〜ｙk-mと、この実績値に対する係数α（＝α1〜αm）と、を用いて、計測データｙの予測値ｙ＾kを表す式１である。決定部９２０は、修正された自己回帰モデルにおける係数αを、一般的に知られている自己回帰モデルにおける係数の決定に用いる式６（ユール・ウォーカー方程式）において、式７で表される自己相関行列Ｒの代わりに、自己相関行列Ｒから物体の異常診断に有用な成分を抽出した第１の行列Ｒ’（式１０）を用いる式である式１２を用いて決定する。式６は、時差が０からｍ−１までの計測データｙの自己相関を成分とする自己相関行列Ｒを係数行列とし、時差が１からｍまでの計測データｙの自己相関を成分とする自己相関ベクトルを定数ベクトルとする方程式である。
式６は、式１で算出される計測データｙの予測値ｙ＾kと、計測データｙの予測値ｙ＾kに対応する時刻ｋにおける計測データｙの実測値ｙkとの二乗誤差を最小化する条件を示す条件式として導出することができる。第１の行列Ｒ’は、自己相関行列Ｒを特異値分解（式８）することで導出される自己相関行列Ｒの固有値を対角成分とする対角行列Σと、自己相関行列Ｒの固有ベクトルを列成分とする直交行列Ｕと、から導出される。第１の行列Ｒ’は、自己相関行列Ｒの固有値のうち１以上且つｍ未満の設定された使用固有値数ｓ個の固有値を用いて、対角行列Σの部分行列であって、使用固有値数ｓ個の固有値を対角成分とする行列Σsと、直交行列Ｕの部分行列であって、使用固有値数ｓ個の固有値に対応する固有ベクトルを列成分ベクトルとする行列Ｕsと、から導出される行列ＵsΣsＵsTである。使用固有値数ｓ個の固有値は、自己相関行列Ｒの固有値のうち、値が最大の固有値を含むようにすればよいが、値が大きいものから順に選ばれるのが好ましい。
後述する判定部９４０により、物体の異常の診断が成功していないと判定されると、決定部９２０は、使用固有値数ｓを設定し直し、修正された自己回帰モデルにおける係数αを決定し直す。
【００６３】
≪診断部９３０≫
診断部９３０は、決定部９２０により決定された係数αに基づいて、物体の異常を診断する。診断部９３０は、修正された自己回帰モデルの周波数の分布を示す周波数特性を、決定部９２０により決定された係数αを用いて式１８により導出し、この周波数特性に基づいて、物体の異常を診断してもよい。更に、診断部９３０は、周波数特性におけるピークを示す周波数に基づいて、物体における異常の有無を診断の結果としてもよいし、物体における異常を有する部位を診断の結果としてもよい。
本実施形態では、軸受１０６における疵の有無、又は、疵を有する部位を診断の結果とする。
≪判定部９４０≫
判定部９４０は、診断部９３０による物体の異常の診断が成功したか否かを判定する。判定部９４０により物体の異常の診断が成功していないと判定されると、前述した決定部９２０は、使用固有値数ｓを設定し直し、修正された自己回帰モデルにおける係数αを決定し直す。
≪出力部９５０≫
出力部９５０は、診断部９３０による診断の結果に係る情報を出力する。
【００６４】
（情報処理装置のハードウェア構成）
図１０は、情報処理装置９００のハードウェア構成の一例を示す図である。
情報処理装置９００は、ＣＰＵ１０００、主記憶装置１００１、補助記憶装置１００２、入出力Ｉ／Ｆ１００３を含む。各構成要素は、システムバス１００４を介して、相互に通信可能に接続されている。
ＣＰＵ１０００は、情報処理装置９００を制御する中央演算装置である。主記憶装置１００１は、ＣＰＵ１０００のワークエリアやデータの一時的な保管場所として機能するＲＡＭ（Ｒａｎｄｏｍ  Ａｃｃｅｓｓ  Ｍｅｍｏｒｙ）等の記憶装置である。補助記憶装置１００２は、各種のプログラム、設定データ、振動計測装置９０１から出力される計測データ、診断情報等を記憶するＲＯＭ（Ｒｅａｄ  Ｏｎｌｙ  Ｍｅｍｏｒｙ）、ＨＤＤ（Ｈａｒｄ  Ｄｉｓｋ  Ｄｒｉｖｅ）、ＳＳＤ（Ｓｏｌｉｄ  Ｓｔａｔｅ  Ｄｒｉｖｅ）等の記憶装置である。入出力Ｉ／Ｆ１００３は、振動計測装置９０１等の外部の装置との間での情報のやり取りに利用されるインターフェースである。
ＣＰＵ１０００が、補助記憶装置１００２等に記憶されたプログラムに基づき処理を実行することによって、図９で説明した情報処理装置９００の機能及び図１１で後述するフローチャートの処理等が実現される。
【００６５】
（異常診断処理）
本実施形態では、診断システムは、図１と同様の状況で位置１１１に設置された振動計測装置９０１により計測された振動の計測データに基づいて、軸受１０６の異常を診断することとする。本実施形態では、診断システムは、軸受１０６のどこに疵が存在するかを特定することで、軸受１０６の異常を診断することとするが、軸受１０６に疵が存在するか否かを特定することで、軸受１０６の異常を診断することとしてもよい。本実施形態では、診断システムは、予め振動計測装置９０１により測定された測定データに基づいて、軸受１０６の異常を診断することとする。
図１１は、情報処理装置９００の処理の一例を示すフローチャートである。
Ｓ１１０１において、取得部９１０は、振動計測装置９０１により計測された振動の計測データｙを取得する。取得部９１０は、計測データｙとして、時刻１〜Ｍまでの計測データであるｙ1〜ｙMを取得する。
Ｓ１１０２において、決定部９２０は、Ｓ１１０１で取得した計測データｙと、予め設定された定数Ｍと、修正された自己回帰モデルにおいて、ある時刻のデータを過去の幾つのデータで近似するかを示す数ｍと、に基づいて、式５と式７とを用いて自己相関行列Ｒを生成する。決定部９２０は、予め記憶されているＭ、ｍの情報を読み出すことで、Ｍとｍとを取得する。本実施形態では、ｍの値は、１５００である。また、Ｍは、ｍよりも大きい整数である。
【００６６】
Ｓ１１０３において、決定部９２０は、Ｓ１１０２で生成した自己相関行列Ｒを特異値分解することで、式８の直交行列Ｕと対角行列Σを取得し、対角行列Σから自己相関行列Ｒの固有値σ11〜σmmを取得する。
Ｓ１１０４において、決定部９２０は、修正された自己回帰モデルの係数αを求めるのに利用する自己相関行列Ｒの固有値の数である使用固有値数ｓを１にする。
Ｓ１１０５において、決定部９２０は、自己相関行列Ｒの複数の固有値であるσ11〜σmmのうち、最大のものから使用固有値数ｓだけの固有値σ11〜σssを、修正された自己回帰モデルの係数αを求めるのに利用する自己相関行列Ｒの固有値として選択する。そして、決定部９２０は、計測データｙと、固有値σ11〜σssと、自己相関行列Ｒの特異値分解により得られた直交行列Ｕと、に基づいて、式１３を用いて、修正された自己回帰モデルの係数αを決定する。
Ｓ１１０６において、診断部９３０は、Ｓ１１０５で決定した係数αに基づいて、式１８を用いて、Ｓ１１０５で決定した係数αで特定される修正された自己回帰モデルの周波数特性を取得する。
【００６７】
Ｓ１１０７において、診断部９３０は、Ｓ１１０６で取得した周波数特性に基づいて、軸受１０６の異常を診断する。
本実施形態では、予め、周波数特性におけるピークが立っている周波数と、軸受１０６の状態と、の対応を示す情報であり、軸受１０６の異常診断に用いられる情報である診断情報が、補助記憶装置に記憶されている。診断情報とは、例えば、周波数特性におけるピークが立つ周波数が設定された範囲であれば、軸受１０６が正常であり、周波数特性におけるピークが立つ周波数が設定された範囲外であれば、軸受１０６に疵が存在することを示す情報等である。また、診断情報は、例えば、軸受１０６の特定の部位に疵が存在する状態に対応する周波数特性におけるピークが立つ周波数の範囲を示す情報であってもよい。また、診断情報は、例えば、設定された閾値以上の周波数においてピークが立っていれば、疵が存在する状態であり、また、その閾値未満の周波数においてピークが立っていれば、正常な状態である場合、その閾値の情報であってもよい。
また、図６、図８に示すように、軸受１０６等の診断対象の物体が同じ状態であったとしても、使用固有値数ｓの値によって、修正された自己回帰モデルの周波数特性におけるピークが立つ場所が異なる。そのため、補助記憶装置１００２は、予め、使用固有値数に応じた診断情報を記憶している。
【００６８】
診断部９３０は、使用固有値数ｓの値に応じた診断情報を補助記憶装置１００２から取得する。本実施形態では、診断部９３０は、診断情報として、軸受１０６の状態と、その状態に応じた周波数特性におけるピークが立つ周波数の範囲と、の対応情報を取得する。
診断部９３０は、Ｓ１１０６で取得した周波数特性に基づいて、ピークが立っている周波数を特定する。診断部９３０は、例えば、Ｓ１１０６で取得した周波数特性から、信号強度が設定された閾値以上となる周波数を、ピークが立っている周波数として特定する。ＣＰＵ１０００は、例えば、ピークが立っていると特定した周波数が連続する場合（例えば、１００Ｈｚ、１０１Ｈｚ、１０２Ｈｚでピークが立っている場合）、それらの連続する周波数のうち、信号強度が最も高いものをピークが立っている周波数として特定してもよい。
診断部９３０は、特定したピークが立っている周波数と、取得した診断情報と、に基づいて、軸受１０６の異常を診断する。診断部９３０は、例えば、診断情報が示す軸受１０６の状態と、周波数の範囲と、の対応に基づいて、特定した周波数が、軸受１０６のどの状態に対応するかを特定する。
【００６９】
Ｓ１１０８において、判定部９４０は、Ｓ１１０７での診断結果に基づいて、異常診断が成功したか否かを判定する。判定部９４０は、例えば、Ｓ１１０７で、特定の部位に疵が存在することを示す状態であること、又は正常なことを示す状態であることを特定した場合、異常診断が成功したとして、Ｓ１１１２の処理に進む。判定部９４０は、例えば、Ｓ１１０７で、部位が特定できない疵が存在することを示す状態であることを特定した、又は疵が存在するか否かを特定できなかった場合、異常診断が失敗したとして、Ｓ１１０９の処理に進む。
Ｓ１１０９において、決定部９２０は、使用固有値数ｓの値を変更する。ＣＰＵ１０００は、例えば、現在の使用固有値数ｓの値に１を加えた値を、新たな使用固有値数ｓの値とする。
Ｓ１１１０において、決定部９２０は、使用固有値数ｓの数が設定された閾値よりも大きいか否かを判定する。決定部９２０は、使用固有値数ｓの数が設定された閾値よりも大きいと判定した場合、Ｓ１１１１の処理に進む。また、決定部９２０は、使用固有値数ｓの数が設定された閾値以下であると判定した場合、Ｓ１１０５の処理に進む。これにより、決定部９２０は、値が変更された使用固有値数ｓを用いて、再度、修正された自己回帰モデルの係数αを決定する。
Ｓ１１１１において、決定部９２０は、軸受１０６の異常の診断が失敗したとして、Ｓ１１１２の処理に進む。
Ｓ１１１２において、出力部９５０は、軸受１０６の異常の診断の結果を示す情報を出力する。出力部９５０は、例えば、表示装置に、軸受１０６の異常の診断の結果を示す情報を表示することで出力する。また、出力部９５０は、例えば、軸受１０６の異常の診断の結果を示す情報を、補助記憶装置１００２に記憶することで出力することとしてもよい。また、出力部９５０は、例えば、外部のサーバ装置等の設定された送信先に、軸受１０６の異常の診断の結果を示す情報を送信することで出力することとしてもよい。
【００７０】
（まとめ）
以上、本実施形態では、診断システムは、物体の周期的な運動から計測された計測データから、自己相関行列Ｒを生成し、自己相関行列Ｒを特異値分解し、得られた固有値のうち、最大のものから設定された数を用いて、計測データを近似する修正された自己回帰モデルの係数αを決定した。そして、診断システムは、決定した係数αから、修正された自己回帰モデルの周波数特性を求め、求めた周波数特性に基づいて、軸受１０６の異常を診断した。診断システムは、自己相関行列Ｒの固有値のうち一部の固有値を用いることで、軸受１０６の異常診断に有用な成分が残り、有用でない成分が残らないように、計測データを近似する修正された自己回帰モデルの係数αを決定できる。これにより、診断システムは、計測データの軸受１０６の異常診断に有用な成分に基づいて、より精度よく異常診断を行うことができる。
また、本実施形態の処理により、診断システムは、予め異常診断に用いる信号がどのような周波数の信号であるか等の想定を要することなく、計測データから異常診断に有用な成分を、抽出できる。
【００７１】
＜実施形態２＞
実施形態１では、軸受１０６を構成する部品のうち、内輪、保持器に疵がある場合についてのみ、修正された自己回帰モデルの周波数特性を求める実験を行った（図６、図８を参照）。このため、Ｓ１１０１において取得部９１０が振動の計測データｙを取得する際のサンプリング周波数が、図６、図８のグラフにおいてピークが立つ周波数よりも高い周波数であっても、転動体、外輪に疵がある場合の修正された自己回帰モデルの周波数特性のグラフにおいてピークが立つ周波数が、このサンプリング周波数よりも高い場合には、転動体、外輪の疵を精度よく検出することができない虞がある。
そこで、発明者らは、更に、軸受１０６として転動体に疵のある軸受、外輪に疵のある軸受それぞれを用いて、実施形態１で説明した軸受異常の診断実験と同様の実験を行った。実験結果について説明する。
図１２は、図３と同様に、自己相関行列Ｒの固有値の分布を示す図である。
図１２（ａ）のグラフは、軸受１０６に転動体に疵がある軸受を用いた場合に計測された計測データを利用した際の自己相関行列Ｒの固有値の分布を示すグラフである。図１２（ｂ）のグラフは、軸受１０６に外輪に疵がある軸受を用いた場合に計測された計測データを利用した際の自己相関行列Ｒの固有値の分布を示すグラフである。
図１２の各グラフは、図３と同様、自己相関行列Ｒを特異値分解して得られた固有値σ11〜σmmを昇順に並べ替え、プロットしたグラフであり、横軸が固有値のインデックス、縦軸が固有値の値を示す。
図１２（ａ）、（ｂ）の何れの場合にも、自己相関行列Ｒの固有値のうち値が他よりも顕著に高い固有値の数は、固有値全体の数よりも顕著に少ないことが分かる。
【００７２】
図１３は、図４と同様に、計測データそれぞれの波形を示す図である。
図１３（ａ）のグラフは、軸受１０６に転動体に疵がある軸受を用いた場合に計測された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの実績値を示す。図１３（ｂ）のグラフは、軸受１０６に外輪に疵がある軸受を用いた場合に計測された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの実績値を示す。
図１４は、図５と同様に、使用固有値数ｓを１として求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示す図である。
図１４（ａ）のグラフは、使用固有値数ｓを１として、軸受１０６に転動体に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図１４（ｂ）のグラフは、使用固有値数ｓを１として、軸受１０６に外輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。
図１４の各グラフは、図１３の各グラフと比べると、ノイズ成分が低減されていることが分かる。
【００７３】
図１５は、図６と同様に、使用固有値数ｓを１として求められた係数αで特定される修正された自己回帰モデルの周波数特性を示す図である。図１５の各グラフは、図１４の各波形についての周波数特性を示すグラフである。
図１５（ａ）のグラフは、使用固有値数ｓを１として、軸受１０６に転動体に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフであり、横軸が周波数、縦軸が信号強度を示す。図１５（ｂ）のグラフは、使用固有値数ｓを１として、軸受１０６に外輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフであり、横軸が周波数、縦軸が信号強度を示す。
【００７４】
図１５（ａ）のグラフを見ると、周波数１９８５３Ｈｚの部分にピークが立っているのが分かる。また、図１５（ｂ）のグラフを見ると、周波数２３００７Ｈｚの部分にピークが立っているのが分かる。
図１５の各グラフにおけるピークが立っている周波数を見ると、図１５（ａ）よりも、図１５（ｂ）の方が高い周波数であることが分かる。また、図１５の各グラフにおけるピークが立っている周波数は、図６の各グラフにおけるピークが立っている周波数よりも顕著に高いことが分かる。そのため、使用固有値数ｓを１として、求められた係数αで特定される修正された自己回帰モデルの周波数特性から、軸受１０６の転動体又は外輪の疵の有無を判断することが可能であることが分かる。
【００７５】
図１６は、図７と同様に、使用固有値数ｓを５として求められた係数αで特定される修正された自己回帰モデルの波形を示す図である。
図１６（ａ）のグラフは、使用固有値数ｓを５として、軸受１０６に転動体に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図１６（ｂ）のグラフは、使用固有値数ｓを５として、軸受１０６に外輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。
【００７６】
図１７は、図８と同様に、使用固有値数ｓを５として求められた係数αで特定される修正された自己回帰モデルの周波数特性を示す図である。図１７の各グラフは、図１６の各波形についての周波数特性を示すグラフである。
図１７（ａ）のグラフは、使用固有値数ｓを５として、軸受１０６に転動体に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフであり、横軸が周波数、縦軸が信号強度を示す。図１７（ｂ）のグラフは、使用固有値数ｓを５として、軸受１０６に外輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフであり、横軸が周波数、縦軸が信号強度を示す。
【００７７】
図１７（ａ）のグラフを見ると、周波数１９８５３Ｈｚの部分にピークが立っているのが分かる。また、図１７（ｂ）のグラフを見ると、周波数２２７８５Ｈｚ、２３０２０Ｈｚの部分にピークが立っているのが分かる。
図６、８、１５、１７の各グラフにおけるピークを見てみると、軸受１０６に疵がある場合に生じるピークのうち、最も周波数の高いものは、外輪に疵がある場合に２３０００Ｈｚ付近に生じるピークであることが分かる。
このことから、発明者らは、周波数２３０００Ｈｚのピークを検出可能なサンプリング周波数で、軸受１０６の振動を計測することで、軸受１０６の各部位に生じる疵を検出可能であるとの知見を得た。また、軸受１０６の振動の計測データから得られる修正された自己回帰モデルの周波数特性でピークを示す周波数は、軸受１０６の内輪の回転数に比例するとみられる。修正された自己回帰モデルの周波数特性でピークを示す周波数のうち最大の周波数が２３０００Ｈｚであり、軸受１０６の内輪の回転数は、２９５４．７ｒｐｍである。そのため、発明者らは、軸受１０６には、内部の疵により、最大で、内輪の回転数の約７．８（≒２３０００／２９５４．７）倍の周波数の信号が生じることとなるとの知見を得た。即ち、軸受１０６の内輪の回転数（ｒｐｍ）の７．８倍の周波数（Ｈｚ）のピークを検出可能なサンプリング周波数で、軸受１０６の内輪の回転運動に係る振動を計測することで、軸受１０６の各部位に生じる疵を検出できる。
そこで、本実施形態では、内輪が回転する軸受１０６から、軸受１０６の内輪の回転数（ｒｐｍ）の７．８倍の周波数（Ｈｚ）の信号を検出可能なサンプリング周波数で、測定された軸受１０６の内輪の回転運動に係る振動に係るデータに対して、実施形態１と同様の処理を行うこととする。
【００７８】
本実施形態の診断システムのシステム構成は、実施形態１と同様である。また、情報処理装置９００のハードウェア構成及び機能構成も、実施形態１と同様である。
本実施形態の診断システムの処理のうち、実施形態１と異なる点について説明する。
振動計測装置９０１は、本実施形態では、軸受１０６の内輪の回転数（ｒｐｍ）の７．８倍の周波数（Ｈｚ）の信号を検出可能なサンプリング周波数で、軸受１０６の振動を計測する。振動計測装置９０１は、例えば、ナイキスト周波数（Ｈｚ）が軸受１０６の内輪の回転数（ｒｐｍ）の７．８倍となるように、軸受１０６の内輪の回転数（ｒｐｍ）の１５．６倍のサンプリング周波数（Ｈｚ）で軸受１０６の内輪の回転運動に係る振動を計測する。また、振動計測装置９０１は、例えば、軸受１０６の内輪の回転数（ｒｐｍ）の１５．６倍以上のサンプリング周波数（Ｈｚ）で軸受１０６の内輪の回転運動に係る振動を計測することとしてもよい。
Ｓ１１０１において、取得部９１０は、振動計測装置９０１により軸受１０６の内輪の回転数（ｒｐｍ）の７．８倍の周波数（Ｈｚ）の信号を検出可能なサンプリング周波数で計測された振動の計測データｙを取得する。
【００７９】
以上、本実施形態では、診断システムは、軸受１０６の内輪の回転数（ｒｐｍ）の７．８倍の周波数（Ｈｚ）の信号を検出可能なサンプリング周波数で測定された測定データに基づいて、軸受１０６の異常診断を行うこととした。これにより、診断システムは、軸受１０６の各部位に生じる疵を検出できる。
【００８０】
＜変形例＞
実施形態１、２では、診断システムは、鉄道台車の軸受１０６の異常診断を行うこととした。しかし、診断システムは、歯車等の回転体、振動子等の振動体、バネ等の伸縮体、生物の心臓等の周期的な運動を行う他の物体についての異常診断を行うこととしてもよい。
実施形態１、２では、ｍは、予め設定された１５００という数であるとした。しかし、診断対象に応じて実験を行う等して適切に計測データを近似できる過去のデータの数を特定し、特定した数をｍの値としてもよい。例えば、情報処理装置９００は、診断対象から計測された計測データについて、ある時刻の計測データを、その時刻よりも過去の計測データを複数用いて自己回帰モデル等を利用して近似して、近似した値とその時刻の計測データとの差分が設定された閾値未満となる場合の近似に用いられた過去の計測データの数を、ｍとして決定してもよい。
実施形態１、２では、情報処理装置９００は、使用固有値数ｓの初期値を１として、異常診断に失敗すると使用固有値数ｓの値を変更することとした。しかし、情報処理装置９００は、以下のようにして、使用固有値数ｓの値を決定してもよい。即ち、情報処理装置９００は、まず、自己相関行列Ｒの固有値全ての合計値を求める。そして、情報処理装置９００は、自己相関行列Ｒの固有値のうち、最大のものからａ個の合計が、求めた合計値の設定された割合（例えば９０％）以上となるａのうち、最小のものを使用固有値数ｓとして決定してもよい。
実施形態１、２では、情報処理装置９００は、自己相関行列Ｒの固有値のうち、最大のものから、使用固有値数ｓ個の固有値を用いて、物体の異常診断を行うこととした。しかし、情報処理装置９００は、自己相関行列Ｒの固有値のうち、任意のｓ個の固有値を用いて、物体の異常診断を行うこととしてもよい。例えば、情報処理装置９００は自己相関行列Ｒの固有値のうち、最大のものと、３つ目に大きいものと、の２つを利用して、異常診断することとしてもよい。
実施形態１、２では、診断システムは、予め振動計測装置９０１により測定された測定データに基づいて、軸受１０６の異常を診断することとした。しかし、診断システムは、リアルタイムで振動計測装置９０１により測定されて出力され続けている測定データを用いて、稼働中の軸受１０６の異常を診断することとしてもよい。即ち、情報処理装置９００は振動計測装置９０１から新たな計測データｙが取得される度に、計測データｙ1〜ｙMを更新して、図１１の処理を行うようにしてもよい。
実施形態１、２では、情報処理装置９００は、補助記憶装置１００２に記憶されたプログラムを実行することで、図１１の処理を実現することとした。しかし、情報処理装置９００は、外付けの記憶媒体や外部の記憶サーバ等に記憶されたプログラムを実行することで、図１１の処理を実現することとしてもよい。
また、例えば、上述した情報処理装置９００の機能の一部又は全てをハードウェアとして情報処理装置９００に実装してもよい。
また、以上説明した本発明の実施形態は、何れも本発明を実施するにあたっての具体化の例を示したものに過ぎず、これらによって本発明の技術的範囲が限定的に解釈されてはならないものである。即ち、本発明はその技術思想、またはその主要な特徴から逸脱することなく、様々な形で実施することができる。
【符号の説明】
【００８１】
９００  情報処理装置
９１０  取得部
９２０  決定部
９３０  診断部
９４０  判定部
９５０  出力部
９０１  振動計測装置
１０００  ＣＰＵ
【図１】
2018163135.tif 000021
【図２】
2018163135.tif 000022
【図３】
2018163135.tif 000023
【図４】
2018163135.tif 000024
【図５】
2018163135.tif 000025
【図６】
2018163135.tif 000026
【図７】
2018163135.tif 000027
【図８】
2018163135.tif 000028
【図９】
2018163135.tif 000029
【図１０】
2018163135.tif 000030
【図１１】
2018163135.tif 000031
【図１２】
2018163135.tif 000032
【図１３】
2018163135.tif 000033
【図１４】
2018163135.tif 000034
【図１５】
2018163135.tif 000035
【図１６】
2018163135.tif 000036
【図１７】
2018163135.tif 000037