鉄道軌道に関する研究

営業車に搭載された軌道検測装置により取得した軌道変位データの活用を想定して，軌道力学状態のモニタリング手法に関する理論的検討を行った．具体的には，レール通り変位の変動に基づく，レール軸力と道床横抵抗力分布の推定法を提案した．その際に，通り変位データに含まれる測定ノイズに対処するため，粒子フィルタの援用を試みた．その結果，通常の通り変位振幅下での適用にはS/N比のさらなる改善が必要であるものの，比較的大きな通り変位の下では絶対軸力が推定可能であることを確認した．一方，道床横抵抗力分布については，レール軸力と道床横剛性とが適切に推定できれば，十分な精度で評価可能であることがわかった．

　左右レールは締結装置を介してまくらぎに接合されている．したがって，左右いずれか一方のレール作用力がわかっても道床横抵抗力を与えず，それらの合力を評価する必要がある．そのため，左右レール・まくらぎ連成系のつり合い状態について考える．
　締結部における横作用力は，レールとまくらぎとの相対変位に依存する．すると，左右レールの水平たわみに関するつり合い式は次式で与えられる．

　ここで，レールはEulerばりによりモデル化するものとし，Eはレールのヤング率，Iはレール弱軸回りの断面二次モーメント，Nは軸力(圧縮を正)，w_L, w_Rは左右レールの弾性たわみ，w_L0, w_R0はN=0における初期通り変位であり，( )'は軌道長手方向座標x に関する微分である．また，締結部のレール拘束力は本来離散的に作用するが，ここではそれを連続支持モデルにより近似表現している．なお，w_sはまくらぎ横変位，k_rは締結部の横剛性(単位長さ当たり)である．本研究では，左右レール軸力は等しいものとする．また，「通り変位」はたわみ波形を意味する量として用いる．
　道床横抵抗力がまくらぎ横変位に関して線形ばねk_sにより近似できるものと仮定すると，まくらぎのつり合い式は次式で与えられる．

式(2)をw_sについて解くと次式を得る．

式(3)を式(1)に代入してw_sを消去すると次式を得る．

さらに式(4)両式の和をとると次式を得る．

　一方，軸力がN+ΔNであるときのつり合い式は，たわみ増分Δw_LRを用い次式により与えられる．

式(6)から式(5)第１式を引くと次式を得る．

さらに式(7)のxに関するFourier変換より次式を得る．

ここでkは波数，(^)はレール長手方向に関するFourier変換を意味する．
　式(8)を整理して次式を得る．

式(9)左辺は，通り変位の測定データとその差分をFourier変換することで求めることができる．式(9)右辺より，これは次のように波数kの4次関数で与えられる．

なお，以下では10m弦正矢ではなく，通り変位源波形データを対象として定式化を進める．
　式(10)右辺の4次関数と，その係数a,b,c，および極小値を与える波数kmとの関係は，概略図1に示す様になる．なお，波数kmと絶対軸力N+ΔN，未知係数a,bらの関係は次式で与えられる．

式(9), (10)および(11)より次式を得る．

以上より，通り変位データのFourier変換より未知係数a,b,cを最小二乗法等により決定し，kmを求めれば，軸力N+ΔN(同様にN)および軌道横剛性を推定することができる．

図1 式(9)の4次関数とその最小値を与える波数km

2.1 数値モデル

　1.に述べた測定原理の正当性を検証するために，軌道モデルによる数値実験を実施した．軌道は弾性連続支持されたEulerばりとしてモデル化する．これをはり要素で離散化し，初期通り変位とレール軸力を設定して変形解析を行う．
　式(5)に対応する離散化求解方程式は次式で与えられる．

ここで，[K₁]，[K₂] は，それぞれEIw_LR''''，w_LR''を離散化して得られる行列である．また，{W}，{W₀}は，それぞれw_LRとw_LR0に関する節点ベクトルである．さらには，の項を離散化して得られる行列である．

2.2 解析結果

　50kgNレールを対象に，数値モデルの軌道長を200mとし，0.25mのはり要素で離散化した．通り変位は0.25m間隔で測定するものとし，はり要素各節点におけるたわみより擬似測定データを作成した．また，左右レールから構成される軌道系を考えた場合の軌道横剛性は2MN/m²と設定した．
　たわみ解のFourier変換よりを求めた一例を図2に示す．なお，図には式(9)右辺の理論曲線を実線で示した．数値実験結果には長波長成分に乱れが認められるが，

図2 数値実験結果(ノイズ無し)

これは有限長で軌道をモデル化したためと考えられる．スペクトル比には多少の変動が認められるものの，概ね理論どおり式(9)右辺の4次曲線が得られている．

　次に，式(12)より軸力を推定した結果を 図3に示す．なお，10m弦正矢データの場合，波数がk≒1.26(1/m)で感度がゼロとなる．このことを念頭に，ここではそれ以下の波数域データを用いている．図より，いずれの軸力においても，前述の理論によって極めて高い精度で軸力推定がなされており，当該理論の正当性が確認できる．

図3 軸力推定結果(ノイズ無し)

　続いて，全通り変位にσ_ε=0.1mmのガウスノイズを加えた場合を対象に，式(10)左辺のスペクトル比を求めた例を図4に示す．なお，図中の曲線はプロット点の最小二乗近似より得られた4次曲線である．当該曲線の最小点は正の波数域には無く，適切な軸力推定がなされていない．式(10)左辺より軸力推定する場合に，上述の問題を生ずる主な原因は，分母が通り変位増分のFourier変換で与えられていることにあると考えられる．この場合，当該量は微小値をとり，ノイズと同オーダーとなり得る．
　実際の軌道で得られる通り変位データには，さらに大きなノイズ(σ_ε≒0.5mm)[1]が含まれているものと考えられる．したがって，本手法を用いる場合，ノイズに対する適切な処理が不可欠である．

図4 数値実験結果σ_ε=0.1mm のノイズ有り)

3.1 軸力推定問題の設定

　ある時刻における通り変位をw_LR1，その際に作用している軸力をN₁とする．また，時間経過後における通り変位をw_LR2，そのときの軸力をN₂とすると，式(9)と同様に次式が成り立つ．

ここで，ΔN=N₂-N₁である．
　測定データは列車が同一軌道を複数回走行して得られる時系列データを計測量としている．よって，この2つのデータを同等に活用する目的で，軸力推定に用いる式を，式(14)両式の総和平均より得られる次式に修正する．

本推定問題を次式により設定する．

ここで，X_tは推定すべき未知量を成分とするベクトルであり，ε_tは通り変位測定データにおけるノイズによる当該Fourierスペクトル比への影響項である．なお，tは時間ステップである．2つの連続測定データに基づく推定を対象とする場合，X_t は次式で構成される．

y_tは通り変位測定データのFourier変換の比(スペクトル比)であり，関数G_t(X_t)は式(15)右辺で与えられるもので，未知量である軸力N_i (i=1, 2)，および軌道横剛性から次式により表すことができる．

粒子フィルタでは，まず推定値の候補となり得るような様々なX₀ = X_t を生成する．これは粒子と呼ばれ，このとき生成した粒子の個数を粒子数と呼ぶ．粒子フィルタはこれら粒子の集合により未知量の確率分布を近似する手法である．具体的には，まず各粒子毎に式(15)右辺を設定し，それと測定データより得られた同式左辺との差異から前者の尤度を評価する．続いて，これに比例する様に各粒子数を再配分する．なお，通常の時間発展問題の場合，この操作を繰り返して事後確率分布を更新するが，本問題では2つの軸力を同時推定することもあり，以降に示す解析例では1回の操作で概ね収束している．

3.2 数値モデルに基づく推定精度の検討

　左右レールの初期通り変位{W_0L}, {W_0R}を期待値ゼロ，標準偏差σのガウス分布に従うものとして生成する．また，その際に距離相関関数を次式により規定した．

　ここで，aはレール初期通り変位の相関長である．以下の解析例では，まくらぎ間隔を一つの目安として0.7mに設定した．
　以上のように作成した初期通り変位の下で弾性変位を求め，さらにノイズを加えて得られた擬似測定データをFourier変換して，前述の粒子フィルタによる軸力N₁, N₂と軌道横剛性の推定を試みた．初期状態の粒子は，軸力と軌道横剛性の推定範囲をそれぞれ0～500kN, 1.5～3.5MN/m²と設定して，各粒子におけるこれらの値を一様乱数により与えた．
　また，推定に用いる波数の範囲を0.5≤k≤ 2 (1/m)とした．ちなみに，式(15)右辺に正解値を設定して得られた曲線と，σ_ε=0.5mmのガウスノイズを含むデータから求めた同式左辺を図5に示す．なお，擬似測定データより求めたスペクトル比は，wavelet変換によりパルス状の変動成分を除去したものである．図4に示した様に，k≤ 1 (1/m)の波数域データではノイズの影響が大きく推定が困難であったが，k≤ 2 (1/m)までの波数範囲を用いれば，概ね適切に再現し得ることが確認できる．

図5 式(15)によるスペクトル比の再現と正解との比較

　前述のとおり，検測装置により入手されたデータは，通り変位を10m弦正矢に変換したものとなっている．この場合，波数が適切な精度で得られる範囲は概ねk≤ 1 (1/m)であり，k≒1.26 (1/m)以上の成分を破棄していることになる．しかし図5より，低波数成分は測定ノイズの影響が大きいため推定での使用は難しく，0.5≤k≤ 2 (1/m)の範囲のデータが有効であるとの結論を得た．以上より，10m弦正矢データではなく，通り変位データ自体を使用することが必要であると考えられる．
　粒子数を50000個，測定ノイズの標準偏差を0.5mmと設定した場合を対象に，初期通り変位の標準偏差を1～5cmまで変動させ，その差異が推定結果に及ぼす影響について調べた．正解値はN₁, N₂=100, 200kN, =2MN/m²と設定した．なお，事前の検討結果より，粒子数を10000個以上に設定すれば推定精度が概ね一定値に収束することを確認している．
　結果を図6～8に示す．初期粒子生成の際に乱数を用いているため，同じ擬似測定データに対して解析を行っても推定値は毎回異なる結果となる．そのため，推定を5回実施し，その平均値と標準偏差を図中に示している．また，図には正解値を水平線で示している．初期通り変位の標準偏差が1～4cmに対しては，明確な精度向上が認められないが，5cmのケースでは弾性たわみが比較的大きいこともあり，絶対軸力も含めた精度が顕著に改善されており，実際に存在し得るレベルのノイズ下であっても良好な推定結果を得ることができている．
　ただし，通常の状態における通り変位の標準偏差が5cmにまで達することは稀であると考えられる．図示はしていないが，実軌道での軸力推定のためにはノイズの標準偏差を0.1mm以下にする必要があり，測定精度のさらなる向上が望まれる．

図6 初期通り変位の標準偏差を変動させた場合のN₁推定値

図7 初期通り変位の標準偏差を変動させた場合のN₂推定値

図8 初期通り変位の標準偏差を変動させた場合の推定値

4.1 問題の設定

　これまでの検討結果より，レールの通り変位振幅が比較的大きな場合においては，現在の測定精度の下でも絶対軸力が推定可能であるとの結論を得た．本節では，絶対軸力と軌道横剛性とが良好な精度の下で推定された場合を前提として，道床横抵抗力分布を評価するための手法について検討する．そのために，まずは対象とする問題について改めて整理する．
　左右レールの横方向つり合い式は式(5)第1式で与えられる．レール締結剛性k_rは道床横剛性k_sに比べ一般に十分大きな値を持つ．この場合，は次式により近似することができる．

よって，軌道横剛性の推定結果より道床横剛性k_sを概ね推定可能である．
　式(20)を(5)に代入すると次式を得る．

式(21)をはり要素により離散化すると，次の連立方程式を得る．

初期通り変位w_LR0({W₀})が得られたならば，式(22)により弾性たわみw_LRを求めることができ，さらに次式より道床横抵抗力分布を算出することができる．

したがって，道床横抵抗力分布の推定は，初期通り変位の推定に帰着する．

4.2 目的関数の設定と勾配の評価

　式(22)に基づいた{W₀}の推定について考える．左右レール通り変位の和より与えられる測定データを{y_w}とおく．{y_w}と通り変位の節点ベクトルとの関係は次式で与えられる．

ここで，[B]は全節点変位ベクトルから測定点のたわみに対応する成分のみ抽出する行列である．{ε}は各測点における測定ノイズを成分とするベクトルであり，各成分は互いに独立であるものとする．
　{y_w}が与えられた時の，未知量{W₀}の事後確率密度関数を，以下の目的関数より求める．

ここで，( )_jは第j番測定データに関する値，{λ_j}は未定乗数ベクトルである．
　式(25)より，は次式で評価できる．

　最小点は式(26)がゼロとなる条件より求められる．これは{ΔW₀}の線形方程式となるが，{λ_j}が初期通り変位の陰的関数で与えられるため，これを直接解くことはむしろ煩雑となり得策ではない．そこで，式(26)より目的関数の勾配を求め，準Newton法により目的関数の最小点を求める．なお，その際にBFGS法(Broyden, Fletcher, Goldfarb, Shanno)により目的関数のHessian逆行列を求め，Armijoのルールを援用して修正ベクトル{ΔW₀}を設定する．

4.3 解析例

　初期通り変位の推定，およびそれを用いた道床横抵抗力分布評価に対して提案法の適用を試みた．軌道条件はこれまでのものと同じである．また，初期通り変位とノイズの標準偏差は，それぞれ1cmおよび0.5mmとしている．作成した擬似測定データを用いて前述の目的関数を設定し，その最小点探索を行い，初期通り変位を推定する．ただしその際に，軸力と道床横剛性は既に推定されているものとし，これらは正しい値に設定した．
　軸力をN=100, 150, 200, 250kNの4段階で設定し，各軸力における擬似測定データを用いて求めた初期通り変位推定結果の例を図9に，道床横抵抗力分布の推定結果を図10に示す．なお，図10は最終軸力(N=250kN)における道床横抵抗力を示したものである．{W₀}={0}を初期値としたが，2ステップ程で概ね正解に一致する初期通り変位を得ることができた．また，道床横抵抗力についても，正解との差異は比較的小さく，分布の特徴も良好に再現できている．
　以上より，少なくとも道床横剛性が線形モデルで近似し得る場合，初期通り変位および道床横抵抗力を良好な精度下で推定可能であることがわかる．

　通り変位の変化から絶対軸力と道床横剛性とを推定するための手法(軌道の力学状態のモニタリング手法)について，その可能性を検討した．その結果，初期通り変位が比較的大きく，よって弾性変位も比較的大きい場合，推定法として粒子フィルタを用いることで，現在の測定精度でも概ね推定可能であることがわかった．一方，初期通り変位がそれ程大きくない場合は，レールの弾性たわみも小さくなるため，良好な推定精度の確保には，測定データに対してより高いS/N比が要求される．また，絶対軸力の感度は決して高くないため，精度確保のためには，粒子数を相当多くとる必要がある．ただし，1粒子当りの計算負荷は非常に少なく，粒子数増加はそれ程負担とはならない．
　現在取得されているデータは，10m弦正矢に変換したものとなっている．この場合，波数k≒1.26(1/m)以上(波長5m以下)の成分を破棄していることになる．しかし本研究では，低波数成分は測定ノイズの影響が大きいため使用は難しく，k=0.5～2(1/m)程度の範囲のデータが有効であるとの結果を得た．したがって，10m弦正矢データでは無く，通り変位データ自体を使用することが必要である．
　軸力と道床横剛性とが良好な精度で推定できた場合を前提として，初期通り変位と道床横抵抗力分布の推定について検討した．その結果，何れも比較的良好に推定可能であることを確認した．

図9 初期通り変位の推定

図10 道床横抵抗力分布の推

参考文献

[1] 坪川洋友，矢沢栄治，小木曽清高，南木聡明 : 車体装架型慣性正矢軌道検測装置の開発，鉄道総研報告，Vol.26(2), 7-12, 2012. W
[2] 樋口知之，上野玄太，中野慎也，中村和幸，吉田　亮：データ同化入門，朝倉書店，2010．