基于双平面镜双相机组合测量全场变形的方法

针对单相机难以用于全场变形测量以及多相机在全场变形测量中复杂性等问题，提出一种利用双平面镜辅助的双相机系统来实现全场变形测量的新方法。该方法基于相机标定的坐标变换方法，即通过双平面镜实现从真实相机坐标系到虚拟相机坐标系之间的转换，进而实现从实物到虚像之间的转换。将该方法用于空心六角铝棒的热变形测量中，并与ANSYS软件的仿真结果进行比较，结果表明，该方法不但可以实现对被测件的三维重建，而且可以实现全场的变形测量，相对误差仅为0.51%，满足精度要求、效果良好。

Abstract

Results and Discussions With a hollow hexagonal aluminum bar as the research object, the proposed method is used to conduct 3D reconstruction and thermal deformation research (Fig. 4). In the 3D reconstruction experiment, the maximum reconstruction deviation between the reconstructed surface and the fitted surface is about 0.004 mm, and it is mainly concentrated in the splicing area, while the deviation in other areas is relatively stable. The stereo reconstruction angle also proves that the system has stable performance and high accuracy (Fig. 6). In the thermal deformation research experiment, the proposed method can be used to obtain the average height of the first and last frame images in the vertical direction of the B and C surfaces of the part, respectively (Table 1). It can be seen that the maximum change in the height of surface B of the part is 0.5467 mm, and that of surface C is 0.5589 mm (Fig. 7). Then, the height change of surface A of the part processed by the system are directly compared with that of surfaces B and C processed by the proposed method during the cooling process, and comparison is carried out every 30 frames from frame 0 to 1290. There are 44 frames of comparison images (Fig. 8). Finally, the finite element analysis software ANSYS is used to simulate the part, and the steady-state simulation is selected. The temperature of the part rises from 20 ℃ to 310 ℃ in simulation experiments. It can be seen that the maximum thermal deformation of the part in the vertical direction is 0.5500 mm, which indicates that the thermal deformation of the part along the z direction is 550.0 μM when the temperature difference is 290 ℃ [Fig. 9 (b)].Objective

The application of digital image correlation (DIC) technology is very extensive, and the technology has high practical value in biotechnology, civil engineering, aerospace, medical application, and other fields. With the continuous advancement of related technologies, the demand for full-field deformation measurement and 3D shape reconstruction of objects has increased accordingly. This requires the DIC technology to not only have higher measurement accuracy but also be more economical and practical, so as to make itself be applied to more fields. In recent years, many scholars have done a lot of research on the full field strain and deformation measurement by DIC, and they have made many valuable research results. Among them, the multi-camera DIC system has been proven to have high measurement accuracy, and it is feasible to achieve full field or double surface deformation measurement. However, in actual use, the system takes a long time to be built and involves a complex operation and high economic costs, and there is interference between multiple cameras. The rotation of a single camera is used to realize the full field deformation measurement of the measured object under different load conditions. In other words, the camera is continuously moved to rotate around the same measured object, and it will shoot and record in multiple different positions, so as to finally cover all the required fields of view. The system has low cost, complex operation, and low accuracy. In view of the problem that a monocular camera is difficult to be used in full-field deformation measurement and the complexity of multiple cameras in full-field deformation measurement, a new method of full-field deformation measurement by using a double-sided mirror-assisted dual-camera system is proposed.

Methods

This research adopts a coordinate transformation method based on camera calibration. In other words, the transformation from the real camera coordinate system to the virtual camera coordinate system is realized through the double-sided mirror, and then the transformation from the real object to the virtual image is realized. It is assumed that the camera coordinate system of L coincides with the world coordinate system of the binocular DIC system (Fig. 3). Therefore, the conversion relationship between the world coordinate system of the virtual binocular DIC system and that of the real binocular DIC system is equivalent to the conversion relationship between the camera coordinate systems of L′ and L. First, it is necessary to find out the positional relationship between the coordinate system o_c-x_cy_cz_c and the O₁-X₁Y₁Z₁. The rotation and translation matrices between the camera coordinate system o_c-x_cy_cz_c and O₁-X₁Y₁Z₁ can be obtained through camera calibration. Then, an intermediate coordinate system is introduced, namely, the coordinate system O₂-X₂Y₂Z₂, which is a rotating coordinate system of O₁-X₁Y₁Z₁. According to the imaging law of the plane mirror and the nature of the Euler angle, the rotation and translation matrices from the coordinate system o_v-x_vy_vz_v to the coordinate system O₂-X₂Y₂Z₂ can be obtained, and then the rotation and translation matrices from the coordinate system o_v-x_vy_vz_v to the coordinate system O₁-X₁Y₁Z₁ can be obtained. Finally, the conversion relationship between the real camera coordinate system and the virtual camera coordinate system can be calculated by synthesizing the above results.

Conclusions

A new method of full field deformation measurement using a double-sided mirror-assisted binocular DIC system is proposed. With a hollow hexagonal aluminum bar as the measuring object, the thermal deformation results of three outer surfaces are measured during the cooling process from 310 ℃ to 20 ℃, and they are compared with the simulation results of the finite element software. The results show that the change curves of the three outer surfaces of the part along the height direction basically coincide with the simulation results, and the absolute error between the calculated average thermal deformation values of A, B, and C surfaces along the height direction of the part by using the proposed method and the simulation results is 2.8 μm. The relative error is only 0.51%. It can be seen that the proposed method not only overcomes the limitation that a monocular camera cannot realize full field measurement but also discards the complexity of a multi-camera DIC system.

1　引言

数字图像相关（DIC）技术已有近四十年的发展史^［1-4］，其应用非常广泛，在生物技术、土木工程、航空航天、医学应用等领域有着很高的实用价值。随着相关技术的不断推进，对物体全场变形测量和三维形貌重建的需求也相应提升。这不仅要求DIC技术有较高的测量精度，也需要该技术更加经济实用，可以推广到更多的应用领域。

近年来，不少学者在DIC实现全场应变和变形测量方面进行了大量研究，也取得了很多有价值的研究结果^［5-12］。其中，多相机DIC系统被证明有较高的测量精度，对实现全场或双表面变形测量有较高的可行性。Orteu等^［13］将这一类技术称为多视角DIC（MV-DIC），并开发了以“主摄像机”配置为基础的MV-DIC系统。该系统包括4个同步摄像机，其中一个摄像机被定义为“主摄像机”，其余的摄像机均与主摄像机有公共视场，通过两两匹配，实现全场测量。由于摄像机之间要有公共视场，因此该系统的测量范围受限。Li等^［14］随后提出以“面对面的相机对”配置为基础的MV-DIC系统，即系统包含4个相机，每2个相机构成一个常规双目DIC系统，且2对相机处于同一高度左、右位置。该方法实现较简单，可实现应力-应变的大范围测量。但是，在以多相机搭建的MV-DIC系统中仍存在视场受限问题，为此一些学者提出使用更多的相机来扩大相机视场^［15-18］。在实际使用时，系统的搭建时间长，操作复杂，经济成本较高，多相机之间时有干扰。针对此类问题，Spera等^［19］提出使用单台相机旋转实现对被测对象在不同载荷情况下的全场变形测量，即连续移动这台相机围绕同一被测对象旋转，并在多个不同位置进行拍摄记录，最终覆盖全部所需视场。Badel等^［20］提出保持单台相机固定，让放在旋转平台上的被测样品按照一定速度旋转，由这台相机分别拍摄被测样品变形前后的图像序列，并对图像序列进行处理，从而实现对被测样品的全景变形测量。上述两种方法利用单台相机构成多个伪相机，简化DIC系统搭建流程，降低经济成本。另外，Barone等^［21］将相机、投影仪、球面镜结合，将投影仪发出的结构光通过球面镜投影到被测对象上，重建了不同被测对象的表面形状。Pan等^［22-23］提出双目DIC测量系统，利用反射镜上预先制作的散斑图案，可以检索两个反射镜的反射变换，并将重建的虚拟曲面转换为其在反射镜前面的真实位置。鉴于MV-DIC系统和单相机系统在测量精度和经济成本上存在的矛盾，有必要设计一种兼顾精度和成本的方案，即使用较少的相机实现多目DIC系统的功能。本文提出一种利用平面镜辅助双目DIC系统实现全场变形测量的新方法。该方法通过相机标定的坐标变换方法，得到镜面虚实像点之间的坐标变换关系，实现被测件的三维重构和全场变形测量，在提高DIC系统测量精度的基础上，不增加CCD相机数量。

2　镜面辅助双目DIC系统

2.1　系统配置和基本原理

搭建的镜面辅助双目DIC系统如图1所示。该系统由一个常规的双目DIC系统和两个平面反射镜M1、M2组成，其中，L、R为置于同一水平面的左、右相机。在被测试样表面人工喷上随机散斑并放置在水平实验台上，为了更好地恢复待测物的三维形貌，在试样后方放置两个互成120°夹角的平面镜。散斑图A、B、C中的框线区域分别表示试样的前表面、左后侧面以及右后侧面，可以被左、右相机同时拍摄到，如图1（a）、（b）所示。在试样上选取一点 $P_{M} (x_{M}, y_{M}, z_{M})$ ，在左、右相机视图中可以分别找到该点的对应点 $P_{1} (x_{1}, y_{1})$ 和 $P_{2} (x_{2}, y_{2})$ ，左、右相机视图可以构成一组立体图像对。

图 1. 镜面辅助双目DIC系统的原理图。（a）左相机视图；（b）右相机视图

Fig. 1. Schematic of mirror-assisted binocular DIC system.(a) Left camera view; (b) right camera view

下载图片查看所有图片

2.2　基于双目DIC的三维重建

在双目DIC系统中，两个相机从不同位置拍摄到被测物体的同一公共区域，类似于人眼的成像原理，其数学模型如图2所示。点 $P_{N} (X_{N}, Y_{N}, Z_{N})$ 为被测物体表面的一点，可以同时被两个相机捕获，并分别在两个相机的图像平面上成像，对应的像点分别为 $p_{l} (x_{l}, y_{l})$ 和 $p_{r} (x_{r}, y_{r})$ 。坐标系 $o_{c} ‑ x_{c} y_{c} z_{c}$ 、 $o_{c r} ‑ x_{c r} y_{c r} z_{c r}$ 、 $O_{W} ‑ X_{W} Y_{W} Z_{W}$ 分别为L、R的相机坐标系以及双目DIC系统的世界坐标系。为计算方便，令L的相机坐标系与双目DIC系统的世界坐标系重合，坐标系 $o_{1} ‑ x_{1} y_{1}$ 和 $o_{2} ‑ x_{2} y_{2}$ 分别为L、R的图像坐标系，有效焦距分别为 $f_{1}$ 、 $f_{2}$ 。

图 2. 双目DIC系统的数学模型

Fig. 2. Mathematical model of binocular DIC system

下载图片查看所有图片

采用棋盘格标定板^［24］可以方便地确定每个相机的内部参数，以及R相对于世界坐标系的旋转矩阵和平移矩阵。由图2可知，双目DIC系统的数学模型和相机的成像原理为

s_{l} [\begin{matrix} x_{l} \\ y_{l} \\ 1 \end{matrix}] = [\begin{matrix} f_{l x} & 0 & u_{l 0} \\ 0 & f_{l y} & v_{l 0} \\ 0 & 0 & 1 \end{matrix}] [\begin{matrix} X_{N} \\ Y_{N} \\ Z_{N} \end{matrix}]

，（1）

s_{r} [\begin{matrix} x_{r} \\ y_{r} \\ 1 \end{matrix}] = [\begin{matrix} f_{r x} & 0 & u_{r 0} \\ 0 & f_{r y} & v_{r 0} \\ 0 & 0 & 1 \end{matrix}] [\begin{matrix} X_{r} \\ Y_{r} \\ Z_{r} \end{matrix}]

，（2）

式中： $s_{l}$ 、 $s_{r}$ 为比例尺度因子； $f_{l x}$ 、 $f_{l y}$ 分别为L相机在x和y方向上的焦距； $f_{r x}$ 、 $f_{r y}$ 分别为R相机在x和y方向上的焦距； $u_{l 0}$ 、 $v_{l 0}$ 分别为L相机图像坐标系原点在其像素坐标系下的横、纵坐标值； $u_{r 0}$ 、 $v_{r 0}$ 分别为R相机图像坐标系原点在其像素坐标系下的横、纵坐标值。

对于一般的工业相机，二阶畸变校正即可满足对图像的校正。然而，所提方法的目的之一是实现双相机的标定，每个视角下的图像质量对系统的测量精度都有影响，有必要考虑更高阶的畸变系数。因此，所提方法采用非线性成像模型，综合考虑畸变校正效率，将径向畸变系数提高为三阶。

径向畸变校正模型为

\{\begin{matrix} x_{d} = x (1 + k_{1} r^{2} + k_{2} r^{4} + k_{3} r^{6}) \\ y_{d} = y (1 + k_{1} r^{2} + k_{2} r^{4} + k_{3} r^{6}) \end{matrix}

，（3）

式中： $(x_{d}, y_{d})$ 为畸变图像点； $(x, y)$ 为校正过的理想图像点； $r^{2} = x^{2} + y^{2}$ ； $k_{1}$ 、 $k_{2}$ 、 $k_{3}$ 为径向畸变系数。

切向畸变校正模型为

\{\begin{matrix} x_{d} = x + 2 p_{1} x y + p_{2} (r^{2} + 2 x^{2}) \\ y_{d} = y + 2 p_{2} x y + p_{1} (r^{2} + 2 y^{2}) \end{matrix}

，（4）

式中： $p_{1} 、 p_{2}$ 为切向畸变系数。

一般工业相机采用一阶校正系数 $k_{1}$ 即可实现对图像的校正，二阶校正系数 $k_{2}$ 主要用于对图像边缘较大畸变的校正，而三阶校正系数 $k_{3}$ 更适用于对成像质量差、图像严重畸变情况的校正。在双相机精密测量中，高阶畸变模型更有利于避免随机误差的产生。

对于L、R的相机坐标系，有如下关系：

[\begin{matrix} X_{r} \\ Y_{r} \\ Z_{r} \end{matrix}] = H_{l r} [\begin{matrix} X_{N} \\ Y_{N} \\ Z_{N} \\ 1 \end{matrix}] = [R | t] [\begin{matrix} X_{N} \\ Y_{N} \\ Z_{N} \\ 1 \end{matrix}]

，（5）

式中： $R = [\begin{matrix} r_{1} & r_{2} & r_{3} \\ r_{4} & r_{5} & r_{6} \\ r_{7} & r_{8} & r_{9} \end{matrix}]$ 、 $t = [\begin{matrix} t_{1} \\ t_{2} \\ t_{3} \end{matrix}]$ 分别为L坐标系到R坐标系的旋转矩阵和平移矩阵，可由双相机外参推导得到； $H_{l r}$ 为空间转移矩阵。

上述过程完成了对所有相机参数的初步估计。双相机DIC系统对标定精度要求较高，单个相机的高质量标定并不能有效提高系统的性能，多坐标系转换过程中会有误差的积累，需要对双相机的全局标定结果做进一步优化来提高系统的标定精度。捆绑调整（BA）是一种常用的全局优化方案，常用的优化算法有最速下降法、牛顿法、Levenberg-Marquardt（LM）法。其中，LM法结合了另外两种方法的优点，计算速度快、不受初始值的限制。所提方法选择LM法对全局标定参数进行优化，以重投影误差表示的最小代价函数形式为

m i n f = \sum_{i = 1}^{I} \sum_{j = 1}^{J} {‖x_{i j} - x_{i j}^{'} (M, k_{1}, k_{2}, k_{3}, p_{1}, p_{2}, H_{i})‖}^{2}

，（6）

式中： $i = 1,2, \dots, I$ ，为靶标图像数量； $j = 1,2, \dots, J$ ，为单张标定板角点数； $M$ 为相机的内参矩阵； $H_{i}$ 为第i幅标定图像的空间转移矩阵。从最小代价函数求解出重投影误差的最小值，从而达到优化全局标定参数的目的。

将式（1）和式（5）代入式（2），即可得到点 $P_{N}$ 的三维坐标。

2.3　基于双相机标定的镜面反射变换

为了实现被测物体的三维重建，需要知道图1（a）、（b）中B和C的信息，即平面镜前后实像和虚像的转换关系。为了提高测量精度以及简化计算过程，提出基于平面镜虚实像的转换关系，即通过平面镜实现从真实相机坐标系到虚拟相机坐标系之间的转换，进而实现从实像到虚像之间的转换，其转换原理如图3所示。L、R分别为真实的左、右相机， $L^{'}$ 、 $R^{'}$ 为L、R关于平面镜M1对称的虚拟相机； $o_{c} ‑ x_{c} y_{c} z_{c}$ 、 $o_{v} ‑ x_{v} y_{v} z_{v}$ 、 $O_{1} ‑ X_{1} Y_{1} Z_{1}$ 分别为相机L、 $L^{'}$ 以及平面镜M1的坐标系； $O_{2} ‑ X_{2} Y_{2} Z_{2}$ 为 $O_{1} ‑ X_{1} Y_{1} Z_{1}$ 的旋转坐标系； $[R_{v}^{(2)}, T_{v}^{(2)}]$ 为坐标系 $o_{v} ‑ x_{v} y_{v} z_{v}$ 到坐标系 $O_{2} ‑ X_{2} Y_{2} Z_{2}$ 的旋转和平移矩阵； $[R_{2}^{(1)}, T_{2}^{(1)}]$ 为坐标系 $O_{2} ‑ X_{2} Y_{2} Z_{2}$ 到坐标系 $O_{1} ‑ X_{1} Y_{1} Z_{1}$ 的旋转和平移矩阵； $[R_{c}^{(1)}, T_{c}^{(1)}]$ 为坐标系 $o_{c} ‑ x_{c} y_{c} z_{c}$ 到坐标系 $O_{1} ‑ X_{1} Y_{1} Z_{1}$ 的旋转和平移矩阵； $[R_{v}^{(c)}, T_{v}^{(c)}]$ 为 $o_{v} ‑ x_{v} y_{v} z_{v}$ 坐标系到坐标系 $o_{c} ‑ x_{c} y_{c} z_{c}$ 的旋转和平移矩阵。

图 3. 虚像和实像之间的转换关系

Fig. 3. Conversion relationship between virtual and real images

下载图片查看所有图片

由图3可知，L的相机坐标系与双目DIC系统的世界坐标系重合，因此，虚拟双目DIC系统的世界坐标系到真实双目DIC系统的世界坐标系的转换关系等价于 $L^{'}$ 和L的相机坐标系的转换关系，即只需计算出 $[R_{v}^{(c)}, T_{v}^{(c)}]$ 即可。

首先需要找出坐标系 $o_{c} ‑ x_{c} y_{c} z_{c}$ 到 $O_{1} ‑ X_{1} Y_{1} Z_{1}$ 的位置关系。可以通过如下标定步骤得到相机坐标系 $o_{c} ‑ x_{c} y_{c} z_{c}$ 到 $O_{1} ‑ X_{1} Y_{1} Z_{1}$ 之间的旋转和平移矩阵 $[R_{c}^{(1)}, T_{c}^{(1)}]$ ：1）将定制的已知厚度为 $d$ 的标定板背面紧贴在左平面镜表面；2）将标定板与左平面镜作为一个整体进行标定，得到相机坐标系 $o_{c} ‑ x_{c} y_{c} z_{c}$ 到标定板的旋转和平移矩阵 $[R_{c}^{(b)}, T_{c}^{(b)}]$ ；3）由于标定板到平面镜的旋转矩阵 $R_{b}^{(1)} = 0$ 和平移矩阵 $T_{b}^{(1)} = {[0,0, - d]}^{T}$ ，进而得到相机坐标系 $o_{c} ‑ x_{c} y_{c} z_{c}$ 到 $O_{1} ‑ X_{1} Y_{1} Z_{1}$ 之间的旋转和平移矩阵 $[R_{c}^{(1)}, T_{c}^{(1)}] = [R_{c}^{(b)}, T_{c}^{(b)} + T_{b}^{(1)}]$ ；4）同理可得到相机坐标系 $o_{c} ‑ x_{c} y_{c} z_{c}$ 到右平面镜的旋转和平移矩阵。

得到相机坐标系 $o_{c} ‑ x_{c} y_{c} z_{c}$ 到 $O_{1} ‑ X_{1} Y_{1} Z_{1}$ 之间的旋转和平移矩阵 $[R_{c}^{(1)}, T_{c}^{(1)}]$ 后，令 $α_{c 1}$ 、 $β_{c 1}$ 、 $γ_{c 1}$ 分别为L绕坐标系 $o_{c} ‑ x_{c} y_{c} z_{c}$ 的 $o_{c} x_{c}$ 轴、 $o_{c} y_{c}$ 轴、 $o_{c} z_{c}$ 轴旋转的欧拉角。每个欧拉角的正方向为绕其轴顺时针旋转（右手定则）的方向。根据旋转欧拉角 $α_{c 1}$ 、 $β_{c 1}$ 、 $γ_{c 1}$ 可分别得到坐标系 $o_{c} ‑ x_{c} y_{c} z_{c}$ 的 $o_{c} x_{c}$ 轴、 $o_{c} y_{c}$ 轴、 $o_{c} z_{c}$ 轴的旋转矩阵 $T (α_{c 1})$ 、 $T (β_{c 1})$ 、 $T (γ_{c 1})$ 。令坐标系 $o_{c} ‑ x_{c} y_{c} z_{c}$ 到坐标系 $O_{1} ‑ X_{1} Y_{1} Z_{1}$ 先绕 $o_{c} z_{c}$ 轴旋转，然后绕 $o_{c} y_{c}$ 轴旋转，最后绕 $o_{c} x_{c}$ 轴旋转，则旋转矩阵 $R_{c}^{(1)}$ 的计算公式为

R_{c}^{(1)} = [\begin{matrix} r_{1} & r_{2} & r_{3} \\ r_{4} & r_{5} & r_{6} \\ r_{7} & r_{8} & r_{9} \end{matrix}]

，（7）

相应地，平移矩阵可表示为

T_{c}^{(1)} = [\begin{matrix} t_{x} \\ t_{y} \\ t_{z} \end{matrix}]

。（8）

在得到坐标系 $o_{c} ‑ x_{c} y_{c} z_{c}$ 与 $O_{1} ‑ X_{1} Y_{1} Z_{1}$ 之间的位置关系后，还需要计算出坐标系 $o_{v} ‑ x_{v} y_{v} z_{v}$ 与 $O_{1} ‑ X_{1} Y_{1} Z_{1}$ 之间的位置关系，在计算该位置关系之前，引入一个中间坐标系，即坐标系 $O_{2} ‑ X_{2} Y_{2} Z_{2}$ 。如图3所示， $O_{2} ‑ X_{2} Y_{2} Z_{2}$ 为 $O_{1} ‑ X_{1} Y_{1} Z_{1}$ 的旋转坐标系。从坐标系 $O_{2} ‑ X_{2} Y_{2} Z_{2}$ 到 $O_{1} ‑ X_{1} Y_{1} Z_{1}$ ，只需将前者绕 $O_{2} Y_{2}$ 轴旋转180°，则两个坐标系之间的关系可表示为

[\begin{matrix} X_{1} \\ Y_{1} \\ Z_{1} \end{matrix}] = R_{2}^{(1)} [\begin{matrix} X_{2} \\ Y_{2} \\ Z_{2} \end{matrix}] = [\begin{matrix} c o s 180^{\circ} & 0 & - s i n 180^{\circ} \\ 0 & 1 & 0 \\ s i n 180^{\circ} & 0 & c o s 180^{\circ} \end{matrix}] [\begin{matrix} X_{2} \\ Y_{2} \\ Z_{2} \end{matrix}] = [\begin{matrix} - 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & - 1 \end{matrix}] [\begin{matrix} X_{2} \\ Y_{2} \\ Z_{2} \end{matrix}]

。（9）

由图3可知，旋转平移矩阵 $[R_{c}^{(1)}, T_{c}^{(1)}]$ 和矩阵 $[R_{v}^{(2)}, T_{v}^{(2)}]$ 关于平面镜M1对称， $α_{v 2}$ 、 $β_{v 2}$ 、 $γ_{v 2}$ 可分别表示坐标系 $o_{v} ‑ x_{v} y_{v} z_{v}$ 到坐标系 $O_{2} ‑ X_{2} Y_{2} Z_{2}$ 的绕 $o_{v} x_{v}$ 轴、 $o_{v} y_{v}$ 轴、 $o_{v} z_{v}$ 轴旋转的欧拉角，根据平面镜的成像规律以及欧拉角的性质，它们具有如下关系：

\{\begin{array}{l} α_{v 2} = α_{c 1} \\ β_{v 2} = - β_{c 1} \\ γ_{v 2} = - γ_{c 1} \end{array}

。（10）

根据式（7）和式（10），从坐标系 $o_{v} ‑ x_{v} y_{v} z_{v}$ 到坐标系 $O_{2} ‑ X_{2} Y_{2} Z_{2}$ 的旋转和平移矩阵 $[R_{v}^{(2)}, T_{v}^{(2)}]$ 可表示为

R_{v}^{(2)} = [\begin{matrix} r_{1} & - r_{2} & - r_{3} \\ - r_{4} & r_{5} & r_{6} \\ - r_{7} & r_{8} & r_{9} \end{matrix}], T_{v}^{(2)} = [\begin{matrix} - t_{x} \\ t_{y} \\ t_{z} \end{matrix}]

。（11）

计算出中间坐标系 $O_{2} ‑ X_{2} Y_{2} Z_{2}$ 与坐标系 $o_{v} ‑ x_{v} y_{v} z_{v}$ 的位置关系后，可以得到从坐标系 $o_{v} ‑ x_{v} y_{v} z_{v}$ 到坐标系 $O_{1} ‑ X_{1} Y_{1} Z_{1}$ 的旋转和平移矩阵：

\{\begin{array}{l} R_{v}^{(1)} = R_{v}^{(2)} R_{2}^{(1)} \\ T_{v}^{(1)} = T_{v}^{(2)} \end{array}

。（12）

联立式（7）~（9）和式（12）可得到坐标系 $o_{v} ‑ x_{v} y_{v} z_{v}$ 到坐标系 $o_{c} ‑ x_{c} y_{c} z_{c}$ 的旋转和平移矩阵：

\{\begin{array}{l} R_{v}^{(c)} = R_{v}^{(1)} R_{1}^{(c)} = R_{v}^{(2)} R_{2}^{(1)} R_{1}^{(c)} \\ T_{v}^{(c)} = T_{v}^{(1)} - R_{v}^{(1)} R_{1}^{(c)} T_{c}^{(1)} = T_{v}^{(2)} - R_{v}^{(2)} R_{2}^{(1)} R_{1}^{(c)} T_{c}^{(1)} \end{array}

。（13）

由于旋转矩阵为单位正交矩阵，故 $R_{1}^{(c)} = {[R_{c}^{(1)}]}^{T}$ ，则式（12）可改写为

\{\begin{array}{l} R_{v}^{(c)} = R_{v}^{(2)} R_{2}^{(1)} {[R_{c}^{(1)}]}^{T} \\ T_{v}^{(c)} = T_{v}^{(2)} - R_{v}^{(2)} R_{2}^{(1)} {[R_{c}^{(1)}]}^{T} T_{c}^{(1)} \end{array}

，（14）

故虚像到实像的旋转和平移矩阵为 $[R_{v}^{(c)}, T_{v}^{(c)}]$ 。

3　实验和结果

3.1　实验装置

实验以空心六角铝棒零件为研究对象，用所提方法对其热变形进行研究。双平面镜辅助双目DIC系统的实验装置如图4所示，包括LED冷光源、加热装置（高温炉）、红外线电子测温仪、两个CCD摄像机与镜头（相机由德国Allied Vision Technologies公司生产，型号为MANTA G-201B，最高帧率为30 frame/s，分辨率为1624 pixel×1234 pixel；镜头由德国Schneider-Kreuznach公司生产，型号为Cinegon-1.4/12，光圈范围为1.4~11.0）、用于图像采集和处理的个人计算机、两个光学反射镜及镜架等。铝棒高为78.5 mm，边长为14.0 mm，内空心直径为11.8 mm。为提高测量精度^［25］，实验时，需对零件喷涂散斑，如图5所示。其中A面为零件前表面，B面为零件左后侧面，C面为零件右后侧面，即为图1（a）、（b）中的A、B、C。

图 4. 实验装置示意图。（a）主要设备；（b）高温炉实物图

Fig. 4. Schematic of experimental setup. (a) Main equipment; (b) physical map of high temperature furnace

下载图片查看所有图片

图 5. 零件的实物图。（a）原型图；（b）三个侧面示意图

Fig. 5. Physical map of the part. (a) Prototype diagram; (b) three side schematic

下载图片查看所有图片

3.2　实验过程

实验前，将加热前的零件置于图4所示的测量系统中，为防止加热后零件热辐射对镜面的影响，将零件与镜面保持一定距离。将两个反射镜的夹角调整为120°，同时对两个摄像机的相对位置进行调整，以使零件的全景可以被镜面辅助双目DIC系统捕捉到，并标记零件位置、完成相机标定。将零件置于恒温炉中加热至350 ℃，然后将零件放置于测量系统的标记位置自然冷却，其间使用测温仪对零件温度进行检测。等零件温度降至310 ℃时，用相机进行拍摄，相机的采样频率设置为1 Hz，待零件温度降至20 ℃时相机停止拍摄。按照设置的采样频率1 Hz共采集到1290帧变形图像，参考图像为310 ℃时采集到的图像，目标图像为20 ℃时采集到的图像。

3.3　实验结果

3.3.1　三维重建结果

通过相机标定，且已知标定板的厚度为2 mm，由所提方法可得到左平面镜中虚像到实像的旋转矩阵为 $[\begin{array}{r} 0.1188 & - 0.6953 & - 0.7088 \\ - 0.6953 & 0.4514 & - 0.5593 \\ - 0.7088 & - 0.5593 & 0.4298 \end{array}]$ ，平移矩阵为 $[\begin{matrix} - 96.3091 \\ - 75.9894 \\ - 77.4696 \end{matrix}]$ 。同理可得到右平面镜中虚像到实像的旋转矩阵为 $[\begin{array}{r} 0.9843 & - 0.0612 & 0.1657 \\ - 0.0612 & 0.7621 & 0.6446 \\ 0.1657 & 0.6446 & - 0.7463 \end{array}]$ ，平移矩阵为 $[\begin{matrix} 10.3276 \\ 40.1662 \\ - 108.8145 \end{matrix}]$ 。最后，将平面镜前的实像点以及由左、右平面镜后的虚像点转化得到的实像点统一到一个世界坐标系中，可得到图6所示的零件外表面轮廓结果。

图 6. 零件外表面的重建过程。（a）三维散点图；（b）三维轮廓效果图；（c）三维重建表面与拟合曲面偏差云图

Fig. 6. Reconstruction process of outer surface of the part. (a) 3D scatter plot; (b) 3D outline sketch; (c) 3D reconstruction of the surface versus the fitted surface deviation cloud map

下载图片查看所有图片

重建表面与拟合表面的偏差如图6（c）所示，最大重建偏差大约为0.004 mm，而且主要集中在拼接区域，其他区域的偏差比较稳定。立体重建角度也证明该系统的性能比较稳定，精度较高。

3.3.2　热变形测量结果

参考图像为310 ℃时采集到的零件图像，目标图像为20 ℃时采集到的零件图像。以y方向上的两个中间边缘点M、N作为零件在y方向上的两个特征点，M、N之间的距离即为零件的高度，如图7（a）所示。用基于DIC的方法找到全部图像对应的特征点，由系统导出零件A面第一帧和最后一帧图像对应的y方向上的高度，结果如表1所示。零件A面y方向上高度的变化量（热变形）为 $Δ l_{y} = 0.5557 m m$ 。

图 7. 零件3个外侧表面的热变形结果。（a）零件A面热变形前（310 ℃）和热变形后（20 ℃）的位移云图；（b）零件B面的热变形后位移云图；（c）零件C面的热变形后位移云图

Fig. 7. Thermal deformation results of three outer surfaces of the part. (a) Displacement cloud diagram of part A surface before (310 ℃) and after (20 ℃) thermal deformation; (b) displacement cloud diagram of part B surface after thermal deformation; (c) displacement cloud diagram of part C surface after thermal deformation

下载图片查看所有图片

表 1. 起止温度下零件A面、B面和C面的高

Table 1. Height of parts A, B, and C at starting and ending temperatures

Measured temperature /℃

Height of A surface /mm

Height of B surface /mm

Height of C surface /mm

310

78.4999

77.9442

78.4987

77.9520

78.4978

77.9389

查看所有表

对于零件的B、C面，由于所拍摄的是经平面镜所成的虚像，从系统导出虚像点的坐标，再用所提方法进行虚实转换，将虚像点转换成双目立体DIC系统的世界坐标系下的实像点，最后进行热变形的计算。用所提方法可以分别得到零件B、C面竖直方向上第一帧和最后一帧图像的高度平均值，结果如表1所示。图7（b）、（c）分别为零件B面热变形后位移云图、C面热变形后位移云图。可以看到，零件B面高度的最大变化量（热变形）为 $Δ {\bar{l}}_{l} = 0.5467 m m$ ，C面高度的最大变化量（热变形）为 $Δ {\bar{l}}_{r} = 0.5589 m m$ 。最后，将用系统处理的零件A面与所提方法处理的零件B、C面在降温过程中的高度变化量进行直观对比，结果如图8所示，从第0帧到第1290帧的每30帧做一次对比，共有44帧对比图像。

图 8. 降温过程零件A面、B面和C面高的变化曲线

Fig. 8. Variation curves of heights of surfaces A, B, and C during the cooling process

下载图片查看所有图片

3.3.3　热变形仿真结果

为验证所提方法在热变形测量中的可行性和准确性，用ANSYS软件对零件的热变形过程进行仿真，空心六角铝棒零件的物理参数如表2所示。

表 2. 铝的物理参数

Table 2. Physical parameters of aluminum

Temperature /℃

Density /

（g·cm^-3）

Thermal conductivity /

［W·（m·K）^-1］

Coefficient of thermal expansion

Specific heat

capacity /

［J·（g·K）^-1］

Elastic modulus /GPa

Possion ration

2.71

217.7

23.6

0.88

68.9

0.33

查看所有表

使用有限元分析软件ANSYS对零件进行网格划分，综合考虑精度和效率后，设置网格单元边长为1 mm，共得到31600个四面体网格，如图9（a）所示。在热变形实验中，零件竖直放置在测量系统中自然冷却，故仿真时对零件的一端进行约束，选择稳态仿真，零件温度从20 ℃升到310 ℃的仿真结果如图9（b）所示。可以看到，零件竖直方向上最大的热变形量为0.5500 mm，表明零件在290 ℃温差时沿z方向的热变形量为 $Δ L = 550.0 μ m$ 。使用所提方法计算得到的零件A、B、C面沿高度方向上的热变形的平均值为 $Δ \bar{l} = 552.8 μ m$ ，该均值与仿真结果的绝对误差为 $e = Δ \bar{l} - Δ L = 2.8 μ m$ ，相对误差仅为 $ε = 0.51 %$ 。可见，利用所提方法计算的空心六角铝棒零件的热变形量与仿真结果接近，证明了所提方法可以用于零件的热变形全场测量。

图 9. 利用ANSYS对零件进行网格划分及热变形仿真。（a）零件的有限元网格划分图；（b）零件热变形仿真结果

Fig. 9. Meshing diagram and thermal deformation simulation on the part by ANSYS. (a) Finite element meshing diagram of the part; (b) simulation result of thermal deformation of the part

下载图片查看所有图片

4　结论

提出一种利用双平面镜辅助双目DIC系统来实现全场变形测量的新方法，并以空心六角铝棒零件为测量对象，测量其从310 ℃冷却至20 ℃过程中3个外侧表面的热变形结果，并与有限元软件仿真结果进行对比。结果表明，零件3个外侧表面沿高度方向的变化曲线基本重合，与仿真结果接近，利用所提方法计算得到的零件A、B、C面沿高度方向上的热变形平均值与仿真结果的绝对误差为 $2.8 μ m$ ，相对误差仅为 $0.51 %$ 。可见，所提方法既克服了单相机不能实现全场测量的局限性，又避开了多目DIC系统的复杂性。

参考文献

[1] 杜鉴昕, 赵加清, 王海涛, 等. 针对裂尖变形场测量的包络单元局部数字图像相关方法[J]. 光学学报, 2022, 42(1): 0112003.

Du J X, Zhao J Q, Wang H T, et al. Envelope element local digital image correlation method for crack tip deformation field measurement[J]. Acta Optica Sinica, 2022, 42(1): 0112003.

[2] 王凡, 赵亮, 吴晓东, 等. 基于CT及数字图像相关法的混凝土变形场测量分析[J]. 激光与光电子学进展, 2020, 57(20): 200401.

Wang F, Zhao L, Wu X D, et al. Measurement and analysis of concrete deformation field based on CT and digital image correlation method[J]. Laser＆Optoelectronics Progress, 2020, 57(20): 200401.

[3] 潘济宇, 张水强, 苏志龙, 等. 基于数字图像相关的水下螺旋桨三维变形测量[J]. 光学学报, 2021, 41(12): 1212001.

Pan J Y, Zhang S Q, Su Z L, et al. Measuring three-dimensional deformation of underwater propellers based on digital image correlation[J]. Acta Optica Sinica, 2021, 41(12): 1212001.

[4] Palanca M, Tozzi G, Cristofolini L. The use of digital image correlation in the biomechanical area: a review[J]. International Biomechanics, 2016, 3(1): 1-21.

[5] Yoneyama S, Ueda H. Bridge deflection measurement using digital image correlation with camera movement correction[J]. Materials Transactions, 2012, 53(2): 285-290.

[6] Ji Y, Li Z, Liu L P, et al. A deformation detection method for aircraft skin on uniform pressure by using speckle image correlation technology[J]. Measurement, 2020, 154: 107525.

[7] Tiwari V, Sutton M A, McNeill S R, et al. Application of 3D image correlation for full-field transient plate deformation measurements during blast loading[J]. International Journal of Impact Engineering, 2009, 36(6): 862-874.

[8] Neri P, Paoli A, Razionale A V, et al. Low-speed cameras system for 3D-DIC vibration measurements in the kHz range[J]. Mechanical Systems and Signal Processing, 2022, 162: 108040.

[9] 吴荣, 刘依, 周建民, 等. 数字图像相关在旋转物体全场应变测量中的应用[J]. 光学学报, 2020, 40(13): 1312005.

Wu R, Liu Y, Zhou J M, et al. Full-field strain measurement of rotating object using digital image correlation[J]. Acta Optica Sinica, 2020, 40(13): 1312005.

[10] 朱飞鹏, 陆润之, 白鹏翔, 等. 双反射镜成像用于二维数字图像相关应变精度提升[J]. 光学学报, 2019, 39(12): 1212002.

Zhu F P, Lu R Z, Bai P X, et al. Enhancement of strain measurement accuracy of two-dimensional digital image correlation based on dual-reflector imaging[J]. Acta Optica Sinica, 2019, 39(12): 1212002.

[11] 单宝华, 霍晓洋, 刘洋. 一种极线约束修正数字图像相关匹配的立体视觉测量方法[J]. 中国激光, 2017, 44(8): 0804003.

Shan B H, Huo X Y, Liu Y. A stereovision measurement method using epipolar constraint to correct digital image correlation matching[J]. Chinese Journal of Lasers, 2017, 44(8): 0804003.

[12] 王永红, 朱奕磊, 高启学, 等. 基于数字图像相关法的空间目标位姿测量[J]. 光学学报, 2022, 42(8): 0812001.

Wang Y H, Zhu Y L, Gao Q X, et al. Position and pose measurement of spatial object based on digital image correlation[J]. Acta Optica Sinica, 2022, 42(8): 0812001.

[13] Orteu J J, Bugarin F, Harvent J, et al. Multiple-camera instrumentation of a single point incremental forming process pilot for shape and 3D displacement measurements: methodology and results[J]. Experimental Mechanics, 2011, 51(4): 625-639.

[14] Li J R, Yang G B, Siebert T, et al. A method of the direct measurement of the true stress-strain curve over a large strain range using multi-camera digital image correlation[J]. Optics and Lasers in Engineering, 2018, 107: 194-201.

[15] Wang Y Q, Lava P, Coppieters S, et al. Application of a multi-camera stereo DIC set-up to assess strain fields in an erichsen test: methodology and validation[J]. Strain, 2013, 49(2): 190-198.

[16] Malowany K, Malesa M, Kujawińska M, et al. Application of 3D digital image correlation for development and validation of FEM model of self-supporting arch structures[J]. Applied Sciences, 2019, 9(7): 1305.

[17] Denys K, Coppieters S, Seefeldt M, et al. Multi-DIC setup for the identification of a 3D anisotropic yield surface of thick high strength steel using a double perforated specimen[J]. Mechanics of Materials, 2016, 100: 96-108.

[18] Malowany K, Malesa M, Kowaluk T, et al. Multi-camera digital image correlation method with distributed fields of view[J]. Optics and Lasers in Engineering, 2017, 98: 198-204.

[19] Spera D, Genovese K, Voloshin A. Application of stereo-digital image correlation to full-field 3-D deformation measurement of intervertebral disc[J]. Strain, 2011, 47: e572-e587.

[20] Badel P, Genovese K, Avril S. 3D residual stress field in arteries: novel inverse method based on optical full-field measurements[J]. Strain, 2012, 48(6): 528-538.

[21] Barone S, Neri P, Paoli A, et al. Structured light stereo catadioptric scanner based on a spherical mirror[J]. Optics and Lasers in Engineering, 2018, 107: 1-12.

[22] Pan B, Chen B. A novel mirror-assisted multi-view digital image correlation for dual-surface shape and deformation measurements of sheet samples[J]. Optics and Lasers in Engineering, 2019, 121: 512-520.

[23] Bin C, Bing P. Calibrating mirror-assisted multi-view digital image correlation system using a speckled planar object[J]. Measurement Science and Technology, 2021, 32(3): 034008.

[24] Zhang Z Y. A flexible new technique for camera calibration[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2000, 22(11): 1330-1334.

[25] Chen Z N, Shao X X, Xu X Y, et al. Optimized digital speckle patterns for digital image correlation by consideration of both accuracy and efficiency[J]. Applied Optics, 2018, 57(4): 884-893.

李桂华, 李涛, 孙卫庆, 龚启诚, 王晓宇, 张梅. 基于双平面镜双相机组合测量全场变形的方法[J]. 光学学报, 2023, 43(2): 0212007. Guihua Li, Tao Li, Weiqing Sun, Qicheng Gong, Xiaoyu Wang, Mei Zhang. Method for Measuring Full-Field Deformation Based on Double-Sided Mirror and Dual-Camera Combination[J]. Acta Optica Sinica, 2023, 43(2): 0212007.