一般財団法人 新潟県建設技術センター
平成30年度 研究助成事業報告書

打音試験により得られた打音応答波形に基づく空洞位置同定手法の高精度化に関する事業

平成31年3月

長岡技術科学大学
長岡工業高等専門学校

機械創造工学専攻
環境都市工学科

倉橋村上貴彦祐貴（研究代表者）（研究分担者）

長岡工業高等専門学校
長岡工業高等専門学校
長岡工業高等専門学校
長岡技術科学大学

機械工学科
機械工学科
電子制御工学科
機械創造工学専攻

池田井山
外山
井原

富士雄
徹郎
茂浩
郁夫

（研究分担者）（研究分担者）（研究分担者）（研究分担者）1

1. はじめに
従来より，コンクリート内における空洞の状態について計測値を用いて非破壊的にまた数値的に評価する方法について検討が行われている．北原ら 1)-3)や廣瀬ら 4)-8)は，コンクリート表面に超音波を送信し，測定した散乱波形から散乱体の形状や物性を決定する方法（線形化逆散乱解析法）について検討を行っている．これらの研究では，二次元問題 1)，三次元問題 2),4)，非均質材料内の欠陥形状を再構成する問題 3)，複数の鉄筋の問題 6)，欠陥形状をリアルタイムで画像化するために超音波リニアアレイ探触子が有する高速画像処理の技術と

図 1 自己組織化マップによる空洞サイズ・位置逆散乱イメージング手法を合わせた方法に関するの分類の一例

研究 7),8)等が行われている．線形化逆散乱解析法とは異なる方法として，西村ら 5)による開口合成法

観測点による欠陥形状のイメージング手法に関する研究や，超音波の計測波形に以外にレーザ干渉計によ打撃点

り対象表面の粒子速度を計測し，その計測値を用いてコンクリート表面におけるクラックの深さおよび角度を決定する方法についても検討を行っている 9),10)．

空洞境界Γcav

本研究では，上記の内容に基づき，計測値を用いた随伴変数法 11)による空洞位置・サイズの同定問題について検討を実施する．検討の導入として，

図 2 打音試験モデルにおける空洞位置同定問題

音圧の測定値を用いた空洞位置の同定問題に取りのイメージ図組んだ．逆解析における空洞位置の初期条件を決定するために，自己組織化マップにより解析事前2. 空洞位置同定に向けた定式化に空洞のサイズや位置の分類（図 1）を作成し，そ2.1 支配方程式の情報をもとにした空洞位置の同定解析を実施し本研究では，インパルスハンマによる打撃によた 12)．検討結果としては良好な結果が得られたが，り得られる加速度センサの応答を用いた空洞位置静穏な環境で実施する必要があるという制約があの同定方法を考えているため，三次元の振動方程り，試験を行う環境によっては，適切な音圧デー式を用いる．総和規約により表示した振動方程式タが取得できず，空洞位置の同定解析が困難になを式(1)～(3)に示す．

るものと想定される．そこで，本検討では加速度センサによる計測値を用いて，空洞の位置を同定

σ𝑖𝑗,𝑗 = 𝜌𝑢𝑖̈(1)

σ𝑖𝑗 = 𝐷𝑖𝑗𝑘𝑙 𝜖𝑘𝑙(2)

する方法について検討を実施し（図 2），定式化および同定解析の流れ，また解析によって得られた結果について報告する．
2

1ϵ𝑖𝑗 = 2 (𝑢𝑖,𝑗 + 𝑢𝑗,𝑖 )較をすることで，空洞位置の評価を行う．このこ(3)とに基づき，評価関数を式(10)のように定義する．
ここに，σij は各応力成分，ui は x,y,z 方向の変位，Dijkl は弾性係数行列，εij は各ひずみ成分を示す．有𝐽=限要素法による定式化に基づき，重み付き残差方

１２({𝑢i } − {𝑢iobs. })𝑇 [𝑄]({𝑢i } − {𝑢iobs. })(10)

程式を誘導し，四面体一次の補間関数を適用すると，式(4)に示す離散化された振動方程式が得られ式(9)において[Q]は重み対角行列，T は転地記号，る．

{ui}は振動解析により求まった変位，{uiobs.}は加速度センサにより測定された加速度応答を変位応答[𝑀e ]{𝑢𝑒𝑖̈ } + [𝐾e ]{𝑢ei} = {𝑓𝑒𝑖 }に変換した値を示す．本研究では式(9)を最小とす(4)

る空洞位置を求める問題を考える．式(9)が最小に添え字 e は各要素番号を示す．質量減衰また剛性なるということは，解析変位{ui}が測定変位{uiobs.}

減衰も考え，[Ce]= α[Me]+ β[Ke]（α，β は比例係数）に近づくということを表している．とし，式(4)に減衰項を導入すると式(5)に示す離散化式が得られる．

2.3 ラグランジュ関数の停留条件の導出
式(9)の最小化問題において，式(6)～式(9)が制約[𝑀e ]{𝑢𝑒𝑖̈ } + [𝐶𝑒 ]{𝑢𝑒𝑖̇ } + [𝐾e ]{𝑢ei } = {𝑓𝑒𝑖 }条件となる．そのため，随伴変数ベクトル{λ}を導(5)入し，評価関数を式(11)のように拡張する．
全要素を足し合わせた式として添え字 e を外し，𝐽∗ = 𝐽 + {𝜆𝑖 }𝑇 ([𝑀]{𝑢̈𝑖 } + [𝐶]{𝑢̇𝑖 } + [𝐾]{𝑢i } − {𝑓𝑖 })式(6)のように表す．

(11)[𝑀]{𝑢̈𝑖 } + [𝐶]{𝑢̇𝑖 } + [𝐾]{𝑢i } = {𝑓𝑖 }(6)式(11)はラグランジュ関数と呼ばれ，この関数の停また，式(6)における初期条件・境界条件を式(7)～留条件（δJ*=0）を満たすように，空洞位置の計算(9)のように定義する．

を行う．停留条件はラグランジュ関数の第一編分を計算することで求めることができ，結果として，{𝑢i (𝑡0)} = {𝑢̂𝑖 } {𝑢̇ i (𝑡0 )} = {𝑢̂𝑖̇ } in Ω

状態変数に対応した随伴変数に対する方程式（式(7)(12) ）が得られる．式(12)は随伴方程式と呼ばれて{𝑢i } = {𝑢̂i } on Γ1(8){𝑓i } = {𝑓̂i } on Γ2

(9)

いる．
[𝑀]𝑇 {𝜆𝑖̈ } − [𝐶]𝑇 {𝜆𝑖̇ } + [𝐾]𝑇 {𝜆i }+[𝑄]𝑇 ({𝑢i } − {𝑢iobs. }) = {0}(12)

ここに t0 は解析初期時刻，Ω は解析領域，Γ1，Γ2は第一種境界，第二種境界を示す．

また，式(13)～式(15)もラグランジュ関数の停留条件から得られ，式(13)～式(15)および振動解析によ

2.2 評価関数の定義り得られた変位応答{ui}を用いることで式(12)を加速度センサにより得られた加速度を変位応答解くことになる．状態変数に対して初期条件が与に変換し，振動解析により求まった変位応答と比えられることに対して，随伴変数に対しては式(13)3に示す終端条件が与えられるため，随伴方程式は4. 各反復回数における評価関数の差の絶対値が解析終端時刻 tf から逆時間方向に解かれることに許容判定値 ε 以内の場合，計算を終了する．そなる．また，状態変数に対応した形で境界条件もうでない場合は，次のステップに進む．

得られ，各境界 Γ1，Γ2 に対して変位，外力に対応

5. 随伴方程式に対する有限要素方定式（式(12)）

した随伴変数{λi}，{si}が{0}として与えられる．

を逆時間方向に計算する．
6. ラグランジュ関数の空洞位置に対する勾配ベ{𝜆i (𝑡𝑓 )} = {0} {𝜆̇i (𝑡𝑓 )} = {0} in Ωクトル{∂J*/∂xi}を計算する．

(13)

7. {xi}(l+1)={xi}(l)-{∂J*/∂xi}(l)により座標値を更新し，{𝜆i } = {0} on Γ1(14){𝑠i } = {0} on Γ2(15)ステップ 2 に戻る．
3. 空洞位置同定法に関する検証
本検討では図 3 に示す有限要素メッシュ・固定

2.4 状態方程式および随伴方程式に対する有限要境界条件を用いて空洞位置の同定解析を行う．観

素方程式の時間方向に対する離散化

測点は打撃点から R=200mm，および R=400mm の状態方程式に対する有限要素法方程式（式(6)）

ところに 90 度間隔に 4 点設定する(図 4)．観測値の時間方向の離散化にはニューマークの β 法を適は応答変位の値を用いる．その他の解析条件は表用する．ニューマーク β 法における速度，変位の1 に示す．本解析では，空洞の位置をコンクリート

時間進展式は式(16)，(17)のように与えられる．

∆𝑡{𝑢̇𝑖 }𝑛+1 = {𝑢̇𝑖 }𝑛 + 2 ({𝑢̈ 𝑖 }𝑛+1 + {𝑢̈ 𝑖 }𝑛 )(16){𝑢𝑖 }𝑛+1 = {𝑢𝑖 }𝑛 + ∆𝑡{𝑢𝑖̇ }𝑛1+Δ𝑡 2 (2 − 𝛽) {𝑢̈ 𝑖 }𝑛 + 𝛽Δ𝑡 2 {𝑢̈ 𝑖 }𝑛+1X-Y-Z Fix(17)

図 3 有限要素メッシュおよび固定境界条件
随伴方程式に対する有限要素法方程式（式(12)）に対しても同様にニューマーク β 法を適用する．
2.5 空洞位置同定解析の流れ
以上に誘導した方程式を用いた空洞位置の同定計算の流れを以下に示す．以下の反復計算においては，最急降下法を適用し，反復回数は(l)とする．
1. 同定解析に必要な計算条件の入力する．
2. 振動方程式に対する有限要素方定式（式(6)）を順時間方向に計算する．
3. 評価関数（式(10)）を計算する．

図 4 打撃点および観測点の位置について
4Table 1

表 1 解析条件いものの，正しい方向に空洞位置の同定解析を正Computational conditions

Computational area
(X, mm×Y, mm×Z, mm)

900×900×180

Diameter of the cavity φ, mm

400

Element typeLineartetrahedron

Average side length, mm50

Mass density ρ, kg/m32300

Young’s modulus E, MPa

25000

Poisson’s ratio

0.16

Damping coefficient κ0100

Damping coefficient κ110-7

Time increment Δt, ms0.0390625

Time step128

Step length α10-11

しい方向に行えていることがわかる．

図 5 観測点の位置の違いに対する表面から 0.06m の位置に仮定し，順解析を行い，評価関数の収束履歴の違いの比較

それぞれの観測点位置で得た変位応答値を観測値として使用し，空洞の初期深さを正解値からずらすことにより，同定解析を行う．
解析結果を以下に示す．図 5 に評価関数の履歴を示す．観測点の位置が打撃点からの距離
R=200mm の場合と R=400mm の場合について示す．
結果として，最初の評価関数の値を 1 とした正規化した値では，打撃点から遠くの距離に観測点を設定した場合（R=400mm）の方が評価関数の値を下げることができた．また，観測点位置を打撃点から等距離に設定しているため，代表点における

図 6 R=200mm の位置に観測点がある場合変位応答の比較を図 6，図 7 に示す．図 6 はの応答変位の比較

R=200mm の場合の結果であり，図 7 は R=400mmの場合の結果である．どちらの結果も空洞位置同定後（After estimation）は空洞位置道程前（Beforeestimation）に比べ，観測値（Target）に近づいていることを確認できる．また，図 8 に空洞位置の同定解析の履歴を示す．正解とする位置（Target）に対 し て ， 観 測 点 の 位 置 が R=200mm お よ び
R=400mm の両方の場合において，どちらも Targetの値に近づいていく様子がわかる．R=400mm の方が若干 Target の値に近い結果を示しているが，両ケースとも同等の結果を得ていると考えられる．

図 7 R=400mm の位置に観測点がある場合正解とする Target の値には完全には一致していなの応答変位の比較
5

図 8 正解値と空洞位置の同定履歴の比較
4. トンネリングアルゴリズムの導入による改良
次に，本応募の目的とした空洞位置が表面から深い箇所に位置する問題に対する検討結果を示す．

図 9 計算モデル図

表 2 に計算条件，図 9 に計算モデルを示す．正解とする空洞の深さは表面から 100mm の位置とし，反復計算における評価関数の履歴および空洞位初期空洞位置を 80mm に設定し，本研究により提置の同定履歴を図 10 に示す．2.5 節に説明した計案する方法による同定可能性について検証する．

算アルゴリズムにトンネリングアルゴリズムを加前章においても，完全に空洞位置を同定できなかえることにより，評価関数も概ねゼロに近づけるったことから，本検討では，2.5 節に説明した計算ことができ，空洞の位置も 100mm に同定可能となアルゴリズムに対して，文献 13)に示すトンネリンった．また，図 11 に示すように観測点においては，グアルゴリズムを加え，検討を実施した．

測定変位（空洞位置を 100mm に設定した場合のシミュレーション結果）と同定後の変位応答は概ね

表 2 解析条件
Number of elements (Target, Initial)等しくなっていることを確認できた．
9300, 9298

Number of nodes (Target, Initial)

2102, 2102

Time increment Δt, s3.90625×10-5

Number of time steps

256

Mass density ρ, kg/m32300

Young’s modulus E, MPa

35096

Damping parameters κ0, κ1100, 10-7

Convergence criterion ε10-5

Weight parameter Q1019

Step length of steepest decent α

1.00

Number of random seeds, irmax

200

Number of maximum iterations, itrmax

200

Step length of ARTA T106

Tunneling phase

図 10 評価関数の履歴と空洞位置の履歴6更なる同定精度の向上について検討を行った．その結果，同定精度の向上を確認できた．
付録 A：空洞なしの実試験体を用いた実際の変位応答と振動解析結果（順解析結果）の比較
本研究では，打撃試験により得られた応答変位を用いた空洞位置の同定問題について検討を実施したが，得られた測定値は，解析結果を人工的な測定値として利用している．実際の応答波形を用(a) 空洞位置が初期位置の場合いて空洞位置の同定をする場合は，先に物性パラメータの条件を合わせておく必要がある．本付録
では，空洞の無い無筋コンクリートを対象とし，シミュレーションによる応答波形の再現性について考察する．
試験体のモデルを図 A-1 に，また境界条件を図A-2 に示す．計算条件として，節点数は 25248，要素数は 5134，ポアソン比は 0.16，時間増分量 Δt は3.90625×10-5 （ 測 定 器 の サ ン プ リ ン グ レ ー ト が25,600Hz のため，その逆数を時間増分量としている．），時間ステップ数は 256 と設定する．また，(b) 空洞位置が同定後の場合

ヤング率は 35096MPa，密度は 1630kg/m3，減衰行

図 11 観測点における変位履歴の比較列を決めるための係数 α は 600，β は 0 と設定する．打撃点での打撃力の経時変化（測定値）を図5. おわりにA-3，図 A-3 の値を用いた振動解析結果（観測点に本研究では，打音試験におけるコンクリート表おける変位応答）を図 A-4，図 A-4 の結果につい面における応答変位を用いた空洞位置の同定解析て周波数解析した結果を図 A-5 に示す．図 A-4 でに対する検討を行った．同定解析においては，打は，実測の変位の値と比較をしており，実測の変音波形の伝播の解析には有限要素法およびニュー位値はハイパスフィルターを介して低周波側の成マーク β 法を適用し，空洞位置の同定における感分をカットして再現した応答波形を示している．

度を計算するために随伴変数法を適用した．

応答変位（解析結果）は，測定値に比べて解析開始

数値実験においては，空洞の同定位置について時におけるピーク変位は若干合わないが，周期的は若干，正解とする値には一致しない結果が得らには良好に一致する波形を示している．また，同れたが，観測点の位置を打撃点からの位置様の傾向は周波数解析の結果からも確認できる．

（R=200mm，R=400mm）と変えたモデルについて，実際の構造物に対して，本研究で検討を行った空両モデルとも初期に設定した空洞位置から正解と洞位置の同定計算を行う場合には，同等の試験片する値に近づいていることを確認できた．しかしに対して，本付録のような事前検証を行い，先にながら，完全に正解値を求めることができていな解析条件を確定した上で実施することが望ましい

いことから，トンネリングアルゴリズムを導入し，と考えられる．
7

Z

YX

(b) Case2
ZY

図 A-1 試験体モデルのメッシュ図

打撃点

観測点

Z

YX

Z

X

Y

200220X

100(c) Case3
図 A-3 打撃点における打撃力の経時変化
10040360400

図 A-2 境界条件について(a) Case1(a) Case1(b) Case2
8

謝辞
本論文の解析結果は，九州大学情報基盤研究開発センターの高性能アプリケーションサーバを使用し計算を行ったものである．計算機の使用について，センター関係者の方へ謝意を表す
参考文献
1) 中畑和之, 大西正浩, 北原道弘, 線形化逆散乱
解析の高速化と欠陥形状再生への適用，応用(c) Case3
図 A-4 観測点における変位の応答履歴

力学論文集, Vol. 5, pp.67-74, 2002.
2) 山田真幸，村上賢治，中畑和之，北原道弘, 円
柱側面からの計測波形による三次元欠陥像の再構成, 応用力学論文集, Vol. 6, pp.69-76, 2003.
3) 山田真幸, 大西正浩, 北原道弘, 線形化逆散乱
解析法による非均質材料内の欠陥形状再構成,
応用力学論文集, Vol. 7, pp.83-90, 2004.
4) 木本和志，廣瀬壮一, 二層体中の散乱体に対する三次元線形化逆散乱解析, 応用力学論文集,(a) Case1

Vol. 5, pp.75-80, 2002.
5) 木本和志, 松江剛士, 廣瀬壮一, 数値シミュレーションを用いた開口合成法の欠陥形状再構
成 能 に関 する 研 究, 応用 力 学論 文集 , Vol.7,pp.91-96, 2004.
6) 中畑和之, 藤澤伸匡, 廣瀬壮一,超音波を用いた鉄筋形状および付着評価のための一解析法,
応用力学論文集, Vol.8, pp.95-102, 2005.
7) 中畑和之, 廣瀬壮一, 超音波リニアアレイ探触子のための逆散乱イメージング法の開発,(b) Case2

応用力学論文集, Vol.9, pp.107-114, 2006.
8) 木本和志, 廣瀬壮一, SH 波超音波リニアアレ
イ探触子の時間領域モデリング, 応用力学論
文集, Vol.9, pp.15-121, 2006.
9) K.J.Langenberg, Introduction to the special issue on
inverse problem, Wave Motion, Vol.11, pp.99-112,
1989.
10) 吉川仁, 大田裕貴, 西村直志, レーザ計測によ(c) Case3
図 A-5 観測点における変位の応答履歴

る水平・鉛直超音波速度波形データを用いたクラック決定解析の精度検証, 応用力学論文
9集, Vol.7, pp.77-82, 2004.
11) T.Kurahashi and H.Oshita, Shape determination of
3-D reinforcement corrosion in concrete based onobserved

temperatureonconcretesurface,

Computers and Concrete, Vol.7, pp.63-81, 2010.
12) 松岡瑛喜，廣瀬陽一，倉橋貴彦，村上祐貴，外山茂浩，池田富士雄，井山徹郎，井原郁夫，打音検査による欠陥定量的評価の高精度化に対する随伴変数法の適用，材料，Vol.67，No.9，pp.869-876，2018．
13) 北山 哲士, 山崎 光悦, 一般化ランダム・トン
ネリング・アルゴリズムによる大域的最適化 :第 1 報,アルゴリズムの提示と数値計算例, 日本 機 械 学 会 論 文 集 A 編 , Vol.69, No.684,
pp.1250-1256, 2003.

10