4. 面形误差的Zernike多项式拟合
任务说明
本案例给定一个镜面面形误差数据文件（离散采样点上的矢高偏差，单位nm），要求使用Zernike圆域多项式对该面形进行最小二乘拟合，提取前15项（对应noll序号1～15，含平移项）系数，并利用拟合系数重构面形，输出重构残差以评估拟合质量。
数据约定
输入数据为文本文件 surface_error.txt，每行三个浮点数：归一化坐标 x、y（圆域内，x²+y²≤1）以及该点面形误差值 z（单位nm）。数据点数量不少于200，均匀分布于单位圆内。本案例附带生成该数据文件的Python脚本，确保可复现。
实现代码
import numpy as np
from numpy.polynomial.legendre import leggauss
# ---------- Zernike 圆域多项式（Noll 序号 1~15） ----------
def zernike_noll(n, m, x, y):
"""返回第 noll 序号对应的 Zernike 多项式值（向量化）"""
r = np.sqrt(x**2 + y**2)
theta = np.arctan2(y, x)
if n == 0 and m == 0: # 1: piston
return np.ones_like(r)
if n == 1 and m == -1: # 2: tilt y
return 2*r*np.sin(theta)
if n == 1 and m == 1: # 3: tilt x
return 2*r*np.cos(theta)
if n == 2 and m == 0: # 4: defocus
return np.sqrt(3)*(2*r**2 - 1)
if n == 2 and m == -2: # 5: astig 45°
return np.sqrt(6)*r**2*np.sin(2*theta)
if n == 2 and m == 2: # 6: astig 0°
return np.sqrt(6)*r**2*np.cos(2*theta)
if n == 3 and m == -1: # 7: coma y
return np.sqrt(8)*(3*r**3 - 2*r)*np.sin(theta)
if n == 3 and m == 1: # 8: coma x
return np.sqrt(8)*(3*r**3 - 2*r)*np.cos(theta)
if n == 3 and m == -3: # 9: trefoil y
return np.sqrt(8)*r**3*np.sin(3*theta)
if n == 3 and m == 3: # 10: trefoil x
return np.sqrt(8)*r**3*np.cos(3*theta)
if n == 4 and m == 0: # 11: primary spherical
return np.sqrt(5)*(6*r**4 - 6*r**2 + 1)
if n == 4 and m == -2: # 12: secondary astig y
return np.sqrt(10)*(4*r**4 - 3*r**2)*np.sin(2*theta)
if n == 4 and m == 2: # 13: secondary astig x
return np.sqrt(10)*(4*r**4 - 3*r**2)*np.cos(2*theta)
if n == 4 and m == -4: # 14: tetrafoil y
return np.sqrt(10)*r**4*np.sin(4*theta)
if n == 4 and m == 4: # 15: tetrafoil x
return np.sqrt(10)*r**4*np.cos(4*theta)
raise ValueError(f"Unsupported Noll index: n={n}, m={m}")
# Noll 序号对应的 (n, m) 列表，前15项
NOLL_LIST = [(0,0),(1,-1),(1,1),(2,0),(2,-2),(2,2),
(3,-1),(3,1),(3,-3),(3,3),(4,0),
(4,-2),(4,2),(4,-4),(4,4)]
def fit_zernike(x, y, z, terms=15):
"""最小二乘拟合，返回系数数组（顺序对应 NOLL_LIST 前 terms 项）"""
A = np.column_stack([zernike_noll(n, m, x, y) for (n, m) in NOLL_LIST[:terms]])
coeffs, _, _, _ = np.linalg.lstsq(A, z, rcond=None)
return coeffs
def reconstruct(x, y, coeffs):
"""利用系数重构面形"""
terms = len(coeffs)
A = np.column_stack([zernike_noll(n, m, x, y) for (n, m) in NOLL_LIST[:terms]])
return A @ coeffs
# ---------- 生成模拟面形数据（用于演示，实际使用时替换为读文件） ----------
def generate_demo_data(n_points=500, seed=42):
"""在单位圆内生成随机点，并叠加已知Zernike面形与噪声"""
rng = np.random.default_rng(seed)
# 均匀圆域采样（拒绝采样）
pts = []
while len(pts) < n_points:
x = rng.uniform(-1, 1, n_points*2)
y = rng.uniform(-1, 1, n_points*2)
mask = x**2 + y**2 <= 1
pts.extend(zip(x[mask], y[mask]))
pts = np.array(pts[:n_points])
x, y = pts[:,0], pts[:,1]
# 真实系数（前15项，单位nm）
true_coeffs = np.array([5.0, -3.0, 2.0, 8.0, -1.5, 4.0,
2.5, -6.0, 1.0, -2.0, 7.5,
0.5, -1.0, 0.8, -0.3])
z_true = reconstruct(x, y, true_coeffs)
noise = rng.normal(0, 0.5, n_points) # 0.5nm 噪声
z = z_true + noise
return x, y, z, true_coeffs
# ---------- 主流程 ----------
if __name__ == "__main__":
# 方式1：读取外部数据文件（取消注释使用）
# data = np.loadtxt("surface_error.txt")
# x, y, z = data[:,0], data[:,1], data[:,2]
# 方式2：生成演示数据（本案例自带，保证可复现）
x, y, z, true_coeffs = generate_demo_data()
# 拟合前15项
coeffs_fit = fit_zernike(x, y, z, terms=15)
# 重构面形与残差
z_recon = reconstruct(x, y, coeffs_fit)
residual = z - z_recon
rms_residual = np.sqrt(np.mean(residual**2))
rms_original = np.sqrt(np.mean(z**2))
# 输出结果
print("拟合系数（Noll 1~15，单位nm）：")
for i, c in enumerate(coeffs_fit, 1):
print(f" Z{i:2d}: {c:8.3f}")
print(f"\n原始面形 RMS: {rms_original:.3f} nm")
print(f"重构残差 RMS: {rms_residual:.3f} nm")
print(f"残差/原始比: {rms_residual/rms_original*100:.2f}%")
# 保存系数到文件
np.savetxt("zernike_coeffs.txt", coeffs_fit, fmt="%.6f", header="Noll 1-15 coefficients (nm)")落地与运行方式
- 将上述代码保存为
zernike_fit.py； - 若使用外部数据，将
surface_error.txt置于同目录，并取消主流程中“方式1”注释、注释掉“方式2”； - 运行
python zernike_fit.py，依赖仅需numpy。
验证方法与判定标准
- 程序输出15个拟合系数、原始RMS、残差RMS及百分比；
- 使用自带演示数据时，拟合系数与真实系数（
true_coeffs）偏差应小于0.1nm（因噪声0.5nm，15项拟合可良好恢复）； - 残差RMS应明显小于原始RMS（通常低于10%），表明前15项已捕获主要面形成分；
- 若使用外部数据，观察残差RMS是否达到预期拟合精度，并检查
zernike_coeffs.txt是否成功生成。