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

news2024/11/23 10:06:48

     41496ea2caa4663d7832a93e7b04568a.png

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

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

进行这种外推的一种特殊方法是克里金法,这是一种以南非采矿工程师 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) 的平均值。

ca705a8b487f9007b5778d3e43ea3c80.png

可视化井的空间位置

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

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

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

6bde3006dc8598a88763d52511ca5959.png

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

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

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()

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

8c4a09f78823caafc76438abd64cb6e0.png

应用克里金法

为了更好地了解我们的数据点以及 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')

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

4cee008d1dec1f9982c0590a3bd4df08.png

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

  • 块金:块金是变差函数的 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) 的数据网格。8d8e469448e456625cc4b6cc9a34cb30.png

我们可以看到,在北纬 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()

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

2e737be8dcb2506a02083fdbf7ba4dab.png

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

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

概括

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

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

使用的数据集

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

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.coloradmin.cn/o/677820.html

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈,一经查实,立即删除!

相关文章

三阶魔方有多少种状态

魔方有 3 种不同的方块,分别为角块(8 个,每个角块有三种颜色),棱块(12 个,每个棱块有两种颜色)与中心块(6 个,每个中心块有一种颜色)。 魔方总共…

每天学一点知识有用吗

在探索如何学习的路上,我注意到了基于微习惯的学习方式,比如每天在用十分钟的时间练习下普通话,或者每天写500字的总结。 我简单回顾一下: 这种方法虽然颇受欢迎,但是它限制了你可以尝试的活动种类,有时候…

深度学习(24)——YOLO系列(4)

深度学习(24)——YOLO系列(4) 文章目录 深度学习(24)——YOLO系列(4)1. dataset准备(1)数据详解(2)dataset(3)…

广告数仓:全流程调度

系列文章目录 广告数仓:采集通道创建 广告数仓:数仓搭建 广告数仓:数仓搭建(二) 广告数仓:全流程调度 文章目录 系列文章目录前言一、ClickHouse安装1.修改环境2.安装依赖3.单机安装4.修改配置文件5.启动clickhouse6.创建需要的数…

012-从零搭建微服务-接口文档(二)

写在最前 如果这个项目让你有所收获,记得 Star 关注哦,这对我是非常不错的鼓励与支持。 源码地址(后端):https://gitee.com/csps/mingyue 源码地址(前端):https://gitee.com/csps…

统一拦截--过滤器Filter

1.过滤器Filter 1. 概述 概念: Filter过滤器,是JavaWeb三大组件(Servlet、Filter、Listener)之一。过滤器可以把对资源的请求拦截下来,从而实现一些特殊的功能。过滤器一般完成一些通用的操作,比如:登录校验、统一编码处理、敏感字符处理等…

Tcp协议的十大特性详解+示例

前言 之前我们简单了解了一下Tcp是什么及它的套接字如何使用:基于UDP和TCP套接字实现简单的回显客户端服务器程序_Crystal_bit的博客-CSDN博客 因为要给大家介绍Tcp的十大特性,所以这里给出Tcp报头结构: 目录 1. 确认应答 2. 超时重传 3. 连接管理 3…

【Android复习笔记】Parcelable 为什么速度优于 Serializable ?

Q:Parcelable 为什么速度优于 Serializable ? 首先,抛开应用场景谈技术方案都是在耍流氓,所以如果你遇到有面试官问这样的题目本身就是在给面试者挖坑。 序列化 将实例的状态转换为可以存储或传输的形式的过程。 Serializable 实现方式: Serializable 是属于 Java 自带的…

Solid Converter PDF v10 安装及使用教程

目录 一、软件介绍二、下载教程三、安装教程四、使用教程1.PDF转Word、Html等2.合并PDF文件 一、软件介绍 Solid Converter PDF是一套专门将PDF文件转换成Word的软件。 能够将PDF转换为Word、Excel、HTML、PowerPoint、纯文本文件从PDF文档中提取数据并以CSV等格式保存能够转…

数仓工程师理解复杂业务的思考方法论

模型设计框架(业务过程驱动)还是在经典的三层数据模型架构下去进行,概念模型、逻辑模型、物理模型 首先概念模型其实是业务过程(流程图),其中需要考虑到几个方面: 1.数据 业务覆盖 业务感知、…

循坏队列CircularQueue

前言 一、CircularQueue 二、特点 三、设计思路 1)判空与判满 2)链表还是数组实现? 四、实现 1).IsEmpty() 2).IsFull() 3)CircularQueueCreate创建 4)CircularQueueEnQueue插入 5)CircularQueueDeQueue删除 6&#xf…

React Hook之useCallback 性能优化

上文 对比之前的组件优化说明React.memo的作用我们说了 React.memo的妙用 但是 它却并非万能 我们来看这个情况 我们子组件代码编写如下 import React from "react";const ChildComponent ({ dom1funt }) > {console.log("ChildComponent 被重新渲染"…

规则引擎--规则集:规则集合的组织和执行

目录 回顾easy-rules的rules执行如何想规则集合的构造 规则集合定义普通规则集和执行定义树形规则集 当弄清楚了一个规则的设计和执行逻辑后,接下来需要考虑的就是许多的规则如何组织了,即规则集的抽象设计。 来看一些例子 回顾easy-rules的rules执行 …

NFCEE Discovery and Mode Set

10.1 NFCEE ID NFCC 动态为 NFCEE 分配 ID(称为“NFCEE ID”)。 DH 通过执行 NFCEE Discovery 来了解 ID 值。 在配置状态为 0x01 的 NFCC 重置之前,NFCEE ID 一直有效。 值为 0x00 的 ID 在本规范中称为 DH-NFCEE ID,并且应代表…

五、Docker本地镜像发布到阿里云/发布到私有库

目录 前言一、本地镜像发布到阿里云1.1 流程图1.2 注册阿里云创建容器服务个人实例1.3 创建命名空间1.4 创建镜像仓库1.5 将镜像推送到阿里云本地仓库 二、从阿里云仓库拉去自己推送的镜像三、本地镜像发布到阿里云总结四、本地镜像发布到私有库4.1 流程图4.2 下载镜像Docker R…

Shell编程从入门到实践——实践篇

欢迎关注 「Android茶话会」 回 「学习之路」 取Android技术路线经典电子书回 「pdf」 取阿里&字节经典面试题、Android、算法、Java等系列武功秘籍。回 「天涯」 取天涯论坛200精彩博文,包括小说、玄学等 背景 之前在搞一些CI/CD,使用到了shell脚本,shell的开…

nvdiffrec在Windows上的配置及使用

nvdiffrec是NVIDIA研究院开源的项目,源代码地址:https://github.com/NVlabs/nvdiffrec ,论文为《Extracting Triangular 3D Models, Materials, and Lighting From Images》,从图像中提取三角形三维(三角网格)模型、空间变化的材质…

uni-app微信小程序获取手机号授权登录(复制即用,js完成敏感数据对称解密,无需走服务端处理)

目录 一、示例 二、具体实现说明 一、示例 获取到的手机号 二、具体实现说明 属性说明 属性名说明生效时机getphonenumber获取用户手机号回调open-type"getPhoneNumber" 按钮写法 <template><view class"login"><view class"content…

为什么要写这个带点玄幻气息的英语单词记忆博客

&#x1f31f;博主&#xff1a;命运之光 ☀️专栏&#xff1a;英之剑法&#x1f5e1; ❤️‍&#x1f525;专栏&#xff1a;英之试炼&#x1f525; ☀️博主的其他文章&#xff1a;点击进入博主的主页 &#x1f433; 开篇想说的话&#xff1a;开学就大三了&#xff0c;命运之光…

DMA详解及应用(嵌入式学习)

DMA 0. 前言1. DMA作用2. DMA特性3. DMA寄存器4. DMA的增量或者循环模式5. 练习 0. 前言 DMA&#xff08;Direct Memory Access&#xff0c;直接内存访问&#xff09;是一种计算机系统中用于高效地实现数据传输的技术。它允许数据在外设和内存之间直接传输&#xff0c;而无需C…