Python 中数据的 3D 插值限制拟合函数仅在 z 上增加

Question

我有以下数据：

x = array([ 0, 0.08885313,  0.05077321,  0.05077321,  0.03807991,  0.03807991,
        0.03807991,  0.02538661,  0.02538661,  0.0126933 ,  0.0126933 ,
       -0.        , -0.        , -0.0126933 , -0.0126933 , -0.02538661,
       -0.02538661, -0.02538661, -0.03807991, -0.05077321, -0.05077321,
       -0.05077321, -0.06346652, -0.07615982, -0.11423973, -0.12693304,
       -0.13962634])
y = array([ 0, -0.15231964, -0.08885313, -0.17770625, -0.08885313, -0.10154643,
       -0.12693304, -0.07615982, -0.08885313, -0.07615982, -0.08885313,
       -0.07615982, -0.10154643, -0.08885313, -0.10154643, -0.07615982,
       -0.08885313, -0.17770625, -0.10154643, -0.08885313, -0.11423973,
       -0.12693304, -0.10154643, -0.08885313, -0.17770625, -0.12693304,
       -0.11423973])
z = array([ 0, 1.21839241, 0.78673339, 1.21839241, 0.70318648, 0.82850684,
       0.96078945, 0.64748854, 0.71014872, 0.63356406, 0.71711096,
       0.65445078, 0.77977115, 0.73103545, 0.77977115, 0.68926199,
       0.73799769, 1.19750569, 0.85635581, 0.84243133, 0.96078945,
       1.02344963, 0.93294048, 0.97471393, 1.24624138, 1.20446793,
       1.22535466])

我使用了以下代码，将 (x, y, z) 空间中的点形式的曲面拟合到我的数据。

# Find the maximum and minimum x values
i_max, = np.where(np.isclose(x, np.max(x)))
i_min, = np.where(np.isclose(x, np.min(x)))
x_max = x[i_max][0]
x_min = x[i_min][0]

# Find the corresponding y values (i.e. y values of those points where x is maximum)
y_max = y[i_max][0]
y_min = y[i_min][0]

# Function finding 2D-line coefficients, based on this answer https://stackoverflow.com/a/21566184/3715182
def line_coeffs(points):
    x_coords, y_coords = zip(*points)
    A = vstack([x_coords, ones(len(x_coords))]).T
    # y = a*x + b
    a, b = lstsq(A, y_coords, rcond=None)[0]
    return (a, b)

# Find coefficients of the two lines "limiting" all points left and right across the x-axis in the XY-plane
k1_max, k2_max = line_coeffs([(x_max, y_max), (0, 0)])
k1_min, k2_min = line_coeffs([(x_min, y_min), (0, 0)])

# Define mathematical function for curve fitting 
def func(xy, a, b, c, d, e, f):
    x, y = xy
    return a + b*x + c*y + d*x**2 + e*y**2 + f*x*y

# Create a grid of points over the surface
x_fill = np.linspace(-0.2, 0.2, 40)
y_fill = np.linspace(-0.3, 0, 40)
X_fill, Y_fill = np.meshgrid(x_fill, y_fill)

# Evaluate the function at each point on the surface
Z_fill = func((X_fill, Y_fill), *popt)

# Plot the data points and the filled surface
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x, y, z, color='blue', label='Data Points')

ax.scatter(np.where(X_fill >= (Y_fill - k2_min)/k1_min, np.where(X_fill <= (Y_fill - k2_max)/k1_max, X_fill, np.nan), np.nan), Y_fill, Z_fill, color='red', alpha=0.5, s= 2)

ax.set_xlim([-0.2, 0.2])
ax.set_ylim([-0.3, 0.2])
ax.set_zlim([0, 1.2])

ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')

plt.show()

这给出了以下以两个不同角度显示的图：

您可以看到红点与我的蓝色数据点非常吻合。然而，我预计数据将继续上升，而不是在 y=-0.3 附近回落。如何限制（或更改）我的拟合函数，以便该函数在 z 上不断增加。我希望有点像下面的图的一部分（仅在两个角度之间）：

Answer 1

您发布的第三张图看起来像围绕 (0, 0) 旋转对称，所以我建议拟合一个具有旋转对称的函数。一个简单的方法是将笛卡尔坐标转换为极坐标，并忽略 theta。

在采用这种方法之前，检查数据是否通过极坐标形式的垂直线测试是有用的，以确保像这样丢弃 theta 能够起作用。

这看起来或多或少还不错。

该图看起来并不完全笔直，所以我们将其转换为极坐标形式，然后拟合二次方程。

定义函数：

def func(xy, a, b, c):
    x, y = xy
    r = np.sqrt(x ** 2 + y ** 2)
    theta = np.arctan2(x, y)
    return a*(r - b)**2 + c

适合它：

popt, pcov = scipy.optimize.curve_fit(func, [x, y], z)

结果非常合适，但未满足您的要求之一。

该函数对于 x 或 y 不是单调的。换句话说，当您从中心移出时，它会上升，然后又下降。不过，这是可以修复的。我们知道抛物线的峰值在哪里。在半径

b

处，抛物线达到

c

的峰值，然后下降。我们可以把这个函数改成只有一半抛物线，剩下的用直线填充。

def func(xy, a, b, c):
    x, y = xy
    r = np.sqrt(x ** 2 + y ** 2)
    theta = np.arctan2(x, y)
    parabola = a*(r - b)**2 + c
    return np.where(r <= b, parabola, c)

这产生以下推断：

这为您提供了相同的拟合质量，并且参数更少，控制更多。

您还可以尝试通过包含 theta 项来降低旋转对称性的版本：

def func(xy, a, b, c, d):
    x, y = xy
    r = np.sqrt(x ** 2 + y ** 2)
    theta = np.arctan2(x, y)
    parabola = a*(r - b)**2 + c + d*np.sin(theta)
    return np.where(r <= b, parabola, c)

我没有发现这增加了多少准确性。