hR e p o r t ，  E a r L h q  uakεResearιh I n s l i t u t e ，  U n i v e r s i t y  o f   1 ' okyo ，  No ，  7 ，  pp ，  6

*;;~大守:地震研究所技術研究報告， 010 ，

6 ‑ 1 / 1   U ，  2 0 t H  + 1

T e c h n i c a l  R e s e a r

hR e p o r t ，  E a r L h q  uakεResearιh I n s l i t u t e ，  U n i v e r s i t y  o f   1 ' okyo ，  No ，  7 ，  pp ，  6

H ， 2 0 0 1  

を

F

‑書士

^{* *}  

Development  ( ) f   a  Fast Three

dimensionalElectromagnetic Induction  Solvcr  on  Equation Method 

Hiroshi MUNEKANE* and  Hi ^{sashi UT  ADA}

Ahstract 

A f a s t   s o l v e r   f o r   a three

dimensional induction problem i s   developed  based on the i n t e g r a l   equation method.  The algorithm us

di n  the s o l v e r  i s   b l

don t h a t  o f  Avdeev  e l  a l .   ( 1 9 告 7 ) ， w h . i c h   u t i l i z e d   I t e r a t i v e  Di s s i p a t i v e  Method (IDM) t o  ensur ・ econvεrgence o f  the r e s u l t a n t  Neumann s e r I e s .   The algorithm i s   modified t o  s o l v e  r e s u l t a n t  large and den

el i n e a r   hy the KI

3

l o vsub

pace method ，  instead o f  the Jacobi method.  This modification enhances the e f f i c i e n c y  o f  the solver up t o   t w i c e .   Several variants o f  the Krylov subspace m

dsare t e s t e d ，  and i t   i s   found that the BiCGSTAB  2  method i s   most ε 血 c i e n tf o r  the bεnchmark model o f  the COMMEMI project (Zhdanov  e t  a l . ，  1 9 9 7 ) .  

Key words :  3D  i n d u c t i o n  s o l v e r

l n t e g r a l  

は じ め に

1 ¥ 1 1 ' 1 、 T e l l u r l c ) 法を用いて地下の比抵抗構造 を推定する躍には，構造の二次元性を1Eしく評価する必要 がある.構濯の三次元性の強い場所で される M T レ スポンスは， Currcnt channcling ( c 芯，J i r a c c k 1 9 9 0 )な どの三次芯特有の彰響を大きく受ける. このような場合，

無理に一次元，二次元の解析を行うと誤った構造を推定し てしまう可能性がある.従って.TF.確な構造を推定するた めには，三次川解析を行うことが不可欠である. しかしな がら，三次忌構造に対する電磁応答の計算手法は未だに発 躍進上であり，決定版が与在しないのが現状である.

二 次 元 構 造 に 対 す る 電 磁 の 計 算 の ア ル ゴ リ ズ ム は 大 きく分けて?三種類ある. すなわち有限差分法 CMackiee t   aL ，  1 9 9 3 ，  1 9 9 4  ;  Smith ，  1 9 9 6 a ，  b ) ，有限要素法 CMogi ， 1 9 9 6  

;  Zyserman  and  Santos ，  2 0 0 0 ) ， そして積分方程式法 2 0 0 1 午 、 8 月 1 6 円受付， 2 0 0 1 年 1 0 月 1 5 日受理.

*東京大学地震研究所火山噴火予知研究推進センター，

料海半球観測研究センター，

十現出土交通省国土地理院.

和

V o l c a n oR e s e a r c h  C e n t e r  

料

Ocean H

m i s p h e r e R e s e a r c h   C e n t e r ，  E a r l h q u a k c   R e ‑ s e a r c h  I n s t i t u l e ，  U n i v c r s i t y  o f  Tokyo 

十

Nowa t  G

o g r a p h i c a lS u r v e y  I n s t i t u t e ，  M i n i s t r y  o f  Land ，  I n f r a s t r u c t u r e  and T r a n s p o r t  

6 

l t e r a t i v e 刀 i s s i ρ a t i v eMethod ， ρ a c e  

H J 9 1  ;  Xiong and Tripp ，  1 9 9 5 ) である.

それぞれのアルゴリズムには利点，欠点がある.有限差分 法，釘限要素法は解くべき線形方程式の係数行列が非常に 疎であり，共役勾配法 (CG 法)などのクリ口フ部分空間社;

が適用しやすいという利点をもっ. しかしながら空気層が そデリング領域に含まれていることから，係数行列の各要 さのばらつきが芳しい. このような場合句クリロ フ部分空間 j 去をそのまま適川しても収束が思し場合に よっては発散してしまう.精度良く解くには前処理や，解 が Maxwell 方程式を厳需に満たすように，強制限 j にポテ ンシャル場を差し引く，などの数値計算上の工夫が必裂で あ る 官t ぱ ， 1 9 9 4 ) . 一方，積分方程式法は，モデ リング領域が背景の

次元構造と異なる屯気伝導度をもっ 電気伝導度異常領域に限られるために，空気層はモデリン グ領域に含まれない.従って係数行列の各要素の大きさの ばらつきは有限差分法や有限要素法の場合に比べて小さ し本来クリロフ部分空間法に向いている.しかしながら，

係数行列が一般的に密である上，密行列に適切な前処理が

あまりないため，特に電気伝導度異常のコントラストが大

きい場合には適用は難しかった.一方，庭接法によって線

形方程式を解く場合はモデリング領域が大きくなると係数

行列の次元が高くなり，計算時間やメモリーの￨台 i i から計算

が著しく闘難になるため，その j 血用範囲は従来電気伝導度

積分方 f 宇式 j 去をや用いた

当速二次 J 乙電鈴応答計算コー i ごの開発

実詰領域が小さい場合に限られていた

Avdeev e tα . l は I t e r a t i v e

Method  (IDM) 去に基づく新し L 、アルゴリズムを提唱し た. IDMl 去は Fainberg

ndSinger  ( 1 9 8 0 )において開発 された方法であり，電場および電気伝導度の代わりに新 L

L  、変数を導入することによって，得られた積分方程式が逐 次 反 復 法 抗o b i 法)によって必ず収束する形になって L 、

る.従って， IDM 法の係数行列i まクリロフ部分空間法にも 有利であることが期待される.そこで，本研究ではまず Avde

e la l .   ( 1 9 9 7 ) に準拠した三次兄電総応容計算コー ドを作成することにし，その上でさらなる高連化をはかる ために， Avdccv

e ta l .   ( 1 9 9 7 )では Jacobi 法によっ れ て い た 線 形 程 式 忠 ク リ ロ フ 部 分 空 間 J よそ適用して解 いた.その結果， ， J a c o b i 法で反復するよりも計算時間を短 縮することに成功した.

本論文では，まず第二平で IDM 法のアルゴリズムの披 まとめる

次に第三主主で，三次記モデリングコードの テストのための標準モデルとなっている COMMEMIpro 

(Zhdanov  e t  a l . .   1 9 9 7 ) のモデルについて今回開発し たコードを用いて電般応答を計算する.さまざまなクリロ フ部分空間法を適用した場合について，反復回数や計算時 間などを Jacobii まなどの他の解法と比較してクリロブ 分空間法の釘効性を示すとともに，タリロブ部分空間法の

中でもどの解法が適しているかを検証する.

I t erative  Di ssipative Method (IDM) 去のまとめ

まず，点 r での電場 E 議場 H および屯流密度 j が exp ( i ω t )   ( ωは角周波数)の時間依存性を持って L  、ると仮定する.同波数領域においてマクスウェル方程式 は以下のように書き去される.

¥ 1 > く i ω μH(r)

マ

_ε 

マ

XH(r)

i ; E ( r ) + j s r c ( r )

¥ 1   .H(r)

O ど = σ

z ω ε

ただし ε は誘電率ヲ μ は透礎率司 ρ は電荷密度， σ

導度， C は変位電流を含めた一般化された電気伝導度，j s r c   ( 1 ' ) はソース電流である . Q を背景の一次元成屑構造に対 するグリーンテンソル，そして E n ( 1 ' )を，ソース電流が背 景の構造に対して作る電場であるとすると，ソース電流が 任意の電気伝導度構造に対して作る電場 E ( r ) は次のよう な積分方程式で表現される (Weidelt ， 1 9 7 5 )

E ( r )  = E n C r ) 十ミ f A C ど と 。 )QE(r 勺 d r ' ( 6 )   ただし，ふ ( 1 ' ) は背長の構造に対する一般化された電気伝 導度であり，また S A d r ' は電気伝導度異常が存在する領域 全体における体積積分を表す.グリーンテンソル Q の具体 的な計算方法は Appendix に示している.

7 

式 ( 6 ) を逐次反復法 次 j !i似級数 (Neuman

られる

I Q I I δ ど I 1 

で解く場合，対応する が収束する条件は次のよう and 1 弘 ng ， 1 9 9 7 ) .

(r1¥  べ1)

一 一

⁝ YL a  h xU   1レ

だ た

である. A  vdccv  e t  a   l . ( 1 9 9 7)は式 ( 6 ) を次のように変形して匂この収東条件が常に成¥T.するよ うにした.対応する積分点謹式は次のように書き表せる:

χ

X o ( r )

κ① (RX) f こだし，

( 8 )  

ι 一 一 乙

一

R 

X o

κ①

l  ⁿ

I‑

j s = =  ( ど乙 !)E n

! 

lik 

演算子 κ はベクト

レ η

χ

して次のように作用する.

(κ@η)

η 1 ' Q  η d r '   ( 1 3 )  

( 2 )  

どっはむの共投複嘉数を表す.Pankratov  e t   a   1 . ( 1 9 9 5 )は地 表での電総場のエネルギーの流れが常に外向きであること を手 ) 1 J ‑ l J し て ， この積分ノんら

i

況

1 1 叩 R i く 司 1 4 羽

が常に i 満荷たされることを示した.従って式矧を J a c o l 刀法

く擦の巡次近似級数は必ず収束する.

Avd

v e t  a l .   ( J に お い て は 式 ( 8 ) を Jacobi 法を用

、ている置 しかしながら， 反復法には Jacobi 法や Ga  u s s ‑ S e . i d  e 1 法などのいわゆる定常反復法の他にも，共投 勾 配 訟 な ど に 代 表 さ れ る い わ ゆ る ク リ ロ フ 部 分 雫 問 訟があるCI3 a r r e t e t  a l . ，  1 9 9 4 ) . クリロフ吉 ; 1 分宅問法とは.

n 回目の反復の更新ベクトルを線 ! f c ; h ' 程式の初期残差ベク ト } v から構成されるクリロブ部分空間の中から選択すると いう五法であり，クリロフ部分空間の構成の{上方.またど のような原理で更新ベクトルを決定するかなどによって CG  此 BiCG 法 ， BiCGSTAB 法，おiCGSTAB2 法 ， GP‑

BiCG 法などの種類がある(藤野・ 3 ， 長 1 9 9 6 ) .   一般的にク リロフ部分空間 i 去は Jacobl 法なと、の定常解法に比べて高 速で，特に係数行列が疎の場合に有効であると言われてい る. しかしながらその収束条件は必ずしも明らかで、なく，

式 ( 8 ) に適用した場合， Jacobi 法を用いた場合と異なり得 られた級数が必ず収束するという保証は存在しない.その と，式 ( 8 ) の係数行列は密であるため，クリロフ部分空間法 がそれほど有効でない可能性もある.そこで本研究では，

従来の Jacobi 法などの定常反復法といくつかのクリロフ 部分空間法について，ある標準モデルに対して電磁応答の 計算を行い，クリロフ部分空間法を用いた場合，定常反復 法に比べて収束回数や収束時間は短縮されるか，またいく

( 3 )  

( 4 )  

8  宗包治志;ー歌回久円j

つ か の ク リ ロ フ 部 分 空 間 法 の 中 で ど ということを評価した.

も いる

標準モデルとしては CO.MMEMI

Modeling  Methods  f o r   l n d u c t i o n   Problcmsl  e t  a 1 . .   1 9 9 7 ) において提案 さ れ た そ デ ル そ 採 し た . 国 lに J f Jlサニモデルを示す. こ のそテゾレに対して，田期 1 0 秒 お よ び 1 0 0 秒における M T

レ ス ポ ン ス を 計 算 し た . 収 束 基 準 は 次 に 示 す に よ っ た

川 [κd]‑ 1 ) χ χ 。

f

… ¹ 

^{χ ￨} 。

が 1 0

に 述 し た と き に 収 束 し た と み な し ， 反 複 を 中

lO(Ohm‑m) 

40km  1(Ohm‑m)II00(Ohm‑m)  一I‑‑...‑X‑‑

20km  20km 

判HI 

K 

l 

lO(Ohm-J~l)

20km  100(Ohm‑m) 

O.l(Ohm‑m) 

図1. COMMEMI p r o j e c t において提案されたベンチマー クモデル目 c u 真 1 : から見た肉(下)

図.図中比抵抗の単位 Ohm‑m は i 2.  m を表す.

。[

止 し た 気 伝 導 度 異 常 領 域 と し て ， 原 点 を 巾 心 と す る 8 0 kmX80km の長五形内の深さ 10km までの領域を取り，

その領域を 8 0X  8 0  X  5 の直方体グリッドに区切った.その 際 必 要 と さ れ た メ モ リ ー は す べ て 倍 精 度 計 慌 で 1 6 1 l ¥ l ! b y ( βちであった.比較に用いた反 f 立法は，定常反復法か ら ぬ c o b l 法および準 G a u s s ‑ S e i d e l t:去，またク 1 ) 1 コフ部分 空間法からはBl CGS 、 I ' AB 法 ， BiCGSTA 別法， GP‑BiCG  法である.各解法の詳細については藤野・張(1 9 9 6 ) に詳

し い . 代 表 的 な ク リ ロ ブ 部 分 空 間 法 で あ る CG 法 お よ び BiCG 法については司 CG 法 は 係 数 行 列 が エ ル ミ ー ト で な いと適}はできないのに対し，本研究の係数行列はそうでは ないこと司また BiCG 法は計算のため

る 必 裂 が あ り 繁 雑 と な る こ と や . 収 束 の 仕 が 不 規 制であること ( B a r r e t e t  a ， . l などの理由で採出しな かった.準 G a u s s ‑ S e i d e l i 去とは.本来の Ga  u s s ‑ S e i d e  1  では式 ( 8 ) において i t e r a t i o n の結果求められた χ の僚を巡 次更新するのに対し，本研究で問先したコードは，式 ( 8 ) の 計 算 を 各 層 毎 i こ FFT を用いた Convolutionf i l t e r により

る 寸 前 v e t   a   ， . l J 方式を採用してい るため， 1 直を逐次更新するのではなく属単位でまとめて更

るという方法である.

い ( [ 2 に，周期 1 0 0 秒の計算の i

l 日

反復回

k の!苅が x 向に、 l 正面波が人射した場合の反復 y  h i r l J に人射した場合の反窪田数を示す.

された. BiCGSTAB ，  BiCGST  AB2 ，  GP‑BiCG の い ず れ の ク リ 口 ブ 部 分 空 間

も ， Jacobi  G a u s s ‑ S e i d e l  

て 半 分 程 虚 の 反 復 回 数 で 請 む こ と が 分 か っ た . Table 1  に，反榎討会算に必要な時間を竺関試行を行って測定した結 果を示す . χ ，y方向の入射波に対して要した時間の手rJを示 している f 吏吊した CPU は 21164A( 6 0 0 で ある.これを見ると，所要時間の i 主 i でも， ク 1 ) ロフ 間 j 去を採JIjすることによりお c o b i 法 や 準 G a u s s ‑ S e i d e l

とされる以復計算の時間の半分以下に押えられる ことが分かる.クリロブ部分空間 j 去の中で i t ， BiCGST  AB  2 お よ び GP‑BiCG 1 . 去 が BiCGSTAB; 去 に 比 べ て 反 復 I I I J

数 は 一 割 ほ ど 少 な く な っ た . BiCGST  AB2  i 去 お よ び GP‑

BiCG 法は必要となる反復回数は全く同じであったが，必 要 と さ れ る 計 算 時 間 は BiCGSTAB2 法 の 方 が わ ず か な が ら少なく済んだ. これは偶数岡目の i t e r a t i o n における演 算 が BiCGSTAB2 法のん命が少なくて済むからである(藤 野・張， 1 9 9 6 )   ちなみに今回は Jacobi 法の } f が準 Gauss S e i d e l 法よりも収束が速かった.これは，準 G a u s s ‑ S e i d e l 法では得られた級数の収束性が成立しなくなる影響である

と考えられる.

次に，本計算で得られた結果をこのモデル計算で基準と

される Wannamaker ( 1 9 9   I)の結果と比較する.図 3 に ，

積分方程式 j 去を汁 i いた高速で

次元電総応答 J I 算コードの J I 市 立 ヨ

x ‑ p o l   陶1 。

0 . 0 1  

0 . 0 0 1

^{、、¥ 、 、 、 ¥ 、 、 、、 、、 呂 iCGSTAB}

包0001 。

1e

05 ~ ^{GP 君 iCG}

F  ^{、 、 ， も}

1e

08 1e

G 窃

20  40  60  80  100  120  i t e r a t i o n  number 

y

p o l

0 . 1   J a c o b i  

0 . 0 1   Gauss

S e i d e l

0 . 0 0 1   ^自 iCGSTA 日

0 . 0 0 0 1   1e 心 5 1e‑06  1e 心 8 1e 舟窃

20  40  告 。 8 合 100  120  i t e r a t i o n  number 

図 2 . 周期 1 0 0 秒において各反復法において必安な以復回数. 横 i~111 に反復回数， 縦 軸l こ を 示 す ー 上の i 苅が x 方向に平而波が 入射した場合， ドの凶に y 方 I ( ' Ji こ平面波が入射した場介会示す. BiCGSTAB ，  I 3 i CGSTAD2 ，  GN3iCG などのケリロフ部分空間法 長用いた場合

l a

b i ; 去や子供 G a u s s ‑ S

i d e l 法などの定取解法を用いた場合に比べて、￨三分程度の反復回数で収束することが分かる.

表1. 反復計算に必要な時間 3 同 の 試 行 の 結 巣 お よ び そ の 平 均 妥 示 し た 方 向 の 人 射 波 に 対 し て 要 し た 時 間 の 手 口 を 示 し て い る . 単 位 は 秒 . BiCGST  AB ，  BiCGST  AB  2 ，  GP‑BiCG などのクリロフ部分空間法を用いた場合， ， J a c o b i 法 や 準 G a u s s ‑ S e i d e l 法を用い た場合に比べて、子分程度の時間で収束することが分かる.

解法

J a c o b i   準 G a u s s ‑ S e i d e l BiCGSTAB  B i CGSTAB2  GP‑BiCG  T r i a l  1  1246  1589  4 2 1   386  403  Trial2  1 2 7 1   1817  482  413  422  Trial3  1305  1826  497  434  462  Average  1274  1744  467  4 1 1   429 

」 一 一 一 一 一 一

ハυー

〉 j i 忽j 岳志・歌凶久司

Mxx ( r e a l )  10s  0 . 8  

0 . 6   0 . 4   0

2 0 

幽

0 . 2

‑ 0

4

‑ 4 0  

³⁰

20

10 9  10  20  30  40  Distance (km) 

陥 xx(ima 思 ) 10s  0 . 2  

0.15  0 . 1   0.05 

0 

幽

0.05

‑ む 圃 1

国

0.15

画 立 2

贋

40

Distance (km) 

図 3 悶 1 のモデルに対する， y‑Q の断商における地織気変換関数の x 成分の言 i 狩 結 果 周 期 1 0 秒の結果を示す.横軸にエ}判票 の{底縦

地に地磁気変換問数の値を示す.上の凶が実剖5 ， ドの凶が虚部を表すー十字がも N

nnamaker ( 1 9 9 1 ) で得られた結巣，

実線が本研究で開発したコードによるもの. I f f i j 存は良い」致そ示す.

周期 10 秒での y 二 O の ~JTI~î ( 間 l の点線)での地織気変換

i 理数の実部および虚部を去す.また， 1 ' ; ( 1   1 1   I こ席 1 0 0 の見掛け比抵抗および佼相のりおよび μ 成分を示す島ど ちらの凶からも，本計での結果と Wannamaker( 1 9 9 1 )   との結果が良 L トー・致を示すことが分かる.

まとめと今後の課題

は ， A  vdeev  e t  a   ( . l 1 9 9 7 ) によって提案された IDM 法に慕づく電感応答計算アルゴリズムそ改良し，得

， J a c o b i 法ではなくクリロツ部分?主問

、た場合の振舞いについて調べた，クリロフ部 分空間法を用いた場合は IDM 法の平Ij点である，どのよう な比抵抗構造に対しでも対応する級数が収束し，原理的に 解が求まるという利点は失われるが，モデル計算ではクリ

ロフ部分空間法を用いても収束は達成され，また収束速度 もほぼ二倍程度に改善されることが分かった.クリロフ部 分空間法の中では，収束回数，計算時間の面でわずかに BiCGSTAB2 法が勝った.

ソースが一般に公開され，良く三次元計算の際に用いら れる Mackie e t   a l .   ( 1 9 9 4 ) による差分法のコードが平面波 入射の応答しか計算できないのに対し，本研究で開発され たコードはあらゆるソースに対する応答を計算できるとい

う 手 J I 点がある.ダイポールソースに対する応容が計算でき ると， Jacobian を相反定理により計算することが可能に なるため，二

次元インパージョンの際に差分で J ζ l c o b i a n

る場合に比べて計算時間を人;帽に短縮することが 可能である.その上，積分方程式法はモデリン夕、領域が電 ある場所だけに眠られるため，先見的な情 報により電気伝導度異常が与在する領域の広さを押えるこ

したコードは，二次元インパージョン L : で有利であると言える.

本 研 究 で は ， 線 形 碍 式 を 解 く 際 前 処 1 寝をJlJl、なかっ た.これは前述、の通り積分) j 程式法における係数行ダ J I は密 であり適切な前処理行列が無いこと，また IDM 法を用い ることで係数行列の性質が改普されると期待されることな どが理由である. しかし，実際さまざまな前処理法を適用 してどの程度改善が見られるのかテストすることは必要で、

あろう. また， 1 7 1 J えぱ Habashy e t  a   . l ( 1 9 9 3 ) ，  T o r r e s ‑

Verdin and Habashy ( 1 9 9 4 )で示された Nonlinears c a t ‑

t e r i n g  approximation のような，式 ( 6 ) の近似解が前処理

行列として利用できる可能性があり，さらに研究を進める

必要がある.

積分 ) j i J ' J 式法を用いた高迷三次元電磁)， t 、終計算コ ドの開発 1 1  

Rhoxy 100s  100 む

100  10 

0 . 1  

。 ^{i s t a n c e( k m )}  

Phase xy 100s 

(印

刷

w

立 ︒ の の こ 仏

30  40  20 

10  O 

胆

10

‑ 2 0  

‑ 3 0  

値

4 0

D i s t a n c e  ( k m )  

( E E C

‑ E O )

Rhoyx 100s 

ハ unu

n v  

dEEE 

100  10 

な 1

‑ 4 0  

30

20

・

υ O 

^U 20  30  40 

D i s t a n c e  ( k m )   Phase yx 100s 

40  20  30 

10 

‑ 1 0   O 

ー 20

‑ 3 0  

剛

4 0

D i s t a n c e  ( k m )  

図 4 関 l のモデルに対する ， y

O の断面におけるインピーダンステンソルの非対角成分の計算結果.周期 1 0 0 秒の結果を示す.

横軸に Z 座標の値，縦剥

には上から順に見掛け比低抗のり成分，位相のり成分，見掛け比抵抗の yx 成分，位相の yx 成分を示

す.記号は図 3 と￨司じ.図中の単位 C O h r n ‑ r n ) は Q'rn を表す 図 3 と同様， Wannarnaker  ( 1 9 9 1)の結果と本研究の結果は良

L  、一致を示す.

1 2   宗包浩志・歌回久可

文 献

Anderson ，  W.  L . ，   1 9 8 2 ，  F a s t  Hankel transforms u s i n g  r e l a t e d   and lagged c o n v o l u t i o n ，  ACM  T r a n s .  on Mαt h .   S o f t ω αr e ， 8 ，  3 4 4 ‑ 3 6 8 .  

Anderson ，  W.  L . ，   1 9 8 9 ，  A hybrid f a s t  Hankel transform a l g o ‑ rithm f o r  electromagn

t i cmodeling ，  G e o p h y s i c s ，  5 4 ，  2 6 3 ‑ 2 6 6 .   Avdeev ，  D .   B . ，   A .  V .   Kuvshinov ，  O .  V .   Pankratov and G .  A .   Newman ，  1 9 9 7 ，  High

performancet h r e e ‑ d i m e n s i o n a l  e l e c ‑ tromagnetic  modelling  u s i n g   m o d i f i e d   Neumann s e r i e s .   Wide‑band numerical s o l u t i o n  and examples ，  J .   Geomag. 

G e o e l e c t r . .  4 9 .  1 5 1 9 ‑ 1 5 3 9  

Barr

t ，R . ，   M. Berry ，  T .  F .   Chan ，  J .   Dammel ，  J .  M. Donato ，  J .   Dongarra ，  V .  E i j k h o u t ，  R .  Poso ，  C .  Romine and H .  Van d e r   V o l s t ，  1 9 9 4 ，  Templates f o r  t h e  s o l u t i o n  o f  l i n e a r  systems :  B u i l d i n g  b l o c k s  f o r  i t e r a t i v e  methods ，  2nd

d i t i o n ， SIAM. 

Chave ，  A .  D . ，  1 9 8 3 ，  Numerical i n t e g r a t i o n  o f  r e l a t e d  Hankel  t r a n s f o r m s  by quadrature and continued f r a c t i o n  expan‑

s i o n ，  G e o p h y s i c s ，  4 8 ，  1 6 7 1 ‑ 1 6 8 6 .  

Fainberg ，  E .  B .  and B .  S h .  S i n g e r ，  1 9 8 0 ，  E l e c t r o m a g n e t i c  induc  t i o n  i n   a  nonuniform s p h e r i c a l  model o f  t h e  Earth ，  Ann  G e o p h y s . ，  3 6 ，  1 2 7 ‑ 1 3 4 .  

藤野清次・張紹良， 1 9 9 6 ，反復法の数理，朝倉書庖， 1 4 0 頁.

Habashy ，  T .  M. ，  R .   W. Groom and  B .   S p i e s ，  1 9 9 3 ，  Beyond t h e   Born and Rytov approximations :  A n o n l i n e a r  approach t o   e l e c t r o m a g n e t i c  s c a t t e r i n g ，  J .   G e o p h y s .  R e s . ，  9 8 ，  1 7 5 9 ‑ 1 7 7 5 .   J i r a c e k ，  G .   R . ，   1 9 9 0 ，  N e a r ‑ s u r f a c e  and t o p o g r a p h i c  d i s t o r t i o n s  

i n  e l

ctromagnetici n d u c t i o n ，  S u r v .  i n  G e o p h y s . ，  1 1 ，  1 6 3 ‑ 2 0 3 .   Macki

R し T . R . Madden and P . E .   Wannamake r .   1 9 9 3  

Three‑dimensional m a g n e t o t e l l u r i c  modeling u s in g  d i f f e ‑ r e n c e   equations‑Theory  and  comparisons  t o   i n t e g r a l   e q u a t i o n  s o l u t i o n s ，  Geo ρ h y s i c s ，  5 8 ，  2 1 5 ‑ 2 2 6 .  

Mackie ，  R .  L . ，   J .  T .   Smith and T .  R .   Madden ，  1 9 9 4 ，  T h r e e ‑ d i m e n s i o n a l  e l e c t r o m a g n e t i c  modeling u si n g  f i n i t e   d i f f e ‑ r e n c e  e q u a t i o n s  :  The magnetot

l l u r i cexample ，  Radio S c i ‑ e n c e ，  2 9 ，  9 2 3 ‑ 9 3 5 .  

Mogi ，  T . ，  1 9 9 6 ，  Three‑dim

n s i o n a lmodeling o f  m a g n e t o t e l l u ‑ r i c  d a t a  u s i n g  f i n i t e  element method ，  J .   A p p l .  G e o p h y s . ，  3 5 ，  1 8 5 ‑ 1 8 9  

Pankratov ，  O .  V . ，  D .   B .   Avdeev and A .  V .   Kuvshinov ，  1 9 9 5 ， 

E l e c t r o m a g n e t i c  f i e l d  s c a t t e r i n g  i n  a  heterogeneous e a r t h  :  A  s o l u t i o n  t o  t h e  forward p r o b l

m ， P h y s .  S o l i d  E a r t h ，  3 1 ， 2 0 1  

‑ 2 0 9 .  

Smith ，  J .  T . ，  1 9 9 6 a ，  C o n s e r v a t i v e  modeling o f  3 ‑ D  e l e c t r o m a g ‑ n e t i c  f i e l d s ，  P a r t  1  :  P r o p e r t i e s  and e r r o r  a n a l y s i s ，  G e o p h y s

i c s ，  6 1 ，  1 3 0 8 ‑ 1 3 1 8 .  

Smith ，  J .  T . ，  1 9 9 6 b ，  C o n s e r v a t i v e  mod

l i n go f  3 ‑ D  e l e c t r o m a g ‑ n e t i c  f i e l d s ，  P a r t  I I   :  B i c o n j u g a t e  g r a d i e n t  s o l u t i o n  and a c ‑ c e l e r a t o r ，  G e o p h y s i c s ，  6 1 ，  1 3 1 9 ‑ 1 3 2 4 .  

T o r r e s ‑ V e r d i n ，  C .   and  T .  M.  Habashy ，  1 9 9 4 ，  Rapid  2 . 5   d i m e n s i o n a l  forward modeling and i n v e r s i o n   v i a  a new  n o n l i n e a r  s c a t t e r i n g  approximation ，  Radio S c i e n c e ，  2 9 ，  1 0 5 1  

1 0 7 9 .  

Wannamaker ，  P . E . ，  1 9 9   ， 1 Advances i n  t h r e e ‑ d i m e n s i o n a l  mag‑

n e t o t e l l u r i c  modeling u s i n g  i n t e g r a l  equations ，  G e o p h y s i c s ，  5 6 ，  1 7 1 6 ‑ 1 7 2 8 .  

Ward ，  S . H .  and G .  W. Hohmann ，  1 9 8 8 ，  Electromagnetic t h e o r y   f o r  g e o p h y s i c a l  a p p l i c a t i o n s ，  i n Electromagnetic methods  i n  a p p l i e d  g e o p h y s i c s " ，  e d i t e d  by Nabighian ，  M. N .  eds ，  S o c i

e t y  o f  E x p l o r a t i o n  G e o p h y s i c i s t s ，  1 8 3 ‑ 3 1 1 .  

Weaver ，  J .  T . ，  1 9 9 4 ，  Mathematical  m

thods f o r   g e o ‑ e l e c t r o m a g n e t i c  i n d u c t i o n ，  John Wiley  &  Sons ，  3 1 6  p p .   Weidelt ，  P . ，  1 9 7 5 ，  E l e c t r o m a g n e t i c   i n d u c t i o n   i n   t h r e e ‑

dimensional s t r u c t u r e s ，  J .   G e o p h y s . ，  4 1 ，  8 5 ‑ 1 0 9  

Xiong ，  Z .  and A .  C .  Tripp ，  1 9 9 5 ，  Electromagnetic s c a t t e r i n g  o f   l a r g e  s t r u c t u r e s  i n  l a y e r

de a r t h  u s i n g  i n t e g r a l  e q u a t i o n s ，  Radio S c i e n c e .  3 0 .  9 2 1 ‑ 9 2 9  

Zhdanov ，  M. S .   and S .  Fang ，  1 9 9 7 ，  Q u a s i ‑ l i n e a r  s

r i e si n  t h r e e ‑ d i m e n s i o n a l  e l e c t r o m a g n e t i c  modeling ，  Radio S c i e n c e ，  3 2 ，  2 1 6 7 ‑ 2 1 8 8 .  

Zhdanov ，  M.S. ，  1 . M. Varentsov ，  J . T .  Weaver ，  N.G. Golubev  and V .  A .  Krylov ，  1 9 9 7 ，  Methods f o r  modeling e l e c t r o m a g ‑ n e t i c   f i e l d s ‑ r e s u l t s   from  COMMEMI‑the i n t e r n a t i o n a l   p r o j e c t  on t h e  comparison o f  modeling m

thodsf o r  e l e c t r o ‑ magnetic i n d u c t i o n ，  J .   A p p l .  G e o p h y s . ，  3 7 ，  1 3 3 ‑ 2 7 1  

Zyserman ，  F . 1 .   and J .  E .   Santos ，  2 0 0 0 ，  P a r a l l e l  f i n i t e  element 

a l g o r i t h m   with  domain  decomposition  f o r   t h r e e ‑

dimensional m a g n e t o t e l l u r i c  modeling ，  J .   App   l . Geo ρ h y s . ， 

4 4 ， 3 3 7 ‑ 3 5 1 .  

積分方程式 j 去を用いた尚速三次

I 主総応答計算コ…ドの i 羽発 1 3  

Appendix 

グ リ ー ン テ ン ソ ル の 計 算 法

取に ; J ( 平成国情 i 左大地におけるグリーンテン 凶の…般解は次のよう ソ

レ の 計 部 方 法 を ま と め る 般 に

層

‑1内部での電場お

よび礎束密度は，ポロイダルポテンシャル I I， およびト 口イダルポテンシャル r を用いて次のように書き表すこと ができる (We

v e r ，1 9 9 4 ) .  

f 子 … ' i 7 x  ' i 7   I I   i w て x r H マ×マ xrlμ t : ' i 7 >

n

f こだし，

川)

r r ( ( )  

である. これらのポテンシャルはヘルムホル ' i 7

nωμど T I 0 

' i 7

r ω μ ど下一一 O を満たす.

w i 境界での1Iと Fの境界条件は，電場および磁場 密度)の水、ド成分が層境界をはさんで連続であるという事

( l l i ) および ( 1 7 ) を書き下してみると， 電 各成分は次のように書き表すことがで きる.

^{: ) ω μ} 山叩

B 2 G 1   ^{i G j 1 1 c (  (} 

従って，電場の水平成分の連続性から a T I ，およびにが述 続であること，また儲場(融束密度)の水平成分の述枕性 から θ i 二およびt;IJ

が連続であることが言える.

次にヘルムホルツ方程式側および凶の一般解を波数領 域で求める.例えば式酬の z 成分は，波数領域で次のよう に表現される.

H  

ャ

j

一

J 1 F  

︐

d

一 {

十 二 Jiw μ 仁 十 十 時 ( 2 5 )   ここで， T I

は T I

の波数成分であり，次のように定義される:

T I

f J I T

( z ) e x p ( i ( k x x

k ， ^{y ) ) d 丸 dk}

凶) ただし，弘および h ^{はそれぞれ}

^{， Y}

^{I j l 白 j の波数を去す.式}

去される

ニご

I I 間領域での T I z 1

) エ変換することによって次のように書き表すことができ る.

包

( H i )  

( 日 1 1 匂二fJ ( i l

十

( 2 8 )  

( 1 8 )  

A および β は対称何より司 λ 二 J(ピ十た~;)のみの関数

である.その場合 式 1 2 印は次のようにハンケル変換の形に

ことカ

で き る 匂a v e r

1 9 9 4 ). 

I I z

f ^え ) e x p ( γ ' Z )

s C A ) εxp( ( 2 9 )  

n叫υ

1  ( 

r  1 3 0 1  

位

。

^{。 1 )}

ただし， < : 1 0 は 0 次のベッセル関数である.式仰の一般解も 同様にしてうえることができる

次にダイボールが存在する出での、/ース項を求める.波 数 領 域 に お い て ダ イ ポ ー ル が 作 る 一 次 電 場 次 密 度 の ? 成 分 は そ れ ぞ れ の ポ ラ ニ ン ン ャ ル の ソ

止を用いて次のように与えられる(式 ( 2 2 ) ，

ニ

am

^{i ω}

^{U ど} m ^{( 3 2 )}  

( 3 3 )   この‑.)，占う}にしでは， トロイ夕、 l V ポテンシャ l V によって 作られる場とボロイ夕、

レポテンシャルによって作られる場 しているので，それぞれの成分を用いて両ポテン シャルのソース E 討を， ダイポールソースが作る一次喝の

を説明するように J 犬めることができる (Ward and Hohmann ，  1 9 8 8 ) 置単位強さの電流ダイポールソース

た時の一次電場， ー次磁京情}支の宅問領域での値 EP ，  sP  ま ( Herlzv c

t o r S を用いて j 欠のようにうえられる

CWeaver ， 1 9 9 4 ) .   EP ニ マ ( ' i 7 ・ . ' 1 ) マ

. ' 1 (

2 )   1

制

倒)

B P

μ t :' i 7   x  . ' 1   ( 3 5 )   ただし，.'1の波数領域での値 s ^{は，ダイポールの￨白 l きの単}

位ベクトルをぬとして波数領域で次のように与えられる.

S=SX;  ( 3 6 )  

S  EEP(‑hztu￨) 

( 3 7 )   4 π t : r 

. z ‑ d i p o l e の場合

式(札 ( 3 5 ) により ム次電場， ふ次磁束密度はそれぞれ

El ご âîS-iωμ~ G 印

1 4   : 4 1 包 j 告;正、 歌田久寸

。

と書き表すことができるー従って、

比較して

( 3 9 )   でうえられる.

( 3 2 ，(坊と ・ ^{y ‑ d i p o l c の場合}

x ‑ d i p o l e の場合のグリーンテンソル成分を時計開りに 9 0 L ， 符 台 を 入 れ 替 え た も の が の 場 合 の グ リーンテンソル成分になるので省略.

ダイポールが存 u : する層でのグリーン関数のソース虫が求 まると，先に求めたポテンシャルの層境界での連続条件，

また地表面，以下隔でのポテンシャルの境界条件を)泣いて 任 意 の 地 点 で の グ リ ー ン 出 数 を 言 す る こ と が で き る 1 9 9 4 :   A ) .   数値三十昨によりグリーン 関 数 を 求 め る 際 に 必 要 と な る ハ ン ケ ル 変 換 に は C h a v ε ( の i 丘陵積分?去を用いた.この方 i 法は， cODvolutioD  f i l t c r を用いた Anderso

( の 方 法 に 比 べ て 計 第 時 間

る代わりに誤ぷが予測できること，ならびにより条

( 1 1 0 )   exp( 示 ‑w i )

4 πど γ

n  0 

・ を得る. ^{x ‑ d i p o l e の場合} 闘により

E~= δ)å4 S

ー

!A

ま/町

一次磁束密度はそれぞれ

である.

( 4 2 )   ( i 1 3 )   ( 3 3 ) と上ヒ較して

( 判

r ;  

を得る.ただし司 t ( ; i

l l Z   ^?と o

1  z‑w<O 

( 4 5 )   の : コ

チンを出用し ディングでは Andcrson ( 1 9 8 9 )  

( 4 6 )