特許 6912487 放射性崩壊による熱産生の計算方法

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】6912487

(24)【登録日】2021年7月12日

(45)【発行日】2021年8月4日

(54)【発明の名称】放射性崩壊による熱産生の計算方法

(51)【国際特許分類】

   G01V 99/00 20090101AFI20210727BHJP

【ＦＩ】

   G01V99/00

【請求項の数】16

【全頁数】24

(21)【出願番号】特願2018-544153(P2018-544153)

(86)(22)【出願日】2017年2月20日

(65)【公表番号】特表2019-507347(P2019-507347A)

(43)【公表日】2019年3月14日

(86)【国際出願番号】NO2017050044

(87)【国際公開番号】WO2017142422

(87)【国際公開日】20170824

【審査請求日】2020年2月18日

(31)【優先権主張番号】1602935.7

(32)【優先日】2016年2月19日

(33)【優先権主張国】GB

(73)【特許権者】

【識別番号】515339653

【氏名又は名称】エクイノール・エナジー・アーエス

【氏名又は名称原語表記】ＥＱＵＩＮＯＲ ＥＮＥＲＧＹ ＡＳ

(74)【代理人】

【識別番号】100110423

【弁理士】

【氏名又は名称】曾我 道治

(74)【代理人】

【識別番号】100111648

【弁理士】

【氏名又は名称】梶並 順

(74)【代理人】

【識別番号】100166235

【弁理士】

【氏名又は名称】大井 一郎

(74)【代理人】

【識別番号】100179936

【弁理士】

【氏名又は名称】金山 明日香

(74)【代理人】

【識別番号】100195006

【弁理士】

【氏名又は名称】加藤 勇蔵

(72)【発明者】

【氏名】ホクスタド、ケティル

(72)【発明者】

【氏名】デュフォート、ケネス

(72)【発明者】

【氏名】フィシュラー、クリスティン

(72)【発明者】

【氏名】キルケボ、ルネ

(72)【発明者】

【氏名】アラソナティ・タサロヴァ、ズザナ

【審査官】 佐野 浩樹

(56)【参考文献】

【文献】 米国特許出願公開第２００４／０２２０７４０（ＵＳ，Ａ１）

【文献】 米国特許出願公開第２０１５／０２４２３６２（ＵＳ，Ａ１）

【文献】 特表２０１１−５０８８７６（ＪＰ，Ａ）

【文献】 特開２００４−２９２８８０（ＪＰ，Ａ）

【文献】 特表２０１１−５０９４１２（ＪＰ，Ａ）

【文献】 米国特許出願公開第２０１０／０１８５４２２（ＵＳ，Ａ１）

【文献】 米国特許出願公開第２０１０／０３２６６６９（ＵＳ，Ａ１）

【文献】 RYBACH, Ladislaus et al.，THE VARIATION OF HEAT GENERATION, DENSITY AND SEISMIC VELOCITY WITH ROCK TYPE IN THE CONTINENTAL LITHOSPHERE，TECTONOPHYSICS，ELSEVIER，１９８４年，Vol. 103, No. 1-4，p. 335-344

【文献】 MUKAI, Makiko, et al.，Measurement of radioactive heat generation in rocks by means of gamma ray spectrometry[online]，Proceedings of the Japan Academy, Series B，日本，The Japan Academy，１９９９年，Vol. 75, No. 7，p. 181-185，[検索日: 2021.1.13], インターネット<https://www.jstage.jst.go.jp/article/pjab1977/75/7/75_7_181/_article>で入手可能。<https://doi.org/10.2183/pjab.75.181>

(58)【調査した分野】（Int.Cl.，ＤＢ名）

Ｇ０１Ｖ１／００−９９／００

(57)【特許請求の範囲】

【請求項1】

地球物理学的構造の放射性崩壊による熱産生（ＲＨＰ）を計算する方法であって、前記地球物理学的構造の少なくとも１つの地球物理学的パラメータが提供され、前記少なくとも一つの地球物理学的パラメータと前記ＲＨＰとの間に統計的依存性が存在し、
前記方法は、前記地球物理学的構造の前記ＲＨＰを推定するために、前記少なくとも１つの地球物理学的パラメータをインバージョンさせることであって、前記地球物理学的パラメータの前記ＲＨＰは、前記少なくとも１つの地球物理学的パラメータの少なくとも１つの原因の１つである、ステップを含み、
前記少なくとも１つの地球物理学的パラメータの前記インバージョンのステップは、前記少なくとも１つの地球物理学的パラメータと、前記地球物理学的構造の前記ＲＨＰとの間の関係を定義する、フォワードモデルを選択することを含み、
前記少なくとも１つの地球物理学的パラメータは、密度、地震波速度、磁化率、電気伝導率、電気抵抗率又は残留磁気の少なくとも１つを含む、地球物理学的構造の放射性崩壊による熱産生（ＲＨＰ）を計算する方法。

【請求項2】

前記地球物理学的構造のサンプルからの前記少なくとも１つの地球物理学的パラメータの少なくとも１つの測定と、前記地球物理学的構造の前記ＲＨＰとを含む、較正データが提供され、
前記方法は、前記較正データに基づいて、前記フォワードモデルを最適化するステップを含む、請求項１に記載の方法。

【請求項3】

前記地球物理学的構造の少なくとも２つの地球物理学的パラメータが提供され、
前記方法は、地球物理学的構造の前記ＲＨＰを推定するために、前記少なくとも２つの地球物理学的パラメータをインバージョンさせるステップを含む、請求項１又は２に記載の方法。

【請求項4】

前記インバージョンさせるステップは、前記少なくとも２つの地球物理学的パラメータ間に統計的独立性が存在し、且つ、各地球物理学的パラメータと前記地球物理学的構造の前記ＲＨＰとの間に統計的依存性が存在する、モデルを使用することを含む、請求項３に記載の方法。

【請求項5】

前記少なくとも２つの地球物理学的パラメータは、少なくとも１つの電磁地球物理学的パラメータと、少なくとも１つの力学的な地球物理学的パラメータとを含むことができる、請求項３又は４に記載の方法。

【請求項6】

前記インバージョンさせるステップは、各地球物理学的パラメータに対してフォワードモデルを用いることを含み、
前記フォワードモデルの各々は、前記各地球物理学的パラメータと、前記地球物理学的構造の前記ＲＨＰとの間の関係を定義する、請求項３、４又は５に記載の方法。

【請求項7】

前記各フォワードモデルのそれぞれにおける唯一の変数は、前記地球物理学的構造の前記ＲＨＰである、請求項６に記載の方法。

【請求項8】

前記地球物理学的構造の１つ又は複数のサンプルからの前記少なくとも２つの地球物理学的パラメータの各々の少なくとも１つの測定と、前記地球物理学的構造の前記ＲＨＰとを含む、較正データが提供され、
前記方法は、前記較正データに基づいて、前記各フォワードモデルを最適化するステップを含む、請求項６又は７に記載の方法。

【請求項9】

前記地球物理学的少なくとも１種類の地球物理学的データが提供され、
前記方法は、前記少なくとも１つの地球物理学的パラメータを計算するために、前記少なくとも１種類の地球物理学的データをインバージョンさせるステップを含む、請求項１〜８のいずれか一項に記載の方法。

【請求項10】

前記地球物理学的データを取得するステップを含む、請求項９に記載の方法。

【請求項11】

ＲＨＰの空間依存関数を計算する方法であって、
前記地球物理学的構造上の複数の点に対して点ごとに、請求項１〜１０のいずれか一項に記載の方法を実行することにより、前記複数の点の各々において前記ＲＨＰを計算するステップと、
ＲＨＰの前記空間依存関数を構築するステップと
を含む、ＲＨＰの空間依存関数を計算する方法。

【請求項12】

地球物理学的構造における基礎熱流量を計算する方法であって、マントルからの熱流量寄与が提供され、前記方法は、
地球物理学的構造の空間にわたって、請求項１〜１０のいずれか一項に記載の地球物理学的構造のＲＨＰを計算する方法又は請求項１１に記載のＲＨＰの空間依存関数を計算する方法を用いて、前記地球物理学的構造のＲＨＰを計算するステップと、
前記地球物理学的構造の空間の少なくとも一部の前記ＲＨＰを合計するステップと、
前記ＲＨＰの前記合計値を前記マントルからの前記熱流量寄与に加算するステップと
を含む、地球物理学的構造における前記基礎熱流量を計算する方法。

【請求項13】

地表上の地表熱流量を計算する方法であって、前記地表熱流量に対する堆積物寄与が提供され、前記方法は、
請求項１２に記載の地球物理学的構造における前記基礎熱流量を計算する方法を用いて、前記基礎熱流量を計算することであって、前記地球物理学的構造は前記地殻又はリソスフェアであり、前記基礎熱流量は前記地殻又はリソスフェアの最上部で計算される、ステップと、
前記堆積物寄与を前記基礎熱流量に加算するステップと
を含む、地表上の地表熱流量を計算する方法。

【請求項14】

請求項１３に記載の方法を用いて計算された前記地表熱流量を用いるステップを含む、前記地球物理学的構造の温度を計算する方法。

【請求項15】

請求項１〜１４のいずれか一項に記載の方法を含む、地球物理学的構造の熱及び／又は温度モデルを生成する方法。

【請求項16】

コンピュータ上で実行された際に、プロセッサに請求項１〜１５のいずれか一項に記載の方法を行わせるように構成されるコンピュータ可読命令を含む、コンピュータプログラムプロダクト。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、地球物理学的構造における、放射性崩壊による熱産生、基礎熱流量及び地表熱流量を計算する方法を提供する。

【背景技術】

【0002】

地質年代による地表下温度の変化は、堆積盆の石油のプロスペクティビティにとって非常に重要である。熱履歴は、根源岩の成熟及び貯留岩の質を制御する。クオーツのセメンテーションに関連した間隙率の減少は、温度に直接依存する。堆積盆の熱レジームは、２つの主要因である、地球のマントル及び地殻からの寄与を有する基礎熱流量と、堆積シーケンスの熱伝導率プロファイルとによって制御される。

【0003】

深さ及び地質年代にともなう温度分布は、主に、これらの一次影響間の相互作用によって決定される。熱的モデリングは、堆積盆モデリング（非特許文献１）の一部である。従来の堆積盆モデリングは、地質学的分野及び地球化学からの入力によって推進される。しかしながら、本方法は、有利に、地球物理学的データによって推進される。

【0004】

先行技術では、放射性崩壊による熱産生（ＲＨＰ：Radiogenic Heat Production）は、一般的に、スペクトルガンマ線記録を用いて、ウラン、トリウム及びカリウム（これらは、ＲＨＰの主因である）などの放射性元素の量を推定することによって求められる。例えば、（非特許文献２）において教示される方法では、放射性崩壊の直接測定を用いてＲＨＰを計算する。この方法は、ボアホールで測定されたスペクトルガンマ線記録に基づいている。これは標準的なアプローチであり、（非特許文献３）によって示される経験式を使用する。この方法において、ウラン、トリウム及びカリウムからの寄与に関する分離した情報を取得するためには、種々の元素が種々の周波数帯で放射することから、ガンマ線のスペクトル（周波数依存）記録が必要とされる。

【0005】

代替的には、手作業によって、堆積盆モデリング中に温度データに適合するようにＲＨＰ値を調整することが知られている。これらの方法は、不正確である傾向があるとともに、ＲＨＰを計算するという特定の目的のためにデータが収集されることを必要とする。

【0006】

先行技術の別の例として、（非特許文献４）は、地震波速度を用いて、地殻の温度を推定すること（圧力が既知であるという仮定を用いると共に、圧力及び温度を地震波速度にリンクさせる状態方程式を用いる）を教示する。この後、熱伝導方程式が解かれ、地殻の温度が地表熱流量及び地殻のＲＨＰに関連付けられる。（非特許文献５）は、概ね類似の方法を教示する。

【先行技術文献】

【非特許文献】

【0007】

【非特許文献1】Ａｌｌｅｎ ａｎｄ Ａｌｌｅｎ，２００５

【非特許文献2】１４ｔｈ Ｉｎｔｅｒｎａｔｉｏｎａｌ Ｃｏｎｇｒｅｓｓ ｏｆ ｔｈｅ Ｂｒａｓｉｌｉａｎ Ｇｅｏｐｈｙｓｉｃａｌ Ｓｏｃｉｅｔｙ，２０１５，Ｏｌｉｖｅｉｒａ ｅｔ ａｌ，“Ｒａｄｉｏｍｅｔｒｉｃ ａｎｄ ｔｈｅｒｍａｌ ｓｉｇｎａｔｕｒｅｓ ｏｆ ｔｕｒｂｉｄｉｔｅ ｆｌｏｗｓ ｉｎ Ｎａｍｏｒａｄｏ ｏｉｌ ｆｉｅｌｄ”，ｐａｇｅｓ ７９９−８０３

【非特許文献3】Ｒｙｂａｃｈ（１９８６）

【非特許文献4】Ｇｅｏｃｈｅｍｉｓｔｒｙ Ｉｎｔｅｒｎａｔｉｏｎａｌ，Ｖｏｌ．４４，Ｎｏ．１０，２００６，Ｋｒｏｎｒｏｄ ｅｔ ａｌ，“Ｄｅｔｅｒｍｉｎｉｎｇ ｈｅａｔ ｆｌｏｗｓ ａｎｄ ｒａｄｉｏｇｅｎｉｃ ｈｅａｔ ｇｅｎｅｒａｔｉｏｎ ｉｎ ｔｈｅ ｃｒｕｓｔ ａｎｄ ｌｉｔｈｏｓｐｈｅｒｅ ｂａｓｅｄ ｏｎ ｓｅｉｓｍｉｃ ｄａｔａ ａｎｄ ｓｕｒｆａｃｅ ｈｅａｔ ｆｌｏｗｓ”，ｐａｇｅｓ １０３５−１０４０

【非特許文献5】Ｐｈｙｓｉｃｓ ｏｆ ｔｈｅ ｓｏｌｉｄ ｅａｒｔｈ，Ｖｏｌ．４３，Ｎｏ．１，２００７，Ｋｒｏｎｒｏｄ ｅｔ ａｌ，“Ｍｏｄｅｌｉｎｇ ｏｆ ｔｈｅ ｔｈｅｒｍａｌ ｓｔｒｕｃｔｕｒｅ ｏｆ ｃｏｎｔｉｎｅｎｔａｌ ｌｉｔｈｏｓｐｈｅｒｅ”，ｐａｇｅｓ ９１−１０１

【発明の概要】

【課題を解決するための手段】

【0008】

ある局面では、本発明は、地球物理学的構造における放射性崩壊による熱産生（ＲＨＰ）を計算する方法であって、地球物理学的構造の少なくとも１つの地球物理学的パラメータが提供され、地球物理学的構造のＲＨＰを推定するために、少なくとも１つの地球物理学的パラメータをインバージョンすることを含む方法を提供する。

【0009】

本発明者らは、地球物理学的構造のＲＨＰは、密度、地震波速度（好ましくは、地震ｐ波速度）又は磁化率などの地球物理学的構造の少なくとも１つの地球物理学的パラメータをインバージョンすることによって求めることができることを発見した。この方法は、それが地球物理学的構造全体にわたってＲＨＰを求めるために地球物理学的構造の詳細なサンプリングを必要としないので、有利である。より正確に言えば、表面測定及び観測により、地球物理学的構造全体にわたって、地球物理学的パラメータ値を取得することができるので、本発明の方法は、表面測定及び観測を用いて、ＲＨＰが計算されることを可能にする。

【0010】

ガンマ線記録の使用を必要とする上述の先行技術の方法（上述のＯｌｉｖｅｉｒａらなど）とは異なり、本発明は、ＲＨＰを計算するために、単に、地球物理学的パラメータ値（これらは既に／簡単に入手可能なものでもよく、またボアホール内ではなく地表から取得されてもよい）に頼ることができる。さらに、地震波速度をＲＨＰ推定に変換する際に多数の中間ステップを伴う上述のＫｒｏｎｒｏｄらの先行技術とは異なり、本方法は、インバージョンを用いて、１つ又は複数の地球物理学的パラメータから直接ＲＨＰを計算する。

【0011】

本方法によって計算されるＲＨＰは、好ましくは、現在のＲＨＰであることができる。これは、地球物理学的パラメータの現在の値が使用される場合に当てはまり得る。もちろん、必要に応じて、履歴値を使用して、履歴ＲＨＰを計算することもできる。地球物理学的構造は、好ましくは、地球の地殻又はリソスフェアであることができる。地球の地殻は、マントルと地表又は堆積物層との間の層である。地球のリソスフェアは、延性マントルと地表又は堆積物層との間の層である（すなわち、リソスフェアは、地殻及び（大部分は）脆性の上部マントルを包含する）。地殻又はリソスフェアのＲＨＰを計算することは、下記にさらに説明されるように、このＲＨＰ値を用いて、地表熱流量を計算することができるので、特に有用である。地表熱流量の計算を用いて、堆積層又は地殻中に、炭化水素の生成にとって適切な条件が存在するか、又は存在したか否かを評価する際に重要な温度分布を計算することができる。

【0012】

インバージョンすること、又はインバージョンは、当該分野においてよく知られている用語である。これは、少なくとも１つの観測／測定パラメータから、パラメータの原因（又はパラメータの複数の原因の少なくとも１つ）を計算するプロセスを表現する。従って、このケースでは、物理的に言って、ＲＨＰは地球物理学的パラメータに影響を与える。しかしながら、測定されるのは、地球物理学的パラメータであって、ＲＨＰではない。従って、地球物理学的パラメータからＲＨＰを計算することは、インバージョンすることであると表現することができる。インバージョンは、１つ又は複数の地球物理学的パラメータから直接ＲＨＰ値を計算するために、１つ又は複数のモデル（下記に説明されるフォワードモデルなどの、１つ又は複数の地球物理学的パラメータをＲＨＰに関連付ける１つ又は複数の岩石物理学モデルなど）を用いる計算であると見なすことができる。

【0013】

本明細書全体を通して、「計算すること」及び「推定すること」などの用語が使用される。これらは、限定することを意図したものではなく、これらは、単に、ＲＨＰなどの実際の物理値（又は少なくとも物理値の（近い）近似）の値を決定すること、又は取得することを意味することが意図されている。

【0014】

地球物理学的パラメータは、密度、磁化率、地震波速度（好ましくは、地震ｐ波速度）、電気伝導率、抵抗率又は残留磁気などの地球物理学的構造の任意の特性であることができる。特に、地球物理学的パラメータは、そのような特性を表す１つ又は複数の値であることができる。特に、地球物理学的パラメータは、ＲＨＰによって影響を受ける任意のそのような特性であることができる。

【0015】

インバージョンのステップは、少なくとも１つの地球物理学的パラメータと、地球物理学的構造のＲＨＰとの間の関係を定義する、フォワードモデルを選択することを含むことができる。

【0016】

インバージョンの計算において、フォワードモデルは、既知／測定パラメータ（この場合、地球物理学的パラメータ）と、未知の数量（この場合、ＲＨＰ）との間の関係である。

【0017】

フォワードモデルは、関連の地球物理学的パラメータをＲＨＰに関連付ける期待される傾向に基づいて選択することができる。例えば、地球物理学的パラメータは、通常、ＲＨＰの増加に伴って、（地球物理学的パラメータに応じて）増加又は減少し得る。地球物理学的パラメータが密度又は地震波速度（好ましくは、地震ｐ波速度）又は伝導率である場合、地球物理学的パラメータは、ＲＨＰの増加に伴って減少し得る。地球物理学的パラメータが磁化率である場合、地球物理学的パラメータは、ＲＨＰの増加に伴って増加し得る。使用される特定のフォワードモデルは、本発明に必須ではなく、幾つかのそのようなフォワードモデルが、当該分野で知られている。また、当業者は、どの１つ又は複数のモデルが使用され得るかに気付くであろう。実際、フォワードモデルが地球物理学的パラメータとＲＨＰとの間の一般傾向をモデリングすることができる限り、種々のフォワードモデルを用いて類似の結果を達成することができる。

【0018】

地球物理学的パラメータがＲＨＰの増加に伴って減少する場合、モデルは任意の崩壊関数であることができる。例えば、地球物理学的パラメータは、ＲＨＰ（Ａ）の逆数又はＲＨＰの逆数の（自然）対数に比例してもよく、すなわち、

【数1】

又は

【数2】

である。ただし、ａは定数である。地球物理学的パラメータがＲＨＰの増加に伴って増加する場合、モデルは任意のシグモイド関数であることができる。

【0019】

地球物理学的パラメータとＲＨＰとの間の関係の説明例として、Ｒｙｂａｃｈ（１９７８）は、ＲＨＰと地震ｐ波速度との間の指数関係

【数3】

を提案している。ただし、ＡはＲＨＰであり、ｖ_ｐは地震ｐ波速度であり、ａ及びｂは定数である。もちろん、これは単なる説明例であり、ＲＨＰと地震波速度との間の他の関係を使用することもできる。

【0020】

従って、上記から理解できるように、正確なフォワードモデルは、岩石物理学関係の知識に基づいて、当業者によって選択されることが可能である。

【0021】

好ましくは、ＲＨＰと地球物理学的パラメータとの間のモデル関係は、その他の地球物理学的パラメータなどのその他の変数に依存しない。もちろん、他の定係数が存在してもよいが、好ましくは、１つの変数のみが存在する。例えば、ＲＨＰと地震ｐ波速度との間の例示的なＲｙｂａｃｈ（１９７８）の関係に再度注目すると、ＲＨＰが依存する唯一の変数は、地震ｐ波速度である。Ｒｙｂａｃｈ（１９７８）の上記の式の他の係数（ａ及びｂ）は、単なる定数である。下記に説明されるように、定係数はデータの較正によって求められてもよい。

【0022】

１つ又は複数のモデルは、システムの完全な複雑さを示さない場合があること、すなわち、地球物理学的パラメータがＲＨＰにのみ依存するようにモデルが意図的に単純化される場合があることを理解されたい。実際には、１つ又は複数の地球物理学的パラメータは、一般に多くの変数に依存する。しかしながら、本方法で使用される１つ又は複数のモデルでは、１つ又は複数の地球物理学的パラメータは、関心のある変数、この場合ＲＨＰにのみ依存することができる。

【0023】

地球物理学的構造のサンプルからの少なくとも１つの地球物理学的パラメータの少なくとも１つの測定と、地球物理学的構造のＲＨＰとを含む、較正データが提供されてもよい。この方法は、較正データを取得することをさらに含むことができる。較正は、好ましくは、地球物理学的構造のサンプルからの少なくとも１つの地球物理学的パラメータの複数の測定と地球物理学的構造のＲＨＰとを含むことができる。少なくとも１つの地球物理学的パラメータの少なくとも１つの測定と、ＲＨＰの測定とは、好ましくは、サンプルの実質的に同じ位置から獲得されたものであることができ、又はサンプル全体としての全体的／平均測定であることができる。

【0024】

本方法におけるインバージョンのステップは、較正データに基づいて、フォワードモデルを最適化することを含むことができる。この最適化は、較正データを使用して、フォワードモデルの定係数の最適値を求めることを含むことができる。一般的に、較正データの量が多いほど、最適化がより良くなる。

【0025】

Ｒｙｂａｃｈ（１９７８）に提示された地球物理学的パラメータとＲＨＰとの間の例示的な関係に再度注目すると、較正データを用いて求めることができるのは、係数ａ及びｂである。

【0026】

フォワードモデルを最適化するために、提供された地球物理学的パラメータに対して、フォワードモデル（これは所与のＲＨＰから地球物理学的パラメータを計算する）が、ある特定の誤差分布を有する（すなわち、提供／観測された地球物理学的パラメータと、各フォワードモデルによって計算された地球物理学的パラメータとの間の差が、誤差分布を与える）と仮定されることができる。好ましくは、誤差分布は、好ましくは平均が０のガウス誤差分布であると仮定される。フォワードモデルは、例えば、可能な限り０に近くなる誤差分布の平均を有することによって、及び可能な限り誤差分布の小さな分散を有することによって、誤差分布が可能な限り小さくなるように誤差分布を減少させることによって、最適化される。最適化は、フォワードモデルを最適化するフォワードモデルの１つ又は複数の定係数（Ｒｙｂａｃｈ（１９７８）の関係におけるａ及びｂなど）の１つ又は複数の値を求めることによって達成されてもよい。

【0027】

最適化されたフォワードモデルは、その後、インバージョンにおいて使用されることにより、より正確なインバージョンを生じさせることができる。

【0028】

ＲＨＰのある特定の値を所与として、地球物理学的パラメータの確率分布（及び／又は平均及び／又は分散値（直接的に））を計算するために、フォワードモデルがインバージョンにおいて使用されてもよい（下記の式（１３）を参照）。この確率分布関数は、地球物理学的パラメータの特定の複数の値を所与として、ＲＨＰの確率分布（及び／又は平均及び／又は分散値（直接的に））を計算するために使用されてもよい（下記の式（４）〜式（７）を参照）。

【0029】

地球物理学的構造における少なくとも２つの地球物理学的パラメータが提供されてもよい。この場合、本方法は、地球物理学的構造のＲＨＰを推定するために、少なくとも２つの地球物理学的パラメータをインバージョンさせることを含むことができる。

【0030】

少なくとも２つの地球物理学的パラメータを使用することは、それを行うことにより、ＲＨＰへの地球物理学的パラメータのインバージョンをかなり制限することができるので好ましい。１つのみの地球物理学的パラメータを用いてＲＨＰを推定することは、計算されたＲＨＰに大きな誤差及び不確実性を残す場合がある。しかしながら、同じＲＨＰを計算するために、同じインバージョンにおいて他の地球物理学的パラメータが使用されるとすぐに、不確実性が劇的に減少する。実際、より多くの地球物理学的パラメータが使用されるほど、計算されるＲＨＰはより正確となり得る。従って、少なくとも３つ、４つ又は５つの地球物理学的パラメータが、インバージョンに使用されてもよい。使用される１つのみ、２つ、３つ、４つ又は５つの地球物理学的パラメータが存在してもよい。

【0031】

インバージョンのステップは、少なくとも２つ（又は３つ、４つ、５つなど）の地球物理学的パラメータ間に統計的（条件付き）独立性が存在し、且つ各地球物理学的パラメータと地球物理学的構造のＲＨＰとの間に統計的依存性が存在する、モデルを使用することを含むことができる。

【0032】

本明細書における「モデル」とは、単に、１つ又は複数のフォワードモデルなどの、インバージョンに使用される数学的関係を意味するものであってよい。

【0033】

種々の地球物理学的パラメータに対するＲＨＰの統計的依存性と、種々の地球物理学的パラメータの互いに対する統計的（条件付き）独立性とは、本発明者らが発見した重要な概念である。インバージョン問題をこのようにモデリングすることによって、複数の地球物理学的パラメータの使用が、ＲＨＰに関して求められる値を制限し、従ってＲＨＰにおける誤差／不確実性を減少させるネットワークとして、地球物理学的パラメータ及びＲＨＰが見なされることを可能にする。

【0034】

この仮定及びこのようなモデルを用いて、地球物理学的パラメータとＲＨＰとの間の関係は、図１などの有向非巡回グラフ（ＤＡＧ）に示すことができるベイジアンネットワークの観点から表現されることが可能である。従って、本インバージョンは、インバージョン問題のベイジアン公式化であることができる。インバージョンは、ベイジアン統計設定において行われることができる。

【0035】

違う言い方をすれば、インバージョンに使用される数学的関係は、ＲＨＰの確率が、個別に各地球物理学的パラメータを条件とし、且つ地球物理学的パラメータ間に条件付き確率が存在しないという仮定に基づいて、選択されることができる。

【0036】

地球物理学的適用において、インバージョン問題にベイジアンアプローチを用いることが知られている。例えば、Ａｆｏｎｓｏ（２０１３ａ、２０１３ｂ）は、ベイジアンアプローチを用いることを教示する。しかしながら、ＲＨＰを計算するためにベイジアンアプローチを用いるという教示はない。例えば、Ａｆｏｎｓｏ（２０１３ａ、２０１３ｂ）は、実際に、マントルの温度及び組成を推定するためのインバージョン計算への入力としてＲＨＰを使用する。Ａｆｏｎｓｏ（２０１３ａ、２０１３ｂ）におけるＲＨＰの入力値は、単に、地球物理学的構造の組成を推定し、推定された組成をその組成に関するＲＨＰの期待公表値に関連付けることによって、求められる。従来の方法では、上記のような統計的ネットワークにおいて複数の地球物理学的パラメータを使用してＲＨＰを計算することはもちろん、地球物理学的パラメータからＲＨＰを計算するという教示又は提案は存在しない。

【0037】

少なくとも２つの地球物理学的パラメータは、少なくとも１つの電磁地球物理学的パラメータ（磁化率、電気伝導率又は抵抗率又は残留磁気など）と、少なくとも１つの力学的な地球物理学的パラメータ（密度又は地震波速度（好ましくは地震ｐ波速度）など）とを、含むことができる。好ましくは、重力及び磁気地球物理学的データ（それらから磁化率及び密度が計算されてもよい）が、通常、地球の広いエリアにわたって、３次元カバレージで入手可能であり／取得が簡単であるので、少なくとも磁化率及び密度が使用される。もちろん、磁化率、電気伝導率、抵抗率、残留磁気、密度又は好ましくは地震ｐ波速度である地震波速度（又はＲＨＰが依存するその他の地球物理学的パラメータ）の任意の組み合わせが使用されることができる。

【0038】

１つのみの地球物理学的パラメータがインバージョンに使用されてもよいケースと同様に、少なくとも２つの地球物理学的パラメータが使用される場合、インバージョンのステップは、各地球物理学的パラメータに対して、フォワードモデルを選択することを含んでもよい。フォワードモデルの各々は、それぞれの地球物理学的パラメータと地球物理学的構造のＲＨＰとの間の関係を定義する。

【0039】

複数のフォワードモデルは、関連の地球物理学的パラメータをＲＨＰに関連付ける期待される傾向に基づいて選択することができる。例えば、複数の地球物理学的パラメータは、通常、ＲＨＰの増加に伴って、（地球物理学的パラメータに応じて）増加又は減少し得る。複数の地球物理学的パラメータの１つが密度又は地震波速度（好ましくは、地震ｐ波速度）である場合、その地球物理学的パラメータは、ＲＨＰの増加に伴って減少し得る。複数の地球物理学的パラメータの１つが磁化率である場合、その地球物理学的パラメータは、ＲＨＰの増加に伴って増加し得る。全く同じフォワードモデルは本発明に必須ではなく、当業者は、１つ又は複数の地球物理学的パラメータとＲＨＰとの間の関係をモデリングするために使用することができる可能なモデルに気付くであろう。実際、種々の複数のフォワードモデルを、それらのフォワードモデルが複数の地球物理学的パラメータとＲＨＰとの間の一般傾向をモデリングすることができる限り、類似の結果を達成するために使用することができる。地球物理学的パラメータとＲＨＰとの間の関係の説明例として、Ｒｙｂａｃｈ（１９７８）は、ＲＨＰと地震ｐ波速度との間の指数関係

【数4】

を提案している。ただし、ＡはＲＨＰであり、ｖ_ｐは地震ｐ波速度であり、ａ及びｂは定数である。もちろん、これは単なる説明例であり、ＲＨＰと地震波速度との間の他の関係を使用することもできる。従って、上記から理解できるように、正確なフォワードモデルは、岩石物理学関係の知識に基づいて、当業者によって選択されることが可能である。

【0040】

また、上記で説明したように、ＲＨＰと各地球物理学的パラメータとの間の関係は、好ましくは、１つ又は複数のその他の地球物理学的パラメータなどのその他の変数に依存しない。もちろん、他の定係数が存在してもよいが、好ましくは、１つの変数のみが存在する。例えば、ＲＨＰと地震ｐ波速度との間の例示的なＲｙｂａｃｈ（１９７８）の関係に再度注目すると、ＲＨＰが依存する唯一の変数は、地震ｐ波速度である。Ｒｙｂａｃｈ（１９７８）の上記の式の他の係数（ａ及びｂ）は、単なる定数である。従って、例示的なＲｙｂａｃｈ（１９７８）の関係は、密度又は磁化率などの他の地球物理学的パラメータから統計的に（条件付きで）独立する。従って、それは、上記で説明したベイジアンネットワークでの使用に適している。別の言い方をすれば、そのような各フォワードモデルにおける唯一の変数は、地球物理学的構造のＲＨＰである。

【0041】

複数のモデルがシステムの完全な複雑さを示さない場合があること、すなわち、地球物理学的パラメータがＲＨＰのみに依存するようにモデルが意図的に単純化される場合があることを理解されたい。実際には、複数の地球物理学的パラメータは、一般に、多くの変数に依存する。しかしながら、本方法で使用される複数のモデルでは、複数の地球物理学的パラメータは、関心のある変数、この場合ＲＨＰにのみ依存することができる。

【0042】

地球物理学的構造のサンプルからの少なくとも２つの地球物理学的パラメータの各々の少なくとも１つの測定と地球物理学的構造のＲＨＰとを含む較正データが提供されてもよい。本方法は、較正データを取得することを含むことができる。較正は、好ましくは、地球物理学的構造のサンプルからの地球物理学的パラメータの各々の複数の測定と地球物理学的構造のＲＨＰとを含有することができる。

【0043】

インバージョンのステップは、較正データに基づいて、１つ又は複数の各フォワードモデルを最適化することを含むことができる。この最適化は、較正データを使用して、フォワードモデルの定係数の最適値を求めることを含むことができる。一般的に、較正データの量が多いほど、最適化がより良くなる。

【0044】

Ｒｙｂａｃｈ（１９７８）に提示された地球物理学的パラメータとＲＨＰとの間の例示的な関係に再度注目すると、較正データを用いて求めることができるのは、係数ａ及びｂである。

【0045】

複数のフォワードモデルを最適化するために、各提供された地球物理学的パラメータに対して、複数のフォワードモデル（これらは、所与のＲＨＰからそれぞれの地球物理学的パラメータを計算する）の各々が、ある特定の誤差分布を有する（すなわち、各提供された地球物理学的パラメータと、各フォワードモデルによって計算された各地球物理学的パラメータとの間の差が、誤差分布を与える）と仮定することができる。好ましくは、各フォワードモデルの誤差分布は、好ましくは平均が０のガウス誤差分布であると仮定される。これら複数のフォワードモデルは、例えば、可能な限り０に近くなる誤差分布の平均を有することによって、及び可能な限り誤差分布の小さな分散を有することによって、誤差分布が可能な限り小さくなるように誤差分布を減少させることによって、最適化されることができる。最適化は、複数のフォワードモデルを最適化するフォワードモデルの１つ又は複数の定係数（Ｒｙｂａｃｈ（１９７８）の関係におけるａ及びｂなど）の１つ又は複数の値を求めることによって達成することができる。

【0046】

１つ又は複数の最適化されたフォワードモデルは、その後、インバージョンにおいて使用されることにより、より正確なインバージョンを生じさせることができる。

【0047】

各地球物理学的パラメータに対して、種々の複数のフォワードモデルが使用されてもよい。

【0048】

ＲＨＰのある特定の値を所与として、複数の地球物理学的パラメータの各々の確率分布（及び／又は平均及び／又は分散値（直接的に））を計算するために、複数のフォワードモデルがインバージョンにおいて使用されることができる（下記の式（１３）を参照）。これらの確率分布関数は、複数の地球物理学的パラメータの特定の値を所与として、ＲＨＰの確率分布（及び／又は平均及び／又は分散値（直接的に））を計算するために、組み合わせることができる（下記の式（４）〜式（７）を参照）。

【0049】

地球物理学的構造の少なくとも１種類の地球物理学的データが提供されてもよく、この方法は、少なくとも１つの地球物理学的パラメータを計算するために、少なくとも１種類の地球物理学的データをインバージョンさせることを含む。同様に、地球物理学的構造の少なくとも２種類の地球物理学的データが提供されてもよく、この方法は、少なくとも２つ、３つ、４つ又は５つの地球物理学的パラメータを計算するために、少なくとも２、３、４又は５種類の地球物理学的データをインバージョンさせることを含む。

【0050】

地球物理学的パラメータが、密度、地震波速度（好ましくは、地震ｐ波速度）、磁化率、電気伝導率、電気抵抗率又は残留磁気である場合、地球物理学的データの種類は、それぞれ、重力データ、地震データ、磁気データ又はマグネトテルリックデータであることができる。データは、地震データ収集などの知られている手法を用いて収集することができる。本方法は、地球物理学的データを収集／取得することを含むことができる。

【0051】

少なくとも１つの地球物理学的パラメータを計算するために地球物理学的データをインバージョンさせることは、単一ドメインインバージョン又は結合インバージョン（これは、２次元又は３次元インバージョンでもよい）などの知られている手法を用いて行うことができる。例えば、重力及び磁気データが密度及び磁化率へとインバージョンされる場合、標準Ｇｒａｖｍａｇインバージョン手法（Ｇｅｏｓｏｆｔソフトウェアによって提供されるものなど）を使用することができる。当業者は、地球物理学的データを地球物理学的パラメータにインバージョンさせる多数のインバージョン方法を知っており、これらは本出願で説明される必要はない。

【0052】

上記で説明したように、少なくとも１つの地球物理学的パラメータは、ＲＨＰに依存するが、好ましくは、密度、地震波速度（好ましくは、地震ｐ波速度）、磁化率、電気伝導率、電気抵抗率若しくは残留磁気などの他の地球物理学的パラメータから独立した、又は他の地球物理学的パラメータから条件付きで独立するとしてモデリングされることが可能な任意の地球物理学的パラメータであることができる。パラメータとして表現されることが可能であり、且つＲＨＰに依存した（及び好ましくは、その他の地球物理学的パラメータ／変数から独立しながらＲＨＰに依存すると表現されることが可能な）地球物理学的構造のその他の特性を使用することができる。幾つものそのようなパラメータの任意の組み合わせを使用することができる。

【0053】

上記の方法は、地球物理学的構造のある特定の点／位置／体積／空間に関してＲＨＰを計算することができ、これらの点／位置／体積／空間は、インバージョンのステップで使用される地球物理学的パラメータの点／位置／体積／空間に対応する（これらの方法で使用される１つ又は複数の地球物理学的パラメータは、地球物理学的構造中のある点／位置／体積／空間におけるそのパラメータの値でもよい）ことを認識されたい。従って、空間依存のＲＨＰ関数Ａ（ｘ，ｙ，ｚ）を取得するために、上記のインバージョン方法は、地球物理学的構造中の各点／位置／体積／空間に関して点ごとに行うことができる。認識されるように、１つ又は複数の地球物理学的パラメータは、地球物理学的構造の空間上で変化する場合があり、これは空間的に変化するＲＨＰに対応する場合がある。

【0054】

従って、本方法は、空間依存のＲＨＰ関数Ａ（ｘ，ｙ，ｚ）を構築することを含むことができる。この関数は、地球物理学的構造中の各点／位置／体積／空間に関してＲＨＰを計算することによって、構築することができる。ＲＨＰは、地球物理学的構造の実質的に全体にわたって、又はある特定のエリア（ｘｙ）及び深さ（ｚ）にわたって、計算することができる。（当該分野で標準であるように、ｘ及びｙ軸は互いに直交する水平方向であり、ｚ軸は垂直方向である。）

【0055】

別の局面では、本発明は、地球物理学的構造における基礎熱流量を計算する方法であって、マントルからの熱流量寄与が提供され、地球物理学的構造の空間にわたって、上記で説明した方法のいずれかを用いて地球物理学的構造のＲＨＰを計算することと、地球物理学的構造の空間の少なくとも一部のＲＨＰを合計することと、ＲＨＰの合計値をマントルからの熱流量寄与に加算することと、を含む方法を提供する。

【0056】

地球物理学的構造の空間にわたるＲＨＰは、上記で説明した空間関数Ａ（ｘ，ｙ，ｚ）であることができる。

【0057】

ＲＨＰを合計することは、ある特定の深さ範囲にわたってＲＨＰを合計することを含むことができる。深さ範囲は、地殻の底部（すなわち、マントルの最上部の深さ）、又はリソスフェアの底部（すなわち、液状マントルの最上部の深さ）から、地殻又はリソスフェアの底部より上のある深さまでの深さであることができる。この地殻又はリソスフェアの底部より上の深さは、好ましくは、地殻又はリソスフェアの最上部であることができる。従って、ＲＨＰの合計は、地殻又はリソスフェア全体にわたるＲＨＰの合計であることができる。この積算は、２つの深さ間のｚ上で、Ａ（ｘ，ｙ，ｚ）を積分することによって達成することができる。

【0058】

ＲＨＰのこのような積算は、積算の上側深さ（例えば、地殻又はリソスフェアの最上部）における２次元熱流量分布を提供する。２次元熱流量分布は、ｘ及びｙに依存してもよい。

【0059】

本方法は、基礎熱流量に対するマントル寄与を取得することを含むことができる。マントル寄与は、マントルからの対流に起因する場合がある。マントル寄与は、知られている手法及び／又は知られているソフトウェアパッケージを用いて、計算することが可能である。

【0060】

基礎熱流量は、現在の基礎熱流量であってもよい。

【0061】

別の局面では、本発明は、地表上の地表熱流量を計算する方法であって、地表熱流量に対する堆積物寄与が提供され、上記の方法の何れかを用いて基礎熱流量を計算することであって、地球物理学的構造が地殻又はリソスフェアであり、基礎熱流量が地殻又はリソスフェアの最上部で計算されることと、堆積物寄与を基礎熱流量に加算することと、を含む方法を提供する。

【0062】

堆積物は、地表上、すなわち地殻／リソスフェアの最上部上に存在する場合又は存在しない場合がある層である。堆積物層が存在しない場合、基礎熱流量の計算方法は、堆積物寄与を考慮する必要なしに（あるいは、この場合、堆積物寄与が０であると見なすことができる）、地表熱流量を計算することができる。

【0063】

堆積物層は、堆積物中の放射性元素によって、熱を産生することができる。堆積物寄与は、一般的に、約１μＷｍ^−３であると考えられる。本方法は、堆積物寄与を計算することを含むことができる。ある地表位置における地表熱流量に対する寄与は、その地表位置における堆積層の深さを１μＷｍ^−３で乗算することによって求めることができる。

【0064】

地表熱流量は、好ましくは、（現在のデータが計算で使用される場合には）現在の地表熱流量でもよいが、（履歴データが使用される場合には）履歴熱流量でもよい。

【0065】

好ましくは定常状態近似において、地球物理学的構造の現在の温度分布及び最大古温度の推定を取得するために、地表熱流量を使用することができる。当業者は、この計算を行うための手法に気付くであろう。

【0066】

定常状態近似とは、熱状態（熱流量及び温度分布）が時間と共に変化しないという仮定である。定常状態近似では、温度は、フーリエの法則（１次元の場合では、ｑ＝ｋｄＴ／ｄｚ）によって与えられる。この場合、熱伝導率（ｋ）及び熱流量（ｑ）を所与として、深さの関数としての温度は、フーリエの法則の積分によって算出することができる。

【0067】

系が時間に依存する（すなわち、定常状態ではない）場合、一般には熱平衡から外れ、温度分布は、時間依存の拡散方程式（これは、フーリエの法則をエネルギー保存則と組み合わせて導出することができる）を解くことによって、算出することができる。

【0068】

本方法は、地球物理学的構造の熱履歴をモデリングすることを含むこともできる。熱履歴モデリングの場合、現在の熱流量は、運動学的復元及び熱流量履歴モデリングの照準点として使用することができる。

【0069】

別の局面では、本発明は、先行する請求項のいずれか一項に記載の方法を含む、地球物理学的構造の熱及び／又は温度モデルを生成する方法を提供する。

【0070】

認識されるように、上記の方法は、炭化水素の探査を行う際に、例えば、（見込みのある）掘削作業を計画して実行する際に使用される。本方法は、炭化水素の探査を行うために、計算されたＲＨＰ、地表熱流量、温度、又は熱／温度モデルを用いることをさらに含むことができる。

【0071】

別の局面では、本発明は、コンピュータ上で実行された際に、プロセッサに上記の方法のいずれかを行わせるように構成されたコンピュータ可読命令を含むコンピュータプログラムプロダクトを提供する。

【図面の簡単な説明】

【0072】

【図1】放射性崩壊による熱産生Ａと地球物理学的パラメータ｛χ，ρ，ｖｐ｝との間のモデル関係を表すベイジアンネットワーク（ＤＡＧ）である。

【図2】岩石物理学フォワードモデル及び較正データである。

【図3】測定ＲＨＰ（一番上）、密度（上から二番目）、磁化率（下から二番目）及び地震波速度（一番下）のデータサンプルである。

【図4】図３のデータに関する、本開示のベイジアン岩石物理学インバージョンの結果である。

【図5】図３のサンプル１に関する、本開示のベイジアン岩石物理学インバージョンから計算されたＲＨＰの確率密度である。

【図6】図３のサンプル７に関する、本開示のベイジアン岩石物理学インバージョンから計算されたＲＨＰの確率密度である。

【発明を実施するための形態】

【0073】

本発明の好適な実施形態が、単なる例として、添付の図面を参照して、これより説明される。

【0074】

図１は、放射性崩壊による熱産生Ａと地球物理学的パラメータ｛χ，ρ，ｖｐ｝との間のモデル関係を表すベイジアンネットワーク（ＤＡＧ）を示している。図１の右側に示されるように、地球物理学的パラメータは、拡張ベイジアンネットワークに含まれ得る地球物理学的データ（重力データ（ｍｚ）、磁気データ（ｇｚ）及び地震データ（Ｐ）のような）などの他のパラメータにも依存し得る。下記に提示される方法では、より単純な統計モデルのみが考慮される。さらに、図１は、ＲＨＰ（Ａ）が、Ｔｈ、Ｕ及びＫなどの放射性元素に依存することを示している。

【0075】

図２は、岩石物理学フォワードモデル及び較正データを示している。フォワードモデルは、地震波速度ｖ_ｐをＲＨＰ（Ａ）に関連付ける式（１０）の対数モデル（図２に示される線）であり、較正データ（図２に示されるデータ点）を用いて較正されている。

【0076】

図３は、測定ＲＨＰ（一番上）、密度（上から二番目）、磁化率（下から二番目）及び地震波速度（一番下）と共に、データサンプルを示している。これらのデータサンプルを用いて、本方法の有効性を示している。

【0077】

図４は、図３のデータに関する、本開示のベイジアン岩石物理学インバージョンの結果を示している。一番上のグラフにおいて、ＲＨＰは、密度のみを用いて計算される。中央のグラフにおいて、ＲＨＰは、密度及び磁化率のみを用いて計算される。一番下のグラフにおいて、ＲＨＰは、密度、磁化率及び地震波速度を用いて計算される。各サンプルの測定ＲＨＰ値は、菱形ドットで与えられる。また、各サンプルの計算ＲＨＰ値は、正方形ドットで与えられる。正方形ドットは、事後平均

【数5】

を示している。また、エラーバーは、

【数6】

によって与えられる。

【0078】

図５は、図３のサンプル１に関する、本開示のベイジアン岩石物理学インバージョンから計算されたＲＨＰの確率密度を示している。一番上のグラフには、事前分布が示されている。中央のグラフには、尤度分布が示されている。一番下のグラフには、事後分布が示されている。中央及び一番下のグラフにおいて、分布は、密度のみ（鎖線）と、密度及び磁化率（実線）と、密度、磁化率及び地震波速度（実線）とのインバージョンによって計算される（サンプル１の場合、密度及び磁化率の分布と、密度、磁化率及び地震波速度の分布とは、互いに実質的に重なるので、合計２本の線のみを見ることができる）。

【0079】

図６は、図３のサンプル７に関する、本開示のベイジアン岩石物理学インバージョンから計算されたＲＨＰの確率密度である。一番上のグラフには、事前分布が示されている。中央のグラフには、尤度分布が示されている。一番下のグラフには、事後分布が示されている。中央及び一番下のグラフでは、分布は、密度のみ（点線）と、密度及び磁化率（鎖線）と、密度、磁化率及び地震波速度（実線）とのインバージョンによって計算される。

【0080】

下記の例示的方法に示されるように、基礎熱流量及び地表熱流量の推定は、好ましくは、少なくとも２つの地球物理学的パラメータのインバージョンを用いて取り扱われる。本方法の一般的ワークフローとして、第１に、地球物理学的パラメータが、単一ドメイン又は結合インバージョン（例えば、重力及び磁気学、任意選択的に地震データのもの）によって取得される。これは２次元又は３次元のインバージョンでもよい。第２に、地殻岩石物理学関係を用いて、放射性崩壊による熱産生（ＲＨＰ）を取得するために、地球物理学的パラメータ（密度、磁化率及び任意選択的に地震波速度など）がインバージョンされる。現在の熱流量は、ＲＨＰを熱流量に対するマントル寄与と合わせることによって取得される。熱流量のマントル部分を算出するために、Ａｆｏｎｓｏら（２００８）及びＦｕｌｌｅａら（２００９）によって開発された、ＬｉｔＭｏｄ３Ｄソフトウェアを使用することができる。推定された基礎熱流量及び地表熱流量は、定常状態近似の範囲内で温度を算出するために、直接適用することができる。熱履歴モデリングの場合、現在の熱流量は、運動学的復元及び熱流量履歴モデリングの照準点として使用することができる。このような運動学的復元を達成するための手法が、当該分野で知られている（ＭｃＫｅｎｚｉｅ、１９７８）。

【0081】

従って、本開示では、ＲＨＰの特徴を利用して、地球物理学的データ及び地球物理学的パラメータのインバージョンから、現在の地殻ＲＨＰの推定値を取得するワークフローが示される。下記の例では、３種類の地球物理学的データである磁気データ、重力データ、及び地震データが考慮される。地域研究にとって最も重要なのは、通常、広いエリアにわたって３次元カバレージで入手可能な重力及び磁気データである。マグネトテルリック（ＭＴ）データ、残留磁気及び抵抗率が含まれてもよい。地球物理学的データから、磁化率、密度及び地震ｐ波速度が、インバージョンによって取得される。地震波の場合、インバージョンという用語は、トモグラフィ又は全波形インバージョンとして理解されるものとする。磁気及び重力データの場合、インバージョンは、ジオメトリで制限されるＧｒａｖｍａｇインバージョンである。

【0082】

下記の例示的方法では、放射性崩壊による熱産生（ＲＨＰ）、現在の基礎熱流量及び地表熱流量が、地球物理学的インバージョンから算出される。

【0083】

ＲＨＰは、地殻岩石中のＴｈ、Ｕ及びＫの長寿命同位体の発生によるものである。本方法では、地球物理学的データとＲＨＰとの間の関係が、図１に示されるベイジアンネットワーク又は有向非巡回グラフ（ＤＡＧ）の観点からモデリングされることが可能であると仮定される。

【0084】

ベイジアンネットワークの同時確率分布は、親ノードの周辺分布及び子の条件付き分布によって定義される。変数（ｘ_ｉ，…，ｘ_ｎ）の同時分布は、

【数7】

によって与えられる。ただし、

【数8】

は、親ノードを示す。ネットワークのトップノードは、親を持たない。一般ベイジアン因数分解ルールである式（１）を図１のＤＡＧに適用すると、ＲＨＰ及び地球物理学的パラメータの同時確率は、

【数9】

として記述することができる。ただし、ｍ＝｛χ，ρ，ｖ_ｐ｝であり、χは磁化率であり、ρは質量密度であり、ｖ_ｐはＰ波速度であり、ＡはＲＨＰである。ＲＨＰの次元は、μＷ／ｍ^３である。事前分布ｐ（Ａ；λ）は、後に説明されるハイパーパラメータλに依存する。ｍには、電気抵抗率が含まれていてもよい。パラメータｍ_ｉの条件付き独立性を用いて、式（２）の同時分布は、

【数10】

として記述することもできる。また、式（２）及び式（３）から、ＲＨＰの事後分布を取得することができ、

【数11】

ｍ_ｉの代わりに実際のパラメータを用いることによって、以下の式：

【数12】

が与えられる。

【0085】

Ａの事後分布が既知である場合、事後期待値及び事後分散は、

【数13】

によって与えられる。

【0086】

所与の地球物理学的パラメータの予想されるＲＨＰの計算に最も有用なのは、式（４）〜式（７）である。しかしながら、式（４）及び式（５）から明らかなように、これを行うためには、地球物理学的パラメータｍ_ｉの各々の尤度関数ｐ（ｍ_ｉ｜Ａ）を知る必要がある。さらに、事前分布ｐ（Ａ；λ）を知る必要がある。これらの計算方法は、下記で与えられる。

【0087】

事前分布ｐ（Ａ；λ）に関して、本方法では、ガウス分布であると仮定される。

【数14】

ただし、μ_Ａ及び

【数15】

は、それぞれ事前期待値及び分散である。事前分布は、例えば、ＲＨＰが通常比較的狭い範囲０＜Ａ＜１０μＷ／ｍ^３内にあるというＲＨＰに関するユーザの事前知識を組み込む。ハイパーパラメータλは、地質学的又は岩石学的設定に関するユーザの事前知識を反映する。事前知識として、例えば、ユーザは、平均ＲＨＰが、常に海洋地殻よりも大陸地殻において高く、一般的に珪長質岩の場合μ_Ａ〜２μＷ／ｍ^３であり、苦鉄質岩の場合μ_Ａ≦１μＷ／ｍ^３であることを知っていることとする。ユーザの事前知識が疎である場合、事前分散σ_Ａは、それに応じて大きくなるはずである。

【0088】

従って、事前分布は、好ましくは統計的分布、好ましくはガウス分布でもよい。好ましくは、事前分布の平均及び分散は、問題の地球物理学的構造のユーザの事前知識（例えば、地球物理学的構造が海洋又は大陸地殻であるか、及びＲＨＰの標準分散の知識）に基づいて、ユーザによって選択される。

【0089】

式（４）及び式（５）の右側の尤度関数ｐ（ｍ_ｉ｜Ａ）に関して、これらは、フォワードモデルＦ_ｉ（Ａ）を用いて計算される。フォワードモデルは、所与のＲＨＰ（Ａ）に対して、関連の地球物理学的パラメータｍ_ｉを算出する数学的関係である。本方法に関して、フォワードモデルＦ_ｉ（Ａ）の各々が、それぞれの測定／観測地球物理学的パラメータと比較された際に、平均が０のそれぞれのガウス誤差分布、すなわち、

【数16】

を有すると仮定される。ただし、σ_ｅｉは誤差分散である。

【0090】

上記に詳細に説明されたように、使用される特定のフォワードモデルは、本発明に必須ではなく、当業者は、使用すべき適切なフォワードモデルに気付くであろう。例えば、地球物理学的パラメータが地震Ｐ波速度である場合、フォワードモデルは、Ｒｙｂａｃｈ（１９７８）に開示されたような、ＲＨＰに対する対数依存性を有するものであってもよい。

【数17】

【0091】

他の地球物理学的パラメータに対して、対応するフォワード関数を求めることができ、各フォワード関数は、ＲＨＰ（Ａ）をそれぞれの地球物理学的パラメータに関連付ける。例えば、密度、地震波速度及び／又は伝導率に関して、

【数18】

又は

【数19】

などの任意の崩壊関数が使用されてもよい。磁化率に関しては、任意のシグモイド関数が使用されてもよい。

【0092】

式（１０）から認識されるように、フォワード関数は、定係数（式（１０）におけるａ及びｂなど）によって制御されてもよい。本方法では、これらの定係数は、地球物理学的構造から採取され得る岩石サンプルに対して測定され得る較正データから決定されてもよい。Ｎ個の岩石サンプル測定

【数20】

のセット（ただし、ｊは地球物理学的パラメータｍ_ｉ及びＲＨＰ（Ａ）のサンプル番号を示す）、並びに、フォワードモデルＦ_ｉ（Ａ）が与えられると、式（９）のガウス誤差の期待値及び分散の標準推定量は、それぞれ、

【数21】

として求めることができる。

【0093】

【数22】

が最小であり、且つμ_ｅｉ≒０であるような、最適パラメータ（式（１０）におけるｖ_∞、ｖ_０及びβなど）。これらは、回帰などの最適化を目的とした任意の知られている数学的手法を用いて求めることができる。図２は、較正データ（データ点）を用いた較正後のＰ波速度及びＲＨＰの対数フォワードモデル（線）の一例を示している。

【0094】

各フォワードモデルの最適パラメータが較正によって求められると、最大尤度関数（すなわち、ＲＨＰの値を所与とした地球物理学的パラメータの最大尤度である、式（４）及び式（５）におけるｐ（ｍ_ｉ｜Ａ））が、

【数23】

によって明確に与えられる。

【0095】

これらの最大尤度関数は、上記で説明された事前分布と共に、式（４）〜式（７）において使用されて、測定／取得地球物理学的パラメータを所与としたＲＨＰの事後分布

【数24】

と、事後期待値及び事後分散を計算する。このようにして、ＲＨＰが計算される。すなわち、これらの数量は、基礎熱流量、地表熱流量及び温度分布の計算に使用することができるＲＨＰの有用な値を提供する。

【0096】

インバージョン問題のベイジアン公式化を用いて、本方法は、提案された地殻岩石物理学モデル（フォワードモデル）が、完全には観測（較正データ）を記述しないという事実を尊重することを認識されたい。この不完全さは、式（１３）で与えられる尤度分布における誤差分散σ_ｅｉによって説明される。これは次に、岩石物理学インバージョンによって取得されるＲＨＰの事後分散

【数25】

の量的推定を提供する。従って、式（３）は、効果的に、事後平均及び事後分散が式（６）及び式（７）によって与えられる、Ａの一変量分布である。

【0097】

上記のステップは、単に、地球物理学的構造中のある特定の点に関して、ＲＨＰ（Ａ）を計算した場合があることも認識されたい。この点は、ＲＨＰの計算に使用された各地球物理学的パラメータの値の位置に対応した位置である。従って、好ましくは、上記のステップで使用される全ての地球物理学的パラメータは、地球物理学的構造中の同じ、又は少なくとも類似した位置から得られる。

【0098】

地球物理学的構造のある領域又は全体にわたって、ＲＨＰの概観を構築するために、Ａを計算するための上記のステップは、地球物理学的構造中の異なる位置に対して実行される。しかしながら、フォワードモデルの較正は、一度のみ実行されてもよい。すなわち、それは異なるそれぞれの位置に対して実行される必要はない。状況によっては、較正は各位置に対して実行されてもよい。

【0099】

従って、上記で説明した岩石物理学インバージョンは、上に述べたように点ごとに適用されて、空間的に変化するＲＨＰのＡ（ｘ，ｙ，ｚ）を取得する。

【0100】

Ａ（ｘ，ｙ，ｚ）が求められると、基礎熱流量を求めることができる。上述の通り、基礎熱流量（これは、地殻の最上部の熱流量でもよい）は、２つの主な部分である、地殻におけるＲＨＰからの寄与と、マントル対流からの寄与とから成る。地殻部分は、地殻の深さｚ（例えば、地殻の最上部Ｚ_Ｔから地殻の底部Ｚ_Ｂまで）にわたって、Ａ（ｘ，ｙ，ｚ）の積分によって近似することができる。基礎熱流量は、基礎熱流量を取得するために地殻（ＲＨＰ）及びマントル（対流）からの寄与を加算することによって求めることができる。従って、基礎熱流量は、

【数26】

によって与えられる。ただし、ｑ_Ｍはマントルからの寄与である。マントルからの寄与は、ＬｉｔＭｏｄ３Ｄソフトウェアを用いるなど、知られている手法を用いて求めることができる。

【0101】

さらに、基礎熱流量が求められると、地表熱流量ｑ_０の近似は、堆積物からの地表熱流量に対する寄与を基礎熱流量に加算することによって取得することができる。例えば、地表熱流量は、平均堆積物熱産生を基礎熱流量に加算することによって求めることができる。

【数27】

ただし、Ｚ_Ｔ−Ｚ_Ｓは、堆積パッケージの厚さであり、Ａ_Ｓ〜１μＷ／ｍ^２である。Ｚ_Ｔは、位置に応じて表面から０〜約２０ｋｍの深さでもよい。Ｚ_Ｓは、陸上では０であってもよく、又は海上では海底の深さであってもよい、堆積層の上側部分の深さである。Ｚ_Ｂは、モホ面まで又は下部地殻の最上部（地殻の最上部とモホ面との間の約２／３であってもよい）までの深さであってもよい。

【0102】

またさらに、地球物理学的構造の（例えば、地殻／リソスフェア、及び任意選択的に堆積物層の）現在の温度及び最大古温度は、地表熱流量から、定常状態近似を用いて求めることができる。

【0103】

従って、上記のステップは、本発明による、地球物理学的構造のＲＨＰ、基礎熱流量、地表熱流量及び温度を計算するための例示的ワークフローを提示する。

【0104】

以下は、上記で説明された方法が測定データに対して実施された数値例である。この例は、本方法の有効性、及び具体的には複数の地球物理学的パラメータの使用がどのようにＲＨＰを求めるインバージョン計算を向上させるかを示す。

【0105】

図３に関して、これは、ＲＨＰ（Ａ）及び幾つかの地球物理学的パラメータ（密度、磁化率、及び地震波速度）の測定値を示している。（ちなみに、サンプルのこのような測定値は、上述のようにフォワードモデルを較正するために使用することもできる）。地球物理学的パラメータの１１個のサンプルが示されている。

【0106】

数値テストの目的は、本インバージョン方法を用いることによって、ＲＨＰが、地球物理学的パラメータから正確に計算されることが可能であることを示すことである。

【0107】

図４は、本推定方法の３つの異なるテストランの結果を示している。第１のテストラン（一番上のグラフ）では、１つのみの地球物理学的パラメータが使用される（密度のみ）。第２のテストラン（中央のグラフ）では、２つの地球物理学的パラメータが、ベイジアン設定で使用される（密度及び磁化率のみ）。第３のテストラン（一番下のグラフ）では、３つの地球物理学的パラメータが、ベイジアン設定で使用される（密度、磁化率、及び地震波速度のみ）。正方形ドットは、事後平均

【数28】

を示しており、エラーバーは、

【数29】

によって与えられる。

【0108】

このように、インバージョンにおいて１つのみの地球物理学的パラメータを用いることにより、一般傾向が正しく、より高い及びより低いＲＨＰを有するサンプルを区別し、並びに、測定ＲＨＰは、インバージョンのエラーバーの範囲内である。しかしながら、エラーバーはかなり大きく、また、平均値の一部は、ＲＨＰサンプルの実際の値（菱形ドット）からかなり離れている。定量的には、インバージョン結果は、大きなＲＨＰを有するサンプルの場合、測定値からかけ離れている。従って、密度のみのインバージョンは、非常に低いＲＨＰ（高密度による）を検出することはできるが、高ＲＨＰを有するサンプルの定量値は、あまり上手くとらえることができない。

【0109】

しかしながら、下側の２つのグラフから分かるように、２つ、又は好ましくは３つの地球物理学的パラメータを使用することによって、全体として、ＲＨＰの平均及び分散のサイズの精度を劇的に向上させることができる。この向上は、複数の地球物理学的パラメータを用いたＲＨＰ分布のより良い制限により生じる。

【0110】

この向上は、図５及び図６にさらに図示されている。図５及び図６は、図３から得たサンプル１（図５）及びサンプル７（図６）からのデータを用いて、本インバージョン方法によって計算されたＲＨＰの事前（一番上のグラフ）、尤度（中央のグラフ）及び事後（一番下のグラフ）分布を示している。これら２つのサンプルは、我々が区別したい２つの主要なケースである、比較的高いＲＨＰ（サンプル１、図５）と、低ＲＨＰ（サンプル７、図６）との典型である。

【0111】

図５及び図６の両方の一番上のグラフから分かるように、事前分布は、広く（無情報）、また反転をあまり誘導しない。

【0112】

図５の中央のグラフから分かるように、１つのみの地球物理学的パラメータ（密度）が使用された場合、ＲＨＰの尤度分布は、最初は非常に広い（鎖線）。しかしながら、２つ又は３つの地球物理学的パラメータ（密度及び磁化率、並びに密度、磁化率及び地震波速度）がインバージョンにおいて使用される場合、分散が減少し、且つ平均がＲＨＰの測定値（サンプル１の場合、これは、Ａ＝２．８０μｍＷ／ｍ^３であった）を中心とする。図５の中央のグラフに関して、これらの分布の両方は、それらが互いに密接に重なるので、実線で示される。

【0113】

同様に、図６の中央のグラフから分かるように、１つのみの地球物理学的パラメータ（密度）が使用された場合、ＲＨＰの尤度分布は、最初は非常に広い（点線）。しかしながら、２つの地球物理学的パラメータ（密度及び磁化率）がインバージョンにおいて使用される場合、分散が減少し、且つ平均がＲＨＰの測定値（サンプル７の場合、これは、Ａ＝０．３７μｍＷ／ｍ^３であった）を中心とする。これは鎖線で示される。さらに、３つの地球物理学的パラメータ（密度、磁化率及び地震波速度）がインバージョンにおいて使用される場合、分散がさらに減少し、且つ平均がさらにＲＨＰの測定値（この場合、これは、Ａ＝０．３７μｍＷ／ｍ^３であった）を中心とする。これは実線で示される。

【0114】

尤度分布と同様に、図５の一番下のグラフから分かるように、１つのみの地球物理学的パラメータ（密度）が使用された場合、ＲＨＰの事後分布は、最初は広い（鎖線）。しかしながら、２つ又は３つの地球物理学的パラメータ（密度及び磁化率、並びに密度、磁化率及び地震波速度）がインバージョンにおいて使用される場合、分散が減少し、且つ平均がＲＨＰの測定値（サンプル１の場合、これは、Ａ＝２．８０μｍＷ／ｍ^３であった）を中心とする。図５の中央のグラフに関して、これらの分布の両方は、それらが互いに密接に重なるので、実線で示される。

【0115】

同様に、図６の中央のグラフから分かるように、１つのみの地球物理学的パラメータ（密度）が使用された場合、ＲＨＰの事後分布は、最初は広い（点線）。しかしながら、２つの地球物理学的パラメータ（密度及び磁化率）がインバージョンにおいて使用される場合、分散が減少し、且つ平均がＲＨＰの測定値（サンプル７の場合、これは、Ａ＝０．３７μｍＷ／ｍ^３であった）を中心とする。これは鎖線で示される。さらに、３つの地球物理学的パラメータ（密度、磁化率及び地震波速度）がインバージョンにおいて使用される場合、分散がさらに減少し、且つ平均がさらにＲＨＰの測定値（この場合、これは、Ａ＝０．３７μｍＷ／ｍ^３であった）を中心とする。これは実線で示される。

【0116】

従って、図３〜図６は、地球物理学的パラメータからＲＨＰを計算する本インバージョン方法が、たった１つの地球物理学的パラメータに対しても効果的に機能することを実証する。しかしながら、それらは、ＲＨＰの確率分布を正しい値に制限する際の、インバージョンにおける複数の地球物理学的パラメータの使用の重要性も明白に実証する。

【0117】

数学的関係、ステップ及び手法が提示される本明細書では、これらは、これらの関係、ステップ及び手法に対する、どのような些細な変更も対象に含めると見なされるべきであることを認識されたい。当業者は、どのような小さな／些細な／形式上の変更（すなわち、例えば単純に式を再配置する、式を組み合わせる、又はステップを些細に並べ替えるといった、本発明に使用される数学的プロセス全体を変えない変更）も、本明細書に提示されたのと同じ関係、ステップ及び手法を使用していると気付くであろう。

【0118】

参照文献
Ａｆｏｎｓｏ，Ｊ．，Ｈ．Ｆｅｒｎａｎｄｅｚ，Ｇ．Ｒａｎａｌｌｉ，Ｗ．Ｇｒｉｆｆｉｎ，ａｎｄ Ｊ．Ｃｏｎｎｏｌｌｙ，２００８，Ｉｎｔｅｇｒａｔｅｄ ｇｅｏｐｈｙｓｉｃａｌ−ｐｅｔｒｏｌｏｇｉｃａｌ ｍｏｄｅｌｌｉｎｇ ｏｆ ｔｈｅ ｌｉｔｈｏｓｐｈｅｒｅ ａｎｄ ｓｕｂｌｉｔｈｏｓｐｈｅｒｅ ｕｐｐｅｒ ｍａｎｔｌｅ：Ｍｅｔｈｏｄｏｌｏｇｙ ａｎｄ ａｐｐｌｉｃａｔｉｏｎｓ：Ｇ３，９．
Ａｆｏｎｓｏ，Ｊ．，Ｆｕｌｌｅａ，Ｊ．，Ｇｒｉｆｆｉｎ，Ｗ．，Ｙａｎｇ，Ｙ．，Ｊｏｎｅｓ，Ａ．，Ｃｏｎｎｏｌｌｙ，Ｊ．ａｎｄ Ｏ’Ｒｅｉｌｌｙ，Ｓ．，２０１３ａ，３−Ｄ ｍｕｌｔｉｏｂｓｅｒｖａｂｌｅ ｐｒｏｂａｂｉｌｉｓｔｉｃ ｉｎｖｅｒｓｉｏｎ ｆｏｒ ｔｈｅ ｃｏｍｐｏｓｉｔｉｏｎａｌ ａｎｄ ｔｈｅｒｍａｌ ｓｔｒｕｃｔｕｒｅ ｏｆ ｔｈｅ ｌｉｔｈｏｓｐｈｅｒｅ ａｎｄ ｕｐｐｅｒ ｍａｎｔｌｅ．Ｉ：ａ ｐｒｉｏｒｉ ｐｅｔｒｏｌｏｇｉｃａｌ ｉｎｆｏｒｍａｔｉｏｎ ａｎｄ ｇｅｏｐｈｙｓｉｃａｌ ｏｂｓｅｒｖａｂｌｅｓ，Ｊｏｕｒｎａｌ ｏｆ Ｇｅｏｐｈｙｓｉｃａｌ Ｒｅｓｅａｒｃｈ：Ｓｏｌｉｄ Ｅａｒｔｈ，Ｖｏｌ，１１８，２５８６−２６１７．
Ａｆｏｎｓｏ，Ｊ．，Ｆｕｌｌｅａ，Ｊ．，Ｙａｎｇ，Ｙ．，Ｊｏｎｅｓ，Ａ．，ａｎｄ Ｃｏｎｎｏｌｌｙ，２０１３ｂ，３−Ｄ ｍｕｌｔｉ−ｏｂｓｅｒｖａｂｌｅ ｐｒｏｂａｂｉｌｉｓｔｉｃ ｉｎｖｅｒｓｉｏｎ ｆｏｒ ｔｈｅ ｃｏｍｐｏｓｉｔｉｏｎａｌ ａｎｄ ｔｈｅｒｍａｌ ｓｔｒｕｃｔｕｒｅ ｏｆ ｔｈｅ ｌｉｔｈｏｓｐｈｅｒｅ ａｎｄ ｕｐｐｅｒ ｍａｎｔｌｅ．ＩＩ：Ｇｅｎｅｒａｌ ｍｅｔｈｏｄｏｌｏｇｙ ａｎｄ ｒｅｓｏｌｕｔｉｏｎ ａｎａｌｙｓｉｓ，Ｊｏｕｒｎａｌ ｏｆ Ｇｅｏｐｈｙｓｉｃａｌ Ｒｅｓｅａｒｃｈ：Ｓｏｌｉｄ Ｅａｒｔｈ，Ｖｏｌ，１１８，１６５０−１６７６．
Ａｌｌｅｎ，Ｐ．，ａｎｄ Ｊ．Ａｌｌｅｎ，２００５，Ｂａｓｉｎ ａｎａｌｙｓｉｓ：Ｂｌａｃｋｗｅｌｌ Ｐｕｂｌｉｓｈｉｎｇ．
Ｆｕｌｌｅａ，Ｊ．，Ｊ．Ａｆｏｎｓｏ，Ｊ．Ｃｏｎｎｏｌｌｙ，Ｈ．Ｆｅｒｎａｎｄｅｚ，Ｄ．Ｇａｒｃｉａ−Ｃａｓｔｅｌｌａｎｏｓ，ａｎｄ Ｈ．Ｚｅｙｅｎ，２００９，Ｌｉｔ−Ｍｏｄ３Ｄ：Ａｎ ｉｎｔｅｒａｃｔｉｖｅ ３−Ｄ ｓｏｆｔｗａｒｅ ｔｏ ｍｏｄｅｌ ｔｈｅ ｔｈｅｒｍａｌ，ｃｏｍｐｏｓｉｔｉｏｎａｌ，ｄｅｎｓｉｔｙ，ｓｅｉｓｍｏｌｏｇｉｃａｌ，ａｎｄ ｒｈｅｏｌｏｇｉｃａｌ ｓｔｒｕｃｔｕｒｅ ｏｆ ｔｈｅ ｌｉｔｈｏｓｐｈｅｒｅ ａｎｄ ｓｕｂｌｉｔｈｏｓｐｈｅｒｉｃ ｕｐｐｅｒ ｍａｎｔｌｅ：Ｇ３，１０．
ＭｃＫｅｎｚｉｅ，Ｄ．Ｐ．，１９７８，Ｓｏｍｅ ｒｅｍａｒｋｓ ｏｎ ｔｈｅ ｄｅｖｅｌｏｐｍｅｎｔ ｏｆ ｓｅｄｉｍｅｎｔａｒｙ ｂａｓｉｎｓ：Ｅａｒｔｈ ａｎｄ Ｐｌａｎｅｔａｒｙ Ｓｃｉｅｎｃｅ Ｌｅｔｔｅｒｓ，４０，２５−３２．
Ｒｙｂａｃｈ，Ｌ．，１９７８，Ｔｈｅ ｒｅｌａｔｉｏｎｓｈｉｐ ｂｅｔｗｅｅｎ ｓｅｉｓｍｉｃ ｖｅｌｏｃｉｔｙ ａｎｄ ｒａｄｉｏａｃｔｉｖｅ ｈｅａｔ ｐｒｏｄｕｃｔｉｏｎ ｉｎ ｃｒｕｓｔａｌ ｒｏｃｋｓ：Ａｎ ｅｘｐｏｎｅｎｔｉａｌ ｌａｗ：Ｐｕｒｅ Ａｐｐｌ．Ｇｅｏｐｈｙｓ，１１７，７５−８２．

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】