地球物理学中的 Wolfram 数学

感谢博客作者 安东·埃基门科 为他的报告

介绍

这篇笔记是在会议结束后写的 Wolfram 俄罗斯技术会议 并包含我提供的报告的摘要。 该活动于六月在圣彼得堡举行。 考虑到我的工作地点距离会场只有一个街区，所以我忍不住参加了这次活动。 2016年、2017年我听了会议报告，今年我做了报告。 首先，出现了一个有趣的（在我看来）主题，我们正在开发该主题 基里尔·别洛夫其次，经过长时间研究俄罗斯联邦有关制裁政策的立法，在我工作的企业，出现了多达两张许可证 Wolfram Mathematica.

在进入演讲主题之前，我想先介绍一下这次活动的良好组织。 会议的访问页面使用了喀山大教堂的图像。 大教堂是圣彼得堡的主要景点之一，从会议举行的大厅可以清晰地看到。

在圣彼得堡国立经济大学入口处，学生们的助手迎接了参与者，他们不让他们迷路。 注册期间，我们会赠送一些小纪念品（一个玩具 - 一根闪光的钉子、一支笔、带有 Wolfram 符号的贴纸）。 午餐和茶歇也包括在会议日程中。 我已经在小组的墙上注意到了美味的咖啡和馅饼 - 厨师很棒。 通过这个介绍部分，我想强调的是，活动本身、它的形式和地点已经带来了积极的情绪。

我和 Kirill Belov 准备的报告名为“使用 Wolfram Mathematica 解决应用地球物理学中的问题”。 地震数据或“古代河流流经的地方”的频谱分析。 报告内容分为两部分：一是现有算法的使用。 Wolfram Mathematica 用于分析地球物理数据，其次，这是如何将地球物理数据放入 Wolfram Mathematica 中。

地震勘探

首先，您需要对地球物理学进行简短的了解。 地球物理学是研究岩石物理性质的科学。 那么，由于岩石具有不同的性质：电、磁、弹性，因此有相应的地球物理学方法：电法勘探、磁法勘探、地震勘探……在本文中，我们将只更详细地讨论地震勘探。 地震勘探是寻找石油和天然气的主要方法。 该方法基于弹性振动的激发以及随后记录构成研究区域的岩石的响应。 在陆地上（使用炸药或弹性振动的非爆炸性振动源）或在海上（使用气枪）激发振动。 弹性振动通过岩体传播，在不同性质的层的边界处折射和反射。 反射波返回表面，并由陆地上的地震检波器（通常基于悬挂在线圈中的磁铁的运动的电动装置）或海洋中的水听器（基于压电效应）记录。 通过波浪到达的时间，可以判断地质层的深度。

地震船拖曳设备


气枪激发弹性振动


波穿过岩体并被水听器记录


地球物理调查研究船“伊万·古布金”号停泊在圣彼得堡布拉戈维申斯基大桥附近的码头

地震信号模型

岩石具有不同的物理性质。 对于地震勘探，弹性特性至关重要——弹性振动的传播速度和密度。 如果两层具有相同或相似的属性，那么波“将不会注意到”它们之间的边界。 如果各层中的波速不同，则在各层的边界处会发生反射。 性质差异越大，反射越强烈。 其强度由反射系数 (rc) 决定：

式中，ρ为岩石密度，ν为波速，1和2表示上层和下层。

最简单且最常用的地震信号模型之一是卷积模型，当记录的地震道表示为反射系数序列与探测脉冲的卷积结果时：

其中 s(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深度处设置两个边界，反射系数将是随机数

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]]

反射系数序列


让我们计算并显示地震道。 由于反射系数具有不同的符号，我们在地震道上得到两个交替的反射。

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 公里、宽 15 公里的地点时，工作结果是每条迹线都表征一个 25x25 米的单元（这样的单元称为 bin），最终的数据阵列将包含 600000 个迹线。 采样时间为1毫秒，记录时间为5秒，最终的数据文件将超过11GB，原始“原始”材料的体积可达数百GB。

如何与他们合作 Wolfram Mathematica?

包 地质学IO

开始开发包 问题 在俄语支持小组的 VK 墙上。 感谢社区的响应，很快就找到了解决方案。 结果，它发展成为一项重大的发展。 相应的 Wolfram 社区墙柱 它甚至被版主标记了。 目前，该软件包支持使用地质行业中积极使用的以下数据类型：

1. 导入 ZMAP 和 IRAP 格式的地图数据
2. 以 LAS 格式孔导入测量值
3. 地震文件格式的输入和输出 赛格

要安装该软件包，您必须按照已组装软件包的下载页面上的说明进行操作，即在任意位置执行以下代码 数学笔记本:

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"]]

“Marmousi 数据集是在研究所生成的……最小速度为 1500 m/s，最大速度为 5500 m/s）”

您可以通过使用“traces”键访问地震道来显示实际的地质模型（该包的功能之一是键不区分大小写）：

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

模型 Marmousi


目前，该软件包还允许您从大文件中分批加载数据，从而可以处理大小可达数十GB的文件。 该包的功能还包括将数据导出到 .segy 以及部分附加到文件末尾的功能。

另外，值得注意的是在处理 .segy 文件的复杂结构时该包的功能。 因为它不仅允许您使用键和索引访问各个跟踪和标头，而且还可以更改它们，然后将它们写入文件。 GeologyIO 实现的许多技术细节超出了本文的范围，可能值得单独描述。

频谱分析在地震勘探中的相关性

将地震数据导入 Wolfram Mathematica 的功能使您可以使用内置的信号处理功能来处理实验数据。 由于每条地震道代表一个时间序列，因此研究它们的主要工具之一是频谱分析。 在分析地震数据频率组成的先决条件中，我们可以列举以下内容：

1. 不同类型的波具有不同频率组成的特征。 这使您可以突出显示有用的波并抑制干扰波。
2. 孔隙度和饱和度等岩石特性会影响频率成分。 这使得识别具有最佳特性的岩石成为可能。
3. 不同厚度的层会导致不同频率范围的异常。

第三点是本文的重点。 下面是一个代码片段，用于在厚度变化的层（楔形模型）的情况下计算地震道。 该模型传统上在地震勘探中研究，以分析当多层反射波相互叠加时的干扰效应。

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、55Hz、85Hz 的源信号


通过计算地震道并绘制反射波振幅图，我们可以看到，对于不同频率，在不同地层厚度处观察到异常。

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]

不同频率下楔块上边缘反射波的振幅图


从地震观测结果得出有关地层厚度的结论的能力非常有用，因为石油勘探的主要任务之一是评估最有希望打井的点（即地层所在的区域）。更厚的）。 此外，在地质剖面中，可能存在其成因导致地层厚度急剧变化的物体。 这使得光谱分析成为研究它们的有效工具。 在本文的下一部分中，我们将更详细地考虑此类地质物体。

实验数据。 您从哪里获得它们以及在其中寻找什么？

文章中分析的材料是在西西伯利亚获得的。 众所周知，该地区是我国的主要石油产区。 该地区上世纪 60 年代开始积极开发矿藏。 寻找石油矿床的主要方法是地震勘探。 看看这个地区的卫星图像很有趣。 小比例尺下，可以看到数量巨大的沼泽、湖泊；放大地图，可以看到丛式钻井场地，将地图放大到极限，还可以分辨出地震剖面上的空地。进行了观察。

Yandex 地图卫星图像 - 诺亚布尔斯克市区


其中一个油田的井场网络


西西伯利亚的含油岩石的深度范围很广——从 1 公里到 5 公里。 含油岩石的主要体积形成于侏罗纪和白垩纪。 侏罗纪时期可能是许多人通过同名电影而了解的。 侏罗纪气候 与现代有很大不同。 大英百科全书有一系列古地图，描绘了每个神学时代的特征。

目前，

侏罗纪时期


请注意，在侏罗纪时期，西西伯利亚的领土是海岸（河流和浅海穿过的土地）。 由于气候舒适，我们可以假设当时的典型景观是这样的：

侏罗纪西伯利亚


在这张照片中，对我们来说重要的不是动物和鸟类，而是背景中河流的图像。 这条河与我们之前停过的地质物体是同一个。 事实上，河流的活动使得分选良好的砂岩积累起来，然后成为石油储层。 这些水库可能具有奇异、复杂的形状（如河床），并且厚度可变 - 靠近河岸的厚度较小，但靠近河道中心或在蜿蜒区域的厚度会增加。 因此，侏罗纪形成的河流现在深度约为三公里，是寻找油藏的目标。

实验数据。 处理和可视化

让我们立即对文章中显示的地震材料进行保留 - 由于用于分析的数据量很大 - 文章文本中仅包含原始地震道集的片段。 这将允许任何人重现上述计算。

在处理地震数据时，地球物理学家通常使用专门的软件（有几个行业领导者的开发成果被积极使用，例如 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]

通过这种方式下载和导入的数据是在10×5公里的区域内记录的路线。 如果使用三维地震勘测方法获得数据（不是沿着单独的地球物理剖面记录波，而是在整个区域同时记录波），则可以获得地震数据立方体。 这些是三维物体，其垂直和水平截面允许对地质环境进行详细研究。 在所考虑的示例中，我们正在处理三维数据。 我们可以从文本标题中获取一些信息，就像这样

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

C 1 这是地质包测试的演示文件
Ç2的
Ç3的
Ç4的
C 5 日期用户名：WOLFRAM 用户
C 6 调查名称：西伯利亚某处
C 7 文件类型 3D 地震体积
Ç8的
Ç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]

地震剖面轨迹之一


知道研究区域的每个方向上有多少条迹线，您可以生成三维数据数组并使用 Image3D[] 函数将其显示

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 语言使用该函数进行小波变换 连续小波变换[]。 并且将该函数应用于整组迹线将使用该函数进行 桌子[]。 这里值得注意的是 Wolfram Mathematica 的优势之一 - 使用并行化的能力 并行表[]。 在上面的例子中，不需要并行化——数据量并不大，但是当处理包含数十万条迹线的实验数据集时，这是必要的。

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

应用该功能后 连续小波变换[] 将出现与所选频率相对应的新数据集。 在上面的示例中，这些频率是：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];

对于其中三个立方体（低频、中频和高频分量），RGB 混合通常用于一起可视化数据。 每个组件都分配有自己的颜色 - 红色、绿色、蓝色。 在 Wolfram Mathematica 中，这可以使用函数来完成 颜色组合[].

结果是可以进行地质解释的图像。 该剖面上记录的曲流使得能够描绘古河道，这些古河道更有可能是储层并含有石油储量。 对这种河流系统的现代类似物的搜索和分析使我们能够确定曲流中最有希望的部分。 这些河道本身的特点是厚层的分选良好的砂岩，是良好的石油储层。 “花边”异常区域之外的区域与现代洪泛区沉积物相似。 洪泛区沉积物主要以粘土岩为代表，在这些区域钻探是无效的。

数据立方体的 RGB 切片。 在中心（中心稍偏左），您可以追踪蜿蜒的河流。

数据立方体的 RGB 切片。 左边可以看到蜿蜒的河流。


在某些情况下，地震数据的质量可以使图像变得更加清晰。 这取决于现场工作方法、降噪算法所使用的设备。 在这种情况下，不仅可以看到河流系统的碎片，还可以看到整个延伸的古河流。

地震数据立方体（水平切片）的三个分量的 RGB 混合。 深度约2公里。

萨拉托夫附近伏尔加河的卫星图像

结论

Wolfram Mathematica 允许您分析地震数据并解决与矿产勘探相关的应用问题，而 GeologyIO 软件包使此过程更加方便。 地震数据的结构是这样的，使用内置方法来加速计算（并行表[], 并行[],...) 非常高效，允许您处理大量数据。 在很大程度上，GeologyIO 包的数据存储功能促进了这一点。 顺便说一句，该软件包不仅可以用于应用地震勘探领域。 探地雷达和地震学中使用的数据类型几乎相同。如果您对如何改进结果有建议，Wolfram Mathematica 武器库中的哪些信号分析算法适用于此类数据，或者您有任何批评意见，请发表评论。

来源： habr.com