多相机网络联合约束优化的高精度三维变形全场测量

张贵阳，霍 炬，杨 明，周 婞，朱子健

（1. 哈尔滨工业大学 航天学院，黑龙江 哈尔滨 150001；2. 哈尔滨工业大学 电气工程及自动化学院，黑龙江 哈尔滨 150001；3. 广州东软医疗有限公司CT 影像物理与重建研发部，广东 广州 510000）

1 引 言

随着计算机技术、图像处理和传感器成像技术的不断发展，视觉测量技术的应用遍及航天器对接、大型飞机制造、医学图像分析等不同领域［1-3］。而由两台或多台相机构成的立体视觉测量系统，因其精度高、测量范围大、重复性测量结果误差稳定而获得了更为广泛的关注。

近年来，基于视觉测量模式的物体三维形貌及其变形信息的测量越来越受到科研人员的高度重视，三维变形信息的测量对试件力学性能评估、结构优化设计以及健康状态检测等方面提供重要的数据依据。数字图像相关法［4-5］是当前实现变形测量应用最为广泛的方法，立体变形测量可以借助数字图像相关方法配准不同相机采集的物体表面图像散斑，然后根据相机标定得到的相机内外参数，从而重构得到空间坐标，并最后计算出三维位移场和应变场［6-7］。近年来清华大学金观昌等［8］、北京航空航天大学潘兵等［9］、东南大学何小元等［10］以及中国科技大学张青川等［3］基于数字图像相关方法对试件变形测量进行了非常系统的深入研究，并取得了令人瞩目的研究成果。

然而在实际测量需求中发现，对于航空器机翼、机械轴承以及试验弹体的应变及形貌测量中，这些被测对象具有的共同特点是长宽比例较大。若只考虑图像宽度方向尺寸，则高度方向的有效分辨率降低；而若只考虑图像高度方向尺寸，则宽度方向尺寸会超出相机的视场范围，不能有效地进行全场测量。为了提高成像的有效分辨率，杨峰等［11］通过拼接的方式获取目标全场图像。由于配准的精度有限，拼接后的全场图像的质量会决定最终的大尺度测量精度。彭聪等［12］基于视频相位的二维振动测量方法，并结合双目立体视觉，设计了一种基于机器视觉的三维振动测量方法，但是振动状态下相机标定必然会存在参数扰动。项大林等［13］所用的方法由于视差的存在，导致试件在左右相机中成像会产生一定的畸变，仅经过相关运算和亚像素定位获得的右参考图像上的对应投影点会存在较大偏差。文档［14］通过控制相邻相机之间的视场成像重叠区域，并利用传统的极线几何进行特征点立体配准，计算出的结果会具有折痕且精度不高等局限性。Dizaji 等［15］通过寻找子区相似度的整个区域较大，在立体配准过程中对目标整幅图像进行全局搜索，消耗大量的运算时间，同时配准结果的稳定性得不到有效保证。刘聪等［16］根据特征点检测及配准算法定位特征点对的亚像素位置，然后通过数字图像相关方法对特征点对进行高精度配准，最后利用发展的逐步优化单应矩阵方法求解图像变换关系，得到变形前后的无缝拼接图像，但是依然没有考虑立体相机间的固定约束关系。单宝华等［17］提出一种基于极线约束修正数字图像相关配准的立体视觉测量方法，应用极线约束对配准计算结果进行修正，取位于右极线上距数字图像相关配准结果最近的点作为新的配准点，并将该种配准修正方法推广到右图像时序配准，但是基础矩阵的求取精度决定着修正的准确性。此外，虽有一些方法［18-19］考虑到同一个空间点在立体相机上的成像满足一定的对应关系，然而对应关系矩阵的求取存在较大的偏差，这样也不能保证立体配准达到一个较高的精度。

另一方面，多相机间相对位姿的高精度求取是实现散斑精确立体配准的关键因素。王安然等［20］先获得小标定板的相对位置关系，然后将小标定板统一到同一坐标系下，最后得到多相机之间的外参数，无需考虑相机间视场是否重叠。张来刚等［21］，将空间控制点的三维坐标与其对应像素点坐标之间的共线方程作为约束条件，建立测量视场中控制点与像点坐标间的约束方程，实现了系统参数的在线校准。文献［22］则致力于双目相机系统中三维重建特征点坐标与真值之间的误差距离，构建包含多参数的最小化目标方程，最后通过Levenberg-Marquardt 迭代算法寻找最优解。张超等［23］提出基于消失点之间互相约束的多摄像机标定方法，通过靶标特征点约束结合摄像机的射影不变性排除杂点的干扰，再利用空间消失点之间夹角一致性，以及靶标特征点所构成的直线和消失点的反向射线平行性来求解相机参数。为了进一步提高相机的标定精度，误差补偿模型［24］被建立并用于减小系统标定误差，通过拟合估计模型位置与实际位置之间的对应关系来补偿偏心误差，并取得了很好的效果。文献［25］则是从标定物的灵活度入手，利用一维标定物可自由运动，不受场地环境约束的优势，实现了基于归一化算法的多相机分层逐步标定法。以上用于多相机网络标定的方法在操作步骤、计算量上存在一定的繁琐性，同时没有从全局的角度去构建多相机参数的联合优化方式。

基于上述分析，本文针对弹体、机翼等几何尺寸比例细长型被测对象，通过建立多相机网络联合约束优化把多个相机绑定为一个相机，并基于数字图像相关法分区域实现细长型或大尺寸试件的立体变形测量。将相机间的固定约束关系引入到相关函数运算中，有效缩减亚像素搜索范围，实现散斑的高精度快速立体配准。同时，多相机间联合约束关系的引入优化了试件与多相机网络之间的位置及姿态参数，进一步提高了散斑点空间三维坐标的重构精度，从而提高变形量的测量精度。

2 三维变形测量原理

2.1 数字图像相关方法

数字图像相关方法（Digital Image Correlation Method，DICM）是通过处理被测试件变形前的参考图像和变形后的目标图像来获取位移和变形量信息。在参考图像中选取以待测点为中心的配准区域，在使用数字图像相关方法进行实际计算时，通常将参考图像中间的待计算区域划分成虚拟网格的形式，通过依次计算每个网格节点的位移从而得到全场的位移信息。

设被测目标变形前后的子区对应关系如图1所示，其中参考子区中心点为P，其像素坐标为(x0，y0)，配 准 区 域 包 含 像 素 数 为(2M＋1)×(2M＋1)，Q(x0＋Δx，y0＋Δy)点为参考子区内任意一点，变形后参考子区在被测物体表面的位置发生移动，形状发生改变，Q点则移动到Q'(x'，y')点，P点变为P'点。将这个子区内任意一点在变形前和变形后的坐标建立关系式，可得：

图1 子区相关运算示意图Fig.1 Schematic diagram of subset correlation operations

对于非均匀应变，还可采用二阶或更高阶形函数来描述变形前后的坐标关系：

2.2 立体配准问题分析

通常情况下，只需使用一阶均匀应变来描述常规材质试件的变形量，从式（1）可以看出，一阶形 函 数 包 含 六 个 分 量：u，ux，uy，v，vx，vy，分 别 表示水平和竖直方向上的位移量以及对应的应变系数。相关函数是评价参考图像子区与变形图像子区相似程度的函数。零均值归一化最小平方距离相关函数具有很好的抗噪特性和适应性，因此本文采用该相关函数进行散斑的二维配准计算：

然而上述的数字图像相关测量方法需要满足被测物体为一个平面或者近似为一个平面、物体表面与光轴垂直且靠近视场中央、加载过程中的物体离面位移分量需非常小等条件，而这些必要条件限制了数字图像相关法在实际三维变形测量中的直接运用。

而在三维变形测量中，如图2 所示，除了使用相关函数将变形前后的图像进行时序配准之外，还需将选用的不同相机采集的图像对进行立体配准，从而得到被测点的空间三维坐标信息。如果依旧使用相关函数作为左右相机两个子区相似度的评判标准，找到最大相关系数时对应的图像点，那么在立体配准过程中对目标图像进行全图搜索，将会消耗大量的计算时间。同时由于视差的存在，导致试件在左右相机中成像会产生一定的畸变，那么经过相关运算和亚像素定位获得的右参考图像上的对应投影点会存在较大偏差，配准结果的偏差将会严重影响后续的三维重建效果。

图2 多相机三维立体变形测量示意图Fig.2 Schematic diagram of three-dimensional deformation measurement with multi-cameras

3 算法设计

3.1 多相机网络联合约束优化

影响基于数字图像相关法位移和应变测量结果的两个关键性因素：一个是相机网络内外参数精度，另一个是立体视觉的配准方法。即使是最好的配准算法也不能克服因立体相机间相对位姿不准确引起的测量误差，因此下面将进行多相机间联合约束关系的优化。

设定多相机网络中共有η＋1 个相机组成，其中一个为主相机ξ，其余η个为从相机。令空间中一点P的三维坐标为(Xw，Yw，Zw)。根据透视投影原理可以得出主相机的共线性方程：

其中：λξ是主相机投影尺度系数；(uξ，vξ)为主相机采集图片中特征点的像素坐标；[Xw Yw Zw]T是空间三维点世界坐标；Mξ为主相机的内部参数；(Rξ，tξ)为主相机的外参矩阵。

通过式（4）分别对相机外参和空间点三维坐标分量求一阶偏导，得到共线性误差方程：

矩阵γξ表示图像上特征点的实际坐标与用共 线 性 方 程 计 算 得 到 的 重 投 影 坐 标(uˆξ，vˆξ)之 间的差值：

多相机网络联合约束优化目的是充分利用从相机与主相机之间的固定约束关系，因此需要将主从相机间的位姿关系进行统一表示。我们通过传统标定法获得主从相机间的相对外参作为联合约束迭代运算的初值，记从相机与主相机之间的固定几何约束关系分别表示为(Rξ1，tξ1)，(Rξ2，tξ2)，…，(Rξη，tξη)，则 从 相 机 的 外 参 数(Rη，tη)可通过（8）式进行转换：

因此，根据主从相机间的相对外参关系，第η个相机的共线性方程可以表示为：

式（9）中：λη是从相机投影尺度系数；(uη，vη)为第η个相机下图像特征点像素坐标；Mη为第η个相机的内参矩阵。

结合式（8）和（9）可以得到第η个从相机的共线性误差方程：

式（10）中参数含义可以与式（5）进行类比。由式（10）可以看出，只需要对主相机外部参数进行求导并且只用计算主相机外参和三维点的改正数即可。这极大程度上减少了对所有相机都进行运算的时间与运算量。

接下来是对法化矩阵进行推导，首先Jacobi矩阵J可以表示为：

从（13）式可以看出多相机联合约束优化算法在处理时，方程总个数不变，但是待优化的参数减少了。并且随着相机数量增多，此算法的优势愈明显。为了直观体现多相机网络联合优化步骤，其流程图如图3 所示。

图3 多相机网络联合优化流程图Fig. 3 Flowchart of multi-camera network joint optimization

由于多相机间的外部参数直接影响变形测量的精度，基于多相机网络联合约束优化的数字图像相关方法能够实现相机间位置及姿态参数的优化，提高了散斑点的空间坐标重构精度。

3.2 散斑立体配准策略

基于数字图像相关法的三维变形测量中，需要将左右相机采集的散斑图像进行立体配准，通过二维像素信息还原被测目标表面各点的三维空间坐标信息，以此来完成目标表面三维形貌及变形信息的全场测量。图像对间相配准的散斑满足极线约束关系，而2.1 节则精确地获得相机间的外参关系，因此可以进行相机网络约束条件与相关函数的结合，这样立体配准中仅需在极线附近的阈值范围内进行对应点搜索，能够节省散斑立体配准的时间。

将通过联合约束优化得到的立体相机间外参的数学关系式代入到相关函数中，设极线方程为：

多相机网络外参相对关系的确定，使得不同相机成像的图像对之间同名散斑点的立体配准搜索区域被限制在极线约束的邻近范围内。因此所组成的双目相机系统中，右相机成像面的散斑搜索范围可表示为：

式（15）中，emin和emax分别表示在右相机图像中在极线附近搜索对应散斑点亚像素的下限和上限。

我们通过[emin，emax]设定了在极线方程的上下界内浮动的搜索区域，寻求最优的相关系数峰值。则新生成的相关函数可改写成式（16）关于x'i和e的 约 束 方 程：

使用加入极线约束后的相关函数进行亚像素搜索能够简化在后续的迭代算法的过程中的Jacobi 矩阵及Hessian 矩阵的求取。基于多相机网络联合约束优化的数字图像相关三维变形测量方法，加快了散斑点立体配准速度的同时，提高了变形测量精度并扩大了有效的测量范围。

3.3 三维变形测量实现步骤

设空间中的散斑点P(Xw，Yw，Zw)在选定的左右两台相机中畸变矫正后的投影点坐标分别为(ul，vl)和(ur，vr)，选取左相机坐标系为世界坐标 系，u0l，v0l，αl，βl为 左 相 机 的 内 参 数，[r1，r2，...，r9]和[tx，ty，tz]分别是构成双目相机间的外参数。在完成左右相机图像对中同名散斑的立体配准后，就可以通过下式求取散斑在世界坐标系下的立体坐标。

基于多相机网络联合约束优化的高精度三维变形全场测量流程如图4 所示，具体实施步骤如下：

图4 基于多相机网络联合约束优化的三维变形测量流程Fig. 4 Schematic of three-dimensional deformation measurement based on multi-camera network joint constraint optimization

Step 1：搭建视觉测量系统硬件架构，通过解算立体相机间联合约束优化方程，精确获得相机网络的内外参数，内参数包括：焦距、主点、畸变系数等；外参数包括两个相机之间的旋转矩阵和平移向量。

Step 2：对附有散斑的被测试件进行图像采集，对各个相机从不同角度采集的图像进行编号。由于受到噪声、环境干扰等因素的影响，需要对图像进行预处理。

Step 3：利用步骤1 中多相机间的内外优化参数，结合式（16）完成选定的左右相机下的散斑图像子区立体配准。

Step 4：进行散斑图像的时序配准，通过选定的左相机获取试件变形前的参考图像与变形后的目标图像，进行时序配准运算。然后再次用步骤3 完成变形后左右图像中同名散斑点配准。

Step 5：对运算子区散斑点坐标进行畸变矫正，然后根据立体相机间的外参数，采用式（17）进行散斑的空间三维坐标重构。

Step 6：在物体表面建立坐标系，然后对变形场数据进行拟合，进而完成变形量的计算。

4 实验结果与分析

4.1 位移测量实验

测量系统由4 台4M140MCX 数字相机和采集传输设备组成。 像幅尺寸设定为2 048×2 048 pixels，图像的采样帧频设为30 Hz，相机镜头型号均为CHIOPT HC1605A 12.5 mm/1∶1.4。多相机采集系统硬件架构如图5 所示。

图5 中的具有99 个圆的平面标定板用于校准多相机间的联合约束关系，平面标定板的圆心间距为10 mm，加工精度是5 μm，能够满足使用要求。

在进行散斑图像采集之前，先对多相机间的联合约束关系进行优化求解。为了方便求解，考虑到相机型号相同，因此4 个相机的内参在一定程度上是相同的。基于如图5 所示的相机分布方式，对视场内不同位姿下的平面标定板进行图像采集，然后提取圆心、编码，得到基于多相机联合约束优化算法的多相机参数如表1 所示：

图5 测量系统相机架构Fig. 5 Camera architecture in measurement system

表1 相机内参以及相对位姿参数Tab.1 Camera internal parameters and fixed pose relationship

通过表1 中的相机参数可以将多相机绑定在一起，重构出的多相机网络空间位置和姿态相对关系如图6 所示。

图6 多相机网络空间分布重构Fig.6 Multi-camera network space scene reconstruction

在完成多相机网络间的相对参数优化后，就可以采用该多相机架构进行目标的位移测量实验。将附有人工散斑图案的测试铝板固定在三坐标测量机（Hexagon Global SR）的测头上，如图7 所示。三坐标机测量范围为900 mm×2000 mm×800 mm，测量分辨率为0.2 μm，测量精度为0.8＋L/400 μm（L为测量距离）。

图7 三坐标测量机实验装置Fig. 7 The CMM experimental device

三坐标机测头能够在X-Y-Z三个方向分别移动。先通过摄像机采集运动之前的图像作为基准图像，然后让铝板依次在X，Y和Z方向上每隔5 mm 做一次点到点的运动，并进行图像采集，等间隔选取150 个测试点，子区尺寸为41×41，步长为15 pixel，利用本文方法进行位移计算与分析，并将计算结果与三坐标机给定的真值数据进行比对，得到的测量误差归一化统计直方图如图8 所示。

从图8 可以看出X与Y方向的误差范围均在在0.036 mm 以内，同时根据直方图发现误差较为集中的分布在0.022 mm～0.030 mm 范围内。对于Z方向面外位移最大误差是0.044 mm，最小误差是0.014 mm，深度方向位移测量误差比面内位移测量误差大。

图8 位移测量误差统计直方图Fig. 8 Statistical histogram of displacement measurement error

接着分析本文算法在不同数量散斑点下的运行时间。为了方便叙述，基于本文联合约束关系优化的配准方法记为RJCO（Registration with Joint Constraint Optimization）算 法，并 与IGGA（Iterative and Gray-Gradient Algorithm）［26］，Normal-SRPG（Normal Stereo Registration Based on Polar Geometry）［14］和IC-GN（Inverse Compositional Gauss–Newton）［27］三 种 具 有 代 表 性 的 方法进行执行效率的对比，经过MATLAB 多次计时的平均耗时结果如图9 所示。

从图9 中可以看出，随着被测散斑点数量的增加，这四种算法的运行时间均呈增大趋势。虽然MATLAB 的tic 和toc 组合计时指令具有一定的误差，但是在散斑配准条件一致的情况下，从耗时的相对性来看，本文RCJO 算法立体配准速度要比Normal-SRPG 与IC-GN 算法快；虽然IGGA 算法耗时较少，但是设置迭代精度较低，不然难以收敛。

图9 同名散斑点立体配准时间Fig. 9 Homologue speckle registration time

由此可见本文算法在保证测量精度的基础上具有效率上的优势，这是因为联合约束关系的优化将图像对间立体配准点的搜索限制在对极线附件，而不是整幅图像上，由此提高了搜索效率。

4.2 变形测量实验

位移测量实验证明了本文方法能够对成功匹配的散斑点进行高精度三维重建，接着进行试件三维变形测量实验。将喷有散斑的测试铝板紧固在三坐标测量机固定位置上，首先采集变形前参考图像作为基准，然后将外力作用于铝板进行不同程度的弯曲变形，采集的图像如图10（a）所示。对采集的图像进行ROI 区域的选定和测量。通过本文方法测量的结果如图10（b）和（c）所示。

图10 计算区域三维形貌重构Fig. 10 Reconstruction of 3D shape of calculated region

从图10（b）可以直观地看出测量结果与试件的实际三维变形是一致的。图10（c）显示试件变形较为均匀，而且通过本文方法测量结果拟合较为平滑，结果较为稳定。此外，图10（b）中三个变形图离面方向的最大位移量分别是3.690 mm，6.985 mm 与10.260 mm，与三坐标机校准的结果是一致的。

为了体现本文方法的变形测量精度，与另外三种方法测量结果进行定量对比。其中，IGGA和IC-GN 算法的三维重建参数获取均采用Normal-SRPG 方法中的立体相机标定方式。用三坐标测量机在如图10（b）中蓝色直线上进行等间距采样点坐标测量并记录，通过相对应的分量相减，获得离面（Z向）变形量曲线，如图11 所示。同时，通过以上四种视觉变形测量方法对蓝色直线上的对应点进行空间三维坐标重建，并绘制出离面位移量曲线，如图12（a）所示；与三坐标测量机提供的真值之间的误差绝对值曲线如图12（b）所示。

图11 三坐标机校准曲线Fig. 11 Deformation curve by CMM

图12 蓝色标识线上点的离面变形测量结果Fig. 12 Measurement of the off-surface deformation of the points on blue marking line

总体来看，图12 表明相机间联合约束关系的优化很好地修正了散斑图像的立体配准，能够将三维变形测量误差控制在0.06 mm 范围内，有效提高了三维变形的测量精度和稳定性。对于Normal-SRPG 算法，其测量结果很稳定，但由于基础矩阵求解的精度不高，导致图像间的极线方程有所偏差，从而影响最终测量精度。相比之下，IC-GN 算法虽然在二维图像相关求解时精度很高，但是直接用于图像对间的同名散斑立体配准时，会受视差的影响导致匹配精度下降。如果在相机光轴角度较大情况下采集的图像对间的同名散斑点进行立体配准，那么IC-GN 和IGGA算法容易会出现去相关现象。而本文的RJCO算法对选定左右相机部分成像的区域进行配准，因此减小了图像畸变的影响，同时也提高了图像的有效分辨率。此外，在这种情况下，相机之间位置关系得到的基础矩阵元素所确定的极线斜率较小，减小了竖直轴的搜索方向范围。综上可知，RJCO 算法的总体性能优于其它三种算法。

5 结 论

本文提出一种基于多相机网络联合约束优化的三维变形实用测量方法。每对相机仅获取试件的部分图像，从而提高了图像的有效分辨率。通过相机间约束关系的优化求解，使图像对间的立体配准点搜索范围限于极线附近，以使得散斑图像对中的同名点可以被准确且快速地搜索。然后利用立体视觉视差原理重构求出变形前后同名空间点的坐标，进而可以进行变形量的拟合。此外，通过优化试件与多相机网络之间的位置及姿态参数，进一步提高了散斑点空间三维坐标的重构精度。最后，实际测量实验验证了本文所提方法的有效性。

下一步工作将考虑散斑图像的对比度对三维变形测量的影响，同时研究适用于环境光照变化下的立体配准策略，以便实现能够在场外环境下进行立体变形测量的视觉应用系统。