引　言

光学平面是光学系统和仪器工程中最常见的光学元件。在光学干涉测量中，平面波前像差标准是其他各种测试的基础。标准参考面的形状会影响传统相对测量方法的精度。当需要较高的测量精度（如精度高于1/20波长）时，必须消除干涉仪的系统误差和参考面误差。为了解决这个问题，提出了绝对平面测量。绝对平面测量是通过消除光学仪器的系统误差和参考平面误差对测量结果的影响，来获得测试平面的绝对表面的方法。目前，使用最广泛的方法是由Fritz^[1]于 1984年提出的，是一种基于传统三面互检的可编程三面互检法。首先，利用 Zernike 多项式性质将平面的表面形状误差分解为正交基函数；然后，再利用最小二乘法对这些基函数的系数进行拟合；最终，得到被测绝对平面。许多研究人员在此方法的基础上进行了改进研究^[2]，并提出了许多新方法^[3-5]。但由于计算复杂或条件苛刻，其中大多数方法并未在工业领域得到广泛应用^[6]。随着测量技术的发展，人们对大型精密光学测量仪器的需求更加迫切，如美国国家点火装置以及中国神光–Ⅲ装置^[7]均需要数千块超大口径光学器件。这些超大口径光学器件对干涉仪器的横向分辨率提出了更高的要求，然而现有算法的测量效率尚不能满足这一要求。因此，本研究采用了图形处理单元（graphics processing unit，GPU）加速技术，旨在提升测量效率。在干涉测量仪器中获取到干涉图像后，先提取相位，然后对提取的相位进行解包裹，最后使用Zernike多项式拟合的方法实现绝对平面测量。在此流程中，针对干涉图的异构并行相位提取和相位解包裹已有诸多学者进行了探索和研究^[8-9]。但是，对Zernike多项式拟合的绝对平面测量步骤进行探索的还比较少。在Zernike多项式的应用方面，大多数研究人员关注的是 Zernike 多项式的 GPU 并行实现及其在图像处理中的应用^[10]。例如，Toharia等^[11]将 Zernike多项式应用于射击边界检测，并通过多GPU 多中央处理单元（central processing unit，CPU） 异构计算实现了超过1 000帧/s的检测率。然而，针对光学波前像差拟合异构并行技术，测量参数求解的并行化和分辨率优化的研究较少。

本文旨在解决传统绝对平面检测Zernike多项式计算算法在处理高分辨率干涉图时效率低下的问题。为了提高计算效率，提出了一种CPU– GPU异构并行加速算法。该算法通过在CPU和GPU之间分配计算任务，优化了笛卡尔坐标系下的Zernike多项式计算、多项式拟合系数求解、理想平面拟合以及测量指标[如峰谷（peak-to-valley，PV）值和均方根（root mean square，RMS）值]计算等关键步骤。此外，本文还采用协作组线程束原语指令和混合精度方法，进一步提升了算法的执行效率。

1 并行算法设计 1.1 异构并行计算

传统的科学计算性能提升依赖于 CPU 的时钟频率。然而，芯片上高密度晶体管的功耗和散热所造成的物理瓶颈极大地延缓了 CPU 频率的提升^[12]。与 CPU 不同，GPU 自诞生之日起就被用于计算与图像相关的数据。

从图1可知，GPU 包含比 CPU 更多的核心数量，因此，GPU比CPU更适合处理大规模数据的并行计算，但它的核心频率、缓存和控制单元较弱，对复杂问题处理能力较弱。为了充分发挥两者优势，有人提出了异构并行计算方式。异构并行计算就是将算法分解为复杂流程控制部分和简单并行处理部分，由CPU处理复杂逻辑运算部分，GPU处理简单高吞吐运算部分，从而达到高性能的诉求。在科学计算领域，越来越多的学者采用异构并行计算来解决相关领域问题：赵锴坤等^[13]将CPU/GPU异构并行计算用于物理领域，计算加速重力场反演，将计算性能提升了380倍；李浩^[14]将异构并行计算用于天文数据处理，在进行相干消色散脉冲星数据处理时，该计算相较于CPU处理速度提升了2倍。

先讨论有关GPU并行算法的设计和应用，图2 展示了GPU内核运行方式。在CPU端启动核函数并根据问题动态分配资源后，GPU会根据核函数中的代码进行运行。在GPU端运行核函数时，每个核函数会启动线程网格，每个线程网格由若干线程块组成，而线程块由多个线程构成，这些线程负责处理数据。同一个线程块中若干个线程组成一个线程束，线程束中的线程以不同数据资源执行相同的指令。

利用CPU/GPU异构技术来实现基于Zernike多项式绝对平面测量算法。在实际异构加速计算过程中，GPU部分主要负责对数据进行运算，CPU部分负责逻辑流程处理以及CPU/GPU之间的数据通信，整体的流程如图3所示。在本研究中，CPU主要负责初始化需要拟合的Zernike项数，加载面型数据，判断当前项数是否满足测量需要，根据待测平面分辨率分配GPU资源（包括GPU内存和线程数量以及线程块数量），启动GPU计算核心以及从GPU获取计算结果。GPU主要负责计算Zernike多项式在笛卡尔坐标系下的映射，拟合Zernike多项式参数，拟合理想面型以及计算PV值和RMS值。其中，整个计算最复杂的部分就是GPU算法的设计，接着将对这一部分的实现方法进行介绍。

1.2 Zernike多项式映射计算

Zernike 多项式是一组定义在单位圆域内的完全正交基^[15]。其数学描述包括2个主要部分

$ {R}_{n}^{m}\left(\rho \right) = \frac{{\sum }_{x=0}^{\left(n - \left|x\right|\right)/2}{\left(- 1\right)}^{s}\left(n - s\right)!}{s!\left[0.5\left(n+\left|m\right|\right) - s\right]!\left[0.5\left(n - \left|m\right|\right)- s\right]!}{\rho }^{n-2s} $

(1)

$ \Theta \left(m\theta \right)\left\{\begin{array}{c}{\mathrm{cos}}\left(m\theta \right),m\geqslant 0\\ {\mathrm{sin}}\left(-m\theta \right),m < 0\end{array}\right. $

(2)

式（1）为独立的径向多项式$ {R}_{n}^{m}\left(\mathrm{\rho }\right) $，式（2）为角多项式$ \mathrm{\Theta }\left(m\mathrm{\theta }\right) $。$ {R}_{n}^{m}\left(\mathrm{\rho }\right) $的各项只与径向有关，n 表示它们的多项式阶数，s 表示多项式次数，ρ 表示弧度，m 表示角频率。$ \mathrm{\Theta }\left(m\mathrm{\theta }\right) $的各项与振幅角度有关，θ表示振幅角度。其极坐标表达式为

$ {Z}_{n}^{m}\left(\rho ,\theta \right)={R}_{n}^{m}\left(\rho \right)\cdot \Theta \left(m\theta \right) $

(3)

Zernike 多项式与光学设计中的各种经典像差相对应，常用于描述光学系统的波前像差。

图4为Zernike多项式映射的异构并行的计算流程在5 × 5分辨率下的示例图。在干涉测量中，通常是对解相位包裹后的图像进行测量，解相位包裹后的数据使用直角坐标表示。因此，需要建立径向多项式$ {R}_{n}^{m}\left(\rho \right) $和角多项式$ \Theta \left(m\theta \right) $单位极坐标与直角坐标的映射关系，即图4中cart2pol计算部分。这种关系的转换可依据待测面型的宽度和高度信息生成。假设待测拟合波面$ W\left(\rho ,\theta \right) $的像素尺寸为ROWs×COLs，在区间 [−1，1] 内生成网格采样点，然后根据式（4）和式（5） 的转换，便可得出这组映射关系。

$ Radial\_map\left(x,y\right)=\sqrt{X{\left(x,y\right)}^{2}+Y{\left(x,y\right)}^{2}} $

(4)

$ Angular\_map\left(x,y\right)={\mathrm{arctan}}\left[\frac{Y\left(x,y\right)}{X\left(x,y\right)}\right] $

(5)

可以发现，该计算中像素之间数据不存在依赖关系。因此，可以借助GPU多核特性来分配像素尺寸的线程数量并进行计算，快速得到该分辨率直角坐标系下的极坐标信息，然后按照式（1）和式（2）便可分别计算$ {R}_{n}^{m}\left(\rho \right) $和$ \Theta \left(m\theta \right) $。最后，通过对应相乘（图4中mult计算部分）可得到当前分辨率下笛卡尔坐标系的Zernike多项式。

1.3 Zernike多项式系数求解

由于 Zernike 多项式与光学像差相对应，因此任何波前像差都可以用 Zernike 多项式作为基函数表示为

$ W\left(x,y\right)={\sum }_{i=1}^{n}{a}_{i}{Z}_{i}\left(x,y\right) $

(6)

式中：$ {Z}_{i}\left(x,y\right) $ 是第i项Zernike多项式；$ {a}_{i} $ 是该多项式对应的 Zernike 系数；$ W\left(x,y\right) $是真实离散面型数据。为了方便计算机运算，将式（6）转为矩阵表示形式

$ {\boldsymbol{W}}={\boldsymbol{Z}}{\boldsymbol{A }}$

(7)

式中：$ {\boldsymbol{A}}={\left({a}_{1},{a}_{2},\cdots ,{a}_{n}\right)}^{{\mathrm{T}}} $；$ {\boldsymbol{Z}}={\left({Z}_{i}\right)}_{R\times C} $，R×C为分辨率。多项式拟合是为了求解其系数，即公式中的系数向量A。如式（7）所示，需要根据已知真实波面数据$ W\left(x,y\right) $计算出相应的系数A。借助GPU可实现LU分解算法，通过三角求解器可求解该问题。

如表1所示，当前异构传统直接算法性能相较于单纯的CPU计算性能提升了4～8倍。考虑到数值算法的发展以及GPU 低精度矩阵乘硬件的性能优势，本文将使用混合精度迭代法作为求解器，进一步提升线性求解器性能^[16]。

如图5所示，混合精度迭代法通过混合计算半精度（FP16）和单浮点精度（FP32）在传统的直接异构并行方法基础上将性能提升了1～5倍。同时可以发现，随着待测面型的分辨率增大，线性求解器获得的性能加速比就越大，在4 096 × 4 096分辨率条件下其性能提升了约5倍。

计算出所有系数后，根据需要用 Zernike 多项式与系数计算拟合波面。在实际计算过程中，借助GPU核心数量多，理论吞吐量大的优势为每个独立像素分配线程来完成计算。在假设求解系数为36项的基础上测得的性能数据如图6所示，可得该方法在超高分辨率（4 096 × 4 096）情况下相较于传统方法的性能提升了71.8倍。

1.4 PV值与RMS值异构并行计算

在绝对平面测量中，通常使用峰谷值（V_pv）和均方根值（V_rms）作为评估指标

$ V_{{\mathrm{pv}}}={\mathrm{max}}\left[H\left(x,y\right)\right]-{\mathrm{min}}\left[H\left(x,y\right)\right] $

(8)

$ V_{{\mathrm{rms}}}=\sqrt{\frac{{\sum }_{x,y}H{\left(x,y\right)}^{2}}{RC}} $

(9)

式中：$ H\left(x,y\right) $为拟合波面与真实波面$ W\left(x,y\right) $的差；R为分辨率宽度；C为分辨率高度。

根据式（8）可知，V_pv是通过对矩阵中的最小和最大值做差求得。该计算数据之间的依赖性较强，无法使用Zernike多项式映射的方式加速，可通过组合规约并行模型来进一步提高求解速度。在Chen等 ^[8]和Jradi等 ^[17]提出的CPU/GPU异构规约模型中，存在较多分支语句，GPU 逻辑处理分支的效率较差，复杂度较高。其具体解决方案的步骤如图7（a）所示，计算步骤需要N/2步，其中N为线程块中线程数量。当前线程块线程束（线程束由多个线程组成）的数量为2。为进一步提升效率，对其进行了性能优化，主要借助协作组线程束规约原语指令，该指令通常用于GPU程序性能优化，能够将每个线程束中的数据进行单指令规约^[18]。图7（b）所示为使用指令优化后的峰/谷值求解算法，优化前需要N个线程迭代N/2次，每次对问题规模进行折半计算，直到求解出最后的元素。优化后，只需要执行2步即可完成计算：首先，通过规约原语指令将每个线程束的最小/最大值筛出；随后，只需比较线程束中筛出元素即可。算法的复杂度和耗时均降低。

图8所示为当前实验环境下求解峰值或者谷值的性能对比图。本文算法相较于未优化的传统异构并行计算峰/谷值的算法快8～20倍，并且更适合处理大分辨率的情况。求峰谷值需要计算峰值和谷值，并对它们进行做差处理，因此需要在上述算法进行2次的基础上增加做差核函数。

根据式（9）可知，求解RMS值的计算中，需要计算波前拟合误差矩阵所有元素和的平方，该算法与峰谷值求解逻辑相似。此外，需要对求得的和进行代数计算，可利用额外的核函数来进行处理。

2 实验与分析 2.1 实验平台与仪器

基于上述并行异构算法搭建实验平台，CPU 为 AMD Ryzen 5900H，加速 GPU 为 RTX3070– Laptop，两者的具体参数见表2，整机的本地内存为 16 GB。GPU 部分的代码根据 NVCC–11.6 来编译，CPU 部分代码用 GCC–7.3.0 来编译，测试环境操作系统内核为Linux–5.17。

所使用的数据为图9（a）所示的平面面型数据，该数据由干涉条纹解相位包裹获得。干涉条纹数据通常由图9（b）所示干涉仪捕获。本文测试数据由苏州慧利仪器有限责任公司提供，共4组数据，每组数据具有不同的分辨率。

2.2 功能验证

在性能验证之前，先对算法进行功能验证。将图9（a）所示样例作为算法的输入，可以得到如图10所示的面型。该面型与常规计算方法算得的面型相同。

2.3 性能验证

在当前实验环境下，分别对Zernike多项式计算、Zernike 系数求解以及PV 值和 RMS 值求解等步骤进行性能测试。

图11所示为不同步骤在4组不同分辨率待测面型情况下异构计算加速比。其中图11（a）为笛卡尔坐标系下的多项式计算，加速比为50～70倍；图11（b）为多项式拟合，加速比为6～25倍；图11（c）为测量面型、PV值和RMS值计算部分，加速比为50～70倍。与常规计算相比，异构计算所获得的性能增益随着像素的增加而增加。

图12所示为在4组不同分辨率待测面型情况下不同步骤的异构计算标准差。其中图12（a）为笛卡尔坐标系下的多项式计算，标准差为0.83～5.28 ms；图12（b）为多项式拟合，标准差为0.90～1.29 ms；图12（c）为测量面型、PV值和RMS值计算部分，标准差为0.08～0.16 ms。3个步骤的标准差均随着分辨率增大而减小，原因在于随着数据量的增加，GPU的调度器可以更有效地分配任务，使每个处理单元的工作负载更加均衡，内存带宽的利用也更充分，从而有助于减小完成上述步骤的耗时差异，降低标准差。

表3所示为不同分辨率下采用常规计算方法和异构计算方法的平均耗时情况。通过分析可知，如果采用传统计算方法，在针对4 096 × 4 096高分辨率测量时的平均耗时约42 s，而采用异构方法平均仅需要602 ms；在针对1 024 × 1 024分辨率时，异构计算法的测量速度能够达20次/s左右，相较于常规计算法其整体加速了40～70倍。在实际的干涉测量研究和应用中，需要对面型进行重复测量，使用异构计算法将极大提高测试效率。

3 结　论

本文实现了CPU/GPU异构并行加速绝对平面检测，包括Zernike 多项式映射计算求解、Zernike 多项式拟合系数求解、绝对平面测量参数并行化算法设计与实现以及CPU流程控制部分。该方法可以提高高分辨率绝对平面测量的效率。测试了多组不同分辨率的光学平面，随着分辨率的提高，GPU 并行加速带来的性能提升也越高。在当前实验环境中，针对不同分辨率，其性能提升了40～70倍，且分辨率越大性能提升越明显，其中在处理4 096×4 096 分辨率的光学平面时，平均只需要 0.6 s左右，比常规CPU计算性能提高了70 倍左右，极大提升了干涉测量效率。未来将进一步探索异构并行在其他光学面型测量领域的应用。