一种非稳腔激光器光瞳位置提取方法

马瑞浩 ，甘永东 ，李新阳 ，梅月 ，斯那卓玛 ，朱枫 ，贾启旺 ，吴天祎

（1 中国科学院自适应光学重点实验室，成都 610209）

（2 中国科学院光电技术研究所，成都 610209）

0 引言

由于共焦非稳腔具有大的可控模体积、横模鉴别力高、近衍射极限平面波输出，可得到高功率和高质量的远场光斑等特点［1］，在化学激光系统中得到了广泛的应用。导引光在非稳腔的两个镜面间来回反射，相互干涉，理论上从非稳腔输出的光束为同心圆环［2］，但实际上，由于非稳腔的准直、空气扰动、平台震动等因素的影响，输出的干涉同心圆环灰度起伏不定，并且处于快速的动态变化之中，增大了导引光光瞳位置提取的难度。

在非稳腔的自动调腔中，也需要提取光斑的位置。2007年，熊木地、贾思楠和张增宝等采用位置传感器检测光斑位置［3］。2010年，胡永军［4］提出使用随机Hough变换椭圆拟合获取同心圆的位置，但是随机Hough变换椭圆拟合方法时间达到分钟级，不具备实时性。2020年，李景润、熊木地等［2］提出使用霍夫圆检测的方法获取同心圆的位置。

非稳腔激光器准直之后，需要将非稳腔出射的导引光与后续光路对准，由于工作环境温度变化，工作环境震动、工作平台发生相对移动等因素的影响，激光器发射出的光束不能准确的到达指定位置。手动调整光路对准需要耗费大量人力和时间资源，并且手动光路调整精度受主观因素影响，重复精度较低。因此使用光路自动准直技术，通过提取导引光光瞳图像的位置和光轴图像的位置，计算提取的位置与基准位置的偏差，通过计算出的偏差量控制伺服电机运动，将激光器发射的光束调整在光路可允许的误差范围内。

1995年，陈庆浩、徐仁芳、彭增云等［5］使用硅四象限管光电探测器提取光路自动对准中的光斑位置。2005年，上海光机所的吕凤年、刘代忠等［6］使用中值滤波、圆拟合的方法完成光路对准中的光斑位置提取。2015年，徐瑞华、何俊华等［7］使用重心法提取近场光斑的位置，完成光路自动准直。2017年李红、朱健强等［8］使用阈值分割、边缘检测和圆拟合的方式提取光路自动准直中近远场图像的位置。2022年，曾沛颖、朱宝强等［9］提出使用蚁群算法提取图像轮廓，增强自准直系统的自适应性，但是该算法的计算时间受蚂蚁爬行步数NC的影响，需要进一步优化。已有的光瞳位置提取方法主要分为：重心法、圆拟合方法、随机Hough椭圆拟合方法。由于导引光经过非稳腔后图像灰度起伏不定，并且光瞳图像处于快速的动态变换之中，重心法无法准确的提取出光瞳的位置，但是对于灰度分布均匀的光瞳图像，可以使用重心法快速提取光瞳位置。因为干涉同心圆并非严格垂直的入射到传感器上，干涉同心圆环经常以干涉同心椭圆的形式在传感器上成像，圆拟合的方法具有一定的局限性。为了更好的提取激光光束的近场图像（光瞳图像）的位置，采用椭圆拟合的方式提取光瞳图像的位置。随机Hough椭圆拟合主要存在实时性的问题，但是当需要检测多个椭圆时，可以使用随机Hough椭圆拟合方法。为减少噪声干扰，准确提取出光瞳位置，本文提出一种随机采样一致性（Random Sample Consensus， RANSAC）椭圆拟合算法的光瞳位置提取方法，可以在保证光路自动准直实时性的基础上，在一定程度上抑制轮廓的噪声对椭圆拟合精度的影响，提高了光瞳位置提取的精度和鲁棒性。

1 非稳腔光瞳位置提取方法

当非稳腔调腔完毕后，从非稳腔输出明暗交替、由内到外逐渐减弱的对称干涉圆环［4］。虽然光瞳图像的内圆环较为清晰，但内环轮廓依然具有噪声，简单的椭圆拟合的精度易受轮廓噪声的影响，为了减少边缘噪声对椭圆拟合的影响，加入随机采样一致性算法（Random Sample Consensus， RANSAC），提高椭圆拟合的精度。

1.1 光瞳位置提取方法

光瞳位置提取的方法如图1所示。输入图像经过图像预处理和阈值分割后，提取光瞳图像的内轮廓，并利用内轮廓的特征将内轮廓从众多轮廓中提取出来，利用椭圆拟合算法拟合内轮廓，使用拟合出的椭圆计算出轮廓的中心，使用轮廓的中心作为光瞳的位置。

图1 光瞳位置提取算法流程Fig.1 Flow chart of Pupil position extraction algorithm

1.2 图像预处理和阈值分割

图像预处理中采用双边滤波［10］的方式去除图像中的噪声，采用最大类间方差法（OTSU）［11］计算图像阈值，将光瞳图像和背景分开，提取出光瞳图像，然后对二值图像进行中值滤波，去除较小的噪声点并填充图像中的孔洞。

双边滤波在去除噪声的同时可以较好的保留图像的边缘信息。双边滤波器公式如式（1）所示。

式中，w(i，j，k，l)由定义域核和值域核的乘积决定，滤波窗口大小为N，i和j表示待滤波点的坐标，k和l表示滤波窗口中心点的坐标，定义域核如式（2）所示，值域核如式（3）所示，w(i，j，k，l)如式（4）所示。

原始图像和滤波后的图像如图2（a）和图2（b）所示。

图2 输入图像和双边滤波后的图像Fig.2 Input image and bilateral filtered image

使用大津算法计算图像阈值T1，对于光瞳图像，前景（即目标）和背景的分割阈值记作T1，属于前景的像素点数占整幅图像的比例记作w0，其平均灰度μ0，属于背景的像素点数占整幅图像的比例记作w1，其平均灰度μ1，类间方差记作g。采用遍历的方法得到使类间方差最大的值T1，就是大津算法所求出的阈值T1。最终图像的阈值T=T1·β。阈值分割后的图像如图3（a）所示。

图3 二值化图像和中值滤波后的二值化图像Fig.3 Binary image and binary image after median filtered

对二值化后的图像进行中值滤波，中值滤波不仅能够有效的滤除图像中的随机噪声和脉冲干扰，而且还能较好的保留图像的边缘信息［13］。滤除噪声同时填充较小的孔洞，滤波后的图像如图3（b）所示。

1.3 图像轮廓提取和轮廓识别

为提取光瞳二值化图像的轮廓，两次扫描光瞳二值化图像，第一次扫描光瞳二值化图像，寻找出图像的轮廓；然后扫描轮廓图像，对每条轮廓进行编号。光瞳二值图像轮廓提取步骤为

假设输入图像为F，建立图像f，用于存放轮廓图像，Fij表示图像灰度值，图像大小为M×N。建立变量NBD，用于标识轮廓；建立一个一维向量V0用于检测轮廓；建立一个二维向量V2，用于存放所有轮廓；建立一个一维向量V1，用于存放当前轮廓。

第一遍从上到下，从左到右扫描图像F。当扫描到某个像素点(i，j)的灰度值Fij=1时，判断(i，j)点的四邻域或八邻域内是否存在0像素点，若存在0像素点，则认为该点位边缘点，将该点灰度值置为1，即fij=1，否则将(i，j)点坐标的灰度值置为0，即fij=0。扫描完图像F后，图像f为光瞳图像的轮廓图像。

从上到下，从左到右扫描图像f，将图像f的轮廓提取出来，步骤为：

1）初始化 NBD=1。

2）从上至下，从左至右扫描图像f，若扫描到像素点fij=1时，执行步骤3），若图像扫描完毕，执行步骤 6）。

3）当扫描到fij=1时，NBD=NBD+1，fij=NBD，将(i，j)点压入向量V1中。

4）判断(i，j)点八邻域内是否存在灰度值为1的像素点，若存在灰度值为1的像素点(i1，j1)，将(i1，j1)压入V0向量和V1向量中。

5）搜索完(i，j)点的八邻域后，判断V0向量是否为空，若V0向量为空，将向量V1压入二维向量V2中，跳转至步骤2），若V0向量不为空，将向量中的最外层的坐标(i2，j2)弹出，使fi2j2=NBD。i=i2，j=j2，重复步骤 4）。

6）图像扫描结束。

光瞳图像的轮廓如图4所示，图像中不同的轮廓用不同的颜色标识出来。轮廓提取完毕后，得到光瞳图像的轮廓图。根据光瞳图像的内轮廓的中心与组成内轮廓的点数，将内轮廓从众多轮廓中提取出来。首先，光瞳的同心椭圆图像经过图像处理后，轮廓图像存在内轮廓、最外围轮廓和干扰轮廓。经过实验，发现干扰轮廓形心距离光瞳位置较远，而最外围和最内围轮廓形心距离光瞳位置较近，根据这个发现，利用各个轮廓中心与整个轮廓图像的中心的距离作为判断条件，从众多轮廓中找到内轮廓和最外围轮廓。中心距离计算公式如公式6所示。经过试验，构成内围轮廓的点数范围在［380，431］，考虑冗余性，将范围扩大一倍，轮廓点数范围在［354，458］之间的，而且轮廓点数最少的轮廓则认为是内围轮廓。

图4 轮廓图像和内轮廓图像Fig.4 Contour image and inner contour image

式中，(ximg，yimg)为轮廓图像的中心，(xNBD，yNBD)为编号为NBD的轮廓的中心。

1.4 椭圆拟合

激光光束一般以非垂直入射的方式入射到传感器靶面，在传感器上成像一般为椭圆，找到内轮廓后，对内轮廓进行椭圆拟合，椭圆拟合方法参考文献［14］。

椭圆的一般方程如式（7）和式（8）所示。

当式（7）中a≠0时，αa代表了相同的圆锥曲线，在一个合适的尺度下，式（8）可以转换成式（9）。

定义A=[a，b，c，d，e，f]T，X=[x2，xy，y2，x，y，1]，则式（7）表示为式（10）。

定义一个N×6的矩阵D和一个5×6的矩阵C，

椭圆拟合问题重新表述为：在约束条件9下，求取式（10），用矩阵的方式表述为

构造拉格朗日函数，得到式（12），对式（12）求偏导数，得到式（13），把椭圆问题转换为特征向量求取问题，如式（14）所示。

为简化计算，将矩阵进行拆分，将矩阵D拆分为D=(D1]D2)，；将矩阵S拆分为。类似的将C矩阵拆分为；将A拆分为。矩阵拆分过后，式（13）转换为式（15）。

经过矩阵拆分，将6阶矩阵运算降阶为3阶矩阵运算，便于显卡进行简单的矩阵求逆和计算特征向量。

为了减少轮廓噪声对椭圆拟合结果的影响，椭圆拟合中加入RANSAC算法，RANSAC基本流程［15］为：

1）在数据集中随机选取n个数据，计算出模型的参数。

2）将所有的点作用于得到的模型，误差在设定阈值TError以内的点为内点，反之，则为外点。

3）重复步骤1）和步骤2），当迭代次数达到设定阈值后，保存内点数目最多的模型参数为最终模型的参数。

假设数据集中内点的比例为η，RANSAC算法迭代K次，Pn表示K次选择中至少一次全是内点的概率，当已知η和Pn，通过式（16）和式（17），可推导出的迭代次数K。

优化RANSAC算法的第一个步骤，在本方法中，将数据集均匀分为n份，从每份中随机采样一个数据，一共采样n个数据。

由于图像中提取的轮廓点均为整数，并不会完全准确的落在椭圆的轮廓上，采用n=6个点拟合出来的椭圆并不能非常精确的拟合出椭圆的参数。因此，将所有判断为内点的轮廓点重新进行椭圆拟合，提高椭圆拟合的精度。

拟合出椭圆后，根据式（18）计算出椭圆的中心，该椭圆的中心即为光瞳的位置，光瞳位置的计算结果如图5所示，白线为提取的光瞳轮廓，红线为拟合的椭圆轮廓，红色的十字为提取的光瞳的位置。

图5 光瞳位置计算结果图像Fig.5 The result image of pupil position calculation

2 算法测试结果分析

算法的准确性、稳定性和实时性是判断算法优劣的重要指标。使用MATLAB模拟生成光瞳图像，验证算法的准确性。通过对模拟的图像轮廓添加噪声和处理实时的光瞳图像，来验证算法的稳定性。

2.1 仿真图像测试结果分析

使用MATLAB模拟生成100帧光瞳图像，图像大小为580×580，每帧图像中心一致，中心坐标为(243.0，289.0)，旋 转 角 度 范 围 为 0°～356.4°，内 椭 圆 短 轴 范 围 为 (65.0，75.0)，外 椭 圆 短 轴 范 围 为(125.0，135.0)。为了使仿真图像更加接近真实图像，在仿真图像的基础上，增加3个处理措施：1）设置合适的图像灰度值；2）增加比例为0.2的椒盐噪声；3）对仿真图像进行高斯滤波，处理后的仿真图像更加接近真实图像。仿真的100帧光瞳图像如图6所示。

图6 第一张光瞳仿真图像和最后一张光瞳仿真图像Fig.6 The first pupil simulation image and the last pupil simulation image

实验确定Td=75时，可以将内轮廓、最外层轮廓和其他轮廓进行区分。然后确定RANSAC算法中TError取值。根据提出的算法提取出光瞳的内轮廓，提取出轮廓后，确定Pn=0.999 999 999，使n=6，TError=10～100，遍历出RANSAC参数TError的范围。使Pn尽可能的大，是为了增加RANSAC计算结果准确性的概率，使n=6是因为拟合椭圆方程中有6个变量。为了验证所提方法对噪声的鲁棒性，对图像的内轮廓添加概率为40%，幅值为0～4的噪声，计算每幅图像的位置偏移量，以光瞳位置偏移量的RMS（Root Mean Square， RMS）作为评判计算结果的判据，遍历结果如图7所示。当TError=72时，RMS值最小。因此确定RANSAC参数中Pn=0.999 999 999，n=6和TError=72。

图7 不同TError取值下的光瞳位置计算结果的RMS值Fig.7 The RMS of pupil position calculation results in different TError

为检测所提方法对噪声的鲁棒性，对提取的仿真图像的内轮廓添加0%～40%的随机噪声，随机噪声幅值分别为0，0～4，0～8，和0～13。噪声幅值为0，表示未添加噪声。以位置偏移量的RMS值作为判据，图8中蓝色曲面为圆拟合的位置计算结果，绿色曲面为直接椭圆拟合的位置计算结果，红色曲面为所提方法的位置计算结果，从图8中可以看出，在噪声较小时，直接椭圆拟合的结果优于圆拟合的结果。随着噪声比例和幅值的增加，圆拟合和直接椭圆拟合的位置计算结果RMS值增大，而所提方法的RMS值基本不变，且RMS值小于0.05。因此所提方法对噪声具有较好的鲁棒性。

图8 圆拟合、直接椭圆和所提方法结果的RMS值Fig.8 The RMS of circle fitting， direct ellipse fitting and proposed the method

图9表示圆拟合、直接椭圆拟合和所提方法在仿真图像的轮廓中添加0～40%幅值为0～13的噪声点的情况下光瞳位置的提取结果。从图9可以看出，圆拟合和直接椭圆拟合的位置偏移量均值随着噪声比例而增加，直接椭圆拟合的位置偏移量小于圆拟合的位置偏移量。但是在噪声的比例为20%时，圆拟合的位置偏移量略小于椭圆拟合。所提方法的位置偏移量基本不受噪声的影响。所提方法的位置偏移量的均值小于0.05个像素。随着噪声比例的增加，圆拟合法和直接椭圆拟合法位置偏移量的平均值（Average， AVG）与真实值之间的误差增加，偏移量最大值（MAX）和RMS值增加。而所提方法的光瞳位置偏移量结果的AVG值、MAX值和RMS值基本不变，在100帧仿真图像的轮廓中添加40%幅值为0～13的噪声点的情况下，所提方法比直接椭圆拟合法光瞳位置偏移量结果的AVG值的误差小0.283个像素、MAX值小0.660个像素、RMS值小0.136个像素。因此所提方法可以有效抑制噪声对光瞳位置提取结果的影响。

图9 添加不同比例噪声后，圆拟合、直接椭圆拟合和所提方法的光瞳位置提取结果Fig.9 The pupil position extraction result of circle fitting， direct ellipse fitting and proposed method after adding different proportions of noise

图10 添加不同比例噪声后，圆拟合、直接椭圆拟合和所提方法光瞳位置提取结果的MAX值Fig.10 The pupil position extraction MAX value of circle fitting， direct ellipse fitting and proposed method after adding different proportions of noise

图11 添加不同比例噪声后，圆拟合、直接椭圆拟合和所提方法光瞳位置提取结果的RMS值Fig.11 The pupil position extraction RMS value of circle fitting， direct ellipse fitting and proposed method after adding different proportions of noise

2.2 真实图像结果分析

根据项目需求，光瞳位置提取结果在标定位置的4.6个像素内，即满足项目需求。如图12所示，红色十字为标定位置，光瞳位置的偏移量在红圈内变化，即偏移量小于4.6个像素即满足项目需求。

图12 光瞳位置提取精度要求Fig.12 Accuracy requirements for pupil position extraction

对于真实图像，首先进行图像预处理，双边滤波窗口N=3，σr=22，σd=1.5，β=0.3，中值滤波窗口大小为15。考虑到RANSAC算法只有一定概率得到可信模型，将RANSAC算法的迭代参数K变为两倍。分别使用圆拟合、直接椭圆拟合和所提方法提取100帧光瞳图像位置，位置提取结果如表1所示，表中Position Offset AVG表示提取的光瞳位置与标定位置的偏移量的均值，图像标定位置为(296.1，277.4)，Position Offset MAX表示提取的光瞳位置的偏移量的最大值，Position Offset RMS表示光瞳位置的偏移量的波动情况，从表1中可以看出，直接椭圆拟合方法在偏移量均值、偏移量最大值方面均优于圆拟合的方法，所提方法在偏移量均值、偏移量最大值和偏移量的波动3方面优于直接椭圆拟合方法。

表1 不同方法光瞳位置提取结果Table 1 Results of pupil position extraction by different methods

表2中为3次提取真实光瞳图像的位置偏移量结果，从表2中可以看出，3次光瞳位置提取结果均满足光路对准的要求。

表2 所提方法光瞳位置提取结果Table 2 Proposed method results of pupil position extraction

直接椭圆拟合和所提方法一帧光瞳位置提取结果如图13所示，图13中，白色点是提取的内轮廓，蓝色点表示RANSAC算法计算出的内点，红色点是使用内点拟合出的椭圆，绿色点是椭圆拟合算法拟合出的椭圆。从图13中可以看出，直接椭圆拟合方法拟合出的椭圆受轮廓中的噪声影响，拟合出的椭圆向上和向右偏移，而使用所提方法拟合出的椭圆，受轮廓中噪声的影响较小。

图13 不同方法提取光瞳位置结果Fig.13 The pupil position results of different meshod

提出的RANSAC椭圆拟合方法分别在CPU平台和GPU平台上运行，算法运行的硬件平台为：联想隐士 X1笔记本上，Inter®Core™ i7-8750H CPU @2.20 GHz 2.21 GHz，16 GB 运行内存，显卡为 NVIDIA GeForce GTX 1050Ti，软件平台为Qt5.12.10。多次在CPU和GPU平台上运行，运行时间如图14所示。RANSAC椭圆拟合算法在CPU平台上运行平均时间为808.2 ms，在GPU平台上运行平均时间为14.2 ms，GPU时间加速比为56.9。所提方法总的运行平均时间为25.9 ms（38.6 Hz），满足实时性要求。

图14 CPU和GPU运行时间Fig.14 The running time of CPU and GPU

3 结论

面对光路自动准直系统中的实际需求，提出了基于RANSAC模型的椭圆拟合的光瞳位置提取方法。对所提方法的准确性和稳定性在仿真光瞳图像上进行验证，随着噪声的增加，所提方法结果的位置偏差AVG误差小于0.05 Pixel，对噪声具有明显的抑制作用。对于实际的光瞳图像，提取100帧真实光瞳图像的位置，所提方法的形心提取结果的最大值小于4.6 Pixel，满足系统对光瞳位置提取的要求。RANSAC算法只有一定的概率得到可信的模型，概率与迭代次数成正比，当噪声点比例明显增多时，需要增加迭代次数，迭代次数越多，得到正确结果的概率越大，相应的计算效率会下降。并且噪声点的比例受阈值的影响，因此要取得较好的效果，需要找到合理的阈值。