简介面向自适应光学与波前校正方向的开发者和研究人员这套Python仿真资源围绕soapy框架提供了大气湍流模拟、波前探测、变形镜建模与控制算法的完整代码实现可用于天文观测、激光传输等场景的算法验证与教学科研。包体内共429个文件以280个Python源码为核心辅以yaml参数配置、ui界面、rst文档、fits数据及bat启动脚本压缩包仅1.39MB结构紧凑且易于部署。目前已有97人学习适合具备Python和NumPy基础、希望深入掌握自适应光学核心算法的中高级读者。通过这套代码读者可以理解分层相位屏湍流模拟、Zernike模式分解与波前重构、卡尔曼预测控制、压电驱动器影响函数建模等关键环节并利用Strehl比、点扩散函数等诊断模块评估校正效果模块化的代码组织和配置文件也便于按需替换传感器或校正器模型为二次开发与科研复现提供了良好基础。1. 自适应光学仿真与波前校正Python能把什么跑通自适应光学Adaptive Optics, AO用于实时补偿大气湍流或热形变带来的波前畸变。一个完整的AO系统包含波前传感器、波前重建器、控制器和变形镜每部分都能用线性代数模型描述。在建造真实硬件之前用Python把这条链路在数值上跑通能显著降低选型、控制器调试和噪声鲁棒性验证的成本。这套仿真不需要光学专用软件只用NumPy和SciPy就能生成湍流相位屏、模拟Shack-Hartmann采样、做波前重建并完成闭环校正。代码层面关键是把连续光场离散成二维相位网格把传感器取样建模为对梯度的测量把波前重建看成一个线性反问题。下文按“相位表示—传感采样—重建—闭环控制”给出可运行的实现片段最后讨论参数标定和稳定性。具备信号处理或光学背景的Python使用者可以按步骤复现再扩展到自己关心的场景。2. 波前相位与湍流屏的Python表示AO仿真第一步是生成待校正的波前误差。相位屏可以用二维数组表示数组中每个元素是某个网格点上的光程差单位通常取弧度。生成方式分成两类一个是叠加Zernike多项式获得低频像差另一个是用FFT频谱滤波法生成符合Kolmogorov统计的湍流屏。前者适合做确定性测试后者适合做统计性能评估。2.1 用Zernike多项式构造平滑相位畸变Zernike多项式是单位圆上正交的一组基函数低阶项对应离焦、像散、慧差等经典像差。实现时把系数乘以基函数再叠加就能得到连续的平滑波前。下面是生成三种常见像差的函数坐标范围按单位圆归一化import numpy as np def zernike_map(n, m, X, Y): r np.sqrt(X**2 Y**2) theta np.arctan2(Y, X) if n 2 and m 0: return 2*r**2 - 1 # 离焦 elif n 2 and m 2: return r**2 * np.cos(2*theta) # 0°像散 elif n 2 and m -2: return r**2 * np.sin(2*theta) # 45°像散 elif n 3 and m 1: return (3*r**3 - 2*r) * np.cos(theta) # 慧差 elif n 3 and m -1: return (3*r**3 - 2*r) * np.sin(theta) else: return np.zeros_like(X) # 未实现项先置零 N 64 x np.linspace(-1, 1, N) X, Y np.meshgrid(x, x) phase 0.8 * zernike_map(2, 0, X, Y) 0.3 * zernike_map(2, 2, X, Y)这段代码用meshgrid生成二维坐标随后按标准Zernike径向多项式计算每个像差面型。离焦项系数0.8弧度、像散项系数0.3弧度叠加后波前峰谷约1.1弧度大约是六分之一波长对中等强度畸变来说是比较合理的起点。注意这里的径向多项式没有乘以角向归一化因子因为后续重建和校正关心的是相对形状与系数关系绝对RMS可以通过系数换算但在仿真中更常用的是直接用相位数组计算统计量。若要生成任意Noll序号的多项式建议引入现成的zernike库自己写容易在归一化因子上出错。2.2 用FFT频谱滤波法生成Kolmogorov湍流屏真实大气湍流相位屏的功率谱密度服从f^(-11/3)。频谱滤波法的思路是生成一个复高斯随机矩阵把它乘上功率谱的平方根再做二维逆傅里叶变换得到随机相位屏。实现代码不长但有几个细节影响结果def kolmogorov_phase_screen(N, L, r0, seed42): rng np.random.default_rng(seed) dx L / N fx np.fft.fftfreq(N, ddx) FX, FY np.meshgrid(fx, fx) f np.hypot(FX, FY) f[0, 0] 1e-6 # 避免零频除零 spectrum 0.023 * (r0 ** (-5.0/3.0)) * np.power(f, -11.0/3.0) gauss rng.normal(size(N, N)) 1j * rng.normal(size(N, N)) screen np.fft.ifft2(np.sqrt(spectrum) * gauss).real screen * N # 补偿FFT离散归一化 return screen - screen.mean()代码中fftfreq生成空间频率坐标单位是1/mspectrum的系数0.023对应经典Kolmogorov湍流功率谱。gauss的实部和虚部都是独立标准正态分布保证相位屏的傅里叶分量有随机振幅和相位。乘以N这一步是为了把离散傅里叶变换的缩放因子抵消实际执行后相位屏的幅度才能落在合理范围。生成后会看到明显的低频大尺度起伏这是湍流能谱中低频能量集中的表现。若要看结构函数是否符合理论值需要保留多帧做统计平均直接观察单帧会以为随机性不足。r0 (cm)湍流强度典型场景10强近地面、白天强热对流20中天文台夜视常见条件50弱高海拔优良台址使用时应根据口径和等效高度选择L和r0。L表示相位屏的空间尺寸单位米工程上取口径的1.2倍以上。r0越小湍流越强生成的相位屏梯度和振幅越大。仿真时如果发现后面闭环控制无法收敛往往不是控制算法出了问题而是r0过小导致相位屏局部梯度超出传感器采样范围。这时应先减小r0的仿真强度或者加密子孔径采样再回头查控制参数。3. 波前传感与斜率计算的Python仿真波前传感器的常见硬件是Shack-Hartmann它用微透镜阵列把波前分割成许多子孔径每个子孔径的焦点在探测器上形成一个光斑。光斑相对参考位置的偏移量与子孔径内平均波前斜率成正比。仿真这一步的关键是把相位网格映射到子孔径斜率测量值并叠加探测器噪声。3.1 Shack-Hartmann子孔径平均斜率计算假设相位屏大小为N×N布置N_sub×N_sub个子孔径每个子孔径覆盖sub_size×sub_size个网格点。子孔径的平均斜率可以用该区域相位梯度的平均值表示。用NumPy的gradient函数就能快速实现def sh_slopes(phase, N_sub, sub_size): N phase.shape[0] assert N N_sub * sub_size, 尺寸不匹配 slopes_x np.zeros((N_sub, N_sub)) slopes_y np.zeros((N_sub, N_sub)) for i in range(N_sub): for j in range(N_sub): block phase[i*sub_size:(i1)*sub_size, j*sub_size:(j1)*sub_size] dy, dx np.gradient(block) slopes_x[i, j] dx.mean() slopes_y[i, j] dy.mean() return slopes_x, slopes_y sx, sy sh_slopes(phase, N_sub8, sub_size8)这里的dx和dy分别表示相位在x和y方向上的差分单位是弧度/网格。由于子孔径内取平均局部噪声会被平滑。N_sub决定波前传感器的空间分辨率子孔径数越多能采样的高频信息越丰富但每个孔径收集的光子数减少信噪比会下降。反过来子孔径太少则斜率信息不足重建出来的波前会忽略高阶像差。实际系统里常用8×8或16×16仿真可以按需求赋值。3.2 探测器像素化与噪声模型真实探测器的光斑位置估计受读出噪声和光子噪声影响因此斜率测量值是带有误差的。仿真时可以在理想斜率上叠加高斯噪声噪声尺度用弧度表示def add_measurement_noise(slopes, sigma0.01): noise np.random.normal(0, sigma, sizeslopes.shape) return slopes noise sx_noisy add_measurement_noise(sx, sigma0.02) sy_noisy add_measurement_noise(sy, sigma0.02)这里的sigma需要结合探测器性能选择。一个经验值是取子孔径像素衍射极限的1/10比如工作波长632nm、子孔径直径2mm时光斑质心误差约0.5μm换算到斜率约0.4μrad。但如果仿真网格尺寸和波长没有严格绑定sigma直接按弧度设成0.010.05即可。注意sigma过大会让后面重建矩阵小型特征值被噪声主导导致重建结果出现棋盘格状抖动。这是AO仿真里最典型的“噪声-分辨率矛盾”后面一章会看到正则化如何缓解。4. 波前重建从斜率到相位网格传感器给出的是斜率数组而不是相位。重建器要把斜率映射回二维相位分布。这里采用Southwell网格重建它假设相位采样点位于子孔径中心x方向相邻相位点的差值等于该方向斜率乘以子孔径间距dy方向同理。把所有方程写成线性系统用最小二乘求解。4.1 Southwell最小二乘重建矩阵设子孔径网格为m×m相位向量长度为m²。对每个水平相邻点对有m×(m-1)个方程对每个垂直相邻点对有(m-1)×m个方程。再增加一个活塞零约束保证方程有唯一解。实现如下def southwell_reconstruct(sx, sy, d): m sx.shape[0] n m * m A [] b [] # x方向约束 for i in range(m): for j in range(m-1): row np.zeros(n) row[i*m j] -1 row[i*m j1] 1 A.append(row) b.append(sx[i, j] * d) # y方向约束 for i in range(m-1): for j in range(m): row np.zeros(n) row[i*m j] -1 row[(i1)*m j] 1 A.append(row) b.append(sy[i, j] * d) A np.array(A) b np.array(b) # 固定第一个相位值为0消除活塞平移自由度 A np.vstack([A, np.zeros(n)]) A[-1, 0] 1 b np.append(b, 0) phi, _, _, _ np.linalg.lstsq(A, b, rcondNone) return phi.reshape(m, m) reconstructed southwell_reconstruct(sx_noisy, sy_noisy, d1.0)代码先构造系数矩阵A每一行代表一个差分方程。x方向约束里-1和1分别对应左右相邻相位点右侧等于斜率乘以间距d。y方向类似。最后加的零约束行用于去掉全活塞项否则A不满秩最小二乘结果会带上任意常数偏移。使用lstsq而不是直接求逆是因为该线性系统是超定的最小二乘解在数值上更稳定。若子孔径数增加到16×16矩阵规模约为512×256计算瞬时完成但换成真实硬件在线运行时需要提前把A的伪逆算好避免每帧都做SVD。4.2 正则化与重建误差评估当斜率噪声sigma较大时直接最小二乘重建的波前会高频振荡。一个常见改进是给最小二乘增加Tikhonov正则项将目标变为min ||Aφ - b||² λ||Lφ||²其中L是二阶差分算子。在代码里可以简单地在反演前给A的协方差矩阵加一个小对角阵AtA A.T A Atb A.T b recon np.linalg.solve(AtA 0.01 * np.eye(n), Atb).reshape(m, m)这里的0.01就是正则化系数它等价于假设噪声方差与信号方差之比为0.01。系数过小无法压制噪声过大会把真实像差抹平。调试方法是用一组已知相位屏生成斜率并加噪声然后改变λ看重建误差RMS曲线选择曲线拐点处的λ。在工程上这一步往往比重建算法本身更影响性能因为真实系统的斜率噪声不可能为零。5. 变形镜响应与控制闭环变形镜Deformable Mirror, DM是校正阶段的执行器通过施加电压使镜面产生面形变化。在仿真里DM可以用一组影响函数的线性组合表示。每个驱动器的电压变化会在镜面上产生一个近似高斯形状的位移多个驱动器叠加得到总相位校正量。控制器根据重建出的残留相位计算电压增量形成闭环。5.1 影响函数与电压到相位的响应矩阵影响函数定义为一个驱动器在镜面上产生的相位分布通常用高斯函数模拟。把驱动器放在均匀网格上每个驱动器占据一个列向量所有列拼成响应矩阵R。生成R的代码如下def build_dm_response(N, drive_positions, width): x np.linspace(-1, 1, N) X, Y np.meshgrid(x, x) R np.zeros((N*N, len(drive_positions))) for k, (xc, yc) in enumerate(drive_positions): r2 (X - xc)**2 (Y - yc)**2 R[:, k] np.exp(-r2 / (width**2)).ravel() return R # 假设 5x5 驱动器间距0.5高斯宽度0.8 gx, gy np.meshgrid(np.linspace(-1, 1, 5), np.linspace(-1, 1, 5)) pos list(zip(gx.ravel(), gy.ravel())) R build_dm_response(N, pos, width0.8)注意影响函数宽度要设置成驱动器间距的0.61.2倍。宽度太小会导致镜面产生波浪形起伏驱动器之间无法平滑过渡宽度太大则相邻驱动器耦合格外强烈矩阵R条件数变大求逆时会把噪声放大。生成R后可以用np.linalg.pinv(R)计算电压到相位的逆映射但实际控制器中更常用的是预先计算控制矩阵把重建相位直接映射到电压增量。5.2 闭环积分控制与仿真发散经典AO闭环是纯积分控制当前控制电压 上一步电压 增益 × 校正电压。仿真每一帧先计算当前残余相位下的斜率重建出相位再乘以响应矩阵的伪逆得到电压增量最后更新残余相位。实现示意如下def closed_loop(init_phase, R, R_inv, N_sub, sub_size, gain0.5, n_iter30): residual init_phase.copy() V np.zeros(R.shape[1]) history [] for k in range(n_iter): sx, sy sh_slopes(residual, N_sub, sub_size) sx np.random.normal(0, 0.01, sx.shape) sy np.random.normal(0, 0.01, sy.shape) phase_rec southwell_reconstruct(sx, sy, d1.0) V_new gain * (R_inv phase_rec.ravel()) V V V_new correction (R V).reshape(init_phase.shape) residual init_phase - correction history.append(np.sqrt(np.mean((residual - residual.mean())**2))) return residual, history res, hist closed_loop(phase, R, np.linalg.pinv(R), 8, 8, gain0.4)这段代码模仿真实CCD采样每次循环都对斜率加噪声。gain是环路增益典型范围0.10.8。增益过大一步校正量超过当前误差会来回震荡并出现仿真发散表现为RMSE曲线先下降后迅速上升增益过小收敛速度慢在动态湍流下永远追不上畸变变化。另一个常见发散原因是phase_rec重建时未消除活塞项导致校正量中混入整体平移不过这不会影响像差RMS。若运行后发现residual的RMS在第10步后反弹优先把gain降到0.3以下再检查噪声sigma是否过大。仿真发散不是代码bug而是采样率、噪声和增益三者失配的物理后果需要在参数层面做平衡。6. 优化验证技巧与在线扩展闭环仿真跑通后下一步是让它给出可信的性能数字而不是只输出一堆二维图。这里给出三个实用技巧能让仿真结果更接近实验室测量。第一个技巧是用相位RMS和Strehl比做单帧评估。假设工作波长与相位屏长度单位一致残余波前方差σ²可以通过np.var(residual)计算Strehl近似用exp(-σ²)。RMS小于0.1弧度时系统已接近衍射极限此时再提升增益收益不大反而容易引起振荡。第二个技巧是扫描增益和正则化系数画出收敛曲线热力图。把两个参数各取10个值运行同一组湍流屏统计收敛后10帧的RMS平均值用matplotlib.imshow画成二维图能看到明显的稳定区域和发散区域。性能指标直接决定后续硬件选型比如控制器带宽要求由湍流截止频率决定而传感器帧率决定了能校正多大的r0。第三个技巧是离线相位屏库的复用。如果每次运行都重新生成随机屏不同参数之间的性能对比不公平。较好做法是保存一批固定种子的屏先用np.save存成.npy文件对比算法时从文件重复读取。这也能让复现实验变得容易审计时只需要提供文件和数据版本。最后把仿真算法从离线阶段迁移到实时系统时注意避免两处性能陷阱一是在每帧循环里重复构造Southwell矩阵应预先计算伪逆运行时只做两次矩阵乘二是scipy.sparse.linalg.lsqr比np.linalg.lstsq快很多在大子孔径数和多帧迭代下差距可以达到十倍。优化后的循环里单帧处理时间可以在普通笔记本上压缩到2毫秒以内这足以支撑数百赫兹的控制带宽模拟。本文还有配套的精品资源点击获取
Python实现自适应光学仿真：波前重建与闭环校正详解
简介面向自适应光学与波前校正方向的开发者和研究人员这套Python仿真资源围绕soapy框架提供了大气湍流模拟、波前探测、变形镜建模与控制算法的完整代码实现可用于天文观测、激光传输等场景的算法验证与教学科研。包体内共429个文件以280个Python源码为核心辅以yaml参数配置、ui界面、rst文档、fits数据及bat启动脚本压缩包仅1.39MB结构紧凑且易于部署。目前已有97人学习适合具备Python和NumPy基础、希望深入掌握自适应光学核心算法的中高级读者。通过这套代码读者可以理解分层相位屏湍流模拟、Zernike模式分解与波前重构、卡尔曼预测控制、压电驱动器影响函数建模等关键环节并利用Strehl比、点扩散函数等诊断模块评估校正效果模块化的代码组织和配置文件也便于按需替换传感器或校正器模型为二次开发与科研复现提供了良好基础。1. 自适应光学仿真与波前校正Python能把什么跑通自适应光学Adaptive Optics, AO用于实时补偿大气湍流或热形变带来的波前畸变。一个完整的AO系统包含波前传感器、波前重建器、控制器和变形镜每部分都能用线性代数模型描述。在建造真实硬件之前用Python把这条链路在数值上跑通能显著降低选型、控制器调试和噪声鲁棒性验证的成本。这套仿真不需要光学专用软件只用NumPy和SciPy就能生成湍流相位屏、模拟Shack-Hartmann采样、做波前重建并完成闭环校正。代码层面关键是把连续光场离散成二维相位网格把传感器取样建模为对梯度的测量把波前重建看成一个线性反问题。下文按“相位表示—传感采样—重建—闭环控制”给出可运行的实现片段最后讨论参数标定和稳定性。具备信号处理或光学背景的Python使用者可以按步骤复现再扩展到自己关心的场景。2. 波前相位与湍流屏的Python表示AO仿真第一步是生成待校正的波前误差。相位屏可以用二维数组表示数组中每个元素是某个网格点上的光程差单位通常取弧度。生成方式分成两类一个是叠加Zernike多项式获得低频像差另一个是用FFT频谱滤波法生成符合Kolmogorov统计的湍流屏。前者适合做确定性测试后者适合做统计性能评估。2.1 用Zernike多项式构造平滑相位畸变Zernike多项式是单位圆上正交的一组基函数低阶项对应离焦、像散、慧差等经典像差。实现时把系数乘以基函数再叠加就能得到连续的平滑波前。下面是生成三种常见像差的函数坐标范围按单位圆归一化import numpy as np def zernike_map(n, m, X, Y): r np.sqrt(X**2 Y**2) theta np.arctan2(Y, X) if n 2 and m 0: return 2*r**2 - 1 # 离焦 elif n 2 and m 2: return r**2 * np.cos(2*theta) # 0°像散 elif n 2 and m -2: return r**2 * np.sin(2*theta) # 45°像散 elif n 3 and m 1: return (3*r**3 - 2*r) * np.cos(theta) # 慧差 elif n 3 and m -1: return (3*r**3 - 2*r) * np.sin(theta) else: return np.zeros_like(X) # 未实现项先置零 N 64 x np.linspace(-1, 1, N) X, Y np.meshgrid(x, x) phase 0.8 * zernike_map(2, 0, X, Y) 0.3 * zernike_map(2, 2, X, Y)这段代码用meshgrid生成二维坐标随后按标准Zernike径向多项式计算每个像差面型。离焦项系数0.8弧度、像散项系数0.3弧度叠加后波前峰谷约1.1弧度大约是六分之一波长对中等强度畸变来说是比较合理的起点。注意这里的径向多项式没有乘以角向归一化因子因为后续重建和校正关心的是相对形状与系数关系绝对RMS可以通过系数换算但在仿真中更常用的是直接用相位数组计算统计量。若要生成任意Noll序号的多项式建议引入现成的zernike库自己写容易在归一化因子上出错。2.2 用FFT频谱滤波法生成Kolmogorov湍流屏真实大气湍流相位屏的功率谱密度服从f^(-11/3)。频谱滤波法的思路是生成一个复高斯随机矩阵把它乘上功率谱的平方根再做二维逆傅里叶变换得到随机相位屏。实现代码不长但有几个细节影响结果def kolmogorov_phase_screen(N, L, r0, seed42): rng np.random.default_rng(seed) dx L / N fx np.fft.fftfreq(N, ddx) FX, FY np.meshgrid(fx, fx) f np.hypot(FX, FY) f[0, 0] 1e-6 # 避免零频除零 spectrum 0.023 * (r0 ** (-5.0/3.0)) * np.power(f, -11.0/3.0) gauss rng.normal(size(N, N)) 1j * rng.normal(size(N, N)) screen np.fft.ifft2(np.sqrt(spectrum) * gauss).real screen * N # 补偿FFT离散归一化 return screen - screen.mean()代码中fftfreq生成空间频率坐标单位是1/mspectrum的系数0.023对应经典Kolmogorov湍流功率谱。gauss的实部和虚部都是独立标准正态分布保证相位屏的傅里叶分量有随机振幅和相位。乘以N这一步是为了把离散傅里叶变换的缩放因子抵消实际执行后相位屏的幅度才能落在合理范围。生成后会看到明显的低频大尺度起伏这是湍流能谱中低频能量集中的表现。若要看结构函数是否符合理论值需要保留多帧做统计平均直接观察单帧会以为随机性不足。r0 (cm)湍流强度典型场景10强近地面、白天强热对流20中天文台夜视常见条件50弱高海拔优良台址使用时应根据口径和等效高度选择L和r0。L表示相位屏的空间尺寸单位米工程上取口径的1.2倍以上。r0越小湍流越强生成的相位屏梯度和振幅越大。仿真时如果发现后面闭环控制无法收敛往往不是控制算法出了问题而是r0过小导致相位屏局部梯度超出传感器采样范围。这时应先减小r0的仿真强度或者加密子孔径采样再回头查控制参数。3. 波前传感与斜率计算的Python仿真波前传感器的常见硬件是Shack-Hartmann它用微透镜阵列把波前分割成许多子孔径每个子孔径的焦点在探测器上形成一个光斑。光斑相对参考位置的偏移量与子孔径内平均波前斜率成正比。仿真这一步的关键是把相位网格映射到子孔径斜率测量值并叠加探测器噪声。3.1 Shack-Hartmann子孔径平均斜率计算假设相位屏大小为N×N布置N_sub×N_sub个子孔径每个子孔径覆盖sub_size×sub_size个网格点。子孔径的平均斜率可以用该区域相位梯度的平均值表示。用NumPy的gradient函数就能快速实现def sh_slopes(phase, N_sub, sub_size): N phase.shape[0] assert N N_sub * sub_size, 尺寸不匹配 slopes_x np.zeros((N_sub, N_sub)) slopes_y np.zeros((N_sub, N_sub)) for i in range(N_sub): for j in range(N_sub): block phase[i*sub_size:(i1)*sub_size, j*sub_size:(j1)*sub_size] dy, dx np.gradient(block) slopes_x[i, j] dx.mean() slopes_y[i, j] dy.mean() return slopes_x, slopes_y sx, sy sh_slopes(phase, N_sub8, sub_size8)这里的dx和dy分别表示相位在x和y方向上的差分单位是弧度/网格。由于子孔径内取平均局部噪声会被平滑。N_sub决定波前传感器的空间分辨率子孔径数越多能采样的高频信息越丰富但每个孔径收集的光子数减少信噪比会下降。反过来子孔径太少则斜率信息不足重建出来的波前会忽略高阶像差。实际系统里常用8×8或16×16仿真可以按需求赋值。3.2 探测器像素化与噪声模型真实探测器的光斑位置估计受读出噪声和光子噪声影响因此斜率测量值是带有误差的。仿真时可以在理想斜率上叠加高斯噪声噪声尺度用弧度表示def add_measurement_noise(slopes, sigma0.01): noise np.random.normal(0, sigma, sizeslopes.shape) return slopes noise sx_noisy add_measurement_noise(sx, sigma0.02) sy_noisy add_measurement_noise(sy, sigma0.02)这里的sigma需要结合探测器性能选择。一个经验值是取子孔径像素衍射极限的1/10比如工作波长632nm、子孔径直径2mm时光斑质心误差约0.5μm换算到斜率约0.4μrad。但如果仿真网格尺寸和波长没有严格绑定sigma直接按弧度设成0.010.05即可。注意sigma过大会让后面重建矩阵小型特征值被噪声主导导致重建结果出现棋盘格状抖动。这是AO仿真里最典型的“噪声-分辨率矛盾”后面一章会看到正则化如何缓解。4. 波前重建从斜率到相位网格传感器给出的是斜率数组而不是相位。重建器要把斜率映射回二维相位分布。这里采用Southwell网格重建它假设相位采样点位于子孔径中心x方向相邻相位点的差值等于该方向斜率乘以子孔径间距dy方向同理。把所有方程写成线性系统用最小二乘求解。4.1 Southwell最小二乘重建矩阵设子孔径网格为m×m相位向量长度为m²。对每个水平相邻点对有m×(m-1)个方程对每个垂直相邻点对有(m-1)×m个方程。再增加一个活塞零约束保证方程有唯一解。实现如下def southwell_reconstruct(sx, sy, d): m sx.shape[0] n m * m A [] b [] # x方向约束 for i in range(m): for j in range(m-1): row np.zeros(n) row[i*m j] -1 row[i*m j1] 1 A.append(row) b.append(sx[i, j] * d) # y方向约束 for i in range(m-1): for j in range(m): row np.zeros(n) row[i*m j] -1 row[(i1)*m j] 1 A.append(row) b.append(sy[i, j] * d) A np.array(A) b np.array(b) # 固定第一个相位值为0消除活塞平移自由度 A np.vstack([A, np.zeros(n)]) A[-1, 0] 1 b np.append(b, 0) phi, _, _, _ np.linalg.lstsq(A, b, rcondNone) return phi.reshape(m, m) reconstructed southwell_reconstruct(sx_noisy, sy_noisy, d1.0)代码先构造系数矩阵A每一行代表一个差分方程。x方向约束里-1和1分别对应左右相邻相位点右侧等于斜率乘以间距d。y方向类似。最后加的零约束行用于去掉全活塞项否则A不满秩最小二乘结果会带上任意常数偏移。使用lstsq而不是直接求逆是因为该线性系统是超定的最小二乘解在数值上更稳定。若子孔径数增加到16×16矩阵规模约为512×256计算瞬时完成但换成真实硬件在线运行时需要提前把A的伪逆算好避免每帧都做SVD。4.2 正则化与重建误差评估当斜率噪声sigma较大时直接最小二乘重建的波前会高频振荡。一个常见改进是给最小二乘增加Tikhonov正则项将目标变为min ||Aφ - b||² λ||Lφ||²其中L是二阶差分算子。在代码里可以简单地在反演前给A的协方差矩阵加一个小对角阵AtA A.T A Atb A.T b recon np.linalg.solve(AtA 0.01 * np.eye(n), Atb).reshape(m, m)这里的0.01就是正则化系数它等价于假设噪声方差与信号方差之比为0.01。系数过小无法压制噪声过大会把真实像差抹平。调试方法是用一组已知相位屏生成斜率并加噪声然后改变λ看重建误差RMS曲线选择曲线拐点处的λ。在工程上这一步往往比重建算法本身更影响性能因为真实系统的斜率噪声不可能为零。5. 变形镜响应与控制闭环变形镜Deformable Mirror, DM是校正阶段的执行器通过施加电压使镜面产生面形变化。在仿真里DM可以用一组影响函数的线性组合表示。每个驱动器的电压变化会在镜面上产生一个近似高斯形状的位移多个驱动器叠加得到总相位校正量。控制器根据重建出的残留相位计算电压增量形成闭环。5.1 影响函数与电压到相位的响应矩阵影响函数定义为一个驱动器在镜面上产生的相位分布通常用高斯函数模拟。把驱动器放在均匀网格上每个驱动器占据一个列向量所有列拼成响应矩阵R。生成R的代码如下def build_dm_response(N, drive_positions, width): x np.linspace(-1, 1, N) X, Y np.meshgrid(x, x) R np.zeros((N*N, len(drive_positions))) for k, (xc, yc) in enumerate(drive_positions): r2 (X - xc)**2 (Y - yc)**2 R[:, k] np.exp(-r2 / (width**2)).ravel() return R # 假设 5x5 驱动器间距0.5高斯宽度0.8 gx, gy np.meshgrid(np.linspace(-1, 1, 5), np.linspace(-1, 1, 5)) pos list(zip(gx.ravel(), gy.ravel())) R build_dm_response(N, pos, width0.8)注意影响函数宽度要设置成驱动器间距的0.61.2倍。宽度太小会导致镜面产生波浪形起伏驱动器之间无法平滑过渡宽度太大则相邻驱动器耦合格外强烈矩阵R条件数变大求逆时会把噪声放大。生成R后可以用np.linalg.pinv(R)计算电压到相位的逆映射但实际控制器中更常用的是预先计算控制矩阵把重建相位直接映射到电压增量。5.2 闭环积分控制与仿真发散经典AO闭环是纯积分控制当前控制电压 上一步电压 增益 × 校正电压。仿真每一帧先计算当前残余相位下的斜率重建出相位再乘以响应矩阵的伪逆得到电压增量最后更新残余相位。实现示意如下def closed_loop(init_phase, R, R_inv, N_sub, sub_size, gain0.5, n_iter30): residual init_phase.copy() V np.zeros(R.shape[1]) history [] for k in range(n_iter): sx, sy sh_slopes(residual, N_sub, sub_size) sx np.random.normal(0, 0.01, sx.shape) sy np.random.normal(0, 0.01, sy.shape) phase_rec southwell_reconstruct(sx, sy, d1.0) V_new gain * (R_inv phase_rec.ravel()) V V V_new correction (R V).reshape(init_phase.shape) residual init_phase - correction history.append(np.sqrt(np.mean((residual - residual.mean())**2))) return residual, history res, hist closed_loop(phase, R, np.linalg.pinv(R), 8, 8, gain0.4)这段代码模仿真实CCD采样每次循环都对斜率加噪声。gain是环路增益典型范围0.10.8。增益过大一步校正量超过当前误差会来回震荡并出现仿真发散表现为RMSE曲线先下降后迅速上升增益过小收敛速度慢在动态湍流下永远追不上畸变变化。另一个常见发散原因是phase_rec重建时未消除活塞项导致校正量中混入整体平移不过这不会影响像差RMS。若运行后发现residual的RMS在第10步后反弹优先把gain降到0.3以下再检查噪声sigma是否过大。仿真发散不是代码bug而是采样率、噪声和增益三者失配的物理后果需要在参数层面做平衡。6. 优化验证技巧与在线扩展闭环仿真跑通后下一步是让它给出可信的性能数字而不是只输出一堆二维图。这里给出三个实用技巧能让仿真结果更接近实验室测量。第一个技巧是用相位RMS和Strehl比做单帧评估。假设工作波长与相位屏长度单位一致残余波前方差σ²可以通过np.var(residual)计算Strehl近似用exp(-σ²)。RMS小于0.1弧度时系统已接近衍射极限此时再提升增益收益不大反而容易引起振荡。第二个技巧是扫描增益和正则化系数画出收敛曲线热力图。把两个参数各取10个值运行同一组湍流屏统计收敛后10帧的RMS平均值用matplotlib.imshow画成二维图能看到明显的稳定区域和发散区域。性能指标直接决定后续硬件选型比如控制器带宽要求由湍流截止频率决定而传感器帧率决定了能校正多大的r0。第三个技巧是离线相位屏库的复用。如果每次运行都重新生成随机屏不同参数之间的性能对比不公平。较好做法是保存一批固定种子的屏先用np.save存成.npy文件对比算法时从文件重复读取。这也能让复现实验变得容易审计时只需要提供文件和数据版本。最后把仿真算法从离线阶段迁移到实时系统时注意避免两处性能陷阱一是在每帧循环里重复构造Southwell矩阵应预先计算伪逆运行时只做两次矩阵乘二是scipy.sparse.linalg.lsqr比np.linalg.lstsq快很多在大子孔径数和多帧迭代下差距可以达到十倍。优化后的循环里单帧处理时间可以在普通笔记本上压缩到2毫秒以内这足以支撑数百赫兹的控制带宽模拟。本文还有配套的精品资源点击获取
