看起来你已经做了很扎实的尝试——用3D插值把非均匀数据转成均匀网格再积分的思路是可行的，但确实会在函数平缓区域浪费大量计算资源，还可能因为网格密度限制影响精度。你想要的不用插值、直接基于非均匀采样数据计算积分的方法是完全可行的，核心思路是换个角度处理：不对整个3D空间做插值，而是针对每个(x,y)位置，直接利用该位置附近的非均匀z采样点完成积分。

核心逻辑 #

我们的目标是得到每个(x,y)处的积分结果 $F(x,y) = \int_a^b f(x,y,z)dz$。由于原始数据是大量的$(x,y,z, f(x,y,z))$样本，我们可以拆分任务：

1. 先定义最终需要的$(x,y)$输出网格（这个网格完全可以按需调整，比如在函数变化剧烈的区域加密，平缓区域稀疏，不必强制均匀）；
2. 对每个目标$(x_0,y_0)$网格点，找到所有邻近的采样点，提取对应的z值和函数值；
3. 对这些非均匀的z采样点直接执行数值积分，得到该$(x_0,y_0)$处的积分结果。

这种方法既保留了原始非均匀采样的优势（在函数变化快的地方有更多样本），又避免了全3D插值带来的资源浪费。

具体实现代码 #

我们基于你提供的示例场景来实现这个方法：

1. 生成示例数据（和你的代码一致）

import numpy as np

def f(x,y,z) :
    "function to be integrated"
    return x**2 + y**2 + z**2

# 生成非均匀采样点
n_points = 10000
x = np.random.random_sample(n_points)
y = np.random.random_sample(n_points)
z = np.random.random_sample(n_points)
vals = f(x,y,z)

2. 基于非均匀采样直接计算积分

from scipy.spatial import KDTree
from scipy.integrate import trapezoid
import matplotlib.pyplot as plt

# 定义输出的(x,y)网格（可按需调整密度）
n_xy_grid = 25
x_grid = np.linspace(0, 1, n_xy_grid)
y_grid = np.linspace(0, 1, n_xy_grid)
X_grid, Y_grid = np.meshgrid(x_grid, y_grid)

# 用KDTree快速查找邻近点：将(x,y)作为特征构建索引
xy_points = np.stack([x, y], axis=1)
kdtree = KDTree(xy_points)

# 每个(x,y)点取邻近的50个采样点（可根据样本量调整）
k_neighbors = 50

# 初始化结果数组
integral_result = np.zeros_like(X_grid)

# 遍历每个(x,y)网格点计算积分
for i, x0 in enumerate(x_grid):
    for j, y0 in enumerate(y_grid):
        # 查找距离最近的k个采样点
        _, neighbor_indices = kdtree.query([x0, y0], k=k_neighbors)
        # 提取对应的z值和函数值
        z_group = z[neighbor_indices]
        vals_group = vals[neighbor_indices]
        
        # 积分要求z值有序，先排序
        sort_idx = np.argsort(z_group)
        z_sorted = z_group[sort_idx]
        vals_sorted = vals_group[sort_idx]
        
        # 用梯形法计算非均匀z轴的积分（兼容性好，精度足够）
        integral = trapezoid(vals_sorted, x=z_sorted)
        integral_result[j, i] = integral

3. 和解析解对比验证

# 计算解析解用于对比
def integral_f_dz(x,y,a,b) :
    "analytical expression of definite integral of f through z axis"
    return 1/3 * ( -a**3 + 3*(b-a)*(x**2+y**2)+b**3 )

analytical_result = integral_f_dz(X_grid, Y_grid, 0, 1)

# 绘图对比
fig, (ax1, ax2) = plt.subplots(ncols=2, constrained_layout=True, figsize=(10, 4))
for ax in [ax1, ax2]:
    ax.set_aspect('equal')

ax1.set_title("解析解结果")
cont1 = ax1.contourf(X_grid, Y_grid, analytical_result, cmap='viridis')

ax2.set_title("非均匀采样直接积分结果")
cont2 = ax2.contourf(X_grid, Y_grid, integral_result, cmap='viridis')

plt.colorbar(cont1, ax=ax1)
plt.colorbar(cont2, ax=ax2)
plt.show()

方法优化与拓展 #

• 分箱替代K近邻：如果你的采样点数量极大，可以将(x,y)空间划分为固定bin，每个bin内的所有z点统一积分，速度会更快；
• 加权积分：为了提升精度，可以根据采样点到$(x_0,y_0)$的距离赋予权重（距离越近权重越高），再进行积分；
• 采样点筛选：如果你的原始数据中存在离群点，可以在提取z_group后先剔除异常值再积分；
• 辛普森法适配：若每个(x,y)附近的z采样点数量足够多且分布相对均匀，也可以使用scipy.integrate.simpson，但注意辛普森法对样本数量的奇偶性有要求。

这种方法的计算效率和精度都远优于全3D插值后积分的方案，完全适配你拥有大量非均匀采样点的场景。

备注：内容来源于stack exchange，提问作者user28616287