Neural network augmented wave-equation simulation
Ali Siahkoohi Mathias Louboutin Felix J. Herrmann

A fifty-fold improvement of thermal noise limited inertial sensitivity by operating at cryogenic temperatures
    J.V. van Heijningen^†,
    December 30, 2023
=============================================================================================================


Accurate forward modeling is important for solving inverse problems. An
inaccurate wave-equation simulation, as a forward operator, will offset
the results obtained via inversion. In this work, we consider the case
where we deal with incomplete physics. One proxy of incomplete physics
is an inaccurate discretization of Laplacian in simulation of wave
equation via finite-difference method. We exploit intrinsic one-to-one
similarities between timestepping algorithm with Convolutional Neural
Networks (CNNs), and propose to intersperse CNNs between low-fidelity
timesteps. Augmenting neural networks with low-fidelity timestepping
algorithms may allow us to take large timesteps while limiting the
numerical dispersion artifacts. While simulating the wave-equation with
low-fidelity timestepping algorithm, by correcting the wavefield several
time during propagation, we hope to limit the numerical dispersion
artifact introduced by a poor discretization of the Laplacian. As a
proof of concept, we demonstrate this principle by correcting for
numerical dispersion by keeping the velocity model fixed, and varying
the source locations to generate training and testing pairs for our
supervised learning algorithm.




§ INTRODUCTION


In inverse problem, we heavily rely on having an accurate forward
modeling operators. Often, we can not afford being physically or
numerically accurate. In other words, being numerically inaccurate can
be due to computationally complexity of accurate methods, or incomplete
knowledge of the underlying data generation process. In either case,
motivated by <cit.>, we propose to intersperse CNNs
between timestepping algorithm for simulating acoustic wave-equation. We
mimic incomplete/inaccurate physics by simulating wave equation with
finite-difference method, while utilizing a poor (second-order)
discretization of Laplacian.

Conventional method for solving partial differential equations (PDEs),
e.g., finite-difference and finite-element method, given enough
computational recourses, are able to simulate high-fidelity solutions to
PDEs. On one hand, as long as the Courant–Friedrichs–Lewy conditions
for stability are satisfied, finite-difference methods are able to
compute solutions to PDE, regardless of medium parameters, with
arbitrary precision. On the other hand, finite-element method requires
careful meshing of the medium in order to carry out the simulation.
Another motivation behind this work is to exploit the fact that in
seismic, wave simulations are usually carried out for specific families
of velocity models and source/receiver distributions. We hope our
proposed method meets halfway between two mentioned extremes—i.e.,
being too generic (finite-difference method) and being too problem
specific (finite-element method).

There are several attempts to exploit learning methods in wave-equation
simulation. <cit.> approximates the solution to a
nonlinear PDE with a neural network. The neural network, given points on
the computational grid as input, computes the solution of PDE. Training
data is obtained by computing the solution of PDE on several points.
<cit.> completely ignore the Laplacian and they solely
rely on predicting the next timestep from the previous two timesteps by
learning the action of the spatially varying velocity and Laplacian.
While possible in principle, their approach needs to train for long
times to provide reasonable simulations on relatively simple models.
<cit.> instead of ignoring the physics, relies on
low-fidelity wave-equation simulation, and by exploiting transfer
learning, they utilize a single CNN to correct wavefield snapshots
simulated on a “nearby” velocity mode for numerical dispersion at any
given timestep. In this work, we extend ideas in
<cit.> and propose using multiple CNNs,
interspersed between low-fidelity timesteps. Finally,
<cit.> propose interspersing Krylov-subspace
iterations and neural nets while inverting the Helmholtz equation. They
show improvement in convergence by “propagating” an approximated
wavefield, obtained from a limited number of iterations, with the aid of
a trained convolutional neural net. This technique can be seen as the
frequency-domain counterpart of our proposed method.

Our paper is organized as follows. First, we describe our approach in
detail by first, describing our formulation for learned wave simulation.
Next, we introduce our training objective function. Due to dependencies
of CNN parameters, we devised an training heuristic that we describe.
Before explaining our numerical experiments, we state used CNN
architecture and training details. Next, we describe our three numerical
experiments we conduct and discuss effectiveness of the proposed method.



§ THEORY


We describe how we augment low-fidelity physics with learning techniques
to handle incomplete and/or inaccurate physics, where the low-fidelity
physics is modeled via finite-difference method with a poor
discretization of the Laplacian. To ensure accuracy, the temporal and
spatial discretization in high-fidelity wave-equation simulations have
to be chosen very fine, typically one to two orders of magnitude smaller
than Nyquist sampling rate. As mentioned earlier, we will utilize a poor
discretization (only second order) for the Laplacian to carry out
low-fidelity wave-equation simulations, but the scheme can be extended
to other proxies of incomplete or inaccurate physics.



 §.§ Simulations by
timestepping


After discretization of the acoustic wave equation, a single timestep of
of scalar wavefields, simulated on 0 ≤ t ≤ T, can be written as
below:


    𝐮_j+1 = 2𝐮_j-𝐮_j-1+δ t^2𝐜^2 Δ𝐮_j,    j=0,1,  … , N-1 ,


 where 𝐮_j is the high-fidelity scalar wavefield at
j^th timestep, δ t is the temporal discretization
(timestep), 𝐜 is the spatially varying velocity in the
medium, and Δ is the high-order discretization of
Laplacian. Similar to Equation <ref>, the low-fidelity
timestepping equation can be formulated as


    𝐮̅_j+1 = 2𝐮̅_j-𝐮̅_j-1+δ T^2𝐜̅^2 Δ̅𝐮̅_j,   j=0, 1, … , M-1,


 where 𝐮̅_j is the low-fidelity scalar wavefield,
δ T is the coarse timestep, 𝐜̅ is the coarse
spatially varying velocity, and Δ̅ is the coarse
(only second second order) discretized Laplacian.

Motivated by <cit.>, we consider every timestep as a
single layer in a neural network, where the discretized Laplacian is a
linear operator, followed by the (nonlinear) action of the spatial
varying velocity. Moreover, noticing the additional terms in the
Equation <ref>, each timestep is similar to a residual
block introduced by <cit.>. Figures <ref>
and <ref> schematically indicate each timestep as a block,
corresponding to high- and low-fidelity discretization of wave equation.
respectively. The similarity of high- and low -fidelity timestepping
method and CNNs can be perceived from Figures <ref>
and low-fidelity-step, respectively, where red and yellow blocks
correspond to high- and low-fidelity timestepping equations,
respectively.





As it can be seen, high-fidelity simulation of the wave equation up to
time t=T requires a lot of high-fidelity timesteps. On the other hand,
Figure <ref> shows that the low-fidelity simulations
can be done with much less low-fidelity timesteps, due to course time
sampling, which each timestep is cheaper than the high-fidelity
timesteps due to the coarse discretization of Laplacian. Although
computationally cheap, the low-fidelity wave-equation simulations suffer
from numerical dispersion artifacts.



 §.§ Learned wave simulations


Depending on the domain of application, we can assume wave simulations
are typically carried out for specific families of velocity models and
source/receiver distributions. This motivates us to deploy a data-driven
wave simulation algorithm which is coupled with low-fidelity and cheap
physics and hope to recover high-fidelity wave simulations on a family
of velocity models. In our method, we propagate the coarse-grained
wavefields according to Equation <ref> with a coarsened
Laplacian. After k timesteps, where k is a hyperparameter, we apply
a correction with a CNN, 𝒢_θ_i, parameterized by
θ_i, to the obtained wavefield at j^th timestep and
proceed with the timestepping. The proposed data-driven timestepping
wave simulation method is formalized in Equation <ref>.


    𝐮̅_j+1 = 𝒢_θ_i( 2𝐮̅_j-𝐮̅_j-1+δ T^2 𝐜̅^2 Δ̅𝐮̅_j),  i = ⌊j/k⌋    if  j≡ k-1  ( k),   
     
    2𝐮̅_j-𝐮̅_j-1+δ T^2 𝐜̅^2 Δ̅𝐮̅_j     else


 where j=0,1, … , M-1. The schematic representation of
Equation <ref> is illustrated in Figure <ref>. Yellow
blocks represent low-fidelity timesteps (see Equation <ref>
and Figure <ref>) and blue blocks correspond to CNNs,
𝒢_θ_i,  i = 0, 1, … , ⌊M-1/k⌋.



The CNNs correct for the effects of inaccurate physics—i.e., numerical
dispersion in our experiments, at every k^th low-fidelity
timestep. In this work, the parameters, θ_i, are not shared among
the CNNs. We have not explored the possibility of shared weights among
the CNNs in different stages of wave propagation.

Note that although parameters of the CNNs are not shared, they are not
independent—i.e., after j^th timestep, CNN
𝒢_θ_i, i = ⌊j/k⌋, corrects for
errors in the wavefield introduced by low-fidelity timestepping and
imperfections present in the output of (i-1)^th CNN, which
have been propagated through timestepping. Therefore, a small
perturbation in the parameters of a CNN in the initial stages of neural
network augmented timestepping causes noticeable differences in the
input of the CNNs in later stages of the wave propagation. The described
dependencies among the CNN parameters introduces difficulties in
optimizing the parameters of the CNNs. Below we describe our heuristic
for training the CNNs.



 §.§ Training objective


During training, we train all the CNNs toward the high-fidelity solution
of wave-equation, at the corresponding timestep, obtained by solving
Equation <ref>. As it can be seen from
Equation <ref>, after j^th timestep, CNN
𝒢_θ_i, i = ⌊j/k⌋, is tasked to
correct the effects of low-fidelity timestepping. During training,
i^th CNN maps its input,
2𝐮̅_j-𝐮̅_j-1+δ T^2 𝐜̅^2 Δ̅𝐮̅_j,
to 𝐮_j+1, result of j^th timestep using
high-fidelity timestepping, obtained by Equation <ref>.

Define function F̅_k(.) as the action of k low-fidelity
timesteps—i.e., F̅_k(.) represents k consecutive low-fidelity
time stepping blocks, depicted in Figure <ref>. Clearly,
F̅_k is only a function of k, δ T, 𝐜̅,
and Δ̅. Using the defined notation, we can write
the input to i^th CNN,
û_i,  i =0, 1, … , ⌊M-1/k⌋,
as follows:


    û_i = F̅_k (𝒢_θ_i-1(û_i-1 )),    i = 1, 2, … , ⌊M-1/k⌋, 
        û_0 = F̅_k(𝐪),


 where 𝐪 is the source. Also let 𝐮_τ_i denote
the wavefield obtained at j^th timestep of high-fidelity
timestepping (Equation <ref>), where
τ_i = j = (k+1)i-1.

The input-output pair of the i^th CNN is
(û_i, 𝐮_τ_i). We can generated multiple
training pairs for CNNs by simulating
(û_i, 𝐮_τ_i) pairs, for various velocity
models and source locations. Assume we have n pairs of of training
data for CNNs, namely,
(û_i^(p), 𝐮_τ_i^(p)),  p = 0,1, … , n-1.
The objective function for optimizing i^th CNN can be written
as follows:


    ℒ_i = 1/n∑_p=0^n-1𝒢_θ_i(û_i^(p)) - 𝐮_τ_i^(p)_1,    i = 0, 1, … , ⌊M-1/k⌋ .


 In the past, in a similar attempt, we used Generative Adversarial
Networks <cit.> to train a CNN in order to
remove numerical dispersion from wavefield snapshots
<cit.>. In this work we choose
to use ℓ_1 norm as the misfit function based on two reasons. First,
training GANs is computationally expensive since it requires training an
additional neural network that discerns between high-fidelity wavefield
snapshots and corrected ones. The computational complexity of the
proposed method in this work is significantly higher than our previous
attempts <cit.>, because it
involves training multiple CNNs. Based on the mentioned facts, for
limiting the computation time we chose to use ℓ_1 misfit function.
Second, motivated by a numerical experiment performed by
<cit.>, ℓ_1 norm misfit function yields
the second best results after the misfit function utilizing a
combination of GANs and ℓ_1 norm misfit. In the next section, we
describe our heuristic for training the CNNs.



 §.§ Training heuristic


To overcome complexities caused by dependencies between parameters of
CNNs, we optimize the objective functions ℒ_i with a
heuristic described below. We minimize
ℒ_i,  i = 0, 1, … , ⌊M-1/k⌋,
with respect to θ_i, respectively. In other words, we minimize
ℒ_i with respect to θ_i, by keeping the rest of the
parameters fixed. We keep updating all the set of parameters, in a
cyclic fashion—i.e., once we updated all the parameters,
ℒ_i,  i = 0, 1, … , ⌊M-1/k⌋, we
start over and update them again, in order, until a stopping criteria is
achieved. We will describe the stopping criteria used in our experiments
later.

We minimize objective functions the above objectives <ref>
with a variant of Stochastic Gradient Descent known as the Adam
optimizer <cit.> with momentum parameter β = 0.9
and a linearly decaying stepsize with initial value
μ = 2 × 10^-4 for both the generator and discriminator
networks. During each iteration of Adam, the gradient ℒ_i is
approximated by a single randomly selected training pair. These pairs
are selected without replacement. Once all the training pairs have been
selected, we start over by randomly picking training pairs, without
replacement from the entire training set.

The optimization carries out for a predetermined number of total
iterations, where each iterations consists of drawing a random
training pair, without replacement, and updating parameters of a CNN.
Additionally, while optimizing θ_i by keeping the rest of the
parameters fixed, before proceeding to the next set of parameters, we
carry out the optimization to update θ_i for number of
iterations, which we refer to it as mini-iterations.
Algorithm <ref> indicates the steps for optimizing objective
functions <ref>.


  1. INPUT: 
MaxItr  
MaxMiniItr  
F̅_k(.) 
𝐪^(p) p = 0,1, … , n-1 
𝐮_τ_i^(p),  p = 0,1, … , n-1,  i = 0, 1, … , ⌊M-1/k⌋ 
θ^0_i,  i = 0, 1, … , ⌊M-1/k⌋ 
2. ItrNum ← 0
3. FOR i = 0:⌊M-1/k⌋ DO
4.      θ_i = θ^0_i
5. FOR p = 0:n-1 DO
6.      û_0^(p) = F̅_k(𝐪^(p))
7. WHILE itrNum < MaxItr DO
8.      FOR i = 0:⌊M-1/k⌋ DO
9.          IF i>0 DO
10.             FOR p = 0:n-1 DO
11.                 û_i^(p) = F̅_k(𝒢_θ_i-1(û_i-1^(p)))
12.         FOR miniItrNum = 1:MaxMiniItr DO
13.             p ← SampleWithoutReplacement({ 0,1, … , n-1 })
14.             θ_i ←min_θ_i𝒢_θ_i(û_i^(p)) - 𝐮_τ_i^(p)_1
15.             ItrNum ← ItrNum + 1
16. RETURN θ_i,  i = 0, 1, … , ⌊M-1/k⌋
Heuristic for optimizing CNNs
𝒢_θ_i,  i = 0, 1, … , ⌊M-1/k⌋.




 §.§ CNN architecture


Motivated by our previous attempts for numerical dispersion removal from
wavefield snapshots <cit.>, we
use the exact architecture provided by <cit.>,
which includes Residual Blocks, the main building block of ResNets,
introduced by <cit.>, for all the CNNs
𝒢_θ_i,  i = 0, 1, … , ⌊M-1/k⌋.



 §.§ Training details and
implementation


While CNNs are known to generalize well—i.e., maintain the quality of
performance when applied to unseen data, they can only be successfully
applied to a data set drawn from the same distribution as the training
data. Because of the Earth's heterogeneity and complex geological
structures present in realistic-looking models, training a neural
network that can generalize well when applied to another velocity model
can become challenging. While we have successfully demonstrated that
transfer learning <cit.> can be used in
situations where the neural network is initially trained on data from a
proximal survey <cit.>, we chose in this
contribution, as a proof of concept, to keep the velocity model fixed,
and vary the source locations to generate different training/testing
pairs.

We use the Marmousi velocity model and out of 401 available shot
locations with 7.5 m spacing, we allocate half of the shot
locations to training and use the rest of the shot locations to generate
testing pairs, for evaluation purposes. The maximum simulation time in
our experiments in 1.1 s.

We designed and implemented our deep architectures in
TensorFlow[<https://www.tensorflow.org/>]. To carry out our
wave-equation simulations with finite differences, we used
Devito[<https://www.devitoproject.org/>]
<cit.>. We used the functionality of
Operator Discretization Library[<https://odlgroup.github.io/odl/>]
to wrap Devito operators into a TensorFlow layers. Our implementation
can be found on GitHub[<https://github.com/alisiahkoohi/NN-augmented-wave-sim>].

We ran our algorithm on Amazon Web Services' g3.4xlarge instance, where
we optimize the CNN parameters on a NVIDIA Tesla M60 GPU and Devito
utilizes 16 CPU cores to perform finite-difference wave-equation
simulations. Initially, we simulate the high-fidelity training wavefield
snapshots,
𝐮_τ_i^(p),  p = 0,1, … , n-1,  i = 0, 1, … , ⌊M-1/k⌋,
only once, in the beginning, and store them. In order to limit CPU-GPU
communication, before utilizing the GPU to to update
θ_i,  i = 0, 1, … , ⌊M-1/k⌋, we
generate the input to i^th CNN,
û_i^(p),  p = 0,1, … , n-1 all at once, and
store them. Afterwards, i^th CNN can be (re)trained using the
stored input/output wavefield snapshot pairs for several
mini-iterations.



§ NUMERICAL EXPERIMENTS


We want to indicate that neural networks, when augmented with inaccurate
physics, e.g., a poor discretization of Laplacian, are able to
approximate the wavefields obtained by an accurate approximation to wave
equation. To demonstrate this, we conduct three numerical experiments in
which we keep the velocity model fixed, and vary the source locations to
generate different training/testing pairs. The experiments differ in the
number of CCNs used throughout learned wave propagation. We use three,
five, and ten CNNs while keeping the total number of iterations fixed.
This implies that an experiment with more CNNs, optimizes each CNN with
a smaller number of iterations per CNN, because,
iterations per CNN×number of CNNs = total number of iterations.

A neural network augmented wave simulator with n_1 CNNs needs more
training iterations and possibly more training data to perform equally
as well as a neural network augmented wave simulator with n_2 CNNs,
when n_1 > n_2. For a fixed number of total iterations, iterations per
CNN is inversely proportional to number of CNNs utilized. Therefore, the
first n_2 CNNs in the neural network augmented wave simulator with
n_1 CNNs will perform worse than the CNNs in the wave propagator with
n_2 CNNs. Consequently, the error accumulated by the poor performance
of first n_2 CNNs, combined with artifacts introduced by low-fidelity
timestepping makes the matters worse for the later CNNs in the more
complex learned wave propagator. Therefore, in our experiments, since
the total number of iterations is fixed, we expect to see the quality of
dispersion removal degrade as the number of CNNs increase in a learned
wave propagator. Table <ref> summarizes the total
number of iterations, iterations per CNN, training pairs, training time,
and number of tunable parameters for the three different experiments.



As described earlier and presented in Table <ref>, the
 variable used in the  condition in line
8 of Algorithm <ref> is set to 500 for all our experiments.
Figures <ref> - <ref> demonstrate the values of
objective function presented in Equation <ref> in orange, and
the wavefield correction signal-to-noise ratio (SNR) in blue, evaluated
on testing data pairs during training, for experiments with three, five,
and ten CNNs, respectively. Note that the SNR curves have not been used
to determine when to stop training and they are only depicted for
demonstration purposes.

Figures <ref>, <ref>, and <ref>
show the wavefield correction SNR for first, second, and third CNN,
respectively, in the neural network augmented wave simulator that
includes three CNNs. Similarly, Figures <ref>,
<ref>, and <ref> depict the training
objective values throughout training for first, second, and third CNNs.
As it can be seen from objective function curves, the raining heuristic
has been effective and the objective function values have a decreasing
trend. Note that CNNs has been trained for 33500 iterations, on
average, with a total of 100500 iterations. Several equispaced spikes
can be noticed on the objective function value curves. For instance, see
the objective value function curve of the third CNN, in
Figure <ref>, at 6030, 8040, 10050, and  12060
iterations. The mini-batch we use in this experiment is 10. Those
spikes occur in moments in training when we have started retraining the
third CNN, after updating the first and second CNNs. As discussed
before, a change in the parameters of the CNNs preceding a CNN causes
changes in the input of the later CNN, and consequently the objective
function becomes large when starting to retrain the CNNs in later stages
again.

Similar objective function value and SNR curves for other two neural
network augmented wave propagators, utilizing five and ten CNNs, can be
found in Figures <ref> and <ref>, respectively.
First column in Figures <ref> and <ref>, from top
to bottom, indicate SNR of wavefield correction obtained by the first to
the last CNN, evaluated on testing pairs while training, respectively.
The second column of Figures <ref> and <ref>
indicate the objective function value curves throughout optimization of
Equation <ref> for training neural network augmented wave
propagators, utilizing five and ten CNNs, respectively. In both columns,
from top to bottom, the objective function value curves correspond to
CNNs from beginning to the end of the learned wave propagators, in
order.

We make two main observations from Figures <ref>
and <ref>. First, the objective function values indicate
overall decreases, validating effectiveness of the introduced heuristic.
Also, the spikes on the objective function value curve can be seen,
which are correlated with the stages in training when
Algorithm <ref> revisits a CNN after updating the rest of the
CNN parameters. As explained before, spikes are caused by change in
parameters of preceding CNNs to a CNN, which in turn alters the input
training wavefields of the CNN. Second, due to decrease in iterations
per CNN as the number of CNNs increases, the SNR curves converge to a
lower value when the number of CNNs increase.







Next, we will demonstrate the the corrected wavefields in three
conducted experiments evaluated over one testing shot location. For each
experiment, we show the high-fidelity wavefield snapshots,
𝐮_τ_i,  i = 0, 1, … , ⌊M-1/k⌋,
where i iterates over the CNNs, numerically dispersed low-fidelity
wavefields, and the corrected wavefield snapshots by the CNNs. To
evaluate the performance of each correction, we also depict the
correction error—i.e., difference between the high-fidelity and
corrected wavefield snapshots. Figure <ref> shows the mentioned
wavefield snapshots for the neural network augmented wave simulator with
three CNNs. First column shows the high-fidelity wavefields by solving
Equation <ref>, second column depicts low-fidelity
simulations by solving Equation <ref>, third column
indicates the result of neural network augmented wavefield simulations,
and the fourth column is the learned wave simulation error—i.e.,
difference between the first and last column in Figure <ref>.
Similarly, Figures <ref> - <ref> show the high- and
low- fidelity and learned wave simulation wavefield snapshots in the
first three columns, in order, for the neural network augmented wave
simulator with five and ten CNNs, respectively.

As expected because of the reasons stated before, we observe that the
quality of neural network augmented wave-equation simulation degrades as
the number of CNNs increases. On the other hands, the high quality of
learned wave simulation with few CNNs (see Figure <ref>) suggests
the quality of the simulation with more CNNs might be improved by
increasing the number of iterations. As it can be seen in the last
column of Figures <ref> and <ref>, the learned wave
simulation with ten CNNs has the lowest quality. It can be seen that the
learned wave simulation has the least accuracy in direct wave, which
happen to be the events with largest amplitudes. Also, it appears that
most of the numerical dispersion has been removed, the phase has been
recovered, and residual is mostly amplitude differences.











 §.§ Performance comparison: Single CNN low-to-high-fidelity
mapping


In order to evaluate the effectiveness of the proposed method we also
train a single CNN similar to our previous attempt to remove numerical
dispersion from wavefield snapshots
<cit.> and compare the result
of numerical dispersion removal with the proposed method. To be more
precise, for each presented neural network augmented wave-equation
simulation experiment, where we use three, five, and ten CNNs, we train
a single CNN, 𝒢_θ, with the same architecture as the
architecture used in learned wave propagators, in order to remove
numerical dispersion from all the low-fidelity wavefield snapshots
simulated by solving Equation <ref> for
j≡ k-1  ( k), on training shot locations. Likewise to
previous examples, here we also use half of the available shot locations
to simulate training pairs, and the rest is used to evaluate the
performance of the trained CNN. The input to 𝒢_θ
during training can be written as follows (compare with
Equation <ref>):


    ũ_i = F̅_k (ũ_i-1 ),    i = 1, 2, … , ⌊M-1/k⌋, 
        ũ_0 = F̅_k(𝐪),


 The desired output for the mentioned CNN is the high-fidelity wavefield
snapshots simulated on training shot locations,
𝐮_τ_i^(p),  p = 0,1, … , n-1,  i = 0, 1, … , ⌊M-1/k⌋.
The objective function for the mentioned CNN can be represented as
follows:


    ℒ = 1/n (⌊M-1/k⌋ + 1)∑_i=0^⌊M-1/k⌋∑_p=0^n-1𝒢_θ(ũ_i^(p)) - 𝐮_τ_i^(p)_1.


 We minimize objective function <ref> over θ with
Adam optimizes, using the same maximum number of iterations as before,
this time by combining all the training pairs associated with different
CNNs in the learned wavefield simulation example. As mentioned before,
in order to compare with the proposed method, we minimize objective
function <ref> over three different set of input-output
pairs, each corresponding to our presented experiments with varying
number of CNNs. Table <ref> summarizes the
total number of iterations, training pairs, training time, and number of
tunable parameters for three different cases, which differ in number of
timesteps which we choose to correct the numerical dispersion. This
selected timesteps are associated with the timesteps that the CNNs
operated on, in our three previous examples.



The slight difference in runtime among three different cases provided in
Table <ref> is partly due to different number
of training pairs needed to be generated. Figure <ref>
depicts the wavefield snapshot correction SNR curves, evaluated on
testing pairs while training, and the value of objective
function <ref>, in single CNN low-to-high-fidelity
mapping experiment, as a function of number of iterations.
Figures <ref>, <ref>,
and <ref> show the SNR curves, when we trained the CNN on
wavefield snapshots correspond to three, five, and ten, timesteps,
respectively. Similarly, Figures <ref>,
<ref>, and <ref> depict the training
objective function value (Equation <ref>), when the CNN
is trained on wavefield snapshots correspond to three, five, and ten,
timesteps, respectively. SNR curves depicted in first column of
Figure <ref> show the evolution of wavefield correction SNR
evaluated on randomly selected testing wavefield wavefields form the
wavefield snapshots combined from different timesteps. Therefore,
Figures <ref>,
<ref>, and <ref> indicate that the three
different CNNs converge to a wavefield correction SNR around 20 dB,
regardless of number of timesteps they are correcting for. Although this
does not suggest that the performance will stay the same as we increase
the number of timesteps needed to be corrected. By comparing
Figure <ref> with first column of Figure <ref>
(SNR curves for neural network augmented wave-equation simulation with
three CNNs), we observe that, on average the two methods are performing
equally well.



Finally, we will show the wavefield corrected by the single CNN
low-to-high-fidelity mapping method, for comparison with our proposed
method. Figures <ref> and <ref> indicate the
corrected wavefields, for cases where three and five timesteps need to
be corrected. Figures <ref> and <ref>
demonstrate the corrected wavefields, for the case where ten timesteps
need to be corrected, in two parts. In Figures <ref> -
<ref>, the first, second, third, and fourth columns depict
the high-fidelity and low fidelity wavefield snapshots, corrected
low-fidelity wavefield snapshots, and the error in numerical dispersion
removal, respectively. In each columns, from top to bottom, the
simulation time increases.

For comparison between our proposed method and single CNN
low-to-high-fidelity mapping, compare Figures <ref>
with <ref>, <ref> with <ref>, <ref>
with <ref>, <ref> with <ref>. As it
can be seen, the single CNN low-to-high-fidelity mapping method
maintains the quality of its performance when the number of timesteps
that need to be corrected increases. On the other hand, as the number of
CNNs in neural network augmented wave-equation simulation increases, the
performance drops, by keeping the maximum number of iterations fixed.
Also, by comparing Tables <ref>
and <ref>, we observe that the training time
needed for single CNN low-to-high-fidelity mapping, when number of
timesteps needed to be corrected increases, for fixed number of maximum
iterations, grows very slowly compared to the training time required for
neural network augmented wave-equation simulation.











§ CONCLUSIONS


Our numerical experiments demonstrate that, given suitable training
data, the well-trained neural network augmented wave-equation simulator
is capable of approximating wavefield snapshots simulated by
high-fidelity simulation. In this work, as a proxy of inaccurate
physics, we simulate wave-equation with finite-difference method, using
a poor discretization of Laplacian. Although not computationally
favorable to high-fidelity wave simulation, we showed that the learned
wave simulator deals with inaccurate physics. An important observation
we made is that training time of the proposed method gets quickly very
long and to achieve high accuracy, it may not be possible to utilize too
many CNNs. On the other hand, training time required for the single CNN
low-to-high-fidelity mapping experiments, conducted for the sake of
comparison, grows very slowly as the number of timesteps to be corrected
increases. In future, we intend to initialize the CNN parameters in the
proposed method with the parameters of a CNN trained by the single CNN
low-to-high-fidelity mapping algorithm. The initialization may
significantly reduce the training time needed for the neural network
augmented wave-equation simulation method, and may give the chance to
fine-tune the CNNs to the specific timestep that each CNN is assigned to
correct.



§ ACKNOWLEDGMENTS


The authors thank Xiaowei Hu for his open-access repository[<https://github.com/xhujoy/CycleGAN-tensorflow>]
on GitHub. Our software implementation built on this work.




