重力および磁気異常解析のための２次元フィルタに関する一考察

重力および磁気異常解析のための２次元フィルタに関する一考察

Discussion on the two-dimensional filter operations for gravity and magnetic anomalies

Ryuji Kubota*

＋2016年11月4日受付，2017年8月30日受理

＊川崎地質株式会社，〒108-0073 東京都港区三田2-11-15

 Kawasaki Geological Engineering Co., Ltd., 2-11-15 Mita, Minato-ku, Tokyo, 108-0073, Japan

ABSTRACT

In this paper, we discuss expressions of potential fields in the frequency domain and filter operations by the two dimensional Fast Fourier Transform method in interpretations of the gravity and the magnetic anomalies. Especially by using the complex Fourier series, the representations in the frequency domain, which correspond to differentiation, integration and continuation in the spatial domain, are simplified, and the meanings of the coefficients used in the two-dimensional filter operation are clarified.

Some useful filter operations are described, which are the pseudo-gravity, reduction to the pole, vector magnetic anomalies derived from the total magnetic field, upward continuation， three dimensional density and magnetization structures confined by horizontal planes and analytic signals.

A simple new iterative method for analyzing double-layer structures from gravity (or pseudo-gravity) data is also proposed, in which an expression in the frequency domain by Parker (1972) is used as a model function. Either of the subsurface structure or the two-dimensional variation of the density (or the magnetization) can be obtained by this method.

１．はじめに  テクトニクスや地震・火山の研究，あるいは海底資源 探査において，重力・地磁気の観測は極めて重要である ことは言うまでもない．特に海上観測や自律型無人潜水 機により均質で大量のデータが蓄積されつつあり，海洋 域の詳細な重力異常や磁気異常の分布が得られるように なってきた．また火山島の研究においても航空磁気測量 が繰り返し実施され，噴火予知のための研究に大きく貢 献している．したがってこれらの観測結果から，面的に あるいは３次元的に地下構造を把握する方法を議論して おくことは極めて重要である．  重力異常と磁気異常はともにラプラスの方程式を満足 するポテンシャル関数から導かれる量であることから， それらの解析法には共通点が多い．これらの観測値から 地下の地質構造を量的に推定するための方法が数多く考 案されている．代表的なものとしては，あらかじめ地下 構造の幾何学的模型を想定して，それから得られる重力 異常や磁気異常の分布を観測値と比較し，最も適合する と思われるモデルを探し出す方法（例えばTalwani and Ewing, 1960）や，地表あるいは海上の観測値から，下 方接続によって直接地下構造を推定する方法などがある （例えばTsuboi and Tomoda, 1958; Saito and Takeuchi,

1964）．また，適当なフィルタによって対象としている 地域の局所的，あるいは広域的構造を把握する方法など がある．  重力異常や磁気異常の解析を目的としたフィルタの種 類は数多くあり，鉛直微分，上・下方接続，磁気異常の 擬似重力および極磁力変換などが代表的である．これら の計算法には大別して二通りあり，格子状に分布した観 測値に，あらかじめ定められた半径あるいは正方領域内 の格子点について，それぞれ異なる重み係数をかけて和 を計算する方法（例えばStencil法，Rosenbach, 1953; Baranov, 1957），ラプラスの方程式の解を有限フーリエ 級数の形に展開して，これらの量を計算する方法である （例えば，加藤，1987）．  特に加藤（1987）の教科書には，２次元の空間フィル タ理論が体系的にまとめられている．これは，観測され た重力異常または磁気異常の分布を２次元の有限フーリ エ級数に展開し，施すべきフィルタの波数領域での係数 を４行４列のマトリックスで表現し，両者の積をとるも のである．フィルタの特性が空間領域での微分に対応し ているとき，そのマトリックスを微分作用子，積分に対 応するときを積分作用子，上・下方接続に対応するとき を接続作用子とよんでいる．また，駒澤（1980）はやは り有限フーリエ級数を用いて，反復法により重力異常か ら３次元の地下構造を解析する方法を示している．  一方，高速フーリエ変換（Fast Fourier Transform: FFT）は，等間隔な標本点の値について，有限複素フ ーリエ級数を用いてその係数を求めるものである（例え ば，大崎，1994）．この有限複素フーリエ級数を用いる ことによって，空間領域での微分，積分および接続に対 応する波数領域での表現が極めて簡潔になり，フィルタ 操 作 で 使 わ れ る 係 数 の 意 味 が 明 確 に な る． す で に Blakely （1995）の教科書には，ポテンシャル場に対す るいくつかの有用なフィルタの表現式とFortran 77によ るプログラムが示されている． 

 本論ではBhattacharyya and Navolio （1976）による 空間領域と波数領域でのポテンシャル場の表現式を基 に，波数表現による擬似重力および極磁力，全磁力異常 の３成分変換，上方接続とパワースペクトル，３次元平 板モデルの磁化と密度分布，複素勾配法について述べ る．さらに地下の起伏分布を波数領域でテーラー展開す るParker （1972）の表現式をモデル関数として用いる， 反復法による構造解析の方法を示す．また参考のため に，２次元高速フーリエ変換と極磁力のFortran77プログ ラムの例を付録に示す．なお本論文の一部は，過去３回 の海洋調査技術学会研究成果発表会（1989, 1997, 2000） にて発表したものである． ２．２次元フーリエ変換のまとめ  いまݔ െ ݕ_{平面上に定義された空間を考え，それぞれの} 方向に対応する波数をݑ, ݒとする．このとき関数݂(ݔ, ݕ)の フーリエ変換をܨ(ݑ, ݒ)とすれば， ݂(ݔ, ݕ) = ඵ ܨ(ݑ, ݒ) ݁ݔ݌[݅(ݑݔ + ݒݕ)] ݀ݑ݀ݒ, ஶ ିஶ ܨ(ݑ, ݒ) = 1 4ߨଶඵ ݂(ݔ, ݕ) ݁ݔ݌[െ݅(ݑݔ + ݒݕ)] ݀ݔ݀ݕ ஶ ିஶ (2.1) である．ここに݅ = ξെ1である．この関係をここでは記 号的に݂(ݔ, ݕ) ՞ ܨ(ݑ, ݒ)と書く．記号՞は空間領域におけ る݂(ݔ, ݕ)と波数領域におけるܨ(ݑ, ݒ)とが対をなし，互い にフーリエ変換あるいはフーリエ逆変換の関係にあるこ とを表す．実際には，ݔ െ ݕ平面上での有限個の離散標 本値を扱うことから，式（2.1）を有限複素フーリエ級 数になおす．いまݔ方向およびݕ方向の標本値の数をܯ, ܰ とすると， ݂(݇, ݈) = ෍ ෍ ܨ(݉, ݊)݁ݔ݌ [2ߨ݅ ൬݇݉ ܯ+ ݈݊ ܰ൰] ேିଵ ௡ୀ଴ ெିଵ ௠ୀ଴ ,

ܨ(݉, ݊) = 1 ܯܰ෍ ෍ ݂(݇, ݈)݁ݔ݌ [െ2ߨ݅ ൬ ݇݉ ܯ+ ݈݊ ܰ൰] ேିଵ ௟ୀ଴ ெିଵ ௞ୀ଴ , (2.2) k,m=0,1,2,…㸪 M-1, l,n=0,1,2,…㸪 N-1 となる．式（2.2）の計算では，ܯ, ܰの数が大きくなる と膨大な計算時間が必要となるが，高速フーリエ変換に よって短時間に計算することができる．  つぎに与えられた関数݂(ݔ, ݕ)の空間領域での微分と積 分，さらにそれに対応する波数領域での表現を整理して おくこととする．すなわち微分に関しては，݂(ݔ, ݕ)が連 続であれば， ߲݂(ݔ, ݕ) ߲ݔ ՞ ݅ݑܨ(ݑ, ݒ), ߲௡_{݂(ݔ, ݕ)} ߲ݔ௡ ՞ (݅ݑ)௡ܨ(ݑ, ݒ), (2.3) ߲ଶ_{݂(ݔ, ݕ)} ߲ݔ߲ݕ ՞ െݑݒܨ(ݑ, ݒ). 積分に関しては， න න ݂(ݔ, ݕ)݀ݔ݀ݕ ՞ െܨ(ݑ, ݒ)_ݑݒ (2.4) 㹷 ିஶ 㹶 ିஶ が成り立つ．   付 録Aに ２ 次 元 高 速 フ ー リ エ 変 換 の プ ロ グ ラ ム 例 FAST2 を示す．この例では大崎（1994）の１次元高 速フーリエ変換のプログラム FAST を用いて，まず x方向にフーリエ変換を施し，その結果を続いてy方向 にフーリエ変換している． ３．空間領域におけるポテンシャル場  ここではまず全磁力異常について考える．いまݖを上 向きとする直交座標系で，ݖ < 0 _{に均一に磁化した物体} を考えると，ある観測点(ݔ, ݕ, ݖ)における磁気ポテンシャ ルܸ(ݔ, ݕ, ݖ)はつぎのようになる（Bhattacharyya and Navolio, 1976）. ܸ(ݔ, ݕ, ݖ) = ܵ(ݔ, ݕ, ݖ) כ ቐ݉௢(ܮݔ + ܯݕ + ܰݖ) (ݔଶ_{+ ݕ}ଶ_{+ ݖ}ଶ₎ଷଶ ቑ . (3.1) ここでܵ(ݔ, ݕ, ݖ)は，磁化している物体の形状関数（Source geometry）である．moは磁化の強さ，L, M, N は磁化 の方向余弦である．*は合積（コンボリューション： Convolution）を表す．式（3.1）の中括弧㹹  㹻内は点源 に対する表現式であり，合積はそれを物体ܵにわたって 積分していることを意味している．  したがって式（3.1）の中括弧㹹  㹻内は，スカラーポ テンシャルに対するグリーン関数である．式（3.1）は また ܸ(ݔ, ݕ, ݖ) = െܵ(ݔ, ݕ, ݖ) כ ቊ݉௢߲݁(ݔ, ݕ, ݖ)_߲ߛ ቋ  (3.2) とも書ける．ここで߲ߛ_{は磁化の方向の微分量に対応し，} ݁(ݔ, ݕ, ݖ) = 1/ඥݔଶ_{+ ݕ}ଶ_{+ ݖ}ଶ_{は点源と観測点との距離の逆数} である．  また全磁力異常ܶ(ݔ, ݕ, ݖ)は，式（3.2）で与えられる磁 気ポテンシャルを地球磁場方向に微分したものであるこ とから， ܶ(ݔ, ݕ, ݖ) = ܵ(ݔ, ݕ, ݖ) כ ቊ݉௢߲ଶ݁(ݔ, ݕ, ݖ) ߲ߛ߲ݐ ቋ (3.3) と表される．ここで߲ݐは，方向余弦݈, ݉, ݊を持つ地球磁 場方向の微分量である．  つぎに，式（3.1）で用いた物体ܵについて，その内部 が 一 様 な 密 度ȡで 構 成 さ れ て い る と き， あ る 観 測 点 (ݔ, ݕ, ݖ) における重力ポテンシャルは, ܷ(ݔ, ݕ, ݖ) = ܩ_௢ߩ[ܵ(ݔ, ݕ, ݖ) כ ݁(ݔ, ݕ, ݖ)]  (3.4) となる．ここにܩ௢は万有引力定数である．したがって重 力異常は， ݃(ݔ, ݕ, ݖ) = ܩ௢ߩ[ܵ(ݔ, ݕ, ݖ) כ ߲݁(ݔ, ݕ, ݖ)/߲ݖ] (3.5) となる． ４．波数領域におけるポテンシャル場  前節において重力ポテンシャルと磁気ポテンシャル は，ともにグリーン関数と物体の形状関数との合積によ って表されることを示した．したがって波数領域では， 重力ポテンシャルと磁気ポテンシャルは，これらの２つ の関数の単なる積となる．いま全磁力異常について，そ の波数領域での表現を考えることにする．  物体の形状関数ܵ(ݔ, ݕ, ݖ)とグリーン関数ܩ(ݔ, ݕ, ݖ)のフー リエ変換を，ݑ, ݒ, ݓを波数，３次元のフーリエ変換の対 を฻で 表 し て，ܵ(ݔ, ݕ, ݖ) ฻ ܣ(ݑ, ݒ, ݓ)㸪 ܩ(ݔ, ݕ, ݖ) ฻ ܤ(ݑ, ݒ, ݓ) とする．このとき式（3.3）について，そのフーリエ変

換ܶி_は

ܶி_{(ݑ, ݒ, ݓ) = ܣ(ݑ, ݒ, ݓ) ή ܤ(ݑ, ݒ, ݓ)   (4.1)}

となる．ただし本章ではBhattacharyya and Navolio （1976）の定義にしたがってフーリエ変換を， ܨ(ݑ) = න ݂(ݔ)݁ି௜௨௫ ஶ ିஶ ݀ݔ, ݂(ݔ) =_2ߨ1 න ܨ(ݑ)݁௜௨௫_݀ݑ ஶ ିஶ と定義する．なおܤ(ݑ, ݒ, ݓ)には磁化の方向，地球磁場方 向，観測点の座標に関する項が含まれている．ܤ(ݑ, ݒ, ݓ) を求めるために，まず݁(ݔ, ݕ, ݖ)についてのフーリエ変換 を考える．݁(ݔ, ݕ, ݖ)のフーリエ変換をܧ(ݑ, ݒ, ݓ)とすれば， ܧ(ݑ, ݒ, ݓ) = ම݁ݔ݌ [െ݅(ݑݔ + ݒݕ + ݓݖ)] (ݔଶ_{+ ݕ}ଶ_{+ ݖ}ଶ₎ଵ/ଶ ஶ ିஶ ݀ݔ݀ݕ݀ݖ = 4ߨ ݑଶ_{+ ݒ}ଶ_{+ ݓ}ଶ (4.2)

となる（Bhattacharyya and Navolio, 1976）.  また式（3.3）の中に現れる方向微分は， ߲ ߲ߛ= ܮ ߲ ߲ݔ+ ܯ ߲ ߲ݕ+ ܰ ߲ ߲ݖ, ߲ ߲ݐ= ݈ ߲ ߲ݔ+ ݉ ߲ ߲ݕ+ ݊ ߲ ߲ݖ (4.3) となる．したがってこれらの方向微分に対応するフーリ エ変換は， ߲݁(ݔ, ݕ, ݖ) ߲ߛ ฻ ݅(ܮݑ + ܯݒ + ܰݓ)ܧ(ݑ, ݒ, ݓ) (4.4) および ߲ଶ_{݁(ݔ, ݕ, ݖ)} ߲ߛ߲ݐ ฻ െ(ܮݑ + ܯݒ + ܰݓ) ή (݈ݑ + ݉ݒ + ݊ݓ) ή ܧ(ݑ, ݒ, ݓ) (4.5) となる．すなわち(ܮݑ + ܯݒ + ܰݓ)および(݈ݑ + ݉ݒ + ݊ݓ)は， それぞれ磁化の方向および地球磁場方向の波数を表して いる．ゆえに 式（3.3）,（4.1）,（4.2） および（4.4）よ り全磁力異常ܶ(ݔ, ݕ, ݖ)のフーリエ変換ܶி_{(ݑ, ݒ, ݓ)は,} ܶி_{(ݑ, ݒ, ݓ) = െ4ߨ݉} ௢(ܮݑ + ܯݒ + ܰݓ) ή (݈ݑ + ݉ݒ + ݊ݓ) ή _ݑଶܣ(ݑ, ݒ, ݓ)_{+ ݒଶ+ ݓଶ}     (4.6) となる．  一方，全磁力異常あるいは重力異常は，ある一定の高 さݖ = ݖଵなる平面上で与えられる．したがって観測され る全磁力異常はܶ(ݔ, ݕ, ݖଵ)であり，ݔおよびݕの２変数につ いての関数となる．またT(x,y,z1)のフーリエ変換も，u とvの２次元平面上で定義される．いまそれをܶ௙_{(ݑ, ݒ)と} すれば，ܶ௙_{(ݑ, ݒ)とܶ}ி_{(ݑ, ݒ, ݓ)の関係は,} ܶ௙_{(ݑ, ݒ) =} 1 2ߨ න ܶி ஶ ିஶ (ݑ, ݒ, ݓ) ݁ݔ݌(െ݅ݓݖଵ) ݀ݓ (4.7) である．またここでݖଵ= 0として結局， ܶ௙_{(ݑ, ݒ) = െ2ߨ݉} ଴(ܮݑ + ܯݒ െ ݅ܰݏ) ή (݈ݑ + ݉ݒ െ ݅݊ݏ) ή ܣ(ݑ, ݒ, െ݅ݏ) ݏ (4.8)

となる（Bhattacharyya and Navolio, 1976）．ただし， s = ξݑଶ_{+ ݒ}ଶ_{である．また同時に重力異常では求めるべき} フーリエ変換を݃௙_{(ݑ, ݒ)とすれば,} ݃௙_{(ݑ, ݒ) = 2ߨܩ} ௢ߩ ή ܣ(ݑ, ݒ, െ݅ݏ) (4.9) となる．ここで導入したݏは，つぎのような意味をもつ． すなわち磁気異常源が観測平面より下方にのみ存在する とき，全磁力異常ܶ(ݔ, ݕ, ݖ)は次のラプラスの方程式を満 たし, ߲ଶ_ܶ ߲ݔଶ+ ߲ଶ_ܶ ߲ݕଶ+ ߲ଶ_ܶ ߲ݖଶ= 0,   (4.10) そのフーリエ変換は， (ݑଶ_{+ ݒ}ଶ_{+ ݓ}ଶ_)ܶி_{(ݑ, ݒ, ݓ) = 0 (4.11)} となる．ここで自明解ܶி_{(ݑ, ݒ, ݓ) = 0以外でも成立するた} め に は，ݑଶ_{+ ݒ}ଶ_{+ ݓ}ଶ_{= 0で な け れ ば な ら な い． こ こ で} ݑଶ_{+ ݒ}ଶ_{= ݏ}ଶ_{とするとw=±݅ݏとなるが，ݖ が大きく（高く）} なるほどポテンシャルは減衰しなければならないから，      ݓ = െ݅ݏ   (4.12)

となる． ５．擬似重力および極磁力  重力ポテンシャルܷ_{と磁気ポテンシャル}ܸ_{の間には，} 次式で示されるポアッソンの関係式が成り立つ． ܸ =݉௢ ܩ௢ߩή ߲ܷ ߲ߛ , (5.1) すなわち，磁気ポテンシャルは重力ポテンシャルの磁化 方 向 へ の 一 次 微 分 に 等 し い（ 例 え ば 物 理 探 査 学 会， 2005）．このことは，得られた全磁力異常に何らかの操 作を施してやれば，同一の物体について期待されるべき 重力を算出できることを示している．このような考えか ら，観測された全磁力異常から重力異常を求める方法 が，Baranov （1957） によって考案され，そのようにし て求められた重力異常を擬似重力（Pseudo-gravity）と よんでいる．また長谷川（1967）は，Baranovの方法を さらに改良した方法を示した．いずれもこれらの方法 は，式（5.1）の積分解から出発する方法であり，適当 な格子状の係数列を求めてStencil 法によって擬似重力 を求めるものであるが，これらの係数列を求めるための 数値計算はきわめて複雑である．  ここでは波数領域で擬似重力を求める方法を考えてみ る．全磁力異常ܶ(ݔ, ݕ, 0)と重力異常݃(ݔ, ݕ, 0)は，式（5.1） をzとtについて微分すると， െμܶ μݐ= ܩ_଴ߩ ݉଴ ߲ଶ_݃ ߲ߛ߲ݐ (5.2) となる．ここで式（4.3）と（4.12）から ߲ ߲ߛ ฻ ݅(ܮݑ + ܯݒ െ ݅ܰݏ), ߲ ߲ݐ ฻ ݅(݈ݑ + ݉ݒ െ ݅݊ݏ) となるから，したがって݃(ݔ, ݕ, 0)のフーリエ変換݃௙_{(ݑ, ݒ)} は， ݃௙_{(ݑ, ݒ) =}ܩ௢ߩ ݉௢ή െݏܶ௙_{(ݑ, ݒ)} (ܮݑ + ܯݒ െ ݅ܰݏ) ή (݈ݑ + ݉ݒ െ ݅݊ݏ) (5.3) となる．  実際の計算は，観測された全磁力異常ܶ(ݔ, ݕ, 0)_を式 （2.2）にしたがって複素フーリエ級数に展開し，ܶ௙_{(ݑ, ݒ)} を求め，式（5.3）の計算を行った後，そのフーリエ逆 変換をとれば，擬似重力を得ることができる．  ここではポアッソンの関係式に現れる空間微分と，そ の波数領域で現れる係数との関係を明確にするために， 以上のような導き出しを試みたが，式（5.3）はポアッ ソンの関係式を経由せずとも，式（4.8）および式（4.9） からも直ちに導かれることがわかる．

 一方Baranov and Naudy （1964）は，擬似重力をさら に鉛直方向に微分することによって，極磁力 （reduction to the magnetic pole） を求める方法を示した．これは 例えば北極では磁場の伏角が90°となることから全磁力 異常の極が一つとなり，擬似重力同様物体の形状あるい は磁化の分布が把握できる．  極磁力ܶோ்௉(ݔ, ݕ, 0)は，重力ポテンシャルを鉛直方向に ２度微分したものとして定義されるので，式（5.3）か ら ܶோ்௉(ݔ, ݕ, 0) ՞ െݏଶ_ܶ௙_{(ݑ, ݒ)} (ܮݑ + ܯݒ െ ݅ܰݏ) ή (݈ݑ + ݉ݒ െ ݅݊ݏ) (5.4) となる．  付録Bに極磁力変換のプログラム RTP を示す．プ ログラム中の２次元データ配列(ݔ, ݕ)は，東方向をx，北 方向をyとし，配列（1, 1）は南西端としている．また 入力データ数は両者とも２の累乗でなくてはならない が，実際にはそのようなことは稀である．通常はデータ 端部にテーパを施し不足分を 0.0 で埋める方法も考え られるが，筆者の場合はフーリエ変換の繰り返し性を考 慮して，データの始点と終点が滑らかに繋がるように適 当な関数を選んで不足分を埋めるようにしている．計算 終了後は元の格子範囲だけを取り出してファイルに格納 する．  以下に計算例を示す．Fig.1にモデル（a）とそれに対す る重力異常（b），全磁力異常（c）を示す．モデルには 50 km四方の中に２個の円錐台が配置されており，いず れも底面の深さを５km，直径10 km，高さ３km，頂上 部の直径を４kmとしている．このモデルを100 mメッ シュの格子点（501×501）で表し，各格子点を角柱で近 似し，同様のメッシュの大きさで観測平面を０mとして それぞれの異常分布を計算した．角柱に対する重力異常 は萩原（1979）を，全磁力異常はBhattacharyya （1964） の表現式をそれぞれ用いた．なお円錐台の密度を1,000 kg/m3_{，磁化の大きさを１A/m，地球磁場および磁化の} 方向をいずれも偏角（Dec.）を０°，伏角（Inc.）を40° とした．Fig.2には，Fig.1（c）に対して式（5.5）およ び（5.6）から求めた擬似重力変換（a），極磁力変換（b） の結果を示す．Fig.2（a）とFig.1（b）を比べると擬似

重力の結果が全体に−５mgalほど下がっていることが わかる．これは波数０での計算が不能となってしまうた め，その複素振幅を０と置いていることによる．すなわ ち，場の平均値を０と仮定しているために起きる現象で ある．またFig.2（b）では地形高まり部の中心に極磁力 分布の極が一致していることがわかる．

Fig. 1. (a) Subsurface structure model with two circular truncated cones of the same form in an area of 50 km x 50 km. The depths of the top and the bottom are 2 km and 5 km, and their diameters are 4km and 10 km,

respectively. (b) Theoretical gravity anomalies by the equation for a prism-shaped body. (c) Theoretical total magnetic anomalies by the equation for a prism-shaped body. ȡ=1,000 kg/m3_{, m}

o=1 A/m, Inc.=40°，Dec.=0° .

The solid lines indicate positive values and the dashed lines, negative values. The symbols H and L indicate higher and lower than surroundings, respectively.

Fig. 2. (a) Pseudo-gravity and (b) Reduction to the pole derived from the total magnetic anomalies of Fig.1(c) by 2D FFT. The symbols H indicate higher than surroundings.

６．全磁力異常から導き出される磁気異常の3成分  全磁力異常から直接磁場の３成分を推定できることを 示したのは，Lourenco and Morrison （1973）である． 全磁力異常の３成分ベクトル൫ܪ௫, ܪ௬, ܪ௭൯は， ܪ௫= െ ߲ܸ ߲ݔ, ܪ௬= െ ߲ܸ ߲ݕ, ܪ௭= െ ߲ܸ ߲ݖ (6.1) で定義され，また式（4.3）および（5.2）から， ܶ(ݔ, ݕ, 0) = െ߲ܸ ߲ݐ= െ ൬݈ ߲ ߲ݔ+ ݉ ߲ ߲ݕ+ ݊ ߲ ߲ݖ൰ ܸ (6.2) なので，全磁力異常ܶ(ݔ, ݕ, 0)は， ܶ(ݔ, ݕ, 0) = ൫݈ܪ_௫+ ݉ܪ_௬+ ݊ܪ_௭൯ (6.3) で表される．ここで݈，݉，݊は地球磁場方向の方向余弦 である．したがって式（6.2）のフーリエ変換表示は， 式（2.3）から ܶ௙_{(ݑ, ݒ) = െ݅(݈ݑ + ݉ݒ െ ݅݊ݏ)ܸ}௙_{(ݑ, ݒ) (6.4)} となる．またܸ௙_{はܸࡢのフーリエ変換である．したがっ} て， ܸ௙_{(ݑ, ݒ) =} െܶ௙(ݑ, ݒ) ݅(݈ݑ + ݉ݒ െ ݅݊ݏ) となるから，求める൫ܪ௫, ܪ௬, ܪ௭൯は ܪ_௫(ݔ, ݕ, 0) ՞ െ݅ݑܸ௙_{(ݑ, ݒ) =} ݑܶ௙(ݑ, ݒ) ݈ݑ + ݉ݒ െ ݅݊ݏ , ܪ௬(ݔ, ݕ, 0) ՞ െ݅ݒܸ௙(ݑ, ݒ) = ݒܶ௙_{(ݑ, ݒ)} ݈ݑ + ݉ݒ െ ݅݊ݏ , ܪ_௭(ݔ, ݕ, 0) ՞ െݏܸ௙_{(ݑ, ݒ) = െ} ݅ݏܶ௙(ݑ, ݒ) ݈ݑ + ݉ݒ െ ݅݊ݏ (6.5) となる．  Fig.3（a）に角柱の式から求めた磁気異常の３成分を， （b）に式（6.5）よりFig.1（c）の全磁力異常から求めた それを示す．矢印は水平２成分によるベクトルを，赤は 鉛直成分が正，青は負を意味する．両者がほぼ一致して いることがわかる．

Fig. 3. (a) Theoretical vector magnetic anomalies by the equation for a prism-shaped body. (b) Vector magnetic anomalies derived from the total magnetic anomalies of Fig.1(c) by 2D FFT. Red and blue arrows indicate that vertical components are positive and negative, respectively.

７．上方接続とパワースペクトル  つぎに高さݖ = 0で与えられる全磁力異常ܶ(ݔ, ݕ, 0)（あ るいは重力異常）と，高さݖ = ݖଶで与えられるܶ(ݔ, ݕ, ݖଶ)の 関係を考えてみる． ݖଶ> 0, ܶ(ݔ, ݕ, 0) ՞ ܶ௙(ݑ, ݒ), ܶ(ݔ, ݕ, λ) = 0 (7.1) とする．いまܶ(ݔ, ݕ, ݖଶ)࡟についてݔ െ ݕ平面に関するフーリ エ変換をܶכ_{(ݑ, ݒ, ݖ} ଶ)とすれば，式（2.1）と同様に ܶ(ݔ, ݕ, ݖଶ) = ඵ ܶכ(ݑ, ݒ, ݖଶ) ݁ݔ݌[݅(ݑݔ + ݒݕ)] ݀ݑ݀ݒ , ஶ ିஶ (7.2) である．したがってܶ(ݔ, ݕ, ݖଶ)のz方向の微分は ߲ܶ(ݔ, ݕ, ݖଶ) ߲ݖ = ඵ ݀ ݀ݖܶכ(ݑ, ݒ, ݖଶ) ݁ݔ݌[݅(ݑݔ + ݒݕ)] ݀ݑ݀ݒ, (7.3) ஶ ିஶ となるから式（2.3）にならって, ߲௡_{ܶ(ݔ, ݕ, ݖ} ଶ) ߲ݖ௡ ՞ ݀௡ ݀ݖ௡ܶכ(ݑ, ݒ, ݖଶ) (7.4) となる．  ここで再び，ラプラスの方程式（4.10）のフーリエ変 換を，式（7.4） を使って式（4.11）を書き改めると, ቆെݑଶ_{െ ݒ}ଶ₊ ݀ଶ ݀ݖଶቇ ܶכ(ݑ, ݒ, ݖଶ) = 0 または ݀ଶ_ܶכ_{(ݑ, ݒ, ݖ} ଶ) ݀ଶ_ݖ െ (ݑଶ+ ݒଶ)ܶכ(ݑ, ݒ, ݖଶ) = 0 (7.5) となる.式（7.1）の境界条件を考慮すると上式は ܶכ_{(ݑ, ݒ, ݖ} ଶ) = ܶ௙(ݑ, ݒ)݁ݔ݌ ቀെݖଶඥݑଶ+ ݒଶቁ , (7.6) と解くことができる．このことは，ܶ(ݔ, ݕ, 0)が既知なら ばܶ(ݔ, ݕ, ݖଶ)も既知となることを示している．したがって ܶ(ݔ, ݕ, ݖଶ) ՞ ܶ௙(ݑ, ݒ)݁ݔ݌ ቀെݖଶඥݑଶ+ ݒଶቁ , (7.7) である．式（7.7）による操作が上方接続とよばれるも のである．

Fig. 4. (a) Upward continuation of 2 km by the equation (7.7) for the total magnetic anomalies of Fig.1(c). (b) Theoretical total magnetic anomalies at the altitude of 2 km for the model of Fig.1(a). The symbols H and L indicate higher and lower than surroundings, respectively.

 つぎにܶ(ݔ, ݕ, ݖ)のパワースペクトルについて考える． 式（7.6）からパワースペクトルは |ܶכ_{(ݑ, ݒ, ݖ} ଶ)|ଶ= หܶ௙(ݑ, ݒ)หଶ݁ݔ݌ ቀെ2ݖଶඥݑଶ+ ݒଶቁ , (7.8) となる．ここで|ܶכ_{(ݑ, ݒ, ݖ} ଶ)|ଶ=P(u,v)とし，またหܶ௙(ݑ, ݒ)หଶを 白色スペクトルとすればหܶ௙_{(ݑ, ݒ)ห}ଶ_{=cܿ݋݊ݏݐ.であるから，定} 数cを用いて式（7.8）はlog௘ܲ(ݑ, ݒ) = ܿଶെ 2ݖଶξݑଶ+ ݒଶとな る．したがってある高さで得られた全磁力異常分布につ い て パ ワ ー ス ペ ク ト ル を 求 め，log_௘ܲ(ݑ, ݒ)_{を 縦 軸 に，} ξݑଶ_{+ ݒ}ଶ_{を横軸にとったグラフの勾配から}ݖଶの値（White Depth）が推定できる．物体上では磁気異常が白色スペ クトルを持つので，ݖଶの値は物体の深さに相当する．  Fig.4（a）はFig.1（c）を+２kmに上方接続した結果 であり，（b）はFig.1（a）のモデルについて，+２km上 方において角柱の式から計算した結果である．両者がほ ぼ一致していることがわかる． ８．３次元平板モデルの磁化と密度分布  地磁気の研究において，過去の地球磁場逆転の問題を 取り扱うために，現在の地球磁場方向の影響を取り除い て，過去の磁場の方向を正負の磁化の大きさとして表現 する方法がある（Schouten and McCamy, 1972）．この 方法は，地下に一定の厚さを持つ２次元平板を仮定し， 線形フィルタ理論によって平板内の磁化の分布を求める 方法である．ここではこの方法を３次元の平板問題に拡 張するとともに，重力に対する同様な３次元平板の密度 分布を求める方法について議論する．  まず磁化の分布の問題を考える．いま磁化している層 は一定の深さおよび厚さを持つと考え，磁気異常はその 層の磁化の強さの変化に起因していると考える．すなわ ち式（3.1）の݉௢をݔとݕ)の関数݉௢(ݔ, ݕ)と考え，そのフ ーリエ変換を݉_௢௙(ݑ, ݒ)とする．また式（4.8）のܣ(ݑ, ݒ, െ݅ݏ) をSchouten and McCamy （1972）にしたがって，磁化 し て い る 層 の 上 面 お よ び 下 面 深 度 をa，bと し て，

ܣ(ݑ, ݒ, െ݅ݏ) = ݁ି௔௦_{െ ݁}ି௕௦_{とすると，式（4.8）は}

ܶ௙_{(ݑ, ݒ) =}െ2ߨ݉଴௙(ݑ, ݒ)

ݏ (ܮݑ + ܯݒ െ ݅ܰݏ) ή (݈ݑ + ݉ݒ െ ݅݊ݏ) ή (݁ି௔௦_{െ ݁}ି௕௦_{) (8.1)}

と な る． 式（8.1） の 各 項 をSchouten and McCamy （1972） の 結 果 と 比 較 す る と，2ߨ(݁ି௔௦_{െ ݁}ି௕௦_)はearth filterに，(ܮݑ + ܯݒ െ ݅ܰݏ)(݈ݑ + ݉ݒ െ ݅݊ݏ)/ݏはphase filterに対 応することがわかる．これより求める݉௢(ݔ, ݕ)はフーリ エ変換を用いて ݉௢(ݔ, ݕ) ՞ െܶ௙_{(ݑ, ݒ) ή ݏ} 2ߨ(ܮݑ + ܯݒ െ ݅ܰݏ) ή (݈ݑ + ݉ݒ െ ݅݊ݏ) ή (݁ି௔௦_{െ ݁}ି௕௦₎ (8.2) となる．  同様に密度分布の問題も同じように扱える．すなわち 密度分布をߩ(ݔ, ݕ)として，そのフーリエ変換をߩ௙_{(ݑ, ݒ)と} すると，式（4.9）は ݃௙_{(ݑ, ݒ) = 2ߨܩ} ௢ߩ௙(ݑ, ݒ) ή (݁ି௔௦െ ݁ି௕௦) (8.3) となる．したがって求める密度分布は， ߩ(ݔ, ݕ) ՞ ݃௙(ݑ, ݒ) 2ߨܩ௢(݁ି௔௦െ ݁ି௕௦) (8.4) となる． ９．複素勾配法  複素勾配法（analytic signal）は，ポテンシャル場デ ータの水平勾配と鉛直勾配がヒルベルト変換の関係で結 ばれていることから，その勾配を複素数で表現し，さら にその振幅（絶対値）を求める方法である（物理探査学 会，2005）．いま全磁力異常の場合を考えると，その複 素勾配Aは， ܣ =߲ܶ(ݔ, ݕ, 0) ߲ݔ ଙԦ + ߲ܶ(ݔ, ݕ, 0) ߲ݕ ଚԦ + ݅ ߲ܶ(ݔ, ݕ, 0) ߲ݖ ࢑ሬሬԦ で表される．ここで，ଙԦ, ଚԦ, ࢑ሬሬԦ はݔ, ݕ, ݖ方向の単位ベクトル である．したがってそれぞれの項は定義から， ߲ܶ(ݔ, ݕ, 0) ߲ݔ ՞ ݅ݑܶ௙(ݑ, ݒ), ߲ܶ(ݔ, ݕ, 0) ߲ݕ ՞ ݅ݒܶ௙(ݑ, ݒ), ߲݅ܶ(ݔ, ݕ, 0) ߲ݖ ՞ ݏܶ௙(ݑ, ݒ) (9.1) となりその振幅は， |ܣ| = ቈ൬߲ܶ(ݔ, ݕ, 0)_߲ݔ ൰ ଶ + ൬߲ܶ(ݔ, ݕ, 0) ߲ݕ ൰ ଶ + ൬߲ܶ(ݔ, ݕ, 0) ߲ݖ ൰ ଶ ቉ ଵ/ଶ (9.2)

で与えられる．  Fig.5はFig.1（c）に対する複素勾配の振幅であるが， 極磁力と同様に，異常源の形に沿った分布となる．この 方法の特徴は，極磁力とは異なり磁化の方向が未知であ っても異常源をイメージング（単極化）できることであ る． 10．地下の起伏分布による重力・磁力分布  Parker （1972）は，地下の起伏分布をテーラー展開し ,高次の項まで取り扱うことにより，より複雑な起伏分 布を持つ構造に対する重力・全磁力異常分布の表現式を 示した．この方法の特徴は，テーラー展開した起伏分布 をさらにフーリエ変換することにより，高速計算に優れ たモデリングの方法を示したことにある．  いまFig.6に示すような深さܦを中心に݄(ݔ, ݕ)で定義さ れる起伏分布を考えると， ݃௙_{(ݑ, ݒ) = 2ߨܩ} ௢݁ݔ݌(െܦݏ) ෍ ݏ௡ିଵ ݊! ஶ ௡ୀଵ ܨ[ߩ(ݔ, ݕ)݄௡_{(ݔ, ݕ)]  (10.1)} となる．ここに， ෍ݏ௡ିଵ_݊! ஶ ௡ୀଵ ܨ[ߩ(ݔ, ݕ)݄௡_{(ݔ, ݕ)] = ܨ[ߩ(ݔ, ݕ)݄(ݔ, ݕ)] +} ݏ 2ܨ[ߩ(ݔ, ݕ)݄ଶ(ݔ, ݕ)] + ݏଶ 6ܨ[ߩ(ݔ, ݕ)݄ଷ(ݔ, ݕ)] ڮ であり，式（4.9）に対比すると物体の形状関数のフー リエ変換である．ここにܨ[ ]はフーリエ変換である． また，ここでは密度ȡも２次元的に変化するものとして いる．したがって式（10.1）のフーリエ逆変換により ݃(ݔ, ݕ)を得る.また基盤の下面を考慮する場合には，やは り平均深度Dからの下面の起伏，すなわちݎ(ݔ, ݕ)を用い て，式（10.1）のフーリエ変換を示すܨ[ ]内の表現を， hĺ(h㸫 r㸧 㸪 h2_ĺ(h2_{㸫 r}2_{㸧 㸪…… (10.2)} とすればよい．また全磁力異常の場合にはߩ(ݔ, ݕ)が磁化強 度݉௢(ݔ, ݕ)となり，2ߨܩ௢をെ(ܮݑ + ܯݒ െ ݅ܰݏ) ή (݈ݑ + ݉ݒ െ ݅݊ݏ)/ݏ で置き換える．ここでݏで割っているのは，式（10.1） を重力ポテンシャルに一度変換し，続いて磁化方向およ び地球磁場方向に微分することを意味している．  以上述べた方法を用いて，Fig.1（a）のモデルについ て計算した重力異常および全磁力異常の結果をFig.7 （a），（b）にそれぞれ示す．Fig.1（b），（c）とほぼ同様 な結果が得られている．なおこれらの結果は式（10.1） のn=4までを考慮しているが，n=3でも大きな差はない． 11．反復法による構造解析法  前節で述べたParker （1972）のモデリングの方法は， その計算の高速性に優れており，これを用いてParker and Huestis （1974）は地磁気の波数領域での反復法に よる逆解析手法を提案した．この方法は，磁化を帯びた 層が地形と平行に一定の厚さを持つと仮定しているこ と，また下方接続の項が含まれていることから，高周波 成分の発散を回避するためのローパス・フィルタが用い られる．本節では式（10.1）のフーリエ逆変換を用いて， 重力異常についての空間領域での反復計算による地下構 造の解析手法について述べる．この方法の特徴は，下方 F i g . 5. A n a l y t i c s i g n a l s o f t h e t o t a l m a g n e t i c

anomalies of Fig.1(c). The symbols H indicate higher than surroundings.

Fig. 6. Schematic model of a subsurface double-layer structure. The function ΰ(ݔ, ݕ) indicates a two-dimensional boundary surface from the depth D.

接続の操作が含まれないことから，ローパス・フィルタ 等を一切施す必要がないことにある．また擬似重力を用 いれば，磁気異常の解析にも適用できる．  いま式（10.1）から得られる重力異常݃௖௔௟_{(ݔ, ݕ)と，観} 測された重力異常（あるいは擬似重力）݃௢௕௦_{(ݔ, ݕ)から次} のような反復計算によって基盤の起伏分布݄(ݔ, ݕ)，ある いは密度分布ߩ(ݔ, ݕ)（あるいは磁化強度分布݉௢(ݔ, ݕ)）を 求める．  すなわち，n回目のモデル(hn,ȡn)に対して， ο݃௡(ݔ, ݕ) = ݃௢௕௦(ݔ, ݕ) െ ݃௖௔௟(ݔ, ݕ) として修正量ǻh，およびǻȡを式（10.1）のDを用いて， ο݄(ݔ, ݕ) = ο݃௡(ݔ, ݕ) 2ߨܩ௢ߩ௡(ݔ, ݕ) , (11.1) あるいは， οߩ(ݔ, ݕ) =ο݃௡(ݔ, ݕ) 2ߨܩ௢ܦ (11.2) で与え，n+１回目の結果を縮小因子ߙ_ΰ, ߙఘ（0.0 ∼ 1.0） を用いて ݄_௡ାଵ(ݔ, ݕ) = ݄௡(ݔ, ݕ) + ߙΰο݄(ݔ, ݕ) , (11.3) あるいは， ߩ_௡ାଵ(ݔ, ݕ) = ߩ௡(ݔ, ݕ) + ߙఘοߩ(ݔ, ݕ) (11.4) とする．これらの修正を繰り返すことにより，基盤の起 伏あるいは密度分布を推定することができる．なお式 （11.2）においてDで割っているが，本来ならば式（11.1） の 類 推 か ら 起 伏 量 分 布݄௡(ݔ, ݕ)が 置 か れ る べ き だ が， ݄௡(ݔ, ݕ)が小さな部分では式（11.2）の結果が発散してし まう．それを避けるためDで置き換えているが，実際に は 縮 小 因 子 と 同 等 な 意 味 を 持 つ． ま たParker and Huestis （1974）のように一定の層厚を仮定する場合は， （10.2）の関係を用いればよい．  Fig.8（a）にFig.1（b）の重力異常から反復法により 求めた地下構造の起伏分布を示す．計算では密度を 1,000 kg/m3_{と固定し，縮小因子を0.5とした．真のモデ} ルよりやや裾が広がった円錐形に近い形となっている． またFig.8（b）は，地下構造を固定してFig.2（a）の擬 似重力より反復法により求めた磁化強度分布である．真 のモデルでは円錐台が一様に１A/mであるが，求めら れた結果は円錐台脚部から頂上に向かって磁化強度が大 Fig. 7. (a) Gravity and (b) total magnetic anomalies calculated by 2D FFT for the model in Fig.1(a). The solid lines indicate positive values and the dashed lines, negative values. The symbols H and L indicate higher and lower than surroundings, respectively.

きくなる．これは式（11.3），（11.4）において反復過程 での修正量をその直下の構造にのみ加えていくために， 地下構造の高周波成分の特徴がやや損なわれるものと考 えられる．これは下方接続による逆解析で用いるローパ ス・フィルタの影響と同様な問題ではあるが，層厚一定 の仮定が必要ないこと，重力異常，磁気異常の問題を共 通のアルゴリズムで解けるという特徴を持っている． 12．おわりに  高速フーリエ変換を利用する前提にたって，重力異 常，磁気異常の解析方法のうち，そのいくつかについて 有限複素フーリエ級数を用いて表現した．こうすること によって，空間領域でとり扱われる微分，積分と波数領 域で現われる各種の係数との関係を明確にすることがで きた．  全磁力異常を擬似重力に変換することによって，重力 による地下構造の解析手法がそのまま利用できる．また 植田（1990）やUeda et al. （2002）は重力異常，磁気異 常の異常源が共通であると仮定して，重力異常と擬似重 力から磁化の方向と密度/磁気モーメント比を求める方 法を示した．

 Lourenco and Morrison （1973）による全磁力異常か ら磁場の３成分を求める方法は，フィルタリング技術と しては極めて興味深く，例えば大久保（1981）はこの方 法を空中磁気データに適用して，残留磁化情報の抽出に ついて議論している．また複素勾配法は，海底ケーブル や危険物探査など磁場の方向が明確でない場合に，磁気 探査による異常源の位置特定手法として有効である．  一方，今回示した反復法による解析手法は，結果が元 の重力異常分布や擬似重力分布とやや相似形になるとい う問題はあるが，水深に比べ十分に広いデータを扱う場 合には計算の高速性に優れており，重力による地下構造 解析，あるいは面的な磁化強度分布を把握する手法とし て極めて有効である．しかし，単体の火山あるいは海山 等を扱う場合には，深さ方向の変化も考慮しなければな らないことから，例えば植田ほか（2001）や久保田ほか （2003），Kubota and Uchiyama （2005）等による３次元

インバージョンによる解析手法が有効である．  なおここに掲載したプログラムを利用した場合には， 本論文によったことを書き添えていただくようお願いす る． 13．謝辞  本研究は川崎地質㈱および海洋情報部の方々と海洋の 地磁気・重力や航空磁気測量のデータ解析を通して議論 を重ねてきたものであり，ここに深く謝意を表する．ま た本論文を執筆するにあたり査読者からは，有益な助言 を得たことに深く謝意を表する．

Fig. 8. (a) Subsurface structure derived from the gravity anomalies of Fig.1(b) by using our new iterative method. (b) Magnetizations derived from the pseudo-gravity anomalies of Fig.2(b) by using our new iterative method. The symbols H indicate higher than surroundings.

引 用 文 献

Baranov, V.: A new method for interpretation of  a e r o m a g n e t i c m a p s : p s e u d o - g r a v i m e t r i c anomalies, Geophysics, 22, p.359-383, （1957） Baranov, V. and H. Naudy: Numerical calculation of the

formula of reduction to the magnetic pole, Geophysics, 29, p.67-79, （1964）

Bhattacharyya, B. K.: Magnetic anomalies due to prism-shaped bodies with arbitrary polarization, Geophysics, 29, p.517-531, （1964）

Bhattacharyya, B. K.: Two-dimensional harmonic analysis as a tool for magnetic interpretation, Geophysics, 30, p.829-857, （1965）

Bhattacharyya, B. K. and M. E. Navolio: A fast Fourier transform method for rapid computation of g r a v i t y a n d m a g n e t i c a n o m a l i e s d u e t o arbitrary bodies, Geophys. Pros., 24, p.633-649, （1976）

Blakely, R. J: Potential theory in gravity and magnetic applications, Cambridge Univ. Press, （1996） 物理探査学会：新版物理探査用語辞典，愛智出版,（2005） 萩原幸男:地球重力論，共立全書，（1979） 長谷川博:擬似重力異常の新計算法，物理探鉱，20, p.22-31, （1967） 加藤元彦:２次元フィルターの理論と重力・磁力分布の 解析，ラテイス，（1987） 駒澤正夫：有限フーリエ級数を用いた重力の３次元解析 について,物理探鉱，33, p.108-115, （1980） 久保田隆二・植田義夫・小野寺健英・瀬尾徳常・加藤剛 ：海域火山の３次元地磁気・重力インバージョ ン，物理探査学会第109回学術講演論文集，p72-75，（2003）

Kubota, R. and A. Uchiyama: Three-dimensional magnetization vector inversion of a seamount, Earth Planets Space, 57, p.691-699, （2005） Lourenco, J. S. and H. F. Morrison: Vector magnetic

anomalies derived from measurements of a single component of the field, Geophysics, 38, p.356-369, （1973）

大久保泰邦:空中磁気図より求まる残留磁気の方向，物 理探査, 34, p.399-413, （1981）

大崎順彦:新・地震動のスペクトル解析入門，鹿島出版 会，（1994）

Parker, R. L.: The rapid calculation of potential anomalies, Geophys. J. Int., 31, p.447-455, （1973）

Parker, R. L. and S. P. Huestis: The inversion of m a g n e t i c a n o m a l i e s i n t h e p r e s e n c e o f topography, J. Geophys. Res., 79, p.1587-1593, （1974）

Rosenbach, O.: A contribution to the computation of the second derivative from gravity data, Geophysics, 18, p.894-909, （1953）

Saito, M. and H. Takeuchi: Gravity anomalies and the c o r r e s p o n d i n g c r u s t a l s t r u c t u r e . P a r t 1） numerical tables useful in three-dimensional gravity interpretations, Bull. Earthq. Res. Inst., 42, p.39-92, （1964）

Schouten, H. and K. McCamy: Filtering marine magnetic anomalies, J. Geophys. Res., 77, p. 7089-7099, （1972）

Talwani, M. and M. Ewing: Rapid computation of gravitational attraction of three-dimensional bodies of arbitrary shape, Geophysics, 25, p.203-225, （1960）

Tsuboi, C. and Y. Tomoda: The relation between the FOURIER series method and sinx/x method for gravity interpretations, Jour. Phys. Earth, 6, p.1-5, （1958） 植田義夫：応答関数法による地磁気・重力異常の総合的 解析法の開発とその応用結果，海洋調査技術，２ , p.49-60，（1990） 植田義夫・中川久穂・小野寺健英・鈴木晃・熊川浩一・ 久保田隆二: 2000年噴火前の三宅島の３次元磁気 構造−地磁気トモグラフィーの試み−，水路部研 究報告，37，p.19-36，（2001）

Ueda, Y., R. Kubota, and J. Segawa: Magneto-gravity response function and its application to the Daito Ridge, Geophysics, 67, p.110-116, （2002）

付録A  ２次元高速フーリエ変換プログラム FAST2 を以 下に示す． 目的  与えられた２次元の複素数型データの高速フーリエ変 換またはフーリエ逆変換を，単精度で計算する． 使用法 （１）接続方法 引数 型* 副プログラムを呼ぶ時の内容 副プログラムから_{戻ったときの内容} NX I 複素数型データおよび_{複素数型変換値のx方向の数} 不変 NY I 複素数型データおよび_{複素数型変換値のy方向の数} 不変 X C 2次元配列の複素数型データ フーリエ変換値のNX×NY倍または 逆フーリエ変換値 NDX I 主プログラムにおけるXの_{x方向の長さ} 不変 NDY I 主プログラムにおけるXの_{y方向の長さ} 不変 IND I フーリエ変換のとき IND=−１ フーリエ逆変換のとき IND=＋１ 不変  * 型 I は整数型，C は複素数型． （２）注意事項 ⅰ）NX，NYは２の累乗でなければならない． ⅱ） フーリエ変換IND=−１のとき，変換値はNX×NY 倍されている． （３）必要なサブルーチンプログラム．大崎（1994）に 記されている１次元データに対する高速フーリエ変換プ ログラム FAST ．

付録B  ２次元高速フーリエ変換による極磁力変換プログラム RTP を以下に示す． 目的  与えられた２次元空間(ݔ, ݕ)の全磁力異常の極磁力変 換を，単精度で計算する． 使用法 （１）接続方法 引数 型* 副プログラムを呼ぶ時の内容 副プログラムから_{戻ったときの内容} NX I 実数型磁気異常データおよび_{極磁力変換値のx （東）方向の数} 不変 NY I 実数型磁気異常データおよび_{極磁力変換値のy （北）方向の数} 不変 A R ２次元配列の_{実数型全磁力異常データ} 極磁力変換値 NDX I 主プログラムにおける_{Aのx方向の長さ} 不変 NDY I 主プログラムにおける_{Aのy方向の長さ} 不変 DX R x方向のサンプリング間隔 不変 DY R y方向のサンプリング間隔 不変 DECG R 地球磁場の偏角 不変 INCG R 地球磁場の伏角 不変 DECM R 磁化方向の偏角 不変 INCM R 磁化方向の伏角 不変  * 型 I は整数型，R は実数型． （２）注意事項 ⅰ）NX，NYは2の累乗でなければならない． ⅱ） 配列Aの単位はnT，DXおよびDYの単位はm，偏 角および伏角の単位は度．A（1,1）は南西端の位置 を示す． （３）必要なサブルーチンプログラム．付録Aの２次元 高速フーリエ変換プログラム“FAST2”．