简介：一份面向激光技术学习者与科研人员的主动调Q固体激光器Matlab仿真文件包，聚焦四能级系统建模与声光调Q原理，通过速率方程描述粒子在基态、亚稳态、上能级与激发态间的转移，帮助理解粒子数反转、受激辐射及短脉冲形成过程，适合理论学习与科研入门。压缩包内共2个文件，均为.m脚本，其中一个为可执行的主程序，另一个为速率方程求解函数，共同完成从粒子浓度到输出波形的建模。可运行并修改泵浦功率、增益系数、AOM开启时间等参数，观察脉冲宽度、峰值功率、能量效率等指标变化，还可对照被动调Q例程，从调制机制上比较两种方案的差异。整包仅919B，轻量便携，便于快速部署与二次开发。目前已有2093人学习下载，其价值在于不仅给出可运行仿真代码，还通过参数敏感性分析引导读者掌握激光器优化思路，对比不同泵浦速率下的输出，直观理解调Q开关对脉冲成形的控制机制，为课题预研、实验验证或技术方案比对提供直接工具，值得反复调试研究。 接到“主动调Q固体激光器MATLAB仿真”这类任务时，我第一反应不是写代码，而是先把扔了很久的《激光原理》翻出来。因为调Q仿真最大的坑不在编程，而在物理模型和参数单位。你搜到的那些示例程序，十个里有八个符号体系和你的教材对不上，直接跑往往要么“负光子数”蹦出来，要么脉冲宽度仿真出飞秒量级的不合理结果。这篇文章我按自己实际跑通顺的路径，把速率方程建立、MATLAB数值求解、结果验证一整套思路拆开讲，最后附上我踩过的几个典型问题。激光原理课设、毕业设计，或者想通过仿真理解调Q物理过程的工程师，都可以直接参考这套框架。
1. 主动调Q仿真第一步不是写代码，而是把速率方程的单位先理顺
1.1 调Q物理过程与哪些参数决定了脉冲形态
主动调Q的原理，说直白点就是“先憋着，再放大”。泵浦阶段，Q开关让谐振腔处于高损耗状态，腔内无法起振，反转粒子数在增益介质里大量积累；等积累到一定程度，外部信号控制Q开关瞬间把损耗降下去，此时增益远大于损耗，腔内光子数像雪崩一样暴涨，受激辐射在极短时间内把储存的能量释放出来，形成巨脉冲。
仿真要做的就是把这个过程变成数学表达。最常用的办法是求解两个耦合的一阶常微分方程：一个描述腔内光子数密度随时间的变化，另一个描述反转粒子数密度随时间的变化。这两个量决定了脉冲的核心特征，包括脉冲宽度、峰值功率、建立时间和脉冲能量。对于典型的灯泵浦或者LD泵浦Nd:YAG固体激光器，增益介质和腔型的参数大致在一个很窄的范围内，仿出来的结果应该落在合理区间，这也是后面验证结果时最重要的依据。
要理解决定脉冲形态的关键参数，首当其冲是初始反转粒子数密度n0，它由泵浦功率和泵浦时间决定，直接决定了储能大小；其次是腔内光子寿命τc，它反映了谐振腔损耗的大小，输出镜透过率越高，损耗越大，τc越短；增益介质本身的受激发射截面σ，比如Nd:YAG的σ大约是2.8e-19 cm²，Nd:YVO₄的σ更大，约2.5e-18 cm²量级，截面大的介质更容易出短脉冲。这三个参数基本决定了脉冲宽度和峰值功率的量级。
1.2 我用的速率方程形式与三项常见简化陷阱
不计空间分布、采用单模均匀加宽近似的速率方程，我习惯写成如下形式：
dφ/dt = φ(2σnl/tr - δ/tr) dn/dt = Rp - n/τf - σcnφ
其中，φ为腔内平均光子数密度（cm⁻³），n为反转粒子数密度（cm⁻³），tr = 2L/c为光子在腔内的往返时间，L是光学腔长，c是光速；σ是受激发射截面，l是增益介质长度；δ是谐振腔单程总损耗，包含输出镜透射损耗、内部散射吸收损耗和Q开关插入损耗；τf是上能级荧光寿命；Rp为泵浦速率，σcnφ项是受激辐射对反转粒子数的消耗。
这套方程看起来简单，真正动手时麻烦全在细节上。第一，受激辐射项的系数写法五花八门，有的教材在dn/dt里写-2σcnφ，有的写-σcnφ，原因在于对φ的定义不同，有的用腔内总光子数，有的用光子数密度，有的只算单向传播的光子。我建议选定一套符号体系后就始终如一，我的做法里：φ是腔内双向往返的平均光子数密度，dn/dt里受激辐射项系数为σcnφ，dφ/dt里的增益项用2σnl/tr，两者在能量上是自洽的。如果你看到别人代码里系数差2倍，不要慌，先确认他的φ是单向还是双向。
第二，δ不是输出镜透过率T那么简单。δ写成-0.5*ln(R1R2)再加内部损耗更准确，R1、R2是两个腔镜的反射率。当R接近1时，δ近似等于T/2加上内部损耗，但这个近似只在透过率很小时才够用。
第三，泵浦项Rp在单脉冲仿真里经常被人为删掉。如果只关注单个Q开关脉冲，泵浦阶段结束后的初始反转粒子数n0可以直接给定，脉冲过程中泵浦的影响很小，删除问题不大。但如果你想做重复频率仿真，脉冲与脉冲之间增益的恢复依赖Rp这项，删掉它整个脉冲串就做不出来。我会在后面的重复频率部分细说。
2. 用MATLAB求解Q开关脉冲：求解器选择和分阶段求解策略
2.1 为什么ode45会在脉冲上升沿翻车
很多人一上来就ode45，跑完发现要么极其缓慢，要么报错说计算失败。原因是Q开关脉冲的物理时间尺度跨度太大：泵浦阶段微秒到毫秒量级，脉冲建立和释放过程只有纳秒量级。举个例子，腔长20cm，光子在腔内往返时间约1.33ns，光子寿命τc大约十几纳秒，而脉冲上升沿可能只有几纳秒。ode45是显式非刚性求解器，在这种陡峭变化面前，为了让误差控制在容差范围内，步长会被压到极小，泵浦阶段几微秒就要算几十万步。
我的做法分两段：泵浦阶段能解析求解就不数值求解。若泵浦速率Rp恒定，dn/dt = Rp - n/τf这个方程有解析解，直接算脉冲触发时刻的n0即可，没必要让求解器耗在这里。真正需要数值求解的只有Q开关触发后的那一段，时间窗口从开关触发到脉冲结束，一般是几百纳秒到几微秒。这一段用ode15s或者ode23s这类刚性求解器，配合低容差，能稳定捕捉脉冲上升沿和下降沿。
2.2 最小可运行代码骨架与参数设置
下面这段是能跑通核心逻辑的骨架代码，参数以Nd:YAG为例，单位统一采用cm-g-s制，注意σ的单位是cm²，对应的长度单位必须是cm。如果改用国际单位制，σ和n的数值会差好几个数量级，很多人就是栽在这里。
% 参数定义（cm-g-s单位制） c = 3e10; % 光速, cm/s h = 6.626e-34; % 普朗克常数, J*s nu = c / 1.064e-4; % 1064nm频率, 1.064e-4 cm sigma = 2.8e-19; % 受激发射截面, cm^2 l = 0.6; % Nd:YAG晶体长度, cm L = 20; % 光学腔长, cm tr = 2*L/c; % 往返时间, s R1 = 1.0; % 全反镜反射率 R2 = 0.85; % 输出镜反射率 delta_i = 0.02; % 内部散射和插入损耗 delta = -0.5*log(R1*R2) + delta_i; % 单程总损耗 tau_c = tr / delta; % 光子寿命, s tau_f = 230e-6; % 荧光寿命, s % 泵浦阶段结束后初始反转粒子数密度 % 假设小信号增益系数 g0 = 0.3 cm^-1 g0 = 0.3; n0 = g0 / sigma; phi0 = 1e-3; % 自发辐射种子光子数密度, 不能太小 % 求解时间窗口：从Q开关触发到脉冲结束 t_span = [0, 5e-6]; % 5微秒足够 % 调用刚性求解器 opt = odeset('RelTol', 1e-8, 'AbsTol', 1e-12, ... 'Events', @(t,y) pulse_end_event(t,y)); p = struct('sigma',sigma,'l',l,'tr',tr,'tau_c',tau_c, ... 'tau_f',tau_f,'nu',nu); [t, y] = ode15s(@(t,y) qswitch_rhs(t,y,p), t_span, [phi0; n0], opt); phi = y(:,1); n = y(:,2);方程函数和事件函数写法如下：
function dydt = qswitch_rhs(~, y, p) phi = y(1); n = y(2); dydt = zeros(2,1); dydt(1) = phi * (2*p.sigma*p.l*n / p.tr - 1/p.tau_c); dydt(2) = -n/p.tau_f - p.sigma*3e10*n*phi; % 单脉冲内忽略泵浦 end function [value, isterminal, direction] = pulse_end_event(~, y) value = y(1) - 1e-4; % 光子数降到峰值千分之一以下时结束 isterminal = 1; % 终止求解 direction = -1; % 下降沿触发 end这里的Rp在单脉冲阶段我直接忽略了，因为脉冲过程只有几十纳秒到几微秒，泵浦项和自发辐射项的影响远小于受激辐射项。如果你要严格一点，把Rp - n/τf保留也行，但单脉冲仿真里差别很小。
2.3 事件函数和保证数值收敛的小技巧
事件函数是这类的核心技巧。脉冲释放完之后，光子数密度会指数衰减到极低，如果不加事件函数终止，求解器还会在这个“无意义”的长尾上耗费大量计算时间。我在事件函数里设置当光子数密度降到峰值的千分之一时停止，这样既能完整覆盖脉冲，又不会浪费算力。
另一个提高稳定性的技巧是合理设置种子光子数密度φ0。理论上Q开关脉冲的起点是自发辐射噪声，对应的光子数密度非常小，量级可能低到10⁻⁶甚至10⁻⁸。但数值上φ0取太小会有问题：求解器为了满足绝对容差，会把步长压到极小，甚至出现光子数为负的情况。我用1e-3到1e-2这个范围作为种子光子密度，不仅计算稳定，而且因为Q开关增益远大于损耗，种子光子的具体初始值对脉冲峰值和宽度影响很小。这个“种子不敏感”特性本身也是判断仿真是否正常的一个标志。
3. 仿完怎么看：用脉宽、峰值功率和提取效率交叉验证结果
3.1 三个“一眼判定”结果是否合理的量级标准
仿真跑通了不代表结果对，还要用物理量级来交叉验证。以Nd:YAG主动调Q为例，典型参数下我的判断标准有三个。
第一，脉冲宽度应该在纳秒量级，通常几个纳秒到几十纳秒。如果仿真出来是皮秒甚至飞秒，先检查是不是把腔长写成了毫米而不是厘米，或者把光子寿命τc算错了一个数量级；如果仿真出来是微秒级，大概率是初始反转粒子数n0太接近阈值，增益不够，脉冲建立过程被拉长。
第二，峰值功率在kW到MW量级。一个mJ级别的脉冲，脉宽10ns，对应的峰值功率就是100kW。仿真算出来的峰值功率如果到了GW量级，那就要怀疑是否把输出镜反射率R2写得太低，导致输出能量和脉冲宽度不匹配。
第三，能量提取效率不能超过初始储能。脉冲释放的能量来自反转粒子数的消耗，初始储能E_stored = hν * n0 * V，其中V是增益介质中光束的有效体积。输出脉冲能量与初始储能的比值就是能量提取效率，物理上必然小于1。如果仿出来的输出能量大于初始储能，说明速率方程里增益项或损耗项写错了，这是最硬性的检查。
3.2 从光子数到输出功率和脉冲能量的换算
速率方程直接解出来的是光子数密度φ，要跟实验对比还得换算成输出功率。换算关系可以这样理解：每个光子每经过一个往返时间tr，有一定概率从输出镜逃逸，这个概率等于输出镜透过率T（T=1-R2）。腔内光子总数为N_ph = φ * V，V是腔内模式体积，那么输出功率就是：
P_out(t) = N_ph(t) * hν * T / tr
用光子数密度表示就是：
P_out(t) = φ(t) * V * hν * T * c / (2L)
这个公式物理含义很清晰：腔内光子越多，输出镜透过率越高，腔长越短，输出功率越大。脉冲能量直接把P_out对时间积分即可。
这里要特别提醒，模式体积V一定要和n0对应的泵浦体积一致。很多人只算了速率方程，最后发现能量对不上，往往就是体积项没匹配。合理做法是用光束在增益介质内的有效截面积A乘上增益介质长度l，如果你用了高斯光束假设，有效截面积取π*w²，w是束腰半径。这样算出来的能量才是可以和实验对比的结果。
4. 我踩过的坑：从振荡发散到结果完全对不上
4.1 光子数初始值太小导致数值发散
我最初做仿真时，为了追求所谓“物理上准确的种子光子数”，把φ0取到1e-8量级。结果是ode15s在刚开始的几步就报错，或者干脆算出负的光子数。原因很简单：绝对容差AbsTol设成1e-12时，求解器会尝试把每个分量的误差控制在这个量级，而种子光子数密度本身就是1e-8，远大于容差，步长被迫缩到极短，数值误差被放大，最终导致振荡。
解决办法就是我前面提到的，把φ0抬高到1e-3左右。我自己测试过，φ0从1e-4到1e-1之间变化，脉冲峰值和宽度的差异不到0.1%，完全可以忽略。这是Q开关增益远大于损耗的特性决定的，激光器自己会把初始条件的微小差异“抹平”。如果你发现φ0对结果影响很大，那说明增益太接近阈值，仿真条件本身就需要重新审视。
4.2 损耗项符号和反射率换算
这个坑我记忆犹新。有一次我把速率方程里的δ符号写反了，写成dφ/dt = φ(2σnl/tr + δ/tr)，结果一运行，脉冲在t=0的时刻就直接爆发，完全没有积累过程。原因很明显，损耗项符号写反等于让谐振腔处于负损耗状态，光子数在没有初始反转粒子数的情况下也能指数增长。
更隐蔽的问题是损耗δ的换算。我第一次做的时候直接用δ = -ln(R2) 而不是 -0.5*ln(R1R2)，结果把损耗高估了一倍，仿真出来的脉冲宽度偏宽，峰值功率偏低。这个问题特别容易发生在电光调Q结构里，因为除了输出镜透过率，还要加上电光晶体的插入损耗和偏振片的损耗，如果不把这些损耗逐项列清楚，仿真结果和实验很难对上。我的建议是先把各项损耗列一张表：输出耦合损耗、散射损耗、插入损耗，再合成总的δ，不要嫌麻烦。
4.3 主动调Q开关时刻与脉冲建立时间
主动调Q和被动调Q在仿真上有个重要区别：主动调Q可以由外部信号控制开关时刻，但开关动作完成后，脉冲并不会立刻出现。从损耗降低到脉冲真正达到峰值，中间有一段建立时间，这个时间是触发信号到收到激光脉冲的延迟，在实验同步（比如触发示波器、同步泵浦源）中非常重要。
我见过不少同学把“开关触发时刻”误当成“脉冲峰值时刻”，在仿真里硬性把Q开关安排在预期输出时刻前，结果发现峰值总是往后偏。正确做法是先粗略跑一次仿真，量出建立时间，再反过来设置触发时刻。对于典型Nd:YAG主动调Q，建立时间从几十纳秒到几百纳秒不等，取决于初始反转粒子数超过阈值的程度。超阈值越高，建立时间越短，但脉冲宽度会更窄，这是调Q仿真中可以直接观察到的一个规律。
5. 让仿真贴近真实激光器：进阶修饰的几个方向
5.1 从单脉冲到重复频率脉冲串
单脉冲仿真跑通之后，很多人下一步会遇到重复频率问题。比如声光主动调Q激光器工作频率在1kHz到100kHz之间，这时候泵浦是连续的，每个脉冲消耗掉一部分反转粒子数，泵浦又在两个脉冲之间把反转粒子数重新“充满”。
这种场景下，不能再用“给定n0然后直接解脉冲”的单脉冲思路。我的做法是保持速率方程完整，把泵浦项Rp加回来：
dn/dt = Rp - n/τf - σcnφ
然后用循环求解。先令φ很小，求解一段时间（比如1/PRF），此时Q开关损耗高，腔内几乎不起振，反转粒子数从残余值回升；接着把损耗突然降低，求解脉冲；然后恢复高损耗，继续下一个周期。MATLAB里可以用事件函数在每个脉冲结束后重置积分区间，也可以简单用固定时间步长循环。注意保留Rp项之后，仿真时间步长覆盖范围更大，计算量会明显增加，建议泵浦阶段用大步长，只在脉冲阶段加密步长。
5.2 热效应、空间分布和电光开关过渡过程
如果想让仿真结果更贴近真实实验，还有几个修饰方向可以尝试。
热效应方面，高平均功率泵浦会让增益介质内部形成温度梯度，产生热透镜效应，等效于在谐振腔里加了一个焦距随泵浦功率变化的透镜。这会改变模式体积V和腔长L，进而影响τc和tr。最简做法是给δ加一个随泵浦功率变化的修正项，或者修正tr的数值，虽然粗糙，但能看出趋势。
空间分布方面，把常数反转粒子数n改为空间分布n(r)，泵浦光斑和振荡光斑都有横向分布，通常假设高斯分布。这样需要把速率方程在空间上离散化，计算量上一个台阶，但可以看到输出光束的空间分布和不同位置的粒子数消耗差异；对研究高功率激光器的热致双折射等问题，这一步是必要的。
电光调Q的开关过渡过程也值得提一句。Q开关的损耗不是理想阶跃变化的，KD*P电光晶体在加压和退压时都有亚纳秒到纳秒级的过渡时间。对脉宽几十纳秒的激光器，这个过渡时间影响很小，可以忽略；但对追求亚纳秒脉冲的高压调Q系统，过渡时间会直接影响脉冲建立时机和峰值功率。我当时用一个smoothstep函数代替阶跃函数表示开关损耗随时间的变化，对比下来确实发现脉宽有可测量的差异。
最后分享一个个人习惯。整套仿真环境稳定后，我会把泵浦功率、输出镜透过率、腔长这几个核心参数都做成可扫描的数组，一次性跑几十组，画出脉宽、峰值功率、脉冲能量随参数变化的曲线。这个操作在我后来设计实验时省了特别多时间，做实验前先在仿真里摸一遍趋势，到了实验台心里有底。另外，仿真文件里一定要把用的单位制和符号体系注释清楚，否则三个月后你自己回来看都会怀疑这里的系数到底是怎么来的。
本文还有配套的精品资源，点击获取
