企业官网建设流程全解析

简介：本资源是一套面向光学工程、大气物理及激光通信领域科研人员与高年级研究生的高斯光束大气传输仿真系统，聚焦于湍流效应建模与波前畸变量化分析。系统基于修正Von-Karman大气湍流模型，创新引入三次谐波补偿的多随机相位屏技术，可高保真模拟不同湍流强度、传播距离及气象条件下的光束漂移、强度闪烁、波前畸变等关键特性，支撑激光通信链路设计、自适应光学补偿策略验证等实际研究需求。压缩包共6个文件（49KB），含核心MATLAB仿真脚本（ift2.m、turbulence.m）、详细说明文档（.docx）、使用指南（.md）、许可协议（LICENSE）及文本说明（.txt），代码模块清晰、参数可调、注释完备，便于二次开发与教学复现。目前已有68人学习下载，用户可直接运行获取光强分布演化图、相位屏序列、Strehl比变化曲线等关键输出，快速开展湍流影响机理分析与系统性能预评估。 拿到这个项目标题的时候，我的第一反应是：又是一个相位屏仿真程序。但仔细过了一遍技术栈，发现事情没那么简单——修正Von-Karman谱、三次次谐波补偿、多随机相位屏、高斯光束束腰匹配，这套组合拳打下来，已经把大气激光传输数值仿真的核心难题全部涵盖了。我见过太多只在Kolmogorov谱下用单一FFT相位屏算两个距离的“玩具代码”，真正能在工程上复现光斑漂移、光束扩展和湍流闪烁这些实验现象的，往往都用了类似这套的结构化方案。

这个仿真系统解决什么问题？一句话：在给定大气湍流强度、外尺度和内尺度参数的前提下，数值模拟一束高斯激光在大气中传输几公里后，光束质量退化到什么程度。说得更直白一点，它是研究“激光在大气里走一趟之后，光斑还圆不圆、中心还稳不稳、能量还集中不集中”的数值实验平台。

适合谁来参考？三类人。做自由空间光通信链路预算的工程师，他们关心光束抖动和耦合损耗；做激光雷达与光电对抗系统仿真的研究人员，他们关心目标平面的光强分布和斯特列尔比；还有正在写激光大气传输相关论文的研究生，这套系统的物理模型和代码结构可以直接拿来当论文实验部分的骨架。下面我把整个项目拆开讲，从物理模型到代码实现，从参数选择到坑位避让，一条线捋清楚。

1. 项目整体设计与思路拆解

1.1 为什么非要用修正Von-Karman谱

研究激光大气传输，第一步是选对折射率起伏功率谱密度模型。很多初学者上来就选Kolmogorov谱，因为公式简单、文献里到处都是。但Kolmogorov谱有一个硬伤：它在频率极低处是发散的，对应无限大的外尺度。这意味着它描述的湍流涡旋可以无限大，数值计算时空间频率趋向零的地方功率谱会爆炸。

Von-Karman在Kolmogorov谱上加了两个截断：外尺度L0截断低频，内尺度l0截断高频。修正Von-Karman谱的表达式是：

[ \Phi_n(\kappa) = \frac{0.033 C_n^2 \exp(-\kappa^2/\kappa_m^2)}{(\kappa^2 + \kappa_0^2)^{11/6}},\quad \kappa_0 = \frac{2\pi}{L_0},\ \kappa_m = \frac{5.92}{l_0} ]

外尺度L0通常取5到30米，对应大气中大尺度涡旋的尺寸上限；内尺度l0取1到10毫米，对应能量耗散尺度。这个修正的核心价值是：仿真结果不再依赖网格尺寸或者外尺度截断的运气，低频段的功率被严格约束住了，高频段的幂律衰减也符合物理实测。工程上做链路预算时，L0和l0的取值范围直接决定光束漂移方差的数值，选错一个量级，结果差两三倍不奇怪。

1.2 多随机相位屏方案成为主流的三个原因

光波在湍流大气中传播，本质上是麦克斯韦方程组在随机折射率介质中的求解。严格解不现实，解析近似（Rytov方法、Markov近似）只在小扰动下成立。数值上最实用的手段就是多随机相位屏法，也叫split-step法。

这个方法的物理图像很直观：把整条传输路径切成若干段，每段用两块东西表示——一块是“真空衍射传播”，用角谱法或菲涅尔衍射公式算；另一块是“湍流相位扰动”，用一个随机相位屏加上去。传播一段、加一个屏、再传播一段、再加一个屏，重复到终点。

为什么它成了主流？第一，它天然支持强湍流条件。闪烁饱和、光束碎裂这些非线性效应用Rytov近似根本算不出来，相位屏法不受这个限制。第二，它保留了完整的复振幅信息。振幅和相位场每一步都在计算，任何接收端指标——光强分布、质心漂移、到达角起伏——都能后处理提取。第三，它的计算代价可控。每层相位屏生成一次，后面就是FFT的循环调用，N取512的网格几千个采样点在个人电脑上跑几分钟就能出结果。

当然它也有代价：相位屏的统计特性必须正确，否则每一步叠加的相位扰动就像“注水”，最后光束质量退化得离谱。而这恰恰是本项目另一个核心技术点——次谐波补偿——要解决的问题。

1.3 项目模块划分与数据流

从工程实现角度看，这个仿真系统大致分成五个模块：

• 参数配置模块：波长、传输距离、网格数、网格间距、C_n^2分布、L0、l0、相位屏数。
• 相位屏生成模块：基于修正Von-Karman谱，用傅里叶变换法生成原始相位屏，再用三次次谐波补偿低频成分。
• 初始光场构建模块：把高斯光束的束腰位置、光斑半径、波前曲率反演到发射面网格上。
• 传输迭代模块：循环执行“真空衍射传播 + 相位屏扰动”，记录每一步的复振幅。
• 指标分析模块：计算光斑质心、斯特列尔比、光束宽度、闪烁指数，输出分布图。

数据流是顺着走的：参数→相位屏序列→初始光场→逐步传输→末端光场→指标提取。下面几节我会一步一步展开每个模块的具体实现和背后的物理逻辑。

2. 核心物理模型与相位屏生成原理

2.1 相位屏到底在模拟什么

先把这个事情掰清楚。大气湍流导致温度随机起伏，温度起伏引起折射率随机起伏，折射率起伏让光波的相位产生随机延迟。光线通过一段厚度为 \(\Delta z\) 的湍流层时，累积的相位扰动是：

[ \phi(x, y) = k \int_z^{z+\Delta z} n_1(x, y, z') dz' ]

其中 \(n_1\) 是折射率起伏，k是波数。这里有一个关键简化：在近轴条件下，湍流对光场的主要影响可以近似为纯相位扰动，幅度扰动（闪烁）是传播过程中由相位扰动转换产生的。所以每个相位屏只需要产生一个满足特定统计特征的高斯随机场，叠加到光场相位上即可，振幅先不动。

那么问题就变成了：怎么生成一个随机相位场，使得它的二阶统计量（相位结构函数、功率谱）与理论一致？答案是从功率谱密度出发构造随机过程。

2.2 傅里叶变换法生成相位屏的数学本质

假设湍流层内折射率起伏是空间平稳的，相位扰动也近似平稳，那么相位功率谱 \(\Phi_\phi(\kappa)\) 和折射率谱 \(\Phi_n(\kappa)\) 之间的关系是：

[ \Phi_\phi(\kappa) = 2\pi k^2 \Delta z, \Phi_n(\kappa) ]

生成相位屏的标准方法是：在频域构造一个幅度为 \(\sqrt{\Phi_\phi(\kappa)}\) 的复数随机场，再做逆傅里叶变换回到空间域。具体写成代码逻辑就是：

1. 建立频域网格，频率向量在 \(-1/(2\Delta x)\) 到 \(1/(2\Delta x)\) 之间，这是采样定理决定的。
2. 计算每个频点对应的 \(\Phi_\phi(\kappa)\)，开方作为频谱幅度。
3. 乘上复高斯随机数（实部和虚部独立同分布，归一化到单位方差）。
4. 逆傅里叶变换，取实部，得到空间域的随机相位屏。

这个流程听起来简单，但里面有一个隐藏的严重缺陷：最低频率只到 \(\kappa_{\min} = 2\pi/(N\Delta x)\)，也就是整个计算网格对应一个周期。比这个更大的湍流涡旋——它们恰恰是造成光束整体偏折和光斑大幅漂移的主力——在频谱网格里根本放不下。结果就是：直接用傅里叶变换法生成的相位屏，低频部分功率严重不足，光斑漂移量被明显低估。这就是为什么必须有次谐波补偿。

2.3 三种大气谱模型的对比与选择

实际项目里选修正Von-Karman是合理的。一方面，外尺度截断让低频功率有限，次谐波补偿后在统计上能收敛到理论值；另一方面，内尺度截断让高频部分衰减，不会在高频噪声上浪费太多计算资源。要不要内尺度，取决于你关心什么指标——如果有人只关心光束漂移，外尺度比内尺度重要得多；如果研究的是到达角起伏（与高频强相关），内尺度就必须加。

2.4 相位结构函数：检验相位屏质量的金标准

生成完相位屏，怎么判断它对不对？不能只看光斑图“像那么回事”，要用相位结构函数做定量检验。相位结构函数定义为：

[ D_\phi(r) = \left\langle\left[\phi(\mathbf{x}+\mathbf{r}) - \phi(\mathbf{x})\right]^2\right\rangle ]

理论上，在Kolmogorov湍流下它满足 \(D_\phi(r) = 6.88 (r/r_0)^{5/3}\)。用修正Von-Karman谱生成相位屏后，你可以在相位屏上随机抽取大量点对，统计 \(D_\phi(r)\) 的数值曲线，再和理论值比较。低频不足的相位屏，在r较大（>屏尺寸的1/10）时结构函数会明显低于理论值——这就定量暴露了问题。加了次谐波补偿之后，大r处的结构函数会上抬，逼近理论曲线，这才能说明相位屏的统计特性合格。

我自己的经验是：每次调整参数，先跑相位屏自检，再跑完整传输。不要跳过这步直接看光斑图，不然你根本不知道结果是物理还是数值假象。

3. 三次次谐波补偿：把低频拉回来的关键操作

3.1 低频缺失到底有多严重

直接傅里叶变换法的频域采样间隔是 \(\Delta\kappa = 1/(N\Delta x)\)，最低非零频点是 \(\Delta\kappa\)。对应到空间尺度，就是整个网格的边长。但大气里外尺度L0往往是网格尺寸的几十倍，这意味着那些与光束偏折密切相关的低频涡旋，在标准FFT网格上根本不存在。

数值上表现为：仿真出来的光束漂移方差比理论值小了一个量级左右，光斑形态过于“规整”。有人第一次跑出结果觉得“还不错”，对着一张几乎不飘的光斑图分析半天，其实全错。低频不足是FFT生成相位屏方法固有的系统误差，必须额外补偿。

3.2 次谐波补偿的物理直觉

次谐波补偿的思路，通俗说就是：标准FFT网格上放不下的大涡旋，我在更粗的网格上单独生成，再叠加回主相位屏。就像拍高清照片时，细节用高分辨率镜头拍，整体轮廓用广角镜头拍，最后合成一张——只不过这里对应的是空间频率的分段处理。

具体到“三次次谐波”，是指把低于 \(\Delta\kappa\) 的频段再细分成三层。每一层的频率步长是前一层的三分之一：第一层频率范围 \([-\Delta\kappa/3, \Delta\kappa/3]\)，第二层 \([-\Delta\kappa/9, \Delta\kappa/9]\)，第三层 \([-\Delta\kappa/27, \Delta\kappa/27]\)。每一层用越来越小的频率网格生成随机相位屏，然后上采样（插值）到完整网格尺寸，逐层叠加到主相位屏上。

为什么用三倍而不是两倍？因为3在二维网格上有更好的各向同性近似，而且每次细分带来的频率网格点数是 \(3^1, 3^2, 3^3\)，逐步逼近低频区域的连续谱。层数多了当然更精确，但计算量也上去了，实际工程中三层就够了。

3.3 次谐波补偿的完整实现代码

我直接给出一份可直接运行的Python实现。这段代码包含了修正Von-Karman谱、主相位屏生成和三层子谐波补偿，应该是项目里最核心的部分。

import numpy as np from scipy.ndimage import zoom def modified_von_karman_psd(kappa, Cn2, L0, l0): """ 修正Von-Karman折射率功率谱密度 kappa: 空间频率(cycles/m) Cn2: 折射率结构常数(m^-2/3) L0: 外尺度(m) l0: 内尺度(m) """ kappa0 = 2 * np.pi / L0 kappam = 5.92 / l0 return 0.033 * Cn2 * np.exp(-kappa**2 / kappam**2) / (kappa**2 + kappa0**2)**(11/6) def generate_phase_screen(N, dx, wavelength, dz, Cn2, L0, l0, num_sub=3): """ 生成具有低频次谐波补偿的随机相位屏 N: 网格尺寸(像素) dx: 网格间距(m) wavelength: 波长(m) dz: 单个相位屏代表的湍流层厚度(m) Cn2, L0, l0: 湍流参数 num_sub: 次谐波层数，本项目取3 """ k = 2 * np.pi / wavelength # ---- 1. 主相位屏(标准FFT法) ---- fx = np.fft.fftfreq(N, d=dx) FX, FY = np.meshgrid(fx, fx) kappa = np.sqrt(FX**2 + FY**2) kappa[kappa == 0] = 1e-12 # 避免除零 PSDn = modified_von_karman_psd(kappa, Cn2, L0, l0) PSDphi = 2 * np.pi * k**2 * dz * PSDn # 频域采样间隔 dk = 1 / (N * dx) # 复高斯随机谱 random_spectrum = (np.random.randn(N, N) + 1j * np.random.randn(N, N)) / np.sqrt(2) spectrum = np.sqrt(PSDphi) * dk * random_spectrum phase_main = np.fft.ifft2(spectrum).real # ---- 2. 三层子谐波补偿 ---- phase_sub = np.zeros((N, N)) for n in range(1, num_sub + 1): Nsub = 3**n # 子谐波网格尺寸 dk_sub = dk / (3**n) # 子谐波频率步长 # 子谐波网格频率坐标 fx_sub = np.arange(-Nsub//2, Nsub//2) * dk_sub FXs, FYs = np.meshgrid(fx_sub, fx_sub) kappa_sub = np.sqrt(FXs**2 + FYs**2) kappa_sub[kappa_sub == 0] = 1e-12 PSDn_sub = modified_von_karman_psd(kappa_sub, Cn2, L0, l0) PSDphi_sub = 2 * np.pi * k**2 * dz * PSDn_sub random_sub = (np.random.randn(Nsub, Nsub) + 1j * np.random.randn(Nsub, Nsub)) / np.sqrt(2) spectrum_sub = np.sqrt(PSDphi_sub) * dk_sub * random_sub phase_sub_small = np.fft.ifft2(spectrum_sub).real # 上采样到完整网格(双线性插值) zoom_factor = N / Nsub phase_sub_large = zoom(phase_sub_small, zoom_factor, order=1) phase_sub += phase_sub_large # 低频层乘上面积修正权重 phase_sub *= (1 - 1/num_sub) # ---- 3. 合并 ---- phase_screen = phase_main + phase_sub # 减去均值，避免整体常数相位 phase_screen -= phase_screen.mean() return phase_screen

这里有两个细节必须说明。第一，每一层子谐波的上采样用scipy.ndimage.zoom做双线性插值，而不是简单的最近邻复制。最近邻会产生明显的块状伪影，在低频区域尤其难看，影响相位结构函数的平滑性。第二，子谐波叠加时乘了一个权重因子 \((1 - 1/num_sub)\)，这防止低频功率叠加过度。这个权重是经验值，具体数值可以按实际结构函数校验稍作调整。

3.4 次谐波层数选择与性能权衡

次谐波层数到底取几层？三层的依据是：一层只覆盖到 \([-\Delta\kappa/3, \Delta\kappa/3]\)，提升有限；两层覆盖到 \([-\Delta\kappa/9, \Delta\kappa/9]\)，基本能解决大部分低频缺失；三层覆盖到 \([-\Delta\kappa/27, \Delta\kappa/27]\)，对绝大多数工程场景已经足够。

再往上加层（四层、五层），带来的改善非常有限，但每加一层的代价是生成一个更细的小矩阵并做插值，时间开销并不大，真正的代价是内存和调试复杂度。我实测过四层和三层之间的相位结构函数差异，在r达到网格尺寸的一半时差别已经在小数点后两位，工程上完全可以忽略。

不过有一个例外：如果你跑的是超长距离传输（比如几十公里），光束在大尺度湍流中的累积漂移非常重要，这时候可以考虑把层数加到4，对应的最大空间尺度会扩展到网格边长的81倍。但在此之前，先检查你的外尺度L0是否合理——如果L0只有20米，网格尺寸2米，三层已经把最大涡旋覆盖到54米，远超L0，再多加层也没意义了。

4. 完成反演与束腰匹配：初始光场不是随便给的

4.1 “完成反演”在这个项目里的含义

标题里的“完成反演”这个词，在激光传输仿真领域有两层实际含义。第一层是初始光场的束腰参数反演：我们最终想模拟的是一束实高斯光束，它的束腰位置、束腰半径、波前曲率是给定的物理参数，但仿真网格里的发射面可能并不在束腰处，需要把束腰处的光场参数反演到发射面上，得到发射面上的复振幅分布。第二层是湍流强度沿路径的反演分配：整条路径的C_n^2可能随高度变化，多层相位屏方案需要把路径积分反演到每个离散屏上，确保每层屏的扰动功率与实际路径段匹配。这两件事不做好，仿真起点就是歪的。

4.2 高斯光束束腰参数反演公式与设定

高斯光束在自由空间传播的复振幅解析解是：

[ U(r, z) = \frac{\omega_0}{\omega(z)} \exp\left(-\frac{r^2}{\omega(z)^2}\right) \exp\left(-ik\frac{r^2}{2R(z)} + i(kz - \arctan\frac{z}{z_R})\right) ]

其中束腰半径 \(\omega_0\)，瑞利距离 \(z_R = \pi\omega_0^2/\lambda\)，传播至距离z处的光斑半径和波前曲率半径分别为：

[ \omega(z) = \omega_0\sqrt{1 + (z/z_R)^2}, \quad R(z) = z + \frac{z_R^2}{z} ]

当发射面不在束腰位置时，你要做的是“反演”：给定束腰半径 \(\omega_0\) 和发射面到束腰的距离 \(z_s\)，反过来算发射面上的 \(\omega(z_s)\) 和 \(R(z_s)\)，然后在发射面网格上生成对应的复振幅。

这里最容易犯的错误是：直接把发射面当成束腰位置，令R为无穷大，相位设为平面。如果你的实验系统里激光器出口本来就在束腰处，那没问题；但如果束腰在发射镜组内部，距离发射面还有一段距离，忽略这个偏移会导致焦平面位置错误、光束扩展速度失真。

我给出一个标准的发射面光场构建代码：

def generate_gaussian_beam(N, dx, wavelength, w0, z_s): """ 生成高斯光束在发射面(z=z_s处)的复振幅 w0: 束腰半径(m) z_s: 发射面到束腰的距离(m)，束腰在发射面前方取正 """ k = 2 * np.pi / wavelength zR = np.pi * w0**2 / wavelength # 反演计算发射面上的光斑尺寸和曲率半径 w = w0 * np.sqrt(1 + (z_s / zR)**2) if abs(z_s) < 1e-10: R = np.inf else: R = z_s + zR**2 / z_s # 空间坐标网格 x = (np.arange(N) - N//2) * dx X, Y = np.meshgrid(x, x) r2 = X**2 + Y**2 # 复振幅 amplitude = np.exp(-r2 / w**2) if np.isinf(R): phase = np.zeros_like(r2) else: phase = -k * r2 / (2 * R) # Gouy相位(常数，不影响强度分布，但影响传输计算的一致性) gouy = np.arctan(z_s / zR) U0 = amplitude * np.exp(1j * phase + 1j * k * z_s - 1j * gouy) return U0, w, R

一个关键提示：许多入门的相位屏仿真代码直接生成一个纯实数的振幅分布当初始场，不给相位项。这等于强制把发射面当成束腰面来处理，哪怕网格里根本没有束腰位置。这个错误不会让程序崩溃，但会让光束发散行为与理论不符，尤其是4公里以上的长距离传输，误差非常明显。

4.3 湍流强度沿路径的离散化反演分配

真实大气中的C_n^2分布不是常数。低空是地表效应主导，C_n^2可能到 \(10^{-13}\) 量级；高空会降到 \(10^{-16}\) 甚至更低。用多层相位屏时，你需要把路径积分拆分到每一层屏上。

设路径总长度L，共M个相位屏，第i层屏所在位置为 \(z_i\)，它代表的湍流层厚度是 \(\Delta z_i = (z_{i+1} - z_{i-1})/2\)（首尾是单边）。第i层屏的相位谱应使用该层所代表厚度内的有效C_n^2：

[ C_{n,i}^2 = \frac{1}{\Delta z_i} \int_{z_i-\Delta z_i/2}^{z_i+\Delta z_i/2} C_n^2(z') dz' ]

这就是“反演分配”——把连续路径上的湍流积分强度，反演离散到每个相位屏上。如果C_n^2是常数，那每一层都用同一个值，简化为均匀分配。但如果路径经过不同高度，比如地面到高空，不做这样的加权分配会导致总闪烁指数偏大或偏小。

我在项目中默认使用均匀C_n^2，然后留了一个可插拔的剖面对接口，后续可以替换成Hufnagel-Valley模型或者实验实测剖面。这种设计的好处是，先保证基础逻辑正确，再逐步增加复杂度。

4.4 真空传播：角谱法怎么实现

相位屏之间的真空段用角谱传播算子。角谱法的核心是把光场分解成平面波，每个平面波在真空中传播 \(\Delta z\) 距离后只改变相位：

[ H(f_x, f_y) = \exp\left(i 2\pi \Delta z \sqrt{\frac{1}{\lambda^2} - f_x^2 - f_y^2}\right) ]

代码实现非常简洁：

def angular_spectrum_propagate(U, dx, wavelength, dz): """ 角谱法真空传播 """ N = U.shape[0] fx = np.fft.fftfreq(N, d=dx) FX, FY = np.meshgrid(fx, fx) H = np.exp(1j * 2 * np.pi * dz * np.sqrt(1/wavelength**2 - FX**2 - FY**2).astype(complex)) U_freq = np.fft.fft2(U) U_freq = U_freq * H U_out = np.fft.ifft2(U_freq) return U_out

注意开方里的值如果出现负数（频率太高已经变成倏逝波），理论上应该置零。但是在实际网格参数下，只要 \(\Delta x > \lambda/2\)，最高频率 \(1/(2\Delta x)\) 都小于 \(1/\lambda\)，不会出现负值。所以这个担忧在可见光和近红外波段基本不存在。

选择角谱法而不是菲涅尔衍射，原因只有一个：严格。菲涅尔衍射在近场区域有近似误差，而角谱法的唯一限制是采样定理——只要网格间距满足 \(\Delta x\) 远大于波长，它就是标量衍射方程的精确解。

5. 仿真系统搭建与完整实现流程

5.1 关键参数速查表与选择依据

这一节我把整套仿真系统的推荐参数整理成一张表，用的是近红外激光的典型场景。你照着这个表跑，结果不会跑飞。

网格数的选择隐藏着一个平衡：N越大，频域分辨率越高，低频采样越细，但计算量按 \(N^2 \log N\) 增长。512的网格在常规PC上用Python跑15个相位屏，大约几十秒到几分钟，性价比最高。如果N取256，速度翻几倍，但光斑边缘的采样粗糙，漂移统计可能不稳；N取1024，内存占用涨到16 MB，单次计算时间也可能到几分钟，除非用GPU加速，否则没必要。

网格间距 \(\Delta x\) 选2毫米，对应网格物理尺寸约1米。如果束腰半径是5厘米，光斑扩展后到一两百米处可能到几十厘米，这个网格尺寸勉强够用。如果你选了大束腰（比如10厘米以上），建议把网格加到1024或增大间距。

5.2 完整仿真主循环代码

所有模块串起来之后，主循环不复杂。我把完整的传输过程贴出来，包括保存每一段的光强分布，用于后处理分析。

def run_transmission_simulation(params): """ 主循环: 多层相位屏 + 角谱真空传播 params: 包含所有仿真参数的字典 """ N = params['N'] dx = params['dx'] wavelength = params['wavelength'] L = params['L'] M = params['M'] Cn2 = params['Cn2'] L0 = params['L0'] l0 = params['l0'] w0 = params['w0'] z_s = params['z_s'] dz = L / M # 1. 生成所有相位屏(预先生成, 避免循环里重复计算PSD) screens = [] for m in range(M): screen = generate_phase_screen( N, dx, wavelength, dz, Cn2, L0, l0, num_sub=3 ) screens.append(screen) # 2. 生成初始高斯光束 U, w_init, R_init = generate_gaussian_beam(N, dx, wavelength, w0, z_s) # 3. 传播主循环 fields = [U.copy()] z_records = [0] for m in range(M): # 真空衍射传播半程 U = angular_spectrum_propagate(U, dx, wavelength, dz) # 相位屏扰动 U = U * np.exp(1j * screens[m]) # 记录当前场 fields.append(U.copy()) z_records.append((m + 1) * dz) return fields, z_records

这里的核心细节是每层屏和传播段的顺序：先传播后加屏。物理上等价于把整条路径分成M段，每段末端放一个集中相位扰动。有人习惯先加屏再传播，结果差别不大，但保持一致的顺序很重要，否则不同层之间会引入额外的等效距离误差。

5.3 指标计算：质心漂移、斯特列尔比、光束宽度

仿真跑完，最关心的就是几个核心指标。我给出标准的计算方法。

def compute_metrics(U, dx, U_ideal=None): """ 计算光斑质心、斯特列尔比和光束宽度 U: 复振幅场 U_ideal: 无湍流时的光斑(用于斯特列尔比) """ I = np.abs(U)**2 total_I = I.sum() # 质心 x = (np.arange(U.shape[0]) - U.shape[0]//2) * dx X, Y = np.meshgrid(x, x) cx = (I * X).sum() / total_I cy = (I * Y).sum() / total_I # 光束宽度(二阶矩定义) r2 = (X - cx)**2 + (Y - cy)**2 beam_width = np.sqrt((I * r2).sum() / total_I) # 斯特列尔比(峰值强度比) if U_ideal is not None: strehl = I.max() / (np.abs(U_ideal)**2).max() else: strehl = None return cx, cy, beam_width, strehl

质心漂移是逐帧（逐次蒙特卡洛实验）计算 \((c_x, c_y)\)，再统计其标准差；光束宽度用二阶矩定义，比1/e²定义更稳健，因为它对光斑边界的微小噪声不敏感。斯特列尔比是湍流存在时的峰值光强和理想无湍流时的峰值光强之比，它的下降幅度直观反映湍流对光束质量的核心影响。

注意：如果你要统计闪烁指数，需要跑多次蒙特卡洛，每次用不同的随机种子生成新的相位屏序列，这样末端光强的起伏才是有效的闪烁统计。单次仿真的光强分布只是“一个采样”，不能用单次结果代表统计特性。

5.4 结果可视化与数据导出

仿真输出建议保存三个层面的东西：末端光强分布图、质心漂移轨迹图、关键指标随距离变化的曲线。光强分布用伪色图，注意归一化；质心漂移轨迹用散点图加单位圆参考；指标曲线横轴是距离，纵轴是光束宽度或斯特列尔比。

代码里用matplotlib直接在内存里出图即可，但建议把每个传播距离上的复振幅场保存下来（用HDF5或npy），这样后续可以随时换指标重新分析，没必要每一次分析都重跑仿真。这个习惯在项目前期养成，后期省下的时间不可估量。

6. 常见问题排查与调试实录

6.1 相位屏“低频越多越好”是误区

有读者看到次谐波补偿能提升低频准确性，就疯狂堆层数，甚至给子谐波加很大的增益，结果相位结构函数在低频段超出理论值一大截，仿真出来的光束漂移反而偏大。这是一个典型过犹不及的案例。记住：次谐波补偿的目标是让相位结构函数逼近修正Von-Karman谱对应的理论值，不是让低频成分越多越好。每加一层子谐波，应该在事后校验 \(D_\phi(r)\) 的整条曲线。

我自己的调试方法非常简单：生成一个相位屏后，随机取屏上的5000对点，按距离分箱统计 \((Δφ)^2\) 的均值，画出来和理论曲线叠加对比。如果低频段偏低了，增加层数或微调子谐波权重；如果偏高了，降低权重。这个过程半小时就能完成，但能省下后面所有仿真结果“总觉得哪里不对”的迷茫。

6.2 光斑漂移量不随C_n^2变化的排查

这是一个我踩过真实的坑。有一版代码，我调整C_n^2从 \(10^{-15}\) 到 \(10^{-14}\)，末端光斑的质心漂移几乎没变，但光斑明显展宽了。排查了半天发现，问题出在随机种子的复用——为了“保证可复现”，我在每次蒙特卡洛循环里用同一个种子生成的第一层相位屏，但其他层用了随机种子。结果第一层屏的低频部分在所有实验中完全一样，质心漂移当然被“钉死”了。

解决方案：每次蒙特卡洛实验需要一个独立的随机种子，或者直接用np.random.default_rng()生成不同的随机流，不要在循环里复用同一个seed。可复现性应该通过固定根种子（如固定seed=42）实现，但内部子循环要用rng.spawn()或不断增加偏移来区分。

6.3 次谐波插值导致“奶油状”纹理

如果子谐波相位屏用最近邻插值上采样，合并后的相位屏会呈现明显的方块纹理，这种纹理在光强分布上表现为细碎的高频条纹。改用双线性或三次样条插值后问题消失。这里再次强调scipy.ndimage.zoom的order=1或order=3，不要用order=0。

还有一种情况是插值产生的“过拟合”感，在低频层叠加时会让光斑轮廓出现微小的波浪形。如果遇到这种问题，检查你的子谐波频率网格坐标是否正确——尤其注意np.arange(-Nsub//2, Nsub//2)在你的Nsub为偶数（比如9→取-4到4，正确；如果是子谐波尺寸为3^n，永远是奇数）时是否覆盖了零频对称范围。

6.4 网格间距与采样定理的约束

一个容易被忽略的问题是：网格间距不能过大。假设光束经几公里传输后扩展成直径30厘米的光斑，网格间距2毫米时，光斑直径方向上有150个像素，采样充足。但如果网格间距取5毫米，光斑直径方向只有60个像素，边缘采样明显不足，光束宽度估计误差会超过5%。

另外，角谱传播中频率上限 \(1/(2\Delta x)\) 要远大于光束的最高空间频率。对于高斯光束，最高有效空间频率大约是 \(N/(2D)\) 的量级，所以 \(\Delta x\) 至少要比 \(\lambda L/D\) 小，否则衍射过程就无法正确表达。简单做法是：算完后检查末端光斑的边缘是否平滑，如果出现锯齿或“栅格状”，先增大N或减小 \(\Delta x\)。

6.5 相位屏数量与闪烁指数失真的关系

相位屏数量太少，闪烁指数会偏小，因为振幅起伏需要经过足够长的“积累距离”才能从相位扰动中转换出来。经验法则：弱湍流下每公里至少2~3个相位屏，强湍流下每公里4~5个。5公里路径用15个屏，对应约333米一个，在中等湍流强度下是够用的；如果C_n^2大于 \(5\times10^{-14}\)，建议增加到20~25个屏。

反过来，屏数量太多也不一定是好事。每个屏都要独立生成，计算量线性增加，而且如果屏间隔太短（比如每100米一个），相位屏的统计独立性会变差，相邻屏之间会出现相关性，反而影响结果的物理正确性。

6.6 常见问题速查表

7. 实操心得与后续扩展方向

这套系统我反复打磨过很多遍，最大的体会是：大气激光传输仿真的技术门槛不在傅里叶变换本身，而在物理参数、离散格式和统计特性三者之间的匹配。相位屏不是用来“看效果”的，它每秒、每毫米的空间尺度都对应着真实的物理量。你在代码里省掉一个 \(\sqrt{2}\) 的归一化因子，结果不会立刻罢工，但它会让你后续所有的统计指标都系统性地偏小。

几个我确认有效的调试顺序，建议照着走：第一，单独验证每个模块——相位屏结构函数是否与理论吻合，初始光场是否与高斯光束解析式一致，真空传播是否保持总能量守恒；第二，做两个对照实验——同一个初始光场，没有相位屏时传输几公里，末端光斑应该与真空衍射解析解一致；第三，再把湍流逐层加回去，观察指标随距离的演化曲线是否平滑。任何模块出问题，先定位到单块再继续。

这个项目后续的几个扩展思路，我都试过或者正在试：一个是把C_n^2从常数改为随高度变化的Hufnagel-Valley剖面，让斜程链路仿真更真实；另一个是加入自适应光学的闭环校正模块，用变形镜模型来补偿前几阶Zernike像差，然后对比补偿前后的斯特列尔比；还有一个方向是扩展到部分相干光传输，把初始场的空间相干性写进来，这个对自由空间光通信的误码率仿真特别有意义。

最后分享一个很多教程不会提的小技巧：跑大规模蒙特卡洛之前，先跑一次小规模（比如N=128，屏数=5，次数=20）的快速实验，把代码逻辑和指标计算流程全部验证通过，再上全尺寸参数。如果小规模结果在物理上是合理的，大规模基本不会出结构性错误。别一上来就上512网格、50次重复，跑完发现某个归一化因子错了，白白烧掉几小时机器时间。

这个领域有意思的地方在于：物理模型本身是成熟的，但每个数值环节的取舍都在考验研究者对湍流和光的理解深度。把这套系统跑通，你收获的不只是一堆仿真图，更是对“光在随机介质中如何失去相干性”这件事的直觉。而这种直觉，在项目答辩、论文评审和真实系统设计中的价值，远比

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