Introduction

The digital image correlation (DIC) method 1, as one of the modern noncontact photomechanics methods, has attracted extensive attention from scientific and engineering circles. According to the gray invariant hypothesis theory, the traditional DIC algorithm measures the displacement of each point in the digital speckle image of the object surface by matching the gray information of the reference image and the distorted image to obtain the deformation field of the object surface.

In the calculation and solution process of the DIC method, the integral pixel displacement is usually solved first, and then the subpixel displacement is solved [1, 2], forming an algorithm for solving the subpixel displacement represented by the Newton–Raphson algorithm [3, 4], gray gradient iterative algorithm [5], noniterative gray gradient algorithm [6], and surface fitting method [7]. Broggiato et al. [8] first proposed introducing the time-continuity constraint into the DIC algorithm and assumed that the acquisition time interval between each image in the image sequence was equal, thus realizing the continuity of displacements in the five images on the time axis. Besnard et al. [9, 10] combined the images collected in the experiment into a volume image with the time axis as the third dimension in the time sequence; this formed a movie-DIC considering time continuity, which improved the calculation accuracy of the DIC method. Volz et al. [11] proposed a novel parameterization method for multi-frame optical flow computation, this method can naturally embed the assumption of temporally coherent spatial flow structure, and the assumption that the optical flow is smooth along motion trajectories. To address the shortcomings of the above algorithms in time, Xian Wang [12] introduced the continuity assumption of solid deformation into the DIC algorithm and proposed a DIC model considering spatiotemporal continuity, namely, the space–time window DIC method. This algorithm averages the deformed images in a short period of time, establishes the displacement representation mode based on the first-order function relationship, and obtains a low signal-to-noise ratio image deformation analysis algorithm (STS-DIC). Starting from the algorithm side, Xuejin Liu [13] improved the space–time window DIC method and proposed an image analysis algorithm based on a quadratic function model to represent continuous deformation. Xiaoyong Liu [14] proposed a nonlinear optical flow equation to describe the variation in speckle gray level before and after deformation under the influence of unstable or nonuniform illumination, excessive or insufficient exposure of the camera, and surface deformation of the measured object.

After continuous basic theoretical research, the DIC method has also been applied in practice, especially in some severe environments. Roux S. and Hild F. [15] used the DIC method to study crack propagation in brittle materials, achieving nanometer-level displacement measurement accuracy and 7% error value in stress intensity factor calculations. Avril S. [16] and Sutton M. A. [17] measured conventional materials using the DIC method. Hwang S. F. [18], Van Paepegem W. [19], and others applied the DIC method in new materials. Lee D. et al. [20] studied fracture behavior of a graphite/epoxy composite using the DIC method and a high-speed camera. Sztefek P. et al. [21] tracked speckle patterns on bone surfaces during loading to determine surface strain. Yoneyama S. [22] measured bridges and Pilch A. et al. [23] measured surfaces of various vehicles. The DIC method was also used in harsh environments, such as fatigue tests of fiber-matrix interfaces [24] by Hao Wenfeng et al. and soil deformation measurements [25] by Liu Changbo. In high-speed rolling mechanics, the DIC method was used to measure automobile tire properties [26] and calculate displacement fields near a precast crack [27]. The DIC method was applied to study fracture characteristics of layered slate [28] by Li Erqiang, underwater propeller deformation [29] by Pan Jiyong, and ice fracture performance [30] by Wang Juan. Weina Guo et al. [31] used the DIC method to measure and calculate the influence of fly ash content on the mechanical properties of cement-based composites. Hartmann [32] proposed a spatiotemporal optical flow method and used it to determine the correlations between the microstructure of the material and macroscopic process results.

However, to date, most scholars have mainly utilized the DIC method in the frequency domain and the spatial domain, and research on time-domain DIC has been relatively limited. The time-domain algorithm is still in use as a displacement mode, and based on the linear deformation of time, the displacement mode has some deficiencies. For example, it is not applicable to complex deformation, only for intensive acquisition of image processing and so on. Therefore, it is necessary to develop a high-order time-domain DIC algorithm model. Aiming at the problems in the above algorithms, a high-order time-domain DIC algorithm based on the nonlinear optical flow equation is proposed by introducing the time-dependent nonlinear displacement mode into the optical flow equation. The solution and precision of the algorithm are analyzed for the subpixel displacement of low-quality images.

High-order Time-domain DIC Algorithm for Nonlinear Optical Flow Equations

High-Order DIC Algorithm in the Time Domain

A high-order time-domain DIC algorithm model is established based on the nonlinear optical flow equation. Suppose the gray matrix of the image before and after deformation is f(x) and g(x). Meanwhile, it is assumed that the gray value of any point in the speckle pattern presents a quadratic nonlinear relationship before and after deformation

$$a{f}^{2}(X)+bf(X)+c=g(X+U(X))$$
(1)

where \(X\text{=}\left[\begin{array}{c}x\\ y\end{array}\right]=\left[\begin{array}{ccccc}{x}_{1}& {x}_{2}& {x}_{3}& ...& {x}_{N}\\ {y}_{1}& {y}_{2}& {y}_{3}& ...& {y}_{N}\end{array}\right]\) is the coordinate of the reference image, \(U(X)=\left[\begin{array}{c}u\\ v\end{array}\right]=\left[\begin{array}{ccccc}{u}_{1}& {u}_{2}& {u}_{3}& ...& {u}_{N}\\ {v}_{1}& {v}_{2}& {v}_{3}& ...& {v}_{N}\end{array}\right]\) is the displacement vector, and N is the number of pixels in the reference image. The undetermined coefficients in equation (1) are a, b, and c. When a = 0, b = 1, and c = 0, equation (1) is reduced to traditional optical flow setup [33], which is also a special case of the model in this Section. Equation (1) degenerates into the basic assumption of an invariable gray level in image motion estimation, which is also a special case in the model.

During the experiment, a charge-coupled device (CCD) camera continuously collects a speckle image sequence. Assuming that the sequence contains n + 1 images, the speckle image collected at time t = 0 is taken as the reference image, and the image collected at time t = ti (i = 1, 2…, n) is a deformed image, the nonlinear optical flow equation of the speckle image at time t can be expressed as:

$$a{f}^{2}(X)+bf(X)+c={g}_{t}(X+U(X,t))$$
(2)

where U (X, t) is the displacement matrix of the deformed image collected at time t relative to the reference image. Then, Eq. (1) has been connected to the time domain.

The image between any time period [ti-m, ti +m] is selected on the time axis of the speckle image sequence, and a total of 2 m + 1 (m = 1,2,3…) speckle images are taken in this time period. According to equation (2), 2 m + 1 nonlinear optical flow equations are established, and they can be obtained by integration:

$$a{f}^{2}(X)+bf(X)+c=\frac{1}{2m+1}{\int }_{{t}_{i-m}}^{{t}_{i+m}}{g}_{t}(X+U(X,t))dt$$
(3)

The right side of equation (3) is the average of the above 2 m + 1 deformation images, which can be regarded as a new deformation speckle image, expressed as:

$${\overline{G} }_{i}{}^{(m)}(X+U(X,t))=\frac{1}{2m+1}{\int }_{{t}_{i-m}}^{{t}_{i+m}}{g}_{t}(X+U(X,t))dt$$
(4)

where U (X, t) is the displacement matrix to be solved, and the independent variables are X and t.

If the deformation images in the experiment are collected continuously and the specimen has not been damaged, then U (X, t) can be considered a continuous function of time t. Therefore, U (X, t) carries out a second-order Taylor expansion on time t and ignores small quantities of higher order, which is obtained:

$$U(X,t)=U(X,t){|}_{t={t}_{i}}+\frac{\partial U(X,t)}{\partial t}{|}_{t={t}_{i}}\Delta t+\frac{1}{2!}\frac{{\partial }^{2}U(X,t)}{\partial {t}^{2}}{|}_{t={t}_{i}}{(\Delta t)}^{2}$$
(5)

where \(\Delta t=t-{t}_{i}\) is the interval between the acquisition time of two images. The first term is the traditional displacement shape function, which usually considers the affine deformation of the relevant window. The first-order Taylor expansion of the central point coordinate, XC, can be written as follows:

$$U(X,t){|}_{t={t}_{i}}=\left(\begin{array}{c}u\\ v\end{array}\right){|}_{\begin{array}{c}t={t}_{i}\\ X={X}_{c}\end{array}}+\left(\begin{array}{cc}\frac{\partial u}{\partial x}& \frac{\partial u}{\partial y}\\ \frac{\partial v}{\partial x}& \frac{\partial v}{\partial y}\end{array}\right){|}_{\begin{array}{c}t={t}_{i}\\ X={X}_{c}\end{array}}\left(\begin{array}{c}\Delta x\\ \Delta y\end{array}\right)$$
(6)

The second term is the deformation velocity term, which can be obtained after first-order Taylor expansion:

$$\frac{\partial U(X,t)}{\partial t}{|}_{t={t}_{i}}=\left(\begin{array}{c}\frac{\partial u}{\partial t}\\ \frac{\partial v}{\partial t}\end{array}\right){|}_{\begin{array}{c}t={t}_{i}\\ X={X}_{C}\end{array}}+\left(\begin{array}{cc}\frac{{\partial }^{2}u}{\partial x\partial t}& \frac{{\partial }^{2}u}{\partial y\partial t}\\ \frac{{\partial }^{2}v}{\partial x\partial t}& \frac{{\partial }^{2}v}{\partial y\partial t}\end{array}\right){|}_{\begin{array}{c}t={t}_{i}\\ X={X}_{c}\end{array}}\left(\begin{array}{c}\Delta x\\ \Delta y\end{array}\right)$$
(7)

The third term is the deformation acceleration term, which is expanded as follows:

$$\frac{{\partial }^{2}U(X,t)}{\partial {t}^{2}}{|}_{t={t}_{i}}=\left(\begin{array}{c}\frac{{\partial }^{2}u}{\partial {t}^{2}}\\ \frac{{\partial }^{2}v}{\partial {t}^{2}}\end{array}\right){|}_{\begin{array}{c}t={t}_{i}\\ X={X}_{C}\end{array}}+\left(\begin{array}{cc}\frac{{\partial }^{3}u}{\partial x\partial {t}^{2}}& \frac{{\partial }^{3}u}{\partial y\partial {t}^{2}}\\ \frac{{\partial }^{3}v}{\partial x\partial {t}^{2}}& \frac{{\partial }^{3}v}{\partial y\partial {t}^{2}}\end{array}\right){|}_{\begin{array}{c}t={t}_{i}\\ X={X}_{C}\end{array}}\left(\begin{array}{c}\Delta x\\ \Delta y\end{array}\right)$$
(8)

Substitute equations (6), (7) and (8) into equation (5), keep the second derivative term and ignore the higher-order term. Then, the displacement of any point in the relevant window at time t is:

$$\begin{array}{c}U(X,t)=U(X,t;{k}_{0})=\left(\begin{array}{c}u\\ v\end{array}\right){|}_{\begin{array}{c}t={t}_{i}\\ X={X}_{C}\end{array}}+\left(\begin{array}{cc}\frac{\partial u}{\partial x}& \frac{\partial u}{\partial y}\\ \frac{\partial v}{\partial x}& \frac{\partial v}{\partial y}\end{array}\right){|}_{\begin{array}{c}t={t}_{i}\\ X={X}_{C}\end{array}}\left(\begin{array}{c}\Delta x\\ \Delta y\end{array}\right)+\left(\begin{array}{c}\frac{\partial u}{\partial t}\\ \frac{\partial v}{\partial t}\end{array}\right){|}_{\begin{array}{c}t={t}_{i}\\ X={X}_{C}\end{array}}\left(\Delta t\right)+\\ \left(\begin{array}{cc}\frac{{\partial }^{2}u}{\partial x\partial t}& \frac{{\partial }^{2}u}{\partial y\partial t}\\ \frac{{\partial }^{2}v}{\partial x\partial t}& \frac{{\partial }^{2}v}{\partial y\partial t}\end{array}\right){|}_{\begin{array}{c}t={t}_{i}\\ X={X}_{C}\end{array}}\left(\begin{array}{c}\Delta x\Delta t\\ \Delta y\Delta t\end{array}\right)+\frac{1}{2!}\left(\begin{array}{c}\frac{{\partial }^{2}u}{\partial {t}^{2}}\\ \frac{{\partial }^{2}v}{\partial {t}^{2}}\end{array}\right){|}_{t={t}_{i}}{(}^{\Delta }={k}_{0}{\left(\begin{array}{ccccccc}1& \Delta x& \Delta y& \Delta t& \Delta x\Delta t& \Delta y\Delta t& \frac{\Delta {t}^{2}}{2!}\end{array}\right)}^{T}\end{array}$$
(9)

and:

$${k}_{0}\text{=}{\left.\left[\begin{array}{ccccccc}u& \frac{\partial u}{\partial x}& \frac{\partial u}{\partial y}& \frac{\partial u}{\partial t}& \frac{{\partial }^{2}u}{\partial x\partial t}& \frac{{\partial }^{2}u}{\partial y\partial t}& \frac{{\partial }^{2}u}{\partial {t}^{2}}\\ v& \frac{\partial v}{\partial x}& \frac{\partial v}{\partial y}& \frac{\partial v}{\partial t}& \frac{{\partial }^{2}v}{\partial x\partial t}& \frac{{\partial }^{2}v}{\partial y\partial t}& \frac{{\partial }^{2}v}{\partial {t}^{2}}\end{array}\right]\right|}_{\begin{array}{c}t={t}_{i}\\ x={x}_{c}\end{array}}$$
(10)

Since in the parameter matrix, ΔxΔt and ΔyΔt are functions of coordinate and time parameters, rather than independent parameters, the final form of equations (3) and (4) is:

$${\overline{G} }_{i}{}^{(m)}(X+U(X,t;k))=a{f}^{2}(X)+bf(X)+c=\frac{1}{2m+1}{\int }_{{t}_{I-m}}^{{t}_{i+m}}{g}_{t}[X+k{\left(\begin{array}{ccccc}1& \Delta x& \Delta y& \Delta t& \frac{{(\Delta t)}^{2}}{2}\end{array}\right)}^{T}]dt$$
(11)
$$k=\left(\begin{array}{ccccc}u& \frac{\partial u}{\partial x}& \frac{\partial u}{\partial y}& \frac{\partial u}{\partial t}& \frac{{\partial }^{2}u}{\partial {t}^{2}}\\ v& \frac{\partial v}{\partial x}& \frac{\partial v}{\partial y}& \frac{\partial v}{\partial t}& \frac{{\partial }^{2}v}{\partial {t}^{2}}\end{array}\right){|}_{\begin{array}{c}t={t}_{i}\\ x={x}_{c}\end{array}}$$
(12)

Solution of the High-Order DIC Algorithm in the Time Domain

Change equation (11) into:

$${\text{Q}}_{i}(P)=\frac{1}{2m+1}{\int }_{{t}_{i-m}}^{{t}_{i+m}}{g}_{t}[X+k{\left(\begin{array}{ccccc}1& \Delta x& \Delta y& \Delta t& \frac{{(\Delta t)}^{2}}{2}\end{array}\right)}^{T}]dt-a{f}^{2}(X)-bf(X)-c\cong 0$$
(13)

The parameter in equation (13) to be solved is:

$$\begin{array}{c}P={\left(\begin{array}{cccc}k& a& b& c\end{array}\right)}^{T}\\ \, \, \, \text{=}{\left(\begin{array}{ccccccccccccc}u& \frac{\partial u}{\partial x}& \frac{\partial u}{\partial y}& \frac{\partial u}{\partial t}& \frac{{\partial }^{2}u}{\partial {t}^{2}}& v& \frac{\partial v}{\partial x}& \frac{\partial v}{\partial y}& \frac{\partial v}{\partial t}& \frac{{\partial }^{2}v}{\partial {t}^{2}}& a& b& c\end{array}\right)}^{T}\end{array}$$
(14)

Each pixel in the image subset (i.e., i = 1, 2…, n) can be obtained, as shown in equation (13). These equations form a statically indeterminate system of equations, and using the least squares method to solve it, the least squares iterative solution of equation (15) can be obtained:

$${P}^{l+1}={P}^{l}-{[\nabla Q{({P}^{l})}^{\text{T}}\nabla Q({P}^{l})]}^{-1}\nabla Q{({P}^{l})}^{\text{T}}\nabla Q({P}^{l})$$
(15)

and

$$\nabla Q({P}^{l})\text{=}\left(\begin{array}{ccccccccccccc}{\overline{G} }_{{x}_{1}}{}^{(m)}& {\overline{G} }_{{x}_{1}}{}^{(m)}\Delta {x}_{1}& {\overline{G} }_{{x}_{1}}{}^{(m)}\Delta {y}_{1}& {\overline{G} }_{{x}_{1}}{}^{(m)}\Delta {t}_{1}& \frac{{\overline{G} }_{{x}_{1}}{}^{(m)}{\left(\Delta {t}_{1}\right)}^{2}}{2}& {\overline{G} }_{{y}_{1}}{}^{(m)}& {\overline{G} }_{{y}_{1}}{}^{(m)}\Delta {x}_{1}& {\overline{G} }_{{y}_{1}}{}^{(m)}\Delta {y}_{1}& {\overline{G} }_{{y}_{1}}{}^{(m)}\Delta {t}_{1}& \frac{{\overline{G} }_{{y}_{1}}{}^{(m)}{\left(\Delta {t}_{1}\right)}^{2}}{2}& -{f}_{1}{}^{2}& -{f}_{1}& -1\\ {\overline{G} }_{{x}_{2}}{}^{(m)}& {\overline{G} }_{{x}_{2}}{}^{(m)}\Delta {x}_{2}& {\overline{G} }_{{x}_{2}}{}^{(m)}\Delta {y}_{2}& {\overline{G} }_{{x}_{2}}{}^{(m)}\Delta {t}_{2}& \frac{{\overline{G} }_{{x}_{2}}{}^{(m)}{\left(\Delta {t}_{2}\right)}^{2}}{2}& {\overline{G} }_{{y}_{2}}{}^{(m)}& {\overline{G} }_{{y}_{2}}{}^{(m)}\Delta {x}_{2}& {\overline{G} }_{{y}_{2}}{}^{(m)}\Delta {y}_{2}& {\overline{G} }_{{y}_{2}}{}^{(m)}\Delta {t}_{2}& \frac{{\overline{G} }_{{y}_{2}}{}^{(m)}{\left(\Delta {t}_{2}\right)}^{2}}{2}& -{f}_{2}{}^{2}& -{f}_{2}& -1\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ {\overline{G} }_{{x}_{n}}{}^{(m)}& {\overline{G} }_{{x}_{n}}{}^{(m)}\Delta {x}_{n}& {\overline{G} }_{{x}_{n}}{}^{(m)}\Delta {y}_{n}& {\overline{G} }_{{x}_{n}}{}^{(m)}\Delta {t}_{n}& \frac{{\overline{G} }_{{x}_{n}}{}^{(m)}{\left(\Delta {t}_{n}\right)}^{2}}{2}& {\overline{G} }_{{y}_{n}}{}^{(m)}& {\overline{G} }_{{y}_{n}}{}^{(m)}\Delta {x}_{n}& {\overline{G} }_{{y}_{n}}{}^{(m)}\Delta {y}_{n}& {\overline{G} }_{{y}_{n}}{}^{(m)}\Delta {t}_{n}& \frac{{\overline{G} }_{{y}_{n}}{}^{(m)}{\left(\Delta {t}_{n}\right)}^{2}}{2}& -{f}_{n}{}^{2}& -{f}_{n}& -1\end{array}\right)$$
(16)

Pl and Pl + 1 are the l th and l + 1 th iterative solutions, respectively.

Since the displacement gradient value and acceleration in practice are much smaller than the displacement value, and in the case of the slow loading process and the firmly mounted specimen, the initial value of the displacement gradient is usually set to 0. In the iterative solution process, the initial value vector is set as:

$${P}^{0}={\left[\begin{array}{ccc}{u}_{0}& 0& \begin{array}{ccc}0& 0& \begin{array}{ccc}0& {v}_{0}& \begin{array}{ccc}0& 0& \begin{array}{ccc}0& 0& \begin{array}{ccc}0& 1& 0\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\right]}^{T}$$
(17)

Before iteration, the bicubic spline interpolation method is used to obtain the gray values and the first-order gray gradients of the subpixel positions on the transformed subsets. After iteration and convergence, three sets of data are obtained, which are the displacement at the center of the reference subset (x0, y0), the displacement gradient in the x and y directions, and the first and second partial derivatives of time t. The above iterative process is adopted for each pixel in the reference image to obtain the entire displacement field of the speckle pattern.

Speckle Pattern Making

The experiment is carried out by manually spraying a speckle pattern (Fig. 1) and applying image displacement by the MATLAB algorithm [33].

Fig. 1
figure 1

Speckle patterns

Full size image

To simulate the noise generated by the CCD camera in the shot scene, Gaussian white noise with standard deviation σg = 5 is added to the images.

To simulate the influence of ambient brightness changes on the shooting scene, seven groups of images with different gray levels are generated from the original images to analyze the influence of ambient brightness changes on the result errors.

By adjusting the contrast and brightness of the translated image, seven groups of target images with different gray levels are obtained to simulate the changes in the speckle pattern gray level caused by the changes in the test site environment and the imaging system aperture. The specific adjustment methods are as follows:

  1. (1)

    Increase the brightness of the translated speckle pattern by 20% and gradually increase the brightness by 20% until it reaches 140% to obtain 7 groups of distorted images with different brightnesses, denoted as the first to the seventh gray levels.

  2. (2)

    After the gray values of the above 7 groups of deformed images are normalized, the contrast stretching is processed according to equation (18).

    $$\begin{array}{cc}{\text{I}}_{\text{out}}\text{=}& \left\{\begin{array}{ccc}\begin{array}{c}0\\ \left[\left({I}_{in}-0.05\right)/0.9\right]\\ 1\end{array}& \begin{array}{c}{\text{if}}\\ {\text{if}}\\ {\text{if}}\end{array}& \begin{array}{c}{I}_{in}\le 0.05\\ 0.05<{I}_{in}<0.95\\ {I}_{in}\ge 0.95\end{array}\end{array}\right.\end{array}$$
    (18)
  3. (3)

    After contrast stretching, the deformed images are reverse-normalized to obtain 7 groups of deformed images with different levels of gray change.

The figure below shows the deformation speckle pattern of the seventh gray level and its gray distribution. Compared with the reference image, most of the pixel points of the deformation speckle pattern of the seventh gray level are concentrated at the point where the gray level is 0 and the gray level is 255, with significant polarization (Figs. 2 and 3).

Fig. 2
figure 2

Reference speckle image’s gray level distribution

Full size image
Fig. 3
figure 3

The 7th grayscale deformation speckle image’s gray level distribution

Full size image

Algorithm Verification Based on the Simulation

The high-order time-domain algorithm is used to select 361 points from the reference image for calculation. The relevant window size is 19 × 19 pixels, and the calculation step size is 13 pixels. Equation (19) is used to calculate the overall standard deviation of the displacement:

$$e=\sqrt{\frac{1}{361}\sum_{k=1}^{361}{\left({u}_{k}-{\widetilde{u}}_{k}\right)}^{2}}$$
(19)

where \({u}_{k}\) and \({\widetilde{u}}_{k}\) are the calculated displacements and theoretical displacements, respectively.

Test of the Linear Displacement Field in the Time Domain

In the experiment, the simulated speckle image is taken as the reference image. Given the displacement field of rigid body translation of the speckle image sequence, 9 deformation images with different displacements are obtained. In this analysis, t = 5 is taken as the central image of the calculation, and t = 5 ± m (m = 1,2,3,4) is taken as the auxiliary calculation image. The central image and the auxiliary image are taken as the calculated deformation image through equation (3). The relation between image displacement u and time t is as follows:

$$\begin{array}{cc}u={a}_{1}t,& ({a}_{1}=0.05\text{pixel/s},t=1{\text{s}},2{\text{s}},...,9{\text{s}})\end{array}$$
(20)

Figure 4 shows the displacement measurement error curve changing with the number of image sequences. As the statistical curves show, the minimum error value of the overall horizontal displacement error when the number of image sequences is 5, so 5 is the optimal value of the window.

Fig. 4
figure 4

The displacement measurement error curve with the change in the number of image sequences

Full size image

Figure 5 represents the displacement error curve with the change in gray level when the number of image sequences is 1, 5, and 9. Under the condition of a certain number of image sequences, the error tends to rise with increasing gray level. The error is on the order of 10–3 with increasing gray level and will not exceed 1.6 × 10–2.

Fig. 5
figure 5

Displacement measurement error curve of the linear displacement with the change in gray level

Full size image

Test of the Quadratic Displacement Field in Time Domain

The quadratic displacement with respect to time as listed in Eq. (21) is applied to the simulated speckle image sequence:

$$\begin{array}{cc}u={a}_{2}{t}^{2},& ({a}_{2}=0.01{\text{pixel/s}}^{2},t=1{\text{s}},2{\text{s}},...,9{\text{s}})\end{array}$$
(21)

The brightness adjustment method in 2.1 is used to obtain 7 groups of distorted images with different brightnesses. The high-order time-domain DIC method is used to calculate the calculation area selected in 2.1, and the statistical errors are obtained, as shown in Fig. 6.

Fig. 6
figure 6

The displacement measurement error curve with the change in the number of image sequences

Full size image

As seen from Fig. 6, with an increase in the number of image sequences, the overall horizontal displacement error is the minimum at 5, so the calculation effect is the best when the number of image sequences is 5, which is the critical size with the best calculation effect. The three curves in Fig. 7 represent the displacement error curves with the change in gray level when the number of image sequences is 1, 5, and 9. Under the condition of a certain number of image sequences, the error tends to rise with increasing gray level.

Fig. 7
figure 7

Displacement measurement error curve of the quadratic displacement with the change in gray level

Full size image

Figure 8 compares the curve of the horizontal displacement error measured by the traditional spatial domain algorithm [34], first-order time-domain algorithm and second-order time-domain algorithm with gray level. It can be seen from the figure that with the increase in the change level of grayscale, the average measurement displacement error of the traditional algorithm increases sharply, while the variation in the measurement displacement error of the high-order time-domain DIC algorithm also increases but very slowly. When the gray change level is 140%, the average measurement error of the traditional algorithm is as high as 0.055 pixels, which is more than 20% of the real displacement (which is 0.25 pixels). However, the measurement error of the high-order time-domain DIC algorithm in linear displacement mode is 0.015 pixels, and the measurement error of the quadratic displacement mode is less than 0.005 pixels, which is far less than the measurement error of the traditional algorithm. The analysis shows that the high-order time-domain DIC algorithm is not sensitive to changes in illumination and other external conditions.

Fig. 8
figure 8

Displacement measurement error curves under the two algorithms

Full size image

Test of Fixed Axial Compression

In the following experiment, the displacement of each point is related to its own coordinates, so a plane Cartesian coordinate system with the pixel point in the lower left corner of the picture as the origin is established on the picture. The schematic diagram of the coordinate system is as follows (Fig. 9):where the coordinate in the lower left corner of the picture is (0,0), and the coordinate in the upper right corner is (255, 255).

Fig. 9
figure 9

Cartesian coordinate system established on the picture

Full size image

Transverse compression is applied to the simulated speckle image sequence as follows:

$$\begin{array}{cc}\Delta d=d'-d=-a_3t,&(a_3=1\text{pixel/s},t=1\text{s},2\text{s},...,9\text{s})\end{array}$$
(22)

where d is the width of the image before compression, d' is the width of the image after compression, and Δd is the change in the width of the image after compression.

Different from the previous two experiments, the displacement of each column of the picture is different in the fixed-axis compression experiment, and its displacement increases with the increase in the abscissa in the original picture. The displacement formula is as follows:

$$u(x,t)=\frac{\Delta dx}{d}=\frac{-tx}{d}$$
(23)

where x is the abscissa of the pixel point in the coordinate system. After displacement, the schematic diagram is shown in the following figure (Fig. 10):

Fig. 10
figure 10

Schematic diagram of fixed axial compression

Full size image

The displacement nephogram when t = 5 is shown in the following figure (Fig. 11):

Fig. 11
figure 11

Displacement nephogram when t = 5

Full size image

The same method is used to obtain 7 groups of deformed images with different brightness for calculation, and the statistical errors are obtained, as shown in the figure below (Figs. 12 and 13):

Fig. 12
figure 12

The displacement measurement error curve with the change in the number of image sequences

Full size image
Fig. 13
figure 13

Displacement measurement error curve of the fixed axial compression experiment with the change in gray level

Full size image

It can be seen that in the fixed-axis compression experiment, different from the traditional algorithm, the displacement error value is stable within 0.050 pixels with increasing gray level. This indicates that this method has a high accuracy in the calculation of the fixed-axis compression displacement.

Test of Shear Deformation

Shear deformation is applied to the simulated speckle image sequence, and the effect is shown in the figure below (Fig. 14):where the displacement of the pixel at the top of the image (y = 255) is:

Fig. 14
figure 14

Schematic diagram of shear deformation

Full size image
$$\begin{array}{cc}u(255,t)=-{a}_{4}t,& ({a}_{4}=0.5\text{pixel/s},t=1{\text{s}},2{\text{s}},...,9{\text{s}})\end{array}$$
(24)

The displacement of any point in the image is related to its y-coordinate, and the displacement formula is:

$$u(y,t)=\frac{-0.5yt}{d}$$
(25)

where d is the width of the image and t is the moment when deformation occurs.

The same method is used to obtain 7 groups of deformed images with different brightnesses, and the DIC method is used for calculation. The statistical errors are obtained, as shown in the figure below (Figs. 15 and 16):

Fig. 15
figure 15

The displacement measurement error curve with the change in the number of image sequences

Full size image
Fig. 16
figure 16

Displacement measurement error curve of the constant axial compression experiment with the change in gray level

Full size image

In the shear deformation experiment, if the number of image sequences is appropriate, then the displacement error value of the solution tends to be stable with increasing gray level. When the number of image sequences is 5, the error can be stabilized within 0.037 pixels, indicating that this method also has a high accuracy in the calculation of shear deformation.

Test of Fixed-Point Rotation

The rotation quantity with respect to time, as listed in equation (26), is applied to the sequence of simulated speckle images.

$$\begin{array}{cc}\theta ={a}_{5}t,& ({a}_{5}=0.24^\circ /{\text{s}},t=1{\text{s}},2{\text{s}},...,9{\text{s}})\end{array}$$
(26)

where the rotation center of the image is the midpoint of the image, and the rotation direction is counterclockwise.

Different from the previous four experiments, each pixel in the fixed-point rotation experiment has not only axial displacement but also longitudinal displacement. Here, the axial displacement is related to the y-coordinate, the longitudinal displacement is related to the x-coordinate, and the relation is as follows:

$$\begin{array}{c}u(y,t)=(\frac{d-1}{2}-y)(\frac{1-\mathrm{cos}\theta }{\mathrm{tan}\theta })=(\frac{d-1}{2}-y)[\frac{1-\mathrm{cos}(0.24^\circ t)}{\mathrm{tan}(0.24^\circ t)}]\\ v(x,t)=(\frac{d-1}{2}-x)(\frac{1-\mathrm{cos}\theta }{\mathrm{tan}\theta })=(\frac{d-1}{2}-x)[\frac{1-\mathrm{cos}(0.24^\circ t)}{\mathrm{tan}(0.24^\circ t)}]\end{array}$$
(27)

It can be seen that the closer the horizontal (or vertical) coordinate of the pixel is to the midpoint in the image, the smaller its vertical (or horizontal) displacement is. The vertical (or horizontal) displacement can be regarded as the first-order function of the horizontal (or vertical) coordinate, which is:

$$\begin{array}{c}u(y,t)={k}_{1}y+{b}_{1}\\ v(x,t)={k}_{2}y+{b}_{2}\end{array}$$
(28)

After calculating the displacement of each point in the image, the rotation angle of each point can be calculated by the following formula [35]:

$$\theta =\frac{1}{2}[\frac{\partial v(x,t)}{\partial x}-\frac{\partial u(y,t)}{\partial y}]=\frac{1}{2}({k}_{2}-{k}_{1})=\frac{90^\circ}{\pi }({k}_{2}-{k}_{1})$$
(29)

where k1 and k2 can be obtained by fitting the calculated displacement.

The same method is used to obtain 7 groups of deformed images with different brightnesses, and the DIC method is used for calculation. The rotation angle of each point is calculated by the above formula, and the error is obtained as shown in the figure (Figs. 17 and 18):

Fig. 17
figure 17

Rotation angle measurement error curve with the change in the number of image sequences

Full size image
Fig. 18
figure 18

Rotation angle error curve of the fixed-point rotation experiment with a change in gray level

Full size image

The error of the algorithm in calculating the fixed-point rotation angle is stable below 0.001° when the error of the traditional algorithm is more than 0.01°.

Chapter Summary

It can be seen from Figs. 4, 6, 12, 15 and 17 that even when Gaussian noise is added, when the brightness of the deformed image changes, the increase of calculation error is still small, indicating that the algorithm can effectively resist the influence caused by environmental changes.

And from Fig. 8, it can be seen that traditional DIC method is more susceptible to brightness changes. Therefore, the light changes and/or other external conditions have little effect on the application of the proposed method.

Algorithm Verification Based on Actual Experiment

Experimental Design

A resin mode I crack experiment is carried out. In this experiment, the same batch of photosensitive resin materials is used to prefabricate the Brazilian half-disk specimen with a 90° incision. The fracture experiment is conducted under three-point bending loading conditions. The specimen is cut into disks with a thickness B = 20 mm on a cutting machine, and the disks are cut into two semidisk specimens along their diameters. Additionally, a 0.5 mm thick circular fine steel stone saw blade is used to cut a 10 mm long precast notch from the center of the half disk along the direction of the upper loading point. The specific size of the specimen is as follows: radius R = 60 mm, crack length a = 10 mm, thickness B = 10 mm, span 2S = 50 mm, and height ratio a/R = 0.17. The specimen size is shown in the figure below (Fig. 19).

Fig. 19
figure 19

Schematic diagram of the mode I crack experiment

Full size image

A 100 kN MTS CMT5105 electronic servo testing machine provides uniform loading through displacement control, with a loading speed of 0.02 mm/min. To ensure close contact between the loading system and the specimen, a 20 N prestress is applied before formal loading. Prior to the experiment, the specimen is numbered, and artificial speckle is sprayed on the front of the specimen. To improve the speckle contrast on the specimen surface, the specimen is first sprayed white with self-spray paint, and then random black spots are sprayed on the surface after the white paint is dry. Speckle images of the crack tip area on the specimen surface in the loading process are captured by a 404 K array black-and-white industrial CCD camera produced by Germany Basler company, and the resolution of the image is 1280 × 1024 pixels. The object plane resolution of the collected speckle field on the specimen surface is 0.050 mm/pixel, and the image acquisition rate is 15 fps. Then, the traditional DIC method and high-order DIC time-domain algorithm are used to process the data.

Figure 20 is the image of the specimen, the highlighted part is the region of interest (ROI).

Fig. 20
figure 20

Image of the specimen taken during the experiment

Full size image

Data Processing

Same as Sect. "Algorithm Verification Based on the Simulation", the relevant window size is 19 × 19 pixels, and the calculation step size is 13 pixels.

The traditional DIC method and high-order time-domain DIC method are used to calculate the displacement fields of sample images at 100, 150 and 200 s. The results are as follows:

In Figs. 21, 22 and 23, (a) and (b) exhibit the nephograms of the displacement calculation results of the traditional DIC algorithm and the higher-order time-domain DIC algorithm, respectively.

Fig. 21
figure 21

Lateral displacement nephogram at an experiment time of 100 s

Full size image
Fig. 22
figure 22

Lateral displacement nephogram at an experiment time of 150 s

Full size image
Fig. 23
figure 23

Lateral displacement nephogram at an experiment time of 200 s

Full size image

Compared with the displacement field results obtained by the traditional DIC method, it can be seen:

  1. 1.

    There are bad points in the displacement nephogram of the traditional DIC algorithm, especially in the place which these place's speckle effect is poor. Meanwhile, the displacement nephogram of higher-order DIC algorithm does not have these bad points.

  2. 2.

    In displacement nephogram of traditional DIC algorithm, contour lines are very dispersed and unsmooth. Meanwhile, in displacement nephogram of higher-order time-domain DIC algorithm, all contour lines are smooth, which conform to displacement of continuum.

Therefore, according to the above results, it can be concluded that the high-order time-domain DIC algorithm has higher accuracy than the traditional DIC method, and is more resistant to the impact of harsh environment.

Conclusion and Discussions

The high-order time-domain DIC algorithm based on the nonlinear optical flow equation is proposed. This method is suitable for solving the deformation field with continuous deformation over time. The least squares method can effectively solve this algorithm and directly measure the subpixel displacement. Simulations showed that this algorithm can calculate the displacement caused by rigid body translation, fixed axis compression, shear deformation, rigid body rotation and other deformation and can suppress the influence of grayscale changes in speckle images before and after deformation. When the image sequence contains 5 images, the error of the displacement calculation results is minimal, and the accuracy is the highest. An experimental study of resin mode I cracks shows that the accuracy of this algorithm is higher than that of the traditional algorithm. This study provides a possible theoretical foundation for analyzing high-precision deformation from low-quality images in extreme testing environments and provides corresponding technical methods.

There are still some considerations in this method.

  1. 1.

    By utilizing a smoothing algorithm, high-order time-domain DIC algorithm can help reduce noise and improve the resolution of the measured data, but it may result in the loss of information. However, with correct parameters (such as step size, number of time window, and size of relative window), it can effectively avoid this problem. This will be taken into account in future studies.

  2. 2.

    In the experiment of Sect. "Algorithm Verification Based on the Simulation", Gaussian noise with standard deviation σg = 5 is added to all deformed images, in order to prove that the proposed method can effectively reduce the noise effect. In future studies, the standard deviation of Gaussian noise will change, so as to analyze the variation of calculation accuracy with noise.

  3. 3.

    In this paper, only small displacement simulation and actual experiment are carried out, and the photosensitive resin used in the actual experiment is brittle material, so the method and simulation proposed in this paper are also established and carried out for small displacement. For large displacement, parameters such as selection of correlation window and step need to be adjusted. In theory, higher-order time-domain DIC algorithm can still adapt to the calculation of such large and complex displacement.