Zernike系数到PSF：光学仿真中zernike_psf原理与MATLAB实现

简介：这是一份面向光学工程与视觉科学研究者的波前光学与Zernike像差分析工具包，围绕点扩散函数（PSF）计算与成像质量评估展开，可帮助理解Zernike系数如何影响系统成像分辨率，在光学设计、视觉模型验证与成像系统优化等场景中具有实用价值。内容基于WavefrontOptics项目，整合了56个MATLAB函数脚本、11个mat数据文件、9个txt参数文档、6个pdf参考论文等共86个文件，压缩包整体约11.36MB，体量小巧且便于快速部署。目前已有247人学习下载。通过Zernike多项式可对球差、彗差、像散等常见像差进行分解和量化，并借助提供的wvfComputePSF、wvfPlot等核心工具模拟人眼光学系统，分析Stiles-Crawford效应、色差与像差对PSF的影响。包内附带Thibos等人眼模型文献与实验数据，便于复现研究结果和验证算法。目录结构清晰，包含tutorial、docs、validate、scripts、utility等模块，适合具备一定光学基础的研究者用于教学演示、算法验证或二次开发。

1. 从 Zernike 系数到 PSF：zernike_psf 这条链在做什么

拿到一组 Zernike 系数，第一件事往往不是看波前图，而是想知道像差叠加之后点光源成像到底糊成什么样。zernike_psf 就是 WavefrontOptics 这类工具包里干这件事的函数：输入系数和入瞳参数，输出点扩散函数 PSF，顺带能算 Strehl 比和 MTF。做过光学装调的人都有体会，干涉仪给的是系数，评审、写报告、判断成像质量，要的却是 PSF。

这套链路在光学设计、系统装调、显微成像和自适应光学里是标配。下面把从系数到 PSF 的完整链路拆开讲，给出可复现的 MATLAB 调用、参数设置依据和验证手段，新手能照做，老手能对照检查自己的常见疏漏。

2. zernike_psf 的计算原理：光瞳相位到 PSF 的四步变换

2.1 为什么用 Zernike 系数而不是直接给出相位分布

Zernike 多项式在单位圆上正交，每一项对应一种典型像差形态，这是它成为波前描述标准的核心原因。第 1 项是平移（piston），第 2、3 项是倾斜，第 4 项是离焦，第 5、6 项是像散，第 7、8 项是彗差，第 11 项是球差。用系数描述波前，每一项可以独立调整、独立评估；干涉仪和夏克-哈特曼波前传感器输出的也正是这套系数。要把系数还原成波前相位，做法是在光瞳每个采样点上求多项式值再线性叠加：

φ(ρ,θ) = Σ cᵢ · Zᵢ(ρ,θ)

其中 ρ 是归一化半径（0 到 1），θ 是极角，cᵢ 是第 i 项系数。两个容易忽略的点：一是系数单位通常约定为波长（waves），干涉仪如果输出微米要先除以工作波长；二是归一化半径以入瞳半径为基准，同一套光学系统换口径时系数不变，变的只是相位分布覆盖的空间范围。常见低阶项的单索引编号如下表，按 OSA/ANSI 系排列，具体实现里可能微调，但离焦、彗差、球差这几项的位置基本一致：

2.2 光瞳函数构造与 FFT 的关键代码

zernike_psf 内部的计算可以拆成四步：生成单位圆网格、按系数叠加相位、构造复光瞳函数、做傅里叶变换取模平方。用代码表达就是：

[x, y] = meshgrid(linspace(-R, R, N), linspace(-R, R, N)); [theta, rho] = cart2pol(x, y); rho = rho / R; % 归一化到单位圆，R 是入瞳半径 pupil = double(rho <= 1); % 圆形孔径振幅掩膜 phase = zeros(size(rho)); for k = 1:numel(c) phase = phase + c(k) .* zernike_term(k, rho, theta); end pupil = pupil .* exp(1i * 2 * pi * phase); psf = abs(fftshift(fft2(ifftshift(pupil)))).^2;

这段代码里有三个决定结果正确性的细节。第一，系数单位是波长时，相位里必须乘 2π，漏掉这一项会让离焦 0.1λ 变成 0.1 radian，PSF 几乎看不出变化，这是新手最常踩的坑。第二，fft2 之前用 ifftshift、取模之后用 fftshift，不能来回都用同一个函数，否则在奇数尺寸网格下 PSF 中心会偏一个像素。第三，振幅掩膜不只包含圆形判定，实际系统里的中心遮挡、支撑桁架阴影都要以乘法形式作用在这里，只能在改掩膜这一步介入，不应该去动系数。zernike_term 可以用工具包自带的多项式生成函数，也可以按 Noll 或 ANSI 编号的递推公式自己实现，两种方式的编号顺序必须和系数向量严格一致。

2.3 傅里叶变换这一步的物理含义

从光瞳面到焦平面的传播，在傍轴近似下就是夫琅禾费衍射，数学形式恰好是一次二维傅里叶变换，所以 PSF 不需要逐点做衍射积分，一次 FFT 就能完成，这也是 zernike_psf 这类函数快的根本原因。反过来，对 PSF 再做一次傅里叶变换得到的是 OTF，取模就是 MTF——同一族函数里顺手就能算出来的量。理解这条互逆关系，排错时才有方向：PSF 里出现规则栅格条纹，问题多半在采样密度或补零方式，而不是系数本身；而低频响应缺失，则要回头看振幅掩膜是不是把不该遮的地方遮掉了。

3. 最小可复现：MATLAB 里跑通 zernike_psf 的完整调用

3.1 解压目录与 addpath 路径检查

WavefrontOptics-master.zip 解压后，典型结构是源码目录、示例脚本和文档说明。写第一行调用之前，先确认 zernike_psf.m 确实存在，再把根目录递归加入 MATLAB 路径：

addpath(genpath('D:\work\WavefrontOptics-master')); which zernike_psf

which 返回实际路径说明加载成功；返回 empty 时，优先怀疑 genpath 没覆盖子目录。这类工具包通常把函数放在 src 或 private 子目录里，只 addpath 根目录会找不到。另外注意 MATLAB 对函数名按路径顺序解析，如果机器上装了多个工具箱存在同名函数，which 会给出第一个命中的路径，调用前确认它指向的是 WavefrontOptics 里的文件，避免同名函数互相遮蔽。

3.2 完整示例代码与逐段说明

% 主脚本: demo_zernike_psf.m clearvars; close all; % 1. 定义 Zernike 系数，单位: 波长 waves，编号按 OSA 单索引 c = zeros(12, 1); c(4) = -0.12; % 离焦 -0.12λ c(5) = 0.06; % 0° 像散 +0.06λ c(8) = 0.08; % 彗差 +0.08λ c(11) = -0.05; % 球差 -0.05λ % 2. 光学系统参数 D = 8e-3; % 入瞳直径 8 mm lambda = 632.8e-9; % He-Ne 波长 f = 50e-3; % 焦距 50 mm N = 512; % 网格点数，取 2 的幂 % 3. 调用 zernike_psf [psf, xgrid] = zernike_psf(c, ... 'aperture', D, 'wavelength', lambda, ... 'focal_length', f, 'grid_size', N); % 4. 显示，横轴换算成微米 imagesc(xgrid*1e6, xgrid*1e6, psf); axis image; colormap(hot); colorbar; xlabel('x (um)'); ylabel('y (um)'); title('PSF: defocus -0.12\lambda, spherical -0.05\lambda');

代码逻辑说明：系数向量 c 的长度决定最多计算到第几项，没有赋值的位置自动按零处理。aperture 传的是直径而不是半径，很多人在这把 8 mm 写成 4 mm，结果 PSF 尺度差一倍，后面和理论值对照时完全对不上。zernike_psf 返回的 psf 已做归一化，xgrid 是像面上每个像素对应的物理坐标，显示时乘 1e6 转成微米更直观。如果调用时报Unrecognized parameter name，打开 zernike_psf.m 看函数头部的注释，以工具包实际定义的参数名为准，不同版本之间参数拼写可能略有差异。

提示：如果算出来的 PSF 和预期差一个数量级，先查 aperture 传的是直径还是半径，这是这类函数最常见的第一个坑。

3.3 zernike_psf 参数速查表

这几个参数里，aperture、focal_length、wavelength 描述的是光学系统本身，grid_size 是数值计算参数，normalize 和 coefficient_unit 是约定参数。调参顺序一般是先固定系统参数，再根据第 4 章的采样规则决定 grid_size，最后核对单位约定。

4. zernike_psf 参数设置：采样密度、孔径遮挡与编号约定

4.1 网格采样密度与混叠的关系

FFT 算 PSF 的隐含条件是光瞳被离散采样，采样密度直接决定像面视场大小。像面像素间隔满足 Δu = λf / (N·Δp)，其中 Δp 是光瞳面采样间隔，N 是网格边长。Δu 太大，艾里斑只占两三个像素，看不出旁瓣结构；Δu 太小，计算区域外的高频成分折返进来，PSF 周围会出现假的周期条纹。

经验法则是让光瞳半径对应 64 到 128 个像素。按第 3 章的参数，入瞳半径 4 mm、N = 512 时，单像素 Δp = 31.25 μm，代入公式得到 Δu ≈ 2.07 μm。艾里斑第一暗环半径约 1.22λf/D = 4.83 μm，对应大约 2.3 个像素，属于能看清细节又不浪费计算量的中间档；如果要精细比较旁瓣或 Strehl 比的小数点后两位，把 N 提到 1024，像素间隔减半。判断有没有混叠，最简单的方法是固定其他参数只改 N，如果 PSF 中央区域的形状随 N 明显变化，说明原来的采样不够。

4.2 中心遮挡与振幅掩膜修改

反射式望远镜、折返镜头都有中心遮挡，直接拿圆形孔径算会高估低频响应，PSF 第一暗环会偏深。zernike_psf 如果支持遮挡参数（类似 obscuration_ratio），直接传入遮挡直径与入瞳直径的比值即可；如果不支持，在拿到 pupil 后逐点修改振幅掩膜：

pupil(rho <= 0.3) = 0; % 30% 线性直径中心遮挡 pupil = pupil .* double(rho <= 1); % 重新确认边界 psf = abs(fftshift(fft2(ifftshift(pupil)))).^2;

注意遮挡不是简单地把能量减去一部分，它改变了光瞳函数的空间形状，能量会从中心主瓣向外围旁瓣转移，Strehl 比下降幅度比能量损失比例更大。改完掩膜之后，建议先用 imagesc 画一遍 pupil 的实部，确认遮挡圆环的位置和比例正确，再继续后续计算。支撑桁架（spider）同理，在掩膜上叠一条细线状黑色区域，PSF 会出现典型的十字衍射条纹。

4.3 系数单位、符号与编号顺序

这是 zernike_psf 类函数最常翻车的地方，三个约定必须同时对齐：单位是 waves 还是微米、符号是凸起为正还是凹陷为正、编号是 Noll 还是 OSA/ANSI。同一组干涉仪数据，单位用微米、编号用 Noll 写进去，PSF 会完全对不上，而且这种错不会报任何错误信息。

我一般这样处理：收到数据先问清楚干涉仪软件的导出设置，再在调用前加一个显式转换，把约定问题留在代码里可见的位置：

% 干涉仪给的是微米，函数内部按 waves 处理 if strcmp(unit, 'um') c_waves = c_um / (lambda * 1e6); end

转换完不要急着算 PSF，先把系数代回波前表达式画一幅波前图，和干涉仪软件里的波前图对比。条纹方向一致、幅度对得上再往下走。这一步虽然多花两分钟，但能把后面一半的排错时间省掉——PSF 算错时，追溯到底往往是编号或符号错位，而不是傅里叶变换本身。

5. 把 zernike_psf 封装成 function 节点：批量扫描与调用校验

5.1 为什么要包一层批量接口

单组系数跑一次没什么问题，但实际工作里很少只算一次。要扫描离焦量从 -0.5λ 到 +0.5λ 的 PSF 序列做景深判断，要对比不同像差组合下的 Strehl 比，还要把结果喂给优化循环。每次手写参数容易出错，常见做法是封装一个批量 function 节点：输入系数矩阵，输出 PSF 数组和关键指标，上层逻辑只关心输入输出，不关心 FFT 和归一化细节。

5.2 批量封装函数代码

function results = zernike_psf_batch(coeff_matrix, D, lambda, f, N) % coeff_matrix: n×m，每行一组 Zernike 系数，m 为项数 % 返回结构体数组，含 psf、strehl、rms_wfe 三个字段 n_cases = size(coeff_matrix, 1); results = repmat(struct('psf',[],'strehl',[],'rms_wfe',[]), n_cases, 1); % 先算无像差参考 PSF，用于 Strehl 比归一化 [ref_psf, ~] = zernike_psf(zeros(size(coeff_matrix(1,:)))', ... 'aperture', D, 'wavelength', lambda, ... 'focal_length', f, 'grid_size', N); ideal_peak = max(ref_psf(:)); for k = 1:n_cases [psf, ~] = zernike_psf(coeff_matrix(k,:)', ... 'aperture', D, 'wavelength', lambda, ... 'focal_length', f, 'grid_size', N); results(k).psf = psf; results(k).strehl = max(psf(:)) / ideal_peak; % 去掉 piston 后按系数平方和近似 RMS，批量筛选够用 results(k).rms_wfe = sqrt(sum(coeff_matrix(k,2:end).^2)); end end

这个封装有三点值得说明。一是输入行向量转成列向量，避免 zernike_psf 内部对系数维度敏感；二是 ideal_peak 单独用无像差 PSF 算出来，而不是写死常量，因为 grid_size 一改峰值就变；三是 RMS 用去掉 piston 后的系数平方和只是近似，精确计算要按每项 Zernike 的方差权重来，但做批量筛选已经足够。

5.3 暴露给外部时的入参校验与 schema 错误

如果这批模块要暴露给上层脚本或 Agent 调用，也就是把它挂成一个 function 节点，输入校验就变得很关键。函数本体是 MATLAB，调用方可能是 Python 或其他语言，常见的传参错误有两类：参数名拼写不一致，比如把 aperture 写成 aperture_size；数组维度不对，比如系数传成了行向量。前者可以用 inputParser 做严格参数名校验，后者在封装入口加一行assert(size(coeff_matrix, 2) >= 1, 'coeff_matrix must be n×m')即可。

实际对接时经常看到的报错是api error: 400 invalid schema for function 'artifact' ... is not a "regex"，这类问题出在调用框架那边的 schema 定义：比如给某个字符串参数写了^(?!.*$)[^\p{Cc}...这种正则，部分环境不支持\p{Cc}写法，schema 校验直接拒绝。解决方式是通读 schema 里每个字段的 pattern，把不支持的转义序列换成普通字符类。回到 zernike_psf 的封装，教训是一样的：对外暴露的每个参数都要写明类型、单位和取值范围，由框架侧做 schema 校验，不要让一个非法字符串运行到第 3 章那张参数表里才报错。

6. zernike_psf 结果验证：能量守恒、衍射极限与系数回代

6.1 能量守恒校核

PSF 是光瞳函数傅里叶变换的模平方，离散 FFT 下满足 Parseval 恒等式：PSF 所有像素求和后乘像素面积，应该等于光瞳面振幅平方和除以 N²。每次跑完新参数我都会做一次：

dx = xgrid(2) - xgrid(1); energy_psf = sum(psf(:)) * dx^2; energy_pupil = sum(abs(pupil(:)).^2) / N^2; fprintf('ratio = %.4f\n', energy_psf / energy_pupil);

如果比值明显偏离 1，先查 fftshift/ifftshift 是否配对，再查 pupil 里有没有 NaN，以及网格边界处 rho 恰好等于 1 造成的异常点。

6.2 与理论艾里斑对照

无像差情况下 PSF 应与艾里斑解析解一致，第一暗环半径 1.22λf/D。用第 3 章参数计算：D = 8 mm、f = 50 mm、λ = 632.8 nm，得到 r = 4.83 μm。取 PSF 过中心一行做截面，找第一个极小值位置，如果偏离理论值超过一个像素，说明采样参数需要复核。

6.3 系数回代与 Strehl 自检

最后一个技巧是系数回代：小像差下 Strehl ≈ exp(-(2πσ)²)，σ 是 RMS 波前误差（waves）。把输入系数算出近似 RMS，再和 PSF 峰值相对无像差峰值的比值互相对照，两者差超过 20% 时，多半是 Zernike 编号顺序或符号里有一个错了。这种自检不需要额外工具，任何一次 zernike_psf 调用后顺手就能做，判断时只看相对误差，不用绝对阈值，因为网格大小和归一化方式都会同时影响两个量的数值。

本文还有配套的精品资源，点击获取