Bayesian seismic tomography based on velocity-space Stein variational
  gradient descent for physics-informed neural network
In this study, we propose a Bayesian seismic tomography inference method
using physics-informed neural networks (PINN). PINN represents a recent advance
in deep learning, offering the possibility to enhance physics-based simulations
and inverse analyses. PINN-based deterministic seismic tomography uses two
separate neural networks (NNs) to predict seismic velocity and travel time.
Naive Bayesian NN (BNN) approaches are unable to handle the high-dimensional
spaces spanned by the weight parameters of these two NNs. Hence, we reformulate
the problem to perform the Bayesian estimation exclusively on the NN predicting
seismic velocity, while the NN predicting travel time is used only for
deterministic travel time calculations, with the help of the adjoint method.
Furthermore, we perform BNN by introducing a function-space Stein variational
gradient descent (SVGD), which performs particle-based variational inference in
the space of the function predicted by the NN (i.e., seismic velocity), instead
of in the traditional weight space. The result is a velocity-space SVGD for the
PINN-based seismic tomography model (vSVGD-PINN-ST) that decreases the
complexity of the problem thus enabling a more accurate and physically
consistent Bayesian estimation, as confirmed by synthetic tests in one- and
two-dimensional tomographic problem settings. The method allows PINN to be
applied to Bayesian seismic tomography practically for the first time. Not only
that, it can be a powerful tool not only for geophysical but also for general
PINN-based Bayesian estimation problems associated with compatible NNs
formulations and similar, or reduced, complexity.
Ryoichiro Agata, Kazuya Shiraishi, Gou Fujie
Bayesian seismic tomography based on velocity-space Stein variation gradient descent for physics-informed neural network

###### Abstract

In this study, we propose a Bayesian seismic tomography inference method using physics-informed neural networks (PINN). PINN represents a recent advance in deep learning, offering the possibility to enhance physics-based simulations and inverse analyses. PINN-based deterministic seismic tomography uses two separate neural networks (NNs) to predict seismic velocity and travel time. Naive Bayesian NN (BNN) approaches are unable to handle the high-dimensional spaces spanned by the weight parameters of these two NNs. Hence, we reformulate the problem to perform the Bayesian estimation exclusively on the NN predicting seismic velocity, while the NN predicting travel time is used only for deterministic travel time calculations, with the help of the adjoint method. Furthermore, we perform BNN by introducing a function-space Stein variational gradient descent (SVGD), which performs particle-based variational inference in the space of the function predicted by the NN (i.e., seismic velocity), instead of in the traditional weight space. The result is a velocity-space SVGD for the PINN-based seismic tomography model (vSVGD-PINN-ST) that decreases the complexity of the problem thus enabling a more accurate and physically consistent Bayesian estimation, as confirmed by synthetic tests in one- and two-dimensional tomographic problem settings. The method allows PINN to be applied to Bayesian seismic tomography practically for the first time. Not only that, it can be a powerful tool not only for geophysical but also for general PINN-based Bayesian estimation problems associated with compatible NNs formulations and similar, or reduced, complexity.

 Seismic tomography, physics-informed neural network (PINN), Bayesian neural network, function-space stein variation gradient descent (fSVGD)

## I Introduction

Seismic tomography is a technique used to determine the interior seismic structure of the Earth, using seismic waves excited by natural earthquakes and artificial sources. This technique is essential for studying the evolution of the Earth, plate tectonics, and mechanisms of earthquake generation. Seismic tomography, similarly to many other physics-based estimation problems in geoscience, represents an under-determined inverse problem because the number of sources and receivers is limited. A priori constraints on the parameters of interest (e.g., seismic velocity) can be used to regularize the problem and provide the uncertainty quantification (UQ) of the estimates, ensuring reliable analytical results. To meet such requirement, recently, a growing body of research uses Bayesian estimation to perform UQ estimation, considering unknown parameters as stochastic variables. In the field of seismic tomography, there are many examples of Bayesian estimations in two-dimensional (2D) surface wave [1, 2, 3, 4], 2D seismic refraction [5], and three-dimensional (3D) tomography [6, 7, 8]. These methods employ a Bayesian estimation of the posterior probability distribution of the seismic velocity structure, combining travel time data, prior information on the target velocity structure, and forward numerical calculations of travel time based on grid or mesh discretization of the target domain.

Recent developments in deep learning techniques provide entirely new options for solving partial differential equations (PDEs) and inversion problems. Physics-informed neural networks (PINN) [9], in particular, have attracted considerable attention owing to their high applicability and flexibility. Raissi et al. [9] presented a PINN application to Navier-Stokes equations. It now has been applied to a variety of fields, including solid earth science, in forward simulations (e.g., seismic travel time calculation [10, 11], seismic wave field simulation [12], and crustal deformation modeling [13]), and inverse problems (e.g., seismic tomography [14, 15] and full-waveform inversion [16]). PINN uses neural networks (NN) to construct functions predicting physical quantities, such as travel time, taking spatial coordinates and time (for dynamic problems) as inputs. The NN is trained based on a loss function consisting of residuals of the governing PDE in some given evaluation (or collocation) points. The advantages brought by the PINN travel time calculation [10, 11] and seismic tomography in general [14, 15], such as mesh-free frameworks not requiring initial model setup, demonstrate a considerable potential for further important developments. The one considered in this study aims at extending PINN-based seismic tomography to include UQ based on Bayes' theorem.

Bayesian inversion based on PINN is an application of Bayesian neural network (BNN) where the posterior probability density function (PDF) of the weight parameters in NN is estimated. Stochastic properties of target physical quantities, such as seismic wave velocity, are obtained from the posterior predictive PDF, calculated based on the posterior PDF. Bayesian sampling methods and variational inference (VI) are the methods most commonly used to estimate the posterior PDF of the weight parameters in BNN [17]. Bayesian sampling is the most accurate method for BNN, enabling sampling from the true posterior PDF. For instance, the Hamiltonian/hybrid Monte Carlo (HMC) method [18] has been used in previous applications of Bayesian PINN [19],
[20]. However, being a sequential sampling method, HMC computation time and scalability are not suitable for large problems. In VI, the posterior PDF estimation problem is replaced with an optimization one, minimizing the Kullback-Leibler (KL) divergence between the target PDF and its parametric approximation. VI are usually based on the mean-field approximation (MFVI), which approximates the target PDF combining multiple independent simple PDFs. However, MFVI displays poor performance in Bayesian PINN, giving an over-simplified PDF approximation [19, 21]. A more recent non-parametric class of VI, called particle-based VI (ParVI) and best known as Stein variational gradient descent (SVGD) [22], has attracted attention for the characteristic high approximation accuracy and computational parallelism. ParVI methods iteratively update a set of particles and use the corresponding empirical probability measure to approximate the target posterior PDF accurately. These methods have already been applied both to BNN and Bayesian PINN problems, using relatively small NN structures (e.g., [23]). Even so, BNNs have a challenging aspect, that is, the multi-modality of the posterior PDF in the high-dimensional weight space resulting from over-parametrization. In fact, optimized networks are able to express the same function values using multiple parameter combinations. Exploring PDFs solutions using ParVI and SVGD in such an over-parametrized space is a difficult task, resulting in degraded approximation performance [24].

However, in physics problems, such as those of seismic tomography, the posterior PDFs of interest are defined in the space of the functions predicted by NNs (e.g., seismic velocities); therefore, they are expected to have a considerably simpler shape than those defined in the weight space. Performing ParVI in function space, instead of weight space with multi-modality, in BNN avoids posterior PDF degraded approximation issues [24]. In addition, considering Bayesian estimations in the function space also enables the efficient incorporation of physically meaningful components of prior information, which would not be possible in the weight space [26].

Another important aspect to consider is the confinement of a NN targeted in Bayesian estimation. In fact, previous PINN-based seismic tomography methods [14, 15] included two NNs: one predicting travel time as a solution of the governing PDE and the other seismic wave velocity. A similar NN formulation was employed in other inverse problems (e.g. [16]). Targeting the weight parameters of two NNs in Bayesian estimation is not efficient, considering we are exclusively interested in the seismic velocity UQ. Moreover, the only requirement for the travel time NN is to satisfy the governing PDE for the given velocity structure. Therefore, if we could confine the target parameters to the weights of the NN predicting seismic velocity (or other relevant functions for UQ), the exploration of the posterior PDF would be drastically simplified.

This study proposes a novel method, based on Bayesian inference and PINN, to address UQ and inverse problems applied to seismic tomography, overcoming existing limitations. The main contributions of this study can be summarized as follows. We reformulated the problem to perform the Bayesian estimation only on the velocity NN, while using the travel time NN for deterministic travel time calculations, with the help of a numerical technique called the adjoint method [27]. We performed Bayesian estimation for the velocity network by introducing SVGD directly in the space of functions predicted by the NN (i.e., seismic velocity), instead of in the whole weight space, leveraging a recent development in BNN techniques [24]. The proposed approach is called "velocity-space SVGD for PINN-based seismic tomography (vSVGD-PINN-ST)" and represents the first practical application of Bayesian estimation in PINN-based seismic tomography. Furthermore, vSVGD-PINN-ST can be generalized to the UQ of PINN-based inverse problems sharing similar NN formulation and problem size.

This paper is organized as follows: In Section II, we present a PINN formulation of seismic tomography from previous studies and explain a naive formulation of Bayesian estimation as a baseline method. In Section III, we present our vSVGD-PINN-ST, introducing the improvements presented earlier. In Section IV and V, we validate the applicability of the proposed method to realistic problems through synthetic seismic tomography tests in one and two dimensions. In Section VI and VII we discuss results and give concluding remarks, respectively.

## II Baseline method

### _PINN-based deterministic seismic tomography_

First, we first introduce the PINN formulation of the deterministic seismic tomographic problem as proposed by [14]. In this work, we follow a slightly modified approach. We start from the eikonal equation relating the spatial derivative of the travel time field to the velocity structure as follows:

\[|\nabla T(\mathbf{x},\mathbf{x}_{s})|^{2} = \frac{1}{v^{2}(\mathbf{x})},\quad\forall\,\mathbf{x}\in\Omega \tag{1}\] \[T(\mathbf{x}_{s},\mathbf{x}_{s}) = 0, \tag{2}\]

where \(\Omega\) is a \(\mathbb{R}^{d}\) domain, with \(d\) as the space dimension, \(T(\mathbf{x},\mathbf{x}_{s})\) is the travel time at the point \(\mathbf{x}\) from the source \(\mathbf{x}_{s}\), \(v(\mathbf{x})\) is the velocity defined on \(\Omega\), and \(\nabla\) denotes the gradient operator. The second equation defines the point source condition. To avoid singularities in this condition, previous studies modeling the travel time using PINN introduced the following factored form [10, 11]:

\[T(\mathbf{x},\mathbf{x}_{s})=T_{0}(\mathbf{x},\mathbf{x}_{s}),\tau(\mathbf{x},\mathbf{x}_{s}) \tag{3}\]

where \(T_{0}(\mathbf{x})\) is defined as

\[T_{0}(\mathbf{x},\mathbf{x}_{s})=\left|\mathbf{x}-\mathbf{x}_{s}\right|. \tag{4}\]

Using this factorization, the point source condition becomes:

\[\tau(\mathbf{x}_{s},\mathbf{x}_{s})=\frac{1}{v(\mathbf{x}_{s})}. \tag{5}\]

Considering \(v\) positive, we can introduce the residuals of the eikonal equation \(r_{\text{EE}}\) and the point source condition \(r_{\text{PC}}\) in terms of velocity:

\[r_{\text{EE}} = v(\mathbf{x})-\frac{1}{|\nabla T(\mathbf{x},\mathbf{x}_{s})|}, \tag{6}\] \[r_{\text{PC}} = v(\mathbf{x}_{s})-\frac{1}{\tau(\mathbf{x}_{s},\mathbf{x}_{s})}. \tag{7}\]
The PINN-based seismic tomography method solves the eikonal equation and estimates the velocity structure simultaneously, by training NNs to predict the travel-time and velocity function. Waheed et al. [14] proposed to use two different NNs to construct the functions \(f_{T}\) and \(f_{v}\), characterized by weight parameters \(\mathbf{\theta}_{T}\) and \(\mathbf{\theta}_{v}\). Optimized NNs are expected to return accurate approximations of the travel time and velocity. Such formulation, which introduces two NNs, one for the solution of the governing PDE and the other for PDE parameters optimization, is generally used in PINN-based inverse analysis [16]. In this study, we define NNs functions as

\[T(\mathbf{x},\mathbf{x}_{s}) \simeq f_{T}(\mathbf{x},\mathbf{x}_{s},\mathbf{\theta}_{T}) \tag{8}\] \[= T_{0}(\mathbf{x},\mathbf{x}_{s})/f_{\tau^{-1}}(\mathbf{x}, \mathbf{x}_{s},\mathbf{\theta}_{T}),\] \[v(\mathbf{x}) \simeq f_{v}(\mathbf{x},\mathbf{\theta}_{v})\] (9) \[= v_{0}(\mathbf{x})+f_{v_{\rm{pth}}}(\mathbf{x},\mathbf{\theta}_{v}),\]

where \(f_{\tau^{-1}}\) is an NN-based function approximating \(1/\tau(\mathbf{x},\mathbf{x}_{s})\), \(v_{0}(\mathbf{x})\) is the reference velocity set by the user, and \(f_{v_{\rm{pth}}}(\mathbf{x},\mathbf{\theta}_{v})\) is the neural network function approximating the velocity perturbation component. We approximated \(1/\tau\) and \(v_{\rm{pth}}\) using NNs, instead of directly computing \(\tau\) and \(v\) as done in previous studies [10, 11], to improve convergence performances. We used fully connected feed-forward networks to implement both \(f_{\tau^{-1}}\) and \(f_{v_{\rm{pth}}}\). The \(f_{v_{\rm{pth}}}\) network is characterized by the parameters \(v_{\rm{pth}}^{\rm{max}}\) and \(v_{\rm{pth}}^{\rm{max}}\), representing the maximum and minimum velocity, set a priori. Additional operations were applied to normalize the input and output of NNs, improving the convergence performance and setting upper and lower limits of the final output values (see Appendix A). For a deterministic tomographic problem, the two NNs are trained simultaneously using the same loss function:

\[L = \alpha_{1}\sum_{i=1}^{N_{T}}\left(T_{\rm{obs}}^{(i)}-f_{T}( \mathbf{x}_{r}^{(i)},\mathbf{x}_{s}^{(i)};\mathbf{\theta}_{T})\right)^{2} \tag{10}\] \[+ \alpha_{2}\sum_{i=1}^{N_{c}}\left(f_{v}(\mathbf{x}_{c}^{(i)}; \mathbf{\theta}_{v})-\frac{1}{|\nabla f_{T}(\mathbf{x}_{c}^{(i)},\mathbf{x}_{s}^ {(i)};\mathbf{\theta}_{T})|}\right)^{2}\] \[+ \alpha_{3}\sum_{i=1}^{N_{s}}\left(f_{v}(\mathbf{x}_{s}^{(i)};\bm {\theta}_{v})-f_{\tau^{-1}}(\mathbf{x}_{s}^{(i)},\mathbf{x}_{s}^{(i)};\mathbf{ \theta}_{T})\right)^{2},\]

where \(N_{T}\), \(N_{c}\) and \(N_{s}\) are the number of travel time data, collocation points, and source points, respectively; \(T_{\rm{obs}}\) represents the travel time data; \(\mathbf{x}_{r}\) and \(\mathbf{x}_{c}\) are the coordinates of receiver and collocation points, respectively; \(\alpha_{1}\), \(\alpha_{2}\), and \(\alpha_{3}\) are loss weights, typically assigned to \(1/N_{T}\), \(1/N_{c}\) and \(1/N_{s}\), respectively, although manual or automated tuning can be used for a stable training [28]. The collocation points, which are usually selected randomly within the target domain \(\Omega\), are set as evaluation points of the PDE residuals [9]. The first term on the right-hand side is the loss due to observational constraints, consisting of the sum of the squared difference of true and predicted travel time values for each source-receiver pair. The second term is the loss due to physics-informed constraints, consisting of the sum of the squared residuals defined in Equation 6 for each pair of source and collocation point. The third term is the one modeling the point source constraints, consisting of the sum of the squared residual defined in Equation 7 for each source point.

Following the above definitions, we defined the training datasets for the observation and PDE loss functions as \(\mathbf{X}_{T}=\{(\mathbf{x}_{r}^{(1)},\mathbf{x}_{s}^{(1)},T_{\rm{obs}}^{(1)},(\mathbf{x}_{r}^{(2)},\mathbf{x}_{s}^{(2)},T_{\rm{obs}}^{(2)}),\ldots,(\mathbf{ x}_{r}^{(N_{T})},\mathbf{x}_{s}^{(N_{T})},T_{\rm{obs}}^{(N_{T})})\}\) and \(\mathbf{X}_{c}=\{(\mathbf{x}_{c}^{(1)},\mathbf{x}_{s}^{(1)}),(\mathbf{x}_{c}^ {(2)},\mathbf{x}_{s}^{(2)}),\ldots,(\mathbf{x}_{c}^{(N_{c})},\mathbf{x}_{s}^ {(N_{c})})\}\), respectively. Tomographic estimation of \(v(\mathbf{x})\) is given by \(f_{v}(\mathbf{x},\mathbf{\theta}_{v}^{*})\), where \(\mathbf{\theta}_{T}^{*},\mathbf{\theta}_{v}^{*}=\underset{\mathbf{\theta}_{T},\mathbf{\theta} _{v}}{\arg\min}L(\mathbf{\theta}_{T},\mathbf{\theta}_{v})\) represent the NNs optimized values. Fig. 1 schematically illustrates settings and constraints used to train the NNs involved in both the tomographic deterministic and Bayesian formulation (described in the following section).

### _Naive Bayesian formulation of PINN-based seismic tomography_

In this study, we considered the weight parameters as stochastic variables, formulating the estimation problem using the Bayes' theorem. In a naive formulation, we consider a Bayesian neural network (BNN) for both \(\mathbf{\theta}_{T}\) and \(\mathbf{\theta}_{v}\) given by:

\[P(\mathbf{\theta}_{T},\mathbf{\theta}_{v}|\mathbf{d}) = \frac{P(\mathbf{d}|\mathbf{\theta}_{T},\mathbf{\theta}_{v})P(\mathbf{\theta}_ {T},\mathbf{\theta}_{v})}{P(\mathbf{d})} \tag{11}\] \[\propto P(\mathbf{d}|\mathbf{\theta}_{T},\mathbf{\theta}_{v})P(\mathbf{\theta}_{T}, \mathbf{\theta}_{v}),\]

where \(\mathbf{d}\) is the data vector, including not only travel time data but also the constraints from the eikonal equation: The likelihood function includes the two following components:

\[P(\mathbf{d}|\mathbf{\theta}_{T},\mathbf{\theta}_{v})=P(\mathbf{T}_{\rm{obs}}|\mathbf{ \theta}_{T},\mathbf{\theta}_{v})P(\mathbf{r}|\mathbf{\theta}_{T},\mathbf{\theta}_{v}). \tag{12}\]

The first term on the right-hand side represents the stochastic observation error of the travel time, for which the following Gaussian distribution is assumed:

\[P(\mathbf{T}_{\rm{obs}}|\mathbf{\theta}_{T},\mathbf{\theta}_{v})\] \[=\frac{1}{Z}\exp\left(-\frac{1}{2}(\mathbf{T}_{\rm{obs}}-\mathbf{ f}_{T}(\mathbf{\theta}_{T}))^{\top}\mathbf{E}_{\rm{obs}}^{-1}(\mathbf{T}_{\rm{obs}}- \mathbf{f}_{T}(\mathbf{\theta}_{T}))\right), \tag{13}\]

where

\[\mathbf{T}_{\rm{obs}} = \left[T_{\rm{obs}}^{(1)},T_{\rm{obs}}^{(2)},\ldots,T_{\rm{obs}}^{ (N_{T})}\right], \tag{14}\] \[\mathbf{f}_{T}(\mathbf{\theta}_{T}) = \left[f_{T}^{(1)},f_{T}^{(2)},\ldots,f_{T}^{(N_{T})}\right],\] (15) \[f_{T}^{(i)} = f_{T}(\mathbf{x}_{r}^{(i)},\mathbf{x}_{s}^{(i)};\mathbf{\theta}_{T}), \tag{16}\]

and \(\mathbf{E}_{\rm{obs}}\) and \(Z\) are the covariance matrix for the observation error and a normalizing coefficient, respectively. We see that \(P(\mathbf{T}_{\rm{obs}}|\mathbf{\theta}_{T},\mathbf{\theta}_{v})\) is an explicit function of \(\mathbf{\theta}_{T}\), implicitly dependent on \(\mathbf{\theta}_{v}\). To materialize the second term of Equation 12, we assume Gaussian properties for the residuals as well:

\[P(\mathbf{r}|\mathbf{\theta}_{T},\mathbf{\theta}_{v})\] \[=\frac{1}{Z}\exp\left(-\frac{1}{2}\mathbf{r}^{\top}(\mathbf{\theta}_{ T},\mathbf{\theta}_{v})\mathbf{E}_{r}^{-1}\mathbf{r}(\mathbf{\theta}_{T},\mathbf{\theta}_{v}) \right), \tag{17}\]
where

\[\mathbf{r}(\mathbf{\theta}_{T},\mathbf{\theta}_{v})=\left[r^{(1)}_{\rm EE},r^{(2)}_{\rm EE },\ldots,r^{(N_{e})}_{\rm EE},r^{(1)}_{\rm PC},r^{(2)}_{\rm PC},\ldots,r^{(N_{s}) }_{\rm PC}\right], \tag{18}\]

\[r^{(i)}_{\rm EE}=f_{v}(\mathbf{x}^{(i)}_{c};\mathbf{\theta}_{v})-\frac{1}{|\nabla f _{T}(\mathbf{x}^{(i)}_{c},\mathbf{x}^{(i)}_{s};\mathbf{\theta}_{T})|}, \tag{19}\]

\[r^{(i)}_{\rm PC}=f_{v}(\mathbf{x}^{(i)}_{s};\mathbf{\theta}_{v})-f_{\tau^{-1}}( \mathbf{x}^{(i)}_{s},\mathbf{x}^{(i)}_{s};\mathbf{\theta}_{T}), \tag{20}\]

and \(\mathbf{E}_{\rm r}\) is the covariance matrix of PDE residuals, which may only have diagonal components given depending on the confidence level of the forward model. We use the notations \(\mathbf{r}(\mathbf{\theta}_{T},\mathbf{\theta}_{v})\) and \(\mathbf{r}\) intermittently. The simplest choices for the prior PDF \(P(\mathbf{\theta}_{T},\mathbf{\theta}_{v})\) are, for instance, an independent and identically distributed (i.i.d.) zero-mean Gaussian distribution [19] or a student's \(t\)-distribution [23]. Once the specific forms of the likelihood function and prior PDF are determined, the approximate posterior PDF can be obtained by several Bayesian estimation methods, introduced in the next section. The stochastic property of the seismic velocity is then obtained as the predictive PDF \(P(\mathbf{v}|\mathbf{d})\) based on the marginal posterior PDF of \(\mathbf{\theta}_{v}\) as follows:

\[P(\mathbf{v}|\mathbf{d}) = \int P(\mathbf{v}|\mathbf{\theta}_{v})P(\mathbf{\theta}_{v}|\mathbf{d})d \mathbf{\theta}_{v}. \tag{21}\]

### _Stein variation gradient descent for naive Bayesian PINNN-based seismic tomography_

The naive Bayesian formulation involves the Bayesian estimation of the weight parameters posterior PDF \(P(\mathbf{\theta}_{T},\mathbf{\theta}_{v}|\mathbf{d})\). ParVI methods, such as SVGD [22], are paradigms of Bayesian estimation, recently gaining more popularity in various fields, including geophysics research (e.g., [4, 29, 30]). SVGD is known for a more efficient approximation ability in the calculation of the posterior distribution compared to HMC, which has degraded performance with high dimensional problems and large datasets. As in ordinary MFVI, in SVGD the Bayesian estimation is replaced with a minimization problem of the KL divergence, defined as follows:

\[KL\left(Q(\mathbf{\theta})\|P(\mathbf{\theta}|\mathbf{d})\right)=\int Q(\mathbf{\theta}) \log\frac{Q(\mathbf{\theta})}{P(\mathbf{\theta}|\mathbf{d})}d\mathbf{\theta}, \tag{22}\]

where \(P(\mathbf{\theta}|\mathbf{d})\) and \(Q(\mathbf{\theta})\) are the target posterior and approximate distributions, respectively. SVGD employs a set of particles \(\{\mathbf{\theta}\}_{i=1}^{n}\) to approximate the target posterior PDF by minimizing the KL divergence. These particles iteratively move towards the posterior distribution following the gradient of the KL divergence \(\mathbf{\varphi}\), which is obtained from the kernelized Stein discrepancy defined in a reproducing kernel Hilbert space (RKHS). The update equations of SVGD are given as follows:

\[\mathbf{\theta}_{i}^{l+1}=\mathbf{\theta}_{i}^{l}+\epsilon_{l}\mathbf{\phi}\left(\mathbf{ \theta}_{i}^{l}\right), \tag{23}\]

where

\[\mathbf{\phi}(\mathbf{\theta})=\frac{1}{n}\sum_{j=1}^{n}\{k(\mathbf{\theta}_{j}^{l},\mathbf{ \theta})\nabla_{\mathbf{\theta}_{j}^{l}}\log P(\mathbf{\theta}_{j}^{l}|\mathbf{d})+ \nabla_{\mathbf{\theta}_{j}^{l}}k(\mathbf{\theta}_{j}^{l},\mathbf{\theta}_{v})\}, \tag{24}\]

\(\epsilon_{l}\) is the step size at each iteration \(l\), which can be determined using various adaptive optimizers (such as Adam [31]), and \(k(\mathbf{x},\cdot)\) represents a positive definite kernel. In particular, we adopt a radial basis function (RBF) kernel with bandwidth determined by the median heuristic as in previous studies. The "driving force" term is a smoothed gradient of the log posterior density that moves the particles toward the high-density regions of the posterior distribution. The "repulsive force" term promotes diversity and prevents particles from concentrating on the mode of the target PDF. This combination of two forces results in an efficient non-parametric approximation of the posterior PDF, using a finite number of particles. In an SVGD optimization, a mini-batch stochastic gradient descent can be used [22] for efficient optimization with large training datasets. In BNN applications, NNs corresponding to the particles \(\{\mathbf{\theta}\}_{i=1}^{n}\) are trained simultaneously by SVGD. SVGD can be applied to the naive Bayesian formulation of the target inverse problem introduced in the previous section by taking \(\mathbf{\theta}=(\mathbf{\theta}_{T}\,\mathbf{\theta}_{v})^{\top}\). We call this naive approach "SVGD for PINN-based seismic tomography" (SVGD-PINN-ST). Algorithm 1 summarizes the above baseline SVGD-PINN-ST algorithm with a mini-batch stochastic descent.

## III Formulation of velocity-space SVGD for PINN-based seismic tomography

### _Formulation of Bayesian estimation exclusive for velocity NN_

Although SVGD is applicable to BNN or Bayesian PINN theoretically, the parameter space \(\mathbf{\theta}=(\mathbf{\theta}_{T}\,\mathbf{\theta}_{v})^{\top}\) of a practical problem setting is usually high-dimensional (above \(10^{4}\) parameters), making the problem intractable for both HMC and SVGD. Consequently, reducing the effective parameter space dimension in the target Bayesian estimation is essential. For this reason, we seek to reduce the target space from \((\mathbf{\theta}_{T}\,\mathbf{\theta}_{v})^{\top}\) to \(\mathbf{\theta}_{v}\), considering that we are exclusively interested in the uncertainty of \(v\). The only requirement for \(\mathbf{\theta}_{T}\) is to minimize the residuals of the eikonal equation. We reformulate the Bayesian estimation defined in Equation 11 by considering exclusively the estimation of the posterior PDF of \(\mathbf{\theta}_{v}\) as follows:

\[P(\mathbf{\theta}_{v}|\mathbf{d}) \propto P(\mathbf{T}_{\rm obs}|\mathbf{\theta}_{v})P(\mathbf{\theta}_{v}), \tag{25}\]

with the condition

\[\mathbf{r}(\mathbf{\theta}_{T},\mathbf{\theta}_{v})=\mathbf{0}. \tag{26}\]

To estimate the posterior probability in Equation 25, we consider an SVGD update for \(\mathbf{\theta}_{v}\) only:

\[\mathbf{\theta}_{v\,i}^{l+1}=\mathbf{\theta}_{v\,i}^{l}+\epsilon_{l}\mathbf{\phi}\left(\mathbf{ \theta}_{v\,i}^{l}\right), \tag{27}\]

where

\[\mathbf{\phi}(\mathbf{\theta}_{v})\] \[=\frac{1}{n}\sum_{j=1}^{n}\{k(\mathbf{\theta}_{v\,j}^{l},\mathbf{\theta}_ {v})\nabla_{\mathbf{\theta}_{v\,j}^{l}}\log P(\mathbf{\theta}_{v\,j}^{l}|\mathbf{d})+ \nabla_{\mathbf{\theta}_{v\,j}^{l}}k(\mathbf{\theta}_{v\,j}^{l},\mathbf{\theta}_{v})\}. \tag{28}\]

To obtain \(\nabla_{\mathbf{\theta}_{v\,j}}\log P(\mathbf{\theta}_{v\,j}|\mathbf{d})\), the calculation of the total derivative is required, because the likelihood function is implicitly dependent on \(\mathbf{\theta}_{v}\) through \(\mathbf{\theta}_{T}\), as seen in Equation 13.
The above calculation can be performed using the Lagrange multiplier method introducing the equality constraints \(\mathbf{r}=\mathbf{0}\). This is a type of the discrete version of the adjoint method [27]. Here, we define the Lagrange function as the sum of the log likelihood function \(J=\log P(\mathbf{T}_{\mathrm{obs}}|\boldsymbol{\theta}_{v})\) and a Lagrange multiplier satisfying the above constraints, obtaining

\[J_{\mathrm{L}}=J+\boldsymbol{\lambda}^{\top}\mathbf{r}, \tag{29}\]

where \(J_{\mathrm{L}}\) and \(\boldsymbol{\lambda}\) are the Lagrange function and multiplier, respectively. The total derivative of \(J_{\mathrm{L}}\) with respect to \(\boldsymbol{\theta}_{v}\) can be calculated using the following chain rule of differentiation:

\[\frac{dJ_{\mathrm{L}}}{d\boldsymbol{\theta}_{v}} =\frac{\partial J_{\mathrm{L}}}{\partial\boldsymbol{\theta}_{v}}+ \frac{\partial J_{\mathrm{L}}}{\partial\boldsymbol{\theta}_{T}}\frac{d \boldsymbol{\theta}_{T}}{d\boldsymbol{\theta}_{v}}+\frac{\partial J_{\mathrm{L }}}{\partial\boldsymbol{\lambda}}\frac{d\boldsymbol{\lambda}}{d\boldsymbol{ \theta}_{v}}\] \[=\boldsymbol{\lambda}^{\top}\frac{\partial\boldsymbol{\tau}}{ \partial\boldsymbol{\theta}_{v}}+(\frac{\partial J}{\partial\boldsymbol{ \theta}_{T}}+\boldsymbol{\lambda}^{\top}\frac{\partial\boldsymbol{\tau}}{ \partial\boldsymbol{\theta}_{T}})\frac{d\boldsymbol{\theta}_{T}}{d\boldsymbol {\theta}_{v}}+\mathbf{r}^{\top}\frac{d\boldsymbol{\lambda}}{d\boldsymbol{ \theta}_{v}}. \tag{30}\]

The second and third terms on the right-hand side still include a total derivative. However, the third term vanishes if \(f_{T}\) is trained for \(f_{v}\) and \(\mathbf{r}\) becomes sufficiently small. The second term can also be eliminated by solving the following equation for \(\boldsymbol{\lambda}\):

\[\frac{\partial J}{\partial\boldsymbol{\theta}_{T}}+\boldsymbol{\lambda}^{\top} \frac{\partial\mathbf{r}}{\partial\boldsymbol{\theta}_{T}}=\mathbf{0}. \tag{31}\]

We approximate the solution of the last equation using the L-BFGS algorithm [32] using a zero vector as the initial solution. Equation 31 represents a discrete adjoint equation. The target total derivative is obtained by evaluating the first term in Equation 30, using the solution \(\boldsymbol{\lambda}^{*}\) satisfying Equation 31:

\[\frac{dJ_{\mathrm{L}}}{d\boldsymbol{\theta}_{v}}=\boldsymbol{\lambda}^{*\top} \frac{\partial\mathbf{r}}{\partial\boldsymbol{\theta}_{v}}. \tag{32}\]

Using this method, the gradient of the logarithm of the posterior PDF can be calculated as follows:

\[\nabla_{\boldsymbol{\theta}_{v}}\log P(\boldsymbol{\theta}_{v}| \mathbf{d}) = \nabla_{\boldsymbol{\theta}_{v}}\log P(\mathbf{T}_{\mathrm{obs}}| \boldsymbol{\theta}_{v})+\nabla_{\boldsymbol{\theta}_{v}}\log P(\boldsymbol{ \theta}_{v}) \tag{33}\] \[= \frac{dJ_{\mathrm{L}}}{d\boldsymbol{\theta}_{v}}+\nabla_{ \boldsymbol{\theta}_{v}}\log P(\boldsymbol{\theta}_{v}).\]

The discrete adjoint method is usually applied to regular systems, such as discretized PDEs. However, our adjoint equation is over-determined, in which the number of constraints (collocation points) is larger than that of the unknowns (weight parameters). Because PINN-based systems of PDEs are usually over-parameterized, the number of the former is significantly larger than the latter. For practical applications, considering that \(\frac{\partial J}{\partial\boldsymbol{\theta}_{T}}\) is a highly sparse vector due to the over-parametrization, Equation 31 is approximately satisfied by a moderate number of collocation points.

After each SVGD update for \(\boldsymbol{\theta}_{v}\) for all particles, \(\boldsymbol{\theta}_{T}\) is trained for a new \(\boldsymbol{\theta}_{v}\) satisfying \(\mathbf{r}=\mathbf{0}\) as follows:

\[\boldsymbol{\theta}_{T}^{*}=\mathop{\mathrm{arg\,min}}_{\boldsymbol{\theta}_{ T}}L(\boldsymbol{\theta}_{T}), \tag{34}\]

where

\[L = \frac{1}{N_{c}}\sum_{i=1}^{N_{c}}\left(f_{v}(\mathbf{x}_{c}^{(i)}, \boldsymbol{\theta}_{v})-\frac{1}{|\nabla f_{T}(\mathbf{x}_{c}^{(i)},\mathbf{x }_{s}^{(i)},\boldsymbol{\theta}_{T})|}\right)^{2} \tag{35}\] \[+ \frac{1}{N_{s}}\sum_{i=1}^{N_{c}}\left(f_{v}(\mathbf{x}_{s}^{(i)} ;\boldsymbol{\theta}_{v})-f_{\tau^{-1}}(\mathbf{x}_{s}^{(i)},\mathbf{x}_{s}^{( i)};\boldsymbol{\theta}_{T})\right)^{2}.\]

The resulting process is equivalent to a PINN-based forward solver of the eikonal equation [10, 11]. Algorithm 2 summarizes the improved SVGD algorithm in which the target parameter space is reduced to that of \(\boldsymbol{\theta}_{v}\) only.

### _Bayesian estimation of velocity NN with particle-based functional variational inference_

In the previous section, the parameter space for the Bayesian estimation was reduce to that of \(\boldsymbol{\theta}_{v}\). However, the NN for \(f_{v}\) is high-dimensional and over-parameterized as well, which result in a posterior PDF with multi-modal features in the weight space that impairs Bayesian estimation. To address this, we introduce a function-space SVGD (fSVGD or, more generally, fParVI) [24] to formulate the Bayesian estimation directly in the function space predicted by the associated NN. This reformulation remarkably improves the estimation accuracy of the posterior PDF in BNN by avoiding the direct exploration of a multi-modal PDF over the weight space [24]. We rewrite the equation for the Bayes' theorem as follows:

\[P(\mathbf{v}|\mathbf{d}) \propto P(\mathbf{T}_{\mathrm{obs}}|\mathbf{v})P(\mathbf{v}), \tag{36}\]

where \(\mathbf{v}\) is evaluated at \(N_{v}\) velocity evaluation points. The evaluation points are randomly chosen, similarly to collocation points:

\[\mathbf{v} = \left[f_{v}(\mathbf{x}_{v}^{(1)},\boldsymbol{\theta}_{v})\,f_{v}( \mathbf{x}_{v}^{(2)},\boldsymbol{\theta}_{v}),\ldots,f_{v}(\mathbf{x}_{v}^{(N_{v })},\boldsymbol{\theta}_{v})\right]^{\top}, \tag{37}\]

where the associated velocity dataset is defined as \(\mathbf{X}_{v}=\{(\mathbf{x}_{v}^{(1)},\mathbf{x}_{s}^{(1)}),(\mathbf{x}_{v}^{(2 )},\mathbf{x}_{s}^{(2)}),\ldots,(\mathbf{x}_{v}^{(N_{v})},\mathbf{x}_{s}^{(N_{v })})\}\). The source locations are used for the evaluation of the residual vector in the adjoint calculation. The update vector in Equation 27 is modified in fSVGD as follows:

\[\boldsymbol{\phi}\left(\boldsymbol{\theta}_{v\,i}^{l+1}\right)=\frac{\partial \mathbf{v}_{i}}{\partial\boldsymbol{\theta}_{v\,i}^{l}}\overset{\top}{\boldsymbol {\psi}}(\mathbf{v}_{i}^{l}), \tag{38}\]

where

\[\boldsymbol{\psi}(\mathbf{v})=\frac{1}{n}\sum_{j=1}^{n}\{k(\mathbf{v}_{j}^{l}, \mathbf{v})\nabla_{\mathbf{v}_{j}^{l}}\log P(\mathbf{v}_{j}^{l}|\mathbf{d})+ \nabla_{\mathbf{v}_{j}^{l}}k(\mathbf{v}_{j}^{l},\mathbf{v})\}. \tag{39}\]

The above rule calculates the update vector in the function (velocity) space based on SVGD and then converts it in the weight space using the Jacobian matrix. As the repulsive force is defined over the function space, the approximation of the posterior PDF does not suffer from the multi-modal feature present in the weight space formulation. Furthermore, as the functional prior \(P(\mathbf{v})\) is explicitly included, more physically meaningful prior information \(P(\boldsymbol{\theta}_{v})\) than that in the weight space can be introduced straightforwardly. Although velocity evaluation points can be changed in each epoch, the same points must be used for all the SVGD particles because the point locations define the PDF evaluated by SVGD. The calculation of \(\nabla_{\mathbf{v}_{j}^{l}}\log P(\mathbf{v}_{j}^{l}|\mathbf{d})\) requires the adjoint formulation \(\nabla_{\boldsymbol{\theta}_{v}\,j}\log P(\boldsymbol{\theta}_{v\,j}|\mathbf{d})\), as in the analysis presented in the previous section. In fact, the same formulation in Equations 29, 30 and 31 is applicable, simply replacing \(\boldsymbol{\theta}_{v}\) with \(\mathbf{v}\). These methods overcome the existing limitations in BNN for PINN-based inverse analyses in realistically complex problem
setting. As we focused on the velocity space to perform fsvgd, the final version of the proposed approach was called "velocity-space svgd for PINN-based seismic tomography" (vsvgd-pinn-st). Algorithm 3 summarizes the steps of vsvgd-pinn-st.

## IV Synthetic test in 1D tomographic problem

To verify the ability of the proposed method to perform Bayesian estimation, we apply vsvgd-pinn-st to synthetic 1D tomography tests in a simple problem setting and compare the results with those obtained using the baseline svgd-pinn-st and linear travel time tomography, which provides an analytical solution.

### _Linear traveltime tomography_

Linear travel time tomography estimates velocity perturbation from a reference model using a Taylor series expansion. We obtain, as a result, a linear inverse problem for the residual travel time. When a conjugate pair of the prior and posterior pdf is adopted, such as Gaussian distributions, Bayesian linear regression for linearized tomography provides an analytical solution for the posterior pdf (see Text S1 in Supporting Material). However, linear tomography neglects the dependence of the ray path (i.e., travel time in 1D problems) on velocity perturbations in the reference model. Consequently, the method accuracy degrades when the reference model does not offer an accurate approximation of the true one.

### _1D Synthetic Test_

In the 1D synthetic test (1DST), we set a simple true velocity model, with a constant velocity of 1 km/s in the 1D domain defined by \(0\leq x\leq 1.2\,\mathrm{km}\), to test the ueq performance (see Fig. 2). We employ such a simple structure because the focus here is the ability of ueq, not the estimation of velocity. In this configuration, the travel time is readily obtained. Ten points, serving both as receivers and sources, are evenly distributed in two regions defined by the intervals \(0.2\leq x\leq 0.4\,\mathrm{km}\) and \(0.8\leq x\leq 1\,\mathrm{km}\) using a 0.05 km spacing. We refer to the five points in each of the above intervals as Group 1 and Group 2, respectively. We only consider rays between points within the same group. Therefore the number of travel time data points is \(5\times 4+5\times 4=40\). No ray paths exist in the intervals \(0\leq x\leq 0.2\,\mathrm{km}\), \(0.4\leq x\leq 0.8\,\mathrm{km}\), and \(1\leq x\leq 1.2\,\mathrm{km}\), in which the uncertainty of velocity estimation is expected to be closer to that given by prior. \(\mathbf{E}_{\mathrm{obs}}\) was set assuming an i.i.d. zero-mean Gaussian noise distribution with standard deviation \(0.005\,s\), although we did not actually add artificial noise to the travel times, which were calculated analytically between points. The prior probability in the velocity space is represented as a stochastic process for Bayesian estimation because \(f_{v}\) is a continuous function of the 1D coordinate \(x\) in PINN, evaluated at arbitrary collocation points within the target domain. Hence, we use a Gaussian process with the same mean \(\mu(x)=1\,\mathrm{km}/\mathrm{s}\) as the true model and kernel function defined by:

\[k(x_{i},x_{j})=\sigma_{1}^{2}\exp\left(-\frac{1}{2\sigma_{2}^{2}}|x_{i}-x_{j} |^{2}\right). \tag{40}\]

This is an RBF kernel, where \(\sigma_{1}\) and \(\sigma_{2}\) are the standard deviations of the marginal probability and correlation length scale, respectively. We set \(\sigma_{1}=0.1\,\mathrm{km}/\mathrm{s}\) and adopt three different \(\sigma_{2}\) values, namely, \(0.25\,\mathrm{km}\), \(0.15\,\mathrm{km}\), and \(0.075\,\mathrm{km}\), for comparison.

In linearized tomography, which we consider to calculate the ground truth of the posterior probability, we set the reference model constant velocity at the true value of 1 km/s to achieve the best accuracy for the linear approximation. We divide the region using a 0.025 km spacing and estimate 48 unknowns parametrizing the velocity perturbation. For these parameters, a prior pdf is generated based on the Gaussian process defined above.

In vsvgd- and svgd-pinn-st, we use fully connected feed-forward neural networks for both \(f_{\tau^{-1}}\) and \(f_{v_{\mathrm{ph}}}\), setting \(v_{0}(\mathbf{x})=1\,\)km/s. The values of \(v_{\mathrm{ph}}^{\mathrm{max}}\) and \(v_{\mathrm{ph}}^{\mathrm{min}}\) are set to -0.4 and 0.4 km/s, giving the upper and lower limits of the velocity prediction \(f_{v}\) as 1.4 and 0.6 km/s, respectively. The Swish activation function [33] is applied in each layer, except in the output one where a linear activation is specified. We use three hidden layers for both \(f_{\tau^{-1}}\) and \(f_{v_{\mathrm{ph}}}\) with 50 and 10 hidden units, respectively. \(\mathbf{\theta}_{v}\) is initialized using the He's method [34]. \(\mathbf{\theta}_{T}\) is trained using the initialized \(\mathbf{\theta}_{v}\) according to Equations 34 and 35, before running the algorithms. For both algorithms, 256 SVgd particles are employed and the Adam optimizer [31] is used to determine \(\epsilon_{l}\). The travel time batch data size for \(\mathbf{X}_{2}^{b}\) (see line 2 in Algorithm 3) is set to 40 (i.e., full batch). For practical convenience, the coordinate data of the velocity evaluation points \(\mathbf{X}_{v}\) are generated in each iteration by random sampling in the target domain. The associated batch data \(\mathbf{X}_{v}^{b}\) are taken as the corresponding \(\mathbf{X}_{v}\) values. The data size is assigned to 200; the number of epochs for each iteration \(l\) of vsvgd-pinn-st is set to 2,000. Each training session of \(\mathbf{\theta}_{T}\) (see line 6 in Algorithm 3) is conducted by using a L-BFGS algorithm [32] for 10 epochs with \(N_{c}=200\). The collocation points coordinates are generated randomly at each iteration, similarly to the \(\mathbf{X}_{v}\) selection process. The initial learning rate of the Adam optimizer is set to \(10^{-2}\). For the baseline naive svgd-pinn-st, introducing a prior probability in the weight space, equivalent to that in the velocity space in Equation 40, represents a challenging task. Hence, we use an i.i.d. zero-mean Gaussian distribution as the prior probability of \(\mathbf{\theta}=(\mathbf{\theta}_{T}\,\mathbf{\theta}_{v})^{\top}\) with variance \(\sigma_{\theta}^{2}\). Considering that a physically meaningful choice of \(\sigma_{\theta}\) values is difficult, we test the standard value \(\sigma_{\theta}^{2}=10^{0}\) (see, for instance, [19]) and a significantly large one \(\sigma_{\theta}^{2}=10^{2}\). The number of epochs for each iteration \(l\) and the initial learning rate of the Adam optimizer are set to 90,000 and \(10^{-3}\), respectively.

When the prior probability has a long correlation length (\(\sigma_{2}=0.25\,\mathrm{km}\)), the uncertainty estimation by vsvgd-pinn-st agrees well with the analytical solution of the linearized tomography (see Fig. 2 (a) and (b)). In the region including the sources and receivers, the standard deviation (\(1\sigma\)) of the posterior probability is small for both methods. This value increases where there is a lack of observations. However, it is always significantly smaller than the standard deviation of the prior probability due to the imposed spatial correlation constraint. With an intermediate and a short correlation length
(\(\sigma_{2}=0.15\,\mathrm{km}\) and \(0.075\,\mathrm{km}\)), the uncertainty estimation accuracy by vSVGD-PINN-ST slightly decreases, as shown by the increasing difference from the analytical solution (see Fig. 2 (c)(e) and (d)(f) for linearized tomography and vSVGD-PINN-ST, respectively). In particular, for \(\sigma_{2}=0.15\,\mathrm{km}\) the uncertainty estimated by vSVGD-PINN-ST is close to the analytical solution at both boundaries of the region. However, simulated results deviate from analytical values in the central region. For \(\sigma_{2}=0.075\,\mathrm{km}\), although the estimated standard deviation in the central region is accurate, the curve is smoother than that of the analytical solution on the whole domain. This feature may be also related to a slightly overestimated uncertainty at the boundaries. The prior probability with a short correlation length assumes the existence of shorter-wavelength components in the target velocity structure. The NN architecture adopted in this analysis fails to learn some of these components possibly due to some learning bias, such as the spectral bias [35] that affects the ordinary PINN formulation based on fully connected feed-forward NNs. A learning bias may also have slightly compromised the uncertainty estimation accuracy of the experiments with the two lowest \(\sigma_{2}\) values. However, we consider the discrepancy in the estimated uncertainty acceptable for practical use at present. Further advances, in fact, would require improved NN architectures. Alternatively, a different configuration of the loss function may reduce the effect of the learning bias. This aspect is discussed in Section VI-A, as a future development.

We examine the relationship between UQ results and two other parameters, namely, the number of particles and epochs, focusing mainly on the experiment with \(\sigma_{2}=0.25\,\mathrm{km}\). Surprisingly, even when the number of SVGD particles (\(n\)) is small (e.g., 32), the estimated posterior mean and standard deviation still capture the basic features of the analytical solution (see Fig. 3(a)). This result demonstrates the advantage of using SVGD in the approximation efficiency over Bayesian sequential sampling methods. However, we observe an asymmetry in the shape of the curve of the standard deviation when \(n\) is small, showing a imperfect convergence near Group 2. The asymmetry disappears when \(n\) is increased. Increasing the number of epochs (or iterations) makes the UQ results converge from the side of larger standard deviations (see Fig. 3 (b)). This means that the SVGD update started with sufficiently diverging particles. Interestingly, the standard deviation curves of the first epochs show asymmetric patterns similar to the one seen in the results where a smaller \(n\) is used. Another interesting finding is that the same analysis applied to the experiment with \(\sigma_{2}=0.15\,\mathrm{km}\) shows faster convergence than the previous case (see Fig. S1 in Supporting Material). These results suggest SVGD requires more iterations to converge for the prior probability distributions with longer correlation length scale.

SVGD-PINN-ST completely failed to reproduce the expected spatial variation features of the posterior probability for both of the assumed prior PDFs in the weight space (the one with \(\sigma_{\theta}^{2}=10^{0}\) and the other with \(\sigma_{\theta}^{2}=10^{2}\)), showing no correspondence between these features and the locations of sources and receivers (see Fig. 4 and Fig. S2 in Supporting Material). These findings suggest two possible reasons explaining the poor performance of SVGD-PINN-ST in the estimation of the posterior probability. The first one is that the parameter space for the Bayesian estimation using the naive SVGD approach is too broad and multi-modal; the second one is that unfavorable effects on the prior probability in the weight space are not physically interpretable. Estimated SVGD-PINN-ST mean values are close to the true on. However, they show a larger discrepancy than those calculated by vSVGD-PINN-ST result.

In summary, vSVGD-PINN-ST can estimate the uncertainty in tomographic velocity estimation problems more accurately than the baseline SVGD-PINN-ST method. We confirmed that improvements in vSVGD-PINN-ST, which avoid direct Bayesian estimation in high-dimensional and multi-modal weight spaces, are essential to conduct an efficient Bayesian estimation in PINN-based inversion analyses. We also found that, even when the spatial correlation length of the velocity is short, vSVGD-PINN-ST returns competitive uncertainty estimates.

## V Synthetic tests in 2D tomographic problems

In this section, we test the applicability of the vSVGD-PINN-ST method to realistic 2D synthetic problems in two scenarios: surface-wave and refraction tomography.

### _2D Synthetic Test 1: surface-wave tomography_

The 2D synthetic test 1 (2DST1) is a synthetic test in 2D surface-wave tomography, in which a low velocity anomaly is surrounded by sources and receivers aligned in a circular arrangement. The purpose of this test is to show that the proposed vSVGD-PINN-ST algorithm gives a reasonable estimation of the velocity and its uncertainty, consistent with characteristic ray paths drawn using this circular configuration.

We use a true velocity model with a homogeneous background set at 2 km/s, containing a circular low velocity anomaly with 1.2 km/s at the center. Sixteen receivers are evenly distributed around the anomaly at a radius of 4 km (see Fig. 5 (a)). Each receiver also serves as a source point, thus, the number of travel time data points is \(16\times 15=240\). This configuration is inspired by the one in [3, 4]. However, we use connected background and anomaly regions to comply with the Gaussian prior distribution introduced later, which supports a NN-predicted continuous velocity function. We calculate the travel time for each source-receiver pair using the fast sweeping method [36] implemented with the tcrpy python package [37]. We add independent zero-mean Gaussian noise with standard deviation 0.01 s to the calculated travel times and use the resulting values as the observation data. \(\mathbf{E}_{\mathrm{obs}}\) is set accordingly, assuming that the error distribution of the travel time observations is known. The ray paths drawn using the numerical solution are absent at the center of the anomaly and outside of the sources and receivers (see Fig. 5 (b)). An accurate Bayesian estimation should infer a larger standard deviation for the estimated velocity in these regions compared to that in the areas where the ray paths go through. We adopted a Gaussian process as prior probability, defined by the kernel function in Equation 40 with the 2D L2 norm as
argument, mean \(\mu(\mathbf{x})=1.75\,\mathrm{km/s}\), and standard deviations \(\sigma_{1}=0.5\,\mathrm{km/s}\), and \(\sigma_{2}=1.5\,\mathrm{km}\).

In vSVGD-PINN-ST, we use fully connected feed-forward neural networks for both \(f_{\tau^{-1}}\) and \(f_{v_{\mathrm{pth}}}\), with \(v_{0}(\mathbf{x})=\mu(\mathbf{x})\). \(v_{\mathrm{pth}}^{\mathrm{max}}\) and \(v_{\mathrm{pth}}^{\mathrm{min}}\) are set to -1.5 km and 1.5 km, giving the upper and lower limits of the velocity predicted by \(f_{v}\) as 3.25 km/s and 0.25 km/s, respectively. For both networks, we apply the Swish activation function [33] in each layer except in the output one where a linear activation is specified. Six hidden layers are used for both \(f_{\tau^{-1}}\) and \(f_{v_{\mathrm{pth}}}\) with 50 and 20 hidden units, respectively. The number of SVGD particles employed for this experiment is 512. The Adam optimizer [31] with an initial learning rate of \(10^{-2}\) is used to determine \(\epsilon_{l}\). Batch sizes for \(\mathbf{X}_{T}^{b}\) and \(\mathbf{X}_{v}^{b}\) (see line 2 in Algorithm 3) are 240 (i.e., full batch) and 400, respectively. \(\mathbf{\theta}_{v}\) is initialized using a prior Gaussian process, obtained by using vSVGD-PINN-ST with zero weight on the travel time observation data. \(\mathbf{\theta}_{T}\) is trained using the initialized \(\mathbf{\theta}_{v}\) according to Equations 34 and 35 before the vSVGD-PINN-ST algorithm is applied. The number of epochs for iteration \(l\) is set to 800. Each training step of \(\mathbf{\theta}_{T}\) (see line 6 in Algorithm 3) is conducted by using the L-BFGS algorithm [32] for 50 epochs with \(N_{c}=1,600\), taking the previous SVGD iteration result as the initial guess. To prevent the solution from being trapped into local minima due to an ineffective initial guess, the training of \(\mathbf{\theta}_{T}\) is conducted from scratch for selected epochs (i.e., epoch number 50, 100, 200, and 400). The coordinates of the velocity evaluation and collocation points are randomly generated within the target domain and \(\mathbf{X}_{v}^{b}\) and \(\mathbf{X}_{c}^{b}\) are set to the corresponding \(\mathbf{X}_{v}\) and \(\mathbf{X}_{c}\) values, respectively.

The mean velocity model estimated by vSVGD-PINN-ST agrees well with the true value in the region where ray paths are present (see Fig. 6 (a) and (b)). In the regions outside the sources and receivers circumference, the mean is close to the true value probably because of the assumption on spatial correlation. The posterior PDF standard deviation values in these regions smoothly approach that of the prior Gaussian process (see Fig. 6 (c)). The estimated values in the central region reflect the absence of ray paths and are, as expected, larger than those in the surrounding areas where the ray paths go through. However, overall values are smaller than those of the prior uncertainty. This is primarily due to the assumption on the spatial correlation of the prior probability. In addition, the absence of ray paths in such a configuration does not allow too large velocity in this region, providing an additional indirect constraint. The mean values in the central area are similar to the true ones, suggesting that the true model used is highly consistent with the prior probability. To further analyze results, we compared the marginal probability distributions of three points (black marks in Fig. 5 (b)). Histograms in Fig. 7 show the difference in the uncertainty of the estimated velocity in the regions with and without ray paths (marked by a circle and a square or an inverse triangle, respectively). These results suggest that tomographic results obtained by vSVGD-PINN-ST are consistent with the true values and prior probability, and the ray paths distribution for the true velocity model.

To prove that vSVGD-PINN-ST improves UQ accuracy in a realistic problem setting, we conducted the same analysis using the naive method, SVGD-PINN-ST. Similarly to the SVGD-PINN-ST analysis in 1DST, we use i.i.d. zero-mean Gaussian prior PDFs as weight parameters (see Text S2 in Supporting Material for other details), taking \(\sigma_{\theta}^{2}=10^{0}\). Results show that the velocity models estimated by the SVGD-PINN-ST are not consistent with the true model and ray paths distributions (see Fig. S3 in Supporting Material), reconfirming that SVGD-PINN-ST could not estimate the posterior probability in seismic tomographic problems effectively.

### _2D Synthetic Test 2: seismic refraction tomography_

In this section, we present the results of the 2D synthetic test 2 (2DST2), in which we consider a refraction tomographic problem where sources and receivers are located on the Earth's surface. We use refracted rays penetrating underground to estimate 2D velocity profiles. Compared to 2DST1, this tomographic problem is ill-posed, given the distribution of sources and receivers [38]. Furthermore, there exist only a few studies on Bayesian seismic tomography for this type of setting [5]. However, in the following, we show that our vSVGD-PINN-ST algorithm is suitable for this type of tomographic problem as well.

We set a relatively simple true velocity model with an almost depth-dependent structure in a domain ranging from 0 to 30 and from 0 to 200 km in depth and horizontal distance, respectively. A velocity bump is included in the region between 70 and 130 km in the horizontal axis. We consider 20 sources and 96 receivers distributed with a uniform spacing between 5 and 195 km in the horizontal distance on the surface of the model (see Fig. 8 (a)). Therefore, the total number of travel time data is \(20\times 96=1920\). As in 2DST1, we create synthetic data by calculating the travel time for each source-receiver pair using the fast sweeping method and add the i.i.d. zero-mean Gaussian noise with a standard deviation of 0.05 s. \(\mathrm{E_{obs}}\) is set accordingly, assuming that the error distribution of traveltime observations is known. The ray paths drawn by using the numerical solution are densely distributed in a shallow portion of the target region (see Fig. 8 (b)). For the prior probability, we set a Gaussian process with an RBF kernel, similarly to the previous tests. In refraction tomography, we assume that velocity increases with increasing depth. Therefore, we adopted a mean velocity model \(\mu(\mathbf{x})\) that is linearly depth-dependent (see Fig. S4 in Supporting Material). Under such an assumption, the correlation length in the horizontal and vertical directions are likely to have significant differences. To reflect this prior information, we redefine the RBF kernel as follows:

\[k(\mathbf{x}_{i},\mathbf{x}_{j})=\sigma_{1}^{2}\exp\left[-\frac{1}{2}\left( \frac{(x_{i}-x_{j})^{2}}{\sigma_{x}^{2}}+\frac{(z_{i}-z_{j})^{2}}{\sigma_{z}^ {2}}\right)\right] \tag{41}\]

where \(\mathbf{x}=(x\,z)^{\top}\), and \(\sigma_{x}\) and \(\sigma_{z}\) are the correlation lengths in the horizontal and vertical direction, respectively. We set \(\sigma_{1}=1\,\mathrm{km/s}\), \(\sigma_{x}=35\,\mathrm{km}\), and \(\sigma_{z}=5\,\mathrm{km}\).

We use the same fully connected feed-forward neural networks as in the previous experiment, with six hidden layers for both \(f_{\tau^{-1}}\) and \(f_{v_{\mathrm{pth}}}\), including 50 and 20 hidden units, respectively. We set \(v_{0}(\mathbf{x})=\mu(\mathbf{x})\) as the mean velocity of the
prior Gaussian process introduced in the following. \(v_{\rm ptb}^{\rm max}\) and \(v_{\rm ptb}^{\rm min}\) are set to -3 km and 3 km, respectively, giving the depth-dependent upper and lower limits of the velocity predicted by \(f_{v}\). The Adam optimizer with an initial learning rate of \(10^{-3}\) is used to determine \(\epsilon_{l}\) in both algorithms. Weight parameters are initialized in the same way as in 2DST1. The batch sizes of \(\mathbf{X}_{T}^{b}\) and \(\mathbf{X}_{v}^{b}\) (see line 2 in Algorithm 3) are 480 and 400, respectively. The number of epochs for each iteration \(l\) is set to 200. Each training session of \(\boldsymbol{\theta}_{T}\) (see line 6 in Algorithm 3)is performed by the L-BFGS algorithm [32] for 50 epochs with \(N_{c}=1,600\). The initial guess was given in the same way as in 2DST1, except for the epochs 12, 25, 50, and 100 that were trained from scratch. The coordinates of the velocity evaluation and collocation points were randomly generated within the target domain, assigning to \(\mathbf{X}_{v}^{b}\) and \(\mathbf{X}_{c}^{b}\) the corresponding values in \(\mathbf{X}_{v}\) and \(\mathbf{X}_{c}\), respectively.

The mean velocity model estimated by vSVGD-PINN-ST agrees well with the true value in the region with a high density of ray paths (see Fig. 9 (a) and (b)). Fig. 9 (c) shows that the standard deviation is generally small in the same regions, with some irregular increases marked in light blue that may reflect complex ray-path patterns and heterogeneous velocity structure. In the bottom and side region not reached by ray paths, it smoothly increases to around 1 km/s, the standard deviation of the prior probability. The results of the experiment are consistent with the true values and the prior probability, proving that our vSVGD-PINN-ST algorithm can estimate the posterior probability successfully, even in highly ill-posed problem setting.

To further analyze the results, we compared the estimated mean model, true model, prior mean, and marginal probability distributions along the two lines marked in Fig. 8 (b) (see Fig. 10). In both profiles, the true velocity model agrees well with the estimated mean or, at least, its values are included within a high frequency region in the shallow portion. The mean and true velocity begin to grow apart at about 15 and 20 km of depth in the 20 (dashed line) and 100 km (dash-dot line) horizontal lines, respectively. We hypothesize that this is due to the reduction in the number of ray paths as the depth increases. At 20 km in the horizontal distance, the uncertainty is small even in the depth between 15 and 20 km that is out of the ray paths coverage, probably due to the spatial correlation assumption in the prior probability. Below these depths, the frequency color maps show a broad distribution, with increasing uncertainty due to the lack of ray paths. We expected the estimated mean model to agree well with the prior one in this depth, since information is not obtained directly from data. However, the former does not approach the latter as the depth increases, it grows apart from it. As we hypothesize in Section IV-B, we attribute this finding to some learning bias affecting our NN architectures, such as the spectral bias [35]. A possible way to address this issue is discussed as a further development in Section VI-A.

## VI Discussion

### _Advantages and future developments of PINN-based Bayesian seismic tomography_

The test problem setup for seismic refraction tomography in 2DST2 is different from the two previous ones, presenting the same observation geometry as that of real subsurface structural exploration. The successful performances of our method in this test confirm its applicability to actual observational data in geophysical exploration and subsurface structural studies in seismogenic zones. Existing seismic tomographic methods used in these research domains consider an arbitrary initial model to calculate the theoretical travel time and, subsequently, they update the velocity model using an iterative method to minimize the residual between the theoretical and observed travel time. The final velocity model is obtained when the residual becomes sufficiently small. A Monte Carlo (MC) analysis with initial model randomization is used to evaluate the uncertainty and reliability of the solution obtained using such tomographic methods [39, 40]. The idea behind the MC analysis is that the main source of uncertainty in the estimation results of traditional tomographic methods lies in the choice of the initial model. The introduction of Bayesian estimation in the proposed method mitigates the dependence of the estimation upon the initial value. The uncertainty evaluation reflects the quality and quantity of the observation data instead. In conventional tomography methods, model parameterization (quantity and arrangement selection) and regularization (e.g., smoothing parameters selection) are often determined subjectively. In this study, model parameterization is performed controlling the continuous functions represented by NNs, using the prior probabilities introduced by Bayes' theorem in the form of stochastic processes. The selection of regularization parameters is performed choosing the ones that characterize prior probabilities in Bayesian statistics, such as the correlation distance in the stochastic process adopted in our method. Although this aspect is not considered in this study, the most suitable parameters can be objectively determined within the framework of Bayesian statistics (e.g., using hierarchical or empirical Bayesian frameworks, see [53]). This aspect represents an important future development.

PINN-based Bayesian seismic tomography benefits from the general advantages of PINN methods for the solution of PDEs and deterministic inversion problems (see [41], for instance). The PINN-based approach we used does not require, in fact, a mesh or grid for numerical calculation because the continuous functions defined by NNs are differentiable and automatically generated in the domain of interest. Moreover, neural networks, automatic differentiation, and optimization algorithms can be introduced easily using existing deep-learning frameworks, such as Pytorch and TensorFlow. These functions also contributed to the efficient implementation of the SVGD algorithm and the adjoint method required in vSVGD-PINN-ST. The extension of this approach to higher dimensional domains (i.e., from 1D to 2D and from 2D to 3D) is easier than in ordinary numerical simulations since the mesh-free framework reduces the dependency of the algorithm on the dimension of the target problem. In fact, our Pytorch
code for vSVGD-PINN-ST required the modification of less than 100 lines of code to switch from 1D to 2D analyses.

The traditional PINN formulation with fully connected feed-forward NNs shows poor performance when the target functions include high-frequency or multi-scale features [42]. This is due to the NN spectral bias [35], which may be also responsible for the slight degradation of the estimation accuracy when the correlation length is small in 1DST and the discrepancy between the estimated and prior mean velocity at the bottom of the target domain in 2DST. Further improvements, such as the introduction of random Fourier feature methods [42, 43] and loss functions that account for physical causality [44], have been reported to tackle the above learning issues. We expect that a combination of these improvements will improve the UQ of our PINN-based seismic tomography method, resulting in increasingly realistic velocity structures.

All calculations performed for vSVGD-PINN-ST were accelerated using full parallelization for each SVGD particle. For instance, performing 2DST1 took 15.9 hours, using 512 CPU cores (64-core AMD EPYC 7742 \(\times\) 8 in Earth Simulator 4, made available by the Japan Agency for Marine-Earth Science and Technology (JAMSTEC)) assigned to each one of the 512 particles. The 2DST1 experiment mimics the setting of the synthetic test by Zhang & Curtis in [4], where the authors analyzed the computational cost required to execute the several Bayesian seismic tomography methods compared in their study. They found that analyses based on the Metropolis-Hastings Markov Chain Monte Carlo (MCMC) [45, 46] and reversible jump MCMC (rj-MCMC) method [1, 47] required 80.05 and 17.1 calculation hours, respectively, using six CPUs. Further parallelization of these sampling algorithms is difficult due to their sequential features. These considerations suggest that vSVGD-PINN-ST is competitive with some existing methods in terms of time-to-solution because of the high parallelism of the algorithm (note that Zhang & Curtis proposed highly efficient variational inference methods [4]). Furthermore, the estimation using our PINN-based method provides mesh-free continuous velocity models (such as those of individual SVGD particles shown in Fig. S5 in Supporting Material), whereas those obtained in [4] are parametrized with a relatively coarse spatial grid, which reduces computation cost. Incorporating finer grids will impose significantly larger cost on these other methods. Most of the vSVGD-PINN-ST computational time is required by the PINN-based solution of the eikonal equation in each SVGD update (in line 6 of Algorithm 3). Hyper parameters affecting this calculation time, such as the number of hidden layers in the NNs, hidden units in each NN layer, epochs, and collocation points, were not optimized in this study. To improve the SVGD efficiency, for instance, reducing the number of particles and iterations, and using improved kernels [48, 49] and second-order methods [50, 51] will be considered for future developments. For large-scale problems, requiring increasing numbers of NNs hidden layers and units, GPUs efficient acceleration might prove more advantageous for our PINN-based method compared with other ones based on ordinary numerical calculation.

To verify the accuracy of the PINN-based travel time calculation in vSVGD-PINN-ST, we compared the final travel time predicted by \(f_{T}\) with the result of the fast sweeping method, which uses the final velocity predicted by \(f_{v}\). Such comparison with reference solutions obtained from ordinary numerical simulation methods is currently the only available method to check the convergence to the true solution. Further theoretical studies are required to understand convergence properties of PINN-based solutions [41].

### _Prior probability_

Compared to the 2D surface wave tomography targeted in 2DST1, 2D seismic refraction tomography in 2DST2 is more likely to lead to a highly ill-posed inverse problem due to the distribution of sources and receivers [5, 38]. Leveraging prior information to provide proper constraints is crucial in this type of seismic tomography. In 2DST2, we introduced a mean function and correlation length dependent on depth and direction, respectively, in the Gaussian process with a RBF kernel used as the prior probability. This choice is based on the solid Earth science knowledge that seismic velocity structures are nearly horizontally stratified. The values of the parameters used in the kernel of the prior Gaussian process are reviewed in Table I. In our vSVGD-PINN-ST, such prior constraints can be incorporated directly because the Bayesian inference is performed in velocity space. In previous studies on Bayesian PINN, the estimation was performed in weight space using a simple prior probability, such as the i.i.d. Gaussian distribution, for the weight parameters [19, 21]. To understand how such a simple prior probability defined in weight space behaves in the physical space, we examine the corresponding Gaussian process in the function space to this prior. We draw 1,000 random samples of the weight parameter set of the NN predicting velocity used in 2DST2, using the i.i.d. Gaussian prior distribution, and we generate 1,000 corresponding predicted velocity structures. Subsequently, we find the best fitting \(\sigma_{1}\), \(\sigma_{x}\) and \(\sigma_{z}\) values in Equation 41 for the 1,000 velocity models (see Text S3 in Supplementary Material for further details). In Table I, two examples of the estimated parameter sets for the kernel function are presented. For instance, when the standard deviations of the i.i.d. Gaussian distributions is set to \(\sigma_{\theta}=10^{-0.6}\), the marginal standard deviation (\(\sigma_{1}\)) has a similar value to that of the prior distribution in 2DST. However, the estimated horizontal correlation length (\(\sigma_{x}\)) is larger than the horizontal domain size, which is physically not appropriate. In contrast, when \(\sigma_{\theta}=10^{-0.4}\), although the correlation length appears to be within an acceptable range, the marginal standard deviation (\(\sigma_{1}\)) is so large that the constraint from the depth-dependent mean function becomes too weak. These two examples demonstrate that introducing a physically interpretable prior probability is not an easy task when Bayesian PINN is performed in the weight space and a simple function such as i.i.d. Gaussian is adopted. A learning algorithm has been proposed in [52] to obtain the prior probability in the weight space from a target one defined in the function space. However, the high computational complexity of the learning algorithm is inevitably problematic when the problem size becomes large. Introducing SVGD in the function space for Bayesian PINN is an effective solution
not only because it overcomes the multi-modality issues in the weight space, but also because it results in a physically interpretable Bayesian inference.

When comparing different methods, we should always consider that a different parametrization leads to a different prior probability, resulting in a consequent discrepancy in the posterior probability. In [4], for instance, the results obtained with methods based on adaptive parametrization (i.e., the rj-MCMC method [1, 47]) were significantly different from those based on a fixed parametrization. Similarly, even if we had considered the exact same problem setting as the one in [4], performing a meaningful comparison would have been difficult because the seismic velocity obtained by the PINN-based tomography is also adaptively parameterized. Therefore, comparison of Bayesian seismic tomography methods in standardized problem settings with equivalent prior probabilities represents an important aspect to be considered for future developments.

Previous studies claim that the PINN-based approach can effectively solve ill-posed inverse problems without introducing prior knowledge for the target parameters estimation [14, 16]. From the Bayesian viewpoint, a PINN-based inverse analysis without prior constraints can be interpreted as a maximum a posteriori (MAP) estimation of the weight parameters, incorporating a prior PDF with uniform distribution of a wide value range (e.g., an improper flat prior). If we focus on the space of the function predicted by a NN, such a prior PDF is not a "non-informative prior" any more, because it imposes implicit prior information in the function space, obtained by a nonlinear transformation from the weight space (see [53] for uniform distributions defined by a nonlinear variable transformation from a different feature space). As a result, a MAP solution obtained without introducing explicit prior constraints may include effects from implicit prior information that are nonnegligible. As proposed in this study, explicitly providing physically interpretable prior information defined in the function space helps in addressing this issue hidden in ordinary PINN-based inverse analyses.

### _Application of vSVGD-PINN-ST to other PINN-based inversion problems_

The proposed vSVGD-PINN-ST algorithm can be applied to general PINN-based inversion problems that incorporate two NNs, one predicting the solution of the governing equation and the other its parameters (e.g., the full waveform inversion [16]). Some of the previous studies on Bayesian PINN targeted simple problems in which the solution of the governing equation predicted by a single NN is the only target of Bayesian estimation (e.g., [20, 23]). The proposed approach can also be used for this type of problems by introducing a simplified method that we call "function-space SVGD for PINN" (fSVGD-PINN). fSVGD-PINN can be derived from vSVGD-PINN-ST simply incorporating fSVGD in PINN as described in Section III-B, removing the procedure used to separate one of the two NNs from the Bayesian estimation described in Section III-A.

## VII Conclusion

In this study, we developed the vSVGD-PINN-ST algorithm, which performs PINN-based Bayesian seismic tomography using SVGD, the best-known particle-based variational inference method, applied only in the velocity space and enhanced with several mathematical and numerical techniques. The vSVGD-PINN-ST performance was tested in one- and two-dimensional Bayesian seismic tomography synthetic tests. Such problems cannot be handled by naive baseline algorithms that perform SVGD in the weight space of the component NNs predicting velocity and travel time. Results show that our method not only allows for accurate UQ but it can also incorporate physically-interpretable prior probability defined in the velocity (function) space, overcoming existing limitations of traditional BNN approaches based on Bayesian estimation in the weight space. To the authors' best knowledge, this is the first success in PINN-based Bayesian seismic tomography with practical estimation accuracy. The success of the last synthetic test adopting a realistic observation geometry, similar to a subsurface structural exploration, suggest that our method can be applied to actual observational data in geophysical exploration and subsurface structural studies. Finally, the proposed method offers a new fundamental Bayesian approach that can be applied to inverse problems sharing the same formulation, in geoscience and other fields, leveraging on the flexibility and extendibility of PINN.

## Acknowledgments

Comments from Dr. Tatsu Kuwatani were valuable for designing the 1D synthetic test. This research was supported by JSPS KAKENHI Grant Number 21K14024. Computational resources of the Earth Simulator 4 provided by JAMSTEC was used.

```
0: A set of initial particles \(\{\mathbf{\theta}^{0}\}_{i=1}^{n}\), where \(\mathbf{\theta}=(\mathbf{\theta}_{T}\,\mathbf{\theta}_{v})^{\top}\).
0: A set of particles \(\{\mathbf{\theta}\}_{i=1}^{n}\), which approximates the target distribution \(P(\mathbf{\theta}|\mathbf{d})\).
1:for iteration \(l\)do
2: Sample a mini batch \(\mathbf{X}_{T}^{b}\) and \(\mathbf{X}_{c}^{b}\) from training set \(\mathbf{X}_{T}\) and \(\mathbf{X}_{c}\).
3: For each \(i\in[n]\), calculate the SVGD update vector for \(\mathbf{X}_{T}^{b}\) and \(\mathbf{X}_{c}^{b}\) according to Equation 24.
4: For each \(i\in[n]\), calculate \(\mathbf{\theta}_{i}^{l+1}\) according to Equation 23.
5: Set \(l\gets l+1\)
6:endfor
```

**Algorithm 1** SVGD-PINN-ST, a naive approach based on SVGD in weight space.

## Appendix

Normalization strategy of NN input and output

In this section, we outline the methods used to normalize the input and output value of the model NNs to improve their convergence performance, using a similar approach to the one in [16]. In order to obtain normalized coordinates in
**Algorithm 2** An updated approach based on SVGD only in the weight space of velocity NN.

```
0: A set of initial particles \(\{\mathbf{\theta}_{v}^{0}\}_{i=1}^{n}\) and \(\{\mathbf{\theta}_{T}^{0}\}_{i=1}^{n}\) that are initially trained for \(\{\mathbf{\theta}_{v}^{0}\}_{i=1}^{n}\).
0: A set of particles \(\{\mathbf{\theta}_{v}\}_{i=1}^{n}\), which approximates the target distribution \(P(\mathbf{\theta}_{v}|\mathbf{d})\).
1:for iteration \(l\)do
2: Sample a mini batch \(\mathbf{X}_{T}^{b}\) and \(\mathbf{X}_{c}^{b}\) from training set \(\mathbf{X}_{T}\) and \(\mathbf{X}_{c}\).
3: For each \(i\in[n]\), calculate the SVGD update vector for \(\mathbf{X}_{T}^{b}\) according to Equation 28.
4: For each \(i\in[n]\), calculate \(\mathbf{\theta}_{v\,i}^{l+1}\) according to Equation 27.
5: For each \(i\in[n]\), update \(\mathbf{\theta}_{T\,i}^{l+1}\) for \(\mathbf{X}_{c}^{b}\) according to Equation 34 and 35.
6: Set \(l\gets l+1\)
7:endfor
```

**Algorithm 3** VSVGD-PINN-ST, the final version of the algorithm based on velocity-space SVGD.

```
0: A set of initial particles \(\{\mathbf{\theta}_{v}^{0}\}_{i=1}^{n}\) and \(\{\mathbf{\theta}_{T}^{0}\}_{i=1}^{n}\) that are initially trained for \(\{\mathbf{\theta}_{v}^{0}\}_{i=1}^{n}\).
0: A set of particles \(\{\mathbf{\theta}_{v}\}_{i=1}^{n}\), such that \(f(\mathbf{x};\mathbf{\theta}_{v\,i})\) approximates the target distribution \(P(v(\mathbf{x})|\mathbf{d})\).
1:for iteration \(l\)do
2: Sample a mini batch \(\mathbf{X}_{T}^{b}\), \(\mathbf{X}_{c}^{b}\) and \(\mathbf{X}_{v}^{b}\) from training set \(\mathbf{X}_{T}\), \(\mathbf{X}_{c}\) and \(\mathbf{X}_{v}\).
3: For each \(i\in[n]\), calculate the SVGD update vector for \(\mathbf{X}_{T}^{b}\) and \(\mathbf{X}_{v}^{b}\) in the velocity space according to Equation 39.
4: For each \(i\in[n]\), calculate the SVGD update vector in the weight space according to Equation 38.
5: For each \(i\in[n]\), calculate \(\mathbf{\theta}_{v\,i}^{l+1}\) according to Equation 27.
6: For each \(i\in[n]\), update \(\mathbf{\theta}_{T\,i}^{l+1}\) for \(\mathbf{X}_{c}^{b}\) according to Equation 34 and 35.
7:endfor
```

**Algorithm 4** VSVGD-PINN-ST, the final version of the algorithm based on velocity-space SVGD.

## Appendix B

Fig. 1: Neural network formulation adopted in this study and a schematic view of PINN-based (deterministic) seismic tomography.
Figure 2: The results of 1DST. (a)(b) Those obtained by using linearized tomography, which we consider as the ground truth, and \(\forall\)SVGD-PINN-ST, respectively, with \(\sigma_{2}=0.25\,\mathrm{km}\). (c)(d) Those with \(\sigma_{2}=0.15\,\mathrm{km}\). (e)(f) Those with \(\sigma_{2}=0.075\,\mathrm{km}\).
where \(a\) is the output of either \(f_{v_{\mathrm{pub}}}\) or \(f_{\tau^{-1}}\) and \(a^{\mathrm{max}}\) and \(a^{\mathrm{min}}\) represent the maximum and minimum range values defined in the following. For \(f_{v_{\mathrm{pub}}}\), \(a^{\mathrm{max}}=v_{\mathrm{pub}}^{\mathrm{max}}\) and \(a^{\mathrm{min}}=v_{\mathrm{pub}}^{\mathrm{min}}\), which are given a priori. For \(f_{\tau^{-1}}\), \(a^{\mathrm{max}}=v^{\mathrm{max}}\) and \(a^{\mathrm{min}}=v^{\mathrm{min}}\), where \(v^{\mathrm{max}}\) and \(v^{\mathrm{min}}\) given by \(v^{\mathrm{max}}=\max(v_{0}(\mathbf{x}))+v_{\mathrm{ptb}}^{\mathrm{max}}\) and \(v^{\mathrm{min}}=\min(v_{0}(\mathbf{x}))+v_{\mathrm{ptb}}^{\mathrm{min}}\), respectively. This operation imposes direct output values of the networks that are included in the interval \([-1,1]\), ensuring that the final upper and lower output limit values are determined by \(a^{\mathrm{max}}\) and \(a^{\mathrm{min}}\).
