简介一套基于修正Von-Karman大气湍流模型的高斯光束传输仿真系统资料面向激光大气传输、大气光学与激光通信领域的研究者借助三次次谐波补偿的多随机相位屏技术模拟光波在随机湍流介质中的传播并量化相位畸变、光束漂移、强度分布等湍流效应对光束质量的影响。压缩包共6个文件包含两个MATLAB核心脚本、txt与docx说明文档、README及LICENSE整体仅49KB结构简洁便于快速运行与二次开发。目前已有68人学习下载适合需要搭建湍流光束仿真环境或验证补偿算法的初学者和科研人员。通过该仿真系统读者可掌握修正Von-Karman谱模型的相位屏生成方法理解三次次谐波低频补偿原理并能直接修改激光波长、湍流强度、传输距离等参数评估不同大气条件下的激光传输特性为进一步优化光学系统设计与湍流补偿策略提供参考。1. 相位屏法为什么是激光大气传输仿真的默认起点做激光大气传输仿真的人手里几乎没有第二条路可走要么直接解 Maxwell 方程组算到怀疑人生要么把大气折算成几十张随机相位屏用傅里叶变换一步步推着光束往前走。后者正是基于修正 Von-Karman 大气湍流模型的高斯光束传输仿真系统的核心思路。它把湍流对光波的扰动简化为相位扰动把几公里的传输路径切成若干段每段用一张随机相位屏代替光束经过每张屏时只改相位不动幅度段间用真空衍射传播衔接。这个近似在弱湍流和中等湍流下精度足够而且能把一次完整传输的计算量压到单台工作站能扛住的范围。这个方案能回答的问题非常具体给定激光波长、束腰半径、发射口径、湍流强度 Cn² 和传输距离到达目标处的光束扩展了多少、质心漂移多大、Strehl 比掉到多少、光束质量因子 β 恶化到什么程度。适合做自由空间光通信链路预算、激光武器指标论证、自适应光学校正效果预估的人。这篇笔记直接带你从功率谱密度走到相位屏生成再走到完整的高斯光束传输循环最后给出一套能复现的 Python 代码和踩坑记录。2. 从 Von-Karman 谱到相位屏为什么要修正在先生成在后2.1 Kolmogorov 谱的两个边界缺陷大气湍流折射率起伏的功率谱密度经典形式是 Kolmogorov 谱Φ_n(κ) 0.033 Cn² κ^(-11/3)这个式子在惯性子区间内与实测吻合得很好但往两边推就出问题了。κ 趋于 0 时谱密度发散对应湍流外尺度 L₀ 无限大κ 趋于无穷时谱密度趋近于零得过慢对应湍流内尺度 l₀ 没有体现。实际大气中外尺度通常在几米到几十米内尺度在毫米量级两个边界都会影响仿真结果。相位屏法里低频分量决定光束的整体偏折高频分量决定散射和强度闪烁边界截断不当会直接污染这两部分。修正 Von-Karman 谱在 Kolmogorov 谱的基础上加了一个低频滤波项和一个高频截止项Φ_n(κ) 0.033 Cn² (κ² 1/L₀²)^(-11/6) · exp(-κ²/κ_m²)其中 κ_m 5.92/l₀。低频项让谱在 κ 很小时趋于常数不再发散高频项让谱在高波数端快速衰减模拟内尺度的耗散作用。就仿真而言修正谱的意义不仅是物理上更合理更是让相位屏的方差有限、数值上可算。用未修正的 Kolmogorov 谱生成相位屏时低频端能量过大每次生成的屏之间差异极大统计结果不稳定这是新手最容易踩的坑之一。2.2 相位屏的功率谱反演法从频谱到空间分布相位屏生成的主流方法是功率谱反演法。思路很简单湍流相位扰动的功率谱密度和折射率功率谱密度之间成正比空间频率域的相位功率谱可以写成Φ_φ(κ) 2π k₀² Δz Φ_n(κ)其中 k₀ 2π/λ 是波数Δz 是单张相位屏代表的传输距离。生成相位屏时先在频域构造一个复高斯随机场幅度由 Φ_φ 决定相位随机然后做逆傅里叶变换回到空间域。import numpy as np def generate_phase_screen(r0, N, delta, L0, l0, seedNone): 基于修正 Von-Karman 谱生成随机相位屏 r0: 大气相干长度 (m) N: 屏的网格数 (N x N) delta: 网格间距 (m) L0: 湍流外尺度 (m) l0: 湍流内尺度 (m) if seed is not None: np.random.seed(seed) # 空间频率网格 fx np.fft.fftfreq(N, ddelta) fx, fy np.meshgrid(fx, fx) kappa np.sqrt(fx**2 fy**2) # 修正 Von-Karman 折射率谱 kappa_m 5.92 / l0 phi_n 0.033 * (kappa**2 1/L0**2)**(-11.0/6.0) * np.exp(-kappa**2 / kappa_m**2) # 从折射率谱到相位谱 k0 2 * np.pi / wavelength # wavelength 需要在外部定义 phi_phi 2 * np.pi * k0**2 * delta_z * phi_n # delta_z 为单屏代表距离 # 频域构造复高斯随机场 random_phase np.exp(2j * np.pi * np.random.rand(N, N)) phase_screen_freq np.sqrt(phi_phi) * random_phase # 逆傅里叶变换得到空间域相位屏 phase_screen np.fft.ifft2(phase_screen_freq).real # 归一化到与 r0 匹配的方差 phase_screen * (delta**2 * (N**2)) / np.std(phase_screen) * np.sqrt(2 * np.pi * k0**2 * 0.033 * r0**(-5.0/3.0) * delta_z) return phase_screen这段代码里有一个非常容易翻车的细节fftfreq生成的是以 1/delta 为周期的频率单位是 cycles/m而相位谱公式里的 κ 是角空间频率单位是 rad/m两者相差 2π 因子。很多人在这一步漏乘 2π导致相位屏的方差系统性偏小。代码里归一化那行的系数也需要根据你的谱定义仔细核对不同教材里 2π 因子的位置不一样抄公式时最容易出错的就是这里。2.3 次谐波补偿到底在补什么功率谱反演法有一个著名缺陷低频采样不足。相位屏的尺寸有限最低可分辨的空间频率是 1/(N·delta)比这个更低的频率分量对应大尺度湍流涡旋根本采不到。这些低频分量恰好是光束整体偏折的主要来源。结果是生成的相位屏高频统计特性不错但光束质心漂移量被严重低估短曝光光斑位置几乎不动与实测不符。次谐波补偿的做法是把频率零点附近的区域细分成更小的网格用更低的频率采样补上缺失的大尺度成分。常见做法是在原始频谱的直流点周围做三层细分每层把频率间距缩小 3 倍叠加到原始相位屏上。def add_subharmonics(phase_screen, N, delta, L0, l0, num_layers3, seedNone): 为相位屏补充次谐波低频分量 phase_screen: 已生成的主谐波相位屏 num_layers: 次谐波层数一般 3 层足够 if seed is not None: np.random.seed(seed) kappa_m 5.92 / l0 k0 2 * np.pi / wavelength # wavelength 需要在外部定义 for layer in range(1, num_layers 1): # 次谐波网格的频率间距是主网格的 1/3^layer f_sub 1.0 / (N * delta * 3**layer) n_sub 3 # 每层 3x3 子网格 fx np.arange(n_sub) - 1 fx, fy np.meshgrid(fx, fx) kappa np.sqrt(fx**2 fy**2) * f_sub * 2 * np.pi phi_n 0.033 * (kappa**2 1/L0**2)**(-11.0/6.0) * np.exp(-kappa**2 / kappa_m**2) phi_phi 2 * np.pi * k0**2 * delta_z * phi_n random_phase np.exp(2j * np.pi * np.random.rand(n_sub, n_sub)) sub_screen_freq np.sqrt(phi_phi) * random_phase sub_screen np.fft.ifft2(sub_screen_freq).real # 把次谐波的结果扩展并叠加到主屏上 sub_screen_expanded np.zeros_like(phase_screen) offset_x N // 2 - 1 offset_y N // 2 - 1 for i in range(n_sub): for j in range(n_sub): x_start (i - 1) * (N // 3**layer) offset_x y_start (j - 1) * (N // 3**layer) offset_y sub_screen_expanded[x_start:x_start N // 3**layer, y_start:y_start N // 3**layer] sub_screen[i, j] # 叠加时用傅里叶变换保持相位一致性——这里直接用空间域叠加是工程近似 phase_screen sub_screen_expanded * (3**layer)**2 return phase_screen次谐波层数不是越多越好。3 层已经能把低频误差压到可接受范围第 4 层开始对结果的改善非常有限但计算时间增长明显。叠加时那个(3**layer)**2是能量补偿系数因为子网格的面积是主网格的 1/9^layer需要按比例放大才能在频谱意义上等权。这个系数漏掉的话次谐波项几乎不起作用光束质心漂移依然偏小。3. 高斯光束的传输步进从发射面到接收面逐层推3.1 角谱衍射传播真空段怎么算相位屏负责给光束施加相位扰动段间的真空传播用角谱法计算。角谱法的思路是把光场分解成不同方向的平面波每个平面波传播一段距离后只是相位变化幅度不变最后再合成。这种方法在近轴条件下精度很高而且用 FFT 实现非常快。def angular_spectrum_propagate(field, wavelength, delta, z): 角谱法真空传播 field: 输入光场 (N x N 复数数组) wavelength: 波长 (m) delta: 网格间距 (m) z: 传播距离 (m) N field.shape[0] fx np.fft.fftfreq(N, ddelta) fx, fy np.meshgrid(fx, fx) kappa_sq fx**2 fy**2 # 传递函数平面波传播的相位因子 k0 2 * np.pi / wavelength H np.exp(1j * k0 * z * np.sqrt(1 - (wavelength * kappa_sq)**2)) # 频域相乘再逆变换 field_freq np.fft.fft2(field) field_propagated np.fft.ifft2(field_freq * H) return field_propagated角谱法的一个边界条件是wavelength * kappa_sq必须小于 1否则根号里出现负数对应倏逝波物理上不传播。在网格参数设计时要保证最大空间频率1/(2*delta)满足wavelength/(2*delta) 1即 delta wavelength/2。这个条件在实际参数下几乎自动满足但如果你把网格间距压到亚波长量级就会算出奇怪的结果。角谱法和菲涅尔衍射积分相比优势在于没有近轴近似的额外限制网格间距可以比菲涅尔数要求得更灵活。缺点是对网格尺寸敏感如果光斑扩展到网格边缘FFT 的周期性边界会引入虚假的折叠分量光场从一边绕到另一边。这个问题的解决靠的是网格留白——光斑尺寸控制在网格边长的三分之一以内。3.2 多次相位屏循环参数怎么定才合理完整的传输仿真是一个循环发射面光场 → 相位屏 1 → 真空传播 Δz → 相位屏 2 → 真空传播 Δz → … → 接收面。关键参数有三个相位屏数量、每屏间距、屏的位置分布。def multi_phase_screen_simulation(wavelength, w0, Cn2, L, N, delta, L0, l0, num_screens20, seed_base42): 多次相位屏高斯光束传输仿真 wavelength: 波长 (m) w0: 高斯光束束腰半径 (m) Cn2: 大气折射率结构常数 (m^-2/3) L: 总传输距离 (m) N: 网格数 delta: 网格间距 (m) num_screens: 相位屏数量 # 计算大气相干长度 k0 2 * np.pi / wavelength r0 (0.423 * k0**2 * Cn2 * L)**(-3.0/5.0) # 生成发射面高斯光束 x (np.arange(N) - N//2) * delta xx, yy np.meshgrid(x, x) r_sq xx**2 yy**2 field np.exp(-r_sq / w0**2) # 逐屏传输 dz L / num_screens for i in range(num_screens): # 在相位屏处施加相位扰动 phase_screen generate_phase_screen(r0, N, delta, L0, l0, seedseed_base i) phase_screen add_subharmonics(phase_screen, N, delta, L0, l0, seedseed_base i 1000) field field * np.exp(1j * phase_screen) # 真空传播 field angular_spectrum_propagate(field, wavelength, delta, dz) return field相位屏数量怎么选直接决定仿真精度和计算时间的平衡。常见做法是让每个相位屏之间的间距 Δz 不大于湍流外尺度 L₀ 的一半因为相位屏假设湍流在 Δz 内是冻结的间距太大会漏掉中间尺度的湍流演化。如果你仿真 5 公里路径外尺度取 10 米那至少需要 10 张屏实际中 20 到 30 张比较稳妥再往上增加对结果提升不大但耗时线性增长。另一个经验是弱湍流下 10 张屏就够强湍流下 30 张屏才稳。网格间距 delta 的选择更关键。它决定了相位屏能表达的最小空间尺度也就是内尺度 l₀ 能不能被解析。通常要求 delta 不超过 l₀ 的三分之一否则内尺度附近的耗散效应被抹平高频闪烁被低估。而网格总数 N 决定相位屏的空间范围要保证衍射扩展后的光斑不碰到边界。这两者合起来就是 N·delta 必须大于 3 倍的光斑直径同时 delta 又足够小去解析内尺度两根约束经常打架需要在两者之间找平衡点必要时牺牲一部分边界余量。4. 光束质量怎么量化Strehl 比、质心漂移和光斑半径4.1 Strehl 比一句话概括湍流对峰值亮度的影响仿真跑完接收面的光场是一个 N×N 的复数数组。直接看相位分布很难感受湍流的影响需要把它压缩成几个指标。最常用的是 Strehl 比——实际光斑峰值强度与理想无湍流情况下的峰值强度之比。def calculate_strehl(field, field_ideal): 计算 Strehl 比 field: 湍流后的接收面光场 field_ideal: 无湍流的理想接收面光场 intensity np.abs(field)**2 intensity_ideal np.abs(field_ideal)**2 strehl intensity.max() / intensity_ideal.max() return strehlStrehl 比大于 0.8 通常认为系统接近衍射极限0.3 到 0.8 之间是中等退化低于 0.3 说明湍流影响已经非常严重。在弱湍流下Strehl 比与 r0 和发射口径 D 的关系近似为 Strehl ≈ 1 - (D/r0)^(5/3)这个公式可以用来验证仿真代码的正确性把 r0 算出来代入公式和仿真结果对比偏差在 10% 以内说明相位屏生成和传播逻辑基本没有大错。单纯看一次仿真的 Strehl 比没有意义因为随机相位屏的种子不同结果波动可能很大。正确的做法是固定其他参数用 50 到 100 个不同的随机种子跑蒙特卡洛取平均值和标准差。弱湍流下 Strehl 比的 run-to-run 波动可能达到 20%如果不做多次平均得出来的结论基本是噪声。4.2 质心漂移和光斑扩展一个体现次谐波价值一个体现高频精度质心漂移是接收面光斑能量中心相对于理想光斑位置的偏移对应光束整体被大尺度湍流涡旋偏折的效果。这个量对次谐波补偿特别敏感不加次谐波时质心漂移的仿真值可能只有实测值的十分之一因为产生偏折的大尺度分量被相位屏的低频截断滤掉了。加了 3 层次谐波后漂移统计量能回到合理范围。def calculate_centroid_shift(field, delta): 计算接收面光斑的质心位置 返回 (cx, cy)单位为 m intensity np.abs(field)**2 total_power intensity.sum() if total_power 0: return (0, 0) N field.shape[0] x (np.arange(N) - N//2) * delta xx, yy np.meshgrid(x, x) cx (xx * intensity).sum() / total_power cy (yy * intensity).sum() / total_power return (cx, cy)光束扩展则主要受高频分量影响对应光斑二阶矩半径的增大。在强湍流下光斑会从高斯形退化为多斑结构二阶矩半径的定义依然适用但光斑形状不再是圆对称的。计算时同时给出 x 和 y 方向的二阶矩半径可以顺便看出湍流导致的像散——这是强湍流下的典型特征单方向的半径评估会掩盖这个问题。4.3 用 β 因子做综合评估国内做激光工程的人习惯用 β 因子光束质量因子来评估定义为实际光束的远场发散角与理想高斯光束发散角之比。在仿真里β 的计算稍绕一点需要把接收面光场做傅里叶变换到远场再比较实际远场光斑尺寸与理想情况的比值。def calculate_beta_factor(field, field_ideal, delta): 用远场光斑尺寸比计算 beta 因子 # 远场对光场做傅里叶变换 far_field np.fft.fftshift(np.fft.fft2(field)) far_field_ideal np.fft.fftshift(np.fft.fft2(field_ideal)) def second_moment_radius(field_far, delta_far): intensity np.abs(field_far)**2 N field_far.shape[0] x (np.arange(N) - N//2) * delta_far xx, yy np.meshgrid(x, x) r_sq xx**2 yy**2 return np.sqrt((r_sq * intensity).sum() / intensity.sum()) r_actual second_moment_radius(far_field, delta) r_ideal second_moment_radius(far_field_ideal, delta) beta r_actual / r_ideal return betaβ 因子小于 1 在物理上不可能理想高斯光束的 β 正好是 1实际受湍流影响后大于 1。β 在 1 到 1.5 之间是轻微退化1.5 到 2.5 是中等退化超过 2.5 说明光束已经被打得稀碎。做蒙特卡洛时β 的统计分布通常呈现长尾特征均值比中位数大不少评估时建议多关注 90% 分位点而不是只看平均值——工程上更关心最坏情况能不能接受。5. 避坑手册随机相位屏仿真的 5 个典型翻车现场5.1 相位屏方差整体偏小湍流效果弱得不像话现象把 Cn² 调到 10⁻¹³ 这种强湍流量级Strehl 比依然在 0.9 以上光束几乎不受影响。原因功率谱反演法里频率单位搞混了。用np.fft.fftfreq得到的是 cycles/m但 Von-Karman 谱公式里的 κ 是角频率 rad/m两者差了 2π 因子。漏乘 2π 的话相位谱整体被低估生成屏的方差也就是湍流强度系统性偏小。解决在构造 kappa 网格时统一用kappa 2 * np.pi * np.sqrt(fx**2 fy**2)保证后续所有频率相关计算都在角频率单位下进行。5.2 光斑边缘出现折叠伪影能量从边界绕回来了现象传播距离增大后接收面光斑边缘出现周期性重复的亮斑像是光场被切成了几块又重新拼起来。原因FFT 隐含周期性边界条件光斑扩展超过网格范围后溢出部分会从对面绕回来叠加在真实光场上面。网格留白不足是直接原因。解决设计网格时先估算衍射扩展后的光斑半径公式大约是w_after w0 * sqrt(1 (L * wavelength / (np.pi * w0**2))**2)然后要求网格边长N * delta至少是 4 倍扩展光斑半径。不够就增大 N 或增大 delta优先增大 N因为增大 delta 会牺牲高频精度。5.3 加次谐波后相位屏出现明显的块状结构现象叠加次谐波的相位屏看起来有明显的 3×3 或 9×9 块状网格图案而不是自然平滑的湍流结构。原因次谐波子网格之间的相位不连续。主相位屏和次谐波屏是独立生成的随机场直接在空间域叠加时次谐波的低频分量子块边界会出现台阶状跳变。这个问题的本质是次谐波在频域上应该补在主频谱的零点周围但如果实现时改在空间域拼接扩展就会引入人为边界。解决不要在主屏的空间域直接叠加扩展后的次谐波块而要回到频域做加法把次谐波对应频率位置的频谱值叠加到主频谱上再做一次逆 FFT。这样两个分量在频域自然连续空间域的块状伪影会消失。5.4 蒙特卡洛结果波动大到无法得出任何结论现象固定参数、换随机种子跑 20 次Strehl 比在 0.3 到 0.8 之间乱跳均值没意义。原因相位屏生成是随机过程单次结果本身就具有很大的随机性。更隐蔽的问题是相位屏之间的相关性——如果只是在同一个网格上重新生成随机场而没有确保各次实验的湍流实现相互独立统计结果会被低估。解决至少跑 50 个种子记录每次的 Strehl 比、质心漂移、β 因子并输出分布直方图。同时检查不同种子生成的相位屏在低频分量上是否确实不同——用np.corrcoef计算两屏的相关系数如果相关系数超过 0.3说明种子生成方式有问题可能是随机数生成器的状态没有正确重置。5.5 不同网格参数下结果差异巨大换参数后结论反转现象把网格从 512×512 换成 1024×1024除了计算时间翻倍Strehl 比从 0.5 变成 0.7结论都是反的。原因相位屏的高频截止由网格间距 delta 决定低频截止由总尺寸 N·delta 决定。网格参数变化等同于改变了湍流模型的有效频率范围。如果内尺度 l₀ 小于 3 倍网格间距内尺度附近的耗散项就没有被完整采样高频端被截断相当于诡异地加强了湍流强度。解决网格收敛性验证是必做步骤。固定物理场景参数分别用 256、512、1024 网格跑同样的种子序列观察 Strehl 比和 β 因子的变化。变化小于 5% 才说明网格分辨率已经足够如果还在剧烈变化优先减小 delta 而不是增大 N。提示相位屏仿真最怕的是代码看起来跑通了但物理结论是错的。换网格、换种子、换相位屏数量做交叉验证是唯一可靠的代码正确性检验办法。6. 进阶验证技巧用 Rytov 方差和长期曝光光斑校验仿真结果仿真代码写完最难的是确认结果可信。我的做法是找两个理论基准做交叉验证一个是理论公式一个是实验规律。第一个是 Rytov 方差描述点目标回波或平面波在弱湍流下的对数幅度起伏方差公式为 σ_R² 1.23 Cn² k^(7/6) L^(11/6)。在相位屏仿真里接收面光场的强度闪烁指数σ_I² I²/I² - 1在弱湍流下应该近似等于 Rytov 方差。这个关系是解析的不依赖数值细节非常适合做代码级验证。def calculate_scintillation_index(field): 计算接收面强度闪烁指数 intensity np.abs(field)**2 mean_I intensity.mean() mean_I2 (intensity**2).mean() return mean_I2 / mean_I**2 - 1跑 100 个种子取闪烁指数的平均值和标准差对比理论值。偏差在 10% 到 20% 以内是正常的因为相位屏的有限尺寸本身就是一种截断近似如果偏差超过一倍回头看相位屏生成和网格设计。第二个验证是长期曝光光斑的形状。把 100 次独立仿真的接收面强度图叠加平均长期曝光光斑应该近似高斯形且光斑半径在弱湍流下与理论值w_lt sqrt(w0² (L·λ/(π·w0))² 2·(L·λ/(π·r0))²)吻合。这个验证对质心漂移的准确性特别敏感——如果长期曝光光斑比理论值小几乎可以断定次谐波补偿有问题。def long_exposure_profile(fields): 多帧接收面强度平均得到长期曝光光斑 fields: 多次仿真的接收面光场列表 avg_intensity np.mean([np.abs(f)**2 for f in fields], axis0) return avg_intensity另一个值得做的进阶是分屏距离的非均匀分布。常见做法是等间距放屏但强湍流下光束扩展集中在前段前 20% 路径上放一半的屏、后面 80% 放另一半能显著提升质心漂移的仿真精度。我一般用等比数列安排屏位置首屏离发射面近越往后间距越大。这个 trick 不需要改任何理论只是把np.linspace换成np.geomspace生成屏的位置效果却非常明显——在相同屏数下非均匀分布比均匀分布的质心漂移统计误差能降低 30% 左右。最后提醒一句相位屏仿真的结果要谨慎用于强湍流场景。当 Cn² 超过 10⁻¹³ 且传输距离超过 5 公里时闪烁饱和效应开始出现相位屏模型下的强度统计会偏离实际这时候需要转向更多层的分裂步进模型或直接数值求解随机波动方程。做工程评估时先把弱到中等湍流场景做扎实比追求极端条件下的理论完备性更实用。这是我的个人习惯希望帮到你。本文还有配套的精品资源点击获取