如何在非正交有序网格上计算数据梯度？NumPy实现咨询

问题背景 #

用np.meshgrid生成梯形xy坐标网格：x方向间距恒定，y坐标随x缩放，左侧边长为2、右侧边长为1。需要计算函数f(x,y)在y方向和每行点的切线方向（xy空间内的行方向）的梯度，但对np.gradient的参数用法存疑，官方文档理解难度大。

示例代码：

import numpy as np

x = np.linspace(0, 1, num=20)
y = np.linspace(-1/2, 1/2, num=10)
x_grid, y_grid = np.meshgrid(x, y)
y_grid = y_grid * (2 - x_grid)  # 生成梯形网格

f = x_grid ** 2 - y_grid ** 3  # 示例函数

# 期望实现类似逻辑，但不确定正确用法
# gradient = np.gradient(f, x_grid, y_grid)
# 希望gradient[0]为行方向梯度，gradient[1]为y方向梯度

现有循环实现 #

发现np.gradient仅支持传入一维坐标数组，目前通过循环逐行逐列计算，但效率较低：

# 计算y方向梯度（逐列处理）
df_dy_loop = np.array([np.gradient(f[:, j], y_grid[:, j], edge_order=2) for j in range(y_grid.shape[1])]).T

# 计算行方向梯度：先求每行点的弧长坐标，再逐行计算梯度
row_distances = np.linalg.norm([x_grid - x_grid[:, [0]], y_grid - y_grid[:, [0]]], axis=0)
df_drow_loop = np.array([np.gradient(f[j], row_distances[j], edge_order=2) for j in range(x_grid.shape[0])])

向量化优化方案 #

1. y方向梯度的向量化计算 #

y方向对应数组的axis=0（列内方向），直接将y_grid作为坐标参数传入np.gradient并指定轴即可，无需循环：

df_dy = np.gradient(f, y_grid, axis=0, edge_order=2)

np.gradient会自动对每个列（axis=0方向）使用对应的y_grid列作为坐标，完成向量化计算。

2. 行方向梯度的向量化计算 #

行方向是xy空间的切线方向，需先计算每行的累积弧长，再沿axis=1（行方向）计算梯度：

步骤1：计算每行的累积弧长

# 计算每行相邻点的x、y差值
dx = np.diff(x_grid, axis=1)
dy = np.diff(y_grid, axis=1)
# 计算相邻点的欧氏距离
seg_dist = np.sqrt(dx**2 + dy**2)
# 累积求和得到弧长，开头补0（起点弧长为0）
row_arclen = np.concatenate([np.zeros((x_grid.shape[0], 1)), np.cumsum(seg_dist, axis=1)], axis=1)

步骤2：向量化计算行方向梯度

df_drow = np.gradient(f, row_arclen, axis=1, edge_order=2)

结果验证 #

可以对比向量化版本和循环版本的结果误差，确认正确性：

print(np.max(np.abs(df_dy - df_dy_loop)))  # 误差接近0
print(np.max(np.abs(df_drow - df_drow_loop)))  # 误差接近0

内容的提问来源于stack exchange，提问作者Alex V.