使用Python数据可视化库可视化地下
在地球科学中，了解地下地质表面至关重要。通过知道地层的确切位置和几何形状，我们可以更容易地识别潜在的石油和天然气勘探目标，以及潜在的碳捕获和储存位置。我们可以使用各种方法来优化这些表面，从地震数据到井记录派生的地层顶部。通常，这些技术会结合使用以优化最终表面。
本文是我之前文章的延续，该文章着重介绍了将井录测量值外推到区域以了解和可视化地理空间变化的过程。在本文中，我们将看看如何使用交互式的Plotly图表创建3D表面。
由于建模地质表面是一个复杂的过程，通常涉及多次迭代和优化，因此本文演示了如何使用Python非常简单的示例来可视化这些数据。
为了了解如何使用Plotly在区域内可视化地质层顶，我们将使用两组数据。
第一组数据源自于与地层相交的28个井眼派生的数据，这些数据用作输入进行克里格插值以产生低分辨率表面。相比之下，第二组数据源自于解释的地震数据，提供了更高分辨率的表面。
这两组数据都来自Equinor Volve数据集。有关详细信息，请参见本文底部。
您也可以在下面的链接中查看本系列中的以下文章：
- Plotly和Python：创建用于岩石物理和地质数据的交互式热度图
- 使用pykrige和matplotlib对地质变化进行空间可视化
- 使用Plotly Express在3D线图上可视化井道
导入库和数据
在处理数据之前，我们首先需要导入所需的库。这些是：
- pandas – 读取我们的数据，数据格式为
csv - matplotlib – 创建我们的可视化
- pykrige – 执行克里格插值
- numpy – 用于一些数值计算
- plotly.graph_objects – 以3D形式可视化表面
import pandas as pdimport matplotlib.pyplot as pltimport plotly.graph_objects as gofrom pykrige import OrdinaryKrigingimport numpy as np接下来，我们可以使用pd.read_csv()加载数据。
由于这些数据包含有关Volve油田中所有井的地质表面信息，因此我们可以使用query()提取我们需要的地层数据。在本例中，我们将查看Hugin Formation。
df = pd.read_csv('Data/Volve/Well_picks_Volve_v1 copy.csv')df_hugin = df.query('SURFACE == "Hugin Fm. VOLVE Top"')df_hugin运行上面的代码后，我们会得到以下表格。您可能会注意到，一些井眼多次遇到了Hugin Formation。这可能是由于井眼或地层几何形状多次穿过地层所致。
利用TVDSS外推生成地质表面
在我的以前的文章中，我专注于如何使用一种称为克里金的过程来“填补”测量点之间的空白。在本文中，我们不会涵盖这个过程的细节；但是，您可以查阅这篇文章获取更多信息。
一旦我们的数据被加载，我们可以通过调用pykrige的OrdinaryKriging方法来运行kriging过程。
在这个调用中，我们传入我们的x和y数据，它们表示井眼在地下岩层中遇到的地点的东向和北向位置。
由于我们想要生成Hugin岩层的表面，我们需要使用TVDSS – True Vertical Depth Subsea – 测量值。这给出了表面在所选数据下面的真实深度。
OK = OrdinaryKriging(x=df_hugin['Easting'], y=df_hugin['Northing'], z=df_hugin['TVDSS'], variogram_model='linear', verbose=True, enable_plotting=True)
一旦模型生成完成，我们可以将其应用于覆盖整个井/穿透点范围的两个数组。
这使我们能够生成一组值的网格，然后将其传递到我们生成的OrdinaryKriging对象中。
gridx = np.arange(433986, 438607, 50, dtype='float64')gridy = np.arange(6477539, 6479393, 50,dtype='float64')zstar, ss = OK.execute('grid', gridx, gridy)最后，我们可以使用matplotlib的imshow方法生成简单的2D地图视图。
fig, ax = plt.subplots(figsize=(15,5))# Create a 2D image plot of the data in 'zstar'# The 'extent' parameter sets the bounds of the image in data coordinates# 'origin' parameter sets the part of the image that corresponds to the origin of the axesimage = ax.imshow(zstar, extent=(433986, 438607, 6477539, 6479393), origin='lower')# Set the labels for the x-axis and y-axisax.set_xlabel('X Location (m)', fontsize=14, fontweight='bold')ax.set_ylabel('Y Location (m)', fontsize=14, fontweight='bold')# Add contourscontours = ax.contour(gridx, gridy, zstar, colors='black')colorbar = fig.colorbar(image)colorbar.set_label('DTC (us/ft)', fontsize=14, fontweight='bold')# Display the plotplt.show()
使用Plotly创建简单的3D表面图
要将我们的2D表面转换为3D，我们需要使用Plotly。我们可以使用matplotlib来做到这一点；但是，从我的经验来看，使用Plotly生成3D可视化更容易、更流畅、更交互。
在下面的代码中，我们首先创建我们的坐标网格。为此，我们使用numpy的linspace函数。该函数将在指定范围内创建一组均匀间隔的数字。对于我们的数据集和示例，范围从xgrid_extent和ygrid_extent的最小值到最大值。
在此范围内使用的值的总数将等于我们在上面创建的zstar网格中存在的x和y元素的数量。
形成网格后，我们随后调用Plotly。
首先，我们创建我们的图形对象，然后使用fig.add_trace将我们的3D表面绘图添加到图形中。
添加完毕后，我们需要调整绘图的布局，使其具有轴标签、定义的宽度和高度以及一些填充。
xgrid_extent = [433986, 438607]ygrid_extent = [6477539, 6479393]x = np.linspace(xgrid_extent[0], xgrid_extent[1], zstar.shape[1])y = np.linspace(ygrid_extent[0], ygrid_extent[1], zstar.shape[0])fig = go.Figure()fig.add_trace(go.Surface(z=zstar, x=x, y=y))fig.update_layout(scene = dict( xaxis_title='X位置', yaxis_title='Y位置', zaxis_title='深度'), width=1000, height=800, margin=dict(r=20, l=10, b=10, t=10))fig.show()运行上面的代码后，我们获得以下交互式绘图，显示了基于钻井井眼的多次遇到的Hugin地层的地质表面。
使用Plotly查看完全解释的表面
Volve数据集具有许多从地质解释（包括地震数据）生成的完全解释的表面。
这些数据包含场中数据点的x和y位置，以及我们的TVDSS数据（z）。
Volve数据门户网站上提供的数据以.dat文件的形式提供，但是，通过在文本编辑器中进行一些操作，可以轻松地将其转换为CSV文件并使用pandas加载。
hugin_formation_surface = pd.read_csv('Data/Volve/Hugin_Fm_Top+ST10010ZC11_Near_190314_adj2_2760_EasyDC+STAT+DEPTH.csv')
加载数据后，我们可以通过提取x、y和z数据到变量中使事情变得更容易。
x = hugin_formation_surface['x']y = hugin_formation_surface['y']z = hugin_formation_surface['z']然后，我们需要在x和y数据位置内的最小和最大位置之间创建一个定期间隔的网格。这可以使用numpy的meshgrid完成。
xi = np.linspace(x.min(), x.max(), 100)yi = np.linspace(y.min(), y.max(), 100)xi, yi = np.meshgrid(xi, yi)有几种方法可以在一系列数据点之间进行插值。所选择的方法将取决于数据的形式（定期采样的数据点与非定期采样的数据点），数据大小和计算能力。
如果我们有像这里一样的大数据集，则使用一些方法（例如径向基函数）将更加计算昂贵。
在此示例中，我们将使用scipy中的LinearNDInterpolator方法构建我们的模型，然后将其应用于我们的z（TVDSS）变量。
为了使我们在点之间进行插值，我们需要将xi，yi重塑为用于插值的1D数组，因为LinearNDInterpolator期望1D数组。
xir = xi.ravel()yir = yi.ravel()interp = LinearNDInterpolator((x, y), z)zi = interp(xir, yir)计算完成后，我们可以使用 Plotly Graph Objects 创建我们的 3D 表面。
fig = go.Figure()fig.add_trace(go.Surface(z=zi, x=xi, y=yi, colorscale='Viridis'))fig.update_layout(scene = dict( xaxis_title='东向距离（米）', yaxis_title='北向距离（米）', zaxis_title='深度', zaxis=dict(autorange='reversed')), width=1000, height=800, margin=dict(r=20, l=10, b=10, t=10))fig.update_traces(contours_z=dict(show=True, usecolormap=True, project_z=True, highlightcolor="white"))fig.show()运行上述代码后，我们得到了以下 Hugin Formation 的 3D 表面图。
当我们将此图与由井孔测量生成的图进行比较时，我们可以明显看到两者的整体形状相似，中间有山谷。但是，由地震数据生成的表面提供了比由井孔测量生成的地层顶面更详细的信息。
总结
在这个简短的教程中，我们看到了如何使用 Plotly 的 3D 表面图生成地质表面的交互式 3D 可视化。使用井孔测量导出的地层顶面，我们可以生成一个非常基本的 3D 表面。这是因为测量受限于已经穿过 Hugin Formation 的井孔，这意味着我们得到的表面分辨率很低。
相比之下，如果我们有更详细的测量点，如地震导出的层位面，我们可以生成一个更加真实的地质表面。
两种方法都是有效的，但是，您必须记住，当从仅由井孔测量导出的地层顶部进行外推时，我们可能无法在整个区域内生成该地层的全面图像。
使用的数据集
本教程使用的数据是 Equinor 在 2018 年发布的 Volve 数据集的子集。数据集的完整详情（包括许可证）可以在下面的链接中找到：
Volve 油田数据集
Equinor 已经发布了 2008-2016 年 Volve 油田的完整数据集。点击此处下载进行研究、研究…
www.equinor.com
Volve 数据许可证基于 CC BY 4.0 许可证。完整的许可协议详情可以在以下链接中找到：
https://cdn.sanity.io/files/h61q9gi9/global/de6532f6134b9a953f6c41bac47a0c055a3712d3.pdf?equinor-hrs-terms-and-conditions-for-licence-to-data-volve.pdf
感谢阅读。在离开之前，您应该订阅我的内容，并在您的收件箱中获取我的文章。 您可以在这里完成订阅！
其次，您可以通过注册会员，获得完整的小猪AI体验，并支持成千上万的其他作家和我。每月仅需5美元，您就可以完全访问所有精彩的小猪AI文章，并有机会通过写作赚钱。
如果您使用我的链接注册，您将直接支持我一部分的费用，而不会增加您的费用。如果您这样做了，非常感谢您的支持。