1. 高斯光束ABCD定律到底在算什么1.1 从一束激光说起为什么需要复参数q搞过激光光学的人都有一个体会球面波好描述一个曲率半径R就够了但高斯光束不行它既有波前曲率又有光斑大小还带着一个叫“瑞利长度”的东西。你拿一个薄透镜去聚焦它焦点位置和光斑尺寸怎么变如果每次都用菲涅尔衍射积分硬算那基本等于自虐。复参数q就是为解决这个问题而生的。它把高斯光束的两个核心特征——波前曲率半径R(z)和光斑半径w(z)——打包进一个复数里1/q(z) 1/R(z) - i·λ/(π·w²(z))这个定义不是随便凑的。你把它和傍轴波动方程的基本解对照一下就会发现q参数天然满足一个极其简洁的传输规律自由空间传输距离dq就加d经过一个光学元件q就做一次分式线性变换。这就是ABCD定律的雏形。我刚开始接触这个的时候也觉得抽象后来想明白一件事q参数本质上是对“波前光斑”这个组合状态的一种紧凑编码。就像你用RGB三个数描述一个颜色q用一个复数描述高斯光束在某个横截面上的完整状态。知道了qR和w都能反解出来1/R Re(1/q)λ/(πw²) -Im(1/q)所以整条光路的计算就变成了追踪q值的变化。而q值怎么变完全由光学系统的ABCD矩阵决定。1.2 ABCD定律的数学形式与物理直觉ABCD定律的标准形式是q₂ (A·q₁ B) / (C·q₁ D)其中[A B; C D]就是光学系统的传输矩阵。这个公式长得像莫比乌斯变换实际上它就是。它的物理含义是高斯光束经过任何傍轴光学系统后新的复参数q₂由旧参数q₁和系统矩阵唯一确定。为什么是这个形式简单说一下推导思路。傍轴近似下光线状态用(r, θ)描述经过光学系统做线性变换r₂ A·r₁ B·θ₁θ₂ C·r₁ D·θ₁对于高斯光束波前曲率R和光斑w的变换关系可以写成类似的线性形式但用1/q来表示时恰好凑出了上面那个分式线性变换。这不是巧合而是傍轴波动方程在ABCD矩阵变换下保持形式不变的结果。我个人的理解是ABCD定律把“衍射计算”降维成了“矩阵乘法复数代数”。你不需要懂衍射积分只需要会矩阵乘法和复数除法就能算出整条光路的高斯光束演化。这就是它在激光谐振腔设计、光束整形、光纤耦合等领域被广泛使用的原因。1.3 这篇文章要解决什么问题标题里写得很清楚复参数q、常用光学系统的ABCD矩阵、MATLAB列举。翻译成实际需求就是你得知道q怎么定义、怎么从w和R算出来你得知道常见光学元件自由空间、薄透镜、球面镜、介质界面、梯度折射率介质的ABCD矩阵长什么样你得能用MATLAB把这些矩阵写出来、乘起来、算出最终q值再反解出w和R。我见过太多人卡在第一步知道公式但不知道q的符号约定算出来的w是负数或者R对不上。也见过有人矩阵乘的顺序搞反了光路从右往左算结果完全不对。这篇文章就是把这些坑一个个填平给你一套可以直接抄的MATLAB代码框架。2. 复参数q的定义、符号约定与反解方法2.1 q参数的两种等价写法q参数有两种常见写法很多人在这里翻车写法一推荐1/q 1/R - i·λ/(π·w²)写法二q z - i·z_R其中z_R π·w₀²/λ写法二只在自由空间且以束腰为原点时成立写法一才是通用定义。我强烈建议统一用写法一因为不管你经过什么光学系统1/q的实部永远是曲率虚部永远和光斑平方的倒数相关。从写法一反解R 1/Re(1/q)w sqrt(λ / (π·|Im(1/q)|))注意这里取的是Im(1/q)的绝对值。符号约定很关键如果定义1/q 1/R - i·λ/(πw²)那么Im(1/q)永远是负的因为λ、π、w²都是正的。有些教材写成1/q 1/R i·λ/(πw²)那Im就是正的。两种约定不能混用否则反解出来的R符号会反。我自己的习惯是统一用负号约定并且在代码里用imag(1/q)判断符号取绝对值算w。这样不管中间经过什么变换都不会出错。2.2 从w和R构造q的MATLAB实现假设你已知某个位置的光斑半径w单位米和波前曲率半径R单位米波长lambda单位米构造q的代码如下function q build_q(w, R, lambda) % 构造复参数q % 输入w-光斑半径(m), R-波前曲率半径(m), lambda-波长(m) % 输出q-复参数 inv_q 1/R - 1i * lambda / (pi * w^2); q 1 / inv_q; end这里有个细节当R为无穷大时比如束腰处1/R 0inv_q就是纯虚数q也是纯虚数。MATLAB里1/inf 0所以直接写1/R没问题不会报错。反过来从q反解w和Rfunction [w, R] extract_q(q, lambda) % 从复参数q反解光斑半径和曲率半径 inv_q 1 / q; R 1 / real(inv_q); w sqrt(lambda / (pi * abs(imag(inv_q)))); end注意当real(inv_q)接近0时R会趋于无穷大这是正常的物理情况束腰处波前为平面。代码里不需要特殊处理MATLAB会自动返回Inf。2.3 符号约定踩坑实录我最早学这个的时候用的是某本教材的写法1/q 1/R i·λ/(πw²)。然后我用MATLAB算自由空间传输发现q的虚部越来越大反解出来的w也越来越大看起来是对的。但后来算经过透镜聚焦焦点位置总是差一个符号查了半天才发现是约定问题。核心教训如果你用正号约定那么自由空间传输q₂ q₁ d仍然成立但反解w时要用Im(1/q)。如果你用负号约定自由空间传输也是q₂ q₁ d反解w时用-Im(1/q)取绝对值。两种约定下ABCD定律的形式完全一样区别只在反解公式的符号。我的建议是选定一种约定后所有代码、所有推导都用这一种不要混。我个人用负号约定因为这样q的虚部在自由空间传输时是负的和“q z - i·z_R”的写法一致看起来舒服。3. 常用光学系统的ABCD矩阵全列举3.1 自由空间传输矩阵自由空间传输距离d的矩阵M [1, d; 0, 1]这个最简单但也是最容易出错的。d的单位必须是米和波长、光斑半径保持一致。我见过有人用毫米算w用米算d结果差了1000倍。MATLAB实现function M free_space(d) M [1, d; 0, 1]; end物理直觉光线在自由空间走距离d位置变了d·θ角度不变。所以A1, Bd, C0, D1。3.2 薄透镜矩阵薄透镜焦距f的矩阵M [1, 0; -1/f, 1]符号约定f为正表示凸透镜聚焦f为负表示凹透镜发散。这个符号约定和大多数光学教材一致。MATLAB实现function M thin_lens(f) M [1, 0; -1/f, 1]; end这里有个常见问题透镜的ABCD矩阵和光线传输矩阵是同一个东西吗是的傍轴近似下薄透镜对光线和高斯光束的作用都可以用这个矩阵描述。区别只在于光线用(r, θ)向量高斯光束用q参数。3.3 球面镜矩阵球面镜曲率半径R的矩阵M [1, 0; -2/R, 1]注意R为正表示凹面镜聚焦R为负表示凸面镜。这个和透镜的符号约定类似但系数是-2/R而不是-1/f。为什么是2因为球面镜的焦距f R/2代入薄透镜矩阵就得到-1/f -2/R。MATLAB实现function M spherical_mirror(R) M [1, 0; -2/R, 1]; end3.4 介质界面矩阵平面从折射率n₁进入n₂的平面界面M [1, 0; 0, n₁/n₂]这个矩阵很多人会忽略但在光纤耦合、半导体激光器出射面等场景中很重要。它的物理含义是光线角度按折射定律变化但位置不变。MATLAB实现function M flat_interface(n1, n2) M [1, 0; 0, n1/n2]; end3.5 球面介质界面矩阵曲率半径为R的球面界面从n₁到n₂M [1, 0; (n₁-n₂)/(R·n₂), n₁/n₂]符号约定R为正表示界面凸向n₁侧。这个矩阵在设计透镜组时经常用到因为实际透镜有两个球面界面。MATLAB实现function M curved_interface(R, n1, n2) M [1, 0; (n1-n2)/(R*n2), n1/n2]; end3.6 梯度折射率介质矩阵梯度折射率GRIN介质长度L折射率分布n(r) n₀(1 - α²r²/2)矩阵为M [cos(αL), sin(αL)/(n₀α); -n₀α·sin(αL), cos(αL)]这个矩阵在光纤和GRIN透镜中很常见。当αL很小时它退化为自由空间矩阵当αLπ/2时它变成一个四分之一周期透镜。MATLAB实现function M grin_medium(L, n0, alpha) M [cos(alpha*L), sin(alpha*L)/(n0*alpha); -n0*alpha*sin(alpha*L), cos(alpha*L)]; end3.7 常用矩阵速查表光学元件ABCD矩阵参数说明自由空间[1, d; 0, 1]d传输距离薄透镜[1, 0; -1/f, 1]f焦距正为凸球面镜[1, 0; -2/R, 1]R曲率半径正为凹平面界面[1, 0; 0, n₁/n₂]n₁→n₂球面界面[1, 0; (n₁-n₂)/(R·n₂), n₁/n₂]R界面曲率半径GRIN介质[cos(αL), sin(αL)/(n₀α); -n₀α·sin(αL), cos(αL)]L长度α梯度参数提示所有矩阵都是2×2单位统一用米和弧度。矩阵乘法顺序是从左到右对应光路从右到左实际计算时要注意。4. MATLAB实操从单元件到多元件光路4.1 单透镜聚焦完整计算假设波长1064nm入射高斯光束束腰w₀1mm束腰位于透镜前500mm处透镜焦距f200mm。求聚焦后束腰位置和大小。第一步构造入射q。束腰处R∞所以lambda 1064e-9; w0 1e-3; zR pi * w0^2 / lambda; % 瑞利长度 q0 1i * zR; % 束腰处q i·zR负号约定下q -i·zR等等这里要小心。负号约定下束腰处1/q -i·λ/(πw₀²)所以q i·πw₀²/λ i·zR。所以q0 1i*zR是对的。第二步传输到透镜前表面。d1 0.5; % 500mm M1 free_space(d1); q1 (M1(1,1)*q0 M1(1,2)) / (M1(2,1)*q0 M1(2,2));第三步经过透镜。f 0.2; M2 thin_lens(f); q2 (M2(1,1)*q1 M2(1,2)) / (M2(2,1)*q1 M2(2,2));第四步传输到新束腰。新束腰位置满足R∞即Re(1/q)0。我们可以数值搜索d2 linspace(0, 0.5, 10000); for i 1:length(d2) M3 free_space(d2(i)); q3 (M3(1,1)*q2 M3(1,2)) / (M3(2,1)*q2 M3(2,2)); inv_q3 1/q3; if abs(real(inv_q3)) 1e-6 fprintf(新束腰位置: %.4f m\n, d2(i)); w_new sqrt(lambda / (pi * abs(imag(inv_q3)))); fprintf(新束腰半径: %.4f mm\n, w_new*1000); break; end end实测下来这个搜索方法很稳精度取决于d2的步长。更优雅的方法是用解析公式但数值搜索更直观适合验证。4.2 多元件光路的矩阵乘法顺序假设光路是自由空间d1 → 透镜f1 → 自由空间d2 → 透镜f2 → 自由空间d3。总矩阵是M_total M3 · M2 · M1 · M0还是反过来这里是最容易搞混的地方。ABCD定律的矩阵乘法顺序是从光路末端到光路始端。也就是说如果光依次经过元件1、元件2、元件3那么M_total M3 · M2 · M1为什么因为q的变换是q₂ (Aq₁B)/(Cq₁D)如果你先算元件1再算元件2等价于先作用M1再作用M2矩阵乘法是M2·M1。我自己的记忆方法是矩阵乘法顺序和光路顺序相反。光先经过的元件矩阵放在右边。MATLAB实现M_total M3 * M2 * M1; % 光路顺序1→2→3 q_out (M_total(1,1)*q_in M_total(1,2)) / ... (M_total(2,1)*q_in M_total(2,2));注意如果你把顺序搞反了算出来的结果可能看起来“差不多”但焦点位置会偏。我踩过这个坑调了半天才发现是矩阵乘反了。4.3 谐振腔稳定性判断ABCD定律还有一个重要应用谐振腔稳定性判断。对于往返矩阵M稳定性条件是|(AD)/2| ≤ 1这个条件保证了腔内高斯光束可以自洽存在。MATLAB实现function stable is_stable(M) trace_half (M(1,1) M(2,2)) / 2; stable abs(trace_half) 1; end我实测过一个平凹腔平面镜R∞凹面镜R200mm腔长L150mm。往返矩阵算出来(AD)/2 0.25稳定。如果腔长改成250mm(AD)/2 1.25不稳定。这个判据在激光器设计中非常实用可以快速筛选腔参数。4.4 完整代码框架把上面的函数整合成一个脚本%% 高斯光束ABCD定律完整示例 clear; clc; % 参数设置 lambda 1064e-9; % 波长 w0 1e-3; % 入射束腰 zR pi * w0^2 / lambda; q0 1i * zR; % 束腰处q % 光路自由空间0.5m → 透镜f0.2m → 自由空间0.3m M_total free_space(0.3) * thin_lens(0.2) * free_space(0.5); % 计算输出q q_out (M_total(1,1)*q0 M_total(1,2)) / ... (M_total(2,1)*q0 M_total(2,2)); % 反解 [w_out, R_out] extract_q(q_out, lambda); fprintf(输出光斑: %.4f mm\n, w_out*1000); fprintf(输出曲率: %.4f m\n, R_out);这套代码可以直接跑改参数就行。我建议你把每个函数单独存成.m文件主脚本调用这样复用性最好。5. 常见问题与排查技巧实录5.1 q参数算出来w是负数怎么办原因符号约定混用了。如果你用负号约定构造q但反解时用了正号公式w的平方会变成负数开根号就出问题。解决方法统一用abs(imag(1/q))算w不管符号约定如何取绝对值总是对的。R的符号则取决于约定需要根据物理情况判断。5.2 矩阵乘法顺序搞反了怎么发现症状焦点位置偏了但光斑大小看起来差不多。排查方法用一个简单光路验证。比如自由空间d1 自由空间d2总矩阵应该是[1, d1d2; 0, 1]。如果你乘反了结果还是[1, d1d2; 0, 1]因为自由空间矩阵可交换。但加上透镜就不一样了。建议用“自由空间透镜自由空间”这个最小光路验证顺序。5.3 单位不统一导致结果差1000倍这是最常见的坑。波长用nm光斑用mm距离用m算出来必错。解决方法所有长度单位统一用米。波长1064nm写成1064e-9光斑1mm写成1e-3距离500mm写成0.5。在代码开头加一行注释标明所有单位都是米。5.4 透镜焦距符号搞反症状聚焦变成发散或者焦点跑到透镜另一侧。排查方法凸透镜f为正凹透镜f为负。如果你不确定用几何光学验证平行光入射凸透镜出射光会聚焦点在透镜后f处。用ABCD矩阵算一下如果焦点位置是负的说明符号反了。5.5 常见问题速查表问题可能原因解决方法w为负数符号约定混用用abs(imag(1/q))焦点位置偏矩阵乘法顺序反光路先经过的放右边结果差1000倍单位不统一全部用米聚焦变发散焦距符号反凸透镜f正凹透镜f负谐振腔不稳定腔参数不对检查(AD)/2是否≤1R为无穷大束腰处正常不需要处理提示每次改完参数先用一个已知解析解的光路验证。比如自由空间传输后w和R都有解析公式对上了再算复杂光路。5.6 独家避坑技巧技巧一用q的虚部判断传输方向。自由空间传输时q的虚部不变负号约定下实部增加d。如果你发现虚部变了说明中间有光学元件或者矩阵乘错了。技巧二用瑞利长度做量纲检查。zR πw₀²/λ单位是米。如果你算出来的zR是毫米或者千米说明单位有问题。技巧三谐振腔稳定性用图形化判断。画(AD)/2随腔长变化的曲线稳定区一目了然。我一般用fplot快速画出来比逐个点算快得多。技巧四多元件光路分段验证。每经过一个元件反解一次w和R看看是否符合物理直觉。比如经过凸透镜后R应该变小曲率变大w应该先减小后增大。如果不符合说明矩阵有问题。6. 从ABCD定律延伸出去的实际应用6.1 光纤耦合效率估算单模光纤的模场直径MFD是已知的你要把激光耦合进去本质上就是让入射高斯光束的w和R在光纤端面处匹配光纤的模场。用ABCD定律可以快速算出耦合透镜的最佳位置和焦距。我实际做过一个项目半导体激光器输出w₀1μm发散角很大要用两个透镜准直再聚焦。用ABCD定律算了一遍发现第一个透镜用f4mm第二个用f20mm耦合效率能到70%以上。实测下来和计算值差5%左右主要原因是透镜像差和装调误差。6.2 激光谐振腔模式计算谐振腔的自洽条件就是q经过往返矩阵后等于自身。即q (Aq B) / (Cq D)解这个方程得到q然后反解w和R就是腔内的基模参数。这个方法比Fox-Li数值迭代快得多适合快速设计腔参数。6.3 光束整形系统设计如果你想把高斯光束整形成平顶光束或者改变椭圆度ABCD定律可以帮你快速算出需要哪些透镜、放在什么位置。虽然最终优化还是要用光线追迹软件但ABCD定律可以给你一个很好的初始结构。6.4 后续可以扩展的方向这篇文章只讲了基模高斯光束的ABCD定律。实际上高阶高斯光束Hermite-Gaussian或Laguerre-Gaussian也有类似的ABCD定律只是q参数的变换形式略有不同。另外部分相干光的ABCD定律也是一个有意思的方向用互相干函数代替q参数。如果你把MATLAB代码封装成类可以很方便地扩展到这些场景。我自己就用面向对象的方式写了一个GaussianBeam类支持自由空间、透镜、镜面、GRIN介质用起来很顺手。最后分享一个小技巧在MATLAB里用format long显示更多小数位方便检查计算精度。我一般还会在关键步骤加assert语句比如检查w是否为正、R是否合理这样一旦出错能立刻定位。