☰
从Zernike系数到MTF曲线:光学像差模拟全链路与Python实现
2026/10/9 13:02:04 网站建设 项目流程

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系数定义了波前相位,波前相位决定了衍射积分的相位分布,衍射积分计算出PSF,PSF再做一次傅里叶变换得到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序把第一项叫piston,Fringe序里第一项却是piston加倾斜的组合。归一化系数也有区别:有的约定里,像散项的峰谷值就是系数的两倍,有的约定里系数的RMS值等于1。两种约定下,同一个物理像差写出来的系数甚至差出两倍多。

我写模拟时统一用标准Zernike(Born & Wolf版),即每项在单位圆上积分的RMS为1(除piston外)。这样做的好处是系数大小直接跟RMS误差挂钩,评价像质时方便换算。代码里我会直接按这个约定实现radial多项式,然后归一化每一项,再把系数乘上去。注意:如果后续要和Zemax或Code V对标,务必查清楚他们导出的Zernike系数用的是哪一套约定,否则结果没法直接比较。

2.2 光瞳网格采样:N=512够不够

FFT计算PSF的精度很大程度上取决于光瞳网格的采样密度。网格太稀,光瞳边缘的圆形边界会变成锯齿状,高频衍射分量被污染;网格太密,内存和计算时间又会飙升。

我常用的经验值是N=512。对于圆形光瞳,这个采样量足以保证PSF中心到第三四个暗环的形态都准确,MTF曲线也能平滑延伸到截止频率。如果是快速验证,N=256也能用,但MTF高频段会开始出现可见的噪声毛刺。N=1024适合做最终出图,特别是在要放大看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,圆外为0),W是波前像差函数(用波长单位)。然后对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接口更快,多核环境下自动并行,对于N=1024这种规模不会卡顿。

安装命令就不多写了,标准的三件套。要是你用的是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, axis=0) 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(N=512): """生成归一化坐标网格和圆形孔径掩膜""" 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, cmap='inferno', vmin=psf_log.max()-6, vmax=psf_log.max())

vmin取最大值往下6个数量级,基本能看到外圈衍射环的整个结构。如果你要量化中心能量,直接从归一化后的psf数组里取中心点的值即可,它其实就是斯特列尔比(Strehl Ratio)的近似。

相位缠绕不需要解包,因为exp(2jπ×0.1)和exp(2jπ×1.1)是完全一样的值。真正需要担心的反而是光瞳外的NaN——如果光瞳掩膜直接把r>1置零,而你在phase_masked里又对r=0做了除零操作,结果会出现nan,沿着孔径边缘扩散。代码里已经规避了这个:aperture乘在最后。

3.4 第三步:提取MTF曲线

MTF的算法已经说过,重点在于提取曲线和坐标换算。先算出二维MTF:

otf = 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, d=2.0/N)) profile_x = mtf[N//2, :] profile_y = mtf[:, N//2] plt.plot(freq, profile_x, label='tangential') plt.plot(freq, profile_y, label='sagittal') 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项去掉(设为0),PSF的中心位置不会移动,但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,光瞳外的"相位噪声"也会参与FFT,PSF就会出现莫名其妙的背景能量。每次构造光瞳函数前,我都建议先把光瞳图和相位图并排打印出来看一眼,确认掩膜对不对齐。

4.2 MTF高频掉太快还是全是毛刺

MTF曲线出现毛刺,大概率是PSF采样不足导致的。PSF在中心附近的峰值非常陡峭,如果网格不够密,FFT结果的高频分量就会出现频谱泄漏。解决办法是提高N,或者在做FFT之前对光瞳数组做零填充(pad)到4倍大小。零填充不增加信息量,但能让频域插值更平滑,曲线上的毛刺会少很多。

如果MTF高频掉得比理论值快很多,可能是PSF没有正确归一化,或者PSF数组里有零值导致log运算出现inf。我调试时会先单独打印psf.min()和psf.max(),如果min=0,取对数时记得加一个小常数1e-9。MTF归一化也有讲究:要用mtf[0,0](零频分量)作为除数,而不是mtf.max(),因为数值误差下两者可能略有差异。我一般直接用中心值,这样曲线严格从1出发。

4.3 相位缠绕和NaN问题

相位本身不会缠绕出问题,因为复数指数是周期函数。真正麻烦的是在极坐标网格里,当r=0时,theta计算的是arctan(0,0),在NumPy里这个值没问题,但Zernike多项式的cos(m·theta)项会变成cos(0)或者cos(π/2),导致中心点数值不确定。

处理办法:要么单独给r=0点赋一个固定值(一般取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三张图对比着看。看多了,对"像差到底怎么影响成像"的直觉就会变得非常准。这在光学设计里比任何公式都管用。

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

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

立即咨询