简介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( nxnx, pixel_scalepixel_scale, r0r0, L0L0, ) phase ps.scrn print(f相位 RMS: {phase.std():.3f} rad)代码逻辑并不复杂构造PhaseScreenKolmogorov时内部完成频谱初始化ps.scrn取出当前相位屏返回值是一个nx × nx的二维数组单位是弧度。打印出来的 RMS 在几弧度到十几弧度的量级是正常的如果发现 RMS 接近 0.x 甚至 1e-3大概率这个版本返回的是 waves 而不是 radians需要自行除以波长再乘 2π 换算。参数配置是这里最容易失控的地方几个关键参数整理出来参数示例值作用与调整方向nx256空间分辨率越高高频越完整FFT 成本按 n² log n 增长pixel_scale0.01 m每像素代表实际尺度应保持在 r0/2 到 r0 之间r00.15 m湍流强度越小相位起伏越大L025 m外尺度截断低频发散几十米量级即可提示纯 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( nxnx, pixel_scalepixel_scale, r0r0, L0L0, wind_speedwind_speed, wind_directionwind_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_subap8): 把相位屏切成 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 成四维前两维是子孔径行列后两维是子孔径内部像素然后对每个子块做中心差分梯度相当于真实系统中光斑质心偏移最后把子孔径内的梯度平均成一个斜率值对应探测器上一个光斑的净偏移。真实系统的质心算法在这里被梯度平均替代物理含义一致。子孔径数的选择要跟变形镜驱动器数和采样分辨率匹配子孔径布局每个子孔径像素适用场景4×464×64低阶像差为主驱动器数少8×832×32和 8×8 变形镜配套最常用16×1616×16空间采样细但斜率噪声上升明显注意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_modes6): A np.stack([zernike_basis(m, N).ravel() for m in range(1, n_modes 1)], axis1) coeff, *_ np.linalg.lstsq(A, phase.ravel(), rcondNone) 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 surfact_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.0H 的形状是 128×6464 个驱动器各标定一路每路产生 128 个斜率输出8×8 子孔径 × x/y 两个方向。维度对上了后面的控制矩阵才能算。这个标定过程在仿真里只是循环 64 次但在真实系统里是每次实验前必须执行的标定流程代码结构完全一致。4.4 控制矩阵伪逆、截断与正则化控制矩阵理论上就是响应矩阵的伪逆control np.linalg.pinv(H)直接求伪逆会在边缘驱动器上出问题。四角驱动器的有些影响函数落在有效子孔径外H 中对应列的能量很小伪逆会给这些列很大的增益测量噪声被放大成很高的电压。常用的处理是截断奇异值或者加 Tikhonov 正则化# 截断伪逆小于最大奇异值 rcond 倍的奇异值直接丢弃 control np.linalg.pinv(H, rcond1e-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)两种方案的参数对闭环行为影响很大具体表现如下方案参数调小了调大了截断伪逆rcond1e-3边缘增益放大噪声明显校正能力下降残差变大Tikhonovλ1.0接近普通伪逆校正变钝低阶像差不干净rcond 和 λ 建议在闭环之前扫一遍以闭环残差方差最小为准则。常见做法是从 rcond1e-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.90.5 rad² 时掉到 0.6 附近。如果闭环残差比开环还大不要急着调增益先回头检查响应矩阵维度、单位换算和斜率排列顺序。把这套开环/闭环对比脚本固定成回归任务之后每次改子孔径数、驱动器数、r0 或风速重跑一遍并对比曲线曲线偏离基线就说明这次改动引入了问题。本文还有配套的精品资源点击获取 aotools大气湍流仿真与自适应光学闭环控制实践：从相位屏到变形镜
简介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( nxnx, pixel_scalepixel_scale, r0r0, L0L0, ) phase ps.scrn print(f相位 RMS: {phase.std():.3f} rad)代码逻辑并不复杂构造PhaseScreenKolmogorov时内部完成频谱初始化ps.scrn取出当前相位屏返回值是一个nx × nx的二维数组单位是弧度。打印出来的 RMS 在几弧度到十几弧度的量级是正常的如果发现 RMS 接近 0.x 甚至 1e-3大概率这个版本返回的是 waves 而不是 radians需要自行除以波长再乘 2π 换算。参数配置是这里最容易失控的地方几个关键参数整理出来参数示例值作用与调整方向nx256空间分辨率越高高频越完整FFT 成本按 n² log n 增长pixel_scale0.01 m每像素代表实际尺度应保持在 r0/2 到 r0 之间r00.15 m湍流强度越小相位起伏越大L025 m外尺度截断低频发散几十米量级即可提示纯 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( nxnx, pixel_scalepixel_scale, r0r0, L0L0, wind_speedwind_speed, wind_directionwind_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_subap8): 把相位屏切成 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 成四维前两维是子孔径行列后两维是子孔径内部像素然后对每个子块做中心差分梯度相当于真实系统中光斑质心偏移最后把子孔径内的梯度平均成一个斜率值对应探测器上一个光斑的净偏移。真实系统的质心算法在这里被梯度平均替代物理含义一致。子孔径数的选择要跟变形镜驱动器数和采样分辨率匹配子孔径布局每个子孔径像素适用场景4×464×64低阶像差为主驱动器数少8×832×32和 8×8 变形镜配套最常用16×1616×16空间采样细但斜率噪声上升明显注意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_modes6): A np.stack([zernike_basis(m, N).ravel() for m in range(1, n_modes 1)], axis1) coeff, *_ np.linalg.lstsq(A, phase.ravel(), rcondNone) 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 surfact_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.0H 的形状是 128×6464 个驱动器各标定一路每路产生 128 个斜率输出8×8 子孔径 × x/y 两个方向。维度对上了后面的控制矩阵才能算。这个标定过程在仿真里只是循环 64 次但在真实系统里是每次实验前必须执行的标定流程代码结构完全一致。4.4 控制矩阵伪逆、截断与正则化控制矩阵理论上就是响应矩阵的伪逆control np.linalg.pinv(H)直接求伪逆会在边缘驱动器上出问题。四角驱动器的有些影响函数落在有效子孔径外H 中对应列的能量很小伪逆会给这些列很大的增益测量噪声被放大成很高的电压。常用的处理是截断奇异值或者加 Tikhonov 正则化# 截断伪逆小于最大奇异值 rcond 倍的奇异值直接丢弃 control np.linalg.pinv(H, rcond1e-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)两种方案的参数对闭环行为影响很大具体表现如下方案参数调小了调大了截断伪逆rcond1e-3边缘增益放大噪声明显校正能力下降残差变大Tikhonovλ1.0接近普通伪逆校正变钝低阶像差不干净rcond 和 λ 建议在闭环之前扫一遍以闭环残差方差最小为准则。常见做法是从 rcond1e-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.90.5 rad² 时掉到 0.6 附近。如果闭环残差比开环还大不要急着调增益先回头检查响应矩阵维度、单位换算和斜率排列顺序。把这套开环/闭环对比脚本固定成回归任务之后每次改子孔径数、驱动器数、r0 或风速重跑一遍并对比曲线曲线偏离基线就说明这次改动引入了问题。本文还有配套的精品资源点击获取