Matlab涡旋光束仿真详解：LG光束、角谱传播与拓扑荷检测

Matlab做涡旋光束仿真这件事我是从研究生阶段啃到现在的。刚接触时对着公式一头雾水等真正把Laguerre-Gaussian光束的螺旋相位在Matlab里画出来的那一刻才明白为什么这种光束能拿诺贝尔奖、能成为光通信OAM复用的核心载具。这篇内容我就围绕“常见涡旋光束仿真”来展开把LG光束、螺旋相位屏、菲涅尔传播、湍流影响和拓扑荷检测这些最常用到的仿真模块一次讲清楚代码直接可跑大家拿去就能复现过程中踩过的坑也会逐一说明。Matlab做涡旋光束仿真的核心是把复振幅的实部和虚部都完整表达出来。涡旋光束的特殊之处在于相位分布是螺旋状的中心存在相位奇点导致光强分布呈中空环形。很多人一开始只画了强度图根本看不出涡旋结构必须同时出强度图和相位图才能真正验证光束的拓扑荷和螺旋波前。这套仿真的应用场景相当广从光镊操控、超分辨成像到自由空间光通信的OAM模式复用全都离不开涡旋光束的数值模拟做光学仿真、通信物理层研究或者光学实验课设的朋友都能用得上。1. 涡旋光束仿真到底在仿什么1.1 涡旋光束的核心物理图像涡旋光束最典型的特征是中心存在相位奇点光场在此处相位不确定因此强度严格为零形成暗核。无论是Laguerre-Gaussian光束、高阶贝塞尔光束还是其他涡旋光场物理本质上都围绕一个核心相位因子( e^{il\phi} )其中 ( l ) 称为拓扑荷是一个整数也是决定涡旋光束行为最关键的参数。( \phi ) 是方位角绕着光轴方向从0到2π旋转一圈相位就改变了 ( 2\pi l )相当于波前被扭成了螺旋阶梯状。用生活化的方式去理解如果把普通高斯光束的波前比作一层平整的纸那么涡旋光束的波前就像旋转楼梯从中心往外走波前高度随方位角连续爬升。拓扑荷 ( l ) 代表爬楼梯时一层跨越了几个台阶高度每跨过一整圈波前就积累了 ( 2\pi l ) 的相位差。这个“螺旋楼梯”结构带来了一个重要结果光场中心处相位无法定义强度必须归零所以在远场或聚焦后观察涡旋光束会呈现中心暗核的环形光斑。在Matlab仿真里我们并不需要像COMSOL那样解完整的电磁场方程而是用标量衍射理论和解析模式公式直接构造复振幅场 ( U(r,\phi,z) )然后用傅里叶光学的方法做传播计算。这样做精度足够覆盖大多数科研和工程验证场景而且速度非常快。1.2 为什么要用Matlab来做这套仿真Matlab做涡旋光束仿真有三个不可替代的优势一是矩阵运算天然适配二维光场。一束光在某个横截面上的分布就是一个二维复数矩阵强度是模方相位是辐角Matlab直接在矩阵层面运算省去了一层一层写循环的麻烦。二是可视化能力极强。imagesc、surf、mesh配合伪彩色图能快速把强度、相位、干涉条纹呈现出来这在物理概念验证阶段太重要了。三是FFT快速傅里叶变换函数成熟稳定。衍射传播计算的核心就是傅里叶变换Matlab的fft2/ifft2加上fftshift/ifftshift组合几行代码就能实现角谱传播法。我实际仿真中绝大多数场景都是用LG模解析公式角谱传播法完成的。无论用户关注的是实验室里的空间光调制器SLM产生的涡旋光场还是自由空间链路中受湍流影响后的光强闪烁和模式串扰这套仿真路径都能覆盖。而且Matlab生态里优化工具箱、深度学习工具箱还支持进一步的数据分析比如把多个拓扑荷模式混叠后用神经网络解调这也是现在的热门方向。2. 常见涡旋光束的生成方法核心实操2.1 LG光束的场分布公式与Matlab实现常见涡旋光束里面Laguerre-Gaussian光束最经典也是绝大多数论文和实验的起点。LG光束在柱坐标下的归一化复振幅分布可以写成[ LG_{p}^{l}(r,\phi,z) \sqrt{\frac{2p!}{\pi (p|l|)!}} \frac{1}{w(z)} \left( \frac{\sqrt{2}r}{w(z)} \right)^{|l|} L_p^{|l|}\left( \frac{2r^2}{w(z)^2} \right) \exp\left( -\frac{r^2}{w(z)^2} \right) \exp\left( -il\phi \right) \exp\left( ikz \right) \exp\left( \frac{-ikr^2}{2R(z)} \right) \exp\left( i(2p|l|1)\zeta(z) \right) ]这个公式看着吓人但仿真中真正每次都要算的关键项其实就几个( w(z) w_0 \sqrt{1(z/z_R)^2} )光束半径( L_p^{|l|}(x) ) 是广义拉盖尔多项式( p ) 是径向指数( \exp(-il\phi) ) 是螺旋相位项( l ) 就是拓扑荷( \exp(i(2p|l|1)\zeta(z)) ) 是Gouy相位( \zeta(z)\arctan(z/z_R) )。对于涡旋光束仿真大多数情况下取 ( p0 ) 就已经能体现全部涡旋特性因为拉盖尔多项式 ( L_0^{|l|}(x)1 )公式瞬间简洁很多。我把写好的生成函数直接放出来这个函数支持任意拓扑荷和径向指数function [U, x, y] lg_beam(N, w0, l, p, lambda, z) % LG光束场生成函数 % N: 网格像素数; w0: 束腰半径(m); l: 拓扑荷; p: 径向指数; % lambda: 波长(m); z: 传播距离(m) L 2e-2; % 物理尺寸, 根据束腰调节, 一般取束腰4-6倍 x linspace(-L/2, L/2, N); y x; [X, Y] meshgrid(x, y); [Phi, R] cart2pol(X, Y); zR pi * w0^2 / lambda; % 瑞利距离 w w0 * sqrt(1 (z/zR)^2); % z处的束腰 Rz z * (1 (zR/z)^2); % 波前曲率半径 (z0时无穷大,代码里特判) zeta atan2(z, zR); % Gouy相位 % 拉盖尔多项式值, 广义版用 laguerreL 需要符号工具箱, 这里手写定义 % 对于p0和p1最常用: if p 0 Lag ones(size(R)); elseif p 1 Lag 1 abs(l) - 2 * R.^2 / w^2; else % 一般场景p3, 用递推关系; 更多阶可定义递归函数 Lag laguerre_gen(p, abs(l), 2*R.^2/w^2); end % 径向归一化因子 C sqrt( 2*factorial(p) / (pi * factorial(pabs(l))) ); norm_factor C / w; % 复振幅 U norm_factor * (sqrt(2)*R/w).^abs(l) .* Lag .* exp(-R.^2/w^2) ... .* exp(-1i*l*Phi) .* exp(1i*k*z) ... .* exp(-1i*pi*R.^2/(lambda*Rz)) .* exp(1i*(2*pabs(l)1)*zeta); % 处理z0处半径无穷大的问题: Rz赋大值 if z 0 U norm_factor * (sqrt(2)*R/w0).^abs(l) .* Lag .* exp(-R.^2/w0^2) ... .* exp(-1i*l*Phi); end % 归一化保证总功率为1 U U / sqrt(sum(sum(abs(U).^2)) * (x(2)-x(1))^2); end注意几个细节网格中心用linspace(-L/2, L/2, N)构造这样中心点严格位于矩阵中心cart2pol得到方位角Phi和径向距离Rz0时波前曲率半径无穷大需要特殊处理归一化通过总功率完成保证不同参数下仿真结果可以横向比较。2.2 用螺旋相位屏快速生成各类涡旋光束实际工程里不总需要完整LG解很多场景只是想在已有光场比如平面波或高斯光上套一个涡旋相位最典型的做法就是用螺旋相位屏Spiral Phase PlateSPP的相位掩模[ U_{out}(x,y) U_{in}(x,y) \cdot \exp(i l \arctan2(y,x)) ]这种方式在仿真空间光调制器SLM时非常实用因为SLM本质上就是一个可编程相位掩模。而且它不局限于高斯光你想让Airy光束带上涡旋相位也可以直接把螺旋相位乘进去就行。代码非常简单function SPP spiral_phase(N, l) % 生成螺旋相位屏 x linspace(-1, 1, N); % 归一化坐标 [X, Y] meshgrid(x, x); SPP exp(1i * l * atan2(Y, X)); end然后生成涡旋光束只需要N 512; w0 1e-3; % 束腰1mm lambda 632.8e-9; % 氦氖激光波长 % 初始高斯光 x linspace(-5e-3, 5e-3, N); [X, Y] meshgrid(x, x); R2 X.^2 Y.^2; G exp(-R2 / w0^2); % 叠加上拓扑荷l2的螺旋相位 l 2; U0 G .* spiral_phase(N, l); % 画强度 I0 abs(U0).^2;用这个方法改变拓扑荷比直接调LG参数更直观。而且拓扑荷的符号也有意义正负代表涡旋旋转方向螺旋相位屏的方法在实验上对应SLM加载的掩模仿真和实验能保持高度一致。2.3 高阶贝塞尔-高斯涡旋光束的补充除了LG光束高阶贝塞尔-高斯光束在“无损传输”“自重建”等场景中非常常见。贝塞尔涡旋光束在传播过程中中心暗斑尺寸几乎不变这在粒子囚禁领域有独特价值。生成方式是高斯光乘以零阶或高阶贝塞尔函数 ( J_l(k_t r) )再乘螺旋相位[ U(r,\phi) J_l(k_t r) \exp(-r^2/w_g^2) \exp(il\phi) ]其中 ( k_t k \sin\theta ) 是横向波矢( \theta ) 是锥角决定了中心亮环的半径。Matlab里调用besselj(l, kt*R)就能生成。和LG一样贝塞尔涡旋的相位图同样是螺旋结构但强度分布不太一样LG是中空高斯-like分布贝塞尔涡旋中心暗核更锐利旁瓣更丰富。3. 涡旋光束的传播仿真与强度、相位重构3.1 角谱法实现涡旋光束的衍射传播生成初始光场后最重要的一步是传播。涡旋光束经过一段距离后强度图样会演化这个衍射过程在Matlab中的标准做法是角谱法Angular Spectrum Method。角谱法的思路是把光场分解为无数不同方向的平面波每个平面波传播一段距离只改变相位不改变振幅分布然后叠加回来。具体表达为[ U(x,y,z) \mathcal{F}^{-1} \left{ \mathcal{F}{U(x,y,0)} \cdot H(f_x,f_y;z) \right} ]传递函数 ( H(f_x,f_y;z) \exp\left( ikz \sqrt{1 - \lambda^2 f_x^2 - \lambda^2 f_y^2} \right) )其中 ( f_x, f_y ) 是空间频率。这个方法的优点是没有近轴近似中的距离限制只要网格采样满足奈奎斯特条件短距离长距离都能算。实现代码function Uz propagate_ASM(U0, lambda, L, z) % 角谱传播 % U0: 起始面复振幅; lambda: 波长; L: 物理尺寸; z: 传播距离 N size(U0, 1); dx L / N; fx (-N/2 : N/2-1) / (N*dx); % 空间频率坐标 [FX, FY] meshgrid(fx, fx); H exp(2i*pi*z/lambda * sqrt(1 - (lambda*FX).^2 - (lambda*FY).^2)); Uf fftshift(fft2(U0)); % 先搬到中心 Uf Uf .* H; Uz ifft2(ifftshift(Uf)); end这一版代码有三个必须提醒的坑一是先用fftshift还是先用fft2的顺序。推荐对初始场执行fftshift(fft2(U))这样频域坐标用(-N/2:N/2-1)这组正负频率坐标对应方便构建传递函数。如果顺序反了频谱中心和传递函数中心错位结果直接错乱。二是传递函数中的倏逝波问题。当 ( 1 - \lambda^2 f_x^2 - \lambda^2 f_y^2 0 ) 时根号内为负传递函数变成衰减项。理论上这是倏逝波在距离远大于波长的场景下迅速衰减可以直接置零或保留复数形式。仿真时默认不会出现太大问题但如果网格很密、频率分量很高要注意数值稳定性。三是网格尺寸和采样间隔匹配。角谱法要求传播距离不能超过一个上限否则会出现混叠。经验公式是 ( z N dx^2 / \lambda )按这个控制仿真参数一般不会出问题。3.2 涡旋光束传播后的强度与相位分析我用一个实际仿真来演示完整流程这个例子也是大家最容易复现的标准场景% 参数设置 N 512; % 网格数 lambda 632.8e-9; % 波长 w0 1e-3; % 束腰 L 1e-2; % 物理尺寸10mm z 2; % 传播2m l 1; % 拓扑荷 p 0; % 生成初始LG光束 [U0, x, y] lg_beam(N, w0, l, p, lambda, 0); % 传播 Uz propagate_ASM(U0, lambda, L, z); % 强度与相位图 figure(Position, [100 100 1200 500]); subplot(1,2,1); imagesc(x*1e3, y*1e3, abs(Uz).^2); axis xy image; colormap hot; colorbar; title(传播2m后的强度分布); xlabel(x (mm)); ylabel(y (mm)); subplot(1,2,2); imagesc(x*1e3, y*1e3, angle(Uz)); axis xy image; colormap(parula); colorbar; title(传播2m后的相位分布); xlabel(x (mm)); ylabel(y (mm));运行后你会看到强度图是典型的甜甜圈形状中心暗核清晰相位图则是围绕中心从蓝色到黄色渐变一圈的螺旋阶梯。拓扑荷 ( l1 ) 时相位刚好转一圈相位图上有1个2π跳跃( l2 ) 时转动两圈。相位图中奇点位置就是强度暗核中心。从这个结果可以引申出一个仿真中的重要经验判断涡旋光束生成对不对不要只看强度一定要看相位图。强度中空结构可能是高阶模式或失焦造成的只有相位图呈现干净的螺旋状才能确认光束确实携带轨道角动量。4. 大气湍流影响与拓扑荷检测4.1 大气湍流相位屏的生成与叠加涡旋光束在自由空间传输时大气湍流会造成相位随机扰动导致模式串扰和光强闪烁。这也是目前光通信OAM复用研究的核心问题。仿真时通常用多层随机相位屏加角谱传播来模拟整个湍流路径这里我给出最常用的Kolmogorov谱相位屏生成方法采用谐波叠加法function phase_screen kolmogorov_phase_screen(N, dx, Cn2, L0) % 生成Kolmogorov湍流相位屏 % Cn2: 折射率结构常数; L0: 外尺度 % 空间频率坐标 fx (-N/2 : N/2-1) / (N*dx); [FX, FY] meshgrid(fx, fx); f sqrt(FX.^2 FY.^2); f(f0) eps; % 避免除零 % Kolmogorov功率谱 (含外尺度修正) r0 0.185 .* (lambda^2 / Cn2 / z_path)^(3/5); % 需要先算Fried参数 % 实际使用时常直接用功率谱: PHI 0.023 * r0^(-5/3) * f .^ (-11/3); % 这里省略外尺度项 % 随机复数高斯场 randn_phase (randn(N) 1i*randn(N)) / sqrt(2); % 滤波 filtered randn_phase .* sqrt(PHI) * N * dx; % 逆变换得到相位屏 phase_screen real(ifft2(ifftshift(filtered))); end这段代码为了演示思路做了简化实际仿真时还需要考虑次谐波补偿低频成分否则相位屏的低频湍流能量会不足。另一个关键点是Fried参数 ( r_0 )[ r_0 0.185 \left( \frac{\lambda^2}{C_n^2 z} \right)^{3/5} ]( r_0 ) 越小湍流越强。在标准大气湍流条件下可见光波段的 ( r_0 ) 一般是几厘米到十几厘米而实验室SLM光束口径一般1mm左右也就是说在实验中湍流影响往往远弱于真实大气链路仿真时要注意参数设定不要仿真得太弱以至于模式串扰根本看不出来。叠加方式更简单每个传播平面上让光场乘以exp(1i * phase_screen)然后继续用角谱法传到下一个相位屏。用了湍流后原本干净的甜甜圈强度分布会出现破碎和闪烁相位图也不再是完美螺旋这会直接影响接收端的拓扑荷识别。4.2 干涉法检测涡旋光束拓扑荷涡旋光束检测是个经典话题而干涉法是理解涡旋结构最直观的方式。实验室中常用的马赫-曾德尔干涉仪是把涡旋光束和平面参考光同轴或离轴叠加形成叉形干涉条纹。仿真同样可以做% 涡旋光束与平面波同轴干涉 l 2; U_vortex G .* spiral_phase(N, l); U_ref ones(N, N); % 平面参考光振幅为1 % 离轴干涉: 参考光倾斜, 相当于叠加一个线性相位 kx 0.2; % 倾斜因子 [X, Y] meshgrid(linspace(-1, 1, N)); U_ref_tilt exp(1i * 2*pi * kx * X); I_interf abs(U_vortex U_ref_tilt).^2; % 干涉条纹会从中心分出l条叉形分叉涡旋光束和倾斜平面波干涉时会形成叉形干涉条纹分叉数等于拓扑荷数。拓扑荷为2时中心会出现一个双叉结构每条干涉条纹从中心分开成两条整体看起来像一把叉子——这就是实验上快速判别涡旋光束拓扑荷的标准方法。这一步的物理解释非常干脆涡旋光束的螺旋波前与参考平面波干涉时等相位面交汇的位置形成条纹。相位围绕奇点旋转2πl干涉条纹就会在这个方向上多出l条分叉。仿真中注意倾斜角不要太大否则条纹太密采样数不够会混叠太小则无法看清分叉结构一般取整个视场内有10~20条条纹比较合适。5. 常见问题排查与实操经验分享5.1 相位图出现噪点与NaN问题我刚开始仿真的时候相位图中心总是一团噪点看起来像“彩色的盐粒”后来检查发现是数值精度问题。涡旋光束中心是相位奇点复振幅实部和虚部都接近0计算幅角angle(U)时微小数值噪声会被放大导致相位随机跳变。这不完全是bug物理上本来如此但会影响展示效果。解决办法是绘图时把中心附近区域强度低于阈值的地方用NaN代替或者固定为某个颜色I abs(U).^2; phase angle(U); phase_thresh phase; phase_thresh(I 0.02*max(I(:))) NaN; % 暗核区域设为NaN imagesc(x, y, phase_thresh);另一个常见现象是相位图出现莫名其妙的横条纹这通常不是涡旋光束的问题而是meshgrid的坐标顺序没搞对。Matlab默认的meshgrid(x, y)生成矩阵时第一维对应y轴如果直接imagesc(angle(U))看起来是旋转了90度的记得加axis xy以及用meshgrid(x, y)而不是meshgrid(y, x)。5.2 采样点数、物理尺寸与传播距离的匹配仿真参数的选择直接决定结果可信度。网格数N、物理尺寸L、波长lambda、传播距离z四者之间必须满足角谱法的采样条件[ z \frac{N \Delta x^2}{\lambda} ]其中 ( \Delta x L/N )。也就是说网格越细、波长越长能传播的距离上限越大。如果传播距离超出该值仿真结果会出现严重的边缘折叠误差。我给出的建议是先用一个自己熟悉的参数组合跑通再逐步调整。比如N512,L1cm,lambda632.8nm,z2m满足条件512 * (1e-2/512)^2 / 632.8e-9 31.6m完全够用。还要注意束腰和物理尺寸的匹配如果w01mm而物理范围L1cm高斯光和涡旋环占整个视场的比例适中如果L取得太小光场会被截断衍射图样出现虚假的方形边缘条纹。5.3 避坑指南从仿真到实验的映射最后分享几个在仿真和实验对照中发现的重要经验相位屏的螺旋方向要和实验一致。SLM加载的相位图通常按8位灰度图编码相位2π对应灰度255。Matlab里exp(1i*angle)和mod(angle, 2*pi)要对应正确否则仿真的螺旋方向跟实验相反。不用刻意把网格取太大。512×512在多数计算机上秒出结果1024×1024也行超过2048×2048速度会明显下降而且对角谱法来说大矩阵的FFT内存开销很高有些老机器会直接内存溢出。角谱法传播多段距离时网格尺寸始终保持不变。这正是角谱法相对单次菲涅尔积分的优点。若用菲涅尔衍射积分传播后网格采样间隔会改变给多屏湍流仿真带来额外麻烦。所以我强烈推荐统一用角谱法。仿真归一化很重要。实验室激光器功率恒定仿真中如果不做总功率归一化传播前后能量可能因数值误差而变化。用我代码里的归一化方式可以确保不同拓扑荷、不同传播距离之间光场达到可比性。对比不同拓扑荷时务必保持初始光场总功率一致。涡旋光束的中心暗核随拓扑荷增大而增大如果初始功率不一致会误判为传播损耗。归一化之后这个坑就避免了。我个人在多次仿真中的体会是涡旋光束仿真的最大价值不是复现教科书结果而是让你建立对相位结构、传播行为和检测方法的直觉。比如拓扑荷越高中心暗核越大传播中受湍流影响越敏感——这些规律看论文印象不深自己跑一遍透射强度图和相位图马上就记住了。后续想深入做可以从两个方向扩展一是把湍流相位屏换成时间动态序列模拟大气闪烁对OAM复用系统误码率的影响二是把深度学习解调加进来利用Matlab的Deep Learning Toolbox训练一个CNN直接从畸变强度图中识别拓扑荷。这些都是在当前代码基础上的自然演进做通之后涡旋光束从生成到检测的整个链路就全部掌握了。