增量数字图像相关方法的误差分析及应用考量1)

池宇希 杨昊天 潘 兵

(北京航空航天大学航空科学与工程学院，强度与结构完整性全国重点实验室，北京 100191)

数字图像相关方法（digital image correlation,DIC）通过跟踪变形前后被测物体表面数字图像中的信息载体（如人工散斑、自然纹理等）的位置变化来获得其表面的全场变形[1-2]。由于具有硬件设备和测量过程简单、测量精度高、适用测量范围广等优点，这种基于数字图像的非接触式、全场光测力学方法已成为实验力学领域中最受欢迎的主流变形测试技术。然而在实际工程应用中，当被测物体发生大变形[3]，或应用于高温环境时，剧烈的光照及散斑变化[4-7]导致被测物表面信息变化较大，变形后图像会出现所谓的“退相关效应（decorrelation effect）”。此时，采用固定参考图像的传统DIC（即序列DIC，下同）计算策略往往会失效。为了解决DIC 在上述场景中的应用问题，研究人员发展了增量计算策略[3,5,8-9]，即增量DIC 方法。通过更新参考图像，增量DIC方法提高了当前参考图像与变形图像之间的相关性，使DIC 计算能够顺利进行。在更新参考图像后，通过将各个阶段的位移场累加以得到相对于初始参考图像的位移场，避免了固定初始参考图像时因“退相关效应”导致的计算失效的问题。

除了在上述易发生图像退相关的应用场景外，增量计算策略在其他DIC 领域还有广泛应用。例如，在3D–DIC 中，左右相机分别采集了一组序列图像，为了将左右相机的计算点关联起来，需要在当前帧中进行左右匹配[9]，并通过左相机进行时序匹配；或通过初始参考帧进行左右匹配，并以左右视图各自的初始帧为参考图像进行时序匹配[10]，这些过程均会发生参考图的更新。再如，在多相机DIC 中，相邻相机之间两两匹配时，其参考图像也各不相同，为了在多个视角中确定同一个计算点，也需要对参考图像进行更新。除此之外，在使用数字体图像相关（digital volume correlation, DVC）测量内部大变形[11-13]时，增量算法也被用来匹配多个加载阶段的体图像。在这些需要更新参考图像的场景中，每次更新参考图像时，新的图像子区中心的坐标往往位于亚像素位置，因此需要对参考图像子区进行插值，会带来一定的插值误差。由于标准的反向组合高斯–牛顿算法中[14]，图像梯度和Hessian 矩阵均在参考图像中计算，这些微分操作会进一步放大插值误差，在更新参考图像后，计算结果会呈现出混叠现象[15-16]，故从测量精准度考虑，插值操作应当尽量避免。

为了减小增量算法中由于对参考图像插值带来的插值误差，Zhou 等[17]提出了自适应参考图像子区平移计算策略。在这一策略中，在更新参考图像时，其使用距离亚像素位置计算点最近的整像素点作为新的图像子区中心，避免了对更新参考图像子区的灰度值的插值计算，这一增量策略也被用于DVC 中[13]。由于平移量较小，且由平移带来的位移差可以通过形函数进行补偿。这一策略不仅消除了更新参考图像时由于参考图像插值导致的系统误差（通过自适应平移将其转换成了随机误差），同时也减少了更新过程中由于插值带来的额外计算量。这一策略虽然减轻了由于参考图像插值导致的系统误差，但由变形图插值带来的插值误差仍不可避免。另外，现有参考图更新的策略有：连续更新、固定步长更新和条件更新（即发生图像退相关再进行更新）等策略。其中，连续更新由于频繁更新参考图会带来额外的计算量，也会带来额外的误差累积，而固定步长更新的步长又往往难以选取。所以，在应用增量DIC 计算时使用者会面临两个令人困惑的问题：（1）更新参考图像会不会引入额外的测量误差，以及这一误差的分布范围及分布规律是怎样的？（2）为尽量减少这一误差和计算量，应该采用何种策略更新参考图像？显然，这两个问题是增量DIC/DVC 算法实际使用时必需研究和讨论的问题。

针对以上问题，本文通过一组亚像素平移模拟实验，分析并证明了图像子区平移策略用于亚像素平移测量时由于亚像素插值导致的系统误差的累积效应。并通过控制图像子区平移的步长，总结了增量算法在不同步长条件下对亚像素插值误差的累积规律。另外，通过实际实验，研究了DIC 在实际实验中亚像素插值误差的分布规律，进一步证明了自适应图像子区平移增量算法在实际实验中将系统误差转化为随机误差的性质。最后，结合亚像素子区平移和经典的高斯低通预滤波方法，极大地降低了增量算法中的累积误差，并给出了应用增量DIC 方法时的一些建议。

1 增量数字图像相关方法

1.1 增量数字图像相关计算

序列DIC 算法[14]可实现高效、精准而鲁棒的位移场测量，但在一些容易出现退相关的场景，序列DIC 算法容易失效，需引入增量计算策略。在条件更新策略中，增量算法在当前图像出现图像退相关时会将当前图像的前一帧图像更新为当前参考图像，由于DIC 匹配后，计算点往往处在亚像素位置，直接的做法是对当前亚像素位置的图像子区像素重新采样，以获得新的参考图像子区的灰度信息。

为避免参考像素在亚像素位置插值，文献[17]提出将图像子区中心平移到最邻近的整像素位置。由于计算量小、实现简单，且可避免插值误差，这一自适应图像子区平移策略在增量DIC 方法普遍采用。图1 为自适应图像子区平移策略示意图。以计算点 (x,y) 为例，在更新参考图像时，其在当前参考图中的匹配点位于 (xref,yref) ，如图1 中绿色箭头所示。在进行图像子区平移时，将其平移到最接近的整像素点，即 (round(xref), round(yref)) ，其中 round(x) 函数将x取整到最接近的整数。则当前计算点在平移后图像子区的局部坐标系中的坐标 (δx,δy) 可以表示为

图1 自适应图像子区平移策略示意图Fig.1 Schematic diagram of the adaptive subset offset strategy

若当前参考图像的匹配点 (xref,yref) 相对初始参考子区的累积位移为 (uref,vref) 。更新参考图像后，当前帧的计算结果为相对当前参考图像的结果，其匹配点位于 (xcurr,ycurr) ，如图1 中红色箭头所示。其相对当前参考子区的位移参数为{ucurr,vcurr,ux,vx,uy,vy}。另外，还需要结合形函数，得到平移子区的位移补偿量，如图1 中蓝色箭头部分所示。故在自适应图像子区平移策略中，当前帧的计算点相对初始参考图像的位移将由3 个部分构成。即当前参考图像的累积位移，当前帧相对当前参考图像的位移，以及图像子区平移的补偿位移。当前计算点相对初始参考图像的位移可由下式给出

1.2 增量DIC 计算中的参考图像更新策略

增量DIC 计算中更新参考图像的策略可以归纳为两类，即定步长更新策略和条件增量策略。其中，经典的连续更新策略可以视作步长为1 的定步长更新策略。如图2 所示，定步长增量策略在经过给定次数的计算后，即将当前变形图作为下一次相关计算的参考图像。而条件增量策略一直使用当前参考图进行计算，仅在后续计算中出现“退相关效应”时，将参考图像更新为上一帧图像。如果相邻的参考图之间相关系数较高，这些增量更新策略都能保证DIC 计算的顺利进行。

图2 增量DIC 方法中的定步长增量与条件增量策略示意图Fig.2 Schematic diagram of the fixed-step and the conditional strategies in incremental DIC

在增量DIC 的实际使用中，当图像变化较为平缓时，各图之间的相关系数较高，应该选取较大的增量步长。而当图像变化较快时，相关系数下降较快，则应选取较小的增量步长。故定步长的增量方法难以适应复杂多变的实际应用场景。而对于条件增量策略，为了确定更新参考图的时机，需要使用一定的判据。为此，Pan 等[3]选择使用计算区域中某个种子点的相关系数作为判据，当种子点的相关系数小于给定的阈值时，即可视为图像发生“退相关效应”，需使用上一帧图像作为新的参考图像进行后续计算。这一条件增量策略可以自适应图像的变化，并尽可能地减少更新次数，仅在需要的时候进行更新。由于更新参考图后需要对新的参考图的各个计算点重新预计算灰度梯度、Hessian 矩阵等，会带来额外的计算量，而各阶段的匹配误差在更新中也会不断累积，故在实际应用中，为了更好地使用增量DIC 方法，本文将通过实验验证来回答前文提出的两个问题，即：（1）更新参考图会不会带来额外的误差，以及这一误差的分布范围和分布规律是怎样的；（2）为了减少这一误差和计算量，应该采用何种策略更新参考图像。

2 实验与分析

2.1 模拟平移实验

在模拟平移实验中，本文选取了DIC challenge[18]中Sample 5 的参考图像作为基准，并使用快速傅里叶变换方法以0.1 像素为步长平移了20 步，共平移了2 个像素。通过一组模拟平移实验来评估增量算法在不同亚像素处进行更新参考图像带来的误差的累积规律与分布特性。

对于定步长增量策略，本文设计了如下的对比实验。其中使用传统序列DIC 的对照组，始终以初始参考图像为参考进行DIC 匹配。而使用增量DIC 的测试组，则以特定的步长进行更新参考图像。同时，选取了不同的计算步长（1～9 像素）以研究更新步长的影响。其中，在每组实验中均采用相同的计算点和计算参数。参考图像大小为512×512 像素，计算区域为边长460 像素的正方形区域，区域左上角像素坐标为 (30, 30)。以10像素为步长进行计算点采样，共有2209 个计算点。由于不同的插值方法会导致不同水平的插值误差，实验统一采用双三次B 样条插值作为图像灰度插值方式。

使用不同的计算方案完成计算后，统计了图像中各计算点的平均位移，并与真实位移做差，得到了各组实验的平均误差。将各组平移实验的平均误差绘制如图3 所示，其中横轴为亚像素平移量，纵轴为各组实验的平均误差。Seq 表示使用序列DIC 方法，当前图像相对初始参考图像的亚像素位移分布，而Inc 表示使用增量DIC 方法时，当前图像相对上一张参考图像的亚像素位移分布。图中黑色粗线为标准序列DIC 方法的误差计算结果，可以看到，序列DIC 方法的平均误差呈现出类似正弦函数的波动，这与已有DIC 插值误差的认识一致。而其余各组使用不同增量步长计算方案的结果均表现出对这一插值误差的累积现象，即每次更新参考图后，该参考图的插值误差会被累积到后续计算中，而最后的总误差即为各阶段的误差之和。这一结论是符合直觉的，增量算法不仅累积了位移，也同时累积了各个阶段的插值误差。这进一步说明了在使用增量DIC 算法时，应尽可能减少更新参考图的数量，仅在必要时更新参考图，以避免插值误差的累积。值得注意的是，从式（2）可知，增量DIC 的误差还包括使用形函数传递的补偿位移导致的误差，但这一误差相对亚像素插值误差的量级而言可以忽略不计，故不再赘述。

图3 采用不同增量步长的增量DIC 方法与传统序列DIC 方法的平均误差对比Fig.3 Comparison of average error between incremental DIC method with different step sizes and traditional sequential DIC method

2.2 真实平移实验

在实际实验中，各计算点的亚像素位移并非如模拟平移实验一样集中分布在给定位移附近，而是会存在一定的分布，故插值误差也会呈现出与模拟平移实验不同的分布。由前述分析可知，插值误差存在累积效应，且可以根据计算点的亚像素位移来估计其对应的插值误差。为研究真实平移实验的误差分布特性，对一固定在平移台上的平板进行平移实验。相机型号为FLIR-GS3-PGE-91QS6M，分辨率为3376×2704 像素，正对试样进行拍摄。实验中，试件沿着水平方向以0.2 mm 的步长进行了5 次平移。对于各组图像，分别采用序列DIC 方法和连续增量DIC 方法进行计算，并统计各计算点的位移分布（计算区域共有74×71 个计算点，计算步长为15 像素，图像子区为51×51 像素）。

在计算完成后，统计了两种方法得到的亚像素位移的分布，如图4 所示。可以看到序列DIC方法的亚像素位移随着平移量的增大，逐渐趋于均匀分布U(0,1)。这是因为在真实实验中，相机光轴与试样平面往往不是完全垂直，且存在镜头的放大率、畸变等多种因素的影响，经过一定的位移后，图像的亚像素位移将趋于随机的均匀分布。然而，在增量DIC 方法中，由于每次均更新参考图像，每一帧相对上一帧的位移较小，其增量位移的亚像素位移分布相对集中。图中可以看出，各次增量的位移主要分布在大于0.5 像素的区间，经过多次累积后会导致所得位移偏小。接下来，将通过前述的分析来估计其累积误差的分布。

图4 使用传统序列及增量DIC 方法在各个平移阶段的亚像素位移分布Fig.4 Subpixel displacement distributions obtained by traditional sequential and incremental DIC methods in real translation experiment

根据前述的模拟平移实验结果可知，亚像素位移的插值误差与亚像素位移之间的关系可用正弦函数来近似，而插值误差的量级与图像的噪声水平和所使用的插值函数有关。根据前述的实验结果，插值误差可以近似表示为

其中 du为亚像素位移部分，e为估计的插值误差。

对于序列DIC 方法，其插值误差可以直接通过式（3）进行估计。而对于增量DIC 方法，则需要将各增量步的插值误差进行累加，使用式（3）估计的插值误差的分布如图5 所示。由于正弦函数的非线性效应，序列DIC 方法中均匀分布的亚像素位移的插值误差呈现出两边高中间低的特性。而对于增量DIC 方法，由于每次更新参考图像均会带来一定的增量误差，故在累加后会得到一个范围更广的累积误差分布。另外，由图4 的结果可知，增量DIC 方法的亚像素位移分布较为集中且偏向于大于0.5 像素的区间。由式（3）计算得出其对应的插值误差会有明显的偏移，这导致累积后的插值误差相对序列DIC 方法呈现出较为明显的偏离。在经过4 次参考图像更新后，累积插值误差的最大值已经达到 0.02 像素，且有50%以上的计算点的插值误差大于0.01 像素。这一累积误差的量级已经足够体现出明显的差异，会进一步影响后续DIC 应变计算的精度。所以，在实际使用增量DIC 时，应尽可能地减少增量次数，以减少连续帧之间亚像素位移分布集中导致的累积误差。

图5 使用传统序列和增量DIC 方法的累积插值误差分布Fig.5 Cumulative error distributions obtained by traditional sequential and incremental DIC methods

为了降低这一累计误差的影响，Pan[19]提出使用高斯低通预滤波对图像进行预处理。这一预处理手段可以将插值的系统误差转换为随机误差，体现为位移结果的平均偏差降低，但标准差变大。文中指出，通过选取合适的高斯核大小，可以在不显著增加标准差的同时大幅降低偏差。为此，本研究对前述实验图片使用5×5 像素高斯核函数进行了高斯低通预滤波处理，并统计了增量DIC与序列DIC 方法之间的相对误差。为了直观展示二者之间的误差分布，图6 展示了误差的累积概率分布曲线，其中横轴为相对误差的分布范围，纵轴为对应的累积概率。其中实线为图像未经高斯滤波的结果，其相对误差较大且分散。另外，随着增量次数的增加，增量算法的最大相对误差也逐渐增大，经过4 次增量计算后最大相对误差约为0.04 像素。虚线为图像经过高斯预滤波之后的结果，经过高斯预滤波后，绝大多数点的相对误差均落在±0.01 像素之间，且增量算法并未导致明显的相对误差增加。结果证明，在使用增量DIC 算法时配合高斯低通预滤波算法，除了显著降低由于亚像素插值带来的误差之外，还可以降低由于增量算法带来的累积误差。

图6 使用传统序列和增量DIC 方法的相对误差的累积概率分布Fig.6 Cumulative probability distributions of the relative displacement error between traditional sequential and incremental DIC methods

2.3 橡胶拉伸大变形的增量DIC 测量

在材料测试实验中，增量DIC 常用于存在大变形的实验场景。为了展示增量DIC 算法在实际实验中的适用性，使用增量DIC 方法计算了一组橡胶材料的大变形拉伸实验图像。橡胶试样按ISO 37:2005 标准切割为标准哑铃型拉伸试件，其尺寸和外形如图7(a)所示，厚度为2 mm，并在表面使用喷漆制作了随机的散斑。如图7(b)所示，实验中使用一台SUNS-UTM5 万能实验机以2.5 mm/min 的速度匀速拉伸试件，最大拉伸位移为68 mm。实验中采用的相机型号为FLIRGS3-PGE-91QS6M（分辨率为3376×2704 像素，装有Schneider 23mm F/1.4 镜头），其固定在试件前方约1.2 m 处，正对试件表面进行采集，按每分钟15 帧的采集速率进行采集，实验全过程一共采集了407 帧图像。

图7 (a) 橡胶试件的尺寸与实际图像和(b) 实验现场装置Fig.7 (a) The dimensions and physical image of the rubber specimen and (b) the experimental setup

在匹配计算中，图像子区大小为41×41 像素，步长为5 像素，所有图像均经过了3×3 像素的高斯低通预滤波处理。由于橡胶试件变形过大，常规DIC 方法在初值估计阶段即会失效。即使采用时序传递策略（使用上一帧的结果作为初值），在变形过大时仍会由于形函数欠匹配、采样失真等原因导致相关系数下降，从而发生“退相关”。而采用增量算法，则可以顺利完成全部计算。在实验中采用了条件增量策略，即当计算区域中相关系数最低的计算点的相关系数低于0.95 时，则采用上一帧的图像作为参考图。这一条件判据相对使用特定种子点的判据[3]，更能体现全场的匹配质量，避免由于部分计算区域已“退相关”而种子点仍保持较高相关系数而带来的误判。所得部分计算结果如图8 所示，可见增量DIC 算法在经过多次增量后能得到实验全过程的变形，克服了传统DIC 方法在大变形时难以继续计算的问题。

图8 实验过程中橡胶试件沿拉伸方向的位移(单位：mm)Fig.8 The displacement of the rubber specimen in the stretching direction during the tensile test (unit: mm)

除此之外，我们还选取了不同的增量步长来计算同一组数据，以展示不同的增量步长在实际实验中的累积效应。在这里，我们设计了7 组定步长的测试组，其步长分别为1, 2, 4, 8, 16, 32,64 帧。图9(a) 展示了对实验的最后一帧，计算区域横向中心线上的数据点（如图8 中第0 帧所示），使用不同增量策略得到的位移结果。其中，中间的黑色线条为采用条件增量策略的增量方法，各计算点位移沿着拉伸方向表现出非常好的线性关系。在这一均匀拉伸实验中，各点位移应与其坐标成正比，故这一结果符合预期。而对于其他步长的结果，则表现出不同水平的累积误差，并沿着拉伸方向有周期性的变化。其中，在序号为50 的计算点处，各种步长的计算结果均相同。将这一点使用不同固定增量步长及使用条件增量步长得到的相关系数绘制如图9(b)所示，这里采用了零均值归一化互相关（zero-mean normalized cross correlation, ZNCC）相关系数。可见，在不更新参考图时，随着变形逐渐增大，相关系数会不断下降，而更新参考图后，相关系数又会回到1 附近。使用条件增量策略时，仅更新了4 次参考图，相比使用不同增量步长的其他实验组，条件增量策略可以尽可能少地进行增量计算，同时也能获得尽可能少的累积误差。除此之外，还统计了使用不同增量步长计算全部406 帧图像的计算时间，列于表1 中。可见，使用条件增量策略由于减少了增量的次数，可以尽可能地减少计算量，提高计算效率。

表1 使用不同固定增量步长及条件增量策略的计算耗时Table 1 Computation time of fixed-step strategy of different steps and conditional strategy

图9 (a) 不同位置的计算点使用不同更新策略下的位移和(b) 70 号计算点使用不同更新策略下全过程的ZNCC 相关系数Fig.9 (a) Displacements of different calculation points with different updating strategy and (b) ZNCC of #70 calculation point with different updating strategy during the experiment

值得注意的是，图9(a)中不同步长的累积误差在空间上呈现出周期性的波动，且其周期与增量步长成反比例关系。这是由于在这一实验中采用了稳定的加载速度和固定的采集帧率，各计算点的增量位移保持相对稳定。图10(a)和图10(b)分别展示了使用条件增量和连续增量策略时不同位置全过程的增量位移，即每一帧图像位移与上一帧图像的位移差。从图10(a)可以观察到，增量位移的大小也与其横向坐标成正比，这一结果符合均匀拉伸的假设。且序号为20 和序号为80的计算点的增量位移分别保持在0.5 像素和1.5像素附近，这表示间距为60 个计算点的两点的平均增量位移差约为1 像素，故这两点应具有相同的误差分布。这一现象与图9(a)中不同步长的增量结果出现周期性变化相符。由于匀速加载和拉伸试验的加持方式，增量位移与加载时间和空间坐标均呈线性关系，导致不同步长的累积误差呈现出周期性变化。其中采用连续增量策略的结果变化周期约为60 个计算点，而对于其他步长的增量算法，由于其间隔时间变长，两点的增量位移差会同比增长，而变化周期会同比缩短。

图10 沿拉伸方向不同位置计算点各帧的增量位移：(a) 使用条件增量策略, (b) 使用连续增量策略Fig.10 Incremental displacement of different calculation points along stretch direction with (a) conditional strategy and (b)sequential strategy

另外，对于序号为50 的计算点而言，其各帧的增量位移恰好为1 个像素。按前述的理论估计，该点的累积误差应该接近于0，这也与图9(a)的结果相符。另外，在图10(b)中也可看出，该点在连续增量中也保持稳定，但对于其他位置的计算点而言，均存在一定的计算误差。图中可以看出对于各位移均有向最近的整像素偏移的系统误差，而经过约400 次连续增量后，累积误差最大将近80 像素，如图9(a)中所示。这一结论再次证明了增量算法中插值误差的累积效应，也证明了条件更新策略的优越性。

3 结论

增量DIC 计算是解决大变形、环境光变化剧烈、散斑退相关等非常规场景中变形测量的有效手段。本文通过一系列的模拟及真实实验，分析了增量DIC 计算误差的理论累积规律，并通过真实实验展示了增量DIC 方法在实际实验中由亚像素位移分布所引起的累积误差的分布和形式，证明了高斯低通预滤波对这一累积误差的抑制效果。另外，还展示了增量DIC 算法在真实大变形实验中的应用，并通过对比实验证明了条件更新策略相比固定步长更新策略的优势。

根据上述研究，本文给出使用增量DIC 的如下建议。

（1）增量DIC 方法由于包含亚像素位移插值过程，存在稳定的系统误差，并会随着增量的过程不断累积。

（2）更新参考图像会不可避免地增加测量误差和计算量，应尽可能少地更新参考图像。在实际使用增量DIC 时，建议采用条件更新策略，如使用种子点或计算区域中相关系数最小的计算点的相关系数作为更新判据。

（3）更新参考图时，应采用随机图像子区平移策略将系统误差转化为随机误差，结合图像高斯低通预滤波处理，能显著降低由于插值带来的系统误差，以获得更好的误差抑制效果。

另外，增量DIC 算法不仅可用在容易发生退相关的测试环境，在3D–DIC、多相机DIC 等需要统一不同参考图中的同一计算点坐标系的场景也适用。在这些场景中，由于存在多个视角，故需要选定一主相机，其他相机的参考图均与主相机之间进行一次增量更新，再各自进行时序匹配。由于各相机时序匹配与2D 匹配过程相同，在多个相机之间的增量更新及后续的时序更新中，本文给出的建议仍然适用。