火山 第 48 巻 (2003)
第 5 号 407ῌ413 頁

論 説

火山地域における震源計算についての提案
西潔῍῎石 原 和 弘῍῍ῐ2003 年 2 月 3 日受付ῌ 2003 年 8 月 14 日受理ῑA Proposal of Hypocenter Calculation for Volcanic Regions
Kiyoshi N>H=>῍ and Kazuhiro IH=>=6G6῍῍A new scheme for hypocenter calculation is proposed for volcanic regions. The scheme ﬁnds the global
minimum of the travel time residual, and the resulting solution is more stable than that from the conventional
Geiger’s method, especially when the number of observations is small or the station distribution is unsuited. In theﬁrst step (Step 1) of this scheme, the target area is parameterized by the node. The node distance depends on theheterogeneity of the velocity structure; usually about 0.3ῌ1 km is su$cient. Travel times between nodes andstations are stored. In the ﬁrst half of the second step (Step 2A), a node is sought that minimizes the sum ofsquares of the travel time residuals for the event. By changing the depth of this node, a set of initial hypocentercandidates is generated. In the following step (Step 2B), with these candidates for the initial hypocenter, precise
locations of the event are obtained by nonlinear calculations using the simplex method. The hypocenter withminimum travel time residual is then selected as the most probable hypocenter. For the travel time calculation, the3-D robust seismic ray tracer known as Fermat (Nishi, 2001) is used. If the velocity structure and the locationsof the observation stations do not change, the results of Step 1 are e#ective for all events. Consequently, repetitionsof only Step 2A and 2B are su$cient for every event. Successful outcomes of calculation using synthetic and actual
data are obtained with practicable CPU times.
Key words: nonlinear hypocenter calculation, initial hypocenter, simplex method

1. は じ め にいては地形上の制約からこれらの装置をもってしても震源分布は応力の集中場所を反映しているのでῌ 火山

デ῏タの取得及び送信が不可能なことがありῌ 地震計を地域の場合はマグマの滞留場所や移動経路を含む火山体

震源決定上望ましい場所ῌ すなわち震源を取り囲む場所の構造に関する情報を含んでいる῍ またῌ 震源分布の時に配置することが困難な場合がある῍ またῌ 十分な数の間的変化はマグマ活動の時間的変化を示しており火山活観測点を設置できない場合もある῍ さらにῌ 火山性地震動の推移を見極めるうえで重要である῍

は構造性の地震と比較すると震源が浅くῌ 規模も小さい

近年ῌ 高度に発達した半導体素子と通信網および GPSものが多い῍ したがってῌ 観測点を多く配置しても震源システムを利用することにより火山地域における地震波

計算に使えるデ῏タ数はすくなくなる῍ このように拘束形取得能力は質 ῎ 量共に飛躍的に向上した ῐ例えばῌ

条件が弱くなるとῌ 線形化により行列演算を伴う Geiger

Iguchi, 1991; 森田῎浜口ῌ 1996ῑ῍ しかしῌ 火山地域にお法による震源計算では解が不安定となり発散や振動が起きる ῐ例えばῌ Sambridge and Kennett, 1986ῑ῍

῍ ῒ606ῌ8223 京都市左京区田中東樋ノ口町 24ῌ2
24ῌ2 Tanaka Higashi Hinokuchi-cho, Sakyo-ku,
Kyoto 606ῌ8223, Japan.
῍῍ ῒ891ῌ1419 鹿児島県鹿児島郡桜島町横山
京都大学防災研究所附属火山活動研究センタ῏
Sakurajima Volcano Research Center, Disaster
Prevention Research Institute, Kyoto University,
Sakurajima, Kagoshima 891ῌ1419, Japan.
Corresponding author: Kiyoshi Nishi
e-mail: nah02621@nifty.com

行列演算を伴わない非線形震源計算法としてシンプ
レックス法がよく用いられるが ῐ例えばῌ Prugger andGendzwill, 1988; 西ῌ 2001ῑῌ 計算には初期震源を必要とする῍ 震源計算の対象領域内に local minimum が存在する場合ῌ 初期震源が適切でないと解が local minimum に収束してしまう危険性がある ῐ例えばῌ Prugger andGendzwill, 1988ῑ῍ この点を解決するためῌ 本稿では複数の最適初期震源候補を用いてシンプレックス法による非線形震源計算を行いῌ その結果の中から最適震源を選び

西

408

潔῎石原和弘出して最終結果とする計算スキ῏ムを提案する῍
震源と観測点の走時計算には Nishi (2001) によるῌ 火

頂点を含めた 2ῐ4 個程度で良い῍ Node の間隔としてはῌ 0.3ῐ1 km 程度である῍

山地域の 3 次元速度構造に対応した波線追跡法を用い

2) 速度構造をブロックの頂点 ῑ格子点ῒ に付与したる῍ この波線追跡法は Fermat と命名されているのでῌ

値あるいは V῕f (x, y, z) で与える῍ ここで V は座標 (x,

本震源計算スキ῏ムを HypoFermat 法と称する῍y, z) における地震波速度である῍
3)観測点と node 間のῌ P 波と S 波の走時と波線経

2. HypoFermat 法による震源決定路をῌ 地震波線追跡法 Fermat (Nishi, 2001) の前半の手

HypoFermat 法は Step 1 と Step 2 で構成される῍ Step法を用いて求める῍ この方法ではまず全ブロックについ1 は震源計算領域内に node を設定しῌ 各 node と観測点てῌ 同一ブロック上にある全ての node 間で走時を計算

間の計算走時ῌ 波線経路を求め格納する῍ この node はする῍ 次にそれらの足し合わせとして観測点と任意の初期震源候補としてῌ また node の点列は波線経路の近node 間の走時を表現しῌ 複数ある経路の中から走時が似としてῌ 用いられる῍ Step 2 はῌ Step 1 の結果と event

最小となるものを初期波線経路として採用する῍

毎の初動デ῏タから初期震源を求める Step 2A とῌ 初期

こうして得られた走時及び波線は保存されῌ 走時は震源を用いて非線形震源計算を行い真の震源を見つける

Step 2A において初期震源の探索にῌ また波線経路は

Step 2B からなる῍

Step 2B において node 以外の任意の点と観測点間の走速度構造と観測点座標が変らない限り Step 1 の結果時を求める際にῌ それぞれ用いられる῍ 波線経路は nodeは全ての event に共通であるからῌ Step 1 は最初に一回の点列で近似されるのでῌ 波線経路の保存はῌ 波線経路だけ実行すればよい῍ 処理の流れを Fig. 1 に示した῍ 各を表す node の接続の順番を保存しておくことで達成さStep の具体的な内容を以下に示す῍れる῍

2ῌ1 領域のパラメῌタ化と走時計算 (Step 1)

速度構造が一次元成層構造で与えられる場合には効率

対象領域をブロックに分割しῌ 稜線上に node をを考慮して波線追跡法 Fermat の代わりに通常のスネル

設定する῍ 最適の node 間隔は速度構造の不均質の程度(Snell) の法則を用いて走時を求める῍ これは Step 2B にによるがῌ 通常はブロックサイズ 1 km 程度ῌ node 数はおける震源計算においても同様である῍

1)

2ῌ2 初期震源候補の自動探索 (Step 2A)
obs1) 観測点 i の P 波初動時刻 T obs
i はῌ 観測走時を t i ῌ

発震時を T とすると
obs
T obs
i ῕t i ΐT(1)である῍ ここでῌ i῕1ῐnp, np: P 波初動デ῏タの個数῍
発震時 T に代わる reference time T obs
0 を P 波初動時刻の平均値で定義すると (Prugger and Gendzwill, 1988),(2) 式が得られる῍
np

np

iῌ1

iῌ1

obs
obs
T obs
0 ῕( ῌ T i )/np῕( ῌ t i )/npΐT

(2)

obs
次に T obs
0 に準拠した観測走時 t i を (3) 式で定義する῍
obs
obs
t obs
i ῕T i ῔T 0

(3)(1), (2), (3) より
np

np

iῌ1

iῌ1

obs
obs
obs
obs
t obs
i ῕T i ῔( ῌ T i )/np῕t i ῔( ῌ t i )/np

(4)

を得る῍
Reference time を用いることにより以後の計算において未知数は (x, y, z) の 3 個である῍ 発震時 T は震源が確定した後ῌ 次式 (5) によって求めることができる῍npcalT῕( ῌ (T obs
i ῔t i ))/np

(5)

iῌ1

ここでῌ t cal
i は震源と観測点 i の間の計算走時である῍
Fig. 1. General ﬂow of calculation scheme.

2)

S 波時刻 S obs
i の読みとり値がある場合は SῌP 時間

火山地域における震源計算についての提案
SP obsi を (6) 式で求める

409

Prugger and Gendzwill, 1988  震源 (x, y, z) を探索するobs
obs
SP obsi S i T i (i1ns)(6)

ここで ns は S 波時刻の観測デタ数である

3)

震源探索に必要となる任意の点 P から観測点までの走時は次のようにして求める 先ず点 P に直近のStep 1 により node j から観測点 i までの計算走時

node を求める この node から観測点に至る最短経路は

t cal
ij が得られる (4) 式は観測値 (obs) を計算値 (cal) とStep 1 で得られ保存されている この経路に P の座標を3)置き換えても成立するので node j と観測点 i 間につい加え 点 P から観測点に至る経路 node 座標の列 を得

calて reference time に準拠した計算走時 t cal
ij を t ij を用いる この経路は node に拘束された点列である この点て表すことができる すなわち (4) 式より次式 (7) 式

列にもシンプレックス法を適用し (Nishi, 2001) 走時がが得られる

最小となる経路と走時を求める
np

cal
cal
t cal
ij t ij ( ῌ t ij )/np

(7)

iῌ1

4)

Step 1 の結果を用い node j と観測点 i 間について 2) に対応する SῌP 時間 SP cal
ij を求める5)4) 震源が得られた後 発震時 T を (5) 式を用いて求める
5)Step 2A から渡された 深さの異なる数個程度の初期震源について 1)4) を行い 最小の走時残差 2 乗

node j について reference time に準拠した P 波走和を与える震源を最終結果とする

時残差の 2 乗と SῌP 時間残差の 2 乗和 Ej を求める
np

ns

iῌ1

iῌ1

cal 2
obs
cal 2
Ej ῌ (t obs
i t ij ) k ῌ (sp i sp ij )(8)

3. HypoFermat 法による計算例

前節で述べた計算スキムを Fortran (F90) を用いて

ここで k は S 波の読みとり精度に対する重みで通常 1.0

コド化し 走時残差に Local minimum がある場合の震以下の値を与える源決定 均質 2 層構造 synthetic デタ  観測点が少な

6)

Ej を最小値にする node j を求めその node の座標

を初期震源とする
震源が node と node の中間に近い場合には 6) で求い場合の震源決定 均質 2 層構造 実デタ  人工地震の震源決定 3 次元速度構造 実デタ に適用し性能を確認した

まる node は震源の深さ方向には最寄の node とはなら

3ῌ1 走時残差に local minimum がある例

ない場合がある これは 観測点間の標高の差が観測点

走時残差に local minimum がある場合の震源計算例の水平方向の広がりに比べて一般に小さいため 深さ方を 震源計算結果と走時残差の root mean square (rms)

向の分解が悪いことによる 走時残差のコンタは深さのコンタとともに Fig. 2 に示した

方向に伸びた形となり 間隔の粗い node 分布では震源

速度構造は標高 0.0 km に不連続面を持つ均質 2 層構に最寄の node の残差が必ずしも最小とはならない場合

造で 第 1 層と 2 層の P 波速度はそれぞれ 2.0 km/s とである このため 6) で得られる node の座標から深さ

4.0 km/s また 2 層とも P 波と S 波の速度比を 1.60 としを変えた複数の初期震源を作成し 次の Step 2B に渡すた 震源位置を (0.5, 0.5, 0.5) に設定し 印  観

Simplex 法による非線形震源計算 (Step 2B)

測点  印 までの走時を求め観測値とした 観測値はStep 2A で得られる初期震源の座標は水平方向に関し

観測点 A が P 波と S 波の到達時刻 他の 3 観測点は P

2ῌ3ては node の座標として得られるので node 間隔以下の分解能はない Step 2B では 初期震源の周りの領域にお波の到達時刻のみとした
真の震源の下方 2.9 km に顕著な local minimum Fig.

いて P 波走時残差の 2 乗と SῌP 時間残差の 2 乗和をさ

2(b)らに減少させて最小とする点 すなわち求める震源を探えば 初期震源の深さが 0.0 km 以深の場合 には local

索する
一般にある関数を最小にするパラメタを探索する方法は最適化 (optimization) 問題として研究され 多種多印 があるため 初期震源が適切でない場合 例

minimum (0.66, 0.39, 2.40) に収束してしまう 一

方 最適な初期震源が与えられる HypoFermat 法では
浅い方の真の震源に収束させることができる

様の方法が考案されている この中で シンプレックス

3ῌ2

法 (Nelder and Mead, 1965; Press et al., 1992) は微分演算

Merapi 火山はインドネシアの火山の中でも活動度のを要せず安定した解が得られるのでよく用いられる 震

高い活火山で VSI (Volcanological Survey of Indonesia)

源決定の概要を次に示す

によって観測と研究が実施されている 1991 年におけるMerapi 火山における例

1) 初期震源の近傍 Step 1 におけるブロックサイズ

VSI の地震観測で明瞭に P 波初動が読み取れた観測点は程度 に初期震源を重心とする任意の 4 点を設定する

Fig. 3 に示すように 4 点であった 山頂火口に最も近い

2) この 4 点にシンプレックス法を適用し 例えばPUS 点が 3 成分地震計で 他の点は上下動地震計が設置

410

西

潔 石原和弘Fig. 3. Locations of stations (solid dots) at Merapi
Volcano in 1991. Triangle indicates the summitof Merapi volcano. Dotted line between P1 andP2 indicates the location of vertical cross section
cutting the determined epicenters (see caption in
Fig. 4). Contour interval of altitude is 500
meter.

うち 収束する震源と振動する震源をそれぞれとで示し HypoFermat 法による震源をで示した
Fig. 4(c) は Geiger 法による震源と HypoFermat 法による震源を通る垂直断面である 断面図の走時残差コン
Fig. 2. Contours of root mean square (rms) oftravel time residual for synthetic data. Contour
interval is 0.025 seconds. Solid square indicatesthe location of the station. Plus mark and blanksquare indicate the location of true hypocenterand the local minimum of rms of travel time
residual respectively.
(a) Horizontal cross section at Z
0.5 km.
(b) Vertical cross section at Y
0.5 km in
Fig. 2(a).

タからは 観測点が少ないために拘束条件が弱く 残差の小さい領域が上下方向に広い範囲に分布していることがわかる このため Geiger 法では震源位置を正しく求められなかった また 正しい震源に近い初期値を与えた場合でも 解は振動し正しい震源は求められなかった この場合 damping を加えることによって振動は抑えられるが damping の値に関しては震源要素の resolution を低下させないように注意深く選ぶ必要があるとされている 例えば Sambridge and Kennett, 1986
されていた したがって 震源決定に用いることができる読み取り値の数は P 波が 4,

S 波が 1 の合計 5 であった 今回は 1991 年 1 月 18 日に発生した火山性地震の着一方 HypoFermat 法では 与えられた計算スキムだけで 走時残差が最小となる位置を正しく求めていることが分る Fig. 4(b), (c) の印

震時デタに Geiger 法と HypoFermat 法を適用し 結

3ῌ3 霧島火山における人工地震の震源決定例

果を比較した3 次元速度構造における計算例を 1994 年に実施され用いた速度構造は標高 0.0 km に不連続面をもつ均質た霧島火山における構造探査 鍵山 他 1995 により

2 層 構 造 で 第 1 層 と 2 層 の P 波 速 度 を そ れ ぞ れ 2.7

得られたデタを用いて示す 震源は人工地震であるか

km/s, 3.5 km/s とし また 2 層とも P 波と S 波の速度ら 計算の結果を震源パラメタの実測値と比較するこ比を 1.86 とした 震源計算結果を走時残差の rms のコとができる 今回の計算には Fig. 5 に示した観測点 発

ンタとともに Fig. 4 に示した Geiger 法による震源の破点のみを用いた 3 次元速度構造は西 (1997) による P

火山地域における震源計算についての提案

411

Fig. 4. Contours of rms of travel time residual for volcanic earthquake at Merapi Volcano. Contour interval is0.05 seconds.
(a) Horizontal cross section at Z῔ῒ1.6 km. Open square indicates the location of epicenter calculated byGeiger’s method with incorrect initial hypocenter. (b) Horizontal cross section at altitude 0.0 km. Cross markindicates the location of epicenter calculated by the HypoFermat method (this study). (c) Vertical crosssection at the dotted line shown in Fig. 3, cutting the epicenters determined by the HypoFermat and Geiger’s
methods. Cross mark indicates the location of epicenter calculated by the HypoFermat method. Open squareindicates the location of epicenter calculated by Geiger’s method with incorrect initial hypocenter. Solidsquares indicate the oscillatory locations of hypocenter calculated by Geiger’s method with correct initial
hypocenter.

波速度構造を用いた῍

走時残差の rms の分布を Fig. 6 に示した῍ Fig. 6(a) はHypoFermat 法 と Geiger 法 の 結 果 を Table 1 に 示 し

rms の深さ 0 km における水平分布である῍ 白丸とΐ印

た῍ Geiger 法としては Thurber (1983) の SIMUL3 に内は Geiger 法で正解が得られた場合と得られなかった場

蔵されている震源計算部分を用いた῍ これは初期値に対合の初期値の位置をそれぞれ示す῍ 発震時の初期値は発する補正値を特異値分解 (SVD) によって求めるもので

破時刻 ῐ鍵山῎他ῌ 1995ῑ とした῍ 初期震源の位置が真成分ごとの固有値が設定値以下の場合はその成分の補正の震源から離れると正解が得られにくい傾向が認められを行わないことにより解を安定化している῍ Table 1 にる῍ これはῌ 深さ方向成分の固有値が許容値未満となる示されるようにῌ Geiger 法では初期震源の位置が適切でことが多くῌ 初期値からの補正が行われないためであないと正解が得られない῍る῍ Fig. 6(b), (c) に示した垂直断面における走時残差のこれに対してῌ HypoFermat 法では初期震源を人の手

rms 分布はῌ 震源周辺を除いてコンタ῏が深さ方向に開で与えることなしにῌ Geiger 法の正解と同等かそれ以上いていて深さ方向に対しての拘束条件が弱いことを示しの精度で震源が求められていることがわかる῍

ている῍

412

西

潔῎石原和弘3ῌ4 計算時間

実デῐタを用いた例についてῌ Step 1 と Step 2 に要する計算時間を Table 2 に示した῍ Step 1 は対象領域内の全ての node と各観測点間の走時を求めるため計算時間を要するがῌ すでに述べたように一度計算をしておけばよい῍ 計算時間は Table 2 に示されるように実用的な範囲にあるといえる῍ 個῏の event ごとに行う Step 2 の計算時間を Geiger 法と比較してみるとῌ 1 次元の均質構造
でおよそ数十倍以下ῌ 3 次元構造でおよそ数百倍以下である῍ 数百倍となるのは Geiger 法が本来高速である為Fig. 5. Locations of stations (solid squares) andcontrolled event (open star) at Kirishima
Volcano in 1994. Open squares A and B
indicate the initial hypocenters for Geiger’s
method. The results of calculations for A and B
are shown in Table 1. Rectangle shows thetarget area for hypocenter calculation.

でありῌ 実時間としては許容の範囲といえる῍
計算機の処理速度は加速度的に向上している῍ 計算時
間の多少はさほど問題ではなくῌ 計算精度が良いことが重要である῍Table 1. Comparative test of Geiger’s method and the HypoFermat method by the data of controlled event at
Kirishima Volcano.

Fig. 6. Contours of rms of travel time residual for the controlled event at Kirishima Volcano in 1994. Contour
interval is 0.05 seconds.
(a) Cross section at Zῒ0.0 km. Solid squares show the locations of the stations. White small circles and crossmarks indicate the locations of initial hypocenters for Geiger’s method. From the initial hypocenter indicatedby white small circle, reasonable result is obtained whereas from the initial hypocenter indicated by cross mark,reasonable solution is not obtained. Dotted lines indicate the location of the vertical cross section shown in
Fig. 5 (b) and (c). (b) Cross section at Xῒῑ0.5 km in Fig. 6(a). (c) Cross section at Yῒ4.5 km in Fig. 6(a).

火山地域における震源計算についての提案
Table 2.

Examples of CPU time for Step 1 and Step 2.

4. ま と め

火山地域における震源計算に適した方法としてῌ 複数の最適初期震源候補を自動探索した後ῌ シンプレックス法による非線形震源計算を行いῌ 最適震源を求める計算
スキ῏ムῌ HypoFermat 法を提案した῍ 走時は波線追跡法に Fermat を用いて求めるのでῌ 高度に不均質な速度構造に対応できる῍ また非線形の震源計算のためῌ マト
リックス演算に伴う不安定性を回避した結果が得られる῍
HypoFermat 法の有効性を確認するためῌ 1 次元速度構造における Synthetic デ῏タと実デ῏タῌ また 3 次元
速度構造における実デ῏タを用いて Geiger 法との比較を行った῍ その結果ῌ HypoFermat 法の解には初期値依存性がないことῌ さらに観測点数が少なく拘束条件が弱い場合にも安定した解が得られることが示された῍
謝

413

辞

メラピ火山のデ῏タ使用に関して Ratdomopurbo 博士(VSI) にお礼申し上げます῍ 二人の匿名査読者と編集担

当の東北大学後藤章夫博士からは本稿の改善に役立つ多くの重要なコメントを頂きました῍ 深く感謝します῍
引用文献
Iguchi, M. (1991) Geophysical data collection using aninteractive personal computer system (part 1)ῌExperimental monitoring at Suwanosejima volcanoῌ. Bull.

Volcanol. Soc. Japan, 36, 335῍343.
鍵山῎他 68 名 (1995) 霧島火山群における人工地震探査
ῌ観測及び初動の読みとりῌ῍ 震研彙報ῌ 70, 33῍60.
森田裕一῎浜口博之 (1996) 火山体構造探査のための高精度小型デ῏タロガ῏の開発῍ 火山ῌ 41, 127῍139.
Nelder, J. A. and Mead, R. (1965) A simplex method forfunction minimization. Computer J., 7, 308῍313.
西 潔 (1997) 構造探査デ῏タを用いた霧島火山体浅部の 3D 速度構造῍ 火山ῌ 42, 165῍170.
西 潔 (2001) 火山体浅部の地震波速度構造解明のための手法ῌ火山地域用 3 次元地震波線追跡法および非線形震源決定法の開発ῌ῍ 京都大学防災研究所研究集会
ῒ火山の浅部構造と火山流体ΐ 報告書 (12K῍3), 61῍69.
Nishi, K. (2001) A three dimensional robust seismic ray
tracer for volcanic regions. Earth Planets Space, 53, 101῍
109.
Press, W. H., Teukolsky, S. A., Vetterling, W. T. and Flannery, B. P. (1992) Numerical Recipes in Fortran,
Second Edition. Cambridge University Press, 402 p.
Prugger, A. F. and Gendzwill, D. J. (1988) Microearthquake location: A nonlinear approach that makes use ofa simplex stepping procedure. Bull. Seism. Soc. Am., 78,
799῍815.
Sambridge, M. S. and Kennett, B. L. N. (1986) A novel
method of hypocenter location. Geophys. J. R. Astr. Soc.,
87, 679῍697.
Thurber, C. (1983) Earthquake location and threedimensional crustal structure in the Coyote lake area,
central Calfornia. J. Geophys. Res., 88, 8226῍8236.
ῐ編集担当 後藤章夫ῑ