如何在Python中使用克里格法对二维空间数据进行插值？

Question

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

\n

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

\n

import numpy as np\ncoordinates = 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]])\nobservations = np.array([1.0,0.5,0.75,-1.0,0.0,1.0])\n

\n

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

\n

x = coordinates[:,0]\ny = coordinates[:,1]\n

\n

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

\n

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

\n

\n

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

\n

Answer 1

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

\n

首先，我们需要将 Numpy 数组转换coordinates为observationsOpenTURNSSample对象：

\n

import openturns as ot\ninput_train = ot.Sample(coordinates)\noutput_train = ot.Sample(observations, 1)\n

\n

该数组coordinates的形状为(6, 2)，因此它被转换为Sample大小为\xc2\xa06 且维度为\xc2\xa02 的数组。该数组observations的 shape(6,)是不明确的：它将是Samplesize\xc2\xa06 和维度\xc2\xa01 的 a，还是Samplesize\xc2\xa01 和维度\xc2\xa06 的 a？为了澄清这一点，我们在类构造函数的调用中指定维度\xc2\xa0(1) Sample。

\n

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

\n

inputDimension = 2\nbasis = ot.ConstantBasisFactory(inputDimension).build()\ncovariance_kernel = ot.SquaredExponential([1.0]*inputDimension, [1.0])\nalgo = ot.KrigingAlgorithm(input_train, output_train,\n                           covariance_kernel, basis)\n

\n

然后我们拟合趋势值和协方差核的参数（幅度参数和尺度参数）并获得元模型：

\n

# Fit\nalgo.run()\nresult = algo.getResult()\nkrigingMetamodel = result.getMetaModel()\n

\n

结果krigingMetamodel是 a Function，它接受 2DPoint作为输入并返回 1D Point。它预测兴趣的数量。为了说明这一点，让我们构建 2D 域 [0,1]\xc3\x97[0,1] 并使用规则网格将其离散化：

\n

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

\n

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

\n

# Predict\nvertices = myMeshBox.getVertices()\npredictions = krigingMetamodel(vertices)\n

\n

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

\n

# Format for plot\nX = np.array(vertices[:, 0]).reshape((ny, nx))\nY = np.array(vertices[:, 1]).reshape((ny, nx))\npredictions_array = np.array(predictions).reshape((ny,nx))\n

\n

以下脚本生成情节：

\n

# Plot\nimport matplotlib.pyplot as plt\nfig = plt.figure()\nplt.pcolor(X, Y, predictions_array)\nplt.colorbar()\nplt.show()\n

\n

\n

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

\n

该元模型是坐标的平滑函数：其平滑度随着协方差核平滑度的增加而增加，并且平方指数协方差核恰好是平滑的。

\n