简介：这份资源是面向激光物理、光学工程方向的学习者与科研人员整理的MATLAB仿真代码包，聚焦高斯光束传输与谐振腔设计这一激光系统建模中的核心问题。包内共95个文件，以67个m脚本为主体，配合17个fig图形结果、5个mat数据文件，以及docx说明、pdf文档、md与txt说明等，压缩包约27.37MB，覆盖从基础计算到结果可视化的完整流程。内容围绕高斯光束传播模型、热效应分析、非线性频率变换、光束调节与匹配算法、稳定驻波谐振腔与行波谐振腔仿真、ABCD矩阵计算等模块展开，可帮助读者理解束宽变化、焦点位置、热透镜效应、倍频混频效率及模式匹配等关键环节。已有122人学习，适合需要借助MATLAB快速搭建光路模型、验证谐振腔稳定性或优化激光输出参数的读者参考使用。

1. 从一束光说起：MATLAB 里把激光高斯光束传输与谐振腔仿真跑通

很多做激光系统的人都有个习惯：先在纸上画光路，再在软件里把光路“跑”一遍。问题是，纸上算高斯光束的束腰、发散角、ABCD 矩阵，手推几页公式还行，一旦涉及热效应、非线性频率变换、驻波腔和行波腔的联合仿真，手算基本就翻车了。这套基于 MATLAB 的激光高斯光束传输与谐振腔仿真代码包，解决的正是这个痛点：用脚本把高斯光束传播模型、ABCD 矩阵计算、谐振腔稳定性判断、热效应分析和非线性频率变换串成一条可复现的链路。它适合做固体激光器、光纤激光器、倍频与和频系统、腔体设计的工程师和研究生，也适合想用 MATLAB 把光学公式变成可调参数工具的人。下面我按自己搭这套东西的顺序，把关键模型、参数和踩过的坑讲清楚。

2. 高斯光束传播模型：从束腰到 ABCD 矩阵的完整链路

2.1 为什么先用 q 参数而不是直接追光线

几何光学追光线在自由空间里是直线，到了透镜就折一下，看起来简单。但激光束不是一条线，它是有束腰、有曲率半径、有发散角的复参数光束。高斯光束的复参数 q 把束腰位置、瑞利长度、波前曲率全打包在一起，传播和经过光学元件时只需要做一次代数变换。常见做法是：

• 自由空间传播距离 L：q₂ = q₁ + L
• 薄透镜焦距 f：1/q₂ = 1/q₁ − 1/f
• 球面镜反射：等效于焦距 f = R/2 的透镜

这套规则的好处是，任何由自由空间、透镜、镜片组成的轴对称光路，都能写成一个 2×2 的 ABCD 矩阵，然后对 q 做分式线性变换。MATLAB 里实现起来就是矩阵乘法和复数运算，不需要解波动方程。

2.2 用 MATLAB 写一个最小可跑的高斯光束传播脚本

下面这段代码是我一般用来验证光路是否搭对的最小示例。它定义了一个束腰在 z=0 处的高斯光束，然后让它经过一段自由空间、一个薄透镜，再传播一段距离，最后输出束腰半径和波前曲率。

% 高斯光束 q 参数传播最小示例 lambda = 1064e-9; % 波长，单位 m w0 = 0.5e-3; % 束腰半径，单位 m zR = pi * w0^2 / lambda; % 瑞利长度 q0 = 1i * zR; % 束腰处 q 参数 % 光学系统：自由空间 L1 -> 薄透镜 f -> 自由空间 L2 L1 = 0.2; % 第一段自由空间，单位 m f = 0.1; % 透镜焦距，单位 m L2 = 0.3; % 第二段自由空间，单位 m % 自由空间 ABCD 矩阵 M_free = @(L) [1, L; 0, 1]; % 薄透镜 ABCD 矩阵 M_lens = @(f) [1, 0; -1/f, 1]; M = M_free(L2) * M_lens(f) * M_free(L1); q_out = (M(1,1)*q0 + M(1,2)) / (M(2,1)*q0 + M(2,2)); % 从 q 参数反推束腰半径和曲率半径 inv_q = 1 / q_out; R = 1 / real(inv_q); % 波前曲率半径 w = sqrt(lambda / (pi * abs(imag(inv_q)))); % 当前光束半径 fprintf('输出 q = %.4f + %.4fi\n', real(q_out), imag(q_out)); fprintf('光束半径 w = %.4f mm\n', w*1e3); fprintf('曲率半径 R = %.4f m\n', R);

逻辑说明：q0 是束腰处的复参数，实部为 0，虚部为瑞利长度。M_free 和 M_lens 分别生成自由空间和薄透镜的 ABCD 矩阵。矩阵乘法顺序是从右到左，对应光先经过 L1，再经过透镜，最后经过 L2。q_out 用分式线性变换得到。最后从 1/q 的实部和虚部分别提取曲率半径和光束半径。

参数说明：lambda 换成 532 nm 就是倍频光，w0 决定瑞利长度和发散角，L1、f、L2 根据实际光路填写。如果 q_out 的虚部变成负数，说明光束已经过了焦点，这是正常现象，不是错误。

2.3 ABCD 矩阵计算里最容易搞混的符号和顺序

ABCD 矩阵本身不难，难的是符号约定和乘法顺序。我见过不少人把自由空间矩阵写成 [1, 0; L, 1]，结果光束越传越收敛，明显不对。自由空间是 [1, L; 0, 1]，透镜是 [1, 0; -1/f, 1]，球面镜反射等效焦距 f = R/2，其中 R 是曲率半径，凹面镜 R 取正，凸面镜 R 取负。乘法顺序是后经过的元件放在左边，先经过的放在右边。如果光路里有反射，要注意传播方向反转后 q 参数的符号处理，常见做法是把反射后的光路重新定义坐标方向，而不是直接套用透射矩阵。

提示：每次改完光路，先用一个已知的简单例子验证，比如束腰在透镜前焦面，经过透镜后束腰应该在透镜后焦面，且大小满足 w0' = lambdaf/(piw0)。如果这个都对不上，后面的谐振腔和热效应分析就不用看了。

3. 稳定驻波谐振腔设计：从 g 参数到自洽场

3.1 驻波腔的稳定性判据和 g 参数

驻波腔由两面镜子组成，光在两面镜子之间来回反射。稳定腔的条件是 0 < g1*g2 < 1，其中 g1 = 1 - L/R1，g2 = 1 - L/R2，L 是腔长，R1、R2 是两面镜子的曲率半径。这个判据的本质是，光在腔内往返一次后，光束参数必须回到自身，否则光会跑出腔外。MATLAB 里判断稳定性只需要几行代码，但实际设计时还要考虑镜面热透镜效应，因为高功率激光的镜面会吸收部分光，产生热变形，等效于改变镜子的曲率半径。

3.2 用往返矩阵求稳定腔的基模参数

稳定腔的基模光束参数可以通过往返矩阵的自洽条件求出。往返矩阵 M_rt 是光在腔内走一圈的 ABCD 矩阵，自洽条件是 q = (Aq + B)/(Cq + D)。解这个方程可以得到 q 的解析表达式，然后反推束腰位置和大小。

% 驻波腔基模参数计算 L = 0.5; % 腔长，单位 m R1 = 1.0; % 镜 1 曲率半径，单位 m R2 = 1.0; % 镜 2 曲率半径，单位 m lambda = 1064e-9; % 波长 % 单程矩阵：从镜 1 到镜 2 M_half = [1, L; 0, 1] * [1, 0; -2/R1, 1]; % 往返矩阵：镜 1 -> 镜 2 -> 镜 1 M_rt = [1, 0; -2/R2, 1] * M_half; A = M_rt(1,1); B = M_rt(1,2); C = M_rt(2,1); D = M_rt(2,2); % 稳定性判据 g1 = 1 - L/R1; g2 = 1 - L/R2; if g1*g2 > 0 && g1*g2 < 1 fprintf('腔是稳定的，g1*g2 = %.4f\n', g1*g2); else fprintf('腔不稳定，g1*g2 = %.4f\n', g1*g2); end % 自洽解 q = (A-D)/(2C) ± sqrt(((A-D)/(2C))^2 - B/C) disc = ((A-D)/(2*C))^2 - B/C; q_sol = (A-D)/(2*C) + sqrt(disc); inv_q = 1/q_sol; w = sqrt(lambda / (pi * abs(imag(inv_q)))); fprintf('基模光束半径 w = %.4f mm\n', w*1e3);

逻辑说明：M_half 是从镜 1 到镜 2 的单程矩阵，包含自由空间传播和镜 1 的反射。M_rt 是往返矩阵，在单程基础上再加镜 2 的反射和反向传播。自洽解用二次方程求根，取物理上合理的那个根。最后从 q 的虚部提取光束半径。

参数说明：L 是腔长，R1、R2 是镜面曲率半径，lambda 是波长。如果 g1*g2 接近 0 或 1，腔处于临界稳定，实际中容易受扰动失稳，一般会留一定余量。镜面热透镜效应可以通过把 R1、R2 替换为等效热透镜焦距来近似。

3.3 热效应分析：镜面变形和热透镜怎么进模型

高功率激光腔里，镜面吸收的光会让镜面局部升温，产生热膨胀和折射率变化，等效于一个热透镜。常见做法是把热透镜当作一个薄透镜放在镜面附近，焦距由吸收功率、镜面材料和冷却条件决定。MATLAB 里可以先用一个经验公式估算热透镜焦距，然后把它代入 ABCD 矩阵重新计算腔的稳定性和基模参数。

% 热透镜焦距估算与腔参数修正 P_abs = 10; % 镜面吸收功率，单位 W kappa = 1.2; % 材料热导率，单位 W/(m*K) dn_dT = 2.5e-6; % 折射率温度系数，单位 1/K alpha = 1e-6; % 热膨胀系数，单位 1/K w_pump = 1e-3; % 泵浦光半径，单位 m % 经验公式：热透镜焦距 f_th = pi*kappa*w_pump^2 / (P_abs*(dn_dT + alpha)) f_th = pi * kappa * w_pump^2 / (P_abs * (dn_dT + alpha)); fprintf('热透镜焦距 f_th = %.4f m\n', f_th); % 把热透镜插入腔矩阵 M_th = [1, 0; -1/f_th, 1]; M_rt_th = M_th * M_rt; % 重新计算稳定性和基模参数，步骤同上

逻辑说明：热透镜焦距的经验公式把吸收功率、材料热导率、折射率温度系数和热膨胀系数放在一起，给出一个量级估算。实际中这个公式需要根据具体材料和冷却条件修正，但作为初步设计足够。把热透镜矩阵乘到往返矩阵里，就能看到热效应对腔稳定性的影响。

参数说明：P_abs 是镜面实际吸收的功率，不是入射功率。kappa、dn_dT、alpha 是材料参数，不同玻璃和晶体差别很大。w_pump 是泵浦光在镜面上的半径，越小热透镜越强。如果算出来的 f_th 和腔长可比，说明热效应不能忽略。

注意：热透镜焦距对吸收功率非常敏感，P_abs 翻倍，f_th 减半。实际中镜面吸收功率很难精确知道，一般会留 20% 到 50% 的余量，或者用可调热透镜模型做灵敏度分析。

4. 行波谐振腔仿真：环形腔里的方向性和模式匹配

4.1 行波腔和驻波腔的本质区别

驻波腔里光来回反射，形成驻波，腔内有波节和波腹。行波腔通常是环形结构，光沿一个方向循环传播，没有驻波效应。行波腔的好处是避免了空间烧孔效应，适合需要单频运转的激光器。但行波腔的稳定性分析比驻波腔复杂，因为光在腔内走一圈要经过多个镜片和可能的光学元件，ABCD 矩阵的维度更高。MATLAB 里处理行波腔的思路和驻波腔类似，只是往返矩阵变成绕腔一圈的单程矩阵。

4.2 环形腔的 ABCD 矩阵和模式匹配

环形腔一般由多个镜片组成，光绕一圈回到起点。模式匹配的目标是让光在腔内循环时，光束参数在每一点都保持合理，尤其是增益介质和非线性晶体所在的位置，光束半径要匹配泵浦光或满足相位匹配条件。

% 环形腔模式匹配示例：四镜环形腔 L1 = 0.3; L2 = 0.2; L3 = 0.3; L4 = 0.2; % 四段自由空间 R1 = 0.5; R2 = 0.5; R3 = 0.5; R4 = 0.5; % 四面镜曲率半径 lambda = 1064e-9; % 绕腔一圈的 ABCD 矩阵 M = [1, 0; -2/R4, 1] * [1, L4; 0, 1] * ... [1, 0; -2/R3, 1] * [1, L3; 0, 1] * ... [1, 0; -2/R2, 1] * [1, L2; 0, 1] * ... [1, 0; -2/R1, 1] * [1, L1; 0, 1]; % 自洽解 A = M(1,1); B = M(1,2); C = M(2,1); D = M(2,2); disc = ((A-D)/(2*C))^2 - B/C; q_sol = (A-D)/(2*C) + sqrt(disc); inv_q = 1/q_sol; w = sqrt(lambda / (pi * abs(imag(inv_q)))); fprintf('环形腔基模光束半径 w = %.4f mm\n', w*1e3); % 检查增益介质位置的光束半径 % 假设增益介质在 L1 中间，传播 L1/2 后计算 M_gain = [1, L1/2; 0, 1]; q_gain = (M_gain(1,1)*q_sol + M_gain(1,2)) / (M_gain(2,1)*q_sol + M_gain(2,2)); w_gain = sqrt(lambda / (pi * abs(imag(1/q_gain)))); fprintf('增益介质处光束半径 w_gain = %.4f mm\n', w_gain*1e3);

逻辑说明：环形腔的 ABCD 矩阵按光路顺序从右到左乘，每个镜片用薄透镜矩阵等效，自由空间用传播矩阵。自洽解和驻波腔一样，用二次方程求根。增益介质处的光束半径通过从参考面传播到增益介质位置再计算。

参数说明：L1 到 L4 是四段自由空间长度，R1 到 R4 是四面镜的曲率半径。如果增益介质处光束半径太小，功率密度过高可能损伤晶体；太大则泵浦光重叠效率低。一般会调整镜片曲率半径和距离，让增益介质处光束半径在 0.2 到 0.5 mm 之间。

4.3 行波腔里的非线性频率变换位置选择

非线性频率变换，比如倍频和和频，对光束半径和相位匹配非常敏感。倍频效率与基频光功率密度的平方成正比，所以通常会把非线性晶体放在光束半径最小的位置，也就是束腰附近。但束腰太小会导致晶体损伤，太大则效率低。MATLAB 里可以扫描晶体位置，计算不同位置的光束半径和对应的倍频效率，找到最优位置。

% 扫描非线性晶体位置，找最优倍频效率 positions = linspace(0, L1, 50); % 在 L1 段内扫描 eta = zeros(size(positions)); for k = 1:length(positions) M_pos = [1, positions(k); 0, 1]; q_pos = (M_pos(1,1)*q_sol + M_pos(1,2)) / (M_pos(2,1)*q_sol + M_pos(2,2)); w_pos = sqrt(lambda / (pi * abs(imag(1/q_pos)))); % 简化倍频效率模型：与 1/w_pos^2 成正比，但受损伤阈值限制 eta(k) = 1 / (w_pos^2 + (0.3e-3)^2); % 0.3 mm 为损伤阈值对应的最小半径 end [~, idx] = max(eta); fprintf('最优晶体位置 = %.4f m，光束半径 = %.4f mm\n', ... positions(idx), w_pos*1e3);

逻辑说明：扫描晶体在 L1 段内的位置，每次计算该位置的光束半径，然后用一个简化模型估算倍频效率。模型里加了一个最小半径限制，模拟损伤阈值。实际中倍频效率还和相位匹配角、晶体长度、走离效应有关，但位置扫描是第一步。

参数说明：positions 是扫描范围，一般覆盖整个自由空间段。损伤阈值对应的最小半径根据晶体材料不同，常见 KTP 和 LBO 的损伤阈值在 GW/cm² 量级，对应光束半径在 0.1 到 0.5 mm 之间。如果最优位置在镜片附近，要考虑镜片像差和热效应的影响。

5. 避坑与排查：高斯光束和谐振腔仿真里最容易翻车的 5 个地方

5.1 现象：q 参数传播后光束半径变成 NaN 或无穷大

原因：ABCD 矩阵里出现了除零，或者 q 参数在某个位置正好等于 −D/C，导致分式线性变换的分母为零。常见于透镜焦距设为零，或者腔长正好等于镜片曲率半径的两倍。

解决：检查所有光学元件的参数，确保焦距和曲率半径不为零。在代码里加一个判断，如果分母绝对值小于 1e-12，就报错并输出当前矩阵和 q 值。另外，腔长等于 2R 时腔处于临界稳定，实际中应该避开这个点。

5.2 现象：驻波腔稳定性判据算出来是稳定的，但实际光路不出光

原因：稳定性判据只考虑了几何光学，没有考虑热透镜、镜面吸收和增益介质的折射率变化。高功率泵浦下，增益介质本身就是一个强热透镜，等效于改变腔的 g 参数。

解决：把增益介质的热透镜矩阵加入往返矩阵，重新计算 g1*g2。如果热透镜焦距和腔长可比，腔可能从稳定变成不稳定。常见做法是先用低功率泵浦测出冷腔参数，再逐步增加功率，观察输出功率和光束质量的变化。

5.3 现象：行波腔模式匹配算出来的光束半径和实验测的差很多

原因：环形腔的镜片不是理想薄透镜，有像差和反射相位延迟。另外，光在腔内循环时，镜片的实际曲率半径可能因为安装应力或热变形而改变。

解决：先用一个已知曲率半径的标准镜片验证 ABCD 矩阵，确认符号和乘法顺序正确。然后测量实际镜片的曲率半径，或者用干涉仪测波前。如果差得不多，可以在模型里加一个等效焦距修正项，通过实验数据反推。

5.4 现象：非线性频率变换效率扫描结果和实验趋势相反

原因：简化模型只考虑了光束半径，忽略了相位匹配和走离效应。如果晶体位置改变了相位匹配角，效率可能反而下降。另外，光束半径太小会导致晶体损伤，实验上不敢用最优位置。

解决：在扫描模型里加入相位匹配角和晶体长度的修正，或者直接用实验测得的效率曲线反推最优位置。如果损伤阈值限制严格，就把最小半径设大一点，牺牲一点理论效率换安全。

5.5 现象：MATLAB 脚本跑得慢，扫描几百个点要等很久

原因：每次循环都重新计算矩阵乘法和自洽解，没有向量化。MATLAB 的循环效率不高，尤其是涉及复数运算时。

解决：把扫描参数写成向量，用矩阵运算一次性计算所有位置。比如把 positions 写成行向量，用 arrayfun 或者直接矩阵乘法批量处理。另外，预分配数组，避免在循环里动态增长。如果还是慢，可以把自洽解写成函数，用 parfor 并行计算。

提示：仿真结果和实验对不上时，先检查单位。波长用米，束腰用米，焦距用米，不要混用毫米和米。我见过有人把 w0 写成 0.5 毫米，但 lambda 用 1064 纳米，结果瑞利长度算出来差三个数量级。

6. 进阶技巧：用 MATLAB 把仿真和实验数据串起来

6.1 用实测光束半径反推 q 参数和热透镜焦距

仿真做得再好，最终还是要和实验对上。我一般会先用 CCD 相机测几个位置的光束半径，然后用最小二乘拟合高斯光束的传播公式，反推束腰位置、束腰半径和热透镜焦距。MATLAB 的 lsqcurvefit 或者 fminsearch 都能做这件事。

% 用实测光束半径拟合高斯光束参数 z_data = [0.1, 0.2, 0.3, 0.4, 0.5]; % 测量位置，单位 m w_data = [0.4, 0.35, 0.32, 0.35, 0.42] * 1e-3; % 实测半径，单位 m lambda = 1064e-9; % 拟合模型：w(z) = w0 * sqrt(1 + ((z - z0)/zR)^2) model = @(p, z) p(1) * sqrt(1 + ((z - p(2))/(pi*p(1)^2/lambda)).^2); p0 = [0.3e-3, 0.3]; % 初始猜测：w0 和 z0 p_fit = lsqcurvefit(model, p0, z_data, w_data); fprintf('拟合束腰半径 w0 = %.4f mm\n', p_fit(1)*1e3); fprintf('拟合束腰位置 z0 = %.4f m\n', p_fit(2));

逻辑说明：model 是高斯光束传播公式，p(1) 是束腰半径，p(2) 是束腰位置。lsqcurvefit 用实测数据拟合这两个参数。拟合完可以反推瑞利长度和发散角。

参数说明：z_data 和 w_data 是实测数据，测量点至少 5 个，覆盖束腰前后。如果拟合残差大，可能是热透镜效应导致光束不是理想高斯，需要在模型里加修正项。

6.2 把热透镜焦距作为拟合参数

如果实测数据偏离理想高斯，可以把热透镜焦距也作为拟合参数，用 ABCD 矩阵模型拟合。这样能同时得到冷腔参数和热透镜焦距。

% 用 ABCD 矩阵模型拟合热透镜焦距 function w_model = abcd_model(p, z_data, lambda, L, R1, R2) f_th = p(1); % 热透镜焦距 w0 = p(2); % 冷腔束腰 z0 = p(3); % 束腰位置 q0 = 1i * pi * w0^2 / lambda; q0 = q0 + z0; % 把束腰移到 z0 w_model = zeros(size(z_data)); for k = 1:length(z_data) M = [1, z_data(k); 0, 1] * [1, 0; -1/f_th, 1]; q = (M(1,1)*q0 + M(1,2)) / (M(2,1)*q0 + M(2,2)); w_model(k) = sqrt(lambda / (pi * abs(imag(1/q)))); end end

逻辑说明：这个模型把热透镜放在腔的起点，然后传播到测量位置。拟合参数包括热透镜焦距、冷腔束腰和束腰位置。实际使用时，可以把腔的 ABCD 矩阵完整写进去，而不是只用一个薄透镜近似。

参数说明：f_th 初始猜测可以根据经验公式给，w0 和 z0 根据冷腔测量给。如果拟合不收敛，先固定 f_th，拟合 w0 和 z0，再放开 f_th。

6.3 一个我常用的习惯：每次改光路先跑一遍 ABCD 矩阵检查

我自己的习惯是，不管光路多复杂，先在 MATLAB 里把 ABCD 矩阵写出来，跑一遍 q 参数传播，看看光束半径在关键位置是否合理。如果某个位置光束半径突然变得很大或很小，说明矩阵符号或顺序有问题。这个检查只需要几秒钟，但能省掉后面几个小时的调试。仿真和实验对不上时，先怀疑单位，再怀疑符号，最后怀疑热效应。希望帮到你。

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