On the Use of Neural Networks for Full Waveform Inversion
Neural networks have recently gained attention in solving inverse problems.
One prominent methodology are Physics-Informed Neural Networks (PINNs) which
can solve both forward and inverse problems. In the paper at hand, full
waveform inversion is the considered inverse problem. The performance of PINNs
is compared against classical adjoint optimization, focusing on three key
aspects: the forward-solver, the neural network Ansatz for the inverse field,
and the sensitivity computation for the gradient-based minimization. Starting
from PINNs, each of these key aspects is adapted individually until the
classical adjoint optimization emerges. It is shown that it is beneficial to
use the neural network only for the discretization of the unknown material
field, where the neural network produces reconstructions without oscillatory
artifacts as typically encountered in classical full waveform inversion
approaches. Due to this finding, a hybrid approach is proposed. It exploits
both the efficient gradient computation with the continuous adjoint method as
well as the neural network Ansatz for the unknown material field. This new
hybrid approach outperforms Physics-Informed Neural Networks and the classical
adjoint optimization in settings of two and three-dimensional examples.
Leon Herrmann, Tim Bürchner, Felix Dietrich, Stefan Kollmannsberger
# On the Use of Neural Networks for Full Waveform Inversion

###### Abstract

Neural networks have recently gained attention in solving inverse problems. One prominent methodology are Physics-Informed Neural Networks (PINNs) which can solve both forward and inverse problems. In the paper at hand, full waveform inversion is the considered inverse problem. The performance of PINNs is compared against classical adjoint optimization, focusing on three key aspects: the forward-solver, the neural network Ansatz for the inverse field, and the sensitivity computation for the gradient-based minimization. Starting from PINNs, each of these key aspects is adapted individually until the classical adjoint optimization emerges. It is shown that it is beneficial to use the neural network only for the discretization of the unknown material field, where the neural network produces reconstructions without oscillatory artifacts as typically encountered in classical full waveform inversion approaches. Due to this finding, a hybrid approach is proposed. It exploits both the efficient gradient computation with the continuous adjoint method as well as the neural network Ansatz for the unknown material field. This new hybrid approach outperforms Physics-Informed Neural Networks and the classical adjoint optimization in settings of two and three-dimensional examples.

keywords: physics-informed neural networks, full waveform inversion, deep learning, adjoint state method, inverse problems, finite difference method +
Footnote †: journal: _Engineering with Computers_.
## List of Symbols

### Motivation

Since the emergence of Physics-Informed Neural Networks (PINNs) [1, 2, 3], they have been of great interest to the scientific community of computational mechanics. Despite their novelty, multiple review articles have already emerged [4, 5, 6]. The essence of PINNs is a variational approach, where a solution field is approximated with a neural network whose weights and biases are learned by minimizing the residual of a differential equation evaluated at collocation points. This general formulation allows for solving a large variety of problems that can at least partially be described by differential equations. At the same time, current machine learning frameworks like PyTorch [7] and TensorFlow [8] provide generic tools to construct PINNs with very little additional and rather simple code compared to the complexity of classical differential equation solvers.

The method shows potential in settings where the amount of data available is insufficient to build supervised surrogate models, as e.g. presented in [9, 10, 11, 12, 13, 14, 15]. If used as a forward solver, PINNs are often computationally more expensive for many common differential equations which are efficiently solvable with established forward solvers [16, 17], as e.g. finite element [18] or finite difference methods [19]. However, in settings where classical methods show difficulties, PINNs are promising. Possible applications are high-dimensional differential equations [20, 21], surrogate modeling [22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32], accelerating traditional solvers [16], and inverse problems [25, 26, 27, 28, 29, 30, 31, 32].

The novelty in this paper is twofold. Firstly, it pinpoints the strengths and weaknesses of PINNs for inverse problems in comparison to classical full waveform inversion where the continuous adjoint state method [33, 34] is considered as the state-of-the-art reference method. To this end, a thorough investigation of the strengths and weaknesses of the individual components of each method is carried out by going from PINNs to the adjoint optimization step by step. The investigation is conducted on the scalar wave equation, popular in non-destructive testing [35, 36], and recently explored with PINNs in [25, 26]. In particular, the following three aspects will be examined extensively
1. the use of PINNs as a forward-solver,
2. neural networks as an Ansatz for the coefficient field of the PDE,
3. gradient computation.

This dissection gives rise to the methods summarized in table 1 which are described and studied in detail in section 2 apart from the more general setting of PINNs described in the introductory section 1.3.

Most of these methods naturally emerge as individual steps moving from a fully PINN-based approach to a classical adjoint optimization. By contrast, and as a second novelty, the paper at hand presents a hybrid method, which combines the strengths of the investigated approaches and delivers more accurate results.

### Full Waveform Inversion

Full waveform inversion aims to non-destructively estimate the material distribution in a specimen. A conceptual illustration is provided in figure 1 with a rectangular domain \(\Omega\) and a circular void \(\Omega_{V}\), that is to be detected. Voids are identified using a fixed number of sources emitting a signal (red circles) and sensors (black circles) recording the image of the signal.

The physics of this problem can be modeled with the scalar wave equation in an isotropic heterogeneous medium, which for the solution field \(u(\mathbf{x},t)\), the density \(\rho(\mathbf{x})\), the wave speed \(c(\mathbf{x})\) and

\begin{table}
\begin{tabular}{l l l} \hline
**Method** & **Ansatz** & **Gradient Computation** \\ \hline PINNs & \(\hat{\gamma}(\mathbf{x})=A_{\gamma}(\mathbf{\theta}_{\mathbf{\gamma}};\mathbf{x})\) & \\ (section 1.3) & \(\hat{u}(\mathbf{x},t)=A_{u}(\mathbf{\theta}_{\mathbf{u}};\mathbf{x},t)\) & \(\nabla_{\mathbf{\theta}_{\mathbf{\gamma}}}\mathcal{L}\) with automatic differentiation \\ \hline PINNs with non-trainable & \(\hat{\gamma}(\mathbf{x})=A_{\gamma}(\mathbf{\theta}_{\mathbf{\gamma}};\mathbf{x})\) & \\ forward operator & \(\hat{u}(\hat{\gamma},\mathbf{x},t)=F[\hat{\gamma},\mathbf{x},t]\) & \(\nabla_{\mathbf{\theta}_{\mathbf{\gamma}}}\mathcal{L}\) with automatic differentiation \\ (section 2.1) & \(\hat{u}(\hat{\gamma},\mathbf{x},t)=F[\hat{\gamma},\mathbf{x},t]\) & \(\nabla_{\mathbf{\theta}_{\mathbf{\gamma}}}\mathcal{L}\) with automatic differentiation \\ \hline Adjoint Method & \(\hat{\gamma}(\mathbf{x})=C_{\gamma}(\mathbf{c},\mathbf{x})\) & \(\nabla_{\mathbf{c}}\mathcal{L}\) with automatic differentiation \\ (section 2.2) & \(\hat{u}(\hat{\gamma},\mathbf{x},t)=F[\hat{\gamma},\mathbf{x},t]\) & or \(\nabla_{\mathbf{c}}\mathcal{L}\) with the adjoint method \\ \hline \hline Hybrid Approach & \(\hat{\gamma}(\mathbf{x})=A_{\gamma}(\mathbf{\theta}_{\mathbf{\gamma}},\mathbf{x})\) & \(\nabla_{\hat{\gamma}}\mathcal{L}\) with the adjoint method \\ (section 2.3) & \(\hat{u}(\hat{\gamma},\mathbf{x},t)=F[\hat{\gamma},\mathbf{x},t]\) & and \(\nabla_{\mathbf{\theta}_{\mathbf{\gamma}}}\hat{\gamma}\) with automatic differentiation \\ \hline \end{tabular}
\end{table}
Table 1: Investigated methods.
the volume force \(f(\mathbf{x},t)\) is given as

\[\rho(\mathbf{x})u_{tt}(\mathbf{x},t)-\nabla\cdot(\rho(\mathbf{x})c(\mathbf{x})^{2}\nabla u(\mathbf{x},t))-f(\mathbf{x},t)=0\qquad\text{on }\Omega\times\mathcal{T}. \tag{1}\]

Additionally, homogeneous Neumann boundary conditions and homogeneous initial conditions are prescribed:

\[u_{x}(\mathbf{x},t) =0 \text{on }\Gamma_{y}\times\mathcal{T} \tag{2}\] \[u_{y}(\mathbf{x},t) =0 \text{on }\Gamma_{x}\times\mathcal{T}\] (3) \[u(\mathbf{x},0) =u_{t}(\mathbf{x},0)=0 \text{on }\Omega \tag{4}\]

Parametrizing the inverse problem with an indicator function \(\gamma\) by scaling the density leads to the most accurate reconstruction of voids and unknown homogeneous Neumann boundaries. The constant intact material has the density \(\rho_{0}\) and \(\gamma\) is constrained between \([\epsilon,1]\), where \(\epsilon\ll 1\) is a non-zero lower bound:

\[\rho(\mathbf{x})=\gamma(\mathbf{x})\rho_{0} \tag{5}\]

In contrast, the wave speed is assumed to be constant \(c(\mathbf{x})=c_{0}\). With this parametrization, the scalar wave equation takes the form:

\[\gamma(\mathbf{x})\rho_{0}u_{tt}(\mathbf{x},t)-\nabla\cdot(\gamma(\mathbf{x})\rho_{0}c_{0 }^{2}\nabla u(\mathbf{x},t))-f(\mathbf{x},t)=0\qquad\text{on }\Omega\times\mathcal{T}. \tag{6}\]

The volume force \(f(\mathbf{x},t)\) emits the signal at the position \(\mathbf{x}_{s}\). This is usually modeled by a spatial Dirac delta using a source term \(\psi(t)\).

\[f(\mathbf{x},t)=\psi(t)\delta(\mathbf{x}-\mathbf{x}_{s}) \tag{7}\]

A typical source in the context of non-destructive testing is the sine burst with \(n_{c}\) cycles, a frequency \(\omega=2\pi f_{\psi}\), and an amplitude \(\psi_{0}\), so that

\[\psi(t)=\begin{cases}\psi_{0}\sin(\omega t)\sin(\frac{\omega t}{2n_{c}}),& \text{for }0\leq t\leq\frac{2\pi n_{c}}{\omega}\\ 0,&\text{for }\frac{2\pi n_{c}}{\omega}<t.\end{cases} \tag{8}\]

This burst type is used in all examples in this paper.

The goal of the inversion is to reconstruct a material distribution, i.e. the indicator \(\hat{\gamma}(\mathbf{x};\mathbf{m})\) parametrized by \(\mathbf{m}\), that best reproduces the \(N_{\mathcal{M}}\) sensor measurements \(u_{\mathcal{M}}(\mathbf{x}_{i},t)\) during a forward simulation which yields \(\hat{u}(\mathbf{x},t;\mathbf{m})\). The inversion is posed as an optimization problem, where the cost function is the mean squared error between the prediction \(\hat{u}\) and the measurements \(u_{\mathcal{M}}\):

\[\mathcal{L}_{\mathcal{M}}(\hat{\mathbf{u}},\mathbf{u}_{\mathcal{M}})=\frac{1}{2}\int_ {\mathcal{T}}\int_{\Omega}\sum_{i=1}^{N_{\mathcal{M}}}\biggl{(}\hat{u}(\mathbf{x},t;\mathbf{m})-u_{\mathcal{M}}(\mathbf{x}_{i},t)\biggr{)}^{2}\delta(\mathbf{x}-\mathbf{x}_{i} )d\Omega d\mathcal{T}. \tag{9}\]

The gradient of \(\mathcal{L}_{\mathcal{M}}\) w.r.t. the parameters \(\mathbf{m}\) is expressed by \(\nabla_{m}\mathcal{L}_{\mathcal{M}}(\hat{\mathbf{u}}(\mathbf{m}),\mathbf{u}_{\mathcal{M}})\) and computed with the adjoint state method. The gradient is then used in a gradient-based optimization procedure to minimize \(\mathcal{L}_{\mathcal{M}}\) by iteratively updating the parameters \(\mathbf{m}\):

\[\mathbf{m}^{(j+1)}=\mathbf{m}^{(j)}-\alpha\nabla_{m^{(j)}}\mathcal{L}_{\mathcal{M}} \tag{10}\]

in which \(\alpha\) is the step size and the superindex \((j)=1,2,\dots\) is the iteration counter. This iterative process then yields an optimized estimate of the indicator \(\hat{\gamma}(\mathbf{x};\mathbf{m})\). In case of multiple sources, an average of the mean squared errors \(\mathcal{L}_{\mathcal{M}}\) and their gradients \(\nabla_{m}\mathcal{L}_{\mathcal{M}}\) must be considered in the iteration update from equation (10). More details on the computation of the gradient using the adjoint method are found in [37]. For a more in depth explanation of the adjoint method for inverse problems, see Givoli [33].
### Physics-Informed Neural Network

The idea of PINNs is to use an artificial neural network \(A\) to approximate either the solution \(\hat{u}\) or and/or the coefficients of the differential equation \(\hat{\gamma}\) using a collocation type approach. In the paper at hand, equation (6) is considered where the solution \(\hat{u}\) is the wave field and \(\rho_{0}\) is scaled by \(\gamma(\mathbf{x})\) as given in equation (5). In the most general case, both \(u\) and \(\gamma\) are approximated by a neural network \(A=(A_{u},A_{\gamma})\):

\[\hat{u} = A_{u}(\mathbf{x},t;\mathbf{\theta}_{u}) \tag{11}\] \[\hat{\gamma} = A_{\gamma}(\mathbf{x},t;\mathbf{\theta}_{\gamma}) \tag{12}\]

where \(\mathbf{\theta}=\{\mathbf{\theta}_{u},\mathbf{\theta}_{\gamma}\}\) are the model parameters of the neural network and the hat symbol indicates that the values of \(u\) and \(\gamma\) are predictions from a neural network. In the remaining paper, in models using only one of the two neural networks, i.e. either \(A_{u}\) or \(A_{\gamma}\), the notation \(\mathbf{\theta}\) without subindex is used as the corresponding neural network parameters.

The initial boundary value problem defined by equations (2) to (4) and (6) is then solved by an optimization procedure which is similar to the one used in the classical full waveform inversion described in the previous section. To this end, the initial boundary problem is first cast into its residual form:

\[\mathcal{N}[u;\gamma;\mathbf{x},t] = \gamma(\mathbf{x})\rho_{0}u_{tt}(\mathbf{x},t)-\nabla\cdot(\gamma(\mathbf{x}) \rho_{0}c_{0}^{2}\nabla u(\mathbf{x},t))-f(\mathbf{x},t) \tag{13}\] \[\mathcal{B}_{x}[u;\mathbf{x},t] = u_{x}(\mathbf{x},t)\] (14) \[\mathcal{B}_{y}[u;\mathbf{x},t] = u_{y}(\mathbf{x},t)\] (15) \[\mathcal{I}_{0}[u;\mathbf{x},t] = u(\mathbf{x},t)\] (16) \[\mathcal{I}_{1}[u;\mathbf{x},t] = u_{t}(\mathbf{x},t). \tag{17}\]

The operators \(\mathcal{N},\mathcal{B}_{x},\mathcal{B}_{y},\mathcal{I}_{0},\mathcal{I}_{1}\) are the residuals of the corresponding equations (2) to (4) and (6). These residuals are evaluated at a number of collocation points \(N=N_{\mathcal{N}}+N_{\mathcal{B}_{x}}+N_{\mathcal{B}_{y}}+N_{\mathcal{I}_{0}} +N_{\mathcal{I}_{1}}\) using the current approximation of the solution coefficient fields given by equations (11) and (12). This requires the computation of the first and second derivatives of equations (11) and (12) w.r.t \(\mathbf{x}\) and \(t\) which is carried out by either using reverse-mode automatic differentiation [38] or numerical differentiation [10, 39].

The loss functional to be minimized by the overall optimization procedure, i.e. the cost function \(\mathcal{L}\), is then defined as the weighted sum of squares of the residual evaluations (the differential equation residual \(\mathcal{L}_{\mathcal{N}}\), the boundary condition residual \(\mathcal{L}_{\mathcal{B}}\), and the initial condition residual \(\mathcal{L}_{\mathcal{I}}\)) at the collocation points:

\[\mathcal{L}_{\mathcal{N}}(\mathbf{u},\mathbf{\gamma}) = \frac{1}{2N_{\mathcal{N}}}\sum_{i=1}^{N_{\mathcal{N}}}\lambda_{ \mathcal{N},i}(\mathcal{N}[u_{i};\gamma_{i}])^{2} \tag{18}\] \[\mathcal{L}_{\mathcal{B}}(\mathbf{u}) = \frac{1}{2N_{\mathcal{B}_{x}}}\sum_{i=1}^{N_{\mathcal{B}_{x}}} \lambda_{\mathcal{B}_{x},i}(\mathcal{B}_{x}[u_{i}])^{2}+\frac{1}{2N_{\mathcal{ B}_{y}}}\sum_{i=1}^{N_{\mathcal{B}_{y}}}\lambda_{\mathcal{B}_{y},i}( \mathcal{B}_{y}[u_{i}])^{2}\] (19) \[\mathcal{L}_{\mathcal{I}}(\mathbf{u}) = \frac{1}{2N_{\mathcal{I}_{0}}}\sum_{i=1}^{N_{\mathcal{I}_{0}}} \lambda_{\mathcal{I}_{0},i}(\mathcal{I}_{0}[u_{i}])^{2}+\frac{1}{2N_{\mathcal{ I}_{1}}}\sum_{i=1}^{N_{\mathcal{I}_{1}}}\lambda_{\mathcal{I}_{1},i}( \mathcal{I}_{1}[u_{i}])^{2}\] (20) \[\mathcal{L}(\mathbf{u},\mathbf{\gamma}) = \lambda_{\mathcal{N}}\mathcal{L}_{\mathcal{N}}(\mathbf{u},\mathbf{\gamma} )+\lambda_{\mathcal{B}}\mathcal{L}_{\mathcal{B}}(\mathbf{u})+\lambda_{\mathcal{I}} \mathcal{L}_{\mathcal{I}}(\mathbf{u}) \tag{21}\]

To counteract imbalances in the magnitude of the loss terms \(\mathcal{L}_{\mathcal{N}},\mathcal{L}_{\mathcal{B}},\mathcal{L}_{\mathcal{I}}\) usually each loss term is multiplied with a weight \(\lambda\) individual to that loss term. The penalty weights can be computed automatically using attention mechanisms from deep vision [40, 41], where the penalty weights are updated via a maximization w.r.t. the cost function \(\mathcal{L}\) resulting in a minimax optimization procedure. Expanding on this idea, as discussed in [42, 43], each evaluation at a collocation point is considered an individual equality constraint, thereby allocating one penalty weight per collocation point, as shown in the loss equations (18), (19), (20). This extension makes it possible to assign a
greater emphasis on more important collocation points, i.e. points which lead to larger residuals. This yields an overall minimax optimization of the form

\[\min_{\mathbf{\theta}}\max_{\mathbf{\lambda}}\mathcal{L}(\mathbf{\theta},\mathbf{\lambda}), \tag{22}\]

This minimax optimization introduces the learning rate of the maximization \(\alpha_{\max}\) as an additional hyperparameter to the learning rate for the minimization \(\alpha_{\min}\).

The optimization problem posed in equation (22) is solved using the gradients \(\nabla_{\mathbf{\theta}}\mathcal{L}\) and \(\nabla_{\mathbf{\lambda}}\mathcal{L}\) obtained via reverse-mode automatic differentiation. The Adam optimizer [44] is used in the paper at hand as an optimization procedure in all examples. However, more advanced optimizers such as combinations of Adam and L-BFGS [45] often yield a faster convergence (see, e.g. [46]).

Important differences exist between the application of PINNs to either forward or inverse problems. These differences are summarized in table 2 using three problem categories: (a) Forward problems, (b) Inverse problems with full domain knowledge and (c) Inverse problems with partial domain knowledge.

**Problem type (a) - Forward Problem**: In the forward problem, only the neural network \(A_{u}\) from equation (11) is used as an approximator of the wave field \(u(\mathbf{x},t)\), as \(\gamma(\mathbf{x},t)\) is known. The gradients of the approximated field \(\hat{u}\) are computed w.r.t. the inputs \(\mathbf{x}\) and \(t\), which enables the computation of the cost function \(\mathcal{L}(\hat{\mathbf{u}},\mathbf{\gamma})\) from equation (21). The wavefield is then solved by performing the minimax optimization from equation (22). Variations of PINNs for the forward problem of the wave equation can be found in literature [47, 48, 49, 50].

With a minor modification, the framework of PINNs can also be applied to the inverse problem, i.e. the identification of the spatially varying indicator function \(\gamma(\mathbf{x})\).

**Problem type (b) - Inverse problems with full domain knowledge of the solution field \(u(\mathbf{x},t)\)**: Since \(u(\mathbf{x},t)\) is known in the form of measurements \(\mathbf{u}_{\mathcal{M}}\) everywhere, a neural network to predict \(A_{u}\) is not needed. Instead, \(\mathbf{u}_{\mathcal{M}}\) can directly be used to evaluate the cost function \(\mathcal{L}\) from equation (21). Since the initial and boundary conditions are already approximately satisfied by the wave field, the cost function \(\mathcal{L}(\mathbf{u}_{\mathcal{M}},\hat{\mathbf{\gamma}})\) is now solely composed of the residual term of the differential equation \(\mathcal{L}_{\mathcal{N}}(\mathbf{u}_{\mathcal{M}},\hat{\mathbf{\gamma}})\). Therein, a neural network \(A_{\gamma}\) is used to predict \(\hat{\mathbf{\gamma}}\) and \(u(\mathbf{x},t)\) is directly plugged into \(\mathcal{L}_{\mathcal{N}}(\mathbf{u}_{\mathcal{M}},\hat{\mathbf{\gamma}})\) as \(\mathbf{u}_{\mathcal{M}}\). A detailed description of solving this inverse problem with PINNs can be found in [25] which documents the state-of-the-art.

**Problem type (c) - Inverse problems with partial knowledge of the solution field \(u(\mathbf{x},t)\)**: Full knowledge of the wave field, i.e. everywhere in \(\Omega\times\mathcal{T}\) is rare in the context of non-destructive testing. Typically receivers can only be placed at parts of the boundary as depicted in figure 1 such that only partial knowledge is given. Yet, to compute the differential equation loss \(\mathcal{L}_{\mathcal{N}}\) of equation (18), the wave field \(u(\mathbf{x},t)\) and the scaling function \(\gamma(\mathbf{x})\) must be known everywhere in \(\Omega\times\mathcal{T}\). Thus, in the inverse problem with partial domain knowledge, the evaluation of the differential equation loss is computed by approximating both quantities with two separate neural networks, i.e. using both equation (11) and equation (12), as visualized in figure 2. The cost \(\mathcal{L}(\hat{\mathbf{u}},\hat{\mathbf{\gamma}})\) is then composed of the boundary condition loss, equation (19), initial condition loss, equation (20), the differential equation loss \(\mathcal{L}_{\mathcal{N}}(\hat{\mathbf{u}},\hat{\mathbf{\gamma}})\), equation (18), and the measurement loss \(\mathcal{L}_{\mathcal{M}}(\hat{\mathbf{u}},\mathbf{u}_{\mathcal{M}})\), equation (9). Note, that the measurement data \(\mathbf{u}_{\mathcal{M}}\) is now incorporated through

\begin{table}
\begin{tabular}{l l l} \hline
**Problem** & **Cost Function** & **Neural Networks** \\ \hline
**(a)** Forward Problem & \(\mathcal{L}(\hat{\mathbf{u}},\mathbf{\gamma})\) & \(\hat{u}_{i}=A_{u}(\mathbf{x}_{i},t_{i};\mathbf{\theta})\) \\  & \(=\mathcal{L}_{\mathcal{N}}(\hat{\mathbf{u}},\mathbf{\gamma})+\mathcal{L}_{\mathcal{B} }(\hat{\mathbf{u}})+\mathcal{L}_{\mathcal{I}}(\hat{\mathbf{u}})\) & \\ \hline
**(b)** Inverse Problem & \(\mathcal{L}(\hat{\mathbf{u}}_{\mathcal{M}},\hat{\mathbf{\gamma}})\) & \(\hat{\gamma}_{i}=A_{\gamma}(\mathbf{x}_{i},t_{i},\mathbf{\theta})\) \\ (Full Domain Knowledge) & \(=\mathcal{L}_{\mathcal{N}}(\mathbf{u}_{\mathcal{M}},\hat{\mathbf{\gamma}})\) & \\ \hline
**(c)** Inverse Problem & \(\mathcal{L}(\hat{\mathbf{u}},\mathbf{u}_{\mathcal{M}},\hat{\mathbf{\gamma}})\) & \(\hat{u}_{i}=A_{u}(\mathbf{x}_{i},t_{i};\mathbf{\theta}_{u})\) \\ (Partial Domain Knowledge) & \(=\mathcal{L}_{\mathcal{M}}(\hat{\mathbf{u}},\mathbf{u}_{\mathcal{M}})+\mathcal{L}_{ \mathcal{N}}(\hat{\mathbf{u}},\hat{\mathbf{\gamma}})\) & and \(\hat{\gamma}_{i}=A_{\gamma}(\mathbf{x}_{i},t_{i};\mathbf{\theta}_{\gamma})\) \\  & \(+\mathcal{L}_{\mathcal{B}}(\hat{\mathbf{u}})+\mathcal{L}_{\mathcal{I}}(\hat{\mathbf{u}})\) & \\ \hline \end{tabular}
\end{table}
Table 2: Overview of differences in the PINN formulation between forward and inverse problems.
the measurement loss, evaluated with the prediction of the forward solution \(\hat{\mathbf{u}}\). The current state-of-the-art of this inverse problem with PINNs is given e.g. in [26] with a representative Python code provided in [51].

## 2 Effect of Neural Networks in the Inversion

To investigate the strengths and weaknesses of PINNs for full waveform inversion, the PINN and its training procedure as illustrated in figure 2 is incrementally modified until it is indistinguishable from the classical adjoint optimization.

The work in [25, 26] discusses the use of PINNs for full waveform inversion, where fully-connected neural networks with inputs \(\mathbf{x}\) and \(t\) were employed. By constrast, in the paper at hand convolutional neural networks are employed, since their spatial equivariance leads to a significant speed-up in training and faster convergence. This has already been reported in a more general context in computational mechanics in e.g. [10, 39] and is not the main focus of the paper at hand. This modification additionally enables a better transition and comparison to the presented methods and is not seen as a detriment to the PINN methodology. Despite this minor difference, we will still consider [26] as the reference case for what is possible within the current state-of-the-art using pure PINNs and only partial domain knowledge.

### Physics Informed Neural Networks with non-trainable Forward Operator

The first modification of PINNs visualized in figure 2 is to replace the neural network \(A_{u}(\mathbf{\theta})\), whose task is to learn the forward solution, with a non-trainable forward operator \(F\) with \(\hat{\gamma}\) as additional input. The new structure of a Physics-Informed Neural Network with a non-trainable Forward Operator is illustrated in figure 3. The forward operator could be any type of surrogate model such as pre-trained neural networks, or classical solution methods such as the finite element method [18] or the finite difference method [19]. An important computational requirement on the forward operator \(F\) is the viability of a GPU-based implementation, which is an advantage of pre-trained neural networks over CPU-based solvers. However, finite difference methods are commonly also accelerated with GPUs, as e.g. in [52]. Therefore, finite differences are used in the sequel as the forward operator \(F\). A basic but relatively efficient implementation in PyTorch is used by exploiting its convolutional network infrastructure as stencils. This results in simulation runs which take about \(500\upmu\)s wall clock time for all 1200 timesteps on a \(252\times 124\) grid1. Details are laid out in Appendix A.

Footnote 1: tested on an Nvidia A100-SXM4-40GB

The introduction of a non-trainable forward operator should, in the ideal case, lead to very small values of the loss components related to the initial-boundary value problem, i.e. \(\mathcal{L}_{\mathcal{N}},\mathcal{L}_{\mathcal{B}},\mathcal{L}_{\mathcal{I}}\). Typically, these are only approximately zero due to modelling and discretization errors. Their

Figure 2: Prediction and training of PINNs for inversion with the scalar wave equation. Predictions of the wavefield \(\hat{u}\) and indicator \(\hat{\gamma}\) are obtained through the neural networks \(A_{u}(\mathbf{\theta}_{u})\) and \(A_{\gamma}(\mathbf{\theta}_{\gamma})\). These are updated during training with the gradients of the cost function \(\mathcal{L}\) obtained by automatic differentiation.
contribution is, however, still negligible in the presence of the measurement loss \(\mathcal{L}_{\mathcal{M}}\). For this reason, they can be eliminated from the cost function \(\mathcal{L}\), such that it is now solely composed of the measurement loss \(\mathcal{L}_{\mathcal{M}}\).

The introduction of a non-trainable forward operator leads to a significantly less complex optimization task. In the case of a time-stepping scheme as forward operator \(F\), the disadvantage of the design depicted in figure 3 is a more expensive sensitivity computation \(\nabla_{\boldsymbol{\theta}}\mathcal{L}(F(\boldsymbol{\hat{\gamma}}))\ =\ \nabla_{ \boldsymbol{\hat{\gamma}}}\mathcal{L}(F(\boldsymbol{\hat{\gamma}}))\cdot \nabla_{\boldsymbol{\theta}}\boldsymbol{\hat{\gamma}}(\boldsymbol{\theta})\) because the backpropagation needs to be piped through each timestep of the forward operator \(F\) by means of automatic differentiation to compute \(\nabla_{\boldsymbol{\hat{\gamma}}}\mathcal{L}(F(\boldsymbol{\hat{\gamma}}))\). The cost of computing the gradient of the inverse field w.r.t. the network parameters \(\nabla_{\boldsymbol{\theta}}\boldsymbol{\hat{\gamma}}(\boldsymbol{\theta})\) is unchanged.

### Adjoint Method

The next step is to replace the neural network Ansatz of the indicator function \(\gamma\) with a piecewise polynomial Ansatz. Here, a piece-wise constant Ansatz is selected and defined in the one-dimensional case as \(\hat{\gamma}(x)=\sum_{i=1}^{N_{c}}c_{i}N_{i}(x)\) with

\[N_{i}(x)=\begin{cases}1\text{ if }x_{i}\leq x<x_{i+1}\\ 0\text{ else.}\end{cases} \tag{23}\]

For the two-dimensional and three-dimensional cases the discretization is handled in the same way for the additional dimensions. This modification is depicted in figure 4. The coefficients \(\boldsymbol{c}=\{c_{i}\}_{i=1}^{N_{c}}\) are still updated with the Adam optimizer using the gradients obtained through automatic differentiation. The reason for choosing a piecewise constant Ansatz over a piecewise higher-order polynomial is that a piecewise constant Ansatz is already capable of representing the jump of going from a constant material to a void, as is the case in all examples presented in section 3. Without the neural network, the Sigmoid activation function is missing and the range of \(\hat{\gamma}\) has to be constrained in a different manner. To this end, an additional weight, i.e. coefficient clipping is used, where the absolute value is clipped to the range of a small value \(\epsilon\) close to zero and an upper bound one or slightly exceeding one [37].

Given equation (23), the sensitivity w.r.t. the coefficients

\[\nabla_{\boldsymbol{c}}\mathcal{L}(F(\boldsymbol{\hat{\gamma}}))\ =\ \nabla_{ \boldsymbol{\hat{\gamma}}}\mathcal{L}(F(\boldsymbol{\hat{\gamma}}))\cdot \nabla_{\boldsymbol{c}}\boldsymbol{\hat{\gamma}}(\boldsymbol{c}) \tag{24}\]

can either be computed via automatic differentiation or the continuous adjoint method, as indicated in figure 4. In the continuous adjoint method, the gradient is derived analytically before its discretization, as e.g. derived in [37]. This third modification has transformed the approach into the classical full waveform inversion based on finite differences. No major implications on the accuracy are expected due to this modification, as the gradients are the same, apart from numerical errors introduced by the finite difference discretization in the adjoint method. These originate in the computation of the Frechet kernel

\[K_{\gamma}=\int_{\mathcal{T}}-\rho_{0}u_{t}^{\dagger}u_{t}+\rho_{0}c_{0}^{2} \nabla u^{\dagger}\cdot\nabla udt \tag{25}\]

Figure 3: PINN with non-trainable forward operator. The modification of the original PINN from figure 2 is marked in red. Instead of learning the forward problem, a non-trainable forward operator \(F\) with \(\hat{\gamma}\) as additional input is used. Simultaneously, this leads to a simplification of the cost function \(\mathcal{L}\), as \(\mathcal{L}_{\mathcal{N}},\mathcal{L}_{\mathcal{B}},\mathcal{L}_{\mathcal{I}}\) are approximately satisfied due to \(F\).
used to compute the gradient \(\nabla_{\hat{\gamma}}\mathcal{L}_{\mathcal{M}}\) via the integration of \(K_{\gamma}\) over the spatial and time domains (c.f. [37] for more details). To compute equation (25), the temporal and spatial gradients of the solution field \(u\) and the adjoint field \(u^{\dagger}\) are computed with finite differences. Note, that the sensitivity w.r.t. the coefficients \(c_{i}\) is computed in a second step, by multiplication with the corresponding shape function \(N_{i}(x)\), obtained by deriving equation (23) w.r.t. \(c_{i}\), i.e. \(\nabla_{\mathbf{c}}\hat{\gamma}(\mathbf{c})\).

In fact, the key difference between the adjoint method and automatic differentiation is the order in which the gradient computation and discretization of the state equation is performed. In the continuous adjoint method, the gradient is first derived analytically, and thereupon the state equation is discretized as described above. By contrast, in the case of automatic differentiation, the differential equation is discretized first, and the gradient is propagated through each timestep. Thus, automatic differentiation bears a great similarity to the discrete adjoint method and is, as stated in [53], fundamentally the same approach. Manually deriving the sensitivities in the discrete adjoint or computing them through reverse mode automatic differentiation as in [53, 54] will not result in a difference in accuracy. Thus, regarding the accuracy, the difference between automatic differentiation and the continuous adjoint method is only affected by the difference between the discrete and continuous adjoint method.

Unlike the accuracy, the computational cost of the automatic differentiation approach compared to the continuous adjoint method is different. Automatic differentiation adds a significant memory cost due to the backpropagation. In addition, the backpropagation along all timesteps is costly with automatic differentiation, as also mentioned in the work on the ordinary differential equation networks [55]. The theoretical complexity of the backpropagation and continuous adjoint method w.r.t. the number of timesteps is in both cases \(O(n)\)[56] and in practice only a factor of roughly two is observed between the computation times. To save memory, a truncated backpropagation through time [56, 57] can be used, as typically carried out for recurrent neural networks. However, a truncation was not used in the examples investigated in section 3 since it was found to be detrimental to the learning behavior of the indicator function.

### Hybrid Approach

This section introduces a hybrid approach which combines the advantages of a neural network Ansatz for the scaling function \(\gamma\) and the more efficient gradient computation of the gradient of the loss functional equation (9) using the continuous adjoint method equation (25). To this end, the neural network is reintroduced as an Ansatz for \(\gamma\). In analogy to equation (24), the gradient of the cost function can be computed with the chain rule but must now be computed with respect to the neural network parameters \(\mathbf{\theta}\):

\[\nabla_{\mathbf{\theta}}\mathcal{L}(\mathbf{\hat{u}}(\mathbf{\hat{\gamma}}(\mathbf{\theta})) )=\nabla_{\mathbf{\hat{\gamma}}}\mathcal{L}(\mathbf{\hat{u}}(\mathbf{\hat{\gamma}}))\cdot \nabla_{\mathbf{\theta}}\mathbf{\hat{\gamma}}(\mathbf{\theta}) \tag{26}\]

The sensitivity of the cost function w.r.t. the neural network outputs, i.e. the indicator field, \(\nabla_{\mathbf{\hat{\gamma}}}\mathcal{L}\) is efficiently computed with the continuous adjoint method equation (25), while the sensitivity of the indicator field w.r.t. the neural network parameters \(\nabla_{\mathbf{\theta}}\mathbf{\hat{\gamma}}\) is obtained via backpropagation through the neural network \(A_{\gamma}(\mathbf{\theta})\). The hybrid approach results in similar computational times as

Figure 4: Adjoint method. The modification of the previously modified PINN from figure 3 is indicated in red. The neural network Ansatz is replaced by a constant Ansatz parametrized by the coefficients \(\mathbf{c}\).
the adjoint method for a sufficient number of timesteps. The scheme is illustrated in figure 5. A similar scheme was recently proposed in the context of topology optimization [58].

In addition to saving memory compared to the scheme discussed in the previous section, the hybrid scheme bears the advantage of separating neural networks from forward solvers such that an efficient GPU implementation of the forward solver is no longer strictly necessary. Instead, the neural network can provide the indicator prediction \(\mathbf{\hat{\gamma}}\) and its sensitivity w.r.t. the parameters \(\nabla_{\mathbf{\theta}}\mathbf{\hat{\gamma}}\) on a CPU while, the forward solver can also be executed independently on a CPU. The increase in computational effort of evaluating and updating the neural network on a CPU instead of a GPU is very small compared to the cost of the forward solver. Furthermore, the forward solver can even be a completely independent black-box solver which provides the wave-field \(\hat{u}\) and the sensitivity of the loss functional \(\nabla_{\hat{\gamma}}\mathcal{L}\).

## 3 Configuration

Two model problems with increasing complexity are considered. First, a two-dimensional plate is used to determine the strengths and weaknesses of neural networks in the inversion process. Afterwards, a more complex three-dimensional case is used to compare the proposed hybrid approach to the adjoint method. The forward operator is always a finite difference solver. To avoid an inverse crime [59], the reference data was generated with a boundary-conforming finite element solver for the partial domain knowledge cases (Voxel-FEM for the three-dimensional case). However, in the case of the two-dimensional full domain knowledge case the finite cell method [60] is employed for the reference data generation. In all studies, the sine burst, presented in equation (8) is applied with \(n_{c}=2\) cycles, a frequency of \(f_{\psi}=500\) kHz, and an amplitude of \(\psi_{0}=10^{12}\) kg(ms)\({}^{-2}\). Special care must be taken when using finite differences during the source application in equation (7). The Dirac delta has to be scaled by a factor of \(1/||d\mathbf{x}||_{2}^{2}\) to account for the grid spacings [61]. Additionally, both cases share the same material with a density of \(\rho_{0}=2700\) kg/m\({}^{3}\) and a wave speed of \(6000\) ms\({}^{-1}\).

### Two-Dimensional Case

The two-dimensional case is strongly inspired by [37], where a circular hole with a diameter of 5 mm is to be detected on a 100 mm \(\times\) 50 mm rectangular plate, illustrated in figure 6(a). We do not provide specific information about the hole to the machine learning algorithms, which implies the entire material distribution \(\gamma(\mathbf{x})\) has to be estimated. While the geometry, material, and source type are identical, the sensor and source locations are slightly different from [37]. In total, 54 sensors and four sources are used on the top edge in four separate measurements. The exact configuration is provided in figure 6, showing the upper left edge of the domain using the domain definitions from figure 1. The setup is symmetric along the vertical midline indicated by the dashed grey

Figure 5: Hybrid approach. Modification of the previously modified PINN from figure 3 is marked in red. The gradient computation is now split into two parts. First, the sensitivity of the cost function \(\mathcal{L}\) w.r.t. the predicted indicator field \(\hat{\gamma}\) is computed with the adjoint method using equation (25), whereupon the sensitivity of the indicator field \(\hat{\gamma}\) w.r.t. the neural network parameters \(\mathbf{\theta}\) is found with the automatic differentiation used to train the neural network. This yields the desired gradient \(\nabla_{\mathbf{\theta}}\mathcal{L}\) by using the chain rule equation (26).
line. The slight peculiarity of the sensor and source locations is due to them being aligned to the finite difference nodes and neural network grid points. This is, however, not a limitation of the presented methods, as it can easily be circumvented by, e.g. using Kaiser windowed sinc functions as interpolation of the sources [61] and linear interpolation for the sensors.

A forward simulation is illustrated in figure 7, where the third source emits a pulse at \(t=0\) s, which then leads to a disturbance of the wavefield in the center due to the circular void. The finite difference simulation for the inversion uses \(252\times 124\) grid points and two ghost cells per dimension, a lower indicator function bound of \(\epsilon=10^{-5}\), and a timestep size of \(\Delta t=3\cdot 10^{-8}\) s\({}^{-1}\) for 1200 timesteps.

### Three-Dimensional Case

The purpose of the three-dimensional case is to showcase the possible strengths of neural networks when the adjoint optimization is used in a setting where it does not perform well anymore. For this, two cuboid slices from a CT-scanned drill core [62] are considered and illustrated in figure 8. The edge length is 20 mm and discretized by 92 grid points in each dimension during the inversion, where 1200 timesteps with a timestep size of \(\Delta t=1.8\cdot 10^{-8}\) s\({}^{-1}\) are considered. Again a lower indicator function bound of \(\epsilon=10^{-5}\) is used, while the upper bound is increased to 1.2 for the adjoint method, as this was found to be beneficial. The sources and sensors are applied to the backside of the cubes, as indicated by the darker shading on the right of figure 8. The four sources are placed at a distance of 5.27 mm from each edge, while the sensors coincide with the grid points of the plane.

The case was selected such that a single forward and adjoint simulation can be performed within the limits of the 40 GB GPU memory. To compute larger three-dimensional, efficient parallelizations of the simulation have to be used, which can either be implemented on GPU clusters or alternatively on CPUs using efficient communication between the GPU and CPU for the neural network.

Figure 6: Symmetrical sensor (black circles) and source (red circles) configuration.

Figure 7: The wavefield is computed with the underlying material distribution, depicted in the leftmost figure, with the third source. When the waves reach the defect, a disturbance in the wavefield can be seen.
### Neural Network Architecture

The neural network considered as Ansatz for the indicator function \(\hat{\gamma}\) has a generator structure with pixel-wise normalization inspired by generative adversarial networks [63]. It upsamples a random three-dimensional matrix of size \(128\times 8\times 4\) with subsequent convolutional filters. The upsampling and convolutions are repeated five times to get the desired resolution for the indicator function \(\hat{\gamma}\), see figure 9 for the neural network used for the two-dimensional case. A similar network is employed in the three-dimensional case using three-dimensional convolutions and a four-dimensional random input matrix of size \(256\times 3\times 3\times 3\). To increase the ratio between the number of indicator voxels and network parameters, twice as many filters are employed in the three-dimensional case. For more details on both networks, c.f. Appendix B. Nearest-Neighbor interpolation was used for upsampling, as it performed better than polynomial interpolation, as it preserves discontinuities needed to represent the voids. Due to their convergence accelerating behavior [64], adaptive activation functions are used. Specifically, PReLU [65] was used for all layers except the last layer, where an adaptive Sigmoid was used to ensure an output in the range \([0,1]\). The output was then rescaled to the range \([\epsilon,1]\) to ensure numerical stability. Optimization is always conducted with Adam [44] to make the comparison more straightforward, although L-BFGS leads to fewer iterations. Gradient clipping [66, 67, 68] was found to be essential to ensure an early convergence. Without it, a plateau was observed, followed by a strongly delayed breakthrough. In addition, a learning rate scheduler with a polynomial decay \((b\cdot\text{epoch}+1)^{a}\) was employed with the values \(a=-0.5\), \(b=0.2\). The initial network weights are set with Glorot initialization [69].

Figure 8: The two cases extracted from a CT-scanned drill core [62]. The sources and sensors are applied on the backside of the cubes indicated by the darker surface on the right.

Figure 9: The convolutional neural network for the two-dimensional inversion case. It performs upsampling and convolutions in a repeating pattern on a random input of size \(128\times 4\times 8\), yielding a prediction \(\hat{\gamma}\) of size \(1\times 126\times 254\). For more details, see Appendix B.
## 4 Numerical Results

### Two-Dimensional Case

In the subsequent sections, the individual components of PINNs will be investigated and compared to the classical adjoint optimization using the case presented in section 3.1. All simulations are performed on an Nvidia A100 SXM4 40GB GPU and the same convolutional neural network presented in figure 9 and Appendix B.

#### 4.1.1 Physics-Informed Neural Networks

The results presented in this section refer to case PINNs, table 1, whose methodology is described in section 1.3. The various inversion tasks presented in [26] ran for about \(400\,000-600\,000\) epochs with computational times of about \(0.2\) s per epoch2. The need for such a large number of epochs illustrates the complexity of the nested optimization task of simultaneously learning the forward and inverse solution for a single problem. Additional complications arise when multiple sources are applied, and a forward operator taking the different sources into account must be learned. In [26], this problem was circumvented by using source stacking [70], such that only one forward simulation had to be learned.

If only one field has to be learned, as in the case of inversion with full domain knowledge, as also presented in [25], the complexity of the optimization task decreases. To illustrate the decrease in complexity and the capabilities of PINNs in specialized settings, the problem presented in section 3.1 is solved with a PINN using additional sensors. The sensors are placed throughout the entire domain, coinciding with the finite difference grid points. As the measurements are only given at discrete points, the residual of the differential equation (18) is evaluated numerically using finite differences for the derivative terms. Additionally, the indicator function is prescribed at the boundaries, i.e. set to one to ensure material at the boundaries, which was essential to retrieve an accurate inversion. The resulting inversion after 100 epochs is provided in figure 10 with a clear identification of the circular void and only minor artifacts. It was essential to prescribe the density \(\rho(\mathbf{x})\) at the boundary and use a relatively large learning rate \(\alpha_{\min}=2\cdot 10^{-3}\). Furthermore the use of the learnable penalty weights using the attention mechanism was helpful with a learning rate of \(\alpha_{\max}=2\cdot 10^{-2}\).

Footnote 2: Based on the Python script provided by [51] and tested on an Nvidia A100-SXM4-40GB GPU.

It is important to emphasize that the improvement from \(400\,000-600\,000\) epochs to 100 epochs is a result of the simplification of the physical problem by using full domain knowledge, which should not be taken as a reference of what is possible with PINNs in the partial domain knowledge case

Figure 10: Prediction and error of the material distribution after 100 epochs using a PINN with full domain knowledge. Average elapsed time per epoch is \(0.51\) s.
presented in section 3.1. Instead, \(400\,000-600\,000\) epochs are to be considered as the reference of the capabilities of PINNs for a similar inversion task.

#### 4.1.2 Physics-Informed Neural Networks with non-trainable Forward Operator

An alternative to reducing the complexity of learning both the forward and inverse solutions is using a non-trainable forward operator classified in table 1 as PINNs with non-trainable forward operator and described in section 2.1. The neural network responsible for the forward solution is now replaced by a finite difference scheme. A reasonable inversion, shown in figure 11, is reached after only 50 epochs. Thus the convergence w.r.t. iterations is even faster than the PINN with full domain knowledge, with the detriment of a greater elapsed time per epoch at 2.53 s. Overall, considering the \(400\,000-600\,000\) epochs observed in [26] given only partial domain knowledge, the PINN is outperformed by using a non-trainable forward operator. For the optimization a learning rate of \(\alpha_{\min}=2\cdot 10^{-3}\) was employed with gradient clipping. The learnable penalty weights improved the results slightly, but were not found to be essential.

#### 4.1.3 Adjoint Method

Next, a performance evaluation of the adjoint method is provided in table 1. Section 2.2 then contains more details and a special focus on the importance of the Ansatz and the gradient computation. A piece-wise constant Ansatz is considered instead of the neural network using automatic differentiation for the gradient computation. The void is clearly identified in the prediction after 50 epochs, shown in figure 12, and achieved with a slightly increased learning rate of \(\alpha_{\min}=6\cdot 10^{-2}\). However, additional oscillatory artefacts are observed. This coincides with the observations in [37], where linear shape functions are used as the Ansatz of the unknown material field. However, interestingly changing the Ansatz of \(\gamma\) to a neural network leads to a very different and better local minimum. Specifically, the mean squared error is reduced from \(3.2\cdot 10^{-3}\) using the piecewise constant Ansatz for \(\gamma\) as in equation (23) to \(9.5\cdot 10^{-4}\) using a neural network as in equation (12). Furthermore, the smoothness of the solution is increased in regions where the true solution is constant, while the sharp gradients around the circular void are retained as depicted in figure 11 and quantified in Appendix C. It is noteworthy that this increase in accuracy is achieved without a significant increase in computational effort.

However, one major difference accompanying the change in Ansatz must be mentioned. The initial field of the constant Ansatz is chosen to be equal to one throughout the entire domain, which in the neural network case is impossible without pre-training. Pre-training to learn a constant field of one with a large regularization on the neural network was not found to be helpful. Thus the neural network starts with a random field, as later shown in the prediction histories in figure 16. However,

Figure 11: Prediction and error of the material distribution after 50 epochs using automatic differentiation and a neural network Ansatz. Average elapsed time per epoch is 2.53 s.
introducing a random field as initialization for the constant Ansatz is found to be detrimental, observable in an unsuccessful convergence toward a meaningful inversion. Thus, the strength of the neural network Ansatz must lie in its non-linear connectivity and strong over-parametrization, \(526\,252\) parameters compared to the \(252\times 124=31\,248\) parameters for the constant Ansatz.

It has to be noted that increasing the number of parameters in the constant Ansatz might be beneficial and yield a better comparison to the neural network Ansatz. However, when using finite differences, increasing the number of parameters, by e.g. increasing the number of voxels or the polynomial degree per finite difference grid point do not yield an improved indicator discretization. These methods are possible in the context of the finite element and cell method, where the added parameters can be taken into account via multiple Gauss points per element. Thus the benefit of neural networks as discretization becomes even clearer, as they provide a universal tool independent of the simulation method.

With this final modification, the change to the classical adjoint optimization as presented in [37] is almost complete, with the minor difference of using automatic differentiation for the gradient computation. However, as previously stated, moving from automatic differentiation to the continuous adjoint method for the gradient computation is not likely to have a significant impact. To assess this numerically, the initial gradient of the cost function w.r.t. the inverse field \(\nabla_{\boldsymbol{c}}\mathcal{L}\) using a constant Ansatz is in both cases compared to a reference solution obtained with finite differences, by varying one parameter \(c_{i}\) at a time. Only minor differences are observable at the boundaries, illustrated in figure 13. The difference can be explained by the numerical evaluation of the Frechet kernel, equation (25), which introduces the most significant numerical errors at the boundaries due to the central finite difference scheme. The absolute error w.r.t. the reference solution is, however, nonetheless negligible. Thus only minor differences in the learned prediction from figure 12 are expected, as confirmed by the numerical comparison. The tiny differences accumulate throughout the iterations leading to slightly different local minima with a similar quality of inversion. Thus, no improvement is expected in the solution quality when using either automatic differentiation or the adjoint method for the gradient computation. However, using the adjoint method to compute the gradient \(\nabla_{\hat{\gamma}}\mathcal{L}_{\mathcal{M}}\) leads to a significant drop in the computational effort from \(2.46\) s to \(1.05\) s per epoch for a simulation comparable to the one presented in figure 12.

#### 4.1.4 Hybrid Approach

The proposed hybrid approach classified in table 1 and described in detail in section 2.3 combines the advantage of the neural network Ansatz with the lower computational effort of the adjoint method for the gradient computation \(\nabla_{\hat{\gamma}}\mathcal{L}_{\mathcal{M}}\). The resulting prediction is presented in figure 14 with only a minor decrease in prediction quality when compared to figure 11, while retaining computation times of \(1.06\) s per epoch. A slight increase of the learning rate \(\alpha_{\min}=4\cdot 10^{-3}\) was

Figure 12: Prediction and error of the material distribution after 50 epochs using automatic differentiation and a constant Ansatz. Average elapsed time per epoch is \(2.46\) s.
Figure 13: Initial gradient using a constant Ansatz computed with automatic differentiation (a) and the adjoint method (b). The corresponding absolute errors in (c) and (d) are computed with a numerical reference solution obtained with finite differences.
beneficial.

#### 4.1.5 Comparison

To conclude the comparison between the investigated methodologies, the training histories are juxtaposed in figure 15. While the constant Ansatz enables the best minimization of the cost function, the mean squared error is lowest for the neural network Ansatz, confirming once again the advantage of the neural network Ansatz. Furthermore, the impact of using automatic differentiation or the adjoint method with finite differences becomes clearer with a consistently lower cost and mean squared error when using automatic differentiation. In view of the computational benefit, the slight loss in accuracy is acceptable. This difference only occurs due to the finite difference evaluations of the Frechet Kernel, equation (25). In the case of finite elements, the Frechet Kernel can be evaluated exactly by analytical differentiation of the shape functions.

To gain an additional understanding of the training behaviour, consider the history of indicator predictions provided in figure 16. Here, automatic differentiation with a neural network Ansatz, the adjoint method with a constant Ansatz and the hybrid approach are juxtaposed as snapshots showing the entire domain at specific instances of the iteration. In all three cases, the void is essentially recovered after 15 epochs. After the identification, only minor improvements in the remaining domain are made, where the effect is strongest for the hybrid approach and weakest for the adjoint method. Interestingly, all three approaches lead to three completely different local minima within a comparable amount of computational effort, despite initial noise in the case of the neural networks.

Thus two critical conclusions arise from this investigation. Firstly, learning the forward and inverse fields simultaneously for a single inversion problem yields a much more complex optimization task than using an non-trainable forward operator such as finite differences. The improvement of using a non-trainable forward operator is about four orders of magnitude when considering the number of iterations. Secondly, a neural network as Ansatz for the inverse fields outperforms classical approaches, such as a piece-wise constant Ansatz, with regard to the recovered field. To once again highlight the difference, a surface plot of the prediction of both the neural network Ansatz and the constant Ansatz after optimization using automatic differentiation is provided in figure 17. The surfaceplot illustrates the difference in smoothness of the two solutions. Although the neural network solution is smoother, the sharpness retained (quantified in Appendix C). The plot indicates a regularization effect to be studied further using the three-dimensional case in the upcoming section 4.2.

Figure 14: Prediction and error of the material distribution after 50 epochs using the combination of automatic differentiation and adjoint method and a neural network Ansatz. Average elapsed time per epoch is 1.06 s.
Figure 16: Prediction histories starting at epoch 0 until epoch 25. From left to right, automatic differentiation with a neural network Ansatz, the adjoint method with a piece-wise constant Ansatz, and the hybrid approach with a neural network Ansatz are depicted.

Figure 15: Comparison of the training histories for each presented approach. The cost function, which consists of the measurement loss \(\mathcal{L}_{\mathcal{M}}\) is shown on the left, and the mean squared error (MSE) w.r.t. the true material distribution \(\gamma(\mathbf{x})\) is on the right. The normalization of cost and mean squared error is with respect to the corresponding values with an indicator function equal to one in the entire domain.
### Three-Dimensional Case

Due to the quality of the inversion results using neural networks as Ansatz, a three-dimensional case is considered. As the additional spatial dimension introduces an additional computational burden, only the most efficient and simultaneously most promising methods are compared, i.e. the hybrid method and the adjoint method with a piece-wise constant Ansatz. The inversion results of the two approaches are illustrated in figure 18 in terms of identified voids for the two CT-scan slices from figure 8. Once again, the neural networks enable a superior reconstruction, quantifiable by considering the mean squared error provided as annotations in figure 18. Center slices are shown in figure 19 to provide an insight into the overall distribution and smoothness of the solution. From this, it can be observed that the reconstruction of the adjoint method suffers from the oscillatory artefacts to such an extent that the voids are no longer clearly identified. The hybrid approach exhibits much less artefacts and is able to recover the voids with superior accuracy.

We attribute this superiority to the different characteristics of the neural network Ansatz which not only leads to an over-parametrization of the indicator field but also has an in-built non-linearity. Specifically, \(6\,306\,764\) parameters were used for the neural network compared to the \(92^{3}=778\,688\) parameters in the classical discretization. Interestingly, the reconstruction is worse when using fewer neural network parameters. By using half and quarter as many filters, yielding only \(1\,577\,452\) and \(394\,748\) parameters, the quality of the reconstruction decreases. At the same time, using double the amount of filters, resulting in \(25\,221\,004\) parameters, leads to a minor decrease in quality, due to the now dominating regularization effect of the over-parameterization. The mean squared error (MSE) w.r.t. the true material distribution presented for the three variations of the network in figure 20 quantifies this observation. Furthermore, a smoothing effect is noticeable, with an increase of parameters indicating a regularization-like behaviour. Thus, the over-parametrization enabled by the neural network is one, yet not the only key ingredient in enhancing the inversion results.

Several other aspects distinguish the neural network from the piece-wise constant Ansatz. The non-linear activations inside the neural network represent an important difference. Also, the neural network acts as a global function, while the constant Ansatz is defined on local patches. Furthermore, other network architectures, such as fully-connected or graph neural networks, may show different behaviour. Currently, it is unclear which aspects are critical aside from the over-parametrization. All aspects can be studied in a straightforward empirical manner and are subject

Figure 17: Surface plot of the predicted material distribution after 50 epochs. The training was performed with a neural network and constant Ansatz using automatic differentiation.
to further research.

## 5 Conclusion

The aim of this work was to compare PINNs with the classical adjoint optimization for full waveform inversion and thereby quantify the strengths of each approach by incremental modification. Firstly, estimating the forward solution via a nested optimization using two neural networks was shown to be disadvantageous. Instead, a classical forward solver, as in the adjoint optimization, is recommended. Secondly, a neural network as Ansatz of the inverse field is highly beneficial with a regularization-like effect without a loss of sharpness. Thirdly, no significant difference in accuracy exists in the sensitivity computation when comparing automatic differentiation and the continuous adjoint method. However, the continuous adjoint method is preferred in the gradient estimation of the cost w.r.t. the inverse field due to the lower computational effort.

As a consequence of these insights, the hybrid method is proposed, utilizing a neural network to discretize the inverse field and computing the cost function gradient w.r.t. the inverse field with the adjoint method, which then is passed on to the backpropagation algorithm to update the neural network parameters. Comparable computational times are achieved, while the inverse reconstructions improved.

The performance of neural networks as the discretization of the inverse fields enables not only

Figure 18: Identified void by the adjoint method and the hybrid approach after 100 epochs in comparison to the true void. The illustrated void is extracted by thresholding indicator values lower than 0.5. Average elapsed time per epoch is about 29.1 s for both methods.
Figure 19: Center slices of the indicator function predictions using the adjoint method and the hybrid approach after 100 epochs. The sliced planes are annotated in the top. Average elapsed time per epoch is 29.1 s for the hybrid approach and 29.1 s for the adjoint method.

Figure 20: Mean squared error after 100 epochs for different number of parameters illustrated on case 1.
further research for other inverse problems but also a thorough investigation into the underlying cause of its success. Insight can be gained by considering the difference between global and local functions, different over-parametrizations, i.e. the number of parameters, the non-linearities in the network, different network architectures, and possibly a thorough investigations of the initialization of the wavefield and its effect on the inversion.

## Acknowledgements

The authors gratefully acknowledge the funding through the joint research project Geothermal-Alliance Bavaria (GAB) by the Bavarian State Ministry of Science and the Arts (StMWK) wich finance Leon Herrmann as well as the Georg Nemetschek Institut (GNI) under the project Deep-Monitor. Furthermore, we gratefully acknowledge funds received by the Deutsche Forschungsgemeinschaft under Grant no. KO 4570/1-1 and RA 624/29-1 which support Tim Burchner. Felix Dietrich would like to acknowledge the funds received by the Deutsche Forschungsgemeinschaft - project no. 468830823.

## Declarations

### Conflict of interest

No potential conflict of interest was reported by the authors.

### Data Availability

We provide a PyTorch [7] implementation of all methods in [71].
## Appendix A Efficient Finite Difference Implementation in Deep Learning Frameworks

The explicit finite difference scheme, shown for the two-dimensional case,

\[u_{i,j}^{n+1}=-u_{i,j}^{n-1}+2u_{i,j}^{n}\] \[+\frac{2}{\gamma_{i,j}}\left(\frac{c_{0}\Delta t}{\Delta x}\right) ^{2}\left(\left(\frac{1}{\gamma_{i,j}}+\frac{1}{\gamma_{i+1,j}}\right)^{-1}(u_{ i+1,j}^{n}-u_{i,j}^{n})-\left(\frac{1}{\gamma_{i-1,j}}+\frac{1}{\gamma_{i,j}} \right)^{-1}(u_{i,j}^{n}-u_{i-1,j}^{n})\right)\] \[+\frac{2}{\gamma_{i,j}}\left(\frac{c_{0}\Delta t}{\Delta y}\right) ^{2}\left(\left(\frac{1}{\gamma_{i,j}}+\frac{1}{\gamma_{i,j+1}}\right)^{-1}(u_{ i,j+1}^{n}-u_{i,j}^{n})-\left(\frac{1}{\gamma_{i,j-1}}+\frac{1}{\gamma_{i,j}} \right)^{-1}(u_{i,j}^{n}-u_{i,j-1}^{n})\right)\] \[+\frac{\Delta t^{2}}{\rho_{0}\gamma_{i,j}}f_{i,j}^{n}\]

was implemented in PyTorch [7] using convolutional kernels \(\mathbf{K}\), thus exploiting PyTorch's GPU capabilities.

\[\mathbf{u}^{n}=-\mathbf{u}^{n}+2\mathbf{u}^{n-1}\] \[+\frac{1}{\mathbf{\gamma}}\Bigg{(}\frac{1}{\frac{1}{\mathbf{\gamma}}*\mathbf{ K}_{\gamma_{2_{0}}}}(\mathbf{u}*\mathbf{K}_{u_{u_{0}}})-\frac{1}{\frac{1}{\mathbf{\gamma}}*\mathbf{ K}_{\gamma_{2_{1}}}}(\mathbf{u}*\mathbf{K}_{u_{u_{1}}})\] \[+\frac{1}{\frac{1}{\mathbf{\gamma}}*\mathbf{K}_{\gamma_{3_{0}}}}(\mathbf{u}* \mathbf{K}_{u_{u_{0}}})-\frac{1}{\frac{1}{\mathbf{\gamma}}*\mathbf{K}_{\gamma_{3_{1}}}}( \mathbf{u}*\mathbf{K}_{u_{y_{1}}})\] \[+\frac{\Delta t^{2}}{2\rho_{0}}\mathbf{f}\Bigg{)}\]

The corresponding kernels are defined as

\[\mathbf{K}_{\gamma_{a_{0}}}=\begin{pmatrix}0&0&0\\ 0&1&1\\ 0&0&0\end{pmatrix},\mathbf{K}_{u_{a_{0}}}=\begin{pmatrix}0&0&0\\ 0&-\left(\frac{c_{0}\Delta t}{\Delta x}\right)^{2}&\left(\frac{c_{0}\Delta t} {\Delta x}\right)^{2}\\ 0&0&0\end{pmatrix},\] \[\mathbf{K}_{\gamma_{x_{1}}}=\begin{pmatrix}0&0&0\\ 1&1&0\\ 0&0&0\end{pmatrix},\mathbf{K}_{u_{x_{1}}}=\begin{pmatrix}0&0&0\\ -\left(\frac{c_{0}\Delta t}{\Delta x}\right)^{2}&\left(\frac{c_{0}\Delta t}{ \Delta x}\right)^{2}&0\\ 0&0&0\end{pmatrix},\] \[\mathbf{K}_{\gamma_{y_{0}}}=\begin{pmatrix}0&0&0\\ 0&1&0\\ 0&1&0\end{pmatrix},\mathbf{K}_{u_{y_{0}}}=\begin{pmatrix}0&0&0\\ 0&-\left(\frac{c_{0}\Delta t}{\Delta y}\right)^{2}&0\\ 0&-\left(\frac{c_{0}\Delta t}{\Delta y}\right)^{2}&0\\ 0&\left(\frac{c_{0}\Delta t}{\Delta y}\right)^{2}&0\end{pmatrix},\] \[\mathbf{K}_{\gamma_{y_{1}}}=\begin{pmatrix}0&1&0\\ 0&1&0\\ 0&0&0\end{pmatrix},\mathbf{K}_{u_{y_{1}}}=\begin{pmatrix}0&-\left(\frac{c_{0} \Delta t}{\Delta y}\right)^{2}&0\\ 0&\left(\frac{c_{0}\Delta t}{\Delta y}\right)^{2}&0\\ 0&0&0\end{pmatrix}\]

The homogeneous Neumann boundary conditions were enforced via ghost cells (see e.g. [19]) and implemented by using pooling layers, and manual adjustments according to the ghost cells.

## Appendix B Network Architecture

A quantitative summary of the presented network from figure 9 is provided in table 3 for the two-dimensional case and table 4 for the three-dimensional case. The choice of network hyperparameters varies slightly for the different approaches and cases. The interested reader is referred to the code made available at [71].
\begin{table}
\begin{tabular}{l l l} \hline
**layer** & **shape after layer** & **learnable parameters** \\ \hline input & \(256\times 3\times 3\times 3\) & 0 \\ \hline \hline upsample & \(256\times 6\times 6\times 6\) & 0 \\ \hline
3D convolution \& PReLU & \(256\times 6\times 6\times 6\) & \(1\,769\,728+1\) \\ \hline
3D convolution \& PReLU & \(256\times 6\times 6\times 6\) & \(1\,769\,728+1\) \\ \hline \hline upsample & \(256\times 12\times 12\times 12\) & 0 \\ \hline
3D convolution \& PReLU & \(128\times 12\times 12\times 12\) & \(884\,864+1\) \\ \hline
3D convolution \& PReLU & \(128\times 12\times 12\times 12\) & \(442\,496+1\) \\ \hline \hline upsample & \(128\times 24\times 24\times 24\) & 0 \\ \hline
3D convolution \& PReLU & \(128\times 24\times 24\times 24\) & \(442\,496+1\) \\ \hline
3D convolution \& PReLU & \(128\times 24\times 24\times 24\) & \(442\,496+1\) \\ \hline \hline upsample & \(128\times 48\times 48\times 48\) & 0 \\ \hline
3D convolution \& PReLU & \(64\times 48\times 48\times 48\) & \(221\,248+1\) \\ \hline
3D convolution \& PReLU & \(64\times 48\times 48\times 48\) & \(110\,656+1\) \\ \hline \hline upsample & \(64\times 96\times 96\times 96\) & \(0\) \\ \hline
3D convolution \& PReLU & \(64\times 96\times 96\times 96\) & \(110\,656+1\) \\ \hline \hline
3D convolution without padding & \(1\times 94\times 94\times 94\) & \(1\,729+1\) \\ \& adaptive Sigmoid & \\ \hline \end{tabular}
\end{table}
Table 4: Three-dimensional convolutional neural network architecture for the prediction of the indicator \(\gamma\). The number of filters per layer are double compared to the two-dimensional convolutional neural network. The total number of parameters is \(6\,306\,764\).

\begin{table}
\begin{tabular}{l l l} \hline
**layer** & **shape after layer** & **learnable parameters** \\ \hline input & \(128\times 8\times 4\) & 0 \\ \hline \hline upsample & \(128\times 16\times 8\) & 0 \\ \hline
2D convolution \& PReLU & \(128\times 16\times 8\) & \(147\,584+1\) \\ \hline
2D convolution \& PReLU & \(128\times 16\times 8\) & \(147\,584+1\) \\ \hline \hline upsample & \(128\times 32\times 16\) & 0 \\ \hline
2D convolution \& PReLU & \(64\times 32\times 16\) & \(73\,792+1\) \\ \hline
2D convolution \& PReLU & \(64\times 32\times 16\) & \(36\,928+1\) \\ \hline \hline upsample & \(64\times 64\times 32\) & 0 \\ \hline
2D convolution \& PReLU & \(64\times 64\times 32\) & \(36\,928+1\) \\ \hline
2D convolution \& PReLU & \(64\times 64\times 32\) & \(36\,928+1\) \\ \hline \hline upsample & \(64\times 128\times 64\) & 0 \\ \hline
2D convolution \& PReLU & \(32\times 128\times 64\) & \(18\,464+1\) \\ \hline
2D convolution \& PReLU & \(32\times 128\times 64\) & \(9\,248+1\) \\ \hline \hline upsample & \(32\times 256\times 128\) & 0 \\ \hline
2D convolution \& PReLU & \(32\times 256\times 128\) & \(9\,248+1\) \\ \hline
2D convolution \& PReLU & \(32\times 256\times 128\) & \(9\,248+1\) \\ \hline \hline
2D convolution without padding & \(1\times 254\times 126\) & \(289+1\) \\ \& adaptive Sigmoid & \\ \hline \end{tabular}
\end{table}
Table 3: Two-dimensional convolutional neural network architecture for the prediction of the indicator \(\gamma\). The total number of parameters is \(526\,252\).
## Appendix C Sharpness Quantification

The sharpness of the identified material distributions can be quantified by the \(L_{2}\)-norm of the spatial gradients of the recovered material, inspired by Sobel filters commonly used in image processing [72]:

\[||\nabla\hat{\gamma}(\mathbf{x})||_{2}=\sqrt{\nabla_{x}\hat{\gamma}(\mathbf{x})+\nabla_{y }\hat{\gamma}(\mathbf{x})} \tag{27}\]

A greater norm indicates a sharper edge. To highlight the improved sharpness recovery by the neural network Ansatz compared to the constant Ansatz, the norm from equation (27) around the circular void is illustrated in figure 21 for the two-dimensional case. At the cirular boundary, the norm is mostly larger in the neural network case showcasing a sharper inversion, which can be quantified by considering the mean of the non-zero values, which are extracted by selecting the points with a norm greater than a positive threshold, with the purpose of excluding the noise present in the constant Ansatz solution. For any threshold, the mean is larger for the neural network case, confirming an overall sharper recovered material distribution.

Figure 21: The \(L_{2}\)-norm of the spatial gradients of the recovered material distribution \(||\nabla\hat{\gamma}(\mathbf{x})||_{2}\).