- Email: [email protected]
An approach to automatic location of regional events
Physics of the Earth and Planetary Interiors 113 Ž1999. 293–302
An approach to automatic location of regional events V.I. Pinsky
)
Geophysical Institute of Israel, PO Box 2286, Holon 58122, Israel Received 18 May 1998; received in revised form 9 July 1998; accepted 9 July 1998
Abstract The conventional network location algorithms for automatic processing are based on P and S first arrivals. These procedures have the following shortcomings: Ž1. it is often difficult to identify the proper phase: ŽPg or Pn, Sg or Sn.; and Ž2. first arrivals are often masked by noise. Both factors may cause significant location errors. An alternative is to look for maximums of seismic energy time–distance distribution, which are less sensitive to these factors. We measure the time-of-maximum of P and S waves envelopes vs. distance for each available station of the Israel Seismic Network ŽISN., thus providing a travel time curve ŽTTC.. The record envelopes are obtained using 1D x 2 optimal detector and 3–6 Hz ‘short-time-average’ time-curves, having enhanced sensitivity for seismic signal arrivals. The corresponding P and S time-of-maximum vs. distance functions are approximated linearly by a least-square method for a set of local earthquakes and quarry blasts. Travel time inversion for location of small events comprises three main steps: Ž1. cleaning of the records from noise bursts and computation of the envelopes, Ž2. triggering and identification of the P and S phases and computation of their energy maximums, and Ž3. maximization of a sum of station residuals as a function of epicenter coordinates. The maximum of the functional is looked for on a 60 = 60 km grid, step of 2 km, covering different parts of Israel and Jordan. The algorithm is not sensitive to the source depth, thus providing epicenter determination only, but shows to be convenient and robust. As a result of the preliminary study for a set of relatively weak 74 local earthquakes and 58 quarry blasts, M L ; 1.5–2.5, we have obtained the accuracy of epicenter estimation "6 km for 80–90% of both types of events, which is satisfactory for automatic location. The accuracy is measured relative to the ISN bulletin locations for earthquakes and the ground truth information for quarries, respectively. q 1999 Elsevier Science B.V. All rights reserved. Keywords: Seismic event location; Envelopes; Travel time curves; Grid search optimization
1. Introduction Routine source location in seismology requires estimation from observed signal arrivals of four coordinates; geographical coordinates of the epicenter Žprojection of the point-source on the surface., depth
)
Fax: q972-3-5502925; E-mail: [email protected]
and source time. For regional events, an analyst marks onsets of Pg, Pn, Sg, Sn phases on the seismogram. These onset times are then converted to the corresponding source coordinates and time through the given travel time model, often determined by the velocity structure. For small and remote events with slow and gradual onsets, the accurate estimation of the first arrival time is problematic due to low signal-to-noise ratio. Thus, more reliable alternative is required. Sometimes, the depth can be considered known, e.g., when the source is classified as a quarry
0031-9201r99r$ - see front matter q 1999 Elsevier Science B.V. All rights reserved. PII: S 0 0 3 1 - 9 2 0 1 Ž 9 9 . 0 0 0 0 8 - 4
294
V.I. Pinskyr Physics of the Earth and Planetary Interiors 113 (1999) 293–302
blast or some other surface event. Depth knowledge is also not critical on a stage of preliminary location. In these cases, a possible solution is to utilize seismic energy distribution in time, space and frequency domains. The physical basis is that energy of the P and S seismic phases obey diffusion law and propagates with constant velocity in the stratified media. On this basis, Ryzhikov et al. Ž1995. derived theoretically a functional, which has a maximum, over the observed envelopes of the seismic phases, pointing to the source epicenter. Shoubik et al. Ž1996. applied principles of the so-called ‘seismic tomography’ based on the calculation of semblance and ‘signalto-noise ratio’ functionals over the envelopes of the seismograms. By the method, using a grid search approach, Shoubik et al. satisfactory estimated source epicenters of the 19 Galilee earthquakes, recorded by the Israel Seismic Network ŽISN.. Husebye et al. Ž1998. applied the short-time-average ŽSTA. procedure for calculation of the signal envelopes with the 3–6 Hz band-pass prefiltering the data of several seismic networks, and showed its effectiveness for epicenter determinations. In this paper, we use a similar approach focusing on the estimation of the seismic phases energy travel time curves ŽTTC.. We measure TTC for the maximums of the seismogram envelopes in the time intervals corresponding to P and S waves and use them for location instead of the usual P and S TTC of the first arrivals. For location, we use envelopes obtained through different types of filtering and develop a simple triggering procedure to extract proper P and S phases. The P and S phases vary in a complex manner and are often misidentified. In routine analysis, this may lead to large errors in source location. We concentrate on this problem and construct a procedure robust enough to phase wrong identifications. For the analysis, we use the expanded ISN database, including earthquakes and quarry blasts from the three different areas of Israel and its surroundings, and find their epicenters using a grid search method.
0.5–12 Hz. of the ISN Ž; 30 stations. telemetered to a Hub via radio FM link. The data include 132 seismic events which: Ž74 earthquakes, 58 quarry blasts. occurred in Galilee Ž30 earthquakes, 36 quarry blasts., Gilad Ž18 earthquakes. and the Dead Sea region Ž26 earthquakes and 22 quarry blasts. with magnitudes of M L s 1.0–2.8 at distances of 15–310 km. Event epicenters, locations of stations and quarries are presented in Fig. 1.
3. Method The whole procedure may be split in three blocks; Ži. preprocessor, Žii. picker and Žiii. locator. The Preprocessor performs despiking and a band-pass filtering of the seismograms in an appropriate frequency band. The records we are operating with are often spoiled by a different kind of spikes and noise bursts, caused by radio impulses, industry and micro-seismic noise that provide false detections. Thus, cleaning of records is crucial to the success of the whole procedure, based on detecting of seismic signals. To further enhance signal-to-noise ratio, we perform prefiltering Žwhitening. of each channel seismograms Vt by the procedure:
ž
Ut s Vt y Ý a k Vtyk k
/
s 2,
k s 1, . . . , M ,
Ž 1.
where a k and s 2 are the coefficients of the autoregression time series model, adapted to the noise preceding the signal. Ut is the whitened time-series being the output of the filter. Time interval for adaptation is selected equal to 200 samples Ž4 s. from the start of the recording. The next step is to calculate and square the sum of autocovariance coefficients Cm of the Ut in a moving window of length S, equivalent to a sort of envelope: 2
Zt s Ž C0,t y 1 . q Ý Cm2 ,t ,
m,s 1, . . . , M ,
m
Cm ,t s Ý Utqsym UtqsrS,
s s 1, . . . ,S.
Ž 2.
s
2. Data The data are presented by the short period seismograms Ž50 Hz sampling rate with efficient bandwidth
This procedure is equivalent to the statistically optimal detector by Kushnir et al. Ž1990., which enhances sensitivity of amplitude and whiteness changes of the Ut time-series due to the incoming signal. For the stationary noise, Vt variance of the
V.I. Pinskyr Physics of the Earth and Planetary Interiors 113 (1999) 293–302
295
Fig. 1. Map of Israel and adjacent areas, showing the three seismic database regions: 1—Galilee, 2—Gilad and 3—Dead Sea region.
corresponding white noise Ut is C0,t f 1, and autocorrelation Cm ,t f 0 for m ) 0, hence, Zt f 0. A
more explicit statement is that the product Zt S asymptotically obeys a x 2 probability distribution
296
V.I. Pinskyr Physics of the Earth and Planetary Interiors 113 (1999) 293–302
with M q 1 degrees of freedom for S ™ `. When the moving window reaches the signal in the seismogram, the time-series Vt is no longer stationary, and, hence, Ut is no longer white, which yields the strong increase of the sequence Zt values. Even for a small number of coefficients wwe choose M s 3, S s 50 Ž1 s window.x and rather weak seismic phases, Zt exhibits abrupt fronts and sharp energy maximums. In most cases, there is a clear separation between P and S phases even at short distances where these signals are close in time. The Eq. Ž1. noise-adapted filter may be efficient also for the suppression of Rg surface waves, excited by explosions and lying in the same frequency band as micro-seismic noise. For computational efficiency, the Zt curve is decimated by 1r10. Examples of the Zt curves are shown in Fig. 2; the Zt curves have local maximums Pmax , S max , at time points Tp and Ts , which are coincidental with P and S arrivals in the original non-transformed records. When M s 0, the function Zt coincides with the usual STA trace.
The Picker automatically provides Tp time tied to Pmax by finding a local maximum of Zt in the time window wTTp ,TTp q Ip x and similarly for Ts in the window wTTs ,TTs q Is x, where TTp and TTs are the two sequential moments of triggering and Ip and Is Ž Ip s Is s 5 s. are the specified time intervals, assumed to contain signal. The first triggering is pronounced if for the majority of points inside a moving window of length LS ŽLS s 50 samples or 10 s. Zt ) TL = LTA t Žthis requirement serves to skip spikes and bursts, left after the ‘preprocessing’., where TL s TL p ŽTL p s 10., LTA t is a long-timeaverage, continuously calculated as: LTA t s Ý Zkrt
k s 1, . . . ,t.
The second triggering is pronounced if the condition above is true for t ) TTp q Ip with a new trigger level TL s TL s ŽTL s s 100.. The Locator. In Fig. 3, we display the time–distance pairs ŽTp , R ., ŽTs , R . covering the distance range
Fig. 2. Optimal detector curves in logarithmic scale for the seismograms for four ISN seismic stations. Time picks Tp and Ts correspond to the local maximums Pma x and S max of the Zt curves Ženvelopes..
V.I. Pinskyr Physics of the Earth and Planetary Interiors 113 (1999) 293–302
297
Fig. 3. Tp and Ts time picks vs. distance for Ža. 30 earthquakes, Žb. 36 quarry blasts from Galilee. Best fitting lines—travel time curves—are also shown.
0–350 km for the 30 Galilee earthquakes and the 36 Galilee explosions. Fig. 3 shows that most of the
pairs are clustered along two straight lines. Therefore, travel times Tp Ž R ., TsŽ R ., for the arrivals of
V.I. Pinskyr Physics of the Earth and Planetary Interiors 113 (1999) 293–302
298
Pmax and S max , were estimated as TaŽ R . s T0a q RrVa , a s P, S. Using a simple least squares procedure Žexcluding outliers., we estimate for earthquakes Vp ( 6.4, Vs ( 3.6 and T0p ( 2.5 s, T0s ( 3.5 s and for explosions Vp ( 6.4, Vs ( 3.5, T0p ( 3.5 s, T0s ( 6. Thus, the slopes of the travel time functions for earthquakes and explosions are close, but the intercepts are different especially for S waves. We choose for our ‘locator’ Vp s 6.4, Vs s 3.5, T0p s 2.5, T0s s 3.5. The usual procedure for location, efficient in many cases, is a least-squares technique, where solution is due to minimization of a sum of squares of the residuals Da i s Ta y TaŽ R i . y T0 , for a number of stations, where T0 is the source time. If the observations include random outliers of the time Ta estimations, this procedure is no more efficient and besides, may cause large location errors. In Fig. 3, we see that there are many outliers due to both false triggering and false phase identification. Those which are below the P line are mostly spikes and bursts, while those over the S line are low velocity Rg surface wave triggers. Many P identifications by mistake give S travel time and vice versa. In such cases, residuals appear like Dac i s Ta y Tc Ž R i . y T0 , a,c s P,S. These features should be taken into account in the specialized location scheme, since the conventional least-squares procedure is not applicable. Next, instead of one station residuals Dac i , we shall consider two station residuals Dabc i j s Dac i y D bc j , i, j s 1, . . . , N, a,b,c s P,S, which do not depend on origin time T0 , thus, reducing the number of unknown variables from three to two. On the basis of the above discussion, we suggest to perform the location of the epicenters by maximization of the following target function: F Ž X 0 ,Y0 . s Ý
Ý Wi j A Ž Dabc i j . ,
ij abc
i s 1, . . . , N y 1, j s i q 1, . . . , N, a,b,c,s S,P 4 , Ž 3. comprising 8 = N = Ž N q 1.r2 observations. From the ‘picker’ above, we receive two time picks Ta and
T b for each channel, but the ‘picker’ is not aware to which phase: P, S, or Rg do they belong. To distinguish between correct and wrong identifications, we use the Eq. Ž3. function AŽ D . s Ž< D < q ´ .y1 of deviations D Žinstead of usual AŽ D . s D 2 ., ‘ voting’ for true pickings by its large values and ‘ voting’ against false pickings by its small values, respectively. Thus, the procedure unites in one the two tasks: the location with the phase identification, so that majority of ‘good’ pickings, close to P or S track, would guarantee successful solution. The control parameter ´ ) 0 helps to avoid possible zeroes of the residuals Dabc i j , thus providing robustness of the Xmax , Ymax estimator, which maximizes F Ž X,Y .. The parameter ´ is selected automatically as ´ s 0.1 q k = 0.1, k s 1, . . . ,30 according to empirical criterion that for more than Nr2 stations, there is at least one combination of a, b, c and j that Ž< Dabc i j < ´ . - 3. The Wi j coefficients are chosen and for the time being, equal 1, though too close and too remote stations, may be given lower weight. We prefer the simple case of two independent variables to obtain X max , Ymax through a simple grid search Žusing grid of 60 km = 60 km, 2 km step., though any other non-linear optimization procedure may be employed.
4. Preliminary location Our non-linear optimization procedure requires an initial value of X 0 , Y0 , preferably as close to the true source as possible. In this paper, we deal with local sources being, at most, some tens of kilometers from the nearest station. Hence, as a simple preliminary estimate, we choose the coordinates of the station which is triggered firstly. Let us call it station N1. Of course, due to the noise bursts preceding seismic signal the estimator may appear to be too crude. A likely improvement appears feasible by using observational experience that for relatively remote stations with R G 100 km or ŽTs y Tp G 10 s., Tp and Ts are determined more reliably than for closer stations. For
Fig. 4. Map and histogram of the 30 Galilee earthquake epicenters relative to ISN bulletin locations wŽ0,0. is true epicenter of each eventx, using Tp , Ts picks, provided by Ža. the optimal detector with 0.2–10 Hz prefiltering, Žb. STArLTA detector with 3–6 Hz prefiltering; and the Žc. best of the OD and STArLTA solutions according to the maximum of the target function Ž3.. Error ellipsoids, 80% and 90%, are provided by the normal distribution fit.
V.I. Pinskyr Physics of the Earth and Planetary Interiors 113 (1999) 293–302
299
V.I. Pinskyr Physics of the Earth and Planetary Interiors 113 (1999) 293–302
300
such remote stations, distance R can be directly found from the equation: Tp y Ts s T0 p y T0 s q R Ž 1rVp y 1rVs . ,
Ž 4.
and then the source origin time T0 is given by: T0 s Ý Ta j y Ž R jrVa . rN.
Ž 5.
j
Then for each station, observed phase time is: ta j s T0 q T0 a q R jrVa .
Ž 6.
Taking station N1 and some remote station N2, which fits Eq. Ž4. best, we solve two equations Ž6.
for N1 and N2 relative to X 0 , Y0 , using an equivalent quadratic equation. From Fig. 1, we see that Israel is clustered by tens of quarries which are close to the natural seismic sources. Taking into account that 90% of seismic events in the region are quarry blasts, it is reasonable to use as preliminary estimator coordinates of the quarry, which provide the largest value to the F Ž X,Y . function. Being applied together with any reliable automatic classification procedure, such as those described in ŽShapira et al., 1996. or ŽPinsky and Shapira, 1998., this decision may be regarded as
Fig. 5. Location of the Ža. 64 earthquakes and Žb. 58 explosions from the Galilee, Gilad and the Dead Sea regions, using the best of the OD and STArLTA solutions Žas for Fig. 4.. See text for discussion.
V.I. Pinskyr Physics of the Earth and Planetary Interiors 113 (1999) 293–302
final in the case when event is classified as a quarry blast, or at least competing with the full optimization procedure. It may be also regarded as methods of event recognition, such as described in ŽJoswig, 1990. and ŽFedorenko et al., 1998.. Finally, what estimator should be chosen as initial value for the optimization procedure? The answer depends on how many stations are triggered. In most cases, we choose the one which provides the largest contribution to the F Ž X,Y . function. For small events Žless than five stations triggering., we utilize the first station ŽN1. coordinates for epicenter determination.
5. Results Our epicenter location results in Figs. 4 and 5 are 1D by histograms of deviations R from X 0 , Y0 and 2D plots of points Ž X y X 0 , Y y Y0 .. Here, X 0 , Y0 are presumed true source coordinates as reported in bulletins in the case of an earthquake and ground truth coordinates of a quarry, in case of explosions. The 80% and 90% error ellipsoids are presented on the corresponding plots. Our picking procedure, described above, is parameterized in order to ensure operational flexibility. It proved to be robust with regards to variations of time window length and trigger level but sensitive to the number of coefficients M in the adaptive filter Ž2. and the frequency band-pass ŽFBP.. An instructive example is given in Fig. 4a,b with the usage of the event analysis procedures for 30 Galilee earthquakes: Ž1. M s 3 and FBP s w0.2–10 Hzx; and Ž2. M s 0 Žthe STA case., FBP s w3–6 Hzx. The final results look promising for both cases Ž95% within 8 km from epicenter and almost 100% within 10 km., though there are some outliers. For each event, one solution from the two cases was selected, corresponding to the larger value of function Ž3. for the fixed value of ´ . The result shown on Fig. 4c improves considerably Žmore than 70% of events are within 2 km range Ž1 grid step. and more than 95% within 4.5 km.. The above analysis procedures were applied to all the events in our dataset. Fig. 5a gives the results for the 74 earthquakes and Fig. 5b—for the 58 explosions. Statistically, the results of location for both types of events look very good and appear to be
301
almost the same: 84% of events within 4.5 km radius, 96% within 7.5 km, and almost 100% within 10 km. For many explosions, the locations were better than those obtained routinely, and this probably applies to the shallow earthquake locations as well.
6. Discussion and conclusions The above results look promising for automatization of local seismicity monitoring. Our novel procedure performed well providing robust and accurate solutions in the different geological regions of Israel, Galilee, Gilad and Dead Sea. Nevertheless, there is some scatter of locations, caused probably by difference in depth of source ranging between 0 and 30 km and by heterogeneity of the wave path. Some of the earthquakes that were located with error of 10–14 km were relatively deep and thus, probably, our assumption of depth independent velocity model of wave energy propagation was invalid. The most adverse factors influencing source location are the low signal-to-noise ratio and abrupt changes of the noise character preceding the phase arrival. The corresponding results are the ‘missing of target’ and the ‘false alarms’, respectively. Both provide the ŽTa , R . ‘outliers’, but really misleading are only those which deviate too much from the P and S curves in Fig. 3. The late P and S phase identifications are mostly due to the high noise level, while too early phase pickings are due to the sudden change in the noise frequency or amplitude. The Locator ŽEq. Ž3.. eliminates the outliers Žsee explanation above. successfully if their number is relatively small, while the rest of the P and S pickings are accurate. Otherwise, the outliers may ensure large but infrequent location errors as seen on Figs. 4 and 5. Nearly half of the events are explosions which occasionally proved problematic to analysis. The reason is that this kind of signal sources cause strong excitation of the Rg surface waves, propagating with low velocities in the range of 1–3 kmrs. This, in turn, may result in too late pickings of the true P and S phases, demonstrated in Fig. 3. The Rg waves are often characterized by signal frequencies, lower than for body waves and thus, filtering in the proper frequency band is the possible way for bettering the
302
V.I. Pinskyr Physics of the Earth and Planetary Interiors 113 (1999) 293–302
P- and S-pickings. Another factor is scattering and multi-pathing due to heterogeneity of the upper-most part of the crust. We chose the 3–6 Hz frequency band suggested by Husebye et al. Ž1998., and saw that filtering by a band-pass Butterworth filter with gradual slopes performed much better. For both earthquakes and explosions, efficient despiking was important. In spite of additional difficulties, final results of location for explosions are almost the same as for earthquakes. We think that some sort of match filters or velocity filters may be helpful for further refinements.
Acknowledgements I would like to thank Eystein Huseby and Avi Shapira for suggesting this research and for helpful advice and Anton Dainty for careful revision and valuable comments. I am grateful to Yuli Zaslavsky, Yefim Gitterman and Ilya Turetsky for fruitful discussions. Many thanks to Boris Gurevich for careful reading of the text and useful comments and to Nitzan Rabinovich for the text corrections. This study was supported partially by the U.S. Department of Defense, Contract No. DSWA01-97-C-0151 and partially by the German Israel Foundation for Scientific Research and Development, Contract No. G 0463245.08r95.
References Fedorenko, Yu.V., Husebye, E.S., Heincke, B., Ruud, B.O., 1998. Recognizing explosion sites without seismogram readings: neural network analysis of envelope-transformed multi-station SP recordings 3–6 Hz. Geophys. J. Int. 133, F1–F6. Husebye, E.S., Ruude, B.O., Dainty, A.M., 1998. Fast, robust and reliable epicenter determinations: envelope processing of local network data. Bull. Seism. Soc. Am. 88, 284–290. Joswig, M., 1990. Pattern recognition for earthquake detection. Bull. Seism. Soc. Am. 80, 170–186. Kushnir, A.F., Lapshin, V.M., Pinsky, V.I., Fyen, J., 1990. Statistically optimal event detection using small array data. Bull. Seism. Soc. Am. 80, 1934–1947. Pinsky, V.I., Shapira, A., 1998. Kinematic multi-station discriminator between local quarry blasts and earthquakes. Geophys. J. Int. 135, 975–987. Ryzhikov, G.A., Birulina, M.S., Husebye, E.S., 1995. Automatic event location using local network data. Proceedings of the 17th Annual Seismic Research Symposium on Monitoring a Comprehensive test Ban Treaty, 12–15 September 1995, pp. 393–399. Shapira, A., Gitterman, Y., Pinsky, V., 1996. Discrimination of seismic sources using the Israel Seismic Network. Proceedings of the 18th Annual Seismic Research Symposium on Monitoring a Comprehensive test Ban Treaty, 4–6 September 1996, pp. 612–621. Shoubik, B.M., Kushnir, A.F., Haikin, L.M., Dainty, A., 1996. SCANLOC: automatic seismic event location based on local seismic network data. Proceedings of the 18th Annual Seismic Research Symposium on Monitoring a Comprehensive Test Ban Treaty, 4–6 September 1996, pp. 774–781.
An approach to automatic location of regional events V.I. Pinsky
)
Geophysical Institute of Israel, PO Box 2286, Holon 58122, Israel Received 18 May 1998; received in revised form 9 July 1998; accepted 9 July 1998
Abstract The conventional network location algorithms for automatic processing are based on P and S first arrivals. These procedures have the following shortcomings: Ž1. it is often difficult to identify the proper phase: ŽPg or Pn, Sg or Sn.; and Ž2. first arrivals are often masked by noise. Both factors may cause significant location errors. An alternative is to look for maximums of seismic energy time–distance distribution, which are less sensitive to these factors. We measure the time-of-maximum of P and S waves envelopes vs. distance for each available station of the Israel Seismic Network ŽISN., thus providing a travel time curve ŽTTC.. The record envelopes are obtained using 1D x 2 optimal detector and 3–6 Hz ‘short-time-average’ time-curves, having enhanced sensitivity for seismic signal arrivals. The corresponding P and S time-of-maximum vs. distance functions are approximated linearly by a least-square method for a set of local earthquakes and quarry blasts. Travel time inversion for location of small events comprises three main steps: Ž1. cleaning of the records from noise bursts and computation of the envelopes, Ž2. triggering and identification of the P and S phases and computation of their energy maximums, and Ž3. maximization of a sum of station residuals as a function of epicenter coordinates. The maximum of the functional is looked for on a 60 = 60 km grid, step of 2 km, covering different parts of Israel and Jordan. The algorithm is not sensitive to the source depth, thus providing epicenter determination only, but shows to be convenient and robust. As a result of the preliminary study for a set of relatively weak 74 local earthquakes and 58 quarry blasts, M L ; 1.5–2.5, we have obtained the accuracy of epicenter estimation "6 km for 80–90% of both types of events, which is satisfactory for automatic location. The accuracy is measured relative to the ISN bulletin locations for earthquakes and the ground truth information for quarries, respectively. q 1999 Elsevier Science B.V. All rights reserved. Keywords: Seismic event location; Envelopes; Travel time curves; Grid search optimization
1. Introduction Routine source location in seismology requires estimation from observed signal arrivals of four coordinates; geographical coordinates of the epicenter Žprojection of the point-source on the surface., depth
)
Fax: q972-3-5502925; E-mail: [email protected]
and source time. For regional events, an analyst marks onsets of Pg, Pn, Sg, Sn phases on the seismogram. These onset times are then converted to the corresponding source coordinates and time through the given travel time model, often determined by the velocity structure. For small and remote events with slow and gradual onsets, the accurate estimation of the first arrival time is problematic due to low signal-to-noise ratio. Thus, more reliable alternative is required. Sometimes, the depth can be considered known, e.g., when the source is classified as a quarry
0031-9201r99r$ - see front matter q 1999 Elsevier Science B.V. All rights reserved. PII: S 0 0 3 1 - 9 2 0 1 Ž 9 9 . 0 0 0 0 8 - 4
294
V.I. Pinskyr Physics of the Earth and Planetary Interiors 113 (1999) 293–302
blast or some other surface event. Depth knowledge is also not critical on a stage of preliminary location. In these cases, a possible solution is to utilize seismic energy distribution in time, space and frequency domains. The physical basis is that energy of the P and S seismic phases obey diffusion law and propagates with constant velocity in the stratified media. On this basis, Ryzhikov et al. Ž1995. derived theoretically a functional, which has a maximum, over the observed envelopes of the seismic phases, pointing to the source epicenter. Shoubik et al. Ž1996. applied principles of the so-called ‘seismic tomography’ based on the calculation of semblance and ‘signalto-noise ratio’ functionals over the envelopes of the seismograms. By the method, using a grid search approach, Shoubik et al. satisfactory estimated source epicenters of the 19 Galilee earthquakes, recorded by the Israel Seismic Network ŽISN.. Husebye et al. Ž1998. applied the short-time-average ŽSTA. procedure for calculation of the signal envelopes with the 3–6 Hz band-pass prefiltering the data of several seismic networks, and showed its effectiveness for epicenter determinations. In this paper, we use a similar approach focusing on the estimation of the seismic phases energy travel time curves ŽTTC.. We measure TTC for the maximums of the seismogram envelopes in the time intervals corresponding to P and S waves and use them for location instead of the usual P and S TTC of the first arrivals. For location, we use envelopes obtained through different types of filtering and develop a simple triggering procedure to extract proper P and S phases. The P and S phases vary in a complex manner and are often misidentified. In routine analysis, this may lead to large errors in source location. We concentrate on this problem and construct a procedure robust enough to phase wrong identifications. For the analysis, we use the expanded ISN database, including earthquakes and quarry blasts from the three different areas of Israel and its surroundings, and find their epicenters using a grid search method.
0.5–12 Hz. of the ISN Ž; 30 stations. telemetered to a Hub via radio FM link. The data include 132 seismic events which: Ž74 earthquakes, 58 quarry blasts. occurred in Galilee Ž30 earthquakes, 36 quarry blasts., Gilad Ž18 earthquakes. and the Dead Sea region Ž26 earthquakes and 22 quarry blasts. with magnitudes of M L s 1.0–2.8 at distances of 15–310 km. Event epicenters, locations of stations and quarries are presented in Fig. 1.
3. Method The whole procedure may be split in three blocks; Ži. preprocessor, Žii. picker and Žiii. locator. The Preprocessor performs despiking and a band-pass filtering of the seismograms in an appropriate frequency band. The records we are operating with are often spoiled by a different kind of spikes and noise bursts, caused by radio impulses, industry and micro-seismic noise that provide false detections. Thus, cleaning of records is crucial to the success of the whole procedure, based on detecting of seismic signals. To further enhance signal-to-noise ratio, we perform prefiltering Žwhitening. of each channel seismograms Vt by the procedure:
ž
Ut s Vt y Ý a k Vtyk k
/
s 2,
k s 1, . . . , M ,
Ž 1.
where a k and s 2 are the coefficients of the autoregression time series model, adapted to the noise preceding the signal. Ut is the whitened time-series being the output of the filter. Time interval for adaptation is selected equal to 200 samples Ž4 s. from the start of the recording. The next step is to calculate and square the sum of autocovariance coefficients Cm of the Ut in a moving window of length S, equivalent to a sort of envelope: 2
Zt s Ž C0,t y 1 . q Ý Cm2 ,t ,
m,s 1, . . . , M ,
m
Cm ,t s Ý Utqsym UtqsrS,
s s 1, . . . ,S.
Ž 2.
s
2. Data The data are presented by the short period seismograms Ž50 Hz sampling rate with efficient bandwidth
This procedure is equivalent to the statistically optimal detector by Kushnir et al. Ž1990., which enhances sensitivity of amplitude and whiteness changes of the Ut time-series due to the incoming signal. For the stationary noise, Vt variance of the
V.I. Pinskyr Physics of the Earth and Planetary Interiors 113 (1999) 293–302
295
Fig. 1. Map of Israel and adjacent areas, showing the three seismic database regions: 1—Galilee, 2—Gilad and 3—Dead Sea region.
corresponding white noise Ut is C0,t f 1, and autocorrelation Cm ,t f 0 for m ) 0, hence, Zt f 0. A
more explicit statement is that the product Zt S asymptotically obeys a x 2 probability distribution
296
V.I. Pinskyr Physics of the Earth and Planetary Interiors 113 (1999) 293–302
with M q 1 degrees of freedom for S ™ `. When the moving window reaches the signal in the seismogram, the time-series Vt is no longer stationary, and, hence, Ut is no longer white, which yields the strong increase of the sequence Zt values. Even for a small number of coefficients wwe choose M s 3, S s 50 Ž1 s window.x and rather weak seismic phases, Zt exhibits abrupt fronts and sharp energy maximums. In most cases, there is a clear separation between P and S phases even at short distances where these signals are close in time. The Eq. Ž1. noise-adapted filter may be efficient also for the suppression of Rg surface waves, excited by explosions and lying in the same frequency band as micro-seismic noise. For computational efficiency, the Zt curve is decimated by 1r10. Examples of the Zt curves are shown in Fig. 2; the Zt curves have local maximums Pmax , S max , at time points Tp and Ts , which are coincidental with P and S arrivals in the original non-transformed records. When M s 0, the function Zt coincides with the usual STA trace.
The Picker automatically provides Tp time tied to Pmax by finding a local maximum of Zt in the time window wTTp ,TTp q Ip x and similarly for Ts in the window wTTs ,TTs q Is x, where TTp and TTs are the two sequential moments of triggering and Ip and Is Ž Ip s Is s 5 s. are the specified time intervals, assumed to contain signal. The first triggering is pronounced if for the majority of points inside a moving window of length LS ŽLS s 50 samples or 10 s. Zt ) TL = LTA t Žthis requirement serves to skip spikes and bursts, left after the ‘preprocessing’., where TL s TL p ŽTL p s 10., LTA t is a long-timeaverage, continuously calculated as: LTA t s Ý Zkrt
k s 1, . . . ,t.
The second triggering is pronounced if the condition above is true for t ) TTp q Ip with a new trigger level TL s TL s ŽTL s s 100.. The Locator. In Fig. 3, we display the time–distance pairs ŽTp , R ., ŽTs , R . covering the distance range
Fig. 2. Optimal detector curves in logarithmic scale for the seismograms for four ISN seismic stations. Time picks Tp and Ts correspond to the local maximums Pma x and S max of the Zt curves Ženvelopes..
V.I. Pinskyr Physics of the Earth and Planetary Interiors 113 (1999) 293–302
297
Fig. 3. Tp and Ts time picks vs. distance for Ža. 30 earthquakes, Žb. 36 quarry blasts from Galilee. Best fitting lines—travel time curves—are also shown.
0–350 km for the 30 Galilee earthquakes and the 36 Galilee explosions. Fig. 3 shows that most of the
pairs are clustered along two straight lines. Therefore, travel times Tp Ž R ., TsŽ R ., for the arrivals of
V.I. Pinskyr Physics of the Earth and Planetary Interiors 113 (1999) 293–302
298
Pmax and S max , were estimated as TaŽ R . s T0a q RrVa , a s P, S. Using a simple least squares procedure Žexcluding outliers., we estimate for earthquakes Vp ( 6.4, Vs ( 3.6 and T0p ( 2.5 s, T0s ( 3.5 s and for explosions Vp ( 6.4, Vs ( 3.5, T0p ( 3.5 s, T0s ( 6. Thus, the slopes of the travel time functions for earthquakes and explosions are close, but the intercepts are different especially for S waves. We choose for our ‘locator’ Vp s 6.4, Vs s 3.5, T0p s 2.5, T0s s 3.5. The usual procedure for location, efficient in many cases, is a least-squares technique, where solution is due to minimization of a sum of squares of the residuals Da i s Ta y TaŽ R i . y T0 , for a number of stations, where T0 is the source time. If the observations include random outliers of the time Ta estimations, this procedure is no more efficient and besides, may cause large location errors. In Fig. 3, we see that there are many outliers due to both false triggering and false phase identification. Those which are below the P line are mostly spikes and bursts, while those over the S line are low velocity Rg surface wave triggers. Many P identifications by mistake give S travel time and vice versa. In such cases, residuals appear like Dac i s Ta y Tc Ž R i . y T0 , a,c s P,S. These features should be taken into account in the specialized location scheme, since the conventional least-squares procedure is not applicable. Next, instead of one station residuals Dac i , we shall consider two station residuals Dabc i j s Dac i y D bc j , i, j s 1, . . . , N, a,b,c s P,S, which do not depend on origin time T0 , thus, reducing the number of unknown variables from three to two. On the basis of the above discussion, we suggest to perform the location of the epicenters by maximization of the following target function: F Ž X 0 ,Y0 . s Ý
Ý Wi j A Ž Dabc i j . ,
ij abc
i s 1, . . . , N y 1, j s i q 1, . . . , N, a,b,c,s S,P 4 , Ž 3. comprising 8 = N = Ž N q 1.r2 observations. From the ‘picker’ above, we receive two time picks Ta and
T b for each channel, but the ‘picker’ is not aware to which phase: P, S, or Rg do they belong. To distinguish between correct and wrong identifications, we use the Eq. Ž3. function AŽ D . s Ž< D < q ´ .y1 of deviations D Žinstead of usual AŽ D . s D 2 ., ‘ voting’ for true pickings by its large values and ‘ voting’ against false pickings by its small values, respectively. Thus, the procedure unites in one the two tasks: the location with the phase identification, so that majority of ‘good’ pickings, close to P or S track, would guarantee successful solution. The control parameter ´ ) 0 helps to avoid possible zeroes of the residuals Dabc i j , thus providing robustness of the Xmax , Ymax estimator, which maximizes F Ž X,Y .. The parameter ´ is selected automatically as ´ s 0.1 q k = 0.1, k s 1, . . . ,30 according to empirical criterion that for more than Nr2 stations, there is at least one combination of a, b, c and j that Ž< Dabc i j < ´ . - 3. The Wi j coefficients are chosen and for the time being, equal 1, though too close and too remote stations, may be given lower weight. We prefer the simple case of two independent variables to obtain X max , Ymax through a simple grid search Žusing grid of 60 km = 60 km, 2 km step., though any other non-linear optimization procedure may be employed.
4. Preliminary location Our non-linear optimization procedure requires an initial value of X 0 , Y0 , preferably as close to the true source as possible. In this paper, we deal with local sources being, at most, some tens of kilometers from the nearest station. Hence, as a simple preliminary estimate, we choose the coordinates of the station which is triggered firstly. Let us call it station N1. Of course, due to the noise bursts preceding seismic signal the estimator may appear to be too crude. A likely improvement appears feasible by using observational experience that for relatively remote stations with R G 100 km or ŽTs y Tp G 10 s., Tp and Ts are determined more reliably than for closer stations. For
Fig. 4. Map and histogram of the 30 Galilee earthquake epicenters relative to ISN bulletin locations wŽ0,0. is true epicenter of each eventx, using Tp , Ts picks, provided by Ža. the optimal detector with 0.2–10 Hz prefiltering, Žb. STArLTA detector with 3–6 Hz prefiltering; and the Žc. best of the OD and STArLTA solutions according to the maximum of the target function Ž3.. Error ellipsoids, 80% and 90%, are provided by the normal distribution fit.
V.I. Pinskyr Physics of the Earth and Planetary Interiors 113 (1999) 293–302
299
V.I. Pinskyr Physics of the Earth and Planetary Interiors 113 (1999) 293–302
300
such remote stations, distance R can be directly found from the equation: Tp y Ts s T0 p y T0 s q R Ž 1rVp y 1rVs . ,
Ž 4.
and then the source origin time T0 is given by: T0 s Ý Ta j y Ž R jrVa . rN.
Ž 5.
j
Then for each station, observed phase time is: ta j s T0 q T0 a q R jrVa .
Ž 6.
Taking station N1 and some remote station N2, which fits Eq. Ž4. best, we solve two equations Ž6.
for N1 and N2 relative to X 0 , Y0 , using an equivalent quadratic equation. From Fig. 1, we see that Israel is clustered by tens of quarries which are close to the natural seismic sources. Taking into account that 90% of seismic events in the region are quarry blasts, it is reasonable to use as preliminary estimator coordinates of the quarry, which provide the largest value to the F Ž X,Y . function. Being applied together with any reliable automatic classification procedure, such as those described in ŽShapira et al., 1996. or ŽPinsky and Shapira, 1998., this decision may be regarded as
Fig. 5. Location of the Ža. 64 earthquakes and Žb. 58 explosions from the Galilee, Gilad and the Dead Sea regions, using the best of the OD and STArLTA solutions Žas for Fig. 4.. See text for discussion.
V.I. Pinskyr Physics of the Earth and Planetary Interiors 113 (1999) 293–302
final in the case when event is classified as a quarry blast, or at least competing with the full optimization procedure. It may be also regarded as methods of event recognition, such as described in ŽJoswig, 1990. and ŽFedorenko et al., 1998.. Finally, what estimator should be chosen as initial value for the optimization procedure? The answer depends on how many stations are triggered. In most cases, we choose the one which provides the largest contribution to the F Ž X,Y . function. For small events Žless than five stations triggering., we utilize the first station ŽN1. coordinates for epicenter determination.
5. Results Our epicenter location results in Figs. 4 and 5 are 1D by histograms of deviations R from X 0 , Y0 and 2D plots of points Ž X y X 0 , Y y Y0 .. Here, X 0 , Y0 are presumed true source coordinates as reported in bulletins in the case of an earthquake and ground truth coordinates of a quarry, in case of explosions. The 80% and 90% error ellipsoids are presented on the corresponding plots. Our picking procedure, described above, is parameterized in order to ensure operational flexibility. It proved to be robust with regards to variations of time window length and trigger level but sensitive to the number of coefficients M in the adaptive filter Ž2. and the frequency band-pass ŽFBP.. An instructive example is given in Fig. 4a,b with the usage of the event analysis procedures for 30 Galilee earthquakes: Ž1. M s 3 and FBP s w0.2–10 Hzx; and Ž2. M s 0 Žthe STA case., FBP s w3–6 Hzx. The final results look promising for both cases Ž95% within 8 km from epicenter and almost 100% within 10 km., though there are some outliers. For each event, one solution from the two cases was selected, corresponding to the larger value of function Ž3. for the fixed value of ´ . The result shown on Fig. 4c improves considerably Žmore than 70% of events are within 2 km range Ž1 grid step. and more than 95% within 4.5 km.. The above analysis procedures were applied to all the events in our dataset. Fig. 5a gives the results for the 74 earthquakes and Fig. 5b—for the 58 explosions. Statistically, the results of location for both types of events look very good and appear to be
301
almost the same: 84% of events within 4.5 km radius, 96% within 7.5 km, and almost 100% within 10 km. For many explosions, the locations were better than those obtained routinely, and this probably applies to the shallow earthquake locations as well.
6. Discussion and conclusions The above results look promising for automatization of local seismicity monitoring. Our novel procedure performed well providing robust and accurate solutions in the different geological regions of Israel, Galilee, Gilad and Dead Sea. Nevertheless, there is some scatter of locations, caused probably by difference in depth of source ranging between 0 and 30 km and by heterogeneity of the wave path. Some of the earthquakes that were located with error of 10–14 km were relatively deep and thus, probably, our assumption of depth independent velocity model of wave energy propagation was invalid. The most adverse factors influencing source location are the low signal-to-noise ratio and abrupt changes of the noise character preceding the phase arrival. The corresponding results are the ‘missing of target’ and the ‘false alarms’, respectively. Both provide the ŽTa , R . ‘outliers’, but really misleading are only those which deviate too much from the P and S curves in Fig. 3. The late P and S phase identifications are mostly due to the high noise level, while too early phase pickings are due to the sudden change in the noise frequency or amplitude. The Locator ŽEq. Ž3.. eliminates the outliers Žsee explanation above. successfully if their number is relatively small, while the rest of the P and S pickings are accurate. Otherwise, the outliers may ensure large but infrequent location errors as seen on Figs. 4 and 5. Nearly half of the events are explosions which occasionally proved problematic to analysis. The reason is that this kind of signal sources cause strong excitation of the Rg surface waves, propagating with low velocities in the range of 1–3 kmrs. This, in turn, may result in too late pickings of the true P and S phases, demonstrated in Fig. 3. The Rg waves are often characterized by signal frequencies, lower than for body waves and thus, filtering in the proper frequency band is the possible way for bettering the
302
V.I. Pinskyr Physics of the Earth and Planetary Interiors 113 (1999) 293–302
P- and S-pickings. Another factor is scattering and multi-pathing due to heterogeneity of the upper-most part of the crust. We chose the 3–6 Hz frequency band suggested by Husebye et al. Ž1998., and saw that filtering by a band-pass Butterworth filter with gradual slopes performed much better. For both earthquakes and explosions, efficient despiking was important. In spite of additional difficulties, final results of location for explosions are almost the same as for earthquakes. We think that some sort of match filters or velocity filters may be helpful for further refinements.
Acknowledgements I would like to thank Eystein Huseby and Avi Shapira for suggesting this research and for helpful advice and Anton Dainty for careful revision and valuable comments. I am grateful to Yuli Zaslavsky, Yefim Gitterman and Ilya Turetsky for fruitful discussions. Many thanks to Boris Gurevich for careful reading of the text and useful comments and to Nitzan Rabinovich for the text corrections. This study was supported partially by the U.S. Department of Defense, Contract No. DSWA01-97-C-0151 and partially by the German Israel Foundation for Scientific Research and Development, Contract No. G 0463245.08r95.
References Fedorenko, Yu.V., Husebye, E.S., Heincke, B., Ruud, B.O., 1998. Recognizing explosion sites without seismogram readings: neural network analysis of envelope-transformed multi-station SP recordings 3–6 Hz. Geophys. J. Int. 133, F1–F6. Husebye, E.S., Ruude, B.O., Dainty, A.M., 1998. Fast, robust and reliable epicenter determinations: envelope processing of local network data. Bull. Seism. Soc. Am. 88, 284–290. Joswig, M., 1990. Pattern recognition for earthquake detection. Bull. Seism. Soc. Am. 80, 170–186. Kushnir, A.F., Lapshin, V.M., Pinsky, V.I., Fyen, J., 1990. Statistically optimal event detection using small array data. Bull. Seism. Soc. Am. 80, 1934–1947. Pinsky, V.I., Shapira, A., 1998. Kinematic multi-station discriminator between local quarry blasts and earthquakes. Geophys. J. Int. 135, 975–987. Ryzhikov, G.A., Birulina, M.S., Husebye, E.S., 1995. Automatic event location using local network data. Proceedings of the 17th Annual Seismic Research Symposium on Monitoring a Comprehensive test Ban Treaty, 12–15 September 1995, pp. 393–399. Shapira, A., Gitterman, Y., Pinsky, V., 1996. Discrimination of seismic sources using the Israel Seismic Network. Proceedings of the 18th Annual Seismic Research Symposium on Monitoring a Comprehensive test Ban Treaty, 4–6 September 1996, pp. 612–621. Shoubik, B.M., Kushnir, A.F., Haikin, L.M., Dainty, A., 1996. SCANLOC: automatic seismic event location based on local seismic network data. Proceedings of the 18th Annual Seismic Research Symposium on Monitoring a Comprehensive Test Ban Treaty, 4–6 September 1996, pp. 774–781.