从Zernike系数计算PSF与MTF:Python光学仿真实现指南
2026/9/20 11:19:18 网站建设 项目流程

简介:在光学成像系统分析中,点扩散函数、调制传递函数与泽尼克多项式是评价像差与成像质量的核心工具。面向光学设计、图像处理和波前模拟学习者,该资源提供了一套从泽尼克多项式拟合到点扩散函数与调制传递函数计算的可运行脚本,并配有讲解用幻灯片和网页讲义,能够帮助使用者对照公式理解代码逻辑、快速掌握光学仿真流程。压缩包共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 和 8,Noll 里是 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,外部为 0),W 是波前像差(单位:波长,直接取 Zernike 展开结果)。MTF 反过来是 PSF 的傅里叶变换模,归一化到零频后取模:

MTF(f) = |FFT{PSF(x_f)}| / |FFT{PSF(x_f)}|_(f=0)

需要注意一个关键点:这里 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 径向多项式用递推公式 # 先处理 m=0 的特殊情况 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(n+1)的差异。

调用时只需构造极坐标网格:

N = 512 # 采样点数,偶数 x = np.linspace(-1, 1, N, endpoint=False) X, Y = np.meshgrid(x, x) rho = np.sqrt(X**2 + Y**2) theta = np.arctan2(Y, X)

这里endpoint=False很关键:配合偶数 N,网格点关于原点对称,FFT 以后的频谱中心不会偏移半个像素。

3.2 从波前相位到 PSF 再到 MTF 的完整计算函数

把相位、孔径、傅里叶变换串起来,核心计算函数可以压缩成下面这样:

def compute_psf_mtf(zernike_coeffs, wavelength=0.55e-3, N=512, pupil_radius=1.0, pixel_scale=1.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, endpoint=False) 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 两步分别用了fftshiftifftshift的组合,这是为了避免光瞳网格的偏移误差:第一步fft2之后用fftshift把零频挪到中心,第二步由于psf已经是中心化过的数据,做 FFT 之前要先ifftshift还原到 FFT 算法的原点布局。这两个函数混用错了一个,MTF 会整体错位,而且很难从视觉上察觉。

3.3 离焦与球差的仿真验证

用一个 20 行以内的驱动脚本把上面函数跑起来:

# 只有 0.25 波长的初级球差 (n=4, m=0, 系数 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, wavelength=0.55e-3, N=512)

如果只加了球差,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 翻倍,速度掉四倍。工程上用 N=512 起步,绝大多数单视场点扩散函数计算都在一秒钟之内完成。真正影响精度的是光瞳边缘落在哪几个像素上:孔径圆边界的锯齿状量化误差会造成 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 系数中离焦项(n=2, m=0)设置从 -1 到 1 波长的扫描序列,观察 PSF 和 MTF 的变化。离焦量与焦面位移 z 的换算关系是:

W_defocus = (N.A.² / (2λ)) · z

其中 N.A. 是数值孔径,λ 是波长。比如一个 N.A.=0.1 的显微物镜,波长为 0.55 μm,1 波长离焦大约对应 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 和光瞳直径 D,PSF 平面像素间隔和光瞳平面像素间隔之间的关系是Δx_psf = λf / (N·Δx_pupil),换算 MTF 横轴为 lp/mm 时用这个关系导出即可。

坐标原点问题经常出现在从外部文件读入 Zernike 系数时。圆孔径的圆心必须落在网格中心,如果网格点数是偶数(比如 512),中心其实在相邻四个像素的交界处。fftshift会把这个位置安排到 (N/2, N/2) 点。如果你从某个工具箱导入网格,它把原点放在 (0,0) 角点,zernike 多项式的奇偶项符号会全部错乱。排查方法:只加倾斜项(如 n=1, m=1),看 PSF 是否严格沿某个方向平移不变形。若 PSF 变成非对称形状而不是整体偏移,原点定义一定错了。

最后留一个实用技巧:在调试阶段,把离焦和球差同时设为 0.2 波长,算出的 PSF 应该呈现“胖中心 + 一环较亮圆环”的结构。这个组合对大部分坐标和排序错误非常敏感——只要看到这个典型形态,说明整套计算链路基本正确,随后再做扫参和批量计算,效率会高很多。

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

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询