Full-waveform Velocity Inversion Based on the Acoustic Wave Equation

Full-waveform velocity inversion based on the acoustic wave equation in the time domain is investigated in this paper. The inversion is the iterative minimization of the misfit between observed data and synthetic data obtained by a numerical solution of the wave equation. Two inversion algorithms in combination with the CG method and the BFGS method are described respectively. Numerical computations for two models including the benchmark Marmousi model with complex structure are implemented. The inversion results show that the BFGS-based algorithm behaves better in inversion than the CG-based algorithm does. Moreover, the good inversion result for Marmousi model with the BFGS-based algorithm suggests the quasi-Newton methods can provide an important tool for large-scale velocity inversion. More computations demonstrate the correctness and effectives of our inversion algorithms and code.


Introduction
The velocity of underground media is an important parameter for resources exploration.It may indicate the structure and position of reservoir.Full-waveform inversion uses travel time, amplitude and phase information simultaneously to inverse media velocity.It is a powerful way in reconstructing complex velocity structures.The inversion can be performed in the time-space domain [1][2][3][4] or in the frequency-space domain [5][6][7][8][9].The frequency-domain inversion approach is equivalent to the time-domain inversion approach if all of the frequency data components are used in the inversion process [9].
Full-waveform inversion can be formulated as a nonlinear optimization problem with an objective function consisting of the difference between observed and synthetic wave field in a suitable norm, usually L 2 norm.This choice results in least-squares formulation.However, solving full-waveform inversion problems using nonlinear optimization methods remains a challenge.This is mostly because of mathematical and numerical difficulties, which include multiple minima and ill-posedness.The objective function is relatively insensitive to wave number components of the model that are shorter than wavelength.The eigenvalues of the linearized forward operator corresponding to those components are nearly zero.Recently, the technique of multifrequency inversion of full-waveform has been investigated.For frequency domain inversion, the important issue is how to choose the various frequencies for inversion [10].Moreover, the wave modeling in the frequency domain requires solving a large-scale system at each iteration which is time-consuming for large-scale problem.Comparing with the frequency-domain inversion methods, the time-domain methods have the advantage of high efficiency in forward modeling.
The algorithms used to solve the optimization problem include linearized approximation, the steepest descent algorithm, the conjugate gradient (CG) [11] and the quasi-Newton methods [12][13][14][15].However, these methods behave differently as their own characters and the essential ill-posedness of the problem.In this paper, we implement full-waveform inversion with the CG method which is a typical gradient-based method, and with the BFGS method which is a representative of the quasi-Newton methods.As to the forward method, the finitedifference method is employed as its high efficiency.The discretization schemes including the absorbing boundary conditions and corner conditions are given.

Finite Difference Simulation
The wave equations can be solved in either the time or the frequency domain.The usual numerical methods such as the finite difference method, the finite volume method and the finite element method can all be employed.Here we use the finite difference method for its high efficiency.The 2-D acoustic wave equation in the time domain can be written as where is the pressure, is the velocity, ( , , ) x z is the source location.Usually the initial conditions are zero, i.e.
Let , be the pressure at the time and at the spatial position , and , n m be the velocity at , then the difference scheme of (1) is with the initial conditions 1 0 , , 0, 0 The numerical stability condition for (3) and ( 4) is min( , ) 2 max( ).
In the wave simulation, the unsuitable choice of time and spatial steps may cause serve dispersion and wave distortion.In order to suppress numerical dispersion, it is usually required that there are 10 sampling points in a wavelength for the second-order difference scheme, i.e., max min( ) min ( , ) . 10 Obviously for the same frequency the numerical dispersion increases as the spatial steps rise.However, the dispersion can be depressed if the high-order accurate computational schemes are used.
In numerical computations, the computational domain is usually truncated to a finite region.So absorbing boundary conditions (ABCs) are required to absorb boundary reflections of outgoing waves.A number of ABCs have been introduced in wave simulation, for example, the paraxial approximation method [16] and the perfectly matched layer (PML) method [17].The PML technique was proposed for electromagnetic fields equations, but can be applied to other wave equations and can fit irregular boundary.
Suppose the computational domain is a rectangular area We set the boundary condition on surface ( z 0  ) be the Dirichlet boundary condition as the data are gathered near the surface in our problem.For the boundary conditions of other three sides and four corners, we select the following Clayton and Engquist [16] conditions.For the top boundary: For the bottom boundary: For the left boundary: For the right boundary: The absorbing conditions for four corners can be derived through the 45°rotation of the ABCs above.The results are the following: For the left-upper corner: For the right-upper corner: For the left-down corner: For the right-down corner: The difference schemes for the absorbing boundary condition ( 7)- (10) can also be obtained by the secondorder central finite difference method, i.e., 1 , 0, ( ) 0 , 4 ( ) 0 , 4 ( ) 0 , 4 where   and   are the forward and backward difference ope ators with first-order accuracy respectively, r 0  is th ral rence operator with secondcuracy, e cent diffe order acx  and z  are the spatial steps in x and directions respectively.Similarly, the difference schemes rs z for four corne conditions ( 12)-( 15) can be obtained in the following: For the left-upper corner: ( , ) (0, 1) .
upper corner: For the right- .
For the left-down corner: .
For the right-down corner: .
Full-waveform inversion is minimization of the following objective function (24)

Two Full-waveform I
n form the where obs u is the observed wave field and cal u is the synthetic data.The functional (25) Different choices of k  lead to different conjugate gradient methods.For example, the FR method [11], the HS method [18] and the PRP method [19,20]: where partial derivative . Here we adopt the PRP m The of with resp known velocity v is ethod.ect to un- where J et derivative matrix with dimension is the Fréch ( ) The gradient g is a vector with dimension M , which can as be written For the line search method we may use the Ai gorithm [21] or the strong Wolfe algorithm [22,23].The A , mijo alimijo algorithm requires 1 ) ( ) ( ) where 1 (0,1) b  and we select 1 0.001 b in computations.And the strong Wolfe algorithm requires the following condition besides condition (33) above 3.2 : Update the velocity model: , The BFGS method modifies matrix and its n by verse is give B  with (36).3.9: Use the Cholesky method to solve the equation to get .
, continue with step 3.In step 3.7 above, the scaling factor k  is introduced which can improve convergence [24].CG method is on nve he BFGS method is supe each iteration.

merica
ary conditions are used.Figure 2 is the 1 but the ABCs are ndary reflections in The ly linearly co rgent while t rlinear convergent but requires solving a system at

Nu l Computations
First of all, we test the validity of the finite difference method for wave modeling.relt after 20 iterations and Figure 6 is that after 200 iterations.For the BFGS-based algorithm the inversion results are shown in Figures 7-9, which are the results after 20 iterations, 100 iterations and 200 iterations, respectively.Comparing Figures 7-9 with Figures 5 and 6, we can see that the BFGS-based algorithm has higher inverse accuracy than the CG-based algorithm.In order to test the ability of the BFGS-based algorithm further, we inverse the Marmousi model next.
The Marmousi model is a benchmark model for testing the ability of migration or inversion methods [25].The velocity model is shown in Figure 10.We can see it has many complex structures.The synthetic data is also generated by the finite difference code.The following computations parameters are adopted: The inversion results with the CG-based algorithm are shown in Figures 5 and 6. Figure 5 is the inversion su

1 : 2 :
Give the initial velocity model and the iterative 2 0.9 b  full-wav computations.N ive the following form inversion algorithm combining with the CG method.Compute the objective function 0 f and the gradient g.If the gradient satisfies the stopping condition || || 0 k  ber g  stop  then stop the iteration.Otherwise set p 0  The algorithm above doesn't use the information of second-order derivate of objective function.The BFGS m

2 : 1 : 3 . 2 :
Compute the objective function 0 f and the gradient g.If the gradient satisfies the stopping rule || || itera Use the strong Wolfe line search algorithm to obtain the step k  ; Update the velocity model:

Figure 1
is the wave filed received near the surface for a fixed source.The boundary reflections are very serious in Figure 1 since only Dirichlet bound wave field corresponding to Figure used.We can see the serious bou Figure 1 are eliminated obviously in Figure 2.

Figure 1 .Figure 2 .
Figure 1.Wavefield received near the surface by the finite difference method with Dirichlet boundary conditions.

Figure 3 .
Figure 3. True model for testing two inversion algorithms.

Figure 4 .
Figure 4.The configuration of sources and receivers for full-waveform inversion.

Figure 11 .
Figure 11.Initial velocity model for inversion.6 x z m     , 0.6 t ms   , 3501 N  .t is initial model and Figure 12 is the inversion resu

Figure 11
Figure 11 lt fter 200 iterations.This result clearly shows that BFGSbased algorithm has obtained very good inversion accuracy.t a he