简介资源围绕艾里光束的角谱衍射传播与MATLAB数值模拟展开适合光学专业学生、科研人员以及需要掌握光束仿真技巧的MATLAB使用者。压缩包内共2个文件一个可直接运行的MATLAB脚本用于生成艾里光束、进行角谱变换并模拟不同传播距离下的场分布一篇PDF文档则梳理了衍射理论对局域空心光束及无衍射光束重建的描述补充相关理论基础。整个包仅468KB轻量易用目前已有2170人学习下载可作为课程设计、毕业设计或科研预研的参考素材。通过调整初始光束尺寸、强度分布和传播距离等参数可系统比较不同条件对光束演化路径的影响直观观察艾里光束的无衍射、自弯曲与自修复特性同时掌握fft2做角谱传播仿真的核心思路为光操控、自由空间光通信等应用研究提供前置验证与理论支持。1. 角谱理论配合艾里光束一个能打通光束仿真全流程的实用组合第一次听到“艾里光束”这个词的人多半会以为它是某种高深的量子光学理论。但当你真正把“基于角谱理论的艾里光束传播”这个标题拆开会发现它其实是一条非常完整、也很适合上手的仿真链路艾里光束是物理模型角谱衍射理论是传播算法MATLAB 是落地工具。三者拼起来就能在电脑上复现自弯曲、无衍射、自愈这些听起来像玄学的光束特性。无论你是刚接触光束传播仿真的学生还是需要在实验前做参数预演的工程师这个组合都值得花一个下午跑通。我最早做这个仿真时走了不少弯路最大的问题不是公式难而是网格坐标、频率坐标、物理单位这三者在 MATLAB 里经常打架。把三件事对齐整套仿真就通透了一大半。下面我会从原理、参数、代码、避坑四个角度把它拆开讲清楚尽量做到你照着敲就能出图。2. 先把要算的东西说清楚艾里光束初始场与角谱传播原理2.1 有限能量艾里光束为什么指数截断是标配艾里函数 Ai(x) 在光学里的地位主要来自它作为薛定谔方程无衍射解的性质。理想情况下一个无限能量的艾里光束可以保持横向形状不变地传播主瓣还会沿着一条抛物线轨迹横向加速。但“无限能量”在物理实验里是不存在的——你没法用实验装置产生一个总能量为无穷大的光场。所以实际仿真和实验中使用的是有限能量艾里光束。做法是在初始场里乘一个指数截断因子最常见的形式是φ(x, y, 0) Ai(x / x₀) · Ai(y / x₀) · exp(a · (x y) / x₀)其中 x₀ 是横向尺度因子控制主瓣和旁瓣的宽度a 是指数截断因子通常取 0.01 到 0.1 之间的值。a 越小光束越接近理想艾里光束但旁瓣拖得越长计算窗口也要开得越大a 越大光束能量越集中但自加速和无衍射特性会被削弱。这里有一个容易搞混的地方Ai(x / x₀) 在 x 0 方向是指数衰减的而在 x 0 方向是振荡的。乘上 exp(a · (x y) / x₀) 之后虽然 x 为正时指数因子在增大但艾里函数本身的衰减速度远快于指数增长所以整体能量仍然是有限的。你在设置网格坐标时把主瓣放在 x 负半轴附近振荡尾巴朝 x 负方向延伸这样最直观。二维艾里光束的初始场不是某个二维解析函数而是两个一维艾里函数的乘积。这样做的好处是便于在 MATLAB 里用meshgrid直接构造。很多教程里会直接给出这个乘积形式并没有解释为什么其实核心原因是艾里光束在 x 和 y 方向的自加速是独立的乘积形式正好保留了这种各向异性。2.2 角谱衍射理论将空间域卷积简化为频域乘法角谱法的基本思想是把初始光场看作无数不同方向传播的平面波的叠加。光场在 z 0 处的空间分布经过二维傅里叶变换后得到的就是各平面波分量的复振幅即角谱。当光传播一段距离 z 后每个平面波分量只会获得一个与传播方向相关的相位延迟。在频域里这表现为初始角谱乘以一个传递函数H(f_x, f_y) exp(i · k · z · sqrt(1 - (λ f_x)² - (λ f_y)²))其中 k 2π / λ 是波数λ 是波长f_x 和 f_y 是空间频率坐标。这个公式没有做近轴近似适用于大角度传播如果传播距离远、角度小也可以简化成 H ≈ exp(i · k · z) · exp(-i · π · λ · z · (f_x² f_y²))后者就是常见的菲涅尔近似形式。用角谱法做仿真本质上是把空间域的卷积运算换成了频域的乘法运算。卷积的数值计算复杂度高而乘法配合 FFT 只需要一次fft2、一次ifft2。与直接对衍射积分做双重求和相比角谱法的速度优势在 1024 × 1024 这种网格规模上非常明显直接求和可能要算几分钟角谱法零点几秒就能出结果。但角谱法也有它的前提条件空间频率分量必须满足 λ²(f_x² f_y²) ≤ 1否则 sqrt 里面会变成负数对应的是倏逝波分量。数值实现时如果不处理这些高频分量计算结果可能发散或者出现噪点。我在后面的避坑章节会专门说这个问题。2.3 参数表波长、尺度因子、截断因子、网格大小做角谱仿真最先要确定的是一组基础参数它们之间互相制约尤其是采样间隔和频率范围的关系。下面是我常用的参考值你先按这个跑通再根据自己课题里的需求调整。参数符号典型值作用与取值建议波长λ632.8 nm决定传输函数的指数项改它直接影响加速轨迹横向尺度x₀50 ~ 200 μm控制主瓣宽度太大需要增大窗口截断因子a0.02 ~ 0.1控制能量衰减和旁瓣拖尾长度计算窗口L4 ~ 10 mm必须足够容纳艾里旁瓣的振荡范围采样点数N512 ~ 2048越大越细腻但内存和速度要权衡传播距离z0.1 ~ 1 m超过窗口后主瓣可能跑出视野这里最关键的是采样间隔 dx L / N 和频域间隔 df 1 / L 的关系。空间域窗口 L 越大频域分辨率越好采样点数 N 越大频域范围越宽。艾里光束的旁瓣结构非常精细如果你的 dx 太大高频细节会被直接抹掉表现出来的就是旁瓣数量变少、主瓣变胖。另外一个容易忽视的点是所有长度单位必须统一。我见过不少人在 MATLAB 里把波长写为 632.8e-9却把 x₀ 写成 100没有标注单位结果算出来的加速轨迹完全对不上理论。建议从一开始就用米作为基本单位变量名里带上单位比如x0_m、lambda_m避免后续换算混乱。3. 用 MATLAB 跑通第一次传播一份可以直接复制的完整代码3.1 生成初始场airy 函数与指数截断打开 MATLAB新建一个脚本把下面这段代码粘贴进去。它会生成一个 1024 × 1024 的二维有限能量艾里光束初始场并保存强度图。%% 参数定义单位米 clear; clc; close all; lambda 632.8e-9; % 波长 632.8 nm x0 100e-6; % 横向尺度因子 100 um a 0.05; % 指数截断因子 L 8e-3; % 计算窗口边长 8 mm N 1024; % 采样点数每边 dx L / N; % 空间采样间隔 %% 构造二维网格 x linspace(-L/2, L/2, N); [X, Y] meshgrid(x, x); %% 初始场有限能量艾里光束 % 二维场 两个一维艾里函数乘积再乘截断因子 phi0 airy(X / x0) .* airy(Y / x0) .* exp(a * (X Y) / x0); phi0 phi0 / max(abs(phi0(:))); % 归一化方便后续观察相对强度 %% 显示初始场强度分布 I0 abs(phi0).^2; figure; imagesc(x*1e3, x*1e3, I0); axis square; colormap hot; xlabel(x (mm)); ylabel(y (mm)); title(有限能量艾里光束初始强度);这段代码的逻辑不复杂先用meshgrid生成二维坐标网格然后调用 MATLAB 内置的airy函数直接计算每个网格点上的艾里函数值。注意airy(X/x0)默认返回第一类艾里函数 Ai这正是我们需要的形式。指数截断因子通过exp(a * (X Y) / x0)实现逐元素乘法直接作用在网格矩阵上。参数说明a取 0.05 是折中值旁瓣不至于太长导致窗口放不下自加速特性也保留得够明显。如果你缩小a到 0.01初始场会多出好几个旁瓣视觉上更接近教科书里的艾里光束图但L 8e-3可能不够宽建议同步把L加到 12 mm。归一化那一步不是必须的但强烈建议保留因为后续观察强度分布时相对值比绝对值直观得多。3.2 角谱传播频域相乘与一次逆变换初始场准备好之后真正的传播计算只需几行代码。继续在同脚本里追加%% 角谱传播 z 0.5; % 传播距离 0.5 m k 2 * pi / lambda; % 波数 % 频域坐标注意范围是 [-N/2, N/2-1]/L fx (-N/2 : N/2-1) / L; [FX, FY] meshgrid(fx, fx); % 角谱传递函数 H exp(1i * k * z .* sqrt(1 - (lambda * FX).^2 - (lambda * FY).^2)); % 频域滤波防止倏逝波项导致数值发散 valid (lambda * FX).^2 (lambda * FY).^2 1; H(~valid) 0; % 传播计算fft2 到频域乘传递函数ifft2 回空间域 PHI0 fft2(phi0); PHI_z PHI0 .* H; phi_z ifft2(PHI_z); %% 显示传播后强度 Iz abs(phi_z).^2; figure; imagesc(x*1e3, x*1e3, Iz); axis square; colormap hot; xlabel(x (mm)); ylabel(y (mm)); title(sprintf(传播 %.2f m 后的强度分布, z));频域坐标fx的构造是这段代码的关键。fft2输出的频率排列是从 0 到 N-1 对应的频率而物理上我们希望频率为 0 的分量在数组中心。这里用(-N/2 : N/2-1) / L构造坐标配合后文不额外做fftshift的做法正好和fft2的默认输出顺序对应。传递函数里的sqrt项如果出现负数物理对应的是倏逝波数值上会导致指数发散。我加了一个valid掩码直接将这类高频分量置零这是做角谱传播时一个非常实用的保护措施也避免了许多莫名其妙的结果。3.3 三张必看的图俯视图、截面图、侧视图强度分布图只能看到二维结果但要判断艾里光束的传播特性至少还要看两张图通过主瓣中心的横向截面图以及传播方向上的侧视图。%% 截面图观察主瓣与旁瓣结构 figure; plot(x*1e3, Iz(N/2, :), LineWidth, 1.2); xlabel(x (mm)); ylabel(相对强度); title(y0 处的横向强度截面); grid on; %% 侧视图取不同传播距离观察主瓣横向偏移 z_list [0.1, 0.3, 0.5, 0.8]; figure; hold on; for zz z_list Hz exp(1i * k * zz .* sqrt(1 - (lambda * FX).^2 - (lambda * FY).^2)); Hz(~valid) 0; phi_zz ifft2(fft2(phi0) .* Hz); I_zz abs(phi_zz).^2; plot(x*1e3, I_zz(N/2, :), LineWidth, 1.2, ... DisplayName, sprintf(z %.1f m, zz)); end xlabel(x (mm)); ylabel(相对强度); title(不同传播距离下的截面演化); legend show; grid on;截面图能直接看出主瓣宽度有没有变化、旁瓣位置是否发生了横向移动。侧视图则把不同传播距离下的截面叠在一起你会看到主瓣整体向 x 正方向偏移这就是自加速现象的直观体现。要注意侧视图里的曲线最好做归一化否则不同距离下峰值强度差异很大不容易对比。4. 角谱仿真里的五类翻车现场现象、原因和解决办法4.1 边缘出现周期复制传播结果像被“拼接”过现象传播一定距离后强度图四个边缘出现明显重复的亮斑像把画布卷起来首尾相接。旁瓣从左边消失又从右边冒出来。原因这是 FFT 周期卷积的典型表现。空间域计算窗口不是无限大FFT 默认把窗口当成周期性重复光场穿过窗口一边后会从另一边再进入。艾里光束的旁瓣很宽传播过程中横向扩散后碰到窗口边界就产生了周期伪影。解决最直接的办法是增大计算窗口 L让光场在目标传播距离内始终远离边界。其次是加吸收边界常见做法是在窗口边缘设置一个平滑衰减掩码但这样会增加代码复杂度。我一般先试增大 L 到原来的 1.5 到 2 倍如果仍然有伪影再考虑加吸收层。4.2 中心一个刺眼亮斑艾里旁瓣完全看不清现象初始场看起来正常但传播后整个图中心出现一个特别亮的光斑周围的艾里旁瓣几乎不可见。不管怎么调色标细节都藏在阴影里。原因最常见的是初始场里混入了零频分量。如果 phi0 矩阵均值不为零或者你在做频域处理时没有把直流项处理好fft2 之后 (0,0) 位置的角谱分量会特别大传播后这部分平面波直直地打在中心盖过了旁瓣信号。解决在生成初始场后检查mean(phi0(:))是否接近 0。如果偏离严重可以对初始场做一次去直流处理phi0 phi0 - mean(phi0(:))。另一个办法是调整指数截断因子 a 不要太大a 过大会让光场整体向正方向倾斜直流分量显著增强。4.3 主瓣明明在左边传播后却跑到右边现象初始场主瓣在 x 负方向按照理论应该向右偏移但仿真结果里主瓣一开始就出现在右边而且方向完全反了。原因频域坐标方向与空间域坐标方向的对应关系搞反了。在构造 FX 和 FY 时如果用meshgrid后没有注意维度方向或者多做了一次fftshift又没有配对ifftshift最终结果就会翻转。解决检查频域坐标构造是否和fft2的输出顺序严格对应。一个稳妥的做法是做完一次零距离传播测试把 z 设为 0看ifft2(fft2(phi0))是否精确恢复初始场。如果恢复出来的场左右翻转说明坐标没对齐。另外不要混用fftshift和ifftshift这两者只在奇数 N 时有区别建议只在频域坐标可视化时使用计算链路里保持fft2→ 相乘 →ifft2的原始顺序。4.4 传播距离增大后出现细微网格状振铃现象z 在 0.1 m 左右图面很干净加到 0.5 m 以上强度分布上出现规则的网格状细纹像蒙了一层纱。原因网格状振铃来自频域采样范围不够也就是高频分量在传播过程中被传递函数截断了。当我们设置valid掩码时大于截止频率的分量被直接置零这在频域形成一个矩形截止边界对应的空间域就表现为振铃。距离越大相位变化越剧烈振铃越明显。解决不要直接硬截断改用平滑过渡。例如给传递函数加一个余弦衰减窗让截止区域附近的幅度逐渐降到零。滤波器过渡带越宽振铃越弱。另外增大采样点数 N 也可以提升频域覆盖范围从 1024 加到 2048 往往能明显改善。4.5 主瓣偏移量比理论值小一大截现象把不同传播距离下的主瓣峰值位置提取出来和理论抛物线轨迹公式对比发现偏移量偏小有时甚至只有理论值的一半。原因大概率是 x₀ 或 λ 的单位写错了。艾里光束自加速的横向偏移约等于 (λ² · z²) / (16π² · x₀³)这是一个对 x₀ 特别敏感的量x₀ 差一点偏移量就差很多。另一个原因是初始场里主瓣峰值位置可能不在原点提取“主瓣位置”时用的是强度质心而不是峰值点导致比较基准不一致。解决把物理单位全部换算成米再代入。提取主瓣位置时先用max找到峰值索引再取该索引对应的 x 坐标值。不要用sum(x .* Iz) / sum(Iz)这种质心公式因为质心会被旁瓣的强度分布拉偏。5. 跑通只是第一步用仿真验证自加速和自愈这两个招牌性质5.1 主瓣轨迹拟合用数值结果反推加速曲线验证自加速最直接的方法是记录不同传播距离下主瓣峰值位置然后拟合成抛物线。下面这段代码会输出一条偏移量随传播距离变化的曲线%% 计算不同距离下的主瓣峰值位置 z_list linspace(0.05, 0.8, 16); peak_pos zeros(size(z_list)); for idx 1:length(z_list) zz z_list(idx); Hz exp(1i * k * zz .* sqrt(1 - (lambda * FX).^2 - (lambda * FY).^2)); Hz(~valid) 0; phi_zz ifft2(fft2(phi0) .* Hz); I_zz abs(phi_zz).^2; % 沿 y0 取截面找峰值位置 profile I_zz(N/2, :); [~, peak_idx] max(profile); peak_pos(idx) x(peak_idx); end %% 和理论抛物线对比 theory_pos x0 (lambda^2 .* z_list.^2) / (16 * pi^2 * x0^3); figure; plot(z_list, (peak_pos - peak_pos(1))*1e3, o-, LineWidth, 1.2); hold on; plot(z_list, (theory_pos - theory_pos(1))*1e3, --, LineWidth, 1.2); xlabel(传播距离 z (m)); ylabel(主瓣横向偏移 (mm)); legend(仿真结果, 理论抛物线, Location, northwest); grid on;理论公式里的x0是初始场主瓣的起始位置不是横向尺度因子。用peak_pos - peak_pos(1)可以消除初始位置的偏差让对比集中在“偏移量”上。仿真结果通常在 z 较小的时候和理论高度一致随着 z 增大由于窗口边界的影响和数值采样误差两者会出现明显偏离。这个偏离本身也是信息它告诉你当前参数下角谱法的适用边界在哪里。5.2 遮挡主瓣验证自愈性自愈是艾里光束一个特别有辨识度的特性挡住主瓣光束在传播一段距离后会自动把主瓣“长回来”。验证方法很简单在初始场上人为加一个掩码%% 构造带遮挡的初始场 mask ones(size(X)); mask(abs(X - x0) 0.5*x0 abs(Y - x0) 0.5*x0) 0; % 挡住主瓣附近 phi0_masked phi0 .* mask; %% 传播 0.5 m 后观察恢复情况 zz 0.5; Hz exp(1i * k * zz .* sqrt(1 - (lambda * FX).^2 - (lambda * FY).^2)); Hz(~valid) 0; phi_masked_z ifft2(fft2(phi0_masked) .* Hz); I_masked abs(phi_masked_z).^2; figure; subplot(1,2,1); imagesc(x*1e3, x*1e3, abs(phi0_masked).^2); axis square; colormap hot; title(遮挡后的初始场); subplot(1,2,2); imagesc(x*1e3, x*1e3, I_masked); axis square; colormap hot; title(传播 0.5 m 后的强度);自愈性的物理解释并不神秘主瓣只是艾里光束的一个组成部分旁瓣的能量会随着传播逐渐向主瓣区域转移从而填补遮挡造成的空缺。仿真里你会看到被遮挡的暗区不会一直存在而是逐渐被周围旁瓣“填满”。这个特性在粒子操控和光片显微中有实际应用做完这个验证会对艾里光束的物理图像有更深一层理解。5.3 参数扫描截断因子和尺度因子如何影响传播行为做一次简单的参数扫描可以让数值实验的说服力强很多。我习惯把 a 分别取 0.02、0.05、0.1再把 x₀ 分别取 50、100、200 μm各跑一次 0.5 m 的传播记录主瓣宽度、旁瓣数量、以及峰值位置偏移量。小结一下规律a 越小旁瓣范围越大自加速轨迹越接近理想抛物线但是计算窗口需求也越大a 越大光束能量越集中但“艾里味”越淡。x₀ 越大主瓣在空间上的尺度越大对应的空间频率带宽越窄相同传播距离下的横向偏移越小这符合公式里 x₀³ 在分母上的直觉。这些参数扫描的结果做成表格放在论文里是很扎实的数值证据。6. 把静态仿真变成动态演示一个值得反复使用的收尾技巧每次给合作者看仿真结果时静态强度图总是不够直观。我会把传播距离做成一串数值循环生成强度图并拼接成动画同时在右下角叠加主瓣位置曲线。这个方法不复杂但演示效果很好能让看的人几秒钟就理解“自加速”到底在讲什么。%% 动态传播演示循环 z 并记录主瓣轨迹 z_video linspace(0.05, 0.8, 60); pos_record zeros(size(z_video)); figure(Position, [100, 100, 560, 500]); for idx 1:length(z_video) zz z_video(idx); Hz exp(1i * k * zz .* sqrt(1 - (lambda * FX).^2 - (lambda * FY).^2)); Hz(~valid) 0; phi_zz ifft2(fft2(phi0) .* Hz); I_zz abs(phi_zz).^2; profile I_zz(N/2, :); [~, peak_idx] max(profile); pos_record(idx) x(peak_idx); imagesc(x*1e3, x*1e3, I_zz); axis square; colormap hot; title(sprintf(z %.2f m, zz)); drawnow; pause(0.05); end动画的代码量不大但有一个我反复踩过的坑主瓣位置记录必须在动画循环里同步进行不能等动画跑完再回头补否则就得把每一帧的强度图全存在内存里60 帧 1024 × 1024 的 double 数组会占掉近 500 MB。预先分配pos_record数组边跑边记录是节省内存的好习惯。做完动画后把pos_record和理论抛物线画在同一张图里收尾你会看到一段漂亮的曲线仿真点沿着理论抛物线稳步爬升最后在窗口边界附近偏离。这个偏离并非“错误”而是数值方法和计算窗口共同作用的必然结果。我做这个仿真最大的教训是不要在脚本开头一股脑把所有参数写死而是把波长、尺度、截断因子、窗口大小这些量全部声明成变量并让窗口大小随截断因子和传播距离自动调整。这样每次改参数重新跑代价只是几秒钟的运行时间而这恰恰是探索光学现象最值得投入的地方——改一个参数、看一次结果、记录一条曲线比闷头调半天代码有意义得多。希望这一套流程能帮你也少走几步弯路。本文还有配套的精品资源点击获取
