BimAnt

探索空间数据的插值

admin -

本文将介绍空间数据的插值算法，包括最近邻算法、高斯过程、克里金插值等。

0、序言

在之前的文章中，我解释了如何使用 plotly Dash 创建一个 GIS webapp，这是一个用于构建 web 分析应用程序的高效 Python 框架。 在这篇文章中，我想展示我们如何使用 dash-vtk 来渲染 3D 数据。 例如，你可以使用 dash-vtk 渲染表示 LiDAR 数据集或地形跟随网格的网格。 此处的目的是比较用于插值日内瓦（瑞士）地区地质单元（莫拉斯）深度的不同算法。

有一个与本教程相关联的 GitHub 存储库，其中包含用于空间分析的笔记本。该应用程序托管在 GitHub 上，并部署在 Heroku 上。 请耐心等待，唤醒应用程序 Heroku dynos 🤞🤞 需要几秒钟。

1、简介

本教程中使用的数据集来自日内瓦州地质调查局 (GESDEC) 管理的井地登记册，它是 SITG (System d’Information du Territoire à Genève) 开放数据目录的一部分。 可以在这里下载。 从这个数据集中，我们想要使用以下插值算法在日内瓦的一个子区域内插值 Molasse 的深度：

  • 最近邻
  • 双谐波样条（使用verde）
  • 高斯过程
  • 克里金法
  • 顺序高斯模拟

2、数据集下载

数据集可以在这里下载或使用 opendata.swiss API 获取，如下所示：

from zipfile import ZipFile
import urllib.request
from io import BytesIO
data_url = ‘https://ge.ch/sitg/geodata/SITG/OPENDATA/4108/CSV_GOL_SONDAGE.zip'
zip_file = ZipFile(BytesIO(urllib.request.urlopen(data_url).read()))
dfs = {text_file.filename: pd.read_csv(zip_file.open(text_file.filename), encoding=’ISO-8859–1', 
sep=’;’,)
 for text_file in zip_file.infolist()
 if text_file.filename.endswith(‘.csv’)}
sondage = dfs[‘GOL_SONDAGE.csv’]

3、插值算法代码

有关本文后续介绍的空间分析的详细说明，请参阅这些笔记。

4、最近邻算法

应用 NN 算法的最简单方法是使用 Verde 开发包 (Uieda L, 2018)¹。 数据集已被处理为 50m x 50m 的规则网格。

import verde as vd
hain = vd.Chain(
    [
        ("mean", vd.BlockReduce(np.mean, spacing=spacing )),
        ("spline", vd.ScipyGridder(method="nearest")),
    ]
)

train, test = vd.train_test_split(
    coordinates,df.MOLASSE, random_state=0
)

chain.fit(*train)

grid = chain.grid(
    region=region,
    spacing=spacing,
    dims=["N", "E"],
    data_names="MOLASSE",
)

5、双谐波样条

接下来，使用已经过测试格林函数的双谐波样条插值，其中两个主要参数是 mindist （点作用力与数据点的最小距离）以及 damping（控制对估计力施加的平滑度）。

chain = vd.Chain(
    [
        ("mean", vd.BlockReduce(np.mean, spacing=spacing )),
        ("spline", vd.Spline()),
    ]
)

chain.fit(coordinates, df.MOLASSE)

grid = chain.grid(
    region=region,
    spacing=spacing,
    # projection=projection,
    dims=["N", "E"],
    data_names="MOLASSE",
)

6、双谐波样条（交叉验证）

verde.SplineCV 类提供了一个交叉验证的 verde.Spline 版本。 它具有几乎相同的接口，但在拟合数据集时会自动完成上述所有操作。 唯一的区别是你必须提供要尝试的 damping和 mindist参数列表，而不是仅提供单个值：

dampings = [None, 1e-4, 1e-3, 1e-2]
mindists = [1e-5, 1, 100, 1000, 1000]

spline = vd.SplineCV(
    dampings=dampings,
    mindists=mindists,
)
spline.fit(coordinates, df.MOLASSE)
print("Highest score:", spline.scores_.max())
print("Best damping:", spline.damping_)
print("Best mindist:", spline.mindist_)

调用 fit 将对所有参数组合进行网格搜索，以找到最大化交叉验证分数的参数组合。

chain = vd.Chain(
    [
        ("mean", vd.BlockReduce(np.mean, spacing=spacing )),
        ("spline", vd.Spline(damping=spline.damping_, mindist=spline.mindist_)),

    ]
)

chain.fit(coordinates,df.MOLASSE)

grid = chain.grid(
    region=region,
    spacing=spacing,
    # projection=projection,
    dims=["N", "E"],
    data_names="MOLASSE",
)

请注意，对于像这样的稀疏数据，更平滑的模型往往是更好的预测器。 这表明你可能不应该相信我们从默认设置中获得的许多短波长特性。

7、高斯过程

高斯过程 (GP) 是一种通用的监督学习方法，旨在解决回归和概率分类问题。

高斯过程的优点是：

  • 预测对观察结果进行插值（至少对于常规内核）。
  • 预测是概率性的（高斯），因此可以计算经验置信区间并根据这些决定是否应该重新拟合（在线拟合、自适应拟合）某个感兴趣区域的预测。
  • 通用性：可以指定不同的内核。 提供了通用内核，但也可以指定自定义内核。
from sklearn.gaussian_process.kernels import RBF
from sklearn.gaussian_process import GaussianProcessRegressor as GPR
from sklearn.gaussian_process.kernels import RBF, ConstantKernel as C

kernel = C(50.0) * RBF([50 ,50])


gp = GPR(normalize_y=True, alpha=0.5, kernel=kernel,)

gp.fit(df[['E', 'N']].values, df.MOLASSE.values)

8、克里金法

克里金法（Kriging）是一种广泛使用的空间插值器。这里我们使用 GeostatPy 包 (Pyrcz et al., 2021)²，它是从原始 FORTRAN GSLIB: Geostatistical Library methods (Deutsch and Journel, 1997)³ 翻译的 Python 版本。以下是克里金法的一些重要性质（Pyrcz，2021）：

  • 精确插值器 - 使用位置的数据值进行克里金估计
  • 可以在获取样本信息之前计算克里金方差，因为克里金估计方差不依赖于数据的值，也不依赖于克里金估计，即克里金估计量是同方差的。
  • 空间背景——克里金法考虑到了空间连续性、接近性和冗余性的陈述，我们可以说克里金法解释了数据的配置和被估计变量的结构连续性。
  • 规模——克里金法可以推广到考虑数据和估计的支持量。我们稍后会介绍。
  • 多变量——克里金法可以推广到使用协同克里金法系统来解释空间估计中的多个辅助数据。我们稍后会介绍。
  • 可以预测克里金法的平滑效果。稍后我们将使用它来构建随机模拟。

获得克里金估计的步骤稍微复杂一些，首先应该计算一个变异函数（Variogram）。如需深入了解变异函数计算和建模，请查看1和2。

9、顺序高斯模拟

SGS （Sequential Gaussian Simulation）的目标是建立尊重感兴趣属性的单变量和空间分布的空间模型。 至于克里金法，可以在此处、此处和此处找到 SGS 的一些优秀资源。

对于本教程，已经生成了10 个可能的实现：

与克里金法和确定性插值器相比，地理统计模拟的一个优势是可以总结实现的不确定性，这对以下方面有用：

  • 量化远离油井的局部不确定性
  • 获取特定本地结果的概率/风险

我们特别关注以下局部摘要：

  • e-type 是局部期望值（因为等权重与平均值相同）
  • 条件标准偏差是实现的局部标准偏差
  • 局部百分位数是与实际相比的局部百分位数

10、E型和条件方差

我们将从 e 类型和条件方差开始。

  • e-type 是局部期望（只是位置 u(α​) 的 L 个实现的平均值，因为我们假设所有实现的可能性相同）。
  • 条件方差是局部方差

e 型模型与克里金模型非常相似，除了：

  • 高斯正向和反向变换可能会改变结果
  • 由于实现太少，结果很嘈杂

条件方差在数据位置最低，远离数据增加

  • 由于实现太少，结果很嘈杂

11、局部百分位数

现在让我们看一下，局部百分位图是具有从局部实现中采样的局部百分位值的地图

我们可以将它们解释如下，在一个位置，如果我们的局部 P10 为 -300m，那么我们有 90% 的概率深度更接近地表。

🚀🚀 让我们通过这个应用程序以 3D 交互方式探索所有插值曲面（请耐心等待，唤醒应用程序 Heroku dynos 需要几秒钟）。 🚀🚀

原文链接：Exploring spatial interpolation

BimAnt翻译整理，转载请标明出处