企业官网建设流程全解析

简介：一份围绕高斯光束Matlab仿真的完整技术文档，适合激光原理、光学工程等课程的学生与科研人员参考。文档从高斯光束数学模型出发，推导强度分布公式，并给出谐振腔中位置处的归一化强度分布仿真方法，覆盖二维、三维强度分布图绘制。同时结合CCD采集的真实激光光斑图像，通过imread读取强度数据，与理论高斯曲线进行对比，展示数学仿真与实际测量的差异。进一步利用角谱衍射法模拟高斯光束在自由空间传播过程中光强、光斑有效截面半径及等相位面的变化，包含核心Matlab代码与参数输入示例，方便读者复现。资源为单份docx文档，体积432KB，内容结构清晰，既有公式推导也有程序实现，可直接用于实验报告撰写或课程设计参考。已有461人学习下载，适合需要快速完成高斯光束仿真作业或理解光场传播原理的读者。

1. 高斯光束的Matlab仿真：先画强度分布，再谈传播

用Matlab画高斯光束，十行代码就能出一张三维图；但把CCD拍到的光斑和理论曲线叠在一起时，峰值位置、半宽、背景噪声都会对不上，三维图还会有一边整体竖起来。问题往往不在公式，而在矩阵读取方向、灰度通道和白边处理。下面这套流程来自激光原理课程里的高斯光束仿真作业，解决两个具体问题：一是用imread读取实际光斑图像，生成二维、三维归一化强度分布并与理论高斯曲线对比；二是用角谱传播模型仿真光束传播不同距离后光斑半径和峰值强度的变化。适合正在做激光实验、光电课程设计，或者想搞懂Matlab图像处理与傅里叶光学仿真边界的读者。

2. 高斯光束数学模型与二维强度矩阵的生成

2.1 基模高斯光束的强度表达式与参数约定

谐振腔中常见的基模场，直角坐标下振幅分布写作

E(x,y)=E0 exp(-(x^2+y^2)/w0^2)

对应的强度分布是 |E|^2，写成

I(x,y)=I0 exp(-2(x^2+y^2)/w0^2)

这里的 w0 是束腰半径，在振幅定义里是振幅下降到中心值 1/e 的位置；在强度定义里，1/e^2 处对应同样半径。很多Matlab示例直接用 exp(-r^2/w0^2) 画“强度分布”，严格说那画的是振幅场。如果只做归一化展示，两者形状相近，但和实测光斑做半宽对比时，强度曲线的半高位置会在 sqrt(2) 倍处出现偏差。做这一题之前，先把要画的是振幅还是强度定下来，后面所有公式、阈值和束宽计算才不会乱。

这份课程作业文档资料里的 M 文件，把理论曲线写成a2=exp(-x2.^2/10)，形式上也是 exp 函数，但它没有按严格的 1/e^2 强度定义写，更像是一个拟合用高斯函数。作业里可以沿用，但在程序注释里最好写明w_fit=sqrt(10)而不是w0，避免后续传播仿真和二维图对比时概念混淆。

2.2 用meshgrid生成二维高斯强度场

二维高斯分布最简单的写法是用 meshgrid 构建网格，再把径向距离平方代入 exp。下面的代码生成 224×224 的归一化强度场，和后面 CCD 图像尺寸保持一致。

N = 224; % 与CCD图像行数保持一致，方便对比 L = 10e-3; % 仿真区域边长，单位m wo = 1e-3; % 束腰半径，单位m dx = L / N; % 空间采样间隔 x = (-N/2 : N/2-1) * dx; % 以0为中心的横坐标 [X, Y] = meshgrid(x, x); % 若按强度分布，指数项为 -2*(X.^2+Y.^2)/wo^2 Gau = exp(-(X.^2 + Y.^2) / wo^2); % 振幅场 Gau_norm = Gau / max(Gau(:)); % 归一化 figure; surf(X*1e3, Y*1e3, Gau_norm, 'EdgeColor', 'none'); xlabel('x (mm)'); ylabel('y (mm)'); zlabel('归一化强度');

这段代码里，N=224不是必须的，后面做 FFT 传播时也会用 500。关键是(-N/2 : N/2-1)这种写法让光斑中心落在网格中央，而不是数组的(1,1)位置；如果直接用1:N，后续fft2和surf都要处理中心偏移，麻烦得多。dx=L/N决定空间分辨率，L取 10mm 时，单个像素相当于约 44.6µm，这个量级和实验中 CCD 的像素尺寸接近，适合做毫米级光斑仿真。绘图时乘 1e3 是为了把米换成毫米，否则横轴刻度会显示成0.001这类小数值。

2.3 参数表：波长、束腰、瑞利长度和采样点

传播仿真中需要提前确定的参数如下。

瑞利长度由zR = k0*wo^2/2计算，代入k0=2*pi/lambda后等于pi*wo^2/lambda。用上面的参数算出来约 5.53m，所以文档中z=100000mm即 100m 时，已经远大于瑞利长度，光斑半径接近线性扩展，峰值强度明显下降。运行结果里峰值从 1 降到 0.21，是这个参数下的正常物理趋势，不是程序错误。

2.4 为什么要做归一化

CCD 读出的灰度值不是光功率绝对值，单位是“数”，和理论公式中的 I0 没有可比性。所以画二维强度分布图前，一定要把数据除以自身最大值。更稳妥的做法是先减背景再归一化：I_norm=(I-min(I))/(max(I)-min(I))。原始光斑照片如果有暗电流或环境光，底部不是 0，直接除以最大值会把背景也纳入比较，导致理论曲线和实验曲线的边缘不贴合。这一行处理在后面的 imread 实验数据中同样适用。

3. 用imread读取CCD光斑：二维/三维强度分布与白边处理

3.1 读取图片前先看矩阵维度

imread是 Matlab 读取位图最常用的命令，但读到的是什么尺寸，取决于图片本身是灰度图还是 RGB 图。

A = imread('D:\documents\作业\激光原理与应用\高斯.bmp'); disp(size(A)); imshow(A); axis off; title('CCD采集的高斯光束光强分布');

如果size(A)返回224 244，说明是单通道灰度图；如果返回224 244 3，说明是 RGB 三通道，后面所有取行、取列的操作都必须指定第三维下标，否则A(:,122)拿到的是 244×3 的矩阵，plot会画出三条线。原文档写的是 224×244 矩阵，但代码里又出现了[high, width, color] = size(A)，说明实际文件很可能是 RGB 格式。常见做法是先用size确认，然后取第一通道或者用rgb2gray转成灰度。imshow对uint8数据会自动按灰度范围映射，直接看即可，不要先用double转换再 imshow，否则显示结果可能变黑。

3.2 提取直径方向数据绘制二维强度分布

原始程序的读取核心是两行：

A1 = A(:, 122); % 取第122列，长度224 x1 = 1:1:224; figure; plot(x1, A1, 'LineWidth', 1, 'Color', 'r'); title('理论高斯曲线');

这里要特别注意“中间一行”和“中间一列”的差别。题目文字说的是“取中间一行（122行）”，但代码A(:,122)取的是第 122 列。如果矩阵是 224 行 × 244 列，那么第 122 列对应的剖面长度是 224，和第 122 行对应的剖面长度 244 并不相同。两个方向理论上都可以反映高斯剖面，但横轴刻度必须和取的方向一致。更好的写法是：

if size(A,3) == 3 A = A(:,:,1); % 取单通道，后续都是二维矩阵 end line_data = double(A(122, :)); % 第122行，长度244 line_data = line_data / max(line_data); % 归一化到[0,1] Nx = length(line_data); xpix = (1:Nx) - (Nx+1)/2; % 以图像中心为原点 figure; plot(xpix, line_data, 'r', 'LineWidth', 1);

double转换是必须的，因为imread返回的是uint8，直接做除法时 Matlab 会按整数运算取整。归一化后，纵轴可以统一到 0～1，方便和理论曲线放在同一张图里。

3.3 理论高斯曲线与实验曲线的叠加

文档里的理论曲线是x2=-100:1:100; a2=exp(-x2.^2/10)。它的 x 范围是 [-100,100]，和实验曲线的像素范围 [1,224] 不在同一坐标基准上，所以原图把两条曲线分开画，看起来各自都像高斯，叠在一起却对不上。叠加时要做两件事：把实验曲线中心移到 0，再让理论曲线的横轴范围和实验曲线匹配。

x_theory = -100:0.1:100; I_theory = exp(-x_theory.^2 / 10); figure; plot(xpix, line_data, 'r', 'LineWidth', 1); hold on; plot(x_theory, I_theory, 'b--', 'LineWidth', 1); legend('实验数据', '理论曲线'); xlabel('相对中心像素'); ylabel('归一化强度');

如果发现实验峰值不在 x=0 附近，说明光斑中心没有落在图像中心，可以用[~, idx] = max(line_data)把峰值位置平移到 0 后再比较。像素坐标到毫米坐标需要标定像素间距，原题没有给出该参数，所以这一步只能做“形”的比较。严格的做法是用光斑的实际直径和像素数算出pixel_size_mm，再把xpix乘上它。

3.4 用mesh命令画三维强度分布并去白边

二维剖面只能看到一条直径上的分布，完整光斑要用整幅矩阵画三维图。原文档里的做法是：

[high, width, color] = size(A); x = 1:width; y = 1:high-1; mesh(x', y', double(A(2:224,:,1))); grid on; xlabel('x'); ylabel('y'); zlabel('z'); title('三维强度分布');

这段代码已经处理了顶部白边：A(2:224,:,1)从第 2 行开始取，跳过了第 1 行的亮线。但它写死了 224，如果图片尺寸变化会越界。我一般会写成：

A_gray = double(A(:,:,1)); A_crop = A_gray(2:end-1, :); % 去掉顶部和底部各一行 [Ny, Nx] = size(A_crop); x = 1:Nx; y = 1:Ny; figure; mesh(x, y, A_crop); grid on; xlabel('x (pixel)'); ylabel('y (pixel)'); zlabel('强度'); title('三维强度分布');

这里2:end-1比2:224更通用。白边在 mesh 图里的表现是 z 轴一整边被拉高，因为白边灰度接近 255，而光斑中心灰度可能不到 200，如果不剔除，三维图会有一条边整体竖起来，光斑本身的形状反而看不清楚。如果白边在侧边而不是顶部，可以通过imshow(A(1:5,:))或size(A)确认，再相应调整列方向的范围。

注意：mesh 的输入必须是 double 类型，uint8 矩阵在三维坐标下可能不按数值直接渲染，会出现颜色被当作索引的问题。

4. 角谱法传播仿真：不同z位置的光斑半径与强度变化

4.1 为什么是角谱法而不是直接积分

高斯光束在自由空间中传播，理论上可以用菲涅尔衍射积分，但离散化之后每个观察点都要做一次二维求和，计算量很大。角谱法的思路是把初始光场分解成平面波，乘上自由空间传递函数，再逆变换回空间域。由于只涉及两次 FFT 和一次矩阵乘法，扫描多个 z 值时速度很快，非常适合画“传播 10m、20m、50m”这组图。

近轴条件下，传递函数可以写成：

H(fx,fy)=exp(jz(kx^2+ky^2)/(2*k0))

其中 kx=2pifx，ky=2pify。这个公式是原始 M 文件的核心，也是高斯光束从束腰传播到远场时振幅和束宽变化的数值基础。适用前提是传播角足够小、介质均匀、没有增益或损耗；激光器谐振腔外检测正好满足这些条件。

4.2 角谱传播代码与fftshift配合

下面这段脚本可以直接替换原 M 文件中的传播部分，解决中心偏移和 ifftshift 缺失的问题。

clear; close all; N = 500; % 采样点数，100~500 L = 10e-3; % 仿真区域宽度，单位m dx = L / N; % 空间采样间隔 x = (-N/2 : N/2-1) * dx; y = x; [X, Y] = meshgrid(x, y); lambda = 0.568e-6; % 波长，单位m k0 = 2*pi / lambda; wo = 1e-3; % 束腰半径，单位m zR = k0 * wo^2 / 2; % 瑞利长度 fprintf('Rayleigh range: %.2f m\n', zR); Gau = exp(-(X.^2 + Y.^2) / wo^2); % 初始振幅分布 fx = (-N/2 : N/2-1) / (N*dx); % 空间频率，单位1/m [FX, FY] = meshgrid(fx, fx); z = input('Propagation distance z (m): '); H = exp(1j * z / (2*k0) * ((2*pi*FX).^2 + (2*pi*FY).^2)); FGau = fftshift(fft2(Gau)); % 零频移到矩阵中心 Gau_pro = ifft2(ifftshift(FGau .* H)); % 乘H后逆shift再ifft2 figure; subplot(1,2,1); mesh(x*1e3, y*1e3, abs(Gau)); title('Initial Gaussian Beam'); xlabel('x (mm)'); ylabel('y (mm)'); zlabel('Amplitude'); subplot(1,2,2); mesh(x*1e3, y*1e3, abs(Gau_pro)); title(['z=', num2str(z), ' m']); xlabel('x (mm)'); ylabel('y (mm)'); zlabel('Amplitude');

代码的要点是频域坐标必须和fft2的输出顺序一致。fft2的零频在矩阵的(1,1)位置，fftshift后零频在中心；同时H是按零频在中心构造的，所以两者可以直接相乘。乘完之后要再ifftshift，把零频挪回(1,1)再送进ifft2。原始 M 文件最容易忽略这一步，结果是光斑在空域发生循环平移，峰值下降，数值束宽也可能从理论上的十几毫米变成一两毫米。

fx = (-N/2 : N/2-1) / (N*dx)是空间频率轴，单位为 1/m。频域间隔等于 1/(N*dx)=1/L，也就是仿真区域越宽，频域越密，角谱法的分辨率越高。H中z/(2*k0)这个系数决定了相位累积的快慢；当 z 很大时，H 在频域边缘振荡很快，N 不够就会发生仿真发散。

4.3 扫描z值并比较数值束宽与理论束宽

要观察光斑随传播距离的变化，不能只跑一次，可以把 z 写成一个数组循环。理论束宽用wz = wo*sqrt(1+(z/zR)^2)。

z_list = [0.1, 1, 5, 10, 20, 50, 100]; w_num = zeros(size(z_list)); w_theory = zeros(size(z_list)); FGau = fftshift(fft2(Gau)); for k = 1:length(z_list) zk = z_list(k); Hk = exp(1j * zk / (2*k0) * ((2*pi*FX).^2 + (2*pi*FY).^2)); Gk = ifft2(ifftshift(FGau .* Hk)); Ik = abs(Gk).^2; Ik = Ik / max(Ik(:)); line_prof = Ik(N/2+1, :); pos1 = find(line_prof > exp(-2), 1, 'first'); pos2 = find(line_prof > exp(-2), 1, 'last'); w_num(k) = (x(pos2) - x(pos)) / 2; w_theory(k) = wo * sqrt(1 + (zk / zR)^2); end figure; semilogy(z_list, w_num*1e3, 'ro-', z_list, w_theory*1e3, 'b^-'); legend('数值束宽', '理论束宽', 'Location', 'northwest'); xlabel('z (m)'); ylabel('束宽 (mm)');

line_prof取的是通过中心的水平强度剖面，阈值exp(-2)对应强度降到峰值 1/e^2 的位置。find找到左右两个越过阈值的点，相减再除以 2 就是光斑半径。如果pos1为空，说明光斑已经扩展出仿真区域，需要增大L或N。原文档运行结果中，100m 处的数值束宽是 1.94mm，理论值 18.107mm，二者差一个数量级，最可能的原因就是fftshift/ifftshift不配对；修正后 100m 处会明显接近理论束宽。

不同传播距离下的束宽变化趋势可以总结成下表，供报告排版参考。

表中的理论值按wz=wo*sqrt(1+(z/zR)^2)计算，z_R=5.53m。数值仿真的峰值强度也会随 z 增大而下降，但要注意abs(Gau_pro)是振幅，画强度图时要取平方再归一化。

5. 调试细节：白边剔除、采样点与束宽计算的一致性

5.1 白边不一定要手动数行数

三维图画出来一边高，先不要急着改数据。用下面几行定位白边位置：

imshow(A(1:5, :, 1)); % 查看顶部5行 figure; imagesc(double(A(:,:,1))); colorbar; % 完整显示强度范围

白边在imagesc下会是一条接近 colorbar 顶部的亮色条，比直接看 mesh 图更直观。确认白边只有顶部一行时，用A(2:end-1, :)去掉首尾；如果白边有十行，就把前面范围改成A(11:end-10, :)。不要写死A(2:224,:,1)，换一张图就容易越界。

5.2 仿真发散时的三个参数调整

高频相位因子在频域边缘振荡时，会出现边缘条纹或能量不衰减，这就是常见的“仿真发散”。按以下顺序调整：

1. 提高采样点数：N从 100 加到 500，频域采样更密；
2. 缩小仿真区域：L太大时光斑只占很少像素，相位变化集中在少数网格上，建议L为光斑直径的 3～5 倍；
3. 分段传播：z 超过几十米时，把 z 拆成若干段循环传播，每段不超过瑞利长度，可以显著减少混叠。

每次调整后检查中心行剖面，如果边缘仍有周期性起伏，说明还没收敛。

5.3 束宽定义统一，才能对比数值和理论

理论束宽公式里的 w(z) 对应振幅降到 1/e 的位置，换算成强度就是 1/e^2。所以代码里数值束宽也应该按强度阈值exp(-2)找边缘，才能和wz=wo*sqrt(1+(z/zR)^2)对齐。如果按半高 FWHM 找，结果会差约 1.18 倍，不要混用。

把像素坐标换成空间坐标时，先标定像素间距。CCD 图像里的一个像素对应多少毫米，应通过光斑实际直径和像素数算出来，而不是对矩阵下标直接取abs。否则二维图、三维图和传播仿真四张图之间的束宽永远对不上。

本文还有配套的精品资源，点击获取