こうえいフォーラム第９号／２００１．１

トモグラフィ的手法を用いた屈折法弾性波探査解析法の開発
DEVELOPMENT OF A SEISMIC REFRACTION ANALYSIS USING THE TOMOGRAPHICMETHOD

小原大輔*・ジャナキラマン カリヤン クマール**
Daisuke OHARA and Janakiraman K. KUMAR

Inversion is bound to become the geophysical interpretation method in the future. For thisreason, research on seismic inversion was undertaken. Computer programs for forward traveltime reconstruction, and ray tracing, were developed in order to study inversion. Linear Traveltime Interpolation (LTI)method was applied for forward calculation. Conventional and latest
inversion techniques were explored, including Iterative Least Square Technique (ILST) andSimultaneous Iterative Reconstruction Technique (SIRT) and Genetic Algorithm(GA). Seismic
inversion was applied to various cases. This paper briefly presents the results of theinvestigation.
Key Words : seismic inversion, linear travel-time interpolation,
genetic algorithm

１．は じ め に
現在、屈折法弾性波探査の解析においては、萩原の解析法もしくはその拡張解析法（以下、はぎ取り法）が一般的の速度構造を推定する調査手法であるが、以下の点において屈折法と異なる。
①

メッシュに分割による地盤のモデル化を行い、数値に用いられている 。この解析法は、地盤を凹凸が無い一計算を行うので、地盤の複雑な速度変化を表現する次元の層構造（水平方向に速度変化がない）と仮定し、基ことが可能である。

1）

盤より上部層の弾性波速度構造を求める方法である。この②的に求め、逆解析手法を用いてその理論走時と観測解析法には下記のような問題点がある。
①

一次元構造を仮定した解析法のため、山岳地のよう走時の差が小さい、観測データを最もよく説明するな地形、地質構造が複雑な箇所では解析結果の信頼速度構造を逐次的に求めることが可能である。

性が低下する（正確に言えば適用できない）。
②

起振点・受振点間の走時や伝播経路（波線）を理論図式解法のため、データにばらつきがある場合、判定に個人差の生ずる可能性がある。

③

解析の際、計算機を用いることで膨大な量のデータを処理することが可能である。

このような特長を屈折法弾性波探査に応用することによ一方、近年の計算機の性能向上に伴って、屈折法弾性波り、はぎ取り法の問題点である、複雑に変化する地盤への探査の解析自動化や高精度化に関する研究が盛んに試みら適用性や図式解法の問題点が改善できることから、屈折法れている。これらの研究に共通することは、弾性波トモグ

弾性波探査解析技術の精度向上のため必要不可欠な技術とラフィ技術を屈折法弾性波探査に応用している点である。

考えられた。

この弾性波トモグラフィは屈折法弾性波探査と同様に地盤

本稿では、弾性波トモグラフィ技術を応用した屈折法弾性波探査解析の理論を説明した後、実データの適用例につ＊＊＊中央研究所 開発研究部
元中央研究所 開発研究部

いて報告する。

7

トモグラフィ的手法を用いた屈折法弾性波探査解析法の開発

２．屈折法弾性波探査の概要
屈折法弾性波探査は、ダム・トンネル等、土木構造物の基礎地盤調査に用いられる調査手法である。この調査は、地表面に起振点（ダイナマイトなどの振源）、受振点（受振器）を設置し、その起振点から受振点までの弾性波の到達時間（走時）を観測することによって、地盤の速度構造を推定するものである（図−１参照）。さらに、この結果

Yesに基づき、地質構造や地盤強度を判定し、構造物の設計条件の設定や施工上の問題点を抽出することができる。

No

起振点において励起された弾性波動は、地盤内を透過し伝播する。一般的に地盤の速度分布は一様でないため、弾性波動は屈折・反射をしながら伝播する。

図−２ 解析フローチャート

例えば、２種の媒質内を弾性波動が通過して行く場合、その入射角が臨界角に達すると、屈折波は媒質の境界面を（1）順解析手法伝播する。このような屈折波は臨界屈折波（Head wave）

1）概要と呼ばれる。屈折法弾性波探査は、一度潜り込んだ弾性波

一般的に、屈折現象はSnellの法則によって記述するこ動が境界面を伝播し戻ってくる屈折現象を利用して探査をとができ、当法則は従来の屈折法弾性波探査の解析に用い行う。

られている。しかし、当法則を順解析手法として、より複雑な速度変化を持つ地盤に適用する場合、起振点からの射出角を試行錯誤的に変化させ、受振点に至る波線を決めなければならないなど、計算に困難さが伴う欠点がある。
一方、最近では、Snellの法則を用いずに、より複雑に速度が変化する地盤に適用できる様々な順解析手法が開発されている。主な順解析手法を表−１に示す。
今回開発した解析システムでは、

図−１ 屈折法弾性波探査の概念図（物理探査適用の手引き、物理探査学会より引用）

・線形近似を用いた高精度手法であること。

・計算時間の短い手法であること。

から、Asakawa and Kawanaka（1993）による線形走時内挿法（Linear Traveltime Interpolation）を順解析手法として採３．屈折法弾性波探査の解析理論用した。

トモグラフィ技術を応用した屈折法弾性波探査解析においては、まず起振点から受振点に至る弾性波の伝播経路

表−１ 主な順解析手法一覧（波線）を再現し、その到達時間を計算によって求めることが必要である。
一般的にこれを順問題（Forward Problem）、この問題を解く過程を順解析と呼ぶ。 特に弾性波トモグ
ラフィでは「波線追跡計算」または「レイトレーシング（Raytracing）」とも呼ばれる。次に、計算にて求められた計算走時と観測された走時の差（残差）を小さくするように地盤内の弾性波伝播速度を修正し、最適な速度構造を求める。つまり、残差が収束するまで、「波線追跡」→「速度構造の修正」の計算を繰り返す。与えられた観測データ

2）線形走時内挿法から最適な速度構造を求めることを逆問題（Inverse Pro-

線形走時内挿法（以下、LTI法）では、図−３のように、blem）、これを解く過程を逆解析（Inversion）と呼ぶ。逆解セル（cell）と呼ばれる小片に分割し、さらにセル境界上析手法を取り入れた屈折法弾性波探査解析でも、このように節点（node）を設け、地盤のモデル化を行う。セルには、な手順で解析を行う。解析の一連の流れを図−２に示す。

一定の速度値を与え（実際には、速度の逆数であるスロー

8こうえいフォーラム第９号／２００１．１

波線を追跡する過程では、受振点から各セル境界での波線の通過点を計算する。基本式を以下に示す。
（2）
ここでrは、線分AB間における波線の通過点である。
3）計 算 手 順
当手法の計算手順は、下記のようなステップからなり、図−５に示した。まず、領域全体の節点での走時を求める
過程では、以下の計算手順を踏む。
（a）起振点（SP）のあるセル（以下、起振セル）に含まれる節点での走時を、セルの速度と起振点−節点
図−３ LTI法における地盤のモデル化概念図間の距離から求める。

ネス（slowness）を与える）、速度構造を表現する。

（b）起振セルと同列にあるセルに含まれる節点での走時

このLTI法は、以下の２つの計算過程からなる。

を（1）式を用いて求める。この際、節点での走時①が最小であるものを採用する。

ある起振点から発生した弾性波動の、領域全体の節（c）起振セルの列の横にあるセル列の節点の走時を（1）点での走時を計算する。
②

その走時を用いて、各受振点まで最短時間で到達す式を用いて計算する。

る波線を追跡し、決定する。
LTI法では、①、②の過程を、線形近似を用いることで、高精度に走時を求めることが出来る。LTI法の概念図を図−４に示す。
図のように、すでに走時が判明している２節点A・B間の走時を線形近似し、隣の節点Dの走時を求める。図内の矢印は波線を示している。点Dの走時TDを求める基本式は以下のように記述される。
（1）
ここで、TAは点Aでの走時、ΔTは点A、Ｂの走時差、sは走時計算過程

各セルでのスローネス（速度の逆数）を表す。

波線追跡過程図−５ LTI法の計算手順
図−４ LTI法の概念図（Asakawa and Kawanaka（1993）より）

9トモグラフィ的手法を用いた屈折法弾性波探査解析法の開発（d）計算した列にある節点の中で最小走時を持つ節点をって計算を展開する。
（e）最終的に右辺上にある節点の走時を計算する。

Elevation(m)選び出す。そこから上辺もしくは下辺の節点に向か

0.0
-10.0
-20.00.5km/s-30.00.2km/s

0.5km/s-40.0-50.00.0（f）これら一連の過程を、領域全体の節点での走時が求

20.0

40.060.0

80.0

100.0120.0

140.0160.0

180.0

200.0

Distance(m)まるまで繰り返す
特に、（d）は臨界屈折波の挙動を計算する過程であり、図−６ 想定した理論速度構造モデル重要である。この計算が終わった後、各受振点に最小走時Elevation(m)で到達する波線を求める。この過程では、以下の計算手順を踏む。
波線を追跡する過程では、（a）受振点（RP）のあるセル上の節点から受振点までDistance(m)の走時を計算し、走時が最小である節点を選び出す。
（b）その節点の両隣のセグメントで最小走時を与えるを、（2）式を用いて求める。
（c）その点を新たな受振点とし、各セル境界での最適な波線通過点を求めていく。最終的に起振点（SP）に到達したら終了とする。
4）数値実験結果
LTI法の精度を検証するため、数値実験を行った。ここでは、２種類の速度から構成される、鉛直に延びる低速度帯モデルを想定した。この理論速度構造を図−６に示す。
当構造においては、Snellの法則を用いて、理論的に走時を得ることが可能であり、LTI法を用いた順解析計算による走時と、Snellの法則を用いた走時とを比較した。そ

図−７ 数値実験結果の結果を、図−７に示す。
図−７のように、Snellの法則とLTI法を用いた走時を比ここで、Skは、k番目のセルのスローネス、likは番目の波較すると、各受振点での平均走時残差は、0.03msであり、線がk番目のセルを横切る長さである。また、mは全セルLTI法により精度よく理論走時が求められることがわか数、nは全観測データ数である。

る。

（3）式の解法は数多くあるが、弾性波トモグラフィの解析には反復最小二乗法（Iterative Least Square Technique、（2）逆解析手法
1）概以下ILST法）や同時反復再構成法（Simultaneous Iterative要

Reconstruction Technique、以下SIRT法）がよく用いられて逆解析の目的は、観測走時を説明できる地下速度構造をいる。また、最近では遺伝的アルゴリズムを適用した研究推定することである。しかし、この問題は高次の非線形問

事例7）も報告されている。

題であり、容易に解くことは出来ない。そこで、v(x,y)と2）反復最小自乗法

して連続的に分布する速度構造を、一定速度を持つセルにILST法は、最小二乗法的にマトリクス計算を行う手法て離散的に表現する。そのため、観測走時と地下速度構造であり、一回の反復計算に伴う速度の修正量ΔSは、下式の関係は、（3）式のような線形方程式にて記述できる。

にて記述される。
（4）（3）

ここで、Bは計算走時と観測走時の差であり、残差を表している。ILST法は解の収束が早いが、問題の性質により計算が不安定になる特徴がある。

10

こうえいフォーラム第９号／２００１．１

3）同時反復再構成法

Select a sub-population for offspring production

SIRT法はマトリクス計算を行わず、残差や各セルを透For example Sub-population = {x3,x2}

過する波線の長さより、直接、速度の修正量を求める。速Produce offsprings by Recombination
For example = {x5,x6}

度の修正量ΔSは下式にて記述される。

Produce offsprings by Mutation
For example = {x7}
Evaluate Offspring Population（5）x5= ; x6= , x7= （図中●で示す。
）
Select the survivors from actual fitness
Best 4 Individuals = { }ここで、OSiはi番目の波線の残差を全波線長で割ったものEndである。このSIRT法は計算を安定的に進めることが出来るが、解の収束が遅いという特徴がある。
上記の解析法においては、はじめに初期速度構造（初期モデル）を与え、その構造に対して順解析を行い、計算走時ならびに各セルを透過する波線長を求める。次に計算走時と観測走時の差が小さくなるように、（4）または（5）式を用いて速度の修正量を求め、速度構造の修正を行う。
求められた速度構造に対して再度順解析を行い、さらに速度の修正を行う。最終的に残差が許容出来るほど小さくなるまで、反復計算を繰り返し、速度構造を推定する。
4）遺伝的アルゴリズム
地下速度構造を求める逆解析は、非線形の最適化問題である。このため、前述のILST法やSIRT法を用いた場合、

図−８ GAの最小解探索概念図

初期モデルの与え方によっては局所解（local optima）に陥いり、真の最適解（global optima）に到達しないことがあ以上の過程を繰り返すことにより、局所解に陥ることなく、最小解を探索する。

る。
一方、このような最適化問題に遺伝的アルゴリズム

GAでは、前述したように子孫を発生させる方法として、（Genetic Alogorithm、以下GA）を適用することも可能であ遺伝子組み替え（Recombination）と突然変異（Mutation）る。GAは生物の進化過程を模した計算手法であり、理論の２手法を用いている。遺伝子組み替えにて子孫｛x5, x6｝的にどのような初期モデルを設定しても、局所解に陥るこを発生させる方法を以下に示す。となしに真の最適解を探索できる。このGAについても逆解析手法として検討中であり、以下にGAによる最適解のDecimal 10Bit Binary form

探索手順、とGAによる速度構造推定例（水平２層構造12x3 =10=[1010101010]セル）を示す。

x2 =0=[0111000111]

X3 and X2 were crossed at randomly selected point.
X5 = [1111000111]= 22.53; Y = -11.69

・GAによる最適解の探索手順
図−８に示した非線形関数y=xsin(x)より、最小解（図中●で示す。
）を求めるGAの計算手順を説明する。
STEP-1 : Initialize a random population of individuals.

X6 =[0010101010]= -12.5; Y = -0.5突然変異は、個体をランダムに変化させる方法である。
突然変異にて子孫{x7}を発生させる方法を以下に示す。

STEP-2 : Evaluate fitness of all the initial individuals ofpopulation
1st best = x3; 2nd best = x2, 3rd best = x1, 4th best = x4
（x1〜x4を図中●で示す。
）

X2 =0

=[0111000111]

X2 was mutated at 4th binary point.
X7 =[0110000111]= 2.8; Y = 0.94

STEP-3 : While NOT DONE

11トモグラフィ的手法を用いた屈折法弾性波探査解析法の開発

・GAによる速度構造推定例
GAの適用例として、水平２層構造（セル数12）にて数値実験を行った。観測データは、その設定した水平２層構造より求めたものとした。ここで、速度構造モデルを定式化すると以下になる。
（6）（7）
ここで、f(x−)は観測値と計算値の差、x−は速度分布、viはi番目のセルの速度である。図−９の上に正しい理論速度構造、図−10 速度分布図、波線追跡図の出力例

下にGAを用いて得られた解析結果を示す。

５．山岳トンネル調査への適用例
1.00km/s

トモグラフィ的手法を用いた屈折法弾性波探査解析の妥当性を検討するため、山岳トンネル調査地にて得られたデータを用いてはぎ取り法と当解析法を適用した場合の解析精度の差違について検討を行った。
（1）はぎ取り法解析結果の検討
調査地は延長約500mの山岳トンネルが計画されている地点であり、地山状況を把握するため、屈折法弾性波探査

図−９ GAを用いた逆解析結果が実施されている。特に、トンネル線形に沿って実施されたのがL-1測線（測線長650m）であり、はぎ取り法によっ４．開発した屈折法弾性波探査解析プログラムて図−11のような速度構造が推定されている。特に図示し今回の開発では、３章で述べた順解析、および逆解析手ていないが、距離程100〜200m区間では、不明瞭ながら２法をプログラム化した。このプログラムの特色は以下の通つの低速度帯（速度値不明）の存在が認められている。

りである。
①

はぎ取り法は、地盤を凹凸の少ない層構造と仮定し解析速度分布図（地盤の推定速度構造）の表現が可能でしているが、実際、層構造と仮定出来ない地盤については、ある（図−10参照）。また、地盤内部における弾性解析結果の精度が低下している可能性がある。したがって、波の伝播状態を直観的に理解するのに有効な等走時順解析手法を用いて、はぎ取り法の解析結果がどの程度観面図や弾性波の経路を描いた波線追跡図も出力可能測走時を説明できているか、シミュレーションを行った。

である。

順解析では、はぎ取り法解析結果をセルに分割し、モデ入力データ（観測走時データおよび初期速度構造モル化し、LTI法にて計算を行った。モデル化に当たっては、デルデータ）の設定のあとは、自動的に反復計算をセルに一定の速度値を与えた。はぎ取り法解析結果では、繰り返し、最適な速度構造を求めることができる。

なめらかな層境界が示されているが、今回は長方形セル本プログラムを屈折法弾性波探査解析に活用する事によ（セルサイズ：5×2m）にて、近似的に層境界を設定した。

り、これまで費していた解析作業時間を大幅に短縮できる

また、距離程100〜200m区間での低速度帯は速度値不明なとともに技理学的に客観的な速度構造評価が可能となる。

ことから、設定しなかった。順解析結果として、波線追跡②

図を図−12に、走時曲線図を図−13に示した。
走時曲線図においては、黒線は観測走時を、赤線が順解析で得られた計算走時を表している。その結果、観測走時と計算走時の差は、平均して 9. 3 msであった。
図−13のように、今回のシミュレーション結果では大きく走時が異なる箇所が認められており、この理由として、

12

こうえいフォーラム第９号／２００１．１①実際の地盤が、はぎ取り法で仮定した凹凸の少ない層構

表−２ 解析データ一覧

造ではなかったこと、②図式解法のため、解析者の主観がより多く介在したことという可能性がある。したがって、トモグラフィ的手法を用いて、より観測データを説明できる速度構造を推定することとした。
においては、白線にてトンネルフォーメーションを図示し（2）トモグラフィ的手法による解析結果
トモグラフィ的解析における解析データの一覧を表−２た。また、走時曲線図においては、黒線が観測走時を、赤線が逆解析後の計算走時を示している。
本解析結果を以下にとりまとめた。に示した。
順解析手法としてLTI法、逆解析手法としてSIRT法を用①

はぎ取り法解析結果に比べ、凹凸に富む速度構造といた。また、RMS残差（残差自乗和の平均とした）の減なっている。これは、はぎ取り法が凹凸の少ない成少がみられなくなった時点での速度構造を最適解とした。

層構造と仮定して解析するのに対し、本解析では層解析結果として、速度分布図、波線追跡図、走時曲線図

構造を仮定せずに、セルに分割して解析しているかをそれぞれ図−14〜16に示した。速度分布図、波線追跡図

らである。

速度(m/s)図−11 はぎ取り法解析結果

図−14 トモグラフィ的手法解析結果
白線：トンネルフォーメーション

図−12 波線追跡図（はぎ取り法）

図−15 波線追跡図（トモグラフィ的手法）
白線：トンネルフォーメーション

図−13 走時曲線図（黒：観測走時、赤：計算走時）

図−16 走時曲線図（黒：観測走時、赤：計算走時）

13トモグラフィ的手法を用いた屈折法弾性波探査解析法の開発②

はぎ取り法解析結果で基盤速度として推定されていより観測データを説明できる速度構造を得ることがる3.4〜3.6km/sの速度境界が、逆解析結果では、距出来た。

離程約150〜500m区間でより標高の低い位置に認め

当解析手法は計算機を用いた解析法であるため、従来のられる。

はぎ取り法と比べ、解析時間やコストを縮減できるものと距離程約50mから200mの区間で、地形に対してお考えられる。また、探査手法（データ取り）も従来と変わおよそ2.0km/sの速度層が凹状の構造を示している。

らないことから、既往データの再解析に適用することがで走時曲線図によると、平均の走時差は2.6msであり、きる。今後はより多くの探査データについて解析を行うこ（1）の順解析結果の9.3msより小さい。とにより、その適用性、ならびに適用限界について検討す③④

図−13と16のそれぞれの走時曲線図を比較すれば分かるとおり、今回の弾性波トモグラフィ技術を応用した屈折る予定である。また、３次元解析への拡張性についても検討する予定である。

法弾性波探査解析結果は、はぎ取り法解析結果より、より観測データを説明できる速度構造であると思われる。ただ謝辞：本研究を進めるにあたり、中央研究所開発研究部高し、本解析では波線の通過領域に対して速度の修正を行う

橋秀博士、今野正雄氏、総合技術センター鈴木弘明氏、仙ため、図−15の波線追跡図を参考にすると、波線分布が密台支店中村一樹氏より多くのご助言をいただきました。こである表層付近は、結果の精度・信頼性が高く、逆に分布こに記して謝意を表します。

が粗である山体内部は精度・信頼性の低い箇所であることに注意する必要がある。当手法は新たな観測データを用いて追加解析を行うことが可能であり、波線分布の粗な領域や施工上問題点となりうる箇所を対象に測定を行うことで、より高精度で信頼性の高い速度構造が推定出来るものと期待される。

参考文献
１）物理探鉱技術協会：土木地質調査における屈折法弾性波探査の調査計画に関する要綱、1979
２）佐々宏一、芦田譲、小澤岳史、山田雅行：効果的なレイトレーシン
グによるサイスミックトモグラフィの精度向上について、物理探査
学会第78回学術講演会論文集、pp.92-93、1988
３）林宏一、斉藤秀樹：高精度屈折法地震探査の開発と適用例、物理探査、Vol.51、No.5、pp.471-492、1998

６．お わ り に
弾性波トモグラフィ技術を屈折法弾性波探査解析に応用し、その解析手法を開発した。その結果を以下に取りまとめた。
①

を改善できる解析手法を開発した。
当手法を用いて、トンネル計画地点における屈折法弾性波探査解析を行い、はぎ取り法の解析結果より、

14

traveltime interpolation，Geophysical prospecting，Vol.41，pp.99-111，
1993５）横田裕、伊藤信一、野田利一、香川敏幸：トモグラフィー的手法を用いた屈折法解析−手解析と比較して−、物理探査学会第91回学術講演会論文集、pp.71-74、1994

弾性波トモグラフィ技術を、屈折法弾性波探査解析に応用することにより、従来のはぎ取り法の問題点②４）Asakawa，E. and Kawanaka，T.：Seismic raytracing using linear６）Qin, F, Y. Luo, K. B. Olsen, W. Cai, G. T. Schuster : Finite difference
solution of the eikonal equation along expanding wavefronts，Geophysics，
Vol.57, pp.478-487, 1992
７）Boschetti, F., Dentith, M., and List, R.:Inversion of seismic refraction datausing Genetic Algorithms, Geophysics，Vol.61，pp.1715-1727, 1996