M. Yamada, T. Heaton, and J. Beck

1

Real-time Estimation of Fault Rupture Extent Using
Near-source versus Far-source Classification
Masumi Yamada1, Thomas Heaton2, and James Beck2
masumi@eqh.dpri.kyoto-u.ac.jp, heaton t@caltech.edu, jimbeck@caltech.edu

Abstract: To estimate the fault dimension of an earthquake in real time, we present
a methodology to classify seismic records into near-source or far-source records. Characteristics of ground motion, such as peak ground acceleration, have a strong correlation withthe distance from a fault rupture for large earthquakes. This study analyzes peak groundmotions and finds the function that best classifies near-source and far-source records basedon these parameters. We perform: (1) Fisher’s linear discriminant analysis and two different Bayesian methods to find the coefficients of the linear discriminant function; and (2)Bayesian model class selection to find the best combination of the peak ground motion parameters. Bayesian model class selection shows that the combination of vertical accelerationand horizontal velocity produces the best performance for the classification. The linear discriminant function produced by the three methods classifies near-source and far-source data,
and, in addition, the Bayesian methods give the probability for a station to be near-source,
based on the ground motion measurements. This discriminant function is useful to estimate
the fault rupture dimension in real time, especially for large earthquakes.

1
2

Kyoto University, Gokasyo, Uji, 611-0011 Japan
California Institute of Technology, MC104-44, 1200 E. California Blvd. Pasadena CA 91125 USA

M. Yamada, T. Heaton, and J. Beck

1

2

Introduction

Recent studies show that earthquake early warning systems, such as the Virtual Seismologist(VS) Method (Cua, 2005; Cua and Heaton, 2006), can accurately estimate the location of
the epicenter a few seconds after the first arrival station records the ground motion of themain shock (Nakamura, 1988; Allen and Kanamori, 2003; Odaka et al., 2003; Wu andKanamori, 2005). The VS method assumes a point source model for the rupture, and it
works well for small to moderate earthquakes (magnitude < 6.5) (Cua, 2005). However,for large earthquakes, the fault rupture length can be on the order of tens to hundreds ofkilometers, and the prediction of ground motion at a site requires approximated knowledgeof the rupture geometry. Early warning information based on a point source model may
underestimate the ground motion at a site, if a station is close to the fault and distant fromthe epicenter. This occurs because, for large earthquakes, the peak characteristics of groundmotion, such as peak ground acceleration, have stronger correlation with the fault rupture
distance rather than with the epicentral or hypocentral distance (Campbell, 1981). (Thedefinition of the fault rupture distance in this paper is the shortest distance between thestation and the surface projection of the fault rupture surface.)
In order to construct an early warning system that is more reliable for large earthquakes,
it is necessary to estimate the fault rupture extent and slip on the fault in real time. The
objective of this paper is to develop a methodology to classify stations into near-source andfar-source since this can be used for identifying the fault geometry if there is a sufficiently
dense seismic network. Peak ground motions recorded in past earthquakes are analyzed to
predict whether a station recording ground motion is close to the earthquake fault area. This
classification problem can be stated as follows: given ground motion data from past earthquake records, what is the probability that a station is near-source when a new observation
is obtained?
To approach this problem, we:1) Collect strong motion data from earthquake strong motion archives and classify thesesamples into two predefined groups: records from near-source stations and far-source stations.
This particular set of data is called the training set.
2) Discover a discriminant function of the samples features (e.g. peak ground acceleration(PGA), velocity (PGV), displacement (PGD)) which provides the best performance in termsof near-source / far-source classification.
3) Allocate new observations when they are obtained to one of the two groups based on
the discriminant function.
The first step is quite straightforward; strong motion data from past earthquakes arecollected based on certain selection criteria. The second step is the main topic of this paper;and we investigate linear discriminant functions by using the traditional Fisher method andtwo Bayesian methods. The third step can then be accomplished in a real-time analysis.
Given a new ground motion observation from on-going rupture, the discriminant function
gives the probability that the observation is located in the near-source.

M. Yamada, T. Heaton, and J. Beck

2

3

Strong Motion Data

We used strong motion datasets from nine earthquakes with magnitude greater than 6.0
and containing records of near-source stations. The selected earthquake dataset is shownin Table 1. Here, we define a near-source station as a station whose fault rupture distanceis less than 10km. 695 three-component strong motion data are used for the classification
analysis and 14% (100 stations) are from near-source stations.

2.1

Data Sources

We obtained the strong motion dataset for the Imperial Valley (October 15, 1979), Loma Prieta (October 18, 1989), Landers (June 28, 1992), Northridge (January 17, 1994), and Denali(November 3, 2002) earthquakes from the COSMOS Virtual Data Center (http://db.cosmoseq.org) which includes data from the California Strong Motion Instrumentation Program(CSMIP) seismic network and the United States Geological Survey (USGS) seismic network. The Northridge earthquake dataset in the COSMOS Virtual Data Center also includes records from seismic networks of the California Institute of Technology, Los Angeles Department of Water and Power, Metropolitan Water District, Southern California
Earthquake Center, and University of Southern California. All these data were recordedby accelerometers and processed appropriately before distribution to users. The correction process may apply baseline corrections, bandpass filters to remove noise contamination,and instrument correction to remove the effects of frequency-dependent instrument response(http://nsmp.wr.usgs.gov/processing.html).
Strong motion data from the Hyogoken-nanbu earthquake (January 16, 1995) are provided by Japan Meteorological Agency (JMA), the Committee of Earthquake Observation
and Research in the Kansai Area (CEORKA) in Japan (Toki et al., 1995), and the Japan
Railway Institute (JR) whose records were scanned and digitized by Wald (1996). Seismometers installed in the CEORKA network record velocity, and those records are differentiatedonce to obtain accelerograms.
The national strong-motion accelerograph network in Turkey recorded the strong motions during the Izmit earthquake (August 17, 1999) (Akkar and Gülkan, 2002). They
can be downloaded from the ftp site of the Earthquake Research Department of General
Directorate of Disaster Affairs, Ministry of Public Works and Settlement, Ankara, Turkey(ftp://angora.deprem.gov.tr/). The COSMOS Virtual Data Center archived the datasetof another network operated by Kandilli Observatory and Earthquake Research Institute,
Earthquake Engineering Department, Bogaziçi University, Istanbul, Turkey. Stations withfault distance greater than 200 km are excluded since ground motion amplitudes of thosestations are quite small which results in a low signal-to-noise ratio. We use four digitaland six analog acceleration records from the national network and eight digital acceleration
records from the Bogaziçi University network.
The Chi-Chi earthquake (September 20, 1999) is one of the best recorded earthquakes
with a large number of stations and a dense station distribution both in the near-source andfar-source. Strong motion records for the Chi-Chi earthquake are available on the attached
CD in the Special Issue of the Bulletin of the Seismological Society of America, Vol. 93,

M. Yamada, T. Heaton, and J. Beck4

No. 5 (Lee et al., 2001). These records were produced by the Central Weather Bureau
Seismic Network (CWBSN) and they are the largest set of strong motion data recorded froma major earthquake (Shin and Teng, 2001). Shin and Teng (2001) classified the recorded
accelerograms into four quality groups based on the existence of absolute timing, pre-events,
and defects. For this analysis, QA-class data (best for any studies) and QB-class data (nextbest but no absolute timing) are used.
Strong motion data from the Niigataken-chuetsu earthquake (October 23, 2004) wererecorded by the K-NET and KiK-net seismic networks operated by the National ResearchInstitute for Earth Science and Disaster Prevention in Japan. Those data are available
at their websites (http://www.k-net.bosai.go.jp/ and http://www.kik.bosai.go.jp/). The
stations with epicentral distance less than 100 km are used for this analysis.

2.2

Data Processing

We processed the accelerograms obtained from the nine earthquakes according to the following method. The DC offset of the accelerograms is corrected by subtracting the meanof the pre-event portion. Because a small DC offset has a large effect when the record is
integrated, this process is applied to all accelerograms.
The peak amplitude of the horizontal components is calculated by the square root of
the sum of the squares of the peaks of NS and EW components. If one of the horizontal
components (NS or EW) of a station has been clipped or is not well-recorded, the squareroot of twice the other well-recorded horizontal component is used for the peak amplitudeof the horizontal component.
The peak amplitude of UD (up-down) component is used directly for the peak verticalcomponent. The station records that have defects in the vertical component are excluded.
The following processes are completed for all the data.
Jerk: The three-component accelerograms are differentiated in the time domain, using
a simple finite-difference approximation. The peak value of each component is selected.
Acceleration: Original accelerograms are used to select the peak value.
Velocity: Some velocity records have a linear trend due to either tilting, the responseof the transducer to strong shaking, or problems in the analog-to-digital converter. The
baseline correction scheme applied to obtain appropriate velocity records is as follows (Iwanet al., 1985; Boore, 2001):
1) Determine the straight line to be subtracted from the velocity trace. The line is givenby the equation vf (t) = a1 t + a2 where coefficients a1 and a2 are determined by least-squares
fitting to the velocity trace after the strong shaking. The segment of the record used forleast-squares fitting is from t1 to t2 (see Figure 1). t1 is the time when the strong shaking
has subsided. The results of baseline correction are not very sensitive to the choice of t1(Boore, 2001). The second cut-off time, t2 , is generally chosen as the end of the record;
2) Remove this linear trend from the velocity record. The initial time to subtract thelinear trend is determined as the intersection between the linear trend and x-axis.
This baseline correction scheme assumes the baseline shift of the acceleration occurs

M. Yamada, T. Heaton, and J. Beck5

only once. There may be records that have more than one baseline shift during strong
shaking. However, our purpose is to get the peak value of each velocity record, and thisdoes not require accurate integration of the entire record. After time-domain integration,
the distortion is not very large in the first portion of the record where the peak value isgenerally recorded.
Displacement: The corrected velocity records are integrated once in the time domain
and high-pass filtered using a fourth-order Butterworth filter with a corner frequency of 0.075
Hz.
The peak features used for the classification analysis are shown in Table 2. Several
combinations of these 8 features are tried to find the best performance of the classification.

2.3

Data Classification

The classification as near-source or far-source in the training set is based on rupture areamodels used for waveform inversions. These rupture area models are typically determinedfrom the aftershock distribution (Sekiguchi et al., 1996), and the shape of the rupture area isapproximated by a rectangular box. Fault models used for classifying stations are shown in
Table 1 and Figure 2. In Figure 2, black solid lines indicate the surface projection of the fault
rupture surface based on the fault models. Stations within 10 km of this fault projection (thewhite area in the figures) are classified as near-source, indicated by solid circles. Far-source
stations are shown in open circles.
High-frequency near-source ground motions have long been researched by engineers andseismologists. High-frequency ground motions depend weakly on magnitude in the nearsource (Hanks and Johnson, 1976; Joyner and Boore, 1981; Hanks and McGuire, 1981).
This helps to analyze ground motions with a wide range of magnitude. Figure 3 shows
horizontal and vertical PGA of near-source records in our training set as a function of
moment magnitude. The slope of a regression line would be almost equal to zero, whichis consistent with past studies. On the other hand, low-frequency motion has a strong
correlation with magnitude. Figure 4 shows horizontal and vertical PGD as a functionof moment magnitude. The PGD are log-proportional to the magnitude. Based on suchobservations, we assume that high-frequency motion does not depend on magnitude forlarge earthquake and that accelerations do not exceed 2g, whereas low-frequency motion ishighly correlated with magnitude, and its amplitude increases as the magnitude becomeslarge.
High-frequency ground motion decays in amplitude more rapidly with distance thanlow-frequency motion (Hanks and McGuire, 1981). Therefore, high-frequency motions (e.g.
acceleration, jerk) have high correlations with the fault distance. We compute the log ofthe ground motion amplitudes and find the means and standard deviations for the nearsource and far-source records. Figure 5 shows the histograms and Gaussian densities givenby the sample means and standard deviations for the near-source and far-source records.
The Gaussian densities are good approximations of the histograms of the log of the groundmotion data. Figure 5 also shows that the distance between means for the near-source
and far-source datasets is larger in high-frequency than low-frequency motions. Therefore,

M. Yamada, T. Heaton, and J. Beck6

we expect that the high-frequency motions is a good measure to classify near-source andfar-source records.

3

Near-source versus Far-source Discriminant Function

We assume the discriminant function to classify records into near-source and far-source isexpressed as a linear combination of the log of ground motion amplitudes:f (Xi |θ) =c1 xi1 + c2 xi2 + ... + cm xim − dmX=
ck xik − d

(1)

k=1

=Xi · c − dwhere
xik = kth feature parameter of the ground motion at the ith stationm = the number of feature parameters
Xi =[xi1 , xi2 , ..., xim ]=[log10 (component1), log10 (component2), ..., log10 (componentm) ]c1 , ..., cm =the regression coefficientsd = decision boundary constantθ =[c1 , c2 , ..., cm , d]T
We may use m components out of the eight ground motion components shown in Table 2.
The coefficients c1 , ..., cm , and d are determined from the training data set by two different
approaches: Fisher’s linear discriminant analysis and Bayesian analysis.
This discriminant function is used to allocate new observations to one of the near-sourceor far-source groups, where f (Xi |θ) = 0 is the boundary between the two groups in thefeature parameter space. The station with observation Xi is classified as near-source if
f (Xi |θ) is positive. If f (Xi |θ) is negative, the station is classified as a far-source station.
Note that the decision boundary may also be expressed using equation (1) as: Xi · c = d.

3.1

Fisher’s Linear Discriminant Analysis

Fisher’s Linear Discriminant Analysis (LDA) is a method to classify data by using a linear
function (1) that best discriminates two or more naturally occurring groups. LDA was
first described by Fisher (1936) to separate two groups optimally. In general, LDA requires
placing objects (e.g. humans) in predefined groups (e.g. Caucasoid, Mongoloid, and Negroid)
based on certain feature parameters (e.g. related to physical characteristics), and finding afunction to distinguish the groups. The parameters ck in the linear function (1) are selected7

M. Yamada, T. Heaton, and J. Beckto minimize the within-group variance (variance of the samples centered on the group mean)
and maximize the between-group variance (variance between group means). The following isa brief discussion about the procedure of linear discriminant analysis (Venables and Ripley,2002):
Consider n × m data matrix X where n is the number of samples and m is thenumber of different features of samples. Each sample is assigned to one of ggroups Nj , j = 1, ..., g, with nj observations in each group. Let G denote thegroup indicator matrix, which indicates the group each sample is assigned to,
and let M denote the group mean matrix, then within-group covariance matrix
W and between-group covariance matrix B are:(X − GM)T (X − GM)n−g(GM − 1µ)T (GM − 1µ)
B=
g−1
W =

where

(2)
(3)

X =[xik ] : n × m data matrix
G =[gij ] : n × g group indicator matrix
M =[mjk ] : g × m group mean matrixµ =[µ1 , µ2 , ..., µm ] : 1 × m mean vector
1 =n × 1 column vector of 1s
xik = kth feature of the ith sample
gij =1 ⇐⇒ Xi = [xi1 , xi2 , ..., xim ] is assigned to group j
1 X
mjk =
xik
nj i∈N
j

µk =

n
1X

n i=1

xik

We would like to find a linear combination X ·c of the data such that the different
groups are maximally separated, that is, maximizing the following separation
ratio λ:λ=cT Bc
between-group variance=Tc Wc
within-group variance(4)= 0. By substituting equation (4)A necessary condition to maximize λ is ∂λ∂cinto this condition, we get:W −1 Bc = λc

(5)assuming W is invertible. This is an eigenvalue problem and the weight vector cand the separation ratio λ are eigenvectors and eigenvalues of W −1 B, respectively.
X · c is called a canonical variate, and the canonical variate of the eigenvector cwhich corresponds to the largest eigenvalue is called the first canonical variate.

M. Yamada, T. Heaton, and J. Beck

8

For the near-source versus far-source classification problem, the data matrix X is thedataset of peak seismic ground motions, where n is the number of stations, and m is numberof the object features (PGA, PGV, PGD, etc.). We have two groups: near-source groupand far-source group (g = 2). LDA finds the optimal set of coefficients of ground motion
amplitudes to classify near-source or far-source records.
Since the traditional LDA does not treat which choice of the ground motion parameters isthe best, Bayesian model class selection is performed (the results are shown later). According
to this analysis, the best selection is (Za and Hv), and their coefficients obtained from LDA
are shown in Table 3.
We choose the decision boundary constant d to maximize the classification performancefor the set of coefficients obtained by the LDA. The classification performance is evaluatedby the following function:
Pc (d) =(P (f (Xi |θ) ≥ 0|Yi = 1) + P (f (Xi |θ) < 0|Yi = −1))/2

(6)

where
f (Xi |θ) =Xi · c − d(1
if near-sourceYi =−1 if far-sourceThis is the average probability between the probability that a near-source station is classified
correctly and the probability that a far-source is classified correctly. The parameter d which
maximizes this function for the given coefficients (Table 3) is 25.903, and the performance
defined by the function above is 93.4%. Another way to compute d is to take the midpoint of
the two group means of the first canonical variate. This method makes it easier to computethe value of d and it gives d = 25.045, a good approximation to d = 25.903 which shows
maximum performance.
As a conclusion, the discriminant function computed from the LDA is:f (Xi |θ) =7.233 log10 Za + 6.813 log10 Hv − 25.903(7)(f (Xi |θ) ≥ 0 near-sourceiff (Xi |θ) < 0 far-sourceThis discriminant function is applied to all the dataset, and the values of f (Xi |θ) are
shown in Figure 6. The figure shows that most of the near-source data lie on the right side ofthe decision boundary, which means the classification performance is very good. Althougha fraction of the far-source records are misclassified, the misclassification of far-source datais less critical than that of near-source data.

3.2

Bayesian Approach

In this section, a Bayesian approach is applied to determine the coefficients of the discriminant function that classifies near-source and far-source data (Sivia, 1996; Jaynes, 2003). The

9

M. Yamada, T. Heaton, and J. Beckprobability density function (pdf) of parameter θ conditioned on data Dn and model class
M can be expressed using Bayes’ theorem:p(θ|Dn , M) ∝ p(Dn |θ, M) × p(θ|M)
posterior

likelihood

∝

n
Y
i=1

prior

P (Yi |Xi , θ) × p(θ|M)

(8)

where
θ =[c1 , c2 , ..., cm , d]T : parameter vector
Dn ={(Xi , Yi ) : i = 1, ..., n} : available set of data
Xi =[xi1 , xi2 , ..., xim ] : ground motion at the station i=[log10 (component1), log10 (component2), ..., log10 (componentm)]
1; if near-sourceat the station iYi =−1 ; if far-sourcem = the number of object featuresn = the number of data

Note that the model class M defines the likelihood for each value of θ in some set of values
and also the prior pdf p(θ).
We determine the parameters c1 , ..., cm , d based on a Bayesian approach using the same
notation as the LDA. The goal of the Bayesian approach is to obtain the posterior pdf of themodel parameters (θ) and determine the most plausible value of θ by maximizing this pdf.
Choice of Prior Distribution
Assume that the model class M is globally identifiable based on Dn (Beck and Katafygiotis, 1998), that is, there is a unique θ maximizing the likelihood p(Dn |θ, M). In this case,given a sufficiently large dataset Dn , the choice of prior pdf does not affect the resulting posterior pdf, and all posteriors with different priors will converge to the same answer (Sivia,
1996). Here, the prior is chosen to cover a wide range of the parameter space by selecting theprior of each model parameter to be a Gaussian pdf with zero mean and standard deviation
σ=100, so:
m
1
1
1 T
1 X 2
p(θ|M) = √c + d2 ))exp(− 2 θ θ) = √exp(− 2 (2σ
2σ k=1 k( 2πσ)m+1( 2πσ)m+1(9)

Choice of Likelihood function
Let the predictive probability that station i is near-source be P (Yi = 1|Xi , θ). The
predictive probability that a station is far-source is then P (Yi = −1|Xi , θ) = 1 − P (Yi =1|Xi , θ). A standard approach in Bayesian classification is to define the predictive probability
by applying the logistic sigmoid function φ(x) = 1/(1 + e−x ) to the linear function f (Xi |θ)that is also used in the traditional LDA (Li et al., 2002). The logistic sigmoid function isa smooth, positive, and monotonically increasing function, as shown in Figure 7. Although
there are other sigmoid functions which have these properties, the logistic sigmoid function

10

M. Yamada, T. Heaton, and J. Beckis mathematically convenient, and the class probability (shown in the Bayesian model classselection) is robust to the choice of sigmoid function for the following reason. Notice fromequation (1) that the location of the separating boundary f (Xi |θ) = 0 is independent of a
uniform scaling of the parameters. The Bayesian updating automatically produces a scaling
appropriate to the separation of the classes in the feature parameter space that is impliedby the data. If the data are well separated, a large scale will be chosen so that there is asteep transition in the class probability as the separating boundary is crossed; on the other
hand, if the data for the classes have significant overlap in the feature parameter space, thena smaller scale will be chosen to give a more gradual transition.
The predictive probability that the ith station is near-source is therefore defined here by:
P (Yi = 1|Xi , θ) = φ(f (Xi |θ)) =

1

(10)

1 + e−f (Xi |θ)As f (Xi |θ) becomes larger, the station is more likely to be near-source, and the probability
that the station is near-source becomes closer to one. Note that the predictive probability
that the station is far-source is then:P (Yi = −1|Xi , θ) = 1 − φ(f (Xi |θ)) = φ(−f (Xi |θ)) =1
1 + ef (Xi |θ)(11)where, from equation (1),f (Xi |θ) =

m
X
k=1

ck xik − d = Xi · c − dFrom equations (10) and (11), the likelihood function can be expressed as:p(Dn |θ, M) =

n
Y
i=1

P (Yi |Xi , θ) =

n
Y
i=1

φ(Yi f (Xi |θ)) =nY

1

1 + e−Yi f (Xi |θ)
i=1

(12)

Posterior Distribution
By substituting equations (9) and (12) into equation (8), the posterior can be expressed
as:
n
1 T Y
1
1
exp(− 2 θ θ)p(θ|Dn , M) ∝ √
−Y
m+1
2σ
1 + e i f (Xi |θ)( 2πσ)
i=1(13)

Both an asymptotic approximation and stochastic simulation are performed to characterize the pdf defined by equation (13). In the asymptotic approach, the posterior isrepresented by a Gaussian distribution for θ with mean θ̂, the most probable value of θ,and a covariance matrix Σ̂ defined later. Stochastic simulation uses the Metropolis algorithm to generate random samples of the parameter vector θ from the posterior pdf. It
is noted that it is computationally challenging to evaluate the proportionality constant inequation (13) that normalizes the posterior pdf because it requires numerical integration
over a high-dimensional parameter space. However, this evaluation can be avoided in boththe asymptotic approximation and stochastic simulation methods.

11

M. Yamada, T. Heaton, and J. Beck
3.2.1

Asymptotic Approximation

We first find the optimal value θ̂ of θ that maximizes the posterior pdf. This multidimensional
optimization problem is solved by a numerical optimization algorithm provided by Matlab.
Using Laplace’s method of asymptotic approximation, Beck and Katafygiotis (1998) showthat the posterior pdf for a set of model parameters θ for a globally identifiable model class
M (which has a unique most probable value) may be approximated accurately by a Gaussian
distribution with mean θ̂ and covariance matrix Σ̂, given a large amount of data. Define
H(θ) by:nY
H(θ) = −∇∇ log[p(Dn |θ, M)p(θ|M)] = −∇∇ log[ P (Yi |Xi , θ)p(θ|M)]

(14)

i=1

then Σ̂ = H(θ̂)−1 . By substituting equations (9) and (12) into equation (14);[H(θ)](α,β) =[−∇∇ log

n
Y
i=1

P (Yi |Xi , θ) − ∇∇ log p(θ|M)](α,β)

n
Y
∂
1
=−
(log
φi ) + 2 δαβ∂cα ∂cβ
σ
i=1
2

=−
=−
=

n
X

∂2
1
(log φi ) + 2 δαβ∂cα ∂cβ
σ
i=1

n
X
∂ 1
∂(Yi f (Xi |θ))1
[ φi (1 − φi )] + 2 δαβ∂cβ φi
∂cα
σ
i=1

n
X
i=1

φi (1 − φi )xiα xiβ +

1
δαβ
σ2

(15)

where φi = φ(Yi f (Xi |θ)), and equation (1), along with Yi2 = 1, has been used. The optimal
parameter values and their standard deviations for the selection of features Za and Hv areshown in Table 3. Note that for large σ, the effect of the prior in equation (15) is negligible.
In order to examine the sensitivity of the Bayesian approach to the training dataset, we
perform a cross-validation analysis. First, the training dataset is randomly divided into two
datasets and the discriminant function is constructed from one dataset (training set). This
discriminant function is applied to the other dataset (validation set) to check its classification
performance. We then switch the testing set and validation set, and repeat this crossvalidation analysis. We set the near-source / far-source boundary so that the probability isa half that the station is near-source, that is, the station is classified as near-source if theprobability that it is near-source is more than 1/2. The confusion matrices of these twoanalysis and the previous analysis which uses all of the dataset are shown in Table 4. The
classification error with half of the dataset is as small as that of the analysis which uses allof the dataset. Therefore, we confirm that the sensitivity to the training dataset is small,giving more confidence that the discriminant function from Bayesian analysis will perform
well for future earthquake data.

12

M. Yamada, T. Heaton, and J. Beck
3.2.2

Stochastic Simulation using Metropolis Algorithm

The asymptotic approximation is valid only if the posterior pdf for the model parameters canbe approximated well with a Gaussian distribution. This requires a large sample size and thatthe class of models M is globally identifiable based on data Dn (Beck and Katafygiotis, 1998).
On the other hand, a stochastic simulation algorithm can be applied to the problem which
generates samples from a Markov Chain whose stationary pdf is the posterior pdf, that is,
the samples are asymptotically distributed according to the posterior pdf for the parameters.
The Metropolis algorithm is used to solve this high-dimensional problem, because it does
not require evaluation of the normalizing constant for sampling the posterior pdf in equation(13).
The Metropolis algorithm is a Markov chain Monte Carlo (MCMC) method proposed byMetropolis et al. (1953). It is a simulation technique for generating random samples fromany given probability distribution. The algorithm uses a proposal pdf Q which depends onthe current sample of parameters, θ(t) at tth iteration (Mackay, 1998). Here, we choose as the
proposal density a Gaussian pdf centered on the current parameters θ(t) with the covariancematrix Σ of the parameters in the asymptotic approximation. The optimal parameters
estimated from direct optimization of the posterior pdf are used as an initial θ(t) . The
expression for Q is:
Q(θ′ |θ(t) ) =1(2π)m′ /2 |Σ|1/2

1
exp(− (θ′ − θ(t) )T Σ−1 (θ′ − θ(t) ))2(16)

where |Σ| is the determinant of the covariance matrix, and m′ is the dimension of theparameter θ, which is m + 1. A candidate sample is drawn from Q(θ′ |θ(t) ). The ratio of theposterior pdf in equation (8) at the current sample θ(t) and the candidate sample θ′ determines
whether to accept the candidate sample, according to:p(θ′ |Dn , M)r=p(θ(t) |Dn , M)
(
θ′
(t+1)
θ
=
θ(t)

with probability min(1, r)
with probability1 − min(1, r)

(17)

(18)

If r ≥ 1 then the candidate is accepted as the next sample in the Markov Chain. Otherwise,
the candidate is accepted with probability r as follows; we generate a random number uniformly distributed between zero and one, and if it is less than r, the candidate is accepted,
that is, θ(t+1) = θ′ . If it is not accepted, the current sample is repeated (θ(t+1) = θ(t) ).
This procedure is repeated until the desired number of samples are generated. There is aburn-in period at the beginning of the MCMC method until the probability distribution of
the current sample θ(t) is sufficiently close to the posterior pdf, which is the stationary pdf
of the Markov chain, so judgment is used to discard initial samples.
Figure 8 shows 5000 parameter samples generated with the Metropolis algorithm forthe optimal selection of features Za and Hv. This selection of the ground motion features
comes from Bayesian model class selection explained later. After discarding the samples

M. Yamada, T. Heaton, and J. Beck13

in the burn-in period (taken as the first 100 samples), the mean and standard deviation of
the samples are shown in Table 3. The average acceptance ratio of the candidate samples
θ′ is 44%, which indicates the method works well (Roberts et al., 1997). The stability ofthe sample mean and standard deviation of each parameter is examined in Figure 9. The
mean and standard deviation of the samples converge after the first 1000 samples are added.
The most probable values of the parameters from maximization of the posterior pdf are also
shown in Figure 9. Note that the means of the marginal pdf’s and the most probable valuesof the joint posterior pdf need not agree if these pdf’s are skewed.
The distribution of sample values for parameters θ and the resulting histogram of probability that a station is near-source calculated by the generated set of parameters are shown in
Figure 10. The distribution of parameter samples agrees well with the Gaussian distribution
defined by the optimal parameters and standard deviations given by the asymptotic approximation. The standard deviations of c1 and c2 are similar to each other and the distribution
is peaked close to the mean of the samples. The distribution of samples for the decisionboundary constant d has a standard deviation almost three times as large as that of c1 andc2 . However, in terms of coefficient of variation, the uncertainty in d is smaller than that ofother parameters (11.7% compared with 14.9% and 15.3% for c1 and c2 , respectively).
Figure 11 shows the correlation of samples of model parameters generated from thesimulation. This is the result of the model class with all parameters corresponding to theeight ground motion parameters, not the result of the optimal model class. The figure
shows that the parameter d is not correlated significantly with any other parameter. The
combinations of parameters which have significant interaction are horizontal and verticaljerk (c1 and c2 ), horizontal and vertical acceleration (c3 and c4 ), and horizontal and verticaldisplacement (c7 and c8 ). Parameters with the same component and similar frequency range(e.g., jerk and acceleration (c1 and c3 , and c2 and c4 ), acceleration and velocity (c3 andc5 , and c4 and c6 ), velocity and displacement (c5 and c7 , and c6 and c8 )) are also strongly
correlated. This result agrees with our intuition; horizontal and vertical components ofthe same quantity are correlated, and records with similar frequency ranges have similarattenuation relationships and so are correlated.

3.3

Comparison between Traditional LDA and Bayesian ApproachP
Parameters for the linear discriminant function f (Xi |θ) = m
k=1 ck xik − d are estimated bytraditional LDA and by the Bayesian approach with two different techniques to characterize
the posterior pdf. The results are shown in Table 3. The parameters for LDA are scaled
such that the norm of the vector c = [c1 , c2 ] is equal to the norm of the vector from theasymptotic approximation. Note that the discriminant function f (Xi |θ) is a linear function,
so for the traditional LDA, multiplying all ck and d by an arbitrary positive constant does
not change the result of classification. However, this is not true for the Bayesian approach,where the modulus of f (Xi |θ) affects the probability that a station is near-source.
The estimated parameters are close for the three methods. The coefficients from LDA
are within one standard deviation of those from both Bayesian methods, except that c1 from
LDA is slightly over one standard deviation from the corresponding mean and most probablevalues from the Bayesian methods.

14

M. Yamada, T. Heaton, and J. Beck

For the asymptotic approximation and Metropolis algorithm, the estimates and standard
deviations for the posterior parameter distribution are very close. If the posterior is a skewed
pdf, the mean is not necessarily equal to the most probable value (e.g. consider log-normal
distribution), as mentioned before. However, Figure 10 suggests that the posterior pdf isalmost symmetric, and the means of the samples and the most probable values should show
very good agreement. In this case, the Gaussian distribution is a good approximation forthe posterior pdf of the parameters.
By using the discriminant functions defined by the values of the parameters in Table 3, we
performed a classification analysis using the whole dataset. The classification performancefor the discriminant function from LDA and two Bayesian approaches are shown in Table 5.
The results for LDA show 100% of near-source data and 86% of far-source data are classified
correctly, and the result of Bayesian approach shows 78% of near-source data and 98% of farsource data are classified correctly. This discriminant function is the function which has thesmallest prediction error. To obtain this function, the misclassification of near-source data
and that of far-source data are considered to be of equal importance. Generally speaking,the misclassification of near-source data is more critical than that of far-source data, and we
may want to decrease the misclassification rate of near-source data. This misclassification
rate can be easily controlled by changing the decision boundary constant d. We also cancontrol this by shifting the near-source / far-source boundary in the Bayesian approach tocorrespond to some other probability than the 1/2 used in this classification analysis.
We performed the leave-one-out cross-validation to compare the misclassification rate
between LDA and the Bayesian method with asymptotic approximation. The idea of thismethod is to predict the probability of a station from the discriminant function constructedfrom the dataset from which that station is excluded. This process is repeated for all 695 data
and the accuracy of prediction is computed. The percentage of misclassified data is shown in
Table 6. It shows the prediction error of the Bayesian approach is much smaller than that of
LDA. In other words, the Bayesian approach is able to construct a more robust discriminantfunction. Therefore, we use the discriminant function obtained from the Bayesian methodwith asymptotic approximation for further analysis.

4

Bayesian Model Class Selection

4.1

MethodBayesian model class selection determines which combination of the eight ground motion
parameters gives the best classification for the near-source and far-source. The essential idea
is to find the most probable model class based on data Dn within a set of candidate model
classes Mj , j = 1, ..., J (Gull, 1988; Beck and Yuen, 2004). Applying Bayes’ theorem, theprobability of model class Mj can be expressed as follows:
evidenceprior

p(Dn |Mj )P (Mj |M)
P (Mj |Dn , M) =p(Dn |M)
normalizing constant

(19)

15

M. Yamada, T. Heaton, and J. Beckwhere
M ={M1 , M2 , ..., MJ } : a set of candidate model classes
J =number of the model classes

The left-side of equation (19) is the probability of a particular model class Mj given thedataset and a set of candidate model classes. On the right hand side, p(Dn |Mj ) is the evidence
for each model class, P (Mj |M) is the prior over the candidate model classes evaluated for
Mj , and p(Dn |M) is a normalizing constant given by:p(Dn |M) =

J
X
j=1

p(Dn |Mj )P (Mj |M)(20)

Assuming a uniform prior for the model class, P (Mj |M) in the numerator and denominatorof equation (19) cancel. By the total probability theorem, the evidence for Mj provided by
the dataset Dn is given as:Zp(Dn |Mj ) =p(Dn |θj , Mj )p(θj |Mj )dθj
(21)
θjThis is simply the integral of the likelihood of the data for a vector of parameters weighted
by its prior probability integrated over the whole parameter set for θj for model class Mj .
An asymptotic approximation for large sample sizes n can be used to compute the evidence of the model (Papadimitriou et al., 1997):
p(Dn |Mj ) ≈2π Nj /2 p(θ̂j |Mj )q
× p(Dn |θ̂j , Mj )likelihood|Hj (θ̂j )|(22)

Ockham f actor

where
Hj (θj ) = − ∇∇ log[p(Dn |θj , Mj )p(θj |Mj )]θ̂j = optimal parameter vector (most probable value) for model class Mj
Nj = number of parameters for model class Mj

Here, Hj (θj ) is given by equation (15) for the choice of parameters θj corresponding to model
class Mj . p(θ̂j |Mj ) is the prior defined in equation (9) and p(Dn |θ̂j , Mj ) is the likelihood
function defined in equation (12), evaluated at the optimal parameter vector for model class
Mj . For the model class selection results, the effect of the standard deviation of the Gaussianprior on the choice of most probable model class is examined later.

4.2

Results of Bayesian Model Class Selection

We used Bayesian model class selection to find the best combination of the eight groundmotion parameters with the same dataset as the previous classification problem. First,

16

M. Yamada, T. Heaton, and J. Beckwe impose the condition that both horizontal and vertical components be included in themodel for any selected ground motion quantity. Under this condition, there are four groupsof ground motion parameters (peak jerk, acceleration, velocity, and filtered displacement)
giving fifteen possible combinations. These fifteen candidate model classes are shown in
Table 7.
The results in Table 7 indicate that the combination of acceleration and velocity is themodel with highest probability, although the jerk and velocity combination also has significant probability. The log of prior (p(θ̂j |Mj )) is simply a function of Nq
j and becomes smaller
as the number of parameters increases. The factor p(θ̂j |Mj )(2π Nj /2 )/ |Hj (θ̂j )| in equation(22) is called the Ockham factor by Gull (Gull, 1988; Beck and Yuen, 2004). It penalizes
a more complicated model and so makes a simpler model preferable.
The Ockham factor isqalso shown in Table 7. Although the coefficient 2π Nj /2 andOckham factor, the effect of prior p(θ̂j |Mj ) is dominant.

|Hj (θ̂j )| are included in the

The log of the likelihood function p(Dn |θ̂j , Mj ) becomes larger as the number of theparameters in the model class increases because a more complicated model class will fit
the data better than a less complicated one. However, the Bayesian model class selection
automatically accounts for the trade-off between the complexity of the model (e.g. numberof parameters) and the fit of the data to find a well-balanced model (Beck and Yuen, 2004).
A useful information-theoretic interpretation of this trade-off is given in Muto and Beck(2007).
To examine the possible model classes further, the constraint that horizontal and verticalcomponents be used together is removed. We test all 255 model classes created from the
combinations of 8 parameters. The results for the best five model classes are shown in Table
8. The sum of the posterior probability of the five model classes is 95% out of all 255 model
classes.
Model class 1, which has the coefficients of the vertical acceleration and horizontal velocity, is the most probable model within the set of 255 model classes. The discriminantfunction for the most probable model in model class 1 is:
f (Xi |θ) =6.046 log10 Za + 7.885 log10 Hv − 27.091

(23)

where
P (Yi = 1|Xi , θ) =1
1 + e−f (Xi |θ)(24)

is the probability that station i is near-source. This result indicates that the amplitudeof high-frequency components is effective in classifying near-source and far-source stations.
Note that the probability that the station is near-source is higher, if f is larger.

4.3

Effect of the Choice of Prior

In this section, we examine the choice of prior for the parameters in the model class selection.
As we stated, for the Gaussian prior distribution, the effect of the number of parameters,

M. Yamada, T. Heaton, and J. Beck17

Nj , is significant if the prior standard deviation, σ, is large. We demonstrate this featureby performing model class selection with a Gaussian prior with different values of σ and a
uniform prior with different widths of boundary b. The posterior probabilities of the model
classes are shown in Table 9.
In the table, we can see the effect of the prior standard deviation in the Gaussian prior.
As we increase σ, it tends to bias the posterior probability towards simpler models (i.e.
models with less parameters). For example, the probability of model jav slightly decreases
as σ increases. The small probability of model jv with Gaussian prior (σ=10) is caused by
the narrow prior range. If σ is too small, it restricts the range of parameters as shown in
Table 10. Also, for the uniform prior case the results are very similar to the Gaussian priorwith σ=100. Based on these results, we judge that the choice of σ=100 for the Gaussianprior is a reasonable one for Bayesian model class selection in our classification application.

5

Results and Discussion

We apply the optimal discriminant function from Bayesian approach (in equations (23) and(24)) to all the stations in the dataset. Figure 12 shows the classification results. The
distribution of stations with a high probability of being in the near-source is consistent withthe fault geometry. As mentioned before, the fault models that are used here are thosefrom the source inversion, and they are not necessarily the best indicator of near-source andfar-source stations.
To examine the application for real-time analysis, the optimal discriminant function inequations (23) and (24) is applied to the Chi-Chi earthquake strong motion records. We
generated snapshots of the probability that a station is near-source from 10 seconds to
40 seconds after the beginning of rupture. Peak ground motions used for this classification
analysis are computed from the observed data every 10 seconds for each station and evaluatedin the discriminant function. The results are shown in Figure 13. A darker mark at a stationin Figure 13 indicates that the station is more likely to be near-source, and a lighter markindicates that the station is more likely to be far-source.
Ten seconds after the rupture initiation, the map shows that stations with high probabilityof being in the near-source are located near the epicenter, and it indicates that the rupturearea is propagating concentrically. At 20 seconds, the probability of being in the nearsource at thirteen stations is computed to be greater than 50 %, but the concentric station
distribution makes it difficult to identify any directivity of rupture propagation. The averageslip velocity is 2 km/s (Ji et al., 2003), and the rupture front propagates 40 km from the
hypocenter at this point. We can see the North-South character of the rupture direction
clearly after 30 seconds of rupture. At 40 seconds, the distribution of stations with high nearsource probability agrees with the fault surface projection, and stations at the near-source
and far-source boundary have around 50 % probability. Even though the fault geometries
used for the wave inversion are not necessarily the actual extent of the fault, to a first-orderapproximation, the classification results are in good agreement with them. The near-source
region at the north of the main rupture is a secondary rupture at the Shihtan fault, which issuggested by Shin and Teng (2001). This event may not be clear in the low frequency ground18

M. Yamada, T. Heaton, and J. Beckmotions, so it is not considered in the waveform inversion. However, the accelerograms atthat region are clearly larger than those of neighbor region, and the classification results
detected this secondary rupture.

6

Conclusion

We presented a methodology to classify seismic records into near-source or far-source records
as a prelude to estimating fault dimensions in an earthquake early warning system. Groundmotion records from some past earthquakes are analyzed to find a linear function that best
discriminates near-source and far-source records. Peak values of jerk, acceleration, velocity,and displacement are used in a traditional LDA and in a Bayesian approach to find the linear
combination of peak values which provides the best performance to classify near-source andfar-source records. All methods gave similar discriminant functions. We also analyzed which
combination of ground motion features had the best performance for classification using
Bayesian model class selection, and the best discriminant function is:
f (Xi |θ) =6.046 log10 Za + 7.885 log10 Hv − 27.091
P (Yi = 1|Xi , θ) =1
1 + e−f (Xi |θ)

(25)

(26)

where Za and Hv denote the peak values of the vertical acceleration and horizontal velocity,respectively, and P (Yi = 1|Xi , θ) is the probability that a station is near-source. Thisfunction indicates that the amplitude of high-frequency components is effective in classifying
near-source and far-source stations.
The probability that a station is near-source obtained using this optimal discriminantfunction for all the earthquakes shows the extent of the near-source area quite well, suggestingthat the approach provides a good indicator of near-source and far-source stations for realtime analyses. Note that this function is constructed by the training dataset with magnitude
greater than 6.5, so it only works for large earthquakes.

Acknowledgments
We thank Dr. David Wald of U.S. Geological Survey and Dr. Hiroo Kanamori of California
Institute of Technology for their help with obtaining data for the Hyogoken-nanbu earthquake. We
also appreciate Dr. Chen Ji for providing the fault model for the Chi-Chi earthquake. Members ofCaltech Earthquake Engineering group, especially Sai Hung Cheung gave us useful comments for
the Bayesian classification problem. Some of the figures are generated by Generic Mapping Tools(Wessel and Smith, 1991).

References
Akkar, S. and P. Gülkan (2002). A critical examination of near-field accelerograms from the sea ofMarmara region earthquakes, Bull. Seism. Soc. Am., 92, 428-447.

M. Yamada, T. Heaton, and J. Beck

19

Allen, R. M., and H. Kanamori (2003). The potential for earthquake early warning in Southern
California, Science 300, 786-789
Beck, J. L., and K. Yuen (2004). Model selection using response measurements: Bayesian probabilistic approach, J. Eng. Mech, 130, 192-203.
Beck, J. L., and L. S. Katafygiotis (1998). Updating models and their uncertainties. I: Bayesianstatistical framework, J. Eng. Mech., 124(4), 455-461.
Boore, D. M. (2001). Effect of baseline correction on displacement and response spectra for several
recordings of the 1999 Chi-Chi, Taiwan, earthquake, Bull. Seism. Soc. Am. 91, 1199-1211.
Campbell, K. W. (1981). Near-source attenuation of peak horizontal acceleration, Bull. Seism.
Soc. Am. 71, 2039-2070.
Cua, G. (2005). Creating the Virtual Seismologist: developments in earthquake early warning andground motion characterization. Ph.D. thesis, Department of Civil Engineering, California
Institute of Technology, 2005.
Cua, G. and T. H. Heaton (2006). The Virtual Seismologist method: a Bayesian approach toseismic early warning. In ”Seismic early warning”, eds: P. Gasparini, J. Zschau, Spring
Heidelberg. (in press)
Fisher, R. A. (1936). The use of multiple measurements in taxonomic problems, Ann Eugen 7,
179-188.
Gull, S. (1988) Bayesian Inductive Reference and Maximum Entropy, in Maximum Entropy andBayesian Methods in Science and Engineering, C.J. Erickson and C.R. Smith eds., 53-74.
Hanks, T. C., and D. A. Johnson (1976). Geophysical assessment of peak accelerations, Bull.
Seism. Soc. Am., 66, 959-968
Hanks, T. C. and R. K. McGuire (1981). The Character of High-frequency Strong Ground Motion,
Bull. Seism. Soc. Am., 71, no.6, 2071-2095.
Hartzell, S., and T. Heaton. (1983). Inversion of strong ground motion and teleseismic waveform
data for the fault rupture history of the 1979 Imperial Valley, California, earthquake, Bull.
Seism. Soc. Am. 73, 1553-1583.
Honda, R., S. Aoi, N. Morikawa, H. Sekiguchi, T. Kunugi, and H. Fujiwara (2005). Ground Motion
and Rupture Process of the 2004 Mid Niigata Prefecture Earthquake Obtained from Strong
Motion Data of K-NET and KiK-net, Earth Planets Space, (submitted).
Iwan, W. D., M. A. Moser, and C.-Y. Peng (1985). Some observations on strong-motion earthquake
measurement using a digital accelerograph, Bull. Seism. Soc. Am. 75, 1225-1246.
Jaynes, E. T. (2003). Probability Theory: The Logic of Science, Cambridge University Press.
Ji, C., D. V. Helmberger, D. J. Wald and K. F. Ma (2003). Slip history and dynamic implicationof 1999 Chi-Chi earthquake. J. Geophys. Res., 108(B9), 2412, doi:10.1029/2002JB001764.
Joyner, W. B. and D. M. Boore (1981). Peak horizontal acceleration and velocity from strongmotion records including records from the 1979 Imperial Valley, California, earthquake Bull.
Seism. Soc. Am., 71, 2011-2038.
Lee, W. H. K., T. C. Shin, K. W. Kuo, K. C. Chen, and C. F. Wu (2001). CWB Free-Field
Strong-Motion Data from the 21 September Chi-Chi, Taiwan, Earthquake, Bull. Seism. Soc.
Am., 91, no.5, 1370-1376.
Li, Y., C. Campbell, and M. Tipping (2002). Bayesian automatic relevance determination algorithms for classifying gene expression data. Bioinformatics, 18, 1332-1339.
MacKay, D. J. C. (1998). Introduction To Monte Carlo Methods. In: Learning in GraphicalModels, MIT Press, 175-204.
Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., and Teller, A. H. (1953). Equations ofstate calculations by fast computing machines. J. Chem. Phys., 21(6), 1087-1092.
M. Muto and J. L. Beck (2007). Bayesian Updating and Model Class Selection for Hysteretic
Structural Models Using Stochastic Simulation. Journal of Vibration and Control (in print).
Nakamura, Y. (1988). On the urgent earthquake detection and alarm system (UrEDAS), Proc. of
the 9th World Conference on Earthquake Engineering, Tokyo-Kyoto, Japan. VII. 673-678.
Odaka, T., K. Ashiya, S. Tsukada, S. Sato, K. Otake, and D. Nozaka (2003). A New Method of
Quickly Estimating Epicentral Distance and Magnitude from a Single Seismic Record, Bull.
Seism. Soc. Am., 93, no.1, 526-532
Papadimitriou, C., Beck, J. L., and Katafygiotis, L. S. (1997). Asymptotic expansions for reliabilities and moments of uncertain dynamic systems, J. Eng. Mech., 123(12), 1219-1229.
Roberts, G. O., A. Gelman, and W.R. Gilks (1997). Weak convergence and optimal scaling of

M. Yamada, T. Heaton, and J. Beck20

random walk Metropolis algorithms, Ann. Appl. Prob. 7, 110-120.
Sekiguchi, H., K. Irikura, T. Iwata, Y. Kakehi, and M. Hoshiba (1996). Minute locating of faulting
beneath Kobe and the waveform inversion of the source process during the 1995 Hyogo-ken
Nanbu, Japan, earthquake using strong ground motion records, J. Phys. Earth, 44, 473-488.
Sekiguchi, H., and T. Iwata (2002). Rupture process of the 1999 Kocaeli, Turkey, earthquake
estimated from strong-motion waveforms, Bull. Seism. Soc. Am., 92, 300-311.
Shin, T.-C. and T.-L. Teng (2001). An Overview of the 1999 Chi-Chi, Taiwan, Earthquake Bull.
Seism. Soc. Am., 91, 895-913.
Sivia, D. S. (1996). Data Analysis: A Bayesian Tutorial, Oxford University Press, Oxford.
Toki, K., K. Irikura, and T. Kagawa (1995). Strong motion data recorded in the source area of the
Hyogo-ken-nanbu earthquake, January 17, 1995, Japan, J. Nat. Disast. Sci., 16, 23-30.
Tsuboi, S., D. Komatitsch, C. Ji, and J. Tromp (2003). Broadband modeling of the 2002 Denali
fault earthquake on the Earth Simulator, Phys. Earth Planet. Interiors 139, 305-312.
Venables, W. N. and B. D. Ripley (2002). Modern Applied Statistics with S. Fourth Edition,
Springer, New York, NY.
Wald, D. J. (1996). Slip history of the 1995 Kobe, Japan, earthquake determined from strong
motion, teleseismic, and geodetic data, J. Phys. Earth, 44, 489-503.
Wald, D. J., T. H. Heaton, and D. V. Helmberger, (1991), Rupture model of the 1989 Loma Prieta
earthquake from the inversion of strong motion and broadband teleseismic data, Bull. Seism.
Soc. Am, 81, 1540-1572.
Wald, D. J., and T. H. Heaton (1994). Spatial and temporal distribution of slip for the 1992
Landers, California, earthquake, Bull. Seism. Soc. Am. 84, 668-691.
Wald, D. J., T. H. Heaton, and K. W. Hudnut (1996). A dislocation model of the 1994 Northridge,
California, earthquake determined from strong-motion, GPS, and leveling-line data, Bull.
Seism. Soc. Am. 86, S49-70.
Wessel, P., and W. H. F. Smith (1991). Free software helps map and display data, EOS. Trans.
AGU 72, 441, 445-446.
Wu, Y. M. and H. Kanamori (2005). Experiment on an onsite early warning method for the Taiwan
early warning system, Bull. Seism. Soc. Am., 95, 347-353. 2005.

21

M. Yamada, T. Heaton, and J. Beck

Figures and Tables
Earthquake
Mw
Imperial Valley (1979)
6.5
6.9
Loma Prieta (1989)
Landers (1992)
7.3
Northridge (1994)
6.6
Hyogoken-Nanbu (1995)
6.9
Izmit (1999)
7.6
Chi-Chi (1999)
7.6
Denali (2002)
7.8
Niigataken-Chuetsu (2004) 6.6
Total

NS
14
8
1
17
4
4
42
1
9
147

FS Total
Fault Model
20
34
Hartzell and Heaton, 1983
39
47
Wald et al., 1991
112 113
Wald and Heaton, 1994
138 155
Wald et al., 1996
14
18
Wald, 1996
13
17
Sekiguchi and Iwata, 2002
172 214Ji et al., 2003
29
30
Tsuboi et al., 2003
58
67
Honda et al., 2004
623 770

Table 1: The earthquake dataset used for the classification analysis. Moment magnitude (Mw )is cited from Harvard CMT solution. The numbers of near-source (NS) and far-source (FS) data
for each earthquake are also shown. The fault models are used as selection criteria to classify
near-source and far-source stations.

Code
Hj
Zj
Ha
Za
Hv
Zv
Hd
Zd

Measurement
Horizontal Peak Ground Jerk
Vertical Peak Ground Jerk
Horizontal Peak Ground Acceleration
Vertical Peak Ground Acceleration
Horizontal Peak Ground Velocity
Vertical Peak Ground Velocity
Horizontal Peak Ground Displacement
Vertical Peak Ground Displacement

Unit(cm/s3 )(cm/s3 )(cm/s2 )(cm/s2 )(cm/s)(cm/s)(cm)(cm)

Table 2: Eight measurements of peak ground motions are calculated from three component accelerograms. Codes and units of the components used in this paper are shown.

22

M. Yamada, T. Heaton, and J. Beck
MethodLDA
Bayesian-Asym.
(σ)
Bayesian-MA(σ)

c1 (Za)c2 (Hv)
d
7.233
6.813
25.903
6.046
7.886
27.090(± 0.903) (± 1.206) (± 3.163)
6.194
8.150
27.872(± 0.946) (± 1.224) (± 3.330)

Table 3: Estimated model parameters by Fisher’s LDA, Bayesian approach with asymptotic approximation, and Bayesian approach with Metropolis algorithm. The standard deviations for eachparameter are shown in brackets.

DatasetAll dataset
Half of dataset
Other half of dataset

NS/FS
NS
FS
NS
FS
NS
FS

Near-source78 (78%)12 (2%)
39 (74%)4 (1%)
37 (79%)8 (3%)

Far-source
22 (22%)583 (98%)
14 (26%)291 (99%)
10 (21%)292 (97%)

Table 4: The confusion matrix for the cross-validation analysis with the Bayesian method withasymptotic approximation. “All dataset” is the analysis which uses the whole dataset as a training
set and a validation set. “Half of dataset” is the analysis which uses half of dataset as a training
set and the other half as a validation set. “Other half of dataset” is the analysis which switches
the training and validation set. NS and FS stand for near-source and far-source, respectively.

DatasetLDA
Bayesian-Assym.
Bayesian-MA

NS/FS
NS
FS
NS
FS
NS
FS

Near-source
100 (100%)
82 (14%)78 (78%)12 (2%)
78 (78%)12 (2%)

Far-source
0 (0%)
513 (86%)
22 (22%)583 (98%)
22 (22%)583 (98%)

Table 5: The confusion matrix for near-source versus far-source classification by the discriminantfunction obtained from LDA, Bayesian approach with asymptotic approximation, and Bayesian
approach with Metropolis algorithm.

MethodLDA
Bayesian approach

Prediction Error
82 / 695 (12%)
36 / 695 (5%)

Table 6: Results of leave-one-out cross-validation for LDA and Bayesian approach.

23

M. Yamada, T. Heaton, and J. Beck

Model
j
a
v
d
ja
jv
jd
av
ad
vd
jav
jad
jvd
avd
javd

Hj
1.53
-2.74
2.57
3.44
1.36
0.55
2.72
1.40

Zj
4.30
2.04
2.79
3.43
1.47
0.43
2.68
1.29

Ha
4.38
6.60
2.54
4.93
1.36
4.35
3.47
2.05

Za
4.37
2.95
4.38
5.02
2.28
4.49
4.50
2.49

Hv
8.57
7.00
7.01
12.55
6.93
6.66
4.58
5.05

Zv
0.87
2.00
0.91
2.30
1.50
2.91
1.06
2.11

Hd
2.49
3.48
3.89
-3.38
3.89
0.66
1.80
1.69

Zd
1.44
0.79
0.22
-0.25
0.27
-1.12
-0.47
-1.02

d
23.84
21.43
16.33
5.76
20.82
36.09
33.17
29.47
29.40
19.99
33.75
30.72
36.66
30.16
34.31

Ock.
-17
-16
-16
-17
-25
-25
-26
-24
-25
-25
-33
-33
-34
-33
-41

Likeli.
-140
-117
-118
-192
-114
-80
-94
-80
-82
-106
-78
-81
-80
-79
-78

Evi.
-156
-133
-134
-209
-139
-105
-120
-104
-106
-131
-111
-115
-113
-112
-119

Prob.
0.00
0.00
0.00
0.00
0.00
0.32
0.00
0.62
0.05
0.00
0.00
0.00
0.00
0.00
0.00

Table 7: Results for Bayesian model class selection when fifteen combinations of the ground motion
parameters are examined under the condition that the horizontal and vertical components are used
together. The most probable value of the decision boundary parameter corresponding to eachground-motion parameter is given first for each model class. The values for the Ockham factor(Ock.), likelihood (likeli.), and evidence (evi.) of each model class are log-scaled. The last column
is the posterior probability that measures how plausible the model class is. It is scaled such thatthe total probability of the fifteen model classes is 1.0.

Model
1
2
3
4
5

Hj
Zj
Ha
Za
Hv Zv Hd
6.05 7.89 1.91
4.41 8.31 1.86 4.88 7.86 1.59
4.31 8.02 4.43
8.52 -

Zd
-

d
Ock.
27.09 -15
31.92 -20
29.17 -20
29.71 -20
32.22 -16

Likeli.
-81
-79
-80
-80
-84

Evi. Prob.
-96
0.81
-99
0.07
-100 0.03-100 0.03-100 0.02

Table 8: The best five model classes in the Bayesian model class selection when 255 combinationsof the ground motion parameters are examined. The columns are in the same format as in 7.

24

M. Yamada, T. Heaton, and J. Beck

Model
j
a
v
d
ja
jv
jd
av
ad
vd
jav
jad
jvd
avd
javd

Gaussian prior
σ=10
σ=100
σ=1000
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
7.2
32.4
33.0
0.0
0.0
0.0
78.9
62.1
61.7
7.3
5.3
5.3
0.0
0.0
0.0
3.3
0.1
0.0
0.1
0.0
0.0
0.1
0.0
0.0
3.0
0.0
0.0
0.1
0.0
0.0

Uniform prior|b| <20|b| <100
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
31.5
32.9
0.0
0.0
59.0
61.6
5.0
5.3
0.0
0.0
3.0
0.1
0.0
0.0
0.3
0.0
1.1
0.0
0.0
0.0

Table 9: The posterior probability of the model class selection with different types of prior distribution for parameters. σ is the standard deviation for the Gaussian distribution and |b| is thewidth of the boundary for the uniform distribution.

Prior
c1 (Za) c2 (Hv)
Gaussian(σ=10)
5.522
7.147
Gaussian(σ=100)
6.046
7.885
Gaussian(σ=1000) 6.053
7.895
6.053
7.895
Uniform Casesd24.68627.091
27.12227.122

Table 10: The estimated parameters from Bayesian approach with different types of prior distribution for parameters.

25

M. Yamada, T. Heaton, and J. Beck
150

t2

100
t

(cm/s)

1

50

0

−50

−100
0

Original velocity
Corrected velocity
Linear trend
50

100

150

(s)

Figure 1: An example of baseline correction for a velocity record from the Chi-Chi earthquake.
The corrected velocity trend is obtained by subtracting the linear trend from the original velocity
record. The portion of the record from t1 to t2 is used for least square fitting to obtain the linear
trend.

26

M. Yamada, T. Heaton, and J. Beck(a) Imperial Valley (1979)(b) Loma Prieta (1989)(c) Landers (1992)(d) Northridge (1994)(e) Hyogoken-Nanbu (1995)(f) Izmit (1999)

Figure 2: Maps the fault projections and station distributions. The fault projections are shown inthe solid lines. The white area around the fault lines indicates the area with distance less than 10

27

M. Yamada, T. Heaton, and J. Beck(g) Chi-Chi (1999)(h) Denali (2002)(i) Niigataken-Chietsu (2004)

Figure 2: Maps of the fault projections and station distributions (continued).

28

M. Yamada, T. Heaton, and J. Beck

4

10

Horizontal
VerticalPGA (cm/s2)

3

10

2

10

1

10

6

6.5

7

7.5

8

Figure 3: Distribution of horizontal and vertical PGA for near-source stations with respect tomagnitude.

3

10

Horizontal
Vertical2

PGD (cm)

10

1

10

0

10

6

6.5

7

7.5

8

Figure 4: Distribution of horizontal and vertical PGD for near-source stations with respect tomagnitude.

29

M. Yamada, T. Heaton, and J. Beck

Horizontal

Vertical22

far−sourcenear−source

1
0.5
0

near−source
far−source

1.5

pdf

pdf

1.5

1
0.5

2

3

4

5

0

6

2

3

log(jerk)

pdf

2
1.5

pdf

2
1.5
1
0.5
0

6

3

4

2

3

2

3

1

0

1

2

3

0

4

0

1

2

log(acc)
2

1.5

1.5

pdf

2

pdf

5

0.5

log(acc)

1
0.5

1
0.5

0
−1

0

1

2

0
−1

3

0

log(vel)1log(vel)
2

1.5

1.5

pdf

2

pdf

4

log(jerk)

1
0.5

1
0.5

0
−1

0

1

2

0
−1

3

0

log(dis)1log(dis)Figure 5: Histograms and Gaussian densities based on the sample means and standard deviationsof the log of ground motions for the near-source and far-source records. These are distributions forjerk, acceleration, velocity, and displacement from the top.

2
1.8

near−source
far−source

1.6
Far−source

Near−source

1.4

pdf

1.2
1
0.8
0.6
0.4
0.2
0
−2

−1.5

−1

−0.5

0
f(X |θ)

0.5

1

1.5

2

i

Figure 6: Histogram of the near-source and far-source data to which the discriminant function
obtained from traditional LDA is applied. The column heights are normalized by the number ofthe data in each group. f (Xi |θ) = 0 is the decision boundary between the two groups. The curves
are the Gaussian distribution with the same mean and standard deviation as the values of f (Xi |θ)for each group.

30

M. Yamada, T. Heaton, and J. Beck1−xφ(x)=1/(1+e )
0.8

φ(x)

0.6

0.4

0.2

0
−10

−5

0
x

5

10

Figure 7: Logistic sigmoid function φ(x) = 1/(1 + e−x ) is used to express the predictive probability
for classification. The function approaches zero as x → -∞, and one as x → ∞. The function is
0.5 when x is zero.

c (Za)

20

1

10

0
0

1000

2000

3000

4000

5000

1000

2000

3000

4000

5000

1000

2000

3000

4000

5000

10

2

c (Hv)

20

0
0
0

d

−20
−40
−60
0

Figure 8: Samples generated by Metropolis algorithm plotted in the parameter space. The x-axis
denotes the sample number. The vertical dotted lines indicate the end of the burn-in period (100samples).

31

M. Yamada, T. Heaton, and J. Beck61c (Za)

8

4
0

1000

2000

3000

4000

5000

1000

2000

3000

4000

5000

1000

2000

3000

4000

5000

c (Hv)

10

2

8
6
4
0
−20

d

−25
−30
−35
0

Figure 9: Mean and standard deviation of samples plotted against the number of samples included(excluding first 100 samples). The solid line is the sample mean, and the dashed lines represent
the mean plus and minus one standard deviation. The small circle is the most probable values ofthe model parameters estimated from optimization.

0.5

1c (Za)

1

5

10

15

20

25

30

5

10

15

20

25

30

15

20

25

30

35

40

0.5

2

c (Hv)

0
0
1

d

0
0
1
0.5
0
10

Figure 10: Distribution of samples for 3 parameters generated by Metropolis algorithm. The
Gaussian distributions obtained from the asymptotic approximation are added in the figure, andfit the histogram well.

32

M. Yamada, T. Heaton, and J. Beckc2(Zj)c3(Ha)

10

10

−0.57
0
−10
−10

c4(Za)

10
−0.76
0
0

−10
10 −1010−10
10 −10
10
0.47

0
−10
−10

0
0

−10
10 −10
10

0

0
−10
−10

0

−10
10 −20
10
−0.69
0
−10
10 −20
10
−0.46
0
−10
10 −2010

0−10
20 −10
10
−0.18
0
−10
20 −10
10
−0.42
0
−10
20 −10
10

0

−10
−20

−10
20 −1020−10
10 −10
10
−0.09
0
−10
10 −10
10

0

0

0

0

−10
10 −10
10
−0.52
0
−10
10 −10
20
−0.39
0
−20
10 −1010−10
10 −6010

0−10
10 −1010−10
10 −60
10

0

0

−10
10 −10
20
−0.71
0
−20
10 −1010−10
10 −10
10

−40

0.36

0

0

0

−20
10 −60
10
−0.77
0
−10
10 −60
10
−0.71
0
−10
10 −60
10

−40

0

−40

c7(Hd)

−20
0.13

−10
−60

c6(Zv)

−20
−0.21

−40

c5(Hv)

−20
−0.25

−40

c4(Za)

−20
−0.06

0

0
−10
−10

−10
10 −6020

0c3(Ha)

−20
0.12

0

0

0

−40

0.13

c2(Zj)

−20
−0.12

0

0

0

−40

0.12

c1(Hj)

−20
−0.31

0

0.13

0

−40

0.11

0.4

−10
−10

−10
10 −60
10

0

0
0

−0.25
0

0.11
0

0
−20
−10

−10
10 −10
10
0.17

0
0

0

0
0

−0.31
0

0.22

0.1
0

−10
10 −10
10

d
10

0.14
0

0

c8(Zd)
10

0.12
0

0

c7(Hd)
10

0.22
0

0

c6(Zv)
10

0.27
0

0

c5(Hv)10c8(Zd)−20

Figure 11: Correlation plot of posterior samples of the model parameters generated by Metropolis
algorithm. The most probable values of the parameters are shown as ”x”. The numbers in the
figure are the correlation coefficient of parameters.

33

M. Yamada, T. Heaton, and J. Beck(a) Imperial Valley (1979)(b) Loma Prieta (1989)(c) Landers (1992)(d) Northridge (1994)(e) Hyogoken-Nanbu (1995)(f) Izmit (1999)

Figure 12: Probabilities of near-source based on the optimal discriminant function obtained by the
Bayesian approach. Darker marks have higher probability that the station is located at near-source.

34

M. Yamada, T. Heaton, and J. Beck(g) Chi-Chi (1999)(h) Denali (2002)(i) Niigataken-Chietsu (2004)

Figure 12: Probabilities of near-source based on the optimal discriminant function obtained by the
Bayesian approach (continued).

35

M. Yamada, T. Heaton, and J. Beck(a) 10 seconds(b) 20 seconds(c) 30 seconds(d) 40 seconds

Figure 13: Snapshots of the probabilities of near-source for the Chi-Chi earthquake, based on
the optimal discriminant function from the Bayesian approach. The large circle is the theoreticalrupture front assuming the rupture velocity 2km/s.