简介本资源是一套基于MATLAB/Simulink开发的热晕相位屏数值仿真程序面向光学工程、高能激光传输、自适应光学等方向的科研人员与研究生用于定量研究大气或介质中热致折射率畸变对光束波前的影响。程序支持调节激光功率、光束尺寸、介质参数、传输距离等关键变量可生成对应条件下的动态相位屏数据为后续波前校正、光束质量评估提供基础输入。压缩包共12个文件7个核心M脚本实现相位计算、FFT/IFFT处理、热效应建模与可视化4幅BMP图像为典型相位屏结果示例1个ASV备份文件总大小541KB结构紧凑、模块清晰便于理解热晕物理模型与代码实现逻辑。目前已有219人学习下载用户可直接运行主函数xuanhuan.m复现完整仿真流程并通过对比不同参数下的相位屏图像如1.bmp、2.bmp等直观掌握热晕演化规律。 做强激光大气传输仿真的人第一眼看到“热晕”两个字往往觉得它是个比湍流更难缠的东西。我最初用MATLAB/Simulink搭热晕相位屏仿真程序时连相位屏长什么样都没概念等把第一个模型跑通才明白这玩意儿跟湍流屏完全是两码事。这篇文章就分享我怎么用MATLAB/Simulink实现不同条件下的热晕相位屏仿真包括物理模型、频域算法、Simulink封装和实际调试中踩过的坑。无论你是做激光通信、光电对抗还是高能激光系统设计只要需要在数值仿真里刻画大气热效应这套思路应该都能直接用。热晕效应最反直觉的一点是它不依赖随机涨落而是激光自己把传输路径上的空气“烤”出了一个大透镜再用这个大透镜去畸变自己。所以仿真热晕相位屏本质上不是在生成随机场而是在解一个流动与光强耦合的方程。理解了这一点后面所有代码和模型都能串起来。1. 热晕相位屏不是“随机屏”它是光给自己挖的坑1.1 激光加热空气到底发生了什么大气对激光并非完全透明尤其在高功率激光波段大气分子和气溶胶会吸收一部分光能。这部分能量变成热量让激光路径上的空气温度升高。空气的折射率随温度变化常温常压下折射率温度系数大约是 dn/dT ≈ -1e-6 K^-1也就是说温度每升高1K折射率大概下降百万分之一。百万分之一听起来很小但激光传输距离动辄几百米到几公里光程差的累积效应足够把一束好端端的高斯光撕得乱七八糟。关键是光强分布不均匀。高斯光束中心光强最强那么中心空气温升最大折射率下降也最多。于是光束中心区域的光速比边缘快相当于经过了一个负透镜。这个效应就是“热晕”或者叫“热透镜效应”。如果风速为零光斑中心持续加热热积累越来越严重如果存在横风被加热的空气团会顺着风往下游漂移形成一个像彗尾一样的温度尾流导致光斑中心向下游偏移同时产生非对称畸变。1.2 相位屏在这种问题里扮演什么角色严格来说热晕是一个三维热流体与光传播的耦合问题。光束在大气中一边传播一边不断加热沿途空气空气的温度场反过来又改变光束相位光束再带着这个新相位继续传播这种相互作用沿传播路径处处发生。要完全数值求解三维NS方程加光传播算力开销非常吓人工程上很少这么干。相位屏的思路是把整条路径上的相位扰动压缩到一个或几个薄平面上。比如从发射端到目标距离L我把中间这段大气对光程的扰动等效折算成一块“相位板”光束经过这块板时只改相位、不改振幅。这样光传播就变成了“真空衍射传播 相位屏”交替进行的标准数值格式。湍流仿真是这么干的热晕仿真也可以这么干。当然热晕相位屏通常不能只放一块因为热晕沿路径持续累积工程上一般会切成多层每层算一个热晕相位屏层与层之间用角谱法或傅里叶法做衍射传播。本文先聚焦单层相位屏怎么算多层扩展最后讲。1.3 热晕相位屏和湍流相位屏最大的区别很多人第一次写热晕仿真会下意识去翻“相位屏生成”的老代码找到Kolmogorov谱、von Kármán谱换一个谱型就开始生成。这是最典型的误区。湍流相位屏是随机场图案不可重复统计特性由功率谱密度决定而热晕相位屏在给定光强分布、风速和吸收系数下是确定性的它的“畸变形状”完全由光强沿风向上的累积积分决定。换句话说湍流屏是“天上下雪落到哪算哪”热晕屏是“水流过石头哪里流速快哪里就冲出一道沟”。随机相位屏可以用频域滤波器加随机数来生成热晕相位屏却必须先算温度场再做光程积分。这两者的物理逻辑完全不同算法框架也完全不同。下面我就从热传导方程出发把整个推导串一遍。2. 从传热方程到可计算的频域传递函数2.1 对流主导的温升方程我做的是稳态近似忽略热传导、忽略浮力对流只保留横向风把热量带走的机制。对于高功率激光连续波传输的典型场景这个近似是合理的。方程写出来很简单ρCp (v·∇) T(x,y) β I(x,y)其中ρ是空气密度Cp是定压比热容v是风速矢量T是温升β是大气吸收系数I是光强分布。这个方程的物理含义是空气微团沿风方向流动时沿途不断吸收激光热量所以温度升高。温度升高的速率取决于局地光强和吸收系数。假设风沿x轴正方向风速大小为V那么方程变成V * ∂T/∂x β I(x,y) / (ρCp)两边对x积分从上游无穷远积分到当前位置T(x,y) β / (ρCp V) * ∫_{-∞}^{x} I(s,y) ds这就是关键结果某一点的温升等于该点上游所有光强沿风向的累积。注意这里的积分只沿风向与光斑垂直方向的光强分布会自然形成温度场的横向轮廓。高斯光束中心光强最强所以温度场在中心处达到峰值边缘低这就是负透镜的来源。2.2 把积分变成频域的除法直接逐点做累加积分不是不行但风速方向一变或者需要跟Simulink的矩阵信号配合效率就低了。我更推荐在频域里做这个积分。对等式两侧做二维傅里叶变换积分运算 ∂^{-1} 在频域对应除以一个虚数项。先定义空间频率 kx 和 ky。如果风向与x轴夹角为θ那么沿风向的导数算子在频域对应i (kx cosθ ky sinθ)因此求温升就是T_hat I_hat / ( i (kx cosθ ky sinθ) )这里的 T_hat 和 I_hat 分别是温度场和光强场的二维傅里叶变换。看到这个式子你可能会担心分母等于零怎么办确实当 kx cosθ ky sinθ 0也就是垂直于风向的空间频率分量这个除法会爆炸。物理上这些分量对应“沿风方向无限长的等值带”它们在稳态下会产生无限累积的温度这属于非物理的直流分量直接在频域里把这一项强制置零就行。最后温度场乘以一个系数就得到相位屏φ(x,y) k * dn/dT * L * T(x,y)其中 k 2π/λ 是波数L是有效传输距离dn/dT是折射率温度系数。把前面的系数整合起来频域实现的核心就是一行传递函数G(kx,ky) 1 / ( i (kx cosθ ky sinθ) )这个G乘上光强的傅里叶变换再做逆变换再乘上系数就得到热晕相位屏。2.3 那个常数项怎么定我见过不少同学把公式推导看懂了结果算出来的相位屏数值大得离谱或小得可以忽略。根因几乎都在系数上。把完整系数写开coef k * dn_dT * β * L * dx / (ρ Cp V)注意这里有个 dx因为傅里叶变换对应的是求和而不是积分离散求和要把网格间距乘回去。很多人容易漏掉这个dx。另外光强I必须用真实的W/m²不能直接用归一化后的灰度值。大气吸收系数β的单位是1/m空气密度ρ约1.225 kg/m³Cp约1005 J/(kg·K)dn/dT取 -1e-6左右。把这些都代进SI单位制出来的相位就是弧度。还有一点容易混淆如果光沿z方向传输的整段距离L内光强近似不变那么相位就是把每个薄层dz上的温升累加一次等效于乘以L。如果你做的是多层相位屏那么每一层的L就是这一层的厚度而不是总距离。这块单独拎出来说是因为我在Simulink里做参数扫描时曾把L误设成总距离同一组数据跑了三遍条纹全对不上最后才发现是这里翻车了。3. MATLAB实现写一个工程可用的热晕相位屏生成函数3.1 函数接口和输入参数设计写这个函数之前我给自己定了几个要求输入输出单位统一参数全部显式传入不依赖全局变量方便后面Simulink批量扫描。最终接口是这样的function phi thermalBloomingPhaseScreen(I, lambda, dx, V, theta, beta, rho, Cp, dn_dT, L) % thermalBloomingPhaseScreen 计算稳态热晕相位屏 % 输入 % I - 二维光强分布单位 W/m^2 % lambda - 激光波长单位 m % dx - 网格间距单位 m % V - 风速大小单位 m/s % theta - 风向角单位 rad0表示沿矩阵列方向x轴吹 % beta - 大气吸收系数单位 1/m % rho - 空气密度单位 kg/m^3 % Cp - 空气定压比热容单位 J/(kg*K) % dn_dT - 折射率温度系数空气约 -1e-6单位 1/K % L - 等效传输距离单位 m % 输出 % phi - 热晕相位屏单位 rad尺寸与 I 相同 k 2 * pi / lambda; [Nr, Nc] size(I); fx (-Nc/2 : Nc/2-1) / (Nc * dx); fy (-Nr/2 : Nr/2-1) / (Nr * dx); [FX, FY] meshgrid(fx, fy); kx 2 * pi * FX; ky 2 * pi * FY; kproj kx * cos(theta) ky * sin(theta); G zeros(size(kproj)); valid abs(kproj) 1e-10; G(valid) 1 ./ (1i * kproj(valid)); Ihat fft2(I); Ihat_shifted fftshift(Ihat); Jhat Ihat_shifted .* G; J real(ifft2(ifftshift(Jhat))); coef k * dn_dT * beta * L * dx / (rho * Cp * V); phi coef * J; end这个函数的核心逻辑很简洁先算空间频率网格再构造沿风向的积分传递函数把光强变换到频域后乘传递函数再逆变换回空间域。注意我用了fftshift和ifftshift确保零频在中心构造传递函数时也是中心零频这个配对的坑后面专题讲。3.2 关键代码逐段解读构造频率网格时fx的长度是列数Ncfy的长度是行数Nr。这里必须和矩阵维度对应好否则最终相位屏会转置。我的习惯是行方向对应y轴列方向对应x轴这样矩阵形状就是直接的标量场不用反复转置。传递函数里最需要注意的就是valid掩膜。我把abs(kproj) 1e-10的项置零这是为了避免分母为零。但这里有个经验值问题阈值取多少合适取太大低频分量被削掉相位屏会出现“缺了骨架”的波浪形取太小数值上容易出NaN。我测试下来在典型的128x128网格下取1e-10这个绝对阈值可行但如果网格尺寸变成1024建议用相对阈值delta 1e-6 * max(abs(kproj(:))); valid abs(kproj) delta;这个做法更稳健因为网格变密后频率单位不变但零频附近的离散频率间隔变小了绝对阈值可能把不该削的成分一起削掉。还有一个细节相位屏只取real(ifft2(...))。理论上频域算子共轭对称逆变换结果应该是实部但由于浮点误差虚部会有一点噪声直接丢弃即可。如果你发现相位屏虚部很大说明传递函数构造不对称多半是fftshift和ifftshift配错了。3.3 任意风向的处理扩展上面的代码已经支持任意风向角 θ。θ0 表示风沿矩阵列方向也就是x轴θπ/2 表示风沿矩阵行方向。实际场景里风向一般不是正交的这时候kproj kx*cos(θ) ky*sin(θ)就自动完成了沿任意方向的线积分。这里有个容易误解的点频域积分是“全平面”积分等价于从风向上游无穷远积分到当前位置这是正确的物理极限但实际热晕只在下游有尾流上游还没有被加热。当光强分布是有限尺寸的高斯光斑时这个区别不大因为上游无穷远处的光强为0积分贡献为0。但如果你的光强矩阵边界上还有非零值周期延拓会让光强从另一侧“绕回来”造成虚假的尾流。解决办法很简单计算相位屏之前把光强矩阵外围补一圈零算完再裁掉边缘。这个操作在工程上非常实用我建议所有用户都加上。3.4 边界效应和正则化频域积分自带周期边界这是FFT方法固有的。直观表现是沿风向的尾流应该从左边界一直延伸到右边界但周期边界会让尾流从右边界消失后又从左边界冒出来。当你用128x128网格仿真一个小尺寸光斑时光斑周围本来就有大块零区域这种绕回效应不明显但如果光斑填满了网格尾流绕回就会形成假的条纹。对策有两个一是扩边推荐扩到光斑直径的3到4倍二是用更大的网格并保持光斑居中。扩边之后边界绕回的假信号被限制在远离光斑的区域而你关心的光斑区域相位是准的。我一般会在仿真初始化脚本里统一处理比如N 256; I_large zeros(N, N); I_large(64:191, 64:191) I_original; dx_large dx; % 网格间距不变 phi_large thermalBloomingPhaseScreen(I_large, lambda, dx_large, V, theta, beta, rho, Cp, dn_dT, L); phi phi_large(64:191, 64:191);这样既避免边界绕回又没有改变网格间距后续和衍射传播模块对接也方便。4. Simulink集成把相位屏变成可批量扫描的仿真模块4.1 MATLAB Function模块封装相位屏计算单纯用MATLAB脚本算相位屏很直接但要把热晕相位屏放进一个更大的仿真链路里比如配合Simulink的光束控制模型、伺服系统模型一起跑就需要把相位屏算法集成到Simulink中。我惯用的做法是用MATLAB Function块调用上面的函数。在Simulink模型里放一个MATLAB Function块双击进去写代码function phi TBPhase(I) coder.extrinsic(thermalBloomingPhaseScreen); p coder.load(thermalBloomingParams.mat); phi zeros(size(I)); phi thermalBloomingPhaseScreen(I, p.lambda, p.dx, p.V, p.theta, ... p.beta, p.rho, p.Cp, p.dn_dT, p.L); end这里coder.extrinsic是关键。它的作用是告诉Simulink这个函数不求代码生成而是在仿真的时候用MATLAB解释器执行。这样MATLAB Function块就可以调用普通m函数不用担心fft2在代码生成环境下的兼容性问题。缺点是仿真速度会变慢但对于以数据处理为主的相位屏模块来说128x128左右的矩阵一次调用也就几毫秒完全能接受。如果模型里已经有实时光强信号比如从一个波前传感器模块输出光强矩阵直接把信号线接到这个MATLAB Function的输入口就行。需要注意Simulink里矩阵信号的尺寸是固定的所以输入I的维度在仿真过程中不能变化。我一般在模型初始化回调里固定I_dim 128然后把所有相关模块的维度都设为这个值。4.2 用sim命令批量扫风速和功率把相位屏模块封装好后最大的好处是可以脚本化批量扫描。比如我想看风速从1m/s到20m/s变化时热晕相位屏的RMS和Strehl比怎么变不需要改Simulink模型只需要在MATLAB脚本里循环调用sim。我用的是这种写法V_list [1 2 5 10 15 20]; for i 1:length(V_list) set_param(thermalBloomingSim/windV, Value, num2str(V_list(i))); simOut sim(thermalBloomingSim, StopTime, 0, ... ReturnWorkspaceOutputs, on, SaveOutput, on); phi squeeze(simOut.yout{1}.Values.Data); rms_phi(i) std(phi(:)); strehl(i) exp(-rms_phi(i)^2); end模型里我把风速参数做成一个Constant模块变量名windV放在MATLAB Function块前面作为参数传给相位屏函数。循环里用set_param修改该Constant的值再用sim运行。这样同一个模型可以扫任意工况不需要为每个工况手动搭模型。如果你的MATLAB版本较新simOut.yout{1}.Values.Data的取法可能略有差异建议先用simOut.who看一下输出对象里有哪些变量。这种脚本化批扫思路比在Simulink里手动改参数再点运行效率高两个数量级也是我做参数敏感性分析的主力工具。4.3 后续接远场衍射的扩展有了热晕相位屏下一个自然的需求是看这个相位屏对光束远场光斑的影响。这时可以把相位屏看成是一个纯相位物体入射复振幅A_in乘上exp(1i*phi)然后做一次夫琅禾费衍射。在Simulink里也可以加一个MATLAB Function块做FFT衍射function I_far farField(U_in, dx, z, lambda) % U_in: 入射复振幅场 % dx: 空间采样间距 % z: 传播距离 % lambda: 波长 N size(U_in, 1); fx (-N/2 : N/2-1) / (N*dx); U_out fftshift(fft2(ifftshift(U_in))); % 夫琅禾费衍射的频率-坐标映射 x_far lambda * z * fx; I_far abs(U_out).^2; I_far I_far / sum(I_far(:)); end调用时入射复振幅可以是A_0 .* exp(1i*phi)其中A_0是初始高斯振幅分布。这样Simulink模型就变成了一条完整链路光强 - 热晕相位屏 - 远场光斑。我在实际项目里还加过自适应光学校正模块在这个链路上做补偿算法验证效果非常好。5. 不同条件下热晕相位屏的对比结果与评价指标5.1 典型工况设置为了直观展示程序的效果我用下面这组参数跑了一个单层热晕相位屏的对照实验参数数值波长 lambda1.064e-6 m网格数 Nx x Ny256 x 256网格间距 dx0.002 m光束半径 w00.05 m峰值光强 I01e7 W/m²对应总功率约75kW吸收系数 beta1e-4 /m空气密度 rho1.225 kg/m³空气比热 Cp1005 J/(kg·K)折射率温度系数 dn_dT-1e-6 /K等效传输距离 L1000 m风速 V1 / 5 / 20 m/s风向角 theta0 rad光强分布直接用高斯光束的公式生成[X, Y] meshgrid((1:Nx)*dx, (1:Nx)*dx); I I0 * exp(-2 * ((X-Nx*dx/2).^2 (Y-Nx*dx/2).^2) / w0^2);5.2 三个典型结果怎么读风速1m/s时相位屏呈现非常明显的“彗尾”结构。沿风向方向相位从光斑中心靠上游的位置开始快速变化到下游逐渐拖出一条长长的尾巴。相位峰值可能到几个弧度量级光束畸变非常严重。此时热晕效应的核心表现是光斑中心区域等效为强负透镜光束在远场扩展非常厉害同时光斑中心向下游偏移。风速5m/s时同样的光强和吸收条件下温度累积的梯度被风拉平了一部分相位屏的峰值明显变小。尾流仍然存在但长度更长、幅度更缓。用专业一点的话说热晕的“相位振幅”减小了但“相位梯度”也降低了因此对光束的偏折作用减弱。此时远场光斑的扩展会比1m/s小很多但也比无热晕时严重。风速20m/s时热量被风快速带走温度场几乎没有什么累积相位屏的PV值可能降到0.1rad以下。这时热晕不再是主要限制因素光束畸变主要由湍流或者系统本身的像差决定。这个趋势也符合实际经验热晕效应在无风或弱风条件下最严重一旦风速起来效应迅速减弱。5.3 用Strehl比和RMS快速评价定量评价热晕相位屏有多严重我常看两个指标相位屏的RMS值和Strehl比近似值。Strehl比这里用近似公式Strehl ≈ exp(-σ_phi²)其中σ_phi是相位屏RMS。这个近似成立的前提是相位均值为0且方差不大工程上够用。计算代码phi_centered phi - mean(phi(:)); rms_phi sqrt(mean(phi_centered(:).^2)); strehl exp(-rms_phi^2);我的经验阈值大概是这样RMSradStrehl 估算主观判断 0.1 0.99可忽略0.1 ~ 0.30.91 ~ 0.99轻度热晕0.3 ~ 10.37 ~ 0.91明显热晕 1 0.37严重热晕当然这只是针对单相位屏的粗略评价实际系统还要看光束质量和环围能量等指标但作为快速判断已经足够了。在批量扫描过程中我会每算一个工况就把RMS和Strehl记下来画成曲线一目了然。6. 我在调试这套程序时踩过的几个坑6.1 第一坑FFT零频和奇点我最早写的版本相位屏一出来就是一片巨大的常数加竖条纹。查了半天发现是传递函数构造时零频项没有处理。由于1/(i*kproj)在零频处是无穷大逆变换后会在整个画面上叠一个极大的常数背景而那条竖条纹则来自kx0那一列没有完全置零。把零频和kproj≈0的位置显式设置为0之后画面立刻正常了。这个小问题折磨了我一个下午所以我强烈建议写完传递函数先surf(real(G))看一眼确认零频处是0而不是NaN或Inf。6.2 第二坑单位换算翻车我之前说系数里有dx这就是踩坑后的教训。第一次实现只算了k*dn_dT*beta*L/(rho*Cp*V)出来的相位屏小到看不出形状。后来把dx加进去量级才合理。为什么必须乘dx因为离散FFT是求和而物理量是积分。MATLAB没有现成的“频域积分算子”可以自动带间距所以这个dx必须手动乘。同样如果光强有单位但数值很大比如峰值1e7 W/m²相位可能会达到几十弧度这时候要注意数值溢出或显示范围问题。6.3 第三坑Simulink矩阵信号尺寸爆炸Simulink里接矩阵信号时如果输入光强矩阵是 256x256输出相位屏也必须是 256x256。但你从MATLAB Function块拉输出线时Simulink有时会推断成“可变大小信号”这时候需要手动指定输出为固定大小。否则一旦某次仿真输入尺寸变化模型直接报错。我自己的习惯是在MATLAB Function块内部一开始就用phi zeros(size(I))这种方式让 Simulink 识别输出尺寸并且不要在模型里通过Signal Specification强制改变矩阵维度。如果实在不行就把输入信号从矩阵改成打包的结构体或者用Matlab System Block但复杂度会高一点。6.4 扩展方向的一点想法这套程序目前是单层稳态热晕相位屏。如果要更贴近实际可以往两个方向扩展一是多层相位屏沿传播路径切N层每层算一个热晕相位屏层之间用角谱法衍射传播这样能模拟热晕沿路径累积的完整过程二是瞬态热晕考虑风场随机性和时间演化在频域传递函数中加一个时间因子配合Simulink的连续采样跑动态过程。我最近的实验是把热晕相位屏和湍流相位屏叠加在一起同一个网格里先加随机湍流屏再加热晕确定性相位这样仿真更接近真实大气环境。如果你也打算写自己的热晕相位屏程序我强烈建议先用一个均匀光强分布做自检均匀光强不会产生横向温差相位屏应该几乎为零。如果程序输出一个很大的相位那一定是积分方向或零频处理出了问题。这个自检虽然简单但能帮你省掉大量排查时间。本文还有配套的精品资源点击获取