强激光大气热晕效应相位屏建模与MATLAB仿真实现-北京尧图网络

简介基于MATLAB的热晕相位屏仿真程序包面向光学系统与大气传输研究场景帮助科研人员和工程师模拟高能激光、望远镜等系统中由大气温度不均匀引发的波前畸变。压缩包共13个文件包括7个.m脚本、4个.bmp图像、1个.asv临时备份和1个.txt说明文档整体仅541KB。7个.m脚本覆盖了从相位屏随机生成、傅里叶变换到光学传递函数(OTF)与点扩散函数(PSF)计算再到反变换成像分析的完整链路.bmp文件直观展示不同环节的仿真输出.txt说明文档便于快速理解程序结构与参数含义。目前已有254人浏览学习。通过调整大气条件、温度梯度等输入参数可定量分析不同热晕场景下的成像退化效果进而为激光传输优化、望远镜设计及自适应补偿技术开发提供参考是一份适合科研入门和工程验证的实用光学仿真工具。1. 热晕相位屏仿真到底在解决什么问题强激光在大气里的“自毁”效应做强激光大气传输的人基本都撞见过这个现象发射功率越往上加靶面上的光斑反而不变亮而是变大、变形、重心飘移甚至在风速方向上被“吹”出一个不对称的尾巴。这不是光学平台没调平也不是自适应光学校得不准而是高功率光束把传播路径上的空气加热了空气折射率随之改变反过来又把光束自己偏折、发散掉——这个效应叫热晕。基于MATLAB的热晕相位屏仿真程序就是把“空气被加热后怎么改变相位”这一环节单独建模用一张或一组相位屏叠到自由空间传播的波前上在数值上复现热晕造成的散焦、偏转和高阶畸变。适合做激光大气传输预判的人、自适应光学方案评估的工程师以及刚接触相位屏法想先跑通一个最小模型的学生。本文不聊复杂的三维流体求解只讲工程上最常用、最容易验证的对流主导相位屏模型。2. 从光强到相位畸变热晕相位屏的物理模型与 MATLAB 生成函数2.1 对流主导假设为什么热晕相位是光强沿风向的积分热晕的物理链条不复杂光束被大气吸收单位体积吸收功率是 α·I其中 α 是吸收系数I 是光强吸收的功率变成热量空气温度上升温度变化引起密度变化进而引起折射率变化折射率变化量近似为 Δn (dn/dT)·ΔT波前经过这段空气时积累了额外相位 φ (2π/λ)·∫Δn dz。问题在于温度场怎么算。严格说需要解含对流的传热方程但在工程上大多数外场场景下横向风速远大于光束自身的扩散速度对流是热量输运的主导机制。这个时候当地温度增量正比于“从风上游吹过来的一路上积累的热量”也就是光强沿风向的积分。一维近似下设风沿 x 方向则ΔT(x,y) ≈ [α / (ρ·cp·v)] · ∫_{-∞}^{x} I(x,y) dxρ 是空气密度cp 是比热容v 是横向风速。这个近似成立的前提是风速导致的空气替换时间远小于导热扩散时间。一般风速 2 m/s 以上、光束尺寸在厘米量级时对流时间约几十毫秒导热时间在秒量级对流主导的假设是合理的。风速趋近于 0 时这个公式会发散需要用含扩散项的模型或者把多帧时间平均后的光强代入。把上面两段合起来就得到一个可以直接写进 MATLAB 的表达式φ(x,y) (2π/λ) · (dn/dT) · [α / (ρ·cp·v)] · ∫_{-∞}^{x} I(x,y) dx · dz_seg注意 dn/dT 在常温大气里约为 -1e-6 K⁻¹是负数。所以加热区域的相位增量整体为负意味着光程缩短、波前局部“提前”等效于一个负透镜这就是热晕把光束发散的来源。这个负号在代码里很容易弄丢丢了以后符号就反了焦斑会偏向下游和实测相反后面避坑章还会强调。2.2 生成函数用 cumsum 一次算出整面相位屏按上面的公式生成一张热晕相位屏只需要一个沿着风向的累积求和。MATLAB 里最稳的写法是 cumsum不要用双重 for 循环去逐点积分网格到 512×512 之后循环慢得让人怀疑人生。function phi thermalBloomingPhaseScreen(I, dx, param) % 生成一张热晕相位屏 % I : 当前光束光强分布 W/m^2, 尺寸 NxN % dx : 网格间距 m % param : 结构体包含 lambda, dn_dT, alpha, rho, cp, v, dz_seg % % 风向约定: x 方向为风向网格最左边是迎风边缘 N size(I, 1); % 沿风向的上游积分: cumsum 默认从左到右刚好是从迎风边到下游边 intI cumsum(I, 2) * dx; % 相位系数: C 本身是负值, 因为 dn_dT 0 C (2 * pi / param.lambda) * param.dn_dT ... * param.alpha / (param.rho * param.cp * param.v) * param.dz_seg; phi C * intI; % 去掉整体常数倾斜, 只保留波前相对起伏 phi phi - mean(phi(:)); end逻辑说明cumsum 是向量化累积求和intI 的每个像素等于该像素左侧所有光强按 dx 累加的结果这正好对应“风把上游热量带到我这里”的物理图像。mean 减掉的是波前整体常数项它只影响常数相位不影响焦斑形态和偏移。C 变量把波长、材料参数、段长全部收进去后续调参数只用改结构体不用动主循环。参数说明alpha 单位是 1/m典型值在 1e-5 到 1e-4 之间rho 取 1.2 kg/m³cp 取 1004 J/(kg·K)这两个是常温常压空气的工程值v 是横向风速单位 m/s最低不要低于 1 m/s否则对流近似就不成立了。dz_seg 是这段相位屏代表的传播距离单位 m与分段数配套下一章会讲。这个函数输出的 phi 是 N×N 的相位阵列单位是弧度量级直接在频域或空间域乘到光场上都行。2.3 光强从哪里来初始光束采样与网格布局相位屏生成需要 I而 I 来自上一段传播结束后的光场。初始光场一般按高斯光束给但要注意采样网格和实验参数匹配。网格宽度 D 和采样数 N 直接影响相位屏的空间分辨率网格左边沿要对应迎风侧否则风向积分就从反方向开始算了。% 初始高斯光束采样 N 512; D 0.4; % 网格总宽度 0.4 m dx D / N; x (-N/2 : N/2-1) * dx; [X, Y] meshgrid(x, x); w0 0.05; % 束腰半径 5 cm lambda 10.6e-6; % CO2 激光波长 P0 5000; % 峰值功率 5 kW A0 sqrt(2 * P0 / (pi * w0^2)); % 中心振幅 E A0 * exp(-(X.^2 Y.^2) / w0^2);逻辑说明meshgrid 生成的 X、Y 是物理坐标包含正负方向这比用数组下标直接算光场要直观得多。振幅平方后就是光强 I。这里用的是高斯光束如果实际是平顶光把振幅包络换成 super-Gaussian 的 exp(-((r/w0)^10)) 即可其余流程不变。参数说明D 要留足余量光束传播到 3 km 后衍射扩展半径至少是初始半径的几倍D 只给 4 倍束腰宽度是远远不够的。一般我会把 D 取到初始束腰的 8 倍以上再配合角谱传播窗口否则边缘衍射伪影会和热晕效应混在一起。N512 是平衡速度和精度的一个常用起点要是相位屏出现明显锯齿再往上加。3. 把相位屏装进分段传播主程序结构、角谱传播与一张参数表3.1 分布效应不能一张屏贴在末端分段传播的物理原因热晕是沿光束路径连续发生的不是集中在某个位置的一次性扰动。每一段路径上的光强都不相同因此每一段的相位屏贡献也不同。常见错误是把整段路径的相位一次性算好然后放在发射端这样会把高功率条件下的像差分布完全算错。分段传播的做法是把 L 分成 M 段每段长度 dz L/M光场先自由传播 dz在段末端采样光强由当前光强生成这一段的热晕相位屏乘到光场上再继续下一段。这样每一段相位屏都使用“该段实际光强”作为加热源等效于把分布效应离散成一串薄透镜。需要叠加大气湍流时也是同样思路把 Kolmogorov 相位屏插在每一段上和热晕相位屏相乘或相加二者互不干扰。分段数的选取直接决定仿真的可信度工程上比较实用的判据是单段产生的峰值相位不超过 λ/10。按这个标准初始功率越高、吸收越强、风速越慢M 就要越大。后面避坑章会给具体的对比方法。3.2 主循环代码传播、采样光强、叠加热晕相位屏主循环里最核心的传播函数是角谱法。它比傅里叶变换法的实现更稳因为只需要一步 FFT 和一步 IFFT不会出现 z 方向步长相关的相位混叠。freePropagate 函数实现如下。function E freePropagate(E, dz, lambda, N, D) % 角谱法自由空间传播 dx D / N; k 2 * pi / lambda; fx (-N/2 : N/2-1) / D; [FX, FY] meshgrid(fx, fx); F fftshift(fft2(fftshift(E))) * dx^2; H exp(1i * pi * lambda * dz * (FX.^2 FY.^2)); E ifftshift(ifft2(ifftshift(F .* H))) / dx^2; end逻辑说明角谱法的核心是把光场变换到频率域乘上一个与传播距离相关的传递函数 H再逆变换回空间域。H 使用的是近轴二次近似对束腰半径在厘米级、传播距离在千米级的场景精度足够。这里 FFT 前后都用了 fftshift为的是让频率零点和空间零点对齐到数组中心很多新手的相位屏错位问题都是漏了这一步。主程序的完整结构可以写成下面这样% thermal_blooming_main.m 分段热晕相位屏仿真主程序 clear; close all; % 物理参数 lambda 10.6e-6; P0 5000; w0 0.05; L0 3000; % 总传播距离 alpha 5e-5; % 吸收系数 1/m v 3; % 横向风速 m/s rho 1.2; cp 1004; dn_dT -0.9e-6; % 数值参数 N 512; D 0.4; dx D / N; M 12; % 分段数 dz L0 / M; % 初始光场见 2.3 节 E initialGaussianBeam(N, D, w0, P0, lambda); % 分段传播 for ii 1 : M E freePropagate(E, dz, lambda, N, D); I abs(E).^2; param struct(lambda, lambda, dn_dT, dn_dT, ... alpha, alpha, rho, rho, cp, cp, v, v, dz_seg, dz); phi thermalBloomingPhaseScreen(I, dx, param); E E .* exp(1i * phi); end % 输出焦平面光强 I_target abs(E).^2; imagesc(I_target); axis image;逻辑说明每个循环周期内光场先走完一段自由空间再按当前光强生成热晕相位屏乘到光场上。由于相位屏的值代表的是这一段距离内的累积效应放在传播之后是物理上正确的顺序。exp(1i·phi) 直接把相位屏转成相位因子这在标量衍射理论里是标准操作。参数说明M12 对 5 kW、3 km、风速 3 m/s 的情况属于保守偏大的设置跑一次大约十秒量级可以接受。调到 M4 也能出结果但热晕高阶畸变会被低估焦斑底部形态明显不同。dz 一定要保持在自由传播能正常工作的范围内角谱法在 dz 很小时数值上不会有问题这点比直接分步傅里叶法要省心。3.3 一张参数表把物理量换算成程序变量新手最容易犯的错不是算法不对而是参数单位混乱。下面这张表是我自己写这类仿真时固定放在脚本头部的对照表每换一台机器、换一组实验条件先改这张表再动代码。参数符号典型值单位备注波长lambda10.6e-6m10.6 μm CO2 激光热晕效应明显发射功率P0100010000WP0 越高分段数需求越大束腰半径w00.030.10m决定初始光强峰值吸收系数alpha1e-51e-41/m10.6 μm 在湿空气中偏大空气密度rho1.2kg/m³常温常压近似比热容cp1004J/(kg·K)干空气折射率温变系数dn_dT-0.9e-61/K恒为负值横向风速v210m/s低于 1 m/s 时模型失效传播距离L010005000m与分段数一起决定 dz网格宽度D0.30.8m至少覆盖衍射扩展后的光束采样数N2561024-512 是常用折中分段数M620-按单段相位 λ/10 校验表的最后一行分段数最容易被忽略。我一般会在程序里加一个探针变量把每一段生成的 phi 最大值打出来如果某段的最大相位超过 0.6 rad就自动提示要增加 M。这个做法相当于给仿真装了一个自检仪表比事后看焦斑畸形再去猜原因高效得多。4. 四个决定结果对错的关键参数吸收系数、风速、网格与分段数4.1 吸收系数 α量纲错一位光斑直接“消失”α 是热晕相位屏里最敏感的参数有个数量级的偏差相位屏的强度就会差十倍。常见错误是把手册里给的 cm⁻¹ 单位直接当成 m⁻¹ 用。比如某本参考资料写“吸收系数 5×10⁻⁷ cm⁻¹”换算成 m⁻¹ 要乘以 100得到 5×10⁻⁵ 1/m。漏掉这个换算热晕效应会小到完全看不见仿真结果和裸光束衍射几乎没有区别。怎么判断 α 给得合不合理一个快速的量级检查是估算峰值相位。取 w00.05 mP05000 Wv3 m/sα5e-5 1/mL03000 m沿风向上的积分长度大约是光束宽度量级算下来峰值相位在零点几个弧度到 1 弧度之间这个量级正好能产生可观察的焦斑畸变。如果算出来的峰值相位是 1e-3 弧度以下大概率是 α 的单位错了如果超过 3 弧度就要怀疑是程序里 I 的单位不对很可能把峰值功率当成了平均功率。4.2 风速 v 和风向相位屏符号与焦斑偏移方向的互检风速的绝对值决定相位屏强度方向决定积分方向。前面代码里把风定义成 x 方向网格最左端是迎风边缘相位屏系数 C 为负下游光程更短。这个设置下焦斑重心应该往 -x 方向偏移也就是朝迎风侧偏移这是热晕的经典特征。我每次跑完仿真第一件事不是看 Strehl 比而是看焦斑重心往哪个方向偏。偏错了不用查别的先去查相位屏符号或者查 cumsum 方向。有人会问既然偏转只是整体倾斜对 Strehl 比没影响为什么还要关注因为在有自适应光学的系统里整体倾斜是最好校正的一阶项真正的麻烦是散焦和高阶像差。如果整体偏移方向都错了整个仿真对自适应光学系统的评估就全部失去了意义。4.3 网格宽度与采样数窗口选多大才不卷绕网格宽度 D 不能拍脑袋定。初始束腰 5 cm传播 3 km 后衍射极限焦斑半径大约在 λL/πw0 ≈ 0.2 m 量级。焦面上的光斑尺寸远比发射端大如果 D 还是取 0.2 m焦斑会超出网格窗口边缘信息被截断imae 里能看到光强硬生生被窗口切掉那后面的重心和 Strehl 比全是错的。经验规则是 D 至少取到初始束腰的 8 倍以上并且要大于自由传播后的衍射扩展直径。N512 时0.4 m 宽的窗口配 5 cm 束腰是比较均衡的组合。如果窗口里出现从一侧绕过传回另一侧的干涉条纹说明 D 不够不是 N 不够盲目加 N 只会让计算更慢而问题照旧。4.4 分段数 M单屏与多屏对比的判定方法分段数 M 是热晕仿真里最像“玄学”的参数不同论文里给出的经验值差距很大。我的做法是直接用判据卡每一段的峰值相位不超过 λ/10。峰值相位可以用上一段传播后的光强和公式算出来写几行 MATLAB 代码就能在循环里检查。如果某段算出的 φ_max 大于 0.6 rad就增加 M 重新跑。判断 M 够不够的另一个方法是做对比实验分别用 M4、M8、M16 跑同一组参数比较目标面焦斑的 Strehl 比。如果 M4 到 M8 差距很大M8 到 M16 差距很小说明 M16 附近已经收敛。如果 M4 和 M16 几乎没差别那说明当前功率和吸收条件下热晕本来就弱M 取 4 也行不需要为了心理安慰把计算量翻倍。5. 热晕相位屏仿真的五个常见坑现象、原因与排查顺序5.1 焦斑偏移方向反了风向与积分方向冲突现象仿真焦斑的偏移方向和文献里热晕的经典结论正好相反往顺风方向偏。原因相位屏符号或 cumsum 方向有一个反了。常见的具体错法是先在代码里把网格左边缘当成了下游边或者把系数 C 前的负号漏掉导致下游光程变长等效成正透镜光束被朝顺风侧偏折。解决把风向约定写死在参数注释里并在生成函数入口处固定网格左边为迎风边。跑完第一步传播后打印相位屏在 x 方向的斜率正常情况斜率应为负值代表下游光程更短。如果斜率为正立刻能定位。5.2 网格边缘出现周期性亮纹卷积卷绕伪影现象焦面上除了主光斑边缘出现一排间距均匀的细干涉条纹像衍射光栅的谱线。原因角谱法的 FFT 隐含着周期性边界条件光场超出网格宽度后会从另一侧卷回。热晕相位屏让光束整体偏转后一部分光强被推到网格边缘卷绕回来以后就和主光斑相干叠加形成周期条纹。解决把 D 增大到光束扩展后的 2 倍以上同时给相位屏加一个软边窗函数比如在相位屏外围乘以 cos² 锥形衰减范围内保持原值。这个操作不影响有效光斑只是把网格边缘的残余能量压掉极大地改善边缘伪影。5.3 相位屏算出来数值巨大相位超过了 π 的处理现象把 phi 数组打出来最大值超过 3、4 弧度甚至到了几十弧度。原因exp(1i·phi) 对相位本身没有幅度限制但过大的相位意味着单段传播内热晕畸变已经强到超出薄屏近似的适用条件。相位屏代表的是薄扰动一次承担几个波长的光程差数值上虽然不出错物理上已经不可信。解决把 M 加大让每一段的峰值相位回落到 λ/10 以内。还有一种情况是 I 的单位给错会导致相位虚高比如光强用了 mW/cm² 忘了换算成 W/m²这时候先查单位表再提 M。5.4 分段数太少导致畸变被低估怎么判断 M 够不够现象M4 和 M8 的结果差别很大但很多人只跑了一次 M4 就开始分析最后得出的“热晕不明显”结论是错的。原因热晕的非线性来自相位屏和光强的耦合。分段太少时光强没有在段间被相位屏充分修改后续段使用的光强过于接近自由传播结果等效于低估了热晕的反馈作用。解决把 M 作为扫描参数跑一组对比。如果 Strehl 比随 M 单调变化且超过 10%说明当前 M 不够。我自己的习惯是至少跑 M8 和 M16 两组结果一致才继续往下做参数分析。5.5 单位不一致从 α 到 ρcp 的换算清单现象仿真结果有时明显偏弱有时强到离谱排查代码结构无误最后发现是 α 或 I 的单位错了。原因这个模型里的单位链条很长α 有 cm⁻¹ 和 m⁻¹光强有 W/cm² 和 W/m²D 有 cm 和 m任何一处不一致都会让系数 C 差出几个数量级。解决把参数表中的单位写成注释贴在自己的脚本头部。每次从论文或手册抄参数先换算成 SI 单位制再填入。我自己翻车次数最多的就是 α所以现在习惯把 α 的换算单独写一行注释查到的值若以 cm⁻¹ 为单位填表时乘以 100 转成 m⁻¹。6. 从静态屏到时间序列用 Strehl 比和光斑重心验证仿真6.1 两个验证指标Strehl 比与焦斑重心坐标静态相位屏跑完不能只看 imagesc 的图觉得“像”就完事要用数字指标卡住结果。两个最实用的指标是 Strehl 比和焦斑重心。Strehl 比直接反映热晕造成的能量集中度损失重心坐标用来核对偏移方向。% 焦面指标计算 I_far abs(E).^2; I_far I_far / sum(I_far(:)); % 归一化能量 % 参考: 无热晕的理想衍射焦斑 I_diff SR max(I_far(:)) / max(I_diff(:)); % Strehl 比, 通常 1 % 重心坐标(像素单位) [Ng, ~] size(I_far); [gx, gy] ndgrid(1:Ng, 1:Ng); cx sum(gx .* I_far) / sum(I_far(:)); % x方向重心 cy sum(gy .* I_far) / sum(I_far(:)); % y方向重心 fprintf(Strehl %.3f, 重心偏移 [%.2f, %.2f] px\n, SR, cx-Ng/2, cy-Ng/2);逻辑说明Strehl 比的定义是把实际焦斑峰值强度除以理想衍射焦斑峰值强度理想值 1实际小于 1。重心偏移的单位是像素要先换算成物理坐标再和风向核对。偏移方向必须是迎风侧这一点无论功率高低都应该成立。6.2 把静态屏变成时间序列风场平移法的取舍静态单帧仿真只能回答“稳态热晕长什么样”回答不了时间演化问题。如果需要看焦斑随时间抖动常见做法是风场平移法把生成的静态相位屏按风向按时间步长平移若干个像素模拟空气团被风吹过光束路径的过程。时间步长取 tc 2w0/v也就是风把整个光束截面吹过去一次所需的时间。这个方法的优点是简单缺点是它强行假设温度场被风整体平移忽略了导热扩散对温度分布的平滑作用。我一般只用它来做初步的动态评估比如看焦斑重心震荡幅度和频率范围。如果项目对动态精度要求高就要上含时间项的分步传热求解那已经不是相位屏程序能覆盖的范围了。从这一段一路写下来我自己最大的感受是热晕相位屏仿真的精度不是靠某一招提升的而是靠分段数、网格宽度和符号约定这些看似琐碎的细节一起堆出来的。前期花十分钟做参数自检比后期对着离谱的焦斑图猜原因省一整天时间。希望这份流程能帮你在自己的项目里少踩几个坑。本文还有配套的精品资源点击获取

强激光大气热晕效应相位屏建模与MATLAB仿真实现

相关推荐

航拍滑坡数据集4315张VOC+YOLO双格式

程序员数学知识地图：概率统计线代离散图论速查与Python验证

拆解Amical的whisper.cpp封装：如何构建带Metal/CUDA/CPU自动回退的C++原生模块

更多资讯

CSS 闪烁效果（Blink）实战：用 You-need-to-know-css 示例彻底搞懂 animation-direction 四个取值

AWS SDK for Go v2 DynamoDB 服务客户端演进全解：从版本历史到 Cortex 落地实践

2019大数据机器学习教学实践：从答案文档反推工程闭环

高校运动会管理系统数据库设计：从ER图到MySQL完整实现

PID温度控制实战：从参数整定到积分饱和的完整避坑指南

国产芯片如何通过工控机实现工业级实时性与环境适配

今日更新

Susi AI API实战指南：susi_alexa_skill如何用Node.js调用chat.json获取智能回答

换新电脑了？KeyStats 恢复码数据找回完全指南，端到端加密统计一键重建

重建二叉树：由前序与中序遍历序列还原二叉树的原理与 Java 实现（YCBlogs 剑指 Offer 精讲）

本周精选

交通拥堵预测毕设：基于Python+LSTM的车辆流量预测系统全流程解析

MATLAB/Simulink搭建10机39节点电力系统暂态稳定仿真实战

无人机航拍三维重建全流程：从SfM到网格生成的避坑指南

本月热门

我发现了一个新思路：用 Remotion + Claude Code 像写代码一样自动化生成短视频

Windows下 Codex 中 Chrome 和 Computer Use 插件不可用问题排查及解决参考方式：TaoToken 统一 Key 配置与验证

2026 大模型集体涨价：用 Python 做企业 Token 成本测算与选型避坑（附配置）

读完文章，想聊聊您的网站？