Development of a new high resolution waveform migration location method and its applications to induced seismicity

1.   Introduction

To monitor induced seismicity, it is essential that seismic events are detected and located accurately. Seismic monitoring has been employed widely in unconventional oil and gas development (e.g. Rutledge et al., 1998; Lakings et al., 2006; Le Calvez et al., 2007; Gharti et al., 2010; Maxwell et al., 2010; Grigoli et al., 2013), geothermal development (e.g. Phillips et al., 2002; Dyer et al., 2008), and mining (e.g. Ge MC, 2005; Tang LZ et al., 2006; Qian JW et al., 2018). These activities can induce considerable seismicity when they create new fractures or reactivate existing fractures/faults. The locations of these induced seismic events can be used to delineate fracture distributions, evaluate stimulated reservoir volume, and identify active faults; such information is essential to improve understanding of the hydraulic fracturing used in unconventional resources development, as well as in forecasting potential hazards in mining (Maxwell and Urbancic, 2002; Rutledge and Phillips, 2003; Maxwell, 2011; Close et al., 2012; Li JL et al., 2013; Grechka et al., 2015).

Conventionally, seismic event location methods used to locate earthquakes have relied on picking the P- and S-wave arrival times manually or automatically (e.g. Douglas, 1967; Pujol, 1992; Waldhauser and Ellsworth, 2000; Zhang HJ and Thurber, 2003; Waldhauser and Schaff, 2008; Zhang HJ et al., 2010; Kwiatek et al., 2013). When seismic signals have low signal-to-noise ratios (SNRs), picking their phase arrivals reliably at different stations — whether manually or automatically — has been challenging, even in the case of downhole microseismic monitoring (e.g. Kersey, 2000; Maxwell et al., 2010; Eisner et al., 2010; Lin Y and Zhang HJ, 2016). For this reason, waveform-based location methods have been proposed to avoid having to pick phase arrivals (e.g., Kao H and Shan SJ, 2004; Pesicek et al., 2014; Zhebel and Eisner, 2015; Beskardes et al., 2018). For a complete review of waveform-based location methods and their applications, refer to Li L et al. (2020). The two principal categories of waveform-based location methods are partial waveform stacking and time reverse imaging (Li L et al., 2020).

Partial waveform stacking sums various waveform attributes that are corrected by the corresponding travel times from potential sources to different receivers (Li L et al., 2020). The trial location associated with the maximum stacking value is assumed to be the hypocenter. Most methods based on partial waveform stacking use the concept of Kirchhoff migration. By employing more waveform information than conventional arrival-based location methods, partial waveform migration-based methods can be more suitable for locating weak seismic events and may be able to improve location accuracy. By stacking actual waveform amplitudes, Duncan (2005) first introduced the waveform migration-based method to locate microseismic events induced by hydraulic fracturing of shale gas reservoirs. However, due to polarity reversals caused by shear-slip-dominated source focal mechanisms, the maximum brightness value that results from stacking actual waveform amplitudes may not correspond to the true event location because of waveform destruction. To reduce the influence of waveform polarity, Kao H and Shan SJ (2007) proposed a source-scanning algorithm (SSA) that searches for the origin time and location simultaneously by maximizing the brightness function that sums the waveform envelopes at all stations. Alternatively, to account for waveform polarity reversals, Grigoli et al. (2013) proposed stacking the short-term average to long-term average ratio (STA/LTA) series of the waveforms, and Liang CT et al. (2016) proposed stacking the focal mechanism corrected waveforms in the joint source scanning algorithm (jSSA).

By stacking different waveform attributes, the maximum brightness spot with largest stacking value can be treated as the most likely seismic event location. However, if the receivers are not well distributed and the event waveforms are greatly contaminated by noise, the brightness functions from stacking different waveform attributes generally are not focused, causing large uncertainties in the event location.

Another type of waveform-based location method is time reverse imaging, which is based on the time reversal theory and technique (Fink, 1999; Li L et al., 2020). Time reverse imaging involves backpropagating the recorded full waveforms at receivers in time to focus the source energy, by taking advantage of the reciprocity of the elastic wave equation. Different location methods based on time reverse imaging have been proposed with different imaging conditions (Gajewski and Tessmer, 2005; Steiner et al., 2008; Artman et al., 2010; Sun JZ et al., 2015; Nakata and Beroza, 2016; Zhu TY et al., 2019). Compared to the conventional arithmetic-mean imaging condition, Nakata and Beroza (2016) proposed a geometric-mean imaging condition by multiplying backpropagated receiver wavefields to increase the spatial resolution of the source location image, which is highly dependent on the accuracy of the velocity model. Sun JZ et al. (2015) proposed a hybrid multiplicative time reverse imaging method that divides the receivers into groups; arithmetic-mean imaging and geometric-mean imaging conditions are used within and between groups, respectively. The application of hybrid multiplicative time reverse imaging methods to field data has shown to be effective in imaging the spatial and temporal distribution of microseismic events (Zhu TY et al., 2019).

In this study, we adapt the hybrid multiplicative imaging condition used for time reverse imaging to a partial waveform migration-based location method. We first introduce the method by dividing receivers into groups and multiplying the summed waveform series between groups. We then compare the hybrid multiplicative imaging condition with the conventional imaging condition that stacks different waveform attributes including the actual waveforms (linear stacking), the waveform envelopes (envelope stacking), and the STA/LTA series of the waveforms (STA/LTA stacking). Finally, we test these methods using both synthetic and real induced seismicity datasets related to conventional and to unconventional gas development.

2.   Migration-Based Seismic Location Methods

Seismic waveform migration-based location methods are derived from the widely used concept of seismic migration in seismic exploration, which involves backpropagating through the medium the received wavefields at different receivers. When seismic records are backpropagated according to the seismic source excitation time, they should be focused around the source location. Several types of seismic migration methods are available (Li L et al., 2020). For simplicity, we use the Kirchhoff-type migration by stacking different waveform attributes that are corrected in time according to different travel times from potential source to different receivers.

The migration-based location method generally involves the following steps: scanning potential source points (x, y, z) and origin time $\tau$, correcting and shifting the waveforms at different receivers by corresponding travel times, then stacking different time-shifted waveform attributes to calculate the brightness functions at each point. If one point is the true event location, the shifted waveforms will be aligned well and the stacking value would be largest or brightest. Otherwise, the stacking value would be relatively small because of destruction of the waveforms. With various brightness functions, the migration-based location method has different accuracy, resolution, and anti-noise ability. For the conventional imaging condition, at each potential point $s$, all the shifted waveforms or attributes are directly summed. At the origin time $\tau$, the brightness function can be defined as,

where $N$ is the number of receivers, ${t}_{\mathrm{p}}$ and ${t}_{\mathrm{s}}$ are P- and S-wave travel times from source point $s$ to $i$th receiver, respectively, and ${l}_{\mathrm{p}}$ and ${l}_{\mathrm{s}}$ represent the lengths of the time windows for P and S, respectively.

For linear stacking, ${u}_{i}^{p}$ and ${u}_{i}^{s}$ are the actual P and S waveforms recorded at $i$th receiver. But due to non-isotropic source focal mechanisms, waveforms may have reversed polarities at different receivers. As a result, linear stacking can result in a brightness function value that is very low or even near zero at the true source location. One way to avoid this issue is to stack the waveform envelopes for different receivers. In this case, ${u}_{i}^{p}$ and ${u}_{i}^{s}$ in Equation (1) represent the waveform envelopes (Gharti et al., 2010), which are computed via the Hilbert transform (Figure 1). Compared with the linear stacking using original waveforms, stacking the waveform envelopes can give a more stable location image but with relatively low resolution.

Figure  1.  Example of a real waveform and its different waveform attributes including envelope, STA/LTA, and characteristic function. Red dashed line denotes P-wave arrival; blue dashed line denotes S-wave arrival.

DownLoad: Full-Size Img PowerPoint

Another way to avoid the waveform polarity issue is to stack the ratio of short-term average (STA) to long-term average (LTA), or STA/LTA series (Allen, 1982). STA/LTA has been widely used for earthquake monitoring in weak-motion seismology to detect events and pick first arrival times. At the $j$th time sample, the STA/LTA for a continuous waveform record can be calculated as follows:

where $j$ is the time index, ${n}_{s}$ is the length of the short-term window, and ${n}_{l}$ is the length of long-term window. The STA/LTA is sensitive to the amplitude change, and thus can be used to detect seismic events and pick phase arrivals (Figure 1). By replacing the term ${u}_{i}$ in Equation (1) with the STA/LTA, all positive time series will be stacked to avoid the waveform polarity issue.

To improve location accuracy and resolution, in this study we adapt the hybrid multiplicative imaging condition to the migration-based location method. This method divides receivers into groups and multiplies instead of summing the waveform attributes within groups. To better characterize amplitude and phase changes in the seismic waveform, we adopt the characteristic function $C$ for each trace as follows,

where x(i) represents the original seismic waveform. Compared to the original waveform, the characteristic function is more sensitive to changes in either amplitude or phase (Figure 1). Furthermore, $C$ is always positive, thus avoiding the potential polarity reversal effect of the source radiation pattern. Instead of stacking the characteristic functions at all receivers, we divide the receivers into several groups. Within the same group, the characteristic functions of each trace are added together, but between groups the sums from individual groups are multiplied. For the grid point $s$ in the image domain, at the origin time $\tau$, the brightness function is defined as,

where $m$ is the number of the groups, ${n}_{j}$ represents the number of receivers in the $j$th group, and ${t}_{\mathrm{p}}$ and ${t}_{\mathrm{s}}$ are P- and S-wave travel times from $s$ to the $i$th receiver, respectively. ${l}_{\mathrm{p}}$, and ${l}_{\mathrm{s}}$ represent the lengths of the time windows for the P and S waves. E is the envelope of the characteristic function C. W is the quality weighting factor related to the SNR of each trace. W is close to 1 if the SNR is high; otherwise W is near 0.

As illustrated in Figure 1, different waveform attributes corresponding to the four stacking methods have different sensitivities to phase arrivals. Compared to the original waveform and its envelope, both STA/LTA and the characteristic function have clearer phase arrivals. However, STA/LTA could distort actual phase arrival times.

3.   Synthetic Test

We first use the synthetic data to test different migration-based location methods with different stacking functions to compare the location accuracy and robustness with noise. We adopt the microseismic monitoring geometry shown in Figure 2, which is based on a real induced seismicity monitoring array that consists of five boreholes, with eight geophones in each borehole (Zhang HJ et al., 2009; Li JL et al., 2011). However, due to the instrumentation issue for some receivers, in each borehole there are only three to eight receivers working normally. As a result, we use just 27 receivers to simulate the synthetic data. The one-dimensional (1D) velocity model used in the synthetic test is derived from the sonic logs (Figure 3). To calculate the synthetic data, we use the discrete wavenumber method (Bouchon, 1981) to generate the theoretical seismograms.

Figure  2.  Microseismic monitoring geometry used in the synthetic test. The blue triangles denote the locations of borehole receivers in five monitoring wells.

DownLoad: Full-Size Img PowerPoint

Figure  3.  One-dimensional (red line) P- and (blue line) S-wave velocity models used in the synthetic test.

DownLoad: Full-Size Img PowerPoint

We set the seismic source at (X = 0, Y = 0, Z = 1.5) km and superimpose different levels of noise (10%, 50% and 200%) on calculated seismograms to test the abilities of different migration-based location methods in the noisy environment. Figure 4 shows the synthetic waveforms with 10%, 50%, and 200% noise, which correspond to SNRs of 10, 2 and 0.5, respectively. The blue lines denote the noisy seismograms added with various levels of noise; the red lines represent the waveforms after bandpass filtering of 10−35 Hz.

Figure  4.  Synthetic waveforms with (a) 10% noise, (b) 50% noise, (c) 200% noise. The blue lines denote the seismograms added with various levels of noise; the red lines represent the waveforms after bandpass filtering of 10−35 Hz.

DownLoad: Full-Size Img PowerPoint

Next, we apply four different migration-based location methods to these noisy waveforms. For the different stacking methods, the stacking window length depends on the length of one cycle of P-waves and S-waves in different waveform attributes. For location methods with linear stacking and envelope stacking, ${l}_{\mathrm{p}}$ is set as 0.035 s and ${l}_{\mathrm{s}}$ is set as 0.075 s, equivalent to 35 and 75 sample points, respectively. For the STA/LTA stacking, we choose the short-term window length to be 50 points, which is determined by the length of the P-wave; the long-term window length is chosen to be 200 points (4 times the length of short-term window) to calculate the STA/LTA series; the stacking window length is 50 points for both P- and S-waves. For the new location method with hybrid multiplicative imaging condition, we treat the receivers in each well as one group, and thus we divide the 27 receivers into 5 groups. The lengths of the stacking windows for P-wave (${l}_{\mathrm{p}}$) and S-wave (${l}_{\mathrm{s}}$) are 30 and 40 points, respectively. We set the quality parameter $W$ of each receiver equal to 1. Considering the location of the synthetic seismic source to be (0, 0, 1.5) km and its origin time to be 1 s, we search X in the range of [−1.1, 1] km, Y in the range of [−2.5, 2.5] km, and Z in the range of [0, 2.5] km with a grid interval of 0.05 km, and $\tau$ from 0.5 to 1.5 s with a time interval of 0.1 s. In total, the number of searching points is 43 by 101 by 81 by 11, which is the same as the real data. Figure 5 shows the comparison of locations for the four stacking location methods. The highest brightness point is treated as the source location for the four location methods. If the image is more focused, then the location is more accurate with smaller uncertainty. In the cases of SNRs equal to 10 and 2, the stacking images based on linear stacking, envelope stacking, and STA/LTA stacking clearly show the footprints of observation systems; the latter two have more focused images. By comparison, the stacking images produced by our proposed method are highly focused around the true location. With the noise level reaches 200%, the images from linear stacking are poorly focused; envelope and STA/LTA stacking methods slightly improve the accuracy; all three of these location methods are greatly affected by high levels of noise, producing poorly determined event locations when noise is high. The new location method we propose, however, remains able to image the event with high resolution even in high noise situations, reliably determining its location. When the SNR is set at 0.5, the new location method with the hybrid multiplicative imaging condition has the highest stacking value (> 0.5), while the linear stacking value is below 0.25 and envelope and STA/LTA stacking values are even lower (Figure 6).

Figure  5.  Location images in the XZ and XY planes from the four stacking methods with different SNRs: (a) SNR = 10; (b) SNR = 2; (c) SNR = 0.5. The highest brightness point is treated as the source location.

DownLoad: Full-Size Img PowerPoint

Figure  6.  The variation of stacking values at the true location along with times for the four stacking methods in cases of different waveform SNRs: (a) linear stacking; (b) envelope stacking; (c) STA/LTA stacking; (d) new method.

DownLoad: Full-Size Img PowerPoint

A further test of the new proposed method is to assess its ability to locate events that are close in time and space. For the spatial resolution test, we set two sources closely located — at (0, 0, 1.25) km and (0, 0, 1.5) km — with the same origin time, and use the same searching grid range as before. The location results from the different stacking methods are shown in Figure 7. It can be seen that none of the three conventional stacking methods correctly resolve the two events, but the new stacking method based on hybrid multiplicative imaging condition is able to do so. For the temporal resolution test, the two events are set at the same location but with close origin times of 1 s and 1.2 s. From the stacking results, it can be seen that envelope stacking and STA/LTA stacking detect only one event, while linear stacking and new method correctly detect two time-separated events (Figure 8). The new method gives the higher stacking values at the corresponding origin times. These two synthetic tests show that the new method has superior ability to locate events that are close in space and time.

Figure  7.  Comparison, in the case of two sources, of stacking values in the vertical direction at the horizontal location of (0, 0) generated by the four stacking methods. The black dashed lines denote the exact depths (1.25 km and 1.5 km, respectively) of the two sources.

DownLoad: Full-Size Img PowerPoint

Figure  8.  Comparison of temporal stacking value variations at the true source location for four stacking methods in the case of two sources located in the same location but separated in time. The black dashed lines denote the origin times of 1 s and 1.2 s for the two sources.

DownLoad: Full-Size Img PowerPoint

4.   Real Data Applications

In this section, the four stacking methods are applied to two real datasets. One is the borehole microseismic monitoring dataset from an oil field in Oman; the other is the induced seismicity monitoring dataset from a surface seismic network in the Changning shale gas development field of the Sichuan basin, China.

4.1   Application to the Microseismic Dataset from an Oilfield in Oman

To monitor induced seismicity caused by oil/gas production in an oilfield in Oman, a borehole network consisting of five closely spaced monitoring wells was established in the most seismically active part of the reservoir; each monitoring well has multiple three-component receivers with depths ranging from 750 to 1250 m (Figure 2; Zhang HJ et al., 2009; Li JL et al., 2011). From the data recorded by those borehole receivers, more than 15800 events have been identified and nearly 5400 events are located and analyzed for reservoir monitoring using an arrival-based location method (Figure 9; Zhang HJ et al., 2009). It can be seen that most seismic events were distributed along the two major faults, oriented in the NE-SW direction. We selected 15 events from each of the two major faults and labeled them as Group A and Group B (Figure 9). In general, the data in Group A have higher quality than those in Group B. For each group, we apply bandpass filtering of 10−35 Hz to the waveform recordings. In this real data application, the searching grid is the same as that used for the synthetic test, i.e. [−1.1, 1] km in X, [−2.5, 2.5] km in Y, and [0, 2.5] km in Z with a grid interval of 0.05 km, and [0.5, 1.5] s in $\tau$ with a time interval of 0.025 s. For linear stacking and envelope stacking, ${l}_{\mathrm{p}}$ is 0.125 s and ${l}_{\mathrm{s}}$ is 0.25 s, equivalent to 250 and 500 samples, respectively. For the STA/LTA stacking, we choose the stacking window length to be 100 points for both P- and S-waves, and the length of the short-term window to be 100 points and the length of the long-term window to be 400 points to calculate the STA/LTA series. For the new stacking method, we treat the receivers in each well as one group so that we have five groups of receivers. The length of the stacking windows for P- and S-waves (${l}_{\mathrm{p}}$ and ${l}_{\mathrm{s}}$) is set at 200 points.

Figure  9.  Distribution of events (brown dots) and monitoring boreholes (blue triangles). Red dots denote the selected events used for the analysis of stacking methods. Black lines indicate fault traces in the study area.

DownLoad: Full-Size Img PowerPoint

We select an induced event as an example to show the capability of the new stacking method over other stacking methods. The data quality for most receivers is high; only the data recorded by receivers in the well YE have relatively poor quality (Figure 10). For this reason, the quality factor W is set as 1 for the receivers in wells YA, YB, YC, and YD, but 0.5 for the receivers in well YE. The stacking images from the four stacking methods applied to this event are shown in Figure 11. It can be seen that linear stacking yields no clear focusing points in the XY and XZ images. The envelope stacking method is slightly better; one clear focusing point can be seen in the XY image. Compared to envelope stacking, the focusing point in the XY image from STA/LTA stacking is more concentrated. For both envelope and STA/LTA stacking, there are apparently two focusing points in the XZ image. In comparison, the new stacking method, based on our hybrid multiplicative imaging condition, greatly improves the focusing effect and yields a single focusing point in the XZ image by eliminating potential artifacts.

Figure  10.  The vertical component waveforms of an induced event recorded by different receivers of the borehole monitoring network. The time sampling interval is 0.5 milliseconds.

DownLoad: Full-Size Img PowerPoint

Figure  11.  Comparison of stacking images in the XZ and XY planes from four stacking methods. The highest brightness point is treated as the source location.

DownLoad: Full-Size Img PowerPoint

We further systematically compare, for 30 selected events, the locations obtained from the four stacking location methods to the reference locations for these events, which we obtain from the double-difference tomography study of Zhang HJ et al. (2009) using 1999 events. From the statistical analysis (Table 1), it can be seen that locations for these 30 events from the new stacking location method are closest to the reference locations, indicating that it is most accurate among the four stacking methods and can provide locations comparable in accuracy to those obtained from arrival-based location methods.

Table  1.  The mean and standard deviation values of location differences between four stacking methods and reference locations from Zhang HJ et al. (2009).

Stacking methods	Mean (km)	Standard deviation (km)
Linear stacking	0.2834	0.2536
Envelope stacking	0.1826	0.1329
STA/LTA stacking	0.1289	0.1035
New method	0.1220	0.0931

 | Show Table

DownLoad: CSV

5.   Application to Induced Seismicity in Changning, Sichuan

Since 2014, earthquake activity has increased abnormally in the southern Sichuan Basin due to the shale gas development (Lei XL et al., 2017; Tan YY et al., 2020). On June 17, 2019, a M6.0 earthquake struck the Changning County, Sichuan, followed by a series of earthquakes occurring in the surrounding area. Using a CNN based event detection method (Yang SB et al., 2021), 8697 earthquakes were detected from continuous waveform recordings by a local seismic array consisting of nine stations (Figure 12). We applied the new stacking location method to the waveform segments consisting of detected events. The nine stations are divided into 3 groups according to their spatial distribution (Figure 12). The 1D V_p and V_s models in the study area are shown in Figure 13. We first test two events with waveforms shown in Figure 14, which have relatively high SNRs after 1 Hz high-pass filtering. Event 1 is located within the station array, but Event 2 has poor station coverage and is located in the northern edge of the station array (Figure 12). The imaging domain is [0, 70] km in X, [0, 70] km in Y, and [−2, 8] km in Z with a grid interval of 0.02 km. For the origin time $\tau$, the search range is [0, 25] with an interval of 0.1 s. The length of stacking windows for P- and S-waves, ${l}_{\mathrm{p}}$ and ${l}_{\mathrm{s}}$, is 50 points (0.5 s).

Figure  12.  Distribution of 8697 events (black dots) located by the proposed new stacking location method. The blue triangles represent stations used in this study. The red solid lines denote the faults and the blue dashed lines denote the anticlines. The red star represents the location of the M6.0 earthquake occurring on June 17, 2019.

DownLoad: Full-Size Img PowerPoint

Figure  13.  1D P- and S-wave velocity models used for location.

DownLoad: Full-Size Img PowerPoint

Figure  14.  Waveforms of two induced seismic events. (a) waveforms of Event 1 located within the station array; (b) waveforms of Event 2 located in the northern edge of the station array. The locations of the two events are shown as white stars in Figure 12.

DownLoad: Full-Size Img PowerPoint

For Event 1, the stacking images clearly show focusing points in the XY, XZ, and YZ planes (Figure 15a), indicating that the event is well located. In comparison, for Event 2 the stacking image is focused only in the XY plane; it is not well focused in the vertical direction, indicating that it is not well constrained in depth. This is mainly because Event 1 is better covered by the stations; Event 2 has poor station azimuthal coverage. Considering the poor constraint on event depth for events outside the station array, we test the two-dimensional (2D) stacking by fixing the event depth at 4 km (Figure 15b). It can be seen that, in the horizontal plane, the image from 2D stacking is very similar to the three-dimensional (3D) stacking image, indicating that the horizontal locations of events are relatively well constrained. We also test the linear stacking method on both Events 1 and 2. The linear stacking method gives comparatively poorly focused location images, especially in the vertical direction, even for Event 1 which was well covered by the stations (Figure 15c).

Figure  15.  Comparison of stacking images from linear stacking and the new stacking method based on hybrid multiplicative imaging condition, for (left) Event 1 and (right) Event 2. (a) 3D stacking image from the new stacking method; (b) 2D stacking image from the new stacking method; (c) 3D stacking image from the linear stacking method. White triangles denote seismic stations.

DownLoad: Full-Size Img PowerPoint

Based on the above tests, to determine the horizontal locations of 8697 detected earthquakes, we fix the event depths at 4 km. Figure 12 shows the location results of these 8697 earthquakes. It can be seen that the distribution of induced earthquakes in June, 2019 is similar to the results from other recent studies (Lei XL et al., 2019; Long F et al., 2020; Jia K et al., 2020). Earthquakes are distributed primarily in two blocks: Block A corresponds to the Changning salt mine area and Block B corresponds to the shale gas development area. The events in Block A are distributed mainly along the Changning–Shuanghe anticlines; most of them are aftershocks of the Changning M6.0 earthquake that occurred on June 17, 2019. In Block B, the seismic events are distributed in clusters around the shale gas well pads within a very short time window (one to two months), suggesting that these events were more likely induced by nearby hydraulic fracturing operations intended to produce the shale gas.

In summary, based on the real data test results, we show that the new stacking method based on hybrid multiplicative imaging condition can better suppress noise and eliminate artifacts in the source location image than other stacking methods. From its application to data from the induced seismicity due to conventional oil/gas production in Oman, it is clear that the new stacking method can more reliably and accurately locate seismic events and produce event locations closer to the arrival-based location method than the other three more conventional methods. For the case in the Changning shale gas development area, even when the number of receivers is small and some events are located outside the station array, the new stacking method can still produce relatively well constrained horizontal event locations.

6.   Discussion and Conclusions

Locating seismic events is a crucial step for earthquake monitoring. Conventional location methods use seismic phase arrival times. However, when seismic waveforms exhibit low SNR or the first arrivals are emergent, it is a challenge to pick first arrivals accurately and efficiently. In comparison, migration-based methods do not need to pick the first arrivals and can determine event locations by stacking the waveforms and searching for the brightest spot in the stacking image.

In this study, we have proposed a new waveform stacking method based on a characteristic function and a hybrid multiplicative imaging condition. The proposed method first divides the stations into groups according to their distances from each other, then sums the characteristic functions of traces within each group and multiplies the summation series among groups (instead of simple summation) to improve the focusing of the stacking image.

Compared with linear, envelope, and STA/LTA stacking approaches, the new method can produce stacking images of significantly higher resolution, thus resulting in the improved determination of locations. Linear stacking images can be degraded by polarity reversals among the original waveforms. In comparison, envelope stacking and STA/LTA stacking avoid the issue of waveform polarity reversals and provide relatively high-quality event locations when waveform SNRs are relatively high. However, in cases of low SNRs and/or poor station coverage, envelope stacking and STA/LTA stacking may deliver event locations with large uncertainties.

We have applied the new location method to two real induced-seismicity datasets that are related to conventional oil/gas production in an oilfield in Oman and unconventional shale gas development in southern Sichuan basin, China. For 30 selected induced seismic events in the Oman oilfield that were monitored by downhole receivers in five boreholes, the locations derived from the proposed new stacking location method are comparable to double-difference locations determined by using first arrivals. For the induced seismicity in the southern Sichuan basin that were monitored by a surface seismic array consisting of nine stations, the new stacking method has provided high-quality event locations in the horizontal plane. Both the synthetic and real data tests suggest that our proposed stacking location method can locate seismic sources reliably without picking their phase arrivals, which indicates that this new method is of potential use in real-time location of earthquake epicenters.

Acknowledgments

This study is supported by National Key R&D Program of China under grant 2018YFC1504102 and the National Natural Science Foundation of China under grant 41961134001.