python也可以使用克里金插值算法吗？

     

挪威大陆架的声学压缩慢度测量的空间变化

在处理地质和岩石物理数据时，我们通常希望了解这些数据在我们的地区是如何变化的。我们可以做到这一点的方法之一是对我们的实际测量值进行网格化，并推断这些值。

进行这种外推的一种特殊方法是克里金法，这是一种以南非采矿工程师 Danie G. Krige 命名的地质统计程序。克里金法背后的思想在于其估计技术：它使用观测数据之间的空间相关性来预测未测量位置的值。

通过衡量变量如何随距离变化，该方法建立了一种统计关系，可用于预测整个区域的值，将分散的数据点转换为连贯的空间地图。

在本教程中，我们将了解一个名为pykrige的 Python 库。该库专为 2D 和 3D 克里金法计算而设计，易于与数据一起使用。

导入库和数据

首先，我们需要导入我们需要的库。对于本文，我们将需要以下库：

• pandas — 读取我们的数据，csv格式

• matplotlib来创建我们的可视化

• pykrige进行克里金法

• numpy用于一些数值计算

import pandas as pd          
import matplotlib.pyplot as plt          
         
from pykrige import OrdinaryKriging          
import numpy as np

导入库后，我们现在可以导入数据。

在本教程中，我们将使用来自 Xeek 和 Force 2020 机器学习竞赛的数据集来根据测井数据预测岩性。竞赛数据集的这个子集包含 65 个井位，其中包含 Balder 地层的平均声学纵波慢度测量值。要读取我们的数据，我们可以使用 pandasread_csv()函数，并传入数据文件的位置。

df=pd.read_csv('Data/Xeek/Force2020/Xeek_2020_Baldr_DTC_AVG.csv')          

当我们查看数据时，我们会看到我们有 65 个井，其中包含 Balder Formation 顶部的位置（X_LOC 和 Y_LOC 用于网格坐标，LAT 和 LON 用于纬度和经度）。我们还有遇到地层的真实海底垂直深度 (TVDSS)，以及声学纵波时差 (DTC) 的平均值。

可视化井的空间位置

现在我们的数据已成功加载到数据框中，我们可以可视化我们的数据以了解我们的井的位置。为此，我们将使用 matplotlib 的散点图并传入经度和纬度列。

plt.scatter(df['Longitude'],df['Latitude'], c=df['DTC'])

当我们运行上面的代码时，我们得到下面的图。

我们可以看到上图非常基本，没有颜色条或轴标签。

让我们通过向其中添加这些功能来稍微修改绘图。

cm = plt.cm.get_cmap('viridis')          
         
plt.figure(figsize=(10,10))          
scatter = plt.scatter(df['LON'], df['LAT'], c=df['DTC_MEAN'], cmap=cm, s=50)          
         
plt.colorbar(scatter)          
         
plt.xlabel('Longitude')          
plt.ylabel('Latitude')          
plt.show()

当我们运行上面的代码时，我们得到下图，它告诉我们更多关于我们的数据。我们可以使用颜色条来估计我们的点值。

应用克里金法

为了更好地了解我们的数据点以及 Balder Formation 整个区域的 DTC 测量值如何变化，我们可以使用克里金法和我们的数据点来填补测量值之间的差距。

为此，我们需要OrdinaryKriging从 pykrige 库创建一个对象。

我们将 x 和 y 的位置数据以及我们要映射到 z 参数的数据传递到该对象中。

我们还需要选择我们想要使用的变差函数模型。在这种情况下，我们将使用指数模型。可以在文档中找到有关模型类型的更多详细信息。

由于我们使用纬度和经度作为 x 和 y 坐标，我们可以将 coordinates_type 参数更改为geographic

OK = OrdinaryKriging(x=df['LON'],          
                      y=df['LAT'],          
                      z=df['DTC_MEAN'],          
                      variogram_model='exponential',          
                      verbose=True, enable_plotting=True,          
                      coordinates_type='geographic')

当我们运行上面的代码时，我们返回以下模型摘要和半变异函数。

以下是返回参数的简要说明：

• 块金：块金是变差函数的 y 轴截距，表示零距离处的方差，通常是由于测量误差或非常小的尺度变化引起的。

• Full Sill：sill 是变差函数达到并开始趋于平稳的最大方差，当点相距很远时会发生这种情况。

• 范围：范围是变差函数到达基台的距离，意思是超过该距离进一步分离点不会增加方差。

• 偏基台：偏基台是基台和块块之间的差异，表示数据中空间结构化的方差量。

这可以让我们根据生成的线和点的形状了解我们的模型对数据的适用程度。

显示克里格结果

要开始显示我们的数据，我们需要创建一个数据网格。

为此，我们首先为我们定义的坐标之间的纬度和经度创建数组。在这种情况下，我们希望图表从 57.5 度 N 扩展到 62 度 N 以及从 1.5 度 E 到 4.5 度 E。

使用np.arange将允许我们以规则的间距创建这些阵列。

grid_lat = np.arange(57.5, 62, 0.01, dtype='float64')          
grid_long = np.arange(1.5, 4.5, 0.01,dtype='float64')

现在我们有了 X 和 Y 坐标，我们可以创建我们的值网格。为此，我们调用OK.execute，并传入我们的纬度和经度数组。

zstar, ss = OK.execute('grid', grid_long, grid_lat)

这将返回两个数组。我们的数据网格 (zstar) 和与之相关的不确定性 (ss)

接下来，我们现在可以使用我们的数据数组并使用 matplotlib 的imshow.

import matplotlib.pyplot as plt          
         
fig, ax = plt.subplots(figsize=(10,10))          
         
image = ax.imshow(zstar, extent=(1.5, 4.5, 57.5, 62), origin='lower')          
         
ax.set_xlabel('Longitude', fontsize=14, fontweight='bold')          
ax.set_ylabel('Latitude', fontsize=14, fontweight='bold')          
         
scatter = ax.scatter(x=df['LON'], y=df['LAT'], color='black')          
         
colorbar = fig.colorbar(image)          
         
colorbar.set_label('DTC (us/ft)', fontsize=14, fontweight='bold')          
         
plt.show()

当我们运行它时，我们得到以下地图，显示了我们 65 口井中 Balder 地层的声学纵波慢度变化。

使用 pykrige 生成的声学压缩慢度 (DTC) 的数据网格。

我们可以看到，在北纬 59 到 60 度附近，岩石速度要快得多，而在东北和西南部地区，岩石要慢得多。

要解释这一点，我们需要了解每口井的地层深度。这将使我们能够确定差异是否与埋藏和压实或其他地质过程有关。

我们将在以后的文章中看到如何做到这一点。

可视化克里格不确定性

查看此类数据的关键之一是了解与克里金法相关的不确定性。

我们可以通过重新调用相同的绘图代码来做到这一点，而不是zstar传入，我们可以将它交换为ss我们之前创建的变量。

fig, ax = plt.subplots(figsize=(10,10))          
         
image = ax.imshow(ss, extent=(1.5, 4.5, 57.5, 62), origin='lower')          
         
ax.set_xlabel('Longitude', fontsize=14, fontweight='bold')          
ax.set_ylabel('Latitude', fontsize=14, fontweight='bold')          
         
scatter = ax.scatter(x=df['LON'], y=df['LAT'], color='black')          
         
colorbar = fig.colorbar(image)          
         
colorbar.set_label('DTC (us/ft)', fontsize=14, fontweight='bold')          
         
plt.show()

通过下图，我们可以看到不确定性高或低的区域。

使用 pykrige 生成的声学压缩慢度 (DTC) 的不确定性数据网格。

在我们的井覆盖较少的地区，我们的不确定性会高得多，而在我们有多口井的地区，我们的不确定性会低得多。

概括

在本教程中，我们了解了如何获取测井测量 (DTC) 的平均值并将它们映射到整个区域。这使我们能够了解某个地理区域的数据趋势。

然而，在查看这些数据时，我们必须记住，我们正在查看的是 2D 表面，而不是我们在地下遇到的更复杂的 3D 结构。因此，测量的变化可归因于深度的变化。

使用的数据集

本文中使用的数据集是训练数据集的一个子集，该数据集用作 Xeek 和 FORCE 2020 （Bormann 等人，2020）举办的机器学习竞赛的一部分。它是根据挪威政府的 NOLD 2.0 许可证发布的，有关详细信息，请参见：挪威开放政府数据许可证 (NLOD) 2.0。可以在此处（https://data.norge.no/nlod/en/2.0）访问完整的数据集。