2. Zernike多项式拟合自由曲面

任务定义

本案例给定一个自由曲面（非球面基底叠加正弦波纹扰动），目标为：使用标准Zernike圆域多项式对该面形进行最小二乘拟合，得到展开系数，并利用系数重构面形，最终通过残差RMS与PV值评估拟合精度。输入为离散采样网格上的矢高数据，输出为Zernike系数、重构面形及精度指标。

实现方案

采用极坐标下的标准Zernike多项式（Noll序号，前36项，含活塞、倾斜、离焦及高阶像差）。将直角坐标采样点映射至单位圆内，构建设计矩阵，通过最小二乘求解系数。重构面形由系数线性组合得到，残差为原始矢高与重构值之差。精度指标取残差面形的RMS（均方根）与PV（峰谷值）。

代码实现

import numpy as np
from numpy.polynomial.legendre import leggauss

def zernike_poly(n, m, rho, theta):
    """标准Zernike多项式（Noll序号对应n,m）"""
    if m > 0:
        R = np.zeros_like(rho)
        for s in range((n - m)//2 + 1):
            coeff = (-1)**s * np.math.factorial(n - s)
            coeff /= (np.math.factorial(s) * 
                      np.math.factorial((n + m)//2 - s) * 
                      np.math.factorial((n - m)//2 - s))
            R += coeff * rho**(n - 2*s)
        return R * np.cos(m * theta)
    elif m < 0:
        m_abs = abs(m)
        R = np.zeros_like(rho)
        for s in range((n - m_abs)//2 + 1):
            coeff = (-1)**s * np.math.factorial(n - s)
            coeff /= (np.math.factorial(s) * 
                      np.math.factorial((n + m_abs)//2 - s) * 
                      np.math.factorial((n - m_abs)//2 - s))
            R += coeff * rho**(n - 2*s)
        return R * np.sin(m_abs * theta)
    else:
        R = np.zeros_like(rho)
        for s in range(n//2 + 1):
            coeff = (-1)**s * np.math.factorial(n - s)
            coeff /= (np.math.factorial(s) * np.math.factorial(n//2 - s)**2)
            R += coeff * rho**(n - 2*s)
        return R

def noll_to_nm(j):
    """Noll序号转(n,m)"""
    n = 0
    while (n + 1)*(n + 2)//2 < j:
        n += 1
    m = 0
    if (j - n*(n+1)//2) % 2 == 0:
        m = 2 * ((j - n*(n+1)//2) // 2)
    else:
        m = 2 * ((j - n*(n+1)//2) // 2) + 1
    if m > n:
        m = n - (m - n)
    # 确定符号
    if m == 0:
        return n, 0
    elif (j - n*(n+1)//2) % 2 == 0:
        return n, m
    else:
        return n, -m

def fit_zernike(x, y, z, max_order=36):
    """最小二乘拟合Zernike系数"""
    # 归一化到单位圆
    r_max = np.sqrt(x**2 + y**2).max()
    x_n = x / r_max
    y_n = y / r_max
    rho = np.sqrt(x_n**2 + y_n**2)
    theta = np.arctan2(y_n, x_n)
    
    # 构建设计矩阵
    A = np.zeros((len(z), max_order))
    for j in range(1, max_order + 1):
        n, m = noll_to_nm(j)
        A[:, j-1] = zernike_poly(n, m, rho, theta)
    
    # 最小二乘求解
    coeffs, _, _, _ = np.linalg.lstsq(A, z, rcond=None)
    z_fit = A @ coeffs
    residual = z - z_fit
    
    rms = np.sqrt(np.mean(residual**2))
    pv = residual.max() - residual.min()
    return coeffs, z_fit, rms, pv

# ============ 测试与验证 ============
if __name__ == "__main__":
    # 生成测试面形：已知系数组合 + 高频波纹
    N = 200
    x = np.linspace(-1, 1, N)
    X, Y = np.meshgrid(x, x)
    mask = X**2 + Y**2 <= 1.0
    Xm, Ym = X[mask], Y[mask]
    
    # 真实面形：离焦(0.5) + 像散(0.3) + 彗差(0.2) + 波纹扰动
    rho = np.sqrt(Xm**2 + Ym**2)
    theta = np.arctan2(Ym, Xm)
    z_true = (0.5 * zernike_poly(2, 0, rho, theta) +
              0.3 * zernike_poly(2, 2, rho, theta) +
              0.2 * zernike_poly(3, 1, rho, theta) +
              0.05 * np.sin(20 * rho) * np.cos(15 * theta))
    
    # 拟合
    coeffs, z_fit, rms, pv = fit_zernike(Xm, Ym, z_true, max_order=36)
    
    # 输出结果
    print("=== 拟合精度 ===")
    print(f"残差RMS: {rms:.6f} mm")
    print(f"残差PV:  {pv:.6f} mm")
    print(f"面形RMS: {np.sqrt(np.mean(z_true**2)):.6f} mm")
    print(f"RMS比:   {rms/np.sqrt(np.mean(z_true**2))*100:.3f}%")
    
    print("\n=== 前10项Zernike系数 ===")
    names = ['活塞', 'X倾斜', 'Y倾斜', '离焦', '斜像散', '垂直像散',
             'X彗差', 'Y彗差', '三叶草', '三叶草']
    for i in range(10):
        print(f"Z{i+1:2d} ({names[i]:4s}): {coeffs[i]:+.6f}")

运行与验证

依赖：Python 3.8+，NumPy（无其他第三方库）。

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

预期输出与判定标准：

• 程序打印残差RMS、PV及面形RMS。由于波纹频率（20周期/半径）超出36项Zernike的拟合能力，残差RMS应约为面形RMS的1%～5%（即RMS比在1～5%之间），PV约为0.1～0.3 mm量级——这验证了截断误差的存在。
• 前几项系数中，Z4（离焦）应接近0.5，Z5/Z6（像散）应接近0.3，Z7/Z8（彗差）应接近0.2，偏差小于1e-6，证明拟合算法正确恢复了低频成分。
• 若将波纹幅度设为0（即纯Zernike面形），残差RMS应接近机器精度（<1e-10），此时可验证算法本身的数值精度。

判定通过条件：低频系数恢复误差<1e-5，且含波纹时残差RMS比在1%～5%区间内，即证明拟合与精度分析流程有效。