在航空工业中，装配连接孔的质量对飞机安全性和使用寿命有至关重要的影响^[1-2]。随着航空数字化制造的不断发展，视觉测量技术在机器人自动钻铆系统中得到广泛应用，一方面通过检测基准孔的偏差进行位置补偿，提高机器人制孔的位置精度^[3-4]；另一方面能够替代传统的手工测量方式，实现连接孔的高效实时检测^[5]。

边缘检测作为圆孔视觉测量过程的核心技术^[6]，直接影响了圆孔的在线检测精度。目前，常用的像素级边缘检测算法包括一阶微分方法中的Sobel算子、Prewitt算子及Roberts算子^[7]等和二阶微分方法中的Laplacian算子及Canny算子^[8]等。但是，传统像素级检测算子的定位精度已不能满足数字化制造中高精度测量的要求。因此，国内外学者在像素细分的理论基础上，不断深入地开展了亚像素边缘检测算法的研究。其中，矩方法相比拟合法和插值法具有对噪声不敏感和效率高等优点^[9]，因此得到了广泛关注。Lyvers ^[10]等提出了利用计算空间矩的方法提取边缘，通过理想阶跃边缘模型的空间矩计算边缘参数。Ghosal ^[11]等首次提出了利用Zernike正交矩(Zernike orthogonal moment, ZOM)来进行边缘检测，该算法只需要3个Zernike矩就可计算出边缘参数，提高了边缘检测的效率。李金泉等^[12]通过考虑模板效应改进了传统的Zernike算法，细化了图像提取的边缘。高世一等^[13]在Ghosal算法中模板的基础上，计算了模板系数并推导了边缘参数的表达式，进而提高了边缘检测的精度。张晓林等^[14]利用几何信息，仅通过Zernike零阶矩获得亚像素边缘，减少了坐标点的计算量。传统Zernike算法在理想阶跃模型的基础上进行亚像素边缘检测，然而实际的检测目标和背景的交界边缘处存在过渡区，理想阶跃边缘模型与实际的边缘灰度分布存在误差。通过分析边缘的灰度过渡区，国内外学者建立了边缘的三灰度过渡模型^[15-16]，推导了边缘参数的原理误差补偿函数。这些研究在利用三灰度过渡模型进行边缘参数推导过程中仍然采用了理想阶跃模型的参数计算公式，虽然通过原理误差补偿的方式获取了实际的参数，但是增加了计算成本。

在综合分析边缘检测相关文献的基础上，本文采用像素级粗定位和亚像素级精定位相结合的方式，准确推导了三灰度过渡模型的边缘参数，优化了边缘检测灰度阈值的条件，提出了一种基于改进Zernike矩的亚像素钻铆圆孔检测算法。仿真和实验的对比表明，本文算法比传统算法具有更高的检测精度，适应于不同尺寸的钻铆圆孔检测，并且对噪声的敏感程度较低，在机器人钻铆视觉测量系统中具有重要的应用意义。

1 Zernike矩边缘检测原理

Teague根据Zernike正交实值多项式，定义了灰度图像n阶m次的Zernike矩^[17]：

$ {G_{nm}} = \frac{{n + 1}}{\pi }\iint\limits_{{x^2} + {y^2} \leqslant 1} {f\left( {x,y} \right)V_{nm}^ * \left( {\rho ,\theta } \right){\text{d}}x{\text{d}}y}. $

(1)

其中：f(x, y)是图像中点(x, y)的灰度值，V_nm^*(ρ, θ)是Zernike多项式V_nm(ρ, θ)的复共轭。在离散条件下，图像的二维ZOM在x²+y²≤1的条件下可表达为

$ {G_{nm}} = \sum\limits_x {\sum\limits_y {f\left( {x,y} \right)V_{nm}^ * \left( {\rho ,\theta } \right)} } . $

(2)

亚像素边缘定位依据是Zernike矩的旋转不变性，如果将图像逆时针旋转α，旋转前后的图像G_nm′和G_nm的关系可表示为

$ {G_{nm}} = {G'_{nm}}{{\text{e}}^{ - {\text{j}}m\alpha }}. $

(3)

则旋转前后图像ZOM的模不变，只是相角发生了改变，利用旋转后的图像ZOM可以相对容易地计算出图像的边缘参数，从而实现对边缘的亚像素定位。

Ghosal ^[11]算法建立了理想的边缘阶跃模型，如图 1所示，h为图像中背景的灰度值，t为背景与目标之间的阶跃灰度，在单位圆中，L为目标的边缘线，d为原点到边缘线的距离，α为d与x轴的夹角。

该算法在计算模板的基础上，利用3个积分核函数推导了图像的ZOM模板系数，移动模板与待检图像进行卷积运算，得到图像的ZOM和边缘参数，根据边缘参数的阈值条件判断该像素点是否为边缘点。在考虑计算模板放大效应^[12]的基础上，检测得到的边缘点坐标(x_m, y_m)修正公式为

$ \left[ {\begin{array}{*{20}{c}} {{x_m}}\\ {{y_m}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{x_{\rm{e}}}}\\ {{y_{\rm{e}}}} \end{array}} \right] + \frac{{Nd}}{2}\left[ {\begin{array}{*{20}{c}} {\cos \alpha }\\ {\sin \alpha } \end{array}} \right]. $

(4)

其中：(x_e, y_e)为原点坐标，N为模板大小。

2 改进的亚像素边缘检测模型 2.1 理想阶跃模型的修正

Ghosal算法基于边缘的理想阶跃模型推导了亚像素坐标的公式，实际中由于非理想的离散采样和光学系统衍射效应的限制^[18]，在待检测的圆孔边缘附近存在过渡区，即图像的灰度值从背景逐渐变为圆孔。为了描述与实际过渡区相符的灰度，国内外学者建立了边缘的三灰度过渡模型^[15-16]，如图 2所示，虚线代表边缘的理想阶跃模型，实线代表三灰度过渡模型，其中h是背景灰度值，h+Δt是过渡区灰度值，(h+t)是目标灰度值，d₁和d₂分别是中心与过渡区左右界限的距离。其中，d₁和d₂的距离差为a，即模板中每个像素的距离，本文采用模板7×7，因此a=2/7。

根据图 2所示的亚像素边缘检测三灰度模型，推导各个ZOM的表达式，本文计算了旋转后的4个ZOM，分别如式(5)—(8)所示。

$ \begin{array}{*{20}{c}} {{G_{00}} = \Delta t\left[ {{d_2}\sqrt {1 - d_2^2} + {{\sin }^{ - 1}}\left( {{d_2}} \right) - } \right.}\\ {\left. {{d_1}\sqrt {1 - d_1^2} - {{\sin }^{ - 1}}\left( {{d_1}} \right)} \right] + }\\ {t\left[ {\frac{{\rm{ \mathsf{ π} }}}{2} - {d_2}\sqrt {1 - d_2^2} - {{\sin }^{ - 1}}\left( {{d_2}} \right)} \right] + h{\rm{ \mathsf{ π} ,}}} \end{array} $

(5)

$ \begin{array}{*{20}{c}} {{G_{11}} = \frac{2}{3}\Delta t\left[ {{{\left( {1 - d_1^2} \right)}^{\frac{3}{2}}} - {{\left( {1 - d_2^2} \right)}^{\frac{3}{2}}}} \right] + }\\ {\frac{2}{3}t{{\left( {1 - d_2^2} \right)}^{\frac{3}{2}}},} \end{array} $

(6)

$ \begin{array}{*{20}{c}} {{G_{20}} = \frac{2}{3}\Delta t\left[ {{d_1}{{\left( {1 - d_1^2} \right)}^{\frac{3}{2}}} + {d_2}{{\left( {1 - d_2^2} \right)}^{\frac{3}{2}}}} \right] + }\\ {\frac{2}{3}t{d_2}{{\left( {1 - d_2^2} \right)}^{\frac{3}{2}}},} \end{array} $

(7)

$ \begin{array}{*{20}{c}} {{G_{31}} = \frac{2}{3}\Delta t\left[ {{{\left( {1 - d_1^2} \right)}^{\frac{3}{2}}} - {{\left( {1 - d_2^2} \right)}^{\frac{3}{2}}}} \right] + }\\ {\frac{4}{5}\Delta t{{\left( {1 - d_2^2} \right)}^{\frac{5}{2}}} - \frac{4}{5}\Delta t{{\left( {1 - d_1^2} \right)}^{\frac{5}{2}}} + }\\ {\frac{2}{{15}}t{{\left( {1 - d_2^2} \right)}^{\frac{5}{2}}}\left( {6d_2^2 - 1} \right).} \end{array} $

(8)

本文在计算G₀₀、G₁₁、G₂₀、G₃₁的基础上，准确推导了该模型的边缘参数d、Δt、t、α，为获得钻铆圆孔的亚像素坐标提供了计算依据。由式(3)推导可得旋转前后G′_n1和G_n1的关系：

$ \begin{array}{*{20}{c}} {{G_{n1}} = \cos \left( {{\alpha _{n1}}} \right){\mathop{\rm Re}\nolimits} \left[ {{{G'}_{n1}}} \right] + \sin \left( {{\alpha _{n1}}} \right){\mathop{\rm Im}\nolimits} \left[ {{{G'}_{n1}}} \right] + }\\ {\left( {\sin \left( {{\alpha _{n1}}} \right){\mathop{\rm Re}\nolimits} \left[ {{{G'}_{n1}}} \right] - \cos \left( {{\alpha _{n1}}} \right){\mathop{\rm Im}\nolimits} \left[ {{{G'}_{n1}}} \right]} \right){\rm{i}}{\rm{.}}} \end{array} $

(9)

其中Re[G′_n1]和Im[G′_n1]分别代表G′_n1的实部和虚部。旋转后图像关于x轴对称，因此可以得到式(10)和边缘线上对应旋转角度的表达式(11)。

$ \sin \left( {{\alpha _{n1}}} \right){\mathop{\rm Re}\nolimits} \left[ {{{G'}_{n1}}} \right] - \cos \left( {{\alpha _{n1}}} \right){\mathop{\rm Im}\nolimits} \left[ {{{G'}_{n1}}} \right] = 0. $

(10)

$ {\alpha _{n1}} = {\tan ^{ - 1}}\left( {\frac{{{\mathop{\rm Im}\nolimits} \left[ {{{G'}_{n1}}} \right]}}{{{\mathop{\rm Re}\nolimits} \left[ {{{G'}_{n1}}} \right]}}} \right). $

(11)

联立式(5)—(8)可以推导得到准确的边缘参数d₁、d₂、Δt、t的表达式：

$ \left\{ \begin{array}{l} {d_1} = \frac{{6{G_{20}} - 3a{G_{11}} + \sqrt { - 6G_{11}^2 - 30{G_{31}}{G_{11}} + 36G_{20}^2 + 9G_{11}^2{a^2}} }}{{6{G_{11}}}},\\ t = \frac{{\frac{3}{2}{G_{11}} - \Delta t\left[ {{{\left( {1 - d_1^2} \right)}^{\frac{3}{2}}} - {{\left( {1 - d_2^2} \right)}^{\frac{3}{2}}}} \right]}}{{{{\left( {1 - d_2^2} \right)}^{\frac{3}{2}}}}},\\ \Delta t = \frac{{3\left( {{d_2}{G_{11}} - {G_{20}}} \right)}}{{2a{{\left( {1 - d_1^2} \right)}^{\frac{3}{2}}}}},\\ {d_2} = {d_1} + a. \end{array} \right. $

(12)

通过图 2三灰度过渡模型的几何关系，边缘参数d与d₁、d₂、Δt、t的关系可以表示为

$ d = {d_2} - \frac{{\Delta t}}{t}\left( {{d_2} - {d_1}} \right). $

(13)

本文与之前边缘灰度模型的不同之处在于：1)考虑了边缘过渡区对亚像素定位精度的影响；2)通过7×7模板的ZOM准确推导了三灰度过渡模型的边缘参数，避免了采用阶跃模型计算三灰度过渡模型的参数所带来的原理误差及补偿计算；3)与传统7×7模板的ZOM计算相比，减少了G₄₀的计算量，提高了检测效率。

2.2 相对灰度阈值条件的优化

传统Ghosal算法的边缘点由边缘参数d和t确定，其阈值条件为

$ t \ge {t_{\rm{z}}} \cap d \le {d_{\rm{z}}}. $

(14)

其中t_z和d_z分别是相对灰度的阈值和距离的阈值。2d_z必须小于一个像素的长度，因此d_z的选择范围较小，传统算法中d_z=1/$\sqrt 2 $，考虑模板效应后距离阈值的条件应为d≤2d_z/N，其中N为模板大小7；但是t_z的选择范围较大，如果相对灰度的阈值太小会获取很多的伪边缘，而阈值过高虽然可以有效地抑制伪边缘，但同时也会丢失圆孔的真实边缘信息。为了保证本文算法对不同圆孔图像的检测精度的稳定性，需要根据图像优化t_z，这样不仅避免手动调节的人为因素影响，而且提高了边缘检测的效率，实现了钻铆圆孔的高精度在线实时检测。因此，本文采用了Otsu阈值法^[19]，利用图像的灰度值，以圆孔和背景的灰度值类间方差为测度准则，当该准则函数取最大值时图像分割效果最好，即可得到最优的相对灰度分割阈值。通过该算法进行钻铆圆孔的拍摄图像分割，从而得到最优的相对灰度阈值t_z。

3 钻铆圆孔的亚像素检测方法

在改进的Zernike矩亚像素边缘模型的基础上，本文提出了钻铆圆孔的亚像素检测方法。飞机部件的钻铆连接孔大多为排孔或群孔，为了减少亚像素定位的计算量，提高边缘检测效率，本文首先采用Canny算子进行圆孔边缘的像素级粗定位，快速获得边缘点坐标信息，生成像素级边缘点集；其次提取待测圆孔所在的感兴趣区域(region of interest，ROI)；然后基于改进的Zernike矩亚像素边缘模型对圆孔边缘进行精定位，获得准确的边缘点亚像素坐标；最后由于视觉测量中圆孔的成像为椭圆孔^[20]，根据最小二乘原理建立椭圆参数的优化目标函数，采用Levenberg-Marquardt非线性优化算法求解钻铆圆孔的中心点亚像素坐标和半径参数，从而实现了基准孔的定位和连接孔的参数在线测量。其算法流程图如图 3所示。

本文将该方法的主要步骤概述如下：

步骤1 将工业相机拍摄的原始待检图像进行灰度化处理，采用Gauss滤波消除噪声影响，以实现待检图像的预处理。

步骤2 计算各像素点的梯度幅值与方向，通过Canny算子实现图像中圆孔的粗定位，得到像素级边缘点集{D}，计算图像中心坐标(x₀, y₀)及圆半径r₀；以(2r₀+m)为边长提取ROI，其中m为某一较小的常数。

步骤3 分别计算Zernike矩的7×7模板系数M₀₀、M₁₁、M₂₀、M₃₁^[13]，将其与{D}中的边缘点卷积运算得到G₀₀、G₁₁、G₂₀、G₃₁；与此同时，在ROI中采用Otsu自适应算法优化相对灰度的阈值t_z。

步骤4 选取坐标点集{D}中的第i个边缘点D_i，根据改进的亚像素检测模型求解其相应的边缘参数d_i、t_i和α_i。

步骤5 判断D_i的边缘参数是否满足阈值条件t_i≥t_z∩d_i≤d_z，若不满足则剔除D_i并检测下一个边缘点；若满足则进一步判断‖D_i－D_is‖≤ε是否成立，不成立则剔除坐标超差的边缘点，其中ε取0.5，满足条件时根据式(15)计算其亚像素坐标D_is(x_is, y_is)并保存至点集{D_s}。

$ \left[ {\begin{array}{*{20}{c}} {{x_{i{\rm{s}}}}}\\ {{y_{i{\rm{s}}}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{x_0}}\\ {{y_0}} \end{array}} \right] + \frac{{N{d_i}}}{2}\left[ {\begin{array}{*{20}{c}} {\cos {\alpha _i}}\\ {\sin {\alpha _i}} \end{array}} \right]. $

(15)

步骤6 按照上述算法遍历{D}，生成亚像素点集{D_s}={(x_1s, y_1s), (x_2s, y_2s), …, (x_is, y_is), …, (x_ns, y_ns)}。根据椭圆的一般方程建立最小二乘的优化目标函数，如式(16)所示，进而求解椭圆参数得到圆心坐标(x, y)及长短轴。通过计算长短半轴的平均值，得到圆孔的半径参数r。

$ \begin{array}{*{20}{c}} {T = \sum\limits_{is = 1}^n {\left( {{T_1}x_{i{\rm{s}}}^2 + {T_2}{x_{i{\rm{s}}}}{y_{i{\rm{s}}}} + } \right.} }\\ {{{\left. {{T_3}y_{i{\rm{s}}}^2 + {T_4}{x_{i{\rm{s}}}} + {T_5}{y_{i{\rm{s}}}} + 1} \right)}^2}.} \end{array} $

(16)

4 仿真与实验 4.1 仿真运算

为了验证本文提出算法的圆孔检测有效性，人工仿真制作了61×61和101×101像素的图像，如图 4所示，对传统Zernike算法和本文算法进行仿真对比。

表 1和2分别为图 4a和图 4b的检测仿真数据，记录了仿真图中圆形的边缘点坐标、圆心坐标和半径。通过对比分析检测数据和相应理想数据的误差，传统算法的边缘点相对误差小于0.5像素，圆心坐标的相对误差在0.02像素范围内，半径的相对误差在0.2像素内；本文算法的边缘点相对误差小于0.25像素，圆心坐标的相对误差在0.01像素范围内，半径的相对误差在0.1像素内。因此，本文算法的边缘点亚像素坐标、圆心坐标和半径的偏差比传统Zernike算法相应参数的偏差小。此外，传统Zernike算法和本文算法的计算时间大致相同，说明在不增加计算时间成本的基础上，本文算法具有更高的检测精度。

4.2 实验测试

为验证本文算法的稳定性和抗噪性，搭建了视觉检测实验装置，由200万像素的工业相机(感光尺寸为1/2.7英寸)、工业相机支架、环形光源、电脑和铝合金试验件构成，如图 5所示。图 6a和6b分别是平面和曲面试验件，其中平面试验件的尺寸为100 mm×100 mm，曲面试验件的尺寸为70 mm×40 mm，试验件上预制了不同尺寸的待检测圆孔，由于待钻铆的航空部件曲率较小，则曲面试验件设计了较大的曲率半径。通过相机和像素当量的标定得到，在该条件下像素当量为每像素0.17 mm，据此计算圆孔的实际中心坐标和半径参数。

第1组实验为算法的稳定性测试，由于飞机部件上待钻铆的连接孔具有不同的尺寸，因此需要测试本文算法对不同连接孔的检测精度稳定性。根据最小二乘法拟合圆后，较难对比出Zernike算法改进前后的效果，因此本文从对比改进和传统Zernike算法的边缘点提取图。如图 7—8中的方框所示，传统Zernike算法虽然可以检测出圆孔的边缘轮廓，但是与本文算法相比出现了断裂，丢失了一些关键边缘点信息，降低了亚像素边缘检测的精度。相比传统Zernike算法，本文算法对不同尺寸的圆孔具有较稳定的提取精度。

第2组实验为该算法的抗噪性测试，由于算法在图像预处理阶段采用了Gauss滤波，因此本文主要针对椒盐噪声进行测试，通过加入不同密度的椒盐噪声，对比亚像素检测的不同圆孔参数，其中椒盐噪声的密度是指噪声所占当前图像全部像素的百分比。以相同理想尺寸的圆孔3和7为例，圆心坐标x、y和圆孔半径r的误差定义为有噪声影响和无噪声的相应参数差值。对于在椒盐噪声的影响下平面试验件，r的误差在0.03 mm范围内，x或y的误差小于0.015 mm，如图 9所示。对于曲面试验件，r的误差在0.03 mm范围内，x的误差小于0.04 mm，y的误差小于0.02 mm，如图 10所示。因此，本文算法适用于平面或曲面待检测部件，对于一般的视觉测量环境具有较稳定的抗噪能力和检测精度。

5 结论

在机器人钻铆视觉测量系统中，针对钻铆圆孔检测精度难以满足高精度检测要求的问题，本文在Zernike矩边缘检测原理的基础上，推导了准确的亚像素边缘参数，从而改进了三灰度边缘过渡模型；通过Otsu算法优化了边缘检测的相对灰度阈值，并提出了一种基于改进Zernike矩的亚像素钻铆圆孔检测方法。仿真结果表明，相比传统Zernike算法，在相同计算效率的前提下，本文算法对钻铆圆孔的中心坐标和半径参数具有更高的检测精度；平面和曲面试验件的检测结果说明对于不同尺寸的钻铆圆孔，本文算法具有较稳定的提取精度，而且检测的圆孔参数对椒盐噪声的敏感程度较低，能够适用于一般的视觉测量环境，满足机器人钻铆过程中快速和高精度的钻铆圆孔检测要求，在机器人钻铆视觉测量系统中具有广阔的应用前景。