令和三年度

筑波大学工学システム学類卒業研究論文

複数車両を用いた VBI システム同定のパラメータ更新方法に関する数値的検討

学籍番号201811171

氏

名

藤原

幸代

指導教員

山本

亨輔

201811171 FUJIWARA - 1

複数車両を用いた VBI システム同定のパラメータ更新方法に関する数値的検討
Numerical Study of Parameter Update Method
in VBI System Identification Using Multiple Vehicles

藤原幸代

Sachiyo Fujiwara
（指導教員

山本亨輔）

This study aims to improve the estimation accuracy of VBISI (Vehicle-Bridge Interaction System
Identification) technology is an on-going method using only vehicle vibration and position data forsimultaneously estimating the vehicle’s and bridge’s mechanical parameters, and road roughness. In thecurrently devised methods, the estimation accuracy is degraded when even a small amount of noise is added
to the vehicle vibration data. As a more accurate method, we attempt to improve the estimation accuracy
when multiple vehicles are moved on the same bridge and shared the bridge’s mechanical parameters. The
Nelder-Mead method is used as the optimization method. In the presence of noise, the vehicle’s and bridge’s
mechanical parameters, and road roughness are estimated with good accuracy. On the other hand, in theabsence of noise, data errors tend to accumulate, and there is no improvement even with multiple vehicles.
The accuracy of the estimation is also affected by differences in vehicle models. The Nelder-Mead method
is computationally efficient and has reasonably good estimation accuracy, but it sometimes fails to reachthe correct solution, and there is room for improvement.

1. はじめに
日本には多くの橋梁が存在しており，その多くった様々な業務・事業が急増し，リソースが慢性的に不足している地方自治体の運営をさらに圧迫していくと予想される．

は高齢化が進んでいる．国土交通省によれば，現在

現在，橋梁の健全性は 5 年に 1 度の詳細点検が日本には道路橋が 72 万基以上存在しており[1]，供義務付けられており[1]，熟練技術者による目視点検用開始から 50 年以上が経過しつつある橋梁の割合

や打音検査等により判定が行われる．しかし，少子

は 2029 年時点で 52％にまで増加すると予測され高齢化に伴う技術者不足は深刻な問題で，少ない

る．経年劣化に伴う社会的コスト，すなわち，維持

人材で効率よく橋梁を点検しなければならない．

管理に必要な技術人材の育成・確保や，劣化・損傷

そのような方法として振動モニタリング [2]が検討箇所の早期発見・監視・修繕，代替橋梁の建設といされている．例えば，加速度センサ等を橋梁に設置

201811171 FUJIWARA - 2

し，実測値に基づいて橋梁の特微量を推定する．点できる．車が直進する場合，前輪と後輪で推定さ検技術者の技能に左右されず客観的に評価可能でれる路面凹凸は一致すると考えられるので，そのあり，技術者不足に対して有効な手段と期待され一致度を目的関数とする．本研究では，このようる．しかし，橋梁の数は膨大であるため，すべてのな手法を，VBISI（VBI System Identification）技術

橋梁で実施するには，大量のセンサや電源の設置，と呼ぶ．VBISI 技術により，車両振動データさえそれに伴う交通整理など，費用や労力が大きく現取得できれば，車両・橋梁の力学的パラメータと実的ではない．そこで，本研究では間接的アプロー

路面凹凸が同時に点検できる可能性が示されたとチを試みる．この手法は車両に加速度センサを設言える．村上モデルでは車両速度振動，車両変位

置し，得られた振動から橋梁の特徴量を推定する．

振動を数値積分により求める．数値積分は誤差を

直接的な橋梁振動モニタリングと比較して，精度蓄積しやすく，わずかでも車両振動データにノイには課題が残るものの，交通整理などを必要としズが加わると推定精度が低下してしまうという課

ない簡易的な健全性診断手法である．

題がある．そのため，ノイズに対しても頑健なモ

既往の研究では，Yang ら[2]が橋梁上を走行する

デルに拡張する必要がある．

車両のパワースペクトルから橋梁の固有振動数のそこで本研究では，同じ橋梁上を複数台の車両抽出が試みた．結果，大まかな推定はできることを走行させ，車両・橋梁の力学的パラメータと路面を実験的に示した．固有振動数は、質量や剛性の凹凸の推定精度の向上を，数値実験によって試み

情報を含んでおり、点検指標としても期待できる．その際，橋梁の力学的パラメータを平均化し，る。但し、健全時データとの比較が必要なため、車両間で共有し，目的関数も一つに集約する．車両

社会コストが大きい。そこで，山本ら[3]は走行車

2 台が同じ橋梁上を走行する場合，まず，前・後輪両の振動応答に基づき橋梁のモード形状を推定す

における計 4 種類の推定路面凹凸の平均求める．

る手法を提案した．一方，長山ら[4]は車両振動か

次に各路面凹凸と平均との差をとり，その二乗誤

ら車両パラメータと路面縦断形状を同時推定する差の総和を目的関数として一つに集約し，パラメ

手法を提案している．また，村上モデル[5][6] でータを更新する．橋梁は健全時と損傷時の 2 種類

は，長山らの手法を拡張し，車両・路面に加えてを用い，損傷は橋の中心に起こる．得られた車両・橋梁も考慮しており，クオーターカーモデルの車

橋梁の力学的パラメータと路面凹凸を正規化し，両 2 台から得られた振動データに粒子群最適化をパラメータ推定精度の変化と橋の損傷の検知の有適用することで，車両・橋梁の力学的パラメータ

無を検証する．

（質量・減衰・剛性）と路面凹凸を同時に推定する．村上モデルでは，ランダムに車両・橋梁パラメータを仮定し，計測された車両振動データ組み

2. 数値実験

合わせて VBI（Vehicle-Bridge Interaction，車両－

VBISI 法では，車両振動から，車両・橋梁の力学

橋梁相互作用）システムの運動方程式を作成す

的パラメータ（質量・減衰・剛性）と路面凹凸を推る．VBI システムの運動方程式からは路面凹凸が定する．今回は数値計算によって橋梁走行中の車推定できるので，最も尤もらしい路面凹凸を推定

両振動データを再現し，提案手法の適用性を検証できたパラメータを採用することで，車両・橋梁

する．まず数値実験について説明する．

の力学的パラメータと路面凹凸の同時推定が実現

201811171 FUJIWARA - 3

2.1.

橋梁は曲げのみを考慮した 1 次元オイラー・ベ

車両モデル

本研究では，車両は図１に示すハーフカーモデ

ルヌーイ梁とし，図２に示す．橋梁の損傷は中心部ルを採用する．このモデルはサスペンション𝑖（剛分（図２の黄色部分，要素 4）とする．橋梁の変位性𝑘𝑠𝑖 ，減衰𝑐𝑠𝑖 ）によって，質量𝑚𝑠 のばね上剛体と振動を𝑦(𝑥, 𝑡)とし，曲げ剛性を𝐸𝐼，単位長さ当たり

質量𝑚𝑢𝑖 のばね下質点が接続されている．ばね上剛

質量を𝜌𝐴，外力を𝑝とすると，橋梁の運動方程式は

体は車体を，ばね下質点はタイヤと車軸をそれぞ式(8)で表される.

れモデル化している．重心から前輪・後輪までの距
離を𝑑𝑖 とする．添字の𝑖は車軸に対応し，前輪が 1，
後輪が 2 となる．車両の運動方程式は，
𝐌𝑣 𝒛̈ (𝑡) + 𝐂𝑣 𝒛̇ (𝑡) + 𝐊𝑣 𝒛(𝑡) = 𝑭𝑣

𝑑1 𝑚𝑠
𝑑1 + 𝑑2𝐼−
𝑑1 + 𝑑2𝑐𝑠2−𝑑2 𝑐𝑠2
0
−𝑐𝑠2

(1)

反力𝑅A ，𝑅B から成る．支点の位置をそれぞれ𝑥A =

すると
2

(2)

𝑖=1

𝑚𝑢1

−𝑐𝑠1
−𝑑1 𝑐𝑠1
𝑐𝑠1
0−𝑘𝑠1−𝑑1 𝑘𝑠1
𝑘𝑠1 + 𝑘𝑢10(9)𝑝 = ∑ 𝛿(𝑥 − 𝑥𝑖 )𝑃𝑖 (𝑡) + 𝛿(𝑥)𝑅A + 𝛿(𝑥 − 𝐿)𝑅B

𝑚𝑢2 ]と書ける．𝛿(𝑥)はディラックのデルタ関数である．

−𝑐𝑠2
𝑑2 𝑐𝑠20 ]𝑐𝑠2(8)式を数値的に解くために，本研究では有限要(3)素法を適用する．(8)式の重み付き残差式は，𝐿∫ 𝜔 (𝜌𝐴
𝑘𝑠1
𝑘𝑠2
𝑑1 𝑘𝑠1 −𝑑2 𝑘𝑠2
𝐊𝑣 = [
−𝑘𝑠1
0
0
−𝑘𝑠2

(8)

0，𝑥B = 𝐿とし，車両の接地力𝑃𝑖 (𝑡)の位置を𝑥𝑖 (𝑡)と[𝑐𝑠1
𝑑1 𝑐𝑠1
𝐂𝑣 = [ −𝑐
𝑠1
0

𝜕2
𝜕2
𝐸𝐼
(
𝑦(𝑥, 𝑡)) = 𝑝(𝑥, 𝑡)
𝜕𝑥 2
𝜕𝑥 2

ここで，外力𝑝(𝑥, 𝑡)は，車両の接地力𝑃𝑖 (𝑡)と支点と表すことができる．ここで，
𝑑2 𝑚𝑠
𝑑1 + 𝑑2𝐼
𝐌𝑣 =𝑑1 + 𝑑2𝜌𝐴𝑦̈ (𝑥, 𝑡) +−𝑘𝑠2
𝑑2 𝑘𝑠2]0
𝑘𝑠2 + 𝑘𝑢2

0

(4)

(5)

𝑭𝑣 = [0 0 𝑘𝑢1 𝑢1 𝑘𝑢2 𝑢2 ]𝑇

(6)

(10)

となり，𝜔は重みを表す．式(10)弱形式化すると，𝐿𝒛(𝑡) = [𝑧𝑠1 𝑧𝑠2 𝑧𝑢1 𝑧𝑢2 ]𝑇

𝜕2𝑦
𝜕4 𝑦
+ 𝐸𝐼 4 − 𝑝) d𝑥 = 0
2𝜕𝑡𝜕𝑥∫ (𝜌𝐴𝜔0𝜕2𝑦𝜕2𝜔 𝜕2𝑦+𝐸𝐼− 𝑝) d𝑥 = 0
𝜕𝑡 2
𝜕𝑥 2 𝜕𝑥 2(11)

橋梁振動𝑦(𝑥, 𝑡)は，Hermite 補間関数を用いた一であり，𝐌𝑣 ,𝐂𝑣 ,𝐊𝑣 はそれぞれ車両の質量, 減衰，次元有限要素梁モデルによって離散化する．ここ

剛性マトリクス，𝒛̇ ，𝒛̈ は𝒛(𝑡)のそれぞれ一階・二階で正規化された要素座標系𝑋（−1 ≤ 𝑋 ≤ 1）に対し時間微分を表しており，速度振動と加速度振動でて，Hermite 基底を以下のように定義する．

ある．ばね上剛体の慣性モーメント𝐼は，回転中心1
𝜙1 (𝑋) = (𝑋 − 1)(𝑋 − 1)(𝑋 + 2)
4
Δ𝑥
(𝑋 − 1)(𝑋 − 1)(𝑋 + 1)
𝜙2 (𝑋) =81
𝜙3 (𝑋) = − (𝑋 + 1)(𝑋 + 1)(𝑋 − 2)4Δ𝑥{𝜙4 (𝑋) = 8 (𝑋 − 1)(𝑋 + 1)(𝑋 + 1)

によって値が異なる．ここでは，重心を回転中心と仮定し，
𝐼 = 𝑚𝑠 𝑑1 𝑑2

(7)

とした．

(12)

本数値実験では，既往研究[4][7]で用いられた車両

全体座標系𝑥が節点𝑥𝑗 と節点𝑥𝑗+1 から構成されるモデルを参考に車両の力学的パラメータを定めた．

梁要素𝑗の内部にある時，𝛥𝑥 = 𝑥𝑗+1 − 𝑥𝑗 である．基

具体的なパラメータを表１に示す．

底関数ベクトル𝑵(𝑥)を，その成分が以下になるよう定義する．

2.2.

橋梁モデル

201811171 FUJIWARA - 4

表 1 車両パラメータ
パラメータ10[t]車両

18[t]車両

車体質量𝑚𝑠 [kg]8.31×103

1.66×104

ばね下質量𝑚𝑢1[kg]4.96×102

7.00×102

ばね下質量𝑚𝑢2[kg]7.51×102

7.00×102減衰𝑐𝑠1 [kg/s]

2.42×104

1.00×104減衰𝑐𝑠2 [kg/s]

2.90×1041.00×104剛性𝑘𝑠1 [N/m]4.56×104

4.00×104

剛性𝑘𝑠2 [N/m]

4.10×105

4.00×104

タイヤ剛性𝑘𝑢1 [N/m]4.79×106

1.75×106

タイヤ剛性𝑘𝑢2 [N/m]

4.31×106

1.75×106

車間距離𝑑[m]

4.40

4.75

車軸間距離𝑑1 [m]

1.215

2.375

車軸間距離𝑑2 [m]

3.185

2.375

𝑑1

𝑑2

𝐷
𝑘𝑠2

𝑐𝑠2

𝑘𝑠1

𝑚 𝑢2

𝑚 𝑢1

𝑘𝑢2

𝑘𝑢 1

𝑐𝑠1

図 1 車両モデル

図２ 橋梁モデル

表２ 橋梁パラメータ
曲げ剛性𝐸𝐼[Nm3]

1.56×1010

損傷部分の曲げ剛性𝐸𝐼[Nm3]

7.8×109単位長さあたり質量𝜌𝐴[kg/m]

4.4×103

スパン長𝐿[𝑚]30

レイリー減衰𝛼0.7024

レイリー減𝛽

0.0052

要素数
𝑁2𝑗−1 (𝑥) = 𝜙1 (𝑋)𝑁2𝑗+0 (𝑥) = 𝜙2 (𝑋)𝑁2𝑗+1 (𝑥) = 𝜙3 (𝑋){𝑁2𝑗+2 (𝑥) = 𝜙4 (𝑋)

図 3 車両橋梁相互作用7
𝝎T (𝐌𝑏 𝒚̈ (𝑡) + 𝐊𝑏 𝒚(𝑡) − 𝑭(𝑡)) = 0
(13)

(16)

𝐌𝑏 ，𝐊𝑏 はそれぞれ橋梁の質量，剛性マトリクスを示している．𝑭(𝑡)は各節点における外力（集中荷重

なお，要素外ではいずれの成分もゼロとする．一方，と力のモーメント）を成分とする外力ベクトルで

変形量ベクトル𝒚(𝑡)は，節点𝑥𝑗 における橋梁たわみ

ある．

𝑦(𝑥𝑗 , 𝑡)とたわみ角𝜃(𝑥𝑗 , 𝑡)を用いて{

𝑦2(𝑗−1)+1 (𝑡) = 𝑦(𝑥𝑗 , 𝑡)𝑦2(𝑗−1)+2 (𝑡) = 𝜃(𝑥𝑗 , 𝑡)

𝐿

𝐌𝑏 = ∫ 𝑵𝑵T d𝑥(14)

𝜕𝑵 𝜕𝑵Td𝑥
0 𝜕𝑥 𝜕𝑥𝐿

𝐊𝑏 = ∫とする．この時，𝑦(𝑥, 𝑡)の近似解は
𝑦(𝑥, 𝑡) = 𝑵(𝑥) ⋅ 𝒚(𝑡)(15)となる．同様に，重みも𝜔(𝑥) = 𝑵(𝑥) ⋅ 𝝎とし，重み付き残差式（(11)式）に代入すると，以下を得る．

(17)

0

(18)

さらに(16)式において，レイリー減衰を考慮し，
𝐂𝑏 = 𝛼𝐌𝑏 + 𝛽𝐊𝑏(19)を導入する．任意の𝝎に対して，式(16)の等号成立201811171 FUJIWARA - 5

条件を求めると，以下に示す有限要素式が得られ

橋梁振動𝑦̃𝑖 (𝑡) = 𝑦(𝑥𝑖 (𝑡), 𝑡)である．つまり，橋梁振る．

動𝒚(𝑡)は橋梁上に固定された各節点の変形量ベク
𝐌𝑏 𝒚̈ (𝑡) + 𝐂𝑏 𝒚̇ (𝑡) + 𝐊𝑏 𝒚(𝑡) = 𝑭(𝑡)

(20)

である．

トルであるから，離散化で用いたのと同じ基底を用いて，車軸位置𝑥𝑖 (𝑡)における橋梁変位𝑦̃𝑖 (𝑡)に変本数値実験で扱う車両・橋梁の力学的パラメー換する必要がある．変換マトリクス𝐋(𝑡)を，タの値を表１，表２に示す

𝐋(𝑡) = [𝑵(𝑥1 (𝑡)) 𝑵(𝑥2 (𝑡))](23)とおくと，橋梁プロファイルは，
2.3.

̃ (𝑡) = 𝑳T (𝑡)𝒚(𝑡)𝒚

車両-橋梁の相互作用

車両と橋梁の応答は互いの出力を入力とする相

(24)

である．

互作用によってモデル化される．VBI（VehicleBridge Interaction，車両－橋梁相互作用）システム

接地力

2.3.2.

は図３のようにモデル化でき，本稿ではその動的

次に接地力について説明する．橋梁への入力で挙動を数値的に再現する．まず，仮定した路面プロ

ある接地力は，タイヤ（剛性𝑘𝑢𝑖 ）に働く復元力にファイルと表 1 に示す車両パラメータを入力とし相当する．ただし，車両の運動方程式（(1)式）ではて Newmark-𝛽法を用いて車両振動を計算する．そ

つり合い位置を基準としているため，重力項が消して，得られた車両振動と路面プロファイルからえていることに注意を要する．復元力の計算では橋梁への接地力を求め，表２に示す橋梁パラメー

自然長を基準としているため重力による影響を考タと Newmark-𝛽法を用いて橋梁振動を計算する．

慮する必要がある．前輪と後輪のそれぞれの接地

この橋梁振動を路面プロファイルに足し合わせて力は，

新たな路面プロファイルとし，表 1 に示す車両パ
ラメータを入力として Newmark-𝛽法を用いて車両振動を求め，再度同じ工程を繰り返すことで必要なデータを得る．

𝑉1 (𝑡) =𝑑2 𝑚𝑠(𝑔 − 𝑧̈𝑠1 ) + 𝑚𝑢1 (𝑔 − 𝑧̈𝑢1 )𝑑1 + 𝑑2(25)

𝑑1 𝑚𝑠(𝑔 − 𝑧̈𝑠2 ) + 𝑚𝑢2 (𝑔 − 𝑧̈𝑢2 )𝑉2 (𝑡) =𝑑1 + 𝑑2である．橋梁に働く外力ベクトルは

2.3.1.

入力プロファイル

𝑭(𝑡) = 𝐋(𝑡)[𝑉1 (𝑡) 𝑉2 (𝑡)] + 𝑯(𝑡)車両ばね下に入力される強制変位を，入力プロ(26)である．𝑯(𝑡)は支点反力を表す．

ファイルと呼ぶこととする．入力プロファイル𝒖(𝑡)̃ (𝑡)は，路面プロファイル𝒓(𝑡)と橋梁プロファイル𝒚の和で与えられ，式(21)で表される．

2.3.3.

Newmark-𝛽法

車両振動𝒛(𝑡)，橋梁振動𝒚(𝑡)はそれぞれの運動方(21)

程式((1)式，(8)式)に Newmark-β 法を用い，数値積

ここで路面プロファイルとは車軸位置における路分をして求めている．以下に Newmark-β 法を示す．

面凹凸を表す．路面凹凸が𝑅(𝑥)，車軸位置が𝑥𝑖 (𝑡)で運動方程式が MCK 系で以下のように表される

̃ (𝑡)𝒖(𝑡) = 𝒓(𝑡) + 𝒚

ある時，

とする．
𝑟𝑗 (𝑡) = 𝑅 (𝑥𝑗 (𝑡))(22)である．𝑟𝑖 (𝑡)は𝒓(𝑡)の成分である．
一方，橋梁プロファイルとは車軸位置𝑥𝑖 (𝑡)での𝐌𝒔̈ (𝑡) + 𝐂𝒔̇ (𝑡) + 𝐊𝒔(𝑡) = 𝑸(𝑡)(27)時刻関数である𝒔(𝑡) を離散化して𝒔𝑘 とする．𝒔𝑘 は車両または橋梁の変位であり𝛥𝑡を時間刻みとして，

201811171 FUJIWARA - 6

表３ 車両の固有振動数

表４ 橋梁の固有振動数

10t 車両18t 車両健全時[Hz]

損傷時[Hz]1.12

6.20×10-1

3.37

2.98

1.80

9.96×10-1

1.26×101

8.84

1.64×101

8.84

1

1.45×10

1.43×101

3.63×101

3.34×101

動，路面凹凸，橋梁振動，車両と橋梁のパワースペ
クトルの結果を図４から図７に示す．車両と橋梁(28)

𝒔𝑘 = 𝒔(𝑘𝛥𝑡)である．Newmark-β 法ではで一定である．

𝒔̇ 𝑘 = 𝒔̇ 𝑘−1 + (1 − 𝛾)𝛥𝑡𝒔̈ 𝑘−1 + 𝛾𝛥𝑡𝒔̈ 𝑘(29)

1
𝒔𝑘 = 𝒔𝑘−1 + 𝛥𝑡𝒔̇ 𝑘−1 + 𝛥𝑡 2 {( − 𝛽) 𝒔̈ 𝑘−1 + 𝛽𝒔̈ 𝑘 }
2

(30)

1

1

2

4

と仮定し，本研究では𝛾 = ，𝛽 = として，(27)式に Newmark-β 法を適用すると，−1

図４より，バネ下の車両振動は路面凹凸と似た波形を示すことがわかる．前・後輪で描いた振動も一致した．
図５より，10t 車両に比べて，18t 車両走行時では，橋梁振動が大きいことが確認できる．これは車(31)𝒔̈ 𝑘 = 𝐀 𝒃𝑘

重の違いによる影響である．
車両の固有振動数は 1.1[Hz]，1.8[Hz]，12.6[Hz]，と書ける．ここで，𝐀 = [𝐌 +の固有振動数を表４と表５に示す．車速は 10[m/s]

(𝛥𝑡)2
𝛥𝑡
𝐂+
𝐊]
2
4

(32)

16.4[Hz]である．図６より，1.24[Hz]，9.91[Hz]，12.1[Hz]，17.7[Hz]で車両パワースペクトルの卓越

𝒃𝑘 = 𝑸𝑘 + 𝐂𝒃1 + 𝐊𝒃2(33)が確認できる．固有振動数付近で卓越しており，ズ

𝑸𝑘 = 𝑸(𝑘𝛥𝑡)(34)

レは 0.1 から 1.3[Hz]である．

𝒃1 = −𝒔̇ 𝑘−1 −𝒔̈ 𝑘−1𝛥𝑡2

𝒃2 = −𝒔𝑘−1 − 𝒔̇ 𝑘−1 𝛥𝑡 −𝒔̈ 𝑘−1
(𝛥𝑡)2
4

(35)

(36)

である．

橋梁の固有振動数は健全時 3.37[Hz]，14.5[Hz]，損傷時 2.98[Hz]，14.3[Hz]である．図７より，健全時 3.24[Hz]，14.9[Hz]，損傷時 2.58[Hz]，14.9[Hz]で橋梁パワースペクトルの卓越が確認できる．固有振動数付近で卓越しており，
ズレは 0.3 から 0.6[Hz]である．橋梁の卓越振動数にはほとんど変化が見2.4.

数値実験の結果

数値シミュレーションを行い，得られた車両振られない．

201811171 FUJIWARA - 7

前輪

18t

後輪

前輪

Disp[m]後輪

Disp[m]

10t

[m]

[m]

(a)バネ上振動（10t車両）

後輪

前輪

18t

後輪

前輪

Disp[m]Disp[m]10t(d)バネ上振動（18t車両）

[m]

[m]

(b)バネ下振動（10t車両）

Disp.[m]Disp.[m](e)バネ下振動（18t車両）

区間

区間

[m]

[m]
Traveling position[m]Traveling position[m](c)路面凹凸(g)路面凹凸

図４ 車両のバネ下振動(a)と路面凹凸(b)10tDisp[m]

前輪進

後輪退出

前輪退出

10t

18t

区間
前輪位置[m]

図５ 橋梁変位振動（15[m]地点）

18t

後輪退出

201811171 FUJIWARA - 8

1.24[Hz]1.12[Hz] 1.8[Hz]16.4[Hz]10t

17.7[Hz]12.2[Hz]

PSD[

PSD[

健全
損傷

16.4[Hz]10t

12.2[Hz]健全

9.91[Hz]損傷

0.62[Hz] 0.996[Hz] 8.84[Hz]0.92[Hz]
0.92[Hz]
8.75[Hz]

18t

健全
損傷

18t

8.75[Hz]
9.23[Hz]

PSD[]

12.6[Hz]17.7[Hz]

0.62[Hz] 0.996[Hz] 8.84[Hz]

PSD[

]

12.6[Hz]]]1.12[Hz] 1.8[Hz]健全
損傷

Frequency [Hz]Frequency [Hz](a)前輪バネ下質点(b)後輪バネ下質点

図６ 車両のパワースペクトル3.24[Hz]

14.9[Hz]
14.9[Hz]

14.3[Hz]
2.98[Hz]
14.5[Hz] 33.4[Hz]3.37[Hz]]

2.58[Hz]36.2[Hz]

3.24[Hz]

健全
損傷

36.2[Hz]

2.58[Hz]

PSD[

PSD[

]

14.3[Hz]
2.98[Hz]
14.5[Hz] 33.4[Hz]3.37[Hz]

健全
損傷

Frequency[Hz]Frequency[Hz](a)7.5[m]地点(b)15[m]地点

図７ 橋梁のパワースペクトル3. 提案手法
3.1.

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

車両振動から路面凹凸を推定する為に，長山らの離散時間拡張状態空間モデル[4] を導入し，定義した状態空間モデルに線形カルマンフィルタを適る．
𝒛(𝑘𝛥𝑡)
𝒛̇ (𝑘𝛥𝑡)
𝒁𝑘 = {}𝒖(𝑘𝛥𝑡)𝒖̇ (𝑘𝛥𝑡)(37)𝑧̈𝑠1 (𝑘𝛥𝑡)𝑧̈ (𝑘𝛥𝑡)𝒔𝑘 = 𝑠2𝑧𝑠1 (𝑘𝛥𝑡){𝑧𝑠2 (𝑘𝛥𝑡)}(38)

用する．カルマンフィルタとは，状態ベクトルの推状態ベクトルは，車両の鉛直変位振動と速度成分，

定を，観測データと力学モデルに基づいて，効率的

入力プロファイルおよび速度成分で与えられる．

に行う方法である．

本研究における離散時間拡張状態空間モデルは次のように書ける．

3.1.1.

状態空間モデルの構築

状態ベクトル𝒁𝑘 と観測ベクトル𝒔𝑘 を次式で与え̅ 𝒁𝑘−1 + 𝝎𝑘
𝒁𝑘 = 𝐕(39)

𝒔𝑘 = 𝐇𝒁𝑘 + 𝝐𝑘(40)

201811171 FUJIWARA - 9

出出 に欠損があるため
全ての
を辿れない

に対して出

がない

図８ 可観測性の概要

̅ は(35)式で与えられ，𝝎𝑘 , 𝝐𝑘 はシステム
行列指数𝐕
ノイズと観測ノイズを表し，𝝎~𝑵(𝟎, 𝐐), 𝝐~𝑵(𝟎, 𝐑)
である．

である．
システムノイズ𝝎𝑘 と観測ノイズ𝝐𝑘 の分散共分散
行列は𝐐，𝐑とする．

̅ = expm[𝐕𝛥𝑡]𝐕(41)= 𝐔 diag(exp(𝐃)) 𝐔−1

𝐐 = E[𝝎𝑘 𝝎T𝑘 ](45)

𝐑 = E[𝝐𝑘 𝝐T𝑘 ](46)

ただし，𝐕は，
𝐎4×4−𝐌v−1𝐊 v
𝐕=[
𝐎2×4
𝐎2×4

𝐈4×4
−𝐌v−1 𝐂v
𝐎2×4
𝐎2×4

𝐎4×2
𝐌v−1𝐅v
𝐈2×2
𝐎2×2

3.1.2.

𝐎2×2
𝐎2×2
]
𝐎2×2
2×2
𝐎

(42)

可観測性分析

出力を観測した状態空間モデルから内部の状態を推定可能であるかの尺度として可観測性分析をであり，𝐔および𝐃は，𝐕𝛥𝑡を対角化したときのモー行う必要がある．構築した状態空間モデルの可観

ド行列および対角行列であり，測性を調べるための必要十分条件は次式となる．

𝐕𝛥𝑡 = 𝐔𝐃𝐔−1

(43)

𝐇
𝐇𝐕
]) = 𝑛⋮
𝐇𝐕 𝑛−1と表される．また，𝐇は，
−

𝑘𝑠1
𝑚𝑠
0

𝑘𝑠1
𝑚𝑢1

0

𝑘𝑠2
𝑚𝑠

0

𝑘𝑠2
𝑚𝑢2

(𝑘𝑠1 + 𝑘𝑢1)
𝑚𝑢1

0

0

0

𝑘𝑠2
𝑚𝑠

𝑐𝑠1
𝑚𝑠
0

𝑐𝑠1
𝑚𝑠

0
−

T

0

𝑘𝑠1
𝑚𝑠

−
𝐇=

−

𝑟𝑎𝑛𝑘 ([

𝑐𝑠2
𝑚𝑠
0

0

𝑐𝑠1
𝑚𝑠

0

0

−

0
𝑐𝑠1
𝑚𝑢1
0
−

𝑐𝑠1
𝑚𝑢1
0

𝑘𝑢1
𝑚𝑢1

0

0

0

0
[ 0

0
0

0
0

−

(47)

𝑛は状態ベクトルの自由度である．本研究の状態空間モデルは𝑛 = 8であり，可観測性は満たしていな
い．可観測性を満たしていないとは，出力に欠損が(𝑘𝑠2 + 𝑘𝑢2)
𝑚𝑢2

あり，全ての入力を辿れないモデルを指す．図８に

0
𝑐𝑠2
𝑚𝑢2

(44)

可観測性を満たしていない場合の概要を示す．

0
−

𝑐𝑠1
𝑚𝑢2

3.1.3.

̂𝑘 とおく．𝒁̂𝑘−1 が得推定された状態ベクトルを𝒁

0
𝑘𝑢2
𝑚𝑢2
0
0

線形カルマンフィルタ

̂𝑘 の候補𝑋は事前分布をられた場合，(39)式から𝒁]

𝐏𝑘−1 として，平均𝜇𝑎 ，分散𝚺𝑎 の正規分布に従う．一201811171 FUJIWARA - 10

̂𝑘 の候方，𝒔𝑘 が得られた時，(40)式から導出される𝒁ータを式(55)に代入して接地力𝑉1 (𝑡)，𝑉2 (𝑡)を求め補𝑌も同様に平均𝜇𝑏 ，分散𝚺𝑏 の正規分布に従う．

る．橋梁パラメータもランダムに仮定し，式(8)に有̅𝒁
̂𝑘−1
𝜇𝑎 = 𝐕(48)

限要素法を適用し離散化して，Newmark-β 法を用

̅ 𝐏𝑘−1 𝐕
̅T + 𝐐
𝚺𝑎,𝑘 = 𝐕(49)

いると，橋梁振動𝒚(𝑡)を求めることができる．式

𝜇𝑏 = 𝐇−1 𝒔𝑘(50)

̃ (𝑡)を得る．
(54)に代入して橋梁プロファイル𝒚𝚺𝑏,𝑘 = 𝐇−1 𝐑𝐇−T

(51)である．初期値として，

ここで，先に求めていた𝒖(𝑡)と橋梁プロファイル
̃ (𝑡)の差をとると，路面プロファイル𝒓(𝑡)が得られ
𝒚̂0 = 𝒁0𝒁(52)

る．得られた路面プロファイルを𝒙(𝑡)で位置同期す𝐏0 = (𝒁0 − 𝐸[𝒁0 ])(𝒁0 − 𝐸[𝒁0 ])T(53)

ると，前輪𝑅1 (𝑡)と後輪𝑅2 (𝑡)の路面凹凸が得られる．

を与える．

車輌が直進していると仮定すれば，前輪と後輪は

̂𝑘 は 2 つカルマンフィルタ適用後の状態ベクトル𝒁

同じ路面凹凸上を走行するため，推定される路面の候補𝑋，𝑌が正規分布に従うと仮定し，最尤値と凹凸は一致するはずである．しかし，力学的パラメ

して推定される．以下にその式を示す．

ータのランダムに与えているため，一致しない．そ

−1
−1
−1
̂𝑘 = (𝚺𝑎,𝑘𝒁+ 𝚺𝑏−1 ) (𝚺𝑎,𝑘
𝜇𝑎 + 𝚺𝑏−1 𝜇𝑏 )こで，前輪と後輪で推定される路面凹凸の二乗誤

̅𝒁̂𝑘−1 } + [𝐆𝑘 𝐇]{𝐇−1 𝒔𝑘 }= [𝐈 − 𝐆𝑘 𝐇]{𝐕

(54)

̅𝒁
̂𝑘−1 } + 𝐆𝑘 𝒔𝑘= [𝐈 − 𝐆𝑘 𝐇]{𝐕

面凹凸が一致するようにパラメータを更新できれここで，𝐆𝑘 = 𝚺𝑎,𝑘 𝐇T (𝐇𝚺𝑎,𝑘 𝐇T + 𝐑)

差を最小化する最適化問題を考える．算出する路

ば，最終的にパラメータも正解値に近づくと期待
−1

−1
−1
𝐏𝑘 = (𝚺𝑎,𝑘+ 𝚺𝑏−1 ) = [𝐈 − 𝐆𝐇]𝚺𝑎(55)

できる．このようなパラメータ同定手法を「VBI シ(56)

ステム同定」と呼ぶこととする．

である．
3.2.1.
3.1.4.

路面凹凸推定非線形最適化問題

車両と橋梁の力学的パラメータは全て正の値でカルマンフィルタを用いた状態ベクトルの推定あるため，非負条件をもつ．また車間距離𝑑は車軸

̃ (𝑡)をにより得られた𝒖(𝑡)から，橋梁プロファイル𝒚

間距離𝑑𝑖 の和で表されるため，等号制約をもつ．そ差し引くことで，路面凹凸𝒓(𝑡)を推定する．

こで本研究では，ペナルティ法を用いて，制約付き

複数車両では車軸間距離が各車両で違うため，

非線形最適化問題を，制約なし非線形最適化問題

推定する路面凹凸の長さは，車軸距離が短い方にに変換し Nelder-Mead 法により最適解を探索する．

合わせたものを採用した．

非線形計画問題において，勾配法を用いたいとき
は，テイラー展開により目的関数の 2 次近似を求

3.2.

車両-橋梁相互作用システム同定

め，近似最適解を探索することが多い．しかし，勾まず，計測データとして得られた車両加速度𝒛̈ (𝑡)配法は目的関数の傾きを元に次の更新点を求めるに Newmark-β 法を適用し，𝒛̇ (𝑡)，𝒛(𝑡)を求める．ラため，傾きの非常に小さい最適化問題においては，

ンダムに車両パラメータを仮定すると，𝐌𝑣 ，𝐂𝑣 ，

効率が悪い．そこで，傾きを必要としない Nelder-

𝐊𝑣 も求められる．ここでカルマンフィルタを用い

Mead 法採用する．

て状態ベクトルを推定し，𝑢1 ，𝑢2 を求める．
次に車両振動データ𝒛̈ (𝑡)と仮定した車両パラメ

3.2.2.

ペナルティ法を用いた目的関数の構成

201811171 FUJIWARA - 11

最適化する力学的パラメータは，非負条件と等号制約を持つため，目的関数𝐽(𝑥)にペナルティ項

Step3：𝐽′(𝒙𝑟𝑒𝑓 ) < 𝐽′(𝒙1 )なら，式(62)を元に拡大点

𝑃(𝑥)に適当な実数𝜇をかけたものを加えたものを，

𝒙𝑒𝑥𝑝を計算する．

新たな目的関数とし，式(57)に示す．これを元問題(62)

𝒙𝑒𝑥𝑝 = (1 + 𝛽)𝒙𝐺 − 𝛽𝒙𝑛+1のペナルティ問題とよぶ．式(57)は制約なし最適化

𝐽′(𝒙𝑒𝑥𝑝 )の計算を行い，𝐽′(𝒙𝑒𝑥𝑝 ) < 𝐽′(𝒙𝑟𝑒𝑓 )を満たす問題となるため，制約なし最適化の解法を使える．

場合，𝒙𝑛+1 に𝒙𝑒𝑥𝑝 を代入する．そうでない場合，

𝐽 ′ (𝑥) = 𝐽(𝑥) + 𝜇𝑃(𝑥)(57)

𝒙𝑛+1 に𝒙𝑟𝑒𝑓 を代入する．

𝑁

𝑃(𝑥) = ∑ 𝑔𝑖 (𝑥)

(58)

𝑖=1

0,
𝑔𝑖 (𝑥) = {
∞,

𝑥∈𝑆
𝑥∉𝑆

(59)

Step4：𝐽′(𝒙𝑛 ) ≤ 𝐽′(𝒙𝑟𝑒𝑓 ) < 𝐽′(𝒙𝑛+1 )なら，式(63)を元に外部収縮点𝒙𝑜𝑐 を計算する．
(63)

𝒙𝑜𝑐 = (1 + 𝛾)𝒙𝐺 − 𝛾𝒙𝑛+1

𝜇を徐々に大きくすると，式(59)の最適解はペナル

𝐽′(𝒙𝑜𝑐 )の計算を行い，𝐽′(𝒙𝑜𝑐 ) < 𝐽′(𝒙𝑟𝑒𝑓 )を満たす場

ティ項を小さく保とうとし，𝑥 ∈ 𝑆を満たそうとす

合，𝒙𝑛+1 に𝒙𝑜𝑐 を代入する．そうでない場合，Step6

る．十分に大きい𝜇に対する最適解が得られれば，に進む．

その最適解は元問題の最適解に近いと考えられる．
本研究では𝜇の初期値を 1.00×105 とし，繰り返し計

Step5：𝐽′(𝒙𝑛+1 ) ≤ 𝐽′(𝒙𝑟𝑒𝑓 )なら，式(64)を元に内部収算ごとに 1.00×105 倍し，2 回で終了する．

縮点𝑥𝑖𝑐 を計算する．
(64)

𝒙𝑖𝑐 = (1 − 𝛾)𝒙𝐺 + 𝛾𝒙𝑛+1
3.2.3.

Nelder-Mead 法

𝐽′(𝒙𝑖𝑐 )の計算を行い，𝐽′(𝒙𝑖𝑐 ) < 𝐽′(𝒙𝑟𝑒𝑓 )を満たす場

Nelder-Mead 法は解を含む𝑛次元解空間を𝑛 + 1個の頂点からなる単体で囲い，反転・拡大・縮小の 3

合，𝒙𝑛+1 に𝒙𝑖𝑐 を代入する．そうでない場合は Step6に進む．

種類を繰り返しながら解空間を探索することで最適化問題の近似解を発見することのできる手法でStep6：2 ≤ 𝑖 ≤ 𝑛 + 1に対して，式(65)を元に次の収ある．以下にそのアルゴリズムを示す．

縮点を計算する．
𝒙′𝑖 = 𝒙1 + 𝛿(𝒙𝑖 − 𝒙1 )(65)

Step 0：ランダムに𝑛 + 1点の変数𝒙𝑖 を生成する．

そして，𝒙′𝑖 を新たな𝒙𝑖 に置き換えて Step1 へ戻る．

Step1：変数𝑥𝑖 それぞれにおいて目的関数の値𝐽′(𝒙)本研究では，{𝛼, 𝛽, 𝛾, 𝛿} = {1,1 + 𝑛 , 0.75 − 𝑛 , 1 −を計算する．得られた値を式(60)の条件を満たすよ

1
𝑛

0.5

}とした．また，値の更新が十分小さくなるもしく

は，100 回𝒙𝑖 の更新を行ったときにアルゴリズムをうに昇順に並べる．
𝐽′(𝒙1 ) ≤ 𝐽′(𝒙2 ) ≤ ⋯ ≤ 𝐽′(𝒙𝑛 ) ≤ 𝐽′(𝒙𝑛+1 )

2

(60)

終了する．

その後，𝒙𝑛+1 を除き重心𝒙𝐺 を計算する．
3.3.
Step2：式(61)を元に鏡像点𝒙𝑟𝑒𝑓 を計算する．

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

作成したシミュレーションモデルの妥当性を確(61)認するために，表２を入力諸元とする単純梁につ𝐽′(𝒙𝑟𝑒𝑓 )の計算を行い，𝐽′(𝒙1 ) ≤ 𝐽′(𝒙𝑟𝑒𝑓 ) < 𝐽′(𝒙𝑛 )をいて 数値解を 求めた． 図 2 の 橋梁の中 央部に満たす場合，𝒙𝑛+1 に𝒙𝑟𝑒𝑓 を代入する．

100[kN]の静的荷重を加え，この時の中央部の変位

𝒙𝑟𝑒𝑓 = (1 + 𝛼)𝒙𝐺 − 𝛼𝒙𝑛+1

201811171 FUJIWARA - 12

Disp.[ ]数値解

解析解Time.[ ]図 9 静的荷重を加えた橋梁モデルのたわみの結果を図９に示す．梁の平衡状態時の変位は究では２車両であるため，各車両の前輪と後輪で3.605×10-3[m]と解析解と一致した．この時，中央部推定される路面凹凸は合計４種類である．平均解析解を求める数式は次式となる．

𝑅̅(𝑥)は，

𝑤𝑚𝑎𝑥 =

𝑃𝐿3
48𝐸𝐼

1 (1)
(1)
(2)
(2)
𝑅̅(𝑥) = (𝑅𝐹 + 𝑅𝑅 + 𝑅𝐹 + 𝑅𝑅 )
4

(66)

(67)

ここで，𝑤𝑚𝑎𝑥 は中央部のたわみ，𝑃は中央部に加えである．(1)，(2)は車両番号，𝐹𝑅はそれぞれ前輪，

た静的荷重，Lは梁の長さ，𝐸𝐼は橋梁の曲げ剛性で後輪を表す．したがって目的関数は，ある．複数車両における VBI システム同定
本研究では，同じ橋梁上を複数台の車両が走行

(1)

2

(1)

2

2

(2)
(2)
+|𝑅𝐹 (𝑥) − 𝑅̅ (𝑥)| + |𝑅𝑅 (𝑥) − 𝑅̅(𝑥)| }した際の車両振動データから，車両・橋梁の力学的パラメータと路面凹凸を推定する手法を提案する．

2

𝐽 = ∑{|𝑅𝐹 (𝑥) − 𝑅̅(𝑥)| + |𝑅𝑅 (𝑥) − 𝑅̅(𝑥)|

(68)

となる．

複数車両での同定フローを図 10 に示す．
まず，数値実験で車両振動および位置データを

4. 提案手法の適用検証得る．次に，VBI システム同定によって，車両・橋梁の力学的パラメータと路面凹凸を推定する．パ

4.1.

概要

ラメータを更新する際，橋梁の力学的パラメータ

複数車両の VBI システム同定における車両・橋はそれぞれの車両で推定された値を平均し，推定梁の力学的パラメータと路面凹凸の推定結果を異値とする．

なるノイズの条件下において確認する．
本研究でのノイズは，エンジン振動や観測ノイ

3.3.1.

複数車両における目的関数

ズを想定しているので，全てのセンサに同じ大き複数車両が同じ橋梁上を直進していると仮定しさのノイズが乗ることが想定する．そこで，1.0%のた場合，各車両は同じ路面上を走行するため，推定

ノイズ(観測データの実効値に対するホワイトノ

される路面凹凸は一致するはずである．そこで，そイズ)を加えた場合をノイズなしの場合と比較す

れぞれの車両で推定された路面凹凸の平均を求め，る．

各路面凹凸との二乗誤差を目的関数とする．本研201811171 FUJIWARA - 13

Measure Vehicle Vibration:
and Position:,,

Randomly Assume Parameters:
,
,
,
,
,
,

Calculate Input Profile:

,

,

,

,

Randomly Assume Parameters:
,
,
,
,
,
,

Calculate Input Profile:,,

Randomly Assume Parameters:
,
,
,
,
,
,

Calculate Input Profile:,

Vehicle (1)Vehicle (1)Vehicle (1)(Kalman Filter)(Kalman Filter)(Kalman Filter)

Vehicle (2)Vehicle (2)Vehicle (2)(Kalman Filter)(Kalman Filter)(Kalman Filter)

Calculate Bridge Profile:,Calculate Bridge Profile:,Calculate Bridge Profile:,

BridgeBridgeBridge(Newmark-beta)(Newmark-beta)(Newmark-beta)

Calculate Road Profile:Calculate Road Profile:Calculate Road Profile:

Transformto:Calculate Objective Function:

Transform

to:Transform

Calculate Objective Function:

to

,

,

:

Calculate Objective Function:

Update -th Parameters based on Nelder-Mead method:
,
,
,
,
,
,
,
,

Repeat 100 times
Calculate Objective Function:

図 10 Nelder-Mead 法を用いたパラメータ更新

車両・橋梁の力学的パラメータの推定結果

車両・橋梁の力学的パラメータの推定精度が悪い．

表１と表２の車両・橋梁モデルを使用し，得られ計測データに含まれるノイズの影響で，推定パラ4.2.

た車両振動から，Nelder-Mead 法を用いて車両・橋メータに誤差が蓄積することが確認できた．

梁パラメータの推定を行った．図 11 に橋梁が健全

10t 車両ではノイズがない場合でも車両・橋梁の時の 10t 車両と 18t 車両におけるノイズなしと 1.0%

力学的パラメータは正解値に一致しない．今回採

ノイズで比較した車両・橋梁の力学的パラメータ

用した最適化手法である Nelder-Mead 法では正しの結果を示す．また図 12 は橋梁が損傷時の結果をい解に到達できなかったと考えられる．一方，18t示す．数値実験は 28 回行い，得られたデータを正

車両ではノイズがない場合，車両パラメータは正規化し，28 回の平均をとった．パラメータの初期解値にほぼ一致する．つまり，車体の重心位置と車探索空間の上下限は，車軸間距離𝑑1 は正解値に対重が推定精度に影響する．

して 10%から 90％までとし，その他のパラメータは 80％から 120％とした．
ノイズありの場合のほうがノイズなしに比べて全体として 10t 車両のほうが，18t 車両に比べて推定精度が悪い．特に𝑑1 の推定精度が悪いことから，車体重心位置が推定精度に影響していると考201811171 FUJIWARA - 14

10

10

×18
10

10
×10

×18
18

10

×18

正解値

𝑑1 𝑐𝑠1 𝑐𝑠2 𝑘𝑠2 𝑘𝑠1 𝑚𝑢1 𝑚𝑢2 𝜌𝐴

𝑑1 𝑐𝑠1 𝑐𝑠2 𝑘𝑠2 𝑘𝑠1 𝑚𝑢1 𝑚𝑢2 𝜌𝐴

損傷無し

𝑑1 𝑐𝑠1 𝑐𝑠2 𝑘𝑠2 𝑘𝑠1 𝑚𝑢1 𝑚𝑢2 𝜌𝐴

𝑑1 𝑐𝑠1 𝑐𝑠2 𝑘𝑠2 𝑘𝑠1 𝑚𝑢1 𝑚𝑢2 𝜌𝐴

損傷無し

損傷無し

損傷無し

正解値

𝐸𝐼1 𝐸𝐼2 𝐸𝐼3 𝐸𝐼4 𝐸𝐼5 𝐸𝐼 𝐸𝐼 𝑘𝑢1 𝑘𝑢2

(a)10

ノイズなし

𝐸𝐼1 𝐸𝐼2 𝐸𝐼3 𝐸𝐼4 𝐸𝐼5 𝐸𝐼 𝐸𝐼 𝑘𝑢1 𝑘𝑢2

(b)10

ノイズあり

𝐸𝐼1 𝐸𝐼2 𝐸𝐼3 𝐸𝐼4 𝐸𝐼5 𝐸𝐼 𝐸𝐼 𝑘𝑢1 𝑘𝑢2(c)18

𝐸𝐼1 𝐸𝐼2 𝐸𝐼3 𝐸𝐼4 𝐸𝐼5 𝐸𝐼 𝐸𝐼 𝑘𝑢1 𝑘𝑢2

ノイズなし

(d)18

ノイズあり

図 11 10t 車両と 18t 車両の推定パラメータ(橋梁：健全)

10
10

×18
10

10
×10

×18
18

10

×18

正解値

𝑑1 𝑐𝑠1 𝑐𝑠2 𝑘𝑠2 𝑘𝑠1 𝑚𝑢1 𝑚𝑢2 𝜌𝐴

損傷有り

𝑑1 𝑐𝑠1 𝑐𝑠2 𝑘𝑠2 𝑘𝑠1 𝑚𝑢1 𝑚𝑢2 𝜌𝐴

損傷有り

𝑑1 𝑐𝑠1 𝑐𝑠2 𝑘𝑠2 𝑘𝑠1 𝑚𝑢1 𝑚𝑢2 𝜌𝐴

𝑑1 𝑐𝑠1 𝑐𝑠2 𝑘𝑠2 𝑘𝑠1 𝑚𝑢1 𝑚𝑢2 𝜌𝐴

損傷有り

損傷有り

正解値

𝐸𝐼1 𝐸𝐼2 𝐸𝐼3 𝐸𝐼4 𝐸𝐼5 𝐸𝐼 𝐸𝐼 𝑘𝑢1 𝑘𝑢2

(a)10

ノイズなし

𝐸𝐼1 𝐸𝐼2 𝐸𝐼3 𝐸𝐼4 𝐸𝐼5 𝐸𝐼 𝐸𝐼 𝑘𝑢1 𝑘𝑢2

(b)10

ノイズあり

𝐸𝐼1 𝐸𝐼2 𝐸𝐼3 𝐸𝐼4 𝐸𝐼5 𝐸𝐼 𝐸𝐼 𝑘𝑢1 𝑘𝑢2(c)18

𝐸𝐼1 𝐸𝐼2 𝐸𝐼3 𝐸𝐼4 𝐸𝐼5 𝐸𝐼 𝐸𝐼 𝑘𝑢1 𝑘𝑢2

ノイズなし

(d)18

ノイズあり

図 12 10t 車両と 18t 車両の推定パラメータ(橋梁：損傷)えられる．また，図５より，18t 車両に比べて 10t 車

棒グラフ)のほうが，誤差が大きいものもあり，車両の橋梁振動が卓越していないことが確認できる．

両 2 台において精度が向上しているとはいえない．

これは車重の影響であり，橋梁パラメータにおい車両を増やすことでノイズ特性を明らかにできれて推定精度低下の原因となる．

ば，ノイズに対してより頑健なモデルにできる可バネ下剛性𝑘𝑢𝑖 や橋梁パラメータ𝐸𝐼𝑖 において車
両 1 台(青色の棒グラフ)より，車両 2 台(赤と黃の能性がある．
損傷時の橋梁の値は，健全時の橋梁の値で正規

10t

前輪

ノ ズ無し

10t

ノ ズ1.0%

後輪
ノ ズ無し

18t

ノ ズ1.0%

Disp.[ ]10-10tノ ズ無し

10-10tノ ズ1.0%

Disp.[ ]10-18tノ ズ無し

10-18tノ ズ1.0%

10-18tノ ズ無し

10-18tノ ズ1.0%18-18tノ ズ無し

18-18tノ ズ1.0%

Disp.[ ]18t

Disp.[ ]正解値

Disp.[ ]Disp.[ ]201811171 FUJIWARA - 15

Traveling position[ ]Traveling position[ ]

ノ ズ無し

10t

ノ ズ1.0%

ノ ズ無し

18t

ノ ズ1.0%

Disp.[ ]10-10tノ ズ無し

10-10tノ ズ1.0%

10-18tノ ズ無し

10-18tノ ズ1.0%

10-18tノ ズ無し

10-18tノ ズ1.0%18-18tノ ズ無し

18-18tノ ズ1.0%

Disp.[ ]18t

Disp.[ ]前輪 後輪

Disp.[ ]正解値

10tDisp.[ ]Disp.[ ]図 13 路面凹凸推定（橋梁：健全）

Traveling position[ ]Traveling position[ ]図 14 路面凹凸推定（橋梁：損傷）

201811171 FUJIWARA - 16化した．しかし，𝐸𝐼4 の誤差が大きくならず，橋梁

解には到達できない．したがって，計算効率の損傷の検知にはいたらなかった．

よく最適解に到達できる最適化手法が求めら
れる．

4.3.

路面凹凸推定(3)

次に 10t 車両と 18t 車両の推定された路面凹凸をじる可能性があり，明確な原因を追求する必正解値と比較し，橋梁が健全時の結果を図 13 に，損傷時の結果を図 14 に示す．

車両モデルの違いによって推定精度に差が生

要がある．
(4)

2 車両では 1 車両のときと比べて推定精度はノイズがある場合はノイズがない場合に比べて向上しなかった．さらに車両を増やすことで推定精度が悪い．計測データに含まれるノイズのノイズ特性の影響を明らかにすることができ影響で，推定パラメータに誤差が蓄積することがれば推定精度が向上する可能性がある．

確認できた．
全体的に 10t 車両のほうが 18t 車両に比べて推定精度が悪く，車両・橋梁の力学的パラメータと同じ

参考文献
[1]

結果となった．したがって，車体の重心位置と車重が推定精度に影響すると考えられる．

国土交通省：道路の老朽化対策，老朽化対策
の取り組み

[2]

Y.B. Yang, C.W. Lin and J.D. Yau：Extracting
bridge frequencies from the dynamic response ofa passing vehicle, Journal of Sound and5. まとめと今後の課題
5.1.

まとめ

Vibration, 272(3-5),pp.471-493,2004.
[3]

本研究では，複数車両が同じ橋梁を走行したと

車両応答に基づく橋梁のモード形状推定法，構造工学論文集 A1，67，pp.242-257,2011.

きの車両・橋梁の力学的（質量・減衰・剛性）パラ
メータと路面凹凸の推定精度の変化を検討した．

山本亨輔，大島義信，杉浦邦征，川野広隆：[4]

長山智則，趙博宇，薛凱：走行時の車体振動

車体中心に重心があり，ノイズがない場合は車両・を利用したハーフカーモデルの同定と路面縦

橋梁パラメータと路面凹凸を精度よく推定できた．

断形状の推定．
土木学会論文集 E1，
75(1)，
pp.1-

しかし，ノイズがある場合は推定パラメータの誤

16,2019.

差の蓄積により推定値が正解値に一致せず，複数[5]

村上翔：粒子群最適化に基づく複数車両の振車両にしても改善が見られなかった．車両モデル

動データを用いた車両・橋梁・路面のパラメの違いによっても推定精度に差が見られた．

ータ同定，筑波大学，学士論文，2019．

Nelder-Mead 法は計算効率が良く，それなりに推定[6]

村上翔：車両振動へのカルマンフィルタ適用
による VBI システム同定の可能性に関する数

精度も良いが，正しい解には到達できなかった．

値的検討，筑波大学，修士論文，2021
5.2.
(1)

(2)

今後の課題

[7]

Ahmed EI-Hattab, Nasim UDDIN, Eugene

計測データに含まれるノイズはパラメータ推

Obrien : Drive-by bridge damage detection using

定精度を低下させる．ノイズの影響による誤apparent profile, Birmingham Alabama USA,

差を小さくすることが必要とされる．

UTC Conference for the Southeastern Region, 14-

Nelder-Mead 法は，計算効率は良いが，正しい

15 December 2015

201811171 FUJIWARA - 17

謝辞
本研究を進めるにあたり，多大な助言をいただいた指導教員の山本亨輔先生と松島亘志先生に深く感謝いたします．また，ゼミ等で助言をいただいた研究室の先輩方，日々の生活を支えてくれた両親に感謝いたします．