Waveform inversion for parameters of microseismic sources in VTI media
Oscar Jarillo Michel* and Ilya Tsvankin, Center for Wave Phenomena, Colorado School of Mines
SUMMARY
Waveform inversion (WI), which has been used primarily for
high-resolution velocity analysis, can also be employed to obtain the source parameters of microseismic events. Here, we
implement WI to estimate the location, origin time, and seismic moment tensor of microseismic sources embedded in heterogeneous VTI (transversely isotropic with a vertical symmetry axis) media. The algorithm operates with multicomponent wavefields modeled using an elastic anisotropic finitedifference code. The gradient of the objective function for the
three classes of parameters is calculated with the adjoint-state
method, whereas the interval VTI parameters are assumed to
be known. Synthetic tests for data recorded by vertical receiver
arrays show that it is possible to tightly constrain all source parameters, if a sufficiently accurate initial model is available. In
particular, the source location in layered VTI media can be
estimated simultaneously with the moment tensor. The resolution of event location, however, somewhat decreases when the
origin time is unknown or there are errors in the VTI parameters.
INVERSE PROBLEM
The wave equation for a dislocation-type source located at
point xs in heterogeneous anisotropic media can be written as
(Aki and Richards, 2002):
„
«
∂[δ(x − xs )]
∂uk
cijkl
= − Mij
S(t) ,
∂xl
∂xj
(1)
where u(x, t) is the displacement field, t is time, cijkl is the
stiffness tensor (i, j, k, l = 1, 2, 3), ρ (x) is the density, Mij
is the seismic moment tensor, S(t) is the source time function,
and δ(x−xs ) is the spatial δ-function. We use finite-difference
(FD) algorithm sfewe in MADAGASCAR to obtain the exact
wavefield from equation 1. P- and SV-waves that propagate
in the [x1 , x3 ]-plane depend on the components M11 , M13 ,
and M33 of the moment tensor, which can be represented as
a function of the fault dip angle θ (Jost and Herrmann, 1989;
Aki and Richards, 2002):
ρ
∂ 2 ui
∂
−
∂t2
∂xj
M11 = −
INTRODUCTION
Applications of waveform inversion in exploration seismology have been mostly limited to velocity tomography. Seismic waveforms, however, contain information that can be used
to constrain other important quantities, such as parameters of
earthquake sources.
Recently, waveform inversion has been extended to elastic and
anisotropic media (Lee et al., 2010; Kamath and Tsvankin,
2013), which makes it appropriate for multicomponent reflection and microseismic data. Also, WI has been employed to
estimate earthquake source parameters in global seismology
(Tromp et al., 2005; Liu and Tromp, 2006; Kim et al., 2011)
and geothermal studies (Morency and Mellors, 2012). These
techniques incorporate the adjoint-state method (Lions, 1972;
Plessix, 2006; Fichtner, 2006, 2009) as a practical way to calculate the gradient of the WI objective function. The main advantage of this method (Talagrand and Courtier, 1987) is that
the gradient for all model parameters is computed using only
two numerical-modeling simulations.
Following the approach described in Kim et al. (2011), Jarillo Michel and Tsvankin (2013) implement WI gradient calculation for the source location xs and moment tensor M of a
microseismic event in a 2D VTI medium. They also show that
the adjoint wavefield can be used to identify the presence of
missed sources and, in general, improve the initial source position for WI. Here, we use the algorithm of Jarillo Michel and
Tsvankin (2013) to develop a 2D waveform-inversion methodology for estimating the microseismic source parameters.
Σu
¯
(c13 − c11 ) sin 2θ ,
2
M13 = Σ u
¯ c55 cos 2θ ,
M33 = −
Σu
¯
(c33 − c13 ) sin 2θ ,
2
(2)
(3)
(4)
where Σ is the fault area, u
¯ is the magnitude of the slip, and
c11 , c13 , c33 , and c55 are the stiffness coefficients in the twoindex Voigt notation.
Here, our goal is to invert for the source coordinates xs1 and xs3 ,
the event origin time t0 , and the moment-tensor elements M11 ,
M13 , and M33 assuming that the velocity model is known.
Hence, the vector of model parameters to be estimated by WI
is defined as:
m = {xs1 , xs3 , t0 , M11 , M13 , M33 } .
(5)
We assume the elastic displacement field u(xs , xrn , t) to be
excited by a single source at xs and recorded by receivers located at xrn (n = 1, 2, ..., N ). The data residuals are measured by the least-squares objective function F, which is minimized by the inversion algorithm:
F(m) =
1
kdpre (m) − dobs k2 ,
2
(6)
where dobs and dpre (m) are the observed and predicted displacements, respectively. We carry out simultaneous inversion
Waveform inversion for parameters of microseismic sources in VTI media
for xs , t0 , and M employing the nondimensionalization approach suggested by Kim et al. (2011). This approach eliminates the difference between the units of different parameters
classes, which makes possible simultaneous parameter updating. Estimation of M, xs , and t0 involves complications typical for WI. In particular, the trial source should be within about
one-half of the predominant wavelength from the actual source
location.
Application of the adjoint-state method
=
N
X
(dobs − dpre )(T − t) δ(x − xrn ) ,
(7)
n=1
where u† is called the “adjoint wavefield.” The “adjoint source”
on the right-hand side of equation 7 consists of the time-reversed
data residuals and is obtained by differentiating F(m) with
respect to the forward wavefield u. The adjoint simulation to
generate u† can be carried out with the forward-modeling code
that solves equation 1 by using the receivers as simultaneous
adjoint sources.
The derivatives of the objective function with respect to the
model parameters are:
Figure 1: Actual source (white dot), trial source (red dot) and
a vertical line of receivers (spacing is 6 m) embedded in a homogeneous VTI medium. The medium parameters are ρ = 2
g/cm3 , VP 0 = 4047 m/s, VS0 = 2638 m/s, = 0.4, and
δ = 0. The actual source is located at x1 = 0.3 km and
x3 = 0.75 km with θ = 0◦ (see equations 2 – 4). For the trial
source, x1 = 0.32 km, x3 = 0.8 km, and θ = 15◦ . Both
events occur at the same time (t0 = 0.049 s).
1
Objective function
The adjoint-state method is designed to efficiently calculate the
derivatives of the objective function with respect to the model
parameters, ∂F(m)/∂m. Computation of the gradient of the
objective function for our problem is discussed in Kim et al.
(2011) and Jarillo Michel and Tsvankin (2013). In addition
to equation 1, the method requires solving the adjoint wave
equation:
!
∂u†k
∂ 2 u†i
∂
ρ
cijkl
−
∂t2 ∂xj
∂xl
1e−02
1e−04
1e−06
g
xsi
∂F
=
=
∂xsi
gt0
T
Z
0
∂F
=
=
∂t0
gMij =
˛
∂[M : † (xts , t)] ˛˛
˛ ts S(T − t) dt , (8)
∂xi
x
Z
0
∂F
=
∂Mij
T
∂S(T − t)
dt ,
M : † (xts , t)
∂t
T
Z
†ij (xts , t) S(T − t) dt ,
0
2
4
6
8
Iteration number
10
12
Figure 2: Change of the normalized objective function F(m)
with iterations for the model in Figure 1.
(9)
(10)
0
where xts is the trial source location, T is the recording time,
† = 21 [∇u† + (∇u† )T ] is the adjoint strain tensor, and M :
† is the double inner product of the tensors M and † . Equations 8 – 10 applied to the 2D problem at hand yield the gradient vector g for the six source parameters:
and t0 (equations 8 and 9) include the double-inner product
M : † (xts , t), which involves summation over all elements
of M. Due to this dependence, stable inversion for xs and t0
requires an accurate initial model for the moment tensor.
NUMERICAL RESULTS
(11)
We test the WI algorithm on wavefields from dislocation-type
sources embedded in homogeneous and layered VTI models.
In all experiments, the observed data represent a single microseismic event recorded by a vertical array of closely spaced
receivers.
For inversion purposes, the derivatives in equations 8 – 10 are
needed only at the trial source position. The derivatives for xs
In the first example, we invert for the parameters xs1 , xs3 , M11 ,
M13 , and M33 with the origin time t0 fixed at the actual value
(Figure 1). The objective function becomes practically negligi-
g=
∂F ∂F ∂F ∂F
∂F
∂F
,
,
,
,
,
∂xs1 ∂xs3 ∂t0 ∂M11 ∂M13 ∂M33
ff
.
Waveform inversion for parameters of microseismic sources in VTI media
0.325
0.32
15
0.31
s
x1 (km)
0.315
10
M11 (GN.m)
0.305
0.3
0.295
0
2
4
6
8
Iteration number
10
12
5
(a)
0.81
0
0
0.8
2
4
s
x3 (km)
0.79
6
8
Iteration number
10
12
10
12
10
12
(a)
0.78
2
0.77
0
0.76
−2
0.75
−4
2
4
6
8
Iteration number
10
M33 (GN.m)
0.74
0
12
(b)
−6
−8
−10
xs1
xs3
Figure 3: Change of the source coordinates (a)
and (b)
with iterations for the model in Figure 1. The actual values are
marked by dashed lines.
−12
−14
−16
0
ble after about 10 iterations (Figure 2). The source coordinates
and all elements of the moment tensor are estimated with high
accuracy (Figures 3 and 4); the errors in xs1 and xs3 are on the
order of centimeters. Evidently, for the 2D problem the direct
P and SV arrivals recorded in a single borehole provide sufficient information to constrain the parameters t0 and M.
In the last test we evaluate the influence of velocity errors on
the inversion results by distorting the anisotropy coefficient for the homogeneous VTI model in Figure 1. The algorithm
4
6
8
Iteration number
(b)
18
16
14
M13 (GN.m)
Next, we use the five-layer model in Figure 5 and assume that
all source parameters (xs , t0 , and M) are unknown. Although
the origin time t0 for the actual and trial sources coincides, it
varies with iterations as WI tries to match the observed and
predicted data. In general, the parameters xs and t0 are decoupled and should not generate trade-offs during inversion. However, when both xs and t0 are unknown, the objective function
oscillates and cannot be reduced as much as in the previous
example. Also, as the trial model approaches the actual one,
P- and SV-traveltime shifts produced by further perturbations
in xs and t0 become too small to guide the inversion toward
the actual model. Hence, the obtained coordinate xs1 and time
t0 are slightly distorted (Figures 6 and 7). The small error in
xs1 is eliminated if t0 is known.
2
12
10
8
6
4
0
2
4
6
8
Iteration number
(c)
Figure 4: Change of the moment-tensor elements (a) M11 , (b)
M33 , and (c) M13 with iterations for the model in Figure 1.
The actual values are marked by dashed lines.
Waveform inversion for parameters of microseismic sources in VTI media
0.038
t0 (s)
0.037
0.036
0.035
0.034
0
10
20
30
Iteration number
40
50
Figure 7: Change of the origin time t0 with iterations for the
model in Figure 5.
Figure 5: Actual (white dot) and trial (red dot) sources and a
vertical receiver array in a five-layer VTI model. The receiver
geometry is the same as in Figure 1. The parameters ρ = 2
kg/m3 , = 0.4, and δ = 0 are the same throughout the model
but the vertical P- and SV-vertical velocities vary from layer to
layer. The actual source is located at x1 = 0.3 km, x3 = 0.75
km with θ = 0◦ . The trial source is located at x1 = 0.32 km,
x3 = 0.8 km with θ = 15◦ . Both events occur at the same
time (t0 = 0.035 s).
0.325
estimates the parameters xs and M, whereas the origin time
t0 is fixed at the actual value. An error in changes the Pwave horizontal velocity Vhor , which should influence estimation of the horizontal source coordinate xs1 . After a fast initial
decrease, the objective function flattens out at a larger value
than that in Figure 2. Still, errors in the source coordinates
are relatively small (about 5 m for xs1 ), and the elements M11
and M33 of the moment tensor are also recovered with sufficient accuracy. However, there is a more significant error (over
20%) in the element M13 , which is most sensitive to the quality of waveform matching.
0.32
CONCLUSIONS
0.31
s
x1 (km)
0.315
0.305
0.3
0.295
0
10
20
30
Iteration number
40
50
(a)
We presented a WI methodology for simultaneous estimation
of the source parameters (location xs , origin time t0 , and moment tensor M) of microseismic events from multicomponent
data. The waveform-inversion algorithm operates with the elastic anisotropic wave equation and is designed for dislocationtype sources embedded in 2D heterogeneous VTI media. The
gradient of the objective function is calculated with the adjointstate method followed by the nondimensionalization approach
designed to handle updating for different parameter classes.
0.81
0.8
s
x3 (km)
0.79
0.78
0.77
0.76
0.75
0.74
0
10
20
30
Iteration number
40
50
Synthetic tests were performed for data recorded by a dense
vertical array of two-component receivers in homogeneous and
horizontally layered VTI media. If the initial model is located
within the basin of convergence, WI accurately estimates the
parameters xs , t0 , and M. To evaluate the influence of errors
in the velocity model, we carried out inversion for a homogeneous VTI medium using an inaccurate anisotropy parameter
(0.3 instead of 0.4). Predictably, a distortion in propagates
into the horizontal source coordinate, but the error in xs1 is not
significant. Also, there is a 20% error in the element M13 of
the moment tensor, whereas M11 and M33 are recovered with
high accuracy.
(b)
Figure 6: Change of the source coordinates (a) xs1 and (b) xs3
with iterations for the model in Figure 5.
Although the VTI parameters were assumed known, they can
be included in the adjoint calculation at little additional cost.
Reconstruction of the velocity model, however, can generate
trade-offs between the source and medium parameters.
Waveform inversion for parameters of microseismic sources in VTI media
REFERENCES
Aki, K., and P. G. Richards, 2002, Quantitative seismology,
second ed.: University Science Books.
Fichtner, A., 2006, The adjoint method in seismology: I. Theory: Physics of the Earth and Planetary Interiors, 157, 86–
104.
——–, 2009, Full seismic waveform inversion for structural
and source parameters: PhD thesis, Ludwig-MaximiliansUniversit¨at M¨unchen.
Jarillo Michel, O., and I. Tsvankin, 2013, Gradient computation for full-waveform inversion of microseismic data in
VTI media: SEG Technical Program Expanded Abstracts,
2238–2242.
Jost, M. L., and R. B. Herrmann, 1989, A students guide to and
review of moment tensors: Seismological Research Letters,
60, 37–57.
Kamath, N., and I. Tsvankin, 2013, Full-waveform inversion
of multicomponent data for horizontally layered VTI media: Geophysics, 78, WC113–WC121.
Kim, Y., Q. Liu, and J. Tromp, 2011, Adjoint centroid-moment
tensor inversions: Geophysical Journal International, 186,
264–278.
Lee, H.-Y., J. M. Koo, D.-J. Min, B.-D. Kwon, and H. S. Yoo,
2010, Frequency-domain elastic full waveform inversion
for VTI media: Geophysical Journal International, 183,
884–904.
Lions, J., 1972, Nonhomogeneous boundary value problems
and applications: Springer Verlag, Berlin.
Liu, Q., and J. Tromp, 2006, Finite-frequency kernels based
on adjoint methods: Bulletin of the Seismological Society
of America, 96, 2383–2397.
Morency, C., and R. J. Mellors, 2012, Full moment tensor and
source location inversion based on full waveform adjoint
inversion: application at the Geysers geothermal field: SEG
Technical Program Expanded Abstracts, 532, 1–5.
Plessix, R.-E., 2006, A review of the adjoint-state method for
computing the gradient of a functional with geophysical
applications: Geophysical Journal International, 167, 495–
503.
Talagrand, O., and P. Courtier, 1987, Variational assimilation
of meteorological observations with the adjoint vorticity
equation. I: Theory: Quarterly Journal of the Royal Meteorological Society, 113, 1311–1328.
Tromp, J., C. Tape, and Q. Liu, 2005, Seismic tomography, adjoint methods, time reversal and banana-doughnut kernels:
Geophysical Journal International, 160, 195–216.
© Copyright 2024 ExpyDoc
