ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

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

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 系数到 PSFzernike_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 系排列具体实现里可能微调但离焦、彗差、球差这几项的位置基本一致单索引名称对应像差典型来源4defocus离焦焦面位置误差5 / 6astigmatism像散镜片应力、柱面误差7 / 8coma彗差偏心、视场离轴9trefoil三叶草镜面加工低阶误差11spherical球差球面镜、平行平板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 radianPSF 几乎看不出变化这是新手最常踩的坑。第二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_psfwhich 返回实际路径说明加载成功返回 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由系统给定入瞳直径或半径直径/半径搞混是头号错误先看函数注释wavelength0.55 / 0.6328 μm工作波长系数单位是 waves 时不影响相位本身focal_length由系统给定焦距决定 PSF 空间缩放影响 PSF 尺寸不改变形状grid_size256 / 512FFT 网格边长取 2 的幂像差大时提到 1024normalizetrue是否做能量归一化算 Strehl 比时保持 truecoefficient_unitwaves系数单位传 micrometer 时内部会除以波长这几个参数里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 调用后顺手就能做判断时只看相对误差不用绝对阈值因为网格大小和归一化方式都会同时影响两个量的数值。本文还有配套的精品资源点击获取
RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进