# Micro-seismic Elastic Reflection Full Waveform Inversion with An Equivalent
Source

Hanchen Wang1
&Qiang Guo2
&Tariq Alkhalifah3
&1Los Alamos National Laboratory, US 2SLB, UK
3King Abdullah University of Science and Technology, KSA
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
###### Abstract

In micro-seismic event measurements, pinpointing the passive source’s exact
spatial and temporal location is paramount. This research advocates for the
combined use of both P- and S-wave data, captured by geophone monitoring
systems, to improve source inversion accuracy. Drawing inspiration from the
secondary source concept in Elastic Reflection Full Waveform Inversion
(ERFWI), we introduce an equivalent source term. This term combines source
functions and source images. Our optimization strategy iteratively refines the
spatial locations of the source, its temporal functions, and associated
velocities using a full waveform inversion framework. Under the premise of an
isotropic medium with consistent density, the source is defined by two spatial
and three temporal components. This offers a nuanced source representation in
contrast to the conventional seismic moment tensor. To address gradient
computation, we employ the adjoint-state method. However, we encountered
pronounced non-linearity in waveform inversion of micro-seismic events,
primarily due to the unknown source origin time, resulting in cycle skipping
challenges. To counteract this, we devised an objective function that is
decoupled from the source origin time. This function is formulated by
convolving reference traces with both observed and predicted data. Through the
concurrent inversion of the source image, source time function, and velocity
model, our method offers precise estimations of these parameters, as validated
by a synthetic 2D example based on a modified Marmousi model. This nested
inversion approach promises enhanced accuracy in determining the source image,
time function, and velocity model.

##  1 Introduction

In recent years, the field of micro-seismic events has been a hot topic in the
exploration geophysics community, due mainly to the fact that hydraulic
fracturing is becoming a widely used process in unconventional exploration in
the oil and gas reservoirs with dense rocks such as shale. Clouds of micro-
seismic events indicate the growth of faults and the occurrence of fractures,
which can help the engineers reactivate failed regions or discover new ones.
The hydraulic stimulation will often trigger or induce micro-seismic events
[15]. A high pressure condition is created in the reservoir area by hydraulic
injections, eventually cracking the dense rock in order for the gas and oil to
flow freely. Such cracks can be monitored by seismic sensors, either in a well
or on the earth’s surface, detected by passive seismic event detection
techniques [10]. Therefore, petroleum engineers are looking for more accurate
methods to locate these micro-seismic events with the recorded signals, as
knowing the accurate locations of micro-seismic events is critical for
engineering the hydraulic process [49].
The most common method of finding micro-seismic source locations or passive
earthquakes is based on traveltime picks of both acoustic and shear wave
arrivals [16, 45]. Such techniques require the manual determination of micro-
seismic event arrival times in the recorded data. With the help of a simple
relationship between traveltime, velocity and distance, each event can be
located by smearing traveltimes over isochrones; ideally, several isochrones
of different traces will intersect at the same point, granted the accurate
velocity, referred to as the source location. This manual picking process,
however, can hardly be done accurately, due to the obvious reason that micro-
seismic data naturally have a poor signal-to-noise ratio (S/N). [8, 26, 37].
Alternatively, one can use the migration based methods [14, 27, 28, 29, 38, 4,
21, 31, 46], among which using the time-reversal imaging (TRI) to migrate the
micro-seismic signals to the source locations is a popular one. With TRI, less
user input is required compared to the traveltime based methods. It also
provides more accuracy than the traveltime based methods by utilizing not only
the traveltime, but also the waveform information. We can back propagate the
time-reversed recorded traces, guided by the wave equation, and focus the
signals at the sources, in both time and space. Not having to pick arrival
times is the biggest advantage of TRI methods, especially when imaging with
low S/N seismic data. If a good enough acquisition system is available,
provided we have an accurate model and physics representation, the correct
source locations will be imaged. If, however, only a limited aperture or a
sparse spatial sampled receiver system is available, TRI is obviously not the
optimal method to locate the passive sources. Much research has been published
on handling poor acquisition systems and sparsely sampled wavefields [17, 43,
44, 5, 23, 34, 19]. All TRI methods to locate the micro-seismic source may
suffer from spatially defocused images, which is caused by the errors between
the true and the imaging velocity models. Migration velocity analysis (MVA)
may be a solution to minimize the velocity error [39, 35].
Another approach of micro-seismic imaging is based on the full waveform
inversion (FWI) technique [24, 6, 41, 25, 36]. This technique makes it
possible to generally update the micro-seismic event locations and the
velocity models by minimizing the objective function that is usually described
by the differences between the true and predicted data as any other FWI. The
update is terminated when the data misfit is sufficiently small. The inverted
source location is considered as the final located source position. Most
existing FWI methods, however, assume an acoustic case and sufficiently good
initial conditions are also essential to overcome the cycle skipping problem
caused by the unknown source origin time, which makes it challenging in
practice.
On the other hand, recent development of machine learning methods shows great
potential of efficient and accurate micro-seismic event relocation and
velocity model estimation [53, 54, 55, 56]. The machine learning methods take
seismic observations as input images and predict corresponding seismic
parameters as outputs, such as velocity models and passive seismic event
locations. However, such data-driven methods require high demanding labeled or
unlabelled datasets for the network training stage, and without enough high
quality data, these methods will not perform at the desired level.
Besides the source location and function, the source mechanism inversion is
also critical for passive seismic data. A frequency domain method is proposed
to localize micro-seismic events blind to the source time function [22].
However, they confine their search to the hypo-centers of the sources and the
normalized amplitudes of the moment tensors. Other methods were introduced
that include simultaneous determination of the micro-seismic event location
and the source mechanism given the input velocity model [13, 51, 3]. The key
requirement for these methods is the knowledge of the correct velocity model.
When the velocity is inaccurate, the source location and mechanism are
affected.
The work presented here builds on the work of [47], based on an acoustic FWI
method, and further introduces an equivalent source term as the force function
in the elastic wave equation. This term is extended from the concept of the
equivalent sources in elastic reflection waveform inversion [20]. The
introduced source term is constructed by separating the source term in space
and time, named as the source images and the source functions, respectively.
An FWI approach measuring the L2-norm data misfit is used to update the
spatial and the temporal components of the source. An FWI scheme with a source
time function independent objective function is used to mitigate the cycle
skipping problem caused by the unknown source origin time while updating the
velocity models. This source function independent FWI approach works well with
any given source function, because it is based on convolved wavefields between
the observed and modeled ones to mitigate the role played by the source
function [11]. After all the parameters are inverted, the proposed method can
accurately estimate the source with initial velocity models, which are quite
different from the true ones so that the key requirement of accurate migration
velocity can be relaxed. We implement the proposed method to a modified
Marmousi model to show the novelty of inverting for the source and the
velocity models with a nested approach.

##  2 Equivalent Source In Elastic Micro-seismic Source Case

In this section, we introduce an equivalent source term for the elastic micro-
seismic case. We borrow the concept of the scatters (or perturbations) in
Elastic Reflection Full Waveform Inversion (ERFWI) to represent the micro-
seismic source term in our inversion method. Thus, we first review the concept
of the ERFWI and then we derive the equivalent source term and its gradients.

###  2.1 Scatters In Elastic Reflection Waveform Inversion

Full waveform inversion is highly non-linear, and the non-linearity increases
in the elastic cases because of the complicated wave propagation in the
subsurface medium. ERFWI is one way to overcome this issue, by separating the
model parameters into low-wavenumber components, which controls the
kinematics, and the high-wavenumber scattering [48, 2].
In an isotropic medium with a constant density, it is assumed that the squared
$V_{p}$ and $V_{s}$ velocities (denoted by
$\alpha^{2}\,\mbox{and}\,\beta^{2}$, respectively) can be divided into two
parts,

$\alpha^{2}=\alpha_{0}^{2}+\delta\alpha^{2}\,\mbox{and}\,\beta^{2}=\beta_{0}^{2}+\delta\beta^{2}$
(1)

where $\alpha_{0}$ and $\beta_{0}$ refer to the low-wavenumber (background)
parts; $\delta\alpha$ and $\delta\beta$ are the high-wavenumber scattering of
the P- and S- velocities, respectively. The total wavefield, denoted by $\bf
u$, can thus be separated correspondingly into two parts,

$\bf{u}=\bf{u}_{0}+\delta\bf{u},$ (2)

where $\bf u_{0}$ is the wavefield propagating in the smooth background
medium, and $\delta\bf u$ is the scattered wavefield simulated by the
scatters. The perturbations in the stress tensors ($\delta\tau$) can be
derived from Hooke’s law, giving us:

$\begin{array}[]{cl}&\displaystyle\delta\tau_{xx}=\left(\alpha^{2}_{0}+\delta\alpha^{2}\right)\frac{\partial\delta
u_{x}}{\partial
x}+\left(\alpha^{2}_{0}+\delta\alpha^{2}-2\beta^{2}_{0}-2\delta\beta^{2}\right)\frac{\partial\delta
u_{z}}{\partial z}+\widetilde{\tau}_{xx};\\\
&\displaystyle\delta\tau_{zz}=\left(\alpha^{2}_{0}+\delta\alpha^{2}-2\beta^{2}_{0}-2\delta\beta^{2}\right)\frac{\partial\delta
u_{x}}{\partial
x}+\left(\alpha^{2}_{0}+\delta\alpha^{2}\right)\frac{\partial\delta
u_{z}}{\partial z}+\widetilde{\tau}_{zz};\\\
&\displaystyle\delta\tau_{xz}=\left(\beta^{2}_{0}+\delta\beta^{2}\right)\left(\frac{\partial\delta
u_{x}}{\partial z}+\frac{\partial\delta u_{z}}{\partial
x}\right)+\widetilde{\tau}_{xz},\end{array}$ (3)

where $\widetilde{\tau}_{ij}$ can be described as

$\begin{array}[]{cl}&\displaystyle\widetilde{\tau}_{xx}=\delta\alpha^{2}\left(\frac{\partial
u_{x}}{\partial x}+\frac{\partial u_{z}}{\partial
z}\right)-2\delta\beta^{2}\frac{\partial u_{z}}{\partial z};\\\
&\displaystyle\widetilde{\tau}_{zz}=\delta\alpha^{2}\left(\frac{\partial
u_{x}}{\partial x}+\frac{\partial u_{z}}{\partial
z}\right)-2\delta\beta^{2}\frac{\partial u_{x}}{\partial x};\\\
&\displaystyle\widetilde{\tau}_{xz}=\delta\beta^{2}\left(\frac{\partial
u_{x}}{\partial z}+\frac{\partial u_{z}}{\partial x}\right).\end{array}$ (4)

The stress functions of equation 4 are substituted into equation 3 as forces
in order to generate the scattered waves at the model perturbation locations.
It is considered as a de-migration procedure. It is critical to estimate the
scatter image and the incident wavefield at the scatter points as source
wavelets in the de-migration procedure [50, 20].
Displacement field in the background medium ($\alpha_{0}$ and $\beta_{0}$) are
then taken into account for estimating the model perturbations ($\delta\alpha$
and $\delta\beta$).

###  2.2 Equivalent Source Term In Micro-seismic Case And Its Gradient

We extend the concept of the scatters in ERFWI to the micro-seismic case. The
passive seismic events triggered from the deep subsurface can be considered as
separated single scatters in the de-migration procedure. Now we re-write the
micro-seismic source term by a form of the equivalent source, which produces
the same micro-seismic seismic events.
If we consider the 2D case, the elastic wave propagation is granted by the
following wave equation:

$\rho\left(\frac{{\partial}^{2}{u}_{i}}{{\partial}t^{2}}\right)-\frac{\partial}{\partial{x}_{j}}\left({c}_{ijkl}\frac{\partial{u}_{k}}{\partial{x}_{l}}\right)=-M_{ij}\frac{\partial[\delta(\textbf{x}-\textbf{x}^{s})]}{\partial
x_{j}}S(t),$ (5)

where ${\bf u}$ is the elastic displacement field with two components,
$c_{ijkl}$ are components of the stiffness tensor; $\rho$ is the density; the
r.h.s. source term is an equivalent force term [1, 12]; $\textbf{x}^{s}$ is
the source location coordinate and $S(t)$ is the source time function.
$M_{ij}$ represents the components of the full seismic moment tensor. One of
the limitations of this description for inversion applications of the source
is that the moment tensor $\bf{M}$, the source location $x^{s}$, and the
source function are separate variables, with accurate estimation of these
variables requiring the estimation of the source location.
We now modify equation 5 using an equivalent source term to mimic the source
mechanism in the micro-seismic case, which is derived from the ERFWI, and such
that the r.h.s. of the elastic wave equation becomes:

${f}_{i}=\partial_{j}\tau^{0}_{ij},$ (6)

where each component of ${\tau}_{ij}^{0}$ can be written as:

$\begin{array}[]{cl}&\displaystyle{\tau}_{11}^{0}=\delta\alpha^{2}\left({w}_{11}+{w}_{22}\right)-2\delta\beta^{2}{w}_{22};\\\
&\displaystyle{\tau}_{22}^{0}=\delta\alpha^{2}\left({w}_{11}+{w}_{22}\right)-2\delta\beta^{2}{w}_{11};\\\
&\displaystyle{\tau}_{12}^{0}={\tau}_{21}^{0}=\delta\beta^{2}\left({w}_{12}+{w}_{21}\right).\end{array}$
(7)

Here, $\delta\alpha$ and $\delta\beta$ are denoted as source images and
$\left[w_{ij}\right]$ is the multi-components source time function. The stress
tensors in equation 7 and equation 4 have the same form, but the incident
wavefield at the scatter points replaced by $\left[w_{ij}\right]$, that
$\frac{\partial u_{x}}{\partial x}$ is replaced by $w_{11}$, $\frac{\partial
u_{z}}{\partial z}$ is replaced by $w_{22}$, $\frac{\partial u_{x}}{\partial
z}$ and $\frac{\partial u_{z}}{\partial x}$ are replaced by $w_{12}$ and
$w_{21}$, respectively. In ERFWI, $\delta\alpha$ and $\delta\beta$ represent
the model perturbations. It is expected that all the possible scatters in the
model are described by the two parameters. The back-propagated wavefield $\bf
u$ is saved at each model point as the source term at the scatter point. In
the micro-seismic case, $\delta\alpha$ and $\delta\beta$ represent the actual
micro-seismic sources in space, expressed in terms of perturbations in $P-$
and $S-$ waves velocities, respectively. Compared to the model perturbations
in the ERFWI case, instead of having scatters all over the model, it is
expected to have only one single scatter at the source location. Thus, the
back-propagated wavefield is only needed at the source location as its source
time function. In this case, the source will also have inherent representation
of the radiation patterns besides the P- and S- wave velocity perturbations,
which we will further discuss in the following sections. The source time
function $\left[w_{ij}\right]$ corresponds to the perturbation locations.
Later we will estimate the coefficients $\left[w_{ij}\right]$ by the back-
propagated wavefield at the source location. Thus, the source spatial location
$\bf x^{s}$ and the seismic moment tensor M are embedded into these source
images. This decomposition of the source assumes that the time and space
dependency are decoupled, which leads to a more manageable space of unknown
that is spatially equal to the model space.
If we ignore the density variable, by considering it as constant for example,
the wave equation then becomes:

$\left(\frac{{\partial}^{2}{u}_{i}}{{\partial}t^{2}}\right)-\frac{\partial}{\partial{x}_{j}}\left({a}_{ijkl}\frac{\partial{u}_{k}}{\partial{x}_{l}}\right)=\nabla\cdot{\tau}_{ij}^{0}.$
(8)

where $a_{ijkl}$ is the density normalized elastic coefficients. In reality
the density may have some influences on the passive seismic data inversion,
but density is not easy to recover. To simplify the inversion, we here choose
to ignore the density in this study. Considering the data corresponding to the
passive events being acquired on the Earth surface (or in a well),
$\textbf{d}(x_{r},t)$, we can devise an iterative approach by defining the
optimization problem using the classic least-square minimization of

${E}\left(m\right)=\frac{1}{2}\sum_{r}{\int{\left\|\textbf{u}\left({x}_{r},t;m\right)-\textbf{d}\left({x}_{r},t\right)\right\|^{2}_{2}dt}},$
(9)

with respect to $\delta\alpha$, $\delta\beta$ and $\left[w_{ij}\right]$, where
$\textbf{u}(x_{r},t;m)$ is the modeled waveform recorded by the receivers at
$x_{r}$ determined by the current parameters, including $\delta\alpha$,
$\delta\beta$, $w_{ij}$ and velocities. The parameters $\delta\alpha$ and
$\delta\beta$ can be retrieved iteratively using the following gradients:

$\nabla_{\delta\alpha}{E}=-{\int_{t}{dt\left({w}_{11}+{w}_{22}\right)\left(\frac{\partial\widehat{{u}}_{x}}{\partial
x}+\frac{\partial\widehat{{u}}_{z}}{\partial z}\right)}},$ (10)
$\nabla_{\delta\beta}{E}=-\int_{t}dt\,2{w}_{11}\frac{\partial\widehat{{u}}_{x}}{\partial
x}+2{w}_{22}\frac{\partial\widehat{{u}}_{z}}{\partial
z}+({w}_{12}+{w}_{21})(\frac{\partial\widehat{{u}}_{x}}{\partial
z}+\frac{\partial\widehat{{u}}_{z}}{\partial x}),$ (11)

where $\hat{\bf u}$ refers to the adjoint wavefield. Besides the source
spatial location, $\delta\alpha$ and $\delta\beta$ represent the source
mechanism in a more compact form other than the seismic moment tensor.
Although the objective function we use to update the velocity models, which
will be introduced in the following section, frees us from knowing the source
time function prior to the inversion, we prefer to apply a simultaneous
inversion for the source time functions, in order to get a more accurate
estimate than the case of ignoring the time functions:

$\begin{array}[]{cl}&\displaystyle\nabla{E}_{{w}_{11}}=\int_{x}\int_{z}\,dxdz\delta\alpha^{2}\left(\frac{\partial\widehat{{u}}_{x}}{\partial
x}+\frac{\partial\widehat{{u}}_{z}}{\partial
z}\right)-2\delta\beta^{2}\frac{\partial\widehat{{u}}_{z}}{\partial z};\\\
&\displaystyle\nabla{E}_{{w}_{22}}=\int_{x}\int_{z}dxdz\,\delta\alpha^{2}\left(\frac{\partial\widehat{{u}}_{x}}{\partial
x}+\frac{\partial\widehat{{u}}_{z}}{\partial
z}\right)-2\delta\beta^{2}\frac{\partial\widehat{{u}}_{x}}{\partial x};\\\
&\displaystyle\nabla{E}_{{w}_{12/21}}=\int_{x}\int_{z}dxdz\,\delta\beta^{2}\left(\frac{\partial\widehat{{u}}_{x}}{\partial
z}+\frac{\partial\widehat{{u}}_{z}}{\partial x}\right).\end{array}$ (12)

More details on the gradient derivations are provided in Appendix $A$.
In multi-parameter optimization, there will always be the potential for trade-
off. In this case, we need to wisely choose the step-length of each parameter
or use a small enough fixed step-length to mitigate the coupling issue.

##  3 Source Image Radiation Patterns Analysis

In this section, we show the radiation pattern analysis of the source images,
$\delta\alpha$ and $\delta\beta$, as well as the relationship between the
source images and the seismic moment tensors, so that we can better understand
the physical meaning of the source image term we introduced.
The most complete form used to describe the source mechanism is the seismic
moment tensor. In 3D, nine force couples are used to fully represent every
possible source [1]. The moment tensor in 2D can be written in the following
symmetric matrix form

$\bf{M}={M_{0}}\left(\begin{matrix}M_{11}&M_{12}\\\
M_{21}&M_{22}\end{matrix}\right),$ (13)

where $M_{12}$ equals to $M_{21}$, as the system would otherwise be unstable
and would rotate. $M_{0}$ refers to the moment magnitude. Thus we have three
independent components in 2D to mathematically represent a seismic source
mechanism.
In our proposed theory, $\delta\alpha$ and $\delta\beta$ can also be used to
describe seismic source mechanisms because they are simply another formulation
of equation 5. Both equivalent source terms should have the same physical
meaning. Since we merge the moment tensor and the space coordinates into a
single term, the source images, we need to find out what exactly the source
images represent. A good way to understand this representation is to find the
connection between the source images and the well-known moment tensor. In
order to understand the relationship between the two, a simple experiment is
performed under ideal conditions. A source is ignited in the center of a 2D
square homogeneous velocity model, with $V_{p}=2km/s$ and
$V_{s}=\frac{2}{\sqrt{3}}km/s$. Every grid point on all four boundaries is set
as a receiver, in order for the full Green’s function to be measured. We use a
moment tensor to describe the source in forward modeling, and then the event
is recorded by the receivers. We subsequently perform an inversion of the
source images using equations 10 and 11, under the assumption that the
original velocities and the original source time function are known. We can
then compare the inverted results of source images with the moment tensor. The
inverted source images are shown in Figure 1. There are other ways of
representing the moment tensor with images. For example, a diffraction stack
migration combines ray-theory Green’s functions with a reverse time imaging
method to simultaneous locate and characterize the seismic source properties
[9]. Although the idea seems quite similar to ours, the theory part is quite
different. In this paper, we use our inverted source images to represent the
source mechanism.
These inverted source images show different features, indicating different
source types. The cases of $M_{1}$ and $M_{2}$ show an explosive and an
implosive sources, which means the moment tensor is either a positive or a
negative diagonal matrix. Both $\delta\alpha$ and $\delta\beta$, Figures 1(a)
to 1(d), show round dots at the source location, which means that the source
is symmetric in all directions in space. In other words, Figures 1(a) to 1(d)
show the cases of isotropic sources, just like what we used for forward
modeling. For cases of $M_{3}$ and $M_{4}$ (Figures 1(e) to 1(h)),
$\delta\alpha$ values equal to zero because according to equation 10,
$(w_{11}+w_{22})$ equals to $0$, and therefore, the whole gradient is zero. In
the source images of double coupled (DC) and compensated-linear-vector-dipole
(CLVD) sources $(M_{3}$ to $M_{9})$, positive (white color) and negative
(black color) values can be directly connected to the source type by looking
at the sign distribution at the source location. This is because positive
values mean that the wave is propagating outward, and negative values mean
that the wave is propagating inward. On the other hand, the source signatures
in these source images exhibit the same features with the beach ball plots of
the same moment tensors. These observations lead us to the conclusion that
source images in Figure 1 represents components of a moment tensor.
For $M_{1}$ and $M_{2}$ (Figures 1(a) to 1(d)), source images exhibit the same
features, no matter whether the source is an implosive or explosive one, which
means that we cannot tell the difference between the expansive and compressive
sources. The same issue also shows in Figure 1(m) to 1(p). Both cases have the
similar features in $\delta\alpha$ and $\delta\beta$, but the moment tensors
for the two sources have opposite signs. This is reasonable as the moment
tensor in 2D case has three independent parameters $(M_{11},M_{22}$ and
$M_{12})$, whereas here, we only have two parameters
$(\delta\alpha,\delta\beta)$. As a result, some information is lost.
In field experiments, we rarely get a perfectly illuminated acquisition
system, unlike in the experimental setup described above. Therefore, here, we
present a more realistic set of tests with surface receivers only. The other
conditions, including the moment tensors, remain the same as in the previous
tests.
In this set of tests, the source images share the same features with the
source signatures and in almost all cases, the resolution along the horizontal
direction is poor and makes it difficult to recognize what source types each
case refers to. This is expected, since we have only surface receivers. We
lose most of the source signature information from left, right and bottom.
This issue could be even more severe in 3D cases with 2D plane receivers,
because we should lose more information and six independent moment tensor
components need to be represented by two source images. There is a limitation
in identifying source types, when using these source images due to the fact
they highly rely on with acquisition systems (illumination).
To understand the sensitivity of source images with respect to illumination, a
radiation pattern analysis in $2D$ case is presented here, with radiation
pattern plots shown in Figure 3.
With these radiation pattern analysis results, we may recognize whether the
source image difference is identifiable or not under certain acquisition
configuration. For example, in cases $3$ and $4$, both P- and S-waves show
different sensitivities from angle $-45(315)$ to $45$ degrees, which reflect
the presence of surface receivers only. Relatively large differences in source
images are observed, which allow us to identify the different source types
(Figures 2(e) to 2(h)). However, we cannot detect any difference between cases
5 and 9 (surface receivers only) or between cases 6, 7 and 8 (vertical well
receivers only), since they share similar sensitivities at certain angles.
Although the source images indicate the source types, we still can hardly rely
on these parameters to identify which types they are, as the source images
highly depend on the acquisition aperture and inversion quality of the other
model parameters than source images, which is the trade-off commonly
experienced in FWI. With these tests and analysis, we now understand the
physical meaning of the source images. The most reliable information provided
by the source images, which is also the most important one, is the source
location in space. Besides the spatial locations, $\delta\alpha$ and
$\delta\beta$ also represent the source mechanism to a certain extent (when
the model parameters are good enough and the acquisition is dense). If
additional information is needed, such as the moment tensor, one can first
locate the source by the proposed method according to the inverted source
images, and then invert the moment tensor with a conventional seismology
method [13].

##  4 Velocity Update With Source Function Independent Objective Function

The velocity models can be updated by minimizing the data misfit in equation
9. The gradients with respect to $\alpha$ and $\beta$ can be written as:

$\nabla E_{\alpha}=-\int_{t}dt(\frac{\partial u_{x}}{\partial
x}+\frac{\partial u_{z}}{\partial z})(\frac{\partial\hat{u}_{x}}{\partial
x}+\frac{\partial\hat{u}_{z}}{\partial z});$ (14) $\nabla
E_{\beta}=-\int_{t}dt(2\frac{\partial u_{x}}{\partial
x}\frac{\partial\hat{u}_{x}}{\partial x})+(2\frac{\partial u_{z}}{\partial
z}\frac{\partial\hat{u}_{z}}{\partial z})+(\frac{\partial u_{x}}{\partial
z}+\frac{\partial u_{z}}{\partial x})(\frac{\partial\hat{u}_{x}}{\partial
z}+\frac{\partial\hat{u}_{z}}{\partial x}).$ (15)

In case of a poor initial guess of the velocity model, the initial source
terms will be far from the true ones in both space and time, especially the
source origin time. Using such an inaccurate source function as the initial
input in the inversion is not acceptable because a slight time shift in the
source origin time leads to severe cycle skipping problems in the data. To
avoid such cycle skipping problems, we can use some advanced objective
functions such as the adaptive waveform inversion introduced by [40]. We
propose to use a convolution based objective function [11, 47], which has
source time function independency, as we mitigate the role of the source
origin time in the objective function. We also add a model total variation
(TV) term as a regularization in the objective function. Thus, the whole
objective function can be written as:

${E_{indp}}\left(m\right)=\frac{1}{2}\sum_{r}{\int{\left\|\textbf{u}\ast{\textbf{d}}_{ref}-\textbf{d}\ast{\textbf{u}}_{ref}\right\|^{2}_{2}dt}}+\lambda\sqrt{\frac{\partial
m}{\partial x}^{2}+\frac{\partial m}{\partial z}^{2}+\epsilon^{2}},$ (16)

where the subscript ${\ast}_{ref}$ represents the reference traces in both the
true and the predicted data. $\lambda$ is a factor that balances the weight of
the regularization term. $\epsilon$ is a factor that controls the smoothness
of the gradient. The larger $\epsilon$ is, the smoother the gradient will be.
The symbol $*$ is the convolution operator. In practice, our experimentation
suggests that a near-offset trace is a better reference than a far-offset
trace. The TV term makes the gradients blocky. However, in our experiment, the
main goal of adding this TV regularization is to obtain smooth velocity
updates in every iteration with a large $\epsilon$. In equation 16, a new
source function that contains both the predicted and the true data source
function information is convolved in both terms equally, which means that the
effect of the delay in the source origin time is muted. Thus, the true source
origin time is no longer necessary. The velocity gradients of the objective
function are calculated by minimizing equation 16 with respect to $\alpha$ and
$\beta$. The gradients can then be given as (see Appendix B):

$\frac{\partial E_{indp}}{\partial
V}=\sum_{i}^{nr}{\left[\left(\frac{\partial{u_{i}}}{\partial
V}*{\textbf{d}}_{ref}\right)\cdot\textbf{r}_{i}-\left({\textbf{d}_{i}}*\frac{\partial{\textbf{u}}_{ref}}{\partial
V}\right)\cdot\textbf{r}_{i}\right]}+\lambda\nabla\cdot\left(\frac{\nabla
V}{\sqrt{\frac{\partial V}{\partial x}^{2}+\frac{\partial V}{\partial
z}^{2}+\epsilon^{2}}}\right).$ (17)

The adjoint wavefields for the gradient calculation are thus simulated by
injecting the adjoint sources:

$\textbf{r}^{(1)}_{i,j}=\textbf{d}_{i,ref}\bigotimes\left({\textbf{u}}_{i,j}\left({x}_{r},t\right)*{\textbf{d}}_{i,ref}-{\textbf{d}}_{i,j}*{\textbf{u}}_{i,ref}\left({x}_{ref},t\right)\right),$
(18)
$\textbf{r}^{(2)}_{i,j}=-\textbf{d}_{i,j}\bigotimes\left({\textbf{u}}_{i,j}\left({x}_{r},t\right)*{\textbf{d}}_{i,ref}-{\textbf{d}}_{i,j}*{\textbf{u}}_{i,ref}\left({x}_{ref},t\right)\right),$
(19)

at the $i^{th}$ receiver locations and the reference trace location,
respectively, where $\bigotimes$ refers to the cross-correlation operation in
time. The two seismograms in equations 18 and 19 are back-propagated at
meanwhile to keep the computational cost equivalent to FWI.
A nested approach is used as the inversion strategy, i.e. $(1)$ updating the
source images with equations 10 and 11 (for five iterations and this was
decided by experimentation); $(2)$ updating the velocities using equation 17
with the adjoint sources equations 18 and 19 (for ten iterations); $(3)$
updating w using equation 12 (for five iterations). We repeat the above steps
until convergence is achieved, which, in our case, means the objective
function is no longer decreasing. Details are discussed in the implementation
section below with the workflow.

##  5 Numerical Examples

In this section, we establish our work on a modified elastic Marmousi model
whose S-wave velocity is modified by adding a half-circle-shaped low velocity
layer (representing the hydraulic injected area) pointed out by an arrow in
Figures 4(a) and 4(b). Five sources are separately triggered at locations
$x(i)=2.4+0.2(i-1)$ km, $i=1,...,5,z=1.55$ km.
The model size is $6$km by $2$km with a $0.01$km spatial interval in the $X$
and $Z$ directions, respectively. Receivers are set on all grid points on the
top surface and in a vertical well, which locates from $(x=1,z=0.5)$ km to
$(x=1,z=1.5)$ km in space with the same $0.01$km receiver interval. Totally
$600$ receivers on the top surface and $100$ receivers in the well work as our
monitoring system. We use a $10$Hz dominant frequency Ricker wavelet as the
true source function. In this case, five grid points represent the minimum
wavelength. The record has been convolved with a high-pass filter with the
minimum frequency of $3$Hz to make the experiment more realistic. The initial
$\alpha$ and $\beta$ are generated by applying strong smoothing to the true
velocities and then making $\alpha$ and $\beta$ $10\%$ slower. As a result,
shown in Figures 4(c) and 4(d), the initial models are almost linearly
increasing with depth. The smoothing operator here is a Gaussian smoothing
window with a length of 300m in the X and 100m in the Z directions and we
repeat the smoothing for 3 times. The $10\%$ deviation makes the starting
model here even harder to invert than the linear increasing velocity model,
which is widely used for Marmousi inversion experiments. We intended to show
the effectiveness of our method using this initial model. When we generate the
initial $V_{s}$ model, we do not include the modified hydraulic injected area
information in the smoothing because it is our target zone, which is the low
velocity zone pointed by an arrow in Figure 4(b). The record length is 5
seconds with a $1$ms temporal sampling interval. The top surface is considered
as a free surface boundary and the other three boundaries have absorbing
boundary conditions applied. The modeling is implemented using a finite
difference algorithm.
It is assumed that the five sources are temporally separated, which allows us
to separate the event from each other. This is commonly satisfied in the field
monitoring. Each source has a unique moment tensor and origin time. Moment
tensors $M_{1}$ to $M_{5}$ are set to be the source types of the five sources,
where the first source on the left uses $M_{3}$, the second source uses
$M_{4}$, the third source uses $M_{1}$, the fourth source uses $M_{2}$ and the
last source uses $M_{5}$.

###  5.1 Implementation

We implement our method on this numerical experiment using a nested approach,
as shown in Figure 6. We need an initial velocity model and source term to
start the inversion. The initial velocities have been discussed already. The
initial source images are completely zero $(\delta\alpha=\delta\beta=0)$,
which means that the source spatial distribution is not available. To obtain
the initial source function, the true record is back-propagated using the
initial velocity model and then the wavefield at the maximum energy point is
extracted from the wavefield, which is similar to the time-reversal imaging
method and the maximum energy imaging condition. Once the initial parameters
are ready, we can then start the inversion and update the three parameters
iteratively until we get a convergent result. To suppress the source
signatures at the five passive source locations, we use the pseudo Hessian to
condition the gradient, which is calculated from the squared values of the
adjoint wavefields to balance the energy contribution. To mitigate the cross-
talk between different parameters, we need to carefully choose the
optimization scheme.
The data residual (Equation 9) reaches its minimum after ten outer-loop
iterations. In each outer loop, five inner loops of source image and source
time function updates and ten inner loops of velocity model updates are
executed. The weighting factor $\lambda$ corresponding to the TV term is set
to $0.5$ in the first iteration and linearly decreases to $0.05$ in the last
iteration as a function of the outer-loop iteration number. Because of the
unreliable initial velocities, we start with higher weights for the
regularization in the beginning. As the velocities improve, we use smaller
weights for the TV term to allow the velocity to have higher resolution
information. We use a constant five percent of the current iteration TV
maximum value as the smoothing factor $\epsilon$. Our objective is to have a
smooth velocity model at the end by using a small $\epsilon$, as we are
focusing on transmission waves. Usually $\epsilon$ is set to be $1\%$ to keep
the model update blocky and approximately $5\%$ to keep it smooth. If the
smoothing parameter is too large, the model update will lose most of the
details, which may have negative effect on the inversion. The weighting factor
$\lambda$ is dependent on to what extent you believe in the initial model.
Since our initial model in the experiment is quite far from the original one,
we need a relatively larger $\lambda$ in the first few iterations. The
velocity updating step lengths are fixed for the whole inversion, that is the
step length for $\alpha$ is $0.05$ and the step length for $\beta$ is $0.03$.
The step length for $\beta$ is approximately $\frac{1}{\sqrt{3}}$ of that of
that for $\alpha$. If the computational resources are sufficient, an even
smaller step-length would be better. However from our experimentation, keep
the maximum update no more than 50 meters per iteration is a good enough step-
length for the modified Marmousi model.

###  5.2 Results and Discussions

The final results of $\alpha$ and $\beta$ are shown in Figures 4(e) and 4(f).
The inverted velocity models generally exhibit low resolution, due to the
frequency band of the data and the limited illumination with only five nearby
sources. Due to the distribution of the micro-seismic events and the
receivers, transmission waves are the dominant recorded energy. Under the FWI
scheme, transmission waves mainly update the low wave-number components of the
model, whereas the reflection waves update the high wave-number parts. Thus,
the inverted velocities in the experiment are mostly controlled by the smooth
background updates. We may also find that the lower part of the model seems to
have a higher resolution than the upper part. This is expected since the parts
of the model below the sources are single scattered waves, whereas the
reflection energy above the source can only be recorded after multi-
scattering. In other words, the recorded reflection energy from the bottom
reflections is higher than that from the upper reflections, which results in
the upper part of the model having fewer high wavenumber updates
(reflections), compared to the lower part, resulting in higher resolution in
the lower part and a poorer resolution in the upper part. If the passive
sources are triggered in a wider range, the illumination of the area would be
much better with the same acquisition system. The passive sources in this
experiment is quite confined since there are only five sources in a very small
area. It leads to the accurate but smooth results for the kinematic components
of the model. We compare the vertical profiles of $\alpha$ and $\beta$ at the
position of $x=3.0$km as a quality control, shown in Figure 5. Apparently, the
initial velocities are lower than the true ones. Moreover, in the initial
$\beta$, no hydraulic injected area information is given at the arrow-pointed
zone in Figure 4(b). After the inversion, we see that the velocities have been
updated in the right direction, and the hydraulic injected area (low velocity
zone) in $\beta$ is also well retrieved. Certainly, the inverted velocities
are more kinematic accurate.
Although the inverted velocity models are imperfect that the high wave-number
part of the model is not well retrieved, they are good enough to locate the
micro-seismic sources in space, which reflects their kinematic accuracy (the
low wave-number part is accurate). The inverted source images, $\delta\alpha$
and $\delta\beta$, are shown in Figures 7(a) to 7(j) for the five sources,
respectively. Extracting the source coordinates from the source images by
searching for the maximum energy point in each source image, where most of the
energy is focused at that point, although some unfocused energy caused
probably by the remaining data misfit is left over. For verification, the
source locations after the first inner loop of the source image update, which
is equivalent to the TRI locations, are shown as the blue stars in the figures
and the true source locations are shown as the red stars. We may find that
both $\delta\alpha$ and $\delta\beta$ focus in the accurate locations, while
the initial (TRI) locations, are totally wrong due to the initial velocity
errors. Note that $\delta\alpha$ of sources 1 and 2 are zero, because of the
reason we discussed in the previous section $(w_{11}+w_{22}=0)$. In the case
of an imperfect acquisition and velocity model, the term $(w_{11}+w_{22})$ is
not exactly zero, but they are at least sufficiently small. When dealing with
the source images, we intended to ignore small values in $\delta\alpha$ as
compared to the value of $\delta\beta$. By comparing $\delta\alpha$ with
$\delta\beta$, we may find that $\delta\alpha$ exhibits lower resolution than
$\delta\beta$. It is reasonable because $\delta\alpha$ simply depends on
$V_{p}$, but $\delta\beta$ depends on both $V_{p}$ and $V_{s}$. The fact that
$V_{p}$ is naturally higher than $V_{s}$ leads to the wavelength (higher
resolution) of $\delta\alpha$ being larger (lower) than that of $\delta\beta$.
Although the source locations are accurately inverted, we can hardly tell the
differences between the source images for the various source types. It is
mostly because of the limited acquisition aperture (surface and vertical well
receivers) and the inverted velocity imperfections. In addition, there is a
leakage of energy in the inverted source images due to the model misfit as
well. One can choose to ignore these weak artifacts or use a regularization,
for example L1-norm regularization, to suppress the unfocused energy.
The inverted source wavelets are shown in Figure 8. For the first and the
second sources on the left (moment tensors $M_{3}$ and $M_{4}$), the
components $w_{11}$ and $w_{22}$ have opposite signs dictated by the source
mechanism. As a result, the summation of $w_{11}$ and $w_{22}$ admits
reasonably small values compared to the value of $w_{12}$, as we discussed
above. In theory, the $w_{11}$ and $w_{22}$ components of $M_{3}$ should be
both zero but the inversion seems to fail in inverting their correct values,
though the result does not affect the source image inversion. For the third
and the fourth sources on the left (moment tensors $M_{1}$ and $M_{2}$), the
components $w_{11}$ and $w_{22}$ show similar features due to the isotropic
source mechanism.
The computational cost for each inner loop of the proposed method is
equivalent to a standard FWI iteration. In our inversion scheme, however,
besides the nonlinear velocity update, a minor extra cost of the linear source
term update in every iteration is required. Thus, the total computational cost
is still challenging. When processing with the field monitored data,
especially with low S/N data, we could apply some pre-processing of the data
to enhance the S/N so that we may obtain clearer subsurface features.

##  6 Conclusions

We introduced an equivalent source term for micro-seismic events, under a
similar scheme to that of an elastic reflection full waveform inversion. We
decoupled the micro-seismic source term into spatial and temporal components.
We then invert for the decoupled source and the velocity models with a nested
approach. In the proposed method, the equivalent decoupled source term
represents the source mechanisms with the different features of the source
images. However, the sensitivity of the source images to the acquisition
system affects the inversion results. To mitigate the cycle skipping problem
caused by an unknown source function, a modified objective function, that is
independent of the source time function, is applied by convolving the ground
truth and the predicted data with reference traces, through which the
necessity of a source origin time information is mitigated. The inversion for
all three parameters yields reasonably good results and the low wavenumber
part of the model is retrieved well enough to accurately locate the passive
sources. The low resolution in the velocity model is due to the nature of the
passive seismic experiment setup. For the modified Marmousi model, the
proposed method is successful in finding all the source locations, as well as
in updating a kinematically equivalent velocity model.

(a) $\delta\alpha,\left(1,1,0\right)$

(b) $\delta\beta,\left(1,1,0\right)$

(c) $\delta\alpha,\left(-1,-1,0\right)$

(d) $\delta\beta,\left(-1,-1,0\right)$

(e) $\delta\alpha,\left(0,0,1\right)$

(f) $\delta\beta,\left(0,0,1\right)$

(g) $\delta\alpha,\left(1,-1,0\right)$

(h) $\delta\beta,\left(1,-1,0\right)$

(i) $\delta\alpha,\left(1,0,0\right)$

(j) $\delta\beta,\left(1,0,0\right)$

(k) $\delta\alpha,\left(0,1,0\right)$

(l) $\delta\beta,\left(0,1,0\right)$

(m) $\delta\alpha,\left(1,-2,0\right)$

(n) $\delta\beta,\left(1,-2,0\right)$

(o) $\delta\alpha,\left(-1,2,0\right)$

(p) $\delta\beta,\left(-1,2,0\right)$

(q) $\delta\alpha,\left(-2,1,0\right)$

(r) $\delta\beta,\left(-2,1,0\right)$

Figure 1: Inverted source images with receivers on all four boundaries. Moment
tensors $M_{i}={M_{11},M_{22},M_{12}}$ used for forward modeling: (a)(b)
$M_{1}={1,1,0}$; (c)(d) $M_{2}={-1,-1,0}$; (e)(f) $M_{3}={0,0,1}$; (g)(h)
$M_{4}={0,-1,0}$; (i)(j) $M_{5}={1,0,0}$; (k)(l) $M_{6}={0,1,0}$; (m)(n)
$M_{7}={1,-2,0}$; (o)(p) $M_{8}={-1,2,0}$; (q)(r) $M_{9}={-2,1,0}$.

(a) $\delta\alpha,\left(1,1,0\right)$

(b) $\delta\beta,\left(1,1,0\right)$

(c) $\delta\alpha,\left(-1,-1,0\right)$

(d) $\delta\beta,\left(-1,-1,0\right)$

(e) $\delta\alpha,\left(0,0,1\right)$

(f) $\delta\beta,\left(0,0,1\right)$

(g) $\delta\alpha,\left(1,-1,0\right)$

(h) $\delta\beta,\left(1,-1,0\right)$

(i) $\delta\alpha,\left(1,0,0\right)$

(j) $\delta\beta,\left(1,0,0\right)$

(k) $\delta\alpha,\left(0,1,0\right)$

(l) $\delta\beta,\left(0,1,0\right)$

(m) $\delta\alpha,\left(1,-2,0\right)$

(n) $\delta\beta,\left(1,-2,0\right)$

(o) $\delta\alpha,\left(-1,2,0\right)$

(p) $\delta\beta,\left(-1,2,0\right)$

(q) $\delta\alpha,\left(-2,1,0\right)$

(r) $\delta\beta,\left(-2,1,0\right)$

Figure 2: Inverted source images with only surface receivers. Moment tensors
$M_{i}={M_{11},M_{22},M_{12}}$ used for forward modeling: (a)(b)
$M_{1}={1,1,0}$; (c)(d) $M_{2}={-1,-1,0}$; (e)(f) $M_{3}={0,0,1}$; (g)(h)
$M_{4}={0,-1,0}$; (i)(j) $M_{5}={1,0,0}$; (k)(l) $M_{6}={0,1,0}$; (m)(n)
$M_{7}={1,-2,0}$; (o)(p) $M_{8}={-1,2,0}$; (q)(r) $M_{9}={-2,1,0}$.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

Figure 3: Radiation pattern analysis of P- and S-wave with respect to aperture
(acquisition angles). Red lines: P-wave sensitivity; blue lines: S-wave
sensitivity. Moment tensors $M_{i}={M_{11},M_{22},M_{12}}$ used for forward
modeling: (a) $M_{1}={1,1,0}$; (b) $M_{2}={-1,-1,0}$; (c) $M_{3}={0,0,1}$; (d)
$M_{4}={1,-1,0}$; (e) $M_{5}={1,0,0}$; (f) $M_{6}={0,1,0}$; (g)
$M_{7}={1,-2,0}$; (h) $M_{8}={-1,2,0}$; (i) $M_{9}={-2,1,0}$.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 4: Original (a) Vp and (b) Vs 2D velocity profiles; initial (c) Vp and
(d) Vs; inverted (e) Vp and (f) Vs 2D velocity profiles.

(a)

(b)

Figure 5: Vertical profiles at $x=3km$ of (a) Vp and (b) Vs. Magenta, blue and
red curves represent the original, initial and retrieved vertical velocity
profiles, respectively. Figure 6: Nested approach workflow with equations for
each step. Magenta: source image update, Equations 10 and 11; blue: velocity
model update, Equation 17; red: source wavelet update, Equations 12.

(a) $\delta\alpha$ of source $No.1$, $M_{3}$

(b) $\delta\beta$ of source $No.1$, $M_{3}$

(c) $\delta\alpha$ of source $No.2$, $M_{4}$

(d) $\delta\beta$ of source $No.2$, $M_{4}$

(e) $\delta\alpha$ of source $No.3$, $M_{1}$

(f) $\delta\beta$ of source $No.3$, $M_{1}$

(g) $\delta\alpha$ of source $No.4$, $M_{2}$

(h) $\delta\beta$ of source $No.4$, $M_{2}$

(i) $\delta\alpha$ of source $No.5$, $M_{5}$

(j) $\delta\beta$ of source $No.5$, $M_{5}$

Figure 7: Inverted source images (left) $\delta\alpha$ (right) $\delta\beta$;
True source locations $(x,z)$ are marked as red stars: (1) $(2.4km,1.55km)$,
(2) $(2.6km,1.55km)$, (3) $(2.8km,1.55km)$, (4) $(3.0km,1.55km)$, (5)
$(3.2km,1.55km)$. Blue stars refers to the initial source locations using TRI
method.

(a)

(b)

(c)

(d)

(e)

Figure 8: Inverted source functions for source (a) No.1 (b) No.2 (c) No.3 (d)
No.4 (e) No.5. Black: $w_{11}$, blue: $w_{22}$, red: $w_{12}$.

## References

  * [1] Keiiti Aki and Paul G. Richards.  Quantitative seismology.  2002\.
  * [2] Tariq Alkhalifah and Zedong Wu.  The natural combination of full and image-based waveform inversion.  Geophysical Prospecting, 64(1):19–30, 2015.
  * [3] D Anikiev, J Valenta, F Staněk, and L Eisner.  Joint location and source mechanism inversion of microseismic events: Benchmarking on seismicity induced by hydraulic fracturing.  Geophysical Journal International, 198(1):249–258, 2014.
  * [4] Brad Artman, Igor Podladtchikov, and Ben Witten.  Source location using time-reverse imaging.  Geophysical Prospecting, 58(5):861–873, 2010.
  * [5] J-F Aubry, M Tanter, J Gerber, J-L Thomas, and M Fink.  Optimal focusing by spatio-temporal inverse filter. II. Experiments. Application to focusing through absorbing and reverberating media.  The Journal of the Acoustical Society of America, 110(1):48–58, 2001.
  * [6] Jyoti Behura.  Expedited FWI of microseismic data for mapping velocity changes and hypocenters.  In SEG Technical Program Expanded Abstracts 2015, pages 2461–2466. Society of Exploration Geophysicists, 2015.
  * [7] Biondo Biondi and Thomas Tisserant.  3D angle-domain common-image gathers for migration velocity analysis.  Geophysical Prospecting, 52(6):575–591, 2004.
  * [8] S Bose, H-P Valero, Q Liu, RG Shenoy, and A Ounadjela.  An automatic procedure to detect microseismic events embedded in high noise.  In SEG Technical Program Expanded Abstracts 2009, pages 1537–1541. Society of Exploration Geophysicists, 2009.
  * [9] Kit Chambers, Ben DE Dando, Glenn A Jones, Raquel Velasco, and Stephen A Wilson.  Moment tensor migration imaging.  Geophysical Prospecting, 62(Vertical Seismic Profiling and Microseismicity Frontiers):879–896, 2014.
  * [10] Yangkang Chen.  Fast waveform detection for microseismic imaging using unsupervised machine learning.  Geophysical Journal International, 215(2):1185–1199, 2018.
  * [11] Yunseok Choi and Tariq Alkhalifah.  Source-independent time-domain waveform inversion using convolved wavefields: Application to the encoded multisource waveform inversion.  Geophysics, 76(5):R125–R134, 2011.
  * [12] FAand Dahlen and Jeroen Tromp.  Theoretical global seismology.  Princeton university press, 1998.
  * [13] AB Droujinine, S Oates, and J Ita.  Elastic full waveform inversion for locations and moment tensors of microseismic events.  In 73rd EAGE Conference and Exhibition incorporating SPE EUROPEC 2011, pages cp–238. European Association of Geoscientists & Engineers, 2011.
  * [14] Alexander Droujinine and Elmar Rothert.  The generalized discrete Radon transform applied to passive source decoupled P/S imaging.  In SEG Technical Program Expanded Abstracts 2005, pages 1284–1287. Society of Exploration Geophysicists, 2005.
  * [15] Peter M Duncan.  Is there a future for passive seismic?  First Break, 23(6), 2005.
  * [16] Leo Eisner, Peter M Duncan, Werner M Heigl, and William R Keller.  Uncertainties in passive seismic monitoring.  The Leading Edge, 28(6):648–655, 2009.
  * [17] Mathias Fink, Didier Cassereau, Arnaud Derode, Claire Prada, Philippe Roux, Mickael Tanter, Jean-Louis Thomas, and François Wu.  Time-reversed acoustics.  Reports on progress in Physics, 63(12):1933, 2000.
  * [18] Dirk Gajewski and Ekkehart Tessmer.  Reverse modelling for seismic event characterization.  Geophysical Journal International, 163(1):276–284, 2005.
  * [19] Thomas Gallot, Stefan Catheline, Philippe Roux, and Michel Campillo.  A passive inverse filter for Green’s function retrieval.  The Journal of the Acoustical Society of America, 131(1):EL21–EL27, 2012.
  * [20] Qiang Guo and Tariq Alkhalifah.  Elastic reflection-based waveform inversion with a nonlinear approach.  Geophysics, 82(6):R309–R321, 2017.
  * [21] JBU Haldorsen, M Milenkovic, N Brooks, C Crowell, and MB Farmani.  A full-waveform, migration-based deconvolution approach to locating micro-seismic events.  In 74th EAGE Conference and Exhibition incorporating EUROPEC 2012, pages cp–293, 2012.
  * [22] Hadi Jamali-Rad, Zijian Tang, Xander Campman, Alexander Droujinine, and Geert Leus.  Research Note: Sparsity-Aware Multiple Microseismic Event Localization Blind to the Source Time-Function.  Geophysical Prospecting, 63(1):70–77, 2014.
  * [23] B Lars G Jonsson, Maarten V de Hoop, Mats Gustafsson, and Vaughan H Weston.  Retrofocusing of acoustic wave fields by iterated time reversal.  SIAM Journal on Applied Mathematics, 64(6):1954–1986, 2004.
  * [24] Jordan Kaderli, Matthew D McChesney, and Susan E Minkoff.  A self-adjoint velocity-stress full-waveform inversion approach to microseismic source estimation.  Geophysics, 83(5):R413–R427, 2018.
  * [25] R Kamei and D Lumley.  Full waveform inversion of repeating seismic events to estimate time-lapse velocity changes.  Geophysical Journal International, 209(2):1239–1264, 2017.
  * [26] J Kummerow.  Using the value of the crosscorrelation coefficient to locate microseismic events.  Geophysics, 75(4):MA47–MA52, 2010.
  * [27] Carene Larmat, Jean-Paul Montagner, Mathias Fink, Yann Capdeville, Arnaud Tourin, and Eric Clévédé.  Time-reversal imaging of seismic sources and application to the great Sumatra earthquake.  Geophysical Research Letters, 33(19), 2006.
  * [28] Carène S Larmat, Robert A Guyer, and Paul A Johnson.  Time-reversal methods in geophysics.  Phys. Today, 63(8):31–35, 2010.
  * [29] Rongrong Lu, M Nafi Toksöz, and Mark E Willis.  Locating microseismic events with time reversed acoustics: a synthetic case study.  In SEG Technical Program Expanded Abstracts 2008, pages 1342–1346, 2008.
  * [30] George A McMechan.  Determination of source parameters by wavefield extrapolation.  Geophysical Journal International, 71(3):613–628, 1982.
  * [31] Nori Nakata and Gregory C Beroza.  Reverse time migration for microseismic sources using the geometric mean as an imaging condition.  Geophysics, 81(2):KS51–KS60, 2016.
  * [32] R-E Plessix.  A review of the adjoint-state method for computing the gradient of a functional with geophysical applications.  Geophysical Journal International, 167(2):495–503, 2006.
  * [33] R Gerhard Pratt, Changsoo Shin, and GJ Hick.  Gauss–Newton and full Newton methods in frequency–space seismic waveform inversion.  Geophysical Journal International, 133(2):341–362, 1998.
  * [34] Paul Sava.  Micro-earthquake monitoring with sparsely sampled data.  Journal of Petroleum Exploration and Production Technology, 1(1):43–49, 2011.
  * [35] Chao Song, Tariq Alkhalifah, and Zedong Wu.  Microseismic event estimation and velocity analysis based on a source-focusing function.  Geophysics, 84(3):KS85–KS94, 2019.
  * [36] Chao Song, Zedong Wu, and Tariq Alkhalifah.  Passive seismic event estimation using multiscattering waveform inversion.  Geophysics, 84(3):KS59–KS69, 2019.
  * [37] Fuxian Song, H Sadi Kuleli, M Nafi Toksöz, Erkan Ay, and Haijiang Zhang.  An improved method for hydrofracture-induced microseismic event detection and phase pickingMicroseismic event detection & phase picking.  Geophysics, 75(6):A47–A52, 2010.
  * [38] Brian Steiner, Erik H Saenger, and Stefan M Schmalholz.  Time reverse modeling of low-frequency microtremors: Application to hydrocarbon reservoir localization.  Geophysical Research Letters, 35(3), 2008.
  * [39] Bingbing Sun and Tariq Alkhalifah.  Automatic wave-equation migration velocity analysis by focusing subsurface virtual sources.  Geophysics, 83(2):U1–U8, 2018.
  * [40] Bingbing Sun and Tariq Alkhalifah.  Adaptive traveltime inversion.  Geophysics, 84(4):U13–U29, 2019.
  * [41] Junzhe Sun, Zhiguang Xue, Tieyuan Zhu, Sergey Fomel, and Norimitsu Nakata.  Full-waveform inversion of passive seismic data for sources and velocities.  In SEG Technical Program Expanded Abstracts 2016, pages 1405–1410. Society of Exploration Geophysicists, 2016.
  * [42] WW Symes and James J Carazzone.  Velocity inversion by differential semblance optimization.  Geophysics, 56(5):654–663, 1991.
  * [43] Mickaël Tanter, Jean-Louis Thomas, and Mathias Fink.  Time reversal and the inverse filter.  The Journal of the Acoustical Society of America, 108(1):223–234, 2000.
  * [44] Michael Tanter, J-F Aubry, J Gerber, J-L Thomas, and M Fink.  Optimal focusing by spatio-temporal inverse filter. I. Basic principles.  The Journal of the Acoustical Society of America, 110(1):37–47, 2001.
  * [45] Felix Waldhauser and William L Ellsworth.  A double-difference earthquake location algorithm: Method and application to the northern Hayward fault, California.  Bulletin of the Seismological Society of America, 90(6):1353–1368, 2000.
  * [46] H Wang and T Alkhalifah.  Time reversal migration for passive sources using a maximum variance imaging condition.  In 79th EAGE Conference and Exhibition 2017, volume 2017, pages 1–5. European Association of Geoscientists & Engineers, 2017.
  * [47] Hanchen Wang and Tariq Alkhalifah.  Microseismic imaging using a source function independent full waveform inversion method.  Geophysical Journal International, 214(1):46–57, 2018.
  * [48] Shen Wang, Feng Chen, Houzhu Zhang, and Yunqing Shen.  Reflection-based full waveform inversion (RFWI) in the frequency domain.  In SEG Technical Program Expanded Abstracts 2013, pages 877–881. Society of Exploration Geophysicists, 2013.
  * [49] Norm Warpinski and others.  Microseismic monitoring: Inside and out.  Journal of Petroleum Technology, 61(11):80–85, 2009.
  * [50] Ru-Shan Wu and Yingcai Zheng.  Non-linear partial derivative and its De Wolf approximation for non-linear seismic inversion.  Geophysical Journal International, 196(3):1827–1843, 2014.
  * [51] O Zhebel and L Eisner.  Simultaneous microseismic event localization and source mechanism determination: 82nd Annual International Meeting, SEG.  Expanded Abstracts, doi, 10:1190, 2012.
  * [52] Tieyuan Zhu.  Time-reverse modelling of acoustic wave propagation in attenuating media.  Geophysical Journal International, 197(1):483–494, 2014.
  * [53] Tariq Alkhalifah, Hanchen Wang, and Oleg Ovcharenko.  MLReal: Bridging the gap between training on synthetic data and real data applications in machine learning.  Artificial Intelligence in Geosciences, 3:101–114, 2022.
  * [54] Hanchen Wang, Tariq Alkhalifah, Umair bin Waheed, and Claire Birnie.  Data-driven microseismic event localization: An application to the Oklahoma Arkoma basin hydraulic fracturing data.  IEEE Transactions on Geoscience and Remote Sensing, 60:1–12, 2021.
  * [55] Chengyuan Deng, Shihang Feng, Hanchen Wang, Xitong Zhang, Peng Jin, Yinan Feng, Qili Zeng, Yinpeng Chen, and Youzuo Lin.  OpenFWI: Large-scale Multi-structural Benchmark Datasets for Full Waveform Inversion.  Advances in Neural Information Processing Systems, 35:6007–6020, 2022.
  * [56] Hanchen Wang and Tariq Alkhalifah.  Direct microseismic event location and characterization from passive seismic data using convolutional neural networks.  Geophysics, 86(6):KS109–KS121, 2021.

##  7 Appendix A: The Gradient calculation with the conventional objective
function

Here we derive the gradients of the objective function, Equation 9, with
respect to the model squared velocity $\bf m(\alpha,\beta)$, source images
$\bf{\delta}\bf m(\delta\alpha,\delta\beta)$ and source function $\bf w$:

$\begin{array}[]{cl}dE(\bf u)&\displaystyle=\frac{\partial E}{\partial\bf
m}d\bf m+\frac{\partial E}{\partial\delta\bf{m}}d\delta\bf{m}+\frac{\partial
E}{\partial\left[w_{ij}\right]}d\bf w\\\ &\displaystyle=\int_{t}dt(\bf u-\bf
d)d\bf u,\end{array}$ (20)

The elastic displacement field when density is constant satisfies

$\bf u_{tt}-\nabla\cdot(\bf m:\nabla u)=\nabla\cdot(\delta\bf m:\bf w),$ (21)

where $:$ is the Frobenius inner product operator. The total derivatives of
the displacement vector $\bf u$ is given by

$d\bf u_{tt}=\nabla\cdot\nabla d\bf u:\bf m+\nabla\cdot d\bf m:\nabla\bf
u+\nabla\cdot d\bf w:\delta\bf m+\nabla d\delta\bf m:\bf w.$ (22)

The perturbation $d\bf u$ can be written as

$d\bf u=(\frac{\partial^{2}}{\partial t^{2}}-\nabla\cdot(\bf
m:\nabla))^{-1}(\nabla\cdot d\bf m:\nabla\bf u+\nabla\cdot d\bf w:\delta\bf
m+\nabla d\delta\bf m:\bf w).$ (23)

Equation 23 contains three terms representing the partial derivatives of the
displacement vector with respect to $d\bf m$, $d\bf w$ and $d\delta\bf m$.
Define the adjoint wavefield for back propagating the data residual into the
media as

$\bf\hat{u}=(\frac{\partial^{2}}{\partial t^{2}}-\nabla\cdot(\bf
m:\nabla))^{-T}(\bf u-\bf d),$ (24)

where $T$ is the adjoint operation.
To calculate the gradient for $\bf m(\alpha,\beta)$, $\delta\bf m$ or $\bf w$,
we substitute $d\bf u$ into Equation 20

$\begin{array}[]{cl}&\displaystyle dE=(\frac{\partial^{2}}{\partial
t^{2}}-\nabla\cdot(\bf m:\nabla))^{-1}(\bf u-\bf d)(\nabla\cdot d\bf
m:\nabla\bf u+\nabla\cdot d\bf w:\delta\bf m+\nabla d\delta\bf m:\bf w)\\\
&\displaystyle=\nabla\cdot(\frac{\partial^{2}}{\partial t^{2}}-\nabla\cdot(\bf
m:\nabla))^{-1}(\bf u-\bf d):\nabla\bf ud\bf
m+\nabla\cdot(\frac{\partial^{2}}{\partial t^{2}}-\nabla\cdot(\bf
m:\nabla))^{-1}(\bf u-\bf d):\delta\bf md\bf w\\\
&\displaystyle+\nabla\cdot(\frac{\partial^{2}}{\partial t^{2}}-\nabla\cdot(\bf
m:\nabla))^{-1}(\bf u-\bf d):\bf wd\delta\bf m.\end{array}$ (25)

To calculate the gradients for each parameter, we set the other two
parameters’ total derivatives to be zero (do not change in the current
iteration). For example, to calculate the gradient for $\bf m(\alpha,\beta)$,
set $d\delta\bf m=0$ and $d\bf w=0$

$\nabla E_{\bf m}=\frac{dE}{d\bf
m}=\int_{t}dt\nabla\cdot(\frac{\partial^{2}}{\partial t^{2}}-\nabla\cdot(\bf
m:\nabla))^{-1}(\bf u-\bf d):\nabla\bf u;$ (26)

Similarly, by setting $d\bf m=0$ and $d\bf w=0$, we can obtain the gradient
for $\delta\bf m$

$\nabla{E}_{\delta\bf m}=\frac{dE}{d\delta\bf
m}=\int_{t}dt\nabla\cdot(\frac{\partial^{2}}{\partial t^{2}}-\nabla\cdot(\bf
m:\nabla))^{-1}(\bf u-\bf d):\bf w,$ (27)

and we get the gradient for $\bf w$ by setting $d\bf m=0$ and $d\bf m=0$

$\nabla{E}_{\bf w}=\frac{dE}{d\bf w}=\int_{\bf x}d{\bf
x}\nabla\cdot(\frac{\partial^{2}}{\partial t^{2}}-\nabla\cdot(\bf
m:\nabla))^{-1}(\bf u-\bf d):\delta\bf m$ (28)

##  8 Appendix B: The Gradient for The Source Independent Objective Function

We derive the gradient of the source function independent objective function,
Equation 16, with respect to elastic velocities $c(x,z)=[V_{p};V_{s}]$, under
the constant density assumption. In this case, the objective function without
the optional total variation term is given by

$\
E_{ind}=\sum_{i}^{nr}{\left\|r_{i}\right\|^{2}}=\sum_{i}^{nr}{\left\|{d}_{i}^{pred}*{d}_{ref}^{obs}-{d}_{i}^{obs}*{d}_{ref}^{pred}\right\|^{2}},$
(29)

Taking the derivative of Equation 29 with respect to $c(x,z)$. The gradient
can be written as

$\frac{\partial E_{ind}}{\partial
c}=\sum_{i}^{nr}{\left[\left(\frac{\partial{d}_{i}^{pred}}{\partial
c}*{d}_{ref}^{obs}\right)\cdot
r_{i}-\left({d}_{i}^{obs}*\frac{\partial{d}_{ref}^{pred}}{\partial
c}\right)\cdot r_{i}\right]}.$ (30)

The first convolution term $(\frac{\partial{d}_{i}^{pred}}{\partial
c}*{d}_{ref}^{obs})\cdot r_{i}$ in Equation 30 can be rewritten in an integral
form,

$\left(\frac{\partial{d}_{i}^{pred}}{\partial
c}*{d}_{ref}^{obs}\right)\cdot{r}_{i}=\int_{-\infty}^{+\infty}{\int_{-\infty}^{+\infty}{\frac{\partial{d}_{i}^{pred}}{\partial
c}\left(t-\tau\right)}{d}_{ref}^{obs}\left(\tau\right)r_{i}\left(t\right)}d\tau
dt.$ (31)

Let $\eta=t-\tau$, so that, equation 31 can be written as,

$\left(\frac{\partial{d}_{i}^{pred}}{\partial
c}*{d}_{ref}^{obs}\right)\cdot{r}_{i}=-\int_{-\infty}^{+\infty}{\int_{-\infty}^{+\infty}{\frac{\partial{d}_{i}^{pred}}{\partial
c}\left(\eta\right)}{d}_{ref}^{obs}\left(t-\eta\right)r_{i}\left(t\right)}d\eta
dt.$ (32)

The partial derivative $\frac{\partial{d}_{i}^{pred}}{\partial
c}\left(\eta\right)$ is now independent of $t$ so we move it outside of the
integration of $t$, which gives

$\left(\frac{\partial{d}_{i}^{pred}}{\partial
c}*{d}_{ref}^{obs}\right)\cdot{r}_{i}=-\int_{-\infty}^{+\infty}{\frac{\partial{d}_{i}^{pred}}{\partial
c}\left(\eta\right)}\left[\int_{-\infty}^{+\infty}{{d}_{ref}^{obs}\left(t-\eta\right)r_{i}\left(t\right)}dt\right]d\eta.$
(33)

This derivative $\frac{\partial{d}_{i}^{pred}}{\partial c}\left(\eta\right)$
represents the convolution between $U_{x,i}$ and $S_{x}$, where $U_{x,i}$ is
the adjoint wavefield of the $ith$ backpropagated residual, $S_{x}$ is the
source wavefield ignited at a virtual source location and $x$ represents the
model space location. Using the convolution form to rewrite Equation 31, we
obtain

$\left(\frac{\partial{d}_{i}^{pred}}{\partial
c}*{d}_{ref}^{obs}\right)\cdot{r}_{i}=-\int_{-\infty}^{+\infty}{\left[\int_{-\infty}^{+\infty}{S_{x}\left(\eta-\tau\right)U_{{x,i}}}\left(\tau\right)d\tau\right]}\left[\int_{-\infty}^{+\infty}{{d}_{ref}^{obs}\left(t-\eta\right)r_{i}\left(t\right)}dt\right]d\eta$
(34)

Again let $\eta-\tau=t$, then equation 34 becomes

$\begin{array}[]{lcl}\left(\frac{\partial{d}_{i}^{pred}}{\partial
c}*{d}_{ref}^{obs}\right)\cdot{r}_{i}&=&\displaystyle\int_{-\infty}^{+\infty}{\left[\int_{-\infty}^{+\infty}{S_{x}\left(t\right)U_{{x,i}}}\left(\eta-t\right)dt\right]}\left[\int_{-\infty}^{+\infty}{{d}_{ref}^{obs}\left(t-\eta\right)r_{i}\left(t\right)}dt\right]d\eta\\\
&=&\displaystyle\int_{-\infty}^{+\infty}{S_{x}\left(t\right)}\left[\int_{-\infty}^{+\infty}{{U_{{x,i}}\left(\eta-t\right)d}_{ref}^{obs}\left(t-\eta\right)r_{i}\left(t\right)d\eta}\right]dt\end{array}$
(35)

Because of the reciprocity of $U_{x,i}$, the inner integration
$\int_{-\infty}^{+\infty}{U_{x,i}\left(\eta-t\right)}d_{ref}^{obs}\left(t-\eta\right)r_{i}\left(t\right)d\eta$
can be considered as a back propagation of the new data residual
$R_{i}^{(1)}=\int_{-\infty}^{+\infty}d_{ref}^{obs}\left(t-\eta\right)r_{i}\left(t\right)dt$.
Thus, the first term of the gradient can be obtained by back-propagating
$R_{i}^{(1)}$ as adjoint sources.
A similar derivation can be applied to the second term of Equation 30,
$-\left({d}_{i}^{obs}*\frac{\partial{d}_{ref}^{pred}}{\partial
c}\right){r}_{i}$. This yields

$-\left({d}_{i}^{obs}*\frac{\partial{d}_{ref}^{pred}}{\partial
c}\right)\cdot{r}_{i}=\int_{-\infty}^{+\infty}{S_{x}\left(t\right)}\left[\int_{-\infty}^{+\infty}{{U_{{x,ref}}\left(\eta-t\right)d}_{i}^{obs}\left(t-\eta\right)r_{i}\left(t\right)d\eta}\right]dt.$
(36)

In Equation 36, the second term of the gradient can be obtained by back-
propagating the second residual
$R_{i}^{(2)}=-\int_{-\infty}^{+\infty}d_{i}^{obs}\left(t-\eta\right)r_{i}\left(t\right)dt$.
Thus, the gradient involves calculating the adjoint wavefields with the
adjoint sources:

$\
R^{(1)}_{i}=d^{obs}_{ref}\bigotimes\left({d}_{i}^{pred}*{d}_{ref}^{obs}-{d}_{i}^{obs}*{d}_{ref}^{pred}\right),$
(37)

at the $i^{th}$ receiver position, and

$\
R^{(2)}_{i}=-d^{obs}_{i}\bigotimes\left({d}_{i}^{pred}*{d}_{ref}^{obs}-{d}_{i}^{obs}*{d}_{ref}^{pred}\right),$
(38)

at the reference receiver position.