基于径向分量限制3D插值范围并约束插值单调性

现有如下3D网格点数据：

import numpy as np
x = np.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 = np.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 = np.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])

已通过Python的curve_fit实现二次多项式曲面拟合，代码如下：

from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

data = np.array([x, y, z]).T 

# 定义拟合函数
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 

# 执行拟合
popt, pcov = curve_fit(func, (x, y), z) 

# 打印优化参数
print(popt) 

# 绘制数据点与拟合曲面
fig = plt.figure() 
ax = fig.add_subplot(111, projection='3d') 
ax.scatter(x, y, z, color='blue') 
x_range = np.linspace(-0.2, 0.2, 50) 
y_range = np.linspace(-0.2, 0.2, 50) 
X, Y = np.meshgrid(x_range, y_range) 
Z = func((X, Y), *popt) 
ax.plot_surface(X, Y, Z, color='red', alpha=0.5) 

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

plt.show()

当前存在两个问题：

• 生成的曲面超出目标径向范围，需添加径向+角度限制，使插值仅在数据点与原点之间的区域内生效
• 调整后发现y<-0.3区域插值下降，需约束该区域内插值结果随y减小（更负）时仅递增

1. 径向与角度范围约束 #

通过极坐标转换，确定每个角度下数据点的最大径向距离，将超出该范围的网格点设为NaN（Matplotlib会自动忽略这些点，不绘制曲面）：

# 计算数据点的极坐标
r_data = np.sqrt(x**2 + y**2)
theta_data = np.arctan2(y, x)

# 建立角度到最大径向距离的映射
theta_unique = np.unique(np.round(theta_data, 4))
max_r_per_theta = {}
for theta in theta_unique:
    mask = np.isclose(np.round(theta_data,4), theta)
    max_r_per_theta[theta] = np.max(r_data[mask])

# 生成网格并转换为极坐标
X, Y = np.meshgrid(np.linspace(-0.2, 0.2, 50), np.linspace(-0.35, 0.2, 50))
r_grid = np.sqrt(X**2 + Y**2)
theta_grid = np.arctan2(Y, X)

# 初始化Z矩阵，超出范围设为NaN
Z = func((X, Y), *popt)
for i in range(Z.shape[0]):
    for j in range(Z.shape[1]):
        theta = np.round(theta_grid[i,j],4)
        # 找到最接近的角度对应的最大r
        closest_theta = min(max_r_per_theta.keys(), key=lambda t: abs(t - theta))
        if r_grid[i,j] > max_r_per_theta[closest_theta]:
            Z[i,j] = np.nan

2. y<-0.3区域递增约束 #

由于curve_fit不支持直接添加导数约束，改用scipy.optimize.minimize进行带约束的拟合，强制y<-0.3时曲面的偏导数∂z/∂y ≥ 0（保证z随y减小而递增）：

from scipy.optimize import minimize

# 定义损失函数：平方误差
def loss(params):
    return np.sum((func((x,y),*params)-z)**2)

# 定义约束：y<-0.3时，∂z/∂y = c + 2e*y + f*x ≥ 0
def constraint(params):
    a,b,c,d,e,f = params
    mask = y < -0.3
    if np.any(mask):
        return c + 2*e*y[mask] + f*x[mask]  # 约束所有值≥0
    return np.array([1])  # 无满足条件的点时返回正数

# 设置约束条件
cons = [{'type': 'ineq', 'fun': constraint}]

# 初始值用curve_fit的结果
initial_guess = popt

# 执行带约束的优化
result = minimize(loss, initial_guess, constraints=cons)
popt_constrained = result.x

完整修改后代码 #

import numpy as np
from scipy.optimize import curve_fit, minimize
import matplotlib.pyplot as plt

# 原始数据
x = np.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 = np.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 = np.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])

# 定义拟合函数
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 

# 先执行无约束拟合得到初始参数
popt, _ = curve_fit(func, (x, y), z)

# 计算数据点的极坐标映射（径向约束）
r_data = np.sqrt(x**2 + y**2)
theta_data = np.arctan2(y, x)
theta_unique = np.unique(np.round(theta_data, 4))
max_r_per_theta = {}
for theta in theta_unique:
    mask = np.isclose(np.round(theta_data,4), theta)
    max_r_per_theta[theta] = np.max(r_data[mask])

# 带约束的拟合（y<-0.3区域递增）
def loss(params):
    return np.sum((func((x,y),*params)-z)**2)

def constraint(params):
    a,b,c,d,e,f = params
    mask = y < -0.3
    if np.any(mask):
        # 偏导数∂z/∂y ≥0
        return c + 2*e*y[mask] + f*x[mask]
    return np.array([1])

cons = [{'type':'ineq','fun':constraint}]
result = minimize(loss, popt, constraints=cons)
popt_constrained = result.x

# 生成网格并应用双重约束
X, Y = np.meshgrid(np.linspace(-0.2,0.2,50), np.linspace(-0.35,0.2,50))
r_grid = np.sqrt(X**2 + Y**2)
theta_grid = np.arctan2(Y,X)
Z = func((X,Y),*popt_constrained)

# 应用径向约束，超出范围设为NaN
for i in range(Z.shape[0]):
    for j in range(Z.shape[1]):
        theta = np.round(theta_grid[i,j],4)
        # 找到最接近的角度
        closest_theta = min(max_r_per_theta.keys(), key=lambda t: abs(t-theta))
        if r_grid[i,j] > max_r_per_theta[closest_theta]:
            Z[i,j] = np.nan

# 绘图
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x,y,z, color='blue', label='原始数据点')
ax.plot_surface(X,Y,Z, color='red', alpha=0.5, label='约束拟合曲面')
ax.set_xlim([-0.2,0.2])
ax.set_ylim([-0.35,0.2])
ax.set_zlim([0,1.4])
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.legend()
plt.show()

关键改动说明 #

• 径向约束：通过极坐标转换，为每个角度区间设定最大允许径向距离，超出部分设为NaN避免绘制
• 递增约束：使用minimize替代curve_fit，添加不等式约束保证y<-0.3区域的偏导数非负，强制z随y减小（更负）时递增

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