Microseismic source imaging using physics-informed neural networks with
  hard constraints
Microseismic source imaging plays a significant role in passive seismic
monitoring. However, such a process is prone to failure due to aliasing when
dealing with sparsely measured data. Thus, we propose a direct microseismic
imaging framework based on physics-informed neural networks (PINNs), which can
generate focused source images, even with very sparse recordings. We use the
PINNs to represent a multi-frequency wavefield and then apply inverse Fourier
transform to extract the source image. To be more specific, we modify the
representation of the frequency-domain wavefield to inherently satisfy the
boundary conditions (the measured data on the surface) by means of a hard
constraint, which helps to avoid the difficulty in balancing the data and PDE
losses in PINNs. Furthermore, we propose the causality loss implementation with
respect to depth to enhance the convergence of PINNs. The numerical experiments
on the Overthrust model show that the method can admit reliable and accurate
source imaging for single- or multiple- sources and even in passive monitoring
settings. Compared with the time-reversal method, the results of the proposed
method are consistent with numerical methods but less noisy. Then, we further
apply our method to hydraulic fracturing monitoring field data, and demonstrate
that our method can correctly image the source with fewer artifacts.
Xinquan Huang, Tariq Alkhalifah
# Microseismic source imaging using physics-informed neural networks with hard constraints

###### Abstract

Microseismic source imaging plays a significant role in passive seismic monitoring. However, such a process is prone to failure due to the aliasing problem when dealing with sparse measured data. Thus, we propose a direct microseismic imaging framework based on physics-informed neural networks (PINNs), which can generate focused source images, even with very sparse recordings. We use the PINNs to represent a multi-frequency wavefield and then apply the inverse Fourier transform to extract the source image. Specially, we modify the representation of the frequency-domain wavefield to inherently satisfy the boundary conditions (the measured data on the surface) by means of the hard constraint, which helps to avoid the difficulty in balancing the data and PDE losses in PINNs. Furthermore, we propose the causality loss implementation with respect to depth to enhance the convergence of PINNs. The numerical experiments on the Overthrust model show that the method can admit reliable and accurate source imaging for single- or multiple- sources and even in passive monitoring settings. Then, we further apply our method on the hydraulic fracturing field data, and demonstrate that our method can correctly image the source.

## Plain Language Summary

Source imaging plays an important role in seismology and exploration geophysics as a reliable tool to locate earthquake events and monitor human-induced microseismic events. The conventional process numerically propagates the waveform recorded, for example, on the Earth surface backward in time so that energy focuses at the source location. However, when dealing with sparse recordings, this numerical implementation fails to focus the energy at the source, as it encounters considerable numerical artifacts and energy leakage. In addition, its computational cost is high. Neural networks have the potential to solve these problems because of their efficient execution and great interpolation abilities. Here, we propose to train a neural network (NN) to learn the observed data from imperfect recordings and then incorporate the prediction of this small NN into another NN, which learns the wavefield solutions, by means of a hard constraint. The functional form of the NN allows for source imaging in irregular geometry and the proposed method could provide reliable and robust source imaging for sparse recordings like those we deal with in global seismology.

## 1 Introduction

Microseismic events location is the essential foundation of passive seismic monitoring. The common way to locate seismic events is based on the travel time inversion (Lienert et al., 1986), which requires the time-consuming first arrival travel time picking. The process of picking arrivals becomes a challenge in low signal-to-noise data and in multiple-events cases. What's more, the resolution and accuracy of the source locations via this type of method are relatively low. Directly making use
of the full waveform information in the seismograms can help us avoid such issues and allow for higher resolution and accuracy in locating the microseismic events. Kao and Shan (2004) proposed an imaging approach based on summing the amplitudes of the measured data at the corresponding estimated arrival times without the need for phase picking, called the source scanning algorithm. It is widely used in seismology, e.g. for detailed earthquake rupture imaging (Ishii et al., 2007). Another type of method based on full waveforms is time-reversal imaging (McMechan, 1982). Its key idea is to backpropagate the waveform recorded on the surface until the energy is focused. Gajewski and Tessmer (2005) utilized the time reversal imaging (TRI) to reverse the waveform in the time domain for source locations. For this category of methods, the accuracy of the location results is highly determined by the imaging condition. In the past decade, many types of imaging conditions, including the maximum energy imaging condition (or the zero-lag cross-correlation imaging condition) (Artman et al., 2010; Yang and Zhu, 2019), the interferometric imaging condition for sparsely sampled data (Sava, 2011), the deconvolution imaging condition (Douma and Snieder, 2014) and the geometric mean imaging condition (Nakata and Beroza, 2016) for higher resolution, the maximum variance imaging condition (Wang and Alkhalifah, 2017) and the grouped seismic imaging condition (Lin et al., 2020; Chen et al., 2021) for noise-robust source location, have been proposed to enhance waveform-based location results. However, the accuracy for TRI is still an issue in case of sparse and irregular observations. The other type of waveform-based method is to invert the observed data for a more accurate velocity model (Kamei and Lumley, 2014; Kaderli et al., 2018; Song et al., 2019; Wang and Alkhalifah, 2021). These methods require massive repeated wavefield computation, and thus, they are costly.

Thanks to recent advances in computing and algorithms, machine learning has achieved a lot of success in many fields, e.g. natural language processing, computer vision, and science. Many computer-vision-based methods have been adapted to microseismic source location tasks (Perol et al., 2018; Kriegerowski et al., 2019; Zhang et al., 2020; Wang and Alkhalifah, 2021; Zhang et al., 2022; Chen et al., 2022). Most of these methods deal with source locations in a supervised manner. That is to say, they train a neural network on simulated data with labels (source locations) available and then minimize the loss using labels. However, the generalization of these methods to the field data is a challenge due to the lack of a sufficiently diverse dataset. Labels for the field data are often attained using conventional methods and human intervention, and thus, they are prone to errors. On the other hand, purely data-driven methods, which do not require labels, are an attractive alternative to the supervised approach. One way to devise such an approach is by incorporating the physics priors or specifically training the neural networks with physics constraints.

Based on the universal function approximation theory (Hornik et al., 1989), neural networks can be utilized to represent functions like the seismic wavefield, which is a key component of the source location process using time reversal imaging. Thus, we represent the seismic wavefield with a function of the spatial coordinates and frequency, yielding significant flexibility of the simulation in irregular domains. The governing equations for seismic modeling, are used as a loss function, to optimize the neural network using physics-informed neural networks (PINN framework (Raissi et al., 2019)). By means of PINNs for source imaging, the subsurface velocity information is well embedded into the neural network, and it can image the source with irregular and sparse receivers.

In this paper, we propose a novel direct microseismic source imaging method by means of PINNs, in which the source image is a snapshot of the time-domain wavefield. We represent the wavefield as a function of the spatial coordinates and frequency via a neural network (NN) and then train the NN with the Helmholtz equation as a loss function. Specifically, we incorporate the observed data on the surface (boundary condition for the Helmholtz equation) into the wavefield via the hard constraint and use the modified Helmholtz equation as the loss function. Besides, inspired by the fact that the information guiding the optimization of the NN parameters is coming from the recorded data on the surface, we propose to impose causality to the loss calculation along the depth axis, yielding an ideal from-surface-down reconstruction of the wavefield, which can improve the convergence speed and accuracy. We will first demonstrate the effectiveness of the proposed method via single- and multiple-passive sources scenarios on the Overthrust model using sparsely sampled data. We further test our proposed method on the Oklahoma Arkoma Basin Hydraulic Fracturing data to highlight the benefits and potential of the proposed method.

In summary, our main contributions are three-fold:

* We develop a direct source imaging framework based on PINNs with hard constraints.
* To ensure a stable and reliable training process, we use a reference frequency loss function incorporated in the hard constraint implementation and impose causality on the loss function with respect to depth.
* We evaluate the proposed method on the synthetic data for two different synthetic cases and field data, and achieve reliable source location results, showing the potential for even global event location.

We first introduce the proposed source imaging framework and key concepts including PINNs with hard constraints, reference frequency loss function, causality implementation, and data fitting. To demonstrate the effectiveness of the proposed method, Sections 3 and 4 present the numerical examples and field tests. Finally, we discuss the potential of the approach and its limitations in Section 5 and conclude in Section 6.

## 2 Theory

### Source imaging in the form of wavefield modeling

The source location imaging problem can be formulated using frequency-domain wavefield modeling, which offers a reduction of dimensionality compared to time-domain wavefield modeling, but more importantly, it will allow us to implement a causal loss function in the depth direction. The wave equation in the frequency domain is described by the Helmholtz equation and for a 2D acoustic isotropic medium, it is, along with the data boundary condition, given by

\[\frac{\omega^{2}}{v^{2}(x,z)}u(x,z,\omega)+\nabla^{2}u(x,z,\omega )=0, \tag{1}\] \[u(x,z=0,\omega)=D(x,\omega), \tag{2}\]

where \(u(x,z,\omega)\) is the wavefield corresponding to the frequency \(\omega\) and location \((x,z)\), \(v(x,z)\) is the velocity, and \(D(x,\omega)\) is the recorded data on the surface, which can be transformed to the frequency domain using Fast Fourier transform (FFT). We use Equations 1 and 2 to solve for the wavefield and then transform it to the time domain in which the source image corresponds to a slice of the time-domain wavefield where the energy focus is the highest.

However, frequency domain modeling for irregular geometry and high frequencies, where the fine grids are needed, are memory and computationally intensive. It also requires solving for every frequency. On the other hand, function learning, by utilizing an NN, is a potential solution in view of its flexibility to irregular mesh and its instant inference speed. Thus, we use PINNs to represent the frequency-domain wavefield.

### Physics-informed neural networks with hard constraints

In seismic wavefield modeling, physics-informed neural networks (PINNs) can be used to find \(\Phi(x,z,\omega,\theta)\) with parameters \(\theta\), that maps the input spatial coordinates and frequency to the value of the complex wavefield that satisfies the Helmholtz equation. Using a vanilla PINN often includes two loss terms for each of Equations 1 and 2. Balancing the contributions of these two terms using a weight affects the convergence and accuracy of the solution. The hard constraints (Berg and Nystrom, 2018; Schiassi et al., 2021) concept can alleviate this issue by modifying the representation of the wavefield \(u(x,z,\omega)\) to a new form where the boundary condition is inherently satisfied. Specifically, in this paper, we propose to use the following form:

\[u(x,z,\omega)=D(x,\omega)+z\Phi(x,z,\omega,\theta). \tag{3}\]

Note that the zero-depth (\(z=0\)) frequency-domain wavefield in this case is equal to the recorded surface data \(D(x,\omega)\). By means of the hard constraint, the governing equations are reduced to a single equation given by inserting Equation 3 into Equation 1, and thus, the loss function of PINNs is reduced to the mean square errors of the resulting single equation. Thus, the loss function with hard
constraints can be written as follows:

\[\begin{split}\mathcal{L}&=\frac{1}{N}\sum_{i=1}^{N} \left|\mathbf{z}^{i}\frac{(\omega^{i})^{2}}{(v^{i})^{2}}\Phi(\mathbf{x}^{i}, \mathbf{z}^{i},\omega^{i},\theta)+\mathbf{z}^{i}\nabla^{2}\Phi(\mathbf{x}^{i}, \mathbf{z}^{i},\omega^{i},\theta)+\right.\\ &\left.2\nabla_{z}\Phi(\mathbf{x}^{i},\mathbf{z}^{i},\omega^{i}, \theta)+\frac{(\omega^{i})^{2}}{(v^{i})^{2}}D(\mathbf{x}^{i},\omega^{i})+2 \nabla_{x}^{2}D(\mathbf{x}^{i},\omega^{i})\right|_{2}^{2},\end{split} \tag{4}\]

where \(N\) is the number of samples used for training, which includes the samples with variable spatial coordinates \((x^{i},z^{i})\) and angular frequencies \(\omega^{i}\).

### A modified loss function

Recall that if we do the source imaging in the frequency domain, then the wavefield belonging to multiple frequencies is needed for the transformation to the time domain wavefield. As shown in Huang and Alkhalifah (2022), direct use of the PINN with the loss (Equation 4) would decrease the accuracy and convergence speed. Similar to their implementation, we modify the loss function with a single reference frequency loss by replacing \(\omega\) to \(\alpha\omega_{ref}\) and \(U(x,z,\omega)\) to \(D(x,\omega)+\alpha z\Phi(x,z,\omega,\theta)\), yielding

\[\begin{split}\mathcal{L}&=\frac{1}{N}\sum_{i=1}^{N} \left|\alpha\mathbf{z}^{i}\frac{\omega_{ref}^{2}}{(v^{i})^{2}}\Phi(\mathbf{x} ^{i},\mathbf{z}^{i},\omega^{i},\theta)+\alpha\mathbf{z}^{i}\frac{\partial^{2} \Phi(\mathbf{x}^{i},\mathbf{z}^{i},\omega^{i},\theta)}{\partial(\alpha x)^{2}} \right.\\ &\left.+\alpha\mathbf{z}^{i}\frac{\partial^{2}\Phi(\mathbf{x}^{i},\mathbf{z}^{i},\omega^{i},\theta)}{\partial(\alpha z)^{2}}+2\frac{\partial \Phi(\mathbf{x}^{i},\mathbf{z}^{i},\omega^{i},\theta)}{\partial(\alpha z)}+ \right.\\ &\left.\frac{\omega_{ref}^{2}}{(v^{i})^{2}}D(\mathbf{x}^{i},\omega ^{i})+2\frac{\partial^{2}D(\mathbf{x}^{i},\omega^{i})}{\partial(\alpha x)^{2} }\right|_{2}^{2},\end{split} \tag{5}\]

where \(\alpha\) is a scaling factor equal to the ratio of the frequency \(\omega\) to the reference frequency \(\omega_{ref}\). With this scaling factor, the frequency variation is compensated by the scaling of the spatial dimensions, so that the wavelength within the wavefield apparently (to the NN) does not change much with frequency. Accordingly, the dimensions of the velocity \(v\) is also scaled, by \(\alpha\), to \(v(\alpha x,\alpha z)\). For details, we refer the reader to Huang and Alkhalifah (2022).

### Causality loss implementation

To improve the convergence of PINN, we replace the conventional loss function with a causality loss implementation. Since the information guiding the wavefield is coming from the data on the surface, we impose causality to the loss function along the depth axis, so that the wavefield is reconstructed from the earth surface down. This has an equivalence in imaging referred to as downward continuation. Like the loss calculation based on the temporal causality Wang et al. (2022), where the causality is with respect to time, here we can apply it with respect to depth. In this paper, we are formulating the solution with respect to an initial condition \(D(x,\omega)\) (for a solution in depth instead of time). The corresponding modified loss function is given by,

\[\mathcal{L}_{c}=\frac{1}{N_{z}}\sum_{i=1}^{N_{z}}\exp\left(-\epsilon\sum_{k=1 }^{i-1}\mathcal{L}\left(z_{k},\mathbf{\theta}\right)\right)\mathcal{L}\left(z_{i},\mathbf{\theta}\right), \tag{6}\]

where \(\mathcal{L}_{c}\) is the PDE loss with causality implementation, \(\mathcal{L}(z_{i}\theta)\) is the loss corresponding to the parameters \(\theta\) (equation 5) at a depth of \(z_{i}\), \(N_{z}\) is the number of samples along the depth axis, and \(\epsilon\) is the causality parameter. It could be a constant hyperparameter, but we modify it to be a varying one that varies with the number of epochs. We use \(\epsilon=\epsilon_{0}+\frac{\lambda}{epoch+1}\) as an alternative, where \(\epsilon_{0}\) is the initial value of \(\epsilon\) and \(\lambda\) is a factor to adjust the change in \(\epsilon\). The causality implies that, for example, the solution at depth 0.5 m depends on the solution at depth 0.3 m, but not the opposite, since the initial condition is at depth zero. With this modification, the convergence is faster and the NN admits better predictions, which we will see later.

### Fitting the observed data with NNs

As suggested in Schiassi et al. (2021), the hard constraint boundary condition could be stored in a neural network to simplify the calculation of the partial derivatives for the hard constraint PINN
loss later. As a result, we fit a small multi-layer perception (MLP) to the surface observed data in a supervised fashion, so it becomes as an NN function of the horizontal position of the receiver and frequency. Such training allows for fitting the NN to irregular acquisition (receivers) geometry, which is common in field applications (e.g. the lack of receivers or the non-uniform receiver intervals). The size of this data NN function will control the smoothness involved in the interpolation (Llanas and Sainz, 2006). More importantly, it allows for robust derivatives calculation for equation 5 using automatic differentiation (AD). Here, the data fitting branch is trained in a supervised manner, where the input of the NN are \(x\) and \(\omega\), and the output are the real and imaginary parts of the frequency-domain data at that location and frequency. Then, the supervised loss function for the data fitting is defined as follows:

\[\mathcal{L}_{d}(\theta)=\lambda\frac{1}{N_{r}\times N_{f}}\sum_{j=1}^{N_{r} \times N_{f}}\left|D\left(x_{j},\omega_{j}\right)-D_{\theta}\left(x_{j},\omega _{j}\right)\right|_{2}^{2}, \tag{7}\]

where \(N_{r}\) and \(N_{f}\) are the numbers of samples along the receivers and frequencies. Using all samples is not always needed because the noise in the samples may decrease the accuracy and we usually only need the use of a random subset of the data points to train the NN. Here, we show a simple demonstration to highlight the versatility of the approach for data fitting in handling irregularly sampled data. We generate data for a small layered model (same as the one used in Huang et al. (2021)) and place 100 receivers on the surface. We randomly drop receivers to imitate irregular data to train the network. As shown in Figure 1, the prediction looks reasonable even when using only 20% of the receivers for training. The second-order derivative w.r.t \(x\) of the NN data function calculated by AD, which is needed for equation 5, also has good accuracy compared with the second-order derivative calculated using a 5-point central difference formulation or AD with the full data.

### The pipeline and the backbone NN architecture

As shown in Figure 2, the preprocessed (e.g., non-local means filtering, NLM) recorded data on the surface are first transformed into the frequency domain via FFT and then used to train a small NN

Figure 1: The comparison between the NN function and the 3Hz observed data. a) shows the prediction of NN function trained with 20% to 90% of the dataset, and b) is the comparison between the corresponding second-order derivates calculated by AD or FD for various coverage percentages.
for data fitting. This trained NN will be used for the evaluation of the data term and calculation of the second-order derivate by AD in the partial differential equation (PDE) fitting. In the PDE fitting branch, the parameters of the data fitting branch are frozen, and another new larger NN is used to fit the PDE. To accelerate the convergence, we use the causality implementation of the PDE loss. After training, we predict the multi-frequency wavefield using the NN and utilize inverse FFT to get the time-domain snapshot. The locations where the energy focus are the source locations.

As for the implementation of the neural networks used in this pipeline, we utilize positional encoding [Huang et al., 2021] to help the network deal with the complex wavefield variations and help the convergence of the NN. The backbone of the NN is MLP (Multilayer Perception) with sine as the activation function. The size of the networks and training details will be shared in the examples.

## 3 Numerical examples

We first evaluate the proposed method on synthetic data including single-source event and multiple-source events cases. The synthetic data are generated for the Overthrust model (Figure 2(a)). The model size is 501\(\times\)161 with a spatial interval of 25 m. During the training for all synthetic examples, we use an MLP with four hidden layers, with 128, 64, 32, and 16 neurons in the layers from shallow to deep respectively, to fit the data, while for the PDE fitting, we use a larger MLP with six hidden layers, with 256, 256, 128, 128, 64, and 64 neurons from shallow to deep, respectively. We also include the positional encoding for \(x\), \(z\), and \(\omega\) with a level of 4. The frequency range for both tasks is from 3 to 12 Hz and we choose 12 Hz as the reference frequency.

### A single-source event case

In this subsection, we first test the proposed method in a toy problem, that has only one source event. The source is placed at (5.0, 3.0) km. We first do the forward modeling on the Overthrust model due to this source and recorded the data with 50 randomly placed receivers on the surface, shown in Figure 2(a). We train the NN for data fitting with positional encoding and train for 20000 epochs using an Adam optimizer with a learning rate of 1e-3. Then we train another NN for PDE fitting with 80000 random samples for 6000 epochs where we freeze the parameters of the NN for data fitting. We evaluate the NN on a uniform grid with 501\(\times\)161 grid points and an interval of 25 m to obtain the multi-frequency wavefields. As we have the source information for this case, we directly obtain the

Figure 2: The pipeline of the proposed method. There are two branches: the data fitting branch (**top**) and the PDE fitting branch (**bottom**). The recorded data \(d(x,t)\) are first transformed to the frequency domain and then be used to train the data NN. And then the data NN is combined into the PDE fitting branch to train the PDE NN. Then the output of the PDE NN are transformed to the time-domain snapshots using inverse FFT, where the red box denotes the source imaging where the energy focus. The module PE denotes the positional encoding with sinusoidal functions [Huang et al., 2021].
source image by the summation of the wavefield over frequencies (equal to the time-domain snapshot at zero time). As shown in Figure 3 b, the source location in the image is consistent with the ground truth (black star in Figure 3a), which demonstrates the accuracy of the proposed method. Next, we test more complex scenarios.

### A more complex case including triple-sources event

In a more realistic scenario, there might be multiple sources distributed at variable locations in the subsurface. Thus, we test the proposed method in this more complex situation. The black stars in Figure 4a shows the source locations, and we use forward modeling from these sources to obtain the data, recorded by a varying numbers of random receivers. With the same hyperparameters and training process mentioned in the last section, we tested the proposed method with data received at 50 (Figure 4 a), 40 (Figure 4 c), 30 (Figure 4 e), and 20 (Figure 4 g) random receivers. As shown in Figures 4b,d,f, and h, the source locations are accurately extracted by the proposed method. Even in the case of only 20 random receivers placed along the 12.5 km distance, we can still get focused source imaging results. This demonstrates the capability of the proposed method in handling sparse and irregular observations, showing good potential for the application to field data.

But before that, we need to test the performance of the proposed method in the passive seismic monitoring scenario as this is a more common case. The above experiments are built on an active seismic scenario in which the source is ignited at zero time. For passive seismic data, we often do not know the source ignition time. Thus, instead of the direct summation of the wavefields over frequencies, we inverse Fourier transform the wavefields to obtain time-domain snapshots. The snapshot that shows the largest focusing of energy is often regarded as the ignition time of the source [Artman et al., 2010]. Thus, with a similar setting to the previous experiments, we execute the forward modeling from the triple sources but ignite the source at 5 s, which we consider as an unknown. We use the proposed pipeline to obtain the time-domain snapshots, and we show several snapshots from 4.75 to 5.5 s in Figure 5. We observe that the energy focuses best at the snapshot at 5.00 s, which is consistent with the ignition time. The locations where the energy focus are also consistent with the ground truth of the source locations.

## 4 The field test to the hydraulic fracturing data

The above synthetic examples have demonstrated the potential for a field data application. Thus, in this section, we have further tested on the field vertical-component data acquired in the Arkoma Basin in North America. The data are recorded through passive seismic monitoring during a hydraulic fracturing stimulation for a shale gas reservoir [Anikiev et al., 2014, Stanek and Eisner, 2017, Wang et al., 2022a]. Figure 6a shows the whole region covering a square area of 5\(\times\)5 km\({}^{2}\) with a maximum depth of 2.45 km. There were 75 events (red dots) captured in the four days of monitoring and recorded by 911 receivers (blue dots) on the surface. The 3D P-wave velocity (Figure 6b) is obtained from an active seismic experiment in the region and adjusted by the well-log information [Anikiev et al., 2014]. In this case, we do not have the exact zero-tag information of the recordings. Here, to demonstrate the performance of the proposed method on the field data, we select two events based on their signal-to-noise ratio of the recorded event and the geometry of the line and the event, so it complies with the 2D implementation of the proposed method. The first event is extracted from

Figure 3: a) shows the Overthrust velocity model, the source event location (denoted by the black star), and the random receivers (denoted by the red triangles); b) is the source image courtesy of the proposed method.
the recordings along line 10, denoted by the black box, shown in Figure 7. The corresponding 2D velocity profile and recorded data for the event are shown in Figure 8. Prior to applying our PINN imaging, we use several pre-processing steps, including the Non-local means filter [Buades et al., 2005] and a bandpass filter on the data.

Figure 4: The source location results using various numbers of random receivers, e.g. 50 (a), 40 (c), 30 (e), and 20 (g) receivers denoted by the red triangles. The sources denoted by the black stars are distributed at different locations: \((5.0,2.5)\), \((7.5,3.0)\), and \((10.0,3.5)\) km. The source location results corresponding to various numbers of receivers by the proposed method are shown in b), d), f), and h).

Figure 5: Time-domain snapshots of wave propagation at various times, e.g., a) 5.50 s, b) 5.25 s, c) 5.00 s, and d) 4.75 s, where the energy of the wavefield focus best at 5.0 s (c).
After preprocessing, we start the proposed pipeline (Figure 2). We transform the recorded data to the frequency domain and train an NN of size {18, 128, 2} with positional encoding where L=4 to fit the 3-12Hz recorded data on the surface. As we claimed before, the data fitting network could deal with a big part of the receivers missing. Here, we randomly pick 70 receivers out of 122 for the training. The benefit of this coarse fitting is to allow for a smooth representation of the data that captures the

Figure 8: The original recorded data for a single event a) and the 2D P-wave velocity profile b), where the source location provided by the data provider is denoted by a red star.

Figure 6: Passive seismic monitoring in Oklahoma, a) is the recording geometry in this region where the blue dots denote the receivers, and red dots denote the locations of the previously estimated events projected on the 2D plane; b) is the corresponding 3-D P-wave velocity.

Figure 7: The selected event 8 and the recording line. a) is the original geometry, where the black box denotes the selected recording line, and the black star denotes the location of the event; b) is the zoom of the corresponding geometry.
key moveout information. After 4000 epochs of supervised training with an Adam optimizer and a learning rate of 1e-3, the predicted data are shown in Figure 10. We could obviously observe that the main features of the original data are well-reconstructed by the neural network, granted we ignore the sharp changes, which could be due to noise (wavefields are inherently smooth).

Then the next step is to perform the PDE fitting branch in Figure 2. Here we use a larger NN of size {27, 256, 256, 128, 128, 64, 64, 2} with positional encoding where L=4 to ensure the capacity of the NN is enough to represent this high-dimensional wavefield. The loss function here is the causality implementation of equation 5. After training with an Adam optimizer using a learning rate of 1e-3 for 1500 epochs, we obtain the frequency-domain wavefield. An inverse Fourier transform provides time-domain snapshots (source imaging). Figure 11 shows these snapshots in which the source image seems to focus on the right source location consistent with the provided location at the estimated time of around 0.7 s.

In the proposed method, we do not need label information for source location during the training as the supervision comes from the governing equations and the data. On the other hand, for supervised based approaches, errors in the labels would result in an error-prone trained model. For example, as for event 62 in line 2 (Figures 12 and 13), the source image result is shown in Figure 14, where the image focuses at a location different than the provided source locatio

Figure 10: A comparison between the data fitting result and recorded data of 3Hz for the first 100 receivers. The blue line denotes the recorded data, while the orange line is the predicted data by NN.

Figure 9: The recording of event 8, a) is the original recording with recording time from 0.8 s to 1.6 s; b) is the corresponding processed data from 3 Hz to 12 Hz after non-local means filtering where the patch size used for denoising is 15, maximal distance in pixels where to search patches used for denoising is 21, and the cut-off distance is 1.0.
supervised learning). We simulate the data from the provided source location label using the 2D P-wave velocity profile and a finite difference wave equation solver. We realize that the simulated records and recorded data are very different in curvature and shape (Figure 15). Using the simulated data, we image the source using our PINN approach and obtain the images in Figure 16. The image focused at the location of the source. This demonstrates that the provided source location is not reliable, but our method is consistent, and forms an adjoint to the forward modeling, which is based on a finite difference method.

Figure 11: The time-domain snapshots of the wavefield at various times. It seems that the wavefield focus best at 0.7 s. The red box denotes the location where the energy is focused, which is consistent with the provided location (red star).

Figure 12: The selected event 62 and the recording line 2. a) is the original geometry, where the black box denotes the selected recording line, and the black star denotes the location of the event; b) is a zoom of the corresponding geometry.
## 5 Discussion

In this paper, we showed the results of our novel direct source imaging method using PINNs with hard constraints on both synthetic and field data. The functional form of NNs offers flexibility for irregular and sparse recordings. The synthetic demos and field examples demonstrate the effectiveness and potential for source imaging. In the following subsections, we will discuss the causality implementation, the drawbacks of the method, and further potential improvements.

### The performance with/without causality implementation

As we claimed earlier, using the PDE loss with causality could accelerate the convergence and improve the final predictions. Here, we take the triple-sources case with 20 random receivers to demonstrate the benefits of the causality implementation. We use the same velocity model to generate the data, the

Figure 14: The time-domain snapshots of the wavefield at various times. The red box denotes the location where the energy is focused but is inconsistent with the given label.

Figure 13: The original shot gathers of line 2 due to event 62, a) and the 2D P-wave velocity profile b), where the red star denotes the source location given by the data provider.
same network configuration but with 3000 epochs in the training, and the same workflow to obtain the source imaging results. For the causality implementation, we set \(\epsilon_{0}=1e-7\) and \(\lambda=1e-5\). Figure 17 shows the loss curves of the training process with and without the causality modification. The convergence with causality implementation is much faster than the vanilla implementation. The prediction after 3000 epochs (Figure 18) with causality implementation is better than the one without causality implementation. The reconstructed source image (Figure 18b) has better energy focusing.

Figure 16: The time domain snapshots of the wavefield at various times with the proposed method applied to the simulated data (Figure15b).

Figure 15: a) is the observed data while b) is the simulated data on the 2D velocity profile (Figure13b) due to the event on the location given by the data provider.
### Resolution

The proposed method still has room for further developments in terms of resolution. This can be achieved by extending the frequency range used in the proposed method. High frequencies have been a thorn for PINNs. However, with the reference frequency approach and the strategy of the frequency extension (Huang and Alkhalifah, 2022; Huang et al., 2022), we can gear the approach to work better. Increasing the frequency range used in the data fitting and PDE fitting would increase the resolution of the final source image.

Figure 17: The loss curves of training with and without the causality implementation. The purple line is the version without causality implementation while the black line is with causality implementation.

Figure 18: The source image a) with a conventional loss, and b) with a causal loss equation 6.
### Efficiency

Although the workflow could deal with irregular and sparse observations and show highly reliable performance, there exists the main drawback of PINNs for now and that is efficiency. The cost is reasonably high, especially considering that the wavefield corresponds to a particular velocity and source location. However, PINNs are new compared to, for example, the finite difference methods for solving the wave equation, so speed up is a matter of time. Even with the slower implementation, the flexibility in PINN in its functional representation of the wavefield, no mesh is required, that allows it to handle irregular geometry might make it all worthwhile to explore and utilize.

## 6 Conclusion

We present a novel direct microseismic imaging framework utilizing physics-informed neural networks with hard constraints. It allows us to image the source due to the use of the neural network functional representation, in addition to its interpolation capabilities. Meanwhile, we propose a loss function with causality with respect to depth to accelerate the convergence and improve the prediction of the NN. With only 20 random receivers on a 12.5 km lateral stretch, the multiple-source event is stably and accurately located. We successfully applied to the Oklahoma Arkoma Basin hydraulic fracturing data, showing the effectiveness of the proposed method and its potential for large-scale source location.

## 7 Acknowledgments

The authors thank KAUST and the DeepWave Consortium sponsors for supporting this research. We thank Microseismic Inc. for the use of the Arkoma data, and Hanchen Wang and Fu Wang for discussing the field data preprocessing. We would also like to thank KAUST for its support and the SWAG group for the collaborative environment. This work utilized the resources of the Supercomputing Laboratory at King Abdullah University of Science and Technology (KAUST) in Thuwal, Saudi Arabia. The datasets in this paper are available at [https://doi.org/10.5281/zenodo.7543192](https://doi.org/10.5281/zenodo.7543192).
