企业官网建设流程全解析

简介：这是一份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) 通常取 20～100 米；(\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 步，则每步Lz=100。n*dx是物理尺寸，必须足够大，否则相位屏边缘的统计特性会失真。randn(n)产生复高斯随机场，乘上功率谱的平方根相当于在频域整形。

2.2 分步傅里叶法：为什么不能直接乘一个相位屏

如果光在真空和湍流中交替传播，可以在一个薄屏处叠加所有扰动？不行。大气湍流是连续分布的三维不均匀体，但工程上经常用“相位屏近似”简化：将传播路径切成若干段，每段长度足够短，使得湍流效应可以集中在该段中心的薄屏上。段间用真空衍射传播，这就是分步傅里叶法。

每一步的操作是：

1. 在真空中传播距离 (L_z)，用角谱法或菲涅尔衍射积分。
2. 将光场乘以相位屏 (\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}; % 保存每步光强 end

record_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^2=10^{-14})，距离 2 km，波长 1.55 μm 的条件下，我跑过一组对比，大致结果如下：

这个结果符合一般规律：拓扑荷越高，初始相位变化越剧烈，湍流引起的相位扰动和光强闪烁越强。但并不是说涡旋光束一定更差，在无湍流时涡旋光束能保持螺旋相位结构，只是统计结果上对湍流更敏感。

4.3 把光斑序列压缩成GIF

gif.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^2=10^{-16})、(z=500) m，(\sigma_R^2) 可能在 0.01 以下，光斑基本不抖。要看到明显效应，至少让 (\sigma_R^2) 在 0.1 以上。我自己常用的“保守起效”组合是：(C_n^2=10^{-14})、距离 1000 米以上。

另一个坑是相位屏的分辨率。相位屏的频域采样间隔是 (1/(n dx))，如果这个值大于湍流内尺度对应的空间频率，小而强的相位涡旋会被平滑掉，表现为光斑只有整体漂移而没有破碎细节。以 (l_0=5) mm 为例，最大空间频率需要达到 (1/l_0=200) /m。取 (n=512)，则 (dx) 必须小于 (1/200=0.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/2+s, :) - phz_all(1+j:n/2, n/2+s, :); 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小，就应缩小步长。这套参数调试逻辑放之四海而皆准，不管用的是高斯还是环形光束，最终核查基于同一个原则：让数值统计量收敛到解析预报值的单调区间内。

本文还有配套的精品资源，点击获取