简介本资源是一套面向光学工程、激光技术及大气物理方向研究者的MATLAB热晕相位屏仿真工具包聚焦高能激光传输与天文观测中由大气温度不均匀性引发的波前畸变建模问题。程序完整实现相位屏生成、二维傅里叶变换含自定义myfft2/myifft2、热晕效应下的光学传递函数OTF与点扩散函数PSF计算、以及成像质量定量分析等核心流程支持用户通过调节大气参数模拟不同热晕强度场景。压缩包共13个文件7个.m主程序脚本承担算法主体4幅.bmp图像用于输入/输出可视化验证1个.txt说明文档提供关键参数注释总大小541KB结构紧凑、模块清晰便于理解热晕物理机制并快速开展二次开发。目前已有252人学习下载适用于具备基础MATLAB编程能力的研究生与工程师开展热晕效应机理研究、补偿算法预研或教学演示。1. 项目背景与热晕效应初探最近在整理一个老项目翻出来一个名为“热晕程序.zip”的压缩包里面是关于大气热晕效应的MATLAB相位屏仿真程序。这让我想起了当年做激光大气传输研究时被这个“热晕”现象折腾得够呛的经历。简单来说当高功率激光束在大气中传播时它会吸收空气中的气溶胶和分子导致光束路径上的空气被加热。被加热的空气密度降低折射率随之减小这就形成了一个天然的“负透镜”。这个透镜效应会使激光束本身发生扩散、畸变甚至导致光束中心能量严重衰减就像人被热晕了头一样视线模糊、行动不稳所以这个现象被形象地称为“热晕”Thermal Blooming。对于从事高能激光、激光通信、自适应光学或者大气物理研究的朋友来说热晕效应是一个必须跨越的坎。你不能只设计一个功率强大的激光器却不管它在实际大气中传播时会不会“自毁前程”。而研究热晕最经济、可控且能深入机理的方法就是数值仿真。MATLAB凭借其强大的矩阵运算能力和丰富的可视化工具成为了实现这类物理场仿真的利器。这个仿真程序的核心就是利用“相位屏”方法来模拟激光光束在受热大气中传播时其波前相位发生的动态、非线性畸变过程。接下来我就结合这个老程序拆解一下如何用MATLAB搭建一个能跑起来、且能说明问题的热晕相位屏仿真模型并分享一些当年调试时踩过的坑和心得。2. 热晕相位屏仿真的核心物理与数学模型要仿真首先得把物理过程用数学语言描述清楚。热晕不是一个简单的瞬时效应它是一个流体力学与电磁传播耦合的时变过程。我们通常基于以下几个关键假设来建立简化模型以便于在计算机上实现2.1 基本假设与控制方程首先我们通常假设激光是连续波CW或者长脉冲这样有时间建立稳定的热场。其次忽略横向风或者假设风场是均匀稳定的这能极大简化模型。最重要的我们使用“薄屏近似”和“线性热晕”理论。薄屏近似意味着我们把大气对激光的影响集中到若干个垂直于传播方向的平面上即相位屏光束在自由空间传播只在经过这些屏时才发生相位突变。线性热晕理论则假设折射率变化与吸收的激光能量密度成正比。核心的耦合过程可以概括为激光吸收生热激光强度I(x, y, z, t)被大气吸收吸收系数为α产生的体热源功率密度为Q α * I。热传导与对流产生的热量通过空气的热传导和风场风速v的对流来散开。在稳态或准稳态假设下这个热扩散过程可以用一个泊松方程类的方程来描述。一个非常经典的简化模型是“等压、稳态、有风”模型它给出折射率扰动Δn与激光强度I的积分关系。对于横向风沿x方向的情况有Δn(x, y, z) ∝ -∫_{-∞}^{x} I(x’, y, z) dx’在特定简化下。 这个公式直观地反映了“上游累积”效应光束上某点的折射率变化正比于该点上游风来的方向所有点吸收的激光能量的累积。这正是热晕导致光束向上风方向弯曲和畸变的根源。相位延迟折射率变化Δn导致光束通过该区域后产生额外的相位延迟φφ(x, y) (2π / λ) * ∫ Δn(x, y, z) dz。 在薄屏近似下这个沿传播方向z的积分被简化为相位屏上的一个二维相位分布φ(x, y)。光束传播携带了相位畸变φ(x, y)的光束通过下一段自由空间传播到下一个相位屏位置。这个传播过程通常用角谱法Angular Spectrum Method或菲涅尔衍射积分来实现。2.2 数值仿真中的离散化与迭代在计算机里我们需要把连续的物理场离散化。将光束横截面x-y平面离散为N×N的网格。时间也被离散为步长Δt。仿真流程是一个典型的“预测-校正”迭代过程步骤一初始给定初始时刻t0的光束复振幅场U0(x, y)包含振幅和相位。步骤二传播U0通过一段自由空间传播到第一个相位屏位置。步骤三加热与相位计算根据当前屏位置的光强I |U|^2结合风速、吸收系数等参数利用离散化的积分公式如梯形法、辛普森法计算该屏上的折射率扰动Δn进而得到相位扰动φ。步骤四相位施加将计算得到的相位屏乘以光束复振幅场U’ U * exp(i * φ)。步骤五继续传播畸变后的光束U’继续传播到下一个相位屏或观测面。步骤六时间推进更新时间为t Δt。此时由于风的存在大气“被加热”的状态是移动的。我们需要更新每个相位屏上的“热记忆”。一种常见方法是根据风速v将上一个时刻的加热图样在空间上进行平移相当于风把热吹走了然后叠加上新时刻在当前位置产生的加热。这构成了下一个时间步的初始热场。循环重复步骤二到步骤六直到达到设定的仿真时间。这个循环清晰地揭示了热晕的动态性和非线性当前时刻的光束形状决定了加热分布加热分布改变了相位屏相位屏又影响了下一时刻光束的传播如此反复耦合。3. MATLAB仿真程序的结构化实现与关键代码解析光有理论不够还得能写成代码。一个结构清晰的热晕仿真程序通常包含以下几个模块。我会结合关键代码片段使用MATLAB语法进行说明并解释为什么这么写。3.1 主程序框架与参数初始化主程序通常是一个脚本负责调用各个函数并管理迭代循环。参数初始化是第一步也是容易出错的地方。务必保证所有物理量使用国际单位制SI并注意量纲一致性。%% 热晕仿真主程序 - 参数设置 clear; close all; clc; % 物理参数 lambda 1.06e-6; % 激光波长 [m] 常见1.06微米Nd:YAG w0 0.05; % 光束初始束腰半径 [m] P0 1e6; % 激光功率 [W] 1MW alpha 1e-3; % 大气吸收系数 [1/m] v_wind 5.0; % 横向风速 [m/s] 沿x正方向 Cn2 1e-15; % 大气湍流折射率结构常数 [m^{-2/3}] 用于对比或叠加湍流效应 L_prop 5000; % 总传播距离 [m] % 计算与大气热晕相关的特征参数用于评估和缩放 % 热晕畸变参数 N_D Distortion Parameter N_D越大热晕越强 % N_D ∝ (α * P * L^2) / (ρ * Cp * v * w0^3) 这里省略了常系数用于定性比较 rho 1.2; Cp 1005; % 空气密度 [kg/m^3] 和定压比热容 [J/(kg·K)] N_D_est (alpha * P0 * L_prop^2) / (rho * Cp * v_wind * w0^3); fprintf(估计的热晕畸变参数 N_D ≈ %.2f\n, N_D_est); if N_D_est 1 fprintf(提示N_D 1 热晕效应预计将非常显著。\n); end % 仿真网格参数 N 256; % 网格点数 (建议为2的幂次便于FFT) L_grid 0.5; % 计算窗口的物理尺寸 [m] 要能容纳扩散后的光束 dx L_grid / N; % 空间采样间隔 [m] x (-N/2:N/2-1) * dx; % 空间坐标向量 [m] [X, Y] meshgrid(x, x); % 传播与相位屏参数 dz 100; % 相位屏间隔 [m] z_positions 0:dz:L_prop; % 相位屏位置向量 num_screens length(z_positions) - 1; % 相位屏数量起点不算屏 dt 0.01; % 时间步长 [s] 需满足 CFL 类条件 v_wind*dt dx num_steps 200; % 总时间步数注意参数N_D是一个无量纲数是评估热晕严重程度的关键。在程序开头进行估算可以帮助你预判仿真结果是否合理。网格大小L_grid需要仔细选择要确保在整个传播过程中光束及其旁瓣不会触及计算窗口边界否则会发生混叠Aliasing。3.2 初始光束场与传播函数生成初始光束通常采用基模高斯光束。角谱传播法是实现自由空间衍射的稳定且准确的方法。%% 生成初始高斯光束场 k 2 * pi / lambda; % 波数 % 在z0面束腰处的高斯场 U0 sqrt(2*P0/(pi*w0^2)) * exp(-(X.^2 Y.^2)/w0^2); % 振幅 已归一化功率 U0 U0 .* exp(1i * 0); % 初始相位设为平面波 %% 角谱传播函数 (用于自由空间传播段) % 注意角谱法在频域操作需要构建空间频率坐标 fx (-N/2:N/2-1) / (N*dx); % 空间频率 [1/m] [FX, FY] meshgrid(fx, fx); H_prop (dz) exp(1i * k * dz * sqrt(1 - (lambda*FX).^2 - (lambda*FY).^2)); % 传递函数 % 为了避免sqrt内出现负数对应倏逝波通常做傍轴近似 H_prop_paraxial (dz) exp(1i * k * dz) .* exp(-1i * pi * lambda * dz * (FX.^2 FY.^2)); % 对于大多数激光大气传输傍轴近似足够精确且计算稳定。这里提供了精确和傍轴近似两种传递函数。对于热晕仿真传播距离通常远大于光束尺寸傍轴近似即菲涅尔衍射的频域形式在计算效率和精度上是一个很好的折中。使用函数句柄(dz)方便了对不同距离的传播进行调用。3.3 热晕相位屏计算的核心函数这是程序的“心脏”。它根据当前的光强分布、风速和吸收系数计算相位屏。我们实现前面提到的“上游累积”模型。function phase_screen calc_thermal_blooming_phase(I, x, y, v, alpha, dt, lambda, history) % 计算热晕引起的相位屏 % 输入 % I - 当前时刻在相位屏位置的光强分布 [W/m^2] % x, y - 网格坐标向量 [m] % v - 风速标量假设沿x方向 [m/s] % alpha - 吸收系数 [1/m] % dt - 时间步长 [s] % lambda - 波长 [m] % history - 上一个时刻的“累积加热图”与I同尺寸首次调用时为零 % 输出 % phase_screen - 相位屏分布 [rad] % 可选更新后的 history dx x(2) - x(1); [X, Y] meshgrid(x, y); % 1. 计算当前时间步产生的“新加热” % 简化模型加热量正比于光强乘以吸收系数和时间步长 dQ alpha * I * dt; % [J/m^3] 量级 这里忽略了密度和比热容的常系数因其最终会归入一个比例常数 % 2. 风场效应将历史热场沿风向-x方向平移 % 风速v0表示风沿x方向吹那么热斑会向-x方向移动。 shift_pixels round(v * dt / dx); % 计算需要平移的像素数 if shift_pixels ~ 0 history_shifted circshift(history, [0, -shift_pixels]); % 注意方向 % 对于被移出计算窗口的部分我们简单丢弃相当于热被吹走了。 % 更精细的模型可能需要考虑边界处理。 else history_shifted history; end % 3. 更新累积加热图旧的热被吹走一部分加上新的加热 Q_cumulative history_shifted dQ; % 4. 计算折射率扰动 (简化线性模型Δn ∝ -Q_cumulative) % 比例常数 K 综合了空气的热光系数等物理常数。 % 对于标准空气一个常用的近似是 Δn ≈ -1e-6 * (ΔT)而 ΔT ∝ Q/(ρ*Cp) rho 1.2; Cp 1005; K -1e-6 / (rho * Cp); % 这是一个示例性的比例常数实际值需根据具体物理模型校准 delta_n K * Q_cumulative; % 5. 计算相位延迟 (薄屏近似相位 k * Δn * 屏厚度) % 屏厚度取为相位屏间隔 dz。注意这里的 dz 是调用此函数时外部决定的。 % 为了模块化我们将dz作为另一个输入参数或者在主循环中计算相位后乘以dz。 % 本例中我们在函数外部处理dz的乘法因此这里只返回与Δn成正比的量。 phase_screen (2*pi/lambda) * delta_n; % 注意这里缺少了 * dz 将在主循环中补上 % 6. 可选返回更新后的累积加热图供下一时间步使用 % 在函数调用中如果需要更新history可以设计为返回两个值。 end这个函数实现了动态热晕的核心circshift操作模拟了风对热分布的平移dQ的累加模拟了持续的激光加热。比例常数K是连接物理模型和仿真结果的关键它需要根据实际的物理参数如空气的热光系数dn/dT来精确计算或通过实验数据标定。在学术研究中这个常数常常被归一化或者直接使用N_D参数来缩放整个效应。3.4 主时间迭代循环与可视化将以上模块组装起来形成动态仿真。%% 主仿真循环 U_current U0; % 初始化光场 % 初始化每个相位屏位置的“累积加热历史” Q_history cell(1, num_screens); for i 1:num_screens Q_history{i} zeros(N, N); end % 预分配数组存储结果例如在目标处的光强随时间变化 I_target_time zeros(N, N, num_steps); figure(Position, [100, 100, 1200, 400]); for t_step 1:num_steps fprintf(时间步 %d / %d\n, t_step, num_steps); U_prop U_current; % 分段传播通过每个相位屏 for screen_idx 1:num_screens % 1. 传播到当前相位屏位置如果是第一个屏就是从起点传播dz prop_dist dz; U_prop fftshift(fft2(U_prop)); % 到频域 U_prop U_prop .* H_prop_paraxial(prop_dist); % 角谱传播 U_prop ifft2(ifftshift(U_prop)); % 回空域 % 2. 在屏位置计算当前光强 I_at_screen abs(U_prop).^2; % 3. 调用函数计算热晕相位屏 phase_pert calc_thermal_blooming_phase(I_at_screen, x, y, v_wind, alpha, dt, lambda, Q_history{screen_idx}); phase_pert phase_pert * dz; % 补上屏厚度带来的相位累积 % 4. 应用相位屏 U_prop U_prop .* exp(1i * phase_pert); % 5. 更新该屏的累积加热历史 (这里需要函数能返回更新后的history) % 假设我们修改了函数使其返回[phase_screen, Q_updated] % [~, Q_history{screen_idx}] calc_thermal_blooming_phase(...); % 为简化演示此处省略更新步骤实际代码必须包含。 end % 传播到最后的目标面最后一段dz U_prop fftshift(fft2(U_prop)); U_prop U_prop .* H_prop_paraxial(dz); U_target ifft2(ifftshift(U_prop)); % 存储目标面光强 I_target_time(:, :, t_step) abs(U_target).^2; % 实时可视化每10步显示一次避免图形刷新过慢 if mod(t_step, 10) 0 subplot(1, 3, 1); imagesc(x, x, log10(I_at_screen 1e-10)); % 显示最后一个相位屏处的光强对数尺度 axis image; colorbar; title(sprintf(屏处光强 (t%.2fs), t_step*dt)); xlabel(x (m)); ylabel(y (m)); subplot(1, 3, 2); phase_to_show angle(exp(1i*phase_pert)); % 包装相位到[-pi, pi] imagesc(x, x, phase_to_show); axis image; colorbar; title(最后一个热晕相位屏); xlabel(x (m)); ylabel(y (m)); subplot(1, 3, 3); imagesc(x, x, log10(I_target_time(:, :, t_step) 1e-10)); axis image; colorbar; title(目标面光强); xlabel(x (m)); ylabel(y (m)); drawnow; end % 为下一个时间步准备将当前目标场作为下一个时间步的初始场这取决于模型 % 对于“准直光束经过静态介质”的常用模型我们通常从光源重新开始传播。 % 但为了模拟风场带来的动态变化我们需要更新每个屏的Q_history。 % U_current U0; % 简单重置为初始光束 end主循环清晰地展示了“传播-加热-相位施加-再传播”的迭代过程。可视化部分非常重要它能让你直观地看到光束如何从圆形高斯光斑逐渐被“吹”向上风方向-x方向并变得不对称、发散。4. 仿真调试中的关键问题与实战心得写完代码只是第一步让仿真结果物理可信、数值稳定才是真正的挑战。以下是我在多次调试中总结的几个核心要点和踩过的坑。4.1 空间与时间采样满足奈奎斯特与CFL条件这是数值计算的基础也是最容易导致诡异结果的地方。空间采样dx你的网格间距dx必须足够小以分辨光束中最小的特征结构。对于高斯光束主要能量集中在束腰w0内但为了准确计算衍射和相位屏的高频成分经验法则是dx ≤ λ/2或更小。同时计算窗口大小L_grid必须足够大确保光束在传播过程中尤其是发生热晕扩散后不会触及边界。一个实用的检查方法是仿真结束后检查光束边缘处的光强是否已经衰减到接近零例如小于峰值光强的1e-6。如果边界处还有显著能量说明窗口太小会发生混叠误差。时间采样dt由于我们引入了风场v_wind这本质上是一个对流问题。时间步长dt必须满足CFLCourant-Friedrichs-Lewy稳定性条件的一个变体v_wind * dt dx。这意味着在一个时间步内风场移动的距离不能超过一个空间网格。否则circshift操作会跳过网格点导致热量“瞬移”破坏物理真实性并可能引发数值不稳定。通常取dt 0.5 * dx / v_wind是比较安全的选择。4.2 相位屏的“包装”问题与处理方法计算出的相位屏φ(x, y)可能包含很大的数值其范围远超过[-π, π]。当使用exp(i*φ)施加相位时MATLAB的exp函数内部会处理复数的相位部分理论上φ多大都可以。但是当你需要可视化相位屏时直接使用imagesc(phi)会得到一个颜色条范围极大的图看不出细节结构。此时需要对相位进行“包装”Wrapping将其映射到[-π, π]区间phi_wrapped angle(exp(1i * phi))。这仅用于可视化不影响实际计算。注意在有些涉及干涉或相位复原的算法中需要处理真实的“包裹相位”但我们的仿真中φ是真实计算的相位延迟直接用于计算exp(i*φ)即可包装只是为了看图方便。4.3 功率守恒校验与能量流失分析一个重要的物理自洽性检查是光束总功率能量是否守恒。在理想的无吸收、无散射的自由空间传播中角谱法是严格功率守恒的。但是我们的模型引入了“吸收”来产生热晕这本身意味着光功率会随着传播而衰减。因此我们需要区分程序性功率流失由于数值误差如FFT的舍入误差、窗口截断导致的非物理功率变化。这应该很小。物理性功率衰减由吸收系数α模型决定的衰减。我们的简化模型Δn ∝ ∫ I dx’隐含了吸收但在代码中我们并没有在光场U的振幅上直接乘以一个衰减因子exp(-α*z/2)。这是因为在我们的模型里吸收的能量已经转化为热并影响了相位但光束振幅的衰减没有被显式模拟。这是一个模型近似。更自洽的模型需要将吸收导致的振幅衰减和加热导致的相位扰动分开处理或者使用更复杂的耦合波方程。一个实用的校验方法是监测光束在传播过程中横截面上积分光强sum(sum(I))*dx*dy的变化。如果除了在施加了显式吸收衰减的步骤外这个值发生剧烈变化那很可能就是数值不稳定或边界处理出了问题。4.4 从仿真到物理参数标定与结果解读仿真程序跑通了出来的图也挺好看但怎么知道它对不对呢这就需要将仿真结果与理论预测或实验现象进行对比。特征量验证计算仿真结果中的光束质量因子如Strehl比下降、光束重心偏移量、环围能量半径等与基于N_D参数的理论公式进行趋势性对比。例如理论表明Strehl比随N_D增大而近似按1/(1N_D^2)下降。你可以固定其他参数改变激光功率P0从而改变N_D运行一系列仿真绘制Strehl比随N_D变化的曲线看是否符合理论趋势。极限情况测试设置风速v_wind为极大值如1e6热晕效应应几乎消失光束应接近衍射极限传播。设置吸收系数α为0应无任何热晕相位屏光束按常规衍射传播。设置功率P0非常小热晕效应应可忽略。与文献对比查找已发表的、使用类似模型的热晕仿真论文对比他们文中展示的光束形态演变图、峰值光强衰减曲线等看是否定性一致。5. 程序优化与功能扩展思路基础模型跑起来后可以考虑以下优化和扩展让仿真更高效、更贴近实际。5.1 计算性能优化向量化与并行化热晕仿真是计算密集型任务尤其是当网格数N较大、时间步数多时。优化是关键。避免循环中的重复计算例如角谱传递函数H_prop_paraxial对于固定的dz是常数可以在循环外计算一次并存储。使用parfor进行并行循环时间迭代步t_step之间通常是独立的取决于模型可以考虑使用parfor并行计算多个时间步。但要注意这需要将每个时间步的仿真完全独立并且Q_history的更新逻辑可能需要调整因为parfor循环顺序是不确定的。更常见的是对多个不同参数如不同风速、功率的仿真进行并行。使用GPU加速MATLAB支持GPU计算。可以将大型矩阵如U_prop,I_at_screen转换为gpuArray这样FFT、点乘等操作会在GPU上执行速度可能有数量级的提升。但要注意GPU内存限制。if gpuDeviceCount 0 X_gpu gpuArray(X); Y_gpu gpuArray(Y); U0_gpu gpuArray(U0); % ... 后续计算使用这些 GPU 数组 % 记得在需要时将结果 gather 回 CPU end5.2 模型进阶引入湍流与动态风场真实大气中热晕和湍流是同时存在的且风场也不是均匀稳定的。叠加湍流相位屏可以在每个相位屏位置在热晕相位φ_thermal上叠加一个由大气湍流生成的随机相位屏φ_turbulence。生成湍流相位屏有成熟的方法如功率谱反演法使用符合Kolmogorov谱的随机相位屏或Zernike多项式展开法。这能研究热晕与湍流的耦合效应比如湍流是否会“打散”热晕形成的稳定热透镜从而部分缓解热晕。时变风场与剪切风将风速v_wind从一个标量变为一个随时间t变化的函数v_wind(t)甚至是一个随高度z变化的剖面v_wind(z)风剪切。这需要修改calc_thermal_blooming_phase函数中的风平移逻辑可能需要对每个屏位置采用不同的风速并更精细地处理热量的三维输运。5.3 仿真结果的后处理与定量分析除了看光强分布图定量分析更能说明问题。光束质量分析% 计算环围能量半径 (桶中功率 PIB) I_target I_target_time(:, :, end); % 取最终时刻目标面光强 total_power sum(I_target(:)) * dx^2; [cx, cy] meshgrid(x, x); centroid_x sum(sum(cx .* I_target)) / total_power; centroid_y sum(sum(cy .* I_target)) / total_power; r_squared (cx - centroid_x).^2 (cy - centroid_y).^2; % 计算包含63.2%能量的半径类似于高斯光束的束宽 sorted_vals sort(I_target(:), descend); cumulative_power cumsum(sorted_vals) * dx^2; idx_63 find(cumulative_power 0.632 * total_power, 1); threshold_val sorted_vals(idx_63); beam_radius sqrt(mean(r_squared(I_target threshold_val)));Strehl比计算Strehl比定义为实际光束的峰值光强与理想衍射极限光束峰值光强之比。你需要先运行一个无热晕α0的仿真得到理想峰值光强I_ideal_peak然后与有热晕时的峰值光强I_actual_peak相比SR I_actual_peak / I_ideal_peak。这是衡量系统性能退化最直接的指标之一。回过头看这个“热晕程序.zip”它不仅仅是一段MATLAB代码更是一个理解复杂物理耦合过程的窗口。从最简单的线性稳态模型出发逐步加入动态风场、湍流、自适应光学补偿这又是一个大话题可以通过在相位屏上减去一个共轭相位来实现你可以构建出一个非常强大且贴近实际的大气激光传输仿真平台。调试这类程序最大的收获不是学会了某个MATLAB函数而是锻炼了一种将连续物理现象离散化、数值化并通过计算实验进行验证和探究的思维方式。每一次参数调整后观察到的光束形态变化都是对理论公式最直观的注解。本文还有配套的精品资源点击获取