筑波大学大学院博士課程
システム情報工学研究科修士論文

車両振動へのカルマンフィルタ適用による
VBI システム同定の可能性に関する
数値的検討

村上

翔

修士（工学）（構造エネルギー工学専攻）

指導教員 松島

亘志

２０２１年

３月

概要
現在，我が国には既設の橋梁が約 70 万基存在するが，2033 年にはその半数以上が建設後50 年を経過する見込みであり，老朽化が懸念されている．地方公共団体がこのような道路構
造物を定期点検するにあたって，点検業務の低コスト化が求められている．
本研究では橋梁上を走行した車両から得られる振動データから前後輪の路面凹凸を推定したのち，粒子群最適化を用いることで車両・橋梁のパラメータを同定する手法を提案する．
本研究により，事前のキャリブレーションを一切必要とせずに VBI システムの同定すること
で，車両・橋梁点検技術の簡易化を目指す．
提案手法の適用性検証として，数値シミュレーション上で車両と橋梁のモデルを作成した
のち，得られた車両振動データから前後輪における路面凹凸を推定し，最適化問題を解くことで車両・橋梁パラメータの同定を試みた．
線形カルマンフィルタを導入したことにより路面凹凸推定のロバスト性が向上したが，車両・橋梁パラメータの同定については，設定したケースにより同定精度は異なり，同定可能なパラメータは限定的であることが分かった．最適化手法の高度化や目的関数のより良い設定が 今後の課題であると考えられる．

目次
第1章
はじめに ······························································································ 4
1.1
研究背景 ································································································· 4
1.2
既往の研究 ······························································································ 4
1.3
研究目的 ································································································· 5
第2章
シミュレーションモデルの概要 ································································ 6
2.1
車両・橋梁モデル ····················································································· 62.2
車両-橋梁相互作用 ···················································································· 92.3
数値シミュレーションの妥当性検証 ··························································· 10
第3章
提案手法 ···························································································· 11
3.1
概要 ······································································································ 11
3.2
線形カルマンフィルタによる路面プロファイル推定 ······································· 11
3.2.1
状態空間モデルの構築 ······································································· 12
3.2.2
可観測性分析 ··················································································· 14
3.2.3
線形カルマンフィルタ ······································································· 15
3.2.4
Rauch-Tung-striebel 平滑化 ······························································· 16
3.3
前後輪の路面凹凸推定 ············································································· 16
3.4パラメータ最適化 ··················································································· 173.4.1
粒子群最適化 ··················································································· 17
3.4.2
本研究への適用 ················································································ 17
第4章
提案手法の適用検証 ············································································· 18
4.1
路面プロファイル・路面凹凸の推定 ··························································· 184.1.1
観測値にノイズがないことが既知である場合 ········································· 184.1.2
観測値へのノイズが未知である場合 ····················································· 20
4.2
車両・橋梁パラメータ同定 ······································································· 24
4.2.1
観測値にノイズがないことが既知である場合 ········································· 244.2.2
観測値へのノイズが未知である場合 ····················································· 25
第5章
まとめと今後の課題 ············································································· 30
謝辞 ················································································································· 31
参考文献 ··········································································································· 32

2

図目次
図 2-1 車両モデル ............................................................................................................ 9
図 2-2 車両-橋梁相互作用 ................................................................................................ 9
図 2-3 橋梁モデル ........................................................................................................... 10
図 2-4 静的荷重を与えた際の橋梁モデルのたわみ ....................................................... 10
図 3-1 パラメータ同定フロー ........................................................................................ 11
図 4-1 前後輪の路面プロファイル推定値（ノイズなしが既知） .................................. 19
図 4-2 前後輪の路面凹凸の推定値（ノイズなしが既知） ............................................. 19
図 4-3 各ノイズを加えたときの観測データ .................................................................. 21
図 4-4 前後輪の路面プロファイル推定値（RMS 値 0%ノイズ） ................................. 22
図 4-5 前後輪の路面凹凸の推定値（RMS 値 0%ノイズ） ............................................ 22
図 4-6 前後輪の路面プロファイル推定値（RMS 値 15%ノイズ） ............................... 23
図 4-7 前後輪の路面凹凸の推定値（RMS 値 15%ノイズ） .......................................... 23
図 4-8 前後輪の路面プロファイル推定値（RMS 値 35%ノイズ） ............................... 23
図 4-9 前後輪の路面凹凸の推定値（RMS 値 35%ノイズ） ...................... 23
図 4-10 前後輪の路面プロファイル推定値（既往の手法[5]・RMS 値 35%ノイズ） ... 24
図 4-11 前後輪の路面凹凸の推定値（既往の手法[5]・RMS 値 35%ノイズ） .............. 24
図 4-12 推定パラメータでの前後輪の路面凹凸推定値（ノイズなしが既知） .............. 25
図 4-13 推定パラメータでの前後輪の路面凹凸推定値（RMS 値 0%ノイズ） .............. 27
図 4-14 推定パラメータでの前後輪の路面凹凸推定値（RMS 値 15%ノイズ） ............ 28
図 4-15 推定パラメータでの前後輪の路面凹凸推定値（RMS 値 35%ノイズ） ............ 293

第1章 はじめに
1.1

研究背景

現在，我が国には既設の橋梁が約 70 万基存在するが，2033 年には約 63%[1]が建設後 50年を経過する見込みであり，老朽化が懸念されている．そこで，5 年に 1 度の近接目視点検が義務付けられている．しかし，熟練技術者や，費用を必要とする．橋梁を管理する各自治体において，財政・人的リソースの不足から補修を必要とする橋梁に満足な点検を行うことができなくなる危険性がある．よって，この問題を解決するため，低コストに実施できる橋
梁点検技術の開発が喫緊の課題である．

1.2

既往の研究

橋梁のような構造物にセンサを設置して得られる振動データなどから，構造物の健全性を評価する構造ヘルスモニタリングが注目を集めている．橋梁においては，損傷の影響が振動特性の変化として現れることから，振動特性を推定し，損傷を検知する手法が有用である．
Vehicle-Bridge Interaction（VBI：車両-橋梁相互作用）モデルを用いて，得られた振動デー
タから橋梁特性を把握する研究が盛んである．しかし，橋梁・車両振動には車両の種類や車両速度，路面凹凸など様々な要因が関与するので，そこから車両や橋梁の特性を直接，把握するのは困難である．また，橋梁振動を計測する際には各橋梁にセンサをそれぞれ直接，設置しなければならず，多大な費用と労力が見込まれる．一方で，車両振動のみからパラメー
タが推定できれば，車両にのみセンサを設置すればよいため，低コスト化と省力化が期待できる．そこで，既往の研究では，車両振動から各特性を推定するための様々な手法が提案され て い る ． Li ら [2] は 観 測 し た 車 両 加 速 度 振 動 デ ー タ か ら GSO(Glowworm Swarm
optimization)を用いて車両のパラメータを推定している．Nagayama ら[3]は車両振動に粒子フィルタを用いて，車両のパラメータを推定している．これらのように，観測データを用いて，最適化問題として構造物のパラメータを推定する研究が行われてきた．しかし，容易には計測できないパラメータが既知と仮定している場合や，そのため事前に車両や橋梁のキャ
リブレーションが必要であることなどが課題として挙げられる．また，山本ら[4]は車両の前後軸振動から空間補正により 2 点の離散点における橋梁振動モード形状を推定し，SSMA という指標を用いて橋梁の損傷を推定している．SSMA は橋梁の損傷の有無を感度良く評価することはできるが， SSMA が反応するメカニズムは力学的モデルで十分に説明できず，統計的妥当性を確認したに過ぎない．そのため，力学的モデルに従う推定手法を構築する必要がある．
そこで，著者らはパラメータ推定手法として，粒子群最適化に基づく複数車両の振動デー
タを用いた車両・橋梁・路面パラメータの推定手法を提案した[5]．本手法の特徴は，1）計測が必要な振動データは車両振動のみで橋梁振動は用いないこと，2）車両パラメータのキャリ
ブレーションを事前に行う必要がないこと 3）車両・橋梁パラメータだけでなく路面凹凸も同時に推定すること，の 3 点が挙げられる．推定された橋梁パラメータから橋梁損傷を検知し、点検に活用できるだけでなく、推定された車両パラメータから車両故障の検出・事前予測や路面凹凸推定値を用いた舗装の状態評価も期待できる．
しかし，本手法の問題の 1 つとして，ロバスト性が低いことが挙げられる．推定する過程の4

中で数値積分を行っており，観測データにノイズが載ることでその積分結果にトレンドが生じてしまうことが原因として考えられる．
Nagayama ら[6][7][8]はスマートフォンで計測した車両振動から路面縦断形状を推定する手法を提案している．ノイズを考慮した状態空間モデルに従うことで，ロバスト性のある路
面縦断形状推定が可能となっている.実際に，前後輪それぞれの推定した路面縦断形状から車両のパラメータを推定している[6]が，舗装上を走行しており，橋梁での走行を対象とはしていない．固有振動数ベースの評価手法として橋梁成分に着目した研究も行われている[7]が，車両にハンプキャリブレーションを行う必要がある．

1.3

研究目的

したがって，Nagayama ら[6]の手法を我々が提案した手法[5]と一部組み合わせることによって，提案手法のロバスト性を向上させることができるのではないかと考えた．そのため，
数値シミュレーションによって本研究の提案手法の可能性を検証することを目的とする．本
研究によって，車両や橋梁のキャリブレーションを一切必要とせずに車両・橋梁・路面を推定することができれば，より低コストのインフラ点検技術として活用されることが期待できる．

5

第2章 シミュレーションモデルの概要
2.1

車両・橋梁モデル

本研究では，車両は図 2-1 に示すハーフカーモデルを採用する．このモデルの運動方程式は次式となる．1 行目は車両のバネ上における鉛直方向，2 行目は回転方向，3，4 行目はそれぞれ前後輪のバネ下における鉛直方向の運動方程式である．
(1)
𝐌 𝒁̈(𝑡) + 𝐂 𝒁̇(𝑡) + 𝐊 𝒁(𝑡) = 𝑷𝒖(𝑡)
v

v

v

𝑚𝐻
𝐼𝑧

𝐌v =

(2)

𝑚𝑓
[

𝑐𝑓 + 𝑐𝑟

𝐿𝑟 𝑐𝑟 − 𝐿𝑓 𝑐𝑓

−𝑐𝑓

−𝑐𝑟

𝐿𝑟 𝑐𝑟 − 𝐿𝑓 𝑐𝑓
−𝑐𝑓
[
−𝑐𝑟

𝐿𝑓2 𝑐𝑓 + 𝐿2𝑟 𝑐𝑟
𝐿𝑓 𝑐𝑓−𝐿𝑟 𝑐𝑟𝐿𝑓 𝑐𝑓
𝑐𝑓0−𝐿𝑟 𝑐𝑟0
𝑐𝑟 ]𝑘𝑓 + 𝑘𝑟

𝐿𝑟 𝑘𝑟 − 𝐿𝑓 𝑘𝑓

−𝑘𝑓

−𝑘𝑟

𝐿𝑟 𝑘𝑟 − 𝐿𝑓 𝑘𝑓
−𝑘𝑓
[
−𝑘𝑟

𝐿𝑓2 𝑘𝑓 + 𝐿2𝑟 𝑘𝑟

𝐿𝑓 𝑘𝑓
𝑘𝑓 + 𝑘𝑡𝑓
0−𝐿𝑟 𝑘𝑟0
𝑘𝑟 + 𝑘𝑡𝑟 ]

𝐂v =

𝐊v =

𝑚𝑟 ]

𝐿𝑓 𝑘𝑓
−𝐿𝑟 𝑘𝑟(3)(4)

𝒁(𝑡) = [𝑧 𝜃 𝑧𝑓 𝑧𝑟 ](5)

0 0 𝑘𝑡𝑓 0 𝑇]
𝑷=[
0 0 0 𝑘𝑡𝑟
𝒖(𝑡) = [𝑢𝑓 𝑢𝑟 ]𝑇

(6)

𝑇

(7)

ここで，𝐌v ,𝐂v ,𝐊 v はそれぞれ車両の質量,減衰，剛性マトリクスである．𝑷𝒖(𝑡)は車両への強制入力であり，𝒁(𝑡)は車両応答である．
橋梁は曲げのみを考慮した 1 次元オイラー・ベルヌーイ梁とする．橋梁の運動方程式は次
式となる．
𝑛

𝜕2
𝜕 2𝑦
𝜌𝐴𝑦̈ (𝑥, 𝑡) + 2 (𝐸𝐼 2 ) = ∑ 𝛿 (𝑥 − 𝑥𝑖 )𝑃𝑖
𝜕𝑥
𝜕𝑥

(8)

𝑖=1

ここで，𝑦(𝑥, 𝑡)は橋梁のたわみであり，𝑥 は位置であり，0から橋梁のスパン長𝐿までとなる．
( ̈ )は時間𝑡における 2 次導関数である．添字𝑖は車両の各車軸であり，𝑛は車軸の数である．
荷重𝑃𝑖 は𝑖番目の車軸による接触力，𝛿 (𝑥 )はディラックのデルタ関数であり，次式となる．
0 (at 𝑥 ≠ 0)(9)
𝛿 (𝑥 ) = {∞ (at 𝑥 = 0)∞

∫ 𝛿 (𝑥 )𝑓 (𝑥 ) d𝑥 = 𝑓 (0)

(10)

−∞

ここで，𝑓 (𝑥 )は任意の関数，(8)式の𝛿 (𝑥 − 𝑥𝑖 )𝑃𝑖 は𝑥 = 𝑥𝑖 での集中荷重として表している．
6

この運動方程式を解くために，重み付き残差法によって有限要素式を導出する．この時，基底関数としてエルミート補間関数を採用する．
1
𝜙1 ( 𝑋 ) = ( 𝑋 − 1 ) 2 ( 𝑋 + 2 )
41
𝜙2 ( 𝑋 ) = ( 𝑋 − 1 ) 2 ( 𝑋 + 1 )
4(11)
1
𝜙3 ( 𝑋 ) = − ( 𝑋 + 1 ) 2 ( 𝑋 − 2 )41
𝜙4 ( 𝑋 ) = ( 𝑋 + 1 ) 2 ( 𝑋 − 1 )4
ここで，𝑋は各要素の正規化したローカル座標である．𝑗番目と(𝑗 + 1)番目のノードが𝑗番目のビーム要素を構成すると仮定すると，𝑋 = −1とは𝑗番目のノードの位置を示し，𝑋 = 1は(𝑗 +1)番目のノードの位置を示している．すると，運動方程式の解である𝑦(𝑥, 𝑡)は次式となる．
𝑦(𝑥, 𝑡) = 𝒚(𝑡) ⋅ 𝑵(𝑥 )𝑦1 (𝑡)
𝑁1 (𝑥 )
𝑁2 (𝑥 )
𝜃1 (𝑡)⋮⋮=⋅
𝑁2𝑛−1 (𝑥 )𝑦𝑛 (𝑡){𝜃𝑛 (𝑡)} { 𝑁2𝑛 (𝑥 ) }(12)

ここで，𝒚(𝑡)は解ベクトル，𝑵(𝑥 )は基底ベクトルである．𝑦𝑗 (𝑡)は𝑦(𝑥𝑗 , 𝑡)を示しており，ノード𝑥𝑗 における橋梁のたわみである．同様に，𝜃𝑗 (𝑡)はノード𝑥𝑗 における橋梁のたわみ角である．
𝑥𝑗 ≤ 𝑥 ≤ 𝑥𝑗+1 の時，基底ベクトルは次式となる．
𝜙𝑘 ( 𝑋 )𝛥𝑥
𝑁2𝑗−2+𝑘 (𝑥 ) = { 𝜙𝑘 (𝑋)20(𝑘 = 1,3)(𝑘 = 2,4)(13)(𝑘 = others)

ここで，グローバル座標𝑥は次式を用いてローカル座標𝑋に変換できる．
𝑥𝑗
𝜙1 ( 𝑋 )𝛥𝑥/2
𝜙2 ( 𝑋 )𝑥 ={𝑥 }⋅
𝜙3 ( 𝑋 )𝑗+1
𝛥𝑥/2 {𝜙4 (𝑋)}(14)

ここで，𝛥𝑥 = 𝑥𝑗+1 − 𝑥𝑗 である．残差式は次式となる．
𝑛

𝜕2
𝜕2𝑦
(
)
(
)
𝑅 𝑥, 𝑡 = 𝜌𝐴𝑦̈ 𝑥, 𝑡 + 2 (𝐸𝐼 2 ) − ∑ 𝛿 (𝑥 − 𝑥𝑖 )𝑃𝑖
𝜕𝑥
𝜕𝑥

(15)

𝑖=1

重み付き残差𝐼の平均が任意の重み関数で常にゼロに等しいと仮定すると，次式で表せる．

7

𝐿∫ 𝑤 (𝑥 )𝑅(𝑥, 𝑡)d𝑥
0
𝐿

𝑛

𝐿

𝐿
𝜕2
𝜕2𝑦
= ∫ 𝑤𝜌𝐴𝑦̈ d𝑥 + ∫ 𝑤 2 (𝐸𝐼 2 ) d𝑥 − ∑ ∫ 𝑤 (𝑥 )𝛿 (𝑥 − 𝑥𝑖 )𝑃𝑖 d𝑥
𝜕𝑥
𝜕𝑥
0
0
0

(16)

𝑖=1

𝐿

𝑛

𝐿

𝜕2𝑤 𝜕2𝑦= ∫ 𝑤𝜌𝐴𝑦̈ d𝑥 + ∫ 𝐸𝐼 2 2 d𝑥 − ∑ 𝑤 (𝑥 − 𝑥𝑖 )𝑃𝑖 = 0.
𝜕𝑥 𝜕𝑥
0
0
𝑖=1

したがって，近似解𝑦 = 𝒚 ⋅ 𝑵および𝑤 = 𝒘 ⋅ 𝑵を代入すると，
𝐿𝑛𝐿

𝜕 2 𝑵 𝜕 2 𝑵T
𝐼 = 𝒘 ([∫ 𝜌𝐴𝑵𝑵 d𝑥 ] 𝒚̈ + [∫ 𝐸𝐼 2
d𝑥 ] 𝒚 − ∑ 𝑵(𝑥 − 𝑥𝑖 )𝑃𝑖 ) = 0
𝜕𝑥 𝜕𝑥 2
0
0
T

T

(17)

𝑖=1

𝐿

𝐿

ここで，𝐌𝑏 = [∫0 𝜌𝐴𝑵𝑵T d𝑥 ], 𝐊 𝑏 = [∫0 𝐸𝐼

𝜕2 𝑵 𝜕2 𝑵T
𝜕𝑥 2 𝜕𝑥 2d𝑥 ], 𝑳 = [𝑵(𝑥 − 𝑥1 ) … 𝑵(𝑥 − 𝑥𝑛 )]，そして𝑃1
𝑷 = { ⋮ }とおくと，有限要素式として次式のように書き直せる．
𝑃𝑛𝐌b 𝒀̈(𝑡) + 𝐊 b 𝒀(𝑡) = 𝑳𝑷(18)

数値安定性のため，レイリー減衰を導入すると最終的に次式となる．
𝐌b 𝒀̈(𝑡) + 𝐂b 𝒀̇(𝑡) + 𝐊 b 𝒀(𝑡) = 𝑭(𝑡)(19)

ここで，𝐌b ,𝐂b ,𝐊 b はそれぞれ橋梁の全体質量，減衰，剛性マトリクスである．また，レイリ
ー減衰を𝐂b = 𝛼𝑏 𝐌b + 𝛽𝑏 𝐊 b とした．
車両と橋梁の振動はそれぞれの運動方程式(1),(19)から Newmark-β 法を用いて，以下のように数値積分して求めている．
𝒖̈ 𝑘 = (𝐌 +

(∆𝑡)2
∆𝑡
𝐂+
𝐊)
2
4

−1

(𝑷𝑘 + 𝐂𝑩𝑐 + 𝐊𝑩𝑘 )(20)

𝒖̈ 𝑘−1 + 𝒖̈ 𝑘
(21)
∆𝑡
2
𝒖̈ 𝑘−1 + 𝒖̈ 𝑘
(22)
(∆𝑡)2
𝒖𝑘 = 𝒖𝑘−1 + 𝒖̇ 𝑘−1 ∆𝑡 +4
ここで，𝒖𝑘 は橋梁または車両の変位応答であり，∆𝑡は時間刻み，添字𝑘はタイムステップである．𝐌，𝐂，𝐊はそれぞれ運動方程式の質量,減衰，剛性マトリクスであり，車両橋梁それぞれにおける𝐌v ,𝐂v ,𝐊 v ，𝐌b ,𝐂b ,𝐊 b に対応する．𝑩𝑐 ，𝑩𝑘 は次式となる．
𝒖̇ 𝑘 = 𝒖̇ 𝑘−1 +𝒖̈ 𝑘−1(23)∆𝑡
2
𝒖̈ 𝑘−1(24)(∆𝑡)2
𝑩𝑘 = −𝒖𝑘−1 − 𝒖̇ 𝑘−1 ∆𝑡 −4
本数値実験で扱う車両・橋梁のパラメータの値を表 2-1，表 2-2 に示す．
車速は 10 [m⁄s] で一定であり，路面上を 6 秒間走行する．車両前輪位置のスタート地点
を 0[m]地点とすると，橋梁は 10[m]~40[m]地点である．
𝑩𝑐 = −𝒖̇ 𝑘−1 −

8

𝐿𝑟

𝜃

𝑚𝐻 , 𝐼𝑧

SensorCenter of Gravity

𝑐𝑟

𝑘𝑟

𝑑
𝐿𝑓

𝑧

𝑐𝑓

𝑘𝑓

𝑧𝑟

𝑚𝑓

𝑚𝑟
𝑢𝑟

𝑘𝑡𝑓

𝑘𝑡𝑟

𝑧𝑓

𝑢𝑓

𝐸𝐼, 𝜌𝐴
図 2-1 車両モデル

2.2

表 2-1 車両のパラメータ𝑚𝐻 [kg]4.155 × 103
𝑚𝑓 [kg]4.69 × 102
𝑚𝑟 [kg]3.755 × 102
𝑐𝑓 [kg⁄s]
1.21 × 104
𝑐𝑟 [kg⁄s]
1.45 × 104
𝑘𝑓 [N⁄m]
2.28 × 105
𝑘𝑟 [N⁄m]
2.05 × 105
𝑘𝑡𝑓 [N⁄m]
2.395 × 106
𝑘𝑡𝑟 [N⁄m]
2.155 × 106
𝐼𝑧 [kg ∙ m2 ]9.011 × 103
𝐿𝑓 [m]
1.215
𝐿𝑟 [m]
1.785
𝑑[m]
0.6
表 2-2 橋梁パラメータ
𝐸𝐼[Nm3 ]𝜌𝐴[kg⁄m]
𝛼𝑏
𝛽𝑏
𝐿[m]

1.56 × 1010
4.4 × 103
0.7024
5.2 × 10−330

車両-橋梁相互作用

車両と橋梁はそれぞれが線形な力学システムとしてモデル化できるが，相互作用システムとした時は，車両移動の影響により非線形システムとなる．この相互作用を考慮した数値シミュレーションを実現するために，反復法を車両-橋梁相互作用に導入した Green らの手法[9]を本研究では用いている．その過程を図 2-2 に示す．はじめに，車両振動を路面凹凸のみを用いて計算する．そして，計算した車両振動と路面プロファイルから橋梁振動を計算する．
この橋梁振動を路面凹凸に足し合わせて路面プロファイルとし，再度同じ工程を繰り返す．
この反復過程は 1 サイクル前との算出結果の差がある規定値より小さくなった時に終了する．

Road Profile

VehicleContact Force

図 2-2 車両-橋梁相互作用9Bridge2.3

数値シミュレーションの妥当性検証

作成したシミュレーションモデルの妥当性を確認するために，表 2-2 を入力諸元とする単純梁について数値計算を行った．FEM の要素分割数は 15 としている．図 2-3 のように単純
梁の中央部 15[m]地点に 10[kN]の静的荷重を加えた．この時の変位の結果を図 2-4 に示す．
数値実験による梁の平衡状態時の変位は3.53 × 10−3 [m]であった．解析解が3.53 × 10−3 [m]なので一致しており，数値実験の精度に問題はないと分かる．この時，中央部の解析解を求める数式は次式となる．
𝑃𝐿3
(25)
48𝐸𝐼
ここで，𝑤max は中央部のたわみ，𝑃は中央部に加えた静的荷重，𝐿は梁の長さ， 𝐸𝐼は橋梁の曲げ剛性である．
𝑤max =

10[kN]

15[m]

30[m]
図 2-3 橋梁モデル

図 2-4 静的荷重を与えた際の橋梁モデルのたわみ
10

第3章 提案手法
3.1

概要

本研究では，車上計測により得られた車両振動および位置データを用いて，車両と橋梁の力学的パラメータを同定する手法を提案する．この同定プロセスでは，計測されていない橋
梁振動および路面凹凸も推定される．同定する力学的パラメータとは，車両の質量・減衰・剛性・水平方向重心位置であり，また，橋梁の質量・減衰・曲げ剛性である．先ず，同定フローの概略図を図 3-1 に示す．数値実験により，橋梁上に車両を走行させる．その際に，車両に設置した計測装置から車両振動データを得る．振動データと推定した車両パラメータから
Nagayama ら[6][7][8]が提案した手法を用いて前後輪それぞれの路面プロファイルを推定する(3.2)．その後，我々が提案した手法を用いる[5]．推定したプロファイルや橋梁パラメータを用いて橋梁振動を数値的に再現し，上記の路面プロファイルから引くことで，前後輪それぞれの推定路面凹凸を得る(3.3)．同一の橋梁を走行しているので，与えたパラメータが正解値
ならば，推定した 2 つの路面凹凸は一致するはずである．したがって，推定した路面凹凸が一致するように車両・橋梁パラメータを最適化することで，パラメータの正解値を同定することができる(3.4)と考えた．

Vehicle vibration
Kalman-Filter

Vehicle Bridge parameter

PSO

Input profile
Analytic solved
evaluation

Road surface

図 3-1 パラメータ同定フロー

3.2

線形カルマンフィルタによる路面プロファイル推定

著者が以前提案した手法[5]では，車両のバネ下加速度振動データから車両への鉛直変位を推定する過程において，観測加速度データを数値積分することで算出した速度・変位を用いて計算を行っていた．あくまで数値実験においてノイズが一切含まれない場合を考えるには影響がなかった．しかしノイズを付加した場合，積分値にはノイズによるトレンドが発生し，車両への鉛直変位を推定することができない．そこで，Nagayama ら[6][7][8]が提案した手法
に従い，一定のノイズを考慮した状態空間モデルを構築し推定を行うことでロバスト性を向上させることができると考えた．以下，提案されている手法の概略を示す．
11

3.2.1

状態空間モデルの構築

初めに，式(1)で示した運動方程式について，連続時間状態空間モデルは次式となる．
(26)
𝒁̇ (𝑡) = 𝑨 𝒁 (𝑡) + 𝑩 𝒖(𝑡) + 𝒘(𝑡)
𝑐

𝑐 𝑐

𝑐

𝒔(𝑡) = 𝑪𝑐 𝒁𝑐 (𝑡) + 𝑫𝑐 𝒖(𝑡) + 𝒗(𝑡)(27)

𝒁𝑐 (𝑡) = [𝒛(𝑡)𝑇𝒛̇ (𝑡)𝑇 ]𝑇

(28)

𝑶4×4
𝑨𝑐 = [
−𝑴−1
𝑣 𝑲

𝑰4×4
]
−𝑴−1
𝑣 𝑪

(29)

𝑩𝑐 = [

𝑶4×2
]
𝑴−1
𝑣 𝑃

(30)

ここで，𝒁𝑐 (𝑡)は車両応答に関する状態ベクトル，𝒔(𝑡)は観測ベクトルである．状態ベクト
ルとは，通常計測することができないが，力学的モデルの状態を説明するために必要である
状態変数を成分としたベクトルである．また，観測ベクトルとは力学的モデルにおいて計測することができる変数を成分としたベクトルであり，その成分はどういう値を計測するかという計測環境によって決められる．𝑪𝑐 , 𝑫𝑐 は観測ベクトル𝒔(𝑡)の成分によって定められ，具体的な成分は後述する．𝒘(𝑡), 𝒗(𝑡)はそれぞれシステムノイズ，観測ノイズである．
次に，路面プロファイルを状態変数に加えた拡張状態ベクトルを考えると，その連続時間
状態空間モデルは次式となる．
(31)
𝒁𝑎 (𝑡) = [𝒁𝑇 𝒖𝑇 ]𝑇𝑐𝒁̇𝑎 (𝑡) = 𝑨𝑐𝑎 𝒁𝑎 (𝑡) + 𝜻(𝑡)(32)

𝒔(𝑡) = 𝑪𝑐𝑎 𝒁𝑎 (𝑡) + 𝒗(𝑡)

(33)

𝑨𝑐
𝑶2×8

𝑩𝑐
]
𝑶2×2

(34)

𝑫𝑐 ]

(35)

𝒘(𝑡)
]
𝜼(𝑡)

(36)

𝑨𝑐𝑎 = [

𝑪𝑐𝑎 = [𝑪𝑐𝜻(𝑡) = [ここで，路面プロファイルに加えて，その微分も状態変数に加えることとする．プロファイ
ルの微分は高周波成分が大きくなるため，両者のプロセスノイズを考慮することにより，プロファイルの周波数特性を柔軟に扱うことができると考えられているためである[6]．すると，連続時間状態空間モデルは次式となる．
𝒁𝑎 (𝑡) = [𝒁𝑇𝑐

𝑢𝑓

𝑢𝑓̇

𝑢𝑟

𝑢𝑟̇ ]

𝑇

(37)

𝒁̇𝑎 (𝑡) = 𝑨𝑐𝑎 𝒁𝑎 (𝑡) + 𝜻(𝑡)(38)

𝒔(𝑡) = 𝑪𝑐𝑎 𝒁𝑎 (𝑡) + 𝒗(𝑡)

(39)

𝑨𝑐
𝑶4×8

𝑩𝑐 ′
]
𝑻4×4

(40)

𝑫𝑐 ]

(41)

𝑨𝑐𝑎 = [

𝑪𝑐𝑎 = [𝑪𝑐12

𝑩′𝑐 = [

𝑶4×4
]
𝑷′

(42)

0
0
𝑻=[
0
0

0
0
0
0

1
0
0
0

0
0
𝑷′ = [𝑘 ⁄𝑚
𝑢𝑓
𝑓
0

0
0
0
0

0
0
0
𝑘𝑢𝑟 ⁄𝑚𝑟𝜻(𝑡) = [

0
0
]
1
0

(43)

0
0
0]
0

(44)

𝒘(𝑡)
]
𝜼(𝑡)

(45)

この時，状態方程式の路面プロファイルの微分成分は，ガウス白色ノイズ𝜼(𝑡)を用いて，次式となる．
𝒚̈ (𝑡) = 𝜼(𝑡)(46)
次に，観測ベクトルを設定する．ここでは，センサにより車両ボディの鉛直加速度及びピッチング角速度を計測する．また，可観測性を満たすためにそれぞれを 2 回・1 回数値積分して得られる鉛直変位とピッチング角を疑似的に観測値として加えると次式のように設定でき
る．
𝒔(𝑡) = [𝑧̈ − (𝐿𝑓 − 𝑑)𝜃̈𝜃̇𝑧 − (𝐿𝑓 − 𝑑)𝜃

𝑇

𝜃]

(47)

なお，鉛直変位とピッチング角にはカットオフ周波数が 0.2[Hz]のハイパスフィルタを適用する[7]．これは観測データの積分時にノイズが低周波のトレンドとして載ってしまうことから，その影響を軽減するためである．この時，観測行列は次式となる．
𝑪 𝑇
𝑪𝑐𝑎 = [ 1 ]
𝑶4×4

13

(48)

𝑘𝑓 + 𝑘𝑟 (𝐿𝑓 − 𝑑)(𝐿𝑟 𝑘𝑟 − 𝐿𝑓 𝑘𝑓 )
+
0
1
𝑚𝐻
𝐼𝑧
𝐿𝑟 𝑘𝑟 − 𝐿𝑓 𝑘𝑓 (𝐿𝑓 − 𝑑)(𝐿𝑓2 𝑘𝑓 + 𝐿2𝑟 𝑘𝑟 )−+0 𝑑 − 𝐿𝑓𝑚𝐻
𝐼𝑧
𝑘𝑓 (𝐿𝑓 − 𝑑)𝐿𝑓 𝑘𝑓
+
0
0
𝑚𝐻
𝐼𝑧
𝑘𝑟 (𝐿𝑓 − 𝑑)𝐿𝑟 𝑘𝑟
−
0
0
𝑚𝐻
𝐼𝑧
𝑪1 =
𝑐𝑓 + 𝑐𝑟 (𝐿𝑓 − 𝑑)(𝐿𝑟 𝑐𝑟 − 𝐿𝑓 𝑐𝑓 )
−
+
0
0
𝑚𝐻
𝐼𝑧
𝐿𝑟 𝑐𝑟 − 𝐿𝑓 𝑐𝑓 (𝐿𝑓 − 𝑑)(𝐿𝑓2 𝑐𝑓 + 𝐿2𝑟 𝑐𝑟 )
−
+
1
0
𝑚𝐻
𝐼𝑧
𝑐𝑓
(𝐿𝑓 − 𝑑)𝐿𝑓 𝑐𝑓
+
0
0
𝑚𝐻
𝐼𝑧
(𝐿𝑓 − 𝑑)𝐿𝑟 𝑐𝑟
𝑐𝑟
−
0
0
[
𝑚𝐻
𝐼𝑧
−

0
1
0
0

(49)

0
0
0
0

]

以上より作成した連続時間状態空間モデルを離散化し，離散時間状態空間モデルへと変換
すると，次式となる．
(50)
𝒁𝑎 = 𝑨𝑎 𝒁𝑎 + 𝜻𝑘𝑘+1𝑘𝒔𝑘 = 𝑪𝑎 𝒁𝑎𝑘 + 𝒗𝑘

(51)

𝒘𝑘
𝜻(𝑡) = [ 𝜼 ](52)

𝑨𝑎 = 𝒆𝑨𝑐𝑎 𝑑𝑡(53)

𝑪𝑎 = 𝑪𝑐𝑎(54)

𝒁𝑎𝑘 = 𝒁𝑎 (𝑘𝑑𝑡)(55)

𝒔𝑘 = 𝒔(𝑘𝑑𝑡)

(56)

𝑘

ここで，

である．ここでは，サンプリング時間間隔を𝑑𝑡として離散化し，タイムステップを𝑘としている．𝒘𝒌 , 𝒗𝒌 , 𝜼𝒌 に関する分散共分散行列は次のように定義する．
(57)
𝐸[𝒘𝑘 𝒘𝑇 ] = 𝑸𝑤,𝑘
𝑘

3.2.2

𝐸[𝒗𝑘 𝒗𝑇𝑘 ] = 𝑹𝑘

(58)

𝐸[𝜼𝑘 𝜼𝑇𝑘 ] = 𝑸𝜂,𝑘(59)

可観測性分析

出力を観測した状態空間モデルから内部の状態を推定可能であるかの尺度として可観測性分析を行う必要がある．構築した状態空間モデルの可観測性を調べるための可観測行列は次
式となる．
14

𝑶 = [𝑪𝑇𝑎2

𝑨𝑇𝑎 𝑪𝑇𝑎𝑨𝑇𝑎 𝑪𝑇𝑎…

𝑨𝑇𝑎

𝑛−1 𝑇 𝑇𝑪𝑎 ](60)

ここで，𝑛は状態ベクトルの自由度であり，本研究では𝑛 = 12である．
システムが可観測であるための必要十分条件として，可観測行列が次式を満たさなければ
ならない．
(61)
rank(𝑶) = 𝑛

3.2.3

線形カルマンフィルタ

以上より定義した状態空間モデルに線形カルマンフィルタを適用する．線形カルマンフィ
ルタとは，正規線形状態空間モデルの状態変数を推定するための代表的な手法の一つであり，状態方程式および観測方程式の線形性，ノイズの正規性を前提とする．本研究において，状態空間モデルは式(38)(39)より線形であり，ノイズには後述する白色ガウス雑音を疑似的に用いるので，この仮定は満たされると考えて良い．
̂𝑎 と事後誤差共分散行列𝑷
̂𝑎 と𝑷
̂𝑘 の初期値である𝑿
̂0 を次のように初めに，事後状態推定値𝑿
0
𝑘
表す．
(62)
̂𝒂 = 𝐸[𝑿𝑎 ]𝑿
𝟎0̂0 = 𝐸[(𝑿𝑎0 − 𝐸 [𝑿𝑎0 ])(𝑿𝑎0 − 𝐸 [𝑿𝑎0 ])𝑇 ]
𝑷
𝑎 −
̂
̂−
事前状態推定値𝑿
𝑘+1 と事前誤差分散共分散行列𝑷𝑘+1 は次式となる．
𝑎 −
̂
̂𝑎
𝑿
= 𝑨𝑎 𝑿
𝑘+1
−

𝑘

̂
̂ 𝑇
𝑷
𝑘+1 = 𝑨𝑎 𝑷𝑘 𝑨𝑎 + 𝑸𝑘

(63)
(64)
(65)

ここで，𝑸𝑘 はプロセスノイズ分散共分散行列である．
𝑎̂
次に，𝑘 + 1ステップにおけるカルマンゲイン𝑮𝑘+1 ，事後状態推定値𝑿
𝑘+1 ，そして事後誤差
̂
分散共分散行列𝑷は次式となる．
𝑘+1
−1
− 𝑇
̂
̂− 𝑇
𝑮𝑘+1 = 𝑷
𝑘+1 𝑪𝑎 (𝑪𝑎 𝑷𝑘+1 𝑪𝑎 + 𝑹𝑘+1 )

(66)

𝑎
𝑎 −
𝑎 −
̂
̂
̂
𝑿
𝑘+1 = 𝑿𝑘+1 + 𝑮𝑘+1 (𝒖𝑘+1 − 𝑪𝑎 𝑿𝑘+1 )

(67)

̂
̂−
𝑷
𝑘+1 = (𝑰 − 𝑮𝑘+1 𝑪𝑎 )𝑷𝑘+1(68)

ここで，𝑹𝑘+1 は観測ノイズ分散共分散行列である．
̂𝑎 と，
カルマンゲイン𝑮𝑘+1 は，タイムステップ𝑘の状態方程式から導かれる状態推定値𝑨𝑎 𝑿
𝑘
タイムステップ𝑘 + 1の観測方程式から導かれる状態推定値𝑪𝑎−1 𝒖𝑘+1 のどちらの推定値を信頼するかの塩梅を決める役割をしている．もし，状態方程式による状態推定値が信用できると−̂
して事前誤差分散共分散行列𝑷
𝑘+1 をゼロとすると，カルマンゲイン𝑮𝑘+1 はゼロとなり，事後𝑎̂
̂𝑎
状態推定値𝑿
𝑘+1 は𝑨𝑎 𝑿𝑘 となるから，これは観測データを一切用いず力学的モデルに従うこととなる．一方，観測方程式による状態推定値が信頼できるとして観測ノイズ分散共分散行列𝑎−1
̂
𝑹𝑘+1 をゼロとすると，カルマンゲイン𝑮𝑘+1 は𝑪−1
𝑎 となり，事後状態推定値𝑿𝑘+1 は𝑪𝑎 𝒖𝑘+1 とな
るから，これは観測データのみから推定していることになる．

15

3.2.4

Rauch-Tung-striebel 平滑化

次に，Rauch-Tung-Striebel(RTS)平滑化[10]を適用する．カルマンフィルタによって𝑘 = 0~𝑁に関して事後状態推定値等が求められている時，まず𝑘 = 𝑁において推定誤差分散共分散行列
̂𝑏 と状態ベクトル𝑿𝑅𝑇𝑆 を初期化する．
𝑷𝑁𝑁
̂𝑏 = 𝑷
̂𝑁
𝑷
𝑁

(69)

𝑅𝑇𝑆
̂
̂𝑎
𝑿
=𝑿
𝑁
𝑁

(70)

𝑅𝑇𝑆
̂
以下の式を繰り返すことで，状態推定値𝑿𝑘 を求める．
− −1
̂𝑘 𝑨𝑇𝑎 [𝑷̂𝚽𝑘 = 𝑷
𝑘+1 ]

(71)

−
T
̂𝑏 = 𝑷̂𝑏̂𝑘 − 𝚽𝑘 (𝑷̂𝑷
𝑘+1 − 𝑷𝑘+1 ) 𝚽k
𝑘

(72)

𝑅𝑇𝑆
𝑅𝑇𝑆
𝑎 −
̂
̂
̂𝑎 + 𝚽𝑘 (𝑿
̂
𝑿
=𝑿
𝑘
𝑘
𝑘+1 − 𝑿𝑘+1 )

(73)

𝑅𝑇𝑆
̂
以降，本論文で言及される路面プロファイルの推定値は，この𝑿
𝑘 のことを指すこととする．

3.3

前後輪の路面凹凸推定

カルマンフィルタにより推定される状態値は以下である．
𝑇

𝒁𝒂 (𝑡) = [𝑧 𝜃 𝑧𝑓 𝑧𝑟 𝑧̇ 𝜃̇ 𝑧𝑓̇ 𝑧𝑟̇ 𝑢𝑓 𝑢𝑓̇ 𝑢𝑟 𝑢𝑟̇ ](74)

車両の前後輪から橋梁への接地力はそれぞれ，次式となるから，推定した状態値のうち
𝑧𝑓 , 𝑧𝑟 , 𝑢𝑓 , 𝑢𝑟 を用いて算出できる．
(75)
𝑃𝑓 = (𝑚𝐻𝑓 + 𝑚𝑓 )𝑔 + 𝑘𝑡𝑓 (𝑧𝑓 − 𝑢𝑓 )𝑃𝑟 = (𝑚𝐻𝑟 + 𝑚𝑟 )𝑔 + 𝑘𝑡𝑟 (𝑧𝑟 − 𝑢𝑟 )(76)

ここで，橋梁の運動方程式は次式となる．上記の式を次式に代入し解くことにより，橋梁振
動を得る．
(77)
𝐌 𝒀̈(𝑡) + 𝐂 𝒀̇(𝑡) + 𝐊 𝐘(𝑡) = 𝑷b𝑏𝑏

最後に，強制入力と橋梁振動から次式より路面凹凸を算出する．
𝒓=𝒖−𝒀
このようにして車両振動データから路面凹凸を解析的に算出した．

16

(78)

3.4

パラメータ最適化

3.4.1

粒子群最適化

粒子群最適化(Particle Swarm Optimization)[11]は，最適化問題を解く際のメタヒューリステ
ィクス手法の 1 つである．1995 年に J.Kennedy と R.C.Eberhart らによって提案された．粒子
群最適化は解空間に粒子を複数配置し，それらの粒子は他の粒子と相互に情報交換を行い，解空間を探索することで最適化問題の近似解を発見することのできる手法であり，群知能の1 つである．各粒子は速度と位置を持っており，それぞれの更新は以下の式で行われる．
𝑡+1
𝑡
𝑡
𝑡
𝑡
(79)
𝑣𝑖𝑗
= w𝑣𝑖𝑗+ 𝑐1 1 (𝑝𝑏𝑖𝑗− 𝑥𝑖𝑗) + 𝑐2 2 (𝑔𝑏𝑗𝑡 − 𝑥𝑖𝑗
)
(80)

𝑡+1
𝑡
𝑡+1
𝑥𝑖𝑗
= 𝑥𝑖𝑗
+ 𝑣𝑖𝑗

ここで、w は慣性速度定数，𝑐1 ,𝑐2 は加速度係数， 1 , 2 は[0,1]の一様乱数である．また，pbは自身の探索過程における最良の位置ベクトル，gb は粒子群の探索過程における最良の位置ベクトルである．
この時に最終的な gb を最適解の近似解とみなす．

3.4.2

本研究への適用

本研究が対象とするシステム同定問題は，ブラックボックス最適化問題の一種であり，目的関数の形状が不明である．したがって，目的関数を微分し，その勾配を用いて最適解を探索する方法は適用できない．そこで，最適化問題の近似解を求めるための代表的なメタヒュ
ーリスティクスアルゴリズムの一つである PSO を適用することとした．さらに，解空間の探索に多様性を持たせるため，
PSO に動的な近傍構造を導入した Independent-minded PSO（IPSO）[12]を採用する．IPSO では，粒子が他粒子と一定確率で情報を交換し，最適解を探索する．
具体的には，あるステップにおいて対象の粒子が，それまでの粒子群の探索過程における最良の位置ベクトル gb を速度更新に反映させるのかどうかを確率的に決める．
粒子群内での i 番目の粒子を𝑥 (𝑖) とすると，
𝒙

(𝒊)

=

{𝑚𝑓 , 𝑚𝑟 , 𝑐𝑓 , 𝑐𝑟 , 𝑘𝑓 , 𝑘𝑟 ,

(𝑖)

(81)

𝑘𝑡𝑓 , 𝑘𝑡𝑟 , 𝐿𝑓 , 𝐸𝐼, 𝜌𝐴, 𝛼𝑏 , 𝛽𝑏 }で表される．
3.2，3.3 より示した手順に従って，観測した車両振動データと𝑥𝑖 から車両の前後輪の路面
凹凸 𝑓 ， 𝑟 を算出する．この 2 つには同時刻においてホイールベース長のズレが生じるが，車両速度からズレを補正し，両輪が通過した部分の路面凹凸のみを用いることとした．
本研究の最適化問題の目的関数F(𝑥𝑖 )は次式のように定める．
2

F(𝑥𝑖 ) = ∑( 𝑓 − 𝑟 )この目的関数が最小となる値を探索する．

17

(82)

第4章 提案手法の適用検証
ハーフカーモデルとベルヌーイ・オイラー梁を 2 章で示したように構築し，提案手法を適用することで検証する．なお，パラメータ推定において車両の総重量(𝑚𝐻 + 𝑚𝑓 + 𝑚𝑟 )と車両のホイールベース長(𝐿𝑓 + 𝐿𝑟 )は簡易に計測できると判断し，それぞれ既知であるとしている．
観測データのサンプリング時間間隔は 0.01[s]とし，路面凹凸のサンプリング間隔は 10[cm]である．
本章で検証するケースを表 4-1 に示す．全ケースについて，路面プロファイル・路面凹凸を順に推定し，最終的に車両橋梁パラメータを推定するという流れに変わりはない．ノイズがないことが既知であるケースでは，ノイズを数値積分することによるトレンドが載らないことが分かっているため，観測値，そして推定したい路面プロファイル・路面凹凸にもハイパスフィルタを適用せずに推定した．一方，ノイズの大きさが未知である場合，ノイズの RMS
値が 0%（つまりノイズは付加されていない）だとしても，ハイパスフィルタを適用し，トレンドを抑える処理を行っている違いがある．実際には，観測値にノイズは当然載るため，ノイズがないことが既知であるケースというのは数値実験上での理想的な条件に過ぎず，ノイズの大きさが未知であるケースが現実的である．
表 4-1
ケース

ノイズの
RMS

ノイズがないことが既知

検証ケース
推定する路面凹凸・
路面プロファイル
真値

0%
ノイズの大きさが未知

推定パラメータ15%

真値に 0.2Hz のハイパスフィルタを適用したもの

車両・橋梁
パラメータ

35%

4.1
4.1.1

路面プロファイル・路面凹凸の推定
観測値にノイズがないことが既知である場合

車両パラメータが既知で，観測ノイズもないことが既知である理想的な状態からカルマン
フィルタによる路面プロファイル推定を考える．数値シミュレーションにて取得した観測値にノイズを付加せず，またその情報が既知であるため，取得した観測値のうち，センサ位置での鉛直加速度及びピッチング角速度を 2 回・1 回数値積分した鉛直変位及びピッチング角には 0.2Hz のハイパスフィルタを適用せずに観測値として定めた．車両のパラメータが既知であり，観測データにノイズがないことも既知であるため，𝑸𝑘 の第 1～8 成分，𝑹𝑘 のいずれの成分も 0 とするのが適切であると考えたが，数値安定性の検討を経て，Nagayama ら[6]を参18

考に適度に小さい値をそれぞれ設定した．また，𝑸𝑘 の第 9～12 成分は，拡張状態ベクトルとした路面プロファイル成分にそれぞれ対応しているため，実験対象とする路面プロファイル
やその微分値のオーダーに基づき決定した．分散共分散行列𝑸𝑘 ，𝑹𝑘 ともに非対角成分は 0 としている．設定したプロセスノイズ・観測ノイズ分散共分散行列の成分を表 4-2 に示す．
以上より設計したカルマンフィルタにより前後輪の路面プロファイル及び路面凹凸を推定した結果を図 4-1，図 4-2 に示す．前後輪ともに路面プロファイル真値とほぼ一致していることが分かる．前輪における推定プロファイルと真値の差は，全長で合計しても 0.037m であり，後輪においても 0.038mm という非常に小さい誤差しかない．ノイズがなく，観測値もそのま𝑎 −̂
ま用いているため，事前状態推定値𝑿
𝑘+1 が状態真値とほとんど誤差なく状態方程式から導
𝑎 −
𝑎
̂
̂
かれ，観測値𝒖𝑘+1 と𝑪𝑎 𝑿𝑘+1 にもほとんど差が出ないことから事後状態推定値𝑿
𝑘+1 も非常に高精度に推定されていると考えられる．また，路面凹凸の推定値も前後輪でほぼ一致していることから橋梁成分が高精度に推定され差し引かれていることが分かる．
表 4-2 ノイズ無しが既知の時の分散共分散行列
行列成分
設定値
行列成分
設定値
𝑸𝑘 (1,1)
𝑸𝑘 (9,9)5.7 × 10−7
1.0 × 10−5
𝑸𝑘 (2,2)
𝑸𝑘 (10,10)
5.8 × 10−7
1.0 × 10−5
𝑸𝑘 (3,3)
𝑸𝑘 (11,11)
1.2 × 10−7
1.0 × 10−5
𝑸𝑘 (4,4)
𝑸𝑘 (12,12)
1.2 × 10−7
1.0 × 10−5
𝑸𝑘 (5,5)
𝑹𝑘 (1,1)
7.6 × 10−7
5.2 × 10−7
𝑸𝑘 (6,6)
𝑹𝑘 (2,2)
2.4 × 10−7
5.2 × 10−7
𝑸𝑘 (7,7)
𝑹𝑘 (3,3)
1.1 × 10−7
5.7 × 10−7
𝑸𝑘 (8,8)
𝑹𝑘 (4,4)
1.1 × 10−7
5.8 × 10−7

図 4-1 前後輪の路面プロファイル推定値（ノイズなしが既知）

図 4-2 前後輪の路面凹凸の推定値（ノイズなしが既知）

19

4.1.2

観測値へのノイズが未知である場合

現実的に得られる観測データには必ず観測ノイズが載ることから，様々なノイズを付加した際の路面プロファイル推定について考える．数値実験により取得した観測値のうち，セン
サ位置での鉛直加速度及びピッチング角速度には，その観測値の RMS 値に基づいた RMS 値を持ったガウス性白色ノイズを， 観測ノイズとしてそれぞれ付加した．N 個のデータを持つ
データ𝑥に対して RMS（二乗平均平方根）値は以下の式から導かれる．
𝑁

1
RMS[𝑥 ] = √ ∑(𝑥𝑖 )2
𝑁

(83)

𝑖=1

その後，前述の通り，それぞれを 2 回・1 回数値積分した鉛直変位及びピッチング角にはカッ
トオフ周波数が 0.2Hz のハイパスフィルタを適用し，疑似的に観測値として加えた．この際，
路面プロファイルも真値に同様のフィルタを適用したものを推定対象としている．未知のノイズ成分に考慮してカルマンフィルタを適用する際には，プロセスノイズ分散共分散行列𝑸𝑘
や観測ノイズ分散共分散行列𝑹𝑘 を適切に設定することが推定精度に大きく影響を及ぼす．プ
ロセスノイズは実車両をハーフカーモデルなどの単純なモデルにモデル化した際の誤差や，走行環境など様々な要因によりオーダーが定められるため，適切に設定するのが困難である．
本ケースでは車両のパラメータを既知としているため，𝑸𝑘 の第 1~8 成分は小さい値を設定しつつも，試行錯誤的に設定した．𝑸𝑘 の第 9~12 成分は先ほどと同様，路面プロファイルやその微分値のオーダーに基づき決定した．観測ノイズ分散共分散行列𝑹𝑘 の各成分は，付加したノイズのレベルごとに試行錯誤的に定めた．
本研究では，観測データであるセンサ位置での鉛直加速度及びピッチング角速度の RMS 値に対して，0%，15%，35%の RMS 値を持つガウス性白色ノイズをそれぞれ付加したケースを検証した．その時のセンサ位置での鉛直加速度とピッチング角速度を図 4-3 に示す．各ケー
スにおいてのプロセスノイズ・観測ノイズ分散共分散行列の成分をそれぞれ表 4-3，表 4-4．
表 4-5 に示す．
以上より設計したカルマンフィルタにより前後輪の路面プロファイル及び路面凹凸を推定した結果をそれぞれ図 4-4，図 4-5，図 4-6，図 4-7，図 4-8，図 4-9 に示す． 前後輪の路面
プロファイル推定値は真値との差は走行個所によるものの，誤差が確認できる．RMS 値 0%
ノイズの場合は，ノイズが実際に付加されていないものの，前輪における真値との誤差は，全長の合計で 0.40m，後輪においては 0.69m と大きく増加している．これは観測値のうち数値積分して得られるものに 0.2Hz のハイパスフィルタを適用していることが理由と考えられる．数値積分による低周波成分のトレンドを抑えるためにハイパスフィルタは必要ではあるのだが，これによって観測値から除くべきでないデータが取れ除かれてしまっている可能性が大いにある．次に，ノイズの RMS 値が大きくなるにつれて誤差も大きくなることが分かる．RMS 値 15%ノイズの場合は，前輪における路面プロファイル真値との誤差は，全長の合計で 0.75m，後輪においては 0.88m であり，RMS 値 35%ノイズの場合は，前輪における路面プロファイル真値との誤差は，全長の合計で 1.04m，後輪においては 1.43m であった．これはノイズが大きくなれば観測データが当てにならなくなり，推定精度が落ちてしまうと考えられる．次に，ノイズの大きさに関わらず後輪よりも前輪で推定精度が良いことが分かる．これは観測値を取得するセンサが前輪軸の方に近いことから推定感度が高いことが理由として20考えられる．これらの結果は前述の通り，システムノイズ・観測ノイズ分散共分散行列𝑸𝑘 , 𝑹𝑘
をどう設定するかによって変化する．今回の結果は試行錯誤的に各成分を決定したが，この設定が最も良いとは決して言えない．この設定をより最適に行うことは推定精度を高めるう
えで非常に重要であるが，それぞれの成分が複雑に関係しており，最適な設定を探すことは困難である．この問題に関して Nagayama ら[6]は，Robbins-Monro Algorithm[13]を導入し各成分を動的に見積もっていることから，本研究への適用により改善が期待できる．次に，著者が以前提案した手法[4]を用いて，観測値である車両加速度振動に 35%の RMS 値を持つ
ガウス性白色ノイズを付加した時の路面プロファイル及び路面凹凸推定結果を図 4-10，図 411 に示す．今回の結果と比較すると，本研究の提案手法が路面プロファイル及び路面凹凸の推定に対するロバスト性が飛躍的に向上していることが見て取れる．これには，ハイパスフ
ィルタによる低周波トレンドの除去やカルマンフィルタ適用によりノイズを適切に考慮し推定できていることが理由として考えられる．

図 4-3 各ノイズを加えたときの観測データ21

表 4-3 RMS 値 0%ノイズの時の分散共分散行列
行列成分
設定値
行列成分
設定値
−7
𝑸𝑘 (1,1)
𝑸𝑘 (9,9)5.7 × 10
1.0 × 10−5
𝑸𝑘 (2,2)
𝑸𝑘 (10,10)
5.8 × 10−7
1.0 × 10−5
𝑸𝑘 (3,3)
𝑸𝑘 (11,11)
1.2 × 10−7
1.0 × 10−5
𝑸𝑘 (4,4)
𝑸𝑘 (12,12)
1.2 × 10−7
1.0 × 10−5
𝑸𝑘 (5,5)
𝑹𝑘 (1,1)
7.6 × 10−7
5.2 × 10−7
𝑸𝑘 (6,6)
𝑹𝑘 (2,2)
2.4 × 10−7
5.2 × 10−7
𝑸𝑘 (7,7)
𝑹𝑘 (3,3)
1.1 × 10−7
5.7 × 10−7
𝑸𝑘 (8,8)
𝑹𝑘 (4,4)
1.1 × 10−7
5.8 × 10−7

図 4-4 前後輪の路面プロファイル推定値（RMS 値 0%ノイズ）

図 4-5 前後輪の路面凹凸の推定値（RMS
値 0%ノイズ）表 4-4 RMS 値 15%ノイズの時の分散共分散行列
行列成分
設定値
行列成分
設定値
−7
𝑸𝑘 (1,1)
𝑸𝑘 (9,9)5.7 × 10
1.0 × 10−5
𝑸𝑘 (2,2)
𝑸𝑘 (10,10)
5.8 × 10−7
1.0 × 10−5
𝑸𝑘 (3,3)
𝑸𝑘 (11,11)
1.2 × 10−7
1.0 × 10−5
𝑸𝑘 (4,4)
𝑸𝑘 (12,12)
1.2 × 10−7
1.0 × 10−5
𝑸𝑘 (5,5)
𝑹𝑘 (1,1)
7.6 × 10−7
3.2 × 10−1
𝑸𝑘 (6,6)
𝑹𝑘 (2,2)
2.4 × 10−7
9.2 × 10−4
𝑸𝑘 (7,7)
𝑹𝑘 (3,3)
1.1 × 10−7
2.7 × 10−1
𝑸𝑘 (8,8)
𝑹𝑘 (4,4)
1.1 × 10−7
9.8 × 10−422

図 4-6 前後輪の路面プロファイル推定値（RMS 値 15%ノイズ）

図 4-7 前後輪の路面凹凸の推定値（RMS
値 15%ノイズ）

表 4-5 RMS 値 35%ノイズの時の分散共分散行列
行列成分
設定値
行列成分
設定値
−6
𝑸𝑘 (1,1)
𝑸𝑘 (9,9)5.7 × 10
1.0 × 10−5
𝑸𝑘 (2,2)
𝑸𝑘 (10,10)
5.8 × 10−6
1.0 × 10−5
𝑸𝑘 (3,3)
𝑸𝑘 (11,11)
1.2 × 10−6
1.0 × 10−5
𝑸𝑘 (4,4)
𝑸𝑘 (12,12)
1.2 × 10−6
1.0 × 10−5
𝑸𝑘 (5,5)
𝑹𝑘 (1,1)
7.6 × 10−6
7.2 × 10−1
𝑸𝑘 (6,6)
𝑹𝑘 (2,2)
5.2 × 10−4
2.4 × 10−6
𝑸𝑘 (7,7)
𝑹𝑘 (3,3)
1.1 × 10−6
7.7 × 10−1
𝑸𝑘 (8,8)
𝑹𝑘 (4,4)
5.8 × 10−4
1.1 × 10−6

図 4-8 前後輪の路面プロファイル推定値（RMS 値 35%ノイズ）

図 4-9 前後輪の路面凹凸の推定値（RMS
値 35%ノイズ）23図 4-10 前後輪の路面プロファイル推定値（既往の手法[5]・RMS 値 35%ノイズ）

4.2

図 4-11 前後輪の路面凹凸の推定値（既往の手法[5]・RMS 値 35%ノイズ）

車両・橋梁パラメータ同定

車両モデルが既知である前節で設計したカルマンフィルタを用い，本来は車両や橋梁が未知であるとして，路面凹凸及びパラメータ推定を行う．観測データから前後輪の路面凹凸を推定し，目的関数である路面凹凸の 2 乗誤差を最小化する最適化問題として粒子群最適化を適用した．粒子群最適化における粒子数は 1000，更新回数は 5 回とした．パラメータごとのの探索空間の上下限は，真値に対して 60%~140%とし，一様分布により初期位置を生成した．
パラメータ同定に使用した粒子群最適化から得られる結果は試行ごとに同様の結果を得られるとは限らないため，各ケースにおいて 10 回試行した．

4.2.1

観測値にノイズがないことが既知である場合

パラメータ推定結果を表 4-6 に示す．𝐹(𝑥)の値はすべて10−6 オーダーである．全パラメー
タが真値の時の目的関数の値は3.25 × 10−8 であることから，10 回の推定すべてで目的関数の値が比較的大きいものの，非常に小さい値を探索できていることが分かる．しかし，パラメ
ータごとに見ると，𝑚𝑓 , 𝐿𝑓 , 𝐸𝐼は比較的精度よく推定できているように考えられるものの，その他のパラメータは推定精度が比較的悪いことが分かる．これは，パラメータごとの推定感
度の差によるものと考えられる．次に，10 回の試行のうち，最も目的関数の値が小さい時のパラメータに基づく前後輪ごとの路面凹凸を推定した結果を図 4-10 に示す．前後輪の路面凹
凸の推定値はほとんど誤差がなく，路面凹凸の真値とも誤差が小さいことが分かる． 以上より，前後輪の推定路面凹凸は高精度に一致しているものの，パラメータ同定の精度が良いとは言えない結果となった．これは初期収束性によって解空間を探索しきれていないことや，目的関数の形状が複雑となっている可能性が考えられる．この問題を解決するために，最適化手法の高度化が今後の課題として挙げられる．

24

𝑚𝑓
1
2
3
4
5
6
7
8
9
10

1.06
1.05
0.94
0.90
1.04
0.90
1.11
0.99
0.93
0.93

𝑚𝑟
1.06
0.97
0.98
1.08
1.14
1.22
0.87
0.97
0.87
1.05

𝑐𝑓
1.02
1.09
1.01
0.92
1.13
0.88
1.01
1.02
1.06
1.09

𝑐𝑟
1.22
0.91
1.14
1.12
1.17
1.21
1.02
1.12
0.93
1.24

表 4-6 パラメータ推定結果
𝑘𝑓
𝑘𝑡𝑓 𝑘𝑡𝑟
𝐿𝑓
𝑘𝑟
0.97
1.12
1.02
0.93
1.05
0.89
1.08
1.00
1.23
1.00

1.13
1.04
1.16
1.12
1.28
1.11
1.08
1.14
1.20
1.16

1.03
1.24
0.88
0.92
1.00
0.88
1.23
1.06
1.14
0.84

1.07
1.22
0.96
1.12
1.05
1.18
0.96
1.04
1.17
0.99

1.08
0.94
1.06
1.12
1.04
1.17
0.97
1.05
0.95
1.08

𝐸𝐼
1.03
1.04
1.08
1.02
1.06
1.03
0.98
1.06
0.96
1.08

𝜌𝐴
1.24
1.18
0.96
0.85
0.97
0.97
1.15
0.77
0.97
0.92

𝛼𝑏
1.16
0.84
1.35
1.14
0.83
0.93
0.99
0.87
0.84
0.98

𝛽𝑏
0.73
1.11
0.76
0.84
0.85
0.89
0.66
0.87
0.98
0.96

𝐹(𝑥)
3.33
3.66
2.13
0.57
4.96
2.12
3.49
2.15
3.74
2.75

図 4-12 推定パラメータでの前後輪の路面凹凸推定値（ノイズなしが既知）

4.2.2

観測値へのノイズが未知である場合

観測データであるセンサ位置での鉛直加速度及びピッチング角速度の RMS 値に対して，0%，15%，35%の RMS 値を持つガウス性白色ノイズをそれぞれ付加したケースについてのパ
ラメータ推定結果及びその時の目的関数値𝐹(𝑥)を表 4-7，表 4-8，表 4-9 に示す．𝐹(𝑥)の値は
表 4-7 で10−5 オーダー，表 4-8，表 4-9 で10−4 オーダーである．全パラメータが真値の時の目的関数の値はそれぞれ，RMS 値 0%ノイズの時4.33 × 10−5 ，RMS 値 15%ノイズの時
3.04 × 10−4 ，RMS 値 35%ノイズの時3.45 × 10−4 である．各ケースにおいて目的関数値𝐹(𝑥)が全パラメータ真値の時の値よりも小さくなっている．よって，本研究で設定した目的関数は25

全パラメータ真値で最小値とならないことが分かる．ここで，目的関数である前後輪の路面
凹凸について考える．前後輪に入力される路面凹凸を等しいため，(84)
𝑢𝑓 (𝑥) − 𝑦𝑓 (𝑥) = 𝑢𝑟 (𝑥) − 𝑦𝑟 (𝑥)ここで，𝑥は位置，𝑢𝑓 (𝑥), 𝑢𝑟 (𝑥)は前後輪の路面プロファイル，𝑦𝑓 (𝑥), 𝑦𝑟 (𝑥)は前後輪へ入力される橋梁変位振動である．前輪の走行開始位置を-10[m]とすると，車両のホイールベース長は3[m]であることから，(85)
𝑢𝑓 (−10: 47) − 𝑦𝑓 (−10: 47) = 𝑢𝑟 (−10: 47) − 𝑦𝑟 (−10: 47)
ハイパスフィルタの線形性から，𝐻 (𝑢𝑓 (−10: 47)) − 𝐻 (𝑦𝑓 (−10: 47)) = 𝐻(𝑢𝑟 (−10: 47)) − 𝐻(𝑦𝑟 (−10: 47))(86)

ここで，𝐻 ( )はハイパスフィルタを適用したことを示す．
本研究ではまず，線形カルマンフィルタによって路面プロファイルを推定しているが，これは走行開始から終了までのハイパスフィルタが適用された路面プロファイルが推定されてるため，出力される前後輪の路面プロファイルはそれぞれ，𝐻 (𝑢𝑓 (−10: 50)), 𝐻(𝑢𝑟 (−13: 47))となっている．したがって，本提案手法によって推定される路面凹凸は，(𝐻 (𝑢𝑓 (−10: 50)) − 𝐻 (𝑦𝑓 (−10: 50))) (−10: 47)(87)≠ (𝐻(𝑢𝑟 (−13: 47)) − 𝐻(𝑦𝑟 (−13: 47))) (−10: 47)
前後輪の路面凹凸をそれぞれ 𝑓 (𝑥), 𝑟 (𝑥)とすると(𝐻 ( 𝑓 (−10: 50))) (−10: 47) ≠ (𝐻( 𝑟 (−13: 47))) (−10: 47)(88)

以上より，不等号となることから，推定している前後輪の路面凹凸に差が生じてしまっている．実際には，橋梁部分は 0~30[m]地点であるため，車両へ入力される橋梁振動は等しく，
路面プロファイルの計測距離による違いがズレの要因となる．これが本ケースにおいて全パ
ラメータ真値で目的関数が最小値とならない理由の一つである．小さなずれであるため，本研究では前後輪で路面凹凸が一致するものとして目的関数を定めたが，パラメータを同定する結果に大きく影響を及ぼしていることが考えられる．
10 回の試行のうち，最も目的関数の値が小さい時のパラメータに基づく前後輪ごとの路面凹凸を推定した結果をそれぞれ図 4-11，図 4-12，図 4-13 に示す． 観測値に付加するノイズが大きくなるにしたがって振動による誤差等が大きくなるものの，真値とそこまで大きな差はないことが確認できる．しかし，この誤差が本研究のパラメータ推定にとっては，特に感度の低いパラメータへの推定精度低下への大きくなると考えられる．これには分散共分散行
列の適切な設定などが課題として考えられるが，真値に限りなく近づけることは現実的に難しいため，パラメータ推定に関してもロバスト性を持つ目的関数の新たな構築が必要である．

26

𝑚𝑓
1
2
3
4
5
6
7
8
9
10

1.21
0.91
0.91
0.95
1.02
1.16
0.98
0.93
0.83
1.16

𝑚𝑟
1.15
0.80
0.83
1.08
1.08
0.95
1.16
0.91
0.95
0.89

表 4-7 パラメータ推定結果(RMS 値 0%ノイズ)
𝑐𝑓
𝑘𝑓
𝑘𝑡𝑓 𝑘𝑡𝑟
𝐿𝑓
𝑐𝑟
𝑘𝑟
𝐸𝐼
𝜌𝐴
1.13 1.11 1.08 1.25 1.11 1.07 1.01 0.82 1.02
1.36 1.24 1.12 1.11 0.97 0.93 0.96 0.88 0.92
1.17 1.17 1.06 1.13 0.93 0.94 1.00 0.78 1.12
0.99 1.11 0.94 1.11 1.03 1.32 1.08 0.78 0.91
1.33 1.00 1.22 1.17 1.07 1.08 0.90 0.77 0.98
1.10 1.00 1.10 1.05 1.16 0.94 0.94 0.76 0.85
1.04 0.93 1.10 1.13 0.97 1.12 1.00 0.86 1.02
1.02 0.99 1.05 1.08 0.91 0.92 1.02 0.82 1.25
1.03 1.05 0.97 1.07 0.99 1.22 1.05 0.87 1.01
1.22 0.83 1.16 0.89 1.17 0.89 0.83 0.74 0.72

𝛼𝑏
0.96
1.11
0.96
0.98
1.06
0.87
0.84
0.92
0.95
0.88

𝛽𝑏
1.01
0.83
1.31
1.08
0.82
1.05
1.01
0.79
0.89
1.04

図 4-13 推定パラメータでの前後輪の路面凹凸推定値（RMS 値 0%ノイズ）

27

𝐹(𝑥)
3.64
3.85
3.80
3.85
3.84
3.70
3.96
3.95
3.92
3.82

𝑚𝑓
1
2
3
4
5
6
7
8
9
10

1.02
0.94
1.03
0.89
0.85
1.06
1.10
0.85
0.78
0.85

𝑚𝑟
1.17
0.95
1.05
1.25
0.85
1.04
1.11
0.89
0.86
0.85

表 4-8 パラメータ推定結果(RMS 値 15%ノイズ)
𝑐𝑓
𝑘𝑓
𝑘𝑡𝑓 𝑘𝑡𝑟
𝐿𝑓
𝑐𝑟
𝑘𝑟
𝐸𝐼
𝜌𝐴
1.01 1.03 1.11 1.34 1.04 1.22 1.06 0.78 0.97
1.14 1.32 0.98 1.30 1.11 1.29 1.12 0.83 0.73
1.32 1.32 1.05 1.15 1.19 1.28 1.01 0.73 0.76
1.00 1.36 0.89 1.27 0.94 1.33 1.20 0.81 1.10
1.23 1.14 1.14 1.26 1.07 1.16 0.98 0.80 1.05
1.12 1.21 0.96 1.22 1.09 1.22 1.08 0.84 0.84
1.11 1.13 1.03 1.21 0.98 1.05 1.03 0.80 1.04
1.20 1.18 1.09 1.24 0.93 1.05 1.02 0.82 1.14
1.08 1.46 0.98 1.38 0.93 1.12 1.16 0.86 0.79
1.23 1.14 1.14 1.26 1.07 1.16 0.98 0.80 1.05

𝛼𝑏
1.20
0.95
1.04
1.01
0.90
0.75
0.88
1.11
0.73
0.90

𝛽𝑏
0.78
1.09
1.14
1.15
1.14
0.86
1.10
1.05
0.78
1.14

図 4-14 推定パラメータでの前後輪の路面凹凸推定値（RMS 値 15%ノイズ）

28

𝐹(𝑥)
2.91
2.87
2.89
2.85
2.90
2.91
2.92
2.90
2.86
2.90

𝑚𝑓
1
2
3
4
5
6
7
8
9
10

1.20
1.27
1.06
1.10
0.91
1.16
0.87
1.03
1.33
0.91

𝑚𝑟
1.26
1.26
1.21
1.21
1.18
1.21
1.11
1.12
1.13
1.18

表 4-9 パラメータ推定結果(RMS 値 35%ノイズ)
𝑐𝑓
𝑘𝑓
𝑘𝑡𝑓 𝑘𝑡𝑟
𝐿𝑓
𝑐𝑟
𝑘𝑟
𝐸𝐼
𝜌𝐴
0.70 0.95 1.06 1.01 0.84 0.80 1.09 0.86 1.05
0.83 0.74 1.40 0.90 0.73 0.66 0.85 0.82 1.29
0.82 1.20 0.98 0.99 0.85 0.88 1.15 0.88 0.69
0.80 1.10 0.99 0.99 0.87 0.84 1.12 0.84 0.93
0.74 1.24 0.92 1.12 0.75 0.91 1.21 0.89 0.99
0.87 0.99 1.23 1.01 0.69 0.65 0.99 0.87 1.02
0.73 1.13 0.95 1.03 0.70 0.88 1.19 0.83 0.89
0.81 1.09 1.03 0.87 0.71 0.64 1.07 0.82 0.80
0.57 0.80 1.00 0.81 1.13 0.86 1.01 0.84 1.19
0.74 1.24 0.92 1.12 0.75 0.91 1.21 0.89 0.99

𝛼𝑏
0.81
1.09
0.88
1.11
0.96
0.85
1.26
1.03
1.04
0.96

𝛽𝑏
1.15
1.06
0.88
1.03
0.88
0.86
1.22
0.76
0.71
0.88

図 4-15 推定パラメータでの前後輪の路面凹凸推定値（RMS 値 35%ノイズ）

29

𝐹(𝑥)
3.02
3.01
3.05
3.06
3.05
3.01
3.01
3.01
3.09
3.05

第5章 まとめと今後の課題
数値実験により，車両振動データを用いた VBI システム同定にカルマンフィルタを適用し検討を行った．
本研究の特徴は，ピッチングを考慮したハーフカーモデルのパラメータと，その車両が走行した路面プロファイル，路面凹凸，そして橋梁パラメータが同時に推定できることである．
本提案は，橋梁をはじめとしたインフラモニタリングとしての実用化を考えた際に，橋梁上を走行する車両のキャリブレーションを事前に行う必要がないという利点が挙げられる．
本研究の提案手法の適用を検証した結果として，
1) ノイズがないことが既知である場合については，前後輪の路面凹凸を高精度に推定でき，また目的関数の形状としては設定に問題はないと考えられるものの，PSO 法の推定精度の問題により，感度の高い一部のパラメータのみを推定できる結果となった．
2) ノイズの大きさが未知である場合については，目的関数の設定や付加した観測ノイズの影響により推定される路面凹凸に前後輪で差が生じ，路面凹凸の推定にロバスト性向上を確認できたものの，パラメータ同定については精度が低い結果となった．
本研究にはいくつかの課題が挙げられる．1 つ目は，システムノイズ・観測ノイズの成分を試行錯誤的に定めていることである．これに関し Nagayama ら[6]は，Robbins-Monro
Algorithm[13]を導入し各成分を動的に見積もっていることから，本研究への適用により改善が期待できる．2 つ目は，目的関数の最尤値が全パラメータ正解値である値からずれてしまうことである．これは前述の通り，推定される路面凹凸に前後輪で差が生じることが原因として考えられるが，これにはカルマンフィルタによる路面プロファイルの推定精度の向上や目的関数の改善が必要である．1 つ目の課題に挙げたノイズ成分の動的見積もりや，周波数成分に着目した目的関数の設計をすることで解決できるのか試すことが考えられる．3 つ目は，各パラメータの推定精度のばらつきが大きいことである．目的関数への感度の高いパラメータに，他のパラメータの推定が引きずられ，初期収束を起こしてしまう可能性が高いと考えられる．ある程度は目的関数の出力結果から判断することができるが，最適化手法をより高度化することで改善を図る必要がある．

30

謝辞
本研究に際して，丁寧なご指導を頂きました山本亨輔先生と松島亘志先生に深謝いたします．加えて，日々の生活から研究活動までお世話になりましたフロンティア工学研究グルー
プの皆様には感謝の念に堪えません．そして生活を支えてくれた家族に心よりお礼申し上げ
ます．

31

参考文献
[1] 国土交通省「インフラメンテナンス情報 社会資本の老朽化の現状と将来」[2] H.L.Li, Z.R.Lu, J.K.Liu,”Identification of moving vehicular parameters based on
glowworm swarm optimization algorithm,” JVE INTERNATIONAL LTD.
VIBROENGINEERING PROCEDIA.,2015[3] Haoqi Wang, Tomonori Nagayama, Boyu Zhao, Di Su,” Identification of moving
parameters using bridge responses and estimated bridge pavement roughness,”Engineering Structures 153(2017), pp.57-70[4] Yuta Takahashi, Kyosuke Yamamoto, Yukihiko Okada, “空間特異モード角度を用いた
比較的軽度な橋梁損傷検知の可能性検討”,構造工学論文集 Vol. 65A, pp283-292, 2019.
[5] “粒子群最適化に基づく複数車両の振動データを用いた車両・橋梁・路面のパラメータ同定”,Kakeru Murakami, 筑波大学卒業論文, 2019[6] Tomonori Nagayama, Boyu Zhao, Kai Xue, “走行時の車体振動を利用したハーフカー
モデルの同定と路面縦断形状の推定”,土木学会論文集 E1(舗装工学), Vol. 75, No.1, 1-16,
2019.
[7] Junki Nakasuka, Wang Haoqi, Tomonori Nagayama, “車両応答計測を利用した路面
プロファイル推定に基づく橋梁固有振動数の抽出”, 構造工学論文集 Vol.64A(2018 年 3月)[8] Boyu Zhao, Tomonori Nagayama, Kai Xue, “Road profile estimation, and its numericaland experimental validation, by smartphone measurement of the dynamic responses
of an ordinary vehicle”, Journal of Sound and Vibration 457(2019)92-117[9] M.F.Green, D.Cebon,”Dynamic interaction between heavy vehicles and highway
bridges,”Computers & Structures 62(1997),pp.253-264[10] Rauch, H.E., Tung, F. and Striebel, C. T.”Maximum likelihood estimates of linear
dynamic systems”, AIAA Journal,Vol.3,No.8,pp.1445-1450,1965[11] J.Kennedy, R.C.Eberhart,”Particle swarm optimization,”in Proc. of IEEE
Int.Conf.on Neural Networks,pp.1942-1948,1995[12] H.Matsushita, Y.Nishio and T.Saito, ”Particle swarm optimization with novel concept
of complex network”,in Proc. of NOLTA,pp.197-200,2010.
[13] Robbins,H. and Monro, S.”A stochastic approximation method”, The Annals of
Mathematical Statistics, Vol.22, No.3, pp.400-407, 1951.

32