尧图网页设计 
简介Zernike多项式是光学波前分析与像差校正中的常用数学工具借助Matlab进行曲面拟合可大幅提升计算效率。这份代码面向计算机、电子信息工程、数学等专业的大学生及科研人员适用于课程设计、期末大作业、毕业设计或光学仿真任务。资源包仅含2个文件主要为1个.m函数文件和1个txt说明文件压缩包大小约3KB结构精简、定位明确。代码采用参数化编程参数可灵活调整思路清晰并附有详细注释同时提供可直接运行的案例数据兼容Matlab 2014至2024a版本。使用者无需从零搭建拟合流程只需修改参数或导入数据即可快速获得Zernike拟合结果并观察曲面重构效果有助于理解多项式拟合原理、缩短调试时间。目前已有204人学习下载适合刚接触Zernike拟合或希望在项目中快速落地该算法的读者。1. Zernike 多项式拟合曲面镜面检测与像差分析都在用的那套函数你手里有一台干涉仪或者共聚焦显微镜的测量结果几千到几十万个 (x, y, z) 散点描述一块镜面、一片晶圆或者一个光学元件的表面起伏。要把这批数据变成可对比、可传递、可反推的结果最通用的做法就是把面形投影到 Zernike 多项式上。这个标题里的 .rar 在光学圈流传了很多年解开后通常就是一个或一组 .m 文件输入离散坐标与面形高度输出一组 Zernike 系数。反直觉的地方在于拟合得到的几十个系数比原始云图更接近“答案”。离焦 0.8 µm、水平散光 0.2 µm——这些物理量直接从系数里读出来而不是从云图里量出来。原因是 Zernike 基函数在单位圆域上正交每个系数只对应一种像差模式互不干扰。正因如此从望远镜主镜检测到眼科波前像差仪比如 iTrace的数据交换都拿 Zernike 系数当标准语言。这篇文章按我平时给光学和量测方向同事梳理的思路展开先讲清 Zernike 多项式的定义、排序与归一化再给出一个可以直接抄进项目的 Matlab 拟合函数然后处理归一化半径、阶数和数据预处理这些最容易翻车的参数最后落到与像差仪数据对接和波前重建的验证技巧上。新手可以一步步跑通熟手可以直接跳到 3.2 和 4.2 对照自己的实现。2. Zernike 多项式拟合的理论基础正交性、排序与单位圆坐标2.1 定义式与两种主流排序ANSI 序和 Noll 序Zernike 多项式定义在单位圆盘上用极坐标 (ρ, θ) 表示ρ ∈ [0,1]θ ∈ [0,2π)。每一项由径向多项式与角向函数相乘Z_n^m(ρ, θ) N_n^m · R_n^{|m|}(ρ) · {cos(mθ), sin(mθ)}其中径向多项式 R_n^{|m|}(ρ) 的闭合形式是R_n^{|m|}(ρ) Σ_{s0}^{(n−|m|)/2} (−1)^s (n−s)! / [s! ((n|m|)/2 − s)! ((n−|m|)/2 − s)!] · ρ^{n−2s}约束条件n 与 |m| 同奇偶且 |m| ≤ n。这里的 n 是径向阶数|m| 是角向频率。归一化系数 N_n^m 按 ANSI Z80.28 标准取m0 时为 √(n1)m≠0 时为 √(2(n1))。这样每项在单位圆上的均方值为 1拟合得到的系数直接等于该项对波前方差的贡献单位与输入面形数据的单位一致。排序是个大坑。同一套 Zernike 多项式ANSI 序眼科行业常用、Noll 序天文望远镜常用和 Born Wolf 序干涉仪厂家常用的编号完全不同。我推荐的实现方式是显式维护一个 (n, m) 项表让系数含义自解释系数位置nm物理名称表达式100平移piston121−1y 方向倾斜2ρ sinθ311x 方向倾斜2ρ cosθ42−245° 散光√6 ρ² sin2θ520离焦defocus√3(2ρ²−1)6220° 散光√6 ρ² cos2θ73−1垂直彗差√8(3ρ³−2ρ) sinθ831水平彗差√8(3ρ³−2ρ) cosθ93−3三叶草30°√8 ρ³ sin3θ1033三叶草0°√8 ρ³ cos3θ提示表格里 m 为负时用 sin、为正时用 cos这只是惯用约定。iTrace、Zygo 等设备可能用相反符号跨设备对比系数前必须先验证方向定义。2.2 为什么必须在圆域上拟合正交性与“核函数”视角普通多项式比如幂级数 x^i y^j在矩形域上做拟合各项之间存在强相关性稍微加点噪声系数就乱跳。Zernike 多项式在连续单位圆盘上满足严格正交性∫∫ Z_i(x,y) · Z_j(x,y) dxdy π · δ_ij也就是说把任意光滑圆孔面形展开成 Zernike 级数时每一项的能量独立分离。这是它成为光学面形“标准语言”的根本原因。但要注意离散测量点不满足连续正交条件。当数据点在圆盘内近似均匀分布时设计矩阵 Z每列是一个基函数在所有数据点上的取值的各列近似正交当数据稀疏、偏心或者有大块缺失时列间会出现相关性。从核函数的角度看最小二乘解依赖 Gram 矩阵 K ZᵀZK 的对角占优程度决定了拟合稳定度。K 的条件数是 Z 的条件数的平方这就是为什么后面要用 QR 分解而不是直接解正规方程。另一个容易被忽略的点Zernike 多项式只在归一化单位圆内定义拟合前必须把物理坐标映射到单位圆并且把圆外数据彻底剔除。映射关系是 ρ √(x²y²)/Rθ atan2(y, x)其中 R 是归一化半径量纲与 x、y 一致。2.3 坐标变换与掩膜的最小实现这段代码是所有 Zernike 拟合函数的第一步先单独拿出来验证% x, y : 列向量物理坐标单位任意但必须一致 % R : 归一化半径通常取有效孔径半径 rho sqrt(x.^2 y.^2) / R; theta atan2(y, x); mask rho 1 1e-12; % 圆内有效点1e-12 是浮点容差逻辑说明ρ 超过 1 的点不在单位圆定义域内必须剔除这里的 mask 同时服务于两项任务——筛选有效数据、保证后续设计矩阵的行数正确。θ 必须用 atan2 而不是 atan(y/x)因为 atan 会把二、三象限的点误判到一、四象限导致倾斜项和彗差项的符号全部反掉。1e-12 的容差是为了把恰好落在边界上的点保留下来这种点通常来自网格采样时 ρ 恰好等于 1 的情况。3. 用 Matlab 写一个 Zernike 曲面拟合函数函数声明、设计矩阵与最小二乘3.1 函数声明输入输出与单位约定常见的 Matlab 拟合函数开头长这样文件名必须与函数名一致保存为zernikeFitSrf.mfunction [coeff, terms, zern, stats] zernikeFitSrf(x, y, z, R, N) % ZERNIKEFITSRF 用 Zernike 多项式拟合圆域面形数据 % 输入 % x, y, z : 面形采样点的坐标与高度同为列向量单位一致 % R : 归一化半径与 x, y 同单位 % N : 需要拟合的 Zernike 项数 % 输出 % coeff : N x 1 系数向量单位与 z 相同 % terms : N x 2 矩阵每行是该项的 (n, m) % zern : 拟合回代值只含有效圆域内的点 % stats : 结构体含残差 RMS、PV 和设计矩阵条件数这里要明确一个单位约定z 用 µm、nm 还是 mm由你的测量设备决定系数单位自动跟随 z。x、y 用 mmR 就用 mm三者保持同一套长度单位即可Zernike 基函数本身是无量纲的。关于 Matlab 函数声明还有一条规则一个 .m 文件里只有第一个函数对外可见其余都是子函数所以radialZernike、zernikeTerm这些辅助函数要么放同一个文件里要么各自建文件。3.2 径向多项式、单顶生成与设计矩阵构造先实现径向多项式。用阶乘公式直接求中等阶数n ≤ 30完全没问题function R radialZernike(n, m, rho) % 计算径向多项式 R_n^|m|(rho)rho 可以是向量 m abs(m); R zeros(size(rho)); for s 0:(n - m)/2 R R (-1)^s * factorial(n - s) ./ ... (factorial(s) * factorial((n m)/2 - s) * factorial((n - m)/2 - s)) ... .* rho.^(n - 2*s); end end逻辑说明求和上限 (n−m)/2 由 n 与 m 同奇偶保证为整数每一项的指数 n−2s 从 n 递减到 m因此 ρ0 处的值只在 m0 时非零这与 piston 项的物理意义一致。阶乘形式在 n 超过约 30 时可能引入较大舍入误差但一般面形拟合用到 7 阶36 项以内这个实现足够。接着是单顶生成与项表生成function Z zernikeTerm(n, m, rho, theta) % 生成第 (n, m) 项 Zernike 基函数按 ANSI 归一化 R radialZernike(n, m, rho); if m 0 Z sqrt(2*(n 1)) * R .* cos(m*theta); elseif m 0 Z sqrt(2*(n 1)) * R .* sin(abs(m)*theta); else Z sqrt(n 1) * R; end end function terms zernikeTerms(N) % 按 n 升序、(m -n:2:n) 生成前 N 项 (n, m) 表 terms zeros(N, 2); k 0; n 0; while k N for m -n:2:n k k 1; if k N break; end terms(k, :) [n, m]; end n n 1; end end系数 √(n1) 与 √(2(n1)) 就是 2.1 节说的 ANSI 归一化因子去掉它们也能拟合但系数含义会变成“未归一化基下的系数”无法直接和其他设备对比。项表生成用 while 循环按 n 逐层展开保证任意 N 都有明确的 (n, m) 对应关系。主函数把这些拼起来function [coeff, terms, zern, stats] zernikeFitSrf(x, y, z, R, N) assert(isequal(size(x), size(y), size(z)), x/y/z 必须同尺寸); x x(:); y y(:); z z(:); % 统一为列向量 rho sqrt(x.^2 y.^2) / R; theta atan2(y, x); mask rho 1 1e-12; % 有效圆域掩膜 if nnz(mask) N error(有效数据点 (%d) 少于拟合项数 (%d)请增大 R 或减小 N, nnz(mask), N); end terms zernikeTerms(N); Zmat zeros(nnz(mask), N); for k 1:N Zmat(:, k) zernikeTerm(terms(k,1), terms(k,2), ... rho(mask), theta(mask)); end coeff Zmat \ z(mask); % 直接最小二乘QR 分解 zern Zmat * coeff; % 回代重构 resid z(mask) - zern; stats struct(rmse, sqrt(mean(resid.^2)), ... pv, max(resid) - min(resid), ... cond, cond(Zmat)); end参数说明Zmat \ z(mask)用的是 Matlab 反斜杠运算符内部走 QR 分解加列主元比显式计算(Z*Z)\(Z*z)数值稳定得多因为正规方程会把条件数平方。stats.cond是诊断用的关键值后面 4.2 节会专门讲。如果有效点数小于项数最小二乘就没有唯一解这里直接报错而不是让结果静默出错。3.3 用已知系数的模拟面形验证函数没有写错拿到新写的拟合函数先别急着喂实测数据。构造一个已知系数组合的仿真面形跑一遍看能不能把系数原样拿回来Npix 401; [xg, yg] meshgrid(linspace(-1, 1, Npix)); rho sqrt(xg.^2 yg.^2); theta atan2(yg, xg); mask rho 1; % 构造0.5 平移 0.8 离焦 0.2 水平散光单位自定义 zs 0.5 0.8*sqrt(3)*(2*rho.^2 - 1) 0.2*sqrt(6)*rho.^2 .* cos(2*theta); zs(~mask) NaN; idx find(mask); [coeff, terms, zern, stats] zernikeFitSrf(xg(idx), yg(idx), zs(idx), 1, 15); disp([terms, coeff]); % 检查 (n, m) 与系数的对应 fprintf(RMSE %.3e, cond %.3e\n, stats.rmse, stats.cond);按 2.1 节的项表期望结果是 coeff(1) ≈ 0.5piston、coeff(5) ≈ 0.8离焦、coeff(6) ≈ 0.2水平散光。如果这三个位置对得上其他位置接近 10⁻¹⁰ 量级RMSE 接近机器精度说明径向多项式、归一化因子和拟合主流程都没问题。对不上时按三个方向查先查 R 与 mask 是否把数据圈对再查zernikeTerm里 m 的符号约定最后查zernikeTerms的排列是否与你的预期一致。4. Zernike 拟合参数怎么设半径、阶数与数据预处理4.1 归一化半径 R 选错是头号问题R 的定义是“物理坐标除以 R 后落在单位圆内”所以 R 必须等于有效孔径半径。常见的错误有两类。第一类R 取小了圆外本应参与拟合的点被 mask 砍掉真实面形的高阶成分会泄漏到低阶系数里最典型的表现是拟合残差沿着孔径边缘出现一条亮环。第二类R 取大了拟合区域只剩中心一小块系数描述的只是局部面形而且中心区域数据密度高会把边缘的真实形状完全稀释掉。更隐蔽的问题是不同的数据集用了不同的 R。改变了 R 等于改变了基函数本身所有系数都会变不存在一个简单的比例换算关系。所以同一批对比实验R 必须固定。对于干涉仪数据R 取镜头设计孔径半径对于轮廓仪扫描数据R 取有效数据点的最大半径或者略小一点保证圆内有足够多的点参与拟合。R_fit max(sqrt(x.^2 y.^2)) * 0.999; % 按数据边界自适应取 R这里乘 0.999 而不是直接用最大值是为了避免恰好落在边界上的散点因为浮点误差被 mask 误杀。确定 R 后把它作为常量传进函数并在输出结果里记录这是工程上防止系数误读的最低要求。4.2 阶数 N 与条件数能用 QR 就不要用正规方程拟合项数 N 决定最高阶 n。取太少会欠拟合面形里的真实像差被低阶项分摊取太多则高阶项开始拟合噪声。判断标准不是“看起来平滑”而是看残差有没有明显结构以及新增系数是否发散。条件数是更硬的指标。设计矩阵 Z 的条件数随 N 增长而用正规方程求解时实际参与运算的是 Gram 矩阵 ZᵀZ它的条件数是 cond(Z)²。以下是我在均匀网格和带缺口数据上看到的典型量级不同采样差异很大关键是自己会读 stats.cond项数 N最高阶 n均匀圆盘网格 cond(Z)稀疏/扇形缺口 cond(Z)153约 1–1010–10036710–10010³–10⁴661010²–10³典型大于 10⁵提示cond 超过 10⁶ 时双精度浮点有效位数已经被吃掉了将近一半此时高阶层级系数不可信。对策按优先级排列第一始终用Zmat \ z而不是显式正规方程第二尽量保持数据点均匀覆盖圆盘避免数据全部集中在边缘或中心第三当 N 超过 36 且数据有缺口时改用带岭惩罚的拟合lambda 1e-6; % 正则化强度按残差量级试 coeff (Zmat*Zmat lambda*eye(N)) \ (Zmat*z(mask));正则化只在高阶系数明显跳动时启用λ 过大会把低阶系数也拉偏。判断方法在 λ0 和 λ1e−6 下分别拟合如果低阶系数前 6 项变化超过 1%说明 λ 太大。4.3 数据预处理从 CSV 导入、去离群点到样本均衡实测数据很少能直接喂给拟合函数。最常见的第一步是从 CSV 导入这也是 Matlab 日常数据处理里频率最高的操作之一T readmatrix(surface.csv); % 三列格式x, y, z x T(:,1); y T(:,2); z T(:,3);如果 CSV 是矩阵格式第 i 行第 j 列是网格高度需要自己生成网格坐标Zgrid readmatrix(heightmap.csv); [ny, nx] size(Zgrid); [xg, yg] meshgrid((0:nx-1) * pitch_x, (0:ny-1) * pitch_y); xg xg(:); yg yg(:); zg Zgrid(:);导入后先剔除无效值再处理离群点。离群点的典型特征是局部突变我用中值滤波做基准检测zmed medfilt1(z, 5); % 一维数据用 medfilt1 bad abs(z - zmed) 5 * std(z - zmed); z(bad) NaN; keep ~isnan(z) ~isnan(x) ~isnan(y); x x(keep); y y(keep); z z(keep);最后一步是样本均衡。如果数据在中心区域密集、边缘稀疏转台扫描特别常见最小二乘会被密集区主导等效于给中心数据加权。常见做法是把数据重采样到均匀网格或者对每个环带分别抽样后再合并拟合。重采样用 scatteredInterpolant 构造插值器插值到目标网格后同样要重新应用圆域掩膜。4.4 几个反复出现的误用对照最容易踩的坑都集中在坐标和定义上。把矩形像素坐标直接当物理坐标拟合不设 R、不掩膜圆外的四个角被当成面形数据低阶系数全面失真。用 atan(y/x) 代替 atan2二、三象限符号反转倾斜项和彗差项的正负完全错误。混用归一化与非归一化基函数系数直接对标设备输出的归一化系数量级对不上还查不出原因。还有一个常见错误是把工件坐标系的原点设置偏离光轴中心导致拟合结果出现虚假的彗差和倾斜——拟合前先做一次几何定心或者把原点设在孔径中心。5. 拟合结果落地iTrace 系数对齐、RMS 计算与波前重建5.1 与 iTrace 等像差仪导出的 Zernike 系数对齐iTrace 这类客观像差仪导出的 Zernike 系数有自己定义的排序、归一化符号。不同版本导出的 CSV 列序可能不一样直接拿过来对比之前必须先做一次“对表”。我一般用一段重排代码把自己的 (n, m) 项表转成 ANSI 单序号形式function c_std reorderToAnsi(coeff, terms) % 把系数从 (n, m) 项表转到 ANSI 单序号向量 c_std zeros(size(coeff)); for k 1:numel(coeff) n terms(k,1); m terms(k,2); j (n*(n2) m)/2 1; % ANSI 单序号公式1-based c_std(j) coeff(k); end end验证方法是构造一个已知离焦面形分别用你的函数和对方的导出格式处理对比 defocus 项的符号和量级。符号相反就翻转 m 的约定量级差 √2 或 2 倍就检查归一化方案。这一步花十分钟能省掉后面几周的“系数对不上”排查。5.2 从系数直接计算 RMS 与斯特列尔比归一化基函数有个直接红利总波前 RMS 就是去掉 piston 后的系数向量的二范数。这是因为各项正交且每项均方值为 1方差可以直接叠加c coeff(2:end); % 去掉 piston rms_wave sqrt(sum(c.^2)); % 波前 RMS单位与 z 相同 lambda 0.6328; % 例如 He-Ne 激光波长单位与 z 统一 strehl exp(-(2*pi*rms_wave/lambda)^2); % 小像差近似Strehl 近似式在 RMS 小于 λ/10 时误差可以接受超过这个范围只能走数值衍射计算。PV 值则需要从重构面形取不能直接从系数算——除非你有全部系数的解析极值表。5.3 波前重建成像的完整收尾拟合结束后把系数回代到细网格上成像是验证质量最直观的一步也可以顺手画出残差云图xi linspace(-R, R, 401); [xg, yg] meshgrid(xi); rho sqrt(xg.^2 yg.^2); theta atan2(yg, xg); mask rho R; Zr zeros(numel(xg), numel(coeff)); for k 1:numel(coeff) Zr(:, k) zernikeTerm(terms(k,1), terms(k,2), rho(:), theta(:)); end rec reshape(Zr * coeff, size(xg)); rec(~mask) NaN; imagesc(xi, xi, rec); axis image; colorbar;最后检查交叉验证误差把原始数据插值到同一网格后减去 rec得到残余面形残余的 PV 应该远小于主面形 PV。如果残余呈现明显的环形条纹说明 R 与真实孔径不贴合如果残余集中在边缘说明 N 取小了。跨设备、跨程序交换 Zernike 系数前永远先核对三件事——归一化基准、R 的取值、排序与符号约定这三个点对不上系数本身再精确也没有可比性。本文还有配套的精品资源点击获取 Zernike多项式拟合曲面：从理论到Matlab实现
简介Zernike多项式是光学波前分析与像差校正中的常用数学工具借助Matlab进行曲面拟合可大幅提升计算效率。这份代码面向计算机、电子信息工程、数学等专业的大学生及科研人员适用于课程设计、期末大作业、毕业设计或光学仿真任务。资源包仅含2个文件主要为1个.m函数文件和1个txt说明文件压缩包大小约3KB结构精简、定位明确。代码采用参数化编程参数可灵活调整思路清晰并附有详细注释同时提供可直接运行的案例数据兼容Matlab 2014至2024a版本。使用者无需从零搭建拟合流程只需修改参数或导入数据即可快速获得Zernike拟合结果并观察曲面重构效果有助于理解多项式拟合原理、缩短调试时间。目前已有204人学习下载适合刚接触Zernike拟合或希望在项目中快速落地该算法的读者。1. Zernike 多项式拟合曲面镜面检测与像差分析都在用的那套函数你手里有一台干涉仪或者共聚焦显微镜的测量结果几千到几十万个 (x, y, z) 散点描述一块镜面、一片晶圆或者一个光学元件的表面起伏。要把这批数据变成可对比、可传递、可反推的结果最通用的做法就是把面形投影到 Zernike 多项式上。这个标题里的 .rar 在光学圈流传了很多年解开后通常就是一个或一组 .m 文件输入离散坐标与面形高度输出一组 Zernike 系数。反直觉的地方在于拟合得到的几十个系数比原始云图更接近“答案”。离焦 0.8 µm、水平散光 0.2 µm——这些物理量直接从系数里读出来而不是从云图里量出来。原因是 Zernike 基函数在单位圆域上正交每个系数只对应一种像差模式互不干扰。正因如此从望远镜主镜检测到眼科波前像差仪比如 iTrace的数据交换都拿 Zernike 系数当标准语言。这篇文章按我平时给光学和量测方向同事梳理的思路展开先讲清 Zernike 多项式的定义、排序与归一化再给出一个可以直接抄进项目的 Matlab 拟合函数然后处理归一化半径、阶数和数据预处理这些最容易翻车的参数最后落到与像差仪数据对接和波前重建的验证技巧上。新手可以一步步跑通熟手可以直接跳到 3.2 和 4.2 对照自己的实现。2. Zernike 多项式拟合的理论基础正交性、排序与单位圆坐标2.1 定义式与两种主流排序ANSI 序和 Noll 序Zernike 多项式定义在单位圆盘上用极坐标 (ρ, θ) 表示ρ ∈ [0,1]θ ∈ [0,2π)。每一项由径向多项式与角向函数相乘Z_n^m(ρ, θ) N_n^m · R_n^{|m|}(ρ) · {cos(mθ), sin(mθ)}其中径向多项式 R_n^{|m|}(ρ) 的闭合形式是R_n^{|m|}(ρ) Σ_{s0}^{(n−|m|)/2} (−1)^s (n−s)! / [s! ((n|m|)/2 − s)! ((n−|m|)/2 − s)!] · ρ^{n−2s}约束条件n 与 |m| 同奇偶且 |m| ≤ n。这里的 n 是径向阶数|m| 是角向频率。归一化系数 N_n^m 按 ANSI Z80.28 标准取m0 时为 √(n1)m≠0 时为 √(2(n1))。这样每项在单位圆上的均方值为 1拟合得到的系数直接等于该项对波前方差的贡献单位与输入面形数据的单位一致。排序是个大坑。同一套 Zernike 多项式ANSI 序眼科行业常用、Noll 序天文望远镜常用和 Born Wolf 序干涉仪厂家常用的编号完全不同。我推荐的实现方式是显式维护一个 (n, m) 项表让系数含义自解释系数位置nm物理名称表达式100平移piston121−1y 方向倾斜2ρ sinθ311x 方向倾斜2ρ cosθ42−245° 散光√6 ρ² sin2θ520离焦defocus√3(2ρ²−1)6220° 散光√6 ρ² cos2θ73−1垂直彗差√8(3ρ³−2ρ) sinθ831水平彗差√8(3ρ³−2ρ) cosθ93−3三叶草30°√8 ρ³ sin3θ1033三叶草0°√8 ρ³ cos3θ提示表格里 m 为负时用 sin、为正时用 cos这只是惯用约定。iTrace、Zygo 等设备可能用相反符号跨设备对比系数前必须先验证方向定义。2.2 为什么必须在圆域上拟合正交性与“核函数”视角普通多项式比如幂级数 x^i y^j在矩形域上做拟合各项之间存在强相关性稍微加点噪声系数就乱跳。Zernike 多项式在连续单位圆盘上满足严格正交性∫∫ Z_i(x,y) · Z_j(x,y) dxdy π · δ_ij也就是说把任意光滑圆孔面形展开成 Zernike 级数时每一项的能量独立分离。这是它成为光学面形“标准语言”的根本原因。但要注意离散测量点不满足连续正交条件。当数据点在圆盘内近似均匀分布时设计矩阵 Z每列是一个基函数在所有数据点上的取值的各列近似正交当数据稀疏、偏心或者有大块缺失时列间会出现相关性。从核函数的角度看最小二乘解依赖 Gram 矩阵 K ZᵀZK 的对角占优程度决定了拟合稳定度。K 的条件数是 Z 的条件数的平方这就是为什么后面要用 QR 分解而不是直接解正规方程。另一个容易被忽略的点Zernike 多项式只在归一化单位圆内定义拟合前必须把物理坐标映射到单位圆并且把圆外数据彻底剔除。映射关系是 ρ √(x²y²)/Rθ atan2(y, x)其中 R 是归一化半径量纲与 x、y 一致。2.3 坐标变换与掩膜的最小实现这段代码是所有 Zernike 拟合函数的第一步先单独拿出来验证% x, y : 列向量物理坐标单位任意但必须一致 % R : 归一化半径通常取有效孔径半径 rho sqrt(x.^2 y.^2) / R; theta atan2(y, x); mask rho 1 1e-12; % 圆内有效点1e-12 是浮点容差逻辑说明ρ 超过 1 的点不在单位圆定义域内必须剔除这里的 mask 同时服务于两项任务——筛选有效数据、保证后续设计矩阵的行数正确。θ 必须用 atan2 而不是 atan(y/x)因为 atan 会把二、三象限的点误判到一、四象限导致倾斜项和彗差项的符号全部反掉。1e-12 的容差是为了把恰好落在边界上的点保留下来这种点通常来自网格采样时 ρ 恰好等于 1 的情况。3. 用 Matlab 写一个 Zernike 曲面拟合函数函数声明、设计矩阵与最小二乘3.1 函数声明输入输出与单位约定常见的 Matlab 拟合函数开头长这样文件名必须与函数名一致保存为zernikeFitSrf.mfunction [coeff, terms, zern, stats] zernikeFitSrf(x, y, z, R, N) % ZERNIKEFITSRF 用 Zernike 多项式拟合圆域面形数据 % 输入 % x, y, z : 面形采样点的坐标与高度同为列向量单位一致 % R : 归一化半径与 x, y 同单位 % N : 需要拟合的 Zernike 项数 % 输出 % coeff : N x 1 系数向量单位与 z 相同 % terms : N x 2 矩阵每行是该项的 (n, m) % zern : 拟合回代值只含有效圆域内的点 % stats : 结构体含残差 RMS、PV 和设计矩阵条件数这里要明确一个单位约定z 用 µm、nm 还是 mm由你的测量设备决定系数单位自动跟随 z。x、y 用 mmR 就用 mm三者保持同一套长度单位即可Zernike 基函数本身是无量纲的。关于 Matlab 函数声明还有一条规则一个 .m 文件里只有第一个函数对外可见其余都是子函数所以radialZernike、zernikeTerm这些辅助函数要么放同一个文件里要么各自建文件。3.2 径向多项式、单顶生成与设计矩阵构造先实现径向多项式。用阶乘公式直接求中等阶数n ≤ 30完全没问题function R radialZernike(n, m, rho) % 计算径向多项式 R_n^|m|(rho)rho 可以是向量 m abs(m); R zeros(size(rho)); for s 0:(n - m)/2 R R (-1)^s * factorial(n - s) ./ ... (factorial(s) * factorial((n m)/2 - s) * factorial((n - m)/2 - s)) ... .* rho.^(n - 2*s); end end逻辑说明求和上限 (n−m)/2 由 n 与 m 同奇偶保证为整数每一项的指数 n−2s 从 n 递减到 m因此 ρ0 处的值只在 m0 时非零这与 piston 项的物理意义一致。阶乘形式在 n 超过约 30 时可能引入较大舍入误差但一般面形拟合用到 7 阶36 项以内这个实现足够。接着是单顶生成与项表生成function Z zernikeTerm(n, m, rho, theta) % 生成第 (n, m) 项 Zernike 基函数按 ANSI 归一化 R radialZernike(n, m, rho); if m 0 Z sqrt(2*(n 1)) * R .* cos(m*theta); elseif m 0 Z sqrt(2*(n 1)) * R .* sin(abs(m)*theta); else Z sqrt(n 1) * R; end end function terms zernikeTerms(N) % 按 n 升序、(m -n:2:n) 生成前 N 项 (n, m) 表 terms zeros(N, 2); k 0; n 0; while k N for m -n:2:n k k 1; if k N break; end terms(k, :) [n, m]; end n n 1; end end系数 √(n1) 与 √(2(n1)) 就是 2.1 节说的 ANSI 归一化因子去掉它们也能拟合但系数含义会变成“未归一化基下的系数”无法直接和其他设备对比。项表生成用 while 循环按 n 逐层展开保证任意 N 都有明确的 (n, m) 对应关系。主函数把这些拼起来function [coeff, terms, zern, stats] zernikeFitSrf(x, y, z, R, N) assert(isequal(size(x), size(y), size(z)), x/y/z 必须同尺寸); x x(:); y y(:); z z(:); % 统一为列向量 rho sqrt(x.^2 y.^2) / R; theta atan2(y, x); mask rho 1 1e-12; % 有效圆域掩膜 if nnz(mask) N error(有效数据点 (%d) 少于拟合项数 (%d)请增大 R 或减小 N, nnz(mask), N); end terms zernikeTerms(N); Zmat zeros(nnz(mask), N); for k 1:N Zmat(:, k) zernikeTerm(terms(k,1), terms(k,2), ... rho(mask), theta(mask)); end coeff Zmat \ z(mask); % 直接最小二乘QR 分解 zern Zmat * coeff; % 回代重构 resid z(mask) - zern; stats struct(rmse, sqrt(mean(resid.^2)), ... pv, max(resid) - min(resid), ... cond, cond(Zmat)); end参数说明Zmat \ z(mask)用的是 Matlab 反斜杠运算符内部走 QR 分解加列主元比显式计算(Z*Z)\(Z*z)数值稳定得多因为正规方程会把条件数平方。stats.cond是诊断用的关键值后面 4.2 节会专门讲。如果有效点数小于项数最小二乘就没有唯一解这里直接报错而不是让结果静默出错。3.3 用已知系数的模拟面形验证函数没有写错拿到新写的拟合函数先别急着喂实测数据。构造一个已知系数组合的仿真面形跑一遍看能不能把系数原样拿回来Npix 401; [xg, yg] meshgrid(linspace(-1, 1, Npix)); rho sqrt(xg.^2 yg.^2); theta atan2(yg, xg); mask rho 1; % 构造0.5 平移 0.8 离焦 0.2 水平散光单位自定义 zs 0.5 0.8*sqrt(3)*(2*rho.^2 - 1) 0.2*sqrt(6)*rho.^2 .* cos(2*theta); zs(~mask) NaN; idx find(mask); [coeff, terms, zern, stats] zernikeFitSrf(xg(idx), yg(idx), zs(idx), 1, 15); disp([terms, coeff]); % 检查 (n, m) 与系数的对应 fprintf(RMSE %.3e, cond %.3e\n, stats.rmse, stats.cond);按 2.1 节的项表期望结果是 coeff(1) ≈ 0.5piston、coeff(5) ≈ 0.8离焦、coeff(6) ≈ 0.2水平散光。如果这三个位置对得上其他位置接近 10⁻¹⁰ 量级RMSE 接近机器精度说明径向多项式、归一化因子和拟合主流程都没问题。对不上时按三个方向查先查 R 与 mask 是否把数据圈对再查zernikeTerm里 m 的符号约定最后查zernikeTerms的排列是否与你的预期一致。4. Zernike 拟合参数怎么设半径、阶数与数据预处理4.1 归一化半径 R 选错是头号问题R 的定义是“物理坐标除以 R 后落在单位圆内”所以 R 必须等于有效孔径半径。常见的错误有两类。第一类R 取小了圆外本应参与拟合的点被 mask 砍掉真实面形的高阶成分会泄漏到低阶系数里最典型的表现是拟合残差沿着孔径边缘出现一条亮环。第二类R 取大了拟合区域只剩中心一小块系数描述的只是局部面形而且中心区域数据密度高会把边缘的真实形状完全稀释掉。更隐蔽的问题是不同的数据集用了不同的 R。改变了 R 等于改变了基函数本身所有系数都会变不存在一个简单的比例换算关系。所以同一批对比实验R 必须固定。对于干涉仪数据R 取镜头设计孔径半径对于轮廓仪扫描数据R 取有效数据点的最大半径或者略小一点保证圆内有足够多的点参与拟合。R_fit max(sqrt(x.^2 y.^2)) * 0.999; % 按数据边界自适应取 R这里乘 0.999 而不是直接用最大值是为了避免恰好落在边界上的散点因为浮点误差被 mask 误杀。确定 R 后把它作为常量传进函数并在输出结果里记录这是工程上防止系数误读的最低要求。4.2 阶数 N 与条件数能用 QR 就不要用正规方程拟合项数 N 决定最高阶 n。取太少会欠拟合面形里的真实像差被低阶项分摊取太多则高阶项开始拟合噪声。判断标准不是“看起来平滑”而是看残差有没有明显结构以及新增系数是否发散。条件数是更硬的指标。设计矩阵 Z 的条件数随 N 增长而用正规方程求解时实际参与运算的是 Gram 矩阵 ZᵀZ它的条件数是 cond(Z)²。以下是我在均匀网格和带缺口数据上看到的典型量级不同采样差异很大关键是自己会读 stats.cond项数 N最高阶 n均匀圆盘网格 cond(Z)稀疏/扇形缺口 cond(Z)153约 1–1010–10036710–10010³–10⁴661010²–10³典型大于 10⁵提示cond 超过 10⁶ 时双精度浮点有效位数已经被吃掉了将近一半此时高阶层级系数不可信。对策按优先级排列第一始终用Zmat \ z而不是显式正规方程第二尽量保持数据点均匀覆盖圆盘避免数据全部集中在边缘或中心第三当 N 超过 36 且数据有缺口时改用带岭惩罚的拟合lambda 1e-6; % 正则化强度按残差量级试 coeff (Zmat*Zmat lambda*eye(N)) \ (Zmat*z(mask));正则化只在高阶系数明显跳动时启用λ 过大会把低阶系数也拉偏。判断方法在 λ0 和 λ1e−6 下分别拟合如果低阶系数前 6 项变化超过 1%说明 λ 太大。4.3 数据预处理从 CSV 导入、去离群点到样本均衡实测数据很少能直接喂给拟合函数。最常见的第一步是从 CSV 导入这也是 Matlab 日常数据处理里频率最高的操作之一T readmatrix(surface.csv); % 三列格式x, y, z x T(:,1); y T(:,2); z T(:,3);如果 CSV 是矩阵格式第 i 行第 j 列是网格高度需要自己生成网格坐标Zgrid readmatrix(heightmap.csv); [ny, nx] size(Zgrid); [xg, yg] meshgrid((0:nx-1) * pitch_x, (0:ny-1) * pitch_y); xg xg(:); yg yg(:); zg Zgrid(:);导入后先剔除无效值再处理离群点。离群点的典型特征是局部突变我用中值滤波做基准检测zmed medfilt1(z, 5); % 一维数据用 medfilt1 bad abs(z - zmed) 5 * std(z - zmed); z(bad) NaN; keep ~isnan(z) ~isnan(x) ~isnan(y); x x(keep); y y(keep); z z(keep);最后一步是样本均衡。如果数据在中心区域密集、边缘稀疏转台扫描特别常见最小二乘会被密集区主导等效于给中心数据加权。常见做法是把数据重采样到均匀网格或者对每个环带分别抽样后再合并拟合。重采样用 scatteredInterpolant 构造插值器插值到目标网格后同样要重新应用圆域掩膜。4.4 几个反复出现的误用对照最容易踩的坑都集中在坐标和定义上。把矩形像素坐标直接当物理坐标拟合不设 R、不掩膜圆外的四个角被当成面形数据低阶系数全面失真。用 atan(y/x) 代替 atan2二、三象限符号反转倾斜项和彗差项的正负完全错误。混用归一化与非归一化基函数系数直接对标设备输出的归一化系数量级对不上还查不出原因。还有一个常见错误是把工件坐标系的原点设置偏离光轴中心导致拟合结果出现虚假的彗差和倾斜——拟合前先做一次几何定心或者把原点设在孔径中心。5. 拟合结果落地iTrace 系数对齐、RMS 计算与波前重建5.1 与 iTrace 等像差仪导出的 Zernike 系数对齐iTrace 这类客观像差仪导出的 Zernike 系数有自己定义的排序、归一化符号。不同版本导出的 CSV 列序可能不一样直接拿过来对比之前必须先做一次“对表”。我一般用一段重排代码把自己的 (n, m) 项表转成 ANSI 单序号形式function c_std reorderToAnsi(coeff, terms) % 把系数从 (n, m) 项表转到 ANSI 单序号向量 c_std zeros(size(coeff)); for k 1:numel(coeff) n terms(k,1); m terms(k,2); j (n*(n2) m)/2 1; % ANSI 单序号公式1-based c_std(j) coeff(k); end end验证方法是构造一个已知离焦面形分别用你的函数和对方的导出格式处理对比 defocus 项的符号和量级。符号相反就翻转 m 的约定量级差 √2 或 2 倍就检查归一化方案。这一步花十分钟能省掉后面几周的“系数对不上”排查。5.2 从系数直接计算 RMS 与斯特列尔比归一化基函数有个直接红利总波前 RMS 就是去掉 piston 后的系数向量的二范数。这是因为各项正交且每项均方值为 1方差可以直接叠加c coeff(2:end); % 去掉 piston rms_wave sqrt(sum(c.^2)); % 波前 RMS单位与 z 相同 lambda 0.6328; % 例如 He-Ne 激光波长单位与 z 统一 strehl exp(-(2*pi*rms_wave/lambda)^2); % 小像差近似Strehl 近似式在 RMS 小于 λ/10 时误差可以接受超过这个范围只能走数值衍射计算。PV 值则需要从重构面形取不能直接从系数算——除非你有全部系数的解析极值表。5.3 波前重建成像的完整收尾拟合结束后把系数回代到细网格上成像是验证质量最直观的一步也可以顺手画出残差云图xi linspace(-R, R, 401); [xg, yg] meshgrid(xi); rho sqrt(xg.^2 yg.^2); theta atan2(yg, xg); mask rho R; Zr zeros(numel(xg), numel(coeff)); for k 1:numel(coeff) Zr(:, k) zernikeTerm(terms(k,1), terms(k,2), rho(:), theta(:)); end rec reshape(Zr * coeff, size(xg)); rec(~mask) NaN; imagesc(xi, xi, rec); axis image; colorbar;最后检查交叉验证误差把原始数据插值到同一网格后减去 rec得到残余面形残余的 PV 应该远小于主面形 PV。如果残余呈现明显的环形条纹说明 R 与真实孔径不贴合如果残余集中在边缘说明 N 取小了。跨设备、跨程序交换 Zernike 系数前永远先核对三件事——归一化基准、R 的取值、排序与符号约定这三个点对不上系数本身再精确也没有可比性。本文还有配套的精品资源点击获取