Research Article  Open Access
Seismic Waveform Inversion Using the FiniteDifference Contrast Source Inversion Method
Abstract
This paper extends the finitedifference contrast source inversion method to reconstruct the mass density for twodimensional elastic wave inversion in the framework of the fullwaveform 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 finitedifference operator is used to represent the differential operator governing the twodimensional 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
Fullwaveform inversion (FWI) is a powerful tool to get highresolution 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 gradientbased 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 gradientbased 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 socalled contrast source inversion (CSI) method, which has two versions: IECSI based on integral equation formulation (see [11, 12]) and FDCSI based on the finitedifference operator (see [13]). The IECSI method is very efficient for the homogeneous or simple layered background medium inverse problems, since the IECSI method needs to calculate Green function of the background medium, while the FDCSI method has a good performance for solving high complex problems in an inhomogeneous background medium. Another highly efficient scattering inversion method called subspacebased 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 IECSI and FDCSI 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 FDCSI method to reconstruct the mass density based on a twodimensional elastic wave equation (in fact, all theories in this paper are suitable for threedimensional problems without any changes) in an isotropic medium, either homogeneous or inhomogeneous background medium, and its IEversion 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 finitedifference operators involving the FDCSI 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 finitedifference responses to improve the efficiency of the forward modeling. We invert these finitedifference 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 FDMRCSI 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 finitedifference 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 twodimensional forward modeling using finitedifference 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 wellknown 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 pwave incident excitation.
The parameters and wave fields are uniformly discretized on a regular twodimensional domain at intervals of and by using the secondorder finitedifferences 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 fivepoint 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 rowmajor rule, we obtain a linear equation system as where is the finitedifference 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 finitedifference matrix only once. However, a new matrix decomposition is required for another frequency.
3. FiniteDifference Contrast Source Inversion Algorithm
3.1. Problem Formulation
We consider the twodimensional 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 finitedifference 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. FiniteDifference 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 finitedifference 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 FDCSI method is to update contrast source vectors by the conjugategradient method with PolakRibiere 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 updatestep size and is the PolakRibiere 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 conjugategradient can be found in [30]. After updating the gradient , the updatesteps 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 FDCSI algorithm, Considering the nonlinearity and illposedness of the inverse problem, and the local convergence of the FDCSI 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 backpropagation 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 finitedifference multiplicative regularization contrast source inversion (FDMRCSI) method can be found in Appendix A of [8].
It is very similar to minimizing functional (29) using conjugategradient 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 conjugategradient method with PolakRibiere 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 timeconsuming 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 FDMRCSI is to calculate these operators and to estimate the conjugate and updatesteps . 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 backsubstitution needs operations for each source position. In the framework of FDMRCSI, backsubstitution is required twice for each of the two operators and per iteration. Hence, the computational complexity of FDMRCSI 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 FDMRCSI method as follows: in which represents the LU decomposition computational cost for forward modeling and denotes the backsubstitution computation for each source during each iteration process. Table 1 indicates the details of computational complexity for FDMRCSI and usual inversion methods.
 
is the iteration number, presents the number of source, and denotes the number of the grid cells. 
From Table 1, we observe that the FDMRCSI 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 FDMRCSI method is a greatly efficient method.
4. Numerical Results
In this section, we present some numerical examples based on the FDMRCSI 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 finitedifference timedomain 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 timedomain forward method are transformed into frequency domain using a fast Fourier transform (FFT) (see [31]). After that, we select some frequencies for the FDMRCSI 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/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.
(a) True model
(b) Reconstruction result
(c) Discrepancy
(d) Model slice
The inversion domain is divided into uniform grids using the secondorder 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 FDMRCSI 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 FDMRCSI 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.
(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 FDMRCSI 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 FDMRCSI method is convergent.
4.2. Reconstruction in an Inhomogeneous Background Medium
In this section, we examine the performance of the FDCSI 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 fourlayer distribution of a mass density of 1.8, 1.6, 1.5, and 1.8 g/cm^{3} 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 twodimensional 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 finitedifference timedomain (FDDT) 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 FDMRCSI 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 finitedifference 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 socalled sequential multifrequency inversion approach, that one usually employs this technique in the GaussNewton 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 FDMRCSI 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 FDMRCSI and FDMRCSI 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 finitedifference contrast source inversion method incorporating a multiplicative regularization factor is successfully employed in seismic fullwaveform 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 FDMRCSI 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.
References
 A. Tarantola, “A strategy for nonlinear elastic inversion of seismic reflection data,” Geophysics, vol. 51, no. 10, pp. 1893–1903, 1986. View at: Google Scholar
 P. Mora, “Nonlinear twodimensional elastic inversion of multi offset seismic data,” Geophysics, vol. 52, no. 9, pp. 1211–1228, 1987. View at: Publisher Site  Google Scholar
 R. G. Pratt and M. H. Worthington, “Inverse theory applied to multisource crosshole tomography. Part 1: acoustic waveequation method,” Geophysical Prospecting, vol. 38, no. 3, pp. 287–310, 1990. View at: Publisher Site  Google Scholar
 R. G. Pratt, “Inverse theory applied to multisource crosshole tomography. Part 2: elastic waveequation method,” Geophysical Prospecting, vol. 38, no. 3, pp. 311–329, 1990. View at: Publisher Site  Google Scholar
 B. Han, H. S. Fu, and Z. Li, “A homotopy method for the inversion of a twodimensional acoustic wave equation,” Inverse Problems in Science and Engineering, vol. 13, no. 4, pp. 411–431, 2005. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 H. BenHadjAli, S. Operto, and J. Virieux, “Velocity modelbuilding by 3D frequencydomain, fullwaveform inversion of wideaperture seismic data,” Geophysics, vol. 73, no. 5, pp. VE101–VE117, 2008. View at: Publisher Site  Google Scholar
 W. W. Symes, “Migration velocity analysis and waveform inversion,” Geophysical Prospecting, vol. 56, no. 6, pp. 765–790, 2008. View at: Publisher Site  Google Scholar
 A. Abubakar, W. Hu, T. M. Habashy, and P. M. van den Berg, “Application of the finitedifference contrastsource inversion algorithm to seismic fullwaveform data,” Geophysics, vol. 74, no. 6, pp. WCC47–WCC58, 2009. View at: Publisher Site  Google Scholar
 J. J. Cao, Y. F. Wang, J. T. Zhao, and C. Yang, “A review on restoration of seismic wavefields based on regularization and compressive sensing,” Inverse Problems in Science and Engineering, vol. 19, no. 5, pp. 679–704, 2011. View at: Publisher Site  Google Scholar  MathSciNet
 A. Fichtner, P. H. Bunge, and H. Igel, “The adjoint method in seismology. I. Theory,” Physics of the Earth and Planetary Interiors, vol. 157, no. 12, pp. 86–104, 2006. View at: Publisher Site  Google Scholar
 P. M. van den Berg and R. E. Kleinman, “A contrast source inversion method,” Inverse Problems, vol. 13, no. 6, pp. 1607–1620, 1997. View at: Publisher Site  Google Scholar  MathSciNet
 P. M. van den Berg, A. L. van Broekhoven, and A. Abubakar, “Extended contrast source inversion,” Inverse Problems, vol. 15, no. 5, pp. 1325–1344, 1999. View at: Publisher Site  Google Scholar  MathSciNet
 A. Abubakar, W. Hu, P. M. van den Berg, and T. M. Habashy, “A finitedifference contrast source inversion method,” Inverse Problems, vol. 24, no. 6, Article ID 065004, pp. 1–17, 2008. View at: Publisher Site  Google Scholar
 X. Chen, “Subspacebased optimization method for solving inversescattering problems,” IEEE Transactions on Geoscience and Remote Sensing, vol. 48, no. 1, pp. 42–49, 2010. View at: Publisher Site  Google Scholar
 X. Chen, “Subspacebased optimization method for inverse scattering problems with an inhomogeneous background medium,” Inverse Problems, vol. 26, Article ID 074007, 2010. View at: Publisher Site  Google Scholar  MathSciNet
 A. Abubakar, G. Pan, M. Li, L. Zhang, T. M. Habashy, and P. M. van den Berg, “Threedimensional seismic fullwaveform inversion using the finitedifference contrast source inversion method,” Geophysical Prospecting, vol. 59, no. 5, pp. 874–888, 2011. View at: Publisher Site  Google Scholar
 R. F. Bloemenkamp, A. Abubakar, and P. M. van den Berg, “Inversion of experimental multifrequency data using the contrast source inversion method,” Inverse Problems, vol. 17, no. 6, pp. 1611–1622, 2001. View at: Publisher Site  Google Scholar  MathSciNet
 A. Abubakar and P. M. van den Berg, “The contrast source inversion method for location and shape reconstructions,” Inverse Problems, vol. 18, no. 2, pp. 495–510, 2002. View at: Publisher Site  Google Scholar  MathSciNet
 L. P. Song, C. Yu, and Q. H. Liu, “Throughwall imaging (TWI) by radar: 2D tomographic results and analyses,” IEEE Transactions on Geoscience and Remote Sensing, vol. 43, no. 12, pp. 2793–2798, 2005. View at: Publisher Site  Google Scholar
 M. Li, A. Abubakar, and P. M. van den Berg, “Application of the multiplicative regularized contrast source inversion method on 3D experimental Fresnel data,” Inverse Problems, vol. 25, no. 2, Article ID 024006, 23 pages, 2009. View at: Publisher Site  Google Scholar  MathSciNet
 A. Zakaria, C. Gilmore, and J. LoVetri, “Finiteelement contrast source inversion method for microwave imaging,” Inverse Problems, vol. 26, no. 11, Article ID 115010, 21 pages, 2010. View at: Publisher Site  Google Scholar  MathSciNet
 G. Pelekanos, A. Abubakar, and P. M. van den Berg, “Contrast source inversion methods in elastodynamics,” Journal of the Acoustical Society of America, vol. 114, no. 5, pp. 2825–2834, 2003. View at: Publisher Site  Google Scholar
 D. Komatitsch and R. Martin, “An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation,” Geophysics, vol. 72, no. 5, pp. SM155–SM167, 2007. View at: Publisher Site  Google Scholar
 R. Martin, D. Komatitsch, and S. D. Gedney, “A variational formulation of a stabilized unsplit convolutional perfectly matched layer for the isotropic or anisotropic seismic wave equation,” Computer Modeling in Engineering & Sciences, vol. 37, no. 3, pp. 274–304, 2008. View at: Google Scholar  MathSciNet
 T. A. Davis, Direct Methods for Sparse Linear Systems, SIAM, Philadelphia, Pa, USA, 2006. View at: Publisher Site  MathSciNet
 P. M. van den Berg, A. Abubakar, and J. T. Fokkema, “Multiplicative regularization for contrast profile inversion,” Radio Science, vol. 38, no. 2, 2003. View at: Publisher Site  Google Scholar
 K. R. Kelly, R. W. Ward, S. Treitel, and R. M. Alford, “Synthetics eismograms: a finite difference approach,” Geophysics, vol. 41, no. 1, pp. 2–27, 1976. View at: Publisher Site  Google Scholar
 R. G. Pratt, “Frequencydomain elastic wave modeling by finite differences: a tool for crosshole seismic imaging,” Geophysics, vol. 55, no. 5, pp. 626–632, 1990. View at: Publisher Site  Google Scholar
 K. W. Brodlie and D. A. H. Jacobs, “Unconstrained minimization,” in State of the Art in Numerical Analysis, pp. 229–268, Academic Press, 1997. View at: Google Scholar
 J. Nocedal and S. J. Wright, Numerical Optimization, Springer, New York, NY, USA, 2006. View at: MathSciNet
 P. Duhamel and M. Vetterli, “Fast Fourier transforms: a tutorial review and a state of the art,” Signal Processing, vol. 19, no. 4, pp. 259–299, 1990. View at: Publisher Site  Google Scholar  MathSciNet
Copyright
Copyright © 2014 Bo Han et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.