(57)【要約】  物体の周期的な運動に係る計測データに基づいて、修正された自己回帰モデルにおける係数である修正係数を、予め定められた個数だけ決定し、決定された修正係数それぞれのパターンを示す非負値のデータが格納されたデータ行列を、予め定められた状態に対応する修正係数のパターンを示すデータを格納する非負値の行列である教師基底行列と、データ行列における教師基底の重みを示す非負値の行列である教師重み行列と、データ行列から抽出される基底を格納する非負値の行列である抽出基底行列と、データ行列における抽出基底の重みを示す非負値の行列である抽出重み行列と、に分解し、教師重み行列と、抽出重み行列と、に基づいて、物体が予め定められた状態であるか否かを診断する。
2020145215.tif 000002
【特許請求の範囲】
【請求項１】
周期的な運動を行う物体の前記周期的な運動に係る計測データを取得する取得手段と、
前記取得手段により取得された前記計測データに基づいて、修正された自己回帰モデルにおける係数である修正係数を、予め定められた個数だけ決定する決定手段と、
前記決定手段により決定された前記個数の前記修正係数それぞれのパターンを示す非負値のデータが格納されたデータ行列を、教師基底行列と、教師重み行列と、抽出基底行列と、抽出重み行列と、に分解する分解手段と、
前記分解手段により前記データ行列が分解された前記教師重み行列と、前記抽出重み行列と、に基づいて、前記物体が前記予め定められた状態であるか否かを診断する診断手段と、
を有し、
前記修正された自己回帰モデルは、前記計測データの実績値と、前記実績値に対する前記修正係数と、を用いて、前記計測データの予測値を表す式であり、
前記決定手段は、前記計測データから導出される自己相関行列を特異値分解することで導出される前記自己相関行列の固有値を対角成分とする対角行列と、前記自己相関行列の固有ベクトルを列成分とする直交行列と、から導出される第１の行列を係数行列とし、前記計測データから導出される自己相関ベクトルを定数ベクトルとする方程式を用いて、前記修正係数を決定し、
前記自己相関ベクトルは、時差が１から前記修正された自己回帰モデルで用いられる前記実績値の数であるｍまでの前記計測データの自己相関を成分とするベクトルであり、
前記自己相関行列は、時差が０からｍ−１までの前記計測データの自己相関を成分とする行列であり、
前記第１の行列は、１以上且つｍ未満の設定された数であるｓに対して、前記自己相関行列のｓ個の固有値と前記対角行列とから導出される第２の行列Σｓと、前記ｓ個の固有値と前記直交行列とから導出される第３の行列Ｕｓと、から導出される行列ＵｓΣｓＵｓＴであり、
前記第２の行列は、前記対角行列の部分行列であって、前記ｓ個の固有値を対角成分とする行列であり、
前記第３の行列は、前記直交行列の部分行列であって、前記ｓ個の固有値に対応する固有ベクトルを列成分ベクトルとする行列であり、
前記教師基底行列は、教師基底を格納する非負値の行列であり、
前記教師重み行列は、前記教師基底の重みを示す非負値の行列であり、
前記抽出基底行列は、抽出基底を格納する非負値の行列であり、
前記抽出重み行列は、前記抽出基底の重みを示す非負値の行列であり、
前記教師基底は、前記物体の予め定められた状態に対応する前記修正係数のパターンを示す基底ベクトルの集合であり、
前記抽出基底は、前記物体の予め定められた状態に対応する前記修正係数のパターン以外の前記修正係数のパターンを示す基底ベクトルの集合であり、
前記分解手段は、前記教師基底行列と前記教師重み行列との積と前記抽出基底行列と前記抽出重み行列との積との和と、前記データ行列と、の差分に関するコスト関数の値を適正化するように、前記教師重み行列と前記抽出基底行列と前記抽出重み行列との値を更新することで、前記データ行列を、前記教師基底行列と前記教師重み行列と前記抽出基底行列と前記抽出重み行列とに分解する情報処理装置。
【請求項２】
前記ｓ個の固有値は、前記自己相関行列の固有値のうち、値が最大の固有値を含む請求項１記載の情報処理装置。
【請求項３】
前記ｓ個の固有値は、前記自己相関行列の固有値の中から、値が大きいものから順に選ばれた固有値である請求項２記載の情報処理装置。
【請求項４】
前記自己相関行列は、前記計測データの予測値と当該計測データの予測値に対応する時刻における前記計測データの実測値との二乗誤差を最小化する条件を示す条件式に含まれる行列である請求項１乃至３何れか１項記載の情報処理装置。
【請求項５】
前記コスト関数は、前記教師基底行列と前記教師重み行列との積と前記抽出基底行列と前記抽出重み行列との積との和と、前記データ行列と、の差分を評価する項と、前記教師重み行列のスパース性の向上に関する制約と、前記抽出重み行列のスパース性の向上に関する制約と、前記抽出基底の正規直交性の向上に関する制約と、前記教師基底と前記抽出基底との直交性の向上に関する制約と、のうちの少なくとも１つの制約を示す項とを含む請求項１乃至４何れか１項記載の情報処理装置。
【請求項６】
前記コスト関数の値の適正化は、前記コスト関数の値の最小化また最大化である請求項１乃至５何れか１項記載の情報処理装置。
【請求項７】
前記計測データは、前記物体の振動に係る計測データである請求項１乃至６何れか１項記載の情報処理装置。
【請求項８】
前記診断手段による診断の結果に係る情報を出力する出力手段を更に有する請求項１乃至７何れか１項記載の情報処理装置。
【請求項９】
前記物体は、鉄道台車に利用される軸受であり、
前記診断手段は、前記教師重み行列と、前記抽出重み行列と、に基づいて、前記物体の状態が、前記予め定められた状態である疵のない正常な状態であるか否かを診断する請求項１乃至８何れか１項記載の情報処理装置。
【請求項１０】
情報処理装置が実行する情報処理方法であって、
周期的な運動を行う物体の前記周期的な運動に係る計測データを取得する取得ステップと、
前記取得ステップで取得された前記計測データに基づいて、修正された自己回帰モデルにおける係数である修正係数を、予め定められた個数だけ決定する決定ステップと、
前記決定ステップで決定された前記個数の前記修正係数それぞれのパターンを示す非負値のデータが格納されたデータ行列を、教師基底行列と、教師重み行列と、抽出基底行列と、抽出重み行列と、に分解する分解ステップと、
前記分解ステップでの分解により前記データ行列が分解された前記教師重み行列と、前記抽出重み行列と、に基づいて、前記物体が前記予め定められた状態であるか否かを診断する診断ステップと、
を含み、
前記修正された自己回帰モデルは、前記計測データの実績値と、前記実績値に対する前記修正係数と、を用いて、前記計測データの予測値を表す式であり、
前記決定ステップでは、前記計測データから導出される自己相関行列を特異値分解することで導出される前記自己相関行列の固有値を対角成分とする対角行列と、前記自己相関行列の固有ベクトルを列成分とする直交行列と、から導出される第１の行列を係数行列とし、前記計測データから導出される自己相関ベクトルを定数ベクトルとする方程式を用いて、前記修正係数を決定し、
前記自己相関ベクトルは、時差が１から前記修正された自己回帰モデルで用いられる前記実績値の数であるｍまでの前記計測データの自己相関を成分とするベクトルであり、
前記自己相関行列は、時差が０からｍ−１までの前記計測データの自己相関を成分とする行列であり、
前記第１の行列は、１以上且つｍ未満の設定された数であるｓに対して、前記自己相関行列のｓ個の固有値と前記対角行列とから導出される第２の行列Σｓと、前記ｓ個の固有値と前記直交行列とから導出される第３の行列Ｕｓと、から導出される行列ＵｓΣｓＵｓＴであり、
前記第２の行列は、前記対角行列の部分行列であって、前記ｓ個の固有値を対角成分とする行列であり、
前記第３の行列は、前記直交行列の部分行列であって、前記ｓ個の固有値に対応する固有ベクトルを列成分ベクトルとする行列であり、
前記教師基底行列は、教師基底を格納する非負値の行列であり、
前記教師重み行列は、前記教師基底の重みを示す非負値の行列であり、
前記抽出基底行列は、抽出基底を格納する非負値の行列であり、
前記抽出重み行列は、前記抽出基底の重みを示す非負値の行列であり、
前記教師基底は、前記物体の予め定められた状態に対応する前記修正係数のパターンを示す基底ベクトルの集合であり、
前記抽出基底は、前記物体の予め定められた状態に対応する前記修正係数のパターン以外の前記修正係数のパターンを示す基底ベクトルの集合であり、
前記分解ステップでは、前記教師基底行列と前記教師重み行列との積と前記抽出基底行列と前記抽出重み行列との積との和と、前記データ行列と、の差分に関するコスト関数の値を適正化するように、前記教師重み行列と前記抽出基底行列と前記抽出重み行列との値を更新することで、前記データ行列を、前記教師基底行列と前記教師重み行列と前記抽出基底行列と前記抽出重み行列とに分解する情報処理方法。
【請求項１１】
周期的な運動を行う物体の前記周期的な運動に係る計測データを取得する取得ステップと、
前記取得ステップで取得された前記計測データに基づいて、修正された自己回帰モデルにおける係数である修正係数を、予め定められた個数だけ決定する決定ステップと、
前記決定ステップで決定された前記個数の前記修正係数それぞれのパターンを示す非負値のデータが格納されたデータ行列を、教師基底行列と、教師重み行列と、前抽出基底行列と、抽出重み行列と、に分解する分解ステップと、
前記分解ステップでの分解により前記データ行列が分解された前記教師重み行列と、前記抽出重み行列と、に基づいて、前記物体が前記予め定められた状態であるか否かを診断する診断ステップと、
を含むステップをコンピュータに実行させ、
前記修正された自己回帰モデルは、前記計測データの実績値と、前記実績値に対する前記修正係数と、を用いて、前記計測データの予測値を表す式であり、
前記決定ステップでは、前記計測データから導出される自己相関行列を特異値分解することで導出される前記自己相関行列の固有値を対角成分とする対角行列と、前記自己相関行列の固有ベクトルを列成分とする直交行列と、から導出される第１の行列を係数行列とし、前記計測データから導出される自己相関ベクトルを定数ベクトルとする方程式を用いて、前記修正係数を決定し、
前記自己相関ベクトルは、時差が１から前記修正された自己回帰モデルで用いられる前記実績値の数であるｍまでの前記計測データの自己相関を成分とするベクトルであり、
前記自己相関行列は、時差が０からｍ−１までの前記計測データの自己相関を成分とする行列であり、
前記第１の行列は、１以上且つｍ未満の設定された数であるｓに対して、前記自己相関行列のｓ個の固有値と前記対角行列とから導出される第２の行列Σｓと、前記ｓ個の固有値と前記直交行列とから導出される第３の行列Ｕｓと、から導出される行列ＵｓΣｓＵｓＴであり、
前記第２の行列は、前記対角行列の部分行列であって、前記ｓ個の固有値を対角成分とする行列であり、
前記第３の行列は、前記直交行列の部分行列であって、前記ｓ個の固有値に対応する固有ベクトルを列成分ベクトルとする行列であり、
前記教師基底行列は、教師基底を格納する非負値の行列であり、
前記教師重み行列は、前記教師基底の重みを示す非負値の行列であり、
前記抽出基底行列は、抽出基底を格納する非負値の行列であり、
前記抽出重み行列は、前記抽出基底の重みを示す非負値の行列であり、
前記教師基底は、前記物体の予め定められた状態に対応する前記修正係数のパターンを示す基底ベクトルの集合であり、
前記抽出基底は、前記物体の予め定められた状態に対応する前記修正係数のパターン以外の前記修正係数のパターンを示す基底ベクトルの集合であり、
前記分解ステップでは、前記教師基底行列と前記教師重み行列との積と前記抽出基底行列と前記抽出重み行列との積との和と、前記データ行列と、の差分に関するコスト関数の値を適正化するように、前記教師重み行列と前記抽出基底行列と前記抽出重み行列との値を更新することで、前記データ行列を、前記教師基底行列と前記教師重み行列と前記抽出基底行列と前記抽出重み行列とに分解するプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、情報処理装置、情報処理方法及びプログラムに関する。本願は、２０１９年１月９日に日本に出願された特願２０１９−２０５５号に基づき優先権を主張し、それらの内容を全てここに援用する。
【背景技術】
【０００２】
周期的な運動を行う物体の診断を行う技術として、運動を行っている物体から得られる信号を計測し、その信号の特徴量を算出し、その特徴量から物体の診断を行う方法がある。特許文献１には、信号の周波数特性を自己回帰モデルで同定することで、より高精度なスペクトル特性を取得する技術が開示されている。特許文献２には、以下の技術が開示されている。まず、物体の周期的な運動から計測された計測データから、自己相関行列を生成する。自己相関行列を特異値分解することにより得られた固有値のうち、最大のものから設定された数を用いて、計測データを近似する修正された自己回帰モデルの係数を決定する。決定した係数から、修正された自己回帰モデルの周波数特性を求める。求めた周波数特性に基づいて、軸受の異常を診断する。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】特開２００３−１１６８１１号公報
【特許文献２】特開２０１８−１６３１３５号公報
【発明の概要】
【発明が解決しようとする課題】
【０００４】
周期的な運動を行う物体が予め定められた状態であるか否かをより精度よく診断したいという要望がある。特許文献１、２等の従来技術では、この要望を満たすことができなかった。
本発明は、以上のような問題点に鑑みてなされたものであり、周期的な運動を行う物体が予め定められた状態であるか否かをより精度よく診断することを目的とする。
【課題を解決するための手段】
【０００５】
本発明の情報処理装置は、周期的な運動を行う物体の前記周期的な運動に係る計測データを取得する取得手段と、前記取得手段により取得された前記計測データに基づいて、修正された自己回帰モデルにおける係数である修正係数を、予め定められた個数だけ決定する決定手段と、前記決定手段により決定された前記個数の前記修正係数それぞれのパターンを示す非負値のデータが格納されたデータ行列を、教師基底行列と、教師重み行列と、抽出基底行列と、抽出重み行列と、に分解する分解手段と、前記分解手段により前記データ行列が分解された前記教師重み行列と、前記抽出重み行列と、に基づいて、前記物体が前記予め定められた状態であるか否かを診断する診断手段と、を有し、前記修正された自己回帰モデルは、前記計測データの実績値と、前記実績値に対する前記修正係数と、を用いて、前記計測データの予測値を表す式であり、前記決定手段は、前記計測データから導出される自己相関行列を特異値分解することで導出される前記自己相関行列の固有値を対角成分とする対角行列と、前記自己相関行列の固有ベクトルを列成分とする直交行列と、から導出される第１の行列を係数行列とし、前記計測データから導出される自己相関ベクトルを定数ベクトルとする方程式を用いて、前記修正係数を決定し、前記自己相関ベクトルは、時差が１から前記修正された自己回帰モデルで用いられる前記実績値の数であるｍまでの前記計測データの自己相関を成分とするベクトルであり、前記自己相関行列は、時差が０からｍ−１までの前記計測データの自己相関を成分とする行列であり、前記第１の行列は、１以上且つｍ未満の設定された数であるｓに対して、前記自己相関行列のｓ個の固有値と前記対角行列とから導出される第２の行列Σｓと、前記ｓ個の固有値と前記直交行列とから導出される第３の行列Ｕｓと、から導出される行列ＵｓΣｓＵｓＴであり、前記第２の行列は、前記対角行列の部分行列であって、前記ｓ個の固有値を対角成分とする行列であり、前記第３の行列は、前記直交行列の部分行列であって、前記ｓ個の固有値に対応する固有ベクトルを列成分ベクトルとする行列であり、前記教師基底行列は、教師基底を格納する非負値の行列であり、前記教師重み行列は、前記教師基底の重みを示す非負値の行列であり、前記抽出基底行列は、抽出基底を格納する非負値の行列であり、前記抽出重み行列は、前記抽出基底の重みを示す非負値の行列であり、前記教師基底は、前記物体の予め定められた状態に対応する前記修正係数のパターンを示す基底ベクトルの集合であり、前記抽出基底は、前記物体の予め定められた状態に対応する前記修正係数のパターン以外の前記修正係数のパターンを示す基底ベクトルの集合であり、前記分解手段は、前記教師基底行列と前記教師重み行列との積と前記抽出基底行列と前記抽出重み行列との積との和と、前記データ行列と、の差分に関するコスト関数の値を適正化するように、前記教師重み行列と前記抽出基底行列と前記抽出重み行列との値を更新することで、前記データ行列を、前記教師基底行列と前記教師重み行列と前記抽出基底行列と前記抽出重み行列とに分解する。【０００６】
本発明の情報処理方法は、情報処理装置が実行する情報処理方法であって、周期的な運動を行う物体の前記周期的な運動に係る計測データを取得する取得ステップと、前記取得ステップで取得された前記計測データに基づいて、修正された自己回帰モデルにおける係数である修正係数を、予め定められた個数だけ決定する決定ステップと、前記決定ステップで決定された前記個数の前記修正係数それぞれのパターンを示す非負値のデータが格納されたデータ行列を、教師基底行列と、教師重み行列と、抽出基底行列と、抽出重み行列と、に分解する分解ステップと、前記分解ステップでの分解により前記データ行列が分解された前記教師重み行列と、前記抽出重み行列と、に基づいて、前記物体が前記予め定められた状態であるか否かを診断する診断ステップと、を含み、前記修正された自己回帰モデルは、前記計測データの実績値と、前記実績値に対する前記修正係数と、を用いて、前記計測データの予測値を表す式であり、前記決定ステップでは、前記計測データから導出される自己相関行列を特異値分解することで導出される前記自己相関行列の固有値を対角成分とする対角行列と、前記自己相関行列の固有ベクトルを列成分とする直交行列と、から導出される第１の行列を係数行列とし、前記計測データから導出される自己相関ベクトルを定数ベクトルとする方程式を用いて、前記修正係数を決定し、前記自己相関ベクトルは、時差が１から前記修正された自己回帰モデルで用いられる前記実績値の数であるｍまでの前記計測データの自己相関を成分とするベクトルであり、前記自己相関行列は、時差が０からｍ−１までの前記計測データの自己相関を成分とする行列であり、前記第１の行列は、１以上且つｍ未満の設定された数であるｓに対して、前記自己相関行列のｓ個の固有値と前記対角行列とから導出される第２の行列Σｓと、前記ｓ個の固有値と前記直交行列とから導出される第３の行列Ｕｓと、から導出される行列ＵｓΣｓＵｓＴであり、前記第２の行列は、前記対角行列の部分行列であって、前記ｓ個の固有値を対角成分とする行列であり、前記第３の行列は、前記直交行列の部分行列であって、前記ｓ個の固有値に対応する固有ベクトルを列成分ベクトルとする行列であり、前記教師基底行列は、教師基底を格納する非負値の行列であり、前記教師重み行列は、前記教師基底の重みを示す非負値の行列であり、前記抽出基底行列は、抽出基底を格納する非負値の行列であり、前記抽出重み行列は、前記抽出基底の重みを示す非負値の行列であり、前記教師基底は、前記物体の予め定められた状態に対応する前記修正係数のパターンを示す基底ベクトルの集合であり、前記抽出基底は、前記物体の予め定められた状態に対応する前記修正係数のパターン以外の前記修正係数のパターンを示す基底ベクトルの集合であり、前記分解ステップでは、前記教師基底行列と前記教師重み行列との積と前記抽出基底行列と前記抽出重み行列との積との和と、前記データ行列と、の差分に関するコスト関数の値を適正化するように、前記教師重み行列と前記抽出基底行列と前記抽出重み行列との値を更新することで、前記データ行列を、前記教師基底行列と前記教師重み行列と前記抽出基底行列と前記抽出重み行列とに分解��る。【０００７】
本発明のプログラムは、周期的な運動を行う物体の前記周期的な運動に係る計測データを取得する取得ステップと、前記取得ステップで取得された前記計測データに基づいて、修正された自己回帰モデルにおける係数である修正係数を、予め定められた個数だけ決定する決定ステップと、前記決定ステップで決定された前記個数の前記修正係数それぞれのパターンを示す非負値のデータが格納されたデータ行列を、教師基底行列と、教師重み行列と、前抽出基底行列と、抽出重み行列と、に分解する分解ステップと、前記分解ステップでの分解により前記データ行列が分解された前記教師重み行列と、前記抽出重み行列と、に基づいて、前記物体が前記予め定められた状態であるか否かを診断する診断ステップと、を含むステップをコンピュータに実行させ、前記修正された自己回帰モデルは、前記計測データの実績値と、前記実績値に対する前記修正係数と、を用いて、前記計測データの予測値を表す式であり、前記決定ステップでは、前記計測データから導出される自己相関行列を特異値分解することで導出される前記自己相関行列の固有値を対角成分とする対角行列と、前記自己相関行列の固有ベクトルを列成分とする直交行列と、から導出される第１の行列を係数行列とし、前記計測データから導出される自己相関ベクトルを定数ベクトルとする方程式を用いて、前記修正係数を決定し、前記自己相関ベクトルは、時差が１から前記修正された自己回帰モデルで用いられる前記実績値の数であるｍまでの前記計測データの自己相関を成分とするベクトルであり、前記自己相関行列は、時差が０からｍ−１までの前記計測データの自己相関を成分とする行列であり、前記第１の行列は、１以上且つｍ未満の設定された数であるｓに対して、前記自己相関行列のｓ個の固有値と前記対角行列とから導出される第２の行列Σｓと、前記ｓ個の固有値と前記直交行列とから導出される第３の行列Ｕｓと、から導出される行列ＵｓΣｓＵｓＴであり、前記第２の行列は、前記対角行列の部分行列であって、前記ｓ個の固有値を対角成分とする行列であり、前記第３の行列は、前記直交行列の部分行列であって、前記ｓ個の固有値に対応する固有ベクトルを列成分ベクトルとする行列であり、前記教師基底行列は、教師基底を格納する非負値の行列であり、前記教師重み行列は、前記教師基底の重みを示す非負値の行列であり、前記抽出基底行列は、抽出基底を格納する非負値の行列であり、前記抽出重み行列は、前記抽出基底の重みを示す非負値の行列であり、前記教師基底は、前記物体の予め定められた状態に対応する前記修正係数のパターンを示す基底ベクトルの集合であり、前記抽出基底は、前記物体の予め定められた状態に対応する前記修正係数のパターン以外の前記修正係数のパターンを示す基底ベクトルの集合であり、前記分解ステップでは、前記教師基底行列と前記教師重み行列との積と前記抽出基底行列と前記抽出重み行列との積との和と、前記データ行列と、の差分に関するコスト関数の値を適正化するように、前記教師重み行列と前記抽出基底行列と前記抽出重み行列との値を更新することで、前記データ行列を、前記教師基底行列と前記教師重み行列と前記抽出基底行列と前記抽出重み行列とに分解する。
【図面の簡単な説明】
【０００８】
【図１Ａ】図１Ａは、実験の状況を示す図である。
【図１Ｂ】図１Ｂは、歯車箱に取り付けられている軸受を示す図である。
【図２】図２は、軸受の詳細の一例を説明する図である。
【図３Ａ】図３Ａは、固有値の分布の一例を示す図である。
【図３Ｂ】図３Ｂは、固有値の分布の一例を示す図である。
【図３Ｃ】図３Ｃは、固有値の分布の一例を示す図である。
【図３Ｄ】図３Ｄは、固有値の分布の一例を示す図である。
【図３Ｅ】図３Ｅは、固有値の分布の一例を示す図である。
【図４Ａ】図４Ａは、計測データの波形の一例を示す図である。
【図４Ｂ】図４Ｂは、計測データの波形の一例を示す図である。
【図４Ｃ】図４Ｃは、計測データの波形の一例を示す図である。
【図４Ｄ】図４Ｄは、計測データの波形の一例を示す図である。
【図４Ｅ】図４Ｅは、計測データの波形の一例を示す図である。
【図５Ａ】図５Ａは、修正された自己回帰モデルの波形の一例を示す図である。
【図５Ｂ】図５Ｂは、修正された自己回帰モデルの波形の一例を示す図である。
【図５Ｃ】図５Ｃは、修正された自己回帰モデルの波形の一例を示す図である。
【図５Ｄ】図５Ｄは、修正された自己回帰モデルの波形の一例を示す図である。
【図５Ｅ】図５Ｅは、修正された自己回帰モデルの波形の一例を示す図である。
【図６Ａ】図６Ａは、修正された自己回帰モデルの周波数特性の一例を示す図である。
【図６Ｂ】図６Ｂは、修正された自己回帰モデルの周波数特性の一例を示す図である。
【図６Ｃ】図６Ｃは、修正された自己回帰モデルの周波数特性の一例を示す図である。
【図６Ｄ】図６Ｄは、修正された自己回帰モデルの周波数特性の一例を示す図である。
【図６Ｅ】図６Ｅは、修正された自己回帰モデルの周波数特性の一例を示す図である。
【図７Ａ】図７Ａは、修正された自己回帰モデルの波形の一例を示す図である。
【図７Ｂ】図７Ｂは、修正された自己回帰モデルの波形の一例を示す図である。
【図７Ｃ】図７Ｃは、修正された自己回帰モデルの波形の一例を示す図である。
【図７Ｄ】図７Ｄは、修正された自己回帰モデルの波形の一例を示す図である。
【図７Ｅ】図７Ｅは、修正された自己回帰モデルの波形の一例を示す図である。
【図８Ａ】図８Ａは、修正された自己回帰モデルの周波数特性の一例を示す図である。
【図８Ｂ】図８Ｂは、修正された自己回帰モデルの周波数特性の一例を示す図である。
【図８Ｃ】図８Ｃは、修正された自己回帰モデルの周波数特性の一例を示す図である。
【図８Ｄ】図８Ｄは、修正された自己回帰モデルの周波数特性の一例を示す図である。
【図８Ｅ】図８Ｅは、修正された自己回帰モデルの周波数特性の一例を示す図である。
【図９】図９は、係数のパターンの一例を示す図である。
【図１０】図１０は、学習により特定された係数のパターンの一例を説明する図である。
【図１１】図１１は、各パターンの重みの一例を説明する図である。
【図１２】図１２は、教師基底の一例を示す図である。
【図１３Ａ】図１３Ａは、実験対象の軸受の一例を説明する図である。
【図１３Ｂ】図１３Ｂは、実験対象の軸受の一例を説明する図である。
【図１３Ｃ】図１３Ｃは、実験対象の軸受の一例を説明する図である。
【図１４Ａ】図１４Ａは、データ行列の分解結果の一例を説明する図である。
【図１４Ｂ】図１４Ｂは、データ行列の分解結果の一例を説明する図である。
【図１５Ａ】図１５Ａは、データ行列の分解結果の一例を説明する図である。
【図１５Ｂ】図１５Ｂは、データ行列の分解結果の一例を説明する図である。
【図１６Ａ】図１６Ａは、データ行列の分解結果の一例を説明する図である。
【図１６Ｂ】図１６Ｂは、データ行列の分解結果の一例を説明する図である。
【図１７Ａ】図１７Ａは、データ行列の分解結果の一例を説明する図である。
【図１７Ｂ】図１７Ｂは、データ行列の分解結果の一例を説明する図である。
【図１８】図１８は、診断システムのシステム構成等の一例を示す図である。
【図１９】図１９は、情報処理装置のハードウェア構成の一例を示す図である。
【図２０】図２０は、学習処理の一例を示すフローチャートである。
【図２１】図２１は、診断処理の一例を示すフローチャートである。
【図２２Ａ】図２２Ａは、教師基底と抽出基底との重みの一例を示す図である。
【図２２Ｂ】図２２Ｂは、教師基底と抽出基底との重みの一例を示す図である。
【図２２Ｃ】図２２Ｃは、教師基底と抽出基底との重みの一例を示す図である。
【図２２Ｄ】図２２Ｄは、教師基底と抽出基底との重みの一例を示す図である。
【図２２Ｅ】図２２Ｅは、教師基底と抽出基底との重みの一例を示す図である。
【図２３Ａ】図２３Ａは、教師基底と抽出基底との重みの一例を示す図である。
【図２３Ｂ】図２３Ｂは、教師基底と抽出基底との重みの一例を示す図である。
【図２３Ｃ】図２３Ｃは、教師基底と抽出基底との重みの一例を示す図である。
【図２３Ｄ】図２３Ｄは、教師基底と抽出基底との重みの一例を示す図である。
【図２３Ｅ】図２３Ｅは、教師基底と抽出基底との重みの一例を示す図である。
【図２４Ａ】図２４Ａは、教師基底と抽出基底との重みの一例を示す図である。
【図２４Ｂ】図２４Ｂは、教師基底と抽出基底との重みの一例を示す図である。
【図２４Ｃ】図２４Ｃは、教師基底と抽出基底との重みの一例を示す図である。
【図２４Ｄ】図２４Ｄは、教師基底と抽出基底との重みの一例を示す図である。
【図２５Ａ】図２５Ａは、教師基底と抽出基底との重みの一例を示す図である。
【図２５Ｂ】図２５Ｂは、教師基底と抽出基底との重みの一例を示す図である。
【図２５Ｃ】図２５Ｃは、教師基底と抽出基底との重みの一例を示す図である。
【図２５Ｄ】図２５Ｄは、教師基底と抽出基底との重みの一例を示す図である。
【図２６Ａ】図２６Ａは、教師基底と抽出基底との重みの一例を示す図である。
【図２６Ｂ】図２６Ｂは、教師基底と抽出基底との重みの一例を示す図である。
【図２６Ｃ】図２６Ｃは、教師基底と抽出基底との重みの一例を示す図である。
【図２６Ｄ】図２６Ｄは、教師基底と抽出基底との重みの一例を示す図である。
【図２７Ａ】図２７Ａは、教師基底と抽出基底との重みの一例を示す図である。
【図２７Ｂ】図２７Ｂは、教師基底と抽出基底との重みの一例を示す図である。
【図２７Ｃ】図２７Ｃは、教師基底と抽出基底との重みの一例を示す図である。
【図２７Ｄ】図２７Ｄは、教師基底と抽出基底との重みの一例を示す図である。
【発明を実施するための形態】
【０００９】
＜実施形態１＞
以下、本発明の一実施形態について図面に基づいて説明する。
【００１０】
（本実施形態の処理の概要）
本実施形態では、診断システムが鉄道台車に利用される軸受から計測された振動の信号に基づいて、その軸受の状態を診断する。
診断システムは、軸受の振動に応じた信号を計測する。診断システムは、計測された信号を自己回帰モデルで表現し、自己回帰モデルの係数に関する条件式における行列を特異値分解し、その行列の固有値を求める。診断システムは、求めた固有値のうち最大のものから設定された数のみを用いて、自己回帰モデルの係数を決定する。そして、診断システムは、決定した係数のパターンに基づいて、軸受の状態を診断する。
【００１１】
（実験１）
鉄道台車に利用されている軸受の状態診断精度の評価のために行った実験について説明する。
図１Ａは、本実験の状況を説明する図である。図１Ａの状況は、実験室に鉄道台車における車輪の駆動機構が設置されている状況である。駆動モータ１０１は、駆動モータ軸１０２、歯車継手１０３、小歯車軸１００を介して、小歯車１０５を回転させる。小歯車１０５は、回転することで、大歯車１０７を回転させる。それに合わせて、車軸１０９が回転し、車軸１０９に接続される車輪が回転することになる。本実験では、車軸１０９には、車輪の代わりにダイナモ１１０が接続されている。ダイナモ１１０は、疑似的な走行負荷を与えるために接続されている発電機である。本実験では、ダイナモ１１０を回転させるための負荷を、走行の際の負荷として想定する。
【００１２】
小歯車１０５が嵌められた小歯車軸１００と大歯車１０７が嵌められた車軸１０９とは、歯車箱１０４に固定されている。そして、歯車箱１０４には、小歯車軸との間に軸受１０６が装着されており、車軸１０９との間に軸受１０８が装着されている。図１Ｂは、歯車箱１０４に取り付けられている軸受１０６、１０８を側面から見た様子を示す図である。
【００１３】
歯車箱１０４上の振動計測位置１１１は、振動計測装置が設置される位置である。振動計測装置は、加速度センサ、レーザ変位計等のセンサを含み、センサを介して物体の振動を検知し、検知した振動に応じた信号を出力する。振動計測位置１１１に設置された振動計測装置が振動計測位置１１１に伝達される振動を計測する。振動計測装置は、計測した振動を示す信号を外部の情報処理装置等に有線又は無線による通信を介して出力する。振動計測装置の計測周波数は、２０４８００Ｈｚである。ただし、他の例として、振動計測装置の計測周波数は、４０９６００Ｈｚ、１０２４００Ｈｚ等の他の周波数であることとしてもよい。本実験では、振動計測装置は、振動計測位置１１１に設置されるとするが、軸受１０６に起因する振動が伝達される位置であるならば、任意の位置に設置することとしてもよい。
【００１４】
ここで、軸受１０６について、説明する。図２は、軸受１０６の一例の構成を説明する図である。軸受１０６は、外輪、内輪、転動体、保持器の４つの部分から構成される。図２には、軸受１０６の外輪、内輪、転動体、保持器の概要が示されている。軸受１０６は、外輪と内輪との間に保持器により保持された転動体が挟まれる構成となっている。
本実験では、小歯車１０５に装着されている軸受１０６として、疵のない正常な軸受、内輪に疵のある軸受、保持器に疵のある軸受のそれぞれを利用した。それぞれの軸受１０６毎に、駆動モータ１０１を駆動させ、振動計測位置１１１に設置された振動計測装置が振動を計測することで、軸受１０６の状態診断を行うための計測データを取得した。本実験では、駆動モータ１０１は、軸受１０６の内輪の回転数が２９５４．７ｒｐｍ（ｒｏｕｎｄ  ｐｅｒ  ｍｉｎｕｔｅ）となるように駆動した。
【００１５】
以上のように、軸受１０６が正常な場合、軸受１０６の内輪に疵がある場合、軸受１０６の保持器に疵がある場合、軸受１０６の転動体に疵がある場合、軸受１０６の外輪に疵がある場合について、振動計測位置１１１における振動の計測データを取得した。
【００１６】
次に、振動計測装置が計測した計測データを用いた軸受１０６の状態診断の方法を説明する。以下では、振動計測装置が計測した計測データを、計測データｙとする。
【００１７】
本実験では、計測データｙを、修正された自己回帰モデルで近似し、この修正された自己回帰モデルの係数から特徴量を算出することとした。
【００１８】
ある時刻ｋ（ｋ：１≦ｋ≦Ｍの整数）における計測データｙの値をｙｋとする。Ｍは、計測データｙがどの時刻までのデータを含むかを示す数であり、予め設定されている。Ｍは、例えば、後述する自己相関（式５で定義されるＲｊｌ）が統計量として十分な信頼性を得るように、事前に試行錯誤して決定すればよい。ｙｋを近似する自己回帰モデルは、例えば、以下の式１のようになる。式１に示すように、自己回帰モデルとは、時系列データにおけるある時刻ｋ（ｍ＋１≦ｋ≦Ｍ）のデータの予測値ｙ＾ｋ（式において＾はｙの上に付けて表記）を、時系列データにおけるその時刻よりも前の時刻ｋ−ｌ（１≦ｌ≦ｍ）のデータの実績値ｙｋ−ｌを用いて表す式である。
【００１９】
【数１】
2020145215.tif 000003
【００２０】
式１におけるαは、自己回帰モデルの係数である。また、ｍは、自己回帰モデルにおいてある時刻ｋにおける計測データｙの値であるｙｋを、その時刻よりも前の過去幾つのデータを用いて近似するかを示すＭ未満の整数である。ここでは、ｍを１５００とする。
【００２１】
続いて、最小二乗法を用いて、自己回帰モデルによる予測値ｙ＾ｋが実測値であるｙｋに近似するための条件式を求める。自己回帰モデルによる予測値ｙ＾ｋが実測値であるｙｋに近似するための条件として、式１で与えられるｙ＾ｋとｙｋとの二乗誤差を最小化することを考える。即ち、本実験では、自己回帰モデルによる予測値ｙ＾ｋをｙｋに近似するために最小二乗法を用いる。以下の式２は、計測データと自己回帰モデルによる予測値との二乗誤差を最小にするための条件式である。
【００２２】
【数２】
2020145215.tif 000004
【００２３】
式２より、以下の式３の関係を満たす。
【００２４】
【数３】
2020145215.tif 000005
【００２５】
また、式３を変形（行列表記）することで、以下の式４のようになる。
【００２６】
【数４】
2020145215.tif 000006
【００２７】
式４におけるＲｊｌは計測データｙの自己相関と呼ばれるもので、以下の式５で定義される値である。このときの｜ｊ−ｌ｜を時差という。
【００２８】
【数５】
2020145215.tif 000007
【００２９】
式４を基に、以下の自己回帰モデルの係数に関する関係式である式６を考える。式６は、自己回帰モデルによる計測データの予測値と、その予測値に対応する時刻における計測データと、の誤差を最小化する条件から導出される方程式である。式６は、ユール・ウォーカー（Ｙｕｌｅ−Ｗａｌｋｅｒ）方程式と呼ばれるものである。また、式６は自己回帰モデルの係数から成るベクトルを変数ベクトルとする線形方程式である。式６における左辺の定数ベクトルは、時差が１からｍまでの計測データの自己相関を成分とするベクトルである。以下では、式６における左辺の定数ベクトルを自己相関ベクトルとする。また、式６における右辺の係数行列は、時差が０からｍ−１までの計測データの自己相関を成分とする行列である。以下では、式６における右辺の係数行列を自己相関行列とする。
【００３０】
【数６】
2020145215.tif 000008
【００３１】
また、式６における右辺の自己相関行列（Ｒｊｌで構成されるｍ×ｍの行列）を、以下の式７のように、自己相関行列Ｒと表記する。
【００３２】
【数７】
2020145215.tif 000009
【００３３】
一般に、自己回帰モデルの係数を求める際には、式６を係数αについて解くという方法が用いられる。式６では、自己回帰モデルで導出されるある時刻ｋにおける計測データの予測値ｙ＾ｋが、その時刻ｋにおける計測データの実績値ｙｋにできるだけ近づくように係数αが導出される。よって、自己回帰モデルの周波数特性に、各時刻における計測データの実績値ｙｋに含まれる多数の周波数成分が含まれる。したがって、例えば、計測データｙに含まれるノイズが多い場合には、軸受１０６の振動に係る信号を抽出できなかったり、軸受１０６の故障の態様に応じた特徴を抽出することができなかったりするという問題が生じる。
【００３４】
そこで、本発明者らは、自己回帰モデルの係数αに乗算される自己相関行列Ｒに着目し、鋭意検討した結果、自己相関行列Ｒの固有値の一部を用いて、計測データｙに含まれるノイズの影響が低減され、軸受の振動に係る信号成分が強調される（ＳＮ比を高める）ように自己相関行列Ｒを書き換えればよいという着想に至った。
【００３５】
以下で、このことの具体例を説明する。
自己相関行列Ｒを特異値分解する。自己相関行列Ｒの成分は、対称であるので、自己相関行列Ｒを特異値分解すると以下の式８のように、直交行列Ｕと、対角行列Σと、直交行列Ｕの転置行列との積となる。
【００３６】
【数８】
2020145215.tif 000010
【００３７】
式８の行列Σは、式９に示すように、対角成分が自己相関行列Ｒの固有値となる対角行列である。対角行列Σの対角成分を、σ１１、σ２２、・・・、σｍｍとする。また、行列Ｕは、各列成分ベクトルが自己相関行列Ｒの固有ベクトルとなる直交行列である。直交行列Ｕの列成分ベクトルを、ｕ１、ｕ２、・・・、ｕｍとする。自己相関行列Ｒの固有ベクトルｕｊに対する固有値がσｊｊという対応関係が有る。自己相関行列Ｒの固有値は、自己回帰モデルによる計測データの予測値の時間波形に含まれる各周波数の成分の強度に反映する変数である。
【００３８】
【数９】
2020145215.tif 000011
【００３９】
自己相関行列Ｒの特異値分解の結果得られる対角行列Σの対角成分であるσ１１、σ２２、・・・、σｍｍの値は、数式の表記を簡略にするために降順とする。これらの自己相関行列Ｒの固有値のうち、最大のものから１以上且つｍ未満の設定された数である使用固有値数ｓ個の固有値を用いて、以下の式１０のように、行列Ｒ’を定義する。行列Ｒ’は、自己相関行列Ｒの固有値のうち使用固有値数ｓ個の固有値用いて自己相関行列Ｒを近似した行列である。
【００４０】
【数１０】
2020145215.tif 000012
【００４１】
式１０における行列Ｕｓは、式８の直交行列Ｕの左からｓ個の列成分ベクトル（使用される固有値に対応する固有ベクトル）により構成されるｍ×ｓ行列である。つまり、行列Ｕｓは、直交行列Ｕから左のｍ×ｓの成分を切り出して構成される部分行列である。また、ＵｓＴはＵｓの転置行列であり、式８の行列ＵＴの上からｓ個の行成分ベクトルにより構成されるｓ×ｍ行列である。式１０における行列Σｓは、式８の対角行列Σの左からｓ個の列と上からｓ個の行により構成されるｓ×ｓ行列である。つまり、行列Σｓは、対角行列Σから左上のｓ×ｓの成分を切り出して構成される部分行列である。
行列Σｓ及び行列Ｕｓを行列成分表現すれば、以下の式１１のようになる。
【００４２】
【数１１】
2020145215.tif 000013
【００４３】
自己相関行列Ｒの代わりに行列Ｒ’を用いることで、式６の関係式を、以下の式１２のように書き換える。
【００４４】
【数１２】
2020145215.tif 000014
【００４５】
式１２を変形することで、係数αを求める式１３が得られる。式１３によって求められた係数αを用いて、式１により予測値ｙ＾ｋを算出するモデルを「修正された自己回帰モデル」とする。
【００４６】
これまで対角行列Σの対角成分であるσ１１、σ２２、・・・、σｍｍの値を降順として説明したが、係数αの算出過程において対角行列Σの対角成分は降順である必要はない。対角行列Σの対角成分は降順でない場合には、行列Ｕｓは直交行列Ｕから左のｍ×ｓの成分を切り出すのではなく、使用される固有値に対応する列成分ベクトル（固有ベクトル）を切り出して構成される部分行列である。また、対角行列Σの対角成分は降順でない場合には、行列Σｓは対角行列Σから左上のｓ×ｓの成分を切り出すのではなく、使用される固有値を対角成分とするように切り出される部分行列である。
【００４７】
式１３は、修正された自己回帰モデルの係数の決定に利用される方程式である。式１２、１３の行列Ｕｓは、自己相関行列Ｒの特異値分解により得られる直交行列の部分行列であって、利用される固有値に対応する固有ベクトルを列成分ベクトルとする行列である第３の行列である。また、式１２、１３の行列Σｓは、自己相関行列Ｒの特異値分解により得られる対角行列の部分行列であって、利用される固有値を対角成分とする行列である第２の行列である。そして、式１２、１３の行列ＵｓΣｓＵｓＴは、行列Σｓと行列Ｕｓとから導出される行列である第１の行列である。
【００４８】
【数１３】
2020145215.tif 000015
【００４９】
式１３の右辺を計算することにより、修正された自己回帰モデルの係数αが求まる。以上に、修正された自己回帰モデルの係数の導出方法の一例について説明してきたが、その基となる自己回帰モデルの係数の導出については直感的に分かり易いように予測値に対して最小二乗法を用いる方法で説明した。しかしながら、一般的には確率過程という概念を用いて自己回帰モデルを定義し、その係数を導出する方法が知られている。その場合に、自己相関は確率過程（母集団）の自己相関で表現される。この確率過程の自己相関は時差の関数として表されるものである。従って、本実施形態における計測データの自己相関は、確率過程の自己相関を近似するものであれば他の計算式で算出した値に代えても良い。例えば、Ｒ２２〜Ｒｍｍは時差が０の自己相関であるが、これらをＲ１１に置き換えてもよい。
【００５０】
本実験では、以下の４つの場合のそれぞれについて、得られた計測データ毎に、式５と式７を用いて自己相関行列Ｒを求め、式８で表される特異値分解を行うことによって自己相関行列Ｒの固有値を求め、自己相関行列Ｒの固有値の分布を求めた。４つの場合は、軸受１０６が正常な場合、軸受１０６の内輪に疵がある場合、軸受１０６の保持器に疵がある場合、および軸受１０６の転動体に疵がある場合、軸受１０６の外輪に疵がある場合である。
【００５１】
次に、求めた係数αで特定される修正された自己回帰モデルの周波数特性を求める方法を説明する。
式１で与えられる計測データの予測値ｙ＾ｋの予測誤差ｘｋが白色雑音になるという性質を使うと、修正された自己回帰モデルは白色雑音である予測誤差ｘｋを入力とし、計測データの実測値ｙｋを出力とする線形時不変システムとみなせる。従って、以下の手順で周波数特性の算出式を導出することができる。以下では、修正された自己回帰モデルを線形時不変モデルとみなしたものを、単にシステムと称することにする。予測誤差ｘｋは、以下の式１４で表される。
【００５２】
【数１４】
2020145215.tif 000016
【００５３】
式１４の両辺にｚ変換を施すと、以下の式１５が得られる。
【００５４】
【数１５】
2020145215.tif 000017
【００５５】
式１５から、システムのインパルス応答のｚ変換である伝達関数Ｈ（ｚ）は、以下の式１６のように求まる。
【００５６】
【数１６】
2020145215.tif 000018
【００５７】
システムの周波数特性は、正弦波入力に対する出力の振幅と位相との変化として現れ、インパルス応答のフーリエ変換で求められる。言い換えると、ｚが複素平面の単位円上を回転した時の伝達関数Ｈ（ｚ）が周波数特性となる。ここで、式１６におけるｚを、以下の式１７のように置くことを考える。
【００５８】
【数１７】
2020145215.tif 000019
ここで、ｊは虚数単位、ωは角周波数、Ｔは標本間隔である。
【００５９】
その場合、Ｈ（ｚ）の振幅特性（システムの周波数特性）は、以下の式１８のように表すことができる。式１８におけるωＴは、０〜２πの範囲で変化するものとする。
【００６０】
【数１８】
2020145215.tif 000020
【００６１】
本実験では、得られた計測データ毎に、式５と式７を用いて自己相関行列Ｒを求め、式８で表される特異値分解を行うことによって自己相関行列Ｒの固有値を求め、自己相関行列Ｒの固有値の分布を求めた。また、得られた計測データ毎に、計測データの波形を求めた。また、得られた計測データ毎に、上記手順で求められた自己相関行列Ｒの特異値分解から式１１を用いて行列Σｓと行列Ｕｓとを求め、式５と式１３を用いて修正された自己回帰モデルの係数αを求め、求めた係数αで特定される修正された自己回帰モデルの波形を求めた。そして、得られた計測データ毎に、式１８を用いて、修正された自己回帰モデルの周波数特性等を求めた。
求められた実験結果を、図３Ａ〜図８Ｅを用いて説明する。
【００６２】
図３Ａ〜図３Ｅは、計測データそれぞれについて、上記手順で求められた自己相関行列Ｒの固有値の分布を示す図である。
図３Ａのグラフは、軸受１０６に疵がない正常な軸受を用いた場合に計測された計測データを利用した際の自己相関行列Ｒの固有値の分布を示すグラフである。図３Ｂのグラフは、軸受１０６に内輪に疵がある軸受を用いた場合に計測された計測データを利用した際の自己相関行列Ｒの固有値の分布を示すグラフである。図３Ｃのグラフは、軸受１０６に保持器に疵がある軸受を用いた場合に計測された計測データを利用した際の自己相関行列Ｒの固有値の分布を示すグラフである。図３Ｄのグラフは、軸受１０６に転動体に疵がある軸受を用いた場合に計測された計測データを利用した際の自己相関行列Ｒの固有値の分布を示すグラフである。図３Ｅのグラフは、軸受１０６に外輪に疵がある軸受を用いた場合に計測された計測データを利用した際の自己相関行列Ｒの固有値の分布を示すグラフである。
【００６３】
図３Ａ〜図３Ｅの各グラフは、自己相関行列Ｒを特異値分解して得られた固有値σ１１〜σｍｍを昇順に並べ替え、プロットしたグラフであり、横軸が固有値のインデックス、縦軸が固有値の値を示す。
【００６４】
図３Ａを見ると、他よりも顕著に高い値をもつ固有値が５つあるのが分かる。その５つの固有値の中でも特に２つの固有値が他の３つよりも高い値となっていることが分かる。
何れの場合にも、自己相関行列Ｒの固有値のうち値が他よりも顕著に高い固有値の数は、固有値全体の数よりも顕著に少ないことが分かる。
【００６５】
また、本実験では、得られた計測データ毎に、計測データの波形を求めた。
図４Ａ〜図４Ｅは、計測データそれぞれの波形を示す図である。
図４Ａのグラフは、軸受１０６に疵がない正常な軸受を用いた場合に計測された計測データの波形を示すグラフである。図４Ａのグラフでは、横軸が時間、縦軸が対応する時間における計測データの実績値を示す。図４Ｂのグラフは、軸受１０６に内輪に疵がある軸受を用いた場合に計測された計測データの波形を示すグラフである。図４Ｂのグラフでは、横軸が時間、縦軸が対応する時間における計測データの実績値を示す。図４Ｃのグラフは、軸受１０６に保持器に疵がある軸受を用いた場合に計測された計測データの波形を示すグラフである。図４Ｃのグラフでは、横軸が時間、縦軸が対応する時間における計測データの実績値を示す。図４Ｄのグラフは、軸受１０６に転動体に疵がある軸受を用いた場合に計測された計測データの波形を示すグラフである。図４Ｄのグラフでは、横軸が時間、縦軸が対応する時間における計測データの実績値を示す。図４Ｅのグラフは、軸受１０６に外輪に疵がある軸受を用いた場合に計測された計測データの波形を示すグラフである。図４Ｅのグラフでは、横軸が時間、縦軸が対応する時間における計測データの実績値を示す。
【００６６】
図４の各グラフを見ると、例えば、図４Ａ及び図４Ｃが似通った波形をしているのが分かる。また、例えば、図４Ｄ及び図４Ｅが似通った波形をしているのも分かる。そのため、図４Ａ〜図４Ｅのグラフから軸受１０６に生じた異常の状態を診断することは困難である。
【００６７】
図５Ａ〜図５Ｅは、使用固有値数ｓを１として求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示す図である。
【００６８】
図５Ａのグラフは、使用固有値数ｓを１として、軸受１０６に疵がない正常な軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフである。図５Ａのグラフでは、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図５Ｂのグラフは、使用固有値数ｓを１として、軸受１０６に内輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフである。図５Ｂのグラフでは、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図５Ｃのグラフは、使用固有値数ｓを１として、軸受１０６に保持器に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフである。図５Ｃのグラフでは、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図５Ｄのグラフは、使用固有値数ｓを１として、軸受１０６に転動体に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフである。図５Ｄのグラフでは、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図５Ｅのグラフは、使用固有値数ｓを１として、軸受１０６に外輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフである。図５Ｅのグラフでは、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。
図５Ａ〜図５Ｅの各グラフは、図４Ａ〜図４Ｅの各グラフと比べると、ノイズ成分が低減されていることが分かる。
【００６９】
図６Ａ〜図６Ｅは、使用固有値数ｓを１として求められた係数αで特定される修正された自己回帰モデルの周波数特性を示す図である。図６Ａ〜図６Ｅの各グラフは、図５Ａ〜図５Ｅの各波形についての周波数特性を示すグラフである。
【００７０】
図６Ａのグラフは、使用固有値数ｓを１として、軸受１０６に疵がない正常な軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフである。図６Ａのグラフでは、横軸が周波数、縦軸が信号強度を示す。図６Ｂのグラフは、使用固有値数ｓを１として、軸受１０６に内輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフである。図６Ｂのグラフでは、横軸が周波数、縦軸が信号強度を示す。図６Ｃのグラフは、使用固有値数ｓを１として、軸受１０６に保持器に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフである。図６Ｃのグラフでは、横軸が周波数、縦軸が信号強度を示す。図６Ｄのグラフは、使用固有値数ｓを１として、軸受１０６に転動体に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフである。図６Ｄのグラフでは、横軸が周波数、縦軸が信号強度を示す。図６Ｅのグラフは、使用固有値数ｓを１として、軸受１０６に外輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフである。る。図６Ｅのグラフでは、横軸が周波数、縦軸が信号強度を示す。
【００７１】
図６Ａのグラフを見ると、周波数８２４Ｈｚの部分にピークが立っているのが分かる。また、図６Ｂのグラフを見ると、周波数１６４９Ｈｚの部分にピークが立っているのが分かる。また、図６Ｃのグラフを見ると、周波数１６４６Ｈｚの部分にピークが立っているのが分かる。
【００７２】
図６Ｄのグラフを見ると、周波数１９８５３Ｈｚの部分にピークが立っているのが分かる。また、図６Ｅのグラフを見ると、周波数２３００７Ｈｚの部分にピークが立っているのが分かる。
【００７３】
図６Ａ〜図６Ｅの各グラフにおけるピークが立っている周波数を見ると、軸受１０６に疵がない場合は、軸受１０６に疵がある場合に比べて低い周波数でピークが立っていることが分かる。即ち、使用固有値数ｓを１として、求められた係数αで特定される修正された自己回帰モデルの周波数特性から、軸受１０６における疵の有無を判断することが可能であることが分かる。また、図６Ｄよりも、図６Ｅの方が高い周波数であることが分かる。また、図６Ｄ〜図６Ｅの各グラフにおけるピークが立っている周波数は、図６Ａ〜図６Ｃの各グラフにおけるピークが立っている周波数よりも顕著に高いことが分かる。そのため、使用固有値数ｓを１として、求められた係数αで特定される修正された自己回帰モデルの周波数特性から、軸受１０６の転動体又は外輪の疵の有無を判断することが可能であることが分かる。
【００７４】
以上より、自己相関行列Ｒの１個の固有値に基づいて求められた係数αで特定される修正された自己回帰モデルの周波数特性に基づいて、軸受１０６の状態診断を行うことができるという知見を得ることができた。
【００７５】
しかし、図６Ｂと図６Ｃとのグラフに示されるように、軸受１０６の内輪に疵がある場合と、軸受１０６の保持器に疵がある場合と、では、似通った値の周波数において周波数特性のピークが立っている。そのため、使用固有値数ｓを１として、求められた係数αで特定される修正された自己回帰モデルの周波数特性から、軸受１０６に内輪の又は保持器の疵の有無を判断することは可能であるが、内輪の疵であるか、保持器の疵であるか、両方の疵であるか、を判断することは困難である。
【００７６】
そこで、本実験では、更に、使用固有値数ｓの値を５に増やして、自己相関行列Ｒの固有値のうち、最大のものから５つを用いて係数αを求め、求めた係数αで特定される修正された自己回帰モデルの波形、周波数特性を求めた。
【００７７】
図７Ａ〜図７Ｅは、使用固有値数ｓを５として求められた係数αで特定される修正された自己回帰モデルの波形を示す図である。
【００７８】
図７Ａのグラフは、使用固有値数ｓを５として、軸受１０６に疵がない正常な軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフである。図７Ａのグラフでは、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図７Ｂのグラフは、使用固有値数ｓを５として、軸受１０６に内輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフである。図７Ｂのグラフでは、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図７Ｃのグラフは、使用固有値数ｓを５として、軸受１０６に保持器に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフである。図７Ｃのグラフでは、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図７Ｄのグラフは、使用固有値数ｓを５として、軸受１０６に転動体に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフである。図７Ｄのグラフでは、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図７Ｅのグラフは、使用固有値数ｓを５として、軸受１０６に外輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフである。図７Ｅのグラフでは、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。
【００７９】
図８Ａ〜図８Ｅは、使用固有値数ｓを５として求められた係数αで特定される修正された自己回帰モデルの周波数特性を示す図である。図８Ａ〜図８Ｅの各グラフは、図７Ａ〜図７Ｅの各波形についての周波数特性を示すグラフである。
【００８０】
図８Ａのグラフは、使用固有値数ｓを５として、軸受１０６に疵がない正常な軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフである。図８Ａのグラフでは、横軸が周波数、縦軸が信号強度を示す。図８Ｂのグラフは、使用固有値数ｓを５として、軸受１０６に内輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフである。図８Ｂのグラフでは、横軸が周波数、縦軸が信号強度を示す。図８Ｃのグラフは、使用固有値数ｓを５として、軸受１０６に保持器に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフである。図８Ｃのグラフでは、横軸が周波数、縦軸が信号強度を示す。図８Ｄのグラフは、使用固有値数ｓを５として、軸受１０６に転動体に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフである。図８Ｄのグラフでは、横軸が周波数、縦軸が信号強度を示す。図８Ｅのグラフは、使用固有値数ｓを５として、軸受１０６に外輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式１８で求められた周波数特性を示すグラフである。図８Ｅのグラフでは、横軸が周波数、縦軸が信号強度を示す。
【００８１】
図８Ａのグラフを見ると、周波数８２１Ｈｚと１０４７Ｈｚの部分にピークが立っているのが分かる。また、図８Ｂのグラフを見ると、周波数１６４８Ｈｚ、２０２１Ｈｚ、６４７４Ｈｚの部分にピークが立っているのが分かる。また、図８Ｃのグラフを見ると、周波数１６４６Ｈｚ、３２９１Ｈｚ、１７４６Ｈｚの部分にピークが立っているのが分かる。図８Ｄのグラフを見ると、周波数１９８５３Ｈｚの部分にピークが立っているのが分かる。また、図８Ｅのグラフを見ると、周波数２２７８５Ｈｚ、２３０２０Ｈｚの部分にピークが立っているのが分かる。
【００８２】
図８Ａ〜図８Ｅの各グラフにおけるピークと、図６Ａ〜図６Ｅの各グラフにおけるピークと、を見比べてみると、図６Ａ〜図６Ｅのグラフでピークが立っている周波数近辺においては、図８Ａ〜図８Ｅのグラフでもピークが立っているのが分かる。図８Ａ〜図８Ｅのグラフでは、更に、図６Ａ〜図６Ｅのグラフで見られなかったピークが見られるようになっている。
【００８３】
このように、係数αを求めるのに利用される自己相関行列Ｒの固有値を増やすことで、修正された自己回帰モデルによる予測値において、図６Ａ〜図６Ｅでピークが立っていた周波数以外の周波数の成分の信号強度が増加したことが分かる。即ち、係数αを求めるのに利用される自己相関行列Ｒの固有値は、求められた係数αで特定される修正された自己回帰モデルによる計測データの予測値の時間波形に含まれる各周波数の成分の強度と相関があることが見出された。
【００８４】
また、図８Ｂ、図８Ｃのグラフを見てみると、図６Ｂ、図６Ｃのグラフと同様に、１６４０Ｈｚ付近にピークが立っているのが分かる。しかし、図８Ｂでは、更に、２０２１Ｈｚ、６４７４Ｈｚにおいてもピークが立っており、図８Ｃでは、更に、３２９１Ｈｚ、１７４６Ｈｚにおいてもピークが立っていることが分かる。このように、図８Ｂ〜図８Ｃでは、内輪に疵のある軸受１０６と保持器に疵のある軸受１０６とで、ピークが立っている周波数に顕著な違いが見て取れる。即ち、使用固有値数ｓを５にすることで、軸受の疵が、内輪の疵であるか、保持器の疵であるか、両方の疵であるか、を見分けることができる。
【００８５】
このように本発明者らは、使用固有値数ｓが１の場合には見分けることが困難であった疵の存在する部位を、使用固有値数ｓを増やすことで、見分けることができるようになるという知見を得た。一方で、本発明者らは、使用固有値数ｓを増やし過ぎると修正された自己回帰モデルが、実際の計測データに近づいていくことによって、軸受１０６の状態診断に有用でないノイズ成分等についても信号強度が増加されてしまい、式１８で得られる周波数特性にもノイズ成分に起因するピークが生じるという知見を得た。即ち、自己相関行列Ｒの固有値のうち、相対的に値の高い一部の固有値を用いて係数αを求めることで、計測データに含まれる成分のうち、物体の状態診断に有用な成分を近似した修正された自己回帰モデルを求めることができると考えられる。
【００８６】
以上より、係数αを求めるのに利用される自己相関行列Ｒの固有値の数である使用固有値ｓを調整することで、より精度よく、軸受１０６の状態診断を行うことができるという知見を得ることができた。
【００８７】
（実験２）
また、発明者らは、実験１と同様の実験状況において、以下のような実験を行った。軸受１０６が疵のない正常な状態である場合、内輪に疵がある状態である場合、保持器に疵がある状態である場合、転動体に疵がある状態である場合、外輪に疵がある状態である場合それぞれについて、以下の作業をｑ回行った。ｑは５０としたが、３０や１００等の他の値としてもよい。即ち、軸受１０６について、連続したＭ個の計測データを取得し、取得した連続したＭ個の計測データ毎に、式１１を用いて行列Σｓと行列Ｕｓとを求め、式５と式１３を用いて修正された自己回帰モデルの係数αを求めた。使用固有値数ｓは、３とした。
【００８８】
これにより、軸受１０６が正常な状態である場合、内輪に疵がある状態である場合、保持器に疵がある状態である場合、転動体に疵がある状態である場合、外輪に疵がある状態である場合の５つの場合それぞれについて、ｑ個の修正された自己回帰モデルの係数αを取得した。
【００８９】
軸受１０６が取り得る状態は、疵のない正常な状態、内輪に疵がある状態、保持器に疵がある状態、転動体に疵がある状態、および外輪に疵がある状態の５つの状態であるとする。以下では、軸受１０６が取り得る状態の数をｐとする。ｐは、これらの５つの状態に対応して５とした。
【００９０】
図９は、軸受１０６が正常な状態である場合、内輪に疵がある状態である場合、保持器に疵がある状態である場合、転動体に疵がある状態である場合、軸受１０６の外輪に疵がある場合それぞれについての、上記手順で求められた係数αのパターンを示す図である。係数αのパターンとは、係数αの成分を並べた際の形状を示す情報である。ｍの値は、５００とした。そのため、係数αは、α１〜α５００の５００個の成分から成る。図９に示されるパターンそれぞれは、αｉ（１＜＝ｉ＜＝５００）をインデックスｉ順に並べたパターンである。このインデックスｉは、式１において係数の成分がかかる計測データｙの計測時点を示すインデックスとみなすことができる。即ち、図９に示されるパターンそれぞれは、係数αの各成分を、式１において各成分がかかる計測データｙの計測時点に基づいて並べたパターンの一例である。
【００９１】
係数αのパターンには、例えば、係数αの成分をインデックス順に並べたベクトル、係数αの成分をインデックス順に予め定められた数（例えば、１等）ずつとばしに並べたベクトル、係数αの成分をインデックスと逆順に並べたベクトル等がある。また、係数αのパターンには、例えば、各成分に同じ値が加算された係数αの成分をインデックス順に並べたベクトル、各成分に同じ値が加算された係数αの成分をインデックス順に予め定められた数（例えば、１等）ずつとばしに並べたベクトル、各成分に同じ値が加算された係数αの成分をインデックスと逆順に並べたベクトル等がある。また、係数αのパターンには、例えば、各成分に同じ値が掛けられた係数αの成分をインデックス順に並べたベクトル、各成分に同じ値が掛けられた係数αの成分をインデックス順に予め定められた数（例えば、１等）ずつとばしに並べたベクトル、各成分に同じ値が掛けられた係数αの成分をインデックスと逆順に並べたベクトル等がある。また、係数αの成分を並べた際の形状を示す情報である係数αのパターンは、ベクトルでなく、係数αの成分を並べたグラフの画像であるとしてもよい。
【００９２】
図９に示されるパターンを見てみると、以下のようなことが見てとれる。軸受１０６の外輪に疵がある場合、ｑ個の係数αのパターンそれぞれは、上下対称な２つの波線で囲まれた領域を示すようなパターンとなっている。また、軸受１０６の転動体に疵がある場合、ｑ個の係数αのパターンそれぞれは、複数の波線が重なり合っているようなパターンとなっている。また、軸受１０６の内輪に疵がある場合、ｑ個の係数αのパターンそれぞれは、４つのピークのある山形のパターンが連続しているパターンとなっている。また、軸受１０６に疵がない場合、ｑ個の係数αのパターンそれぞれは、徐々に減衰していく波型のパターンとなっている。また、軸受１０６の保持器に疵がある場合、ｑ個の係数αのパターンそれぞれは、２つのピークのある山形のパターンが連続しているパターンとなっている。
【００９３】
このように、軸受１０６が正常な場合、軸受１０６の内輪に疵がある場合、軸受１０６の保持器に疵がある場合、軸受１０６の転動体に疵がある場合、軸受１０６の外輪に疵がある場合それぞれについて、修正された自己回帰モデルの係数αが、特徴的なパターンを示すことが見て取れる。
【００９４】
そこで、発明者らは、修正された自己回帰モデルの係数αのパターンが、どのようなパターンを示すのかを特定することで、軸受１０６の診断を行うことができるという知見を得た。
【００９５】
また、本実験では、以下のようにして、軸受１０６が正常な場合、軸受１０６の内輪に疵がある場合、軸受１０６の保持器に疵がある場合、軸受１０６の転動体に疵がある場合、軸受１０６の外輪に疵がある場合それぞれについて、修正された自己回帰モデルの係数αの代表的なパターンを学習により特定した。
【００９６】
まず、各列が図９に示される係数αそれぞれを示すｍ×ｄのサイズの行列Ｙを、以下の式１９のように生成する。ｄは、ｐ×ｑであり、２５０となる。
【００９７】
【数１９】
2020145215.tif 000021
【００９８】
式１９のαｋ，（ｉ，ｊ）は、ｉに対応する状態の軸受１０６から計測されたｑ個の計測データのうちｊ番目の計測データから決定された係数αの第ｋ成分αｋを示す。本実験では、外輪に疵がある状態の軸受１０６から計測された計測データから決定された係数αの成分は、式１９におけるαｋ，（０，ｊ）（１＜＝ｋ＜＝ｍ、１＜＝ｊ＜＝ｑ）である。また、転動体に疵がある状態の軸受１０６から計測された計測データから決定された係数αの成分は、式１９におけるαｋ，（１，ｊ）（１＜＝ｋ＜＝ｍ、１＜＝ｊ＜＝ｑ）である。また、内輪に疵がある状態の軸受１０６から計測された計測データから決定された係数αの成分は、式１９におけるαｋ，（２，ｊ）（１＜＝ｋ＜＝ｍ、１＜＝ｊ＜＝ｑ）である。また、疵がない正常な状態の軸受１０６から計測された計測データから決定された係数αの成分は、式１９におけるαｋ，（３，ｊ）（１＜＝ｋ＜＝ｍ、１＜＝ｊ＜＝ｑ）である。また、保持器に疵がある状態の軸受１０６から計測された計測データから決定された係数αの成分は、式１９におけるαｋ，（４，ｊ）（１＜＝ｋ＜＝ｍ、１＜＝ｊ＜＝ｑ）である。ここでは、ｐは５である。式１９では、表記の都合で、αｋ，（４，ｊ）をαｋ，（ｐ−１，ｊ）とする。また、式１９では、αｋ，（２，ｊ）およびαｋ，（３，ｊ）の表記を省略する。
【００９９】
本実験では、この行列Ｙを非負値行列因子分解することで、軸受１０６の状態それぞれについて、修正された自己回帰モデルの係数αの代表的なパターンを特定する。行列Ｙを非負値行列因子分解するためには、行列Ｙの各列の係数αのパターンを変えずに、行列Ｙ自体の各成分を、非負値にするよう調整する必要がある。そこで、例えば、行列Ｙの０未満の成分のうち、最小のものの絶対値を、行列Ｙの各成分に加算することで、行列Ｙの各成分が非負値となるようにする。以下では、各成分が非負値である行列を、非負値行列とする。
【０１００】
以下の式２０に示すように、この行列Ｙを、非負値行列である行列Ａと行列Ｐとに非負値行列因子分解する。行列Ａは、ｍ×ｐのサイズの行列である。行列Ｐは、ｐ×ｄのサイズの行列である。行列Ａの列数は、非負値行列因子分解における基底ベクトルの数であり、本実験では、軸受１０６が取り得る状態の数としてｐ（５）とする。これにより、行列Ａの各列には、行列Ｙから導出されたｐ（５）個の係数αのパターン（基底ベクトルともいう）が格納されることとなる。
【０１０１】
【数２０】
2020145215.tif 000022
【０１０２】
以下の式２１のようにコスト関数Ｊを定義する。
【０１０３】
【数２１】
2020145215.tif 000023
【０１０４】
式２１の右辺第１項は、行列Ａと行列Ｐとの積と、行列Ｙと、が、これらの残差の二乗和を最小にするという意味で近似的に一致するための条件を示す。式２１の右辺第１項の右下のＦｒ記号は、フロベニウスノルムであることを示す記号である。即ち、式２１の右辺第１項内の｜｜ＡＰ−Ｙ｜｜は、フロベニウスノルムであり、行列（ＡＰ−Ｙ）の各成分の絶対値の２乗の合計の平方根となる。式２１の右辺第２項は、行列Ｐのスパース性を要求する条件を示す。この条件を付加することによって、非負値行列因子分解が過学習となることを防ぐことができる。
【０１０５】
式２１の右辺第２項におけるインデックスｉは、行列の行を示すインデックスである。また、式２１の右辺第２項におけるインデックスｋは、行列の列を示すインデックスである。式２１の右辺第３項は、行列Ａの異なる列に格納されるパターン同士の直交性を近似的に要求する条件を示す。式２１の右辺第３項によって、行列Ａのｐ（５）個の列に格納されるパターンが可能な限り線形従属となることを防ぐ効果が得られる。行列Ａの列に格納されるパターンが線形従属となれば、あるパターンを行列Ａの列に格納されるパターンの線形結合で近似するときの係数（重みともいう）が一通りではなくなるため、パターンの同定が難しくなる。また、式２１の右辺第３項における行列ＡＴは、行列Ａの転置行列であり、行列Ｉは、単位行列である。また、式２１の右辺第３項におけるインデックスｋ、ｋ−（式中では、上にバーが引かれたｋとして表現）は、それぞれ行列の行と、列と、を示すインデックスである。
【０１０６】
ここで、コスト関数Ｊとして、式２１の右辺第１項のみを用いても、非負値行列因子分解は機能するので、コスト関数Ｊとして、式２１の右辺第１項のみを用いてもよい。また、上述の過学習や線形従属といった問題が生じる場合に、コスト関数Ｊとして、式２１の右辺第１項に式２１の右辺第２項、第３項を追加した式を用いることとしてもよい。そして、右辺第２項、第３項を追加する場合、軸受１０６の各状態と、図１４Ａ〜図１４Ｂ等で後述する方法で決定される学習パターンと、を良好に対応付けることが出来るように、例えば、ケーススタディにより係数βと係数γとを予め決定しておく。本実施形態では、βを０．０３５、γを０とした。
【０１０７】
このコスト関数Ｊを最小化するように、行列Ｙを、行列Ａと行列Ｐと分解する。即ち、以下の式２２で定義されるコスト関数Ｊの最小化問題を解くこととなる。
【０１０８】
【数２２】
2020145215.tif 000024
【０１０９】
式２２における第１の制約条件は、行列Ａと行列Ｐとが非負値行列であることを示す条件である。式２２における第２の制約条件は、行列Ａを基準化（ｎｏｒｍａｌｉｚｅ）するための条件である。また、式２２の第２の制約条件におけるインデックスｋは、行列の行と、列と、を示すインデックスである。
本実験では、コスト関数Ｊを最小化するために、以下のような処理を行う。即ち、以下の式２３と式２４とを用いて行列Ａと行列Ｐとの更新を繰り返すことで、コスト関数Ｊを最小化する。また、コスト関数Ｊを最小化するために、以下に記載する方法以外の方法（例えば、最急降下法、ニュートン法、確率的勾配降下法等の最適化手法）を用いてもよい。
【０１１０】
【数２３】
2020145215.tif 000025
【０１１１】
【数２４】
2020145215.tif 000026
【０１１２】
式２３のηＰは、行列Ｐの更新に用いられる緩和係数である。式２４のηＡは、行列Ａの更新に用いられる緩和係数である。
式２３におけるコスト関数Ｊの行列Ｐによる偏微分は、以下の式２５の通りとなる。
【０１１３】
【数２５】
2020145215.tif 000027
【０１１４】
式２５の行列Ｅ１は、サイズがｐ×ｄの全ての成分が１の行列である。
式２４におけるコスト関数Ｊの行列Ａによる偏微分は、以下の式２６の通りとなる。
【０１１５】
【数２６】
2020145215.tif 000028
【０１１６】
式２６の行列Ｅ２は、サイズがｐ×ｐの全ての成分が１の行列である。
式２５を用いて、式２３を書き換えると以下の式２７が求まる。
【０１１７】
【数２７】
2020145215.tif 000029
【０１１８】
式２７のインデックスｉは、行列の行を示すインデックスである。式２７のインデックスｋは、行列の列を示すインデックスである。
【０１１９】
行列Ｐは、非負値であるため、式２７の値は、非負値である必要がある。式２７の右辺の第３項は、必ず非負値となる項である。即ち、式２７の右辺から第３項を除いた式の値が０であるならば、必ず式２７は、非負値となる。そのため、以下の式２８に示すように、式２７の右辺から第３項を除いた式が０であるという条件を満たすように、行列Ｐの更新を行うこととする。
【０１２０】
【数２８】
2020145215.tif 000030
【０１２１】
式２８のインデックスｉは、行列の行を示すインデックスである。式２８のインデックスｋは、行列の列を示すインデックスである。
式２８から、式２３における緩和係数ηＰの値が、以下の式２９のように求まる。
【０１２２】
【数２９】
2020145215.tif 000031
【０１２３】
式２９のインデックスｉは、行列の行を示すインデックスである。式２９のインデックスｋは、行列の列を示すインデックスである。
式２９を用いて、式２３を書き換えると、以下の式３０のようになる。
【０１２４】
【数３０】
2020145215.tif 000032
【０１２５】
式３０のインデックスｉは、行列の行を示すインデックスである。式３０のインデックスｋは、行列の列を示すインデックスである。この式３０を用いて、行列Ｐの各成分を更新することで、行列Ｐを更新していくこととなる。
また、式２６を用いて、式２４を書き換えると以下の式３１が求まる。
【０１２６】
【数３１】
2020145215.tif 000033
【０１２７】
式３１のインデックスｋは、行列の行を示すインデックスである。式２７のインデックスｊは、行列の列を示すインデックスである。
行列Ａは、非負値であるため、式３１の値は、非負値である必要がある。式３１の右辺の第３項は、必ず非負値となる項である。即ち、式３１の右辺から第３項を除いた式の値が０であるならば、必ず式２７は、非負値となる。そのため、以下の式３２に示すように、式３１の右辺から第３項を除いた式が０であるという条件を満たすように、行列Ａの更新を行うこととする。
【０１２８】
【数３２】
2020145215.tif 000034
【０１２９】
式３２のインデックスｋは、行列の行を示すインデックスである。式３２のインデックスｊは、行列の列を示すインデックスである。
式３２から、式２４における緩和係数ηＡの値が、以下の式３３のように求まる。
【０１３０】
【数３３】
2020145215.tif 000035
【０１３１】
式３３のインデックスｋは、行列の行を示すインデックスである。式３３のインデックスｊは、行列の列を示すインデックスである。
式３３を用いて、式２４を書き換えると、以下の式３４のようになる。
【０１３２】
【数３４】
2020145215.tif 000036
【０１３３】
式３４のインデックスｋは、行列の行を示すインデックスである。式３４のインデックスｊは、行列の列を示すインデックスである。この式３４を用いて、行列Ａの各成分を更新することで、行列Ａを更新していくこととなる。
また、式２２における第２の制約条件を満たすようにするため、以下の式３５に示すように、行列Ａを更新する。
【０１３４】
【数３５】
2020145215.tif 000037
【０１３５】
式３５のｋは、行列の行を示すインデックスである。式３５のｊは、行列の行又は列を示すインデックスである。
コスト関数Ｊの値が、収束するまで、式３０を用いて行列Ｐを更新し、式３４、式３５を用いて行列Ａを更新することを繰り返す。例えば、各更新段階において、更新後の行列Ｐ及び行列Ａから、それぞれ更新前の行列Ｐ及び行列Ａを引いた差である行列△Ｐ及び行列△Ａを求める。そして、行列△Ｐ及び行列△Ａそれぞれのフロベニウスノルムを計算する。各更新段階における行列△Ｐ及び行列△Ａのフロベニウスノルムの、それぞれ最初の更新段階における行列△Ｐ及び行列△Ａのフロベニウスノルムに対する比が、予め定められた閾値以下となる更新段階をもって、コスト関数Ｊの値が収束したとする。そして、その際の行列Ａと行列Ｐとをコスト関数Ｊを最小化させる行列とする。
【０１３６】
本実験では、以上のように、コスト関数Ｊを最小化させる行列Ａと行列Ｐとを求めることとした。ただし、他の例として、行列Ｐと行列Ａとの更新前後で、コスト関数Ｊの値の変動が予め定められた閾値以下となることが、予め定められた回数連続したら、コスト関数Ｊの値が収束したとすることとしてもよい。そして、その際の行列Ａと行列Ｐとをコスト関数Ｊを最小化させる行列とすることとしてもよい。
【０１３７】
コスト関数Ｊが収束した際の行列Ａの各列が示すパターンを確認した。図１０に、コスト関数Ｊが収束した際の行列Ａの各列が示すパターンの一例を示す。行列Ａ列には、列毎に、ｍ（５００）個の値が含まれ、係数αのパターンを示す。
図１０の一番上のグラフは、行列Ａの１列目のｍ個の値（Ａｉ，１（１＜＝ｉ＜＝ｍ））を、行を示すインデックス順に並べたパターンを示す。図１０の上から２つ目のグラフは、行列Ａの２列目のｍ個の値（Ａｉ，２（１＜＝ｉ＜＝ｍ））を、行を示すインデックス（ｉ）順に並べたパターンを示す。図１０の上から３つ目のグラフは、行列Ａの３列目のｍ個の値（Ａｉ，３（１＜＝ｉ＜＝ｍ））を、行を示すインデックス順に並べたパターンを示す。図１０の下から２つ目のグラフは、行列Ａの４列目のｍ個の値（Ａｉ，４（１＜＝ｉ＜＝ｍ））を、行を示すインデックス順に並べたパターンを示す。図１０の一番下のグラフは、行列Ａの５列目（ｐ列目）のｍ個の値（Ａｉ，５（１＜＝ｉ＜＝ｍ））を、行を示すインデックス順に並べたパターンを示す。
【０１３８】
また、コスト関数Ｊが収束した際の行列Ｐの各成分の値を確認した。式２０から得られる以下の式３６に示されるように、行列Ｐの各列（第ｊ列Ｐｋ，ｊ（１＜＝ｋ＜＝ｐ））は、行列Ｙの各列（第ｊ列Ｙｉ，ｊ（１＜＝ｉ＜＝ｍ））に示されるパターンを行列Ａの各列に示されるパターン（基底ベクトル）の線形結合で近似するときの係数（重み）を示している。
【０１３９】
【数３６】
2020145215.tif 000038
【０１４０】
図１１に、コスト関数Ｊが収束した際の行列Ｐの各行の値を示すグラフを示す。図１１のグラフの横軸は、行列Ｐの列のインデックスを示す。図１１のグラフの縦軸は、行列Ｐの各成分の値を示す。
【０１４１】
図１１のライン１１０１は、行列Ｐの１行目のｄ（２５０）個の値を結んだラインである。ライン１１０１は、横軸が０以上、５０未満の範囲（１列目からｑ（５０）列目までの範囲）では、約１２となり、それ以外の範囲では、１２と比べて非常に微小な値（約０）となっている。
【０１４２】
図１１のライン１１０２は、行列Ｐの２行目のｄ個の値を結んだラインである。ライン１１０２は、横軸が５０以上、１００未満の範囲（ｑ＋１列目から２ｑ列目までの範囲）では、約１２となり、それ以外の範囲では、１２と比べて非常に微小な値（約０）となっている。
【０１４３】
図１１のライン１１０３は、行列Ｐの４行目のｄ個の値を結んだラインである。ライン１１０３は、横軸が１００以上、１５０未満の範囲（２ｑ＋１列目から３ｑ列目までの範囲）では、約１２となり、それ以外の範囲では、１２と比べて非常に微小な値（約０）となっている。
【０１４４】
図１１のライン１１０４は、行列Ｐの５行目のｄ個の値を結んだラインである。ライン１１０４は、横軸が１５０以上、２００未満の範囲（３ｑ＋１列目から４ｑ列目までの範囲）では、値が約１２となり、それ以外の範囲では、１２と比べて非常に微小な値（約０）となっている。
【０１４５】
図１１のライン１１０５は、行列Ｐの３行目のｄ個の値を結んだラインである。ライン１１０５は、横軸が２００以上、２５０未満の範囲（４ｑ＋１列目から５ｑ列目までの範囲）では、約１２となり、それ以外の範囲では、１２と比べて非常に微小な値（約０）となっている。
【０１４６】
即ち、行列Ｐの各列の値は、１列目からｑ（５０）列目までの範囲では、１行目の値が約１２となり、２〜５行目の値が約０となる。また、行列Ｐの各列の値は、ｑ＋１列目から２ｑ列目までの範囲では、２行目の値が約１２となり、１行目、３〜５行目の値が約０となる。また、行列Ｐの各列の値は、２ｑ＋１列目から３ｑ列目までの範囲では、４行目の値が約１２となり、１〜３行目、５行目の値が約０となる。また、行列Ｐの各列の値は、３ｑ＋１列目から４ｑ列目までの範囲では、５行目の値が約１２となり、１〜４行目の値が約０となる。また、行列Ｐの各列の値は、４ｑ＋１列目から５ｑ列目までの範囲では、３行目の値が約１２となり、１、２、４、５行目の値が約０となる。
そのため、行列ＡＰの１列目からｑ列目までの各列のデータが示すパターンは、図１０の一番上のグラフに示されるパターンに約１２の重み、図１０のそれ以外のグラフに示されるパターンに約０の重み、で図１０の各パターンが重ね合わされたパターンであり、図１０の一番上のグラフに示されるパターンに類似するパターンとなる。また、行列ＡＰのｑ＋１列目から２ｑ列目までの各列のデータが示すパターンは、図１０の上から２つ目のグラフに示されるパターンに約１２の重み、図１０のそれ以外のグラフに示されるパターンに約０の重み、で図１０の各パターンが重ね合わされたパターンであり、図１０の上から２つ目のグラフに示されるパターンに類似するパターンとなる。また、行列ＡＰの２ｑ＋１列目から３ｑ列目までの各列のデータが示すパターンは、図１０の下から２つ目のグラフに示されるパターンに約１２の重み、図１０のそれ以外のグラフに示されるパターンに約０の重み、で図１０の各パターンが重ね合わされたパターンであり、図１０の下から２つ目のグラフに示されるパターンに類似するパターンとなる。また、行列ＡＰの３ｑ＋１列目から４ｑ列目までの各列のデータが示すパターンは、図１０の一番下のグラフに示されるパターンに約１２の重み、図１０のそれ以外のグラフに示されるパターンに約０の重み、で図１０の各パターンが重ね合わされたパターンであり、図１０の一番下のグラフに示されるパターンに類似するパターンとなる。また、行列ＡＰの４ｑ＋１列目から５ｑ列目までの各列のデータが示すパターンは、図１０の上から３番目のグラフに示されるパターンに約１２の重み、図１０のそれ以外のグラフに示されるパターンに約０の重み、で図１０の各パターンが重ね合わされたパターンであり、図１０の上から３つ目のグラフに示されるパターンに類似するパターンとなる。
【０１４７】
行列Ｙの１列目からｑ列目までの各列のデータが示すパターンは、図９に示す外輪に疵がある軸受１０６に対応するパターンである。また、行列Ｙのｑ＋１列目から２ｑ列目までの各列のデータが示すパターンは、図９に示す転動体に疵がある軸受１０６に対応するパターンである。また、行列Ｙの２ｑ＋１列目から３ｑ列目までの各列のデータが示すパターンは、図９に示す内輪に疵がある軸受１０６に対応するパターンである。また、行列Ｙの３ｑ＋１列目から４ｑ列目までの各列のデータが示すパターンは、図９に示す疵のない正常な軸受１０６に対応するパターンである。また、行列Ｙの４ｑ＋１列目から５ｑ列目までの各列のデータが示すパターンは、図９に示す保持器に疵がある軸受１０６に対応するパターンである。
【０１４８】
行列Ｙと行列Ａとを見比べてみると、双方の各列のデータが示すパターンは、類似するパターンとなっていることが分かる。即ち、行列Ａの各列のデータが示すパターンは、それぞれ、各状態の軸受１０６に対応する係数αのパターンとなっていることが分かる。行列Ａの１列目のデータは、外輪に疵がある状態の軸受１０６に対応する係数αの代表的なパターンを示す。行列Ａの２列目のデータは、転動体に疵がある状態の軸受１０６に対応する係数αの代表的なパターンを示す。行列Ａの３列目のデータは、保持器に疵がある状態の軸受１０６に対応する係数αの代表的なパターンを示す。行列Ａの４列目のデータは、内輪に疵のある状態の軸受１０６に対応する係数αの代表的なパターンを示す。行列Ａの５列目のデータは、疵のない正常な状態の軸受１０６に対応する係数αの代表的なパターンを示す。
【０１４９】
以上より、行列Ｙを生成し、生成した行列Ｙを、式２２を満たすように、非負値行列因子分解により学習することで、行列Ａの各列（基底ベクトル）として、各状態の軸受１０６に対応する係数αのパターンが特定されたことが確認された。
【０１５０】
（実験３）
また、発明者は、軸受１０６の状態を診断する方法の着想を得た。この着想の概要について説明する。状態が既知な軸受１０６についての計測データから、予めその状態に対応する修正された自己回帰モデルの係数αのパターンを示す非負値の行列を、教師基底行列として特定する。そして、状態が不知な軸受１０６についての計測データから修正された自己回帰モデルの係数αを予め定められた個数だけ求めて、求めたその個数の係数αそれぞれのパターンを示す非負値のデータを各列に格納する行列を、教師基底行列と、教師基底行列の重み行列と、この非負値のデータを各列に格納する行列から抽出される基底を示す抽出基底行列と、抽出基底行列の重み行列と、に分解する。そして、教師基底行列の重み行列と、抽出基底行列の重み行列と、を比較することで、状態が不知な軸受１０６の状態が、教師基底行列に対応する状態であるか否かを診断する。以上が、発明者が得た着想の概要である。
【０１５１】
そこで、発明者らは、この着想で得られた方法の実効性を検証するため、実験１と同様の実験状況において、以下のような実験を行った。本実験では、使用固有値数ｓは、３とした。また、ｍは、５００とした。
【０１５２】
実験２で説明した方法で、疵のない正常な状態の軸受１０６に対応する係数αのパターンを示すデータ（式２０の行列Ａの正常な軸受１０６に対応する列のデータ）を、教師基底に属する基底ベクトルとして取得した。取得した正常な状態の軸受１０６に対応する教師基底に属するパターンは、図１２に示すようなパターンとなった。図１２のグラフは、非負値となるように修正された正常な状態の軸受１０６に対応する係数αを示す。図１２のグラフの横軸は、係数αのインデックスを示し、縦軸は、対応するデータの値を示す。
【０１５３】
そして、発明者らは、軸受１０６が疵のない正常な状態である場合、図１３Ａ〜図１３Ｃに示す各状態にある場合それぞれと、について、以下の作業をｑ’回行った。図１３Ａ〜図１３Ｃに示す各状態とは、図１３Ａに示す外輪に幅３ｍｍ・深さ約０．１ｍｍの疵がある状態と、図１３Ｂに示す外輪に幅０．５ｍｍ・深さ０．１ｍｍ以下の疵がある状態と、図１３Ｃに示す外輪に予め定められた大きさ（径が１ｍｍ等）以下の微小な疵がある状態とのそれぞれである。ｑ’は、２００としたが、１００等の他の値としてもよい。各場合について、軸受１０６についての連続したＭ個の計測データを取得し、取得した連続したＭ個の計測データに基づいて、式１１を用いて行列Σｓと行列Ｕｓとを求め、式５と式１３を用いて修正された自己回帰モデルの係数αを求める作業をｑ’回行った。
【０１５４】
ｑ’が大きいほど診断結果の信頼性が向上するため、ｑ’は、大きい方が望ましい。ｑ’の適正な値は、診断の対象、診断結果に要求される信頼性の程度、計測データの種類等に応じ、実験等の手段を駆使して定めればよい。本実験３では、ｑ’は、１０以上とするのが望ましく、５０以上とするのがより望ましく、１００以上とするのがさらに望ましい。
【０１５５】
これにより、疵のない正常な状態の軸受１０６と、図１３Ａ〜図１３Ｃに示す各状態にある軸受１０６それぞれと、について、ｑ’個の修正された自己回帰モデルの係数αを取得した。
【０１５６】
そして、軸受１０６が正常な場合、図１３Ａ〜図１３Ｃに示す各状態にある場合それぞれについて、以下の式３７で示すようなデータ行列Ｙｄを生成した。
【０１５７】
【数３７】
2020145215.tif 000039
【０１５８】
式３７の行列の各列のデータは、取得されたｑ’個の修正された自己回帰モデルの係数αそれぞれを示すデータである。即ち、αｉ，ｊは、ｑ’個の修正された自己回帰モデルの係数αのうちのｊ個目の係数αにおけるｉ番目の要素を示す。
本実験では、このデータ行列Ｙｄを、非負値の複数の行列に分解する。データ行列Ｙｄを非負値の複数の行列に分解するためには、行列Ｙの各列の係数αのパターンを変えずに、行列Ｙ自体の各成分を、非負値にするよう調整する必要がある。そこで、例えば、行列Ｙの０未満の成分のうち、最小のものの絶対値を、データ行列Ｙｄの各成分に加算することで、データ行列Ｙｄの各成分が非負値となるようにする。データ行列Ｙｄの各列には、修正された自己回帰モデルにおける係数αのパターンを示す非負値のデータが格納されていることとなる。
【０１５９】
以下の式３８に示すように、非負値の行列に修正したデータ行列Ｙｄを、複数の非負値行列に分解する。より具体的には、データ行列Ｙｄを、教師基底行列Ｆと、教師基底に対する重みを示す教師重み行列Ｇと、データ行列Ｙｄから抽出される基底である抽出基底を示す抽出基底行列Ｑと、抽出基底に対する重みを示す抽出基底重み行列Ｗと、に分解する。教師基底行列Ｆと、教師重み行列Ｇと、抽出基底行列Ｑと、抽出基底重み行列Ｗは、何れも非負値行列である。本実験では、教師基底は、疵のない正常な状態の軸受１０６に対応する係数αのパターンを示す基底ベクトルの集合である。抽出基底は、疵のない正常な状態の軸受１０６に対応する係数αのパターン以外の係数αのパターンを示す基底ベクトルの集合である。疵のない正常な状態の軸受１０６に対応する係数αのパターン以外の係数αのパターンは、例えば、状態が不知の軸受１０６に対応する係数αのパターンを含む。
【０１６０】
ＹｄをＦとＧとＱとＷとに分解する処理は、以下のような処理とみなすことができる。即ち、Ｙｄの各列に格納された係数αのパターンのデータを、教師基底に属するパターンと、それ以外のパターンと、の重み付けの足し合わせとして分解する処理とみなすことができる。
【０１６１】
Ｆは、ｍ×１のサイズの行列である。Ｇは、１×ｑ’のサイズの行列である。Ｑは、ｍ×（データ行列から抽出される基底ベクトルの数を示す数ｎ（本実験では、１とした））のサイズの行列である。Ｗは、ｎ×ｑ’のサイズの行列である。
Ｇにおける各列の要素の値は、データ行列Ｙｄのその列に格納された係数αのパターンを示すデータに含まれる教師基底の成分の重みを示す。Ｗにおけるｎ行目の各列の要素の値は、データ行列Ｙｄのその列に格納された係数αのパターンを示すデータに含まれるＱのその列に格納された抽出基底の成分の重みを示す。
【０１６２】
【数３８】
2020145215.tif 000040
【０１６３】
以下の式３９のようにコスト関数Ｊｄを定義する。
【０１６４】
【数３９】
2020145215.tif 000041
【０１６５】
式３９の右辺第１項のＦｒは、フロベニウスノルムであることを示す添え字である。即ち、式３９の右辺第１項は、フロベニウスノルムである。式３９の右辺第１項は、ＦとＧとの積とＱとＷとの積との和と、データ行列Ｙと、の差分の大きさを示す。式３９の右辺第２項は、教師重み行列のスパース制約を示す項である。式３９の右辺第２項のδ１は、予め定められた実数である。式３９の右辺第２項のｋ、ｔは、それぞれ行列の行、列を示すインデックスである。
【０１６６】
式３９の右辺第３項は、抽出重み行列のスパース制約を示す項である。式３９の右辺第３項のδ２は、予め定められた実数である。式３９の右辺第３項のｌ、ｔは、それぞれ行列の行、列を示すインデックスである。式３９の右辺第４項は、抽出基底の正規直交制約（即ち、抽出基底行列Ｑの各列に格納されているパターンの全体が近似的に正規直交系をなすとの制約）を示す項である。即ち、式３９の右辺第４項は、正規直交性の向上に関する項である。式３９の右辺第４項のε１は、予め定められた実数である。式３９の右辺第４項のｉ、ｉ’は、それぞれ行列の行、列を示すインデックスである。式３９の右辺第４項のＩは、単位行列である。式３９の右辺第５項は、教師基底と抽出基底との間の直交制約（即ち、教師基底行列Ｆの各列に格納されている全てのパターンと抽出基底行列Ｑの各列に格納されている全てのパターンとが近似的に直交するとの制約）を示す項である。式３９の右辺第５項のε２は、予め定められた実数である。式３９の右辺第４項のｉ、ｉ’は、それぞれ行列の行、列を示すインデックスである。
【０１６７】
本実験では、教師重み行列のスパース制約、抽出重み行列のスパース制約、抽出基底の正規直交制約、教師基底と抽出基底との間の直交制約の何れも付加しない。即ち、δ１、δ２、ε１、ε２それぞれは、０とする。そのため、コスト関数Ｊｄは、式３９の右辺第１項のみの式となる。
【０１６８】
Ｊｄの値を適正化するように、Ｇ、Ｑ、Ｗの値を更新することで、Ｙｄを、ＦとＧとＱとＷとに分解する。以下の式４０で定義されるコスト関数Ｊｄの最小化問題を解くこととなる。
【０１６９】
【数４０】
2020145215.tif 000042
【０１７０】
本実験では、コスト関数Ｊｄを最小化するために、以下のような処理を行う。即ち、ＧとＱとＷとの更新を繰り返すことで、コスト関数Ｊｄを最小化する。また、コスト関数Ｊｄを最小化するために、以下に記載する方法以外の方法（例えば、ニュートン法、確率的勾配降下法等の最適化手法）を用いてもよい。
Ｇの更新について説明する。本実験では、以下の式４１に示すようにＧの更新を行った。
【０１７１】
【数４１】
2020145215.tif 000043
【０１７２】
式４１のｋ、ｔは、それぞれ行列の行、列を示すインデックスである。ηＧｋｔは、Ｇのｋ、ｔ要素についての更新のステップサイズを示す非負値の値である。
コスト関数ＪｄをＧで偏微分することで、以下の式４２が得られる。
【０１７３】
【数４２】
2020145215.tif 000044
【０１７４】
式４２のＥ３は、サイズが１×ｑ’の全ての成分が１の行列である。
式４１に式４２を代入することで、更新後のＧを示す以下の式４３が得られる。
【０１７５】
【数４３】
2020145215.tif 000045
【０１７６】
行列Ｇは、非負値であるため、式４３の値は、非負値である必要がある。式４３の最右辺の末尾の項は、必ず非負値となる項である。即ち、式４３の最右辺から末尾の項を除いた式の値が０であるならば、必ず更新後のＧは、非負値となる。そのため、以下の式４４に示すように、式４３の最右辺から末尾の項を除いた式が０であるという条件を満たすように、行列Ｇの更新を行うこととする。
【０１７７】
【数４４】
2020145215.tif 000046
【０１７８】
式４４からηＧｋｔは、以下の式４５のようになる。
【０１７９】
【数４５】
2020145215.tif 000047
【０１８０】
式４１に式４２と式４５とを代入することで、更新後のＧを示す以下の式４６が得られる。
【０１８１】
【数４６】
2020145215.tif 000048
【０１８２】
即ち、Ｙｄと、Ｆと、処理時点における更新前のＧ、Ｑ、Ｗと、予め定められた係数δ１（本実験では、０）と、に基づいて、式４６を用いて、Ｇを更新した。
次に、Ｑの更新について説明する。本実験では、以下の式４７に示すようにＱを更新した。
【０１８３】
【数４７】
2020145215.tif 000049
【０１８４】
式４７のｗ、ｔは、それぞれ行列の行、列を示すインデックスである。ηＱｗｔは、Ｑのｗ、ｔ要素についての更新のステップサイズを示す非負値の値である。
コスト関数ＪｄをＱで偏微分することで、以下の式４８が得られる。
【０１８５】
【数４８】
2020145215.tif 000050
【０１８６】
式４８のＥ４は、サイズがｎ×ｎの全ての成分が１の行列である。式４８のＥ５は、サイズが１×１の全ての成分が１の行列である。
式４７に式４８を代入することで、更新後のＱを示す以下の式４９が得られる。
【０１８７】
【数４９】
2020145215.tif 000051
【０１８８】
行列Ｑは、非負値であるため、式４９の値は、非負値である必要がある。式４９の最右辺の末尾の項は、必ず非負値となる項である。即ち、式４９の最右辺から末尾の項を除いた式の値が０であるならば、必ず更新後のＱは、非負値となる。そのため、以下の式５０に示すように、式４９の最右辺から末尾の項を除いた式が０であるという条件を満たすように、行列Ｑを更新することとする。
【０１８９】
【数５０】
2020145215.tif 000052
【０１９０】
式５０からηＱｗｔは、以下の式５１のようになる。
【０１９１】
【数５１】
2020145215.tif 000053
【０１９２】
４７式に４８式と５１式とを代入することで、更新後のＧを示す以下の式５２が得られる。
【０１９３】
【数５２】
2020145215.tif 000054
【０１９４】
即ち、Ｙｄと、Ｆと、処理時点における更新前のＧ、Ｑ、Ｗと、予め定められた係数ε１、ε２（本実験では、共に０）と、に基づいて、式５２を用いて、Ｑを更新した。
次に、Ｗの更新について説明する。本実験では、以下の式５３に示すようにＷを更新した。
【０１９５】
【数５３】
2020145215.tif 000055
【０１９６】
式５３のｌ、ｔは、それぞれ行列の行、列を示すインデックスである。ηＷｌｔは、Ｗのｌ、ｔ要素についての更新のステップサイズを示す非負値の値である。
コスト関数ＪｄをＷで偏微分することで、以下の式５４が得られる。
【０１９７】
【数５４】
2020145215.tif 000056
【０１９８】
式５４のＥ６は、サイズが１×ｑ’の全ての成分が１の行列である。
式５３に式５４を代入することで、更新後のＷを示す以下の式５５が得られる。
【０１９９】
【数５５】
2020145215.tif 000057
【０２００】
行列Ｗは、非負値であるため、式５５の値は、非負値である必要がある。式５５の最右辺の末尾の項は、必ず非負値となる項である。即ち、式５５の最右辺から末尾の項を除いた式の値が０であるならば、必ず更新後のＷは、非負値となる。そのため、以下の式５６に示すように、式５５の最右辺から末尾の項を除いた式が０であるという条件を満たすように、行列Ｗの更新を行うこととする。
【０２０１】
【数５６】
2020145215.tif 000058
【０２０２】
式５６からηＷｌｔは、以下の式５７のようになる。
【０２０３】
【数５７】
2020145215.tif 000059
【０２０４】
５３式に５４式と５７式とを代入することで、更新後のＷを示す以下の式５８が得られる。
【０２０５】
【数５８】
2020145215.tif 000060
【０２０６】
即ち、Ｙｄと、Ｆと、処理時点における更新前のＧ、Ｑ、Ｗと、予め定められた係数δ２（本実験では、０）と、に基づいて、式５８を用いて、Ｇを更新した。
以上のように、本実験では、式４６を用いたＧの更新と、式５２を用いたＱの更新と、式５８を用いたＷの更新と、を繰り返す。各更新段階において、更新後のＧ、Ｑ、Ｗから、それぞれ更新前のＧ、Ｑ、Ｗを引いた差であるΔＧ、ΔＱ、ΔＷについて、これらの行列のフロベニウスノルムを計算する。そして、各更新段階におけるΔＧ、ΔＱ、ΔＷのフロベニウスノルムの、それぞれ最初の更新段階におけるΔＧ、ΔＱ、ΔＷのフロベニウスノルムに対する比が、予め定められた閾値以下となる更新段階をもって、コスト関数Ｊｄの値が収束したとした。そして、Ｊｄが収束した際のＧ、Ｑ、Ｗを、コスト関数Ｊｄを最小化させる行列とした。これにより、データ行列Ｙｄを、式３８に示すように、それぞれ非負値の行列であるＦとＧとＱとＷとに分解した。
【０２０７】
本実験では、以上のようにして、コスト関数Ｊｄを最小化させるＧとＱとＷとを求めることとした。ただし、他の例として、式４６を用いたＧの更新と、式５２を用いたＱの更新と、式５８を用いたＷの更新と、を繰り返し行い、Ｇ、Ｑ、Ｗの更新前と更新後とでコスト関数Ｊｄの値の変動が予め定められた閾値以下となることが、予め定められた回数連続したら、コスト関数Ｊｄの値が収束したこととしてもよい。そして、Ｊｄが収束した際のＧ、Ｑ、Ｗを、コスト関数Ｊｄを最小化させる行列としてもよい。
【０２０８】
図１４Ａ〜図１４Ｂを用いて、疵のない正常な状態にある軸受１０６についての本実験の実験結果について説明する。
図１４Ａのグラフは、Ｆが示す教師基底（に属する係数αのパターン）と、正常な状態にある軸受１０６についての計測データから本実験で求めたＱが示す抽出基底（に属する係数αのパターン）と、を示すグラフである。図１４Ａのグラフの横軸は、係数αのインデックスを示し、縦軸は、対応するデータの値を示す。
【０２０９】
図１４Ｂのグラフは、正常な状態にある軸受１０６についての計測データから本実験で求めたＧとＷとを示すグラフである。図１４Ｂのグラフの横軸は、ＧとＷとの列のインデックスを示し、縦軸は、対応する列のＧとＷとの要素の値を示す。
ＧとＷとの各要素は、データ行列Ｙｄに格納されている係数αのパターンを示すデータそれぞれの中に、教師基底の成分と抽出基底の成分とのそれぞれが、どの程度の重みで存在しているかを示す重みを示す。そのため、図１４Ｂのグラフは、データ行列Ｙｄ内に、教師基底の成分と抽出基底の成分とのそれぞれが、どの程度の重みで存在しているかを示すグラフである。
【０２１０】
図１４Ａのグラフから、教師基底と抽出基底とは、似通っていることが見て取れる。また、図１４Ｂのグラフから、データ行列Ｙｄ内における教師基底の成分と抽出基底の成分との重みには、顕著な違いがみられないことが分かる。図１４Ｂのグラフにおいて、ＧとＷとを、２００個の列の要素それぞれについて比較すると、１１３点においてＧの方が大きくなっており、８７点においてＷの方が大きくなっている。
【０２１１】
図１５Ａ〜図１５Ｂを用いて、図１３Ａに示すように幅３ｍｍ・深さ約０．１ｍｍの線上の疵が外輪に存在する状態にある軸受１０６についての本実験の実験結果について説明する。
【０２１２】
図１５Ａのグラフは、Ｆが示す教師基底（に属する係数αのパターン）と、幅３ｍｍ・深さ約０．１ｍｍの線上の疵が外輪に存在する状態の軸受１０６についての計測データから本実験で求めたＱが示す抽出基底（に属する係数αのパターン）と、を示すグラフである。図１５Ａのグラフの横軸は、係数αのインデックスを示し、縦軸は、対応するデータの値を示す。
【０２１３】
図１５Ｂのグラフは、幅３ｍｍ・深さ約０．１ｍｍの線上の疵が外輪に存在する状態の軸受１０６についての計測データから本実験で求めたＧとＷとを示すグラフである。図１５Ｂのグラフの横軸は、ＧとＷとの列のインデックスを示し、縦軸は、対応する列のＧとＷとの要素の値を示す。
【０２１４】
図１５Ａのグラフから、教師基底と抽出基底とは、顕著に異なることが見て取れる。また、図１５Ｂのグラフから、データ行列Ｙｄ内における教師基底の成分と抽出基底の成分との重みには、顕著な違いがみられることが分かる。抽出基底の成分の重みが、教師基底の成分の重みよりも顕著に高くなっている。図１５Ｂのグラフにおいて、ＧとＷとを、２００個の列の要素それぞれについて比較すると、２００個全てにおいてＷの方が大きくなる。
【０２１５】
図１６Ａ〜図１６Ｂを用いて、図１３Ｂに示すように幅０．５ｍｍ・深さ約０．１ｍｍ以下の線上の疵が外輪に存在する状態にある軸受１０６についての本実験の実験結果について説明する。
【０２１６】
図１６Ａのグラフは、Ｆが示す教師基底と、幅０．５ｍｍ・深さ０．１ｍｍ以下の線上の疵が外輪に存在する状態の軸受１０６についての計測データから本実験で求めたＱが示す抽出基底と、を示すグラフである。図１６Ａのグラフの横軸は、係数αのインデックスを示し、縦軸は、対応するデータの値を示す。
【０２１７】
図１６Ｂのグラフは、幅０．５ｍｍ・深さ０．１ｍｍ以下の線上の疵が外輪に存在する状態の軸受１０６についての計測データから本実験で求めたＧとＷとを示すグラフである。図１６Ｂのグラフの横軸は、ＧとＷとの列のインデックスを示し、縦軸は、対応する列のＧとＷとの要素の値を示す。
【０２１８】
図１６Ａのグラフから、教師基底と抽出基底とは、大まかな波形の動きとしては、似通っているともとれる。しかし、抽出基底に属する基本ベクトルの波形は、教師基底に属する基本ベクトルの波形に比べてより細かな振動が多数存在することが分かる。また、図１６Ｂのグラフから、データ行列Ｙｄ内における教師基底の成分と抽出基底の成分との重みには、顕著な違いがみられることが分かる。図１６Ｂのグラフにおいて、ＧとＷとを、２００個の列の要素それぞれについて比較すると、１９９個においてＷの方が大きくなっており、１個においてＧの方が大きくなっている。
【０２１９】
図１７Ａ〜図１７Ｂを用いて、図１３Ｃに示すように微小な疵が外輪に存在する状態にある軸受１０６についての本実験の実験結果について説明する。
【０２２０】
図１７Ａのグラフは、Ｆが示す教師基底と、微小な疵が外輪に存在する状態の軸受１０６についての計測データから本実験で求めたＱが示す抽出基底と、を示すグラフである。図１７Ａのグラフの横軸は、係数αのインデックスを示し、縦軸は、対応するデータの値を示す。
【０２２１】
図１７Ｂのグラフは、微小な疵が外輪に存在する状態の軸受１０６についての計測データから本実験で求めたＧとＷとを示すグラフである。図１７Ｂのグラフの横軸は、ＧとＷとの列のインデックスを示し、縦軸は、対応する列のＧとＷとの要素の値を示す。
【０２２２】
図１７Ａのグラフから、教師基底と抽出基底とは、似通っているともとれる。また、図１７Ｂのグラフから、データ行列Ｙｄ内における教師基底の成分と抽出基底の成分との重みには、一見、顕著な違いがないようにもとれる。しかし、図１７Ｂのグラフを見ると、Ｗの値は、大半が０．１２付近に存在しており、Ｇの値は、大半が０．０３付近に存在していることが分かる。図１７Ｂのグラフにおいて、ＧとＷとを、２００個の列の要素それぞれについて比較すると、１５９個においてＷの方が大きくなっており、４１個においてＧの方が大きくなっている。図１４Ｂのグラフと比べて、ＧがＷよりも大きくなっている箇所が顕著に多いのが分かる。
【０２２３】
図１４Ａ〜図１７Ｂに示した実験結果から、データ行列Ｙｄに対応する軸受１０６の状態と、教師基底に対応する軸受１０６の状態と、が異なっている場合、データ行列Ｙｄの分解により得られたＷとＧとも顕著に異なる傾向が確認された。
【０２２４】
発明者らは、本実験の実験結果から、データ行列Ｙｄの分解により得られたＧとＷとを比較することで、データ行列Ｙｄに対応する軸受１０６の状態が、教師基底が示す状態であるか否かを診断することができるとの知見を得た。
【０２２５】
そこで、本実施形態の処理は、以下のような処理である。即ち、予め定められた状態の軸受１０６に対応する修正された自己回帰モデルの係数αのパターンを示す非負値のデータを、教師基底に属する基底ベクトルとして求める。そして、診断対象の軸受１０６から計測された計測データからデータ行列Ｙｄを生成し、生成したデータ行列Ｙｄを、教師基底を示す行列と、教師基底の重みを示す教師重み行列と、抽出基底を示す抽出基底行列と、抽出基底の重みを示す抽出重み行列と、に分解する。そして、データ行列Ｙｄの分解により得られた教師重み行列と、抽出重み行列と、に基づいて、軸受１０６の状態が、教師基底に対応する状態であるか否かを診断する処理である。
【０２２６】
（システム構成）
図１８は、本実施形態の診断システムのシステム構成の一例を示す図である。診断システムは、周期的な運動を行う物体についての状態診断を行うシステムである。診断システムは、情報処理装置１８００、振動計測装置１８０１を含む。診断システムは、鉄道台車に利用されている軸受の状態診断を行う。本実施形態では、診断システムは、図１Ａおよび図１Ｂと同様の状況で振動計測位置１１１に設置された振動計測装置１８０１により計測された振動の計測データに基づいて、軸受１０６の状態を診断する。
【０２２７】
情報処理装置１８００は、振動計測装置１８０１により計測された信号に基づいて、軸受の状態診断を行うパーソナルコンピュータ（ＰＣ）、サーバ装置、タブレット装置等の情報処理装置である。また、情報処理装置１８００は、電車に組み込まれたコンピュータ等であってもよい。
【０２２８】
振動計測装置１８０１は、加速度センサ等のセンサを含み、センサを介して物体の振動を検知し、検知した振動に応じた信号を有線又は無線で情報処理装置１８００等の外部の装置に出力する計測装置である。
【０２２９】
（情報処理装置の機能構成）
図１８を参照しながら、情報処理装置１８００が有する機能の一例を説明する。
情報処理装置１８００は、取得部１８１０、決定部１８２０、学習部１８３０、分解部１８４０、診断部１８５０、出力部１８６０を含む。
【０２３０】
≪取得部１８１０≫
取得部１８１０は、周期的な運動を行う物体について、この周期的な運動に係る計測データｙを取得する。
本実施形態では、取得部１８１０は、軸受１０６の内輪を回転させ、振動計測位置１１１に設置された振動計測装置１８０１により計測された振動の計測データｙを取得する。
【０２３１】
≪決定部１８２０≫
決定部１８２０は、取得部１８１０により取得された計測データｙに基づいて、修正された自己回帰モデルにおける係数αを決定する。修正された自己回帰モデルは、計測データｙの実績値ｙｋ−１〜ｙｋ−ｍと、この実績値に対する係数α（＝α１〜αｍ）と、を用いて、計測データｙの予測値ｙ＾ｋを表す式１である。決定部１８２０は、修正された自己回帰モデルにおける係数αを、式１２を用いて決定する。式１２は、一般的に知られている自己回帰モデルにおける係数の決定に用いる式６（ユール・ウォーカー方程式）において、式７で表される自己相関行列Ｒの代わりに、自己相関行列Ｒから物体の状態診断に有用な成分を抽出した第１の行列Ｒ’（式１０）を用いる式である。式６は、時差が０からｍ−１までの計測データｙの自己相関を成分とする自己相関行列Ｒを係数行列とし、時差が１からｍまでの計測データｙの自己相関を成分とする自己相関ベクトルを定数ベクトルとする方程式である。
【０２３２】
式６は、式１で算出される計測データｙの予測値ｙ＾ｋと、計測データｙの予測値ｙ＾ｋに対応する時刻ｋにおける計測データｙの実測値ｙｋとの二乗誤差を最小化する条件を示す条件式として導出することができる。第１の行列Ｒ’は、自己相関行列Ｒの固有値を対角成分とする対角行列Σと、自己相関行列Ｒの固有ベクトルを列成分とする直交行列Ｕと、から導出される。自己相関行列Ｒの固有値は、自己相関行列Ｒを特異値分解（式８）することで導出される。第１の行列Ｒ’は、自己相関行列Ｒの固有値のうち１以上且つｍ未満の設定された使用固有値数ｓ個の固有値を用いて、対角行列Σの部分行列であって、使用固有値数ｓ個の固有値を対角成分とする行列Σｓと、直交行列Ｕの部分行列であって、使用固有値数ｓ個の固有値に対応する固有ベクトルを列成分ベクトルとする行列Ｕｓと、から導出される行列ＵｓΣｓＵｓＴである。使用固有値数ｓ個の固有値は、自己相関行列Ｒの固有値のうち、値が最大の固有値を含むようにすればよいが、値が大きいものから順に選ばれるのが好ましい。
【０２３３】
≪学習部１８３０≫
学習部１８３０は、複数の状態の物体について、予め取得部１８１０により取得された計測データから決定部１８２０により決定された修正された自己回帰モデルにおける係数αを学習データとして、予め定められた状態の物体に対応する修正された自己回帰モデルにおける係数αのパターンを学習により、教師基底に属する基底ベクトルとして特定する。
【０２３４】
≪分解部１８４０≫
分解部１８４０は、診断対象の物体から計測された計測データに基づいて決定部１８２０により決定された係数αを予め定められた個数だけ取得する。分解部１８４０は、取得した係数αそれぞれのパターンを示す非負値のデータを格納するデータ行列を生成する。そして、分解部１８４０は、生成したデータ行列を、学習部１８３０により学習された教師基底を格納する教師基底行列と、データ行列内の教師基底の成分の重みを示す教師重み行列と、データ行列から抽出される基底が格納される抽出基底行列と、データ行列内の抽出基底の成分の重みを示す抽出重み行列と、に分解する。
【０２３５】
≪診断部１８５０≫
診断部１８５０は、分解部１８４０により求められた教師重み行列と抽出重み行列とに基づいて、物体の状態が教師基底に対応する状態であるか否かを診断する。
本実施形態では、軸受１０６の状態が教師基底に対応する状態であるか否かを診断の結果とする。
≪出力部１８６０≫
出力部１８６０は、診断部１８５０による診断の結果に係る情報を出力する。
【０２３６】
（情報処理装置のハードウェア構成）
図１９は、情報処理装置１８００のハードウェア構成の一例を示す図である。
情報処理装置１８００は、ＣＰＵ１９００、主記憶装置１９０１、補助記憶装置１９０２、入出力Ｉ／Ｆ１９０３を含む。各構成要素は、システムバス１９０４を介して、相互に通信可能に接続されている。
【０２３７】
ＣＰＵ１９００は、情報処理装置１８００を制御する中央演算装置である。主記憶装置１９０１は、ＣＰＵ１９００のワークエリアやデータの一時的な保管場所として機能するＲＡＭ（Ｒａｎｄｏｍ  Ａｃｃｅｓｓ  Ｍｅｍｏｒｙ）等の記憶装置である。補助記憶装置１９０２は、各種のプログラム、設定データ、振動計測装置１８０１から出力される計測データ、診断情報等を記憶するＲＯＭ（Ｒｅａｄ  Ｏｎｌｙ  Ｍｅｍｏｒｙ）、ＨＤＤ（Ｈａｒｄ  Ｄｉｓｋ  Ｄｒｉｖｅ）、ＳＳＤ（Ｓｏｌｉｄ  Ｓｔａｔｅ  Ｄｒｉｖｅ）等の記憶装置である。入出力Ｉ／Ｆ１９０３は、振動計測装置１８０１等の外部の装置との間での情報のやり取りに利用されるインターフェースである。
ＣＰＵ１９００が、補助記憶装置１９０２等に記憶されたプログラムに基づき処理を実行することによって、図１８で説明した情報処理装置１８００の機能及び図２０、２１で後述するフローチャートの処理等が実現される。
【０２３８】
（学習処理）
図２０は、診断システムが実行する学習処理の一例を示すフローチャートである。図２０を用いて、教師基底を学習により特定する処理を説明する。
本実施形態では、情報処理装置１８００は、軸受１０６が外輪に疵のある状態である場合、軸受１０６が転動体に疵のある状態である場合、軸受１０６が内輪に疵のある状態である場合、軸受１０６に疵がなく正常な状態である場合、軸受１０６が外輪に疵のある状態である場合のそれぞれについて、予め、ｑ（５０）回、以下の処理を繰り返す。
【０２３９】
即ち、取得部１８１０は、振動計測装置１８０１により計測された振動の計測データｙを取得する。取得部１８１０は、計測データｙとして、ランダムに設定した時刻を時刻１とした場合の時刻１〜Ｍまでの計測データであるｙ１〜ｙＭを取得する。そして、決定部１８２０は、取得した計測データｙと、予め設定された定数Ｍと、修正された自己回帰モデルにおいて、ある時刻のデータを過去の幾つのデータで近似するかを示す数ｍと、に基づいて、式５と式７とを用いて自己相関行列Ｒを生成する。決定部１８２０は、予め記憶されているＭ、ｍの情報を読み出すことで、Ｍとｍとを取得する。本実施形態では、ｍの値は、５００である。また、Ｍは、ｍよりも大きい整数である。
【０２４０】
そして、決定部１８２０は、生成した自己相関行列Ｒを特異値分解することで、式８の直交行列Ｕと対角行列Σを取得し、対角行列Σから自己相関行列Ｒの固有値σ１１〜σｍｍを取得する。決定部１８２０は、自己相関行列Ｒの複数の固有値であるσ１１〜σｍｍのうち、最大のものから使用固有値数ｓだけの固有値σ１１〜σｓｓを、修正された自己回帰モデルの係数αを求めるのに利用する自己相関行列Ｒの固有値として選択する。本実施形態では、使用固有値数ｓは、３とするが、２以下としてもよいし、４以上としてもよい。そして、決定部１８２０は、計測データｙと、固有値σ１１〜σｓｓと、自己相関行列Ｒの特異値分解により得られた直交行列Ｕと、に基づいて、式１３を用いて、修正された自己回帰モデルの係数αを決定する。そして、決定部１８２０は、決定した係数αを、非負値化して、学習データとして補助記憶装置１９０２に記憶する。非負値化は、例えば、係数αのそれぞれの成分に、係数αに含まれる０未満の最低の成分の値に−１をかけた値を足すことにより行われる。
これにより、情報処理装置１８００は、軸受１０６が外輪に疵のある状態である場合、軸受１０６が転動体に疵のある状態である場合、軸受１０６が内輪に疵のある状態である場合、軸受１０６に疵がなく正常な状態である場合、軸受１０６が外輪に疵のある状態である場合のそれぞれについて、ｑ（５０）個ずつ修正された自己回帰モデルにおける係数αのデータを、学習データとして取得し、記憶することができる。
【０２４１】
本実施形態では、軸受１０６が取り得る状態は、疵のない正常な状態、内輪に疵がある状態、保持器に疵がある状態、転動体に疵がある状態、および外輪に疵がある状態の５つであるとする。即ち、軸受１０６が取り得る状態の数ｐは、５となる。
【０２４２】
Ｓ２００１において、学習部１８３０は、補助記憶装置１９０２に記憶されている学習データを取得する。
【０２４３】
Ｓ２００２において、学習部１８３０は、Ｓ２００１で取得した学習データに基づいて、式１９の行列Ｙを生成する。学習部１８３０は、ｍ（５００）×ｄ（２５０）のサイズの行列を生成し、各列のデータを、学習データに含まれる各係数αの値とすることで、行列Ｙを生成する。また、ｍ×ｐのサイズの行列Ａを生成し、行列Ａの各成分に、予め定められた初期値を設定する。また、ｐ×ｄのサイズ行列Ｐを生成し、行列Ｐの各成分に、予め定められた初期値を設定する。
【０２４４】
Ｓ２００３において、学習部１８３０は、Ｓ２００２で生成した行列Ｙと、行列Ａと、行列Ｐと、予め定められた係数βと、予め定められた係数γと、に基づいて、式２１を用いて、コスト関数Ｊの値を、コストとして求める。学習部１８３０は、求めたコストを、更新前コストとして、主記憶装置１９０１等に記憶する。
【０２４５】
Ｓ２００４において、学習部１８３０は、行列Ｙと、行列Ａと、行列Ｐと、係数γと、サイズがｐ×ｐの全ての成分が１の行列である行列Ｅ２と、に基づいて、式３４を用いて、行列Ａを更新する。また、学習部１８３０は、更新した行列Ａを、式３５を用いて、更に更新する。また、学習部１８３０は、行列Ｙと、行列Ａと、行列Ｐと、係数βと、サイズがｍ×ｐの全ての成分が１の行列である行列Ｅ１と、に基づいて、式３０を用いて、行列Ｐを更新する。
【０２４６】
Ｓ２００５において、学習部１８３０は、行列Ｙと、Ｓ２００３で更新した行列Ａと、Ｓ２００３で更新した行列Ｐと、係数βと、係数γと、に基づいて、式２１を用いて、コスト関数Ｊの値を、コストとして求める。学習部１８３０は、求めたコストを、更新後コストとして、主記憶装置１９０１等に記憶する。
【０２４７】
Ｓ２００６において、学習部１８３０は、更新後コストと、更新前コストと、の差分の絶対値が、予め定められた閾値以下であるか否かを判定する。学習部１８３０は、更新後コストと、更新前コストと、の差分が、予め定められた閾値以下であると判定した場合、Ｓ２００７の処理に進む。また、学習部１８３０は、更新後コストと、更新前コストと、の差分の絶対値が、予め定められた閾値よりも大きいと判定した場合、Ｓ２００８の処理に進む。学習部１８３０は、Ｓ２００６の処理の後、更新後コストを、新たな更新前コストとして、主記憶装置１９０１に記憶する。
【０２４８】
Ｓ２００７において、学習部１８３０は、主記憶装置１９０１等に記憶されたカウンタの値に１を加える。学習部１８３０は、このカウンタの値を、図２０の処理の開始前に０に初期化している。
Ｓ２００８において、学習部１８３０は、主記憶装置１９０１等に記憶されたカウンタの値を、０に初期化して、Ｓ２００４の処理に進む。
【０２４９】
Ｓ２００９において、学習部１８３０は、主記憶装置１９０１等に記憶されたカウンタの値が、予め定められた閾値以上であるか否かを判定する。学習部１８３０は、主記憶装置１９０１等に記憶されたカウンタの値が、予め定められた閾値以上であると判定した場合、コスト関数Ｊの値が極小値に収束したとして、Ｓ２０１０の処理に進む。学習部１８３０は、主記憶装置１９０１等に記憶されたカウンタの値が、予め定められた閾値未満であると判定した場合、Ｓ２００４の処理に進む。
【０２５０】
Ｓ２０１０において、学習部１８３０は、行列Ｐに基づいて、行列Ａの各列のデータが示すパターンが、どの状態である軸受１０６に対応する修正された自己回帰モデルにおける係数αのパターンであるかを特定する。
【０２５１】
行列Ｙのある列のデータは、行列Ａと、行列Ｐにおけるその列と、をかけたデータとして求まる。即ち、行列Ｙのある列のデータは、行列Ａの各列のデータを、行列Ｐのその列における各行の値を重みとして、重ね合せたデータとなる。そこで、学習部１８３０は、Ｓ２０１０で以下のようにして、行列Ａの各列のデータが示すパターンが、どの状態である軸受１０６に対応する修正された自己回帰モデルにおける係数αのパターンであるかを特定する。
【０２５２】
学習部１８３０は、行列Ｙの各列のうち、ある１つの状態である軸受１０６に対応する係数αを示すデータを格納する列を特定する。例えば、学習部１８３０は、行列Ｙの各列のうち、外輪に疵がある状態の軸受１０６に対応する係数αを示すデータを格納する列である１列目〜ｑ列目を特定する。そして、学習部１８３０は、行列Ｐにおける特定した列（１列目〜ｑ列目）について、行毎に値を集計する。そして、学習部１８３０は、集計した値のうち、最も大きいものに対応する行を特定する。学習部１８３０は、特定した行列Ｐの行とかけ合わされる行列Ａの列を特定する。例えば、学習部１８３０は、行列Ｐの１行目を特定した場合、行列Ｐの１行目とかけ合わされる行列Ａの列は、１列目なので、１列目を特定する。そして、学習部１８３０は、行列Ａの特定した列が示すデータを、その状態である軸受１０６に対応する係数αの代表的なパターンを示すデータとして特定する。
【０２５３】
学習部１８３０は、軸受１０６が取り得る各状態について、以上の処理を行うことで、行列Ａの各列のデータが示すパターンが、どの状態である軸受１０６に対応する修正された自己回帰モデルにおける係数αのパターンであるかを特定する。
【０２５４】
学習部１８３０は、特定した各状態の軸受１０６に対応する修正された自己回帰モデルにおける係数αの代表的なパターンを示す非負値のデータのうち、予め定められた状態の軸受１０６に対応するデータを、教師基底に属する基底ベクトルとして、補助記憶装置１９０２等に記憶する。本実施形態では、この予め定められた状態は、疵のない正常な状態であるとする。そのため、学習部１８３０は、正常な状態に対応する係数αのパターンを示すデータを、教師基底に属する基底ベクトルとして、補助記憶装置１９０２等に記憶する。
【０２５５】
（診断処理）
本実施形態では、診断システムは、軸受１０６の状態が教師基底に対応する状態か否かを特定することで、軸受１０６の状態を診断することとする。ただし、他の例として、診断システムは、予め振動計測装置１８０１により計測された計測データに基づいて、軸受１０６の状態を診断することとしてもよいし、リアルタイムで振動計測装置１８０１により計測されて出力され続けている計測データを用いて、稼働中の軸受１０６の状態を診断することとしてもよい。
【０２５６】
図２１は、診断処理の一例を示すフローチャートである。
Ｓ２１０１において、取得部１８１０は、振動計測装置１８０１により計測された振動の計測データｙを取得する。取得部１８１０は、計測データｙとして、振動の計測期間内からランダムに時刻を選択する。そして、取得部１８１０は、選択した時刻を時刻１として、時刻１〜Ｍまでの計測データを、ｙ１〜ｙＭとして取得する。
【０２５７】
Ｓ２１０２において、決定部１８２０は、Ｓ２１０１で取得した計測データｙと、予め設定された定数Ｍと、修正された自己回帰モデルにおいて、ある時刻のデータを過去の幾つのデータで近似するかを示す数ｍと、に基づいて、式５と式７とを用いて自己相関行列Ｒを生成する。決定部１８２０は、予め記憶されているＭ、ｍの情報を読み出すことで、Ｍとｍとを取得する。本実施形態では、ｍの値は、５００である。また、Ｍは、ｍよりも大きい整数である。
【０２５８】
Ｓ２１０３において、決定部１８２０は、Ｓ２１０２で生成した自己相関行列Ｒを特異値分解することで、式８の直交行列Ｕと対角行列Σを取得し、対角行列Σから自己相関行列Ｒの固有値σ１１〜σｍｍを取得する。
【０２５９】
Ｓ２１０４において、決定部１８２０は、自己相関行列Ｒの複数の固有値であるσ１１〜σｍｍのうち、最大のものから使用固有値数ｓだけの固有値σ１１〜σｓｓを、修正された自己回帰モデルの係数αを求めるのに利用する自己相関行列Ｒの固有値として選択する。本実施形態では、使用固有値数ｓは、３とするが、２以下としてもよいし、４以上としてもよい。そして、決定部１８２０は、計測データｙと、固有値σ１１〜σｓｓと、自己相関行列Ｒの特異値分解により得られた直交行列Ｕと、に基づいて、式１３を用いて、修正された自己回帰モデルの係数αを決定する。Ｓ２１０４で決定された係数αは、計測データから決定された修正された自己回帰モデルにおける係数である修正係数の一例である。Ｓ２１０４で決定された係数αのパターンは、計測パターンの一例である。
【０２６０】
Ｓ２１０５において、決定部１８２０は、予め定められた数ｑ’個の係数αを決定したか否かを判定する。決定部１８２０は、予め定められた数ｑ’個の係数αを決定したと判定した場合、処理をＳ２１０５に進め、予め定められた数ｑ’個以下の係数αしか決定していないと判定した場合、処理をＳ２１０１に進め、係数αの決定処理を繰り返す。
【０２６１】
Ｓ２１０６において、分解部１８４０は、補助記憶装置１９０２から、Ｓ２０１０で記憶された教師基底を取得する。
【０２６２】
Ｓ２１０７において、分解部１８４０は、決定部１８２０により決定されたｑ’個の係数αを用いて、式３７のように、ｍ×ｑ’のサイズのデータ行列Ｙｄを生成する。即ち、データ行列Ｙｄの各列が、Ｓ２１０４で決定された係数αのパターンのデータを格納することとなる。また、分解部１８４０は、データ行列Ｙｄの各成分を、非負値化する。非負値化は、例えば、０未満の行列Ｙの成分のうち最低のものに−１をかけた値を、行列Ｙの各成分に足すことにより行われる。
【０２６３】
Ｓ２１０８において、分解部１８４０は、教師基底を格納するｍ×１のサイズの教師基底行列Ｆを生成する。そして、分解部１８４０は、Ｆの各要素に、教師基底に属するデータを設定する。また、分解部１８４０は、データ行列Ｙｄに含まれる教師基底の成分の重みを格納する１×ｑ’のサイズの教師重み行列Ｇを生成する。また、分解部１８４０は、データ行列Ｙｄから抽出されるパターンの集合である抽出基底を格納するｍ×予め定められた数ｎのサイズの行列Ｑを生成する。本実施形態では、ｎは１とする。また、分解部１８４０は、データ行列Ｙｄに含まれる抽出基底の成分の重みを格納するｎ×ｑ’のサイズの教師重み行列Ｗを生成する。
そして、分解部１８４０は、Ｇ、Ｑ、Ｗそれぞれの各成分に予め定められた初期値を設定することで、Ｇ、Ｑ、Ｗそれぞれを初期化する。
【０２６４】
Ｓ２１０９において、分解部１８４０は、ＹｄとＦとＧとＱとＷと予め定められた係数δ１、δ２、ε１、ε２と、に基づいて、式３９を用いて、コスト関数Ｊｄの値を、コストとして求める。本実施形態では、δ１、δ２、ε１、ε２それぞれの値は、０とする。分解部１８４０は、求めたコストを、更新前コストとして、主記憶装置１９０１等に記憶する。
【０２６５】
Ｓ２１１０において、分解部１８４０は、Ｇ、Ｑ、Ｗそれぞれを更新する。より具体的には、分解部１８４０は、Ｙｄと、Ｆと、更新前のＧ、Ｑ、Ｗと、予め定められた係数δ１（本実施形態では、０）と、サイズが１×ｑ’の全ての成分が１の行列であるＥ３と、に基づいて、式４６を用いて、Ｇを更新する。ただし、分解部１８４０は、本実施形態ではδ１が０なので、Ｅ３を用いずに、Ｇを更新してもよい。
【０２６６】
また、分解部１８４０は、Ｙｄと、Ｆと、更新前のＧ、Ｑ、Ｗと、予め定められた係数ε１、ε２（本実施形態では、共に０）と、サイズがｎ×ｎの全ての成分が１の行列であるＥ４と、サイズが１×１の全ての成分が１の行列であるＥ５と、に基づいて、式５２を用いて、Ｑを更新する。ただし、分解部１８４０は、本実施形態ではε１、ε２それぞれが０なので、Ｅ４とＥ５とを用いずに、Ｑを更新してもよい。
【０２６７】
また、分解部１８４０は、Ｙｄと、Ｆと、更新前のＧ、Ｑ、Ｗと、予め定められた係数δ２（本実施形態では、０）と、サイズが１×ｑ’の全ての成分が１の行列であるＥ６と、に基づいて、式５８を用いて、Ｗを更新する。ただし、分解部１８４０は、本実施形態ではδ２が０なので、Ｅ６を用いずに、Ｑを更新してもよい。
【０２６８】
Ｓ２１１１において、分解部１８４０は、Ｙｄと、Ｆと、Ｓ２１１０で更新したＧ、Ｑ、Ｗと、係数δ１、δ２、ε１、ε２と、に基づいて、式３９を用いて、コスト関数Ｊｄの値を、コストとして求める。分解部１８４０は、求めたコストを、更新後コストとして、主記憶装置１９０１等に記憶する。
【０２６９】
Ｓ２１１２において、分解部１８４０は、更新後コストと、更新前コストと、の差分の絶対値が、予め定められた閾値以下であるか否かを判定する。分解部１８４０は、更新後コストと、更新前コストと、の差分の絶対値が、予め定められた閾値以下であると判定した場合、処理をＳ２１１３に進める。また、分解部１８４０は、更新後コストと、更新前コストと、の差分の絶対値が、予め定められた閾値よりも大きいと判定した場合、処理をＳ２１１４に進める。分解部１８４０は、Ｓ２１１２の処理の後、更新後コストを、新たな更新前コストとして、主記憶装置１９０１に記憶する。
【０２７０】
Ｓ２１１３において、分解部１８４０は、主記憶装置１９０１等に記憶されたカウンタの値に１を加える。なお、分解部１８４０は、このカウンタの値を、図２１の処理の開始前に０に初期化している。
【０２７１】
Ｓ２１１４において、分解部１８４０は、主記憶装置１９０１等に記憶されたカウンタの値を、０に初期化して、処理をＳ２１１０に進める。
【０２７２】
Ｓ２１１５において、分解部１８４０は、主記憶装置１９０１等に記憶されたカウンタの値が、予め定められた閾値以上であるか否かを判定する。分解部１８４０は、主記憶装置１９０１等に記憶されたカウンタの値が、予め定められた閾値以上であると判定した場合、コスト関数Ｊｄの値が極小値に収束したとして、処理をＳ２１１６に進める。分解部１８４０は、カウンタの値が、予め定められた閾値未満であると判定した場合、処理をＳ２１１０に進める。
【０２７３】
Ｓ２１１６において、診断部１８５０は、処理時点までに更新されたＧとＷとに基づいて、診断対象の軸受１０６の状態が教師基底に対応する状態であるか否かを診断する。診断部１８５０は、ＧとＷとの各列について、要素同士を比較して、何れが大きいかを特定する。そして、診断部１８５０は、Ｇの要素の方がＷの要素よりも大きい列の数が、予め定められた閾値以上である場合、診断対象の軸受１０６の状態が、教師基底に対応する状態（本実施形態では、正常な状態）であると診断する。予め定められた閾値は、例えば、ＧとＷとの列数の７５％の個数（本実施形態では１５０）である。また、診断部１８５０は、Ｇの要素の方がＷの要素よりも大きい列の数が、予め定められた閾値未満である場合、診断対象の軸受１０６の状態が、教師基底に対応する状態ではないと診断する。
【０２７４】
Ｓ２１１７において、出力部１８６０は、Ｓ２１１６での軸受１０６の診断の結果を示す情報を出力する。出力部１８６０は、例えば、表示装置に、軸受１０６の診断の結果を示す情報を表示することで出力する。また、他の例としては、出力部１８６０は、例えば、軸受１０６の診断の結果を示す情報を、補助記憶装置１９０２に記憶することで出力することとしてもよい。また、出力部１８６０は、例えば、外部のサーバ装置等の設定された送信先に、軸受１０６の診断の結果を示す情報を送信することで出力することとしてもよい。
【０２７５】
（まとめ）
診断対象の軸受１０６の計測データから得られた係数αのパターンに含まれる、予め定められた状態に対応するパターンの成分の重みと、それ以外のパターンの成分の重みと、に基づいて、軸受１０６の状態が予め定められた状態か否かを適切に診断できるとの知見が実験３から見出された。本実施形態では、診断システムは、この知見に基づいた軸受１０６の診断処理を行う。
【０２７６】
即ち、診断システムは、予め、疵のない正常な状態の軸受１０６に対応する修正された自己回帰モデルの係数αのパターンを、教師基底に属する基底ベクトルとして求める。そして、診断システムは、診断対象の軸受１０６から計測された計測データからｑ’個の修正された自己回帰モデルの係数αを求め、求めたｑ’個の係数αそれぞれのパターンを示すデータを格納するデータ行列Ｙｄを生成する。そして、診断システムは、Ｙｄを、Ｆと、教師基底の重みを示す教師重み行列Ｇと、抽出基底を示す抽出基底行列Ｑと、抽出基底の重みを示す抽出重み行列Ｗと、に分解する。診断システムは、Ｙｄを分解して得られたＧとＷとに基づいて、診断対象の軸受１０６が教師基底に対応する状態であるか否かを診断する。これにより、診断システムは、軸受１０６が予め定められた状態であるか否かを適切に診断できる。
【０２７７】
また、診断システムは、自己相関行列Ｒの固有値のうち一部の固有値を用いることで、軸受１０６の診断に有用な成分がより多く残り、有用でない成分がより残らないように、計測データを近似する修正された自己回帰モデルの係数αを決定する。診断システムは、このような係数αを用いて、軸受１０６の診断を行うことで、様々な周波数域の信号成分を含んだ信号を用いて診断を行う場合に比べて、より精度よく軸受１０６の診断を行うことができる。
【０２７８】
（変形例１）
本実施形態では、診断システムは、鉄道台車の軸受１０６の状態診断を行うこととした。しかし、診断システムは、歯車等の回転体、振動子等の振動体、バネ等の伸縮体、生物の心臓等の周期的な運動を行う他の物体についての状態を診断することとしてもよい。
【０２７９】
（変形例２）
本実施形態では、ｍは、予め設定された５００という数であるとした。しかし、診断対象に応じて実験を行う等して適切に計測データを近似できる過去のデータの数を特定し、特定した数をｍの値としてもよい。例えば、情報処理装置１８００は、診断対象から計測された計測データについて、ある時刻の計測データを、その時刻よりも過去の計測データを複数用いて自己回帰モデル等を利用して近似して、近似した値とその時刻の計測データとの差分が設定された閾値未満となる場合の近似に用いられた過去の計測データの数を、ｍとして決定してもよい。
【０２８０】
（変形例３）
本実施形態では、情報処理装置１８００は、使用固有値数ｓを３とした。しかし、情報処理装置１８００は、以下のようにして、使用固有値数ｓの値を決定してもよい。即ち、情報処理装置１８００は、まず、自己相関行列Ｒの固有値全ての合計値を求める。そして、情報処理装置１８００は、自己相関行列Ｒの固有値のうち、最大のものからａ個の合計が、求めた合計値の設定された割合（例えば９０％）以上となるａのうち、最小のものを使用固有値数ｓとして決定してもよい。
【０２８１】
（変形例４）
本実施形態では、情報処理装置１８００は、自己相関行列Ｒの固有値のうち、最大のものから、使用固有値数ｓ個の固有値を用いて、修正された自己回帰モデルにおける係数αを決定することとした。しかし、情報処理装置１８００は、自己相関行列Ｒの固有値のうち、任意のｓ個の固有値を用いて、修正された自己回帰モデルにおける係数αを決定することとしてもよい。例えば、情報処理装置１８００は自己相関行列Ｒの固有値のうち、最大のものと、３つ目に大きいものと、の２つを利用して、修正された自己回帰モデルにおける係数αを決定することとしてもよい。
【０２８２】
（変形例５）
本実施形態では、式１９に示す行列Ｙの各列には、複数の状態の軸受１０６に対応する計測データから決定された修正された自己回帰モデルにおける係数αが、状態毎に、ｑ個ずつ格納されることとした。しかし、行列Ｙの各列には、複数の状態の軸受１０６に対応する計測データから決定された修正された自己回帰モデルにおける係数αが、状態毎に、異なる数ずつ格納されることとしてもよい。
【０２８３】
本実施形態では、診断システムは、行列Ｙを、同じ状態の軸受１０６に対応する係数αを格納する列が連続するように、生成した。しかし、診断システムは、同じ状態の軸受１０６に対応する係数αを格納する列が連続しないように、行列Ｙを生成してもよい。
【０２８４】
本実施形態では、診断システムは、各列に各状態の軸受１０６に対応する係数αが格納された行列Ｙを生成し、行列Ｙを非負値行列因子分解により学習することで、各状態の軸受１０６に対応する係数αの代表的なパターンを特定した。しかし、診断システムは、例えば、学習データを取得する際の処理と同様の処理で、ある状態の軸受１０６に対応する係数αを複数決定し、決定した複数の係数αの平均を求めることにより学習することで、求めた平均の係数が示すパターンを、その状態の軸受１０６に対応する係数αの代表的なパターンとして特定してもよい。
【０２８５】
（変形例６）
本実施形態では、診断部１８５０は、Ｙｄの分解により求まったＧとＷとに基づいて、図２１のＳ２１１６で説明した方法で軸受１０６の診断を行うこととした。しかし、診断部１８５０は、他の方法で、診断対象の軸受１０６の状態が教師基底に対応する状態であるか否かを診断してもよい。
【０２８６】
例えば、診断部１８５０は、以下のようにしてもよい。即ち、診断部１８５０は、Ｇの各要素について、偏差を求めて、偏差が小さいものから予め定められた数を特定し、特定したＧの要素の平均値を求める。予め定められた数は、例えば、Ｇの全要素の７０％の数（本実施形態では、１４０）である。また、診断部１８５０は、Ｗの各要素について、偏差を求めて、偏差が小さいものから予め定められた数を特定し、特定したＷの要素の平均値を求める。予め定められた数は、例えば、Ｇの全要素の７０％の数（本実施形態では、１４０）である。そして、診断部１８５０は、求めたＧの要素の平均値がＷの要素の平均値よりも予め定められた閾値以上大きい場合、診断対象の軸受１０６が教師基底に対応する状態であると診断する。一方、診断部１８５０は、求めたＧの要素の平均値がＷの要素の平均値よりも予め定められた閾値以上大きくない場合、診断対象の軸受１０６が教師基底に対応する状態でないと診断する。
【０２８７】
また、診断部１８５０は、以下のようにしてもよい。即ち、診断部１８５０は、Ｇの各要素のうち値が中間の予め定められた数の要素の平均値を求める。予め定められた数は、例えば、Ｇの全要素の７０％の数（本実施形態では、１４０）である。また、診断部１８５０は、Ｗの各要素のうち値が中間の予め定められた数の要素の平均値を求める。そして、診断部１８５０は、求めたＧの要素の平均値がＷの要素の平均値よりも予め定められた閾値以上大きい場合、診断対象の軸受１０６が教師基底に対応する状態であると診断する。一方、診断部１８５０は、求めたＧの要素の平均値がＷの要素の平均値よりも予め定められた閾値以上大きくない場合、診断対象の軸受１０６が教師基底に対応する状態でないと診断する。
【０２８８】
（変形例７）
本実施形態では、データ行列Ｙｄから抽出される、抽出基底に属する基底ベクトルの数（抽出基底行列Ｑの列数、抽出重み行列Ｗの行数）を１とした。ただし、データ行列Ｙｄから抽出される、抽出基底に属する基底ベクトルの数を２以上としてもよい。その場合、Ｓ２１１６では、診断部１８５０は、以下のようにしてもよい。
【０２８９】
即ち、診断部１８５０は、診断部１８５０は、Ｇと、Ｗの各行それぞれと、の各列について、要素同士を比較して、何れが大きいかを特定する。そして、診断部１８５０は、Ｇの要素が最も大きい列の数が、予め定められた閾値以上である場合、診断対象の軸受１０６の状態が、教師基底に対応する状態であると診断する。一方、診断部１８５０は、Ｇの要素が最も大きい列の数が、予め定められた閾値未満である場合、診断対象の軸受１０６の状態が、教師基底に対応する状態ではないと診断する。
【０２９０】
（変形例８）
本実施形態では、診断システムは、各列に係数αのパターンを示す非負値のデータを格納するように行列Ｙｄを生成した。ただし、他の例として、診断システムは、各列ではなく各行に係数αのパターンを示す非負値のデータを格納するように行列Ｙｄを生成してもよい。その場合、診断システムは、行列Ｙｄを、Ｙｄ＝ＧＴＦＴ＋ＷＴＱＴのように、分解すればよい。その場合、診断システムは、コスト関数Ｊｄとして、式３９における「ＦＧ＋ＱＷ−Ｙｄ」の部分を「ＧＴＦＴ＋ＷＴＱＴ−Ｙｄ」に置き換えた関数を用いればよい。
【０２９１】
（変形例９）
本実施形態では、診断システムは、診断対象の軸受１０６が疵のない正常な状態であるか否かを診断することとした。ただし、診断システムは、診断対象の軸受１０６が他の状態（例えば、外輪に特定の疵の有る状態、内輪に疵の有る状態、転動体に疵の有る状態等）であるか否かを診断することとしてもよい。
【０２９２】
（変形例１０）
本実施形態では、診断システムは、式３９で表される関数Ｊｄを、コスト関数として用いることとした。ただし、診断システムは、Ｇ、Ｑ、Ｗの更新の指標となる関数であれば、他の関数をコスト関数として用いてもよい。例えば、診断システムは、式３９で表される関数Ｊｄの逆数で表される関数をコスト関数として用いてもよい。その場合、診断システムは、コスト関数の値を適正化（最大化）するため、Ｇ、Ｑ、Ｗの値をコスト関数の値が増大するように更新していくこととなる。そして、診断システムは、コスト関数が収束したら、コスト関数の値が適正化（最大化）できたとして、その際のＧ、Ｑ、Ｗを、Ｙｄの分解結果として取得する。
【０２９３】
ここで、最小化とは、コスト関数Ｊｄの値を最小化する解を導出する最適化アルゴリズムにおいてコスト関数Ｊｄを最小化することを指す。即ち、最小化とは、コスト関数の値の厳密な最小値を求めるものに限られない。このことは、最大化についても同じである。
【０２９４】
（変形例１１）
本実施形態では、診断システムは、図２０の学習処理において、行列Ｐと行列Ａとの更新前後で、コスト関数Ｊの値の変動が予め定められた閾値以下となることが、予め定められた回数連続するまで、行列Ｐと行列Ａとの更新を繰り返した。これにより、診断システムは、コスト関数Ｊが最小化した際の行列Ｐと行列Ａとを求めることとした。ただし、他の例として、診断システムは、以下のようにしてもよい。
【０２９５】
即ち、診断システムは、式３０を用いて行列Ｐを更新し、式３４、式３５を用いて行列Ａを更新することを繰り返す。そして、診断システムは、各更新段階において、更新後の行列Ｐ及び行列Ａから、それぞれ更新前の行列Ｐ及び行列Ａを引いた差である行列△Ｐ及び行列△Ａを求める。診断システムは、行列△Ｐ及び行列△Ａそれぞれのフロベニウスノルムを計算する。診断システムは、各更新段階における行列△Ｐ及び行列△Ａのフロベニウスノルムの、それぞれ最初の更新段階における行列△Ｐ及び行列△Ａのフロベニウスノルムに対する比が、予め定められた閾値以下となる更新段階をもって、コスト関数Ｊの値が収束したとする。そして、診断システムは、その際の行列Ａと行列Ｐとをコスト関数Ｊを最小化させる行列として求める。
【０２９６】
（変形例１２）
本実施形態では、診断システムは、図２１の診断処理において、Ｇ、Ｑ、Ｗの更新前と更新後とでコスト関数Ｊｄの値の変動が予め定められた閾値以下となることが、予め定められた回数連続するまで、Ｇ、Ｑ、Ｗそれぞれの更新を繰り返すことで、コスト関数Ｊｄの値が収束した際のＧ、Ｑ、Ｗを求めることとした。ただし、他の例として、診断システムは、以下のようにしてもよい。
【０２９７】
即ち、診断システムは、式４６を用いたＧの更新と、式５２を用いたＱの更新と、式５８を用いたＷの更新と、を繰り返し行う。診断システムは、各更新段階において、更新後のＧ、Ｑ、Ｗから、それぞれ更新前のＧ、Ｑ、Ｗを引いた差であるΔＧ、ΔＱ、ΔＷについて、これらの行列のフロベニウスノルムを計算する。診断システムは、各更新段階におけるΔＧ、ΔＱ、ΔＷのフロベニウスノルムの、それぞれ最初の更新段階におけるΔＧ、ΔＱ、ΔＷのフロベニウスノルムに対する比が、予め定められた閾値以下となる更新段階をもって、コスト関数Ｊｄの値が収束したとしてもよい。そして、診断システムは、Ｊｄが収束した際のＧ、Ｑ、Ｗを、コスト関数Ｊｄを最小化させる行列として求めてもよい。
【０２９８】
＜実施形態２＞
本実施形態では、コスト関数Ｊｄにおける係数δ１、δ２、ε１、ε２のうちの少なくとも１つが０以外の値である場合の処理について説明する。
【０２９９】
（実験４）
発明者らは、実験３と同様の実験状況において、以下のような実験を行った。即ち、疵のない正常な状態にある軸受１０６と、図１３Ｃに示す外輪に微小な疵のある状態の軸受１０６と、について実験３で説明した実験と同様の実験を、式３９に示すコスト関数Ｊｄのδ１、δ２、ε２それぞれについて、値を変更しつつ行った。本実験では、使用固有値数ｓは、３とした。また、ｍは、５００とした。
【０３００】
本実験では、δ２、ε１、ε２それぞれを０として、δ１の値を、０、０．０００２５、０．０００５、０．０００７５、０．００１それぞれに変更しつつ、実験３で説明した実験と同様の実験を行った。
【０３０１】
また、本実験では、δ１、ε１、ε２それぞれを０として、δ２の値を、０、０．００００５、０．０００２５、０．０００５それぞれに変更しつつ、実験３で説明した実験と同様の実験を行った。
【０３０２】
また、本実験では、δ１、δ２、ε１それぞれを０として、ε２の値を、０、０．０２５、０．０５、０．１それぞれに変更しつつ、実験３で説明した実験と同様の実験を行った。
【０３０３】
δ２、ε１、ε２それぞれを０として、δ１の値を変更しつつ、実験３で説明した実験と同様の実験を行った結果について説明する。δ１は、コスト関数Ｊｄにおける教師重み行列Ｇのスパース性の向上に関する制約の項に係る係数である。以下では、スパース性の向上に関する制約を、スパース制約とする。即ち、δ１の値は、Ｙｄの分解におけるＧのスパース制約の度合を示す係数である。
【０３０４】
図２２Ａ〜図２２Ｅに、外輪に微小な疵のある状態の軸受１０６についての実験結果を示す。また、図２３Ａ〜図２３Ｅに、正常な状態の軸受１０６についての実験結果を示す。図２２Ａ〜図２２Ｅ、および図２３Ａ〜図２３Ｅの各グラフは、データ行列Ｙｄの分解の結果求まったＧとＷとを示すグラフである。横軸は、ＧとＷとの各列のインデックスを示す。縦軸は、ＧとＷとの対応する列の要素の値を示す。
【０３０５】
図２２Ａ、図２３Ａのグラフは、δ１が０の場合の実験結果を示すグラフである。図２２Ｂ、図２３Ｂのグラフは、δ１が０．０００２５の場合の実験結果を示すグラフである。図２２Ｃ、図２３Ｃのグラフは、δ１が０．０００５の場合の実験結果を示すグラフである。図２２Ｄ、図２３Ｄのグラフは、δ１が０．０００７５の場合の実験結果を示すグラフである。図２２Ｅ、図２３Ｅのグラフは、δ１が０．００１の場合の実験結果を示すグラフである。
【０３０６】
図２２Ａ〜図２２Ｅの各グラフを見ると、δ１が０．０００２５、０．０００５、０．０００７５それぞれの場合、δ１が０の場合に比べて、Ｇの値とＷの値とは、より分離しており、正常な状態と微小な疵のある状態との違いがより表れており、診断により有利になっていることが見て取れる。ただし、δ１が０．００１の場合、δ１が０．０００７５の場合に比べて、Ｇの値とＷの値とは、より近くなっており、正常な状態と微小な疵のある状態との違いが小さくなり、診断に不利になっていしまっていることが見て取れる。
【０３０７】
図２３Ａ〜図２３Ｅの各グラフを見ると、δ１が０．０００２５、０．０００５、０．０００７５、０．００１それぞれの場合で、Ｇの値とＷの値との分離の度合に顕著な違いが表れていないことが見て取れる。
【０３０８】
発明者らは、図２２Ａ〜図２２Ｅ、図２３Ａ〜図２３Ｅの各グラフから、δ１が０．０００２５、０．０００５、０．０００７５の場合において、軸受１０６の診断により有利であることを見出した。これにより、発明者らは、δ１の適切な値は、０．０００２５以上０．０００７５以下の範囲の値であるとの知見を得た。
【０３０９】
次に、δ１、ε１、ε２それぞれを０として、δ２の値を変更しつつ、求まった実験結果について説明する。δ２は、コスト関数Ｊｄにおける抽出重み行列Ｗのスパース制約の項に係る係数である。即ち、δ２の値は、Ｙｄの分解におけるＷのスパース制約の度合を示す係数である。
【０３１０】
図２４Ａ〜図２４Ｄに、外輪に微小な疵のある状態の軸受１０６についての実験結果を示す。また、図２５Ａ〜図２５Ｄに、正常な状態の軸受１０６についての実験結果を示す。図２４Ａ〜図２４Ｄ、図２５Ａ〜図２５Ｄの各グラフは、データ行列Ｙｄの分解の結果求まったＧとＷとを示すグラフである。横軸は、ＧとＷとの各列のインデックスを示す。縦軸は、ＧとＷとの対応する列の要素の値を示す。
【０３１１】
図２４Ａ、図２５Ａのグラフは、δ２が０の場合の実験結果を示すグラフである。図２４Ｂ、図２５Ｂのグラフは、δ２が０．００００５の場合の実験結果を示すグラフである。図２４Ｃ、図２５Ｃのグラフは、δ２が０．０００２５の場合の実験結果を示すグラフである。図２４Ｄ、図２５Ｄのグラフは、δ２が０．０００５の場合の実験結果を示すグラフである。
【０３１２】
図２４Ａ〜図２４Ｄの各グラフを見ると、δ２が０．００００５の場合、δ２が０の場合と同等に、Ｇの値とＷの値とが分離していることが見て取れる。また、δ２が０．０００２５、０．０００５の場合、δ２が０、０．００００５の場合に比べて、Ｇの値とＷの値とがより近くなっており、正常な状態と微小な疵のある状態との違いが小さくなり、診断に不利になっていしまっていることが見て取れる。
【０３１３】
図２５Ａ〜図２５Ｄの各グラフを見ると、δ２が０、０．００００５それぞれの場合で、Ｇの値とＷの値との分離の度合に顕著な違いが表れていないことが見て取れる。また、δ２が０．０００２５、０．０００５それぞれの場合で、δ２が０、０．００００５の場合に比べて、Ｇの値とＷの値とがより分離してしまっており、正常な状態同士で違いが大きくなり、診断に不利になっていしまっていることが見て取れる。
【０３１４】
発明者らは、図２４Ａ〜図２４Ｂ、図２５Ａ〜図２５Ｂの各グラフから、δ２が０．０００２５、０．０００５それぞれの場合において、軸受１０６の診断により不利であることを見出した。これにより、発明者らは、δ２の適切な値は、０以上０．００００５以下（好ましくは０以上０．００００５未満）の範囲の値であるとの知見を得た。
【０３１５】
次に、δ１、δ２、ε１それぞれを０として、ε２の値を変更しつつ、求まった実験結果について説明する。ε２は、コスト関数Ｊｄにおける教師基底と抽出基底との直交性の向上に関する制約の項に係る係数である。以下では、直交性の向上に関する制約を、直交制約とする。即ち、ε２の値は、Ｙｄの分解におけるＦとＱとの直交制約の度合を示す係数である。
【０３１６】
図２６Ａ〜図２６Ｄに、外輪に微小な疵のある状態の軸受１０６についての実験結果を示す。また、図２７Ａ〜図２７Ｄに、正常な状態の軸受１０６についての実験結果を示す。図２６Ａ〜図２６Ｄ、図２７Ａ〜図２７Ｄの各グラフは、データ行列Ｙｄの分解の結果求まったＧとＷとを示すグラフである。横軸は、ＧとＷとの各列のインデックスを示す。縦軸は、ＧとＷとの対応する列の要素の値を示す。
【０３１７】
図２６Ａ、図２７Ａのグラフは、ε２が０の場合の実験結果を示すグラフである。図２６Ｂ、図２７Ｂのグラフは、ε２が０．０２５の場合の実験結果を示すグラフである。図２６Ｃ、図２６Ｃのグラフは、ε２が０．０５の場合の実験結果を示すグラフである。図２６Ｄ、図２７Ｄのグラフは、ε２が０．１の場合の実験結果を示すグラフである。
【０３１８】
図２６Ａ〜図２６Ｄの各グラフを見ると、ε２が０．０５、０．１の場合、ε２が０、０．０２５の場合に比べて、Ｇの値とＷの値とがより近くなっており、正常な状態と微小な疵のある状態との違いが小さくなり、診断に不利になっていしまっていることが見て取れる。
【０３１９】
図２７Ａ〜図２７Ｄの各グラフを見ると、ε２が０．１の場合、ε２が０、０．０２５、０．０５の場合に比べて、Ｇの値とＷの値とがより分離してしまっており、正常な状態同士で違いが大きくなり、診断に不利になっていしまっていることが見て取れる。
【０３２０】
発明者らは、図２６Ａ〜図２６Ｄ、図２７Ａ〜図２７Ｄの各グラフから、ε２が０．０５、０．１それぞれの場合において、軸受１０６の診断により不利であることを見出した。これにより、発明者らは、ε２の適切な値は、０以上０．０２５以下の範囲の値であるとの知見を得た。
【０３２１】
（診断処理）
本実施形態の診断システムのシステム構成は、実施形態１と同様である。また、情報処理装置１８００のハードウェ構成及び機能構成についても、実施形態１と同様である。
【０３２２】
本実施形態の処理は、診断システムがδ１、δ２、ε１、ε２の全てが０のコスト関数Ｊｄの代わりにδ１、δ２、ε１、ε２のうちの少なくとも１つが０でないコスト関数Ｊｄを用いること以外は、図２０、図２１で説明した実施形態１の処理と同様である。
【０３２３】
本実施形態では、δ１の値は、実験により見出された適切な値の範囲（０．０００２５以上０．０００７５以下の範囲）から予め定められた値（例えば、０．０００２５、０．０００５、０．０００７５等）である。また、δ２の値は、実験により見出された適切な値の範囲（０以上０．００００５以下の範囲）から予め定められた値（例えば、０、０．００００２５、０．００００５等）である。また、ε２の値は、実験により見出された適切な値の範囲（０以上０．０２５以下の範囲）から予め定められた値（例えば、０、０．０１、０．０２５等）である。また、ε１の値は、予め定められた値（例えば、０等）である。
【０３２４】
δ１、δ２、ε１、ε２のうちの少なくとも１つが０でないコスト関数Ｊｄは、δ１、δ２、ε１、ε２のうちの０でない係数に対応する制約を付加した関数となる。そのため、δ１、δ２、ε１、ε２のうちの少なくとも１つが０でないコスト関数Ｊｄは、教師重み行列Ｇについてのスパース制約、抽出重み行列についてのスパース制約、抽出基底の正規直交制約、および教師基底と抽出基底との間の直交制約のうちの少なくとも１つを示す。
【０３２５】
（まとめ）
本実施形態の処理により、診断システムは、診断対象の軸受１０６が予め定められた状態であるか否かをより精度よく診断できる。尚、理論上は、式３８に示す分解の自由度を低減させるために、コスト関数に抽出基底の正規直交制約と教師基底と抽出基底との間の直交制約とを付加するのが好ましい。本実施形態においても、実施形態１で説明した変形例を採用することができる。
【０３２６】
＜その他の実施形態＞
本実施形態では、情報処理装置１８００は、補助記憶装置１９０２に記憶されたプログラムを実行することで、図２０、図２１の処理を実現することとした。しかし、情報処理装置１８００は、外付けの記憶媒体や外部の記憶サーバ等に記憶されたプログラムを実行することで、図２０、図２１の処理を実現することとしてもよい。
【０３２７】
また、例えば、上述した情報処理装置１８００の機能の一部又は全てをハードウェアとして情報処理装置１８００に実装してもよい。
また、以上説明した本発明の実施形態は、何れも本発明を実施するにあたっての具体化の例を示したものに過ぎず、これらによって本発明の技術的範囲が限定的に解釈されてはならないものである。即ち、本発明はその技術思想、またはその主要な特徴から逸脱することなく、様々な形で実施することができる。
【産業上の利用可能性】
【０３２８】
本発明は、例えば、周期的な運動を行う物体の状態を診断することに利用することができる。
【図１Ａ】
2020145215.tif 000061
【図１Ｂ】
2020145215.tif 000062
【図２】
2020145215.tif 000063
【図３Ａ】
2020145215.tif 000064
【図３Ｂ】
2020145215.tif 000065
【図３Ｃ】
2020145215.tif 000066
【図３Ｄ】
2020145215.tif 000067
【図３Ｅ】
2020145215.tif 000068
【図４Ａ】
2020145215.tif 000069
【図４Ｂ】
2020145215.tif 000070
【図４Ｃ】
2020145215.tif 000071
【図４Ｄ】
2020145215.tif 000072
【図４Ｅ】
2020145215.tif 000073
【図５Ａ】
2020145215.tif 000074
【図５Ｂ】
2020145215.tif 000075
【図５Ｃ】
2020145215.tif 000076
【図５Ｄ】
2020145215.tif 000077
【図５Ｅ】
2020145215.tif 000078
【図６Ａ】
2020145215.tif 000079
【図６Ｂ】
2020145215.tif 000080
【図６Ｃ】
2020145215.tif 000081
【図６Ｄ】
2020145215.tif 000082
【図６Ｅ】
2020145215.tif 000083
【図７Ａ】
2020145215.tif 000084
【図７Ｂ】
2020145215.tif 000085
【図７Ｃ】
2020145215.tif 000086
【図７Ｄ】
2020145215.tif 000087
【図７Ｅ】
2020145215.tif 000088
【図８Ａ】
2020145215.tif 000089
【図８Ｂ】
2020145215.tif 000090
【図８Ｃ】
2020145215.tif 000091
【図８Ｄ】
2020145215.tif 000092
【図８Ｅ】
2020145215.tif 000093
【図９】
2020145215.tif 000094
【図１０】
2020145215.tif 000095
【図１１】
2020145215.tif 000096
【図１２】
2020145215.tif 000097
【図１３Ａ】
2020145215.tif 000098
【図１３Ｂ】
2020145215.tif 000099
【図１３Ｃ】
2020145215.tif 000100
【図１４Ａ】
2020145215.tif 000101
【図１４Ｂ】
2020145215.tif 000102
【図１５Ａ】
2020145215.tif 000103
【図１５Ｂ】
2020145215.tif 000104
【図１６Ａ】
2020145215.tif 000105
【図１６Ｂ】
2020145215.tif 000106
【図１７Ａ】
2020145215.tif 000107
【図１７Ｂ】
2020145215.tif 000108
【図１８】
2020145215.tif 000109
【図１９】
2020145215.tif 000110
【図２０】
2020145215.tif 000111
【図２１】
2020145215.tif 000112
【図２２Ａ】
2020145215.tif 000113
【図２２Ｂ】
2020145215.tif 000114
【図２２Ｃ】
2020145215.tif 000115
【図２２Ｄ】
2020145215.tif 000116
【図２２Ｅ】
2020145215.tif 000117
【図２３Ａ】
2020145215.tif 000118
【図２３Ｂ】
2020145215.tif 000119
【図２３Ｃ】
2020145215.tif 000120
【図２３Ｄ】
2020145215.tif 000121
【図２３Ｅ】
2020145215.tif 000122
【図２４Ａ】
2020145215.tif 000123
【図２４Ｂ】
2020145215.tif 000124
【図２４Ｃ】
2020145215.tif 000125
【図２４Ｄ】
2020145215.tif 000126
【図２５Ａ】
2020145215.tif 000127
【図２５Ｂ】
2020145215.tif 000128
【図２５Ｃ】
2020145215.tif 000129
【図２５Ｄ】
2020145215.tif 000130
【図２６Ａ】
2020145215.tif 000131
【図２６Ｂ】
2020145215.tif 000132
【図２６Ｃ】
2020145215.tif 000133
【図２６Ｄ】
2020145215.tif 000134
【図２７Ａ】
2020145215.tif 000135
【図２７Ｂ】
2020145215.tif 000136
【図２７Ｃ】
2020145215.tif 000137
【図２７Ｄ】
2020145215.tif 000138
【国際調査報告】
2020145215.tif 000139
2020145215.tif 000140
2020145215.tif 000141
2020145215.tif 000142