轴棱锥MATLAB仿真：贝塞尔高斯光束与环形光束的相位与网格采样实战

这阵子一直在调一个轴棱锥的MATLAB仿真目标很明确用轴棱锥产生贝塞尔高斯光束同时还要拿到环形光束的光强图像。听起来就是两个光强图的事实际跑起来我却发现文献里一行公式背后藏着不少网格、采样和相位处理上的坑。如果你正在做光学仿真或者准备把轴棱锥写进课程设计、课题验证里这篇文章应该能帮你少走几趟弯路。我会从轴棱锥的相位模型讲起把贝塞尔高斯光束和环形光束的物理来源说明白再给出一份可以直接运行的MATLAB仿真代码最后聊一聊我在排查“光强图不对”时积累的几条经验。文章里的参数都给出具体数值和单位照着改就能复现。1. 轴棱锥的光学“性格”锥角、相位与光束的底层逻辑1.1 轴棱锥为什么能把高斯光束变成“无衍射光束”轴棱锥axicon本质上是一个绕光轴旋转对称的圆锥面光学元件。普通透镜是球面波前被掰成一个球面轴棱锥则是把入射平面波的波前掰成一个锥面。用数学语言讲轴棱锥的相位调制近似为φ(r) -k (n - 1) α r其中 k 2π/λ 是波数n 是折射率α 是轴棱锥的底角r 是离光轴的径向距离。这个公式说明轴棱锥的相位不是透镜那样的二次项而是线性项。也就是说r 越大相位滞后越多等相位面从平面被掰成圆锥面。为什么圆锥面波前能形成贝塞尔光束可以这样直观理解锥面波可以看作无数个方向一致、绕光轴排列的平面波分量每个分量的波矢都与光轴保持同一个夹角 θ ≈ (n-1)α。这些平面波在锥面重叠区域发生干涉干涉图样在横截面上满足零阶贝塞尔函数 J0(k_t r)在纵向上则呈现“无衍射”特性。换句话说轴棱锥把一个有限尺寸的高斯光束转化为一个在有限距离内近似无衍射的贝塞尔光束。1.2 无衍射距离和中心光斑尺寸先算一笔账仿真之前先把两个关键参数算出来中心主瓣半径 r0 和无衍射距离 z_max。近似公式是横向波矢 k_t ≈ k (n - 1) α中心主瓣半径 r0 ≈ 2.405 / k_t无衍射距离 z_max ≈ w0 / (n - 1) α我用的典型参数是入射高斯光束束腰 w0 1 mm波长 λ 0.6328 μm轴棱锥折射率 n 1.5底角 α 1°0.01745 rad。计算得到k 9.93 × 10^6 m^-1k_t ≈ 8.66 × 10^4 m^-1r0 ≈ 27.8 μmz_max ≈ 0.115 m也就是说中心亮斑直径大约55 μm无衍射区大约从轴棱锥后方一直延伸到11.5 cm。这个数字决定了仿真该在哪个位置观察贝塞尔高斯光束也决定了网格范围和采样点数。很多同学一上来就把传播距离设到1 m出来的图当然是一圈一圈的大环那不是贝塞尔高斯光束而是远场环形图样。1.3 为什么仿真里通常把轴棱锥当纯相位屏实际轴棱锥有厚度有反射、吸收和色散但在单色光、近轴近似下仿真里最常用的做法是把它当成一个纯相位屏只改相位不改振幅。也就是入射复振幅 E0(r) 乘以透过率函数 t(r) exp(-i k (n-1) α r)就完成了轴棱锥的建模。这个做法的前提是锥角足够小相位近似成立同时入射光束要完全覆盖在轴棱锥的有效口径内否则边缘衍射会非常明显。我在仿真里对比过带真实口径的模型发现纯相位屏对光强分布的主特征几乎没有影响却能省掉大量几何建模的麻烦。所以第一次跑通、看趋势用纯相位屏完全够用。2. 贝塞尔高斯光束与环形光束同源的两个面孔2.1 理想贝塞尔光束只是数学理想实际得到的是贝塞尔高斯理想的无衍射光束在数学上由零阶贝塞尔函数描述E(r) ∝ J0(k_t r)这个场在横向始终不扩散能量无限显然无法在现实里实现。轴棱锥加高斯入射光的实际产物学界叫贝塞尔高斯光束。它的表达式可以近似看成高斯包络乘以贝塞尔函数中心亮、旁瓣弱、总能量有限。所以在轴棱锥后的无衍射距离内观察横截面会看到中心一个特别亮的小圆斑周围是几圈逐渐变暗的同心环这就是贝塞尔高斯光束光强图。如果切一条经过中心的横向强度曲线理想贝塞尔光束的第一零点在 2.405/k_t 处贝塞尔高斯光束的第一零点位置会略微偏移但基本接近中心亮斑半径 28 μm。这个特征很适合做验证中心亮斑越小说明锥角越大(n-1)α 越大。2.2 环形光束从哪里来远场和焦平面上的“同一回事”在轴棱锥的无衍射区里光轴上是亮点但如果继续往后传播越过 z_max 之后中心强度会逐渐下降光强分布慢慢变成一个亮环。这个亮环的中心是暗的也就是环形光束。为什么会出现这个变化因为轴棱锥的本质作用是给不同横向分量赋予线性相位在远场或透镜焦平面这些分量被重新摊开形成与锥角严格对应的一个环。环形光束的半径近似为自由传播远场r_ring ≈ z (n - 1) α透镜焦平面r_ring ≈ f (n - 1) α拿 0.6 m 的自由传播距离算环半径大约 0.6 × 0.00873 ≈ 5.2 mm拿焦距 100 mm 的透镜算焦平面环半径大约是 0.873 mm。只要观察位置足够远轴棱锥输出的光场就会呈现环形光强图像。理解这一点很重要贝塞尔高斯光束和环形光束不是两套互不相干的物理现象而是在不同纵向位置看同一个锥面波前的结果。2.3 仿真任务里两种光强图像的正确对应既然标题要求同时得到两类光强图像我的做法是把同一份场分别传播到两个距离z_bessel 50 mm位于 z_max ≈ 115 mm 内横截面为贝塞尔高斯光束z_ring 600 mm远大于 z_max横截面为环形光束只改传播距离不改其他参数就能在一套代码里拿到两张图。如果你真的拿透镜做实验环形光束对应的是透镜焦平面的光强图仿真里也可以用傅里叶变换实现但角谱传播到足够远的效果更直观也更容易和实验里挪动观察屏的习惯对上。3. 仿真参数标定与网格采样很多错误从这里开始3.1 一张参数清单先照着设我建议初学者先照抄下面这组参数跑通之后再根据需求改。所有单位都统一成国际单位制否则 MATLAB 里写错单位是最常见的翻车现场。参数符号数值波长λ0.6328e-6 m入射束腰w01e-3 m轴棱锥折射率n1.5底角α1° 0.0174533 rad网格半宽L10e-3 m采样数N1024观察距离1z_bessel0.05 m观察距离2z_ring0.6 m为什么网格半宽取 10 mm因为传播到 0.6 m 之后环形光束的半径大约 5.2 mm加上高斯包络本身的尺寸10 mm 的半宽不会让光场碰到边界。如果网格半宽太小角谱传播会因截断产生严重边缘衍射条纹。3.2 相位屏的奈奎斯特条件一个容易被忽略的硬约束轴棱锥相位是 φ(r) -k γ r其中 γ (n-1)α。径向方向相邻两个采样点之间的相位变化为Δφ k γ Δr要正确表示这个快速变化的相位必须满足 Δφ π也就是 Δr π / (k γ)。对这组参数k γ ≈ 8.66×10^4 m^-1所以 Δr 必须小于约 36 μm。网格半宽 L 10 mm、N 1024 时Δr 2L/N ≈ 19.5 μm满足条件。如果你把 L 直接改成 20 mm 而不增大 NΔr 变成 39 μm超出约束光强图上会出现一圈一圈的伪周期条纹而且不是贝塞尔旁瓣。我最初就在这上面栽过跟头后来把 N 加到 2048或者减小 L 才解决。这个计算大概是整个仿真里最容易被忽略、又最影响图像正确性的地方。3.3 角谱传播的频率坐标安排自由空间角谱传播的传递函数是H(f_x, f_y) exp[ i k z sqrt(1 - λ²(f_x² f_y²)) ]MATLAB 里做这一步时最需要注意的是频率坐标和 fftshift/ifftshift 的顺序。常见错误是用 linspace 生成频率坐标后不配合 fftshift结果图像整体平移、相位混乱。我的做法是空间坐标用 x (-N/2:N/2-1) * dxdx 2L/N频率坐标用 fx (-N/2:N/2-1) / (N*dx)对输入场先 fftshift再 fft2再乘以 H再 ifft2 ifftshift 还原这个顺序一旦乱了出的图要么中心对不上要么旋转了半个格子。后面我会在代码里按这个顺序写死直接复制就行。4. MATLAB代码实现一次完整的轴棱锥光束仿真4.1 主程序从高斯光到轴棱锥相位下面这段代码在大多数 MATLAB 版本上都能跑核心就是基础函数和一次二维 FFT不需要额外工具箱。clear; clc; close all; % 物理参数 lambda 0.6328e-6; k 2*pi/lambda; w0 1e-3; n_ax 1.5; alpha 1*pi/180; gamma (n_ax - 1) * alpha; % 网格参数 N 1024; L 10e-3; dx 2*L/N; x (-N/2:N/2-1) * dx; [X, Y] meshgrid(x, x); R sqrt(X.^2 Y.^2); % 入射高斯光束 E0 exp(-R.^2 / w0^2); % 轴棱锥相位屏 phase_axicon -k * gamma * R; T_axicon exp(1i * phase_axicon); % 透过轴棱锥后的复振幅 E1 E0 .* T_axicon; % 观察贝塞尔高斯光束无衍射区内 z_bessel 0.05; E2_bessel ang_spectrum(E1, z_bessel, lambda, dx); % 观察环形光束远场 z_ring 0.60; E2_ring ang_spectrum(E1, z_ring, lambda, dx); % 绘图 figure(Name,Bessel-Gauss Beam Intensity); I_bessel abs(E2_bessel).^2; imagesc(x*1e3, x*1e3, I_bessel); axis image; axis xy; colormap(hot); colorbar; xlabel(x / mm); ylabel(y / mm); title(Bessel-Gauss Beam Intensity at z50mm); figure(Name,Ring Beam Intensity); I_ring abs(E2_ring).^2; imagesc(x*1e3, x*1e3, I_ring); axis image; axis xy; colormap(hot); colorbar; xlabel(x / mm); ylabel(y / mm); title(Ring Beam Intensity at z600mm); % 中心横向截面曲线 figure(Name,Cross-section Curves); subplot(1,2,1); Ic I_bessel(N/21, :); plot(x*1e3, Ic/max(Ic), LineWidth, 1.2); xlim([-0.1, 0.1]); xlabel(x / mm); ylabel(Normalized Intensity); title(Bessel-Gauss center section); subplot(1,2,2); Ir I_ring(N/21, :); plot(x*1e3, Ir/max(Ir), LineWidth, 1.2); xlabel(x / mm); ylabel(Normalized Intensity); title(Ring beam center section);4.2 角谱传播函数一套稳定写法角谱传播在多个位置都要调用所以我单独写成一个函数放在同目录下的 ang_spectrum.m 里。function E_out ang_spectrum(E_in, z, lambda, dx) % 自由空间角谱传播 % E_in : 输入复振幅 % z : 传播距离 % lambda: 波长 % dx : 空间采样间隔 N size(E_in, 1); fx (-N/2:N/2-1) / (N*dx); [FX, FY] meshgrid(fx, fx); k 2 * pi / lambda; f2 FX.^2 FY.^2; % 传递函数 H exp(1i * k * z * sqrt(1 - (lambda^2) * f2)); % 倏逝波分量置零避免数值发散 H(lambda * sqrt(f2) 1) 0; % 频域运算 U_in fftshift(fft2(E_in)); U_out U_in .* H; E_out ifft2(ifftshift(U_out)); end这里有几个细节值得说。第一H 中 sqrt 里的负数一定要处理否则 exp(实部×复数) 会让结果指数爆炸。频域里超过 1/λ 的分量属于倏逝波不参与传播直接置零是最稳妥的做法。第二整个 fft2、fftshift、ifftshift 的顺序要保持一致坐标都定义成从负到正的对称顺序这是经过多次比对确认没问题的写法。4.3 运行结果应该长什么样跑完之后第一张图在 z50 mm 处应该是中心一个直径约55 μm的亮斑外围几圈同心条纹条纹间距越靠外越密、亮度越低这就是贝塞尔高斯光束。横向截面曲线在主瓣两侧有明显的第一、第二旁瓣但幅度衰减非常快。第二张图在 z600 mm 处应该是一个半径约5.2 mm的亮环中心强度接近零环本身有一定宽度这就是轴棱锥的远场环形光束。如果第一张图中心没有亮斑或者第二张图中心不是暗的先别急着改物理参数大概率是网格或采样的问题。下一节专门讲排查链路。5. 结果不对时的排查链路光强图常见的几种“翻车现场”5.1 中心不是亮斑而是彩条或条纹错位这种症状最常见的原因是 fftshift 和 ifftshift 的顺序写反了或者频率坐标与 MATLAB 的 FFT 索引没对齐。排查方法很简单先不传播直接画 E1 E0 .* T_axicon 的强度图。如果这步就不对说明相位屏或网格有问题如果这步对说明问题出在角谱传播函数里。我个人的建议是不要在同一个脚本里来回试把角谱传播封装成函数固定输入输出格式排查起来边界清晰。5.2 环形半径和理论值对不上自由传播远场环形光斑的中心半径理论上是 z * gamma。拿 0.6 m 距离算r_ring ≈ 5.2 mm。如果你算出来的半径明显偏大或偏小先检查单位角度有没有从度换成弧度波长和束腰是不是都转成了米。这是我在各类仿真问题里见过最多的参数单位问题。另一个容易出问题的点是网格半宽太小环被边界切掉视觉上半径看起来就不对。5.3 图像出现辐射状条纹或同心硬边辐射状条纹一般对应奈奎斯特条件不满足也就是 Δr 太大相位折叠。回到第3.2节的公式把 N 增大或把 L 减小。同心硬边则多半是边界截断造成的入射光或传播后的光碰到了网格边缘导致矩形孔径衍射。适当加大 L 或者加密网格可以解决。5.4 传播距离一长就发散或者全乱角谱传递函数在使用大 z 值时对网格采样要求更高。如果你的 z 从 0.05 m 改成 0.6 m图像突然出现大量高频伪影建议先检查 H 中倏逝波处理再看频率范围是否覆盖得足够。脚本里已经加了 H(lambda * sqrt(f2) 1) 0 这行如果删掉了长距离传播很快会出问题。此外远场观察要保证光场没有扩散到网格边界必要时把 L 同步增大同时按奈奎斯特约束增大 N。5.5 光强图灰度范围不对看不清条纹贝塞尔高斯光束的中心亮斑强度远大于旁瓣直接 imagesc 默认色标会让旁瓣全部看不见。解决的办法是画横向截面曲线或者在作图时对光强做对数变换比如 log(1 I/max(I)) 拉伸。我平时喜欢把二维强度图和横向截面曲线放在同一张图里两个视角互补信息完整得多。6. 从仿真到实验几个值得提前知道的事6.1 元件参数与标称值的偏差仿真跑通之后如果后面要进实验室有几条经验值得提前写出来。轴棱锥的锥角精度非常敏感。底角 1° 的轴棱锥实际加工误差哪怕只有 0.01°小到肉眼完全看不出来中心光斑半径也会产生百分之几的变化。买元件的时候一定要跟商家要实测数据不要只信标称值。入射光束的束腰要调准。贝塞尔高斯光束中心亮斑的大小几乎由锥角单独决定但无衍射距离 z_max 正比于入射束腰 w0。如果实验室的准直光斑实际只有 0.5 mm无衍射距离会缩到 57 mm观察位置必须相应调整。6.2 对准、饱和与测量细节对准是轴棱锥实验最大的坑。轴棱锥顶点只要横向偏半个毫米光轴就会偏移中心亮斑可能直接被切掉一块。仿真里默认严格共轴实验里要做到这一点往往要花很多时间。用 CCD 拍摄环形光束时切记不要让感光面饱和。贝塞尔高斯光束的中心峰值能量密度远高于平均光强低功率激光也可能打出一个亮到发白的中心点看起来很好看但量化分析就全废了。这些点说起来都不复杂但实际操作里优先级往往比调代码更高。仿真能告诉你理想情况下的光强分布实验里多出来的每一项误差都会在二维图像上留下痕迹。我自己的习惯是每改一次参数就重新算一遍第3.2节的 Δφ 约束并顺手把观察距离除以 z_max 打印出来确认自己到底是在无衍射区、过渡区还是远场。仿真脚本里我会同时保留 z_bessel 和 z_ring 两个距离跑完对照理论半径再用横向截面确认中心暗/亮状态。这套流程看起来笨但确实帮我挡掉过很多次“以为结果错了其实是采样错了”的情况。希望这份代码和排查经验能帮你在轴棱锥的两张光强图上少走几趟弯路。