aotools在大气湍流模拟与自适应光学实验中的应用指南

简介aotools是针对大气湍流模拟与自适应光学研究而设计的Python库主要解决光在大气中传播时因湍流导致的波前畸变问题提供哈特曼探测器模拟与变形镜控制算法支持适合天文望远镜、激光通信和遥感成像等领域的科研与工程人员使用也适合相关专业高年级学生作为课题研究参考。压缩包共70个文件体积约103KB以45个Python源文件与12个RST文档为主另含单元测试、配置和构建脚本完整覆盖湍流相位屏生成、波前重建及自适应校正仿真链路。已有520人学习浏览。内容包含各模块测试与使用说明用户可通过Zernike多项式或高斯随机过程生成相位扰动模拟哈特曼传感器对波前局部曲率的测量。在此基础上还能测试变形镜闭环校正效果为不同湍流条件下的光学系统设计与参数优化提供数据基础与代码参考。1. 为什么说 aotools 是大气湍流模拟的实验拼装台同样是 256×256 的相位屏手写傅里叶变换法时最容易错的是功率谱归一化——差一个系数结构函数就从 5/3 次方变成 2 次方肉眼根本看不出来。aotools 把这类容易踩坑的函数全部收敛到aotools.functions和aotools.wfs下并且带了一整套test_*.py做验证。它不是给做产品的人用的黑盒而是给做自适应光学实验的人一个能拆开改的源码包。这个包覆盖从 Kolmogorov 湍流相位屏、哈特曼子孔径斜率测量到变形镜控制回路所需的主要环节适合刚把 Python 当作仿真主力、又不想从零重写波前重构算法的研究生和光学工程师。2. 大气湍流模拟的第一站傅里叶相位屏与协方差校验2.1 先搞清 r0、L0、l0 三个参数再调代码大气湍流对波前的影响不是白噪声而是服从 Kolmogorov 统计相位结构函数与(r/r0)^(5/3)成正比r0 是 Fried 相干长度代表波前空间尺度超过 r0 后相位差会剧烈变化。aotools 里最常用的是ft_phase_screen它用频谱法生成一个二维随机相位屏先生成高斯白噪声再乘以模拟湍流功率谱的传递函数最后做一次傅里叶逆变换。典型调用如下import numpy as np from aotools.functions import turbulence # r0: Fried相干长度(米), N: 相位屏边长(像素) # delta: 每个像素代表的实际尺寸(米/像素) phase turbulence.ft_phase_screen( r00.1, N256, delta0.002, L050.0, l00.01 ) print(phase.shape) # (256, 256), 单位是弧度这里的L0是湍流外尺度控制低频能量截止l0是内尺度截断高频。delta决定了空间采样率也直接决定模拟孔径的真实大小模拟直径是N * delta。如果delta取太大子孔径尺寸会失真取太小同一 r0 下相位起伏被过度离散导致能量泄漏。我一般按口径 1~8 米、像素数 128~512 来回调保证N * delta落在目标口径上。怎么验证生成的相位对不对不要只看图像不像湍流要算结构函数。aotools 的测试文件里给了思路统计相位屏上两点相位差的方差和理论值6.88*(d/r0)**(5/3)对比。下面是我常用的快速校验def structure_function_error(phase, delta, r0, max_sep64): errors [] for sep in range(2, max_sep): diff phase[::sep, ::sep] - phase[sep::sep, sep::sep] var np.mean(diff**2) theory 6.88 * (sep * delta / r0) ** (5/3) errors.append((sep, var, theory)) return np.array(errors)校验时注意两点一是只统计比口径小得多的间距间距接近口径后采样数骤减方差估计本身就不可靠二是6.88这个常数只在L0趋近无穷时严格成立你给L050的有限外尺度后低频部分会有偏差所以把max_sep压到口径的 1/4 以内再比较更稳。注意ft_phase_screen返回的相位单位是弧度并且默认没有去掉整体活塞项。做闭环仿真前先减去全屏平均值否则每一帧都需要重新对齐控制信号的零点。2.2 无穷相位屏和 Zernike 系数分解单帧相位屏只能做静态分析做闭环仿真需要时间上连续的无穷相位屏。aotools 的inf_phase_screen用协方差法生成符合时间相关性的相位屏序列。它的实现是通过预置随机种子和时域滤波生成下一帧让前一帧与后一帧共享低频成分从而模拟风移效应。调用时除了 r0、N、delta还需要给定风速向量和帧数。我的经验是风速向量要按像素尺寸归一化比如风速 10 m/s、delta0.002 m那么每帧移动量是10*dt/delta个像素dt是帧时间间隔。另一个常用操作是把相位屏分解成 Zernike 模式用于分析湍流的低阶像差分布。aotools 的zernike模块可以生成模式图在离散采样下做最小二乘拟合。例如生成离焦模式from aotools.functions import zernike # zernike_nm(n, m, size) 生成单个模式 defocus_map zernike.zernike_nm(2, 0, 256)如果你想拟合相位屏的前 15 阶模式系数从源码角度看比较稳的是把各模式图铺成矩阵再对相位屏做最小二乘。注意 Zernike 模式在像素网格上并不完全正交拟合前建议先做一次 Gram-Schmidt 正交化否则拟合系数会和 Noll 理论值对不上。这也是为什么仓库里单独有个test_KL.py——Karhunen-Loève 模式从协方差矩阵特征分解来比硬套 Zernike 更贴近真实湍流统计。实验场景里我一般用 Zernike 做控制基底用 KL 做物理统计。Zernike 阶名称物理含义对成像影响2离焦轴向焦点偏移整体模糊3像散两个方向焦点不同星点拉伸4彗差非对称像差拖尾有方向性5三叶草低频不规则星点三角扩展3. 哈特曼波前传感器从子孔径图像到斜率3.1 模拟哈特曼子孔径的完整图像哈特曼传感器在光路上加一块微透镜阵列把入射波前切成几十个子孔径每个子孔径在探测器上形成一个光斑。波前倾斜会让光斑从子孔径中心偏移偏移量正比于该子孔径内的平均斜率。aotools 的wfs模块提供算法部分但光学追迹要自己拼。我的做法是先用 numpy 构造一个 8×8 子孔径的光斑图再用质心函数提取偏移量def make_spot_array(subap_size16, nx8): img np.zeros((subap_size*nx, subap_size*nx)) y, x np.mgrid[0:subap_size*nx, 0:subap_size*nx] for i in range(nx): for j in range(nx): cy i*subap_size subap_size/2 0.3*(i1) cx j*subap_size subap_size/2 0.2*(j1) img np.exp(-((x - cx)**2 (y - cy)**2) / (2*2.5**2)) return img这里每个光斑都是高斯形状cx、cy里加了随 i、j 变化的偏移来模拟重力方向像差和安装误差。实际仿真中可以换成随机倾斜每个子孔径的偏移量从高斯分布里采样方差由该子孔径对应的大气斜率决定。注意subap_size必须大于光斑直径的 2 倍否则光斑被窗口截断质心会系统性偏向子孔径中心。3.2 质心算法别直接用灰度峰值很多从图像处理转过来的同事会下意识找光斑最大值位置这在哈特曼场景是不对的。峰值像素受衍射环和探测器噪声干扰很大正确做法是计算整个光斑的能量重心也就是一阶矩。aotools 提供了centre_of_gravity和thresholded_centre_of_gravity后者会先扣掉低于阈值的能量抑制暗电流from aotools.functions import centroiders def measure_slopes(img, subap_size, nx, thresholdNone): slopes_x np.zeros((nx, nx)) slopes_y np.zeros((nx, nx)) for i in range(nx): for j in range(nx): sub img[i*subap_size:(i1)*subap_size, j*subap_size:(j1)*subap_size] cx, cy centroiders.thresholded_centre_of_gravity( sub, thresholdthreshold or sub.max()*0.3 ) slopes_x[i, j] cx - subap_size/2 slopes_y[i, j] cy - subap_size/2 return slopes_x, slopes_ythreshold取sub.max()*0.3是我常用的经验值对高斯光斑能滤掉读出噪声又保留足够能量。如果光斑太暗阈值偏高会把重心拉偏阈值太低则引入噪声。aotools 的test_centroiders.py里有现成对比可以跑一遍看不同算法在特定噪声下的 RMS 误差比自己拍脑袋选算法可靠得多。质心误差有一个经典近似sigma_c sigma_spot * sigma_noise / (N_photons * sigma_spot)所以低信噪比时增大光斑尺寸反而会让质心抖动变小直到光子噪声占主导。这块调参是哈特曼仿真里最耗时间的部分。算法鲁棒性适用场景重心法 CoG中等受噪声偏置影响高信噪比快读阈值重心法较好能压暗电流中低 SNR 实验平方加权重心强抑制旁瓣多光斑干扰3.3 从斜率重建波前Zernike 导数矩阵有了子孔径斜率之后下一步是把斜率场重建为波前相位。常见做法是把波前展开成 Zernike 多项式然后求多项式对 x、y 的偏导数再用子孔径平均斜率去拟合。关键在构造矩阵A它的每列对应一个模式每两行分别对应 x、y 方向的斜率from numpy.linalg import lstsq def zernike_slope_matrix(n_modes, subap_grid, radius): n_rows subap_grid.size * 2 A np.zeros((n_rows, n_modes)) for mode in range(n_modes): phase_map zernike.zernike_nm(*noll_index[mode1], sizeradius) dy, dx np.gradient(phase_map) A[0::2, mode] average_in_subaps(dx, subap_grid) A[1::2, mode] average_in_subaps(dy, subap_grid) return Aaverage_in_subaps把每个子孔径区域内的像素平均得到该子孔径的平均斜率。构造完A后用coeffs, _, _, _ lstsq(A, slopes, rcondNone)拟合。注意slopes的排布顺序必须和A的行顺序一致否则拟合出的波前会忽高忽低。我见过至少三次因为 x、y 斜率交替顺序写反而拟合出奇怪波前的坑。子孔阵列可稳定重建的模式数说明4×46~10只适合低阶像差8×820~40常见实验配置16×1660~120注意矩阵病态需截断奇异值模式数不要超过子孔径数量的 1/2否则矩阵病态系数开始震荡。判断模式数是否合适的标准是看拟合残差分布残差集中在边缘说明模式基不够均匀分布则说明模式数已经够了。4. 变形镜控制把哈特曼斜率变成可收敛的闭环4.1 先有影响矩阵再有控制信号变形镜在 aotools 里没有显式类但整个自适应光学闭环可以拆成模式基底来组合。大多数变形镜可以近似为一组高斯响应驱动器每个驱动器在镜面上形成一个鼓包。对这个响应函数做采样就得到驱动器的影响矩阵M。控制目标是把波前拟合误差分配到各个驱动器上def dm_influence(x, y, cx, cy, coupling0.2): r np.hypot(x - cx, y - cy) return np.exp(-(r / (0.3 * coupling))**2)coupling是驱动器间距与高斯宽度的比值实际 DM 的耦合度通常在 0.1~0.4 之间。仿真时把驱动器排成 8×8 网格影响矩阵M的每一列是一个驱动器对子孔径斜率的响应。控制回路里先测量当前残余波前的斜率再用影响矩阵的伪逆生成驱动信号。4.2 串起来的闭环仿真把连续相位屏、哈特曼斜率测量和 DM 响应组合在一起就是最简的自适应光学闭环phase_screens inf_phase_screen(...) # 连续大气相位序列 dm np.zeros((n_act, 1)) # 初始驱动电压 gain 0.4 for frame in phase_screens: residual frame - apply_dm(dm) # 残余波前 slopes measure(residual) # 哈特曼测斜率 coeffs reconstruct(slopes) # Zernike 拟合 dm dm gain * M_pinv coeffs # 积分控制这里的M_pinv是把影响矩阵和 Zernike 模式耦合后的伪逆不是直接对斜率矩阵求逆。积分增益gain取 0.3~0.6 通常比较稳妥太大回路会振荡太小收敛慢。判断是否收敛不要只看一帧残余波前的 RMS要看它随时间是否振荡。把每帧residual的 RMS 存下来画曲线如果 RMS 下降到平稳且没有周期起伏就说明控制带宽够了。如果想定量评估闭环带宽可以给相位屏加一个频率逐渐升高的正弦倾斜模式观察残余 RMS 从小于 0.1 rad 变成大于 1 rad 的转折点该频率就是系统的 -3 dB 控制带宽。对于 8×8 子孔径、500 帧/秒的仿真配置积分增益 0.4 通常能到 4~6 Hz。带宽受斜率计算延迟影响很大质心算法越复杂延迟越高这也是为什么工程上倾向用简单重心法做实时控制把更耗时的降噪留给离线处理。aotools 里的test_slopecovariance.py会校验斜率协方差与相位结构函数的关系这在调L0和子孔径大小时很有用。控制带宽不足时残余误差主要来自高频分量此时能在斜率协方差图里看到明显的截止效应帮助判断瓶颈在探测器帧率还是在控制环路增益。5. 把 test_*.py 当排错手册aotools 的高效调参方式这份源码最有价值的其实不在aotools目录而在test/下那一堆test_*.py。很多库把测试当摆设但 aotools 的测试既是回归检查也是浓缩的 API 文档test_infinitephasescreen.py告诉你连续相位屏怎么调参test_wfs.py告诉你各波前传感算法输入输出的 shape 到底是什么。最直接的用法是跑一遍 pytest先建立置信区间cd aotools python -m pytest test/ -k turbulence or wfs -v如果测试失败多数情况不是逻辑错而是环境版本问题。aotools 的setup.py依赖 numpy 和 scipyPython 3.8 以上搭配 numpy 1.20 左右表现稳定。如果你发现ft_phase_screen生成的相位屏有明显网格状条纹先检查是不是 FFT 归一化被环境里的 pyfftw 干扰或者delta取了整数导致采样频率恰好与网格重合。另一个排错技巧是把test_wfs.py里的测试图像直接打印出来看光斑有没有被截断。子孔径边界如果切到光斑一半质心计算会整体偏向网格中心表现为闭环残余波前始终有一个固定方向分量。这种系统性误差靠滤波消不掉只能回头把subap_size加大到光斑直径的 2 倍以上。调试变形镜参数时我会把增益、模式数、子孔径数做成三张趋势表固定两个变一个看 RMS 收敛曲线。跑一组参数只需几百帧12×12 子孔径在 aotools 下大约几十毫秒一帧完全可以做参数扫描。最容易被忽略的是 DM 驱动器布局要和子孔径网格对齐错开半个网格就会引入无法解释的高频误差。最后分享一个常用技巧做像差分解时把 Zernike 模式和 Karhunen-Loève 模式都算一遍如果两者给出的低阶系数差超过 10%说明相位屏采样不足先加大 N 再谈控制优化。本文还有配套的精品资源点击获取