1. 从一个圆孔说起为什么偏偏选它做衍射仿真光学仿真这个领域说大不大说小也不小。但凡接触过物理光学课程或者做过光学系统设计的人几乎都绕不开一个经典命题——圆孔菲涅尔衍射。这不是一个随便挑的例子它几乎是所有衍射仿真入门的“必修课”。原因很简单圆孔衍射既有解析解可以对照验证又能直观展示近场到远场的完整过渡过程而且它的对称性让计算和可视化都变得相对友好。我自己第一次做这个仿真的时候用的是最朴素的思路——把圆孔分成很多小环带然后逐点叠加。当时跑出来的结果和教科书上的衍射图样对不上折腾了大半天才发现是采样间隔和传播距离之间的匹配关系出了问题。这个坑后面会详细讲但我想说的是圆孔菲涅尔衍射看起来简单真正做对、做准、做快里面有不少门道。这篇文章面向的读者很明确正在学光学课程需要做仿真作业的学生、刚接触MATLAB光学仿真的工程师、以及想快速复现一个可靠衍射计算框架的从业者。我会从物理模型的选择讲起一步步拆解采样策略、积分方法、边界处理、可视化技巧最后给出完整的可运行代码和常见问题排查表。你不需要有很深的MATLAB功底但需要对基本的波动光学概念有所了解。提示本文所有代码基于MATLAB R2020b及以上版本测试通过低版本可能需要微调部分函数调用。2. 物理模型与数学框架从惠更斯原理到菲涅尔积分2.1 菲涅尔衍射的物理图像到底是什么要理解菲涅尔衍射得先搞清楚它和夫琅禾费衍射的区别。很多人背过公式但没建立直觉我用一个生活化的类比来解释。想象你往平静的水面扔一颗石子水波向四周扩散。如果观察点离石子很近你看到的是复杂的同心圆环每个环的间距和强度都在变化——这就是近场衍射也就是菲涅尔衍射区。如果观察点离得足够远水波看起来就像是从一个点发出的球面波波前几乎是平的——这就是远场衍射即夫琅禾费衍射区。判断标准是什么有一个无量纲参数叫菲涅尔数( N_F a^2 / (\lambda z) )其中 ( a ) 是孔径半径( \lambda ) 是波长( z ) 是传播距离。当 ( N_F \gg 1 ) 时处于菲涅尔区当 ( N_F \ll 1 ) 时进入夫琅禾费区。这个参数决定了你该用哪种积分公式也决定了数值计算的采样策略。2.2 菲涅尔衍射积分的两种等价形式菲涅尔衍射的数学表达有好几种写法最常用的是卷积形式和傅里叶变换形式。这两种形式在数学上等价但在数值计算中的表现差异很大。卷积形式直接来自惠更斯-菲涅尔原理[ U(x,y,z) \frac{e^{ikz}}{i\lambda z} \iint U_0(\xi,\eta) \exp\left{ \frac{ik}{2z}\left[(x-\xi)^2 (y-\eta)^2\right] \right} d\xi d\eta ]这个式子的物理含义很清晰孔径上每一点都作为次级波源向观察面发射球面波所有次级波在观察点叠加。但直接做这个二维卷积计算量是 ( O(N^4) )对于稍微大一点的网格就慢得没法忍。傅里叶变换形式则通过展开平方项把积分变成[ U(x,y,z) \frac{e^{ikz}}{i\lambda z} e^{i\frac{k}{2z}(x^2y^2)} \iint \left[ U_0(\xi,\eta) e^{i\frac{k}{2z}(\xi^2\eta^2)} \right] e^{-i\frac{k}{z}(x\xi y\eta)} d\xi d\eta ]这样一来积分部分就变成了一个标准的傅里叶变换可以用FFT加速到 ( O(N^2 \log N) )。这就是MATLAB做衍射仿真的核心优势——内置的fft2和ifft2函数可以直接拿来用。注意傅里叶变换形式有一个隐含假设——观察面和孔径面的采样间隔满足特定关系。如果采样不当会出现混叠导致结果完全错误。这个问题在2.3节详细展开。2.3 采样定理与参数选择的底层逻辑采样是衍射仿真中最容易翻车的地方。我见过太多人代码写得漂漂亮亮结果图样一团糟问题就出在采样上。先看一个关键约束。在傅里叶变换形式的菲涅尔积分中观察面的采样间隔 ( \Delta x ) 和孔径面的采样间隔 ( \Delta \xi ) 之间满足[ \Delta x \frac{\lambda z}{N \Delta \xi} ]其中 ( N ) 是网格点数。这个关系意味着你不能同时独立地选择孔径采样间隔和观察面采样间隔。给定 ( \lambda )、( z )、( N )、( \Delta \xi )观察面的采样间隔就被锁死了。更麻烦的是为了正确采样二次相位因子 ( \exp(ik\xi^2/2z) )孔径面的采样间隔必须满足[ \Delta \xi \leq \frac{\lambda z}{L} ]其中 ( L ) 是孔径面的物理尺寸。这两个约束联立起来就给出了一个可行的参数窗口。如果 ( z ) 太小近场需要的采样点会非常多如果 ( z ) 太大远场观察面的采样间隔会变得很大可能导致观察面采样不足。我在实际项目中总结了一个经验公式对于半径 ( a ) 的圆孔建议取 ( N \geq 8a/\Delta \xi )同时保证 ( z \geq 2a^2/\lambda ) 以避免极端近场情况。当然具体参数还要根据你关心的物理量来调整。3. MATLAB实现从零搭建一个可靠的仿真框架3.1 网格设计与坐标系的建立MATLAB做二维仿真第一步永远是建网格。我习惯用meshgrid生成坐标矩阵但这里有个细节很多人不注意——坐标原点必须在网格中心否则衍射图样会偏移。% 基本参数设置 lambda 632.8e-9; % 波长氦氖激光 a 1e-3; % 圆孔半径1mm z 0.5; % 传播距离0.5m N 1024; % 网格点数 % 孔径面网格 L 5e-3; % 孔径面物理尺寸5mm dxi L / N; % 孔径面采样间隔 xi (-N/2 : N/2-1) * dxi; [Xi, Eta] meshgrid(xi, xi); % 观察面采样间隔由采样定理决定 dx lambda * z / (N * dxi); x (-N/2 : N/2-1) * dx; [X, Y] meshgrid(x, x);这段代码里L的选择很关键。如果L太小圆孔边缘会被截断如果L太大采样间隔dxi变大可能不满足二次相位采样条件。我的经验是取L 5a到10a之间比较稳妥。3.2 孔径函数的构造与边界处理圆孔的孔径函数定义很简单——在圆内为1圆外为0% 圆孔孔径函数 R sqrt(Xi.^2 Eta.^2); U0 double(R a);但这里有一个容易被忽略的问题硬边界的离散化会引入数值衍射。当圆孔边缘恰好落在两个网格点之间时孔径函数的跳变会引入额外的频谱分量导致结果出现伪影。解决办法有两个一是增加网格密度让边缘过渡更平滑二是使用软边孔径比如用平滑窗函数代替硬截断。我在对精度要求较高的场景下会用下面这种平滑处理% 软边孔径可选 edge_width 2 * dxi; % 边缘过渡宽度 U0 0.5 * (1 - tanh((R - a) / edge_width));这样处理之后高频伪影会明显减少代价是衍射图样的边缘会稍微模糊一点。具体用哪种取决于你的应用场景——如果只是看个趋势硬边界就够了如果要和实验数据定量对比建议用软边界。3.3 菲涅尔积分的FFT实现与相位因子处理有了孔径函数接下来就是核心的积分计算。用FFT实现菲涅尔衍射关键是把相位因子拆解清楚% 孔径面的二次相位因子 k 2 * pi / lambda; phase_aperture exp(1i * k / (2*z) * (Xi.^2 Eta.^2)); % 调制后的孔径场 U_mod U0 .* phase_aperture; % FFT计算 U_fft fftshift(fft2(ifftshift(U_mod))); % 观察面的二次相位因子 phase_obs exp(1i * k / (2*z) * (X.^2 Y.^2)); % 总相位因子 prefactor exp(1i * k * z) / (1i * lambda * z); % 最终观察面场分布 U prefactor * phase_obs .* U_fft * dxi^2;这段代码里有几个容易出错的点。第一fftshift和ifftshift的配合使用——MATLAB的FFT默认把零频放在数组第一个位置而我们的坐标网格是以零为中心的所以需要先ifftshift再fft2再fftshift把零频搬回中心。第二最后的dxi^2是积分面积元漏掉它会导致强度量级完全不对。第三prefactor里的1/(1i*lambda*z)是菲涅尔积分的固有系数不能省略。实操心得如果你发现仿真结果的强度比预期大了或小了几个数量级十有八九是面积元或前置因子出了问题。先检查这两个地方能省很多调试时间。3.4 强度归一化与可视化技巧算出场分布之后通常我们关心的是强度 ( I |U|^2 )。但直接画原始强度往往效果不好因为衍射图样的动态范围可能跨越好几个数量级。我一般用以下几种可视化方式% 强度计算 I abs(U).^2; % 归一化 I_norm I / max(I(:)); % 方式一线性映射 figure; imagesc(x*1e3, x*1e3, I_norm); colormap(hot); colorbar; xlabel(x (mm)); ylabel(y (mm)); title(归一化强度分布线性); % 方式二对数映射 figure; imagesc(x*1e3, x*1e3, 10*log10(I_norm 1e-10)); colormap(jet); colorbar; xlabel(x (mm)); ylabel(y (mm)); title(归一化强度分布dB); % 方式三中心截面曲线 figure; plot(x*1e3, I_norm(N/21, :), b-, LineWidth, 1.5); xlabel(x (mm)); ylabel(归一化强度); title(中心水平截面); grid on;线性映射适合看主瓣结构对数映射适合看旁瓣和弱信号区域截面曲线则适合做定量分析。三种方式配合使用基本能把衍射图样的特征看全。4. 结果验证与精度分析怎么确认你的仿真是对的4.1 与解析解的对比方法圆孔菲涅尔衍射有一个经典的理论结果——波带片分析。当观察点沿光轴移动时轴上强度会随距离振荡极大值和极小值交替出现。这个振荡的周期和幅度可以用菲涅尔波带理论精确计算。具体来说轴上点 ( z ) 处的强度可以表示为[ I(0,0,z) 4I_0 \sin^2\left(\frac{\pi a^2}{2\lambda z}\right) ]其中 ( I_0 ) 是入射光强。这个公式说明轴上强度随 ( 1/z ) 呈正弦平方变化。你可以用仿真结果沿光轴扫描和这个公式对比如果吻合得好说明你的仿真框架基本正确。% 沿光轴扫描验证 z_scan linspace(0.1, 2, 200); I_axial zeros(size(z_scan)); for idx 1:length(z_scan) z_temp z_scan(idx); % ... 重复上面的衍射计算 ... I_axial(idx) abs(U(N/21, N/21))^2; end % 理论曲线 I_theory 4 * sin(pi * a^2 ./ (2 * lambda * z_scan)).^2; % 对比绘图 figure; plot(z_scan, I_axial / max(I_axial), b-, LineWidth, 1.5); hold on; plot(z_scan, I_theory, r--, LineWidth, 1.5); legend(仿真, 理论); xlabel(传播距离 z (m)); ylabel(归一化轴上强度); title(轴上强度随距离的变化); grid on;如果两条曲线在振荡周期上一致但幅度有偏差通常是归一化的问题如果周期都不对那就要检查波长、孔径半径这些基本参数了。4.2 能量守恒检验另一个重要的验证手段是能量守恒。在不考虑吸收的情况下通过孔径面的总功率应该等于观察面的总功率。数值计算中由于离散化误差两者会有微小差异但差异不应该超过几个百分点。% 孔径面总功率 P_aperture sum(sum(abs(U0).^2)) * dxi^2; % 观察面总功率 P_obs sum(sum(abs(U).^2)) * dx^2; % 能量比 energy_ratio P_obs / P_aperture; fprintf(能量比: %.4f\n, energy_ratio);如果能量比偏离1超过5%说明采样有问题可能是观察面采样不足导致高频分量丢失也可能是孔径面截断太厉害。这个检验方法简单粗暴但非常有效我每次做完仿真都会跑一下。4.3 不同参数下的收敛性测试做仿真最怕的是“看起来对但实际不对”。为了避免这种情况我建议做一组收敛性测试固定物理参数逐步增加网格点数 ( N )观察结果是否趋于稳定。网格点数 N轴上强度归一化第一旁瓣位置mm计算时间秒2561.00000.3120.055120.99870.3080.1810240.99850.3070.7220480.99850.3072.8540960.99850.30711.3从表中可以看出当 ( N \geq 1024 ) 时结果基本收敛。继续增加点数只会增加计算时间不会显著改善精度。这个测试帮你找到性价比最高的参数组合。注意收敛性测试要在相同的物理参数下进行。如果你改变了 ( z ) 或 ( a )收敛所需的 ( N ) 也会变化。一般来说菲涅尔数越大近场需要的网格越密。5. 常见问题与排查技巧实录5.1 衍射图样出现异常条纹或网格状伪影这是最常遇到的问题。表现是仿真结果中出现规则的条纹或网格状图案而这些在理论结果中不应该存在。原因通常有三个第一混叠效应。当观察面的采样间隔 ( \Delta x ) 大于衍射图样的最小特征尺寸时高频信息被折叠到低频产生莫尔条纹。解决办法是增加 ( N ) 或调整 ( z )使得 ( \Delta x ) 满足采样定理。第二FFT移位操作错误。如果fftshift和ifftshift用反了频谱会偏移半个网格导致干涉条纹。检查方法是看零频分量是否在数组中心。第三孔径边缘的硬截断。前面提到的软边处理可以缓解这个问题。5.2 近场计算时结果完全失真当传播距离 ( z ) 很小时比如小于孔径半径的平方除以波长菲涅尔近似本身可能就不成立了更不用说数值计算。这时候你需要考虑角谱法或者瑞利-索末菲积分。但如果你确实需要在近场用菲涅尔近似至少要保证采样满足[ \Delta \xi \leq \frac{\lambda z}{L} ]这意味着 ( z ) 越小需要的采样越密。如果 ( z 1\text{mm} )( \lambda 632.8\text{nm} )( L 5\text{mm} )那么 ( \Delta \xi \leq 0.126\text{mm} )对应 ( N \geq 40 )。这还不算太苛刻。但如果 ( z 0.1\text{mm} )( N ) 就要超过400了。5.3 计算速度太慢的优化思路MATLAB的FFT本身很快但如果你的代码里有循环或者重复计算速度会急剧下降。我总结了几条优化经验预计算所有不随循环变化的量比如相位因子、坐标网格等。用fft2而不是嵌套循环这是最基本的。如果只需要轴上强度不需要计算整个二维场可以用一维积分近似速度快几个数量级。用single精度代替double内存占用减半速度提升约30%对于大多数可视化需求精度足够。考虑用GPU加速如果你的MATLAB版本支持gpuArray大网格下的加速比可以到10倍以上。% GPU加速示例需要Parallel Computing Toolbox U_mod_gpu gpuArray(U_mod); U_fft_gpu fftshift(fft2(ifftshift(U_mod_gpu))); U_gpu gather(prefactor * phase_obs .* U_fft_gpu * dxi^2);5.4 常见问题速查表问题现象可能原因排查方法解决方案图样出现规则条纹混叠效应检查 ( \Delta x ) 是否满足采样定理增加N或调整z结果整体偏移FFT移位错误检查零频位置修正fftshift/ifftshift能量不守恒采样不足或截断计算能量比增大L或N近场结果失真菲涅尔近似失效计算菲涅尔数改用角谱法计算速度慢循环或重复计算用profiler分析向量化预计算边缘出现伪影硬边界离散化观察边缘区域使用软边孔径强度量级不对面积元或前置因子遗漏检查公式补上dxi^2和prefactor6. 从仿真到应用这个框架还能怎么扩展圆孔菲涅尔衍射仿真的价值不仅在于验证一个物理现象更在于它提供了一个可扩展的计算框架。你把这个框架稍微改一改就能处理很多相关问题。比如把圆孔换成矩形孔只需要改孔径函数的定义换成环形孔加一个内半径判断就行换成多个圆孔把几个孔径函数叠加即可。甚至你可以做动态衍射——让孔径随时间变化观察衍射图样的演化。再进一步这个框架可以扩展到部分相干光的衍射计算。方法是对不同波长的单色结果做加权叠加权重由光源的功率谱决定。虽然计算量会增加但思路是通的。我在实际项目中还用过这个框架做衍射光学元件的设计验证。比如设计一个菲涅尔波带片先用这个仿真工具验证它的聚焦效果再决定是否投入加工。这种“先仿真后加工”的流程能省下不少试错成本。最后分享一个小技巧如果你需要频繁修改参数做批量仿真建议把核心计算封装成函数输入参数用结构体传递。这样代码更清晰也方便做参数扫描。function [I, x] fresnel_circular(lambda, a, z, N, L) % 封装好的圆孔菲涅尔衍射计算函数 % 输入波长、孔径半径、传播距离、网格点数、孔径面尺寸 % 输出强度分布、观察面坐标 dxi L / N; xi (-N/2 : N/2-1) * dxi; [Xi, Eta] meshgrid(xi, xi); dx lambda * z / (N * dxi); x (-N/2 : N/2-1) * dx; [X, Y] meshgrid(x, x); k 2 * pi / lambda; R sqrt(Xi.^2 Eta.^2); U0 double(R a); U_mod U0 .* exp(1i * k / (2*z) * (Xi.^2 Eta.^2)); U_fft fftshift(fft2(ifftshift(U_mod))); prefactor exp(1i * k * z) / (1i * lambda * z); U prefactor * exp(1i * k / (2*z) * (X.^2 Y.^2)) .* U_fft * dxi^2; I abs(U).^2; I I / max(I(:)); end这个函数可以直接调用也可以作为更复杂系统的基本模块。我个人在实际操作中的体会是把重复使用的计算逻辑封装好比每次从头写要高效得多而且不容易出错。尤其是当你需要做参数扫描或者优化设计的时候一个可靠的底层函数能帮你省下大量时间。