简介这是一份面向光学与激光通信研究的高斯光束大气湍流传播仿真MATLAB源码聚焦光斑仿真、大气湍流屏建模与大气传输效应分析适合科研人员、工程师及光学相关专业学生使用用于理解湍流对光束质量的影响并辅助后续补偿设计。压缩包仅含1个m源文件整包仅1KB代码精简无需额外依赖即可在MATLAB中直接运行便于快速复现和参数修改资源紧密围绕大气湍流中光束传播的核心问题利用随机相位屏模拟大气折射率起伏可设定不同湍流强度与传输距离计算接收面光强分布直观展示光斑扩展、扭曲与闪烁等现象。目前已有502人学习下载。通过运行gauss.m可评估湍流引起的光束质量退化量化相位扰动导致的波前畸变为激光通信、天文观测等提供理论依据与数据支持相关专业学生可通过修改代码参数深入探究不同湍流条件下的光斑演化规律代码结构清晰适合作为仿真学习与二次开发基础。1. 光斑仿真撞上大气湍流先做一个会“呼吸”的相位屏激光从发射端出发在大气里走过几百米甚至几公里后接收端的探测器上很难再看到教科书里那种边缘光滑的高斯光斑。它可能整体漂移可能边缘抖动可能亮度快速闪烁湍流重的时候甚至会裂成几个不连续的小斑。做自由空间光通信、激光雷达、光学遥感和地面天文观测的人都会被同一个问题卡住这套“高斯光束大气传输”过程到底能不能在电脑里一帧一帧复现出来而且复现得足够可信。能但路径和大多数人第一反应完全不同。解析公式给得出远场发散角的平均值却给不了一帧“现场照片”要同时还原光斑的形变、质心抖动和闪烁实际工程里最通用的办法是放弃一次性积分改用“传输叠加”把光束沿路径一步一步往前推每一步经过一个随机相位屏用这个屏等效代替这一段大气湍流对波前的作用。这个屏就是标题里的“大气湍流屏”整个光斑仿真工作基本就是围绕它的生成、校准、插入和统计展开的。下文从零开始把这条链路整个走通先讲湍流屏的数值生成再讲高斯光束的分步传播模型再把两者拼成完整的大气传输仿真系统最后说怎么判断结果可信以及批量仿真时值得关注的提速手段。2. 大气湍流屏生成从Kolmogorov谱到FFT功率谱反演2.1 Cn²、内尺度和外尺度相位屏背后是哪几个物理量大气湍流的本质是温度、湿度微扰引起的折射率随机起伏。工程上不需要追踪每一个涡旋的瞬时位置而是用统计量描述它。折射率结构常数 Cn² 是湍流强度的核心指标弱湍流大约在 1e-17~1e-16 m^(-2/3)中等强度约 1e-15强湍流能到 1e-14 甚至更高。另两个关键参数是尺度范围外尺度 L0 对应能量注入的大涡旋典型值 10~100 m内尺度 l0 对应粘性耗散的小涡旋通常几毫米到一厘米。L0 和 l0 直接决定了湍流功率谱的低频和高频形态。理想化的 Kolmogorov 谱写成 Φn(κ) 0.033 Cn² κ^(-11/3)只在中频段成立实际仿真时要在低频段加外尺度截断、在高频段加内尺度平滑这就是 von Karman 谱存在的意义。相位屏的生成本质上就是按这套功率谱生产一个二维随机场让它的空间统计特性和真实大气尽量一致。2.2 从功率谱到随机相位屏FFT谱反演的原理与限制谱反演法基于一个非常直接的思路把湍流看成大量随机傅里叶分量的叠加每个分量的幅度由功率谱密度决定相位是随机的。具体到相位屏先计算相位功率谱密度 Φφ(κ)它与折射率谱的关系在薄屏近似下是Φφ(κ) 2π k² Δz Φn(κ)其中 k 2π/λΔz 是这个屏等效代替的传输距离。拿到 Φφ 后用频域复高斯随机数乘上它的平方根再做一次二维傅里叶逆变换取实部就得到一帧弧度单位的相位屏。这个流程的局限也很明显FFT 的低频分辨率受屏尺寸限制不足 1/L0 的空间频率信息几乎全丢。结果就是直接生成的屏只有高频毛刺缺少大尺度波前倾斜和低阶像差光斑仿真时“漂移”效果出不来。这个问题放在 2.4 细说先看能直接跑的代码。2.3 可直接落地的Python相位屏生成函数下面这段代码按 von Karman 谱生成一帧相位屏参数全部显式传入方便逐个调试。import numpy as np def turbulent_phase_screen(N, dx, wavelength, Cn2, L0, l0, dz, seedNone): 通过 FFT 谱反演生成 von Karman 谱相位屏返回弧度。 rng np.random.default_rng(seed) fx np.fft.fftfreq(N, ddx) fxx, fyy np.meshgrid(fx, fx) kappa_sq fxx**2 fyy**2 kappa_sq np.maximum(kappa_sq, 1e-12) # 避免零频除零 k0 1.0 / L0 # 外尺度对应频率 km 5.92 / l0 # 内尺度对应频率 # von Karman 折射率功率谱 phin 0.033 * Cn2 * np.exp(-kappa_sq / km**2) / (kappa_sq k0**2)**(11.0 / 6.0) k 2 * np.pi / wavelength phi_phase 2 * np.pi * k**2 * dz * phin # 相位功率谱薄屏近似 gamma rng.standard_normal((N, N)) 1j * rng.standard_normal((N, N)) phase np.fft.ifft2(gamma * np.sqrt(phi_phase)).real return phase - phase.mean()N是网格边长建议取 2 的幂dx是网格间距wavelength用米做单位Cn2是湍流强度L0、l0是外尺度和内尺度dz是当前屏等效代表的大气厚度seed用于产生可复现的相位屏。返回的phase单位是弧度幅值通常在几个弧度到几十弧度之间不要期望它被限制在 ±π。提示不同教材里 km 的系数有 5.92 和 5.91 两种写法差异不到 0.2%对仿真结果几乎没有影响。外尺度 L0 如果设得太大低频分量增多需要更大的屏才能采样完整网格内存会明显上涨。2.4 低频缺失先认症状再决定要不要补偿跑通上面的函数后建议先做一件事把相位屏显示出来看它是不是只有细碎的颗粒感。如果是说明低频能量不足这是 FFT 谱反演的天然缺陷不一定是参数错了。要确认这一点可以粗略检查相位结构函数它与理论值的偏离在间距较大时会非常明显。def phase_structure_1d(phase, dx): 简化的一维相位结构函数用于低频缺失趋势判断。 N phase.shape[0] r np.arange(1, N // 4) D [] for rr in r: diff phase[rr:, :] - phase[:-rr, :] D.append(np.mean(diff**2)) return r * dx, np.array(D)一段大气引起的相位结构函数理论上是间距 r 的 5/3 次方增长。如果实测曲线在小间距段斜率偏陡、大间距段明显下弯基本可以断定外尺度取小了或者网格不够大。常见修正手段有三种第一加大屏尺寸让 L0 真正被采样到第二对低频子网格做次谐波补偿把低频能量补回来第三在链路上额外叠加前几阶 Zernike 多项式模拟倾斜和离焦等低阶像差。提示不要一上来就同时上三种修正。先固定 L050m、屏尺寸大于 20 倍光束直径跑完一帧看光斑如果质心抖动仍然接近零再考虑次谐波或 Zernike 叠加。3. 高斯光束传播模型角谱法里的网格、采样和混叠约束3.1 从解析高斯公式到数值分步的必要性基模高斯光束的解析表达式很漂亮给定束腰 w0 和波长 λ任意 z 处的光斑半径 w(z)、等相位面曲率半径 R(z) 都有闭式解瑞利距离 zR π w0² / λ 决定近远场分界。真空里这套公式精确、快速、零成本。可一旦中间插入随机相位屏光束不再是理想高斯模式解析工具立即失效。数值分步的思路是把整段传输路径切成许多小段每段内部当作真空传播段末叠加上一个相位屏循环往复。真空传播部分最常用两种实现菲涅尔衍射积分和角谱法。我一般直接用角谱法它在近场和远场都适用不需要像单次菲涅尔变换那样频繁判断距离范围唯一要小心的是频率坐标和网格大小。3.2 角谱法的传递函数与一段可复用的传播代码角谱法的本质是平面波展开先把光场做傅里叶变换到空间频率域乘上自由空间传递函数再逆变换回空间域。传递函数形式为H(fx, fy) exp(i k Δz √(1 - λ² fx² - λ² fy²))这里的 fx、fy 是空间频率单位是 cycles/m和 numpy.fft.fftfreq 的输出正好对应。下面是可直接复制的传播函数。def angular_spectrum(U, wavelength, dx, dz): 角谱法真空传播U 为复振幅场dz 为传播距离米。 N U.shape[0] k 2 * np.pi / wavelength fx np.fft.fftfreq(N, ddx) fxx, fyy np.meshgrid(fx, fx) term 1 - (wavelength * fxx)**2 - (wavelength * fyy)**2 term[term 0] 0 H np.exp(1j * k * dz * np.sqrt(term)) Uf np.fft.fft2(U) return np.fft.ifft2(Uf * H)高频超出物理范围的部分被直接截断实际工程中只要网格选得合理term基本不会出现负值。dx和N决定了空间频率的采样范围如果光斑扩散到超过网格边框边界处会折回混叠输出的强度图会出现对称伪影这是最常见的问题。3.3 网格、束腰与混叠一组经验参数网格选型不是拍脑袋。束腰 w0 决定了初始场的宽度衍射扩展决定了传播后的宽度两者加在一起必须小于网格边长的一半否则边界折返。场景波长网格 N网格间距 dx束腰 w0传播距离 L用途桌面验证633 nm2560.1 mm1 mm20 m快速跑通近地面链路1550 nm5120.5 mm2 cm1 km常规单模长距离激光通信1550 nm10241 mm5 cm5 km需要GPU经验准则是束腰至少占 20~40 个网格传播后的光斑半径不超过网格边长的四分之一。比如 633 nm、束腰 1 mm 的桌面场景衍射扩展后光斑约 4 mmN256、dx0.1 mm 对应网格宽度 25.6 mm余量充足。网格太小最先出现的现象不是误差变大而是光斑从一侧折回另一侧看起来像多了一个对称的“鬼像”。提示把angular_spectrum单独测一遍。让它传播一个理想高斯光束到真空中远场再把数值光斑半径与解析公式 w(z) 对比误差小于 0.5% 再往下走。这个步骤只要几分钟能省掉后面所有怀疑。4. 把湍流屏装进传播链路大气传输仿真系统与参数调优4.1 多屏布局几个屏、间距多少、放在哪单相位屏只能模拟很短的路径长距离必须用一串屏。把总路径 L 切成 M 段每段距离 Δz L/M生成 M 个独立的相位屏按顺序插入真空传播段之间就构成立体的大气传输模型。屏的数量不是越多越好。屏太稀每段路径太长等效薄屏近似失效屏太密相邻屏之间的随机场会引入不必要的计算量而且短距离内多次施加随机相位会让相位方差增长与实际路径积分不一致。我一般取屏间距不超过外尺度 L0屏数量控制在总路径除以 L0 到两倍这个范围。例如 L030 m1 km 路径取 20~40 个屏比较稳妥。4.2 完整的高斯光束大气传输代码下面这段把前面所有函数串成一条可执行链路输入参数输出接收面的复振幅场。def simulate_atmospheric_link(N, dx, wavelength, w0, Cn2, L0, l0, L, num_screens, seed0): 高斯光束经过大气湍流的端到端仿真返回接收面复振幅。 U gaussian_beam(N, dx, wavelength, w0) dz L / num_screens for i in range(num_screens): screen turbulent_phase_screen( N, dx, wavelength, Cn2, L0, l0, dz, seed1000 i ) U U * np.exp(1j * screen) U angular_spectrum(U, wavelength, dx, dz) return Ugaussian_beam生成束腰位置的基模高斯场循环体内每次先乘相位屏再真空传播。seed1000 i保证每次生成的相位屏不同且可复现。这个顺序不能反叠加相位屏必须发生在该段真空传播之前因为相位屏等效的是这一段大气累积的波前畸变。初始光场的归一化也要注意。比较稳妥的做法是把总功率归一化成 1 W这样接收面的|U|²直接就是强度分布后续算闪烁指数、Strehl 比都不用再标定功率。4.3 接收面光斑分析与湍流强度的对应关系跑完一次仿真接收面光斑无法用肉眼可靠判断结果通常要算三个量质心位置、RMS 束宽和闪烁指数。它们的定义分别是光强加权重心、光强二阶矩半径以及σ_I² I²/I² - 1。代码实现如下。def beam_metrics(I, dx): 计算光斑质心、RMS半径I 是二维强度分布。 N I.shape[0] x (np.arange(N) - N // 2) * dx xx, yy np.meshgrid(x, x) total I.sum() cx (xx * I).sum() / total cy (yy * I).sum() / total r2 (xx - cx)**2 (yy - cy)**2 rms np.sqrt((r2 * I).sum() / total) return cx, cy, rms不同湍流强度下光斑表现差异很大下表给的是典型定性特征单位一致时可直接用Cn² 量级1550 nm, 1 km光斑表现闪烁指数量级1e-16接近理想高斯质心轻微抖动0.01 左右1e-15边缘明显抖动光斑半径扩大0.1~0.21e-14出现破碎和暗区中心强度剧烈变化0.5 附近甚至更高如果看到的是“光斑形状几乎不变但整体位置来回摆动”多半是高频毛刺型相位屏造成的假象低频成分没有被正确补偿也就是第 2 章说的那种情况。如果看到光斑大面积扩展但质心稳定则要怀疑网格混叠或传播方向符号问题。提示批量仿真时不要每次重新生成同一个位置的相位屏。把多帧结果存成(frames, N, N)的三维数组后续统计只要一次np.mean(axis0)就够不要用 Python 列表逐帧求和。5. 验证仿真可信度与批量加速的实用技巧5.1 用统计量对齐理论值而不是盯着单帧光斑单帧光斑什么都证明不了。湍流是随机过程必须有足够的独立帧数做统计平均结果才可能与理论匹配。弱湍流条件下最容易对比的量是闪烁指数平面波 Rytov 方差给出σ_I² ≈ 1.23 Cn² k^(7/6) L^(11/6)虽然高斯光束的修正系数略有不同但量级和随 Cn² 的变化趋势可以作为重要参考。批量收集指标的做法是循环不同随机种子分别调用simulate_atmospheric_link对每组强度求闪烁指数和质心抖动标准差。帧数建议不少于 50 帧否则闪烁指数的统计涨落会比你想验证的误差还大。如果统计均值和理论值差出两倍以上优先检查 Cn² 的量纲和 k 的计算是否用了角波数 2π/λ 而不是 λ 本身。5.2 大样本仿真的提速手段批量 FFT 与显存管理一帧 1024² 的复数场约 16 MB链路里同时存在的场、相位屏和临时数组很容易超过 500 MB。跑 100 帧统计时最直接的提速办法是把循环改成批量模式一次性生成(frames, N, N)的强度数组用numpy.fft.fft2的 axes 参数在全部帧上并行处理在 GPU 上则用cupy.fft.fft2一次提交整个 batch。这样做不光省 Python 循环开销更重要的是把 FFT 调用变成少数几次巨大的矩阵运算显存带宽利用率高得多。这里要提一个类似 fast gauss sums via flash attention 的思路flash attention 之所以快核心不是少做浮点运算而是把数据分块组织减少对 HBM 的重复读写光斑仿真批处理面对的问题几乎一样——FFT 算得再快数据搬来搬去也快不起来。所以实际仿真时我会把 50 帧分成几个 block每次处理 8 帧保证显存占满但不溢出如果单帧已经用到 2048²建议把场类型从 complex128 降到 complex64内存直接砍半精度损失对光斑统计影响极小。最后留一个实用习惯先用 256² 网格把整条链路跑通像素级的统计指标确认正常后再一次性放大到 512² 或 1024²。放大后固定随机种子对比放大前后的质心抖动若差异超过 10%说明原来的网格就有混叠而不是湍流屏参数的问题。本文还有配套的精品资源点击获取 
长河君 · 2026/9/14 12:33:12 · 阅读 1.2k 大气湍流相位屏生成与高斯光束传播仿真实战
简介这是一份面向光学与激光通信研究的高斯光束大气湍流传播仿真MATLAB源码聚焦光斑仿真、大气湍流屏建模与大气传输效应分析适合科研人员、工程师及光学相关专业学生使用用于理解湍流对光束质量的影响并辅助后续补偿设计。压缩包仅含1个m源文件整包仅1KB代码精简无需额外依赖即可在MATLAB中直接运行便于快速复现和参数修改资源紧密围绕大气湍流中光束传播的核心问题利用随机相位屏模拟大气折射率起伏可设定不同湍流强度与传输距离计算接收面光强分布直观展示光斑扩展、扭曲与闪烁等现象。目前已有502人学习下载。通过运行gauss.m可评估湍流引起的光束质量退化量化相位扰动导致的波前畸变为激光通信、天文观测等提供理论依据与数据支持相关专业学生可通过修改代码参数深入探究不同湍流条件下的光斑演化规律代码结构清晰适合作为仿真学习与二次开发基础。1. 光斑仿真撞上大气湍流先做一个会“呼吸”的相位屏激光从发射端出发在大气里走过几百米甚至几公里后接收端的探测器上很难再看到教科书里那种边缘光滑的高斯光斑。它可能整体漂移可能边缘抖动可能亮度快速闪烁湍流重的时候甚至会裂成几个不连续的小斑。做自由空间光通信、激光雷达、光学遥感和地面天文观测的人都会被同一个问题卡住这套“高斯光束大气传输”过程到底能不能在电脑里一帧一帧复现出来而且复现得足够可信。能但路径和大多数人第一反应完全不同。解析公式给得出远场发散角的平均值却给不了一帧“现场照片”要同时还原光斑的形变、质心抖动和闪烁实际工程里最通用的办法是放弃一次性积分改用“传输叠加”把光束沿路径一步一步往前推每一步经过一个随机相位屏用这个屏等效代替这一段大气湍流对波前的作用。这个屏就是标题里的“大气湍流屏”整个光斑仿真工作基本就是围绕它的生成、校准、插入和统计展开的。下文从零开始把这条链路整个走通先讲湍流屏的数值生成再讲高斯光束的分步传播模型再把两者拼成完整的大气传输仿真系统最后说怎么判断结果可信以及批量仿真时值得关注的提速手段。2. 大气湍流屏生成从Kolmogorov谱到FFT功率谱反演2.1 Cn²、内尺度和外尺度相位屏背后是哪几个物理量大气湍流的本质是温度、湿度微扰引起的折射率随机起伏。工程上不需要追踪每一个涡旋的瞬时位置而是用统计量描述它。折射率结构常数 Cn² 是湍流强度的核心指标弱湍流大约在 1e-17~1e-16 m^(-2/3)中等强度约 1e-15强湍流能到 1e-14 甚至更高。另两个关键参数是尺度范围外尺度 L0 对应能量注入的大涡旋典型值 10~100 m内尺度 l0 对应粘性耗散的小涡旋通常几毫米到一厘米。L0 和 l0 直接决定了湍流功率谱的低频和高频形态。理想化的 Kolmogorov 谱写成 Φn(κ) 0.033 Cn² κ^(-11/3)只在中频段成立实际仿真时要在低频段加外尺度截断、在高频段加内尺度平滑这就是 von Karman 谱存在的意义。相位屏的生成本质上就是按这套功率谱生产一个二维随机场让它的空间统计特性和真实大气尽量一致。2.2 从功率谱到随机相位屏FFT谱反演的原理与限制谱反演法基于一个非常直接的思路把湍流看成大量随机傅里叶分量的叠加每个分量的幅度由功率谱密度决定相位是随机的。具体到相位屏先计算相位功率谱密度 Φφ(κ)它与折射率谱的关系在薄屏近似下是Φφ(κ) 2π k² Δz Φn(κ)其中 k 2π/λΔz 是这个屏等效代替的传输距离。拿到 Φφ 后用频域复高斯随机数乘上它的平方根再做一次二维傅里叶逆变换取实部就得到一帧弧度单位的相位屏。这个流程的局限也很明显FFT 的低频分辨率受屏尺寸限制不足 1/L0 的空间频率信息几乎全丢。结果就是直接生成的屏只有高频毛刺缺少大尺度波前倾斜和低阶像差光斑仿真时“漂移”效果出不来。这个问题放在 2.4 细说先看能直接跑的代码。2.3 可直接落地的Python相位屏生成函数下面这段代码按 von Karman 谱生成一帧相位屏参数全部显式传入方便逐个调试。import numpy as np def turbulent_phase_screen(N, dx, wavelength, Cn2, L0, l0, dz, seedNone): 通过 FFT 谱反演生成 von Karman 谱相位屏返回弧度。 rng np.random.default_rng(seed) fx np.fft.fftfreq(N, ddx) fxx, fyy np.meshgrid(fx, fx) kappa_sq fxx**2 fyy**2 kappa_sq np.maximum(kappa_sq, 1e-12) # 避免零频除零 k0 1.0 / L0 # 外尺度对应频率 km 5.92 / l0 # 内尺度对应频率 # von Karman 折射率功率谱 phin 0.033 * Cn2 * np.exp(-kappa_sq / km**2) / (kappa_sq k0**2)**(11.0 / 6.0) k 2 * np.pi / wavelength phi_phase 2 * np.pi * k**2 * dz * phin # 相位功率谱薄屏近似 gamma rng.standard_normal((N, N)) 1j * rng.standard_normal((N, N)) phase np.fft.ifft2(gamma * np.sqrt(phi_phase)).real return phase - phase.mean()N是网格边长建议取 2 的幂dx是网格间距wavelength用米做单位Cn2是湍流强度L0、l0是外尺度和内尺度dz是当前屏等效代表的大气厚度seed用于产生可复现的相位屏。返回的phase单位是弧度幅值通常在几个弧度到几十弧度之间不要期望它被限制在 ±π。提示不同教材里 km 的系数有 5.92 和 5.91 两种写法差异不到 0.2%对仿真结果几乎没有影响。外尺度 L0 如果设得太大低频分量增多需要更大的屏才能采样完整网格内存会明显上涨。2.4 低频缺失先认症状再决定要不要补偿跑通上面的函数后建议先做一件事把相位屏显示出来看它是不是只有细碎的颗粒感。如果是说明低频能量不足这是 FFT 谱反演的天然缺陷不一定是参数错了。要确认这一点可以粗略检查相位结构函数它与理论值的偏离在间距较大时会非常明显。def phase_structure_1d(phase, dx): 简化的一维相位结构函数用于低频缺失趋势判断。 N phase.shape[0] r np.arange(1, N // 4) D [] for rr in r: diff phase[rr:, :] - phase[:-rr, :] D.append(np.mean(diff**2)) return r * dx, np.array(D)一段大气引起的相位结构函数理论上是间距 r 的 5/3 次方增长。如果实测曲线在小间距段斜率偏陡、大间距段明显下弯基本可以断定外尺度取小了或者网格不够大。常见修正手段有三种第一加大屏尺寸让 L0 真正被采样到第二对低频子网格做次谐波补偿把低频能量补回来第三在链路上额外叠加前几阶 Zernike 多项式模拟倾斜和离焦等低阶像差。提示不要一上来就同时上三种修正。先固定 L050m、屏尺寸大于 20 倍光束直径跑完一帧看光斑如果质心抖动仍然接近零再考虑次谐波或 Zernike 叠加。3. 高斯光束传播模型角谱法里的网格、采样和混叠约束3.1 从解析高斯公式到数值分步的必要性基模高斯光束的解析表达式很漂亮给定束腰 w0 和波长 λ任意 z 处的光斑半径 w(z)、等相位面曲率半径 R(z) 都有闭式解瑞利距离 zR π w0² / λ 决定近远场分界。真空里这套公式精确、快速、零成本。可一旦中间插入随机相位屏光束不再是理想高斯模式解析工具立即失效。数值分步的思路是把整段传输路径切成许多小段每段内部当作真空传播段末叠加上一个相位屏循环往复。真空传播部分最常用两种实现菲涅尔衍射积分和角谱法。我一般直接用角谱法它在近场和远场都适用不需要像单次菲涅尔变换那样频繁判断距离范围唯一要小心的是频率坐标和网格大小。3.2 角谱法的传递函数与一段可复用的传播代码角谱法的本质是平面波展开先把光场做傅里叶变换到空间频率域乘上自由空间传递函数再逆变换回空间域。传递函数形式为H(fx, fy) exp(i k Δz √(1 - λ² fx² - λ² fy²))这里的 fx、fy 是空间频率单位是 cycles/m和 numpy.fft.fftfreq 的输出正好对应。下面是可直接复制的传播函数。def angular_spectrum(U, wavelength, dx, dz): 角谱法真空传播U 为复振幅场dz 为传播距离米。 N U.shape[0] k 2 * np.pi / wavelength fx np.fft.fftfreq(N, ddx) fxx, fyy np.meshgrid(fx, fx) term 1 - (wavelength * fxx)**2 - (wavelength * fyy)**2 term[term 0] 0 H np.exp(1j * k * dz * np.sqrt(term)) Uf np.fft.fft2(U) return np.fft.ifft2(Uf * H)高频超出物理范围的部分被直接截断实际工程中只要网格选得合理term基本不会出现负值。dx和N决定了空间频率的采样范围如果光斑扩散到超过网格边框边界处会折回混叠输出的强度图会出现对称伪影这是最常见的问题。3.3 网格、束腰与混叠一组经验参数网格选型不是拍脑袋。束腰 w0 决定了初始场的宽度衍射扩展决定了传播后的宽度两者加在一起必须小于网格边长的一半否则边界折返。场景波长网格 N网格间距 dx束腰 w0传播距离 L用途桌面验证633 nm2560.1 mm1 mm20 m快速跑通近地面链路1550 nm5120.5 mm2 cm1 km常规单模长距离激光通信1550 nm10241 mm5 cm5 km需要GPU经验准则是束腰至少占 20~40 个网格传播后的光斑半径不超过网格边长的四分之一。比如 633 nm、束腰 1 mm 的桌面场景衍射扩展后光斑约 4 mmN256、dx0.1 mm 对应网格宽度 25.6 mm余量充足。网格太小最先出现的现象不是误差变大而是光斑从一侧折回另一侧看起来像多了一个对称的“鬼像”。提示把angular_spectrum单独测一遍。让它传播一个理想高斯光束到真空中远场再把数值光斑半径与解析公式 w(z) 对比误差小于 0.5% 再往下走。这个步骤只要几分钟能省掉后面所有怀疑。4. 把湍流屏装进传播链路大气传输仿真系统与参数调优4.1 多屏布局几个屏、间距多少、放在哪单相位屏只能模拟很短的路径长距离必须用一串屏。把总路径 L 切成 M 段每段距离 Δz L/M生成 M 个独立的相位屏按顺序插入真空传播段之间就构成立体的大气传输模型。屏的数量不是越多越好。屏太稀每段路径太长等效薄屏近似失效屏太密相邻屏之间的随机场会引入不必要的计算量而且短距离内多次施加随机相位会让相位方差增长与实际路径积分不一致。我一般取屏间距不超过外尺度 L0屏数量控制在总路径除以 L0 到两倍这个范围。例如 L030 m1 km 路径取 20~40 个屏比较稳妥。4.2 完整的高斯光束大气传输代码下面这段把前面所有函数串成一条可执行链路输入参数输出接收面的复振幅场。def simulate_atmospheric_link(N, dx, wavelength, w0, Cn2, L0, l0, L, num_screens, seed0): 高斯光束经过大气湍流的端到端仿真返回接收面复振幅。 U gaussian_beam(N, dx, wavelength, w0) dz L / num_screens for i in range(num_screens): screen turbulent_phase_screen( N, dx, wavelength, Cn2, L0, l0, dz, seed1000 i ) U U * np.exp(1j * screen) U angular_spectrum(U, wavelength, dx, dz) return Ugaussian_beam生成束腰位置的基模高斯场循环体内每次先乘相位屏再真空传播。seed1000 i保证每次生成的相位屏不同且可复现。这个顺序不能反叠加相位屏必须发生在该段真空传播之前因为相位屏等效的是这一段大气累积的波前畸变。初始光场的归一化也要注意。比较稳妥的做法是把总功率归一化成 1 W这样接收面的|U|²直接就是强度分布后续算闪烁指数、Strehl 比都不用再标定功率。4.3 接收面光斑分析与湍流强度的对应关系跑完一次仿真接收面光斑无法用肉眼可靠判断结果通常要算三个量质心位置、RMS 束宽和闪烁指数。它们的定义分别是光强加权重心、光强二阶矩半径以及σ_I² I²/I² - 1。代码实现如下。def beam_metrics(I, dx): 计算光斑质心、RMS半径I 是二维强度分布。 N I.shape[0] x (np.arange(N) - N // 2) * dx xx, yy np.meshgrid(x, x) total I.sum() cx (xx * I).sum() / total cy (yy * I).sum() / total r2 (xx - cx)**2 (yy - cy)**2 rms np.sqrt((r2 * I).sum() / total) return cx, cy, rms不同湍流强度下光斑表现差异很大下表给的是典型定性特征单位一致时可直接用Cn² 量级1550 nm, 1 km光斑表现闪烁指数量级1e-16接近理想高斯质心轻微抖动0.01 左右1e-15边缘明显抖动光斑半径扩大0.1~0.21e-14出现破碎和暗区中心强度剧烈变化0.5 附近甚至更高如果看到的是“光斑形状几乎不变但整体位置来回摆动”多半是高频毛刺型相位屏造成的假象低频成分没有被正确补偿也就是第 2 章说的那种情况。如果看到光斑大面积扩展但质心稳定则要怀疑网格混叠或传播方向符号问题。提示批量仿真时不要每次重新生成同一个位置的相位屏。把多帧结果存成(frames, N, N)的三维数组后续统计只要一次np.mean(axis0)就够不要用 Python 列表逐帧求和。5. 验证仿真可信度与批量加速的实用技巧5.1 用统计量对齐理论值而不是盯着单帧光斑单帧光斑什么都证明不了。湍流是随机过程必须有足够的独立帧数做统计平均结果才可能与理论匹配。弱湍流条件下最容易对比的量是闪烁指数平面波 Rytov 方差给出σ_I² ≈ 1.23 Cn² k^(7/6) L^(11/6)虽然高斯光束的修正系数略有不同但量级和随 Cn² 的变化趋势可以作为重要参考。批量收集指标的做法是循环不同随机种子分别调用simulate_atmospheric_link对每组强度求闪烁指数和质心抖动标准差。帧数建议不少于 50 帧否则闪烁指数的统计涨落会比你想验证的误差还大。如果统计均值和理论值差出两倍以上优先检查 Cn² 的量纲和 k 的计算是否用了角波数 2π/λ 而不是 λ 本身。5.2 大样本仿真的提速手段批量 FFT 与显存管理一帧 1024² 的复数场约 16 MB链路里同时存在的场、相位屏和临时数组很容易超过 500 MB。跑 100 帧统计时最直接的提速办法是把循环改成批量模式一次性生成(frames, N, N)的强度数组用numpy.fft.fft2的 axes 参数在全部帧上并行处理在 GPU 上则用cupy.fft.fft2一次提交整个 batch。这样做不光省 Python 循环开销更重要的是把 FFT 调用变成少数几次巨大的矩阵运算显存带宽利用率高得多。这里要提一个类似 fast gauss sums via flash attention 的思路flash attention 之所以快核心不是少做浮点运算而是把数据分块组织减少对 HBM 的重复读写光斑仿真批处理面对的问题几乎一样——FFT 算得再快数据搬来搬去也快不起来。所以实际仿真时我会把 50 帧分成几个 block每次处理 8 帧保证显存占满但不溢出如果单帧已经用到 2048²建议把场类型从 complex128 降到 complex64内存直接砍半精度损失对光斑统计影响极小。最后留一个实用习惯先用 256² 网格把整条链路跑通像素级的统计指标确认正常后再一次性放大到 512² 或 1024²。放大后固定随机种子对比放大前后的质心抖动若差异超过 10%说明原来的网格就有混叠而不是湍流屏参数的问题。本文还有配套的精品资源点击获取