Three-dimensional frequency-domain full waveform inversion based on the nearly-analytic discrete method

DeYao Zhang; WenYong Pan; DingHui Yang; LingYun Qiu; XingPeng Dong; WeiJuan Meng

doi:10.26464/epp2021022

Earth and Planetary Physics > 2021 > 5(2): 149-157

Zhang, D. Y., Pan, W. Y., Yang, D. H., Qiu, L. Y., Dong, X. P. and Meng, W. J. (2021). Three-dimensional frequency-domain full waveform inversion based on the nearly-analytic discrete method. Earth Planet. Phys., 5(2), 149–157. DOI: 10.26464/epp2021022

Citation:

RESEARCH ARTICLE | SOLID EARTH: COMPUTATIONAL GEOPHYSICS Open Access

PDF (910 KB)

Three-dimensional frequency-domain full waveform inversion based on the nearly-analytic discrete method

1.
Department of Mathematical Sciences, Tsinghua University, Beijing 100084, China
2.
Key Laboratory of Petroleum Resource Research, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China
3.
Yau Mathematical Sciences Center, Tsinghua University, Beijing 100084, China

More Information

Corresponding author:
D. H. Yang, ydh@tsinghua.edu.cn

Publication History:

Issue Online: February 28, 2021
First Published online: February 28, 2021
Accepted article online: March 01, 2021
Article accepted: January 26, 2021
Article received: September 29, 2020

Key Points:

Nearly analytic discrete method is applied to solving 3D frequency-domain acoustic wave equation.

Compared to traditional finite difference methods, the nearly analytic discrete method can solve the wave equation more accurately and efficiently.

Nearly analytic discrete method is more suitable than traditional methods for 3D frequency-domain full waveform inversion.

Abstract

Abstract

The nearly analytic discrete (NAD) method is a kind of finite difference method with advantages of high accuracy and stability. Previous studies have investigated the NAD method for simulating wave propagation in the time-domain. This study applies the NAD method to solving three-dimensional (3D) acoustic wave equations in the frequency-domain. This forward modeling approach is then used as the “engine” for implementing 3D frequency-domain full waveform inversion (FWI). In the numerical modeling experiments, synthetic examples are first given to show the superiority of the NAD method in forward modeling compared with traditional finite difference methods. Synthetic 3D frequency-domain FWI experiments are then carried out to examine the effectiveness of the proposed methods. The inversion results show that the NAD method is more suitable than traditional methods, in terms of computational cost and stability, for 3D frequency-domain FWI, and represents an effective approach for inversion of subsurface model structures.
- three-dimension,
- frequency-domain,
- NAD method,
- forward modeling,
- full waveform inversion

FullText(HTML)

1. Introduction

In recent decades, full waveform inversion (FWI) has become increasingly popular in global and exploration seismology due to its ability to provide high-resolution subsurface velocity structures (Lailly, 1983; Tarantola, 1986; Virieux and Operto, 2009). FWI can be performed both in the time domain (Mora, 1987; Tromp et al., 2005; Liu QY and Tromp, 2006; Tape et al., 2007; Bozdağ et al., 2016; Wang J et al., 2019; Dong XP et al., 2019) and the frequency-domain (Pratt and Worthington, 1990; Song and Williamson, 1995; Pratt, 1999; Operto et al., 2004; Sirgue and Pratt, 2004; Brossier et al., 2009). Compared with time-domain FWI, frequency-domain FWI is more flexible in selecting specific frequency components in a multi-scale inversion strategy. In addition, frequency-domain FWI is computationally more efficient, thanks to solving different sources simultaneously with one matrix factorization. However, in 3D frequency-domain FWI the matrix size grows exponentially (Operto et al., 2007; Plessix, 2009). The iterative solution of the seismic inverse problem requires a large number of forward simulations. Forward modeling of the seismic waves is considered as the “engine” of and foundation for the implementation of FWI. Therefore, determining an accurate and efficient forward modeling approach is necessary and essential if FWI experiments are to be practical.

In the past decades, researchers have studied a series of forward modeling methods, including finite difference methods (Kelly et al., 1976; Day, 1982; Igel et al., 1995; Yang DH et al., 2003, 2004, 2006; Liu SL et al., 2015, Ma X et al., 2018), finite element methods (Lysmer and Drake, 1972; Marfurt, 1984; Liu SL et al., 2014; He XJ et al., 2020), a pseudo-spectral method (Kosloff and Baysal, 1982), and a spectral element method (Kosloff and Baysal, 1982; Seriani and Priolo, 1994; Komatitsch et al., 2005). These methods are designed to solve different types of questions. Each may have superiorities in some aspects but limitations in others. Among these methods, finite difference methods are widely used because of the advantages of easy implementation, low computational cost, etc.

The nearly analytic discrete (NAD) method is a finite difference method that uses displacement and its gradient to approximate the high-order partial derivatives in the wave equation (Yang DH et al., 2003, 2004, 2006; Tong P et al., 2011, 2013). As this method captures more detailed wave equation information, dispersion analysis and numerical experiments have shown that it can solve the wave equation more accurately and efficiently than traditional finite difference methods (Yang DH et al., 2010, 2012). Previous studies have implemented the NAD method for forward modeling and reverse time migration in time-domain applications (Yang DH et al., 2003; Li JS et al., 2013). This study extends the NAD method to solve the three-dimension (3D) acoustic wave equation in the frequency-domain, which is then used as the “engine” for conducting 3D frequency-domain FWI. The paper is organized as follows. First, we introduce the 3D frequency-domain acoustic wave equation's discretization with a perfect matched layer (PML) absorbing boundary condition, which gives the linear system equations. According to the characteristics of the coefficient matrix (or “impedance matrix”), we choose the Krylov subspace iteration with incomplete LU decomposition to solve the linear system. We then introduce the theoretical basis of 3D FWI in the frequency-domain. The non-linear conjugate gradient method is adopted for model updating in the iterative inversion process. Finally, we give numerical examples of forward modeling and 3D frequency-domain FWI. The forward modeling experiments show that the NAD method for the 3D frequency-domain acoustic wave equation works well and can suppress the numerical dispersion effects more effectively than traditional finite difference methods. The FWI experiments show that the frequency-domain FWI with the NAD method can efficiently and stably reconstruct the subsurface velocity structures.

2. Forward Modeling with the NAD Method in Frequency-Domain

In the frequency-domain, the 3D acoustic wave equation with constant density can be written as

$\Delta u\left({x,y,z} \right) + \frac{{{\omega ^2}}}{{{c^2}}}u\left({x,y,z} \right) = - \frac{1}{{{c^2}}}s\left({x,y,z} \right),\ \ \ \left({x,y,z} \right) \in D,$

(1)

where $u$ is pressure wavefield, $\Delta$ is Laplace operator, $\omega$ indicates angular frequency, $c = c\left({x,y,z} \right)$ is the compressional velocity of the medium, $s$ is the seismic source term, and $D$ denotes the whole 3D computing area.

To discretize the acoustic wave equation with the NAD method, we need to derive the partial derivative form of the frequency-domain equation with respect to x, y and z, respectively, as illustrated in the following equations:

$\left\{ \begin{aligned} &{\frac{{{\partial ^2}u}}{{\partial {x^2}}} + \frac{{{\partial ^2}u}}{{\partial {y^2}}} + \frac{{{\partial ^2}u}}{{\partial {z^2}}} + \frac{{{\omega ^2}}}{{{c^2}}}u = - \frac{1}{{{c^2}}}s}, \\ &{\frac{{{\partial ^2}{u_x}}}{{\partial {x^2}}} + \frac{{{\partial ^2}{u_x}}}{{\partial {y^2}}} + \frac{{{\partial ^2}{u_x}}}{{\partial {z^2}}} + \frac{{{\omega ^2}}}{{{c^2}}}{u_x} = - \frac{1}{{{c^2}}}\frac{{\partial s}}{{\partial x}}}, \\ & {\frac{{{\partial ^2}{u_y}}}{{\partial {x^2}}} + \frac{{{\partial ^2}{u_y}}}{{\partial {y^2}}} + \frac{{{\partial ^2}{u_y}}}{{\partial {z^2}}} + \frac{{{\omega ^2}}}{{{c^2}}}{u_y} = - \frac{1}{{{c^2}}}\frac{{\partial s}}{{\partial y}}}, \\ & {\frac{{{\partial ^2}{u_z}}}{{\partial {x^2}}} + \frac{{{\partial ^2}{u_z}}}{{\partial {y^2}}} + \frac{{{\partial ^2}{u_z}}}{{\partial {z^2}}} + \frac{{{\omega ^2}}}{{{c^2}}}{u_z} = - \frac{1}{{{c^2}}}\frac{{\partial s}}{{\partial z}}}. \end{aligned} \right.$

(2)

Second, we apply the PML condition (Komatitsch and Tromp, 2003), which is widely used in the area of seismic wave propagation, in Equation (2) for eliminating the reflections from artificial boundaries. The real coordinate $x$ in the original differential equations is replaced with the complex coordinate $\tilde x$ in the PML region. The relation between $x$ and $\tilde x$ is:

$\tilde x = x - \frac{\rm i}{w}\int\limits_0^x {{d_x}(s){\rm d}s,}$

(3)

where ${d_x}(s) > 0$ is an attenuation function, and $i$ denotes the imaginary unit. We choose ${\rm d}s: = {d_x}(s)$ as (Komatitsch and Tromp, 2003):

${d_x}(s) = - \frac{{3c}}{{2l}}\log (R){\left({\frac{s}{l}} \right)^2}.$

(4)

In Equation (4), $c$ is the wave propagation velocity between the PML region and real physical region, $l$ indicates the depth of the PML region, $s$ is the position in the PML layer, and $R$ is theoretical reflection parameter which can be very small (in our models, it is set as 0.001).

In the PML region, after replacing x, y and z with $\tilde x$ , $\tilde y$ and $\tilde z$ , Equation (2) can be rewritten as:

$\left\{ \begin{aligned} & {\frac{{{\partial ^2}u}}{{\partial {{\tilde x}^2}}} + \frac{{{\partial ^2}u}}{{\partial {{\tilde y}^2}}} + \frac{{{\partial ^2}u}}{{\partial {{\tilde z}^2}}} + \frac{{{\omega ^2}}}{{{c^2}}}u = - \frac{1}{{{c^2}}}s}, \\ & {\frac{{{\partial ^2}{u_{\tilde x}}}}{{\partial {{\tilde x}^2}}} + \frac{{{\partial ^2}{u_{\tilde x}}}}{{\partial {{\tilde y}^2}}} + \frac{{{\partial ^2}{u_{\tilde x}}}}{{\partial {{\tilde z}^2}}} + \frac{{{\omega ^2}}}{{{c^2}}}{u_{\tilde x}} = - \frac{1}{{{c^2}}}\frac{{\partial s}}{{\partial \tilde x}}}, \\ & {\frac{{{\partial ^2}{u_{\tilde y}}}}{{\partial {{\tilde x}^2}}} + \frac{{{\partial ^2}{u_{\tilde y}}}}{{\partial {{\tilde y}^2}}} + \frac{{{\partial ^2}{u_{\tilde y}}}}{{\partial {{\tilde z}^2}}} + \frac{{{\omega ^2}}}{{{c^2}}}{u_{\tilde y}} = - \frac{1}{{{c^2}}}\frac{{\partial s}}{{\partial \tilde y}}}, \\ &{\frac{{{\partial ^2}{u_{\tilde z}}}}{{\partial {{\tilde x}^2}}} + \frac{{{\partial ^2}{u_{\tilde z}}}}{{\partial {{\tilde y}^2}}} + \frac{{{\partial ^2}{u_{\tilde z}}}}{{\partial {{\tilde z}^2}}} + \frac{{{\omega ^2}}}{{{c^2}}}{u_{\tilde z}} = - \frac{1}{{{c^2}}}\frac{{\partial s}}{{\partial \tilde z}}}. \end{aligned} \right.$

(5)

Then we turn the variables $\tilde x$ , $\tilde y$ and $\tilde z$ in Equation (5) to x, y and z. By referring to Equation (3), we see that the derivative can be computed as

$\frac{{\partial x}}{{\partial \tilde x}} = \frac{{{\rm i}\omega }}{{{\rm i}\omega + {d_x}}}.$

(6)

Applying chain rule to Equation (5), taking $\tilde x$ as an example, the result is given in the following equation:

$\left\{ \begin{aligned} & {\frac{{\partial u}}{{\partial \tilde x}} = \frac{{{\rm i}\omega }}{{{\rm i}\omega + {d_x}}}\frac{{\partial u}}{{\partial x}}}, \\ & {\frac{{{\partial ^2}u}}{{\partial {{\tilde x}^2}}} = {{\left({\frac{{{\rm i}\omega }}{{{\rm i}\omega + {d_x}}}} \right)}^2}\frac{{{\partial ^2}u}}{{\partial {x^2}}} - \frac{{{{({\rm i}\omega)}^2}d'_{xp}}}{{{{({\rm i}\omega + {d_x})}^3}}}\frac{{\partial u}}{{\partial x}}}, \\ &{\frac{{{\partial ^3}u}}{{\partial {{\tilde x}^3}}} = {{\left({\frac{{{\rm i}\omega }}{{{\rm i}\omega + {d_x}}}} \right)}^3}\frac{{{\partial ^3}u}}{{\partial {x^3}}} - \frac{{3{{({\rm i}\omega)}^3}d'_{xp}}}{{{{({\rm i}\omega + {d_x})}^4}}}\frac{{{\partial ^2}u}}{{\partial {x^2}}} + \left({\frac{{3{{({\rm i}\omega)}^3}{d'^2_{xp}}}}{{{{({\rm i}\omega + {d_x})}^5}}} - \frac{{{{({\rm i}\omega)}^3}d''_{xpp}}}{{{{({\rm i}\omega + {d_x})}^4}}}} \right)\frac{{\partial u}}{{\partial x}}}, \\ & {\frac{{{\partial ^3}u}}{{\partial {{\tilde y}^2}\partial \tilde x}} = \frac{{{{({\rm i}\omega)}^3}}}{{({\rm i}\omega + {d_x}){{({\rm i}\omega + {d_y})}^2}}}\frac{{{\partial ^3}u}}{{\partial {y^2}\partial x}} - \frac{{{{({\rm i}\omega)}^3}d'_{yp}}}{{({\rm i}\omega + {d_x}){{({\rm i}\omega + {d_y})}^3}}}\frac{{{\partial ^2}u}}{{\partial y\partial x}}}, \\[-10pt]\end{aligned}\right.$

(7)

where $d'_{xp}$ and $d''_{xpp}$ are the first-order and second-order derivatives of ${d_x}$ . We can replace $\tilde x$ , $\tilde y$ and $\tilde z$ to obtain the corresponding partial derivatives, then apply all of the partial derivatives to Equation (5). The final 3D acoustic wave equation of the NAD method with PML boundary condition is obtained. The whole discretized wave equation is illustrated in the Appendix. In the end, we apply the fourth-order NAD method to discretize the 3D wave equation following Lang C and Yang DH (2017). After discretization with the fourth-order NAD method, the 3D frequency-domain acoustic wave equation can be written in compact matrix form, known as a linear system equation:

${\boldsymbol{AU}}={\boldsymbol{s}},$

(8)

where ${\boldsymbol{A}}$ is the impedance matrix, ${\boldsymbol{U}}$ is the pressure wavefield vector, and ${\boldsymbol{s}}$ is the source term vector. In the NAD method, the wavefield includes its spatial derivatives ${\boldsymbol{U}} = (u,{u_x},{u_y},{u_z})$ . Lang C and Yang DH (2017) have carried out many studies concerning solutions of Equation (8). In this study, to solve the linear system we use the Krylov subspace iteration (GMRES and BiCGSTAB iteration) with incomplete LU (ILU) decomposition precondition. A Ricker wavelet is used as the source time function:

$w(t) = A\left[ {1 - 2{{(\pi {f_0}t)}^2}} \right]{e^{ - {{(\pi {f_0}t)}^2}}},$

(9)

where $A$ is the amplitude of seismic source, ${f_0}$ denotes the central frequency, and $e$ is the Euler number. In the frequency-domain, the source wavelet is obtained as:

$F\left[ {w(t)} \right] = \frac{{\sqrt 2 A}}{{\pi {f_0}}}{\left({\frac{f}{{{f_0}}}} \right)^2}{e^{ - {{\left({\frac{f}{{{f_0}}}} \right)}^2}}},$

(10)

where $f$ denotes frequency.

3. Full Waveform Inversion

In this section, the theory of 3D frequency-domain FWI based on the NAD method is introduced. In FWI, the model properties are updated iteratively by minimizing an objective function that measures the observed and synthetic seismic data differences. The L-2 norm objective function is formulated as:

$E(c) = \frac{1}{2}\sum\limits_{\omega,i} {{{\left[ {\delta {d^{\,i}}(c,\omega)} \right]}^{\rm H}}\delta {d^{\,i}}(c,\omega),}$

(11)

where ${(\cdot)^{\rm H}}$ denotes the conjugate transpose of the matrix, and $\delta {d^{\,i}}$ is the data residual at the i-th source, which is obtained by

$\delta d_j^{\,i}(c,\omega) = u_j^i(c,\omega) - d_j^{\,i}({c_{\rm{true}}},\omega),\ \ \ j = 1,2, \cdot \cdot \cdot,n,$

(12)

where ${c_{\rm{true}}}$ means the real model (or target model), $n$ denotes the total number of receivers, $d_j^{\,i}({c_{\rm{true}}},\omega)$ is the observed data recorded on the j-th receiver of the i-th source, and $u_j^{\,i}(c,\omega)$ is the corresponding synthetic data.

The full waveform inverse problem is commonly solved within the gradient-based optimization framework. Different local optimization methods, including steepest descent, non-linear conjugate gradient (NCG), quasi-Newton methods, etc. (Nocedal and Wright, 2006; Pan WY et al., 2017), can be used for the model updating. Within these optimization methods, the first-order derivative of the misfit function, namely the gradient, needs to be calculated first. Following Pratt (1999), the gradient of the misfit function $E(c)$ can be written as:

$\nabla E(c) = \frac{{\partial E}}{{\partial c}} = {\rm{Re}} ({{\boldsymbol{J}}^{\rm T}}\delta {{\boldsymbol{d}}^{*}}),$

(13)

where ${\boldsymbol{J}}$ is the partial derivative of the synthetic wavefield data ${\boldsymbol{U}} = {({u_1},{u_2}, \cdot \cdot \cdot,{u_n})^{\rm T}}$ ( ${(\cdot)^{\rm T}}$ and ${(\cdot)^{*}}$ mean the transpose and conjugate of vectors):

$\;\;\;\;\;\;{\boldsymbol{J}} = \dfrac{{\partial {\boldsymbol{U}}}}{{\partial c}} = \left({\begin{array}{*{20}{c}} {\dfrac{{\partial {u_1}}}{{\partial {c_1}}}}&{\dfrac{{\partial {u_1}}}{{\partial {c_2}}}}&{ \cdot \cdot \cdot }&{\dfrac{{\partial {u_1}}}{{\partial {c_m}}}} \\ \vdots & \vdots & \vdots & \vdots \\ {\dfrac{{\partial {u_n}}}{{\partial {c_1}}}}&{\dfrac{{\partial {u_n}}}{{\partial {c_2}}}}&{ \cdot \cdot \cdot }&{\dfrac{{\partial {u_n}}}{{\partial {c_m}}}} \end{array}} \right) = \left[ {\begin{array}{*{20}{c}} {\dfrac{{\partial u}}{{\partial {c_1}}}}&{\dfrac{{\partial u}}{{\partial {c_2}}}}& \cdots &{\dfrac{{\partial u}}{{\partial {c_m}}}} \end{array}} \right].$

(14)

However, it is computationally unaffordable to calculate the partial derivative wavefield explicitly for large-scale inverse problems. With the help of the linear system in Equation (8), the partial derivative wavefield can be derived as (Pratt et al., 1998; Pan WY et al., 2016):

$\frac{{\partial {\boldsymbol{U}}}}{{\partial c}} = - {{\boldsymbol{A}}^{ - 1}}\frac{{\partial {\boldsymbol{A}}}}{{\partial c}}{\boldsymbol{U}}.$

(15)

Substituting Equation (15) into Equation (14) then into Equation (13), we can obtain the gradient as:

$\nabla E(c) = - {\rm{Re}} \left[ {{{\boldsymbol{U}}^{\rm T}}{{\left({\frac{{\partial {\boldsymbol{A}}}}{{\partial c}}} \right)}^{\rm T}}{{\boldsymbol{A}}^{ -{\rm T}}}\delta {{\boldsymbol{d}}^{*}}} \right],$

(16)

where $\dfrac{{\partial {\boldsymbol{A}}}}{{\partial c}}$ , ${\boldsymbol{U}}$ and $\delta {{\boldsymbol{d}}^*}$ can be obtained easily, but for large models the cost of calculating ${{\boldsymbol{A}}^{ - {\rm T}}}$ directly is impracticable. Instead, the gradient can be computed efficiently by cross-correlating the forward and adjoint wavefields:

$\nabla E(c) = - {\rm{Re}} \left[ {{{\boldsymbol{U}}^{\rm T}}{{\left({\frac{{\partial {\boldsymbol{A}}}}{{\partial c}}} \right)}^{\rm T}}{\boldsymbol{v}}} \right],$

(17)

where ${\boldsymbol{v}}$ is known as the adjoint wavefield obtained by solving the following adjoint wave equation:

${{\boldsymbol{A}}^{\rm T}}{\boldsymbol{v}} = \delta {{\boldsymbol{d}}^{*}}.$

(18)

The LU decompositions of ${\boldsymbol{A}}$ and ${{\boldsymbol{A}}^{\rm T}}$ are similar. In the previous section, we have obtained the ILU decomposition of ${\boldsymbol{A}}$ . Thus, Equation (19) can be solved using the same approach with the fourth-order NAD method. At each non-linear iteration, the model parameter $m$ can be updated by

${m_{k + 1}} = {m_k} + {\alpha _k}\delta {m_k},$

(19)

where ${\alpha _k}$ and $\delta {m_k}$ are the step length and search direction at the k-th iteration. The step length is calculated using the line search method. In this study, we adopt the NCG method for solving the inverse problem. The search direction is calculated by

$\delta {m_k} = - \nabla {E_k} + {\beta _k}\delta {m_{k - 1}},$

(20)

where ${\beta _k}$ is a scalar such that $\delta {m_k}$ and $\delta {m_{k-1}}$ are conjugate. It can be obtained with the “Fletcher-Reeves”, “Polak-Ribiere”, “Hestenes-Stiefel” or “Dai-Yuan” methods.

4. Numerical Experiments

This section first gives two forward modeling examples to illustrate the NAD method's advantages compared to traditional finite difference methods. We then present numerical examples of 3D frequency-domain FWI based on the NAD method.

4.1 Homogeneous Model Example

The 3D homogeneous model has a size of 4 km, 0.5 km, and 4 km in x, y and z directions. The grid spacing is 25 m. PML layers with a thickness of 6 are applied on all boundaries of the model for absorbing the reflections of artificial boundaries. The compressional velocity of the medium is 2 km/s. A Ricker wavelet with a dominant frequency of 30 Hz, located at the center of the model, is used as the source for forward modeling.

Figure 1 shows the frequency-domain wavefield snapshot at the slice of y = 0.25 km with a frequency of 30 Hz. Figure 2 shows the frequency-domain wavefield snapshots in 3D with a frequency of 11 Hz. Figures 3, 4 and 5 are the wavefield snapshots in the time domain. In Figure 3, at the time of t₂ = 1.233 s, the wavefronts have spread out of the computation region without evident reflections, indicating that the PML boundary condition works effectively. However, it can be seen in Figure 5 that the fourth-order NAD method works much better in suppressing the numerical dispersion effects, as indicated by the red boxes, than the traditional fourth-order finite difference method. Theory regarding suppression of dispersion effects is discussed in Lang C and Yang DH (2017).