简介本资源是一套基于MATLAB/Simulink开发的热晕相位屏仿真程序面向光学工程、激光大气传输及自适应光学领域的科研人员与高年级本科生用于定量研究高能激光在大气中传播时由热晕效应引起的波前畸变。程序支持调节激光功率、光束参数、大气湍流强度及传输距离等关键变量可生成对应条件下的动态相位屏数据为后续波前校正算法验证提供基础输入。压缩包共含12个文件7个核心M脚本实现相位屏建模、FFT/IFFT运算与热晕迭代计算4幅BMP格式中间结果图像用于可视化验证1个ASV备份文件总容量541KB结构紧凑、模块职责明确便于二次开发与参数调试。目前已有219人学习下载配套代码注释清晰包含热晕积分模型reyunjisuanzz.m、相位屏重构xuanhuan.m、真空传输模拟zhenkongchuanshu.m等关键功能可直接运行并拓展至闭环补偿系统仿真。 做激光大气传输的人应该都躲不开“热晕”这两个字。无论是高能激光武器、自由空间光通信还是激光测距只要激光功率上去了大气吸收带来的热效应就会跳出来捣乱。我之前用MATLAB/Simulink写了一套热晕相位屏仿真程序可以模拟不同风速、不同吸收系数、不同功率条件下的大气热晕相位屏用来评估光束质量退化程度。今天把这个程序的核心思路、实现细节、参数选取和踩过的坑都整理出来给同样做相关方向的朋友参考。这套程序解决的核心问题是如何在仿真中把“大气被激光加热后造成的折射率不均匀”定量地表达出来。有了热晕相位屏就能把它直接嵌入到已有的光束传播仿真链路里配合真空衍射传播、湍流相位屏形成一套完整的大气激光传输仿真工具。适合做光传输仿真、自适应光学算法验证、系统参数论证的工程师和研究生使用。1. 热晕相位屏是什么先搞清楚要仿真什么1.1 热晕效应的物理本质热晕效应英文叫thermal blooming是指高功率激光在大气中传输时大气分子和气溶胶粒子吸收了一部分激光能量导致光路附近的空气温度升高。空气温度升高后密度就会下降折射率也随之改变。这样一来原本均匀的大气变成了一个“负透镜”——光束中心区域折射率低边缘折射率高光束被向外扩散焦斑能量不再集中。这个现象在连续激光、高功率脉冲串激光中表现尤其明显。低功率情况下大气对激光的吸收几乎不会引起可观测的折射率变化一切都可以按照线性光学来处理。但一旦功率密度达到一定阈值热晕效应就会显著影响光束质量斯特列尔比Strehl ratio会急剧下降光束中心可能出现空心化、破碎化的现象。这也是高能激光系统设计时绕不开的一个问题。1.2 为什么用相位屏模型大气对激光的影响可以分成两个维度振幅变化和相位变化。对于宏观尺度的大气吸收引起的热晕在单次通过大气的条件下振幅变化相对较小主要是折射率不均匀导致相位畸变。这在物理上非常适合用相位屏来等效。所谓相位屏就是把连续分布的大气折射率扰动压缩到一个极薄的平面上。光波在这个平面上传播时只改变相位、不改变振幅相当于乘以一个复振幅因子 exp(i·φ(x,y))。这其实是一个数学上的“薄透镜近似”——只要传播距离足够短折射率扰动足够弱这种近似就是可靠的。在实践中我们通常把整条大气路径分成若干段每段用一个相位屏代替该段的所有折射率扰动段与段之间用真空角谱传播来衔接。这就是经典的“分步传播法”split-step method。热晕相位屏就是其中描述热晕效应的那一部分相位扰动。1.3 这套程序能做什么我写的这套程序的定位是生成“纯热晕”的相位屏即只考虑大气吸收激光能量后引起的折射率变化。它跟常见的Kolmogorov湍流相位屏是分开的——湍流相位屏描述的是大气本身随机的温度/密度起伏而热晕相位屏描述的是激光“自己加热自己”造成的、空间分布有规律的热畸变。程序可以输出特定传播距离处的热晕相位分布、光强分布也可以输出等效的相位屏数据。这套程序的核心价值在于“参数可调”——改变激光功率、光束直径、大气吸收系数、横向风速、离焦量、传播距离等可以反复计算比较评估不同工况下热晕效应的严重程度为后续的光束控制策略设计提供依据。2. 工具选型为什么用MATLAB加Simulink2.1 MATLAB承担核心计算Simulink承担流程编排很多做物理仿真的朋友可能会问热晕相位屏是纯数值计算Simulink是动态系统建模工具两者有什么关系实际上在这个方案里MATLAB脚本负责所有核心计算——包括光束生成、热晕相位计算、分步传播、数据后处理Simulink模型则负责整体仿真流程的编排、参数扫描控制、数据采集和结果可视化。我在实际使用中发现这个分工非常合理。热晕相位屏的物理计算涉及多层循环和FFT运算这些在Simulink的模块框图里硬搭会非常痛苦而且可读性差、调试困难。而把核心计算封装成MATLAB Function模块放到Simulink里调用既能保留Simulink在系统级建模、参数管理方面的优势又能让物理计算部分保持代码的灵活性。2.2 整体方案架构整个仿真流程大概是这样的Simulink模型作为顶层框架包含参数初始化模块、传播控制模块、结果存储模块。其中传播控制模块按传播层数进行循环每一层内调用MATLAB Function完成相位屏生成和光场传播最终结果输出到MATLAB工作区再做后处理。具体来说Simulink端负责三件事一是管理所有输入参数的配置防止脚本参数杂乱二是通过计时器和使能信号控制多层循环的执行顺序三是收集每一层的中间结果方便观察热晕随传播距离的演化。MATLAB端负责热晕相位屏的计算细节包括温度场更新、折射率场计算、相位积分等。这种架构的另一个好处是后期扩展方便。如果你后续想加入湍流相位屏、跟踪误差、平台抖动等模型只需要在Simulink的传播链路里再串联一个函数模块即可不需要改动核心热晕代码。3. 核心公式与参数设计3.1 热晕相位屏的数学表达热晕相位屏的核心公式是基于流体力学和热力学方程的近似解。在连续激光、匀速横向风、准稳态假设下由热晕引起的折射率变化Δn(x, y, z)可以写成Δn(x, y, z) (dn/dT) · (α·I(x, y, z) · t_c) / (ρ·C_p)其中dn/dT是空气折射率温度系数大约为-1×10⁻⁶ K⁻¹量级α是大气吸收系数单位是m⁻¹I是激光光强分布t_c是特征时间在横向风存在时t_c D/vD是光束直径v是横向风速ρ是空气密度C_p是定压比热容。这里要说明一下这个公式是在“准稳态、对流主导”条件下的简化模型。它假设横向风把被加热的空气吹走达到一个动态平衡而不是无限加热下去。如果没有横向风或者风速极低就必须求解含时热传导方程计算量会大很多。得到Δn之后相位屏的计算就很简单了Δφ(x, y) (2π/λ) · Δn(x, y) · Δzλ是激光波长Δz是这一层大气的厚度。把这个相位叠加到光场上就完成了热晕相位屏的核心功能。3.2 关键无量纲参数在调试程序之前我强烈建议你先手算一下几个无量纲参数判断你关心的工况处于什么区间。第一个是热畸变数N_D它定义为N_D (2√2 · k · L · (dn/dT) · α · P) / (π · n_0 · ρ · C_p · v · D³)其中k是波数2π/λL是传播距离P是激光总功率n_0是大气初始折射率D是光束直径。这个参数衡量热晕效应的强度大致公式。N_D越大热晕越严重。根据经验N_D小于1的工况热晕效应可忽略N_D在1到10之间时热晕明显大于10时光束质量会严重退化。第二个是相变参数或者叫Streck参数它描述了热晕相位畸变的空间尺度与光束直径的比值。这个参数影响相位屏功率谱的形状。在实际写代码时我习惯先算N_D再决定是否需要跑完整仿真。如果N_D远小于1直接跳过热晕仿真只做真空传播和湍流仿真就够了能省很多计算时间。3.3 网格与采样参数选择网格参数选得好不好直接决定仿真结果可不可信。我常用的经验是网格点数N取256或512网格间距Δx根据光束直径D来定保证光束直径占到网格宽度的1/3到1/2。也就是说如果D0.1 m网格宽度L_grid取0.2到0.3 mΔx就是L_grid/N。为什么不能把网格开得太大或者太小网格太小时光场的边缘会被截断产生明显的衍射伪影网格太大时远场采样的角度分辨率不够焦平面光斑细节出不来。经验法则是网格宽度至少是光束直径的2倍留足衍射扩展的空间。另一个关键参数是传播层数。总传播距离L除以每层厚度Δz就是层数。热晕相位屏每层厚度怎么选我的经验是每层中N_D的增量控制在0.1到0.3之间。这样既能捕捉热晕随距离的累积过程又不至于因为层数太多导致计算时间爆炸。比如总N_D在3左右时取10到20层是合理的。4. 程序实现细节4.1 热晕相位屏生成主函数我把热晕相位屏的生成封装成一个独立的函数输入参数包括光场分布、大气参数、传播层厚度输出是相位屏矩阵。函数内部按照3.1节的公式一步步计算。function phase_screen thermal_blooming_screen(E_in, param) % 输入 % E_in: 当前层入射光场NxN复数矩阵 % param: 结构体包含所有大气与光束参数 % 输出 % phase_screen: 热晕相位屏NxN实数矩阵单位rad N param.N; I abs(E_in).^2; % 光强分布 I_mean mean(I(:)); % 平均光强用于归一化 % 对流热扩散的简化模型 % 温度增量正比于光强分布响应函数是低通滤波 delta_T param.alpha * I / (param.rho * param.Cp) * param.t_Char; delta_T imgaussfilt(delta_T, param.filter_sigma); % 考虑热扩散平滑 % 折射率变化 delta_n param.dn_dT * delta_T; % 相位累积2π/λ * Δn * Δz phase_screen (2*pi/param.lambda) * delta_n * param.dz; end这里的imgaussfilt是MATLAB自带的高斯滤波函数用来模拟热扩散的平滑效果。这个平滑的半径不能随意取它对应热扩散的特征尺度。在风主导流动时扩散尺度通常远小于光束直径所以平滑半径一般取1到2个网格点相当于一个数值稳定性的处理而不是物理建模的核心。4.2 多层传播循环有了相位屏生成函数接下来就是多层循环传播。每一层内做两件事先乘相位屏再做角谱传播。% 初始化光场 [x, y] meshgrid(param.x, param.y); r2 x.^2 y.^2; E param.E0 .* exp(-r2/param.w0^2); % 初始高斯光束 for k 1:param.n_layers % 1. 当前层热晕相位屏 phi_tb thermal_blooming_screen(E, param); % 2. 叠加相位屏 E E .* exp(1i * phi_tb); % 3. 真空角谱传播一层 E angular_spectrum_propagate(E, param.dz, param.lambda, param.dx); % 记录中间结果可选 record_intensity{k} abs(E).^2; end角谱传播函数可以用MATLAB的fft2和ifft2实现function E_out angular_spectrum_propagate(E_in, dz, lambda, dx) N size(E_in, 1); fx (-N/2 : N/2-1) / (N*dx); [FX, FY] meshgrid(fx, fx); H exp(1i * 2*pi/lambda * dz .* sqrt(1 - (lambda*FX).^2 - (lambda*FY).^2)); E_out ifft2(fft2(E_in) .* ifftshift(H)); end这个角谱传播函数注意两点一是频域坐标必须用fftshift/ifftshift配对好否则相位会错位二是对于倏逝波分量lambda·FX)² (lambda·FY)² 1的部分需要把H设为0否则会出现数值爆炸。4.3 Simulink模型搭建思路Simulink端我用的是“脚本驱动模型配置”方式而不是把所有计算都拖成模块。模型本身比较简洁主要有三层第一层是参数配置层用Simulink的常量模块或者通过set_param从MATLAB工作区读取参数。这样改参数不需要打开模型直接在脚本里改一个结构体就行。第二层是核心循环层用MATLAB Function模块承载多层传播循环。整个循环放在一个函数里Simulink只需要调用一次这个函数就完成了整个传播过程的计算。第三层是结果输出层把仿真结果写到工作区、保存成.mat文件或者用Scope模块显示光斑演化过程。有人可能会说这样Simulink参与的部分太少了有点形式主义。但实际上Simulink在这里的价值是提供一个“可视化装配界面”让你可以方便地把热晕模块、湍流模块、瞄准误差模块、接收端模块等串联起来形成完整的系统级仿真链路。这在做方案论证时非常有用——直接改连线就能改变系统拓扑。5. 不同工况仿真结果与规律5.1 不同横向风速的影响我跑的第一组对照实验是改变横向风速。风速从2 m/s到10 m/s变化功率固定为10 kW吸收系数取1e-4 m⁻¹传播距离2 km光束直径10 cm。风速对热晕的影响非常直接。低风速2 m/s时热畸变数N_D大焦斑严重扩展光斑呈现明显的“新月形”或“拖尾”形态。这是热晕的标志性特征——被加热的空气没有及时被吹走光束下风向一侧的畸变更严重。风速提高到8 m/s以上时N_D降到3以下焦斑质量明显改善斯特列尔比从0.2左右提升到0.6以上。这个规律说明一个工程问题如果平台本身有相对运动比如机载或车载激光器前进方向的空速本身就能抑制热晕。这对系统设计很有利。值得注意的一个细节是风的方向不能忽略。如果风速方向与光束传播方向垂直抑制作用最明显如果近似平行则要考虑光束路径上风速分量的变化。我的仿真程序里风方向是通过在温度场计算时引入偏移来模拟的简单而有效。5.2 不同功率与吸收系数的耦合效应第二组实验是改变激光功率和大气吸收系数。功率从5 kW到50 kW变化吸收系数取两个典型值1e-4 m⁻¹晴天干净大气和5e-4 m⁻¹有霾或沙尘。结果很清楚功率翻倍N_D近似翻倍热晕几乎线性加重。吸收系数提高5倍N_D也提高5倍而且相位屏的峰值相位也会成比例增大。这表明在吸收严重的大气条件下不太适合用太高的连续功率否则能量都会浪费在加热空气上而不是到达目标。这组实验让我养成一个习惯跑仿真前先查该波长、该大气条件下的典型吸收系数。对1.064 μm波长晴天海平面大气吸收系数大约1e-4到3e-4 m⁻¹对3.8 μm中红外波段吸收系数明显更低对10.6 μm CO₂激光吸收系数略高。选不同波长时程序里的α参数要对应修改。5.3 有无横向风、有无湍流叠加的对比第三组实验是组合工况。第一种工况无风、无湍流只考虑热晕第二种工况有风、无湍流第三种工况有风、有中等强度湍流。第一种工况下热晕相位屏轴向对称焦斑均匀扩展类似一个大的负透镜加球差。光束中心强度下降边缘出现环状结构但整体对称性保持较好。第二种工况下相位屏出现不对称的马鞍形分布焦斑移向风的反方向出现明显的彗差特征。第三种工况叠加了湍流情况复杂得多。湍流相位屏引入随机的高频畸变热晕引入低频的系统性畸变两者的叠加效果不是简单相加。在数值仿真中必须先用湍流相位屏更新光场再用热晕相位屏更新两者的先后顺序会影响结果。我在程序中按“湍流加在每层传播的前半段热晕加在后半段”的方式来处理这是对实际物理过程的一种近似——湍流是环境固有的热晕是光束自致的。6. 常见问题与排查技巧6.1 相位屏出现“网格状”伪影这是初期最容易遇到的问题之一。相位屏在高频部分出现规律性的网格纹路通常对应频域坐标处理错误尤其是fftshift和ifftshift的配对问题。排查方法很简单打印出角谱传播函数H的前几个点检查中心点的值是否接近1。如果H的中心不在数组中心说明fftshift没有配对好。另一个常见原因是网格宽度选取不当导致高频分量被截断产生Gibbs现象。解决方案是统一在频域函数生成时用fftshift在乘法之前用ifftshift(FX, FY)。我习惯把这四行写成固定模板不再改动避免每次重写踩坑。6.2 仿真结果发散或出现NaN出现NaN通常是光场某些点强度为0但相位却巨大经指数运算后溢出。热晕相位屏在光强接近0的边缘区域理论上相位增量也趋近0但如果数值上出现0除以某个小量就可能导致NaN。解决方法是给光强加一个很小的地板值floor比如I_safe max(I, 1e-12)或者在计算δT前对光场做一次高斯滤波。这个操作物理上是合理的——实际光束不会严格到0强度背景噪声和散射光总会提供一个下限。6.3 不同分辨率下结果不一致如果你的仿真结果强烈依赖于网格点数那说明网格参数没有收敛。N从128增加到256时结果变化很大再增加到512才基本稳定。这种情况我建议你用“双网格验证”先用粗网格快速扫描参数再对关键工况用细网格精确计算。还有一个容易被忽略的因素是光束直径占网格宽度的比例。如果网格宽度开得太大比如光束直径只占网格宽度的1/10虽然光场不会被截断但远场的角度采样间隔太大焦点处的光斑细节会丢失。这时需要增加网格点数来同时保证空间分辨率和角度分辨率。6.4 常见问题速查表现象可能原因排查与解决相位屏网格状伪影fftshift配对错误、频域坐标错位检查H矩阵中心统一ifftshift模板结果发散或NaN光强过小导致相位溢出给光强加地板值或做平滑处理不同网格点数误差大网格未收敛、光束占比不匹配做双网格验证调整网格宽度焦斑不对称且不可复现随机噪声源未固定设置随机种子固定初始条件热晕效果被低估层数太少、N_D增量过大增加层数确保单层N_D增量小于0.3Simulink调用速度慢每步都在与MATLAB工作区交互批量传入参数减少set_param调用6.5 Simulink仿真速度优化Simulink调用MATLAB Function时如果参数传递方式不当速度会慢得让人崩溃。我实测的经验是在模型里尽量不要用From Workspace模块去读每个中间变量而是把所有参数打包成一个结构体在模型初始化时一次性传入。另一个速度瓶颈是代码生成。如果你安装了MATLAB Coder可以把热晕相位屏函数单独生成MEX文件Simulink里直接调用MEX版本。我的实际测试里MEX化后速度提升了5到10倍256×256网格、20层传播的计算时间从秒级降到亚秒级。7. 这套程序的扩展方向与实际使用建议7.1 扩展方向叠加湍流、指向抖动与自适应光学热晕相位屏程序本身只是一个模块但它能扩展出很多实用功能。我在项目中主要做了三个扩展第一个是叠加湍流相位屏。把热晕相位屏和Kolmogorov湍流相位屏都叠加到同一个传播层里可以考察热晕效应和湍流效应的耦合作用。这在高功率激光传输中很关键因为湍流会改变光束的局域强度分布进而影响热晕的分布而热晕又会反过来影响湍流的统计特征两者存在非线性耦合。第二个是叠加平台指向抖动。在Simulink里加入一个随机指向误差模块每次传播前对光场施加一个小角度的倾斜相位模拟跟踪系统的残余误差。这能评估系统对瞄准稳定度的需求。第三个是接入自适应光学模块。如果热晕导致的低阶像差是主要问题Shack-Hartmann波前传感器加变形镜的控制环路可以显著改善焦斑质量。Simulink非常适合做这种“光学控制”联合仿真这也是我坚持用Simulink而不是纯脚本来搭框架的核心原因。7.2 实际使用建议如果你的机器配置一般我建议先用256×256网格做快速预研确定参数范围后再用512×512网格精算。不要一开始就上1024×1024否则一次仿真跑几分钟参数扫描根本做不动。参数的保存和追溯也很重要。我习惯给每次仿真建立一个独立的参数结构体存成.mat文件文件名带时间戳和关键参数摘要比如TB_10kW_v5ms_alpha1e-4_N256.mat。这样回看结果时不会搞混哪组参数对应哪个结果。7.3 一个容易犯的错误忽略瞬态过程很多简化模型的公式都基于准稳态假设也就是认为温度场已经稳定了。但实际中激光发射的初期有一个瞬态过程热晕还没有完全建立此时的光斑质量和稳态时有明显差别。如果你的需求是仿真脉冲激光或短时间发射的场景必须在程序中加入时间累积项用时间步进的方式求解温度场而不是直接用稳态公式。我在程序中加了一个可选的时间演化模式使用显式欧拉法逐步更新温度场可以按需开启。7.4 最后说些真心话做了这么多组仿真之后我个人最深的体会是热晕仿真的难点不在“算出一个相位屏”而在于“验证这个相位屏物理上合理”。我在初步完成程序后对照了公开发表的实验数据发现单纯用解析近似公式计算的相位屏往往高估了热晕效应——原因在于真实大气中还存在湍流混合、气流起伏等附加机制会部分缓解热聚焦效应。所以如果你要用这套程序做工程论证千万不要拿单一工况的仿真结果作为最终结论。至少要跑一个参数矩阵看看系统在好天气、坏天气、不同风速下的表现包络再结合统计学方法给出概率性的评估。仿真程序的价值在于帮我们理解趋势和机理而不在于提供一个精确到小数点的绝对数值。这些年我用它做过方案对比、系统参数优化、甚至写了小半篇论文的分析部分希望这些实现细节也能让你少走几步弯路。本文还有配套的精品资源点击获取 基于MATLAB/Simulink的大气热晕相位屏仿真实现与参数分析
简介本资源是一套基于MATLAB/Simulink开发的热晕相位屏仿真程序面向光学工程、激光大气传输及自适应光学领域的科研人员与高年级本科生用于定量研究高能激光在大气中传播时由热晕效应引起的波前畸变。程序支持调节激光功率、光束参数、大气湍流强度及传输距离等关键变量可生成对应条件下的动态相位屏数据为后续波前校正算法验证提供基础输入。压缩包共含12个文件7个核心M脚本实现相位屏建模、FFT/IFFT运算与热晕迭代计算4幅BMP格式中间结果图像用于可视化验证1个ASV备份文件总容量541KB结构紧凑、模块职责明确便于二次开发与参数调试。目前已有219人学习下载配套代码注释清晰包含热晕积分模型reyunjisuanzz.m、相位屏重构xuanhuan.m、真空传输模拟zhenkongchuanshu.m等关键功能可直接运行并拓展至闭环补偿系统仿真。 做激光大气传输的人应该都躲不开“热晕”这两个字。无论是高能激光武器、自由空间光通信还是激光测距只要激光功率上去了大气吸收带来的热效应就会跳出来捣乱。我之前用MATLAB/Simulink写了一套热晕相位屏仿真程序可以模拟不同风速、不同吸收系数、不同功率条件下的大气热晕相位屏用来评估光束质量退化程度。今天把这个程序的核心思路、实现细节、参数选取和踩过的坑都整理出来给同样做相关方向的朋友参考。这套程序解决的核心问题是如何在仿真中把“大气被激光加热后造成的折射率不均匀”定量地表达出来。有了热晕相位屏就能把它直接嵌入到已有的光束传播仿真链路里配合真空衍射传播、湍流相位屏形成一套完整的大气激光传输仿真工具。适合做光传输仿真、自适应光学算法验证、系统参数论证的工程师和研究生使用。1. 热晕相位屏是什么先搞清楚要仿真什么1.1 热晕效应的物理本质热晕效应英文叫thermal blooming是指高功率激光在大气中传输时大气分子和气溶胶粒子吸收了一部分激光能量导致光路附近的空气温度升高。空气温度升高后密度就会下降折射率也随之改变。这样一来原本均匀的大气变成了一个“负透镜”——光束中心区域折射率低边缘折射率高光束被向外扩散焦斑能量不再集中。这个现象在连续激光、高功率脉冲串激光中表现尤其明显。低功率情况下大气对激光的吸收几乎不会引起可观测的折射率变化一切都可以按照线性光学来处理。但一旦功率密度达到一定阈值热晕效应就会显著影响光束质量斯特列尔比Strehl ratio会急剧下降光束中心可能出现空心化、破碎化的现象。这也是高能激光系统设计时绕不开的一个问题。1.2 为什么用相位屏模型大气对激光的影响可以分成两个维度振幅变化和相位变化。对于宏观尺度的大气吸收引起的热晕在单次通过大气的条件下振幅变化相对较小主要是折射率不均匀导致相位畸变。这在物理上非常适合用相位屏来等效。所谓相位屏就是把连续分布的大气折射率扰动压缩到一个极薄的平面上。光波在这个平面上传播时只改变相位、不改变振幅相当于乘以一个复振幅因子 exp(i·φ(x,y))。这其实是一个数学上的“薄透镜近似”——只要传播距离足够短折射率扰动足够弱这种近似就是可靠的。在实践中我们通常把整条大气路径分成若干段每段用一个相位屏代替该段的所有折射率扰动段与段之间用真空角谱传播来衔接。这就是经典的“分步传播法”split-step method。热晕相位屏就是其中描述热晕效应的那一部分相位扰动。1.3 这套程序能做什么我写的这套程序的定位是生成“纯热晕”的相位屏即只考虑大气吸收激光能量后引起的折射率变化。它跟常见的Kolmogorov湍流相位屏是分开的——湍流相位屏描述的是大气本身随机的温度/密度起伏而热晕相位屏描述的是激光“自己加热自己”造成的、空间分布有规律的热畸变。程序可以输出特定传播距离处的热晕相位分布、光强分布也可以输出等效的相位屏数据。这套程序的核心价值在于“参数可调”——改变激光功率、光束直径、大气吸收系数、横向风速、离焦量、传播距离等可以反复计算比较评估不同工况下热晕效应的严重程度为后续的光束控制策略设计提供依据。2. 工具选型为什么用MATLAB加Simulink2.1 MATLAB承担核心计算Simulink承担流程编排很多做物理仿真的朋友可能会问热晕相位屏是纯数值计算Simulink是动态系统建模工具两者有什么关系实际上在这个方案里MATLAB脚本负责所有核心计算——包括光束生成、热晕相位计算、分步传播、数据后处理Simulink模型则负责整体仿真流程的编排、参数扫描控制、数据采集和结果可视化。我在实际使用中发现这个分工非常合理。热晕相位屏的物理计算涉及多层循环和FFT运算这些在Simulink的模块框图里硬搭会非常痛苦而且可读性差、调试困难。而把核心计算封装成MATLAB Function模块放到Simulink里调用既能保留Simulink在系统级建模、参数管理方面的优势又能让物理计算部分保持代码的灵活性。2.2 整体方案架构整个仿真流程大概是这样的Simulink模型作为顶层框架包含参数初始化模块、传播控制模块、结果存储模块。其中传播控制模块按传播层数进行循环每一层内调用MATLAB Function完成相位屏生成和光场传播最终结果输出到MATLAB工作区再做后处理。具体来说Simulink端负责三件事一是管理所有输入参数的配置防止脚本参数杂乱二是通过计时器和使能信号控制多层循环的执行顺序三是收集每一层的中间结果方便观察热晕随传播距离的演化。MATLAB端负责热晕相位屏的计算细节包括温度场更新、折射率场计算、相位积分等。这种架构的另一个好处是后期扩展方便。如果你后续想加入湍流相位屏、跟踪误差、平台抖动等模型只需要在Simulink的传播链路里再串联一个函数模块即可不需要改动核心热晕代码。3. 核心公式与参数设计3.1 热晕相位屏的数学表达热晕相位屏的核心公式是基于流体力学和热力学方程的近似解。在连续激光、匀速横向风、准稳态假设下由热晕引起的折射率变化Δn(x, y, z)可以写成Δn(x, y, z) (dn/dT) · (α·I(x, y, z) · t_c) / (ρ·C_p)其中dn/dT是空气折射率温度系数大约为-1×10⁻⁶ K⁻¹量级α是大气吸收系数单位是m⁻¹I是激光光强分布t_c是特征时间在横向风存在时t_c D/vD是光束直径v是横向风速ρ是空气密度C_p是定压比热容。这里要说明一下这个公式是在“准稳态、对流主导”条件下的简化模型。它假设横向风把被加热的空气吹走达到一个动态平衡而不是无限加热下去。如果没有横向风或者风速极低就必须求解含时热传导方程计算量会大很多。得到Δn之后相位屏的计算就很简单了Δφ(x, y) (2π/λ) · Δn(x, y) · Δzλ是激光波长Δz是这一层大气的厚度。把这个相位叠加到光场上就完成了热晕相位屏的核心功能。3.2 关键无量纲参数在调试程序之前我强烈建议你先手算一下几个无量纲参数判断你关心的工况处于什么区间。第一个是热畸变数N_D它定义为N_D (2√2 · k · L · (dn/dT) · α · P) / (π · n_0 · ρ · C_p · v · D³)其中k是波数2π/λL是传播距离P是激光总功率n_0是大气初始折射率D是光束直径。这个参数衡量热晕效应的强度大致公式。N_D越大热晕越严重。根据经验N_D小于1的工况热晕效应可忽略N_D在1到10之间时热晕明显大于10时光束质量会严重退化。第二个是相变参数或者叫Streck参数它描述了热晕相位畸变的空间尺度与光束直径的比值。这个参数影响相位屏功率谱的形状。在实际写代码时我习惯先算N_D再决定是否需要跑完整仿真。如果N_D远小于1直接跳过热晕仿真只做真空传播和湍流仿真就够了能省很多计算时间。3.3 网格与采样参数选择网格参数选得好不好直接决定仿真结果可不可信。我常用的经验是网格点数N取256或512网格间距Δx根据光束直径D来定保证光束直径占到网格宽度的1/3到1/2。也就是说如果D0.1 m网格宽度L_grid取0.2到0.3 mΔx就是L_grid/N。为什么不能把网格开得太大或者太小网格太小时光场的边缘会被截断产生明显的衍射伪影网格太大时远场采样的角度分辨率不够焦平面光斑细节出不来。经验法则是网格宽度至少是光束直径的2倍留足衍射扩展的空间。另一个关键参数是传播层数。总传播距离L除以每层厚度Δz就是层数。热晕相位屏每层厚度怎么选我的经验是每层中N_D的增量控制在0.1到0.3之间。这样既能捕捉热晕随距离的累积过程又不至于因为层数太多导致计算时间爆炸。比如总N_D在3左右时取10到20层是合理的。4. 程序实现细节4.1 热晕相位屏生成主函数我把热晕相位屏的生成封装成一个独立的函数输入参数包括光场分布、大气参数、传播层厚度输出是相位屏矩阵。函数内部按照3.1节的公式一步步计算。function phase_screen thermal_blooming_screen(E_in, param) % 输入 % E_in: 当前层入射光场NxN复数矩阵 % param: 结构体包含所有大气与光束参数 % 输出 % phase_screen: 热晕相位屏NxN实数矩阵单位rad N param.N; I abs(E_in).^2; % 光强分布 I_mean mean(I(:)); % 平均光强用于归一化 % 对流热扩散的简化模型 % 温度增量正比于光强分布响应函数是低通滤波 delta_T param.alpha * I / (param.rho * param.Cp) * param.t_Char; delta_T imgaussfilt(delta_T, param.filter_sigma); % 考虑热扩散平滑 % 折射率变化 delta_n param.dn_dT * delta_T; % 相位累积2π/λ * Δn * Δz phase_screen (2*pi/param.lambda) * delta_n * param.dz; end这里的imgaussfilt是MATLAB自带的高斯滤波函数用来模拟热扩散的平滑效果。这个平滑的半径不能随意取它对应热扩散的特征尺度。在风主导流动时扩散尺度通常远小于光束直径所以平滑半径一般取1到2个网格点相当于一个数值稳定性的处理而不是物理建模的核心。4.2 多层传播循环有了相位屏生成函数接下来就是多层循环传播。每一层内做两件事先乘相位屏再做角谱传播。% 初始化光场 [x, y] meshgrid(param.x, param.y); r2 x.^2 y.^2; E param.E0 .* exp(-r2/param.w0^2); % 初始高斯光束 for k 1:param.n_layers % 1. 当前层热晕相位屏 phi_tb thermal_blooming_screen(E, param); % 2. 叠加相位屏 E E .* exp(1i * phi_tb); % 3. 真空角谱传播一层 E angular_spectrum_propagate(E, param.dz, param.lambda, param.dx); % 记录中间结果可选 record_intensity{k} abs(E).^2; end角谱传播函数可以用MATLAB的fft2和ifft2实现function E_out angular_spectrum_propagate(E_in, dz, lambda, dx) N size(E_in, 1); fx (-N/2 : N/2-1) / (N*dx); [FX, FY] meshgrid(fx, fx); H exp(1i * 2*pi/lambda * dz .* sqrt(1 - (lambda*FX).^2 - (lambda*FY).^2)); E_out ifft2(fft2(E_in) .* ifftshift(H)); end这个角谱传播函数注意两点一是频域坐标必须用fftshift/ifftshift配对好否则相位会错位二是对于倏逝波分量lambda·FX)² (lambda·FY)² 1的部分需要把H设为0否则会出现数值爆炸。4.3 Simulink模型搭建思路Simulink端我用的是“脚本驱动模型配置”方式而不是把所有计算都拖成模块。模型本身比较简洁主要有三层第一层是参数配置层用Simulink的常量模块或者通过set_param从MATLAB工作区读取参数。这样改参数不需要打开模型直接在脚本里改一个结构体就行。第二层是核心循环层用MATLAB Function模块承载多层传播循环。整个循环放在一个函数里Simulink只需要调用一次这个函数就完成了整个传播过程的计算。第三层是结果输出层把仿真结果写到工作区、保存成.mat文件或者用Scope模块显示光斑演化过程。有人可能会说这样Simulink参与的部分太少了有点形式主义。但实际上Simulink在这里的价值是提供一个“可视化装配界面”让你可以方便地把热晕模块、湍流模块、瞄准误差模块、接收端模块等串联起来形成完整的系统级仿真链路。这在做方案论证时非常有用——直接改连线就能改变系统拓扑。5. 不同工况仿真结果与规律5.1 不同横向风速的影响我跑的第一组对照实验是改变横向风速。风速从2 m/s到10 m/s变化功率固定为10 kW吸收系数取1e-4 m⁻¹传播距离2 km光束直径10 cm。风速对热晕的影响非常直接。低风速2 m/s时热畸变数N_D大焦斑严重扩展光斑呈现明显的“新月形”或“拖尾”形态。这是热晕的标志性特征——被加热的空气没有及时被吹走光束下风向一侧的畸变更严重。风速提高到8 m/s以上时N_D降到3以下焦斑质量明显改善斯特列尔比从0.2左右提升到0.6以上。这个规律说明一个工程问题如果平台本身有相对运动比如机载或车载激光器前进方向的空速本身就能抑制热晕。这对系统设计很有利。值得注意的一个细节是风的方向不能忽略。如果风速方向与光束传播方向垂直抑制作用最明显如果近似平行则要考虑光束路径上风速分量的变化。我的仿真程序里风方向是通过在温度场计算时引入偏移来模拟的简单而有效。5.2 不同功率与吸收系数的耦合效应第二组实验是改变激光功率和大气吸收系数。功率从5 kW到50 kW变化吸收系数取两个典型值1e-4 m⁻¹晴天干净大气和5e-4 m⁻¹有霾或沙尘。结果很清楚功率翻倍N_D近似翻倍热晕几乎线性加重。吸收系数提高5倍N_D也提高5倍而且相位屏的峰值相位也会成比例增大。这表明在吸收严重的大气条件下不太适合用太高的连续功率否则能量都会浪费在加热空气上而不是到达目标。这组实验让我养成一个习惯跑仿真前先查该波长、该大气条件下的典型吸收系数。对1.064 μm波长晴天海平面大气吸收系数大约1e-4到3e-4 m⁻¹对3.8 μm中红外波段吸收系数明显更低对10.6 μm CO₂激光吸收系数略高。选不同波长时程序里的α参数要对应修改。5.3 有无横向风、有无湍流叠加的对比第三组实验是组合工况。第一种工况无风、无湍流只考虑热晕第二种工况有风、无湍流第三种工况有风、有中等强度湍流。第一种工况下热晕相位屏轴向对称焦斑均匀扩展类似一个大的负透镜加球差。光束中心强度下降边缘出现环状结构但整体对称性保持较好。第二种工况下相位屏出现不对称的马鞍形分布焦斑移向风的反方向出现明显的彗差特征。第三种工况叠加了湍流情况复杂得多。湍流相位屏引入随机的高频畸变热晕引入低频的系统性畸变两者的叠加效果不是简单相加。在数值仿真中必须先用湍流相位屏更新光场再用热晕相位屏更新两者的先后顺序会影响结果。我在程序中按“湍流加在每层传播的前半段热晕加在后半段”的方式来处理这是对实际物理过程的一种近似——湍流是环境固有的热晕是光束自致的。6. 常见问题与排查技巧6.1 相位屏出现“网格状”伪影这是初期最容易遇到的问题之一。相位屏在高频部分出现规律性的网格纹路通常对应频域坐标处理错误尤其是fftshift和ifftshift的配对问题。排查方法很简单打印出角谱传播函数H的前几个点检查中心点的值是否接近1。如果H的中心不在数组中心说明fftshift没有配对好。另一个常见原因是网格宽度选取不当导致高频分量被截断产生Gibbs现象。解决方案是统一在频域函数生成时用fftshift在乘法之前用ifftshift(FX, FY)。我习惯把这四行写成固定模板不再改动避免每次重写踩坑。6.2 仿真结果发散或出现NaN出现NaN通常是光场某些点强度为0但相位却巨大经指数运算后溢出。热晕相位屏在光强接近0的边缘区域理论上相位增量也趋近0但如果数值上出现0除以某个小量就可能导致NaN。解决方法是给光强加一个很小的地板值floor比如I_safe max(I, 1e-12)或者在计算δT前对光场做一次高斯滤波。这个操作物理上是合理的——实际光束不会严格到0强度背景噪声和散射光总会提供一个下限。6.3 不同分辨率下结果不一致如果你的仿真结果强烈依赖于网格点数那说明网格参数没有收敛。N从128增加到256时结果变化很大再增加到512才基本稳定。这种情况我建议你用“双网格验证”先用粗网格快速扫描参数再对关键工况用细网格精确计算。还有一个容易被忽略的因素是光束直径占网格宽度的比例。如果网格宽度开得太大比如光束直径只占网格宽度的1/10虽然光场不会被截断但远场的角度采样间隔太大焦点处的光斑细节会丢失。这时需要增加网格点数来同时保证空间分辨率和角度分辨率。6.4 常见问题速查表现象可能原因排查与解决相位屏网格状伪影fftshift配对错误、频域坐标错位检查H矩阵中心统一ifftshift模板结果发散或NaN光强过小导致相位溢出给光强加地板值或做平滑处理不同网格点数误差大网格未收敛、光束占比不匹配做双网格验证调整网格宽度焦斑不对称且不可复现随机噪声源未固定设置随机种子固定初始条件热晕效果被低估层数太少、N_D增量过大增加层数确保单层N_D增量小于0.3Simulink调用速度慢每步都在与MATLAB工作区交互批量传入参数减少set_param调用6.5 Simulink仿真速度优化Simulink调用MATLAB Function时如果参数传递方式不当速度会慢得让人崩溃。我实测的经验是在模型里尽量不要用From Workspace模块去读每个中间变量而是把所有参数打包成一个结构体在模型初始化时一次性传入。另一个速度瓶颈是代码生成。如果你安装了MATLAB Coder可以把热晕相位屏函数单独生成MEX文件Simulink里直接调用MEX版本。我的实际测试里MEX化后速度提升了5到10倍256×256网格、20层传播的计算时间从秒级降到亚秒级。7. 这套程序的扩展方向与实际使用建议7.1 扩展方向叠加湍流、指向抖动与自适应光学热晕相位屏程序本身只是一个模块但它能扩展出很多实用功能。我在项目中主要做了三个扩展第一个是叠加湍流相位屏。把热晕相位屏和Kolmogorov湍流相位屏都叠加到同一个传播层里可以考察热晕效应和湍流效应的耦合作用。这在高功率激光传输中很关键因为湍流会改变光束的局域强度分布进而影响热晕的分布而热晕又会反过来影响湍流的统计特征两者存在非线性耦合。第二个是叠加平台指向抖动。在Simulink里加入一个随机指向误差模块每次传播前对光场施加一个小角度的倾斜相位模拟跟踪系统的残余误差。这能评估系统对瞄准稳定度的需求。第三个是接入自适应光学模块。如果热晕导致的低阶像差是主要问题Shack-Hartmann波前传感器加变形镜的控制环路可以显著改善焦斑质量。Simulink非常适合做这种“光学控制”联合仿真这也是我坚持用Simulink而不是纯脚本来搭框架的核心原因。7.2 实际使用建议如果你的机器配置一般我建议先用256×256网格做快速预研确定参数范围后再用512×512网格精算。不要一开始就上1024×1024否则一次仿真跑几分钟参数扫描根本做不动。参数的保存和追溯也很重要。我习惯给每次仿真建立一个独立的参数结构体存成.mat文件文件名带时间戳和关键参数摘要比如TB_10kW_v5ms_alpha1e-4_N256.mat。这样回看结果时不会搞混哪组参数对应哪个结果。7.3 一个容易犯的错误忽略瞬态过程很多简化模型的公式都基于准稳态假设也就是认为温度场已经稳定了。但实际中激光发射的初期有一个瞬态过程热晕还没有完全建立此时的光斑质量和稳态时有明显差别。如果你的需求是仿真脉冲激光或短时间发射的场景必须在程序中加入时间累积项用时间步进的方式求解温度场而不是直接用稳态公式。我在程序中加了一个可选的时间演化模式使用显式欧拉法逐步更新温度场可以按需开启。7.4 最后说些真心话做了这么多组仿真之后我个人最深的体会是热晕仿真的难点不在“算出一个相位屏”而在于“验证这个相位屏物理上合理”。我在初步完成程序后对照了公开发表的实验数据发现单纯用解析近似公式计算的相位屏往往高估了热晕效应——原因在于真实大气中还存在湍流混合、气流起伏等附加机制会部分缓解热聚焦效应。所以如果你要用这套程序做工程论证千万不要拿单一工况的仿真结果作为最终结论。至少要跑一个参数矩阵看看系统在好天气、坏天气、不同风速下的表现包络再结合统计学方法给出概率性的评估。仿真程序的价值在于帮我们理解趋势和机理而不在于提供一个精确到小数点的绝对数值。这些年我用它做过方案对比、系统参数优化、甚至写了小半篇论文的分析部分希望这些实现细节也能让你少走几步弯路。本文还有配套的精品资源点击获取