JP6919397B2 - 情報処理装置、情報処理方法及びプログラム - Google Patents

情報処理装置、情報処理方法及びプログラム Download PDF

Info

Publication number
JP6919397B2
JP6919397B2 JP2017151727A JP2017151727A JP6919397B2 JP 6919397 B2 JP6919397 B2 JP 6919397B2 JP 2017151727 A JP2017151727 A JP 2017151727A JP 2017151727 A JP2017151727 A JP 2017151727A JP 6919397 B2 JP6919397 B2 JP 6919397B2
Authority
JP
Japan
Prior art keywords
matrix
measurement data
eigenvalues
autocorrelation
bearing
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
JP2017151727A
Other languages
English (en)
Other versions
JP2018163135A (ja
Inventor
中川 淳一
淳一 中川
秀樹 南
秀樹 南
啓之 上島
啓之 上島
章典 竹下
章典 竹下
谷 泰徳
泰徳 谷
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Nippon Steel Corp
Original Assignee
Nippon Steel Corp
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Nippon Steel Corp filed Critical Nippon Steel Corp
Priority to PCT/JP2018/028942 priority Critical patent/WO2019026980A1/ja
Priority to EP18840775.3A priority patent/EP3663741A4/en
Priority to CN201880052067.6A priority patent/CN110998272A/zh
Publication of JP2018163135A publication Critical patent/JP2018163135A/ja
Application granted granted Critical
Publication of JP6919397B2 publication Critical patent/JP6919397B2/ja
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Landscapes

  • Testing Of Devices, Machine Parts, Or Other Structures Thereof (AREA)

Description

本発明は、情報処理装置、情報処理方法及びプログラムに関する。
周期的な運動を行う物体の異常診断を行う技術として、運動を行っている物体から得られる信号を測定し、その信号の特徴量を算出し、その特徴量から正常であるか異常であるかの判定を行う方法がある。このとき、正常時と異常時との特徴量の差異を明確にすることによって異常診断を精度よく行うために、その物体の運動の基本周波数の整数倍に対応する狭帯域周波数領域に属する成分を測定された信号から抽出し、抽出した成分の信号から得られる自己回帰モデルの係数を算出し、この自己回帰モデルの係数から特徴量を算出する技術がある(特許文献1)。
特開2003−50158号公報
特許文献1に開示された方法では、測定された信号に対して、予めフィルタ処理等を施すことで、特徴量を算出するための信号を抽出する。しかしながら、特許文献1に開示された方法では、通過させるフィルタの周波数領域を予め設定しておく必要があるため、物体の異常を反映している信号成分の周波数が想定外で、予め設定した周波数領域から外れた場合は、正常時と異常時とで特徴量に差異が現れないため、異常を検知出来ない可能性があるという問題がある。
本発明は、以上のような問題点に鑑みてなされたものであり、周期的な運動を行う物体の異常診断をより精度よく行うことができるようにすることを目的とする。
本発明の情報処理装置は、周期的な運動を行う物体の前記周期的な運動に係る計測データを取得する取得手段と、前記取得手段により取得された前記計測データに基づいて、修正された自己回帰モデルにおける係数を決定する決定手段と、前記決定手段により決定された前記係数に基づいて、前記物体の異常を診断する診断手段と、を有し、前記修正された自己回帰モデルは、前記計測データの実績値と、前記実績値に対する前記係数と、を用いて、前記計測データの予測値を表す式であり、前記決定手段は、前記係数から成るベクトルを未知ベクトルとし、前記計測データから導出される自己相関行列を特異値分解することで導出される前記自己相関行列の固有値を対角成分とする対角行列と、前記自己相関行列の固有ベクトルを列成分とする直交行列と、から導出される第1の行列を係数行列とし、前記計測データから導出される自己相関ベクトルを定数ベクトルとする線形方程式を用いて、前記係数を決定し、前記自己相関ベクトルは、時差が1から前記修正された自己回帰モデルで用いられる前記実績値の数であるmまでの前記計測データの自己相関を成分とするベクトルであり、前記自己相関行列は、時差が0からm−1までの前記計測データの自己相関を成分とする行列であり、前記第1の行列は、1以上且つm未満の設定された数であるsに対して、前記自己相関行列のs個の固有値と前記対角行列とから導出される第2の行列Σsと、前記s個の固有値と前記直交行列とから導出される第3の行列Usと、から導出される行列UsΣss Tであり、前記第2の行列は、前記対角行列の部分行列であって、前記s個の固有値を対角成分とする行列であり、前記第3の行列は、前記直交行列の部分行列であって、前記s個の固有値に対応する固有ベクトルを列成分ベクトルとする行列であり、前記診断手段は、前記修正された自己回帰モデルの周波数の分布を示す周波数特性を、前記決定手段により決定された係数を用いて導出し、前記周波数特性に基づいて、前記物体の異常を診断する。
本発明によれば、周期的な運動を行う物体の異常診断をより精度よく行うことができる技術を提供することができる。
図1は、実験の状況を示す図である。 図2は、軸受の詳細の一例を説明する図である。 図3は、固有値の分布の一例を示す図である。 図4は、計測データの波形の一例を示す図である。 図5は、修正された自己回帰モデルの波形の一例を示す図である。 図6は、修正された自己回帰モデルの周波数特性の一例を示す図である。 図7は、修正された自己回帰モデルの波形の一例を示す図である。 図8は、修正された自己回帰モデルの周波数特性の一例を示す図である。 図9は、診断システムのシステム構成の一例を示す図である。 図10は、情報処理装置のハードウェア構成の一例を示す図である。 図11は、情報処理装置の処理の一例を示すフローチャートである。 図12は、固有値の分布の一例を示す図である。 図13は、計測データの波形の一例を示す図である。 図14は、修正された自己回帰モデルの波形の一例を示す図である。 図15は、修正された自己回帰モデルの周波数特性の一例を示す図である。 図16は、修正された自己回帰モデルの波形の一例を示す図である。 図17は、修正された自己回帰モデルの周波数特性の一例を示す図である。
<実施形態1>
以下、本発明の一実施形態について図面に基づいて説明する。
(本実施形態の処理の概要)
本実施形態では、診断システムが鉄道台車に利用される軸受から測定された振動の信号に基づいて、その軸受の異常を診断する。
診断システムは、軸受の振動に応じた信号を測定し、測定された信号を自己回帰モデルで表現し、自己回帰モデルの係数に関する条件式における行列を特異値分解し、その行列の固有値を求め、求めた固有値のうち最大のものから設定された数のみを用いて、自己回帰モデルの係数を決定する。そして、診断システムは、決定した係数から、自己回帰モデルの周波数特性を求め、求めた周波数特性に基づいて、軸受の異常を診断する。
(軸受異常の診断実験)
鉄道台車に利用されている軸受の異常診断精度の評価のために行った実験について説明する。
図1(a)は、本実験の状況を説明する図である。図1(a)の状況は、実験室に鉄道台車における車輪の駆動機構が設置されている状況である。駆動モータ101は、駆動モータ軸102、歯車継手103、小歯車軸100を介して、小歯車105を回転させる。小歯車105は、回転することで、大歯車107を回転させる。それに合わせて、車軸109が回転し、車軸109に接続される車輪が回転することになる。本実験では、車軸109には、車輪の代わりにダイナモ110が接続されている。ダイナモ110は、疑似的な走行負荷を与えるために接続されている発電機である。本実験では、ダイナモ110を回転させるための負荷を、走行の際の負荷として想定する。
小歯車105が嵌められた小歯車軸100と大歯車107が嵌められた車軸109とは、歯車箱104に固定されている。そして、歯車箱104には、小歯車軸との間に軸受106が装着されており、車軸109との間に軸受108が装着されている。図1(b)は、歯車箱104に取り付けられている軸受を側面から見た様子を示す図である。
歯車箱104上の位置111は、振動計測装置が設置される位置である。振動計測装置は、加速度センサ、レーザ変位計等のセンサを含み、センサを介して物体の振動を検知し、検知した振動に応じた信号を出力する。位置111に設置された振動計測装置が位置111に伝達される振動を計測し、計測した振動を示す信号を外部の情報処理装置等に有線又は無線による通信を介して出力する。本実験では、振動計測装置は、位置111に設置されるとするが、軸受106に起因する振動が伝達される位置であるならば、任意の位置に設置することとしてもよい。
ここで、軸受106について、説明する。図2は、軸受106の一例の構成を説明する図である。軸受106は、外輪、内輪、転動体、保持器の4つの部分から構成される。図2には、軸受106の外輪、内輪、転動体、保持器の概要が示されている。軸受106は、外輪と内輪との間に保持器により保持された転動体が挟まれる構成となっている。
本実験では、小歯車105に装着されている軸受106として、疵のない正常な軸受、内輪に疵のある軸受、保持器に疵のある軸受のそれぞれを利用した。それぞれの軸受106毎に、駆動モータ101を駆動させ、位置111に設置された振動計測装置が振動を計測することで、軸受106の異常診断を行うための計測データを取得した。本実験では、駆動モータ101は、軸受106の内輪の回転数が2954.7rpm(round per minute)となるように駆動した。
以上のように、軸受106が正常な場合、軸受106の内輪に疵がある場合、軸受106の保持器に疵がある場合について、位置111における振動の計測データを取得した。
次に、振動計測装置が計測した計測データを用いた軸受106の異常診断の方法を説明する。以下では、振動計測装置が計測した計測データを、計測データyとする。
本実験では、計測データyを、修正された自己回帰モデルで近似し、この修正された自己回帰モデルの係数から特徴量を算出することとした。
ある時刻k(1≦k≦M)における計測データyの値をykとする。Mは、計測データyがどの時刻までのデータを含むかを示す数であり、予め設定されている。ykを近似する自己回帰モデルは、例えば、以下の式1のようになる。式1に示すように、自己回帰モデルとは、時系列データにおけるある時刻k(m+1≦k≦M)のデータの予測値y^k(式において^はyの上に付けて表記)を、時系列データにおけるその時刻よりも前の時刻k−l(1≦l≦m)のデータの実績値yk-lを用いて表す式である。
Figure 0006919397
式1におけるαは、自己回帰モデルの係数である。また、mは、自己回帰モデルにおいてある時刻kにおける計測データyの値であるykを、その時刻よりも前の過去幾つのデータを用いて近似するかを示すM未満の整数であり、本実験では、1500とする。
続いて、最小二乗法を用いて、自己回帰モデルによる予測値y^kが実測値であるykに近似するための条件式を求める。自己回帰モデルによる予測値y^kが実測値であるykに近似するための条件として、式1で与えられるy^kとykとの二乗誤差を最小化することを考える。即ち、本実験では、自己回帰モデルによる予測値y^kをykに近似するために最小二乗法を用いる。以下の式2は、計測データと自己回帰モデルによる予測値との二乗誤差を最小にするための条件式である。
Figure 0006919397
式2より、以下の式3の関係を満たす。
Figure 0006919397
また、式3を変形(行列表記)することで、以下の式4のようになる。
Figure 0006919397
式4におけるRjlは計測データyの自己相関と呼ばれるもので、以下の式5で定義される値である。このときの|j−l|を時差という。
Figure 0006919397
式4を基に、以下の自己回帰モデルの係数に関する関係式である式6を考える。式6は、自己回帰モデルによる計測データの予測値と、その予測値に対応する時刻における計測データと、の誤差を最小化する条件から導出される方程式で、ユール・ウォーカー(Yule−Walker)方程式と呼ばれるものである。また、式6は自己回帰モデルの係数から成るベクトルを変数ベクトルとする線形方程式で、式6における左辺の定数ベクトルは、時差が1からmまでの計測データの自己相関を成分とするベクトルで、以下では、自己相関ベクトルとする。また、式6における右辺の係数行列は、時差が0からm−1までの計測データの自己相関を成分とする行列であり、以下では、自己相関行列とする。
Figure 0006919397
また、式6における右辺の自己相関行列(Rjlで構成されるm×mの行列)を、以下の式7のように、自己相関行列Rと表記する。
Figure 0006919397
一般に、自己回帰モデルの係数を求める際には、式6を係数αについて解くという方法が用いられる。式6では、自己回帰モデルで導出されるある時刻kにおける計測データの予測値y^kが、その時刻kにおける計測データの実績値ykにできるだけ近づくように係数αを導出する。よって、自己回帰モデルの周波数特性に、各時刻における計測データの実績値ykに含まれる多数の周波数成分が含まれる。したがって、例えば、計測データyに含まれるノイズが多い場合には、軸受106の振動に係る信号を抽出できなかったり、軸受106の故障の態様に応じた特徴を抽出することができなかったりするという問題が生じる。
そこで、本発明者らは、自己回帰モデルの係数αに乗算される自己相関行列Rに着目し、鋭意検討した結果、自己相関行列Rの固有値の一部を用いて、計測データyに含まれるノイズの影響が低減され、軸受の振動に係る信号成分が強調される(SN比を高める)ように自己相関行列Rを書き換えればよいという着想に至った。
以下で、このことの具体例を説明する。
自己相関行列Rを特異値分解する。自己相関行列Rの要素は、対称であるので、自己相関行列Rを特異値分解すると以下の式8のように、直交行列Uと、対角行列Σと、直交行列Uの転置行列との積となる。
Figure 0006919397
式8の行列Σは、式9に示すように、対角成分が自己相関行列Rの固有値となる対角行列である。対角行列Σの対角成分を、σ11、σ22、・・・、σmmとする。また、行列Uは、各列成分ベクトルが自己相関行列Rの固有ベクトルとなる直交行列である。直交行列Uの列成分ベクトルを、u1、u2、・・・、umとする。自己相関行列Rの固有ベクトルujに対する固有値がσjjという対応関係が有る。自己相関行列Rの固有値は、自己回帰モデルによる計測データの予測値の時間波形に含まれる各周波数の成分の強度に反映する変数である。
Figure 0006919397
自己相関行列Rの特異値分解の結果得られる対角行列Σの対角成分であるσ11、σ22、・・・、σmmの値は、数式の表記を簡略にするために降順とする。これらの自己相関行列Rの固有値のうち、最大のものから1以上且つm未満の設定された数である使用固有値数s個の固有値を用いて、以下の式10のように、行列R’を定義する。行列R’は、自己相関行列Rの固有値のうち使用固有値数s個の固有値用いて自己相関行列Rを近似した行列である。
Figure 0006919397
式10における行列Usは、式8の直交行列Uの左からs個の列成分ベクトル(使用される固有値に対応する固有ベクトル)により構成されるm×s行列である。つまり、行列Usは、直交行列Uから左のm×sの要素を切り出して構成される部分行列である。また、Us TはUsの転置行列であり、式8の行列UTの上からs個の行成分ベクトルにより構成されるs×m行列である。式10における行列Σsは、式8の対角行列Σの左からs個の列と上からs個の行により構成されるs×s行列である。つまり、行列Σsは、対角行列Σから左上のs×sの要素を切り出して構成される部分行列である。
行列Σs及び行列Usを行列要素表現すれば、以下の式11のようになる。
Figure 0006919397
自己相関行列Rの代わりに行列R’を用いることで、式6の関係式を、以下の式12のように書き換える。
Figure 0006919397
式12を変形することで、係数αを求める式13が得られる。式13によって求められた係数αを用いて、式1により予測値y^kを算出するモデルを「修正された自己回帰モデル」とする。
これまで対角行列Σの対角成分であるσ11、σ22、・・・、σmmの値を降順として説明したが、係数αの算出過程において対角行列Σの対角成分は降順である必要はなく、その場合には、行列Usは直交行列Uから左のm×sの要素を切り出すのではなく、使用される固有値に対応する列成分ベクトル(固有ベクトル)を切り出して構成される部分行列であり、行列Σsは対角行列Σから左上のs×sの要素を切り出すのではなく、使用される固有値を対角成分とするように切り出される部分行列である。
式13は、修正された自己回帰モデルの係数の決定に利用される方程式である。式13の行列Usは、自己相関行列Rの特異値分解により得られる直交行列の部分行列であって、利用される固有値に対応する固有ベクトルを列成分ベクトルとする行列である第3の行列である。また、式13の行列Σsは、自己相関行列Rの特異値分解により得られる対角行列の部分行列であって、利用される固有値を対角成分とする行列である第2の行列である。そして、式13の行列UsΣss Tは、行列Σsと行列Usとから導出される行列である第1の行列である。
Figure 0006919397
式13の右辺を計算することにより、修正された自己回帰モデルの係数αが求まる。以上に、修正された自己回帰モデルの係数の導出方法の一例について説明してきたが、その基となる自己回帰モデルの係数の導出については直感的に分かり易いように予測値に対して最小二乗法を用いる方法で説明した。しかしながら、一般的には確率過程という概念を用いて自己回帰モデルを定義し、その係数を導出する方法が知られている。その場合に、自己相関は確率過程(母集団)の自己相関で表現されており、この確率過程の自己相関は時差の関数として表されるものである。従って、本実施形態における計測データの自己相関は、確率過程の自己相関を近似するものであれば他の計算式で算出した値に代えても良く、例えば、R22〜Rmmは時差が0の自己相関であるが、これらをR11に置き換えても良い。
次に、求めた係数αで特定される修正された自己回帰モデルの周波数特性を求める方法を説明する。
式1で与えられる計測データの予測値y^kの予測誤差xkが白色雑音になるという性質を使うと、修正された自己回帰モデルは白色雑音である予測誤差xkを入力とし、計測データの実測値ykを出力とする線形時不変システムとみなせるので、以下の手順で周波数特性の算出式を導出することができる。以下では、修正された自己回帰モデルを線形時不変モデルとみなしたものを、単にシステムと称することにする。予測誤差xkは、以下の式14で表される。
Figure 0006919397
式14の両辺にz変換を施すと、以下の式15が得られる。
Figure 0006919397
式15から、システムのインパルス応答のz変換である伝達関数H(z)は、以下の式16のように求まる。
Figure 0006919397
システムの周波数特性は、正弦波入力に対する出力の振幅と位相との変化として現れ、インパルス応答のフーリエ変換で求められる。言い換えると、zが複素平面の単位円上を回転した時の伝達関数H(z)が周波数特性となる。ここで、式16におけるzを、以下の式17のように置くことを考える。
Figure 0006919397
ここで、jは虚数単位、ωは角周波数、Tは標本間隔である。
その場合、H(z)の振幅特性(システムの周波数特性)は、以下の式18のように表すことができる。式18におけるωTは、0〜2πの範囲で変化するものとする。
Figure 0006919397
本実験では、得られた計測データ毎に、式5と式7を用いて自己相関行列Rを求め、式8で表される特異値分解を行うことによって自己相関行列Rの固有値を求め、自己相関行列Rの固有値の分布を求めた。また、得られた計測データ毎に、計測データの波形を求めた。また、得られた計測データ毎に、上記手順で求められた自己相関行列Rの特異値分解から式11を用いて行列Σsと行列Usとを求め、式5と式13を用いて修正された自己回帰モデルの係数αを求め、求めた係数αで特定される修正された自己回帰モデルの波形を求めた。そして、得られた計測データ毎に、式18を用いて、修正された自己回帰モデルの周波数特性等を求めた。
求められた実験結果を、図3〜9を用いて説明する。
図3は、計測データそれぞれについて、上記手順で求められた自己相関行列Rの固有値の分布を示す図である。
図3(a)のグラフは、軸受106に疵がない正常な軸受を用いた場合に計測された計測データを利用した際の自己相関行列Rの固有値の分布を示すグラフである。図3(b)のグラフは、軸受106に内輪に疵がある軸受を用いた場合に計測された計測データを利用した際の自己相関行列Rの固有値の分布を示すグラフである。図3(c)のグラフは、軸受106に保持器に疵がある軸受を用いた場合に計測された計測データを利用した際の自己相関行列Rの固有値の分布を示すグラフである。
図3の各グラフは、自己相関行列Rを特異値分解して得られた固有値σ11〜σmmを昇順に並べ替え、プロットしたグラフであり、横軸が固有値のインデックス、縦軸が固有値の値を示す。
図3(a)を見ると、他よりも顕著に高い値をもつ固有値が5つあるのが分かる。その5つの固有値の中でも特に2つの固有値が他の3つよりも高い値となっていることが分かる。
何れの場合にも、自己相関行列Rの固有値のうち値が他よりも顕著に高い固有値の数は、固有値全体の数よりも顕著に少ないことが分かる。
図4は、計測データそれぞれの波形を示す図である。
図4(a)のグラフは、軸受106に疵がない正常な軸受を用いた場合に計測された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの実績値を示す。図4(b)のグラフは、軸受106に内輪に疵がある軸受を用いた場合に計測された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの実績値を示す。図4(c)のグラフは、軸受106に保持器に疵がある軸受を用いた場合に計測された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの実績値を示す。
図4の各グラフを見ると、例えば、図4(a)及び図4(c)が似通った波形をしているのが分かる。そのため、図4のグラフから軸受106に生じた異常を診断することは困難である。
図5は、使用固有値数sを1として求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示す図である。
図5(a)のグラフは、使用固有値数sを1として、軸受106に疵がない正常な軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図5(b)のグラフは、使用固有値数sを1として、軸受106に内輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図5(c)のグラフは、使用固有値数sを1として、軸受106に保持器に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。
図5の各グラフは、図4の各グラフと比べると、ノイズ成分が低減されていることが分かる。
図6は、使用固有値数sを1として求められた係数αで特定される修正された自己回帰モデルの周波数特性を示す図である。図6の各グラフは、図5の各波形についての周波数特性を示すグラフである。
図6(a)のグラフは、使用固有値数sを1として、軸受106に疵がない正常な軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式18で求められた周波数特性を示すグラフであり、横軸が周波数、縦軸が信号強度を示す。図6(b)のグラフは、使用固有値数sを1として、軸受106に内輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式18で求められた周波数特性を示すグラフであり、横軸が周波数、縦軸が信号強度を示す。図6(c)のグラフは、使用固有値数sを1として、軸受106に保持器に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式18で求められた周波数特性を示すグラフであり、横軸が周波数、縦軸が信号強度を示す。
図6(a)のグラフを見ると、周波数824Hzの部分にピークが立っているのが分かる。また、図6(b)のグラフを見ると、周波数1649Hzの部分にピークが立っているのが分かる。また、図6(c)のグラフを見ると、周波数1646Hzの部分にピークが立っているのが分かる。
図6の各グラフにおけるピークが立っている周波数を見ると、軸受106に疵がない場合は、軸受106に疵がある場合に比べて小さい周波数でピークが立っていることが分かる。即ち、使用固有値数sを1として、求められた係数αで特定される修正された自己回帰モデルの周波数特性から、軸受106における疵の有無を判断することが可能であることが分かる。
以上より、自己相関行列Rの1個の固有値に基づいて求められた係数αで特定される修正された自己回帰モデルの周波数特性に基づいて、軸受106の異常診断を行うことができるという知見を得ることができた。
しかし、図6(b)と図6(c)とのグラフに示されるように、軸受106の内輪に疵がある場合と、軸受106の保持器に疵がある場合と、では、似通った値の周波数において周波数特性のピークが立っている。そのため、使用固有値数sを1として、求められた係数αで特定される修正された自己回帰モデルの周波数特性から、軸受106に内輪の又は保持器の疵の有無を判断することは可能であるが、内輪の疵であるか、保持器の疵であるか、両方の疵であるか、を判断することは困難である。
そこで、本実験では、更に、使用固有値数sの値を5に増やして、自己相関行列Rの固有値のうち、最大のものから5つを用いて係数αを求め、求めた係数αで特定される修正された自己回帰モデルの波形、周波数特性を求めた。
図7は、使用固有値数sを5として求められた係数αで特定される修正された自己回帰モデルの波形を示す図である。
図7(a)のグラフは、使用固有値数sを5として、軸受106に疵がない正常な軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図7(b)のグラフは、使用固有値数sを5として、軸受106に内輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図7(c)のグラフは、使用固有値数sを5として、軸受106に保持器に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。
図8は、使用固有値数sを5として求められた係数αで特定される修正された自己回帰モデルの周波数特性を示す図である。図8の各グラフは、図7の各波形についての周波数特性を示すグラフである。
図8(a)のグラフは、使用固有値数sを5として、軸受106に疵がない正常な軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式18で求められた周波数特性を示すグラフであり、横軸が周波数、縦軸が信号強度を示す。図8(b)のグラフは、使用固有値数sを5として、軸受106に内輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式18で求められた周波数特性を示すグラフであり、横軸が周波数、縦軸が信号強度を示す。図8(3)のグラフは、使用固有値数sを5として、軸受106に保持器に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式18で求められた周波数特性を示すグラフであり、横軸が周波数、縦軸が信号強度を示す。
図8(a)のグラフを見ると、周波数821Hzと1047Hzの部分にピークが立っているのが分かる。また、図8(b)のグラフを見ると、周波数1648Hz、2021Hz、6474Hzの部分にピークが立っているのが分かる。また、図8(c)のグラフを見ると、周波数1646Hz、3291Hz、1746Hzの部分にピークが立っているのが分かる。
図8の各グラフにおけるピークと、図6の各グラフにおけるピークと、を見比べてみると、図6のグラフでピークが立っている周波数近辺においては、図8のグラフでもピークが立っているのが分かる。図8のグラフでは、更に、図6のグラフで見られなかったピークが見られるようになっている。
このように、係数αを求めるのに利用される自己相関行列Rの固有値を増やすことで、修正された自己回帰モデルによる予測値において、図6でピークが立っていた周波数以外の周波数の成分の信号強度が増加したことが分かる。即ち、係数αを求めるのに利用される自己相関行列Rの固有値は、求められた係数αで特定される修正された自己回帰モデルによる計測データの予測値の時間波形に含まれる各周波数の成分の強度と相関があることが見出された。
また、図8(b)、図8(c)のグラフを見てみると、図6(b)、図6(c)のグラフと同様に、1640Hz付近にピークが立っているのが分かる。しかし、図8(b)では、更に、2021Hz、6474Hzにおいてもピークが立っており、図8(c)では、更に、3291Hz、1746Hzにおいてもピークが立っていることが分かる。このように、図8では、内輪に疵のある軸受106と保持器に疵のある軸受106とで、ピークが立っている周波数に顕著な違いが見て取れる。即ち、使用固有値数sを5にすることで、軸受の疵が、内輪の疵であるか、保持器の疵であるか、両方の疵であるか、を見分けることができる。
このように本発明者らは、使用固有値数sが1の場合には見分けることが困難であった疵の存在する部位を、使用固有値数sを増やすことで、見分けることができるようになるという知見を得た。一方で、本発明者らは、使用固有値数sを増やし過ぎると修正された自己回帰モデルは、実際の計測データに近づいていき、軸受106の異常診断に有用でないノイズ成分等についても信号強度が増加されてしまい、式18で得られる周波数特性にもノイズ成分に起因するピークが生じるという知見を得た。即ち、自己相関行列Rの固有値のうち、相対的に値の高い一部の固有値を用いて係数αを求めることで、計測データに含まれる成分のうち、物体の異常診断に有用な成分を近似した修正された自己回帰モデルを求めることができると考えられる。
以上より、係数αを求めるのに利用される自己相関行列Rの固有値の数である使用固有値sを調整することで、より精度よく、軸受106の異常診断を行うことができるという知見を得ることができた。
本実施形態の処理は、以下のような処理である。即ち、本実験で得られた知見を基づいて、物体の振動の計測データから求められる自己相関行列Rを特異値分解し、得られる固有値の一部を用いて、計測データを近似する修正された自己回帰モデルの係数を求める。次に、求めた係数で特定される修正された自己回帰モデルから、式1の右辺の計算を行って計測データの予測値を求める訳ではなく、求めた係数を用いて、修正された自己回帰モデルの周波数特性を取得する。そして、取得した周波数特性に基づいて、物体の異常を診断する処理である。ただし、本実施形態の処理に、計測データの予測値(式1の左辺の値)を求める処理を加えてもよい。
(システム構成)
図9は、本実施形態の診断システムのシステム構成の一例を示す図である。診断システムは、周期的な運動を行う物体についての異常診断を行うシステムである。診断システムは、情報処理装置900、振動計測装置901を含む。本実施形態では、診断システムは、鉄道台車に利用されている軸受の異常診断を行う。
情報処理装置900は、振動計測装置901により計測された信号に基づいて、軸受の異常診断を行うパーソナルコンピュータ(PC)、サーバ装置、タブレット装置等の情報処理装置である。また、情報処理装置900は、電車に組み込まれたコンピュータ等であってもよい。
振動計測装置901は、加速度センサ等のセンサを含み、センサを介して物体の振動を検知し、検知した振動に応じた信号を有線又は無線で情報処理装置900等の外部の装置に出力する計測装置である。
(情報処理装置の機能構成)
図9を参照しながら、情報処理装置900が有する機能の一例を説明する。
情報処理装置900は、取得部910、決定部920、診断部930、判定部940、出力部950を含む。
≪取得部910≫
取得部910は、周期的な運動を行う物体について、この周期的な運動に係る計測データyを取得する。
本実施形態では、取得部910は、軸受106を回転させ、位置111に設置された振動計測装置901により計測された振動の計測データyを取得する。
≪決定部920≫
決定部920は、取得部910により取得された計測データyに基づいて、修正された自己回帰モデルにおける係数αを決定する。修正された自己回帰モデルは、計測データyの実績値yk-1〜yk-mと、この実績値に対する係数α(=α1〜αm)と、を用いて、計測データyの予測値y^kを表す式1である。決定部920は、修正された自己回帰モデルにおける係数αを、一般的に知られている自己回帰モデルにおける係数の決定に用いる式6(ユール・ウォーカー方程式)において、式7で表される自己相関行列Rの代わりに、自己相関行列Rから物体の異常診断に有用な成分を抽出した第1の行列R’(式10)を用いる式である式12を用いて決定する。式6は、時差が0からm−1までの計測データyの自己相関を成分とする自己相関行列Rを係数行列とし、時差が1からmまでの計測データyの自己相関を成分とする自己相関ベクトルを定数ベクトルとする方程式である。
式6は、式1で算出される計測データyの予測値y^kと、計測データyの予測値y^kに対応する時刻kにおける計測データyの実測値ykとの二乗誤差を最小化する条件を示す条件式として導出することができる。第1の行列R’は、自己相関行列Rを特異値分解(式8)することで導出される自己相関行列Rの固有値を対角成分とする対角行列Σと、自己相関行列Rの固有ベクトルを列成分とする直交行列Uと、から導出される。第1の行列R’は、自己相関行列Rの固有値のうち1以上且つm未満の設定された使用固有値数s個の固有値を用いて、対角行列Σの部分行列であって、使用固有値数s個の固有値を対角成分とする行列Σsと、直交行列Uの部分行列であって、使用固有値数s個の固有値に対応する固有ベクトルを列成分ベクトルとする行列Usと、から導出される行列UsΣss Tである。使用固有値数s個の固有値は、自己相関行列Rの固有値のうち、値が最大の固有値を含むようにすればよいが、値が大きいものから順に選ばれるのが好ましい。
後述する判定部940により、物体の異常の診断が成功していないと判定されると、決定部920は、使用固有値数sを設定し直し、修正された自己回帰モデルにおける係数αを決定し直す。
≪診断部930≫
診断部930は、決定部920により決定された係数αに基づいて、物体の異常を診断する。診断部930は、修正された自己回帰モデルの周波数の分布を示す周波数特性を、決定部920により決定された係数αを用いて式18により導出し、この周波数特性に基づいて、物体の異常を診断してもよい。更に、診断部930は、周波数特性におけるピークを示す周波数に基づいて、物体における異常の有無を診断の結果としてもよいし、物体における異常を有する部位を診断の結果としてもよい。
本実施形態では、軸受106における疵の有無、又は、疵を有する部位を診断の結果とする。
≪判定部940≫
判定部940は、診断部930による物体の異常の診断が成功したか否かを判定する。判定部940により物体の異常の診断が成功していないと判定されると、前述した決定部920は、使用固有値数sを設定し直し、修正された自己回帰モデルにおける係数αを決定し直す。
≪出力部950≫
出力部950は、診断部930による診断の結果に係る情報を出力する。
(情報処理装置のハードウェア構成)
図10は、情報処理装置900のハードウェア構成の一例を示す図である。
情報処理装置900は、CPU1000、主記憶装置1001、補助記憶装置1002、入出力I/F1003を含む。各構成要素は、システムバス1004を介して、相互に通信可能に接続されている。
CPU1000は、情報処理装置900を制御する中央演算装置である。主記憶装置1001は、CPU1000のワークエリアやデータの一時的な保管場所として機能するRAM(Random Access Memory)等の記憶装置である。補助記憶装置1002は、各種のプログラム、設定データ、振動計測装置901から出力される計測データ、診断情報等を記憶するROM(Read Only Memory)、HDD(Hard Disk Drive)、SSD(Solid State Drive)等の記憶装置である。入出力I/F1003は、振動計測装置901等の外部の装置との間での情報のやり取りに利用されるインターフェースである。
CPU1000が、補助記憶装置1002等に記憶されたプログラムに基づき処理を実行することによって、図9で説明した情報処理装置900の機能及び図11で後述するフローチャートの処理等が実現される。
(異常診断処理)
本実施形態では、診断システムは、図1と同様の状況で位置111に設置された振動計測装置901により計測された振動の計測データに基づいて、軸受106の異常を診断することとする。本実施形態では、診断システムは、軸受106のどこに疵が存在するかを特定することで、軸受106の異常を診断することとするが、軸受106に疵が存在するか否かを特定することで、軸受106の異常を診断することとしてもよい。本実施形態では、診断システムは、予め振動計測装置901により測定された測定データに基づいて、軸受106の異常を診断することとする。
図11は、情報処理装置900の処理の一例を示すフローチャートである。
S1101において、取得部910は、振動計測装置901により計測された振動の計測データyを取得する。取得部910は、計測データyとして、時刻1〜Mまでの計測データであるy1〜yMを取得する。
S1102において、決定部920は、S1101で取得した計測データyと、予め設定された定数Mと、修正された自己回帰モデルにおいて、ある時刻のデータを過去の幾つのデータで近似するかを示す数mと、に基づいて、式5と式7とを用いて自己相関行列Rを生成する。決定部920は、予め記憶されているM、mの情報を読み出すことで、Mとmとを取得する。本実施形態では、mの値は、1500である。また、Mは、mよりも大きい整数である。
S1103において、決定部920は、S1102で生成した自己相関行列Rを特異値分解することで、式8の直交行列Uと対角行列Σを取得し、対角行列Σから自己相関行列Rの固有値σ11〜σmmを取得する。
S1104において、決定部920は、修正された自己回帰モデルの係数αを求めるのに利用する自己相関行列Rの固有値の数である使用固有値数sを1にする。
S1105において、決定部920は、自己相関行列Rの複数の固有値であるσ11〜σmmのうち、最大のものから使用固有値数sだけの固有値σ11〜σssを、修正された自己回帰モデルの係数αを求めるのに利用する自己相関行列Rの固有値として選択する。そして、決定部920は、計測データyと、固有値σ11〜σssと、自己相関行列Rの特異値分解により得られた直交行列Uと、に基づいて、式13を用いて、修正された自己回帰モデルの係数αを決定する。
S1106において、診断部930は、S1105で決定した係数αに基づいて、式18を用いて、S1105で決定した係数αで特定される修正された自己回帰モデルの周波数特性を取得する。
S1107において、診断部930は、S1106で取得した周波数特性に基づいて、軸受106の異常を診断する。
本実施形態では、予め、周波数特性におけるピークが立っている周波数と、軸受106の状態と、の対応を示す情報であり、軸受106の異常診断に用いられる情報である診断情報が、補助記憶装置に記憶されている。診断情報とは、例えば、周波数特性におけるピークが立つ周波数が設定された範囲であれば、軸受106が正常であり、周波数特性におけるピークが立つ周波数が設定された範囲外であれば、軸受106に疵が存在することを示す情報等である。また、診断情報は、例えば、軸受106の特定の部位に疵が存在する状態に対応する周波数特性におけるピークが立つ周波数の範囲を示す情報であってもよい。また、診断情報は、例えば、設定された閾値以上の周波数においてピークが立っていれば、疵が存在する状態であり、また、その閾値未満の周波数においてピークが立っていれば、正常な状態である場合、その閾値の情報であってもよい。
また、図6、図8に示すように、軸受106等の診断対象の物体が同じ状態であったとしても、使用固有値数sの値によって、修正された自己回帰モデルの周波数特性におけるピークが立つ場所が異なる。そのため、補助記憶装置1002は、予め、使用固有値数に応じた診断情報を記憶している。
診断部930は、使用固有値数sの値に応じた診断情報を補助記憶装置1002から取得する。本実施形態では、診断部930は、診断情報として、軸受106の状態と、その状態に応じた周波数特性におけるピークが立つ周波数の範囲と、の対応情報を取得する。
診断部930は、S1106で取得した周波数特性に基づいて、ピークが立っている周波数を特定する。診断部930は、例えば、S1106で取得した周波数特性から、信号強度が設定された閾値以上となる周波数を、ピークが立っている周波数として特定する。CPU1000は、例えば、ピークが立っていると特定した周波数が連続する場合(例えば、100Hz、101Hz、102Hzでピークが立っている場合)、それらの連続する周波数のうち、信号強度が最も高いものをピークが立っている周波数として特定してもよい。
診断部930は、特定したピークが立っている周波数と、取得した診断情報と、に基づいて、軸受106の異常を診断する。診断部930は、例えば、診断情報が示す軸受106の状態と、周波数の範囲と、の対応に基づいて、特定した周波数が、軸受106のどの状態に対応するかを特定する。
S1108において、判定部940は、S1107での診断結果に基づいて、異常診断が成功したか否かを判定する。判定部940は、例えば、S1107で、特定の部位に疵が存在することを示す状態であること、又は正常なことを示す状態であることを特定した場合、異常診断が成功したとして、S1112の処理に進む。判定部940は、例えば、S1107で、部位が特定できない疵が存在することを示す状態であることを特定した、又は疵が存在するか否かを特定できなかった場合、異常診断が失敗したとして、S1109の処理に進む。
S1109において、決定部920は、使用固有値数sの値を変更する。CPU1000は、例えば、現在の使用固有値数sの値に1を加えた値を、新たな使用固有値数sの値とする。
S1110において、決定部920は、使用固有値数sの数が設定された閾値よりも大きいか否かを判定する。決定部920は、使用固有値数sの数が設定された閾値よりも大きいと判定した場合、S1111の処理に進む。また、決定部920は、使用固有値数sの数が設定された閾値以下であると判定した場合、S1105の処理に進む。これにより、決定部920は、値が変更された使用固有値数sを用いて、再度、修正された自己回帰モデルの係数αを決定する。
S1111において、決定部920は、軸受106の異常の診断が失敗したとして、S1112の処理に進む。
S1112において、出力部950は、軸受106の異常の診断の結果を示す情報を出力する。出力部950は、例えば、表示装置に、軸受106の異常の診断の結果を示す情報を表示することで出力する。また、出力部950は、例えば、軸受106の異常の診断の結果を示す情報を、補助記憶装置1002に記憶することで出力することとしてもよい。また、出力部950は、例えば、外部のサーバ装置等の設定された送信先に、軸受106の異常の診断の結果を示す情報を送信することで出力することとしてもよい。
(まとめ)
以上、本実施形態では、診断システムは、物体の周期的な運動から計測された計測データから、自己相関行列Rを生成し、自己相関行列Rを特異値分解し、得られた固有値のうち、最大のものから設定された数を用いて、計測データを近似する修正された自己回帰モデルの係数αを決定した。そして、診断システムは、決定した係数αから、修正された自己回帰モデルの周波数特性を求め、求めた周波数特性に基づいて、軸受106の異常を診断した。診断システムは、自己相関行列Rの固有値のうち一部の固有値を用いることで、軸受106の異常診断に有用な成分が残り、有用でない成分が残らないように、計測データを近似する修正された自己回帰モデルの係数αを決定できる。これにより、診断システムは、計測データの軸受106の異常診断に有用な成分に基づいて、より精度よく異常診断を行うことができる。
また、本実施形態の処理により、診断システムは、予め異常診断に用いる信号がどのような周波数の信号であるか等の想定を要することなく、計測データから異常診断に有用な成分を、抽出できる。
<実施形態2>
実施形態1では、軸受106を構成する部品のうち、内輪、保持器に疵がある場合についてのみ、修正された自己回帰モデルの周波数特性を求める実験を行った(図6、図8を参照)。このため、S1101において取得部910が振動の計測データyを取得する際のサンプリング周波数が、図6、図8のグラフにおいてピークが立つ周波数よりも高い周波数であっても、転動体、外輪に疵がある場合の修正された自己回帰モデルの周波数特性のグラフにおいてピークが立つ周波数が、このサンプリング周波数よりも高い場合には、転動体、外輪の疵を精度よく検出することができない虞がある。
そこで、発明者らは、更に、軸受106として転動体に疵のある軸受、外輪に疵のある軸受それぞれを用いて、実施形態1で説明した軸受異常の診断実験と同様の実験を行った。実験結果について説明する。
図12は、図3と同様に、自己相関行列Rの固有値の分布を示す図である。
図12(a)のグラフは、軸受106に転動体に疵がある軸受を用いた場合に計測された計測データを利用した際の自己相関行列Rの固有値の分布を示すグラフである。図12(b)のグラフは、軸受106に外輪に疵がある軸受を用いた場合に計測された計測データを利用した際の自己相関行列Rの固有値の分布を示すグラフである。
図12の各グラフは、図3と同様、自己相関行列Rを特異値分解して得られた固有値σ11〜σmmを昇順に並べ替え、プロットしたグラフであり、横軸が固有値のインデックス、縦軸が固有値の値を示す。
図12(a)、(b)の何れの場合にも、自己相関行列Rの固有値のうち値が他よりも顕著に高い固有値の数は、固有値全体の数よりも顕著に少ないことが分かる。
図13は、図4と同様に、計測データそれぞれの波形を示す図である。
図13(a)のグラフは、軸受106に転動体に疵がある軸受を用いた場合に計測された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの実績値を示す。図13(b)のグラフは、軸受106に外輪に疵がある軸受を用いた場合に計測された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの実績値を示す。
図14は、図5と同様に、使用固有値数sを1として求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示す図である。
図14(a)のグラフは、使用固有値数sを1として、軸受106に転動体に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図14(b)のグラフは、使用固有値数sを1として、軸受106に外輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。
図14の各グラフは、図13の各グラフと比べると、ノイズ成分が低減されていることが分かる。
図15は、図6と同様に、使用固有値数sを1として求められた係数αで特定される修正された自己回帰モデルの周波数特性を示す図である。図15の各グラフは、図14の各波形についての周波数特性を示すグラフである。
図15(a)のグラフは、使用固有値数sを1として、軸受106に転動体に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式18で求められた周波数特性を示すグラフであり、横軸が周波数、縦軸が信号強度を示す。図15(b)のグラフは、使用固有値数sを1として、軸受106に外輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式18で求められた周波数特性を示すグラフであり、横軸が周波数、縦軸が信号強度を示す。
図15(a)のグラフを見ると、周波数19853Hzの部分にピークが立っているのが分かる。また、図15(b)のグラフを見ると、周波数23007Hzの部分にピークが立っているのが分かる。
図15の各グラフにおけるピークが立っている周波数を見ると、図15(a)よりも、図15(b)の方が高い周波数であることが分かる。また、図15の各グラフにおけるピークが立っている周波数は、図6の各グラフにおけるピークが立っている周波数よりも顕著に高いことが分かる。そのため、使用固有値数sを1として、求められた係数αで特定される修正された自己回帰モデルの周波数特性から、軸受106の転動体又は外輪の疵の有無を判断することが可能であることが分かる。
図16は、図7と同様に、使用固有値数sを5として求められた係数αで特定される修正された自己回帰モデルの波形を示す図である。
図16(a)のグラフは、使用固有値数sを5として、軸受106に転動体に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。図16(b)のグラフは、使用固有値数sを5として、軸受106に外輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルで近似された計測データの波形を示すグラフであり、横軸が時間、縦軸が対応する時間における計測データの予測値を示す。
図17は、図8と同様に、使用固有値数sを5として求められた係数αで特定される修正された自己回帰モデルの周波数特性を示す図である。図17の各グラフは、図16の各波形についての周波数特性を示すグラフである。
図17(a)のグラフは、使用固有値数sを5として、軸受106に転動体に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式18で求められた周波数特性を示すグラフであり、横軸が周波数、縦軸が信号強度を示す。図17(b)のグラフは、使用固有値数sを5として、軸受106に外輪に疵がある軸受を用いた場合に計測された計測データから求められた係数αで特定される修正された自己回帰モデルについて、式18で求められた周波数特性を示すグラフであり、横軸が周波数、縦軸が信号強度を示す。
図17(a)のグラフを見ると、周波数19853Hzの部分にピークが立っているのが分かる。また、図17(b)のグラフを見ると、周波数22785Hz、23020Hzの部分にピークが立っているのが分かる。
図6、8、15、17の各グラフにおけるピークを見てみると、軸受106に疵がある場合に生じるピークのうち、最も周波数の高いものは、外輪に疵がある場合に23000Hz付近に生じるピークであることが分かる。
このことから、発明者らは、周波数23000Hzのピークを検出可能なサンプリング周波数で、軸受106の振動を計測することで、軸受106の各部位に生じる疵を検出可能であるとの知見を得た。また、軸受106の振動の計測データから得られる修正された自己回帰モデルの周波数特性でピークを示す周波数は、軸受106の内輪の回転数に比例するとみられる。修正された自己回帰モデルの周波数特性でピークを示す周波数のうち最大の周波数が23000Hzであり、軸受106の内輪の回転数は、2954.7rpmである。そのため、発明者らは、軸受106には、内部の疵により、最大で、内輪の回転数の約7.8(≒23000/2954.7)倍の周波数の信号が生じることとなるとの知見を得た。即ち、軸受106の内輪の回転数(rpm)の7.8倍の周波数(Hz)のピークを検出可能なサンプリング周波数で、軸受106の内輪の回転運動に係る振動を計測することで、軸受106の各部位に生じる疵を検出できる。
そこで、本実施形態では、内輪が回転する軸受106から、軸受106の内輪の回転数(rpm)の7.8倍の周波数(Hz)の信号を検出可能なサンプリング周波数で、測定された軸受106の内輪の回転運動に係る振動に係るデータに対して、実施形態1と同様の処理を行うこととする。
本実施形態の診断システムのシステム構成は、実施形態1と同様である。また、情報処理装置900のハードウェア構成及び機能構成も、実施形態1と同様である。
本実施形態の診断システムの処理のうち、実施形態1と異なる点について説明する。
振動計測装置901は、本実施形態では、軸受106の内輪の回転数(rpm)の7.8倍の周波数(Hz)の信号を検出可能なサンプリング周波数で、軸受106の振動を計測する。振動計測装置901は、例えば、ナイキスト周波数(Hz)が軸受106の内輪の回転数(rpm)の7.8倍となるように、軸受106の内輪の回転数(rpm)の15.6倍のサンプリング周波数(Hz)で軸受106の内輪の回転運動に係る振動を計測する。また、振動計測装置901は、例えば、軸受106の内輪の回転数(rpm)の15.6倍以上のサンプリング周波数(Hz)で軸受106の内輪の回転運動に係る振動を計測することとしてもよい。
S1101において、取得部910は、振動計測装置901により軸受106の内輪の回転数(rpm)の7.8倍の周波数(Hz)の信号を検出可能なサンプリング周波数で計測された振動の計測データyを取得する。
以上、本実施形態では、診断システムは、軸受106の内輪の回転数(rpm)の7.8倍の周波数(Hz)の信号を検出可能なサンプリング周波数で測定された測定データに基づいて、軸受106の異常診断を行うこととした。これにより、診断システムは、軸受106の各部位に生じる疵を検出できる。
<変形例>
実施形態1、2では、診断システムは、鉄道台車の軸受106の異常診断を行うこととした。しかし、診断システムは、歯車等の回転体、振動子等の振動体、バネ等の伸縮体、生物の心臓等の周期的な運動を行う他の物体についての異常診断を行うこととしてもよい。
実施形態1、2では、mは、予め設定された1500という数であるとした。しかし、診断対象に応じて実験を行う等して適切に計測データを近似できる過去のデータの数を特定し、特定した数をmの値としてもよい。例えば、情報処理装置900は、診断対象から計測された計測データについて、ある時刻の計測データを、その時刻よりも過去の計測データを複数用いて自己回帰モデル等を利用して近似して、近似した値とその時刻の計測データとの差分が設定された閾値未満となる場合の近似に用いられた過去の計測データの数を、mとして決定してもよい。
実施形態1、2では、情報処理装置900は、使用固有値数sの初期値を1として、異常診断に失敗すると使用固有値数sの値を変更することとした。しかし、情報処理装置900は、以下のようにして、使用固有値数sの値を決定してもよい。即ち、情報処理装置900は、まず、自己相関行列Rの固有値全ての合計値を求める。そして、情報処理装置900は、自己相関行列Rの固有値のうち、最大のものからa個の合計が、求めた合計値の設定された割合(例えば90%)以上となるaのうち、最小のものを使用固有値数sとして決定してもよい。
実施形態1、2では、情報処理装置900は、自己相関行列Rの固有値のうち、最大のものから、使用固有値数s個の固有値を用いて、物体の異常診断を行うこととした。しかし、情報処理装置900は、自己相関行列Rの固有値のうち、任意のs個の固有値を用いて、物体の異常診断を行うこととしてもよい。例えば、情報処理装置900は自己相関行列Rの固有値のうち、最大のものと、3つ目に大きいものと、の2つを利用して、異常診断することとしてもよい。
実施形態1、2では、診断システムは、予め振動計測装置901により測定された測定データに基づいて、軸受106の異常を診断することとした。しかし、診断システムは、リアルタイムで振動計測装置901により測定されて出力され続けている測定データを用いて、稼働中の軸受106の異常を診断することとしてもよい。即ち、情報処理装置900は振動計測装置901から新たな計測データyが取得される度に、計測データy1〜yMを更新して、図11の処理を行うようにしてもよい。
実施形態1、2では、情報処理装置900は、補助記憶装置1002に記憶されたプログラムを実行することで、図11の処理を実現することとした。しかし、情報処理装置900は、外付けの記憶媒体や外部の記憶サーバ等に記憶されたプログラムを実行することで、図11の処理を実現することとしてもよい。
また、例えば、上述した情報処理装置900の機能の一部又は全てをハードウェアとして情報処理装置900に実装してもよい。
また、以上説明した本発明の実施形態は、何れも本発明を実施するにあたっての具体化の例を示したものに過ぎず、これらによって本発明の技術的範囲が限定的に解釈されてはならないものである。即ち、本発明はその技術思想、またはその主要な特徴から逸脱することなく、様々な形で実施することができる。
900 情報処理装置
910 取得部
920 決定部
930 診断部
940 判定部
950 出力部
901 振動計測装置
1000 CPU

Claims (13)

  1. 周期的な運動を行う物体の前記周期的な運動に係る計測データを取得する取得手段と、
    前記取得手段により取得された前記計測データに基づいて、修正された自己回帰モデルにおける係数を決定する決定手段と、
    前記決定手段により決定された前記係数に基づいて、前記物体の異常を診断する診断手段と、
    を有し、
    前記修正された自己回帰モデルは、前記計測データの実績値と、前記実績値に対する前記係数と、を用いて、前記計測データの予測値を表す式であり、
    前記決定手段は、前記係数から成るベクトルを未知ベクトルとし、前記計測データから導出される自己相関行列を特異値分解することで導出される前記自己相関行列の固有値を対角成分とする対角行列と、前記自己相関行列の固有ベクトルを列成分とする直交行列と、から導出される第1の行列を係数行列とし、前記計測データから導出される自己相関ベクトルを定数ベクトルとする線形方程式を用いて、前記係数を決定し、
    前記自己相関ベクトルは、時差が1から前記修正された自己回帰モデルで用いられる前記実績値の数であるmまでの前記計測データの自己相関を成分とするベクトルであり、
    前記自己相関行列は、時差が0からm−1までの前記計測データの自己相関を成分とする行列であり、
    前記第1の行列は、1以上且つm未満の設定された数であるsに対して、前記自己相関行列のs個の固有値と前記対角行列とから導出される第2の行列Σsと、前記s個の固有値と前記直交行列とから導出される第3の行列Usと、から導出される行列UsΣss Tであり、
    前記第2の行列は、前記対角行列の部分行列であって、前記s個の固有値を対角成分とする行列であり、
    前記第3の行列は、前記直交行列の部分行列であって、前記s個の固有値に対応する固有ベクトルを列成分ベクトルとする行列であり、
    前記診断手段は、前記修正された自己回帰モデルの周波数の分布を示す周波数特性を、前記決定手段により決定された係数を用いて導出し、前記周波数特性に基づいて、前記物体の異常を診断する情報処理装置。
  2. 前記s個の固有値は、前記自己相関行列の固有値のうち、値が最大の固有値を含む請求項1記載の情報処理装置。
  3. 前記s個の固有値は、前記自己相関行列の固有値の中から、値が大きいものから順に選ばれた固有値である請求項2記載の情報処理装置。
  4. 前記診断手段による前記物体の異常の診断が成功したか否かを判定する判定手段を更に有し、
    前記決定手段は、前記判定手段により、前記診断手段による前記物体の異常の診断が成功していないと判定されると、sを設定し直して、前記係数を決定し直し、
    前記診断手段は、更に、前記判定手段により、前記診断手段による前記物体の異常の診断が成功していないと判定されると、前記決定手段により決定し直された係数に基づいて、前記周期的な運動を行う物体の異常を診断する請求項1乃至3何れか1項記載の情報処理装置。
  5. 前記自己相関行列は、前記計測データの予測値と当該計測データの予測値に対応する時刻における前記計測データの実測値との二乗誤差を最小化する条件を示す条件式に含まれる行列である請求項1乃至4何れか1項記載の情報処理装置。
  6. 前記診断手段は、前記周波数特性におけるピークを示す周波数に基づいて、前記物体における異常の有無を特定し、特定した結果を診断結果とすることで、前記物体の異常を診断する請求項1乃至5何れか1項記載の情報処理装置。
  7. 前記診断手段は、前記周波数特性においてピークを示す周波数に基づいて、前記物体における異常を有する部位を特定し、特定した結果を診断結果とすることで、前記物体の異常を診断する請求項1乃至6何れか1項記載の情報処理装置。
  8. 前記計測データは、前記物体の振動に係る計測データである請求項1乃至何れか1項記載の情報処理装置。
  9. 前記診断手段による診断の結果に係る情報を出力する出力手段を更に有する請求項1乃至何れか1項記載の情報処理装置。
  10. 前記物体は、鉄道台車に利用される軸受であり、
    前記診断手段は、前記決定手段により決定された前記係数に基づいて、前記物体における疵の存在を含む異常を診断する請求項1乃至何れか1項記載の情報処理装置。
  11. 前記周期的な運動は、前記軸受における内輪の回転運動であり、
    前記計測データは、前記軸受の内輪の回転数(rpm)の15.6倍以上のサンプリング周波数(Hz)で計測された前記軸受の振動に係る計測データである請求項1記載の情報処理装置。
  12. 情報処理装置が実行する情報処理方法であって、
    周期的な運動を行う物体の前記周期的な運動に係る計測データを取得する取得ステップと、
    前記取得ステップで取得された前記計測データに基づいて、修正された自己回帰モデルにおける係数を決定する決定ステップと、
    前記決定ステップで決定された前記係数に基づいて、前記物体の異常を診断する診断ステップと、
    を含み、
    前記修正された自己回帰モデルは、前記計測データの実績値と、前記実績値に対する前記係数と、を用いて、前記計測データの予測値を表す式であり、
    前記決定ステップでは、前記係数から成るベクトルを未知ベクトルとし、前記計測データから導出される自己相関行列を特異値分解することで導出される前記自己相関行列の固有値を対角成分とする対角行列と、前記自己相関行列の固有ベクトルを列成分とする直交行列と、から導出される第1の行列を係数行列とし、前記計測データから導出される自己相関ベクトルを定数ベクトルとする線形方程式を用いて、前記係数を決定し、
    前記自己相関ベクトルは、時差が1から前記修正された自己回帰モデルで用いられる前記実績値の数であるmまでの前記計測データの自己相関を成分とするベクトルであり、
    前記自己相関行列は、時差が0からm−1までの前記計測データの自己相関を成分とする行列であり、
    前記第1の行列は、1以上且つm未満の設定された数であるsに対して、前記自己相関行列のs個の固有値と前記対角行列とから導出される第2の行列Σsと、前記s個の固有値と前記直交行列とから導出される第3の行列Usと、から導出される行列UsΣss Tであり、
    前記第2の行列は、前記対角行列の部分行列であって、前記s個の固有値を対角成分とする行列であり、
    前記第3の行列は、前記直交行列の部分行列であって、前記s個の固有値に対応する固
    有ベクトルを列成分ベクトルとする行列であり、
    前記診断ステップでは、前記修正された自己回帰モデルの周波数の分布を示す周波数特性を、前記決定ステップで決定された係数を用いて導出し、前記周波数特性に基づいて、前記物体の異常を診断する情報処理方法。
  13. コンピュータを、請求項1乃至1何れか1項記載の情報処理装置の各手段として機能させるためのプログラム。
JP2017151727A 2017-03-24 2017-08-04 情報処理装置、情報処理方法及びプログラム Active JP6919397B2 (ja)

Priority Applications (3)

Application Number Priority Date Filing Date Title
PCT/JP2018/028942 WO2019026980A1 (ja) 2017-08-04 2018-08-01 情報処理装置、情報処理方法及びプログラム
EP18840775.3A EP3663741A4 (en) 2017-08-04 2018-08-01 INFORMATION PROCESSING DEVICE, INFORMATION PROCESSING METHOD AND PROGRAM
CN201880052067.6A CN110998272A (zh) 2017-08-04 2018-08-01 信息处理装置、信息处理方法及程序

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2017058936 2017-03-24
JP2017058936 2017-03-24

Publications (2)

Publication Number Publication Date
JP2018163135A JP2018163135A (ja) 2018-10-18
JP6919397B2 true JP6919397B2 (ja) 2021-08-18

Family

ID=63859924

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2017151727A Active JP6919397B2 (ja) 2017-03-24 2017-08-04 情報処理装置、情報処理方法及びプログラム

Country Status (1)

Country Link
JP (1) JP6919397B2 (ja)

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP7036233B2 (ja) * 2019-01-09 2022-03-15 日本製鉄株式会社 情報処理装置、情報処理方法及びプログラム
US11334414B2 (en) 2019-03-04 2022-05-17 Mitsubishi Heavy Industries, Ltd. Abnormality detecting apparatus, rotating machine, abnormality detection method, and non- transitory computer readable medium
CN110057588B (zh) * 2019-05-09 2020-07-03 山东大学 基于奇异值与图论特征融合的轴承早期故障检测与诊断方法及系统
CN110926811A (zh) * 2019-12-30 2020-03-27 黄石邦柯科技股份有限公司 一种轴承检测系统
CN111721528B (zh) * 2020-05-18 2022-04-05 浙江工业大学 基于cms系统大数据的风力发电机组齿轮箱故障预警方法
CN112474395B (zh) * 2020-12-16 2022-04-15 大连贝林轴承仪器有限公司 一种机械手式全自动三代轮毂轴承振动测量仪及测量方法

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2517300B2 (ja) * 1987-07-21 1996-07-24 コニカ株式会社 高感度かつ生保存性の改良されたハロゲン化銀写真感光材料
JP3630312B2 (ja) * 2002-06-27 2005-03-16 日本ビクター株式会社 ディスク再生方法及びオーディオ信号の記録再生方法
BR0203292A (pt) * 2002-08-14 2004-05-18 Michitoshi Oishi Detetor de defeitos das máquinas rotativas
JP5535954B2 (ja) * 2011-01-28 2014-07-02 三菱重工業株式会社 健全性評価装置及びその方法並びにプログラム

Also Published As

Publication number Publication date
JP2018163135A (ja) 2018-10-18

Similar Documents

Publication Publication Date Title
JP6919397B2 (ja) 情報処理装置、情報処理方法及びプログラム
JP6848813B2 (ja) 情報処理装置、情報処理方法及びプログラム
US10520397B2 (en) Methods and apparatuses for defect diagnosis in a mechanical system
KR101748559B1 (ko) 회전체 진단 장치 및 방법
JP6237774B2 (ja) 情報処理システム、情報処理方法及びプログラム
US11209487B2 (en) Rotor diagnostic apparatus, rotor diagnostic method, and rotor diagnostic program
KR20100094452A (ko) 롤링 베어링 대미지 검측 및 자동 식별 방법
JP6835251B2 (ja) 損傷診断装置、損傷診断方法、及び、損傷診断プログラム
JP7470849B2 (ja) 異常検出装置、異常検出方法、及びプログラム
WO2013180727A1 (en) A method and a system for testing operational integrity of a drilling rig
Jana et al. Computer vision‐based real‐time cable tension estimation algorithm using complexity pursuit from video and its application in Fred‐Hartman cable‐stayed bridge
JPH10258974A (ja) 非定常信号解析装置及び非定常信号解析プログラムを記録した媒体
JP7014080B2 (ja) 情報処理装置、情報処理方法及びプログラム
JP6904418B2 (ja) 情報処理装置、情報処理システム、情報処理方法、及び、プログラム
WO1998011417A1 (fr) Analyseur de signal instable et support d'enregistrement de programme d'analyse de signal instable
KR20150090950A (ko) 진동 장치의 상태 판단 방법
JP2020123229A (ja) 異常予兆検出システム、異常予兆検出方法
TW201633025A (zh) 工具機主軸故障形式的診斷方法及其系統
JP2015161506A (ja) 信号解析装置、信号解析方法、劣化診断装置及び劣化診断方法
JP7036233B2 (ja) 情報処理装置、情報処理方法及びプログラム
JP2020144111A (ja) 異常検出装置、回転機械、異常検出方法、及びプログラム
JP7099296B2 (ja) 評価方法、システム構築方法、及び評価システム
Bapir et al. A comparative analysis of 1D convolutional neural networks for bearing fault diagnosis
WO2024116637A1 (ja) 分類装置、学習モデル生成装置、分類方法、および学習モデル生成方法
US20220163425A1 (en) Method for acquiring data for detecting damage to a bearing

Legal Events

Date Code Title Description
A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20200513

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20210202

A521 Written amendment

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20210301

TRDD Decision of grant or rejection written
A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20210622

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20210705

R151 Written notification of patent or utility model registration

Ref document number: 6919397

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R151