Abstract
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 inversion method is a nonlinear iterative method that alternatively reconstructs contrast sources and contrast function. One of the most outstanding advantages of this inversion method 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.
1. 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–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]). 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–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.
2. 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 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 . 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, 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 denotes the spatial variable, represents the source position, is the angular frequency, the source term 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 .
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 is the finite-difference matrix, or stiffness matrix, which is asymmetric and extremely sparse, is the unknowns vector, and 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 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 only once. However, a new matrix decomposition is required for another frequency.
3. 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 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 and incident field vector 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 , incident filed , and the scattered filed , the scattered fields satisfy where are the contrast source vectors, defined as in which the contrast function is given by
According to (7), the differential operator with the known Lame parameters and only depends on the background mass density . The inverse of the operator 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 location can be formally presented by the following equation: and the object equation or state equation satisfies where is an operator selecting the scattered fields on the measurement surface among the scattered fields in the total computation domain, and the operator 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 and the object (or state) error are denoted by the following equations, respectively:
To establish the objective functional, we should firstly give the specific expressions of the inner product and -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 and is dropped in the remainder of this paper.
3.2. 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 in (12) and (13). On the contrary, it alternatively reconstructs contrast sources 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 represents the data equation error and 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 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, and , by minimizing the objective functional (18). The first step of FD-CSI method is to update contrast source vectors by the conjugate-gradient method with Polak-Ribiere search directions. While updating the contrast source, the contrast function is assumed to be . At the th step, the contrast sources are updated by the following formula: where is an update-step size and is the Polak-Ribiere search directions (see [29]) given by in which is the gradient of the cost function with respect to evaluated at the 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 with conjugate gradient method, the contrast function is updated by finding the minimizer of the functional . We assume the contrast source vectors 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), is the adjoint operator of the operator , which is defined as where and the symbol has the same meaning with , 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 -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 , since the regularization factor is only dependent on contrast function and . 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].
3.3. Computational Complexity Analysis
The highest computational cost part of FD-MRCSI is to calculate these operators and 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 and is done only once at the beginning of the inversion process and then reused at each subsequent step. The decomposition process takes operations, where is the number of the unknowns of the discrete system (, is the number of the grid cells) (see [25]). The back-substitution needs operations for each source position. In the framework of FD-MRCSI, back-substitution is required twice for each of the two operators and per iteration. Hence, the computational complexity of FD-MRCSI for each frequency is approximately given by where 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 and scattered synthetic data 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 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 represents the LU decomposition computational cost for forward modeling and 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.
4. 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 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 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 and are the field data vectors with noise and the noiseless data, respectively, and are uniformly distributed random numbers on the interval , is a two- dimensional vector whose components are equal to one, and is the desired fraction of the noise.
4.1. Reconstruction in a Homogeneous Background Medium
The first example is concerned with recovering a square scatterer from noise observed data 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/cm3, 1.5 g/cm3 outside the scatterer. Hence, the contrast for the scatterer is 0.2. During the whole inversion procedure, the angular frequency is 8 Hz.
(a) True model
(b) Reconstruction result
(c) Discrepancy
(d) Model slice
The inversion domain is divided into uniform grids using the second-order scheme, adding the grids of PML; so, the size of the total computation domain is . For this experiment, and ; therefore, the grid dimension of the total computation domain is . 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/cm3. The large obstacle has a mass density of 2.7 g/cm3, and the small one has a mass density of 2.1 g/cm3. Hence, the contrasts for the two abnormities are 0.8 and 0.4, respectively.
(a) True model
(b) Reconstruction result
(c) Discrepancy
(d) Variation of the data equation error
In the inversion process, we discretize the total computation domain into 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 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.
4.2. 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, 1.6, 1.5, and 1.8 g/cm3 from top to bottom. Figure 3(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 Figure 3(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.
(a) Background medium
(b) True model
(c) Reconstruction result
(d) Variation of the data equation error
In the inversion, the contrast function , contrast sources , and state variables , , and 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.
4.3. 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 al. (see [8]). Firstly, use the previous inversion result as ; then, solve the forward problem with to get the total field vectors ; finally, the initial contrast sources are obtained with the equation .
We employ the strategy to rebuild the two test examples given in Section 4.1, and all relative parameters are the same as those present in Section 4.1. For the one scatterer model case, we successively select three frequencies (4, 6, and 8 Hz). For the two scatterers model case, the sequential frequency applied during the multiple frequencies inversion is 5 Hz, 7 Hz, and 10 Hz. Figures 4(a) and 4(b) show the corresponding reconstruction results. The relative errors are 8.13% and 12.60%, respectively, which are smaller than those in Section 4.1. Table 2 indicates the reconstruction relative errors of the two models using FD-MRCSI and FD-MRCSI with multiple frequencies technique.
(a) Multifrequency inversion result
(b) Multifrequency inversion result
From Figure 4 and Table 2, we can observe that the use of the sequential multifrequency approach can improve reconstructed model to some extent. This can be explained by the fact that a good initial guess can result in a better computation result for a locally convergent method.
5. Conclusions
In this paper, the finite-difference contrast source inversion method incorporating a multiplicative regularization factor is successfully employed in seismic full-waveform inversion to reconstruct the mass density. The method is very efficient because it does not require solving a full forward problem in each inversion iteration, and the highly efficient calculation of the forward modeling makes it have a potential for dealing with large scale computational inverse problems.
From the results of numerical experiments, the algorithm can well recognize the shapes of the scatterers and precisely identify their locations, no matter the background mediums are homogeneous or inhomogeneous medium. The algorithm is very robust when the measurement data is polluted by random noise. Furthermore, the FD-MRCSI using sequential multifrequency technique can successfully improve the reconstruction results, which is very hopeful for getting higher resolution in exploration seismology. Estimating other subsurface parameters using this method, such as Lame coefficients and velocity, is our future work.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors are grateful to Professor van den Berg who gives us many helpful suggestions on both theory and numerical implementation of the contrast source inversion method. This work is supported by the National Natural Science Foundation of China, under Grants no. 91230119, no. 41474102 and no. 11301119.