Downloaded 09/07/15 to 24.114.102.107. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

GEOPHYSICS, VOL. 80, NO. 6 (NOVEMBER-DECEMBER 2015); P. WC61–WC72, 11 FIGS. 10.1190/GEO2015-0046.1

Extended wave-equation imaging conditions for passive seismic data

Ben Witten1 and Jeffrey Shragge1

algorithms, and we have found that they represent a key step toward verifying and updating elastic velocity models. By extending imaging conditions away from zero lag in time and space, we can better evaluate the focusing of a given event based on a self-consistency principle that direct arrival waves focus at zero lag only when the velocity models were correct. We have determined that given an elastic medium and multicomponent recordings, we can propagate and correlate microseismic Pand S-wavefield modes to compute extended imaging conditions that are sensitive to P- and S-wave velocity perturbations. The extended image volume is robust to sparsely sampled data and high levels of noise and is sensitive to velocity model errors. The extended image volume can be used to determine relative P- and S-wave velocity updates to improve focusing.

ABSTRACT The information obtained from seismic monitoring at injection sites is often limited to seismic event properties (e.g., location, origin time, moment tensor), the accuracy of which greatly depends on the assumed or estimated elastic velocity models. Current traveltime-based methods to calibrate velocity models rely on having distinct arrivals. In surface seismic monitoring, however, the signal-to-noise ratio of arrivals is often too poor for pick-based methods, and migration techniques are needed. Wave-equation migration can image microseismic sources without the need for picking or knowledge of the onset time by reconstructing the wavefield through the model space and applying an imaging condition. We have devised extended imaging conditions for passive seismic wave-equation imaging

lies on an accurate model of the physical parameters, particularly the P- and S-wave velocity fields through which seismic waves propagate. Using incorrect velocity models can lead to significant errors in the location estimates (Billings et al., 1994; Husen and Hardebeck, 2010; Gesret et al., 2015) that may produce an incorrect assessment of the project and over- or undercompensation for potential risks. Most methods to quality control (QC) velocity models use events with known spatial and temporal origin, such as perforation shots. For high signal-to-noise recordings, these procedures match forward-modeled traveltimes to picked arrivals in the observed data, the residuals of which can be used to invert for velocity model updates (Warpinski et al., 2003). For low signal-to-noise data, which is often the case for surface-recorded microseismic data, migration methods are generally used to locate events (McMechan, 1982; Artman et al., 2010; Chambers et al., 2010). If focal locations do not match the known locations, it is common to update the velocity

INTRODUCTION Recent years have seen an increase in the injection of fluids into subsurface geologic intervals. This fluid injection has led, directly or indirectly, to more measurable seismicity around the injection site (Ellsworth, 2013) that, in turn, has led to increased seismic monitoring. One of the goals of acquiring seismic monitoring data is to determine the properties of measurable seismic events, such as location, origin time, magnitude, fracture type, and orientation (Dziewonski et al., 1981). These properties are critical for evaluating project success and for mitigating associated potential risks. Obtaining more accurate locations could provide a more robust assessment of project successes or shortcomings (Grechka and Yaskevich, 2014). Determining the location and orientation of smaller events may help to determine the likelihood of inducing larger earthquakes, in which case the injection program can be altered (Royal Society, 2012). Obtaining reliable earthquake properties re-

Manuscript received by the Editor 26 January 2015; revised manuscript received 23 June 2015; published online 26 August 2015. 1 The University of Western Australia, School of Earth and Environment and School of Physics, Centre for Energy Geoscience, Crawley, Australia. E-mail: [email protected]; [email protected]. © 2015 Society of Exploration Geophysicists. All rights reserved. WC61

Witten and Shragge

Downloaded 09/07/15 to 24.114.102.107. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

WC62

models by applying bulk shifts to velocity values or event locations, adjusting individual layers, or introducing anisotropic parameters (Maxwell, 2014). Velocity models, though, are only calibrated along wave/raypaths from the injection well to the receiver array. Currently, the only process to calibrate away from the known shot locations is first-break traveltime tomography (Bardainne and Gaucher, 2010; Grechka and Yaskevich, 2013), which is often not viable for surface data. Artman et al. (2010) introduce multiple zero-lag passive imaging conditions and note that multiple images of the same event (i.e., separate images generated from the P- and Swave data) imaged at inconsistent locations indicate incorrect velocity profiles. We propose using a new passive extended imaging condition (PeIC) to supplement the self-consistency of multiple zero-lag images to QC velocity models for events with unknown source locations and without the need for picking arrivals. We extend a frequency-domain implementation of zero-lag passive waveequation imaging conditions (Artman et al., 2010). For a single event, we invoke a self-consistency principle: All images should have the maximum amplitude at the same spatial location and zero lag in the extended volume. For incorrect velocity models, the maximum amplitude may occur at differing locations in the zero-lag images and away from zero lag in the extended image volume. We first review active-source wave-equation ICs and then introduce the PeICs. We demonstrate a method for QCing an event location by examining the focal location in multiple complementary images, and we demonstrate the sensitivity of the image volume to velocity errors. We examine three reasonable incorrect velocity scenarios. The first model is a 1D assumption of the true laterally varying velocity profile. The second and third models are the correct P-wave model and “good” and “bad” constant V P ∕V S approximations of the true depth-dependent V P ∕V S . These synthetic examples show the utility of PeICs for QCing common approximations of microseismic velocity models. Finally, we discuss practical considerations, including spatial sampling, noise, efficient implementation, and use of the PeIC volume to improve focusing.

WAVE-EQUATION IMAGING Wave-equation migration is a two-step process of propagating two wavefields through a model and combining the wavefield information through an IC to produce an image volume. In active source imaging, one wavefield is constructed from a source function with known spatial location and temporal characteristics and the other from recorded seismograms. Wave-equation propagation reconstructs the wavefields in the model space in either the time domain as a function of space and time, x ¼ fz; x; yg and t, or in the frequency domain as a function of space and frequency, x and ω. Herein, we use a frequency-domain one-way acoustic propagator, and we will restrict our discussion to ωdomain operations. Because we are primarily interested in reconstructing the direct-arrival wavefields, not multiples or mode conversions, this approach represents an efficient and sufficient approximation. Using one-way wave-equation operators, a ω-domain wavefield W can be downward continued from a given depth z to depth step, z þ Δz, according to the following equation:

W zþΔz ¼ W z eikz Δz ;

(1)

where Δz is the extrapolation interval and the vertical wavenumber kz is given by

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ω2 kz ¼ − k2x − k2y ; v2

(2)

where v is the velocity at the current depth step and kx and ky are the horizontal wavenumbers. Further approximations are necessary when there is laterally varying velocity (i.e., v ¼ vðz; x; yÞ). We use a mixed Fourier spatial domain high-order split-step migration scheme (Shan and Biondi, 2005) to solve equations 1 and 2. The sign of the exponential determines whether the data are propagated causally (−) or acausally (+), which is equivalent to forward or backward in time. In active-source imaging scenarios, source wavefields are propagated causally, whereas receiver wavefields are propagated acausally. Once wavefields have been computed, one may apply a correlation IC to form an image. A standard active-source wave-equation IC is constructed by extracting the zero lag of the crosscorrelation between the source and receiver wavefields and summation over all frequencies (Claerbout, 1985):

I a ðxÞ ¼

X ω

W s ðx; ωÞ W r ðx; ωÞ;

(3)

where I a ðxÞ is the resulting active-source image for a single source, W s and W r are the reconstructed source and receiver wavefields, respectively, and the overbar indicates complex conjugation. Conventional IC computation can be extended beyond zero lag by correlating the wavefields shifted symmetrically about an image point. The shifts can be done along the spatial axes (Rickett and Sava, 2002) and/or the time/frequency axis (Sava and Fomel, 2006). The eIC is defined as

I a ðx; λ; τÞ ¼

X ω

e2iωτ W s ðx − λ; ωÞW r ðx þ λ; ωÞ;

(4)

where λ and τ are the spatial and temporal shifts, respectively (Sava and Vasconcelos, 2011). Note that equation 3 is a special case of equation 4, when λ ¼ 0 and τ ¼ 0. Energy from primary reflections focusing away from zero spatial or temporal lag in the extended image volume indicates inaccurate wavefield reconstruction, which we assume to be solely generated from incorrect migration velocities, although this could be related to other factors, such as incomplete wavefield sampling. Misfocusing in extended image volumes has been effectively used to update velocity models in active-source imaging (Sava and Biondi, 2004; Albertin et al., 2005) and teleseismic imaging (Burdick et al., 2013; Shabelansky et al., 2015).

Passive seismic imaging The passive seismic source imaging problem can be treated similarly to the active experiment, except for a few subtle modifica-

Downloaded 09/07/15 to 24.114.102.107. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

Passive eIC tions. Sources of passively recorded seismic energy, such as earthquakes and induced/triggered microseisms, do not have known source functions. Artman et al. (2010) introduce timedomain ICs for passive seismic data using only receiver wavefield components. Contrary to the active-source case, the two wavefields in passive imaging are formed from the recorded wavefield and the sign in the exponent of equation 1 is the same (+) because both wavefields are propagated acausally. Passive source waveequation ICs require recording of the P- and S-wavefields, which necessitates multicomponent data and models for the Pand S-wave velocity profiles. Single-event passive ICs are computed by extracting the zero lag of the autocorrelation of the P- or S-wavefields:

I PP ðxÞ ¼

I SS ðxÞ ¼

X ω

X ω

W P ðx; ωÞ W P ðx; ωÞ;

W S ðx; ωÞ W S ðx; ωÞ;

(5)

(6)

or the crosscorrelation of the P- and S-wavefields

I PS ðxÞ ¼

X ω

W S ðx; ωÞ W P ðx; ωÞ;

(7)

where I PP ðxÞ, I SS ðxÞ, and I PS ðxÞ are the resulting source images and W P and W S are the reconstructed P- and S-wavefields, respectively. We will refer to the images produced by equations 5–7 as the PP, SS, and PS images, respectively. Because we reconstruct the P- and S-wavefields independently through acoustic propagation, the computational cost is greatly reduced with respect to the fully elastic approach. Equation 7 exploits the fact that the P- and S-wavefields are generated by a single event simultaneously and, therefore, they should collapse to the source location cotemporally. In addition, applying equations 5–7 requires separating the vector wavefield into its constituent P- and S-components. Although there are many ways to achieve (or approximate) wave-mode separation (Dankbaar, 1985; Amundsen and Reitan, 1995; Aki and Richards, 2002), we do not elaborate on these methods herein because they are beyond the scope of this paper. If the P- and S-wavefields are separated in the data domain, they can be propagated independently using one-way wave equations. We recognize that this pseudoacoustic scalar approximation of the true wave propagation phenomenon ignores elastic effects; however, we are predominantly interested in developing kinematic constraints and will disregard these effects. Therefore, the v in equation 2 can refer to either the P- or S-wave velocity. Analogous to equation 4, we introduce a PeIC:

I PS ðx; λ; τÞ ¼

X e2iωτ W S ðx − λ; ωÞW P ðx þ λ; ωÞ: (8) ω

Note that equation 7 is a special case of equation 8, where λ ¼ 0 and τ ¼ 0. Like the active source eICs, the PeICs should have en-

WC63

ergy focused around zero lag when using the correct migration velocities and away from zero lag when the velocities are incorrect. As we show in the following sections, extended image volumes are sensitive to P- and S-wave velocity errors and can identify seismic events mislocated due to incorrect velocity models. One may also compute PeICs for the PP and SS images; however, there is no sensitivity to velocity error due to the fact that autocorrelations always have a maximum at zero lag, although the imaged maximum focus may vary spatially with the velocity model. Thus, we will show that combining the PS PeIC with the PP and SS autocorrelation ICs helps to determine the accuracy of the P- and S-wave velocity models as required by the aforementioned self-consistency principle.

SYNTHETIC EXAMPLES We use several synthetic experiments to validate the extended imaging conditions to both QC the velocity models and to demonstrate their sensitivity to velocity. We start by describing the velocity models and forward-modeling operations used to generate the synthetic data. We then describe the simulated data and the resulting images computed using the correct velocity models. Finally, we show the image-domain sensitivity to velocity errors using the same data, but using three different P- and S-wave velocity models: (1) a 1D approximation, (2) a correct P-wave velocity and a representative constant estimation of the V P ∕V S ratio, and (3) and a correct P-wave velocity and a poor constant V P ∕V S approximation.

Forward modeling This section details the 2D velocity models and the forward-modeling procedure used to generate the synthetic data. Figure 1a shows the P-wave model that we use to generate the synthetic data. We create this model by modifying a slice extracted from the SEG/ EAGE 3D overthrust model (Aminzadeh et al., 1997). The right line graph of Figure 1a is the V P ∕V S ratio that we use to specify the S-wave velocity model. The model has minimum and maximum V P ∕V S ratios of 1.83 and 3.17, respectively. The minimum S-wave velocity is 0.55 km∕s at the surface. Note that V P ∕V S is laterally invariant and does not follow the structure. The velocity model has some complexity, including a fault and local heterogeneities, but it is fairly simple overall. The grid size for the forward modeling is 2.5 m in the x- and zdimensions. We place a single source at ½x; z ¼ ½5; 1.5 km, indicated by the asterisk in Figure 1a. We generate multicomponent seismic data by means of an elastic modeling algorithm (Weiss and Shragge, 2013) using a thrust fault source with a dip of 45° (i.e., nonzero stress components Sxx ¼ 1 and Szz ¼ −1) with a Ricker wavelet of central frequency 15 Hz, time steps of 0.2 ms, and no free surface. The resulting wavefields are recorded at each surface grid location on synthetic vertical and horizontal seismograms. Figure 2a and 2b shows the vertical and horizontal data components, respectively. The recorded wavefield contains numerous arrivals with the first P-wave arrival at t ¼ 0.5 s and the first direct S-wave at t ¼ 1.4 s. All other arrivals are reflections, multiples, and/or mode conversions. There is a clear polarity reversal about the horizontal source location (x ¼ 5 km) associated with the thrust fault source mechanism, which also results in a maximum S-wave

Witten and Shragge

WC64

amplitude that is four times larger than the maximum P-wave amplitude.

Downloaded 09/07/15 to 24.114.102.107. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

Imaging We image the passive seismic data through the models as shown in Figure 1 with 5-m grid spacing using the acoustic frequency domain one-way operators described above. We approximately separate the P- and S-components and independently propagate the wavefields with their respective velocities. We estimate the P-wave component by muting the data around the P-wave arrival on the vertical component. Similarly, we estimate the S-wave component by muting before the S-wave arrival on the horizontal component. Because mode conversions and multiples cannot be accurately reconstructed, these wave modes contribute as coherent noise in our resulting image volumes. Figure 3a–3c shows PS, PP, and SS images as grayscale plots with contours overlain constructed using the correct velocities and applying the zero-lag imaging conditions. The contour lines begin at 23% of the maximum value in the panel and in 15% increments, such that the last contour line delineates 98%. All subsequent grayscale plots with contours use the same relative contour values unless otherwise noted. From Figure 3, we note that the maximum amplitude in each image does not occur at the true location, which we indicate by the crosshairs. In fact, the image is a function of the radiation pattern (Artman et al., 2010). Because we are trying to obtain better source location estimates, we need to remove the radiation pattern and ensure the image focuses at the correct location. There are several ways to account for the radiation pattern including crosscorrelation with a master trace in the data domain (Eisner et al., 2008; Ay et al., 2012) or in the image domain by accounting for the radiation pattern in the IC (Saenger, 2011; Anikiev et al., 2014; Chambers et al., 2014). Discussion of the merits of these methods is again beyond the scope of this article. We have chosen to take the envelope of data followed by a band-pass filter (K. Cieslik, personal communication, 2015) prior to imaging to remove the radiation pattern. Note that Figure 1. (a) P-wave velocity model and V P ∕V S ratio used for 2D elastic forward modeling. The asterisk indicates the source location. (b) P-wave model and V P ∕V S ratio for the 1D velocity approximation by extrapolating the velocities at x ¼ 8.5 km, indicated above by the vertical line, throughout the entire model.

a)

b)

this is an essential procedure that is not necessary in active-source imaging. Figure 2c and 2d shows the estimated P- and S-wave data from Figure 2a and 2b after taking the envelope and band-pass filtering. The processed data serve as the input for all subsequent imaging experiments discussed below. Figure 4a–4c shows the zero-lag image using the processed data in Figure 2c and 2d. When comparing Figures 3a–3c and 4a–4c, it is clear that accounting for the radiation pattern collapses the image to a single focus at the true source location in all zero-lag images, and it is thus an essential step in the full-wavefield source imaging process. We now examine PeIC results for the same processed data and the true velocity models. Applying the PeIC transforms 3D wavefields (z; x; ω) to a 5D hypercube (z; x; λz ; λx ; τ). Figure 4d–4f shows slices through the extended image volume with 81 lags in the hypercube space and time dimensions with spatial lag sampling of 5 m and temporal lags of 0.0125 s. The slices in Figure 4d–4f are the fz ¼ zmax ;x ¼ xmax ;λz ¼ 0;λx ;τg, fz ¼ zmax ;x ¼ xmax ;λz ;λx ;τ ¼ 0g, and fz ¼ zmax ; x ¼ xmax ; λz ; λx ¼ 0; τg planes, respectively, where zmax and xmax are the spatial locations of the maximum focal location in the zero-lag PS image (Figure 4a). For the remainder of the paper, when referring to planes extracted from the PeIC volume, we will use the compact notation fλx ; τg. The omitted dimensions, unless stated otherwise, are extracted at zero lag and the maximum focal location in the respective PS image (e.g., fλx ; τg ¼ fz ¼ zmax ; x ¼ xmax ; λz ¼ 0; λx ; τg). All image volume slices show that the maximum image location is at zero lag in all lag dimensions, which is expected when imaging with the correct P- and S-velocities. Note that the maximum focus locations in all zero-lag images and the extended image are internally consistent.

Velocity sensitivity The PP, SS, and PS PeIC volumes shown in Figure 4 were constructed using the correct P- and S-wave migration velocity models. However, to use the various ICs for the purposes of velocity QC,

Passive eIC

Downloaded 09/07/15 to 24.114.102.107. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

these volumes must be sensitive to velocity errors. Thus, we now investigate the sensitivity of the zero-lag images and passive extended image volumes to velocity model errors. We examine two incorrect velocity scenarios that are not uncommon in microseismic monitoring. 1D approximation Figure 1b shows the synthetic single “well-log” velocity model generated by extracting the P- and S-wave velocities at x ¼ 8.5 km from the correct velocity model (the vertical line in Figure 1) and extrapolating this profile throughout the model domain. Therefore, the P- and S-wave velocities are represented by 1D models. Although the V P ∕V S ratio is correct, the S-wave model is not. The 1D approximation has cumulative vertical traveltime errors of 7% for P-waves and 6% for S-waves from the source location to the surface. Although the well log is 3.5 km away from the source location, this model is still quite accurate. We have exactly extracted the true P- and S-wave velocities from depth to the surface, which is seldom the case because well logs are often limited to the reservoir interval. Figure 5a–5c shows the PS, PP, and SS images, respectively, constructed using the 1D velocity models. As is evident from these images, the various zero-lag images have maximum focal locations that are inconsistent in depth. This indicates that either the P- and/or S-wave velocity models are incorrect. We apply the PeIC around the maximum focal location in Figure 5a. Figure 5d–5f shows the fτ; λx g, fλz ; λx g, and fλz ; τg planes extracted from the 5D PeIC volume. The PeIC slices illustrate that the maximum value is shifted away from zero lag along the spatial and temporal lag dimensions, and the focus is smeared due to the a)

b)

c)

d)

WC65

incorrect imaging velocities. Collectively, this provides a significant amount of evidence of velocity inaccuracy. Constant V P ∕V S We now conduct imaging tests using the correct P-wave but incorrect S-wave velocity models. These scenarios are similar to situations in which one has an accurate P-wave velocity from a reflection seismic survey, but no S-wave velocity information is available. The first scenario has V P ∕V S ¼ 2.5, which is approximately equal to the true rms value of V P ∕V S ¼ 2.47. The P-wavepaths are exactly correct, and the S-wavepaths have a cumulative vertical traveltime error of <1%, although nonzero-offsets suffer from greater error. Figure 6a–6c shows the PS, PP, and SS zerolag images, respectively. The PS and PP images focus at the same correct location with modest additional smearing; however, the SS image focuses 200 m deeper indicating errors in at least one of the velocity models. The PS image focuses at the correct location because a significant amount of the S-wave energy in the near offsets arrive at the source location simultaneously with the P-wave energy. The far offsets coalesce with the near-offset energy deeper to produce the maximum focus observed in the SS panel. Figure 6d–6f shows slices extracted from the extended image volume through the fτ; λx g, fλz ; λx g, and fλz ; τg planes. Although the S-wave vertical traveltime is almost correct, the subvertical wavepaths have greater errors that cause lateral smearing. The shift in the focal location along the τ-axis is most noticeable in the fλz ; τg plane, which indicates that even though the zero-lag PS image (Figure 6a) focuses at the correct location, there is additional information about the velocity inaccuracy in the extended images. Figure 2. Modeled 2D elastic data: (a) vertical and (b) horizontal components recorded at the surface from elastic forward modeling of a thrust fault source at ½z; x ¼ ½1.5; 5 km through the model as shown in the top of Figure 1. Panels (c and d) as in (a and b) but after applying a mask around the direct arrivals and computing the wavefield envelope and band-passing the result.

Witten and Shragge

WC66

Downloaded 09/07/15 to 24.114.102.107. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

We next examine the case in which we use a constant V P ∕V S ¼ 2.2, which corresponds to a cumulative 11% vertical S-wave traveltime error. Figure 7a–7c shows the PS, PP, and SS zero-lag images, respectively. Although the maximum focus in

a)

b)

the PP and SS images is very close the true location, the zerolag PS image is maximum within the deeper “blob” with a slightly smaller local maximum near the true source location. Thus, not knowing the correct source location, one may assume that the veloc-

c)

Figure 3. Zero-lag IC images using the correct velocities and the data shown in Figure 2a and 2b: (a) PS, (b) PP, and (c) SS.

a)

b)

c)

d)

e)

f)

Figure 4. Zero-lag IC images using the correct velocities and the data shown in Figure 2c and 2d: (a) PS, (b) PP, and (c) SS. Slices from the PeIC volume at the maximum focal location from panel (a). Panels (d–f) are fτ; λx g, fλz ; λx g, and fλz ; τg, respectively. Note that the focus is at the correct location and zero lag in all dimensions.

Downloaded 09/07/15 to 24.114.102.107. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

Passive eIC ity models are fairly accurate and the global maximum in the PS image in Figure 7a is related to the constructive interference of multiples and/or mode conversions. Figure 7d–7f presents the fτ; λx g, fλz ; λx g, and fλz ; τg planes, respectively. In these slices, there is significant smearing and the maximum amplitude is shifted away from zero lag. Figure 7g–7i shows fτ; λx g, fλz ; λx g, and fλz ; τg planes, respectively, extracted at fz; xg ¼ f1.45; 5g km, near the maximum in the PP and SS images, and zero lag along the omitted lag dimension. This extended image volume would be examined if it were assumed that the maximum generated within the deeper focal area was due to noise. In these panels, we note that there are large shifts away from zero lag in all dimensions. Because the shifts from zero lag are smaller in Figure 7d–7f than Figure 7g–7i, we infer that the deeper focal area is likely related to poorly focused direct arrivals. From the extended image slices, we see that despite all of the zero-lag images exhibiting local maxima near the true source location, there is no internal consistency, and thus velocity errors still remain.

DISCUSSION Figures 4–7 show that without making any assumptions about the source, one may recover information about the accuracy of the mi-

WC67

gration velocities in the image domain. Moreover, fully generating the zero lag and extended image volumes is automated and requires no arrival picking. However, there are practical considerations that must be taken into account, such as spatial sampling, noise, using the PeICs to improve focusing, and the computational expense of performing this analysis. We discuss each of these aspects in turn.

Spatial sampling We first test the effects of sparse spatial sampling of the wavefield by reducing the number of surface receivers. Figure 8a–8d shows the zero-lag PS images, and Figure 8e–8h shows the respective fλz ; λx g slices from the PeIC volume for spatial sampling of 50, 125, 250, and 500 m using the 1D velocity model. The contours for Figure 8a–8d begin at 15% of maximum within the panel and increment by 20%. The contours for Figure 8e–8h begin at 23% of maximum within the panel and step by 15%. When comparing the top row of Figure 8 to Figure 5a, we note that images are very stable with the strength of the background noise in the image increasing relative to the focal strength as the number of receivers decreases. At 500-m spatial sampling, a strong potential false positive is produced above the true source. The decrease in the signal to noise with fewer stations is expected because the signal-to-noise

a)

b)

c)

d)

e)

f)

Figure 5. Zero-lag IC images using the 1D velocity approximation and the data shown in Figure 2c and 2d: (a) PS, (b) PP, and (c) SS. Slices from the PeIC volume at the maximum focal location from panel (a). Panels (d–f) are fτ; λx g, fλz ; λx g, and fλz ; τg, respectively. Note that the focus is shifted from correct location in all panels and away from zero lag in all extended dimensions.

Witten and Shragge

Downloaded 09/07/15 to 24.114.102.107. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

WC68

ratio is proportional to the square root of the number of measurements. The PeIC volumes, though, remain remarkably consistent at the maximum zero-lag image point because the sampling in this space is unchanged. This shows that the PeIC is useful even with sparse spatial sampling given high signal-to-noise data. We next examine how the zero-lag images and PeIC volumes sampled at 125 m respond to noise in the data.

incrementing by 30% using the data from Figure 9 as input. Figure 10 shows the zero-lag images and fτ; λx g, fλz ; λx g, and fλz ; τg the slices through the PeIC volume. Comparing Figures 5 and 10, we note that the PS and PP zero-lag images are severely degraded due to sampling the data every 125 m and the addition of noise. However, in the PeIC slices, there is still a clear focal area that is consistent with those in Figure 5d–5f with a shift toward negative τ and λz .

Noise To generate a noise data set, we propagate single forces at random orientations, times, and locations along the surface of the model. The source wavelet for the noise is the same as the signal. We then add random Gaussian noise, scale, and add this field to the signal. The generated data have a signal-to-noise ratio of 0.33, which we define by the maximum value in the noise-free traces relative to the maximum value in the noise-only data. Figure 9a and 9b shows the vertical and horizontal data components with a signal-to-noise ratio of 0.33 spatially sampled at 125 m. Figure 9c and 9d shows the data after taking the envelope, band-passing, and windowing out the data of interest. Note that the P-wave arrival is only visible in the very near offsets and picked with sufficient accuracy. Figure 10a–10c shows PS, PP, and SS zero-lag images, respectively, with contours beginning at 15% of the panel maximum and

Focal improvement Even without using the PeIC volume for inversion, this analysis can provide a qualitative measure of the direction in which to update the velocity models to improve focusing. To examine how the velocity models can be updated based on the PeIC volume, we forward model data through a homogeneous medium with V P ∕V S ¼ 2.0 and then image with perturbations to the P- and S-wave velocity models. We use 10% maximum velocity perturbation when imaging. Figure 11 shows the τ lag in which the maximum value in the PeIC volume occurs and at zero lag on the remaining axes for all velocity perturbations. Using λz produces a similar pattern. There is a line in which the maximum value occurs at τ ¼ 0, the slope of which is a function of the true V P ∕V S . As is evident from Figure 11, there is more sensitivity to S- than P-wave velocity. For a given S-wave

a)

b)

c)

d)

e)

f)

Figure 6. Zero-lag IC images using the correct P-wave velocity, a constant V P ∕V S ¼ 2.5, and the data shown in Figure 2c and 2d: (a) PS, (b) PP, and (c) SS. Slices from the PeIC volume at the maximum focal location from panel (a): (d) fτ; λx g, (e) fλz ; λx g, and (f) fλz ; τg. Note that the focus is shifted from correct location and zero lag in the λz and τ dimensions.

Downloaded 09/07/15 to 24.114.102.107. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

Passive eIC velocity perturbation, the sign of the τ lag remains relatively constant for varying P-wave velocity. For a given P-wave velocity perturbation, the sign of the τ-lag changes with S-wave velocities. Based on Figure 11, a general rule can be devised to improve focusing using PeIC volumes. If the PeIC volume has a positive τ shift, the S-wave velocity should be increased relative to the P-wave velocity. Conversely, if there is a negative τ shift, the S-wave velocity should be decreased relative to the P-wave velocity.

WC69

Computational expense One potential drawback is that computing PeICs can lead to very large image volumes and require significant resources. For the 2D example presented here, using an admittedly large 81 lags in all dimensions about each model point leads to a full PeIC volume greater than 1 TB. However, instead of calculating the entire extended image, we only need to produce the volume in the neighborhood of the fo-

a)

b)

c)

d)

e)

f)

g)

h)

i)

Figure 7. Zero-lag IC images using the correct P-wave velocity, a constant V P ∕V S ¼ 2.2, and the data shown in Figure 2a and 2b: (a) PS, (b) PP, and (c) SS. Slices from the PeIC volume at the maximum focal location and fz; xg ¼ f1.45; 5g in panel (a): (d and g) fτ; λx g, (e and h) fλz ; λx g, and (f and i) fλz ; τg.

Witten and Shragge

Downloaded 09/07/15 to 24.114.102.107. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

WC70

a)

b)

c)

d)

e)

f)

g)

h)

Figure 8. (a–d) Zero-lag PS images for surface receiver sampling of 50, 125, 250, and 500 m; (e–h) fλz ; λx g slice from the PeIC volume above. The zero-lag image quality degrades as the station density decreases, with the background noise strength increasing relative to the focus. The PeIC volume remains stable because the sampling in this space is not affected.

Figure 9. Data as in Figure 2 with propagated and random noise added. (a) Vertical and (b) horizontal components recorded at the surface. Panels (c and d) as in panels (a and b) but after applying a mask around the direct arrivals, computing the wavefield envelope and band-passing the result.

a)

b)

c)

d)

Downloaded 09/07/15 to 24.114.102.107. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

Passive eIC

WC71

a)

b)

c)

d)

e)

f)

Figure 10. Zero-lag IC images using the 1D velocity approximation and the noisy data subsampled at 125 m shown in Figure 9c and 9d: (a) PS, (b) PP, and (c) SS. Slices from the PeIC volume at the maximum focal location from panel (a): (d) fτ; λx g, (e) fλz ; λx g, and (f) fλz ; τg. Note that PS and PP images are much more degraded than the corresponding fully sampled noise-free images in Figure 5, whereas the PeIC slices show a similar shift in focal position.

Figure 11. The τ lag in which the maximum PeIC value occurs for P- and S-wave velocity perturbations in a homogeneous medium. The trend clearly demonstrates that a qualitative sense of how to update the relative velocities can be derived from the PeIC volume.

cus. We advocate first generating the zero-lag image and automatically extracting the location of the maximum focus or set of foci (Artman and Witten, 2012). One then calculates the PeIC in a small neighborhood about the zero-lag focal location(s) for a user-defined number of lags similar to active-source common image points (Sava and Vasconcelos, 2011). This greatly reduces the computational time and size of the output volume. In this example, computing the PeIC volume approximately 20 points in the spatial dimension reduces the volume to 3.2 GB. If we only computed the PeIC at a single point (e.g., the maximum focal location) the output volume would be reduced from 5D to 3D and require only 2 MB of memory. Although not shown in this article, we assert that the application to 3D is a trivial algorithmic extension. We note that a wavefield volume Wðz; x; y; ωÞ would expand into a 7D PeIC volume ðz; x; y; λx ; λy ; λz ; τÞ making calculation of the entire extended image volume prohibitively expensive. Similar to the 2D case, we can restrict the number of spatial points about which we calculate lags. If we again take 20 points along the spatial dimensions and 81 samples in the lag dimensions, the 7D volume would be 10 TB. Therefore, it may be necessary to calculate the lags around a single point in the 3D case reducing the data to a 4D volume and to a file size of 165 MB for 81 lags in the ðλx ; λy ; λz ; τÞ dimensions, which is small enough to be held in the computer’s memory.

Witten and Shragge

WC72

Downloaded 09/07/15 to 24.114.102.107. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

CONCLUSIONS We introduce and analyze PeICs based on the P- and S-wave receiver wavefields that are analogous to active source-receiver eICs. The PeICs can be applied to any unknown source without any prior information about the source timing or location assuming multicomponent data. The volumes generated by PeICs are sensitive to velocity model error. Imaging with correct velocities produces a PeIC volume that has the maximum focus at zero lag, whereas incorrect velocities give rise to smearing and a shift of the focus away from zero lag. This information can be used for evaluating the velocity model accuracy in conjunction with PP and SS zero-lag images based on a self-consistency principle. The volumes can be constructed efficiently by performing an initial zero-lag imaging step and then calculating the PeIC about the maximal focal location(s). We assert that the PeIC volumes are a key step toward image-domain tomographic updates from passive seismic data. Because PeIC generation relies on wave-equation migration, the wavefield can be approximated with a sparse number of sensors and low signal-to-noise data making it practical to acquire the necessary data. Beyond seismic monitoring, this method may have applications in other research fields, such as nondestructive testing and photoacoustics.

ACKNOWLEDGMENTS We would like to thank the sponsors of the UWA:RM consortium. We are also grateful for an Australian Society of Exploration Geophysicists Research Foundation grant and a University of Western Australia IPRS scholarship that supported this work. We also thank D. Lumley and other members of the UWA Center of Energy Geoscience for useful discussions. We acknowledge the help of K. Cieslik of Spectraseis, Inc., for helpful discussions on removing radiation patterns in our images. We would also like to thank B. Artman, the three anonymous reviewers, and the special section assistant editor L. Eisner for suggestions on improving this manuscript. This work was supported by resources provided by the Pawsey Supercomputing Centre with funding from the Australian government and the government of Western Australia. The reproducible examples were generated using the Madagascar package (www.ahay.org).

REFERENCES Aki, K., and P. Richards, 2002, Quantitative seismology: University Science Books. Albertin, U., P. Sava, J. Etgen, and M. Maharramov, 2005, Adjoint waveequation velocity analysis: 76th Annual International Meeting, SEG, Expanded Abstracts, 3345–3349. Aminzadeh, F., J. Brac, and T. Kunz, 1997, 3-D salt and overthrust models: SEG. Amundsen, L., and A. Reitan, 1995, Decomposition of multicomponent seafloor data into upgoing and downgoing P- and S-waves: Geophysics, 60, 563–572, doi: 10.1190/1.1443794. 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. Artman, B., and B. Witten, 2012, Wave equation microseismic imaging and event selection in the image domain: 82nd Annual International Meeting, SEG, Expanded Abstracts, 1699–1703. Ay, E., H. S. Kuleli, F. Song, and M. N. Toksöz, 2012, Detection and enhancement of microseismic signals with correlation operators: Cross product and master event methods: Presented at the International Geo-

physical Conference and Oil & Gas Exhibition, Expanded Abstracts, 1–4. 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. 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. Burdick, S., M. V. de Hoop, S. Wang, and R. D. van der Hilst, 2013, Reverse-time migration-based reflection tomography using teleseismic free surface multiples: Geophysical Journal International, 196, 996–1017, doi: 10.1093/gji/ggt428. 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. Chambers, K., J.-M. Kendall, S. Brandsberg-Dahl, and J. Rueda, 2010, Testing the ability of surface arrays to monitor microseismic activity: Geophysical Prospecting, 58, 821–830, doi: 10.1111/j.1365-2478.2010 .00893.x. Claerbout, J., 1985, Imaging the earths interior: Blackwell Scientific. 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. Dziewonski, A. M., T.-A. Chou, and J. H. Woodhouse, 1981, Determination of earthquake source parameters from waveform data for studies of global and regional seismicity: Journal of Geophysical Research: Solid Earth, 86, 2825–2852, doi: 10.1029/JB086iB04p02825. 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. Ellsworth, W. L., 2013, Injection-induced earthquakes: Science, 341, 6142, doi: 10.1126/science.1225942. Gesret, A., N. Desassis, M. Noble, T. Romary, and C. Maisons, 2015, Propagation of the velocity model uncertainties to the seismic event location: Geophysical Journal International, 200, 52–66, doi: 10.1093/gji/ggu374. Grechka, V., and S. Yaskevich, 2013, Inversion of microseismic data for triclinic velocity models: Geophysical Prospecting, 61, 1159–1170, doi: 10.1111/1365-2478.12042. Grechka, V., and S. Yaskevich, 2014, Azimuthal anisotropy in microseismic monitoring: A Bakken case study: Geophysics, 79, no. 1, KS1–KS12, doi: 10.1190/geo2013-0211.1. Husen, S., and J. Hardebeck, 2010, Earthquake location accuracy, community online resource for statistical seismicity analysis, http://www.corssa .org, accessed 10 January 2015. Maxwell, S., 2014, Microseismic imaging of hydraulic fracturing: Improved engineering of unconventional reservoirs: SEG. 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. 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. Royal Society, 2012, Shale gas extraction in the UK: A review of hydraulic fracturing, https://royalsociety.org/policy/projects/shale-gas-extraction/ report, accessed 15 January 2015. 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 B. Biondi, 2004, Wave-equation migration velocity analysis. II. Subsalt imaging examples: Geophysical Prospecting, 52, 607–623, doi: 10.1111/j.1365-2478.2004.00448.x. Sava, P., and S. Fomel, 2006, Time-shift imaging condition in seismic migration: Geophysics, 71, no. 6, S209–S217, doi: 10.1190/1.2338824. 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. Shan, G., and B. Biondi, 2005, 3D wavefield extrapolation in laterally-varying tilted TI media: 75th Annual International Meeting, SEG, Expanded Abstracts, 104–107. Warpinski, N. R., R. B. Sullivan, J. E. Uhl, C. K. Waltman, and S. R. Machovoe, 2003, Improved microseismic fracture mapping using perforation timing measurements for velocity calibration: SPE Journal, 10, 14–23, doi: 10.2118/84488-MS. 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.