基于立体视觉的板料成形极限应变测量.pdf

第3 1 卷第1 0 期航空学报 V 0 1 3 1N o 1 0 2 0 1 0 年1 0 月 A C T AA E R O N A U T I C AE TA S T R ( ) N A U T l C AS I N I C AO c t 2 0 1 0 文章编号：1 0 0 0 6 8 9 3 ( 2 0 1 0 ) 1 0 2 0 9 3 1 0 基于立体视觉的板料成形极限应变测量 关键技术及其系统 叶南，张丽艳 ( 南京航空航天大学航空宇航制造工程系。江苏南京 2 1 0 0 1 6 ) K e yT e c h n i q u e sa n dS y s t e mo fS h e e tM e t a lF o r m i n gL i m i tS t r a i n M e a s u r e m e n tB a s e do nS t e r e oV i s i o n Y eN a n ，Z h a n gL i y a n ( A e r o s p a c eM a n u f a c t u r eE n g i n e e r i n g ，N a m i n gU n i v e r s i t yo f A e r o n a u t i c sa n dA s t r o n a u t i c s ，N a n ji n g21 0 016 ，C h i n a ) 摘要：设计并实现了一种基于双目立体视觉和数字图像相关方法的板料成形极限应变测量系统- - B O S A S ( 双目视觉应变测量分析系统) ，对其中的关键技术作了深入讨论。该系统可以克服传统坐标网格方法很难进 行变形过程的动态监控、自动化程度不高等局限，并可计算试件各变形阶段全场应变分布，以及重建不同变形时 刻下试件的几何外形。将分段位移传递法和有限应变理论相结合，可计算大变形下的极限应变。对深冲铝板 6 0 1 6 以及航空铝板2 A 1 2 一T 4 、2 A 1 2 一O 、7 1 3 0 4 一O 等材料的应变测量结果显示本文方法适用于板料成形三维变形 及极限应变的测量。并将测得数据与由坐标网格方法测得的数据进行比较，结果表明本文方法测量准确、可靠。 关键词：立体视觉；成形极限图；应变；数字图像相关，大变形 中图分类号：V 2 6 1 2 + 8 ；T P 3 9 1 4 1文献标识码：A A b s t r a c t ：As h e e tm e t a lf o r m i n gl i m i ts t r a i nm e a s u r e m e n ts y s t e m B ( ) S A S ( b i n o c u l a rs t r a i na n a l y s i sm e a s u r e m e n ts y s t e m ) i sd e s i g n e da n di m p l e m e n t e d ，w h i c hi sb a s e dO nb i n o c u l a rs t e r e ov i s i o na n dd i g i t a li m a g ec o r r e l a t i o nm e t h o d T h ek e yt e c h n i q u e sc o n c e r n e da r ed i s c u s s e di nd e t a i l T h ep r o p o s e dm e t h o dc a nb r e a kt h el i m i t a t i o n so ft h ec l a s s i c a lc o o r d i n a t eg r i da n a l y s i sm e t h o d ，s u c ha sd i f f i c u l td y n a m i ct r a c k i n go fd e f o r m a t i o np r o c e s s a n d l o wl e v e la u t o m a t i o n ，a n dc a nc a l c u l a t ew h o l e - f i e l ds t r a i nd i s t r i b u t i o na n dr e c o n s t r u c tt h r e e - d i m e n s i o n a l ( 3 D ) s h a p eo ft h es p e c i m e n ss u r f a c ei ne v e r ys t a g ed u r i n gt h ew h o l ed e f o r m a t i o np r o c e s s L i m i ts t r a i ni sa b l e t Ob ec a l c u l a t e db yt h ec o m b i n a t i o no ft h ef r a c t i o n i z e dd i s p l a c e m e n tt r a n s f e r ( F D T ) a l g o r i t h ma n dt h ef i n i t e s t r a i nt h e o r y S t r a i nm e a s u r e m e n tr e s u l t so fm a t e r i a lv a r i e t i e ss u c ha sd e e pd r a w i n ga l u m i n u ms h e e t6 0 1 6 ，a n d a e r o n a u t i c a la l u m i n u ms h e e t2 A 1 2 - T 4 。2 A 12 一Oa n d7 8 0 4 一Od e m o n s t r a t et h a tt h ep r o p o s e dm e t h o di sa p p l i c a b l e t ot h em e a s u r e m e n to f3 Dd e f o r m a t i o na n dl i m i ts t r a i no ft h es h e e tm e t a lf o r m i n g C o m p a r i s o no fs t r a i nd a t a o b t a i n e db yB O S A Sa n db yc o o r d i n a t eg r i da n a l y s i sm e t h o ds h o w st h a tt h es u g g e s t e dm e t h o di sa c c u r a t ea n d r e l i a b l e K e yw o r d s ：s t e r e ov i s i o n ；f o r m i n gl i m i td i a g r a m ；s t r a i n ；d i g i t a li m a g ec o r r e l a t i o n ；l a r g ed e f o r m a t i o n 板料成形是材料加T 技术中的重要分支之 一，在航天、航空、船舶、汽车、电器等国民经济的 各部门得到了广泛的应用。以极限应变构成的成 形极限图( F o r m i n gL i m i tD i a g r a m ，F L D ) 卜2 0 可 以为评价板料成形性能以及改进成形工艺提供技 术基础和实用判据【3 J 。 确定F L D 的方法通常有以下几种 4 3 ：根据 塑性理论计算得到；通过有限元数值模拟来分 析板料的成形极限；在实验室条件下由标准实 收稿日期：2 0 0 9 1 0 - 3 0 ；修订日期：2 0 1 0 0 4 0 5 基金项目：国家自然科学基金( 5 0 8 7 5 1 3 0 ) ；教育部博士点基金 ( 2 0 0 8 0 2 8 7 0 0 1 6 ) ；江苏省科技支撑计划( B E 2 0 0 8 1 3 6 ) 通讯作者：张丽艳E m a i l ：z h a n 9 1 ) 7 n u a a e d u c n 验测定。根据塑性理论计算成形极限图主要是通 过采用不同的屈服准则和塑性本构关系，利用不 同的拉伸失稳准则作为判断发生颈缩与破裂的条 件来进行解析的。但每种准则适用范围有限，而 且在对具有复杂应变路径、复杂加载历史的板料 成形极限预测上与实验结果还有一定差距；将韧 性断裂准则与有限元数值计算相结合模拟材料成 形极限方法的计算精度受到诸多因素的影响，如 何恰当地进行单元的划分、单元形函数的选取、迭 代参数的选定和材料模型的选择，以及如何给出 合理的初值和边值条件并选取适当的韧性断裂准 则作为成形极限判据，是决定模拟精度的关键。 万方数据 航空学报 第3 1 卷 基于实验手段的成形极限图是在实验室条件 下，通过改变试件宽度和润滑条件，并采用一定的 应变测量技术直接获得极限应变数据以建立材料 的F L D 。坐标网格技术是目前研究板料成形表 面应变大小和分布的主要实验方法 5 。通常的做 法是在板料表面印制规则的方形网格或圆形阵 列，通过测量对比这些图案在变形前后形状参数， 来确定板料表面的应变。传统的应变测量方法是 利用单个圆形图案在板料产生变形时近似变成椭 圆的性质，采用工具显微镜直接测量椭圆的长轴 和短轴长度近似获取单个圆所在区域的最大、最 小主应变 3 。这种方法对操作者的依赖性较大、 效率不高、精度较低。后来，文献 6 结合数字图 像处理方法自动计算单个圆变形参数，但由于该 方法处理的是在二维平台上经光学投影后的网格 图像，所以存在一定的测量误差。文献 7 中提到 采用单数码相机对板料上的网格进行图像采集， 通过数字图像处理方法得到变形后的方形网格交 点，再由计算机视觉理论解算出网格交点在空间 的三维坐标用以求出相对位移，从而得到板料表 面应变分布。由于该方法不能进行变形过程的动 态监控，只能对成形的最终状态进行测量，而实验 中又很难保证成形最终状态试件表面刚好是出现 局部颈缩的极限状态，因此该方法只能从已破裂 的试件表面非裂纹处近似估计极限应变。另外上 述方法在图像处理中需要有较多的交互操作，这 一方面会影响自动化程度，同时也会给测量结果 的可靠性带来一定影响。 由W H P e t e r s 等 8 1 提出的数字图像相关 ( D i g i t a lI m a g eC o r r e l a t i o n ，D I C ) 方法经过2 0 多 年的发展，已广泛应用于固体力学领域中涉及物 体表面变形和应变测量的场合 9 。文献 1 0 3 报道 了采用立体视觉技术并结合D I C 方法进行动态 应变测量的基本方法，但论文中没有对大变形条 件下应变的测量等关键问题进行讨论，也没有提 及极限应变的计算。文献 1 1 尝试将相关方法应 用于成形极限应变的测定，并利用应变速率的变 化来判断颈缩的出现，但该方法仅使用单个相机 采集图像，只局限于获取板料各点的面内位移。 在对测量系统建模理论、高精度相机标定、立 体图像及时序图像匹配、三维重建、大变形条件下 位移跟踪和极限应变计算等关键问题进行系统研 究和综合集成的基础上，笔者所在课题组自主研发 了基于双目立体视觉的应变测量分析系统( B O - S A S ) ，以提供板料成形动态全场应变和极限应变 的测量手段，满足航空制造中板料成形研究的需 要。B O S A S 系统目前已应用于航空铝板2 A 1 2 - T 4 、2 A 1 2 一O 、7 1 3 0 4 一O 等航空材料极限应变的测定， 并已进一步应用于相应板料成形工艺的优化。本 文将对其中的关键技术及系统实现方法进行讨论。 1 双目立体视觉系统模型 双目立体视觉是一种通过从两个不同视角观 察同一场景以获得其三维描述的技术。图1 是双 目立体视觉原理的简单示意图，其中D 是待测量 的三维点，J l 和J r 分别是点D 在左、右相机像平 面上的投影点。o l 和0 r 分别是左、右两个相机 的光心。在下面两个前提成立的条件下，就可以 利用三角法的原理，从一对图像上( 称为立体图像 对) 解算出实际空间点的三维坐标： ( 1 ) 双目立体视觉装置的系统参数已知，包 括两个相机的相对位置参数以及相机的焦距、主 点等内部参数。这些参数一般可以通过离线系统 标定解决。 ( 2 ) 能够正确地建立图像点J l 和J ，之间的 对应，也就是说对应同一空间点D 在左右两幅图 像上的像点能被正确识别和匹配。这个过程也被 称为立体匹配。 图i 双目立体视觉模型示意图 F i g 1 S c h e m a t i co ts t e r e ov i s i o nm o d e l 1 1 双摄像机成像模型 在理想的针孔成像模型下，三维空间点D 在 左、右相机成像平面上的投影可表示为 即江吼，h 门 陋I ，r 0 A l r210 风， 【oo X 。 y - Z 。 1 ( 1 ) 万方数据 第1 0 期叶南等：基于立体视觉的板料成形极限应变测量关键技术及其系统 式中：下标l 和r 分别表示左相机和右相机( 下 同) ；5 为比例系数；( X ，y 。，Z 。) 为点D 在世界坐 标系下的坐标；( “，铆) 为点D 在左右图像平面上 的像点转换成像素表示的坐标；R 和f 分别为从 世界坐标系到左、右相机坐标系的旋转变换矩阵 和平移变换向量；A 为相机的固有成像参数矩阵； ( 扩，V o ) 为图像平面的主点坐标；口、口为图像在“ 轴和口轴的等效焦距参数。 由于镜头畸变的影响，实际成像点J a 并不严 格是三维点D 和光学中心o 。的连线与图像平面 的交点J u ，而是有了一定的镜头畸变偏移，如图2 所示( 这里f 。可以代表图1 中的J I 和J r ) 。镜头 畸变呈复杂的非线性，一般可用高阶多项式模型 近似。B O S A S 系统综合考虑了径向畸变和切向 畸变，将畸变模型表示为 蝴= “。+ 比。( 五1 ，+ k 2 4 ) +1 z a l l t a 仇+ d 2 ( ，+ 2 u ：) I V d = 仇+ v 。( 志l ，工+ 愚2 ，) + d 。( 产+ 2 诟) + 2 d ：U u V 。 J 式中：下标U 和d 分别表示无畸变和有畸变的情 况；( “。，钆) 和( ，钆) 分别为I 。和J d 的像素坐 标；r 一 “：+ 口：；志。、k 2 为径向畸变系数；d 。、d 2 为 切向畸变参数。式( 1 ) 中的矩阵A 和式( 2 ) 中的 畸变系数k 。、k ：、d 。、d ：统称为相机的内部参数。 图2 畸变模型示意图 F i g 2S c h e m a t i co fl e n sd i s t o r t i o nm o d e l 对于双目相机系统，除了两个相机的内参数 需要标定之外，还需要确定它们之间的相对位置 和姿态。如图1 所示，用R 畦。表示右相机坐标系 相对于左相机坐标系的旋转矩阵，t 吨表示右相机 坐标系相对于左相机坐标系的平移向量，则有如 下变换关系： D i = R r 2 l D ，+ t a l ( 3 ) 式中：D l 和D ，分别表示三维点D 在左相机和右 相机坐标系下的坐标。如果已知从世界坐标系到 左相机坐标系和右相机坐标系的旋转矩阵R 。和 R f 以及平移向量t - 和t ，则很容易推得： R r 2 l = R l 肝1( 4 ) t ，2 l = t l R l R i l t ， ( 5 ) 1 2 双相机内外参数确定 相机内外参数标定对于双目立体测量至关重 要，对此人们作了大量研究。其中基于平面模板 的标定方法在使用便捷性和标定精度方面具有综 合优势。B O S A S 系统使用圃点平面标靶，采用在 文献 1 2 、文献 1 3 基础上加以改进的三步标定 方法 1 胡进行系统标定。 三步标定方法的第1 步是利用平面标定板的 射影变换性质，分别由线性模型初步确定两个相 机的内外参数，即确定A l 、R I 、t i 和A ，、彤、t ；的值 ( 上标i 表示第i 次拍摄，下同) ，再通过式( 4 ) 、式 ( 5 ) 便可以求得两个相机的相对位置关系。 理论上说，各次拍摄时的R 乞。和如- 都应该相 等。但是，由于有噪声和计算误差的影响，而且第 1 步中没有考虑非线性畸变，因此由各个图像求 得的R ；。和如- 的结果有较大的偏差，所以方法的 第2 步是使用L e v e n b e r g - M a r q u a r d t ( L M ) 算法 对两个摄像机的内部参数和相对位置及姿态参数 进行优化。L M 算法最小化的目标函数为 。刀- ； 0 砰- i m ( A ，是，k ：，d 。- ，d 2 。，嗣，f ，Y ) I I2 + i - - IJ 喜1 f m f 0 砖一j d r ( A k ，d d z ， i _ 1j = 1 R i ，t i ，R 1 2 I ，t r 2 I ，) I I 2 ( 6 ) 式中：，z 为共摄取的图像对数；上标f 为第J 个特 征圆点；m i 、m ；分别为左、右相机第i 幅图像上识 别出的特征圆点个数；为第个特征圆点三维 坐标，可由靶点的空间相对位置予以确定；j 。- ( A - ， k l l ，k z l ，d l l ，d 2 l ，R ，t i ，驴) 、j d r ( A l ，足1 ，k 2 ，d l ，d 2 ， 冈，t i ，R 坨，t 吐I ，) 分别为标定板上的第歹个特征 圆点使用线性模型式计算出理想像点位置，再由 式( 2 ) 的非线型成像模型计算出来的在左相机和 右相机第i 幅图像上的图像坐标；开、印分别为 第_ 个特征圆点在左、右相机第i 幅图像上识别 出的真实图像坐标。本步骤中L M 算法优化的 参数包括A I 、惫l I 、志2 l 、d l I 、d 2 I 、闻、t 、A ，、惫l r 、惫2 ，、d l ，、 d 2 ，、R 心l 、t f 2 l 。其中A I 、A ，、扁、t i 的初始值由第1 步标定得到，k l l 矗2 I 、d l I 、d 2 I 和志1 ，、惫2 ，、d l ，、d 2 ，的初 始值设为0 ，Rr 2 l 和t 心的初始值分别取R ：2 - 和如的 万方数据 航空学报 第3 1 卷 平均值。这些初始值已经比较接近精确结果，因 此L M 算法可以快速收敛。 经过前面两个步骤，摄像机所有内外参数包 括非线性畸变因子都可以得到。传统两步标定法 到上一步即结束了。然而，由于平面标靶的平面 度以及各靶点的三维坐标可能存在误差，因此会 影响到系统标定结果的精度。为此在第2 步优化 结果的基础上进行第3 步的优化，即将靶点的三 维坐标也松弛为优化变量进行优化。其中内外参 数的初始值可取为第2 步标定得到的结果，各个 靶点的三维坐标D 的初始值取为靶点坐标的设 计值( 世界坐标系即建立在平面模板上) 。对式 ( 6 ) 再次进行L M 优化( 又称光束平差) ，可获得 B O S A S 系统所需的标定参数。由于本步骤中靶 点的三维坐标点也作为变量参与优化，因此，如同 自标定方法一样，标定后得到的三维坐标与实际 真实的三维坐标可能相差一个比例系数，对左右 摄像机间的平移向量进行简单的缩放即可得到实 际尺度的标定结果，具体为：t 畦l = ( L L ) t 吨l ，其 中L 为标定优化之后得到的模板上两个特征点 的间距，L 为事先精确测量该两个特征点间距的 实际值。通常的标定算法都假定靶标的三维坐标 均完全精确已知，在标定过程中作为固定不变的 量处理，因此有一定误差的标定靶标必然导致系 统标定精度的降低。三步法中的最后一步就是充 分考虑了用于标定的靶标点三维坐标可能存在的 误差，借鉴自标定的思想，在优化过程中将三维坐 标点也加入到优化参数当中，从而得到实际的三 维靶标点，消减靶标误差对标定精度的影响。这 种方法不需要标定板所有点坐标精确已知，而仅 仅需要精确已知一个标准距离即可。这不仅大大 降低了标定板的制作和计量校准要求，而且能够 显著提高系统的标定精度。 1 3 三维重建 对双目立体测量系统进行标定之后，就可以 根据左右图像上提取得到的对应同名像点求得相 应的三维点坐标。其中，左右图像的对应问题将 在第2 节讨论，这里讨论已知同名点的图像坐标 反求三维点的方法。 设图1 中左右图像上点J l 和J f 在各自的相 机坐标系中的三维物理坐标为J 。一、J 仃，下标c 表 示在相机坐标系下。记原点0 I 和J l 连成的射线 为L ：n J d ( 口R ) ；同样，o r 和J ，连成的射线R 在 左相机坐标系( 也是本文所取的测量系统坐标系) 中可表示为R ：勰f 2 l L + t f 2 I ( 6 R ) 。这里，a 和6 分别为两条射线的参数。理论上，射线L 和R 应 该相交，它们的交点即为所求的三维坐标点。但 是由于图像匹配以及数值计算等误差因素的影 响，它们一般不相交，因此本文求取连接L 和R 的公垂线的中点作为重建的三维坐标。设w 表 示与这两条射线都垂直的向量，即： ，= J 。l R 删L ，那么三角法变为求取a 、b 和C ( f R ) ， 满足： d 。I 一硪，2 I j r c ，+ 删= t r 2 l ( 7 ) 根据式( 7 ) 得到参数a 、b 和C 后，则待求的三 维坐标点为 n J 。I + c w z ( 8 ) 2 空间位移定位和应变计算 物体表面应变是由物体表面对应点在空间上 的相对位移得到的，而板料成形极限应变是属于 局部应变，需要根据物体表面局部位移计算得到。 但由于物体在变形过程中其表面形貌已经发生变 化，所谓单个点的跟踪没有明显的物理意义，如果 通过使用子区的概念，来跟踪局部区域的变形程 度，并以相邻子区中心连成的三角片在变形前后 空间上的相对变形来量化局部区域的变形程度可 以近似表达物体表面的局部变形。 数字图像相关是一种基于局部子区的相关方 法，它利用具有随机纹理的被测物体表面在数字 图像上会形成随机灰度分布以致在局部区域内可 以形成某种统计意义上的独一性的特点，以进行 立体匹配和时序图像上的像素跟踪。 2 1 数字图像相关方法 基本原理如下：在参考图像( 可称为待匹配图 像) 上取一待求点( “，口) 为中心的( 2 M + 1 ) 像素 ( 2 M + 1 ) 像素的矩形参考图像子区，在目标图像 ( - - I 称为匹配图像) 寻找与选定参考图像子区相似 程度最大的的以( 乱7 ，口7 ) 为中心的( 2 M + 1 ) 像 素( 2 M + 1 ) 像素矩形区域，这样就形成了一对 匹配点( 甜，口) 和( “，7 3 7 ) 。由于实际物体表面的局 部变形，反映在图像上，即为图像子区的变形。为 了提高匹配精度，需要量化目标图像子区的变形， 并采用一阶或二阶位移模式( 又称形函数 1 5 ) 来 近似表达子区的变形。图3 给出了图像子区中的 点在变形前后的映射关系。其中点P ( U 。，v o ) 为 参考图像子区的中心，点P 7 ( 7 。，。) 为P 点在变 形图像子区中的对应点，并且满足越7 。= U ( p ) ， 万方数据 第1 0 期叶南等；基于立体视觉的板科成形极限应变测量关键技术及其系统 2 0 9 7 口7 。- - - V ( p ) ，其中U ( ) ，y ( - ) 表示p 到P7 的映 射关系( 又称位移模式函数) 。点口( “，D ) 为p 点 附近的一点，且满足U = + A u ，= 饥+ A v 。根 据变形的连续性，口点在变形子区的对应点 q7 ( “7 ，v 7 ) 也应该在P 7 附近，并且满足“7 = U ( 口) ， 口7 = V ( 口) 。若使用一阶位移模式函数来表示 【，( ) ，V ( ) ，则有： U t = U O + A “+ “。+ 警比+ 豢u ? 甜? “ ( 9 ) 口7 = “ 0 0 + A v + 弘+ 婺越+ 警口I a 研a 可J 式中：U ，m 分别为参考图像子区中心在，口方向 上的位移值；A u 、A v 为点( “，口) 到计算窗口中心 ( “o ，仇) 的距离；1 3 u 。a “ 、a “。a 口、B y , B u 、a 翟。a 口 为图像子区的位移梯度。 图3 图像子区中的点在变形前后映射示意图 F i g 3 S c h e m a t i co fs u b s e tm a p p i n gb e f o r ea n da f t e r d e f o r m a t i o n 本文使用零均值归一化互相关函数1 6 3 来量 化子区间的相似程度： c ，。( s ) = M f I ( ，( 剐) 一7 ) ( g ( “7 ，口7 ) 一；) I ! 三二丝! = = 丝： 厂田1 c r 一厂丽玎一 ( ，( 剐) 寸y ( g ( “7 ，) 百) 2 VH = M 俨：一 fY a 一 f口I M ( 1 0 ) 当函数值为0 时，表示两个子区相似程度最 低；当函数值为1 时，相似程度最高。式中：s 为位 移模式中的待求参数向量，对一阶位移模式，有s = “ t a t t S ua 牡。3 v。a 口。a ua 砂，a 印 ； f ( u ，口) ，g ( u 7 ，可7 ) 分别为参考图像子区和目标图 像子区各点的灰度值；7 ，孑分别为参考图像子区 和目标图像子区的灰度平均值。对式( 1 0 ) 使用 N e w t o n - R a p h s o n ( N - R ) 优化进行迭代求解 17 ，并 采用近似H e s s a i n 矩阵的方法 1 8 3 减少N R 优化 过程的复杂度，最终可得亚像素级别的最优匹 配点。 由于用式( 9 ) 计算出的变形子区的像素坐标 一般都不是整数，因此需要对目标子区图像进行 插值以获得亚像素坐标位置的灰度值。而插值方 法的选择将会很大程度上影响N R 优化的效率。 本文采用双三次样条基插值表示以点( U7 ，秽7 ) 为 中心的N N 像素的子区灰度分布曲面。实验 表明N = 6 时，可在时间复杂度和优化效率的平 衡上达到最优。 2 2 空间位移定位 如图4 所示，在初始参考阶段( 在未变形时刻 对应的阶段) 左图像中设定待计算的图像子区，并 通过相关方法找到其在同一阶段右图像中的相对 应的最优匹配子区，再根据1 - 3 节可重建得到子 区中心点对应的空间位置。同理，在变形阶段左、 右图像中通过相关方法也追踪到其相应位置，则 同样可得到子区中心点变形后的空间位置。 图4 立体和时序匹配示意图 F i g 。4 S c h e m a t i co fs t e r e oa n dt e m p o r a lm a t c h i n g 在板料极限应变的测量中，物体表面经常出 现大变形情况，而常用的数字相关方法因参考图 像和目标图像间变形过大造成相关性减低，不能 精确匹配子区中心。因此提出一个建立在数字图 像相关原理基础上的分段位移传递方法进行变形 阶段的立体匹配以及空间位移定位【l 引。其基本思 想是将两个阶段的变形量分散至这两个阶段之间 的若干阶段序列中，再利用序列图像进行相关匹配 获取最终的子区中心匹配。通过控制电荷耦合器 件( C C D ) 的采集频率，可以保证相邻序列图像间有 万方数据 2 0 9 8 航空学报 第3 1 卷 较高的相关性，运用分段位移方法可以保证最后的 位移定位精度。具体过程如图5 所示，如果初始参 考子区与变形子区工进行了正确的相关匹配，但因 变形过大而无法与变形子区进行相关匹配，则将 变形子区I 设置为新参考子区。由于与初始参考子 区一点p ( u ，。) 相匹配的点p 7 ( “7 ，) 一般为亚像素 坐标，因此首先将点p 圆整到整像素位置，即点 而( 么，磊) 。通过D I C 方法，得到矗在变形子区上 的对应点硝( 名，稿) ，同时得到位移模式中的参数 向量s 。根据变形的连续性，可使用式( 9 ) 位移模式 函数得到P 7 的匹配点矿，最终可得点p 与点，之 间的像素位移“。= 一甜，V t = 一掣。值得说明的 是，该方法同时适用于待匹配点为亚像素坐标的相 关匹配，如图4 中右图像时序匹配。 变形子区 新参考子区 图5 分段位移传递法( F D T ) 示意图 F i g 5 S c h e m a t i co ff r a c t i o n i z e dd i s p l a c e m e n tt r a n s f e r ( F D T ) 2 3 应变计算 由于板料成形经历的是塑性变形，因此需要 用有限应变理论瞳o 来计算物体表面的应变分布。 文献 7 提到了在局部各向同性和网格内部应变 均匀的假设下，利用柯西一格林张量与变形前后网 格边长的拉形比关系以及柯西一格林张量与拉格 朗日应变的关系可计算网格内部的真实应变，但 其前提是初始三角网格在进行局部坐标对齐后必 须有两边正交。虽然本文在参考阶段左图像上生 成一系列已知拓扑关系的虚拟三角网格，而且在 图像上这些网格有两边是正交的( 如图6 ( a ) 所 示) ，但射影变换并不是保角性的，重建后的空间 三角网格并不能保证有两边正交。因此本文采用 面内形变梯度张量和柯西一格林张量 1 8 以及柯 西一格林张量与拉格朗日应变的关系计算网格内 部的真实应变，具体过程如下。 图6 ( b ) 为变形阶段虚拟网格变形示意图。 已知变形前后相对应的空间三角网格各顶点坐 标，如图7 ( a ) 所示。定义面内形变梯度张量为 F ，它是二阶张量，则根据图7 ( b ) 所示的坐标对应 关系，有 z ：= F 1 12 C l + F 1 2 y 1 ，y ：= F z l - r I + F 2 2 Y l ，1 、 z ：= F 1 1 2 2 + F 1 2 y 2 ，y ：= F 2 1 2 2 + F z 2 y 2J 图6 虚拟三角网格示意图 F i g 6 S c h e m a t i co fv i r t u a lt r i a n g u l a rm e s h e s ( 五，K ，互) ( 掣，球z p ( ”啦标变换后变形前后空间三角网格 图7 应变计算示意图 F i g 7S c h e m a t i co fs t r a i nc a l c u l a t i o n 事实上，F 可进一步分解为：F - - - - - R U ，其中R 为刚体旋转矩阵，【，为形变张量，因此F 并不具有 对刚体运动的不变性。而柯西一格林张量可表示 为：c = F T F ，可见C 对刚体运动具有不变性。因 此，采用对空间三角形网格进行坐标变换的方式消 除刚体运动分量，则图7 ( b ) 中的坐标计算为 z l I ，lI ，Y l = 0 ，z ：= I1 ，：I ，Y ：= 01 z 2 = ( ，l 1 ，2 ) I ，lI ，Y 2 = I ，l 屹I I1 ，lI z ：= ( ，。耽) l ，：I ，y ：= M 1 ，：I lJ ( 1 2 ) 式中：y 。、 ，：、耽、 ，：分别为三角形各边所对应的 向量。 万方数据 第1 0 期叶南等：基于屯体视觉的板料成形极限应变测量关键技术及其系统 而柯西一格林张量与拉格朗日应变的关系为 C - - J 4 - 2 E ，其中E 为拉格朗日应变张量，J 为单 位矩阵。定义： E m 。x m i n = E 。，+ E 2 。以瓦1 i F F 囹瓦丁 2 ( 1 3 ) 则最大最小真实应变可表达为 。砸= I n 14 - 2 E 。；m i 。 ( 1 4 ) 另外，B O S A S 系统使用节点应变来平滑用上 述方法所得到的应变值。节点应变可以通过共节 点的三角网格上的托格朗日应变得到。一般有3 种常用的节点应变类型，如图8 所示。 泌| 、 弋 ( a ) 权重型( b ) 局部化型 X 哌 ( c ) 对角型 图8 节点应变类型 F i g 8T y p eo fn o d es t r a i n 图8 中O 点即为节点应变点，数字表示所在的 三角网格在计算节点应变时所占的权重。本文使用 图8 ( b ) 局部化型的节点应变。具体过程为：由式 ( 1 1 ) 式( 1 3 ) 计算节点O 相邻的4 个三角网格各自 的拉格朗日应变。由于拉格朗日应变各分量与坐标 的选择无关，因此有 4 44 毋。= 日。4 ，E o ：一鹾：4 ，E o ：= 鹾z 4 I = 1I = 1 = 1 再由式( 1 4 ) 可计算节点O 处的最大、最小真 实应变值。 读入 图像 序列 测量结束 3系统硬件组成和软件系统流程 3 1 系统硬件组成 B O S A S 系统在实物上主要由两个数字摄像 机、一个光源、平面标定板、相机同步采集控制器、 电脑以及显示设备组成，如图9 所示。实验使用的 相机为D H l 3 1 0 F M 数字摄像机，光学镜头使用 M 5 0 1 8 一M P 镜头，焦距为5 0m m 。平面标定板上图 案为1 1 9 圆形阵列，相邻圆心的标称距离为 1 。5m m ，5 个大圆用于建立标定板坐标系( 亦为世 界坐标系) ，其中距离最远的2 个较大圆圆心经过 精密仪器测定距离，用于三步标定法中的比例因子 恢复( 详见1 2 节) 。在双相机和光源装置固定好 后对双目相机进行离线标定，然后启动板料成型 机，再由计算机经由控制器控制两个C C D 摄像机 同步拍摄以获取同一时刻的立体图像对和时间轴 上的图像序列。最后在软件系统中对采集到的数 据进行处理以及显示结果。 图9 系统硬件结构简图 F i g 9 S c h e m a t i co fs y s t e mh a r d w a r es t r u c t u r e 3 2 软件系统流程 软件系统流程如图1 0 所示。除需要手动拾 取一个初始点，选定图像上计算区域以及设定相 荛餮蓥娶翼绂H 碧嘉叠肇譬H 斋藩舞H 主篓并敬定计算参数Il 初始匹配点l 。l时序匹配l 宦建 仁缵澎躲嘴，H 自笺蚕存H 篙萋H 网慧嘉扑 图1 0 软件系统流程图 F i g 1 0 F l o w c h a r to fs o f t w a r es y s t e m 万方数据 2 1 0 0 航空学报第3 1 卷 应的计算参数外，从匹配到结果显示的各个过程 都可以自动完成( 计算参数包括：子区大小、各子 区中心间距、位移模式以及像点匹配成功的判据， 通常这些参数可以取系统缺省设置) 。其中灰色 框图代表了用户需要参与交互操作的步骤。 4 系统测量实例 应用B O S A S 系统已经完成了深冲铝板6 0 1 6 、 航空铝板2 A 1 2 一T 4 、2 A 1 2 一O 、7 8 0 4 一O 等材料的成形 极限测定，并已应用于相应材料板料成形件的数字 模拟分析和成形工艺优化。本文给出在相同实验 环境下针对深冲铝板6 0 1 6 ，用B ( ) S A S 系统测得的 F L D 与用坐标网格法测得F L D 的结果对比与分 析，以及航空铝板3 种材料的F L D 测量结果。 4 1 试件外形和测量现场 实验用的深冲铝板6 0 1 6 厚度为0 9m l n ，参考 I S 01 2 0 0 4 2 ：2 0 0 8 ( 即G B T2 4 1 7 1 2 2 0 0 9 ) 标准 制备试件和安排实验环境。其中试件的外形及几 何尺寸如图1 1 所示，为增加试件表面的纹理对比 度以便相关计算顺利进行并避免金属高光造成数 据缺失，可将试件表面经过去油处理后用强附着力 的哑光黑白漆交替喷涂，形成人工散斑。实验采用 N a k a j i m a 试验，实验工具剖面图如图1 2 所示。润 滑条件采用表面涂有凡士林的硅胶进行润滑。 媳 I S O 一3 0I S O 一7 0I S 0 - 9 0 I S O - 1 3 0I S O - 1 8 0 单位：m i l l 图1 1试件尺寸和几何外形示意图 F i g 11 D i m e n s i o na n dg e o m e t r yo fs p e c i m e n R 5 - 1 0纩 、 r 彤彩勿钐黝 隧惑阑蕊忒 、木乙 1 0 0 2 图1 2N a k a j i m a 试验的工具剖面 F i g 1 2 S e c t i o no fN a k a j i m ae x p e r i m e n t st o o l 4 2 系统输入和输出 实验在如图1 3 所示的现场环境下实时采集 不同宽度的板料变形过程的图像序列。在后续数 据处理中，对应不同宽度的图像序列分别选取若 干阶段作为软件系统输入，在计算参数设定步骤设 定参数如下：子区大小为2 5 x 2 5 ，各子区中心间距 7 个像素，一阶位移模式，像点匹配成功的判断条 件为：相关函数值大于0 8 且迭代次数小于3 0 次。 图1 4 为破裂前最后一个时刻各种宽度下试件表面 各采样点的应变分布图。结果输出后，选取破裂前 最后一个时刻的最大的1 0 个主应变点的平均值来 绘制成形极限曲线( F L C ) ，如图1 5 所示。 图1 3 实验现场硬件装置 F i g 13H a r d w a r ec o n f i g u r a t i o ni nl a b o r a t o r ys p o t 岛矗蕊00 0 50 1 00 1 5 0 2 0 0 2 5 最小k - l t t 变 图1 4 破裂前最后时刻各宽度下试件表面应变分布图 F i g 1 4S t r a i nd i s t r i b u t i o no fs p e c i m e n ss u r f a c ea t t h el a s tm o m e n tb e f o r ef r a c t u r e 斟 型 “ 二 嚼 一兰 图1 5 深冲铝板6 0 1 6 成形极限曲线( F L C ) F i g 1 5F o r m i n gl i m i tc u r v e s ( F L C s ) o fA A 6 0 1 6 丽囊- 一。僦一燮i 万方数据 第1 0 期叶南等：基于立体视觉的板料成形极限应变测量关键技术及其系统 为验证本文方法的可行性和准确性，在相同 的实验条件和环境下，采用坐标网格方法 7 j 测量 相同批次深冲铝板6 0 1 6 的另一条成形极限曲线， 如图1 5 所示。从图中可见两条曲线的应变分布 趋势比较接近，但坐标网格法一般在已破裂的试 件上进行测量，因此测量值比临界状态时的真实 值偏小，而本文方法可以动态跟踪到破裂前的临 界状态，因此测量值更接近真实值，从图1 5 中可 见B O S A S 系统测量曲线整体在坐标网格法测量 曲线的上方。另外，从图中可以观察到前者中各 处的主应变点相对集中，方差小，而后者受到变形 制 强 “ 嘣 前后网格质量和测量方法等的影响，数据较为 分散。为显示试件动态变形过程中的应变分布 情况，从I S ( ) 一7 0 试件变形过程中选取1 5 个阶 段，分别生成各变形时刻下试件表面三维外形 和最大主应变分布，如图1 6 所示。从图中可 见，最大主应变的极值首先出现在试件中间部 位并逐渐向两边发展，且呈现对称分布，这与实 验观察相一致。 图1 7 和图1 8 分别给出了使用B O S A S 系统 测得的2 A 1 2 - T 4 和2 A 1 2 一O 材料3 种不同厚度 ( 1 0 ，1 2 ，1 5r a m ) 的F L C 。 图1 6I S O 一7 0 不同阶段试件几何外形和表面应变分布图 F i g 16 3 Ds u r f a c ea n ds t r a i nd i s t r i b u t i o no fI S O - 7 0 w i t hd i f f e r e n td e f o r m a t i o ns t a g e s 、 - - - _ _ _ 鹣 ? 萨 一 一：， 浮 二 一2 A 1 2 一I 41 O m m 2 A 1 2 1 “ 41 2 m m 2 A 1 2 一T 41 5m m 宽煌3 0 m m 宽，2 7 0 m m 宽度9 0 m m 宽艘1 1 0 m m 。者府13 0 m m 宽馊15 0 m m 宽度1 8 0 m m 图1 73 个厚度2 A 1 2 一T 4 材料的F L C F i g 17 F L C so f2 A 1 2 一T 4m a t e r i a lw i t ht h r e ed i f f e r - e n tt h i c k n e s s e s 5结论 为时序跟踪板料成形各时刻表面三维外形和 应变分布情况并测定板料成形极限应变，综合运 用双目立体相机标定、基于图像相关的立体匹配 方法、大变形条件下的位移定位和应变计算方法 等，自主研发了基于立体视觉的板料应变测量系 统B O S A S 。该系统已用于深冲铝板6 0 1 6 ，航空 制 翻 州 堪 图1 83 个厚度2 A 1 2 一O 材料的F L C F i g 1 8F L C so f2 A 1 2 一Om a t e r i a lw i t ht h r e ed i f f e r e n t t h i c k n e s s e s 铝板2 A 1 2 一T 4 、2 A 1 2 一O 、7 8 0 4 - O 等材料的成形极 限测定以及进一步的成形工艺优化。结果表明： B O S A S 系统有以下优点： ( 1 ) 可以克服传统坐标网格方法难以进行