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

微弱打撃により取得した加速度応答データを援用した
検査場所非依存型の空洞位置同定解析手法の構築に関する事業

2020年3月

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

機械創造工学専攻
環境都市工学科
機械工学科
機械工学科
電子制御工学科
機械創造工学専攻

1

倉橋
村上
池田
井山
外山
井原貴彦
祐貴
富士雄
徹郎
茂浩
郁夫

（研究代表者）（研究分担者）（研究分担者）（研究分担者）（研究分担者）（研究分担者）

1. はじめに2.1 3 次元領域に対する振動方程式

2012 年の笹子トンネルの崩落事故を契機に，打本研究では，式(1)～(3)に示す 3 次元の振動方程撃による非破壊試験に注目が集められている．本式を用いる．式(1)～(3)は総和規約により記載して事業関係者においても，打撃による非破壊検査にいる．

おいて定量的に検査することを目指し，有限要素σ𝑖𝑖𝑖𝑖,𝑗𝑗 = 𝜌𝜌𝑢𝑢̈𝚤𝚤

法と随伴変数法による打撃検査システムの提案を行った 1)．しかしながら，打撃の環境が変わる場合は十分な打撃ができず，微弱な応答しか得られない場合もある．例えば，高所における打撃試験で

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

ϵ𝑖𝑖𝑖𝑖 = �𝑢𝑢𝑖𝑖,𝑗𝑗 + 𝑢𝑢𝑗𝑗,𝑖𝑖 �

(3)

1

は，近年，ドローンによる非破壊検査も考えられ2ており，研究が進められている 2)．しかし，飛行中(1)

では打撃時における反力が取れないこともあり，ここに，σij は各応力成分，ui は x,y,z 方向の変位，打撃方法についても工夫が必要である．また，狭Dijkl は弾性係数行列，εij は各ひずみ成分を示す．四小区間における打撃検査においても，十分に打撃面体一次の補間関数を適用した有限要素法によりが実施できない場合もあり，上記のような高所や定式化を実施すると，式(4)に示す離散化された有狭小な領域における打撃試験では，弱い打撃力に限要素方程式が得られる．

おいて，空洞の有無や空洞の大きさ・深さ等を探[𝑀𝑀e ]{𝑢𝑢𝑒𝑒𝑒𝑒̈ } + [𝐾𝐾e ]{𝑢𝑢ei} = {𝑓𝑓𝑒𝑒𝑒𝑒 }査することになる．このようなことを踏まえ，昨年度の事業内容を拡張し，微弱打撃の条件下にお(4)いて，コンクリート表面における応答変位量を用

添え字 e は各要素番号を示す．質量減衰また剛性

いた空洞位置同定の可否について検証を行う．解減衰も考え，[Ce]= κ0[Me]+ κ1[Ke]（κ0，κ1 は減衰パ

析モデルとしては図 1 に示すモデルを用いる．

ラメータ）とし，式(4)に減衰項を導入すると式(5)に示す離散化式が得られる．

応答変位量の測定点

打撃点[𝑀𝑀e ]{𝑢𝑢𝑒𝑒𝑒𝑒̈ } + [𝐶𝐶𝑒𝑒 ]{𝑢𝑢𝑒𝑒𝑒𝑒̇ } + [𝐾𝐾e ]{𝑢𝑢ei } = {𝑓𝑓𝑒𝑒𝑒𝑒 }空洞(5)全要素において重ね合わせを行った式として添え供試体字 e を外し，式(6)のように表す．
支持具[𝑀𝑀]{𝑢𝑢̈𝚤𝚤 } + [𝐶𝐶]{𝑢𝑢̇𝚤𝚤 } + [𝐾𝐾]{𝑢𝑢i } = {𝑓𝑓𝑖𝑖 }図 1 打撃試験のモデル図(6)

また，式(6)における初期条件・境界条件を式(7)～2. 空洞位置同定手法について
本検討では，空洞位置の同定問題に際し，コン(9)のように定義する．

クリート表面における観測変位とシミュレーショ
{𝑢𝑢i (𝑡𝑡0)} = {𝑢𝑢�𝚤𝚤 } {𝑢𝑢̇ i (𝑡𝑡0 )} = �𝑢𝑢�𝚤𝚤̇ � in Ω

ンにより求まった変位の差の二乗積分値を評価関数として定義し，支配方程式(3 次元領域に対する{𝑢𝑢i } = {𝑢𝑢�i } on Γ1振動方程式)の制約条件のもと評価関数の最小化を行う．定式化の詳細は以下に整理する．
2(7)
(8)

{𝑓𝑓i } = �𝑓𝑓̂i � on Γ2(9)�𝜆𝜆i (𝑡𝑡𝑓𝑓 )� = {0} �𝜆𝜆̇i (𝑡𝑡𝑓𝑓 )� = {0} in Ω

ここに t0 は解析初期時刻，Ω は解析領域，Γ1，Γ2{𝜆𝜆i } = {0} on Γ1

は第一種境界，第二種境界を示す．

{𝑠𝑠i } = {0} on Γ22.2 評価関数の定義
加速度センサにより得られた加速度から変位応(13)(14)(15)答を求め，振動解析により求まった変位応答と比2.4 状態方程式および随伴方程式に対する有限要較をすることで，空洞位置の評価を行う．ここで，素方程式の時間方向に対する離散化状態方程式に対する有限要素法方程式（式(6)），

評価関数を式(10)のように定義する．

随伴方程式に対する有限要素法方程式（式(12)）の𝐽𝐽 =１２({𝑢𝑢i } − {𝑢𝑢iobs. })𝑇𝑇 [𝑄𝑄]({𝑢𝑢i } − {𝑢𝑢iobs. })時間方向の離散化にはニューマークの β 法を適用(10)する．

式(9)において[Q]は重み対角行列，T は転地記号，

2.5 空洞位置同定解析の流れ{ui}は振動解析により求まった変位，{uiobs.}は加速以上の定式化により誘導した方程式を用いて空度センサにより測定された加速度応答を変位応答

洞位置の同定計算を行う．本解析では局所最適化に変換した値を示す．本研究では式(10)を最小とす

アルゴリズムとして最急降下法を適用し，大域的る空洞位置を求める問題に対して解析を行う．

最適化アルゴリズムとしてランダム・トンネリング・アルゴリズム（以下，RTA）を適用する．まず2.3 ラグランジュ関数の停留条件の導出

最急降下法の計算の流れを以下に示す．反復回数式(10)の最小化問題において，式(6)～式(9)が制は(l)としている．

約条件となる．そのため，随伴変数ベクトル{λ}を導入し，評価関数を式(11)のように拡張する．

1. 同定解析に必要な計算条件の入力する．
2. 振動解析を実施する．
（式(6)）𝐽𝐽∗= 𝐽𝐽 + {𝜆𝜆𝑖𝑖 }𝑇𝑇 ([𝑀𝑀]{𝑢𝑢̈𝚤𝚤 } + [𝐶𝐶]{𝑢𝑢̇𝚤𝚤 } + [𝐾𝐾]{𝑢𝑢i } − {𝑓𝑓𝑖𝑖 })3. 評価関数（式(10)）を計算する．
4. 各反復回数における評価関数の差の絶対値が(11)許容判定値 ε 以内の場合，計算を終了する．そ式(11)はラグランジュ関数と呼ばれ，式(11)に対す

うでない場合は，次のステップに進む．

る停留条件より式(12)に示す随伴方程式が得られ

5. 随伴方程式の解析を実施する．
（式(12)）

る．

6. ラグランジュ関数の空洞位置に対する勾配ベクトル{∂J*/∂xi}を計算する．
[𝑀𝑀]𝑇𝑇 �𝜆𝜆𝚤𝚤̈ � − [𝐶𝐶]𝑇𝑇 �𝜆𝜆𝚤𝚤̇ � + [𝐾𝐾]𝑇𝑇 {𝜆𝜆i }+[𝑄𝑄]𝑇𝑇 ({𝑢𝑢i } − {𝑢𝑢iobs. }) = {0}7. {xi}(l+1)={xi}(l)-{∂J*/∂xi}(l)により座標値を更新し，ステップ 2 に戻る．

(12)

また，式(13)～式(15)もラグランジュ関数の停留条

本検討では，上記に説明した局所最適化のアルゴ件から得られ，振動解析により得られた変位応答リズムの計算後，文献 3)に示す RTA による最適化{ui}を用いることにより式(12)を解くことになる．

を行う．RTA は図 2 に示すように多峰性を有する
3

評価関数を対象とする場合は非常に有効な手法で900あり，局所最適解に収束した場合も，大域的最適900解を探索するフェーズに移行し，解析を実施する．

空洞y

図 2 RTA のイメージ図
xh

3.1 空洞位置の同定解析

z180x3. 数値実験による検討
数値実験における打撃位置を図 3，図 4（図 4 は
図 3 の赤枠拡大図）に示す．正解とする空洞深さ50(50)800をコンクリート上表面から 50mm と設定し，打撃荷重のピーク値を 2000N，1000N，500N と変え空図 3 打撃試験の打撃・観測点領域（赤枠部）洞位置の同定結果について考察を行う．振動解析における物性値，時間方向の離散化に対する解析
条件，空間方向の離散化に対する解析条件，空洞位置の同定解析に対する解析条件を表 1 に整理す
る．また，打撃による境界条件（時間応答）の一例を図 5 に示す．上記の解析条件に基づき，大域的最適化アルゴリズムである RTA に基づき空洞深さの同定解析を行う．今回の検討では，図 4 に示すように空洞の中心と，打撃点の位置をずらして設

図 4 打撃試験の打撃・観測位置（赤枠拡大図）定しており，空洞の平面位置が不明瞭な場合についても近傍を打撃することにより，空洞深さが同表１ 解析における条件定可能であるか否かについて検討を実施する．コ

総要素数2087

ンクリート上面において，空洞中心が総節点数9197(x,y)=(450mm,450mm)であることに対して，打撃点

タイムステップ数256の位置は(x,y)=(550mm,550mm)とする．空洞の初期時間刻み

39.0625

位置は 70.5mm として空洞位置同定解析を実施す

ヤング率35.096

る．

ポアソン比0.16

質量密度2300

減衰パラメータ

90, 10-6評価関数の重みQ

41090.380

図 5 打撃荷重応答（Case1：ピーク値 2000N）
0.381

同定解析を行った結果を以下に整理する．図 60.380は評価関数の履歴を示す．また図 7 に空洞深さの同定履歴を示す．結果より，図 6 において，評価

図 6 各打撃荷重における評価関数の比較関数が不連続的に変わっている箇所は RTA を適用することにより，局所最適解（空洞位置を誤って同定した場合の結果）を回避し，大域的最適解のフェーズへシフトし解析を実施していることを示しており，大域的最適解を求めるための方法論が必要であることを示している．また図 6 より，どの荷重状況下においても，正解とする空洞深さ

46.36 mm

50mm 近傍へ解（空洞位置）が収束していることが46.35 mm

46.34 mmわかる．
上記の解析結果の妥当性に関する考察として，各荷重下の解析において空洞の深さと評価関数の値と，空洞位置同定解析により得られた結果を同図上にプロットする．図 8～10 は，各荷重下にお図 7 各打撃荷重における空洞深さの同定履歴

いて空洞の深さと評価関数の値の関係を図化したものであり，赤線は空洞位置を少しずつ変え順解析を実施した際に得られた評価関数の値，青のプ
ロットは空洞位置同定解析（逆解析）により得られた結果を示している．図 8～10 はそれぞれ，荷重のピーク値が 2000N，1000N，500N の際の結果である．結果より 70.5mm（各図の右側 x=70.5mmのところ）からスタートし，除々に右側へ推移し，局所最適解へ収束した時点で，左側の凹部における探索へ移行していることを確認できる．以上に示したように，多峰性を有する評価関数を対象とする場合には，RTA 等の大域的最適化のアルゴリ

図 8 空洞深さと評価関数の関係ズムの適用が必須となる．

（荷重ピーク値 2000N の場合）
5

図 13 は評価関数の履歴を示している．局所最適化と大域的最適化のフェーズを繰り返しながら最終的に収束をしていることがわかる．また，形状同定後の結果を図 14 に示す．同定後の形状（青線）は完全に目的形状（赤線）一致はしてないものの，初期形状からスタートし，目的形状に近づいていることを確認できる．また各観測点における目的形状時における変位の応答波形（赤線）と同定後の形状における変位の応答波形（青線）の比較を
図 15 に示す．観測点全てにおいて良好な一致を示
図 9 空洞深さと評価関数の関係していることを確認できる．

（荷重ピーク値 1000N の場合）

(a) 鳥瞰図

図 10 空洞深さと評価関数の関係（荷重ピーク値 500N の場合）
3.2 空洞形状の同定解析
上記の内容をもとにして，形状の同定解析へ展開し検討を実施した内容を以下に示す．試験体モデルのサイズは 900mm×900mm であり，高さは180mm である．空洞は薄い円柱を仮定し，直径300mm，コンクリート表面から空洞中心までの垂直距離は 100mm としている．図 11 にターゲットとする空洞の形状（赤線）と，逆解析の際に設定する初期形状（青線）を示す．また図 12 にメッシュ図，固定境界条件，打撃点，観測点を示す．上記の(b) 平面図条件のもと空洞の平面形状を同定した結果を以下

図 11 空洞の目的形状（赤線）に整理する．

および初期形状（青線）
6

(b) 平面図
図 14 空洞の目的形状（赤線）
および同定後の形状（青線）
図 12 メッシュ図，固定境界条件，打撃点および観測点について(a) 平面図(図 12)における左下の観測点
図 13 評価関数の収束履歴(b) 平面図(図 12)における右下の観測点(a) 鳥瞰図7また，空洞位置の同定問題を空洞形状の同定問題に拡張し検討を行ったところ，形状の一致性はまだ十分ではないが，初期形状に比べて目的形状に近い同定結果を得ることができた．各観測点における応答波形を見ても，十分に一致していることから，本解析結果以上の結果（より目的形状に近い形状）を求めることは困難であると考えられる．そのため，より観測点を増やす等，目的形状をより正確に表すために必要な観測点の位置等につ(c) 平面図(図 12)における右上の観測点

いて将来的には検討を要すると考えられる．
謝辞
本報告は，2019 年度の一般財団法人 新潟県建設技術センターの「微弱打撃により取得した加速度応答データを援用した検査場所非依存型の空洞位置同定解析手法の構築に関する事業」による研究助成により実施したものである．本報告内における解析結果は九州大学情報基盤研究開発センターの高性能アプリケーションサーバを使用し計算を行ったものである．新潟県建設技術センターな(d) 平面図(図 12)における左上の観測点

らびに九州大学情報基盤研究開発センターの関係

図 15 各観測点における目的形状時と同定形状時者の方へ謝意を表す．

における応答変位量の比較
参考文献
4. おわりに

1) 吉原 健太，倉橋 貴彦，村上 祐貴，外山 茂浩，

本事業では，打撃試験における空洞位置の同定

池田 富士雄，井山 徹郎，井原 郁夫， 随伴

解析に対して RTA を適用し，荷重量を変えた場合変数法に基づくコンクリート構造内部におけについて検討を行った．各荷重の条件に対して数る欠陥深度の推定解析 （打音検査による観測値実験による比較を行ったところ，どの荷重条件変位量を想定した数値実験）, 日本機械学会論下においても，RTA による探索フェーズを介する文集, 第 85 巻，第 869 号 pp.1-13, 2019.

ことで，正解とする解（空洞位置）近傍の値を得る

2) DiaTec Corporation,ことができることを確認できた．実際の打撃検査

http://www.diatec-net.co.jp/service/list.html

においても，打撃箇所の中心直上部を打撃できる

3) 北山 哲士, 山崎 光悦, 一般化ランダム・トンとは限らず，本検討で実施したように，打撃箇所ネリング・アルゴリズムによる大域的最適化 :と空洞中心位置がずれた場合についても考える必第 1 報, 日本機械学会論文集 A 編, Vol.69,

要がある．この様な場合（局所最適解に落ち着い

No.684, pp.1250-1256, 2003.てしまう場合）には，少なくても本解析法（RTA）

4) 竹内則雄，寺田賢二郎，樫山和男，計算力学―が必須であると考えられる．

有限要素法の基礎，森北出版(株)，2003．
8