EAGE-HAGI 1st Asia Pacific meeting on Near Surface Geoscience & Engineering

11-12 April 2018, Yogyakarta, Indonesia

Introduction

Full-Waveform Inversion (FWI) can be defined as an iterative fitting data procedure for obtaining

physical properties of the Earth based on the full wave-field data simulation. Hence, FWI is widely

known as a comprehensive way to solve the complex structure below the earth surface, it performed a

high-resolution image, and very powerful when it combined with the good prior model. Generally,

FWI can be split up into two main problems, first is the forward modeling problem that consist of

computing the solution of the seismic wave for each source of the seismic experiment. Second,

Inverse problem that consist of estimating the model parameters that describing the subsurface from

the recorded data through the resolution of the optimization model (Virieux, 2009).

One of the most challenging steps of this methodology is the cost of computation. Comprehensive

calculation of the Gradient and the Hessian as the first and second derivatives of the cost function

must be the most challenging step on the whole process of FWI. To approach the computation of the

gradient, we used a backward propagation from the adjoint source that generated by the residual data.

However, it needs a high computational cost to store the huge amount of wavefield data. To get the

gradient computation done, the trajectory of the regular wavefield for each time and each shot are

needed to be stored, and then we convolved them with the adjoint wavefield.

In this research, we especially focusing on how we worked with the 2-D Full-Waveform Inversion

(FWI) using Gauss-Newton approach in elastic media. The steps included the forward modeling

problem based on the finite-difference and staggered grid scheme that bounded by free-surface

boundary condition (on the top of model) and Perfectly Matched Layer (PML) for the rest, also, we

applied the Gauss-Newton inversion that exploit the approximate-Hessian into this methodology. For

the result, we tested the inversion modeling in simple layer cake model and complex marmousi2

model, both of those models are located in the shallow area.

The Forward Problem

The forward modeling machine is greatly important to solving the inverse problem, it is needed to

calculate the predicted data based on our prior model. If the predicted data is failed to represent the

real behavior of the seismic wave, high possibility the objective function will produce the wrong

interpretation of the gradient of inversion. To anticipate the failure of gradient development, many

researchers have developed tons of forward modeling engines, one that we used for this paper is

called as a finite-difference approximation.

Consider a two-dimensional plane wave in homogeneous continuum. Using the stress-strain equation,

equation of continuity, and Cauchy’s equations of motion, we can write the resulting system of two

differential equations as,

Where is density, is particle displacement, is stress tensor, is body force, and is stiffness

tensor. The full staggered-grid scheme is applied in this research; it is very useful to maintain the

elastic parameters. In 1989, Virieux demonstrated in an elegant way how to put the density, lame

parameters, normal stress, shear stress, and velocity components in different position of grid.

Free-surface and Perfectly Matched Layer (PML)

Free-surface boundary condition is simply defined as a condition that maintain the seismic wave that

propagate through the top of solid interface. In the real condition, this medium can be as an air

condition or water condition. The implementation of this concept is assuming all normal and shear

stress that propagate perpendicular with z-plane in the edge of interface are zero. PML is slightly

EAGE-HAGI 1st Asia Pacific meeting on Near Surface Geoscience & Engineering

11-12 April 2018, Yogyakarta, Indonesia

different with other boundary condition, it is more similar with an absorbing layer than an absorbing

condition. Absorbing boundary layer is a layer of artificial absorbing material that is placed adjacent

to the edges of the grid, completely independent of the boundary condition. PML has been developed

by Berenger in 1994, it first applicable for electromagnetic wave and solve the free-space simulation.

Then, in 2001, Collino and Tsogka proposed the same absorbing layer, but it addressed for seismic

wave simulation in elastic media. We illustrate the damping profile in our computational domain,

simply, we can write the function of damping factor as,

Where, and are damping factor along horizontal and vertical axis, is depth of the barycenter

element inside the PML, and is reflection coefficient. We define the factor is zero inside the

interest area (physical model domain), and nonzero inside the PML area (non-physical model

domain). An absorbing layer is placed adjacent to the edges of the computational region.

Consequently, a perfect absorbing layer would absorb outgoing waves and suppress the undesirable

reflection waves into our interest domain.

The Inverse Problem

Instead being too much spend the cost computation for the Frechet, Adjoint State method that

proposed by Tarantola and Gauthier (1986) leads us to perform the simple cross-correlation between

the forward propagation with its pair (we called it as a backward propagation) to obtain the gradient.

Theoretically, seismic events are defined as the result of encounter event between up-going wave and

down-going wave. Based on these knowledge, we manipulate and produce those waves using the

forward propagation and backward propagation simulation. Thus, the principle formula for the

gradient can be written as,

Where is the constant dependent on its parameter. The arrow toward the right direction indicate the

model parameter of forward propagation, and the arrow toward the left direction means the backward

propagation. Consequently, we can directly calculate the forward propagation and the

backward propagation to obtain the gradient using finite-difference method.

Generally, Gauss-Newton methodology neglects the second derivative part of the objective function,

it simplified the Full-Hessian calculation into the inverse product between the first-derivative with its

conjugate. Both gradient and hessian are equipped by the damping constant factor, it proven effective

to stabilize the value.

Examples

We tested the full-waveform inversion (FWI) using two different synthetic models, the layer cake and

marmousi2 model. In the first model, it contains four different horizontal layers, 1500 m/s, 2000 m/s,

2500 m/s, and 3000 m/s, the second model contains the marmousi2 model with 3000-4500 m/s in

range. We also used 0.0002s of time step in our computation. PML (the left, right, and bottom of

model) and free-surface boundary condition (the top of model) are applied for this simulation to

mimic the real condition of earth surface. Figure 3 showed the result of layer cake model, it used a

EAGE-HAGI 1st Asia Pacific meeting on Near Surface Geoscience & Engineering

11-12 April 2018, Yogyakarta, Indonesia

homogeneous model as an initial model. The result successful to delineate the four different layers of

layer cake model while it also succeeds to produce the good result and show the same trend of P-wave

and S-wave velocity in complex marmousi2 model (see Figure 4).

Concluding Remarks

Like the other FWI methodology that exploit the Hessian, the advantage of the FWI using GaussNewton

on this research is capable to produce the high quality of resolution of the updated model, the

approximate Hessian leads the model to avoid the local minima efficiently (7 iterations in marmousi2

model and 10 iterations in layer cake model). However, it needs a high-end machine to accommodate

the huge amount storage of all wave trajectory for gradient and hessian calculation. Nevertheless, this

research used very simple configuration (9 sources, 9 receivers, and 50×50 of grids) to maximize the

potential of our laptop computer (8GB of RAM), but the technology is still promising to be developed

in the future, especially for the near surface investigation technology development.

References

Aki, K. and Richards, P. G., 2002, Quantitative seismology, second ed.: University Science Books.

Berenger, J. P., 1994, A perfectly matched layer for absorption of electromagnetic waves, Journal of

Computational Physics, 114, 185–200.

Collino, F. and Tsogka, C., 2001, Application of the perfectly matched absorbing layer model to the

linear elastodynamic problem in anisotropic heterogeneous media, Geophysics 66: 294–307.

Gauthier, O., Virieux, J., and Tarantola, Albert., 1986, Two-dimensional nonlinear inversion of

seismic waveforms: Numerical results, Geophysics, 51, 1387–1403.

Tarantola, A., 1984, Inversion of seismic reflection data in the acoustic approximation, Geophysics,

49, 1259–1266.

Virieux, J., 1986, P-SV wave propagation in heterogeneous media: velocity stress finite difference

method, Geophysics, 51, 889–901.

Virieux, J. and Operto, S., 2009, An overview of full-waveform inversion in exploration geophysics,

Geophysics, 74.

Figures

Figure 1 Examples of forward modeling simulation: Snapshot of normalize velocity wave field in

marmousi2 model. The snapshots captured in 0.15s, 0.30s, and 0.45s.

EAGE-HAGI 1st Asia Pacific meeting on Near Surface Geoscience & Engineering

11-12 April 2018, Yogyakarta, Indonesia

Figure 2 Snapshot of seismic wave field, boundary problem simulation, interest area is surrounded by

four-different boundary condition, Free-surface boundary condition in the top of model, Perfectly

Matched Layer (PML) in left-side of model, Absorbing Boundary Condition (ABC) in right-side of

model, and rigid-wall Dirichlet boundary condition with zero value along the edge in the bottom of

model.

Figure 3 The result of layer cake model (Lx=1000m, Lz=1000m, nx=50, nz=50, dt=0.0002s,

freq=2.5-15Hz, damping Hessian for lambda, mu, and rho= 5e-5, 2.5e-5, 5e-4, initial

model=homogeneous model, and convergent in 10 iteration)

Figure 4 The result of Marmousi2 model (Lx=1000m, Lz=1000m, nx=50, nz=50, dt=0.0002s,

freq=2.5-15Hz, initial model=gradation model, and convergent in 7 iteration)