3. 自由曲面斜率计算与曲率分析

任务定义

本案例针对自由曲面光学设计中面形质量评估的核心环节：给定一组离散面形数据（x, y, z），计算各点处的斜率（一阶导数）与曲率（二阶导数），并可视化曲率分布以识别局部面形异常。输入数据采用旋转对称非球面叠加正弦波纹的解析面形，便于与理论值对比验证算法正确性。

数值方法

采用五点中心差分格式计算偏导数。对于均匀网格间距 h，一阶偏导为：

∂z/∂x ≈ (-z(x+2h) + 8z(x+h) - 8z(x-h) + z(x-2h)) / (12h)

二阶偏导及混合偏导类似。平均曲率 H 和高斯曲率 K 由下式得到：

H = ((1+zx²)zyy - 2zx·zy·zxy + (1+zy²)zxx) / (2(1+zx²+zy²)^(3/2))
K = (zxx·zyy - zxy²) / (1+zx²+zy²)²

边界点采用单侧差分，内部点全部使用中心差分以保证精度。

自由曲面斜率与曲率计算流程 输入面形 z = f(x,y) 离散网格 五点中心差分 zx, zy, zxx, zyy, zxy 边界：单侧差分 曲率合成 平均曲率 H 高斯曲率 K 可视化 曲率分布 验证路径 1. 解析面形：旋转对称非球面 + 正弦波纹 2. 理论斜率/曲率：对解析式求导 3. 对比数值解与理论解，计算最大误差 判定标准 • 内部点斜率误差 < 1e-6（双精度浮点下五点差分精度 O(h⁴)） • 曲率误差 < 1e-4（二阶差分精度 O(h²)） • 曲率分布热力图可清晰分辨波纹周期与幅值变化

代码实现


import numpy as np
import matplotlib.pyplot as plt

# ========== 1. 解析面形定义 ==========
def surface(x, y):
    """旋转对称非球面(二次曲面) + 正弦波纹扰动"""
    c = 0.02          # 曲率半径倒数 (1/mm)
    k = -0.5          # 圆锥系数
    r2 = x**2 + y**2
    # 非球面 sagitta
    z_base = c * r2 / (1 + np.sqrt(1 - (1+k) * c**2 * r2))
    # 正弦波纹: 幅值 2μm, 空间频率 0.5/mm
    z_ripple = 0.002 * np.sin(2 * np.pi * 0.5 * x) * np.cos(2 * np.pi * 0.5 * y)
    return z_base + z_ripple

# ========== 2. 理论导数（解析求导） ==========
def surface_derivatives(x, y):
    """返回 (zx, zy, zxx, zyy, zxy) 理论值"""
    r2 = x**2 + y**2
    c, k = 0.02, -0.5
    denom = np.sqrt(1 - (1+k) * c**2 * r2)
    # 非球面部分一阶导
    zx_base = c * x / denom
    zy_base = c * y / denom
    # 非球面部分二阶导
    zxx_base = c / denom + (1+k) * c**3 * x**2 / denom**3
    zyy_base = c / denom + (1+k) * c**3 * y**2 / denom**3
    zxy_base = (1+k) * c**3 * x * y / denom**3
    # 波纹部分
    wx, wy = 2*np.pi*0.5, 2*np.pi*0.5
    A = 0.002
    zx_rip = A * wx * np.cos(wx*x) * np.cos(wy*y)
    zy_rip = -A * wy * np.sin(wx*x) * np.sin(wy*y)
    zxx_rip = -A * wx**2 * np.sin(wx*x) * np.cos(wy*y)
    zyy_rip = -A * wy**2 * np.sin(wx*x) * np.cos(wy*y)
    zxy_rip = -A * wx * wy * np.cos(wx*x) * np.sin(wy*y)
    return (zx_base+zx_rip, zy_base+zy_rip, zxx_base+zxx_rip,
            zyy_base+zyy_rip, zxy_base+zxy_rip)

# ========== 3. 数值差分计算 ==========
def numerical_derivatives(Z, h):
    """五点中心差分计算一阶/二阶偏导，边界用单侧差分"""
    ny, nx = Z.shape
    zx = np.zeros_like(Z); zy = np.zeros_like(Z)
    zxx = np.zeros_like(Z); zyy = np.zeros_like(Z); zxy = np.zeros_like(Z)

    # 内部点（五点中心差分）
    for i in range(2, ny-2):
        for j in range(2, nx-2):
            zx[i,j] = (-Z[i,j+2] + 8*Z[i,j+1] - 8*Z[i,j-1] + Z[i,j-2]) / (12*h)
            zy[i,j] = (-Z[i+2,j] + 8*Z[i+1,j] - 8*Z[i-1,j] + Z[i-2,j]) / (12*h)
            zxx[i,j] = (-Z[i,j+2] + 16*Z[i,j+1] - 30*Z[i,j] + 16*Z[i,j-1] - Z[i,j-2]) / (12*h**2)
            zyy[i,j] = (-Z[i+2,j] + 16*Z[i+1,j] - 30*Z[i,j] + 16*Z[i-1,j] - Z[i-2,j]) / (12*h**2)
            zxy[i,j] = (Z[i+1,j+1] - Z[i+1,j-1] - Z[i-1,j+1] + Z[i-1,j-1]) / (4*h**2)

    # 边界点（单侧差分，二阶精度）
    # 上下边界 (i=0,1,ny-2,ny-1)
    for i in [0, 1, ny-2, ny-1]:
        for j in range(2, nx-2):
            if i == 0:
                zy[i,j] = (-3*Z[i,j] + 4*Z[i+1,j] - Z[i+2,j]) / (2*h)
                zyy[i,j] = (2*Z[i,j] - 5*Z[i+1,j] + 4*Z[i+2,j] - Z[i+3,j]) / h**2
            elif i == 1:
                zy[i,j] = (Z[i+1,j] - Z[i-1,j]) / (2*h)
                zyy[i,j] = (Z[i+1,j] - 2*Z[i,j] + Z[i-1,j]) / h**2
            elif i == ny-2:
                zy[i,j] = (Z[i+1,j] - Z[i-1,j]) / (2*h)
                zyy[i,j] = (Z[i+1,j] - 2*Z[i,j] + Z[i-1,j]) / h**2
            else:  # i == ny-1
                zy[i,j] = (3*Z[i,j] - 4*Z[i-1,j] + Z[i-2,j]) / (2*h)
                zyy[i,j] = (2*Z[i,j] - 5*Z[i-1,j] + 4*Z[i-2,j] - Z[i-3,j]) / h**2
            # 边界处 x 方向仍用中心差分
            zx[i,j] = (-Z[i,j+2] + 8*Z[i,j+1] - 8*Z[i,j-1] + Z[i,j-2]) / (12*h)
            zxx[i,j] = (-Z[i,j+2] + 16*Z[i,j+1] - 30*Z[i,j] + 16*Z[i,j-1] - Z[i,j-2]) / (12*h**2)
            zxy[i,j] = (Z[i+1,j+1] - Z[i+1,j-1] - Z[i-1,j+1] + Z[i-1,j-1]) / (4*h**2)

    # 左右边界 (j=0,1,nx-2,nx-1)，内部 i 范围
    for j in [0, 1, nx-2, nx-1]:
        for i in range(2, ny-2):
            if j == 0:
                zx[i,j] = (-3*Z[i,j] + 4*Z[i,j+1] - Z[i,j+2]) / (2*h)
                zxx[i,j] = (2*Z[i,j] - 5*Z[i,j+1] + 4*Z[i,j+2] - Z[i,j+3]) / h**2
            elif j == 1:
                zx[i,j] = (Z[i,j+1] - Z[i,j-1]) / (2*h)
                zxx[i,j] = (Z[i,j+1] - 2*Z[i,j] + Z[i,j-1]) / h**2
            elif j == nx-2:
                zx[i,j] = (Z[i,j+1] - Z[i,j-1]) / (2*h)
                zxx[i,j] = (Z[i,j+1] - 2*Z[i,j] + Z[i,j-1]) / h**2
            else:  # j == nx-1
                zx[i,j] = (3*Z[i,j] - 4*Z[i,j-1] + Z[i,j-2]) / (2*h)
                zxx[i,j] = (2*Z[i,j] - 5*Z[i,j-1] + 4*Z[i,j-2] - Z[i,j-3]) / h**2
            zy[i,j] = (-Z[i+2,j] + 8*Z[i+1,j] - 8*Z[i-1,j] + Z[i-2,j]) / (12*h)
            zyy[i,j] = (-Z[i+2,j] + 16*Z[i+1,j] - 30*Z[i,j] + 16*Z[i-1,j] - Z[i-2,j]) / (12*h**2)
            zxy[i,j] = (Z[i+1,j+1] - Z[i+1,j-1] - Z[i-1,j+1] + Z[i-1,j-1]) / (4*h**2)

    # 四个角点（双线性外推近似）
    for (i,j) in [(0,0),(0,nx-1),(ny-1,0),(ny-1,nx-1)]:
        zx[i,j] = zx[i, 2] if j == 0 else zx[i, nx-3]
        zy[i,j] = zy[2, j] if i == 0 else zy[ny-3, j]
        zxx[i,j] = zxx[i, 2] if j == 0 else zxx[i, nx-3]
        zyy[i,j] = zyy[2, j] if i == 0 else zyy[ny-3, j]
        zxy[i,j] = 0.0

    return zx, zy, zxx, zyy, zxy

# ========== 4. 曲率计算 ==========
def compute_curvature(zx, zy, zxx, zyy, zxy):
    """计算平均曲率 H 和高斯曲率 K"""
    E = 1 + zx**2 + zy**2
    H = ((1+zx**2)*zyy - 2*zx*zy*zxy + (1+zy**2)*zxx) / (2 * E**1.5)
    K = (zxx*zyy - zxy**2) / E**2
    return H, K

# ========== 5. 主程序 ==========
# 网格参数
L = 10.0          # 口径 10mm × 10mm
N = 101           # 网格点数
h = L / (N-1)     # 网格间距
x = np.linspace(-L/2, L/2, N)
y = np.linspace(-L/2, L/2, N)
X, Y = np.meshgrid(x, y)

# 采样面形
Z = surface(X, Y)

# 数值导数
zx_num, zy_num, zxx_num, zyy_num, zxy_num = numerical_derivatives(Z, h)

# 数值曲率
H_num, K_num = compute_curvature(zx_num, zy_num, zxx_num, zyy_num, zxy_num)

# 理论导数与曲率（用于验证）
zx_theo, zy_theo, zxx_theo, zyy_theo, zxy_theo = surface_derivatives(X, Y)
H_theo, K_theo = compute_curvature(zx_theo, zy_theo, zxx_theo, zyy_theo, zxy_theo)

# ========== 6. 误差评估 ==========
# 排除边界两圈（边界差分精度较低）
inner = slice(2, -2)
err_zx = np.max(np.abs(zx_num[inner,inner] - zx_theo[inner,inner]))
err_zy = np.max(np.abs(zy_num[inner,inner] - zy_theo[inner,inner]))
err_H = np.max(np.abs(H_num[inner,inner] - H_theo[inner,inner]))
err_K = np.max(np.abs(K_num[inner,inner] - K_theo[inner,inner]))

print(f"内部区域最大斜率误差 zx: {err_zx:.2e}")
print(f"内部区域最大斜率误差 zy: {err_zy:.2e}")
print(f"内部区域最大平均曲率误差 H: {err_H:.2e}")
print(f"内部区域最大高斯曲率误差 K: {err_K:.2e}")

# ========== 7. 可视化 ==========
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# 面形
im0 = axes[0,0].imshow(Z, extent=[-L/2,L/2,-L/2,L/2], origin='lower', cmap='viridis')
axes[0,0].set_title('Surface Sagitta z(x,y) [mm]')
plt.colorbar(im0, ax=axes[0,0])

# 平均曲率
im1 = axes[0,1].imshow(H_num, extent=[-L/2,L/2,-L/2,L/2], origin='lower', cmap='RdBu_r')
axes[0,1].set_title('Mean Curvature H [1/mm]')
plt.colorbar(im1, ax=axes[0,1])

# 高斯曲率
im2 = axes[1,0].imshow(K_num, extent=[-L/2,L/2,-L/2,L/2], origin='lower', cmap='RdBu_r')
axes[1,0].set_title('Gaussian Curvature K [1/mm²]')
plt.colorbar(im2, ax=axes[1,0])

# 曲率误差分布
im3 = axes[1,1].imshow(np.abs(H_num - H_theo), extent=[-L/2,L/2,-L/2,L/2],
                       origin='lower', cmap='hot')
axes[1,1].set_title('|ΔH| Numerical vs Theoretical')
plt.colorbar(im3, ax=axes[1,1])

plt.tight_layout()
plt.savefig('curvature_analysis.png', dpi=150)
plt.show()

运行与验证

依赖环境：Python 3.8+，NumPy，Matplotlib。运行前确认已安装：pip install numpy matplotlib。

执行方式：将上述代码保存为 curvature_analysis.py，在终端运行 python curvature_analysis.py。

预期输出：程序打印内部区域（排除边界两圈）的数值误差，并生成 curvature_analysis.png 四联图。由于网格间距 h=0.1mm，五点中心差分截断误差为 O(h⁴)=1e-4 量级，二阶差分误差为 O(h²)=1e-2 量级，因此预期误差范围：斜率误差约 1e-6~1e-5，曲率误差约 1e-4~1e-3。具体数值取决于浮点舍入与边界处理，但应满足：斜率误差 < 1e-4，曲率误差 < 1e-2。

判定标准：

  • 程序无异常退出，四联图正常显示并保存；
  • 打印的斜率误差与曲率误差在预期量级范围内；
  • 平均曲率图中可清晰看到正弦波纹引起的周期性曲率变化（波纹频率 0.5/mm，对应图上约 5 个周期）；
  • 误差分布图显示最大误差集中在边界附近，内部区域颜色接近黑色（误差极小），说明中心差分精度良好。

若误差显著偏大（如斜率误差 > 1e-3），检查网格间距 h 是否与 x/y 数组一致，或边界处理代码是否存在索引越界。

← 上一章 📖 返回目录 下一章 →