三维均匀网格中随机散乱数据的插值与外推

Question

我有一个256x256x32的网格，由x、y和z上的规则间隔点组成，并有一个相关的变量“a”。我还有一组随机分散的点，在一个更为有限的x，y，z空间中，有一个相关的变量“b”。我基本上想做的是插值和外推我的随机数据到一个与“a”立方体匹配的规则间隔网格中，如下所示：

到目前为止，我已经使用了scipy的griddata来实现插值，它似乎工作得很好，但是它无法处理外推（据我所知），并且输出急剧地截断为“nan”值。在研究这个问题时，我遇到了一些人第二次使用griddata，用“nearest”作为插值方法来填充“nan”值。我试过了，但结果似乎不可靠。如果我在“线性”模式下使用fill_值，会得到更合适的结果，但目前它更像是一个骗局，因为fill_值必须是一个常量。在

我注意到MATLAB有一个ScatteredInterpolant类，它似乎能做我想做的事情，但是我在Python中找不到一个等价的类，也不知道如何在3D中高效地实现这样一个例程，任何帮助都将不胜感激。在

我用于插值的代码如下：

x, y, z, b = np.loadtxt(scatteredfile, unpack = True)

# Create cube to match aCube dimensions
xi = np.linspace(-xmax_aCube, xmax_aCube, 256)
yi = np.linspace(-ymax_aCube, ymax_aCube, 256)
zi = np.linspace(zmin_aCube, zmax_aCube, 32)

# Interpolate scattered points
X, Y, Z = np.meshgrid(xi, yi, zi)
bCube = griddata((x, y, z), b, (X, Y, Z), method = 'linear')    

Answer 1

这种讨论适用于任何维度。对于你的3D例子，让我们先谈谈计算几何，来理解为什么部分区域从griddata给出NaN。在

体积中的分散点构成了一个具有以下特性的几何图形：

• 表面总是凸的（顾名思义）
• 在不破坏凸性的情况下，形状的体积是尽可能小的
• 曲面（在3d中）是三角剖分和闭合的

不那么正式，凸壳（which you can compute easily with scipy）就像是在一个框架上拉伸一个气球，框架的角点是分散的簇的最外面的点。在

在气球内部的规则网格位置，你被已知的点包围着。可以对这些位置进行插值。在它之外，你必须推断。在

外推法是很难的。没有一般的规则来解释。。。这是具体问题。在这个区域，像griddata这样的算法选择返回NaN-这是通知科学家必须选择一种合理的推断方式的最安全的方法。在

让我们来看看做那件事的一些方法。在

1。[最差]搞砸

在外壳外部指定一些标量值。In the numpy docs您将看到这是通过以下方式完成的： s=平均值（b） bCube=网格数据（（x，y，z），b，（x，y，z），方法='线性'，填充值=s）

缺点：这会在赫尔边界的插值场中产生尖锐的不连续性，严重偏离平均标量场值，并且不尊重数据的函数形式。在

2。[下一个最差]“混搭”

假设在你的域的角落，你应用了一些值。这可能是与分散点相关联的标量场的平均值。在

抱歉，这是伪代码，因为我根本不使用numpy，但可能会很清楚

# With a unit cube, and selected scalar value
x, y, z, b = np.loadtxt(scatteredfile, unpack = True)
s = mean(b)
x.append([0 0 0 0 1 1 1 1])
y.append([0 0 1 1 0 0 1 1])
z.append([0 1 0 1 0 1 0 1])
b.append([s s s s s s s s])
# drop in the rest of your code

缺点：这会在赫尔边界的插值场的梯度中产生尖锐的不连续性，相当严重地偏差平均标量场值，并且不尊重数据的函数形式。在

3。[仍然相当糟糕]最近的邻居

对于每个常规的NaN点，找到最近的非NaN点并分配该值。这是有效和稳定的，但粗糙的，因为你的领域可以结束图案特征（如条纹或光束从船体辐射出去），往往视觉上不吸引人，或更糟糕的是，在数据平滑方面不可接受

根据数据的密度，可以使用最近的分散数据点，而不是最近的非NaN规则点。这可以通过（再次，伪代码）来实现：

^{pr2}$

使用MATLAB基于delaunay的方法将揭示在一行代码中实现类似功能的更强大的方法，但是numpy在这里看起来有点有限。在

4。[不总是很糟糕]自然权重

很抱歉在这一节中解释得不好，我从来没有写过算法，但我相信对自然邻域技术的一些研究会让你远远地

使用带有某个参数D的距离加权函数，该参数可能与长方体的长度相似，或者是长方体长度的两倍。你可以调整。对于每个NaN位置，计算到每个分散点的距离。在

# Don't do it this way for anything but small matrices - this is O(NM) 
# and it can be done much more effectively (e.g. MATLAB has a quick 
# natural weighting option), but for illustrative purposes:
for each NaN point 1:N
    for each scattered point 1:M
        calculate a basis function using inverse distance from NaN to point, normalised on D, and store in a [1 x M] vector of weights
Multiply weights by the b value, summate and divide by M   

你基本上想得到一个函数，它在离船体的距离D处，平滑地达到B的平均强度，但在边界处与船体重合。在远离边界的地方，它在最近的点上的权重最大。在

优点：非常稳定，并且相当连续。由于加权，在单个数据点比最近邻点更能抵抗噪声。在

5。【英雄摇滚明星】功能形式假设

你对物理学了解多少？假设函数一种形式，它代表了你期望物理所做的事情，然后对散乱的数据进行最小二乘（或类似的）拟合。使用函数稳定外推法。在

一些有助于构建函数的好主意：

• 你期望对称性还是周期性？在
• b是向量场的一个分量吗？它有一些像零散度的性质？在
• 方向性：你希望所有的角都一样吗？或者是一个方向上的线性变化？在
• 字段b是在某个时间点上吗？也许可以使用平滑的时间序列来获得一个基本函数？在
• 有没有已知的形式，像高斯或二次方？在

一些例子：

• b代表通过一个物体的激光束的强度。你期望入口侧名义上与出口相同，其他四个边界为零强度。强度将具有同心高斯分布。

• b是不可压缩流体中速度场的一个分量。流体必须是无散度的，所以在NaN区域产生的任何场也必须是无散度的，所以你应用这个条件。

• b代表房间的温度。你期望顶部温度更高，因为热空气上升。

• b代表机翼上的升力，通过三个独立变量进行测试。你可以很容易地在停机坪查看电梯，所以要确切地知道它在空间的某些部分会是什么。

优点/缺点：把这件事做好，那就太棒了。如果弄错了，特别是非线性函数形式，它将变得非常错误，并可能导致非常不稳定的结果。在

健康警告你不能假设一个函数形式，得到很好的结果，然后用它们来证明函数形式是正确的。那真是糟糕的科学。表单需要表现良好并且与数据分析无关。在

Answer 2

如果你的散点很好地符合立方体形状，一种方法可以是使用^{}插值到适合点云的规则数据网格上（因此避免了nans），然后使用这个规则的值网格作为^{}的输入，这确实有助于线性外推（但需要一个规则的网格作为输入）。在

通过这种方式，您可以像以前一样使用^{}来处理分散点的凸壳内的所有点，并且可以使用^{}来估计返回为nan的点。在

这还远远不够完美，但我认为它更接近于实现你所期待的。在

优点：

• 避免尖锐的不连续性。在
• 捕捉数据集边缘的基本线性趋势，而不必知道函数形式。在
• 尊重数据中的不对称性（例如，在远距离时，数据集的一方的值可能比另一方大，因此在远距离时，数据集的一方的值可能比另一方大。）

缺点：

• 这种方法的有效性在很大程度上取决于一个立方体的大小，你可以在你最初散布的点的凸包内容纳多大。如果你的数据是尖峰的/不规则的，那么即使是凸包边缘上的点也可能被外推到离嵌套立方体边缘相当远的地方，这会导致错误，因为外推法不会考虑到位于立方体外部的更近的数据点。在
• 线性外推法将受到数据中噪声的严重影响 在点云的边缘。在
• 做两组插值的计算成本。在