简介本资源是一份面向光学工程、大气物理及激光通信领域初学者与科研人员的MATLAB仿真脚本聚焦高斯光束在真实大气湍流环境下的传输特性建模重点解决光强闪烁与相位畸变的定量分析问题。适用于自由空间光通信系统设计、激光雷达性能预估、天文自适应光学算法验证等实际场景。压缩包仅含1个核心文件gauss.m1KB为轻量级可执行脚本完整实现了高斯光束初始场构建、基于Kolmogorov谱的大气湍流相位屏生成、多层湍流传输的傅里叶光传播模拟以及光强分布演化与闪烁指数计算等关键功能。已有1387人学习下载用户可直接运行该脚本复现典型湍流条件下的光束展宽、散斑形成与强度起伏现象并基于源码快速修改参数如波长、束腰半径、结构常数Cn²开展对比实验是理解Rytov近似与湍流光学效应的高效入门工具。 搞激光通信或者激光雷达外场实验的人多半见过这个场面接收屏上的光斑一会儿聚成一个亮芯一会儿弥散成一大团模糊亮斑能量读数像抽风一样上下乱跳。这台“看不见的手”就是大气湍流。我去年在做一套1550nm自由空间激光链路的预研时想在开实验之前就知道一束高斯光束经过1公里大气后光强会演变成什么样子于是花了两周时间搭了一套“高斯光束在大气湍流中的传输仿真”。这套仿真后来成了链路预算、接收口径选取的重要参考也让我把相位屏、闪烁指数这些概念从“书上见过”变成了“心里有数”。今天把这套仿真从物理模型到代码实现完整写一遍适合正在做大气光学仿真、激光传输评估的研究生和工程师参考。1. 湍流大气里的高斯光束仿真到底在算一个什么问题1.1 激光在大气中会经历哪些事为什么湍流最难缠激光在大气里传播通常会同时碰到三件事吸收、散射和折射率随机起伏。吸收和散射是大气分子和气溶胶干的活儿它们的效应相对稳定可以用大气透过率模型算个大概。真正让人头疼的是第三件——大气湍流引起的折射率随机起伏。地面上方的空气并不是均匀介质。太阳加热地面热空气团不断上升与冷空气混合形成了大大小小的漩涡。这些漩涡的温度和密度跟周围不一样折射率也就有细微差异。光波穿过这些大小不一的“折射率斑块”时不同位置的相位被随机调制光波前变得歪歪扭扭。等光继续往前传播这种相位畸变会通过衍射演变成光强起伏、光束扩展、光束漂移和到达角抖动这就是实验里看到的那些现象。这套仿真要做的事就是把这一串物理过程数字化从光源发出一束高斯光束让它在带有湍流扰动的大气中传播一段距离最后在接收面上统计光强的分布和抖动程度。1.2 高斯光束为什么是光学仿真的“标准起点”实际激光器输出的光斑并不是理想均匀平面波更接近基模高斯光束——光强从中心向外按高斯函数衰减。选择高斯光束作为仿真对象有三个很实际的原因绝大部分单模光纤激光器、半导体激光器经过准直后的输出都可近似为高斯光束仿真结果能直接对应工程问题。高斯光束有解析的真空传播公式方便做无湍流情况下的自检。很多更复杂的光束如平顶光束、涡旋光束都可以看作高斯光束的变形或叠加先把高斯情况吃透后续扩展就有基础。高斯光束的复振幅可写成U0(r) exp(-r^2 / w0^2)其中 r 是离光轴的距离w0 是束腰半径。在束腰处等相位面是平面所以这里不需要加二次相位因子。波长、束腰、传播距离这几个参数共同决定光束的衍射特性。1.3 仿真的输出指标平均光强、闪烁指数、光束漂移与扩展仿真结束之后我们关心的不是某一次随机实现长什么样而是一组统计指标。我实际用到的有四个平均光强分布多次独立湍流实现下接收面光强系综平均。它反映光斑在长时间曝光下的整体轮廓能看到中心能量被摊开、峰值下降。闪烁指数σI^2 I^2/^2 - 1。这个指标量化光强抖动的剧烈程度是激光通信里衡量误码率的重要参数。光束质心漂移光强加权质心在接收面上的随机偏移量统计其标准差。它反映湍流大尺度涡旋对光束整体的“掰弯”作用。二阶矩光束半径按光强的二阶矩计算等效光斑半径用来量化光束扩展量评估接收端探测器口径够不够。这四个指标从不同角度描述一束光被湍流“折磨”后的状态也是仿真是否可信的关键判据。2. 相位屏法把一整条大气路径拆成一摞“随机透镜”2.1 薄屏近似的物理依据传播与扰动可以分离大气湍流是连续分布在整条路径上的严格求解光波在其中传播需要解随机介质的波动方程这几乎不可能在普通计算机上完成。工程上常用的处理方法是“相位屏法”把连续的随机介质离散成一系列薄相位屏。物理依据是这样的在厚度足够小的一段路径上可以近似认为湍流只改变光波的相位不改变振幅而在这段路径之外光波按真空衍射规律传播。于是整条大气路径被切成一张张“随机透镜”光每走一段就先在真空中衍射传播再被叠加上一个随机相位扰动。这张“透镜”不会真的聚焦或发散光束它只是让波前变得高低不平。需要强调的是薄屏近似成立的条件是单段路径长度 dz 远小于湍流外尺度 L0同时整个路径上的湍流起伏不太强。在仿真中一般取 dz 在50米到200米之间路径总长1公里取10~20个屏足够。2.2 折射率功率谱怎么选Kolmogorov谱与Von Karman谱相位屏的“随机”不是瞎随机的它在频域上必须满足大气湍流的功率谱特性。湍流的折射率起伏常用功率谱密度来描述最经典的是Kolmogorov谱但它在低频和高频两端都有发散问题实际仿真里几乎都用修正后的Modified Von Karman谱Φn(kappa) 0.033 Cn^2 · (kappa^2 k0^2)^(-11/6) · exp(-kappa^2 / km^2)其中 kappa 是空间频率角波数Cn^2 是折射率结构常数是湍流强弱的度量k0 2π/L0L0 是湍流外尺度km 5.92/l0l0 是湍流内尺度。外尺度 L0 决定了湍流最大漩涡的尺寸通常取几米到几十米内尺度 l0 对应能量耗散的尺度典型值在毫米到厘米量级。两者不直接参与传播计算但决定了相位屏频谱的形状和能量分布。两种谱的选择可以这样理解功率谱类型适用范围问题Kolmogorov谱理论分析中频段低频和高频发散仿真中直接使用会导致能量异常Von Karman谱实际数值仿真需要指定L0和l0但更符合物理Modified Von Karman谱工程仿真通用在Von Karman基础上加高频截断项数值稳定我一般直接取 L020m、l01cm 做近地面水平链路。如果仿真的是高空或星地链路L0会更大l0也可能不同需要按实际场景调整。2.3 谱反演生成相位屏的完整步骤与次谐波补偿用谱反演法生成相位屏的思路很直接在频域构造一个满足湍流相位谱的随机复数场再傅里叶逆变换回空间域。具体可以写成function phz phase_screen_fft(N, dx, lambda, dz, Cn2, L0, l0) % 谱反演法生成一个湍流相位屏不含次谐波补偿 % N : 网格点数 % dx : 网格间距 (m) % lambda : 波长 (m) % dz : 单层湍流厚度 (m) % Cn2 : 折射率结构常数 (m^(-2/3)) % L0 : 外尺度 (m) % l0 : 内尺度 (m) fx (-N/2 : N/2-1) / (N*dx); [fx, fy] meshgrid(fx); kappa sqrt(fx.^2 fy.^2); kappa(kappa 0) 1e-10; k0 2*pi/L0; km 5.92/l0; PS 0.033 * Cn2 .* exp(-kappa.^2 / km^2) ./ (kappa.^2 k0^2).^(11/6); S_phi 2*pi * (2*pi/lambda)^2 * dz .* PS; cn (randn(N) 1i*randn(N)) .* sqrt(S_phi) * (N*dx)^2; phz real(ifft2(ifftshift(cn))); end这里有一个新手很容易踩的坑直接生成的相位屏低频分量严重不足因为频谱在低频端能量高但FFT网格上低频采样点太少导致大尺度涡旋的起伏被“饿死”。表现为仿真出的光束漂移量远小于理论值。解决办法是叠加“次谐波补偿”在低频区域用更细的网格补充采样。我一般叠加三层每层将频域采样间隔缩小3倍。这一步加上之后光斑漂移量才跟理论对得上。业内常说“不加次谐波的相位屏只有高频细节没有大尺度骨架”就是这个道理。3. 分步传播的数值模型角谱法为什么是首选3.1 傅里叶光学视角下的单步传播有了相位屏接下来要解决的是光波如何在两个相位屏之间传播。在真空里光波传播可以看作一个线性系统用角谱法处理非常顺手。角谱法的核心思想把光场分解成一系列朝不同方向传播的平面波分量每个分量在空间传播一段距离后只改变相位不改变振幅。数学上就是在频域乘一个传递函数H(fx, fy) exp(i · 2π · z · sqrt(1/λ^2 - fx^2 - fy^2))具体实现时光场先做二维FFT变换到频域乘上传递函数再做逆FFT回到空间域。这套操作每个传播步要算两次FFT对N512的网格来说速度很快。3.2 角谱传播与菲涅尔近似的取舍分步传播也可以用菲涅尔衍射积分它在近轴条件下成立形式更简单FFT实现也快。但菲涅尔近似要求传播距离满足一定条件而且当网格尺寸和传播距离的搭配不理想时容易出现采样问题。我做对比之后固定用角谱法原因有两个。第一角谱法没有近轴限制对某些大发散角的光束场景也能处理。第二角谱法对传播步长的适应范围更宽不容易因为 dz 取得不合适而出现数值震荡。实际仿真中角谱传播的代码是这样一段% 角谱法传播一步 function U ang_spectrum_prop(U, L, lambda, z) N size(U, 1); dx L / N; fx (-N/2 : N/2-1) / (N*dx); [FX, FY] meshgrid(fx); KX 2*pi*FX; KY 2*pi*FY; k 2*pi/lambda; Kz sqrt(k^2 - KX.^2 - KY.^2); H exp(1i * Kz * z); UF fftshift(fft2(U)); UF UF .* H; U ifft2(ifftshift(UF)); end这里注意 fftshift 和 ifftshift 的配对使用。频域坐标的原点在数组中心fft2 之后直接乘传递函数会因为频域偏移不对而产生整体平移的伪影很多人第一次跑出来光斑整体错位就是忘了处理这个。3.3 网格、采样点数与传播步长的匹配规则分步传播看着简单真正决定仿真质量的往往是网格参数。我自己调试时总结了几条规则网格总宽度 L 必须大于光束到达接收面时的最大扩展宽度否则光场被边界截断会产生严重的衍射伪影。空间采样间隔 dx 要满足奈奎斯特条件保证高频分量不混叠。dx 太小会增大FFT计算量太大则丢失高频成分。相位屏数量要足够通常每个屏之间距离不超过200米。也不能一味增加屏数每一层都会引入相位噪声屏太多反而会让统计结果偏向“过湍流化”。拿我常用的例子来看波长1550nm束腰2cm总距离1公里束腰处的瑞利距离大概在811米所以1公里处的真空光斑半径约3.2cm不到0.1米。取网格宽度 L0.2mN512dx约0.39mm满足采样要求同时保留了足够大的边界余量。如果你仿真的距离达到5公里以上光斑扩散可能到几十厘米量级网格宽度要相应加大。经验是先跑一次无湍流情况看光斑在接收面的理论半径然后让 L 至少是光斑半径的5倍。4. 主程序实现与光强统计指标的计算细节4.1 仿真实例参数设置波长1550nm、1公里链路把前面的模块拼起来一个完整的仿真主程序如下。这套参数对应近地面1公里水平链路Cn2取典型的白天值1e-15 m^(-2/3)湍流中等偏强。% 主程序高斯光束经过湍流大气传输1km clear; clc; % 基本参数 lambda 1550e-9; % 波长 1550nm w0 0.02; % 束腰半径 2cm L 0.2; % 网格总宽度 0.2m N 512; % 网格点数 dx L / N; % 空间采样间隔 k 2*pi/lambda; % 湍流参数 Cn2 1e-15; % 折射率结构常数 L0 20; % 外尺度 20m l0 0.01; % 内尺度 1cm z_total 1000; % 总传播距离 1km nscr 10; % 相位屏数量 dz z_total / nscr; % 单步距离 % 初始高斯光束 [x, y] meshgrid((-N/2 : N/2-1) * dx); r2 x.^2 y.^2; U0 exp(-r2 / w0^2); % 蒙特卡洛次数 nMC 50; % 统计量累积 I_accum zeros(N, N); I2_accum 0; I_mean_accum 0; cx_accum 0; cy_accum 0; W2_accum 0; for mc 1:nMC U U0; for i 1:nscr % 真空传播 U ang_spectrum_prop(U, L, lambda, dz); % 叠加相位屏 phz phase_screen_fft(N, dx, lambda, dz, Cn2, L0, l0); U U .* exp(1i * phz); end % 接收面强度 I abs(U).^2; I_accum I_accum I / nMC; Itot sum(I(:)) * dx^2; % 总功率 I_mean_accum I_mean_accum Itot / nMC; I2_accum I2_accum (sum(I(:).^2) * dx^2) / nMC; % 质心 cx_accum cx_accum (sum(sum(I.*x)) * dx^2) / Itot / nMC; cy_accum cy_accum (sum(sum(I.*y)) * dx^2) / Itot / nMC; % 二阶矩半径平方 W2_accum W2_accum (2 * sum(sum(I.*r2)) * dx^2) / Itot / nMC; end % 输出统计指标 I_mean I_accum; % 平均光强分布 sigma_I2 I2_accum / (I_mean_accum^2) - 1; % 闪烁指数 W_rms sqrt(W2_accum); % 等效光斑半径这段代码可以直接跑注意 nMC 不要取太少否则统计误差很大。我自己一般至少跑30次以上看闪烁指数波动趋于平稳再加次数。4.2 相位屏生成与传播主循环的常见实现顺序问题相位屏和传播循环的配合顺序看起来是小细节但很影响结果。我见过不少初学朋友把相位屏放在传播之前即先乘屏再传播。这样做不是完全不行但物理含义有偏差。更合理的做法是每一段路径先按真空衍射传播到这段的终点再叠加这个位置处的相位屏。也就是“先走一段再被打一巴掌”。这样每个相位屏都对应一段大气路径的累积效应和薄屏近似的物理图像一致。还有一个容易被忽略的点每个相位屏在每次蒙特卡洛运行里必须重新生成否则整个统计结果就只有一种湍流实现算出来的闪烁指数会严重偏离真实值。相位屏生成函数里的 randn 每次调用都会产生不同随机数所以只要把 phz 放进循环内就行。4.3 光强统计指标的计算闪烁指数、质心漂移、光斑半径统计指标的计算藏在主循环里这里单独提出来解释因为很多朋友拿到了正确光强分布却输出了错误的指标。闪烁指数 σI^2 I^2/^2 - 1这里的尖括号表示系综平均实际操作中对多次蒙特卡洛结果做平均。要注意的是I 是接收面上所有像素的光强值加起来后的整体统计还是取某个特定点的光强随时间起伏两者含义不同。我这里的实现用的是接收面总功率的起伏对应通信系统里“接收总能量抖动”的场景。如果研究的是单个探测器小孔处的闪烁应该固定一个空间位置例如光斑中心收集这个点在多次湍流实现下的光强再计算闪烁指数。这个区别在写论文时容易被审稿人追问务必弄清楚。质心漂移和光斑半径用的是光强加权矩。质心就是光强的一阶矩位置每次随机实现质心都会在接收面上抖动光斑半径用的是二阶矩定义的 rms 半径公式里的系数2对应高斯光束的二阶矩半径约等于 1/e^2 强度半径这是光学里的通用约定不要漏掉。5. 仿真自检与常见陷阱怎么确认结果没跑偏5.1 无湍流时的衍射极限自检第一次跑完仿真先别急着上湍流。把 Cn2 设成0跑一遍看结果能不能对得上高斯光束的真空衍射公式。这是成本最低的自检手段。高斯光束在真空里传播距离 z 后光斑半径满足w(z) w0 · sqrt(1 (z / zR)^2)其中 zR π·w0^2/λ 是瑞利距离。按上面那组参数zR约811米1公里处的理论光斑半径约3.17cm。在仿真里算出的 W_rms 应该和这个值接近误差通常在1%以内。如果偏差超过5%先检查网格总宽度是不是太小导致边界截断再检查角谱传递函数的频域坐标是否正确。我遇到过一种情况传播步数多了之后光斑半径逐次增加明显不符合衍射规律。后来排查发现是角谱法中某一步把 ifftshift 写成了 fftshift导致每个传播步都叠加了一个微小平移积累起来就让光斑看起来在“膨胀”。这种问题只有无湍流自检才能暴露出来。5.2 弱湍流条件下与Rytov理论的对比自检完无湍流情况就把湍流打开但此时选择弱湍流参数Cn2取1e-16左右与Rytov理论对比。弱起伏条件下平面波的闪烁指数约等于Rytov方差σR^2 1.23 · Cn2 · k^(7/6) · z^(11/6)高斯光束的闪烁指数比平面波略复杂但数量级应该与 σR^2 接近。仿真中取同样的 Cn2、1公里路径σR^2 算出来约0.34。仿真得到的 σI^2 落在0.2~0.5之间就说明相位屏的幅度和传播算法基本正确。如果算出来的闪烁指数明显偏大比如超过1多半是相位屏数量不足单层屏携带的能量太大等效于把湍流“压扁”成了强扰动偏离了弱起伏假设。这时候增加 nscr 到20或30闪烁指数会降下来并趋于一个稳定值。5.3 常见的三个坑混叠、次谐波缺失、相位屏尺寸不足最后集中说三个我在调试过程中踩得最深的坑。第一个是频域混叠。角谱法里频域传递函数 H 对高频分量有截断效应如果网格采样间隔 dx 太大高频分量会混叠回低频产生奇怪的周期性格子条纹。判断方法是看接收面光强分布边缘是否出现规则的条纹伪影。解决办法是减小 dx也就是在总宽度不变的前提下增大 N。但 N 太大会显著增加计算时间512×512点的单次蒙特卡洛已经要跑十几秒更大就要考虑GPU加速了。第二个是次谐波缺失导致的漂移量偏小。这个前面提过症状是仿真光斑在接收面上基本不“乱飘”质心标准差远小于理论估算。解法就是加次谐波补偿。如果你用现成的相位屏函数最好确认一下它是否包含次谐波。第三个是相位屏尺寸不足。相位屏的物理尺寸就是网格总宽度 L它必须覆盖光束经过湍流后可能偏转到的范围。强湍流下光束漂移可以达到厘米级甚至更高如果 L 只比光斑大一点点统计结果会因为“光束撞墙”而失真。我习惯在强湍流下把 L 调大到0.5m以上再跑一次对比两次结果如果统计指标变化超过10%说明原网格宽度不够。最后一个个人经验仿真参数不是越多越好。相位屏数量、采样点数、蒙特卡洛次数这三者的合理组合是靠“先跑一遍对照理论值”来确定的。每换一组大气条件都应该重新做一次5.1和5.2的验证。等你有了一套自己的验证流程这套仿真就能很轻松地扩展到不同波长、不同距离、不同湍流强度的场景甚至为后续做自适应光学校正、多光束合成提供仿真基础。说到底数值仿真不是追求代码好看而是要让结果能在实验之前帮你做出靠谱的判断。本文还有配套的精品资源点击获取