Physics-informed neural network for seismic wave inversion in layered
  semi-infinite domain
Estimating the material distribution of Earth's subsurface is a challenging
task in seismology and earthquake engineering. The recent development of
physics-informed neural network (PINN) has shed new light on seismic inversion.
In this paper, we present a PINN framework for seismic wave inversion in
layered (1D) semi-infinite domain. The absorbing boundary condition is
incorporated into the network as a soft regularizer for avoiding excessive
computation. In specific, we design a lightweight network to learn the unknown
material distribution and a deep neural network to approximate solution
variables. The entire network is end-to-end and constrained by both sparse
measurement data and the underlying physical laws (i.e., governing equations
and initial/boundary conditions). Various experiments have been conducted to
validate the effectiveness of our proposed approach for inverse modeling of
seismic wave propagation in 1D semi-infinite domain.
Pu Ren, Chengping Rao, Hao Sun, Yang Liu
# Physics-informed neural network for seismic wave inversion in layered semi-infinite domain

###### Abstract

Estimating the material distribution of Earth's subsurface is a challenging task in seismology and earthquake engineering. The recent development of physics-informed neural network (PINN) has shed new light on seismic inversion. In this paper, we present a PINN framework for seismic wave inversion in layered (1D) semi-infinite domain. The absorbing boundary condition is incorporated into the network as a soft regularizer for avoiding excessive computation. In specific, we design a lightweight network to learn the unknown material distribution and a deep neural network to approximate solution variables. The entire network is end-to-end and constrained by both sparse measurement data and the underlying physical laws (i.e., governing equations and initial/boundary conditions). Various experiments have been conducted to validate the effectiveness of our proposed approach for inverse modeling of seismic wave propagation in 1D semi-infinite domain.

keywords: Physics-informed neural network, Seismic wave propagation, Absorbing boundary conditions, Full waveform inversion, Material distribution +
Footnote †: journal: Journal name

## 1 Introduction

Seismic imaging is essential to a wide range of scientific and engineering tasks. Thanks to the development of a theoretical foundation for wave propagation, researchers have been dedicated to tackling seismic inversion problems by combining physical simulation and field measurement data. One of the most important applications is full waveform inversion (FWI), which determines the high-resolution and high-fidelity velocity model (or material distribution) of the subsurface from the seismic waveforms collected on the field surface. This task is usually solved as an inverse problem of minimizing the misfit error between the observed seismic signals and the synthetic solution simulated based on an estimated velocity model. However, the traditional numerical methods suffer from three significant limitations for inverse analysis (e.g., FWI). First, there is a harsh requirement of computational resources and memory for calculating the Frchet derivatives in the optimization procedure. Second, it limits the resolution of the inverted velocity model by using specific smoothing regularization approaches for alleviating the ill-posedness of inverse problems. The last one is that the final solution heavily depends on the initial guess of the model.

Due to the great advancement in machine learning (ML) and artificial intelligence (AI), scientists have embraced a new research direction for forward and inverse analysis of complex systems in the domain of seismology. For instance, deep neural networks (DNNs) and auto differentiation
are utilized for solving the Eikonal equation [1] and seismic inversion [2; 3]. [4] also explore the potential of Fourier Neural Operator [5] for acoustic wave propagation and inversion. In addition to inferring velocity structures, [6] develop a generalized Expectation-Maximization (DeepGEM) algorithm for estimating the wave sources simultaneously. Nevertheless, these methods rely on large amounts of high-fidelity data, which are usually inaccessible in real-world applications. Therefore, there is an urgent need to develop novel ML approaches to address the sparsity and noise issues of measurement data. It is worthwhile to mention that Graph Network (GN) [7] has been proven to be an effective tool for estimating the underlying velocity model and the corresponding initial condition of the specific wave propagation based on sparse measurement data [8].

Moreover, the recent physics-informed neural networks (PINNs) [9; 10] have also succeeded in forward and inverse analysis of various complex systems with limited sensor measurements [11; 12; 13; 14; 15] or even no labeled data [16; 17; 18; 19; 20]. The general principle of PINNs is to combine the known physical laws and the measured data to optimize the DNNs simultaneously. In particular, we observe a growing interest in applying PINNs for seismic wave modeling [21; 22; 23; 24], which can overcome the obstacles of traditional methods. Note that the current PINN implementations for seismic wave propagation primarily focus on studying acoustic wave propagation.

In this paper, we aim to design PINN frameworks for elastic wave modeling and inversion in 1D domain, which is inspired by the previous work [25]. To be more concrete, we first design a specific PINN architecture to approximate the solution variables of interest. Next, for the sake of seismic inversion, we introduce a new part of shallow NN to learn the material distribution (i.e., Young's modulus). The synthesized networks work simultaneously to model the specific seismic wave propagation. Furthermore, the absorbing boundary condition is considered here to truncate the domain and avoid unnecessary computation, which works as a soft regularizer in the loss function to constrain the entire network. Finally, we evaluate the performance of our proposed PINN frameworks on various numerical experiments, including solving wave equations and FWI. More specifically, we consider different material configurations when validating the method.

The rest of the paper is organized as follows. The Methods Section formulates the seismic inversion task (i.e., FWI) and introduces the proposed PINN framework. The Experiments Section presents various numerical experiments for validating the performance of our proposed method on seismic wave inversion. Moreover, we conclude the entire paper in the Conclusions Section.

## 2 Method

### Problem setup

In this paper, we focus on the inverse analysis of 1D seismic wave propagation in elastic media. First of all, the momentum and constitutive equations of elasticity in 1D are defined as

\[\rho\frac{\partial u^{2}}{\partial t^{2}} =\frac{\partial\sigma}{\partial x}, \tag{1}\] \[E\frac{\partial u}{\partial x} =\sigma,\]

where \(t\) and \(x\) are the time and spatial dimensions; \(\sigma(t,x)\) and \(u(t,x)\) are the stress and displacement fields, respectively. \(E(x)\) denotes Young's modulus with respect to \(x\). \(\rho\) represents the density. The initial condition of the displacement field is set as \(u(0,:)=0\). In addition, we apply the absorbing boundary condition for the truncated boundary, which is given by

\[\frac{\partial u}{\partial t}+\sqrt{\frac{E}{\rho}}\cdot\frac{\partial u}{ \partial x}=0. \tag{2}\]
Our specific goal is to design PINN frameworks to take the spatiotemporal coordinate \((x,t)\) as inputs and approximate the displacement and stress variables \(\{u,\sigma\}\) based on Eq. (1) and Eq. (2). For the FWI problem, we aim to reconstruct the velocity model (in the form of Young's modulus) of the subsurface given the displacement waveform collected on the surface (i.e., \(x=0\)).

### PINN framework

In this part, we introduce the strategy of combining DNNs and physical laws (i.e., governing equations and boundary conditions) for learning seismic wave propagation. Various network architectures have been explored in the area of scientific ML, such as fully-connected NNs [9], convolutional NNs [19], graph NNs [26; 27], and transformer [28]. Considering the complexity of our task, we employ the fully-connected feed-forward NNs to represent the seismic dynamics in this paper. The common implementation of fully-connected NNs utilizes an input layer, multiple hidden layers, and an output layer. In specific, the connection between \((i-1)\)-th layer and \(i\)-th layer is given by

\[\mathbf{X}^{i}=\sigma\left(\mathbf{W}^{i}\mathbf{X}^{i-1}+\mathbf{b}^{i} \right),\;\;i\in[1,n+1], \tag{3}\]

where \(\mathbf{X}^{i-1}\) and \(\mathbf{X}^{i}\) are input and output features at \(i\)-th layer; \(\mathbf{W}^{i}\) and \(\mathbf{b}^{i}\) denote the weight tensor and bias vector of \(i\)-th layer; \(n\) refers to the number of hidden layers; \(\sigma(\cdot)\) represents the nonlinear activation function (e.g., Hyperbolic Tangent Function).

The major difference between pure data-driven ML and physics-informed learning lies in the introduction of physical constraints apart from the traditional data loss. The physics constraints are achieved by constructing the residual loss between the predicted dynamics and the underlying true physics. In PINNs, automatic differentiation [29] is applied to obtain the derivative terms of interest. More details about the specific design of PINN architecture for seismic inversion can be found in Section 2.3.

### Seismic inversion

The goal of seismic inversion is to infer the high-fidelity velocity model (i.e., material distribution) of the subsurface based on the observed seismic signals (e.g., seismic waveforms), which are measured from the field surface. The challenge here is that the distribution of seismic wavespeed in the medium is unknown. To this end, we introduce another shallow NN to approximate the material variable \(\widetilde{E}\) independently apart from the PINN architecture designed for approximating the solution variables \(\widehat{u}\) and \(\widetilde{\sigma}\). The overall framework of our proposed method is presented in Fig. 1. Note that the shallow NN part, which simulates the material distribution, holds a smaller network size than the deep NN module due to less spatiotemporal complexity [14; 24] compared with learning solution variables. Herein, this new synergetic PINN framework serves as a model pipeline for FWI.

In addition, to avoid unnecessary computation, we consider applying the absorbing boundary condition for handling the semi-infinite domain in seismic wave propagation. The truncated boundary condition is defined as Eq. (2). Finally, the loss function for inferring the velocity model is defined as

\[J=J_{g}+\alpha\cdot J_{d}+\beta\cdot J_{ibc} \tag{4}\]

where \(\alpha\) and \(\beta\) denote weighting coefficients. \(J_{g}\), \(J_{ibc}\), and \(J_{d}\) denote the governing equation loss, initial/boundary condition loss, and data loss, respectively. The loss function guarantees the optimization of both the misfit between the predicted and the measured waveform, as well as the physical consistency described by Eq. (1). Moreover, as shown in Fig. 1, Young's modulus \(\widetilde{E}\) is approximated by the shallow NNs part, which is constrained by \(J_{g}\) and \(J_{ibc}\).
## 3 Experiments

To validate the effectiveness of our proposed method, we conduct two numerical experiments on evaluating the performance of FWI in 1D domain. Specifically, two velocity models are designed: (1) a mixed velocity model with both linear and Gaussian functions; (2) a four-layer material distribution. All of the experiments are coded with TensorFlow [30] and implemented on an NVIDIA Tesla V100 GPU card (32G) in a standard workstation.

### Setup

In these numerical examples, we consider a 1D subsurface under the elasticity assumption, as shown in Fig. (2). To ensure the computational efficiency of PINN, we truncate the original computational semi-infinite domain with a finite depth, i.e., \(x=2.0\). As a result, the absorbing boundary condition would be applied to the truncated boundary to avoid wave reflection. Note that the waveform we used for inversion is synthesized through finite element simulation in an enlarged domain \(x\in[0,10]\) though a truncated domain (i.e., \(x\in[0,2]\)) is considered in PINN.

Furthermore, we also define an evaluation metric, the mean absolute relative error (MARE), to measure the prediction discrepancy of material distribution in the form of Young's modulus \(E\). MARE is given by

\[\text{MARE}=\frac{1}{N}\sum_{i=1}^{N}\frac{|\widetilde{E}_{i}-E_{i}|}{|E_{i}| }\times 100\%, \tag{5}\]

where \(N\) is the number of collocation points.

Figure 1: Architecture of the PINN for 1D full waveform inversion. A shallow neural network is used to approximate the unknown distribution of \(E\).
### Linear and Gaussian material distribution

The first case considers a velocity model of the subsurface with both linear and Gaussian functions, which is defined in the form of Young's modulus

\[E=2.0+0.9x+4.0\exp\left[-\frac{(x-1.1)^{2}}{0.2}\right]. \tag{6}\]

Moreover, in the simulation, we prescribe a time-varying force on the surface, which is given by

\[\sigma(t)=\sigma_{0}\cdot\exp\left[-\frac{(t-t_{s})^{2}}{0.03}\right], \tag{7}\]

where the amplitude \(\sigma_{0}\) is set as 1.0 and the offset \(t_{s}\) is 0.5. The entire simulation duration is 3.

We design the deep NNs of \(50\times 3\) (i.e., depth of 3 and width of 50) for approximating the solution variables \(u\) and \(\sigma\). The shallow NN has a network size of \(4\times 2\) (i.e., depth of 2 and width of 4) for learning the material distribution \(E\). In addition, the Hyperbolic Tangent activation function (i.e., Tanh) is utilized in the framework. To train the composite PINN jointly, we adopt a routine with 2,000 iterations of Adam optimizer [31] followed by 40,000 iterations of L-BFGS-B optimizer [32].

In addition, 30,000 collocation points are randomly sampled within the spatiotemporal domain for constructing the \(J_{g}\). The synthetic wave signal includes 301 data points collected in \(t\in[0,3]\). When the optimization convergences, we can directly infer the material distribution \(E\) and predict the waveform on the surface. The comparison between the predicted velocity model and the ground truth is shown in Fig. 3. It is obvious that the learned velocity model agrees well with the ground truth despite the minor discrepancy near the truncated boundary. In specific, the evaluation metric MARE shows that the prediction error of material distribution is only 0.47%. Furthermore, we present the comparison between the predicted surface waveform (i.e., displacement responses) and the synthetic measurements in Fig. 4. The result also proves an excellent performance of our proposed PINN framework in approximating the solution variables.

Figure 2: Computational domain of the 1D full waveform inversion problem. \(F(t)\) denotes the time-varying load excited on the surface. The truncated depth is \(L=2.0\).
### Four-layer material distribution

The second numerical example is a four-layer velocity model, which is a more complicated subsurface compared with the first case. The mathematical formulation is given by

\[E=\begin{cases}1.5,&0\leq x\leq 0.5\\ 6.0,&0.5\leq x\leq 1.0\\ 4.0,&1.0\leq x\leq 1.5\\ 6.5,&\text{others}\end{cases}. \tag{8}\]

To avoid the discontinuity of Young's modulus \(E\), we smooth the transition between two adjacent layers with the Sine function. Note that the derivatives of \(E\) are still discontinuous. The transition length is set to \(0.1\). This problem becomes more complicated since the distribution of \(E\) would dissipate the seismic wave during propagation and lead to a faint waveform. Therefore, the excitation of surface load applied here features a higher frequency which is described by

\[\sigma(t)=\sigma_{0}\cdot\exp\left[-\frac{(t-t_{s})^{2}}{0.008}\right], \tag{9}\]

Figure 4: Comparison between the learned surface waveform (i.e., displacement) and the ground truth.

Figure 3: Comparison between the predicted material distribution in the form of Young’s modulus \(E\) and the true velocity model.
\(\sigma_{0}\), \(t_{s}\), and time duration are still defined as 1.0, 0.5, and 3, respectively.

The NNs with more powerful expressiveness are employed to approximate the solutions, i.e., deep NN with a network size of \(80\times 4\) and shallow NN with the size of \(10\times 2\) for learning \(E\). Additionally, 100,000 collocation points are sampled within the spatiotemporal domain. 5,000 iterations of Adam optimizer and 400,000 iterations of L-BFGS-B optimizer are used to train the composite PINNs. The inverted distribution of \(E\) is inferred from the trained model for comparison with the ground truth, as shown in Fig. 5. It can be seen that the PINN is capable of identifying the subsurface with four layers. The prediction of four-layer material distribution has a MARE of 1.70%, which also exhibits great performance for the complex velocity model. Nevertheless, we also observe that the discrepancy increases as the depth increases. This is owing to the inherent nature of this problem - as the wave propagates deep, the reflected wave that is observable on the surface carries fewer features. It reflects the poor identifiability of the deeper subsurface and the ill-posedness of such tasks. Moreover, the discontinuity (i.e., derivatives of \(E\)) near the interface of adjacent layers also makes it challenging to accurately reconstruct the subsurface material as the PINN is a continuous approximator. Besides, the comparison of the surface waveform between the prediction from PINN and the ground truth also validates the effectiveness of our proposed method. To summarize, our designed PINN framework shows promising results on seismic inversion problems.

## 4 Conclusions

This paper introduces a novel PINN framework for learning seismic wave propagation in elastic media. In specific, a deep NN is designed to approximate the physical laws (governing equations and the initial/boundary condition), and a shallow NN is employed to learn the underlying velocity structure. Moreover, the absorbing boundary is incorporated into the network to truncate the semi-infinite domain and avoid excessive computation. We validate the effectiveness of our proposed approach in two numerical cases in the context of seismic inversion (i.e., FWI). In the future, we would like to explore the potential of physics-informed learning for high-dimensional seismic inversion problems (e.g., in 2D/3D domains) but are not limited to employing fully-connected NNs that are typically used in PINNs. This is due to the observation that the physics-informed discrete learning scheme [15; 18; 19] is promising to estimate more complex velocity models.

Figure 5: Comparison between the predicted material distribution in the form of Young’s modulus \(E\) and the true velocity model.
## Acknowledgement

P. Ren would like to gratefully acknowledge the Vilas Mujumdar Fellowship at Northeastern University. Y. Liu would like to thank the support from the Fundamental Research Funds for the Central Universities. The work is also supported by the National Natural Science Foundation of China (No. 62276269).
