# Simultaneous Inversion of Layered Velocity and Density Profiles Using Direct Waveform Inversion (DWI): 1D Case Zhonghan Liu, et al.

Jan 13, 2022

## Introduction

Seismic full waveform inversion (FWI), formulated originally by (Lailly, 1983; Tarantola, 1984), is a powerful process in subsurface velocity model building. The goal of FWI is to find a model such that the model-predicted waveforms fit the observed waveforms. Since FWI is an iterative gradient-based method, its success depends on how much the initial model differs from the true model (Virieux and Operto, 2009). The limitation of the iterative FWI scheme was recognized early on by many authors (Gauthier et al., 1986; Tarantola, 1986; Mora, 1987; Bourgeois et al., 1989). Tarantola (2005, p.128) pointed out that the local Fréchet gradient used in FWI was equivalent to the single scattering Born approximation. Therefore the performance of FWI relies on an accurate and long-wavelength initial velocity model in which case the Born approximation is more accurate. As the correspondence between low-frequency seismic data and low-wavenumber/large-scale structures is linear in the Born single scattering (Wu and Zheng, 2014), due to the lack of low-frequency content (<5 Hz) in most reflection seismic data, most developments in FWI have been focusing on how to recover large-scale structural information when low-frequency data are not available. These developments include, for example, the Laplace FWI (Shin and Cha, 2008; Shin and Ha, 2008; Kim et al., 2013), envelope inversion (Wu et al., 2014; Luo and Wu, 2015; Chen et al., 2018), intensity inversion (Liu et al., 2018; Liu et al., 2020), and the FWI using deep learning techniques (Richardson, 2018).

To circumvent the challenges in FWI, we proposed an alternative waveform inversion scheme (Liu and Zheng, 2015; 2017), called the direct waveform inversion (DWI), to invert for subsurface models without the need for a global initial model. DWI combines seismic imaging and velocity model building into one single process. In order to use DWI, it is necessary for the input seismic data to include both free-surface and inter-bed multiples. Using surface recorded reflection seismic data, DWI is able to deliver accurate P-wave velocity inversion results without using a global initial model, for both 1D and 2D layered scalar (i.e., no density variation) models (Zheng and Liu, 2020). Without using a global model, DWI inverts the model from shallow to deep depths. In this regard, DWI is similar to the layer-stripping methods (Claerbout, 1976) and the approach by Goupillaud (1961). However, there are important differences in the methods, in particular the explicit use of the time-space causality in DWI and local inversion in both space and time.

In the current industry, simultaneous inversion of multiple rock physical parameters is one of the state-of-art workflows that can much directly benefit the subsurface reservoir characterization and improve production (Brossier et al., 2009), using surface seismic data. To address the demand of multi-parameter inversion, many FWI methods are developed to invert for parameters like the P-wave velocity, S-wave velocity, density, and seismic anisotropy (Sears et al., 2008; Brossier et al., 2009; Jeong et al., 2012; Warner et al., 2013; Alkhalifah and Plessix, 2014). However, for these multi-parameter inversion methods based on FWI, the increased number of the model parameters means higher computational cost, and increased ill-posedness of the inverse problem (Virieux and Operto, 2009). In contrast, benefited from the localized inversion and explicit use of the time-space causality, our DWI method can be implemented for multi-parameter inversion without increasing much computational cost and numerical instability.

In this paper, we start from a 1D acoustic layered medium to demonstrate the ability of DWI in the simultaneous inversion of multiple parameters. Such 1D layered examples can fundamentally validate the feasibility of this state-of-art DWI method. The numerical examples demonstrate that the DWI method could play an important role in the multi-parameter seismic inversion.

### 1D Scalar DWI With Constant Density Throughout the Model

In this section, we briefly summarize the scalar DWI procedure for inverting the sound wave velocity in a horizontally stratified layered medium that has a constant density throughout the model (Figure 1), i.e.,

$ρi=ρ0$

. In the next section, we will consider the case where densities may be different in different layers.

FIGURE 1. Wavefields in a 1D layered model.

$c0$

,

$c1$

, … are layer velocities.

,

are densities in layers.

,

$H2$

are the layer thicknesses.

$Pi$

,

$vzi$

are pressure and vertical component of the particle velocity at the depth

$zi$

, respectively.

,

$Di$

are the up-going and down-going pressure fields respectively at depth

$zi$

, where “-” indicates the wavefields on the top side of

$zi$

, “+” indicates the wavefields on the bottom side of

$zi$

. Figure modified from Zheng and Liu (2020).

DWI explicitly uses the time-space causality property of the wavefield in the inversion. Starting from the source-receiver layer near the surface, we recursively (not iteratively) build the model downward by fitting the earliest parts of the waveforms of pairs of source-receivers of short (or zero) offsets. We then extrapolate the sources and receivers downward to the bottom of the inverted region, and repeat the process.

To illustrate the DWI process, we assume that both the pressure waveform, P, and particle velocity,

$Vz$

, are recorded on the surface (

$z0$

). The incident plane wave is vertical or at zero incident angle. We further assume that

$c0$

, the velocity of the first layer, is known. We decompose the wavefields (P and

$Vz$

) into up-going, U, and down-going, D, pressure wavefields, respectively (Liu and Zheng, 2015) as follows

where

$ρ$

and

$c$

are density and wave velocity, respectively. Conversely, if we know U and D, we can compose P and

$Vz$

. In this section, we assume a constant density profile

$ρ$

.

DWI consists of four steps:

• Step 1. At the acquisition plane depth

${z}_{0}$

, we use the recorded waveforms,

${P}_{0}\left(t\right)$

and

${V}_{z0}\left(t\right)$

, and the first layer’s velocity

${c}_{0}$

(assumed to be known), to calculate the up-going and down-going pressure fields at depth

${z}_{0}^{+}$

, denoted as

${U}_{0}^{+}$

and

${D}_{0}^{+}$

, respectively, using Eq. 1. From the point of view of the input/output system, we can view the medium below

${z}_{0}$

as a linear system.

${D}_{0}^{+}$

is the incident wave (or input), and

${U}_{0}^{+}$

is the reflection response (or output) of the system. Following the time-space causality of the wavefield, the earliest up-going impulse in

${U}_{0}^{+}$

must be generated from the first (earliest) down-going impulse in

${D}_{0}^{+}$

, reflected by the reflector at depth

${z}_{1}$

to be determined. A causal time-domain deconvolution between

${D}_{0}^{+}$

and

${U}_{0}^{+}$

can generate a response consisting of a series of impulses. The time of the first impulse gives a time difference,

$2\tau$

. Hence the depth of the reflector

${z}_{1}$

, or the thickness of the first layer

${H}_{1}$

, can be calculated by multiplying

${c}_{0}$

with the one-way traveltime

$\tau$

.

• Step 2. We then extrapolate fields

${U}_{0}^{+}$

and

${D}_{0}^{+}$

to the bottom of the first layer (depth

${z}_{1}^{-}$

) in the frequency (

$\omega$

) domain

$U1−=U0+exp(−iωτ)(3)$

This wavefield extrapolation can be done by many different methods. In 1d, we choose the phase shift approach shown here because it is easy to implement.

• Step 3. After extrapolation, the first impulse in

${U}_{0}^{-}$

and the first impulse in

${D}_{0}^{-}$

should be time-shifted to the same time as if the incidence and reflection occur right above

${z}_{1}$

. We use their amplitude ratio,

${R}_{0}$

, which is the reflectivity in a constant-density medium

To determine the velocity,

$c1$

, of the next (or the second) layer since

$c0$

is known.

• Step 4. Finally, we can use

${c}_{1}$

,

${U}_{1}^{-}$

, and

${D}_{1}^{-}$

, to obtain the pressure and particle velocity fields,

${P}_{1}$

and

${v}_{z1}$

, respectively at depth

${\text{z}}_{1}^{-}$

, using Eqs 1, 2. Because the pressure and the particle velocity fields should be continuous across a boundary, we can get their values at

${z}_{1}^{+}$

in layer 2. At this point, we also know

${c}_{1}$

, so our situation is the same as in Step 1.

The aforementioned process, using the recorded fields,

$P0$

and

$Vz0$

, and

$c0$

, to obtain the other parameters of the second layer (

$z1$

,

$P1$

,

$Vz1$

, and

$c1$

), can be recursively repeated downward layer by layer. As the inversion process goes deeper, there will be fewer and fewer remaining seismic events in both the up-going data. Eventually, DWI stops when there are no seismic events in the extrapolated upgoing fields due to finite recording time of the seismic traces. At this point, all the layers have been inverted or the inverted model has expanded downward to its maximum extent and converged to the final model. In this sense, DWI always converges, unconditionally.

### Simultaneous DWI for Both Velocity and Density

In the previous 1D DWI scheme, we assume the density is constant throughout the model. For the 1D inversion on models of depth dependent density profiles, there were some relevant work by Coen in 1980s (Coen, 1981a; Coen, 1981b; Coen, 1981c). In Coen’s work, the density and velocity are inverted separately using a dataset from oblique incident plane waves based on the Gel’fand-Levitan-Marchenko (GLM) theory (Agranovich and Marchenko, 1963; Berryman and Greene, 1980). In our study, instead of applying the GLM theory, we directly use the incident angle (

$θ$

)-dependence of the reflectivity,

$R(θ)$

, to invert for both velocities and densities of a layered model. To achieve simultaneous inversion of velocities and densities, we show how to modify the steps in the previous section respectively.

Assuming the wave is incident from medium-1 (

$ρ1,c1$

) at an angle

$θ$

to medium-2 (

$ρ2,c2)$

, we have the angle-dependent reflectivity

$R(θ)=ρ2c2cosθ−ρ1c12−c22sin2θρ2c2cosθ+ρ1c12−c22sin2θ(6)$

If we have two plane waves of two different incident angles

$θ1$

and

$θ2$

and two measured amplitudes,

$R1=R(θ1)$

and

$R2=R(θ2)$

, we can in principle determine

$c2$

, and

$ρ2$

, simultaneously.

To increase robustness of the inversion, we can make use of waves of multiple incident angles (

$n≥2$

), and minimize the objective function

$J[c2,ρ2]=∑i=1n|R(θi;ρ2,c2)−ri|2(7)$

In Eq. 7, for a plane wave at the incident angle

$θi$

,

$R(θi)$

is theoretically modeled reflectivity using Eq. 6, and

$ri$

is the measured reflectivity. As we only have two unknowns (

$c2$

and

$ρ2$

), using a grid search method can quickly get the results.

Assuming the incident wave angle is

$θ$

and in order to invert for the density profile, the steps of the scalar DWI need to be modified as follows:

In Step 1, Eq. 2, the relationship between the pressure and vertical component of the particle velocity, should be changed to

In Step 2, the extrapolation of up-going and down-going pressure fields should be modified as

$U2−=U1+exp(−iωτcosθ)(9)$
$D2−=D1+exp(+iωτcosθ)(10)$

In Step 3, using the amplitude ratio

$R(θi)$

and the measured data,

$ri$

, from multiple incident angles,

$θi$

, we can obtain

$c2$

and

$ρ2$

by either solving Eq. 6 directly or fitting Eq. 7.

In Step 4, we need to use Eqs 1, 8 to compose P and

$Vz$

from the up-going and the down-going fields.

### Numerical Examples

In this section, we will present two synthetic examples to demonstrate the effectiveness of our proposed method: a simple layered model with six layers, and a more complex layered model with thirty-one layers. Both models are horizontally stratified. Within each layer, the velocity and density are constant. However, different layers have different properties. Both the top and bottom boundaries of the model are set up as half-space boundary conditions.

The synthetic data (pressure and particle velocity) in both examples are generated by a propagator matrix method (e.g., Eftekhar et al., 2018). The plane wave is injected at a depth of 0 m and propagated downward into the model. The receivers are placed at the same depth. Both the pressure and particle velocity wavefields are recorded at a time sampling interval of 1 ms.

Example 1 In the first example, there are six layers (Figure 2) and the velocity contrast is up to 200%. Here we use a 15 Hz Ricker wavelet as the incident plane wave for the model (Figure 2). We conduct the modeling for four plane wave sources at different incident angles: 0, 5, 10, and 15°. The waveform records of the pressure and the vertical component of particle velocity are shown in Figure 3. In Figure 3, the recorded waveforms contain full information of the wavefields, including the primary reflections and multiples.Using the recorded data shown in Figure 3 and following the DWI steps in the previous section, we inverted for both the velocity and density profiles, shown in Figure 4. We also calculated the misfits in velocities and densities between the DWI results and the true models shown in Table 1, where most of them are less than 1%, except the velocity in the last layer (layer 6).To check the validity of the inverted model in the data space, we conduct a forward synthetic modeling using the DWI inverted model (Figure 4). The modeled waveforms fit the data very well (Figure 5). Both the primary reflections and the internal multiples can be reproduced by the DWI inverted model.

FIGURE 2. Velocity (A) and density (B) profiles of the true model.

FIGURE 3. Recorded waveforms of pressure (A–D) and vertical component of particle velocity (E–H) in response to four different plane waves.

FIGURE 4. Comparisons between DWI inversion result and true model for both velocity (A) and density (B) structure.

TABLE 1. The misfit of velocities and densities between the inverted model and true model (the velocity and density information of the first layer are known).

FIGURE 5. Comparisons of the recorded data (red) and synthetics (black) modeled using the DWI inverted model. (A) pressure waveforms; (B) particle velocity waveforms at the 0-degree incidence. Note the waveform amplitudes of the multiples in the red dashed box are amplified by 300 times.

Example 2 In the second example, we build a model with 31 velocity and density layers (Figure 6). An 80 Hz Ricker wavelet was used as the incident plane wave source wavelet and we modeled the data for four plane waves at angles: 0, 5, 9, and 16°. The recorded waveform data of the pressure and the vertical particle velocity are shown in Figure 7.Compared with the recorded waveforms in the first example (Figure 3), both the primary reflections and the internal multiples (Figure 7) are much more complicated. Using these recorded data, we applied our DWI scheme and obtained the inversion results of velocity, density, and impedance models shown in Figure 8.From Figure 8 we can see that the DWI scheme almost exactly recovers the impedance model. For the velocity and density models, although there are some small misfits (less than 2%), the inverted models still agree well with the true model. To further examine the influences of these modeled misfits, here we also conduct a forward synthetic modeling based on the inverted model (Figure 8). The results are shown in Figure 9. The modeled waveforms using the DWI model fit the data very well including not only the primary reflections but also the internal multiples.

FIGURE 6. The velocity (A) and density (B) profiles of the true model.

FIGURE 7. Recorded waveforms of pressure (A–D) and vertical component of particle velocity (E–H) in response to four different plane waves.

FIGURE 8. Comparisons between DWI inversion result and true model on velocity (A), density (B), and impedance (C) models.

FIGURE 9. Comparisons of data (red) and synthetics (black) modeled using the DWI inverted model for pressure at the 0-degree incidence (A), the events in the red dashed box are amplified by 10 times. A zoom-in view of the events in the red dashed box in (A) is shown in (B).

## Discussion

We remark on limitations of DWI. DWI depends on reflection events in data to invert for the subsurface model. If the true model is smooth and does not have many reflectors, DWI may fail to find the true model. If the incident angle is too large and the total reflection occurs, DWI is not able to invert for model parameters below the total reflection depth.

The performance of the recursive DWI scheme may suffer from the accumulation of errors as the inversion process goes from shallow to deep depths. Data redundancy can help. In order to resolve reflectivities into P-wave velocity contrast and density contrast, respectively, we need to use several distinct incident angles. If the range of incident angles is narrow, DWI may not be able to resolve Vp and density correctly.

It is also worthwhile to provide some general remarks on the differences between DWI and FWI. Based on the single-scattering approximation (Tarantola, 2005), FWI linearizes the global non-linear seismic inversion problem by iteratively minimizing the misfit between the recorded data and the model predicted data. For FWI, if the intial model is far from the true model, the FWI convergence can be a problem (Sirgue and Pratt, 2004). Different from FWI, DWI does not rely on an initial global model to start the waveform inversion process. It only needs the local velocity around the surface receivers. The DWI scheme converts the FWI global optimization problem into many localized reflectivity inversions by explicitly invoking the causality principle. Hence it reduces the nonlinearity significantly which is a strength over the FWI method. Another advantage for DWI is that it does not need low-frequency data as seen in our example 2. On the other hand, if the true model is smooth and data have only a few reflection events, FWI may perform better. In cases where DWI can be applicable, DWI can be significantly faster numerically to obtain a model.

There are also marked differences between DWI and the 1D inversion using the GLM theory. Most developments of the GLM theory are aiming at imaging and redatuming methods (e.g., Broggini et al., 2014; Wapenaar et al., 2014; van der Neut et al., 2015; Nowack and Kiraz, 2018). Recently, Wu and He (2020) used the GLM theory to invert for 1D impedance profile. For 1D GLM problem, a time to depth conversion is needed and is usually carried out by the Liouville transform. But for 2D and 3D spatial problems, a macro-velocity model need to be used and should be obtained a prior. However, for DWI, the inversion is localized in a shallow to deep fashion and the inverted model is automatically obtained in the depth domain and there is no need to use a global velocity model. In future we will show 2D inversion results in which the wavefield extrapolation is much more involved using integral equations.

## Conclusion

We extend the scalar DWI scheme to invert for the subsurface density and velocity simultaneously, using multiple plane waves. Using recorded seismic data on the surface, our method inverts for the model parameters locally by explicitly employing time-space causality principle of the wavefield and recursively from shallow to deep depths. The new DWI scheme makes use of the angle dependent property of the reflectivity to solve for density and velocity simultaneously. The input seismic data to DWI must include all types of data including multiples. Numerical examples demonstrate the feasibility of the DWI approach to invert for both velocity and density using four plane wave sources. We find that the acoustic impedance profile is better resolved than the P-wave velocity or density owing to slight tradeoff (<1%) between the two parameters.

## Data Availability Statement

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

## Author Contributions

YZ, H-WZ, and ZL contributed to idea and design of the study. ZL realized the idea and generated numerical examples for testing. ZL wrote the first draft of the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.

## Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.