基于GPU 的961 单元自适应光学系统波前处理器设计及其动态性能

杨乐强，王建立，姚凯男，李宏壮，陈 璐，邵 蒙

（中国科学院 长春光学精密机械与物理研究所，吉林 长春 130033）

1 引 言

地基大口径光学望远镜的成像分辨率受到大气湍流引起的波前畸变影响会下降。自适应光学系统通过探测并补偿入射波前的相位畸变，能够有效克服大气湍流引起的相位畸变影响，目前已成为提高地基大口径光学望远镜成像分辨率的必要手段［1-3］。在自适应光学系统中，波前处理器负责解算入射波前的畸变信息，根据控制算法计算产生相应的驱动控制信号，控制波前校正器补偿入射波前畸变，是自适应光学系统的运算核心，其处理能力与计算延时直接决定了自适应光学系统的闭环动态性能。

随着地基望远镜口径的逐步增大，自适应光学系统波前校正单元数也不断增加。为满足4 m级光学望远镜的波前探测需求，自适应光学系统的规模需要提高至千单元量级，相应的波前处理器的计算规模随之呈平方倍增加。另一方面，波前探测器的采样频率要超过1 000 Hz，才能保证自适应光学系统具有足够的误差抑制带宽，从而抑制大气湍流引起的动态波前畸变。因此，大规模自适应光学系统对波前处理器的数据吞吐规模、计算速度提出了很高的要求［4-6］。

为了降低系统延迟，提高数据处理速度，传统的自适应光学系统大多采用基于数字信号处理器（Digital Signal Processor，DSP）阵列构成的波前处理器以及基于现场可编程门阵列（Field Programmable Gate Array，FPGA）的波前处理器［7-8］。美国星火靶场的SOR3.5 m 望远镜自适应光学系统，系统规模为941 单元，采用9 块FPGA 级联作为波前处理器，计算延时约为0.297 ms［9-10］。2012 年，Veran 针对三十米望远镜（Thirty Meter Telescope，TMT）在近红外波段观测的多层共轭自适应光学系统的计算需求，提出了采用16 块FPGA 刀片的ATCA 架构的处理方案［11］。中国科学院成都光电所的61 单元自适应光学系统，其波前处理器采用17 个多核DSP 组成，峰值计算速度可达8.5 亿次/秒，实验结果表明，该波前处理器的运算延时为330 μs，误差抑制带宽可达54 Hz［12］。长春光机所研制的基于FPGA 的千单元波前处理板卡，采用10 块板卡组成波前处理器，采用多线并行流水算法缩短波前处理延时，提高系统控制带宽［4-5］。2016 年，成都光电所报道了1 m 太阳望远镜自适应光学系统，该系统采用151 单元变形镜，波前处理器采用FPGA 与DSP 协同工作的方式，能够在系统采样频率为3 100 Hz 的条件下实现稳定闭环，误差抑制带宽可达170 Hz［13-14］。

基于DSP 及FPGA 的波前处理器，具有并行程度高、计算数据快以及集成度高等优点，但其多芯片分布式的计算架构调试困难，扩展性较差。随着商业化图形处理器（Graphics Processing Unit，GPU）与多核中央处理器（Central Processing Unit，CPU）的并行计算能力的逐步增强，自2010 年 以 来，GPU 与 多 核CPU 逐 渐 成 为大口径望远镜自适应光学系统波前处理计算架构的研究热点［15-17］。2012 年，Basden 等针对欧洲南方天文台计划建设的极大型望远镜系统中的自适应光学系统波前处理要求，采用CPU 及GPU 方案对其波前处理能力进行了分析［18］。对于规模为40×40 的自适应光学系统，采用E5645的双CPU 系统的波前处理延时约为323 μs。Basden 等还分析了基于GPU 的波前处理架构能力及其适应性，并提出了一系列优化措施。2015年，Basden 等针对规模数为80×80 的自适应光学系统，采用基于GPU 的波前处理方案，实现了550 Hz 的波前处理速度［19-20］。2018 年，该课题组采用多核多节点CPU 处理器作为波前处理器处理同样规模的自适应光学系统，整个链路延时为1.35 ms，处理延时测量结果的方差为20 μs［21-22］。

上述GPU 波前处理架构研究大多基于对计算架构的波前处理计算速度及其适应性进行分析，仅停留在分析和模拟的阶段，实际的大规模自适应光学系统的系统级性能分析结果相对较少。本文在基于GPU 的波前处理架构的基础上，通过桌面961 单元自适应光学系统的动态校正实验，验证了基于GPU 的波前处理架构的动态像差抑制能力，通过湍流屏动态模拟分析了不同格林伍德频率下自适应光学系统的成像结果，并对其动态性能进行了分析。

2 桌面961 单元自适应光学系统组成

自适应光学系统的原理示意图如图1 所示。该系统主要由波前传感器、波前处理器以及波前校正器3 个子系统构成。其中，波前传感器负责测量入射波前的畸变信息，并将测量结果发送给波前处理器；波前处理器接收波前传感器的测量结果，并计算得到需要施加到波前校正器上的控制量；波前校正器根据控制结果产生相应波前补偿入射的畸变波前，完成自适应光学系统波前控制闭环。经过自适应光学系统校正后的远场成像，可以达到近衍射极限的成像分辨率。

图1 自适应光学系统原理Fig.1 Principle diagram of adaptive optics system

图2 为桌面961 单元自适应光学系统光路示意图。白光光源发出的光束经过湍流模拟器后通过准直透镜准直为平行光束，平行光束入射到两级快速反射镜构成的倾斜像差校正系统，进行倾斜像差校正，经过快速反射镜的光束再经过扩束元件扩大后入射到961 单元变形镜上进行高阶像差的闭环校正。经过变形镜校正后的光束通过分色镜将光分为两束，其中波段为700～900 nm 的光经过透镜会聚成像到高分辨率成像相机中，而另外一束500～700 nm 的光则经过缩束系统将光束口径缩小后入射到夏克-哈特曼波前传感器进行波前探测。夏克哈特曼波前传感器与变形镜构成闭环控制系统进行系统的入射波前残差校正。系统中961 单元变形镜的有效光瞳口径为235 mm，变形镜的促动器采用PZT 型压电陶瓷促动器，单个促动器行程为5 μm，促动器间距为7 mm，呈矩形排布方式。夏克-哈特曼波前传感器的有效子孔径数为1 020，微透镜阵列排布为37×37，单个微透镜口径为200 μm，焦距为7 mm，微透镜阵列成像点斑经过中继匹配镜组缩放后入射到波前传感器相机靶面中，以满足子孔径成像艾里斑的尺寸要求；波前传感器相机采用FirstLight 公司Ocam2 相机，相机靶面分辨率为240×240，像元尺寸为24 μm，最高采样帧频可达2 000 Hz。成像相机采用滨松公司的ORCA Flash4.0 科学级成像相机，相机分辨率为2 048×2 048，像元尺寸为6.5 μm。

图2 961 单元自适应光学系统测试光路Fig.2 Optical layout of 961-element adaptive optics testbed system

3 基于GPU 的波前处理器的设计与优化

3.1 波前处理算法

自适应光学系统波前处理算法主要由波前斜率解算、波前复原运算以及波前控制量解算三部分构成，其计算流程如图3 所示。波前传感器获取的带有入射波前畸变信息的图像，并将图像传输给波前处理器，波前处理器接收波前传感器的图像首先进行波前斜率计算，夏克哈特曼波前传感器第k个子孔径波前斜率的计算公式如下：

其中：（xi，yj）是像素在子孔径的x和y方向上的坐标，Ii，j则是在（xi，yj）坐标点上的像素灰度值。（xc，yc）为实际入射波前经过夏克哈特曼波前传感器单个微透镜后聚焦在探测器靶面上的光斑实际质心位置，（xref，yref）则为入射波前为理想平面波时探测器靶面形成的光斑质心位置。由实际质心位置与理想光斑质心位置的偏差，就可以获得夏克-哈特曼波前传感器内子孔径的二维波前斜率信息（Δxi，Δyi）。将每个子孔径内的斜率向量按x和y方向重新排列成一维列向量，就可以获得波前斜率向量：

基于波前斜率向量，本文通过直接斜率法进行波前复原计算。该方法通过矩阵向量乘法的方式得到相应的复原电压，可以表示为：

上述形式的矩阵表示为：

其中：E是误差向量，是复原运算的结果；D是响应矩阵的广义逆矩阵，一般称为控制矩阵。

根据波前误差向量进行波前控制运算，本文采用自适应光学系统中常用的PI 控制运算［12，23］。

3.2 基于GPU 的波前处理架构

基于GPU 的波前处理架构如图3 所示，主要由图像采集卡、负责波前处理的GPU 以及负责协同控制与数据传输的CPU 构成。在波前处理过程中，哈特曼波前探测器作为数据源，将微透镜阵列所成的光斑阵列图像通过CameraLink 线缆发送到GPU 所在的波前处理计算机的图像采集卡中。图像采集卡在接收完一帧图像信息后，通过DMA 的方式将波前探测器的图像传输到上位机内存当中，同时以中断的形式通知CPU 波前探测器；CPU 响应该中断，将接收到的图像数据通过PCIe 总线传输到GPU 显存部分。GPU接收到图像后，首先进行波前斜率计算即子孔径质心运算，得到质心偏移向量；然后进行波前复原运算，将得到的质心偏移向量和初始化过程中就已经拷贝到显存当中的波前复原矩阵相乘，得到波前误差向量；最后根据波前误差向量，进行波前控制运算，再根据PID 控制参数计算得到变形镜驱动电压控制量。GPU 同样通过PCIe 总线传输的方式将显存中的驱动电压控制向量传回CPU 内存当中，CPU 将驱动电压控制量打包并使用UDP 通信协议通过千兆网口将它发送到数字模拟转换卡（DA），转换后的模拟控制量经过放大器放大后驱动变形镜对应的压电陶瓷促动器产生位移形变，生成共轭波前，完成自适应光学系统闭环校正的一次迭代过程。

图3 基于GPU 的自适应光学系统波前处理算法流程Fig.3 Flow chart of wavefront processor of GPU-based adaptive optics system

3.3 基于GPU 的波前处理算法优化

在GPU 并行计算框架中，CUDA 和OPENCL 是两种常见的GPU 并行计算框架。由于CUDA 计算框架广泛应用于深度学习、大规模并行计算中，本文选择CUDA 计算框架作为GPU的并行计算框架。

为了提高基于GPU 波前处理架构的波前处理速度，本文在基于CUDA 计算框架的波前处理算法中进行了一系列的算法映射与优化设计，尽可能地提高GPU 计算资源的利用效率，进而提高波前处理速度。主要的算法映射过程与优化处理步骤如下：

（1）减少无效图像带来的传输开销。在基于GPU 的波前处理机构中，图像数据在内存与显存之间的交换导致GPU 计算延时，需要尽可能降低由传输数据带来的传输延时。夏克-哈特曼波前传感器在波前斜率解算过程中，其有效数据仅为子孔径图像信息，而其余部分为无效图像数据，因此可以通过去除图像中无效像素信息的方式，降低内存与GPU 之间传输的通信数据量。根据961 单元自适应光学系统夏克-哈特曼波前传感器的设计参数，有效子孔径占据的像素数为1020×6×6=36 700，而相机全靶面像素数为240×240=57 600，通过无效像素排除，可以降低36%的数据传输量。

（2）波前斜率计算算法映射与优化。由于夏克-哈特曼子孔径之间的斜率解算具有天然的并行计算特点，同时每个子孔径仅有36 个像素参与计算，计算规模较小，因此本文将每个子孔径的斜率计算分配到GPU 的每个计算线程当中，1 020 个子孔径的斜率计算对应在并行的1 020个GPU 线程当中。由于计算规模较小，在质心计算过程中，GPU 的内存访问开销不可忽略，为了提高质心计算速度，本文将夏克-哈特曼波前传感器的图像按照子孔径排序的方式进行了重排与向量化操作，使每个并行计算的线程能够快速连续访问子孔径的图像信息，提高显存的读取效率。在得到子孔径质心向量后，需要对质心向量解算平均值，将倾斜像差数据单独发送到快速反射镜驱动器进行倾斜像差校正，减小变形镜的行程负担。GPU 计算得到的质心误差向量是中间结果，位于访问速度最快的寄存器当中，因此质心误差的平均值计算采用CUDA 架构中的__shuffle__指令，通过该指令线程间可以进行信息交互，提高内存访问效率进而提高计算速度。

（3）波前复原运算算法映射。波前复原向量通过质心误差向量与复原矩阵相乘得到，其运算实质是向量矩阵乘法，因此本文使用CUDA 架构中的cublas 库，采用列主序的数据排列方式进行矩阵向量乘法运算。

（4）波前控制算法映射。得到波前复原误差向量后，波前控制运算可以在GPU 中映射为向量的乘加运算，同样采用cublas 库进行基于PID的控制计算得到波前控制向量。为了减少CPU与GPU 之间的通信开销，波前控制向量的存储采用CUDA 架构中的零拷贝内存，CPU 可以直接访问该内存，减少了通信开销。

4 实验结果与分析

根据上述波前处理器架构与算法设计，将基于GPU 的波前处理器应用到961 单元变形镜自适应光学系统中。其中，波前处理器采用的GPU 为GeForce RTX2080Ti 显卡，CPU 采用Intel i9-9900，8 核心16 线程，主频为3.6 GHz。

4.1 系统计算时延

为了验证GPU 实现波前处理的有效性，本文还进行了基于CPU 的波前处理实验，与GPU波前处理延时结果进行对比，其中CPU 波前处理中波前复原运算中最耗时的矩阵乘法部分采用OpenCV 库函数进行矩阵乘法，统计10 000 帧波前处理延时后，得到波前处理延时均值，如表1所示。由表1 可以看出，基于CPU 的波前处理延时平均值为1.53 ms，当系统采样频率超过1 000 Hz时，无法满足波前处理器的实时处理要求；而基于GPU的波前处理延时平均值为0.3 ms，满足系统采样频率为2 000 Hz时的波前处理速度要求。

表1 波前处理延时结果Tab.1 Results of wavefront process delay

4.2 静态像差校正实验

首先对桌面自适应光学系统进行了静态像差展平实验。图4 为一组变形镜静态像差校正实验校正前后夏克-哈特曼波前传感器测量得到的波前像差结果，由图可知，校正前系统像差峰峰值（Peak Valley，PV）为5.30λ（λ=600 nm），均方根值（Root Mean Square，RMS）为0.81λ；校正后，波 前 残 差PV 值 为0.24λ，RMS 值 为0.05λ。校正前后远场光斑图像如图5 所示，可以看出，校正前远场光斑图像的灰度峰值为3 347，校正后的灰度峰值为45 876，远场光斑的能量集中度明显提高。

图4 自适应光学系统校正前后系统波前像差Fig.4 System wavefront errors before and after optical adaptive corrections

图5 校正前后光源远场光斑图像Fig.5 Far-field images of light source before and after correction

4.3 动态像差校正实验

动态像差校正实验通过湍流模拟器模拟大气扰动。湍流模拟器通过转动刻蚀在位相屏表面特定的光程差部分形成动态像差，动态像差能够模拟大气湍流的变化情况。其中，湍流模拟器所模拟的大气格林伍德频率与风速之间的关系为：

其中：r0代表大气相干长度，vw代表风速，通过改变湍流模拟器转速即可改变模拟的大气格林伍德频率。实验中，湍流模拟器等效模拟的大气相干长度为11 cm，湍流屏转速为1 rad/s，等效模拟的格林伍德频率为60 Hz。

为验证变形镜的动态性能，实验中采用变形镜进行像差校正，快速反射镜不进行倾斜校正，仅作为普通反射镜使用（实际应用中通常采用快速反射镜校正行程较大的倾斜像差）。由湍流模拟器产生的倾斜像差全部由变形镜进行校正，根据校正前后倾斜像差与波前RMS 值的动态变化，衡量波前处理系统的动态性能。首先，测试动态条件下自适应光学系统对倾斜像差的校正结果，如图6 和图7 所示。可以看出，自适应光学系统在校正前，由湍流屏引入的系统X轴倾斜量均方根值为0.18″，Y轴倾斜量均方根值为0.26″；经过变形镜校正后，X轴倾斜量均方根值下降到0.02″，Y轴倾斜量均方根值下降为0.03″。校正前后系统入射波前（去除倾斜像差）的均方根值变化如图8 所示，校正前系统入射波前为0.5λ～1.3λ（λ=600 nm），均值为0.83λ；校正后系统波前残差明显降低，波前残差均值为0.16λ，约为100 nm。

图6 校正前后X 轴倾斜误差曲线Fig.6 Track errors of X-axis before and after correction

图7 校正前后Y 轴倾斜误差曲线Fig.7 Track errors of Y-axis before and after correction

图8 校正前后系统波前残差曲线Fig.8 Residual wavefront error before and after correction

图9～图11 分别为X轴倾斜、Y轴倾斜以及高阶像差均方根值校正前后的功率谱变化情况。由图可以看出，系统的0 dB 误差抑制带宽为100 Hz，高频部分存在校正后功率谱高于校正前的情况，这是因为在工程实际过程中，为了保证系统具有较高的动态特性，控制参数选取较大，相位裕度较小，存在放大高频噪声的情况［12，23］。

图9 校正前后X 轴倾斜误差功率谱Fig.9 Power spectra density of X-tilt before and after correction

图11 校正前后系统波前残差功率谱Fig.11 Power spectra density of residual wavefront error before and after correction

图10 校正前后Y 轴倾斜误差功率谱Fig.10 Power spectra density of Y-tilt before and after correction

图12 是不同格林伍德频率下，961 单元自适应光学系统校正前后成像相机远场能量的三维分布。可以看出，当自适应光学系统闭环时，远场能量峰值得到很大提高，但随着格林伍德频率的提高，自适应光学系统校正后的效果逐渐变差，其原因是波前处理器的闭环控制带宽有限，随着格林伍德频率的提高，时域误差也随之提高，并逐渐成为波前残余误差的主要成分，从而降低自适应光学系统的闭环控制性能。

图12 不同格林伍德频率下校正后远场光斑能量分布Fig.12 Far-field images of light source before and after correction at different Greenwood frequencies

5 结 论

本文根据961 单元自适应光学系统波前处理器的性能要求，采用NVIDIA 公司GeForce RTX2080Ti 显卡作为波前处理运算核心，实现了961 单元自适应光学系统的闭环波前处理运算。实验结果表明，对于961 单元自适应光学系统，基于GPU 的波前处理器在系统采样频率为1 500 Hz 时，系统误差抑制带宽可达100 Hz，满足千单元级自适应光学系统的实时波前处理需求。