Michel Bornert

François Hild

Jean-José Orteu

Stéphane Roux

Digital Image Correlation

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Digital Image Correlation

Background

Digital Image Correlation or DIC, which appeared in the early 1980s [LUC 81, BUR 82, SUT 83, SUT 86] has a major impact in the field of the Mechanics of Solids and Structures. Nowadays, it is still undergoing very spectacular developments.

The challenge is to measure displacement fields of surfaces (or in volumes) of stressed specimens and structures from images acquired at different stages of loading. A specific advantage of this tool is that it exploits numerical images, which are usually acquired by optical means. Imaging devices have made significant progress not only in terms of quality and definition, but also (lower) cost. These imaging means are inherently contactless, non intrusive, tolerant to aggressive conditions (e.g., temperature, chemical environment), easy to use, efficient and cheap, many of these features that can only be appealing in the context of mechanical tests. DIC can easily be used at different scales of space and time to the extent that it relies on principles applicable to pictures obtained by very different imaging systems. It is nowadays possible to use images shot by fast and ultra-fast cameras at time scales down to the microsecond or less [SCH 03a, SCH 03b, SIE 07, TIW 07, BES 08, BES 10], acquired by a scanning electron microscope (SEM) [DOU 00a, DOU 00b, SOP 01, DOU 03, TAT 05, SUT 06, SUT 07a] or an atomic force microscope (AFM) [CHA 02, CHO 05b, CHO 07a, CHO 07b, HAN 10] at nanometric scales, but also satellite images at geophysical scales [SCA 92, LEP 07]. Multicamera systems give access to 3D shapes and displacement fields of surfaces of an observed object (see Subsection 6.4). Three-dimensional images obtained by computed (micro)tomography [BAY 99, BOR 04, LEN 07, BAY 08, RAN 10] or magnetic resonance imaging [NEU 08, BEN 09] can also be utilized to measure 3D displacement fields in the bulk of various (optically opaque) materials.

The wealth of kinematic data obtained allows not only for a quantitative exploitation through identification techniques (see following chapters), but also to validate this identification, to enrich or degrade it according to the needs. By specifying appropriate forms of the displacement fields to be measured, image correlation can directly address this identification step (see Subsection 6.2.8).

Surface and Volume Digital Image Correlation

DIC techniques can be applied indifferently to classical two-dimensional or to volumetric images to which recent imaging tools give access. As a consequence, the term "pixel" will refer to in the following as the elementary discrete datum of a digital image defined in a two-or in a three-dimensional space. We detail in this section the algorithms providing evaluations of the apparent mechanical transformation Φ a that links two images of the same mechanical system under two different configurations.

After a first subsection devoted to guiding principles and ingredients common to all presented algorithms, the discussion then will focus on so-called "local" image correlation techniques, which evaluate the transformation Φ a piecewise, through a large number of independent analyses on sub-images, called "correlation windows" or "domains." Various flavors of the more recent "global" approaches, which may be enriched by some a priori mechanical information suited to each problem dealt with, will subsequently be described. The main sources of uncertainties and their quantification will be discussed in Section 6.3.

Images

The input data of the analysis are positive integers, called gray levels, of the image of the first configuration, known as the "reference image", noted f I , and those of the second one, called the "deformed image", noted g I , where the subscript I refers to a pair of integers (column, line) in the case of 2D images or to a triplet (column, line, plane) if 3D images are considered. These integers vary between 0 and some maximum value related to the digitization or encoding depth (or dynamic range) of the image sensor (e.g., 256 levels for 8-bit images, 4096 for 12-bit images). Indices I vary in planar or volumetric domains, whose extension characterizes the image definition (e.g., [0; 1023] × [0; 1023] for a one-megapixel 2D image), and in which it is possible to define continuous positions, with real-valued coordinates x 1 . The two images can be extracted from a temporal sequence comprising a large number of images, of which any can be selected as the reference image or the deformed image. However we do not address in the following the algorithms that process such a sequence globally.

The images f I and g I result from a complex acquisition chain, today evermore various. It is essential to consider this acquisition process in the final interpretation of the estimated transformation, and to wonder about the presence of biases that could result from it. In particular, the analysis of this chain would make it possible to specify the link between some physical quantity continuously varying in space, denoted by f , and the discrete gray level f I that this quantity induces through the imaging system at position I in the image. Because of space limitation, this aspect will not be developed herein.

DIC principles rest on the essential assumption that the physical quantity that leads to the image is associated with some physical property of the matter that constitutes the analyzed system, and, as a consequence, is transported by the sought apparent mechanical transformation Φ a . Noting g the continuous physical quantity associated with the deformed image, this assumption writes

g(Φ a (x)) = f (x).
(6.1)

Texture of images

As will be made more specific later on, DIC exploits a texture that must be the signature of each surface (or volume) element, simply transported, according to Equation (6.1), by the displacement field, without any other deterioration. Moreover, the gray levels representing this texture in the images need to exhibit a broad dynamic range that covers as much as possible the available encoding depth of the images, i.e., from 8 to 16 bits, without however showing saturation. Lastly, it is desirable to have strong contrasts from one pixel to the next to be able to be sensitive to small displacement amplitudes.

The simplicity of the implementation of the method would recommend to rely on the natural texture of the studied material. This is sometimes possible, as illustrated in Figure 6.1. However, the correlation function of the texture will then dictate the performances of the analysis. Thus, for the example of the steel shown in the Figure 6.1a, the various phases revealed by a selective chemical etching have an orientation that will allow a good estimate of the displacement along the horizontal direction but a worse one in the perpendicular direction. In this context, the surface topography of silica observed by AFM (Figure 6.1c) does not lead to a strong sensitivity. Conversely, the phase image of the same zone (Figure 6.1d) can lead to a much better resolution. The selection of an appropriate mode of observation, associated, when required, with an adequate surface preparation, is thus an essential element for the use of natural contrast in DIC. The existence of such a contrast is not guaranteed. A possible way to circumvent this difficulty is to deposit an artificial texture on the surface of the studied samples as illustrated in Figure 6.2. The pulverization of fine droplets of black paint on a white background (or conversely) is a process called "speckle painting" that can be used on many materials. The airbrush allows to form an aerosol whose droplets size can be adjusted by an adapted nozzle. For the largest scales (for example for civil engineering structures) marking can be done with a stencil key set. A limitation can emerge from the strong deformability of materials as in the example of Figure 6.2c which is an elastomer. It is difficult for a painting to follow deformations of large amplitude. In this case, a powder deposit (talc in this example) allows for a good stability of texture up to several hundreds percent of strain.

The deposition of particles on the surface can also be used in the context of scanning electron microscope (SEM) imaging [SCR 07]. Electron lithography techniques, in particular in their simplified version that can be implemented without expensive equipment other that the SEM itself [ALL 94], make it possible to mark surfaces by means of microgrids with a step varying from a few tens of micrometers to less than one micrometer, by metal deposition or chemical engraving. The periodicity of the obtained local contrast requires however some minor adaptations of the correlation algorithms to avoid offsets of one or several grip steps in the evaluation of the displacement field. An advantage if this type of marking is that it does not hide the microstructure and thus facilitates the micromechanical analyses of heterogeneous materials [DOU 00a, HER 07]. However the analysis cannot be as dense as with a fine speckle painting. Lastly, when three-dimensional images of computed tomography are to be exploited, a texture in the bulk of the studied material needs to be relied on (Figure 6.3). Some materials are then privileged, like most rocks [LEN 07], stonewool [HIL 09b], nodular graphite cast iron [RAN 10], or foams [ROU 08], whose micro-pores, inclusions or heterogeneous microstructure behave as a speckle painting. In the absence of such contrast or within the framework of an analysis at a finer scale, it is sometimes possible to reinforce it artificially (for example by diffusion of a heavy element to mark grain boundaries, or by addition of fine markers during the processing of the material [BOR 04]) but the difficulty is then to make sure that this marking does not introduce any modification of the mechanical properties, or to take into account the modified microstructure of the material in the interpretation of the results.

Guiding principles

The aim of the exposed methods is to evaluate the apparent transformation Φ a on a Region Of Interest (ROI) R of the reference image starting from the knowledge of the gray levels f I and g I . This essentially is an ill-posed problem insofar as the available information (i.e., the gray levels of the pixels) is insufficient to uniquely determine a vectorial displacement at each pixel. It is thus imperative to regularize the problem by restricting this determination to a particular family of transformations Φ 0 (α, •), parameterized by N scalars, α, written here in vector form. Image correlation methods can then be given the following general formulation

α min = Argmin α∈V α C(Φ 0 (α, •),R,[f ], [g]) (6.2)
where the correlation coefficient C is a scalar that measures the similarity2 between the ROI R of the image [f ] (i.e., the set of the pixels f I ) transformed by Φ 0 (α, •) and the corresponding region of the image [g]. The parameters α belong to a certain domain V α of R N . The numerical resolution of this minimization problem leads to an evaluation Φ 0 (α min , •) of the sought transformation Φ a (•).

The various available image correlation algorithms derive from specific choices of the expressions of the correlation coefficient C (section 6.2.4), of the considered family of transformations and its parameter setting, which distinguish in particular the local methods (Section 6.2.6) from the global ones (Section 6.2.8), as well as of the minimization algorithms in use. Moreover, the evaluation of a continuous transformation from the knowledge of discrete values of gray levels provided by the input images, requires to recourse to some interpolation method at some stage of the minimization problem. The various ways to do so define additional alternatives to specify image correlation algorithms (Section 6.2.5). These interpolation algorithms make it possible to evaluate the local displacement field with a resolution notably lower than one pixel, known as subpixel accuracy. For the sake of simplicity, the discrete nature of the input data will be ignored in the following presentation, except in Section 6.2.5 where this question is explicitly addressed. This simplification consists in assuming the pixels infinitely small with respect to the spatial gray level fluctuations, so that discrete gray levels f I and g I can be represented by continuous functions f (x) and g(x).

Correlation coefficients

The simplest and most classical similarity criterion is the sum of squared differences between the images3

C 1 = R [f (x) -g(Φ 0 (x))] 2 dx = R [f (x) -g(x + u 0 (x))] 2 dx (6.3)
This quantity is clearly positive. It vanishes when the transformation Φ 0 coincides with the apparent transformation Φ a and when the conservation condition (6.1) is exactly satisfied not only for the physical quantities f and g at the origin of the gray levels in the images, but also for the continuous gray levels f and g. The choice of a squared difference induces some operational advantages but may be replaced by an absolute value or any other measure of discrepancy.

Some imaging systems (e.g., SEM) do not make it possible to guarantee the stability of the conversion into gray levels f I of the physical quantity f at the origin of the image contrast and convected by motion according to Equation (6.1). One can then choose a less demanding similarity measure such as

C 2 = min a,b R [f (x) -(ag(Φ 0 (x)) + b)]
2 dx (6.4) that consists in seeking the best linear regression between the gray levels of the two paired images. The optimized coefficients a and b can be interpreted as changes in contrast and brightness, respectively, of the image between the two configurations 4 . The optimization in Equation (6.4) is straightforward and leads to a = R (f -f )(gḡ)dx/ R (g -ḡ) 2 dx and b = f -aḡ, where f and ḡ are the averages of f and g • Φ 0 over R.A s R (f -f ) 2 dx does not depend on Φ 0 , optimizing C 2 is then equivalent to minimizing

C 3 =1- R (f -f ).(g -ḡ)dx R (g -ḡ) 2 dx R (f -f ) 2 dx (6.5)
which can be interpreted as the correlation coefficient of a linear regression and varies from 0 (perfect similarity) to 1 (no link between images) 5 .

The definition of these criteria leads immediately to the definition of the image of residuals, whose gray levels defined on R are given by f (x) -g(Φ 0 (x)) or f (x) -(ag(Φ 0 (x)) + b), depending on which correlation criterion is considered. The qualitative or quantitative analysis of these residuals is a means of evaluating the relevance of the optimal transformation with respect to the actual unknown transformation. In case the registration would be perfect, and the residuals would contain nothing but image noise. It is worth noting that in this case and under the assumption of a moderate noise level, the criterion C 3 is equal at its optimum to (1/2)(σ(f ′ )/σ(f )) 2 where σ(f ) et σ(f ′ ) are the standard deviations of the noiseless image f and of its noise f ′ on R. With an optimal contrast occupying the full dynamic range of the gray levels and a noise level of about 1% of the latter, as often observed with images of conventional optical cameras, the minimum of the C 3 coefficient might become as small as 10 -3 . This concept of residuals will again be discussed in Section 6.2.8.

The above list of criteria is not exhaustive [CHA 11], but gathers those most commonly used for quantitative analyses in mechanics. Other similarity measures might be used within other image registration applications (e.g., satellite or medical imaging), in particular when the condition of convection of the gray levels is not well satisfied, because for instance of high noise levels [BOR 04] or of a too strong change of image contrast between the various configurations.

Sub-pixel interpolation

All DIC algorithms need to specify the interpolation method required to access a resolution in displacement notably less than one pixel. One should be aware of the fact that the associated technical choice might have strong consequences on the measurement accuracy.

The computation of the correlation coefficient, whatever the formulation kept among those presented in Section 6.2.4, needs to face the discrete nature of the images. The integration in equations (6.3)-(6.5) is in practice replaced by a summation over a 5. When suppressing the absolute value of the numerator of this expression, the criterion varies between 0 and 2, the maximum value corresponding to a "perfect inversion" of the contrast. set of pixels in the reference image and requires to specify how the gray level in the deformed image at point Φ 0 (x(I)) is calculated. This point has in general non integer coordinates and its gray level can only be evaluated by means of some interpolation scheme making use of the gray levels of surrounding pixels. There are numerous possible choices to do so, namely, evaluation by the gray level of the nearest neighbor, bilinear interpolation making use of the 4 first neighbors (8 neighbors for a trilinear interpolation in 3D), bi-cubic interpolation (various possible formulations) based on the 16 first neighbors (64 in 3D), spline functions to various orders, Fourier transforms on windows of adjustable size, or even wavelets. A high order interpolation will ensure the continuity of the correlation coefficient and of its derivative to various orders with respect to the kinematic parameters α, required by some optimization techniques, but will be greedier in computation time.

The errors induced by this interpolation takes the form of an over-or underestimation of the displacement components, which depends, on average, on its fractional part, expressed in pixels, as a consequence of the periodicity of the sensor behavior. The average of these errors is thus a function of this fractional part, which exhibits a symmetry with respect to the half-pixel displacement, and is often referred to as the "S-shaped" systematic error curve [CH0 97, SCH 00]. The precise shape and amplitude of this curve express the more or less good efficiency of the used interpolation scheme to restore the evolution of gray levels in an image induced by a real sub-pixel translation. This error needs to be quantified and, if possible controlled when small strains are to be analyzed, by means of the selection of an interpolation scheme adapted to the actual texture of the image, or conversely, by means of a modification of the latter [YAN 10], which depends on the marking of the sample but also on the parameters of the optical system.

Various procedures have been proposed to determine this S-shaped curve. One can for instance generate an artificial translated image from some real reference image and compare the results of the DIC analysis to the prescribed translation. However the representativity of the virtually translated image with respect to a real one itself depends of the used numerical translation procedure. A Fourier-based translation [SCH 00, HIL 06], consisting in a multiplication in Fourier space by a phase term, e ikδ , for each wave number k, is a particular procedure that leads to a C ∞ interpolation of the reference image, which exhibits some satisfactory properties with respect to this question. Another possibility is to numerically simulate the integration of the optical signal by the sensor [BOR 09, DUP 10] with now the difficulty to define a good numerical model of the sample texture. A last option consists in analyzing a pair or a sequence of real images corresponding to a well known transformation. This can be done but is delicate in practice for a pure translation [WAN 09], but can easily be achieved with out-of-plane motions that induce homogeneous apparent deformation gradients [YAN 10, DAU 11]. Recent analysis have also shown that this systematic error is strongly sensitive to image noise [WAN 09], at least for some algorithms [DUP 10].

It would ideally be advisable to reduce the amplitude of the systematic errors to a level less than that of the random errors, which correspond to the fluctuations of the measured displacements. The amplitude of the latter may also depend on the fractional part of the displacement field. When systematic errors are appropriately monitored, random errors constitute the actual limiting factor. This question will be addressed again in Section 6.3.2 to discuss measurement uncertainty.

Local approaches

In traditional approaches the transformation Φ 0 over R is decomposed into a multitude of independent and local transformations,orshape-functions, parameterized by the coefficients of their local expansion near centers x 0 and used in the neighborhood D x 0 of these centers

∀x ∈ D x 0 , Φ 0 (α, x) -x = u 0 (x)= i α 0 i + j α 1 ij (x j -x 0 j ) •••+ jk α 2 ijk (x j -x 0 j )(x k -x 0 k )+.... e i (6.6)
where e i is the unit vector along direction i in the image and x i the coordinates of point x (i ∈{ 1, 2} in 2D and i ∈{ 1, 2, 3} in 3D). Zeroth-order methods are limited to the translation vector α 0 and were historically the first ones. First-order methods are the most common ones and take into account the components of the first deformation gradient α i j 1 . There are thus six parameters to locally describe Φ 0 in 2D and twelve in 3D. Some variants include in addition the second order cross-terms (α 2 112 and α 2 212 in 2D, leading to a total of eight coefficients). One may also choose to restrict the local transformation to rigid motions (i.e., two translation components and one rotation angle in 2D, three components and three angles in 3D). Higher order expansions are seldom used.

With such a choice of local parametric description, the standard processing procedure consists in defining in the region of interest R a set of points x k , often organized into rows and columns (and planes in 3D) and usually regularly spaced, and to associate to each of them a small local domain D x k , in general square (or cubic in 3D), called correlation window, subset or domain. Problem (6.2) is then decomposed into as many independent optimizations leading to local parameters α k for each point x k , which define local estimations of the transformation Φ a over the domain D x k , from which in general only the value Φ 0 (α k , x k ) at the center is retained. The final result is then a discrete collection of displacement vectors on a set of points x k ∈ R.

Such a description of the mechanical transformation by means of a set of local shape-functions requires thus essentially three choices, namely, the order of the local transformation, the size of the correlation window and the sampling step of the ROI. These choices are usually left to the user in classical image correlation softwares, but they have some implications on the relevance of the kinematic measurements (displacement or strain). While questions relative to experimental errors will be addressed in detail in Section 6.3, the following aspects briefly are emphasized:

-the selection of the order of the local transformation is strongly linked with the window size. The systematic analysis performed in reference [BOR 09] on virtual images exhibiting non uniform deformation gradients with variable characteristic lengths demonstrates the existence of two error regimes. The first one is linked to the inaptitude of the shape function to describe the real transformation over the whole correlation window. This error is all the more important as the window is large, and that the order of the transformation is low. The second error regime, referred to as the ultimate error, coincides with the one encountered when the real transformation is a pure translation. It is on the contrary characterized by a decrease with window size and an increase with the order of the shape-function. This latter phenomenon is discussed in detail in Section 6.3.2; -sampling step and window size can be chosen independently. One often selects a step less than or equal to the window size, so as to make use of the whole information contained in the image. An arbitrarily fine sampling of the image can even be chosen, at the price of an increased computation time. However one should be aware of the difference between this measurement step and the spatial resolution of the measurement of the displacement field, which remains controlled by the correlation window size. The displacement measurement associated with each position of the ROI sampling is the result of some averaging process over the whole correlation window;

-the sampling step may be a parameter entering implicitly or explicitly in the procedures used to compute strain components from the discrete values of the displacement components evaluated at each measurement point. We refer to Chapter 7 for a deeper discussion of this question. Note however that the accuracy of the measurement of the components of the deformation gradient relative to some gauge length is inversely proportional to the latter. When low strain levels are to be quantified, one should then privilege larger sampling steps or make use of differentiation procedures based on displacements evaluated over several sampling steps [ALL 94]. Moreover, in case of overlapping correlation windows, measurements of displacement at neighboring positions are no longer independent, which may induce some bias in the computation of the gradient components.

-the DIC analysis results in a discrete sampling of the displacement vectors (and their gradients) over a regular mesh of the ROI, the value at each sampling point being itself some average of the local displacement field over the correlation window. In order to define a displacement field at each position x in the ROI, these values need to be somehow interpolated. There is however no systematic way to evaluate the distance between the interpolated displacement field and the experimentally investigated field.

To conclude this short description of local formulations, let us note that they can also deal with an arbitrary selection of measurement points, adapted to the problem under consideration. These points may for instance be defined as the nodes of the finite element mesh of a structural calculation modeling the experiment [HER 07], as some isolated points of interest when the analysis does not require a dense evaluation of the kinematics -in such a case, the DIC techniques would operate as a markers tracking technique -or as the only points with an appropriate local image contrast, in case of a DIC analysis of a system with a non uniformly distributed natural local contrast. Local formulations can also make use of arbitrarily shaped correlation windows, for instance associated with some discrete entities of the considered system (such as the grains in a sand sample [HAL 10]), or delimited by characteristic boundaries of the system (such as for instance the interface between constitutive phases of a heterogeneous material [RUP 07]). The use of image masks is useful for the practical implementation of such alternatives of the general algorithm.

Optimization algorithms

At this stage, it is possible to compute C and its gradients for any set of parameters α. The optimisation, Equation (6.2), can then in principle be performed with any iterative numerical algorithm. A minimization algorithm making use of gradients (e.g., gradient descent, Newton-Raphson, Levenberg-Marquardt) or not (e.g., simplex, Powell) will allow to converge more or less quickly towards a local minimum close to the initializing set of parameters. In order to avoid convergence to an erroneous local minimum, various strategies might be employed, such as the preliminary systematic scanning of a set of displacement fields sufficiently large to contain the actual one, the use of known displacement fields relative to a previous stage of the mechanical test, or of the displacement of already analyzed neighboring positions, in the context of a "propagating front" type algorithm, or even, in the most difficult cases, the assistance of an operator. In this context, the preliminary systematic scanning is often restricted to integer components of the translation vector α 0 i , with other possible parameters set to zero or frozen to some other value determined by some other mean. To speed up this calculation stage, the gray level interpolation scheme may be set to a simple nearest neighbor algorithm; an alternative is to compute correlation coefficients in Fourier space, by means of efficient fast Fourier transform algorithms.

It is also possible to make use of these discrete evaluations of the correlation coefficient to interpolate this quantity by a bi-quadratic (or tri-quadratic in 3D) polynomial function near its discrete optimum, making use of the 9 (or 27) discrete values surrounding this optimum, and to analytically seek the continuous optimum of this interpolating function. Such an algorithm allows to reach a sub-pixel accuracy at a particularly low computation cost [WAT 01], but is restricted to the optimization of a limited number of components of the vector α.

Another proposed strategy to avoid local minima consists in a multi-scale algorithm [GAR 01a, HIL 02]. The idea is to start with strongly low-pass filtered images in order to artificially enlarge the well associated with the global minimum and to avoid local trappings far away from the solution. The obtained estimate is however not accurate as details of the image have been filtered out. The strategy consists then in gradually restoring these details and to perform the optimization again with these richer images, starting the algorithm with the displacement obtained on the coarser images, until the original image is recovered. This procedure (which also applies to the global approaches described in following section) turns out to be extremely robust and makes it possible to benefit from an optimization by descent, while limiting the effects of local minima. This wide variety of options, which can be combined with each other, leads to a very vast range of DIC algorithms, which in addition will be enriched by the global approaches described in the following section. The selection of a particular combination of options is a compromise, whose definition depends on the available local image contrast, the awaited level of strain, the required measurement resolution and its spatial resolution, as well as on the numerical cost and the robustness of the algorithms.

Global Approaches

The main difference between local and global approaches is the choice of the displacement basis used to account for the transformation Φ 0 and its spatial definition. The local approach results in an independent calculation of the transformation on each correlation (small) window. However, no specific regularity of the displacement field is exploited. Conversely, it is possible to make the choice of a continuous transformation basis in a global approach. The displacement field is expressed as

u 0 (x)=Φ 0 (α, x) -x = i α i ψ i (x) (6.7)
where the fields ψ i form the chosen kinematic basis that can be defined over the whole region of interest. There is no restriction at this stage for the definition of this basis. The superposition of different displacement fields throughout domain R makes interdependent the sought amplitudes α, hence the term global to address this problem. The determination of α is performed as above, namely, via the minimization of the sum of squared differences, C 1 (Φ 0 (α, •),R,[f ], [g]). A Newton-Raphson scheme, which resorts to successive linearizations, is implemented to determine the solution, provided the amplitudes of the displacement fields are sufficiently small. As mentioned above, a multi-scale strategy can significantly relax this requirement [BES 06], and large displacements can be measured.

Unlike interpolated fields obtained with a local approach, each pixel is directly confronted with the value of the displacement field at this point and no interpolating step is required. The correlation residual normalized by the dynamic range of gray levels is comparable to that encountered in local approaches, i.e., of the order of one percent. As mentioned in Section 6.2.4, the residual field allows the correlation errors to be spatially located. This property allows the user to quickly check for the good convergence of the algorithm, and especially to possibly correct / extend the chosen basis of fields. For example, if a fractured medium is observed and the chosen displacement basis is continuous, then the residual error will focus on the support of the kinematic discontinuity. This is a very effective way of precisely locating the crack front in all its geometric complexity [RAN 10] as shown in Figure 6.4.

The nature of the kinematic basis was not specified beyond its overall character. Great latitude exists for this choice. Sometimes some specific information (or assumptions) on the expected displacement field are available, and this is a valuable guide in the choice of this basis. Let us cite some examples:

-Without any specific information on the kinematics, the representation of the displacement field by shape functions of finite elements is a convenient possibility that will enable for an easy and direct interface with numerical simulations. So-called Q4elements in two dimensions [BES 06], or C8-elements in three dimensions [ROU 08] (i.e., four-noded quadrilaterals with bilinear interpolation (2D) or 8-noded cubes (3D) with trilinear interpolations) are very well adapted to the discretization of the image. An example of a displacement field obtained with Q4 elements is shown in Figure 6.5(left). It is also possible to work with an unstructured mesh [LEC 09] or to use different shape functions [RET 09a]. Specific finite elements can then be considered with an enriched kinematics as in the extended finite element method (X-FEM) framework [RET 08] to add discontinuities in a given mesh (see Chapter 14). Last, the fact that discretizations of the displacement field are shared with numerical modeling makes the coupling direct and seamless between image correlation on the one hand with e.g., an elastic calculation on the other hand (see [RET 09b] for such an example).

-The test can be analyzed with an analytical expression for the displacement field in the context of an elastic behavior. This solution can be an element of the selected kinematic basis. It is generally supplemented by additional degrees of freedom (such as rigid body motions). An example of this approach is the analysis of diametral compression (i.e., 'Brazilian test' [HIL 06]) or a point force loading (i.e., Flamant's problem [ROU 05]). Note however that the solution involves the Poisson's ratio ν of the material. If it is not known, the general solution can be split into two terms, one dependent on ν, and the other not. Similarly, the amplitude of the displacement field (except rigid body motions) is proportional to the ratio of the mechanical loading level to the shear modulus of the material. Thus by measuring the amplitudes of the basic fields, associated with the knowledge of the load level, allow to directly assess the elastic properties of considered medium [ROU 05]. Beyond the pure kinematic measurement, we begin to see the identification of mechanical properties (elastic in that case) without additional work via global digital image correlation.

-A particularly important case concerns cracks [ROU 06, ROU 09]. A family of analytical solutions exists, namely, Williams' series [WIL 57], which satisfies the cancelation of the tractions on the crack mouth. An appropriate selection of some terms of this infinite series can account for the sought displacement field in the vicinity of the crack tip. Among the selected fields, it is natural to include those associated with a 1/ √ r singularity for the stresses and strains. The amplitudes of these fields are related (via the physical size of the pixel and elastic properties of the medium) to the stress intensity factors in modes I and II. An example of such a determination is shown in Figure 6.5(right) for the same pair of images as that used in the left part of the figure. However, this family of fields requires a priori knowledge of the geometry of the crack (here assumed to be straight). The support of the crack is generally easy to appreciate as apparent in the displacement field or in the residual field (see Figure 6.4). The position of the crack tip is yet another issue. One solution is to use specific terms of Williams' series (i.e., strain and stress fields with r -3/2 singularity [ROU 09]). The amplitudes of these fields provide information about the actual position of the crack tip.

-When considering slender objects, beam-type models are often used. This description assumes a particular (e.g., Euler-Bernoulli) kinematics that can be taken into account as in the displacement basis. This approach then allows to assess by image correlation the motion of a beam or a frame in a language directly suitable for modeling purposes [HIL 09a]. In this context, it is possible to finely characterize nonlinear modes (e.g., plastic hinge formation) in the framework of the strength of materials (i.e., constitutive law of a localized mode, equivalent position of a point corresponding to a localized failure mechanism, see Chapter 12). Moreover, the measurement with an elastic kinematics over a given part of the beam gives access to the actual loading conditions.

-Beyond elasticity, it is sometimes possible to have a priori information to define a suitable displacement basis. The formation of shear bands in a uniaxial tensile test results in specific fields obeying Hadamard's conditions that can be accounted for (see Chapter 14).

-Last, the basis fields can be computed numerically. If one wishes to consider several unknowns (e.g., boundary conditions, material parameters) it suffices to calculate the various fields associated with changes in the parameters associated with the unknowns of the problem and consider them as a kinematics basis. The value of the amplitudes of these fields correspond to the best evaluation reflecting the observed situation [LEC 09].

Errors et uncertainties

Image correlation gives access to quantitative measurements of displacement fields. It is therefore crucial to assess the sources and levels of error and uncertainty associated with its implementation.

Main Error Sources

The first source of error is extrinsic to DIC because it results from image acquisition. Among them, let us cite:

-those related to the 3D-2D projection, which will be described in more details in Section 6.4. In particular the out-of-plane motions can induce artificial yet quantifiable displacements and strains [SUT 08b]. The use of telecentric lenses largely restricts some of these artifacts. It is also possible to correct for variations in optical magnification induced by global motions of an object with respect to the imaging system if one has an additional means for quantifying them [YAN 10].

-those related to the imperfect positioning of pixel coordinates. For optical images, the distortions induced by objective lenses fall within this category. For images obtained by scanning (e.g., SEM, AFM) the imperfect positioning of the locally probed area (by the electronic spot or cantilever tip) may also introduce significant spurious strains.

-those due to intrinsic noise in images. Two images acquired in identical conditions have gray levels that differ by a random amount whose standard deviation varies typically between a few thousandths (very good camera) and a few tenths (SEM image acquired very quickly) of the dynamic range of the gray levels. Its impact on image correlation [ROU 06] is detailed in Section 6.3.3.

A second category of errors is due to DIC:

-The selected basis of displacement fields is unable to describe the apparent transformation Φ a . This error regime has already been mentioned in Section 6.2.5 and has been investigated in detail in Reference [BOR 09] for various softwares based on local approaches.

-The second type of error is related to the way the gray levels (or correlation coefficient) are interpolated in the correlation procedure when moving images for noninteger values (Section 6.2.5).

One could also mention the errors induced by possible changes of local contrast and brightness invalidating the basic assumption of gray level conservation (6.1). If their modeling is accessible, it is possible to take them into account [HIL 12b] in the correlation procedure (e.g., for the analysis of AFM images [HAN 10]). Last, when talking about errors or uncertainties, it is important to remember that the correlation residuals are very good indicators. They provide information on the local quality of the registration between images. They also allow for the modification or enrichment of the chosen kinematic basis [RET 08, RAN 10].

Uncertainty and Spatial Resolution

The choice of the discretization of the displacement field is crucial. Systematic studies on virtual images show two different regimes of error occur [BOR 09]. The first is related to the inability of the shape function to describe the actual transformation of the entire correlation window (for a local approach) or the discretized kinematics (for a finite-element based global approach). The second type of error, referred to as ultimate, is similar to that encountered for pure translation transformations. It is characterized by a decrease of the error with the window size and increased order of the function shapes. It is a direct consequence of the ill-posedness of the correlation problem (6.2). The information contained in a small correlation window will only identify a limited number of parameters and the increase in their number will be accompanied by an increasing uncertainty of the value of each of them. The user is faced with a dilemma, which is comparable in some ways to Heisenberg's uncertainty principle, where position and displacement cannot be fully resolved simultaneously. Good spatial resolution (i.e., small window or element size) will be accompanied by a high displacement uncertainty and vice versa. A way of breaking this limit is to resort to regularized DIC [ROU 12] or Digital Volume Correlation [LEC 11].

To quantify a priori the measurement uncertainty, we restrict ourselves to the second case by choosing a simple kinematics (i.e., uniform translation) that belongs to the space of shape functions. We proceed as in the evaluation of the interpolation sub-pixel error (Section 6.2.5). From the displacement field evaluated on an image pair made of the reference picture and its translated copy, the standard deviation of the displacement field gives an a priori estimate of the measurement uncertainty. This uncertainty is usually much higher than the systematic error. It is customary to observe a dependence of the standard displacement uncertainty, σ u , with the size, ℓ in pixels of the interrogation window (local approach) or elements (global approach using finite elements) as a power law

σ u = Aℓ -η (6.8)
where A is of the order of unity, and η varies between 1 and 2 as appropriate. This allows measurement uncertainties to be evaluated (excluding other factors related to noise or complexity of the observed displacement field) that commonly reach the hundredth or the even thousandth of one pixel in the most favorable cases.

From this estimate it is also possible to evaluate the strain uncertainty when they are obtained, say, by simple finite differences from the displacement field. It should be noted however that the fluctuations will be observed on strains are strongly correlated (or anti-correlated) at small scales. This is especially important for the subsequent use of those kinematic measurements for the identification of mechanical properties (see second part of this book) when these correlations are not taken into account. It is therefore advised to prefer formulations of inverse problems that are based on displacements rather than strains, or to be aware of the implicit filtering often linked to the strain calculation (see Chapter 7).

Noise Sensitivity

Another limiting factor is due to the presence of noise in images. This noise can be of different nature and therefore has very different statistical features. It can sometimes be reduced by image averaging, provided that the sample does not evolve (i.e., move) between acquisitions. A good approximation is often given by white noise (i.e., uncorrelated from pixel to pixel). Its effect can be followed up to the estimate of the displacement fluctuations. In the case of a local approach, when the correlation windows are not overlapping, then the displacement fluctuations exhibit no correlation between different discrete estimates. However when these windows overlap, or in the case of a global approach, the correlation matrix of the kinematic degrees of freedom is no longer diagonal but its calculation is possible for a low noise level using a linearized operator for the final stage of the displacement estimation [HIL 12a].

In general, the standard displacement uncertainty (i.e., its standard deviation σ u ) is proportional to that of the noise images ǫ (by linearity). The amplitude will be inversely proportional to the mean gradient of the image, and will decrease as a power law with the number of pixels that are used to measure the considered motion. Thus, in d dimensions, for both local and global approaches

σ u ∝ ǫ ℓ d/2 |∇f | 2 1/2
(6.9)

A more detailed discussion of this evaluation is provided in References [BES 06, ROU 06, HIL 12a].

It is possible to take this expression as a practical recommendation on image quality and correlation parameters to use. Regarding the acquisition, it should minimize image noise (i.e., ǫ), and to acquire images with a large depth of gray level (as long as it is not effectively truncated by the noise level ǫ). The texture of the analyzed images also affects the displacement uncertainty. The dependence on the gray level gradient calls for textures having a high contrast, and secondly having a short correlation length. The limit to this suggestion is the ability to interpolate gray levels at a sub-pixel level, which instead requires a regular variation at the pixel level. A good compromise is to deal with textures whose correlation length (radius) is of the order of 2 to 5 pixels. Last, for the parameters of DIC itself, the size of correlation windows should increase to reduce its sensitivity to noise, but at the cost of lower spatial resolution. Again, a compromise must be made. One can also note that the sensitivity of different correlation algorithms vis-à-vis this image noise is not identical [DUP 10].

Stereo-Correlation or 3D-DIC

The main advantage of 2D DIC techniques described in Section 6.2 is their (apparent) simplicity and versatility, namely, a single camera is sufficient and the object surface preparation (when required) can be easily done (spraying paint is an easy way to create an appropriate speckle pattern at the surface of the object). However, some important points need to be kept in mind:

-by using a single camera only in-plane displacements/strains can be measured on a planar object, -from a practical point of view the camera and the object need to be aligned so that the camera image plane and the planar object be parallel and remain parallel all along the experiment, -if the object undergoes some out-of-plane displacements they will not be properly detected and they will produce in the images false apparent strains [SUT 08b].

With the stereo-correlation technique (also called 3D-DIC) that is described hereafter all these problems can be tackled. This technique allows to measure 3D displacements (in-plane and out-of-plane) and 2D strains undergone by any 3D object (not necessarily planar) while the experimental constraints are minimized [ORT 09a].

It should be noted that stereo-correlation procedures combine two techniques: the DIC technique, which is an image matching technique that allows to compute 2D displacement of pixels between two images. With this technique, several images taken at different instants with a single camera can be registered (see Section 6.2), or two images (or more) taken at the same instant by two cameras (or more) can be matched (see below).

the stereovision technique, which is a triangulation-based 3D reconstruction technique, that allows to compute the 3D position of a scene point from its stereo projections in two images (or more).

The stereovision technique

Binocular stereovision is a technique for recovering the three-dimensional structure of a scene from two (or more) different viewpoints (Figure 6.6a where P is the 3-D point to be measured, p 1 and p 2 are its stereoprojections in the images, C 1 and C 2 are the optical centers of the two cameras).

From Figure 6.6a, we can see that it is possible to compute the 3D coordinates of point P provided that: a) the two image-points p 1 and p 2 , which correspond to the projection onto the images of the same physical point P, can be identified. This step is called stereo image matching or search for stereo-correspondents. It is the critical step of the stereovision technique (see Section 6.4.1.3). b) the two lines C 1 p 1 and C 2 p 2 that intersect at P can be computed. This step requires that the intrinsic parameters (focal length, size of the pixels, distortion coefficients) of each camera and the extrinsic parameters of the stereo rig (relative position and orientation of the two cameras) be known. The intrinsic and extrinsic parameters are obtained by using an off-line calibration step (see Section 6.4.1.2).

. Geometry of a stereovision sensor

A binocular stereovision sensor is made up of two cameras positioned in such a way that their field of view intersects. Let us write R c the reference frame associated to the left camera, R ′ c the reference frame associated to the right camera, and R w the so-called world reference frame in which the 3D measurements will be expressed. T, T ′ and T s are the rigid-body transformations (a rotation matrix and a translation vector) that link these reference frames (Figure 6.7). These three transformations are linked by: T s = T ′ T -1 . T s is often called the stereoscopic transformation.I t represents the orientation and translation of a camera with respect to the other one.

It should be noted that the 3D measurements are often expressed in the reference frame associated to the left or the right camera. In that case T or T ′ are the identity matrix.

m′ = K ′ T ′ M = K ′ T s T M (6.11)
The image registration problem mentioned in Section 6.4.1 consists in finding the image-points m and m ′ that are stereo-correspondent.

In stereovision, the epipolar geometry [HOR 95] illustrated in Figure 6.8 provides the following important geometrical property (also called epipolar constraint):

Given a point m in the left image, its corresponding point m ′ in the right image appears to be always lying along a line of the right image entirely defined by the coordinates of m. This line is called the epipolar line associated with m.

6. The pinhole model is based on the perspective projection. In such a model the distortions (which can be more or less important) induced by the lenses are not taken into account. For metrology applications, which require a high measurement accuracy, distortions need to be taken into account and more sophisticated camera models need to be used [GAR 01a, GAR 01b, COR 05b, ORT 09b, ORT 09a, SUT 09]. Thanks to this important geometrical property, inherent to any stereo imaging system, the search for the stereo-correspondent of a given point in the left image is simplified from a 2D search across the entire image to a 1D search along its epipolar line. Using this geometrical constraint the image registration can be quicker and more robust (since possible false registration that could be found away from the epipolar line are avoided).

Calibration of a stereovision sensor

Camera calibration is an important task in 3D computer vision, particularly when metric data are required for applications involving accurate dimensional measurements. Calibrating a camera involves determining its intrinsic parameters (matrix K, and possibly distortion parameters). Calibrating a stereovision sensor made up of two cameras involves determining the intrinsic parameters of each camera and the relative position and orientation between the two cameras (T s transformation). These calibration data are required to compute, by triangulation, the 3D coordinates of a point corresponding to matched pixels on the two images. For more details on the calibration of a stereovision sensor, see References [GAR 01a, GAR 01b, ORT 09b].

Matching of stereo images by DIC

This section is dedicated to the matching of stereo images by DIC (which leads to the so-called stereo-correlation technique). It should be noted that stereovision can be used without DIC-based matching. It is for instance the case when stereo images are matched using feature matching [ORT 02, ORT 09a].

DIC-based image matching simply consists in using the DIC technique described in Section 6.2 in order to find stereo-correspondents. The only significant difference is that stereo-correlation can use (it is not compulsory) the epipolar constraint in order that the search for correspondent pixels be restricted to the epipolar line (or to a pixel band around the epipolar line).

3D reconstruction by triangulation

If the image-points m =( u, v) and m ′ =( u ′ ,v ′ ) (provided by the registration) and matrices K, K ′ , T and T s (provided by the calibration of the stereovision sensor) are known, Equations (6.10) and (6.11) lead to an over-determined system of four equations and three unknowns that are the three coordinates of the searched 3D point M =(X, Y, Z) [HAR 97]:

⎛ ⎜ ⎜ ⎝ u v u ′ v ′ ⎞ ⎟ ⎟ ⎠ = ⎛ ⎝ H 4×3 ⎞ ⎠ ⎛ ⎝ X Y Z ⎞ ⎠ (6.12)
Equation (6.12) can be written as b = H M and the 3D coordinates of point M are determined analytically using the pseudo-inverse method

M = (H T H) -1 H T b
The triangulation problem can also be solved using non-linear optimization by minimizing the distance between the measured image-points (u, v) m and the imagepoints (u, v) p predicted by the camera model [SUT 09]

χ 2 = (u m -u p ) 2 +(v m -v p ) 2 +(u ′ m -u ′ p ) 2 +(v ′ m -v ′ p ) 2 (6.13)
Note that Equation (6.12) implies a linear model of the camera (pinhole model). When lens distortions are taken into account, the camera model is no longer linear. To compute the 3D coordinates of point M using Equation (6.12), the image-points need to be corrected for their distortion beforehand. If the distortion is not corrected beforehand, then Equation (6.13) needs to be used.

3D displacement measurement by stereo-correlation

Using the stereovision technique, the shape variation of an object can be measured by analyzing a sequence of pair of stereo images. However, in experimental mechanics, we are generally interested in the surface strain field, which can be obtained by tracking the displacement of some points at the surface of an object undergoing some mechanical or thermal loading as illustrated in Figure 6.9. The 2D technique described in Section 6.2 allows the 2D displacement field to be computed at the surface of a planar object. By combining the stereo-correlation technique (DIC-based spatial registration between two cameras and triangulation, to get a cloud of 3D points) and the DIC-based temporal matching between the images acquired by the left or the right camera at different instants, the 3D displacement field of any object (with any complex shape) can be computed. This is illustrated in Figure 6.10. This measurement technique provides the 3D displacements, from which the strains can be computed as described in the next section.

Computation of surface strains from 3D displacements

The stereo-correlation technique combined with temporal image matching leads to the 3D displacement field of any object with any complex shape (not necessarily planar). The surface strain field can be computed from the 3D displacement field by differentiation of the 3D displacement at each point of the surface, which can be a complex task.

From a practical point of view, the computation of the strain at a given point P of the surface is performed through several steps [COR 05a, SUT 09]:

-first, a small contiguous collection of 3D surface points (e.g., a 7 × 7 array with center-point P) is selected from the initial surface shape, -next, a least squares plane is fitted to this 3D position data set providing the tangent plane, -a local coordinate system (n, b, t) is associated with the tangent plane and centered in P, -the array of data points for each component of displacement is then fitted with a least squares functions of the form g(X, Y ), -differentiating the functional fit g(X, Y ) for each displacement component in both the X and Y directions leads to the ε XX (P), ε YY (P) and ε XY (P) components of the strain tensor at point P.

Applications

The stereo-correlation technique has proven to be a powerful non-contact technique for measuring 3D displacement / 2D strain fields on any 3D object surface and is now widely used for industrial applications in different fields of mechanics and materials [ (small or large structures), it can measure both small (a few hundreds of micro-strains) and large strains (>200%), it can be used for high-speed dynamic tests by using high-speed cameras [SCH 03a, SCH 03b, SIE 07, TIW 07, BES 08, BES 10].

Conclusions

Digital image correlation, as we have seen through its many variants, in two or three dimensions, surface and even volume, at different scales of space and time allows for the measurement of displacement fields with a rich description and remarkable accuracy especially in light of its ease of implementation. Image correlation can exploit with great flexibility pictures acquired by different imaging means every day better and more informative. Many tools are available to incorporate in this measurement procedure any a priori knowledge that is available on the kinematics to assess the measurement uncertainty, or the covariance matrix of the measured degrees of freedom induced by image noise, and finally to visualize the residual field, often carrying very valuable information. Today's experiments provide a number of measured data comparable or sometimes even higher than what traditionally was the privilege of numerical simulations. These measurements are used, and currently more often than never before, to monitor mechanical tests, whether in the context of conventional characterization, or when used in an industrial context. The latter distinction is increasingly blurred as kinematic measurements, in particular by stereo-correlation, can account for complex loadings and geometries. This merger opens up outstanding perspectives in terms of quantitative use of tests, and therefore relevant responses to the use of materials, structures, or assemblies under realistic conditions, representative of their in-service life.

Last, what explains the spectacular development of this technique is probably its coupling with numerical simulations for identification and validation purposes. The end of the book is devoted to this broad issue. Again the classical boundary between experimental and numerical fields gradually vanishes, thanks to measured fields. They thus appear in different dimensions as the focal point to a particularly fruitful dialogue between all actors interested in the mechanical behavior of materials and structures, experimentalists and modelers, designers and users, researchers and producers.

Figure 6 . 1 .

 61 Figure 6.1. Examples of natural image textures suitable for digital image correlation. a) 304L stainless steel (1 pix. = 3.2 µm); b) Silicon carbide (1 pix. = 1.85 µm); c) Topographic AFM image of a silica glass (1 pix. = 2 nm); d) Phase AFM image (same area as in c)

Figure 6 . 2 .

 62 Figure 6.2. Examples of artificial image textures suitable for digital image correlation

Figure 6 . 3 .

 63 Figure 6.3. Examples of sections through 3D X-rays computed tomography images; a) Cast iron with spheroidal graphite modules; b) Polymeric foam; c) Stonewool

Figure 6 . 4 .

 64 Figure 6.4. Three-dimensional residue field showing the concentration of errors on the support of the crack (after [RAN 10])

Figure 6 . 5 .

 65 Figure 6.5. Horizontal displacement component obtained on a silicon carbide specimen with (left) Q4 finite elements (of size 32 pixels) or an integrated approach with few closed-form solutions (right). The physical size of one pixel is 1.85 µm. The displacements are expressed in pixels (after [ROU 09])

Figure 6

 6 Figure 6.6. (a) The principle of stereovision. (b) Picture of a stereo rig

  Figure 6.7. The three reference frames associated to a stereovision sensor

  7. To be able to write the pinhole model in a matrix form, homogeneous coordinates are used [FAU 93, FAU 01]. m are the homogeneous coordinates associated to m.

Figure 6 . 8 .

 68 Figure 6.8. The epipolar geometry: the epipolar lines are defined by the intersection of the image plane of each camera with the plane defined by points M, C and C ′ , where C and C ′ are the optical centers of the left and right camera respectively

Figure 6

 6 Figure 6.9. 3D displacement field corresponding to two deformation states of the object

Figure 6

 6 Figure 6.10. 3-D displacement field computation. Stereo-correlation (DIC-based spatial matching + triangulation) provides the shape of the object at a given instant. DIC-based temporal matching between the images taken by the left or right camera at different instants allows to temporally link the 3D points of the shapes

  SUT 09], namely, fracture mechanics [SUT 92, LUO 94, LUO 00, SUT 07b], biomechanics [SUT 08a, TYS 02], for the study of different types of materials, namely, composites [MCG 01, LOP 04, MUL 05, KAR 06], ceramics / concrete [ORT 07, ROB 07], foams [GUA 07], elastomers [MIS 05, JON 06], tyre [MOS 07], soft membranes [VIA 05]. The technique is very versatile and can be used for a large scale of experimental mechanics problems. It was applied at micro [SCH 04] or macro scales [HEL 03]

By convention, one considers here that the similarity is larger when C is smaller, although the reverse could be adopted (by changing the sign of C for example).

For the sake of ease of notation, the various arguments of C are omitted in the following definitions.

If these changes can be quantified independently, it is possible not to carry out optimization, and to apply the criterion C 1 to the image f and the rescaled image ag+ b.


Jérome Van Der Woerd

Chapitre 11 Déformation cumulée, vitesse long-terme et cycle sismique des décrochements intra-continentaux

Introduction

L'un des objectifs de la morphotectonique est de caractériser la récurrence des séismes le long d'une faille afin de contraindre des modèles de risque sismique et d'estimer les probabilités de taille et d'occurrence des futurs séismes. Connaître le passé du fonctionnement sismique d'une faille c'est un des moyens de prévoir son comportement futur (cf. Chapitre 12). Ainsi, classiquement on associe la vitesse moyenne de déplacement d'une faille sur le temps long et une certaine taille de séisme pour élaborer des modèles de glissement qui ont un certain caractère prédictif en terme de séquence possible de séismes. La vitesse d'une faille est la mesure du taux de déplacement d'un côté de la faille par rapport à l'autre. Une vitesse géologique est en général déterminée par la mesure à travers une faille du décalage d'une forme du relief, accumulé sur une longue période de temps (plusieurs milliers d'années à plusieurs millions d'années). Une vitesse géodésique est déterminée par la mesure du déplacement du sol de part et d'autre d'une faille, et à une certaine distance de celle-ci, sur une période de temps courte (quelques années à quelques dizaines d'années) en utilisant le plus communément le système GPS ou l'interférométrie radar (Figure 11.1, cf. Chapitres 1, 14).

Si nombre de failles restent peu étudiées, l'avancée des connaissances a permis de progresser dans leur caractérisation. La question se déplace de : « Quelle est la vitesse de cette faille ? » vers : « Cette vitesse est-elle constante dans le temps ? Varie-t-elle au cours du cycle sismique ou à l'échelle de plusieurs cycles ? » (cf. Chapitre 1). Comparer des vitesses issues d'approches différentes (géodésie, géologie) nécessite de bien définir ce qui est effectivement mesuré. A cette condition on pourra discuter de l'origine possible des différences ou variations observées. Celles-ci peuvent être liées à des propriétés physiques de la faille, à l'évolution de la géométrie des failles et leurs interactions au cours du temps, ou bien encore à des changements des conditions aux limites lointaines imposées par la tectonique des plaques.

Ce progrès des connaissances est essentiellement un résultat de l'accumulation de données de vitesses obtenues avec des méthodes variées et en constante évolution de datation des marqueurs géomorphologiques, de paléosismologie, et de géodésie (GPS, InSAR) en de nombreux sites le long de failles à différents stades de leur cycle sismique. Pour les vitesses géologiques, les progrès en matière de datations par les méthodes basées sur les isotopes cosmogéniques (10Be/26Al/36Cl/21Ne) mais aussi par thermoluminescence (OSL) sont indéniables. Il faut également souligner ces dernières années l'apport considérable que représente l'imagerie à haute résolution (issue des satellites Ikonos 1m, Quickbird 0.60m, SPOT5 2.5m [LI 05, MER 05] et la topographie radar à haute résolution (Lidar au sol ou aéroporté [ZIE 10, ZIE 15]) ou encore l'automatisation de méthodes photogrammétriques pour l'acquisition rapide de modèles numériques de terrain par drone en particulier.

Dans la suite nous nous concentrerons sur la définition et la mesure d'une vitesse géologique à partir de la déformation cumulée de marqueurs géomorphologiques. Nous aborderons certaines méthodologies couramment utilisées et discuterons certaines de leurs limites. Nous discuterons de la relation entre le décalage cumulé et le séisme, le décalage cumulé résultant de l'addition de déplacements sismiques au cours du temps. Enfin, nous avons choisi de discuter seulement des failles décrochantes.

Du décalage géomorphologique à la vitesse d'une faille

Décalage tectonique des rivières

L'idée que l'on puisse faire une mesure quantitativement significative de certaines propriétés des séismes et des failles à partir de ruptures de surface est née de l'observation des déformations suite à certain grands séismes, notamment le séisme de San Andreas du 18 avril 1906 (cf. Chapitre 1, [ZOB 06]). Non seulement il a été possible de mesurer une longueur de faille (rupture), mais également de caractériser la distribution du glissement co-sismique le long de la rupture et de quantifier la déformation des blocs séparés par la faille. Bien que certains segments aient été cartographiés au préalable, il a fallu le séisme de 1906 pour faire la démonstration de la continuité de la faille de San Andreas sur des centaines de kilomètres. Un autre aspect clairement révélé par la juxtaposition de la rupture récente avec des morphologies héritées de déformations antérieures établissait le caractère récurrent des grands événements sismiques. Ces conclusions confortaient des observations similaires faites à la fin du XIXème siècle le long de la faille normale de Wasatch (Utah) du décalage répété de lignes de rivages abandonnées du paléo-lac de Bonneville par Gilbert [GIL 84].

Un exemple classique de marqueurs géomorphologiques décalés est représenté par le site de Wallace Creek le long de la faille de San Andreas (Figure 11.2) [START_REF] Wallace | Notes on stream channels offset by the San Andreas fault, southern Coast Ranges, California[END_REF][START_REF] Sieh | The san Andreas fault at Wallace Creek, San Luis Obispo County, California[END_REF][START_REF] Sieh | Holocene activity of the San Andreas fault at Wallace Creek, California[END_REF]. A cet endroit situé dans la plaine relativement aride de Carrizo, une rivière intermittente, le Wallace Creek, prend sa source dans les montagnes de Temblor au nord-est. Cette rivière incise le piedmont en traversant la faille de San Andreas formant un chenal de plusieurs mètres de profondeurs. La vitesse moyenne de mouvement de la faille de San Andreas à ce site est contrainte à 33.9±2.9 mm/an pour les derniers 3700 ans et 35.8±5.4 mm/an pour les derniers 13250 ans, en bon accord avec des mesures de géodésies de nivellement et avec des mesures de géodésie spatiale GPS (35.9±0.5 mm/an [MEA 05]).

Les bords de terrasses alluviales (risers) pour la détermination des vitesses moyennes

Dans la suite nous nous focalisons sur les bords de terrasses comme marqueurs de la déformation le long des décrochements. En effet, ce sont des marqueurs très répandus et relativement nombreux le long des décrochements. Dans de nombreuses régions de tectonique active en convergence oblique la déformation se partitionne entre décrochement et chevauchement sur des failles distinctes, le décrochement longeant une chaîne de montagne. Cette géométrie des structures entraine presque automatiquement la formation d'un drainage perpendiculaire à la direction du décrochement et c'est dans ces environnements que l'on peut espérer trouver de nombreux marqueurs morphologiques décalés.

Ces environnements sont idéaux pour quantifier les déplacements récents le long des décrochements pour plusieurs raisons : -le drainage perpendiculaire (ou quasi) permet de reconstruire avec une bonne précision les marqueurs géomorphologiques naturels décalés -ces marqueurs sont associés au dépôt de sédiments et sont donc potentiellement datables -les phases d'incision et de déposition des rivières entrainent la formation de marqueurs d'âges variés -la longueur des chaînes de montagne et des failles décrochantes (quelques centaines au millier de kilomètres) augmentent d'autant le nombre d'observations possibles.

Ces caractéristiques sont les points clés de l'approche morphotectonique (ou tectonic geomorphology) : mesure du décalage, datation du décalage, représentativité temporelle et spatiale des mesures.

La figure 11.4 montre schématiquement la formation d'une séquence de terrasses concomitante du déplacement le long du décrochement inspiré de situations rencontrées le long du segment de Xidatan de la faille du Kunlun [VAN 98, VAN 02, TAP 01, RYE 06]. La précision de la mesure des décalages horizontaux des bords de terrasses et les relations temporelles entre l'incision et l'abandon des différentes terrasses nécessitent de reconstituer la mise en place des différents dépôts et phase d'incisions. Ceci se fait par une approche en 3 dimensions à partir de données topographiques de précision. En effet, selon les situations, il faut préciser la relation entre l'incision du bord de terrasse et l'âge des terrasses adjacentes. Ceci est parfois rendu complexe lorsque certains niveaux de terrasses ne sont pas présents de part et d'autre de la faille en raison de l'érosion ou quand l'incision est progressive et les bords de terrasses diachroniques.

Variation spatio-temporelle de la vitesse long-terme des failles

Pour comprendre le rôle que joue les failles dans la déformation des continents il est nécessaire de connaître leur vitesse de mouvement en en particulier si cette vitesse change le long de la faille ou au cours du temps. Si à petite échelle une faille est un objet linéaire simple, à grande échelle un décrochement est une succession de segments avec des jonctions et des bifurcations parfois complexes. Le lien entre un segment et le système de failles entier n'est pas toujours simple. De plus les endroits où on peut déterminer une vitesse géologique ne sont pas répartit uniformément, les situations qui associent des décalages géomorphologiques bien contraints avec des opportunités de datation précise sont rares. Si à cela on ajoute que la vitesse est susceptible de varier dans le temps et spatialement le long de la faille, le bilan cinématique d'une faille devient complexe.

Figure 11.7. Deux approches de détermination de la vitesse moyenne d'une faille. a) Le décalage est bien déterminé, mais l'âge de son accumulation est mal connu et borné par l'âge des terrasses inférieures et supérieures. b) Cas différent où l'âge est bien déterminé (ici âge d'une terrasse) mais dont le décalage est moins bien déterminé et compris entre des bornes minimale et maximale [MER 12].

La vitesse de la faille du Karakorum dans l'ouest du Tibet a été déterminée à une dizaine d'endroit le long de la faille par la datation de décalages de dépôts alluviaux, de bords de terrasses et de moraines abandonnées [CHE 15]. A certains endroits la faille est double et nécessite d'intégrer des vitesses sur plusieurs segments. Selon les sites, la vitesse est déterminée pour des périodes variées allant de 10-15 ka jusqu'à 150 ka. Les vitesses déterminées varient de 2 à 8 mm/an, avec une majorité entre 3 et 6 mm/an. Si l'on additionne les vitesses aux endroits où la faille est double, en faisant abstraction du fait que les vitesses sont déterminées pour des périodes de temps différentes, alors la limite haute de la vitesse à 8 mm/an est validée. Si au contraire, on fait l'hypothèse de variation temporelle alors la vitesse présente de la faille serait plutôt autour de 3 à 6 mm/an avec une augmentation vers le sud. La comparaison des vitesses obtenues par des méthodes différents comme la géodésie, la géomorphologie et la géologie montre une certaine constance temporelle entre 3 et 10 mm/an de l'actuel à 30 Ma, mais la précision des méthodes et les incertitudes associées à la mesure des objets considérés autorise à y voir également des variations temporelles de la vitesse et peut-être une diminution de la vitesse au cours du temps.

Figure 11.8. Variation le long de la faille du Karakorum de la vitesse long-terme pour des périodes de temps différentes. Sommer des vitesses déterminées sur des périodes différentes présuppose que les vitesses ne varient pas temporellement.

Les limites de blocs continentaux ou des plaques tectoniques sont rarement des zones de faille simple. Si la faille de San Andreas accommode une large part du déplacement entre la plaque Pacifique et la plaque Nord Américaine en Californie du sud, d'autres failles jouent également un rôle. Ces failles peuvent se répartir jusqu'à plusieurs centaines de kilomètres de distance de la faille principale, définissant un système de failles complexe dont la géométrie évolue dans le temps avec le déplacement progressif des micro-blocs que ces failles délimitent. La vitesse faible sur certaines failles n'en fait pas moins représenter un risque sismique potentiel en particulier dans les zones très peuplées.

En Californie du sud, la vitesse relative des plaques Pacifique et Nord Américaine est d'environ 5 cm/an. La faille de San Andreas, puis plus au sud, la faille de San Andreas et la faille de San Jacinto, accommodent plus de la moitié de ce mouvement relatif des plaques. Le reste de la déformation est accommodée par d'autres failles, en général plus petites, mais capables de produire des séismes de magnitude supérieure à 7. Connaître la vitesse de mouvement respective de l'ensemble des failles est donc un moyen de comprendre comment la zone se déforme mais également de faire face aux aléas que représentent chacune des failles. Dans cette situation avec de nombreuses failles sub-parallèles proches, les méthodes géodésiques qui nécessitent de se placer à distance des failles ne résolvent pas correctement la manière dont la déformation se répartit sur chacune des failles. Lorsqu'on additionne le mouvement sur chaque faille le long de trois coupes perpendiculaires à la limite de plaque (Figure 11.9) on retrouve le déplacement relatif de 5 cm/an des plaques tectoniques. Dans le détail ces reconstructions permettent de discuter des incertitudes sur le mouvement sur certaines des failles. En particulier, à l'est de Los Angeles, les vitesses respectives de la faille de San Jacinto et de la faille de San Andreas sont débattues [ VAN 06]. S'il ne fait aucun doute que la faille principale de San Andreas a accommodé le plus de déplacement cumulé au cours des temps géologiques, il semble qu'à l'heure actuelle, la faille de San Jacinto serait plus rapide. Ainsi le système de limite de plaque évoluerait avec un changement des zones actives, la faille de San Jacinto prenant le relais et devenant la faille principale. Il n'est pas possible de dire si cette évolution actuelle est irréversible du point de vue géologique ou si elle n'est que le résultat de l'alternance de l'activité sismique le long des deux failles. Les observations de terrain suite au séisme mettent en évidence des variations rapides et d'un facteur 2 à 4 du glissement co-sismique lors de la rupture (Figure 11.11). Bien que la faille présente une géométrie relativement simple constituée d'une faille unique à l'échelle régionale, une partie de ces variations résultent de la complexité de la rupture sismique de la surface qui affecte une zone large autour du plan de faille et de la difficulté de mesurer le déplacement précisément (Figure 11.10). Cependant, il est probable que ces variations se résorbent dans le temps au cours de séismes successifs. Il est légitime de se demander si ces variations peuvent s'accumuler au cours des cycles sismiques et si elles pourraient biaiser la détermination de la vitesse de la faille au long-terme. Les datations par les isotopes cosmogéniques de certaines des terrasses décalées permettent de proposer un scénario de glissement sur les derniers 5000 ans qui implique, compte tenu des mesures variables des décalages co-sismiques et cumulés, une variation de la vitesse de glissement au cours du temps [WEL 04].

Répétition caractéristique des ruptures et des séismes

La compilation des mesures des décalages de la morphologie alluviale le long des décrochements permet d'accéder à l'histoire sismique passée de la faille. A condition que l'érosion soit modérée mais que des événements d'incision soit relativement fréquents par rapport à l'occurrence des séismes majeurs, alors les chenaux incisés au travers de la faille sont autant de marqueurs possibles de la déformation au cours du temps. Les plus petits décalages de la morphologie correspondent au décalage co-sismique du dernier événement sismique, et les décalages plus grands résultent de l'accumulation de déplacement lors d'événements antérieurs. Le long de segment de faille déterminé l'analyse de la distribution des décalages de manière statistique permet de mettre en évidence la répétition d'événements dont le glissement moyen est de même ampleur. Par exemple, le long du segment de faille de la Karakax sur le décrochement de l'Altyn Tagh, un ensemble de décalages de chenaux et de bords de terrasse est décalé de l'ordre de 6-7 m, déterminant ainsi le déplacement moyen du dernier grand séisme sur cette faille. Les décalages cumulés plus grand sont de moins en moins nombreux, comme on peut s'y attendre en raison de l'érosion, mais suggèrent des ensembles multiples Ainsi, les déformations cumulées le long des failles résultent d'incréments sismiques ponctuels, l'ensemble du déplacement des blocs a lieu sur un plan de faille. Il n'en est pas toujours ainsi lorsque plusieurs failles définissent une zone d'endommagement complexe, la déformation est alors accommodée sur une plus large avec des déplacements réduits sur chacune des failles.

Conclusion

La variabilité de la vitesse géologique d'une faille à différentes échelles est directement liée à la variabilité des mesures de déplacements le long d'une faille. Pour des périodes de temps longues la vitesse moyenne de déplacement est mieux déterminée mais le nombre de séisme et leurs tailles deviennent incertain, alors que pour les périodes plus courtes et récentes, la vitesse moyenne est moins bien déterminée et le nombre et la taille des séismes sont mieux déterminés. La combinaison des mesures à différentes échelles de temps, ainsi que l'augmentation du nombre d'observations par l'apport de techniques et d'outils en constante progression permettront d'obtenir des jeux de données de plus en plus fiable et significatif pour la compréhension du fonctionnement des failles.

Figure 11 . 1 .

 111 Figure 11.1. Répartition à la surface de la Terre des décrochements continentaux actifs d'au moins 1000 km de long, avec longueur et bornes des vitesses géologiques/géodésiques d'après la littérature.

Figure 11 .

 11 Figure 11.2. a) photo aérienne du site de Wallace Creek. b) vue oblique d'un modèle numérique de terrain à partir de Lidar aérien (https://opentopography.org/) révélant les détails du relief (données B4 25cm de résolution). La faille de San Andreas suit un escarpement topographique sur lequel naissent des chenaux régressifs. Wallace Creek, le chenal le plus large et profond en amont est décalé de manière dextre au passage de la faille. En aval (sur la gauche), un ancien chenal abandonné, plus large et aux bords plus émoussés représente la position antérieure de Wallace Creek.

Figure 11. 3 .

 3 Figure 11.3. Représentation schématique progressive de l'évolution du site de Wallace Creek depuis la mise en place des dépôts colluviaux de piedmont, l'incision du chenal de Wallace Creek et les phases d'accumulation de décalages et de réincision de la rivière intermittente. Les contraintes temporelles sont issues de datations au Carbone-14 des sédiments abandonnés dans les chenaux [SIE 84].

Figure 11 . 4 .

 114 Figure 11.4. Blocs diagrammes montrant une séquence de formation de terrasses alluviales recoupées par une faille décrochante. 1. Mise en place d'un cône alluvial T0 lors d'une décharge sédimentaire importante. La faille est enfouie. 2. La rivière creuse un lit T1. La surface T0 est abandonnée et commence à enregistrer les mouvements sur la faille, mais son bord est constamment rafraîchi par érosion latérale. 3. Durant une nouvelle phase d'incision, T1 est abandonnée et le bord de terrasse T0/T1, devenu marqueur passif, commence à enregistrer les déplacements horizontaux. L'âge de l'abandon de T1 date le décalage du bord de terrasse. 4. Des épisodes successifs de dépôts de terrasses et d'incision conduisent à la formation de terrasses dont les bords sont décalés variablement et de manière incrémentielle. 5. Situation similaire, mais avec une composante de mouvement vertical. Le décalage vertical ne s'accumule qu'au moment de l'abandon de la surface. Ainsi, dans cette séquence, le décalage vertical de T1 (ou T2) est corrélé avec le décalage horizontal de T0/T1 (ou T1/T2) et avec l'âge de T1 (ou T2) [VAN 02].

Figure 11 . 5 .

 115 Figure 11.5. Image Quickbird et interprétation géomorphologique de terrasses vdécalées le long de la partie centrale de la faille de l'Altyn Tagh [MER 12]. Dans l'exemple de la faille de l'Altyn Tagh, certaines terrasses ont été érodées complètement en amont de la faille (Figure 11.5), ainsi la géométrie des terrasses change drastiquement entre l'amont et l'aval de la faille et enfin, certaines phases d'incision entre les niveaux de terrasses sont importantes (> 5 m). La reconstruction 3D de la mise en place des dépôts permet de mieux comprendre les relations entre l'incision et le mouvement sénestre sur la faille (Figure 11.6). En particulier, on peut remarquer que la largeur du lit actif de la rivière change de 120 m à moins de 20 m au cours de l'Holocène en s'incisant de plus de 15 m. Dans ce cas où les terrasses ne sont préservées que d'un seul côté de la faille, seuls des bornes du décalage du chenal principal sont déterminées pour estimer la vitesse de mouvement sénestre plutôt que les bords de terrasses dont l'âge est moins bien contraint (Figure 11.6).L'interprétation des relations entre incision et décalage n'est pas simple et les objets géologiques sont complexes. Des études contradictoires ont conduit à discuter la part d'interprétation dans les incertitudes liées aux déterminations de vitesse longterme (epistemic uncertainty [SAL 15]). Il est clair que considérer le décalage d'un bord de terrasse qui n'est jamais daté directement est une source d'incertitude pour la détermination des vitesses long-terme. En particulier, la possibilité qu'un bord de terrasse, bien que séparant des terrasses de même âge puisse avoir un âge différent de part et d'autre de la faille a été considéré comme une limite à l'utilisation des bords de terrasse pour caractériser un décalage tectonique [LEN 64, LEN 68, SUG 60], la préférence étant mis sur le décalage de chenaux incisés à l'intérieur des terrasses. La combinaison de méthodes de datations variées, datation au Carbone-14,

Figure 11 . 6 .

 116 Figure 11.6. Evolution complexe irrégulière de l'interaction entre les mouvements tectoniques horizontaux et la dynamique de la rivière, qui incise son lit en ajustant la largeur de son lit. Ici, 4 terrasses majeures sont formées lors des variations climatiques des derniers 14 000 ans [MER 12].

Figure 11

 11 Figure 11.9. A l'échelle de la limite de plaque, le système de faille de San Andreas avec ses failles associées accommode le déplacement lointain de l'ordre de 5 cm/an des plaques tectoniques Pacifique et Amérique du Nord. On peut voir comment se répartit la vitesse sur les différentes failles le long de 3 coupes perpendiculaires à la limite de plaque. Le long de la coupe 3, il a été proposé que la variabilité des vitesses sur les failles de San Jacinto (SJF) et San Andreas (SAF) résultent de l'alternance de l'activité sismique le long de celles-ci.

Figure 11 .

 11 Figure 11.10. Vue vers l'ouest le long de la rupture du séisme de magnitude Mw7.9 du 14 novembre 2001 de Kokoxili le long du décrochement du Kunlun (Tibet). La silhouette de Paul Tapponnier rend compte de la largeur de la rupture et de la taille des déformations co-sismiques.

Figure 11 . 11 .

 1111 Figure 11.11. Répartition du glissement co-sismique le long de la faille du Kunlun lors du séisme de magnitude Mw7.9 du 14 novembre 2001, tel qu'il a été déterminé par différentes méthodes. En rouge, position du site d'étude décrit plus loin.Quel que soit le modèle de glissement sismique considéré (glissement caractéristique, prédictible en temps, etc… [SCH 84, SIE 96, TAP 01]), il est généralement admis que le déplacement des blocs loin de la faille finit par se traduire par la même quantité de déplacement sur la faille. Alors que la détermination de la vitesse long-terme sur des périodes de temps largement supérieures au temps de retour des séismes (plusieurs milliers d'années) montre une relative constance sur des sites distants le long de la faille du Kunlun [VAN 02], on peut donc supposer que la faille accumule le même déplacement en n'importe quel point de la faille. Si cela est vrai pour les déplacements cumulés, il semble que cela ne soit plus vrai pour des mesures de décalages et temporelles qui se rapprochent des valeurs du déplacement co-sismique et du temps de récurrence des séismes : il y aurait donc un lissage des disparités de glissement au cours du temps et au cours des

Figure 11 . 12 .

 1112 Figure 11.12. Cartographie des terrasses décalées et mesures des décalages à partir d'observations de terrain, d'analyse d'images satellitaires de haute précision et de levés topographiques de terrain. Projections de profils topographiques levés parallèlement à la direction de la faille, avec indication de la hauteur des terrasses et

Figure 11 . 13 .

 1113 Figure 11.13. Compilation des mesures de décalages géomorphologiques les plus petits et estimation du glissement co-sismique moyen.

  de 6-7 m, à 14, 19 et 24 m. Des observations similaires sont faites sur d'autres décrochements, comme par exemple la faille de San Andreas [ZIE 10, ZIE 15] ou la faille de Fuyun [KLI 11]. Ces observations des décalages suggèrent donc un comportement sismique où la quantité de glissement à chaque séisme majeur est semblable, caractérisant ainsi le comportement d'un segment de faille. Seuls l'addition d'un travail de datation permet ensuite de préciser la date d'occurrence des événements sismiques dont rien ne dit, a priori, qu'ils interviennent à intervalle régulier.

Figure 11 .

 11 Figure 11.14. Exemple de distribution des décalages morphologiques de chenaux et terrasses le long de la faille de la Karakax -Altyn Tagh, Tibet occidental. La courbe à droite ou fonction de densité de probabilité des décalages fait apparaître des déplacements co-sismiques moyens caractéristiques (ici multiple de 6-7 m).

  Il est possible de déterminer des limites maximales et minimales de l'âge du bord de terrasse. En effet, l'incision ne commence que lorsque le niveau supérieur est abandonné, et son décalage peut avoir commencé bien avant l'abandon de la terrasse inférieure. Cette approche reste valable si les terrasses ont des âges proches. On peut aussi ne considérer que l'âge maximum possible, c'est-à-dire l'âge de la terrasse supérieure conduisant à estimer des vitesses minimales [COW 07].ENCART : METHODE DE DATATION PAR LES ISOTOPES COSMOGENIQUESLes roches de la surface de le Terre sont bombardées par des particules cosmiques, comme des protons de haute énergie. Les interactions nucléaires entre ces particules et les atomes des minéraux des roches provoquent la formation d'isotopes radioactifs, dit cosmogéniques, qui vont rester piégés dans les minéraux. Ainsi, ces réactions provoquent la formation de Beryllium-10 ou d'Aluminum-26 à partir des atomes de Silicium et de l'Oxygène du quartz (SiO2). Le taux de production des isotopes cosmogéniques dépend du flux de particules cosmiques chargées à travers

l'atmosphère terrestre et d'une modulation par le champ magnétique terrestre, il dépend donc de l'altitude et de la latitude et doit être connu pour chaque site de mesure. Ceci est l'objet de programmes internationaux de modélisation et de calibration permettant une constante amélioration des mesures. Ce taux de production des isotopes cosmogéniques diminue également avec l'épaisseur de roche traversée et donc avec la profondeur. Quantifier l'accumulation de ces isotopes dans les roches de la surface et la sub-surface est donc un moyen de mesurer à la fois le temps d'exposition de la roche au rayonnement cosmique et de quantifier l'érosion subie par la surface analysée. On peut donc après une isolation chimique des éléments et l'utilisation de spectromètre à accélération mesurer la quantité d'isotopes accumulées dans des galets ou des blocs abandonnées par les rivières ou les glaciers dans des terrasses ou des moraines pour en connaître l'âge [GOS 01, BEN 14].