# Using explainability to design physics-aware CNNs for solving subsurface
inverse problems

J. Crocker *, K. Kumar  B. Cox

###### Abstract

We present a novel method of using explainability techniques to design
physics-aware neural networks. We demonstrate our approach by developing a
convolutional neural network (CNN) for solving an inverse problem for shallow
subsurface imaging. Although CNNs have gained popularity in recent years
across many fields, the development of CNNs remains an art, as there are no
clear guidelines regarding the selection of hyperparameters that will yield
the best network. While optimization algorithms may be used to select
hyperparameters automatically, these methods focus on developing networks with
high predictive accuracy while disregarding model explainability (descriptive
accuracy). However, the field of Explainable Artificial Intelligence (XAI)
addresses the absence of model explainability by providing tools that allow
developers to evaluate the internal logic of neural networks. In this study,
we use the explainability methods Score-CAM and Deep SHAP to select
hyperparameters, such as kernel sizes and network depth, to develop a physics-
aware CNN for shallow subsurface imaging. We begin with a relatively deep
Encoder-Decoder network, which uses surface wave dispersion images as inputs
and generates 2D shear wave velocity subsurface images as outputs. Through
model explanations, we ultimately find that a shallow CNN using two
convolutional layers with an atypical kernel size of 3$\times$1 yields
comparable predictive accuracy but with increased descriptive accuracy. We
also show that explainability methods can be used to evaluate the network’s
complexity and decision-making. We believe this method can be used to develop
neural networks with high predictive accuracy while also providing inherent
explainability.

###### keywords:

machine learning , CNN , explainability , XAI , subsurface imaging , MASW

††journal: Computers and Geotechnics

[1]The University of Texas at Austin, Department of Civil, Architectural, and
Environmental Engineering, Austin, TX, USA, 78712 [2]Utah State University,
Department of Civil and Environmental Engineering, Logan, UT, USA, 84322

##  1 Introduction

In recent years, non-invasive techniques like surface wave testing have become
increasingly popular for subsurface imaging due to their computational
efficiency and ease of field data acquisition. Surface wave techniques involve
collecting experimental wavefields recorded at the ground surface and solving
an inverse problem to identify a subsurface profile that matches the
characteristics of the recorded wavefield measurements. The inverse problem
involves iteratively modifying a 1D subsurface model and solving a forward
wave propagation problem (e.g., the Thomson-Haskell transfer matrix approach;
Haskell, 1953; Thomson, 1950) until a subsurface model that yields theoretical
data comparable to the experimental field data is found. This inverse problem,
referred to as surface wave inversion, is complex and challenging, as it is
non-unique and ill-posed (Cox and Teague, 2016; Foti et al., 2014, 2018;
Vantassel and Cox, 2022). Traditional surface wave inversions also lose the
true 2D/3D nature of the subsurface contained in the recorded waveforms by
only matching against theoretical 1D subsurface models. Although full-waveform
inversion (FWI) can provide true 2D/3D subsurface images, it is much more
time-consuming in terms of both field data acquisition and inversion compared
to traditional surface wave methods. Hence, there is a need for alternative
subsurface imaging methods to surface wave testing and FWI that are
computationally efficient and accurate. A promising alternative to traditional
surface wave inversion for developing subsurface images is using deep neural
networks, such as convolutional neural networks (CNNs). An appropriately
trained CNN can quickly and accurately generate subsurface images directly
from experimental data without needing to solve the forward problem. Although
the adoption of deep learning for seismic imaging is relatively recent
(Dramsch, 2020), many authors have successfully used CNNs to image the deep
subsurface (Li et al., 2019; Sun and Demanet, 2020; Wang et al., 2018; Wu et
al., 2018; Zheng et al., 2019). However, limited research has been conducted
on using CNNs for shallow subsurface imaging (depths less than 30 m). Shallow
subsurface imaging is particularly challenging because material properties
vary significantly over short distances (vertically and laterally) and the
various components of the elastic wavefield (i.e., compression, shear, and
surface waves) are mixed, having not yet propagated far enough to spread out
from one another. Vantassel et al. (2022) developed a CNN for shallow
subsurface imaging of soil-over-rock profiles that accepts active-source
seismic wavefields recorded by a linear array of receivers as input to rapidly
generate 2D subsurface images of shear wave velocity (Vs). Abbas et al. (2022)
expanded on this network by demonstrating the possibility of generalizing such
CNNs to a variety of field-testing configurations using a frequency-domain
approach rather than a time-domain approach. The input of the frequency-domain
CNN is dispersion images calculated from active-source wavefields recorded by
a linear array of receivers placed on the ground surface, exactly like those
collected from Multi-channel Analysis of Surface Waves (MASW) testing, as
illustrated in the top half of Fig. 1. While their CNN was trained on
synthetic data, it yields a high predictive accuracy on field data collected
at a site with corroborating subsurface ground truth. Despite the promising
results, CNNs are traditionally designed using a standard trial-and-error
approach by varying the network architecture and performing hyperparameter
tuning until the network yields a high predictive accuracy. This blind trial-
and-error approach ignores the logic, or explanations, behind model
predictions, limiting their trustworthiness and applicability. In this work,
we propose a framework for designing a physics-aware neural network
architecture using model explainability to select the network’s
hyperparameters, leading to a network that has high predictive accuracy and
inherent explainability.

Deep learning models vary in accuracy, complexity, and interpretability.
Simple models, such as linear regression, classification rules, and decision
trees, use linear functions that are highly interpretable, but lack predictive
accuracy. Highly non-smooth, non-linear, or multilayer models, such as CNNs,
trade this interpretability for high accuracy. CNNs are multilayer neural
network models capable of solving complex problems such as image recognition,
but offer no explanations regarding their internal logic or decision making
process (i.e., they are black box in nature). Hence, it is difficult to assess
the physics behind their predictions. Typically, models are designed to have
high predictive accuracy (i.e., the accuracy of the CNN’s predictions versus
known data) but low descriptive accuracy (i.e., how accurately the model
explains why it made that prediction; Murdoch et al., 2019). Therefore, model
explainers are required to address the low descriptive accuracy of these black
box models. An alternative to training black box models is to incorporate a
physics-informed loss function. Although Physics Informed Neural Networks
(PINNs) are trained to conserve physical properties, they do not guarantee a
physically-consistent response, and the reasoning behind their predictions
still remains a black box (Raissi et al., 2019). The field of Explainable
Artificial Intelligence (XAI) offers tools and techniques to comprehend the
results of deep learning models. These tools include model explainers such as
LIME (Ribeiro et al., 2016), class activation maps (CAMs; Zhou et al., 2016),
eXplanation with Ranked Area Integrals (XRAI; Kapishnikov et al., 2019), and
integrated gradients (IG; Sundararajan et al., 2017), which provide post hoc
explanations of model predictions without the need for retraining. These post
hoc explanations interpret the model’s decision-making process after its
development, thus improving its trust without compromising the model’s
predictive accuracy (Bien et al., 2018; Lauritsen et al., 2020; Lou et al.,
2019; Rajpurkar et al., 2017, 2018; Mitani et al., 2020; Zech et al., 2018).
However, the post hoc nature of these explainers also means the explanations
are not rigorous and are fallible, and may only capture a specific aspect of a
model (Linardatos et al., 2020). Moreover, model explanations are limited
mainly to understanding and interpreting predictions rather than designing and
optimizing the deep learning model, and thus are disregarded during the design
of neural networks.

Despite advances in neural network algorithms, the process of designing a
neural network remains an art. Developing a CNN requires careful tuning of
more than 100,000 hyperparameters, such as the kernel size, depth of the
network, type of layers, learning rate, or optimization algorithm, through a
trial-and-error iterative process of minimizing its predictive error. Without
a clear set of guidelines for selecting hyperparameters, model tuning uses a
brute-force approach to optimize the network until it yields a high predictive
accuracy without considering its descriptive accuracy. Approaches such as
Bayesian optimization or search algorithms offer the automatic selection of
hyperparameters in CNNs that yield a low error; however, these approaches also
disregard model explanations (Hinz et al., 2018; Lee et al., 2018; Loshchilov
and Hutter, 2016; Snoek et al., 2012; Wu et al., 2019). The design of CNNs by
minimizing the predictive error does not necessarily yield models that are
logically consistent. The network may not capture the physics or explain the
logic behind its predictions, potentially yielding untrustworthy results.
Although descriptive accuracy metrics may offer post hoc justification of
model behavior, it is often considered separately from the model’s predictive
accuracy, leading to a disconnect between model design and use. This
disconnect is particularly noticeable in the field of geophysics, as most
models are designed only using predictive accuracy while offering no
explanations (Abbas et al., 2022; Araya-Polo et al., 2018; Çaylak and Kaftan,
2014; Hu et al., 2021; Li et al., 2019; Lu et al., 2019; Vantassel et al.,
2022; Yang and Ma, 2019; Zhu et al., 2022).

Figure 1: The proposed framework for developing a convolution neural network
(CNN) for solving an inverse problem for shallow subsurface imaging. The CNN
takes surface wave dispersion images as inputs to generate 2D shear wave
velocity images. The development process begins with the training of a
proposed CNN. The network’s predictive and descriptive accuracies are
evaluated and hyperparameter tuning is performed accordingly until a CNN is
found that yields logical predictions while maintaining a high predictive
accuracy.

In this work, we address the disconnect between predictive and descriptive
accuracy by developing a framework for creating physics-aware neural networks
by integrating model explainability into the network development process. We
create a CNN for solving an inverse problem associated with shallow 2D
subsurface imaging, where model explainability is used to iteratively select
network hyperparameters, such as the depth of the network and the size of the
convolutional layers. An overview of our approach is shown in Fig. 1. The CNN
is trained on 80,000 synthetic subsurface images representing the common
geological condition of soil with varying thickness and stiffness overlying
undulating rock of varying stiffness. We use finite difference elastic wave
propagation simulations and MASW-style wavefield transformation methods to
obtain surface wave dispersion images that capture key features of each
subsurface model. These key features of the subsurface model are expected to
be best predicted by the high-power regions (i.e., red colors) in the
dispersion images, which represent the fundamental and potentially higher-
modes of surface wave propagation. Thus, our dataset consists of pairs of
known 2D subsurface velocity images and their corresponding surface wave
dispersion images. The network is first trained with respect to its predictive
accuracy, but we include an additional step that evaluates the descriptive
accuracy of the trained network on the validation dataset. By doing so, we can
interpret the CNN’s decision-making and determine if its predictions are
logical. If we determine the network makes illogical predictions, the
development process is revisited, and various hyperparameters are iteratively
tuned to improve the network.

In this work, we propose a framework for performing hyperparameter tuning and
developing a neural network architecture that improves its descriptive
accuracy while maintaining a high predictive accuracy. The aim of this work is
not to generalize the network to geologic conditions outside the bounds of the
training dataset or to improve the accuracy of its predicted 2D shear wave
velocity profiles; we instead guide the network to make logical decisions
based on the current problem. We believe that progress made in this area will
benefit the development of CNNs for solving inverse problems related to
subsurface imaging.

###  1.1 Convolutional neural networks for solving an inverse problem for
shallow subsurface imaging

Convolutional neural networks are a type of deep learning model that excels at
image classification and generation. CNNs are typically composed of one or
more convolutional layers, each followed by a pooling layer, and commonly
ending with a fully connected layer. Convolutional layers consist of a set of
filters defined by a kernel size. These filters are used to perform
convolution on the input, and an activation function is applied to obtain
feature maps, which contain information about the spatial features of the
input. Pooling layers are used to reduce the dimensionality of each feature
map representation, thus reducing the number of learnable parameters and
computational costs. Finally, after multiple cycles of convolution and pooling
layers, fully connected layers may be used to learn non-linear combinations of
the feature maps.

With the addition of each layer, CNNs improve their predictive accuracy as
each additional layer extracts more complex information from the input image.
For example, the first layer typically identifies the layer boundaries within
an image, and subsequent layers identify complex features such as gradients
and curvatures. Many different metrics, such as mean square error (MSE), mean
absolute error (MAE), or mean absolute percentage error (MAPE), are used to
evaluate the accuracy of the network’s predictions, and typically networks are
developed to have high predictive accuracy. An example of a CNN architecture
is provided as Fig. 2, which is the final network developed in this study, as
described in later sections. Note that we use only two 2D convolutional layers
separated by a single max pooling layer ending with a final fully connected
layer. For more information on CNNs, please refer to Goodfellow et al. (2016).

Figure 2: The architecture of a convolutional neural network used to generate
subsurface images using MASW-style data. The structure, kernel size of the
convolutional layers, and depth are selected considering the network’s
predictive and descriptive accuracy.

##  2 Explainability metrics

###  2.1 Predictive accuracy metrics

In this study, we define the global accuracy of our network’s predictions
based on the entire subsurface image rather than a specific feature (e.g.,
only the Vs of the overlying soil layer or depth to the rock layer).
Therefore, we measure the accuracy of our networks across 20,000 testing image
pairs using two global error metrics: (1) MAPE and (2) mean structural
similarity index (MSSIM; Wang et al., 2004). MAPE measures the mean of the
absolute value of the pixel-by-pixel percent error of each predicted Vs image
compared to the true image. We use MAPE to relate to the physical units of
each image (i.e., Vs magnitude), while other common metrics such as MSE or MAE
only correspond to normalized images. A limitation of MAPE is that it does not
consider the structure of an image. For example, if a true image and its
predicted counterpart are reordered pixel-by-pixel in the same way, then the
MAPE does not change despite the changes in images. Therefore, we also
consider the SSIM index. The SSIM index compares a predicted image to its true
counterpart using three criteria: luminance, contrast, and structure. This
index is computed on a pixel-by-pixel basis to create an SSIM index map, and a
mean value may be taken as the MSSIM to provide a single overall measure of
image similarity for each prediction. The MSSIM ranges from zero for
predictions with poor similarity to their respective true images to one for a
prediction image that is identical to the true image. In this study, we follow
the parameters recommended by Wang et al. (2004) to calculate MSSIM with a
dynamic range equal to the difference between the maximum and minimum Vs
values of the true subsurface profiles in the testing dataset.

We choose to use both MAPE and the MSSIM index to provide two perspectives
regarding the predictive accuracy of our networks. As MAPE provides an error
in physical units, we propose that this measures the overall accuracy of the
predicted Vs values in each subsurface image despite not considering the
overall image structure. In contrast, the MSSIM index considers the predicted
Vs values in combination with the image structure, and therefore we
hypothesize that it reflects the accuracy of the predicted layer boundary in
each subsurface image. We illustrate the use of both metrics using a subset of
eight true subsurface Vs images from the testing dataset (Fig. 4), and the
respective predictions of the Transpose Convolution CNN, as described in a
later section, with MAPE and MSSIM values (Fig. 4).

Figure 3: A selected set of true synthetic subsurface soil-over-rock Vs images
ranging from simple images (a-d) to complex images (e-h).

Figure 4: The Transpose Convolution network’s predictions of the subsurface
images depicted in Fig. 4.

Note that the images shown in Fig. 4 and Fig. 4 were specifically selected to
better illustrate the following trends, which are consistent for all the
images analyzed throughout this study. Comparing Figs. 4a and b to the true
images in Figs. 4a and b, we see the predictions are visually similar, which
is reflected by their low MAPE values and high MSSIM indices. Specifically,
the Vs values and the layer boundaries, which are simply the depths to the
rock layer, are predicted well in these images. In Figs. 4c and d, we see the
Vs values are generally predicted well (relatively low MAPEs), while the
undulations of the layer boundary are relatively smoothed (relatively low
MSSIM indices that indicate structural differences at the interface). In Figs.
4e and f, the predicted Vs values are higher than the true images (high
MAPEs), and the layer boundaries are overly smoothed, although the average
depths to the rock layer are captured (relatively low MSSIM indices). Finally,
in Figs. 4g and h, we see the predicted Vs values are generally lower than the
true images (high MAPEs) and undulations in the layer boundary are lost (Fig.
4g) or inverted (Fig. 4h; low MSSIM indices). From these results, we conclude
that the magnitude of MAPE is strongly affected by the overall predicted Vs
values, while the magnitude of the MSSIM index strongly correlates with the
prediction of layer boundaries, and thus we use both metrics to evaluate the
predictive accuracy of our networks. Note that the CNN is trained only using
the mean absolute error (MAE) loss function, as described later.

###  2.2 Descriptive accuracy metrics

We use post hoc explainability methods to comprehend the black box CNN model
used in this study. As stated previously, there is no single post hoc
explainability method that can cover all the complexities of a CNN, as they
are not rigorous. Therefore, we use two different descriptive accuracy metrics
to review the internal logic of our CNNs: Score-weighted Class Activation
Mapping (Score-CAM; Wang et al., 2020) and Deep SHAP (Lundberg and Lee, 2017).
We present a brief overview of these methods, but the reader is encouraged to
review the literature for a more in-depth discussion of the algorithms and
implementations.

####  2.2.1 Score-weighted Class Activation Mapping (Score-CAM)

Heatmaps are often used as a model explainer, as they provide a visual
representation of a network’s internal logic. One popular method is the
Gradient-weighted Class Activation Mapping (Grad-CAM; Selvaraju et al., 2017)
algorithm, which is a post hoc explainer that uses the gradients of a class,
or the classification of an image, flowing into a CNN’s final convolutional
layer to create heatmaps associated with that class. It is a generalization of
the method used to create class activation maps (CAMs) by Zhou et al. (2016),
as it does not require using a particular architecture (i.e., a final
convolutional layer followed by a global average pooling layer and fully
connected layer) to develop the heatmaps

Although a popular explainability method, recent work shows that Grad-CAM may
generate heatmaps that highlight irrelevant regions in the input image that
are not used by a network for prediction (Chattopadhay et al., 2018; Wang et
al., 2020). This is caused by the gradient averaging step, as the computed
weights used to create a heatmap do not adequately capture each feature map’s
importance. For example, two feature maps may have similar weights, despite
one feature map having a lower contribution to the output image, as
demonstrated by Wang et al. (2020). Hence, we use Score-CAM, proposed by Wang
et al. (2020), which uses score-based weights to obtain a final heatmap rather
than using gradients to weight feature maps. Score-CAM extracts the feature
maps ($A_{k}$) generated by the last convolutional layer in a network for a
given input image, where _k_ is the number of feature maps. Each feature map
is then treated as a mask on the original input image, and the masked image is
passed through the network to obtain its score of a target class. For image
generation applications, each pixel is assigned a class rather than an
individual class for the whole image. The scoring process is repeated _k_
times to obtain the score-based weights ($\alpha_{k}^{c}$) and the final
heatmap is generated using a linear combination of $\alpha_{k}^{c}$ with
$A_{k}$. Score-CAM uses ReLU to preserve positive values in the heatmap while
setting negative values to zero, so the heatmap only highlights parts of the
input image that contribute to the output. Additionally, the authors note that
while the last convolutional layer in a network is the preferable choice for
use with Score-CAM, this process can be applied to any convolutional layer to
obtain its associated heatmap. An example of a Score-CAM heatmap is shown in
Fig. 1 under the descriptive accuracy header, while its interpretation is
discussed in a later section.

####  2.2.2 Deep Learning Important FeaTures (DeepLIFT) with Shapley values
(Deep SHAP)

Shapley values are used to quantify the relative contributions of various
input features to an output. Lundberg and Lee (2017) used this concept to
create SHAP (SHapley Additive exPlanation) values, which quantify the relative
change in the expected model prediction caused by the introduction of each new
feature to the model. Deep SHAP (Lundberg and Lee, 2017) combines SHAP values
with the Deep Learning Important FeaTures method (DeepLIFT; Shrikumar et al.,
2017). The DeepLIFT method explains the difference in an output from a set
“reference” output by using the differences between perturbed inputs and a set
“reference” input. Specifically, the relative contribution caused by the
differences between neurons for an input image and its reference input must
equal the difference between an output and its reference output. These
relative contributions are separated into positive and negative components,
and a linear composition is used to assign contribution scores for each neuron
to its immediate inputs. Although similar to other backpropagation methods,
SHAP addresses the issue of saturation that is inherent in gradient-based
methods.

SHAP values provide a method of measuring additive feature importance, which
allows it to be combined with other additive feature attribution methods, such
as DeepLIFT. As DeepLIFT can linearize the non-linear components of a deep
neural network (i.e., through its linear composition rule) and its input
images are assumed to be independent of one another, DeepLIFT can approximate
SHAP values, leading to Deep SHAP. In this case, SHAP values measure the
relative contributions of perturbed inputs versus the unperturbed input. Like
Shapley values, the magnitude of SHAP values refers to the degree of influence
an input has on output, while positive and negative values correspond to
positive and negative contributions to the output. We compute Deep SHAP values
using the open-source package shap (Lundberg and Lee, 2017). An example Deep
SHAP result is shown in Fig. 1 under the descriptive accuracy header, while
its interpretation is discussed in a later section.

##  3 Dataset creation

Our dataset contains two components: (1) a set of 2D subsurface Vs images and
(2) a corresponding set of surface wave dispersion images (one dispersion
image for each 2D subsurface model, explained in the following sections). We
begin by developing 2D subsurface Vs images, then use these Vs images as the
basis for generating corresponding compression wave velocity (Vp) and mass
density images using common relationships. The Vs, Vp, and mass density images
together form the subsurface models needed for elastic wave propagation
simulations. Following subsurface model development, we use 2D finite
difference numerical wave propagation simulations with these subsurface models
to generate waveforms recorded by a linear array of receivers at the ground
surface. These waveforms are then transformed into dispersion images that
capture key features of each subsurface model. In this study, only the pairs
of subsurface Vs images and their corresponding dispersion images (refer to
Fig. 5) are considered for CNN development. We discuss additional details of
the dataset development in the following sections.

###  3.1 Development of subsurface Vs images

We design each subsurface model to mimic a relatively simple but common
subsurface geological condition: soil with varying thickness and stiffness
overlying undulating rock of varying stiffness (i.e., soil-over-rock with an
irregular interface). These subsurface models are programmatically generated
using an in-house code. We create the subsurface models using a range of
typical values for subsurface material characteristics, such as Vs, density,
and Poisson’s ratio, as per Foti et al. (2018). Additionally, we also vary the
undulations of the layer interface, the depth to the rock/weathered rock
layer, and the depth of the groundwater table. The development process is
summarized as follows: we generate a 104 m-wide and 24 m-deep subsurface Vs
image using a range of typical Vs values for the subsurface materials. We
separate the upper soil layer from the lower rock layer with a spatially
variable and undulating interface, which is created by summing three
sinusoidal waves, each controlled by a respective spatial undulation frequency
and amplitude pair. Three types of layer interfaces are considered: (1)
“linear”, which has nearly linear interfaces with very little to no
undulations, but may be inclined (Fig. 5a), (2) “slightly undulating”, which
has slight undulations (Fig. 5b), and (3) “highly undulating”, which has more
frequent undulations (Fig. 5c). Our dataset contains 10% linear, 60% slightly
undulating, and 30% highly undulating subsurface models, such that the network
is trained primarily on more realistic subsurface model geometries (i.e.,
linear and slightly undulating layer interfaces) while still accounting for
more extreme or complex layer interface geometries. We place the layer
interface at a variable depth between 5–12 m, which completes the base
subsurface Vs image. However, to create relatively more realistic subsurface
images, we use random sampling to introduce lateral and vertical perturbations
to each image.

We next create subsurface Vp images by relating each subsurface Vs image to Vp
using a Poisson’s ratio between 0.15–0.35 for the soil layers and between
0.20–0.25 for the rock layers. A groundwater table is randomly placed between
the surface and bottom of the Vs image, and any soil materials beneath the
water table are assigned a saturated Poisson’s ratio between 0.47–0.49 instead
of the default range noted above. Finally, we create subsurface density images
by assigning a value between 1,650–2,000 kg/m3 for soil and 2,100–2,400 kg/m3
for rock/weathered rock.

Following these procedures, we develop 100,000 synthetic subsurface models in
total, with 80,000 for training and 20,000 for testing. Note that we only
consider the subsurface Vs images in our dataset, but the Vs, Vp, and density
images combined create the subsurface models on which we perform the wave
propagation simulations described in the following section. To prepare the Vs
images for our networks, we normalize them by the maximum Vs value in the
training set. Additionally, the Vs images are trimmed to include only the 48 m
beneath the center of each image so that they only represent the subsurface
conditions directly beneath the linear array of receivers used to record the
wave propagation simulations, as per Abbas et al. (2022). Therefore, each
input 2D Vs image is 24 m deep and 48 m wide, corresponding to an array of Vs
values of size 24$\times$48$\times$1.

###  3.2 Development of dispersion images

Each dispersion image is calculated from an active-source seismic wavefield
recorded by a linear array of surface sensors. We develop seismic wavefields
by performing elastic wave propagation simulations with the 2D finite-
difference open-source software DENISE (Köhn, 2011; Köhn et al., 2012). For
brevity, we summarize our experimental configuration setup, which is similar
to Abbas et al. (2022). The surface array used to record the wavefields
consists of 48 receivers equally spaced at 1 m intervals, and the source used
in our simulations was a 30 Hz Ricker wavelet, which is meant to simulate a
broadband impact source like a sledgehammer or drop weight. Additionally, we
use the same simulation parameters (i.e., finite difference operator, boundary
conditions, etc.) and model two seconds of wave propagations through each of
the 100,000 models. These simulations are performed on the Texas Advanced
Computing Center’s (TACC’s) high-performance cluster Stampede2 using a single
compute node with Skylake processors.

We generate dispersion images from the seismic wavefields recorded at the
surface of each model by using the frequency-domain beamforming (FDBF)
wavefield transformation (Zywicki, 1999). Additional details about the
dispersion image computations may be found in Abbas et al. (2022). The FDBF
transformation was used to generate dispersion images over a common frequency
range of 5–80 Hz with a 1 Hz interval, resulting in 75 frequencies, and a
common phase velocity range of 50–1,250 m/s with a velocity interval of 3 m/s,
resulting in 400 velocities. Therefore, each dispersion image corresponds to
an array of phase velocity values of size 75$\times$400$\times$1\. To
normalize each dispersion image, we use frequency-dependent normalization,
which simplifies the images for the networks compared to using absolute
maximum normalization. These dispersion images are generated using the open-
source Python package _swprocess_ (Vantassel, 2021).

Examples of pairs of subsurface Vs images and their respective dispersion
input images are shown in Figs. 5a-c and Figs. 5d-f, respectively. As noted
above, and documented in many publications (e.g., Foti et al., 2018), surface
wave dispersion images capture key features of the subsurface (e.g., the
velocity and layering variations). These key features are expected to be best
predicted by the high-power regions (i.e., red colors) in the dispersion
images, which represent the fundamental and potentially higher-modes of
surface wave propagation that are excited by an active source operated on the
surface. In traditional surface wave testing, it is important to discern which
mode(s) are present in the dispersion image, such that appropriate, mode-
specific inversions can be performed. While mode jumps from fundamental to
higher-modes can be noted in some of the dispersion images (refer to Figs. 5d
and e), it is not necessary to isolate/delineate the fundamental-mode prior to
network training. Rather, the entire dispersion image is used with an
understanding that the network should ideally focus on learning what
subsurface characteristics tend to generate the high-power trends in the
dispersion image, regardless of the exact mode of surface wave propagation.
Thus, for simplicity herein, we do not try to distinguish which modes are
present in the dispersion images, and rather we simply refer to these
fundamental and possibly higher-mode regions as the _high-power modal trend_.

Figure 5: Examples of (a-c) subsurface Vs images with linear, slightly
undulating, and highly undulating layer boundaries, respectively, and (d-f)
their associated dispersion input images.

##  4 Developing a convolutional neural network using explainability

We now illustrate our proposed framework for developing a physics-aware CNN
using model explainability. Note that we focus on tuning a select group of
hyperparameters, such as the network’s decoder structure, the kernel size of
each convolution layer, and the depth of the network, as the network’s
performance was most sensitive to these hyperparameters.

###  4.1 Selecting a decoder structure

We begin with a simple Encoder-Decoder CNN architecture, as this is one of the
most widely used networks for solving inverse problems in imaging (Genzel et
al., 2022). An Encoder-Decoder CNN processes an input image $x\in
X\subset\mathbb{R}^{n}$ and maps it to a feature space $z\in
Z\subset\mathbb{R}^{d}$. The decoder produces an output $y\in
Y\subset\mathbb{R}^{l}$ from the feature map $z$. Traditionally, an Encoder-
Decoder CNN is symmetric with skipped connections. We develop an asymmetric
Encoder-Decoder CNN architecture to match the shape of the input dispersion
images (400$\times$76$\times$1) to the shape of the output Vs subsurface
images (24$\times$48$\times$1). The encoder contains five convolutional layers
using a conventional 3$\times$3 kernel size, each followed by a rectified
linear unit (ReLU) and a max pooling operation.

We evaluate two commonly used decoder structures: (1) a decoder structure
consisting of 2D Transpose Convolution and UpSampling layers and (2) a
Dense/fully connected layer. The first trial decoder structure used in the
”Transpose Convolution” network contains five 2D Transpose Convolution layers,
each followed by an UpSampling layer to mirror the proposed encoder structure.
We vary the kernel size for the transpose convolutional layers to match the
shape of the output Vs subsurface images to the input dispersion images. The
second trial decoder structure used in the ”Dense” network consists of a
single Dense layer containing 1152 nodes, which correspond to each pixel in
our flattened Vs subsurface image (i.e., 24$\times$48$\times$1 flattened).
Additionally, we choose a learning rate of 5E-4, batch size of 16, and 40
training epochs for the training of both networks. The Adam optimizer (Kingma
and Ba, 2014) and mean absolute error (MAE) are used as the optimizer and loss
function. Table 1 provides a detailed description of each network, and the
difference between the decoder structures is highlighted. We note that
validation was performed using 20% of the training dataset, and the loss
function for validation did not indicate any overfitting by the network.

Table 1: Detailed descriptions of the Transpose Convolution Network (left) and the Dense Network (right). All convolution, pooling, and upsampling layers listed are 2D. Note that the only difference between the networks is the decoder structure, which is highlighted. Transpose Convolution Network | Dense Network
---|---
Network Layer Type | Filters |  Kernel/Pool Size |  Output Shape |  Network Layer Type | Filters |  Kernel/Pool Size |  Output Shape
Input |  |  |  400$\times$76$\times$1 (Dispersion Image) |  Input |  |  |  400$\times$76$\times$1 (Dispersion Image)
Convolution | 32 |  (3, 3) |  398$\times$74$\times$32 |  Convolution | 32 |  (3, 3) |  398$\times$74$\times$32
Max Pooling |  |  (3, 1) |  132$\times$74$\times$32 |  Max Pooling |  |  (3, 1) |  132$\times$74$\times$32
Convolution | 32 |  (3, 3) |  130$\times$72$\times$32 |  Convolution | 32 |  (3, 3) |  130$\times$72$\times$32
Max Pooling |  |  (3, 1) |  43$\times$72$\times$32 |  Max Pooling |  |  (3, 1) |  43$\times$72$\times$32
Convolution | 64 |  (3, 3) |  41$\times$70$\times$64 |  Convolution | 64 |  (3, 3) |  41$\times$70$\times$64
Max Pooling |  |  (1, 3) |  41$\times$23$\times$64 |  Max Pooling |  |  (1, 3) |  41$\times$23$\times$64
Convolution | 128 |  (3, 3) |  39$\times$21$\times$128 |  Convolution | 128 |  (3, 3) |  39$\times$21$\times$128
Max Pooling |  |  (3, 3) |  13$\times$7$\times$128 |  Max Pooling |  |  (3, 3) |  13$\times$7$\times$128
Convolution | 128 |  (3, 3) |  11$\times$5$\times$128 |  Convolution | 128 |  (3, 3) |  11$\times$5$\times$128
Transpose Convolution | 128 |  (1, 1) |  11$\times$5$\times$128 |  |  |  |
UpSampling |  |  (1, 2) |  11$\times$10$\times$128 |  |  |  |
Transpose Convolution | 64 |  (1, 2) |  11$\times$11$\times$64 |  |  |  |
UpSampling |  |  (1, 2) |  11$\times$22$\times$64 |  Flatten |  |  |  7040
Transpose Convolution | 32 |  (1, 2) |  11$\times$23$\times$32 |  Dense |  |  |  1152 (Flattened Vs Image)
UpSampling |  |  (1, 1) |  11$\times$23$\times$32 |  |  |  |
Transpose Convolution | 32 |  (2, 2) |  12$\times$24$\times$32 |  |  |  |
UpSampling |  |  (2, 2) |  24$\times$48$\times$32 |  |  |  |
Convolution | 1 |  (1, 1) |  24$\times$48$\times$1 (Vs Image) |  |  |  |

We now review the predictive accuracy of both networks. Fig. 6 presents a set
of randomly selected subsurface Vs images from the testing dataset (left
column), the Transpose Convolution network’s predictions for these images
(middle column), and the Dense network’s predictions for these images (right
column). Note that while these images were randomly selected, they have been
plotted in order of decreasing layer boundary complexity for ease of
discussion and visualization.

Overall, across all 20,000 testing models, the Transpose Convolution network
has an MAPE of 11.3% and MSSIM index of 0.54. Recall that lower values of MAPE
correspond to more accurate Vs predictions, while higher values of MSSIM
correspond to more accurate layer boundary predictions. Given the complexity
of predicting Vs subsurface profiles from dispersion images, an MAPE of 11.3%
is considered adequate and only a few percentage points higher than the best
MAPEs predicted by Abbas et al. (2022) for their frequency-velocity CNN.
Additionally, a lower MSSIM index is expected; because Transpose Convolution
layers followed by UpSampling layers typically cause pixelation in predicted
images, the boundary between soil layers is poorly predicted, causing a lower
MSSIM index (compare this network’s prediction to the true image in row c of
Fig. 6).

Figure 6: Comparisons between five randomly selected subsurface shear-wave
velocity (Vs) images from the testing dataset (left column) and the
predictions made by the Transpose Convolution convolutional neural network
(CNN, middle column) and the Dense CNN (right column). Each row has been
labeled (a-e) for ease of discussion. The mean absolute percent error (MAPE)
and mean structural similarity (MSSIM) index of each predicted image is listed
above the image. Note that lower values of MAPE correspond to better Vs value
predictions, while higher values of MSSIM correspond to better layer boundary
predictions.

In comparison, we find that the performance of the Dense CNN is superior to
the Transpose Convolution CNN, as it has an overall lower MAPE of 8.5% and a
higher MSSIM index of 0.65. This network likewise predicts Vs values of the
subsurface relatively well but performs poorly when predicting some of the
more complex layer boundaries. For example, compare the true image in row b of
Fig. 6 to the Dense CNN’s prediction. The highly undulating layer boundary in
the true image is not captured by the network’s prediction and is instead
predicted as a horizontal line. This smoothing of complex layer boundaries is
somewhat expected because convolutional layers inherently blur edge and
boundary features, thus resulting in smoother outputs. Nonetheless, the Dense
network performs reasonably well in predicting subsurface Vs values and less
complex layer boundaries and should be selected over the Transpose Convolution
network due to its higher predictive accuracy. However, this network selection
ignores the internal workings of the CNNs, meaning the selection of the CNN is
made without understanding the explainability of the network.

Recall that in the input dispersion image, the high-power modal dispersion
trend provides a significant amount of information about the subsurface.
Therefore, we hypothesize that the network must learn the most from this high-
power modal trend (i.e., regions with the highest normalized power). To assess
this hypothesis, we evaluate the descriptive accuracy of the networks using
the Score-CAM algorithm to generate heatmaps that show where the CNN is
“looking” when it learns from the training dataset. These heatmaps highlight
regions of the input dispersion image that are most influential to the
predictions. To illustrate the results of Score-CAM for both networks, an
input dispersion image was randomly selected as an example (Fig. 8a). Figs. 8b
and c show the heatmaps generated using Score-CAM for this image for the
Transpose Convolution and Dense networks, respectively, and these heatmaps are
superimposed onto a translucent grayscale version of the input dispersion
image in Figs. 8d and e. Regions of the dispersion image that should
theoretically contain less important information (i.e., low-power non-modal
regions, highlighted as i and ii) are outlined in red in Fig. 8a and these
regions are shown again in the heatmaps in Figs. 8b and c.

The warmer colors in Score-CAM heatmaps indicate regions of positive influence
on the CNN’s decisions, while cooler colors indicate regions of no influence.
In Figs. 8b and c, the portion of the dispersion image containing high
normalized power is clearly visible as a warmer region, meaning both networks
strongly consider this part of the input image when learning to make
predictions. Additionally, a large region of the Transpose Convolution
network’s heatmap (Fig. 8b-i) is highlighted as a warm region, which indicates
that this region affects its predictions. This explanation of the prediction
is non-physical, as the highlighted region corresponds to dispersion data not
associated with the high-power modal trend and not used during typical surface
wave inversions (Fig. 8a-i). Conversely, very little of this region is
highlighted in the Dense network’s heatmap (Fig. 8c-i), indicating that it
does not consider this portion of the input image when making predictions.
Additionally, although the bottom region of the original dispersion image
(Fig. 8a-ii) contains no significant information, the Transpose Convolution
network strongly considers the bottom right corner of the image to be
important (Fig. 8b-ii), while the Dense network considers regions in the
bottom left and right corners as significant (Fig. 8c-ii). These results are
typical for many other input images reviewed throughout this study, although
we present only one image in Fig. 8 for brevity.

Figure 7: (a) A randomly selected dispersion input image from the training
dataset, the average heatmaps for all layers in the (b) Transpose Convolution
network and the (c) Dense network generated using Score-CAM with respect to
the selected dispersion input image, and (d, e) the respective heatmaps
superimposed on the original input image (shown in translucent grayscale for
clarity). Note that warmer colors indicate regions that positively influence
the CNN’s decisions, while cooler colors indicate regions that do not
influence the CNN’s decisions. Regions containing less important information
in the dispersion image (a-i, a-ii) are outlined in red and are shown again in
the associated heatmaps.

Figure 8: The SHAP values calculated using Deep SHAP for the selected input
image shown as Fig. 8a for the (a) Transpose Convolution network and (b) Dense
network. Red colors indicate regions that positively influence the CNN’s
decisions, blue colors indicate regions that do not influence the CNN’s
decisions, and the magnitudes of the SHAP values indicate the significance of
the influence.

We now use another method, Deep SHAP, to provide another perspective on the
CNN’s logic. Using the same input dispersion image as Fig. 8a, the SHAP values
measuring the importance of regional differences relative to the first layer
in the network are shown in Figs. 8a and b for the Transpose Convolution and
Dense networks, respectively. Note that in Fig. 8, positive SHAP values are
red in color, while negative SHAP values are blue in color. The SHAP values
reveal that both networks strongly consider portions of the high-power modal
trend at low frequencies relative to the rest of the input when making
predictions. However, the mix of negative and positive SHAP values indicates
the Transpose Convolution network has difficulties identifying the important
features in this region. Fig. 8i further shows that the Transpose Convolution
network incorrectly considers unimportant parts of the input image (recall
Fig. 8a-i) as important to its decision making. In contrast, the Dense network
considers more of the high-power modal trend as a positive influence on its
decision making, despite indicating some of the lower left portions of the
dispersion image as erroneously important. Therefore, we find that the Dense
network follows our hypothesis more closely than the Transpose Convolution
network.

When reviewing the predictive accuracy of the Transpose Convolution and Dense
networks, it seems that the Dense network is superior, as it has a higher
overall MAPE and MSSIM index. Additionally, the descriptive accuracy of the
Dense network is relatively higher than the Transpose Convolution network,
although its logic could likely be improved through hyperparameter tuning.
Therefore, we choose a Dense layer as the decoder structure for the CNN given
its relatively higher predictive and descriptive accuracy, and we perform
additional hyperparameter tuning to further improve the network’s descriptive
accuracy. We next focus on tuning the kernel size of the convolutional layers
used in the network’s encoder structure.

###  4.2 Selecting a convolution kernel size in the encoder

We next evaluate two different kernel sizes for the encoder in the Dense
network: (1) a traditional square 3$\times$3 kernel for all convolutional
layers and (2) an atypical rectangular 3$\times$1 kernel for the first three
convolutional layers. We compare the performance of the two kernel sizes in
preserving the frequency-dependent information (i.e., velocity dispersion
occurring at each frequency) before convolution with other frequencies. To
compare kernel sizes, we use the previous Dense network as the baseline (now
referred to as the 3$\times$3 network) and compare this to a similar network
in which only the kernel size of the first three convolutional layers is
changed to 3$\times$1 (referred to as the 3$\times$1 network). Both networks
are trained and tested using the same datasets described previously. We use
the same learning rate, batch size, number of training epochs, optimizer, and
loss function as before. Table 2 provides a detailed description of the
3$\times$3 and 3$\times$1 CNNs, and the kernel sizes of the first three
convolutional layers are highlighted to emphasize the difference between the
networks. Figs. 9a and b show the MAE loss per epoch during training for the
3$\times$3 network and the 3$\times$1 network, respectively. Both networks
performed comparably in terms of MAE loss reduction during validation, and we
emphasize that neither network shows signs of overfit.

Table 2: A detailed description of the convolutional neural networks (CNNs) designed with 3$\times$3 kernels (left) and 3$\times$1 kernels (right). Note that the only difference between the networks is the kernel size of the first three convolutional layers, which is highlighted. |  3$\times$3 Network |  3$\times$1 Network
---|---|---
Network Layer Type |  Filters |  Kernel/Pool Size |  Output Shape |  Kernel/Pool Size |  Output Shape
Input |  |  |  400$\times$76$\times$1 (Dispersion Image) |  |  400$\times$76$\times$1 (Dispersion Image)
2D Convolution |  32 |  (3, 3) |  398$\times$74$\times$32 |  (3, 1) |  398$\times$76$\times$32
2D Max Pooling |  |  (3, 1) |  132$\times$74$\times$32 |  (3, 1) |  132$\times$76$\times$32
2D Convolution |  32 |  (3, 3) |  130$\times$72$\times$32 |  (3, 1) |  130$\times$76$\times$32
2D Max Pooling |  |  (3, 1) |  43$\times$72$\times$32 |  (3, 1) |  43$\times$76$\times$32
2D Convolution |  64 |  (3, 3) |  41$\times$70$\times$64 |  (3, 1) |  41$\times$76$\times$64
2D Max Pooling |  |  (1, 3) |  41$\times$23$\times$64 |  (1, 3) |  41$\times$25$\times$64
2D Convolution |  128 |  (3, 3) |  39$\times$21$\times$128 |  (3, 3) |  39$\times$23$\times$128
2D Max Pooling |  |  (3, 3) |  13$\times$7$\times$128 |  (3, 3) |  13$\times$7$\times$128
2D Convolution |  128 |  (3, 3) |  11$\times$5$\times$128 |  (3, 3) |  11$\times$5$\times$128
Flatten |  |  |  7040 |  |  7040
Dense |  |  |  1152 (Flattened Vs Image) |  |  1152 (Flattened Vs Image)
Figure 9: Mean absolute error (MAE) loss reduction observed for the training
and validation datasets during convolutional neural network (CNN) optimization
for (a) the 3$\times$kernel network and (b) the 3$\times$kernel network. Note
that neither network is overfit. Figure 10: Comparisons between five randomly
selected subsurface shear-wave velocity (Vs) images from the testing dataset
(left column) and the predictions made by the 3$\times$3 convolutional neural
network (CNN, middle column) and the 3$\times$1 CNN (right column). Each row
has been labeled (a-e) for ease of discussion. The mean absolute percent error
(MAPE) and mean structural similarity (MSSIM) index of each predicted image is
listed above the image. Note that lower values of MAPE correspond to better Vs
value predictions, while higher values of MSSIM correspond to better layer
boundary predictions.

We now compare the predictions of both networks for the same randomly selected
set of images from the testing dataset shown in Fig. 6 for consistency. These
true subsurface images are shown again as the left column in Fig. 10, while
the corresponding predictions of the 3$\times$3 and 3$\times$1 networks are
shown as the middle and right columns, respectively. In Fig. 10, each
prediction made by the 3$\times$3 network (middle column) is visually similar
to each respective prediction made by the 3$\times$1 network (right column),
and the MAPE and MSSIM index values calculated for these predictions are
nearly identical. Overall, the 3$\times$3 network has an MAPE of 8.5% and
MSSIM index of 0.65, and the 3$\times$1 network has similar values of 8.8% and
0.65, respectively. These networks predict Vs values of the subsurface quite
well, but perform poorly when predicting the more complex layer boundaries.
For example, compare the true image in row b of Fig. 10 to each network’s
prediction. The highly undulating layer boundary in the true image is not
captured by the predictions of either network and is instead predicted as a
horizontal line. This is again expected because convolutional layers
inherently blur edge and boundary features as the network learns, thus
resulting in smoother outputs. Despite this blurring of complex layer
boundaries, the predictions capture the average depth of the boundaries quite
well.

Because both networks perform similarly, the challenge is now choosing between
the 3$\times$3 network and the 3$\times$1 network (i.e., which kernel size to
use in the first three convolutional layers). The simplest choice is to select
the 3$\times$3 kernel, as it provides the marginally lower MAPE and the same
MSSIM index. However, this choice would again ignore the descriptive accuracy
of each CNN, and so we propose using post hoc explainers to provide more
insight regarding the choice of our final network.

We again evaluate each network’s descriptive accuracy by developing heatmaps
that provide insight into what the network considers important in each
dispersion input image. We require that the network learns the most
information from the high-power modal trend in the input image. The results of
Score-CAM averaged across all layers in either network are shown in Fig. 11,
where the same input dispersion image shown in Fig. 8a is shown again as Fig.
11a, the heatmaps generated using Score-CAM for this image are shown as Figs.
11b and c for the 3$\times$3 and 3$\times$1 networks, respectively, and these
heatmaps are superimposed onto a translucent grayscale version of the input
dispersion image in Figs. 11d and e. We also used Deep SHAP to quantify the
relative importance of each part of the input image and found that the Deep
SHAP results confirmed the trends present in the Score-CAM heatmaps as before.
Therefore, we omit the Deep SHAP plots for brevity.

Figure 11: (a) A randomly selected dispersion input image from the training
dataset, the average heatmaps for all layers in the (b) 3$\times$3 network and
the (c) 3$\times$1 network generated using Score-CAM with respect to the
selected dispersion input image, and (d, e) the respective heatmaps
superimposed on the original input image (shown in translucent grayscale for
clarity). Note that warmer colors indicate regions that positively influence
the CNN’s decisions, while cooler colors indicate regions that do not
influence the CNN’s decisions. The region outlined in red in (a) illustrates a
portion of the input image containing less important information.

Although both networks seem to perform similarly when considering their
predictions, MAPEs, and MSSIM indices, the results of Score-CAM show that the
networks’ learned features differ. In Fig. 11b, we see that while the
3$\times$3 network’s heatmap highlights the general shape of the high-power
modal trend in the input image, it also strongly considers the entire left
portion of the image. Fig. 11c shows that the 3$\times$1 network can better
distinguish the general shape of the high-power modal trend from surrounding
features. Despite this, both networks still fixate on some unimportant
features, such as regions in the bottom left and right corners of the input
image where only low-power dispersion data are present in the input image, as
illustrated by the red outlined region in Fig. 11a. Note that while we present
a single input image in Fig. 11, this finding is typical across other input
images analyzed in this study.

Both explainability methods reveal a fixation on regions that are typically
unimportant in dispersion images used for standard surface wave testing, and
we believe both networks are excessively complex, leading to learning from
noise or insignificant features. Therefore, we now use Score-CAM to develop
averaged heatmaps for each convolutional layer in both networks to visualize
when the network begins to learn from these insignificant features. To do so,
we first calculate a heatmap for each individual layer in the network, then
simply average the results at each convolutional layer. The results are shown
in Fig. 12, where the original input image is repeated as Fig. 12a and the
heatmaps for the 3$\times$3 network and 3$\times$1 network are shown in the
left and right columns, respectively. Each convolutional layer is labeled to
the right of each row and increasing row depth corresponds to increasing
network depth. Note that these heatmaps are shown at a higher intensity and
are not superimposed on the original input image for better visualization.

Figure 12: (a) A selected dispersion input image from the training dataset and
the average Score-CAM heatmaps for each convolutional layer in the 3$\times$3
network (left column) and the 3$\times$1 network (right column). Note that
warmer colors indicate regions that positively influence the CNN’s decisions,
while cooler colors indicate regions that do not influence the CNN’s
decisions. Increasing row depth corresponds to increasing CNN depth, as shown
to the right of the figure.

In Fig. 12, we see that both networks learn a significant amount of
information regarding the high-power modal trend in the input image within the
first two convolutional layers (Figs. 12b-e). However, beginning at the third
convolutional layer, we see the networks begin to focus on noise in the bottom
right corner (outlined in red). At the fourth convolutional layer, the
3$\times$3 network begins focusing on additional noisy features at the bottom
left corner. By the final convolutional layer, both networks focus on noisy
features (outlined in red), with the 3$\times$3 network considering relatively
more noisy features as relevant to its decision making. We believe this
fixation on low-power noise indicates a type of network overfit, although
neither network shows signs of conventional overfit when considering their
predictive accuracies (recall Fig. 9). While the 3$\times$1 network has a
relatively higher descriptive accuracy compared to the 3$\times$3 network, we
find that both networks begin capturing noisy features, and thus we create
shallower architectures as a means to improve their descriptive accuracy.

###  4.3 Selecting the network depth

As the CNN learns most of the high-power modal trend features in the first two
layers, we shorten the previous two architectures to have only two 2D
convolutional layers instead of five. Additionally, our only pooling layer now
uses a 3$\times$3 pooling area instead of a 3$\times$1 area, as this
dramatically reduces the number of trainable parameters from approximately 350
million to 100 million, thus reducing the training time from approximately 26
hours to 10 hours. We retain the final Dense layer as our decoder structure
given its superior performance compared to the decoder structure of the
Transpose Convolution CNN. We also use the same learning rate, batch size,
number of training epochs, optimizer, and loss function as the previous
networks. A detailed description of both the shallow 3$\times$3 CNN and the
shallow 3$\times$1 CNN are provided in Table 3. We note that the loss
reduction during training and validation showed that neither network overfits.

Table 3: A detailed description of the shallow convolutional neural networks (CNNs) designed with 3$\times$3 kernels (left) and 3$\times$1 kernels (right). |  Shallow 3$\times$3 Network |  Shallow 3$\times$1 Network
---|---|---
Network Layer Type |  Filters |  Kernel/Pool Size |  Output Shape |  Kernel/Pool Size |  Output Shape
Input |  |  |  400$\times$76$\times$1 (Dispersion Image) |  |  400$\times$76$\times$1 (Dispersion Image)
2D Convolution |  32 |  (3, 3) |  398$\times$74$\times$32 |  (3, 1) |  398$\times$76$\times$32
2D Max Pooling |  |  (3, 3) |  132$\times$42$\times$32 |  (3, 3) |  132$\times$25$\times$32
2D Convolution |  32 |  (3, 3) |  130$\times$22$\times$32 |  (3, 1) |  130$\times$25$\times$32
Flatten |  |  |  91520 |  |  104000
Dense |  |  |  1152 (Flattened Vs Image) |  |  1152 (Flattened Vs Image)

We now compare the predictions of each revised, shallow network. The same
subset of images selected from the testing dataset shown in the left column in
Fig. 6 is presented again as the left column in Fig. 13. The predictions of
these images using the shallow 3$\times$3 network and the shallow 3$\times$1
network are shown as the middle and right columns, respectively, in Fig. 13.
When reviewing the results shown in Fig. 13, we find the shallow networks
follow similar trends to the deeper networks shown previously. Specifically,
we find that both shallow networks have: (1) relatively low MAPE values,
meaning Vs values are generally well-predicted; (2) relatively low MSSIM
values for more complex models, which correspond to poor layer boundary
predictions; and (3) smoothed layer boundaries resulting from the
convolutional layers in the network. Overall, the shallow 3$\times$3 network
has an MAPE of 11.2% and an MSSIM index of 0.58, while the shallow 3$\times$1
network has an MAPE of 12.1% and an MSSIM index of 0.57. In general, the
networks provide similar predictions to one another, although the shallow
3$\times$1 network typically has slightly higher errors. Despite their
differences, both shallow networks perform poorly relative to the deeper
networks shown in the previous section, as shown in Table 4. Recall, however,
that the choice to use a shallow network was based on the network’s fixation
on unimportant features seen in the explanations of the deeper CNNs (Fig. 12).
Therefore, we attempt to explain the predictions of the shallow networks.

Table 4: A comparison between the MAPE (shown as a percentage) and MSSIM index (shown in parentheses) between the deep and shallow 3$\times$3 and 3$\times$1 networks. |  3$\times$3 Network |  3$\times$1 Network
---|---|---
Deep |  |
(5 convolutional layers) |  8.5% (0.65) |  8.8% (0.65)
Shallow |  |
(2 convolutional layers) |  11.2% (0.58) |  12.1% (0.57)

Figure 13: Comparisons between five randomly selected subsurface shear-wave
velocity (Vs) images from the testing dataset (left column) and the
predictions made by the shallow 3$\times$3 convolutional neural network (CNN,
middle column) and the shallow 3$\times$1 CNN (right column). Each row has
been labeled (a-e) for ease of discussion. The mean absolute percent error
(MAPE) and mean structural similarity (MSSIM) index of each predicted image is
listed above the image. Note that lower values of MAPE correspond to better Vs
value predictions, while higher values of MSSIM correspond to better layer
boundary predictions.

To evaluate the descriptive accuracy of our networks, we again use Score-CAM
to explain our networks’ decisions. The results of Score-CAM averaged across
all layers in either network are shown in Fig. 14 for the same input
dispersion image in Fig. 8a. The heatmaps generated using Score-CAM for this
image are shown as Figs. 14b and c for the shallow 3$\times$3 and 3$\times$1
networks, respectively, and these heatmaps are superimposed onto a translucent
grayscale version of the input image in Figs. 14d and e.

The Score-CAM results shown in Figs. 14b and d reveal that the shallow
3$\times$3 network has difficulty distinguishing features in the input image,
while Figs. 14c and e show the shallow 3$\times$1 network considers nearly the
entire high-power modal trend as important to its decision making. Further,
unlike the deep 3$\times$1 network in Fig. 11e, the shallow 3$\times$1 network
does not fixate on other low-power features present in the input dispersion
image (Fig. 14e), which we believe indicates that this network is not
excessively complex. As with the previous hyperparameter tuning, these methods
were used to analyze a variety of randomly selected input dispersion images,
and the trends presented here are representative of the training dataset.
Although the shallow 3$\times$1 network has a relatively low predictive
accuracy compared to other dense-layer networks considered in this study
(recall Table 4), we find that it has the highest descriptive accuracy and
therefore makes the most logical predictions. Therefore, we recommend the CNN
with two convolutional layers using a 3$\times$1 kernel as our final network,
as it yields a high descriptive accuracy without a considerable loss to
predictive accuracy.

Figure 14: (a) A selected dispersion input image from the training dataset,
the average heatmaps for all layers in the (b) shallow 3$\times$3 kernel
network and the (c) shallow 3$\times$1 kernel network generated using Score-
CAM with respect to the selected dispersion input image, and (d, e) the
respective heatmaps superimposed on the original input image (shown in
translucent grayscale for clarity). Warmer colors indicate regions that
positively influence the CNN’s decisions, while cooler colors indicate regions
that do not influence the CNN’s decisions.

We recommend a shallow two-layer CNN that yields high descriptive accuracy but
a relatively lower predictive accuracy, as we emphasize creating a physics-
aware network that makes logical predictions. Although black box models such
as CNNs are typically designed to have high predictive accuracy, they may
sometimes learn incorrect physical behaviors, which may not be noticeable when
tested on the local dataset, i.e., the network is interpolating (Kumar et al.,
2022). This leads to issues when trying to generalize the model to conditions
absent from the local dataset (i.e., different subsurface conditions, type of
active source, etc.). By focusing on improving our network’s descriptive
accuracy, we create a more robust model with respect to modeling the physics
behind our problem, thus leading to more trustworthy predictions. However, we
note that shallower networks are typically sensitive to noisy inputs (recall
that deeper layers in CNNs capture more complex information in the input
image). Therefore, while we recommend beginning with the shallow 3$\times$1
network proposed herein, we further recommend that researchers looking to
develop networks for field applications should consider adding additional
denoising layers to the model. Alternatively, researchers may add noise to
their training dataset, resulting in a network that is more robust to noisy
field data while still being physics-aware.

##  5 Conclusion

As deep learning models, such as CNNs, have increased in popularity in recent
years, the demand for explainable models has grown significantly. Although
post hoc explanations have offered insights into the decision-making process
of CNNs, the neural network design largely relies on brute-force optimization
of the architecture that yields the lowest predictive error. We propose a
framework for designing deep neural networks using explainability. We
demonstrate that the descriptive accuracy of a network offers insights and
guidelines for selecting its hyperparameters. We present the XAI-guided design
of a CNN for predicting shallow subsurface 2D Vs profiles based on input
dispersion images obtained from active-source wavefields recorded by a linear
array of surface sensors. We develop initial CNNs trained on subsurface Vs and
surface wave dispersion image pairs. We define the predictive accuracy of each
network using two metrics: the Mean Absolute Percentage Error (MAPE), which
captures the overall prediction of Vs, and the Mean Structural Similarity
Index (MSSIM), which captures the interface boundary. The intial Transpose
Convolution and Dense networks show a relatively low overall MAPE and high
MSSIM value, indicating good performance. However, post hoc explainability
methods, such as Score-CAM and Deep SHAP, show that these networks provide
predictions based on some of the low-power regions in the input dispersion
images containing no meaningful information. We select the Dense network as a
baseline due to its relatively higher descriptive accuracy compared to the
Transpose Convolution network, and further explore the predictive and
descriptive accuracies of using different kernel sizes: a conventional square
3$\times$3 kernel and a non-conventional rectangular 3$\times$1 kernel. The
new 3$\times$3 and 3$\times$1 networks result in relatively similar predictive
accuracies. However, post hoc explanations reveal that these networks use
regions of low power in the input images for their predictions, while also
focusing on non-physical, unimportant, and noisy features. We minimize this
effect by altering the network depth to be relatively shallow (i.e., fewer
convolutional layers). Although the shallow 3$\times$3 and 3$\times$1 networks
perform slightly worse than the previous deeper networks in predictive
accuracy, their model explanations showed that the 3$\times$1 network made
logical decisions focused primarily on the high-power modal trends and avoided
overfit (i.e., had higher descriptive accuracy). Through this iterative
process, we demonstrate the usefulness and necessity of using model
explanations to influence the development and design of physics-aware deep
learning models.

## Acknowledgments

This material is based upon work supported by the National Science Foundation
Graduate Research Fellowship under Grant No. DGE – 1610403. However, any
opinions, findings, and conclusions or recommendations expressed in this
material are those of the authors and do not necessarily reflect the views of
the National Science Foundation.

The open-source software DENISE (Köhn, 2011; Köhn et al., 2012) was used to
perform all wave propagation simulations conducted for this study. The Texas
Advanced Computing Center’s (TACC’s) clusters Stampede2 and Frontera were used
during the construction of our dataset and the training and testing of our
CNNs, with an allocation provided by the DesignSafe-CyberInfrastructure
(Rathje et al., 2017). The open-source machine learning libraries TensorFlow
(Abadi et al., 2016) and Keras (Chollet et al., 2015) were used to develop,
train, and test the CNNs presented in this study. The wavefield
transformations used to develop the dispersion images were performed using the
open-source Python package _swprocess_ (Vantassel et al., 2022).

## Data availability

Further information regarding the networks and dataset (including subsurface
models, wavefields, and dispersion images) used in this study is available
from the corresponding author upon reasonable request.

## References

  * Haskell (1953) N. A. Haskell,  The dispersion of surface waves on multilayered media,  Bulletin of the Seismological Society of America 43 (1953) 17–34. doi:10.1785/BSSA0430010017.
  * Thomson (1950) W. T. Thomson,  Transmission of elastic waves through a stratified solid medium,  Journal of Applied Physics 21 (1950) 89–93. doi:10.1063/1.1699629.
  * Cox and Teague (2016) B. R. Cox, D. P. Teague,  Layering ratios: a systematic approach to the inversion of surface wave data in the absence of a priori information,  Geophysical Journal International 207 (2016) 422–438. doi:10.1093/gji/ggw282.
  * Foti et al. (2014) S. Foti, C. G. Lai, G. J. Rix, C. Strobbia, Surface wave methods for near-surface site characterization, CRC press, 2014\. doi:10.1201/b17268.
  * Foti et al. (2018) S. Foti, F. Hollender, F. Garofalo, D. Albarello, M. Asten, P.-Y. Bard, C. Comina, C. Cornou, B. Cox, G. Di Giulio, et al.,  Guidelines for the good practice of surface wave analysis: a product of the interpacific project,  Bulletin of Earthquake Engineering 16 (2018) 2367–2420. doi:10.1007/s10518-017-0206-7.
  * Vantassel and Cox (2022) J. P. Vantassel, B. R. Cox,  Swprocess: a workflow for developing robust estimates of surface wave dispersion uncertainty,  Journal of Seismology (2022) 1–26. doi:10.1007/s10950-021-10035-y.
  * Dramsch (2020) J. S. Dramsch,  70 years of machine learning in geoscience in review,  Advances in Geophysics 61 (2020) 1–55. doi:10.48550/arXiv.2006.13311.
  * Li et al. (2019) S. Li, B. Liu, Y. Ren, Y. Chen, S. Yang, Y. Wang, P. Jiang,  Deep-learning inversion of seismic data,  arXiv preprint arXiv:1901.07733 (2019). doi:10.1109/TGRS.2019.2953473.
  * Sun and Demanet (2020) H. Sun, L. Demanet,  Extrapolated full-waveform inversion with deep learningefwi-cnn,  Geophysics 85 (2020) R275–R288. doi:10.1190/geo2019-0195.1.
  * Wang et al. (2018) W. Wang, F. Yang, J. Ma,  Velocity model building with a modified fully convolutional network,  in: 2018 SEG International Exposition and Annual Meeting, OnePetro, 2018. doi:10.1190/segam2018-2997566.1.
  * Wu et al. (2018) Y. Wu, Y. Lin, Z. Zhou,  Inversionnet: Accurate and efficient seismic waveform inversion with convolutional neural networks,  in: 2018 SEG International Exposition and Annual Meeting, OnePetro, 2018. doi:10.1190/segam2018-2998603.1.
  * Zheng et al. (2019) Y. Zheng, Q. Zhang, A. Yusifov, Y. Shi,  Applications of supervised deep learning for seismic interpretation and inversion,  The Leading Edge 38 (2019) 526–533. doi:10.1190/tle38070526.1.
  * Vantassel et al. (2022) J. P. Vantassel, K. Kumar, B. R. Cox,  Using convolutional neural networks to develop starting models for near-surface 2-d full waveform inversion,  Geophysical Journal International 231 (2022) 72–90. doi:10.1093/gji/ggac179.
  * Abbas et al. (2022) A. Abbas, J. P. Vantassel, B. R. Cox, K. Kumar, J. Crocker,  A frequency-velocity cnn for developing near-surface 2d vs images from linear-array, active-source wavefield measurements,  arXiv preprint arXiv:2207.09580 (2022).
  * Murdoch et al. (2019) W. J. Murdoch, C. Singh, K. Kumbier, R. Abbasi-Asl, B. Yu,  Definitions, methods, and applications in interpretable machine learning,  Proceedings of the National Academy of Sciences 116 (2019) 22071–22080. doi:10.1073/pnas.1900654116.
  * Raissi et al. (2019) M. Raissi, P. Perdikaris, G. Karniadakis,  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations,  Journal of Computational Physics 378 (2019) 686–707. URL: https://www.sciencedirect.com/science/article/pii/S0021999118307125. doi:10.1016/j.jcp.2018.10.045.
  * Ribeiro et al. (2016) M. T. Ribeiro, S. Singh, C. Guestrin,  ”why should i trust you?” explaining the predictions of any classifier,  in: Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, 2016, pp. 1135–1144. doi:10.18653/v1/N16-3020.
  * Zhou et al. (2016) B. Zhou, A. Khosla, A. Lapedriza, A. Oliva, A. Torralba,  Learning deep features for discriminative localization,  in: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2016, pp. 2921–2929. doi:10.1109/CVPR.2016.319.
  * Kapishnikov et al. (2019) A. Kapishnikov, T. Bolukbasi, F. Viégas, M. Terry,  Xrai: Better attributions through regions,  in: Proceedings of the IEEE/CVF International Conference on Computer Vision, 2019, pp. 4948–4957. doi:10.48550/arXiv.1906.02825.
  * Sundararajan et al. (2017) M. Sundararajan, A. Taly, Q. Yan,  Axiomatic attribution for deep networks,  in: International conference on machine learning, PMLR, 2017, pp. 3319–3328. doi:10.48550/arXiv.1703.01365.
  * Bien et al. (2018) N. Bien, P. Rajpurkar, R. L. Ball, J. Irvin, A. Park, E. Jones, M. Bereket, B. N. Patel, K. W. Yeom, K. Shpanskaya, et al.,  Deep-learning-assisted diagnosis for knee magnetic resonance imaging: development and retrospective validation of mrnet,  PLoS Medicine 15 (2018) e1002699. doi:10.1371/journal.pmed.1002699.
  * Lauritsen et al. (2020) S. M. Lauritsen, M. Kristensen, M. V. Olsen, M. S. Larsen, K. M. Lauritsen, M. J. Jørgensen, J. Lange, B. Thiesson,  Explainable artificial intelligence model to predict acute critical illness from electronic health records,  Nature Communications 11 (2020) 1–11. doi:10.48550/arXiv.1912.01266.
  * Lou et al. (2019) B. Lou, S. Doken, T. Zhuang, D. Wingerter, M. Gidwani, N. Mistry, L. Ladic, A. Kamen, M. E. Abazeed,  An image-based deep learning framework for individualising radiotherapy dose: a retrospective analysis of outcome prediction,  The Lancet Digital Health 1 (2019) e136–e147. doi:10.1016/S2589-7500(19)30058-5.
  * Rajpurkar et al. (2017) P. Rajpurkar, J. Irvin, K. Zhu, B. Yang, H. Mehta, T. Duan, D. Ding, A. Bagul, C. Langlotz, K. Shpanskaya, et al.,  Chexnet: Radiologist-level pneumonia detection on chest x-rays with deep learning,  arXiv preprint arXiv:1711.05225 (2017).
  * Rajpurkar et al. (2018) P. Rajpurkar, J. Irvin, R. L. Ball, K. Zhu, B. Yang, H. Mehta, T. Duan, D. Ding, A. Bagul, C. P. Langlotz, et al.,  Deep learning for chest radiograph diagnosis: A retrospective comparison of the chexnext algorithm to practicing radiologists,  PLoS Medicine 15 (2018) e1002686. doi:10.1371/journal.pmed.1002686.
  * Mitani et al. (2020) A. Mitani, A. Huang, S. Venugopalan, G. S. Corrado, L. Peng, D. R. Webster, N. Hammel, Y. Liu, A. V. Varadarajan,  Detection of anaemia from retinal fundus images via deep learning,  Nature Biomedical Engineering 4 (2020) 18–27. doi:10.1038/s41551-019-0487-z.
  * Zech et al. (2018) J. R. Zech, M. A. Badgeley, M. Liu, A. B. Costa, J. J. Titano, E. K. Oermann,  Variable generalization performance of a deep learning model to detect pneumonia in chest radiographs: a cross-sectional study,  PLoS Medicine 15 (2018) e1002683. doi:10.1371/journal.pmed.1002683.
  * Linardatos et al. (2020) P. Linardatos, V. Papastefanopoulos, S. Kotsiantis,  Explainable ai: A review of machine learning interpretability methods,  Entropy 23 (2020) 18. doi:10.3390/e23010018.
  * Hinz et al. (2018) T. Hinz, N. Navarro-Guerrero, S. Magg, S. Wermter,  Speeding up the hyperparameter optimization of deep convolutional neural networks,  International Journal of Computational Intelligence and Applications 17 (2018) 1850008\. doi:10.1142/S1469026818500086.
  * Lee et al. (2018) W.-Y. Lee, S.-M. Park, K.-B. Sim,  Optimal hyperparameter tuning of convolutional neural networks based on the parameter-setting-free harmony search algorithm,  Optik 172 (2018) 359–367. doi:10.1016/j.ijleo.2018.07.044.
  * Loshchilov and Hutter (2016) I. Loshchilov, F. Hutter,  Cma-es for hyperparameter optimization of deep neural networks,  arXiv preprint arXiv:1604.07269 (2016).
  * Snoek et al. (2012) J. Snoek, H. Larochelle, R. P. Adams,  Practical bayesian optimization of machine learning algorithms,  Advances in Neural Information Processing Systems 25 (2012). doi:10.48550/arXiv.1206.2944.
  * Wu et al. (2019) J. Wu, X.-Y. Chen, H. Zhang, L.-D. Xiong, H. Lei, S.-H. Deng,  Hyperparameter optimization for machine learning models based on bayesian optimization,  Journal of Electronic Science and Technology 17 (2019) 26–40. doi:10.11989/JEST.1674-862X.80904120.
  * Araya-Polo et al. (2018) M. Araya-Polo, J. Jennings, A. Adler, T. Dahlke,  Deep-learning tomography,  The Leading Edge 37 (2018) 58–66. doi:10.1190/tle37010058.1.
  * Çaylak and Kaftan (2014) Ç. Çaylak, İ. Kaftan,  Determination of near-surface structures from multi-channel surface wave data using multi-layer perceptron neural network (mlpnn) algorithm,  Acta Geophysica 62 (2014) 1310–1327. doi:10.2478/s11600-014-0207-8.
  * Hu et al. (2021) W. Hu, Y. Jin, X. Wu, J. Chen,  Progressive transfer learning for low-frequency data prediction in full-waveform inversion,  Geophysics 86 (2021) R369–R382. doi:10.1190/geo2020-0598.1.
  * Lu et al. (2019) P. Lu, Y. Zhang, J. Chen, Y. Xiao, G. Zhao,  Enhanced seismic imaging with predictive neural networks for geophysics,  arXiv preprint arXiv:1908.03973 (2019).
  * Yang and Ma (2019) F. Yang, J. Ma,  Deep-learning inversion: A next-generation seismic velocity model building method,  Geophysics 84 (2019) R583–R599. doi:10.1190/geo2018-0249.1.
  * Zhu et al. (2022) W. Zhu, K. Xu, E. Darve, B. Biondi, G. C. Beroza,  Integrating deep neural networks with full-waveform inversion: Reparameterization, regularization, and uncertainty quantification,  Geophysics 87 (2022) R93–R109. doi:10.1190/geo2020-0933.1.
  * Goodfellow et al. (2016) I. Goodfellow, Y. Bengio, A. Courville, Deep learning, MIT press, 2016.
  * Wang et al. (2004) Z. Wang, A. C. Bovik, H. R. Sheikh, E. P. Simoncelli,  Image quality assessment: from error visibility to structural similarity,  IEEE transactions on image processing 13 (2004) 600–612. doi:10.1109/TIP.2003.819861.
  * Wang et al. (2020) H. Wang, Z. Wang, M. Du, F. Yang, Z. Zhang, S. Ding, P. Mardziel, X. Hu,  Score-cam: Score-weighted visual explanations for convolutional neural networks,  in: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops, 2020, pp. 24–25. doi:10.1109/CVPRW50498.2020.00020.
  * Lundberg and Lee (2017) S. M. Lundberg, S.-I. Lee,  A unified approach to interpreting model predictions,  Advances in Neural Information Processing Systems 30 (2017). doi:arXiv:1705.07874.
  * Selvaraju et al. (2017) R. R. Selvaraju, M. Cogswell, A. Das, R. Vedantam, D. Parikh, D. Batra,  Grad-cam: Visual explanations from deep networks via gradient-based localization,  in: Proceedings of the IEEE International Conference on Computer Vision (ICCV), 2017, pp. 618–626.
  * Chattopadhay et al. (2018) A. Chattopadhay, A. Sarkar, P. Howlader, V. N. Balasubramanian,  Grad-cam++: Generalized gradient-based visual explanations for deep convolutional networks,  in: 2018 IEEE Winter Conference on Applications of Computer Vision (WACV), IEEE, 2018, pp. 839–847. doi:10.1109/ICCV.2017.74.
  * Shrikumar et al. (2017) A. Shrikumar, P. Greenside, A. Kundaje,  Learning important features through propagating activation differences,  in: International conference on machine learning, PMLR, 2017, pp. 3145–3153.
  * Köhn (2011) D. Köhn, Time domain 2D elastic full waveform tomography, Ph.D. thesis, 2011.
  * Köhn et al. (2012) D. Köhn, D. De Nil, A. Kurzmann, A. Przebindowska, T. Bohlen,  On the influence of model parametrization in elastic full waveform tomography,  Geophysical Journal International 191 (2012) 325–345. doi:10.1111/j.1365-246X.2012.05633.x.
  * Zywicki (1999) D. J. Zywicki, Advanced signal processing methods applied to engineering analysis of seismic surface waves, Georgia Institute of Technology, 1999.
  * Vantassel (2021) J. Vantassel,  jpvantassel/swprocess:v0.1.0b0,  Zenodo 4584129 (2021). doi:10.5281/zenodo.4584129.
  * Genzel et al. (2022) M. Genzel, J. Macdonald, M. Marz,  Solving inverse problems with deep neural networks-robustness included,  IEEE Transactions on Pattern Analysis and Machine Intelligence (2022). doi:10.1109/TPAMI.2022.3148324.
  * Kingma and Ba (2014) D. P. Kingma, J. Ba,  Adam: A method for stochastic optimization,  arXiv preprint arXiv:1412.6980 (2014).
  * Kumar et al. (2022) K. Kumar, E. Rathje, J. Crocker, Can we trust ai? - a case for explainable ai (x-ai), 2022. URL: https://doi.org/10.5281/zenodo.7308507. doi:10.5281/zenodo.7308507.
  * Rathje et al. (2017) E. M. Rathje, C. Dawson, J. E. Padgett, J.-P. Pinelli, D. Stanzione, A. Adair, P. Arduino, S. J. Brandenberg, T. Cockerill, C. Dey, et al.,  Designsafe: New cyberinfrastructure for natural hazards engineering,  Natural Hazards Review 18 (2017) 06017001. doi:10.1061/(ASCE)NH.1527-6996.0000246.
  * Abadi et al. (2016) M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, et al.,  Tensorflow: Large-scale machine learning on heterogeneous distributed systems,  arXiv preprint arXiv:1603.04467 (2016).
  * Chollet et al. (2015) F. Chollet, et al., Keras, https://keras.io, 2015.