Downloaded 07/06/17 to 172.219.108.51. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

GEOPHYSICS, VOL. 82, NO. 5 (SEPTEMBER-OCTOBER 2017); P. KS71–KS83, 10 FIGS. 10.1190/GEO2016-0561.1

Image-domain velocity inversion and event location for microseismic monitoring

Ben Witten1 and Jeffrey Shragge1

of the arrivals. We have developed a full-wavefield adjoint-state method for locating seismic events while inverting for P- and Swave velocity models that optimally focus multiple complementary images of recorded seismic events. This method requires neither picking nor initial estimates of event location or origin time. Because the inversion relies on (image domain) residuals that satisfy the differential semblance criterion, there is no requirement that the starting model be close to the true velocity. We determine synthetic results derived from a model with conditions similar to a field-acquisition scenario in terms of the number and spatial sampling of receivers and recorded coherent and random noise levels. The results indicate the effectiveness of the methodology by demonstrating a significantly enhanced focusing of event images and a reduction of 95% in event location error from a reasonable initial model.

ABSTRACT Microseismic event locations obtained from seismic monitoring data sets are often a primary means of determining the success of fluid-injection programs, such as hydraulic fracturing for oil and gas extraction, geothermal projects, and wastewater injection. Event locations help the decision makers to evaluate whether operations conform to expectations or parameters need to be changed and may be used to help assess and reduce the risk of induced seismicity. However, obtaining accurate event location estimates requires an accurate velocity model, which is not available at most injection sites. Common velocity updating techniques require picking arrivals on individual seismograms. This can be problematic in microseismic monitoring, particularly for surface acquisition, due to the low signal-to-noise ratio

There are currently two classes of methods for locating microseismic events: one coming from the earthquake seismology and the other from the exploration reflection seismology communities. The first are the traditional earthquake location methods that require picking event arrivals on individual traces (House, 1987; Rutledge et al., 1994). Modeled traveltimes are calculated by ray tracing from every point in a given velocity model to each sensor. One then solves an optimization problem, whereby the event location is estimated as the point in the model space having the minimum residual between picked and computed arrival times. This can be done using either a deterministic (Thurber and Engdahl, 2000) or a probabilistic (Myers et al., 2007) framework. Although these methods can locate events with a limited number of recordings, the signal-to-noise ratios of most surface microseismic monitoring data are too low for accurate arrival picking, which can lead to large location errors (Pavlis, 1986). The second class of microseismic location methods uses fullwavefield migration techniques adapted from reflection seismic

INTRODUCTION Microseismic monitoring is an increasingly important tool for a variety of subsurface fluid-injection programs, and it is now commonly applied during hydraulic fracturing and waste-water injection. In most cases, the main seismic monitoring objective is to determine accurate locations and magnitudes of induced or triggered (micro) earthquakes (Batchelor et al., 1983; Maxwell and Urbancic, 2001). For hydraulic fracturing, event locations are used to calculate frac length and height (Maxwell, 2014), which help improve the efficiency of an injection program by optimizing well and stage spacing. Microseismic monitoring in hydraulic fracturing and wastewater injection programs is important for ensuring that fluids do not leak into potentially preexisting faults, which can trigger larger and sometimes felt earthquakes (Ellsworth, 2013). Therefore, generating accurate event location estimates is paramount for comprehensively assessing and derisking a fluid-injection program.

Manuscript received by the Editor 1 November 2016; revised manuscript received 30 January 2017; published ahead of production 18 May 2017; published online 05 July 2017. 1 The University of Western Australia, Centre for Energy Geoscience School of Earth Sciences and School of Physics, Perth, Australia. E-mail: bwitten@gmail .com; [email protected]. © 2017 Society of Exploration Geophysicists. All rights reserved. KS71

Downloaded 07/06/17 to 172.219.108.51. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

KS72

Witten and Shragge

imaging. The migration process transforms, or focuses, the recorded wavefield from the acquisition surface to the subsurface location from which the energy originated based on the velocity-dependent Green’s function solution of the governing wave equation. One uses this information to form an image at each subsurface model location by stacking the wavefield energy. In an ideal scenario, in which high-quality microseismic data and an accurate velocity model are used, the imaged energy (i.e., signal) will correctly focus at the true event location at an amplitude significantly higher than any accompanying recorded noise component. Although migration methods can be applied to either single- or multicomponent recordings, the use of multicomponent data can lead to the generation of a series of complementary images from different wavefield components (i.e., P- and S-waves). Kirchhoff migration (also called back-projection) is one approach for generating images of microseismic events (e.g., Kao and Shan, 2004; Baker et al., 2005; Ishii et al., 2005). Similar to pick-based methods, this method requires calculating the traveltimes from each receiver to each point in the model. The wavefield data are then stacked along these isosurfaces at each model point to produce an image. This stacking is done for sliding time windows or for multiple time samples in the data. The possible source locations are represented by model space locations, where energy stacks coherently to produce an image-domain amplitude greater than the background level. An alternate microseismic event location approach is to use wave-equation migration (or back-propagation) methods. These techniques use the recorded data to approximately reconstruct the microseismic wavefield, which is propagated through a velocity model by recursively computing a numerical solution of the governing wave equation. One then applies an imaging condition (e.g., McMechan, 1982; Artman et al., 2010; Nakata and Beroza, 2016) to collapse the time axis to extract energy from the propagated wavefield to form an image of the microseismic event. Although migration techniques are more computationally expensive and require additional sensors than pick-based methods (i.e., tens of sensors per square kilometer usually are required for microseismic monitoring; Eisner et al., 2010; Birkelo et al., 2012), they are less sensitive to noise due to the signal-to-noise enhancement provided by effectively stacking across multiple traces. All seismic event location techniques are sensitive to velocity model error (Billings et al., 1994; Eisner et al., 2009; Husen and Hardebeck, 2010), even for relative location techniques (e.g., Michelini and Lomax, 2004). The pick- and migration-based techniques may have location errors due to incomplete or inaccurate physics used in modeling, either through traveltime calculation or through propagation algorithms. However, although a suitably accurate velocity model is known to be critical for accurate location estimates, all too often, little to no velocity information is available prior to seismic monitoring at injection sites. A common scenario is that only a single or a small number of 1D well logs is available. Moreover, measurements are frequently restricted to only the reservoir interval (i.e., the oil-/gas-bearing geologic interval) and therefore place no constraints on the overburden velocities. Minimal velocity knowledge may lead to unreliable event locations and necessitate introducing a velocity model updating procedure. Velocity updating can be approximately accomplished in a straightforward way, such as applying a bulk shift to the velocity until the event locations conform to expectations (Maxwell,

2014). Alternatively, traditional tomographic inversion provides a quantitative means of simultaneously improving the velocity model and hypocenter estimate (Thurber, 1992; Bardainne and Gaucher, 2010; Grechka et al., 2011). However, inversion methods that rely on picked arrivals may not be suitable for microseismic monitoring applications, particularly surface-based ones, due to the low data signal-to-noise ratios. Thus, we look toward wave-equation migration techniques to provide velocity information. One common approach in reflection seismology is to use information from migrated images to update the velocity model and improve the focal quality of the images. One may generate an image by applying a zero-lag imaging condition (Claerbout, 1971), which correlates the source and receiver wavefields (i.e., the forwardpropagating explosive source wavefield and the back-propagating recorded wavefield reconstructed from the data, respectively) at each point in the model over the full recording time. Although it is known that the quality of the image is a function of the velocity model accuracy, zero-lag reflection images do not provide information for velocity improvement. Subsequent work on imaging conditions has demonstrated that extending the correlation of these two wavefields beyond zero-lag produces an extended image volume (also called an extended image) that is useful for assessing the accuracy of migration velocity models (Rickett and Sava, 2002; Sava and Vasconcelos, 2011). In particular, extended images demonstrably have velocity sensitivity and can be used to invert for velocity models that optimize image focusing through the application of adjoint-state tomography (Tarantola, 1984; Symes and Carazzone, 1991; Shen and Symes, 2008). For microseismic monitoring applications, Artman et al. (2010) develop similar zero-lag imaging conditions to produce multiple images of a single event by correlating different wave modes (i.e., P- and S-wave autocorrelation and crosscorrelation). Witten and Shragge (2015) extend the microseismic imaging conditions and show that the combination of multiple complementary zero-lag images and the P-S extended image provide a velocity quality-control (QC) measure. This suggests that, like in reflection seismology, the inclusion of these migrated images in an adjoint-state inversion formalism would lead to a suitable approach for velocity modeling updating in microseismic monitoring scenarios. There are four key criteria that any surface-based microseismic velocity inversion algorithm should satisfy: (1) be adequate to provide accurate event locations, (2) converge even with a poor starting model, (3) have minimal sensitivity to initial estimates of event locations and origin time, and (4) be robust to surface-field scenarios (i.e., sparse spatial sampling and strong noise). The goal of this paper is to demonstrate that the combination of multicomponent seismic recordings, a suite of microseismic (extended) images, and an adjoint-state tomography formalism leads to a powerful framework for velocity model updating and accurate event location that meets all of these criteria. This paper extends the QC measure of Witten and Shragge (2015) to a formal optimization problem that meets the aforementioned requirements. We use wave-equation imaging and adjointstate tomography to simultaneously invert for P- and S-wave velocity models that produce optimally focused and internally consistent imaged event locations. We use multiple redundant images of events constructed through autocorrelation and crosscorrelation of the Pand S-wavefields to generate image-domain residuals based on the consistency of imaged events. The residuals are derived from what

Downloaded 07/06/17 to 172.219.108.51. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

Microseismic joint location and inversion we call the self-consistency principle: Given the correct velocity model, all images must produce a focal maximum at the same spatial location and the extended image must achieve its maximum at zero lag. The energy at spatial locations in the images where this does not occur forms the image-domain residuals, which we use to jointly invert for multiparameter velocity model updates that optimally focus the suite of images, thereby optimizing event location estimates. We begin by detailing the wave-equation imaging conditions and presenting a multiterm objective function to optimize the consistency amongst a series of images through adjoint-state tomography. Next, we derive the P- and S-wave velocity model gradients for each term in the objective function. To demonstrate the applicability of the inversion procedure, we present 2D synthetic data examples. We first show a noise-free inversion example using sparse surface acquisition, and then we present the results when using a data set comprised of the same signal but contaminated with significant coherent and random noise that mimic field-like conditions. We conclude with a discussion on parameter choice for the inversion algorithm and a comparison of advantages and shortcomings of the proposed methodology with the well-known full-waveform inversion (FWI) methodology (Virieux and Operto, 2009; Wang and Alkhalifah, 2016).

1 ∂2t ui ðx; t; eÞ 2 V i ðxÞ

KS73

− ∇2 ui ðx; t; eÞ ¼ Di ðxr ; t; eÞ;

(1) where ui is the acoustic wavefield, V i is the velocity model, x is the spatial location, e is the event index, Di is the data recorded at sensor locations xr , and subscripts i indicate the P- or S-wave mode. Although not formally correct for particle motion, the acoustic assumption maintains kinematic traveltimes for the direct arrivals, making the assumption suitable for this application. As indicted by the data Di it is necessary to separate wave modes prior to propagation. In the following section, we discuss our method of preprocessing to derive P- and S-wave data estimates from multicomponent recordings. For a given event, e, we produce multiple complementary images by autocorrelating and crosscorrelating the P- and S-wavefield components. We generate zero-lag image volumes I PP and I SS by autocorrelating the P- and S-wavefields, respectively, and an extended image volume I PS by crosscorrelating the P- and S-wavefields spatially and/or temporally shifted with respect to each other and summing over time. The imaging conditions are written as

Z I PP ðx; eÞ ¼

THEORY The adjoint-state inversion methodology developed herein is based on an iterative application of five main steps. First, we inject the recorded P- and S-wave data of isolated events to form wavefields (state variables), which we back propagate through the current P- and Swave velocity model estimates using the acoustic wave equation. Second, we apply a series of imaging conditions that stack over the time axis to produce image volumes of the event. The image amplitudes are a measure of the auto focusing and cross focusing of the P- and S-wave data throughout the model. Third, we generate imagedomain residuals by applying penalty functions to the image volumes that remove energy consistently focused amongst the suite of images for the zero-lag images or around zero lag in the extended image. Fourth, we use the residual energy to develop the adjoint-state variables, which are correlated with the back-propagated wavefield data to produce P- and S-wave velocity gradients. Finally, we use a linesearch method to find an appropriate step length to determine the magnitude of the multiparameter model update. The P- and S-wave model gradients are scaled by the step length and added to the current velocity model. We detail each phase of the inversion process below. Although it is possible to derive similar equations using the elastic wave equation, as done by Shabelansky et al. (2015) for teleseismic converted waves using extended images, we have found that the acoustic approximation is sufficient because we are only interested in determining the source location, and not perfectly reconstructing the wavefield, and this approach is considerably less computationally expensive. Finally, although the examples we present below are in 2D, we emphasize that the methodology as written is equally applicable for 3D applications.

Imaging We begin by looking at the pseudoacoustic imaging conditions. We back-propagate recorded P- and S-wave data with their corresponding time-domain scalar acoustic wave equations:

i ¼ P; S;

T

Z I SS ðx; eÞ ¼

0

T

0

uP ðx; t; eÞuP ðx; t; eÞdt;

(2)

uS ðx; t; eÞuS ðx; t; eÞdt;

(3)

and

I PS ðx; λ; τ; eÞ  Z 0 uP ðx − λ; t − τ; eÞuS ðx þ λ; t þ τ; eÞdt; 0 ; ¼ max T

(4) where x ¼ ½x; y; z are the spatial dimensions, T is the maximum time of the extracted trace data, and λ ¼ ½λx ; λy ; λz  and τ are the spatial and temporal shifts between correlated wavefield components, respectively. Because I PP and I SS are the autocorrelations, there is no velocity sensitivity in the extended images when the origin time is unknown. We elect to remove all negative values from the I PS volume, which has the effect of filtering out strong out-ofphase correlations that disrupt our residuals and inversion procedure. The integrals begin at time T to indicate time reversal or back-propagation. We select the maximum of each zero-lag image as an estimated event location. Therefore, for each event, we have three location estimates, one from each image. Of all the images, I PS affords the highest resolution due to the P- and S-wavefields propagating at different velocities and focusing over a narrower zone. An important advantage of this imaging approach is that no prior estimates of the event location or origin time are needed. When λ ¼ 0 m and τ ¼ 0 s, I PS reduces to the zero-lag image, similar to I PP and I SS . Equation 4 is an extended imaging condition using the acoustic approximation that is an extension of the zero-lag elastic imaging condition of Artman et al. (2010). For a sufficiently accurate velocity model and by compensating for the radiation pattern, the zero-lag I PP , I SS , and I PS images should

Witten and Shragge

KS74

Downloaded 07/06/17 to 172.219.108.51. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

focus at the same correct spatial location and the extended I PS at zero lag (i.e., the image suite is self-consistent). The imaging conditions exploit the fact that for correct velocity models the P- and S-waves simultaneously collapse to the source location during back-propagation, but they will diverge elsewhere.

J¼

1X 2 e

Z  Z τ Z ϵ1 −τ

þ

ϵ2 F2SS I 2PP

þ

ϵ5 F2PP I 2PS

λ

−λ

 F2λ ðλ; τÞI 2PS ðλ; τÞdλdτ

þ

ϵ3 F2PP I 2SS

þ ϵ4 F2SS I 2PS

þ

ϵ6 F2PS I 2PP

ϵ7 F2PS I 2SS

þ

 dx;

(7)

Image-domain residuals The lack of tightness in focus around zero-lag in the extended image-volume and inconsistent focal locations in the zero-lag images are indicators of velocity model error along source-receiver wavepaths for any given event. We wish to use this poorly and/ or inconsistently focused energy in an adjoint-state formalism by defining a set of residuals that are useful for velocity updating. We generate image-domain residuals by use of penalty functions Fj which remove what we consider to be “well-focused” energy in the image volumes. For extended images, this type of penalty function is a differential semblance operator (Symes and Carazzone, 1991; Shen, 2008), which is frequently used in image-domain inversion strategies (Mulder and ten Kroode, 2002; Tang and Biondi, 2011). It is usually defined as an absolute value function in λ- and τ-spaces (Shen and Symes, 2008; Shabelansky et al., 2015; Yang and Sava, 2015), which may be thought of as the inverse of an idealized image. However, for low signal-to-noise data, this can put too much weight on noise at large lags, and therefore we choose a multidimensional Gaussian function centered at zero lag to penalize the extended image volume. This removes the energy at and around zero lag, but it allows us to independently adjust the variances, σ 2i , in the shift dimensions:

 Fλ ¼ 1 − exp

λ2 − x2 2σ x

 λ2y λ2 τ2 − 2 − z2 − 2 : 2σ y 2σ z 2σ τ

(5)

To calculate residuals for the zero-lag images, we create a penalty function to annihilate energy co-located between a pair of images; thus, only inconsistent energy remains. For each zero-lag image, we generate a penalty function Fj



 αj I j Fj ¼ sech ; maxðI j Þ

(6)

where I j is a zero-lag image (i.e., j ¼ PP; SS, or PS) and αj is a dimensionless parameter that determines the broadness of the penalty operator. This penalty function is near zero around the maximum focus in each image and tends to unity away from this location. Shragge et al. (2013) introduce a similar penalty function for 4D seismic inversion. In our case, we create a penalty for each zero-lag image and apply it to the other zero-lag images to remove energy where the respective image pairs are consistent.

Objective function We can now formulate an optimization problem based on minimizing the residuals in the extended image volume and the crosspenalized zero-lag images through differential semblance optimization (Symes, 1993). Thus, we wish to minimize the following objective functional J:

where ϵi are the weighting functions that determine the relative importance of each term in the overall objective function. We omit the dependence on spatial coordinates x and event e for brevity. Because the formulation of equation 7 satisfies the differential semblance criterion (Symes, 1993), the initial velocity need not be close to the true model to achieve model optimization.

Inversion We solve the optimization problem expressed in equation 7 by means of the adjoint-state formalism (Chevant, 1974; Tarantola, 1984), which is based on the computation of “state variables.” In the passive monitoring scenarios, the state variables are the back-propagated P- and S-wavefields, uP and uS (from equation 1). Perturbing the objective function (equation 7) with respect to the state variables yields

Z δJ ¼

½δV P K P ðxÞ þ δV S K S ðxÞdx;

(8)

where K P and K S are the P- and S-wave velocity model gradients, which are obtained by solving

Ki

2 X ¼ 3 Vi e

Z

0 T

∂2t ui ðx; tÞai ðx; T − tÞdtdx;

i ¼ P; S; (9)

where ai is an “adjoint-state variable” obtained by forward propagation of adjoint sources (see Appendix A). As seen in equation 9, the gradients are calculated by correlating the wavefields reconstructed from the back-propagated data and the forward-propagated adjoint sources. Appendix A presents the full derivation (again independent of the dimensionality of x) and definitions of ai . After calculating the P- and S-wave gradients, we determine Pand S-wave step lengths through a multiparameter line-search method (Tang and Ayeni, 2015), whereby the step lengths for Pand S-velocity updates are calculated first independently then jointly. The negative of the gradients scaled by the determined step length are added to the current velocity model: ðlþ1Þ

Vi

ðlÞ

ðlÞ

¼ V i − mi K iðlÞ ;

(10)

where l is the iteration number, V i is the velocity model, and mi is the step length for the P- and S-velocity model updates. To summarize, our algorithm is an iterative five-step process: (1) independently propagate the acoustic P- and S-wavefields (equation 1), (2) apply imaging conditions (equations 2–4), (3) apply penalty functions, Fj and Fλ , to construct residuals (equations 5 and 6), (4) calculate P- and S-wave velocity gradients (equation 9), and (5) determine a step length and update the model (equation 10). We repeat the process until a stopping criterion is met (e.g., residual below a required threshold, no further change in residual, maximum number of iterations reached). Note that although we have derived the formu-

Microseismic joint location and inversion lation with nonzero τ-lags, for computational efficiency, we do not calculate temporal lags (τ ¼ 0 s) for the remainder of the paper.

experimental condition. Figure 1 shows the V P model and V P ∕V S ratio that we use to simulate data from 10 double-couple sources with random orientations, shown as Xs, with a 2D isotropic elastic finitedifference algorithm (Weiss and Shragge, 2013) on a Δx ¼ Δz ¼ 10 m grid. The source locations are constrained to a narrow geologic interval of 90 m. The simulated wavefield data are extracted at synthetic multicomponent (vertical and horizontal motions) seismometers placed on the surface at 300 m intervals. Figure 2a and 2b shows the recorded vertical and horizontal component data. As noted above, the acoustic assumption used in equations 2 and 3 requires the separation of the data into P- and S-wave components. We choose to estimate this by straightforwardly taking the vertical component as the P-wave data and the horizontal component as the S-wave data. However, we recognize that other methods, including polarization steering filters (Galperin, 1984) and Fourierdomain techniques (Dankbaar, 1985), could be applied that may lead to modest improvements. Because the sources are not explosive, the radiation pattern of the event manifests in the image (Artman et al., 2010), when no corrective measures are implemented. Numerous authors have devised methods for removing radiation pattern effects in either the data domain (Eisner et al., 2008; Ay et al., 2012) or the image domain (Saenger, 2011; Anikiev et al., 2014; Chambers et al., 2014). Following Witten and Shragge (2015), we remove the radiation pattern by taking the envelope of data followed by a band-pass filter prior to back-propagation. This ensures that the P- and S-wave data have the same polarity across the array. The only other data preprocessing step is the application of a time-window mask about the arrivals and a trace balance to compensate for amplitude differences related to the radiation pattern. In field-acquisition scenarios, masking or windowing the arrival data would require a detection step prior to extracting the appropriate time window. However, short timewindowed data are necessary because the image quality is a

Forward modeling To demonstrate the effectiveness of the method above, we present a synthetic example based upon a microseismic monitoring scenario. We first present a noise-free example and then add coherent and random Gaussian noise to the data set to create a more realistic Distance (km) 0

1

2

4

5

6

VP (km/s)

4

1 3 2

b)

Distance (km) 0

Depth (km)

3

0

1

2

3

4

5

6 2.2

VP/VS

Depth (km)

0

1

2

2

1.8

Figure 1. (a) True V P and (b) V P ∕V S model with 10 source locations, Xs.

a) 0

2

4

6

8

Trace 10 12

c) 14

16

18

20

0

2

4

6

8

Trace 10 12

e) 14

16

18

20

2

4

6

8

Trace 10 12

14

16

18

20

2

4

6

8

Trace 10 12

14

16

18

20

0 0.5

0.5

1

1 t (s)

t (s)

t (s)

0.5

1.5

1.5

2

2

2.5

0

2

4

6

8

Trace 10 12

2 2.5

d) 14

16

18

20

0

1 1.5

2.5

b)

2

4

6

8

Trace 10 12

f) 14

16

18

20 0

0.5

0.5 0.5

t (s)

1 t (s)

1 t (s)

Downloaded 07/06/17 to 172.219.108.51. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

SYNTHETIC EXAMPLE

a)

KS75

1

1.5

1.5

2

2

2

2.5

2.5

2.5

1.5

Figure 2. Forward modeled (a) vertical and (b) horizontal component data. Data after adding coherent and random noise for (c) vertical and (d) horizontal components. Estimated (e) P- and (f) S-wave data after preprocessing.

Witten and Shragge

ϵ ¼ ½2.0; 0.1; 0.1; 1.0; 1.0; 0.01; 0.01;

(11)

and penalty parameterizations to effectively remove well-focused energy

½αP ; αS ; αPS  ¼ ½1.0; 1.0; 1.0;

and

½σ x ; σ z  ¼ ½0.1; 0.1 km:

(12)

In the “Discussion,” section, we further discuss our strategies for choosing these parameters. Figure 4 shows the inverted V P and V P ∕V S velocity models after applying five iterations of the adjoint-state tomography formalism

0

0

1

Distance (km) 2 3 4

5

6

3.5

1

3 2

VP (km/s)

a)

2.5

b) 0

0

1

Distance (km) 2 3 4

5

6 2.2

1

2

2

1.8

VP/VS

To generate an initial velocity model, we extract a 1D trace from the center of the V P model (Figure 1a) and assume a constant V P ∕V S ratio of 1.7, which is approximately 10% too low on average. We use this to model a single borehole scenario, in which we have logged only V P along the entire depth of the well bore. Figure 3a–3c shows the zero-lag I PP , I SS , and I PS images for the multicomponent data shown in Figure 2a and 2b using the 1D velocity model assumption. Independent event locations are estimated for each imaging condition and event as the maximum value in each image. For these images, and similar ones in the sections below, the estimated event locations are indicated by magenta dots, whereas the true event locations are denoted by green Xs. For reference, the white dots and Xs indicate the estimated and true event locations for the other nine events, respectively. Figure 3d shows a slice through the extended image volume computed for ðjλx ; jλz jÞ ≤ 0.2 km at the I PS estimated event location in Figure 3d (i.e., the magenta dot). As observed in Figure 3, the various images produce different event location estimates and a maximum value in the extended I PS image shifted away from zero lag. The event location in the I PP image is nearly correct because the P-wave velocity model is the most accurate. The I SS event locations exhibit a consistent pattern, but they are imaged too shallow due to the S-wave velocity being too fast. The zero-lag I PS imaging condition produces event locations that are very inconsistent with the other images due to a large error in the V P ∕V S ratio. When viewed collectively, we note that large velocity errors must be present due to the inconsistency of the event locations amongst the images and the poor focusing around zero lag in the extended image volume. We calculate the spatial error for each event as the Euclidean distance from the true source location to the estimated location. The average spatial root-meansquare (rms) error for I PS is 620 m. Although, the I PP and I SS event locations have a lower average error (25 and 230 m, respectively), we take the I PS locations as a reference because these images should produce more highly resolved location estimates and ensure suitable P- and S-wave velocity models.

Inversion results There are many ways to parameterize the inversion, including by changing the spatial locations where one computes the extended

Figure 4. Inverted (a) V P and (b) V P ∕V S model with final estimated I PS event locations using the noise-free data shown in Figure 2a and 2b.

a)

c)

b)

3

3.5

4

1

2

3

3.5

4

1.5 2

d)

Distance (km) 2.5

1 Depth (km)

1.5

Distance (km) 2.5

3

3.5

4

1 Depth (km)

Distance (km) 2.5

Depth (km)

Downloaded 07/06/17 to 172.219.108.51. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

Initial imaging

image volume, the number of spatial and temporal lags, the weights ϵi , and the broadness of the penalty functions, αi and σ i , in equation 7. In the following examples, we calculate the extended image volume at a single spatial location (the estimated event location) for each event and iteration. We use 20 spatial lags in the x- and z-directions sampled at 10 m intervals and zero in τ. We choose ϵi values to favor the optimization of the I PS image at the expense of complete agreement between all the zero-lag images

Depth (km)

function of the overall trace signal-to-noise ratio, and including additional samples not containing arrivals degrades the image.

Depth (km)

KS76

1.5 2

Figure 3. Zero-lag (a) I PP , (b) I SS , and (c) I PS images using the initial 1D velocity model and noise-free data. The magenta dot and green X indicate the estimated and true locations for the shown event images. The white dots and Xs show the estimated and true locations for the other nine events. (d) Slice through an extended image volume at the magenta point in panel (c) for −0.2 ≤ λ ≤ 0.2 km.

developed above. Most of the update is applied to slow the S-wave velocity. The velocity update is very smooth, as expected from the low spatial frequency of the images. The white dots indicate the imaged locations extracted from the zero-lag I PS image, whereas the Xs show the true location. It is evident when comparing Figures 1 and 4 that the inversion procedure has greatly improved the estimated event locations. The average spatial rms error is reduced from 620 m using the initial velocity model to 17 m with the velocity model obtained through the inversion process (Figure 4). Figure 5a–5c shows the zero-lag I PP , I SS , and I PS images when using the inverted velocity model shown in Figure 4. Figure 5d shows a slice through the extended image volume taken at the location of the magenta dot in Figure 5c. Comparing Figures 3 and 5, we observe that the inversion procedure has greatly improved the I SS and I PS images and event location estimate. The spatial errors in I PP and I SS are reduced from 25 to 12 m and 230 to 79 m, respectively. All imaging conditions produce smaller location errors and, when taken as a group, are more internally self-consistent. In particular, the I PS image is much better focused and has a maximum in the extended image closer to zero lag. Figure 6a shows the convergence curve after five iterations of the inversion procedure with the contributions for each term in the objective function (equation 7) in a stacked bar graph. We observe that most of the improvement comes from bringing the zero-lag I PS into better focus and agreement with the I PP and I SS images (terms 4 and 5). The second largest improvement is in term 1, focusing the extended image around zero lag. Figure 7a presents a crossplot between the initial and final I PS rms errors for all 10 events. A point falling along the dashed line would indicate no improvement in location, whereas those below the line are where errors have decreased. This plot clearly shows a significant decrease in errors between the initial and final models.

KS77

and horizontal data components and apply the same preprocessing procedure as the noise-free data example. Figure 2e and 2f shows the preprocessed, noisy P- and S-wave data used in the imaging and inversion process for a single representative event. As expected, the S-wave arrival (Figure 2f) has a significantly higher signal-to-noise ratio than its P-wave counterpart (Figure 2e) due to the double-couple source mechanism. Although it would be relatively easy to pick the S-wave arrivals, doing so with the P-wave data poses a greater challenge and arguably necessitates a migration approach. We follow the same wave-equation migration and tomography procedure for the noisy data as the noise-free example presented above. Figure 8a–8c shows the initial I PP , I SS , and I PS images using the 1D velocity model and noisy data shown in Figure 2e and 2f. Figure 8d presents a slice through the extended image at the estimate event location shown in Figure 8c. Comparing the initial noise-free

Noisy data To produce a more realistic synthetic data example, we now add noise to the data shown in Figure 2a and 2b. We generate correlated noise by exciting sources randomly in time and space along the model’s surface. We then add the correlated noise along with random Gaussian noise to the source arrivals, such that the signal-tonoise, defined as the ratio of the maximum value of the modeled signal and noise for each event, is set at 0.5. Figure 2c and 2d shows the vertical and horizontal components of the noisy data, respectively. We again associate the P- and S-wave data with the vertical

2.5

Distance (km) 3 3.5

b) 4

2.5

1

2

c) 4

2.5

1 Depth (km)

1.5

Distance (km) 3 3.5

Figure 6. Convergence curve as a stacked bar graph of each term in equation 7 using (a) noise-free and (b) noisy data.

1.5 2

Distance (km) 3 3.5

d) 4

1 Depth (km)

a)

Depth (km)

Downloaded 07/06/17 to 172.219.108.51. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

Microseismic joint location and inversion

1.5 2

Figure 5. The (a) I PP , (b) I SS , and (c) I PS zero-lag images using the inverted model shown in Figure 4 and noise-free data. The magenta dot and green X indicate the estimated and true locations for the shown event images. The white dots and Xs show the estimated and true locations for the other nine events. (d) Slice through an extended image volume at the magenta point in panel (c) for −0.2 ≤ λ ≤ 0.2 km.

Witten and Shragge

and noisy data images (Figures 3 and 8), we note the presence of additional scatter in the noisy images due to the lower signal-tonoise data. The noisy I PP event locations have an error of 110 m compared with only 25 m for the noise-free example. The I SS locations do not change substantially with errors of 243 m compared with 230 m for the noise-free case. However, the I PS event locations have a smaller error in the noisy data, 537−620 m. We apply the same inversion procedure with the same parameters as the noise-free example. Figure 9 shows the resulting V P model and V P ∕V S ratio. The resulting velocity models are similar to the noise-free case. Figure 10 shows imaging results as in Figure 8 except using the inverted velocity model (Figure 9). Although we have not improved the I PP event locations much (110 m error decreased to 93 m), the I SS and I PS locations have substantially improved. The average rms error in the I SS locations decreases from 243 to 99 m. The average I PS rms location error has dropped from 541 to 29 m, a decrease of 95%. Figure 7b shows a crossplot of the I PS initial and final estimated location errors for all 10 events. Again, we see a large decrease in error with the events plotting close to the horizontal axis. Figure 6b shows the convergence curve for the noisy data as a stacked bar graph for each term in equation 7. Similar to the noisefree case (Figure 6a), the greatest decrease is from improved focusing of I PS . However, in the noisy case, the residual cannot be driven

as close to zero largely due to the effects of coherent noise in the images.

DISCUSSION

Final error (m)

Final error (m)

The inversion procedure and corresponding results are dependent on the parameterization of the extended images and objective function. In our experience, we find that the particular choice depends largely on the specific fitting goal and the perceived quality of the initial P- and S-wave velocity models. Because the I PS image provides the highest resolution and is dependent on the P- and S-wave velocity models, we suggest that optimizing this image be the primary focus of the inversion procedure. To improve the I PS image, the crosspenalized autocorrelation terms (terms 2 and 3) should not be weighted strongly because they ensure only co-spatial but not cotemporal focusing between I PP and I SS . Conversely, the extended image contribution in term 1 should be weighted more heavily to optimize co-spatial and co-temporal focusing. Terms 4−7 should be weighted based upon the perceived quality of the initial images. Although we could use only the extended image (i.e., term 1), there are drawbacks to this approach. Calculating the extended images and adjoint sources is a computationally expensive operation. The computational cost is a function of the number of spatial locations at which the extended image is calculated and the number of lags. A full analysis of the computational cost is detailed in Sava and Vasconcelos (2011). If the initial imaged event is poorly focused, one b) a) should include additional spatial locations and/ 600 or a greater number of spatial lags. We have 600 found that using a multitude of terms is beneficial for achieving optimal inversion results with al400 400 most no additional computational expense. In addition, one may wish to inspect the initial images for quality to determine weighting for the 200 200 terms along with assumptions about the quality of the initial velocity models. For example, we know the P-wave model is better constrained, 0 0 and therefore we choose to upweight the terms 0 200 400 600 0 200 400 600 penalized by the I PP image. Conversely, we obInitial error (m) Initial error (m) serve very poor focusing in the I PS image and thus opt to weight these very weakly. Figure 7. Crossplots of the I PS initial and final rms errors for (a) noise-free and (b) noisy Comparing the procedure presented herein data for 10 forward modeled events. The dashed line indicates no change in error. Points falling below the line indicate improved event location estimates. All events observe a with another popular technique, FWI, we note decreased error through the inversion procedure. certain advantages of the image-domain methods

** **

2.5

Distance (km) 3 3.5

b) 4

2.5

Depth (km)

2

Distance (km) 3 3.5

c) 2.5

4

1.5 2

Distance (km) 3 3.5

d) 4

1

1

1 1.5

* *** *

*

Depth (km)

a)

Depth (km)

Downloaded 07/06/17 to 172.219.108.51. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

KS78

1.5 2

Figure 8. The (a) I PP , (b) I SS , and (c) I PS zero-lag images using the initial 1D velocity model and noisy data shown in Figure 2e and 2f. The magenta dot and green X indicate the estimated and true locations for the shown event images. The white dots and Xs show the estimated and true locations for the other nine events. (d) Slice through an extended image volume at the magenta point in panel (c) for −0.2 ≤ λ ≤ 0.2 km.

for surface microseismic monitoring. Although both methods can handle sparse data and, depending on the parameterization, can be tuned to prioritize event location, FWI results may be biased based on the initial estimation of the source location and origin time. Using wave-equation migration imaging means that no initial source estimate is required. It is a well-known limitation of FWI that the initial model must be sufficiently accurate (depending on the frequency content of the data) to avoid local minima in the inversion (Virieux and Operto, 2009). In contrast, image-domain inversion is significantly less constrained by this limitation. Additionally, FWI may suffer from low signal-to-noise data because the algorithm will try to fit the noise on individual traces, whereas the image-domain method enhances the signal-to-noise through stacking-based migration principles. The main drawback of the image-domain method, though, is the limited resolution of the resulting models. The method described above cannot produce the same quality resolution as FWI, but we believe that generating a model that produces accurate event locations is sufficient for most microseismic monitoring applications given the constraints of data quality and acquisition geometry commonly observed in surface seismic monitoring.

KS79

CONCLUSION We derive an adjoint-state tomography method to jointly invert for P- and S-wave velocity models by optimizing image-space focusing for multicomponent data recorded on a passive seismic array. We construct an objective function that incorporates multiple images of the seismic event to ensure optimal focusing and consistency between the various images. Our method is based upon wave-equation migration and adjoint-state tomography, and it is thus applicable for events recorded on large numbers of sensors. It does not require arrival picking, estimation of source location or onset time. Because there is no picking, it is suitable for low signal-to-noise data, such as surface-based microseismic or induced seismic monitoring. Using realistic synthetic examples, we have shown that our methodology can generate elastic velocity models that result in location errors of the order of tens of meters. Relative to initial model estimates, the location errors are reduced by 95%. Beyond microseismic monitoring, this technique can be applied in any emission monitoring program from engineering applications, such as nondestructive testing through global-scale seismology.

ACKNOWLEDGMENTS

a)

The authors thank the staff and students of the Centre for Energy Geosciences for useful discussions. In particular, we thank T. Potter for help implementing boundary conditions and R. Kamei for insightful comments on the manuscript. We also thank W. Debski and an anonymous reviewer for their comments. This work was supported by an Australian Government Research Training Program Scholarship for Ben Witten and an Australian Society of Exploration Research Foundation Grant. We acknowledge the support of the University of Western Australia Reservoir Management Consortium. Jeffrey Shragge acknowledges the support of Woodside Petroleum Ltd. The reproducible numeric examples use the Madagascar open-source package (http://www.reproducibility.org).

Distance (km) 0

1

2

3

4

5

6 VP (km/s)

3.5

1

3 2 2.5

b)

Distance (km) 2 3 4

1

5

6 2.2

1

2

2

1.8

APPENDIX A

VP/VS

Depth (km)

0

0

P- AND S-WAVE VELOCITIES GRADIENT DERIVATION

Figure 9. Inverted (a) V P and (b) V P ∕V S model with final estimated I PS event locations using the noisy data shown in Figure 2e and 2f.

a)

b)

Distance (km) 2.5

3

3.5

4

1

2

3

3.5

4

1.5 2

d)

Distance (km) 2.5

1 Depth (km)

1.5

c)

Distance (km) 2.5

We follow a similar strategy to that outlined in Plessix (2006) and Shabelansky et al. (2015) to derive P- and S-wave velocity model gradients, K P and K S , for each term in the objective functional

3

3.5

4

1 Depth (km)

Depth (km)

0

Depth (km)

Downloaded 07/06/17 to 172.219.108.51. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

Microseismic joint location and inversion

1.5 2

Figure 10. The (a) I PP , (b) I SS , and (c) I PS zero-lag images using the final inverted model shown in Figure 9 and noisy data shown in Figure 2e and 2f. The magenta dot and green X indicate the estimated and true locations for the shown event images. The white dots and Xs show the estimated and true locations for the other nine events. (d) Slice through an extended image volume at the magenta point in panel (c) for −0.2 ≤ λ ≤ 0.2 km.

Witten and Shragge

KS80

Downloaded 07/06/17 to 172.219.108.51. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

(equation 7), for a single event e. For each of seven terms, the derivation follows a similar procedure with multiple terms being identical save for a switch of subscripts P and S. Therefore, the full derivation for every term will not be demonstrated, but it can be easily verified by the reader.

where

δLP ¼ −



δuP ¼ ðLP Þ−1

We introduce J1 as the first term of equation 7

ϵ1 2

ZZZ

F2λ ðλ; τÞI 2PS ðλ; τÞdλdτdx;

(A-1)

and begin by perturbing equation A-1 with respect to our P- and S-wavefields, uP and uS , to obtain

ZZZ

δJ 1 ¼ ϵ1

R1 ðx; λ; τÞδI PS ðx; λ; τÞdλdτdx;

(A-2)

 2 2u δV ∂ P t P ; V 3P

δI PS ðx;λ;τÞ ¼

0

δJP1 ¼ ϵ1 ⨌ R1 ðx; λ; τÞuS ðx þ λ; t þ τÞ    2  −1 2 δV P ðxÞ∂t uP ðx − λ; t − τÞ dλdτdtdx: × ðLP Þ V 3P (A-13)

T

Using the inner product rule, we rearrange the integral in equation A-13 to remove the operator dependence on δV P :

½uS ðxþ λ;tþ τÞδuP ðx − λ;t− τÞ

þ uP ðx− λ;t− τÞδuS ðx þ λ;t þ τÞdt:

(A-3)



 uS R1 ;ðLP Þ−1

Substituting equation A-3 into equation A-2 gives

δJ1 ¼ δJP1 þ δJS1 ;

(A-4)

where δJP1 is a function of δuP and δJS1 is a function of δuS

δJ P1 ≔ϵ1 ⨌ R1 ðx;λ;τÞuS ðxþλ;tþτÞδuP ðx−λ;t−τÞdλdτdtdx

and

δJS1 ≔ϵ1 ⨌ R1 ðx;λ;τÞuP ðx−λ;t−τÞδuS ðxþλ;tþτÞdλdτdtdx: (A-6) From here, we only present the derivation for the P-wave model gradient K P from equation A-5. The expression for the S-wave gradient may be derived following the same approach. To find δuP , we rewrite equation 1 using linear operator notation:

(A-7)

where LP is the back-propagation acoustic-wave-equation operator:

1 LP ¼ 2 ∂2t − ∇2 : VP

Z δJP1 ¼

2ϵ1 δV P ðxÞ V 3P

(A-9)

   2 −1 2 ¼ LP ðuS R1 Þ; 3 δV P ∂t uP : VP (A-14)

ZZZ

L−1 P ðuS ðx þ λ; t þ τÞ

× R1 ðx; λ; τÞÞ∂2t uP ðx − λ; t − τÞdλdτdtdx:

(A-15)

We apply a shift in the spatial coordinates to simplify the calculation of equation A-15 (Shen and Symes, 2008)

Z δJP1 ¼

2ϵ1 δV P ðxÞ V 3P

ZZZ

L−1 P ½uS ðx þ 2λ; t þ 2τÞ

× R1 ðx þ λ; λ; τÞ∂2t uP ðx; tÞdλdτdtdx:

(A-16)

This can be written more succinctly as

Z δJP1

¼

2ϵ1 δV P ðxÞ V 3P

Z T

0

aP ðx; T − tÞ∂2t uP ðx; tÞdtdx; (A-17)

where

(A-8)

We perturb equation A-7 with respect to model variable V P to obtain

δLP uP þ LP δuP ¼ 0;

2 δV P ∂2t uP V 3P

Therefore, equation A-13 becomes

(A-5)

LP uP ¼ DP ;

(A-12)

and then further substituting this expression into equation A-5 produces

where R1 ðx; λ; τÞ ¼ F2λ ðx; λ; τÞI PS ðx; λ; τÞ, and

Z

(A-11)

Introducing δLP into equation A-10

Term 1: Extended image

J1 ¼

2 δV P ∂2t : V 3P

LP aP1 ¼ wP1 ;

(A-18)

where LP is the forward-propagation operator and aP are the adjoint wavefields through propagation of the adjoint source wP1

ZZ

wP1

¼2

R1 ðx þ λ; λ; τÞuS ðx þ 2λ; t þ 2τÞdλdτ:

and solve for δuP

(A-19) δuP ¼

−ðLP Þ−1 δLP uP ;

(A-10)

Finally, we arrive at

Microseismic joint location and inversion

KS81

Z δJ1 ¼

Z

x

δV P K P1 dx;

(A-20)

Downloaded 07/06/17 to 172.219.108.51. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

where K P1 is the convolution of the forward- and back-propagated wavefields:

K P1 ¼

ϵ1 V 3P

Z

0 T

aP1 ðx; T − tÞ∂2t uP ðx; tÞdt:

(A-21)

By following the above derivation for K P1, we find the expression for K S1 to be

K S1

ϵ ¼ 13 VS

Z

0 T

aS1 ðx; T − tÞ∂2t uS ðx; tÞdt;

δJ P2

¼ −2ϵ2 αS

ZZ

R1 ðx−λ;λ;τÞuP ðx−2λ;t−2τÞdλdτ:

(A-24)

The derivation of the gradients for the remaining terms in equation 7 follow the steps above, ultimately obtaining the same form as equations A-21 and A-22. Therefore, it is sufficient to show only the derivation of the adjoint sources wji , where j is the component (i.e., P or S) and i is the term number in equation 7. In addition because the remaining terms form pairs where P and S are reversed, we will derive only the even terms as the subsequent derivation is identical. For the following derivations, we omit the explicit dependency on space and time for brevity because these terms are based solely on zero-lag images.

Terms 2 and 3: Cross-penalized autocorrelation images

Z F2SS I 2PP dx:

wP2 ¼ 4RP2 uP ;

(A-31)

wS2 ¼ −4αS RS2 uS :

(A-32)

wP3 ¼ −4αP RP3 uP ;

(A-33)

wS3 ¼ 4RS3 uS ;

(A-34)

Next, we derive the gradients for the zero-lag crosscorrelation image I PS penalized by the autocorrelation images I SS and I PP . We only derive the gradient for term 4 of the objective function

ϵ J4 ¼ 4 2

Z F2SS I 2PS dx:

(A-35)

(A-25)

Z

(A-26)

ðF2SS I PS δI PS þ FSS δFSS I 2PS Þdx:

δJ 4 ¼ ϵ4

(A-36)

Substituting the definitions of δI PS from equation A-3 and δFSS from equation A-28 and separating terms based on δuP and δuS , we obtain

Z δJ P4 Z δJ S4 ¼ ϵ4

¼ ϵ4

Z R4

δuP uS dtdx;

Z R4

Z δuS uP dtdx − 2ϵ4 αS

(A-37) Z

RS4

where

uS δuS dtdx; (A-38)

Z δI PP ¼ 2

(A-30)

Terms 4 and 5: Crosscorrelation images penalized by autocorrelation

Z

ðF2SS I PP δI PP þ FSS δFSS I 2PP Þdx;

uS δuS dtdx;

where RP3 ¼ F2PP T PP I 2SS and RS3 ¼ F2PP I SS .

We again perturb equation A-25 with respect to the state variables to obtain

δJ 2 ¼ ϵ2

Z RS2

Following the pattern from above, we perturb equation A-35 to yield

We now derive the gradients for the cross-penalized autocorrelation images I PP and I SS . As Terms 2 and 3 are identical with P and S swapped, we only examine the former:

ϵ J2 ¼ 2 2

(A-29)

Similarly, we can derive the adjoint sources for term 3 of equation 7 to get

and wS1 is the adjoint source

wS1 ¼2

uP δuP dtdx;

where RP2 ¼ F2SS I PP and RS2 ¼ F2SS T SS I 2PP . Following the steps above, we replace δuP by equation A-12 and δuS by its equivalent, apply the inner product rule (similar to equation A-14) and simplify. As stated above, we obtain a version of equations A-21 and A-22, where

(A-22)

(A-23)

¼ 2ϵ2 Z

δJ S2

where

LS aS1 ¼ wS1 ;

Z RP2

uP δuP dt;

δFSS ¼ −αS FSS T SS δI SS ;

(A-27) (A-28)

where T SS ¼ tanhðαS I SS Þ and δI SS is equivalent to A-27 when replacing uP by uS. We again separate J2 into component functions of δuP and δuS

where R4 ¼ F2SS I PS and RS4 ¼ F2SS T SS I 2PS . After substitution and rearrangement (analogous to the procedure followed in equations A-13 and A-14), we find

wP4 ¼ 2R4 uS ;

(A-39)

wS4 ¼ 2ðR4 uP − 2αS RS4 uS Þ:

(A-40)

Witten and Shragge

KS82

δJ ¼ δJP þ δJ S ;

The adjoint sources for term 5, wP5 and wS5 , can be derived similarly as where

Downloaded 07/06/17 to 172.219.108.51. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

wP5 ¼ 2ðR5 uS − 2αP RP5 uP Þ; wS5

(A-41)

¼ 2R5 uP ;

δV P K P dx Z Z δV P 0 P ¼ a ðx; T − tÞ∂2t uP ðx; tÞdtdx V 3P T

where R5 ¼ F2PP I PS and RP5 ¼ F2PP T PP I 2PS .

Z δJ S ¼

δV S K S dx Z Z δV S 0 S a ðx; T − tÞ∂2t uS ðx; tÞdtdx; ¼ V 3S T

Z

F2PS I 2PP dx

(A-43)

Li ai ¼ wi ¼

Z δJ6 ¼ ϵ6

þ

FPS δFPS I 2PP Þdx;

(A-54)

where aP and aS are the adjoint sources generated by forward propagation of the total adjoint sources wP and wS , i.e.

to get

ðF2PS I PP δI PP

(A-53)

and

Finally, we derive the adjoint sources for the autocorrelation images I PP and I SS penalized by the zero-lag crosscorrelation image I PS . We perturb the functional

ϵ J6 ¼ 6 2

Z

δJ P ¼

(A-42)

Terms 6 and 7: Autocorrelation images penalized by crosscorrelation image

(A-52)

7 X

wij ;

(A-55)

j¼1

(A-44)

where i indicates either the P- or S-wave. The total adjoint sources are summations of the individual adjoint sources, wj s, derived above. The total adjoint sources can be written as

where

δFPS ¼ −αPS FPS T PS δI PS ;

(A-45)

and δI PP and δI PS are defined in equations A-27 and A-3 with λ ¼ 0 m, τ ¼ 0 s, respectively, and T PS ¼ tanhðαPS I PS Þ. Separating equation A-44, we find

Z δJP6 ¼ 2ϵ6

Z RP6

Z uP δuP dtdx − ϵ6 αPS

δJ S6 ¼ −ϵ6 αPS

R6

uS δuP dtdx;

Z R6

uP δuS dtdx;

(A-47)

where R6 ¼ F2PS T PS I 2PP and RP6 ¼ F2PS I PP . After substitutions and rearranging (similar to equations A-12–A-14, we obtain the adjoint sources

wP6 ¼ 2ð2RP6 uP − αPS R6 uS Þ;

(A-48)

wS6 ¼ −2αPS R6 uP :

(A-49)

The adjoint sources for term 7 (J7 ¼

wP ¼ 2ðϵ1

ðϵ7 ∕2Þ∫ I 2SS F2PS dx)

R1 ðx þ λ; λ; τÞuS ðx þ 2λ; t þ 2τÞdλdτ

þ 2uP ½ϵ2 RP2 − ϵ3 αP RP3 − ϵ5 αP RP5 þ ϵ6 RP6  þ uS ½ϵ4 R4 þ ϵ5 R5 − ϵ6 αPS R6 − ϵ7 αPS R7 Þ

Z (A-46)

Z

ZZ

(A-56)

and

ZZ wS ¼ 2ðϵ1

R1 ðx − λ; λ; τÞuP ðx − 2λ; t − 2τÞdλdτ

þ 2uS ½ϵ3 RS3 − ϵ2 αS RS2 − ϵ4 αS RS4 þ ϵ7 RS7  þ uP ½ϵ4 R4 þ ϵ5 R5 − ϵ7 αPS R7 − ϵ6 αPS R6 Þ:

(A-57)

The R functions are defined above as functions of the images, penalty functions, and hyperbolic tangents, which come about from the derivative of the zero-lag penalty functions.

REFERENCES are

wP7 ¼ −2αPS R7 uS ;

(A-50)

wS7 ¼ 2ð2RS7 uS − αPS R7 uP Þ;

(A-51)

where R7 ¼ F2PS T PS I 2SS and RS7 ¼ F2PS I SS .

Total gradient From the terms above, we construct the total perturbation of the objective functional

Anikiev, D., J. Valenta, F. Stank, and L. Eisner, 2014, Joint location and source mechanism inversion of microseismic events: Benchmarking on seismicity induced by hydraulic fracturing: Geophysical Journal International, 198, 249–258, doi: 10.1093/gji/ggu126. Artman, B., I. Podladtchikov, and B. Witten, 2010, Source location using time-reverse imaging: Geophysical Prospecting, 58, 861–873, doi: 10 .1111/j.1365-2478.2010.00911.x. Ay, E., H. S. Kuleli, F. Song, and M. N. Toksz, 2012, Detection and enhancement of microseismic signals with correlation operators: Cross product and master event methods: International Geophysical Conference and Oil & Gas Exhibition, Expanded Abstracts, 1–4. Baker, T., R. Granat, and R. W. Clayton, 2005, Real-time earthquake location using Kirchhoff reconstruction: Bulletin of the Seismological Society of America, 95, 699–707, doi: 10.1785/0120040123. Bardainne, T., and E. Gaucher, 2010, Constrained tomography of realistic velocity models in microseismic monitoring using calibration shots: Geophysical Prospecting, 58, 739–753, doi: 10.1111/j.1365-2478.2010 .00912.x.

Downloaded 07/06/17 to 172.219.108.51. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

Microseismic joint location and inversion Batchelor, A., R. Baria, and K. Hearn, 1983, Monitoring the effects of hydraulic stimulation by microseismic event location: A case study: Presented at the SPE 58th Annual Technical Conference and Exhibition, SPE12109, doi: 10.2118/12109-MS. Billings, S. D., M. S. Sambridge, and B. L. N. Kennett, 1994, Errors in hypocenter location: Picking, model, and magnitude dependence: Bulletin of the Seismological Society of America, 84, 1978–1990. Birkelo, B., K. Cieslik, B. Witten, S. Montgomery, B. Artman, D. Miller, and M. Norton, 2012, High-quality surface microseismic data illuminates fracture treatments: A case study in the Montney: The Leading Edge, 31, 1318–1325, doi: 10.1190/tle31111318.1. Chambers, K., B. D. Dando, G. A. Jones, R. Velasco, and S. A. Wilson, 2014, Moment tensor migration imaging: Geophysical Prospecting, 62, 879–896, doi: 10.1111/1365-2478.12108. Chevant, G., 1974, Identification of function parameters in partial differential equations, in Goodson, R. E., and M. Polis, eds., Identification of parameter distributed systems: American Society of Mechanical Engineers, 31–48. Claerbout, J., 1971, Toward a unified theory of reflector mapping: Geophysics, 36, 467–481, doi: 10.1190/1.1440185. Dankbaar, J. W. M., 1985, Separation of P- and S-waves: Geophysical Prospecting, 33, 970–986, doi: 10.1111/j.1365-2478.1985.tb00792.x. Eisner, L., D. Abbott, W. B. Barker, J. Lakings, and M. P. Thornton, 2008, Noise suppression for detection and location of microseismic events using a matched filter: 78th Annual International Meeting, SEG, Expanded Abstracts, 1431–1435. Eisner, L., P. Duncan, W. Heigl, and W. R. Keller, 2009, Uncertainties in passive seismic monitoring: The Leading Edge, 28, 648–655, doi: 10 .1190/1.3148403. Eisner, L., B. J. Hulsey, P. Duncan, D. Jurick, H. Werner, and W. Keller, 2010, Comparison of surface and borehole locations of induced seismicity: Geophysical Prospecting, 58, 809–820, doi: 10.1111/j.1365-2478.2010.00867.x. Ellsworth, W., 2013, Injection-induced earthquakes: Science, 341, 1225942, doi: 10.1126/science.1225942. Galperin, E., 1984, The polarization method of seismic exploration: First Break, 2, 11–14. Grechka, V., P. Singh, and I. Das, 2011, Estimation of effective anisotropy simultaneously with locations of microseismic events: Geophysics, 76, no. 6, WC143–WC155, doi: 10.1190/geo2010-0409.1. House, L., 1987, Locating microearthquakes induced by hydraulic fracturing in crystalline rock: Geophysical Research Letters, 14, 919–921, doi: 10 .1029/GL014i009p00919. Husen, S., and J. Hardebeck, 2010, Earthquake location accuracy, community online resource for statistical seismicity analysis, http://www.corssa .org, accessed 10 October 2016. Ishii, M., P. M. Shearer, H. Houston, and J. E. Vidale, 2005, Extent, duration and speed of the 2004 Sumatra-Andaman earthquake imaged by the HiNet array: Nature, 435, 933–936. Kao, H., and S.-J. Shan, 2004, The source-scanning algorithm: Mapping the distribution of seismic sources in time and space: Geophysical Journal International, 157, 589–594, doi: 10.1111/j.1365-246X.2004.02276.x. Maxwell, S., 2014, Microseismic imaging of hydraulic fracturing: Improved engineering of unconventional reservoirs: SEG. Maxwell, S. C., and T. I. Urbancic, 2001, The role of passive microseismic monitoring in the instrumented oil field: The Leading Edge, 20, 636–639, doi: 10.1190/1.1439012. McMechan, G. A., 1982, Determination of source parameters by wavefield extrapolation: Geophysical Journal International, 71, 613–628, doi: 10 .1111/j.1365-246X.1982.tb02788.x. Michelini, A., and A. Lomax, 2004, The effect of velocity structure errors on double difference earthquake location: Geophysical Research Letters, 31, L09602. Mulder, W., and A. ten Kroode, 2002, Automatic velocity analysis by differential semblance optimization: Geophysics, 67, 1184–1191, doi: 10.1190/ 1.1500380. Myers, S. C., G. Johannesson, and W. Hanley, 2007, A Bayesian hierarchical method for multiple-event seismic location: Geophysical Journal International, 171, 1049–1063, doi: 10.1111/j.1365-246X.2007.03555.x.

KS83

Nakata, N., and G. Beroza, 2016, Reverse time migration for microseismic sources using the geometric mean as an imaging condition: Geophysics, 81, no. 2, KS51–KS60, doi: 10.1190/geo2015-0278.1. Pavlis, G. L., 1986, Appraising earthquake hypocenter location errors: A complete, practical approach for single-event locations: Bulletin of the Seismological Society of America, 76, 1699–1717. 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, doi: 10.1111/j.1365-246X.2006 .02978.x. Rickett, J., and P. Sava, 2002, Offset and angle-domain common imagepoint gathers for shot-profile migration: Geophysics, 67, 883–889, doi: 10.1190/1.1484531. Rutledge, J. T., W. S. Phillips, A. Roff, J. N. Albright, T. Hamilton-Smith, K. C. Jones, K. Steven, and Kimmich, 1994, Subsurface fracture mapping using microearthquakes detected during primary oil production, Clinton County, Kentucky: Presented at the SPE 69th Annual Technical Conference and Exhibition, SPE-28384–MS, doi: 110.2118/28384-MS. Saenger, E. H., 2011, Time reverse characterization of sources in heterogeneous media: NDT & E International, 44, 751–759, doi: 10.1016/j .ndteint.2011.07.011. Sava, P., and I. Vasconcelos, 2011, Extended imaging conditions for waveequation migration: Geophysical Prospecting, 59, 35–55, doi: 10.1111/j .1365-2478.2010.00888.x. Shabelansky, A. H., A. E. Malcolm, M. C. Fehler, X. Shang, and W. L. Rodi, 2015, Source independent full wavefield converted-phase elastic migration velocity analysis: Geophysical Journal International, 200, 952–966, doi: 10.1093/gji/ggu450. Shen, P., 2008, Wave equation migration velocity analysis by differential semblance optimization: Ph.D. thesis, Rice University. Shen, P., and W. Symes, 2008, Automatic velocity analysis via shot profile migration: Geophysics, 73, no. 5, VE49–VE59, doi: 10.1190/1.2972021. Shragge, J., T. Yang, and P. Sava, 2013, Time-lapse image-domain tomography using adjoint-state methods: Geophysics, 78, no. 4, A29–A33, doi: 10.1190/geo2013-0044.1. Symes, W., 1993, A differential semblance criterion for inversion of multioffset seismic reflection data: Journal of Geophysical Research: Solid Earth, 98, 2061–2073, doi: 10.1029/92JB01304. Symes, W., and J. Carazzone, 1991, Velocity inversion by differential semblance optimization: Geophysics, 56, 654–663, doi: 10.1190/1 .1443082. Tang, Y., and G. Ayeni, 2015, Efficient line search methods for multi-parameter full wave-field inversion: U.S. Patent 20, 150, 323, 689. Tang, Y., and B. Biondi, 2011, Target-oriented wavefield tomography using synthesized Born data: Geophysics, 76, no. 5, WB191–WB207, doi: 10 .1190/geo2010-0383.1. Tarantola, A., 1984, Inversion of seismic reflection data in the acoustic approximation: Geophysics, 49, 1259–1266, doi: 10.1190/1.1441754. Thurber, C., 1992, Hypocenter-velocity structure coupling in local earthquake tomography: Physics of the Earth and Planetary Interiors, 75, 55–62, doi: 10.1016/0031-9201(92)90117-E. Thurber, C. H., and E. R. Engdahl, 2000, Advances in global seismic event location: Springer Netherlands, 3–22. Virieux, J., and S. Operto, 2009, An overview of full-waveform inversion in exploration geophysics: Geophysics, 74, no. 6, WCC1–WCC26, doi: 10 .1190/1.3238367. Wang, H., and T. Alkhalifah, 2016, Microseismic imaging using a sourceindependent full-waveform inversion method: 86th Annual International Meeting, SEG, Expanded Abstracts, 2596–2600. Weiss, R., and J. Shragge, 2013, Solving 3D anisotropic elastic wave equations on parallel GPU devices: Geophysics, 78, no. 2, F7–F15, doi: 10 .1190/geo2012-0063.1. Witten, B., and J. Shragge, 2015, Extended imaging conditions for passive seismic data: Geophysics, 80, no. 6, WC61–WC72, doi: 10.1190/ geo2015-0046.1. Yang, T., and P. Sava, 2015, Image-domain wavefield tomography with extended common-image-point gathers: Geophysical Prospecting, 63, 1086–1096, doi: 10.1111/1365-2478.12204.