简介:本资源是一套基于MATLAB/Simulink开发的热晕相位屏仿真程序,面向光学工程、激光大气传输及自适应光学领域的科研人员与高年级本科生/研究生,用于定量研究高能激光在大气中传播时由热晕效应引起的波前畸变。程序支持灵活设置激光功率、大气参数、光束特性等条件,可生成对应热晕相位屏并完成频域/空域转换、相位重建与可视化分析。压缩包共12个文件(7个核心M脚本实现相位计算、FFT/IFFT处理、温度场建模与屏生成;4幅BMP图像为典型相位屏或中间结果示例;1个ASV备份文件),总容量541KB,结构紧凑、模块清晰,便于理解热晕物理机制与算法实现逻辑。目前已有219人学习下载,用户可直接运行主程序xuanhuan.m,结合reynjisuanzz.m、myfft2.m等关键函数掌握热晕建模全流程,快速复现实验结果并开展参数敏感性分析。
1. 项目概述与整体设计思路
1.1 热晕相位屏仿真到底解决了什么问题
做激光大气传输的人应该都绕不开一个现象——热晕。高能激光在空气中传输时,大气分子和气溶胶会吸收一部分光能,这部分能量转化为热量,导致光束路径上的空气被局部加热。空气被加热后密度下降,折射率跟着变化,而折射率的变化又会反过来改变光束的相位分布,造成光束扩展、畸变、弯曲,严重时还会出现光束分裂。这个非线性过程就是热晕效应(thermal blooming)。
单靠解析方法来分析热晕几乎不可能,因为这是一个光束与介质互相耦合的非线性过程。工程上最常用的做法就是数值仿真,而相位屏(phase screen)方法是其中效率最高、最灵活的手段之一。相位屏的本质是把连续介质的相位畸变等效为一系列离散的薄屏,光束每经过一个屏就累积一次相位扰动,从而在大大降低计算量的同时保留主要物理特征。
这个Matlab/Simulink程序做的就是这件事,核心功能是在给定激光功率、波长、光束半径、风速、大气吸收系数等条件下,生成对应的热晕相位屏。程序支持不同参数工况的切换,可以用一组仿真代码批量研究不同条件下的热晕强度变化。对于做激光传输评估、自适应光学校正算法验证、光通信链路仿真的工程师和科研人员来说,这类工具能省下大量重复造轮子的时间。
1.2 为什么选Matlab/Simulink做热晕仿真
要是把热晕仿真放到大型C++工程里做,光处理数组操作和可视化就得写不少代码,调试周期往往不短。Matlab在矩阵运算和快速原型验证上有天然优势,二维复振幅数组就是天然的网格数据,傅里叶变换调用fft2一行搞定,相位屏的生成、叠加、传播都能用矩阵化编程实现,代码量能压缩到C++版本的十分之一。
Simulink在这个项目中的定位不是替代Matlab的数值计算,而是负责仿真流程管理和多工况调度。热晕仿真往往需要跑不同风速、不同功率、不同湍流强度的组合,用Simulink搭一个模块化框架,可以把相位屏生成、传输计算、数据记录封装成独立模块,切换工况时不用改脚本,直接改模块参数就能批量跑。这种组合方式在实际项目中很实用,尤其是后续要接入控制系统仿真或者半实物仿真时,Simulink的接口优势就更明显了。
1.3 程序整体框架
程序围绕“参数输入—相位屏生成—光束传输计算—结果可视化”这条主线设计:
% 主入口:run_thermal_blooming.m % 1. 设置物理参数 % 2. 生成热晕相位屏 % 3. 模拟光束经过相位屏后的远场光斑 % 4. 输出对比图和评价指标Simulink模型则负责把主循环包装成模块化架构,用Parameter Configuration模块统一管理参数。整体上Matlab负责核心计算,Simulink负责流程调度,两者通过Interpreted MATLAB Function或MATLAB Function模块进行数据交互。
2. 热晕相位屏的核心物理模型与数学表达
2.1 热晕的基本物理图景
要生成一个靠谱的热晕相位屏,首先得把热晕的物理过程用数学表达出来。热晕效应可以粗略地分成两步:第一步是激光加热空气引起的折射率变化,第二步是折射率变化对光束相位的影响。
从流体力学角度看,激光加热导致的空气密度变化满足等压近似下的状态方程:
[ \frac{\Delta n}{n_0} = -(n_0 - 1)\frac{\Delta T}{T_0} ]
其中n₀是未扰动空气折射率,T₀是环境温度,ΔT是温升。而温升由吸收激光能量引起,在连续波激光且风速横向吹过光束的情况下,温度分布会呈现出上游冷、下游热的不对称形态。
基于这个物理解释,工程上最常用的热晕相位屏近似模型把热晕相位分解成离焦项、像散项和三次畸变项的组合。离焦项对应于热晕引起的等效负透镜效应,光束通过后发散;像散项对应于横向风导致的不对称性;三次畸变项则描述了更高阶的波前畸变。这三项叠加后,用Zernike多项式或者直接计算相位分布,就能构造出热晕相位屏。
2.2 Bradley-Hermann近似模型
实际项目中用得最多的热晕解析模型是Bradley-Hermann模型,它给出了热晕导致的远场峰值光强退化因子。对于单层相位屏近似,热晕相位可以写成无量纲参数的形式:
[ \phi_{th}(x, y) = -N_D \cdot \frac{\sqrt{\pi}}{8} \cdot \frac{\int_0^z (x - vt')^2 , dt'}{...} ]
这个表达式看着很复杂,但物理上有明确的含义。N_D是热畸变数,它综合了激光功率P₀、吸收系数α、风速v、光束半径a和传输距离z等参数,是一个无量纲数,用来衡量热晕效应的强度。N_D越大,热晕越强。
在程序中,热晕相位屏的最终表达式可以简化为:
[ \phi_{th}(x,y) = -N_D \cdot \text{norm_factor} \cdot \left[ \frac{x}{a} + \frac{1}{2}\left(\frac{x^2 + y^2}{a^2}\right) + \frac{1}{3}\left(\frac{x^3 + 3xy^2}{a^3}\right) \right] ]
其中x方向是横向风方向。这就是热晕相位屏的核心生成公式,程序中只需要确定N_D的数值和光束半径a,就能生成对应的相位屏。
2.3 随机湍流相位屏与热晕相位屏的叠加
实际大气中,热晕和湍流总是同时存在的。程序里的热晕相位屏,准确的说是“热晕+湍流”总相位屏,或者提供选项让用户选择只生成热晕分量、只生成湍流分量、还是两者叠加。
湍流相位屏的生成使用功率谱反演法。基本原理是:大气折射率起伏的功率谱密度符合Kolmogorov或von Kármán谱,通过对随机频谱进行滤波,再做逆傅里叶变换,就可以得到空间域的随机相位分布。核心代码如下:
function phase = generate_turbulence_phase(N, delta, L0, l0, Cn2) % 生成湍流相位屏 % N: 网格数 % delta: 网格间距(m) % L0: 外尺度(m) % l0: 内尺度(m) % Cn2: 折射率结构常数(m^-2/3) k = 2*pi*(-N/2:N/2-1)/(N*delta); [Kx, Ky] = meshgrid(k, k); K = sqrt(Kx.^2 + Ky.^2); K(N/2+1, N/2+1) = 1e-10; % 避免零频 % von Karman谱 phi_phi = 0.49*r0^(-5/3) ./ (K.^2 + 1/L0^2).^(11/6) ... .* exp(-K.^2 ./ (2*pi/l0)^2); % 随机复高斯滤波 random_phase = (randn(N) + 1i*randn(N)) / sqrt(2); phase = real(ifft2(sqrt(phi_phi) .* fftshift(random_phase) * N^2)); end热晕相位屏生成的主函数则是基于Bradley-Hermann近似模型:
function phi_th = generate_thermal_phase_screen(N, delta, x_wind, N_D, a) % 生成热晕相位屏 % N_D: 热畸变数 % a: 光束半径(m) % x_wind: 横向风速方向 [x, y] = meshgrid((-N/2:N/2-1)*delta, (-N/2:N/2-1)*delta); r2 = (x.^2 + y.^2) / a^2; x_norm = x / a; phi_th = -N_D * (x_norm + 0.5*r2 + (1/3)*(x_norm.^3 + 3*x_norm.*(y/a).^2)); end需要提醒的是,这里的N_D是热畸变数,通常定义为一个正数,前面这个负号表示热晕产生的是负透镜效应——光束中心区域相位滞后,等效于发散透镜。相位符号的判断直接影响后续远场计算结果,很多人在这里容易搞错符号导致光斑聚焦反而增强的假象。
3. Simulink环境下的模块化集成
3.1 为什么要在Simulink里搭框架
如果用Matlab脚本直接跑,对于单次仿真确实足够了。但实际工程中往往需要扫描几十组参数,比如考察不同风速下热晕强度变化,或者比较不同功率下的光斑畸变程度。每次手动修改脚本参数再运行,效率低还容易出错。
Simulink在这时候的价值就体现出来了。通过把热晕仿真流程封装成Simulink模型,可以用MATLAB Function模块嵌入相位屏生成的算法,用Constant模块设置参数,用To Workspace模块记录结果,再用Simulink的batch simulation功能一次性跑完所有工况组合。还可以把Simulink模型打包成引用模型(Model Reference),嵌入到更大的激光系统仿真框架中。
3.2 Simulink模型结构
模型的顶层结构可以这样搭:
- 参数输入层:用Simulink.Parameter对象定义Nd_base、wind_speed、Cn2、wavelength等参数
- 相位屏生成层:一个MATLAB Function模块,输入归一化坐标和物理参数,输出相位屏矩阵
- 光束传播层:计算光束通过相位屏后的远场分布
- 数据输出层:To Workspace输出不同位置的光斑和相位数据
MATLAB Function模块里的代码可以直接复用脚本函数,接口用coder.extrinsic声明外部函数,确保代码生成兼容性。如果后续要做嵌入式部署,还可以用Simulink Coder把光束传播部分生成C代码。
3.3 多工况批处理实现
在Simulink中做多工况批处理,有几种方式。最简单的是用仿真输入对象(Simulink.SimulationInput),在循环里给不同工况赋值:
for i = 1:length(Nd_array) simInput(i) = Simulink.SimulationInput('thermal_model'); simInput(i) = simInput(i).setVariable('N_D', Nd_array(i)); simInput(i) = simInput(i).setVariable('wind_speed', wind_array(i)); end out = sim(simInput, 'ShowProgress', 'on', 'UseParallel', 1);用UseParallel参数配合并行计算工具箱,能把批处理时间缩短到原来的四分之一。我实测下来,在8核机器上跑15组工况,从28分钟压缩到7分钟左右,效率提升非常明显。
另一个方案是用Simulink的batch simulation manager(在Simulink工具条里直接打开),把参数组合填到表格里一键运行,跑完还能自动生成对比报告。这个方式更适合不太熟悉脚本操作的人。
4. 不同条件下的仿真结果分析
4.1 参数设置与基准工况
仿真程序需要定义一组基准物理参数,后续所有对比都基于这组参数进行单变量调整。典型的基准工况如下:
| 参数 | 符号 | 数值 | 单位 |
|---|---|---|---|
| 激光波长 | λ | 1.064e-6 | m |
| 激光功率 | P | 100 | kW |
| 光束半径 | a | 0.1 | m |
| 大气吸收系数 | α | 0.05 | km⁻¹ |
| 横向风速 | v | 5 | m/s |
| 传输距离 | z | 5 | km |
| 折射率结构常数 | Cn² | 1e-15 | m⁻²/³ |
| 网格数 | N | 256 | - |
| 网格间距 | δ | 0.002 | m |
根据这些参数先计算热畸变数N_D,公式为:
[ N_D = \frac{2\sqrt{2}(-dn/dT)\alpha P z^2}{\pi a^2 n_0 \rho c_p v a} ]
具体数值计算在程序中自动完成。对于基准工况,N_D大约在几到十几之间,属于中等强度热晕。
4.2 不同功率下的热晕相位屏
把激光功率从20 kW逐步提高到200 kW,观察热晕相位屏的变化。功率为20 kW时,热畸变数N_D小,相位屏的离焦分量较小,相位起伏平缓,光斑畸变不明显。功率达到200 kW时,N_D增大到原来的10倍左右,相位屏中心相位差显著增大,离焦项占主导,等效负透镜效应明显增强,远场光斑扩展显著。
程序输出的相位屏对比图能清晰看出,低功率时相位分布趋于平缓,高功率时相位屏呈现明显的碗状凹陷,颜色从中心到边缘的梯度变化剧烈。这就是热晕“自散焦”的数值体现。
4.3 横向风速的影响
风速是影响热晕的重要因素。风速越大,热量被带走越快,空气温升越小,热晕越弱。程序在风速为1 m/s、5 m/s、20 m/s三种条件下生成的热晕相位屏存在明显差异。
1 m/s风速时,热晕相位屏的像散项非常突出,相位分布沿风向呈不对称形态,光束畸变严重;20 m/s风速时,由于对流冷却作用增强,相位屏整体起伏显著减小,近似接近纯湍流情形。这组对比特别能解释为什么实际激光系统会关注风速风向——条件允许时会选择迎风发射或者大风天作业。
4.4 湍流与热晕叠加效应
程序的一个关键特性是支持湍流热晕叠加仿真。只加湍流时,相位屏呈现高频随机起伏;只加热晕时,相位屏呈现低频大尺度畸变;两者叠加后,高频和低频分量同时存在,远场光斑同时出现小尺度破碎和大尺度扩展。
这一结果对自适应光学校正系统的研究特别有价值。校正系统通常能有效校正低频像差,但对高频湍流分量的校正能力受限于子孔径数和波前传感器采样率。用这个叠加相位屏就能定量评估校正系统在不同热晕强度下的性能边界。
5. 程序设计中的几个关键细节
5.1 网格分辨率与采样间隔的选择
相位屏仿真的精度严重依赖网格参数。采样间隔δ决定了相位屏的最高空间频率,网格数N决定了频率分辨率。热晕相位屏主要是低频分量,δ的选择可以相对宽松;但湍流相位屏包含高频分量,δ必须足够小才能捕捉内尺度以内的湍流涡旋。
经验法则是:网格间距至少小于内尺度l₀的一半,即δ < l₀/2;仿真区域边长N·δ至少要大于光束直径的两到三倍。如果用256×256网格、δ=2mm,那么仿真区域边长约为0.5m,是光束直径的2.5倍,满足要求。如果出现明显的光束能量“溢出”计算区域的问题,优先增加网格数而不是增大间距。
5.2 相位屏生成后要做低通滤波吗
热晕相位屏本身是平滑的解析函数,不需要滤波。但湍流相位屏用功率谱反演法生成时,由于FFT的周期性假设,会出现低频分量不足的问题。解决办法是用次谐波(subharmonic)补偿低频段,或者直接用Zernike多项式叠加低频项。
在叠加热晕和湍流相位屏时,建议先分别生成两个屏,再进行相位加法。直接混合生成会导致高频和低频相互干扰。相位相加时要注意单位统一,两个相位屏的波长必须一致。
5.3 傅里叶变换传播的数值细节
光束传播计算中,程序使用了角谱法或者菲涅尔衍射积分的FFT实现。这里有一个常见坑:使用fftshift和ifftshift时容易搞混顺序。标准流程是:输入复振幅场先fftshift,再做fft2,然后乘以传输函数,再做ifft2,最后ifftshift回来。如果中间传递函数也做了fftshift,那会重复移位一次,结果出现倒像。我在调试程序时花了不少时间排查这个问题,后来干脆封装成函数,固定好每一步的shift操作。
5.4 Simulink仿真速度优化
Simulink模型跑一次仿真涉及Matlab引擎调用,如果直接在MATLAB Function模块里写循环,速度会很慢。优化手段有三个:
第一,把可并行的数组运算全部向量化,避免for循环。相位屏生成的核心计算全部用矩阵操作,没写一个for循环。
第二,用coder.extrinsic声明不需要代码生成的函数,比如绘图函数和文件读写函数,这样在Simulink中调用时不至于触发不必要的代码生成检查。
第三,如果做批处理,用‘UseParallel’开启并行池,配合Fast Restart模式,能大幅缩短循环仿真时间。不要用sim函数里的StopTime参数去做多次稳态仿真再拼接,那种方式反而更慢。
6. 工程应用中的扩展方向
做完基础的热晕相位屏仿真,程序还可以向多个方向扩展。
一是把单层相位屏推广为多层相位屏。对于长距离传输,热晕沿路径是不断累积的,用单层相位屏是粗略近似,更严谨的做法是把传输路径切成若干段,每段分别生成相位屏,光束在段与段之间做真空衍射传播。这个扩展在程序中预留了接口,只要把传输路径向量作为参数传入即可。
二是瞬态热晕仿真。目前程序默认是稳态热晕,即热平衡已经建立。但实际上激光开启后需要几十毫秒才能建立稳定的热晕相位分布,这个时间尺度对脉冲激光或者快速指向控制系统非常重要。瞬态仿真需要在原模型上引入时间变量,把相位屏的时间演化方程联立起来,计算量会比稳态大很多,但这类仿真能回答很多实际工程问题。
三是与控制仿真联合。热晕相位屏模型可以作为自适应光学仿真链路中的被控对象,用Simulink输出波前畸变,接入控制系统做闭环分析,评估校正算法的稳定裕度。这个方向在实际的项目里非常有价值,也是我当时做这个程序时最核心的后续用途。
关于相位屏的验证,我自己的习惯是拿程序跑出来的远场光斑与理论公式对比,尤其是低功率极限下,相位屏接近纯湍流,远场光斑要能恢复到近衍射极限的情况,否则就说明热晕项的处理有误。用这个方式能快速定位相位屏生成代码里的bug。
7. 实操中的几点体会
这套热晕相位屏仿真程序用Matlab和Simulink搭建,核心计算代码不长,但涉及到的物理模型和数值技巧相当密集。
我的体会是,相位屏方法真正难的不是生成相位屏本身,而是怎么把物理参数正确映射到相位屏系数上。N_D的计算结果对不对、符号取没取反、风向定义是否与坐标轴一致,任何一个环节出错,得到的仿真结果都会偏离物理真实,而且往往看起来还挺合理、不容易发现问题。所以做这类仿真,建议始终保留一组已发表文献中的验证工况,每次修改代码后先跑基准对比,确认结果一致再继续往下做。
另外,Matlab和Simulink联合仿真时,我比较推荐把物理计算尽量放在纯Matlab中完成,Simulink只做流程控制。原因很简单,Simulink的MATLAB Function模块虽然有代码生成能力,但调试体验远不如脚本环境,变量查看和断点设置都不够灵活。核心计算放脚本,流程管理放Simulink,两边都能发挥优势。
如果你也在做激光大气传输相关的仿真工作,希望能从这套程序的框架里得到一些参考。还有一个小技巧是,在生成热晕相位屏之前,先关闭湍流分量,单独看热晕的相位分布,确认物理形态合理了再打开湍流叠加,这样能少走不少弯路。
本文还有配套的精品资源,点击获取