(57)【要約】  本発明は、ある領域の地下の特性分析を行う方法であって、前記領域内の受動的地震信号が記録される少なくとも１つの地点において前記領域の地表付近で取得される受動的地震信号のスペクトル密度を表す複数のスペクトルであって、それぞれが運動を表す信号から作成される複数のスペクトルを作成するステップ（Ｅ０１）と、記録および周波数に関連付けられたスペクトル属性のセットを取得するため、前記スペクトルそれぞれに現れる周波数ごとに少なくとも１つのスペクトル属性を決定するステップ（Ｅ０２）と、前記スペクトル属性のセットを、各行が１つの記録に関連付けられた行列に編成するステップ（Ｅ０３）と、主成分を決定して前記主成分から前記地下の特性を推定するために、主成分分析法を前記行列に適用するステップ（Ｅ０４）と、を含む方法を提供する。【選択図】図１
2019519799.tif 000002
【特許請求の範囲】
【請求項１】
ある領域の地下の特性分析を行う方法であって、
前記領域内の受動的地震信号が記録される少なくとも１つの地点において前記領域の地表付近で取得される受動的地震信号のスペクトル密度を表す複数のスペクトルであって、それぞれが運動を表す信号から作成される複数のスペクトルを作成するステップ（Ｅ０１）と、
記録および周波数に関連付けられたスペクトル属性のセットを取得するため、前記スペクトルそれぞれに現れる周波数ごとに少なくとも１つのスペクトル属性を決定するステップ（Ｅ０２）と、
前記スペクトル属性のセットを、各行が１つの記録に関連付けられた行列に編成するステップ（Ｅ０３）と、
主成分を決定して前記主成分から前記地下の特性を推定するために、主成分分析法を前記行列に適用するステップ（Ｅ０４）と、
を含む方法。
【請求項２】
前記運動は鉛直運動および／または水平運動であり、
周波数ごとの前記スペクトル属性は、鉛直地震動のスペクトル密度と水平地震動のスペクトル密度の比、前記水平地震動の周波数に関するスペクトル密度の導関数、および、前記鉛直地震動の周波数に関するスペクトル密度の導関数からなる群より選択される種類のものである、請求項１に記載の方法。
【請求項３】
前記水平地震動の周波数に関する前記スペクトル密度の前記導関数、および／または、前記鉛直地震動の周波数に関する前記スペクトル密度の前記導関数は、選択した数のスペクトルの点に関して線形回帰を適用することで算出される、請求項２に記載の方法。
【請求項４】
信号から各スペクトルを作成することは、
前記信号（ＳＩＧ）を、期間がすべて同じである複数の連続するサブ信号に分割することと、
前記サブ信号ごとにスペクトル密度のサブスペクトル（ＰＳＤ＿１、ＰＳＤ＿２、ＰＳＤ＿３）を作成することと、
前記サブスペクトルの各周波数について、前記サブスペクトルそれぞれにおける前記周波数についてのスペクトル密度の値から、スペクトル密度の統計学的属性を求めることと、
すべての前記周波数のすべての前記統計学的属性から作成される前記スペクトル（ＰＳＤ＿ｍ）を取得することと、を含む、請求項１〜３のいずれか１項に記載の方法。
【請求項５】
前記サブ信号はそれぞれ、先行するサブ信号の０を超える期間にわたり前記先行するサブ信号に重なる、請求項４に記載の方法。
【請求項６】
前記記録は前記領域の異なる地点で作られ、
前記スペクトル属性のセットの各属性は地点にも関連付けられている、請求項１〜５のいずれか１項に記載の方法。
【請求項７】
前記記録は所定の時点から所定の期間にわたり作られる、請求項１〜６のいずれか１項に記載の方法。
【請求項８】
前記記録は、同時に作られた記録ごとのグループとして作られ、
前記グループはそれぞれ、前記グループの前記記録が作られた日にちに対応し、
前記記録は前記領域の異なる地点で、および／または、異なる時点から作られる、請求項７に記載の方法。
【請求項９】
前記行列における所定の属性（Ａｔｔｒ＿ｉ）に関連付けられた列は、すべて隣接している、請求項１〜８のいずれか１項に記載の方法。
【請求項１０】
前記記録のグループはそれぞれ、前記行列の複数行からなる１つのグループに関連付けられ、
前記属性の値は、複数行からなるグループごとに正規化される、請求項８に記載の方法。
【請求項１１】
前記主成分は射影作用素であり、
前記行列を前記射影作用素それぞれに射影して、記録ごとに前記行列を射影した結果を示す前記領域のグラフィックによる表現（ＰＲＪ１、ＰＲＪ２、ＰＲＪ３、ＰＲＪ４）を、射影作用素ごとに取得する、請求項１〜１０のいずれか１項に記載の方法。
【請求項１２】
Ｋ個の射影作用素を前記射影作用素の中から決定する、請求項１１に記載の方法。
【請求項１３】
長さＫのベクトルを前記行列の行ごとに得るため、前記行列を前記Ｋ個の前記射影作用素に射影することと、
長さＫの前記ベクトルから第２の行列を得ることと、
Ｎ個の、複数行のクラスを得るため、一次元または二次元で編成された分類法を前記第２の行列に適用することと、
前記領域の前記地下のアノマリーの規模を表す前記行列の各行に少なくとも１つの値を割り当てることと、
複数行のクラスごとにクラスのヘッドを作成することと、
Ｎ×Ｋ次元の第３の行列を前記クラスのヘッドから得ることと、
スペクトル属性の第１の行列の各行に現れる周波数の数のＮ倍の次元の第４の行列を得るため、前記第３の行列に擬似逆行列法を適用することと、
をさらに含む、請求項１２に記載の方法。
【請求項１４】
長さＫのベクトルを前記行列の行ごとに得るため、前記行列を前記Ｋ個の射影作用素に射影することと、
長さＫの前記ベクトルから第２の行列を得ることと、
スペクトル属性の第１の行列と同じ次元を有する第３の行列を得るため、前記第２の行列に擬似逆行列法を適用することと、
Ｎ個の、複数行のクラスを得るため、一次元または二次元で編成された分類法を前記第３の行列に適用することと、
前記領域の前記地下のアノマリーの規模を表す前記行列の各行に少なくとも１つのクラス番号を割り当てることと、
複数行のクラスごとにクラスのヘッドを作成することと、
スペクトル属性の前記第１の行列の各行に現れる周波数の数のＮ倍の次元の第４の行列を前記クラスのヘッドから得ることと、
をさらに含む、請求項１２に記載の方法。
【請求項１５】
ある領域の地下の特性分析を行うシステムであって、
前記領域内の受動的地震信号が記録される少なくとも１つの地点において前記領域の地表付近で取得される受動的地震信号のスペクトル密度を表す複数のスペクトルであって、それぞれが運動を表す信号から作成される複数のスペクトルを作成するモジュール（４１）と、
記録および周波数に関連付けられたスペクトル属性のセットを取得するのに適した、前記スペクトルそれぞれに現れる周波数ごとに少なくとも１つのスペクトル属性を決定するモジュール（４２）と、
前記スペクトル属性のセットを、各行が１つの記録に関連付けられた行列に編成するモジュール（４３）と、
主成分を決定して前記主成分から前記地下の特性を推定するために、主成分分析法を前記行列に適用するモジュール（４４）と、
を含む、システム。
【請求項１６】
コンピュータプログラムであって、前記コンピュータプログラムがコンピュータにより実行されたときに請求項１〜１４のいずれか１項に記載の方法のステップを実行するための命令を含んだコンピュータプログラム。
【請求項１７】
請求項１〜１４のいずれか１項に記載の方法のステップを実行するための命令を含んだコンピュータプログラムを格納している、コンピュータ可読のデータ媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、ある領域の地下の特性分析を行う分野全般に関し、詳細には、受動的地震信号を調査することによって、さらに詳細には、低周波の受動的地震信号を調査することによって、ある領域の地下の特性分析を行う分野全般に関する。
【背景技術】
【０００２】
低周波の受動的地震信号とは、０．１ヘルツ（Ｈｚ）〜１０Ｈｚ、さらには０．１Ｈｚ〜４Ｈｚまたは５Ｈｚの範囲の周波数で自然発生した大地の運動を表す信号である。これらの周波数範囲では、大地の運動は、海水面より高い陸地の内部を伝搬する波動を生み出す海の波の動きに関係している場合がある。周波数が１Ｈｚを超える波動は、人間の活動（道路、産業）によっても起こりうる。
【０００３】
本明細書において「受動的」という用語は、ユーザーまたはユーザーが制御するツールによって地震信号が生成されないという意味で使用するものと理解されたい。
【０００４】
また、前述の周波数範囲は、１０Ｈｚ〜１５０Ｈｚの範囲に存在しうる高周波の地震信号であって、特に誘発地震の観察に関係する高周波の地震信号に関する周波数範囲とは明確に異なることも理解されたい。
【０００５】
文献“Phenomenology of tremor-like signals observed over hydrocarbon reservoirs”（S． Dangel et al．，Journal of Volcanology and Geothermal Research）の記載によると、信号のスペクトルにおいて１．５Ｈｚ〜４Ｈｚの範囲にある周波数ピークは、炭化水素（ひいては複数種類の流体、または同一種類の流体の複数の相）を含む地下貯留層の上で観察される。その結果、そのような信号を観察することで、複数種類の流体または単一種類の流体の複数の相が領域の地下にあるかどうかを判断することが可能となる。
【０００６】
そのような状況下において、低周波の受動的地震信号を調査することにより、炭化水素（例えば天然ガス）、蒸気、または様々な種類の気体（例えばＣＯ2、Ｈ2）の貯蔵に利用される貯留層における変化を追跡でき、さらには炭化水素または地熱地帯における水および蒸気を含む帯水層について探査する作業を実施できる。
【０００７】
従来技術では、以下の文献が知られている。
【０００８】
文献ＵＳ２００８／００２１６５６には、受動的地震信号を得て、それらの信号のスペクトルを取得し、スペクトルの鉛直成分とスペクトルの水平成分の比を算出し、その比を積分するという、地震データを処理する方法が記載されている。この文献に記載された方法を実施するには、処理が実施される周波数範囲を選択する必要がある。この解決手段は融通性が十分でなく、この方法で得られた結果の正確性が十分でない。
【０００９】
文献ＷＯ２０１０／０８０３６６では、受動的地震データを別の種類の地球物理学的データ（例えば能動的地震データ）と組み合わせて使用することで炭化水素を検知する方法が開示されている。
【００１０】
文献ＷＯ２００９／０２７８２２では、地震データを取得し、信号からスペクトルを取得し、それらのスペクトルの極大値を求めてマップを得る、という、炭化水素貯留層の位置を特定する方法が提案されている。
【００１１】
また、文献ＥＰ２０３００４６には、スペクトルの振幅の比を調査するという方法が記載されている。また、スペクトルの平滑化についても記載されている。
【００１２】
文献ＥＰ１９６０８１２および文献ＥＰ２０３００４６には、地震データを取得し、該地震データからスペクトルを取得し、該スペクトルを平滑化する、という方法が記載されている。
【００１３】
文献ＷＯ２００９／０８１２１０には、地表下における運動を測定することで取得したスペクトルの周波数帯においてエネルギーを特定するという方法が記載されている。
【００１４】
文献ＵＳ２０１４／０２５４３１９には、受動的地震信号の取得について記載されている。この文献では０〜５Ｈｚの範囲の周波数用のローパスフィルタを使用することが提案されているが、その理由は、受動的地震信号に基づいて現象が観察されるのがそれらの周波数においてであるためである。この文献では、信号は周波数―波数領域に変換される。
【００１５】
文献ＷＯ２０１４／１０８８４３には、微小地震信号を得て畳み込んでからフィルターを適用する、という方法が記載されている。
【００１６】
これらの文献に記載された解決手段は満足のいくものではない。これらの解決手段では、領域の地下に異なる種類の流体または相があるかどうかを十分正確に特定することができない。具体的には、これらの方法は人間が生み出すノイズに非常に影響されやすく、流体が存在する領域と他の領域とを正確に区別するのが難しいことがある。
【００１７】
一部の方法では、観察の対象となる周波数範囲をあらかじめ決定しておく必要がある。その解決手段は、観察の対象となる周波数範囲外に現れるアノマリーを考慮していないため、満足のいくものではない。
【発明の概要】
【発明が解決しようとする課題】
【００１８】
本発明は、特に、これらの欠点を軽減することを目的とする。
【課題を解決するための手段】
【００１９】
本発明は、ある領域の地下の特性分析を行う方法であって、
前記領域内の受動的地震信号が記録される少なくとも１つの地点において前記領域の地表付近で取得される受動的地震信号のスペクトル密度を表す複数のスペクトルであって、それぞれが運動を表す信号から作成される複数のスペクトルを作成するステップと、
記録および周波数に関連付けられたスペクトル属性のセットを取得するため、前記スペクトルそれぞれに現れる周波数ごとに少なくとも１つのスペクトル属性を決定するステップと、
前記スペクトル属性のセットを、各行が１つの記録に関連付けられた行列に編成するステップと、
主成分を決定して前記主成分から前記地下の特性を推定するために、主成分分析法を前記行列に適用するステップと、
を含む方法を提案することによって上記の要求を満たす。
【００２０】
運動を運動速度信号として、または運動加速度信号として感知してもよいことがわかるであろう。好ましくは、運動速度信号を使用する。また、３つの方向（１つの鉛直方向と２つの水平方向）に運動が生じうること、および、鉛直運動のみに関連付けられたスペクトルおよび／または水平運動のみに関連付けられたスペクトルを作成できることもわかるであろう。
【００２１】
結果として得られる行列の次元数は、記録の数に属性ごとの周波数の数（属性１つ当たりの周波数の数は属性により異なってもよい）を掛けた数である。
【００２２】
本発明者らは、主成分分析法を使うと、別々の記録に対応するスペクトル間の差異をより容易に示す主成分が得られることを確認した。具体的には、主成分空間はスペクトル間の差異を表す最良の空間である。結果として、この空間は上記領域の地下の特性を推定するのに良い空間となる。
【００２３】
複数の地点を調査する場合、主成分を使用すると、領域のグラフィックによる表現を得ることができ、そのグラフィックによる表現は地中にある流体の存在を示す。１つの地点のみを異なる時点に取得した複数の信号を用いて調査する場合は、流体の量が経時的に変化する限りにおいて、経時的な変化を観察することができる。
【００２４】
通常、主成分分析は当業者によりＰＣＡという略称で呼ばれる。
【００２５】
よって、取得されるスペクトルはサンプルであり、有限な数の周波数を呈するものと理解できる。これらの周波数は、例えば０．１Ｈｚ〜４Ｈｚまたは５Ｈｚといった、広い周波数範囲において選択されてもよい。文献ＵＳ２００８／００２１６５６に記載の方法ではより狭い周波数範囲を定義する必要があるが、本発明を実施するのにその必要はない。
【００２６】
属性は、スペクトルの周波数を特徴づけるパラメーターであって、かつ、流体を含んでいる範囲の鉛直上方に上記領域の上記地点がある場合に変化する特定のパラメーターであってもよいこともわかるであろう。よって、本発明は、複数の属性を使用している点で、使用する属性が鉛直運動と水平運動の比のみである文献ＵＳ２００８／００２１６５６に記載の解決手段と異なる。本発明者らは、複数の属性を使用することで流体の存在を効果的に検知できることを確認した。
【００２７】
さらに本発明では、カナダのサプライヤーＮａｎｏｍｅｔｒｉｃｓにより商品名「Ｔ−４０」として販売されている装置などの地震計を用いて受動的地震信号を取得してもよい。そのような装置は領域の地表付近、例えば約５０センチメートルの深さに埋めてもよい。あるいは、地表が大地と良好につながっているならば、地表に置いてもよい。
【００２８】
ある特定の実施形態では、前記運動は鉛直運動および／または水平運動であり、周波数ごとの前記スペクトル属性は、鉛直地震動のスペクトル密度と水平地震動のスペクトル密度の比、前記水平地震動の周波数に関するスペクトル密度の導関数、および、前記鉛直地震動の周波数に関するスペクトル密度の導関数からなる群より選択される種類のものである。
【００２９】
本発明者らは、これら複数のスペクトル属性を組み合わせることで流体の存在の効果的な検知を実現できることを確認した。
【００３０】
ある特定の実施形態では、前記水平地震動の周波数に関する前記スペクトル密度の前記導関数、および／または、前記鉛直地震動の周波数に関する前記スペクトル密度の前記導関数は、選択した数のスペクトルの点に関して線形回帰を適用することで算出される。
【００３１】
選択される点は、周波数軸を０．５Ｈｚの幅に分割することで得てもよい。
【００３２】
ある特定の実施形態では、信号から各スペクトルを作成することは、
前記信号を、期間がすべて同じである複数の連続するサブ信号に分割することと、
前記サブ信号ごとにスペクトル密度のサブスペクトルを作成することと、
前記サブスペクトルの各周波数について、前記サブスペクトルそれぞれにおける前記周波数についてのスペクトル密度の値から、スペクトル密度の統計学的属性（例えば中央値）を求めることと、
すべての前記周波数のすべての前記統計学的属性から作成される前記スペクトルを取得することと、を含む。
【００３３】
この特定の実施形態は、中央値など統計学的属性が活用されるので、スペクトルを良好に平滑化するのに役立つ。
【００３４】
ある特定の実施形態では、前記サブ信号はそれぞれ、先行するサブ信号の０を超える期間、例えば前記先行するサブ信号の期間の５０％、にわたり前記先行するサブ信号に重なる。
【００３５】
ある特定の実施形態では、前記記録は、または場合によってはすべての前記記録は、前記領域の異なる地点で作られ、前記スペクトル属性のセットの各属性は地点にも関連付けられている。換言すると、記録自体が記録ごとに異なりうる地点に関連付けられているので、属性のセットの各属性は周波数および地点に関連付けられている。
【００３６】
ある特定の実施形態では、前記記録は所定の時点から所定の期間にわたり作られる。
【００３７】
開始時間を真夜中（例えば、午前０時に開始）とし、期間は約４時間とすることが好ましい。これにより、人間が生み出すノイズが最も低い期間に信号を得ることが可能になる。
【００３８】
ある特定の実施形態では、同時に作られた記録ごとのグループとして作られ、前記グループはそれぞれ、前記グループの前記記録が作られた日にちに対応し、前記記録は前記領域の異なる地点で、および／または、異なる時点から作られる。
【００３９】
この特定の実施形態は、限られた数の測定装置を使って領域を調査するのに役立つ。この特定の実施形態では、良好な解像度で領域をカバーできるよう、装置を日ごとに移動させる。
【００４０】
ある特定の実施形態では、前記行列における所定の属性に関連付けられた列は、すべて隣接している。この特定の実施形態では、主成分分析を利用することで、ある属性の影響を別の属性よりも観測しやすくなる。
【００４１】
ある特定の実施形態では、前記記録のグループはそれぞれ、前記行列の複数行からなる１つのグループに関連付けられ、前記属性の値は、複数行からなるグループごとに正規化される。
【００４２】
行列において、これらの複数行からなるグループは、それらを共通のグループに連続的に置くことでグループ化してもよいこともわかるであろう。
【００４３】
したがって、グループはそれぞれ１つの日にちに関連付けられているため、センサーの位置を変更することにより記録で得た信号の振幅の違いが生じる場合でも、正規化によって同じ数値範囲内にある属性値を取得できる。正規化は低次元化と中心化による正規化（une normalisation centree reduite）であってもよい。
【００４４】
ある特定の実施形態では、前記主成分は射影作用素であり、前記行列を射影作用素それぞれに射影して、地点ごとに前記行列を射影した結果を示す前記領域のグラフィック画像を、射影作用素ごとに取得する。
【００４５】
ある特定の実施形態では、Ｋ個の射影作用素を前記射影作用素の中から決定する。換言すると、射影作用素の数は少なくとも１つの基準に応じて選択される。例として、この基準は、領域の地下または地下における経時的変化についての良好なグラフィックによる表現がある射影作用素により得られると判断できるようにするものであってもよい。当業者であれば、厳密には本発明の一部となる手段を除いた手段で取得される地下に関係するデータに応じて、つまり、領域の大地の構造に関する知識に応じて、どのようにこの基準を評価するかがわかる。
【００４６】
例えば、グラフィックによる表現が領域のグレースケールのマップである場合、Ｋ個の射影作用素を、変化が予期される位置にアノマリーが現れるかどうかを判断することで決定してもよい。なお、アノマリーは地質学上のトラップの構造形状に対応していてもよい。
【００４７】
この特定の実施形態は、他の手段（坑井など）により実施した初期の分析から範囲が分かっている貯留層の変化の調査に特に適している。
【００４８】
ある特定の実施形態では、本方法は、
長さＫのベクトルを前記行列の行ごとに得るため、前記行列を前記Ｋ個の射影作用素に射影することと、
長さＫの前記ベクトルから第２の行列を得ることと、
Ｎ個の、複数行のクラスを得るため、一次元または二次元で編成された分類法を前記第２の行列に適用することと、
前記領域の前記地下のアノマリーの規模を表す前記行列の各行に少なくとも１つのクラス番号を割り当てることと、
複数行のクラスごとにクラスのヘッドを作成することと、
Ｎ×Ｋ次元の第３の行列を前記クラスのヘッドから得ることと、
スペクトル属性の第１の行列の各行に現れる周波数の数のＮ倍の次元の第４の行列を得るため、前記第３の行列に擬似逆行列法を適用することと、
をさらに含む。
【００４９】
上記編成された分類法は、自己組織化写像（ＳＯＭ）または生成位相写像（ＧＴＭ）など、当業者によく知られている方法であってもよい。
【００５０】
Ｎは厳密に１を超えるよう選択され、約１０程度、または数十程度であってもよい。
【００５１】
一次元分類では、アノマリー（すなわちスペクトルの変化）の規模を表すため、クラス番号が各行に割り当てられる。これらのクラス番号はグレースケールに表してもよい。異なる行が同じクラス番号を有していてもよい。
【００５２】
二次元分類を使用する場合、例えば補色のスケールを２つ組み合わせるなど、他の表現方法を使用してもよい。
【００５３】
各クラス番号について、そのクラス番号が割り当てられたすべての行の重心としてクラスのヘッドを定義できることがわかるであろう。あるいは、クラスのヘッドはその重心に最も近い行であってもよく、例として、ユークリッド計量によって求めてもよい。
【００５４】
あるいは、本方法は、
長さＫのベクトルを前記行列の行ごとに得るため、前記行列を前記Ｋ個の射影作用素に射影することと、
長さＫの前記ベクトルから第２の行列を得ることと、
スペクトル属性の第１の行列と同じ次元を有する第３の行列を得るため、前記第２の行列に擬似逆行列法を適用することと、
Ｎ個の、複数行のクラスを得るため、一次元または二次元で編成された分類法を前記第３の行列に適用することと、
前記領域の前記地下のアノマリーの規模を表す前記行列の各行に少なくとも１つのクラス番号を割り当てることと、
複数行のクラスごとにクラスのヘッドを作成することと、
スペクトル属性の前記第１の行列の各行に現れる周波数の数のＮ倍の次元の第４の行列を前記クラスのヘッドから得ることと、
をさらに含む。
【００５５】
この変形例では、擬似逆行列のステップを分類法の適用前に実施する。
【００５６】
また、本発明は、ある領域の地下の特性分析を行うシステムであって、
前記領域内の受動的地震信号が記録される少なくとも１つの地点において前記領域の地表付近で取得される受動的地震信号のスペクトル密度を表す複数のスペクトルであって、それぞれが運動を表す信号から作成される複数のスペクトルを作成するモジュールと、
記録および周波数に関連付けられたスペクトル属性のセットを取得するのに適した、前記スペクトルそれぞれに現れる周波数ごとに少なくとも１つのスペクトル属性を決定するモジュールと、
前記スペクトル属性のセットを、各行が１つの記録に関連付けられた行列に編成するモジュールと、
主成分を決定して前記主成分から前記地下の特性を推定するために、主成分分析法を前記行列に適用するモジュールと、
を含む、システムを提供する。
【００５７】
本システムは、本方法の上述した実施形態をすべて実施するように構成されていてもよい。
【００５８】
また、本発明は、コンピュータプログラムであって、前記コンピュータプログラムがコンピュータにより実行されたときに先に定義したような本方法のステップを実行するための命令を含んだコンピュータプログラムを提供する。
【００５９】
本開示に記載されるコンピュータプログラムは、任意のプログラミング言語を使用してよく、ソースコード、オブジェクトコード、または、部分的にコンパイルした形式もしくは他の望ましい形式など、ソースコードとオブジェクトコードの中間のコードでもよいことがわかるであろう。
【００６０】
また、本発明は、先に定義したような本方法のステップを実行するための命令を含んだコンピュータプログラムを格納している、コンピュータ可読のデータ媒体を提供する。
【００６１】
本開示に記載されるデータ（または記録）媒体は、上記プログラムを格納可能な任意のエンティティまたはデバイスでもよい。媒体としては、例えば、読み出し専用メモリ（ＲＯＭ）（例えばコンパクトディスク（ＣＤ）ＲＯＭまたは超小型電子回路ＲＯＭ）、または磁気記録手段（例えばフロッピーディスクまたはハードディスク）などの、格納手段が含まれていてもよい。
【００６２】
さらにデータ媒体は、電気ケーブルまたは光ケーブルを介して、無線通信またはその他の手段によって伝達するのに適した、電気信号または光信号などといった、伝送可能な媒体に該当してもよい。本発明のプログラムは特に、インターネット型のネットワークからダウンロードされてもよい。
【００６３】
あるいは、データ媒体は、上記プログラムが組み込まれた集積回路であって、本方法を実行するのに適した、または本方法を実行するときに使用するのに適した、集積回路に該当してもよい。
【００６４】
本発明の他の特徴および利点は、限定的ではない例を示している添付の図面を参照してなされる以下の説明から明らかになる。
【図面の簡単な説明】
【００６５】
【図１】図１は、本発明の一実施形態における方法のステップを示す図である。
【図２】図２は、本発明の一実施形態におけるシステムを示す図である。
【図３】図３は、ある領域の地下の断面図である。
【図４】図４は、信号からスペクトルをどのように取得するかを示す図である。
【図５】図５は、属性が編成された行列を示す図である。
【図６】図６は、射影作用素への行列の射影を示す図である。
【図７】図７は、再射影したベクトルの分類を示すグラフである。
【図８】図８は、再射影した各ベクトルにパラメーターを関連付けた後に得られるグラフィックによる表現を示す図である。
【発明を実施するための形態】
【００６６】
以下、ある領域の地下の特性分析を行う方法およびシステムについて、本発明の特定の実施形態に即して説明する。
【００６７】
図１は、ある領域の地下の特性分析を行う方法における各ステップを示す図である。
【００６８】
本方法は、複数種類の流体または単一種類の流体の複数の相が領域の地下に存在するかを判断するために実施されうる。
【００６９】
本方法を使用する典型的な用途は、例えば、炭化水素（例えば天然ガス）、蒸気、および様々な種類の気体（例えばＣＯ2、Ｈ2）を含んだ貯留層のモニタリング、炭化水素の探査、および地熱地帯での探査に関する。
【００７０】
本方法の第１ステップＥ０１では、領域内の受動的地震信号が記録される少なくとも１つの地点において該領域の地表付近で取得される受動的地震信号のスペクトル密度を表す複数のスペクトルを作成し、該スペクトルは運動を表す信号からそれぞれ作成される。
【００７１】
すなわち、水平運動を表す信号（異なる方向における２つの信号であり得る）および／または鉛直運動を表す信号を取得するステップをあらかじめ実施する。これらの信号は、カナダのサプライヤーＮａｎｏｍｅｔｒｉｃｓにより商品名「Ｔ−４０」として販売されている装置などの地震計を使用することで得てもよい。そのような装置は、以下にて図３を参照しつつ説明するように、領域の地表付近または該領域の地表に複数個規則的に配置できる。これらの装置は、人間が生み出すノイズを減らすため、夜間に使用することが好ましい。信号はすべて、個々の時点および／または調査対象の領域の個々の位置もしくは地点に関連している。
【００７２】
前述のような装置は、運動を表す運動速度信号を供給することがわかるであろう。
【００７３】
各スペクトルは対応する信号から、その信号のパワースペクトル密度（ＰＳＤ）を求めることで取得してもよい。
【００７４】
また、以下で図４を参照しつつ説明するように、取得したスペクトルを平滑化するための処理を実施できる。
【００７５】
ステップＥ０１の後、複数のスペクトルを取得するが、これらスペクトルは、水平運動に関連したスペクトルと鉛直運動に関連したスペクトルの両方を含んでもよく、これらのスペクトルは記録に関連付けられ、ひいては、領域の地点および取得した時点または日にちで構成される該記録のプロパティに関連付けられる。
【００７６】
水平運動の１つのスペクトルを、異なる２つの方向における運動を表す２つの信号から、各方向に対応するスペクトルの幾何平均を算出することで取得してもよいことがわかるであろう。具体的には、以下の式によって算出する。
【００７７】
【数１】
2019519799.tif 000003
【００７８】
式中、ｆは周波数であり、Ｆｈ（ｆ）は水平運動に対応するスペクトルであり、ＰＳＤEは第１の方向（具体的には東）の水平運動に対応するスペクトルであり、ＰＳＤNは第２の方向（具体的には北）の水平運動に対応するスペクトルである。
【００７９】
よって、スペクトルはサンプリングされ、予め選択された広い範囲に含まれる有限な数の周波数に関係する。
【００８０】
第２のステップＥ０２では、スペクトル属性を決定する。これらの属性は、鉛直地震動のスペクトル密度と水平地震動のスペクトル密度の比、該水平地震動の周波数に関するスペクトル密度の導関数、および、該鉛直地震動の周波数に関するスペクトル密度の導関数からなる群より選ばれてもよい。
【００８１】
これら３つの属性は以下のように書き表せる。
【００８２】
【数２】
2019519799.tif 000004
【００８３】
式中、Ｆｖ（ｆ）は鉛直運動に対応するスペクトルである。
【００８４】
該水平地震動の周波数に関する該スペクトル密度の該導関数、および／または、該鉛直地震動の周波数に関する該スペクトル密度の該導関数は、選択した数のスペクトルの点に関して線形回帰を適用することで算出できることがわかるであろう。
【００８５】
なお、選択される点は、周波数軸を０．５Ｈｚ幅で分割することで得てもよい。
【００８６】
ステップＥ０２では、複数の属性のセットであって、周波数（属性の種類によって異なってもよい）および記録に関連付けられ、ひいては領域の地点および信号を得た時点および／または日にも関連付けられた複数の属性のセットを取得する。
【００８７】
ステップＥ０３では、この属性のセットを、行それぞれが１つの記録（つまり、領域の１つの地点、および、記録が作られた１つの時点および／または１つの日にち）に関連付けられた行列に編成する。
【００８８】
この編成については、図５を参照しつつより詳細に説明する。
【００８９】
ステップＥ０４では、主成分を決定して該主成分から地下の特性を推定するために、主成分分析法を上記成分に適用する。
【００９０】
図２は、図１を参照しつつ説明したステップＥ０１〜Ｅ０４を実施するのに好適なシステム１を示す図である。
【００９１】
システム１はコンピュータシステムでもよく、プロセッサー２およびメモリ３を備える。
【００９２】
コンピュータプログラム４の命令がメモリ３に格納されている。コンピュータプログラム４は、ステップＥ０１を実施するための命令４１、ステップＥ０２を実施するための命令４２、ステップＥ０３を実施するための命令４３、およびステップＥ０４を実施するための命令４４を含んでいる。
【００９３】
命令４１〜４４とプロセッサーとは、組み合わさることでシステム１の複数のモジュールを形成し、該複数のモジュールはステップＥ０１〜Ｅ０４の実施にそれぞれが適合している。
【００９４】
図３は、本発明の方法を実施することで特性分析を行うことが望まれる領域の地下の断面図である。
【００９５】
この目的のために、地震計１００が該領域の地表付近に埋められており、断面の平面にグループ１０１の一部をなす地震計１００を見ることができる。例として、地震計１００は約５０センチメートルの深さに埋めた。こうした設置は、専門家にとって非常に簡単なことである。
【００９６】
あるいは、地震計が地面にしっかり据え付けできる構成ならば、地表に地震計を置くこともできる。当業者であれば、しっかり据え付けるためにはどのように地震計を置いたらよいかがわかる。
【００９７】
この例において、領域の地下は、気体を含む範囲２００および水を含む範囲３００を含んでいる。この領域は貯留層であってもよい。相が異なるこれらの２種類の流体が存在することにより、本発明の方法を実施することが可能となる。
【００９８】
図示した例では、グループ１０１において、図の中央に配置された地震計１００は、該中央の地震計のみが貯留層の鉛直上方にあるので、右側および左側に配置された地震計１００のスペクトルとは異なるスペクトルを示すことになる。
【００９９】
少数の装置を用いて領域内部で記録を作れるようにするため、グループ単位で測定することが可能である。
【０１００】
例えば、１日目の午前０時から午前４時までは、地震計１００を配置してグループ１０１としデータを得る。２日目の午前０時から午前４時までは、地震計１００を配置してグループ１０２としデータを得る。３日目の午前０時から午前４時までは、地震計１００を配置してグループ１０３としデータを得る。４日目の午前０時から午前４時までは、地震計１００を配置してグループ１０４としデータを得る。
【０１０１】
図４は、例えば図３を参照しつつ説明した地震計１００によって取得した信号などの信号からスペクトルをどのように取得するかを示す図である。
【０１０２】
この図には、運動、具体的には１つの方向に沿った運動速度、を表す信号ＳＩＧが見られる。この信号は、午前１２時から開始した４時間の取得期間に得たもので、こうすることで人間が生み出すノイズの出現が減少する。
【０１０３】
信号ＳＩＧは、期間がすべて同じで連続した複数のサブ信号に分割してもよく、この例における複数のサブ信号はそれぞれ、先行するサブ信号の期間の少なくとも半分にわたり該先行するサブ信号に重なる（この重なりは必須ではない）。図において、サブ信号は信号ＳＩＧの下の波括弧によって表されている。
【０１０４】
過度のノイズを示すサブ信号がある場合、それらのサブ信号を省略し、以後処理しなくてもよいことがわかるであろう。例として、閾値を超える絶対値を有するサブ信号を排除することができる。例えば、信号ＳＩＧ全体に対して定義した９９％分位点を超えた絶対値を有するサブ信号を除外できる。
【０１０５】
その後、サブスペクトルをサブ信号ごとに作成する。図にはＰＳＤ＿１、ＰＳＤ＿２、およびＰＳＤ＿３の３つのスペクトル密度サブスペクトルが示されている。
【０１０６】
サブスペクトルの各周波数について、各サブスペクトルの該周波数に対するスペクトル密度の値に基づき、スペクトル密度の中央値を決定する。次に、すべての中央値からスペクトルＰＳＤ＿ｍを取得する。すなわち、スペクトルはこれらの中央値からできている。
【０１０７】
図５は、異なる地点においてある日にち（複数の異なる日にちであってもよい）に取得した信号について得た属性を編成するために、どのように属性を行列Ｍに編成するかを示している。
【０１０８】
行列Ｍでは、以下の表記を使用する。
Ａｔｔｒ＿ｉ：タイプｉの属性
ｘ＿ｊ：領域の地点ｊ（この地点は測定日に関連付けられている）
ｆ＿ｋ：スペクトル属性の周波数ｋ
【０１０９】
行列Ｍにおいて、各行は１つの記録および領域の１つの地点ｘ＿ｊに関連付けられており、各列は属性タイプＡｔｔｒ＿ｉおよびスペクトル属性周波数ｆ＿ｋに関連付けられている。
【０１１０】
行列における所定の属性に関連付けられた列は、すべて隣接している。よって、同時に作られた記録のグループ（例えば、図３を参照しつつ説明したグループ１０１〜１０４で作られた記録のグループ）に対応する行列の行はすべて行列においてグループ化され、複数の、複数行からなるグループを形成する。この例では、各グループは１つの日にちに関連付けられている。
【０１１１】
日にちごとまたは行のグループごとに、属性の値を正規化することが好ましい。
【０１１２】
行列Ｍに属性を編成することで、各行を個体とし各列を変数とする主成分分析法を実施することが可能になる。主成分分析により、「射影作用素」と呼ばれる主成分が取得される。よって、射影作用素は、ｐと表記されるベクトルであり、異なる属性タイプの数に属性タイプごとに存在する周波数の数を掛けた積に等しい長さＬのベクトルである。
【０１１３】
行列Ｍの添え字ｉ（１〜ｍの範囲にあり、ｍは行列Ｍの行の数である）の行の射影作用素ｐへの射影は、以下のように算出される。
【０１１４】
【数３】
2019519799.tif 000005
【０１１５】
この射影の結果は領域の１つの地点に対応するので、信号の取得が実施された領域のすべての地点で、行列Ｍの射影のグラフィックによる表現を取得できる。
【０１１６】
図６にそのようなグラフィックによる表現を示す。この図には、射影ＰＲＪ１、ＰＲＪ２、ＰＲＪ３、およびＰＲＪ４に対応する４つのグラフィックによる表現が示されている。これらのグラフィックによる表現上において、既知の貯留層ＲＥＳの輪郭が強調表示されている。
【０１１７】
予測された空間表現に最もよく一致するグレースケールの変化を示すグラフィックによる表現が、良好な射影作用素に関連付けられているとみなされる。当業者であれば、これらのマップの解釈の仕方がわかる。
【０１１８】
領域の各グラフィックによる表現の下には、周波数の関数として射影作用素自体も示されている。
【０１１９】
図示した例の場合、射影作用素ＰＲＪ１およびＰＲＪ２が良い射影作用素であるとみなされる。よって、Ｋ個、ここでは２個の射影作用素が決定され、これら射影作用素はｐ１およびｐ２と表記される。
【０１２０】
その後、下記式を使用することで、行列を２つの射影作用素ｐ１およびｐ２に射影できる。
【０１２１】
【数４】
2019519799.tif 000006
【０１２２】
これにより、二次元空間に属するｍ個（行列中の行の数、または記録の数）のベクトルを得る。これらのベクトルを編成して第２の行列の行を形成してもよい。
【０１２３】
図７は、この第２の行列に存在する個体を、軸ｘ１、ｘ２によって表される第１の基準座標系に示している。個体はそれぞれ図中の×印に対応する。
【０１２４】
したがって、１次元に編成された分類法をＮ個のクラスを得るために適用する場合、個体を近似するξと表記される曲線を求めた後、円で表されるクラスであって、曲線ξに沿った曲線状の横座標がそのクラス番号（この例では?１から１まで変化する）に対応するクラスを決定することで、個体を分類できる。円の半径が円の共分散の半径に対応することがわかるであろう。
【０１２５】
また、この図は、保持されている２つの射影作用素に対応する軸も示している。該軸はｅ１およびｅ２として表記されるが、アノマリーを十分よく表すことはできない。
【０１２６】
この例において、クラス番号は領域の地下のアノマリーの規模を表す。
【０１２７】
よって、円の中心はクラスのヘッドと見なされる。
【０１２８】
次いで、Ｎ×Ｋ次元の第３の行列をクラスのヘッドから得ることができる。その後、スペクトル属性の第１の行列の各行に現れる周波数の数のＮ倍の次元の第４の行列を得るために、擬似逆行列を求めることができる。
【０１２９】
図８は、領域の各地点のクラス番号によってアノマリーの規模の値を表示することで得られるマップを示し、このマップは擬似逆行列を求めた後に得られる。
【０１３０】
この図では、各クラスのヘッドについて、クラスのヘッドに対応する属性の周波数に応じた変化を示す曲線を２つの属性に関して示しており、この２つの属性とは、鉛直運動に対応するスペクトルの周波数に関する導関数、および、鉛直運動と水平運動の比である。
【０１３１】
これらの曲線が、１Ｈｚ未満の周波数範囲内も含めて、大いに異なる挙動を見せることがわかるであろう。
【０１３２】
本発明者らは、貯留層の鉛直上方で測定する場合、実際には、文献“Phenomenology of tremor-like signals observed over hydrocarbon reservoirs”（S． Dangel et al．，Journal of Volcanology and Geothermal Research）で特定されたピークに先行してスペクトルの減少があることを確認した。本発明は、広い周波数範囲を用い、かつ主成分分析を用いることにより、先行技術の解決手段よりもアノマリーをはっきり提示することができる。
【０１３３】
また、本明細書に示される様々なグラフは、既知の貯留層の上方および周辺で行った測定により得られたこともわかるであろう。
【図１】
2019519799.tif 000007
【図２】
2019519799.tif 000008
【図３】
2019519799.tif 000009
【図４】
2019519799.tif 000010
【図５】
2019519799.tif 000011
【図６】
2019519799.tif 000012
【図７】
2019519799.tif 000013
【図８】
2019519799.tif 000014
【国際調査報告】
2019519799.tif 000015
2019519799.tif 000016
2019519799.tif 000017
2019519799.tif 000018