从Zernike波前到PSF/MTF：光学像质仿真完整指南

做光学系统设计的人应该都有过这种体验：镜头草图画完、光线追迹也跑通了，但别人问"那成像到底清不清楚、能分辨多小的细节"，你一时却拿不出一个直观的答复。光线追迹给的是光路，不是像质。真正用来回答"成像多锐利"的，是点扩散函数（PSF）和调制传递函数（MTF）。而要搞清楚像差是怎么把清晰的点像抹成模糊光斑的，最趁手的分析工具就是Zernike多项式——它能像拆解零件一样，把波前像差拆成离焦、像散、球差等一个个标准项，再通过标量衍射计算把波前误差映射成PSF和MTF曲线。这篇文章就是一套可以直接拿来做仿真的完整路径：从Zernike波前重建开始，到PSF计算，再到MTF提取，最后一个一个看常见像差对成像质量的影响，顺便把我踩过的采样和归一化的坑也一起讲了。

1. 从波前到光斑：PSF和MTF到底在刻画什么

1.1 波前像差：所有成像质量问题的最底层来源

先统一一下坐标系里的概念。一个完美成像系统，从物点发出的球面波经过镜头后，应该汇聚成一个以理想像点为中心的球面波。实际镜头做不到，出射波前相对理想球面波会有偏差，这个偏差就是波前像差。波前像差通常用光程差来表示，单位是纳米或微米，折算成相位就是 (2π/λ)×OPD。

为什么我强调要从波前而不是直接从光斑入手？因为波前是"因"，光斑和传递函数是"果"。同一种光斑形态可以对应完全不同的波前误差组合，但反推却不行。Zernike多项式的作用就是把这个连续分布的波前误差，展开成一组互相正交的模式。每个模式都对应一种经典像差，系数大小就是这种像差的强度。有了这套系数，整个系统的像差状态就变成了一个可以量化、可以传递、可以比较的列表。

1.2 PSF解决"一个点变成什么样"，MTF解决"细节还剩多少"

PSF的定义很直观：物平面上一个理想点光源，经过光学系统后，在像面上形成的强度分布。理想的衍射极限系统，PSF是艾里斑，中心亮斑加一圈圈暗环。有像差之后，能量不再集中在中心，可能偏到旁边、拉成条、或者散到外围，主瓣变宽，旁瓣变高。

MTF则把问题换了个问法：成像系统对不同空间频率的"响应能力"。它描述的是目标的正弦条纹经过系统后，调制对比度从1掉到了多少。MTF越高，说明这个频率的细节保留得越好；MTF掉到0的那个频率，就是系统的分辨率极限。

注意PSF和MTF不是两个独立的信息。MTF就是PSF的傅里叶变换的模，两者是同一枚硬币的两面。很多人在看仿真结果时只盯着PSF光斑图，这是个坏习惯。PSF好看不代表细节能传出去，低频高、中频塌掉的MTF曲线，用眼睛看PSF很难马上察觉，但曲线一眼就能看出来有问题。

1.3 从波前到MTF的完整计算链路

整条链路其实只有四步：

• 用Zernike系数在光瞳坐标系（单位圆）上重建波前像差 W(ρ,θ)
• 构造光瞳复振幅函数：光瞳内的振幅为1，相位是 (2π/λ)·W
• 对光瞳函数做二维傅里叶变换，取模平方得到PSF
• 对PSF再做一次傅里叶变换，取模并归一化，得到MTF

这里每一步都有物理前提：第三步用的是菲涅尔近似和夫琅禾费衍射，模拟的是焦平面上的光场分布；第四步用到的是非相干成像的OTF（光学传递函数），它的模就是MTF。整个链路在Zemax、Code V里也有对应的快速工具，但自己写一遍的好处是——你能知道每一步到底在算什么，系数变了之后曲线怎么动，你会比用黑盒工具清楚得多。

2. Zernike多项式：把任意像差拆成标准零件清单

2.1 为什么偏偏是Zernike

描述波前误差可以用很多基函数，比如幂级数、切比雪夫多项式，但只有Zernike是定义在单位圆域上、并且在圆域内满足正交性的。这带来两个非常实用的性质。

第一，各项之间线性独立。我把波前写成 Σc_j Z_j，每一项的系数 c_j 可以单独提取、单独分析，不会因为加了别的项而改变这一项的"身份"。检测出某个系数，就说明系统里存在对应程度的该种像差，可以直接指导装调动作。

第二，低阶项和光学设计里的经典初级像差一一对应。Zernike的离焦项和光学设计里的离焦是同一种波前形态；Zernike的像散项和赛德尔像散本质也是同一个物理量，只是归一化方式不同。这让你在Zernike系数和镜头设计直觉之间建立起一座桥，不至于对着数学符号发呆。

2.2 标准Zernike项的物理对应

在实际仿真里，我不需要写一个能算任意n、m阶的完整Zernike库——大部分像质评估关注到初级像差和少数高级项就足够了。下面这张表是标准Zernike项（也就是常说的Noll/标准归一化排列）前11项，系数这么写的时候，每一项对整个波前的RMS贡献正好等于系数值本身。

我强烈建议仿真时用标准Zernike排列而不是Zemax里的Fringe Zernike排列。因为标准排列做了RMS归一化，系数直接就是波前RMS误差（单位取λ时），不同项之间可以统一比较。Fringe排列的系数是峰谷归一化，离焦项和像散项的系数不能直接比，容易造成误判。

2.3 从系数到波前重建

重建波前的公式非常简单：

W(ρ,θ) = Σ c_j · Z_j(ρ,θ)

其中ρ是光瞳上的归一化半径，瞳孔边缘处ρ=1；θ是极角。仿真中生成ρ和θ的方式是先在直角网格上计算 Xd=2r/D−1（r是物理径向坐标），再转换到极坐标。这样整个光瞳就落在单位圆内，圆外区域在后面的瞳函数构造里直接置零。

这里有一个很多新手会犯的错误：Zernike多项式定义中的ρ必须是归一化半径，单位圆内最大为1。如果你直接用物理坐标代入，离焦项里的ρ²在边缘会是某个很大的数值，波前形状会完全错掉。所以第一步永远是先做好归一化坐标网格。

3. 从波前到PSF的仿真：标量衍射到底在算什么

3.1 瞳函数和夫琅禾费衍射

那步傅里叶变换的物理依据是：在夫琅禾费衍射条件下，焦平面上的复振幅分布正比于光瞳函数的二维傅里叶变换。具体来说：

E(x,y) ∝ F { P(ρ,θ) }

其中P是光瞳函数，写成P = A·exp(iφ)。A是光瞳振幅分布，通常均匀取1；φ是相位，由波前像差决定。这里有一个理解上的关键点：相位里的(2π/λ)·W是一个很大的数，W即便是0.1λ，相位也有0.63弧度，足以让PSF出现可观察的变化。所以这个傅里叶变换不是随便打个包求个模，相位分布才是决定PSF形态的灵魂。

3.2 一套可以直接跑的Python实现

下面这版代码是我在实际分析时反复用过的骨架。不追求库的完整性，注重的是把每个环节掰开看清楚。

import numpy as np import matplotlib.pyplot as plt # ---------- 系统参数 ---------- diameter = 4e-3 # 入瞳直径，4 mm wavelength = 550e-9 # 波长，550 nm focal = 50e-3 # 焦距，50 mm N = 512 # 采样点数 # ---------- 光瞳归一化坐标 ---------- lin = np.linspace(-1, 1, N, endpoint=False) X, Y = np.meshgrid(lin, lin) Rho = np.sqrt(X**2 + Y**2) Theta = np.arctan2(Y, X) mask = (Rho <= 1).astype(float) # ---------- 标准Zernike前几项 ---------- def zernike(j, Rho, Theta): if j == 1: return np.ones_like(Rho) * 1 elif j == 2: return 2 * Rho * np.sin(Theta) elif j == 3: return 2 * Rho * np.cos(Theta) elif j == 4: return np.sqrt(3) * (2 * Rho**2 - 1) elif j == 5: return np.sqrt(6) * Rho**2 * np.sin(2 * Theta) elif j == 6: return np.sqrt(6) * Rho**2 * np.cos(2 * Theta) elif j == 7: return np.sqrt(8) * (3 * Rho**3 - 2 * Rho) * np.sin(Theta) elif j == 8: return np.sqrt(8) * (3 * Rho**3 - 2 * Rho) * np.cos(Theta) elif j == 11: return np.sqrt(5) * (6 * Rho**4 - 6 * Rho**2 + 1) else: raise ValueError("本骨架只实现了常用项，你需要自行补全其他阶") def wavefront_from_coeffs(coeffs, Rho, Theta): # coeffs是一个字典，键为Zernike序号，值为以波长为单位的系数 W = np.zeros_like(Rho) for j, c in coeffs.items(): W = W + c * zernike(j, Rho, Theta) return W def calc_psf(coeffs, N, wavelength, diameter, focal): # 重建波前，单位：波长 W = wavefront_from_coeffs(coeffs, Rho, Theta) phase = 2 * np.pi * W # 光瞳函数 pupil = mask * np.exp(1j * phase) # 夫琅禾费衍射：傅里叶变换后取模平方 field = np.fft.fftshift(np.fft.fft2(np.fft.fftshift(pupil))) psf = np.abs(field)**2 psf = psf / psf.sum() # PSF平面坐标间隔 dx = wavelength * focal / diameter axis = np.fft.fftshift(np.fft.fftfreq(N, d=1.0)) # 适当处理，见下方说明 return psf, axis * dx * diameter / (wavelength * focal)

这个代码里有几个细节值得解释。

第一，np.fft.fftfreq(N, d=1.0)这里我是故意先不写真实的物理间隔。因为如果你想在PSF图上直接标出"多少个艾里斑半径"，最方便的是把PSF坐标归一化成无量纲参数 u = x·D/(λf)。在这个归一化坐标下，艾里斑第一暗环半径约等于1.22。我习惯直接把坐标轴换成这个无量纲单位来出图，物理含义一目了然。

第二，对pupil做fftshift再fft2，是因为坐标原点在数组中心。fftshift先把原点搬到DFT要求的(0,0)位置，变换后再搬回来，这样PSF的中心亮斑就位于图像正中间，方便后续提取径向曲线。

3.3 PSF归一化的一点经验

代码里我做了psf/psf.sum()，这是把PSF总能量归一化为1。这么做的好处是MTF计算里归一化会更干净，而且不同像差条件下对比PSF的旁瓣比例时，可以直接看强度的相对分布。

但有另一个细节容易踩：如果你直接在PSF图像里比较"峰值亮度"，归一化到总能量为1后，像差越大、光斑越散，峰值必然越低。这其实是合理的物理表现——能量守恒，摊薄了。可如果你拿这个归一化PSF去反推"系统透过率"，那就错了，因为这种归一化掩盖了真实的能量传递效率。要评估实际集光能力，应该做的是"中心艾里斑内的能量占比"这类指标，而不是简单看PSF峰值。

3.4 分辨率不够？先检查采样密度

采样数N的选择直接决定仿真的可信度。我一般的经验值是：模拟衍射极限系统时N=256勉强能用；一旦加入像差，尤其是球差这类高阶项，相位变化变得剧烈，建议N=512或更高。

判断标准很简单：检查最大相位梯度是否满足每个采样点上的相位变化小于π。换句话说，波前像差的表面如果出现"锯齿"，那FFT出来的PSF必然带着假的高频纹路，看起来像是真实的衍射环，实际是采样不足的伪迹。我做球差仿真的时候遇到过：Z11给到1λ时，N=256的结果里PSF边缘出现规则网格状图案，当时还以为是数值问题，后来把N提到1024就完全消失了——就是采样不够。

4. 从PSF到MTF：一条傅里叶变换的捷径

4.1 两种本质等价的计算方法

MTF有两条计算路径。

路径一是直接对PSF做傅里叶变换，取模，再除以零频值。路径二是利用光瞳函数的自相关：OTF的数学定义等于光瞳函数与其自身的自相关，再归一化。两条路径算出来的结果在离散实现里会有细微差别，但物理上等价。我倾向于路径一，因为从PSF出发更直观，而且你会发现，做MTF的过程其实就是又做了一次FFT，代码上省不了太多，但概念上特别干净。

def calc_mtf(psf, wavelength, diameter, focal): # 对PSF做FFT得到OTF otf = np.fft.fft2(np.fft.fftshift(psf)) # 峰值位于索引(0,0) mtf = np.abs(otf) mtf = mtf / mtf[0, 0] # 空间频率轴，单位：线对/毫米 dx = wavelength * focal / diameter freq_cycles_per_m = np.fft.fftshift(np.fft.fftfreq(N, d=dx)) freq_cycles_per_mm = freq_cycles_per_m / 1e3 # 归一化频率：用衍射极限截止频率归一 fc = diameter / (wavelength * focal) # 单位：cycles/m # 注意这是相干截止频率，非相干MTF截断也在同一尺度 norm_freq = freq_cycles_per_m / fc return mtf, freq_cycles_per_mm, norm_freq

这里有个容易混淆的点：空间频率轴我用了两个坐标。绘MTF曲线时，我一般把横轴画成归一化频率（cycles除以截止频率），这样所有系统都可以放在同一张图上比较；但在工程报告里，频率轴应该标成线对每毫米（lp/mm），因为传感器像素尺寸对应的就是物理分辨率。两种坐标之间的换算很简单：真实空间频率 = 归一化频率 × 截止频率。

4.2 MTF曲线的判读顺序

拿到MTF曲线后，我建议按照一个固定顺序来读，避免漏掉信息。

先看零频，MTF(0)必须等于1，这是归一化保证的。再看截止频率，曲线第一次掉到0的频率在哪，它告诉你系统的理论极限分辨率。然后看曲线形态：衍射极限的MTF是平滑单调下降的，像差会让曲线"提前塌陷"——可能在低频就掉得快，也可能在中频出现凹陷甚至零点。

最值得注意的是MTF曲线出现零点的情况。零点意味着该频率的调制度完全丢失，像面上这个频率的正弦条纹对比度是0。这在离焦过大的系统里很典型：中频出现一个零点，然后随着频率升高再有小幅回升，看起来像"蛇形"。你如果只看PSF，看到的只是"光斑变大了"，完全看不出这种周期性的信息丢失。

4.3 衍射极限MTF作为判据

判断像差容忍度的最快捷方法，是把当前系统的MTF和相同孔径/焦距下的衍射极限MTF叠在一张图上。衍射极限MTF有解析公式：

MTF(v) = (2/π) × [ arccos(v) − v·√(1−v²) ]

其中v是归一化空间频率，取值0到1。这个曲线我几乎每次仿真都会画一条作为基准线。比如在评价离焦0.2λ时，你会看到低频段可能还能维持90%以上的对比度，但到0.4归一化频率附近，已经掉到衍射极限曲线的下面一大截。把偏差最大的频率段标记出来，比盯着一个RMS波前误差数字有用得多。

5. 像差诊断实验：离焦、像散、球差在PSF和MTF上长什么样

5.1 离焦：能量摊平，中频塌陷

把Z4系数设为0.25λ（这是常见的光学公差预算量级），PSF中心峰值下降，第一暗环的环状结构开始模糊，能量向周围扩散。但真正戏剧性的变化在MTF上：低频段影响相对有限，中频段出现明显的MTF下降，0.3~0.6归一化频率之间的数值可能从接近0.8掉到0.4以下。继续加大离焦到0.5λ，MTF会在中频出现零点，这是不可修复的信息丢失，后期图像处理也很难救回来，因为该频率的调制度是物理性的零。

这里我要强调一个经验：很多人以为"所谓离焦就是整体变糊"，其实离焦模糊更像是在某个中间尺度上把对比度抽走了。小细节（高频）和大尺度轮廓（低频）都还在，就是中等尺度的纹理没了。所以判断离焦，看中频MTF非常敏感。

5.2 像散：子午和弧矢分道扬镳

像散的特殊之处在于它的轴对称性被打破了。Z5或Z6项会让PSF呈十字形的散焦特征——两个方向的焦点不重合，一个方向清晰的同时，另一个方向是模糊的。在PSF图上你会看到像散特有的椭圆或"双线"结构，方向取决于像散轴。

MTF上像散的标志姿态是：过中心横切和纵切的两条MTF曲线不重合。一条在中频还能维持较高对比度，另一条已经大幅塌陷。这个特征在诊断时极其有用——如果你测到一台镜头的MTF在水平方向和垂直方向差异明显，那基本可以断定这个系统有明显的非对称像差，像散或彗差是头号嫌疑。

5.3 球差：中心能聚焦，旁瓣拖着尾

Z11的球差有一个很反直觉的表现：光斑中心依然是个尖峰，但周围拖着一圈一圈亮度不低的旁瓣，能量并不集中在中心。这跟离焦"整体摊平"完全不一样。所以单看PSF主瓣大小，球差的影响容易被低估。

MTF对球差的暴露非常直接：低频部分尚可，但中频从0.2到0.5之间出现一个持续的凹陷，且不像离焦那样会出现明显的零点回落，而是整体"矮"了一截。球差还有一个特点，就是它的表现和焦面位置有关，你把像面在焦点前后移动，会发现某一侧PSF突然变得集中（即焦前/焦后不对称），这是区分球差和离焦的经典判据。

这些特征描述看起来简单，但做仿真的价值正在于此：你可以精确控制单一像差系数，观察它对PSF和MTF的独立贡献，建立"像差形态到曲线特征"的直觉。有了这套直觉，再去看实测的干涉仪数据和镜头测试报告，就不会一头雾水。

6. 这些坑我踩过：采样、归一化、边缘伪迹与仿真实战建议

6.1 采样不足的振铃是最大的假象来源

我在3.4里提过采样问题，这里展开说一个坑。孔径的圆形边缘在掩膜里是硬边界，直接用二值mask就会在FFT结果里引入振铃，表现为PSF外围有规则的同心圆环纹路。这个环纹和真实的衍射环长得非常像，但强度衰减模式和真实衍射不一样——真实艾里斑的旁瓣能量递减很有规律，振铃伪迹则会在顶部出现"平台状"的异常亮带。

缓解的办法有两层。第一层是提高采样密度，把瞳孔边缘的每个像素覆盖得更细；第二层是对掩膜边缘做亚像素平滑，比如用超采样：在一个像素内部判断圆边覆盖比例，把mask值设成0~1之间的小数。我的经验是，粗糙判断时用简单二值mask+高采样N已经够用；要做发表级结果，建议上超采样平滑，否则审稿人一眼就能看出边缘伪迹。

6.2 归一化约定：波前RMS和相位别混为一谈

这是新手最容易翻车的地方。Zernike系数的单位必须明确。我用系数 c_j 乘以 Z_j(ρ,θ)得到的 W，单位是波长。相位则是2πW。如果你看到别人的代码里相位直接就是c_j乘以Z_j，没有乘2π，说明他用的系数单位是弧度，不是波长。

还有一个隐蔽的问题：标准Zernike的正交性是指向"单位圆内面积分"的RMS归一化。系数0.25λ代表该项波前误差的RMS值为0.25λ。而系统总波前RMS是所有项系数平方和开根号（前提是各项正交且都在完整圆域内）。我经常看到有人把Z4和Z11的系数直接加起来当总像差，这是不对的——应该平方相加再开方。原因就是正交性，各项贡献的能量是叠加的，但振幅是不直接相加的。

6.3 FFT的周期性带来的边界伪迹

np.fft.fft2自带周期性边界条件。光瞳函数在单位圆外是0，但在FFT看来，数组的边缘会"卷绕"。当PSF的扩展范围很宽时（大像差），衍射场的能量跑到数组边缘又被卷回，形成对称的鬼影。这个问题的直接表现是：PSF图像四角出现不应该出现的亮点。

解决思路通常是"补零充填"：在计算FFT之前，把光瞳矩阵嵌入一个更大的全零矩阵。比如原矩阵N=512，外面套一圈零到N=1024再做FFT，相当于在原光瞳四周留出空白区，让衍射能量有地方扩散。这样得到的PSF主瓣更干净，代价是计算量变大。我在做Z11=1λ的球差仿真时，N=512直接算会出现四角伪迹，扩到N=1024并补零后就解决了。

6.4 从真实系统到Zernike系数的衔接

仿真做得再漂亮，系数总得有来源。实际项目中我拿到Zernike系数的途径主要有三种：干涉仪测量结果、光学设计软件的拟合输出、公差分析导出的敏感度列表。无论哪种，都要注意一件事：对方给你的编号约定是Fringe还是标准Zernike，极坐标方向定义是什么，有没有包含倾斜项。一个常见的尴尬是：Zemax里的Fringe Z4是离焦没问题，但Fringe Z11和标准Zernike的Z11球差系数符号可能相反，直接用会导致仿真结果和实测对不上。

我的建议是写一个转换辅助函数，输入原始系数列表、注明约定类型，输出统一到标准Zernike。这段代码不复杂，但能省下大量排查时间。

6.5 仿真实战检查清单

最后整理一份我在每次做波前到PSF/MTF仿真前都会过的检查表：

• 坐标系：ρ是否归一化到单位圆，θ的起点和方向是否与像差约定一致
• 单位链：Zernike系数是波长还是弧度，相位计算是否重复乘了2π
• 采样量：波前最大相位梯度下是否满足每像素相位变化小于π
• 边缘伪迹：光瞳是否补零，PSF四角是否有卷绕能量
• 归一化：PSF是否做了总能量归一化，MTF零频是否严格等于1
• 基准对照：是否叠画了理论衍射极限MTF曲线
• 交叉验证：换一个更大的N，确认PSF关键特征（旁瓣位置、MTF零点）不随N漂移

这套检查表是我吃过几次亏之后固化下来的。特别是最后一条，交叉验证N，花不了几十秒，但能帮你挡住绝大多数数值伪迹。做仿真最怕的不是算错，而是算出一个看起来无比合理、实际上是假的图案，还拿着它去指导下一次光学设计。习惯了先怀疑采样、再怀疑参数、最后才怀疑物理模型，你的仿真结果才能真正拿来作为决策依据。