求助：如何在Python中为3D数据点拟合二次曲面

刚好之前做过类似的二次曲面拟合任务，来给你拆解下怎么实现这个vertical_distance函数，顺便结合你的场景说清楚逻辑～

不管你最终确定哪种二次曲面，通用的表达式可以写成：

$F(x,y,z) = ax² + by² + cz² + dxy + exz + fyz + gx + hy + iz + j = 0$
这里的[a,b,c,d,e,f,g,h,i,j]就是我们要拟合的核心参数，后续的距离计算都基于这个形式。

点$p_0=(x0,y0,z0)$到曲面的垂直距离，本质是找曲面上的点$p=(x,y,z)$，满足两个关键条件：

• $p$严格在曲面上：$F(x,y,z)=0$
• 向量$\overrightarrow{p_0p}$与曲面在$p$处的法向量平行（因为垂直于曲面的方向就是法向量方向）

曲面在$p$处的法向量是$F$的梯度：
$\nabla F = (2ax + dy + ez + g, 2by + dx + fz + h, 2cz + ex + fy + i)$

因为$\overrightarrow{p_0p}$和法向量平行，我们可以引入参数$t$，令：
$x = x0 + t \cdot n_x$
$y = y0 + t \cdot n_y$
$z = z0 + t \cdot n_z$
其中$n=(n_x,n_y,n_z)$是$p$处的法向量。把这三个式子代入曲面方程，解出$t$后，对应的$p$就是曲面上离$p0$最近的垂直点，两点的欧氏距离就是我们要的垂直距离。

因为代入后得到的是关于$t$的高次方程，没有通用解析解，所以用数值方法求解最方便。下面是完整的实现示例：

首先导入依赖库：

import numpy as np
from scipy.optimize import fsolve

定义二次曲面的方程函数：

def surface_equation(params, x, y, z):
    """计算二次曲面方程的值，params是10个系数的列表"""
    a, b, c, d, e, f, g, h, i, j = params
    return a * x**2 + b * y**2 + c * z**2 + d * x * y + e * x * z + f * y * z + g * x + h * y + i * z + j

定义曲面的梯度（法向量）函数：

def surface_gradient(params, x, y, z):
    """计算曲面在(x,y,z)处的法向量"""
    a, b, c, d, e, f, g, h, i, _ = params
    dx = 2 * a * x + d * y + e * z + g
    dy = 2 * b * y + d * x + f * z + h
    dz = 2 * c * z + e * x + f * y + i
    return np.array([dx, dy, dz])

核心的vertical_distance函数实现：

def vertical_distance(p0, params):
    """
    计算点p0到二次曲面的垂直距离，以及对应的曲面上的点
    :param p0: 三维点，格式为(x0, y0, z0)
    :param params: 二次曲面的10个系数列表
    :return: (垂直距离, 曲面上的对应点)
    """
    x0, y0, z0 = p0

    # 定义关于t的隐式方程：代入曲面方程后等于0
    def solve_for_t(t):
        # 先计算当前t对应的曲面上的点的坐标
        x_candidate = x0 + t * surface_gradient(params, x0 + t, y0 + t, z0 + t)[0]
        y_candidate = y0 + t * surface_gradient(params, x0 + t, y0 + t, z0 + t)[1]
        z_candidate = z0 + t * surface_gradient(params, x0 + t, y0 + t, z0 + t)[2]
        # 返回曲面方程的值，目标是让它等于0
        return surface_equation(params, x_candidate, y_candidate, z_candidate)

    # 用t=0作为初始猜测（也就是p0本身），解方程得到t
    t_opt = fsolve(solve_for_t, 0.0)[0]

    # 计算最终的曲面上的点
    x_surf = x0 + t_opt * surface_gradient(params, x0 + t_opt, y0 + t_opt, z0 + t_opt)[0]
    y_surf = y0 + t_opt * surface_gradient(params, x0 + t_opt, y0 + t_opt, z0 + t_opt)[1]
    z_surf = z0 + t_opt * surface_gradient(params, x0 + t_opt, y0 + t_opt, z0 + t_opt)[2]

    # 计算欧氏距离
    dist = np.sqrt((x0 - x_surf)**2 + (y0 - y_surf)**2 + (z0 - z_surf)**2)
    return dist, (x_surf, y_surf, z_surf)

• 如果之后要做曲面拟合，你可以把所有数据点的垂直距离平方和作为损失函数，用scipy.optimize.minimize来优化参数params
• 如果你的曲面形式可以简化（比如旋转抛物面、椭球面等），可以减少参数数量，这样解方程和拟合都会更高效
• 关于误差假设：如果假设数据点有高斯误差，用最小二乘（距离平方和）是合理的；如果是鲁棒性需求，可以用绝对值和作为损失

内容的提问来源于stack exchange，提问作者Annika