python - 随机散布数据的内插和外推到 3D 均匀网格

Question

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



到目前为止，我已经使用 scipy 的 griddata 来实现插值，这似乎工作正常，但它无法处理外推(据我所知)并且输出急剧截断为“nan”值。在研究这个问题时，我遇到了几个人第二次使用 griddata 使用“最近”作为填充“nan”值的插值方法。我试过这个，但结果似乎不可靠。如果我使用带有“线性”模式的 fill_Value 会获得更合适的外观结果，但目前它更像是一种软糖，因为 fill_Value 必须是一个常量。

我注意到 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 案例，让我们先谈谈计算几何，以了解为什么部分区域给出 NaN来自 griddata .

体积中的分散点构成一个凸包；具有以下属性的几何形状:

表面总是凸的(顾名思义)

在不违反凸性的情况下，形状的体积尽可能小

表面(在 3d 中)被三角化并闭合

不太正式地，凸包 ( which you can compute easily with scipy ) 就像在框架上拉伸(stretch)气球，其中框架角是分散簇的最外点。

在气球内的常规网格位置，您被已知点包围。您可以插入到这些位置。在它之外，你必须推断。

外推很难。 没有关于如何去做的一般规则......它是特定于问题的。在那个区域，像 griddata 这样的算法选择返回 NaN - 这是告知科学家他/她必须选择合理的外推方式的最安全方式。

让我们通过一些方法来做到这一点。

1. [最糟糕] 搞砸了

在船体外分配一些标量值。 In the numpy docs你会看到这是通过以下方式完成的:
s = 平均值(b)
bCube = griddata((x, y, z), b, (X, Y, Z), method = 'linear', fill_value=s)

缺点:这会在船体边界的内插场中产生明显的不连续性，严重偏置平均标量场值并且不尊重数据的函数形式。

2. [NEXT WORST]“混合拙劣”

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

抱歉，这是伪代码，因为我根本不使用 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. [STILL PRETTY BAD] 最近的邻居

对于每个常规 NaN 点，找到最近的非 NaN 并分配该值。这是有效且稳定的，但很粗糙，因为您的场最终可能会出现图案特征(例如从船体辐射出的条纹或光束)，通常在视觉上没有吸引力，或者更糟糕的是，在数据平滑度方面是 Not Acceptable

根据数据的密度，您可以使用最近的分散数据点而不是最近的非 NaN 常规点。这可以通过(再次，伪代码)简单地完成:

bCube = griddata((x, y, z), b, (X, Y, Z), method = 'linear', fill_value=nan)
bCubeNearest = griddata((x, y, z), b, (X, Y, Z), method = 'nearest')
indicesMask = isNan(bCube)
# Use nearest interpolation outside the hull, keeping linear interpolation inside.
bCube(indicesMask) = bCubeNearest(indicesMask)

使用 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.【HEROIC ROCKSTAR】功能形式假设

你对物理学了解多少？假设一个函数形式代表您期望物理所做的事情，然后对该形式进行最小二乘法(或某种等效的)拟合散布数据。使用该函数来稳定外推。

一些可以帮助您构建函数的好主意:

你期望对称性还是周期性？

b 是具有零散度等属性的向量场的一个分量吗？

方向性:您是否希望所有角落都相同？或者可能是一个方向的线性变化？

场 b 是否在某个时间点 - 也许可以使用平滑的测量时间序列来提出基本函数？

是否已经存在像高斯或二次这样的已知形式？

一些例子:

b 表示通过体积的激光束的强度。您希望入口侧在名义上与导出相同，其他四个边界为零强度。强度将具有同心高斯分布。

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

b 代表房间内的温度。您预计顶部的温度会更高，因为热空气会上升。

b 代表机翼上的升力，在三个独立变量上进行测试。您可以很容易地在摊位上查找升降机，因此您可以确切地知道它在空间的某些部分会是什么。

优点/缺点:做对了，它会很棒。弄错了，尤其是非线性函数形式，它会出错并且会导致非常不稳定的结果。

健康警示 你不能假设一个函数形式，得到漂亮的结果，然后用它们来证明函数形式是正确的。那只是糟糕的科学。该表单需要表现良好，并且独立于您的数据分析。