大变形场下高速数字图像相关并行计算研究.pdf

由于大变形场下图像去相关效应的影响，数字图像相关（Digital image correlation，DIC）始终无法完成图像间的并行计算。为了突破这一瓶颈，本文提出了一种基于 AcceleratedKAZE（AKAZE）的参考图更新方法，可在 DIC 正式计算之前完成参考图更新工作，为并行计算提供独立数据。并构建了一种图形处理器（Graphics processing unit，GPU）并行计算架构，可对所有子区独立估值，完成图像间和子区间的并行计算。最后对丁腈橡胶进行了拉伸测试，结果表明相比于传统的串行 DIC 计算方法，运用本文的并行方法速度可提升两个数量级。关键词：数字图像相关；并行计算；散斑；大变形；AcceleratedKAZE（AKAZE）中图分类号：TP391.41 文献标志码：AParallel Computing of HighSpeed DIC Under Large Deformation FieldCHEN Houchuang，MA Kun，XUE Yuxuan，MENG Zhi(Faculty of Science,Kunming University of Science and Technology,Kunming 650500,China)Abstract：Due to the effect of image decorrelation under large deformation fields，digital image correlation（DIC）has never been able to complete parallel computation between images.In order to break through this bottleneck，this paper proposes an acceleratedKAZE（AKAZE）-based reference image update method，which can complete the reference image update before DIC is officially calculated，and provide independent data for parallel computing.A graphics processing unit（GPU）parallel computing architecture is constructed，which can independently estimate all subsizes and complete the parallel computation between images and subsizes.Finally，tensile tests are performed on the nitrile butadiene rubber（NBR），and the results show that compared with the traditional serial DIC calculation method，the proposed parallel method can be increased by two orders of magnitude.Key words:digital image correlation;parallel computation;speckle;large deformation;acceleratedKAZE（AKAZE)引 言数字图像相关（Digital image correlation，DIC）13是一种非接触式光学测量方法，通过电荷耦合元件（Charge coupled device，CCD）相机记录试件变形前后信息，可获得其位移、变形等力学性能。DIC 技术装置简单，具有很强的可重复分析的特点，在科学研究和工程实践中得到了广泛应用。收稿日期：20220503；修订日期：20220628陈厚创 等：大变形场下高速数字图像相关并行计算研究随着高速摄像技术的发展，DIC 方法越来越多地用于碰撞4和材料的瞬态变形测量5。这些实验一般都以大于 100 frame/s频率采集图像，在一些高速运动中甚至超过 1 000 frame/s。在此过程中会记录下大量的图像，这些图像记录了材料变形直至破坏的全过程。由于图像的数据量大，导致了相关计算和分析非常耗时，所以计算速度的提升有利于数据分析的时效性。邵新星等6提出了一种多种子点扩散并行计算方法用于土木的准静态测量。文献 7 提出了一种基于 CPU 与图形处理器（Graphics processing unit，GPU）的异构框架，实现了子区内和子区间的并行计算。但是到目前为止，DIC并行算法的研究大部分都集中于子区间并行，未能达到图像间并行的原因是由于在碰撞、快速变形的情况下，试件会出现区域变形较大的情况，导致 DIC方法出现去相关现象，造成较大的测量误差。根据现状，本文基于参考图更新方法和 GPU 并行计算能力，实现了大变形场下高速 DIC 并行计算测量。提出了一种基于 AcceleratedKAZE（AKAZE）的参考图更新方法，可在 DIC 正式计算之前完成参考图更新和分组工作，为图像间并行计算提供独立数据；并构建了 GPU 并行计算架构，可独立获取各子区初值，完成多张变形图的并行计算，提升检测效率。1 基本原理 DIC基本原理是在参考图中定义一个大小为(2M+1)像素(2M+1)像素的参考子区，由相关搜索算法在变形图中找到相对应的目标子区。计算目标子区中心点(x0，y0)和参考子区中心点(x0，y0)的差值初步得到子区的位移信息。再通过优化函数得到更精确的位移分量u，v。计算感兴趣区域（Region of interest，ROI）内的全部感兴趣点就可得到全场变形参数。1.1AKAZE特征点检测算法为了提高特征点的检测速度，文献 89 在 2013年提出了 AKAZE 特征点检测算法，该算法实现了快速显示扩散（Fast explicit diffusion，FED）10来构建非线性尺度空间。在对特征点描述上采用局部差分二值化（Modifiedlocal difference binary，MLDB）11构建信息描述子，使特征点具备旋转不变性且占用更少的存储空间。不同于传统的基于线性高斯金字塔的尺度不变特征变换（Scaleinvariant feature transform，SIFT）12和加速健壮特征（Speeded up robust features，SURF）13，AKAZE可局部自适应小细节，保留目标边界，让尺度空间具备更多的特征信息，提升算法精度9。1.2参考图更新方法对“大变形”的定义是传统 DIC 计算时产生去相关而误匹配的变形，量化来看是参考图和变形图变形量超过 20%或者旋转角度大于 7的变形。更新参考图可降低大变形引起的去相关效应，这里使用基于 AKAZE的参考图更新方法。首先提取参考图和变形图的 3个最强匹配点，再计算以这 3对匹配点为中心构建的子区相关系数大小，来判断是否需要更新参考图。为了减少参考图和变形图的计算次数，采用区间二分法对采集到的所有序列图像 N进行分组计算，具体步骤如下：（1）以 12为单位划分为k个区间。（2）令Dki=In=1，2，N，i 1，12，k 1，N12+1，Dki为第k个区间的第i张变形图，用二分法选取变形图和参考图进行 AKAZE特征点检测，构建待测子区SR（参考图子区），Ski（变形图子区）。（3）定义函数F(Dki)=C(IR，I12(k-1)+i)（1）式中：C(IR，I12(k-1)+i)为参考图IR与第k个区间的第i张变形图I12(k-1)+i待测子区的归一化互相关（Normalized crosscorrelation，NCC）系数。979数据采集与处理 Journal of Data Acquisition and Processing Vol.38,No.4,2023若F(Dki)CTh且F(Dki+1)CTh（CTh为设定阈值），则确定第 1 分组G1 I1，I2，I12(k-1)+i。令I12(k-1)+i为第2组参考图，以I12(k-1)+i+1为第1张变形图重新划分区间，重复步骤（13），直至分组完毕。1.3基于 GPU的 DIC并行计算GPU 具有数量众多的计算逻辑单元，通过计算机统一设备架构（Computer unified device architecture，CUDA）编程平台对其调用可以完成复杂的并行计算任务。如图 1为大变形场的 DIC计算架构，共分为两部分，第 1部分为参考图变换与分组，第 2部分为 2DDIC 计算。经过参考图变换后，选定一组图像进行 2DDIC 多图像间并行计算，计算包括预计算部分和迭代部分。预计算部分主要在 CPU 内完成，包括子区的初值估计，各像素空间位置索引，Hessian 矩阵和图像梯度的计算。迭代部分主要在GPU 内完成，包括形函数计算，双三次插值，迭代参量p的计算等，通过调用线程模块，开启子区间和图像间的并行计算模式。图中 IiR（i 1，n）为每个分组的参考图，IijD、IijDsubset（i 1，n，j1，N）分别为为每个分组的变形图和子区，xn subset为变形图中对应参考图的亚像素位置，Gsubset为亚像素位置的像素值，pij和 pijconverged分别为各子区的迭代参量和收敛值，u和 v分别为收敛后获得的 x，y方向的位移信息。图 1大变形场 2D-DIC计算架构Fig.1Large deformation field 2D-DIC computing architecture980陈厚创 等：大变形场下高速数字图像相关并行计算研究1.3.1初值估计为了打破子区间初值传播的依赖性，文献 7 用 3层金字塔对子区进行初值估计，该方法虽然能够独立估计出每个子区的初值，但是需要连续的 3次 DIC计算才能得到两张图像的变形信息，花费了过多的时间成本。基于此，这里提出一种基于 AKAZE 的全局子区估值方法。如图 2所示，第 1步对参考图和变形图的 ROI区域进行 AKAZE特征点检测，得到可两两匹配的特征点对，第 2步对 ROI区域进行网格划分，网格大小与子区大小相同，第 3步对落在每一个网格内的特征点求位移均值u、v，得到子区初值u=(xn-xn)/n（2）v=(yn-yn)/n（3）式中：xn、yn和xn、yn分别为参考图与变形图相对应的特征点；n为单个子区内特征点存在个数。该方法可减少过多的冗余计算，还能得到更精确的初值估计。1.3.2并行计算形函数本文采用一阶形函数14来跟踪试件 6个自由度的空间位置，从 Algorithm 1可看出，程序执行了 2个数据无关的操作，将一阶形函数作用在参考图的独立子区，得到一个新的坐标值，即变形图中对应子区位置。使用 GPU 计算时，传入形函数参数p和参考图的像素位置(x，y)。启动 Ntg个并行线程对不同变形图多子区进行并行计算。其中，N 为子区像素个数，t为 ROI包含的子区个数，g为同时计算的变形图张数。Algorithm 1 Shape function computation（1）Input：x，y for Ntg pixels in ROI，p of size 6tg（2）for Ntg parallel threads do（3）tid=getThreadId()（4）xn tid=p0+(1+p2)x tid+p3y tid（5）yn tid=p1+p4x tid+(1+p5)y tid（6）cuda.synchronize（）（7）end for图 3显示了 3种不同架构的形函数计算过程，图 3（a）为传统串行方法，每个子区逐个计算，计算次数为 tg次；图 3（b）为 Mullai的子区间并行方法，每张图的所有子区同时计算，计算次数 g次；图 3（c）为本文图像间与子区间并行方法，g张变形图所有子区同时计算，计算次数 1次。1.3.3双三次插值这里采用双三次插值法，通过对待测点周围最近的 16个采样点加权平均得到插值结果f(x，y)=i=03j=03f(xi，yj)W(x-xi)W(y-yj)（4）式中：W(x)为权重函数，如 Algorithm 2 所示，输入形函数的计算结果和双三次插值权重的查找表，输出亚像素插值结果。图 2各子区初值估计Fig.2Estimated initial values for each subsize981数据采集与处理 Journal of Data Acquisition and Processing Vol.38,No.4,2023Algorithm 2 Bicubic interpolation computation（1）Input：xn，yn for Ntg pixels，lookuptable（2）Output：G of size Ntg（3）for Ntg parallel threads do（4）tid=getThreadId()（5）ai tid=BasisWeight(lookup_table)（6）bi tid=Get(xn tid，yn tid)（7）G tid=i=116bi tid ai tid（8）cuda.synchronize（）（9）end for1.3.4并行计算迭代参量p本文采用反向组合高斯牛顿（Inverse compositional GaussNewton，ICGN）算法15对迭代参量进行优化，其形式为CZNCC(p)=f(x+w(；p)-f f-g(x+w(；p)-g g2（5）为方便并行计算，对函数优化如下：（1）计算参考图和变形图各子区像素均值，开启线程数tid=t g。-F tid=n=1NFnN（6）-G tid=n=1NGnN（7）（2）计算参考图和变形图各像素值与所在子区像素均值的强度差，开启线程数tid=N t g。ftemp tid=F tid-F t（8）gtemp tid=G tid-G t（9）（3）计算f和g，开启线程数tid=g。f tid=n=1N tf2temp（10）图 3形函数计算过程Fig.3Shape function calculation process982陈厚创 等：大变形场下高速数字图像相关并行计算研究g tid=n=1N tg2temp（11）（4）计算乘积分量Qtemp，开启线程数tid=N t g，t=tid/N，t=tid/(N t)。Qtemp tid=fwp(F tid-F t)-(f t/g t)(G tid-G t)（12）（5）计算乘积累加项Q，开启线程数tid=t g。Q tid=nNQtemp n（13）（6）计算迭代增量p，更新迭代参量p，开启线程数tid=t g。p tid=H-1 tid Q tid （14）p tid p tid p tid-1（15）以上每个分部计算完成都会执行一个线程同步操作，来确保当前所有线程计算完毕之后再开启下一个任务。2 分析与讨论 2.1参考图更新仿真测试由计算机生成的数字散斑可以有效地控制变形信息，便于评估方法性能。本小节根据文献 16 的模拟方法生成一系列拉伸变形的仿真散斑。实验测试了两组图像，尺寸大小分别为512像素 512像素和1 000像素 500像素。实验结果如图 4所示，随着变形图计算序列的增加，3个测试子区的CNCC系数逐渐降低，到第 10张变形图时，相关系数均低于所设阈值 0.9，说明去相关效应严重，匹配结果不佳。对参考图进行更新，将第 9 张变形图作为新的参考图与后续的变形图继续比较。更新参考图后，第 10张变形图与新参考图的CNCC相关系数再次增大到 0.99 范围，说明匹配可靠，保真度高。基于 AKAZE特征点检测的参考图变换，能够迅速精准地找出参考图和变形图的匹配点，再配以区间二分法对整个序列图像分组，可以有效地减少参考图和变形图的计算次数。2.2GPU并行测试为验证所提出并行算法的准确性和鲁棒性，本文对非线性材料丁腈橡胶（Acrylonitrilebutadiene rubber，NBR）标准试件进行测量。按照 GB/T528 2009标准 硫化橡胶或热塑性橡胶拉伸应力应变性能的测定 中规定尺寸，制成宽度为 25.00 mm，长度为 100.00 mm，拉伸段长度为 20.00 mm，厚度为2.00 mm。实验采用 RGM6000万能试验机对标准试件进行拉伸实验，测试之前对试件表面喷涂白色图 43个对应待测子区CNCC系数Fig.4Three CNCC coefficients corresponding to the subsize to be tested 983数据采集与处理 Journal of Data Acquisition and Processing Vol.38,No.4,2023涂料，形成随机散斑图案。以 1 mm/s速度向下拉伸，拉伸过程中丁腈橡胶逐渐变长直至断裂。采集过程中用 CCD相机（大恒 MER5007UC，分辨率 1 944像素2 592像素）对试件以 400 ms/帧的速度拍摄记录。实验过程中始终保持室温 25。用于计算的硬件配置包括 Inter i510400 6 核 CPU，NVIDIA GeForce GTX760 GPU。软件配置包括 CUDA10 Toolkit等。通过实验采集到了一系列丁腈橡胶变形实验图，图像序列记录了试件沿 y轴单向拉伸下的变形，经过 2.2 节方法更新参考图后，得到了可连续计算的多组图像，每组图像包含了一张参考图和多张变形图。这里将第 1组第 1帧作为参考图像，第 25帧作为变形图像，ROI为840像素 315像素，子区大小为21像素 21像素，进行 DIC并行计算。图 5 给出了沿 y轴的位移分布图，随着拉力的增大，丁腈橡胶逐渐下移，并且越靠下半部分颜色越深，位移越大，符合实际拉动的预设值。为了说明本文并行计算方法的优越性，将其与相同的串行 DIC计算方法做了比较，方便起见分别称之为 paDIC 和 seDIC。沿 y方向上的位移误差如图 6所示，平均误差在2 10-3像素左右。表 1 统计了 4 张变形图的计算时间。seDIC 和 paDIC 计算总时长分别为44 721 ms和 726 ms，计算速度总体提升了 62 倍左右。从表 1 中还可看出双三次插值操作和计算迭代参量p的速度提升效果远大于形函数计算，这是由于这两部分计算复杂度较高，利用 GPU 进行并行计算效果更明显。3 结束语 本文就大变形场下 DIC 计算产生的去相关效应进行了讨论，提出了一种基于 AKAZE 的参考图更新方法，有效地解决了变形图与参考图匹配度低的问题，并打破了图像间的依赖性，为图像间的并行计算提供了独立数据。在 GPU 并行计算的支持下，评估了高精度亚像素级的 DIC 计算性能，并在商用计算机上 DIC 计算速度提升了两个数量级。值得指出的是，通过优化代码，或者在更专业的 GPU 设备上可同时计算更多的图像，仍可进一步提升计算效率。图 5随拉力增大试件表面 y方向位移场分布图Fig.5Distribution of the displacement field in the y direction on the surface of the specimen with increasing tension图 6y方向位移误差Fig.6y-direction displacement error表 1计算时间对比表Table 1Comparison of computation timems项目OverallShape function computationBicubic interpolation computationComputate pseDIC44 72117011 80527 954paDIC7269.225118.2Speedup61.618.5472.2236.5984陈厚创 等：大变形场下高速数字图像相关并行计算研究参考文献：1PAN B,LI K,TONG W.Fast,robust and accurate digital image correlation calculation without redundant computationsJ.Experimental Mechanics,2013,53(7):1277-1289.2SAKAI K,ZHANG Y L,YONEYAMA S,et al.Evaluating distribution and variation of strains near eyes in blinking using digital image correlationJ.Skin Research and Technology,2020,26(5):749-756.3XU Z,DELA C J,FTHENAKIS C,et al.A novel method to measure skin mechanical properties with three-dimensional digital image correlation J.Skin Research&Technology,2019,25(1):60-67.4王永红,但西佐,胡悦,等.基于高速数字图像相关的人车碰撞伤害实验研究J.光电子激光,2017,28(1):81-86.WANG Yonghong,DAN Xizuo,HU Yue,et al.Experimental research on human-vehicle collision injury based on high-speed digital image correlationJ.Journal of Optoelectronics Laser,2017,28(1):81-86.5申海艇,蒋招绣,王贝壳,等.基于超高速相机的数字图像相关性全场应变分析在 SHTB 实验中的应用J.爆炸与冲击,2017,37(1):15-20.SHEN Haiting,JIANG Zhaoxiu,WANG Beike,et al.Application of digital image correlation full-field strain analysis based on ultra-high-speed camera in SHTB experimentJ.Explosion and Shock Waves,2017,37(1):15-20.6邵新星,戴云彤,何小元,等.实时数字图像相关用于土木准静态实验测量J.光学学报,2015,35(10):133-141.SHAO Xinxing,DAI Yuntong,HE Xiaoyuan,et al.Real-time digital image correlation for civil quasi-static experimental measurementJ.Acta Optica Sinica,2015,35(10):133-141.7THIAGU M,SUBRAMANIAN S J,NASRE R.High-speed,two-dimensional digital image correlation algorithm using heterogeneous(CPU-GPU)frameworkJ.Strain,2020,56(3):e12342.8ALCANTARILLA P F,BARTOLI A,DAVISON A J.KAZE featuresC/Proceedings of European Conference on Computer Vision.Berlin,Germany:Springer-Verlag,2012:214227.9ALCANTARILLA P F,NUEVO J,BARTOLI A.Fast explicit diffusion for accelerated features in nonlinear scale spacesC/Proceedings of the British Machine Vision Conference.S.l.:BMVA,2013.10SCHREIER H W,SUTTON M A.Systematic errors in digital image correlation due to undermatched subset shape functionsJ.Experimental Mechanics,2002,42:303-310.11YANG X,CHENG K T.LDB:An ultra-fast feature for scalable augmented realityC/Proceedings of IEEE International Symposium on Mixed and Augmented Reality(ISMAR).Atlanta,USA:IEEE,2012:49-57.12LOWE D G.Distinctive image features from scale-invariant keypointsJ.International Journal of Computer Vision,2004,60(2):91-110.13BAY H,TUYTELAARS T,GOOL L V.SURF:Speeded up robust featuresC/Proceedings of the 9th European Conference on Computer VisionVolume Part I.S.l.:Springer-Verlag,2006:346-359.14XU X H,SU Y,CAI Y,et al.Effects of various shape functions and subset size in local deformation measurements using DICJ.Experimental Mechanics,2015,55:15751590.15GAO Y,CHENG T,SU Yong,et al.High-efficiency and high-accuracy digital image correlation for three-dimensional measurementJ.Optics and Lasers in Engineering,2015,65:73-80.16SU Y,GAO Z R,FANG Z,et al.Theoretical analysis on performance of digital speckle pattern:Uniqueness,accuracy,precision,and spatial resolutionJ.Optics Express,2019,27(16):2243922474.作者简介:陈厚创（1995），男，硕士研究生，研究方向：DIC、图像处理、光电检测、光学信息处理方面的研究，Email：。马琨（1966），通信作者，男，博士，教授，博士生导师，研究方向：光测实验力学和结构损伤识别研究工作，Email：。薛宇轩（1995），男，硕士研究生，研究方向：指纹图像识别、机器学习。孟志（1993），男，硕士研究生，研究方向：指纹图像处理、光学信息处理。(编辑：陈珺）985