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 2018 ExploreDoc