1. 项目概述：从理论公式到可视化光束
高斯光束，这大概是光学和激光领域里最基础也最重要的概念之一了。但凡接触过激光原理、光纤通信或者光学设计，都绕不开它。但说实话，光看教科书上那一堆关于束腰、瑞利长度、发散角的公式，总感觉隔着一层纱，不够直观。我记得自己刚开始学的时候，就特别想“看见”这个光束到底长什么样，它在空间中是怎么传播和演变的。
这就是为什么我觉得用MATLAB来做高斯光束的仿真特别有价值。它不是一个复杂的科研项目，而是一个将抽象理论具象化的绝佳工具。通过几行代码，你就能亲手“创造”出一束光，观察它的强度分布，看着它从束腰处最细的状态慢慢发散开来，所有公式里的参数都变成了屏幕上可以调节的滑块和可以直观看到的图像变化。这个过程不仅能帮你牢牢记住那些关键概念，更能培养一种“计算光学”的思维——用编程和数值计算来解决和验证光学问题。
这篇内容，我就想和你一起，从最基础的高斯光束公式出发，一步步在MATLAB里把它仿真出来。我们会涵盖一维、二维的强度分布，三维的空间传播，以及一些常用的分析技巧。无论你是正在学习相关课程的学生，还是需要快速验证光学设计的工程师，这个简单直接的仿真过程都能给你带来实实在在的帮助。
2. 高斯光束理论基础与MATLAB实现思路
在动手写代码之前，我们得先统一一下“语言”，明确我们要仿真的对象到底是什么。高斯光束，顾名思义，其横截面上的光强分布遵循高斯函数。它的数学描述是解析的，这为我们的仿真提供了清晰的路径。
2.1 核心公式拆解
我们通常从最基本的傍轴近似下的标量波方程解说起。一个沿z轴传播的高斯光束，其复振幅可以表示为：
U(x, y, z) = A0 * (w0 / w(z)) * exp(-(x^2 + y^2) / w(z)^2) * exp(-i*k*z - i*k*(x^2+y^2)/(2*R(z)) + i*ζ(z))
这个公式看起来有点唬人，但我们把它拆开看，每个部分都有明确的物理意义：
振幅部分：
A0 * (w0 / w(z)) * exp(-(x^2 + y^2) / w(z)^2)A0是常数，代表振幅。w0是束腰半径，即光束最细处的半径（光强下降到中心值的1/e^2处的半径）。w(z)是z位置处的光束半径，它随传播距离变化：w(z) = w0 * sqrt(1 + (z / zR)^2)。这里引入了瑞利长度zR = π * w0^2 / λ，它是衡量光束准直程度的重要参数，可以理解为光束大致保持平行的传播距离。
相位部分：
exp(-i*k*z - i*k*(x^2+y^2)/(2*R(z)) + i*ζ(z))exp(-i*k*z)是平面波的相位延迟。exp(-i*k*(x^2+y^2)/(2*R(z)))是波前曲率带来的相位因子，R(z)是波前曲率半径，R(z) = z * (1 + (zR / z)^2)。在束腰处(z=0)，R(z)为无穷大，波前是平面；远离束腰，波前变为球面。exp(i*ζ(z))是古依相位，ζ(z) = atan(z / zR)，这是一个额外的相位移动，在光束通过焦点时会发生π的相位跃变。
我们最常关心的是光强分布I(x, y, z) = |U(x, y, z)|^2。对于许多仿真，我们可以暂时忽略相位部分，专注于强度仿真，公式简化为：I(x, y, z) = I0 * (w0 / w(z))^2 * exp(-2*(x^2 + y^2) / w(z)^2)其中I0是束腰中心处的峰值光强。
2.2 MATLAB仿真策略设计
有了公式，如何在MATLAB中构建它？我的思路是分层次、模块化地进行：
- 参数化定义：将所有物理参数（
λ,w0,zR等）设为变量。这样，我们只需要改变几个参数，就能立刻看到不同波长、不同束腰大小的光束行为，仿真就活了。 - 网格化空间：使用
meshgrid函数创建代表x, y, z坐标的数值网格。这是数值计算的基础，我们的光束将在这个网格上被“计算”出来。 - 向量化计算：利用MATLAB强大的矩阵运算能力，避免使用低效的循环。直接对整个坐标网格应用上述公式，一次性计算出整个空间的光场或光强分布。这是MATLAB仿真的性能关键。
- 可视化呈现：这是点睛之笔。我们将使用：
imagesc或pcolor来显示二维截面图。surf或mesh来绘制三维强度分布。plot来绘制一维光强剖面，用于精确测量束宽。subplot将多角度视图组合在一起，形成全面的分析面板。
实操心得：在开始写代码前，最好在纸上或注释里把公式和计算步骤列清楚。尤其是
w(z)和R(z)这些z的函数，确保你在计算每个z平面时都用对了对应的值。一个常见的错误是在计算三维体数据时，错误地广播了数组维度，导致结果不对。
3. 从二维截面到三维传播：完整仿真实现
理论清晰了，策略定好了，现在打开MATLAB，我们开始动手实现。我会从最简单的静态二维截面开始，逐步扩展到完整的三维空间传播仿真。
3.1 基础二维束腰光强分布仿真
我们从仿真的起点——束腰位置（z=0）的二维光强分布开始。这是最直观的一步。
%% 1. 基础参数设置 clear; close all; clc; % 清空环境，好习惯 lambda = 632.8e-9; % 波长，单位米，这里用常见的He-Ne激光波长 w0 = 1e-3; % 束腰半径，1毫米 I0 = 1; % 中心峰值光强，归一化为1 % 设置观察平面范围和采样点数 x = linspace(-3*w0, 3*w0, 500); % x轴范围：-3倍束腰到+3倍束腰 y = linspace(-3*w0, 3*w0, 500); % y轴范围 [X, Y] = meshgrid(x, y); % 生成二维网格坐标 %% 2. 计算束腰处光强分布 (z=0, 此时 w(z)=w0) % 二维高斯分布公式 I_xy = I0 * exp(-2*(X.^2 + Y.^2) / w0^2); %% 3. 可视化 figure('Position', [100, 100, 1200, 400]); % 设置一个宽幅图窗 % 子图1：二维伪彩色图 subplot(1, 3, 1); imagesc(x*1e3, y*1e3, I_xy); % 坐标转换为毫米显示 axis image; % 保持纵横比相等 xlabel('x (mm)'); ylabel('y (mm)'); title('束腰处光强分布 (二维视图)'); colorbar; colormap('hot'); % 使用‘hot’色图模拟光强 % 子图2：三维曲面图 subplot(1, 3, 2); surf(X*1e3, Y*1e3, I_xy, 'EdgeColor', 'none'); xlabel('x (mm)'); ylabel('y (mm)'); zlabel('相对光强'); title('束腰处光强分布 (三维视图)'); colormap('hot'); view(30, 30); % 调整视角 % 子图3：通过中心的一维剖面 subplot(1, 3, 3); plot(x*1e3, I_xy(ceil(end/2), :), 'LineWidth', 2); % 取中间一行数据 xlabel('x (mm)'); ylabel('相对光强'); title('束腰处水平中心光强剖面'); grid on; hold on; % 标记1/e^2强度点（对应束腰半径w0） plot([-w0, w0]*1e3, [exp(-2), exp(-2)], 'r--', 'LineWidth', 1.5); legend('光强分布', '1/e^2 阈值', 'Location', 'best');这段代码运行后，你会得到一个包含三个视图的图形。二维图让你看到经典的高斯光斑，三维图让你感受强度的起伏，而一维剖面图则是定量分析的关键。你可以从剖面图上精确读出光强下降到中心值1/e^2（约0.135）的位置，那就是±w0。
注意事项：
linspace生成的坐标点数量决定了图像的分辨率和计算量。500个点对于演示是清晰的，但如果你的计算区域很大或者需要非常精细，可以增加到1000或更多。同时，注意meshgrid生成的X和Y是矩阵，这确保了后续的点乘(.^2)和矩阵运算能正确进行。
3.2 三维空间传播仿真
静态截面看完了，我们让光束“动”起来，看看它怎么传播。这意味着我们要计算一系列不同z位置上的二维光强分布，然后堆叠或选取视图来观察。
%% 1. 扩展参数：定义传播轴 z = linspace(-3*zR, 3*zR, 150); % 传播距离范围：-3倍瑞利长度到+3倍瑞利长度 % 计算瑞利长度 zR = pi * w0^2 / lambda; fprintf('束腰半径 w0 = %.2f mm, 瑞利长度 zR = %.2f mm\n', w0*1e3, zR*1e3); % 创建三维空间网格（注意：全三维网格数据量巨大，这里采用循环计算每个z面） % 我们依然使用之前的x-y网格 I_3D = zeros(length(y), length(x), length(z)); % 预分配内存，提升效率 %% 2. 循环计算每个z平面的光强分布 for k = 1:length(z) z_current = z(k); % 计算当前z处的光束半径w(z)和波前曲率半径R(z)（此处仅用于强度，R(z)未用） w_z = w0 * sqrt(1 + (z_current / zR)^2); % 计算当前z平面的光强分布 I_3D(:, :, k) = I0 * (w0 / w_z)^2 .* exp(-2*(X.^2 + Y.^2) / w_z^2); end %% 3. 三维传播可视化 figure('Position', [100, 100, 1400, 500]); % 子图1：X-Z平面（y=0）的切片图 subplot(1, 3, 1); % 提取y=0中心线上的数据 [~, idx_y0] = min(abs(y - 0)); % 找到y坐标最接近0的索引 I_xz = squeeze(I_3D(idx_y0, :, :)); % 挤压维度，得到(x, z)矩阵 imagesc(z*1e3, x*1e3, I_xz'); xlabel('传播距离 z (mm)'); ylabel('x (mm)'); title('X-Z平面 (y=0) 光强分布'); axis xy; % 确保y轴方向正确 colorbar; colormap('hot'); % 标记束腰位置(z=0)和瑞利长度位置 hold on; plot([0,0], [x(1), x(end)]*1e3, 'w--', 'LineWidth', 1); plot([-zR, zR]*1e3, [0,0], 'c--', 'LineWidth', 1); legend('', '束腰位置 z=0', '瑞利长度 ±zR', 'Location', 'best'); % 子图2：不同z位置的光束半径变化曲线 subplot(1, 3, 2); w_z_array = w0 * sqrt(1 + (z / zR).^2); % 理论值 plot(z*1e3, w_z_array*1e3, 'b-', 'LineWidth', 2); xlabel('传播距离 z (mm)'); ylabel('光束半径 w(z) (mm)'); title('光束半径随传播距离变化'); grid on; hold on; % 标记束腰和渐近线 plot([z(1), z(end)]*1e3, [w0, w0]*1e3, 'r--', 'LineWidth', 1); plot([z(1), z(end)]*1e3, [w0/zR * abs(z(end)), w0/zR * abs(z(end))]*1e3, 'g--', 'LineWidth', 1); legend('w(z) 理论曲线', '束腰半径 w0', '渐近线 (发散角)', 'Location', 'best'); % 子图3：选取三个特征z位置，显示其二维光斑 subplot(1, 3, 3); z_plots = [-zR, 0, zR]; % 选择 z = -zR, 0, +zR 三个面 plot_titles = {'z = -zR (束腰上游)', 'z = 0 (束腰处)', 'z = +zR (束腰下游)'}; for p = 1:3 [~, idx_z] = min(abs(z - z_plots(p))); % 找到最接近的索引 I_slice = I_3D(:, :, idx_z); % 创建一个临时坐标（为了在子图中显示） ax = subplot(3, 3, 6+p); % 定位到右列的下三个子图位置 imagesc(ax, x*1e3, y*1e3, I_slice); axis(ax, 'image'); xlabel(ax, 'x (mm)'); ylabel(ax, 'y (mm)'); title(ax, plot_titles{p}); colormap(ax, 'hot'); if p == 3 colorbar('eastoutside'); end end这段代码生成了更丰富的视图。X-Z平面切片图清晰地展示了光束从汇聚到束腰再到发散的全过程，你能看到光束宽度如何平滑变化。光束半径变化曲线则定量验证了w(z)的双曲线公式。最后对比三个特征位置的光斑，直观展示了“束腰处最细，离束腰越远光斑越大”的规律。
实操心得：计算三维数据
I_3D时，预分配数组（zeros）至关重要。如果是在循环内动态扩展数组，当z的点数很多时，MATLAB会反复申请内存，导致速度极慢甚至内存不足。另外，对于非常大的三维网格，直接存储所有数据可能内存吃不消。这时可以采用“按需计算”的策略，即只在需要画图的时候计算当前切片的数据，而不是一次性计算全部。
4. 进阶分析与常见问题调试
基础仿真跑通后，我们可以玩点更深入的，同时也会遇到一些典型问题。这部分是区分“会写代码”和“理解仿真”的关键。
4.1 相位信息仿真与波前可视化
之前我们只关注了强度。要完整描述光场，相位必不可少，尤其是在涉及干涉、衍射或光束整形时。
%% 1. 计算包含相位的复振幅场 (在束腰处z=0) z_target = 0; % 选择束腰位置 k = 2 * pi / lambda; % 波数 % 注意：在束腰处，R(z) -> inf，因此波前曲率相位项为0。古依相位 ζ(z)=0。 U_xy = sqrt(I0) * exp(-(X.^2 + Y.^2) / w0^2); % 复振幅，此时相位均匀 %% 2. 计算离开束腰后的复振幅场 (例如 z = zR) z_target = zR; w_z = w0 * sqrt(1 + (z_target / zR)^2); R_z = z_target * (1 + (zR / z_target)^2); % 波前曲率半径 zeta_z = atan(z_target / zR); % 古依相位 % 完整的复振幅计算 U_xy_z = sqrt(I0) * (w0 / w_z) .* exp(-(X.^2 + Y.^2) / w_z^2) .* ... exp(-1i * k * z_target) .* ... exp(-1i * k * (X.^2 + Y.^2) / (2 * R_z)) .* ... exp(1i * zeta_z); %% 3. 可视化相位 phase_z = angle(U_xy_z); % 提取相位，范围 [-π, π] figure; subplot(1,2,1); imagesc(x*1e3, y*1e3, phase_z); axis image; colorbar; xlabel('x (mm)'); ylabel('y (mm)'); title(sprintf('z = %.1f mm 处的相位分布', z_target*1e3)); colormap('hsv'); % HSV色图非常适合表示周期性相位 subplot(1,2,2); % 绘制通过中心的相位剖面 plot(x*1e3, phase_z(ceil(end/2), :), 'LineWidth', 2); xlabel('x (mm)'); ylabel('相位 (弧度)'); title('相位中心剖面'); grid on; % 可以看到相位随x^2变化的抛物线形状，这正是球面波前的特征。4.2 常见仿真问题与排查技巧
在实际操作中，你可能会遇到下面这些问题。这里是我的排查清单：
| 问题现象 | 可能原因 | 排查与解决方法 |
|---|---|---|
| 图形显示一片空白或全黑/全白 | 1. 数据值超出显示范围或全部为0/NaN。 2. imagesc的坐标轴设置错误，数据画到看不见的地方。3. 色图范围不合适。 | 1. 在imagesc后使用caxis([0, 1])手动设置颜色轴范围。或用max(I_xy(:))检查数据最大值。2. 检查 imagesc(x, y, I)中的x, y向量是否与矩阵I的维度匹配（size(I,2)==length(x),size(I,1)==length(y)）。3. 尝试 colormap(gray)或imagesc(..., 'CDataMapping', 'scaled')。 |
| 三维曲面图非常粗糙或畸变 | 网格点太少 (linspace点数不足)。 | 增加linspace的采样点数，例如从100增加到500。对于曲面图，也可以使用shading interp命令进行平滑插值。 |
| 计算速度极慢 | 1. 使用了多重嵌套for循环。2. 未预分配大数组。 3. 网格分辨率过高。 | 1.务必使用向量化操作。将exp,sqrt等函数直接作用于整个矩阵X,Y。2. 对于像 I_3D这样的结果数组，使用zeros预分配。3. 在保证可视效果的前提下，适当降低采样点数。 |
| 光束形状不对称或奇怪 | 1. x和y的坐标范围或采样点数不一致。 2. 公式中输入了错误的坐标变量。 | 1. 确保linspace(-range, range, N)中的range和N在x和y方向一致，除非你特意需要各向异性的光束。2. 检查公式中是 X.^2 + Y.^2而不是x.^2 + y.^2。前者是矩阵，后者是向量。 |
| “发散角”与理论值对不上 | 1. 在图中测量发散角的方法不对。 2. w0或λ的单位不一致。 | 1. 理论发散角（远场）θ = λ / (π * w0)。在仿真中，可通过拟合远场`w(z) ≈ θ * |
| 相位图出现剧烈跳变（锯齿） | 相位包裹问题。真实相位变化超过[-π, π]范围，angle函数将其折叠到此区间。 | 这是正常现象，表示相位变化剧烈。若要分析连续的相位变化，可以使用unwrap函数（对一维数据效果好），但对于二维相位解包裹需要更复杂的算法，这属于另一个专题。可视化时用hsv色图即可清晰观察相位周期。 |
避坑技巧：在编写每个关键公式后，用一两行代码做个快速验证。比如计算完
w_z后，fprintf('在z=%.2f处，w_z=%.2e\n', z_target, w_z)，看看输出是否合理。对于光束半径，在束腰处应该等于w0，在z=zR处应该等于sqrt(2)*w0。这种即时的小验证能帮你快速定位是公式写错了还是参数单位错了。
5. 仿真扩展与应用实例
一个扎实的基础仿真框架建立后，你可以像搭积木一样，扩展出很多有趣且实用的应用。
5.1 模拟光束通过透镜的变换
这是非常经典的应用。一个薄透镜可以改变高斯光束的波前曲率，从而变换其束腰位置和大小。
%% 模拟透镜对高斯光束的变换 f = 0.1; % 透镜焦距，100mm z_from_lens = 0.2; % 透镜到输入光束束腰的距离，200mm % 假设入射光束参数（在透镜位置处） lambda = 632.8e-9; w0_in = 0.5e-3; % 入射光束束腰半径 zR_in = pi * w0_in^2 / lambda; % 计算入射光束在透镜位置的光斑大小w_in和波前曲率半径R_in w_in = w0_in * sqrt(1 + (z_from_lens / zR_in)^2); R_in = z_from_lens * (1 + (zR_in / z_from_lens)^2); % 应用ABCD定律（薄透镜变换） % 对于薄透镜，传输矩阵为 [1, 0; -1/f, 1] % 复数q参数满足 1/q2 = 1/q1 - 1/f % 其中 q = z + i*zR， 或者用更常用的形式：1/q = 1/R - i*λ/(π*w^2) % 计算入射光束的q参数 q1_inv = 1/R_in - 1i*lambda/(pi*w_in^2); % 经过透镜后的q参数 q2_inv = q1_inv - 1/f; % 从q2参数提取出射光束的R_out和w_out q2 = 1 / q2_inv; R_out = 1 / real(1/q2); w_out = sqrt(-lambda / (pi * imag(1/q2))); % 进一步，可以计算出射光束的新束腰位置（相对于透镜）和束腰大小 z_out = R_out / (1 + (lambda*R_out/(pi*w_out^2))^2); w0_out = w_out / sqrt(1 + (pi*w_out^2/(lambda*R_out))^2); fprintf('经过焦距%.3f m的透镜后：\n', f); fprintf(' 出射光束在透镜处的光斑半径: %.3f mm\n', w_out*1e3); fprintf(' 出射光束在透镜处的波前曲率半径: %.3f m\n', R_out); fprintf(' 新的束腰位置（相对于透镜）: %.3f m\n', z_out); fprintf(' 新的束腰半径: %.3f mm\n', w0_out*1e3);这个例子展示了如何将高斯光束的q参数与ABCD矩阵理论结合，用代码实现复杂光学系统的仿真。你可以改变焦距f和入射距离z_from_lens，观察输出光束参数如何变化，这对于光学系统设计中的光路计算非常有帮助。
5.2 生成动态传播GIF动画
静态图片虽然好，但动态图更能直观展示传播过程。我们可以用循环生成一系列图片，然后合成GIF。
%% 生成光束传播动态GIF filename = 'gaussian_beam_propagation.gif'; % 沿用之前的三维数据 I_3D， 或者重新计算一段更短的z轴用于动画 z_anim = linspace(-2*zR, 2*zR, 50); figure('Position', [100, 100, 800, 600], 'Color', 'white'); for n = 1:length(z_anim) % 计算或提取第n帧的数据 (这里假设重新计算) w_z = w0 * sqrt(1 + (z_anim(n) / zR)^2); I_frame = I0 * (w0 / w_z)^2 .* exp(-2*(X.^2 + Y.^2) / w_z^2); % 绘图 imagesc(x*1e3, y*1e3, I_frame); axis image; xlabel('x (mm)'); ylabel('y (mm)'); title(sprintf('高斯光束传播 | z = %.2f mm', z_anim(n)*1e3)); colormap('hot'); colorbar; caxis([0, 1]); % 固定颜色轴，让动画稳定 drawnow; % 捕获帧并写入GIF frame = getframe(gcf); im = frame2im(frame); [imind, cm] = rgb2ind(im, 256); if n == 1 imwrite(imind, cm, filename, 'gif', 'Loopcount', inf, 'DelayTime', 0.1); else imwrite(imind, cm, filename, 'gif', 'WriteMode', 'append', 'DelayTime', 0.1); end end fprintf('GIF动画已保存为: %s\n', filename);这个动画能非常生动地展示光束从汇聚到发散，束腰处最细的整个过程。生成的GIF可以用于报告、演示或教学材料中。
5.3 像散与椭圆高斯光束仿真
现实中的激光器可能输出非圆对称的，即椭圆高斯光束。它的两个正交方向（x和y）有不同的束腰和发散角。
%% 仿真椭圆高斯光束 w0x = 1.0e-3; % x方向束腰半径 w0y = 0.5e-3; % y方向束腰半径 zRx = pi * w0x^2 / lambda; zRy = pi * w0y^2 / lambda; % 计算特定z位置（例如z=0）的光强 z_astig = 0; wx_z = w0x * sqrt(1 + (z_astig / zRx)^2); wy_z = w0y * sqrt(1 + (z_astig / zRy)^2); % 光强分布为两个方向高斯函数的乘积 I_ellipse = I0 * exp(-2*(X.^2)/wx_z^2 - 2*(Y.^2)/wy_z^2); figure; imagesc(x*1e3, y*1e3, I_ellipse); axis image; xlabel('x (mm)'); ylabel('y (mm)'); title('椭圆高斯光束 (束腰处)'); colorbar; colormap('hot'); % 可以看到光斑在x和y方向宽度不同，呈椭圆形。通过调整w0x和w0y的比例，你可以仿真出从细线到扁椭圆的各种光束形状。更进一步，你还可以让两个方向的束腰位置不在同一z点，这就是像散光束的典型特征，仿真方法类似，只是分别计算x和y方向的光束半径w(z)时使用不同的z偏移量。
这些扩展实例只是抛砖引玉。基于这个核心框架，你还可以去仿真高阶模（如厄米-高斯或拉盖尔-高斯光束），模拟光束在湍流大气中的传播，或者与傅里叶光学结合计算光束的衍射。关键是理解物理模型，然后用MATLAB这个强大的计算工具把它实现出来。
