企业官网建设流程全解析

简介：aotools是面向光学与天文科研场景的Python自适应光学工具库，专用于大气湍流相位扰动模拟、哈特曼波前传感器响应仿真及变形镜闭环控制算法测试。资源共70个文件，以45个py源码模块为核心，涵盖湍流生成、Zernike拟合、光学传播、斜率协方差等实现；另含12个rst文档及配置文件，便于快速了解API与工程组织，压缩包仅103KB。已有520人学习下载。借助该工具包，研究者可在不同湍流强度和望远镜口径条件下生成仿真波前，评估哈特曼传感器采样与重构效果，并设计变形镜控制策略以补偿像差；同时可结合天文观测、图像处理等扩展模块，完成从湍流模拟到成像质量优化的完整链路验证。代码结构清晰、依赖精简，目录按功能划分，适合自适应光学方向的科研预研、教学演示与算法二次开发。

1. 大气湍流模拟为什么是自适应光学仿真的第一道坎

做自由空间光通信或天文 AO 终端的人，大概率经历过这个尴尬：样机装好了，控制代码在跑，唯独没有一段能代表真实大气的波前。实验室里用分划板和透镜组造不出大尺度翘曲加高频毛刺的空间结构，aotools 就是补这个缺口的 Python 工具包，pip install aotools就能拿到，vscode 里把解释器指到虚拟环境后import aotools就通了。这套 adaptive tools 覆盖湍流相位屏、哈特曼斜率到 Zernike 重构，适合做自由空间光通信、天文观测和激光传输仿真的人。

一个反直觉的事实是：AO 仿真成败大多由湍流屏的统计正确性决定，而不是由控制算法决定。相位屏的频谱和时间演化不对，后面算出来的残差再漂亮也没有参考价值。下面把链路拆成四段，从相位屏一路走到闭环残差，每段都给参数和能直接跑的代码。

2. 用 aotools 生成大气湍流相位屏：Kolmogorov 谱、外尺度与动态演化

2.1 相位屏的物理基础：为什么白噪声叠加不出湍流

大气折射率起伏服从 Kolmogorov 统计，对应到波前相位，功率谱密度近似为 Φ(k) ∝ k^(-11/3)。这意味着能量高度集中在低频端：大气湍流的整体形态是米级的大尺度弯折，毫米级的小尺度起伏只是叠加在上面的毛刺。如果直接生成随机白噪声再平滑，得到的屏会太平滑、太均匀，低阶像差占比和真实大气差一个量级。

常见的相位屏生成做法是在频域滤波：生成复高斯随机场，乘上振幅正比于 √Φ(k) 的传递函数，再做一次逆傅里叶变换回到空间域。自己写这个流程有三个坑：FFT 归一化因子容易错、低频采样不足导致大尺度能量缺失、外尺度截断方式不对会引入虚假的周期边界。aotools 把相位屏封装成类，内部补齐了低频补偿（subharmonic 方法），这几个坑基本不用自己处理。

这里唯一必须理解的物理参数是弗里德参数 r0，它同时决定湍流的强度尺度。r0 越小湍流越强，相位屏起伏越大；0.15 m 对应中等强度大气，地面光通信链路经常在 0.05 m 到 0.2 m 之间取值。

2.2 生成静态相位屏：参数配置与最小代码

静态相位屏用于单帧像差分析，比如验证哈特曼的斜率计算和 Zernike 拟合。最小代码是这样的：

import numpy as np from aotools.turbulence.infinitephasescreen import PhaseScreenKolmogorov nx = 256 pixel_scale = 0.01 # 每像素对应的实际尺寸，单位米 r0 = 0.15 # 弗里德参数，单位米 L0 = 25.0 # 外尺度，单位米 ps = PhaseScreenKolmogorov( nx=nx, pixel_scale=pixel_scale, r0=r0, L0=L0, ) phase = ps.scrn print(f"相位 RMS: {phase.std():.3f} rad")

代码逻辑并不复杂：构造PhaseScreenKolmogorov时内部完成频谱初始化，ps.scrn取出当前相位屏，返回值是一个nx × nx的二维数组，单位是弧度。打印出来的 RMS 在几弧度到十几弧度的量级是正常的；如果发现 RMS 接近 0.x 甚至 1e-3，大概率这个版本返回的是 waves 而不是 radians，需要自行除以波长再乘 2π 换算。

参数配置是这里最容易失控的地方，几个关键参数整理出来：

提示：纯 Kolmogorov 谱在零频附近能量趋向无穷，实际大气存在有限外尺度，所以构造时传入L0。如果关心低频段的准确性，可以换用PhaseScreenVonKarman，参数接口基本一致。aotools 不同版本里导入路径调整过，在infinitephasescreen里找不到类时，直接看 aotools/turbulence 目录下 phase 相关模块，类名和参数差异不大。

2.3 动态相位屏：Taylor 冻结流假设与时间步进

静态屏只能做单帧分析，闭环仿真必须处理时间演化。工程上最常用的动态模型是 Taylor 冻结流假设：把湍流看作一整段以恒定风速平移的冻结结构，某个观测点看到的时间变化来自湍流屏的空间平移。aotools 的动态相位屏在傅里叶域做相位增量，近似满足这个假设，比每一帧重新生成整个频谱省得多。

wind_speed = 10.0 # 风速，单位米/秒 wind_direction = 0.0 # 风向，单位弧度，0 表示沿 x 方向 ps = PhaseScreenKolmogorov( nx=nx, pixel_scale=pixel_scale, r0=r0, L0=L0, wind_speed=wind_speed, wind_direction=wind_direction, ) for k in range(100): ps.add_phasescreen() phase = ps.scrn # 取当前帧相位

每次调用add_phasescreen，相位屏向前推进一个时间步。推进的实际空间距离由风速和步长时间换算得到：每帧平移像素数 = wind_speed × dt / pixel_scale。这一步换算会在第 5 章回到哈特曼采样问题上，现在只需要记住这个公式。

3. 哈特曼波前传感器建模：从相位屏到斜率向量的最小实现

3.1 哈特曼为什么只输出斜率向量

Shack-Hartmann 波前传感器由微透镜阵列和探测器组成，每个微透镜把对应子孔径的光聚焦成一个光斑。当子孔径内波前是理想平面波时，光斑落在焦点中心；波前出现局部倾斜时，光斑横向偏移，偏移量正比于子孔径内波前梯度的均值。所以哈特曼输出的是 x/y 两个方向的斜率向量，并不直接给出相位。

这个设计的实际意义在速度。探测器上所有光斑的质心可以并行计算，斜率向量维度被压缩到子孔径数 × 2，控制回路才能在千赫兹量级跑完。在仿真里模拟哈特曼，最直接的做法不是生成光斑图像再做质心，而是对相位屏按子孔径切块、求梯度平均，两步下来结果完全等价。

3.2 最小实现：切块、梯度、区域平均

下面这段代码不依赖 aotools 的具体版本，因为哈特曼测量的本质就是这三步：

def measure_slopes(phase, n_subap=8): """把相位屏切成 n_subap x n_subap 个子孔径，返回 x/y 方向斜率。""" ny, nx = phase.shape sh, sw = ny // n_subap, nx // n_subap sub = phase[:n_subap * sh, :n_subap * sw].reshape( n_subap, sh, n_subap, sw ).transpose(0, 2, 1, 3) gy, gx = np.gradient(sub, axis=(2, 3)) sx = gx.mean(axis=(2, 3)) sy = gy.mean(axis=(2, 3)) return sx, sy

逻辑分三步。先裁掉不能整除的边角，再把相位矩阵 reshape 成四维，前两维是子孔径行列，后两维是子孔径内部像素；然后对每个子块做中心差分梯度，相当于真实系统中光斑质心偏移；最后把子孔径内的梯度平均成一个斜率值，对应探测器上一个光斑的净偏移。真实系统的质心算法在这里被梯度平均替代，物理含义一致。

子孔径数的选择要跟变形镜驱动器数和采样分辨率匹配：

注意sx、sy的形状是[n_subap, n_subap]，闭环前需要拼接成一维向量。拼接顺序一旦定下来，后面响应矩阵和控制矩阵全部按这个顺序排，中间不能改。

提示：忘掉裁边那一行，相位屏尺寸不能被整除时reshape会直接抛错。换用实际探测器光斑数据时也一样，边角上光斑不全的子孔径要么剔除，要么单独标记。

3.3 zonal 与 modal：斜率怎么回到波前

得到斜率向量之后有两条路线。zonal 重构把斜率当作相位梯度，通过求解泊松方程积分出完整相位分布，适合波前诊断和成像补偿；modal 重构把相位展开成一组正交模式，最常用的是 Zernike 多项式，只要拟合几十个系数就能描述波前。闭环控制更常用 modal，因为变形镜驱动器数量远小于哈特曼子孔径数，控制自由度必须截断到能驱动的模式数量。

Zernike 拟合的最小实现可以自己写，前几项基函数就够理解整个流程：

def zernike_basis(mode, N): """生成前 6 项 Zernike 模式，N 为网格大小，归一化坐标在 [-1, 1]。""" yy, xx = np.mgrid[-1:1:1j*N, -1:1:1j*N] r = np.hypot(xx, yy) if mode == 1: # piston return np.ones((N, N)) if mode == 2: # x tilt return 2 * xx if mode == 3: # y tilt return 2 * yy if mode == 4: # defocus return np.sqrt(3) * (2 * r**2 - 1) if mode == 5: # astigmatism 45° return 2 * np.sqrt(6) * xx * yy if mode == 6: # astigmatism 0° return np.sqrt(6) * (xx**2 - yy**2) def fit_zernike(phase, N, n_modes=6): A = np.stack([zernike_basis(m, N).ravel() for m in range(1, n_modes + 1)], axis=1) coeff, *_ = np.linalg.lstsq(A, phase.ravel(), rcond=None) recon = (A @ coeff).reshape(N, N) return coeff, recon

设计矩阵 A 的每一列是一个模式展平后的像素向量，最小二乘求解得到每个模式的系数。这里在相位域做拟合是为了展示原理，闭环里只能拿到斜率，通常是把 Zernike 的偏导数建到斜率空间再求逆，核心的lstsq没有变化。aotools 的湍流模块下也带 Zernike 工具，能生成更高阶和 Noll 序排列的基函数，但不同版本接口差异较大，手写这套当兜底不会有兼容问题。

4. 变形镜建模与控制矩阵：响应矩阵标定与最小二乘求解

4.1 变形镜仿真：高斯影响函数叠加

变形镜由几十到上千个驱动器构成，每个驱动器单独加电压时产生的镜面局部形变叫影响函数，工程仿真里近似为高斯型。总面形是所有驱动器影响函数按电压的线性叠加，这个近似的精度足够支撑控制算法验证，没有必要上有限元模型。

影响函数有两个关键参数：驱动器间距 d 和耦合宽度 sigma。sigma/d 在 0.3 到 0.5 之间时，相邻驱动器有一定重叠，镜面连续性好；比值太小会看到明显的电极格子痕迹，比值太大会让驱动器间串扰过大，控制矩阵条件数恶化。实际建模时，sigma 通常取 0.4 倍驱动器间距。

aotools 没有内置变形镜类，常见做法是自己用 numpy 搭一个。几十行代码就能得到一个可以插值到任意采样网格的面形模型。

4.2 用 numpy 搭一个 8×8 变形镜

N = 64 # 面形网格大小，和湍流屏一致 n_act = 8 # 驱动器按 8x8 排列 d = N / n_act # 驱动器间距，单位像素 sigma = 0.4 * d # 影响函数耦合宽度 yy, xx = np.mgrid[0:N, 0:N] act_positions = [ ((i + 0.5) * d, (j + 0.5) * d) for i in range(n_act) for j in range(n_act) ] def build_dm_surface(voltage): """电压向量 -> 变形镜面形，单位与湍流屏保持一致。""" surf = np.zeros((N, N)) for v, (cx, cy) in zip(voltage, act_positions): surf += v * np.exp( -((xx - cx)**2 + (yy - cy)**2) / (2 * sigma**2) ) return surf

act_positions是 64 个驱动器的中心坐标，build_dm_surface把长度 64 的电压向量映射成N × N面形。面形网格和湍流相位屏必须完全一致，否则后面做phase_turb - dm_surface时会出现像素级错位，闭环残差怎么调都压不下去。

4.3 响应矩阵：push-pull 标定

响应矩阵 H 描述的是给第 k 路驱动器加单位电压后，哈特曼上看到什么斜率变化。真实系统里的标定方法是 push-pull：对第 k 路加 +V 测一组斜率，再加 -V 测一组，相减除以 2V。差分的意义在于消掉影响函数中的常数项和传感器偏置，仿真里做同样处理，可以避免面形里掺入静态像差。

n_subap = 8 n_slopes = 2 * n_subap * n_subap H = np.zeros((n_slopes, n_act * n_act)) for k in range(n_act * n_act): v = np.zeros(n_act * n_act) v[k] = 1.0 s_plus = np.concatenate(measure_slopes(build_dm_surface(v))) s_minus = np.concatenate(measure_slopes(build_dm_surface(-v))) H[:, k] = (s_plus - s_minus) / 2.0

H 的形状是 128×64：64 个驱动器各标定一路，每路产生 128 个斜率输出（8×8 子孔径 × x/y 两个方向）。维度对上了，后面的控制矩阵才能算。这个标定过程在仿真里只是循环 64 次，但在真实系统里是每次实验前必须执行的标定流程，代码结构完全一致。

4.4 控制矩阵：伪逆、截断与正则化

控制矩阵理论上就是响应矩阵的伪逆：

control = np.linalg.pinv(H)

直接求伪逆会在边缘驱动器上出问题。四角驱动器的有些影响函数落在有效子孔径外，H 中对应列的能量很小，伪逆会给这些列很大的增益，测量噪声被放大成很高的电压。常用的处理是截断奇异值，或者加 Tikhonov 正则化：

# 截断伪逆：小于最大奇异值 rcond 倍的奇异值直接丢弃 control = np.linalg.pinv(H, rcond=1e-3) # Tikhonov 正则化版本，lambda 一般从 0.1 到 10 之间试探 lam = 1.0 HtH = H.T @ H + lam * np.eye(n_act * n_act) control = np.linalg.solve(HtH, H.T)

两种方案的参数对闭环行为影响很大，具体表现如下：

rcond 和 λ 建议在闭环之前扫一遍，以闭环残差方差最小为准则。常见做法是从 rcond=1e-4 开始按 10 倍步长递增，画出残差随 rcond 的变化曲线找谷底。

提示：控制矩阵只在初始化阶段算一次。闭环循环里一帧斜率只是一次 64×128 矩阵乘法，再把求逆写进循环既拖慢帧率，也没有任何数值上的必要。

5. 闭环仿真里的三个必调参数：增益、采样与斜率标定

5.1 先跑一个最小闭环循环

把前面几段拼起来，一个 300 帧的闭环循环长这样。需要注意 DM 面形用的是上一帧电压，这对应真实控制系统中的一拍延迟。

gain = 0.4 dm_voltage = np.zeros(n_act * n_act) dm_surface = np.zeros((N, N)) rv_open, rv_closed = [], [] for _ in range(300): ps.add_phasescreen() phase_turb = ps.scrn # 残差波前：DM 面形已经进入光路，再被哈特曼测量 residual = phase_turb - dm_surface slopes = np.concatenate(measure_slopes(residual)) # 积分控制器：电压增量 = -gain * control @ slopes dm_voltage += -gain * control @ slopes dm_surface = build_dm_surface(dm_voltage) rv_open.append((phase_turb**2).mean()) rv_closed.append((residual**2).mean())

rv_open记录没有校正时的湍流相位方差，rv_closed记录闭环后的残差方差，两者对比能直接判断闭环是否把波前压下去了。

5.2 三个必调参数的判断标准

第一个是增益。gain 通常在 0.1 到 0.5 之间。残差随时间震荡，说明增益太大；单调缓慢下降说明太小。从 0.3 起步，每次加减 0.1，观察rv_closed中后段的均值即可。

第二个是采样与风速的匹配。每帧平移像素数 = wind_speed × dt / pixel_scale，这个值要小于子孔径尺寸的十分之一。8×8 子孔径、单个子孔径 32 像素时，阈值约 3.2 像素/帧。超过这个值，哈特曼测到的空间混叠会直接把闭环增益吃掉，残差会出现类似噪声抬升的底噪。

第三个是单位一致性。相位屏是弧度，DM 面形也是弧度，两边才能在残差里直接相减。如果 aotools 版本返回的是 waves，先统一换算再做闭环。斜率向量的排列顺序从响应矩阵到控制矩阵一路绑死，中间改过一次顺序就要重新标定 H。

5.3 用开环/闭环对比验证闭环是否真的在工作

要验证变形镜和控制矩阵没有搭错，最快的方式是在同一段动态相位屏上对比开环方差和闭环残差方差。闭环收敛后残差方差应该是开环的 1/10 到 1/100，对应 Strehl 比：

strehl = np.exp(-np.var(phase_turb - dm_surface))

残差方差 0.1 rad² 时 Strehl 约 0.9，0.5 rad² 时掉到 0.6 附近。如果闭环残差比开环还大，不要急着调增益，先回头检查响应矩阵维度、单位换算和斜率排列顺序。把这套开环/闭环对比脚本固定成回归任务，之后每次改子孔径数、驱动器数、r0 或风速，重跑一遍并对比曲线，曲线偏离基线就说明这次改动引入了问题。

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