A Bayesian approach to the real-time estimation of magnitude from the early P and S wave displacement peaks

of 17

Please download to get full document.

View again

All materials on our website are shared by users. If you have any questions about copyright issues, please report us to resolve them. We are always happy to assist you.
17 pages
0 downs
A Bayesian approach to the real-time estimation of magnitude from the early P and S wave displacement peaks
  A Bayesian approach to the real-time estimation of magnitude fromthe early  P   and  S   wave displacement peaks M. Lancieri 1 and A. Zollo 2 Received 18 September 2007; revised 2 August 2008; accepted 7 August 2008; published 3 December 2008. [ 1 ]  It has been shown that the initial portion of   P   and  S   wave signals can provideinformation about the final earthquake magnitude in a wide magnitude range. Thisobservation opens the perspective for the real-time determination of source parameters. Inthis paper we describe a probabilistic evolutionary approach for the real-time magnitudeestimation which can have a potential use in earthquake early warning. The techniqueis based on empirical prediction laws correlating the low-frequency peak grounddisplacement measured in a few seconds after the  P   and/or   S   phase arrival and the finalevent magnitude. The evidence for such a correlation has been found through the analysisof 256 shallow crustal events in the magnitude range  M   jma   4–7.1 located over the entireJapanese archipelago. The peak displacement measured in a 2-s window from thefirst   P   phase arrival correlates with magnitude in the range  M   = [4–6.5]. While a possiblesaturation effect above  M  ’ 6.5 is observed, it is less evident in an enlarged window of 4 s.The scaling of   S   peaks with magnitude is instead also observed at smaller time lapses(i.e., 1 s) after the first   S   arrival. The different scaling of   P   and  S   peaks with magnitudewhen measured in a 2-s window is explained in terms of different imaged rupture surface by the early portion of the body wave signals. We developed a technique to estimatethe probability density function (PDF) of magnitude, at each time step after the event srcin. The predicted magnitude value corresponds to the maximum of PDF, while itsuncertainty is given by the 95% confidence bound. The method has been applied to the2007 (  M   jma   = 6.9) Noto Hanto and 1995 (  M   jma   = 7.3) Kobe earthquakes. The results of this study can be summarized as follows: (1) The probabilistic algorithm founded onthe predictive model of peak displacement versus final magnitude is able to provide a fast and robust estimation of the final magnitude. (2) The information available after a fewseconds from the first detection of the  P   phase at the network can be used to predict the peak ground motion at a given regional target with uncertainties which are comparable tothose derived from the attenuation law. (3) The near-source  S   phase data can be used jointly with  P   data for regional early warning purposes, thus increasing the accuracy andreliability of magnitude estimation. Citation:  Lancieri, M., and A. Zollo (2008), A Bayesian approach to the real-time estimation of magnitude from the early  P   and  S  wave displacement peaks,  J. Geophys. Res. ,  113 , B12302, doi:10.1029/2007JB005386. 1. Introduction [ 2 ] The final aim of an earthquake early warning system(EEWS) is to evaluate in real-time the probability of lossesand to activate automatic and/or manual procedure tomitigate earthquake effects. An EEWS is a modern real-time information system able to provide rapid notification of an occurring earthquake through the fast telemetry and processing of data from dense instrument arrays deployedin the source region of the event (regional) or surroundingthe target infrastructure (site-specific) [  Kanamori , 2005;  Nakamura , 1988].[ 3 ] A site-specific EEWS consists of a single sensor or anarray deployed in the proximity of the target site, and whosemeasurements of amplitude and predominant period on theinitial  P   wave motion [ Wu and Kanamori , 2005a, 2005b]are used to predict the ensuing peak ground motion (mainlyrelated to the arrival of late  S   and surface waves) at the samesite. A regional EEWS is based on a dense sensor network covering a part of or the entire seismic source area. Therelevant source parameters (event location and magnitude)are estimated from the early portion of the recorded signalsand are used to predict, with a quantified confidence, aground motion intensity measure at a distant location wherea target site of interest is located. The success of a regionalEEWS critically depends on the capability to provide in a JOURNAL OF GEOPHYSICAL RESEARCH, VOL. 113, B12302, doi:10.1029/2007JB005386, 2008 Click Here for Full A rticle 1 Istituto Nazionale di Geofisica e Vulcanologia, Sezione di Napoli, Naples, Italy. 2 Dipartimento di Scienze Fisiche, Universita` di Napoli Federico II, Naples, Italy.Copyright 2008 by the American Geophysical Union.0148-0227/08/2007JB005386$09.00 B12302  1 of 17  fast, secure and reliable way the size and location of anearthquake while it is still occurring [  Kanamori , 2005].[ 4 ] Standard techniques used to determine the source parameters of an earthquake are generally not suited for early warning applications, since they typically use thecomplete event waveforms from all the stations of a localseismic network, with a consequent large time delay for datatransmission and analysis. For this reason, a different strategy is required, where the computation starts when afew seconds of data from a small number of recordingstations is available, and the source parameters estimates are progressively updated with time.[ 5 ] Assuming that a dense seismic network is deployedaround the earthquake source zone, an evolutionary estima-tion of source parameters requires that earthquake locationand magnitude, and their associated uncertainty, are deter-mined as a function of time, i.e., of the number of recordingstations and of the length of ground motion signal recordedat each station.[ 6 ] Recently, robust and effective algorithms for the realtime and evolutionary event location have been developed[  Horiuchi et al. , 2005;  Satriano et al. , 2008;  Cua and  Heaton , 2007]. In particular,  Cua and Heaton  [2007],generalized the concept of not yet arrived data, described by  Horiuchi et al.  [2005] in order to start the location withone single triggering station.[ 7 ]  Satriano et al.  [2008] further extended and general-ized the  Horiuchi et al.  [2005] approach by developing afully probabilistic technique where the hypocenter isdescribed as a PDF instead of as a point. Their techniquemakes use of the equal differential time formulation [  Font et al. , 2004] to incorporate both triggered arrivals and not yet triggered stations and applies a nonlinearized, global searchfor each update of the location process.[ 8 ] On the other hand the feasibility of determining inreal time the final earthquake magnitude using the informa-tion contained in the initial portion of the recorded  P   and  S  signals is matter of investigation and debate [ Olson and  Allen , 2005;  Rydelek and Horiuchi , 2006;  Rydelek et al. ,2007;  Hellemans , 2006;  Zollo et al. , 2006, 2007].[ 9 ] In the ElarmS methodology [  Allen , 2007] the magni-tude is rapidly estimated using the frequency content in awindow of 4 s of signal after the first   P   wave arrival.Founded on a method described first by  Nakamura  [1988],the predominant period,  t   p  of the vertical component waveform, is based upon a recursive calculation of squaredground velocities and displacement  1 t   p ¼  12 p   ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R  t  0 0  _ u 2 t  ð Þ dt  R  t  0 0  u 2 t  ð Þ dt  v uut  ð 1 Þ where  u  is the vertical component of displacement; and  t  0  isthe time window duration (i.e., 4 s). The maximum  t   p  valuewithin 4 s is found to scale with event magnitude [  Allen and  Kanamori , 2003;  Olson and Allen , 2005;  Wu and  Kanamori , 2005a].[ 10 ] An alternative to the predominant period is theanalysis of low-frequency peak displacement measured ina few seconds after the  P   phase arrival, PD.  Wu et al.  [2006]investigated the attenuation relationship between PD, thehypocentral distance and the magnitude (estimated as localmagnitude) on southern California [ Wu et al. , 2006;  Wu and  Zhao , 2006] and Taiwan events.[ 11 ] On the basis of the analysis of the Mediterranean,near-source strong motion database,  Zollo et al.  [2006]showed that the peak displacement amplitude of initial  P  and  S   scales with magnitude in the moment magnituderange 4    M  w     7. This relationship may be of potential usefor EEWS application since the earthquake size can be,therefore, estimated using only a couple of seconds of signalfrom the  P   or   S   wave onset, i.e., while the rupture itself isstill propagating and before the final rupture dimension isachieved. The potential use of   S   waves for early warning purposes is one of main implications of the results of   Zolloet al.  [2006, 2007], despite that their use for early warning purposes is usually excluded [e.g.,  Rydelek et al. , 2007].[ 12 ] In this work we develop a Bayesian real-time mag-nitude algorithm based on empirical prediction laws corre-lating the low-frequency peak displacement measured onthe initial portion of   P   and  S   wave of strong motion signalswith the final earthquake magnitude.[ 13 ] This correlation has been pointed out by  Zollo et al. [2006] for the Euro-Mediterranean earthquake records andthen confirmed on Japanese earthquakes recorded by theK-Net and KiK-Net strong motion network [  Zollo et al. ,2007].[ 14 ] At any time  t  , after the first event detection theconditional probability density function (PDF) of magnitudegiven the observed data vector is expressed, through theBayes’ theorem, as the product between the conditionalPDF of data, given the magnitude, and a priori PDF basedon the prior information available before the time  t  . In the proposed methodology the a priori distribution is derived bythe Gutenberg-Richter occurrence law, similarly to  Iervolinoet al.  [2006] and  Cua and Heaton  [2007].[ 15 ] In order to investigate the feasibility of early warningin Japan using the proposed approach, we have applied it totwo large events that occurred recently in Japan (Kobe 1995and Noto Hanto 2007 [  ABS Consulting  , 2007]) and whoserecords have not been used to determine the prediction law. 2. Peak Ground Displacement and EarthquakeMagnitude 2.1. Empirical Correlation Law for Japan [ 16 ] With the aim to determine an empirical correlationlaw between early peak displacement and magnitude for shallow crustal Japanese earthquakes, we analyzed thestrong motion records from 256 Japanese events, with amaximum depth of 50 km and magnitude ranging from 4to 7.1 distributed over the entire Japanese archipelago(Figure 1).[ 17 ] The data are preliminary grouped into magnitude bins of width 0.3, which approximates the standard error onmagnitude determination [  Katsumata , 1996]. The magni-tude bin grouping has the additional advantage to averageout possible effects on peak amplitudes due to the fault mechanism and rupture directivity.[ 18 ] For earthquake waveforms extracted from K-Net andKiK-Net catalogues, the information concerning the event locations and magnitudes (  M   jma  ) are provided by the JapanMeteorological Agency (JMA). For earthquakes shallower  B12302  LANCIERI AND ZOLLO: REAL-TIME MAGNITUDE2 of 17 B12302  than 60 km  M   jma   is determined by  Tsuboi ’s [1954]relationship.[ 19 ] By a comparison between the  M   jma   and moment magnitude  M  w   retrieved from the CMT solution [  Dziewoskiet al. , 1981],  Katsumata  [1996] found negligible differences between the two units for shallow earthquakes in themagnitude range 5 to 7. Hereinafter we simply indicatemagnitude as  M  .[ 20 ] A total of 2640 records, with hypocentral distancesHD < 60 km have been extracted from the K-Net seismicnetwork. We integrated this waveform database for eventswith  M   > 6.1 using 68 additional records from the KiK-Net  borehole strong motion network.[ 21 ] The distance distribution of analyzed records for Japan earthquakes [  Zollo et al. , 2007, Figure 1] is rather different from that for the Euro-Mediterranean (MED)earthquakes in the work by  Zollo et al.  [2006, Figure 1].Whereas most of MED records occur at an HD < 25 km,most Japanese records are at HD > 40 km. Despite a muchlarger number of records for Japan, we would, therefore,expect a stronger distance attenuation effect on peak ampli-tudes than for Euro-Mediterranean data. 2.2. Data Processing [ 22 ] The procedure of data analysis involved sensitivitycorrection,  P   and  S   phase identification, double integrationandfilteringinafrequencyband[0.075–3]Hz.  P  and S   peaksare read on the modulus of displacement, expressed inmeters, defined as  D  ¼  ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   NS 2 þ  EW 2 þ  UD 2 p   ð 2 Þ where NS, EW, and UD are the north-south, east–west, andvertical components, respectively.[ 23 ] A preliminary  P   arrival identification has been per-formed through an automatic first-  P   picker [  Allen , 1978], Figure 1.  Geographic distribution of the selected event on the Japanese archipelago. The grey linesrefer to plate boundaries. The Pacific plate is subducting from the east beneath the North American andEurasian plates in eastern Japan; the Philippine Sea plate is descending from the south beneath theEurasian plate in southwest Japan. Large intraplate earthquakes occur along the plate boundaries off the Pacific coast of the Japan islands. Intraplate earthquakes within the continental plate take place in theupper crust beneath the Japan islands and along the coast of the Sea of Japan. Although the crustalintraplate earthquakes do not occur as frequently as the intraplate earthquakes, they generally inflict greater damage because they are shallow and near densely populated areas [  Zhao et al. , 2000]. B12302  LANCIERI AND ZOLLO: REAL-TIME MAGNITUDE3 of 17 B12302  then each pick has been manually reviewed. The first   S  arrivals from the horizontal components of all the selectedstrong motion records have been identified and pickedmanually.[ 24 ] Starting from the  P   wave and  S   wave picked arrivals,we measure the maximum amplitude on increasing timewindows, with duration of 2–4 s for the  P   phase and 1–2 sfor the  S   phase.[ 25 ] In order to correct the  P   and  S   wave peak amplitudefor the distance attenuation effect, we assume an attenuationrelationship of the form [ Campbell  , 1985;  Wu and Zhao ,2006;  Zollo et al. , 2006]: log PD ð Þ ¼  A  þ  BM   þ  C   log  R ð Þ ð 3 Þ where PD is the observed peak displacement amplitude,  M  is the magnitude, and  R  is hypocentral distance. The parameters  A ,  B , and  C   have been estimated by a standard best fit regression analysis and are reported along with their uncertainties in Table 1. Using the attenuation relationship(3) retrieved for shallow Japanese earthquakes, the  P   and  S   peaks have been normalized to a reference distance of 10 km[e.g.,  Zollo et al. , 2007]. Figure 2 shows the relation between the logarithm of distance normalized, peak displacement for   P   and  S   waves in the different timewindows and the final magnitude of the analyzed events.[ 26 ] The hypothesis of a Gaussian distribution of peak readings in each magnitude bin has been verified using achi-square test with confidence level of 95%. We can,therefore, consider that the mean and the standard deviationof peak values are adequate estimators of the statisticaldistribution.[ 27 ] Looking in detail at   P   data distribution for the 2-swindow, one can note that peak amplitudes are log linearlycorrelated with magnitude up to about   M   = 6.5, while at larger magnitudes the peaks scatter around a nearly flat trend indicating a saturation effect. We point out the rather  poor   P   data sampling for magnitude ranging between  M   =6.5 and  M   = 7.1 (only 8 events with  M   > 6.5).[ 28 ] Using a larger time window (4 s) the  P   peak displacement appears to correlate with magnitude all over the whole explored magnitude range. Interestingly, the  S   peaks show an excellent correlation with magnitude in thefull considered range for 1- and 2-s windows after the first   S  arrival, while no saturation effect is visible at   M   > 6.5 as for the  P   data in the 2-s time window. 3. Real-Time Bayesian Estimation of EarthquakeMagnitude [ 29 ] On the basis of the observed correlation relation-ships, we implemented a probabilistic, evolutionary algo-rithm for real-time magnitude evaluation (RTMag). 3.1. A Bayesian, Evolutionary Approach [ 30 ] At any time step  t   from the first event detection, onecan define the conditional probability density function(PDF) of magnitude, given the observed data vector   d  t  (i.e., the distance-normalized peak displacement read inthe early portion of   P   and  S   phase):  f m j d  t  ð Þ ¼  f d  t  j m ð Þ  f   t   m ð Þ R   M  max  M  min  f d  t  j m ð Þ  f   t   m ð Þ dm ð 4 Þ where  f   ( m j d  t  ) is the conditional PDF of   M   given data  d  t  ,  f   ( d  t  j m ) is the conditional PDF of   P   and/or   S   data  d  t   givenmagnitude  M  , and  f   t  ( m ) is the ‘‘prior’’ probability densityfunction of magnitude  M  , which expresses the probability of earthquake magnitude based on the available prior informa-tion before the current time  t  . Table 1.  Coefficients of Estimated Curves Along With RetrievedStandard Errors Phase Seconds  A c 0 D  A c 0  B c 0 D  B c 0 SE c  C  c  D C  c  P   2   6.93 0.48 0.75 0.09 0.32   1.13 0.06  P   4   6.46 0.48 0.70 0.09 0.40   1.05 0.10 S   1   6.03 0.54 0.71 0.10 0.38   1.40 0.05 S   2   6.34 0.48 0.81 0.09 0.37   1.33 0.05 Figure 2.  Correlation between low-pass-filtered peak ground motion value and magnitude. Thelogarithm of peak ground displacement normalized at a reference distance of 10 km as a function of   M   jma  is shown. Black diamonds represent the peak average on each magnitude bin with the associated standarddeviation, while the grey dots are the peak values read on each record. The black solid lines represent the best fit line evaluated through a linear regression. The dashed lines show the prediction bounds, i.e., the95%confidenceintervalforanewobservationwhichwillbecenteredon  b Y   =log(  PD 10 km ) k  whoselengthisgiven by  V  [  b Y  ] =  s {(1/  n ) + [(  M  k       M  ) 2 /(  M  k       M  ) 2 ]} 2 , where  s  is the root of the mean square error,  n  is thenumber of observation, and    M   is the mean value of magnitude. B12302  LANCIERI AND ZOLLO: REAL-TIME MAGNITUDE4 of 17 B12302  [ 31 ] The time step is defined as the minimum windowlength of transmitted data packets, here assumed as 1 s.[ 32 ] Assuming that at a given time  t  , the observed data  d  t  ,from an event of given magnitude, have a lognormaldistribution and are statistically independent, the PDF  f   ( d  t  j m ) is expressed by the likelihood  f d  t  j m ð Þ¼ Y  N k  ¼ 1 1  ffiffiffiffiffiffi 2 p  p   s  log  d  t  ; k  ð Þ d  t  ; k   exp   log  d  t  ; k     m  m ð Þ log  d  t  ; k  ð Þ   2 2 s  2log  d  t  ; k  ð Þ 264375 ð 5 Þ where  N   is the number of stations whose data are availableat time  t  .[ 33 ] Following  Iervolino et al.  [2006], the prior informa-tion on magnitude is given by the Gutenberg and Richter relationship for times before the first   P   arrival ( t    first  ) isdetected at the seismic network. As new data becomeavailable from the seismic network ( t   >  t    first  ), the prior PDF must change according to the new available informa-tion on earthquake magnitude provided by the recordedwaveform:  f   t  > t   first   m ð Þ¼  f   t   1  m j d  t   1 ð Þ ð 6 Þ 3.2. Magnitude Estimation From Peak Measurements [ 34 ] At each time step after the first   P   detection, theseismic network outputs a set of peak amplitude measure-ments:  P   peak in windows of 2 and 4 s after first   P  ;  S   peak in windows of 1 or 2 s after first   S  .[ 35 ] While the  P   wave onset is identified by an automatic picking procedure [  Allen , 1978], the  S   onset can be esti-mated from an automatic  S   picking or from a theoretical prediction based on the hypocentral distance given by theactual earthquake location.[ 36 ] To verify the quality of the estimated  S   arrivals, onthe basis of the hypocentral location, we extracted fromeach magnitude bin a set of events, with at least six records,and we located them using the 1-D velocity model by  Fukuyama et al.  [2003]. The obtained hypocentral estimatesare then used to calculate the  S   arrival times which arecompared to the manual picking. The distribution of thedifferences between the theoretical and the manual  S   pick-ing is centered on zero and have a standard deviation of ±0.5 s.[ 37 ] At a given time step after the first   P   wave detectionat the network, progressively refined estimates of magnitudeare obtained from  P   and  S   peak displacement data. Theseare preliminarily normalized to a reference distance  R  =10 km: log PD 10km c   ¼ log PD  Rc    C  c  log  R 10    ð 7 Þ where the subscript   c  is for   P   and  S   phase, constant   C  c  isdetermined through the best fit regression expressed inequation (3), and  R  is the hypocentral distance, expressed inkilometers. Then the distance corrected peaks are correlatedwith the final magnitude through the relationship log PD 10km c   ¼  A 0 c þ  B 0 c  M   ð 8 Þ [ 38 ] Assuming a standard error of SE c  on peak displace-ments regression (8), the mean values and standard devia-tion of the quantity log (PD c ), can be written as m  m ð Þ log PD c ð Þ ¼  B 0 c m þ  A 0 c þ C  c  log  R 10   s  log PD c ð Þ  ¼ SE c þ log  R 10   D C  c þ C  c 1  R D  R 8>>>>><>>>>>: ð 9 Þ where D  R  and D C  c  are the errors on the distance  R  and onthe  C  c  coefficient in equation (7). In this framework it is possible to evaluate the PDF using several strategiesincluding the use of 2-s  P   peak regression laws jointlywith the  S   information.[ 39 ] As discussed in section 3.3, the saturation effect doesnot prevent the use of the regression law for 2-s  P   dataduring the real-time magnitude estimation, so that theregression laws for both 2- and 4-s  P   data (in addition to2-s  S   wave data) can be usefully included in a Bayesian probabilistic approach.[ 40 ] By defining  F  t  ,  M  ( m ) as the cumulative PDF at time  t  ,it is possible to estimate a range of magnitude variation[  M  min ;  M  max ], whose limits are based on the shape of the  F  t  ,  M  (  M  ) function:  M  min  :  F  t  ;  M  min  m ð Þ¼ R   M  min 1  P  t   m j d  t  ð Þ dm ¼ a  M  max  :  F  t  ;  M  max  m ð Þ¼ R   M  max 1  P  t   m j d  t  ð Þ dm ¼ 1  a ð 10 Þ For example, if we assume  a  = 5%, then  M  min  and  M  max will be the  F  t  ( m ) evaluated at 0.05 and 0.95, respectively.Alternatively, the PDF can be used to predict the probabilityof exceeding a critical magnitude threshold (  M  cth ) given theestimated magnitude (  M  est  ) defined as  P M   >  M  cth j  M  est  ð Þ¼ Z   þ1  M  cth  P  t   m j d  t  ð Þ dm  ð 11 Þ 3.3. A Statistical Model to Account for the 2-s  P   Peak Saturation Effect [ 41 ] In this section we discuss how to include in our Bayesian method the empirical correlation law for   P   peak displacement measured in a window of 2 s. The use of 2-s  P  data (PD  P  2 s ) for magnitude estimation is important for several reasons: (1) the PD  P  2 s regression law can predict the magnitude of events below the saturation threshold(  M  sat  ), (2) PD  P  2s data are the only  P   wave information fromthe near-source station with  S  -  P   time window shorter than4 s, and (3) although more uncertain, the PD  P  2s data can provide a preliminary magnitude evaluation, and its associ-ated PDF can be used as prior information for further magnitude estimates.[ 42 ] Let us hypothesize that the distribution for PD  P  2s datafollows the lognormal distribution for   M   <  M  sat   while it is B12302  LANCIERI AND ZOLLO: REAL-TIME MAGNITUDE5 of 17 B12302
Related Search
We Need Your Support
Thank you for visiting our website and your interest in our free products and services. We are nonprofit website to share and download documents. To the running of this website, we need your help to support us.

Thanks to everyone for your continued support.

No, Thanks