简介基于修正Von-Karman大气湍流模型的高斯光束传输仿真系统面向大气光学、激光通信与遥感探测领域的科研人员。系统利用三次次谐波补偿的多随机相位屏技术模拟光波在湍流介质中的传播可观察光强分布、波前畸变和光束漂移等参数随湍流强度的变化。资源包共6个文件核心为两个Matlab程序实现相位屏生成与光场传播计算另含说明文档、docx附赠资料、开源许可证和README整体仅49KB便于直接运行和二次开发目前已有68人学习下载。使用者可调整激光波长、湍流强度、传输距离等条件对比不同谐波补偿阶数的模拟效果从而评估激光系统在复杂大气环境下的大气传输性能。该仿真工具能为激光通信、遥感探测及自适应光学研究提供理论验证与参数优化支持。1. 高斯光束在湍流里跑偏多少取决于低频相位屏的还原程度同样的传输距离、同样的Cn2有人仿真出光斑抖动5微弧度有人仿真出11微弧度差距几乎全在相位屏的低频段。基于修正Von-Karman大气湍流模型的高斯光束传输仿真系统核心干的就是这件事用多随机相位屏技术把光波经过湍流介质的相位扰动拆成一帧一帧的屏再用三次次谐波补偿把离散网格丢失的低频功率补回来最终统计Strehl比、光斑半径和到达角起伏。这套资源适合正在做激光大气传输特性研究、需要复现对照数据或者想快速搭一套光束质量评估脚本的从业者。下面从模型原理开始拆一步步看到它是怎么把湍流“塞”进屏幕里的。2. 从Kolmogorov谱到修正Von-Karman先想清楚为什么是这个模型2.1 相位屏近似激光大气传输仿真绕不开的底衬大气湍流对光束的作用本质上是折射率随机起伏引起的波前相位畸变。光在路径上每走一小段相位就叠加一个随机扰动扰动攒到一定程度光强分布才跟着变化。严格求解这类问题要么数值解麦克斯韦方程组要么走广义Huygens-Fresnel积分路径计算量都不小。相位屏近似则把连续路径切成若干段每段湍流效应浓缩成一块垂直于光轴的屏屏幕上只修改相位、不修改振幅段与段之间的真空部分用衍射传播算子衔接。只要网格数够、屏数够弱起伏和中等起伏下的光斑特征都能复现所以它成了激光大气传输数值仿真的地基。多随机相位屏技术的常规做法是把总路径L分成M段每段长度dz在每段末端放一个相位屏起点不放或只放一个代表整段影响的屏。每个屏上的相位扰动按该段的湍流强度独立生成。这个环节最容易出问题的是“屏的质量”屏的低频分量不准后续所有统计量都会偏。做传输之前先花半小时把屏本身折腾明白比直接跑一百帧蒙特卡洛更有价值。2.2 Von-Karman和Kolmogorov的差别外尺度与高频截断Kolmogorov湍流谱是最经典的模型Φ(f) 0.033·Cn2·f^(-11/3)Cn2是大气折射率结构常数单位m^(-2/3)代表湍流强度f是空间频率。这个谱在低频端有天然缺陷f趋于0时谱值发散。实际大气里的湍涡尺度不可能无限大外尺度L0会截断最大涡旋的尺寸。数值仿真里网格只能表示有限频率范围低频缺失的部分恰恰决定光束整体漂移和光斑大尺度变形。Von-Karman谱在分母上加了一项Φ(f) 0.033·Cn2·(f² 1/L0²)^(-11/6)当f远大于1/L0时它退化成Kolmogorov的-11/3次幂当f趋于0时谱值收敛到有限值不再发散。修正Von-Karman谱再乘一个高频截断项exp(-f²/kl²)其中kl3.3/l0用来模拟内尺度l0对高频湍流的耗散。很多现成代码只做中间那步把Kolmogorov谱往网格里一扔就开始跑传播结果就是低频功率虚高、光斑偏移量偏大而且每次运行的随机性特别大。这个系统把模型选定为修正Von-Karman等于先从物理层面把低频和高频两端都约束住了后面相位屏才谈得上可信。2.3 三次次谐波补偿低频段的功率补回来相位屏是离散网格最低可表示的空间频率是Δf1/(N·dx)。在频谱矩阵里这个分量只对应中心那个网格。Von-Karman谱的能量集中在低频中心网格离散化后f在0到Δf之间的贡献几乎全部丢掉结果就是相位屏看起来太平、缺少大尺度起伏。次谐波补偿的思路很直接把主网格的中心低频单元用更细的子网格重新采样。常见做法是每阶把步长再除以3第一阶用3×3子网格第二阶9×9第三阶27×27。这样原来无法表达的最低频段按几何级数细分p取到3通常足够。取三次而不是更多是因为第四阶开始子网格面积变化对功率谱的增量已经小于插值误差继续加阶只增加计算时间结构函数几乎不再变。系统标题里特别强调的“三次次谐波补偿”落地就在这个循环里。3. 多随机相位屏生成从频谱构造到次谐波补偿的实现3.1 六个参数先定下来网格、外尺度和结构常数生成相位屏之前先定参数。以下是我反复使用的初始值参数符号建议初值作用与边界网格边长N256或512决定最高可表示频率与内存占用网格间距dx0.0020.01 m空间分辨率至少比内尺度小3倍折射率结构常数Cn21e-161e-14 m^(-2/3)湍流强度平方根关系进入相位屏方差湍流外尺度L01050 m决定低频功率拐点位置湍流内尺度l00.0010.01 m决定高频截断位置次谐波阶数num_sub3低频补偿深度N256对大多数光束质量统计够用N512用于看光斑细节或者算闪烁指数。dx取决于内尺度l0经验要求dx小于l0/3否则高频分量会折叠回低频光斑边缘出现伪纹理。Cn2按大气条件取白天近地面晴朗条件下1e-15这个量级比较典型。L0取大一些比如50 m低频拐点更接近真实边界层。这套系统主要面向水平路径和近地斜程路径如果算高空平台之间的链路外尺度可以取更大。提示如果只关心光斑统计而不看单帧细节N256、三级次谐波全开跑100帧的速度大约是512网格的4倍准确度差异不大。3.2 生成修正Von-Karman相位屏的主函数下面这段是相位屏生成的核心函数我在MATLAB里按这种方式组织输入输出都可以直接复用function ph vonkarman_screen(N, dx, Cn2, L0, l0, num_sub) % 生成修正Von-Karman谱随机相位屏 % N : 网格边长(像素) % dx : 网格间距(m) % Cn2 : 折射率结构常数(m^-2/3) % L0 : 外尺度(m) % l0 : 内尺度(m) % num_sub: 次谐波补偿阶数, 3即三次次谐波 fx (-N/2:N/2-1) / (N*dx); [FX, FY] meshgrid(fx); f sqrt(FX.^2 FY.^2); f(1,1) 1e-9; % 中心点置小值, 直流分量由次谐波共同表达 kl 3.3 / l0; % 高频截止波数 Phi 0.033 * Cn2 * (f.^2 1/L0^2).^(-11/6) .* exp(-f.^2/kl^2); df 1/(N*dx); H (randn(N) 1i*randn(N)) .* sqrt(Phi) * df; % 复高斯随机频谱 ph real(ifft2(ifftshift(H))) * N * N; % 逆傅里叶回空间域 % 三次次谐波补偿: 逐级用3x3, 9x9, 27x27子网格细分低频 for p 1:num_sub n_sub 3^p; fx_s (-1:1) * (df / n_sub); [FXs, FYs] meshgrid(fx_s); fs sqrt(FXs.^2 FYs.^2); fs(2,2) 1e-9; % 子网格中心同样置小值 Phis 0.033 * Cn2 * (fs.^2 1/L0^2).^(-11/6) .* exp(-fs.^2/kl^2); Hs (randn(3) 1i*randn(3)) .* sqrt(Phis) * (df/n_sub); ps real(ifft2(ifftshift(Hs))) * 9; % 3x3子网格相位 ps_up imresize(ps, [N, N], bilinear); % 插值放大回主网格 ph ph ps_up / (3^p); % 逐阶累加低频分量 end end先讲主函数部分。频率轴fx按(-N/2:N/2-1)/(N·dx)构造对MATLAB的ifftshift排版是匹配的零频在矩阵中心附近逆变换前需要ifftshift把零频搬回矩阵左上角。f(1,1)置成1e-9而不是0是因为0会让谱值变成无穷大后续ifft2没法处理这个中心点对应直流分量它不影响相位屏的起伏结构只影响整体活塞相位统计光斑时不关心。生成随机频谱时用了复高斯随机数实部和虚部各占一个自由度乘以sqrt(Phi)·df是对连续功率谱做离散采样的标准幅度映射df1/(N·dx)是频域网格间距。ifft2自带1/N²归一化所以前面乘N·N把幅度恢复回来。最后取real不取abs是因为复频谱的对称性约束下逆变换的虚部只是数值误差。3.3 次谐波补偿的权重、插值与标定次谐波权重是这套系统里最需要手工调的地方。不同代码对低频补偿系数的写法不统一有的分母用3^p有的用9^p有的根本不除。原因是前置频谱采样幅度的写法不一致有的乘了df有的乘了df²导致同样一段代码在不同版本之间表现差异很大。直接照搬别人的权重很可能低频过冲或不足这部分做成“可调参数”更稳妥。我拿到这类代码的第一件事不是跑传输而是生成20帧相位屏统计相位结构函数和理论值对比看低频段斜率和幅值对不对。插值方式也有讲究imresize用bilinear对3×3或9×9的相位做放大会引入一些平滑这个平滑能去除子网格边缘的台阶感但也会轻微压低低频幅度。补偿不足时可以把插值改成cubic或者先在子网格上做零填充再傅里叶变换效果更锐利。最后的相位方差标定是收尾工作先不乘任何缩放直接统计多帧相位屏的方差和期望方差比较后统一乘一个全局校正系数。这个方法不优雅但实用能把屏幕总能量先保证住再做低频段细节微调。4. 高斯光束多屏传输初始化、传播算子和光束质量评估4.1 高斯光束初始化和角谱传播算子传输仿真从基模高斯光束开始。假设z0处是束腰光场振幅分布写作Eexp(-r²/w0²)w0是束腰半径。初始化代码如下% 高斯光束初始化 x (-N/2:N/2-1) * dx; [X, Y] meshgrid(x); r2 X.^2 Y.^2; E exp(-r2 / w0^2); % 束腰处振幅, 暂不引入初始波前曲率 % 角谱传播算子(傍轴近似) fx (-N/2:N/2-1) / (N*dx); [FX, FY] meshgrid(fx); H exp(1i * pi * lambda * dz * (FX.^2 FY.^2)); % 单位: lambda和dz都换算成米w0的选择要照顾网格尺寸。w0太小光斑只占几个像素后面统计光斑半径没有意义w0太大网格边缘截断产生衍射环又会污染光斑。我一般让w0大概占网格物理长度的1/10到1/4。波长lambda常见取值有532 nm、1.064 μm、1.55 μm这套系统都可以跑只需要保持lambda、dx、dz单位一致。角谱传播算子H是自由空间传播在傍轴近似下的频域形式代表一段距离dz的真空衍射传播。在频域做乘法比时域卷积快得多而且单步传播本身不引入额外的奈奎斯特限制。但有个前提条件实际使用时要检查网格物理尺寸LxN·dx与dz之间需要满足lambda·dz小于Lx·dx的量级否则高频分量混叠光斑周围出现稳定散点。这个条件不满足时优先把N翻倍或者缩短单步传播距离。4.2 多屏循环与Cn2分配让每一层屏都贡献正确的扰动相位屏和传播算子的组合循环是整个仿真的主体。总路径分成M段每段dz z_total/M屏幕依次摆在各段末端。M不是越大越好屏幕太多单屏相位方差太小数值噪声占比上升屏幕太少单屏扰动过大背离薄屏假设。我一般用经验标准控制单屏相位方差在0.11 rad²量级用M反推路径分段。循环代码如下% screens为N x N x M的三维数组, 由vonkarman_screen批量生成 Es E; for k 1:M Es Es .* exp(1i * screens(:,:,k)); % 相位屏调制 Es ifft2(fft2(Es) .* H); % 真空传播一段dz end I abs(Es).^2; % 接收面光强每屏生成必须与dz匹配。常见的翻车做法是先把一条长路径的强湍流整体压进一个相位屏再把这个屏当成薄屏直接用结果到达角抖动比理论值大好几倍。正确做法是按段更新湍流强度。如果整条路径的Cn2视为常数就按该段路径单独生成屏如果某一段湍流特别强把该段单独细分弱湍流段合并而不是平均分配。screens三维数组的内存要提前算一下。N512、M10时一个screens数组约20 MB可以接受M50就该按段逐屏生成、逐屏传播不要把所有屏都先存在内存里再循环。很多机器跑512网格卡死不是算法问题是三维数组把内存占了。4.3 光束质量评估光斑、Strehl比和质心抖动传输结束后统计量集中在几个指标上光斑半径、Strehl比、质心位置和长曝光光斑。Strehl比定义为有湍流和无湍流时接收面峰值光强的比值是光束质量最直观的指标。质心抖动对应到达角起伏用多帧统计标准差表达。Iturb I; Iideal abs(ifft2(fft2(E0) .* H_total)).^2; % 无湍流传播结果, H_total为多段总算子 Strehl max(Iturb(:)) / max(Iideal(:)); % 光斑质心 cx sum(sum(I .* X)) / sum(I(:)); cy sum(sum(I .* Y)) / sum(I(:)); % 多帧重复后统计cx/cy的标准差即为到达角抖动注意Strehl的定义在不同文献里略有差别有的用峰值光强比有的用桶内功率比。我习惯把峰值比和63%环围能量半径一起报至少口径不会被质疑。质心抖动算出来后如果波长和传播距离已知可以换算成到达角方差再用理论公式对照数量级。数量级不对大概率是相位屏低频段没补好而不是传播代码写错了。5. 复现避坑这套仿真系统里最容易翻车的五个环节5.1 相位屏出现规则网格纹理现象生成的相位屏放大后能看到明显的棋盘状或条纹状纹理光斑图像边缘出现伪周期结构。原因主网格频域采样到的高频部分超过奈奎斯特频率高频能量折叠回低频。常见触发条件是dx相对内尺度l0太大或者N太小导致最大可表示频率不够覆盖l0对应的频率范围。解决保持dx不超过l0的1/3。如果内尺度很小、网格数又不能增大可以在生成相位屏时适当增大kl对应的截止值人为加快高频耗散让频谱在高频端自然衰减后再进入逆傅里叶变换。5.2 单屏相位方差和理论值对不上现象统计大量相位屏的方差均值要么是理论值的零点几倍要么是好几倍重复运行结果稳定偏大或偏小。原因频域抽样幅度写法不一致是最常见的原因有的代码里sqrt(Phi)·df有的乘df²再叠加ifft2的归一化最终结果天差地别。另一个原因是f(1,1)被置零后直流附近的低频能量没有参与统计屏幕总方差偏低。解决先不纠结公式写法生成30帧屏幕统计方差用一个全局校正系数把方差标定到理论值。然后把结构函数曲线拉出来看低频段斜率再决定是否调整次谐波权重。标定系数只解决总能量问题解决不了低频分布问题所以两步都要做。5.3 光斑抖动偏大或偏小现象质心抖动标准差和理论公式差两三倍而且不是随机波动重复跑都稳定偏移。原因多半是相位屏没有按段匹配Cn2把整段路径的湍流全部压到一个屏上或者M取值太小单屏扰动超出薄屏假设的适用范围。解决让M满足单屏相位方差0.11 rad²。路径上只有某些高度段湍流强时强湍流段单独细分弱湍流段可以合并。调整后再看到达角起伏数量级基本能对上。如果还差检查是否引入初始波前曲率束腰不在起始面时高斯光本身有发散角质心抖动统计会混入几何扩散的贡献。5.4 长距离传播后光斑出现高频散点现象传输距离一长光斑外围出现规则散点看起来像噪声但换随机种子后位置不变。原因角谱法的空间带宽积不够。网格物理尺寸LxN·dx角谱传播的空间频率范围受限于N/(2Lx)当lambda·dz超过Lx·dx量级高频分量折叠产生散点。解决增大N、增大dx或者缩短单步dz。如果dz受屏数限制无法缩短改用菲涅尔衍射积分加单次FFT的方式或者减少每段真空传播长度。注意这个坑和5.1的网格纹理视觉上很像但成因完全不同5.1是相位屏生成时的高频混叠这里是传播算子的高频混叠排查入口一个在屏生成端一个在传播端。5.5 次谐波叠加后低频过冲现象相位屏方差已经标定正确但光斑整体偏移量还是明显大于实测结构函数在低频段高于理论趋势线。原因次谐波每阶子网格的权重偏大。不同实现里分母写成3^p、9^p、3^(p/2)都有可能取决于插值方式和频域抽样写法。bilinear插值本身会压低高频细节有人为了拟合曲线把权重调大实际是把低频放过了。解决把次谐波权重改成逐阶可调的数组比如[1/3, 1/9, 1/27]先跑如果低频上翘继续压缩高阶权重。每改一次就看一次结构函数直到在r大于几个dx之后与理论差5%以内。这个参数不要一次性调到感觉上一定要用结构函数做判据。6. 进阶验证用结构函数和统计收敛性给结果“体检”6.1 单屏结构函数校验结构函数是相位屏质量最直接的体检指标定义是Dφ(r)⟨[φ(ρ)-φ(ρr)]²⟩表示距离r两点的相位差方差。对随机相位屏这个量在小r区间的斜率对应谱模型的特征。用20帧相位屏做估计% ph为vonkarman_screen生成的相位屏, 多帧循环后统计 dphi zeros(1, N/2 - 1); for r 1:N/2 - 1 diffv ph(:, 1r:end) - ph(:, 1:end-r); dphi(r) mean(diffv(:).^2, omitnan); end loglog((1:N/2-1)*dx, dphi);低频段r较大的区域斜率如果低于5/3说明次谐波补得不够如果高于5/3说明低频过冲。这个检查跑一次只要几秒钟却能省下整条传输链路返工的时间。我拿到这套系统时第一件事就是把N512、dx0.005 m、L020 m这组参数做成默认用上面的脚本标定确认无误后再开始跑蒙特卡洛。6.2 蒙特卡洛次数和统计量的收敛判断相位屏的随机性意味着单次仿真没有意义至少几十帧才能让质心抖动和Strehl比稳定下来。跑之前先设定总次数比如100帧然后边跑边看统计量的滑动平均是否收敛sr_avg movmean(Strehl_all, 10); plot(1:num_runs, sr_avg);Strehl比前20帧波动大很正常如果到80帧还在单调上升或下降说明某些帧的相位屏出现低频异常这时候回头检查结构函数而不是继续加帧。我还习惯同时监控质心抖动的累计标准差曲线它比Strehl收敛得更慢通常要到60帧以后才平稳。从那以后我每次把Cn2、L0或者传输距离一改都会先跑一遍这组体检脚本再去看光斑有没有翻车。这个习惯帮我把至少三次返工挡在了出图之前。希望帮到你。本文还有配套的精品资源点击获取