夏瀚笙,沈 峘,王 莹,刘敦强
(南京航空航天大学 能源与动力学院, 南京 210016)
三维数字图像相关法的匹配策略和应变场计算
夏瀚笙,沈 峘,王 莹,刘敦强
(南京航空航天大学 能源与动力学院, 南京 210016)
针对在三维数字图像相关法的三维重建过程中存在着两种不同的匹配策略,利用三维仿真散斑图模拟变形,分析研究二者在不同场景下应用的品质。现有的三维应变场计算有直接对三维位移场计算和对位移场投影后再进行计算两种方法,这两种方法存在忽略局部特征、未去除噪声的缺陷。为了得到精度较高的应变场,采用基于最小二乘拟合的微单元体应变计算方法。仿真分析结果表明:利用该计算方法求得的应变场精度明显提高。
三维数字图像相关法;双目视觉;匹配策略;三维应变场
Abstract: In three-dimensional Digital Image Correlation (DIC), there are two different matching strategies in the process of three-dimensional reconstruction. This paper summarizes the application of these two matching strategies in different scenarios by using the three-dimensional simulation speckle maps to simulate the deformation. There are two methods to calculate the three-dimensional displacement field and to calculate the displacement field directly. The two methods have the disadvantages of ignoring the local feature and removing the noise. In order to obtain the strain field with high precision, the strain calculation method of the micro unit body based on the least squares fitting is used in the calculation of three-dimensional strain fields. The Simulation analysis results show that the strain field precision obtained by using this strain calculation method is reliable.
Keywords: three-dimensional DIC; binocular stereo vision; match strategy; three-dimensional strain computation
随着计算机技术和相机产业的发展,数字图像相关法(DIC)应运而生[1]。与传统的激光干涉原理测量法相比,数字图像相关法具有设备简单、非接触式、全场测量等优点。数字图像相关法在二维测量[2](即单目视觉)中已相对成熟,但平面应变测量在实际工程应用中有较大的局限性,二维数字图像相关法仅适用于平面变形测量,测量过程需要相机光心与被测件表面垂直,且对离面位移敏感[3]。近年来,三维数字图像相关法(3D-DIC)的相关研究成为数字图像相关法领域内的热点。三维数字图像相关法[4]是随着双目视觉[5]发展而广泛应用的。该方法利用双目视觉来测量三维表面应变,且对相机的位置没有太高的要求,能克服二维数字图像相关法无法测量离面位移的缺陷。
三维数字图像相关法主要包括相机标定、图像匹配、三维重建以及三维位移、应变计算等过程。众多学者在三维数字图像相关法的立体匹配策略和应变场计算精度方面提出了许多行之有效的方法。在立体匹配方面,Helm等[6]在小平面假设条件下通过搜索两幅图像子区相关系数最大值的候选空间平面参数完成匹配任务,并通过铝板拉伸试验进行了分析验证。潘兵等[7]提出以左相机采集的第1张图像作为参考图像,该图既是左相机匹配序列的参考图像,也是右相机匹配系列的参考图像;Tang Z Z等[8]提出以左相机采集的一系列图像作为参考图像,对每个左右相机采集对应的两张图像之间做匹配计算,然后利用相机标定参数三维重建。本文就上述两种匹配策略在不同场景下的适用性进行分析研究。在应变场计算方面,Lu H等[9]引入二阶形函数进行物体表面变形的描述,提高非均匀变形位移场和应变场计算精度,提高位移计算精度。Sun Y等[10]分析匹配子区大小对位移场精度的影响。2003年,Schreier[11]首次直接采用二阶形函数拟合左右相机所采集图像之间的变形,说明形函数拟合变形的方法优于基于小平面假设的投影求解方法。2006年,Siebert T[12]分析三维数字图像相关法中标定所产生误差的影响,其计算的位移场误差在0.02 pixel以内,应变误差在500 με左右,证明位移量在50 pixel以内,位移场测量误差呈线性增加,故建议增加径向畸变模型来减小位移场的测量误差。2011年,Helfrick等[13]对比有限元仿真结果和数字图像相关法的计算结果,说明三维数字图像相关法的准确性得到验证,证明三维数字图像相关法结果和有限元结果相符合。清华大学的潘兵在散斑质量评价[14]中提出基于最小二乘拟合[15]的全场应变计算方法,减小位移场引起的粗大误差,并且提出应变场降噪的策略[16]。在三维应变计算过程[17-18]中,直接采用三维位移场通过最小二乘法拟合方法计算应变场,非均匀变形的细节信息将被平滑。高越[19]在计算三维应变时先将位移场进行投影,然后进行柯西应变计算,位移场未进行平滑处理,应变场计算过程会将位移场噪声恶意放大。
本文针对三维数字图像相关法中的立体匹配策略,分析在大变形和小变形两种情形下的适用性。为提升三维数字图像相关法应变场计算精度,提出将微单元体投影到拟合平面转化成二维应变场,结合局部最小二乘拟合计算三维表面应变的方法。
1.1 双目视觉模型
计算机双目视觉技术属于光学测量中的一种被动式测量方法,利用两个相机从不同的视角观察同一物体,获取在不同视角下的图像,然后根据图像间的匹配关系,利用三角测量原理计算出像素间的偏移量,从而获取物体三维信息。双目视觉在建立三维模型的过程中,需要利用坐标系来描述模型中测量点的位置。图1表示了双目视觉系统中多个坐标系系统。o-uv为图像像素坐标系,o-XuYv为图像物理坐标系,Ocl-xclycl、Ocr-xcrycr分别为左、右相机坐标系,Ow-XwYwZw为世界坐标系。物体上点P(xw,yw,zw)对应左、右相机拍摄图像上点P1、P2。由P1、P2的图像坐标系坐标结合相机的内、外参数,即可求出P的世界坐标(xw,yw,zw)。
图1 双目视觉系统坐标系
1.2 相机标定
本文采用张正友[20]的平面标定法,相机模型运用传统的针孔成像模型,公式如下:
(1)
式中:s为任意的比例系数;R、t分别为从世界坐标系到摄像机坐标系的旋转矩阵和平移向量,即R和t为相机的外部参数,A为相机的内部参数。
实际相机镜头的加工制造误差以及安装误差会导致镜头畸变,因此简单的针孔模型无法准确地描述镜头成像关系。当只考虑径向畸变时,利用Tsai[21]两步标定法:
δx=(u-u0)[k1r2+k2r4]
δy=(v-v0)[k1r2+k2r4]
(2)
其中:k1,k2分别为1阶、2阶径向畸变系数;r为像素点到主点(u,v)的距离。
2.1 匹配策略
3D-DIC采用双目视觉技术,匹配过程涉及自匹配和互匹配的问题。三维重建阶段中存在两种匹配策略,一种是左相机所采集的系列图像之间的匹配;另一种是左相机和右相机所采集的图像之间的匹配。
图2(a)策略1中,以左相机采集的第一张图像作为参考图像,第一张图像既是左相机匹配序列的参考图像,也是右相机匹配系列的参考图像。此策略即为潘兵[7]所使用的计算策略。图2(b)策略2中,在左相机采集的图像匹配过程中,以左相机采集的第一幅图像为参考图像。在左、右相机匹配的过程中,是以左相机采集的一系列图像作为参考图像,左、右相机采集对应的每两张图像之间做匹配计算,然后利用相机标定参数三维重建。此策略即为西安交通大学[8]所使用的策略。
2.2 仿真分析
仿真分析硬件配置采用Intel(R) Core(TM) i7-6700 CPU,内存16G。采用Matlab编程软件。
采用文献[22]中的三维仿真模型,图像为随机生成的散斑图像,分辨率为700×700 pixels,散斑颗粒数为19 600,散斑颗粒半径为3 pixel。假设物体是一个ZW=1 200 mm的平面。其中相机的内、外参数为:
在图像匹配计算中采用的网格点间距为10 pixel,图像匹配半径15 pixel,牛顿迭代最大收敛次数20次。
图2 两种匹配策略
2.2.1 平移测量
沿x向生成间隔为0.05pixel的平移图像共100张,在三维空间里对应的实际位移为0~0.6 mm,间隔为0.006 mm。仿真分析结果如图3所示。
图3 策略1和策略2在小变形情况下的测量误差
如图3所示,在位移量较小的情况下,策略1精度高于策略2。这是由于在小变形的情况下,策略1采用参考图像是用左相机图像的第一张图像和右相机图像进行匹配,其视差和变形量较小;而方策略2中,只有第一张图像和第二张图像的匹配,没有引进计算点的误差,在后续的计算中,因为左相机图像作为参考图像,参与匹配的点要求是整像素的,此时在选取计算点时就有取整的步骤,此时就会引起误差增加。
2.2.2 旋转分析
采用仿真图为30°的散斑图,其中计算点为504个,网格间距为10 pixel,计算子区半径为15 pixel。根据仿真分析结果可知:在图像旋转30°时,匹配策略1有130个点计算不收敛,且收敛点的位置也不合理。而经过计算,匹配策略2中所有的初值估计点都收敛,即使在30°旋转时,也可以经过去除异常点得到比较理想的数据。
因此,策略1在大旋转的情况下不适合用做立体匹配。主要原因在于:参考图像采用的是左相机第一张图像,旋转角度越大,参考图像和目标图像的差异越大。其中,既包括本身的变形误差,还包括双目相机的视差,因此导致收敛状况差。旋转仿真分析中,策略1在旋转30°时有较大误差,因此不做计算。在三维重建结果误差图上可以明显看出:策略1不适用于变形较大的立体匹配过程;而策略2因为采用对应的左、右相机做匹配,在进行左、右相机匹配分析时,相机的视差较小,因此可以完全迭代收敛,适用于大旋转大位移的计算情况。
图4 旋转30°时,策略1和策略2的匹配计算对比
图5 两种策略三维重建结果
3.1 微单元体投影定理
弹性体的变形一般用微六面体单元表述。对于本文关注的小变形问题,为简化分析,将微单元体分别投影到oxy,oyz,ozx平面,如图6所示。
图6 微单元体
图6中:MA、MB、MC变形前分别为x、y、z坐标轴平行的棱边;MA′、MB′、MC′为MA、MB、MC变形后对应的棱边。分别用εx、εy、εz表示正应变,γxy、γyz、γzx表示切应变,则有
(3)
本文的应变计算针对物体表面的应变状态,oyz、ozx平面上的应变量相对于oxy面上的应变量很小,可以忽略不计。通过研究oxy面上的投影来分析物体表面的三维应变状态。
设ma、mb分别为MA、MB的投影,m′a′、m′b′分别为M′A′、M′B′的投影。微单元体的棱边长为dx、dy、dz,M点的坐标为(x,y,z)、u(x,y,z)、v(x,y,z)分别表示M点x,y方向的位移分量。则A点的位移为u(x+dx,y,z),v(x+dx,y,z);B点的位移为u(x,y+dy,z),v(x,y+dy,z)。用泰勒级数形式将A、B两点的位移展开,并且略去2阶以上的小量,则A、B点的位移分别为:
(4)
(5)
(6)
同理
(7)
由此可以计算出弹性体内任意一点的微线段的相对伸长度,即正应变。若微线段伸长,则正应变εx、εy、εz为正,缩短则为负。
图7 微单元体在oxy面上的投影
3.2 微单元体最小二乘拟合应变法
最小二乘拟合方法计算应变即在计算应变时采用最小二乘拟合的计算方法求解形函数的过程。通常位移场数据包含系统误差、随机误差和粗大误差,数字图像相关法直接计算的位移场是有误差的,若直接对位移场进行差分应变计算会放大误差。最小二乘拟合的过程会对位移场数据进行平滑处理[18],计算示意图如图8所示。
图8中,蓝色标记点是匹配计算得到的位移场,红色点是大小为(2M+1)×(2M+1)的计算窗口。窗口的中心点即要计算的网格点,其余点只做辅助计算。窗口在位移场中的步长为1pixel,对二者重合部分进行拟合计算。假设形函数采用1阶形函数[15],那么用二元一次函数进行拟合:
(8)
式中:u和v为位移场;x、y的变化范围均为[-m,m];a0,…,b2为拟合系数,与应变场之间的关系是
(9)
求解出拟合方程的系数,即可得到应变值。
图8 局部最小二乘法应变计算示意图
在局部拟合过程中,拟合方程的矩阵形式如下:
(10)
其中只有6个系数为未知量,可利用最小二乘拟合方法解出系数值。
3.3 仿真分析
参考图像为随机生成的散斑图像,分辨率为700 pixel×700 pixel,散斑颗粒数为19 600,散斑颗粒半径为3pixel。假设物体是一个ZW=1 200 mm的平面。相机的内外参数为:
3.3.1 均匀变形仿真分析
本实验的散斑图沿x方向应变为2 000 με,此应变计算的区域半径为20 pixel。
如图9所示,x方向的位移场是一个平滑的斜面,应变场的真值为2 000 με,而实际计算得到的应变场围绕2 000 με上下波动,x方向的应变场计算结果为1 960 με,方差为58 με。y方向的位移场为0.5 μm以内,应变场的均值为50 με。
3.3.2 沿x方向呈sin应变仿真分析
利用反向映射法生成沿x方向非均匀变形图,x方向的位移场符合sin分布,即v=Asin(2πx/T),A=0.12 mm,T=200,则其应变的表达式为
ε=2πA/Tcos(2πx/T)
(11)
计算得到的位移场和应变场如图10所示。
图9 最小二乘拟合法
图10 最小二乘拟合位移场和应变场
仿真分析结果表明:对于三维重建中的两种匹配策略,在小变形情况下,匹配策略1优于匹配策略2;在大变形情况下,匹配策略2优于匹配策略1。
利用均匀应变和非均匀应变的三维仿真图,证明本文提出的三维结构表面应变的计算方法提高了应变计算的精度。
[1] 韩昌元.高分辨力空间相机的光学系统研究[J].光学精密工程,2008,16(11):2164-2172.
[2] SUTTON M A,WOLTERS W J,PETERS W H,et al.Determination of displacement using an improved digital correlation method[M].Brighton:The New International Economic Order and UNCTAD IV Bangladesh Institute of Development Studies,1977:133-139.
[3] 孟利波,金观昌,姚学锋.DSCM中摄像机光轴与物面不垂直引起的误差分析[J].清华大学学报:自然科学版,2006,46(11):1930-1932.
[4] LUO P F,CHAO Y J,SUTTON M A.Application of stereo vision to three-dimensional deformation analyses in fracture experiments[J].Optical Engineering,1994,33:3(3):990.
[5] BURIAN H M,NOORDEN G K V.Binocular Vision and Ocular Motility[J].Journal of Neuro-Ophthalmology,1975,52(9):614.
[6] HELM J D.Deformations in wide,center-notched,thin panels,part I:three-dimensional shape and deformation measurements by computer vision[J].Optical Engineering,2003,42(5):1293-1305.
[7] 潘兵,谢惠民,李艳杰.用于物体表面形貌和变形测量的三维数字图像相关方法[J].实验力学,2007,22(6):556-567.
[8] TANG Z Z,LIANG J,XIAO Z Z,et al.Three-dimensional digital image correlation system for deformation measurement in experimental mechanics[J].Optical Enginee-ring,2010,49(10):1298-1298.
[9] LU H,CARY P D.Deformation measurements by digital image correlation:Implementation of a second-order displacement gradient[J].Experimental Mechanics,2000,40(4):393-400.
[10] SUN Y,PANG J H L.Study of optimal subset size in digital image correlation of speckle pattern images[J].Optics & Lasers in Engineering,2007,45(9):967-974.
[11] SCHREIER H W.Investigation of two and three-dimensional image correlation techniques with applications in experimental mechanics[D].USA:University of South Carolina,2003.
[12] SIEBERT T.Error estimations of 3D digital image correlation measurements[J].Proceedings of SPIE -The International Society for Optical Engineering,2006:63410F-63410F-6.
[13] HELFRICK M N,NIEZRECKI C,AVITABILE P,et al.3D digital image correlation methods for full-field vibration measurement[J].Mechanical Systems & Signal Processing,2011,25(3):917-927.
[14] PAN B,LU Z,XIE H.Mean intensity gradient:An effective global parameter for quality assessment of the speckle patterns used in digital image correlation[J].Optics & Lasers in Engineering,2010,48(4):469- 477.
[15] XIE H,PAN B.Full-field strain measurement using a two-dimensional Savitzky-Golay digital differentiator in digital image correlation[J].Optical Engineering,2007,46(3):033601.
[16] PAN B,YUAN J,Xia Y.Strain field denoising for digital image correlation using a regularized cost-function[J].Optics & Lasers in Engineering,2015,65:9-17.
[17] YAN H,PAN B.Three-dimensional displacement measurement based on the combination of digital holography and digital image correlation[J].Optics Letters,2014,39(17):5166.
[18] PAN B,WU D,YU L.Optimization of a three-dimensional digital image correlation system for deformation measurements in extreme environments[J].Applied Optics,2012,51(19):4409- 4419.
[19] 高越.三维数字图像相关法的关键技术及应用研究[D].合肥:中国科学技术大学,2014.
[20] ZHANG Z Y.A flexible new technique for camera calibration[J].IEEE Transaction son Pattern Analysis and Machine Intelligence,2000,22(11):1330-1334.
[21] TAAI R Y.A versatile camera calibration technique for high-accuracy 3D machine vision metrology using off-the-shelf TV cameras and lenses[J].IEEE Journal on Robotics & Automation,1987,3(4):323-344.
[22] 毛建国,张佩泽,程百顺,等.基于反向映射法逆向描述数字散斑变形的方法[J].光电子·激光,2015(12):2433-2439.
(责任编辑陈 艳)
StudyonStereoMatchingStrategyand3-DStrainComputationinThree-DimensionalDigitalImageCorrelation
XIA Hansheng, SHEN Huan, WANG Ying, LIU Dunqiang
(College of Energy and Power Engineering, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China)
2017-05-03
航空科学基金资助项目(20120952022)
夏瀚笙(1993—),安徽芜湖人,硕士研究生,主要从事计算机视觉、深度学习在人脸识别检测方面的研究,E-mail:hanson_cha@163.com;通讯作者 沈峘,男,博士,主要从事图像分析与检测技术研究,E-mail:huan_shen@nuaa.edu.cn。
夏瀚笙,沈峘,王莹,等.三维数字图像相关法的匹配策略和应变场计算[J].重庆理工大学学报(自然科学),2017(9):110-118.
formatXIA Hansheng,SHEN Huan,WANG Ying,et al.Study on Stereo Matching Strategy and 3-D Strain Computation in Three-Dimensional Digital Image Correlation[J].Journal of Chongqing University of Technology(Natural Science),2017(9):110-118.
10.3969/j.issn.1674-8425(z).2017.09.018
TP391;O348
A
1674-8425(2017)09-0110-09