Measuring acceleration fields via regularized digital image correlation

Digital image correlation (DIC) is a full-field measurement technique. In instantaneous approaches (i.e., registering two images), DIC only gives access to displacement (or velocity) fields. Consequently, acceleration fields are not one of the primary measured variables. To evaluate acceleration fields, a regularization scheme has to be used. The latter may be either heuristic or mechanically motivated. The key idea of the paper is to use spatiotemporal analyses in order to explicitly measure acceleration fields. Various regularization schemes will be assessed, and their relative merits will be studied when performing uncertainty quantifications. Pyrotechnic cutting simulations will provide a set of artificial pictures to be studied in order to validate the new implementations. This analysis enables the measurement performances to be evaluated for the new implementations.


Introduction
Pyrotechnic devices are currently used in the aerospace industry to cut different parts of a launcher or activate on-board systems. However, it is hard to quantify the impact of detonations on such complex structures [1,2]. This observation calls for experimental validations of the numerical predictions. The present paper deals with the issue of acceleration measurements as a way of assessing the reliability of pay-loads against such loadings. Near the source of the shock, the actual loads are hard to quantify. A solution to circumvent the problem consists in measuring acceleration fields induced by shock waves.
Experiments on mock-ups typically consist in the separation of two bolted plates using a pyrotechnic device. Today, the kinematic measurements are performed with accelerometers and strain gauges. These point measurements have limited frequency and amplitude ranges. The purpose of the present study is to evaluate the feasibility of acceleration measurements via full-field measurements, namely, digital image correlation (DIC [3][4][5]).
With the development of digital high speed cameras [6], studies focusing on high strain rate experiments have been published. The first paper to report displacement, velocity, acceleration, strain and strain rate fields [7] was based on instantaneous stereocorrelation with a commercial code. Since then, very few publications presented acceleration fields measured via DIC [8]. Virtually all results required post-processing and filtering procedures from instantaneous DIC analyses [7,[9][10][11][12]. A first alternative consists in performing local time interpolations to measure directly, say, strain rate fields [13] and possibly acceleration fields. In global spatiotemporal DIC [14,15], the acceleration field becomes a direct output via exact differentiation of temporal shape functions [8]. Exact differentiation may also be performed after a post-processing (smoothing) of the data, but the reliability of such an data processing is never checked with respect to the original input images, in contrast to global spatiotemporal DIC. Other types of temporal regularization were proposed within the proper generalized decomposition (PGD) framework [16,17]. Although second order differentiability is not directly ensured, supplementing PGD with a suited additional regularization may provide the needed regularity.
In the following, it is proposed to extend the concept of mechanical regularization [18][19][20], which was applied to static problems, to dynamic analyses. The outline of the paper is as follows. First, spatiotemporal DIC is introduced, and its declination with the use of orthonormal modes is briefly summarized. Uncertainty quantifications are performed to illustrate the gain of temporal regularization accounting for piece-wise linear and continuous acceleration fields. Second, dynamic regularization is discussed as an extension of static regularization. Uncertainty quantifications are also performed. Last, a numerical test case illustrates the relative merits of spatiotemporal and dynamic regularizations.

Spatiotemporal digital image correlation
Instantaneous DIC is based on the registration of two frames [3]. Conversely, spatiotemporal DIC analyses the full movie of an experiment [14]. Let us consider a time interval [t 0 ; t n ], where t 0 corresponds to the first frame and t n the last frame of the movie. Similarly, the spatial domain is defined by the so-called region of interest (ROI). Spatiotemporal analyses will seek displacements minimizing the gap to gray level conservation over the whole set of frames where f denotes the image of the reference configuration, which is either the picture f t 0 (x) of the reference configuration [14] or its denoised estimate [15,16]. Equation (1) accounts for the apparent motion of any point of the ROI over time and corresponds to the optical flow [21]. It is worth noting that in the present setting is Lagrangian since the sought spatiotemporal displacement field is determined with respect to the reference configuration.
The global DIC residual, C , reads The spatiotemporal displacement field υ(x, t) is parameterized as where θ i are spatial shape functions, φ j temporal shape functions, and α ij the sought nodal (spatiotemporal) displacements when a finite element setting is chosen. Thus, high order temporal continuity can be assumed between images by using appropriate temporal shape functions φ j [8].
In the following, a Gauss-Newton scheme is selected to iteratively minimize Eq. (1). The corrections to the spatiotemporal amplitudes δα ij are computed using the spatiotemporal Hessian matrix [H st ] and residual vector {B} with and whereυ denotes the current estimate of the displacement field, and ρ the corresponding gray level residuals evaluated at any pixel position x within the ROI, and any time t between t 0 and t n at which an image was acquired. By constructing an orthonormal temporal basis [14], the previous system reduces to analyses that are formally identical to instantaneous DIC problems, and can be implemented in a non intrusive way in any standard DIC code [14,15,17]. For instance, considering the temporal mode k, the modal correction vector {δα k } becomes where [H s ] denotes the spatial Hessian matrix (just as if the temporal shape function was a Dirac, i.e., instantaneous DIC) and {β k } the residual vector associated with the temporal mode ψ k where the temporal shape function ψ k weights {b(t)}, the column vector gathering all instantaneous contributions In the following analyses, cubic B-splines were used as temporal shape functions φ j . By construction, the displacements and velocities can be made continuous. However, the acceleration field is only piece-wise continuous with such a setting. Additional constraints were added to also enforce full temporal continuity, thereby reducing the total number of kinematic unknowns (Appendix A). For the spatial discretization, 3-noded (T3) elements (with linear interpolation) were selected and regular meshes were constructed. All DIC analyses reported herein were performed within the Correli 3.0 framework developed at LMT [22] in a Matlab environment. The most computation-intensive operations of the Gauss-Newton scheme are related to the computations of the instantaneous Hessian and residual vector. The latter is updated at each iteration as the gray level residuals change. They were performed by calling binary MEX files generated from optimized C++ kernels.
It is proposed to evaluate measurement uncertainties thanks to a 500 frame film of a speckle pattern with no external load acquired with a Photron camera (SA-Z 2100 k, 20,000 fps, definition: 1024 × 1024 px, the physical size of one pixel is 200 µm). Then DIC analyses were carried out for different spatial discretizations. The latter ones were characterized by a typical length scale chosen to be the square root of the mean element surface, which is referred as element size x throughout the paper. It is worth noting that since an FE-based kinematic basis was used, the element size does not correspond to the spatial resolution of the registration scheme. With the uniform mesh used herein, the spatial resolution of inner nodes was equal to √ 6 x (i.e., each node is shared by 6 T3 elements), of edge nodes is √ 3 x (i.e., each node is shared by 3 elements), √ 2 x (i.e., each node is shared by 2 elements).
The variance was evaluated in time for each nodal displacement Then, the standard displacement uncertainty was calculated for all spatial degrees of freedom The same type of quantities was computed for nodal velocities The latter will be referred to as standard acceleration uncertainty. Let us stress that although these quantities are important, they only provide a partial view of the global uncertainty that would require the full covariance matrix. The above variances are the restriction of this matrix to its diagonal elements. The proper use of the acceleration data (e.g., for comparison to models) should make use of the entire covariance matrix to be optimal [23]. In Fig. 1, the standard displacement and acceleration uncertainties are reported as functions of the element size x . Two different sets of results are compared, namely, (i) instantaneous analyses in which no temporal regularization was used, and (ii) spatiotemporal analyses with 20 pictures per temporal element (i.e., t = 1 ms). As for the spatial resolution, t does not correspond to the temporal resolution. A first order approximation is √ 2 t . For any reported quantity, the lower the element size, the higher the measurement uncertainty. This trend corresponds to the classical compromise between The standard compromise between measurement uncertainty and spatial resolution is illustrated for both kinematic quantities. The benefit of temporal regularization is substantial for accelerations, and more limited for displacements spatial resolution and measurement uncertainty [3,4,24]. It is equally valid for displacement and acceleration uncertainties.
The standard displacement uncertainties were compared to the levels obtained from the inverse of the instantaneous Hessian [H s ] multiplied by the variance of acquisition noise (since the reference image was denoised), which is equal to the covariance matrix [14]. The square root of the mean diagonal terms (i.e., variances) is reported in Fig. 1a. A good agreement was observed between these two ways of assessing the displacement uncertainties. For spatiotemporal analyses, the spatiotemporal Hessian [H st ] was considered instead in the covariance matrix. The use of the spline functions allowed the variances of the spatiotemporal degrees of freedom to be reduced by a factor of more than 10, thereby leading to a decrease of the standard displacement uncertainties by a multiplicative factor of 0.3.
For the standard acceleration uncertainties, the gain is significantly higher and from the previous results, it leads to a gain of 0.3 5 ≈ 2.7 × 10 −3 . This prediction is close to what is observed from the a posteriori estimates shown in Fig. 1b. It is worth noting that the gain for velocity uncertainties is expected to be of the order of 0.3 3 ≈ 3 × 10 −2 . When comparing the performance in terms of temporal regularization, there is a very significant gain for the acceleration uncertainties, and a more limited one for the displacement uncertainties [8].
This uncertainty quantification allows us to conclude that only the spatiotemporal analysis provides a good acceleration measurement at higher acquisition rates. Similar trends were observed when an Euler Bernoulli kinematics was selected and uncertainties in terms of displacement and curvature were compared with standard discretizations [25]. The present regularization can be seen as an ad hoc choice of the time shape functions as called for by the needed regularity. Yet, it does not refer to the physics at play in the observed phenomena. The following section aims at adding such prior knowledge in the acceleration measurement.

Dynamic regularization
The interest of implementing mechanical regularization lies in the fact that artifactual high spatial frequency noise produced by DIC measurements can be filtered out as shown for static analyses [19,20]. There, the DIC functional was regularized with a weighted equilibrium gap functional [18], thereby dampening out displacement fluctuations that were not mechanically admissible. Such static regularization may not be compatible with dynamic experiments for which the kinetic energy is not negligible with respect to the strain energy.
The linear momentum balance equation reads where σ denotes the Cauchy stress tensor, ρ the mass density,ü the acceleration field, and f body forces. In the absence of body forces, and assuming linear elasticity, the discretized equation becomes [26] [ The quantity to be minimized in dynamically regularized DIC 2 R corresponds to the DIC functional 2 C penalized by the dynamic functional D 2 where w D is the weight put on the dynamic residual, namely, the higher the weight, the more regularized the measured displacement (and acceleration) fields. The value of this weight will be discussed later on. With the present spatiotemporal setting, the computation of dynamic residuals is straightforward. The nodal displacements {u k (t)} and accelerations {ü k (t)} for the temporal mode k read and the corresponding dynamic residuals become 2 where the penalization matrix [ k ] for mode k is expressed as The iterative minimization (8) becomes where {α k } is the current estimate of the modal amplitudes, and {δα k } their corrections. As for static regularizations, the weight w D is proportional to a regularization length D (in pixels) raised to the power 4, as can be shown by a simple power counting in Fourier space [23]. In order to set the prefactor so that this regularization length scale is given a direct meaning, a trial displacement field u * is chosen to normalize the regularization cost function. A convenient choice is a plane wave u * = u 0 sin(x/L) whose characteristic scale is defined as L (i.e., half the largest length in the ROI). The DIC and mechanical residual norms are computed from this trial field where {α * } is the FE discretization of u * . With these notations, when the regularization at scale D is sought, the weight reads The effect of dynamic (and static) regularization on uncertainty quantification was analyzed by using the same type calculations as in the previous section. For the static regularization, the contribution of the kinetic energy in the estimation of the mechanical residuals (Eq. 18) was discarded. The corresponding residual is denoted by S (instead of D ). In the following, only one (fine) spatial discretization was considered (i.e., x = 3 mm or 15 px) and one temporal element size (i.e., t = 1 ms). The regularization lengths corresponding to the dynamic and static regularizations were varied and their effect on the standard displacement and acceleration uncertainties are shown in Fig. 2. The same tradeoff between measurement uncertainty and regularization lengths is observed, namely, the higher the regularization length, the higher the penalization weight, and the lower the uncertainties. Even though the spatial discretization was small, the uncertainty level was controlled by the regularization length (once it becomes greater than the element size).
Further, the static and dynamic regularizations led to virtually the same uncertainty levels. In both cases they act as low pass filters, and the higher the regularization lengths, the lower the measurement uncertainties. When compared to spatiotemporal DIC, the order of magnitude of the measurement uncertainties was identical. This observation shows that, if discretization errors are not considered, the effect of dynamic regularization is similar to temporal discretization changes. Such a trend was also observed for static regularizations [18,20]. However, the regularization strategy allows the mesh size to be decreased at will, provided the regularization length remains sufficiently large. In the extreme case, pixel-scale discretizations may be considered as was proven for digital volume correlation [27,28].

Numerical experiment
The following proof of concept is a numerical simulation of pyrotechnic cutting. It corresponds to the dynamic separation of two assembled aluminium alloy plates with a density ρ = 2710 kg/m 3 , Poisson's ratio ν = 0.33, and Young's modulus E = 69 GPa (see Fig. 3).
The acceleration amplitude reaches 50,000g, which corresponds to the so-called zone Z1 of pyrotechnic cutting [1,29]. The boundary conditions consist in applying a load history on the plate (at rest at t = 0), which models linear pyrotechnic cutting [30]. The corresponding acceleration history is shown in Fig. 4a for the cutting signal due to ignition in the cord, and the induced shock wave in the plate (Fig. 4b). The maximum acceleration amplitude is 600,000g for the former and 300,000g for the latter.
Its history was assumed to be acquired at a 1 Mfps acquisition rate with a definition of 320 × 192 pixels (Fig. 12a). The resolution was to equal to 500 µm/pixel, which was selected to measure displacements of the order of 20 µm over an area of 160 × 96 mm 2 (Fig. 3). This size was a compromise between the low displacement amplitude of the observed shock wave and a reasonably large zone of the plate. It is worth noting that such definition and temporal sampling are compatible with existing hardware (e.g., Shimadzu HPV-X ultra-high speed camera [8,31]). To create artificial frames, the simulated displacement field was linearly interpolated at every time t of the acquisition procedure. Then, the reference speckle pattern was deformed (see Appendix B). If not otherwise stated, the gray level interpolation scheme is based on cubic splines. To be representative of an actual acquisition, the series of reference images that were deformed were corrupted by acquisition noise (whose mean standard deviation is equal to 3.1 gray levels for 8-bit pictures, see Fig. 12b). The number of frames was determined according to the camera acquisition parameters and the duration of the film (i.e., 180 µs in the present case).

Qualitative analysis
The following analyses are based on a spatial discretization of the displacement field with elements of average size x = 7.5 mm (or 15 pixels). For the temporal interpolation, B-spline elements of size t = 10 µs (or 10 frames) were considered. For the dynamic regularization, different regularization lengths were selected. Figure 5a shows the displacement history of the mesh center. The maximum amplitude was equal to 120 µm (i.e., less than 0.25 pixel). The displacement levels were reasonably well captured, as soon as some dynamic regularization was applied. Conversely, when no regularization was selected, the agreement was less good. This trend was even more pronounced in the acceleration history (Fig. 5b) for which the unregularized case led to very high discrepancies that were less important for spatiotemporal DIC, and even lower for the two regularization lengths. For the instantaneous analysis, unrealistic levels were measured even though the displacements were rather well captured.
The dynamic (mechanical) nodal residuals are shown in Fig. 6 for different regularization lengths at a given time t = 85 µs. Their levels decrease as the regularization length increases in accordance with Eq. (17).

Quantitative results
Let us now compare spatial fields and temporal histories obtained with the different DIC implementations, namely, instantaneous analyses, spatiotemporal registrations without or with dynamic regularization ( D = 75 pixels) with respect to the reference solution   Fig. 7 at time t = 85 µs (i.e., for the 86th frame). It is representative of the first propagation of the shock wave (Fig. 3). All analyses captured reasonably well the spatial distribution. The instantaneous and unregularized spatiotemporal analyses led to RMS differences of 8-9 µm (i.e., less than 0.02 pixel), which is very small given the noise level (Fig. 12b). The dynamically regularized solution lowered the RMS difference by a factor of 2 in comparison with the previous levels.
The various acceleration fields are reported in Fig. 8. For instantaneous DIC, they were computed via centered finite differences and were not filtered. The corresponding field is not representative of the shock wave simulation. It is too noisy and displays amplitudes greater than 3 × 10 6 m/s 2 , while the expected level is less than 5 × 10 5 m/s 2 . The RMS difference with respect to the reference solution is greater than the maximum amplitude. The spatiotemporal analysis with no dynamic regularization leads to lower nonphysical fluctuations (RMS difference ten times lower than in the previous case). A clear gain is observed with such approach. Even closer results The acceleration history of the central node is plotted in Fig. 9a to highlight the interest of dynamic regularization. As for the previous results, the acceleration from the instantaneous analysis is primarily corrupted by noise. If no filtering is applied, this result disqualifies this type of approach in the present case. All spatiotemporal analyses better captured the accelerations. Further, when a small regularization length is chosen, the results are close to those with no dynamic regularization. Conversely, for a larger regularization length, the results are very close to the reference solution (Fig. 9b). It is worth noting that the regularization was voluntarily turned off for the first 10 and last 10 frames and leads to  (Fig. 9c), the previous conclusion also applies.
The change of the mechanical residuals with the regularization length (Fig. 10a) proves the efficiency of the dynamic regularization in comparison with the static regularization (i.e., the contribution of the mass matrix was turned off). The two gray level interpolation schemes (Appendix B) do not yield the same results in the case of static regularization. Conversely, this effect is not observed for the dynamic regularization. Figure 10b focuses on the gray level residuals. In the present case, the dynamic regularization does not make them vary much with an increase of the regularization length D .  They are always less than 1.7% of the dynamic range for a spline interpolation of the gray levels, and 2.2% for the linear interpolation. In the present case, the gray level interpolation scheme has a more significant influence.
Last, Fig. 10c shows the global residuals R as functions of the regularization length. When the latter increases, static residuals increase while dynamic residuals remain essentially constant. This trend proves that the dynamic model describes more faithfully a shock wave, because, even if the penalization with the mechanical model increases, the total residuals do not vary much. Conversely, for the static regularization, the regulariza-tion length should not be too large as it would cut-off physical phenomena. This observation does not hold for the linear gray level interpolation scheme. These last results validate the dynamic regularization and its implementation. The effect of the gray level interpolation scheme is also observed with identical trends as those in the gray level residuals.

Conclusion
Spatiotemporal DIC analyses can be used to measure acceleration fields that are, say, continuous (in time and space when global registrations are performed). Such temporal regularization enables displacement and more importantly acceleration uncertainties to be lowered in comparison with instantaneous approaches. Such expectations were supported by uncertainty quantifications based on actual acquisitions with a high-speed camera.
Within such spatiotemporal framework, it is straightforward to add a dynamic penalization, which dampens out displacement fluctuations that are not mechanically admissible. Such approach corresponds to an enhancement of the static penalization based upon the equilibrium gap functional [32], which was introduced in DIC analyses [18][19][20]. Uncertainty quantifications showed that similar levels (as previously observed) could be achieved for displacements and accelerations with fine spatiotemporal meshes when dynamic regularization was enforced.
The benefit of dynamic regularization was also shown on an artificial test case in which deformed pictures were constructed from simulation results of the pyrotechnic cutting of a plate. Adequate choices of temporal shape functions and of the regularization model provide acceleration measurements with the expected magnitude and with mitigated influence of acquisition noise. It was shown that measurements via DIC were possible with the spatiotemporal code when dynamically regularized to quantify acceleration amplitudes of the order of 10,000g.
The next step of such analyses would be to test them under real experimental conditions to further validate the analyses preformed herein. Further, the mechanical residuals were computed within the small perturbation framework and in linear elasticity. Geometric and material nonlinearities may also be accounted for. Last, the same type of procedures could be used to analyze vibration tests.

Appendix B: Virtual deformation of images
In the present study, the reference images f t (x) were constructed by using actual pictures of the uncertainty analyses. Consequently, they included the effect of acquisition noise (and its dependence with the gray levels) so that each of them was different. The denoised reference image f was constructed by considering the mean (over time) of all pictures f t (x), see Fig. 12. Acquisition noise η(x, t) was then computed as the difference, for each pixel and each time, of f t (x) − f (x). Figure 12b shows that η(x, t) is Poissonian since its variance σ 2 η is proportional to the gray level of f . The Anscombe transform could be used to account for such effect [33,34].
With the selected resolution (i.e., 500 µm per pixel), the numerical mesh was overlaid and the displacements were interpolated for any position according to the shape functions of 3-noded (linear) elements. Similarly, the temporal changes were interpolated on the time basis of the acquisition device, which is not necessarily identical to that of the numerical simulations. In the present case, 250 pictures were generated each microsecond with a definition of 320 × 192 pixels (Fig. 12a).
The construction of the picture in the deformed configuration g t requires the gray levels to be evaluated at pixel positions ξ. Therefore, an inverse mapping is required to determine the position x in the reference configuration such that x + υ(x, t) = ξ [35]. In Fig. 12 a Denoised picture f in the reference configuration of the virtual test case. b Variance of acquisition noise σ 2 η as a function of the gray level of the denoised picture. f is constructed by using an actual picture series of the uncertainty analyses. Its definition is 320 × 192 pixels. The acquisition noise is Poissonian since its variance is proportional (solid red line) to the considered gray level the present case, a Newton's method was used to determine x. Once the position x(ξ) was known, a gray level interpolation scheme enables f t (x(ξ)) to be computed, and invoking gray level conservation, it is equal to g t (ξ). This operation was repeated for all pixels and all frames considered in the analysis. In the studies presented herein two gray level interpolation schemes were considered, namely, bilinear and cubic splines as implemented in Matlab [36].