如何在 Python 中使用克里金法对 2D 空间数据进行插值？

Posted 2023-03-12

【中文标题】如何在 Python 中使用克里金法对 2D 空间数据进行插值？

【英文标题】：How to interpolate 2D spatial data with kriging in Python?

【发布时间】：2020-04-07 08:59:03

【问题描述】：

我有一个空间二维域，例如 [0,1]×[0,1]。在该域中，有 6 个点观察到了一些感兴趣的标量（例如，温度、机械应力、流体密度等）。如何预测未观察点的兴趣数量？换句话说，我如何在 Python 中插入空间数据？

例如，考虑 2D 域中点的以下坐标（输入）和感兴趣数量的相应观察值（输出）。

import numpy as np
coordinates = np.array([[0.0,0.0],[0.5,0.0],[1.0,0.0],[0.0,1.0],[0.5,1.],[1.0,1.0]])
observations = np.array([1.0,0.5,0.75,-1.0,0.0,1.0])

X 和 Y 坐标可以通过以下方式提取：

x = coordinates[:,0]
y = coordinates[:,1]

以下脚本创建一个散点图，其中黄色（分别为蓝色）代表高（分别为低）输出值。

import matplotlib.pyplot as plt
fig = plt.figure()
plt.scatter(x, y, c=observations, cmap='viridis')
plt.colorbar()
plt.show()

我想使用克里金法预测二维输入域内规则网格上的感兴趣标量。如何在 Python 中做到这一点？

【参考方案1】：

在OpenTURNS 中，KrigingAlgorithm 类可以根据特定输入点的已知输出值估计高斯过程模型的超参数。然后KrigingAlgorithm 的getMetamodel 方法返回一个插入数据的函数。

首先，我们需要将 Numpy 数组 coordinates 和 observations 转换为 OpenTURNS Sample 对象：

import openturns as ot
input_train = ot.Sample(coordinates)
output_train = ot.Sample(observations, 1)

数组coordinates 的形状为(6, 2)，因此变成了大小为6、维度为2 的Sample。数组observations 的形状为(6,)，这是不明确的：大小为 6 且维度为 1 的 Sample，还是大小为 1 且维度为 6 的 Sample？为了澄清这一点，我们在对 Sample 类构造函数的调用中指定维度 (1)。

下面，我们定义一个具有恒定趋势函数和平方指数协方差核的高斯过程模型：

inputDimension = 2
basis = ot.ConstantBasisFactory(inputDimension).build()
covariance_kernel = ot.SquaredExponential([1.0]*inputDimension, [1.0])
algo = ot.KrigingAlgorithm(input_train, output_train,
                           covariance_kernel, basis)

然后我们将趋势的值和协方差核的参数（幅度参数和尺度参数）进行拟合，得到一个元模型：

# Fit
algo.run()
result = algo.getResult()
krigingMetamodel = result.getMetaModel()

生成的krigingMetamodel 是一个Function，它接受一个二维Point 作为输入并返回一个一维Point。它预测感兴趣的数量。为了说明这一点，让我们构建二维域 [0,1]×[0,1] 并用规则网格对其进行离散化：

# Create the 2D domain
myInterval = ot.Interval([0.0, 0.0], [1.0, 1.0])
# Define the number of interval in each direction of the box
nx = 20
ny = 10
myIndices = [nx - 1, ny - 1]
myMesher = ot.IntervalMesher(myIndices)
myMeshBox = myMesher.build(myInterval)

使用我们的krigingMetamodel 来预测此网格上感兴趣的数量所取的值可以通过以下语句完成。我们首先将网格的vertices 获取为Sample，然后通过对元模型的一次调用来评估预测（这里不需要for 循环）：

# Predict
vertices = myMeshBox.getVertices()
predictions = krigingMetamodel(vertices)

为了使用 Matplotlib 查看结果，我们首先必须创建 pcolor 函数所需的数据：

# Format for plot
X = np.array(vertices[:, 0]).reshape((ny, nx))
Y = np.array(vertices[:, 1]).reshape((ny, nx))
predictions_array = np.array(predictions).reshape((ny,nx))

下面的脚本产生了情节：

# Plot
import matplotlib.pyplot as plt
fig = plt.figure()
plt.pcolor(X, Y, predictions_array)
plt.colorbar()
plt.show()

我们看到元模型的预测等于在观察到的输入点的观察。

这个元模型是坐标的平滑函数：它的平滑度随着协方差核平滑度的增加而增加，平方指数协方差核恰好是平滑的。

【讨论】：

看起来output_train = ot.Sample(observations,1) 是一个语法错误。对我来说，唯一可行的解​​决方案是将观察结果设为假二维数组，如下所示：[[0]、[1]、[5]、[4]]