Seismic Velocity Inverse via Physics Embedded RNN

1. Introduction

Seismic waveform inversion is an essential geophysical method that enables the inversion of physical parameters of the subsurface medium, such as velocity, density, and elastic moduli, from seismic data[1]. These parameters are of great significance in fields such as oil and gas exploration, earthquake hazard assessment, and geological structure analysis[2]. However, seismic waveform inversion is also a non-linear, multi-modal, non-convex optimization problem and faces various challenges, such as dependency on initial models, susceptibility to local minima traps, and high computational complexity[3].

To overcome these challenges, deep learning techniques have been introduced to seismic inversion. Röth and Tarantola were the first to propose the use of neural networks to transform time-domain seismic data into acoustic velocity models in 1994[4]. Nath et al., in 1999, trained neural networks with synthetic data, resulting in networks capable of estimating formation velocities using inter-well data[5,6]. In 2019, Yang and Ma developed a fully convolutional neural network for inverting P-wave velocity models from raw seismic data[7]. A commonality among these networks is the need for extensive labeled data for training, and once trained, the neural networks can be used as inversion solvers. However, in the field of seismic exploration, acquiring a large labeled dataset is challenging.

Physics-Embedded Neural Networks (PINNs) offer a novel approach to seismic velocity modeling[8]. The fundamental idea behind PINNs is to use neural networks to approximate unknown functions or parameters and incorporate the residuals of physical equations as part of the loss function. The neural networks are trained by minimizing this loss function, resulting in solutions that satisfy physical constraints[9]. Moseley et al. (2020) utilized PINNs as solvers for acoustic wave forward propagation[10,11], while Smith et al.[12] and Waheed et al.[13] applied PINNs to the Eikonal equation for first-arrival prediction and travel-time tomography, respectively. Song et al. employed PINNs to solve the frequency-domain anisotropic acoustic wave equation[14]. Additionally, PINNs have been successfully used for forward and inverse analyses of various complex systems in cases with limited sensor measurements[15,16,17,18,19,20] or without labeled data[21,22,23,24,25].

This paper utilizes Recurrent Neural Networks (RNNs) to simulate the 2D time-domain finite difference acoustic wave equation as the forward process and treats velocity as a network parameter. The velocity model is updated during the backward propagation of the network. This approach enables the integration of information from the acoustic wave equation into the network, guiding the optimization direction of the network, and avoiding explicit solutions for wavefields and gradients, thereby reducing computational costs. Furthermore, it leverages the non-linear fitting capabilities of neural networks to enhance inversion accuracy and robustness. Additionally, given the temporal causality in wavefields, where the current computed wavefield is related to previous wavefields, this paper adopts a time-lapse inversion approach. This ensures that the residuals of the wavefields are minimized at each timestep, guaranteeing the accuracy of subsequent wavefield computations and alleviating the non-uniqueness issue in velocity inversion. The effectiveness and superiority of the method proposed in this paper are verified through numerical experiments.

2. Materials and Methods

2.1. Finite-difference Acoustic Wave Equation

The numerical simulation of wave equation forward modeling is fundamental to inversion. The time-space domain finite difference algorithm is widely used in geophysical forward modeling due to its simplicity and high computational efficiency, and it is employed to simulate the propagation of seismic waves. In the process of solving partial differential wave equations using the finite difference method, the computation area is first discretized into a grid, then grid differences are utilized to approximate the derivative operators, and the wave equation is converted into a difference equation.

In a 2D isotropic acoustic medium, the time-domain wave equation is:

$$\begin{array}{c}{\nabla }^2\mathrm{u}\left(\mathrm{r},\mathrm{t}\right)=\frac{1}{{\mathrm{v}}^2\left(\mathrm{r}\right)}\cdot \frac{\partial {\mathrm{u}}^2\left(\mathrm{r},\mathrm{t}\right)}{\partial {\mathrm{t}}^2}+s\left(\mathrm{r},\mathrm{t}\right)\mathsf{\delta }\left(\mathrm{r}-{\mathrm{r}}_{\mathrm{s}}\right)\end{array}$$

(1)

where, $\nabla {\mathrm{u}}^2\left(\mathrm{r},\mathrm{t}\right)$ represents the pressure or displacement at coordinate r at time t, ${\mathrm{V}}^2\left(\mathrm{r}\right)$represents the speed of sound at coordinate r, $\mathrm{s}\left(\mathrm{r},\mathrm{t}\right)$ represents the source function at coordinate r at time t, and $\mathsf{\delta }\left(\mathrm{r}-{\mathrm{r}}_{\mathrm{s}}\right)$ is a Dirac (Delta) function, indicating that energy is released only at the source position.

For the 2D medium, spatial discretization is done by converting coordinates (x, y) into grid points (i, j), where i represents the index in the x direction, and j represents the index in the y direction. Uniform grid division is usually adopted, dividing the space into a series of equal-sized cells. Time t is discretized into time steps Δt, allowing time to progress in fixed intervals. Based on discrete spatial grids and time steps, difference approximations are used to replace derivative terms. A centered difference approximation is used to approximate the second-order spatial derivative $\nabla {\mathrm{u}}^2\left(\mathrm{r},\mathrm{t}\right)$:

$$\begin{array}{c}\nabla {\mathrm{u}}^2\left(\mathrm{i},\mathrm{j},\mathrm{t}\right)\approx \frac{\mathrm{u}\left(\mathrm{i}+1,\mathrm{j},\mathrm{t}\right)-2\mathrm{u}\left(\mathrm{i},\mathrm{j},\mathrm{t}\right)+\mathrm{u}\left(\mathrm{i}-1,\mathrm{j},\mathrm{t}\right)}{\Delta {\mathrm{x}}^2}+\frac{\mathrm{u}\left(\mathrm{i},\mathrm{j}+1,\mathrm{t}\right)-2\mathrm{u}\left(\mathrm{i},\mathrm{j},\mathrm{t}\right)+\mathrm{u}\left(\mathrm{i},\mathrm{j}-1,\mathrm{t}\right)}{\Delta {\mathrm{y}}^2}\end{array}$$

(2-2)

where Δx and Δy represent the grid spacing in the x and y directions, respectively.

A second-order centered difference approximation is used to approximate the second-order time derivative $\partial {\mathrm{u}}^2\left(\mathrm{r},\mathrm{t}\right)/\partial {\mathrm{t}}^2$:

$$\begin{array}{c}\frac{\partial {\mathrm{u}}^2\left(\mathrm{r},\mathrm{t}\right)}{\partial {\mathrm{t}}^2}\approx \frac{\mathrm{u}\left(\mathrm{i},\mathrm{j},\mathrm{t}+\Delta \mathrm{t}\right)-2\mathrm{u}\left(\mathrm{i},\mathrm{j},\mathrm{t}\right)+\mathrm{u}\left(\mathrm{i},\mathrm{j},\mathrm{t}-\Delta \mathrm{t}\right)}{\Delta {\mathrm{t}}^2}\end{array}$$

(2-3)

Substituting the difference approximations into the original wave equation yields the 2D time-domain finite difference acoustic wave equation:

$$\begin{array}{c}\mathrm{u}\left(\mathrm{i},\mathrm{j},\mathrm{t}+\Delta \mathrm{t}\right)={\mathrm{V}}^2\left(\mathrm{i},\mathrm{j}\right)\Delta {\mathrm{t}}^2\left[\nabla {\mathrm{u}}^2\left(\mathrm{i},\mathrm{j},\mathrm{t}\right)-\mathrm{s}\left(\mathrm{i},\mathrm{j},\mathrm{t}\right)\mathsf{\delta }\left(\mathrm{i},\mathrm{j}-{\mathrm{r}}_{\mathrm{s}}\right)\right]+2\mathrm{u}\left(\mathrm{i},\mathrm{j},\mathrm{t}\right)-\mathrm{u}\left(\mathrm{i},\mathrm{j},\mathrm{t}-\Delta \mathrm{t}\right)\end{array}$$

(2-4)

Through the difference acoustic wave equation, the wavefield at the current time can be calculated using the wavefields of the past two time steps, the source, and the velocity model. This relationship can be seen as a function g:

$$\begin{array}{c}\mathrm{u}\left(\mathrm{r},\mathrm{t}+\Delta \mathrm{t}\right)=\mathrm{g}\left(\mathrm{u}\left(\mathrm{r},\mathrm{t}\right),\mathrm{u}\left(\mathrm{r},\mathrm{t}-\Delta \mathrm{t}\right),\mathrm{s}\left(\mathrm{r},\mathrm{t}\right),{\mathrm{v}}^2\left(\mathrm{r}\right)\right)\end{array}$$

(2-5)

The wavefield at each moment can be considered as part of the time series data. The calculation process of the difference acoustic wave equation is shown in Figure 1 and Figure 2.

2.2. Recurrent Neural Network with Acoustic Wave Equation Embedding

The structure of the RNN (Recurrent Neural Network) is shown in Figure 3. By comparison, it can be observed that the propagation structure of the difference wave equation has high similarity with the network structure of the RNN, and existing experiments indicate that the dynamics of the wave equation are conceptually equivalent to the RNN. Therefore, the acoustic difference equation can be embedded into the RNN structure, which enables a better understanding of the propagation characteristics of waves and maps them into the neural network structure [26].

Considering the dependence of the acoustic difference equation on the wavefields of the two previous time steps, these wavefields can be used as hidden states in the RNN. The source function at the current time step can be used as input, and the signal at the detector position for each time step 't' can be used as the shot gather at time 't'. The velocity model in the difference equation can be treated as a trainable parameter within the network. By using the error between the observed and predicted values for backward correction, this process can be viewed as the inversion of the velocity model.

This approach allows the information from the acoustic wave equation to be incorporated into the network, guiding the optimization direction of the network, and avoiding the explicit calculation of gradients, which reduces computational costs. Additionally, the non-linear fitting capabilities of neural networks can be leveraged to improve inversion accuracy and robustness.

The structure of the RNN with embedded difference acoustic wave equation and the structure of operators are shown in Figure 4 and Figure 5, respectively.

2.3. Time-by-Time Inversion Algorithm

In RNNs, two common network updating algorithms are generally employed: (1) Compute the loss and update the network at each time step: At each time step, the loss is computed and used to update the network parameters. This is a traditional time-based loss update method. (2) Compute the loss at each time step, sum them, and then update the network: At each time step, the current loss is computed and accumulated with losses from previous time steps. The network parameters are then updated using the accumulated sum of losses. This method integrates information from previous time steps into the current time step loss to provide global information.

In the RNN embedded with the acoustic wave equation, updating the network by computing the loss at each time step would result in changing the velocity model at different time steps. According to the calculation formula of the difference wave equation (5), the velocity needs to remain consistent throughout the process of solving the acoustic wave equation. Therefore, computing the loss and updating the parameters separately at each time step would violate the constraints of the acoustic wave equation, leading to the inversion not converging. The article introduces a new network loss update method called the step-by-step inversion algorithm. Considering the causality of wavefield propagation, the calculation of the wavefield at the current time step is related to the previous wavefields. If errors are present in the calculation of the preceding wavefields, it will introduce additional errors in the subsequent wavefield calculations. Therefore, it is necessary to ensure the correctness of the preceding wavefield calculations before calculating subsequent wavefields. The step-by-step inversion adopts the idea of ensuring the correctness of the preceding wavefields for inversion. Specifically, starting from the beginning of wavefield propagation, the sum of the losses of the wavefield at each time step and its preceding time steps is computed to calculate the total loss. The model parameters are updated according to the total loss to ensure the correctness of the wavefield at the current time step. Compared to the traditional method of independently calculating the loss and updating the network at each time step or calculating the sum of the losses for all time steps, this method has significant advantages. Firstly, step-by-step inversion not only minimizes the overall residual but also ensures that the residual of the wavefield at each time step is as small as possible. In addition, by ensuring that the residual of the wavefield at the previous time step is small before propagating to the next time step, this method incorporates causal information of wavefield propagation, satisfying the propagation conditions of the acoustic wave equation. This avoids errors at the current time step that could be caused by errors from the previous time step, improving the accuracy of the inversion. By using a velocity model that has been partially corrected through inversion as the network parameters for the next time step, the initial model for the next time step inversion is altered, which can alleviate the dependency on the initial model in velocity inversion.

The mean squared error loss is calculated at each time step by comparing the output with the label:

$$\begin{array}{c}\mathsf{\delta }{\mathrm{d}}_{\mathrm{t}}={{(\mathrm{d}}_{\mathrm{t}}-{{\mathrm{d}}^{\prime }}_t)}^2\end{array}$$

(2-8)

where ${\mathrm{d}}_{\mathrm{t}}$ is the output of the network at time t, and ${{\mathrm{d}}^{\prime }}_t$ is the label sequence at time t. The residual at time t used for computing the gradient through backpropagation is:

$$\begin{array}{c}{\mathrm{loss}}_{\mathrm{t}}={\displaystyle \sum _{r_s}}{\displaystyle \sum _0^{\mathrm{t}}}\mathsf{\delta }{\mathrm{d}}_{\mathrm{t}}\end{array}$$

(2-9)

where $r_s$​ represents different shot gathers. The gradients are calculated for each time step to compute the update for the model:

$$\begin{array}{c}\mathsf{\delta }{\mathrm{v}}_{\mathrm{t}}\left(\mathrm{r}\right)=-\frac{\partial {\mathrm{l}\mathrm{o}\mathrm{s}\mathrm{s}}_{\mathrm{t}}}{\partial {\mathrm{v}}_{\mathrm{t}}\left(\mathrm{r}\right)}\end{array}$$

(2-10)

The gradient is propagated through the chain rule, so the velocity correction is related to ${\mathrm{d}}_{\mathrm{t}}$​, which causes the residual. Through the two-dimensional finite difference acoustic wave equation in the time domain (5), it is known that ${\mathrm{d}}_{\mathrm{t}}$ is related to the velocity model and the wavefields of the previous two time steps. The spatial domain variations at different time steps are caused by the difference in the wavefields of the previous time step. Therefore, the range of velocity correction gradually expands with the range of the wavefield. Through the step-by-step correction algorithm in the time domain, the velocity model can be gradually corrected in the spatial domain, starting from the source and progressively correcting the velocity model in all directions. The gradual correction in the spatial domain makes the inversion parameters smaller, thereby reducing the non-uniqueness of the inversion problem. Moreover, a more accurate shallow velocity can improve the accuracy and stability of deep velocity corrections. The schematic diagram of the step-by-step inversion process is shown in Figure 6.

3. Results

To validate the effectiveness of the Step-by-Step Inversion Algorithm and its advantages over the Sum-of-Losses Update method, numerical experiments were conducted using both a Homogeneous Layer Model and a Sloping Layer Model for comparison and analysis. Ultimately, the Step-by-Step Inversion Algorithm was applied to synthetic data computed using the Marmousi model to verify the feasibility of this algorithm.

3.1. Homogeneous Layer Model

In order to validate the efficacy of the Step-by-Step Inversion Algorithm, a Homogeneous Layer Model was constructed. In this model, the velocities from the shallow to deep regions are set as [3, 3.5, 4] km/s. Synthetic datasets were computed through forward propagation using a Vertical Seismic Profile (VSP) observation system. The initial model has velocities set as [3, 3.4, 3.9] km/s from shallow to deep regions, simulating a scenario where the initial model is relatively accurate. The inversion was performed using the Adam optimization algorithm, with identical parameter settings for all tests. The Homogeneous Layer Model and its inversion results are depicted in Figure 7.

As can be seen from Figure 7, the Step-by-Step Inversion Algorithm yields smoother results both in shallow and deep areas, especially near the source points close to the surface and the lateral regions. Slices were taken at depths of 500m and 800m for velocity profile comparisons, with the results shown in Figure 8 and Figure 9 respectively.

It can be deduced from Figure 8 and Figure 9 that, in terms of inversion accuracy, the Step-by-Step Inversion generally achieves higher accuracy compared to the Sum-of-Losses Inversion. The inversion results are closer to the true model. Moreover, by employing a time-stepped inversion, which takes into account the propagation of waves, the accuracy of the preceding wavefield is ensured. This, in turn, leads to more precise calculations in subsequent wavefield propagations, manifesting in the inversion results as more accurate representations on the lateral sides and deeper regions of the model.

3.2. Sloping Layer Model

To evaluate the performance of the Step-by-Step Inversion Algorithm in inclined geological structures, as well as its weak dependence on the initial model, a sloping layer model was constructed. In this second model, the velocities increase from the shallow to deeper layers, following [2.5, 3.5, 4.5] km/s respectively. The synthetic dataset was computed using a VSP (Vertical Seismic Profile) observation system via forward modeling. The initial model had velocities [2, 3, 4] km/s, with significant discrepancies between the initial and true models. This was to validate the weak dependency of the Step-by-Step Inversion Algorithm on the initial model, demonstrating that accurate inversion is still achievable despite large errors in the initial model. The inversion employed the Adam optimization algorithm with identical parameter settings. The results of the second model and its inversion are shown in Figure 10.

As can be observed in Figure 10, despite the substantial error in the initial model in the context of a sloping layer model, the Step-by-Step Inversion Algorithm exhibited superior performance overall. This suggests that the algorithm has a lower dependency on the initial model and is capable of achieving reasonably accurate inversions even when the initial model has large discrepancies. In contrast, the traditional Sum-of-Losses Inversion performed poorly, indicating a higher dependency on the initial model. Velocity profiles at depths of 150m, 500m, and 800m were selected for comparison and are displayed in Figure 11, Figure 12, and Figure 13 respectively.

Based on the cross-sectional velocity comparison plots, it can be observed that, even with significant errors in the initial model, the Step-by-Step Inversion Algorithm was able to more accurately invert velocities at various depths, including 150m, 500m, and 800m, demonstrating superior precision in the results.

3.3. Marmousi Model

In order to validate the application effect of the algorithm under more complex actual conditions, the Marmousi model was used for experiments. A relatively accurate and smooth velocity model was used as the initial model for inversion. The real model and initial model are shown in Figure 13. Synthetic data were generated from the real model using the VSP observation system, with the RNN parameters set consistently with the inversion parameters.

The Marmousi model and inversion results are shown in Figure 14. As can be seen from the inversion results, both the overall inversion and time-by-time inversion results are relatively close to the real model. A comparison of horizontal velocity values at depths of 500m and 850m is shown in Figure 15, and a comparison of depth velocity values at a horizontal distance of 1000m is shown in Figure 16. From the velocity comparison curve, it can be seen that both the time-by-time inversion and the overall inversion are relatively close to the real model. Moreover, the velocity of the time-by-time inversion results is smoother, with fewer velocity anomalies, making the inversion more stable.

Figure 14. Marmousi model and initial model. (a)True model. (b)Initial model.

Figure 14. Marmousi model and initial model. (a)True model. (b)Initial model.

Figure 15. Marmousi model and inversion results. (a)True model. (b) Total loss inversion (c) Time by time inversion.

Figure 15. Marmousi model and inversion results. (a)True model. (b) Total loss inversion (c) Time by time inversion.

Figure 16. Comparison of 500m deep velocity profiles.

Figure 16. Comparison of 500m deep velocity profiles.

Figure 17. Comparison of 850m deep velocity profiles.

Figure 17. Comparison of 850m deep velocity profiles.

Figure 18. Comparison of 1000m distance velocity profiles.

Figure 18. Comparison of 1000m distance velocity profiles.

4. Discussion

This paper introduced a Step-by-Step Inversion Algorithm based on a physics-embedded Recurrent Neural Network (RNN). Through the Step-by-Step Inversion Algorithm, the causal information in the temporal evolution of wave equations is mapped to spatial information for velocity model correction. The velocity model is corrected in spatial regions in accordance with the wavefield propagation range. The corrected model then serves as the initial model for subsequent time steps, covering a larger temporal range, which reduces the inversion's dependency on the initial model. Additionally, controlling the correction area results in fewer inversion parameters, effectively reducing the non-uniqueness of inversion solutions within a single time step, and ultimately decreasing the overall non-uniqueness while enhancing the accuracy and stability of the inversion.

5. Patents

This section is not mandatory but may be added if there are patents resulting from the work reported in this manuscript.