Abstract
In this paper, the acoustic full waveform inversion based on the wavelet method is proposed. It allows to determine density and velocity parameters simultaneously in the time domain. The forward problem is solved based on the wavelet method. The numerical schemes for modelling are constructed in detail. The stability condition is derived for the first time. The inversion is an optimization process for minimizing the mismatch between the synthetic data and the observed data. The computational scheme for the gradient of the objective function is constructed based on the matrix analysis technique. The inversion is implemented from low frequency to high frequency successively. Furthermore, the inversion result of current low frequency is served as the initial guess of next high frequency. This hierarchic strategy includes the smaller wavelengths progressively and improves the inversion robustness. Numerical computations for the benchmark Marmousi model demonstrate that the method has good ability to reconstruct media parameters in complex structures.
1. Introduction
Seismic full waveform inversion (FWI) can quantitatively determine model parameters of the media by using all information of wavefield including phase, amplitude and traveltime. It is a feasible and essential technology to build a good subsurface model to explore oil and gas reservoir in geophysics. FWI was initially proposed by Lailly and Tarantola in 1980s for the acoustic wave equation in the time domain [Citation1,Citation2]. It is in fact a data-fitting procedure based on numerical optimization. In [Citation1,Citation2], the gradient method is applied and the gradient direction is constructed by crosscorrelating forward wavefield from the sources with backward wavefield from the residuals. This gradient computational way requires only several forward computations and avoids intensive computational amount which the explicit approaches meet. This key idea was quickly applied to elastic media in the time domain using multicomponent data [Citation3,Citation4]. In 1990, the formulation of FWI in the frequency domain was developed by Pratt [Citation5,Citation6]. Each type of method has its advantages. For the time-domain method there is no need to transform data into the Fourier domain, whereas the frequency-domain needs. So the wavefield in the time domain is real valued whereas it is complex valued in the frequency domain. In addition, the parallel code implementation of FWI in the time domain is relatively easy. An overview of FWI in exploration geophysics may refer to [Citation7]. FWI also has been developed in other domains such as the Laplace domain [Citation8] and the Laplace-Fourier domain [Citation9]. Recently, FWI has been explored in viso media [Citation10] and anisotropic media [Citation11]. Due to the the significant advance of computer power in last ten years, FWI has been used both in academics and geophysical application.
FWI is a non-linear and ill-posed inverse problem and the solution is non-unique. Instead of the global optimization methods such as the simulated annealing algorithm and the genetic algorithms, the gradient-based minimization techniques are usually adopted (e.g. see [Citation7,Citation12,Citation13]). They include the steepest descent method, the Gauss-Newton and full Newton methods [Citation14], the Newton-CG method [Citation15], the truncated Newton method [Citation10,Citation12] and so on. The Newton method is quadratically convergent but requires the computation of the second-order derivatives of the objective function (Hessian matrix) with respect to the model parameters, which is not adopted generally for its huge computational cost. The main problem in FWI is that the inversion usually converges to a local minimum solution when the initial model deviates significantly from the true model. Besides the nonlinearity and non-convexity of the misfit function, the cycle skipping is another major reason to cause the phenomenon of multiple solutions. In particular, each local minimum of the misfit function corresponding to fitting the data up to one phase or several shifts. Several types of techniques have been proposed to overcome this problem. The first is the regularization strategy such as the Tikhonov regularization [Citation16,Citation17] and the total variation (TV) regularization [Citation11,Citation18], which have been applied to stabilize the ill-posedness in the nonlinear inversion process. The second is the multiscale method promoted by Bunk [Citation19], which overcomes the ill-posedness of FWI by incorporating progressively shorter wavelengths of the model. In this method, FWI is performed by starting low frequencies and moving to higher frequencies successively. This strategy includes the smaller wavelengths progressively and guarantees the inversion convergence utmost since the objective function has few local minima at low frequencies. It has been realized it is more effective than the Tikhonov regularization and TV regularization. Recently, the optimal transport method, which increases the convexity of the misfit function by the Wasserstein metric rather than the metric, has been applied to FWI [Citation20,Citation21]. The multiscale method can be considered as a macro-regularization strategy.
In FWI, the forward problem is required to solve repeatedly. There are several numerical methods to solve the wave equation, for example, the finite difference method [Citation22,Citation23], the pseudo-spectral method [Citation24,Citation25], the finite volume method [Citation26], the finite element method [Citation27–29]. Each method has its own advantages and disadvantages. Among of them, the finite difference method is very popular in FWI since its high efficiency and relative easiness of code implementing. In comparison with the traditional numerical methods, the wavelet-based method for solving partial differential equations including wave equations numerically is relatively new [Citation30–34]. In wavelet methods, the wavelet basis has the good potential to capture the complex structures since its multiresolution. In this paper, we investigate and develop the wavelet-based FWI for determining the velocity and density coefficients in the acoustic wave equation.
In this paper, the time-domain acoustic FWI based on the wavelet method is investigated. The rest of the paper is organized as follows. In Section 2, the forward method by the wavelet method is investigated in detail. In Section 3, the inverse method for the density and velocity parameters in the acoustic wave equation is described. In Section 4, numerical computations are presented to illustrate the performance of the method. Finally, the conclusion is given in Section 5.
2. Forward method
The two dimensional acoustic wave equation with density and velocity can be written as (1) (1) where u is the pressure, ρ is the density, v is the velocity, is the source functions at the point position . Introducing the auxiliary variables and , we can rewrite (Equation1(1) (1) ) as (2) (2) (3) (3) (4) (4) The initial conditions are (5) (5) Equation (Equation1(1) (1) ) or the equivalent (Equation2(2) (2) )–(Equation4(4) (4) ) with (Equation5(5) (5) ) and suitable boundary conditions form an initial boundary value problem, which the uniqueness and existence of the solution to this problem is known theoretically [Citation35]. This forward problem is also well posed. Numerically, we have several methods to solve it. For example, the finite difference method [Citation22,Citation23], the finite element method [Citation27,Citation28] and so on. If the staggered-grid finite difference method is applied to (Equation2(2) (2) )–(Equation4(4) (4) ) and (Equation5(5) (5) ), we get the following schemes (6) (6) (7) (7) (8) (8) with initial conditions (9) (9) where is the time step, and and are the spatial steps in the x and z directions respectively. In the following, we apply the wavelet method to solve it. We first introduce the orthonormal wavelet bases and some related properties. Then we derive the wavelet-based computational schemes for solving (Equation2(2) (2) )–(Equation4(4) (4) ).
2.1. Orthonormal wavelet basis
In this subsection, we briefly review the orthonormal bases of compactly supported wavelets. For the details we refer to [Citation36–39]. The orthonormal basis of compactly supported wavelets in is constructed by the dilation and translation of a single function , (10) (10) where . The function has a companion, the scaling function . The functions and satisfy the following relations (11) (11) (12) (12) where (13) (13) and (14) (14) In addition, the function has M vanishing moments (15) (15) Where L in (Equation11(11) (11) )–(Equation12(12) (12) ) is related to M. For Daubechies' wavelets [Citation37] applied in this paper, L = 2M. The coefficient sets and in (Equation11(11) (11) )–(Equation12(12) (12) ) are quadrature mirror filters. Once the filter H has been chosen, the wavelet function ψ and the scaling function φ are determined. The wavelet coefficients of Daubechies' wavelets can be computed by multiresolution analysis algorithm [Citation39].
In the multiresolution analysis, the Hilbert space is decomposed into a chain of closed subspaces: such that Let be the orthogonal complement of in : then can be expressed as a direct sum of , i.e. . For a fixed scale j, the scale functions and the wavelets form the orthonormal bases of and respectively, i.e. Thus , there is the following decomposition (16) (16) where (17) (17) with where is defined as for any functions on I, and An effective way to construct two-dimensional wavelets is by using the tensor product of one-dimensional wavelets. Define . Here denotes the tensor product. Then similarly, is decomposed into a chain of closed subspaces . Making the decomposition of the direct sum, we konw where Thus . Furthermore, form the orthonomal bases of ; form the orthonomal bases of .
2.2. Wavelet-based forward scheme
In this subsection, we construct the forward computational scheme on staggered grids based on the wavelet method. We employ the scaling function φ to express the wavefields u, and in (Equation2(2) (2) )–(Equation4(4) (4) ). Let (18) (18) (19) (19) (20) (20) Applying φ as the test function to the source-free system (Equation2(2) (2) )–(Equation4(4) (4) ), we have (21) (21) (22) (22) (23) (23) We need to compute the integrals in (Equation21(21) (21) )–(Equation23(23) (23) ). Let (24) (24) (25) (25) (26) (26) (27) (27) Obviously, (28) (28)
So by the property in multiresolution analysis, we have (29) (29) (30) (30) (31) (31) Since , then we have (32) (32) (33) (33) (34) (34) By the properties of translation and symmetry, we have (35) (35) (36) (36) (37) (37) (38) (38) Therefore, we only need to compute and . By the orthogonality of φ, we have (39) (39) In this paper, we employ Daubechies wavelet [Citation37] with M = 4 since this wavelet is smooth. The coefficients can be computed numerically. The details are given in Appendix 1. The results are the following: (40) (40) In numerical computations the computational domain is finite and absorbing boundary condition (ABC) is necessary to eliminate the boundary reflections. Several approaches can be used, for example, the paraxial approximation [Citation40], the exact nonreflecting conditions [Citation41,Citation42] and the perfectly matched layer (PML) method [Citation43]. The PML is the most widely used method and we adopt it in this paper. As my best knowledge, there is no work to incorporate the PML to solve the wave equation based on the wavelet method. The idea of PML is to design a particular layer around the computational domain in which the waves will be attenuated satisfactory. The details are omitted for saving space.
2.3. Stability analysis
We use the Fourier method to analyse the stability. Suppose the computational domain is . Note that the wavelet computations map the computational domain to . So the first derivative will produce a factor . For convenience, denote the wavelet coefficients , , for a fixed k by , , , , respectively. We rewrite (Equation21(21) (21) )–(Equation23(23) (23) ) as the following recursive form in the time index n (41) (41) (42) (42) (43) (43) By applying the Fourier stability analysis method [Citation44], we have (44) (44) where , and are the the Fourier transform of , and u respectively, and (45) (45) Here . Note that Eliminating and in (Equation44(44) (44) ), we have (46) (46) Let and . The characteristic equation of (Equation46(46) (46) ) is The stability requires the two characteristic roots are in the unit circle, which yields the condition (47) (47) Solving (Equation47(47) (47) ) we get the sufficient and necessary stability condition for the wavelet schemes (Equation41(41) (41) )–(Equation43(43) (43) ) (48) (48) where is the maximal velocity.
3. Inversion method
To develop the FWI method, we firstly define the objective function as the norm of the residuals between the synthetic data and the observed data : (49) (49) where is the shot position and the receiver position, is the model parameters, and T is the total recording time. The summations in (Equation49(49) (49) ) are performed over the number of all sources and receivers. The FWI is an iterative process for minimizing the objective function. The iterative formulation can be written as (50) (50) where the subscript k denotes the kth iteration, is the initial model, indicates the model update, and is the stepsize and a scalar estimated via a line search method. In this paper, the strong Wolfe line search [Citation45] is used. Generally, the model update has the following form (51) (51) where is a symmetric and nonsingular matrix. In the steepest descent method is simply the identity matrix, whereas in Newton method is the exact Hessian . In quasi-Newton methods, is an approximation to the Hessian that is updated at every iteration. Methods that use the Newton direction have a fast rate of local convergence. The main drawback of the Newton direction is the need for the Hessian. For large-scale FWI, the computation cost of the Hessian and its inverse is beyond current computational capability. So the Newton method is seldom applied in FWI. The quasi-Newton methods do not require computation of the Hessian and yet still attain a superlinear rate of convergence [Citation45]. In this paper, a quasi-Newton method named limited-memory BFGS (L-BFGS) method [Citation45] is applied for its good behaviour in the time-domain FWI. Recently, the truncated Newton method has been applied to FWI (e.g. see [Citation10,Citation12]). The truncated Newton method was initially developed by Dembo and Steihaug [Citation46] in 1983. It solves the Newton equation by the conjugate gradient method iteratively. Thus it has much computational cost to obtain Hessian than the L-BFGS method. The L-BFGS method consists in approximating the inverse Hessian by several previous values of the gradient.
We remark that all optimization-based approaches mentioned above are converged locally and suffer from the local-minima trapping problem. So some techniques such as the successive inversion strategy is necessary to improve the robustness and accuracy of the inversion. The inversion is implemented from low frequency to high frequency successively. The inversion at the stage of using low frequency date has more probability to escape the local-minima trap. We note that a so-called convexification method initially proposed by Klibnanov [Citation47] is a globally convergent method which can overcome the local-minima trapping problem. The application of this method is beyond the scope of this paper. For the details, we refer to the series of works of Klibanov with coauthors [Citation48,Citation49] and the references therein.
As FWI is an ill-posed problem, a regularization term is necessary in the objective function to stable the solution. So (Equation49(49) (49) ) is revised as (52) (52) where α is the regularization parameter which is estimated based on the amplitude of the objective function in this paper, and is the regularization function with two typical approaches. One is the Tikhonov regularization (53) (53) The other is the TV regularization (54) (54) Generally, the TV regularization improves inversion accuracy at discontinuous positions of model parameters whereas the Tikhonov regularization smoothes the inversion solutions.
Note that the gradient of the objective function is required in FWI. The gradient of with respect to model parameters can be computed by the standard finite difference method directly. For example, for the Tikhonov regularization (Equation53(53) (53) ), the result is (55) (55) As for the TV regularization (Equation54(54) (54) ), we have (56) (56) where is a small parameter to prevent the singularity. In the following, we will focus on the gradient of or with respect to the model parameters v and ρ, which is required in the L-BFGS method.
3.1. Computational schemes for the gradient
In this subsection, we will derive the numerical schemes for the gradient of the objective function with respect to the model parameters v and ρ. Let the discrete points number in x and z directions are and respectively. Then the difference scheme for (Equation1(1) (1) ) with the second-order accuracy is (57) (57) where Arranging the spatial indexes in x direction first and then in z direction, we rewrite (Equation57(57) (57) ) as the following matrix form: (58) (58) where (59) (59) (60) (60) with (61) (61) The nonzero elements of in (Equation61(61) (61) ) are the following Note that are the diagonal matrices and . We know is symmetric since . Thus .
Furthermore, the system (Equation58(58) (58) ) can be rewritten as the following matrix form (62) (62) where (63) (63) Note that A and K are all symmetric matrices.
Now making derivatives of the objection with respect to and , we get (64) (64) where and And making derivatives of (Equation62(62) (62) ) with respect to and , we have (65) (65) Inserting (Equation65(65) (65) ) into (Equation64(64) (64) ) yields (66) (66) By (Equation63(63) (63) ) we have (67) (67) where the elements are given in Appendix 2. Then for a given vector : (68) (68) we have (69) (69) We now consider to compute . First by (Equation63(63) (63) ), we have (70) (70) where the elements and are also given in Appendix 2. Thus (71) (71) where , and () Since A and K are symmetric matrices, we have the following lemma.
Lemma 3.1
Suppose is solution of then is the solution of .
Proof.
By the definition we rewrite for as (72) (72) Similarly, we rewrite for as (73) (73) By (Equation72(72) (72) ) and (Equation73(73) (73) ) we know the conclusion of this lemma is true. The proof is complete.
From Lemma 3.1, we know that is the numerical solution to the backward problem corresponding to (Equation1(1) (1) ) with the source term replaced by the residuals.
Therefore, by (Equation66(66) (66) ) we obtain the computational schemes for the gradient of the objection: (74) (74)
and (75) (75) where
4. Numerical computations
In this section we first consider the numerical computations of forward problem. Then we consider the inversion for a complicated benchmark model named Marmousi model.
4.1. Forward computations
We now present two examples to show the numerical results of the forward problem. The first example is a homogeneous triangular model with constant velocity and constant density . The physical domain is . In computations, the spatial step is and the time step is . The source is located at , which is the Ricker wavelet with the time history (76) (76) where is the central frequency. Figure is the snapshot of wavefield at propagation time . Figure (a) is the result by the staggered-grid method without absorbing boundary conditions. In Figure (a), we can see the serious reflections from the top and bottom boundaries. Figure (b) is the result by the wavelet method with the PML boundary conditions. The PML is around the computational model. In Figure (b) we can see the boundary reflections are eliminated obviously.
The second example is the Marmousi model [Citation50] with complex structures. The exact velocity is shown in Figure and it will be used in the inversion computations. The computational domain is . The spatial step is and the time step is . The thickness of PML is . The shots and receivers are arranged under the surface . The sketch map of the data acquisition system is shown in Figure . This data observation geometry is analogous to the acquisition system of reflection seismology in geophysical exploration. In Figure , the blue points denote the shots and the red crosses denote the receivers. Figure (a) is a shot gather with the shot located at the left side of the model. Figure (b) is another shot gather with the shot located at the middle of the model. In both Figure (a,b), we can see the boundary reflections are eliminated obviously. For comparison with Figure , the corresponding shot gather by the staggered-grid method is shown in Figure . Comparing Figure with Figure , we can clearly see that the waveform in Figure by the wavelet method has much less dispersion.
4.2. Inversion computations
We now consider the FWI for the Marmousi model. This model is a benchmark model which is usually applied to test the ability of an inversion method. The exact model of velocity and density for the Marmousi model is shown in Figure . The density shown in Figure (b) varies from to . The density structure is similar with the velocity because high velocity corresponds to large density physically. The data acquisition system is shown in Figure . The number of sampling points in space is and with the step . The number of time points is with the step . There are 80 shots and each shot 40 receivers totally. The space of receivers is and the space of shots is . The initial model for FWI is shown in Figure , which is built by a linear increasing function from the surface to the bottom. Figure (a) is the initial velocity model and Figure (b) is the initial density model. As we can see, the initial model is quite different from the exact model and we can not see the useful information of complex structures such as the faults in the exact model Figure . The frequency multiscale strategy is applied in FWI. Particularly, we implement FWI from low frequency to high frequency step by step. And the inversion result at current frequency band is served as the initial model for the next frequency band. In order to obtain the data in a certain frequency band, we apply the following Blackman-Harris window to filter the data: where is the length of the filter and In computations, we set four different frequency bands. Moreover, the data in the current frequency band is contained the data in the next frequency band so the FWI is more robust. Figure is the inversion result of velocity and density with the data from frequency 0 to 5 Hz. Figure is the inversion result with the data from frequency 0 to 15 Hz. Figure is the inversion result with the 0–25 Hz data. The iteration number is fifty. From Figures –, we can see the inversion accuracy is improved gradually. Figure is the final inversion result with the whole frequency band data. Comparing Figure with Figure , we can see most of the structures of the exact model are inverted accurately.
FWI is a typical large-scale computational problem. We implement the parallel computations based on the message passing interface (MPI). The physical domain is separated into different subdomains. And each subdomain is assigned to a different processer. The data exchanges only occur at the boundary of subdomains. In the L-BFGS algorithm, the computations for searching direction only involve vector inner product in a matrix-free fashion and the computations do not repeat in the overlapping area. The number of discretized parameters in inversion is 245,514. The CPU time for the whole FWI is about 14 h with 128 processes. The processor is the 3.6 GHz Intel. The stoping criteria for the inversion is when the is less than a small quantity while the iteration times are not beyond a fixed number.
We remark that the multi-parameter FWI including the density is harder than the single parameter FWI because of the coupling influence between the density and velocity. The single parameter inversion has better accuracy than the multi-parameter inversion. Figure shows the results of the single parameter inversion for velocity only by the wavelet method with the data in four different frequency bands, i.e. 0–2.5 Hz, 0–5 Hz, 0–15 Hz and all frequencies. From Figure (a,d), we can clearly see that the inversion accuracy improves gradually. Comparing Figure (a) with Figure (a), we also see that the final single parameter inversion result is better than the multi-parameter inversion, which is as we expected. In the computations, the observed data has been added the Gaussian noises.
5. Conclusions
The FWI is an attractive tool to obtain high resolution subsurface parameters in complex geological structures. We have investigated the wavelet-based FWI for determining the density and velocity coefficients in the acoustic wave equation. The numerical methods and schemes are established. The method solves the forward problem by applying Daubechies-4 wavelet to approximate the wavefield on staggered grids. The discrete schemes for forward modelling and the corresponding stability condition are derived in detail. The inverse problem is solved iteratively by minimizing the waveform misfit between the synthetic data and the observed wafefield. The discrete schemes for the gradient are derived based on the matrix technique. The frequency multiscale method implements FWI from low frequency to high frequency successively. Moreover, the inversion result of current frequency band is served as the initial guesses of the next frequency band, which improves robustness and accuracy of FWI. The multi-parameter FWI with density has lower inversion accuracy than the single parameter inversion. In FWI computations, the initial model is far from the exact model, which shows the method is less dependent on the initial model. Numerical computations for a benchmark model show the effectiveness and robustness of the method in this paper. The method can be extended to elastic model FWI.
Acknowledgements
I appreciate J. Jiang for his valuable helps in this research. Thanks also give to the reviewers and the editor for their helps.
Disclosure statement
No potential conflict of interest was reported by the author(s).
Additional information
Funding
References
- Lailly P. The seismic inverse problem as a sequence of before stack migrations. Conference on Inverse Scattering, Theory and Applications. Philadelphia (PA): SIAM; 1983. p. 206–220.
- Tarantola A. Inversion of seismic reflection data in the acoustic approximation. Geophysics. 1984;49(8):1259–1266.
- Mora P. Nonlinear two-dimensional elastic inversion of multioffset seismic data. Geophysics. 1987;52(9):1211–1228.
- Tarantola A. A strategy for nonlinear elastic inversion of seismic reflection data. Geophysics. 1986;51(10):1893–1903.
- Pratt RG, Worthington MH. Inverse theory applied to multi-source cross-hole tomography, part 1: acoustic wave equation method. Geophys Prospect. 1990;38(3):287–310.
- Pratt RG. Inverse theory applied to multi-source cross-hole tomography, part 2: elastic wave equation method. Geophys Prospect. 1990;38(3):311–329.
- Virieux J, Operto S. An overview of full-waveform inversion in exploration geophysics. Geophysics. 2009;74(6):WCC1–WCC26.
- Shin C, Cha YH. Waveform inversion in the Laplace-Fourier domain. Geophys J Int. 2009;177(3):1067–1079.
- Shin C, Cha YH. Waveform inversion in the Laplace domain. Geophys J Int. 2008;173(3):922–931.
- Yang P, Brossier R, Métivier L, et al. A time-domain preconditioned truncated Newton approach to viso-acoostic multiparameter full waveform inversion. SIAM J Sci Comput. 2018;40(4):B1101–B1130.
- Aghamiry HS, Gholami A, Operto S. ADMM-based multi-parameter wavefield reconstruction inversion in VTI acoustic media with TV regularization. Geophys J Int. 2019;219(2):1316–1333.
- Métivier L, Brossier R, Operto S, et al. Full waveform inversion and the truncated Newton method. SIAM Rev. 2017;59(1):153–195.
- Plessix RE. A review of the adjoint-state method for computing the gradient of a functional with geophysical applications. Geophys J Int. 2006;167(2):495–503.
- Pratt RG, Shin C, Hicks GJ. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion. Geophys J Internat. 1998;133(2):341–362.
- Epanomeritakis I, Akelik V, Ghattas O, et al. A Newton-CG method for large-scale three-dimensional elastic full waveform seismic inversion. Inverse Problems. 2008;24(3):1–26.
- Engl H, Hanke M, Neubauer A. Regularization of inverse problems. Dordrecht: Kluwer Academic Publishers Group; 1996.
- Guitton A. Blocky regularization schemes for full-waveform inversion. Geophys Prospect. 2012;60(5):870–884.
- Lin Y, Huang L. Acoustic- and elastic-waveform inversion using a modified total-variation regularization scheme. Geophys J Int. 2015;200(1):489–502.
- Bunks C, Saleck F, Zaleski S, et al. Multiscale seismic waveform inversion. Geophysics. 1995;60(5):1457–1473.
- Engquist E, Froese BD. Application of the Wasserstein metric to seismic signals. Commun Math Sci. 2014;12(5):979–988.
- Métivier L, Brossier R, Mérigot Q, et al. Measuring the misfit between seismograms using an optimal transport distance: application to full waveform inversion. Geophys J Int. 2016;205(1):345–377.
- Virieux J. P-SV wave propagation in heterogeneous media: velocity-stress finite-difference method. Geophysics. 1986;51(4):889–901.
- Zhang W, Tong L, Chung ET. A new high accuracy locally one-dimensional scheme for the wave equation. J Comput Appl Math. 2011;236(6):1343–1353.
- Fornberg B. A practical guide to pseudospectral methods. Cambridge: Cambridge University Press; 1996.
- Zhang W. Stability conditions for wave simulation in 3-D anisotropic media with the pseudospectral method. Commun Comput Phys. 2012;12(3):703–720.
- Zhang W, Zhuang Y, Zhang L. A new high-order finite volume method for 3D elastic wave simulation on unstructured meshes. J Comput Phys. 2017;340:534–555.
- Cohen G. Higher-order numerical methods for transient wave equations. Berlin: Springer-Verlag; 2002.
- Zhang W, Chung ET, Wang C. Stability for imposing absorbing boundary conditions in the finite element simulation of acoustic wave propagation. J Comput Math. 2014;32(1):1–20.
- Zhang W. Elastic full waveform inversion on unstructured meshes by the finite element method. Phys Scr. 2019;94: 115002.
- Beylkin G, Keiser JM. On the adaptive numerical solution of nonlinear partial differential equations in wavelet bases. J Comput Phys. 1997;132(2):233–259.
- Dahmen W. Wavelet and multiscale methods for operator equations. Acta Numer. 1997;6:55–228.
- Fröhlich J, Schneider K. An adaptive wavelet-vaguelette algorithm for the solution of PDEs. J Comput Phys. 1997;130(2):174–190.
- Hong TK, Kennett BLN. On a wavelet-based method for the numerical simulation of wave propagation. J Comput Phys. 2002;183(2):577–622.
- Qian S, Weiss J. Wavelets and the numerical solution of partial differential equations. J Comput Phys. 1993;106(1):155–175.
- Lions LJ, Magenes E. Non-homogeneous boundary value problems and application. Vol. 1. New York (NY): Springer-Verlag; 1972.
- Chui C. Wavelets: a tutorial in theory and applications. San Diego (CA): Academic Press; 1992.
- Daubechies I. Orthonormal bases of compactly supported wavelets. Commun Pure Appl Math. 1988;41:901–906.
- Daubechies I. Ten lectures on wavelet. Philadelphia (PA): Society for Industrial and Applied Mathematics; 1992.
- Mallat S. Multiresolution approximations and wavelet orthogonal bases of L2(R). Trans Amer Math Soc. 1989;315:69–88.
- Clayton R, Engquist B. Absorbing boundary conditions for acoustic and elastic wave equations. Bull Seismol Soc Am. 1977;67(6):1529–1540.
- Grote MJ, Keller JB. Exact nonreflecting boundary condition for elastic waves. SIAM J Appl Math. 2000;60(3):803–819.
- Zhang W, Tong L, Chung ET. Exact nonreflecting boundary conditions for three dimensional poroelastic wave equations. Commun Math Sci. 2014;12(1):61–98.
- Berenger JP. A perfectly matched layer for the absorption of electromagnetic waves. J Comput Phys. 1994;114(2):185–200.
- Thomas JW. Numerical partial differential equations: finite difference methods. New York (NY): Springer Science Business Media; 1995.
- Nocedal J, Wright ST. Numerical optimization. New York (NY): Springer Science Business Media; 1990.
- Dembo RS, Steihaug T. Truncated Newton algorithms for large-scale unconstrained optimization. Math Program. 1983;26(2):190–212.
- Klibanov MV. Global convexity in a three-dimensional inverse acoustic problem. SIAM J Math Anal. 1997;28(6):1371–1388.
- Khoa VA, Klibanov MV, Nguyen LH. Convexification for a three-dimensional inverse scattering problem with the moving point source. SIAM J Imaging Sci. 2020;13(2):871–904.
- Klibanov MV, Kolesov AE, Nguyen LH. Convexification method for an inverse scattering problem and its performance for experimental backscatter data fir buried targets. SIAM J Imaging Sci. 2019;12(1):576–603.
- Martin GS, Wiley R, Marfurt KJ. Marmousi2: an elastic upgrade for Marmousi. Leading Edge. 2006;25:156–166.
- Beylkin G. On the representation of operators in bases of compactly supported wavelet. SIAM J Numer Anal. 1992;29(6):1716–1740.
Appendices
Appendix 1. The computations for coefficients in (40).
By the definition of (Equation26(26) (26) ), we know (A1) (A1) For brevity, we set . Then we can rewrite (EquationA1(A1) (A1) ) as (A2) (A2) and by (Equation11(11) (11) ) we have (A3) (A3) Substituting m = k−n, we rewrite (EquationA3(A3) (A3) ) as (A4) (A4) Using the fact , we obtain (A5) (A5) where are the autocorrelation coefficients of : (A6) (A6) In fact, for Daubechies' wavelets with M vanishing moments, can be expressed as [Citation51] where Note that coefficients with even indices are zero. We remark on integer grids can be computed by solving a linear system. For example, for the Daubechies wavelet with M = 4, we have So by (EquationA5(A5) (A5) ) and setting we can obtain the results on staggered grids in (Equation40(40) (40) ).
Appendix 2. The elements of matrix in (67) and the elements of matrix in (70).
The elements of matrix are . By the expression of A in (Equation59(59) (59) ), we have (A7) (A7) Note that .
The elements of matrix in (Equation70(70) (70) ) are and . By the expression of A in (Equation59(59) (59) ), we have (A8) (A8) And by the expression of K in (Equation60(60) (60) ), we have (A9) (A9) where