简介在光学成像系统分析中点扩散函数、调制传递函数与泽尼克多项式是评价像差与成像质量的核心工具。面向光学设计、图像处理和波前模拟学习者该资源提供了一套从泽尼克多项式拟合到点扩散函数与调制传递函数计算的可运行脚本并配有讲解用幻灯片和网页讲义能够帮助使用者对照公式理解代码逻辑、快速掌握光学仿真流程。压缩包共18个文件以9个M脚本为主涵盖波前像差计算、点扩散函数生成、调制传递函数曲线绘制等核心模块另有若干示意图、幻灯片及说明性文档辅助学习整体仅1.09MB轻量但完整适合直接运行调试。已有2210人学习下载对于希望在光学仿真中减少从零摸索成本、系统理解PSF与MTF理论的研究者和工程师来说是一份实用的上手资料。1. 计算点扩散函数是在算什么从PSF到MTF再到波前像差计算点扩散函数PSF这件事表面上是一个傅里叶变换实际上是把一个光学系统的全部脾气摸清楚。一个理想光学系统点光源成像后应该还是一个点但衍射和像差会让这个点摊开成一团光斑——这团光斑的强度分布就是 PSF对它做傅里叶变换取模得到的就是 MTF。而 Zernike 多项式在这里扮演的角色是用一组正交基函数把波前像差拆成看得懂的成分离焦、球差、彗差、像散每一项都有明确的物理含义和可调系数。做机器视觉、工业镜头选型、光刻照明系统设计或者自适应光学的人都绕不开这套计算流程。市面上很多商业软件帮你把 PSF 和 MTF 算好了但如果你手里只有一堆实测的 Zernike 系数或者需要把自定义孔径、非圆对称像差加进系统里评估自己动手算一遍就是唯一的可靠路径。这篇文章从理论公式一路推到可运行的 Python 代码和参数调试聚焦在“怎么从 Zernike 系数算出 PSF 和 MTF以及算了以后怎么确信算对了”。2. Zernike 多项式的物理意义与波前表示从光圈坐标到像差分量用 Zernike 多项式来描述波前像差最早是 Frits Zernike 在 1934 年提出的核心理念是在单位圆内定义一组互相正交的二维多项式任何连续波前都可以表示成这些多项式的线性组合。这个思路跟傅里叶级数类似区别在于 Zernike 定义在圆域上而光学系统的光瞳天然就是圆的。波前 W(x, y) 可以表示为W(x, y) Σ c_i · Z_i(x, y)其中 Z_i 是第 i 项 Zernike 多项式c_i 是该项的系数单位通常取“波长”waves。每一项 Z_i 由角向频率 m 和径向阶数 n 决定展开来写就是径向函数和角向函数的乘积Z_i(r, θ) R_n^m(r) · sin(mθ) 或 cos(mθ)r 是归一化径向坐标光瞳边缘为 1θ 是极角。径向函数 R_n^m(r) 的形式是固定的多项式组合网上随处可以查到闭式表达式这里不再贴冗长的公式。真正容易混淆的是排序方式Zernike 多项式常用的排序有 Noll 序号、Fringe 序号、OSA 标准序号三种同一项在不同排序里的编号完全不同混用是新手最常见的翻车点。比如离焦项Fringe 序号是 4对应 z[4]Noll 序号是 4一致但彗差两项就不一样了——Fringe 里是 7 和 8Noll 里是 7 和 8也有差别像散项差得更多。如果从实验设备直接导出的 Zernike 系数务必先确认设备用的是哪一种排序。2.1 前几阶 Zernike 项与像差的对应关系实际光学设计里低阶像差解释了大半的成像问题高阶项通常只在强离轴系统或自由曲面系统里才需要关注。下表列出 Fringe 排序下前 15 项对应的像差名称和典型物理来源写计算程序之前先把这张表放在手边序号nm像差名称物理来源典型符号约定100活塞常数相位不影响 PSF可忽略211倾斜 X光斑位置偏移符号决定偏移方向31-1倾斜 Y光斑位置偏移符号决定偏移方向420离焦对焦不准正负对应焦点前后52-2像散 45°柱面/倾斜面与第6项正交622像散 0°柱面元件/装配应力参考轴选择731彗差 X偏心或光轴不重合与倾斜方向相关83-1彗差 Y偏心或光轴不重合同上933三叶草镜片制造误差三项对称103-3三叶草镜片支撑应力与第9项方向正交1140初级球差球面本身固有像差正为边缘聚焦更近1242高阶像散高次曲面残差少见134-2高阶像散同上少见1444四叶草车削残留高次面154-4四叶草装配应力高次面工程上有个非常常见的做法把 Zernike 系数理解为“对应的 Seidel 像差量”即某一项系数值 0.25 就意味着该像差的波前 RMS 贡献约 0.25 波长。不过严格说只有把波前展开后按项做 RMS 统计才算准确因为多项式虽正交但同阶项叠加后的 RMS 与分项系数之间存在勾股关系RMS_total sqrt(Σ c_i²)这个性质在很多相位反演算法里被直接用来约束解空间。做仿真时如果只想调一种像差直接把对应项的系数设成非零值即可其他项保持 0。2.2 为什么用 Zernike 而不是直接搞一个波前函数实践里有人会把波前直接写成一个简单曲面比如球面波前用 W A·ρ²然后去算 PSF这样算出来的结果并不严谨。Polynomial 系数一旦换成物理坐标就丢掉了两个关键性质——正交性和旋转对称性。Zernike 的每一项在圆域内正交意味着调整某一项系数不会改变其他项对波前的贡献这在优化过程里极其重要你调彗差系数不会鬼使神差地影响离焦量的评估。另外Zernike 多项式的旋转对称性角向频率 m对应光学系统的对称特征。m 0 的项旋转对称m 1 的项是“蝴蝶形”m 2 的像散在旋转 90° 后取反号。利用这个性质可以快速判断仿真结果是否符合物理直觉——比如旋转对称的球差项算出来的 PSF 必须是旋转对称的如果不对称那一定是采样网格或坐标定义出了错。这个特性后面章节里做验证时会反复用到。3. 用 Python 从 Zernike 系数算出 PSF 和 MTF 的最小可运行实现光学系统的 PSF 计算基于标量衍射理论核心公式是光瞳函数的傅里叶变换模平方PSF(x_f) |FFT{P(x_p)}|²其中 P(x_p) 是光瞳平面上的复振幅分布由振幅透过率 A(x_p) 和相位项组合而成P(x_p) A(x_p) · exp(j · (2π/λ) · W(x_p))A 代表孔径形状圆孔内部为 1外部为 0W 是波前像差单位波长直接取 Zernike 展开结果。MTF 反过来是 PSF 的傅里叶变换模归一化到零频后取模MTF(f) |FFT{PSF(x_f)}| / |FFT{PSF(x_f)}|_(f0)需要注意一个关键点这里 PSF 的单位网格要和频率空间网格匹配即最后输出的 MTF 横轴是用像素表示的采样频率要换算成物理频率lp/mm必须知道实际系统里的像素缩放关系。3.1 网格生成与 Zernike 多项式函数实现网格生成这一步直接影响计算结果的精度。工程上最常用的是“奇偶数采样”方案采样点数 N 取偶数坐标范围从 -1 到 1 取 N 个点这样避免了原点落在网格中心时 FFT 常见的偏移问题。Zernike 多项式函数采样的传统写法是定义于单位圆内——也就是把物理光瞳半径缩放到 1采样完以后把 r 1 的位置掩膜掉。下面这段实现把多项式生成从 ANSI C 常见写法改成 numpy 向量化实现效率高且不容易错import numpy as np def zernike_poly(n, m, rho, theta): 计算单个 Zernike 多项式在极坐标网格上的值 n: 径向阶数, m: 角向频率(带符号) rho: 归一化径向坐标 (0~1), theta: 方位角 返回: 与 rho 同形状的浮点数组 # R_n^m 径向多项式用递推公式 # 先处理 m0 的特殊情况 if m 0: # p 为多项式的半阶 s_max n // 2 else: s_max (n - abs(m)) // 2 radial np.zeros_like(rho) for s in range(s_max 1): # 组合数直接用阶乘计算避免引入 scipy 依赖 coef ((-1) ** s) * np.math.factorial(n - s) / ( np.math.factorial(s) * np.math.factorial((n abs(m)) // 2 - s) * np.math.factorial((n - abs(m)) // 2 - s) ) radial coef * (rho ** (n - 2 * s)) # 角向分量: 根据 m 的符号选择 cos 或 sin if m 0: angular np.cos(m * theta) elif m 0: angular np.sin(-m * theta) else: angular np.ones_like(theta) norm np.sqrt(2 * (n 1)) if m ! 0 else np.sqrt(n 1) return radial * angular * norm代码逻辑说明s_max决定多项式迭代的项数radial部分由 R 多项式的标准求和式算出angular按 m 符号选择三角基函数——这个符号约定对应“径向多项式为正、角向用 sin/cos 区分方向”的常用形式与 Zemax 导出的数据符号规则一致。norm系数使得多项式在单位圆上 RMS 等于 1这样 Zernike 系数就可以直接代表 RMS 值实测数据里导出的系数如果不做 RMS 归一化这一步会额外引入一个sqrt(2)或sqrt(n1)的差异。调用时只需构造极坐标网格N 512 # 采样点数偶数 x np.linspace(-1, 1, N, endpointFalse) X, Y np.meshgrid(x, x) rho np.sqrt(X**2 Y**2) theta np.arctan2(Y, X)这里endpointFalse很关键配合偶数 N网格点关于原点对称FFT 以后的频谱中心不会偏移半个像素。3.2 从波前相位到 PSF 再到 MTF 的完整计算函数把相位、孔径、傅里叶变换串起来核心计算函数可以压缩成下面这样def compute_psf_mtf(zernike_coeffs, wavelength0.55e-3, N512, pupil_radius1.0, pixel_scale1.0): 输入: zernike_coeffs: 列表, 按单位圆内 RMS 归一化的 Zernike 系数 长度自适应, 缺失项按 0 处理 wavelength: 波长, 单位 mm (可见光约 0.00055) N: 采样网格边长 (偶数) pupil_radius: 光瞳半径, 决定孔径边缘落在哪个像素 pixel_scale: 输出 PSF 的像素缩放, 方便后续换算 MTF 返回: psf: N×N 强度分布, 已归一化到总能量 1 mtf: N×N 调制传递函数, 零频归一化为 1 W: 波前图 (单位: 波长), 用于调试可视化 # 1. 极坐标网格 x np.linspace(-1, 1, N, endpointFalse) X, Y np.meshgrid(x, x) rho np.sqrt(X**2 Y**2) theta np.arctan2(Y, X) # 2. 振幅孔径: 半径为 pupil_radius 的圆, 外部置零 aperture (rho pupil_radius).astype(float) # 3. 累加 Zernike 波前 (单位: 波长) W np.zeros_like(rho) for idx, coef in enumerate(zernike_coeffs): if coef 0: continue n int(np.sqrt(idx 1)) - 1 # 这个索引是近似, 见下文说明 # 实际使用时建议传入 (n, m) 对, 这里简化 m 0 # 占位, 需要按排序映射实际 m Z zernike_poly(n, m, rho, theta) W coef * Z # 4. 构造复振幅光瞳函数 phase (2 * np.pi / wavelength) * W pupil aperture * np.exp(1j * phase) # 5. FFT 算 PSF: 光瞳函数投影到焦平面 psf np.abs(np.fft.fftshift(np.fft.fft2(pupil))) ** 2 # 归一化总能量为 1 psf / psf.sum() # 6. MTF: PSF 再做一次 FFT, 模归一化 mtf np.abs(np.fft.fftshift(np.fft.fft2(np.fft.ifftshift(psf)))) mtf / mtf[0, 0] # 零频归一化 return psf, mtf, W代码里第 3 步有一个关键的简化从数组索引倒推 (n, m) 并不是一个可靠的做法依赖特定排序表。工程上最稳的方式是在函数入口处直接传入(idx, n, m, coef)四元组列表或者用一个标准 Zernike 排序表映射。上面代码只是展示整体数据流结构真用的时候建议把(n, m, coef)作为显式输入。FFT 两步分别用了fftshift和ifftshift的组合这是为了避免光瞳网格的偏移误差第一步fft2之后用fftshift把零频挪到中心第二步由于psf已经是中心化过的数据做 FFT 之前要先ifftshift还原到 FFT 算法的原点布局。这两个函数混用错了一个MTF 会整体错位而且很难从视觉上察觉。3.3 离焦与球差的仿真验证用一个 20 行以内的驱动脚本把上面函数跑起来# 只有 0.25 波长的初级球差 (n4, m0, 系数 0.25) z_coeffs [(4, 0, 0.25)] W np.zeros((N, N)) for n, m, c in z_coeffs: Z zernike_poly(n, m, rho, theta) W c * Z psf, mtf, W compute_psf_mtf_from_phase(W, wavelength0.55e-3, N512)如果只加了球差PSF 结果应该是中心亮斑周围带对称的同心环结构沿着光轴两侧离焦环结构会出现非对称的明暗变化这是球差“焦点位移导致模糊不对称”的典型特征。MTF 在某个中间频率处可能出现凹陷甚至落到零再回升这就是经典的“焦点外 MTF 带有频带缺口”现象做机器视觉镜头评测时如果看到这种 MTF 形状可以直接判断系统有残留球差。4. 像差参数怎么设系数单位、采样率与 MTF 曲线验证跑通第一版代码之后下一步是把仿真结果调到跟实验或设计数据对得上。这个阶段最磨人的不是代码逻辑而是参数设定的一致性。下面三个参数是最常出问题的。4.1 Zernike 系数单位与符号约定Zernike 系数的物理单位有两种习惯一种是以波长为单位的波前值RMS另一种是直接给出“光程差”OPD的单位是微米或纳米。两者换算关系是OPD coeff × wavelength。很多从商用软件导出的系数写的是“waves RMS”拿到手直接用就行如果是干涉仪采出的原始数据单位很可能是微米级 OPD需要除波长。符号约定方面各软件用的坐标系定义不同Zemax 和 Code V 对同一项 Zernike 的符号就可能相反。节省时间的做法是先用单一大像差如 1 波长离焦跑一次仿真确认 PSF 的模糊方向和像差符号的对应关系再拿实验数据对照一次后面大批量数据就按这个映射关系校准。4.2 采样网格点数与孔径边缘像素对精度的影响网格 N 的选择是一个典型的精度-速度权衡。FFT 计算的时间复杂度是 O(N² log N²)N 翻倍速度掉四倍。工程上用 N512 起步绝大多数单视场点扩散函数计算都在一秒钟之内完成。真正影响精度的是光瞳边缘落在哪几个像素上孔径圆边界的锯齿状量化误差会造成 PSF 高频部分的虚假能量表现是 MTF 在高频区域出现“裙边”上翘。缓解办法是给孔径边缘加过渡带使用超采样反走样。一段常用的小技巧# 孔径边缘用 2 像素平滑过渡, 减少振铃 edge 2.0 # 过渡带宽度 (像素) aperture 0.5 * (1 - np.tanh((np.abs(rho) - pupil_radius) / edge))把硬边缘改成软边缘之后PSF 外围的伪振荡幅度明显下降代价是中心强度略微下降、斯特列尔比Strehl Ratio计算值会比真实值低约 0.5% 以下可以接受。高频 MTF 的精度收益通常远大于这个损失。4.3 离焦扫描验证MTF 曲线随离焦量的变化规律一个典型的工程验证场景是“离焦扫描”给 Zernike 系数中离焦项n2, m0设置从 -1 到 1 波长的扫描序列观察 PSF 和 MTF 的变化。离焦量与焦面位移 z 的换算关系是W_defocus (N.A.² / (2λ)) · z其中 N.A. 是数值孔径λ 是波长。比如一个 N.A.0.1 的显微物镜波长为 0.55 μm1 波长离焦大约对应 110 μm 物理位移。跑完离焦扫描后看 MTF 曲线族在零频处所有曲线都归一为 1低频段随离焦量增大迅速跌落中频段出现零点或凹陷。通过对这些 MTF 曲线的包络做进一步处理可以近似还原出系统的“离焦容限”参数。离焦扫描还能验证代码内部的正负符号是否正确理想球面波在焦点前后对称离焦时PSF 严格对称仅中心亮暗变化位置一致如果发现正负离焦的 PSF 不对称那么离焦项符号或者说相位计算里正负号定义就有错误。5. 孔径形状、采样不足与坐标原点三处影响精度的细节这一节把几个工程中反复踩到的坑集中处理同时给出每次算完 PSF 后必须做的三项验证。这三项验证不花时间但能拦住大部分低级错误。第一项验证是孔径边缘检查——画出光瞳函数的实部图确认圆孔径边缘落在预期的像素位置且没有出现伪影。如果孔径边缘出现细密的干涉条纹说明 FFT 之前的孔径边缘过陡可以直接把边缘过度带加宽。第二项是总能量守恒——修改像差量前后PSF 的总能量即psf.sum()必须保持不变归一化前由 Parseval 定理保证。第三项是MTF 横轴单位换算——MTF 图像的像素间隔对应的是空间频率分辨率具体公式是Δf 1 / (N · Δx_eff)其中 Δx_eff 是 PSF 输出平面像素对应的实际物理尺寸在设定 pupil_radius 时已经隐式决定。如果你的系统用了焦距 f 和光瞳直径 DPSF 平面像素间隔和光瞳平面像素间隔之间的关系是Δx_psf λf / (N·Δx_pupil)换算 MTF 横轴为 lp/mm 时用这个关系导出即可。坐标原点问题经常出现在从外部文件读入 Zernike 系数时。圆孔径的圆心必须落在网格中心如果网格点数是偶数比如 512中心其实在相邻四个像素的交界处。fftshift会把这个位置安排到 (N/2, N/2) 点。如果你从某个工具箱导入网格它把原点放在 (0,0) 角点zernike 多项式的奇偶项符号会全部错乱。排查方法只加倾斜项如 n1, m1看 PSF 是否严格沿某个方向平移不变形。若 PSF 变成非对称形状而不是整体偏移原点定义一定错了。最后留一个实用技巧在调试阶段把离焦和球差同时设为 0.2 波长算出的 PSF 应该呈现“胖中心 一环较亮圆环”的结构。这个组合对大部分坐标和排序错误非常敏感——只要看到这个典型形态说明整套计算链路基本正确随后再做扫参和批量计算效率会高很多。本文还有配套的精品资源点击获取
从Zernike系数计算PSF与MTF：Python光学仿真实现指南
简介在光学成像系统分析中点扩散函数、调制传递函数与泽尼克多项式是评价像差与成像质量的核心工具。面向光学设计、图像处理和波前模拟学习者该资源提供了一套从泽尼克多项式拟合到点扩散函数与调制传递函数计算的可运行脚本并配有讲解用幻灯片和网页讲义能够帮助使用者对照公式理解代码逻辑、快速掌握光学仿真流程。压缩包共18个文件以9个M脚本为主涵盖波前像差计算、点扩散函数生成、调制传递函数曲线绘制等核心模块另有若干示意图、幻灯片及说明性文档辅助学习整体仅1.09MB轻量但完整适合直接运行调试。已有2210人学习下载对于希望在光学仿真中减少从零摸索成本、系统理解PSF与MTF理论的研究者和工程师来说是一份实用的上手资料。1. 计算点扩散函数是在算什么从PSF到MTF再到波前像差计算点扩散函数PSF这件事表面上是一个傅里叶变换实际上是把一个光学系统的全部脾气摸清楚。一个理想光学系统点光源成像后应该还是一个点但衍射和像差会让这个点摊开成一团光斑——这团光斑的强度分布就是 PSF对它做傅里叶变换取模得到的就是 MTF。而 Zernike 多项式在这里扮演的角色是用一组正交基函数把波前像差拆成看得懂的成分离焦、球差、彗差、像散每一项都有明确的物理含义和可调系数。做机器视觉、工业镜头选型、光刻照明系统设计或者自适应光学的人都绕不开这套计算流程。市面上很多商业软件帮你把 PSF 和 MTF 算好了但如果你手里只有一堆实测的 Zernike 系数或者需要把自定义孔径、非圆对称像差加进系统里评估自己动手算一遍就是唯一的可靠路径。这篇文章从理论公式一路推到可运行的 Python 代码和参数调试聚焦在“怎么从 Zernike 系数算出 PSF 和 MTF以及算了以后怎么确信算对了”。2. Zernike 多项式的物理意义与波前表示从光圈坐标到像差分量用 Zernike 多项式来描述波前像差最早是 Frits Zernike 在 1934 年提出的核心理念是在单位圆内定义一组互相正交的二维多项式任何连续波前都可以表示成这些多项式的线性组合。这个思路跟傅里叶级数类似区别在于 Zernike 定义在圆域上而光学系统的光瞳天然就是圆的。波前 W(x, y) 可以表示为W(x, y) Σ c_i · Z_i(x, y)其中 Z_i 是第 i 项 Zernike 多项式c_i 是该项的系数单位通常取“波长”waves。每一项 Z_i 由角向频率 m 和径向阶数 n 决定展开来写就是径向函数和角向函数的乘积Z_i(r, θ) R_n^m(r) · sin(mθ) 或 cos(mθ)r 是归一化径向坐标光瞳边缘为 1θ 是极角。径向函数 R_n^m(r) 的形式是固定的多项式组合网上随处可以查到闭式表达式这里不再贴冗长的公式。真正容易混淆的是排序方式Zernike 多项式常用的排序有 Noll 序号、Fringe 序号、OSA 标准序号三种同一项在不同排序里的编号完全不同混用是新手最常见的翻车点。比如离焦项Fringe 序号是 4对应 z[4]Noll 序号是 4一致但彗差两项就不一样了——Fringe 里是 7 和 8Noll 里是 7 和 8也有差别像散项差得更多。如果从实验设备直接导出的 Zernike 系数务必先确认设备用的是哪一种排序。2.1 前几阶 Zernike 项与像差的对应关系实际光学设计里低阶像差解释了大半的成像问题高阶项通常只在强离轴系统或自由曲面系统里才需要关注。下表列出 Fringe 排序下前 15 项对应的像差名称和典型物理来源写计算程序之前先把这张表放在手边序号nm像差名称物理来源典型符号约定100活塞常数相位不影响 PSF可忽略211倾斜 X光斑位置偏移符号决定偏移方向31-1倾斜 Y光斑位置偏移符号决定偏移方向420离焦对焦不准正负对应焦点前后52-2像散 45°柱面/倾斜面与第6项正交622像散 0°柱面元件/装配应力参考轴选择731彗差 X偏心或光轴不重合与倾斜方向相关83-1彗差 Y偏心或光轴不重合同上933三叶草镜片制造误差三项对称103-3三叶草镜片支撑应力与第9项方向正交1140初级球差球面本身固有像差正为边缘聚焦更近1242高阶像散高次曲面残差少见134-2高阶像散同上少见1444四叶草车削残留高次面154-4四叶草装配应力高次面工程上有个非常常见的做法把 Zernike 系数理解为“对应的 Seidel 像差量”即某一项系数值 0.25 就意味着该像差的波前 RMS 贡献约 0.25 波长。不过严格说只有把波前展开后按项做 RMS 统计才算准确因为多项式虽正交但同阶项叠加后的 RMS 与分项系数之间存在勾股关系RMS_total sqrt(Σ c_i²)这个性质在很多相位反演算法里被直接用来约束解空间。做仿真时如果只想调一种像差直接把对应项的系数设成非零值即可其他项保持 0。2.2 为什么用 Zernike 而不是直接搞一个波前函数实践里有人会把波前直接写成一个简单曲面比如球面波前用 W A·ρ²然后去算 PSF这样算出来的结果并不严谨。Polynomial 系数一旦换成物理坐标就丢掉了两个关键性质——正交性和旋转对称性。Zernike 的每一项在圆域内正交意味着调整某一项系数不会改变其他项对波前的贡献这在优化过程里极其重要你调彗差系数不会鬼使神差地影响离焦量的评估。另外Zernike 多项式的旋转对称性角向频率 m对应光学系统的对称特征。m 0 的项旋转对称m 1 的项是“蝴蝶形”m 2 的像散在旋转 90° 后取反号。利用这个性质可以快速判断仿真结果是否符合物理直觉——比如旋转对称的球差项算出来的 PSF 必须是旋转对称的如果不对称那一定是采样网格或坐标定义出了错。这个特性后面章节里做验证时会反复用到。3. 用 Python 从 Zernike 系数算出 PSF 和 MTF 的最小可运行实现光学系统的 PSF 计算基于标量衍射理论核心公式是光瞳函数的傅里叶变换模平方PSF(x_f) |FFT{P(x_p)}|²其中 P(x_p) 是光瞳平面上的复振幅分布由振幅透过率 A(x_p) 和相位项组合而成P(x_p) A(x_p) · exp(j · (2π/λ) · W(x_p))A 代表孔径形状圆孔内部为 1外部为 0W 是波前像差单位波长直接取 Zernike 展开结果。MTF 反过来是 PSF 的傅里叶变换模归一化到零频后取模MTF(f) |FFT{PSF(x_f)}| / |FFT{PSF(x_f)}|_(f0)需要注意一个关键点这里 PSF 的单位网格要和频率空间网格匹配即最后输出的 MTF 横轴是用像素表示的采样频率要换算成物理频率lp/mm必须知道实际系统里的像素缩放关系。3.1 网格生成与 Zernike 多项式函数实现网格生成这一步直接影响计算结果的精度。工程上最常用的是“奇偶数采样”方案采样点数 N 取偶数坐标范围从 -1 到 1 取 N 个点这样避免了原点落在网格中心时 FFT 常见的偏移问题。Zernike 多项式函数采样的传统写法是定义于单位圆内——也就是把物理光瞳半径缩放到 1采样完以后把 r 1 的位置掩膜掉。下面这段实现把多项式生成从 ANSI C 常见写法改成 numpy 向量化实现效率高且不容易错import numpy as np def zernike_poly(n, m, rho, theta): 计算单个 Zernike 多项式在极坐标网格上的值 n: 径向阶数, m: 角向频率(带符号) rho: 归一化径向坐标 (0~1), theta: 方位角 返回: 与 rho 同形状的浮点数组 # R_n^m 径向多项式用递推公式 # 先处理 m0 的特殊情况 if m 0: # p 为多项式的半阶 s_max n // 2 else: s_max (n - abs(m)) // 2 radial np.zeros_like(rho) for s in range(s_max 1): # 组合数直接用阶乘计算避免引入 scipy 依赖 coef ((-1) ** s) * np.math.factorial(n - s) / ( np.math.factorial(s) * np.math.factorial((n abs(m)) // 2 - s) * np.math.factorial((n - abs(m)) // 2 - s) ) radial coef * (rho ** (n - 2 * s)) # 角向分量: 根据 m 的符号选择 cos 或 sin if m 0: angular np.cos(m * theta) elif m 0: angular np.sin(-m * theta) else: angular np.ones_like(theta) norm np.sqrt(2 * (n 1)) if m ! 0 else np.sqrt(n 1) return radial * angular * norm代码逻辑说明s_max决定多项式迭代的项数radial部分由 R 多项式的标准求和式算出angular按 m 符号选择三角基函数——这个符号约定对应“径向多项式为正、角向用 sin/cos 区分方向”的常用形式与 Zemax 导出的数据符号规则一致。norm系数使得多项式在单位圆上 RMS 等于 1这样 Zernike 系数就可以直接代表 RMS 值实测数据里导出的系数如果不做 RMS 归一化这一步会额外引入一个sqrt(2)或sqrt(n1)的差异。调用时只需构造极坐标网格N 512 # 采样点数偶数 x np.linspace(-1, 1, N, endpointFalse) X, Y np.meshgrid(x, x) rho np.sqrt(X**2 Y**2) theta np.arctan2(Y, X)这里endpointFalse很关键配合偶数 N网格点关于原点对称FFT 以后的频谱中心不会偏移半个像素。3.2 从波前相位到 PSF 再到 MTF 的完整计算函数把相位、孔径、傅里叶变换串起来核心计算函数可以压缩成下面这样def compute_psf_mtf(zernike_coeffs, wavelength0.55e-3, N512, pupil_radius1.0, pixel_scale1.0): 输入: zernike_coeffs: 列表, 按单位圆内 RMS 归一化的 Zernike 系数 长度自适应, 缺失项按 0 处理 wavelength: 波长, 单位 mm (可见光约 0.00055) N: 采样网格边长 (偶数) pupil_radius: 光瞳半径, 决定孔径边缘落在哪个像素 pixel_scale: 输出 PSF 的像素缩放, 方便后续换算 MTF 返回: psf: N×N 强度分布, 已归一化到总能量 1 mtf: N×N 调制传递函数, 零频归一化为 1 W: 波前图 (单位: 波长), 用于调试可视化 # 1. 极坐标网格 x np.linspace(-1, 1, N, endpointFalse) X, Y np.meshgrid(x, x) rho np.sqrt(X**2 Y**2) theta np.arctan2(Y, X) # 2. 振幅孔径: 半径为 pupil_radius 的圆, 外部置零 aperture (rho pupil_radius).astype(float) # 3. 累加 Zernike 波前 (单位: 波长) W np.zeros_like(rho) for idx, coef in enumerate(zernike_coeffs): if coef 0: continue n int(np.sqrt(idx 1)) - 1 # 这个索引是近似, 见下文说明 # 实际使用时建议传入 (n, m) 对, 这里简化 m 0 # 占位, 需要按排序映射实际 m Z zernike_poly(n, m, rho, theta) W coef * Z # 4. 构造复振幅光瞳函数 phase (2 * np.pi / wavelength) * W pupil aperture * np.exp(1j * phase) # 5. FFT 算 PSF: 光瞳函数投影到焦平面 psf np.abs(np.fft.fftshift(np.fft.fft2(pupil))) ** 2 # 归一化总能量为 1 psf / psf.sum() # 6. MTF: PSF 再做一次 FFT, 模归一化 mtf np.abs(np.fft.fftshift(np.fft.fft2(np.fft.ifftshift(psf)))) mtf / mtf[0, 0] # 零频归一化 return psf, mtf, W代码里第 3 步有一个关键的简化从数组索引倒推 (n, m) 并不是一个可靠的做法依赖特定排序表。工程上最稳的方式是在函数入口处直接传入(idx, n, m, coef)四元组列表或者用一个标准 Zernike 排序表映射。上面代码只是展示整体数据流结构真用的时候建议把(n, m, coef)作为显式输入。FFT 两步分别用了fftshift和ifftshift的组合这是为了避免光瞳网格的偏移误差第一步fft2之后用fftshift把零频挪到中心第二步由于psf已经是中心化过的数据做 FFT 之前要先ifftshift还原到 FFT 算法的原点布局。这两个函数混用错了一个MTF 会整体错位而且很难从视觉上察觉。3.3 离焦与球差的仿真验证用一个 20 行以内的驱动脚本把上面函数跑起来# 只有 0.25 波长的初级球差 (n4, m0, 系数 0.25) z_coeffs [(4, 0, 0.25)] W np.zeros((N, N)) for n, m, c in z_coeffs: Z zernike_poly(n, m, rho, theta) W c * Z psf, mtf, W compute_psf_mtf_from_phase(W, wavelength0.55e-3, N512)如果只加了球差PSF 结果应该是中心亮斑周围带对称的同心环结构沿着光轴两侧离焦环结构会出现非对称的明暗变化这是球差“焦点位移导致模糊不对称”的典型特征。MTF 在某个中间频率处可能出现凹陷甚至落到零再回升这就是经典的“焦点外 MTF 带有频带缺口”现象做机器视觉镜头评测时如果看到这种 MTF 形状可以直接判断系统有残留球差。4. 像差参数怎么设系数单位、采样率与 MTF 曲线验证跑通第一版代码之后下一步是把仿真结果调到跟实验或设计数据对得上。这个阶段最磨人的不是代码逻辑而是参数设定的一致性。下面三个参数是最常出问题的。4.1 Zernike 系数单位与符号约定Zernike 系数的物理单位有两种习惯一种是以波长为单位的波前值RMS另一种是直接给出“光程差”OPD的单位是微米或纳米。两者换算关系是OPD coeff × wavelength。很多从商用软件导出的系数写的是“waves RMS”拿到手直接用就行如果是干涉仪采出的原始数据单位很可能是微米级 OPD需要除波长。符号约定方面各软件用的坐标系定义不同Zemax 和 Code V 对同一项 Zernike 的符号就可能相反。节省时间的做法是先用单一大像差如 1 波长离焦跑一次仿真确认 PSF 的模糊方向和像差符号的对应关系再拿实验数据对照一次后面大批量数据就按这个映射关系校准。4.2 采样网格点数与孔径边缘像素对精度的影响网格 N 的选择是一个典型的精度-速度权衡。FFT 计算的时间复杂度是 O(N² log N²)N 翻倍速度掉四倍。工程上用 N512 起步绝大多数单视场点扩散函数计算都在一秒钟之内完成。真正影响精度的是光瞳边缘落在哪几个像素上孔径圆边界的锯齿状量化误差会造成 PSF 高频部分的虚假能量表现是 MTF 在高频区域出现“裙边”上翘。缓解办法是给孔径边缘加过渡带使用超采样反走样。一段常用的小技巧# 孔径边缘用 2 像素平滑过渡, 减少振铃 edge 2.0 # 过渡带宽度 (像素) aperture 0.5 * (1 - np.tanh((np.abs(rho) - pupil_radius) / edge))把硬边缘改成软边缘之后PSF 外围的伪振荡幅度明显下降代价是中心强度略微下降、斯特列尔比Strehl Ratio计算值会比真实值低约 0.5% 以下可以接受。高频 MTF 的精度收益通常远大于这个损失。4.3 离焦扫描验证MTF 曲线随离焦量的变化规律一个典型的工程验证场景是“离焦扫描”给 Zernike 系数中离焦项n2, m0设置从 -1 到 1 波长的扫描序列观察 PSF 和 MTF 的变化。离焦量与焦面位移 z 的换算关系是W_defocus (N.A.² / (2λ)) · z其中 N.A. 是数值孔径λ 是波长。比如一个 N.A.0.1 的显微物镜波长为 0.55 μm1 波长离焦大约对应 110 μm 物理位移。跑完离焦扫描后看 MTF 曲线族在零频处所有曲线都归一为 1低频段随离焦量增大迅速跌落中频段出现零点或凹陷。通过对这些 MTF 曲线的包络做进一步处理可以近似还原出系统的“离焦容限”参数。离焦扫描还能验证代码内部的正负符号是否正确理想球面波在焦点前后对称离焦时PSF 严格对称仅中心亮暗变化位置一致如果发现正负离焦的 PSF 不对称那么离焦项符号或者说相位计算里正负号定义就有错误。5. 孔径形状、采样不足与坐标原点三处影响精度的细节这一节把几个工程中反复踩到的坑集中处理同时给出每次算完 PSF 后必须做的三项验证。这三项验证不花时间但能拦住大部分低级错误。第一项验证是孔径边缘检查——画出光瞳函数的实部图确认圆孔径边缘落在预期的像素位置且没有出现伪影。如果孔径边缘出现细密的干涉条纹说明 FFT 之前的孔径边缘过陡可以直接把边缘过度带加宽。第二项是总能量守恒——修改像差量前后PSF 的总能量即psf.sum()必须保持不变归一化前由 Parseval 定理保证。第三项是MTF 横轴单位换算——MTF 图像的像素间隔对应的是空间频率分辨率具体公式是Δf 1 / (N · Δx_eff)其中 Δx_eff 是 PSF 输出平面像素对应的实际物理尺寸在设定 pupil_radius 时已经隐式决定。如果你的系统用了焦距 f 和光瞳直径 DPSF 平面像素间隔和光瞳平面像素间隔之间的关系是Δx_psf λf / (N·Δx_pupil)换算 MTF 横轴为 lp/mm 时用这个关系导出即可。坐标原点问题经常出现在从外部文件读入 Zernike 系数时。圆孔径的圆心必须落在网格中心如果网格点数是偶数比如 512中心其实在相邻四个像素的交界处。fftshift会把这个位置安排到 (N/2, N/2) 点。如果你从某个工具箱导入网格它把原点放在 (0,0) 角点zernike 多项式的奇偶项符号会全部错乱。排查方法只加倾斜项如 n1, m1看 PSF 是否严格沿某个方向平移不变形。若 PSF 变成非对称形状而不是整体偏移原点定义一定错了。最后留一个实用技巧在调试阶段把离焦和球差同时设为 0.2 波长算出的 PSF 应该呈现“胖中心 一环较亮圆环”的结构。这个组合对大部分坐标和排序错误非常敏感——只要看到这个典型形态说明整套计算链路基本正确随后再做扫参和批量计算效率会高很多。本文还有配套的精品资源点击获取