JP7184191B2 - 推定装置、推定方法、およびプログラム - Google Patents

Description

本発明は、推定装置、推定方法、およびプログラムに関し、特に、鉄道車両の状態を推定するために用いて好適なものである。本願は、２０１９年６月２８日に日本に出願された特願２０１９－１２１７５１号に基づき優先権を主張し、特願２０１９－１２１７５１号の内容を全てここに援用する。

鉄道車両の台車枠に亀裂が発生する等、台車枠の状態が異常になると、鉄道車両の走行に支障をきたす虞がある。そこで、台車枠の状態を推定する技術が求められる。この種の技術として特許文献１に記載の技術がある。
特許文献１では、台車の負荷箇所に荷重が与えられた際の車体の加速度検出箇所における単位荷重当たりの加速度の周波数特性（周波数第一特性）と、前記負荷箇所に荷重が与えられた際の台車の検査箇所における単位荷重当たりの応力の周波数特性（周波数第二特性）と、を記憶する。そして、鉄道車両の走行中に前記加速度検出箇所において検出した加速度の周波数特性（周波数第三特性）と、前記周波数第一特性と、に基づいて、走行中に台車に与えられた荷重を算出する。このようにして算出した走行中に台車に与えられた荷重と、前記周波数第二特性と、に基づいて、前記台車の検査箇所における車両の走行中の応力を算出する。そして、時間経過に応じた前記車両の走行中の応力の解析結果に基づいて、前記検査箇所の監視を行う。

特開２０１８－１５５５１７号公報

鷲津久一郎 他４名 編、「有限要素法ハンドブック II 応用編」、初版、１９８３年１月２５日、株式会社培風館、ｐ．６５－６７

しかしながら、特許文献１に記載の技術では、台車の検査箇所は、疲労の発生によって損傷し易い箇所として予め決められている。従って、予め決められた箇所における応力しか導出することができない。また、特許文献１に記載の技術では、台車の検査箇所における応力を、歪みゲージ等を用いて測定する必要がある。このため、台車の検査箇所を増やすことは容易ではない。従って、特許文献１に記載の技術では、鉄道車両が走行しているときの台車枠の状態を正確に推定することが容易ではない。

本発明は、以上のような問題点に鑑みてなされたものであり、鉄道車両が走行しているときの台車枠の状態を正確に推定することができるようにすることを目的とする。

本発明の推定装置は、台車枠と、前記台車枠に直接または他の部品を介して接続された接続部品と、を有する鉄道車両の状態を推定する推定装置であって、前記接続部品に取り付けられた第１のセンサで、前記鉄道車両の走行中に測定された第１の測定データを取得するデータ取得手段と、前記第１の測定データに基づいて、前記台車枠の着力箇所に作用する外力を導出し、前記台車枠の着力箇所に作用する外力に基づいて、前記台車枠に作用する外力の分布を導出する外力導出手段と、前記台車枠に作用する外力の分布に基づいて、前記台車枠の状態を表す情報を導出する状態導出手段と、を有し、前記第１のセンサは、前記台車枠の着力箇所に作用する外力を導出することが可能な物理量を測定するためのセンサであり、前記台車枠の着力箇所は、前記台車枠において、前記接続部品により発生する外力が作用する箇所であり、前記状態導出手段は、前記台車枠の振動を表す運動方程式を解くことにより、前記台車枠の変位の分布を前記台車枠の状態を表す情報として導出する変位分布導出手段を少なくとも有し、前記台車枠の振動を表す運動方程式における外力の分布には、前記外力導出手段により導出された前記台車枠に作用する外力の分布が与えられることを特徴とする。

本発明の推定方法は、台車枠と、前記台車枠に直接または他の部品を介して接続された接続部品と、を有する鉄道車両の状態を推定する推定方法であって、前記接続部品に取り付けられた第１のセンサで前記鉄道車両の走行中に測定された第１の測定データを取得するデータ取得工程と、前記第１の測定データに基づいて、前記台車枠の着力箇所に作用する外力を導出し、前記台車枠の着力箇所に作用する外力に基づいて、前記台車枠に作用する外力の分布を導出する外力導出工程と、前記台車枠に作用する外力の分布に基づいて、前記台車枠の状態を表す情報を導出する状態導出工程と、を有し、前記第１のセンサは、前記台車枠の着力箇所に作用する外力を導出することが可能な物理量を測定するためのセンサであり、前記台車枠の着力箇所は、前記台車枠において、前記接続部品により発生する外力が作用する箇所であり、前記状態導出工程は、前記台車枠の振動を表す運動方程式を解くことにより、前記台車枠の変位の分布を前記台車枠の状態を表す情報として導出する変位分布導出工程を少なくとも有し、前記台車枠の振動を表す運動方程式における外力の分布には、前記外力導出工程により導出された前記台車枠に作用する外力の分布が与えられることを特徴とする。

本発明のプログラムは、台車枠と、前記台車枠に直接または他の部品を介して接続された接続部品と、を有する鉄道車両の状態を推定するための処理をコンピュータに実行させるためのプログラムであって、前記接続部品に取り付けられた第１のセンサで前記鉄道車両の走行中に測定された第１の測定データを取得するデータ取得工程と、前記第１の測定データに基づいて、前記台車枠の着力箇所に作用する外力を導出し、前記台車枠の着力箇所に作用する外力に基づいて、前記台車枠に作用する外力の分布を導出する外力導出工程と、前記台車枠に作用する外力の分布に基づいて、前記台車枠の状態を表す情報を導出する状態導出工程と、をコンピュータに実行させ、前記第１のセンサは、前記台車枠の着力箇所に作用する外力を導出することが可能な物理量を測定するためのセンサであり、前記台車枠の着力箇所は、前記台車枠において、前記接続部品により発生する外力が作用する箇所であり、前記状態導出工程は、前記台車枠の振動を表す運動方程式を解くことにより、前記台車枠の変位の分布を前記台車枠の状態を表す情報として導出する変位分布導出工程を少なくとも有し、前記台車枠の振動を表す運動方程式における外力の分布には、前記外力導出工程により導出された前記台車枠に作用する外力の分布が与えられることを特徴とする。

図１は、鉄道車両の概略の一例を示す図である。 図２は、台車枠およびその周辺の部品の構成の第１の例を示す図である。 図３Ａは、結合要素の第１の例をモデル化して示す図である。 図３Ｂは、結合要素の第２の例をモデル化して示す図である。 図４は、推定装置の機能的な構成の第１の例を示す図である。 図５は、応力限界図の一例を概念的に示す図である。 図６は、推定方法の第１の例を説明するフローチャートである。 図７は、台車枠およびその周辺の部品の構成の第２の例を示す図である。 図８は、推定装置の機能的な構成の第２の例を示す図である。 図９は、推定装置の機能的な構成の第３の例を示す図である。 図１０は、推定方法の第２の例を説明するフローチャートである。 図１１は、推定装置のハードウェアの構成の一例を示す図である。 図１２は、台車枠の空気バネ座における変位と時間との関係の一例を示す図である。 図１３Ａは、比較的大きな応力が発生した台車枠の或る位置における変位と時間との関係の第１の例を示す図である。 図１３Ｂは、比較的大きな応力が発生した台車枠の或る位置における変位と時間との関係の第２の例を示す図である。 図１４は、モード行列の成分とモード行列の成分の近似値との関係の一例を示す図である。 図１５は、推定対象領域における最大主歪みと時間との関係の一例を示す図である。 図１６は、推定対象領域における最大主応力と時間との関係の一例を示す図である。 図１７は、推定対象領域における最大主応力の推定値および測定値の関係の一例を示す図である。

以下、図面を参照しながら、本発明の実施形態を説明する。
（第１の実施形態）
まず、第１の実施形態を説明する。
＜鉄道車両の概略構成＞
まず、本実施形態で例示する鉄道車両について説明する。図１は、鉄道車両の概略の一例を示す図である。図２は、台車枠およびその周辺の部品の構成の一例を示す図である。尚、図１、図２において、鉄道車両は、ｘ_１軸の正の方向に進むものとする（ｘ_１軸は、鉄道車両の走行方向に沿う軸である）。また、ｘ_３軸は、軌道２０（地面）に対し垂直方向（鉄道車両の高さ方向）であるものとする。ｘ_２軸は、鉄道車両の走行方向に対して垂直な水平方向（鉄道車両の走行方向と高さ方向との双方に垂直な方向）であるものとする。また、鉄道車両は、営業車両であるものとする。尚、各図において、○の中に●が付されているものは、紙面の奥側から手前側に向かう方向を示す。

図１、図２に示すように本実施形態では、鉄道車両は、車体１１と、台車１２ａ、１２ｂと、輪軸１３ａ～１３ｄと、を有する。このように本実施形態では、１つの車体１１に、２つの台車１２ａ、１２ｂと、４組の輪軸１３ａ～１３ｄと、が備わる鉄道車両を例に挙げて説明する。輪軸１３ａ～１３ｄは、車軸１５ａ～１５ｄと、その両端に設けられた車輪１４ａ～１４ｄと、を有する。本実施形態では、台車１２ａ、１２ｂが、ボルスタレス台車である場合を例に挙げて説明する。

図１では、表記の都合上、輪軸１３ａ～１３ｄの一方の車輪１４ａ～１４ｄのみを示す。輪軸１３ａ～１３ｄの他方にも車輪が配置されている（図１に示す例では、車輪は合計８個ある）。
図２において、輪軸１３ａ、１３ｂのｘ_２軸に沿う方向の両側には、軸箱１７ａ、１７ｂが配置される。軸箱１７ａ、１７ｂは、モノリンク１８ａ、１８ｂを介して台車枠１６と接続されている。また、軸箱１７ａ、１７ｂは、軸バネ１９ａ、１９ｂを介して台車枠１６と接続されている。尚、鉄道車両は、図１、図２に示す構成要素以外の構成要素を有する。表記および説明の都合上、図１、図２では、当該構成要素の図示を省略する。例えば、図２において、鉄道車両が軸ダンパを有するものである場合、軸箱１７ａ、１７ｂは、軸ダンパを介して台車枠１６と接続されている場合もある。

図２において、１つの台車１２ａに、１つの台車枠１６が配置される。軸箱１７ａ、１７ｂ、モノリンク１８ａ、１８ｂ、および軸バネ１９ａ、１９ｂは、１つの車輪に対して１つずつ配置される。前述したように１つの台車１２ａには、４つの車輪が配置される。従って、１つの台車１２ａには、軸箱、モノリンク、および軸バネが、それぞれ４つずつ配置される。ここで、台車枠１６と一体になっているものも台車枠１６に含まれるものとする。例えば、台車枠１６の本体に溶接されている排障器も、台車枠１６に含まれるものとする。尚、台車枠１６と、台車枠１６と一体となっているもととの境界は定まらない。また、台車枠１６と一体になっているものは、台車枠１６と同じ運動を行う。

図２では、台車１２ａにおける台車枠１６、軸箱１７ａ、１７ｂ、モノリンク１８ａ、１８ｂ、および軸バネ１９ａ、１９ｂのみを示す。台車１２ｂにおける台車枠、軸箱、モノリンク、および軸バネも図２に示すものと同じもので実現される。尚、鉄道車両自体は公知の技術で実現できるので、ここでは、その詳細な説明を省略する。
以下の説明では、台車枠１６と、軸箱１７ａ、１７ｂとを連結する部品を、必要に応じて結合要素と総称する。

結合要素は、台車枠１６と軸箱１７ａ、１７ｂとに接続される。従って、軸箱１７ａ、１７ｂの振動と、結合要素の振動とは連動する。これらの振動が結合要素を介して台車枠１６に伝搬する。台車枠１６に作用する外力は、結合要素の振動における粘性減衰力と剛性力との和で表される。粘性減衰力は、粘性減衰係数と速度との積で表される。剛性力は、剛性と変位との積で表される。以上のことから、本発明者らは、接続部品で測定される加速度（加速度から導出される速度および変位）に基づいて、台車枠１６に作用する外力を導出することにより、台車枠１６に作用する外力を正確に導出することができることを見出した。ここで、接続部品は、台車枠１６と一体となっていない部品であって、台車枠１６に直接または他の部品を介して接続される部品である。

このような接続部品として、例えば、結合要素を構成する部品と、結合要素に直接または部品を介して接続される部品との少なくとも一方の部品が挙げられる。より具体的に、本実施形態では、このような接続部品が、軸箱１７ａ、１７ｂである場合を例に挙げて説明する。そこで、本実施形態では、軸箱１７ａ、１７ｂのそれぞれに加速度センサ２１ａ、２１ｂを取り付ける。本実施形態では、加速度センサ２１ａ、２１ｂが、３次元加速度センサである場合を例に挙げて説明する。加速度センサ２１ａ、２１ｂで測定される加速度のデータから、加速度のｘ_１軸方向成分、ｘ_２軸方向成分、ｘ_３軸方向成分が得られる。尚、加速度センサは、軸箱１７ａ、１７ｂ以外の軸箱にも取り付けられる。本実施形態では、例えば、加速度センサ２１ａ、２１ｂが第１のセンサの一例である。また、加速度センサ２１ａ、２１ｂにより測定される加速度（軸箱１７ａ、１７ｂの加速度）が台車枠の着力箇所に作用する外力を導出することが可能な物理量の一例である。

以上のようにして台車枠１６に作用する外力が得られれば、台車枠１６の振動（運動）を表す運動方程式を解くことにより、台車枠１６の変位分布（台車枠１６の各部の位置の変化量）を導出することができる。台車枠１６の変位分布が得られれば、台車枠１６の応力分布が得られる。また、本実施形態では、計算の負荷を軽減するため、モード解析法を用いて台車枠１６の振動を表す運動方程式を解く。以下に、台車枠１６の応力分布を導出する方法の一例を説明する。

＜台車枠１６の応力分布の導出＞
＜＜台車枠１６の運動方程式＞＞
台車枠１６の振動を表す運動方程式は、以下の（１）式で表される。

ここで、［Ｍ］（∈Ｒ^{３ｌ×３ｌ}）は、台車枠１６の質量行列である。［Ｃ］（∈Ｒ^{３ｌ×３ｌ}）は、台車枠１６の粘性行列である。粘性行列は、減衰行列ともいう。［Ｋ］（∈Ｒ^{３ｌ×３ｌ}）は、台車枠１６の剛性行列である。尚、Ｒは、実数全体の集合を表す（このことは以降の説明でも同じである）
｛ｕ｝（∈Ｒ^３ｌ）は、台車枠１６の変位ベクトルである。｛ｆ｝（∈Ｒ^３ｌ）は、台車枠１６の外力ベクトルである。

本実施形態では、台車枠１６の振動には、ｘ_１軸方向、ｘ_２軸方向、およびｘ_３軸方向の３つの成分があるものとする。従って、｛ｕ｝を構成する各節点の変位および｛ｆ｝を構成する各節点の外力は、３つの自由度を有する。尚、（１）式において、・は、ｄ／ｄｔ（時間の一階微分）を表し、・・は、ｄ^２／ｄｔ^２（時間の二階微分）を表す（このことは、以降の式でも同じである）。

ｌは、数値解析における変位分布の近似解の自由度に対応する。本実施形態では、数値解析として、有限要素法（ＦＥＭ；Finite Element Method）を用いる場合を例に挙げて説明する。従って、ｌは、例えば、有限要素法のメッシュの節点の数である。この場合、台車枠１６の質量行列［Ｍ］、台車枠１６の粘性行列［Ｃ］、および台車枠１６の剛性行列［Ｋ］の成分には、それぞれ、各成分に対応する密度から導出される値、各成分に対応する粘性減衰係数から導出される値、および各成分に対応する剛性から導出される値が与えられる。密度、粘性減衰係数および剛性は、位置によらずに同じ値としても異なる値としてもよい。台車枠１６の質量行列［Ｍ］、粘性行列［Ｃ］、および剛性行列［Ｋ］の成分は、例えば、有限要素法による数値解析を行う公知のソルバーにおいて、有限要素法のメッシュと、台車枠１６全体の密度、粘性減衰係数、および剛性と、を用いて導出される。

（１）式の左辺第１項は、台車枠１６に作用する重力を表す慣性項である。（１）式の左辺第２項は、台車枠１６に作用する粘性力を表す減衰項である。（１）式の左辺第３項は、台車枠１６に作用する剛性力を表す剛性項である。

＜＜台車枠１６に作用する外力＞＞
（１）式の右辺の外力ベクトル｛ｆ｝は、結合要素に作用する外力を導出することによって与えられる。
そこで、図３Ａおよび図３Ｂを参照しながら、結合要素に作用する外力の導出方法の一例について説明する。図３Ａおよび図３Ｂは、結合要素の一例をモデル化して示す図である。図３Ａは、台車枠１６と軸箱１７ａとに接続される結合要素をモデル化した図を示す。図３Ｂは、台車枠１６と軸箱１７ｂとに接続される結合要素をモデル化した図を示す。台車枠１６とその他の軸箱（軸箱１７ｂ等）とに接続される結合要素をモデル化した図も、図３Ａおよび図３Ｂと同じようにして表されるので、ここでは、その詳細な説明を省略する。台車枠１６と結合要素との接続箇所は、台車枠１６と結合要素とが相互に接触する領域（の全体）としても、台車枠１６と結合要素とが接触する領域の代表点（例えば、重心の位置）としてもよい。図３Ａおよび図３Ｂに示す例では、説明を簡単にするため、台車枠１６と結合要素との接続箇所が点であるものとする。以下の説明では、台車枠１６と結合要素との接続箇所を、必要に応じて、着力点と称する。図３Ａおよび図３Ｂに示す例では、バネおよびダンパを並列に接続したモデルで、モノリンク１８ａおよび軸バネ１９ａを表す。モデル化したモノリンク１８ａと台車枠１６とは着力点３１ａで接続される。モデル化した軸バネ１９ａと台車枠１６とは着力点３２ａで接続される。また、軸箱１７ｂについても同様に、モデル化したモノリンク１８ｂと台車枠１６とは着力点３１ｂで接続され、モデル化した軸バネ１９ｂと台車枠１６とは着力点３２ｂで接続される。

結合要素の振動を表す運動方程式は、以下の（２）式で表される。

ここで、［Ｃ_ｂｃ］は、結合要素の粘性行列である。［Ｋ_ｂｃ］は、結合要素の剛性行列である。｛ｕ_０｝は、軸箱１７ａとの接続箇所における結合要素の変位と、着力点における結合要素の変位とで構成される結合要素の変位ベクトルである。｛ｆ_０｝は、結合要素の軸箱１７ａとの接続箇所に作用する外力と、結合要素の着力点に作用する外力とで構成される結合要素の外力ベクトルである。結合要素の粘性行列［Ｃ_ｂｃ］、および結合要素の剛性行列［Ｋ_ｂｃ］の成分には、それぞれ、各成分に対応する粘性減衰係数から導出される値、および各成分に対応する剛性から導出される値が与えられる。結合要素の粘性行列［Ｃ］、および剛性行列［Ｋ］の成分は、例えば、結合要素の軸箱１７ａとの接続箇所の位置と、結合要素の着力点の位置と、結合要素の粘性減衰係数と、結合要素の剛性と、を用いて導出される。

（２）式において、結合要素の着力点におけるｘ_１軸方向、ｘ_２軸方向、およびｘ_３軸方向の変位を０（ゼロ）とする。また、軸箱１７ａとの接続箇所における結合要素の変位を軸箱１７ａの変位とする。以上のようにすることによって、結合要素の外力ベクトル｛ｆ_０｝を構成する外力として、結合要素の着力点に作用する外力を導出することができる。つまり、モノリンク１８ａの着力点３１ａおよび軸バネ１９ａの着力点３２ａに作用する外力は、モノリンク１８ａおよび軸バネ１９ａのそれぞれについて構成した（２）式を用いて導出される。また、軸箱１７ｂについても同様に、軸箱１７ｂの変位から、モノリンク１８ｂの着力点３１ｂおよび軸バネ１９ｂの着力点３２ｂに作用する外力が導出される。そして、台車枠１６の着力点３１ａ、３１ｂおよび着力点３２ａ、３２ｂに作用する外力は、それぞれ、モノリンク１８ａ、１８ｂの着力点３１ａ、３１ｂおよび軸バネ１９ａ、１９ａの着力点３２ａ、３２ｂに作用する外力の反作用力として導出される。

＜＜変位ベクトル｛ｕ｝の導出＞＞
（１）式の外力ベクトル｛ｆ｝の（台車枠１６の）着力点３１ａ、３１ｂ、３２ａ、３２ｂに作用する外力の成分には、前述の方法で導出した値を与え、その他の成分には、０（ゼロ）を与えて、台車枠１６の変位ベクトル｛ｕ｝を導出することにより、台車枠１６の変位分布を導出することができる。本実施形態では、計算時間を短縮するため、（１）式で表される物理座標（現実空間の位置を表す座標）系の運動方程式を、モード座標系の運動方程式で表現し、モード座標系の運動方程式を解く場合を例に挙げて説明する。以下に、モード座標系の運動方程式を解く手法の一例を説明する。

（１）式の運動方程式に固有値解析を適用することにより、固有振動数ωおよび固有ベクトル｛φ｝を導出し、固有ベクトル｛φ｝＝｛φ^（１）｝、・・・、｛φ^（ｎ）｝で構成されるモード行列［φ］（∈Ｒ^３ｌ×ｎ）を導出する。尚、固有ベクトルは、固有モードベクトルとも称される。ここで、ｎ（∈Ｎ）は、モード数である。低い固有振動数から順にモード数ｎだけ固有振動数が選択される。固有値解析の手法は、例えば、非特許文献１に記載のモード解析法を用いることにより実現することができ、公知の技術であるので、ここでは、その詳細な説明を省略する。尚、固有値解析に際しては、非特許文献１に記載のように、（１）式の右辺を０（ゼロ）とおく。また、減衰項［Ｃ］｛ｕ・｝を０（ゼロ）とおいてもよい（ｕ・は、（１）式においてｕの上に・が付されているものに対応する）。また、固有ベクトル｛φ^（１）｝、・・・、｛φ^（ｎ）｝を基底とする変位ベクトルの座標系のことをモード座標系と称する。

モード座標系における台車枠１６の質量行列［Ｍ_ξ］（∈Ｒ^ｎ×ｎ）、モード座標系における台車枠１６の粘性行列［Ｃ_ξ］（∈Ｒ^ｎ×ｎ）、モード座標系における台車枠１６の剛性行列［Ｋ_ξ］（∈Ｒ^ｎ×ｎ）は、それぞれ、以下の（３）式、（４）式、（５）式で表される。

ここで、Ｔは、転置行列であることを示す（このことは、以降の式でも同じである）。また、モード座標系における台車枠１６の質量行列［Ｍ_ξ］、粘性行列［Ｃ_ξ］、剛性行列［Ｋ_ξ］は、それぞれ、以下の（６）式、（７）式、（８）式のように、対角行列である。尚、モード座標系における台車枠１６の質量行列［Ｍ_ξ］および粘性行列［Ｃ_ξ］において、対角成分以外の成分が０（ゼロ）でない場合には、当該成分は０（ゼロ）に近似されるものとする。

ここで、（１）、・・・、（ｎ）は、それぞれ、一次の固有振動モード、・・・、ｎ次の固有振動モードに対応する成分であることを示す（このことは、以降の式でも同じである）。
モード座標系における台車枠１６の質量行列［Ｍ_ξ］、粘性行列［Ｃ_ξ］、剛性行列［Ｋ_ξ］は、対角行列である。従って、各固有振動モードは、相互に独立したものとして扱うことができる。よって、計算時間を短くすることができる。
物理座標系における台車枠１６の変位ベクトル｛ｕ｝は、以下の（９）式のように、モード座標系における変位ベクトル｛ξ｝（∈Ｒ^ｎ）に変換される。物理座標系における台車枠１６の外力ベクトル｛ｆ｝は、以下の（１０）式のように、モード座標系における外力ベクトル｛ｆ_ξ｝（∈Ｒ^ｎ）に変換される。

（１）式の両辺に、左から［φ］^Ｔを乗じる。（１）式の両辺に、左から［φ］^Ｔを乗じた式に、（９）式と（１０）式とを代入する。そうすると、以下の（１１）式が得られる。

ここで、ξ・は、（１１）式においてξの上に・が付されているものに対応する。ξ・・は、（１１）式においてξの上に・・が付されているものに対応する（このことは、以降の式でも同じである）。
（３）式～（５）式を（１１）式に代入すると、［φ］^Ｔ［φ］は単位行列であるので、以下の（１２）式が得られる。

以上のように、（１）式に示す物理座標系の運動方程式を、（１２）式のようにモード座標系の運動方程式で表現することができる。
ここで、前述の方法で導出した（１）式における外力ベクトル｛ｆ｝を、モード座標系における外力ベクトル｛ｆ_ξ｝として与えられる必要がある。
そこで、（１０）式を変形して得られる（１３）式を用いる。

本実施形態では、（１２）式の運動方程式を解くことにより、モード座標系における台車枠１６の変位ベクトル｛ξ｝（＝［ξ^（１） ξ^（２）・・・ξ^（ｎ）］^Ｔ）を導出する。
そのために、まず、モード行列［φ］を、固有値解析により導出する。モード行列［φ］と、台車枠１６の質量行列［Ｍ］、粘性行列［Ｃ］、および剛性行列［Ｋ］と、を用いて、（３）式～（５）式により、モード座標系における台車枠１６の質量行列［Ｍ_ξ］、粘性行列［Ｃ_ξ］、および剛性行列［Ｋ_ξ］を導出する。

また、加速度センサ２１ａ、２１ｂで測定された加速度のデータから、時間に関する一階積分、二階積分を行うことにより、（２）式の左辺の変位ベクトル｛ｕ_０｝、速度ベクトル｛ｕ_０・｝をそれぞれ導出する。モノリンク１８ａ、１８ｂの粘性減衰係数および軸バネ１９ａ、１９ｂの粘性減衰係数から、モノリンク１８ａ、１８ｂの粘性行列［Ｃ_ｂｃ］および軸バネ１９ａ、１９ｂの粘性行列［Ｃ_ｂｃ］をそれぞれ導出する。モノリンク１８ａ、１８ｂの剛性および軸バネ１９ａ、１９ｂの剛性から、モノリンク１８ａ、１８ｂの剛性行列［Ｋ_ｂｃ］および軸バネ１９ａ、１９ｂの剛性行列［Ｋ_ｂｃ］をそれぞれ導出する。このようにして得られた情報を（２）式に与えることにより、モノリンク１８ａ、１８ｂの着力点３１ａ、３１ｂに作用する外力および軸バネ１９ａ、１９ｂの着力点３２ａ、３２ｂに作用する外力をそれぞれ導出する。そして、台車枠１６の着力点３１ａ、３１ｂに作用する外力および着力点３２ａ、３２ｂに作用する外力を、それぞれ、モノリンク１８ａ、１８ｂの着力点３１ａ、３１ｂに作用する外力の反作用力および軸バネ１９ａ、１９ａの着力点３２ａ、３２ｂに作用する外力の反作用力として導出する。（１３）式の外力ベクトル｛ｆ｝の着力点３１ａ、３１ｂ、３２ａ、３２ｂに作用する外力の成分には、このようにして導出した値を与え、その他の成分には、０（ゼロ）を与える。そして、モード行列［φ］を（１３）式の右辺に与える。そして、（１３）式により、モード座標系における台車枠１６の外力ベクトル｛ｆ_ξ｝を導出する。

そして、以上のようにして導出された、モード座標系における台車枠１６の質量行列［Ｍ_ξ］、粘性行列［Ｃ_ξ］、剛性行列［Ｋ_ξ］、外力ベクトル｛ｆ_ξ｝を（１２）式に与える。モード座標系における台車枠１６の質量行列［Ｍ_ξ］、粘性行列［Ｃ_ξ］、剛性行列［Ｋ_ξ］、外力ベクトル｛ｆ_ξ｝が与えられた（１２）式を解くことにより、モード座標系における台車枠１６の変位ベクトル｛ξ｝を導出することができる。状態変数モデルを用いて（１２）式を解く場合、状態変数モデルは、以下のように表される。
まず、状態変数Ξ（∈Ｒ^２ｎ）を以下の（１４）式のように定義する。また、（１２）式を以下の（１５）式の状態方程式で記述する。

ここで、状態遷移行列Ａ（∈Ｒ^{２ｎ×２ｎ}）は、以下の（１６）式で表される。

また、ベクトルＦ（∈Ｒ^ｎ）は、モード座標系における台車枠１６の外力ベクトル｛ｆ_ξ｝を格納するベクトルであり、以下の（１７）式で表される。

また、行列Ｇ（∈Ｒ^２ｎ×ｎ）は、以下の（１８）式で表される。

（１６）式～（１８）式は、状態変数Ξが（１４）式で表現されるものとして、（１２）式を（１５）式の形で表記することにより得られる。
尚、（１２）式は、常微分方程式であり、（１２）式を解く手法は、状態変数モデルを用いる手法に限定されない。（１２）式を解く手法として、常微分方程式を解くための公知の手法を採用することができる。

以上のようにしてモード座標系における台車枠１６の変位ベクトル｛ξ｝（＝［ξ^（１） ξ^（２）・・・ξ^（ｎ）］^Ｔ）が導出される。そして、モード座標系における台車枠１６の変位ベクトル｛ξ｝を、物理座標系における台車枠１６の変位ベクトル｛ｕ｝に変換する。例えば、以下の（１９）式に示すようにして、メッシュの節点ｑにおけるｘ_ｉ軸方向（ｉ∈｛１，２，３｝）の変位ｕ_ｑ，ｉを導出することができる。

ここで、φ_ｑ，ｉ ^（ｎ）は、ｎ次の固有振動モード、メッシュの節点ｑ、およびｘ_ｉ軸方向に対応するモード行列［φ］の成分である。以上のようにして、モード座標系における台車枠１６の変位ベクトル｛ξ｝は、物理座標系における台車枠１６の変位ベクトル｛ｕ｝に変換される。これにより、台車枠１６の変位分布（各メッシュにおける変位）が得られる。

＜＜歪み－変位関係式＞＞
台車枠１６の変位分布が得られれば、メッシュの各要素ｅの内部における歪みテンソルの各成分ε_ｅ,１１、ε_ｅ,１２、ε_ｅ,１３、ε_ｅ,２１、ε_ｅ,２２、ε_ｅ,２３、ε_ｅ,３１、ε_ｅ,３２、ε_ｅ,３３は、以下の（２０ａ）式～（２０ｆ）式で表される。

ここで、ｘ_１、ｘ_２、ｘ_３は、変位の３つの自由度を記述する座標成分（ここでは、ｘ_１軸方向成分、ｘ_２軸方向成分、ｘ_３軸方向成分）を表す。また、ｕ_ｅ，１、ｕ_ｅ，２、ｕ_ｅ，３は、それぞれ各要素ｅの内部における変位のｘ_１軸方向成分、ｘ_２軸方向成分、ｘ_３軸方向成分を表す。

＜＜応力－歪み関係式＞＞
応力と歪みとの関係は、以下の（２１）式で表される。尚、（２１）式は、弾性体の構成則と呼ばれる。

ここで、σ_ｅ,１１、σ_ｅ,１２、σ_ｅ,１３の１１、１２、１３の下２桁目の１は、ｘ_１軸に垂直な微小面に作用する応力であることを表し、σ_ｅ,１１、σ_ｅ,２１、σ_ｅ,３１の１１、２１、３１の下２桁目の１は、ｘ_１軸方向成分の値であることを表す。σ_ｅ,２１、σ_ｅ,２２、σ_ｅ,２３の２１、２２、２３の下２桁目の２は、ｘ_２軸に垂直な微小面に作用する応力であることを表し、σ_ｅ,１２、σ_ｅ,２２、σ_ｅ,３２の１２、２２、３２の下１桁目の２は、ｘ_２軸方向成分の値であることを表す。σ_ｅ,３１、σ_ｅ,３２、σ_ｅ,３３の３１、３２、３３の下２桁目の３は、ｘ_３軸に垂直な微小面に作用する応力であることを表し、σ_ｅ,１３、σ_ｅ,２３、σ_ｅ,３３の１３、２３、３３の下１桁目の３は、ｘ_３軸方向成分の値であることを表す。
λ、μは、それぞれ、ラメ定数、剛性率である。ラメ定数λ、剛性率μは、それぞれ、以下の（２２）式、（２３）式で表される。

ここで、Ｅは、ヤング率である。νは、ポアソン比である。
以上のようにして、台車枠１６の応力分布（各メッシュにおける応力）が得られる。

＜推定装置４００の構成＞
図１に示すように、本実施形態では、推定装置４００は、鉄道車両の車体１１内に配置される。ただし、推定装置４００を車体１１内に配置せず、鉄道車両の外部に配置してもよい。このようにする場合、鉄道車両で測定されるデータは、例えば、無線通信により、鉄道車両から推定装置４００に送信されるようにする。

推定装置４００は、＜台車枠１６の応力分布の導出＞の項で説明したようにして、鉄道車両が走行中の台車枠１６の応力分布を導出する。推定装置４００は、当該台車枠１６の応力分布に基づく指標値が、所定の条件を満たすか否かを判定する。推定装置４００は、台車枠１６の応力分布に基づく指標値が、所定の条件を満たさない場合、当該所定の条件を満たさない台車枠１６の位置（箇所）と、当該所定の条件を満たさないと判定したタイミングにおける鉄道車両の走行位置とを導出する。推定装置４００は、当該所定の条件を満たさない台車枠１６の位置と、当該所定の条件を満たさないと判定したタイミングにおける鉄道車両の走行位置とを相互に関連付けて記憶する。このようにして記憶された情報は、例えば、台車枠１６の点検箇所および軌道２０の点検箇所を特定するための情報として使用することができる。軌道２０の点検箇所とは、具体的には、軌道不整などが生じている虞のある箇所である。台車枠１６の応力分布に基づく指標値は、台車枠１６および軌道２０の状態に関する指標値の一例である。

図４は、推定装置４００の機能的な構成の一例を示す図である。以下に、図４を参照しながら、推定装置４００が有する機能の一例を説明する。推定装置４００のハードウェアは、例えば、ＣＰＵ、ＲＯＭ、ＲＡＭ、ＨＤＤ、および各種のインターフェースを備える情報処理装置、または、専用のハードウェアを用いることにより実現される。尚、以下に説明する処理において、事前に設定する必要があるデータは、測定開始前に推定装置４００に記憶されているものとする。

＜＜データ取得部４０１＞＞
データ取得部４０１は、一定周期のサンプリング時刻が到来する度に、加速度センサ２１ａ、２１ｂで測定された加速度を示すデータと、鉄道車両の走行位置を示すデータと、を含むデータを取得する。鉄道車両の走行位置を示すデータを取得する方法は、特に限定されない。鉄道車両の走行位置を示すデータは、例えば、鉄道車両に配置されたＧＰＳ（Global Positioning System）を用いることにより得ることができる。また、鉄道車両の走行速度と、測定の開始時からの経過時間とに基づいて、各時刻における鉄道車両の、測定の開始時点からの移動距離を導出し、当該移動距離と、当該鉄道車両が走行する軌道２０の配置とに基づいて、鉄道車両の各時刻における走行位置を導出してもよい。このようにする場合、測定の開始時点における鉄道車両の位置を示すデータは、測定開始前に推定装置４００に記憶される。

固有値解析部４０２、質量導出部４０３、粘性導出部４０４、剛性導出部４０５、外力導出部４０６、状態導出部４０７、および応力－位置関係導出部４０８は、サンプリング時刻が到来する度に、以下に説明する処理を実行する。これらは、同一のサンプリング時刻に得られた情報を用いて、以下に説明する処理を実行する。例えば、各部が、或るサンプリング時刻の処理において、データ取得部４０１で取得されたデータを用いる場合、当該データは、当該サンプリング時刻においてデータ取得部４０１で取得されたデータになる。

＜＜固有値解析部４０２＞＞
固有値解析部４０２は、（１）式の運動方程式に対する固有値解析を行うことにより、固有振動数ωおよび固有ベクトル｛φ｝を導出し、固有ベクトル｛φ｝で構成されるモード行列［φ］を導出する。尚、有限要素法による運動方程式の定式化に基づいて、台車枠１６の質量行列［Ｍ］、粘性行列［Ｃ］、および剛性行列［Ｋ］が導出される。

＜＜質量導出部４０３＞＞
質量導出部４０３は、モード行列［φ］と、台車枠１６の質量行列［Ｍ］と、を用いて、（３）式により、モード座標系における台車枠１６の質量行列［Ｍ_ξ］を導出する。以下の説明では、モード座標系における台車枠１６の質量行列を、必要に応じて、モード質量行列と称する。

＜＜粘性導出部４０４＞＞
粘性導出部４０４は、モード行列［φ］と、台車枠１６の粘性行列［Ｃ］と、を用いて、（４）式により、モード座標系における台車枠１６の粘性行列［Ｃ_ξ］を導出する。以下の説明では、モード座標系における台車枠１６の粘性行列を、必要に応じて、モード粘性行列と称する。

＜＜剛性導出部４０５＞＞
剛性導出部４０５は、モード行列［φ］と、台車枠１６の剛性行列［Ｋ］と、を用いて、（５）式により、モード座標系における台車枠１６の剛性行列［Ｋ_ξ］を導出する。以下の説明では、モード座標系における台車枠１６の剛性行列を、必要に応じて、モード剛性行列と称する。

＜＜外力導出部４０６＞＞
外力導出部４０６は、加速度センサ２１ａ、２１ｂで測定された加速度のデータに基づいて、モノリンク１８ａ、１８ｂの変位ベクトル｛ｕ_０｝および速度ベクトル｛ｕ_０・｝および軸バネ１９ａ、１９ｂの変位ベクトル｛ｕ_０｝および速度ベクトル｛ｕ_０・｝をそれぞれ導出する。外力導出部４０６は、モノリンク１８ａ、１８ｂの粘性減衰係数および軸バネ１９ａ、１９ｂの粘性減衰係数から、モノリンク１８ａ、１８ｂの粘性行列［Ｃ_ｂｃ］および軸バネ１９ａ、１９ｂの粘性行列［Ｃ_ｂｃ］をそれぞれ導出する。また、外力導出部４０６は、モノリンク１８ａ、１８ｂの剛性および軸バネ１９ａ、１９ｂの剛性から、モノリンク１８ａ、１８ｂの剛性行列［Ｋ_ｂｃ］および軸バネ１９ａ、１９ｂの剛性行列［Ｋ_ｂｃ］をそれぞれ導出する。

外力導出部４０６は、モノリンク１８ａ、１８ｂの変位ベクトル｛ｕ_０｝、速度ベクトル｛ｕ_０・｝、粘性行列［Ｃ_ｂｃ］、剛性行列［Ｋ_ｂｃ］および軸バネ１９ａ、１９ｂの変位ベクトル｛ｕ_０｝、速度ベクトル｛ｕ_０・｝、粘性行列［Ｃ_ｂｃ］、剛性行列［Ｋ_ｂｃ］を用いて、（２）式により、モノリンク１８ａ、１８ｂの着力点３１ａ、３１ｂに作用する外力および軸バネ１９ａ、１９ｂの着力点３２ａ、３２ｂに作用する外力を導出する。外力導出部４０６は、これらの反作用力として、台車枠１６の着力点３１ａ、３１ｂに作用する外力および着力点３２ａ、３２ｂに作用する外力を導出する。
外力導出部４０６は、台車枠１６の外力ベクトル｛ｆ｝の成分の内、着力点３１ａ、３１ｂ、３２ａ、３２ｂに作用する外力の成分には、このようにして導出した値を与え、その他の成分には、０（ゼロ）を与えることにより、台車枠１６の外力ベクトル｛ｆ｝を、台車枠１６の外力分布として導出する。外力導出部４０６は、このようにして導出された台車枠１６の外力ベクトル｛ｆ｝と、モード行列［φ］と、を用いて、（１３）式により、モード座標系における台車枠１６の外力ベクトル｛ｆ_ξ｝を導出する。以下の説明では、モード座標系における台車枠１６の外力ベクトルを、必要に応じて、モード外力ベクトルと称する。本実施形態では、例えば、台車枠１６の外力ベクトル｛ｆ｝が、台車枠に作用する外力の分布の一例である。また、モード外力ベクトル｛ｆ_ξ｝が、モード座標系における台車枠に作用する外力の分布の一例である。

＜＜状態導出部４０７＞＞
状態導出部４０７は、外力導出部４０６により導出された台車枠１６の外力ベクトル｛ｆ｝に基づいて、台車枠１６の状態を表す情報を導出する。本実施形態では、状態導出部４０７は、台車枠１６の状態を表す情報として、台車枠１６の変位の分布と台車枠１６の応力の分布とを含む情報を導出する。本実施形態では、状態導出部４０７は、変位分布導出部４０７ａと、応力分布導出部４０７ｂとを有する。

＜＜変位分布導出部４０７ａ＞＞
変位分布導出部４０７ａは、モード質量行列［Ｍ_ξ］、モード粘性行列［Ｃ_ξ］、モード剛性行列［Ｋ_ξ］、およびモード外力ベクトル｛ｆ_ξ｝を用いて、（１２）式により、モード座標系における台車枠１６の変位ベクトル｛ξ｝を、モード座標系における台車枠１６の変位分布として導出する。前述したように本実施形態では、変位分布導出部４０７ａは、（１４）式～（１８）式に示す状態変数モデルを用いて、モード座標系における台車枠１６の変位ベクトル｛ξ｝を導出する。以下の説明では、モード座標系における台車枠１６の変位ベクトルを、必要に応じて、モード変位ベクトルと称する。本実施形態では、例えば、モード変位ベクトル｛ξ｝が、モード座標系における台車枠の変位の分布の一例である。
変位分布導出部４０７ａは、モード変位ベクトル｛ξ｝と、モード行列［φ］と、を用いて、（１９）式により、物理座標系における台車枠１６の変位ベクトル｛ｕ｝を、物理座標系における台車枠１６の変位分布として導出する。以下の説明では、物理座標系における台車枠１６の変位ベクトルを、必要に応じて、実変位ベクトルと称する。本実施形態では、例えば、実変位ベクトル｛ｕ｝が、物理座標系における台車枠の変位の分布の一例である。

＜＜応力分布導出部４０７ｂ＞＞
応力分布導出部４０７ｂは、実変位ベクトル｛ｕ｝を用いて、（２０ａ）式～（２０ｆ）式により、メッシュの各要素ｅの内部における歪みテンソルを、台車枠１６の歪み分布として導出する。応力分布導出部４０７ｂは、メッシュの各要素ｅの内部における歪みテンソルと、台車枠１６のヤング率Ｅおよびポアソン比νとを用いて、（２１）式～（２３）式により、メッシュの各要素ｅの内部における応力テンソルを、台車枠１６の応力分布として導出する。

＜＜応力－位置関係導出部４０８＞＞
応力－位置関係導出部４０８は、台車枠１６の応力分布と、鉄道車両の走行位置と、を相互に関連付けて記憶する。以下の説明では、このようにして相互に関連付けられた、台車枠１６の応力分布および鉄道車両の走行位置を、必要に応じて、応力－位置関係情報と称する。以上のようにして各サンプリング時刻における応力－位置関係情報が得られる。

＜＜判定部４０９＞＞
判定部４０９は、応力－位置関係情報（応力の時間変化）に基づく指標値が、所定の条件を満たすか否かを判定する。本実施形態では、判定部４０９は、応力－位置関係情報（応力の時間変化）に基づいて定められる平均応力および変動応力を、ＪＩＳ規格Ｅ４２０７の応力限界図にプロットする。そして、判定部４０９は、応力限界図にプロットした点が、応力限界図の応力許容域内にあるか否かを判定する。例えば、判定部４０９は、或る計測区間の台車枠１６の応力と鉄道車両の走行位置（時間）との関係を示すグラフにおいて、応力の極大値を示す領域周辺に変動区間を設定する。そして、判定部４０９は、当該変動区間内での応力の平均値を平均応力とし、当該変動区間内での応力の変動幅（振幅）を変動応力とする。

図５は、応力限界図の一例を概念的に示す図である。図５では、応力許容域を斜線で示す。
図５において、σ_Ｂは、台車枠１６の材料の引張強さである。σ_０は、台車枠１６の材料の降伏に対する許容応力である。σ_ｗ１は、台車枠１６の材料の疲れ許容応力である。σ_ｗ２は、台車枠１６の溶接止端部を仕上げない場合の疲れ許容応力である。σ_ｗ３は、台車枠１６の溶接止端部を仕上げる場合の疲れ許容応力である。

このように、本実施形態では、平均応力および変動応力から定まる点が、応力分布に基づく指標値の一例となる。また、平均応力および変動応力から定まる点が、応力限界図の応力許容域内にあることが所定の条件を満たすことの一例である。

判定部４０９は、サンプリング時刻が到来する度に当該判定を行っても、予め定められたタイミングで当該判定を行っても、オペレータによる推定装置４００に対する所定の操作があったときに当該判定を行ってもよい。予め定められたタイミングとしては、例えば、鉄道車両の所定の区間における走行が終了したタイミングが挙げられる。

尚、応力分布に基づく指標値および所定の条件は、このようなものに限定されない。例えば、最大主応力を応力分布に基づく指標値の一例として採用し、最大主応力の絶対値が閾値を上回らないことを所定の条件を満たすことの一例として採用してもよい。

＜＜点検箇所特定情報導出部４１０＞＞
点検箇所特定情報導出部４１０は、応力－位置関係情報に基づいて定められる平均応力および変動応力を応力限界図にプロットした点が、応力限界図の応力許容域内にない場合、当該プロットが得られた台車枠１６の位置と、当該プロットした点が応力限界図の応力許容域内にないときの鉄道車両の走行位置と、を導出する。点検箇所特定情報導出部４１０は、当該導出を、応力－位置関係情報から設定される複数の変動区間のそれぞれにおいて行う。これにより、鉄道車両がどの位置を走行しているときに、台車枠１６のどの位置で、応力－位置関係情報に基づいて定められる平均応力および変動応力を応力限界図にプロットした点が、応力限界図の応力許容域内にないのかを特定することができる。以下の説明では、応力限界図の応力許容域内にないプロットが得られた台車枠１６の位置と、当該応力－位置関係情報に示される鉄道車両の走行位置とを示す情報を、必要に応じて、点検箇所特定情報と称する。尚、点検箇所特定情報は、応力限界図の応力許容域内にないプロットが得られた台車枠１６の位置と、当該応力－位置関係情報に示される鉄道車両の走行位置との何れか一方を示す情報であってもよい。また、台車枠１６全体の応力に基づいて、応力限界図にプロットを行ってもよい。

＜＜出力部４１１＞＞
出力部４１１は、点検箇所特定情報を出力する。出力の形態としては、例えば、コンピュータディスプレイへの表示、外部装置への送信、および推定装置４００の内部または外部の記憶媒体への記憶の少なくとも１つを採用することができる。
オペレータは、点検箇所特定情報に示される台車枠１６の位置に基づいて、台車枠１６の点検箇所を特定することができる。また、オペレータは、点検箇所特定情報に示される鉄道車両の走行位置に基づいて、軌道２０の点検箇所を特定することができる。

尚、点検箇所特定情報導出部４１０がこのような特定を行い、特定した点検箇所の情報を点検箇所特定情報としてよい。例えば、点検箇所特定情報導出部４１０は、応力限界図の応力許容域内にないプロットに基づいて特定される台車枠１６の位置を中心とする所定の範囲を、台車枠１６の点検箇所として特定する。また、例えば、推定装置４００は、応力－位置関係情報に基づいて定められる平均応力および変動応力を応力限界図にプロットした点が、応力限界図の応力許容域内にないときの鉄道車両の走行位置を中心とする所定の範囲を、軌道２０の点検箇所として特定する。

＜フローチャート＞
図６のフローチャートを参照しながら、本実施形態の推定装置４００を用いた推定方法の一例を説明する。図６のフローチャートのステップＳ６０１～Ｓ６１３の繰り返し処理は、サンプリング時刻が到来する度に実行されるものとする。
まず、ステップＳ６０１において、データ取得部４０１は、加速度センサ２１ａ、２１ｂで測定された加速度を示すデータと、鉄道車両の走行位置を示すデータと、を含むデータを取得する。
次に、ステップＳ６０２において、固有値解析部４０２は、（１）式の運動方程式に対する固有値解析を行うことにより、固有振動数ωおよび固有ベクトル｛φ｝を導出し、固有ベクトル｛φ｝で構成されるモード行列［φ］を導出する。

次に、ステップＳ６０３において、質量導出部４０３は、モード行列［φ］と、台車枠１６の質量行列［Ｍ］と、を用いて、（３）式により、モード質量行列［Ｍ_ξ］を導出する。
次に、ステップＳ６０４において、粘性導出部４０４は、モード行列［φ］と、台車枠１６の粘性行列［Ｃ］と、を用いて、（４）式により、モード粘性行列［Ｃ_ξ］を導出する。

次に、ステップＳ６０５において、剛性導出部４０５は、モード行列［φ］と、台車枠１６の剛性行列［Ｋ］と、を用いて、（５）式により、モード剛性行列［Ｋ_ξ］を導出する。
次に、ステップＳ６０６において、外力導出部４０６は、加速度センサ２１ａ、２１ｂで測定された加速度のデータに基づいて、モノリンク１８ａ、１８ｂの変位ベクトル｛ｕ_０｝および速度ベクトル｛ｕ_０・｝と、軸バネ１９ａ、１９ｂの変位ベクトル｛ｕ_０｝および速度ベクトル｛ｕ_０・｝とをそれぞれ導出する。外力導出部４０６は、モノリンク１８ａ、１８ｂの変位ベクトル｛ｕ_０｝、速度ベクトル｛ｕ_０・｝、粘性行列［Ｃ_ｂｃ］、および剛性行列［Ｋ_ｂｃ］と、モード行列［φ］と、軸バネ１９ａ、１９ｂの変位ベクトル｛ｕ_０｝、速度ベクトル｛ｕ_０・｝、粘性行列［Ｃ_ｂｃ］、および剛性行列［Ｋ_ｂｃ］と、モード行列［φ］と、を用いて、（２）式および（１３）式により、モード外力ベクトル｛ｆ_ξ｝を導出する。
尚、ステップＳ６０３～Ｓ６０６の順番は、順不同である。

次に、ステップＳ６０７において、変位分布導出部４０７ａは、モード質量行列［Ｍ_ξ］、モード粘性行列［Ｃ_ξ］、モード剛性行列［Ｋ_ξ］、およびモード外力ベクトル｛ｆ_ξ｝を用いて、（１２）式により、モード変位ベクトル｛ξ｝を導出する。変位分布導出部４０７ａは、モード変位ベクトル｛ξ｝と、モード行列［φ］と、を用いて、（１９）式により、実変位ベクトル｛ｕ｝を、台車枠１６の変位分布として導出する。

次に、ステップＳ６０８において、応力分布導出部４０７ｂは、実変位ベクトル｛ｕ｝を用いて、（２０ａ）式～（２０ｆ）式により、メッシュの各要素ｅの内部における歪みテンソルを、台車枠１６の歪み分布として導出する。
次に、ステップＳ６０９において、応力分布導出部４０７ｂは、メッシュの各要素ｅの内部における歪みテンソルと、台車枠１６のヤング率Ｅおよびポアソン比νと、を用いて、（２１）式～（２３）式により、メッシュの各要素ｅの内部における応力テンソルを、台車枠１６の応力分布として導出する。

次に、ステップＳ６１０において、応力－位置関係導出部４０８は、台車枠１６の応力分布と、鉄道車両の走行位置と、を相互に関連付けて記憶する（応力－位置関係情報を記憶する）。
次に、ステップＳ６１１において、判定部４０９は、応力－位置関係情報（応力の時間変化）に基づいて定められる平均応力および変動応力を応力限界図にプロットする。そして、判定部４０９は、プロットした点が、応力限界図の応力許容域内にあるか否かを判定する。この判定の結果、応力－位置関係情報（応力の時間変化）に基づいて定められる平均応力および変動応力を応力限界図にプロットした点が、応力限界図の応力許容域内にある場合、処理は、ステップＳ６１２を省略して後述するステップＳ６１３に進む。

一方、応力－位置関係情報（応力の時間変化）に基づいて定められる平均応力および変動応力を応力限界図にプロットした点が、応力限界図の応力許容範囲内にない場合、処理は、ステップＳ６１２に進む。処理がステップＳ６１２に進むと、点検箇所特定情報導出部４１０は、当該プロットが得られた台車枠１６の位置と、当該プロットした点が応力限界図の応力許容域内にないときの鉄道車両の走行位置と、を、点検箇所特定情報として導出する。
次に、ステップＳ６１３において、推定装置４００は、測定を終了するか否かを判定する。この判定は、例えば、オペレータによる推定装置４００に対する所定の操作があったか否かに基づいて行われる。また、この判定は、予め定められたタイミングになったか否かに基づいて行われてもよい。予め定められたタイミングとしては、例えば、鉄道車両の所定の区間における走行が終了したタイミングが挙げられる。

ステップＳ６１３の判定の結果、測定を終了しない場合、処理は、ステップＳ６０１に戻る。そして、次のサンプリング時刻におけるステップＳ６０１～Ｓ６１３の処理が実行される。
一方、ステップＳ６１３の判定の結果、測定を終了する場合、処理は、ステップＳ６１４に進む。ステップＳ６１４に進むと、出力部４１１は、点検箇所特定情報を出力する。そして、図６のフローチャートによる処理が終了する。

尚、ステップＳ６１４の処理を、ステップＳ６１２の後に行ってもよい。この場合、点検箇所特定情報は、サンプリング時刻の単位で出力される。この場合、ステップＳ６１１の判定の結果、応力－位置関係情報（応力の時間変化）に基づいて定められる平均応力および変動応力を応力限界図にプロットした点が、応力限界図の応力許容域内にある場合、出力部４１１は、そのことを示す情報を点検箇所特定情報として出力してもよい。この場合の点検箇所特定情報は、例えば、応力－位置関係情報（応力の時間変化）に基づいて定められる平均応力および変動応力を応力限界図にプロットした点が、応力限界図の応力許容域内にあることを示す情報と、当該プロットした点が応力限界図の応力許容域内にあるときの鉄道車両の走行位置とを示す情報と、を含む。点検箇所特定情報に含まれる鉄道車両の走行位置は点検が不要であることが特定される。

＜まとめ＞
以上のように本実施形態では、推定装置４００は、軸箱１７ａ、１７ｂに取り付けられた加速度センサ２１ａ、２１ｂで測定された加速度のデータに基づいて、台車枠１６の着力点３１ａ、３１ｂ、３２ａ、３２ｂに作用する外力を導出する。推定装置４００は、台車枠１６に作用する外力の分布を、台車枠１６の振動を表す運動方程式に与えて当該運動方程式を解くことにより、台車枠１６の変位分布を導出する。推定装置４００は、台車枠１６の変位分布を用いて台車枠１６の応力分布を導出する。従って、台車枠の検査箇所を予め指定し、台車枠の検査箇所にセンサを配置しなくても、鉄道車両が走行しているときの台車枠における応力分布を正確に推定することができる。よって、鉄道車両が走行しているときの台車枠の状態や軌道２０の状態を正確に推定することができる。

また、本実施形態では、推定装置４００は、台車枠１６の振動を表す運動方程式に対する固有値解析を行うことにより、固有ベクトルを導出し、固有ベクトルを格納したモード行列［φ］を導出する。推定装置４００は、加速度センサ２１ａ、２１ｂで測定された加速度のデータに基づいて、モノリンク１８ａ、１８ｂの変位ベクトル｛ｕ_０｝および速度ベクトル｛ｕ_０・｝と、軸バネ１９ａ、１９ｂの変位ベクトル｛ｕ_０｝および速度ベクトル｛ｕ_０・｝とを導出する。推定装置４００は、モノリンク１８ａ、１８ｂの粘性行列［Ｃ_ｂｃ］、剛性行列［Ｍ_ｂｃ］、変位ベクトル｛ｕ_０｝、および速度ベクトル｛ｕ_０・｝と、軸バネ１９ａ、１９ｂの粘性行列［Ｃ_ｂｃ］、剛性行列［Ｍ_ｂｃ］、変位ベクトル｛ｕ_０｝、および速度ベクトル｛ｕ_０・｝と、モード行列［φ］と、に基づいて、モード座標系における台車枠１６の外力ベクトル｛ｆ_ξ｝を導出する。推定装置４００は、モード座標系における台車枠１６の振動を表す運動方程式に、モード座標系における台車枠１６の外力ベクトル｛ｆ_ξ｝を与えて当該運動方程式を解くことにより、モード座標系における台車枠１６の変位分布を導出する。従って、鉄道車両の走行中に、台車枠における応力分布を正確に推定するための計算負荷を軽減することができる。

また、本実施形態では、軸箱１７ａ、１７ｂに取り付けられた加速度センサ２１ａ、２１ｂで測定された加速度のデータを用いる。従って、結合要素を構成する部品（モノリンク１８ａ、１８ｂおよび軸バネ１９ａ、１９ｂ）毎に加速度を測定する必要がなくなる。また、軸箱１７ａ、１７ｂは、力を吸収したり減衰したりすることを目的とするものではないので、加速度のデータに外乱が含まれる等、加速度の測定値の信頼性が低下することを抑制することができる。

＜変形例＞
本実施形態では、軸箱１７ａ、１７ｂに加速度センサ２１ａ、２１ｂを取り付ける場合を例に挙げて説明した。しかしながら、台車枠１６の着力点に作用する外力を導出することが可能な物理量を測定するためのセンサであれば、軸箱１７ａ、１７ｂに取り付けられるセンサは、加速度センサに限定されない。例えば、歪みゲージであってもよい。ただし、歪みゲージは、長期間の使用により、正常に動作しなくなる虞がある。このため、本実施形態では、台車枠１６の着力点に作用する外力を導出することが可能な物理量を測定するためのセンサとして、加速度センサを用いる。

前述したように、加速度センサ２１ａ、２１ｂを軸箱１７ａ、１７ｂに取り付けるのが好ましい。しかしながら、加速度センサ２１ａ、２１ｂの取り付け箇所は、軸箱１７ａ、１７ｂに限定されない。加速度センサ２１ａ、２１ｂの取り付け箇所は、台車枠１６と直接または他の部材を介して接続される部品であればよい。例えば、加速度センサ２１ａ、２１ｂの取り付け箇所は、加速度センサを取り付けるための専用の部品であって、結合要素に接続される部品であってもよい。

また、本実施形態では、台車枠１６の応力分布を導出する場合を例に挙げて説明した。しかしながら、鉄道車両が走行しているときの台車枠の状態を反映する物理量を導出していれば、必ずしもこのようにする必要はない。例えば、台車枠１６の応力分布を導出せずに、台車枠１６の変位分布を導出して出力してもよい。

また、前述したように、本実施形態のように、モード座標系における台車枠１６の変位分布を導出し、物理座標系における台車枠１６の変位分布に変換するのが好ましい。しかしながら、物理座標系からモード座標系への変換を行わずに、物理座標系における台車枠１６の変位分布を導出してもよい。

また、本実施形態では、推定装置４００は、軸箱１７ａ、１７ｂに取り付けられた加速度センサ２１ａ、２１ｂで測定された加速度のデータに基づいて、モノリンク１８ａ、１８ｂおよび軸バネ１９ａ、１９ｂを介して当該軸箱１７ａ、１７ｂに接続されている台車枠１６の着力点３１ａ、３１ｂ、３２ａ、３２ｂに作用する外力を導出する。しかしながら、必ずしもこのようにする必要はない。例えば、列車が複数の鉄道車両を有する場合、特定の１つの鉄道車両の軸箱１７ａ、１７ｂに加速度センサ２１ａ、２１ｂを取り付ける。当該特定の鉄道車両については、本実施形態で説明したように、当該特定の鉄道車両の軸箱１７ａ、１７ｂに取り付けられた加速度センサ２１ａ、２１ｂで測定された加速度のデータに基づいて、当該軸箱１７ａ、１７ｂに接続されている台車枠１６の着力点３１ａ、３１ｂ、３２ａ、３２ｂに作用する外力を導出する。一方、当該特定の鉄道車両以外の鉄道車両であって、当該特定の鉄道車両と同一の編成（列車）の鉄道車両については、当該特定の鉄道車両の軸箱１７ａ、１７ｂに取り付けられた加速度センサ２１ａ、２１ｂで測定された加速度のデータに基づいて、当該特定の鉄道車両以外の鉄道車両の軸箱１７ａ、１７ｂに接続されている台車枠１６の着力点３１ａ、３１ｂ、３２ａ、３２ｂに作用する外力を導出する。このように全ての鉄道車両（全ての台車１２ａ、１２ｂ）の軸箱１７ａ、１７ｂに加速度センサを取り付ける必要はない。このようにする場合、特定の鉄道車両以外の鉄道車両においては、少なくとも軌道２０の点検箇所を特定することができる。

（第２の実施形態）
次に、第２の実施形態を説明する。本実施形態では、台車枠１６の所定の位置の変位を導出することが可能な物理量を測定するためのセンサを台車枠１６に取り付ける。そして、本実施形態では、台車枠１６の振動を表す運動方程式を解くことにより導出される変位分布のうち、前記所定の位置（当該センサが取り付けられた箇所に対応する位置）の変位と、当該センサで測定された変位と、に基づいて、台車枠１６の振動を表す運動方程式を解くことにより導出される変位分布を補正するための補正パラメータを導出する。そして、台車枠１６の振動を表す運動方程式を解くことにより導出される変位分布と、補正パラメータと、を用いて、台車枠１６の変位分布を導出する。このように本実施形態は、第１の実施形態に対し、補正パラメータを用いることによる構成および処理が追加されたものである。従って、本実施形態の説明において、第１の実施形態と同一の部分については、図１～図６に付した符号と同一の符号を付す等して詳細な説明を省略する。

図７は、台車枠およびその周辺の部品の構成の一例を示す図である。図７は、図２に対し、加速度センサ２２ａ、２２ｂが追加されたものである。加速度センサ２２ａ、２２ｂは、例えば、加速度センサ２１ａ、２１ｂと同じもので実現することができる。加速度センサ２２ａ、２２ｂは、台車枠１６の所定の位置に取り付けられる。図７では、台車枠１６に取り付けられる加速度センサ２２ａ、２２ｂの数が２である場合を例に挙げて示す。しかしながら、台車枠１６に取り付けられる加速度センサの数は、１以上であれば、幾つであってもよい。また、台車枠１６の所定の位置の変位を導出することが可能な物理量を測定するためのセンサであれば、台車枠１６に取り付けられるセンサは、加速度センサに限定されない。例えば、歪みゲージであってもよい。ただし、歪みゲージは、長期間の使用により、正常に動作しなくなる虞がある。このため、本実施形態では、台車枠１６の所定の位置の変位を導出することが可能な物理量を測定するためのセンサとして、加速度センサを用いる。本実施形態では、例えば、加速度センサ２２ａ、２２ｂが第２のセンサの一例である。また、加速度センサ２２ａ、２２ｂにより測定される加速度（台車枠１６の加速度）が、台車枠の所定の位置の変位を導出することが可能な物理量の一例である。

＜推定装置４００の構成＞
以下に、第１の実施形態の推定装置４００が有する機能と異なる部分を中心に、本実施形態の推定装置４００が有する機能の一例を説明する。本実施形態の推定装置４００と第１の実施形態の推定装置４００とでは、データ取得部４０１および変位分布導出部４０７ａの機能の一部が異なる。

＜＜データ取得部４０１＞＞
第１の実施形態では、データ取得部４０１は、サンプリング時刻が到来する度に、加速度センサ２１ａ、２１ｂで測定された加速度を示すデータと、鉄道車両の走行位置を示すデータと、を含むデータを取得する。これに対し、本実施形態のデータ取得部４０１は、加速度センサ２１ａ、２１ｂで測定された加速度を示すデータと、鉄道車両の走行位置を示すデータと、を含むデータに加えて、加速度センサ２２ａ、２２ｂで測定されたデータを取得する。

＜＜状態導出部４０７＞＞
状態導出部４０７は、外力導出部４０６により導出された台車枠１６の外力ベクトル｛ｆ｝に基づいて、台車枠１６の状態を表す情報を導出する。本実施形態の状態導出部４０７は、第１の実施形態の状態導出部４０７と同様に、変位分布導出部４０７ａと、応力分布導出部４０７ｂとを有する。ただし、本実施形態の変位分布導出部４０７ａおよび応力分布導出部４０７ｂが有する機能の一部は、第１の実施形態の変位分布導出部４０７ａおよび応力分布導出部４０７ｂが有する機能と異なる。

＜＜変位分布導出部４０７ａ＞＞
変位分布導出部４０７ａは、第１の実施形態の変位分布導出部４０７ａと同様に、モード座標系における台車枠１６の変位ベクトル｛ξ｝（モード変位ベクトル｛ξ｝）を導出する。第１の実施形態の変位分布導出部４０７ａでは、（１９）式の計算を行う。これに対し、本実施形態の変位分布導出部４０７ａは、モード座標系における台車枠１６の変位ベクトル｛ξ｝を導出した後、修正係数ｃ_ｋを導出する。本実施形態では、修正係数ｃ_ｋが、モード座標系における台車枠１６の変位ベクトル｛ξ｝を補正するための補正パラメータの一例である。そして、本実施形態の変位分布導出部４０７ａは、修正係数ｃ_ｋと、モード座標系における台車枠１６の変位ベクトル｛ξ｝と、モード行列［φ］と、を用いて、物理座標系における台車枠１６の変位ベクトル｛ｕ｝を、台車枠１６の変位分布として導出する。以下に、修正係数ｃ_ｋを導出する方法の一例と、修正係数ｃ_ｋを用いて、物理座標系における台車枠１６の変位ベクトル｛ｕ｝を導出する方法の一例とを説明する。

本実施形態の説明では、ｊにより識別される有限要素法のメッシュの節点をｑ_ｊとする。メッシュの節点ｑ_ｊ（ｊ＝１，・・・，ｍ）におけるｘ_ｉ軸方向の変位の測定値をｕ～_ｑｊ，ｉとする。尚、ｕ～_ｑｊ，ｉは、各式において、ｕの上に～が付され、ｑｊは、ｑ_ｊと表記されるものである。また、モード座標系における台車枠１６の変位ベクトル｛ξ｝は［ｃ_１ξ^（１） ｃ_２ξ^（２）・・・ｃ_ｎξ^（ｎ）］^Ｔに補正される。そうすると、メッシュの節点ｑにおけるｘ_ｉ軸方向の変位ｕ_ｑｊ，ｉを導出する式は、（１９）式より、以下の（２４）式で表される。

ここで、コスト関数Ｊを以下の（２５）式で定義する。コスト関数Ｊの値を最小にする修正係数ｃ_ｋは、以下の（２６）式を解くことにより導出される。

ここで、αは、正則化係数である。正則化係数αは、正則化係数αが大きい（小さい）ほど、メッシュの節点ｑ_ｊにおけるｘ_ｉ軸方向の変位の測定値ｕ～_ｑｊ，ｉの信頼度が低くなる（高くなる）ように、予め定められる定数である。
コスト関数Ｊの値を最小にする修正係数ｃ_ｋは、例えば、以下の（２７）式および（２８）式で表される。尚、（２６）式は、公知の手法で解くことができるので、ここでは、コスト関数Ｊの値を最小にする修正係数ｃ_ｋの導出過程の詳細な説明を省略する。

変位分布導出部４０７ａは、以上のようにして、コスト関数Ｊの値を最小にする修正係数ｃ_ｋを導出する。変位分布導出部４０７ａは、コスト関数Ｊの値を最小にする修正係数ｃ_ｋを用いて、（２４）式により、メッシュの節点ｑ_ｊにおけるｘ_ｉ軸方向（ｉ∈｛１，２，３｝）の変位ｕ_ｑ，ｉを導出する。

以上のようにして物理座標系における台車枠１６の変位ベクトル（実変位ベクトル）｛ｕ｝が導出される。
応力分布導出部４０７ｂは、（１９）式を用いて導出される変位ｕ_ｑ，ｉを成分とする実変位ベクトル｛ｕ｝に代えて、（２４）式を用いて導出される変位ｕ_ｑ，ｉを成分とする実変位ベクトル｛ｕ｝を用いて、第１の実施形態で説明した処理を実行する。

＜フローチャート＞
以下に、第１の実施形態の推定装置４００を用いた推定方法と異なる部分を中心に、本実施形態の推定装置４００を用いた推定方法一例を説明する。本実施形態の推定装置４００を用いた推定方法と、第１の実施形態の推定装置４００を用いた推定方法とでは、ステップＳ６０１、Ｓ６０７の処理の一部が異なる。

ステップＳ６０１において、データ取得部４０１は、加速度センサ２１ａ、２１ｂで測定された加速度を示すデータと、鉄道車両の走行位置を示すデータと、加速度センサ２２ａ、２２ｂで測定されたデータと、を含むデータを取得する。
ステップＳ６０２～Ｓ６０６の処理は、第１の実施形態で説明した図６のフローチャートの処理と同じである。

ステップＳ６０７において、変位分布導出部４０７ａは、モード質量行列［Ｍ_ξ］、モード粘性行列［Ｃ_ξ］、モード剛性行列［Ｋ_ξ］、およびモード外力ベクトル｛ｆ_ξ｝を用いて、（１２）式により、モード変位ベクトル｛ξ｝を導出する。そして、変位分布導出部４０７ａは、（２７）式および（２８）式により、修正係数ｃ_ｋを導出する。修正係数ｃ_ｋの導出には、加速度センサ２２ａ、２２ｂで測定されたデータ（メッシュの節点ｑ_ｊにおけるｘ_ｉ軸方向の変位の測定値ｕ～_ｑｊ，ｉ）が用いられる。また、モード変位ベクトル｛ξ｝、モード行列［φ］のうち、加速度センサ２２ａ、２２ｂの位置に対応する成分が用いられる。そして、変位分布導出部４０７ａは、モード変位ベクトル｛ξ｝と、モード行列［φ］と、修正係数ｃ_ｋとを用いて、（２４）式により、実変位ベクトル｛ｕ｝を、台車枠１６の変位分布として導出する。
ステップＳ６０８以降の処理は、第１の実施形態で説明した図６のフローチャートの処理と同じである。

＜まとめ＞
以上のように本実施形態では、推定装置４００は、加速度センサ２２ａ、２２ｂにより測定された台車枠１６の所定の位置の変位と、台車枠１６の運動を記述するモード座標系の運動方程式を解くことにより導出されるモード変位ベクトル｛ξ｝に対応する実変位ベクトル｛ｕ｝のうち、前記所定の位置（加速度センサ２２ａ、２２ｂが取り付けられた箇所に対応する位置）の変位との差に基づいて、台車枠１６の運動を記述するモード座標系の運動方程式における変位ベクトル｛ξ｝を補正するための修正係数ｃ_ｋを導出する。そして、推定装置４００は、モード変位ベクトル｛ξ｝と、モード行列［φ］と、修正係数ｃ_ｋと、を用いて、台車枠１６の変位分布を導出する。従って、台車枠の変位分布を、より一層正確に推定することができる。これにより、台車枠における応力分布を正確に推定することができる。
本実施形態においても、第１の実施形態で説明した種々の変形例を採用することができる。

（第３の実施形態）
次に、第３の実施形態を説明する。第１の実施形態および第２の実施形態では、有限要素法のメッシュの全ての節点ｑに対応するモード行列［φ］の成分φ_ｑ，ｉ ^（ｎ）を用いて、メッシュの節点ｑにおけるｘ_ｉ軸方向の変位ｕ_ｑ，ｉを導出する場合を例に挙げて説明した。しかしながら、このようにする場合、メッシュの全ての節点ｑの数と、自由度の数と、モード数ｎとの積の数だけ、モード行列［φ］の成分φ_ｑ，ｉ ^（ｎ）が必要になる。また、所望の位置におけるモード行列［φ］の成分φ_ｑ，ｉ ^（ｎ）を探索することが容易ではない。このため、計算時間が長くなる虞がある。そこで、本実施形態では、モード行列［φ］（固有ベクトル｛φ｝）の成分の近似値を導出し、モード行列［φ］の成分の近似値を用いて、メッシュの節点ｑにおけるｘ_ｉ軸方向の変位ｕ_ｑ，ｉを導出する。このように本実施形態と第１の実施形態および第２の実施形態とは、モード行列［φ］の成分の近似値を導出することによる構成および処理が主として異なる。従って、本実施形態の説明において、第１の実施形態および第２の実施形態と同一の部分については、図１～図７に付した符号と同一の符号を付す等して詳細な説明を省略する。

＜モード行列［φ］の近似＞
本実施形態では、（１）式の運動方程式に固有値解析を適用することにより導出されるモード行列［φ］（固有ベクトル｛φ｝）の成分の一部を用いて、モード行列［φ］の成分の近似式を導出する。モード行列［φ］の成分の近似式は、モード行列［φ］の成分の近似値を計算するための式である。本実施形態では、カーネルリッジ回帰を行うことにより、モード行列［φ］の成分の近似式を導出する場合を例に挙げて説明する。更に、本実施形態では、カーネル関数としてガウスカーネルを用いる場合を例に挙げて説明する。

そこで、モード行列［φ］（∈Ｒ^３ｌ×ｎ）を、ガウス関数φ：Ｒ^３→Ｒで近似する。ここで、ｐ次の固有振動モード、位置ｓ、およびｘ_ｉ軸方向に対応するモード行列［φ］の成分の近似値φ’_ｓ，ｉ ^（ｐ）（ｘ_ｓ，１，ｘ_ｓ，２，ｘ_ｓ，３）を、以下の（２９）式で表す。尚、ｘ_ｉ軸方向は、位置ｓを定める座標系の自由度（ｉ＝１～３）の何れか１つを指す。

ここで、ｐは、モード数を識別する変数である（１≦ｐ≦ｎ）。ｊは、有限要素法のメッシュの節点ｑを識別する変数である（１≦ｊ≦ｍ）。本実施形態の説明では、ｊにより識別される有限要素法のメッシュの節点をｑ_ｊとする。ｘ_ｓ，ｉは、位置ｓの物理座標系におけるｘ_ｉ軸方向の座標である。ｘ_ｊ，ｉは、メッシュの節点ｑ_ｊの物理座標系におけるｘ_ｉ軸方向の座標である。位置ｓは、メッシュの節点ｑ_ｊと一致していても一致していなくてもよい。ｍは、（１）式の運動方程式に固有値解析を適用することにより導出されるモード行列［φ］の成分φ_ｑ，ｉ ^（ｎ）のうち、予め定められたメッシュの節点ｑ_ｊの数である。ｍの数は、有限要素法による数値解析で用いられるメッシュの節点ｑの数よりも少ない数である。

ψ_ｉ ^（ｐ）（ｘ_ｓ，１，ｘ_ｓ，２，ｘ_ｓ，３）は、ガウスカーネルであり、以下の（３０）式で表される。

ここで、σ_{ｉ，（ｐ）} ^２（∈Ｒ）は、ガウスカーネルの幅を調整するパラメータである。σ_{ｉ，（ｐ）} ^２は、３つの自由度（ｘ_ｉ軸）毎およびモード数ｐ毎に設定される。
（２９）式においてλ_ｊ，ｉ ^（ｐ）は、ガウスカーネルψ_ｉ ^（ｐ）（ｘ_ｓ，１，ｘ_ｓ，２，ｘ_ｓ，３）に対する重み係数である。本実施形態では、（３１）式の行列方程式を用いてカーネルリッジ回帰を行うことにより重み係数λ_ｊ，ｉ ^（ｐ）が導出される。

（３１）式の左辺の行列の各成分は、ｐ次の固有振動モード、メッシュの節点ｑ_ｊ（ｑ_ｊ＝１’，・・・，ｍ’）、およびｘ_ｉ軸方向に対応するモード行列［φ］の成分φ_ｑ，ｉ ^（ｐ）（ｘ_ｊ，１，ｘ_ｊ，２，ｘ_ｊ，３）である。

（３１）式において、ｘ_ｊ’，ｉは、メッシュの節点ｑ_ｊ’の物理座標系におけるｘ_ｉ軸方向の座標である。また、ｘ_ｊ，ｉは、メッシュの節点ｑ_ｊの物理座標系におけるｘ_ｉ軸方向の座標である。座標ｘ_ｊ’，ｉは、通常、座標ｘ_ｊ，ｉと同じである。例えば、座標ｘ_１’，１は、通常、座標ｘ_１，１と同じである。即ち、通常、（３１）式におけるｘ_ｊ’，ｉのアポストロフィ（’）はなくてよい。ただし、座標ｘ_ｊ’，ｉは座標ｘ_ｊ，ｉと異なっていてもよい。そこで、（３１）式において、メッシュの節点ｑ_ｊを識別する変数ｊとして、アポストロフィ（’）を付している変数とアポストロフィ（’）を付していない変数とを示す。従って、１’～ｍ’の個数と１～ｍの個数は同じである。即ち、（３１）式の左辺の行列は、ｍ行１列の行列（＝ｍ’行１列の行列）である。

（３１）式の左辺の行列をΦ（∈Ｒ^ｍ’×１）と表記する。（３１）式の右辺の１つ目の行列をΨ（∈Ｒ^ｍ’×ｍ）と表記する。（３１）式の右辺の２つ目の行列をΛ（∈Ｒ^ｍ×１）と表記する。そうすると、（３１）式は、以下の（３２）式のように表される。

行列Λに関する（３２）式の最小二乗解を得るため、以下の（３３）式に示すコスト関数Ｊ：Ｒ^ｍ→Ｒを用いる。

ここで、α（∈Ｒ）は正則化係数である。（３３）式の右辺第２項は、正則化項（ペナルティ項）である。行列Λは、コスト関数ＪをΛで偏微分して極小化する問題を解くことにより導出される。即ち、以下の（３４）式により行列Λが導出される。

ここで、Ｉ（∈Ｚ^ｍ×ｍ）は単位行列である。
一次の固有振動モード、・・・、ｎ次の固有振動モードのそれぞれに対して（３１）式を作成する。即ち、（３１）式のｐの値に、１，２，・・・，ｎを与えて、モード数ｐの数だけ（３１）式を作成する。そして、一次の固有振動モード、・・・、ｎ次の固有振動モードのそれぞれに対する（３１）式に基づいて、（３４）式の計算をそれぞれ実行する。これにより行列Λ（重み係数λ_ｊ，ｉ ^（ｐ））が導出される。行列Λ（重み係数λ_ｊ，ｉ ^（ｐ））を（２９）式に代入することにより、φ’_ｓ，ｉ ^（１）（ｘ_ｓ，１，ｘ_ｓ，２，ｘ_ｓ，３）、・・・、φ’_ｓ，ｉ ^（ｎ）（ｘ_ｓ，１，ｘ_ｓ，２，ｘ_ｓ，３）を計算する近似式（（２９）式）が導出される。

＜推定装置８００の構成＞
図８は、推定装置８００の機能的な構成の一例を示す図である。以下に、図８を参照しながら、第１の実施形態および第２の実施形態の推定装置４００が有する機能と異なる部分を中心に、本実施形態の推定装置８００が有する機能の一例を説明する。推定装置８００のハードウェアは、例えば、推定装置４００のハードウェアと同じもので実現される。尚、以下に説明する処理において、事前に設定する必要があるデータは、推定装置８００に事前に記憶されているものとする。

第１の実施形態および第２の実施形態の推定装置４００では、固有値解析部４０２で導出されたモード行列［φ］が、質量導出部４０３、粘性導出部４０４、剛性導出部４０５、外力導出部４０６、および状態導出部４０７（変位分布導出部４０７ａ）で用いられる。

これに対し、本実施形態では、近似式導出部８０１は、質量導出部４０３、粘性導出部４０４、剛性導出部４０５、外力導出部４０６、および状態導出部４０７の処理が開始される前に、固有値解析部４０２で導出されたモード行列［φ］を用いて、モード行列［φ］の成分の近似式（（２９）式）を導出する。その後、近似値導出部８０２は、モード行列［φ］の成分の近似式を用いて、ｐ次の固有振動モード、位置ｓ、およびｘ_ｉ軸方向に対応するモード行列［φ］の成分の近似値φ’_ｓ，ｉ ^（ｐ）（ｘ_ｓ，１，ｘ_ｓ，２，ｘ_ｓ，３）を導出する。位置ｓは、物理座標系における台車枠１６の変位分布を導出するために必要な位置である。位置ｓは、任意の位置である。位置ｓは、有限要素法のメッシュの節点ｑの位置と異なっていてもよい。位置ｓの数は、有限要素法のメッシュの節点ｑの数より少なくてもよい。ただし、位置ｓは、有限要素法のメッシュの節点ｑの位置と同じであってもよい。また、位置ｓの数は、有限要素法のメッシュの節点ｑの位置の数と同じであってもよい。

質量導出部４０３、粘性導出部４０４、剛性導出部４０５、外力導出部４０６、および変位分布導出部４０７ａは、固有値解析部４０２で導出されたモード行列［φ］に代えて、モード行列［φ］の成分の近似値φ’_ｓ，ｉ ^（ｐ）（ｘ_ｓ，１，ｘ_ｓ，２，ｘ_ｓ，３）を成分として含むモード行列［φ］を用いて、第１の実施形態および第２の実施形態で説明した処理を行う。本実施形態で用いられるモード行列［φ］の成分のうち、位置ｓに対応しない成分は０（ゼロ）になる。以下に、近似式導出部８０１および近似値導出部８０２が有する機能の一例を詳細に説明する。

＜＜近似式導出部８０１＞＞
近似式導出部８０１は、有限要素法のメッシュの節点ｑのうち、ｍ’個の節点とｍ個の節点とを選択する。この選択は、例えば、乱数を用いてランダムに行うことができる。また、ｍ’個の節点とｍ個の節点の位置は、オペレータによって指定されてもよい。尚、前述したように、ｍ個の節点の位置とｍ’個の節点の位置とを同じにしてもよい。

近似式導出部８０１は、固有値解析部４０２により導出されたモード行列［φ］の成分φ_ｑ，ｉ ^（ｐ）のうち、ｍ’個の節点ｑ_ｊに対応する成分φ_１’，ｉ ^（ｐ）（ｘ_１’，１，ｘ_１’，２，ｘ_１’，３）、・・・、φ_ｍ’，ｉ ^（ｐ）（ｘ_ｍ’，１，ｘ_ｍ’，２，ｘ_ｍ’，３）を抽出する。そして、近似式導出部８０１は、モード行列［φ］の成分φ_１’，ｉ ^（ｐ）（ｘ_１’，１，ｘ_１’，２，ｘ_１’，３）、・・・、φ_ｍ’，ｉ ^（ｐ）（ｘ_ｍ’，１，ｘ_ｍ’，２，ｘ_ｍ’，３）を（３１）式の左辺の行列（（３２）式の左辺の行列Φ）に代入する。また、近似式導出部８０１は、（３０）式に、位置ｓの物理座標系におけるｘ_ｉ軸方向の座標ｘ_ｓ，ｉと、メッシュの節点ｑ_ｊの物理座標系におけるｘ_ｉ軸方向の座標ｘ_ｊ，ｉとを代入する。そして、近似式導出部８０１は、（３０）式により、ガウスカーネルψ_ｉ ^（ｐ）（ｘ_ｓ，１，ｘ_ｓ，２，ｘ_ｓ，３）を導出する。近似式導出部８０１は、ガウスカーネルψ_ｉ ^（ｐ）（ｘ_ｓ，１，ｘ_ｓ，２，ｘ_ｓ，３）を、（３１）式の右辺の１つ目の行列（（３２）式の行列Ψ）に代入する。

そして、近似式導出部８０１は、（３４）式の計算を行うことにより、行列Λ（重み係数λ_ｊ，ｉ ^（ｐ））を導出する。
近似式導出部８０１は、以上の行列Λ（重み係数λ_ｊ，ｉ ^（ｐ））の導出を、一次の固有振動モード、・・・、ｎ次の固有振動モードのそれぞれに対して行う。
近似式導出部８０１は、以上のようにして、モード行列［φ］の成分の近似値φ’_ｓ，ｉ ^（１）（ｘ_ｓ，１，ｘ_ｓ，２，ｘ_ｓ，３）、・・・、φ’_ｓ，ｉ ^（ｎ）（ｘ_ｓ，１，ｘ_ｓ，２，ｘ_ｓ，３）を計算する近似式（（２９）式）としてｎ個の近似式を導出する。

＜＜近似値導出部８０２＞＞
近似値導出部８０２は、（２９）式のｘ_ｓ，１，ｘ_ｓ，２，ｘ_ｓ，３に、物理座標系における台車枠１６の変位分布を導出するために必要な位置ｓの物理座標系における座標を与える。そして、近似値導出部８０２は、（２９）式により、ｐ次の固有振動モード、位置ｓ、およびｘ_ｉ軸方向に対応するモード行列［φ］の成分の近似値φ’_ｓ，ｉ ^（ｐ）（ｘ_ｓ，１，ｘ_ｓ，２，ｘ_ｓ，３）を導出する。近似値導出部８０２は、ｐ次の固有振動モード、位置ｓ、およびｘ_ｉ軸方向に対応するモード行列［φ］の成分の近似値φ’_ｓ，ｉ ^（ｐ）（ｘ_ｓ，１，ｘ_ｓ，２，ｘ_ｓ，３）の導出を、一次の固有振動モード、・・・、ｎ次の固有振動モードのそれぞれについて行う。そして、近似値導出部８０２は、モード行列［φ］の成分の近似値φ’_ｓ，ｉ ^（ｐ）（ｘ_ｓ，１，ｘ_ｓ，２，ｘ_ｓ，３）を、モード行列［φ］の成分としてモード行列［φ］に格納する。このようにして近似値φ’_ｓ，ｉ ^（ｐ）（ｘ_ｓ，１，ｘ_ｓ，２，ｘ_ｓ，３）が格納されたモード行列［φ］は、固有ベクトル｛φ｝が格納されたモード行列［φ］の近似値となる。

＜フローチャート＞
以下に、第１の実施形態および第２の実施形態の推定装置４００を用いた推定方法と異なる部分を中心に、本実施形態の推定装置８００を用いた推定方法一例を説明する。本実施形態の推定装置８００を用いた推定方法と、第１の実施形態および第２の実施形態の推定装置４００を用いた推定方法とでは、ステップＳ６０２の処理の一部が異なる。また、図６のフローチャートによる処理が開始する前に、モード行列［φ］の成分の近似式（重み係数λ_ｊ，ｉ ^（ｐ））が近似式導出部８０１により導出されている必要がある。

ステップＳ６０１において、データ取得部４０１は、加速度センサ２１ａ、２１ｂで測定された加速度を示すデータと、鉄道車両の走行位置を示すデータと、を含むデータを取得する。本実施形態を第２の実施形態に適用する場合、データ取得部４０１は、これらのデータに加えて、加速度センサ２２ａ、２２ｂで測定されたデータを取得する。

次に、ステップＳ６０２において、近似値導出部８０２は、（２９）式により、ｐ次の固有振動モード、位置ｓ、およびｘ_ｉ軸方向に対応するモード行列［φ］の成分の近似値φ’_ｓ，ｉ ^（ｐ）（ｘ_ｓ，１，ｘ_ｓ，２，ｘ_ｓ，３）を導出する。そして、近似値導出部８０２は、モード行列［φ］の成分の近似値φ’_ｓ，ｉ ^（ｐ）（ｘ_ｓ，１，ｘ_ｓ，２，ｘ_ｓ，３）を、モード行列［φ］の成分としてモード行列［φ］に格納する。

ステップＳ６０３以降の処理は、第１の実施形態および第２の実施形態で説明した図６のフローチャートの処理と同じである。ただし、第１の実施形態および第２の実施形態では、固有ベクトル｛φ｝を含むモード行列［φ］が用いられる。これに対し、本実施形態では、モード行列［φ］の成分の近似値φ’_ｓ，ｉ ^（ｐ）（ｘ_ｓ，１，ｘ_ｓ，２，ｘ_ｓ，３）を含むモード行列［φ］が用いられる。

＜まとめ＞
以上のように本実施形態では、推定装置８００は、固有値解析部４０２で導出されたモード行列［φ］を用いて、カーネルリッジ回帰を行うことにより、モード行列［φ］の成分の近似式（（２９）式）を導出する。モード行列［φ］の成分の近似式は、モード行列［φ］の成分の近似値として、固有振動モードの次数ｐ、位置ｓ、位置ｓを定める座標系の自由度ｉに応じた近似値を計算する式である。推定装置８００は、モード行列［φ］の成分の近似式を用いて、モード行列［φ］の成分の近似値φ’_ｓ，ｉ ^（ｐ）（ｘ_ｓ，１，ｘ_ｓ，２，ｘ_ｓ，３）を導出する。従って、モード行列［φ］の成分として、任意の位置ｓにおける成分を導出することができる。よって、計算時間をより短くすることができる。

＜変形例＞
本実施形態では、カーネルリッジ回帰を行う際に、カーネル関数としてガウスカーネルを用いる場合を例に挙げて説明した。しかしながら、カーネル関数はガウスカーネルに限定されない。また、モード行列［φ］の成分の近似値φ’_ｓ，ｉ ^（ｐ）（ｘ_ｓ，１，ｘ_ｓ，２，ｘ_ｓ，３）を、カーネルリッジ回帰以外の回帰分析の手法を用いて導出してもよい。例えば、カーネル回帰を用いてもよい。

また、固有値解析部４０２および近似式導出部８０１は、推定装置８００に含まれていなくてもよい。即ち、固有値解析部４０２および近似式導出部８０１は、推定装置８００と異なる装置に含まれていてもよい。この場合、推定装置８００は、推定装置８００と異なる装置から、モード行列［φ］の成分の近似式（（２９）式）の情報を取得する。
その他、本実施形態においても、第１の実施形態および第２の実施形態で説明した種々の変形例を採用することができる。

（第４の実施形態）
次に、第４の実施形態を説明する。第１の実施形態、第２の実施形態、および第３の実施形態では、実変位ベクトル｛ｕ｝を用いて、（２０ａ）式～（２０ｆ）式により、メッシュの各要素ｅの内部における歪みテンソル（台車枠１６の歪み分布）を導出する。また、メッシュの各要素ｅの内部における歪みテンソルを用いて、（２１）式～（２３）式により、メッシュの各要素ｅの内部における応力テンソル（台車枠１６の応力分布）を導出する。（２０ａ）式～（２０ｆ）により、メッシュの各要素ｅの内部における歪テンソルを導出する際に、ＦＥＭの形状関数が必要になる。そこで、本実施形態では、実変位ベクトル｛ｕ｝に基づいて台車枠１６の主歪みおよび主応力を導出することを、ＦＥＭの形状関数を用いずに実行する。このように本実施形態と、第１の実施形態、第２の実施形態、および第３の実施形態とは、実変位ベクトル｛ｕ｝が導出された後の構成および処理が主として異なる。従って、本実施形態の説明において、第１の実施形態、第２の実施形態、および第３の実施形態と同一の部分については、図１～図８に付した符号と同一の符号を付す等して詳細な説明を省略する。本実施形態では、第１の実施形態の変位分布導出部４０７ａにより実変位ベクトル｛ｕ｝が導出される場合を例に挙げて説明する。

後述するように本実施形態では、台車枠１６の推定対象領域における主歪みを導出（推定）する。本実施形態では、台車枠１６の基準点を含む領域を、推定対象領域とする。基準点は、例えば、変位分布導出部４０７ａにより導出された実変位ベクトル｛ｕ｝に含まれる変位ｕ_ｑ，ｉを示すメッシュの節点ｑのうち、任意のメッシュの節点ｑである。基準点は、例えば、変位分布導出部４０７ａにより導出された実変位ベクトル｛ｕ｝に含まれる変位ｕ_ｑ，ｉを示すメッシュの節点ｑのうち、台車枠１６において亀裂が発生しやすい点（位置）に対応する。例えば、鉄道車両を走行させた場合に台車枠１６に発生する応力を有限要素法等の数値解析により導出し、導出した応力に基づいて基準点を決定することができる。例えば、最大主応力の絶対値が最大となる点を基準点として決定する。また、鉄道車両を実際に走行させることにより基準点を決定してもよい。例えば、過去に発生した台車枠１６の亀裂の発生箇所や、台車枠１６に取り付けた振動計の測定値に基づいて、基準点を決定してもよい。

推定対象領域は、基準点と、実変位ベクトル｛ｕ｝に含まれる変位ｕ_ｑ，ｉを示すメッシュの節点ｑのうち、基準点を除く複数のメッシュの節点ｑと、を含む大きさを有する。推定対象領域の境界（外縁）は、例えば、当該基準点を除く複数のメッシュの節点ｑの位置に基づいて定められる。例えば、当該基準点を除く複数のメッシュの節点ｑの位置を通るように、当該基準点を除く複数のメッシュの節点ｑの位置の間を線分で結ぶことにより定まる線を境界線とする領域を推定対象領域とすることができる。

一般的に、基準点における最大主応力と推定対象領域における最大主応力とは一致しない。そこで、推定対象領域における最大主応力で基準点における最大主応力を近似する場合には、基準点を除く複数のメッシュの節点ｑの位置は、例えば、以下のようにして決定することができる。まず、推定対象領域の候補に含まれる基準点に歪みゲージを取り付けて鉄道車両を走行させ、歪みゲージの測定値に基づいて推定対象領域の候補における最大主応力を導出する。また、後述するようにして推定対象領域の候補における最大主応力を導出する。このようにして推定対象領域の候補における最大主応力を導出する場合には、歪ゲージの測定値を用いない。推定対象領域の候補を異ならせて、以上のような推定対象領域の候補における最大主応力の導出を、歪みゲージを用いる場合と歪ゲージを用いない場合とのそれぞれについて実行する。そして、歪みゲージの測定値に基づいて導出される最大主応力と、歪みゲージの測定値を用いずに後述するようにして導出される最大主応力とが近くなる（好ましくは一致する）推定対象領域の候補を探索する。このようにして探索された推定対象領域の候補における基準点を除く複数のメッシュの節点ｑの位置を、基準点を除く複数のメッシュの節点ｑの位置として決定する。また、以上のようにして基準点を除く複数のメッシュの節点ｑの位置を決定するに際し、例えば、鉄道車両を走行させた場合の台車枠１６の歪みを有限要素法等の数値解析により導出し、導出した歪みを、歪みゲージの測定値に代えて用いてもよい。

本実施形態では、３次元の変位（ｘ_１軸方向成分の変位、ｘ_２軸方向成分の変位、およびｘ_３軸方向成分の変位）を成分として含むベクトルが実変位ベクトル｛ｕ｝として変位分布導出部４０７ａにより導出される場合を例に挙げて説明する。即ち、物理座標系は３軸の座標系である。尚、後述するように３次元の成分の変位が変位分布導出部４０７ａにより導出される場合には、実変位ベクトル｛ｕ｝に含まれる変位ｕ_ｑ，ｉを示すメッシュの節点ｑのうち、推定対象領域に含めるメッシュの節点ｑの数は４以上である。２次元の成分の変位が変位分布導出部４０７ａにより導出される場合には、実変位ベクトル｛ｕ｝に含まれる変位ｕ_ｑ，ｉを示すメッシュの節点ｑのうち、推定対象領域に含めるメッシュの節点ｑの数は３以上である。また、同一の台車枠１６に複数の基準点を設定する場合には、それぞれの基準点に対する推定対象領域が設定される。以下の説明では、実変位ベクトル｛ｕ｝に含まれる変位ｕ_ｑ，ｉを示すメッシュの節点ｑのうち、推定対象領域に含まれるメッシュの節点ｑを、必要に応じて、推定対象領域に含まれるメッシュの節点ｑと略称する。

＜変位に基づく歪みの導出＞
本発明者らは、推定対象領域における変位前の位置を示す座標をアフィン変換することにより変位後の位置を示す座標が表現されることに着目した。尚、当該変位前の位置、変位後の位置を示す座標は、物理座標系における座標（実空間座標）である。以下の説明では、物理座標系における座標を必要に応じて位置座標と称する。
変位前の位置座標を格納する縦ベクトルをＸ_１とし、変位後の位置座標を格納する縦ベクトルをＸ_２とする。縦ベクトルＸ_１、Ｘ_２は、以下の（３５）式で表される。

ここで、Ｍは、変位前の位置座標を格納する縦ベクトルＸ_１に対して乗算される行列である。変位前の位置座標を格納する縦ベクトルＸ_１と行列Ｍとの積により、変位前の位置座標を格納する縦ベクトルＸ_１は線型変換される。ｖは、変位前の位置座標を格納する縦ベクトルＸ_１と行列Ｍとの積に加算される縦ベクトルである。変位前の位置座標を格納する縦ベクトルＸ_１と行列Ｍとの積に、縦ベクトルｖを加算することにより、縦ベクトルＭＸ_１は平行移動される。縦ベクトルＭＸ_１は、変位前の位置座標を格納する縦ベクトルＸ_１と行列Ｍとの積である。

（１）式において、縦ベクトルＸ_１、Ｘ_２、ｖの下に余分な行の成分として１を加えると共に、行列Ｍの全ての列の下に０のみからなる余分な行を加えることにより、縦ベクトル［Ｘ_１ ^Ｔ，１］^Ｔ、［Ｘ_２ ^Ｔ，１］^Ｔ、［ｖ^Ｔ，１］^Ｔと、行列［Ｍ^Ｔ，０］^Ｔとが構成される。Ｔは転置行列であることを表す（このことは以降の説明でも同じである）。そして、平行移動を表す縦ベクトルｖを含む縦ベクトル（行列）［ｖ^Ｔ，１］^Ｔを、行列Ｍを含む行列［Ｍ^Ｔ，０］^Ｔの右に加える。すると、以下の（３６）式のように（３５）式を行列の積で表すことができる。

（３５）式および（３６）式はアフィン変換を表す式であり等価な式である。行列Ｍが正則行列であれば、（３６）式の右辺の左側の行列は群の構造を有する。ここでは、（３６）式の右辺の左側の行列をアフィン変換群と称することとする。アフィン変換群（群の構造を有する行列）は行列の積で以て群演算ができる。そして、アフィン変換群に属する行列は、物体の変形（回転，せん断，引張，圧縮など）や平行移動といった挙動を表現する行列の積という形に分解することが可能である。

変位前のｘ_１軸方向の位置座標、ｘ_２軸方向の位置座標、ｘ_３軸方向の位置座標をそれぞれ、ｘ_１1、ｘ_２1、ｘ_３1とする。変位後のｘ_１軸方向の位置座標、ｘ_２軸方向の位置座標、ｘ_３軸方向の位置座標をそれぞれ、ｘ_１２、ｘ_２２、ｘ_３２とする。すると、変位前の位置座標を表す縦ベクトルＸ_１と、変位後の位置座標を表す縦ベクトルＸ_２はそれぞれ、以下の（３７ａ）式、（３７ｂ）式で表される。尚、ｘ_１1、ｘ_２1、ｘ_３1、ｘ_１２、ｘ_２２、ｘ_３２の１１、２１、３１、１２、２２、３２の下１桁目の１、２は、それぞれ、変位前、変位後を表す。ｘ_１1、ｘ_２1、ｘ_３1、ｘ_１２、ｘ_２２、ｘ_３２の１１、２１、３１、１２、２２、３２の下２桁目の１、２、３は、それぞれ、ｘ_１軸方向成分の値、ｘ_２軸方向成分の値、ｘ_３軸方向成分の値であることを表す。

また、行列Ｍの成分をａ_１、ａ_２、ａ_３、ａ_５、ａ_６、ａ_７、ａ_９、ａ_１０、ａ_１１とし、縦ベクトルｖの成分をａ_４、ａ_８、ａ_１２とする。すると、（３６）式は、以下の（３８）式のようになる。

（３８）式を連立方程式に書き直すと以下の（３９）式のようになる。（３９）式の連立方程式の１２個の未知係数ａ_１～ａ_１２を決定するためには、推定対象領域に含めるメッシュの節点ｑが４箇所以上あればよい。推定対象領域に含めるメッシュの節点ｑが４箇所以上あれば、ｘ_１軸方向、ｘ_２軸方向、ｘ_３軸方向の変位の組が４つ以上得られるからである。
変位分布導出部４０７ａにより導出された実変位ベクトル｛ｕ｝に含まれる変位ｕ_ｑ，ｉのうち、推定対象領域に含まれるメッシュの節点ｑにおける変位を、変位前の位置座標に加算することにより、変位後のｘ_１軸方向、ｘ_２軸方向、ｘ_３軸方向の位置座標の組が得られる。

以上のようにして得られた変位前および変位後のｘ_１軸方向、ｘ_２軸方向、ｘ_３軸方向の位置座標の組として４つ以上の組を用いることによって、未知係数ａ_１～ａ_１２に関する１２以上の数の方程式が得られることから、例えば、最小二乗法等の手段を使い未知係数ａ_１～ａ_１２を決定することができる。

ここで、（３８）式より（３６）式の右辺の左側の行列を以下の（４０）式のように表す。即ち、（３６）式の右辺の左側の行列を、以下の（４１）式、（４２）式に示す行列Ｍ、縦ベクトルｖをブロックとする区分行列に分解する。

行列Ｍ∈Ｒ^３×３は回転を含む一般線形群ＧＬ_３であり、縦ベクトルｖ∈Ｒ^３は平行移動を表す正規部分群である。Ｒは（３６）式と同様に実数全体の集合を表す（このことは以降の説明でも同じである）。以下の（４３）式の行列（（３６）式の右辺の左側の行列）は、アフィン変換群のひとつの元である。アフィン変換群の性質から、積演算は、以下の（４４）式のように群の半直積演算で記述される。

すると、Ｉ∈Ｒ^３×３を単位行列として半直積演算のルールを用いれば、アフィン変換群の元は、以下の（４５）式のように分解される。

以下の（４６）式のように行列Ｍを特異値分解する。行列Ｍを特異値分解した結果を（４５）式に代入すると、以下の（４７）式が得られる。

ここで、Ｕ∈Ｒ^３×３は、直交行列である。Σ∈Ｒ^３×３は、特異値を対角成分として有する対角行列である。Ｖ^Ｔは、直交行列Ｖ∈Ｒ^３×３の転置行列を表す。
半直積の演算のルールを用いて（４７）式を、以下の（４８）式のように分解する。

（４８）式は、（４３）式で表されるアフィン変換が、以下の（１）～（４）の一連の変形・平行移動操作を、（１）、（２）、（３）、（４）の順に実行して得られるものであることを意味する。
（１） 転置行列Ｖ^Ｔの中から回転行列を選択し、選択した回転行列に基づく回転角で物体を回転させる。
（２） 特異値を対角成分として有する対角行列Σに基づいて物体を歪み変形させる。
（３） 直交行列Ｕの中から回転行列を選択し、選択した回転行列に基づく回転角で物体を再び回転させる。
（４） 縦ベクトルｖに基づいて物体を平行移動する。

本発明者らは、アフィン変換が（４８）式で表されることに着目し、行列Ｍを特異値分解することにより導出される対角行列Σの対角成分が主歪みに対応することを見出した。
前述したように、変位前のｘ_１軸方向の位置座標ｘ_１1、ｘ_２軸方向の位置座標ｘ_２1、ｘ_３軸方向の位置座標ｘ_３1と変位とに基づいて、変位後のｘ_１軸方向の位置座標ｘ_１２、ｘ_２軸方向の位置座標ｘ_２２、ｘ_３軸方向の位置座標ｘ_３２が導出される。そして、変位前のｘ_１軸方向の位置座標ｘ_１1、ｘ_２軸方向の位置座標ｘ_２1、ｘ_３軸方向の位置座標ｘ_３1と、変位後のｘ_１軸方向の位置座標ｘ_１２、ｘ_２軸方向の位置座標ｘ_２２、ｘ_３軸方向の位置座標ｘ_３２とに基づいて、未知係数ａ_１～ａ_１２が導出される。即ち、行列Ｍの成分ａ_１、ａ_２、ａ_３、ａ_５、ａ_６、ａ_７、ａ_９、ａ_１０、ａ_１１が導出される。そして、行列Ｍを特異値分解することにより、特異値を対角成分として有する対角行列Σが導出される。

対角行列Σの対角成分をε_１１、ε_２２、ε_３３とする。対角行列Σの対角成分ε_１１、ε_２２、ε_３３は、未知係数ａ_１～ａ_１２を導出する際に用いたｈ個の位置座標の集合｛（ｘ_１1，ｘ_２1，ｘ_３1）^Ｔ _１、・・・、（ｘ_１1，ｘ_２1，ｘ_３1）^Ｔ _ｈ｝で構成される基準領域（推定対象領域）が、弾性変形により移動し、新たな位置座標の集合｛（ｘ_１２，ｘ_２２，ｘ_３２）^Ｔ _１、・・・、（ｘ_１２，ｘ_２２，ｘ_３２）^Ｔ _ｈ｝に変化したときの、前記基準領域（推定対象領域）の主歪み成分に該当する。ここで、ｈ≧４である。このような対角行列Σの対角成分ε_１１、ε_２２、ε_３３のうち絶対値が最大となるものが最大主歪みに対応する。
本実施形態は、以上の知見に基づいてなされたものである。

本実施形態では、３次元の成分の変位が変位分布導出部４０７ａにより導出される場合を例に挙げて説明する。２次元の成分の変位が変位分布導出部４０７ａにより導出される場合には、３次元の成分の変位が変位分布導出部４０７ａにより導出される場合に対し、変位分布導出部４０７ａにより導出されない成分（軸方向）の値を０（ゼロ）とすることにより実現することができる。従って、ここでは、２次元の成分の変位が変位分布導出部４０７ａにより導出される場合の詳細な説明を省略する。例えば、推定対象領域がモノリンク１８ａ、１８ｂの場合には、ｘ_２軸方向に生じる変位が無視できるほど小さい。このため、ｘ_１軸方向とｘ_３軸方向との２次元の変位を変位分布導出部４０７ａにより導出すればよい。尚、２次元の成分の変位が変位分布導出部４０７ａにより導出される場合には、（３９）式の連立方程式において未知係数の数は６個になる。推定対象領域に含めるメッシュの節点ｑが３箇所以上あれば、ｘ_１軸方向、ｘ_２軸方向、ｘ_３軸方向のうち、変位分布導出部４０７ａにより導出される２つの軸方向の変位の組が３個以上得られる。従って、６個の未知係数に関する６以上の数の方程式が得られる。よって、６個の未知係数を決定することができる。

＜推定装置９００の構成＞

図９は、推定装置９００の機能的な構成の一例を示す図である。以下に、図９を参照しながら、第１の実施形態、第２の実施形態、および第３の実施形態の推定装置４００、８００が有する機能と異なる部分を中心に、本実施形態の推定装置９００が有する機能の一例を説明する。推定装置９００のハードウェアは、例えば、推定装置４００、８００のハードウェアと同じもので実現される。尚、以下に説明する処理において、事前に設定する必要があるデータは、測定開始前に推定装置９００に記憶されているものとする。

第１の実施形態、第２の実施形態、および第３の実施形態の推定装置４００、８００では、変位分布導出部４０７ａにより導出された実変位ベクトル｛ｕ｝を用いて、応力分布導出部４０７ｂにより、台車枠１６の応力分布（メッシュの各要素ｅの内部における応力テンソル）が導出される。

これに対し、本実施形態では、状態推定部４０７は、変位分布導出部４０７ａにより導出された実変位ベクトル｛ｕ｝を用いて、推定対象領域における主歪み、最大主歪み、主応力、および最大主応力を、台車枠１６の状態を表す情報の一例として導出する。本実施形態では、状態推定部４０７は、変位分布導出部４０７ａ、位置座標導出部４０７ｃ、行列導出部４０７ｄ、特異値分解部４０７ｅ、歪み導出部４０７ｆ、および応力導出部４０７ｇを有する。

＜＜位置座標導出部４０７ｃ＞＞
位置座標導出部４０７ｃは、推定対象領域に含まれるメッシュの節点ｑ（位置）における変位前の位置座標ｘ_１1、ｘ_２1、ｘ_３1に、推定対象領域に含まれるメッシュの節点ｑ（位置）における変位ｕ_ｑ，１、ｕ_ｑ，２、ｕ_ｑ，３を加算することにより、推定対象領域に含まれるメッシュの節点ｑ（位置）における変位後の位置座標ｘ_１２、ｘ_２２、ｘ_３２を導出する。ｕ_ｑ，１、ｕ_ｑ，２、ｕ_ｑ，３の１、２、３は、それぞれ、ｘ_１軸方向成分の値、ｘ_２軸方向成分の値、ｘ_３軸方向成分の値であることを表す。

尚、変位前の位置座標ｘ_１1、ｘ_２1、ｘ_３1は、例えば、静止時または静止時と等価な外力が作用する時刻に導出された位置座標ｘ_１２、ｘ_２２、ｘ_３２である。また、推定対象領域に含まれるメッシュの節点ｑにおける変位前の位置座標ｘ_１1、ｘ_２1、ｘ_３1の初期値は、例えば、台車枠１６の仕様に基づいて定められ、推定装置９００に予め設定される。

＜＜行列導出部４０７ｄ＞＞
行列導出部４０７ｄは、推定対象領域に含まれるメッシュの節点ｑ（位置）における変位前の位置座標ｘ_１1、ｘ_２1、ｘ_３1と、推定対象領域に含まれるメッシュの節点ｑ（位置）における変位後の位置座標ｘ_１２、ｘ_２２、ｘ_３２とに基づいて、行列Ｍの成分ａ_１、ａ_２、ａ_３、ａ_５、ａ_６、ａ_７、ａ_９、ａ_１０、ａ_１１を導出する。尚、このとき、縦ベクトルｖの成分ａ_４、ａ_８、ａ_１２を導出することもできる。本実施形態では、行列Ｍにより変換行列の一例が実現される。

＜＜特異値分解部４０７ｅ＞＞
特異値分解部４０７ｅは、行列導出部４０７ｄにより導出された成分ａ_１、ａ_２、ａ_３、ａ_５、ａ_６、ａ_７、ａ_９、ａ_１０、ａ_１１を有する行列Ｍを特異値分解し、特異値を対角成分として有する対角行列Σを導出する。尚、このとき、直交行列Ｕ、直交行列Ｖの転置行列Ｖ^Ｔも導出される。

＜＜歪み導出部４０７ｆ＞＞
歪み導出部４０７ｆは、特異値分解部４０７ｅにより導出された対角行列Σに基づいて、推定対象領域における主歪みを導出する。本実施形態では、歪み導出部４０７ｆは、特異値分解部４０７ｅにより導出された対角行列Σの対角成分ε_１１、ε_２２、ε_３３を推定対象領域における主歪みとして導出する。以下の説明を簡単にするため、対角行列Σの対角成分ε_１１、ε_２２、ε_３３は、ε_１１≧ε_２２≧ε_３３となるように並べ変えられているものとする。また、歪み導出部４０７ｆは、以下の（４９）式により、推定対象領域における最大主歪みε_ｍａｘを導出する。ここで、特異値分解部４０７ｅにより導出された対角行列Σの対角成分ε_１１、ε_２２、ε_３３のうち、符号を含めて最大の対角成分をε_１１とし、最小の対角成分をε_３３とする。

（４９）式は、対角行列Σの対角成分ε_１１、ε_２２、ε_３３のうち絶対値が最大となるものを、最大主歪みε_ｍａｘとして選択することと等価である。ただし、主歪みは１であるときに歪みが生じていないことを表す（主歪みが１を超えれば引張を表し、１を下回れば圧縮を表す）。そこで、（４９）式においては、直感的に分かりやすいように、対角行列Σの対角成分ε_１１、ε_３３を１から減算した値の絶対値（｜１．０－ε_１１｜、｜１．０－ε_３３｜）に基づいて、対角行列Σの対角成分ε_１１またはε_３３を最大主歪みε_ｍａｘとして選択する場合を例に挙げて示す。

＜＜応力導出部４０７ｇ＞＞
応力導出部４０７ｇは、歪み導出部４０７ｆにより導出された推定対象領域における主歪みに基づいて、推定対象領域における応力を導出する。本実施形態では、応力導出部４０７ｇは、推定対象領域における主歪み（対角行列Σの対角成分）ε_１１、ε_２２、ε_３３に基づいて、推定対象領域における主応力σ_１１、σ_２２、σ_３３を、以下の（５０）式により導出する。

ここで、σ_１１、σ_２２、σ_３３の１１、２２、３３は、それぞれ、せん断応力が０（ゼロ）となる座標系の相互に直交する３つの軸方向成分の値であることを表す。
尚、（５０）式において、対角行列Σの対角成分ε_１１、ε_２２、ε_３３から１を減算しているのは、本実施形態では、主歪みは１を基準とする値であるのに対し、主応力は０（ゼロ）を基準とする値だからである。即ち、本実施形態では、主歪みは１のときに歪みが生じていないことを示すのに対し、主応力は０（ゼロ）のときに応力が生じていないことを示す。尚、主応力は、０（ゼロ）を上回る場合、引張応力であることを示し、０（ゼロ）を下回る場合、圧縮応力であることを示す。
また、λ、μは、それぞれ、ラメ定数、剛性率である（（２２）式、（２３）式を参照）。
応力導出部４０７ｇは、以上のようにして導出した推定対象領域における主応力σ_１１、σ_２２、σ_３３に基づいて、推定対象領域における最大主応力σ_ｍａｘを、以下の（５１）式により導出する。

＜＜応力－位置関係導出部４０８＞＞
第１の実施形態、第２の実施形態、および第３の実施形態では、応力－位置関係導出部４０８は、台車枠１６の応力分布と、鉄道車両の走行位置と、を相互に関連付けて記憶する。本実施形態では、応力－位置関係導出部４０８は、応力導出部４０７ｇにより導出された推定対象領域における主応力σ_１１、σ_２２、σ_３３および最大主応力σ_ｍａｘの少なくとも一方と、データ取得部４０１で取得されたデータに含まれる鉄道車両の走行位置と、を相互に関連付けて記憶する。尚、以下の説明では、応力－位置関係導出部４０８が、応力導出部４０７ｇにより導出された推定対象領域における最大主応力σ_ｍａｘと、鉄道車両の走行位置と、を相互に関連付けて記憶する場合を例に挙げて説明する。本実施形態の説明では、このようにして相互に関連付けられた、推定対象領域における最大主応力σ_ｍａｘと鉄道車両の走行位置を、必要に応じて、応力－位置関係情報と称する。以上のようにして各サンプリング時刻における応力－位置関係情報が得られる。

＜＜判定部４０９＞＞
第１の実施形態、第２の実施形態、および第３の実施形態では、判定部４０９は、応力－位置関係情報（応力の時間変化）に基づいて定められる平均応力および変動応力が、応力限界図の応力許容域内にあるか否かを判定する。これに対し、本実施形態では、判定部４０９は、推定対象領域における最大主応力σ_ｍａｘの絶対値が閾値を上回るか否かを判定する。

＜＜点検箇所特定情報導出部４１０＞＞
第１の実施形態、第２の実施形態、および第３の実施形態では、点検箇所特定情報導出部４１０は、応力限界図にプロットした点が、応力限界図の応力許容域内にない場合に、当該プロットが得られた台車枠１６の位置と、当該プロットした点が応力限界図の応力許容域内にないときの鉄道車両の走行位置とを導出する。これに対し、本実施形態では、点検箇所特定情報導出部４１０は、最大主応力σ_ｍａｘの絶対値が閾値を上回る推定対象領域と、当該推定対象領域における最大主応力σ_ｍａｘの絶対値が閾値を上回るときの鉄道車両の走行位置とを導出する。

＜＜出力部４１１＞＞
出力部４１１は、点検箇所特定情報を出力する。本実施形態では、最大主応力σ_ｍａｘの絶対値が閾値を上回る推定対象領域と、当該推定対象領域における最大主応力σ_ｍａｘの絶対値が閾値を上回るときの鉄道車両の走行位置とを示す情報が点検箇所特定情報となる。

また、出力部４１１は、応力－位置関係情報に代えてまたは加えて、歪み－位置関係情報を出力してもよい。歪み－位置関係情報は、歪み導出部４０７ｆにより導出された主歪みε_１１、ε_２２、ε_３３および最大主歪みε_ｍａｘの少なくとも一方と、鉄道車両の走行位置と、を相互に関連付けた情報である。

＜フローチャート＞
図１０のフローチャートを参照しながら、本実施形態の推定装置９００を用いた推定方法の一例を説明する。図１０のフローチャートのステップＳ６０１～Ｓ６０７、Ｓ１００１～Ｓ１０１０の繰り返し処理は、サンプリング時刻が到来する度に実行されるものとする。
まず、ステップＳ６０１～Ｓ６０７の処理は、第１の実施形態で説明した図６のフローチャートによる処理と同じである。ステップＳ６０７において、変位分布導出部４０７ａにより、実変位ベクトル｛ｕ｝が、台車枠１６の変位分布として導出される。ステップＳ６０７の処理の後、ステップＳ１００１の処理が実行される。

ステップＳ１００１において、位置座標導出部４０７ｃは、推定対象領域に含まれるメッシュの節点ｑ（位置）における変位前の位置座標ｘ_１1、ｘ_２1、ｘ_３1に、推定対象領域に含まれるメッシュの節点ｑにおける変位ｕ_ｑ，１、ｕ_ｑ，２、ｕ_ｑ，３を加算することにより、推定対象領域に含まれるメッシュの節点（位置）における変位後の位置座標ｘ_１２、ｘ_２２、ｘ_３２を導出する。尚、変位前の位置座標ｘ_１1、ｘ_２1、ｘ_３1は、静止時の位置座標ｘ_１２、ｘ_２２、ｘ_３２である。

次に、ステップＳ１００２において、行列導出部４０７ｄは、推定対象領域に含まれるメッシュの節点ｑ（位置）における変位前の位置座標ｘ_１1、ｘ_２1、ｘ_３1と、ステップＳ１００１で導出された推定対象領域に含まれるメッシュの節点ｑ位置）における変位後の位置座標ｘ_１２、ｘ_２２、ｘ_３２とに基づいて、行列Ｍの成分ａ_１、ａ_２、ａ_３、ａ_５、ａ_６、ａ_７、ａ_９、ａ_１０、ａ_１１を導出する。

次に、ステップＳ１００３において、特異値分解部４０７ｅは、ステップＳ１００２で導出された成分ａ_１、ａ_２、ａ_３、ａ_５、ａ_６、ａ_７、ａ_９、ａ_１０、ａ_１１を有する行列Ｍを特異値分解し、特異値を対角成分として有する対角行列Σを導出する。
次に、ステップＳ１００４において、歪み導出部４０７ｆは、ステップＳ１００３で導出された対角行列Σの対角成分ε_１１、ε_２２、ε_３３を推定対象領域における主歪みとして導出する。そして、歪み導出部４０７ｆは、（４９）式により、推定対象領域における最大主歪みε_ｍａｘを導出する。

次に、ステップＳ１００５において、応力導出部４０７ｇは、ステップＳ１００４で導出された推定対象領域における主歪み（対角行列Σの対角成分）ε_１１、ε_２２、ε_３３に基づいて、推定対象領域における主応力σ_１１、σ_２２、σ_３３を、（５０）式により導出する。そして、応力導出部４０７ｇは、推定対象領域における主応力σ_１１、σ_２２、σ_３３に基づいて、推定対象領域における最大主応力σ_ｍａｘを、（５１）式により導出する。

次に、ステップＳ１００６において、応力－位置関係導出部４０８は、推定対象領域における主応力σ₁₁、σ₂₂、σ₃₃および最大主応力σ_maxと、鉄道車両の走行位置と、を相互に関連付けて記憶する（応力－位置関係情報を記憶する）。
次に、ステップＳ１００７において、判定部４０９は、推定対象領域における最大主応力σ_maxの絶対値（｜σ_max｜）が閾値を上回るか否かを判定する。この判定の結果、推定対象領域における最大主応力σ_maxの絶対値が閾値を上回らない場合、処理は、ステップＳ１００８を省略して後述するステップＳ１００９に進む。

一方、推定対象領域における最大主応力σ_maxの絶対値が閾値を上回る場合、処理は、ステップＳ１００８に進む。処理がステップＳ１００８に進むと、点検箇所特定情報導出部４１０は、最大主応力σ_maxの絶対値が閾値を上回る推定対象領域と、当該推定対象領域における最大主応力σ_maxの絶対値が閾値を上回るときの鉄道車両の走行位置と、を、点検箇所特定情報として導出する。
次に、ステップＳ１００９において、推定装置９００は、測定を終了するか否かを判定する。ステップＳ１００９の処理は、例えば、第１の実施形態で説明した図６のステップＳ６１３の処理と同じである。

ステップＳ１００９の判定の結果、測定を終了しない場合、処理は、ステップＳ６０１に戻る。そして、次のサンプリング時刻におけるステップＳ６０１～Ｓ６０６、およびＳ１００１～Ｓ１００９の処理が実行される。
一方、ステップＳ１００９の判定の結果、測定を終了する場合、処理は、ステップＳ１０１０に進む。処理がステップＳ１０１０に進むと、出力部４１１は、点検箇所特定情報を出力する。そして、図１０のフローチャートによる処理が終了する。
尚、ステップＳ１０１０の処理を、ステップＳ１００８の次に行ってもよい。この場合、点検箇所特定情報は、サンプリング時刻の単位で出力される。

＜まとめ＞
以上のように本実施形態では、推定装置９００は、推定対象領域における変位前の位置座標ｘ_１1、ｘ_２1、ｘ_３1と、推定対象領域における変位後の位置座標ｘ_１２、ｘ_２２、ｘ_３２とに基づいて、行列Ｍの成分ａ_１、ａ_２、ａ_３、ａ_５、ａ_６、ａ_７、ａ_９、ａ_１０、ａ_１１を導出する。行列Ｍは、アフィン変換の際に変位前の位置座標ｘ_１1、ｘ_２1、ｘ_３1に乗算される行列である。推定装置９００は、行列Ｍを特異値分解し、特異値を対角成分として有する対角行列Σを導出する。推定装置９００は、対角行列Σに基づいて、推定対象領域における主歪みε_１1、ε_２２、ε_３３を導出する。従って、実変位ベクトル｛ｕ｝に基づいて台車枠１６の主歪みおよび主応力を導出することを、ＦＥＭの形状関数を用いずに実行することができる。また、歪ゲージで歪みを測定しなくても（即ち、歪みそのものの測定値を得なくても）、歪みを導出することができる。よって、歪ゲージで歪みを測定しなくても、外力が作用する物体の状態を推定することができる。尚、本実施形態では、実変位ベクトル｛ｕ｝に含まれる変位ｕ_ｑ，ｉを示すメッシュの節点ｑに基づいて推定対象領域を定める。従って、推定対象領域として様々な領域を設定することができる。

また、本実施形態では、推定装置９００は、対角行列Σの成分ε_１1、ε_２２、ε_３３に基づいて推定対象領域における最大主歪みε_ｍａｘを導出する。従って、外力が作用する物体の状態を示す情報としてより多くの情報を導出することができる。

また、本実施形態では、推定装置９００は、対角行列Σの成分ε_１1、ε_２２、ε_３３に基づいて推定対象領域における主応力σ_１1、σ_２２、σ_３３を導出し、推定対象領域における主応力σ_１1、σ_２２、σ_３３に基づいて、最大主応力σ_ｍａｘを導出する。従って、外力が作用する物体の状態を示す情報としてより一層多くの情報を導出することができる。

＜変形例＞
本実施形態においても、第１の実施形態、第２の実施形態、および第３の実施形態で説明した種々の変形例を採用することができる。また、第２の実施形態または第３の実施形態に対して本実施形態の手法を適用してもよい。

（ハードウェア）
推定装置４００のハードウェアの一例について説明する。推定装置８００、９００のハードウェアは、推定装置４００のハードウェアと同じもので実現することができる。従って、ここでは、推定装置８００、９００のハードウェアの詳細な説明を省略する。図１１において、推定装置４００は、ＣＰＵ１１０１、主記憶装置１１０２、補助記憶装置１１０３、通信回路１１０４、信号処理回路１１０５、画像処理回路１１０６、Ｉ／Ｆ回路１１０７、ユーザインターフェース１１０８、ディスプレイ１１０９、およびバス１１１０を有する。

ＣＰＵ１１０１は、推定装置４００の全体を統括制御する。ＣＰＵ１１０１は、主記憶装置１１０２をワークエリアとして用いて、補助記憶装置１１０３に記憶されているプログラムを実行する。主記憶装置１１０２は、データを一時的に格納する。補助記憶装置１１０３は、ＣＰＵ１１０１によって実行されるプログラムの他、各種のデータを記憶する。

通信回路１１０４は、推定装置４００の外部との通信を行うための回路である。通信回路１１０４は、推定装置４００の外部と無線通信を行っても有線通信を行ってもよい。通信回路１１０４は、無線通信を行う場合、鉄道車両に設けられるアンテナに接続される。

信号処理回路１１０５は、通信回路１１０４で受信された信号や、ＣＰＵ１１０１による制御に従って入力した信号に対し、各種の信号処理を行う。
画像処理回路１１０６は、ＣＰＵ１１０１による制御に従って入力した信号に対し、各種の画像処理を行う。この画像処理が行われた信号は、ディスプレイ１１０９に出力される。
ユーザインターフェース１１０８は、オペレータが推定装置４００に対して指示を行う部分である。ユーザインターフェース１１０８は、例えば、ボタン、スイッチ、およびダイヤル等を有する。また、ユーザインターフェース１１０８は、ディスプレイ１１０９を用いたグラフィカルユーザインターフェースを有していてもよい。

ディスプレイ１１０９は、画像処理回路１１０６から出力された信号に基づく画像を表示する。Ｉ／Ｆ回路１１０７は、Ｉ／Ｆ回路１１０７に接続される装置との間でデータのやり取りを行う。図１１では、Ｉ／Ｆ回路１１０７に接続される装置として、ユーザインターフェース１１０８およびディスプレイ１１０９を示す。しかしながら、Ｉ／Ｆ回路１１０７に接続される装置は、これらに限定されない。例えば、可搬型の記憶媒体がＩ／Ｆ回路１１０７に接続されてもよい。また、ユーザインターフェース１１０８の少なくとも一部およびディスプレイ１１０９は、推定装置４００の外部にあってもよい。
出力部４１１は、例えば、通信回路１１０４および信号処理回路１１０５と、画像処理回路１１０６、Ｉ／Ｆ回路１１０７、およびディスプレイ１１０９との少なくとも何れか一方を用いることにより実現される。

尚、ＣＰＵ１１０１、主記憶装置１１０２、補助記憶装置１１０３、信号処理回路１１０５、画像処理回路１１０６、およびＩ／Ｆ回路１１０７は、バス１１１０に接続される。これらの構成要素間の通信は、バス１１１０を介して行われる。また、推定装置４００、８００、９００のハードウェアは、前述した推定装置４００、８００、９００の機能を実現することができれば、図１１に示すものに限定されない。

（実施例）
次に、実施例を説明する。本実施例では、実際の鉄道車両の台車枠をモデル化し、第１の実施形態および第２の実施形態で説明した手法のそれぞれにおいて、台車枠の応力分布を導出した。本実施例では、台車枠は、ばね帽およびモノリンクから外力を受けるものとした。また、第２の実施形態で説明した手法では、台車枠に複数の加速度センサを取り付けた。導出した応力分布を調べた結果、平均応力および変動応力を応力限界図にプロットした点が応力限界図の応力許容域内に入らない箇所が存在した。当該プロットした点が応力限界図の応力許容域内にないときの鉄道車両の走行位置に対応する軌道を点検したところ、軌道不整が発生していることが判明した。

図１２は、台車枠の空気バネ座の或る位置における変位と時間との関係を示す図である。図１２において、グラフ１２０１は、第１の実施形態で説明した手法で導出した実変位ベクトル｛ｕ｝の当該位置に対応する成分から得られたものである。グラフ１２０２は、当該位置に取り付けた加速度センサの測定値から得られたものである。図１２に示すように、第１の実施形態で説明した手法により、台車枠の変位を正確に導出することができることが分かる。

図１３Ａおよび図１３Ｂは、比較的大きな応力が発生した台車枠の或る位置における変位と時間との関係を示す図である。図１３Ａにおいて、グラフ１３０１は、第１の実施形態で説明した手法で導出した実変位ベクトル｛ｕ｝の当該位置に対応する成分から得られたものである。グラフ１３０２は、当該位置に取り付けた加速度センサの測定値から得られたものである。尚、図１３Ａにおいて、グラフ１３０１は、濃度の薄いグラフである。図１３Ａにおいて、グラフ１３０２は、濃度の濃いグラフである。図１３Ｂにおいて、グラフ１３０３は、第２の実施形態で説明した手法で導出した実変位ベクトル｛ｕ｝の当該位置に対応する成分から得られたものである。図１３Ｂにおいて、グラフ１３０３は、濃度の薄いグラフである。図１３Ａに示すように、第１の実施形態で説明した手法で導出した変位の時間変化は、絶対値が実測値と乖離している部分があるものの、変化の傾向は、実測値と対応する。第１の実施形態の手法では図１３Ａに示すような変位が導出される場合、図１３Ｂに示すように、第２の実施形態の手法を採用することにより、台車枠の変位をより正確に導出することができることが分かる。尚、第２の実施形態の手法で変位を導出するに際し、当該位置に取り付けた加速度センサの測定値は使用しなかった。

また、本実施例では、実際の鉄道車両の台車枠をモデル化し、第１の実施形態および第２の実施形態で説明したモード行列［φ］と、第３の実施形態で説明したモード行列［φ］とを比較した。
本実施例では、有限要素法のメッシュの節点ｑ_ｊの数を４０００００とした。４０００００個のメッシュの節点ｑ_ｊのうち、モード行列［φ］の成分の近似式に用いるメッシュの節点ｑ_ｊの数ｍを１０００とした。ｍ（＝１０００）個の節点の位置を、台車枠に対してランダムに設定した。ｍ個の節点の位置とｍ’個の節点の位置とを全て同じにした。

図１４は、モード行列［φ］の成分φ_ｑ，ｉ ^（ｐ）とモード行列［φ］の成分の近似値φ’_ｓ，ｉ ^（ｐ）（ｘ_ｓ，１，ｘ_ｓ，２，ｘ_ｓ，３）との関係を示す図である。図１４では、１０次の固有振動モード（ｐ＝１０、６０Ｈｚ）およびｘ_１軸方向（ｉ＝１）に対応するモード行列［φ］の成分を例に挙げて示す。図１４では、モード行列［φ］の成分φ_ｑ，ｉ ^（ｐ）をモード行列の値と表記する。また、モード行列［φ］の成分の近似値φ’_ｓ，ｉ ^（ｐ）（ｘ_ｓ，１，ｘ_ｓ，２，ｘ_ｓ，３）をモード行列の近似値と表記する。図１４では、モード行列の値とモード行列の近似値とから定まる点を薄い濃度で示す。

図１４に示すように、第３の実施形態で説明したようにして導出されるモード行列［φ］の成分の近似値と、第１の実施形態および第２の実施形態で説明したモード行列［φ］の成分との平均誤差は１．２％であり、標準偏差は０．００２７（モード行列の値の標準偏差に対し２.０％）である。従って、第３の実施形態の手法により、モード行列［φ］を高精度に近似することができることが分かる。

また、本実施例では、実際の鉄道車両の台車枠１６をモデル化し、第４の実施形態で説明した手法により、台車枠１６の推定対象領域における最大主歪みε_ｍａｘおよび最大主応力σ_ｍａｘを導出した。

鉄道車両を走行させた場合に台車枠１６に発生する応力を有限要素法により導出した。台車枠１６に発生する応力に基づいて、最大主応力が発生する箇所を特定し、最大主応力が発生する箇所を基準点とした。鉄道車両を走行させた場合の基準点における変位と、基準点の周囲の１２箇所の点における変位とを、実変位ベクトル｛ｕ｝に基づいて導出した。ここでは、第１の実施形態の手法で実変位ベクトル｛ｕ｝を導出した。このようにして導出した変位に基づいて、行列Ｍの成分をａ₁、ａ₂、ａ₃、ａ₅、ａ₆、ａ₇、ａ₉、ａ₁₀、ａ₁₁を導出し、行列Ｍを特異値分解して対角行列Σを導出した。そして、対角行列Σに基づいて、基準点と基準点の周囲の１２箇所の点とを含む推定対象領域における最大主歪みε_maxおよび最大主応力σ _maxを導出した。尚、本実施例では、ヤング率Ｅを２０５．９ＧＰａとし、ポアソン比νを０．３とした。
また、鉄道車両を走行させた場合の基準点における最大主応力を基準点に取り付けた歪ゲージの測定値に基づいて導出した。

図１５は、推定対象領域における最大主歪みε_ｍａｘと時間との関係の一例を示す図である。図１６は、推定対象領域における最大主応力σ_ｍａｘと時間との関係の一例を示す図である。図１６において、推定値は、第４の実施形態の手法で導出した推定対象領域における最大主応力σ_ｍａｘと時間との関係を示す。測定値は、基準点に取り付けた歪ゲージの測定値に基づく最大主応力σ_ｍａｘと時間との関係を示す。図１６では、推定値は濃い濃度で示され、測定値は薄い濃度で示される。図１７は、第４の実施形態の手法で導出した推定対象領域における最大主応力σ_ｍａｘと、基準点に取り付けた歪ゲージの測定値に基づく最大主応力σ_ｍａｘとの関係の一例を示す図である。図１７は、図１６において、同一時刻の推定値および測定値の組をプロットすることにより得られる。図１７では、第４の実施形態の手法で導出した推定対象領域における最大主応力σ_ｍａｘを最大主応力の推定値と表記する。基準点に取り付けた歪ゲージの測定値に基づく最大主応力σ_ｍａｘを最大主応力の測定値と表記する。図１７では、最大主応力の推定値と最大主応力の測定値とから定まる点を薄い濃度で示す。

図１７において、基準点に取り付けた歪ゲージの測定値に基づく最大主応力σ_ｍａｘ（測定値）と、第４の実施形態の手法で導出した推定対象領域における最大主応力σ_ｍａｘ（推定値）との平均誤差は１．４ＭＰａであり、標準偏差は１．５ＭＰａである。従って、第４の実施形態の手法により、歪ゲージを用いた場合と同等の結果が得られることが分かる。

（その他の変形例）
尚、以上説明した本発明の実施形態は、コンピュータがプログラムを実行することによって実現することができる。また、前記プログラムを記録したコンピュータ読み取り可能な記録媒体及び前記プログラム等のコンピュータプログラムプロダクトも本発明の実施形態として適用することができる。記録媒体としては、例えば、フレキシブルディスク、ハードディスク、光ディスク、光磁気ディスク、ＣＤ－ＲＯＭ、磁気テープ、不揮発性のメモリカード、ＲＯＭ等を用いることができる。
また、以上説明した本発明の実施形態は、何れも本発明を実施するにあたっての具体化の例を示したものに過ぎず、これらによって本発明の技術的範囲が限定的に解釈されてはならないものである。すなわち、本発明はその技術思想、またはその主要な特徴から逸脱することなく、様々な形で実施することができる。

本発明は、例えば、鉄道車両の状態を推定することに利用することができる。

Claims (20)

1. 台車枠と、前記台車枠に直接または他の部品を介して接続された接続部品と、を有する鉄道車両の状態を推定する推定装置であって、
   前記接続部品に取り付けられた第１のセンサで、前記鉄道車両の走行中に測定された第１の測定データを取得するデータ取得手段と、
   前記第１の測定データに基づいて、前記台車枠の着力箇所に作用する外力を導出し、前記台車枠の着力箇所に作用する外力に基づいて、前記台車枠に作用する外力の分布を導出する外力導出手段と、
   前記台車枠に作用する外力の分布に基づいて、前記台車枠の状態を表す情報を導出する状態導出手段と、を有し、
   前記第１のセンサは、前記台車枠の着力箇所に作用する外力を導出することが可能な物理量を測定するためのセンサであり、
   前記台車枠の着力箇所は、前記台車枠において、前記接続部品により発生する外力が作用する箇所であり、
   前記状態導出手段は、前記台車枠の振動を表す運動方程式を解くことにより、前記台車枠の変位の分布を前記台車枠の状態を表す情報として導出する変位分布導出手段を少なくとも有し、
   前記台車枠の振動を表す運動方程式における外力の分布には、前記外力導出手段により導出された前記台車枠に作用する外力の分布が与えられることを特徴とする推定装置。
2. 前記データ取得手段は、前記台車枠に取り付けられた第２のセンサで測定された第２の測定データを更に取得し、
   前記変位分布導出手段は、前記台車枠の振動を表す運動方程式を解くことにより導出した変位の分布を補正するための補正パラメータを導出し、前記台車枠の振動を表す運動方程式を解くことにより導出した前記台車枠の変位の分布と、前記補正パラメータとを用いて、前記台車枠の変位の分布を導出し、
   前記第２のセンサは、前記台車枠の所定の位置の変位を導出することが可能な物理量を測定するためのセンサであり、
   前記変位分布導出手段は、前記第２の測定データに基づいて導出された前記台車枠の所定の位置の変位と、前記台車枠の振動を表す運動方程式を解くことにより導出した前記台車枠の変位の分布のうち前記所定の位置の変位と、を用いて、前記補正パラメータを導出することを特徴とする請求項１に記載の推定装置。
3. 前記外力導出手段は、前記第１の測定データに基づいて、物理座標系において前記台車枠に作用する外力の分布を導出し、当該台車枠に作用する外力の分布を、モード行列を用いて、モード座標系における前記台車枠に作用する外力の分布に変換し、
   前記モード行列は、物理座標系の前記台車枠の振動を表す運動方程式に対する固有値解析の結果に基づいて導出され、
   前記変位分布導出手段は、モード座標系における前記台車枠の振動を表す運動方程式を解くことにより、モード座標系における前記台車枠の変位の分布を導出し、当該モード座標系における前記台車枠の変位の分布を、前記モード行列を用いて、物理座標系における前記台車枠の変位の分布に変換し、
   モード座標系における前記台車枠の振動を表す運動方程式における外力の分布には、前記外力導出手段により導出されたモード座標系における前記台車枠に作用する外力の分布が与えられることを特徴とする請求項１または２に記載の推定装置。
4. 前記データ取得手段は、前記台車枠に取り付けられた第２のセンサで測定された第２の測定データを更に取得し、
   前記変位分布導出手段は、モード座標系における前記台車枠の振動を表す運動方程式を解くことにより導出した変位の分布を補正するための補正パラメータを導出し、モード座標系における前記台車枠の振動を表す運動方程式を解くことにより導出した前記台車枠の変位の分布を、前記補正パラメータと、前記モード行列とを用いて、物理座標系における前記台車枠の変位の分布に変換し、
   前記第２のセンサは、前記台車枠の所定の位置の変位を導出することが可能な物理量を測定するためのセンサであり、
   前記変位分布導出手段は、前記第２の測定データに基づいて導出された前記台車枠の所定の位置の変位と、モード座標系における前記台車枠の振動を表す運動方程式を解くことにより導出したモード座標系における前記台車枠の変位の分布のうち前記所定の位置の変位と、前記モード行列と、を用いて、前記補正パラメータを導出することを特徴とする請求項３に記載の推定装置。
5. 物理座標系の前記台車枠の振動を表す運動方程式に対する固有値解析を行うことにより前記モード行列を導出する固有値解析手段を更に有することを特徴とする請求項３または４に記載の推定装置。
6. 前記モード行列の成分の近似式を用いて、前記モード行列の成分の近似値を導出する近似値導出手段を更に有し、
   前記モード行列の成分の近似式は、前記固有値解析により導出された固有ベクトルに基づいて導出される式であり、
   前記モード行列の成分の近似式は、前記モード行列の成分の近似値として、固有振動モードの次数、位置、および位置を定める座標系の自由度に応じた近似値を計算する式であり、
   前記モード行列の成分には、前記近似値導出手段により導出された前記モード行列の成分の近似値が含まれることを特徴とする請求項３～５の何れか１項に記載の推定装置。
7. 前記モード行列の成分の近似式は、前記固有値解析により導出された固有ベクトルに基づいて、カーネルリッジ回帰分析を行うことにより導出される式であることを特徴とする請求項６に記載の推定装置。
8. 前記状態導出手段は、前記変位分布導出手段により導出された前記台車枠の変位の分布に基づいて、前記台車枠の応力の分布を前記台車枠の状態を表す情報として導出する応力分布導出手段を更に有することを特徴とする請求項１～７の何れか１項に記載の推定装置。
9. 前記状態導出手段は、前記台車枠の少なくとも３箇所の点における変位前の位置座標と、前記変位分布導出手段により導出された前記台車枠の変位のうち、前記台車枠の少なくとも３箇所の点における変位とに基づいて、前記台車枠の少なくとも３箇所の点における変位後の位置座標を導出する位置座標導出手段と、
   前記台車枠の少なくとも３箇所の点における変位前の位置座標と、前記位置座標導出手段により導出された前記台車枠の少なくとも３箇所の点における変位後の位置座標とに基づいて、変換行列の成分を導出する行列導出手段と、
   前記行列導出手段により導出された前記成分を有する前記変換行列を特異値分解する特異値分解手段と、
   前記特異値分解手段により特異値分解を実行することにより導出された対角行列の成分である特異値を、前記台車枠の主歪みとして導出する歪み導出手段と、を更に有し、
   前記変換行列は、前記台車枠の少なくとも３箇所の点における変位前の位置座標をアフィン変換して前記台車枠の少なくとも３箇所の点における変位後の位置座標を導出する際に前記台車枠の少なくとも３箇所の点における変位前の位置座標に乗算される行列であることを特徴とする請求項１～７の何れか１項に記載の推定装置。
10. 前記歪み導出手段は、前記対角行列に基づいて、前記台車枠の最大主歪みを更に導出することを特徴とする請求項９に記載の推定装置。
11. 前記状態導出手段は、前記歪み導出手段により導出された前記台車枠の主歪みに基づいて、前記台車枠の主応力を導出する応力導出手段を更に有することを特徴とする請求項９または１０に記載の推定装置。
12. 前記応力導出手段は、前記台車枠の主応力に基づいて、前記台車枠の最大主応力を更に導出することを特徴とする請求項１１に記載の推定装置。
13. 前記少なくとも３箇所の点は、少なくとも４箇所の点であり、
    前記位置座標は、３軸の座標系における位置座標であることを特徴とする請求項９～１２の何れか１項に記載の推定装置。
14. 前記台車枠および軌道の少なくとも一方の点検箇所を特定するための情報である点検箇所特定情報を導出する点検箇所特定情報導出手段を更に有し、
    前記データ取得手段は、前記鉄道車両の走行位置を示す走行位置データを更に取得し、
    前記点検箇所特定情報導出手段は、前記状態導出手段により導出された前記台車枠の状態を表す情報に基づく所定の指標値と、前記走行位置データと、に基づいて、前記点検箇所特定情報を導出することを特徴とする請求項１～１３の何れか１項に記載の推定装置。
15. 前記鉄道車両は、軸箱を更に有し、
    前記接続部品は、前記軸箱であることを特徴とする請求項１～１４の何れか１項に記載の推定装置。
16. 前記第１のセンサは、加速度センサであることを特徴とする請求項１～１５の何れか１項に記載の推定装置。
17. 前記外力導出手段は、前記第１のセンサが取り付けられた前記鉄道車両と同一の前記鉄道車両の前記台車枠の着力箇所に作用する外力を、当該第１のセンサで前記鉄道車両の走行中に測定された第１の測定データに基づいて導出することを特徴とする請求項１～１６の何れか１項に記載の推定装置。
18. 前記外力導出手段は、前記第１のセンサが取り付けられた前記鉄道車両と異なる前記鉄道車両の前記台車枠の着力箇所に作用する外力を、当該第１のセンサで前記鉄道車両の走行中に測定された第１の測定データに基づいて導出し、
    前記第１のセンサが取り付けられた前記鉄道車両と異なる前記鉄道車両は、前記第１のセンサが取り付けられた前記鉄道車両と同一の編成の鉄道車両であることを特徴とする請求項１～１７の何れか１項に記載の推定装置。
19. 台車枠と、前記台車枠に直接または他の部品を介して接続された接続部品と、を有する鉄道車両の状態を推定する推定方法であって、
    前記接続部品に取り付けられた第１のセンサで前記鉄道車両の走行中に測定された第１の測定データを取得するデータ取得工程と、
    前記第１の測定データに基づいて、前記台車枠の着力箇所に作用する外力を導出し、前記台車枠の着力箇所に作用する外力に基づいて、前記台車枠に作用する外力の分布を導出する外力導出工程と、
    前記台車枠に作用する外力の分布に基づいて、前記台車枠の状態を表す情報を導出する状態導出工程と、を有し、
    前記第１のセンサは、前記台車枠の着力箇所に作用する外力を導出することが可能な物理量を測定するためのセンサであり、
    前記台車枠の着力箇所は、前記台車枠において、前記接続部品により発生する外力が作用する箇所であり、
    前記状態導出工程は、前記台車枠の振動を表す運動方程式を解くことにより、前記台車枠の変位の分布を前記台車枠の状態を表す情報として導出する変位分布導出工程を少なくとも有し、
    前記台車枠の振動を表す運動方程式における外力の分布には、前記外力導出工程により導出された前記台車枠に作用する外力の分布が与えられることを特徴とする推定方法。
20. 台車枠と、前記台車枠に直接または他の部品を介して接続された接続部品と、を有する鉄道車両の状態を推定するための処理をコンピュータに実行させるためのプログラムであって、
    前記接続部品に取り付けられた第１のセンサで前記鉄道車両の走行中に測定された第１の測定データを取得するデータ取得工程と、
    前記第１の測定データに基づいて、前記台車枠の着力箇所に作用する外力を導出し、前記台車枠の着力箇所に作用する外力に基づいて、前記台車枠に作用する外力の分布を導出する外力導出工程と、
    前記台車枠に作用する外力の分布に基づいて、前記台車枠の状態を表す情報を導出する状態導出工程と、をコンピュータに実行させ、
    前記第１のセンサは、前記台車枠の着力箇所に作用する外力を導出することが可能な物理量を測定するためのセンサであり、
    前記台車枠の着力箇所は、前記台車枠において、前記接続部品により発生する外力が作用する箇所であり、
    前記状態導出工程は、前記台車枠の振動を表す運動方程式を解くことにより、前記台車枠の変位の分布を前記台車枠の状態を表す情報として導出する変位分布導出工程を少なくとも有し、
    前記台車枠の振動を表す運動方程式における外力の分布には、前記外力導出工程により導出された前記台車枠に作用する外力の分布が与えられることを特徴とするプログラム。

Images