大视场远距离视觉测量系统的分步标定

胡 浩，魏 斌，梁 晋，王惠刚，张永庆

（1. 西北工业大学 航海学院 青岛研究院，陕西 西安 710072 2. 河北科技大学 机械工程学院，河北 石家庄 050091 3. 西安交通大学 机械工程学院，陕西 西安 710049）

1 引 言

视觉测量系统正被广泛用于科学研究和工业现场，主要任务是从拍摄的二维图像中解析出对应场景的三维信息，如表征物体几何形状的关键位置坐标，由此识别或重建物体［1-2］。然而受工作距离、视场范围及成像分辨之间的相互制约，视觉系统多用于自身所限定视场范围内的尺寸测量和检测。随着工业制造水平的不断提高，工业制造及军用场景对大范围、大尺寸高精度测量的需求越来越多，如大型装备空间位置的定位、风洞试验模型的运动姿态跟踪、飞机蒙皮的损伤监测、机翼飞行过程中的变形测量、重型直升机旋翼桨叶的动态分析以及水下远距离光学三维成像与探测等［3-9］。

相机标定是三维视觉系统的关键。在建立相机成像模型后，一般将通过实验与计算求解其参数的过程称之为相机标定，其本质是建立二维图像与三维空间的坐标映射关系［2］，标定的精度直接影响三维重建的准确性，尤其是对于大视场大尺寸测量。经典的标定方法一般是通过采集具有已知欧式空间结构的标靶图像来实现的。最常用的R. Y. Tsai 提出的两步法［10］和张正友提出的棋盘网格平面标定法［11］都需要制作较高精度的标靶来保证标定的最终精度。李文龙等提出的用于涡轮叶片检测的高精度标定方法标定精度可达0.02 mm，主要用于小视场测量，且需要精确测量所用三个靶球的中心位置［12］。王向军等提出的利用三坐标设备和发光二极管对相机进行分区域标定方法能够克服大型标靶制造困难和小型标靶在大视场中标定精度低的缺点［13］，但由于需要借助三坐标测量设备，难以在实际工程现场使用。Faugeras 等［14］提出的相机自标定的思想，不需要已知标靶欧式空间结构和其他的空间三维信息，仅利用多次成像之间的约束关系来计算相机参数，并将相机内参数、外参数和标靶三维信息都作为未知量进行整体光束平差解算，目前在工业测量中最为常用［15-20］。武汉大学冯文灏等研究者［21-24］提出并进一步发展了基于摄影测量的相机柔性标定方法。Luhmann 等发展了自动化自标定方法［25］。这些方法在一定条件下能够实现高精度的相机标定，但要求标靶尽可能覆盖整个视场区域，并需要制造与测量视场尺寸匹配的标靶。这对于十米及以上更大尺寸的测量，甚至几十米的超大尺寸测量是不现实的，制造和维护均存在问题［19］，如加工难度大、精度低、自身发生变形、现场移动操作繁琐等。

由此，本文提出一种大视场视觉系统的分步标定方法，并进行了实验验证。近距离采用小尺寸标靶进行相机内参数标定，远距离视场上布置全局标志点组成控制场，实现相机外参数标定，最后进行光束平差优化，从而获得相机精确的内外参数。

2 标定方法

相机成像可用小孔成像模型近似表示。如图1 所示，空间三维点P(X，Y，Z)经过镜头投影中心S(XS，YS，ZS)后，投影到像平面上的像点为p(x，y)，o 为 主 点(x0，y0)，So为 焦 距f（单 位：mm）。理想情况下，P、S和p三点共线，但由于镜头存在各种畸变，使得各像点在像平面上相对其理论位置存在偏差(Δx，Δy)（单位：mm），则相机成像的实际数学模型为：

图1 相机成像模型Fig.1 Camera imaging model

坐标系相对于相机坐标系的变换矩阵。

本文标定算法将所有可能的畸变都考虑在内，采用的畸变模型为：

其中：K1，K2，K3是径向畸变系数；B1和B2是切向畸变系数；E1和E2是像平面畸变系数；r2=x2+y2；r0是畸变曲线的二次零点。

对于十米及以上大视场的标定，传统标定方法由于需定制与测量视场尺寸匹配的标靶，导致标定过程较为复杂且精度不高。本文利用相机透视成像模型截面相似性质，提出将相机的内参数和外参数分离进行标定的新方法。整个标定流程如图2 所示，图中右侧的圆点代表每一步引入计算的数据类型。在近距离采用小幅面标靶首先对相机内参数进行标定；在被测空间布置若干个编码标志点（编码点），利用多片后方交会原理，计算得到相机的外方位参数。

图2 分步标定流程Fig.2 Flowchart of two-step calibration

2.1 前截面内参数标定

将方程（1）线性化可得如下误差方程：

其中：V为像点坐标残差；X1，X2，X3分别为内参数，外方参数以及物方点坐标的改正数；A，B，C分别为内参数，外参数和空间点坐标对应的偏导数矩阵；L为观测值，即图像点坐标。若以图像坐标作为观测值，把相机的内参数、外参数和空间点坐标作为未知数来求解，此过程称为光束平差。

2.1.1 求解标靶标志点的三维坐标及图像坐标

所用标靶的标志点三维坐标为已知，利用西安交通大学研制的XJTUDP 工业摄影测量系统测定，该系统测量精度为0.1 mm/4 m，测量过程详见文献［20］。将标靶分别在两个相机各自的视场内摆放多个不同的方位，利用相机采集标靶的序列图像。然后，对采集图像进行处理，识别并定位标志点，求出其中心的图像坐标。

2.1.2 角锥体解算外参数初值

根据标志点的三维空间坐标和识别出的标志点的二维图像坐标，采用角锥体法计算相机的外参数，即旋转矩阵R和平移矩阵T。如图3 所示，该方法利用成像数学模型中光线顶角相等的原理求解外参数。假设Q1，Q2和Q3分别是世界坐标系下的三个已知空间点，q1，q2和q3分别为Q1，Q2和Q3所对应的图像点，根据相似三角原理：

图3 角锥体法原理Fig.3 Principle of pyramid method

对于角锥体o-q1q2q3，在图像坐标系下：

世界坐标系O-XYZ下，对于角锥体o-Q1Q2Q3存在类似关系，

由 于 空 间 点Qi(i=1，2，3)在 世 界 坐 标 系O-XYZ下的坐标和相应的像点qi在像空间坐标系o-xyz中的坐标均为已知，所以oq1、oq2、oq3、q1q2、q2q3、q1q3、Q1Q2、Q2Q3和Q1Q3的 长 度为已知。其中，oQ1、oQ2和oQ3三个长度是未知数，式（7）为非线性方程，采用最小二乘迭代法求解。求得三边长度oQ1、oQ2和oQ3之后，可得三个空间点在像空间坐标系下的对应坐标。

则旋转矩阵R为：

其中：Qw和Qc分别为三个空间点在世界坐标系和像空间坐标系下的重心。

则平移矩阵T为：

2.1.3 内参数精确求解

得到标志点三维坐标和相机外参数后，对应的误差方程（4）可简化为：

以像点坐标残差优化目标，将相机内参数的初值设为其出厂的参数，然后代入方程（12）采用最小二乘算法解算得到相机内参数X1。

2.2 后截面外参标定

在前截面通过小尺寸标靶分别标定得到每个相机的内参数后，在后截面则通过布置若干编码标志点来标定相机的外参数。

利用XJTUDP 工业摄影测量系统重建后截面编码标志点的三维坐标。如图4 所示，第一步，采用单相机从多个角度拍摄被测区域布置的编码标志点，得到一组图像。第二步，进行标志点检测，得到标志点亚像素级的图像坐标。第三步，选用两幅图像进行相对定向，计算出五个外参数，并重建出至少5 个编码点的三维坐标作为控制点。第四步，根据这些图像包含的控制点数目，依次循环处理剩余的图片；首先利用控制点定向该图片，然后与已经定向好的图片分别搜索公共的未重建的编码标志点，然后重建出来。所有图片都定向完毕后，利用光束平差算法同时调整内、外参数和物体点的三维坐标。最后加入比例尺和温度补偿，最终得到物体点的三维坐标。

图4 摄影测量三维重建处理流程Fig.4 Flow chart of 3D reconstruction process of the photogrammetry

本文采用图5 所示的编码标志点，白色实心圆点周围带有同心编码环带，每个点都拥有唯一的编码ID，有利于降低同一图像中多个标志点的误匹配。首先，采用Canny 算法搜索标志点中间实心圆的边缘，并利用椭圆最小二乘拟合算法来求解其中心图像坐标。然后，从中心向编码环径向作内外边界的连线，在连线上等距离采样5 次，以5 次采样的中值与标志点的灰度阀值（统计）作比较，大于阀值则本环带的编码为1，否则为0；每36°（对应于10 位的编码标志点）重复上述操作，旋转一周后得到形如“0101101111”的二进制编号。

图5 编码标志点Fig.5 Coding mark point

由于编码标志点在图像中的角度不同，得到的序列也会不同，所以先将其转换为十进制、得到一个确定数值，然后将这个序列循环向右移位，每移动一位得到一个新序列和对应的十进制数值，N 位编码会得到N 个十进制数值，取其中最小值作为编码的ID 值。

以测得的编码标志点三维坐标为控制点，对两个测量相机的位姿分别进行绝对定向求解其外参数。已知标志点的三维坐标和内参数，求解外参数的过程称为后方交会，对应的误差方程为：

得到外参数X2后，将相机内参数、外参数以及后截面上标志点的三维坐标一起代入方程（4），并利用列文伯格-马夸尔特（LM）方法［17］求解R，T的最优解，并引入比例系数，从而完成相机内外参数的整体优化标定。

最后，将标定得到的外参数转换到以相机1为基准的坐标系下，从而完成对大视场视觉测量系统相机的整体标定。

3 实验结果与分析

3.1 标定实验

为验证所述标定方法的可行性和精度，设计了如下标定实验。首先搭建如图6 所示的大视场实验环境，在长10 m、高3 m 的墙面上随机粘贴若干编码点（黑底白点，内圆直径40 mm）。构建双目视觉测量系统，硬件组成：BFLY-PGE-50H5M 相机（分辨率2 448×2 048 pixels，像元尺寸3.45 μm）2 个、Schneider 8 mm 定焦镜头2个、工作站1 台、相机触发控制箱1 台、三脚架2个、标定十字架（十字型标靶，表面粘贴直径为15 mm 的编码标志点）1 套、500 W 高亮度LED光源1 套、设备线缆若干，直径200 mm 的10 位编码点1 套。

图6 标定实验现场Fig.6 Experiment environment

3.1.1 内参标定

如图7 所示，将标靶分别在两个相机的视场内摆放8 个不同的方位，利用相机采集标定图像并识别编码点的图像坐标和编码ID 值，如图8。根据光束平差方程求解编码点三维坐标和相机内参数值。相机内参数计算结果见表1。

图7 内参标定示意图Fig.7 Schematic diagram of internal parameter calibration

图8 内参标定图像采集及编码点识别Fig.8 Image acquisition and coding point recognition for internal parameter calibration

表1 相机内参数标定结果Tab.1 Result of internal parameter calibration

3.1.2 外参标定

如图9 所示，拍摄墙面上编码点的多组图像（5～9 幅即可），并利用工业近景摄影测量系统重建出这些编码点的三维坐标，然后将两个相机对准编码点同时采集1 幅图像，基于多片后方交会原理解算得到相机的外参数，结果见表2。

表2 双目相机的外参数Tab.2 Result of external parameter calibration

图9 外参标定图像采集及编码点识别Fig.9 Image acquisition and coding point recognition for internal parameter calibration

为评估标定方法精度，根据标定得到的两相机的内外参数，将重建出的标靶上编码点的三维坐标反投影到相机的像平面上，计算其重投影误差，通过重投影误差来定量评估标定精度。

图10 所示为其中3 个姿态下标靶上编码点在左相机（相机1）像平面上的重投影误差分布，可以看出测得的编码点三维坐标在左图像上的投影点与实际图像点偏差基本一致。重投影误差小于0.08 个像素，表明本文标定方法能够在十米大视场下保持较好的标定精度。

图10 重投影误差分布图Fig.10 Error distribution of reprojection

3.2 精度验证

为进一步验证本文标定方法在大视场三维测量中的精度表现，参照德国VDI-VDE 2634 part 1 三维光学测量精度验证标准，采用图11 中所示标称长度为6 100.056 5 mm（出厂校检精度为0.003 mm）的因瓦合金高精度标尺作为长度基准，利用标定后的双目视觉测量系统对标尺进行多次测量并重建出其三维长度，然后与标尺的标称长度进行对比，根据得到的相对误差来评价该系统的三维测量精度。

图11 因瓦合金标尺Fig.11 Invar alloy standard ruler

测量时，手持标尺在测量视场内移动，尽量使运动范围充满整个测量空间，运动过程中保持标尺端两个标志点面向相机方向，重复拍摄并测量标尺两端点间的长度，然后将测量结果与其标称值比较，得到标尺长度测量的相对误差。

两相机拍摄的标尺部分位姿图像见图12，测量结果见表3。从实验结果可以看出，使用本文标定方法，在测量幅面宽度为10 m 时，标尺长度多次测量的相对误差最大值为0.007%。

图12 左右相机拍摄的标尺部分运动位姿的图像Fig.12 Images of the ruler in different positions

表3 长度测量结果Tab.3 Result of length measurement

综上，本文方法仅需采用小尺寸标靶、并在测量视场后截面布置若干编码点（10～30 个），即可快速方便的完成十米测量视场相机的准确标定。

3.3 十米高空旋翼测量

在某型直升机旋翼试验塔上开展了10 m 直径桨叶高空旋翼测量试验，旋翼高速运转时无地面效应，实验现场如图13 所示。首先利用小尺寸的标靶在室内完成相机内参数的标定，然后将双目视觉测量系统安装到旋翼的正下方，直接在旋翼表面随机布置若干编码点，并以此为标定的特征点进行相机外参数的快速标定和优化。

图13 旋翼塔测量实验现场Fig.13 Experiment site of the rotor tower

在桨叶根部沿着桨叶运行方向做1 条垂直于桨轴A 的异面直线L，旋翼运动时该直线变为L'，角度相对变化量即为总距角，如图14 所示。旋翼转速为514 r/min，当桨叶转入相机视场时，连续采集旋翼运动图像，测量得到的三维变形结果如图15 所示。当输入总距角为0°时，拟合旋翼平飞基准平面，桨轴方向即为基准面的垂线方向。异面直线L和L'则通过测得的桨叶根部三维数据拟合得到。

图14 旋翼总距角Fig.14 Total pitch angle of the rotor

图15 旋翼三维全场变形测量结果Fig.15 Result of three-dimensional deformation of the rotor

以拟合的基准平面为参考平面，机械操控输入总距角从0°、2°、4°、6°度依次变化，测得旋翼桨叶根部总距角如图16 所示，测得桨叶根部总距角与机械操控输入总距角误差在0.1°内，符合实际运动规律。本试验得到了旋翼塔上直升机桨叶的三维全场变形数据，解决了大视场下高速旋翼桨叶全场变形的三维测量难题，为分析直升机动力学及其相互耦合问题分析提供了科学数据。

图16 总距角测量结果Fig.16 Result of total pitch angle of the rotor

4 结 论

本文提出了一种大视场视觉测量系统多相机内、外参数的分步标定方法，并设计实验进行了可行性和精度验证。对于高空旋翼变形测量，近距离首先用小幅面标靶完成内参数的标定，然后在远距离被测空间布置若干个编码标志点，利用多片后方交会原理解算外参数，最后通过光束平差优化进一步得到相机准确的内外参数。实验表明，该方法标定后高速旋翼总距角的测量误差小于0.1°。该方法可实现相机内外参数的分离标定，即内参数标定在实验室内进行，外参数标定在现场实现，因此具备较好的工程价值。

后续工作将重点围绕如何去除后截面标定所需的编码点进行，即直接利用被测目标表面已有的特征（如随机斑点、纹理线、角点等），与前截面识别的标志点融合，在保证标定精度的同时进一步简化大视场视觉测量系统的标定流程。