陈仁虹，梁 晋 ，叶美图，任茂栋，张继耀

（1. 西安交通大学 机械工程学院 机械制造系统工程国家重点实验室，陕西 西安 710049；2. 新拓三维技术（深圳）有限公司 创新实验室，广东 深圳 518060）

1 引 言

在材料领域，柔性薄膜是一种重要的材料，其地位越来越重要，各种不同功能的薄膜都有广泛的应用。薄膜材料的种类繁多，应用广泛，目前常用的有：半导体薄膜复合材料、金刚石薄膜、铁电薄膜、氮化碳薄膜和多层薄膜材料等。随着成膜技术的飞速发展，各种材料的薄膜化已经成为了一种普遍趋势[1]。由于薄膜材料很薄且其刚性不足极易变形，在成形过程中很难控制其承载大小，因此需要测量薄膜材料的成形极限曲线（Forming Limit Curves, FLC），从而能够更加方便有效地生产制造薄膜材料。为了测量出薄膜材料的成形极限曲线，设计制备了由PET、Nylon、Al foil和PP 4种材料通过一定的比例复合而成的柔性复合薄膜[2]。

材料成形技术是由不同应变路径条件下极限应变描绘构成的成形极限曲线，是区分材料成形过程中安全和破裂失效之间的界限，是当前最为广泛的评价材料成形性能的判据。传统的试验测量方法一般利用单个圆形图案在材料产生变形时近似变成椭圆的性质，采用工业软尺或工业显微镜直接测量椭圆的长轴和短轴长度近似获取单个圆所在区域的最大主应变和最大次应变。传统的试验测量方法不仅费时费力，而且不能准确地判断材料成形极限状态发生的时刻，不能保证成形极限的测量精度和成形极限曲线的拟合准确度等[3-4]。

由于材料成形极限试验中传统应变测量方法的不足，随着计算机技术和图像处理技术的发展，人们开始广泛采用坐标网格方法。坐标网格技术是对材料成形表面的应变分布和大小进行计算，是目前研究材料成形方面的主要方法[5-6]，2011年，西安交通大学梁晋[7]等发明了一种三维网格应变测量方法，解算出每个网格节点处的变形梯度张量，从而计算得到最大主应变和最大次应变；2012年，福州大学的梁伟[8]等提出一种适用于金属薄板的网格印制方法；2015年北京航空航天大学的朱凯博[9]等在网格应变测量中提取网格，先对图像进行边缘检测，然后求出边缘图像中每个网格的质心，利用质心重组应变网格计算网格应变。

尽管坐标网格方法在一定程度上满足材料表面应变测量的要求，但仍然存在一些不足：（1）坐标网格方法的网格绘制工作量大且精度一般；（2）只能对材料成形的最终状态进行测量，不能在材料成形过程中实现动态监测，且很难保证材料成形的最终状态正好是材料出现颈缩时的极限状态；（3）薄膜材料刚性差，极易变形，且不易固定，采用坐标网格方法不能准确地测量出薄膜材料的极限应变。

为了克服上述极限应变测量方法出现的不足，提出了一种基于双目立体视觉结合数字图像相关法（Digital Image Correlation, DIC）的材料成形极限应变测量方法。DIC是一种先进的光学非接触式测量方法[10-11]，具有对试验环境适应性好、材料测试范围大、操作简单等特点，在材料力学性能测量方面具有广泛的应用，能实时得到柔性薄膜材料变形全过程的应变值[12-13]。另外，经过多年的发展，计算机双目立体视觉技术和数字相关方法结合的研究和应用也比较成熟，从而间接地推动了数字散斑相关方法在三维变形测量领域的发展，DIC技术已经能够在实际生产中得到应用。

2 DIC原理与成形极限

2.1 基于摄影测量的双目相机标定

相机在拍摄物体时，物体点经过相机光心成像在像平面上，理想的投影成像模型是几何光学中的小孔成像模型，物方点、镜头光心和实际像点在一条直线上，在摄影测量中被表述为共线方程，但是在实际情况下还需要考虑镜头的各种畸变，因此共线方程[14-15]为：

式中 (X,Y,Z)为物方点在世界坐标系下的三维坐标，(XS,YS,ZS)为相机光学中心S的三维坐标，(x,y)为 图像点在像平面上的坐标，(x0,y0)为主点坐标，(dx,dy)为相机的镜头畸变量。具体见图1。

图1 摄影测量原理示意图Fig. 1 Principle diagram of photogrammetry

相机的内方位参数包括主点(x0,y0)、焦距f和相机的各种镜头畸变，镜头畸变包括径向畸变、切向畸变和薄棱镜畸变。本文的相机标定算法将所有可能的三阶径向畸变以及切向畸变和薄棱镜畸变都考虑在内，所采用的的畸变模型为：

其中K1、K2、K3为 径向畸变参数，B1、B2为切向畸变参数，E1、E2为薄棱镜畸变参数，再加上主点坐标(x0,y0)和 焦距f共10个参数，称之为10参数相机畸变模型，从而确定相机的内方位参数。

2.2 数字图像相关匹配

数字图像相关法是对CCD相机拍摄的前后两幅散斑图像进行相关计算，即变形前的参考图像和变形后的变形图像，在两幅图像中找出相对应的点。如图2所示，在参考图像中，取以待匹配点F为中心的( 2M+1)×(2M+1)大小的矩形子图像作为参考子图像，在变形图像中，通过基于种子点的图像匹配搜索算法，并按照预先定义的相关系数进行相关计算，可用引入光强系数的最小距离平方和系数函数[16]，其形式为：

数字图像相关的过程就是通过相关函数，寻找与参考子图像相似度最大的以F′为中心的目标子图像的过程。如图2所示，参考子图像中点F在目标子图像中的对应点为F′。

图2 数字图像相关原理图Fig. 2 Schematic diagram of the digital image correlation method

假设物体参考图像上有一点P(x,y)，变形图像上对应的点为P′(x′,y′)，在物体只发生刚体位移的情况下，两点之间的映射关系可用零阶映射函数描述：

式中u、v分别为点位移的水平分量和竖直分量。考虑到实际情况中变形子区也可能发生旋转和形变，因此参考子区和变形子区之间有一阶映射关系，包括旋转、平移、正应变和剪应变，即：

式（5）中 Δx和 Δy表示A点x、y方向上的距离。

数字图像相关法常被用于均匀变形的情况，当柔性薄膜材料在成形过程中变形太大时，变形图像与参考图像之前的相似性程度大大降低，不利于图像的相关性计算，匹配精度很低，甚至不能匹配。根据系列图像相邻状态变形的连续性，提出了一种针对柔性薄膜材料大变形与裂纹情况的弱相关分步匹配方法。DIC散斑图像匹配技术已经应用于实际社会生产实践中。

（1）大变形匹配

在数字图像相关法中，对于变形前后散斑图像上的某一点，为了寻找该点在变形后的散斑图像上的对应点，通常利用该点周围矩形子区内的散斑灰度信息，在变形图像上通过求取相关系数的最大值（或最小值）来寻找。在变形测量过程中，通常会得到连续变形的多个状态的图像（变形量不断增大）。现有数字图像相关法的匹配方法是每张图像都基于第一张图片进行匹配，从而计算出每张图像的应变值，这种在小变形时可以精确匹配，但在薄膜材料成形试验中，试件的变形量非常大，传统的匹配方法不适用，匹配精度很差[17]。

针对这种柔性薄膜材料大变形或裂纹的情况（图3），在原有的种子点匹配方法的基础上，进行相应的算法改进，提出了一种图像匹配基准自适应更新的弱相关分布匹配方法，如图4所示[18]，匹配步骤如下：

图3 薄膜材料大变形或出现裂缝Fig. 3 Large deformation or cracks of the thin film materials

图4 顺序逐帧基准匹配方法Fig. 4 Sequential frame-by-frame benchmark matching method

①将第0个状态图像作为优先匹配；

②若与第0个状态的散斑图像匹配失败，则与第1个状态的散斑图像匹配；

③若匹配失败，则依次类推，直到第n−1个状态的散斑图像匹配。

初始未变形状态的散斑图像优先匹配，极大程度上减少了单纯以上一帧为参考图像计算带来的累计误差，保证了匹配精度；对于未成功匹配的散斑图像进行顺序逐帧基准匹配，能够提高薄膜材料表面变形场的完整度。

（2）全局匹配

在柔性复合薄膜成形试验中试件会出现裂纹，此时散斑图像变形太大，不能完全使用传统的DIC方法进行图像匹配。含有裂纹的散斑图像的计算区域可分为两个部分，裂纹影响区和非裂纹影响区。在实际散斑图像上，非裂纹影响区域占据了大部分的图像幅面。在网格单元的边长、数目以及计算区域大小相等的情况下，局部DIC方法不管是在计算速度还是计算精度都要高于全局DIC方法，因此提出一种联合匹配方案，在非裂纹影响区域采用局部DIC常规方法进行匹配，在裂纹影响区域采用弱相关分步匹配方法进行匹配，从而提高匹配精度和计算效率。

2.3 成形极限曲线拟合

将XJTUDIC三维数字应变测量软件对采集的图像进行处理，得到试件在变形过程中各个状态的表面应变后，对这些应变值进行处理得到材料的极限应变曲线，步骤如下：

(1)确定试件在变形过程中发生颈缩的位置，找到首张破裂图的前一张（失效位置）；

(2)在首张破裂图的前一张上确定截线的位置，自动确定截线的最佳位置；

(3)利用截线上的各个点的应变值来进行反二次拟合，确定该试件一条截线的极限应变；

(4)利用同种材料各种不同形状试件上截线的极限应变值来确定薄膜材料的极限应变曲线。

2.3.1 紧缩区域位置确定

该截线上的最大主应变是在试件已失稳的情况下计算出的，并不是实际试件发生颈缩时的极限应变，若要找出失效位置，则需重新拟合：找出截线上最大主应变的最大值，并分别向截线左右两端取5个点，用二次曲线f(x)=ax2+bx+c拟合出试件的实际发生紧缩时的最大主应变，即为试件紧缩区域的位置。

2.3.2 最佳拟合区域确定

截线上最大主应变最大的几个点已失稳，不能用于参与极限应变拟合，将这部分点舍弃，用剩下的点参与拟合，实际用来参与拟合的点只有其中一部分，需要根据数据来确定：左半部分的右端点和右半部分的左端点。下面来介绍左半部分用来参与拟合的点的算法。

（1）拐点的确定

从2.3.1可得颈缩区域，由颈缩区域位置向左取20 mm范围内的点，找出其二阶导数最大的地方（最陡的地方），此位置即为左端点。为了确定二阶导数，取3个连续点 (x1,y1)、(x2,y2)、(x3,y3)；通过这3个点确定二次抛物线并使用两倍的二次项系数，这个值就是中间点(x2,y2)的二阶导数，即

这种方法使用超过3点得到的数据被称为“过滤的”二阶导数，为了得到相同的中间点，使用奇数的点数。在这种情况下使用最小平方抛物线拟合，中间点的二阶导数是平方项系数的两倍。这种方法对二阶导数有“过滤”作用，但不会对初始曲线产生影响。

（2）拟合区域宽度确定

由（1）也可以确定右部分的拐点位置，从而得到左右两个部分的最大主应变和最大次应变：左边内部的最大主应变 ε1,BL、左边内部的最大次应变 ε2,BL、右边内部的最大主应变 ε1,Br、右边内部的最大次应变 ε2,Br。拟合区域的宽度 ω应最少包含5个点，并按照下式计算：

式中

由以上方法可得到最佳拟合范围用来拟合，从而确定材料的极限应变，用于极限曲线拟合。

2.3.3 极限应变拟合

对最佳拟合范围中的 ε1和 ε2随距离的变化关系进行拟合，使用两个区域最大主应变的最佳拟合二阶反抛物线f(x)=1/(ax2+bx+c)计算应变ε1。 裂纹位置就是需要的极限应变ε1，同时计算出极限应变ε2，从而得到一个极限应变对( ε1,ε2)。应保证：①每条截线上最多只能得到一个极限应变对 (ε1,ε2)；②要剔出成形试验过程中不满足极限应变标准的点；③对应于最大主应变和最大次应变极值的两个节点之间的距离应小于试件的3倍厚度。

2.3.4 极限应变曲线拟合

根据成形试验，将所有同种材料不同形状的试件的极限应变对( ε1,ε2)进行处理：

(1)将属于相同试件（一个试件上取3条截线，每条截线的数据拟合出一个极限应变对(ε1,ε2)）和具有相同形状试件的极限应变对(ε1,ε2)归类；

(2)计算出相同试件和每种形状试件的极限应变的平均值，将ε1作 为Y轴，ε2作为X轴；

(3)按照X轴从小到大的顺序，将各个极限应变对连起来，则连起来的折线就是该种材料的极限应变曲线。

3 系统组成及测量流程

3.1 系统组成

试验设备为材料成形机及配套计算机和XJTUDIC三维数字应变测量软件，相机的型号为Basler Usb3.0，分辨率为2448 pixel×2048 pixel，镜头采用Schneider焦距为25 mm的定焦镜头，型号为RICOH，2个蓝光LED光源，一台同步控制箱，相机三脚架，一个力传感器及其他辅助器件等，系统组成如图5所示。

图5 实验设备及测量软件Fig. 5 Experimental device and software

3.2 测量流程

将制备好的试件放于成形试验机上用夹具夹紧固定，开启成形试验机，同时打开相机采集试件在成形过程中各个状态的变形图像。

用XJTUDIC三维数字应变测量软件对采集到的图像进行处理，经过散斑区域选取、图像种子点选择、相关匹配、散斑区域位移和应变计算等，得到试件在整个成形过程各个状态的表面应变场，测量流程如图6所示。

图6 测量流程Fig. 6 Measurement process

4 试验及分析

4.1 精度验证

为了验证本文方法的可行性和准确性，在相同的试验条件下，将采用坐标网格方法测量相同材料的Q235钢试件的极限应变描绘的极限应变曲线，与本文方法进行对比。

为了描述完整的成形极限曲线至少需要5组不同几何形状的试件（目的是得到从单向到等双向均匀拉伸的极限应变曲线）。Q235金属试件制备的尺寸由ISO 12004附录A确定，制备的试件如图7所示。

图7 Q235钢试件Fig. 7 Q235 steel specimen

本文方法可以计算出Q235钢试件的极限应变，实验中实际测得的应变数据很多，将同一试件上不同截线或相同形状试件的极限应变进行归类合并，算出其平均值，得到Q235金属试件的极限应变对( ε1,ε2)，如表1所示，然后将各个极限应变对在二维坐标轴上连起来，即为极限应变曲线，如图8所示。

图8 两种方法结果对比Fig. 8 Comparison with two different methods

从图8中可见，本文方法和坐标网格方法测得的两条极限应变曲线的应变分布很接近，但是坐标网格方法是在已经破裂的试件上进行极限应变测量，测量值比实际临界值大；目前，FLC曲线的获取主要是以网格法为标准，而本文方法测得的极限应变比网格方法的结果小，更接近于试件的实际极限应变，精度能够提高0.02%，由此可知采用本文方法测量材料的极限应变曲线更加准确，可用于测量柔性复合薄膜材料的极限应变曲线。

4.2 测量案例

为了测量出薄膜材料的成形极限曲线，设计制备了由PET、Nylon、Al foil和PP 4种材料通过一定的比例复合而成的柔性复合薄膜，参考ISO 12004标准制备了柔性复合薄膜试件并喷涂散斑，如图9所示。

图9 制备柔性复合薄膜并喷涂散斑Fig. 9 Pliable composite film and speckle spray

采用本文方法，将喷有散斑的柔性复合薄膜试件进行成形试验，采用八步法[12]进行双目相机标定，标定图片如图10所示,标定参数如表2所示。

表2 相机标定内参数Tab. 2 Camera calibration internal parameters

图10 相机标定Fig. 10 Camera calibration

得到物方控制点的三维坐标及相机的标定结果后，通过后方交会计算即可得到左右相机的外方位参数[R|T]。

采集柔性复合薄膜成形试验各个状态的变形图像，然后用XJTUDIC三维数字应变测量软件对采集到的图像进行处理，得到柔性复合薄膜试件在整个成形过程各个状态的表面应变场，如图11所示。

图11 不同状态的表面应变分布Fig. 11 Surface strain distribution in different states

4.2.1 失效位置确定

在成形试验中，利用力传感器，可以输出试件受到的力随时间变化的数据；由XJTUDIC软件可以得到薄膜试件各个状态的应变值，由应变变化来确定失效位置，而传感器有延时性，失效位置主要用应变来确定，力值来辅助确定薄膜试件的失效位置。

随着成形试验的进行，柔性薄膜试件的变形会越来越大，在失效位置前后，柔性薄膜试件的应变会发生骤变，如图12所示，只要找出其前后两个状态变化最大的位置即可。为了避免应变计算时出现误差，观察数据可发现薄膜试件破裂之后，其应变值基本不变，找出其中前后变化差值最大的5个位置，然后再计算出其后面5个点的方差，计算方差最小的位置，即为失效位置。在计算试件表面应变时，试件因变形太大而匹配不出来的位置，应变应设置为0，做特殊处理；对比出现0的状态前是不是整个过程的应变的最大值，若是最大值，则此位置即为破裂位置，若不是，则将0之前的所有状态再进行正常处理，从而得到失效位置。

图12 柔性复合薄膜应变-状态图Fig. 12 Strain-state diagram of flexible composite film

由力传感器输出找出力值的最大位置，判断由应变确定的失效位置是否在由力值确定的位置的5个状态范围内（传感器具有延时性），从而判断失效位置，如果不在此范围内，则说明此次试验不规范，剔除该组试验数据。

4.2.2 截线位置确定

由国标ISO12004可知：为了成形极限测定的可重复性，裂纹两边截线的长度均不能小于20 mm，两边测量点均不能少于10个，裂纹中间大约长10 mm的虚线确定了裂纹的起始位置。

对于不同形状的薄膜试件，截线的位置确定的要求不同：当最大次应变≥0时，截线应尽可能垂直于裂纹（在±25°内）；当最大次应变<0时（窄试件），截线应平行于薄膜试件平行部方向，这样可以避免由于裂纹宽度大造成的截线左右两边的平行位置偏差。

由XJTUDIC三维数字应变测量软件计算发现：复合薄膜试件上的最大主应变方向与试件裂纹平行，由试件上的最大次应变最大值的位置和其方向即可确定一条截线的位置，由第一条截线穿过裂纹中心，另外在截线两边间隔约2 mm分别画出一条截线，如图13所示。

图13 截线位置确定Fig. 13 Determination of cut-off position

4.2.3 颈缩区域位置确定

将XJTUDIC三维数字应变测量软件计算出的薄膜试件截线上各点的应变数据用于二次曲线f(x)=ax2+bx+c拟合，确定薄膜试件的发生颈缩时的最大主应变，就是薄膜试件的实际颈缩位置，如图14(a)（彩图见期刊电子版）所示。

4.2.4 最佳拟合区域确定

由最佳拟合区域算法可以确定颈缩区域两端的左右拐点位置以及最佳拟合区域的宽度ω，从而确定了薄膜裂纹左右两部分最佳拟合区域内各点的应变值，如图14(b)（彩图见期刊电子版）所 示。

图14 试件颈缩区域Fig. 14 Necking area of specimen

4.2.5 极限应变拟合

依据极限应变拟合算法，对4.2.4中确定最佳拟合区域的各点的应变值进行反二次曲线拟合，拟合出截线上的最大主应变 ε1和相对应的最大次应变 ε2，从而确定该试件上截线的极限应变对(ε1,ε2)，如图15所示。

图15 极限应变拟合Fig. 15 Limit strain fitting

4.2.6 极限应变曲线拟合

将同一柔性复合薄膜试件上或相同形状薄膜试件上截线的极限应变进行归类合并，分别算出最大主应变和最大次应变的平均值，得到柔性薄膜试件的极限应变对 (ε1,ε2)，将 ε1作 为Y轴，ε2作为X轴，按照X轴从小到大的顺序，将各个极限应变点连起来，则连起来的折线就是该种材料的极限应变曲线，如图16所示。

图16 柔性复合材料的极限应变曲线Fig. 16 Limit strain curve of pliable composite film

4.2.7 讨论

在柔性复合薄膜进行裂纹位置拟合时，由XJTUDIC三维数字应变测量软件计算薄膜表面应变，会遇到裂纹位置太偏（正常情况下，裂纹在试件的中间位置）或是出现多个裂纹的情况，如图17所示，这是因为薄膜试件本身制作问题或者是在做柔性复合薄膜成形试验时不规范，如没有润滑等，致使柔性复合薄膜试件受力不均匀而导致裂纹位置太偏，将拟合出的裂纹位置偏离试件中间15%的范围和多裂纹试件的数据剔除。

图17 裂纹位置太偏或多裂纹Fig. 17 Crack position is too biased or multiple cracks

在柔性复合薄膜进行极限应变拟合，进行反二次函数拟合时拟合出的曲线会发生突变，如图18(a)所示。由于最大次应变太小，接近于X轴，进行反二次拟合，将最大次应变取倒数，变成无穷大，拟合出的曲线就会发生突变，所以不能利用反二次函数进行拟合。根据数据在产生突变位置的特点，改成进行线性拟合，结果如图18(b)所示。

5 结 论

本文提出一种基于双目立体视觉结合数字图像相关法的成形极限曲线测量方法，通过试验验证其可行性，极限应变测定精度能够提高0.02%；搭建出用于材料成形极限应变测量的软件，为测定柔性薄膜材料成形极限应变曲线提供了一种高可靠性、高精度手段。

本文利用光学测量的方法，采用无接触的数字图像相关法对薄膜材料进行极限应变测量，克服了传统应变方法的缺点，具有精度高，操作方便等优点。

该方法不仅能测量出材料最终的成形状态，还能够快速地确定各个状态时的材料表面应变分布，可一次性获得大量的试验数据。

对于薄膜材料成形过程中产生大变形或裂纹时难以匹配的问题，根据系列图像相邻状态变形的连续性，提出了一种图像匹配基准自适应更新的弱相关分步匹配方法。

本文方法不仅可以测量刚性很大的金属材料的成形极限曲线，还可以测定刚性很小的柔性薄膜材料的成形极限曲线。