简介一套基于修正Von-Karman大气湍流模型的高斯光束传输仿真系统面向激光大气传输、湍流效应与光束质量评估领域的研究人员、工程师及高年级学生。系统以三次次谐波补偿的多随机相位屏技术为核心能够较精确地模拟激光在湍流介质中传播时由温度、风速扰动引起的波前随机变化并展示光束强度分布、波前畸变、光束漂移等随传播距离和湍流强度演变的过程。压缩包共含6个文件主体为MATLAB源码包括ift2.m、turbulence.m等核心脚本另有txt、docx、md说明材料及LICENSE文件整体仅49KB结构简洁、便于快速运行与二次开发。该资源已有68人学习适用于激光通信、遥感探测、激光雷达及激光武器等系统的性能仿真还可通过修改波长、湍流强度、传输距离等参数为优化光学系统设计和制定湍流补偿策略提供数值依据。1. 大气湍流仿真为什么要选修正Von-Karman模型做激光大气传输实验的工程师大概都有过这样的经历实验室台架上高斯光束漂漂亮亮一到外场测试光斑抖动、能量扩散、接收端信号衰落得厉害。要解释这些现象不能靠拍脑袋得先把大气湍流对光波的作用量化出来。最通用的做法就是用随机相位屏模拟湍流介质的相位扰动而相位屏的统计特性完全由大气湍流模型决定。这篇要讲的这套仿真系统核心是修正Von-Karman大气湍流模型与高斯光束传输的数值仿真通过三次次谐波补偿的多随机相位屏技术模拟光波在湍流介质中的传播过程用于研究激光大气传输特性以及湍流效应对光束质量的影响。对正在做自由空间光通信链路预算、激光雷达或高能激光传输评估的工程师来说这套方法基本是绕不过去的基础工具。2. 修正Von-Karman湍流谱与相位屏生成原理为什么低频截断和三次次谐波缺一不可2.1 从Kolmogorov到修正Von-Karman低频截断决定仿真可信度大气折射率起伏的功率谱密度经典描述是Kolmogorov谱。它在惯性子区间成立得很好但有两个麻烦频率趋于零时谱密度发散导致相位屏的方差在数值上被虚假放大频率很高时也缺少截止网格采样不足就在高频出现混淆。实际大气中湍流的能量注入发生在尺度L0外尺度耗散发生在尺度l0内尺度以内所以修正Von-Karman模型在原有-11/3幂律两侧加上截断——低频用外尺度L0抑制发散高频用内尺度l0做指数衰减。数值仿真里截断不是可选项。傅里叶变换生成随机相位屏时低频分量决定了光斑的整体偏移和波前低阶畸变。如果低频发散生成出的相位屏在空间上会出现大范围缓慢起伏叠加高斯光束后等效于随机加了一个离焦或倾斜波前这不符合外场实测的统计结果。修正Von-Karman的相位功率谱写成Φ_φ(k) 0.023 r0^(-5/3) (k² k0²)^(-11/6) exp(-k²/km²)其中k02π/L0km5.92/l0r0是大气相干长度。这一形式直接进相位屏的傅里叶域生成流程不用先算C_n²再换算工程上最省事。选择修正Von-Karman模型而不是Kolmogorov模型最直接的原因就在这个低频截断项它保证相位结构函数在尺度接近L0时饱和和实测数据的吻合度更高。提示k0和km的单位必须与网格空间频率一致。常见错误是把L0当特征长度直接代入k空间导致低频截断不起作用仿真结果退化回Kolmogorov特性。2.2 三次次谐波补偿补回FFT网格丢失的低频尾巴直接用FFT在离散网格上生成相位屏最低可表示的空间频率受相位屏边长D限制Δf1/D。外尺度L0通常远大于D低频部分落在Δf以下被整体截掉了。这部分丢失能量恰好对应光束的倾斜、离焦等低阶像差。缺失的后果很直观仿真光束的质心漂移偏小长曝光光斑半径被低估闪烁指数偏低于实际。三次次谐波补偿是对低频段的分层再采样。常见做法是把频率零点附近的区域按3倍、9倍、27倍细分网格对每个子网格重新生成随机傅里叶系数再将这些系数通过插值叠加回主相位屏。这样等效地把最低可表示频率向下扩展了3^p倍p取到3时低频扩展了27倍。三次这个三指的就是p1、2、3这三层迭代每层把频率分辨率提高3倍。具体实现时每一层子网格的大小是3^p×3^p子网格内的频率步长为Δf/3^p。子网格内的每个频点仍按修正Von-Karman谱计算幅度相位取均匀随机分布。生成后对子网格做逆傅里叶变换再插值到主网格尺寸叠加到主相位屏上。叠加时注意幅度归一化否则低频分量被重复计数结构函数会在大尺度处鼓包。def generate_vk_phase_screen(D, N, r0, L0, l0, seedNone): import numpy as np rng np.random.default_rng(seed) dx D / N fx np.fft.fftfreq(N, ddx) FX, FY np.meshgrid(fx, fx) F np.sqrt(FX**2 FY**2) F[F 0] 1e-12 k0 2 * np.pi / L0 km 5.92 / l0 K 2 * np.pi * F # 修正Von-Karman相位功率谱r0直接控制湍流强度 Phi 0.023 * r0**(-5/3) * np.exp(-K**2 / km**2) / (K**2 k0**2)**(11/6) # 随机复高斯系数保证相位屏实部统计正确 noise rng.standard_normal((N, N)) 1j * rng.standard_normal((N, N)) S np.sqrt(Phi) * noise phi np.fft.ifft2(S).real * (N**2) / (D**2) * dx**2 phi phi - phi.mean() # 去掉均值项消除整体相位偏移 # 三次次谐波补偿低频分量p1,2,3对应3倍/9倍/27倍细分 for p in range(1, 4): Ns 3**p dF 1.0 / D fs (np.arange(Ns) - Ns / 2.0) * dF / Ns FSX, FSY np.meshgrid(fs, fs) FSm np.sqrt(FSX**2 FSY**2) FSm[FSm 0] 1e-12 Km 2 * np.pi * FSm Phi_sub 0.023 * r0**(-5/3) * np.exp(-Km**2 / km**2) / (Km**2 k0**2)**(11/6) noise_sub rng.standard_normal((Ns, Ns)) 1j * rng.standard_normal((Ns, Ns)) S_sub np.sqrt(Phi_sub) * noise_sub phi_sub np.fft.ifft2(S_sub).real * (Ns**2) / (D**2) * dF**2 # 上采样到主网格要求N能被Ns整除否则尺寸对不上 phi_sub np.kron(phi_sub, np.ones((N // Ns, N // Ns))) phi phi_sub return phi这段代码说明了两个关键技术点。第一谱采用的是修正Von-Karman的相位功率谱而非折射率谱因此可直接对相位屏做逆傅里叶变换。第二次谐波补偿循环里每一层子网格尺寸Ns按3倍递增频率步长随之细化最后用np.kron将子网格相位上采样到主网格尺寸叠加。这里的上采样因子是N//Ns要求N能被3的幂整除所以N建议取432、540这类同时含27因子的网格数否则上采样尺寸对不齐。2.3 多随机相位屏的排布距离分段决定仿真分辨率单张相位屏只能模拟一次集中扰动真实大气湍流沿传播路径连续分布所以需要用多张相位屏分段逼近。常见做法是把传输距离L分成M段每段长度ΔzL/M在每段末尾放一张相位屏。每张相位屏的r0按该段距离折算如果整条路径Cn²均匀则每段的r0_segr0_total*M^(-3/5)。如果路径不均匀更稳妥的做法是按每段单独算出r0再生成对应强度相位屏。段数M的选择直接影响仿真时间与精度。M过少时湍流被当作几个薄透镜集中作用光束的衍射效应在屏间被过度放大M过多时单段相位屏的统计方差很小却要付出M倍的FFT和传播计算代价。工程上我一般按Δz0.2*r0²/λ来估段数但这只是粗准则具体还要配合弱湍流条件单段的Rytov方差应小于0.1~0.3否则相位屏间交互作用被低估。3. 高斯光束传输仿真实现从网格设计到分步传播代码3.1 网格与波长参数怎么定先定边界再定采样高斯光束传输仿真第一件事不是写代码而是把网格尺寸D、网格数N、波长λ三者绑在一起定死。基本原则有三条D要至少是束腰半径的6~8倍保证衍射扩展后光斑不撞边界N要满足dxD/N小于衍射特征尺度否则高频相位细节被采样丢失λ决定传播相位步长段数Δz要小于衍射尺度dz_dπ*w0²/λ的若干分之一。举个例子波长λ1.064μm束腰w02cmD取20cmN取540。那么dx约0.37mm衍射尺度约为π*(0.02)²/1.064e-6≈1.18km。若总距离L4km分16段每段Δz250m大约是衍射尺度的五分之一满足分步传播的采样经验。网格数N选540而不是512是为了让次谐波上采样N//27能被整除这对3层次谐波格外重要。3.2 高斯光束初始化与相位屏加载核心代码与参数说明光束初始化时束腰处的电场写成高斯型。若束腰不在z0处还需要乘一个二次相位因子来给出初始波前曲率。下面这段初始化代码把束腰放在相位屏阵列的第一屏之前用平面波前近似def initial_beam(N, D, w0, wavelength, seed_phase0.0): dx D / N x (np.arange(N) - N / 2.0) * dx X, Y np.meshgrid(x, x) r2 X**2 Y**2 # 束腰位置的高斯场w0单位米与D保持同一坐标系 u np.exp(-r2 / w0**2) * np.exp(1j * seed_phase) return u这段代码中w0的单位是米必须与D、dx保持同一坐标系。seed_phase用来给整束光加一个全局相位对光束质量评估没有影响但如果要做相干合成或多光束仿真这个相位就不能随意置零。相位屏加载到光束上是逐段进行的。每段传播先用角谱法走真空距离再乘上该段的相位屏。角谱传播的传递函数写成H(f) exp( i2π/λ * sqrt(1 - (λfx)² - (λ*fy)²) * Δz )注意fx、fy是空间频率1/mlambda*f是一个无量纲比值这能保证近轴条件下退化为菲涅尔传播。def propagate_angular_spectrum(u, dz, wavelength, D, N): dx D / N fx np.fft.fftfreq(N, ddx) FX, FY np.meshgrid(fx, fx) k 2 * np.pi / wavelength # 角谱传递函数单位统一后lambda*f无量纲 exp_arg 1.0 - (wavelength * FX)**2 - (wavelength * FY)**2 exp_arg np.clip(exp_arg, 1e-12, None) H np.exp(1j * k * dz * np.sqrt(exp_arg)) U np.fft.fft2(u) u_out np.fft.ifft2(U * H) return u_out def multi_screen_propagate(u, phase_screens, dz, wavelength, D, N): for phi in phase_screens: u propagate_angular_spectrum(u, dz, wavelength, D, N) u u * np.exp(1j * phi) return u传递函数里exp_arg被clip到不小于1e-12是因为高频分量对应的λf可能超过1此时对应的倏逝波项在自由空间实际不存在。clip处理让指数项得到一个很大的负虚部等效衰减对于传输距离远小于数值容差的系统足够。若把仿真距离拉长到几十公里更稳妥的做法是彻底判断exp_arg0时置H为0避免数值噪声累积。multi_screen_propagate里先传播再乘相位屏顺序不能反过来否则湍流作用点就落在了每段的起点而非终点物理图景与实际不符。3.3 多相位屏传输循环步长与相位屏数量的联动调整多屏循环里最容易忽略的是最后一段没有相位屏。物理上接收端之前的那段距离同样有湍流所以循环结束后还应补一次角谱传播到接收面。另外每段相位屏的强度可以不同路径上如果有关键高度层比如近地面边界层可以把段距加密但每条段距内的r0需要重新换算。我的习惯是先在脚本里把所有段的Δz和r0_seg打出来核对再进入循环避免段距不一致导致相位屏方差集中。z_points np.linspace(0, L, M1) dz z_points[1] - z_points[0] r0_seg r0_total * M**(-3/5) # 均匀路径下的单屏折算 screens [generate_vk_phase_screen(D, N, r0_seg, L0, l0, seed100i) for i in range(M)] u initial_beam(N, D, w0, wavelength) # 先走完所有带相位屏的段 u multi_screen_propagate(u, screens, dz, wavelength, D, N) # 最后一段补到接收面如果L不是M的整数倍 u propagate_angular_spectrum(u, L - M*dz, wavelength, D, N)这段代码把均匀路径的折算写得很直观总r0按M^(-3/5)分到每一段。seed100i保证每个相位屏不同但可复现是后续做参数对比的基础。注意r0_seg是分段距离Δz对应的相干长度不是整条路径的r0两者差着M的幂次搞混是新手最容易犯的错误。如果L不是M的整数倍最后一小段需要单独补传否则光束相当于少走了一段距离Strehl比会偏高。4. 仿真结果分析湍流对光束质量的影响如何量化评估4.1 光束质量评价指标从强度分布提炼三个核心数字拿到接收面光场u后第一件事是算强度I|u|²。光束质量评估最常用的三个指标是Strehl比、二阶矩光束宽度和闪烁指数。Strehl比定义为实际峰值强度与无湍流时峰值强度之比弱湍流下它直接反映像差对聚焦性能的劣化。二阶矩光束宽度定义在强度加权上对光斑形状不敏感适合描述能量扩散趋势。闪烁指数σ_I²刻画光强起伏计算公式为⟨I²⟩/⟨I⟩² − 1。这三个指标都要对多次独立实现取平均单次仿真的强度分布只能看形态不能下统计结论。下面这段计算代码把三个指标一起算出来def beam_quality_metrics(u, u_ideal, N, D): I np.abs(u)**2 I_ideal np.abs(u_ideal)**2 P I.sum() dx D / N x (np.arange(N) - N/2.0) * dx X, Y np.meshgrid(x, x) r2 X**2 Y**2 # 二阶矩半径强度加权均方根半径反映能量扩散程度 rms_r np.sqrt((r2 * I).sum() / P) # Strehl比峰值强度之比需保证两束总能量一致 strehl I.max() / I_ideal.max() # 闪烁指数无量纲光强起伏用于和Rytov理论对比 mean_I I.mean() scint (I**2).mean() / mean_I**2 - 1.0 return {rms_r: rms_r, strehl: strehl, scint: scint}这段代码的Strehl比用峰值比避免了绝对强度归一化的麻烦但要保证u_ideal与u的总能量相同否则比值失真。二阶矩半径用的坐标x从中心负半轴到正半轴数组默认原点在索引N//2所以生成坐标时必须先减N/2再乘dx否则所有半径都算偏一个格点。闪烁指数直接对整张强度图计算对接收孔径内的平均有意义若接收口径小于光斑则要先做孔径截取再计算不能拿全图平均去和实验对比。4.2 不同湍流强度下的仿真对比r0、L0、l0三个参数的影响顺序三个谱参数里r0的变化对结果影响最直接它决定相位屏的方差幅度。r0小则湍流强Strehl比迅速下降闪烁指数上升光斑质心波动显著。L0影响的是低频倾斜和离焦成分的含量L0越大低频能量越多光束质心的漂移越大但L0对短曝光强度分布的干扰弱于对长曝光累积分布的干扰。l0的影响主要集中在高频端对光束质量指标影响较小但对闪烁指数的相位调制项有贡献。做参数扫描时建议固定其他参数单独扫描r0和L0。常用扫描式是保持同一组随机种子依次改变r0这样每次变化只来自参数本身排除随机实现差异。扫L0时L0取10m、50m、100m观察质心漂移方差和Strehl比的变化趋势可以发现L0大于一个阈值后影响趋缓这对应大气外尺度饱和效应——这个拐点直观展示了为什么要用修正Von-Karman而不是无截断的Kolmogorov谱。实际处理中我把每组参数跑30次实现并保存均值与标准差画成误差棒图再对比理论曲线比单次运行更有说服力。5. 避坑指南随机相位屏仿真的7个常见问题排查5.1 低频次谐波补偿量级不对光斑质心漂移偏小现象仿真输出的光斑抖动范围明显小于外场实测长曝光光斑中心基本不动。原因主FFT网格能表示的最低空间频率是1/D次谐波补偿层数不够或子网格尺寸错误。子网格的Ns用3^p后最低频率扩展为1/(3^p D)p只到2时扩展9倍低频倾斜分量仍然缺失。解决把p最大值设为3并确认主网格数N能被27整除否则子网格上采样后尺寸对不上补偿会失效。5.2 修正Von-Karman谱高频项写错导致闪烁指数异常现象仿真闪烁指数随距离增长的速度快于理论值且l0缩小后指标几乎不变。原因功率谱里高频截断形式用错写成exp(-k²/l0²)而不是exp(-k²/km²)其中km5.92/l0截断尺度差了一个数量级高频能量被错误保留。解决统一用km5.92/l0并确认代码里指数项是exp(-K²/km²)不要直接拿l0当波长代入。5.3 相位屏周期性伪影污染光斑边缘现象接收面强度图在四角出现对称条纹光斑外围有微弱网格状结构。原因FFT假设相位屏是周期的边界不匹配会在传播中产生虚假衍射。修正Von-Karman的低频截断缓解了一部分但D不够大时光束衍射后仍会碰到边界。解决将D加大到至少8倍束腰半径或在相位屏边界乘Hanning窗函数但要接受窗函数会压低边缘相位幅度对质心漂移统计有影响只能作为临时手段。5.4 随机种子管理混乱导致对比实验失真现象同一参数扫描下Strehl比曲线剧烈跳动看不出单调趋势。原因每次调用生成函数都用了新的随机种子不同参数对应的相位屏实现完全独立噪声淹没了参数效应。解决把所有仿真共用一个种子列表例如seed1000param_index确保扫描时每个参数点复用同一组湍流实现。同时把生成器和传播过程的随机状态记录到文件便于复盘。5.5 段数M与Δz不匹配导致能量守恒破坏现象传输总能量逐渐增大或衰减与理论不符。原因相位屏只乘复指数理论上不改变能量但角谱传播的H在exp_arg趋近于0时数值不稳定段数太多会把数值误差逐段放大。解决检查M段内是否存在exp_arg接近0的分量若存在把N减小或把网格频率上限控制在奈奎斯特极限的0.8倍以内并输出每两步之间的能量差做断言。# 在传播循环里加能量断言及时暴露数值不稳定 energy_before np.abs(u)**2 u propagate_angular_spectrum(u, dz, wavelength, D, N) u u * np.exp(1j * phi) energy_after np.abs(u)**2 assert abs(energy_after.sum() - energy_before.sum()) / energy_before.sum() 1e-6, \ 能量不守恒段距或网格有问题传播函数本应保持能量守恒乘相位屏也只是改变相位不改变幅度。能量偏差超过1e-6的绝对量级基本可以断定是角谱传递函数数值溢出。加上这个断言后参数扫描的每一组都能自动暴露问题不用等到结果图出来才发现曲线乱跳。5.6 网格数N与次谐波上采样尺寸冲突现象程序在np.kron处报错shape对不上。原因N取512这类2的幂时N//Ns不是整数。解决N选540、648这类同时含27因子的数。如果又需要2的幂网格做快速FFT可以把次谐波做在独立的小网格上再用双线性插值叠加到主屏而不是用kron直接铺。5.7 相位屏与传播顺序颠倒导致湍流作用方向反了现象光束先过强湍流再慢慢聚焦和预期的先传播后扰动顺序不符。原因循环里先乘相位屏再调传播函数导致湍流在段首作用。解决保持先角谱传播、再乘相位屏的顺序只有最后一段接收面不需要乘屏。这种细节不会报错只能靠逐屏检查相位屏生成顺序和传播距离来发现。6. 验证仿真结果的进阶方法统计一致性校验与参数敏感性分析多相位屏仿真跑通只是第一步结果能不能信要看统计性质是否收敛到理论预期。我最常用的一组校验是固定r0和L0重复运行30次独立实现计算Strehl比、闪烁指数和质心漂移的均值与方差然后把均值与弱湍流理论公式对照。闪烁指数的理论值是σ_I²≈1.23C_n²k^(7/6)*L^(11/6)可以把r0换算成C_n²后对拍。偏差超过20%通常指向谱形式或归一化错误而不是统计学波动这时回去检查代码远比加跑次数更有用。第二组校验是结构函数校验。生成相位屏后直接计算相位差结构函数D_φ(r)⟨(φ(xr)−φ(x))²⟩与修正Von-Karman公式的理论解画在同一条对数坐标轴上。前几个采样间隔的点如果偏差小于5%说明网格与谱形式正确大尺度r接近L0时应当出现饱和平台如果没有饱和说明低频截断被绕过。这个方法能直接检验次谐波补偿是否真的在起作用比看光斑形态可靠得多。实际验证中我会把r作为横轴、D_φ(r)作为纵轴同一个图上画出自定义判定曲线程序输出是否落在包络内的布尔量这样能自动跑完整组参数扫描。第三组是参数敏感性分析。把r0按0.5、1、2、4倍变化观察Strehl比的短曝光与长曝光曲线。如果曲线在窄区间内出现平台或跳变多半是某一步Δz超过了衍射尺度需要减少段长重跑。这个习惯帮我排掉了至少两版代码里的隐藏bug——不是算法错了而是段距和网格尺寸没对齐。希望帮到你先别急着换更复杂的湍流谱把修正Von-Karman和三次次谐波这一套吃透做绝大多数工程评估已经足够。本文还有配套的精品资源点击获取