MATLAB热晕相位屏仿真：物理模型、代码实现与避坑指南

简介这份资源是一套基于MATLAB的热晕相位屏仿真程序面向光学工程、激光传输与望远镜系统方向的研究人员和学生用于定量分析大气温度不均匀引起的光波前畸变问题。压缩包共13个文件约541KB以7个m脚本为核心配合4个bmp图像、1个asv备份与1个txt说明脚本覆盖相位屏生成、傅里叶变换、光学传递函数与点扩散函数计算及成像分析等关键环节图像文件可用于对照仿真输出结果。程序通过调整大气条件、温度梯度等输入参数帮助读者观察不同场景下热晕对成像质量与分辨率的影响并理解高斯分布、泰勒级数展开在相位屏建模中的作用。已有254人学习下载适合作为高能激光与光学系统热晕效应研究的入门与验证工具也可为后续补偿技术探索提供可复用的仿真基础。1. 热晕相位屏仿真从一束光在介质里的“自毁”说起高功率激光在吸收介质中传输时会经历一个反直觉的过程光越强自己把自己越推越散。介质吸收少量光能后折射率随温度升高而改变形成类似负透镜的相位分布光束中心被“推”向外围远场光斑反而膨胀、畸变。这就是热晕效应。做激光定向能、高能激光大气传输、甚至某些精密干涉测量的同行迟早要面对它。而“相位屏”是把连续介质拆成若干薄片、每片只施加一个相位扰动的数值手段MATLAB 因为矩阵操作顺手、可视化直接成了很多人搭热晕仿真程序的第一选择。这篇笔记就围绕一个典型的热晕相位屏仿真程序把物理模型、代码骨架、参数怎么调、哪里容易翻车一次讲透。2. 热晕相位屏的物理模型与数值离散为什么不能只靠一个相位公式热晕仿真的核心矛盾在于光场强度分布决定介质温度分布温度分布又反过来改变折射率折射率再作用于光场相位。这是一个闭环。如果只写一个静态相位公式相当于默认“光场不变、热晕固定”那只能算线性近似对高功率长距离场景误差很大。真正要复现实验趋势必须做交替迭代光场传播一步算吸收生热更新折射率再施加相位再传播。相位屏方法就是把这段连续过程切成 N 段每段内认为介质参数均匀只把累积效应折算成一个薄相位屏。2.1 从热传导方程到折射率扰动三个必须写进代码的物理关系介质吸收的光功率密度 $Q(r,z)\alpha I(r,z)$其中 $\alpha$ 是吸收系数$I$ 是当地光强。温度场 $T$ 满足热传导方程 $\rho c_p \partial T/\partial t \kappa \nabla^2 T Q$。在准稳态、忽略横向热扩散的常见简化下温升近似正比于累积吸收能量。折射率变化 $\Delta n (\mathrm{d}n/\mathrm{d}T) \Delta T$对多数液体和气体$\mathrm{d}n/\mathrm{d}T$ 为负所以中心温升最高、折射率最低等效于负透镜。相位扰动 $\Delta \phi k \Delta n \Delta z$$k2\pi/\lambda$。这三步关系决定了代码里至少要有三个数组光强矩阵、温度增量矩阵、相位屏矩阵。很多初学者只写了相位屏却忘了温度场需要沿传播方向累积导致热晕强度不随距离增长仿真结果永远偏弱。2.2 相位屏该放多少片、每片多厚采样规则与经验取值相位屏数量 $N_{screen}$ 不是越多越好。每片厚度 $\Delta z L/N_{screen}$$L$ 是总传输距离。如果 $\Delta z$ 太大单次相位扰动过强会引入非物理的衍射伪影如果太小计算量线性增长而且相邻屏之间光场变化微弱浪费算力。经验做法是让每片引起的最大相位变化不超过 $\pi$即 $\max|\Delta \phi| \le \pi$。对于典型高能激光大气传输参数$N_{screen}$ 取 20 到 50 之间比较稳妥。另外相位屏的横向网格必须与光场网格一致且满足采样定理网格间距 $\Delta x \le \lambda z/(2D)$$D$ 是光束直径。如果网格太粗高频相位信息丢失热晕畸变会显得过于平滑远场环状结构出不来。2.3 吸收系数与热光系数两个最容易被低估的输入参数吸收系数 $\alpha$ 的单位是 $\mathrm{m}^{-1}$但很多文献给的是 $\mathrm{cm}^{-1}$ 或 $\mathrm{dB/km}$换算时差两个数量级热晕强度直接差百倍。热光系数 $\mathrm{d}n/\mathrm{d}T$ 对空气约 $-1\times10^{-6}\ \mathrm{K}^{-1}$对某些液体可到 $-4\times10^{-4}\ \mathrm{K}^{-1}$差两个数量级。这两个参数如果从不同来源拼凑仿真结果可能完全不可信。我一般会先做量纲检查算一个稳态温升 $\Delta T \approx \alpha I L/(\rho c_p v)$再乘 $\mathrm{d}n/\mathrm{d}T$ 和 $kL$看总相位是否在合理范围通常几弧度到几十弧度。如果算出来是 0.001 弧度那热晕根本不可见要么参数错了要么场景不对。3. 用 MATLAB 搭一个可跑的热晕相位屏程序从网格到迭代循环下面给出一套最小可运行骨架包含光场初始化、分段传播、热致相位更新和远场评价。代码不依赖任何工具箱纯矩阵运算方便你直接改成自己的场景。注意这里用分步傅里叶法做衍射传播热相位用简化准稳态模型适合观察趋势和做参数扫描不适合替代严谨的流体-光学耦合求解。3.1 初始化光场与网格避免 aliasing 的两种实用手段% 热晕相位屏仿真 - 初始化 lambda 1.06e-6; % 波长单位 m k 2*pi/lambda; D 0.1; % 光束直径单位 m L 1000; % 总传输距离单位 m Nx 512; % 网格点数 dx 1.5e-3; % 网格间距单位 m x (-Nx/2:Nx/2-1)*dx; [X, Y] meshgrid(x, x); R2 X.^2 Y.^2; % 高斯光束初始场 w0 D/2; E0 exp(-R2/w0^2); I0 abs(E0).^2; P0 sum(I0(:))*dx^2; % 总功率用于后续标定 % 频域坐标用于分步傅里叶 df 1/(Nx*dx); fx (-Nx/2:Nx/2-1)*df; [FX, FY] meshgrid(fx, fx); K2 (2*pi)^2*(FX.^2 FY.^2);这段代码做了三件事定义物理单位和网格、生成高斯光场、准备频域坐标。关键参数是dx和Nx。dx必须足够小让光束边缘至少落在 3 到 4 个网格内否则截断会引入高频振荡。如果dx太大光场在边界突然降到零相当于一个硬边光阑远场会出现强烈的衍射环和热晕混在一起分不清。另一个实用手段是加超高斯吸收边界比如在边缘 10% 区域乘一个平滑窗函数把能量缓慢衰减到零而不是直接截断。我一般会检查初始光场的频谱如果高频分量在 Nyquist 频率附近还有明显幅度说明网格不够细需要减小dx或增大Nx。3.2 分步傅里叶传播与热相位施加一个循环里的四个动作Nscreen 30; % 相位屏数量 dz L/Nscreen; % 每段长度 alpha 1e-4; % 吸收系数1/m dndT -1e-6; % 热光系数1/K rho_cp 1.2e3; % 密度乘比热J/(m^3 K) E E0; Tacc zeros(Nx); % 累积温升 for n 1:Nscreen % 1) 频域传播 dz/2 E ifft2(fft2(E).*exp(-1i*K2*dz/4/k)); % 2) 计算吸收生热并累积温升 I abs(E).^2; dT alpha * I * dz / rho_cp; Tacc Tacc dT; % 3) 施加热致相位屏 dphi k * dndT * Tacc * dz; E E .* exp(1i * dphi); % 4) 再传播 dz/2 E ifft2(fft2(E).*exp(-1i*K2*dz/4/k)); end % 远场评价 E_far fftshift(fft2(E)); I_far abs(E_far).^2;循环里四个动作对应物理过程衍射传播、吸收生热、相位调制、再传播。这里用了对称分步傅里叶精度比单步高。注意dphi的计算里乘了dz因为Tacc已经是累积温升而相位是沿传播方向累积的。如果漏掉这个dz热晕强度会随屏数线性增长而不是平方增长趋势就错了。另一个细节是Tacc没有考虑横向热扩散。在短脉冲或低重复频率场景下热扩散时间远大于脉冲间隔这个近似成立但在连续波长时间传输中横向扩散会显著降低中心温升热晕饱和这时候需要加扩散项或者用更复杂的模型。我一般会先跑这个简化版看趋势如果远场光斑膨胀到 2 到 3 倍就停因为再往后简化模型已经不可信了。3.3 远场光斑评价用环围能量和峰值强度两个指标% 环围能量曲线 P_total sum(I_far(:)); P_enc zeros(100,1); r_max Nx*dx/2; for m 1:100 r_th m/100 * r_max; mask (X.^2 Y.^2) r_th^2; P_enc(m) sum(I_far(mask))/P_total; end % 峰值强度比 Strehl max(I_far(:)) / max(I_far(:)) ; % 需要与无热晕对比实际使用时应该跑两次一次dndT0作为基准一次正常参数。然后比较峰值强度比和环围能量曲线。热晕的典型特征是峰值下降、能量向外环转移。如果只看远场图像容易把衍射环误认为热晕。用环围能量曲线更定量无热晕时 80% 能量在某个半径内热晕后同样半径内只剩 40%这个差异很直观。注意Strehl那行需要和无热晕结果对比不能自己除自己。我一般会把两次结果存成结构体最后画在一张图上。4. 参数扫描与结果验证怎么判断你的热晕仿真没跑偏仿真跑通只是第一步更常见的问题是结果看起来“有热晕”但强度不对、趋势不对。这一章给出一套自检流程帮你判断参数是否合理、代码有没有隐藏错误。4.1 用解析近似做量纲校验三个快速检查点第一个检查点总相位量级。用 $\Delta \phi_{max} \approx k |\mathrm{d}n/\mathrm{d}T| \alpha I_0 L^2/(2 \rho c_p)$ 估算。如果算出来小于 0.1 弧度热晕可忽略如果在 1 到 10 弧度远场会有明显畸变如果大于 100 弧度可能出现相位包裹需要检查网格是否够细。第二个检查点热晕特征时间 $t_{th} \rho c_p/(\alpha I_0 |\mathrm{d}n/\mathrm{d}T| k)$如果仿真时间远小于它热晕来不及建立。第三个检查点菲涅尔数 $N_F D^2/(\lambda L)$它决定衍射和热晕的相对强弱。如果 $N_F$ 很小衍射本身就把光斑扩得很大热晕效应被掩盖。我一般会先算这三个数再决定要不要跑仿真。4.2 网格收敛性测试改变 Nx 和 Nscreen 看结果是否稳定固定其他参数把Nx从 256 加到 512 再加到 1024看远场峰值强度比是否收敛。如果 256 和 512 差很多说明网格不够。同样把Nscreen从 10 加到 30 再加到 50看环围能量曲线是否重合。如果 10 和 30 差很多说明分段太粗。收敛性测试是判断仿真可信度的硬标准比任何理论分析都直接。我一般会跑一个 2x2 表格Nx 取 256/512Nscreen 取 15/30四个结果画在一起如果两条曲线基本重合就可以用较小那组做后续扫描省时间。4.3 与文献趋势对比热晕畸变随功率和距离的典型曲线在简化模型下热晕引起的峰值强度下降大致随功率平方或距离三次方增长直到饱和。如果你扫描功率发现峰值强度线性下降那可能漏掉了温度累积的非线性如果随距离指数下降那可能相位计算里多乘了不该乘的因子。常见做法是画双对数图横轴功率或距离纵轴峰值强度比看斜率是否在 2 到 3 之间。如果斜率是 1回去检查Tacc的累积方式如果斜率是 5 以上检查是否有数值不稳定导致能量非物理扩散。5. 避坑与排查热晕相位屏仿真里最容易翻车的五件事这一章记录我在搭这类程序时踩过的坑每条按现象、原因、解决写。你如果遇到类似症状可以直接对照。5.1 远场光斑始终对称热晕环出不来现象跑完仿真远场光斑只是变大但没有出现典型的热晕环或不对称畸变。原因最常见的是dndT符号搞反了。如果设成正数中心折射率升高等效正透镜光束会聚焦而不是发散远场反而变小。另一个原因是Tacc没有沿传播方向累积每段都用同一个温升热晕强度不随距离增长。解决检查dndT是否为负检查Tacc Tacc dT是否在循环内。还有一个隐蔽原因光强I没有归一化导致dT量级错误热晕太弱看不出来。用总功率P0标定一下确保峰值光强在合理范围。5.2 增加相位屏数量后结果反而变差现象Nscreen从 20 加到 50远场峰值强度比不升反降甚至出现高频振荡。原因每片相位扰动过强时分步傅里叶的对称拆分假设失效相位屏之间的衍射不能正确抵消。另外如果dz太小频域传播因子exp(-1i*K2*dz/4/k)在数值上接近 1但浮点误差累积可能导致相位噪声。解决先检查单屏最大相位是否超过 $\pi$如果超过要么减少Nscreen要么降低功率或吸收系数。如果单屏相位很小但结果仍振荡检查频域坐标K2是否计算正确特别是df和dx的对应关系。我一般会固定Nscreen30先调其他参数最后再微调屏数。5.3 功率扫描时热晕阈值不明显现象改变输入功率远场峰值强度比缓慢下降没有明显的阈值或拐点。原因简化模型里热晕是渐变的没有真实实验中的对流阈值。如果你期待看到突变需要加入横向对流项或时间依赖。另一个原因是吸收系数太小热晕效应被衍射掩盖。解决先算菲涅尔数和热晕特征时间确认参数落在热晕主导区。如果确实需要阈值行为可以引入一个经验对流速度把Tacc乘以一个衰减因子模拟热对流带走的能量。但要注意这已经超出简化模型范围结果只能定性看。5.4 内存溢出或跑一次要几十分钟现象Nx1024、Nscreen100时MATLAB 报内存不足或者循环慢到无法接受。原因每个屏都做两次fft2和ifft2复杂度 $O(N_x^2 \log N_x)$屏数线性增加。另外Tacc、E、I都是复数或双精度矩阵1024x1024 下每个约 16 MB多个变量叠加容易爆。解决用single精度存光场和温度场内存减半对趋势观察足够。把fft2换成fft2(..., symmetric)如果光场有对称性但一般不建议。更实用的做法是先用 256 网格做参数扫描找到感兴趣的区域后再用 512 或 1024 精算。我一般会在循环里加if mod(n,10)0, fprintf(...), end看进度避免盲等。5.5 远场评价时环围能量曲线不单调现象环围能量曲线随半径增加出现下降或平台。原因远场光斑没有居中fftshift后峰值不在矩阵中心导致圆形掩膜统计错误。另一个原因是远场网格的物理尺寸没有正确对应r_max算错。解决在fftshift后找到峰值位置把矩阵滚动到中心或者直接用regionprops类方法但不用工具箱。更简单的是检查sum(I_far(:))是否等于sum(I0(:))如果不等说明传播过程中能量不守恒可能是边界截断或相位屏吸收。我一般会在每次传播后算总能量画一条能量随距离的曲线如果下降超过 1%就要检查边界处理。6. 进阶技巧把热晕仿真从“能跑”推到“敢用”最后一章聊几个让仿真结果更可信、更接近实际工程判断的技巧。这些不是必须但如果你要把热晕相位屏程序用于方案对比或参数优化它们能帮你少走弯路。第一个技巧是引入时间步进。简化模型是准稳态的但实际热晕有建立过程。你可以把Tacc写成时间导数形式Tacc Tacc (alpha*I/rho_cp - Tacc/tau)*dt其中tau是热弛豫时间。这样能模拟脉冲串或扫描光束下的热晕累积和消退。dt取热晕特征时间的 1/10 到 1/20tau根据对流速度估算。这个改动不大但能让你的程序从“静态相位屏”升级到“动态热晕”对重复频率激光场景特别有用。第二个技巧是用多物理场耦合验证。如果你有条件把简化模型算出的温升分布导出和某个商业流体软件或实验测量的折射率分布对比。不需要完全一致看趋势和量级是否对得上。我一般会取中心线上的温升曲线简化模型通常高估中心温升因为忽略横向扩散。如果高估在 2 倍以内可以接受如果高估 10 倍说明吸收系数或热光系数需要重新标定。第三个技巧是参数敏感性分析。不要只跑一组参数用拉丁超立方或简单网格扫描看哪些参数对远场峰值影响最大。通常吸收系数和热光系数是主导网格间距和屏数影响较小但也不能忽略。把敏感性结果画成条形图你就能快速判断哪个参数需要精确测量哪个可以粗略估计。这个习惯让我在几次方案对比中避免了把时间浪费在次要参数上。最后一个习惯每次跑完仿真先问自己三个问题——能量守恒吗网格收敛吗量级合理吗这三个问题只要有一个答不上来结果就不能用。热晕仿真没有后悔药前期多花十分钟检查比后期推翻重来划算得多。希望帮到你。本文还有配套的精品资源点击获取