数値座標変換による比抵抗法 3 次元データ解析に関する研究

(1)

数値座標変換による

比抵抗法 3 次元データ解析に関する研究

On the 3-D Data Interpretation in Resistivity Method Using Numerical Coordinate Transformation

2003 年 3 月

汪振洋

(2)

第 1 章序論

1.1 研究の背景

電気探査法は、地下情報を取得する技術として、産業革命の勃興とともに発展してきた。

とくに、19世紀後半から20世紀前半には、欧米諸国を中心として近代化が加速し、物質文明の進歩を支える基盤ともいえる様々な鉱物資源の需要が大幅に増加した。地下資源探査の有力な手段である電気探査技術はこの時期に著しく発展し、近代電気探査の基礎はほぼ完成された。

しかし、これら鉱物資源探査を主な目的として開発されてきた電気探査技術は、20 世紀後半からその応用分野を急速に広げ、地球科学研究を初め、土木建設、防災、遺跡調査、

環境保全などの多くの分野で適用されるようになった。このような適用の拡大に伴い、より複雑な地下構造に対して、より精密な地下情報の取得が要求されるようになった。具体的には、起伏地形への対応や、地下浅部から深部までの詳細な 3 次元可視化技術などの実用化が急務となった。

電気探査法には様々な手法があるが、大別すれば、自然的なソースを用いる受動的な探査手法と人工的なソースを用いる能動的な探査手法との 2 種類に分けることができる。能動的な探査手法である比抵抗法は、そのソースの幾何的、物理的な特徴を完全に把握できるため、S/N比が高く、最も多く適用されている。

比抵抗法は人工電位法の一種であり、大地に接地した電極を通して電流を地下へ注入し、

他の接地電極から電位あるいは電位差を測定し、その電場の変化から地下比抵抗分布を調べる物理探査手法である。比抵抗法は、その歴史が古く、金属鉱物資源探査、地下水調査などを目的として開発され、多くの実績を残してきている。とくに、その測定法ならびに解析理論の発展は、近代電気探査の全般に亘って大きく貢献してきている。

比抵抗法の発展には、近年のエレクトロニクス技術革新の恩恵を多大に受けており、データ取得技術の面ではディジタル化・マルチチャンネル化など、解析のハードウェアの面ではコンピュータの小型・高速・大容量化などにより、飛躍的な進歩を遂げた。しかしながら、地下の 2 次元断面探査については、ソフトウェアと解析理論の面でほぼ確立されてきているのに対し、3次元探査については、最近の 10年間に多くの研究が行われてきてはいるものの、とくに数値解析における解析の不安定性、解像度と効率の低下という難題が残されたままであり、いまだに実用段階には至っていない。

(6)

1.2 研究の目的

比抵抗法の3次元逆解析法の研究は、最小2乗法を中心とした1次元・2次元解析法の拡張という形で進められてきた。しかし、現在、比抵抗法 3 次元データ解析においては、解析結果の解像度が低いという問題に直面している。この問題は、数値モデル化における離散化誤差だけに起因する問題ではなく、むしろ、地下パラメータに対する応答の感度特性という、比抵抗法それ自体の特性・能力に大きく起因するため、コンピュータの計算能力の向上だけでは解決されない問題であることが認識されてきている。

従来の逆解析法による比抵抗法3次元解析では、次の問題点が指摘されている。

① 反復計算のたびに感度行列を計算すると非常にコストの高い計算となる。

② 感度行列の近似計算は平坦な地形でなければ成立しないため、起伏地形の解析を行うことはできない。

③ データに含まれる情報量と逆解析で推定しようとするモデル情報とが互いに矛盾関係にあるため、解析は非常に不安定となる。また、安定化操作は解像度の低下を招く。

これらの問題の解決に向けて、最近では、主に次の 3 つのテーマについて研究が行われている。すなわち、解析の効率と精度に深く関わる地下空間の離散化方法についての研究、

解析の安定性・効率性と深く関わる感度マトリックスの計算方法の研究、および逆問題における解の存在性・一意性と深く関わるデータの取得方法についての研究である。

とくに、順問題と逆問題の解の唯一性の検討は、安定的な解析手法を開発するための重要な指針を導くものであり、それは境界値問題の検討と密接な関係がある。すなわち、

① 全データ空間がわからなくても、調査したい領域をなるべく囲むような境界で測点を配置すれば、逆解析の安定性と唯一性が向上する。

② それが達成できない場合は、既知のすべての情報を用いて、データの欠けた境界部分にできるだけ良い近似条件を与える。

③ 既知データではコントロールできないモデルの部分空間に対して無駄な推定計算をしない。

である。

本研究では、信頼性、効率性と解像度を追求した 3 次元数値解析法を開発することを目的として、次の事項を検討課題として、比抵抗法3次元データ解析法の開発研究を行った。

(7)

① 解析の精度と効率を向上させるための空間離散化の方法

② 任意地形、任意地下比抵抗構造に対応できる順解析法

③ 解析の安定性と唯一性を向上させるためのデータ取得法

④ 信頼性、効率性と解像度を追求した逆解析法

なお、本研究では､2 極法(Pole-Pole 配置)の測定に対応した電位のみについて検討を進めた。その理由は､2 極法測定に対応した電位が求められていれば、相加定理などを用いて、

任意の電極配置による電位(電位差)あるいは見掛比抵抗を容易に求められるからである。

1.3 本研究の概要

本研究では、まず従来の 3 次元数値解析法の問題点を分析し、最適な空間離散化の方法と最適な計算空間理論を検討した上で、数値座標変換による新たな 3 次元数値解析法を導いた。

各章の概要は以下の通りである。

第1章は序論であり、研究の背景、目的と各章の内容について述べた。

第2章では、比抵抗法およびその数値解析法の背景と基礎概念を述べ、3次元データ解析における問題点を明確にした。まず、比抵抗法発展の歴史と最近の研究動向から本研究の位置づけと範囲を明示し、比抵抗法における支配方程式およびその境界条件の導出と基礎概念を紹介した。次に、順解析と逆解析を符号化により定義し、逆解析の不安定性の原因は、測定ノイズを含むデータ部分空間からモデル全空間を推定しようという、比抵抗法における逆解析の本質にあることを指摘した。また、従来の逆解析法について検討し、その解析法での解像度と安定性との矛盾関係は、コンピュータの能力向上だけでは解決されないことを明らかにした。さらに、順解析と逆解析の唯一定理を紹介し、解析の安定性と解像度を向上させるための重要な課題は、空間のより適切な離散化と、その自動化にあることを明らかにした。

第3 章では、比抵抗法の3 次元解析における問題点を具体的に絞り、座標変換による問題解決の理論的な基礎を導出した。まず、2次元解析から3次元解析へ移行する場合に直面する 2 つの問題、すなわち、未知変数の増大による計算量増大の問題と、パラメータに対する応答感度の距離減衰が速くなることによる解析の不安定性の問題を指摘した。この 2 つの問題は、それぞれ順解析と逆解析に相関しているため、それらを解決するためには空

(8)

間の離散化を最適化する必要があることを明らかにした。通常、デカルト座標で表される物理的な空間（以下、物理空間）の離散化を最適化することは、計算を実行する空間（以下、計算空間）を最適化することに等しいことを導いた。物理空間における空間離散化の最適化の問題は、順解析と逆解析とで判断基準が異なるが、最適な計算空間を導入することによって、それらを統一できることを明らかにした。この最適な計算空間とは、物理空間での支配方程式と任意地形を、単純な方程式と平坦地形で表すことのできる座標空間である。それは、計算空間の基本計量テンソルの選択と適切な境界条件の適用により達成され、導出された双曲型変換方程式と楕円型変換方程式およびそれぞれの境界条件により、

理論的に両空間を一対一の対応関係で結びつけることができた。

この、物理空間および計算空間それぞれのグリッド系の座標対応関係を求める過程は、

最適な空間離散化であると同時に、比抵抗法の数値解析そのものである。とくに順問題の解は、最適な計算空間においては単純な解析式で求められることから、物理空間と計算空間との数値座標変換を行うことによって、比抵抗法の順解析が行われることを明らかにした。

以上により、最適な空間離散化と密接な関係を有する最適な計算空間の概念を用いて、3 次元的な任意地形と任意地下比抵抗構造に対応する、効率的な数値解析法の理論的基礎を確立した。

第 4 章では、座標変換問題を数値的に解くためのアルゴリズムを構築した。この数値座標変換は、物理空間および計算空間それぞれのグリッド系節点座標の対応関係を求めることに相当するが、物理空間の地表グリッドは、送信点、受信点に対応し固定されるため、

物理空間のグリッド系を再構築する段階で支障が生じる。そこで、両方向への座標変換方程式を一つの座標空間で解くことができるような実質的な計算空間（以下、媒介空間）を導入し、この問題を解決した。そのために、双曲型変換方程式およびその境界条件について媒介空間座標を引数とする方程式へ書き直し、方程式の因子化と分散項の追加を行った。

また、楕円型変換方程式およびその境界条件について媒介空間座標を引数とする方程式へ書き直し、有限体積法により、Picard反復法の適用可能な差分方程式を導出した。

以上により、数値座標変換システムの構成と変換アルゴリズムを明らかにした。また、

開発した数値座標変換システムを用いて起伏地形モデルに対する変換例を示した。

第 5 章では、比抵抗法の順問題の解は、計算空間での単純な境界値問題の解に相当することを理論的に証明し、数値座標変換による順解析のアルゴリズムを開発した。このアルゴリズムの精度を検討するために、半無限均質大地モデルと不均質半球モデルの 2 例について、計算結果と理論解とを比較し、両モデルとも、極めて高い精度で解が得られることを示した。とくに、特異点となる電流ソース付近でも精度よく求めることができ、この点で、開発した数値座標変換法は、従来の数値計算法よりはるかに優れていることが明らか

(9)

となった。また、起伏地形問題および異常体と媒質との比抵抗コントラストの問題に対して、アルゴリズムの安定性を確認した。さらに、従来、膨大な計算量を必要とした感度解析が、この数値座標変換を用いると、解析式で極めて容易に計算できることを明らかにした。

以上により、この数値座標変換法は、3次元的な任意地形と任意地下比抵抗構造に関する比抵抗法の数値計算に対応し、有限規模のグリッド系で、地下空間の離散化を自動的に最適化することができる。また、電流ソースの位置が変化しても、新たに微分方程式を解く必要がないため、効率と精度の両面で有利であることを明らかにした。

第 6 章では、比抵抗法における逆問題の解は、計算空間での境界値問題の解に相当することを理論的に証明し、効率的なデータ取得法ならびに新たな逆解析法を開発した。ここでは、数値座標変換の利点を十分に活かすために、感度行列を陽的に計算しない新たな逆解析法を開発した。まず、比抵抗法の測定（電位）データと計算空間地表グリッド系との関係を明らかにし、測定データを計算空間グリッド系節点座標に変換する方法を導出した。

また、与えられた地表電極系から、地下情報を最大限に引き出すための送受信配列およびデータ取得の方法について考察した。さらに、本逆解析手法は、計算空間地表グリッド系を適合させる反復計算法であることを理論的に証明した上で、数値座標変換による比抵抗法 3 次元逆解析アルゴリズムを導いた。このアルゴリズムを用いて逆解析計算を行い、計算精度の検討、解析の安定性、解像度ならびに効率性を評価し、さらに、比抵抗 3 次元探査法に関する今後の研究課題を提言した。

以上により、数値座標変換による 3 次元数値解析法は、解析に必要となる測定データセットを推定できる点、解析過程の中で自動的に地下の離散化が達成される点で、従来の解析法とは一線を画す優れた数値解析法であることを明らかにした。

第7 章は結論であり、本研究の成果をまとめた。また、比抵抗3次元探査法の実用化に向けた今後の研究展望を述べ、論文を総括した。

(10)

第 2 章比抵抗法の基礎理論と解析

2.1 はじめに

電気探査比抵抗法は、人工電位法の一種であり、大地に接地した電極を通して電流を地下へ注入し、他の接地電極から電位あるいは電位差を測定し、その電場の変化から地下比抵抗分布を調べる物理探査手法である。比抵抗法は、その歴史が古く、金属鉱物資源探査、

地下水調査などを目的として開発され、多くの実績を残してきている。とくに、その測定法ならびに解析理論の発展は、近代電気探査の全般に亘って大きく貢献してきている。

比抵抗法の3次元逆解析法の研究は、最小2乗法を中心とした1次元・2次元解析法の拡張という形で進められてきた。しかし、現在、比抵抗法 3 次元データ解析においては、解析結果の解像度が低いという問題に直面しており、これは、コンピュータの計算能力の向上だけでは解決されない問題であることが認識されてきている。

比抵抗法は、近年のエレクトロニクス技術革新の恩恵を受けて飛躍的な進歩を遂げてきており、地下の 2 次元断面探査については、ソフトウェアと解析理論の面でほぼ確立されてきている。それに対して3次元探査については、最近の10年間に多くの研究開発が行われてきてはいるものの、とくに数値解析における解析の不安定性、解像度と効率の低下という難題が残されたままであり、いまだに実用段階には至っていない。

本章では、電気探査比抵抗法およびその数値解析法の発展の経緯・現状と基礎概念を整理し、3次元データ解析における問題点を明確にする。

2.2 比抵抗解析法の発展経緯

電気探査法を用いて地下を調べようとする努力は、1830年代の Fox R. W.によるPenzance

Mine of Cornwallでの自然電位調査までさかのぼるが、人工ソースを用いて地下の比抵抗分

布と鉱床の関係を調べる試みは、1883~1891年のBrown Fred H.の研究に始まる。1912年、

Schlumberger (1920)は、初めて電気探査比抵抗法を実際の調査に適用し、同時に、比抵抗法の測定中に地下のIP効果を発見している。

見掛比抵抗の概念は、Wenner(1912) とSchlumberger(1920)のそれぞれによって導入された。

Wenner は等間隔の電流、電位電極配置を使用し、Schlumbergerは電流電極の間隔に比べ極

めて小さな間隔の電位電極を電流電極間の中央に置く配置を使用した。現在、比抵抗法で用いられる電極配置の種類は非常に多く、2 極法配置（Pole-Pole 配置）、3 極法配置

(11)

(Pole-Dipole)、4極法配置がある。4極法配置には、上述のWennerとSchlumbergerの配置に加えて、Dipole-Dipole 配置、Eltran配置、Staggered配置などがある。

データ解析法は、Stefanesco et. al (1930)の水平層状大地(1次元)モデルに対する理論解に関する研究から発展した。彼らは実測データと理論曲線との適合(Fitting)により、地下比抵抗構造を推定する方法を開発した。この解析法が開発されたことにより、それまで定性的解釈にとどまっていた比抵抗法は、定量的な探査法へと大きく発展した。現在の解析法、すなわち、順・逆解析の考え方の基礎はその当時に形成されたと言える。この解析法は 1960 年代まで標準的な解析法として用いられた。

2次元数値解析では、Tagg (1930)が垂直断層モデルの理論解を、Frank and Von Mises (1935) が点電流ソースによる埋没水平円柱モデルの理論解を導いたのが始まりである。その後、

各種の2・3次元モデルに対する解析解が導かれたが、それらは、実測データに対する定性的な解釈には有効となるものの、定量的な解析にはあまり役に立たなかった。この状況は、

1960 年代末、偏微分方程式数値解法の電気探査法順問題への適用により一変した。Jepsen

(1969)は、差分法を比抵抗法ならびIP法の順問題に適用した。Coggon (1971)は、有限要素

法を用いた電気・電磁探査法の順問題の数値解析法を明らかにした。彼は、電磁エネルギーの保存則から導出した汎関数を用いて理論を展開し、3次元ソース、2次元比抵抗構造モデルの場合について検討した。ここでは、支配方程式の比抵抗走行方向の空間座標に対してフーリエ-変換・逆変換を行うことにより、3次元電場の問題を純 2次元問題に変換して解く方法が明らかにされた。この数値モデリング法は、現在の電気・電磁探査法 2 次元モデリングの標準的な方法となっている。

3次元数値モデリングについては、Stefanesco (1950)によって近似計算法(Alpha-Center 法) が提案され、さらに、Hohmann (1975)は、積分方程式法を比抵抗法、IP 法ならびに電磁探査法の数値計算に適用した。その後、Pridmore (1978)は有限要素法を、Dey and Morrison (1979) は差分法をそれぞれ比抵抗法3次元モデリングに適用した。

逆解析の電気探査法への適用は、Slichter (1933)が層状大地モデルについて試みたのが始まりである。しかし、一般逆解析(Generalizes Inversion)は、Ghosh (1971)がLinear Filter 理論を水平多層構造モデル計算に適用したことにより、順解析の効率と精度が大幅に向上した後のことである。Inman et. al (1973) は、最小2乗法による逆解析を層状構造モデルデータ解析(1次元解析)に適用させた。さらに, Inman (1975)は、制約付き最小2乗法反復法(Ridge

Regression Inversion)を1次元解析に適用し、ノイズを含むデータであっても安定に解が収束

することを示し、比抵抗法逆解析の実用化に大きく貢献した。

比抵抗法の2次元逆解析は、Coggonの順解析アルゴリズムとInmanの逆解析アルゴリズムとを結合させることによって可能となった。しかし、当初はコンピュータの計算能力不足で多くの制約を受けたが、その後のコンピュータ分野の飛躍的な発展を背景にして、佐々木(1981)は、Coggonと Inman のアルゴリズムを完全に結合させた。それにより、比抵抗法データ解析法の主要部分はほぼ完成されたといえる。以後の研究は、アルゴリズムの安定

(12)

化、効率化および解析精度の向上に集中した。1980年代後半から、比抵抗法 2次元逆解析法の適用が急速に広がり、実用性が向上した。その成功は比抵抗法だけではなく、電気・

電磁探査の逆解析研究の全般に亘って大きな刺激を与えた。それと同時に、比抵抗法 2 次元探査の限界が次第に明らかとなり、3次元探査とその解析の実用化への要求が高まった。

比抵抗法3次元データ解析では、Park and Van (1991)は、順解析にはDey and Morrisonの差分法を改良して用い、逆解析には平滑化制約付き最小 2 乗法反復法を用いて、実際の調査データの解析を行い、比抵抗法3次元逆解析適用の先鞭をつけた。佐々木(1992)は、3次元有限要素法を順解析に用い、また計算効率を目的として、均質大地の感度行列を近似的に利用した3次元逆解析近似法を提案した。さらにSasaki (1994)は、感度行列を精確に求める 3 次元逆解析アルゴリズムを示し、近似逆解析法との比較を行った。Li and Oldenburg(1994)は、差分法の順解析に近似逆解析法を加えるアルゴリズムを提案した。このように比抵抗法3次元データ解析法は、計算コストを考慮した、2次元解析法の拡張として短期間に開発された。

しかし、現在、比抵抗法 3 次元解析においては、解析結果の解像度が低いという問題に直面している。この問題は、順解析における離散化誤差だけに起因する問題ではなく、

むしろ、地下パラメータに対する応答の感度特性という、比抵抗法それ自体の特性・能力に大きく起因するため、コンピュータの計算能力の向上だけでは解決されない問題であることが認識されてきている。

2.3 比抵抗法の基礎理論

2.3.1 支配方程式と境界条件：

比抵抗法は、地下物質の電気特性を利用する物理探査法の1つである。その基礎理論は、

次のMaxwellの方程式から導出される。

t

∇× = − ∂

∂

E B (2.1a)

s t

∇× = + + ∂

∂

H J J D (2.1b)

ρ

e

∇ ⋅ = D

(2.1c)

∇⋅ = B 0

(2.1d)

ここで、

E：電場強度

(13)

H：磁場強度 B：磁束密度 D：電束密度

J

：伝導電流密度

ρ

e^{：自由電荷密度}

σ

^：導電率

ε

^：誘電率

µ

^：透磁率

真空中であれば

[ ]

σ = 0 S m

[ ]

ε ε =

₀

= π 1

₂

×

⁷

4 10

c F m

[ ]

µ µ =

0

= π ×

⁻

4 10 7 H m

c

^{：真空中の光の速度}

である。

(2.1b)式の t

∂

D項は変位電流と呼ばれる。

J

_sの項はソース電流密度であり、伝導電流密度

σ

J = E

(2.2)

には依存しない電流密度の部分を示す。式(2.2)は微分形式のオームの法則である。伝導電流密度

J

と電場Eは、ともにベクトルであるから、一般的に、導電率

σ

^{はテンソル関数であ}

ることがわかる。本研究では、導電率

σ

は空間任意一点で等方性であると仮定する。この仮定を用いると、導電率

σ

は空間座標のスカラー連続関数に相当する。

ここで導出する比抵抗法の基礎理論は、次の仮定のもとで成立する。

① 低周波数近似が成立する。

② 導電率

σ

は電場と独立である。

③ 電流ソースは、意図的に大地に注入する人工ソース以外には存在しない。

仮定①により、式(2.1b)の変位電流を無視でき、さらに、式(2.1a)の磁束密度の時間微分も無視できる。

仮定②により、導電率

σ

は、印加外部電場の時間(周波数)と強度変化によらず、常に定数である。

仮定③により、電流ソースの幾何的、物理的な特徴はすべて把握できることになる。

(14)

これらの仮定を満足させるために、実際の測定では、とくにデータ取得の際、細心の注意を払わなければならない。すなわち、仮定①を満足させるためには、測定中の電磁カップリングや容量カップリングの影響を取り除く必要がある。仮定②を満足させるためには、

地下物質のIP効果が測定データに含まれていてはならない。仮定③を満足させるためには、

電極の接触電位や地下の自然分極効果、地電流などの影響を取り除いておく必要がある。

これらの仮定のもとで、式(2.1a)、(2.1b)はそれぞれ次の2つの方程式に書き換えられる。

∇× = E 0

(2.3)

σ

s

∇× H = E + J

(2.4)

式(2.3)から、スカラーポテンシャル関数の電位

U

を次のように定義することができる。

=−∇ U

E

(2.5)

ここで、式(2.4)の両辺で発散を取ると、次式が得られる。

( σ )

_s

∇⋅ E =−∇⋅ J

地下空間の

r

₀の点電流ソースから強度Iの電流を注入した場合、上式は次のようになる。

( σ ) I δ (

0

)

∇⋅ E = r − r

(2.6)

式(2.6)は微分形式の電流保存則である。

Γ

_∆_Vを外境界面とする任意体積

∆ V

^で式(2.6)^に対

して体積積分を行い、Gaussの積分定理を適用すると、次の積分形式の電流保存則が得られる。

0 0

0 ( )

( )

V

r V

d I r V

σ

Γ∆

 ∉∆

⋅ =   ∈∆

∫

Ñ ^E ^Γ

^(2.7)

一方、閉曲線

L

_∆_Aを外境界曲線とする任意曲面

∆

Aで式(2.3)に対して面積積分を行い、

Ampereの環路積分定理を適用すると、次式が得られる。

0

LA

d

∆

⋅ =

∫

Ñ

^E ^L ^(2.8)

式(2.5)を式(2.6)に代入すると、比抵抗法の支配方程式が得られる。

( σ U ) I δ (

0

)

∇⋅ ∇ = − r − r

(2.9)

(15)

いま、任意境界面両側の導電率をそれぞれ

σ

₁^と

σ

2^{、電位をそれぞれ}

U

1と

U

₂、境界面の法線方向微分と接線方向微分をそれぞれ

n

∂

^と t

∂

で表す。式(2.5)により、その境界面上の電流ソースが存在しないところでは電場が有限でなければならないため、電位の連続性が要求される。すなわち、

1 2

U = U

^(2.10)

また、式(2.7)により、電流ソースが存在しないところでは、電位法線微分について次式を満足することが要求される。

1 2

U U

n n

σ

∂ =

σ

∂

∂ ∂

^(2.11)

電位接線微分の連続性は式(2.8)により要求され、次式で表される。

1 2

U U

t t

∂ = ∂

∂ ∂

^(2.12)

式(2.9)は楕円型の偏微分方程式である。偏微分方程式理論の唯一定理から、支配方程式の式(2.9)が領域内で、方程式(2.10)~(2.12)のいずれかあるいはその組み合わせが領域の閉合境界で与えてあれば、領域内任意一点の電位関数は唯一に決定される。

2.3.2 相加定理と相反定理について

いま、地下空間任意

N

点r₀ⁱ (i

=

1,2,...,N)^で同時に

N

個の電流ソース

I

_i

( i = 1,2,..., N )

が存在する場合の総合電位

U

を考える。これらの電流ソースがそれぞれ単独に存在する場合の電位を

U

ⁱ

( i = 1,2,..., N )

^{とすれば、総合電位}

U

は次式により求められる。

0

1 1

( , )

N N

i i i

i i

U U I G

= =

= ∑ = ∑

^{r r} ^(2.13)

ただし、

G

はGreen関数である。

この結論は、式(2.9)~式(2.12)の微分演算子の線形性と境界値問題の唯一性から明白であり、

相加定理と呼ばれる。

(16)

一般的には、次の積分形式の相加定理が成立する（Morse,1953）。

'

' ' '

( ) ( ) ( , )

U I G d

Ω

= ∫ Ω

r r r r

(2.14)

このとき、Green関数

G

は、次の微分方程式と式(2.10)~式(2.12)の境界条件を満足する。

( σ G ) δ (

'

)

∇⋅ ∇ = − r − r

(2.15)

ただし、

r

^'とrはそれぞれ点電流ソースの位置と電位測定点の位置を表す。

Green関数

G ( , ) r r

^' の性質として、任意の比抵抗分布と境界形状で、点電流ソース位置

r

^'と

電位測定位置rとを交換しても Green 関数

G

の値は変わらないという対称性があり、これは次式で表される。

' '

( , ) ( , ) G r r = G r r

このGreen関数の性質により、比抵抗法では相反定理が成立する。

2.3.3 電極配置と見掛比抵抗について

導電率

σ

₀を有する半無限均質大地の地表

r

_Aに点電流ソースを置いた場合、任意一点

r

_M

での電位は次式により求められる。

0

( , ) 1

M A

2

M A

U IG I

∆ = = πσ

r r −

r r

^(2.16)

式(2.16)で求められる電位は、無限遠点の電位を0と与えた場合であり、比抵抗法の2極法配置(Fig. 2.1)の測定電位差に相当する。式(2.16)と相反定理、相加定理を用いれば、任意電極配置による測定電位値が求められる。例えば、任意2点

r

_Mと

r

_N との電位差は次式で表される。

0

1 1

2 _M _A _N _A

U I

πσ

 

∆ =   

r

−

r

−

r

−

r

  

^(2.17)

この式(2.17)で与えられる電位差は3極法配置の測定値に相当する。

このことから、4極法配置の測定値も容易に求められ、次式で表される。

(17)

0

1 1 1 1

2 _M _A _N _A _M _B _N _B

U I

πσ

 

∆ =   

r

−

r

−

r

−

r

−

r

−

r

+

r

−

r

  

^(2.18)

見掛比抵抗ρ_aは電極配置により異なるが、次式により統一して定義することができる。

a

K U

ρ

= ∆

I ^(2.19)

ただし、Kは電極係数である。

4極法配置の電極係数は次式で表される。

1 1 1 1

2

M A N A M B N B

K π

 

=   

r

−

r

−

r

−

r

−

r

−

r

+

r

−

r

  

^(2.20)

明らかに、均質大地の場合は

ρ

_a

= 1 σ

₀が成立し、電極配置の影響は完全に取り除かれる。

しかしながら、一般的な3次元大地の場合には、式(2.19)で求められた見掛比抵抗には地下比抵抗分布、地形および電極配置などの影響がすべて含まれている。

本研究では､2 極法(Pole-Pole 配置)の測定に対応した電位のみについて検討を進める。その理由は､上述のように、2極法測定に対応した電位を求めることは方程式(2.9)〜(2.12)の順解析と直接に対応し、その電位が求められていれば、相加定理などを用いて、任意の電極配置による電位(電位差)あるいは見掛比抵抗は容易に求められるからである。

∞ ∞

r

A

r

^M

M A

I

O

σ

0

M

−

A

r r

∆ U

Fig.2.1 Illustration of pole-pole measurement configuration in resistivity method

2.4 順解析と逆解析

(18)

2.4.1 基礎概念

順問題(Forward Problem)、逆問題(Inverse Problem)は、どちらもモデルと応答との関係を追究する問題である。モデルが既知であり、それによって引き起こされる応答を求める問題が順問題である。それに対し、応答が既知であり、それを引き起こす原因となるモデルを推定する問題が逆問題である。順問題、逆問題を解くことは、それぞれ順解析、逆解析と呼ばれる。比抵抗法では地下比抵抗の空間分布がモデルであり、大地に電流ソースを与えた場合の電位分布が応答である。

上述の概念は、次のように符号化して表されることが多い(例えば、Li and Oldenburg(1994))。

順解析は次のように表される。

( )

D

=

F M ^(2.21)

逆解析は次のように表される。

1

( )

M = F

⁻

D

^(2.22)

ただし、Dはデータ空間を、Mはモデル空間を表す。F、

F

⁻¹はそれぞれ順解析演算子、

逆解析演算子と呼ばれる。

式(2.21)と(2.22)で示される関係は、Fig. 2.2のように表現することができる。

D M

−1

F

Fig.2.2 Illustration of forward and inverse mapping: the whole model space and data space are supposed available.

しかし、このような符号化は、実際の比抵抗法順解析と逆解析におけるいくつかの重要な前提条件が無視されている。すなわち、

① 応答は必ずソース位置と受信点位置に相関する。ソースがない場合にはモデルの応答は存在しない。

② 実際の測定は限られた部分空間でしか行われない。全データ空間を用いて逆解析する

(19)

ことはない。

③ 測定データには必ず誤差が含まれる。その誤差レベル以上になるモデル部分空間による応答だけが意味のある測定データとして反映される。

したがって、より精確には、次の符号化で順解析と逆解析が表現される。

順解析は：

*

( )

D = F M I X

(2.23)

逆解析は：

1 *

( )

F

⁻

D ⇒ M

^(2.24)

ただし、X は電極を配置することが可能な物理空間(相反定理により、送受信位置は等価性がある)であり、

D

^*は部分データ空間である。式(2.23)と(2.24)で表される関係は、Fig.2.3 のように表現できる。

Fig. 2.3には、順解析と逆解析の重要な特徴が示されている。すなわち、順解析は安定な

解析であり、モデル空間と電極配置空間とが重なる部分を人為的にコントロールすることができ、必要となる境界条件はすべて既知であるため、計算時間などの制限がなければ全データ空間Dを求めることも不可能ではない。一方、逆解析は、ノイズが含まれた部分データ空間から全モデル空間を推定しようとするため、不安定な解析となる。これは、コンピュータの能力を高くすることだけでは解決されない問題である。

?

D

X

M F

⁻¹

∩

M X

F

D

*

Fig. 2.3 Illustration of forward and inverse mapping: unstable process in the inverse mapping due to attempt to reconstruct the whole model space with only limited and noise contaminated data space being available.

(20)

2.4.2 最小2乗法を中心とする逆解析

比抵抗法において、代表的な逆解析のアルゴリズムは線形化最小 2 乗法である。その順解析計算に用いられるアルゴリズムには、電流保存則(式(2.9))を直接的に解く差分法と、電気エネルギー保存則のもとで変分原理を用いる有限要素法、および相加定理を用いる積分方程式法がある。どのアルゴリズムを用いる場合でも、あらかじめ地下空間を離散化し、

各セル(あるいはブロック、節点)の比抵抗をモデルベクトルの成分に設定する点では同様である。

比抵抗法では、モデルと応答との関係(式(2.23))が非線形の関係にあるため、最小 2 乗法による逆解析を行う場合には、局所の線形化と反復計算が必要となる。

具体的には、まず、式(2.23)を次の関数関係に書き換える。

( )

=

F

d m (2.25)

ただし、

(d_j) (j 1,2,..., )n

= =

d

( m

_i

) ( i 1,2,..., ) m

= =

m

d

と

m

は、それぞれデータベクトルとモデルベクトルである。今、非線形関数の式(2.25) をモデルベクトル初期値

m

₀の近傍、領域

( m

₀

+ ∆ m )

でTaylor展開する。

∆ m

が十分小さければ、

F ( m

₀

+ ∆ m )

の

∆ m

に関する線形部分は全体関数値の絶対優位の割合を締め、非線形部分を無視できる。すなわち、

(

0

)

F A

= + ∆

d m m

(2.26)

が近似的に成立する。

式(2.26)において

0

1 1 1

1 2

2 2 2

1 2

j( )

m

j x m

i

n n n

m

a a a

a a a F

A a

m

a a a

=

 

   ∂ 

 

=         =    = ∂   

L

L M M O M

L

m m m m

m

はF( )m の

m

₀でのヤコビアン、あるいは、感度行列である。

(21)

目的関数を次のように設定する。

( ^- )

²

S = W ∆ d A ∆ m

^* (2.27)

ただし、

∆ = − d d F ( m

₀

)

は残差ベクトルであり、

W

は

d

に含まれるノイズ

ε

^{のレベルに応}

じて決定される重み行列である。

最小2乗法の理論により、

ε = ( ε

_i

) ( i = 1,2,..., ) n

^{が偏りのない}^Gauss^{分布に従えば、目}

的関数Sが最小となるような

∆

m^*は、モデル修正量

∆ m

の偏りのない推定値である(中川・

小柳、1982)。

このことから

∆

m^*に関する次の正規方程式が導かれる。

( WA WA )

^T

∆ m

*

= ( WA W )

^T

d

(2.28)

感度行列Aの条件数が小さい場合は、式(2.28)から安定に

∆

m^*を求めることができる。しかし、もともと行列 Aは非線形問題を線形近似して導かれたものであり、モデルベクトル各成分の応答に対する感度は、地下空間の離散化の仕方と相関して大きな差が存在することが多い。よって、行列Aは悪条件になりやすく、式(2.28)を直接解くと安定な収束解は求められないことが多い。比抵抗法において、実際に用いられている逆解析法は、平滑化拘束最小2乗法である(例えば、Sasaki, 1989)。すなわち、最小にする目的関数を次のように設定する。

( ^- )

² ² ²

S = W ∆ d A ∆ m

^*

+ α C ∆ m

^*

ただし、行列

C

は、比抵抗分布を平坦化(1 次差分の場合)あるいは平滑化(2 次差分の場合) するための行列である。

α

は、平滑の程度を決める平滑化パラメータ(ラグランジェ係数) である。

結局、各反復計算の過程で解く方程式は次式で表される。

2 *

[( WA WA )

^T

+ α C C

^T

] ∆ m = ( WA A )

^T

∆ d

(2.29)

この従来の逆解析法による比抵抗法 3 次元解析では、次の問題について、まだ解決には至ってない。

① 反復計算のたびに感度行列を計算すると非常にコストの高い計算となる。

② 感度行列の近似計算は平坦な地形でなければ成立しないため、起伏地形の解析を行うことはできない。

③ データに含まれる情報量と逆解析で推定しようとするモデル情報とが矛盾関係にあるため、解析は非常に不安定となる。安定化操作は解像度の低下を招く。

(22)

2.4.3 順解析・逆解析と境界値問題

順問題と逆問題の唯一性の検討は、安定的な解析手法を開発するための重要な指針を導くものである。それは境界値問題の検討と密接な関係がある。

比抵抗法順問題の唯一性については、Morse(1953)の詳細な検討がある。それによると、

外境界面を

Γ

^{とする任意}³^{次元区間領域}

Ω

で楕円型偏微分方程式(式(2.9)を含む)が成立し、

閉合境界面

Γ

^{上で、分段的に}^Dilichlet条件(式(2.10)を含む)あるいはNeumann条件(式(2.11) を含む)あるいはそれらの組み合わせが与えられていれば、領域

Ω

^{内(境界面}

Γ

^{を含む)の未}

知関数

U

は、1 個の定数Kの差が許される以外は唯一に決定される。もし、境界

Γ

^上少な

くとも1点でDilichlet条件が与えてあれば、その定数は

K = 0

となる。この結論は、楕円型

偏微分方程式境界値問題の唯一性定理と呼ばれる。

ここで、逆問題と相関する唯一性問題を、次のように定義する。

外境界面を

Γ

^{とする任意}3次元区間領域

Ω

上で問題を考える。境界面

Γ

^上、

U

についてある境界条件

Λ

¹_Γ^と

Λ

_Γ²がそれぞれ与えられた場合、

Ω

上でそれぞれ式(2.9)を満足するような

σ

₁^と

σ

2が求められるとする。このとき、唯一性問題は次の2つ関係が成立するか否かを問う問題である。

1 2

σ σ

Γ Γ

Λ = Λ ⇒ =

^(2.30a)

1 2

σ σ

Γ Γ

Λ ≠ Λ ⇒ ≠

^(2.30b)

この問題について、Sylvester and Uhlmann (1986)は、もし

σ

^が

Ω

^{上で円滑(空間的}²^階微係数が存在する)であれば、式(2.30)が成立することを証明した。Nachman(1988)はさらに進んで、もし

σ

^が

Ω

上で有限であれば、式(2.30)が成立するという境界値逆問題の唯一定理を証明した。

これらの唯一定理は、逆解析法に関して、いくつかの重要なポイントを提示している。

例えば、

① 全データ空間が分からなくても、調査したい領域をなるべく囲むような境界に測定点を配置すれば、逆解析の安定性と唯一性が向上する。

② それが達成できない場合は、既知のすべての情報を用いて、データの欠けた境界部分にできるだけ良い近似条件を与える。

③ 既知データではコントロールできないモデルの部分空間に対して、無駄な推定計算をしない。

(23)

2.5 最新の研究動向

最近、孔井と地表を同時に利用した比抵抗トモグラフィ測定がよく行われるようになった。そのための高密度測定に対応できる探査機器の開発や解析理論の研究が注目される(例えば、Herwanger et al(2002)、Lesur et al(1999)、Zhou and Greenhalgh(2000))。

一方、医療分野では、人体内部の精密な比抵抗構造を調べるための、電気インピーダンストモグラフィ(EIT) の研究進展が著しい(例えば、Cheney et al (1999)、Molinari et al (2001))。

EIT法では、有限な電極数で高解像度を獲得するために、最適な電極位置についての研究に力が注がれ、とくに、逆解析の唯一性問題についての研究が進展している。

比抵抗法 3 次元逆解析の研究においては、感度行列の効率的な計算方法の研究が依然として大きな割合を占めている。Li and Oldenburg (2000)は、感度計算において、近似計算と完全に求める計算とを組み合わせたアルゴリズムが計算効率と解析精度の両面で有利であることを示した。Zhou and Greenhalgh(2002))は、感度計算に解析式を用い、まず近似法により初期モデルを与え、続いて感度を精確に求めていく逆解析アルゴリズムを示した。

一方、測定データから解析領域内部の比抵抗境界を推定する研究も行われている。とくに注目されるのは、EIT法におけるMolinari et al (2001)の、解析領域の離散化格子を解析の進行に伴って修正していくアルゴリズムである。彼らのアルゴリズムは、解析のスタート時点では粗い分割密度を用い、解析過程の中で、比抵抗の傾度に相関する基準関数を用いて、部分領域をさらに細かく分割するか否かを決定していく。そのような変動的な離散化方法は、流体力学で非線形問題を解く場合によく使われる手法であり、彼らは、計算効率と解像度の向上に有利であると報告している。

D

X

−1

F

≈ ∩

M M X *

D

Fig. 2.4 Illustration of forward and inverse mapping: increasing the robustness of the inverse mapping by increasing product space of model and physical space related to the partial data space.

このように、最新の研究動向には、2つの傾向がある。1つは、測点の配置を適切化することにより、データの部分空間にできるだけ全モデル空間情報を引き込もうとする、測定

(24)

法の開発(Fig.2.4で示す)であり、他は、解析の過程中に自動的に、適切な解析領域の設定および離散化を果すことによって、限られたデータから推定可能なモデル空間情報だけを取り出そうとする、解析法の研究(Fig.2.5)である。

D X

M

₋₁

F

∩

M X

F

D

*

Fig.2.5 Illustration of forward and inverse mapping: increasing the robustness of the inverse mapping by trying to reconstruct the partial model space that is sufficient sensitive to the partial data space.

2.6 おわりに

これまで述べてきたように、比抵抗 3 次元探査法を実用化するために、信頼性、効率性と解像度を追求した 3 次元数値解析法の研究と開発が要求されている。それらの背景と現状を踏まえて、次の事項について検討する必要がある。

① 解析の精度と効率を向上させるための空間離散化の方法

② 任意地形、任意地下比抵抗構造に対応できる順解析法

③ 解析の安定性と唯一性を向上させるためのデータの取得法

④ 信頼性、効率性と解像度を追求した逆解析法

本研究では、これらの事項について検討するために、まず従来の 3 次元数値解析法の問題点を分析し、最適な空間離散化の方法と最適な計算空間理論を検討した上で、数値座標変換による新たな3次元数値解析法を導いた。

(25)

参考文献

Cheney, M. Isaacson, D. and Newell, J. C. (1999): Electrical Impedance Tomography, SIAM Review, 41, 85-101.

Coggon, J. H. (1971): Electromagnetic and electrical modeling by finite-element method, Geophysics, 36, 132-155.

Dakhnov, V. N. (1941): Electrical well logging interpretation of electric log, Moscow.

Dey, A. and Morrison, H. F., (1979): Resistivity modeling for arbitrarily shaped three-dimensional structures, Geophysics, 44, 753-780.

Frank, P. and Von Mises, R. (1935): Die differential and integralgliechung der mechanik und physik, Braunschweig, Friedn., Vieweg u. Sohn., 2, 1106

Ghosh, D. P. (1971): The application linear filter theory to the direct interpretation of geo-electrical resistivity sounding measurements, Geophys. Prosp., 19, 192-217

Herwanger, J. V., Pain, C. C. and de Oliveira, C. R. E. (2002): Electic and Seismic Inversion in Anisotropic Inhomogeneous Media, Proceeding of 4^th international conference on Inverse Problems in Engineering.

Hohmann, G. W. (1975): Three-dimensional induced polarization and electromagnetic modeling, Geophysics, 40, 309-324.

Inman, J. R., Ryu, J. and Ward, S. H. (1973): Resistivity inversion, Geophysics, 38, 1088-1108.

Inman, J. R. (1975): Resistivity inversion with ridge regression, Geophysics, 40, 798-817.

Jepson, A. F. (1969): Resistivity and induced polarization modeling, Ph.D. thesis, University of California, 5, 333-352.

Lesur V., Cuer M. and Straub, A. (1999): 2-D and 3-D interpretation of electrical tomography measurements, Part 2: The inverse problem.

Li, Y. and Oldenburg, D. W. (1994): Inversion of 3-D resistivity data using an approximate inverse mapping, Geophys. J. Int., 116, 527-537.

Li, Y. and Oldenburg, D. W. (2000): 3-D inversion of induced polarization data, Geophysics, 65, 1931-1945.

Molinari, M., Cox, S. J., Blott, B. H. and Daniell, G. J. (2001): Efficient Non-Linear 3D Electrical Tomography Reconstruction, Proceeding of 2nd World Congress on Industrial Process Tomography, 424-432.

Morse, P. M. and Feshbach, H. (1953): Method of theoretical physics, McGRAW-HILL Book Comp.

Nachman, A. (1988): The inverse reconstructions from boundary measurements, Ann. of Math., 128, 531-587.

中川徹・小柳義夫 (1982): 最小 2 乗法による実験データ解析, 東京大学出版会, 29-53.

(26)

Park, S. K. and Van, G. P. (1991): Inversion of pole-pole data for 3-D resistivity structure beneath arrays of electrodes, 56, 951-960.

Pridmore, D. F. (1978): Three-dimentional modeling of electric and electromagnetic data using the finite-element method, Ph.D. thesis, Univ. of Utah.

Sasaki, Y. (1989): Two-dimensional joint inversion of magnetotelluric and dipole-dipole resistivity data, Geophysics, 54, 254-262.

Sasaki, Y. (1994): 3-D resistivity inversion using the finite-element method, Geophysics, 59, 1839-1848.

佐々木裕 (1981): 比抵抗垂直探査における 2 次元構造の自動解析(I), 物理探鉱, 34, 341-350.

佐々木裕 (1992): 比抵抗法とIP法の3次元インバージョン, 物理探査, 45, 3-9.

Schlumberger, C. (1920): Etudes sur la prospection Erectrique du sours sol, Gauthiers-Villars, Paris.

Slichter, L. B. (1933): The interpretation of the resistivity prospecting method for horizontal structures, Physics, 4, 307-322.

Stefanesco, S. (1950): Modeles theoreques de milieux heteirogenes pour les methods de prospection electrique a courants stationnaires, Comitetul Geologic, Studii technice si economice, seria D., Nr. 2, and Nationola Bucaresti

Stefanesco, S., Schlumberger, C., and Schlumberger, M. (1930): Sur la distribution electrique potentielle autour d’une prise de terre ponctuelle dans un terrain a couches horizontals homogenes et isotropes, J. Phys. Et Red., Ser. 7, 132-141

Sylvester, J. and Uhlmann, G. (1986): A uniqueness theorem for an inverse boundary value problem in electrical prospection, Comm. Pure Appl. Math., 39, 91-112.

Tagg, G. F. (1930): The earth resistivity method of geophysical prospecting: Some theoretical considerations, Mining Mag., 43, 150-158.

Wenner, F. (1912): The four-terminal conductor and the Thomson bridge, U. S. Bulletin of Bureau of Standards, V8, 559-610.

Zhou, B. and Greenhalgh, S. A. (2000): Crosshole resistivity tomography using different electrode configuration, Geophys. Prosp., 48,443-468.

Zhou, B. and Greenhalgh, S. A. (2002): Rapid 2-D/3-D crosshole resistivity imaging using the sensitivity function, Geophysics, 67, 755-765.

(27)

第 3 章比抵抗法における数値座標変換の基礎理論

3.1 はじめに

この章では、比抵抗法の 3 次元解析における主要な問題点を分析し、最適な空間離散化および最適な計算空間について検討する。

2次元解析から3次元解析へ移行すると、次の2つの状況について大きな変化が生じる。

すなわち、未知変数の数が飛躍的に増大することと、パラメータに対する応答感度の距離減衰が著しく速くなることである。そのため、3 次元解析においては、とくにその効率性、

安定性について十分に検討し、確認しておかなければならない。

解析精度や解析効率が空間離散化の仕方と密接な関係を持つことは、多くの研究者によって指摘されてきた。例えば、Fox(1979)は、差分法による楕円型偏微分方程式数値解法について、特異点付近での離散化誤差を検討し、特異点に近づくにつれて、セルのサイズを急激に小さくすることを提案した。一方、Dey and Morrison (1979) は、非一様なセル分割による差分法を用いて、比抵抗法 3 次元モデリングの効率化をはかった。しかしながら、

Pridmore et al. (1981)は、そのような非対称セルによる離散化の方法は、誤差を増大させる要因となることを指摘した。

明らかに、フォワード計算(順解析)では、地下を必要以上に高い密度で分割すると計算に無駄が生じる。一方、比抵抗変化の激しい箇所や電流ソースのような特異点付近で分割密度が不足すると、計算精度の低下をまねくことになる。インバージョン(逆解析)では、分割密度の足りない部分でのセルの感度は必然的に大きくなり、逆に必要以上にセル密度の高い部分の感度は小さくなるため、悪条件の感度行列を生み出す要因となる。

このように、従来の解析アルゴリズムをそのまま 3 次元解析に適用した場合には、計算に無駄が生じやすく、モデルが大規模になるほど安定性が低下する傾向にある。これらの原因はすべて、その物理問題に最適な空間離散化を考慮していないためであると考えられる。

空間を最適に離散化する方法は、支配方程式と分けて考えることはできず、また、モデルごとに異なってくるため、そのモデルの解が得られるまでは分からないというのが実情である。そのため、空間の最適な離散化については、これまであまり検討されてこなかった。

すでに明らかなように、空間離散化の最適化は、順問題と逆問題とでは異なる意味にあるが、最終的には、統一されなくてはならない。順問題における最適化とは、数値解の離散化誤差を抑えながら空間領域のセル分割数を最小に抑えることである。すなわち、場の

数値座標変換による 比抵抗法 3 次元データ解析に関する研究