邹 翔 张轩豪 王延珺 潘 兵,2)

* (北京航空航天大学固体力学所,北京 100191)

† (中国工程物理研究院化工材料研究所,四川绵阳 621999)

引言

数字体图像相关方法(digital volume correlation,DVC)是二维数字图像相关方法(two-dimensional digital image correlation,2D DIC)在三维体图像上的拓展.通过比较体成像设备获取的被测试样变形前后数字体图像,该方法可测量物体内部三维全场变形[1-2].由于可以提供被测试样内部(而不仅仅是表面)更为丰富的变形信息,DVC 在多孔材料(如骨[3]、金属泡沫[4]、木材[5])、半透明生物材料(如细胞[6]等)、建筑材料(如混凝土[7]等)、复杂结构(如3D 打印结构[8]、生物植入物[9]等)的内部变形场测量、力学行为表征和有限元模拟互证等方面都有着重要应用[10-12].

DVC 的基本原理是从参考体图像中围绕感兴趣的计算点选择一个立方体的图像子体块(subvolume),通过跟踪该图像子体块在变形后图像中的位置以获得图像子体块中心点的三维位移矢量.在实际分析计算时,被测试样体图像需呈现随机的灰度分布(通常称作内部散斑),该散斑场作为变形信息的载体随试样一起变形.因此被测试样内部散斑质量和DVC 位移测量精度密切相关.为对散斑图质量进行评价,在DIC 方法研究中已提出多个参数[13],这些评价参数可分为局部参数和全局参数[14-16].局部参数主要用于单个图像子区内散斑质量的评价,如Pan 等[17]通过理论分析和模拟实验证明了DIC 图像子区的灰度梯度平方和(sum of square subset intensity gradient,SSSIG)参数与其位移测量结果的随机误差成反比,因此图像子区SSSIG 参数越大,其内部散斑质量越好;Bomarito 等[18]将SSSIG 参数和二次自相关峰高参数结合提出了一种组合式的散斑评价参数,该参数能够解决原有SSSIG 参数不适用于规则散斑图像(如棋盘格图案)的情况.此外,如假设散斑质量在整个图像内的分布是均匀的,则可在散斑质量局部评价参数的基础上提出可评价整幅散斑图质量的全局参数,如Pan 等[19]在SSSIG 局部参数的基础上提出了用于整幅散斑图评价的平均灰度梯度(mean intensity gradient,MIG)全局参数等.

类似地,在利用DVC 测量内部变形时,由于被测对象、体成像设备类型和参数以及实验条件等的差异[20],获取的被测物体内部散斑图像千差万别.因此,如何定量评价图像内部散斑的优劣也是DVC 中一个重要的基本问题.将DIC 中已有的散斑评价参数直接推广到DVC 中是最有效的方法之一.在DIC 中的诸多散斑评价参数中,SSSIG 参数是从DIC 位移测量精度的理论模型中推导出来,具有严格的理论基础,具有普适性.并且该参数是一种典型的局部评价参数,能够克服全局参数在遇到材料内部出现空洞等非连续和非均匀分布时失效的问题[21],因此更适合于材料内部散斑质量的评价.然而,在将SSSIG 参数在DVC 内部散斑评价中推广应用需要解决以下几个问题:(1) SSSIG 参数在DVC 中最新的三维反向组合高斯−牛顿亚体素配准算法[22](3D inverse compositional Gauss-Newton,3D IC-GN)中应用的可行性问题.SSSIG 最初被提出来时是基于DIC 经典亚像素配准算法牛顿−拉普森算法[17],因具有更优的抗噪性能和更高的计算效率,IC-GN 算法已成为DIC/DVC 方法中的标准算法,因此系统地研究SSSIG 参数与该算法测量精度的关系是有必要的;(2)由于DVC 内部散斑类型的多样性,这就要求评价指标应具有广泛的适用性,因此需要研究SSSIG参数在不同类型内部散斑中的适用性;(3)与DIC 通过人工制作散斑不同,DVC 主要通过被测物体自身内部纹理的灰度对比度来形成散斑场.由于被测物体内部组成和结构难以改变,无法通过人工优化DVC 测量对象的内部散斑质量,因此从计算方法层面研究减弱低质量散斑对DVC 位移测量精度的影响无疑是有价值的.

基于以上考虑,本文研究了DVC 中体图像的内部散斑质量评价参数及其适用性,并通过优化选取子体块内的体素点来实现增大子体块尺寸以减小测量误差的同时使计算量没有显著增加.论文首先从DVC 方法的位移测量精度理论分析出发,系统介绍了SSSIG 评价参数的理论基础;然后,通过数值模拟平移实验,检验了该参数在不同类型内部散斑中的有效性;最后,为了在算法层面上改善DVC 子体块的SSSIG 参数,提出了基于灰度梯度的体素点选择DVC 方法,并通过模拟实验和实际实验验证了该方法的有效性.

1 DVC 位移测量误差的理论分析及模拟实验验证

本节首先对最新的3D IC-GN 亚体素配准算法[22]进行简单介绍,然后给出基于3D IC-GN 算法的DVC 位移测量误差的理论公式推导[23-25].该公式显示位移测量误差与SSSIG 参数、体图像噪声、亚体素灰度插值算法直接相关.最后通过模拟数值平移实验验证了SSSIG 参数用于不同类型内部散斑质量评价的有效性.

1.1 3D IC-GN 算法

DVC 方法在具体实现过程中需执行两个步骤:初值估计和亚体素位移测量.对任一计算点,首先需要通过简单的整体素搜索(种子点)或可靠性导向的位移初值传递策略(非种子点)来确定其位移迭代初值.为进一步提高位移测量精度,需进行亚体素配准计算[26-27].在众多DVC 亚体素配准算法中,3D ICGN 算法因具有高精度、高效率的优点而获得广泛使用,该算法利用迭代运算使非线性相关函数由初始估计值逐渐收敛至局部最佳极值(即使参考和变形子体块获得最佳匹配的变形参数矢量)以获得具有亚体素精度的位移计算结果.

3D IC-GN 算法计算流程图如图1 所示,该算法通过优化如式(1)所示的零均值归一化最小平方距离相关函数(zero-mean normalized sum-of-square difference,ZNSSD)获得亚体素位移

图1 3D IC-GN 算法计算流程图[22]Fig.1 Flow chart of the 3D IC-GN algorithm[22]

式中f(x)和g(x) 分别表示参考和变形子体块中点x=(x,y,z)T处的灰度值; ξ=(Δx,Δy,Δz)T表示参考子体块中各体素点的局部坐标;fm和gm分别表示参考和变形子体块的灰度平均值;Δf和Δg分别是参考和变形子体块中所有体素点的灰度标准差;p=(u,ux,uy,uz,v,vx,vy,vz,w,wx,wy,wz)是描述目标子体块的位置和形状变化的变形矢量;Δp=(Δu,Δux,Δuy,Δuz,Δv,Δvx,Δvy,Δvz,Δw,Δwx,Δwy,Δwz)表示参考子体块的增量变形矢量;W(ξ;p) 为目标子体块的翘曲函数;W(ξ;Δp)为参考子体块的增量翘曲函数.

采用高斯牛顿算法对ZNSSD 相关函数进行迭代优化,则变形矢量增量可写为

其中 ∇f是参考体图像中各点的灰度梯度,H为Hessian 矩阵.在确定变形矢量增量Δp之后,就可以确定子体块的增量翘曲函数W(ξ;Δp).然后将该增量翘曲函数与现有的翘曲函数W(ξ;p) 反向组合,以更新目标子体块的翘曲函数.对于每个参考子体块,重复进行式(1)的迭代运算,直至满足预设的收敛条件,从而获取最终优化的变形矢量.此外,为表征目标子体块中更为复杂的局部变形模式,也可采用包含二阶变形分量的翘曲函数[28].

1.2 位移测量误差的理论分析

实际实验中获得的体图像不可避免含有随机噪声[3,28-29],且各体素点的图像噪声相互独立,设参考体图像和变形体图像实际的灰度值分别为(x) 和(x),真实的图像灰度值为f(x) 和g(x),则

式中 δf(x) 和 δg(x) 分别表示在参考体图像和变形体图像中的图像噪声(假设标准差分别为 σf(x) 和σg(x))引起的灰度偏差.因此,在实际实验中,式(2)的实际形式为

根据高斯白噪声的统计特性,变形矢量增量Δp的期望和方差为

式中h(x1)=f(x1)−g(x2) 表示灰度插值偏差.

为了更加直观地展示式(6),将其应用到沿x方向的平移工况中,由上节分析可知,如果u为子体块真实变形量,当采用3D IC-GN 算法迭代1 次时,则位移结果为u′≈u+Δu+uxΔu+uyΔv+uzΔw≈u+Δu,此时移u测量误差ue=u′−u≈Δu,因此由式(6)可知

由式(7)可知,当不考虑形函数欠匹配问题时,在3D IC-GN 算法中位移u的系统误差E(Δu) 和随机误差Var(Δu) 均与子体块灰度梯度和(即SSSIG 参数)大小有关.但系统误差主要受灰度插值误差h(x1) 影响,当采用精度较高的插值算法,并对图像进行预滤波,3D IC-GN 算法的系统误差可控制在较低水平,而随机误差Var(ue) 与子体块SSSIG 参数呈反相关关系.也就是说,在DVC 计算中,某方向的SSSIG 参数越大,其位移测量结果越精确.

1.3 内部散斑质量评价参数的模拟实验验证

为验证上节推导的理论公式的正确性,本文对3 种不同类型内部散斑进行了数值模拟平移实验.与实际实验相比,数值模拟实验不但可以精确控制散斑图的位移和应变信息,还能够排除由体图像采集系统(如设备自热)、非理想加载条件等非理想实验因素对位移测量的影响[15],最大程度地保证了计算结果的可靠性.验证实验选取的参考体图像是由不同X-ray CT 获取的3 种具有代表性的材料(铜颗粒填充环氧树脂基复合材料、柚子皮多孔材料和多聚物粘结糖)的体图像,图2 给出了3 种材料体图像的切片图、以及选取的感兴趣区域(volume of interest,VOI)及其灰度分布直方图.

图2 3 种不同类型内部散斑的切片图、VOI 及其灰度分布直方图Fig.2 Slice figures,VOI and histograms of three speckle patterns

在平移实验中,利用傅里叶变换对3 种材料的体图像做精确的位置平移.具体为:将原体图像作为参考体图像,参考体图像沿x方向平移0.2 体素而得到平移变形体图像.为模拟真实情况下的图像噪声,在参考和变形体图像中分别加入标准差为3 灰度值的高斯白噪声.DVC 分析中均选取大小为200×200×200 体素的体图像中心区域VOI 进行研究,计算采用3D IC-GN 亚体素配准算法和一阶形函数,子体块尺寸为41×41×41 体素,计算步长为10 体素.

图3(a)给出了3 种材料位移u测量结果的标准差(standard deviation (SD) error)与子体块图像x方向SSSIG (简记为SSSIGx)参数平均值的关系的理论结果和测量结果,由表可知,在不同材料中,位移测量结果的标准差与理论预测值吻合良好;不同材料间,散斑图像SSSIGx参数越大的材料,位移测量结果标准差也较小,这和上面的理论公式相一致.需要说明的是,由式(7) 可知,均值误差不仅与图像SSSIG 有关还与子体块内各体素点灰度插值偏差有关,而灰度插值偏差在分析中难以统计,因此本文不再对均值偏差进行重点分析.但是为了降低位移测量结果的均值误差,在DVC 计算前需要对体图像进行低通预滤波[31].

图3 3 种不同类型内部散斑测量结果Fig.3 Measurement results of three speckle patterns

图3(b)还给出了3 种材料中子体块测量值u与其子体块散斑图像参数SSSIGx之间的关系,从图3中可以看出:尽管各材料之间测量结果差别较大,然而对同一散斑材料中,测量结果的变化幅度均整体上与散斑图像SSSIGx参数呈现负相关关系,即随子体块散斑图像参数SSSIGx的增大,其测量结果变化幅度越小;此外,由图3 可知,3 种材料中,柚子皮多孔材料的位移标准差最大,其次是铜颗粒填充环氧树脂基复合材料,而多聚物粘结糖最小.

以上的研究结果表明,散斑图像SSSIGx参数能够反映出材料内部散斑质量,对同一种材料而言,其散斑图像SSSIGx值越大,则DVC 位移测量精度越高.

2 DVC 方法中计算体素点的优化选择

2.1 SSSIG 参数的优化

为提高位移测量精度,在DIC 中通常采用增大子区尺寸来优化子区SSSIG 参数[17,21],该方法也可以推广到DVC 中.但由于DVC 处理的是三维数据,简单地通过增大子体块尺寸来优化子体块SSSIG 参数将大幅度增加计算量,进而影响计算效率.通常,在DVC 计算时子体块内所有的体素点均参与相关性计算,这将耗费大量的计算时间.另外一方面,通过式(7)可以看出子体块的SSSIG 参数实质上其内部所有体素点灰度梯度平方的总和,子体块内灰度梯度小的体素点对提高SSSIG 值影响较小,即这些点对子体块计算精度影响较小,因此可在增大子体块尺寸的同时将子体块内灰度梯度平方偏小的体素点剔除出DVC 计算,以实现在增大子体块SSSIG的同时尽可能的降低计算量.在实际实验需要对3 个方向均进行优化,因此选择了3 个方向的SSSIG参数的均值(简记为SSSIGmean)作为目标函数,见式(8).

为了更加直观地展示材料中灰度梯度的分布情况,图4 给出了多孔材料中一子体块内部灰度梯度和的分布情况,由图可知,在子体块内部有较大一部分体素点的灰度梯度平方和处于较低的水平,因此将这些体素点排除出DVC 计算将能够有效减轻计算量.

图4 多孔材料灰度梯度分布图Fig.4 Intensity gradient distribution of porous materials

2.2 体素点优化选择方法

基于上述的SSSIG 参数优化思想,提出了一种基于灰度梯度的计算体素点优化选择方法,图5 给出了该方法在计算一子体块位移时采取的计算流程图.该方法的具体计算步骤为:

图5 体素点优化选择方法计算流程图Fig.5 Calculation flow chart of voxel selection DVC method

(1) 在每一个子体块开始DVC 计算前首先计算其内部所有体素点3 个方向的灰度梯度平方的总和∇f,只有总和大于所设定灰度梯度阈值 ∇f0的体素点(有效点)才会参与到该子体块SSSIGmean参数的计算和之后的DVC 计算;

(2) 计算该子体块SSSIGmean参数,将该子体块SSSIGmean参数与所设定的SSSIGmean参数阈值T0对比,如果低于该阈值,则将该子体块半边长M增大1 体素;

(3) 重复以上步骤,直到子体块SSSIGmean参数大于所设定的阈值T0,然后基于选定的体素点,完成该子体块的DVC 计算.

在该方法中,体素点灰度梯度阈值 ∇f0和SSSIGmean参数阈值T0的设定是重点,在子体块中,3 个方向的灰度梯度平方的总和小于灰度梯度阈值∇f0的体素点将被删除,而只有高于灰度梯度阈值∇f0的体素点才会参与到DVC 计算.灰度梯度阈值∇f0过小,将起不到提高测量精度的作用,灰度梯度阈值 ∇f0过大,将会使子体块有效体素点大大减小,进而影响测量结果.同时,SSSIGmean参数阈值T0设置过大,则使计算量增加,太小又起不到优化作用,因此实际应用中要根据不同的情况选择合适的灰度梯度阈值 ∇f0和SSSIGmean参数阈值T0.

2.3 模拟实验验证与讨论

为验证体素点优化选择方法的有效性,将该方法应用于1.3 节中模拟平移实验,通过3 种材料原有的∇f的分布情况确定了其灰度梯度阈值 ∇f0和SSSIGmean参数阈值T0,其中子体块灰度梯度阈值 ∇f0的取法为:将子体块所有体素点(记为N)按灰度梯度总和∇f大小递增排列,选用编号N1=int(N/2) 体素点对应的灰度梯度总和为阈值 ∇f0,SSSIGmean参数阈值T0设为该子体块原有SSSIGmean的1.5 倍.

DVC 计算中,子体块尺寸设为41×41×41 体素,计算步长10 体素.图6 给出了常规DVC 和体素点选择DVC 位移测量结果的对比情况.

图6 两种方法测量结果对比Fig.6 Comparison of the results of two methods

由图可知,常规DVC 中每个子体块中体素点的数目是一定的,具体为41×41×41=68 921 个体素点,而在体素点选择DVC 中,尽管子体块尺寸增加了,但3 种材料子体块中有效体素点的平均数量依次为58 619,61 016 和64 521,比常规DVC 略低,这说明优化后的DVC 方法计算量不仅没有增加,还存在一定程度的降低.而在测量误差方面,优化后3 种材料散斑图像的测量标准差分别下降了15.63%,25.00%和18.75%.通过以上的分析结果可知,体素点优化选择方法能够在保证计算精度的前提下有效降低由于子体块尺寸增大所引起的计算量的急剧增加.

2.4 泡沫铜材料重扫描实验

为进一步检验体素点选择DVC 方法的有效性,将该方法应用到泡沫铜材料重扫描实验.

重扫描实验是DVC 中一种常用的测量精度评价方法,在重扫描实验中,被测试件静止,且不对其施加任何载荷,利用体成像设备对试件进行两次连续扫描,通过DVC 分析两次扫描获得的体图像可以获得试件在静止状态下的变形量测量值,通过此实验可以检验DVC 方法的测试精准度.实验所用试件是孔隙率为97.5%泡沫铜材料,试件尺寸为11 mm×11 mm×8 mm.体成像设备选用的是 Skyscan1172 台式 X 射线 CT 系统(Bruker Corporation,美国),空间采样率选取为 2000×1336 像素,体素尺寸设定为9 μm,旋转步长为 0.4°.每次CT 扫描持续 45 min,共获取试件2 组灰度体图像.

图7 给出了CT 实验装置图和试件切片图.在DVC 分析中,采用 3D IC-GN 亚体素配准算法和一阶形函数.感兴趣计算区域选在图像中心,大小为300×300×500 体素,子体块尺寸预设为41×41×41 体素,计算步长为20 体素,共有15×15×25=5625 个规则分布的计算点.体素点选择采用了和3.3 节一致的灰度梯度阈值 ∇f0和SSSIGmean参数阈值T0.

图8 对比给了采用常规和体素点选择DVC 计算获得的位移场,由图可以很直观地看出,采用体素点选择DVC 的测量误差要小于常规DVC.表1 还给出了两种方法计算量和位移测量结果的详细对比情况.由表1 可知,优化后位移u,v,w的标准差分别下降了22.80%,27.73%和28.55%;与此同时,优化前后各位移分量的均值误差变化较小.由此可见,在泡沫铜重扫描实验中,本文提出的体素点选择DVC可有效降低DVC 的位移测量误差,但没有明显增加其计算量.

图8 重扫描实验结果Fig.8 Rescan results

图8 重扫描实验结果 (续)Fig.8 Rescan results (continued)

表1 重扫描实验位移测量结果对比Table 1 Comparison of displacement measurement results in rescan experiment

3 结论

提出了适用于DVC 的内部散斑质量评价参数−SSSIG,并通过模拟平移实验验证了SSSIG 参数在3 种(铜颗粒填充环氧树脂基复合材料、柚子皮多孔材料和多聚物粘结糖)不同类型内部散斑图像中的有效性.结果表明,在DVC 中,子体块散斑图像SSSIG 参数越大,测量误差越小;更进一步,为了在增大子体块SSSIG 参数的同时尽可能的降低计算量,提出了一种体素点优化选择方法,该方法通过将子体块中灰度梯度较小的体素点剔除出DVC 计算来实现在增大子体块尺寸的同时而不增加计算量,从而实现高精度、高效率的DVC 分析.3 种不同类型材料散斑图像的模拟平移实验和泡沫铜的重扫描实验均验证了体素点优化选择方法在有效降低DVC方法的测量误差同时提高计算效率.