第41号

国立防災科学技術センター研究報告

1988年3月

550．34＝681．3

震源計算・発震機構解計算プログラムの改良岡田義光＊
国立防災科学技術センター

Improvememt

Hypo㏄11terofDetermimation

RoutimandPmgmms

Mechamism

forCa1cu1ation

By

Yoshimitsu

Okada

M肋〃地∫〃肋α〃召プ〃1）ゴ∫α∫〃1〕κ舳〃ケo〃，ノ；ψ舳Abstmcttion

The

formerroutinewerecheckedandsystem，theAs

Ana1yzingtotheSystem

hypocenter

traveltime

calculation

times

later

than18h00h00m．Improvement
Asto

the

algorithm

usedadjustment

after

ang1e．Using

asolution

withprograms

revisedthetotwoand

theawasprogram，amechanism

that
did

the
not

inputwhich

CPU

time

by

seriousspherical

theandintroduged

data

solution，Thisestimate

fomd

iteration

reducedcoarsenew1ycalcu1a−data

processing

Earthquakes．

thebeginning，the

a1so

the

arriva1time

best

using

determination

inof

that

problems

the

from

them

program，itandcalcu1ation

for

score

low

inadequate

searchgrid

best

Precursors

thesesearch

fine

for

quite

mechanism

to

the

due

hypocenter

installdetermination

was

of

for

to

efficiency
we11forwere

countedprob1em

was

fomd

in

the

out

by

fine

grid

interva1andextent

of

the

from

carriedbeing

routinetheoriginha1f．

wasrevisedofwork

wascoarseable

ambiguity

toofrotation

obtainthetheP−andT■aXeS．

1．はじめに
国立防災科学技術センターの地震データ処理ルーチン作業で従来より用いられてきた震源言十算及び発震機構解計算のプログラムを，新しく導入された地震前兆解析システムに移植するにあたって，そのアルゴリズムを見直し，若干の改良を行った．

‡第2研究部地震活動研究室

153一

国立防災科学技術センター研究報告第41号

1988年3月

本稿では，旧プログラムに内在した問題点及びその改良結果について述べる．

2．震源計算プログラム

従来のノレーチン作業で便用されてきた震源計算プログラムは鵜川ら（1984）によるものであり，その中核部分のフローは図1に示すようになっている．この各部分における処理内容は以下の通りである．

①震源計算の基礎データである各観測点におけるP波，S波の着震時分秒（IH，IM，SS）を，①D・t・tti・gT．b。

②1・itiωg・…

XO．tO

i＝O
i＝i．1i〉10？

y電③Tm舳im砒・1・ti・・T。。ユ1藷，纂④1・・…i・・

Xi，ti⑤C・㎜哩・㏄yes
図1⑥．誰畿鴫。峰⑦M

dj・・t㎜・tX，tFig．1

震源計算プログラムの中核部分の処理フロー
Flowofthehypocenter

mainpartof

thedeterrnination

prog「am・その日の午前零時からの積算秒丁。。。＝IH＊3600＋IM＊60＋SSに変換する．

②3観測点のP時刻と1観測点のS−P時問を基に，二次方程式を解いて第0近似解（震源位置X。，発震時t。）を求める．

③理論走時丁、。1及びその震央距離△，震源の深さZに関する微分値を計算する．理論走時は△の5次式丁。。1＝A。十A1△十………十A。△5によっており，係数A。〜A。は震源の深さ別に震央距離区問毎に表として与えられているので，任意のZについてはこれを狭む数点の深さでのT。。1及び∂T。。1／∂△を求め，多項式によるラグランジュ補間を行う．
④T。。。一丁。。1を入力として，∂T。、1／∂△，∂T。、1／∂Zを係数行列とする重み付き最小二乗法の計算を行い，次の逐次近似解（Xi，ti）に進むための修正量（δX，δt）を求める．この過程で震源の深さが負になった場合は，以下，深さを35，25，15，5kmの順に固定し，収束時点で最も解の標準偏差の小さいものを選択する．

一154一

震源計算・発震機構解計算プログラムの改良一岡田⑤解の収束を判定する．収束条件は，緯度，経度，深さ，発震時の修正量を各々δ1at（deg），δ1ong（deg），δZ（km），批（sec）として，STPF＝110＊1δlat

l＋90＊■δ1ong■

十〇．5＊1δZl＋5＊1δtlの値がO．05より小さくなった場合とする．これは，各量単独なら，各々50m，50m，100m，O．01secの修正量に相当する．解が収束するまで，上記③〜⑤は最大10回まで繰り返される．

⑥O−Cの大きいデータ（Pについては2秒，Sについては4秒以上）を除外して再び③からやり直す．除外するデータのない場合，または除外作業が3回繰返された場合，震源計算を終了し⑦へ進む．

⑦積算秒で得られた発震時を，歴年の時分秒に戻す．

これらの各部分について詳細な再検討を行った結果，以下のような問題点が判明した．
（1）②の第O近似解を得る2次方程式の解法にコーディングミスがあった．
（2）④の逐次近似解を得るための逆行列計算の過程で，対角要素が1O−9を下まわった際になされる列置換にコーディングミスがあった．

（3）深さを固定とする時，最初の35kmについては，逐次近似の回数が最大10回ではなく，深さが負になった時点からの残り回数しか許されない．
（4）⑤の収束判定でδ1at，δ1ongの単位に誤りがあった．

（5）⑥のデータ調整で発震時の積算秒が負になった場合，月の変り目のところでは正しい結語が得られない．

上記（1）〜（5）はごく特別な場合にしか問題にならないバグであり，全体への影響は殆どない種類のものであるが，完壁を期すため修正を行った．

（6）震源計算全体で最も負荷率の高い③の走時計算部分の能率があまりよくない．△の5次式の計算が直接A。＊△＊＊5＋A。＊△＊＊4＋………となっていたのを（A。＊△十A。）＊△十………のような形に改め，ラグランジュ補間の定式化についても割算の回数を大幅に減らす等して簡潔な形に改めた．この結果，1984年の1年分約15000件の震源計算処理に，従来プログラムでは38分44秒のCPU時間を要していたのが27分24秒に短縮された．
（7）解が収束するまでの逐次近似の回数の分布を，1984年分について発震時別に調べた結果を図2の実線で示す．通常は4回前後で収束しているが，18時以降になると，最大の10回まで繰返されて計算が打切られる例が急増する．これは，①で震源計算の入カデータをその日の午前零時からの積算秒としたことに問題がある．⑤の収束判定で，発震時単独ならその修正量はO．01秒以下になる必要がある．しかし1日は86400秒であるため，遅い時刻では，発震時の修正量が32ビット単精度浮動小数点データの精度である10進7桁から外れ，逐次近似が働かなくなってしまう．そこで，時間の積算基準をその日の午前零時ではなく最寄りの時分とするように改良して，同様に発震時別の逐次近似回数分布

一155一

国立防災科学技術センター研究報告ト舳

1988年3月

図2

解収束に要した逐次近似回数の発震時別分布．実線は従来ルーチン，点線は改良ルーチンによる（データは

Fig．2

Hour1yl1・鮒

h

12h

01h

13h

02h

14h

03h

15h

04h

16h

05h

17h

06h

18h

07h

19h

08h

20h

09h

21h

10h

22h

OO

第41号

1984年の1年分）
of

23h

11h

12345678910
Nu1nb6rof

iterotion

1needed2

3

4

toconverge

5678910

distribution

iterationofthenumberneededtoconverge

inlinesareduethe1984．

Solid

formerroutinewhile

dottedare

torevisedroutine．

duethetolines

を求めてみると，図2の点線のようになった．今度は発震時による系統的なエラーがなくなると同時に，以前に存在した，やはり何らかの計算不安定性に起因したと思われる，収束回数が2回のところのピークも消えている．

以上の処置に加え，P波，S波について同じ処理が重複してなされていた部分の整理等，細かい改艮を行った結果，前記1年分の震源計算処理に要するCPU時問は更に18分29秒に短縮された上，震源決定数の増大，深さ固定となった数の減少等，若干の精度向上も得られた．
従来のノレーチンと最終的な改良ルーチンとの処理性能の比較を表1に示す．

3．発震機構解計算プログラム

従来のルーチン作業で使用されてきた発震機構解計算プログラムは，堀内ら（1972）と同様のグリードサーチ方式によるものであり，図3に示すようなフローとなっている．この各部分における処理内容は以下の通りである．

①極性読取数Nが1O個以上ある場合のみを計算の対象とする．

一156一

震源言十算・発震機構解計算プログラムの改良一岡田

表1

震源計算プログラムの従来ルーチンと改良ルーチンとの処理性能比較

Table1Comparisonoftheperformance

theformerand

therevised

Datain1984

CPUtime（usingnot

determinedof

hypocenter

determination

Routine

RevisedFormerbetween

M360）

Routine38m44s

18m29s5262

5223

hypocenter

determined

10275

10314

depth

among

1258

1245

fixed

programroutines．

above①N・mb・・佃k②・脈舳刺d・・1・ti。・③S榊i・dg・舳耐i・・Ng．i・
1−Ng・id④舳い舳i州・。。…〕⑤S・。爬舳山ti。・⑥S・1・・t5t帥岬1・・1・ti。・・

図3

so11−5⑦Fi…dj・・tm・・川・・1・ti。・・

発震機構解計算プログラムの中核部分の処理
フロー（従来ルーチン）

Fig．3

F1owofca1cu1ation

themainpartprogram（theofformet

themechanismroutine）．

②震源から各観測点への射出角を，震源の深さと震央距離別に与えられた表からの1次内挿により求める．

③半球面上に格子点網を設定する．極性読取数Nが45未満なら15。メッシュの85格子点，45以上なら10。メッシュの145格子点をとる．

④上で設定された各格子点にB軸を置き，そこを中心にP・丁軸を回転させた時の極性一致数の変化を見る．回転の刻みはデータ数が45未満なら10。，45以上なら5。とし，極性一致数の列を作る．

⑤上記の列のうち，最初に現れる極性一致数最大の並びの中問点にあたる回転角に相当する解をその格子点での最良解とみなし，そこでのスコアSC二（極性一致数）十Σ1iri／Nを求める．ここで，1iは極性一致なら十1，不一致なら一1をとり，riはradiation

一157一

pattemに国立防災科学技術センター研究報告第41号

1988年3月

よる係数である．

⑥全格子点について上記の④⑤を行った結果を基に，次の様に5つの仮解を選択する．第1の仮解は上記スコアの最大値を与える格子点，第2の仮解は極からサーチしていって最初に見つけた最大極性一致数の点，第3の仮解は赤道からサーチしていって最初に見つけた最大極性一致数の点，第4，第5の仮解は再び極及び赤道から同様にサーチを行う．

⑦上で得られた各仮解のまわりで解の微調整を行う．データ数が45未満の場合，まず仮解の周囲に7．5。離れた6つの副格子点を設け，そこをB軸として5。刻みの回転を行わせ，⑤と同じくスコアを算出する．仮解よりも大きいスコアを与える副格子点があった場合，仮解をそこに移す．次に，その周囲に3．75。離れた6つの副格子点を設け，2．5。刻み回転による同様の微調整作業を行って最終解を確定する．データ数が45以上の場合も全く同様に〔r離れた副格子点，5。刻み回転〕による第1次微調整，〔2．5。離れた副格子

18げ

Focoユdepth

170o

樽薫葦雫；薫1111葦11ガ

160
150

140

＿＿

一一、一＿200130．

2o、・仁120o3｝110
0否

Φ100

わ．

、一一＿100．一一一＿＿80．

−o
← 90．

5．一．

32．01

80
70．．．…1・996030
莞

申≡堰

5040

0102030

40 50 60 70 80 90 10◎ 110 120 130 140 150km
Epicentro工図4Fig．46050
4035distonce

震源の深さをパラメータとした，震央距離と射出角との関係．実線は従来ルーチンのもの，点線は改良ルーチンで付け加えられた部分
Re1ation

betweenSolid1ines

revisedaretake・off

definedinangletheandformer

epicentra1distanceroutineroutine．

一158一

and

atdifferent

dotted1inesarefoca1depths．

addedinthe震源計算・発震機構解計算プログラムの改良

岡田

点，2．5。刻み回転〕による第2次微調整を実施する．
これらの各部分について詳細な再検討を行った結果，以下のような問題点が判明した．
（1）①で極性読取数は10個以上に限っているが，もっとデータが少なくても，観測点配置が良好なら発震機構解をかなりよく推定できる場合がある．

（2）②で射出角を計算する基として用いられていた表は，図4の実線で示す如く，深さについては40kmまで5km刻み，それ以深は50，60，80，100，200，400，600kmの各深度で与えられている．ここから射出角を内挿する場合，震源の深さが30〜35kmで震央距離が約70km以上の場合，及び震源の深さが5km未満で震央距離も数km以内の場合には，いずれも内挿値が不安定となる．

（3）③で作る半球面グリッドのうち，赤道におかれたものの半数（データ数45未満のときは12個，45以上のときは16個）は，発震機構解の問題の球対称性により，余分なものである．

このため，計算の重複がみられると同時に，⑥で5つの解を選択する際に，第3，第5の仮解が全く同一のものになる場合が多い．

（4）⑤で計算されるスコアは，必ずしもその格子点にB軸を置いた際の最良のものではない．

このため，⑥で選択を行う第1の仮解は，各格子点の代表者でないもの同士を比較していることになる．

（5）③④で粗い格子点間隔と回転角刻みでサーチしたのち，最後に⑦で解の微調整を行っているが，このような方式で最艮解にたどりつける保証はない．1984年1月分のデータで極性読取数15以上の27件について，格子点問隔を2。，5。，10。，15。の4通り，回転角刻みを0，2。，0．5。，1。，2。，5。，10。の6通りにとって，その各組合せでサーチを行った場合の最小極性不一致を調査し，図5に示した．ここで，従来のルーチン作業で使用されていた組合せは各図における1番右下に相当するが，半数近い例で最良の極性不一致数を見つけられないことがわかる．また，その格子点のまわりの微調整によって，必ず最良の所にめぐりあえるとは限らない．

以上のように，従来のルーチンプログラムについては色々と問題点の多いこと，そして1つの地震について5つの可能な解が与えられるだけでその信頼度の指標がないのは，データベース化する上で具合が悪いこと等から，発震機構解計算についてはプログラムを全面的に改訂することとした．改良ルーチンのフローを図6に示すが，その各部分の処理内容は以下の通りである．

①半球面上に格子点網を設定する．格子点の粗さについては，図5に示した調査結果，及び5件のサンプル実行に要したCPU時問が，5。問隔だとO分54秒，3。問隔だと3分11秒，2。問隔だと7分10秒になったことを勘案して，3。問隔を採用し，赤道上での対称点60個を除いた1801個よりなる格子点網を設定する．

②極性読取数Nが6個以上ある場合を計算の対象とする．

一159一

国立防災科学技術センター研究報告

oooooooooooooooooooooooo第41号

1988年3月

ooooo6守oo守
oooooo守oo守○うoooo060つooNNN000000

い㎝NNN㎝㎝㎝NNαNnNNN〔NNα㎝αNNooooooooooooooooooooo◎oooo◎ooo㎝㎝N㎝（NNoooooo
ooooooooooooooooooooooooooooooOOOOOOooo◎ooooooooOOOOOOooooooooooooooooooooo〇一一〇〇〇〇〇〇〇〇〇〇〇〇

oooo〇一ooooo−

oooooooooooo
oooooo〇〇〇〇〇口ooooooooooooαい㎝N000∋㎝㎝㎝㎝㎝㎝αNN㎝N㎝副N㎝nα

OOOOOOoooooooooooooooooooooo〇一ooo〇一一〇〇〇〇〇一〇〇〇〇〇〇〇〇○ooo〇〇〇〇〇〇

oooooo〇〇〇〇〇〇o◎ooooOOOOOOoooooooooooo
図5Fig．5

ooooooOOOOOOoooooooooooo

寸守守守守守〔αNN㎝㎝

守守守守守守α㎝ぺNn㎝

αα○，㎝N㎝

N㎝NN㎝㎝o〜㎝elo］㎝N

ぺα㎝N靹㎝

格子点間隔および回転角刻みを変えた場合の，最小極性不一致数の変化（1984年1月の例）．N：極性読取数，M：マグニチュード，H：震源の深さ
Changesphericalofthegridminimum

mmber

interva1andof

rotationpolarity−inconsistentmitinreporteddata，M：magnitude，H：focaldepth．

一160一

dataasJan．，1984−N：mmber

aof

fmctionofpolarity

震源計算・発震機構解計算プログラムの改良一岡田①Sph・・i・・lg・i・f・m・tb・Ng，id②N・mb…h・・k③T・k・一。ff㎝gl…1・ω・㈹・
1−Ng・id④p一，・一・m・di・t・ib・ti。・⑤A・…g・P一，T一・・・…dth・ir・・liαbility⑥R・・k・沽丁一㈱・
1−Ng・id

図6⑦S・。…α1・・1・li・・f。・・一眺

Fig．6⑧S・1…1・・b・・t・山tb・

発震機構解言十算プログラムの中核部分の処理フロー（改良ルーチン）
F1owofnism

ca1culationthemainrevisedroutine）．

ofmecha−programpart（the③震源から各観測点への射出角を求める．内挿に用いる表は，図4の点線で示すように，深さ31．99km及び32．01kmに対応する射出角データを加えると同時に，震源が深い場合のデータも補充した．また，震源の深さ及び震央距離がともに5km末満の場合には，表からの内挿ではなく，震源と観測点の位置から直接に射出角を言十算する方式に改めた．

④各格子点にP軸または丁軸を置きO．2□刻みの回転を行った際に得られる最良の極性不一致数を求める．

⑤P軸，丁軸の各々について，最小の極性不一致数を与える格子点の分布の平均方向及びその標準偏差σ二〉（Σαi2）／（N−1）

を求める．ここで，αiは，平均方向と該当格子点方向とのなす空問角である．

⑥上記のσが5。，10。，20□未満の場合，そのP軸または丁軸方向の信頼性をランクA，B，Cと定義し，20。以上の場合はDと格付けする．なお，σ。，σ丁が共に40。を超える場合は，解不能として発震機構解言十算を打切る．

⑦各格子点にB軸を置きO．2。刻みの回転を行った際に得られる最良のスコアを計算する．

⑧全格子点のうちで最大のスコアを与える格子点の方向をB軸とし，該当する回転角に対応したP軸・丁軸を有する解を最良の発震機構解として採用し，上記⑤で得られたσ・，所を解の信頼度として付与する．

一161一

国立防災科学技術センター研究報告第41号

1988年3月

9．おわりに改良された震源計算プログラムについては，現方式でもっとも能率の良いアルゴリズムのものが出来上がった．今後は，観測点補正値の導入，3次元速度構造の配慮といった質的向上が望まれるところである．一方，発震機構解計算プログラムについては，信頼性の高いアルゴリズムに基づいた改艮ルーチンが得られたものの，計算時問が非常に多くかかる難点がある．極性読取データの位置がアナログ的に分布する以上，グリッドサーチの格子点間隔はいくら密にしてもきりがない．今後は，青木（1986）の試みのように，何か違ったアプローチから，よりスマートに発震機構解を求めるためのアルゴリズムを考案する必要があろう．

本改良作業にあたっては，当センター第2研究部鵜川元雄，石田瑞穂，松村正三の各氏に色々と御教示を戴いた．ここに記して謝意を表する．

参考文献

1）青木治三（1986）：フーリエ解析による地震メカニズム解．地震学会講演予稿集，No1，283．
2）堀内茂木・江村欣也・平沢朋郎（1972）：気象庁観測網を用いる場合のP波による発震機構決定の精度．地
震，25，92−104．

3）鵜川元雄・石田瑞穂・松村正三・笠原敬司（1984）：関東・東海地域地殼活動観測網による震源決定方法について．国立防災科学技術センター研究速報，53，88pp．

（1987年11月13日

一162一

原稿受理）