# Microseismic source locations with deconvolution migration

Microseismic source locations with deconvolution migration Summary Identifying and locating microseismic events are critical problems in hydraulic fracturing monitoring for unconventional resources exploration. In contrast to active seismic data, microseismic data are usually recorded with unknown source excitation time and source location. In this study, we introduce deconvolution migration by combining deconvolution interferometry with interferometric cross-correlation migration (CCM). This method avoids the need for the source excitation time and enhances both the spatial resolution and robustness by eliminating the square term of the source wavelets from CCM. The proposed algorithm is divided into the following three steps: (1) generate the virtual gathers by deconvolving the master trace with all other traces in the microseismic gather to remove the unknown excitation time; (2) migrate the virtual gather to obtain a single image of the source location and (3) stack all of these images together to get the final estimation image of the source location. We test the proposed method on complex synthetic and field data set from the surface hydraulic fracturing monitoring, and compare the results with those obtained by interferometric CCM. The results demonstrate that the proposed method can obtain a 50 per cent higher spatial resolution image of the source location, and more robust estimation with smaller errors of the localization especially in the presence of velocity model errors. This method is also beneficial for source mechanism inversion and global seismology applications. Image processing, Interferometry, Induced seismicity, Seismic interferometry INTRODUCTION Determining the subsurface seismic source location is important in a variety of geophysical experiments. Identifying and locating the microseismic events are critical for passive monitoring, which can help one understand the efficiency of the treatment and further optimize production in unconventional resources (Maxwell 2014). Unlike in the processing of active seismic data, microseismic data are recorded with an unknown source excitation time and source location, similar to tectonic earthquakes. The arrival-time inversion method and its variations, as common earthquake location algorithms (Lee & Stewart 1981), are typically used in microseismic monitoring to retrieve both the source location and excitation time. These techniques aim to minimize the difference between the theoretical and observed P- and/or S-wave arrival times (Geiger 1912). The double-difference method redefines the residuals between theoretical and observed arrival times for pairs of earthquakes at each receiver and can remove the common errors found in both velocity models and records by incorporating the relative information among events (Deichmann & Garciafernandez 1992; Waldhauser & Ellsworth 2000). The paraxial ray-based relative location is modified for microseismicity (Grechka et al.2015). For these schemes, picking the corresponding arrival times is an inevitable step and can be determined manually or by using automatic algorithms (Schimmel & Paulssen 1997; Withers et al.1998; Wu et al.2016). However, obtaining it remains tedious and time-consuming. Furthermore, the chosen algorithm may struggle to produce reliable pickings for a low signal-to-noise ratio (SNR) wavefield or complicated phases. In contrast to the above inversion-based algorithms, migration-based location methods make direct use of recorded wavefield information without identifying or measuring phase arrivals. Array seismology methods, often called ‘delay-and-stack’ methods, are widely used in global seismology to avoid picking (Rost & Thomas 2002). The source scanning algorithm (Kao & Shan 2007), multiscale matched-field algorithms in frequency (Corciulo et al.2012), combining network and array waveform coherence method (Sick & Joswig 2017) and the source mechanism corrected semblance of amplitudes (Anikiev et al.2014; Stanek et al.2015) are the practical implementations of this method that searches the maximum sum of the absolute amplitudes of the receivers with trial location and source excitation time, similar to the Kirchhoff reconstruction (Baker 2005). Given a reasonably accurate velocity model, the time-reversed extrapolation or backpropagation method can focus the observed wavefields onto the original source location (Fink 1997; Gajewski & Tessmer 2005). This algorithm is complemented by some physically meaningful imaging conditions, such as the zero lag of autocorrelation and cross-correlation (Artman et al.2010) and its extension into the frequency domain (Witten & Shragge 2015) and interferometric imaging condition (Wang et al.2013). The deconvolution can be extended to improve the imaging of a microseismic event in an elastic medium (Douma & Snieder 2015). The source locations and excitation times can be determined by scanning the entire extrapolated wavefields (general 2-/3-D in space plus 1-D time) for general passive data (Nakata & Beroza 2016). Further, there are some methods proposed to avoid the need of the source excitation time, such as the semblance technique, which finds the point in space that maximizes a semblance measure of the arrival of specific phases (Duncan & Eisner 2010); the semblance-weighted deconvolution image condition between two or more reconstructed source signatures (Haldorsen et al.2013); and the geometric-mean Reverse Time Migration method, which cross-correlates the extrapolated wavefields (Nakata & Beroza 2016). However, the storage and calculation costs of the extrapolated wavefields are tremendous (Nakata & Beroza 2016). Seismic interferometry is a technique that is commonly used to reconstruct the impulse response between receivers without knowing any information about the source. The most typical interferometry in seismology is implemented by utilizing the cross-correlation of any pair of records at the receivers (Campillo & Paul 2003; Schuster et al.2004; Schuster 2009). There are many applications of seismic interferometry in a wide range of contexts. The uses of active interferometry for exploration include the transformation and migration of vertical seismic profiling and inverse vertical seismic profiling (Yu & Schuster 2006), noise attenuation (Halliday et al.2010) and seismic interpolation (Wang et al.2009, 2010). Passive seismic interferometry aims to turn passive seismic measurements into deterministic seismic responses (Wapenaar et al.2010). It shows potential application in hydraulic fracturing (Poliannikov et al.2011). Interferometric cross-correlation migration (CCM) was thus proposed (Schuster et al.2004; Yu & Schuster 2006) to determine the source distribution by using an asymptotic migration scheme along differential traveltime curves, also for tremor location (Li et al.2017). The method can eliminate the excitation source time by generating the virtual gather. However, the method will degrade the spatial resolution due to the additional square term of the source wavelets. To enhance the resolution of the microseismic source locations, we introduce deconvolution migration (DCM) by combining deconvolution interferometry (Vasconcelos & Snieder 2008) with conventional interferometric CCM (Schuster et al.2004). First, we outline the DCM methodology for the localization of hydraulic fracturing induced seismicity, demonstrate its features by using synthetic data and field data from hydraulic fracturing project to find the source location, and test the sensitivity the velocity model errors and data distortions. METHOD Conventionally, frequency-domain passive data u(k, s, ω) at receiver k are represented as the multiplication of the Green's function G(k, s, ω) convolved with a source function W and a term associated with excitation time τs of the source s:   $$u\ \left( {k,s,\omega } \right) = \ W\left( \omega \right)G\left( {k,s,\omega } \right){e^{i\omega {\tau _s}}},$$ (1) where ω represents the circular frequency. The microseismic data, such as from the hydraulic fracturing monitoring, are usually recorded without knowing the source excitation time τs. This makes the microseismic data processing much more difficult than processing the active seismic data. Seismic interferometry is a useful method for obtaining subsurface structures or source distribution (Schuster et al.2004). To determine the source distribution, the CCM (Schuster et al.2004; Schuster 2009) is treated as an asymptotic migration scheme based on an interferometric imaging method. The algorithm can eliminate the unknown excitation source time by generating the virtual gather with cross-correlation operator. Thus, the virtual data $${\tilde{\Phi }_{{\rm{corr}}}}( {k,j} )$$ of any two traces recorded at k (as the master trace) and j can be expressed generally as   $${\tilde{\Phi }_{{\rm{corr}}}}\left( {k,j} \right) = {\rm{\ }}u\left( {k,s,\omega } \right){u^*}\left( {j,s,\omega } \right).$$ (2) In the context of classical scattering theory (Rodberg & Thaler 1967; Vasconcelos & Snieder 2008), the Green's function  G(s, ω) is expressed as the superposition of the direct wave G0 and all of the following scattered waves Gs. To locate the buried source, only the cross-correlation term of the direct wavefield is needed (Schuster et al.2004). Because of its shortest traveltime (first arrival) from the buried source to the receivers, direct wave is unambiguous and avoids multipath wave propagation in time domain, which is decoupled with the later scattered waves from the complex wavefield. Consequently, the direct wave G0 can be extracted with a mask function M:   $$D\left( {s,\omega } \right) = Mu\left( {s,\omega } \right) = W\left( \omega \right){G_0}\left( {s,\omega } \right){e^{i\omega {\tau _s}}}.$$ (3) After inserting eq. (3) into the general cross-correlation eq. (2), it is more convenient to describe the physics of CCM for source location (see  Appendix) and the new virtual gather $${\tilde{\Phi }_{{\rm{corr}}}}$$ as   $${\tilde{\Phi }_{{\rm{corr}}}}\left( {k,j} \right) = {\rm{\ }}D\left( k \right){D^*}\left( j \right) = {\left| {W\left( \omega \right)} \right|^2}{G_0}\left( {k,\omega } \right)G_0^*\left( {j,\omega } \right).$$ (4) Here, the term |W(ω)|2is the square of the source wavelets, which is unknown and results in lower resolution images of the source location, as demonstrated in the following examples. In the interferometry literature, most authors suggest eliminating the term |W(ω)|2 (Vasconcelos & Snieder 2008), but independently estimation of the source function is not a viable option in many cases. To eliminate the square term |W(ω)|2 and enhance the resolution of the source image, we introduce the deconvolution type method to generate the virtual gather. This virtual gather can be expressed as $${\tilde{\Phi }_{{\rm{de}}}}( {k,j} )$$ (see  Appendix):   $${\tilde{\Phi }_{{\rm{de}}}}\left( {k,j} \right) = \frac{{u\left( k \right)}}{{u\left( j \right)}} = \frac{{u\left( k \right){u^*}\left( j \right)}}{{{{\left| {u\left( j \right)} \right|}^2}}} \approx \frac{{{G_0}\left( {k,\omega } \right)G_0^*\left( {j,\omega } \right)}}{{|{G_0}\left( {j,\omega } \right){|^2}}}$$ (5) Here, the deconvolution operator is similar to the cross-correlation operator in eq. (3), except for an eliminated term |W(ω)|2 and an additional term |G0(k, ω)|2. The numerical stability is always compromised by the presence of the denominator, especially because the whole wavefield results in strong oscillatory effects due to the cross terms between the direct and scattered wavefields (Bleistein & Handelsman 1975). Consider using only the direct wavefield with the mask function in eq. (3), the denominator term |G0(k, ω)|2 will contain no cross terms between the direct and the scattered wavefields (see  Appendix), which is relatively stable than the method using the whole wavefield (Bleistein & Handelsman 1975; Vasconcelos & Snieder 2008). Further, a dynamic damping term (Clayton & Wiggins 1976) is added in the denominator to improve the stability:   $${\tilde{\Phi }_{{\rm{de}}}}\left( {k,j} \right) \approx \frac{{{G_0}\left( {k,\omega } \right)G_0^*\left( {j,\omega } \right)}}{{|{G_0}\left( {k,\omega } \right){|^2} + \varepsilon \left\langle {|{G_0}\left( {j,\omega } \right){|^2}} \right\rangle }}.$$ (6) Here, ε is the scaling factor, which is a function of the source wavelet. In this study, it is set to 0.00001. For active seismic data, the standard diffraction-stack migration can find the image points by smearing the event energy along an ellipsoidal zone in the model space, which is determined by the total reflection traveltime from the source to the receiver on the image point (reflector).For passive seismic data, the standard diffraction-stack migration algorithm only uses the one way traveltime (τkx or τxj) for the buried source location point as shown in Fig. 1, but needs the excitation time τs. Further, it can eliminate the unknown excitation time by restructuring the virtual gather, which reassembles the original one-wave path gather (from the source to receiver k or j) into a two-wave path gather (traveltime difference τxj − τxk of the master trace k and trace j). Then, the interferometric migration (CCM or DCM) can image the buried source location by smearing the energy of the virtual gather into the model space along an asymptotic zone, which is determined by the differential traveltime (subtraction by the traveltimes from the source to the two receivers). Thus, the migration kernel can be expressed as Figure 1. View largeDownload slide Illustration of interferometric cross-correlation migration. The left-hand picture denotes the localization algorithm with the standard diffraction-stack migration algorithm, which needs the source excitation time. The right-hand picture shows the deconvolution method, which can eliminate the source excitation time by generating the virtual gather, as shown at the top of the image. The red dashed line is the migration ray path for source localization, where the energy in trace smears over the true source location. Figure 1. View largeDownload slide Illustration of interferometric cross-correlation migration. The left-hand picture denotes the localization algorithm with the standard diffraction-stack migration algorithm, which needs the source excitation time. The right-hand picture shows the deconvolution method, which can eliminate the source excitation time by generating the virtual gather, as shown at the top of the image. The red dashed line is the migration ray path for source localization, where the energy in trace smears over the true source location.  $${e^{ - i\omega \left( {{\tau _{xj}} - {\tau _{xk}}} \right)}}.$$ (7) Here, the term (τxj − τxk) means to remove the traveltime difference of the virtual gather. The idea is of utilizing differential traveltimes on pairwise receivers from common source is analogous to double-difference inversion method (Waldhauser & Ellsworth 2000; Zhang et al.2010). Both two types of methods can reduce the effect of the velocity uncertainties along the very similar ray paths near the source region when the two stations are close to each other (Zhang et al.2010). Then, the image can be obtained by multiplying the virtual data $${\tilde{\Phi }_{{\rm{de}}}}( {k,j} )$$ or $${\tilde{\Phi }_{{\rm{corr}}}}( {k,j} )$$ with the migration kernel. Consequently, the possible locations of microseismic sources are discretized at each grid point of the model. In this way, the final migration equation for estimating the buried source locations can be expressed with a common form:   $${m_{{\rm{dem}}}}\left( x \right){\rm{\ }} = \mathop \sum \limits_{k = 1}^N \mathop \sum \limits_{j = 1}^N \mathop \sum \limits_{\omega = 1}^W {\tilde{\Phi }_{de}}\left( {k,j,\omega } \right){e^{ - i\omega \left( {{\tau _{xj}} - {\tau _{xk}}} \right)}}$$ (8) where N and W denote the number of the gather and frequency band, respectively. The innermost summation over frequency leads to the constructive interference of the migrated single time trace in the virtual gather. The outer summation over receiver j by the virtual gather yields a single image and is further summated over receiver k as a stacking of each single image. The equation is very similar to the standard diffraction-stack migration, except for the input data of the virtual gather and the time-difference-based migration kernel. As a processing algorithm for locating the microseismic sources, first, the virtual gathers are generated by correlating the master trace k with all other traces using the deconvolution operator. Then, the single image mdem(x, k) of the source location is obtained by using the time-difference-based migrating equation. Finally, the final migration image mdem(x) is achieved by stacking the single image results. Synthetic example In the first example, the proposed DCM is tested with a synthetic seismic data set. The synthetic velocity model is extracted from the Marmousi model in Fig. 2(a). The model contains inclined stratified formation and a high-velocity interbedded to simulate the shale reservoir. The model has 401 gridpoints in the horizontal direction and 251 gridpoints in the vertical direction. The grid interval is 10 m in both directions. In total, 100 receivers with a 40-m interval are placed along the surface, and a point source is set at the central location (1500, 2000). According to previous research, the events observed on surface networks are usually dominated by amplitudes with peak frequencies below 50 Hz (Eisner et al.2013). And also based on the data statistics in Fig. 24(b), the peak frequencies are between 30 and 60 Hz. Thus, a Ricker wavelet with a peak frequencyof 50 Hz is considered representative of an actual microseismic source. The recording length of the trace is 0.7 s, with a temporal sampling interval of 1 ms. To simulate the unknown excitation source time in practice, the passive gather in Fig. 2(b) is randomly shifted with excitation time of 50 ms, which only contains the direct wave. Figure 2. View largeDownload slide The synthetic velocity model and its microseismic gather. The model (a) is from part of the Marmousi model; the data are generated by a point source located at the red star (1500, 2000). The source excitation time τs is set to 0.05 s, which is considered an unknown parameter. The model contains inclined stratified formation and some high-velocity interbed to simulate the shale reservoir. Figure 2. View largeDownload slide The synthetic velocity model and its microseismic gather. The model (a) is from part of the Marmousi model; the data are generated by a point source located at the red star (1500, 2000). The source excitation time τs is set to 0.05 s, which is considered an unknown parameter. The model contains inclined stratified formation and some high-velocity interbed to simulate the shale reservoir. The interferometric migration methods utilize differential traveltimes in generated virtual gather, which is positively correlated with the offset from the master trace. Therefore, the traces in vicinity of zero offset (master trace) of the virtual gather have very small differential traveltimes. Due to the zero differential traveltimes of the master trace itself, it will achieve migration results of constant value. And the small differential traveltime of the nearest two traces results serious smearing throughout wider ray path in migration image shown in Figs 3(a) and (b). In other words, those traces have poor sensitivities to the event locations in both horizontal and vertical directions, and also degrade the SNR and resolution of the final stacking image of the source location. Therefore, the traces in vicinity of zero offset should be muted to obtain a higher resolution image. However, the near-offset substantially reduces the effect of the velocity uncertainties along the very similar ray paths near the source region. Having more muted traces means that there is less smear, but weaker source enhancement with stacking. The optimal parameter of range of the vicinity depends on the case. The range of nearest two or three traces is recommended for our algorithm. Besides, the results also show that the proposed DCM migration can weakens the smearing significantly. Figure 3. View largeDownload slide Virtual trace selection in vicinity of zero offset. (a) and (b) are the single trace migration images at the one-time receiver interval offset with the master trace at 1.3 km, using the CCM and DCM methods, respectively. The traces in vicinity of zero offset (master trace) of the virtual gather have serious smearing throughout wider ray path in its migration image. The DCM can weaken the smearing significantly. Figure 3. View largeDownload slide Virtual trace selection in vicinity of zero offset. (a) and (b) are the single trace migration images at the one-time receiver interval offset with the master trace at 1.3 km, using the CCM and DCM methods, respectively. The traces in vicinity of zero offset (master trace) of the virtual gather have serious smearing throughout wider ray path in its migration image. The DCM can weaken the smearing significantly. The single migration image of the virtual cross-correlation and deconvolution gathers with the master trace at 1.3 km are shown in Figs 4(a) and (b), respectively. The single migration images of both methods smear over the true source location. Both methods show obvious migration noise of a certain bandwidth that crosses over the source. Figs 5(a) and (b) show the final result of stacking the single migration images of the virtual gathers using the CCM and DCM, respectively. For magnified views and analysis of migration images around the source location shown in Fig. 6, the stacked image can focus energy on the source accurately with distinguishable amplitude and suppress the migration noise efficiently. Fig. 6(d) indicates the difference of these two migration images in the source region. Obviously, the deconvolution method shows higher resolution imaging of the source location than does CCM, especially near the true source location, which approximatively eliminates the square term of the wavelets. For more detail, Fig. 6(c) shows the vertical and horizontal traces crossing at the maximum amplitudes of each image. The energy is well focused by both methods, where the peaks of the amplitude are both in the correct location. Figure 4. View largeDownload slide Single migration images of the source location. (a) and (b) are the single migration images with the master trace at 1.3 km using the cross-correlation and deconvolution operator, respectively. The images generated with both methods smear over the true source location. The deconvolution migration method shows higher resolution than does cross-correlation migration. Figure 4. View largeDownload slide Single migration images of the source location. (a) and (b) are the single migration images with the master trace at 1.3 km using the cross-correlation and deconvolution operator, respectively. The images generated with both methods smear over the true source location. The deconvolution migration method shows higher resolution than does cross-correlation migration. Figure 5. View largeDownload slide Stacked migration images of the source location with true velocity. (a) and (b) are the stacked images using the cross-correlation migration and deconvolution migration obtained by stacking all 100 single migration images. Compared to the single migration image, the stacked image can focus energy on the source accurately with distinguishable amplitude and suppress the migration noise efficiently. Figure 5. View largeDownload slide Stacked migration images of the source location with true velocity. (a) and (b) are the stacked images using the cross-correlation migration and deconvolution migration obtained by stacking all 100 single migration images. Compared to the single migration image, the stacked image can focus energy on the source accurately with distinguishable amplitude and suppress the migration noise efficiently. Figure 6. View largeDownload slide Magnified views and analysis of migration images around the source location in Figs 5(a) and (b) are correspond to the images in Figs 4(a) and (b), respectively. Its vertical and horizontal traces (d) cross at each maximum amplitudes (denoted with white star) extracted from (a) and (b). (d) is the difference of these two migration images in the source. The DCM result offers a sharper image, with approximately twice the resolution in the depth direction and a relatively higher resolution in the horizontal direction. Figure 6. View largeDownload slide Magnified views and analysis of migration images around the source location in Figs 5(a) and (b) are correspond to the images in Figs 4(a) and (b), respectively. Its vertical and horizontal traces (d) cross at each maximum amplitudes (denoted with white star) extracted from (a) and (b). (d) is the difference of these two migration images in the source. The DCM result offers a sharper image, with approximately twice the resolution in the depth direction and a relatively higher resolution in the horizontal direction. Sensitivity study We address the effects of velocity model errors and data distortions for microseismic source location and improvement in spatial resolution. All of these factors are important for applications that use field data. Velocity model is critical for migration-based location algorithms. However, velocity model error is inevitable due to the limitations in velocity estimation. First, we use a smoothed mean velocity model (smoothed with a 100 × 100m2 2-D rectangle filter) for migration, which is more reasonable than using the correct velocity model. The smoothed velocity along the ray paths is nearly correct. Besides, utilizing differential traveltimes on pairwise receivers from common source can reduce the effect of the velocity uncertainties along the very similar ray paths near the source region. Therefore, the energy is well focused at the correct source location with either methods in Fig. 7. Similarly, the DCM method shows a higher resolution image of the source location than the cross-correlation method. Then, we also add fixed ratio ± 5 and ± 10 per cent velocity errors for the smoothed migration velocities shown in Figs 8–11. The DCM shows a smaller shift in the maximum amplitude of the source localization (error) than does the CCM in the presence of velocity errors, where the surface geometry is more sensitive in the vertical direction than the horizontal direction. As the performance statistics in Fig. 12, the result from DCM will be more reasonable in the presence of velocity errors, which can achieve small sources location errors for more correct localization especially in the vertical direction and narrower error bars of the resolution for sharper peak values at the focused source location. Figure 7. View largeDownload slide Performance statistics in the presence of velocity model error. This example uses a smoothed mean velocity model with a 100 × 100 m2 2-D rectangle filter. Figure 7. View largeDownload slide Performance statistics in the presence of velocity model error. This example uses a smoothed mean velocity model with a 100 × 100 m2 2-D rectangle filter. Figure 8. View largeDownload slide Performance statistics in the presence of velocity model error. This example uses a smoothed mean velocity model with a 100 × 100 m2 2-D rectangle filter, and add fixed ratio + 10 per cent velocity errors. Figure 8. View largeDownload slide Performance statistics in the presence of velocity model error. This example uses a smoothed mean velocity model with a 100 × 100 m2 2-D rectangle filter, and add fixed ratio + 10 per cent velocity errors. Figure 9. View largeDownload slide Performance statistics in the presence of velocity model error. This example uses a smoothed mean velocity model with a 100 × 100 m2 2-D rectangle filter, and add fixed ratio + 5 per cent velocity errors. Figure 9. View largeDownload slide Performance statistics in the presence of velocity model error. This example uses a smoothed mean velocity model with a 100 × 100 m2 2-D rectangle filter, and add fixed ratio + 5 per cent velocity errors. Figure 10. View largeDownload slide Performance statistics in the presence of velocity model error. Performance statistics in the presence of velocity model error. This example uses a smoothed mean velocity model with a 100 × 100 m2 2-D rectangle filter, and add fixed ratio − 5 per cent velocity errors. Figure 10. View largeDownload slide Performance statistics in the presence of velocity model error. Performance statistics in the presence of velocity model error. This example uses a smoothed mean velocity model with a 100 × 100 m2 2-D rectangle filter, and add fixed ratio − 5 per cent velocity errors. Figure 11. View largeDownload slide Performance statistics in the presence of velocity model error. Performance statistics in the presence of velocity model error. This is a smoothed mean velocity model with a 100 × 100m2 2-D rectangle filter, and add fixed ratio − 10 per cent velocity errors. Figure 11. View largeDownload slide Performance statistics in the presence of velocity model error. Performance statistics in the presence of velocity model error. This is a smoothed mean velocity model with a 100 × 100m2 2-D rectangle filter, and add fixed ratio − 10 per cent velocity errors. Figure 12. View largeDownload slide Performance statistics in the presence of velocity model error. (a) and (b) are vertical and horizontal traces cross at the maximum amplitudes in Figs 7–10, respectively. The dashed lines indicate the location errors of the maximum, and the error bars show a common confidence interval with 85 per cent of the peak values as the resolution. Both methods show incorrect locations of the focused energy in the presence of velocity errors. However, DCM tends to be more reasonable and correct, because it tends to be closer to the true source position (dashed line) and achieves higher resolution with smaller error bars by reducing the interference from the square term of the source wavelet. Figure 12. View largeDownload slide Performance statistics in the presence of velocity model error. (a) and (b) are vertical and horizontal traces cross at the maximum amplitudes in Figs 7–10, respectively. The dashed lines indicate the location errors of the maximum, and the error bars show a common confidence interval with 85 per cent of the peak values as the resolution. Both methods show incorrect locations of the focused energy in the presence of velocity errors. However, DCM tends to be more reasonable and correct, because it tends to be closer to the true source position (dashed line) and achieves higher resolution with smaller error bars by reducing the interference from the square term of the source wavelet. Second, we use a smoothed random errors velocity model. We create a spatially structures correlated Gaussian error on velocities, and then smooth it with a 100 × 100m2 2-D rectangle filter. We set two-scale distribution interval of Gaussian error with approximately ± 5 and ± 10 per cent of the original true velocities values; and have multiple tests with random combination of the velocity errors. We show one general possible scenario with ± 5 per cent in Fig. 13 and with ± 10 per cent in Fig. 14. The DCM tends to be more robust in the presence of velocity errors, which can achieve smaller sources location errors and higher resolution image of the source location than the CCM method. Theoretically, the error of the traveltimes is always mitigated by the summation of the random distribution velocity errors along the ray paths. And multiple tests and calculated results have clearly shown that the random Gaussian errors model achieve smaller sources location errors than the corresponding fixed ratio velocity errors model. Figure 13. View largeDownload slide Performance statistics in the presence of velocity model error. This example uses a 5 per cent Gaussian error on velocities with a 100 × 100 m2 2-D rectangle filter smoothing. Figure 13. View largeDownload slide Performance statistics in the presence of velocity model error. This example uses a 5 per cent Gaussian error on velocities with a 100 × 100 m2 2-D rectangle filter smoothing. Figure 14. View largeDownload slide Performance statistics in the presence of velocity model error. This example uses a 10 per cent Gaussian error on velocities with a 100 × 100 m2 2-D rectangle filter smoothing. Figure 14. View largeDownload slide Performance statistics in the presence of velocity model error. This example uses a 10 per cent Gaussian error on velocities with a 100 × 100 m2 2-D rectangle filter smoothing. Essentially, near-surface medium is consists of continuous velocity field and density, which has no apparent layer structures and reflectivity interfaces. And it can be considered as some velocity anomalies in a large-scale background field. So, we introduces a smooth random model example to simulate the near-surface velocity anomalies. The model in Fig. 15(a) is created by adding smooth strong random velocity anomalies and removing layered structures. And we generate the synthetic event gather in Fig. 15(b) with the same modeling parameters and the same geometry of the previous complex model. Obviously, the arrival time of the direct wave has perturbation due to the strong random velocity anomalies along the ray paths. Then, we also use a smoothed mean velocity model (smoothed with a 100 × 100 m2 2-D rectangle filter) as the migration velocity. The results in Fig. 16 show that both methods focused energy at the correct source location, but degraded with more migration noise than the previous migration images in Fig. 7. Similarly, the DCM method shows a higher resolution image of the source location than the cross-correlation method. To some extent, it gives a preliminary evidence that the proposed method is effective for complex random velocity model. Figure 15. View largeDownload slide The smooth random velocity model and its microseismic gather. The model (a) has no layered structures and some strong random velocity anomalies. Its synthetic event gather (b) has the same parameters and the same geometry of the previous complex model. The arrival time of the direct wave has obvious perturbation due to the strong random velocity anomalies. Figure 15. View largeDownload slide The smooth random velocity model and its microseismic gather. The model (a) has no layered structures and some strong random velocity anomalies. Its synthetic event gather (b) has the same parameters and the same geometry of the previous complex model. The arrival time of the direct wave has obvious perturbation due to the strong random velocity anomalies. Figure 16. View largeDownload slide Performance with smooth random model. This example uses a smoothed mean velocity model with a 100 × 100 m2 2-D rectangle filter. (a) and (b) are the stacked images using the cross-correlation migration and deconvolution migration. Figure 16. View largeDownload slide Performance with smooth random model. This example uses a smoothed mean velocity model with a 100 × 100 m2 2-D rectangle filter. (a) and (b) are the stacked images using the cross-correlation migration and deconvolution migration. Based on the results, the existence of a coincidence localization of maximum amplitude by DCM and CCM indicates a well-estimated migration velocity as shown in Fig. 6, whereas the shift of the localization by DCM and CCM indicates the presence of velocity model errors as shown in Figs 7–13. Waveform distortions test Noise and distortions are inevitable features of seismograms. We add two types of noise model to the synthetic records, contains a Gaussian random noise and actual noise extracted from our field data. First, we use bandlimited random Gaussian noise to the synthetic records as background noise. The random Gaussian noise and its frequency contents of each trace are shown in Figs 17(a) and (b), respectively. It has the random attribute in time and frequency domain. We add two-level amplitude of the random noise, approximately 20 and 50 per cent of the maximum amplitude of the noise-free data. The results in Figs 18 and 19 show that both methods still provide an amplitude peak around the true source location. In addition, the DCM shows higher resolution of the focused source. However, the deconvolution operator is more sensitive, as evidenced by a significantly larger artefact surrounding the focused source. Its resolution is negatively correlated with the noise level. Figure 17. View largeDownload slide The random Gaussian noise model. (a) and (b) are the random Gaussian noise in time and frequency domains, respectively. Figure 17. View largeDownload slide The random Gaussian noise model. (a) and (b) are the random Gaussian noise in time and frequency domains, respectively. Figure 18. View largeDownload slide Performance statistics in the presence of data distortions. This example is with bandlimited random Gaussian noise of 20 per cent of the maximum amplitude of noise-free data. Figure 18. View largeDownload slide Performance statistics in the presence of data distortions. This example is with bandlimited random Gaussian noise of 20 per cent of the maximum amplitude of noise-free data. Figure 19. View largeDownload slide Performance statistics in the presence of data distortions. This example is with bandlimited random Gaussian noise of 50 per cent of the maximum amplitude of noise-free data. Both methods provide an energy-focused amplitude peak around the true source location with two-level random Gaussian noise. DCM has a relatively higher resolution than does CCM in the presence of data distortions. Figure 19. View largeDownload slide Performance statistics in the presence of data distortions. This example is with bandlimited random Gaussian noise of 50 per cent of the maximum amplitude of noise-free data. Both methods provide an energy-focused amplitude peak around the true source location with two-level random Gaussian noise. DCM has a relatively higher resolution than does CCM in the presence of data distortions. Second, we use an actual noise model as the coherent noise from producing environments. The actual model is extracted and modified from our field data. The actual noise model and its frequency contents of each trace are shown in Figs 20(a) and (b), respectively. It is characterized by a burst of energy in time domain, which has clustered frequency content not random distribution. This kind of noise may originate from traffic activities or producting environments (Birnie et al. 2016). We use the actual noise model with approximately 20 per cent of the maximum amplitude of the noise-free data. The results in Fig. 21 show that both methods still provide an amplitude peak around the true source location. In addition, the DCM shows higher resolution of the focused source. Figure 20. View largeDownload slide The actual noise model. (a) is the actual noise extracted from our field data and (b) its frequency contents of each trace. Figure 20. View largeDownload slide The actual noise model. (a) is the actual noise extracted from our field data and (b) its frequency contents of each trace. Figure 21. View largeDownload slide Performance statistics in the presence of data distortions. This example is with actual noise. Both methods provide an energy-focused amplitude peak around the true source location with two-level random Gaussian noise. DCM has a relatively higher resolution than does CCM in the presence of data distortions. Figure 21. View largeDownload slide Performance statistics in the presence of data distortions. This example is with actual noise. Both methods provide an energy-focused amplitude peak around the true source location with two-level random Gaussian noise. DCM has a relatively higher resolution than does CCM in the presence of data distortions. To a certain extent, the presence of additional noise (both random and coherent noise) will make the deconvolution operator unstable and yield a degraded performance close to that of the CCM. Therefore, a more accurate operator for the denominator will improve the resolution of the deconvolution, which is meaningful for improving the stability and resolution of the proposed method. Field example The method is also tested using a field data set from a hydraulic fracturing project in North China. The geometry of microseismic monitoring is shown in Fig. 22. Six lines of a star-like surface array of geophones are used to gather passive seismic data with recording the vertical particle velocity. The surface topography of the geophones is varied between 300 and 850 m. The number of geophones per arm varied between 237 and 346, with 1780 geophones in total. The average receiver geophone spacing is approximately 20.2 m (varied between 19 and 22 m). The continuous microseismic data is acquired at 1 kHz. One recording segment of contain strong event is selected with the length of 0.5 s. The vertical production well is located near the centre of the array. According to the field record, the burial depth of shale reservoir is about 4000–4200 m, where the majority of events are expected to be located. The velocity model assumes a horizontal-layer medium and is built based on the experiment sonic log at the injection well, as shown in Fig. 23. Figure 22. View largeDownload slide Map view of the surface survey geometry. The circles denote the geophones within six survey lines, and the red line is used for the following migration using a 2-D approximation. The injection well is in the centre of the star-like geometry. The shale reservoir is at a depth of 4000–4200 m. Figure 22. View largeDownload slide Map view of the surface survey geometry. The circles denote the geophones within six survey lines, and the red line is used for the following migration using a 2-D approximation. The injection well is in the centre of the star-like geometry. The shale reservoir is at a depth of 4000–4200 m. Figure 23. View largeDownload slide The built velocity model. The velocity model assumes a horizontal-layer medium and is built based on the experiment sonic log at the injection well. Figure 23. View largeDownload slide The built velocity model. The velocity model assumes a horizontal-layer medium and is built based on the experiment sonic log at the injection well. Without loss of generality, one line with 294 geophones is selected in this study (the red line shown in Fig. 22), which has an approximate average length of the all arm lines. The selected raw data have low SNR shown in Fig. 24(a), which is contaminated with noise and distortions. However, the direct wave is relative unambiguous, which is decoupled with the later complex scattered wavefield. As discussed in the previous section, the data distortions will make the deconvolution operator unstable and performance with degraded resolution. The pre-processing is the most vital step to find the source location accurately. First, the direct wave is isolated from the wavefield with a mask function as the abovementioned examples. Then, a bandpass filter is applied for mainly removing the noise in the site acquisition. The bandwidth parameters are chosen based on the frequency contents as shown in Fig. 24(b), where the most of the signal is concentrated between 30 and 60 Hz. The final processed record is shown in Fig. 24(c). The SNR is effectively improved. Figure 24. View largeDownload slide Data pre-processing. (a) is the full recorded wavefield and the extracted direct wavefield. (b) is its frequency contents of each trace. (c) is the final processed record.The signal-to-noise ratio (SNR) is effectively improved. Figure 24. View largeDownload slide Data pre-processing. (a) is the full recorded wavefield and the extracted direct wavefield. (b) is its frequency contents of each trace. (c) is the final processed record.The signal-to-noise ratio (SNR) is effectively improved. The single migration images generated using the cross-correlated virtual gather and deconvolved virtual gather with the master trace at 2.2 km are shown in Figs 25(a) and (b), respectively. The single migration image shows a low-resolution image with a strong band-like migration ray path crossing the source. The final results with the two methods are obtained by stacking all 294 migration images of the virtual gathers, shown in Fig. 26. Note that the two final stacked images can enhance the source location with a distinguishable focused amplitude and suppress the migration noise efficiently. As the indication in the vertical and horizontal traces in Fig. 27(c), the deconvolution type method performs a higher resolution image than does the cross-correlation type method, especially in the depth direction. The migration image reveals the location to be at (2820, 4090) by CCM and (2820, 4100) by DCM, both of which are in the shale reservoir of 4000–4200 m. Based on the performance statistics sensitivity study, the existence of shift localization indicates the presence of velocity model errors. Actually, an accurate estimation of the field velocity model from sonic log is very difficult to determine by the current methods, especially for the strong attenuation and anisotropy of the complex near-surface. However, the shift localization is relatively small, which means that using the estimated velocity model is acceptable and reasonable for the migration scheme. Further from the sensitivity study, we know that the localization by DCM at (2820, 4100) should be more correct and closer to the true source location, though the true source location for field data is unknown. Figure 25. View largeDownload slide Single migration image comparison. (a) and (b) are the cross-correlation and deconvolution gathers with a master trace at 2.2 km, respectively. Both images smear over the true source location with relatively strong migration noise. Figure 25. View largeDownload slide Single migration image comparison. (a) and (b) are the cross-correlation and deconvolution gathers with a master trace at 2.2 km, respectively. Both images smear over the true source location with relatively strong migration noise. Figure 26. View largeDownload slide Stacked migration images of the source location. (a) and (b) are the stacked images using cross-correlation and deconvolution migration obtained by stacking all 294 single migration images. Both methods are better focused at the location point than are the single migration images. Figure 26. View largeDownload slide Stacked migration images of the source location. (a) and (b) are the stacked images using cross-correlation and deconvolution migration obtained by stacking all 294 single migration images. Both methods are better focused at the location point than are the single migration images. Figure 27. View largeDownload slide Magnified views of migration images around the source location. (a) and (b) correspond to the images in Figs 18(a) and (b), respectively. (c) shows the vertical and horizontal traces crossing at the maximum amplitudes extracted from (a) and (b). The DCM shows a sharper peak value than the CCM. Figure 27. View largeDownload slide Magnified views of migration images around the source location. (a) and (b) correspond to the images in Figs 18(a) and (b), respectively. (c) shows the vertical and horizontal traces crossing at the maximum amplitudes extracted from (a) and (b). The DCM shows a sharper peak value than the CCM. DISCUSSION The DCM method provides reasonable microseismic source localization results using both synthetic data and a field data from hydraulic fracturing monitoring. The results show more robust and higher spatial resolution images of the source localization than does CCM. For computational efficiency, DCM is very similar to the CCM in the processing algorithm, except for the operator that generates the virtual gather. When generating the virtual gather, the proposed method uses the same method to calculate the cross-correlation term D(k)D*(j) and the additional numerator |D(j)|2 in the deconvolution operator in frequency, which requires the same calculation. In addition, the deconvolution operator requires the vector division of complex numbers by the length of the numerator. Thus, the deconvolution operator costs more than twice as much as the cross-correlation operator overall. For the migration scheme, both methods have the same complexity of calculation, which is based primarily on the number of the traces and the size of the model. In addition, both methods share a traveltime table and pre-processing in the study. Although the proposed method doubles the calculation time of the cross-correlation method, it is still far less computationally intensive than are time-reversed extrapolation-based methods. Because every single migration permits a high degree of parallelism with different grains, some Graphics Processing Unit acceleration and multithreading may be good choices for hardware acceleration. For the performance of the source localization, the star-like geometry in hydraulic fracturing is considerably more irregular and sparse than the dense and regular acquisition geometry in active seismic exploration. Consequently, it results in an uneven distribution of wave paths, where the region near the centre of star-like geometry will provide high fold number and get better image illumination. Analogous to standard diffraction-stack method, the kernel of the DCM needs a sufficient number of folds to suppress the migration noise. Therefore, the 2-D approximation is used for the proposed method in this study. However, it is still meaningful to extend the theory and its performance of the DCM to 3-D with the irregular and sparse star-like geometry. In this study, the current DCM method only deals with direct wave. The following the scattering wavefield will provide more information and illustrations for the hydraulic fracturing of the reservoir (Reshetnikov et al.2015). Furthermore, locating a single event is not equivalent to locating a hydraulic fracturing. Applications to 3-D data sets and extension to the scattering wavefield remain future research topics. CONCLUSIONS We present a high-resolution microseismic sources localization method. The DCM method combines deconvolution interferometry with interferometric CCM. It can eliminate the unknown excitation source time by generating the virtual gather and can enhance the spatial resolution and robustness of the source localization by removing the square term of the source wavelet in CCM. The results of applying our method to synthetic and field cases prove the validity of the proposed method. Compared to the conventional CCM, the DCM can enhance the spatial resolution of source location with a sharper value and improve the robustness with a smaller error of the source localization especially in the presence of velocity errors. The method can provide valuable information for further source mechanism inversion and can be beneficial for global seismology applications. Acknowledgements The research was funded by National Natural Science Foundation of China (grant no. 41230317) and Strategic Priority Research Programme (grant no. XDB10030500). We also acknowledge the constructive suggestions and manuscript improvement of the editor Herve Chauris, reviewer Pierre-Francios Roux and an anonymous reviewer. REFERENCES Anikiev D., Valenta J., Staněk F., Eisner L., 2014. Joint location and source mechanism inversion of microseismic events: benchmarking on seismicity induced by hydraulic fracturing, Geophys. J. Int. , 198( 1), 249– 258. https://doi.org/10.1093/gji/ggu126 Google Scholar CrossRef Search ADS   Artman B., Podladtchikov I., Witten B., 2010. Source location using time-reverse imaging, Geophys. Prospect. , 58( 5), 861– 873. https://doi.org/10.1111/j.1365-2478.2010.00911.x Google Scholar CrossRef Search ADS   Baker T., 2005. Real-time earthquake location using Kirchhoff reconstruction, Bull. seism. Soc. Am. , 95( 2), 699– 707. https://doi.org/10.1785/0120040123 Google Scholar CrossRef Search ADS   Birnie C., Chambers K., Angus D.A., Stork L., 2016. Analysis and models of pre-injection surface seismic array noise recorded at the Aquistore carbon storage site, Geophys. J. Int. , 206( 2), 1246– 1260. https://doi.org/10.1093/gji/ggw203 Google Scholar CrossRef Search ADS   Bleistein N., Handelsman R.A., 1975. Asymptotic Expansions of Integrals , Dover Publishing, Inc. Campillo M., Paul A., 2003. Long-range correlations in the diffuse seismic coda. Science , 299( 5606), 547– 549. https://doi.org/10.1126/science.1078551 Google Scholar CrossRef Search ADS PubMed  Corciulo M., Roux P., Campillo M., Dubucq D., Kuperman W.A., 2012. Multiscale matched-field processing for noise-source localization in exploration geophysics, Geophysics , 77( 5), KS33– KS41. https://doi.org/10.1190/geo2011-0438.1 Google Scholar CrossRef Search ADS   Clayton R. W., Wiggins R.A., 1976, Source shape estimation and deconvolution of teleseismic body waves, Geophys. J. Int. , 47( 1), 151– 177. https://doi.org/10.1111/j.1365-246X.1976.tb01267.x Google Scholar CrossRef Search ADS   Deichmann N., Garcia-Fernandez M., 1992. Rupture geometry from high-precision relative hypocenter locations of microearthquake clusters, Geophys. J. Int. , 110, 501– 517. Google Scholar CrossRef Search ADS   Douma J., Snieder R., 2014. Focusing of elastic waves for microseismic imaging, Geophys. J. Int. , 200( 1), 390– 401. https://doi.org/10.1093/gji/ggu398 Google Scholar CrossRef Search ADS   Duncan P., Eisner L., 2010. Reservoir characterization using surface microseismic monitoring, Geophysics , 75( 5), 139– 146. https://doi.org/10.1190/1.3467760 Google Scholar CrossRef Search ADS   Eisner L., Gei D., Hallo M., Oprs̆al I., Ali M.Y., 2013. The peak frequency of direct waves for microseismic events, Geophysics , 78( 6), 45– 49. https://doi.org/10.1190/geo2013-0197.1 Google Scholar CrossRef Search ADS   Fink M., 1997. Time reversed acoustics, Phys. Today , 50( 3), 34– 40. https://doi.org/10.1063/1.881692 Google Scholar CrossRef Search ADS   Gajewski D., Tessmer E., 2005. Reverse modelling for seismic event characterization, Geophys. J. Int. , 163( 1), 276– 284. https://doi.org/10.1111/j.1365-246X.2005.02732.x Google Scholar CrossRef Search ADS   Geiger L., 1912. Probability method for the determination of earthquake epicentres from the arrival time only (translated from Geiger's 1910 German article), Bull. St. Louis Univ. , 60– 71. Grechka V., Li Z., Howell B., 2015. Relative location of microseismicity, Geophysics , 80( 6), 149– 158. https://doi.org/10.1190/geo2014-0617.1 Haldorsen J.B.U., Brooks N.J., Milenkovic M., 2013. Locating microseismic sources using migration-based deconvolution, Geophysics , 78( 5), KS73– KS84. https://doi.org/10.1190/geo2013-0086.1 Google Scholar CrossRef Search ADS   Halliday D.F., Curtis A., Vermeer P., Strobbia C., Glushchenko A., van Manen D.-J., Robertsson O.A.J., 2010. Interferometric ground-roll removal: attenuation of scattered surface waves in single-sensor data, Geophysics , 75( 2), SA15– SA25. https://doi.org/10.1190/1.3360948 Google Scholar CrossRef Search ADS   Lee W.H.K., Stewart W.S., 1981. Principles and Applications of Microearthquake Networks , Academic Press, Inc. Li K.L., Sgattoni G., Sadeghisorkhani H., Gudmundsson O., 2017. A double-correlation tremor-location method, Geophys. J. Int. , 208( 2), 1231– 1236. https://doi.org/10.1093/gji/ggw453 Google Scholar CrossRef Search ADS   Kao H., Shan S.J., 2007. Rapid identification of earthquake rupture plane using Source-Scanning Algorithm, Geophys. J. Int. , 168( 3), 1011– 1020. https://doi.org/10.1111/j.1365-246X.2006.03271.x Google Scholar CrossRef Search ADS   Maxwell S.C., 2014. Microseismic imaging of hydraulic fracturing: improved engineering of unconventional shale reservoirs, in SEG Distinguished Instructor Short Course , doi:10.1190/1.9781560803164. Nakata N., Beroza G.C., 2016. Reverse time migration for microseismic sources using the geometric mean as an imaging condition, Geophysics , 81( 2), 51– 60. https://doi.org/10.1190/geo2015-0278.1 Google Scholar CrossRef Search ADS   Poliannikov O.V., Malcolm A.E., Djikpesse H., Prange M., 2011. Interferometric hydrofracture microseism localization using neighboring fracture. Geophysics , 76( 6), 27– 36. https://doi.org/10.1190/geo2010-0325.1 Google Scholar CrossRef Search ADS   Reshetnikov A., Kummerow J., Asanuma H., Häring M., Shapiro S.A., 2015. Microseismic reflection imaging and its application to the Basel geothermal reservoir, Geophysics , 80( 6), WC39– WC49. https://doi.org/10.1190/geo2014-0593.1 Google Scholar CrossRef Search ADS   Rodberg L.S., Thaler R.M., 1967. Introduction to the Quantum Theory of Scattering , Academic Press, Inc. Rost S., Thomas C., 2002. Array seismology: methods and applications, Rev. Geophys. , 40( 3), 2-1–2-27. https://doi.org/10.1029/2000RG000100 Google Scholar CrossRef Search ADS   Schuster G.T., Yu J., Sheng J., Rickett J., 2004. Interferometric/daylight seismic imaging, Geophys. J. Int. , 157( 2): 838– 852. https://doi.org/10.1111/j.1365-246X.2004.02251.x Google Scholar CrossRef Search ADS   Schuster G.T., 2009. Seismic Interferometry, Cambridge Univ. Press. Schimmel M., Paulssen H., 1997. Noise reduction and detection of weak, coherent signals through phase-weighted stacks, Geophys. J. Int. , 130( 2), 497– 505. https://doi.org/10.1111/j.1365-246X.1997.tb05664.x Google Scholar CrossRef Search ADS   Sick B, Joswig M., 2017. Combining network and array waveform coherence for automatic location: examples from induced seismicity monitoring. Geophys. J. Int. , 208( 3), 1373– 1388. https://doi.org/10.1093/gji/ggw468 Google Scholar CrossRef Search ADS   Stanek F., Anikiev D., Valenta J., Eisner J., 2015. Semblance for microseismic event detection. Geophys. J. Int. , 201( 3), 1362– 1369. https://doi.org/10.1093/gji/ggv070 Google Scholar CrossRef Search ADS   Wang C., Cheng J., Yin C., Hong L., 2013. Microseismic event location using reverse-time focusing and interferometric techniques for surface and borehole observations[J], Chin. J. Geophys. , 56( 56), 584– 597. Wang Y., Luo Y., Schuster G., 2009. Interferometric interpolation of missing seismic data, Geophysics , 74( 3), SI37– SI45. https://doi.org/10.1190/1.3110072 Google Scholar CrossRef Search ADS   Wang Y., Dong S., Luo Y., 2010. Model-based interferometric interpolation method, Geophysics , 75( 6), WB211– WB217. https://doi.org/10.1190/1.3505816 Google Scholar CrossRef Search ADS   Waldhauser F., Ellsworth W.L., 2000. A double-difference earthquake location algorithm: method and application to the Northern Hayward Fault, California, Bull. seism. Soc. Am. , 90( 6), 1353– 1368. https://doi.org/10.1785/0120000006 Google Scholar CrossRef Search ADS   Wapenaar K., Draganov D., Snieder R., Campman X., Verdel A., 2010. Tutorial on seismic interferometry: Part 1—Basic principles and applications, Geophysics , 75( 5), 75A195– 75A209. Google Scholar CrossRef Search ADS   Withers M., Aster R., Young C., Beiriger J., Harris M., Moore S., Trujillo J., 1998. A comparison of select trigger algorithms for automated global seismic phase and event detection, Bull. seism. Soc. Am. , 88( 1), 95– 106. Witten B., Shragge J., 2015. Extended wave-equation imaging conditions for passive seismic data, Geophysics , 80( 6), WC61– WC72. https://doi.org/10.1190/GEO2015-0046.1 Google Scholar CrossRef Search ADS   Wu S., Wang Y., Chang X., 2016. Automatic microseismic event detection by band-limited phase-only correlation, Phys. Earth planet. Inter. , 261, 3– 16. https://doi.org/10.1016/j.pepi.2016.09.005 Google Scholar CrossRef Search ADS   Vasconcelos I., Snieder R., 2008. Interferometry by deconvolution: Part 2—theory for elastic waves and application to drill-bit seismic imaging, Geophysics , 73( 3), S129– S141. https://doi.org/10.1190/1.2904985 Google Scholar CrossRef Search ADS   Yu J., Schuster G.T., 2006. Crosscorrelogram migration of inverse vertical seismic profile data, Geophysics , 71( 1), S1– S11. https://doi.org/10.1190/1.2159056 Google Scholar CrossRef Search ADS   Zhang H., Nadeau R.M., Toksoz M.N., 2010. Locating nonvolcanic tremors beneath the San Andreas Fault using a station-pair double-difference location method, Geophys. Res. Lett. , 37( 13), doi:10.1029/2010GL043577. https://doi.org/10.1029/2010GL043577 APPENDIX: DECONVOLUTION INTERFEROMETRY For general cross-correlation, we can expand $${\tilde{\Phi }_{{\rm{corr}}}}( {k,j} )$$ into four terms by inserting eq. (1) into (2) with direct wave G0 and scattered waves Gs:   \begin{eqnarray} && {\tilde {\Phi }_{{\rm{corr}}}}\left( {k, j} \right) = {G_0}\left( k \right)G_0^*\left( j \right) + {G_s}\left( k \right)G_0^*\left(j \right)\\ && + \, {G_0}\left( k \right)G_s^*\left( j \right) + {G_s}\left( k \right)G_s^*\left( j \right). \end{eqnarray} (A1) Each of the four terms on the right-hand side has a different physical meaning (Schuster et al.2004; Vasconcelos & Snieder 2008). For imaging purposes, we want to use only the Gs terms, which are the dominant contribution to the causal scattered wavefield that comes from the correlation between the direct wavefield at G0 and the scattered wavefield Gs. For source location imaging, only the first cross-correlation term of the direct wavefield is needed (Schuster et al.2004). Consequently, the latter three terms in eq. (A1) are dropped by using the muting processing in eq. (3). When considering a full waveform u(i), the deconvolution operator of eq. (5) can be expressed as   $${\tilde{\Phi }_{{\rm{de}}}}\left( {k,j} \right) = \frac{{u\left( k \right){u^*}\left( j \right)}}{{{{\left| {u\left( j \right)} \right|}^2}}}\ = \frac{{{{\tilde{\Phi }}_{{\rm{corr}}}}\left( {k,j} \right)}}{{{{\left| {\tilde{u}\left( j \right)} \right|}^2}}}$$ (A2) where the relationship between deconvolution $${\tilde{\Phi }_{{\rm{de}}}}( {k,j} )$$ and cross-correlation $${\tilde{\Phi }_{{\rm{corr}}}}( {k,j} )$$ can be explicitly identified. The denominator can also be expanded into four terms, as eq. (A1):   \begin{eqnarray} {{\left | {\tilde{u}\left( j \right)} \right|^2} }\\ && = \frac{1}{{{G_0}\left( j \right)G_0^*\left( j \right) \!+\! {G_0}\left( j \right)G_s^*\left( j \right) \!+\! {G_s}\left( j \right)G_0^*\left( j \right) \!+\! {G_s}\left( j \right)G_s^*\left( j \right)}}\\ \end{eqnarray} (A3) Assuming that the terms of the direct wavefield G0 make the most prominent contributions and the wavefield perturbations Gs is small (|Gs|2 ≪ |G0|2), the last term in the denominator can be dropped (Vasconcelos & Snieder 2008):   $${\left| {\tilde{u}\left( j \right)} \right|^{ - 2}} \approx \frac{1}{{{{\left| {{G_0}\left( j \right)} \right|}^2}\left[ {1 + \frac{{{G_0}\left( j \right)G_s^*\left( j \right)}}{{{{\left| {{G_0}\left( j \right)} \right|}^2}}} + \frac{{{G_s}\left( j \right)G_0^*\left( j \right)}}{{{{\left| {{G_0}\left( j \right)} \right|}^2}}}} \right]}}$$ (A4) and the latter two terms can be ignored, because of the muting of the scatter wavefield Gs(s) in eq. (2), which also satisfies the condition of |Gs|2 ≪ |G0|2. Combining eqs (A1) and (A4), we can recover the expression of the deconvolution operator in eq. (5). All of the cross terms between direct and scattered wavefields in the denominator and numerator are dropped. © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# Microseismic source locations with deconvolution migration

, Volume 212 (3) – Mar 1, 2018
28 pages

/lp/ou_press/microseismic-source-locations-with-deconvolution-migration-SUsMb6oicY
Publisher
The Royal Astronomical Society
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx518
Publisher site
See Article on Publisher Site

### Abstract

Summary Identifying and locating microseismic events are critical problems in hydraulic fracturing monitoring for unconventional resources exploration. In contrast to active seismic data, microseismic data are usually recorded with unknown source excitation time and source location. In this study, we introduce deconvolution migration by combining deconvolution interferometry with interferometric cross-correlation migration (CCM). This method avoids the need for the source excitation time and enhances both the spatial resolution and robustness by eliminating the square term of the source wavelets from CCM. The proposed algorithm is divided into the following three steps: (1) generate the virtual gathers by deconvolving the master trace with all other traces in the microseismic gather to remove the unknown excitation time; (2) migrate the virtual gather to obtain a single image of the source location and (3) stack all of these images together to get the final estimation image of the source location. We test the proposed method on complex synthetic and field data set from the surface hydraulic fracturing monitoring, and compare the results with those obtained by interferometric CCM. The results demonstrate that the proposed method can obtain a 50 per cent higher spatial resolution image of the source location, and more robust estimation with smaller errors of the localization especially in the presence of velocity model errors. This method is also beneficial for source mechanism inversion and global seismology applications. Image processing, Interferometry, Induced seismicity, Seismic interferometry INTRODUCTION Determining the subsurface seismic source location is important in a variety of geophysical experiments. Identifying and locating the microseismic events are critical for passive monitoring, which can help one understand the efficiency of the treatment and further optimize production in unconventional resources (Maxwell 2014). Unlike in the processing of active seismic data, microseismic data are recorded with an unknown source excitation time and source location, similar to tectonic earthquakes. The arrival-time inversion method and its variations, as common earthquake location algorithms (Lee & Stewart 1981), are typically used in microseismic monitoring to retrieve both the source location and excitation time. These techniques aim to minimize the difference between the theoretical and observed P- and/or S-wave arrival times (Geiger 1912). The double-difference method redefines the residuals between theoretical and observed arrival times for pairs of earthquakes at each receiver and can remove the common errors found in both velocity models and records by incorporating the relative information among events (Deichmann & Garciafernandez 1992; Waldhauser & Ellsworth 2000). The paraxial ray-based relative location is modified for microseismicity (Grechka et al.2015). For these schemes, picking the corresponding arrival times is an inevitable step and can be determined manually or by using automatic algorithms (Schimmel & Paulssen 1997; Withers et al.1998; Wu et al.2016). However, obtaining it remains tedious and time-consuming. Furthermore, the chosen algorithm may struggle to produce reliable pickings for a low signal-to-noise ratio (SNR) wavefield or complicated phases. In contrast to the above inversion-based algorithms, migration-based location methods make direct use of recorded wavefield information without identifying or measuring phase arrivals. Array seismology methods, often called ‘delay-and-stack’ methods, are widely used in global seismology to avoid picking (Rost & Thomas 2002). The source scanning algorithm (Kao & Shan 2007), multiscale matched-field algorithms in frequency (Corciulo et al.2012), combining network and array waveform coherence method (Sick & Joswig 2017) and the source mechanism corrected semblance of amplitudes (Anikiev et al.2014; Stanek et al.2015) are the practical implementations of this method that searches the maximum sum of the absolute amplitudes of the receivers with trial location and source excitation time, similar to the Kirchhoff reconstruction (Baker 2005). Given a reasonably accurate velocity model, the time-reversed extrapolation or backpropagation method can focus the observed wavefields onto the original source location (Fink 1997; Gajewski & Tessmer 2005). This algorithm is complemented by some physically meaningful imaging conditions, such as the zero lag of autocorrelation and cross-correlation (Artman et al.2010) and its extension into the frequency domain (Witten & Shragge 2015) and interferometric imaging condition (Wang et al.2013). The deconvolution can be extended to improve the imaging of a microseismic event in an elastic medium (Douma & Snieder 2015). The source locations and excitation times can be determined by scanning the entire extrapolated wavefields (general 2-/3-D in space plus 1-D time) for general passive data (Nakata & Beroza 2016). Further, there are some methods proposed to avoid the need of the source excitation time, such as the semblance technique, which finds the point in space that maximizes a semblance measure of the arrival of specific phases (Duncan & Eisner 2010); the semblance-weighted deconvolution image condition between two or more reconstructed source signatures (Haldorsen et al.2013); and the geometric-mean Reverse Time Migration method, which cross-correlates the extrapolated wavefields (Nakata & Beroza 2016). However, the storage and calculation costs of the extrapolated wavefields are tremendous (Nakata & Beroza 2016). Seismic interferometry is a technique that is commonly used to reconstruct the impulse response between receivers without knowing any information about the source. The most typical interferometry in seismology is implemented by utilizing the cross-correlation of any pair of records at the receivers (Campillo & Paul 2003; Schuster et al.2004; Schuster 2009). There are many applications of seismic interferometry in a wide range of contexts. The uses of active interferometry for exploration include the transformation and migration of vertical seismic profiling and inverse vertical seismic profiling (Yu & Schuster 2006), noise attenuation (Halliday et al.2010) and seismic interpolation (Wang et al.2009, 2010). Passive seismic interferometry aims to turn passive seismic measurements into deterministic seismic responses (Wapenaar et al.2010). It shows potential application in hydraulic fracturing (Poliannikov et al.2011). Interferometric cross-correlation migration (CCM) was thus proposed (Schuster et al.2004; Yu & Schuster 2006) to determine the source distribution by using an asymptotic migration scheme along differential traveltime curves, also for tremor location (Li et al.2017). The method can eliminate the excitation source time by generating the virtual gather. However, the method will degrade the spatial resolution due to the additional square term of the source wavelets. To enhance the resolution of the microseismic source locations, we introduce deconvolution migration (DCM) by combining deconvolution interferometry (Vasconcelos & Snieder 2008) with conventional interferometric CCM (Schuster et al.2004). First, we outline the DCM methodology for the localization of hydraulic fracturing induced seismicity, demonstrate its features by using synthetic data and field data from hydraulic fracturing project to find the source location, and test the sensitivity the velocity model errors and data distortions. METHOD Conventionally, frequency-domain passive data u(k, s, ω) at receiver k are represented as the multiplication of the Green's function G(k, s, ω) convolved with a source function W and a term associated with excitation time τs of the source s:   $$u\ \left( {k,s,\omega } \right) = \ W\left( \omega \right)G\left( {k,s,\omega } \right){e^{i\omega {\tau _s}}},$$ (1) where ω represents the circular frequency. The microseismic data, such as from the hydraulic fracturing monitoring, are usually recorded without knowing the source excitation time τs. This makes the microseismic data processing much more difficult than processing the active seismic data. Seismic interferometry is a useful method for obtaining subsurface structures or source distribution (Schuster et al.2004). To determine the source distribution, the CCM (Schuster et al.2004; Schuster 2009) is treated as an asymptotic migration scheme based on an interferometric imaging method. The algorithm can eliminate the unknown excitation source time by generating the virtual gather with cross-correlation operator. Thus, the virtual data $${\tilde{\Phi }_{{\rm{corr}}}}( {k,j} )$$ of any two traces recorded at k (as the master trace) and j can be expressed generally as   $${\tilde{\Phi }_{{\rm{corr}}}}\left( {k,j} \right) = {\rm{\ }}u\left( {k,s,\omega } \right){u^*}\left( {j,s,\omega } \right).$$ (2) In the context of classical scattering theory (Rodberg & Thaler 1967; Vasconcelos & Snieder 2008), the Green's function  G(s, ω) is expressed as the superposition of the direct wave G0 and all of the following scattered waves Gs. To locate the buried source, only the cross-correlation term of the direct wavefield is needed (Schuster et al.2004). Because of its shortest traveltime (first arrival) from the buried source to the receivers, direct wave is unambiguous and avoids multipath wave propagation in time domain, which is decoupled with the later scattered waves from the complex wavefield. Consequently, the direct wave G0 can be extracted with a mask function M:   $$D\left( {s,\omega } \right) = Mu\left( {s,\omega } \right) = W\left( \omega \right){G_0}\left( {s,\omega } \right){e^{i\omega {\tau _s}}}.$$ (3) After inserting eq. (3) into the general cross-correlation eq. (2), it is more convenient to describe the physics of CCM for source location (see  Appendix) and the new virtual gather $${\tilde{\Phi }_{{\rm{corr}}}}$$ as   $${\tilde{\Phi }_{{\rm{corr}}}}\left( {k,j} \right) = {\rm{\ }}D\left( k \right){D^*}\left( j \right) = {\left| {W\left( \omega \right)} \right|^2}{G_0}\left( {k,\omega } \right)G_0^*\left( {j,\omega } \right).$$ (4) Here, the term |W(ω)|2is the square of the source wavelets, which is unknown and results in lower resolution images of the source location, as demonstrated in the following examples. In the interferometry literature, most authors suggest eliminating the term |W(ω)|2 (Vasconcelos & Snieder 2008), but independently estimation of the source function is not a viable option in many cases. To eliminate the square term |W(ω)|2 and enhance the resolution of the source image, we introduce the deconvolution type method to generate the virtual gather. This virtual gather can be expressed as $${\tilde{\Phi }_{{\rm{de}}}}( {k,j} )$$ (see  Appendix):   $${\tilde{\Phi }_{{\rm{de}}}}\left( {k,j} \right) = \frac{{u\left( k \right)}}{{u\left( j \right)}} = \frac{{u\left( k \right){u^*}\left( j \right)}}{{{{\left| {u\left( j \right)} \right|}^2}}} \approx \frac{{{G_0}\left( {k,\omega } \right)G_0^*\left( {j,\omega } \right)}}{{|{G_0}\left( {j,\omega } \right){|^2}}}$$ (5) Here, the deconvolution operator is similar to the cross-correlation operator in eq. (3), except for an eliminated term |W(ω)|2 and an additional term |G0(k, ω)|2. The numerical stability is always compromised by the presence of the denominator, especially because the whole wavefield results in strong oscillatory effects due to the cross terms between the direct and scattered wavefields (Bleistein & Handelsman 1975). Consider using only the direct wavefield with the mask function in eq. (3), the denominator term |G0(k, ω)|2 will contain no cross terms between the direct and the scattered wavefields (see  Appendix), which is relatively stable than the method using the whole wavefield (Bleistein & Handelsman 1975; Vasconcelos & Snieder 2008). Further, a dynamic damping term (Clayton & Wiggins 1976) is added in the denominator to improve the stability:   $${\tilde{\Phi }_{{\rm{de}}}}\left( {k,j} \right) \approx \frac{{{G_0}\left( {k,\omega } \right)G_0^*\left( {j,\omega } \right)}}{{|{G_0}\left( {k,\omega } \right){|^2} + \varepsilon \left\langle {|{G_0}\left( {j,\omega } \right){|^2}} \right\rangle }}.$$ (6) Here, ε is the scaling factor, which is a function of the source wavelet. In this study, it is set to 0.00001. For active seismic data, the standard diffraction-stack migration can find the image points by smearing the event energy along an ellipsoidal zone in the model space, which is determined by the total reflection traveltime from the source to the receiver on the image point (reflector).For passive seismic data, the standard diffraction-stack migration algorithm only uses the one way traveltime (τkx or τxj) for the buried source location point as shown in Fig. 1, but needs the excitation time τs. Further, it can eliminate the unknown excitation time by restructuring the virtual gather, which reassembles the original one-wave path gather (from the source to receiver k or j) into a two-wave path gather (traveltime difference τxj − τxk of the master trace k and trace j). Then, the interferometric migration (CCM or DCM) can image the buried source location by smearing the energy of the virtual gather into the model space along an asymptotic zone, which is determined by the differential traveltime (subtraction by the traveltimes from the source to the two receivers). Thus, the migration kernel can be expressed as Figure 1. View largeDownload slide Illustration of interferometric cross-correlation migration. The left-hand picture denotes the localization algorithm with the standard diffraction-stack migration algorithm, which needs the source excitation time. The right-hand picture shows the deconvolution method, which can eliminate the source excitation time by generating the virtual gather, as shown at the top of the image. The red dashed line is the migration ray path for source localization, where the energy in trace smears over the true source location. Figure 1. View largeDownload slide Illustration of interferometric cross-correlation migration. The left-hand picture denotes the localization algorithm with the standard diffraction-stack migration algorithm, which needs the source excitation time. The right-hand picture shows the deconvolution method, which can eliminate the source excitation time by generating the virtual gather, as shown at the top of the image. The red dashed line is the migration ray path for source localization, where the energy in trace smears over the true source location.  $${e^{ - i\omega \left( {{\tau _{xj}} - {\tau _{xk}}} \right)}}.$$ (7) Here, the term (τxj − τxk) means to remove the traveltime difference of the virtual gather. The idea is of utilizing differential traveltimes on pairwise receivers from common source is analogous to double-difference inversion method (Waldhauser & Ellsworth 2000; Zhang et al.2010). Both two types of methods can reduce the effect of the velocity uncertainties along the very similar ray paths near the source region when the two stations are close to each other (Zhang et al.2010). Then, the image can be obtained by multiplying the virtual data $${\tilde{\Phi }_{{\rm{de}}}}( {k,j} )$$ or $${\tilde{\Phi }_{{\rm{corr}}}}( {k,j} )$$ with the migration kernel. Consequently, the possible locations of microseismic sources are discretized at each grid point of the model. In this way, the final migration equation for estimating the buried source locations can be expressed with a common form:   $${m_{{\rm{dem}}}}\left( x \right){\rm{\ }} = \mathop \sum \limits_{k = 1}^N \mathop \sum \limits_{j = 1}^N \mathop \sum \limits_{\omega = 1}^W {\tilde{\Phi }_{de}}\left( {k,j,\omega } \right){e^{ - i\omega \left( {{\tau _{xj}} - {\tau _{xk}}} \right)}}$$ (8) where N and W denote the number of the gather and frequency band, respectively. The innermost summation over frequency leads to the constructive interference of the migrated single time trace in the virtual gather. The outer summation over receiver j by the virtual gather yields a single image and is further summated over receiver k as a stacking of each single image. The equation is very similar to the standard diffraction-stack migration, except for the input data of the virtual gather and the time-difference-based migration kernel. As a processing algorithm for locating the microseismic sources, first, the virtual gathers are generated by correlating the master trace k with all other traces using the deconvolution operator. Then, the single image mdem(x, k) of the source location is obtained by using the time-difference-based migrating equation. Finally, the final migration image mdem(x) is achieved by stacking the single image results. Synthetic example In the first example, the proposed DCM is tested with a synthetic seismic data set. The synthetic velocity model is extracted from the Marmousi model in Fig. 2(a). The model contains inclined stratified formation and a high-velocity interbedded to simulate the shale reservoir. The model has 401 gridpoints in the horizontal direction and 251 gridpoints in the vertical direction. The grid interval is 10 m in both directions. In total, 100 receivers with a 40-m interval are placed along the surface, and a point source is set at the central location (1500, 2000). According to previous research, the events observed on surface networks are usually dominated by amplitudes with peak frequencies below 50 Hz (Eisner et al.2013). And also based on the data statistics in Fig. 24(b), the peak frequencies are between 30 and 60 Hz. Thus, a Ricker wavelet with a peak frequencyof 50 Hz is considered representative of an actual microseismic source. The recording length of the trace is 0.7 s, with a temporal sampling interval of 1 ms. To simulate the unknown excitation source time in practice, the passive gather in Fig. 2(b) is randomly shifted with excitation time of 50 ms, which only contains the direct wave. Figure 2. View largeDownload slide The synthetic velocity model and its microseismic gather. The model (a) is from part of the Marmousi model; the data are generated by a point source located at the red star (1500, 2000). The source excitation time τs is set to 0.05 s, which is considered an unknown parameter. The model contains inclined stratified formation and some high-velocity interbed to simulate the shale reservoir. Figure 2. View largeDownload slide The synthetic velocity model and its microseismic gather. The model (a) is from part of the Marmousi model; the data are generated by a point source located at the red star (1500, 2000). The source excitation time τs is set to 0.05 s, which is considered an unknown parameter. The model contains inclined stratified formation and some high-velocity interbed to simulate the shale reservoir. The interferometric migration methods utilize differential traveltimes in generated virtual gather, which is positively correlated with the offset from the master trace. Therefore, the traces in vicinity of zero offset (master trace) of the virtual gather have very small differential traveltimes. Due to the zero differential traveltimes of the master trace itself, it will achieve migration results of constant value. And the small differential traveltime of the nearest two traces results serious smearing throughout wider ray path in migration image shown in Figs 3(a) and (b). In other words, those traces have poor sensitivities to the event locations in both horizontal and vertical directions, and also degrade the SNR and resolution of the final stacking image of the source location. Therefore, the traces in vicinity of zero offset should be muted to obtain a higher resolution image. However, the near-offset substantially reduces the effect of the velocity uncertainties along the very similar ray paths near the source region. Having more muted traces means that there is less smear, but weaker source enhancement with stacking. The optimal parameter of range of the vicinity depends on the case. The range of nearest two or three traces is recommended for our algorithm. Besides, the results also show that the proposed DCM migration can weakens the smearing significantly. Figure 3. View largeDownload slide Virtual trace selection in vicinity of zero offset. (a) and (b) are the single trace migration images at the one-time receiver interval offset with the master trace at 1.3 km, using the CCM and DCM methods, respectively. The traces in vicinity of zero offset (master trace) of the virtual gather have serious smearing throughout wider ray path in its migration image. The DCM can weaken the smearing significantly. Figure 3. View largeDownload slide Virtual trace selection in vicinity of zero offset. (a) and (b) are the single trace migration images at the one-time receiver interval offset with the master trace at 1.3 km, using the CCM and DCM methods, respectively. The traces in vicinity of zero offset (master trace) of the virtual gather have serious smearing throughout wider ray path in its migration image. The DCM can weaken the smearing significantly. The single migration image of the virtual cross-correlation and deconvolution gathers with the master trace at 1.3 km are shown in Figs 4(a) and (b), respectively. The single migration images of both methods smear over the true source location. Both methods show obvious migration noise of a certain bandwidth that crosses over the source. Figs 5(a) and (b) show the final result of stacking the single migration images of the virtual gathers using the CCM and DCM, respectively. For magnified views and analysis of migration images around the source location shown in Fig. 6, the stacked image can focus energy on the source accurately with distinguishable amplitude and suppress the migration noise efficiently. Fig. 6(d) indicates the difference of these two migration images in the source region. Obviously, the deconvolution method shows higher resolution imaging of the source location than does CCM, especially near the true source location, which approximatively eliminates the square term of the wavelets. For more detail, Fig. 6(c) shows the vertical and horizontal traces crossing at the maximum amplitudes of each image. The energy is well focused by both methods, where the peaks of the amplitude are both in the correct location. Figure 4. View largeDownload slide Single migration images of the source location. (a) and (b) are the single migration images with the master trace at 1.3 km using the cross-correlation and deconvolution operator, respectively. The images generated with both methods smear over the true source location. The deconvolution migration method shows higher resolution than does cross-correlation migration. Figure 4. View largeDownload slide Single migration images of the source location. (a) and (b) are the single migration images with the master trace at 1.3 km using the cross-correlation and deconvolution operator, respectively. The images generated with both methods smear over the true source location. The deconvolution migration method shows higher resolution than does cross-correlation migration. Figure 5. View largeDownload slide Stacked migration images of the source location with true velocity. (a) and (b) are the stacked images using the cross-correlation migration and deconvolution migration obtained by stacking all 100 single migration images. Compared to the single migration image, the stacked image can focus energy on the source accurately with distinguishable amplitude and suppress the migration noise efficiently. Figure 5. View largeDownload slide Stacked migration images of the source location with true velocity. (a) and (b) are the stacked images using the cross-correlation migration and deconvolution migration obtained by stacking all 100 single migration images. Compared to the single migration image, the stacked image can focus energy on the source accurately with distinguishable amplitude and suppress the migration noise efficiently. Figure 6. View largeDownload slide Magnified views and analysis of migration images around the source location in Figs 5(a) and (b) are correspond to the images in Figs 4(a) and (b), respectively. Its vertical and horizontal traces (d) cross at each maximum amplitudes (denoted with white star) extracted from (a) and (b). (d) is the difference of these two migration images in the source. The DCM result offers a sharper image, with approximately twice the resolution in the depth direction and a relatively higher resolution in the horizontal direction. Figure 6. View largeDownload slide Magnified views and analysis of migration images around the source location in Figs 5(a) and (b) are correspond to the images in Figs 4(a) and (b), respectively. Its vertical and horizontal traces (d) cross at each maximum amplitudes (denoted with white star) extracted from (a) and (b). (d) is the difference of these two migration images in the source. The DCM result offers a sharper image, with approximately twice the resolution in the depth direction and a relatively higher resolution in the horizontal direction. Sensitivity study We address the effects of velocity model errors and data distortions for microseismic source location and improvement in spatial resolution. All of these factors are important for applications that use field data. Velocity model is critical for migration-based location algorithms. However, velocity model error is inevitable due to the limitations in velocity estimation. First, we use a smoothed mean velocity model (smoothed with a 100 × 100m2 2-D rectangle filter) for migration, which is more reasonable than using the correct velocity model. The smoothed velocity along the ray paths is nearly correct. Besides, utilizing differential traveltimes on pairwise receivers from common source can reduce the effect of the velocity uncertainties along the very similar ray paths near the source region. Therefore, the energy is well focused at the correct source location with either methods in Fig. 7. Similarly, the DCM method shows a higher resolution image of the source location than the cross-correlation method. Then, we also add fixed ratio ± 5 and ± 10 per cent velocity errors for the smoothed migration velocities shown in Figs 8–11. The DCM shows a smaller shift in the maximum amplitude of the source localization (error) than does the CCM in the presence of velocity errors, where the surface geometry is more sensitive in the vertical direction than the horizontal direction. As the performance statistics in Fig. 12, the result from DCM will be more reasonable in the presence of velocity errors, which can achieve small sources location errors for more correct localization especially in the vertical direction and narrower error bars of the resolution for sharper peak values at the focused source location. Figure 7. View largeDownload slide Performance statistics in the presence of velocity model error. This example uses a smoothed mean velocity model with a 100 × 100 m2 2-D rectangle filter. Figure 7. View largeDownload slide Performance statistics in the presence of velocity model error. This example uses a smoothed mean velocity model with a 100 × 100 m2 2-D rectangle filter. Figure 8. View largeDownload slide Performance statistics in the presence of velocity model error. This example uses a smoothed mean velocity model with a 100 × 100 m2 2-D rectangle filter, and add fixed ratio + 10 per cent velocity errors. Figure 8. View largeDownload slide Performance statistics in the presence of velocity model error. This example uses a smoothed mean velocity model with a 100 × 100 m2 2-D rectangle filter, and add fixed ratio + 10 per cent velocity errors. Figure 9. View largeDownload slide Performance statistics in the presence of velocity model error. This example uses a smoothed mean velocity model with a 100 × 100 m2 2-D rectangle filter, and add fixed ratio + 5 per cent velocity errors. Figure 9. View largeDownload slide Performance statistics in the presence of velocity model error. This example uses a smoothed mean velocity model with a 100 × 100 m2 2-D rectangle filter, and add fixed ratio + 5 per cent velocity errors. Figure 10. View largeDownload slide Performance statistics in the presence of velocity model error. Performance statistics in the presence of velocity model error. This example uses a smoothed mean velocity model with a 100 × 100 m2 2-D rectangle filter, and add fixed ratio − 5 per cent velocity errors. Figure 10. View largeDownload slide Performance statistics in the presence of velocity model error. Performance statistics in the presence of velocity model error. This example uses a smoothed mean velocity model with a 100 × 100 m2 2-D rectangle filter, and add fixed ratio − 5 per cent velocity errors. Figure 11. View largeDownload slide Performance statistics in the presence of velocity model error. Performance statistics in the presence of velocity model error. This is a smoothed mean velocity model with a 100 × 100m2 2-D rectangle filter, and add fixed ratio − 10 per cent velocity errors. Figure 11. View largeDownload slide Performance statistics in the presence of velocity model error. Performance statistics in the presence of velocity model error. This is a smoothed mean velocity model with a 100 × 100m2 2-D rectangle filter, and add fixed ratio − 10 per cent velocity errors. Figure 12. View largeDownload slide Performance statistics in the presence of velocity model error. (a) and (b) are vertical and horizontal traces cross at the maximum amplitudes in Figs 7–10, respectively. The dashed lines indicate the location errors of the maximum, and the error bars show a common confidence interval with 85 per cent of the peak values as the resolution. Both methods show incorrect locations of the focused energy in the presence of velocity errors. However, DCM tends to be more reasonable and correct, because it tends to be closer to the true source position (dashed line) and achieves higher resolution with smaller error bars by reducing the interference from the square term of the source wavelet. Figure 12. View largeDownload slide Performance statistics in the presence of velocity model error. (a) and (b) are vertical and horizontal traces cross at the maximum amplitudes in Figs 7–10, respectively. The dashed lines indicate the location errors of the maximum, and the error bars show a common confidence interval with 85 per cent of the peak values as the resolution. Both methods show incorrect locations of the focused energy in the presence of velocity errors. However, DCM tends to be more reasonable and correct, because it tends to be closer to the true source position (dashed line) and achieves higher resolution with smaller error bars by reducing the interference from the square term of the source wavelet. Second, we use a smoothed random errors velocity model. We create a spatially structures correlated Gaussian error on velocities, and then smooth it with a 100 × 100m2 2-D rectangle filter. We set two-scale distribution interval of Gaussian error with approximately ± 5 and ± 10 per cent of the original true velocities values; and have multiple tests with random combination of the velocity errors. We show one general possible scenario with ± 5 per cent in Fig. 13 and with ± 10 per cent in Fig. 14. The DCM tends to be more robust in the presence of velocity errors, which can achieve smaller sources location errors and higher resolution image of the source location than the CCM method. Theoretically, the error of the traveltimes is always mitigated by the summation of the random distribution velocity errors along the ray paths. And multiple tests and calculated results have clearly shown that the random Gaussian errors model achieve smaller sources location errors than the corresponding fixed ratio velocity errors model. Figure 13. View largeDownload slide Performance statistics in the presence of velocity model error. This example uses a 5 per cent Gaussian error on velocities with a 100 × 100 m2 2-D rectangle filter smoothing. Figure 13. View largeDownload slide Performance statistics in the presence of velocity model error. This example uses a 5 per cent Gaussian error on velocities with a 100 × 100 m2 2-D rectangle filter smoothing. Figure 14. View largeDownload slide Performance statistics in the presence of velocity model error. This example uses a 10 per cent Gaussian error on velocities with a 100 × 100 m2 2-D rectangle filter smoothing. Figure 14. View largeDownload slide Performance statistics in the presence of velocity model error. This example uses a 10 per cent Gaussian error on velocities with a 100 × 100 m2 2-D rectangle filter smoothing. Essentially, near-surface medium is consists of continuous velocity field and density, which has no apparent layer structures and reflectivity interfaces. And it can be considered as some velocity anomalies in a large-scale background field. So, we introduces a smooth random model example to simulate the near-surface velocity anomalies. The model in Fig. 15(a) is created by adding smooth strong random velocity anomalies and removing layered structures. And we generate the synthetic event gather in Fig. 15(b) with the same modeling parameters and the same geometry of the previous complex model. Obviously, the arrival time of the direct wave has perturbation due to the strong random velocity anomalies along the ray paths. Then, we also use a smoothed mean velocity model (smoothed with a 100 × 100 m2 2-D rectangle filter) as the migration velocity. The results in Fig. 16 show that both methods focused energy at the correct source location, but degraded with more migration noise than the previous migration images in Fig. 7. Similarly, the DCM method shows a higher resolution image of the source location than the cross-correlation method. To some extent, it gives a preliminary evidence that the proposed method is effective for complex random velocity model. Figure 15. View largeDownload slide The smooth random velocity model and its microseismic gather. The model (a) has no layered structures and some strong random velocity anomalies. Its synthetic event gather (b) has the same parameters and the same geometry of the previous complex model. The arrival time of the direct wave has obvious perturbation due to the strong random velocity anomalies. Figure 15. View largeDownload slide The smooth random velocity model and its microseismic gather. The model (a) has no layered structures and some strong random velocity anomalies. Its synthetic event gather (b) has the same parameters and the same geometry of the previous complex model. The arrival time of the direct wave has obvious perturbation due to the strong random velocity anomalies. Figure 16. View largeDownload slide Performance with smooth random model. This example uses a smoothed mean velocity model with a 100 × 100 m2 2-D rectangle filter. (a) and (b) are the stacked images using the cross-correlation migration and deconvolution migration. Figure 16. View largeDownload slide Performance with smooth random model. This example uses a smoothed mean velocity model with a 100 × 100 m2 2-D rectangle filter. (a) and (b) are the stacked images using the cross-correlation migration and deconvolution migration. Based on the results, the existence of a coincidence localization of maximum amplitude by DCM and CCM indicates a well-estimated migration velocity as shown in Fig. 6, whereas the shift of the localization by DCM and CCM indicates the presence of velocity model errors as shown in Figs 7–13. Waveform distortions test Noise and distortions are inevitable features of seismograms. We add two types of noise model to the synthetic records, contains a Gaussian random noise and actual noise extracted from our field data. First, we use bandlimited random Gaussian noise to the synthetic records as background noise. The random Gaussian noise and its frequency contents of each trace are shown in Figs 17(a) and (b), respectively. It has the random attribute in time and frequency domain. We add two-level amplitude of the random noise, approximately 20 and 50 per cent of the maximum amplitude of the noise-free data. The results in Figs 18 and 19 show that both methods still provide an amplitude peak around the true source location. In addition, the DCM shows higher resolution of the focused source. However, the deconvolution operator is more sensitive, as evidenced by a significantly larger artefact surrounding the focused source. Its resolution is negatively correlated with the noise level. Figure 17. View largeDownload slide The random Gaussian noise model. (a) and (b) are the random Gaussian noise in time and frequency domains, respectively. Figure 17. View largeDownload slide The random Gaussian noise model. (a) and (b) are the random Gaussian noise in time and frequency domains, respectively. Figure 18. View largeDownload slide Performance statistics in the presence of data distortions. This example is with bandlimited random Gaussian noise of 20 per cent of the maximum amplitude of noise-free data. Figure 18. View largeDownload slide Performance statistics in the presence of data distortions. This example is with bandlimited random Gaussian noise of 20 per cent of the maximum amplitude of noise-free data. Figure 19. View largeDownload slide Performance statistics in the presence of data distortions. This example is with bandlimited random Gaussian noise of 50 per cent of the maximum amplitude of noise-free data. Both methods provide an energy-focused amplitude peak around the true source location with two-level random Gaussian noise. DCM has a relatively higher resolution than does CCM in the presence of data distortions. Figure 19. View largeDownload slide Performance statistics in the presence of data distortions. This example is with bandlimited random Gaussian noise of 50 per cent of the maximum amplitude of noise-free data. Both methods provide an energy-focused amplitude peak around the true source location with two-level random Gaussian noise. DCM has a relatively higher resolution than does CCM in the presence of data distortions. Second, we use an actual noise model as the coherent noise from producing environments. The actual model is extracted and modified from our field data. The actual noise model and its frequency contents of each trace are shown in Figs 20(a) and (b), respectively. It is characterized by a burst of energy in time domain, which has clustered frequency content not random distribution. This kind of noise may originate from traffic activities or producting environments (Birnie et al. 2016). We use the actual noise model with approximately 20 per cent of the maximum amplitude of the noise-free data. The results in Fig. 21 show that both methods still provide an amplitude peak around the true source location. In addition, the DCM shows higher resolution of the focused source. Figure 20. View largeDownload slide The actual noise model. (a) is the actual noise extracted from our field data and (b) its frequency contents of each trace. Figure 20. View largeDownload slide The actual noise model. (a) is the actual noise extracted from our field data and (b) its frequency contents of each trace. Figure 21. View largeDownload slide Performance statistics in the presence of data distortions. This example is with actual noise. Both methods provide an energy-focused amplitude peak around the true source location with two-level random Gaussian noise. DCM has a relatively higher resolution than does CCM in the presence of data distortions. Figure 21. View largeDownload slide Performance statistics in the presence of data distortions. This example is with actual noise. Both methods provide an energy-focused amplitude peak around the true source location with two-level random Gaussian noise. DCM has a relatively higher resolution than does CCM in the presence of data distortions. To a certain extent, the presence of additional noise (both random and coherent noise) will make the deconvolution operator unstable and yield a degraded performance close to that of the CCM. Therefore, a more accurate operator for the denominator will improve the resolution of the deconvolution, which is meaningful for improving the stability and resolution of the proposed method. Field example The method is also tested using a field data set from a hydraulic fracturing project in North China. The geometry of microseismic monitoring is shown in Fig. 22. Six lines of a star-like surface array of geophones are used to gather passive seismic data with recording the vertical particle velocity. The surface topography of the geophones is varied between 300 and 850 m. The number of geophones per arm varied between 237 and 346, with 1780 geophones in total. The average receiver geophone spacing is approximately 20.2 m (varied between 19 and 22 m). The continuous microseismic data is acquired at 1 kHz. One recording segment of contain strong event is selected with the length of 0.5 s. The vertical production well is located near the centre of the array. According to the field record, the burial depth of shale reservoir is about 4000–4200 m, where the majority of events are expected to be located. The velocity model assumes a horizontal-layer medium and is built based on the experiment sonic log at the injection well, as shown in Fig. 23. Figure 22. View largeDownload slide Map view of the surface survey geometry. The circles denote the geophones within six survey lines, and the red line is used for the following migration using a 2-D approximation. The injection well is in the centre of the star-like geometry. The shale reservoir is at a depth of 4000–4200 m. Figure 22. View largeDownload slide Map view of the surface survey geometry. The circles denote the geophones within six survey lines, and the red line is used for the following migration using a 2-D approximation. The injection well is in the centre of the star-like geometry. The shale reservoir is at a depth of 4000–4200 m. Figure 23. View largeDownload slide The built velocity model. The velocity model assumes a horizontal-layer medium and is built based on the experiment sonic log at the injection well. Figure 23. View largeDownload slide The built velocity model. The velocity model assumes a horizontal-layer medium and is built based on the experiment sonic log at the injection well. Without loss of generality, one line with 294 geophones is selected in this study (the red line shown in Fig. 22), which has an approximate average length of the all arm lines. The selected raw data have low SNR shown in Fig. 24(a), which is contaminated with noise and distortions. However, the direct wave is relative unambiguous, which is decoupled with the later complex scattered wavefield. As discussed in the previous section, the data distortions will make the deconvolution operator unstable and performance with degraded resolution. The pre-processing is the most vital step to find the source location accurately. First, the direct wave is isolated from the wavefield with a mask function as the abovementioned examples. Then, a bandpass filter is applied for mainly removing the noise in the site acquisition. The bandwidth parameters are chosen based on the frequency contents as shown in Fig. 24(b), where the most of the signal is concentrated between 30 and 60 Hz. The final processed record is shown in Fig. 24(c). The SNR is effectively improved. Figure 24. View largeDownload slide Data pre-processing. (a) is the full recorded wavefield and the extracted direct wavefield. (b) is its frequency contents of each trace. (c) is the final processed record.The signal-to-noise ratio (SNR) is effectively improved. Figure 24. View largeDownload slide Data pre-processing. (a) is the full recorded wavefield and the extracted direct wavefield. (b) is its frequency contents of each trace. (c) is the final processed record.The signal-to-noise ratio (SNR) is effectively improved. The single migration images generated using the cross-correlated virtual gather and deconvolved virtual gather with the master trace at 2.2 km are shown in Figs 25(a) and (b), respectively. The single migration image shows a low-resolution image with a strong band-like migration ray path crossing the source. The final results with the two methods are obtained by stacking all 294 migration images of the virtual gathers, shown in Fig. 26. Note that the two final stacked images can enhance the source location with a distinguishable focused amplitude and suppress the migration noise efficiently. As the indication in the vertical and horizontal traces in Fig. 27(c), the deconvolution type method performs a higher resolution image than does the cross-correlation type method, especially in the depth direction. The migration image reveals the location to be at (2820, 4090) by CCM and (2820, 4100) by DCM, both of which are in the shale reservoir of 4000–4200 m. Based on the performance statistics sensitivity study, the existence of shift localization indicates the presence of velocity model errors. Actually, an accurate estimation of the field velocity model from sonic log is very difficult to determine by the current methods, especially for the strong attenuation and anisotropy of the complex near-surface. However, the shift localization is relatively small, which means that using the estimated velocity model is acceptable and reasonable for the migration scheme. Further from the sensitivity study, we know that the localization by DCM at (2820, 4100) should be more correct and closer to the true source location, though the true source location for field data is unknown. Figure 25. View largeDownload slide Single migration image comparison. (a) and (b) are the cross-correlation and deconvolution gathers with a master trace at 2.2 km, respectively. Both images smear over the true source location with relatively strong migration noise. Figure 25. View largeDownload slide Single migration image comparison. (a) and (b) are the cross-correlation and deconvolution gathers with a master trace at 2.2 km, respectively. Both images smear over the true source location with relatively strong migration noise. Figure 26. View largeDownload slide Stacked migration images of the source location. (a) and (b) are the stacked images using cross-correlation and deconvolution migration obtained by stacking all 294 single migration images. Both methods are better focused at the location point than are the single migration images. Figure 26. View largeDownload slide Stacked migration images of the source location. (a) and (b) are the stacked images using cross-correlation and deconvolution migration obtained by stacking all 294 single migration images. Both methods are better focused at the location point than are the single migration images. Figure 27. View largeDownload slide Magnified views of migration images around the source location. (a) and (b) correspond to the images in Figs 18(a) and (b), respectively. (c) shows the vertical and horizontal traces crossing at the maximum amplitudes extracted from (a) and (b). The DCM shows a sharper peak value than the CCM. Figure 27. View largeDownload slide Magnified views of migration images around the source location. (a) and (b) correspond to the images in Figs 18(a) and (b), respectively. (c) shows the vertical and horizontal traces crossing at the maximum amplitudes extracted from (a) and (b). The DCM shows a sharper peak value than the CCM. DISCUSSION The DCM method provides reasonable microseismic source localization results using both synthetic data and a field data from hydraulic fracturing monitoring. The results show more robust and higher spatial resolution images of the source localization than does CCM. For computational efficiency, DCM is very similar to the CCM in the processing algorithm, except for the operator that generates the virtual gather. When generating the virtual gather, the proposed method uses the same method to calculate the cross-correlation term D(k)D*(j) and the additional numerator |D(j)|2 in the deconvolution operator in frequency, which requires the same calculation. In addition, the deconvolution operator requires the vector division of complex numbers by the length of the numerator. Thus, the deconvolution operator costs more than twice as much as the cross-correlation operator overall. For the migration scheme, both methods have the same complexity of calculation, which is based primarily on the number of the traces and the size of the model. In addition, both methods share a traveltime table and pre-processing in the study. Although the proposed method doubles the calculation time of the cross-correlation method, it is still far less computationally intensive than are time-reversed extrapolation-based methods. Because every single migration permits a high degree of parallelism with different grains, some Graphics Processing Unit acceleration and multithreading may be good choices for hardware acceleration. For the performance of the source localization, the star-like geometry in hydraulic fracturing is considerably more irregular and sparse than the dense and regular acquisition geometry in active seismic exploration. Consequently, it results in an uneven distribution of wave paths, where the region near the centre of star-like geometry will provide high fold number and get better image illumination. Analogous to standard diffraction-stack method, the kernel of the DCM needs a sufficient number of folds to suppress the migration noise. Therefore, the 2-D approximation is used for the proposed method in this study. However, it is still meaningful to extend the theory and its performance of the DCM to 3-D with the irregular and sparse star-like geometry. In this study, the current DCM method only deals with direct wave. The following the scattering wavefield will provide more information and illustrations for the hydraulic fracturing of the reservoir (Reshetnikov et al.2015). Furthermore, locating a single event is not equivalent to locating a hydraulic fracturing. Applications to 3-D data sets and extension to the scattering wavefield remain future research topics. CONCLUSIONS We present a high-resolution microseismic sources localization method. The DCM method combines deconvolution interferometry with interferometric CCM. It can eliminate the unknown excitation source time by generating the virtual gather and can enhance the spatial resolution and robustness of the source localization by removing the square term of the source wavelet in CCM. The results of applying our method to synthetic and field cases prove the validity of the proposed method. Compared to the conventional CCM, the DCM can enhance the spatial resolution of source location with a sharper value and improve the robustness with a smaller error of the source localization especially in the presence of velocity errors. The method can provide valuable information for further source mechanism inversion and can be beneficial for global seismology applications. Acknowledgements The research was funded by National Natural Science Foundation of China (grant no. 41230317) and Strategic Priority Research Programme (grant no. XDB10030500). We also acknowledge the constructive suggestions and manuscript improvement of the editor Herve Chauris, reviewer Pierre-Francios Roux and an anonymous reviewer. REFERENCES Anikiev D., Valenta J., Staněk F., Eisner L., 2014. Joint location and source mechanism inversion of microseismic events: benchmarking on seismicity induced by hydraulic fracturing, Geophys. J. Int. , 198( 1), 249– 258. https://doi.org/10.1093/gji/ggu126 Google Scholar CrossRef Search ADS   Artman B., Podladtchikov I., Witten B., 2010. Source location using time-reverse imaging, Geophys. Prospect. , 58( 5), 861– 873. https://doi.org/10.1111/j.1365-2478.2010.00911.x Google Scholar CrossRef Search ADS   Baker T., 2005. Real-time earthquake location using Kirchhoff reconstruction, Bull. seism. Soc. Am. , 95( 2), 699– 707. https://doi.org/10.1785/0120040123 Google Scholar CrossRef Search ADS   Birnie C., Chambers K., Angus D.A., Stork L., 2016. Analysis and models of pre-injection surface seismic array noise recorded at the Aquistore carbon storage site, Geophys. J. Int. , 206( 2), 1246– 1260. https://doi.org/10.1093/gji/ggw203 Google Scholar CrossRef Search ADS   Bleistein N., Handelsman R.A., 1975. Asymptotic Expansions of Integrals , Dover Publishing, Inc. Campillo M., Paul A., 2003. Long-range correlations in the diffuse seismic coda. Science , 299( 5606), 547– 549. https://doi.org/10.1126/science.1078551 Google Scholar CrossRef Search ADS PubMed  Corciulo M., Roux P., Campillo M., Dubucq D., Kuperman W.A., 2012. Multiscale matched-field processing for noise-source localization in exploration geophysics, Geophysics , 77( 5), KS33– KS41. https://doi.org/10.1190/geo2011-0438.1 Google Scholar CrossRef Search ADS   Clayton R. W., Wiggins R.A., 1976, Source shape estimation and deconvolution of teleseismic body waves, Geophys. J. Int. , 47( 1), 151– 177. https://doi.org/10.1111/j.1365-246X.1976.tb01267.x Google Scholar CrossRef Search ADS   Deichmann N., Garcia-Fernandez M., 1992. Rupture geometry from high-precision relative hypocenter locations of microearthquake clusters, Geophys. J. Int. , 110, 501– 517. Google Scholar CrossRef Search ADS   Douma J., Snieder R., 2014. Focusing of elastic waves for microseismic imaging, Geophys. J. Int. , 200( 1), 390– 401. https://doi.org/10.1093/gji/ggu398 Google Scholar CrossRef Search ADS   Duncan P., Eisner L., 2010. Reservoir characterization using surface microseismic monitoring, Geophysics , 75( 5), 139– 146. https://doi.org/10.1190/1.3467760 Google Scholar CrossRef Search ADS   Eisner L., Gei D., Hallo M., Oprs̆al I., Ali M.Y., 2013. The peak frequency of direct waves for microseismic events, Geophysics , 78( 6), 45– 49. https://doi.org/10.1190/geo2013-0197.1 Google Scholar CrossRef Search ADS   Fink M., 1997. Time reversed acoustics, Phys. Today , 50( 3), 34– 40. https://doi.org/10.1063/1.881692 Google Scholar CrossRef Search ADS   Gajewski D., Tessmer E., 2005. Reverse modelling for seismic event characterization, Geophys. J. Int. , 163( 1), 276– 284. https://doi.org/10.1111/j.1365-246X.2005.02732.x Google Scholar CrossRef Search ADS   Geiger L., 1912. Probability method for the determination of earthquake epicentres from the arrival time only (translated from Geiger's 1910 German article), Bull. St. Louis Univ. , 60– 71. Grechka V., Li Z., Howell B., 2015. Relative location of microseismicity, Geophysics , 80( 6), 149– 158. https://doi.org/10.1190/geo2014-0617.1 Haldorsen J.B.U., Brooks N.J., Milenkovic M., 2013. Locating microseismic sources using migration-based deconvolution, Geophysics , 78( 5), KS73– KS84. https://doi.org/10.1190/geo2013-0086.1 Google Scholar CrossRef Search ADS   Halliday D.F., Curtis A., Vermeer P., Strobbia C., Glushchenko A., van Manen D.-J., Robertsson O.A.J., 2010. Interferometric ground-roll removal: attenuation of scattered surface waves in single-sensor data, Geophysics , 75( 2), SA15– SA25. https://doi.org/10.1190/1.3360948 Google Scholar CrossRef Search ADS   Lee W.H.K., Stewart W.S., 1981. Principles and Applications of Microearthquake Networks , Academic Press, Inc. Li K.L., Sgattoni G., Sadeghisorkhani H., Gudmundsson O., 2017. A double-correlation tremor-location method, Geophys. J. Int. , 208( 2), 1231– 1236. https://doi.org/10.1093/gji/ggw453 Google Scholar CrossRef Search ADS   Kao H., Shan S.J., 2007. Rapid identification of earthquake rupture plane using Source-Scanning Algorithm, Geophys. J. Int. , 168( 3), 1011– 1020. https://doi.org/10.1111/j.1365-246X.2006.03271.x Google Scholar CrossRef Search ADS   Maxwell S.C., 2014. Microseismic imaging of hydraulic fracturing: improved engineering of unconventional shale reservoirs, in SEG Distinguished Instructor Short Course , doi:10.1190/1.9781560803164. Nakata N., Beroza G.C., 2016. Reverse time migration for microseismic sources using the geometric mean as an imaging condition, Geophysics , 81( 2), 51– 60. https://doi.org/10.1190/geo2015-0278.1 Google Scholar CrossRef Search ADS   Poliannikov O.V., Malcolm A.E., Djikpesse H., Prange M., 2011. Interferometric hydrofracture microseism localization using neighboring fracture. Geophysics , 76( 6), 27– 36. https://doi.org/10.1190/geo2010-0325.1 Google Scholar CrossRef Search ADS   Reshetnikov A., Kummerow J., Asanuma H., Häring M., Shapiro S.A., 2015. Microseismic reflection imaging and its application to the Basel geothermal reservoir, Geophysics , 80( 6), WC39– WC49. https://doi.org/10.1190/geo2014-0593.1 Google Scholar CrossRef Search ADS   Rodberg L.S., Thaler R.M., 1967. Introduction to the Quantum Theory of Scattering , Academic Press, Inc. Rost S., Thomas C., 2002. Array seismology: methods and applications, Rev. Geophys. , 40( 3), 2-1–2-27. https://doi.org/10.1029/2000RG000100 Google Scholar CrossRef Search ADS   Schuster G.T., Yu J., Sheng J., Rickett J., 2004. Interferometric/daylight seismic imaging, Geophys. J. Int. , 157( 2): 838– 852. https://doi.org/10.1111/j.1365-246X.2004.02251.x Google Scholar CrossRef Search ADS   Schuster G.T., 2009. Seismic Interferometry, Cambridge Univ. Press. Schimmel M., Paulssen H., 1997. Noise reduction and detection of weak, coherent signals through phase-weighted stacks, Geophys. J. Int. , 130( 2), 497– 505. https://doi.org/10.1111/j.1365-246X.1997.tb05664.x Google Scholar CrossRef Search ADS   Sick B, Joswig M., 2017. Combining network and array waveform coherence for automatic location: examples from induced seismicity monitoring. Geophys. J. Int. , 208( 3), 1373– 1388. https://doi.org/10.1093/gji/ggw468 Google Scholar CrossRef Search ADS   Stanek F., Anikiev D., Valenta J., Eisner J., 2015. Semblance for microseismic event detection. Geophys. J. Int. , 201( 3), 1362– 1369. https://doi.org/10.1093/gji/ggv070 Google Scholar CrossRef Search ADS   Wang C., Cheng J., Yin C., Hong L., 2013. Microseismic event location using reverse-time focusing and interferometric techniques for surface and borehole observations[J], Chin. J. Geophys. , 56( 56), 584– 597. Wang Y., Luo Y., Schuster G., 2009. Interferometric interpolation of missing seismic data, Geophysics , 74( 3), SI37– SI45. https://doi.org/10.1190/1.3110072 Google Scholar CrossRef Search ADS   Wang Y., Dong S., Luo Y., 2010. Model-based interferometric interpolation method, Geophysics , 75( 6), WB211– WB217. https://doi.org/10.1190/1.3505816 Google Scholar CrossRef Search ADS   Waldhauser F., Ellsworth W.L., 2000. A double-difference earthquake location algorithm: method and application to the Northern Hayward Fault, California, Bull. seism. Soc. Am. , 90( 6), 1353– 1368. https://doi.org/10.1785/0120000006 Google Scholar CrossRef Search ADS   Wapenaar K., Draganov D., Snieder R., Campman X., Verdel A., 2010. Tutorial on seismic interferometry: Part 1—Basic principles and applications, Geophysics , 75( 5), 75A195– 75A209. Google Scholar CrossRef Search ADS   Withers M., Aster R., Young C., Beiriger J., Harris M., Moore S., Trujillo J., 1998. A comparison of select trigger algorithms for automated global seismic phase and event detection, Bull. seism. Soc. Am. , 88( 1), 95– 106. Witten B., Shragge J., 2015. Extended wave-equation imaging conditions for passive seismic data, Geophysics , 80( 6), WC61– WC72. https://doi.org/10.1190/GEO2015-0046.1 Google Scholar CrossRef Search ADS   Wu S., Wang Y., Chang X., 2016. Automatic microseismic event detection by band-limited phase-only correlation, Phys. Earth planet. Inter. , 261, 3– 16. https://doi.org/10.1016/j.pepi.2016.09.005 Google Scholar CrossRef Search ADS   Vasconcelos I., Snieder R., 2008. Interferometry by deconvolution: Part 2—theory for elastic waves and application to drill-bit seismic imaging, Geophysics , 73( 3), S129– S141. https://doi.org/10.1190/1.2904985 Google Scholar CrossRef Search ADS   Yu J., Schuster G.T., 2006. Crosscorrelogram migration of inverse vertical seismic profile data, Geophysics , 71( 1), S1– S11. https://doi.org/10.1190/1.2159056 Google Scholar CrossRef Search ADS   Zhang H., Nadeau R.M., Toksoz M.N., 2010. Locating nonvolcanic tremors beneath the San Andreas Fault using a station-pair double-difference location method, Geophys. Res. Lett. , 37( 13), doi:10.1029/2010GL043577. https://doi.org/10.1029/2010GL043577 APPENDIX: DECONVOLUTION INTERFEROMETRY For general cross-correlation, we can expand $${\tilde{\Phi }_{{\rm{corr}}}}( {k,j} )$$ into four terms by inserting eq. (1) into (2) with direct wave G0 and scattered waves Gs:   \begin{eqnarray} && {\tilde {\Phi }_{{\rm{corr}}}}\left( {k, j} \right) = {G_0}\left( k \right)G_0^*\left( j \right) + {G_s}\left( k \right)G_0^*\left(j \right)\\ && + \, {G_0}\left( k \right)G_s^*\left( j \right) + {G_s}\left( k \right)G_s^*\left( j \right). \end{eqnarray} (A1) Each of the four terms on the right-hand side has a different physical meaning (Schuster et al.2004; Vasconcelos & Snieder 2008). For imaging purposes, we want to use only the Gs terms, which are the dominant contribution to the causal scattered wavefield that comes from the correlation between the direct wavefield at G0 and the scattered wavefield Gs. For source location imaging, only the first cross-correlation term of the direct wavefield is needed (Schuster et al.2004). Consequently, the latter three terms in eq. (A1) are dropped by using the muting processing in eq. (3). When considering a full waveform u(i), the deconvolution operator of eq. (5) can be expressed as   $${\tilde{\Phi }_{{\rm{de}}}}\left( {k,j} \right) = \frac{{u\left( k \right){u^*}\left( j \right)}}{{{{\left| {u\left( j \right)} \right|}^2}}}\ = \frac{{{{\tilde{\Phi }}_{{\rm{corr}}}}\left( {k,j} \right)}}{{{{\left| {\tilde{u}\left( j \right)} \right|}^2}}}$$ (A2) where the relationship between deconvolution $${\tilde{\Phi }_{{\rm{de}}}}( {k,j} )$$ and cross-correlation $${\tilde{\Phi }_{{\rm{corr}}}}( {k,j} )$$ can be explicitly identified. The denominator can also be expanded into four terms, as eq. (A1):   \begin{eqnarray} {{\left | {\tilde{u}\left( j \right)} \right|^2} }\\ && = \frac{1}{{{G_0}\left( j \right)G_0^*\left( j \right) \!+\! {G_0}\left( j \right)G_s^*\left( j \right) \!+\! {G_s}\left( j \right)G_0^*\left( j \right) \!+\! {G_s}\left( j \right)G_s^*\left( j \right)}}\\ \end{eqnarray} (A3) Assuming that the terms of the direct wavefield G0 make the most prominent contributions and the wavefield perturbations Gs is small (|Gs|2 ≪ |G0|2), the last term in the denominator can be dropped (Vasconcelos & Snieder 2008):   $${\left| {\tilde{u}\left( j \right)} \right|^{ - 2}} \approx \frac{1}{{{{\left| {{G_0}\left( j \right)} \right|}^2}\left[ {1 + \frac{{{G_0}\left( j \right)G_s^*\left( j \right)}}{{{{\left| {{G_0}\left( j \right)} \right|}^2}}} + \frac{{{G_s}\left( j \right)G_0^*\left( j \right)}}{{{{\left| {{G_0}\left( j \right)} \right|}^2}}}} \right]}}$$ (A4) and the latter two terms can be ignored, because of the muting of the scatter wavefield Gs(s) in eq. (2), which also satisfies the condition of |Gs|2 ≪ |G0|2. Combining eqs (A1) and (A4), we can recover the expression of the deconvolution operator in eq. (5). All of the cross terms between direct and scattered wavefields in the denominator and numerator are dropped. © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: Mar 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations