Near subsurface density reconstruction by full waveform inversion in the frequency domain

. The work proposed is part of a global project dealing with the characterization of heterogeneous media by both electromagnetic and mechanical full waveform inversions. Indeed Full Waveform Inversion of seismic reﬂec-tion or Ground Penetrating Radar data is an e ﬃ cient approach to reconstruct subsurface physical parameters with high resolution. This paper focuses on the mechanical part, and more speciﬁcally on quantitative imaging of nearsurface density. Processing ﬁeld data is challenging because the nature of the source and the sensors used impact the signal-to-noise ratio as well as the frequency range appearing in the recorded data. From then it becomes interesting to process the data in the frequency domain and work on a few representative frequencies of the recorded temporal signal. In this article, ﬁeld data are simulated by noisy synthetic data. Di ﬀ erent frequency strategies are used and their results are compared with each other. The inverse problem consists in assessing the density in the probed medium from the data on the displacement ﬁeld measured at the detectors. Such a problem is known to be nonlinear and ill-posed. It is solved iteratively by a regularized Gauss-Newton algorithm, which relies on the Fréchet derivatives obtained through the generalized reciprocity principle equivalent to the well-known adjoint method. The numerical results show an optimal strategy, for which the convergence rate and the computation time are reasonable, the spatial resolution is improved and the density is well reconstructed.


Introduction
Seismic reflection (SR) aims to spatially reconstruct density and Lamé parameters or Compressionnal and Shear elastic wave speeds of an heterogeneous subsurface.The physical configuration consists of an incident wave radiated by a transmitter positioned on the air/ground interface [1].This wave interacts with the interface and the heterogeneities within the subsurface.The resulting displacement field is measured at a receiver also positioned on the air/ground interface but laterally offset from the transmitter.
From the knowledge of the excitation and the measurements at several receivers, the aim is to establish as precise as possible a mapping of the mechanical properties of the investigated medium.This requires a realistic forward model that takes into account all the interaction between the incident wave and the probed medium.Here it is built by the finite element method.The latter is nested in a robust full waveform inversion (FWI) strategy that processes the data to iteratively reconstruct the mechanical properties.There are different methods to exploit the data in order to find the mechanical parameters of the medium which is seen as a perfect elastic material: we can mention in particular the global [2] or local [3] optimization methods and the methods resulting from the recent development of deeplearning [4].However, the approaches related to deep-learning and global optimization are not really adapted to our goal.Indeed, the deep-learning ones require a simple representation of the medium to be reconstructed, and a training phase that is costly in terms of computing time and resources.In the same way, the global optimization is suitable when there are only few unknowns which is not the case here since density and Lamé parameters are evaluated at each node point of the medium mapping.Among the local optimization methods FWI is of major interest for the case under consideration.Although initially developed in the 80's [5], FWI has known a huge development relatively recently, namely since the 2010's with the multi-parameter inversion [6], mainly due to the fast increase in processing and storage capacities of computers.Compared to other approaches, FWI is not costly in terms of computing time and in memory resources which is a great advantage.By exploiting all the information embedded in the data, FWI allows to get a fine reconstruction of the ground without knowing a priori neither the shape nor the parameters of the burried objects.Obviously, it is still needed to have an initial guess of the reconstructed medium to start the inversion process.In our study we have chosen a homogeneous medium.FWI can be formulated equivalently in time and frequency regimes.The second one is chosen here since it better handles the data according to their frequencies.Moreover, it can be applied to both frequency and time data after being Fourier transformed.Inverting SR data is a non linear and ill-posed inverse problem.This difficulty can be overcome by local optimization approaches such as the Quasi-Newton (L-BFGS) [7][8][9] and the regularized Gauss-Newton (RGN) methods.In this work we select the RGN algorithm, because it is accurate, simple to implement and yields a reasonable RAM management, as shown in quantitative imaging of the relative dielectric permittivity [10].There are different ways to use the data according to their frequencies in the inversion process, which we refer to as "strategies".The contribution of this work is thus twofold: firstly, to show the ability of RGN algorithm for density reconstruction; secondly, to compare four strategies using frequency data in different ways, and then to discuss the strengths and weaknesses of each one in order to determine an optimal one.This paper is organized as follows: after the introduction, the study configuration and the forward model are presented in Section 2. Section 3 is devoted to the description of the inversion process.The numerical results are shown and discussed in Section 4. Finally, the conclusion is given in Section 5.

Study configuration and forward model 2.1 Forward model
In a fixed Cartesian coordinate system (O,x,y,z), a two-dimensional medium with an invariance axis along Oz is considered, as shown in Figure 1.The medium is heterogeneous in the S domain and homogeneous around it.It is described by its density ρ(r) depending on the spatial position r = x x+y y (x and y spanning the entire S domain).Lamé parameters λ and μ are assumed to be constant everywhere.These three parameters are assumed to be frequency independent in the frequency range of SR.We consider an harmonic time dependence for the displacement field.A ponctual force density generates an incident wave at radial frequency ω k : it is modeled by the source term f k m δ(r − r m ) = f m e −j ω k t δ(r − r m ) y of amplitude f m = 1 oriented along the Oy axis and positioned at point r m = x m x + y m y (schematized with a triangle in Figure 1), with j the imaginary number such that j 2 = −1.Here δ corresponds to the Dirac delta function.The study is thus led in the context of plane deformation configuration, for which the displacement field u is entirely represented in the plane (O,x,y) and satisfies the following equation [11]: For the sake of clarity, the time dependence is omitted in equation ( 1).Symbols • and ∧ respectively designate the dot and cross products, notation ∇ refers to the nabla operator.Solution u(r) = u(r, r m , ω k ) calculated at point r depends implicitly on the position r m and the radial frequency ω k of the source term.It is obtained by solving equation ( 1) by the finite element method, with the module Structural mechanics of the commercial software Comsol Multiphysics®.

Study configuration
Figure 1 illustrates the specific configuration under study: the light grey area corresponds to the heterogeneous zone which extends over a 5 m wide by 4 m deep rectangle.The outer homogeneous part uses perfectly matched layers to avoid reflections at the boundaries of the domain.The same forward model is used to generate the data and build the inversion process.However, to avoid as much as possible the inverse crime [12], different meshes are involved: one for the inversion process and one for the reference data.A white Gaussian noise is added to the latter with a signal to noise ratio of 30 dB.The sources are placed on the air/ground interface and the common mid-point and common offset acquisition modes are considered for the measurements.Twelve transceivers (schematized with a cross in Figure 1) are used, since the displacement field cannot be measured at the position where the force is applied (l m) and due to reciprocity theorem, 66 separate measurements per frequency are obtained.
The actual medium to be reconstructed is composed of two inclusions, as shown in Figure 1.One is a disk with density ρ = 2000 kg.m −3 corresponding to sandstone material.The other is a square made of clay of density ρ = 2200 kg.m −3 .The karstic medium surrounding these two inclusions has a density of 2600 kg.m −3 .The Lamé parameters are set to λ = 2.1437 × 10 10 N.m −2 and μ = 1.56065 × 10 10 N.m −2 in the whole inversion domain including the inclusions.All these values are chosen thanks to the literature but also thanks to travel time experimental data measured at EMMAH Physics laboratory.Such a configuration aims to demonstrate the ability of the inversion algorithm to deal with the low sensitivity of the data with respect to density and to achieve a good spatial resolution.

Inversion scheme
At each step n of the regularized Gauss-Newton iterative scheme, the correction on the density parameter is calculated in order to minimize the following normalized cost function: with the notation ||.|| 2  2 representing the square of the L 2 norm, u obs (r l , r m , ω k ) the displacement field generated by the actual medium, u n (r l , r m , ω k ) the displacement field simulated at the same iteration, N S the number of transceivers and Ω the set of frequencies involved.This correction updates the density as follows ρ n+1 = ρ n + δρ n , and ensures the decay of the cost function between two successive iterations F (ρ n+1 ) < F (ρ n ).It is obtained by solving in the least squares sense the following linear system [13]: The above integral is computed using the adjoint method and corresponds to the application of the Fréchet derivatives on the density correction.It implies the displacement field u n (r, r m , ω k ) solution of (1), but also the displacement field u n (r, r l , ω k ) solution of the adjoint problem, defined by considering the same medium but by positioning fictitiously the transmitter at the receiver position r l .
Consequently the whole linear system to be solved is formed by concatenating the different Fréchet derivatives taken for each measurement point and for each frequency as well as the differences between the simulated and measured fields.The LSQR algorithm [14] is applied to solve the above linear system, with its maximum number of iterations fixed at 10 for regularization purpose.The inversion process is driven by a homemade Matlab code which repetitively calls upon the software Comsol Multiphysics® to manage the forward model part.

Frequency strategies
To get an optimal reconstruction, a comparison of 4 strategies involving 7 different frequencies is considered.The choice of the frequencies (400, 445, 508, 574, 645, 750, 835, 900 Hz) based on [15] avoids the so-called "cycle skipping" phenomenon.These strategies are described as follows : 1. the Bunks Strategy consists in using simultaneously frequency packets according to the following sequence [16] : ( f 1 ), ( f 1 , f 2 ), ..., ( f 1 , .. , f 8 ).The result of the first packet is used as initialization of the new one.
2. the Group Strategy consists in the use of frequencies by packets of two according to the sequence ( f 1 , f 2 ), ( f 2 , f 3 ), ..., ( f 7 , f 8 ).The result of the previous frequency group is used to initialize the new one.
3. the Sequential Strategy was proposed by Pratt and Worthington [17] and uses the frequencies one at a time in increasing order.Again, the result of the previous frequency is used as initialization of the new frequency.
4. the Simultaneous Strategy uses all the frequencies at once.
The convergence is reached when the difference of the cost function between two iterations becomes smaller than a quantity C : |F n − F n−1 | < C. Choosing C = 10 −4 ensures that the cost function does not evolve anymore, which means that all the information has been taken into account.

Numerical results
The medium under consideration is meshed with a 1 cm characteristic length, which leads to around 200,000 unknowns for the density map.Density is discretized with respect to its position in space ρ(r i j ) = ρ i j with i ∈ [1, N] and j ∈ [1, M].The comparative study of the 4 strategies focuses on the accuracy of the reconstruction, the convergence rate and the computation time.The accuracy of the reconstruction is evaluated by the correlation coefficient written as follows [18] : with ρend and ρtrue the average values of the densities of the final iteration ρ end and of the true model ρ true .Their associated standard deviations are respectively noted σ ρ end and σ ρ true .The correlation coefficient assesses the similarity rate between the true model and the reconstructed one at the end of the inversion process.The initialization step of the inversion algorithm involves a perfectly homogeneous medium of density 2600 kg.m −3 .We assume that the Lamé parameters remain constant and unchanged during the whole inversion process.The correlation coefficient is thus the same for all strategies at the initialization step.With such a medium, the smallest wavelength is given for 900 Hz and is equal to 4 m.All the results presented in this article were obtained on a computer with 128 GB of RAM and AMD Ryzen Threadripper PRO 3955WX 16-Core CPU.
Figure 2 shows the different density maps obtained by the four strategies.To qualitatively appreciate the likelihood of the reconstructions compared to the actual inclusions, we  2b and 2c).The other strategies seem to have more oscillations and blurred objects.The visual and thus qualitative impression provided by Figure 2 is confirmed by Table 1.Indeed, the Group Strategy leads to the maximum value of the density correlation coefficient.However, it is very expensive in computation time.The Sequential Strategy has a CC slightly lower than the one of the Group Strategy but is less expensive in computation time.As expected, the Simultaneous Strategy yields the lowest number of iterations to converge but suffers from a long computation time and a poor reconstruction.The Bunk Strategy has the largest computation time and is not very appropriate in our case.But if the comparison is restricted to a smaller area closely surrounding each inclusion, the Simultaneous Strategy shows the highest CC.However the reliable comparison remains in the whole inversion domain and not just nearby the two inclusions.We also tested the inversion process with a higher (25 dB) Signal to Noise Ratio, for which the Group Strategy is more robust than the others and remain the best strategy to use.However the results have to be taken in their context and nothing ensures that the results will be the same for a mechanical multiparameter inversion.Indeed in an analogous electromagnetic approach, a multiparameter inversion on both dielectric permittivity and electrical conductivity gave the Bunks Strategy as the optimal one.A specific study will be led to investigate which is the best strategy in a mechanical multiparameter inversion.Moreover these results could be geometry dependent so this study has to be made in some new configuration.

Conclusion
In this work a regularized Gauss-Newton algorithm was used to invert data at several frequencies.Different frequency strategies have been tested and compared with each other.Numerical results of reconstruction of the actual density show that the Group Strategy gives the best spatial resolution i.e. with the highest correlation coefficient.Then the Sequential Strategy reconstructs the parameter with a slightly lower resolution but with a much shorter computation time.Both strategies show that the full waveform inversion leads to a good reconstruction of the parameter of the two inclusions, letting guess their shapes, although they have a characteristic length smaller than the smallest wavelength in the homogeneous medium.As a perspective for this work, we plan to apply this full waveform inversion method on SR experimental data at laboratory (50 cm & 50 kHz) and geophysical (100 m & 200 Hz) scales.

¡¢Figure 1 .
Figure 1.Configuration of the numerical study

Table 1 .
Summary table of the results