1.   引言

共光路斐索型动态干涉仪^[1]光机结构装调完毕后，由于在安装时不能保证空间相对位置绝对的精确，存在微量的旋转、上下平移、左右平移等问题；并且由于相机的摆放位置不准确，还有可能出现采集到的条纹图存在镜像、倒像等情况；另外在成像过程中不同的CMOS感光特性也有细微的差别；共光路动态干涉仪在进行面形重构时，相位信息需要从同步采集的4幅条纹图中获得，4幅条纹图中相同位置的像素点所对应的必须是光学面形上同一处细小区域反射光的干涉光强^[2]。解决以上问题的经典方法是进行图像配准，通过对图像进行旋转、平移等操作，使4幅不同相位调制量的条纹图中，相同位置的像素点对应同一空间位置。由于条纹图中的特征信息较少，单纯对采集到的条纹图进行配准，难度较大，并且配准精度很低。因此需要提出一种合适的配准方法来实现4幅条纹图的像素点点对应。

目前能够解决空间一致性问题的方法主要有两种：一种是基于条纹图的图像配准技术，通过对条纹图的特征点进行提取，然后通过配准条纹图像来实现4幅条纹图的点点位置精确对应，此种方法的缺点是，由于采集到的条纹图像几乎没有相同的特征点，特征点提取精度不高导致图像的配准精度很低；第二种是通过CMOS相机的标定技术，通过特殊的配准装置，使4台CMOS相机同时采集含有十字丝的刻划板，通过肉眼检测的方式来确定十字丝的位置重合标定，此方法属于人工检测装调技术，具有很大的操作误差，配准精度依然不高。

南京理工大学左芬^[3]、徐晨^[4]等人提出两种位置匹配方法：(1) 采用在条纹图中加入一个标志点来对4幅条纹进行配准，此方法的缺点是在条纹中人为加入标志点，会更改条纹像素值，最后会影响面型计算结果，并且将条纹与标志点混在一起也加大了提取标志点的难度。(2) 利用计算条纹图区域重心法来对4幅条纹图像进行位置一致性确定，此方法的缺点是有效条纹区域的确定本身就有误差，在有误差的区域内再次进行重心法计算，会继续增加误差，使得中心点计算结果更加不精确，并且此方法无法标定条纹图像的旋转量。内蒙古工业大学田枫等人^[5]提出一种基于统计分析方法的圆形域同步移相干涉图位置配准技术，其精度优势在于条纹图轮廓的精确划分，而通常条纹图对比度不高，影响条纹图轮廓提取，最终影响精度计算。

本文将综合以上方法的优缺点，提出一种新的配准思路，将条纹图与位置配准在关联度上人为分开，搭建了基于十字丝的专用图像配准装置，避免将条纹与标志点混在一幅图像内，然后利用最小二乘法对十字丝的两条直线进行拟合，配准得到CMOS相机的统一坐标系。此种方法的配准精度高，可以实现4幅条纹图的精确配准。

2.   基于整体最小二乘法的十字丝配准方法

2.1   基于刻划板十字丝图像的物理配准

由于相移分光产生的4个不同相位的干涉条纹由4个CMOS相机分别接收，因此在图像处理时要确保处理的图像像素要在严格的相同位置上，这就需要对CMOS相机进行位置配准，为此本文搭建了基于刻划板十字丝图像的物理配准装置。

本文搭建的干涉条纹配准结构图如图 1所示。调整过程如下：在镜头的目视位置安装十字刻划板，在镜头前端安装法线与镜头光轴平行的平面反射镜，通过均匀光源照亮刻划板后并通过镜头放大成像在4个CMOS相机上，使其中一个像在接近靶面中心位置后固定，调整其他3个CMOS相机上的像与第一个像位置重合，即保证了4个CMOS相机靶面位置空间位置的配准。在配准完成的CMOS相机进行空间位置固定之后，CMOS相机依然会出现不同影响的位置偏移及旋转，此时需要通过图像配准技术对十字丝图像进行精确配准。

图  1  干涉条纹配准结构三维图

Figure  1.  3D picture of inference fringe registration structure

下载: 全尺寸图片 幻灯片

利用4个CMOS相机同时采集十字刻划板图像，采集到的原始图像如图 2所示，然后通过算法计算十字丝交点坐标的偏移量以及十字丝的旋转量来实现四幅干涉条纹图的位置配准。

图  2  调整前刻划板采集图

Figure  2.  Collecting picture of depiction board before regulating depiction

下载: 全尺寸图片 幻灯片

2.2   基于sobel算子的图像预处理

由于原始图像十字丝特征不明显，首先我们需要对图 2中的原始图像进行边缘处理，提高十字丝的对比度，常见的边缘提取有Prewitt算子^[6]、拉普拉斯算子^[7]、candy算子^[8]等。由于这些算法局限性较高，不能保证提取目标的完整边缘，本文将利用Sobel算子对采集到的刻划板图像进行预处理，然后在其基础上利用整体最小二乘拟合对十字丝进行直线提取操作，最后联立直线方程，实现求出两直线的交点坐标。

Sobel算子本质上是一种梯度幅度检测算子^[6]，其原理是：以像素点为中心的3×3窗口内进行灰度计算。然后根据该点的极值状态来进行边缘提取检测。令f(x,y)为像素点的灰度值，f_x为水平方向上的梯度，f_y为垂直方向上的梯度。

Sobel边缘检测算子定义为：

Sobel边缘检测算子综合了图像每个像素点的上下左右邻点灰度的和，接近模板中心的权值。适当选取阈值T就可以实现图像边缘检测。

对0°相位的十字丝图像预处理结果如图 3所示。由结果可以看出，预处理后的十字丝图像对比度高，十字丝更加明显，对于之后的直线拟合起到关键性作用。

图  3  预处理结果及十字丝放大显示图

Figure  3.  Result of pretreatment and the enlarged image of cross wire

下载: 全尺寸图片 幻灯片

2.3   基于整体最小二乘法的直线拟合

在图像预处理之后，需要对图 3中的十字丝进行直线拟合，求取两直线交点的坐标，为此我们将求交点坐标转换为直线拟合问题。直线拟合问题可以描述为：假设给定n个坐标点，直线拟合的最终目的是通过计算得到一条最佳拟合直线，此直线能够尽量多的通过或者接近这些坐标点。其本质是求直线斜率和截距的最佳估计，最经典最常用的拟合方法是最小二乘法。普通最小二乘法^[9]的缺陷在于自变量和因变量的选择对直线拟合的结果影响很大，造成拟合精度差异较大。经典最小二乘法要求误差仅来源于模型中的因变量，解释变量无误差或相比因变量误差而言其误差大小可以忽略不计，而在很多工程应用中，模型中解释变量有误差且是不可忽略的，因此20世纪80年代提出了整体最小二乘法^[10](Total Least Squares)，其思想是建立在变量含误差EIV(Errors-In-Variables)模型上，能够很好地解决所有数据被随机误差影响的问题，其精度比经典最小二乘法更高。本文在最小二乘直线拟合方法的基础上，提出了基于整体最小二乘直线拟合方法，完成了十字丝两直线交点坐标的计算。

整体最小二乘法将自变量和因变量的误差加入到直线拟合方程中，其表达式为：

相应的误差方程可以按照EIV模型描述，如式(3) 所示:

式中，设计矩阵A与观测向量l的误差分别由E_A、E_l表示。

在设计矩阵A中，可以看到有一列元素的值恒为1，这是混合最小二乘法的求解问题^[11]，因此，应按混合最小二乘求解。

令：A=[A₁ A₂]，其中，A₁=[1 1 … 1]^T，A₂=[x₁ x₂ … x_m]^T，此时构造一个增广矩阵C，并对其进行QR三角分解：

式中，Q为正交矩阵，R为上三角矩阵，则:

式中，R₁₁、R₁₂、R₂₂、R_1l、R_2l为标量，可将上式写成如下的形式：

在对上述参数进行求解时，先利用整体最小二乘法对方程(7) 进行求解以获取参数δa，把δa带入方程(6) ，然后再按照普通的最小二乘法进行求解参数δb。

构造增广矩阵C_R=[R₂₂ R_2l]，并进行奇异值分解：

式中，U=[U₁ U₂]，V=[N₁ N₂]，Σ=diag(σ₁,σ₂)，σ₁＞σ₂。

可得参数δa的整体最小二乘解如下式所示：

将δ$\hat a$带入方程(6) 可得δ$\hat b$如下式：

得到了参数δa，δb的解，需要对其残差和精度进行理论上的评价。

把增广矩阵C_R写成下述形式：

R的改正量即为：

设计矩阵和观测向量的改正数为：

此时可得出残差：

在正交方向上的距离残差为：

距离残差在几何意义上就是最小的正交距离的残差平方和，所以，样本精度的评定应该按照正交距离的残差的平方和来计算:

由于整体最小二乘法同时考虑了两个坐标轴的误差，并且在两个坐标轴方向上的正交距离残差都比普通最小二乘法小，其单位权重误差比普通的最小二乘法更小。

3.   实验结果与分析

本文搭建了共光路动态干涉仪的CMOS相机图像配准试验平台，平台选择KEYENCE公司的VH-Z50L长距离高性能变焦镜头，其观察距离在85 mm时能够达到500的放大率，该装置能够辅助配准干涉条纹。在缩短CMOS相机像面和PBS距离后，能够选取放大率更大的镜头。利用4个CMOS相机同时拍摄刻划板实物，刻划板上包含一个十字丝标志。表 1为利用整体最小二乘法得到的配准结果，包含平移量和旋转量。图 4为经过整体最小二乘配准后的十字丝刻划板图，图 5为配准完成后利用配准得到的旋转平移公式变换后的4幅干涉条纹图。

图  4  整体最小二乘法配准后刻划板图

Figure  4.  Depiction board picture after registration using total least square

下载: 全尺寸图片 幻灯片

图  5  按配准得到的公式变换后的条纹图

Figure  5.  Fringe picture of transform using the registration result

下载: 全尺寸图片 幻灯片

表  1  整体最小二乘法得到配准结果

Table  1.  Registration result using total least square

图像序列	十字丝中心点	平移量	直线斜率	旋转量/(°)
0相位	(1 030,1 045)	(6,21)	0.013 6	0.779 2
180相位	(1 016,1 021)	(-8,-3)	0.001 3	0.074 5
90相位	(1 023,1 025)	(-1,1)	0.000 8	0.045 8
270相位	(1 019,1 027)	(-5,3)	0.001 1	0.063 0

下载: 导出CSV 

| 显示表格

由图 2可以看出，初始采集到的4幅十字丝图像具有明显的位置差异，利用本文提出的整体最小二乘法得到4个十字丝中心点位置坐标，并将四幅配准图的中心点平移到(1 024 ,1 024) 坐标处，然后利用整体最小二乘法得到四幅配准图的旋转量，以0°相位配准图为标准，将其他3幅配准图转换到0°相位配准图的位置，表 2计算出配准后的4幅条纹图两两互相关度结果表明，经过整体最小二乘方法进行配准后，两两条纹图之间的互相关度结果达到了0.97左右。结果表明，本算法达到了预期的处理目的。

表  2  互相关度计算结果

Table  2.  Result of cross-correlation value

相机标号	配准前	整体最小二乘法	重心法
(1,2)	0.15	0.96	0.83
(1,3)	0.17	0.97	0.84
(1,4)	0.14	0.97	0.82
(2,3)	0.12	0.98	0.92
(2,4)	0.10	0.96	0.91
(3,4)	0.15	0.97	0.90

下载: 导出CSV 

| 显示表格

为了更好的验证本算法的优越性，利用重心法配准结果^[12-15]与本文提出的配准算法进行比较。图 6为重心法条纹图手动提取模板，图 7为利用提取模板得到的有效条纹区域，图 8为配准结果，由表 2可以看出，重心法由于只考虑了条纹中心配准结果，未考虑图像的旋转量，其互相关度结果只在0.92以下，较本文算法有一定的差距。

图  6  手动提取模板

Figure  6.  Cyclostyle by manpower

下载: 全尺寸图片 幻灯片

图  7  有效区域提取结果

Figure  7.  Results collection of valid region

下载: 全尺寸图片 幻灯片

图  8  重心法配准结果

Figure  8.  Registration result of barycenter method

下载: 全尺寸图片 幻灯片

4.   结论

本文首先介绍了条纹图位置关系对检测结果造成的影响，在此基础上提出采用硬件装调与图像处理相结合的方法完成共光路动态干涉仪的条纹图配准工作。以硬件装调为主，使得对4幅条纹图进行的数字信号处理操作尽可能少，以保证其配准的同时减少信息丢失。在目前国内的光机加工水平下，在安装时，难以保证空间相对位置绝对的精确，仍然会存在微量的旋转、俯仰、横摆、上下平移、左右平移等问题，本文采用整体最小二乘法来对十字丝进行配准，提高条纹图的配准精度。最后通过一系列的试验验证了本文算法的正确性，并与重心配准法进行了性能比较，验证了本文算法的优越性。