地球物理学におけるWolfram Mathematica

ブログの作者に感謝 アントン・エキメンコ 彼の報告のために

導入

このメモは会議後に書かれたものです Wolframロシア技術カンファレンス そして私が行ったレポートの要約が含まれています。 このイベントは2016月にサンクトペテルブルクで開催された。 私がカンファレンス会場から 2017 ブロック離れたところで働いていることを考えると、このイベントに参加せずにはいられませんでした。 XNUMX年、XNUMX年は学会報告を聞き、今年は発表をさせていただきました。 まず、興味深い（ように見える）トピックが登場しました。 キリル・ベロフそして第二に、制裁政策に関するロシア連邦の法律を長い間研究した結果、私が働いている企業では、XNUMXつものライセンスが発行されました。 Wolfram Mathematica.

私のスピーチの主題に移る前に、このイベントの構成が優れていたことに注目したいと思います。 会議の訪問ページにはカザン大聖堂の画像が使用されています。 大聖堂はサンクトペテルブルクの主要な観光スポットの XNUMX つであり、会議が行われたホールから非常にはっきりと見えます。

サンクトペテルブルク州立経済大学の入り口で、参加者は学生の中から選ばれたアシスタントに出迎えられ、道に迷うことは許されませんでした。 登録中に小さなお土産が配られました（おもちゃ - 点滅するスパイク、ペン、Wolframのシンボルが入ったステッカー）。 昼食とコーヒーブレイクも会議のスケジュールに含まれていました。 グループの壁にあるおいしいコーヒーとパイについてはすでに書きましたが、シェフは素晴らしいです。 この導入部では、イベント自体、その形式、場所がすでにポジティブな感情をもたらしていることを強調したいと思います。

私とキリル・ベロフが作成したレポートは「Wolfram Mathematica を使用して応用地球物理学の問題を解決する」と呼ばれています。 地震データまたは「古代の川が流れていた場所」のスペクトル分析。 レポートの内容は XNUMX つの部分で構成されています。まず、次の XNUMX つの部分で利用可能なアルゴリズムの使用について説明します。 Wolfram Mathematica XNUMXつ目は地球物理データを分析するためのもので、XNUMXつ目は地球物理データをWolfram Mathematica に入力する方法です。

地震探査

まず、地球物理学について簡単に説明する必要があります。 地球物理学は、岩石の物理的性質を研究する科学です。 さて、岩石には電気、磁気、弾性などのさまざまな特性があるため、電気探査、磁気探査、地震探査など、対応する地球物理学の方法があります...この記事では、地震探査についてのみ詳しく説明します。 地震探査は石油とガスを探す主な方法です。 この方法は、弾性振動の励起とその後の調査地域を構成する岩石からの応答の記録に基づいています。 振動は、陸上 (ダイナマイトまたは弾性振動の非爆発性振動源を使用) または海上 (エアガン) で励起されます。 弾性振動は岩盤中を伝播し、異なる性質を持つ層の境界で屈折および反射されます。 反射波は地表に戻り、陸上の受振器 (通常はコイル内に吊り下げられた磁石の動きに基づく電気力学的装置) または海の水中聴音器 (圧電効果に基づく) によって記録されます。 波が到達する時間によって、地層の深さを判断できます。

地震探査船曳航装置


エアガンは弾性振動を励起します


波は岩塊を通過し、水中聴音器によって記録されます


サンクトペテルブルクのブラゴヴェシチェンスキー橋近くの桟橋にある地球物理探査研究船「イワン・グブキン」

地震信号モデル

岩石にはさまざまな物理的性質があります。 地震探査では、弾性特性、つまり弾性振動の伝播速度と密度が主に重要です。 XNUMX つの層が同じまたは類似の特性を持っている場合、波はそれらの層の間の境界を「認識しません」。 層内の波の速度が異なる場合、層の境界で反射が発生します。 特性の差が大きいほど、反射は強くなります。 その強度は反射率係数 (rc) によって決まります。

ここで、ρ は岩石の密度、ν は波の速度、1 と 2 は上層と下層を示します。

最も単純で最も頻繁に使用される地震信号モデルの XNUMX つは畳み込みモデルです。このモデルでは、記録された地震痕跡が一連の反射係数とプローブ パルスの畳み込みの結果として表されます。

ここで、s(t) — 地震の痕跡、つまり一定の録音時間内に水中聴音器または受振器によって録音されたすべてのもの、 w(t) - エアガンによって生成される信号、 n(t) - ランダムなノイズ。

例として合成地震痕跡を計算してみましょう。 初期信号として、地震探査で広く使用されているリッカー パルスを使用します。

length=0.050; (*Signal lenght*)
dt=0.001;(*Sample rate of signal*)
t=Range[-length/2,(length)/2,dt];(*Signal time*)
f=35;(*Central frequency*)
wavelet=(1.0-2.0*(Pi^2)*(f^2)*(t^2))*Exp[-(Pi^2)*(f^2)*(t^2)];
ListLinePlot[wavelet, Frame->True,PlotRange->Full,Filling->Axis,PlotStyle->Black,
PlotLabel->Style["Initial wavelet",Black,20],
LabelStyle->Directive[Black,Italic],
FillingStyle->{White,Black},ImageSize->Large,InterpolationOrder->2]

初期地震推力


300 ms と 600 ms の深さに XNUMX つの境界を設定し、反射係数は乱数になります。

rcExample=ConstantArray[0,1000];
rcExample[[300]]=RandomReal[{-1,0}];
rcExample[[600]]=RandomReal[{0,1}];
ListPlot[rcExample,Filling->0,Frame->True,Axes->False,PlotStyle->Black,
PlotLabel->Style["Reflection Coefficients",Black,20],
LabelStyle->Directive[Black,Italic]]

反射係数の系列


地震痕跡を計算して表示してみましょう。 反射係数の符号が異なるため、地震軌跡上で XNUMX つの交互の反射が得られます。

traceExamle=ListConvolve[wavelet[[1;;;;1]],rcExample];
ListPlot[traceExamle,
PlotStyle->Black,Filling->0,Frame->True,Axes->False,
PlotLabel->Style["Seismic trace",Black,20],
LabelStyle->Directive[Black,Italic]]

シミュレートされたトラック


この例では、予約する必要があります。実際には、層の深さはもちろんメートル単位で決定され、地震痕跡の計算は時間領域で行われます。 深さをメートル単位で設定し、層内の速度を知って到着時間を計算する方が正確です。 今回は早速時間軸上にレイヤーを設定してみました。

現地調査について言えば、そのような観測の結果、膨大な数の同様の時系列（地震痕跡）が記録されます。 たとえば、長さ 25 km、幅 15 km のサイトを調査する場合、作業の結果、各トレースが 25x25 メートルのセル (このようなセルはビンと呼ばれます) を特徴付けるため、最終的なデータ配列には 600000 のトレースが含まれることになります。 サンプリング時間が 1 ミリ秒、録音時間が 5 秒の場合、最終的なデータ ファイルは 11 GB を超え、元の「生」素材のボリュームは数百ギガバイトになる可能性があります。

彼らとどのように協力するか Wolfram Mathematica?

パッケージ 地質学IO

パッケージの開発が始まりました 質問 ロシア語を話す支援グループのVKの壁に。 コミュニティの反応のおかげで、解決策はすぐに見つかりました。 そして結果的にはシリアスな展開に発展してしまいました。 対応する Wolframコミュニティのウォールポスト 司会者からもマークされました。 現在、このパッケージは、地質業界で積極的に使用されている次のデータ タイプの操作をサポートしています。

1. ZMAPおよびIRAP形式の地図データのインポート
2. LAS 形式のウェルでの測定値のインポート
3. 地震ファ​​イル形式の入出力 セギ

パッケージをインストールするには、組み立てられたパッケージのダウンロード ページの指示に従う必要があります。 任意の場所で次のコードを実行します Mathematica ノートブック:

If[PacletInformation["GeologyIO"] === {}, PacletInstall[URLDownload[
    "https://wolfr.am/FiQ5oFih", 
    FileNameJoin[{CreateDirectory[], "GeologyIO-0.2.2.paclet"}]
]]]

その後、パッケージはデフォルトのフォルダーにインストールされます。そのパスは次のように取得できます。

FileNameJoin[{$UserBasePacletsDirectory, "Repository"}]

例として、パッケージの主な機能を示します。 この呼び出しは伝統的にWolfram言語のパッケージに対して行われます:

Get["GeologyIO`"]

パッケージは次を使用して開発されています Wolframワークベンチ。 これにより、パッケージの主な機能にドキュメントを添付することができます。プレゼンテーション形式の点では、Wolfram Mathematica 自体のドキュメントと変わらないため、最初に知っている人向けにテストファイルをパッケージに提供できます。

このようなファイルは、特にファイル「Marmousi.segy」です。これは、フランス石油協会によって開発された、地質セクションの合成モデルです。 このモデルを使用して、開発者は波動場のモデリング、データ処理、地震軌跡の反転などの独自のアルゴリズムをテストします。 Marmousi モデル自体は、パッケージ自体がダウンロードされたリポジトリに保存されます。 ファイルを取得するには、次のコードを実行します。

If[Not[FileExistsQ["Marmousi.segy"]], 
URLDownload["https://wolfr.am/FiQGh7rk", "Marmousi.segy"];]
marmousi = SEGYImport["Marmousi.segy"]

インポート結果 - SEGYData オブジェクト


SEGY 形式には、観測に関するさまざまな情報が保存されます。 まず、これらはテキストコメントです。 これには、作業場所、測定を行った会社名などの情報が含まれます。 この例では、このヘッダーは TextHeader キーを使用したリクエストによって呼び出されます。 短縮されたテキストヘッダーは次のとおりです。

Short[marmousi["TextHeader"]]

「マルムーシのデータセットは研究所で生成されました...最小速度は 1500 m/s、最大速度は 5500 m/s)」

「traces」キーを使用して地震痕跡にアクセスすると、実際の地質モデルを表示できます (このパッケージの特徴の XNUMX つは、キーの大文字と小文字が区別されないことです)。

ArrayPlot[Transpose[marmousi["traces"]], PlotTheme -> "Detailed"]

モデル マルムーシ


現在、このパッケージでは、大きなファイルからデータを部分的にロードすることもできるため、サイズが数十ギガバイトに達するファイルを処理できるようになります。 このパッケージの関数には、データを .segy にエクスポートし、ファイルの末尾に部分的に追加する関数も含まれています。

これとは別に、.segy ファイルの複雑な構造を扱う場合は、パッケージの機能に注目する価値があります。 これにより、キーとインデックスを使用して個々のトレースとヘッダーにアクセスできるだけでなく、それらを変更してファイルに書き込むこともできるためです。 GeologyIO の実装の技術的な詳細の多くはこの記事の範囲を超えており、おそらく別の説明が必要です。

地震探査におけるスペクトル解析の関連性

地震データをWolfram Mathematica にインポートする機能により，実験データに組み込みの信号処理機能を使用できるようになります。 各地震痕跡は時系列を表すため、それらを研究するための主なツールの XNUMX つはスペクトル分析です。 地震データの周波数構成を分析するための前提条件としては、たとえば次のものを挙げることができます。

1. 異なるタイプの波は、異なる周波数構成によって特徴付けられます。 これにより、有用な電波を強調し、干渉波を抑制することができます。
2. 多孔性や飽和度などの岩石の特性は、周波数構成に影響を与える可能性があります。 これにより、最良の特性を持つ岩石を識別することが可能になります。
3. 層の厚さが異なると、異なる周波数範囲で異常が発生します。

XNUMX 番目のポイントは、この記事の文脈における主要なポイントです。 以下は、厚さが異なる層 (ウェッジ モデル) の場合の地震痕跡を計算するためのコードの一部です。 このモデルは、多くの層から反射された波が互いに重なり合うときの干渉効果を分析するために、地震探査で伝統的に研究されてきました。

nx=200;(* Number of grid points in X direction*)
ny=200;(* Number of grid points in Y direction*)
T=2;(*Total propagation time*)
(*Velocity and density*)
modellv=Table[4000,{i,1,ny},{j,1,nx}];(* P-wave velocity in m/s*)
rho=Table[2200,{i,1,ny},{j,1,nx}];(* Density in g/cm^3, used constant density*)
Table[modellv[[150-Round[i*0.5];;,i]]=4500;,{i,1,200}];
Table[modellv[[;;70,i]]=4500;,{i,1,200}];
(*Plotting model*)
MatrixPlot[modellv,PlotLabel->Style["Model of layer",Black,20],
LabelStyle->Directive[Black,Italic]]

ピンチアウトフォーメーションのモデル


くさびの内側の波の速度は 4500 m/s、くさびの外側は 4000 m/s、密度は一定の 2200 g/cmXNUMX と仮定されます。 このようなモデルについては、反射係数と地震痕跡を計算します。

rc=Table[N[(modellv[[All,i]]-PadLeft[modellv[[All,i]],201,4000][[1;;200]])/(modellv[[All,i]]+PadLeft[modellv[[All,i]],201,4500][[1;;200]])],{i,1,200}];
traces=Table[ListConvolve[wavelet[[1;;;;1]],rc[[i]]],{i,1,200}];
starttrace=10;
endtrace=200;
steptrace=10;
trasenum=Range[starttrace,endtrace,steptrace];
traserenum=Range[Length@trasenum];
tracedist=0.5;
Rotate[Show[
Reverse[Table[
	ListLinePlot[traces[[trasenum[[i]]]]*50+trasenum[[i]]*tracedist,Filling->{1->{trasenum[[i]]*tracedist,{RGBColor[0.97,0.93,0.68],Black}}},PlotStyle->Directive[Gray,Thin],PlotRange->Full,InterpolationOrder->2,Axes->False,Background->RGBColor[0.97,0.93,0.68]],
		{i,1,Length@trasenum}]],ListLinePlot[Transpose[{ConstantArray[45,80],Range[80]}],PlotStyle->Red],PlotRange->All,Frame->True],270Degree]

ウェッジモデルの地震痕跡


この図に示された一連の地震痕跡を地震断面と呼びます。 ご覧のとおり、反射波の形状は前に指定したモデルに明確に対応しているため、その解釈は直感的なレベルでも実行できます。 痕跡をより詳細に分析すると、1 から約 30 までの痕跡に違いがないことがわかります。地層の屋根からの反射と底部からの反射は互いに重なりません。 31 番目のトレースから、反射が干渉し始めます。 また、モデルでは反射係数は水平方向には変化しませんが、地層の厚さが変化すると地震痕跡の強度も変化します。

地層の上部境界からの反射の振幅を考えてみましょう。 60 番目のルートから反射の強度が増加し始め、70 番目のルートで最大になります。 このようにして、層の屋根と底部からの波の干渉が現れ、場合によっては地震記録に重大な異常を引き起こすことになります。

ListLinePlot[GaussianFilter[Abs[traces[[All,46]]],3][[;;;;2]],
InterpolationOrder->2,Frame->True,PlotStyle->Black,
PlotLabel->Style["Amplitude of reflection",Black,20],
LabelStyle->Directive[Black,Italic],
PlotRange->All]

ウェッジの上端からの反射波の振幅のグラフ


信号が低周波の場合、地層の厚い部分で干渉が発生し始め、高周波信号の場合、より薄い層で干渉が発生するのは論理的です。 次のコード スニペットは、35 Hz、55 Hz、および 85 Hz の周波数の信号を作成します。

waveletSet=Table[(1.0-2.0*(Pi^2)*(f^2)*(t^2))*Exp[-(Pi^2)*(f^2)*(t^2)],
{f,{35,55,85}}];
ListLinePlot[waveletSet,PlotRange->Full,PlotStyle->Black,Frame->True,
PlotLabel->Style["Set of wavelets",Black,20],
LabelStyle->Directive[Black,Italic],
ImageSize->Large,InterpolationOrder->2]

周波数 35 Hz、55 Hz、85 Hz のソース信号のセット


地震痕跡を計算し、反射波振幅のグラフをプロットすることにより、さまざまな周波数で、さまざまな地層の厚さで異常が観察されることがわかります。

tracesSet=Table[ListConvolve[waveletSet[[j]][[1;;;;1]],rc[[i]]],{j,1,3},{i,1,200}];

lowFreq=ListLinePlot[GaussianFilter[Abs[tracesSet[[1]][[All,46]]],3][[;;;;2]],InterpolationOrder->2,PlotStyle->Black,PlotRange->All];
medFreq=ListLinePlot[GaussianFilter[Abs[tracesSet[[2]][[All,46]]],3][[;;;;2]],InterpolationOrder->2,PlotStyle->Black,PlotRange->All];
highFreq=ListLinePlot[GaussianFilter[Abs[tracesSet[[3]][[All,46]]],3][[;;;;2]],InterpolationOrder->2,PlotStyle->Black,PlotRange->All];

Show[lowFreq,medFreq,highFreq,PlotRange->{{0,100},All},
PlotLabel->Style["Amplitudes of reflection",Black,20],
LabelStyle->Directive[Black,Italic],
Frame->True]

異なる周波数におけるウェッジの上端からの反射波の振幅のグラフ


石油探査の主なタスクの XNUMX つは、井戸を敷設するのに最も有望な地点 (つまり、地層が存在する領域) を評価することであるため、地震観測の結果から地層の厚さについて結論を導き出す能力は非常に役立ちます。厚い）。 さらに、地質学的セクションには、その起源によって地層の厚さが急激に変化する物体が存在する可能性があります。 このため、スペクトル分析はそれらを研究するための効果的なツールになります。 記事の次の部分では、そのような地質学的オブジェクトについてさらに詳しく検討します。

実験データ。 どこでそれらを入手しましたか?そしてその中で何を探すべきですか?

記事で分析された資料は西シベリアで入手されたものです。 おそらく誰もが例外なく知っているように、この地域は我が国の主要な石油生産地域です。 この地域では、前世紀の 60 年代に鉱床の積極的な開発が始まりました。 石油鉱床を探す主な方法は地震探査です。 この地域の衛星写真を見るのは興味深いです。 小さな縮尺では、膨大な数の沼地や湖に注目することができます。地図を拡大すると、クラスター井戸掘削現場が表示され、地図を限界まで拡大すると、地震が発生したプロファイルの空き地を区別することもできます。観察が行われました。

Yandex マップの衛星画像 - ノヤブリスク市街地


田畑の XNUMX つにある井戸パッドのネットワーク


西シベリアの油を含む岩石は、深さ 1 km から 5 km までの広範囲に存在します。 石油を含む岩石の主な量は、ジュラ紀と白亜紀に形成されました。 ジュラ紀については、おそらく同名の映画で多くの人に知られているでしょう。 ジュラ紀の気候 現代のものとは大きく異なっていました。 ブリタニカ百科事典には、各時代を特徴づける一連の古地図が掲載されています。

現在、

ジュラ紀


ジュラ紀の西シベリアの領土は海岸（川と浅い海が交差する土地）であったことに注意してください。 気候が快適だったので、当時の典型的な風景は次のようなものであったと推測されます。

ジュラ紀のシベリア


この絵で私たちにとって重要なのは、動物や鳥ではなく、背景にある川のイメージです。 川は、先ほど立ち寄ったものと同じ地質物です。 実際のところ、川の活動により、よく選別された砂岩が堆積し、石油の貯留層となります。 これらの貯水池は、奇妙で複雑な形状 (川床のような) を持つことがあり、その厚さは変化します。堤防近くの厚さは薄いですが、水路の中心に近づくか蛇行した領域では厚さが増します。 したがって、ジュラ紀に形成された川は現在、深さ約XNUMXキロメートルにあり、石油貯留層の探索の対象となっています。

実験データ。 処理と可視化

記事に示されている地震資料については、すぐに保留しておきます。解析に使用されるデータの量が膨大であるため、記事の本文には元の地震痕跡の一部のみが含まれています。 これにより、誰でも上記の計算を再現できるようになります。

地震データを扱う場合、地球物理学者は通常、さまざまな種類のデータを分析でき、便利なグラフィカル インターフェイスを備えた特殊なソフトウェアを使用します (Petrel や Paradigm など、その開発が積極的に使用されている業界リーダーがいくつかあります)。 非常に便利であるにもかかわらず、この種のソフトウェアには欠点もあります。たとえば、最新のアルゴリズムを安定バージョンで実装するには時間がかかり、計算を自動化する可能性は通常限られています。 このような状況では、コンピュータ数学システムや高級プログラミング言語を使用すると非常に便利になります。これらを使用すると、幅広いアルゴリズム ベースを使用できると同時に、多くのルーチンを引き受けることができます。 これはWolfram Mathematica で地震データを扱うために使用される原理です。 データを操作するインタラクティブな作業のための豊富な機能を記述することは不適切です。一般に受け入れられている形式から確実にロードし、目的のアルゴリズムを適用して、外部形式にアップロードし直すことがより重要です。

提案されたスキームに従って、元の地震データをロードし、それらを表示します。 Wolfram Mathematica:

Get["GeologyIO`"]
seismic3DZipPath = "seismic3D.zip";
seismic3DSEGYPath = "seismic3D.sgy";
If[FileExistsQ[seismic3DZipPath], DeleteFile[seismic3DZipPath]];
If[FileExistsQ[seismic3DSEGYPath], DeleteFile[seismic3DSEGYPath]];
URLDownload["https://wolfr.am/FiQIuZuH", seismic3DZipPath];
ExtractArchive[seismic3DZipPath];
seismic3DSEGY = SEGYImport[seismic3DSEGYPath]

このようにしてダウンロードおよびインポートされたデータは、10km×5kmのエリアに記録されたルートです。 三次元地震探査法（個々の物理的形状に沿って波を記録するのではなく、全域を同時に記録する方法）を用いてデータを取得すると、地震データキューブを取得することが可能になります。 これらは三次元のオブジェクトであり、その垂直断面と水平断面により、地質環境の詳細な研究が可能になります。 この例では、XNUMX 次元データを扱っています。 次のようにテキストヘッダーから情報を取得できます。

StringPartition[seismic3DSEGY["textheader"], 80] // TableForm

C 1 これは GEOLOGYIO パッケージ テストのデモ ファイルです
C 2
C 3
C 4
C 5 日付ユーザー名: WOLFRAM ユーザー
C 6 調査名: シベリアのどこか
C 7 ファイルタイプ 3D 地震ボリューム
C 8
C 9
C10 Z 範囲: 最初の 2200M 最後の 2400M

このデータセットは、データ分析の主要な段階を実証するには十分です。 ファイル内のトレースは順番に記録され、それぞれは次の図のようになります。これは、垂直軸 (深さ軸) に沿った反射波の振幅の分布です。

ListLinePlot[seismic3DSEGY["traces"][[100]], InterpolationOrder -> 2, 
 PlotStyle -> Black, PlotLabel -> Style["Seismic trace", Black, 20],
 LabelStyle -> Directive[Black, Italic], PlotRange -> All, 
 Frame -> True, ImageSize -> 1200, AspectRatio -> 1/5]

地震断面痕跡の一​​つ


調査対象領域の各方向にトレースの数がわかれば、3 次元データ配列を生成し、ImageXNUMXD[] 関数を使用して表示できます。

traces=seismic3DSEGY["traces"];
startIL=1050;EndIL=2000;stepIL=2; (*координата Х начала и конца съёмки и шаг трасс*)
startXL=1165;EndXL=1615;stepXL=2; (*координата Y начала и конца съёмки и шаг трасс*)
numIL=(EndIL-startIL)/stepIL+1;   (*количество трасс по оис Х*)
numXL=(EndXL-startXL)/stepIL+1;   (*количество трасс по оис Y*)
Image3D[ArrayReshape[Abs[traces/Max[Abs[traces[[All,1;;;;4]]]]],{numIL,numXL,101}],ViewPoint->{-1, 0, 0},Background->RGBColor[0,0,0]]

地震データキューブの XNUMXD 画像 (縦軸 - 深さ)


対象の地質学的特徴が激しい地震異常を引き起こす場合は、透明性のある視覚化ツールを使用できます。 記録の「重要でない」領域を非表示にして、異常部分だけを残すことができます。 Wolfram Mathematica ではこれを次のように行うことができます． 不透明度[] и ラスター3D[].

data = ArrayReshape[Abs[traces/Max[Abs[traces[[All,1;;;;4]]]]],{numIL,numXL,101}];
Graphics3D[{Opacity[0.1], Raster3D[data, ColorFunction->"RainbowOpacity"]}, 
Boxed->False, SphericalRegion->True, ImageSize->840, Background->None]

Opacity[] および Raster3D[] 関数を使用した地震データ キューブ画像

合成例と同様に、元の立方体の断面で、さまざまな起伏を持ついくつかの地質学的境界 (層) を識別できます。

スペクトル解析の主なツールはフーリエ変換です。 これを利用すると、各トレースまたはトレースのグループの振幅周波数スペクトルを評価できます。 ただし、データを周波数領域に転送した後は、周波数がいつ変化するか (どの深度で読み取られるか) に関する情報が失われます。 時間 (深さ) 軸上で信号の変化の位置を特定できるようにするために、ウィンドウ フーリエ変換とウェーブレット分解が使用されます。 この記事ではウェーブレット分解を使用します。 ウェーブレット解析技術は 90 年代から地震探査に積極的に使用され始めました。 ウィンドウフーリエ変換に対する利点は、より優れた時間分解能であると考えられます。

次のコード部分を使用すると、地震痕跡の XNUMX つを個々のコンポーネントに分解できます。

cwd=ContinuousWaveletTransform[seismicSection["traces"][[100]]]
Show[
ListLinePlot[Re[cwd[[1]]],PlotRange->All],
ListLinePlot[seismicSection["traces"][[100]],
PlotStyle->Black,PlotRange->All],ImageSize->{1500,500},AspectRatio->Full,
PlotLabel->Style["Wavelet decomposition",Black,32],
LabelStyle->Directive[Black,Italic],
PlotRange->All,
Frame->True]

トレースのコンポーネントへの分解


さまざまな波の到達時間で反射エネルギーがどのように分布するかを評価するには、スカログラム (スペクトログラムに類似) が使用されます。 原則として、実際にはすべてのコンポーネントを分析する必要はありません。 通常、低周波、中周波、高周波成分が選択されます。

freq=(500/(#*contWD["Wavelet"]["FourierFactor"]))&/@(Thread[{Range[contWD["Octaves"]],1}]/.contWD["Scales"])//Round;
ticks=Transpose[{Range[Length[freq]],freq}];
WaveletScalogram[contWD,Frame->True,FrameTicks->{{ticks,Automatic},Automatic},FrameTicksStyle->Directive[Orange,12],
FrameLabel->{"Time","Frequency(Hz)"},LabelStyle->Directive[Black,Bold,14],
ColorFunction->"RustTones",ImageSize->Large]

スカログラム。 関数の結果 ウェーブレットスカログラム[]


Wolfram言語はウェーブレット変換に関数を使用する ContinuousWaveletTransform[]。 そして、この関数のトレースセット全体への適用は、関数を使用して実行されます。 テーブル[]。 ここでWolfram Mathematica の強みのXNUMXつである並列化を使用できることは注目に値します パラレルテーブル[]。 上記の例では、並列化の必要はありません。データの量は大きくありませんが、数十万のトレースを含む実験データ セットを扱う場合、これは必要です。

tracesCWD=Table[Map[Hilbert[#,0]&,Re[ContinuousWaveletTransform[traces[[i]]][[1]]][[{13,15,18}]]],{i,1,Length@traces}]; 

機能適用後 ContinuousWaveletTransform[] 選択した周波数に対応する新しいデータ セットが表示されます。 上の例では、これらの周波数は 38Hz、33Hz、27Hz です。 周波数の選択は、ほとんどの場合、テストに基づいて行われます。さまざまな周波数の組み合わせに対して効果的なマップが取得され、地質学者の観点から最も有益なマップが選択されます。

結果を同僚と共有したり、顧客に提供したりする必要がある場合は、GeologyIO パッケージの SEGYExport[] 関数を使用できます。

outputdata=seismic3DSEGY;
outputdata["traces",1;;-1]=tracesCWD[[All,3]];
outputdata["textheader"]="Wavelet Decomposition Result";
outputdata["binaryheader","NumberDataTraces"]=Length[tracesCWD[[All,3]]];
SEGYExport["D:result.segy",outputdata];

これらの XNUMX つの立方体 (低周波、中周波、および高周波コンポーネント) では、通常、RGB ブレンディングを使用してデータを一緒に視覚化します。 各コンポーネントには、赤、緑、青という独自の色が割り当てられます。 Wolfram Mathematica では，次の関数を使用してこれを行うことができます カラーコンバイン[].

その結果、地質学的解釈が可能な画像が得られます。 この断面に記録されている蛇行により、貯留層である可能性が高く、石油埋蔵量を含む古水路の輪郭を描くことが可能になります。 このような川系の現代の類似物の検索と分析により、蛇行の最も有望な部分を決定することができます。 この水路自体は、よく選別された砂岩の厚い層によって特徴付けられており、石油の優れた貯留層となっています。 「レース」異常の外側の領域は、現代の氾濫原の堆積物に似ています。 氾濫原の堆積物は主に粘土質の岩石で構成されており、これらのゾーンへの掘削は効果がありません。

データ キューブの RGB スライス。 中央（中央やや左）では蛇行する川の流れを辿ることができます。

データ キューブの RGB スライス。 左側には蛇行する川の流れが見えます。


場合によっては、地震データの品質により、非常に鮮明な画像が得られることがあります。 これは、現場での作業方法、ノイズ低減アルゴリズムで使用される機器によって異なります。 このような場合、河川系の断片だけでなく、広がった古川全体も表示されます。

地震データ キューブ (水平スライス) の 2 つのコンポーネントの RGB 混合。 深さ約XNUMXkm。

サラトフ近郊のヴォルガ川の衛星画像

まとめ

Wolfram Mathematica を使用すると，地震データを分析し，鉱物探査に関連する応用問題を解決できます。GeologyIO パッケージを使用すると，このプロセスがより便利になります。 地震データの構造は、組み込みメソッドを使用して計算を高速化するようなものです (パラレルテーブル[], ParallelDo[],…) は非常に効率的で、大量のデータを処理できます。 これは、GeologyIO パッケージのデータ ストレージ機能によって大幅に促進されます。 ちなみに、このパッケージは応用地震探査の分野だけでなく使用することもできます。 地中レーダーと地震学ではほぼ同じ種類のデータが使用されています. 結果を改善する方法についての提案がある場合, Wolfram Mathematica アーセナルのどの信号解析アルゴリズムがそのようなデータに適用できるか, または批判がある場合は, 連絡​​してください.コメント。

出所： habr.com