简介这是一份MATLAB气动光学仿真程序包主题涵盖高斯光束、涡旋光束与环形光束在大气湍流等条件下的传输特性适合光学工程、大气科学及相关课题的本科生、研究生和科研人员用于理论验证与数值实验。压缩包内共有3个文件全部为m脚本包括光束传播主程序、屏幕系列分析和GIF动态可视化程序整包仅2KB代码结构简洁便于快速运行与二次开发。目前已有299人学习下载。该程序可在多个屏幕位置记录光强分布并生成动态图直观展示不同光束在大气中传输时的扩展、畸变与涡旋稳定过程同时为瑞利散射、大气吸收及光涡旋演化等关键问题提供可复现的仿真框架可作为课程设计或科研预研的起步模板。1. 气动光学效应光束穿过大气时到底发生了什么一束原本发散角很小的激光在几公里外用CCD接收光斑不再是一个漂亮的高斯斑而是出现随机抖动、破碎甚至分裂成多个小块。很多人第一反应是大气散射或吸收但真正的元凶往往是折射率随机起伏带来的相位畸变这就是气动光学的核心问题。大气湍流像一堆大小不同的透镜乱放在光路上每时每刻都在改变光束的波前。要定量评估这种影响光做静态光学设计不够需要对随机过程做统计模拟。接下来要拆的这套MATLAB程序就是用来做这类模拟的构造高斯、涡旋、环形光束让它们穿过湍流相位屏观察光斑演化并统计闪烁指数、质心漂移等参数。适合激光通信链路预算、自适应光学系统仿真、大气光学实验预研的人参考。2. 湍流相位屏与分步傅里叶传播从Kolmogorov谱到可执行的模型气动光学效应的物理根源是大气折射率随机变化而描述这种随机性的标准模型是Kolmogorov湍流理论。工程上不用直接解Navier-Stokes而是用功率谱密度来构造等效相位扰动。2.1 折射率功率谱哪些空间尺度决定光束畸变大气折射率起伏的空间统计特性常用Von Karman谱表示[ \Phi_n(\kappa) 0.033 C_n^2 \frac{\exp(-\kappa^2/\kappa_m^2)}{(\kappa^2\kappa_0^2)^{11/6}} ]其中 (\kappa) 是空间波数(C_n^2) 是折射率结构常数典型值从 (10^{-17})弱湍流到 (10^{-13})强湍流不等。(\kappa_0 2\pi/L_0) 是外尺度对应的波数(L_0) 通常取 20100 米(\kappa_m 5.92/l_0)(l_0) 是内尺度通常在毫米量级。相位屏生成的核心思路是把一个随机高斯场通过频域滤波使其幅度满足上述谱再变换到空间域从而得到一组相位样本。我的习惯是生成大量独立相位屏每次传播用其中一张最后对结果做系综平均。只跑一次得到的光斑形态没有统计意义。function phz phase_screen_km(n, dx, Cn2, Lz, L0, l0, seed) % n: 网格点数; dx: 网格间距; Lz: 相位屏等效厚度; Cn2: 湍流强度 if nargin 6, rng(seed); end k (0:n-1) - n/2; % 频率坐标中心为0 [kx, ky] meshgrid(k, k); kr2 (kx.^2 ky.^2) / (n*dx)^2; % 空间波数平方 k0 2*pi/L0; km 5.92/l0; p 0.023 * Cn2 * Lz * (4*pi^2) * ... % 常数项换算 (1 kr2/k0^2).^(-11/6) .* exp(-kr2/km^2); phz real(ifft2( sqrt(p) .* ... (randn(n) 1i*randn(n)) )) * n^2; % 逆变换后的相位 phz phz - mean(phz(:)); end这段代码生成一张尺寸为 (n \times n) 的相位屏。关键参数是Lz相位屏代表的传播厚度如果模拟传播距离 1000 米并分成 10 步则每步Lz100。n*dx是物理尺寸必须足够大否则相位屏边缘的统计特性会失真。randn(n)产生复高斯随机场乘上功率谱的平方根相当于在频域整形。2.2 分步傅里叶法为什么不能直接乘一个相位屏如果光在真空和湍流中交替传播可以在一个薄屏处叠加所有扰动不行。大气湍流是连续分布的三维不均匀体但工程上经常用“相位屏近似”简化将传播路径切成若干段每段长度足够短使得湍流效应可以集中在该段中心的薄屏上。段间用真空衍射传播这就是分步傅里叶法。每一步的操作是在真空中传播距离 (L_z)用角谱法或菲涅尔衍射积分。将光场乘以相位屏 (\exp(i\theta))其中 (\theta) 是相位屏产生的相位。MATLAB里最常用的真空传播是角谱法function u angular_spectrum_prop(u0, dx, lambda, z) % u0: 输入光场; dx: 采样间距; lambda: 波长; z: 传播距离 [n, ~] size(u0); fx (-n/2 : n/2-1) / (n*dx); [FX, FY] meshgrid(fx, fx); H exp(1i * (2*pi/lambda) * z .* ... sqrt(1 - (lambda*FX).^2 - (lambda*FY).^2)); U fftshift(fft2(u0)); u ifft2(ifftshift(U .* H)); end这段代码用角谱传递函数做衍射计算比直接调用数值积分快得多。注意sqrt里的值若出现负数说明高频分量超出传播条件需要降低dx或增大网格尺寸。fftshift和ifftshift的配对顺序错了会导致光场翻转这是最常见的低级错误。传播路径分段时每段的距离要满足相位屏间的自由传播不会引入采样混叠常见做法是使每段传播的菲涅尔数 (F w^2/(\lambda z)) 大于 1其中 (w) 是光束半径。分段数多则每段湍流弱统计结果更平缓但计算量线性增长。3. 高斯、涡旋与环形光束的MATLAB构造与传输脚本exercise_beam_propagation_95.m 这类脚本核心逻辑并不复杂生成初始光场循环执行“真空传播 相位屏”最后记录输出光强。关键是不同光束的初始场表达式以及统计口径。3.1 三种典型光束的初始场高斯光束的基础形式是[ u_0(r) \exp\left(-\frac{r^2}{w_0^2}\right) \exp\left(-i \frac{k r^2}{2R}\right) ]其中 (w_0) 是束腰半径(R) 是波前曲率半径聚焦情形下 (R) 取正值发散时为负。在MATLAB里构造n 512; dx 0.002; wavelength 1.55e-6; [x, y] meshgrid((-n/2:n/2-1)*dx); r2 x.^2 y.^2; w0 0.02; R -500; k 2*pi/wavelength; u0_gauss exp(-r2/w0^2) .* exp(-1i*k*r2/(2*R));R -500表示入射波前呈发散状相当于光束在传播 500 米处虚聚焦。如果改成R 500光束先聚焦再发散光斑演化过程完全不同。涡旋光束带有螺旋相位结构最简单的模型是phi atan2(y,x); % 方位角 l 2; % 拓扑荷数 u0_vortex u0_gauss .* exp(1i*l*phi);拓扑荷 (l) 为 2 时相位沿方位角变化 4π中心处相位奇异性导致强度为零。涡旋光束在湍流中会发生模式串扰拓扑荷越高抵抗湍流导致的强度衰减能力也不一样这正好适合做对比。环形光束更常用空心高斯或拉盖尔-高斯模式。空心高斯的一种简单实现是ring_radius 0.03; ring_width 0.005; u0_ring exp(-(sqrt(r2)-ring_radius).^2/ring_width^2);上式构造一个中心暗、半径约 3 cm 的环形亮斑。注意sqrt(r2)会在坐标原点产生一个“尖角”如果网格分辨率不够环形光束的中心暗斑会变得不均匀。3.2 传播循环与数据保存把初始场、相位屏、传播距离封装进主循环z_total 2000; dz 200; num_steps z_total / dz; u u0; for s 1:num_steps phz phase_screen_km(n, dx, Cn2, dz, L0, l0, s); u u .* exp(1i*phz); % 湍流相位扰动 u angular_spectrum_prop(u, dx, wavelength, dz); % 真空衍射 I abs(u).^2; record_intensity(s) {I}; % 保存每步光强 endrecord_intensity用元胞数组保存每一步的光强方便后续生成动态图。如果内存紧张可以只保存指定的几个“屏幕”位置。3.3 屏幕序列采样何时保存光斑图像exercise_screen_series_95.m 这个脚本的关键不是光场传播本身而是采样的时机和位置。比如每 100 米保存一次共 20 帧。但要注意保存稀疏帧会丢失光斑抖动的中间细节保存太密则GIF文件会非常大。我一般先跑一次粗采样例如 50 米一步看到光斑特征后再针对感兴趣区间做细采样。4. 闪烁指数、质心漂移与GIF动态可视化光强分布随机起伏单看一帧很难评价。必须用统计量描述再把多帧压缩成GIF直观观察光斑抖动趋势。4.1 统计量计算闪烁指数和质心漂移闪烁指数定义为[ \sigma_I^2 \frac{\langle I^2 \rangle}{\langle I \rangle^2} - 1 ]质心漂移则用光强的能量加权位置function [scint, centroid] beam_stats(I, dx) % 输入光强矩阵和采样间距 total sum(I(:)); I_norm I / total; sx (1:size(I,2)) * dx; % x方向网格坐标 sy (1:size(I,1)) * dx; cx sum(sx .* sum(I_norm,1)); % 光斑质心x cy sum(sy .* sum(I_norm,2)); centroid [cx, cy]; scint mean(I(:).^2) / mean(I(:)).^2 - 1; end注意这个质心计算没有减去中心坐标偏移你需要在调用时自行减去初始质心才能得到抖动偏移量。4.2 不同光束在湍流中的表现对比在 (C_n^210^{-14})距离 2 km波长 1.55 μm 的条件下我跑过一组对比大致结果如下光束类型闪烁指数质心漂移 RMS中心强度保持率高斯光束0.422.3 mm0.28涡旋光束 (l2)0.512.7 mm0.21涡旋光束 (l4)0.573.1 mm0.17环形空心光束0.472.5 mm0.24这个结果符合一般规律拓扑荷越高初始相位变化越剧烈湍流引起的相位扰动和光强闪烁越强。但并不是说涡旋光束一定更差在无湍流时涡旋光束能保持螺旋相位结构只是统计结果上对湍流更敏感。4.3 把光斑序列压缩成GIFgif.m 里最核心的是用imwrite生成GIF文件。注意MATLAB原生GIF动画是通过循环写入实现的figure(Visible,off); for i 1:num_steps imagesc(abs(U{i}).^2); axis image; colormap hot; clim([0, 1]); % 固定颜色范围避免闪烁 set(gca,XTick,[]); set(gca,YTick,[]); frame getframe(gcf); [A, map] rgb2ind(frame.cdata, 256); if i 1 imwrite(A, map, beam_propagation.gif, gif, ... LoopCount, inf, DelayTime, 0.1); else imwrite(A, map, beam_propagation.gif, gif, ... WriteMode, append, DelayTime, 0.1); end end固定clim是最容易忽略的细节。如果每帧都用独立的 colorbar 范围光斑强度涨落会被掩盖看起来整个屏幕都在“闪”而不是光斑在抖。延迟时间 0.1 秒表示10 FPS太长会卡顿太短看不清抖动轨迹。生成的GIF可以直接嵌入演示文稿里比贴一堆光斑图直观得多。5. 参数调优与常见坑从代码能跑到结果可信很多拿到这套程序的人第一反应是“跑起来了但光斑一点不抖”。多数不是代码bug而是参数设置没有处于湍流有效作用区间。5.1 关键参数的量级估计影响结果是否“看得出湍流”的两个核心参数是传播距离 (z) 和 (C_n^2)。弱湍流下Rytov方差为[ \sigma_R^2 1.23 C_n^2 k^{7/6} z^{11/6} ]闪烁指数近似等于 (\sigma_R^2)当 (\sigma_R^2 0.3) 时。如果 (C_n^210^{-16})、(z500) m(\sigma_R^2) 可能在 0.01 以下光斑基本不抖。要看到明显效应至少让 (\sigma_R^2) 在 0.1 以上。我自己常用的“保守起效”组合是(C_n^210^{-14})、距离 1000 米以上。另一个坑是相位屏的分辨率。相位屏的频域采样间隔是 (1/(n dx))如果这个值大于湍流内尺度对应的空间频率小而强的相位涡旋会被平滑掉表现为光斑只有整体漂移而没有破碎细节。以 (l_05) mm 为例最大空间频率需要达到 (1/l_0200) /m。取 (n512)则 (dx) 必须小于 (1/2000.005) m也就是网格物理尺寸小于 2.56 m。如果光束初始半径是 2 cm这个尺寸完全够用但如果模拟 10 cm 宽的光束网格就要到 1024 或 2048。5.2 相位屏统计一致性的验证生成相位屏后先别急着跑全传播。可以单独验证相位屏的统计性质计算结构函数 (D_\phi(r))并与理论值比较。用100张相位屏做系综平均得到的结果应当趋于直线在Kolmogorov区间。phz_all zeros(n,n,100); for s 1:100 phz_all(:,:,s) phase_screen_km(n, dx, Cn2, dz, L0, l0, s); end % 沿特定方向求一维结构函数 r (1:n/2-1)*dx; for j 1:length(r) diff_phz phz_all(1:n/2-j, n/2s, :) - phz_all(1j:n/2, n/2s, :); D_emp(j) mean(diff_phz(:).^2); end loglog(r, D_emp); hold on; D_theory 6.88 * (r / r0).^(5/3); % r0 为大气相干长度如果实测结构函数在短距离开端偏离理论值通常意味着内尺度设置得太大或dx不够细。如果长距离处饱和则是外尺度L0取值小于网格物理尺寸所致。5.3 随机种子与可复现实验相位屏用randn生成如果不固定种子每次运行结果都不同这在做参数扫描时是灾难。我会把随机种子写进函数参数并单独维护一个rng_setting.mat文件记录每个实验的种子。这样调一个参数时其他随机条件不变唯一变量就是被调参数结果差异才能归因。对于强湍流情况相位屏乘在光场上可能产生强度突然到零的区域此时分步傅里叶法里真空传播步长要减小否则衍射会把相位不连续点放大成非物理的强度尖峰。经验法则是每步的Rytov方差变化量不超过 0.05即dz_max (0.05 / (1.23*Cn2*k^(7/6)))^(6/11);算出来如果dz_max比预设的dz小就应缩小步长。这套参数调试逻辑放之四海而皆准不管用的是高斯还是环形光束最终核查基于同一个原则让数值统计量收敛到解析预报值的单调区间内。本文还有配套的精品资源点击获取