在气象、海洋、遥感等领域的数值数据处理中，一个高频出现的需求是：把规则网格数据（例如经纬度网格上的温度场）插值到一组任意分布的观测点上。如果网格数据完全干净，scipy.interpolate.griddata 配合 method='cubic' 就能直接完成任务。但现实往往是网格里散落着若干 NaN，比如陆地掩膜、传感器缺测、云遮挡等场景。NaN 一旦混进输入数组，cubic 方法要么直接抛出 QhullError，要么悄悄把结果污染成大片无效值，排查起来非常头疼。本文将从原因分析入手，给出几种完整可运行的解决方案。
一、为什么 NaN 会让 griddata 的 cubic 方法失效
先看一段典型的失败代码。假设我们有一个 20x20 的规则网格，网格点坐标由 meshgrid 生成，数据中人为放置了若干 NaN：
import numpy as np
from scipy.interpolate import griddata
# 构造 20x20 规则网格
x = np.linspace(0, 10, 20)
y = np.linspace(0, 10, 20)
xx, yy = np.meshgrid(x, y)
values = np.sin(xx) * np.cos(yy)
# 模拟缺测，随机撒入 NaN
rng = np.random.default_rng(42)
mask = rng.random(values.shape) < 0.1
values[mask] = np.nan
# 目标插值点：非结构化散点
pts_x = rng.uniform(0, 10, 50)
pts_y = rng.uniform(0, 10, 50)
target = np.column_stack([pts_x, pts_y])
# 直接调用 cubic，问题就出在这里
result = griddata(
np.column_stack([xx.ravel(), yy.ravel()]),
values.ravel(),
target,
method='cubic'
)运行后大概率会遇到两种情况之一：一是直接抛出 scipy.spatial.QhullError，二是返回的结果中大量元素为 nan。根本原因在于 method='cubic' 底层依赖 Delaunay 三角剖分和 CloughTocher 分片三次样条，NaN 参与 Qhull 的几何计算时，任何涉及 NaN 的比较结果都是 False，导致剖分失败或三角形面片数据非法。即使剖分勉强完成，任何一个顶点为 NaN 的三角形内部的插值结果也会被传染为 NaN。
另一个容易忽视的坑是 NaN 的判断方式。一定要用 np.isnan() 而不是 np.isfinite() 的反向逻辑随意混用，前者只筛 NaN，后者还会把正负无穷一起处理。如果数据来自二进制文件或网络接口，正负无穷同样会破坏 cubic 插值，建议统一用 np.isfinite() 做过滤，一步到位。
二、方案一：先剔除无效点再做 cubic 插值
最直接的思路是，在调用 griddata 之前把 NaN 对应的网格点从输入中删掉。规则网格删掉若干点后，剩余点就变成了事实上的非结构化散点，而这恰好是 griddata 擅长处理的输入形式。代码如下：
import numpy as np from scipy.interpolate import griddata # 展平后过滤掉无效点 points = np.column_stack([xx.ravel(), yy.ravel()]) vals_flat = values.ravel() valid = np.isfinite(vals_flat) points_valid = points[valid] vals_valid = vals_flat[valid] # 用干净的散点做 cubic 插值 result = griddata(points_valid, vals_valid, target, method='cubic') # 检查结果，位于凸包之外的点会得到 nan print(np.isnan(result).sum(), '个目标点落在凸包外或插值失败')
这个方案的关键点在于过滤逻辑 valid = np.isfinite(vals_flat)，它同时处理了 NaN 和 inf 两类脏数据。过滤之后剩余点集仍然要满足一定密度才能保证 cubic 的精度——如果 NaN 恰好连成一片大区域（比如整块陆地区域被掩膜），剔除后该区域的三角形会非常狭长，CloughTocher 样条在这些位置的导数估计会明显失真，插值结果可能出现震荡。
还需要注意凸包问题。剔除点之后，有效点集的凸包可能比原始网格小，靠近边界的插值点会得到 NaN。如果业务上必须给这些点一个值，可以改用 method='linear' 兜底，或者在后处理阶段用最近邻方法补值：
# 用最近邻结果兜底填补凸包外的点 nn_result = griddata(points_valid, vals_valid, target, method='nearest') final = np.where(np.isnan(result), nn_result, result)
这种两级策略在工程上非常常用，核心区域享受 cubic 的高精度，边缘区域用最近邻保底，兼顾了精度和覆盖率。
三、方案二：用 RectBivariateSpline 处理规则网格
既然数据本身来自规则网格，还有一条路可走：先把 NaN 在网格层面填补掉，然后用专门面向规则网格的 RectBivariateSpline 做三次样条插值。填补 NaN 可以用简单的邻域平均，也可以用 scipy.interpolate.NearestNDInterpolator 只对缺失位置插值：
import numpy as np
from scipy.interpolate import NearestNDInterpolator, RectBivariateSpline
# 第一步：只对缺失位置做最近邻填补
finite = np.isfinite(values)
interp_nearest = NearestNDInterpolator(
np.column_stack([xx[finite], yy[finite]]), values[finite]
)
filled = values.copy()
filled[~finite] = interp_nearest(xx[~finite], yy[~finite])
# 第二步：规则网格上的三次样条
spline = RectBivariateSpline(y, x, filled, kx=3, ky=3)
result2 = spline.ev(pts_y, pts_x) # ev 接受任意散点坐标这里有两个细节值得展开。第一，RectBivariateSpline 要求输入网格坐标必须严格单调递增，所以传入的是一维的 x、y 向量而不是二维的 xx、yy；第二，查询散点时要用 ev 方法，它接受两组坐标数组并按元素配对，正好匹配非结构化插值点的需求。
这种方案的优势是性能极佳。规则网格样条的系数只需要计算一次，之后每次查询都是小开销运算，当目标点数量达到十万级甚至百万级时，比 griddata 每次都重新做三角剖分快一个数量级以上。缺点是填补步骤会引入人为数据，如果缺失区域面积大，填补值的误差会通过三次样条的光滑性扩散到周边区域，这一点在精度敏感的场景必须评估。
四、方案三：RBFInterpolator 径向基插值替代
当 NaN 分布很不规则、既不想剔除点也不想填补网格时，径向基函数插值是另一条可行路线。scipy.interpolate.RBFInterpolator 原生支持任意散点，而且对输入点是否规则没有任何要求：
import numpy as np
from scipy.interpolate import RBFInterpolator
rbf = RBFInterpolator(
points_valid, vals_valid,
neighbors=32, # 只用邻近点，控制计算量
smoothing=1e-3 # 少量平滑，抑制噪声放大
)
result3 = rbf(target)neighbors 参数是性能的关键。不设置它时，插值器会构建全量核矩阵，内存和时间的复杂度都是点数的平方；设置成 32 或 64 之后，每个目标点只用最近的邻居参与计算，速度大幅提升，精度损失通常可以忽略。smoothing 参数则允许插值不严格穿过数据点，对含观测噪声的数据尤其有价值，这是 cubic 样条不具备的能力。
三种方案如何选？给一个简单的判断标准：NaN 稀疏且需要严格过点时用方案一的 cubic 剔除法；数据源是规则网格且查询点量大时用方案二的样条路线；缺失模式复杂、数据带噪声时用方案三的 RBF。实际项目中也可以先用小样本对三种方案做交叉验证，量出各自在有效数据点上的误差，再决定正式方案，这比凭直觉选择可靠得多。