1. 链路拆解Zernike系数到MTF曲线到底经历了什么做光学的人大概都绕不开这三个词Zernike多项式、PSF、MTF。搞镜头设计的时候要抠Zernike系数做图像算法的要问PSF长什么样做系统验收的要盯MTF曲线。但真要把给一组Zernike系数变成屏幕上一条能评估像质的MTF曲线中间那几步涉及什么物理和数值细节很多朋友其实没完整跑通过。这篇文章就把这条链路完全拆开从波前像差开始一步一步算到PSF再算到MTF代码和参数都给到位。先说这套模拟能干什么。给定一套光瞳上的Zernike像差系数比如离焦0.25个波长、彗差0.3个波长它就能告诉你这个系统的光斑能量分布变成什么样了、不同空间频率的对比度还剩下多少。不管是做光学设计前期的像差预算分析还是做图像复原时构造PSF模板还是上课做傅里叶光学实验这套链路都是最基础的工具。适合的人也很明确光学工程方向的学生、刚入门的光学设计工程师、以及做成像算法的朋友都可以从中拿到一套可以直接抄的作业。我把整个模拟过程分成几个环节波前相位生成、光瞳函数构造、傅里叶变换求PSF、再傅里叶变换求MTF、最后是曲线提取与评价。但其实里面每个环节都有一些坑比如坐标网格怎么对齐、FFT之后坐标轴怎么换算、MTF要不要归一化。这些细节网上的教程很少讲全踩过坑的都知道往往一个像素没对齐结果就完全不对。1.1 为什么要以Zernike系数为输入Zernike多项式是定义在单位圆上的一组正交基函数它最大的特点是把复杂的波前像差分解成一个个具有明确物理意义的模式。低阶项对应平移、倾斜、离焦高阶项对应球差、彗差、像散等等。光学设计软件里经常直接把镜面面形误差或系统像差展开成Zernike系数因为每一项跟赛德尔像差都有明确的对应关系又能直接跟干涉仪测出来的波前图对上。做模拟时以Zernike系数为输入最实际的好处是参数直观。你跟设计师说这个镜头有0.3波长的三阶彗差他脑子里立刻有画面。你要是直接塞一个二维相位数组给他他没法判断这个相位分布对应什么像差。所以无论是验收还是返工Zernike系数都是行业里通用的语言。不过要注意一个关键问题Zernike归一化约定是有讲究的。常用的有Noll序、Fringe序和标准Zernike序不同约定下同一个系数对应的多项式形式和归一化系数不一样。很多人在这一步栽过跟头以为拿到的是同一个系数结果算出来的PSF完全对不上。后面我会专门讲这块怎么避坑。1.2 PSF和MTF在这个链路里各自扮演什么角色PSF的全称是点扩散函数可以理解为一个理想点光源经过光学系统后在像面上弥散成什么样的光强分布。衍射受限系统再怎么理想光斑也不会是无穷小它会形成一个以艾里斑为核心的分布。像差存在的意义就是把这个理想光斑弄得更复杂中心能量降低、能量往外部环带扩散、甚至出现不对称的结构。MTF则是调制传递函数它衡量的是系统对不同空间频率的正弦光栅的对比度传递能力。MTF曲线的低频段决定了大目标的对比度高频段决定了细节能不能分辨曲线下的面积在某种程度上反映了整体成像质量。MTF跟PSF在数学上是一对傅里叶变换关系PSF做一次傅里叶变换取模就得到了OTF光学传递函数再取模就是MTF。所以这条链路的本质就是Zernike系数定义了波前相位波前相位决定了衍射积分的相位分布衍射积分计算出PSFPSF再做一次傅里叶变换得到MTF。每一步都是上一步的输入没有中间环节可以跳过去。理解了这条因果链再去写代码就顺了。1.3 整个模拟方案的选型逻辑计算PSF的标准方法是用傅里叶光学里的夫琅禾费衍射公式。对于圆形光瞳理论上可以直接做贝塞尔函数积分但那只适用于无像差的情况。一旦引入任意Zernike像差解析解基本不存在数值方法几乎是唯一选择。最常见的做法是把光瞳采样成N×N的网格光瞳函数写成振幅掩膜乘以相位因子然后做一次二维FFT取模平方得到PSF。这套方案的优势很明显实现简单、速度快、可以任意扩展像差模式。FFT的本质是离散傅里叶变换它计算的是周期性延拓后的结果因此采样参数直接决定了计算准确性。N取多少、网格范围怎么定义、光瞳圆域怎么用掩膜切出来这些细节都要提前想清楚不然FFT一跑出来的PSF可能全是混叠和振铃。我没有选更复杂的角谱法或严格衍射积分原因很简单对于光瞳尺寸远大于波长的普通成像系统这种基于FFT的衍射计算精度已经足够了。只有当系统进入深亚波长尺度或需要考虑偏振效应时才需要上严格矢量衍射模型。初学者先跑通这套标量模型是性价比最高的路径。2. 核心细节像差数值化的几个关键点2.1 Zernike多项式的归一化约定与单位Zernike多项式由径向函数和角向函数相乘构成。标准形式里每一项用三个参数标记径向阶数n、角向频率m、以及系数a_nm。相位分布Φ(r,θ)就是这些项的线性叠加。单位通常用波长表示系数a_nm为1就代表该项引入1个波长的光程差。如果你自己写实现可以不用理会复杂的索引表直接按递推关系生成径向多项式。但前提是心里要清楚你在用哪套定义。比如Noll序把第一项叫pistonFringe序里第一项却是piston加倾斜的组合。归一化系数也有区别有的约定里像散项的峰谷值就是系数的两倍有的约定里系数的RMS值等于1。两种约定下同一个物理像差写出来的系数甚至差出两倍多。我写模拟时统一用标准ZernikeBorn Wolf版即每项在单位圆上积分的RMS为1除piston外。这样做的好处是系数大小直接跟RMS误差挂钩评价像质时方便换算。代码里我会直接按这个约定实现radial多项式然后归一化每一项再把系数乘上去。注意如果后续要和Zemax或Code V对标务必查清楚他们导出的Zernike系数用的是哪一套约定否则结果没法直接比较。2.2 光瞳网格采样N512够不够FFT计算PSF的精度很大程度上取决于光瞳网格的采样密度。网格太稀光瞳边缘的圆形边界会变成锯齿状高频衍射分量被污染网格太密内存和计算时间又会飙升。我常用的经验值是N512。对于圆形光瞳这个采样量足以保证PSF中心到第三四个暗环的形态都准确MTF曲线也能平滑延伸到截止频率。如果是快速验证N256也能用但MTF高频段会开始出现可见的噪声毛刺。N1024适合做最终出图特别是在要放大看PSF细节的时候。还有一个容易忽略的点网格的物理尺寸定义。通常把光瞳半径归一化为1网格坐标从-1到1步长2/N。这样FFT之后输出平面的角度尺度是1/(2Δ)即N/4个衍射单位每像素。很多教程只给代码不给换算关系导致PSF图像的横坐标到底对应多少微米完全对不上。这个问题我在经验分享里还会展开。圆域掩膜的处理要特别注意最稳妥的方式是用阈值判断r≤1再乘进去。不要直接用条件索引去修改FFT结果那样会破坏数组形状。掩膜边缘如果出现0.5级别的灰度过渡即抗锯齿处理PSF的振铃会小很多但也会稍微平滑真实衍射环看需求取舍。2.3 从波前到PSF的FFT原理与参数换算夫琅禾费衍射的数值实现可以写成先构造光瞳函数P(x,y) A(x,y)·exp(i·2π·W(x,y)/λ)其中A是孔径振幅圆内为1圆外为0W是波前像差函数用波长单位。然后对P做二维FFT得到像面上的复振幅分布。PSF |FFT(P)|^2再归一化到总能量为1。为什么可以直接用FFT因为FFT计算的是离散傅里叶变换而夫琅禾费衍射本质上就是光瞳函数的傅里叶变换。只要采样满足奈奎斯特条件离散结果就能很好地近似连续积分。这里有一个隐含假设系统满足远场条件且像差导致的相位变化在光瞳面上缓慢变化。对于普通镜头设计这个假设完全成立。相位因子写成2πW/λ还是直接写成2π·系数之和取决于W的单位。我在代码里统一把Zernike系数定义为波长数所以相位因子直接乘2π。比如离焦系数0.25λ等效相位幅度就是1.57弧度。这样从物理到代码单位链条是闭合的。MTF的计算同样简单把PSF再做一次FFT取模并归一化到零频为1就得到MTF。注意MTF的x轴频率单位是每弧度lp还是归一化频率要和PSF的像素尺度联动换算。一般的规则是如果PSF横轴是像素索引那么MTF横轴归一化空间频率的最大值对应奈奎斯特频率截止频率通常出现在约0.5/ciclo附近具体位置由NA和波长决定。3. 实操过程一套可以复现的Python实现3.1 环境准备与工具选型工具组合很简单Python NumPy Matplotlib。这三个库是科学计算标配不需要额外安装光学专用包。如果你想快速验证但不写代码MATLAB的fft2也是一样的操作逻辑但后期做批量和自动化不如Python方便。我用Python还有一个原因是生态完整后续要做图像复原、深度学习PSF估计直接在同一套环境里衔接。另外NumPy 2.0以后的FFT接口更快多核环境下自动并行对于N1024这种规模不会卡顿。安装命令就不多写了标准的三件套。要是你用的是Anaconda环境直接pip install numpy matplotlib就行。代码兼容Python 3.9。3.2 第一步生成Zernike波前相位图先写一个生成Zernike多项式的函数。径向多项式部分用组合数实现递推角向部分分别用cos和sin处理正负m。代码里我按标准Zernike约定归一化了RMS每一项在单位圆上积分的平方均值归一为1。import numpy as np import matplotlib.pyplot as plt def zernike_radial(n, m, r): 计算Zernike径向多项式 R_n^m(r), r 是归一化半径数组 m abs(m) if (n - m) % 2 ! 0: return np.zeros_like(r) terms [] for s in range((n - m) // 2 1): coeff ((-1) ** s * np.math.comb(n - s, s) * np.math.comb(n - 2 * s, (n - m) // 2 - s)) terms.append(coeff * r ** (n - 2 * s)) return np.sum(terms, axis0) def zernike(n, m, r, theta): 标准Zernike多项式返回在(r,theta)网格上的值 radial_val zernike_radial(n, m, r) if m 0: return radial_val * np.cos(m * theta) elif m 0: return radial_val * np.sin(-m * theta) else: return radial_val def make_grid(N512): 生成归一化坐标网格和圆形孔径掩膜 x np.linspace(-1, 1, N) xx, yy np.meshgrid(x, x) r np.sqrt(xx**2 yy**2) theta np.arctan2(yy, xx) aperture (r 1.0).astype(float) return xx, yy, r, theta, aperture有了这三个函数生成波前相位就很直白了。设定一组系数比如离焦项(2,0)系数0.25彗差项(3,1)系数0.2像散项(2,2)系数0.1然后线性叠加xx, yy, r, theta, aperture make_grid(512) coeffs [ (0, 0, 0.1), # piston一般可以忽略 (2, 0, 0.25), # 离焦 0.25λ (2, 2, 0.1), # 像散 0.1λ (3, 1, 0.2), # 彗差 0.2λ ] phase np.zeros_like(r) for n, m, c in coeffs: phase c * zernike(n, m, r, theta) # 光瞳内有效相位 phase_masked phase * aperture跑完这一小段用plt.imshow输出相位图你应该能看到一个既有四分之一弧形的离焦轮廓、又有不对称彗差特征的图案。这一步跑通了后面就是纯粹的FFT变换。3.3 第二步计算PSF并正确显示构造复光瞳函数做FFT取模平方。核心就四行pupil aperture * np.exp(2j * np.pi * phase) psf np.abs(np.fft.fftshift(np.fft.fft2(np.fft.fftshift(pupil))))**2 psf / psf.sum()这里面fftshift的用法有点讲究我习惯在FFT前后各做一次shift这样得到的PSF中心在数组中央方便观察。你也可以只在外侧做一次但那样中心在角落看图会非常别扭。PSF的动态范围很大中心峰值可能比边缘高出好几个数量级直接imshow会变成一团白点。正确的显示方式是取对数并做饱和度裁剪psf_log np.log10(psf 1e-9) plt.imshow(psf_log, cmapinferno, vminpsf_log.max()-6, vmaxpsf_log.max())vmin取最大值往下6个数量级基本能看到外圈衍射环的整个结构。如果你要量化中心能量直接从归一化后的psf数组里取中心点的值即可它其实就是斯特列尔比Strehl Ratio的近似。相位缠绕不需要解包因为exp(2jπ×0.1)和exp(2jπ×1.1)是完全一样的值。真正需要担心的反而是光瞳外的NaN——如果光瞳掩膜直接把r1置零而你在phase_masked里又对r0做了除零操作结果会出现nan沿着孔径边缘扩散。代码里已经规避了这个aperture乘在最后。3.4 第三步提取MTF曲线MTF的算法已经说过重点在于提取曲线和坐标换算。先算出二维MTFotf np.fft.fftshift(np.fft.fft2(np.fft.fftshift(psf))) mtf np.abs(otf) mtf / mtf.max() # 零频归一化到 1二维MTF图可以直接imshow看但工程上更常用一维剖面。对于圆对称系统直接取通过中心的水平切面或垂直切面都行。对于有像散或彗差的系统不同方位角的MTF不一样所以一般会分别取两个正交方向的切面。freq np.fft.fftshift(np.fft.fftfreq(N, d2.0/N)) profile_x mtf[N//2, :] profile_y mtf[:, N//2] plt.plot(freq, profile_x, labeltangential) plt.plot(freq, profile_y, labelsagittal) plt.xlim(0, 0.5)注意fftfreq的d参数即采样间距。因为光瞳网格从-1到1共N个点所以实际步长是2/N。FFT后的频率轴的刻度单位是1/像素对应的空间频率上限是0.5 cycle/pixel奈奎斯特极限。如果你要把横轴换算成实际空间频率lp/mm需要知道光瞳的实际尺寸、波长和焦距换算公式是f_real freq_normalized × (1/(λ×F#))到时候按系统参数乘就行。我给的示例代码里没有做严谨的实际单位换算因为归一化频率已经能反映相对趋势。做工程评估时务必根据你的光学系统参数把横轴换成lp/mm不然光看曲线形状很难判断系统到底达不达标。3.5 实例输出与结果解读拿上面那组系数跑完整流程你会看到两个关键现象一是PSF中心的能量明显下降中心亮斑变小、周围产生不对称的旁瓣二是MTF曲线在中低频段对比无像差的衍射受限曲线有显著下降高频截止位置不变但曲线形状被压低。无像差时MTF是一条接近直线的下降曲线到归一化频率0.5处截止。加入离焦后低频段掉得最快中频段可能出现非单调的凹陷。这是因为离焦的相位是二次型对低频成分的对比度影响特别大。彗差和像散则会让两个方向的MTF不一致这也是为什么实际镜头评测要分切线方向和弧矢方向分别测。如果你把piston项去掉设为0PSF的中心位置不会移动但MTF几乎不变。这说明piston只是整体相位平移不影响成像质量。跑完这组案例后建议自己多调几次系数比如把离焦调到0.5λ、球差调到1λ看看PSF怎么从艾里斑变成甜甜圈形状MTF怎么逐渐失去高频细节。4. 常见问题与排查技巧实录4.1 PSF中心不亮还带条纹先查光瞳掩膜我最初跑这套模拟时遇到最多的就是PSF中心突然变暗而且外面还有一圈竖直的条纹。排查了一圈最后发现是光瞳坐标网格没对中网格的眼神线在像素边界上导致圆形掩膜的左半圈和右半圈差了一列像素。这种情况下PSF肯定出大问题。解决办法很简单生成坐标时用np.linspace(-1, 1, N)并且确保N是偶数。用偶数N时中心点在两列像素之间圆形掩膜对称性最好。用奇数N时中心落在某个像素上虽然有轻微不对称但算PSF也基本可用。我的习惯是N取512或1024都是偶数。另一个bug是掩膜和相位没有同步。比如aperture里光线不经过的区域应该是0但你如果忘了乘aperture光瞳外的相位噪声也会参与FFTPSF就会出现莫名其妙的背景能量。每次构造光瞳函数前我都建议先把光瞳图和相位图并排打印出来看一眼确认掩膜对不对齐。4.2 MTF高频掉太快还是全是毛刺MTF曲线出现毛刺大概率是PSF采样不足导致的。PSF在中心附近的峰值非常陡峭如果网格不够密FFT结果的高频分量就会出现频谱泄漏。解决办法是提高N或者在做FFT之前对光瞳数组做零填充pad到4倍大小。零填充不增加信息量但能让频域插值更平滑曲线上的毛刺会少很多。如果MTF高频掉得比理论值快很多可能是PSF没有正确归一化或者PSF数组里有零值导致log运算出现inf。我调试时会先单独打印psf.min()和psf.max()如果min0取对数时记得加一个小常数1e-9。MTF归一化也有讲究要用mtf[0,0]零频分量作为除数而不是mtf.max()因为数值误差下两者可能略有差异。我一般直接用中心值这样曲线严格从1出发。4.3 相位缠绕和NaN问题相位本身不会缠绕出问题因为复数指数是周期函数。真正麻烦的是在极坐标网格里当r0时theta计算的是arctan(0,0)在NumPy里这个值没问题但Zernike多项式的cos(m·theta)项会变成cos(0)或者cos(π/2)导致中心点数值不确定。处理办法要么单独给r0点赋一个固定值一般取0即可要么在计算theta之前给r加上一个极小值epsilon。我倾向于后一种因为代码简单且不影响物理结果。做完相位计算后再乘上aperture掩膜然后检查数组里是否还有nan。如果发现nan一定是某个Zernike径向多项式的递推在r非常接近1时溢出这时检查径向公式里有没有除以(r-1)之类的项。4.4 Zernike索引与阈值定义混乱这是跟别人对数据时最容易出问题的地方。同一项像差不同软件里的索引号能差出一大截。比如Noll序的第7项是垂直彗差Fringe序里对应的可能是第8项而标准Zernike序里它又是另一项。你拿别人给的Zernike系数表上来就按自己的索引代入结果必然对不上。我的建议是在代码开头明确写入一个索引-物理项对照表用单项重构波前验证核对。比如先生成只有离焦(2,0)的波前看看是不是二次旋转对称的环带分布生成只有彗差(3,1)的波前看看是不是不对称的典型彗形图案。这样即使索引错了一眼就能看出来。另外所有系数的单位统一用波长且要说明RMS还是PV值。工程上常用RMS值因为评价函数更好算但有些供应商给的是PV。5. 像差模拟还能往哪走5.1 用评价函数量化像差跑通PSF和MTF之后下一步自然是量化像质。最常用的三个指标斯特列尔比Strehl Ratio、RMS波前误差、MTF面积分。斯特列尔比直接取PSF中心最大值就行大于0.8通常认为是衍射受限。RMS波前误差可以从Zernike系数算出所有非piston项的系数平方和再开根号就是RMS值前提是系数已按RMS归一化。这三个指标之间是有关联的。经验上RMS误差在0.1λ以下时Strehl比近似等于exp(-(2π·RMS/λ)^2)也就是马雷夏尔判据。把这三个指标加进你的模拟脚本以后做像差预算时就多了一套参考坐标系不用每次只看PSF图猜质量。5.2 从单点扩展到多点与部分相干这套模拟最直接的扩展方向是做离焦扫描、景深分析或全场成像模拟。把离焦系数做成扫描变量得到一系列PSF进一步还能算出系统的离焦MTF曲线族这对自动对焦算法的仿真很有用。另一个方向是把单点的PSF扩展到多个场角每个场角取不同的Zernike系数组合从而模拟真实镜头的像差场依赖性。对于照明条件复杂的情况比如部分相干成像PSF不再是简单的相干叠加需要引入TC和交叉谱密度复杂度会高一个量级。但从Zernike系数到波前相位这一步完全不变。先把这个基础链路牢牢掌握后面无论做多远都有底气。我个人建议初学者第一次跑这套模拟时不要急着加各种花哨功能先把四种典型像差离焦、球差、彗差、像散分别跑一遍把相位图、PSF和MTF三张图对比着看。看多了对像差到底怎么影响成像的直觉就会变得非常准。这在光学设计里比任何公式都管用。