Acoustic Based Crosshole Full Waveform Slowness Inversion in the Time Domain

We develop a new full waveform inversion (FWI) method for slowness with the crosshole data based on the acoustic wave equation in the time domain. The method combines the total variation (TV) regularization with the constrained optimization together which can inverse the slowness effectively. One advantage of slowness inversion is that there is no further approximation in the gradient derivation. Moreover, a new algorithm named the skip method for solving the constrained optimization problem is proposed. The TV regularization has good ability to inverse slowness at its discontinuities while the constrained optimization can keep the inversion converging in the right direction. Numerical computations both for noise free data and noisy data show the robustness and effectiveness of our method and good inversion results are yielded.


Introduction
Seismic exploration is one of the ways of identifying media properties and structures by propagation of waves.The wave is ignited by sources on the surface and propagates into underground.Due to different properties of media, various waves such as reflection wave and refraction wave are produced.The nature of wave reflection and refraction gives the physical information of media.
So the inverse of seismic waves can give geological properties of underground tion methods are impractical.Instead, the gradient-based algorithms are applied to find the local minimum of the misfit function.In 1986, Tarantola [3] extended the theory to elastic wave inversion.The idea was applied to two dimensional crosshole data inversion successfully in the frequency domain [6] [7] [8] [9].Later much more attentions are paid to the frequency domain method, and the frequency domain method has the same inversion ability like time domain method.Its advantage is that it can be implemented with one selected frequency.
However, it needs to solve a large-scale system unlike the time domain method which the forward problem can be implemented effectively in the time domain.
The FWI is a typical ill-posed problem for its nonlinearity and cycle-skipping problem.To overcome the ill-posedness, various methods have been developed, for example the multiscale method [10] [11] [12] and the regularization method [13] [14].Other methods such as the Laplace-domain method [15] [16] [17] and the envelop FWI method [18] are proposed to improve the inversion robustness.
Both the Laplace-domain FWI and envelop FWI are designed to make best of low frequency to enhance the robustness of the inversion.An overview of full waveform inversion can be found in [19].Driven by potential application, the development of FWI is very rapidly, for example, see [20] [21] [22].
In this paper, we develop an effective full waveform slowness inversion method for crosshole data by combing the total variation regularization with the constrained optimization together.We also propose the new skip algorithm to implement the constrained optimization problem.The rest of this paper is arranged as follows.The theoretical method is described in Section 2 in detail.It includes the forward method and the inverse method.In Section 3, numerical computations both for the noise free data and noisy date are implemented.Finally, the conclusion is drawn in Section 4.

Theory
There are three subsections in this section.In Section 2.1, the forward problem and the staggered-grid scheme with the perfectly matched layer are described.In Section 2.2, the inversion method and the corresponding algorithm by combing TV regularization with bound constraints are described.Section 2.3 involves the computations for the gradient of objective function and the step length, which is important to the success of inversion.

Forward Method
We consider the acoustic wave equation excited by the source ( ) x z , which the mathematical form of the pressure ( ) , in 0, , , , , 0, 0, on where is the area of computational domain, T is the final observation time, is the square of slowness parameter, v is the wave velocity.The system (1)-( 2) is the forward problem.In this paper, we solve the pressure u numerically with the finite difference method for its high computational efficiency.Among various of finite difference schemes, the staggered-grid scheme is a typical one [23].In order to apply the staggered-grid method, we introduce new variables x w and z w defined by u x ∂ ∂ and u z ∂ ∂ respectively.Then system (1)-( 2) can be written as the following hyperbolic form , Now we construct the difference scheme of (3)-( 5) on the staggered grids.The schematic map of grid points for different variables is shown in Figure 1.Let , n i j u be the wavfield at time n t ∆ and position ( ) , i x j z ∆ ∆ , and ( ) , ( ) ( ) ( ) ,  ( ) ( ) ( ) ( ) , with initial conditions 1 0 , , 0, 0, ( ) ( ) x z i j i j w w The sufficient and necessary stability for the scheme ( 6)-( 8) is (see Appendix Since the computational domain is limited, we need to use absorbing boundary conditions to eliminate the boundary reflections.Several methods can be applied, for example, the paraxial approximation [24], the exact nonreflecting conditions [25] [26] and the perfectly matched layer (PML) method [27].Here, we adopt the PML method in its split formulation.The idea of PML is to design a particular layer around the computational domain in which the waves will be attenuated satisfactory.In split PML formulation, the wavefield u is split into d z is introduced.Then system (3)-( 5) in the case of the source free can be written as the following PML form ( ) ( ) ( ) ( ) where ( ) d x and ( ) d z are designated to attenuate the refractions in absorbing zone.
We use the following model for the damping parameter ( ) ( ) ( ) where L is the thickness of PML absorbing layer, v is the media velocity, x or z is the distance from current position to PML inner boundary, and R is a parameter chosen as d p is zero in the original computational domain.The staggered-grid scheme for system ( 12)-( 15) is ( ) ( ) ( ) ( ) ( ) .
The forward computations can be extrapolated explicitly in the time direction based on ( 17)-( 21).In Figure 2, the snapshot of wave propagation in a homogeneous media with velocity 3000 m/s at time 0.8 s is shown.

Inversion with Constraints
Consider the inversion of slowness fields in the vector ( ) Define the objective function ( ) where the summation is for all the shots or sources in the well.In Equation ( 22), the first term comes from the misfit between the forward data u and the observed data obs u , the second term is the TV regularization with positive para- meter 2 ε , and 0 η > is the regularization parameter.The TV regularization is applied because it can yield better inversion result at jump discontinuities theoretically [30] [31].For more theory about Tikhonov theory, reader may refer to references, for example [32] [33].The inverse problem is to minimize the objective function, i.e., ( ) It is a high dimensional complex optimization problem.The dimensions are x z N N , where x N and z N are the discretization points along x and z directions respectively.In this paper, the iterative decent method named the global Barzilai-Borwein (GBB) method is applied because it is suitable for solving high-dimensional complex minimization problem [34] [35] [36].Generally, for a given nonlinear function be the initial slowness vector, then the ( ) where the superscript m denotes the iteration number, m α is the step length and f ∇ is the gradient of f at the point m σ .The solution of constraint problem for nonlinear optimization is challenging.
Usually, the penalty function technique is used to change the constrainted optimization to unconstrained optimization problem [37].In this paper, we develop a new effective algorithm called the skip method to solve it.The concept of the skip method is to skip the current inverse result if it goes outside of the bound constraints and simply store previous results.For our experience, the skip method performs better than the penalty method.Numerical computations in Section 3 will show the effectiveness of the skip method.

Computations for the Gradient and Step Length
For minimizing the problem (23), the gradient and the step length are required.
The gradient of the objection functional (22) with respect to σ is ( ) where φ is the backward wavefield satisfying which can be solved effectively like the forward problem in Section 2.1.The backward problem ( 26)-( 27) is introduced in the gradient derivation naturally in Appendix B. There are two terms on the right hand side of Equation (25).The first is the gradient of the misfit functional between the forward data u and the observed data obs u , which is derived in Appendix B. We remark that there is no Born approximation in the gradient derivation in Appendix B. The second is the gradient of TV regularization term and its discrete form can be obtained by the approximation of the gradient operator, i.e., ( ) ( ) , , , where The parameter ε in Equation ( 25) is a small number to prevent the singularity of the denominator.Now we consider the step length which is important to the success of inversion.For a general objective function φ : the one dimensional line search is applied to find the step length, where k p is the search direction and k α is the step length.Line search condition gives the step length k α which guarantees the sufficient decrease in the objective func- tion f by inequality ( ) ( ) for some constant ( ) . The reduction of the objective function is proportional to both the search direction and the step length.The sufficient condition states that k α is the step length only if ( ) ( ) The decrease is not enough if it only satisfies the sufficient condition because it can satisfy for the very small value of α .For this reason it requires another essential condition called the curvature condition.If k α satisfies ( ) for some constant ( ) , it is called to satisfy the curvature condition.The sufficient decrease condition and the curvature condition together are called Wolfe condition [37].It can also be possible a step length satisfies the Wolfe condition but it closes to a minimizer of φ .Due to this reason we use the mod- ified curvature condition to enlarge the step length which lies in a broad neighbourhood of local minimizer.Therefore, k α satisfies the strong Wolfe condi- tion if it satisfies the following two conditions ( ) where 1 2 0 1 c c < < < .For the importance of line search in inversion, the algo- rithm for line search is given by Algorithm 1 and Algorithm 2. The FWI algorithm with both TV regularization and bound constraints is given by Algorithm 3.

Numerical Computations
We implement numerical computations with C-language code.The computational domain is where 0 25 Hz f = is the peak frequency, A is the maximum amplitude and where , : 0.125 0.125 0.000851 As shown in Figure 3, this exact model contains a roughly spherical body with lower slowness in the middle on the background media with high slowness.We will consider the inversion in two cases: the noise free data and nosy data.

Noise Free Data
First we consider the inverse for the noise free data.The inverse result after 54 iterations without TV regularization and bound constraints in Algorithm 3 is shown in Figure 4. We remark that more iterations in Figure 4 is impossible because no step length can be found and the iteration is automatically stopped.
We can see the inverse result differs from the exact model very much, which

Noisy Data
The recorded data are usually affected by noise due to error in measurements or the influence of natural environment.So it is necessary to consider noise analysis.The Gaussian white noise is one of the effective way to add noise.It has the normal distribution and standard deviation of the data.Let the discrete signal be { } , m s Ν denotes a normal distribution with mean m and standard derivation 2  s , and ⋅ is the Euclidean mean.We first present some noisy data Here z N is the number of discretization points in z direction.In these figures, the red line denotes the noise free data while the blue line represents the noisy data and we can see clearly that the lower SNR is the higher noise in the data.Generally, the inversion for noisy data is challenge because the problem is nonlinear and ill-posed.However, numerical inverse results by our proposed method are still good.In Figure 11, the inverse results by three different methods for noisy data with SNR = 100 are given.13 respectively.After comparing these results, we can see that even in the case of lower SNR with SNR = 1, the inverse result with both TV regularization and bound constraints are still good, which shows the robustness of our proposed method.

Conclusion
Full waveform inversion is an effective method for recovering the media parameter.It is an optimization iterative process by minimizing the misfit between the recorded and synthesized data.In this paper, we have developed a new numerical method for crosshole waveform slowness inversion.This method combines TV regularization with constrained optimization together.The TV regularization method can yield better images of slowness at its discontinuities while the bound constraints from logging data can give more reliable constraints to keep the inversion converging to the right result.The combination is necessary to improve the robustness and inversion precision.The novel skip method is proposed to implement the constrained optimization problem effectively.Numerical computations both for noise free data and noisy data with lower SNR show the effectiveness and robustness of the proposed method.In the future, we will design preconditioners to improve the computational efficiency further.
where ( ) and The sufficient and necessary condition for stability is that the eigenvalues of matrix A is within or on the unit circle.The characteristic equation of matrix A is The stability requires which is just the stability (11).

Appendix B. The Derivation for the Gradient
For convenience, we denote the first term in Equation ( 22) by and the forward operator A by where A is given by The operator A is a self-conjugate operator.See the following Lemma 1.
Lemma 1.For any ( ) ( ) where H is the Sobolev space, then operator A is the self-adjoint operator, i.e., T Proof.Based on the definition of the operator A, we have Applying the Green formula in time and using the conditions (52), the first term in Equation (54) becomes Similarly, applying the Green formula in space and noting boundary conditions of u and φ , the second term in Equation (54) becomes The proof is completed.The result (69) is just the first term on the right hand side of Equation (25).We can see that there is no Born approximation in the derivation above.

7 )
Journal of Applied Mathematics and Physics

Figure 1 .
Figure 1.The schematic map of grid points for the staggered-grid scheme.
that the coefficient ( )

Figure 2 (
a) is the result without PML while Figure 2(b) is the result with PML.We can see that the serious boundary reflections in Figure 2(a) are eliminated effectively in Figure 2(b).
Now we propose to solve the minimization problem(23) as a constrained optimization problem subject to the constraints a b σ ≤ ≤ .Here a and b are the lower bound and upper bound for each component of the vector σ respectively.
is shown in Figure 3.The sources are located in the left well marked as * while the receivers are located in the right well marked as o.There are 27 sources in the left well and 29 receivers in the right well.The source function is the Ricker wavelet given by

Figure 3 .
Figure 3.The exact slowness model.The sources are in the left well marked as * and the receivers in the right well marked as o.
shows the inversion does not converge correctly.The inverse results with three different methods after 200 iterations are shown in Figure 5.For contrast, the exact model is shown in Figure 5(a) again.Figures 5(b)-(d) are the inverse results with TV regularization, with bound constraints, and with both TV regularization and bound constraints, respectively.Comparisons show that Figures 5(b)-(d) are much better than Figure 4.Moreover, Figure 5(d) with both TV regularization and bound constraints is best.Figure 6 is the convergence history W. S. Zhang, A. K. Joardar DOI: 10.4236/jamp.2018.650941096 Journal of Applied Mathematics and Physics

Figure 4 .
Figure 4. Inverse result without TV regularization and bound constraints.The inversion is stopped automatically after 54 iterations.

Figure 5 .Figure 6 .
Figure 5. Exact model and inverse results with three different methods after 200 iterations.(a) Exact model; (b) TV regularization; (c) Bound constraints; (d) Both TV regularization and bound constraints.
Gaussian noise which is random number with the mean zero and standard deviation under different signal-to-noise ratio (SNR)

Figure 7 .
Figure 7. Exact model and inverse results with three different methods after 1000 iterations.(a) Exact model; (b) TV regularization; (c) Bound constraints; (d) Both TV regularization and bound constraints.
Figure 11(a) is the exact model.Figures 11(b)-(d) are the inverse results after 1000 iterations with TV regularization, with bound constraints, and with both TV regularization and bound constraints, respectively.Comparing the results in Figure 11, we can clearly see that Figure 11(d) obtained with both TV regularization and bound constraints performs best.Correspondingly, the inverse results for the noisy data with SNR = 10 and SNR = 1 are shown in Figure 12 and Figure

Figure 8 .
Figure 8. Noisy data with SNR = 100 received in the right well at three different positions r z .(a) 2 r z z N = , (b)

Figure 9 .
Figure 9. Noisy data with SNR = 10 received in the right well at three different positions r z .(a) 2 r z z N = ; (b)

Figure 10 .
Figure 10.Noisy data with SNR = 1 received in the right well at three different positions r z .(a) 2 r z z N = ; (b)

Figure 11 .
Figure 11.Exact model and inverse results with three different methods after 1000 iterations for noisy data with SNR = 100.(a) Exact model; (b) TV regularization; (c) Bound constraints; (d) Both TV regularization and bound constraints.

Figure 12 .
Figure 12.Exact model and inverse results with three different methods after 1000 iterations f for noisy data with SNR = 10.(a) Exact model; (b) TV regularization; (c) Bound constraints; (d) Both TV regularization and bound constraints.

Figure 13 .
Figure 13.Exact model and inverse results with three different methods after 1000 iterations for noisy data with SNR = 1.(a) Exact model; (b) TV regularization; (c) Bound constraints; (d) Both TV regularization and bound constraints.

Let 1 Ι
= − and suppose x k and z k are the wavenumber in the x and z directions respectively.By applying , The proof is complete.Making variation for the objective function (49), we obtain

Theorem 2
If φ is the solution of the following backward problem