简介本资源是一份面向光学工程、大气物理及自由空间光通信方向研究者与高年级本科生的仿真工具包聚焦高斯光束在大气湍流环境中的传输建模问题解决光强闪烁、相位畸变等关键效应的定量分析需求。压缩包仅含1个MATLAB脚本文件gauss.m大小1KB代码实现了基于Rytov近似与Kolmogorov湍流模型的光束传播仿真涵盖高斯光束参数初始化、Cn²谱设定、多层相位屏叠加、光强分布演化及闪烁指数计算等核心模块可直接运行并可视化光强衰减与相位扰动过程。目前已有1387人学习下载适用于快速复现大气湍流对激光传输影响的经典仿真流程为自由空间光通信系统设计、自适应光学预研及激光雷达误差建模提供轻量级可调试代码基础。1. 项目到底在算什么高斯光束、光强与大气湍流仿真的目标拆解1.1 一句话说清这个仿真在做什么大气里的激光传输听起来像是物理教科书上的内容但真正做起来你会发现它其实是一个先有数学假设、再做数值计算、最后用统计量说话的工程问题。我最近就在做一个这样的仿真一束波长为1.064微米的高斯光束从发射端出发穿过由大气湍流引起的折射率随机起伏区域经过10公里传播最后落到接收面上。整个过程要在计算机里还原出来同时把接收口径处的光强分布、光束扩展、能量集中度和随机抖动这些指标都量化出来。做过自由空间光通信、激光雷达或者激光定向能系统的人都会遇到同样的问题你没法在实验室里搭一条真正的10公里外场链路即便能搭大气条件也在不断变化没法重复实验。数值仿真就成了研究这类问题的最关键手段。所谓“高斯光束在大气湍流传输仿真”本质上就是用折射率随机场模拟大气的不均匀性再利用标量衍射理论把光场一帧一帧地“推”到目标距离。1.2 为什么要用数值仿真而不是直接套解析公式很多人第一反应是湍流里光传播不是有Rytov近似吗直接套闪烁指数的解析公式不就行了这个问题我在实际项目里被问过很多次。解析公式当然存在弱起伏条件下的平面波闪烁指数有一个非常简洁的表达式高斯光束也有对应的修正形式。但公式能算的东西非常有限它只能给出接收面上某个统计量的期望值给不出光强的二维空间分布更给不出单次实现的随机光斑图样。而工程上恰恰最关心后者。比如接收端用的自适应光学系统需要知道的就是某一帧波前的具体形状通信系统关心的则是某一瞬间落在探测器上的功率是深衰落还是深增强。这些信息只能通过数值仿真来获得。做法是跑蒙特卡罗仿真生成几百乃至上千个独立的湍流实现对每个实现都完成一次光场传播最后统计出光强闪烁、质心漂移、光束扩展等指标这才是实际工程里能直接指导设计的做法。1.3 仿真结果能解决哪些实际问题这个项目做完以后能输出的东西非常具体。第一个是给定发射参数下的接收面平均光强分布用来估计耦合进光纤的效率第二个是质心抖动的时间序列对应你跟瞄系统的跟踪误差预估第三个是闪烁指数的统计结果对应误码率的估算第四个是不同湍流强度下光束扩展的对比曲线。这些数据对链路预算、发射光学系统设计、接收口径尺寸选型都特别有价值。所以这个仿真不是一个纯理论作业它直接关系到一套系统的性能底线。接下来我会把我用的理论模型、实现流程以及调试过程中踩过的坑都整理出来给打算做类似工作的人一个可以直接上手的参考。2. 理论基础高斯光束与大气湍流里的几个关键物理量2.1 高斯光束的基本参数和传播公式先明确一件事高斯光束不代表光强是高斯分布就够了关键还是它复振幅的数学形式。在束腰处假设波前曲率半径为无穷大初始光场可以写成一个纯高斯振幅函数E0(r) sqrt(2P/(π·w0²)) · exp(-r²/w0²)其中P是总功率w0是束腰半径。这里的w0不是光强分布的标准差而是光强下降到峰值的1/e²处对应的半径。很多新手第一次就错在这里——用1/e的点来定义后面所有瑞利距离和发散角的计算都会偏。沿z传播以后理想情况下光场的解析表达式仍然保持高斯形式只是多了三个量光斑半径w(z)、波前曲率半径R(z)和古伊相移ζ(z)。瑞利距离zR π·w0²/λ决定了光束行为的转折点在z远小于zR时光束基本保持发射尺寸在z大于zR之后光束开始以线性发散角扩展。远场发散角大约是λ/(π·w0)这在设计链路的时候直接决定接收端光斑有多大。在例子里λ1064nm、w05cm时zR π×0.05²/1.064e-6 ≈ 7383m。10公里的链路距离是瑞利距离的1.35倍所以即使没有湍流理想光束到达接收面时光斑就已经从5cm扩展到了约8.2cm。这一点必须最先算清楚否则后面所有湍流导致的扩展你都会归因错误。参数数值说明波长 λ1064 nm常用近红外激光束腰半径 w05 cm发射端光斑半径瑞利距离 zR7383 m光束行为分界点链路距离 L10000 m约1.35倍zR无湍流光斑半径 w(L)约8.2 cm理想衍射叠加结果2.2 大气湍流的统计描述Cn²与Kolmogorov谱大气湍流对光的影响是折射率随机起伏引起的相位扰动在光学上可以等价成沿着路径存在很多随机薄相位屏。描述湍流强度的核心物理量是折射率结构常数Cn²单位是m^(-2/3)。在大气近地面、白天强日照条件下Cn²通常在10^-14到10^-13量级到了夜间或者高空可以降到10^-16甚至更低。这个量变化范围极大做链路预算时不能只取一个定值要按不同时段和不同高度分段处理。湍流的统计特性用折射率功率谱密度来描述。最经典的是Kolmogorov谱Φn(κ) 0.033·Cn²·κ^(-11/3)它假设湍流在惯性子区间内各向同性并且从内尺度l0到外尺度L0之间都满足幂律关系。实际数值仿真里我通常用von Kármán谱因为它在低频端有一个截止避免了功率谱在κ→0时发散。von Kármán谱的具体形式是Φn(κ) 0.033·Cn²·exp(-κ²/κm²) / (κ² κ0²)^(11/6)其中κm 5.92/l0κ0 2π/L0。这个修正对仿真特别重要因为相位屏的能量主要集中在低频部分如果你直接用Kolmogorov谱且网格尺寸不够大低频成分的截断会让仿真的光束漂移量严重失真。另外一个常用参数是弗里德相干长度r0它综合反映了路径上所有湍流强度对光学系统的影响r0 (0.423·k²·Cn²·L)^(-3/5)r0越小湍流越强。在λ1064nm、L10km下Cn²1e-15时r0约为1.3cmCn²5e-14时r0只有约1.1mm。这个量直接决定了仿真网格要取多细后面我会结合网格设置再讨论。2.3 从相位谱到相位屏为什么要转换成相位扰动光在湍流中的传播通常被处理成两个过程的循环真空中的自由衍射和薄层相位扰动。一个厚度为Δz的薄大气层对光场的相位调制可以在频域里写成相位功率谱Φφ(κ) 2π·k²·Δz·Φn(κ)这里的k 2π/λ是波数。有了这个相位谱就可以用随机数生成一个满足该谱统计特性的二维随机场这个随机场就是相位屏。数值上最常用的办法是频谱反演法在频域生成复高斯随机数乘上sqrt(Φφ(κ))再做逆傅里叶变换得到空间域的相位分布。这个方法实现简单、计算效率高但有一个已知问题——低频分量不足也就是说只靠原始网格分辨率生成的相位屏会在低频段明显偏离理论谱造成光束整体漂移量被严重低估。解决这个问题最常用的办法是次谐波补偿法也叫低频补偿法。具体做法是把频域网格的零频附近区域再细分比如分成3层或5层每层放大3倍或5倍把更精细的低频采样补进去。这样得到的相位屏既保留了高频细节又把低频大尺度相位起伏补全了。这个细节直接决定仿真结果的准确性我后面会再展开讲。3. 分步法仿真的实现3.1 第一步网格和物理参数初始化开始写代码之前先把所有物理参数和数值参数列清楚。物理参数包括波长λ、束腰半径w0、发射功率P、传输距离L、湍流强度Cn²、内尺度l0、外尺度L0。数值参数包括网格点数N、网格物理尺寸D、相位屏数量M。这几个数值参数的选取有讲究。网格尺寸D要至少是接收端光斑尺寸的3到4倍否则周期性边界效应会让光场从一侧绕到另一侧产生伪影网格点数N决定了空间采样间隔Δx D/N这个间隔要足够小能分辨湍流相位屏的细小结构和光场的干涉条纹相位屏数量M则跟传播步长有关每一步传播距离不能太大否则每一步内部的湍流起伏被平均掉仿真结果会过于乐观。以我的10公里链路为例D取0.5mN取512那么Δx≈0.98mm。这个分辨率对捕捉毫米级以上的湍流结构是够的但对Cn²5e-14的强湍流当r0降到1.1mm左右时这个分辨率就接近极限了。所以在强湍流条件下我会把网格提升到1024×1024或者减小单步Δz不让分辨率成为仿真的瓶颈。提示有一个经验法则网格间距Δx一般要小于0.5倍的弗里德相干长度r0才能比较可靠地还原强湍流下的相位细节。如果Δx比r0还大相位屏的高频信息就直接被网格抹掉了仿真结果的闪烁指数会比理论值偏小。3.2 第二步生成符合湍流统计的随机相位屏相位屏生成是整个仿真里最容易出错的一步。我建议按下面的流程走一遍。第一步生成频域网格频率坐标从-1/(2Δx)到1/(2Δx)步长1/(N·Δx)对二维网格计算每个频点的径向波数κ sqrt(κx² κy²)。第二步计算相位功率谱密度。按von Kármán谱计算Φφ(κ)内尺度和外尺度按实际情况给。注意κ0那个点分母里因为有κ0所以不会真的发散这个点对应的物理意义是无限大尺度的相位起伏实际中直接取零即可。第三步生成复高斯随机数。在频域生成一个复随机场每个点的实部和虚部都服从零均值、方差为1的高斯分布然后乘上sqrt(Φφ(κ))。这里的系数在各文献里经常差一个2π因子和你采用的谱密度定义有关最容易让人踩坑建议先用已知的理论闪烁指数去反验证一遍。第四步做逆傅里叶变换取实部作为相位屏。这里要注意由于傅里叶变换的周期性相位屏在空间上也是周期延拓的所以低频分量会缺失必须做次谐波补偿。次谐波补偿我实际用的是三层的做法。针对原始频域网格的零频周围区域分别用3倍、9倍、27倍细分的小网格重新采样低频相位谱再把结果叠加回原相位屏。这样原始网格最低频段0到1/L之间的信息就被补上了。核心部分可以用下面的伪代码表达# 频域随机场 random_field rng.normal(0, 1, (N, N)) 1j * rng.normal(0, 1, (N, N)) field_freq random_field * np.sqrt(Phi_phi) * delta_x # 逆傅里叶变换得到基础相位屏 phase_screen np.real(np.fft.ifft2(field_freq)) # 次谐波补偿对零频附近细分网格叠加低频分量 for n in range(1, num_subharmonics 1): scale 3**n df_sub 1 / (scale * N * delta_x) # 构建更细的低频网格计算对应相位谱 # 生成随机场叠加到 phase_screen 上注意次谐波补偿的层数不是越多越好。工程上3到5层、每层3到5倍细分已经足够覆盖低频端我实测下来5层比3层在光束漂移统计上的差别已经可以忽略但计算量却增加不少。3.3 第三步分步传播分步传播的实现代码结构很清晰核心就是一个循环。先把初始光场定义在网格上然后每一步执行两个操作自由空间衍射传播和相位屏调制。比较推荐的做法是“对称分裂步”即每一步先加半个相位屏然后是自由传播再加另半个相位屏。对称分裂步的精度比单边分裂步高在相位屏数量较多时误差可以控制在二阶小量几乎可以忽略。自由空间传播在频谱域的传递函数严格解是H(fx, fy) exp(i·k·Δz·sqrt(1 - (λfx)² - (λfy)²))。在小角度近似下可以用抛物近似H exp(-i·π·λ·Δz·(fx² fy²))。做过光场仿真的人都知道后者运算快、实现简单在光束发散角不大时精度也足够。但对于发散角非常大的系统或者传播距离特别长的情况要回到严格形式否则会引入明显的相位误差。我实际项目里的做法是先用角谱严格形式算一遍再用抛物近似算一遍对比接收面光强分布的差异。如果差异小于0.5%就可以放心用抛物近似如果差异大就必须用严格形式。这个验证流程花的时间不多但能避免后面所有结果都建立在有偏差的传播模型上。整个传播循环的骨架是这样的for i in range(M): # 先加半个相位屏 E E * np.exp(1j * phase_screen_half[i]) # 自由空间角谱传播 E np.fft.ifft2(np.fft.fft2(E) * H_prop) # 再加另外半个相位屏 E E * np.exp(1j * phase_screen_half[i])3.4 第四步输出指标计算仿真的最终输出不是光场本身而是从光场里提取的统计量。最常用的是以下几个。第一个是光强分布I |E|²。可以直接绘制二维热图也可以沿径向做平均光强曲线看是否还保持高斯轮廓。第二个是光斑半径w_eff。定义方法可以是二阶矩半径也可以按包含86.5%功率的半径来算。二阶矩半径对光强尾部比较敏感而86.5%功率半径更贴近工程上的“桶中功率”概念。两种都算对比着看更有意义。第三个是质心坐标与光束漂移。单次实现中光斑质心位置(x_c, y_c)可由光强加权计算x_c ∫x·I dxdy / ∫I dxdy。多次独立统计可以得到质心漂移的方差和标准差这个量直接对应你跟瞄系统的跟踪误差。第四个是闪烁指数σ_I²定义为σ_I² I²/I² - 1尖括号表示多次独立实现平均。它用来衡量接收功率的随机起伏强度对通信系统的误码率影响最大。输出指标定义工程意义光强分布 I|E|²接收面能量分布光斑半径 w_eff二阶矩或86.5%功率半径接收口径选型质心漂移光强加权质心的统计方差跟瞄系统误差闪烁指数 σ_I²I²/² - 1误码率估算4. 例程数据从弱湍流到强湍流的几个仿真案例4.1 无湍流情况下的自校验无论后边怎么加湍流第一步都必须先跑一遍无湍流传播。这一步有两层作用第一确认传播算法本身没写错第二得到理想情况下接收面上的光强分布作为后续湍流仿真的基准。我在10公里传输、λ1064nm、w05cm这个参数下的无湍流结果是接收端光斑半径约8.2cm跟理论值w(L) w0·sqrt(1 (L/zR)²)完全吻合。光斑中心峰值强度比发射端下降了约4.1倍这对应着同样的功率分布在更大的面积上。如果这一步出现明显偏差千万别急着加湍流先回头查传递函数的相位符号、网格物理尺寸和能量归一化。4.2 弱湍流Cn²1e-15的仿真结果取Cn²1e-15这是比较典型的夜间大气条件。用512×512网格、传播距离分成10段每段1000m跑200个独立实现。单次实现的接收面光强分布已经能看出明显的不规则光斑不再是平滑的高斯圆边缘出现随机畸变中心位置也有轻微偏移但整体上光斑仍然接近高斯形状。对200个实现做平均后结果是质心漂移的均方根约1.8mm闪烁指数约0.12平均光斑半径仍在8到9cm量级比理想值略有增加但增加幅度小于5%。这些数字都处在一个弱起伏区间说明在弱湍流下发射端的光束质量和接收端的性能主要由理想衍射和轻微湍流共同决定链路设计余量不需要太大。4.3 强湍流Cn²5e-14的仿真结果同样的参数换成白天强日照条件Cn²5e-14情况就完全不一样了。此时按前面公式估计r0≈1.1mm而我的网格间距Δx≈0.98mm刚好接近分辨极限。单次实现的接收面光强分布变得支离破碎光斑出现了多个亮斑和暗斑中心位置每帧都大幅跳动峰值光强也有很大起伏。200个独立实现统计下来质心漂移的均方根达到约25mm比弱湍流放大了一个数量级以上闪烁指数达到0.6以上平均光斑半径扩展到约15到18cm是理想值的两倍左右。这些数字对链路设计是很有冲击力的。比如通信接收端如果只用5cm口径的透镜耦合进光纤这个条件下光斑已经完全发散到没法有效耦合必须考虑使用大孔径接收、自适应光学补偿或者多孔径分集接收。湍流条件r0质心漂移RMS闪烁指数平均光斑半径无湍流—00约8.2cmCn²1e-15约1.3cm约1.8mm约0.12约8.5cmCn²5e-14约1.1mm约25mm约0.6约15-18cm4.4 统计结果怎么“看”才不踩坑有一点我必须提醒单次实现的光强分布虽然直观醒目但只能给你一个感觉不能作为设计依据。统计结果才是有意义的。我在项目里通常的做法是每次仿真跑200到500个独立实现每个实现用不同的随机种子最后统计平均。这里的关键是独立实现的样本量要大因为闪烁指数、质心漂移这类二阶统计量天然收敛慢样本太小时方差非常大你可能跑了好几百帧都看不出趋势。另外要盯住的一个坑是仿真中的光场能量守恒。在无吸收损耗的假设下每一步传播之后全场总能量应当保持恒定。如果发现能量在几个step后明显减小几乎可以肯定是FFT中引入了数值耗散或者相位屏边界处产生了大量数值衍射这时要优先检查网格尺寸和边界处理。5. 实际调试中最容易踩的坑5.1 网格间距与采样率在光学仿真圈子里采样问题大概是出现频率最高的错误来源。网格间距Δx过大会导致高频成分混叠低频成分被截断网格间距过小又会让计算量膨胀特别是N1024时每一步FFT都要耗掉相当多的内存和时间。我自己的经验是先根据光束尺寸和湍流的r0确定Δx再根据需要的视场确定N。具体来说至少要求Δx min(r0/2, λL/(2D))其中λL/(2D)是角谱传播的奈奎斯特条件防止高频分量出现周期性混叠。在强湍流、长距离链路里这个条件往往比r0条件更严格。5.2 相位屏周期性伪影相位屏是通过FFT生成的天然具有周期性自由空间传播用的角谱传递函数也是周期的。这两个周期性叠加在一起会让光场在边界处出现“卷绕”也就是光从左边出去、从右边绕回来。我在做10公里仿真时第一次就碰到了这个现象光斑周围出现了明显的高斯“鬼影”后来检查发现是网格尺寸D只取了光斑尺寸的1.5倍导致大角度衍射分量触到了边界。修正方法其实很简单把D放大到光斑尺寸的3倍以上。如果是强湍流导致光斑剧烈扩展还需要再加大。我一般会先用解析公式粗估接收端光斑尺寸再把它乘上3到4作为D。这个代价是网格点数可能翻倍但换来的是干净可靠的结果。5.3 随机种子与样本数量相位屏生成用的随机数要设置好随机种子。做蒙特卡罗统计时每个独立实现必须用不同种子做算法验证时固定的种子又特别方便因为同样条件下输出完全可以复现。我在代码里预留了一个seed参数验证和正式统计可以一键切换。另一个小建议是不要在循环里用系统时钟作种子否则并行计算时会出问题。样本数量方面我发现200个独立实现对于质心漂移的统计基本够用但闪烁指数的收敛明显更慢至少需要500个以上才能比较稳定地逼近理论值。如果时间紧张可以用一个“先快后慢”的策略先在网格数小一倍的快速模式下定级判断趋势再对关键工况用完整网格跑精细统计。5.4 和理论解析对比验证数值仿真最怕的事情是代码没写对但结果看起来还挺“合理”。唯一有效的校验方法就是拿已知的理论结果对标。我比较常用的校验标准有三个。第一无湍流时光斑半径和理论值误差应在1%以内。第二弱湍流时闪烁指数应与Rytov近似给出的解析值接近比如平面波弱湍流下的σ_I² 1.23·Cn²·k^(7/6)·L^(11/6)对高斯光束还要乘上对应的修正因子具体值可以查经典的Andrews著作。第三光束质心漂移的方差应随传播距离增加而增加而且最终结果不随网格尺寸变化而剧烈变化也就是要做网格收敛性测试。检查项方法通过标准无湍流光斑与理论公式对比误差小于1%弱湍流闪烁与Rytov近似对比相对误差小于10%网格收敛性改变网格数对比结果关键统计量变化小于5%能量守恒监控每一步总能量变化小于0.1%我曾遇到过这样一种情况结果完全符合理论解析公式但一改网格点数就剧烈变化。这说明结果高度依赖数值采样物理上并不可靠。遇到这种情况先降低传播步长、提高网格分辨率看结果是否收敛如果不收敛多半是相位屏生成或者传递函数出了问题必须先解决再谈下一步。6. 一些在实际项目中沉淀下来的经验6.1 参数文件与复现性管理做这类数值实验非常需要参数管理的习惯。每个仿真工况涉及的参数太多——波长、束腰、距离、Cn²、网格数、相位屏数、次谐波层数、随机种子——任何一个改动都可能导致结论变化。我现在的做法是把所有参数写在一个独立的参数文件里而不是散落在代码各处的赋值语句中。文件名带上工况摘要比如Cn15_L10km_w5cm_case1这样跑完之后回看结果永远知道当时用的是什么参数。这一步很土但对长周期项目来说节省的时间非常可观。6.2 关于代码性能的几点建议这个仿真看起来算法简单但性能瓶颈还是很明显的512×512×200次实现的传播循环即使只用抛物近似也要执行几十次FFT和IFFT。我一般会做两件优化。第一尽可能减少不必要的复数数组拷贝在循环里复用临时数组。第二把单个实现的传播循环封装成一个函数方便用并行循环或者GPU数组来加速。这套优化做完单次200实现的仿真时间能从十几分钟降到了两三分钟对于需要扫参数的场景特别有用。GPU加速方面市面上的开源光场仿真库已经不少也可以直接用支持GPU的FFT库。我实测下来GPU在本文还有配套的精品资源点击获取