您可以将Delaunay三角剖分与Z值集一起提供给scipy.interpolate.LinearNDInterpolator，它应该可以为您完成这项工作。
如果你真的想自己做插值，你可以从find_simplex和transform构建它。

我们可以使用scipy.special.Delaunay()从给定的样本点返回的对象计算重心坐标，并按如下方式插值函数（可以在不使用scipy.interpolate.LinearNDInterpolator的情况下完成）：

from scipy.spatial import Delaunay
import matplotlib.pylab as plt

nx, ny = 160, 300

xv, yv = np.meshgrid(np.arange(nx), np.arange(ny))
xv = xv.flatten()
yv = yv.flatten()
pp = np.vstack((xv, yv)).T
    
n = 100

ix = np.random.choice(nx, n).tolist() + [0, 0, nx-1, nx-1]
iy = np.random.choice(ny, n).tolist() + [0, ny-1, 0, ny-1]
f = np.zeros((nx, ny))
for x, y in zip(ix, iy):
    f[x, y] = (x/160)**2 + (y/300)**2 + np.random.rand()
points = np.array(list(zip(ix, iy)))
tri = Delaunay(points)
ss = tri.find_simplex(pp)
ndim = tri.transform.shape[2]
print(n,len(np.unique(ss)))

out = np.zeros((nx, ny))
for i in np.unique(ss): # for all simplices (triangles)
    p = pp[ss == i] # all points in the simplex
    # compute the barycentric coordinates of the points
    b = tri.transform[i, :ndim].dot(np.transpose(p) - tri.transform[i, ndim].reshape(-1,1))
    αβγ = np.c_[np.transpose(b), 1 - b.sum(axis=0)] 
    sp = points[tri.simplices[i]]
    if len(αβγ) > 0:
        out[p[:,0], p[:,1]] = αβγ@f[sp[:,0], sp[:,1]]     
out = out / out.max()

plt.figure(figsize=(10,12))
plt.imshow(out, aspect='auto', cmap='Oranges'), plt.axis(), plt.title('interpolated output with\n barycentric coordinates', size=20)
plt.triplot(points[:,1], points[:,0], tri.simplices)
for p in points:
    plt.scatter(p[1], p[0], c=f[p[0],p[1]], s=50)
plt.title(f'Interpolation with Delaunay Triangulation', size=20)
plt.show()