A Bayesian approach to the realtime 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 realtime determination of source parameters. Inthis paper we describe a probabilistic evolutionary approach for the realtime magnitudeestimation which can have a potential use in earthquake early warning. The techniqueis based on empirical prediction laws correlating the lowfrequency 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 2s 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 2s 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 nearsource
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 realtime 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 realtime the probability of lossesand to activate automatic and/or manual procedure tomitigate earthquake effects. An EEWS is a modern realtime 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 (sitespecific) [
Kanamori
, 2005;
Nakamura
, 1988].[
3
] A sitespecific 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.01480227/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 estimation of source parameters requires that earthquake locationand magnitude, and their associated uncertainty, are determined 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 generalized 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 information 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 magnitude 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
ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 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 lowfrequency 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,nearsource 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 realtime magnitude algorithm based on empirical prediction laws correlating the lowfrequency 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 EuroMediterranean earthquake records andthen confirmed on Japanese earthquakes recorded by theKNet and KiKNet 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 GutenbergRichter 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 magnitude 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 KNet andKiKNet 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: REALTIME 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 KNet seismicnetwork. We integrated this waveform database for eventswith
M
> 6.1 using 68 additional records from the KiKNet 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 EuroMediterranean (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 amplitudes than for EuroMediterranean 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
¼
ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
NS
2
þ
EW
2
þ
UD
2
p
ð
2
Þ
where NS, EW, and UD are the northsouth, east–west, andvertical components, respectively.[
23
] A preliminary
P
arrival identification has been performed 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: REALTIME 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 achisquare 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 2swindow, 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 2s windows after the first
S
arrival, while no saturation effect is visible at
M
> 6.5 as for the
P
data in the 2s time window.
3. RealTime Bayesian Estimation of EarthquakeMagnitude
[
29
] On the basis of the observed correlation relationships, we implemented a probabilistic, evolutionary algorithm for realtime 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 distancenormalized 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 information 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 lowpassfiltered 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: REALTIME 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
ﬃﬃﬃﬃﬃﬃ
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 information 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 information 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 measurements:
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 estimated 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 1D 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
picking 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 displacements regression (8), the mean values and standard deviation 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 2s
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 2s
P
dataduring the realtime magnitude estimation, so that theregression laws for both 2 and 4s
P
data (in addition to2s
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 2s
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 2s
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 nearsource 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 associated 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: REALTIME MAGNITUDE5 of 17
B12302