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) )。
对于涡旋光束仿真，大多数情况下取 ( p=0 ) 就已经能体现全部涡旋特性，因为拉盖尔多项式 ( 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); % 波前曲率半径 (z=0时无穷大,代码里特判) zeta = atan2(z, zR); % Gouy相位 % 拉盖尔多项式值, 广义版用 laguerreL 需要符号工具箱, 这里手写定义 % 对于p=0和p=1最常用: if p == 0 Lag = ones(size(R)); elseif p == 1 Lag = 1 + abs(l) - 2 * R.^2 / w^2; else % 一般场景p<=3, 用递推关系; 更多阶可定义递归函数 Lag = laguerre_gen(p, abs(l), 2*R.^2/w^2); end % 径向归一化因子 C = sqrt( 2*factorial(p) / (pi * factorial(p+abs(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*p+abs(l)+1)*zeta); % 处理z=0处半径无穷大的问题: 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和径向距离R；z=0时波前曲率半径无穷大，需要特殊处理；归一化通过总功率完成，保证不同参数下仿真结果可以横向比较。
2.2 用螺旋相位屏快速生成各类涡旋光束
实际工程里不总需要完整LG解，很多场景只是想在已有光场（比如平面波或高斯光）上套一个涡旋相位，最典型的做法就是用螺旋相位屏（Spiral Phase Plate，SPP）的相位掩模：
[ 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); % 叠加上拓扑荷l=2的螺旋相位 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)');运行后你会看到：强度图是典型的甜甜圈形状，中心暗核清晰；相位图则是围绕中心从蓝色到黄色渐变一圈的螺旋阶梯。拓扑荷 ( l=1 ) 时相位刚好转一圈，相位图上有1个2π跳跃；( l=2 ) 时转动两圈。相位图中奇点位置就是强度暗核中心。
从这个结果可以引申出一个仿真中的重要经验：判断涡旋光束生成对不对，不要只看强度，一定要看相位图。强度中空结构可能是高阶模式或失焦造成的，只有相位图呈现干净的螺旋状，才能确认光束确实携带轨道角动量。
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(f==0) = 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 )。也就是说网格越细、波长越长，能传播的距离上限越大。如果传播距离超出该值，仿真结果会出现严重的边缘折叠误差。我给出的建议是先用一个自己熟悉的参数组合跑通，再逐步调整。比如N=512,L=1cm,lambda=632.8nm,z=2m满足条件：512 * (1e-2/512)^2 / 632.8e-9 = 31.6m，完全够用。
还要注意束腰和物理尺寸的匹配：如果w0=1mm而物理范围L=1cm，高斯光和涡旋环占整个视场的比例适中；如果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直接从畸变强度图中识别拓扑荷。这些都是在当前代码基础上的自然演进，做通之后，涡旋光束从生成到检测的整个链路就全部掌握了。
