JP6939540B2 - 接触角推定システム、接触角推定方法、およびプログラム - Google Patents

Description

本発明は、接触角推定システム、接触角推定方法、およびプログラムに関し、特に、鉄道車両における接触角を推定するために用いて好適なものである。

鉄道車両の安全性を評価する指標として、脱線係数がある。脱線係数は、輪重Ｐに対する横圧Ｑの比で表される。脱線係数と所定の評価値とを比較することで、鉄道車両が乗り上がり脱線する可能性を見積もることができる。脱線係数は、鉄道車両の脱線現象を判定する重要な指標となるものの、車輪と軌条（レール）との接触状態をより詳細に評価するための指標として、脱線係数に大きく影響を与える接触角を求めるニーズがある。接触角は、車輪と軌条との接触位置での接平面と水平面とのなす角度で表される。尚、車輪と軌条とが点接触している場合には、接触位置は、当該点の位置であり、車輪と軌条とが面接触している場合には、接触位置は、当該面の重心の位置とすればよい。以下の説明では、接触角は、車輪と軌条との接触位置での接平面と水平面とのなす角度のうち小さい方の角度（鋭角）であるものとして説明を行う。

接触角を求める技術として特許文献１に記載の技術がある。特許文献１には、１つの輪軸を構成する左右それぞれの車輪の踏面の形状データと、軌条の断面の形状データとを離散値データに変換し、それらに基づいて、接触角を算出することが記載されている。

特許第４６５７７６７号公報 国際公開第２０１７／１６４１３３号公報

しかしながら、特許文献１に記載の技術では、車輪と軌条の形状を計測しなければならない。車輪の形状と、当該車輪が接触している軌条の形状とを同時に計測することは容易ではない。この点について特許文献１には、この計測の具体的な手法についての開示はない。このように、特許文献１に記載の技術では、鉄道車両における接触角を推定することは容易ではない。

本発明は、以上のような問題点に鑑みてなされたものであり、鉄道車両における接触角を容易に推定することができるようにすることを目的とする。

本発明の接触角推定システムは、車体と台車と輪軸とを有する鉄道車両の車輪と軌条との接触位置における接平面と左右方向とのなす角度である接触角を推定する接触角推定システムであって、前記鉄道車両を軌道上で走行させることにより測定される測定値のデータである計測データを取得するデータ取得手段と、前記データ取得手段により取得された前記計測データと、前記接触角を未知数とする方程式からなる第１のフィルタとを用いて、前記接触角を算出する接触角算出手段と、を有し、前記計測データは、前記左右方向の一方側および他方側の前記車輪のそれぞれにおける輪重および横圧の測定値のデータと、前記鉄道車両の走行速度の測定値のデータと、前後方向力の測定値のデータとを含み、前記前後方向力は、前記輪軸と、当該輪軸が設けられる前記台車との間に配置される部材に生じる前後方向の力であり、前記部材は、軸箱を支持するための部材であり、前記前後方向は、前記鉄道車両の走行方向に沿う方向であり、前記左右方向は、前記前後方向と、前記軌道に対し垂直な方向である上下方向との双方に垂直な方向であり、前記方程式は、前記接触角の関数であるコスト関数の更新前の前記接触角における微分値と、前記コスト関数の更新前の前記接触角における値と、緩和係数と、更新前の前記接触角と、更新後の前記接触角との関係を表す式を含み、前記コスト関数は、前記車輪と前記軌条との接触位置での法線荷重の前記左右方向の成分の、前記左右方向の一方側および他方側の前記車輪のそれぞれにおける値と、前記車輪と前記軌条との接触位置での横クリープ力の前記左右方向の成分の、前記左右方向の一方側および他方側の前記車輪のそれぞれにおける値とを加算した値と、前記車輪と前記軌条との接触位置での前記左右方向の一方側および他方側の前記車輪のそれぞれにおける前記横圧を加算した値と、の差を含む第１のコスト関数を含み、前記方程式における前記接触角は、前記左右方向の一方側の前記車輪と前記軌条との接触位置における接平面と前記左右方向とのなす角度である第１の接触角と、前記左右方向の他方側の前記車輪と前記軌条との接触位置における接平面と前記左右方向とのなす角度である第２の接触角とを含み、前記方程式において未知数となる前記接触角は、前記第１の接触角と前記第２の接触角との少なくとも何れか一方であり、前記横圧、前記輪重、前記法線荷重、および前記横クリープ力は、前記車輪と前記軌条との接触位置における値であり、前記第１のコスト関数における前記横クリープ力は、前記横クリープ力を用いずに、前記輪軸の前記左右方向の速度と、前記台車のヨーイング方向における角変位と、前記前後方向力と、前記鉄道車両の走行速度と、を用いて表現され、前記ヨーイング方向は、前記上下方向を回動軸とする回動方向であり、前記接触角算出手段は、前記第１のフィルタに含まれる前記方程式に前記更新前の接触角を与えて前記更新後の前記接触角を算出することを、所定の条件が成立するまで行うことを特徴とする。

本発明の接触角推定方法は、車体と台車と輪軸とを有する鉄道車両の車輪と軌条との接触位置における接平面と左右方向とのなす角度である接触角を推定する接触角推定方法であって、前記鉄道車両を軌道上で走行させることにより測定される測定値のデータである計測データを取得するデータ取得工程と、前記データ取得工程により取得された前記計測データと、前記接触角を未知数とする方程式からなる第１のフィルタとを用いて、前記接触角を算出する接触角算出工程と、を有し、前記計測データは、前記左右方向の一方側および他方側の前記車輪のそれぞれにおける輪重および横圧の測定値のデータと、前記鉄道車両の走行速度の測定値のデータと、前後方向力の測定値のデータとを含み、前記前後方向力は、前記輪軸と、当該輪軸が設けられる前記台車との間に配置される部材に生じる前後方向の力であり、前記部材は、軸箱を支持するための部材であり、前記前後方向は、前記鉄道車両の走行方向に沿う方向であり、前記左右方向は、前記前後方向と、前記軌道に対し垂直な方向である上下方向との双方に垂直な方向であり、前記方程式は、前記接触角の関数であるコスト関数の更新前の前記接触角における微分値と、前記コスト関数の更新前の前記接触角における値と、緩和係数と、更新前の前記接触角と、更新後の前記接触角との関係を表す式を含み、前記コスト関数は、前記車輪と前記軌条との接触位置での法線荷重の前記左右方向の成分の、前記左右方向の一方側および他方側の前記車輪のそれぞれにおける値と、前記車輪と前記軌条との接触位置での横クリープ力の前記左右方向の成分の、前記左右方向の一方側および他方側の前記車輪のそれぞれにおける値とを加算した値と、前記車輪と前記軌条との接触位置での前記左右方向の一方側および他方側の前記車輪のそれぞれにおける前記横圧を加算した値と、の差を含む第１のコスト関数を含み、前記方程式における前記接触角は、前記左右方向の一方側の前記車輪と前記軌条との接触位置における接平面と前記左右方向とのなす角度である第１の接触角と、前記左右方向の他方側の前記車輪と前記軌条との接触位置における接平面と前記左右方向とのなす角度である第２の接触角とを含み、前記方程式において未知数となる前記接触角は、前記第１の接触角と前記第２の接触角との少なくとも何れか一方であり、前記横圧、前記輪重、前記法線荷重、および前記横クリープ力は、前記車輪と前記軌条との接触位置における値であり、前記第１のコスト関数における前記横クリープ力は、前記横クリープ力を用いずに、前記輪軸の前記左右方向の速度と、前記台車のヨーイング方向における角変位と、前記前後方向力と、前記鉄道車両の走行速度と、を用いて表現され、前記ヨーイング方向は、前記上下方向を回動軸とする回動方向であり、前記接触角算出工程は、前記第１のフィルタに含まれる前記方程式に前記更新前の接触角を与えて前記更新後の前記接触角を算出することを、所定の条件が成立するまで行うことを特徴とする。

本発明のプログラムは、前記接触角推定システムの各手段としてコンピュータを機能させることを特徴とする。

本発明によれば、鉄道車両における接触角を容易に推定することができる。

図１は、鉄道車両の概略の一例を示す図である。 図２は、鉄道車両の構成要素の主な運動の方向を概念的に示す図である。 図３は、通り狂い量の一例を示す図である。 図４は、接触角推定装置の機能的な構成の一例を示す図である。 図５は、接触角推定装置のハードウェアの構成の一例を示す図である。 図６は、接触角推定装置における事前処理の一例を示すフローチャートである。 図７は、接触角推定装置における本処理の一例を示すフローチャートである。 図８は、自己相関行列の固有値の分布の一例を示す図である。 図９は、前後方向力の測定値の時系列データ（測定値）と、前後方向力の予測値の時系列データ（計算値）の一例を示す図である。 図１０は、前後方向力の高周波成分の時系列データの一例を示す図である。 図１１は、接触角およびフランジ角の一例を説明する図である。 図１２は、フランジ接触の一例を説明する図である。 図１３は、左右の車輪に作用する力と接触角の一例を示す図である。 図１４は、接触角の時系列データ、横圧の時系列データ、輪重の測定値の時系列データ、および軌道データの一例を示す図である。 図１５は、接触角推定システムの構成の一例を示す図である。

以下、図面を参照しながら、本発明の実施形態を説明する。
（概要）
まず、本発明の実施形態の概要を説明する。
図１は、鉄道車両の概略の一例を示す図である。尚、図１において、鉄道車両は、ｘ軸の正の方向に進むものとする（ｘ軸は、鉄道車両の走行方向に沿う軸である）。また、ｚ軸は、軌道１６（地面）に対し垂直方向（鉄道車両の高さ方向）であるものとする。ｙ軸は、鉄道車両の走行方向に対して垂直な水平方向（鉄道車両の走行方向と高さ方向との双方に垂直な方向）であるものとする。また、鉄道車両は、営業車両であるものとする。尚、各図において、○の中に●が付されているものは、紙面の奥側から手前側に向かう方向を示し、○の中に×が付されているものは、紙面の手前側から奥側に向かう方向を示す。す。

図１に示すように本実施形態では、鉄道車両は、車体１１と、台車１２ａ、１２ｂと、輪軸１３ａ〜１３ｄとを有する。このように本実施形態では、１つの車体１１に、２つの台車１２ａ、１２ｂと４組の輪軸１３ａ〜１３ｄとが備わる鉄道車両を例に挙げて説明する。輪軸１３ａ〜１３ｄは、車軸１５ａ〜１５ｄとその両端に設けられた車輪１４ａ〜１４ｄとを有する。本実施形態では、台車１２ａ、１２ｂが、ボルスタレス台車である場合を例に挙げて説明する。尚、図１では、表記の都合上、輪軸１３ａ〜１３ｄの一方の車輪１４ａ〜１４ｄのみを示すが、輪軸１３ａ〜１３ｄの他方にも車輪が設けられている（図１に示す例では、車輪は合計８個ある）。また、鉄道車両は、図１に示す構成要素以外の構成要素（後述する運動方程式で説明する構成要素等）を有するが、表記の都合上、図１では、当該構成要素の図示を省略する。例えば、台車１２ａ、１２ｂは、台車枠および枕バネなどを有する。また、各輪軸１３ａ〜１３ｄのｙ軸に沿う方向の両側には、軸箱が配置される。また、台車枠と軸箱は、軸箱支持装置により相互に結合される。軸箱支持装置は、軸箱および台車枠との間に配置される装置（サスペンション）である。軸箱支持装置は、軌道１６から鉄道車両に伝わる振動を吸収する。また、軸箱支持装置は、軸箱が台車枠に対してｘ軸に沿う方向およびｙ軸に沿う方向に移動することを抑制するように（好ましくは、当該移動が生じないように）軸箱の台車枠に対する位置を規制した状態で軸箱を支持する。軸箱支持装置は、各輪軸１３ａ〜１３ｄのｙ軸に沿う方向の両側に配置される。尚、鉄道車両自体は公知の技術で実現できるので、ここでは、その詳細な説明を省略する。

鉄道車両が軌道１６上を走行すると、車輪１４ａ〜１４ｄと軌道１６との間の作用力（クリープ力）が振動源となり、輪軸１３ａ〜１３ｄ、台車１２ａ、１２ｂ、車体１１に振動が順次伝搬する。図２は、鉄道車両の構成要素（輪軸１３ａ〜１３ｄ、台車１２ａ、１２ｂ、車体１１）の主な運動の方向を概念的に示す図である。図２に示すｘ軸、ｙ軸、ｚ軸は、それぞれ、図１に示したｘ軸、ｙ軸、ｚ軸に対応する。

図２に示すように本実施形態では、輪軸１３ａ〜１３ｄ、台車１２ａ、１２ｂ、および車体１１が、ｘ軸を回動軸として回動する運動と、ｚ軸を回動軸として回動する運動と、ｙ軸に沿う方向の運動とを行う場合を例に挙げて説明する。以下の説明では、ｘ軸を回動軸として回動する運動を必要に応じてローリングと称し、ｘ軸を回動軸とする回動方向を必要に応じてローリング方向と称し、ｘ軸に沿う方向を必要に応じて前後方向と称する。尚、前後方向は、鉄道車両の走行方向である。本実施形態では、ｘ軸に沿う方向が鉄道車両の走行方向であるものとする。また、ｚ軸を回動軸として回動する運動を必要に応じてヨーイングと称し、ｚ軸を回動軸とする回動方向を必要に応じてヨーイング方向と称し、ｚ軸に沿う方向を必要に応じて上下方向と称する。尚、上下方向は、軌道１６に対し垂直な方向である。また、ｙ軸に沿う方向の運動を必要に応じて横振動と称し、ｙ軸に沿う方向を必要に応じて左右方向と称する。尚、左右方向は、前後方向（鉄道車両の走行方向）と上下方向（軌道１６に対し垂直な方向）との双方に垂直な方向である。また、左方向、右方向は、それぞれ鉄道車両の走行方向（図１に示す例ではｘ軸の正の方向）に向かって見た場合の左方向、右方向であるものとする。以下の説明において、（車輪や軌条の曲がり方向における）左、右も、それぞれ鉄道車両の走行方向に向かって見た場合の左、右であるものとする。尚、鉄道車両は、この他の運動も行うが、各実施形態では説明を簡単にするため、これらの運動については考慮しないものとする。しかしながら、これらの運動を考慮してもよい。

ここで、本発明者らは、鉄道車両の各輪軸１３ａ〜１３ｄに属する左右の車輪１４ａ〜１４ｄにおける接触角を、当該左右の車輪が受ける横圧および輪重と、当該左右の車輪における横クリープ力とに基づいて導出することができることを見出した。横クリープ力は、輪軸１３ａ〜１３ｄのヨーイング方向における回動量（角変位）と、鉄道車両の走行速度と、輪軸１３ａ〜１３ｄの左右方向の速度と、を用いて導出される。

これらのうち、輪軸１３ａ〜１３ｄのヨーイングを記述する運動方程式には、通り狂い量が含まれる。輪軸１３ａ〜１３ｄのヨーイングを記述する運動方程式を状態方程式に含めて、後述するデータ同化を行うフィルタ（カルマンフィルタ）によるフィルタリングを行う場合には、軌道１６のモデルが必要になる。また、通り狂い量は物理法則に則って記述できるものではない。従って、通り狂い量の時間微分が例えばホワイトノイズとなるように、軌道１６のモデルを作成する必要がある。そうすると、軌道１６のモデルの不確かさが、後述するデータ同化を行うフィルタ（カルマンフィルタ）によるフィルタリングの結果に影響を与える虞がある。また、状態方程式を少なくし、状態変数を減らすことにより、後述するデータ同化を行うフィルタ（カルマンフィルタ）の動作を安定させることができる。

そこで、本発明者らは、横クリープ力を、輪軸１３ａ〜１３ｄのヨーイング方向における回動量（角変位）を用いずに、輪軸１３ａ〜１３ｂ（１３ｃ〜１３ｄ）と、当該輪軸１３ａ〜１３ｂ（１３ｃ〜１３ｄ）が設けられる台車１２ａ（１２ｂ）との間に配置される部材に生じる前後方向の力と、台車１２ａ、１２ｂのヨーイング方向における回動量（角変位）とを用いて表現することを見出した。以下の説明では、この部材に生じる前後方向の力を、必要に応じて前後方向力と称する。

そうすると、各輪軸１３ａ〜１３ｄに属する左右の車輪における接触角は、当該輪軸に属する左右の車輪が受ける横圧および輪重と、当該輪軸における前後方向力および左右方向における速度と、当該輪軸が属する台車のヨーイング方向における回動量（角変位）と、鉄道車両の走行速度とを用いて１つの方程式で表現される。

鉄道車両の走行速度、横圧、輪重、および前後方向力は、測定可能な物理量である。また、輪軸１３ａ〜１３ｄの左右方向の速度および台車のヨーイング方向における回動量（角変位）も測定可能ではあるが、これらを状態変数としてデータ同化を行うフィルタ（カルマンフィルタ）によるフィルタリングを行う方が簡便である。

従って、前述した接触角を表現する方程式において未知数は、各輪軸１３ａ〜１３ｄに属する左右の車輪１４ａ〜１４ｄにおける接触角（未知数の数は２つ）になる。このように、各輪軸１３ａ〜１３ｄに属する左右の車輪１４ａ〜１４ｄにおける接触角は、１つの方程式で表現されるので、当該方程式だけでは解くことができない。そこで、本発明者らは、鉄道車両の各種の走行実績から、車輪１４ａ〜１４ｄの踏面が軌条に接触しているとき（所謂フランジ接触していないとき）には、各輪軸１３ａ〜１３ｄに属する左右の車輪１４ａ〜１４ｄにおける接触角の和は一定値とみなしてよいという知見を得た。また、本発明者らは、鉄道車両が軌道１６の曲線部として曲率が大きい領域（即ち、急カーブ）を走行しており、車輪１４ａ〜１４ｄのフランジが軌条に接触（所謂フランジ接触）しているときには、左右方向のうち、軌道１６の曲線部の曲がっている方向の車輪１４ａ〜１４ｄにおける接触角は下限値（０（ゼロ）または０（ゼロ）に近い定数値）に達するという知見を得た（例えば、鉄道車両が右回りの急カーブを走行している場合には、右側の車輪における接触角が下限値に達する）。尚、以下の説明では、軌道１６の直線部を必要に応じて直線軌道と称し、軌道１６の曲線部を必要に応じて曲線軌道と称する。

以上のようにすれば、フランジ接触していない場合には、２つの未知数（各輪軸１３ａ〜１３ｄに属する左右の車輪１４ａ〜１４ｄにおける接触角）に対し２つの方程式を構成することができる。一方、フランジ接触している場合には、未知数は、各輪軸１３ａ〜１３ｄに属する左右の車輪１４ａ〜１４ｄにおける接触角のうちの何れか１つになる。従って、各輪軸１３ａ〜１３ｄに属する左右の車輪１４ａ〜１４ｄにおける接触角を推定することができる。

ところで、後述するデータ同化を行うフィルタ（カルマンフィルタ）によるフィルタリングを行う場合に、曲線軌道を走行する鉄道車両の運動を記述する運動方程式を用いて状態方程式を構成すると、状態変数が発散する虞がある。そこで、データ同化を行うフィルタ（カルマンフィルタ）によるフィルタリングを行う場合の状態方程式を、直線軌道を走行する鉄道車両の運動を記述する運動方程式を用いて構成する。
曲線軌道を走行する鉄道車両の運動を記述する運動方程式では、走行時に鉄道車両が受ける遠心力等を考慮する必要がある。従って、曲線軌道を走行する鉄道車両の運動を記述する運動方程式には、軌条の曲率半径を含む項が含まれる。よって、鉄道車両が曲線軌道を走行しているときに、直線軌道を走行する鉄道車両の運動を記述する運動方程式を用いて構成したデータ同化を行うフィルタ（カルマンフィルタ）を用いて状態変数を導出すると、状態変数を高精度に導出することができなくなる虞がある。

本発明者らは、鉄道車両が曲線軌道を走行する場合には、直線軌道を走行しているときに対し、前後方向力の測定値が或るバイアスを持つことに着目した。そこで、本発明者らは、前後方向力の測定値の時系列データから低周波成分（前記バイアスの挙動）を低減させることにより、鉄道車両の直線軌道の走行時における運動を記述する運動方程式に基づく式を用いて後述するデータ同化を行うフィルタ（カルマンフィルタ）を構成しても、状態変数の推定値から、鉄道車両が曲線軌道を走行することに起因する低周波成分を低減することができると考えた。

（運動方程式）
次に、鉄道車両の直線軌道の走行時における運動を記述する運動方程式の一例を説明する。本実施形態では、特許文献２に記載の運動方程式を例に挙げ、鉄道車両が２１自由度を有する場合を例に挙げて説明する。即ち、輪軸１３ａ〜１３ｄが、左右方向における運動（横振動）とヨーイング方向における運動（ヨーイング）とを行うものとする（２×４セット＝８自由度）。また、台車１２ａ、１２ｂが、左右方向における運動（横振動）とヨーイング方向における運動（ヨーイング）とローリング方向における運動（ローリング）とを行うものとする（３×２セット＝６自由度）。また、車体１１が、左右方向における運動（横振動）とヨーイング方向における運動（ヨーイング）とローリング方向における運動（ローリング）とを行うものとする（３×１セット＝３自由度）。また、台車１２ａ、１２ｂに対してそれぞれ設けられている空気バネ（枕バネ）が、ローリング方向における運動（ローリング）を行うものとする（１×２セット＝２自由度）。また、台車１２ａ、１２ｂに対してそれぞれ設けられているヨーダンパが、ヨーイング方向における運動（ヨーイング）を行うものとする（１×２セット＝２自由度）。

尚、自由度は、２１自由度に限定されない。自由度を大きくすれば計算精度は向上するが、計算負荷が高くなる。また、後述するカルマンフィルタの動作が安定しなくなる虞がある。これらの点を考慮して自由度を適宜決定することができる。また、以下の運動方程式は、それぞれの構成要素（車体１１、台車１２ａ、１２ｂ、輪軸１３ａ〜１３ｄ）のそれぞれの方向（左右方向、ヨーイング方向、ローリング方向）の動作を、例えば、非特許文献１、２の記載に基づいて表すことにより実現することができる。従って、ここでは、それぞれの運動方程式の概要を説明し、詳細な説明を省略する。尚、以下の各式には、軌道１６（軌条）の曲率半径（曲率）を含む項が存在しない。即ち、以下の各式は、鉄道車両が直線軌道を走行することを表現する式となる。鉄道車両が曲線軌道を走行することを表現する式において、軌道１６（軌条）の曲率半径を無限大（曲率を０（ゼロ））とすることにより、鉄道車両が直線軌道を走行することを表現する式が得られる。

以下の各式において、添え字ｗは、輪軸１３ａ〜１３ｄを表す。添え字ｗ（のみ）が付されている変数は、輪軸１３ａ〜１３ｄで共通であることを表す。添え字ｗ１、ｗ２、ｗ３、ｗ４はそれぞれ、輪軸１３ａ、１３ｂ、１３ｃ、１３ｄを表す。
添え字ｔ、Ｔは、台車１２ａ、１２ｂを表す。添え字ｔ、Ｔ（のみ）が付されている変数は、台車１２ａ、１２ｂで共通であることを表す。添え字ｔ１、ｔ２はそれぞれ、台車１２ａ、１２ｂを表す。
添え字ｂ、Ｂは、車体１１であることを表す。

添え字ｘは、前後方向またはローリング方向を表し、添え字ｙは、左右方向を表し、添え字ｚは、上下方向またはヨーイング方向を表す。
また、変数の上に付している「・・」、「・」はそれぞれ、２階時間微分、１階時間微分を表す。
尚、以下の運動方程式の説明に際し、必要に応じて、既出の変数の説明を省略する。また、運動方程式自体は、特許文献２に記載されているものと同じである。

［輪軸の横振動］
輪軸１３ａ〜１３ｄの横振動（左右方向における運動）を記述する運動方程式は、以下の（１）式〜（４）式で表される。

ｍ_wは、輪軸１３ａ〜１３ｄの質量である。ｙ_w1・・（式において・・はｙ_w1の上に付される（以下、その他の変数についても同様））は、輪軸１３ａの左右方向における加速度である。ｆ₂は、横クリープ係数である（尚、横クリープ係数ｆ₂は、輪軸１３ａ〜１３ｄ毎に与えられてもよい）。ｖは、鉄道車両の走行速度である。ｙ_w1・（式において・はｙ_w1の上に付される（以下、その他の変数についても同様））は、輪軸１３ａの左右方向における速度である。Ｃ_wyは、軸箱と輪軸とをつなぐ軸箱支持装置の左右方向におけるダンピング定数である。ｙ_t1・は、台車１２ａの左右方向における速度である。ａは、台車１２ａ、１２ｂに設けられている輪軸１３ａ・１３ｂ、１３ｃ・１３ｄ間の前後方向における距離の１／２を表す（台車１２ａ、１２ｂに設けられている輪軸１３ａ・１３ｂ、１３ｃ・１３ｄ間の距離は２ａになる）。ψ_t1・は、台車１２ａのヨーイング方向における角速度である。ｈ₁は、車軸の中心と台車１２ａの重心との上下方向における距離である。φ_t1・は、台車１２ａのローリング方向における角速度である。ψ_w1は、輪軸１３ａのヨーイング方向における回動量（角変位）である。Ｋ_wyは、軸箱支持装置の左右方向のバネ定数である。ｙ_w1は、輪軸１３ａの左右方向における変位である。ｙ_t1は、台車１２ａの左右方向における変位である。ψ_t1は、台車１２ａのヨーイング方向における回動量（角変位）である。φ_t1は、台車１２ａのローリング方向における回動量（角変位）である。尚、（２）式〜（４）式の各変数は、前述した添え字の意味に従って（１）式の変数を読み替えることにより表される。

［輪軸のヨーイング］
輪軸１３ａ〜１３ｄのヨーイングを記述する運動方程式は、以下の（５）式〜（８）式で表される。

Ｉ_wzは、輪軸１３ａ〜１３ｄのヨーイング方向における慣性モーメントである。ψ_w1・・は、輪軸１３ａのヨーイング方向における角加速度である。ｆ₁は、縦クリープ係数である。ｂは、輪軸１３ａ〜１３ｄに取り付けられている２つの車輪と軌道１６（軌条）との接点の間の左右方向における距離である。ψ_w1・は、輪軸１３ａのヨーイング方向における角速度である。Ｃ_wxは、軸箱支持装置の前後方向のダンピング定数である。ｂ₁は、軸箱支持装置の左右方向における間隔の１／２の長さを表す（１つの輪軸に対して左右に設けられている２つの軸箱支持装置の左右方向における間隔は２ｂ₁になる）。γは、踏面勾配である。ｒは、車輪１４ａ〜１４ｄの半径である。ｙ_R1は、輪軸１３ａの位置での通り狂い量である。ｓ_aは、車軸１５ａ〜１５ｄの中心から軸箱支持バネまでの前後方向におけるオフセット量である。ｙ_t1は、台車１２ａの左右方向における変位である。Ｋ_wxは、軸箱支持装置の前後方向のバネ定数である。尚、（６）式〜（８）式の各変数は、前述した添え字の意味に従って（５）式の変数を読み替えることにより表される。ただし、ｙ_R2、ｙ_R3、ｙ_R4はそれぞれ、輪軸１３ｂ、１３ｃ、１３ｄの位置での通り狂い量である。

ここで、通り狂いとは、日本工業規格（ＪＩＳ Ｅ １００１：２００１）に記載されているように、軌条の長手方向の左右の変位である。通り狂い量は、その変位の量である。図３に、輪軸１３ａの位置での通り狂い量ｙ_R1の一例を示す。図３（ａ）では、軌道１６が、直線軌道である場合を例に挙げて説明する。図３（ｂ）では、軌道１６が、曲線軌道である場合を例に挙げて説明する。図３（ａ）および図３（ｂ）において、１６ａは、軌条を示し、１６ｂは、枕木を示す。図３（ａ）では、輪軸１３ａの車輪１４ａが位置３０１で軌条１６ａと接触しているとする。図３（ｂ）では、輪軸１３ａの車輪１４ａが位置３０２で軌条１６ａと接触しているとする。輪軸１３ａの位置での通り狂い量ｙ_R1は、輪軸１３ａの車輪１４ａと軌条１６ａとの接触位置と、正規の状態であると仮定した場合の軌条１６ａの位置と、の左右方向の距離である。輪軸１３ａの位置とは、輪軸１３ａの車輪１４ａと軌条１６ａとの接触位置である。輪軸１３ｂ、１３ｃ、１３ｄの位置での通り狂い量ｙ_R2、ｙ_R3、ｙ_R4も、輪軸１３ａの位置での通り狂い量ｙ_R1と同様に定義される。

［台車の横振動］
台車１２ａ、１２ｂの横振動（左右方向における運動）を記述する運動方程式は、以下の（９）式、（１０）式で表される。

ｍ_Tは、台車１２ａ、１２ｂの質量である。ｙ_t1・・は、台車１２ａの左右方向における加速度である。ｃ´₂は、左右動ダンパのダンピング定数である。ｈ₄は、台車１２ａの重心と左右動ダンパとの上下方向における距離である。ｙ_b・は、車体１１の左右方向における速度である。Ｌは、台車１２ａ、１２ｂの中心間の前後方向における間隔の１／２を表す（台車１２ａ、１２ｂの中心間の前後方向における間隔は２Ｌになる）。ψ_b・は、車体１１のヨーイング方向における角速度である。ｈ₅は、左右動ダンパと車体１１の重心との間の上下方向における距離である。φ_b・は、車体１１のローリング方向における角速度である。ｙ_w2・は、輪軸１３ｂの左右方向における速度である。ｋ´₂は、空気バネ（枕バネ）の左右方向のバネ定数である。ｈ₂は、台車１２ａ、１２ｂの重心と空気バネ（枕バネ）の中心との間の上下方向における距離である。ｙ_bは、車体１１の左右方向における変位である。ψ_bは、車体１１のヨーイング方向における回動量（角変位）である。ｈ₃は、空気バネ（枕バネ）の中心と車体１１の重心との間の上下方向における距離である。φ_bは、車体１１のローリング方向における回動量（角変位）である。尚、（１０）式の各変数は、前述した添え字の意味に従って（９）式の変数を読み替えることにより表される。

［台車のヨーイング］
台車１２ａ、１２ｂのヨーイングを記述する運動方程式は、以下の（１１）式、（１２）式で表される。

Ｉ_Tzは、台車１２ａ、１２ｂのヨーイング方向における慣性モーメントである。ψ_t1・・は、台車１２ａのヨーイング方向における角加速度である。ψ_w2・は、輪軸１３ｂのヨーイング方向における角速度である。ψ_w2は、輪軸１３ｂのヨーイング方向における回動量（角変位）である。ｙ_w2は、輪軸１３ｂの左右方向における変位である。ｋ´₀は、ヨーダンパのゴムブッシュ剛性である。ｂ´₀は、台車１２ａ、１２ｂに対し左右に配置された２つのヨーダンパの左右方向における間隔の１／２を表す（台車１２ａ、１２ｂに対し左右に配置された２つのヨーダンパの左右方向における間隔は２ｂ´₀になる）。ψ_y1は、台車１２ａに配置されたヨーダンパのヨーイング方向における回動量（角変位）である。ｋ´´₂は、空気バネ（枕バネ）の左右方向のバネ定数である。ｂ₂は、台車１２ａ、１２ｂに対し左右に配置された２つの空気バネ（枕バネ）の左右方向における間隔の１／２を表す（台車１２ａ、１２ｂに対し左右に配置された２つの空気バネ（枕バネ）の左右方向における間隔は２ｂ₂になる）。尚、（１２）式の各変数は、前述した添え字の意味に従って（１１）式の変数を読み替えることにより表される。

［台車のローリング］
台車１２ａ、１２ｂのローリングを記述する運動方程式は、以下の（１３）式、（１４）式で表される。

Ｉ_Txは、台車１２ａ、１２ｂのローリング方向における慣性モーメントである。φ_t1・・は、台車１２ａのローリング方向における角加速度である。ｃ₁は、軸ダンパの上下方向のダンピング定数である。ｂ´₁は、台車１２ａ、１２ｂに対し左右に配置された２つの軸ダンパの左右方向における間隔の１／２を表す（台車１２ａ、１２ｂに対し左右に配置された２つの軸ダンパの左右方向における間隔は２ｂ´₁になる）。ｃ₂は、空気バネ（枕バネ）の上下方向のダンピング定数である。φ_a1・は、台車１２ａに配置された空気バネ（枕バネ）のローリング方向における角速度である。ｋ₁は、軸バネの上下方向のバネ定数である。λは、空気バネ（枕バネ）の本体の容積を補助空気室の容積で割った値である。ｋ₂は、空気バネ（枕バネ）の上下方向のバネ定数である。φ_a1は、台車１２ａに配置された空気バネ（枕バネ）のローリング方向における回動量（角変位）である。ｋ₃は、空気バネ（枕バネ）の有効受圧面積の変化による等価剛性である。尚、（１４）式の各変数は、前述した添え字の意味に従って（１３）式の変数を読み替えることにより表される。ただし、φ_a2は、台車１２ｂに配置された空気バネ（枕バネ）のローリング方向における回動量（角変位）である。

［車体の横振動］
車体１１の横振動（左右方向における運動）を記述する運動方程式は、以下の（１５）式で表される。

ｍ_Bは、台車１２ａ、１２ｂの質量である。ｙ_b・・は、車体１１の左右方向における加速度である。ｙ_t2・は、台車１２ｂの左右方向における速度である。φ_t2・は、台車１２ｂのローリング方向における角速度である。ｙ_t2は、台車１２ｂの左右方向における変位である。φ_t2は、台車１２ｂのローリング方向における回動量（角変位）である。

［車体のヨーイング］
車体１１のヨーイングを記述する運動方程式は、以下の（１６）式で表される。

Ｉ_Bzは、車体１１のヨーイング方向における慣性モーメントである。ψ_b・・は、車体１１のヨーイング方向における角加速度である。ｃ₀は、ヨーダンパの前後方向のダンピング定数である。ψ_y1・は、台車１２ａに配置されたヨーダンパのヨーイング方向における角速度である。ψ_y2・は、台車１２ｂに配置されたヨーダンパのヨーイング方向における角速度である。ψ_t2は、台車１２ｂのヨーイング方向における回動量（角変位）である。

［車体のローリング］
車体１１のローリングを記述する運動方程式は、以下の（１７）式で表される。

Ｉ_Bxは、車体１１のローリング方向における慣性モーメントである。φ_b・・は、車体１１のローリング方向における角加速度である。

［ヨーダンパのヨーイング］
台車１２ａに配置されたヨーダンパ、台車１２ｂに配置されたヨーダンパのヨーイングを記述する運動方程式は、それぞれ以下の（１８）式、（１９）式で表される。

ψ_y2は、台車１２ｂに配置されたヨーダンパのヨーイング方向における回動量（角変位）である。

［空気バネ（枕バネ）のローリング］
台車１２ａに配置された空気バネ（枕バネ）、台車１２ｂに配置された空気バネ（枕バネ）のローリングを記述する運動方程式は、それぞれ以下の（２０）式、（２１）式で表される。

φ_a2・は、台車１２ｂに配置された空気バネ（枕バネ）のローリング方向における角速度である。

（前後方向力）
次に、前後方向力について説明する。尚、前後方向力自体は、特許文献２に記載されているものと同じである。
１つの輪軸における左右の車輪のうち一方の車輪における縦クリープ力と他方の車輪における縦クリープ力との同相の成分は、ブレーキ力や駆動力に対応する成分である。従って、縦クリープ力の逆相成分に対応するように前後方向力を定めるのが好ましい。縦クリープ力の逆相成分とは、１つの輪軸における左右の車輪のうち一方の車輪における縦クリープ力と他方の車輪における縦クリープ力との相互に逆位相となる成分である。即ち、縦クリープ力の逆相成分とは、縦クリープ力の、車軸をねじる方向の成分である。この場合、前後方向力は、１つの輪軸の左右方向の両側に取り付けられた２つの前記部材に生じる力の前後方向の成分のうち、相互に逆位相となる成分となる。

以下に、縦クリープ力の逆相成分に対応するように前後方向力を定める場合の前後方向力の具体例について説明する。
軸箱支持装置が、モノリンク式の軸箱支持装置である場合、軸箱支持装置は、リンクを備えており、軸箱と台車枠とがリンクにより連結されている。このリンクの両端にはゴムブッシュが取り付けられる。この場合、前後方向力は、１つの輪軸の左右方向の端にそれぞれ１つずつ取り付けられる２つのリンクのそれぞれが受ける荷重の前後方向の成分のうち、相互に逆位相となる成分になる。また、リンクの配置および構成により、リンクは、前後方向、左右方向、前後方向の荷重のうち主に前後方向の荷重を受ける。従って、例えば、各リンクに歪ゲージを１つ取り付ければよい。この歪ゲージの測定値を用いて、当該リンクが受ける荷重の前後方向の成分を導出することにより、前後方向力の測定値を得る。また、このようにすることに替えて、リンクに取り付けられたゴムブッシュの前後方向の変位を変位計で測定してもよい。この場合、測定した変位と当該ゴムブッシュのバネ定数との積を、前後方向力の測定値とする。軸箱支持装置が、モノリンク式の軸箱支持装置である場合、前述した、軸箱を支持するための部材は、リンクまたはゴムブッシュになる。

尚、リンクに取り付けられる歪ゲージにより測定される荷重には、前後方向の成分だけでなく、左右方向の成分および上下方向の成分のうち少なくとも何れか一方の成分が含まれる場合がある。しかしながら、このような場合であっても、軸箱支持装置の構造上、リンクが受ける左右方向の成分の荷重および上下方向の成分の荷重は、前後方向の成分の荷重に比べて十分に小さい。従って、各リンクに歪ゲージを１つ取り付けるだけで、実用上要求される精度を有する前後方向力の測定値を得ることができる。このように、前後方向力の測定値には、前後方向の成分以外の成分が含まれることがある。従って、上下方向および左右方向の歪みがキャンセルされるように３つ以上の歪ゲージを各リンクに取り付けてもよい。このようにすれば、前後方向力の測定値の精度を向上させることができる。

軸箱支持装置が、軸はり式の軸箱支持装置である場合、軸箱支持装置は、軸はりを備えており、軸箱と台車枠とが、軸はりにより連結されている。軸はりは、軸箱と一体に構成されていてもよい。この軸はりの台車枠側の端にはゴムブッシュが取り付けられる。この場合、前後方向力は、１つの輪軸の左右方向の端にそれぞれ１つずつ取り付けられる２つの軸はりのそれぞれが受ける荷重の前後方向の成分のうち、相互に逆位相となる成分になる。また、軸はりの配置構成により、軸はりは、前後方向、左右方向、上下方向の荷重のうち前後方向の荷重に加えて、左右方向の荷重も受けやすい。従って、例えば、左右方向の歪みがキャンセルされるように２つ以上の歪ゲージを各軸はりに取り付ける。これらの歪ゲージの測定値を用いて、軸はりが受ける荷重の前後方向の成分を導出することにより、前後方向力の測定値を得る。また、このようにすることに替えて、軸はりに取り付けられたゴムブッシュの前後方向の変位を変位計で測定してもよい。この場合、測定した変位と当該ゴムブッシュのバネ定数との積を、前後方向力の測定値とする。軸箱支持装置が、軸はり式の軸箱支持装置である場合、前述した、軸箱を支持するための部材は、軸はりまたはゴムブッシュになる。

尚、軸はりに取り付けられる歪ゲージにより測定される荷重には、前後方向および左右方向の成分だけでなく、上下方向の成分が含まれる場合がある。しかしながら、このような場合であっても、軸箱支持装置の構造上、軸はりが受ける上下方向の成分の荷重は、前後方向の成分の荷重および左右方向の成分の荷重に比べて十分に小さい。従って、軸はりが受ける上下方向の成分の荷重をキャンセルするように歪ゲージを取り付けなくても、実用上要求される精度を有する前後方向力の測定値を得ることができる。このように、計測された前後方向力には、前後方向の成分以外の成分が含まれることがあり、左右方向の歪みに加えて上下方向の歪みもキャンセルされるように３つ以上の歪ゲージを各軸はりに取り付けてもよい。このようにすれば、前後方向力の測定値の精度を向上させることができる。

軸箱支持装置が、板バネ式の軸箱支持装置である場合、軸箱支持装置は、板バネを備えており、軸箱と台車枠とが、板バネにより連結されている。この板バネの端にはゴムブッシュが取り付けられる。この場合、前後方向力は、１つの輪軸の左右方向の端にそれぞれ１つずつ取り付けられる２つの板バネのそれぞれが受ける荷重の前後方向の成分のうち、相互に逆位相となる成分になる。また、板バネの配置構成により、板バネは、前後方向、左右方向、上下方向の荷重のうち前後方向の荷重に加えて、左右方向の荷重および上下方向の荷重も受けやすい。従って、例えば、左右方向および上下方向の歪みがキャンセルされるように３つ以上の歪ゲージを各板バネに取り付ける。これらの歪ゲージの測定値を用いて、板バネが受ける荷重の前後方向の成分を導出することにより、前後方向力の測定値を得る。また、このようにすることに替えて、板バネに取り付けられたゴムブッシュの前後方向の変位を変位計で測定してもよい。この場合、測定した変位と当該ゴムブッシュのバネ定数との積を、前後方向力の測定値とする。軸箱支持装置が、板バネ式の軸箱支持装置である場合、前述した、軸箱を支持するための部材は、板バネまたはゴムブッシュになる。

尚、前述した変位計としては、公知のレーザ変位計や渦電流式の変位計を用いることができる。
また、ここでは、軸箱支持装置の方式が、モノリンク式、軸はり式、および板バネ式である場合を例に挙げて、前後方向力を説明した。しかしながら、軸箱支持装置の方式は、モノリンク式、軸はり式、および板バネ式に限定されない。軸箱支持装置の方式に合わせて、モノリンク式、軸はり式、および板バネ式と同様に、前後方向力を定めることができる。
また、以下では、説明を簡単にするために、１つの輪軸について１つの前後方向力の測定値が得られる場合を例に挙げて説明する。即ち、図１に示す鉄道車両は、４つの輪軸１３ａ〜１３ｄを有する。従って、４つの前後方向力Ｔ₁〜Ｔ₄の測定値が得られる。

（第１の実施形態）
次に、本発明の第１の実施形態を説明する。

＜接触角推定装置４００＞
図４は、接触角推定装置４００の機能的な構成の一例を示す図である。図５は、接触角推定装置４００のハードウェアの構成の一例を示す図である。図６は、接触角推定装置４００における事前処理の一例を示すフローチャートである。図７は、接触角推定装置４００における本処理の一例を示すフローチャートである。本実施形態では、図１に示すように、接触角推定装置４００が、鉄道車両に搭載される場合を例に挙げて示す。

図４において、接触角推定装置４００は、その機能として、状態方程式記憶部４０１、観測方程式記憶部４０２、データ取得部４０３、第１の周波数調整部４０４、フィルタ演算部４０５、第２の周波数調整部４０６、接触角算出部４０７、および出力部４０８を有する。

図５において、接触角推定装置４００は、ＣＰＵ５０１、主記憶装置５０２、補助記憶装置５０３、通信回路５０４、信号処理回路５０５、画像処理回路５０６、Ｉ／Ｆ回路５０７、ユーザインターフェース５０８、ディスプレイ５０９、およびバス５１０を有する。

ＣＰＵ５０１は、接触角推定装置４００の全体を統括制御する。ＣＰＵ５０１は、主記憶装置５０２をワークエリアとして用いて、補助記憶装置５０３に記憶されているプログラムを実行する。主記憶装置５０２は、データを一時的に格納する。補助記憶装置５０３は、ＣＰＵ５０１によって実行されるプログラムの他、各種のデータを記憶する。補助記憶装置５０３は、後述する状態方程式および観測方程式を記憶する。状態方程式記憶部４０１および観測方程式記憶部４０２は、例えば、ＣＰＵ５０１および補助記憶装置５０３を用いることにより実現される。

通信回路５０４は、接触角推定装置４００の外部との通信を行うための回路である。通信回路５０４は、例えば、前後方向力の測定値や、車体１１、台車１２ａ、１２ｂ、および輪軸１３ａ〜１３ｄの左右方向における加速度の測定値の情報を受信する。通信回路５０４は、接触角推定装置４００の外部と無線通信を行っても有線通信を行ってもよい。通信回路５０４は、無線通信を行う場合、鉄道車両に設けられるアンテナに接続される。

信号処理回路５０５は、通信回路５０４で受信された信号や、ＣＰＵ５０１による制御に従って入力した信号に対し、各種の信号処理を行う。データ取得部４０３は、例えば、ＣＰＵ５０１、通信回路５０４、および信号処理回路５０５を用いることにより実現される。また、第１の周波数調整部４０４、フィルタ演算部４０５、第２の周波数調整部４０６および接触角算出部４０７は、例えば、ＣＰＵ５０１および信号処理回路５０５を用いることにより実現される。

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

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

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

［状態方程式記憶部４０１、Ｓ６０１］
状態方程式記憶部４０１は、状態方程式を記憶する。本実施形態では、特許文献２に記載の状態方程式を用いる場合を例に挙げて説明する。前述したように、本実施形態では、（５）式〜（８）式の輪軸１３ａ〜１３ｄのヨーイングを記述する運動方程式を状態方程式に含めずに、以下のようにして状態方程式を構成する。
まず、（９）式、（１０）式の台車１２ａ、１２ｂの横振動（左右方向における運動）を記述する運動方程式と、（１３）式、（１４）式の台車１２ａ、１２ｂのローリングを記述する運動方程式と、（１５）式の車体１１の横振動（左右方向における運動）を記述する運動方程式と、（１６）式の車体１１のヨーイングを記述する運動方程式と、（１７）式の車体１１のローリングを記述する運動方程式と、（１８）式、（１９）式の台車１２ａに配置されたヨーダンパ、台車１２ｂに配置されたヨーダンパのヨーイングを記述する運動方程式と、（２０）式、（２１）式の台車１２ａに配置された空気バネ（枕バネ）、台車１２ｂに配置された空気バネ（枕バネ）のローリングを記述する運動方程式については、これらをそのまま用いて状態方程式を構成する。

一方、（１）式〜（４）式の輪軸１３ａ〜１３ｄの横振動（左右方向における運動）を記述する運動方程式と、（１１）式、（１２）式の台車１２ａ、１２ｂのヨーイングを記述する運動方程式には、輪軸１３ａ〜１３ｄのヨーイング方向における回動量（角変位）ψ_w1〜ψ_w4や角速度ψ_w1・〜ψ_w4・が含まれる。（１）式〜（４）式および（１１）式、（１２）式からこれらの変数を消去したものを用いて状態方程式を構成する。

まず、輪軸１３ａ〜１３ｄにおける前後方向力Ｔ₁〜Ｔ₄は、以下の（２２）式〜（２５）式で表される。このように、前後方向力Ｔ₁〜Ｔ₄は、輪軸のヨーイング方向の角変位ψ_w1〜ψ_w4と、当該輪軸が設けられる台車のヨーイング方向の角変位ψ_t1〜ψ_t2との差に応じて定まる。

以下の（２６）式〜（２９）式のように、変換変数ｅ₁〜ｅ₄を定義する。このように、変換変数ｅ₁〜ｅ₄は、台車のヨーイング方向の角変位ψ_t1〜ψ_t2と輪軸のヨーイング方向の角変位ψ_w1〜ψ_w4との差で定義される。変換変数ｅ₁〜ｅ₄は、台車のヨーイング方向の角変位ψ_t1〜ψ_t2と輪軸のヨーイング方向の角変位ψ_w1〜ψ_w4とを相互に変換するための変数である。

（２６）式〜（２９）式を式変形すると、以下の（３０）式〜（３３）式が得られる。

（３０）式〜（３３）式を、（１）式〜（４）式の輪軸１３ａ〜１３ｄの横振動（左右方向における運動）を記述する運動方程式に代入すると、以下の（３４）式〜（３７）式が得られる。

このように、（１）式〜（４）式の輪軸１３ａ〜１３ｄの横振動（左右方向における運動）を記述する運動方程式を、変換変数ｅ₁〜ｅ₄を用いて表現することで、当該運動方程式に含まれていた輪軸１３ａ〜１３ｄのヨーイング方向における回動量（角変位）ψ_w1〜ψ_w4を消去することができる。

（２２）式〜（２５）式を、（１１）式、（１２）式の台車１２ａ、１２ｂのヨーイングを記述する運動方程式に代入すると、以下の（３８）式、（３９）式が得られる。

このように、（１１）式、（１２）式の台車１２ａ、１２ｂのヨーイングを記述する運動方程式を、前後方向力Ｔ₁〜Ｔ₄を用いて表現することで、当該運動方程式に含まれていた輪軸１３ａ〜１３ｄのヨーイング方向における角変位ψ_w1〜ψ_w4および角速度ψ_w1・〜ψ_w4・を消去できる。

また、（２６）式〜（２９）式を、（２２）式〜（２５）式に代入すると、以下の（４０）式〜（４３）式が得られる。

以上のように本実施形態では、（３４）式〜（３７）式のようにして輪軸１３ａ〜１３ｄの横振動（左右方向における運動）を記述する運動方程式を表すと共に、（３８）式、（３９）式のようにして台車１２ａ、１２ｂのヨーイングを記述する運動方程式を表し、これらを用いて状態方程式を構成する。また、（４０）式〜（４３）式は、常微分方程式であり、その解である変換変数ｅ₁〜ｅ₄の実績値は、輪軸１３ａ〜１３ｄにおける前後方向力Ｔ₁〜Ｔ₄の値を用いることにより求めることができる。ここで、前後方向力Ｔ₁〜Ｔ₄の値は、後述する第１の周波数調整部４０４により、前後方向力の測定値の時系列データから鉄道車両が軌道の曲線部を走行することに起因して生じる低周波成分の信号強度が低減されたものである。

このようにして求めた変換変数ｅ₁〜ｅ₄の実績値を、（３４）式〜（３７）式に与える。また、輪軸１３ａ〜１３ｄにおける前後方向力Ｔ₁〜Ｔ₄の値を（３８）式、（３９）式に与える。ここで、前後方向力Ｔ₁〜Ｔ₄の値は、後述する第１の周波数調整部４０４により、前後方向力の測定値の時系列データから鉄道車両が軌道の曲線部を走行することに起因して生じる低周波成分の信号強度が低減されたものである。

本実施形態では、以下の（４４）式に示す変数を状態変数とし、（９）式、（１０）式、（１３）式〜（２１）式、（３４）式〜（３９）式の運動方程式を用いて状態方程式を構成する。

状態方程式記憶部４０１は、例えば、以上のようにして構成される状態方程式を、オペレータによるユーザインターフェース５０８の操作に基づいて入力し、記憶する。

［観測方程式記憶部４０２、Ｓ６０２］
観測方程式記憶部４０２は、観測方程式を記憶する。本実施形態では、車体１１の左右方向における加速度、台車１２ａ、１２ｂの左右方向における加速度、および輪軸１３ａ〜１３ｄの左右方向における加速度を観測変数とする。この観測変数は、後述するカルマンフィルタによるフィルタリングの観測変数である。本実施形態では、（３４）式〜（３７）式、（９）式、（１０）式、および（１５）式の横振動を記述する運動方程式を用いて観測方程式を構成する。観測方程式記憶部４０２は、例えば、このようにして構成される観測方程式を、オペレータによるユーザインターフェース５０８の操作に基づいて入力し、記憶する。

以上のようにして、状態方程式および観測方程式が接触角推定装置４００に記憶された後に、データ取得部４０３、第１の周波数調整部４０４、フィルタ演算部４０５、第２の周波数調整部４０６、接触角算出部４０７、および出力部４０８が起動する。即ち、図６のフローチャートによる事前処理が終了した後に、図７のフローチャートによる本処理が開始する。

［データ取得部４０３、Ｓ７０１］
データ取得部４０３は、計測データを所定のサンプリング周期で取得する。
本実施形態では、データ取得部４０３は、計測データとして、車体１１の左右方向における加速度の測定値の時系列データ、台車１２ａ、１２ｂの左右方向における加速度の測定値の時系列データ、輪軸１３ａ〜１３ｄの左右方向における加速度の測定値の時系列データ、鉄道車両の走行速度ｖの時系列データ、車輪１４ａ〜１４ｄにおける輪重の測定値の時系列データ、および車輪１４ａ〜１４ｄにおける横圧の測定値の時系列データを取得する。
各加速度は、例えば、車体１１、台車１２ａ、１２ｂ、および輪軸１３ａ〜１３ｄにそれぞれ取り付けられた歪ゲージと、当該歪ゲージの測定値を用いて加速度を演算する演算装置とを用いることにより測定される。また、輪重および横圧は、例えば、車輪１４ａ〜１４ｄに取り付けられた歪ゲージと、当該歪ゲージの測定値を用いて輪重および横圧を演算する演算装置とを用いることにより測定される。尚、加速度、輪重、および横圧の測定は、公知の技術で実現することができるので、その詳細な説明を省略する。

また、データ取得部４０３は、計測データとして、前後方向力の測定値の時系列データを取得する。前後方向力の測定の方法は、前述した通りである。
データ取得部４０３は、例えば、前述した演算装置との通信を行うことにより、計測データを取得することができる。

［第１の周波数調整部４０４、Ｓ７０２］
第１の周波数調整部４０４は、データ取得部４０３により取得された計測データのうち、前後方向力の測定値の時系列データに含まれる低周波成分の信号強度を低減（好ましくは除去）する。この低周波成分の信号は、鉄道車両が直線軌道を走行している場合には計測されないが、鉄道車両が曲線軌道を走行している場合に計測される信号である。即ち、鉄道車両が曲線軌道を走行している場合に計測される信号は、鉄道車両が直線軌道を走行している場合に計測される信号に、この低周波成分の信号が重畳された信号と見なすことができる。

本発明者らは、自己回帰モデル（ＡＲ（Auto-regressive）モデル）を修正したモデルを考案した。そして、本発明者らは、このモデルを用いて、前後方向力の測定値の時系列データに含まれる低周波成分の信号強度を低減することに想到した。以下の説明では、本発明者らが考案したモデルを、修正自己回帰モデルと称する。これに対し、公知の自己回帰モデルを、単に自己回帰モデルと称する。以下、修正自己回帰モデルの一例について説明する。

時刻ｋ（１≦ｋ≦Ｍ）における物理量の時系列データｙの値をｙ_kとする。Ｍは、物理量の時系列データｙがどの時刻までのデータを含むかを示す数であり、予め設定されている。以下の説明では、物理量の時系列データを必要に応じてデータｙと略称する。データｙの値ｙ_kを近似する自己回帰モデルは、例えば、以下の（４５）式のようになる。（４５）式に示すように、自己回帰モデルとは、データｙにおける時刻ｋ（ｍ＋１≦ｋ≦Ｍ）の物理量の予測値ｙ＾_kを、データｙにおけるその時刻ｋよりも前の時刻ｋ−ｌ（１≦ｌ≦ｍ）の物理量の実績値ｙ_k-lを用いて表す式である。尚、ｙ＾_kは、（４５）式において、ｙ_kの上に＾を付けて表記したものである。

（４５）式におけるαは、自己回帰モデルの係数である。ｍは、自己回帰モデルにおいて時刻ｋにおけるデータｙの値ｙ_kを近似するために用いられるデータｙの値の数であって、その時刻ｋよりも前の連続する時刻ｋ−１〜ｋ−ｍにおけるデータｙの値ｙ_k-1〜ｙ_k-mの数である。ｍは、Ｍ未満の整数である。ｍとして、例えば、１５００を用いることができる。

続いて、最小二乗法を用いて、自己回帰モデルによる時刻ｋにおける物理量の予測値ｙ＾_kが、値ｙ_kに近似するための条件式を求める。自己回帰モデルによる時刻ｋにおける物理量の予測値ｙ＾_kが値ｙ_kに近似するための条件として、例えば、自己回帰モデルによる時刻ｋにおける物理量の予測値ｙ＾_kと値ｙ_kとの二乗誤差を最小化するとする条件を採用することができる。即ち、自己回帰モデルによる時刻ｋにおける物理量の予測値ｙ＾_kを値ｙ_kに近似するために最小二乗法を用いる。以下の（４６）式は、自己回帰モデルによる時刻ｋにおける物理量の予測値ｙ＾_kを値ｙ_kとの二乗誤差を最小にするための条件式である。

（４６）式より、以下の（４７）式の関係が成り立つ。

また、（４７）式を変形（行列表記）することで、以下の（４８）式が得られる。

（４８）式におけるＲ_jlはデータｙの自己相関と呼ばれるもので、以下の（４９）式で定義される値である。このときの｜ｊ−ｌ｜を時差という。

（４８）式を基に、以下の（５０）式を考える。（５０）式は、自己回帰モデルによる時刻ｋにおける物理量の予測値ｙ＾_kと、その予測値ｙ＾_kに対応する時刻ｋにおける物理量の値ｙ_kと、の誤差を最小化する条件から導出される方程式である。（５０）式は、ユール・ウォーカー（Ｙｕｌｅ−Ｗａｌｋｅｒ）方程式と呼ばれる。また、（５０）式は、自己回帰モデルの係数から成るベクトルを変数ベクトルとする線形方程式である。（５０）式における左辺の定数ベクトルは、時差が１からｍまでのデータｙの自己相関を成分とするベクトルである。以下の説明では、（５０）式における左辺の定数ベクトルを必要に応じて自己相関ベクトルと称する。また、（５０）式における右辺の係数行列は、時差が０からｍ−１までのデータｙの自己相関を成分とする行列である。以下の説明では、（５０）式における右辺の係数行列を必要に応じて自己相関行列と称する。

また、（５０）式における右辺の自己相関行列（Ｒ_jlで構成されるｍ×ｍの行列）を、以下の（５１）式のように、自己相関行列Ｒと表記する。

一般に、自己回帰モデルの係数を求める際には、（５０）式を係数αについて解くという方法が用いられる。（５０）式では、自己回帰モデルで導出される時刻ｋにおける物理量の予測値ｙ＾_kが、その時刻ｋにおける物理量の値ｙ_kにできるだけ近づくように係数αを導出する。よって、自己回帰モデルの周波数特性には、各時刻におけるデータｙの値ｙ_kに含まれる多数の周波数成分が含まれる。

そこで、本発明者らは、自己回帰モデルの係数αに乗算される自己相関行列Ｒに着目し、鋭意検討した。その結果、本発明者らは、自己相関行列Ｒの固有値の一部を用いて、データｙに含まれる高周波成分の影響を低減することができることを見出した。即ち、本発明者らは、低周波成分が強調されるように自己相関行列Ｒを書き換えることができることを見出した。

以下に、このことの具体例を説明する。
自己相関行列Ｒを特異値分解する。自己相関行列Ｒの要素は、対称である。従って、自己相関行列Ｒを特異値分解すると以下の（５２）式のように、直交行列Ｕと、対角行列Σと、直交行列Ｕの転置行列との積となる。

（５２）式の対角行列Σは、以下の（５３）式に示すように、対角成分が自己相関行列Ｒの固有値となる行列である。対角行列Σの対角成分を、σ₁₁、σ₂₂、・・・、σ_mmとする。また、直交行列Ｕは、各列成分ベクトルが自己相関行列Ｒの固有ベクトルとなる行列である。直交行列Ｕの列成分ベクトルを、ｕ₁、ｕ₂、・・・、ｕ_mとする。自己相関行列Ｒの固有ベクトルｕ_jに対する固有値がσ_jjという対応関係がある。自己相関行列Ｒの固有値は、自己回帰モデルによる時刻ｋにおける物理量の予測値ｙ＾_kの時間波形に含まれる各周波数の成分の強度を反映する変数である。

自己相関行列Ｒの特異値分解の結果から得られる対角行列Σの対角成分であるσ₁₁、σ₂₂、・・・、σ_mmの値は、数式の表記を簡略にするために降順とする。（５３）式に示す自己相関行列Ｒの固有値のうち、最大のものからｓ個の固有値を用いて、以下の（５４）式のように、行列Ｒ'を定義する。ｓは、１以上且つｍ未満の数である。本実施形態では、ｓは、予め定められる。行列Ｒ'は、自己相関行列Ｒの固有値のうちｓ個の固有値を用いて自己相関行列Ｒを近似した行列である。

（５４）式における行列Ｕ_sは、（５２）式の直交行列Ｕの左からｓ個の列成分ベクトル（使用される固有値に対応する固有ベクトル）により構成されるｍ×ｓ行列である。つまり、行列Ｕ_sは、直交行列Ｕから左のｍ×ｓの要素を切り出して構成される部分行列である。また、（５４）式におけるＵ_s ^Tは、Ｕ_sの転置行列である。Ｕ_s ^Tは、（５２）式の行列Ｕ^Tの上からｓ個の行成分ベクトルにより構成されるｓ×ｍ行列である。（５４）式における行列Σ_sは、（５２）式の対角行列Σの左からｓ個の列と、上からｓ個の行により構成されるｓ×ｓ行列である。つまり、行列Σ_sは、対角行列Σから左上のｓ×ｓの要素を切り出して構成される部分行列である。
行列Σ_sおよび行列Ｕ_sを行列要素で表現すれば、以下の（５５）式のようになる。

自己相関行列Ｒの代わりに行列Ｒ'を用いることで、（５０）式の関係式を、以下の（５６）式のように書き換える。

（５６）式を変形することで、係数αを求める式として、以下の（５７）式が得られる。（５７）式によって求められる係数αを用いて、（４５）式により、時刻ｋにおける物理量の予測値ｙ＾_kを算出するモデルが「修正自己回帰モデル」である。

ここでは、対角行列Σの対角成分であるσ₁₁、σ₂₂、・・・、σ_mmの値を降順とする場合を例に挙げて説明した。しかしながら、係数αの算出過程において対角行列Σの対角成分は降順である必要はない。その場合には、行列Ｕ_sは、直交行列Ｕから左のｍ×ｓの要素を切り出して構成される部分行列ではなく、使用される固有値に対応する列成分ベクトル（固有ベクトル）を切り出して構成される部分行列になる。また、行列Σ_sは、対角行列Σから左上のｓ×ｓの要素を切り出して構成される部分行列ではなく、修正自己回帰モデルの係数の決定に利用される固有値を対角成分とするように切り出される部分行列になる。

（５７）式は、修正自己回帰モデルの係数の決定に利用される方程式である。（５７）式の行列Ｕ_sは、自己相関行列Ｒの特異値分解により得られる直交行列Ｕの部分行列であって、修正自己回帰モデルの係数の決定に利用される固有値に対応する固有ベクトルを列成分ベクトルとする行列（第３の行列）である。また、（５７）式の行列Σ_sは、自己相関行列Ｒの特異値分解により得られる対角行列の部分行列であって、修正自己回帰モデルの係数の決定に利用される固有値を対角成分とする行列（第２の行列）である。（５７）式の行列Ｕ_sΣ_sＵ_s ^Tは、行列Σ_sと行列Ｕ_sとから導出される行列（第１の行列）である。

（５７）式の右辺を計算することにより、修正自己回帰モデルの係数αが求まる。以上、修正自己回帰モデルの係数αの導出方法の一例について説明した。ここでは、修正自己回帰モデルの基となる自己回帰モデルの係数の導出方法を、直感的に分かり易いように、時刻ｋにおける物理量の予測値ｙ＾_kに対して最小二乗法を用いる方法とした。しかしながら、一般的には確率過程という概念を用いて自己回帰モデルを定義し、その係数を導出する方法が知られている。その場合に、自己相関は、確率過程（母集団）の自己相関で表現される。この確率過程の自己相関は、時差の関数として表される。従って、本実施形態におけるデータｙの自己相関は、確率過程の自己相関を近似するものであれば他の計算式で算出した値に代えてもよい。例えば、Ｒ₂₂〜Ｒ_mmは、時差が０（ゼロ）の自己相関であるが、これらをＲ₁₁に置き換えてもよい。

（５３）式に示す自己相関行列Ｒから抽出する固有値の数ｓは、例えば、自己相関行列Ｒの固有値の分布から決定することができる。
ここでは、前述した修正自己回帰モデルの説明における物理量は、前後方向力になる。前後方向力の値は、鉄道車両の状態に応じて変動する。
そこで、まず、鉄道車両を軌道１６上で走行させて、前後方向力の測定値についてのデータｙを得る。得られたデータｙ毎に、（４９）式と（５１）式とを用いて自己相関行列Ｒを求める。この自己相関行列Ｒについて（５２）式で表される特異値分解を行うことによって自己相関行列Ｒの固有値を求める。図８は、自己相関行列Ｒの固有値の分布の一例を示す図である。図８では、輪軸１３ａにおける前後方向力Ｔ₁の測定値のデータｙのそれぞれについての自己相関行列Ｒを特異値分解して得られた固有値σ₁₁〜σ_mmを昇順に並べ替えて、プロットしている。図８の横軸は、固有値のインデックスであり、縦軸は、固有値の値である。

図８に示す例では、他よりも顕著に高い値をもつ固有値が１つある。また、前記顕著に高い値もつ固有値ほどではないが、他と比べると比較的大きな値を持ち０（ゼロ）と見なせない固有値が２つある。このことから、（５３）式に示す自己相関行列Ｒから抽出する固有値の数ｓとして、例えば、２または３を採用することができる。どちらを採用しても結果に顕著な差異は生じない。

第１の周波数調整部４０４は、データ取得部４０３で前後方向力の測定値のデータｙの時刻ｋにおける値ｙ_kが所定のサンプリング周期で取得されるたびに以下の処理を行う。
まず、第１の周波数調整部４０４は、前後方向力の測定値のデータｙと、予め設定されている数Ｍ、ｍと、に基づいて、（４９）式と（５１）式とを用いて自己相関行列Ｒを生成する。

次に、第１の周波数調整部４０４は、自己相関行列Ｒを特異値分解することで、（５２）式の直交行列Ｕおよび対角行列Σを導出し、対角行列Σから自己相関行列Ｒの固有値σ₁₁〜σ_mmを導出する。
次に、第１の周波数調整部４０４は、自己相関行列Ｒの複数の固有値σ₁₁〜σ_mmのうち、最大のものからｓ個の固有値σ₁₁〜σ_ssを、修正自己回帰モデルの係数αを求めるのに利用する自己相関行列Ｒの固有値として選択する。
次に、第１の周波数調整部４０４は、前後方向力の測定値のデータｙと、固有値σ₁₁〜σ_ssと、自己相関行列Ｒの特異値分解により得られた直交行列Ｕと、に基づいて、（５７）式を用いて、修正自己回帰モデルの係数αを決定する。

そして、第１の周波数調整部４０４は、修正自己回帰モデルの係数αと、前後方向力の測定値のデータｙと、に基づいて、（４５）式により、前後方向力の測定値のデータｙの時刻ｋにおける予測値ｙ＾_kを導出する。前後方向力の予測値ｙ＾_kの時系列データは、前後方向力の測定値のデータｙに含まれる低周波成分を抽出した時系列データになる。

図９は、前後方向力の測定値の時系列データ（測定値）と、前後方向力の予測値の時系列データ（計算値）の一例を示す図である。尚、本実施形態では、４つの前後方向力Ｔ₁〜Ｔ₄の測定値が得られる。即ち、前後方向力について４つのデータｙが得られる。図９では、これら４つのデータｙのそれぞれにおける測定値および計算値を示す。図９の横軸は、基準の時刻を０（ゼロ）とした場合の当該基準の時刻からの経過時間（秒）で前後方向力Ｔ₁〜Ｔ₄の測定時刻・計算時刻を表す。縦軸は、前後方向力Ｔ₁〜Ｔ₄（Ｎｍ）である。

図９において、輪軸１３ａにおける前後方向力Ｔ₁の計算値は、概ね１５秒〜３５秒において、バイアスがかかっている（即ち、他の時間よりも大きな値を示す）。この期間は、輪軸１３ａが曲線軌道を通過する期間に対応する。輪軸１３ｂにおける前後方向力Ｔ₂の計算値、輪軸１３ｃにおける前後方向力Ｔ₃の計算値、および輪軸１３ｄにおける前後方向力Ｔ₄の計算値についても、輪軸１３ａにおける前後方向力Ｔ₁の計算値と同様に、輪軸１３ｂ、１３ｃ、１３ｄが曲線軌道を通過する期間にバイアスがかかっている。

従って、図９において、輪軸１３ａ〜１３ｄにおける前後方向力Ｔ₁〜Ｔ₄の測定値から計算値を除けば、前後方向力Ｔ₁〜Ｔ₄の信号のうち、輪軸１３ａ〜１３ｄが曲線軌道を通過することに起因する低周波成分を除くことができる。即ち、図９において、輪軸１３ａ〜１３ｄにおける前後方向力Ｔ₁〜Ｔ₄の測定値から計算値を除けば、輪軸１３ａ〜１３ｄが曲線軌道を通過した場合の前後方向力Ｔ₁〜Ｔ₄として、輪軸１３ａ〜１３ｄが直線軌道を通過した場合と同等の前後方向力を得ることができる。

そこで、第１の周波数調整部４０４は、前後方向力の測定値ｙ_kの時系列データ（データｙ）から、前後方向力の予測値ｙ＾_kの時系列データを減算する。以下の説明では、前後方向力の測定値ｙ_kの時系列データ（データｙ）から、前後方向力の予測値ｙ＾_kの時系列データを減算した時系列データを、必要に応じて前後方向力の高周波成分の時系列データと称する。また、前後方向力の高周波成分の時系列データの各サンプリング時刻における値を、必要に応じて前後方向力の高周波成分の値と称する。

図１０は、前後方向力の高周波成分の時系列データの一例を示す図である。図１０の縦軸は、前後方向力Ｔ₁、Ｔ₂、Ｔ₃、Ｔ₄の高周波成分の時系列データを示す。即ち、図１０の縦軸に示す前後方向力Ｔ₁、Ｔ₂、Ｔ₃、Ｔ₄の高周波成分は、それぞれ、図９に示した、輪軸１３ａ、１３ｂ、１３ｃ、１３ｄにおける前後方向力Ｔ₁、Ｔ₂、Ｔ₃、Ｔ₄の測定値から計算値を減算することにより得られるものである。また、図１０の横軸は、図９の横軸と同様に、基準の時刻を０（ゼロ）とした場合の当該基準の時刻からの経過時間（秒）で前後方向力Ｔ₁〜Ｔ₄の測定時刻・計算時刻を表す。
第１の周波数調整部４０４は、以上のようにして、前後方向力Ｔ₁〜Ｔ₄の高周波成分の時系列データを導出する。

［フィルタ演算部４０５、Ｓ７０３］
フィルタ演算部４０５は、観測方程式を観測方程式記憶部４０２により記憶された観測方程式とし、状態方程式を状態方程式記憶部４０１により記憶された状態方程式として、カルマンフィルタにより、（４４）式に示す状態変数の推定値を決定する。このとき、フィルタ演算部４０５は、データ取得部４０３で取得された計測データのうち、前後方向力Ｔ₁〜Ｔ₄を除く計測データと、第１の周波数調整部４０４で生成された前後方向力Ｔ₁〜Ｔ₄の高周波成分の時系列データとを用いる。前述したように本実施形態では、計測データには、車体１１の左右方向における加速度の測定値、台車１２ａ、１２ｂの左右方向における加速度の測定値、および輪軸１３ａ〜１３ｄの左右方向における加速度の測定値が含まれる。輪軸１３ａ〜１３ｄにおける前後方向力Ｔ₁〜Ｔ₄については、データ取得部４０３で取得された計測データ（測定値）を用いずに、第１の周波数調整部４０４で生成された前後方向力Ｔ₁〜Ｔ₄の高周波成分の時系列データを用いる。

カルマンフィルタは、データ同化を行う手法の一つである。即ち、カルマンフィルタは、観測できる変数（観測変数）の測定値と推定値との差異が小さく（最小に）なるように、未観測の変数（状態変数）の推定値を決定する手法の一例である。フィルタ演算部４０５は、観測変数の測定値と推定値との差異が小さく（最小に）なるカルマンゲインを求め、そのときの未観測の変数（状態変数）の推定値を求める。カルマンフィルタにおいては、以下の（５８）式の観測方程式と、以下の（５９）式の状態方程式を用いる。
Ｙ＝ＨＸ＋Ｖ ・・・（５８）
Ｘ・＝ΦＸ＋Ｗ ・・・（５９）

（５８）式において、Ｙは、観測変数の測定値を格納するベクトルである。Ｈは、観測モデルである。Ｘは、状態変数を格納するベクトルである。Ｖは、観測ノイズである。（５９）式において、Ｘ・は、Ｘの時間微分を示す。Φは、線形モデルである。Ｗは、システムノイズである。尚、カルマンフィルタ自体は、公知の技術で実現できるので、その詳細な説明を省略する。
フィルタ演算部４０５は、（４４）式に示す状態変数の推定値を所定のサンプリング周期で決定することにより、（４４）式に示す状態変数の推定値の時系列データを生成する。

［第２の周波数調整部４０６］
第１の周波数調整部４０４により、前後方向力の測定値の時系列データに含まれる低周波成分の信号強度が十分に除去されていないと、フィルタ演算部４０５により生成される状態変数の推定値の時系列データに、鉄道車両が曲線軌道を走行することに起因する低周波成分の信号が残る虞がある。そこで、第２の周波数調整部４０６は、フィルタ演算部４０５により生成された状態変数の推定値の時系列データに含まれる低周波成分の信号強度を低減（好ましくは除去）する。尚、第１の周波数調整部４０４により、前後方向力の測定値の時系列データに含まれる低周波成分の信号強度が十分に除去されるように、（５３）式に示す自己相関行列Ｒから抽出する固有値の数ｓを定めることができる場合には、第２の周波数調整部４０６の処理は不要になる。

本実施形態では、第２の周波数調整部４０６は、第１の周波数調整部４０４と同様に、修正自己回帰モデルを用いて、状態変数の推定値の時系列データに含まれる低周波成分の信号強度を低減する。

第２の周波数調整部４０６は、所定のサンプリング周期で状態変数毎に以下の処理を行う。
ここでは、前述した修正自己回帰モデルの説明における物理量は、状態変数になる。即ち、状態変数のデータｙは、フィルタ演算部４０５により生成された状態変数の推定値の時系列データになる。状態変数の推定値は、何れも鉄道車両の状態に応じて変動する。
まず、第２の周波数調整部４０６は、状態変数の推定値のデータｙと、予め設定されている数Ｍ、ｍと、に基づいて、（４９）式と（５１）式とを用いて自己相関行列Ｒを生成する。

次に、第２の周波数調整部４０６は、自己相関行列Ｒを特異値分解することで、（５２）式の直交行列Ｕおよび対角行列Σを導出し、対角行列Σから自己相関行列Ｒの固有値σ₁₁〜σ_mmを導出する。
次に、第２の周波数調整部４０６は、自己相関行列Ｒの複数の固有値σ₁₁〜σ_mmのうち、最大のものからｓ個の固有値σ₁₁〜σ_ssを、修正自己回帰モデルの係数αを求めるのに利用する自己相関行列Ｒの固有値として選択する。ｓは、状態変数毎に予め定められる。例えば、鉄道車両を軌道１６上で走行させて、これまで説明してきたようにして各状態変数の推定値のデータｙを得る。そして、自己相関行列Ｒの固有値の分布を状態変数毎に個別に作成する。この自己相関行列Ｒの固有値の分布から、（５３）式に示す自己相関行列Ｒから抽出する固有値の数ｓを状態変数のそれぞれについて決定する。

次に、第２の周波数調整部４０６は、状態変数の推定値のデータｙと、固有値σ₁₁〜σ_ssと、自己相関行列Ｒの特異値分解により得られた直交行列Ｕと、に基づいて、（５７）式を用いて、修正自己回帰モデルの係数αを決定する。

そして、第２の周波数調整部４０６は、修正自己回帰モデルの係数αと、状態変数の推定値のデータｙと、に基づいて、（４５）式により、状態変数の推定値のデータｙの時刻ｋにおける予測値ｙ＾_kを導出する。状態変数の予測値ｙ＾_kの時系列データは、状態変数の推定値のデータｙに含まれる低周波成分を抽出した時系列データになる。

そして、第２の周波数調整部４０６は、状態変数の推定値のデータｙから、状態変数の予測値ｙ＾_kの時系列データを減算する。以下の説明では、状態変数の推定値のデータｙから、状態変数の予測値ｙ＾_kの時系列データを減算した時系列データを、必要に応じて状態変数の高周波成分の時系列データと称する。

［接触角算出部４０７、Ｓ７０５］
接触角算出部４０７は、データ取得部４０３で取得された、車輪１４ａ〜１４ｄにおける輪重の測定値の時系列データ、車輪１４ａ〜１４ｄにおける横圧の測定値の時系列データ、前後方向力Ｔ₁〜Ｔ₄の時系列データ、および鉄道車両の走行速度ｖの時系列データと、第２の周波数調整部４０６により導出された状態変数の高周波成分の時系列データのうち、輪軸１３ａ〜１３ｄの左右方向の速度ｙ_w1・〜ｙ_w4・および台車１２ａ、１２ｂのヨーイング方向における回動量（角変位）ψ_t1〜ψ_t2と、横クリープ係数ｆ₂と、軸箱支持装置の前後方向のバネ定数Ｋ_wxと、軸箱支持装置の左右方向における間隔の１／２の長さｂ₁とに基づいて、各輪軸１３ａ〜１３ｄの各車輪１４ａ〜１４ｄにおける接触角α^L ₁〜α^L ₄、α^R ₁〜α^R ₄を導出する。

ここで、添え字Ｒは、右側の車輪１４ａ〜１４ｄを表し、添え字Ｌは、左側の車輪を表す。α^L ₁、α^L ₂、α^L ₃、α^L ₄は、それぞれ、輪軸１３ａ、１３ｂ、１３ｃ、１３ｄに属する左側の車輪（車輪１４ａ、１４ｂ、１４ｃ、１４ｄの反対側にある車輪）を表し、α^R ₁、α^R ₂、α^R ₃、α^R ₄は、それぞれ、輪軸１３ａ、１３ｂ、１３ｃ、１３ｄに属する右側の車輪１４ａ、１４ｂ、１４ｃ、１４ｄを表す。

また、ここでは、表記の都合上、必要に応じて、輪軸１３ａ、１３ｂ、１３ｃ、１３ｄを表す添え字１、２、３、４を纏めてｉと表す。同様に、必要に応じて、台車１２ａ、１２ｂを表す添え字１、２を纏めてｊと表す。例えば、速度ｙ_w1・〜ｙ_w4・は、ｙ_wi・と表し、台車１２ａ、１２ｂのヨーイング方向における回動量（角変位）ψ_t1〜ψ_t2は、ψ_tjと表し、前後方向力Ｔ₁〜Ｔ₄は、Ｔ_iと表し、接触角α^L ₁〜α^L ₄、α^R ₁〜α^R ₄は、α^L _i、α^R _iと表す。

ここで、接触角αについて説明する。図１１は、接触角α（図１１（ａ））およびフランジ角ａ（図１１（ｂ））の一例を説明する図である。図１１では、軌道１６の右側の軌条と、右側の車輪１４Ｒとの接触部分を鉄道車両の走行方向（ｘ軸方向）に垂直に切った場合の断面を示す。
図１１（ａ）において、法線荷重Ｎは、車輪１４Ｒが軌条から受ける荷重の、車輪１４Ｒと軌条との接触位置での接平面に対する法線方向の成分である。横クリープ力Ｆ_yは、クリープ力の、鉄道車両の走行方向（ｘ軸方向）に垂直な方向のうち車輪１４Ｒと軌条との接触位置での接平面に沿う方向の成分である。尚、縦クリープ力は、クリープ力の、鉄道車両の走行方向（ｘ軸方向）の成分である。法線荷重Ｎと横クリープ力Ｆ_yの合力が、車輪１４Ｒと軌条との接触位置で軌条から車輪１４に作用する力である。即ち、輪重Ｐ、横圧Ｑは、それぞれ、車輪１４と軌条との接触位置で軌条から車輪１４Ｒに作用する力の、高さ方向（ｚ軸方向）の成分、左右方向（ｙ軸方向）の成分である。接触角αは、車輪１４Ｒと軌条との接触位置での接平面と水平面（左右方向（ｙ軸方向）の面）とのなす角度のうち小さい方の角度（鋭角）である。図１１（ｂ）において、フランジ角ａは、車輪のフランジの踏面側（軌条側）の領域と水平面（左右方向（ｙ軸方向）の面）とのなす角度のうち小さい方の角度（鋭角）である。尚、フランジ角は、フランジ角度とも称される。

次に、フランジ接触について説明する。図１２は、フランジ接触の一例を説明する図である。図１２では、軌道１６の左右の軌条と、１つの輪軸１３とを鉄道車両の走行方向（ｘ軸方向）に垂直に切った場合の断面を示す。また、図１２では、軌道１６（軌条）が右方向（ｙ軸の負の方向）に曲がっており、鉄道車両が右方向に曲がりながら走行している場合における輪軸１３の様子を示す。尚、図１２では、左側の車輪１４Ｌにおける横クリープ力Ｆ_y ^L _iおよび法線荷重Ｎ^L _iと、右側の車輪１４Ｒにおける横クリープ力Ｆ_y ^R _iおよび法線荷重Ｎ^R _iを併せて示す。

図１２に示すように、鉄道車両が右方向に曲がっている軌条を走行する場合に、鉄道車両は、左方向（ｙ軸の正の方向）の作用力を受け、輪軸１３が左方向に移動することにより、車輪１４Ｌ、１４Ｒと軌条との接触位置からの左右方向の反力が大きくなり力の釣り合い点に達する。この作用力がさらに大きくなると、輪軸１３は左へさらに移動し、接触角α^Lが左側の車輪１４Ｌのフランジ角ａ^Lと同じになると、図１２に示すように、左側の車輪１４Ｌは、軌条とフランジで接触することになる。このような接触をフランジ接触という。一方、この状態では、右側の車輪１４Ｒは、軌条と踏面で接触する。

本実施形態では、接触角算出部４０７は、時々刻々と変化する接触角α^L _i、α^R _iを、輪重Ｐ、横圧Ｑ、輪軸１３ａ〜１３ｄの左右方向の速度ｙ_wi・、台車１２ａ、１２ｂのヨーイング方向における回動量（角変位）ψ_tj、前後方向力Ｔ_i、および鉄道車両の走行速度ｖの測定値を用いて算出（推定）する。

以下に本実施形態における接触角α^L _i、α^R _iの推定方法の一例を説明する。図１３は、左右の車輪１４Ｌ、１４Ｒに作用する力と接触角α^L _i、α^R _iの一例を示す図である。図１３でも、図１２と同様に、軌道１６の左右の軌条と、１つの輪軸１３とを鉄道車両の走行方向（ｘ軸方向）に垂直に切った場合の断面を示す。また、図１３では、説明の都合上、力の方向を分かりやすくするために、車輪１４Ｌ、１４Ｒと軌条の形状を単純化している。尚、図１３において、車軸１５の中央に示す矢印線は、車軸１５に働く力の方向を示す。

横圧Ｑの方向に関しては、ｙ軸の正の方向を正の方向とし、輪重Ｐの方向に関してはｚ軸の正の方向を正の方向とし、法線荷重Ｎの方向に関しては、当該法線荷重Ｎが作用する方向を正の方向とし、横クリープ力Ｆ_yの方向に関しては、鉄道車両の進行方向（ｘ軸の正の方向）を見て法線荷重Ｎの正の方向を時計回りに９０°回転させた方向を正の方向とする。尚、図１３において、左側の車輪１４Ｌにおける横圧Ｑ^L _iは、右側の車輪１４Ｒの方に働くが、図１３には横圧Ｑの正の方向を示すように描いている。そうすると、左側の車輪１４Ｌにおける横圧Ｑ^L _iは、左側の車輪１４Ｌにおける法線荷重Ｎ^L _i、輪重Ｐ^L _i、横クリープ力Ｆ_y ^L _i、および接触角α^L _iとの関係式として、以下の（６０）式および（６１）式で表される。同様に、右側の車輪１４Ｒにおける横圧Ｑ^R _iは、右側の車輪１４Ｒにおける法線荷重Ｎ^R _i、輪重Ｐ^R _i、横クリープ力Ｆ_y ^R _i、および接触角α^R _iとの関係式として、以下の（６２）式および（６３）式で表される。

尚、（６０）式〜（６３）式は、車輪１４Ｌ、１４Ｒと軌条との接触位置における力のバランスをとるための式である。
本実施形態では、鉄道車両が２１自由度を有する。この場合、各輪軸１３ａ、１３ｂ、１３ｃ、１３ｄにおける横圧Ｑ_iは、当該輪軸に属する左右の車輪１４Ｌ、１４Ｒにおける横圧Ｑ^L _i、Ｑ^R _iの和として表現できるので、以下の（６４）式が成り立つ。

ここで、横圧Ｑ^L _i、Ｑ^R _iの和を、フィルタ演算部４０５により導出されない状態変数である横クリープ力Ｆ_y ^L _i、Ｆ_y ^R _iおよび法線荷重Ｎ^L _i、Ｎ^R _iを用いずに表現するために、本実施形態では、（６４）式を以下のように変形する。
まず、横クリープ力Ｆ_y ^L _i、Ｆ_y ^R _iに関しては、近似的に以下の（６５）式で表されることが知られている。

前述したように、前後方向力Ｔ_iは、輪軸のヨーイング方向の角変位ψ_wiと、当該輪軸が設けられる台車のヨーイング方向の角変位ψ_tiとの差に応じて定まり、以下の（６６）式および（６７）式で表される。尚、（６６）式、（６７）式は、それぞれ、（４０）式〜（４３）式、（２６）式〜（２９）式を、ｉ、ｊを用いて表現したものである。

ここで、通常は、軸箱支持装置の前後方向のダンピング定数Ｃ_wxは小さく無視できるので、（６６）式および（６７）式から、以下の（６８）式が得られる。

（６８）式を（６５）式に代入すると、以下の（６９）式が得られる。

（６９）式を（６４）式に代入すると、以下の（７０）式が得られる。

（７０）式の法線荷重Ｎ^L _i、Ｎ^R _iに、（６０）式および（６２）式を代入すると、以下の（７１）式が得られる。

（７１）式を展開し、三角関数の積和公式を適用すると、以下の（７２）式が得られる。

以上のように、横圧Ｑ^L _i、Ｑ^R _iは、輪重Ｐ^L _i、Ｐ^R _i、接触角α^L _i、α^R _i、前後方向力Ｔ_i、輪軸１３ａ〜１３ｄの左右方向における速度ｙ_wi・、台車１２ａ、１２ｂのヨーイング方向における角回動量（角変位）ψ_tj、および鉄道車両の走行速度ｖを用いて表現される。
（７２）式において、横クリープ係数ｆ₂、軸箱支持装置の前後方向のバネ定数Ｋ_wxおよび軸箱支持装置の左右方向における間隔の１／２の長さｂ₁は定数である。尚、前述したように、横クリープ係数ｆ₂は、輪軸１３ａ〜１３ｄ毎に与えられるものとしてもよい。そうすると、（７２）式において、未知数は、接触角α^L _i、α^R _iになる。従って、未知数が２つになり、（７２）式だけでは、接触角α^L _i、α^R _iを算出することはできない。そこで、本発明者らは、車輪１４ａ〜１４ｄが軌条に対しフランジ接触している場合とフランジ接触していない場合とに分けて、以下のようにして接触角α^L _i、α^R _iの算出式を定めることを見出した。

≪フランジ接触していない場合の接触角α^L _i、α^R _iの算出式≫
鉄道車両が直線軌道や緩やかな曲線軌道を走行しており、車輪１４ａ〜１４ｄが軌条に対しフランジ接触していない場合、（概要）の項で説明したように、同一の輪軸１３ａ〜１３ｄに属する左右の車輪における接触角α^L _i、α^R _iの和は一定値ｃとみなす。一定値ｃは、例えば、３°程度の小さな正値に定めればよい。

接触角α^L _i、α^R _iの推定誤差が０（ゼロ）である場合、（７２）式の右辺の値と左辺の値との差は０（ゼロ）になる。同様に、（７３）式の右辺の値と左辺の値との差も０（ゼロ）になる。そこで、（７２）式の右辺の値と左辺の値との差を含むコスト関数Ｊ₁（α^L _i，α^R _i）を、以下の（７４）式とすると共に、（７３）式の右辺の値と左辺の値との差を含むコスト関数Ｊ₂（α^L _i，α^R _i）を、以下の（７５）式とするコスト関数Ｊ（α）＝（Ｊ₁，Ｊ₂）を考える。ここで、αは、α＝（α^L _i，α^R _i）である。

コスト関数Ｊ（α）の最小値は０（ゼロ）であるから（７６）式のように表現される。（７６）式において、コスト関数Ｊ（α）を或るα_kの周囲でテーラー展開すると、以下の（７７）式が得られる。従って、以下の（７８）式が得られる。

（７８）式のαをα_k+1（α＝α_k+1）とおき、（７８）式のＪ´（α_k）^-1Ｊ（α_k）に緩和係数βをかけて（７４）式および（７５）式を代入すると、以下の（７９）式が得られる。

ここで、緩和係数βは、予め与えられるものである。緩和係数βとして、例えば０．５を採用することができる。尚、緩和係数βを１としてもよい（緩和係数βは、式において、あらわになっていなくてもよい）。

≪フランジ接触している場合の接触角α^L _i、α^R _iの算出式≫
まず、鉄道車両が右方向に曲がりながら走行している場合の接触角α^L _iの算出式を説明する。
（概要）の項で説明したように、鉄道車両が急な曲線軌道（急カーブ）を走行しており、車輪１４ａ〜１４ｄが軌条に対しフランジ接触している場合、左右方向のうち、軌道１６の曲線部の曲がっている方向の車輪１４ａ〜１４ｄにおける接触角α^L _iまたはα^R _iは下限値（０（ゼロ）または０（ゼロ）に近い定数値）とみなす。従って、鉄道車両が右方向に曲がりながら走行している場合には、以下の（８０）式を満たすものとする。

ここで、α＾（式において＾はαの上に付される）は、接触角α^L _i、α^R _iの下限値（０（ゼロ）または０（ゼロ）に近い定数値）である。

まず、（７２）式において、接触角α^R _iをα＾とすると、（７２）式は、以下の（８１）のように、未知数α^L _iの非線形方程式に書き換えられる。

そうすると、（７４）式のコスト関数Ｊ₁（α^L _i，α^R _i）は、以下の（８２）式のコスト関数Ｊ₁（α^L _i）ように書き換えられる。

（８２）式をα^L _iで微分することにより以下の（８３）式が得られる。

（７４）式〜（７９）式を参照しながら説明したのと同様に、（８２）式の右辺を０（ゼロ）として、コスト関数Ｊ₁（α^L _i）を或る（α^L _i）_kの周囲でテーラー展開した式において、（α^L _i）を（α^L _i）_k+1とし、緩和係数をβとすると、以下の（８４）式が得られる。

次に、鉄道車両が左方向に曲がりながら走行している場合の接触角α^R _iの算出式を示す。
鉄道車両が左方向に曲がりながら走行している場合には、以下の（８５）式を満たすものとする。従って、（８２）式、（８３）式、（８４）式に対応する式として、それぞれ、以下の（８６）式、（８７）式、（８８）式が得られる。

≪接触角α^L _i、α^R _iの算出方法≫
（７９）式、（８４）式、（８８）式のようにして非線形方程式の解（ここではα^L _i、α^R _i）を求める手法は、ニュートン法と呼ばれている。（７９）式を例に挙げて説明すると、ニュートン法では、α_k（＝（α^L _i，α^R _i）_k）の初期値を（７９）式に与えて、（７９）式よりα_k+1（＝（α^L _i，α^R _i）_k+1）を求め、求めたα_k+1（＝（α^L _i，α^R _i）_k+1）を（７９）式のα_k（＝（α^L _i，α^R _i）_k）に与えて、（７９）式よりα_k+1（＝（α^L _i，α^R _i）_k+1）を求めることを、所定の収束条件を満たす（解が収束する）と判定されるまで繰り返す。

本発明者らは、このような収束計算を行うと、（７９）式、（８４）式、（８８）式の解が発散することがあることを見出した。特に、（８４）式、（８８）式の解が発散しやすくなる。これは、コスト関数Ｊ₁（α^L _i，α^R _i）、Ｊ₁（α^L _i）、Ｊ₁（α^R _i）において、フィルタ演算部４０５（第２の周波数調整部４０６）により導出された状態変数を用いていることと、（運動方程式）の項で説明した運動方程式が、車輪１４ａ〜１４ｄが軌条に対しフランジ接触していない場合の鉄道車両の直線軌道の走行時における運動を記述する運動方程式であること等のモデリング上の誤差要因があり、系の状態を一貫性をもって正確に記述することが困難であることに起因すると考えられる。

そこで、本発明者らは、（７９）式、（８４）式、（８８）式をフィルタとして用いることを見出した。即ち、（７９）式、（８４）式、（８８）式のα_k（＝（α^L _i，α^R _i）_k）、（α^L _i）_k、（α^R _i）_kを、前回のサンプリング時刻における接触角α^L _i、α^R _iとして、（７９）式、（８４）式、（８８）式により、α_k+1（＝（α^L _i，α^R _i）_k+1）、（α^L _i）_k+1、（α^R _i）_k+1を、今回のサンプリング時刻における接触角α^L _i、α^R _iとして算出する。このように、ニュートン法における解の更新の回数を１回に制限する。このようにすれば、（７９）式、（８４）式、（８８）式の解が発散することを抑制することができる。

また、このようにして算出される接触角α^L _i、α^R _iは、厳密な解ではない。しかしながら、本実施形態では、（４０）式〜（４３）式等を参照しながら説明したように、前後方向力Ｔ₁〜Ｔ₄を用いて運動方程式を構成する。前後方向力Ｔ₁〜Ｔ₄は、鉄道車両が走行している軌道１６（直線軌道および曲線軌道）のそれぞれの位置で精度よく測定される。従って、前後方向力Ｔ₁〜Ｔ₄に対応する状態変数（台車１２ａ、１２ｂのヨーイング方向における角変位ψ_ti）の精度はそれ程低下しないと考えられる。
以上のことから、本実施形態では、（７９）式、（８４）式、（８８）式をフィルタとして用いることとした。

本実施形態では、接触角算出部４０７は、鉄道車両の走行位置に関わらず、（７９）式による接触角α^L _i、α^R _iの算出と、（８４）式または（８８）式による接触角α^L _i、α^R _iの算出との双方を行う。また、接触角算出部４０７は、（７９）式により算出された接触角α^L _i、α^R _iの少なくとも何れか一方が、下限値α＾以下になる場合に、車輪１４ａ〜１４ｄが軌条に対しフランジ接触していると見なせると判定し、そうでない場合に、車輪１４ａ〜１４ｄが軌条に対しフランジ接触していると見なせないと判定する。

この判定の結果、車輪１４ａ〜１４ｄが軌条に対しフランジ接触していると見なせない場合、接触角算出部４０７は、（７９）式を用いて算出した接触角α^L _i、α^R _iを採用し、出力部４０８に出力する。一方、接触角算出部４０７は、車輪１４ａ〜１４ｄが軌条に対しフランジ接触していると見なせる場合であって、鉄道車両が右方向に曲がる曲線軌道を走行している場合には、（８４）式を用いて算出した接触角α^L _iを採用し、出力部４０８に出力する。また、接触角算出部４０７は、車輪１４ａ〜１４ｄが軌条に対しフランジ接触していると見なせる場合であって、鉄道車両が左方向に曲がる曲線軌道を走行している場合には、（８８）式を用いて算出した接触角α^R _iを採用し、出力部４０８に出力する。

ここで、接触角算出部４０７は、曲線軌道の左方向および右方向の何れの方向に曲がりながら走行しているのかを、例えば、軌道データを用いて判定することができる。軌道データは、鉄道車両が走行する軌道の位置と当該位置における軌条の曲率１／Ｒとの関係を示すデータを含む。曲率１／Ｒが０（ゼロ）である位置は直線軌道である。曲率１／Ｒが０（ゼロ）でない位置は曲線軌道である。このようなデータとして、例えば、鉄道会社が路線ごとに管理しているデータを用いることができる。

例えば、軌道データにおいて、右方向に曲がる曲線軌道における曲率１／Ｒの値を正とし、左方向に曲がる曲線軌道における曲率１／Ｒの値を負とする場合、接触角算出部４０７は、鉄道車両の現在の位置に対応する曲率１／Ｒを軌道データから読み出し、読み出した曲率１／Ｒの値に基づいて、鉄道車両が曲線軌道の左方向および右方向の何れの方向に走行しているのかを判定する。この他、接触角算出部４０７は、横圧Ｑ^L _i、Ｑ^R _iや輪重Ｐ^L _i、Ｐ^R _i等の計測データに基づいて、鉄道車両が曲線軌道の左方向および右方向の何れの方向に走行しているのかを判定してもよい。

また、接触角算出部４０７は、軌道データにおける曲率１／Ｒの絶対値が予め設定されている値以上になる場合に、車輪１４ａ〜１４ｄが軌条に対しフランジ接触していると見なせると判定し、そうでない場合に、車輪１４ａ〜１４ｄが軌条に対しフランジ接触していると見なせないと判定してもよい。その他、接触角算出部４０７は、横圧Ｑ^L _i、Ｑ^R _iや輪重Ｐ^L _i、Ｐ^R _i等の計測データに基づいて、車輪１４ａ〜１４ｄが軌条に対しフランジ接触していると見なせるか否かを判定してもよい。

また、接触角算出部４０７は、（７９）式、（８４）式、（８８）式の計算に際し、輪軸１３ａ〜１３ｄの左右方向の速度ｙ_wi・、台車１２ａ、１２ｂのヨーイング方向における回動量（角変位）ψ_tjについては、第２の周波数調整部４０６で低周波成分の信号強度が低減されたものを用いる。また、接触角算出部４０７は、（７９）式、（８４）式、（８８）式の計算に際し、前後方向力Ｔ_iについては、第１の周波数調整部４０４で低周波成分の信号強度が低減されたものではなく、データ取得部４０３で取得されたものを用いる。

接触角算出部４０７は、以上の接触角α^L _i、α^R _iの算出を所定のサンプリング周期で行うことにより、接触角α^L _i、α^R _iの時系列データを生成する。尚、（８４）式、（８８）式を用いて、左方向、右方向に曲がる曲線軌道を鉄道車両が走行している場合、当該方向の車輪１４ａ〜１４ｄにおける接触角α^L _i、α^R _iは、α＾になる。

図１４は、以上のようにして算出される接触角α^L ₁、α^R ₁の時系列データと、車輪１４ａ〜１４ｄにおける横圧Ｑ^L ₁、Ｑ^R ₁の測定値の時系列データと、車輪１４ａ〜１４ｄにおける輪重Ｐ^L ₁、Ｐ^R ₁の測定値の時系列データと、軌道データの一例を示す図である。尚、前述したように軌道データは、鉄道車両が走行する軌道の位置と当該位置における軌条の曲率１／Ｒとの関係を示すデータを含むが、図１４では、鉄道車両が走行する軌道の位置を、鉄道車両が当該位置を走行した時間に変換して示す。また、図１４の横軸の時間は、基準の時刻を０（ゼロ）とした場合の当該基準の時刻からの経過時間（秒）を示す。

図１４では、曲率半径Ｒが１６０ｍであり、曲線軌道が左方向に曲がっている場合を例に挙げて示す。図１４に示す例では、軌条のフランジ角は６７．８°である。
図１４に示すように、曲線軌道において、曲率１／Ｒの絶対値が大きくなるにつれて接触角α^R ₁はフランジ角に近づき、概ね４８５秒〜４９３秒において、車輪１４ａがフランジ接触の状態であることが推定される。その後、曲率１／Ｒの絶対値が小さくなるにつれて接触角α^R ₁は０（ゼロ）に近づく。また、直線軌道においては、接触角α^L ₁、α^R ₁は概ね一定の小さな値を示す。このように、本実施形態の手法によれば、軌条と軌道１６の状態を反映するように接触角α^L _i、α^R _iをリアルタイムで算出することができることが分かる。これにより、脱線に至る予兆を接触角α^L _i、α^R _iはから判定できることや、車輪１４ａ〜１４ｄと軌道１６の摩耗状態の履歴を時系列で監視できること等の効果を得ることが期待できる。

［出力部４０８、Ｓ７０６］
出力部４０８は、接触角算出部４０７により算出（採用）された接触角α^L _i、α^R _iの情報を出力する。出力の形態としては、例えば、コンピュータディスプレイへの表示、外部装置への送信、およびの内部または外部の記憶媒体への記憶の少なくとも何れか１つを採用することができる。

＜まとめ＞
以上のように本実施形態では、接触角推定装置４００は、更新前の接触角α_k（＝（α^L _i、α^R _i）_k）においてコスト関数Ｊ（α）を接触角αで微分した値の逆数J´^-1（α_k）と、更新前の接触角α_kにおけるコスト関数の値Ｊ（α_k）と、緩和係数βとの積を、更新前の接触角α_kから減算した値で、更新後の接触角α_k+1を表す式であるフィルタを用いて、接触角α^L _i、α^R _iを算出する。コスト関数Ｊ（α）は、横圧Ｑ^L _i、Ｑ^R _iの和と、（接触角α^L _i、α^R _iを用いて表現される）法線荷重Ｎ^L _i、Ｎ^R _iおよび横クリープ力Ｆ_y ^L _i、Ｆ_y ^R _iの左右方向の成分の和との差を含む関数である。また、横クリープ力Ｆ_y ^L _i、Ｆ_y ^R _iの左右方向の成分については、横クリープ力Ｆ_y ^L _i、Ｆ_y ^R _iを用いずに、輪軸１３ａ〜１３ｄの左右方向の速度ｙ_wi・、台車１２ａ、１２ｂのヨーイング方向における回動量（角変位）ψ_tj、前後方向力Ｔ_i、鉄道車両の走行速度ｖを用いて表現する。

接触角推定装置４００は、軌条がフランジ接触していない場合には、接触角α^L _i、α^R _iの和が一定値ｃであるとし、鉄道車両が曲線軌道を走行し軌条がフランジ接触している場合には、左右方向のうち曲線軌道の曲がっている方向の接触角α^L _iまたはα^R _iをα＾（０（ゼロ）または０（ゼロ）に近い一定値）とする。接触角推定装置４００は、更新前の接触角α_kとして、前回のサンプリング時刻において算出した（更新後の）接触角α_k+1を与え、更新回数を１回として、今回のサンプリング時刻における（更新後の）接触角α_k+1を導出する。従って、鉄道車両における接触角α^L _i、α^R _iを特別な測定装置を用いることなく容易に推定することができる。

即ち、本実施形態では、センサとして歪ゲージ等の汎用のセンサを用いることができるので、特別なセンサを必要としない。従って、鉄道車両における接触角α^L _i、α^R _iを大きなコストをかけることなく推定することができる。また、特別なセンサを用いる必要がないので、営業車両に歪ゲージを取り付け、営業車に接触角推定装置４００を搭載することにより、営業車両の走行中に、鉄道車両における接触角α^L _i、α^R _iをリアルタイムで推定することができる。従って、検測車を走行させなくても、軌道１６の不整を検知することができる。ただし、検測車に歪ゲージを取り付け、検測車に接触角推定装置４００を搭載してもよい。

また、本実施形態では、接触角推定装置４００は、前後方向力の測定値のデータｙから、自己相関行列Ｒを生成し、自己相関行列Ｒを特異値分解して得られた固有値のうち、最大のものからｓ個の固有値を用いて、前後方向力の測定値のデータｙを近似する修正自己回帰モデルの係数αを決定する。従って、前後方向力の測定値のデータｙに含まれる低周波成分の信号が残り、高周波成分が残らないように、係数αを決定することができる。接触角推定装置４００は、時刻ｋにおける前後方向力の予測値ｙ＾_kを、このようにして係数αが定められた修正自己回帰モデルに、その時刻よりも前の時刻ｋ−ｌ（１≦ｌ≦ｍ）の前後方向力の測定値のデータｙを与えることにより算出する。従って、カットオフ周波数を予め想定することなく、前後方向力の測定値のデータｙから、鉄道車両の曲線軌道の走行に起因する低周波成分の信号を低減することができる。

また、本実施形態では、接触角推定装置４００は、データ取得部４０３で取得された計測データのうち、前後方向力Ｔ₁〜Ｔ₄を除く計測データと、第１の周波数調整部４０４で生成された前後方向力Ｔ₁〜Ｔ₄の高周波成分の時系列データと、をカルマンフィルタに与えて、状態変数（ｙ_w1・〜ｙ_w4・、ｙ_w1〜ｙ_w4、ｙ_t1・〜ｙ_t2・、ｙ_t1〜ｙ_t2、ψ_t1・〜ψ_t2・、ψ_t1〜ψ_t2、φ_t1・〜φ_t2・、φ_t1〜φ_t2、ｙ_b・、ｙ_b、ψ_b・、ψ_b、φ_b・、φ_b、ψ_y1、ψ_y2、φ_a1、φ_a2）を導出する。次に、接触角推定装置４００は、状態変数の推定値の時系列データに含まれる低周波成分の信号強度を低減（好ましくは除去）することにより、状態変数の高周波成分の値を算出する。従って、輪軸１３ａ〜１３ｄのヨーイングを記述する運動方程式として、輪軸１３ａ〜１３ｄの位置での通り狂い量ｙ_R1〜ｙ_R4を変数として含む運動方程式を用いて状態方程式を構成する必要がなくなる。これにより、軌道１６のモデルを作成する必要がなくなると共に状態変数の数を減らすことができる。本実施形態では、モデルの自由度を２１自由度から１７自由度に減らすことができると共に、状態変数の数を３８から３０に減らすことができる。また、前後方向力Ｔ₁〜Ｔ₄の分だけ、カルマンフィルタで用いる測定値が増える。

一方、前後方向力Ｔ₁〜Ｔ₄を使用しないで、（５）式〜（８）式の輪軸１３ａ〜１３ｄのヨーイングを記述する運動方程式を状態方程式に含めると、計算が不安定となり、推定結果が得られない場合がある。即ち、状態変数を選定しないと、計算が不安定となり、推定結果が得られない場合がある。また、仮に推定結果が得られたとしても、本実施形態の方法の方が、状態変数を選定しない方法に比べ、軌道１６の不整の検知精度が高くなる。本実施形態では、輪軸１３ａ〜１３ｄのヨーイングを記述する運動方程式を状態方程式に含めないことと、前後方向力の測定値を用いることとを実現しているからである。

＜変形例＞
本実施形態では、フィルタ演算部４０５においてカルマンフィルタを用いる場合を例に挙げて説明した。しかしながら、観測変数の測定値と推定値との誤差が最小または当該誤差の期待値が最小になるように状態変数の推定値を導出するフィルタ（即ち、データ同化を行うフィルタ）を用いていれば、必ずしもカルマンフィルタを用いる必要はない。例えば、粒子フィルタを用いてもよい。尚、観測変数の測定値と推定値との誤差としては、例えば、観測変数の測定値と推定値との二乗誤差が挙げられる。

また、本実施形態では、車体１１の状態を表す状態変数を含める場合を例に挙げて説明した。しかしながら、車体１１は、車輪１４ａ〜１４ｄと軌道１６との間の作用力（クリープ力）による振動の伝搬が最後に伝わる部分である。従って、例えば、車体１１においてその伝搬による影響が小さいと判断される場合には、車体１１の状態を表す状態変数を含めなくてもよい。このようにする場合、（１）式〜（２１）式の運動方程式のうち、（１５）式〜（１７）式の車体１１の横振動、ヨーイング、ローリングを記述する運動方程式と、（１８）式、（１９）式の台車１２ａに配置されたヨーダンパ、台車１２ｂに配置されたヨーダンパのヨーイングを記述する運動方程式は不要になる。また、（１）式〜（２１）式の運動方程式において、車体に関する状態量（添え字ｂを含む状態量）と、車体に関する状態量（添え字ｂを含む状態量）を含む｛｝内の値（例えば（２１）式の左辺第３項の｛φ_a2−φ_b｝）を０（ゼロ）にする。

また、本実施形態では、フィルタ演算部４０５で導出された状態変数に含まれる（輪軸１３ａ〜１３ｄの左右方向の速度ｙ_wi・、台車１２ａ、１２ｂのヨーイング方向における回動量（角変位）ψ_tj）を用いて、接触角α^L _i、α^R _iを算出する場合を例に挙げて説明した。しかしながら、必ずしもこのようにする必要はなく、例えば、これらを、歪ゲージやジャイロスコープを用いて測定してもよい。

また、本実施形態では、修正自己回帰モデルを用いる場合を例に挙げて説明した。しかしながら、必ずしも、修正自己回帰モデルを用いて、前後方向力測定値のデータｙから、鉄道車両の曲線軌道の走行に起因する低周波成分の信号を低減する必要はない。例えば、鉄道車両の曲線軌道の走行に起因する周波数帯を特定することができる場合には、ハイパスフィルタを用いて、前後方向力の測定値のデータｙから、鉄道車両の曲線軌道の走行に起因する低周波成分の信号を低減してもよい。

また、本実施形態では、（７９）式、（８４）式、（８８）式の計算における更新の回数が１回である場合を例に挙げて説明したが、この更新の回数は、２回であってもよい。また、輪軸１３ａ〜１３ｄの左右方向の速度ｙ_wi・、台車１２ａ、１２ｂのヨーイング方向における回動量（角変位）ψ_tjを高精度に求めることができれば、解が収束するまで（収束条件を満たすまで）更新してもよい。この場合、更新前の接触角α_k（＝（α^L _i、α^R _i）_k）に与える初期値は、前回のサンプリング時刻における（更新後の）接触角α_k+1でなくてもよい。

また、本実施形態では、台車１２ａ、１２ｂがボルスタレス台車である場合を例に挙げて説明した。しかしながら、台車１２ａ、１２ｂは、ボルスタレス台車に限定されない。この他、鉄道車両の構成要素、鉄道車両が受ける力、および鉄道車両の運動の方向等に応じて、運動方程式は、適宜書き換えられる。即ち、運動方程式は、本実施形態で例示したものに限定されない。

（第２の実施形態）
次に、第２の実施形態を説明する。
第１の実施形態では、鉄道車両に搭載した接触角推定装置４００が接触角α^L _i、α^R _iを算出する場合を例に挙げて説明した。これに対し、本実施形態では、接触角推定装置４００の一部の機能が実装されたデータ処理装置が、指令所に配置される。このデータ処理装置は、鉄道車両から送信される計測データを受信し、受信した計測データを用いて接触角α^L _i、α^R _iを算出する。このように、本実施形態では、第１の実施形態の接触角推定装置４００が有する機能を、鉄道車両と指令所とで分担して実行する。本実施形態と第１の実施形態とは、このことによる構成および処理が主として異なる。従って、本実施形態の説明において、第１の実施形態と同一の部分については、図１〜図１４に付した符号と同一の符号を付す等して詳細な説明を省略する。

図１５は、接触角推定システムの構成の一例を示す図である。図１５において、接触角推定システムは、データ収集装置１５１０ａ、１５１０ｂと、データ処理装置１５２０とを有する。図１５には、データ収集装置１５１０ａ、１５１０ｂおよびデータ処理装置１５２０の機能的な構成の一例も示す。尚、データ収集装置１５１０ａ、１５１０ｂおよびデータ処理装置１５２０のハードウェアは、例えば、図５に示すもので実現することができる。従って、データ収集装置１５１０ａ、１５１０ｂおよびデータ処理装置１５２０のハードウェアの構成の詳細な説明を省略する。

鉄道車両のそれぞれには、データ収集装置１５１０ａ、１５１０ｂが１つずつ搭載される。データ処理装置１５２０は、指令所に配置される。指令所は、例えば、複数の鉄道車両の運行を集中管理する。

＜データ収集装置１５１０ａ、１５１０ｂ＞
データ収集装置１５１０ａ、１５１０ｂは、同じもので実現することができる。データ収集装置１５１０ａ、１５１０ｂは、データ取得部１５１１ａ、１５１１ｂと、データ送信部１５１２ａ、１５１２ｂとを有する。

［データ取得部１５１１ａ、１５１１ｂ］
データ取得部１５１１ａ、１５１１ｂは、データ取得部４０３と同じ機能を有する。即ち、データ取得部１５１１ａ、１５１１ｂは、データ取得部４０３で取得する計測データと同じ計測データを取得する。具体的にデータ取得部１５１１ａ、１５１１ｂは、計測データとして、車体１１の左右方向における加速度の測定値の時系列データ、台車１２ａ、１２ｂの左右方向における加速度の測定値の時系列データ、輪軸１３ａ〜１３ｄの左右方向における加速度の測定値の時系列データ、鉄道車両の走行速度ｖの時系列データ、車輪１４ａ〜１４ｄにおける輪重の測定値の時系列データ、および車輪１４ａ〜１４ｄにおける横圧の測定値の時系列データを取得する。これらの測定値を得るための歪ゲージおよび演算装置は、第１の実施形態で説明したものと同じである。

［データ送信部１５１２ａ、１５１２ｂ］
データ送信部１５１２ａ、１５１２ｂは、データ取得部１５１１ａ、１５１１ｂで取得された計測データを、データ処理装置１５２０に送信する。本実施形態では、データ送信部１５１２ａ、１５１２ｂは、データ取得部１５１１ａ、１５１１ｂで取得された計測データを、無線通信により、データ処理装置１５２０に送信する。このとき、データ送信部１５１２ａ、１５１２ｂは、データ収集装置１５１０ａ、１５１０ｂが搭載されている鉄道車両の識別番号を、データ取得部１５１１ａ、１５１１ｂで取得された計測データに付加する。このようにデータ送信部１５１２ａ、１５１２ｂは、鉄道車両の識別番号が付加された計測データを送信する。

＜データ処理装置１５２０＞
［データ受信部１５２１］
データ受信部１５２１は、データ送信部１５１２ａ、１５１２ｂにより送信された計測データを受信する。この計測データには、当該計測データの送信元である鉄道車両の識別番号が付加されている。

［データ記憶部１５２２］
データ記憶部１５２２は、データ受信部１５２１で受信された計測データを記憶する。データ記憶部１５２２は、鉄道車両の識別番号ごとに計測データを記憶する。データ記憶部１５２２は、鉄道車両の現在の運行状況と、計測データの受信時刻とに基づいて、当該計測データの受信時刻における鉄道車両の位置を特定し、特定した位置の情報と当該計測データとを相互に関連付けて記憶する。尚、データ収集装置１５１０ａ、１５１０ｂが、鉄道車両の現在の位置の情報を収集し、取集した情報を計測データに含めてもよい。

［データ読み出し部１５２３］
データ読み出し部１５２３は、データ記憶部１５２２により記憶された計測データを読み出す。データ読み出し部１５２３は、データ記憶部１５２２により記憶された計測データのうち、オペレータにより指定された計測データを読み出すことができる。また、データ読み出し部１５２３は、予め定められたタイミングで、予め定められた条件に合致する計測データを読み出すこともできる。本実施形態では、データ読み出し部１５２３により読み出される計測データは、例えば、鉄道車両の識別番号および位置の少なくとも何れか１つに基づいて決定される。

状態方程式記憶部４０１、観測方程式記憶部４０２、第１の周波数調整部４０４、フィルタ演算部４０５、第２の周波数調整部４０６、接触角算出部４０７、および出力部４０８は、第１の実施形態と説明したものと同じである。従って、ここでは、これらの詳細な説明を省略する。尚、フィルタ演算部４０５は、データ取得部４０３で取得された計測データに代えてデータ読み出し部１５２３で読み出された計測データを用いて、（４４）式に示す状態変数の推定値を決定する。

＜まとめ＞
以上のように本実施形態では、鉄道車両に搭載されたデータ収集装置１５１０ａ、１５１０ｂは、計測データを収集してデータ処理装置１５２０に送信する。指令所に配置されたデータ処理装置１５２０は、データ収集装置１５１０ａ、１５１０ｂから受信した計測データを記憶し、記憶した計測データを用いて、接触角α^L _i、α^R _iを算出する。従って、第１の実施形態で説明した効果に加え、例えば、以下の効果を奏する。即ち、データ処理装置１５２０は、計測データを任意のタイミングで読み出すことにより、任意のタイミングで接触角α^L _i、α^R _iを算出することができる。また、データ処理装置１５２０は、同じ位置における接触角α^L _i、α^R _iの時系列的な変化を出力することができる。

＜変形例＞
本実施形態では、データ収集装置１５１０ａ、１５１０ｂからデータ処理装置１５２０に計測データを直接送信する場合を例に挙げて説明した。しかしながら、必ずしもこのようにする必要はない。例えば、クラウドコンピューティングを利用して接触角推定システムを構築してもよい。
その他、本実施形態においても、第１の実施形態で説明した種々の変形例を採用することができる。

また、第１の実施形態では、状態方程式記憶部４０１、観測方程式記憶部４０２、データ取得部４０３、第１の周波数調整部４０４、フィルタ演算部４０５、第２の周波数調整部４０６、接触角算出部４０７、および出力部４０８が１つの装置に含まれる場合を例に挙げて説明した。しかしながら、必ずしもこのようにする必要はない。状態方程式記憶部４０１、観測方程式記憶部４０２、第１の周波数調整部４０４、フィルタ演算部４０５、第２の周波数調整部４０６、接触角算出部４０７、および出力部４０８の機能を複数の装置で実現してもよい。

（その他の実施形態）
以上説明した本発明の実施形態は、コンピュータがプログラムを実行することによって実現することができる。また、前記プログラムを記録したコンピュータ読み取り可能な記録媒体及び前記プログラム等のコンピュータプログラムプロダクトも本発明の実施形態として適用することができる。記録媒体としては、例えば、フレキシブルディスク、ハードディスク、光ディスク、光磁気ディスク、ＣＤ−ＲＯＭ、磁気テープ、不揮発性のメモリカード、ＲＯＭ等を用いることができる。

また、以上説明した本発明の実施形態は、何れも本発明を実施するにあたっての具体化の例を示したものに過ぎず、これらによって本発明の技術的範囲が限定的に解釈されてはならないものである。すなわち、本発明はその技術思想、またはその主要な特徴から逸脱することなく、様々な形で実施することができる。

（請求項と実施形態との関係）
以下に、請求項の記載と実施形態の記載との関係の一例を説明する。尚、請求項の記載が実施形態の記載に限定されないことは、変形例等に記載した通りである。
＜請求項１、２＞
接触角は、例えば、接触角α^L _i、α^R _i（α^L ₁〜α^L ₄、α^R ₁〜α^R ₄）を用いることにより実現される。
データ取得手段は、例えば、データ取得部４０３を用いることにより実現される。
接触角算出手段は、例えば、接触角算出部４０７を用いることにより実現される。
接触角を未知数とする方程式を含む式からなる第１のフィルタは、例えば、（７９）式、（８４）式、（８８）式を用いることにより実現される。
前記左右方向の一方側および他方側の前記車輪のそれぞれにおける輪重および横圧は、輪重Ｐ^L _i、Ｐ^R _i（Ｐ^L ₁〜Ｐ^L ₄、Ｐ^R ₁〜Ｐ^R ₄）、横圧Ｑ^L _i、Ｑ^R _i（Ｑ^L ₁〜Ｑ^L ₄、Ｑ^R ₁〜Ｑ^R ₄）を用いることにより実現される。
前記鉄道車両の走行速度は、例えば、鉄道車両の走行速度ｖを用いることにより実現される。
前後方向力は、例えば、前後方向力Ｔ_i（Ｔ₁〜Ｔ₄）を用いることにより実現される。
前記接触角の関数であるコスト関数の更新前の前記接触角における微分値は、例えば、（７９）式の右辺第２項の１つ目の行列（（７８）式のＪ´（α_k））、（８４）式の右辺第２項のＪ₁´（α^L _i）_k、（８８）式の右辺第２項のＪ₁´（α^R _i）_kを用いることにより実現される。
前記コスト関数の更新前の前記接触角における値は、例えば、（７９）式の右辺第２項の２つ目の行列（（７８）式のＪ（α_k））、（８４）式の右辺第２項のＪ₁（α^L _i）_k、（８８）式の右辺第２項のＪ₁´（α^R _i）_kを用いることにより実現される。
緩和係数は、例えば、（７９）式、（８４）式、（８８）式のβを用いることにより実現される。
更新前の前記接触角は、例えば、（７９）式の右辺第１項の行列（（７８）式のα_k）、（８４）式の右辺第１項の（α^L _i）_k、（８８）式の右辺第１項の（α^R _i）_kを用いることにより実現される。
更新後の前記接触角は、例えば、（７９）式の左辺の行列（（７８）式のα）、（８４）式の左辺の（α^L _i）_k+1、（８８）式の左辺の（α^R _i）_k+1を用いることにより実現される。
第１のコスト関数は、例えば、（７４）式、（８２）式、（８６）式を用いることにより実現される（（６４）式から（７４）式、（８２）式、（８６）式が得られることも参照）。
前記車輪と前記軌条との接触位置での法線荷重の前記左右方向の成分の、前記左右方向の一方側および他方側の前記車輪のそれぞれにおける値は、例えば、（６４）式の右辺第１項および第２項を用いることにより実現される。
前記車輪と前記軌条との接触位置での横クリープ力の前記左右方向の成分の、前記左右方向の一方側および他方側の前記車輪のそれぞれにおける値は、例えば、（６４）式の右辺第３項および第４項を用いることにより実現される。
前記車輪と前記軌条との接触位置での前記左右方向の一方側および他方側の前記車輪のそれぞれにおける前記横圧を加算した値は、例えば、（６４）式の左辺（＝Ｑ^L _i＋Ｑ^R _i）を用いることにより実現される。
第１の接触角および第２の接触角は、例えば、接触角α^L _i、α^R _i（α^L ₁〜α^L ₄、α^R ₁〜α^R ₄）を用いることにより実現される（例えば、第１の接触角が接触角α^L _iである場合、第２の接触角は接触角α^R _iである）。
前記第１のコスト関数における前記横クリープ力は、前記横クリープ力を用いずに、前記輪軸の前記左右方向の速度と、前記台車のヨーイング方向における角変位と、前記前後方向力と、前記鉄道車両の走行速度と、を用いて表現されることは、例えば、（６９）式を用いることにより実現される。ここで、前記輪軸の前記左右方向の速度は、例えば、輪軸１３ａ〜１３ｄの左右方向の速度ｙ_wi・（ｙ_w1・〜ｙ_w4・）を用いることにより実現され、前記台車のヨーイング方向における角変位は、例えば、台車１２ａ、１２ｂのヨーイング方向における角速度ψ_tj（ψ_t1〜ψ_t2）を用いることにより実現される。
＜請求項３、４＞
前記コスト関数として前記第１のコスト関数のみを用いて構成される前記方程式は、例えば、（８４）式、（８８）式を用いることにより実現される。
前記鉄道車両が前記左右方向の前記他方側に曲がる前記軌道を走行している場合は、例えば、鉄道車両が右側に曲がる曲線軌道を走行している場合に対応する。この場合、前記第１の接触角は、例えば、接触角α^L _iに対応する（（８４）式を参照）。また、予め定められた下限値は、例えば、下限値α＾を用いることにより実現される（（８０）式を参照）。
前記鉄道車両が前記左右方向の前記一方側に曲がる前記軌道を走行している場合は、例えば、鉄道車両が左側に曲がる曲線軌道を走行している場合に対応する。この場合、前記第２の接触角は、例えば、接触角α^R _iに対応する（（８８）式を参照）。また、予め定められた下限値は、例えば、下限値α＾を用いることにより実現される（（８５）式を参照）。
＜請求項５＞
第２のコスト関数は、例えば、（７５）式を用いることにより実現される。
前記コスト関数として前記第１のコスト関数および前記第２のコスト関数を用いる場合には、前記方程式で未知数となる前記接触角を、前記第１の接触角および前記第２の接触角として、前記第１の接触角および前記第２の接触角を算出することは、例えば、（７９）式を用いて、α_k+1（＝（α^L _i，α^R _i）_k+1）を算出することにより実現される。
＜請求項６＞
前記コスト関数として前記第１のコスト関数および前記第２のコスト関数を用いて構成される前記方程式を用いて算出した前記第１の接触角および前記第２の接触角の双方が、予め定められた下限値以下にならない場合は、例えば、（７９）式により算出された接触角α^L _i、α^R _iの双方が、下限値α＾以下にならない場合（即ち、（７９）式により算出された接触角α^L _i、α^R _iの双方が、下限値α＾よりも大きい場合）に対応する。
当該第１の接触角および当該第２の接触角を採用することは、例えば、接触角算出部４０７が、（７９）式を用いて算出した接触角α^L _i、α^R _iを採用し、出力部４０８に出力することにより実現される。
そうでない場合には、前記コスト関数として前記第１のコスト関数のみを用いて構成される前記方程式を用いて算出した前記第１の接触角または前記第２の接触角を採用することは、例えば、接触角算出部４０７が、（８４）式または（８８）式を用いて算出した接触角α^L _i、α^R _iを採用し、出力部４０８に出力することにより実現される。
＜請求項７＞
フィルタ演算手段は、例えば、フィルタ演算部４０５を用いることにより実現される。
状態方程式は、例えば、（９）式、（１０）式、（１３）式〜（２１）式、（３４）式〜（３９）式の運動方程式を用いることにより実現される。
観測方程式は、例えば、（３４）式〜（３７）式、（９）式、（１０）式、および（１５）式の横振動を記述する運動方程式を用いることにより実現される。
前記台車および前記輪軸の左右方向の加速度は、例えば、台車１２ａ〜１２ｂの左右方向における加速度ｙ_t1・・〜ｙ_t2・・、輪軸１３ａ〜１４ｄの左右方向における加速度ｙ_w1・・〜ｙ_w4・・を用いることにより実現される。
前記前後方向力が、前記輪軸のヨーイング方向の角変位と、当該輪軸が設けられる前記台車のヨーイング方向の角変位との差に応じて定まる力であることは、例えば、前後方向力Ｔ₁〜Ｔ₄が、（２２）式〜（２４）式を用いて表されることにより実現される。
前記台車の前記左右方向の変位および速度は、例えば、台車１２ａ〜１２ｂの左右方向における変位ｙ_t1〜ｙ_t2、速度ｙ_t1・〜ｙ_t2・を用いることにより実現される。
前記台車のヨーイング方向の角変位および角速度は、例えば、台車１２ａ〜１２ｂのヨーイング方向における角変位ψ_t1〜ψ_t2および角速度ψ_t1・〜ψ_t2・を用いることにより実現される。
前記台車のローリング方向の角変位および角速度は、例えば、台車１２ａ〜１２ｂのローリング方向における角変位φ_t1〜φ_t2および角速度φ_t1・〜φ_t2・を用いることにより実現される。
前記輪軸の前記左右方向の変位および速度は、例えば、輪軸１３ａ〜１３ｄの左右方向における変位ｙ_w1〜ｙ_w4および速度ｙ_w1・〜ｙ_w4・を用いることにより実現される。
前記鉄道車両に取り付けられている空気バネのローリング方向の角変位は、例えば、台車１２ａ〜１２ｂに配置された空気バネ（枕バネ）のローリング方向における回動量（角変位）φ_a1〜φ_a2を用いることにより実現される。
変換変数は、例えば、変換変数ｅ₁〜ｅ₄を用いることにより実現される（（２６）式〜（２９）式を参照）。
＜請求項８＞
第１の周波数調整手段は、例えば、第１の周波数調整部４０４を用いることにより実現される。
＜請求項９＞
第２の周波数調整手段は、例えば、第２の周波数調整部４０６を用いることにより実現される。

１１：車体、１２ａ〜１２ｂ：台車、１３ａ〜１３ｄ：輪軸、１４ａ〜１４ｄ：車輪、１５ａ〜１５ｄ：車軸、１６：軌道、４００：接触角推定装置、４０１：状態方程式記憶部、４０２：観測方程式記憶部、４０３：データ取得部、４０４：第１の周波数調整部、４０５：フィルタ演算部、４０６：第２の周波数調整部、４０７：接触角算出部、４０８：出力部

Claims (11)

1. 車体と台車と輪軸とを有する鉄道車両の車輪と軌条との接触位置における接平面と左右方向とのなす角度である接触角を推定する接触角推定システムであって、
   前記鉄道車両を軌道上で走行させることにより測定される測定値のデータである計測データを取得するデータ取得手段と、
   前記データ取得手段により取得された前記計測データと、前記接触角を未知数とする方程式からなる第１のフィルタとを用いて、前記接触角を算出する接触角算出手段と、を有し、
   前記計測データは、前記左右方向の一方側および他方側の前記車輪のそれぞれにおける輪重および横圧の測定値のデータと、前記鉄道車両の走行速度の測定値のデータと、前後方向力の測定値のデータとを含み、
   前記前後方向力は、前記輪軸と、当該輪軸が設けられる前記台車との間に配置される部材に生じる前後方向の力であり、
   前記部材は、軸箱を支持するための部材であり、
   前記前後方向は、前記鉄道車両の走行方向に沿う方向であり、
   前記左右方向は、前記前後方向と、前記軌道に対し垂直な方向である上下方向との双方に垂直な方向であり、
   前記方程式は、前記接触角の関数であるコスト関数の更新前の前記接触角における微分値と、前記コスト関数の更新前の前記接触角における値と、緩和係数と、更新前の前記接触角と、更新後の前記接触角との関係を表す式を含み、
   前記コスト関数は、前記車輪と前記軌条との接触位置での法線荷重の前記左右方向の成分の、前記左右方向の一方側および他方側の前記車輪のそれぞれにおける値と、前記車輪と前記軌条との接触位置での横クリープ力の前記左右方向の成分の、前記左右方向の一方側および他方側の前記車輪のそれぞれにおける値とを加算した値と、前記車輪と前記軌条との接触位置での前記左右方向の一方側および他方側の前記車輪のそれぞれにおける前記横圧を加算した値と、の差を含む第１のコスト関数を含み、
   前記方程式における前記接触角は、前記左右方向の一方側の前記車輪と前記軌条との接触位置における接平面と前記左右方向とのなす角度である第１の接触角と、前記左右方向の他方側の前記車輪と前記軌条との接触位置における接平面と前記左右方向とのなす角度である第２の接触角とを含み、
   前記方程式において未知数となる前記接触角は、前記第１の接触角と前記第２の接触角との少なくとも何れか一方であり、
   前記横圧、前記輪重、前記法線荷重、および前記横クリープ力は、前記車輪と前記軌条との接触位置における値であり、
   前記第１のコスト関数における前記横クリープ力は、前記横クリープ力を用いずに、前記輪軸の前記左右方向の速度と、前記台車のヨーイング方向における角変位と、前記前後方向力と、前記鉄道車両の走行速度と、を用いて表現され、
   前記ヨーイング方向は、前記上下方向を回動軸とする回動方向であり、
   前記接触角算出手段は、前記第１のフィルタに含まれる前記方程式に前記更新前の接触角を与えて前記更新後の前記接触角を算出することを、所定の条件が成立するまで行うことを特徴とする接触角推定システム。
2. 前記所定の条件は、更新の回数が１回であるという条件であり、
   前記接触角算出手段は、前回のタイミングで前記第１のフィルタを用いて算出した前記更新後の接触角を前記更新前の接触角として前記第１のフィルタに含まれる前記方程式に与えることにより前記更新後の接触角を今回のタイミングにおける前記接触角として算出することを特徴とする請求項１に記載の接触角推定システム。
3. 前記接触角算出手段は、前記コスト関数として前記第１のコスト関数のみを用いて構成される前記方程式を用いる場合であって、前記鉄道車両が前記左右方向の前記他方側に曲がる前記軌道を走行している場合には、前記方程式で未知数となる前記接触角を、前記第１の接触角とし、前記方程式における前記第２の接触角を、予め定められた下限値で一定として、前記第１の接触角を算出し、
   前記コスト関数として前記第１のコスト関数のみを用いて構成される前記方程式を用いる場合であって、前記鉄道車両が前記左右方向の前記一方側に曲がる前記軌道を走行している場合には、前記方程式で未知数となる前記接触角を、前記第２の接触角とし、前記方程式における前記第１の接触角を、予め定められた下限値で一定として、前記第２の接触角を算出することを特徴とする請求項１または２に記載の接触角推定システム。
4. 前記接触角算出手段は、前記車輪が前記軌条にフランジ接触していると見なせる場合に、前記第１の接触角および前記第２の接触角のうちの一方を、前記方程式で未知数となる前記接触角として算出し、他方を、前記予め定められた下限値で一定とすることを特徴とする請求項３に記載の接触角推定システム。
5. 前記コスト関数は、前記第１の接触角と前記第２の接触角との和と、予め定められた一定値との差を含む関数である第２のコスト関数を更に含み、
   前記接触角算出手段は、前記コスト関数として前記第１のコスト関数および前記第２のコスト関数を用いる場合には、前記方程式で未知数となる前記接触角を、前記第１の接触角および前記第２の接触角として、前記第１の接触角および前記第２の接触角を算出することを特徴とする請求項１〜４の何れか１項に記載の接触角推定システム。
6. 前記接触角算出手段は、前記コスト関数として前記第１のコスト関数のみを用いて構成される前記方程式を用いる場合であって、前記鉄道車両が前記左右方向の前記他方側に曲がる曲線軌道を走行している場合には、前記方程式で未知数となる前記接触角を、前記第１の接触角とし、前記方程式における前記第２の接触角を、予め定められた下限値で一定として、前記第１の接触角を算出し、
   前記コスト関数として前記第１のコスト関数のみを用いて構成される前記方程式を用いる場合であって、前記鉄道車両が前記左右方向の前記一方側に曲がる前記軌道を走行している場合には、前記方程式で未知数となる前記接触角を、前記第２の接触角とし、前記方程式における前記第１の接触角を、予め定められた下限値で一定として、前記第２の接触角を算出し、
   前記コスト関数として前記第１のコスト関数および前記第２のコスト関数を用いて構成される前記方程式を用いて算出した前記第１の接触角および前記第２の接触角の双方が、予め定められた下限値以下にならない場合には、当該第１の接触角および当該第２の接触角を採用し、そうでない場合には、前記コスト関数として前記第１のコスト関数のみを用いて構成される前記方程式を用いて算出した前記第１の接触角または前記第２の接触角を採用することを特徴とする請求項５に記載の接触角推定システム。
7. 前記計測データと、状態方程式と、観測方程式と、を用いて、データ同化を行う第２のフィルタを用いた演算を行うことにより、前記状態方程式で推定値を決定すべき変数である状態変数の推定値を決定するフィルタ演算手段を更に有し、
   前記計測データは、前記台車および前記輪軸の左右方向の加速度の測定値を更に含み、
   前記前後方向力は、前記輪軸のヨーイング方向の角変位と、当該輪軸が設けられる前記台車のヨーイング方向の角変位との差に応じて定まる力であり、
   前記ヨーイング方向は、前記上下方向を回動軸とする回動方向であり、
   前記状態方程式は、前記状態変数と、前記前後方向力と、変換変数と、を用いて記述される方程式であり、
   前記状態変数は、前記台車の前記左右方向の変位および速度と、前記台車のヨーイング方向の角変位および角速度と、前記台車のローリング方向の角変位および角速度と、前記輪軸の前記左右方向の変位および速度と、前記鉄道車両に取り付けられている空気バネのローリング方向の角変位と、を含み、前記輪軸の前記ヨーイング方向の角変位および角速度を含まず、
   前記ローリング方向は、前記前後方向を回動軸とする回動方向であり、
   前記変換変数は、前記輪軸の前記ヨーイング方向の角変位と前記台車の前記ヨーイング方向の角変位とを相互に変換する変数であり、
   前記観測方程式は、観測変数と、前記変換変数と、を用いて記述される方程式であり、
   前記観測変数は、前記台車および前記輪軸の左右方向の加速度を含み、
   前記フィルタ演算手段は、前記観測変数の測定値と、前記前後方向力の値および前記変換変数の実績値を代入した前記状態方程式と、前記変換変数の実績値を代入した前記観測方程式と、を用いて、前記観測変数の測定値と推定値との誤差または当該誤差の期待値が最小になるときの前記状態変数の推定値を決定し、
   前記変換変数の実績値は、前記前後方向力の測定値を用いて導出され、
   前記接触角算出手段は、前記第１のコスト関数における、前記輪軸の前記左右方向の速度および前記台車のヨーイング方向の角変位として、前記フィルタ演算手段により決定された前記状態変数の推定値を用いることを特徴とする請求項１〜６の何れか１項に記載の接触角推定システム。
8. 前記前後方向力の測定値の時系列データから、前記鉄道車両が前記軌道の曲線部を走行することに起因して生じる低周波成分の信号強度を低減する第１の周波数調整手段を更に有し、
   前記変換変数の実績値は、前記第１の周波数調整手段により低周波成分の信号強度が低減された前記前後方向力の値を用いて導出され、
   前記第１のフィルタを構成する前記方程式、前記状態方程式、および前記観測方程式は、前記軌条の曲率半径を含まない式であることを特徴とする請求項７に記載の接触角推定システム。
9. 前記状態変数の推定値から、前記鉄道車両が前記軌道の曲線部を走行することに起因して生じる低周波成分の信号強度を低減する第２の周波数調整手段を更に有し、
   前記接触角算出手段は、前記第１のコスト関数における、前記輪軸の前記左右方向の速度および前記台車のヨーイング方向の角変位として、前記第２の周波数調整手段により低周波成分の信号強度が低減された前記状態変数の推定値を用いることを特徴とする請求項７または８に記載の接触角推定システム。
10. 車体と台車と輪軸とを有する鉄道車両の車輪と軌条との接触位置における接平面と左右方向とのなす角度である接触角を推定する接触角推定方法であって、
    前記鉄道車両を軌道上で走行させることにより測定される測定値のデータである計測データを取得するデータ取得工程と、
    前記データ取得工程により取得された前記計測データと、前記接触角を未知数とする方程式からなる第１のフィルタとを用いて、前記接触角を算出する接触角算出工程と、を有し、
    前記計測データは、前記左右方向の一方側および他方側の前記車輪のそれぞれにおける輪重および横圧の測定値のデータと、前記鉄道車両の走行速度の測定値のデータと、前後方向力の測定値のデータとを含み、
    前記前後方向力は、前記輪軸と、当該輪軸が設けられる前記台車との間に配置される部材に生じる前後方向の力であり、
    前記部材は、軸箱を支持するための部材であり、
    前記前後方向は、前記鉄道車両の走行方向に沿う方向であり、
    前記左右方向は、前記前後方向と、前記軌道に対し垂直な方向である上下方向との双方に垂直な方向であり、
    前記方程式は、前記接触角の関数であるコスト関数の更新前の前記接触角における微分値と、前記コスト関数の更新前の前記接触角における値と、緩和係数と、更新前の前記接触角と、更新後の前記接触角との関係を表す式を含み、
    前記コスト関数は、前記車輪と前記軌条との接触位置での法線荷重の前記左右方向の成分の、前記左右方向の一方側および他方側の前記車輪のそれぞれにおける値と、前記車輪と前記軌条との接触位置での横クリープ力の前記左右方向の成分の、前記左右方向の一方側および他方側の前記車輪のそれぞれにおける値とを加算した値と、前記車輪と前記軌条との接触位置での前記左右方向の一方側および他方側の前記車輪のそれぞれにおける前記横圧を加算した値と、の差を含む第１のコスト関数を含み、
    前記方程式における前記接触角は、前記左右方向の一方側の前記車輪と前記軌条との接触位置における接平面と前記左右方向とのなす角度である第１の接触角と、前記左右方向の他方側の前記車輪と前記軌条との接触位置における接平面と前記左右方向とのなす角度である第２の接触角とを含み、
    前記方程式において未知数となる前記接触角は、前記第１の接触角と前記第２の接触角との少なくとも何れか一方であり、
    前記横圧、前記輪重、前記法線荷重、および前記横クリープ力は、前記車輪と前記軌条との接触位置における値であり、
    前記第１のコスト関数における前記横クリープ力は、前記横クリープ力を用いずに、前記輪軸の前記左右方向の速度と、前記台車のヨーイング方向における角変位と、前記前後方向力と、前記鉄道車両の走行速度と、を用いて表現され、
    前記ヨーイング方向は、前記上下方向を回動軸とする回動方向であり、
    前記接触角算出工程は、前記第１のフィルタに含まれる前記方程式に前記更新前の接触角を与えて前記更新後の前記接触角を算出することを、所定の条件が成立するまで行うことを特徴とする接触角推定方法。
11. 請求項１〜９の何れか１項に記載の接触角推定システムの各手段としてコンピュータを機能させることを特徴とするプログラム。

Images