某型器物量测系统原型通过两个平移轴和一个旋转轴来控制变化摄像机位置与姿态以对器物拍照，并依据同名像点还原器物表面位置关系进行相关证据量测，系统结构如图 1所示。系统设计要求是1.2 m作业范围内的点位误差小于1 mm，因此，需要对系统进行误差检测和建模补偿。国内外学者已对误差量测方法与误差建模进行了大量研究，并取得了可喜的成果^[1-5]。在误差量测方面主要有以下方法：①直接测量法。用千分表、标准平尺等直接测量并给出误差值，该方法对专业性要求较高且测量时间长。②综合测量辨识法。使用激光干涉仪^[6]、球杆仪和激光跟踪仪^[7]等检测误差，利用误差传递原理^[8]分别将运动轴的各个部分所产生的误差联合起来。目前应用较多的是激光干涉仪，并在此基础上发展了9线法^[9]、22线法^[10]等方法，然而其误差辨识项较多，且成本太高。对于误差建模，通常是以多体系统理论为基础，采用刚体运动学与齐次坐标矩阵描述空间误差的传递关系^[11-12]，模型参数较多(以三轴为例，其参数多于20个)是阻碍其实际应用的主要问题^[13-14]。

研究表明，几何误差和热误差约占量测系统中运动轴总体误差的70%，当温度变化较小时，可以认为主要误差源为几何误差，其相对稳定，易于补偿^[15]。同时顾及系统自身的量测属性，本文未考虑上述常规方法，而是基于相机标定原理^[16-18]与多项式误差拟合模型提出了一种利用同名像点位姿变化建立系统误差补偿的新方法。该方法利用实际位姿变化后棋盘格标定板图像的角点像素坐标与其对应的拟合值建立同名像点，通过最小化同名像点间的差值，迭代求解得各轴运动误差。类似地，可获取多组运动误差数据，利用非线性最小二乘求得其多项式运动误差补偿模型。该方法以自校验方式替代测量仪器，成本较低，操作简易。此外，以几何误差作为最主要误差源，采用多项式模型，参数少，模型简单，适用性强，且弱化了各轴间的误差传递。

1 误差补偿模型新方法 1.1 误差补偿模型建立

系统结构示意图如图 1所示。将摄像机放置于旋转轴上，使其可以通过联合各个轴做位姿变化，用S-X_CY_CZ_C表示像空间坐标系，其中摄影中心S为像空间坐标系原点，摄像机主光轴为像空间坐标系的Z_C轴。将棋盘格标定板固定在与运动轴平面中心相对位置处，以其左上角内角点作为物空间坐标系的原点，以原点为起点的棋盘格横向线和纵向线分别为物空间坐标系的X轴和Y轴。虚线处为位置1，实线处为位置2，将摄像机在位置1、2处所拍包含标定板的像片分别记为左、右像片，保证至少有一个位置处的像空间坐标系的各个轴分别与各运动轴平行，两位置间的旋转和平移分别用R、T表示。其中R由绕竖直方向旋转的角度φ表示，T由两平移轴方向的运动值(x，y)表示。位置1与位置2处像空间坐标系与物空间坐标系之间的位姿关系分别用R_l、T_l与R_r、T_r表示。

相较于分析系统结构，利用严格几何成像模型及仿射变换模型可间接得到各运动轴之间的误差传递关系。本文直接利用非线性最小二乘拟合建立各运动值与相应误差值之间的多项式模型，模型参数无实际物理意义，可实现参数透明化，形式简单且可回避误差的传递过程。

根据测量系统中运动轴形式，可将误差补偿模型表示为旋转轴引起的二阶非线性误差、平移轴引起的一阶线性误差以及干扰误差(常数项)，表达式如下：

$ \left[ \begin{array}{l} \Delta \varphi \\ \Delta x\\ \Delta y \end{array} \right] = \left[ \begin{array}{l} {s_{00}}\ {s_{01}}\ {s_{02}}\\ {s_{10}}\ {s_{11}}\ {s_{12}}\\ {s_{20}}\ {s_{21}}\ {s_{22}} \end{array} \right]\left[ \begin{array}{l} {\varphi ^2}\\ {x^2}\\ {y^2} \end{array} \right] + \left[ \begin{array}{l} {f_{00}}\ {f_{01}}\ {f_{02}}\\ {f_{10}}\ {f_{11}}\ {f_{12}}\\ {f_{20}}\ {f_{21}}\ {f_{22}} \end{array} \right]\left[ \begin{array}{l} \varphi \\ x\\ y \end{array} \right] + \left[ \begin{array}{l} {e_0}\\ {e_1}\\ {e_2} \end{array} \right] $

(1)

式中，(Δφ，Δx, Δy)分别表示旋转轴、水平轴和竖直轴的运动误差；(φ，x，y)分别表示系统在旋转轴、水平轴和竖直轴的位姿变化量；s₀₀~s₂₂、f₀₀~f₂₂、e₀~e₂分别表示模型待求二阶、一阶系数及其常数项。

在确定误差补偿模型后，可利用图 1、图 2中的流程确定模型中的未知系数，具体过程可参见后述内容。

1.2 左像片内外方位元素的求解

如图 1所示，在位置1处拍摄左像片，将共线条件方程式^[19]表示为矩阵形式：

$ \left[ \begin{array}{l} x\\ y\\ 1 \end{array} \right] = s\left[ \begin{array}{l} - {f_x}\ 0\ {x_0}\\ 0\ - {f_y}\ {y_0}\\ 0\ \ \ \ 0\ \ \ \ 1 \end{array} \right]\left[ \begin{array}{l} {a_1}\ {b_1}\ {c_1}\ {t_x}\\ {a_2}\ {b_2}\ {c_2}\ {t_y}\\ {a_3}\ {b_3}\ {c_3}\ {t_z} \end{array} \right]\left[ \begin{array}{l} X\\ Y\\ Z\\ 1 \end{array} \right] $

(2)

式中，(x，y)表示左像片中棋盘格角点的像素坐标；s为尺度因子；(X，Y，Z)表示像点所对应的物像空间坐标；a_i、b_i、c_i(i=1，2，3)是像片外方位角元素φ、ω、κ的方向余弦，表示从物空间坐标系到像空间坐标系的旋转矩阵；t_x、t_y、t_z是像片外方位线元素，表示摄影中心S在像空间坐标系里的坐标。

考虑到求解像片内方位元素的过程中关注点仅为所需的平面坐标，对式(2)进行简化，令所选取物像空间坐标(即棋盘格标定板平面)的Z=0，即：

$ \left[ \begin{array}{l} x\\ y\\ 1 \end{array} \right] = s\mathit{\boldsymbol{A}}\left[ {{\mathit{\boldsymbol{r}}_1}\ \ {\mathit{\boldsymbol{r}}_2}\ \ {\mathit{\boldsymbol{r}}_3}\ \ \mathit{\boldsymbol{t}}} \right]\left[ \begin{array}{l} X\\ Y\\ 0\\ 1 \end{array} \right] = s\mathit{\boldsymbol{A}}\left[ {{\mathit{\boldsymbol{r}}_1}\ \ {\mathit{\boldsymbol{r}}_2}\ \ \mathit{\boldsymbol{t}}} \right]\left[ \begin{array}{l} X\\ Y\\ 1 \end{array} \right] $

(3)

随后利用张正友标定法^[18]，即可求得左像片的内外方位元素。

若假定标定板上角点附近的噪声服从高斯分布，则利用最小二乘优化对式(4)求最小值，即得待求参数的极大似然估计值：

$ \min \sum\limits_{i = 1}^n {\sum\limits_{j = 1}^m {{{\left\| {{m_{ij}} - \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over m} (\mathit{\boldsymbol{A}}, {\mathit{\boldsymbol{R}}_i}, {\mathit{\boldsymbol{t}}_i}, {\mathit{\boldsymbol{M}}_j})} \right\|}^2}} } $

(4)

式中，m_ij表示第i幅图像上第j个点的像素坐标；${\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over m} (\mathit{\boldsymbol{A}}, {\mathit{\boldsymbol{R}}_i}, {\mathit{\boldsymbol{t}}_i}, {\mathit{\boldsymbol{M}}_j})}$表示对物像点M_j通过所求单应性矩阵得到的在第i幅图像上的投影。至此可得到优化后左像片的内、外方位元素。

1.3 观测-拟合法求左右像片间位姿参数

如图 1所示，摄像机随运动轴从位置1到位置2，使位置2处于像空间坐标系的X_C轴、Y_C轴，分别与水平、竖直方向的运动轴平行，位姿关系用R、T表示。将控制轴运动所输脉冲数分别作为其初始值R₀、T₀，上节中已求出从物空间坐标系到位置1处像空间坐标系的位姿关系R_l、T_l，则根据坐标转换原理可得：

$ {\mathit{\boldsymbol{P}}_r} = {\mathit{\boldsymbol{R}}_0}{\mathit{\boldsymbol{P}}_l} + {\mathit{\boldsymbol{T}}_0} $

(5)

根据式(2)得：

$ \left\{ \begin{array}{l} {\mathit{\boldsymbol{P}}_l}{\rm{ = }}{\mathit{\boldsymbol{R}}_l}{\mathit{\boldsymbol{P}}_w}{\rm{ + }}{\mathit{\boldsymbol{T}}_l}\\ {\mathit{\boldsymbol{P}}_r}{\rm{ = }}{\mathit{\boldsymbol{R}}_r}{\mathit{\boldsymbol{P}}_w} + {\mathit{\boldsymbol{T}}_r} \end{array} \right. $

(6)

式中，P_l、P_r分别表示物像点在左、右像空间坐标系中的坐标；P_w表示物像空间坐标；R_r、T_r表示从物空间坐标系到位置2处像空间坐标系的位姿关系。

将式(6)代入式(5)可得：

$ \left\{ \begin{array}{l} {\mathit{\boldsymbol{R}}_r} = {\mathit{\boldsymbol{R}}_0}{\mathit{\boldsymbol{R}}_l}\\ {\mathit{\boldsymbol{T}}_r} = {\mathit{\boldsymbol{R}}_0}{\mathit{\boldsymbol{T}}_l} + {\mathit{\boldsymbol{T}}_0} \end{array} \right. $

(7)

可得到R_r、T_r的初始值。仅考虑绕竖直轴的旋转，有：

$ {\mathit{\boldsymbol{R}}_0} = \left[ \begin{array}{l} \cos \varphi\ \ 0\ \ sin \varphi \\ 0\ \ \ \ \ \ \ \ 1\ \ \ \ \ \ \ \ 0\\ - \sin \varphi\ \ 0\ \ cos \varphi \end{array} \right] $

(8)

又因平移轴为水平和竖直方向，故令${\mathit{\boldsymbol{T}}_0} = {\left[ {X\ Y\ 0} \right]^{\rm T}}$，分别令§1.1中所得A、R_l、T_l为：

$ \left\{ \begin{array}{l} \mathit{\boldsymbol{A = }}\left[ \begin{array}{l} {a_{00}}\ 0\ {a_{02}}\\ 0\ \ \ {a_{11}}\ {a_{12}}\\ 0\ \ \ \ 0\ \ \ \ \ 1 \end{array} \right]\\ {\mathit{\boldsymbol{R}}_l} = \left[ \begin{array}{l} {r_{00}}\ {r_{01}}\ {r_{02}}\\ {r_{10}}\ {r_{11}}\ {r_{12}}\\ {r_{20}}\ {r_{21}}\ {r_{022}} \end{array} \right]\\ {\mathit{\boldsymbol{T}}_l} = \left[ \begin{array}{l} {t_x}\\ {t_y}\\ {t_z} \end{array} \right] \end{array} \right. $

(9)

则式(3)可写为如下形式，以此表示物空间坐标系到右像空间坐标系之间的位姿关系：

$ \left[ \begin{array}{l} x\\ y\\ 1 \end{array} \right] = s\mathit{\boldsymbol{A}}\left( {{\mathit{\boldsymbol{R}}_r}\mathit{\boldsymbol{W}} + {T_r}} \right) $

(10)

将式(7)~(9)代入式(10)，可得:

$ \left[ \begin{array}{l} x\\ y\\ 1 \end{array} \right] = s\mathit{\boldsymbol{A}}\left[ \begin{array}{l} X(\cos \varphi {r_{00}} + \sin \varphi {r_{20}}) + Y(\cos \varphi {r_{01}} + \sin \varphi {r_{21}}) + \cos \varphi {t_x} + \sin \varphi {t_z} + x\\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ X{r_{10}} + Y{r_{11}} + {t_y} + y\\ X(\cos \varphi {r_{20}} - \sin \varphi {r_{00}}) + Y(\cos \varphi {r_{21}} - \sin \varphi {r_{01}}) + \cos \varphi {t_z} - \sin \varphi {t_x} \end{array} \right] $

(11)

令：

则式(11)可写为：

$ \left[ \begin{array}{l} x\\ y \end{array} \right] = \left[ \begin{array}{l} \frac{{{a_{00}}{t_1} + {a_{02}}{t_3}}}{{{t_3}}}\\ \frac{{{a_{11}}{t_2} + {a_{12}}{t_3}}}{{{t_3}}} \end{array} \right] $

(12)

至此，已根据求得的物空间坐标系到右像空间坐标系之间的旋转平移量的初值得到了右像片上棋盘格各角点的像素坐标模拟值。之后将其与同名像点的观测像素坐标值比较，利用最小二乘迭代更新左右像片间的旋转平移量，即为两位置间的位姿关系。本文以一个棋盘格角点为例进行最小二乘迭代的计算过程，误差方程为：

$ \left[ \begin{array}{l} \Delta {x_1}\\ \Delta {y_1} \end{array} \right] = \left[ \begin{array}{l} {x_1}\\ {y_1} \end{array} \right] - \left[ \begin{array}{l} {x_1}^0\\ {y_1}^0 \end{array} \right] $

(13)

式中，(x₁，y₁)为经过上述运算得到的棋盘格角点的拟合像素坐标；(x₁⁰，y₁⁰)为直接从像片上读取的同名像点的观测像素坐标，将其作为真值。

将式(12)代入式(13)得：

$ \left[ \begin{array}{l} \Delta {x_1}\\ \Delta {y_1} \end{array} \right] = \left[ \begin{array}{l} \frac{{{a_{00}}{t_1} + {a_{02}}{t_3}}}{{{t_3}}}\\ \frac{{{a_{11}}{t_2} + {a_{12}}{t_3}}}{{{t_3}}} \end{array} \right] - \left[ \begin{array}{l} {x_1}^0\\ {y_1}^0 \end{array} \right] $

(14)

式(14)中待求未知数为φ、x、y，对式(14)进行泰勒展开，并将角点数量扩增为n个角点，写成矩阵形式为：

$ \left[ \begin{array}{l} \frac{{\partial \Delta {x_1}}}{{\partial \varphi }}\ \frac{{\partial \Delta {x_1}}}{{\partial x}}\ \frac{{\partial \Delta {x_1}}}{{\partial y}}\\ \frac{{\partial \Delta {y_1}}}{{\partial \varphi }}\ \frac{{\partial \Delta {y_1}}}{{\partial x}}\ \frac{{\partial \Delta {y_1}}}{{\partial y}}\\ \ \ \ \vdots\ \ \ \ \ \ \ \ \vdots\ \ \ \ \ \ \ \ \vdots \\ \frac{{\partial \Delta {x_n}}}{{\partial \varphi }}\ \frac{{\partial \Delta {x_n}}}{{\partial x}}\ \frac{{\partial \Delta {x_n}}}{{\partial y}}\\ \frac{{\partial \Delta {y_n}}}{{\partial \varphi }}\ \frac{{\partial \Delta {y_n}}}{{\partial x}}\ \frac{{\partial \Delta {y_n}}}{{\partial y}} \end{array} \right]\left[ \begin{array}{l} \Delta \varphi \\ \Delta x\\ \Delta y \end{array} \right] + \left[ \begin{array}{l} \Delta {x_1}\\ \Delta {y_1}\\ \ \ \ \vdots \\ \Delta {x_n}\\ \Delta {y_n} \end{array} \right] = 0 $

(15)

通过式(15)可求出Δφ、Δx、Δy，将之前得到的φ、x、y加上求得的改正数再重复进行最小二乘，直到所求出的改正数小于某一限差即可停止迭代，最终求得的φ、x、y即为两位置处的左右像片之间的位姿关系。因为待求未知数的真值为控制轴运动时所输入的脉冲数，所以最终得到的所有改正数Δφ、Δx、Δy相应项之和即为两运动位置间的运动误差，分别为绕竖直轴旋转的角度、水平轴运动的距离与竖直轴运动的距离。

1.4 误差补偿模型的建立

分别对多个运动位置处拍摄像片，得到多组运动误差值，结合式(1)利用非线性最小二乘算法^[20]和全局最优化算法对误差数据进行处理，最终得到各个运动轴的误差补偿公式。具体求解过程如图 3所示。

2 实验和分析

首先用摄像头从不同角度对棋盘格标定板拍摄25张像片，求出其内方位元素A，且每张像片的平均角点误差为0.398 7个像素。

然后根据摄像头在位置1处拍摄到的标定板像片，按照§1.3中所述步骤计算出其外方位元素，对3个运动轴分别输入脉冲数，使摄像头从位置1运动到位置2，其中两个平移轴均是25 000个脉冲对应10 cm，旋转轴为6 400个脉冲对应360°。此处直接输入的运动单位即为脉冲数，不需转换为长度或角度，方便以后误差补偿模型的建立。接着根据式(5)得到位置2处所拍摄右像片的外方位元素，进而得到该像片角点的拟合值，与对应角点的观测值构成同名像点，列出式(15)所示误差方程，并通过不断迭代得到运动误差，至此即可得到一组误差值。

通过对3个轴输入不同的脉冲数并在不同时间正反向多次测量，采用上述方式共得到300组误差数据，对这些误差数据求均值整合后得到96组误差数据，对其随机排列进行非线性最小二乘拟合建模，得到各个轴的误差补偿公式如下。

旋转轴误差补偿公式为：

$ \begin{array}{c} \Delta \varphi = 1.825 \times {10^{ - 5}}{\varphi ^2} - 1.274 \times {10^{ - 9}}{x^2} + \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 1.726 \times {10^{ - 10}}{y^2} - 1.201 \times {10^{ - 3}}\varphi - 3.290 \times \\ \ \ \ \ \ \ \ \ \ \ \ \ {10^{ - 6}}x - 8.783 \times {10^{ - 7}}y + 1.506 \times {10^{ - 4}} \end{array} $

(16)

式中，Δφ为角度误差值，单位为rad。利用模型求得其拟合均方根误差(root mean square error，RMSE)为0.001 6，相关系数R为0.982 4。

水平轴误差补偿公式为：

$ \begin{array}{c} \Delta x = - 9.881 \times {10^{ - 4}}{\varphi ^2} + 2.813 \times {10^{ - 6}}{x^2} + \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 3.705 \times {10^{ - 6}}{y^2} + 1.490 \times {10^{ - 1}}\varphi + 2.727 \times \\ \ \ \ \ \ \ \ \ \ \ \ \ {10^{ - 4}}x - 1.134 \times {10^{ - 2}}y - 5.363 \times {10^{ - 2}} \end{array} $

(17)

式中，Δx即为水平轴运动误差，单位为mm。利用模型求得其RMSE为0.500 0，相关系数R为0.996 0。

竖直轴补偿公式为：

$ \begin{array}{c} \Delta y = 5.363 \times {10^{ - 5}}{\varphi ^2} - 1.862 \times {10^{ - 6}}{x^2} - \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 2.087 \times {10^{ - 6}}{y^2} + 1.229 \times {10^{ - 1}}\varphi - 1.353 \times \\ \ \ \ \ \ \ \ \ \ \ \ \ {10^{ - 2}}x + 7.721 \times {10^{ - 4}}y + 1.135 \times {10^{ - 1}} \end{array} $

(18)

式中，Δy即为竖直轴运动误差，单位为mm。利用模型求得其RMSE为0.551 7，相关系数R为0.997 1。

需要注意的是，以上误差补偿公式中，φ为其绕竖直轴的旋转角度，单位为控制器中输入的脉冲数；x和y分别为水平方向和竖直方向平移量，单位也为输入的脉冲数。

图 4(a)、4(c)、4(e)分别列出了旋转轴、水平轴和竖直轴的模型拟合图，其中蓝线为对应轴所测的误差数据，红线为拟合模型。图 4(b)、4(d)、4(f)分别表示旋转轴、水平轴和竖直轴模型拟合残差图。通过各轴的模型拟合图可知，该模型可以对各轴所测误差数据进行较好的拟合。分析其对应的模型残差图也可得，经该模型补偿后各轴的误差绝对值分别为：旋转轴小于0.003 rad，水平轴和竖直轴小于1 mm。

另从建模数据源外取10组误差数据进行补偿，其结果见表 1。表 1中Δφ'、Δx'、Δy'分别表示旋转轴、水平轴和竖直轴运动误差的真值，分析表 1可得，利用该模型进行误差补偿，补偿后旋转轴、水平轴、竖直轴的运动误差至少可以补偿80%。

式(16)~(18)即为利用观测-拟合的方法得到的各运动轴误差补偿公式，上述实验结果表明：

1) 利用同名像点位姿变化的方法得到运动误差补偿模型，绕轴旋转的角度对平动轴误差影响较大；相反地，平动轴对旋转角度误差影响较小，两平动轴之间的影响也较小。

2) 实验数据输入值均为控制器中所输入的脉冲数，因为较小的脉冲数对应的角度或距离较小，导致所建模型未知数系数较小，可通过改变输入的数据单位提高精度。

3) 本文提出的方法建立的误差补偿模型操作方便，且能很好地减小误差，经模型补偿后精度可达亚毫米级。

3 结语

考虑到系统自身的量测属性和已有条件，本文从相机标定原理出发，提出了一种利用同名像点位姿变化建立某型器物量测系统误差补偿新方法。与传统的误差检测方法相比较，其无需专业的量测仪器和复杂的操作，效率高，成本低，同时利用二阶多项式建立误差补偿模型，模型参数较少，回避了传统建模中的误差项辨识，简化了建模流程。结果表明，本文方法能对系统误差进行有效的补偿，同时也为类似系统的误差补偿提供了一种新的思路，具有广泛的应用价值。