2. 自由曲面面型参数化建模入门

任务说明

本案例用“径向基函数（RBF）+ 多项式”混合模型对自由曲面面型进行参数化建模：给定一组离散采样点 (x, y, z)，拟合出连续矢高函数 z(x, y)，并输出矢高图与斜率图。核心是构造混合基矩阵并求解最小二乘系数。

环境与依赖

• Python 3.8+，依赖：numpy、matplotlib
• 安装：pip install numpy matplotlib

模型与实现

混合模型：z(x,y) = Σ wᵢ·φ(rᵢ) + (a₀ + a₁x + a₂y + a₃x² + a₄xy + a₅y²)，其中 φ(r)=exp(-(εr)²) 为高斯 RBF，中心取采样点。未知量为 RBF 权重 w 与多项式系数 a，通过最小二乘一次求解。斜率由解析偏导得到。

完整代码

import numpy as np
import matplotlib.pyplot as plt

# 1. 生成自由曲面采样点（含轻微非对称，模拟真实面型）
rng = np.random.default_rng(0)
n = 120
x = rng.uniform(-1, 1, n)
y = rng.uniform(-1, 1, n)
z = (0.8*x**2 + 0.5*y**2 + 0.3*x*y
     + 0.15*np.exp(-((x-0.4)**2 + (y+0.3)**2)/0.05)
     + 0.01*rng.standard_normal(n))

# 2. 构造混合基设计矩阵
eps = 1.2  # RBF 形状参数
def rbf_matrix(x, y, cx, cy, eps):
    r2 = (x[:, None]-cx[None, :])**2 + (y[:, None]-cy[None, :])**2
    return np.exp(-(eps**2) * r2)

def poly_matrix(x, y):
    return np.column_stack([np.ones_like(x), x, y, x**2, x*y, y**2])

A = np.hstack([rbf_matrix(x, y, x, y, eps), poly_matrix(x, y)])
coef, *_ = np.linalg.lstsq(A, z, rcond=None)
w, a = coef[:n], coef[n:]

# 3. 面型与斜率函数（解析偏导）
def surface(xq, yq):
    R = rbf_matrix(xq, yq, x, y, eps)
    P = poly_matrix(xq, yq)
    return R @ w + P @ a

def slopes(xq, yq):
    dx = xq[:, None]-x[None, :]
    dy = yq[:, None]-y[None, :]
    R = np.exp(-(eps**2)*(dx**2+dy**2))
    dzdx = (-2*eps**2*dx*R) @ w + (a[1] + 2*a[3]*xq + a[4]*yq)
    dzdy = (-2*eps**2*dy*R) @ w + (a[2] + a[4]*xq + 2*a[5]*yq)
    return dzdx, dzdy

# 4. 网格评估并绘图
g = np.linspace(-1, 1, 80)
X, Y = np.meshgrid(g, g)
xf, yf = X.ravel(), Y.ravel()
Z = surface(xf, yf).reshape(X.shape)
SX, SY = slopes(xf, yf)
S = np.sqrt(SX**2 + SY**2).reshape(X.shape)

fig, ax = plt.subplots(1, 2, figsize=(11, 4.5))
c0 = ax[0].contourf(X, Y, Z, 30, cmap="viridis")
ax[0].set_title("Sag map z(x,y)")
fig.colorbar(c0, ax=ax[0])
c1 = ax[1].contourf(X, Y, S, 30, cmap="magma")
ax[1].set_title("Slope magnitude |grad z|")
fig.colorbar(c1, ax=ax[1])
plt.tight_layout()
plt.savefig("freeform_sag_slope.png", dpi=120)

# 5. 拟合残差验证
res = z - surface(x, y)
print("RMS residual:", np.sqrt(np.mean(res**2)))

运行与验证

• 运行：python freeform_model.py，生成 freeform_sag_slope.png。
• 通过条件：左图矢高呈平滑连续曲面，右图斜率幅值分布合理；终端打印的 RMS 残差应接近采样噪声量级（约 0.01），说明混合模型能拟合面型。
• 若残差明显偏大，可减小 eps 或增加采样点数量。