- Email: [email protected]
Microseismic event location using multi-scale time reversed imaging
Journal of Petroleum Science and Engineering 174 (2019) 144–160
Contents lists available at ScienceDirect
Journal of Petroleum Science and Engineering journal homepage: www.elsevier.com/locate/petrol
Microseismic event location using multi-scale time reversed imaging a,∗
a
b
b
Meng Li , Huifeng Li , Guo Tao , Mohammed Ali , Yuhua Guo a b c
T
c
School of Earth Sciences and Engineering, Xi'an Shiyou University, 710065, Xi'an, China Department of Petroleum Geosciences, Khalifa University of Science and Technology, P.O. Box 2533, Abu Dhabi, United Arab Emirates Xi'an Changqing Technology Engineering Co.,Ltd., 710021, Xi'an, China
A R T I C LE I N FO
A B S T R A C T
Portions of this work were presented in “Locating microseismic events with multi-scale time reversed imaging: a synthetic case study,” SEG 85th Annual International Meeting held in New Orleans, USA, Oct 18–23, 2015
Locating passive seismic sources in a reservoir monitoring system has attracted considerable interest recently due to its ability to image the induced fracture geometry that results from fluid injection and fracturing processes. Time reversed imaging (TRI) techniques are well recognized localization tools, which back-propagates the received microseismic recordings to focus at its real source location without picking various arrival phases. However, the time reversed images are often contaminated due to the strong noise, coherent noise and poor spatial constraints especially in surface microseismic measurements, leading to unreliable location estimations. To minimize these interferences, we propose a multi-scale TRI technique using the shift-invariant dual-tree complex wavelet transform (DTCWT) to decompose the original waveforms into multiple time-frequency domains (different levels). The Birge-Massart threshold is applied to the wavelet coefficients at each level to further attenuate the noise amplitude. The TRI are then applied to the reconstructed waveform component at each level. Both the synthetic and field studies show that the multi-scale cross-correlation of TR images at effective levels substantially improved the quality of the final image of subsurface microseismic events with a much sharper and brighter focus. On the other hand, the images at noise levels may reveal the velocity structure, which is helpful in event location. In addition, the noise components recognized by the multi-scale TRI can be applied to estimate the background noise level, which can be further used in microseismic sensitivity and location uncertain analysis.
Keywords: Microseismic Event location Multi-scale Time reversed imaging
1. Introduction Microseismic monitoring has been widely applied in reservoir monitoring, hydraulic fracturing stimulation, mining rock-burst monitoring and CO2 geological sequestration (Rutledge et al., 2004; Maxwell et al., 2010) especially as the booming development of unconventional resources since its first case study in a waste water injection well in Denver, USA (Healy et al., 1968). The spatial and temporal distribution of events provide invaluable knowledge of geomechanical deformation, subsurface stress changes and fluid movement that results from massive fluid injection or depletion within the reservoirs (Miyazawa et al., 2008). This information is of great significance in enhanced oil recovery, estimating of simulated rock volumes, adjusting of stimulation design, and reservoir risk assessment (Li et al., 2017a). The microseismic events are usually located based on first arrivals of direct P and S waves using ray tracing techniques to minimize the difference between observed and theoretical arrival times in an assumed velocity model (Eisner et al., 2010; Wang et al., 2016a). The
∗
locations can further be refined by using relative location methods like double difference to lower the effect of inaccurate velocity model (Waldhauser and Ellsworth, 2000; Li et al., 2013, 2014). Although advanced location methods are commonly used, the strong background noise resulting in low-quality picks sometimes may dramatically lower the reliability and precision of source locations. In addition, arrival picking is considerably time-consuming especially for massive data processing such as distributed acoustic sensor (DAS) measurements (Li et al., 2015; Wilks et al., 2017). With the improvement of computational capacity and DAS measurements, the migration-based event location techniques that employ both of the kinematic and dynamic characteristics of the data are developed in the past few years (Mcmechan, 1982; Hu and Mcmechan, 2010; Fink, 2006). Time reversed imaging (TRI) is one of the migration based location methods. It reversely propagates the emitted miroseismic energy back to the origin without the need to pick first arrivals (Gajewski and Tessmer, 2005). Lu (2008) investigated the potential of time reversed imaging technique to locate earthquakes and image the salt dome flank. Artman
Corresponding author. E-mail address: [email protected] (M. Li).
https://doi.org/10.1016/j.petrol.2018.11.015 Received 16 July 2018; Received in revised form 12 October 2018; Accepted 6 November 2018 Available online 11 November 2018 0920-4105/ © 2018 Elsevier B.V. All rights reserved.
Journal of Petroleum Science and Engineering 174 (2019) 144–160
M. Li et al.
Fig. 1. Schematic diagram of DTCWT decomposition. DTCWT employs two different real DWTs, the first (tree a) gives the real part of the complex coefficients while the second (tree b) gives the imaginary part. s(n) refers to the original signal, g0(n), g1(n) refers to low-pass/high pass filter pair for tree a and h0(n) and h1(n) refers to low-pass/high pass filter pair for tree b, the downward arrow represents the down-sample process.
Fig. 2. Wavelet components of DTCWT (left) and DWT (right) at each level for decomposition of Ricker wavelet with sequential time shift.
and strong coherent noise. Moreover, determining event locations from TR images with strong artifacts and unclear focus can be difficult. Therefore, developing an improved TRI technique that attenuates the artifacts to provide a TR image with a clear and high-resolution is of great significance for continuous data processing. In this paper, we propose a multi-scale TRI method to minimize the interference of strong noise and improve the spatial resolution of the location image. The method employs the shift-invariant dual-tree complex wavelet transform (DTCWT) to decompose the original waveforms into multiple time-frequency domains (different levels) for the purpose of attenuating the strong background noise and enhancing the energy of effective events (Kingsbury, 2001). TRI are then applied to the waveform component at each level and the final location image is given by the multiplication of the images at effective levels. The 3D synthetic and field data examples show that the multi-scale TRI works well in extremely low S/N with clear energy focus and improved spatial solution while conventional TRI fails to locate the events.
et al. (2010) located the microseismic events using auto-correlations and cross-correlations of separated P and S potentials based on Helmholtz decomposition (Aki and Richards, 2002; Li et al., 2017b) in a 2D isotropic media. Witten and Artman (2010) mitigated the imaging artifacts by dividing the signal TRI by noise TRI and employed a confidence threshold to set the minimum displayed image values. Zhu (2014) used a method to compensate the attenuation and dispersion to make TRI effective in attenuating media. Douma and Snieder (2015) improved the spatial and temporal resolutions of TRI using deconvolution. Nakata and Beroza (2016) alleviated the imaging artifacts by crosscorrelating the wavefields extrapolated at each receiver. Li and Van der Baan (2017c) combined the 3C particle velocities and pressure wavefield in TRI to remove the ghost focus on the opposite side of the well for borehole recordings. However, although TRI techniques have been applied in earthquake and downhole data processing to estimate the source locations, its successful application in surface microseismic monitoring is rare due to the extremely low S/N. The high-amplitude ambient noises severely degrade the quality of TR images with strong artifacts, resulting in unreliable location estimations. In addition, the imaging artifacts could arise due to the velocity model uncertainty, sparse recording geometry 145
Journal of Petroleum Science and Engineering 174 (2019) 144–160
M. Li et al.
Fig. 3. Snapshots of time reversed wavefield in various time steps.
data redundancy to analyze the characteristics of non-stationary signals. (Kingsbury, 2001). This technique has been successfully applied in cross-dipole logging data processing to remove the leaky P wave in multiple time-frequency domains (Li et al., 2016) and seismic signal random noise attenuation (Meng et al., 2016; Wang et al., 2016b). This paper employs DTCWT to decompose the noisy microseismic data into multiple time-frequency domains for the purpose of S/N improvement and effective identification of true micro-seismic events without changing the original phases of array microseismic signals. The DTCWT ψc is composed of two separate real wavelets as below:
ψc (t ) = ψr (t ) + jψj (t )
(1)
Here, ψr and ψj are real and imaginary part of DTCWT, respectively. They form a Hilbert transform pair to ensure the linear phase of the wavelet filter. Fig. 1 shows a four-level decomposition process of DTCWT, which comprises two separate trees of real filters. Tree a and b produce the real and imaginary parts of complex coefficients, respectively; s(n) refers to the original signal, g0(n), g1(n) refers to low-pass/ high-pass filter pair for tree a; h0(n) and h1(n) refers to low-pass/highpass filter pair for tree b; the downward arrow represents the downsample process. x0a and x1a are scaling and wavelet coefficients at level 1 for tree a, x0b and x1b are scaling and wavelet coefficients at level 1 for tree b, the rest is expressed in the same manner. The two sets of low filters used to construct the DTCWT satisfy a half-sample delay condition to provide the perfect reconstruction of the original signal. More details about the background knowledge on DWT and DTCWT can be found in Mallat (1989, 1999), Kingsbury (2001) and Selesnick et al. (2005), respectively. In general, the reconstructed signals from scaling (e.g. x0a) and wavelet coefficients (e.g. x1a) at each level are named as scaling and wavelet components, respectively. According to Kingsbury (2001), the original signal can be recovered from either the real or imaginary parts. Thus, the original signal equals to the sum of wavelet components at all levels plus the scaling component at level L (decomposition layer). For simplicity, we name the scaling component at level L as the wavelet component at level L+1. Hence the microseismic signal s(n) can be expressed as:
Fig. 4. Workflow chart for locating microseismic events using multi-scale time reversed image.
2. Method 2.1. Fundamentals of dual tree complex wavelet transform The discrete wavelet transform (DWT) has been widely applied in microseismic S/N improvement and phase identification by decomposing the original signals into multiple time-frequency domains and identify the main features (Yomogida, 1994; Fedi et al., 2004; Capilla, 2006; Mousavi et al., 2016). However, due to the lack of shift invariance of DWT, a small time lapse of the input signal will cause a dramatic change in the wavelet components at different levels (Selesnick et al., 2005). Therefore, the phase of original waveforms cannot be retained. Compared with DWT, Dual-tree complex wavelet transform (DTCWT) is a multi-scale analysis tool with shift invariance and limited
L+1
s (n) =
∑ i=1
(x i ),
(2)
where, xi refers to the wavelet component at level i and L stands for the decomposition layer. To further demonstrate the shift invariance of DTCWT, a set of array waveforms composed of 8 Ricker wavelets with sequential shifts in time domain are decomposed by DTCWT and DWT, respectively. Fig. 2 displays the wavelet components at each level. Here, the wavelet components are the average of reconstructed signals from wavelet coefficients of tree a and tree b. Good shift invariance is shown by 146
Journal of Petroleum Science and Engineering 174 (2019) 144–160
M. Li et al.
Fig. 5. (a) 3D view of velocity model and geometry configuration of receiver arrays used in synthetic study. Red triangular stand for receivers and blue points refer to microseismic sources. (b) Synthetic waveforms (Z-components) of inline receivers. (c) Synthetic waveforms added with random noise following Gaussian distribution (S/N = 0.01 dB). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)
signal at each level.
DTCWT where the shapes of the wavelet components at each level are the same and independent with a phase shift as the signals move forward (the linear moveout of the input signal has been retained), while those components of DWT are distorted severely especially at level 4 and level 5. The very weak amplitude (nearly flat lines) of 2nd, 4th, 6th and 8th traces in level 4 from DWT suggests that the shapes of these waveforms are distorted due to the lack of shift invariance of DWT and the linear moveout is no longer existed. The distortion effect is more obvious in level 5 with more flat lines. This distortion will destroy the inherent phase characteristics of the original microseismic waveforms, leading to inaccurate location results using TRI. Several de-noising techniques like Donoho soft thresholding (Donoho, 1995) can then be applied to the wavelet coefficients to further enhance the S/N. Here, we utilize the adaptive Birge-Massart threshold model (Birgé and Massart, 1997) to suppress the random noise:
c (t ) =
∑ m2 (k ) + 2σ 2t ⎛α + ln k≤t
⎝
n ⎞, t = 1,2, ..., n t⎠
2.2. Microseismic event location using time reversed imaging Locating microseismic events using time reversed imaging technique is based on source-receiver reciprocity and time invariance of the wave equation in non-dissipative media (Larmat et al., 2006), propagating the emitted seismic energy back to its origin without the need of picking arrival time. For multi-component data d at each level, the time reversed wavefield u(x,y,z,t) are used as initial conditions at the receiver positions (Artman et al., 2010):
ρu ¨ − (λ + 2μ) ∇∇⋅u + μ∇ × ∇ × u = d (x, z = 0, T − t )
where u is the backward propagated wavefields obtained by using numerical wave propagators like Fourier method, pseudo-spectral method, finite difference or phase shift method (Gajewski and Tessmer, 2005). Although the information around the source is not perfectly sampled, it is expected that the energy will focus to its real position as the wavefield reversely propagates given an accurate velocity model. The conventional time reversed location image is obtained by the integral of wavefield u along a reasonable time window (Artman et al., 2010) or the snapshot of wavefield u where maximum focusing is observed.
(3)
where, c(t) is the threshold function, m(k) is the wavelet coefficient at level k, σ is the noise variance, α is the penalty factor, n is the total number of wavelet coefficients. The threshold function c(t) reaches its minimum when t = tmin and the final threshold is:
T = |c (tmin )|
(4a)
(4)
Iu =
The wavelet component is then reconstructed from the wavelet coefficient after thresholding to obtain the distribution of the de-noised
∫ |u (v, t )|dt
or Iu = max(|u (v , t )|)
(5)
where v is the space model. The event location (x0, y0, z0) is then 147
Journal of Petroleum Science and Engineering 174 (2019) 144–160
M. Li et al.
Fig. 6. Wavelet components of synthetic waveforms (Z-components) at inline receivers at each level. n
estimated by choosing the grid with the greatest value of Iu. Fig. 3 shows the back-propagating snapshots of wavefield at various time steps in a simple isotropic homogenous velocity model where the source and receivers are located in the center and surface of the model, respectively. As expected, the wave fronts converge as reversely propagating and collapsing at the true source location when time equals to 0 ms, which is the origin time. The focus of the high energy in the location image Iu corresponds to the event location. The radius of the focus implies the location resolution, which depends on the acquisition geometry and bandwidth of the measured microseismic waveforms. The origin time tt is then estimated by scanning the entire time axis to obtain the maximum energy around the imaged focus (Kremers et al., 2011):
⎞ ⎛ tt = argmax ⎜ |u (v, t )|2 d v⎟ ⎠ ⎝W
∫
I˜ =
n
T
∏ Ii = ∏ ∫ |ui |dt, i=1
i=1
i = 1,2…n
0
(7)
where, n refers to the total number of effective images. The effective image is defined by peak signal noise ratio (PSNR), which is a measurement to indicate how easy and confident one can identify and locate the event from the energy focus (Lu, 2008):
⎡ PSNR = 20 log10 ⎢ ⎢ ⎣
⎤ ⎥ ∑x ∈ D n2 (x) d x ⎥ ⎦ max[I (x)]
1 N
(8)
where, D is a sub domain of I(x) which is composed of points whose amplitude is less than a pre-defined threshold. Here the threshold is chosen as the half of the peak amplitude in I(x). n (x) = I (x|x ∈ Dn ) and N is the total number of samples in image domain I(x). When the PSNR at a certain level exceeds a threshold, we consider the image of this level as an effective image. Fig. 4 shows the workflow chart for locating the microseismic events using multi-scale TRI technique. First, the original microseismic recordings are input and decomposed into multiple levels using DTCWT. Second, the wavelet components at all levels are reconstructed after applying the Birge-Massart threshold model to the wavelet coefficients. Third, the TRI is applied to the wavelet component to obtain the location image at each level. Fourth, the effective ratio is calculated to identify the valid images. Fifth, the final location image is obtained by the multiplication of all effective images. Finally, the spatial position and origin time of the event is estimated.
(6)
where, W corresponds to the region around the estimated location (x0, y0, z0). Usually, the length of W is the dominant wavelength.
2.3. Workflow for multi-scale time reversed imaging After the decomposition of original microseismic recordings using DTCWT, TRI is applied to the wavelet component to obtain the location image at each level. The final location image I˜ is formed by the multiplication of several effective images Ii: 148
Journal of Petroleum Science and Engineering 174 (2019) 144–160
M. Li et al.
Fig. 7. Time reversed image (TRI) in East-Depth view of wavelet component at each level. The white stair and circles are velocity model used in location and theoretical source locations, respectively.
The P-wave velocity ranges from 3 to 4 km/s and the VP/VS is a constant of 1.8. One receiver array in the inline direction the other array in the perpendicular direction (red triangles) are deployed at the surface with spacing of 10 m, as shown in Fig. 5a. The source function is Ricker wavelet with the dominant frequency of 5 Hz. The sample rate and total recording time is 0.5 ms and 6s, respectively. The synthetic data are added with random noise following Gaussian distribution and S/N equals to 0.01 dB. Fig. 5b shows the pure synthetic waveform gathered at the inline receivers. For simplicity, only Z components are displayed. Fig. 5c shows the same data but added with strong random noise from which we can hardly identify any effective events. We follow the processing procedure of multi-scale TRI and apply DTCWT and Birge-Massart threshold with the decomposition layer of 5 to reconstruct the wavelet components at various levels, as shown in Fig. 6. To clearly display the difference between wavelet components at each level, only Z-component of the waveforms are displayed. It is obvious that the wavelet components from level 1 to 4 mainly contain random noise while the wavelet components at level 5 and 6 show the effective events. The noises are dispersed and attenuated along the scale axis after multi-
Images constructed by different effective levels provide similar location images with different artifacts. Therefore, the multiplication of all effective levels will highlight the source focus and attenuate the artifacts. The product of effective images can be interpreted as a multidimensional cross-correlation, which has been applied in active seismic processing (Arfken et al., 2013). 3. Synthetic example The proposed approach is first checked with a synthetic data set to illustrate the benefits of multi-scale TRI. Although we use a simple layered velocity model, the validity of the proposed method will not be affected by the complexity of the velocity structure. A 3D staggered grid finite difference code (Wang et al., 2013) is applied to simulate the forward and backward elastic wave propagation (Wang et al., 2015) where the non-splitting perfectly matched layer (Wang and Tao, 2011; Meng and Wang, 2018) is used surrounding the model to eliminate the boundary artificial reflections. A 4-layer velocity model with three point sources (blue dots) of different origin time located at the middle depth of the model is tested. 149
Journal of Petroleum Science and Engineering 174 (2019) 144–160
M. Li et al.
Fig. 8. Time reversed image (TRI) in East-North view of wavelet component at each level.
East-North plane because the velocity in East-North section is homogenous and isotropic. Fig. 9 shows the North - Depth view of time reversed image at each level. Similar to Fig. 7, the images at level 1 to level 4 show the velocity structure and the images at level 5 and level 6 show clear energy focuses. Fig. 10 shows the PSNR values of location image at each level. It is obvious that the PSNR values at level 5 and 6 are much higher than those at level 1 to level 4. Therefore, the location images of level 5 and 6 are multiplied to form the final location imaging result. Fig. 11 show the conventional TRI location result and multi-scale TRI location result in East-North (top), East-Depth (middle) and North-Depth (bottom) view. The image of conventional TRI has no obvious energy focus especially in East-Depth slice while the final image obtained from multiscale TRI results in a much sharper focus. Although the wavelet components of level 1 to level 4 contribute nothing to the final location images, the energy of these useless levels can be used to estimate the background noise level, which are calculated by sum of the root mean square (RMS) of samples in these waveforms:
scale decomposition and de-noising using DTCWT. In addition, the phase lags between the original waveforms are reserved precisely. Then TRI is applied to the wavelet components to generate the location image for each level. Here, a smoothed model of the layered velocity used in forward modelling is employed in multi-scale TRI. Fig. 7 shows the East-Depth view of time reversed image at each level. The TRI from level 1 to level 4 can hardly give an energy focus area while the TRI at level 5 and level 6 shows obvious location ellipse with bright focus. In addition, the time reversed images at level 1 to level 4 seem to show the interface of the velocity model where the upper layer with lighter colors corresponds to a smaller velocity zone, while the lower layer with dark colors corresponds to a high velocity zone. The depth of each interface matches with the velocity model (white stairs in Figs. 7 and 9) in forward modelling very well. This result may not give an accurate velocity model but it is still useful in microseismic event location. Detailed velocity inversion analysis from the noise levels are beyond the scope of this study. Fig. 8 shows the East - North view of time reversed image at each level. The images from level 1 to 4 are smeared along the image domains showing that the waveforms at these levels are mainly composed of random noise. These images do not show any velocity profiles in 150
Journal of Petroleum Science and Engineering 174 (2019) 144–160
M. Li et al.
Fig. 9. Time reversed image (TRI) in North-Depth view of wavelet component at each level. The white stair and circles are velocity model used in location and theoretical source locations, respectively.
Fig. 10. PSNR values of time reversed location images at East-North view, East-Depth view and North-Depth view for all levels.
151
Journal of Petroleum Science and Engineering 174 (2019) 144–160
M. Li et al.
Fig. 11. Comparison of conventional and multi-scale time reversed images in East-North (top), East-Depth (middle) and North-Depth (bottom) views.
the conventional and multi-scale TRI are tested in a field data set, respectively.
Ns
RMS =
∑i = 1 x i2 Ns
(9)
where Ns is the number of samples and xi is the value of sample i and the values of RMS were calculated for each of the three components of the eight receivers. The average RMS of three components is considered as the final RMS for each receiver. The estimated background noise level can be further used in microseismic sensitivity as well as the location uncertain analysis (Halló, 2012). Fig. 12 shows the RMS values for each component per trace of inline receiver array and its average value of three components. The red line is the noise RMS in theory which varies in different trace to give a constant S/N ratio and the blue bars are noise RMS calculated from wavelet components at level 1 and level 2. The noise levels are very similar and slightly lower than those we used in the simulation. This is reasonable because we didn't consider the noise amplitude in level 5 and 6.
4.1. Study area and data acquisition The pilot study was conducted on one of China's giant onshore gas field located in the southwest. The main structure of the field was formed by organic-rich shale within the Longmaxi Formation in Lower Silurian (Tenger et al., 2017). Several stages of hydraulic fracturing were performed in the horizontal well to generate the commercial gas flow. The microseismic data was acquired continuously for 2 days using nine 3-component geophones deployed at the surface surrounding the fracturing well with a sample rate of 4 ms. Fig. 13a shows the geometry configuration of geophone array applied in the field study where the red triangles are the geophones. The layered P-wave velocity model built for location is shown in Fig. 13b which is based on sonic log. The Vp/Vs ratio is a constant of 1.7 due to the absence of S-wave sonic log. Since the near-surface velocity is unknown, the velocity model starts from D235m and we assume the velocities above D235m equal to those in the first layer.
4. Field example To further demonstrate the advantages of the proposed technique, 152
Journal of Petroleum Science and Engineering 174 (2019) 144–160
M. Li et al.
Fig. 12. RMS at each component per trace and its average value of three components.
Fig. 13. (a) 3D view of geometry configuration for geophones (red triangular) in the field study. (b) Horizontal layered velocity used in event location. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)
component at level 3 (Fig. 15) for channel 1. It is evident that the original waveform contains a colored noise with the frequency band of 7–9Hz (highlighted with red ellipse in Fig. 16a). In contrast, the wavelet component at level 3 attenuates this colored noise, demonstrating the proposed method is able to suppress both of the white and colored noise. Then TRI is applied to the wavelet components to generate the location image for each level. Fig. 17 shows the East - North view of time reversed image at each level where a brighter color shows a higher energy. The TRI at level 1 to level 2 are smeared along the image domains while the TRI from level 3 to level 6 shows bright focuses. Fig. 18 shows the East - Depth view of time reversed image at each level. The TRI energy fails to focus at level 1 and level 2 while the images from level 3 to level 6 show a bright eclipse although the background energy
4.2. Comparison of TRI and multi-scale TRI locations Fig. 14 shows an example of microseismic event with strong noise and no obvious onsets. The red, blue and black lines show the X, Y and Z components of the original waveform, respectively. Although some features of the effective event can be identified especially in channel 5 (around 6s) and 27 (around 8s), the main energy of this microseismic event is still hidden in strong background noise. Fig. 15 shows the reconstructed wavelet components at each level after multi-scale decomposition and denoising using DTCWT and Birge-Massart threshold model. It is obvious that the wavelet components at level 1 and level 2 are mainly composed of strong random noise while the wavelet components from level 3 to level 6 show clear arrivals. Fig. 16 compares the frequency spectrum of the original signal (Fig. 14) and wavelet
153
Journal of Petroleum Science and Engineering 174 (2019) 144–160
M. Li et al.
for location with TRI in East – Depth view and find that the velocity interface matches with the velocity model very well. Different with the synthetic example where the velocity structures only displayed in the noise components (wavelet components from level 1 to level 4), the strong velocity interface almost displayed on the wavelet components at all the levels. Fig. 19 shows the North - Depth view of TRI at each level. Similar to Fig. 18, the images at level 1 and level 2 show velocity interface only while the images from level 3 to level 6 show clear energy focuses overlaid with a strong velocity interface. Therefore, the location images from level 3 to 6 are combined to form the final location image. Fig. 20 shows the conventional TRI location result and multi-scale TRI location result in East-North (top), East-Depth (middle) and North-Depth (bottom) view. It is evident that conventional TRI fails to converge the back propagated energy to its real origin while multi-scale TRI is able to give a brighter and sharper focus due to its ability to strengthen the effective signal after multi-scale decomposition and wavelet denoising. The migration-based location methods are valid in Gaussian-noise data and the DTCWT inherits the virtue of wavelet transform that effectively attenuates the white noise, making the proposed method superior to the conventional migration-based methods in strong white noise environment. In addition, since the DTCWT is able to analyze the waveform's characters both in time and frequency domain, analyzing its reversed image in each level (frequency domain) can reveal the frequency band of useful signal, thus, can attenuate the noises with banded spectrum.
Fig. 14. An example of microseismic event with strong noise and no obvious onsets. The red, blue and black lines show the X, Y and Z components of the original waveform, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)
above D360m is strong which may be due to the velocity uncertainty of the shallow formation or the poor constraints of the geophone array in depth direction. The red circles are inferred event locations. It is interesting that a velocity interface around D360m exists in East – Depth view. We overlay the layered velocity model (red curve in Fig. 18) built
Fig. 15. Wavelet components of the original waveforms at each level. The red, blue and black lines show the X, Y and Z components of the original waveform, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) 154
Journal of Petroleum Science and Engineering 174 (2019) 144–160
M. Li et al.
Fig. 16. Frequency spectrum of the original microseismic waveform (a) and wavelet component at level 3 (b) for channel 1.
Fig. 17. Time reversed image (TRI) in East-North view of wavelet component at each level.
5. Discussion
5.1. Revealed velocity structure from TRI
Here, we address the issues of velocity structures revealed in the time reversed images, the improved spatial resolution due to multiscale multiplication and the influence of source mechanism on the imaged focus.
Instead of using a layered velocity model as shown in Fig. 13b, we use a homogenous velocity model in multi-scale TRI. The P-wave velocity is set to be 3411 m/s and the other parameters remained the same. Fig. 21a compares the layered velocity model (blue stairs) and homogeneous velocity of P wave (red straight line). The homogeneous velocity can be viewed as an equivalent medium of the layered
155
Journal of Petroleum Science and Engineering 174 (2019) 144–160
M. Li et al.
Fig. 18. Time reversed image (TRI) in East-Depth view of wavelet component at each level.
Fig. 19. Time reversed image (TRI) in North-Depth view of wavelet component at each level.
156
Journal of Petroleum Science and Engineering 174 (2019) 144–160
M. Li et al.
Fig. 20. Comparison of conventional and multi-scale time reversed images in East-North (top), East-Depth (middle) and North-Depth (bottom) views.
Fig. 21. (a) Comparison of layered (blue) and homogenous (red) P-wave velocity model and (b) East-Depth view of revealed velocity structure of wavelet components at level 1 for homogeneous velocity model. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)
velocity model. The multiple-scattering imaging or noise tomography might be the reason that the multi-scale TRI can generate a velocity interface at the depth where the velocity has huge difference. Although the revealed velocity interface is not adequate for establishing an elaborated velocity model to improve the event locations, it still can be used to calibrate the depth of initial velocity model.
formation. Fig. 21b shows the East-Depth of time reversed image obtained from wavelet component at level 1 overlapped with the layered velocity model. Compared with Fig. 18, the velocity interface is still apparent and matches with the layered formation well although there is small misfit in depth direction. This result implies that the revealed velocity structure of multi-scale TRI is not severely affected by the input 157
Journal of Petroleum Science and Engineering 174 (2019) 144–160
M. Li et al.
Fig. 22. Spatial resolution of conventional TRI and multi-scales TRI in East (a), North (b) and Depth (c) directions.
Fig. 23. Acquisition geometry (a) and imaged focuses for horizontal (b), vertical (c) and 45° point force (d) source in a 2D synthetic case study.
directions. It is obvious that the spatial resolutions of multi-scale TRI are higher than those of conventional TRI and TRI at level 3. The reason is that the strong random noise in the original data is dispersed and the effective signal is strengthened when it is decomposed into multiple time-frequency domains with appropriate noise suppression techniques. In addition, the multiplication of several effective TRI certainly sharpens the focused area leading to a higher resolution than that of any single level. The improved spatial resolution is of great significance for multiple sources location using time reversed techniques since the
5.2. Spatial resolution The spatial resolution of multi-scale TRI is also investigated. The amplitude of the time reversed images at the depth and north coordinate of the inferred source location can be viewed as the spatial resolution in the east direction. The spatial resolutions in depth and north directions are estimated in a similar manner. Fig. 22 compares the spatial resolution of conventional TRI (blue circle), multi-scale TRI (red square) and TRI at level 3 (black diamond) in east, north and depth 158
Journal of Petroleum Science and Engineering 174 (2019) 144–160
M. Li et al.
References
microseismic events are often distributed in a close area within a short time during hydraulic fracturing.
Aki, K., Richards, P., 2002. Quantitative Seismology. University Science Books. Arfken, G.B., Weber, H.J., Harris, F.E., 2013. Mathematical Methods for Physicists, seventh ed. Elsevier. Artman, B., Podladtchikov, I., Witten, B., 2010. Source location using time-reverse imaging. Geophys. Prospect. 58, 861–873. Birgé, L., Massart, P., 1997. From Model Selection to Adaptive Estimation. Springer, New York, pp. 55–87. Capilla, C., 2006. Application of the haar wavelet transform to detect microseismic signal arrivals. J. Appl. Geophys. 59 (1), 36–46. Donoho, D.L., 1995. Denoising by soft-thresholding. IEEE Trans. Inf. Theor. 41 (3), 613–627. Douma, J., Snieder, R., 2015. Focusing of elastic waves for microseismic imaging. Geophys. J. Int. 200 (1), 390–401. Eisner, L., Hulsey, B.J., Duncan, P., Jurick, D., Werner, H., Keller, W., 2010. Comparison of surface and borehole locations of induced seismicity. Geophys. Prospect. 58 (5), 809–820. Fedi, M., Primiceri, R., Quarta, T., Villani, A.V., 2004. Joint application of continuous and discrete wavelet transform on gravity data to identify shallow and deep sources. Geophys. J. Roy. Astron. Soc. 156 (1), 7–21. Fink, M., 2006. Time-reversal acoustics in complex environments. Geophysics 71 (4), S1151–S1164. Gajewski, D., Tessmer, E., 2005. Reverse modelling for seismic even characterization. Geophys. J. Int. 163, 276–284. Halló, M., 2012. Microseismic Surface Monitoring Network Design - Sensitivity and Accuracy: 74th Conference and Exhibition. EAGE Extended Abstracts, pp. P021. Healy, J.H., Rubey, W.W., Griggs, D.T., Raleigh, C.B., 1968. The Denver earthquakes. Science 161, 1301–1310. Hu, L.Z., Mcmechan, G.A., 2010. Elastic finite-difference modelling and imaging for earthquake sources. Geophys. J. Roy. Astron. Soc. 95 (2), 303–313. Kingsbury, N.G., 2001. Complex wavelets for shift invariant analysis and filtering of signals. J. Appl. Comput. Harmon. Anal. 10 (3), 234–253. Kremers, S., Fichtner, A., Brietzke, G.B., Igel, H., Larmat, C., Huang, L., 2011. Exploring the potentials and limitations of the time-reversal imaging of finite seismic sources. Solid Earth Discuss. 2 (1), 95–105. Larmat, C., Montagner, J.P., Fink, M., Capdeville, Y., Tourin, A., Eric, C., 2006. Timereversal imaging of seismic sources and application to the great Sumatra earthquake. Geophys. Res. Lett. 33 (19), 1–4. Li, J., Li, C., Morton, S.A., Dohmen, T., Katahara, K., Toksoz, M.N., 2014. Microseismic joint location and anisotropic velocity inversion for hydraulic fracturing in a tight Bakken reservoir. Geophysics 79 (5), C111–C122. Li, J., Zhang, H., Rodi, W.L., Toksoz, M.N., 2013. Joint microseismic location and anisotropic tomography using differential arrival times and differential backazimuths. Geophys. J. Int. 195 (3), 1917–1931. Li, M., Ali, M.Y., Tao, G., Alakberi, W., Alnuaimi, H., 2017a. Monitoring of induced microseismicity in an onshore oilfield from Abu Dhabi, United Arab Emirates: implications for carbonate reservoir monitoring. J. Petrol. Sci. Eng. 152, 33–48. Li, M., Tao, G., Wang, H., Zhang, K., Vega, S., 2016. An improved multiscale and leaky pwave removal analysis for shear-wave anisotropy inversion with crossed-dipole logs. Petrophysics 57 (3), 270–293. Li, M., Wang, H., Tao, G., 2015. Current and future applications of distributed acoustic sensing as a new reservoir geophysics tool. Open Petrol. Eng. J. 8 (1), 272–281. Li, Y., Wang, H., Fehler, M.C., 2017b. Wavefield characterization of perforation shot signals from a shale gas reservoir. Phys. Earth Planet. In. 267, 31–40. Li, Z., Van der Baan, M., 2017c. Elastic passive source localization using rotational motion. Geophys. J. Int. 211 (2), 1206–1222. Lu, R., 2008. Time Reversed Acoustics and Applications to Earthquake Location and Salt Dome Flank Imaging. Massachusetts Institute of Technology. Mallat, S., 1999. A Wavelet Tour of Signal Processing. Elsevier. Mallat, S., 1989. Theory for multiresolution signal decomposition: the wavelet representation. IEEE Trans. Pattern Anal. Mach. Intell. 11 (7), 674–693. Maxwell, S.C., Rutledge, J., Jones, R., Fehler, M., 2010. Petroleum reservoir characterization using downhole microseismic monitoring. Geophysics 75 (5) 75A129–75A137. Mcmechan, G.A., 1982. Determination of source parameters by wavefield extrapolation. Geophys. J. Int. 71 (3), 613–628. Meng, C., Wang, H., 2018. A finite element and finite difference mixed approach for modeling fault rupture and ground motion. Comput. Geosci. 113, 54–69. Meng, K., Li, J., Li, Y., 2016. Noise suppression in the dual-tree complex wavelet domain for seismic signal. J. Petrol. Explor. Prod. Technol. 7 (2), 353–359. Miyazawa, M., Venkataraman, A., Snieder, R., Payne, A., 2008. Analysis of microearthquake data at Cold Lake and its applications to reservoir monitoring. Geophysics 73 (3), O15–O21. Mousavi, S.M., Langston, C.A., Horton, S.P., 2016. Automatic microseismic denoising and onset detection using the synchrosqueezed continuous wavelet transform. Geophysics 81 (4), V341–V355. Nakata, N., Beroza, G.C., 2016. Reverse time migration for microseismic sources using the geometric mean as an imaging condition. Geophysics 81 (2). Rutledge, J.T., Phillips, W.S., Mayerhofer, M.J., 2004. Faulting induced by forced fluid injection and fluid flow forced by faulting: an interpretation of hydraulic-fracture microseismicity, Carthage Cotton valley gas field. Bull. Seismol. Soc. Am. 94 (5), 1817–1830. Selesnick, I.W., Baraniuk, R.G., Kingsbury, N.C., 2005. The dual-tree complex wavelet transform. IEEE Signal Process. Mag. 22 (6), 123–151.
5.3. Influence of source mechanisms on imaged focuses Based on equations (5) and (7), we use the integration of L2 norm of the reversed wavefield along the time axis to image the focus. Therefore, the reversed imaging based on equations (5) and (7), where the multiplication of imaged focuses at different levels does not change the physical meaning, shows the averaged radiation pattern of the located microseismic source if the wavefield around the event is perfectly sampled and the three-component data is pure without noise. Fig. 23 shows the imaged focuses for horizontal, vertical and 45° point forces in a simple 2D homogenous model with a perfect recording geometry. The configuration of receiver arrays (blue triangles) and source (red circle) is shown in Fig. 23a. It is clear that the averaged radiation patterns for the located event with various source mechanisms are perfectly presented. However, in real cases, sparse recording geometries are frequently deployed due to the limitations of field conditions. As a result, the characteristics of P- or S-wave radiation pattern may not be sampled by the sparse recording geometries. For example, the horizontal nodal plane in the P-wave radiation pattern of a vertical single point force is not sampled by the surface arrays. Therefore, the shape of the imaged focus, using the proposed method in this paper, is a system response of the recording geometry and source mechanism. 6. Conclusions The images of surface micro-seismic measurements may give unreliable location results in extremely strong ambient noise environments. This paper has shown, with a synthetic study of multiple microseismic sources in 3D layered velocity model and a field study during hydraulic fracturing in southwest of China, that conventional time reversed imaging (TRI) techniques fail to locate the events in a situation of extremely low S/N ratio. By contrast, the multi-scale TRI technique proposed in this paper minimizes the interference of noise and identifies effective events in multiple time-frequency domains. Time reversed images of effective components substantially improved the quality of the final image of subsurface microseismic sources with much sharper focuses. On the other hand, the images of noise components may reveal the velocity interface which is helpful in event location. In addition, the noise components recognized by the multi-scale TRI can be applied to estimate the background noise level. Multi-scale TRI effectively utilizes signals scattered in multiple timefrequency domains to construct a source image, therefore, improving the ability to converge the events of weak energy. However, it requires solving the wave equation independently for each level, which can be computationally demanding. Nevertheless, it is naturally suitable for parallel computation to improve the efficiency. A detailed analysis to estimate elaborated velocity model using multi-scale TRI will be an interesting research direction in future. Acknowledgements We would like to thank Xi’an Changqing Technology Engineering Co.,Ltd. for providing the field data. The authors also appreciate the technical support and help from Dr. Quanfeng Wang, Chengdu University of Technology. Appendix A. Supplementary data Supplementary data to this article can be found online at https:// doi.org/10.1016/j.petrol.2018.11.015. 159
Journal of Petroleum Science and Engineering 174 (2019) 144–160
M. Li et al.
Wang, H., Tao, G., 2011. Wavefield simulation and data-acquisition-scheme analysis for LWD acoustic tools in very slow formations. Geophysics 76 (3), 59–68. Wang, H., Li, M., Shang, X.F., 2016a. Current developments on micro-seismic data processing. J. Nat. Gas Sci. Eng. 32, 521–537. Wang, J.J., Li, Y., Liu, W.R., Xiao-Hong, X.U., Mathematics, S.O., 2016b. Dual-tree complex wavelet domain bivariate method for seismic signal random noise attenuation. Chin. J. Geophys. 59 (8), 3046–3055. Wilks, M., Wuestefeld, A., Thomas, P., Kolltveit, E., Oye, V., 2017. Initial Results from the Field on the Use of DAS as a Viable Microseismic Monitoring System of CCS Sites. EAGE Research Workshop. Witten, B., Artman, B., 2010. Signal-to-noise estimates of time-reverse images. Geophysics 76 (2), 2176–2180. Yomogida, K., 1994. Detection of anomalous seismic phases by the wavelet transform. Geophys. J. Roy. Astron. Soc. 116 (1), 119–130.
Tenger, B., Sheng, B., Yu, L., Yang, Y., Zhang, W., Tao, C., Xi, B., Zhang, Q., Bao, F., Qin, J., 2017. Mechanisms of shale gas generation and accumulation in the ordovician Wufeng-Longmaxi Formation, Sichuan basin, SW China. Petrol. Explor. Dev. 44 (1), 69–78. Zhu, Tieyuan, 2014. Time-reverse modelling of acoustic wave propagation in attenuating media. Geophys. J. Int. 197 (1), 483–494. Waldhauser, F., Ellsworth, W.L., 2000. A double-difference earthquake location algorithm: method and application to the Northern Hayward Fault, California. Bull. Seismol. Soc. Am. 90, 1353–1368. Wang, H., Tao, G., Fehler, M.C., 2015. Investigation of the high-frequency wavefield of an off center monopole acoustic Logging-While-Drilling tool. Geophysics 80 (4), D329–D341. Wang, H., Tao, G., Zhang, K., 2013. Wavefield simulation and analysis with the finiteelement method for acoustic logging while drilling in horizontal and deviated wells. Geophysics 78 (6), D525–D543.
160
Contents lists available at ScienceDirect
Journal of Petroleum Science and Engineering journal homepage: www.elsevier.com/locate/petrol
Microseismic event location using multi-scale time reversed imaging a,∗
a
b
b
Meng Li , Huifeng Li , Guo Tao , Mohammed Ali , Yuhua Guo a b c
T
c
School of Earth Sciences and Engineering, Xi'an Shiyou University, 710065, Xi'an, China Department of Petroleum Geosciences, Khalifa University of Science and Technology, P.O. Box 2533, Abu Dhabi, United Arab Emirates Xi'an Changqing Technology Engineering Co.,Ltd., 710021, Xi'an, China
A R T I C LE I N FO
A B S T R A C T
Portions of this work were presented in “Locating microseismic events with multi-scale time reversed imaging: a synthetic case study,” SEG 85th Annual International Meeting held in New Orleans, USA, Oct 18–23, 2015
Locating passive seismic sources in a reservoir monitoring system has attracted considerable interest recently due to its ability to image the induced fracture geometry that results from fluid injection and fracturing processes. Time reversed imaging (TRI) techniques are well recognized localization tools, which back-propagates the received microseismic recordings to focus at its real source location without picking various arrival phases. However, the time reversed images are often contaminated due to the strong noise, coherent noise and poor spatial constraints especially in surface microseismic measurements, leading to unreliable location estimations. To minimize these interferences, we propose a multi-scale TRI technique using the shift-invariant dual-tree complex wavelet transform (DTCWT) to decompose the original waveforms into multiple time-frequency domains (different levels). The Birge-Massart threshold is applied to the wavelet coefficients at each level to further attenuate the noise amplitude. The TRI are then applied to the reconstructed waveform component at each level. Both the synthetic and field studies show that the multi-scale cross-correlation of TR images at effective levels substantially improved the quality of the final image of subsurface microseismic events with a much sharper and brighter focus. On the other hand, the images at noise levels may reveal the velocity structure, which is helpful in event location. In addition, the noise components recognized by the multi-scale TRI can be applied to estimate the background noise level, which can be further used in microseismic sensitivity and location uncertain analysis.
Keywords: Microseismic Event location Multi-scale Time reversed imaging
1. Introduction Microseismic monitoring has been widely applied in reservoir monitoring, hydraulic fracturing stimulation, mining rock-burst monitoring and CO2 geological sequestration (Rutledge et al., 2004; Maxwell et al., 2010) especially as the booming development of unconventional resources since its first case study in a waste water injection well in Denver, USA (Healy et al., 1968). The spatial and temporal distribution of events provide invaluable knowledge of geomechanical deformation, subsurface stress changes and fluid movement that results from massive fluid injection or depletion within the reservoirs (Miyazawa et al., 2008). This information is of great significance in enhanced oil recovery, estimating of simulated rock volumes, adjusting of stimulation design, and reservoir risk assessment (Li et al., 2017a). The microseismic events are usually located based on first arrivals of direct P and S waves using ray tracing techniques to minimize the difference between observed and theoretical arrival times in an assumed velocity model (Eisner et al., 2010; Wang et al., 2016a). The
∗
locations can further be refined by using relative location methods like double difference to lower the effect of inaccurate velocity model (Waldhauser and Ellsworth, 2000; Li et al., 2013, 2014). Although advanced location methods are commonly used, the strong background noise resulting in low-quality picks sometimes may dramatically lower the reliability and precision of source locations. In addition, arrival picking is considerably time-consuming especially for massive data processing such as distributed acoustic sensor (DAS) measurements (Li et al., 2015; Wilks et al., 2017). With the improvement of computational capacity and DAS measurements, the migration-based event location techniques that employ both of the kinematic and dynamic characteristics of the data are developed in the past few years (Mcmechan, 1982; Hu and Mcmechan, 2010; Fink, 2006). Time reversed imaging (TRI) is one of the migration based location methods. It reversely propagates the emitted miroseismic energy back to the origin without the need to pick first arrivals (Gajewski and Tessmer, 2005). Lu (2008) investigated the potential of time reversed imaging technique to locate earthquakes and image the salt dome flank. Artman
Corresponding author. E-mail address: [email protected] (M. Li).
https://doi.org/10.1016/j.petrol.2018.11.015 Received 16 July 2018; Received in revised form 12 October 2018; Accepted 6 November 2018 Available online 11 November 2018 0920-4105/ © 2018 Elsevier B.V. All rights reserved.
Journal of Petroleum Science and Engineering 174 (2019) 144–160
M. Li et al.
Fig. 1. Schematic diagram of DTCWT decomposition. DTCWT employs two different real DWTs, the first (tree a) gives the real part of the complex coefficients while the second (tree b) gives the imaginary part. s(n) refers to the original signal, g0(n), g1(n) refers to low-pass/high pass filter pair for tree a and h0(n) and h1(n) refers to low-pass/high pass filter pair for tree b, the downward arrow represents the down-sample process.
Fig. 2. Wavelet components of DTCWT (left) and DWT (right) at each level for decomposition of Ricker wavelet with sequential time shift.
and strong coherent noise. Moreover, determining event locations from TR images with strong artifacts and unclear focus can be difficult. Therefore, developing an improved TRI technique that attenuates the artifacts to provide a TR image with a clear and high-resolution is of great significance for continuous data processing. In this paper, we propose a multi-scale TRI method to minimize the interference of strong noise and improve the spatial resolution of the location image. The method employs the shift-invariant dual-tree complex wavelet transform (DTCWT) to decompose the original waveforms into multiple time-frequency domains (different levels) for the purpose of attenuating the strong background noise and enhancing the energy of effective events (Kingsbury, 2001). TRI are then applied to the waveform component at each level and the final location image is given by the multiplication of the images at effective levels. The 3D synthetic and field data examples show that the multi-scale TRI works well in extremely low S/N with clear energy focus and improved spatial solution while conventional TRI fails to locate the events.
et al. (2010) located the microseismic events using auto-correlations and cross-correlations of separated P and S potentials based on Helmholtz decomposition (Aki and Richards, 2002; Li et al., 2017b) in a 2D isotropic media. Witten and Artman (2010) mitigated the imaging artifacts by dividing the signal TRI by noise TRI and employed a confidence threshold to set the minimum displayed image values. Zhu (2014) used a method to compensate the attenuation and dispersion to make TRI effective in attenuating media. Douma and Snieder (2015) improved the spatial and temporal resolutions of TRI using deconvolution. Nakata and Beroza (2016) alleviated the imaging artifacts by crosscorrelating the wavefields extrapolated at each receiver. Li and Van der Baan (2017c) combined the 3C particle velocities and pressure wavefield in TRI to remove the ghost focus on the opposite side of the well for borehole recordings. However, although TRI techniques have been applied in earthquake and downhole data processing to estimate the source locations, its successful application in surface microseismic monitoring is rare due to the extremely low S/N. The high-amplitude ambient noises severely degrade the quality of TR images with strong artifacts, resulting in unreliable location estimations. In addition, the imaging artifacts could arise due to the velocity model uncertainty, sparse recording geometry 145
Journal of Petroleum Science and Engineering 174 (2019) 144–160
M. Li et al.
Fig. 3. Snapshots of time reversed wavefield in various time steps.
data redundancy to analyze the characteristics of non-stationary signals. (Kingsbury, 2001). This technique has been successfully applied in cross-dipole logging data processing to remove the leaky P wave in multiple time-frequency domains (Li et al., 2016) and seismic signal random noise attenuation (Meng et al., 2016; Wang et al., 2016b). This paper employs DTCWT to decompose the noisy microseismic data into multiple time-frequency domains for the purpose of S/N improvement and effective identification of true micro-seismic events without changing the original phases of array microseismic signals. The DTCWT ψc is composed of two separate real wavelets as below:
ψc (t ) = ψr (t ) + jψj (t )
(1)
Here, ψr and ψj are real and imaginary part of DTCWT, respectively. They form a Hilbert transform pair to ensure the linear phase of the wavelet filter. Fig. 1 shows a four-level decomposition process of DTCWT, which comprises two separate trees of real filters. Tree a and b produce the real and imaginary parts of complex coefficients, respectively; s(n) refers to the original signal, g0(n), g1(n) refers to low-pass/ high-pass filter pair for tree a; h0(n) and h1(n) refers to low-pass/highpass filter pair for tree b; the downward arrow represents the downsample process. x0a and x1a are scaling and wavelet coefficients at level 1 for tree a, x0b and x1b are scaling and wavelet coefficients at level 1 for tree b, the rest is expressed in the same manner. The two sets of low filters used to construct the DTCWT satisfy a half-sample delay condition to provide the perfect reconstruction of the original signal. More details about the background knowledge on DWT and DTCWT can be found in Mallat (1989, 1999), Kingsbury (2001) and Selesnick et al. (2005), respectively. In general, the reconstructed signals from scaling (e.g. x0a) and wavelet coefficients (e.g. x1a) at each level are named as scaling and wavelet components, respectively. According to Kingsbury (2001), the original signal can be recovered from either the real or imaginary parts. Thus, the original signal equals to the sum of wavelet components at all levels plus the scaling component at level L (decomposition layer). For simplicity, we name the scaling component at level L as the wavelet component at level L+1. Hence the microseismic signal s(n) can be expressed as:
Fig. 4. Workflow chart for locating microseismic events using multi-scale time reversed image.
2. Method 2.1. Fundamentals of dual tree complex wavelet transform The discrete wavelet transform (DWT) has been widely applied in microseismic S/N improvement and phase identification by decomposing the original signals into multiple time-frequency domains and identify the main features (Yomogida, 1994; Fedi et al., 2004; Capilla, 2006; Mousavi et al., 2016). However, due to the lack of shift invariance of DWT, a small time lapse of the input signal will cause a dramatic change in the wavelet components at different levels (Selesnick et al., 2005). Therefore, the phase of original waveforms cannot be retained. Compared with DWT, Dual-tree complex wavelet transform (DTCWT) is a multi-scale analysis tool with shift invariance and limited
L+1
s (n) =
∑ i=1
(x i ),
(2)
where, xi refers to the wavelet component at level i and L stands for the decomposition layer. To further demonstrate the shift invariance of DTCWT, a set of array waveforms composed of 8 Ricker wavelets with sequential shifts in time domain are decomposed by DTCWT and DWT, respectively. Fig. 2 displays the wavelet components at each level. Here, the wavelet components are the average of reconstructed signals from wavelet coefficients of tree a and tree b. Good shift invariance is shown by 146
Journal of Petroleum Science and Engineering 174 (2019) 144–160
M. Li et al.
Fig. 5. (a) 3D view of velocity model and geometry configuration of receiver arrays used in synthetic study. Red triangular stand for receivers and blue points refer to microseismic sources. (b) Synthetic waveforms (Z-components) of inline receivers. (c) Synthetic waveforms added with random noise following Gaussian distribution (S/N = 0.01 dB). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)
signal at each level.
DTCWT where the shapes of the wavelet components at each level are the same and independent with a phase shift as the signals move forward (the linear moveout of the input signal has been retained), while those components of DWT are distorted severely especially at level 4 and level 5. The very weak amplitude (nearly flat lines) of 2nd, 4th, 6th and 8th traces in level 4 from DWT suggests that the shapes of these waveforms are distorted due to the lack of shift invariance of DWT and the linear moveout is no longer existed. The distortion effect is more obvious in level 5 with more flat lines. This distortion will destroy the inherent phase characteristics of the original microseismic waveforms, leading to inaccurate location results using TRI. Several de-noising techniques like Donoho soft thresholding (Donoho, 1995) can then be applied to the wavelet coefficients to further enhance the S/N. Here, we utilize the adaptive Birge-Massart threshold model (Birgé and Massart, 1997) to suppress the random noise:
c (t ) =
∑ m2 (k ) + 2σ 2t ⎛α + ln k≤t
⎝
n ⎞, t = 1,2, ..., n t⎠
2.2. Microseismic event location using time reversed imaging Locating microseismic events using time reversed imaging technique is based on source-receiver reciprocity and time invariance of the wave equation in non-dissipative media (Larmat et al., 2006), propagating the emitted seismic energy back to its origin without the need of picking arrival time. For multi-component data d at each level, the time reversed wavefield u(x,y,z,t) are used as initial conditions at the receiver positions (Artman et al., 2010):
ρu ¨ − (λ + 2μ) ∇∇⋅u + μ∇ × ∇ × u = d (x, z = 0, T − t )
where u is the backward propagated wavefields obtained by using numerical wave propagators like Fourier method, pseudo-spectral method, finite difference or phase shift method (Gajewski and Tessmer, 2005). Although the information around the source is not perfectly sampled, it is expected that the energy will focus to its real position as the wavefield reversely propagates given an accurate velocity model. The conventional time reversed location image is obtained by the integral of wavefield u along a reasonable time window (Artman et al., 2010) or the snapshot of wavefield u where maximum focusing is observed.
(3)
where, c(t) is the threshold function, m(k) is the wavelet coefficient at level k, σ is the noise variance, α is the penalty factor, n is the total number of wavelet coefficients. The threshold function c(t) reaches its minimum when t = tmin and the final threshold is:
T = |c (tmin )|
(4a)
(4)
Iu =
The wavelet component is then reconstructed from the wavelet coefficient after thresholding to obtain the distribution of the de-noised
∫ |u (v, t )|dt
or Iu = max(|u (v , t )|)
(5)
where v is the space model. The event location (x0, y0, z0) is then 147
Journal of Petroleum Science and Engineering 174 (2019) 144–160
M. Li et al.
Fig. 6. Wavelet components of synthetic waveforms (Z-components) at inline receivers at each level. n
estimated by choosing the grid with the greatest value of Iu. Fig. 3 shows the back-propagating snapshots of wavefield at various time steps in a simple isotropic homogenous velocity model where the source and receivers are located in the center and surface of the model, respectively. As expected, the wave fronts converge as reversely propagating and collapsing at the true source location when time equals to 0 ms, which is the origin time. The focus of the high energy in the location image Iu corresponds to the event location. The radius of the focus implies the location resolution, which depends on the acquisition geometry and bandwidth of the measured microseismic waveforms. The origin time tt is then estimated by scanning the entire time axis to obtain the maximum energy around the imaged focus (Kremers et al., 2011):
⎞ ⎛ tt = argmax ⎜ |u (v, t )|2 d v⎟ ⎠ ⎝W
∫
I˜ =
n
T
∏ Ii = ∏ ∫ |ui |dt, i=1
i=1
i = 1,2…n
0
(7)
where, n refers to the total number of effective images. The effective image is defined by peak signal noise ratio (PSNR), which is a measurement to indicate how easy and confident one can identify and locate the event from the energy focus (Lu, 2008):
⎡ PSNR = 20 log10 ⎢ ⎢ ⎣
⎤ ⎥ ∑x ∈ D n2 (x) d x ⎥ ⎦ max[I (x)]
1 N
(8)
where, D is a sub domain of I(x) which is composed of points whose amplitude is less than a pre-defined threshold. Here the threshold is chosen as the half of the peak amplitude in I(x). n (x) = I (x|x ∈ Dn ) and N is the total number of samples in image domain I(x). When the PSNR at a certain level exceeds a threshold, we consider the image of this level as an effective image. Fig. 4 shows the workflow chart for locating the microseismic events using multi-scale TRI technique. First, the original microseismic recordings are input and decomposed into multiple levels using DTCWT. Second, the wavelet components at all levels are reconstructed after applying the Birge-Massart threshold model to the wavelet coefficients. Third, the TRI is applied to the wavelet component to obtain the location image at each level. Fourth, the effective ratio is calculated to identify the valid images. Fifth, the final location image is obtained by the multiplication of all effective images. Finally, the spatial position and origin time of the event is estimated.
(6)
where, W corresponds to the region around the estimated location (x0, y0, z0). Usually, the length of W is the dominant wavelength.
2.3. Workflow for multi-scale time reversed imaging After the decomposition of original microseismic recordings using DTCWT, TRI is applied to the wavelet component to obtain the location image at each level. The final location image I˜ is formed by the multiplication of several effective images Ii: 148
Journal of Petroleum Science and Engineering 174 (2019) 144–160
M. Li et al.
Fig. 7. Time reversed image (TRI) in East-Depth view of wavelet component at each level. The white stair and circles are velocity model used in location and theoretical source locations, respectively.
The P-wave velocity ranges from 3 to 4 km/s and the VP/VS is a constant of 1.8. One receiver array in the inline direction the other array in the perpendicular direction (red triangles) are deployed at the surface with spacing of 10 m, as shown in Fig. 5a. The source function is Ricker wavelet with the dominant frequency of 5 Hz. The sample rate and total recording time is 0.5 ms and 6s, respectively. The synthetic data are added with random noise following Gaussian distribution and S/N equals to 0.01 dB. Fig. 5b shows the pure synthetic waveform gathered at the inline receivers. For simplicity, only Z components are displayed. Fig. 5c shows the same data but added with strong random noise from which we can hardly identify any effective events. We follow the processing procedure of multi-scale TRI and apply DTCWT and Birge-Massart threshold with the decomposition layer of 5 to reconstruct the wavelet components at various levels, as shown in Fig. 6. To clearly display the difference between wavelet components at each level, only Z-component of the waveforms are displayed. It is obvious that the wavelet components from level 1 to 4 mainly contain random noise while the wavelet components at level 5 and 6 show the effective events. The noises are dispersed and attenuated along the scale axis after multi-
Images constructed by different effective levels provide similar location images with different artifacts. Therefore, the multiplication of all effective levels will highlight the source focus and attenuate the artifacts. The product of effective images can be interpreted as a multidimensional cross-correlation, which has been applied in active seismic processing (Arfken et al., 2013). 3. Synthetic example The proposed approach is first checked with a synthetic data set to illustrate the benefits of multi-scale TRI. Although we use a simple layered velocity model, the validity of the proposed method will not be affected by the complexity of the velocity structure. A 3D staggered grid finite difference code (Wang et al., 2013) is applied to simulate the forward and backward elastic wave propagation (Wang et al., 2015) where the non-splitting perfectly matched layer (Wang and Tao, 2011; Meng and Wang, 2018) is used surrounding the model to eliminate the boundary artificial reflections. A 4-layer velocity model with three point sources (blue dots) of different origin time located at the middle depth of the model is tested. 149
Journal of Petroleum Science and Engineering 174 (2019) 144–160
M. Li et al.
Fig. 8. Time reversed image (TRI) in East-North view of wavelet component at each level.
East-North plane because the velocity in East-North section is homogenous and isotropic. Fig. 9 shows the North - Depth view of time reversed image at each level. Similar to Fig. 7, the images at level 1 to level 4 show the velocity structure and the images at level 5 and level 6 show clear energy focuses. Fig. 10 shows the PSNR values of location image at each level. It is obvious that the PSNR values at level 5 and 6 are much higher than those at level 1 to level 4. Therefore, the location images of level 5 and 6 are multiplied to form the final location imaging result. Fig. 11 show the conventional TRI location result and multi-scale TRI location result in East-North (top), East-Depth (middle) and North-Depth (bottom) view. The image of conventional TRI has no obvious energy focus especially in East-Depth slice while the final image obtained from multiscale TRI results in a much sharper focus. Although the wavelet components of level 1 to level 4 contribute nothing to the final location images, the energy of these useless levels can be used to estimate the background noise level, which are calculated by sum of the root mean square (RMS) of samples in these waveforms:
scale decomposition and de-noising using DTCWT. In addition, the phase lags between the original waveforms are reserved precisely. Then TRI is applied to the wavelet components to generate the location image for each level. Here, a smoothed model of the layered velocity used in forward modelling is employed in multi-scale TRI. Fig. 7 shows the East-Depth view of time reversed image at each level. The TRI from level 1 to level 4 can hardly give an energy focus area while the TRI at level 5 and level 6 shows obvious location ellipse with bright focus. In addition, the time reversed images at level 1 to level 4 seem to show the interface of the velocity model where the upper layer with lighter colors corresponds to a smaller velocity zone, while the lower layer with dark colors corresponds to a high velocity zone. The depth of each interface matches with the velocity model (white stairs in Figs. 7 and 9) in forward modelling very well. This result may not give an accurate velocity model but it is still useful in microseismic event location. Detailed velocity inversion analysis from the noise levels are beyond the scope of this study. Fig. 8 shows the East - North view of time reversed image at each level. The images from level 1 to 4 are smeared along the image domains showing that the waveforms at these levels are mainly composed of random noise. These images do not show any velocity profiles in 150
Journal of Petroleum Science and Engineering 174 (2019) 144–160
M. Li et al.
Fig. 9. Time reversed image (TRI) in North-Depth view of wavelet component at each level. The white stair and circles are velocity model used in location and theoretical source locations, respectively.
Fig. 10. PSNR values of time reversed location images at East-North view, East-Depth view and North-Depth view for all levels.
151
Journal of Petroleum Science and Engineering 174 (2019) 144–160
M. Li et al.
Fig. 11. Comparison of conventional and multi-scale time reversed images in East-North (top), East-Depth (middle) and North-Depth (bottom) views.
the conventional and multi-scale TRI are tested in a field data set, respectively.
Ns
RMS =
∑i = 1 x i2 Ns
(9)
where Ns is the number of samples and xi is the value of sample i and the values of RMS were calculated for each of the three components of the eight receivers. The average RMS of three components is considered as the final RMS for each receiver. The estimated background noise level can be further used in microseismic sensitivity as well as the location uncertain analysis (Halló, 2012). Fig. 12 shows the RMS values for each component per trace of inline receiver array and its average value of three components. The red line is the noise RMS in theory which varies in different trace to give a constant S/N ratio and the blue bars are noise RMS calculated from wavelet components at level 1 and level 2. The noise levels are very similar and slightly lower than those we used in the simulation. This is reasonable because we didn't consider the noise amplitude in level 5 and 6.
4.1. Study area and data acquisition The pilot study was conducted on one of China's giant onshore gas field located in the southwest. The main structure of the field was formed by organic-rich shale within the Longmaxi Formation in Lower Silurian (Tenger et al., 2017). Several stages of hydraulic fracturing were performed in the horizontal well to generate the commercial gas flow. The microseismic data was acquired continuously for 2 days using nine 3-component geophones deployed at the surface surrounding the fracturing well with a sample rate of 4 ms. Fig. 13a shows the geometry configuration of geophone array applied in the field study where the red triangles are the geophones. The layered P-wave velocity model built for location is shown in Fig. 13b which is based on sonic log. The Vp/Vs ratio is a constant of 1.7 due to the absence of S-wave sonic log. Since the near-surface velocity is unknown, the velocity model starts from D235m and we assume the velocities above D235m equal to those in the first layer.
4. Field example To further demonstrate the advantages of the proposed technique, 152
Journal of Petroleum Science and Engineering 174 (2019) 144–160
M. Li et al.
Fig. 12. RMS at each component per trace and its average value of three components.
Fig. 13. (a) 3D view of geometry configuration for geophones (red triangular) in the field study. (b) Horizontal layered velocity used in event location. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)
component at level 3 (Fig. 15) for channel 1. It is evident that the original waveform contains a colored noise with the frequency band of 7–9Hz (highlighted with red ellipse in Fig. 16a). In contrast, the wavelet component at level 3 attenuates this colored noise, demonstrating the proposed method is able to suppress both of the white and colored noise. Then TRI is applied to the wavelet components to generate the location image for each level. Fig. 17 shows the East - North view of time reversed image at each level where a brighter color shows a higher energy. The TRI at level 1 to level 2 are smeared along the image domains while the TRI from level 3 to level 6 shows bright focuses. Fig. 18 shows the East - Depth view of time reversed image at each level. The TRI energy fails to focus at level 1 and level 2 while the images from level 3 to level 6 show a bright eclipse although the background energy
4.2. Comparison of TRI and multi-scale TRI locations Fig. 14 shows an example of microseismic event with strong noise and no obvious onsets. The red, blue and black lines show the X, Y and Z components of the original waveform, respectively. Although some features of the effective event can be identified especially in channel 5 (around 6s) and 27 (around 8s), the main energy of this microseismic event is still hidden in strong background noise. Fig. 15 shows the reconstructed wavelet components at each level after multi-scale decomposition and denoising using DTCWT and Birge-Massart threshold model. It is obvious that the wavelet components at level 1 and level 2 are mainly composed of strong random noise while the wavelet components from level 3 to level 6 show clear arrivals. Fig. 16 compares the frequency spectrum of the original signal (Fig. 14) and wavelet
153
Journal of Petroleum Science and Engineering 174 (2019) 144–160
M. Li et al.
for location with TRI in East – Depth view and find that the velocity interface matches with the velocity model very well. Different with the synthetic example where the velocity structures only displayed in the noise components (wavelet components from level 1 to level 4), the strong velocity interface almost displayed on the wavelet components at all the levels. Fig. 19 shows the North - Depth view of TRI at each level. Similar to Fig. 18, the images at level 1 and level 2 show velocity interface only while the images from level 3 to level 6 show clear energy focuses overlaid with a strong velocity interface. Therefore, the location images from level 3 to 6 are combined to form the final location image. Fig. 20 shows the conventional TRI location result and multi-scale TRI location result in East-North (top), East-Depth (middle) and North-Depth (bottom) view. It is evident that conventional TRI fails to converge the back propagated energy to its real origin while multi-scale TRI is able to give a brighter and sharper focus due to its ability to strengthen the effective signal after multi-scale decomposition and wavelet denoising. The migration-based location methods are valid in Gaussian-noise data and the DTCWT inherits the virtue of wavelet transform that effectively attenuates the white noise, making the proposed method superior to the conventional migration-based methods in strong white noise environment. In addition, since the DTCWT is able to analyze the waveform's characters both in time and frequency domain, analyzing its reversed image in each level (frequency domain) can reveal the frequency band of useful signal, thus, can attenuate the noises with banded spectrum.
Fig. 14. An example of microseismic event with strong noise and no obvious onsets. The red, blue and black lines show the X, Y and Z components of the original waveform, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)
above D360m is strong which may be due to the velocity uncertainty of the shallow formation or the poor constraints of the geophone array in depth direction. The red circles are inferred event locations. It is interesting that a velocity interface around D360m exists in East – Depth view. We overlay the layered velocity model (red curve in Fig. 18) built
Fig. 15. Wavelet components of the original waveforms at each level. The red, blue and black lines show the X, Y and Z components of the original waveform, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) 154
Journal of Petroleum Science and Engineering 174 (2019) 144–160
M. Li et al.
Fig. 16. Frequency spectrum of the original microseismic waveform (a) and wavelet component at level 3 (b) for channel 1.
Fig. 17. Time reversed image (TRI) in East-North view of wavelet component at each level.
5. Discussion
5.1. Revealed velocity structure from TRI
Here, we address the issues of velocity structures revealed in the time reversed images, the improved spatial resolution due to multiscale multiplication and the influence of source mechanism on the imaged focus.
Instead of using a layered velocity model as shown in Fig. 13b, we use a homogenous velocity model in multi-scale TRI. The P-wave velocity is set to be 3411 m/s and the other parameters remained the same. Fig. 21a compares the layered velocity model (blue stairs) and homogeneous velocity of P wave (red straight line). The homogeneous velocity can be viewed as an equivalent medium of the layered
155
Journal of Petroleum Science and Engineering 174 (2019) 144–160
M. Li et al.
Fig. 18. Time reversed image (TRI) in East-Depth view of wavelet component at each level.
Fig. 19. Time reversed image (TRI) in North-Depth view of wavelet component at each level.
156
Journal of Petroleum Science and Engineering 174 (2019) 144–160
M. Li et al.
Fig. 20. Comparison of conventional and multi-scale time reversed images in East-North (top), East-Depth (middle) and North-Depth (bottom) views.
Fig. 21. (a) Comparison of layered (blue) and homogenous (red) P-wave velocity model and (b) East-Depth view of revealed velocity structure of wavelet components at level 1 for homogeneous velocity model. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)
velocity model. The multiple-scattering imaging or noise tomography might be the reason that the multi-scale TRI can generate a velocity interface at the depth where the velocity has huge difference. Although the revealed velocity interface is not adequate for establishing an elaborated velocity model to improve the event locations, it still can be used to calibrate the depth of initial velocity model.
formation. Fig. 21b shows the East-Depth of time reversed image obtained from wavelet component at level 1 overlapped with the layered velocity model. Compared with Fig. 18, the velocity interface is still apparent and matches with the layered formation well although there is small misfit in depth direction. This result implies that the revealed velocity structure of multi-scale TRI is not severely affected by the input 157
Journal of Petroleum Science and Engineering 174 (2019) 144–160
M. Li et al.
Fig. 22. Spatial resolution of conventional TRI and multi-scales TRI in East (a), North (b) and Depth (c) directions.
Fig. 23. Acquisition geometry (a) and imaged focuses for horizontal (b), vertical (c) and 45° point force (d) source in a 2D synthetic case study.
directions. It is obvious that the spatial resolutions of multi-scale TRI are higher than those of conventional TRI and TRI at level 3. The reason is that the strong random noise in the original data is dispersed and the effective signal is strengthened when it is decomposed into multiple time-frequency domains with appropriate noise suppression techniques. In addition, the multiplication of several effective TRI certainly sharpens the focused area leading to a higher resolution than that of any single level. The improved spatial resolution is of great significance for multiple sources location using time reversed techniques since the
5.2. Spatial resolution The spatial resolution of multi-scale TRI is also investigated. The amplitude of the time reversed images at the depth and north coordinate of the inferred source location can be viewed as the spatial resolution in the east direction. The spatial resolutions in depth and north directions are estimated in a similar manner. Fig. 22 compares the spatial resolution of conventional TRI (blue circle), multi-scale TRI (red square) and TRI at level 3 (black diamond) in east, north and depth 158
Journal of Petroleum Science and Engineering 174 (2019) 144–160
M. Li et al.
References
microseismic events are often distributed in a close area within a short time during hydraulic fracturing.
Aki, K., Richards, P., 2002. Quantitative Seismology. University Science Books. Arfken, G.B., Weber, H.J., Harris, F.E., 2013. Mathematical Methods for Physicists, seventh ed. Elsevier. Artman, B., Podladtchikov, I., Witten, B., 2010. Source location using time-reverse imaging. Geophys. Prospect. 58, 861–873. Birgé, L., Massart, P., 1997. From Model Selection to Adaptive Estimation. Springer, New York, pp. 55–87. Capilla, C., 2006. Application of the haar wavelet transform to detect microseismic signal arrivals. J. Appl. Geophys. 59 (1), 36–46. Donoho, D.L., 1995. Denoising by soft-thresholding. IEEE Trans. Inf. Theor. 41 (3), 613–627. Douma, J., Snieder, R., 2015. Focusing of elastic waves for microseismic imaging. Geophys. J. Int. 200 (1), 390–401. Eisner, L., Hulsey, B.J., Duncan, P., Jurick, D., Werner, H., Keller, W., 2010. Comparison of surface and borehole locations of induced seismicity. Geophys. Prospect. 58 (5), 809–820. Fedi, M., Primiceri, R., Quarta, T., Villani, A.V., 2004. Joint application of continuous and discrete wavelet transform on gravity data to identify shallow and deep sources. Geophys. J. Roy. Astron. Soc. 156 (1), 7–21. Fink, M., 2006. Time-reversal acoustics in complex environments. Geophysics 71 (4), S1151–S1164. Gajewski, D., Tessmer, E., 2005. Reverse modelling for seismic even characterization. Geophys. J. Int. 163, 276–284. Halló, M., 2012. Microseismic Surface Monitoring Network Design - Sensitivity and Accuracy: 74th Conference and Exhibition. EAGE Extended Abstracts, pp. P021. Healy, J.H., Rubey, W.W., Griggs, D.T., Raleigh, C.B., 1968. The Denver earthquakes. Science 161, 1301–1310. Hu, L.Z., Mcmechan, G.A., 2010. Elastic finite-difference modelling and imaging for earthquake sources. Geophys. J. Roy. Astron. Soc. 95 (2), 303–313. Kingsbury, N.G., 2001. Complex wavelets for shift invariant analysis and filtering of signals. J. Appl. Comput. Harmon. Anal. 10 (3), 234–253. Kremers, S., Fichtner, A., Brietzke, G.B., Igel, H., Larmat, C., Huang, L., 2011. Exploring the potentials and limitations of the time-reversal imaging of finite seismic sources. Solid Earth Discuss. 2 (1), 95–105. Larmat, C., Montagner, J.P., Fink, M., Capdeville, Y., Tourin, A., Eric, C., 2006. Timereversal imaging of seismic sources and application to the great Sumatra earthquake. Geophys. Res. Lett. 33 (19), 1–4. Li, J., Li, C., Morton, S.A., Dohmen, T., Katahara, K., Toksoz, M.N., 2014. Microseismic joint location and anisotropic velocity inversion for hydraulic fracturing in a tight Bakken reservoir. Geophysics 79 (5), C111–C122. Li, J., Zhang, H., Rodi, W.L., Toksoz, M.N., 2013. Joint microseismic location and anisotropic tomography using differential arrival times and differential backazimuths. Geophys. J. Int. 195 (3), 1917–1931. Li, M., Ali, M.Y., Tao, G., Alakberi, W., Alnuaimi, H., 2017a. Monitoring of induced microseismicity in an onshore oilfield from Abu Dhabi, United Arab Emirates: implications for carbonate reservoir monitoring. J. Petrol. Sci. Eng. 152, 33–48. Li, M., Tao, G., Wang, H., Zhang, K., Vega, S., 2016. An improved multiscale and leaky pwave removal analysis for shear-wave anisotropy inversion with crossed-dipole logs. Petrophysics 57 (3), 270–293. Li, M., Wang, H., Tao, G., 2015. Current and future applications of distributed acoustic sensing as a new reservoir geophysics tool. Open Petrol. Eng. J. 8 (1), 272–281. Li, Y., Wang, H., Fehler, M.C., 2017b. Wavefield characterization of perforation shot signals from a shale gas reservoir. Phys. Earth Planet. In. 267, 31–40. Li, Z., Van der Baan, M., 2017c. Elastic passive source localization using rotational motion. Geophys. J. Int. 211 (2), 1206–1222. Lu, R., 2008. Time Reversed Acoustics and Applications to Earthquake Location and Salt Dome Flank Imaging. Massachusetts Institute of Technology. Mallat, S., 1999. A Wavelet Tour of Signal Processing. Elsevier. Mallat, S., 1989. Theory for multiresolution signal decomposition: the wavelet representation. IEEE Trans. Pattern Anal. Mach. Intell. 11 (7), 674–693. Maxwell, S.C., Rutledge, J., Jones, R., Fehler, M., 2010. Petroleum reservoir characterization using downhole microseismic monitoring. Geophysics 75 (5) 75A129–75A137. Mcmechan, G.A., 1982. Determination of source parameters by wavefield extrapolation. Geophys. J. Int. 71 (3), 613–628. Meng, C., Wang, H., 2018. A finite element and finite difference mixed approach for modeling fault rupture and ground motion. Comput. Geosci. 113, 54–69. Meng, K., Li, J., Li, Y., 2016. Noise suppression in the dual-tree complex wavelet domain for seismic signal. J. Petrol. Explor. Prod. Technol. 7 (2), 353–359. Miyazawa, M., Venkataraman, A., Snieder, R., Payne, A., 2008. Analysis of microearthquake data at Cold Lake and its applications to reservoir monitoring. Geophysics 73 (3), O15–O21. Mousavi, S.M., Langston, C.A., Horton, S.P., 2016. Automatic microseismic denoising and onset detection using the synchrosqueezed continuous wavelet transform. Geophysics 81 (4), V341–V355. Nakata, N., Beroza, G.C., 2016. Reverse time migration for microseismic sources using the geometric mean as an imaging condition. Geophysics 81 (2). Rutledge, J.T., Phillips, W.S., Mayerhofer, M.J., 2004. Faulting induced by forced fluid injection and fluid flow forced by faulting: an interpretation of hydraulic-fracture microseismicity, Carthage Cotton valley gas field. Bull. Seismol. Soc. Am. 94 (5), 1817–1830. Selesnick, I.W., Baraniuk, R.G., Kingsbury, N.C., 2005. The dual-tree complex wavelet transform. IEEE Signal Process. Mag. 22 (6), 123–151.
5.3. Influence of source mechanisms on imaged focuses Based on equations (5) and (7), we use the integration of L2 norm of the reversed wavefield along the time axis to image the focus. Therefore, the reversed imaging based on equations (5) and (7), where the multiplication of imaged focuses at different levels does not change the physical meaning, shows the averaged radiation pattern of the located microseismic source if the wavefield around the event is perfectly sampled and the three-component data is pure without noise. Fig. 23 shows the imaged focuses for horizontal, vertical and 45° point forces in a simple 2D homogenous model with a perfect recording geometry. The configuration of receiver arrays (blue triangles) and source (red circle) is shown in Fig. 23a. It is clear that the averaged radiation patterns for the located event with various source mechanisms are perfectly presented. However, in real cases, sparse recording geometries are frequently deployed due to the limitations of field conditions. As a result, the characteristics of P- or S-wave radiation pattern may not be sampled by the sparse recording geometries. For example, the horizontal nodal plane in the P-wave radiation pattern of a vertical single point force is not sampled by the surface arrays. Therefore, the shape of the imaged focus, using the proposed method in this paper, is a system response of the recording geometry and source mechanism. 6. Conclusions The images of surface micro-seismic measurements may give unreliable location results in extremely strong ambient noise environments. This paper has shown, with a synthetic study of multiple microseismic sources in 3D layered velocity model and a field study during hydraulic fracturing in southwest of China, that conventional time reversed imaging (TRI) techniques fail to locate the events in a situation of extremely low S/N ratio. By contrast, the multi-scale TRI technique proposed in this paper minimizes the interference of noise and identifies effective events in multiple time-frequency domains. Time reversed images of effective components substantially improved the quality of the final image of subsurface microseismic sources with much sharper focuses. On the other hand, the images of noise components may reveal the velocity interface which is helpful in event location. In addition, the noise components recognized by the multi-scale TRI can be applied to estimate the background noise level. Multi-scale TRI effectively utilizes signals scattered in multiple timefrequency domains to construct a source image, therefore, improving the ability to converge the events of weak energy. However, it requires solving the wave equation independently for each level, which can be computationally demanding. Nevertheless, it is naturally suitable for parallel computation to improve the efficiency. A detailed analysis to estimate elaborated velocity model using multi-scale TRI will be an interesting research direction in future. Acknowledgements We would like to thank Xi’an Changqing Technology Engineering Co.,Ltd. for providing the field data. The authors also appreciate the technical support and help from Dr. Quanfeng Wang, Chengdu University of Technology. Appendix A. Supplementary data Supplementary data to this article can be found online at https:// doi.org/10.1016/j.petrol.2018.11.015. 159
Journal of Petroleum Science and Engineering 174 (2019) 144–160
M. Li et al.
Wang, H., Tao, G., 2011. Wavefield simulation and data-acquisition-scheme analysis for LWD acoustic tools in very slow formations. Geophysics 76 (3), 59–68. Wang, H., Li, M., Shang, X.F., 2016a. Current developments on micro-seismic data processing. J. Nat. Gas Sci. Eng. 32, 521–537. Wang, J.J., Li, Y., Liu, W.R., Xiao-Hong, X.U., Mathematics, S.O., 2016b. Dual-tree complex wavelet domain bivariate method for seismic signal random noise attenuation. Chin. J. Geophys. 59 (8), 3046–3055. Wilks, M., Wuestefeld, A., Thomas, P., Kolltveit, E., Oye, V., 2017. Initial Results from the Field on the Use of DAS as a Viable Microseismic Monitoring System of CCS Sites. EAGE Research Workshop. Witten, B., Artman, B., 2010. Signal-to-noise estimates of time-reverse images. Geophysics 76 (2), 2176–2180. Yomogida, K., 1994. Detection of anomalous seismic phases by the wavelet transform. Geophys. J. Roy. Astron. Soc. 116 (1), 119–130.
Tenger, B., Sheng, B., Yu, L., Yang, Y., Zhang, W., Tao, C., Xi, B., Zhang, Q., Bao, F., Qin, J., 2017. Mechanisms of shale gas generation and accumulation in the ordovician Wufeng-Longmaxi Formation, Sichuan basin, SW China. Petrol. Explor. Dev. 44 (1), 69–78. Zhu, Tieyuan, 2014. Time-reverse modelling of acoustic wave propagation in attenuating media. Geophys. J. Int. 197 (1), 483–494. Waldhauser, F., Ellsworth, W.L., 2000. A double-difference earthquake location algorithm: method and application to the Northern Hayward Fault, California. Bull. Seismol. Soc. Am. 90, 1353–1368. Wang, H., Tao, G., Fehler, M.C., 2015. Investigation of the high-frequency wavefield of an off center monopole acoustic Logging-While-Drilling tool. Geophysics 80 (4), D329–D341. Wang, H., Tao, G., Zhang, K., 2013. Wavefield simulation and analysis with the finiteelement method for acoustic logging while drilling in horizontal and deviated wells. Geophysics 78 (6), D525–D543.
160