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

陈厚创，马琨，薛宇轩，孟志; Chen Houchuang; Ma Kun; Xue Yuxuan; Meng Zhi

网刊加载中。。。

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

摘要

由于大变形场下图像去相关效应的影响，数字图像相关（Digital image correlation，DIC）始终无法完成图像间的并行计算。为了突破这一瓶颈，本文提出了一种基于Accelerated‑KAZE（AKAZE）的参考图更新方法，可在DIC正式计算之前完成参考图更新工作，为并行计算提供独立数据。并构建了一种图形处理器（Graphics processing unit，GPU）并行计算架构，可对所有子区独立估值，完成图像间和子区间的并行计算。最后对丁腈橡胶进行了拉伸测试，结果表明相比于传统的串行DIC计算方法，运用本文的并行方法速度可提升两个数量级。

关键词

引言

数字图像相关（Digital image correlation，DIC）^［

1‑3］是一种非接触式光学测量方法，通过电荷耦合元件（Charge coupled device，CCD）相机记录试件变形前后信息，可获得其位移、变形等力学性能。DIC技术装置简单，具有很强的可重复分析的特点，在科学研究和工程实践中得到了广泛应用。

随着高速摄像技术的发展，DIC方法越来越多地用于碰撞^［

4］和材料的瞬态变形测量^{［参考文献 5

百度学术}5］。这些实验一般都以大于100 frame/s频率采集图像，在一些高速运动中甚至超过1 000 frame/s。在此过程中会记录下大量的图像，这些图像记录了材料变形直至破坏的全过程。由于图像的数据量大，导致了相关计算和分析非常耗时，所以计算速度的提升有利于数据分析的时效性。邵新星等^{［参考文献 6

百度学术}6］提出了一种多种子点扩散并行计算方法用于土木的准静态测量。文献［7］提出了一种基于CPU与图形处理器（Graphics processing unit，GPU）的异构框架，实现了子区内和子区间的并行计算。但是到目前为止，DIC并行算法的研究大部分都集中于子区间并行，未能达到图像间并行的原因是由于在碰撞、快速变形的情况下，试件会出现区域变形较大的情况，导致DIC方法出现去相关现象，造成较大的测量误差。

根据现状，本文基于参考图更新方法和GPU并行计算能力，实现了大变形场下高速DIC并行计算测量。提出了一种基于Accelerated‑KAZE（AKAZE）的参考图更新方法，可在DIC正式计算之前完成参考图更新和分组工作，为图像间并行计算提供独立数据；并构建了GPU并行计算架构，可独立获取各子区初值，完成多张变形图的并行计算，提升检测效率。

1 基本原理

DIC基本原理是在参考图中定义一个大小为 $(2 M + 1) 像素 \times (2 M + 1)$ 像素的参考子区，由相关搜索算法在变形图中找到相对应的目标子区。计算目标子区中心点 $(x_{0}^{'}, y_{0}^{'})$ 和参考子区中心点 $(x_{0}, y_{0})$ 的差值初步得到子区的位移信息。再通过优化函数得到更精确的位移分量 $u, v$ 。计算感兴趣区域（Region of interest， ROI）内的全部感兴趣点就可得到全场变形参数。

1.1　AKAZE特征点检测算法

为了提高特征点的检测速度，文献［

8‑9］在2013年提出了AKAZE特征点检测算法，该算法实现了快速显示扩散（Fast explicit diffusion，FED）^{［参考文献 10

百度学术}10］来构建非线性尺度空间。在对特征点描述上采用局部差分二值化（Modified‑local difference binary，M‑LDB）^{［参考文献 11

百度学术}11］构建信息描述子，使特征点具备旋转不变性且占用更少的存储空间。不同于传统的基于线性高斯金字塔的尺度不变特征变换（Scale‑invariant feature transform，SIFT）^{［参考文献 12

百度学术}12］和加速健壮特征（Speeded up robust features，SURF）^{［参考文献 13

百度学术}13］，AKAZE可局部自适应小细节，保留目标边界，让尺度空间具备更多的特征信息，提升算法精度^{［参考文献 9

百度学术}9］。

1.2　参考图更新方法

对“大变形”的定义是传统DIC计算时产生去相关而误匹配的变形，量化来看是参考图和变形图变形量超过20%或者旋转角度大于7°的变形。更新参考图可降低大变形引起的去相关效应，这里使用基于AKAZE的参考图更新方法。首先提取参考图和变形图的3个最强匹配点，再计算以这3对匹配点为中心构建的子区相关系数大小，来判断是否需要更新参考图。为了减少参考图和变形图的计算次数，采用区间二分法对采集到的所有序列图像N进行分组计算，具体步骤如下：

（1）以12为单位划分为 $k$ 个区间。

（2）令 $D_{i}^{k} = {I_{n = 1,2, \dots, N}}, i \in [1,12], k \in [1, \frac{N}{12} + 1]$ ， $D_{i}^{k}$ 为第 $k$ 个区间的第 $i$ 张变形图，用二分法选取变形图和参考图进行AKAZE特征点检测，构建待测子区 $S_{R}$ （参考图子区）， $S_{i}^{k}$ （变形图子区）。

（3）定义函数

F (D_{i}^{k}) = C (I_{R}, I_{12 (k - 1) + i})

（1）

式中： $C (I_{R}, I_{12 (k - 1) + i})$ 为参考图 $I_{R}$ 与第 $k$ 个区间的第 $i$ 张变形图 $I_{12 (k - 1) + i}$ 待测子区的归一化互相关（Normalized cross‑correlation，NCC）系数。

若 $F (D_{i}^{k}) \geq C_{T h}$ 且 $F (D_{i + 1}^{k}) < C_{T h}$ （ $C_{T h}$ 为设定阈值），则确定第1分组 $G_{1} {I_{1}, I_{2}, \dots, I_{12 (k - 1) + i}}$ 。令 $I_{12 (k - 1) + i}$ 为第2组参考图，以 $I_{12 (k - 1) + i + 1}$ 为第1张变形图重新划分区间，重复步骤（1~3），直至分组完毕。

1.3　基于GPU的DIC并行计算

GPU具有数量众多的计算逻辑单元，通过计算机统一设备架构（Computer unified device architecture，CUDA）编程平台对其调用可以完成复杂的并行计算任务。如图1为大变形场的DIC计算架构，共分为两部分，第1部分为参考图变换与分组，第2部分为2D‑DIC计算。经过参考图变换后，选定一组图像进行2D‑DIC多图像间并行计算，计算包括预计算部分和迭代部分。预计算部分主要在CPU内完成，包括子区的初值估计，各像素空间位置索引，Hessian矩阵和图像梯度的计算。迭代部分主要在GPU内完成，包括形函数计算，双三次插值，迭代参量 $p$ 的计算等，通过调用线程模块，开启子区间和图像间的并行计算模式。

图1 大变形场2D-DIC计算架构

Fig.1 Large deformation field 2D-DIC computing architecture

图中I_i^R（i∈［1，n］）为每个分组的参考图，I_i_j^D、I_i_j^D_subset（i∈［1，n］，j $\in$ ［1，N］）分别为为每个分组的变形图和子区，x_n_subset为变形图中对应参考图的亚像素位置，G_subset为亚像素位置的像素值，p_ij和p_i_j^converged分别为各子区的迭代参量和收敛值，u和v分别为收敛后获得的x，y方向的位移信息。

1.3.1　初值估计

为了打破子区间初值传播的依赖性，文献［

7］用3层金字塔对子区进行初值估计，该方法虽然能够独立估计出每个子区的初值，但是需要连续的3次DIC计算才能得到两张图像的变形信息，花费了过多的时间成本。基于此，这里提出一种基于AKAZE的全局子区估值方法。如图2所示，第1步对参考图和变形图的ROI区域进行AKAZE特征点检测，得到可两两匹配的特征点对，第2步对ROI区域进行网格划分，网格大小与子区大小相同，第3步对落在每一个网格内的特征点求位移均值

u 、 v

，得到子区初值

u = (\sum x_{n}^{'} - x_{n}) / n

（2）

v = (\sum y_{n}^{'} - y_{n}) / n

（3）

式中： $x_{n}$ 、 $y_{n}$ 和 $x_{n}^{'}$ 、 $y_{n}^{'}$ 分别为参考图与变形图相对应的特征点； $n$ 为单个子区内特征点存在个数。该方法可减少过多的冗余计算，还能得到更精确的初值估计。

图2 各子区初值估计

Fig.2 Estimated initial values for each subsize

1.3.2　并行计算形函数

本文采用一阶形函数^［

14］来跟踪试件6个自由度的空间位置，从Algorithm 1可看出，程序执行了2个数据无关的操作，将一阶形函数作用在参考图的独立子区，得到一个新的坐标值，即变形图中对应子区位置。使用GPU计算时，传入形函数参数

p^{'}

和参考图的像素位置

(x, y)

。启动N×t×g个并行线程对不同变形图多子区进行并行计算。其中，N为子区像素个数，t为ROI包含的子区个数，g为同时计算的变形图张数。

Algorithm 1 Shape function computation

（1） Input： $x, y$ for N×t×g pixels in ROI， $p^{'}$ of size 6×t×g

（2） for N×t×g parallel threads do

（3） $t i d = g e t T h r e a d I d ()$

（4） $x_{n} [t i d] = p_{0} + (1 + p_{2}) \cdot x [t i d] + p_{3} \cdot y [t i d]$

（5） $y_{n} [t i d] = p_{1} + p_{4} \cdot x [t i d] + (1 + p_{5}) \cdot y [t i d]$

（6） cuda.synchronize（）

（7） end for

图3显示了3种不同架构的形函数计算过程，图3（a）为传统串行方法，每个子区逐个计算，计算次数为t×g次；图3（b）为Mullai的子区间并行方法，每张图的所有子区同时计算，计算次数g次；图3（c）为本文图像间与子区间并行方法，g张变形图所有子区同时计算，计算次数1次。

图3 形函数计算过程

Fig.3 Shape function calculation process

1.3.3　双三次插值

这里采用双三次插值法，通过对待测点周围最近的16个采样点加权平均得到插值结果

f (x, y) = \sum_{i = 0}^{3} \sum_{j = 0}^{3} f (x_{i}, y_{j}) W (x - x_{i}) W (y - y_{j})

（4）

式中： $W (x)$ 为权重函数，如Algorithm 2所示，输入形函数的计算结果和双三次插值权重的查找表，输出亚像素插值结果。

Algorithm 2 Bicubic interpolation computation

（1） Input： $x_{n}, y_{n}$ for N×t×g pixels， lookup‑table

（2） Output： G of size N×t×g

（3） for N×t×g parallel threads do

（4） $t i d = g e t T h r e a d I d ()$

（5） $a_{i} [t i d] = B a s i s W e i g h t (l o o k u p_t a b l e)$

（6） $b_{i} [t i d] = G e t (x_{n} [t i d], y_{n} [t i d])$

（7） $G [t i d] = \sum_{i = 1}^{16} b_{i} [t i d] a_{i} [t i d]$

（8） cuda.synchronize（）

（9） end for

1.3.4　并行计算迭代参量 $p$

本文采用反向组合高斯牛顿（Inverse compositional Gauss‑Newton，IC‑GN）算法^［

15］对迭代参量进行优化，其形式为

C_{Z N C C} (Δ p) = \sum_{ξ} {\{\frac{[f (x + w (ξ; Δ p)) - \bar{f}]}{Δ f} - \frac{[g (x + w (ξ; p)) - \bar{g}]}{Δ g}\}}^{2}

（5）

为方便并行计算，对函数优化如下：

（1）计算参考图和变形图各子区像素均值，开启线程数 $t i d = t \times g$ 。

\bar{F} [t i d] = \frac{\sum_{n = 1}^{N} F_{n}}{N}

（6）

\bar{G} [t i d] = \frac{\sum_{n = 1}^{N} G_{n}}{N}

（7）

（2）计算参考图和变形图各像素值与所在子区像素均值的强度差，开启线程数 $t i d = N \times t \times g$ 。

Δ f_{t e m p} [t i d] = F [t i d] - \bar{F} [t]

（8）

Δ g_{t e m p} [t i d] = G [t i d] - \bar{G} [t]

（9）

（3）计算 $Δ f$ 和 $Δ g$ ，开启线程数 $t i d = g$ 。

Δ f [t i d] = \sqrt[]{\sum_{n = 1}^{N \times t} Δ f_{t e m p}^{2}}

（10）

Δ g [t i d] = \sqrt[]{\sum_{n = 1}^{N \times t} Δ g_{t e m p}^{2}}

（11）

（4）计算乘积分量 $Q_{t e m p}$ ，开启线程数 $t i d = N \times t \times g$ ， $t = t i d / N$ ， $t^{'} = t i d / (N \cdot t)$ 。

Q_{t e m p} [t i d] = \nabla f \frac{\partial w}{\partial p} \times [(F [t i d] - \bar{F} [t]) - (Δ f [t^{'}] / Δ g [t^{'}]) \times (G [t i d] - \bar{G} [t])]

（12）

（5）计算乘积累加项 $Q$ ，开启线程数 $t i d = t \times g$ 。

Q [t i d] = \sum_{n}^{N} Q_{t e m p} [n]

（13）

（6）计算迭代增量 $Δ p$ ，更新迭代参量 $p^{'}$ ，开启线程数 $t i d = t \times g$ 。

Δ p [t i d] = H^{- 1} [t i d] \times Q [t i d]

（14）

p^{'} [t i d] \leftarrow p [t i d] \cdot Δ p [t i d] -

1

（15）

以上每个分部计算完成都会执行一个线程同步操作，来确保当前所有线程计算完毕之后再开启下一个任务。

2 分析与讨论

2.1　参考图更新仿真测试

由计算机生成的数字散斑可以有效地控制变形信息，便于评估方法性能。本小节根据文献［

16］的模拟方法生成一系列拉伸变形的仿真散斑。实验测试了两组图像，尺寸大小分别为

512 像素 \times 512

像素和

1 000 像素 \times 500

像素。实验结果如图4所示，随着变形图计算序列的增加，3个测试子区的

C_{N C C}

系数逐渐降低，到第10张变形图时，相关系数均低于所设阈值0.9，说明去相关效应严重，匹配结果不佳。对参考图进行更新，将第9张变形图作为新的参考图与后续的变形图继续比较。更新参考图后，第10张变形图与新参考图的

C_{N C C}

相关系数再次增大到0.99范围，说明匹配可靠，保真度高。基于AKAZE特征点检测的参考图变换，能够迅速精准地找出参考图和变形图的匹配点，再配以区间二分法对整个序列图像分组，可以有效地减少参考图和变形图的计算次数。

图4 3个对应待测子区 $C_{N C C}$ 系数

Fig.4 Three $C_{N C C}$ coefficients corresponding to the subsize to be tested

2.2　GPU并行测试

为验证所提出并行算法的准确性和鲁棒性，本文对非线性材料丁腈橡胶（Acrylonitrile‑butadiene rubber， NBR ）标准试件进行测量。按照GB/T528 2009标准《硫化橡胶或热塑性橡胶拉伸应力应变性能的测定》中规定尺寸，制成宽度为25.00 mm，长度为100.00 mm，拉伸段长度为20.00 mm，厚度为2.00 mm。实验采用RGM‑6000万能试验机对标准试件进行拉伸实验，测试之前对试件表面喷涂白色涂料，形成随机散斑图案。以1 mm/s速度向下拉伸，拉伸过程中丁腈橡胶逐渐变长直至断裂。采集过程中用CCD相机（大恒MER‑500‑7UC，分辨率1 944像素 $\times$ 2 592像素）对试件以400 ms/帧的速度拍摄记录。实验过程中始终保持室温25℃。用于计算的硬件配置包括Inter i5‑10400 6核CPU，NVIDIA GeForce GTX760 GPU。软件配置包括CUDA10 Toolkit等。

通过实验采集到了一系列丁腈橡胶变形实验图，图像序列记录了试件沿y轴单向拉伸下的变形，经过2.2节方法更新参考图后，得到了可连续计算的多组图像，每组图像包含了一张参考图和多张变形图。这里将第1组第1帧作为参考图像，第2~5帧作为变形图像，ROI为 $840 像素 \times 315$ 像素，子区大小为 $21 像素 \times 21$ 像素，进行DIC并行计算。

图5给出了沿y轴的位移分布图，随着拉力的增大，丁腈橡胶逐渐下移，并且越靠下半部分颜色越深，位移越大，符合实际拉动的预设值。为了说明本文并行计算方法的优越性，将其与相同的串行DIC计算方法做了比较，方便起见分别称之为paDIC和seDIC。沿y方向上的位移误差如图6所示，平均误差在 $2 \times 10^{- 3}$ 像素左右。表1统计了4张变形图的计算时间。seDIC和paDIC计算总时长分别为44 721 ms和726 ms，计算速度总体提升了62倍左右。从表1中还可看出双三次插值操作和计算迭代参量 $p$ 的速度提升效果远大于形函数计算，这是由于这两部分计算复杂度较高，利用GPU进行并行计算效果更明显。

图5 随拉力增大试件表面y方向位移场分布图

Fig.5 Distribution of the displacement field in the y direction on the surface of the specimen with increasing tension

图6 y方向位移误差

Fig.6 y-direction displacement error

表1 计算时间对比表

Table 1 Comparison of computation time ( ms )

项目	seDIC	paDIC	Speed‑up
Overall	44 721	726	61.6
Shape function computation	170	9.2	18.5
Shape function computation	11 805	25	472.2
Bicubic interpolation computation	11 805	25	472.2
Bicubic interpolation computation	27 954	118.2	236.5
Computate $p$	27 954	118.2	236.5

3 结束语

本文就大变形场下DIC计算产生的去相关效应进行了讨论，提出了一种基于AKAZE的参考图更新方法，有效地解决了变形图与参考图匹配度低的问题，并打破了图像间的依赖性，为图像间的并行计算提供了独立数据。在GPU并行计算的支持下，评估了高精度亚像素级的DIC计算性能，并在商用计算机上DIC计算速度提升了两个数量级。值得指出的是，通过优化代码，或者在更专业的GPU设备上可同时计算更多的图像，仍可进一步提升计算效率。

参考文献

Pan B, Li K, Tong W. Fast, robust and accurate digital image correlation calculation without redundant computations[J]. Experimental Mechanics, 2013, 53(7): 1277-1289. [百度学术]

Sakai K, Zhang Y L, Yoneyama S, et al. Evaluating distribution and variation of strains near eyes in blinking using digital image correlation[J]. Skin Research and Technology, 2020, 26(5): 749-756. [百度学术]

Xu 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. [百度学术]

王永红,但西佐,胡悦,等.基于高速数字图像相关的人车碰撞伤害实验研究[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 correlation[J].Journal of Optoelectronics Laser, 2017,28(1): 81-86. [百度学术]

申海艇,蒋招绣,王贝壳,等.基于超高速相机的数字图像相关性全场应变分析在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 experiment[J].Explosion and Shock Waves,2017,37(1): 15-20. [百度学术]

邵新星,戴云彤,何小元,等.实时数字图像相关用于土木准静态实验测量[J].光学学报,2015,35(10): 133-141. [百度学术]

Shao Xinxing, Dai Yuntong, He Xiaoyuan, et al.Real-time digital image correlation for civil quasi-static experimental measurement[J].Acta Optica Sinica,2015,35(10): 133-141. [百度学术]

Thiagu M, Subramanian S J, Nasre R. High-speed, two-dimensional digital image correlation algorithm using heterogeneous (CPU-GPU) framework[J]. Strain, 2020, 56(3): e12342. [百度学术]

Alcantarilla P F, Bartoli A, Davison A J. KAZE features[C]//Proceedings of European Conference on Computer Vision. Berlin, Germany: Springer-Verlag, 2012:214‑227. [百度学术]

Alcantarilla P F, Nuevo J, Bartoli A. Fast explicit diffusion for accelerated features in nonlinear scale spaces[C]//Proceedings of the British Machine Vision Conference. [S.l.]: BMVA,2013. [百度学术]

Schreier H W, Sutton M A. Systematic errors in digital image correlation due to undermatched subset shape functions[J]. Experimental Mechanics, 2002, 42: 303-310 . [百度学术]

Yang X, Cheng K T. LDB: An ultra-fast feature for scalable augmented reality[C]//Proceedings of IEEE International Symposium on Mixed and Augmented Reality (ISMAR). Atlanta, USA: IEEE, 2012: 49-57. [百度学术]

Lowe D G. Distinctive image features from scale-invariant keypoints[J]. International Journal of Computer Vision, 2004, 60(2): 91-110. [百度学术]

Bay H, Tuytelaars T, Gool L V. SURF: Speeded up robust features[C]//Proceedings of the 9th European Conference on Computer Vision—Volume Part I.[S.l.]: Springer-Verlag, 2006: 346-359. [百度学术]

Xu X H, Su Y, Cai Y, et al. Effects of various shape functions and subset size in local deformation measurements using DIC[J]. Experimental Mechanics, 2015, 55: 1575‑1590 . [百度学术]

GAO Y, CHENG T, SU Yong, et al. High-efficiency and high-accuracy digital image correlation for three-dimensional measurement[J]. Optics and Lasers in Engineering,2015,65: 73-80. [百度学术]

Su Y, Gao Z R, Fang Z, et al. Theoretical analysis on performance of digital speckle pattern: Uniqueness, accuracy, precision, and spatial resolution[J]. Optics Express, 2019, 27(16): 22439‑22474. [百度学术]

您是本站第访问者

通信地址：南京市御道街29号南京航空航天大学(明故宫校区)

邮编：210016 传真：025-84892742

电话：025-84892742 E-mail：sjcj@nuaa.edu.cn

技术支持：北京勤云科技发展有限公司

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

摘要

关键词

引 言

1 基本原理

1.1 AKAZE特征点检测算法

1.2 参考图更新方法

1.3 基于GPU的DIC并行计算