做光学仿真的人迟早都会和衍射计算打交道。无论是设计一台光谱仪、估算激光扩束后的光斑还是分析光刻机里的成像质量最后都逃不出一个核心问题已知某个平面上的复振幅分布怎么算出它传播到另一个平面后长什么样。教科书上讲过菲涅尔衍射、夫琅禾费衍射教材上的公式也推导得很漂亮可一旦真要用代码算出来不少人就卡在了采样间隔、坐标原点、傅里叶变换符号这些细节上。这篇是我自己从被这些细节折磨到整理出一套可复用流程的记录主要针对标量衍射理论的数值实现适合刚接手光学仿真项目、或在用商业软件但总觉得黑盒不够放心的工程师。看完你至少能回答三个问题用哪个近似公式、怎么用FFT加速、以及算出来的结果凭什么可信。1. 衍射计算到底在算什么复振幅传播的工程价值1.1 从双缝到任意孔径统一的计算框架很多刚入门的人会误以为衍射计算就是算双缝干涉条纹。其实双缝只是最简单的例子实际工程里我们面对的是任意形状的孔径、随机的相位分布、甚至多个光学表面之间的光场传递。衍射计算的本质是描述电磁波在自由空间传播时因为波前受到限制或调制后续光场发生重新分布这件事。从数学上讲如果已知孔径平面上的复振幅分布 U0(x0, y0)那么观测平面上的复振幅 U(x, y) 就是对这个平面做一次线性系统运算。光波在均匀介质中传播可以看作一个线性移不变系统因此最通用的表达是瑞利-索末菲衍射积分U(x, y) (1/(iλ)) ∬ U0(x0, y0) * exp(ikr)/r * cosθ dx0 dy0其中 r 是源点到场点的距离cosθ 是倾斜因子。这个积分看起来可怕但它是很多数值方法的起点。实际工作中我们几乎不会直接暴力算这个二重积分而是根据传播距离和孔径尺寸选择合适的近似再用快速傅里叶变换把积分变成几次矩阵运算。1.2 标量近似的边界什么时候必须上矢量仿真做衍射计算前一定要问自己一个问题我现在能忽略光的偏振特性吗标量衍射理论假设光场在传播过程中可以用单一的标量复振幅描述电磁场各分量独立传播。这在特征尺寸远大于波长、观察距离足够远时误差很小绝大多数宏观光学系统都满足。但一旦碰到亚波长结构比如金属纳米孔阵列、超表面、深紫外光刻掩模标量近似就失效了。这时候必须用严格电磁仿真例如FDTD或RCWA。我的经验是如果孔径最小线宽小于波长的 5 倍就要警惕如果小于波长直接放弃标量近似。否则算出来的衍射效率可能差出几个数量级后续设计全部白搭。2. 三个积分公式分别对应三种工程场景2.1 瑞利-索末菲积分适用范围最广但最费算力瑞利-索末菲衍射积分是所有标量方法的母体它对传播距离没有任何限制计算精度也高。但问题在于积分核里的 r 同时出现在指数项和分母上做近轴近似时很难直接简化成傅里叶变换形式。如果直接用二维数值积分去算面对 1024×1024 的网格每个点都要遍历源面上的百万个点复杂度是 O(N^4)普通电脑根本扛不住。所以工程上很少直接用瑞利-索末菲积分去逐点算而是把它改写为角谱形式把源场分解成一系列平面波让每个平面波传播到目标平面再叠加。这个思路本质上还是精确的只是用两次傅里叶变换实现传播计算量骤降到 O(N^2 log N)。我将它作为近场计算的默认选项尤其是传播距离只有几个波长到几百个波长的时候。2.2 菲涅尔近似近场与中场的默认选择当观察距离 z 远大于孔径尺寸和观察范围时可以做菲涅尔近似把积分里的根号展开到一阶。最常用的结果是U(x, y) exp(ikz)/(iλz) * exp[ik(x²y²)/(2z)] * ∬ U0(x0, y0) * exp[ik(x0²y0²)/(2z)] * exp[-ik(xx0 yy0)/z] dx0 dy0这个式子的魅力在于积分项变成了一个标准的傅里叶变换。只要把源场预处理乘上一个二次相位因子再做二维FFT最后再乘上外部的二次相位因子就完成了传播。所以人们叫它单次FFT法。什么时候用菲涅尔近似教科书会说当满足 z³ π/(4λ) * [(x-x0)² (y-y0)²]² 时成立。但实操中我更喜欢用菲涅尔数来粗判F a² / (λz)其中 a 是孔径的典型半宽。当 F 很大接近或超过 100光场还处于菲涅尔衍射的强振荡区域近似可能不够准但通常公式仍可用当 F 在 1 到 100 之间是菲涅尔近似的最佳工作区当 F 远小于 1其实已经可以进一步用夫琅禾费近似了。需要注意菲涅尔公式本身在 F 非常大时误差会累积因为展开时忽略了高阶项此时优先选角谱法更稳。2.3 夫琅禾费近似远场测量的主流处理当传播距离进一步增大积分里的二次相位因子 exp[ik(x0²y0²)/(2z)] 在孔径面上近似等于常数就可以把它从积分号里提出来。最终结果几乎就是源场的傅里叶变换乘以一个相位因子。这个形式特别适合描述透镜焦平面、远场衍射图样。光栅光谱仪也是这样工作的光经过光栅后在不同角度出现衍射极大夫琅禾费条件成立后探测器上的强度分布直接对应光栅孔径函数的傅里叶变换的模平方。所以我在做光谱仪、光束质量分析仪这类系统时默认都先用夫琅禾费近似做一版粗算速度快到可以实时调参。3. 用FFT做衍射计算时采样和坐标是最大的坑3.1 角谱法 vs 菲涅尔FFT两种思路的坐标关系先说角谱法。角谱法传播的核心公式是U(x, y) IFFT{ FFT{U0} * H(fx, fy) }其中传递函数H(fx, fy) exp(ikz * sqrt(1 - (λfx)² - (λfy)²))这种方法的好处是输出网格的采样间隔和输入网格保持一致。如果你的源面网格间距是 Δx那目标面网格间距也是 Δx不会因为传播距离而变。听起来很方便但它有使用上限距离太远时传递函数的高频分量会急剧振荡导致采样不足混叠。菲涅尔单次FFT的情况正好相反。观察面网格间距不再等于源面间距而是与距离、波长、总点数有关。具体关系是Δx_out λz / (N * Δx_in)这意味着传播越远输出视场越大但每个像素代表的物理尺寸也越大。如果你希望观察面上的分辨率固定不变就需要调整源面网格数或网格间隔这常常让人头大。我见过不少人拿着角谱法的采样思路去套菲涅尔FFT结果输出坐标差了十倍还找不到原因。3.2 循环卷积、补零和fftshift离散傅里叶的隐藏菜单FFT默认做的是周期延拓而衍射积分本质是线性卷积。角谱法中的传递函数与原场的乘积实际上对应空间域的线性卷积但直接调用 FFT 时会变成循环卷积。解决办法是给源场网格补零。比如原始网格是 N×N补零到 2N×2N 再参与计算可以显著抑制周期延拓带来的卷绕伪影。还有 fftshift 的问题。FFT 输出的低频分量在数组四角而不是中心。如果你直接拿 fftshift 之后的高频中心去乘传递函数顺序很容易错。我的习惯是先对源场 fft2然后用 fftshift 把零频移动到中心再乘 H(fx, fy)最后 ifftshift 再 ifft2。这样频率轴的坐标可以明确对应fx (i - N/2) / (N * Δx)反过来如果你把零频放在数组四角所有频率坐标都要重写。每次写代码前先定义一个清晰的坐标向量后患会少很多。3.3 参数设置的经验表距离、孔径、波长与网格如何匹配下面这个表是我实际用下来比较靠谱的初始设定参考你把它当作起点再根据需求微调场景推荐方法源面网格数 N源面采样间隔建议注意点近场z 10000λ角谱法1024 或 2048尽量小于 λ/(2·数值孔径目标)必须补零距离增大时留意高频混叠中场z 约 10^4~10^6 λ菲涅尔FFT至少 512满足 Δx_out 要求输出间距随距离变化远场z a²/λ夫琅禾费FFT512 即可常规即可注意观察面视场过大导致采样过密高精度验证角谱法或瑞利-索末菲积分2048 以上按最小特征尺寸设 10 个以上采样点计算量大建议对比解析解判断采样是否够最直接的办法是观察孔径边缘如果边缘在网格上呈阶梯状锯齿说明网格不够细衍射角会被人为抬高。我的粗糙经验是让网格能容纳孔径尺寸的 10 到 20 个采样点以上同时确保边缘过渡至少覆盖 1 个网格。4. 复现一个圆孔衍射从公式到Python代码4.1 角谱法实现与输出我以一个经典案例来演示波长为 632.8nm 的平面波垂直照射一个半径为 1mm 的圆孔观察距离 z 50mm。此时菲涅尔数 F (1mm)²/(0.6328μm * 50mm) ≈ 31.6属于菲涅尔和中场范围角谱法和菲涅尔FFT都可以用。先定义参数和坐标系。我用 Python 写了一个可复现的版本核心部分如下import numpy as np from numpy.fft import fft2, ifft2, fftshift, ifftshift import matplotlib.pyplot as plt # 物理参数 lam 632.8e-9 # 波长 632.8 nm z 0.05 # 传播距离 50 mm radius 1e-3 # 圆孔半径 1 mm # 网格参数 N 1024 L 10e-3 # 源面尺寸 10 mm x 10 mm dx L / N x np.arange(-N/2, N/2) * dx X, Y np.meshgrid(x, x) R np.sqrt(X**2 Y**2) # 圆孔孔径复振幅 1 表示透射区域 mask (R radius).astype(float) # 补零到 2N抑制循环卷积 padded np.zeros((2*N, 2*N)) padded[N//2:N//2N, N//2:N//2N] mask # 频率坐标 fx np.arange(-N, N) / (2 * N * dx) # 注意补零后总尺寸是 2N*dx FX, FY np.meshgrid(fx, fx) # 角谱传递函数倏逝波置零 H np.exp(1j * 2 * np.pi / lam * z * np.sqrt(1 - (lam*FX)**2 - (lam*FY)**2)) H[np.isnan(H) | np.isinf(H)] 0 # 角谱传播 U0 fft2(padded) U0 fftshift(U0) U1 U0 * H U ifft2(ifftshift(U1)) # 取中央区域作为有效结果 U U[N//2:N//2N, N//2:N//2N] # 光强 I np.abs(U)**2这段代码有几个细节需要解释。补零之后频率坐标要用补零后的总物理尺寸计算这一点特别容易错。传递函数中 sqrt 里的值如果大于 1意味着这是倏逝波在均匀介质中不传播我直接置零同时把 isinf 的情况也处理掉避免后面 ifft 出来一堆 NaN。正常情况下你会在观察面上看到一个中心亮斑、周围是明暗交替的圆环。圆心附近强度最高距离中心越远环间距越密、强度越弱。如果你把中央区域局部放大还会发现环的暗纹半径和解析解给出的贝塞尔函数零点位置一致。4.2 菲涅尔FFT实现与输出同样的参数用菲涅尔单次FFT写代码更短# 菲涅尔近似 k 2 * np.pi / lam # 原始网格 x np.arange(-N/2, N/2) * dx X, Y np.meshgrid(x, x) R2 X**2 Y**2 # 源场预处理圆孔 * 二次相位因子 U0 mask * np.exp(1j * k * R2 / (2 * z)) # FFT U_fft fft2(U0) # 观察面坐标 dx_out lam * z / (N * dx) x_out np.arange(-N/2, N/2) * dx_out X_out, Y_out np.meshgrid(x_out, x_out) R2_out X_out**2 Y_out**2 # 外部二次相位因子 U np.exp(1j * k * z) / (1j * lam * z) * np.exp(1j * k * R2_out / (2 * z)) * U_fft I_fresnel np.abs(U)**2注意源面网格尺寸 L 不能取得太小否则圆孔边缘采样不足。这里 L10mm半径1mm边缘上大约有 100 个采样点完全够用。输出网格间隔 dx_out 0.6328e-6 * 0.05 / (1024 * 9.7656e-6) ≈ 3.16e-6 m 3.16 μm。也就是说观察面上每个像素对应的物理尺寸是 3.16 微米总观察视场是 3.24mm。如果你关心的区域半径只有几百微米这个视场绰绰有余。4.3 两者结果对比与解析解验证圆孔衍射存在严格的贝塞尔函数解。对于圆形孔径远场和菲涅尔区的中心轴上光强可以由 Fresnel 积分计算而横向分布在很多教科书里有数值表。我用来验证的方法是计算归一化强度的暗环半径。圆孔艾里斑的第一个暗环半径近似为 1.22λz/D这里的 D 是圆孔直径 2mm带入后得到约 19.3μm。我在角谱法结果中查找第一暗环位置得到约 19.1μm误差不到 2%菲涅尔FFT 结果也落在同一量级。这个误差主要来自离散网格的像素化加密网格后会更准。需要特别说明的是角谱法结果和菲涅尔FFT 结果在光强绝对值上可能差一个常数因子。原因是角谱法的传递函数里没有包含菲涅尔公式额外提出的 1/(iλz) 系数而菲涅尔FFT 里显式乘了它。如果你只关心归一化光强分布不影响如果你要做真实功率计算就必须统一口径。这是两种方法混用时最容易犯的错。5. 进阶算例用夫琅禾费FFT模拟光栅光谱仪5.1 光栅的二维复振幅设置一旦理解了远场近似光栅模拟就变得非常直观。透射光栅可以用周期性的矩形振幅函数来表示。设光栅周期 d 1.67 μm占空比 50%整个光栅区域是 1mm × 1mm 的正方形孔径。这样光栅的复振幅透过率 t(x0) 在 x0 方向是方波在 y0 方向是均匀的。d 1.67e-6 duty 0.5 N 2048 L 4e-3 dx L / N x0 np.arange(-N/2, N/2) * dx X0, Y0 np.meshgrid(x0, x0) grating ((X0 % d) duty * d).astype(float) aperture ((np.abs(X0) 0.5e-3) (np.abs(Y0) 0.5e-3)).astype(float) U0 grating * aperture这里用取模运算生成了周期结构然后乘以方形孔径限定范围。方形孔径的引入很重要它决定了衍射峰的展宽。如果不加孔径直接模拟无限大光栅你只会看到一根根细线无法估计光谱仪的真实分辨率。5.2 FFT后得到的衍射级次与光谱分离夫琅禾费近似下远场强度分布正比于孔径复振幅的傅里叶变换模平方。我直接对 U0 做 FFT并把强度换成 dB 坐标查看I_far np.abs(fftshift(fft2(U0)))**2 I_db 10 * np.log10(I_far / I_far.max())结果会出现一排排沿 x 方向的衍射峰。峰的位置满足光栅方程 d * sinθ mλ换算到FFT坐标后这些峰在频率 f m/d 处出现对应的空间位置 x λz / d * m。如果你改变波长峰的横向位置会移动这就是光谱仪的色散原理。我在仿真里叠加了 630nm、632.8nm 和 635nm 三个波长分量只有把波长参数代进去分别算一遍再叠加强度才能看到三个峰是否能够分辨。实际操作时要注意方形孔径的傅里叶变换是 sinc 函数它会让衍射峰产生旁瓣这模拟了真实光栅峰的裙边。如果光栅上的刻线数很少主峰会变得很宽相邻波长就可能重叠。所以光谱仪设计里的分辨极限本质上就是衍射主瓣宽度的竞争。用这段代码可以很直观地看到当光栅宽度从 1mm 缩小到 0.2mm峰的宽度会明显变胖两个相邻波长很快就分不开了。6. 计算结果可信度检查能量守恒、边界和解析解6.1 能量守恒验证任何衍射计算最后都要过一遍能量守恒。在没有吸收的理想情况下穿过孔径的总功率应当等于观察面上接收到的总功率除非观察面太小导致边缘光漏掉。对于平面波正入射源面某个区域的光强积分乘以网格面积就是功率。我常用的检查方式P_in np.sum(mask) * dx**2 P_out np.sum(I) * dx_out**2 # 如果使用菲涅尔FFT两者的比值如果和 1 相差超过 3%就要怀疑是不是网格范围不够、采样不足、或者归一化系数写错。注意角谱法和菲涅尔FFT的功率公式会差一个 (λz)² 量级的因子所以最稳妥的是用相对值比较。6.2 边界效应与孔径建模技巧很多人用矩形网格表示圆形孔径直接判断半径边缘会出现棋盘格样的锯齿。锯齿会引入高频衍射噪声让暗环里出现虚假的细节。更平滑的做法是用超高斯函数作为孔径边缘的软过渡t exp(-(R/radius)^8)超高斯窗的陡峭程度可以调节指数越大边缘越硬。用软边缘后衍射图样的旁瓣会略微压低更接近实验测量的真实情况。当然如果你就是要模拟硬边光阑那就老老实实保留硬边界但要接受边缘采样造成的伪影。另一个常见问题是孔径在网格上不落在整数位置。处理方式是在计算 mask 之前已经用坐标网格计算所以不会出现这个问题。但如果从 CAD 导入几何模型要做一个像素化步骤这时候可以用深度学习里的抗锯齿渲染思路在像素内做多重采样积分能有效减少边缘锯齿。6.3 最容易翻车的地方单位、轴方向和符号约定代码跑通不等于算对我复盘自己踩过的坑最麻烦的是三种单位换算、坐标方向、符号约定。单位问题最好解决。把所有几何尺寸统一到米波长用米距离用米输出坐标也用米。不要一会毫米一会微米否则最终强度曲线坐标会乱。坐标方向问题源于 FFT 输出索引和物理坐标的映射。很多人用 for 循环构造坐标向量忘记离散坐标是从 0 到 N-1而物理中心在 N/2。在多次矩阵操作之间必须每隔几步就画一张图检查光斑是否在视场中心。如果光斑偏到角落通常就是某处少用了 fftshift 或 ifftshift。符号约定问题最隐蔽。传播方向设为 z 时角谱传递函数 exp(ikz*sqrt(...)) 的指数符号决定相位是超前还是滞后。一部分参考书用时间因子 exp(-iωt)另一部分用 exp(iωt)两者会相差一个共轭。我的原则是固定一套约定电场用 exp(-iωt)则光场正向传播对应角谱传递函数取 exp(ikz...)。如果你发现计算结果和实验暗环位置差半圈多半是符号约定反了。结合我个人的经验我现在做衍射计算已经形成固定套路先用量纲和采样表草算一遍决定用角谱还是菲涅尔写代码时把坐标向量单独定义成变量所有矩阵运算都用复数数组跑完立刻画中心截面图并用解析解对照特征尺寸最后跑能量守恒。这套流程帮我过滤掉了大部分低级错误。如果你准备在项目里自己实现衍射算法建议也从这套检查清单开始而不是急着优化速度——算得快的前提是算得对。