这个讨论适用于任何维度。对于您的 3D 案例,让我们先谈谈计算几何,以了解为什么部分区域NaN
从griddata
.
体积中的散点构成一个凸包;具有以下属性的几何形状:
- 表面总是凸的(顾名思义)
- 在不违反凸性的情况下,形状的体积尽可能小
- 表面(在 3d 中)被三角剖分并闭合
不太正式,凸包(您可以使用 scipy 轻松计算)就像在框架上拉伸气球,其中框架角是分散集群的最外点。
在气球内的常规网格位置,您被已知点包围。您可以对这些位置进行插值。在它之外,你必须推断。
外推很难。没有关于如何做到这一点的一般规则......它是针对特定问题的。在那个区域,算法griddata
选择返回NaN
——这是告知科学家他/她必须选择一种合理的推断方式的最安全方式。
让我们通过一些方法来做到这一点。
1. [最糟糕] 搞砸了
在船体外部分配一些标量值。在 numpy 文档中,您将看到这是通过以下方式完成的: s = mean(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 并分配该值。这是有效且稳定的,但很粗糙,因为您的字段最终可能会出现图案化特征(如从船体向外辐射的条纹或光束),通常在视觉上不吸引人,或者更糟糕的是,在数据平滑度方面不可接受
根据数据的密度,您可以使用最近的分散数据点而不是最近的非 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] 函数形式假设
你对物理学了解多少?假设一个函数形式代表您期望物理学做的事情,然后对该形式进行最小二乘(或某种等效)拟合到分散的数据。使用函数来稳定外推。
一些可以帮助您构建函数的好主意:
- 你期待对称性还是周期性?
- 向量场的 ba 分量是否具有诸如零散度之类的特性?
- 方向性:您是否希望所有角落都相同?或者可能是一个方向的线性变化?
- 字段 b 是否在某个时间点 - 也许可以使用平滑的测量时间序列来提出基本函数?
- 是否已经有已知的形式,例如高斯或二次?
一些例子:
b 表示通过体积的激光束的强度。您希望入口侧名义上与出口相同,其他四个边界为零强度。强度将具有同心高斯分布。
b 是不可压缩流体中速度场的一个分量。流体必须是无发散的,因此在 NaN 区域产生的任何场也必须是无发散的,因此您可以应用此条件。
b 代表房间内的温度。您预计顶部的温度会更高,因为热空气会上升。
b 代表机翼上的升力,在三个独立变量上进行了测试。您可以轻松地在失速时查看电梯,因此确切地知道它在空间的某些部分会是什么。
优点/缺点:做对了,它会很棒。弄错了,尤其是非线性函数形式,它会出错,并可能导致非常不稳定的结果。
健康警告你不能假设一个函数形式,得到漂亮的结果,然后用它们来证明函数形式是正确的。那只是糟糕的科学。表单必须是行为良好且已知独立于您的数据分析的东西。