Seismic Waveform Inversion Using the Finite-Difference Contrast Source Inversion Method

This paper extends the finite-difference contrast source inversion method to reconstruct the mass density for two-dimensional elastic wave inversion in the framework of the full-waveform inversion.The contrast source inversionmethod is a nonlinear iterative method that alternatively reconstructs contrast sources and contrast function. One of the most outstanding advantages of this inversionmethod is the highly computational efficiency, since it does not need to simulate a full forward problem for each inversion iteration. Another attractive feature of the inversion method is that it is of strong capability in dealing with nonlinear inverse problems in an inhomogeneous background medium, because a finite-difference operator is used to represent the differential operator governing the two-dimensional elastic wave propagation. Additionally, the techniques of a multiplicative regularization and a sequential multifrequency inversion are employed to enhance the quality of reconstructions for this inversion method. Numerical reconstruction results show that the inversion method has an excellent performance for reconstructing the objects embedded inside a homogeneous or an inhomogeneous background medium.


Introduction
Full-waveform inversion (FWI) is a powerful tool to get high-resolution parameter models in the context of seismic tomography.This feature of FWI makes it very suitable for estimating subsurface parameters, which is the ultimate goal in exploration seismology.In recent years, with the development of computer technology and the requirement of higher quality of seismic imaging, FWI becomes a very promising technique for exploration seismology as acquisition improves in size and density.The process of FWI is to find an optimal parameter model by minimizing an objective functional that measures the discrepancy between the observed seismograms and the synthetic seismograms.In other words, the specific implementation of FWI is an iterative procedure that is usually based on gradient-based optimization, either in time domain or frequency domain (see [1][2][3][4][5][6][7][8][9]).In this paper, we only focus on frequency domain case.Simulation in frequency domain has its inherent advantages: it can obtain a relatively good reconstruction result with only a few frequencies, which means that it can significantly reduce the data redundancy; in addition, a parallel computation programming can be easily implemented for different frequency inversions, which results from the independence of each frequency.Most methods in the framework of FWI are based on gradient or Newton method.A common feature of these gradient-based optimization methods is the highly expensive forwarding calculation for estimating the gradient of the object functional with respect to the interested parameter.One special method to estimate the gradient for a minimization problem based on gradient method is the adjoint method (see [10]).Additionally, seismic inverse problems usually are nonlinear and illposed, which make the seismic inversion difficult to handle.In short, FWI incorporated both traveltime and amplitude information is a highly computational process, no matter it is implemented in time domain or frequency domain, and this feature greatly prevents from its widespread applications in exploration seismology.
An alternative method is proposed by van den Berg and Kleinman (see [11]), extended by van den Berg and Abubakar to effectively reconstruct the interested parameters in the field of the microwave imaging (see [12]).

Journal of Applied Mathematics
The method is referred to as the so-called contrast source inversion (CSI) method, which has two versions: IE-CSI based on integral equation formulation (see [11,12]) and FD-CSI based on the finite-difference operator (see [13]).The IE-CSI method is very efficient for the homogeneous or simple layered background medium inverse problems, since the IE-CSI method needs to calculate Green function of the background medium, while the FD-CSI method has a good performance for solving high complex problems in an inhomogeneous background medium.Another highly efficient scattering inversion method called subspace-based optimization method (SOM) shares this good property of CSI method, which also has two variants adopted in homogeneous and inhomogeneous background medium, respectively (see [14,15]).The comparisons in detail between SOM and CSI are covered in [14], and this reference also concludes that the CSI is close in spirit to a special case of SOM and the SOM could be convergent for a few special problems.The interesting advantage of CSI method is that it does not require solving any full forward problem in each inversion iteration; therefore, it can significantly reduce the computational cost, and it may be very efficient for large scale or high dimension inverse problems (see [16]).Both the IE-CSI and FD-CSI methods incorporating with multiplicative regularization constraint have been successfully applied in electromagnetic wave imaging and acoustic wave imaging (see [17][18][19][20][21]).
In this paper, we extend the FD-CSI method to reconstruct the mass density based on a two-dimensional elastic wave equation (in fact, all theories in this paper are suitable for three-dimensional problems without any changes) in an isotropic medium, either homogeneous or inhomogeneous background medium, and its IE-version is presented in [22], which is successfully employed in elastodynamics when the background medium is homogeneous.For our forward modeling, a perfectly matched layer (PML) absorbing boundary condition is employed to reduce the boundary reflections (see [23,24]), and, considering the fact that the stiffness matrix usually is large and highly sparse, the techniques of matrix compression store and LU decomposition of a sparse matrix are employed to decrease the storage space and increase the computation efficiency (see [25]).For the inversion processing, since all the finite-difference operators involving the FD-CSI algorithm are only dependent on the background medium and angular frequency that both of them are invariant throughout the inversion process, we use a direct solver such as an LU decomposition method to solve these finite-difference responses to improve the efficiency of the forward modeling.We invert these finite-difference operators only once at the start of the inversion processing, and the results can be reused for multiple source positions and successive iterations of the inversion.
In order to get a better reconstruction profile and some geometrical information of an abnormality, a multiplicative regularization factor called FD-MRCSI is incorporated, which is first introduced by van den Berg et al. (see [26]).The multiplicative regularization factor has the same function as the total variation (TV) regularization.The most promising advantage of the multiplicative regularization is that the regularization parameter is automatically calculated during each inversion.Furthermore, in order to mitigate the nonlinearity and enhance the reconstruction results, the sequential multifrequency inversion approach is employed; that is, the seismic inversion starts from low frequency and moves upward to high frequency.To illustrate the performance of the finite-difference contrast source inversion method in seismic inversion problems, we present several numerical experiments, either in a homogeneous or an inhomogeneous background medium.

Forward Modeling
In this section, we present the two-dimensional forward modeling using finite-difference method, which is of great importance for the inversion method.The forward modeling simulates the propagation process of the elastic wave in an isotropic inhomogeneous medium and drives the inversion algorithm.In time domain, the governing PDE system, called elastodynamic system, is written as the following differential equation system: where u is the displacement vector,  is the mass density,  denotes the stress tensor of the elastic medium, ∇⋅ represents the divergence operator with respect to spatial variables, and the external force is denoted by the term f.In an isotropic medium, the stress tensor  is given by the following formulas: where the indices  and  can be 1 or 2,   is the well-known Kronecker delta symbol,   denotes the strain tensor, and  and  are the Lame parameters of the isotropic medium.
For simplicity, all physical variables in frequency domain are also denoted by the same symbols as in time domain unless it makes an ambiguity.For example, u also represents the displacement vector in frequency domain.By substituting (2) and ( 3) into (1) and employing the Fourier transformation on both sides of (1), we get the Navier equations as where r denotes the spatial variable, r  represents the source position,  is the angular frequency, the source term f is a Ricker wavelet represented in frequency domain throughout this paper,  is a matrix transpose, and ∇ is a gradient operator with respect to the spatial variable r.System (1) and ( 4) defines the wave propagation in an infinite medium in time and frequency domains, respectively.However, the wave field is numerically modeled in a finite domain.Therefore, waves will reflect from the computation edges of the domain, and the artificial waves will be recorded.In time domain, these reflections can be eliminated by the time windowing; however, such reflections cannot be eliminated in frequency domain.Therefore, we should attenuate or suppress the unwanted waves to simulate a modeling in a finite domain, and one type of absorbing boundary conditions should be applied.For our forward modeling, a stabilized unsplit convolutional perfectly matched layer (PML) is employed (see [24]), and the source term f we added is in term of p-wave incident excitation.
The parameters and wave fields are uniformly discretized on a regular two-dimensional domain at intervals of Δ and Δ by using the second-order finite-differences described in [27,28].Note that the spatial sampling intervals Δ and Δ must be small enough so that the numerical solutions of the discretization linear equation system of ( 4) are sufficiently accurate.Generally, this five-point stencil requires about 10 grid points per shear wave wavelength in the lowest velocity zone for an accurate modeling (see Pratt [28]).After the unknowns at each grid point are ordered by row-major rule, we obtain a linear equation system as where A is the finite-difference matrix, or stiffness matrix, which is asymmetric and extremely sparse, U is the unknowns vector, and F denotes the discrete vector of the source term.To solve the linear equation system (5), we choose a direct solver method based on the sparse matrix LU decomposition rather than an iterative solver.The main reason we employ a direct solver method is that the stiffness matrix A is only dependent on the background medium and an angular frequency .Hence, we can get all solutions for all source excitations simultaneously by factoring the finite-difference matrix A only once.However, a new matrix decomposition is required for another frequency.

Finite-Difference Contrast Source Inversion Algorithm
3.1.Problem Formulation.We consider the two-dimensional elastic wave inversion problem in a homogeneous or an inhomogeneous background medium.Throughout the paper, the inversion domain or object domain in which the scatterers are embedded is denoted by symbol , and the data domain  is the surface where the sources and receivers are located.We assume that both of the domains are contained within the total domain  1 called finite-difference computation domain.By this configuration of this inverse problem, we can slightly reduce the inversion domain.Assume that there are   sources successively illuminating the inversion domain and the index for the th source is denoted by subscript .
For each source , the total field vector u  (r, ) and incident field vector u inc  (r, ) satisfy the following Navier equations, respectively: In ( 7), the parameter involving the subscript  denotes the relative parameters when the equation is employed to describe the seismic wave propagation in a known background medium.Subtracting ( 7) from ( 6) and using the relationship between the total filed u  , incident filed u inc  , and the scattered filed u sct  , the scattered fields satisfy where w  (r, ) are the contrast source vectors, defined as in which the contrast function (r) is given by According to (7), the differential operator N  with the known Lame parameters  and  only depends on the background mass density   (r).The inverse of the operator N  is formally given by the linear operator Therefore, the solution vectors of ( 8) can be symbolically denoted by the following equation: By using (11), the measured scattered field data vectors at the receiver r  location can be formally presented by the following equation: and the object equation or state equation satisfies where M  is an operator selecting the scattered fields on the measurement surface  among the scattered fields in the total computation domain, and the operator M  selects the fields inside the inversion domain .Note that ( 12) and ( 13) are the fundamental equations for the CSI method, which are used to establish the cost functional.From ( 12) and ( 13), the data error r   (r, ) and the object (or state) error r   (r, ) are denoted by the following equations, respectively: To establish the objective functional, we should firstly give the specific expressions of the inner product and  2 -norm as where the overbar denotes the complex conjugate for each component of a complex vector and symbol ⋅ represents the product of two vectors.Ω is the integral domain, which can be either  or .
For simplicity, the explicit dependence of the field quantities on r and r  is dropped in the remainder of this paper.

Finite-Difference Contrast Source Inversion Method.
In the framework of the contrast source inversion method (see [11]), it does not eliminate the contrast source quantities w  in ( 12) and ( 13).On the contrary, it alternatively reconstructs contrast sources w  and contrast function  by minimizing a cost functional.In order to avoid minimizing a nonquadratic objective function caused by the presence of contrast function  in both numerators and denominators, we borrow the idea from [14] and slightly modify the cost functional of CSI method.Therefore, the modified cost functional to be minimized at the th step is defined as the following form: where   (w  ) represents the data equation error and    (, w  ) represents the object equation error.The scattered field data is obtained by subtracting the simulated background field from the simulated model field.From ( 14), (15), and ( 18), we note that the background medium does not change throughout the inversion process; therefore, the computation cost of the inversion of the finite-difference operator N  is relatively cheap for each frequency, which results from the LU decomposition required only once.
The main idea of the CSI method is to reconstruct two interlaced sequences, w , and   , by minimizing the objective functional (18).The first step of FD-CSI method is to update contrast source vectors w , by the conjugate-gradient method with Polak-Ribiere search directions.While updating the contrast source, the contrast function  is assumed to be  −1 .At the th step, the contrast sources w , are updated by the following formula: where   , is an update-step size and k , is the Polak-Ribiere search directions (see [29]) given by in which   , is the gradient of the cost function with respect to w , evaluated at the ( − 1)th iteration.The details of this type of the conjugate-gradient can be found in [30].After updating the gradient   , , the update-steps   , can be obtained by minimizing the cost functional   ; that is, and the explicit expression of the minimizer   , of ( 21) can be founded by setting the derivative with respect to  being equal to zero; that is, where the normalization factors   and    are The second step of the CSI method is to update contrast function .Once we update the contrast source vectors w  with conjugate gradient method, the contrast function  is updated by finding the minimizer of the functional    (, w  ).We assume the contrast source vectors w  are constant during this process.It can be shown that the updating formula has the closed form as follows: where In (24), the symbol Re indicates getting the real part of a complex number.The details of finding the minimizer with respect to  of the objective functional (18) can be found in [11].
To start the FD-CSI algorithm, Considering the nonlinearity and illposedness of the inverse problem, and the local convergence of the FD-CSI method, a good initial guess is very important to get a better reconstruction result.We initialize the algorithm by the following way: and the contrast function is calculated by using ( 24) and (25).In (26), L *  is the adjoint operator of the operator L  , which is defined as where and the symbol M  * has the same meaning with L *  , which maps the field values on the measurement surface into the total computation domain.This initialization method is called back-propagation method, and it has been proved that it can provide a good initial guess (see [8]).
To improve the quality of reconstruction results, the multiplicative regularization technique introduced by van den Berg et al. (see [26]) is applied to this inverse problem, which has the same effect as the total variational regularized method.After incorporating the multiplicative regularization factor, the cost functional is rewritten as where the regularization cost function    () for the weighted  2 -norm regularization is defined as where where  is the area of the inversion domain and ΔΔ is the cell area.More details on finite-difference multiplicative regularization contrast source inversion (FD-MRCSI) method can be found in Appendix A of [8].
It is very similar to minimizing functional (29) using conjugate-gradient method.There is no change while updating the contrast source vectors w  , since the regularization factor    () is only dependent on contrast function  and    ( −1 ) ≡ 1.Nevertheless, the contrast function  is updated by using the conjugate-gradient method with Polak-Ribiere search directions.
Compared with an additive regularization method, it is not necessary to determine the artificial weighting regularization parameter because it is automatically calculated at each iteration.As we know that selecting the regularization parameter is a time-consuming and difficult work, especially there is no prior information about the measurement data, multiplicative regularization method is very suitable to invert experimental data as shown in [20].

Computational Complexity Analysis.
The highest computational cost part of FD-MRCSI is to calculate these operators L  and L *  to estimate the conjugate   and update-steps    .Hence, an efficient method to compute the action of these operators is important to increase the efficiency of the inversion algorithm.Considering the fact that these operators only depend on the background medium, which does not change throughout the inversion process for a given frequency , we choose an efficient direct method, that is, the LU decomposition, and the decomposition for each of the two discretization operators L  and L *  is done only once at the beginning of the inversion process and then reused at each subsequent step.The decomposition process takes ( 1.5 ) operations, where  is the number of the unknowns of the discrete system ( = 2,  is the number of the grid cells) (see [25]).The back-substitution needs ( log ) operations for each source position.In the framework of FD-MRCSI, back-substitution is required twice for each of the two operators L  and L *  per iteration.Hence, the computational complexity of FD-MRCSI for each frequency is approximately given by  (,  iter ,   ) ≈  (2 1.5 ) +  iter  (2   log ) , (32) where  iter and   are the number of inversion iteration and the number of sources, respectively.
For usual inversion methods (UIM), the objective functional involving the scattered observed data d  and scattered synthetic data u sct  is given in the following formula: The first term of (33) ensures that the minimizer   of  will be an optimal parameter model, while the second is the regularized term, which guarantees the stabilization of the problem and forces the minimizer   of (33) to satisfy certain regularity properties, and  denotes the regularization  iter is the iteration number,  S presents the number of source, and  denotes the number of the grid cells.
parameter that balances the data error and regularized term.
If functional (33) is also minimized by using a nonlinear conjugate gradient method, then the computational complexity of this algorithm for this inverse problem can be approximately estimated in analogy to FD-MRCSI method as follows: in which ( 1.5 ) represents the LU decomposition computational cost for forward modeling and ( log ) denotes the back-substitution computation for each source during each iteration process.Table 1 indicates the details of computational complexity for FD-MRCSI and usual inversion methods.
From Table 1, we observe that the FD-MRCSI method does not need to do an LU decomposition for each iteration, while usual methods require it, and this is the highest computational cost part.Therefore, the FD-MRCSI method is a greatly efficient method.

Numerical Results
In this section, we present some numerical examples based on the FD-MRCSI algorithm.All inversion domains of these tests are inside a square with side length .In order to measure the discrepancy between the exact profile and reconstructed profile, the relative error with respect to contrast quantity  is introduced, which is given by the following formula: where   is the reconstructed result after  iterations and  exact is the exact solution in which we are interested.Since the "inverse crime" phenomenon may happen, the synthetic data are generated by using a finite-difference time-domain forward method.A Ricker wavelet f with a dominant frequency 30 Hz is employed for generating the synthetic data.However, the inversion algorithm is carried out in frequency domain, and hence the datasets generated by using time-domain forward method are transformed into frequency domain using a fast Fourier transform (FFT) (see [31]).After that, we select some frequencies for the FD-MRCSI algorithm.Once we get the synthetic data, a random noise is added to each frequency synthetic data using this formula as follows: where d noise  and d  are the field data vectors with noise and the noiseless data, respectively,  Re and  Im are uniformly distributed random numbers on the interval (−1, 1), I is a two-dimensional vector whose components are equal to one, and  is the desired fraction of the noise.

Reconstruction in a Homogeneous Background Medium.
The first example is concerned with recovering a square scatterer from noise observed data d noise  in a homogeneous background medium, and the exact model is given by Figure 1(a).The inversion domain  is set as a square with side length 6 km, which does not include PML thickness.For this square scatterer, the second Lame parameter  is equal to 2.0 GPa, and a Poisson's ratio is 0.25.The square with the side length 1.5 km is located in the inversion domain, and its centre is at (3.3) km; the mass density of the scatterer is 1.8 g/cm 3 , 1.5 g/cm 3 outside the scatterer.Hence, the contrast  for the scatterer is 0.2.During the whole inversion procedure, the angular frequency  is 8 Hz.
The inversion domain  is divided into  ×  uniform grids using the second-order scheme, adding the  1 grids of PML; so, the size of the total computation domain is ( + 2 1 ) × ( + 2 1 ).For this experiment,  = 61 and  1 = 15; therefore, the grid dimension of the total computation domain is 91×91.Both the sources and receivers are assigned to 16, which are equally distributed on the four sides of the inversion domain, and all sources are added in the form of wave incident excitation.To show the denoising performance of the FD-MRCSI method, the numerical inversion results are presented when the measurement datasets are added 5% random noise.Figures 1(b) and 1(c) are the reconstruction profile of the experiment and discrepancy between the true model and recovering model after 1000 iterations.Figure 1(d) shows the slices along the depth at horizontal distances of 3 km for the true model and inversion model.The relative error is 12.26%.
In this type of numerical tests of a homogeneous background medium, we will test that the FD-MRCSI method possesses the ability of reconstructing multiple scatterers with different geometrical information and different physical property parameters.For this purpose, we employ a model consisting of two rectangular scatterers as shown in Figure 2(a).The configuration consists of two rectangular objects, which have different dimensions.One is 1.2 km × 1.2 km, and another is 0.8 km × 1 km.We also assume that the inversion domain is a square domain with the dimension of 6 km × 6 km and has a second Lame coefficient  of 5.2 GPa, a Poisson's ratio of 0.25, and a mass density of 1.5 g/cm 3 .The large obstacle has a mass density of 2.7 g/cm 3 , and the small one has a mass density of 2.1 g/cm 3 .Hence, the contrasts for the two abnormities are 0.8 and 0.4, respectively.
In the inversion process, we discretize the total computation domain into 151 × 151 uniform grids, including 15 PML grids, and 24 sources and 24 receivers are used for this example.These sources and receivers are distributed along the inversion domain, six sources or receivers for each side.After the measurement data is generated using the same way as the previous example, we select them at frequency 8 Hz to drive the inversion algorithm, and 5% noise is added.After 1000 iterations, the reconstruction profile is given in Figure 2(b).Figure 2(c) denotes the discrepancy between the true model and recovering model after 1000 iterations, and Figure 2(d) shows the variation of the data equation term error with respect to iterations.The relative error is 16.05%.
At first glance, we should solve the forward modeling for many times for these 1000 inversion iterations.Since the contrast source w and contrast function  are incorporated into this algorithm, which makes the forward stiffness matrix only involve the background medium for each given frequency, the LU decomposition is required only once throughout the whole inversion procedure; therefore, the forward modeling computation cost is very cheap.However, the price that we have to pay is the storing of the LU decomposition arrays of the stiffness matrix of the background medium.
From these numerical results, we observe that the FD-MRCSI method can effectively reconstruct the scatterers including their locations and shapes using the measurement scattered data, and the method is of strong denoising ability.Additionally, Figure 2(d) demonstrates that the FD-MRCSI method is convergent.

Reconstruction in an Inhomogeneous Background
Medium.In this section, we examine the performance of the FD-CSI using an inhomogeneous medium as a background medium.The configuration of the background consists of a compressional wave speed of 2500 m/s, a shear wave speed of 1443 m/s, and a four-layer distribution of a mass density of 1.8,For the elliptical scatterer, it has a semimajor axis of 0.7 km. and a semiminor axis of 0.36 km, and its centre is (2.3, 2) km.We also assume that the inversion domain is a square with dimension of  × .In our numerical experiments,  is assigned to 6 km.In the inversion, the contrast function , contrast sources w, and state variables u, u inc , and u sct in the inversion domain are uniformly discretized on a regular two-dimensional grid at intervals of Δ and Δ.Both Δ and Δ are equal to 0.05 km; in other words, the grid size is 121 × 121.The measurement domain surface  where 24 receivers and 24 sources are distributed is a square around the inversion domain.After the synthetic datasets are generated by using FFT on the finite-difference time-domain (FD-DT) data at frequency 10 Hz, 5% noise is added by using (35).After running 1000 iterations, the reconstruction results are given in Figure 3(c).The relative errors are 18.83%.
As shown in Figure 3, the FD-MRCSI method can successfully reconstruct the objects including their locations and shapes as those present in a homogeneous background medium.The feature of the inversion method results from the flexibility of the finite-difference operator of the governing differential equations.

Reconstruction Improvement Using Multiple Frequencies.
In these examples, we use the previous lower frequency inversion reconstruction results as an initial guess for another frequency inversion, which is the so-called sequential multifrequency inversion approach, that one usually employs this technique in the Gauss-Newton or nonlinear conjugate gradient method.This technique applied to improve the inversion profiles is guaranteed by the following two observations: (1) each frequency is independent of FD-MRCSI algorithm; (2) the inversion algorithm is locally convergent, and hence a good initial guess will result in a better inversion results.Additionally, multiple frequencies inversion technique can migrate the nonlinearity of the elastic inverse problem.
To drive the inversion algorithm, we should use another initialization method, which was introduced by Abubakar et

Figure 1 :
Figure 1: Applying the FD-MRCSI method to one abnormity inversion.(a) True model.(b) Reconstruction result after 1000 inversion iterations.(c) Discrepancy between true model and reconstruction result.(d) Slices at horizontal distances of 3 km for the true model and reconstruction result.
Variation of the data equation error

Figure 2 :
Figure 2: Applying the FD-CSI method to two abnormities inversion.(a) True model.(b) Reconstruction result after 1000 inversion iterations.(c) Discrepancy between true model and reconstruction result.(d) Variation of the data term error with respect to iterations.1.6,1.5, and 1.8 g/cm 3 from top to bottom.Figure3(a)shows the mass density distribution of the background medium.The inversion model consisting of a circle and an ellipse is presented as shown in Figure3(b).The circular abnormity with a contrast of 0.8 has a radius of 0.550 km, which is centered at (3.4, 3.6) km.For the elliptical scatterer, it has a semimajor axis of 0.7 km. and a semiminor axis of 0.36 km, and its centre is (2.3, 2) km.We also assume that the inversion domain is a square with dimension of  × .In our numerical experiments,  is assigned to 6 km.In the inversion, the contrast function , contrast sources w, and state variables u, u inc , and u sct in the inversion domain are uniformly discretized on a regular two-dimensional grid at intervals of Δ and Δ.Both Δ and Δ are equal to 0.05 km; in other words, the grid size is 121 × 121.The measurement domain surface  where 24 receivers and 24 sources are distributed is a square around the inversion domain.After the synthetic datasets are generated by using FFT on the finite-difference time-domain (FD-DT) data at frequency 10 Hz, 5% noise is added by using (35).After running 1000 iterations, the reconstruction results are given in Figure3(c).The relative errors are 18.83%.