简介：这份基于MATLAB的热晕相位屏仿真程序包，面向光学工程、激光传输及望远镜设计领域的研究人员与工程师，用于模拟大气温度不均导致的光波前畸变并评估成像质量。压缩包共13个文件，约541KB，包含7个m脚本、4个bmp图像、1个asv备份文件及1个txt说明文档；m脚本实现了相位屏随机生成、傅里叶变换、OTF与PSF计算等核心流程，bmp图像为仿真结果示例，说明文档便于快速上手。程序支持调整大气条件、温度梯度等输入参数，可对比不同场景下的点扩散函数与成像效果，从而量化热晕对分辨率的影响。已有254人学习下载，适合需要开展高能激光传输或大气光学仿真分析的MATLAB用户，也可作为教学与课题研究的参考工具。
1. 热晕相位屏仿真：一套能直接跑出真实热晕形态的 MATLAB 程序
做高能激光传输链路的人，对热晕这个词应该都不陌生。高功率激光在稠密大气里走几百米，大气吸收能量后升温、折射率下降，光束像穿过一块自生的负透镜，光斑从圆斑被压成新月形甚至中间凹陷的环——这个现象在 MATLAB 里要建出来，靠的就是热晕相位屏仿真程序。这份资源我拆完之后给同事复现过，核心逻辑并不玄：把大气切成若干层，每层里按“吸收升温→折射率变化→相位扰动”算出一个屏，再用 FFT 角谱法让光束逐层往前传，最后就能看到典型的热晕畸变。适合做高能激光链路计算、自适应光学算法验证或者大气光学课程设计的工程师和研究生，拿到之后在 256 像素网格上就能直接验证光斑形态变化。
2. 相位屏的物理基础：热晕为什么不能当成普通湍流屏替代
2.1 热晕的“负透镜”机制与时间尺度
热晕的本质是大气对激光能量的吸收。光强 (I) 穿过吸收系数为 (\alpha) 的大气时，单位体积吸收的功率是 (\alpha I)，这个热量被气流以速度 (v) 带走。稳态下，沿风向的温度变化满足对流方程 (v \frac{\partial T}{\partial x} = \frac{\alpha I}{\rho C_p})，积分后得到横向温度分布：
[ \Delta T(x,y) = \frac{\alpha}{\rho C_p v} \int_{-\infty}^{x} I(x',y),dx' ]
其中 (\rho) 是大气密度，(C_p) 是定压比热。这个式子说明一个关键特性：热晕温度场不是随机的，它由光强分布和风速方向决定，是一个确定性的低阶相位结构。空气的折射率温度系数 (dn/dT) 约为 (-1 \times 10^{-6}) K(^{-1})，温度升高导致折射率下降，光束经过长度为 (dz) 的加热层后，相位扰动为：
[ \phi(x,y) = \frac{2\pi}{\lambda} \frac{dn}{dT} \Delta T(x,y) dz ]
这个相位是负的、凹的，等效于负透镜，光束被向外推。这里有一个新手最容易忽略的点：热晕相位屏和时间有关，风速决定了热量堆积的方向和大小。如果仿真里没有风速项，或者方向设置反了，光斑形态会完全不对。
2.2 谱密度模型里的取舍：Kolmogorov 谱和热晕屏的关系
普通大气湍流相位屏一般用 Kolmogorov 或 von Karman 折射率功率谱密度来生成，谱形式为 (\Phi_n(f) \propto C_n^2 (f^2 + 1/L_0^2)^{-11/6})，它描述的是折射率起伏的随机统计特性，产生的是高频随机相位，典型表现是光斑边缘出现随机破碎。热晕相位屏不一样，它的能量集中在低频——大尺度的散焦和偏向，而不是高频随机抖动。
实际仿真中，正确做法是两类屏叠加：确定性热晕屏负责“光斑整体形态畸变”，随机湍流屏负责“散斑细节和抖动”。我在拆这份资源时特别确认了这一点，程序里两个分量是分函数生成的，可以关闭其中任意一个来对比效果。如果不做区分、只用一个随机谱模型硬套热晕，跑出来的结果是高频噪声掩盖了大尺度变形，光斑看起来像碎裂的散斑，完全没有典型的新月形热晕特征。
| 模型 | 物理来源 | 相位特征 | 生成方式 |
|---|---|---|---|
| 普通湍流屏 | 温度/密度随机起伏 | 高频随机、PSD 满足负幂律 | FFT 谱反演 + 次谐波 |
| 热晕屏 | 激光加热大气的确定性过程 | 低频平滑、负透镜形状 | 对光强沿风积分 + 乘 dn/dT |
| 实际链路 | 两者同时存在 | 热晕主导形态，湍流叠加细节 | 分别生成后直接叠加相位 |
谱模型选择直接决定计算结果是否可信。只做热晕定性分析时，可以关掉随机湍流项，屏会更干净，便于观察形态；要做 Strehl 比或者光束质量 (\beta) 因子这类定量指标，就必须两者都开，且随机项要跑多次做系综平均。
3. 核心模块拆解：温度场积分与谱反演生成热晕相位屏
3.1 分层大气模型和整体传播链路
相位屏法把连续大气等效成一系列离散的薄层。每层厚度 (dz) 内做两件事：先让光束自由传播 (dz)，再乘一个相位屏 (\exp(i\phi))。传播部分用角谱法，在频域乘传播因子：
[ H(f_x,f_y) = \exp\left(-i\pi\lambda dz (f_x^2 + f_y^2)\right) ]
分层数一般取 2 到 10 段，太薄则计算量翻倍，太厚则单层内的温升假设不再成立。常见做法是等间距分层，每层 200 到 500 米。程序整体链路就是“高斯光束初始化→循环各层：角谱传播 + 相位屏叠加→输出焦面或近场光强”。
3.2 热晕相位主函数：光强积分与相位换算
% thermal_phase_screen.m % 输入光强分布，输出一帧确定性热晕相位屏 % I : 网格上的光强分布 [W/m^2] % dxy : 网格间距 [m] % v : 横向风速 [m/s]，定义沿 x 正方向 % alpha : 大气吸收系数 [1/m] % rho : 大气密度 [kg/m^3] % Cp : 定压比热 [J/(kg*K)] % dn_dT : 折射率温度系数 [1/K] % lambda : 波长 [m] % dz : 当前层厚度 [m] function pha = thermal_phase_screen(I, dxy, v, alpha, rho, Cp, dn_dT, lambda, dz) % 沿 x 方向积分光强，得到稳态温度场 dT = (alpha / (rho * Cp * v)) * cumtrapz(I, 2) * dxy; % 相位 = 波数 * dn/dT * 温度变化 * 层厚 k = 2 * pi / lambda; pha = k * dn_dT .* dT * dz; % 去掉整体活塞项，只保留空间起伏 pha = pha - mean(pha(:)); end这个函数核心就一行积分：cumtrapz(I, 2)沿矩阵第二维也就是 x 方向做累积梯形积分。为什么是从左往右积分？因为假设风从 x 负方向吹来，上游光强积分贡献小，越往下游热量堆积越多，温度逐渐升高，形成一边陡一边缓的新月形温度场。参数里alpha是体积吸收系数，单位必须用 1/m；很多资料给的是 dB/km，换算时要除以 (4.34 \times 10^3) 才能得到 1/m，这是最常见的单位坑。rho * Cp * v组合在一起是热对流项，风速越大温度升越小，热晕越弱，物理上直观。
3.3 高频随机项：谱反演生成湍流扰动屏
% random_phase_screen.m % 生成一帧高频湍流扰动相位屏（谱反演法） % N : 网格点数 % dx : 网格间距 [m] % Cn2 : 折射率结构常数 [m^-2/3] % L0 : 湍流外尺度 [m] % l0 : 湍流内尺度 [m] % lambda: 波长 [m] % dz : 层厚 [m] function phz = random_phase_screen(N, dx, Cn2, L0, l0, lambda, dz) df = 1 / (N * dx); fx = (-N/2 : N/2-1) * df; [FX, FY] = meshgrid(fx, fx); f = sqrt(FX.^2 + FY.^2); f(N/2+1, N/2+1) = 1e-10; % 避开零频除零 % von Karman 折射率功率谱密度 PSD = 0.033 * Cn2 .* (f.^2 + 1/L0^2).^(-11/6) .* exp(-(f*l0).^2); % 复高斯随机场激励，频域滤波后做逆傅里叶变换 g = randn(N) + 1i * randn(N); dn = real(ifft2(sqrt(PSD) .* g)); % 折射率起伏换算成相位起伏 phz = (2*pi/lambda) * dn * dz; % 方差标定：用理论相位方差对齐随机场幅度 theory_var = 2.895 * (2*pi/lambda)^2 * Cn2 * dz; % 近似平面波相位方差 phz = phz * sqrt(theory_var / var(phz(:))); end谱反演法的思路是在频域构造符合功率谱密度的随机场：随机复高斯阵g乘上sqrt(PSD)，逆傅里叶变换得到折射率起伏dn。乘以波数 (2\pi/\lambda) 和层厚 (dz) 把它换算成相位。代码最后一步做方差标定，这一步很重要，因为ifft2的归一化在不同 MATLAB 版本里一致，但谱幅度常因为网格分辨率不足偏低，直接生成会损失高频能量，用理论相位方差反推一个比例系数把幅度拉回来，是工程上稳妥的做法。
3.4 主程序整合：分层传播循环
% thermal_blooming_main.m % 主程序：两层大气，先传播后叠屏 N = 256; % 网格点数，越大高频越足 D = 0.5; % 屏物理尺寸 [m] dx = D / N; % 网格间距 [m] lambda = 1.06e-6; % 波长 [m]，典型固体激光 w0 = 0.03; % 高斯光束 1/e 半径 [m] I0_peak = 2e6; % 峰值光强 [W/m^2] rho = 1.225; % 空气密度 [kg/m^3] Cp = 1005; % 定压比热 [J/(kg*K)] alpha = 5e-5; % 吸收系数 [1/m] v = 3; % 横向风速 [m/s] dn_dT = -1e-6; % 折射率温度系数 [1/K] Cn2 = 1e-15; % 折射率结构常数 L0 = 50; l0 = 0.005; % 外内尺度 [m] layers = [250, 250]; % 两段各 250 m [X, Y] = meshgrid((0:N-1)*dx - D/2, (0:N-1)*dx - D/2); E = exp(-(X.^2 + Y.^2) / w0^2); % 归一化高斯场 for iz = 1:length(layers) dz = layers(iz); I = abs(E).^2 * I0_peak; pha_thermal = thermal_phase_screen(I, dx, v, alpha, rho, Cp, dn_dT, lambda, dz); pha_rand = random_phase_screen(N, dx, Cn2, L0, l0, lambda, dz); phz = pha_thermal + pha_rand; % 两个分量叠加 % 角谱传播 dz fx = (-N/2 : N/2-1) / (N*dx); [FX, FY] = meshgrid(fx, fx); H = exp(-1i * pi * lambda * dz * (FX.^2 + FY.^2)); E = ifft2(fft2(E) .* H) .* exp(1i * phz); % 先自由传播，再上相位屏 end I_out = abs(E).^2; imagesc((0:N-1)*dx - D/2, (0:N-1)*dx - D/2, I_out); axis image; colormap hot; title('热晕传播后光斑');主循环里顺序值得注意：先ifft2(fft2(E).*H)做真空传播，再乘exp(1i*phz)叠相位屏，这对应“先走一段路，再在该层位置受到相位扰动”的物理顺序。每层的光强I是用当前层的光斑算出来的，所以第二层热晕屏会根据第一层畸变后的光强重新积分，反应了热晕的自耦合特性。I0_peak和w0一起决定了初始热晕强度，这个值调大或调小时，光斑形态变化非常敏感。
4. 参数设置与验证：Strehl 比和光斑形态的自检方法
4.1 关键参数对照表
| 参数 | 符号 | 典型值 | 单位 | 设置不当的后果 |
|---|---|---|---|---|
| 网格数 | N | 256 / 512 | px | 太小则高频分量缺失，光斑边缘过于平滑 |
| 屏尺寸 | D | 0.3 ~ 0.6 | m | 小于光束半径 3 倍时边界效应严重 |
| 波长 | lambda | 1.06e-6 | m | 波长越短，同等温度场相位扰动越大 |
| 吸收系数 | alpha | 1e-5 ~ 1e-4 | 1/m | 取错数量级，热晕直接消失 |
| 风速 | v | 1 ~ 10 | m/s | 过小则温度积分过大，相位发散 |
| 层厚 | layers | 100 ~ 500 | m | 太厚时温升假设失真 |
| 峰值光强 | I0_peak | 1e6 ~ 1e7 | W/m^2 | 过低看不出热晕，过高导致数值溢出 |
最敏感的往往是alpha和I0_peak的组合。工程上可以先手算一个参考温度：取典型值 (\alpha=5\times10^{-5}) /m，(I_0=2\times10^6) W/m²，光束半径 (w_0=0.03) m，风速 3 m/s，则温度升大约 (\Delta T \approx \alpha I_0 w_0 / (\rho C_p v) \approx 0.008) K。这个量级看起来很小，但乘以波数 (2\pi/\lambda) 后相位就有几十弧度，足以造成强烈散焦。如果跑出来的屏相位最大值只有 0.1 rad 量级，说明参数数量级肯定不对。
4.2 传播验证：看形态和 Strehl 比
% 验证脚本片段：计算远场光斑半径和 Strehl 比 % 对输出近场做一次真空傅里叶变换得到远场分布 E_far = fftshift(fft2(fftshift(E))) / N; I_far = abs(E_far).^2; I_far = I_far / sum(I_far(:)) * sum(abs(E).^2); % 能量归一 % Strehl 比：实际峰值光强 / 无畸变衍射极限峰值光强 I_close = abs(ifft2(fft2(exp(-(X.^2+Y.^2)/w0^2)))).^2; % 无屏近场 strehl = max(I_out(:)) / max(I_close(:)); % 光斑质心偏移，热晕会整体向风下游偏移 [gx, gy] = meshgrid(1:N, 1:N); cx = sum(sum(gx .* I_out)) / sum(I_out(:)); fprintf('Strehl = %.3f, 光斑质心 x = %.2f px\n', strehl, cx);这段验证做三件事：一是远场强度归一，判断能量是否守恒；二是用 Strehl 比评估热晕导致的峰值下降；三是算光斑质心位置，热晕后光斑会整体往风下游偏，质心偏移是热晕的典型特征。无畸变条件下 Strehl 接近 1，强热晕条件下会掉到 0.1 以下，这本身就是最直观的判据。
4.3 从结果反推参数是否合理
如果 Strehl 低于 0.05 且光斑完全碎成多个小块，先不要急着怪热晕模型，很可能是随机湍流项太强或者网格分辨率不够。关闭pha_rand只留热晕屏，此时光斑应该是连续、圆滑的新月形；如果仍然碎裂，问题在传播步长或网格尺寸。我一般把这一招当成第一自检手段：分层关掉组分，看哪部分导致形态异常。另一条经验是观察光斑拖尾方向——热晕的形态是迎风面陡峭、背风面拖尾，如果拖尾方向反了，把风速定义成负值即可。
5. 避坑与排查：热晕相位屏仿真的五个翻车点
5.1 相位屏全是高频雪花点，看不出热晕轮廓
现象：生成的屏展开后全是密密麻麻的颗粒，完全看不到平滑的新月形温度分布。
原因：只调用了随机谱反演函数，没有叠加热晕温度积分项；或者热晕函数里alpha单位给错，导致确定性分量小到可以忽略。
解决：检查主循环是否同时调用了thermal_phase_screen和random_phase_screen并把两者相加。先单独输出pha_thermal画图确认它是平滑的大尺度分布，再叠加随机项。还要确认alpha的单位是 1/m，不是 dB/km。
5.2 Strehl 比高得离谱，热晕好像不存在
现象：初始光强和吸收系数都设了，但跑出来 Strehl 接近 0.98，光斑仍然是完美高斯圆斑。
原因：吸收系数、峰值光强或风速三个参数组合以后，温度升太小，相位扰动不到 1 rad。
解决：先用前面给的公式手算 (\Delta T) 和相位幅度，确认数量级。相位应该有几十 rad 才是强热晕，几 rad 是弱热晕，0.01 rad 就是没热晕。一般先调高alpha或I0_peak一两个数量级看响应是否出现。
5.3 光斑不对称，但拖尾方向反了
现象：光斑确实变成了新月形，但拖尾方向和实验观测相反，且质心偏移方向与设定风向相反。
原因：cumtrapz(I, 2)沿矩阵第二维积分，默认从左往右积累，但主程序里风速方向如果没有和网格 x 坐标统一起来，就会出现符号错误。
解决：把风速方向明确写死在代码注释里，约定v > 0表示风沿 x 正方向，光斑拖尾应在 x 正方向。如果反了，把v改成负值，或者把cumtrapz(I,2)改成对翻转后的矩阵积分再翻回来。
5.4 屏边缘出现亮条纹和鬼像
现象：光斑四角或边缘出现规则亮条纹，像镜像复制。
原因：FFT 假设信号周期延拓。光斑半径接近屏边界时，边界切断产生衍射，角谱传播时边界反射回来形成鬼像。
解决：屏尺寸至少取光束半径的 5 倍以上，让光斑峰值远离边缘 20 个像素以上。更稳妥的办法是给光场加边缘窗函数，比如余弦窗，把边缘光强压到接近零。还有一种做法是扩大网格但只对中心区域做温度积分，边缘区域不叠屏，避免周期延拓。
5.5 同一参数跑两次，Strehl 差出五倍
现象：不改任何参数，只换一次随机种子，输出 Strehl 从 0.05 变成 0.5。
原因：热晕屏是确定性的，但随机湍流分量单帧波动很大，尤其是低频大尺度成分对 Strehl 影响极强。单帧结果形同抽签。
解决：固定rng种子做回归对比，正式评估时跑 50 到 100 帧取平均。热晕的确定性分量在多次平均后保持新月形，随机项会被平均掉，平均后的形态才是稳态热晕特征。交付给别人时，把种子写进参数文件，否则别人无法复现你的结果。
6. 进阶用法：批量生成相位屏并接进光学链路
6.1 批量生成与统计平均
一帧屏只能说明单次随机相位实现，热晕仿真要出统计结论，批量生成是必经步骤。常见做法是把热晕屏和随机屏的生成放进循环，把 50 到 100 帧结果存成三维数组，再统一做平均和方差统计：
% batch_screens.m rng(2025); % 固定种子，保证可以复现 numFrames = 80; screenArray = zeros(N, N, numFrames); strehl_log = zeros(numFrames, 1); for k = 1:numFrames pha_thermal = thermal_phase_screen(I, dx, v, alpha, rho, Cp, dn_dT, lambda, dz); pha_rand = random_phase_screen(N, dx, Cn2, L0, l0, lambda, dz); phz = pha_thermal + pha_rand; screenArray(:, :, k) = phz; E = ifft2(fft2(E0) .* H) .* exp(1i * phz); I_out = abs(E).^2; strehl_log(k) = max(I_out(:)) / max(I_close(:)); end mean_screen = mean(screenArray, 3); std_screen = std(screenArray, 0, 3); save('thermal_screens_batch.mat', 'screenArray', 'mean_screen', 'std_screen', ... 'strehl_log', 'dx', 'lambda', 'N', 'v', 'alpha', '-v7.3');批量保存时把参数结构体一起存进去，这一点很重要。我接手过别人的屏文件，光一个.mat没有任何参数记录，网格间距、层厚、风速全是未知数，根本没法复用。现在习惯是把N, dx, lambda, layers, v, alpha, rng种子全部打包进同一个文件。
6.2 量化评估与交付建议
平均屏mean_screen展示的是确定性热晕形态；方差屏std_screen展示的是随机湍流的不确定性范围。对自适应光学工程师来说，这两个量比单帧光斑图像有用得多：平均形态决定波前校正器需要多大的低频行程，方差分布决定高频校正器的动态范围。如果要做光束质量评估，建议同时输出 63% 环绕能量直径、质心偏移量和 Strehl 比三个量，并给出这批屏对应的均值和标准差，粗糙交付一个光斑图意义不大。
每个项目我都强烈建议保留一次无热晕的基准屏做对照。同一条传播链路，关掉pha_thermal只保留湍流跑一组，再开启热晕跑一组，两者差值就是纯热晕贡献。这个差值在写报告和判断系统余量时，比绝对 Strehl 值更有说服力，也更容易向非光学背景的人解释清楚。
以前我在给一套光束控制链路供屏时，拿单帧屏就跑闭环算法，Strehl 忽高忽低，折腾了两天查不出原因，后来才意识到是忘了做系综平均。从那以后我每次交付屏文件都会强制走一遍固定种子、批量生成、参数打包三步流程，确认平均圈数和相位方差都记录在案才发出去。希望帮到你。
本文还有配套的精品资源，点击获取
