用于物体内部变形测量的数字图像相关技术

1. 引言

二维DIC使用一个相机，可以测量平面物体的面内位移及应变，三维DIC使用两个相机结合立体视觉与DIC技术，可以测量任意形状物体的三维位移及表面应变。然而，无论是二维DIC还是三维DIC，都是通过相机拍摄被测物体的表面图像来测量物体表面的变形信息，无法获取物体内部的变形信息。由于物体材质的非均匀性和多向异性，仅靠测量物体表面的表面变形信息无法全面反映物体的变形情况，因此如何获取物体内部的变形信息是力学研究领域的一个重要问题。

近年来，随着各种新型三维图像采集设备（如显微CT、共聚焦显微镜以及磁共振成像设备等）的发展，物体内部信息的获取也变的越来越容易，为物体内部变形的观察提供了更好的条件，也使得数字图像相关法有了更广阔的应用前景。B.K.Bay等人首先提出了数字体图像相关法的概念[125-129]，并将其应用于骨试样压缩实验的内部变形场测量。作为二维DIC的一种扩展，VDIC的概念非常简单明了，但是在实际的使用过程中却存在着许多包括计算复杂性、测量精度以及计算效率等诸多难题[130]。虽然与二维DIC相比，VDIC仅仅是增加一个维度，但是计算量和计算的复杂度却大大增加。因此对于VDIC来说，计算精度和速度是两个核心问题。

西安交通大学机械工程学院先进制造技术研究所对VDIC的关键技术进行了研究。首先引入基于最小二乘迭代法的亚像素匹配方法，将这种方法与基于种子点的图像匹配方法相结合，在保证图像匹配精度的同时可以获得较高的计算效率。然后通过对模拟变形体图像的测量和微米X射线三维成像系统测量泡沫铝平移实验来证明该方法的有效性。

2.数字体图像相关技术

与DIC数字图像相关技术不同，VDIC处理的不是普通相机采集的二维图片，而是通过体成像设备（如显微CT、共聚焦显微镜以及磁共振成像设备等）获取的物体内部变形前后的体图像。如图1所示，VDIC相关匹配的过程中涉及到变形前的参考体图像（用 表示）和变形后的变形体图像（用 表示）。

通过精确跟踪参考体图像中的离散点在变形体图像中的位置来获取各离散点的三维位移量。以待匹配点为例，在计算过程中，在参考体图像中取以点为中心的大小的立方体子体块，在变形体图像中，通过一定的搜索算法并按照某一相关函数进行相关计算，寻找与选定的参考子体块相关系数最大的以点为中心的变形子体块。则点即为在变形体图像中的对应点。

图5 1   数字体图像相关方法测量物体内部三维位移原理示意图

在进行相关匹配时，参考子体块与变形子体块之间需要按照一定的映射函数进行体素三维坐标的对应，考虑到变形引起的参考子体块的形状变化，本文采用一阶映射函数，其中包含旋转、平移、正应变、剪切应变等分量，假设参考子体块中的任意一点在变形子体块中的对应点为，则二者之间存在如下关系：

（5-1）式中：——参考子体块中任意一点 到中心点 的距离； ——x，y，z三个方向的位移分量；

                      ——子体块内部的位移梯度。

对于变形前后子体块中的对应点和，考虑到噪声和光照引起的灰度差异，二者的灰度值和存在如下关系：

（5-2）式中：

                    E——噪声；

                    ——用来补偿光照引起灰度差异的参数。

设参考子体块内共有n个像素，且像素的灰度值受理想的随机分布噪声影响，那么通过求取相关系数的极值，可以在变形体图像中找到与参考子体块具有最大相似度的目标子区，即：（5-3）式中：为相关参数向量。对方程（5-3）可以通过最小二乘迭代算法进行求解。

3 基于最小二乘迭代的VDIC亚像素匹配方法

针对二维DIC的最小二乘迭代算法在2.4节已有论述，本节在此基础上介绍针对VIDC的最小二乘迭代算法。对公式（5-2）进行泰勒分解并去掉高于一阶的项，可得误差方程：

（5-4）设未知数的改正数向量为x，偏导数矩阵为A，观察值的向量（参考图片与待匹配图像像素值的差）为l，假设子图像的像素数目为n，则对应的误差方程组为：

（5-5）式中：
利用最小二乘的方法求解方程可得：
其中P为权值，一般情况下可以令P为单位矩阵，也即所有的观测量的权值相等。求得x之后，对未知数进行更新，并重复上述过程，直到改正数x小于一定的阀值。

针对VDIC的最小二乘迭代算法与二维DIC原理相同，但是未知数参数数量增加了，计算复杂度也相应的提高了。最小二乘迭代算法是针对一个子体块的优化算法，如前说述，这种算法要求有比较准确的未知数初值，初值不准确会导致迭代次数增加甚至迭代不收敛。为了求取其相关参数初值，传统的方法是，先利用整数像素的粗匹配得到位移的初值，其余参数则取固定的初值，粗匹配计算量很大、效率很低，对于VDIC则更加明显，因此与二维DIC一样也需要配合基于种子点的图像匹配算法，不同的是，对于DIC来说，基于种子点的扩散匹配是在种子点的4临域进行的，而对于VDIC来说，基于种子点的扩散匹配是在种子点的6临域进行的。

4 灰度插值算法

与二维DIC一样，VDIC同样需要对非整数像素进行插值计算，所不同的是插值维度变成三维，计算难度也相应增加不少。如果出于减少计算量的考虑，可以选用三维线性插值。三维线性插值相当于2次双线程插值和1次线性插值的组合。如图2所示，设 点为待插值点， 点的非整数像素坐标为 ，三维线性插值使用周围临近的八个点进行插值计算。首先使用 和 分别进行双线性插值（见1）节）计算获得 和 的灰度及梯度值，然后对 和 进行一次线性插值，即可得到 点的灰度及梯度值：

（5-6）三维线性插值的优势是计算速度高，但是插值精度较低，为了提高插值精度，则可以选用三维三次插值算法。如图3所示，三维三次插值使用周围的64个点进行插值计算，三维三次插值算法可以分解为4次双三次插值和1次一维三次插值。设待插值点的非整数像素坐标为 ，首先在z方向的4个整数像素位置上，分别进行双三次插值，求取四个点的灰度及梯度值，然后在z方向上进行一次一维三次插值，得到待插值点的灰度及梯度值。

由于VDIC的灰度插值计算量非常巨大，无法使用查表法进行灰度插值。为了提高计算效率，本文选用的是三维线性插值算法。

图 2   三维线性插值示意图

图 3   三维三次插值示意图

5 实验验证

5.1 三维散斑模拟体图像生成

本文在2.5节介绍了二维数字散斑模拟图像的生成方法，模拟数字散斑图像是通过算法计算出来的模拟散斑。模拟数字散斑图像的位移、应变量，散斑的大小、密度、灰度以及噪声都可以人为控制，因此可以广泛应用于数字图像相关法的算法研究中。

三维散斑模拟体图像是从二维数字散斑模拟图像扩展而来的，三维散斑模拟体图像是由多个位置随机的空间散斑颗粒的作用叠加而成，在三维空间中每个散斑颗粒可以看作是一个点光源，每个空间体素的灰度是由所有这些点光源的对它的影响叠加而成，具体计算公式见式（5-7），按照这个公式可以计算生成变形前后的两张三维散斑模拟体图像，其中 为变形前的参考图像， 为变形后的图像。

（5-7）式中：

s——散斑颗粒的总个数；

R——散斑颗粒的尺寸；

——每个散斑颗粒的随机光强；

——变形前第k个空间散斑颗粒的随机分布位置；

——变形后第k个空间散斑颗粒的随机分布位置。

模拟过程中可以使用公式（5-8）模拟实际变形过程中可能发生的旋转、平移、拉压变形及剪切变形。

（5-8）
为了验证本文算法的有效性及精度，下面利用模拟算法实施两个实验：刚体平移实验和单向拉伸实验。对于刚体平移实验，设定的模拟参数是：体图像尺寸大小为150×150×150体素，s=40000，R=2， ，其余参数均为0。对于单向拉伸实验，设定的模拟参数是：体图像尺寸大小为150×150×150体素，s=40000，R=2， ，其他参数均为0。图4显示的是模拟单向拉伸实验的150×150×150大小的两幅三维散斑模拟体图像，与参考体图像相比，变形体图像在z方向发生了均匀拉伸变形。

（a）参考体图像                   （b）变形体图像
图4   三维散斑模拟体图像

5.2 模拟刚体平移实验

首先利用前一节模拟得到的刚体平移模拟体图像测试本身的匹配算法。对于平移实验的第一部分模拟图像，在计算过程中，分别设置子体块尺寸为15×15×15体素、25×25×25体素、35×35×35体素，步长设置为15体素，选定计算区域为中间的135×135×135大小的散斑体图像，因此需要匹配的点个数为9×9×9=729个。

在匹配过程中，为了提高匹配效率，仍然使用基于种子点的图像匹配方法（见3.2.1节），首先选取一个子体块作为种子点并对种子点进行匹配计算，对于种子的的匹配使用整像素搜索加最小二乘迭代优化的方法。种子点匹配成功后，利用种子点为其周围临近点提供匹配初值并完成这些临近点的匹配，如此不断向外扩散，直至所有的点匹配完成。通过这种方法，可以省去大量整像素搜索的繁琐工作量，提高匹配速度。

针对不同的子体块设置，共得到三组不同的计算结果，如表 5 1所示。从表中的结果可以看出，随着子体块尺寸的增加，位移平均偏差没有显著变化，但是位移标准差显著减低，同时计算时间也有较大增加。

表 5 1   针对不同子体块设置的计算结果

子体素尺寸/体素 位移X平均偏差/体素 位移Y平均偏差/体素 位移Z平均偏差/体素 位移X平标准差/体素 位移Y标准差/体素 位移Z标准差/体素 平均计算时间/秒

15×15×15 0.00416 -0.00536 0.00674 0.00897 0.00905 0.01381 0.01899

25×25×25 0.00349 -0.00493 0.00772   0.00303 0.00335 0.00772 0.0499

35×35×35 0.00338 -0.00484 0.00698   0.00141 0.00169 0.00546   0.1058

5.3 模拟单向拉伸实验

利用前一节模拟得到的单向拉伸三维散斑模拟体图像进一步测试本文的匹配算法。在计算过程中，设置子体块尺寸为15×15×15体素，步长设置为15体素，选定计算区域为中间的135×135×135大小的散斑体图像，因此需要匹配的点个数为9×9×9=729个。图5（a）显示的是变形前划分好的子体块线框，图5（b）显示的变形后的子体块线框，图6（a）为变形后的三维点，图6（b）则是以三维色谱图形式显示的位移测量结果。从图中可以看出，位移场沿z轴方向线性增大，范围为-0.327~-3.893，测量结果与预先施加的单向拉伸变形完全符合，最大偏差小于0.02个体素。

          （a）变形前                （b）变形后
图5   变形前后的子体块线框

（a）三维点           （b）三维位移色谱图
图6   变形后的三维点及三维位移色谱图

5.4 泡沫铝平移测量实验

通过真实实验验证DVC方法的有效性，测量对象是泡沫铝。泡沫铝是在纯铝或铝合金中加入添加剂后，经过发泡工艺而成，同时兼有金属和气泡特征。它密度小、高吸收冲击能力强、耐高温、防火性能强、抗腐蚀、隔音降噪、导热率低、电磁屏蔽性高、耐候性强、有过滤能力、易加工、易安装、成形精度高、可进行表面涂装。泡沫铝具有优异的物理性能、化学性能和力学性能以及可回收性。泡沫铝的这些优异性能使其在当今的材料领域具有广阔的应用前景，是很有开发前途的工程材料，特别是在交通运输工业，航天事业和建筑结构工业等方面。

实验通过微米X射线三维成像系统测量泡沫铝的内部结构，如图7所示。多功能微焦点X 射线检测系统，结合了Feinfocus的专利技术与先进的高速平板检测技术，可以通过X 射线检测样品三维立体形貌及内部结构。可进行电子产品、微系统、产品装配、原材料及成品的检测，满足新产品研发、原型机设计、失效分板及生产工艺流程控制需求。其性能参数如下：

（1）最大电压：160 KV；最大靶功率：10W；三维分辨率：1.5 μm。

（2）最大样品重量：5 kg；最大样品尺寸：440 mm×550 mm。

（3）操作软件：FGUI 3.0（二维），VG Studio 2.1（三维）。

（4）标准操作轴：X，Y，Z，探测头倾斜角度：145 度 (±72.5°)。

图8是扫描获得的泡沫铝三维立体形貌。通过切片操作可以得到泡沫铝的内部结构图片，见图9。 X射线三维成像系统在X，Y方向的分辨率为512voxel×512voxel，设置切片距离为0.078mm，一共获得内部结构图片512张。第一次测量完成后，沿Z方向将泡沫铝平移一段距离，并进行第二次测量。利用DVC方法对两组图像进行匹配分析，得到的泡沫铝在Z方向的位移色谱图见图10，Z方向位移的平均值为15.13voxel（约1.19mm），标准偏差为0.053voxel（约0.004mm）。虽然实际的平移距离是未知数，但是测量结果与实际情况相符，且标准偏差较小，证明DVC方法可以用于真实固体材料的内部变形测量。

图7   泡沫铝测量实验现场

图9   泡沫铝内部结构

图10   三维位移色谱图

6 试验总结

针对用于物体内部变形测量的VDIC进行了深入研究，完成的主要工作有：

（1）提出了一种基于最小二乘迭代的VDIC亚像素匹配算法。将这种方法与基于种子点的图像匹配方法相结合，在保证图像匹配精度的同时可以获得较高的计算效率。

（2）研究了VDIC中的三维灰度插值算法。VDIC的是插值维度变成三维，计算难度也相应增加不少。出于减少计算量的考虑，本文采用的是三维线性插值。

（3）实验验证了基于最小二乘迭代的VDIC亚像素匹配算法的有效性，分析了不同子体块尺寸对位移测量精度的影响。通过微米X射线三维成像系统测量泡沫铝的刚体平移，进一步验证了本文算法的有效性，也证明DVC方法可以用于真实固体材料的内部变形测量。