简介本资源面向光学工程、激光物理及计算光学方向的研究生与科研工程师聚焦艾里光束这一典型无衍射光束的角谱理论建模与MATLAB数值仿真。资源包共2个文件1个PDF理论文献 1个MATLAB源码总大小468KB轻量紧凑便于快速上手。PDF文献系统阐述局域空心光束与无衍射光束的衍射理论重建方法为理解艾里光束的自我恢复与曲率传播提供理论支撑MATLAB脚本完整实现基于角谱法的艾里光束远场传播模拟涵盖初始场构造、频谱演化计算、逆傅里叶变换重构及多距离动态可视化可直接用于光束操控、微粒牵引或光学通信中的非高斯光场设计验证。已有1417人学习下载是兼顾理论深度与工程实践的入门级光学仿真工具包。 艾里光束近几年在光学仿真里快成“网红”了。我第一次看到它的传播动图时确实被那种整体横向偏转、主瓣几乎不发散的画面震住了后来自己用MATLAB复现时才发现这东西看着神奇但背后就是一套经典的傅里叶光学流程角谱理论加一次FFT。这篇博文就完整走一遍基于角谱理论的艾里光束传播模拟从物理图像、参数设计到MATLAB代码实现再到我实际调试中踩过的坑一并写出来。内容适合光学方向的学生、做光束传播仿真的科研人员以及想深入理解角谱法本质的MATLAB用户不求代码一跑就完而是把每一步为什么这么写讲清楚。1. 艾里光束的物理图像与方法选型1.1 艾里光束为什么是“不走直线”的无衍射光艾里光束这个名字来自艾里函数Ai(x)它是自由空间近轴波动方程的一组特殊解。注意它和我们平时见的高斯光束、拉盖尔高斯光束不太一样那几个解在传播中都会慢慢展宽而艾里函数解最大的特点是主瓣在传播过程中保持不变形也就是“无衍射”特性。更直观的体现是它的能量包络整体沿一条抛物线轨迹横向移动仿佛光束在自由空间里自己做了一个匀加速运动这就是所谓的“自加速”。我最早看文献时觉得“自加速”很玄后来用粒子打比方就理解了想象一个溜冰选手在冰面上不是沿直线滑而是带着一个稳定的侧向推力在滑推力大小恒定、方向恒定它的轨迹就是抛物线。艾里光束的“横向加速”就类似这个效果只不过这个横向推力来自光场内部的干涉自弯曲并不需要介质提供梯度折射率。另一个让人印象深刻的特性是“自愈”如果拿一个小障碍物挡住艾里光束的主瓣它在继续传播一段距离后能重新自己“长回来”。这在实际应用中很实用比如微粒操控或光片成像时少量遮挡不会让光束整个报废。1.2 理想艾里函数能量无限仿真必须做指数截断这里有个非常关键的物理细节严格来说理想艾里函数在整个空间上不是平方可积的它的能量是发散的。也就是说理论上那个无限延伸的振荡尾巴把能量带走了根本无法在物理上实现。经典的解决思路是用指数衰减因子给初始场“加窗”把尾巴削掉得到近似的艾里光束。典型的初始场构造方式是E(x, y, 0) Ai(x/x0) * Ai(y/y0) * exp(a * (x/x0 y/y0))其中x0和y0决定主瓣的横向尺度a是衰减系数通常在0.01到0.1之间取值。a越小光束越接近理想艾里函数尾巴越长能量截断越少但数值模拟时边缘伪影会越明显a越大尾部衰减越快主瓣更干净但艾里特性也会变弱。很多论文里喜欢取0.05我自己在仿真中也是从0.05起步再根据结果微调。1.3 为什么用角谱法而不是直接积分传播方法的选择直接决定代码复杂度和运算速度。常见做法有三种空间域直接积分、分步光束传播法、角谱法。我在仿真自由空间传播时基本无脑选角谱法原因很简单它把偏微分方程问题变成了一个简单的“频域相乘”一次FFT加一次IFFT就能把整段距离的光场算出来。方法核心思想优点短板空间域直接积分菲涅尔/瑞利-索末菲衍射积分物理意义直观二维卷积计算量大距离越远越慢分步光束传播法(BPM)把传播路径切成小段逐段衍射相位修正适合折射率渐变介质需要小步长均匀介质中效率低角谱法分解为平面波频域乘以传播相位单次FFT完成大步长精确无近轴强制限制要求均匀介质网格需规则角谱法的数学本质是把初始光场看成无数不同方向传播的平面波的叠加每个平面波在自由空间传播一段距离z后只是累积一个相位exp(ikzz)其中kz是纵向波数由横向空间频率决定。整个算法一句话概括对初始场做FFT乘以传播传递函数再做IFFT回来。这里没有任何近轴近似只要初始场本身是近轴近似下的艾里解用它来传播就已经足够精确。2. 动笔写码前先把参数设计搞清楚2.1 波长、横向尺度、衰减系数怎么选很多人拿到代码第一件事是复制粘贴结果换个参数就“满屏条纹”问题往往出在参数没有闭环。我先列一组我常用的起始参数波长lambda 532e-9绿光532nm实验室常见光源也方便和文献对比。横向尺度x0 50e-6也就是50微米。这个值决定主瓣宽度经验上主瓣半高宽大概在2到3倍x0。x0取得太小光束加速太快模拟窗口根本装不下x0取得太大主瓣太宽无衍射距离虽然长但视觉效果不明显。衰减系数alpha 0.05需要在“尾巴长度”和“边缘干净度”之间平衡。这三个量不是孤立的。横向位移的估算公式是x_shift z^2 / (4 * k^2 * x0^3)其中k就是波数2π/λ。这个公式来自艾里光束峰值位置的解析表达式也就是艾里函数自变量括号中与z平方相关的那一项。我每次仿真前都会先算一遍确保最大传播距离上的横向位移不会超过模拟窗口的一半否则光场跑出边界后还会从另一边折回来造成严重的混叠干扰。举一个具体例子lambda 532nmx0 50μm传播1米。先算k 2π / 532e-9 ≈ 1.18e7k的平方约1.40e14x0的三次方为1.25e-10于是k^2 * x0^3 ≈ 1.75e4。z^2除以4倍这个值结果约为1.43e-5米也就是14.3微米。看起来不大但如果你把x0缩小到10μm同样传播2米位移会暴涨到几米根本模拟不了。这就是为什么很多环境下的艾里光束实验要用透镜补偿或把尺度做大。2.2 网格与采样这条比物理参数更容易坑人角谱法在MATLAB里的本质是离散傅里叶变换所以采样约束绕不开。设模拟窗口尺寸为L采样点数为N那么空间采样间隔dx L/N频域间隔df 1/L。这里有个容易出现混叠的条件角谱法中可准确表示的横向空间频率范围是[-N/(2L), N/(2L)]对应的最大横向波数kx_max π * N / L。如果初始场的高频分量超过了这个范围就会被折叠产生虚假干涉条纹。有个经验法则N必须是2的幂次至少1024有条件就直接上2048。L的选择则要看光束横向范围。艾里光束的尾部虽然被指数衰减削掉了但如果alpha选得太小尾部振荡延伸很远L不够大就会产生截断伪影。我的经验是L至少取到x0的80到100倍以上配合N2048基本能覆盖大多数常见参数。另外一个容易忽略的细节是频域坐标的顺序。MATLAB的fft2输出频谱顺序是零频在数组的左上角不是中心。如果你按自然顺序用meshgrid方式构造频率轴比如fx (-N/2 : N/2-1) / L那么构造出的传递函数H必须做一次ifftshift再和fft2的结果相乘否则频谱对不齐出来的光场混乱不堪。这个点我在第3部分代码里会明确写出来。2.3 传播距离的分段策略和能量校验角谱法虽然支持大步长但不是“无限大步长”都行。原因在于传递函数exp(ikzz)是相位高度振荡的函数当z非常大时H在频域的变化非常快如果频域采样不够密乘出来的结果会严重失真。这有点像用离散点采样一个高频正弦波采样率不够就出现假的低频信号。所以在实际模拟中我通常不会直接用一次FFT传播几十米而是把传播路径切成若干段每段调用同一个角谱传播函数逐段累加。这样做有两个好处一是每段的相位变化量可控不容易触发混叠二是可以顺便输出每个截面的光强得到传播动画。能量守恒是一个特别好的自检指标传播前后sum(abs(E).^2)*dx^2应该几乎不变如果能量出现了明显跳动说明网格或参数有问题。3. MATLAB核心代码实现与逐段拆解3.1 参数初始化代码这里我直接给出一个完整可跑的版本。先说明环境MATLAB R2019b及以上版本都行不需要额外工具箱。% 基础参数 lambda 532e-9; % 波长 532nm k 2 * pi / lambda; % 波数 x0 50e-6; % 艾里光束横向尺度 y0 50e-6; alpha 0.05; % 截断系数 % 网格参数 L 8e-3; % 模拟窗口尺寸 8mm N 2048; % 采样点数 dx L / N; x (-N/2 : N/2-1) * dx; [X, Y] meshgrid(x, x);这里有个习惯问题我建议无论初始场是什么都把空间坐标原点放在窗口中心也就是用(-N/2 : N/2-1)来构造坐标。这样艾里函数Ai(x/x0)的自变量在原点附近对称后续看光强分布时也比较直观。如果你用(0:N-1)构造初始场中心会偏到窗口角落给后面的调试添乱。3.2 初始场构造与可视化自检构造初始场这步看似简单但很容易写错。二维艾里光束是x和y方向一维场相乘的形式不是把Ai(x y)这样的耦合形式。我见过有人直接写成airy(X Y)出来的图案完全是斜向条纹和真正的艾里光束差很远。正确写法如下% 初始艾里光场 E0 airy(X / x0) .* airy(Y / y0) .* exp(alpha * (X / x0 Y / y0)); % 检查初始光强分布 figure; imagesc(x * 1e3, x * 1e3, abs(E0).^2); axis xy; axis equal tight; colormap(parula); colorbar; title(初始光强分布); xlabel(x / mm); ylabel(y / mm);跑完这段你应该能看到一个沿左上到右下对角线方向排列的明亮主瓣周围带有一串扇形展开的次瓣。如果你看到主瓣不在中心区域或者图案不对称优先检查x0、y0和坐标构造。关于airy函数MATLAB中airy(X)就是Ai(X)等价于airy(0, X)这个返回的是实轴上的艾里函数值。注意不要写成airy(1, X)那是Ai的导数用于其他场景。另外艾里函数在负自变量区域是振荡的MATLAB处理没问题但如果你需要精确控制精度建议核对一下负半轴那段振荡是否符合预期。3.3 频域传递函数的构造与频谱对齐这是角谱法最核心的一步。标准传递函数是H(fx, fy) exp(i * z * sqrt(k^2 - (2πfx)^2 - (2πfy)^2))如果严格按照这个写需要注意频域坐标的顺序。第一种做法是按fft输出顺序直接构造频率避免用fftshift类函数% 按fft输出顺序构造频域坐标 fx_fft (0 : N-1) / L; [FX_fft, FY_fft] meshgrid(fx_fft, fx_fft); kz sqrt(k^2 - (2 * pi * FX_fft).^2 - (2 * pi * FY_fft).^2); H_fft exp(1i * kz * z);这样直接E1 ifft2(fft2(E0) .* H_fft)就完事但缺点是频率轴不是从负到正排的不容易一眼看出哪是零频。我更推荐第二种做法先按自然顺序构造再用ifftshift对齐% 频域坐标从负到正零频在中心 fx (-N/2 : N/2-1) / L; [FX, FY] meshgrid(fx, fx); % 严格角谱传递函数 kz sqrt(k^2 - (2 * pi * FX).^2 - (2 * pi * FY).^2); H exp(1i * kz * z); % 对齐到fft2输出顺序 H ifftshift(H); % 一次传播 E1 ifft2(fft2(E0) .* H);简单解释一下为什么必须加ifftshift。fft2(E0)返回的频谱零频位于数组的(1,1)位置而按自然顺序构造的H零频位于数组中心(N/21, N/21)。ifftshift把H的原点从左下和右上区域整体搬到数组角落正好匹配fft2的输出顺序。这个操作如果漏掉模拟结果就会变成一堆斜纹而且怎么调参数都救不回来。我在初学阶段被这个坑折磨过整整一下午。如果你更习惯用近轴近似也可以把kz换成抛物线近似% 近轴近似传递函数 H exp(-1i * z * (2*pi*FX).^2 / (2*k)) .* exp(-1i * z * (2*pi*FY).^2 / (2*k)); H ifftshift(H);不过既然都用严格角谱了我建议直接用第一版。艾里光束本身是近轴方程的解理论上两种版本差别非常小但严格版在网格很大时更稳。3.4 多距离传播与动态可视化单次传播只是热身实际研究总需要看光束在连续距离上的演化。我会把整个传播过程写成一个循环每次输出一帧图像方便导出动画z_list linspace(0, 1.5, 31); % 从0到1.5米共31个截面 for idx 1 : numel(z_list) z z_list(idx); H exp(1i * kz * z); H ifftshift(H); Ez ifft2(fft2(E0) .* H); % 光强 Iz abs(Ez).^2; % 绘制 imagesc(x * 1e3, x * 1e3, Iz); axis xy; axis equal tight; clim([0, max(Iz(:)) * 0.8]); % 压一下动态范围次瓣更明显 colormap(parula); title(sprintf(z %.3f m, z)); xlabel(x / mm); ylabel(y / mm); drawnow; % 如果想保存视频 % frame getframe(gcf); % writeVideo(v, frame); end注意clim那行我刻意把显示上限压到最大值的80%因为艾里光束的次瓣强度比主瓣低不少不压缩动态范围会显得主瓣过曝、次瓣一片黑。MATLAB R2022a之后的版本建议用clim替代已经废弃的caxis。3.5 能量守恒校验一个动作排查一半问题无论参数怎么改我建议在循环外先跑一次能量校验E0_energy sum(abs(E0(:)).^2) * dx^2; Ez_energy sum(abs(E1(:)).^2) * dx^2; fprintf(初始能量: %.6f, 传播后能量: %.6f, 相对误差: %.3e\n, ... E0_energy, Ez_energy, abs(E0_energy - Ez_energy) / E0_energy);角谱法在均匀介质中理论上严格能量守恒所以你测出来的相对误差通常应该小于1e-10。如果误差到了1e-2级别基本可以断定网格或传递函数有问题。这个检查做起来只需要两秒但它能帮你快速区分是参数问题还是算法问题。4. 常见问题与排查技巧实录4.1 满屏斜条纹和”折叠图像“多半是频谱对齐或混叠我见过太多人把仿真结果贴到论坛上问“为什么我的艾里光束像一坨乱线”十有八九是频谱对齐问题。你在代码里写了fftshift还是ifftshift是加在频域还是空间域一个符号的差别就能毁掉整张图。我的排查顺序是先算能量守恒如果能量完全不变但图案错乱那就只剩对齐问题把传递函数中的ifftshift去掉或改成fftshift再试试。另外一种是混叠表现为图案边缘出现规则波纹这种往往是因为网格L太小艾里尾巴被硬截断导致的。解决办法是把L增大或者把alpha从0.05提到0.08让尾巴更快衰减。4.2 主瓣位置和理论值对不上先检查你把位移公式算对没有艾里光束主瓣位置的解析结果是x_peak z^2 / (4 * k^2 * x0^3)。有次我用x050μm、z1m算出来的理论位移只有14微米结果在图上几乎看不出来还以为代码错了。后来在网格里量了一下实际峰值位置对得相当好。这里提醒一句如果你的x0取的是微米量级位移也是微米量级在毫米级窗口里当然看不出明显移动。想看明显的抛物线轨迹要把传播距离拉长或者用更小的x0但也要保证窗口够大。还有一种情况是二维时x和y方向都移动如果你只看x方向剖面容易把对角线方向位移误判成没移动。4.3 关于AIRY函数在MATLAB中的几个细节airy(X)计算的是Ai(X)支持向量和矩阵输入计算速度很快。但有几个边界情况需要注意。第一airy对NaN和Inf的处理和一般函数一样但如果你把X/x0算出来是Inf那airy(Inf)会返回0可能导致整个场变成黑屏这种基本是坐标或x0的单位搞错了。第二MATLAB的airy函数在负实轴上是实数值不用特殊处理。但如果你不小心传了复数结果就会变成复数光场会莫名出现相位所以尽量确保X/x0是纯实数。第三二维构造时我用的是airy(X/x0) .* airy(Y/y0)两个矩阵点乘。要是误写成矩阵乘法*MATLAB会报维度错误这也算一个常见的“隐形提醒”。4.4 一张速查表搞定大部分问题现象可能原因解决办法能量严重不守恒传递函数对齐错误或kz虚部处理不当检查ifftshift改用严格角谱版画面全是斜条纹忽略频谱对齐H构造后加ifftshift边缘出现规则波纹窗口L太小增大L或增大alpha主瓣看不出横向位移x0偏大或z太短增大z或减小x0重新估算位移量图案中心不在原点坐标轴没从(-N/2:N/2-1)开始调整空间坐标构造光场出现NaN或Infairy输入为Inf检查x0的单位和数量级这些经验都是从实际调试里一条条攒出来的尤其是ifftshift这个问题属于角谱法新手必踩的坑今天写出来希望帮你省掉几个小时的暴力试错。5. 几个值得尝试的进阶方向5.1 实验版产生方式用高斯光束加三次相位做傅里叶变换模拟归模拟真正在实验室里要产生艾里光束常用的办法是对高斯光束做“整形”。核心原理是艾里函数的傅里叶变换是三次相位因子exp(i * u^3 / 3)所以把高斯光束通过一个立方相位板再用透镜做傅里叶变换就能近似得到艾里光束。我后来在模拟里验证过这条链路代码上也很容易复现先给高斯初始场加上三次相位exp(1i * beta * (X.^3 Y.^3))然后用透镜傅里叶变换公式传播到焦平面得到的场和直接用airy函数构造的初始场非常接近只是参数之间要做换算。想深入理解的话建议把这个链路加入你的仿真中对比“直接构造”和“实验模拟”两种方式的异同对提升物理直觉非常有帮助。5.2 加损耗介质或线性势场看看艾里光束如何被调教艾里光束最出圈的应用之一是在光敏折变晶体中产生弯曲的路径。在均匀介质里艾里的自加速曲线是固定的抛物线但如果在数值模型中加一个线性折射率梯度或增益损耗项就能人为控制这个加速方向。做法也不难在角谱传播的每一步之间插入一个吸收或相位屏等价于在传播方程里加上势能项。不过这就超出纯角谱的范畴了需要用分步法的思想把角谱传播和相位调制交替进行。这个方向水挺深但可玩性极高想深入的同学可以从“线性势中的艾里波包”这个关键词入手。5.3 非傍轴、涡旋艾里光束和其他扩展艾里光束家族远不止一个“主瓣加尾巴”这么简单。可以把艾里函数和光涡旋结合得到艾里涡旋光束它同时携带轨道角动量横向自加速的空间模式更丰富。也可以考虑两个艾里光束对撞的干涉场景这在微粒操控中很有意思。这些扩展在模拟层面并不复杂基本都是改一改初始场传播部分完全复用角谱法代码。我个人的体会是角谱法这套工具就像一把瑞士军刀学会了以后不只是模拟艾里光束任何自由空间中光场的传播都可以用它快速搞定。每次看到有人还在用双重循环做衍射积分我都想劝他先把FFT和频谱对齐吃透后面能省下大把时间。最后再分享一个小习惯我每次改参数都会把z_list、x0、alpha这些关键值连同结果图文件名一起记录下来时间久了这就是一笔很宝贵的参数调试资产比到处翻聊天记录强多了。本文还有配套的精品资源点击获取