Liu and Yamada

1

Bayesian Approach for Identification of Multiple Events in an
Early Warning System
Annie Liu∗

Masumi Yamada†

1

Abstract2

The 2011 Tohoku earthquake (Mw9.0) was followed by a large number of aftershocks that3

resulted in 70 early warning messages in the first month after the main shock. Of these4

warnings, a non-negligible fraction (63%) were false warnings where the largest expected5

seismic intensities were overestimated by at least two intensities or larger. These errors can

6

be largely attributed to multiple concurrent aftershocks from distant origins that occur within

7

a short period of time. Based on a Bayesian formulation that considers the possibility of8having more than one event present at any given time, we propose a novel likelihood function

9suitable for classifying multiple concurrent earthquakes, which uses amplitude information.

10

We use a sequential Monte Carlo heuristic whose complexity grows linearly with the number

11

of events. We further provide a particle filter implementation and empirically verify its

12

performance in with the aftershock records after the Tohoku earthquake. The initial case∗†aliu@cms.caltech.edu, Computer Science, California Institute of Technology, Pasadena, CA, USA
masumi@eqh.dpri.kyoto-u.ac.jp, Disaster Prevention Research Institute, Kyoto University, Uji, Japan

Liu and Yamada

213

studies suggest promising performance of this method in classifying multiple seismic events

14

that occur closely in time.

Liu and Yamada

315

Introduction

16

During the highly seismically active period after a major earthquake, multiple earthquakes

17

can occur almost concurrently at different locations. In this case, the seismic waves measured

18

by the ground sensors contain mixed signals from more than one source. If the detection

19

algorithm assumes only one quake, the estimated quake parameters (e.g. location and mag-

20

nitude) will not be accurate. These inaccurate estimates can lead to false warnings that are21

often observed after large earthquakes.

22

The 2011 off the Pacific Coast of Tohoku Earthquake (hereafter called Tohoku earth-23

quake) caused significant damage over a large area of northeastern Honshu. An earthquake

24

early warning (EEW) was issued to the public in the Tohoku region about 8 s after the25

first P-arrival, which is 31 s after the origin time (Hoshiba et al., 2011; Sagiya et al., 2011;26

Hoshiba and Iwakiri, 2011). There was no blind zone, i.e., warnings were received at all

27

locations before the S-arrivals, since the earthquake was fairly far offshore.

28

The main earthquake was followed by a large number of aftershocks that resulted in 70

29early warnings issued in the first month after the main shock (JMA, 2011). Among these,30

63% of the warnings contained significant errors where the estimated seismic intensities were31

at least two scales larger than the observed ones. As a comparison, only 29% of the warnings

32

contain such errors prior to the Tohoku earthquake. Post-event analysis revealed that 73% ofLiu and Yamada

433

these errors could be attributed to failure to classify multiple concurrent quakes either from34

the same hypocenter separated by a short amount of time or from spatially distant origins

35(JMA, 2011). One of the main reasons of these false alarms is that the current approach36

uses mainly P-wave arrival time to estimate the hypocenter.

37

In this paper, we propose a novel approach to detect and classify multiple concurrent

38

earthquakes in the current Japan Meteorological Agency (JMA) system framework. We

39

introduce an approximate Bayesian method that estimates the location, magnitude, and40

origin time of multiple concurrent aftershocks. In contrast to the current JMA system,

41

this approach produces multiple sets of estimation for earthquakes that occur closely in

42time. The experimental results from several case studies suggest that this approach can

43

successfully detect and estimate the parameters of multiple concurrent earthquakes.

44

Data and Processing

45

This paper includes strong motion data observed by the JMA seismic stations during and46

after the Tohoku earthquake. We evaluate the new classification approach on three sections

47of these records as summarized below. For each record, the values included in the JMA

48

EEW are compared to the values that appear in the JMA catalog in Table 1.

Liu and Yamada

549

Data set

50

Case 1: 15 March 2011, 1:36:00 - 1:38:00 (two small earthquakes). Early warnings

51

were issued to the public based on an estimated JMA magnitude of 5.9 at 21 s after the first52

P-wave detection (see Data and Resources). However, the largest observed seismic intensity53

was only 2 in the JMA seismic intensity scale. As shown in Table 1, at least two events

54

about 200 km apart of magnitude 2.5 and 3.3 occurred within 15 s. Since the second event

55

started close in time to the wave arrivals of the first event, the EEW system treated these56

separate events as one single earthquake and as a result, overestimated the magnitude.

57

Case 2: 20 March 2011, 14:19:00 - 14:21:00 (two small earthquakes).

58

EEW system estimated a magnitude of 7.6 at 6.6 s after the initial detection of P-wave and59

issued a warning to the public (see Data and Resources). However, the largest observed

60

seismic intensity was only 3. Again, as shown in Table 1, the overestimation can possibly61

be attributed to mistaking two smaller quakes about 150 km apart (M3.0 and M4.7) that62

occurred within 5 s for one large quake, since the occurrence of the second event was close

63

in time to the wave arrival of the first event.

64

Case 3: 11 March 2011, 14:46:00 - 14:49:00 (Tohoku earthquake). To demonstrate

65

that the method can also handle the classification of a single event, we also include the

The JMA

Liu and Yamada

666

analysis of the Tohoku earthquake (Mw9.0). An early warning was issued to the public in67

the Tohoku region about 8 s after the first P-arrival, which is 31 s after the origin time (see68

Data and Resources).

69

Processing

70

This paper uses the three component acceleration data with a sampling rate of 100 Hz

71

from about 200 stations. The acceleration data was first converted to SAC format and72

decimated by a factor of 100, reducing the sampling frequency to 1 Hz. The decimation was73

not necessary but was used to reduce computation time. Each component of the decimated74

acceleration k(t) was then converted to displacement A(t). The conversion was done by twice75

integration of k(t) using a recursive digital filter with the frequency response of a mechanical76

seismometer (Katsumata, 2008).

A(t) = gn × [k(t) + h0 · k(t − 1) + h1 · k(t − 2)] − h2 · A(t − 1) − h3 · A(t − 2),

(1)

77

where the function gain factor gn and filter constants h0 , h1 , h2 , h3 depend on the sampling78

frequency, damping constant, and natural period of the seismometer. For a JMA seismometer

79

with 100Hz sampling, 0.55 damping constant, and 6 s natural period, the values correspond

80

to:

Liu and Yamada

7

gn = 0.0000248691025, h0 = 1.0, h1 = 1.0, h2 = −1.9889474, h3 = 0.9895828.

(2)

81

The following approach to classification uses both the vector sum of the three component

82

displacement A(t) as well as the vertical component of acceleration k(t). The picking is83

done with STA/LTA of k(t) with a short term window of 1 s and long term window of 10

84

s. The method also computes expected P- and S-wave arrival times (tp and ts ) to determine85

whether a station should have observed P-wave or S-wave or neither. These arrival times

86

are computed with the JMA 1D layered velocity structure (Ueno et al., 2002).

87

Bayesian Method88

The problem of continuous parameter estimation for multiple events can be formulated as89

a Bayesian inference problem. Let θ be the vector of parameters that characterizes an

90

event and Θ be a set of events that are parametrized by θ’s, Θ = {ø, {θ1 }, . . . , {θ1 , θ2 , . . .}}.

91

Suppose z1:t is the complete history of observations from all the stations till the current time

92

t, the posterior P (Θt |z1:t ) reveals the distribution of information of current ongoing events

93

at time t given the evidence and prior information.

P (Θt |z1:t ) =P (zt |Θt ) P (Θt |z1:t−1 ),P (zt |z1:t−1 )(3)Liu and Yamada

894

where P (zt |Θt ) is the likelihood function and is typically denoted as L, L(zt |Θt ) = P (zt |Θt ).

95

P (Θt |z1:t−1 ) is the updated prior at time t,

Z
P (Θt |z1:t−1 ) =P (Θt |Θt−1 )P (Θt−1 |z1:t−1 ) dΘt−1 ,

(4)

96

and P (Θ0 |z0 ) ≡ P (Θ0 ) is the prior distribution of Θ.

97

Particle Filter

98

In general, Equation (3) does not have a closed-form solution, and there exists several sub-99

optimal solutions to approximate the posterior distribution (Arulampalam et al., 2002), one

100of which is grid search. Grid search, though simple to implement, suffers a few problems.

101

First of all, when the parameters are continuous and not sufficiently restricted, the method

102

cannot cover the complete parameter space since there can only be a finite number of grids.

103

Secondly, the grid size is predefined, and as a result, it requires a large number of grids to104

achieve good coverage at a desired resolution.

105

Another solution is the Particle Filter (PF), which is a sequential Monte Carlo method

106that approximates the posterior distribution with a set of weighted particles (Doucet et al.,1072001). As the number of particles goes to infinity, the solution from PF approaches the108

optimal solution. There is a rich literature on PF and its variation (Doucet et al., 2001;

109

Arulampalam et al., 2002; Liu and Chen, 1998). The basic procedure is summarized below

Liu and Yamada

9110

for reference.

111

Sampling. At the beginning of each iteration, the value of each particle is drawn from an

112

important density function q(Θit |Θit−1 , zt ). For i = 1, . . . , NΘit ∼ q(Θit |Θit−1 , zt ).

(5)

113

where ∼ denotes that the sample Θit is drawn according to the distribution q(·).

114

Weight update. PF approximates the posterior with a collection of weighted particles.

P (Θt |z1:t ) ≈NX

wti · δ(Θt − Θit ),

(6)

i=1
115

where wti is the weight for particle i at time t. The sum of total weights are normalized to 1.

NX

wti = 1.

(7)

i=1

116

117

The weights for all particles are updated as new evidence zt comes in and renormalized
at the end of each update.

L(zt |Θit )P (Θit |Θit−1 )iiwt ∝ wt−1,
q(Θit |Θit−1 , zt )
118

(8)

where q(·) is the same important density that appears in the sampling step. To simplify theLiu and Yamada

10

119

calculation, q(·) is often chosen to be the transition prior P (Θit |Θit−1 ). Since the terms cancel

120

out in the right hand side, the new weight is directly proportional to the likelihood L(zt |Θit ).

121

Resampling.

122

suffers sample degeneracy after a few update iterations when the weight is concentrated on

123

a very small number of particles. The decrease in weight variance determines the degree of

124[
degeneracy that can be approximated with N
ef f (Arulampalam et al., 2002),

Because the posterior is approximated with discrete particles, the system

1[N.
ef f = PN
i 2
(w
)
t
i=1

(9)

125

[
Small Nef f indicates severe degeneracy in which case resampling is required. Resampling126

essentially eliminates particles with negligible weight by generating a new set of N equally127

weighted particles according to current distribution P (Θt |z1:t ). There exists many methods

128

for sampling from a discrete distribution, which we will not discuss here.

129

Each iteration typically involves one sampling and one weight update. Resampling only130[
happens when N
ef f drops below a certain threshold.

131

Model

132

In the rest of the section, we discuss the practical implementation details of a PF-based real-

133

time parameter estimation system for multiple earthquakes. The parameters we would like134

to estimate are θ = [x, y, D, M, t0 ], where x = longitude (deg), y = latitude (deg), D = depthLiu and Yamada

11

135

(km), M = JMA magnitude, and t0 = origin time. Complete pseudo code (Algorithm 1) is

136

included in the end of this section.

137

Prior distribution.

138

good prior encodes geographical information such as the location of nearby fault lines to the139

station that first triggered, and the most common magnitudes generated at the fault lines.

140

This information can be compiled from historical earthquake catalog for each station and141

used in real time when initializing the PF. If prior information is absent, then a flat prior142

can be used instead. The choice of prior distribution affect the quality of the estimates and143

the convergence rate. Prior distribution of large coverage may cause the initial estimates to

144

be unstable because little evidence is present. Priors of small coverage may result in slow

145

convergence or false convergence (converging at the wrong values). These tradeoffs can be146evaluated empirically. In this paper, we use a uniform flat prior of ±100 km for location,

147±10 km for depth, ±1 magnitude for event magnitude, and ±10 s for event origin time.

148

Likelihood function.

149

depends largely upon the design of the likelihood function. In addition to the arrival time150

and measured amplitude from the triggered stations that current JMA approach uses, our151

likelihood function also utilize the same information from non-triggered stations as well152

because they also convey important information about the event.

The prior P (θ) determines how the particles are initialized. A

The performance of the particle filter for parameter estimation

Liu and Yamada

12

153

In this paper, we use the attenuation relationship developed by JMA for magnitude

154

estimation. The relationship is stated as follows (Hoshiba and Ozaki, 2013). Let Amax be155

the maximum displacement measured by a seismometer after the onset of an event. The

156

earthquake P-wave and S-wave magnitude Mp and Ms can be expressed as a function of the157

linear distance from the station to the hypocenter (R), the depth of the hypocenter (D), and158

the maximum displacement for P-wave (Apmax ) or the maximum displacement of the entire

159

duration (Ap+smax ).

0.72Mp = log Apmax + 1.2 log R + 5 × 10−4 R − 5.0 × 10−3 D + 0.46,

(10)

−3
−3
0.87Ms = log Ap+s
max + log R + 1.9 × 10 R − 5.0 × 10 D + 0.98.

(11)

160

The relationship between the parameters is illustrated in Figure 1. These formulae are

161

specifically tailored for the geological compositions in Japan (see Data and Resources). The

162

P-wave and S-wave magnitudes are expressed in terms of the maximum displacement Amax

163

rather than the maximum acceleration or velocity because the scatter of displacement is164

smaller.

165

166

Given Equation (10) and Equation (11) and that the displacement is log-normally distributed A ∼ ln N (µ, σ 2 ), we propose the following likelihood function for a single station,

Liu and Yamada

13

2

Aexp )
exp −(log Amax2σ−log
2√
L(z|x, y, D, M, t0 ) =.
Amax · σ 2π

(12)

167

Here Aexp is the expected Amax and σ is the standard deviation of displacement measurement.

168

Depending on whether the station has observed P-wave, S-wave, or neither, the expected

169

maximum displacement and its standard deviation are different. For convenience, by rear-

170

ranging Equation (10) and Equation (11), we can compute Aexp and σ for the following three171

cases.

172

Note that Equation (12) is based on amplitude which departs from standard arrival-time

173

based methods. The main reason for adopting this approach is the observation that the174

information of no shaking is critical in separating and classifying multiple earthquakes that175

occur close in space and time. This will be further discussed in Discussion.

• Has not observed any seismic wave:

log Aexp = log Anoise , σ = σnoise .

(13)• Has observed P-wave:
log Aexp = 0.72Mp − 1.2 log R − 5 × 10−4 R + 5.0 × 10−3 D − 0.46, σ = σp .

(14)

• Has observed S-wave:
log Aexp = 0.87Ms − log R − 1.9 × 10−3 R + 5.0 × 10−3 D − 0.98, σ = σs .

(15)

Liu and Yamada

14

176

Anoise and σnoise are the noise in displacement measurement due to recent environmental

177

noise and can be computed independently for each station by keeping a running window. σp178

and σs can be precomputed from historical earthquake data.

179

The decision of which Aexp to compute for a station depends on whether P-wave, S-wave,

180or neither has arrived at the station. The expected travel time of P-wave and S-wave (tp and181ts ) can be computed with ray theory, given the relative location of the station to a hypocenter

182(x, y, D). Comparison between tp , ts , the absolute current time t, and the absolute event

183

start time t0 gives direct estimation of which Aexp to compute for a station. Figure 2 provides

184

a illustrative summary of these design ideas.
This design of the likelihood function is based on the maximum displacement Amax thata seismometer observes during the shaking of P- or S-wave. However, a seismometer may
not observe the maximum displacement immediately after the wave arrival. In this case,the initial estimates can be highly incorrect using this likelihood function. A simple delay
function g(·) can be included to approximate the instantaneous displacement before themaximum is observed,

Aexp = g(t − t0 − tp )Amax , 0 ≤ g(·) ≤ 1.

(16)

185

where t and t0 are the absolute current time and the absolute event origin time. tp is the186

expected P-wave travel time. An example of g(t) is a left shifted sigmoid function.

Liu and Yamada

15

187

The likelihood L(·|·) is applied in each time step to update the weight of each particle.

188

Assuming that each station makes independent observation and the collection of observations

189

from all stations is z, the complete likelihood function becomesL(z|x, y, D, M, t0 ) =nYL(zi |x, y, D, M, t0 ),

(17)

i=1
190

where n is the number of stations. Note that the independence assumption is a minor

191

simplification since nearby stations may have correlated observations.

192

Generalized Particle Filter

193

Particles are initialized according to a prior distribution on the parameters. Since we are194

approximating an unbounded and continuous 5-dimensional space with a bounded and dis-

195

crete one, care must be taken to ensure that the particles have sufficient coverage and the196

number of required particles stays bounded. This is especially important for the seismic

197

application since both the number of parameters and the range of values they can take are198

large. One way to ensure particle diversity with a limited number of particles is to adopt

199

the Regularized Particle Filter (RPF) approach (Arulampalam et al., 2002).

200

RPF differs from common particle filter only in the resampling stage. Rather than201

sampling from a discrete approximation of the posterior density P (·|z) as in Equation (6),202

RPF samples from a continuous approximation (Musso et al., 2001). More specifically, RPF

Liu and Yamada

203

16

draws samples from the approximation,

P (θ|z) ≈NX

wi · Kh (θ − θi ),

(18)

i=1
204

where Kh (θ) = h1 K(θ/h), h > 0 is the rescaled kernel density of K(·). h is the bandwidth,205

and wi is the normalized weight for particle i. As a comparison, Kh (θ) is the Dirac delta

206

function δ(θ) in the regular particle filter. Special care is given to the design of kernels to207minimize the error between approximated and actual distribution. Under the assumption

208

that all particles are equally-weighted and the density is Gaussian, the optimal kernel is the209

Epanechnikov kernel (Musso et al., 2001).

(
Kopt (x) =nx +2(1 − kxk)
2Cnx

kxk < 10otherwise

,

(19)

210

where nx is the dimension of the parameter space, Cnx is the volume of the unit hypersphere211

in Rnx . Figure 3 lists a few popular kernels in the literature.

212

The bandwidth vector h can be chosen proportionally to the variance in the particle

213

population by computing the Cholesky decomposition of the empirical covariance matrix

214(Bickel and Levina, 2008).

Liu and Yamada

17

215

Approximate Method for Multiple Concurrent Quakes

216

PF allows for solving the Bayesian inference problem when exact inference is intractable;

217

however, for the estimates to approach the optimal solution, the number of required particles

218

must grow exponentially with the number of events.

219

Fortunately, as shown in historical records, the probability of having n concurrent earth-220

quakes within a time window of 60 seconds is exponentially small for large n (n > 3).

221

Incorporating this information into the prior distribution can significantly reduce the size of222

the state space. But the state space may still be too large for efficient real-time computation

223

even with this information. For example, suppose that the quake can be parameterized by224

a 5-parameter vector θ, θ = [x y D M t0 ]T where [x y D]T is the [longitude, latitude, depth]

225

coordinate, M is the event magnitude, and t0 is the event starting time. In the presence of

226

n = 3 quakes, the states to be searched reside in a 5 × 3 = 15-dimensional space.

227

This amount of computation may be executable in reasonable time on a supercomputer228

or a networked system of computers with parallel implementation of particle filter (Durham

229

and Geweke, 2013; Miao et al., 2010). In this paper, however, we propose a simple heuristics

230

to keep track of multiple quakes. The heuristics has the desired property such that the231

complexity grows linearly with the number of the events.

232

As a first approximation, the heurisitcs initializes separate particle filters pf1 (θ1 ), pf2 (θ2 ), . . .

Liu and Yamada

18

233

for all possible quakes rather than keeping track of all events within one particle filter234

pf (Θ = {θ1 , θ2 , . . .}). Each particle filter communicates its current estimate θ̂ at the end of235

each update step to all other particle filters. Specifically, each particle filter pfi computes236

the following posterior at time t,

ˆ , j 6= i}).
P (θit |z, {θjt−1

(20)

237

This approximation breaks down the 5n state space where n is the number of concurrent

238

quakes, and dramatically reduces the required computations to keep all events estimation

239

up to date. It is suboptimal, however, since all the particles from pf1 , pf2 , . . . combined only

240

cover a small fraction of the complete parameter space.

241

The heuristic initializes a new particle filter with each single station P-wave pick, using a242high enough threshold such that noisy detections are filtered out. Since local detection can be243

due to an existing event that is being tracked by another particle filter, it is necessary to con-

244

h
i
dition new initialization on a separate metric. A natural choice of metric is P z|θˆ1 , θˆ2 , . . . ,245

i.e., the probability that the triggered measurement can be explained by existing events.

246

Computation of this metric can follow directly from the single station likelihood calculation247

as in Equation (12); however, determining Aexp is nontrivial in this case since it involves

248

computing the additive effect of the interference of multiple wavefronts. We propose an al-

Liu and Yamada

19

249

ternative metric which allows for rapid computation; the metric is the probability of shaking

250

due to any of the existing events and threshold on the highest probability:(< τ, initialize new pf
max P z|θˆi = max L(z|θˆi ).
ii≥ τ, do nothing
h

251

i

(21)

By tuning the threshold τ , we adjust how conservative the system is in declaring new252

events. The complete algorithm is outlined in Algorithm 1 in the appendix for reference.

253

Results254

We carried out the particle filter parameter estimation approach on the data described in

255

Data and Processing, using a flat prior around the first triggered station and 1,000 particles

256

for each particle filter. The algorithm updates at a one-second interval and all experiments

257

were run in simulated real-time.

258

Case 1: 15 March 2011, 1:36:00 - 1:38:00 (two small earthquakes).

259

were performed during this period of time. Snapshots of the particle distribution for one of260

the runs are shown in Figure 4. The averaged time histories of the estimated parameters

261

across all 20 runs were compared against the JMA unified catalog (marked as dotted lines)

262

in Figure 5. The standard deviations across all runs are included as the error bars. The

263

labeled x-axis corresponds to seconds since the first detection of the first event. As the

20 trials

Liu and Yamada

20

264

results demonstrate, the first particle filter was initiated at the first P-wave arrivals, and265

15 s later, another particle filter was created. This approach successfully identified the two

266

separate events. In addition, all estimates converge within 10 s after the initializations. On

267

average, the method is able to localize the epicenters to within 20 km and produce magnitude268

estimates with an error of ±1, relative to the JMA unified hypocenter catalog (Table 1).

269

Case 2: 20 March 2011, 14:19:00 - 14:21:00 (two small earthquakes). We repeated270

the analyses for the dataset of Case 2, where two small earthquakes occurred 5 s apart. The

271

snapshots of particle distributions and time series of estimated parameters are included in

272

Figure 6 and Figure 7. Note that in this example, because the first event occurred offshore

273

and there were fewer near-source recordings, localization and estimation of other parameters274

are more challenging than for Case 1. Indeed, the results showed that the estimates converge

275

slower (about 30 s for event A), and the averaged localization error was relatively large (about

276

80 km for event A, relative to the JMA unified hypocenter catalog (Table 1). However, the277

algorithm was still able to identify and separate the two events and provide accurate estimates278

of their magnitudes to within ±0.5.

279

Case 3: 11 March 2011, 14:46:00 - 14:49:00 (Tohoku earthquake).

280

dataset of the Tohoku earthquake to show that the approach also works for a single event.

281

The snapshots of particle distributions and time series of estimated parameters are included

We used theLiu and Yamada

21

282

in Figure 8 and Figure 9. Since the event was originated offshore, there was substantial

283

localization error in the initial estimates. However, the averaged error decreased with time284

and converged at less than 40 km at 40 s after the initial P-wave arrival. The magnitude285

estimate grew from 6.0 to 8.4 as the earthquake rupture propagated, which is consistent with286

the earthquake rupture physics. At convergence, all five estimated parameters were close to287

the values in the JMA catalog.

288

Discussion

289

Current JMA methods to detecting and associating multiple quakes perform well when the290

events are far apart in space or time. However, they have been shown to generate many false

291

alarms when events are close in space or time (Sagiya et al., 2011). The empirical studies292

suggest that the particle based heuristic can successfully separate multiple concurrent seismic

293

events and provide reasonable estimates of their parameters. And the speed of convergence294

may be improved by incorporating P-wave arrival time in the likelihood, i.e., the residual

295

between observed and predicted P-wave arrival times. The results show that estimated296

parameters converge in less than 10 s for inland earthquakes. For offshore earthquakes, the297

estimates converge in 20-30 s. In terms of localization error, we observed less than 20 km

298

for inland earthquakes, and 20-80 km for offshore events.

299

In order to classify multiple concurrent earthquakes, the use of non-triggered stations isLiu and Yamada

22

300

important. The current JMA EEW system uses arrival times of waves at only the triggered

301

stations in the hypocenter calculation. As a result, when multiple earthquakes occur around302

the same time and the later event occurs close to the wave arrival times of the earlier

303

event, the EEW system treats these events as one single earthquake. If this is the case304

and the stations around the later event observes non-negligible shakings, the current system305

may overestimate the magnitude because these stations are far away from the estimated306

hypocenter (i.e., the location of the earlier event). In our approach, the likelihood function307

uses information from not only the triggered stations but also the non-triggered ones. This

308

design together with the adaptive measure of Anoise allow the algorithm to identify unaffected309

regions between events and is therefore crucial in separating multiple concurrent earthquakes.

310

Another advantage of our approach is the use of regularized particle filter to circumvent311

the need for intensive computation that traditional grid search requires. Although a prior

312

distribution is still required as mentioned in Model, such a distribution can be compiled from313

historical records. Alternatively, initial measurements can be used to “select” the appropriate

314

priors to achieve better performance (Liu et al., 2011).

315

This approach is also subject to several weaknesses. For example, the algorithm is sen-316

sitive to the choice of prior distribution, the number of particles, the values of Anoise , σnoise ,317

σp and σs . While these values can be adjusted and adapted in real time, it requires extensive

318

empirical studies and analyses of historical records for the algorithm to be robust. Some ofLiu and Yamada

23

319

the slow convergence and high variance results in Results may be attributed to suboptimal

320

choices in these parameters.

321

In this paper we use only three cases to test the proposed method, so we are currently

322

carrying out more extensive evaluations of our method using the many examples of multiple

323

earthquake sequences that have occurred over the last several years.

324

As a side note, the performance of parameter estimation for multiple seismic events is

325

limited by how well one can model the ground motion when multiple wavefronts overlap. In

326

the algorithm proposed in Model, this model is not considered. While the omission makes

327

little difference in the case studies where the events are spatially far apart (greater than 100328km), if we want to apply the same technique to separate aftershocks from mainshock that329

occur close in time, then such model should be considered.

330

Conclusion

331

In the seismically active period, multiple earthquakes of similar or distant origins can take332

place at almost the same time. Failure to identify them as separate events leads to poor

333

estimates of their parameters. The error in estimates can in turn cause false warnings. In

334

this paper, we study the problem of detecting and classifying multiple earthquakes that335

occur close in time. Based on a Bayesian formulation that considers the possibility of having

336

more than one event present at any given time, we propose a novel likelihood function

Liu and Yamada

24

337

suitable for classifying multiple concurrent earthquakes and present a sequential Monte Carlo338

heuristic whose complexity grows linearly with the number of events. The performance339

of the heuristic is empirically validated with three sets of JMA seismic records after the340

2011 Tohoku earthquake. The initial studies show that the approach is able to successfully

341

separate multiple events that occur close in space and time and estimate their parameters342

in realtime to a reasonable degree of precision in comparison to official values determined

343

by JMA in the post event analyses. Although complete validation and characterization are344

required before this method applied in realtime detection, the initial results show that our

345

approach can reduce the chance of overestimation of earthquake magnitude and, as a result,346

contribute to the design of a better EEW system.

347

Data and Resources

348

Waveform data used in the present study were extracted from continuous recordings of349

the stations within the JMA strong motion network.

350

three cases is available at; http://www.seisvol.kishou.go.jp/eq/EEW/kaisetsu/joho/351

20110315013605/content/contentout.html, http://www.seisvol.kishou.go.jp/eq/EEW/

352

kaisetsu/joho/20110320141959/content/contentout.html, and http://www.seisvol.

353

kishou.go.jp/eq/EEW/kaisetsu/joho/20110311144640/content/contentout.html (last354

accessed July 2013). We use Seismic Analysis Code (http://www.iris.washington.edu/

The JMA EEW performance ofLiu and Yamada

25

355

software/sac/manual/fileformat.html, last accessed July 2013) for the data processing.

356

The JMA attenuation relationship are available in the report of the second JMA EEW

357

evaluation committee (http://www.seisvol.kishou.go.jp/eq/EEW/MeetingHYOUKA/t02/

358

shiryou.pdf, last accessed July 2013).

359

Acknowledgements

360

The first author would like to thank professor Jim Mori at the Kyoto University for the361

generous hosting and the warm welcome she received from all the group members. The

362

authors would also like to thank JMA for providing the strong motion seismic data. This

363

work was generously funded by the 2011 National Science Foundation (NSF) East Asian and364

Pacific Summer Institute (EAPSI) Fellowship and the Funding Program for Next Generation

365

World-Leading Researchers (NEXT Program).

366

References

367

Arulampalam, M., Maskell, S., Gordon, N., and Clapp, T. (2002). A tutorial on particle

368

filters for online nonlinear/non-Gaussian Bayesian tracking. IEEE Transactions on Signal

369

Processing, 50(2):174–188.

370

371

Bickel, P. J. and Levina, E. (2008). Regularized estimation of large covariance matrices.
arXiv.org.

Liu and Yamada

372

373

374

375

26

Doucet, A., de Freitas, N., Gordon, N., and Smith, A. (2001). Sequential Monte Carlo
Methods in Practice. Statistics for Engineering and Information Science. Springer.

Durham, G. and Geweke, J. (2013). Adaptive Sequential Posterior Simulators for Massively
Parallel Computing Environments. arXiv.org.

376

Hoshiba, M. and Iwakiri, K. (2011). Initial 30 seconds of the 2011 off the Pacific coast377

of Tohoku Earthquake(M (w) 9. 0)-amplitude and tau(c) for magnitude estimation for

378

Earthquake Early Warning -. Earth Planets Space, 63(7):553–557.

379

Hoshiba, M., Iwakiri, K., and Hayashimoto, N. (2011). Outline of the 2011 off the Pacific380

coast of Tohoku Earthquake (M w 9.0)—Seismicity: foreshocks, mainshock, aftershocks,

381

and induced activity—. Earth Planets Space, 63(7):513–518.

382

Hoshiba, M. and Ozaki, T. (2013). Earthquake early warning and tsunami warning of the383

japan meteorological agency, and theirperformance in the 2011 off the pacific coast of384

tohoku earthquake (mw9.0). In Wenzel, F. and Zschau, J., editors, Early Warning for

385

Geological Disasters. Springer-Verlag Berlin Heidelberg.

386

387

388

JMA (2011). Earthquake early warning report for tohoku 2011. http://www.jma.go.jp/
jma/press/1104/28b/eew_hyouka_2.pdf.

Katsumata, A. (2008). Recursive digital filter with frequency response of a mechanicalLiu and Yamada

27

389

seismograph. Technical report, Seismology and Volcanology Research Department, Mete-390

orological Research Institute.

391

Liu, A. H., Bunn, J., and Chandy, K. M. (2011). Sensor networks for the detection and392

tracking of radiation and other threats in cities. In The 10th International Conference on

393

Information Processing in Sensor Networks (IPSN), pages 1–12.

394

395

Liu, J. S. and Chen, R. (1998). Sequential monte carlo methods for dynamic systems. Journalof the American Statistical Association, 93:1032–1044.

396

Miao, L., Zhang, J. J., Chakrabarti, C., and Papandreou-Suppappola, A. (2010). A new

397

parallel implementation for particle filters and its application to adaptive waveform design.

398

Signal Processing Systems (SIPS), 2010 IEEE Workshop on, pages 19–24.

399

Musso, C., Oudjane, N., and Legland, F. (2001). Improving regularized particle filters. In

400

Doucet, A., de Freitas, N., and Gordon, N., editors, Sequential Monte Carlo Methods in401

Practice. New York, number 12, pages 247–271. Statistics for Engineering and Information

402

Science.

403

404

405

Sagiya, T., Kanamori, H., Yagi, Y., Yamada, M., and Mori, J. (2011). Rebuilding seismology.
Nature, (473):146–148.

Ueno, H., Hatakeyama, S., Aketagawa, T., Funasaki, J., and Hamada, N. (2002). Im-

Liu and Yamada

28

406

provement of hypocenter determination procedures in the japan meteorological agency.

407

Quarterly Journal of Seismology, 65(1-4):123–134.

Liu and Yamada

29

Table 1: Summary of the earthquake information studied in this paper.
Mest Lon
Lat Dep Date Time
M Lon
Lat
Dep
03-15 01:35:57.35 2.5 138.610 36.938 3.4
Case 1 5.9
138.6 36.9 10
03-15 01:36:12.72 3.3 139.879 35.526 20.5
03-20 14:19:38.27 3.0 141.935 38.286 42.3
Case 2 7.6
142.1 38.2 30
03-20 14:19:58.06 4.7 140.794 37.082 7.2Case 3 8.6
142.7 38.2 10
03-11 14:46:48.08 9.0 142.861 38.103 23.7
The first four columns correspond to the real-time JMA EEW records. The last six columns
are the values documented in the JMA unified hypocenter catalog. Both Case 1 and 2
contain two events.
(a)Station(b)

Ap+smax

Apmax

D

R

tp
ts

Hypocentert0

P-waveS-waveFigure 1: Illustrations of the parameters used in Model. (a) Hypocenter and seismic stationand (b) amplitude and arrival times. tp and ts mark the arrival time of the P-wave andS-wave since the start of the earthquake at t0 , tp ≤ ts .

Liu and Yamada

30

EpicenterP-wave front

S-wave front

Stations thathave bothP- and Swave
Eq. (15)

Stations thathave observed
only the Pwave
Eq. (14)

Stations thatseismic shaking
Eq. (13)

Figure 2: Illustrative summary of the design of a single station likelihood function. The
expected observation made by a station depends on whether it should have observed Pwave, S-wave, or neither, given an hypocenter estimate.

1

Epanechnikov
Uniform
Trianglar

0.8
0.6
0.4
0.2
0
−1

−0.5

0

0.5

1

Figure 3: Some popular smoothing kernels used in regularized particle filter. Each kernel
integrates to 1 to ensure that the resulting density is still a probability density function.

Liu and Yamada

38

31

(a)

38

37.5

37.5

37

37

36.5

36.5

36

36

35.5

35.5

35

35

34.5

34.5

34136
38

137

138

139

140

141

(c)

(b)

34136

137

138

139

140

141

(d)

37.5
37

A

36.5
36

B

35.5
35
34.5
34136

137

138

139

140

141

Figure 4: Distributions of 2,000 Particles visualized on the map at (a)1 s, (b)2 s, (c)14 s,
and (d)17 s after 1:36:07 on 15 March, 2011. The time correspond to seconds elapsed since
the first P-wave detection. The official epicenters for the two events as appeared in the JMA
catalog are marked as stars and labeled in (d) for reference.

Liu and Yamada

(a)

25

Magnitude estimate

Event AEvent B

30

20
15
10
5
0
0

25

depth estimate (km)

5

5

10

(c)

15

Time (s)20

25

Event AEvent B

4
3
2
1

5

10

(d)

15

20

25

30

15

20

25

30

Time (s)15
Event AEvent B

20
15
10
5
0
0

(b)

0
0

30

t0 estimate (sec)

Localization error (km)

35

32

5

10

15

Time (s)

20

25

30

Event AEvent B

10
5
0
−5
−10
0

5

10

Time (s)Figure 5: Results compiled from 20 independent runs for the period between 1:36:07 and1:36:37 on 15 March, 2011. Time histories of the (a) localization error, (b) magnitude,(c) depth of the hypocenter and (d) origin time of the event. The two events are labeledaccording to Figure 4(d). Averaged time histories across all 20 runs are marked as solidlines, and the official values that appear in the JMA catalog are marked as dashed lines.
The standard deviations across all runs are shown as error bars. The time displayed on thex-axis is relative to the first pick from the earliest event.

Liu and Yamada

41

33

(a)

41

40

40

39

39

38

38

37

37

36

36

35

35

34
139
41

139.5 140(c)

140.5 141141.5 142142.5 143

34
139
41

40

40

39

39

38

38

37

37

36

36

35

35

34
139

139.5 140

140.5 141141.5 142142.5 143(b)

139.5 140

(d)

34
139

140.5 141B

139.5 140

140.5 141141.5 142142.5 143A

141.5 142142.5 143

Figure 6: Distributions of 2,000 Particles visualized on the map at (a)2 s, (b)7 s, (c)17 s,
and (d)37 s after 14:19:56 on 20 March, 2011. The symbols are defined the same way as in
Figure 4.

Liu and Yamada

34(a)Magnitude estimate

Event AEvent B

100
80
60
40
20
0
0

50

depth estimate (km)

5

10

(c)

40

20

Time (s)

30

4
3
2
1

Event AEvent B
10

(d)

20

Time (s)

30

40

60
Event AEvent B

30
20
10
0
0

(b)

0
0

40

t0 estimate (sec)

Localization error (km)

120

10

20

Time (s)

30

40

55

50
Event AEvent B
45
0

10

20

Time (s)30

40

Figure 7: Results compiled from 20 independent runs in the period between 14:19:56 and14:20:36 on 20 March, 2011. The subfigures and included symbols are defined the same way
as in Figure 5.

Liu and Yamada

40

35

(a)

40

39.5

39.5

39

39

38.5

38.5

38

38

37.5

37.5

37
141
40

(c)

141.5

142

142.5

143

143.5

144

37
141
40

39.5

39.5

39

39

38.5

38.5

38

38

37.5

37.5

37
141
40

141.5

142

142.5

143

143.5

144

(e)

40
39.5

39

39

38.5

38.5

38

38

37.5

37.5

141.5

142

142.5

143

143.5

144

(d)

37
141

39.5

37
141

(b)

141.5

142

142.5

143

143.5

144

141.5

142

142.5

143

143.5

144

141.5

142

142.5

143

143.5

144

(f)

37
141

Figure 8: Distributions of 1,000 Particles visualized on the map at (a)2 s, (b)7 s, (c)13 s,(d)22 s, (e)32 s, and (f)62 s after 14:46:46 on 11 March, 2011. The symbols are defined thesame way as in Figure 4.

Liu and Yamada

36(a)Magnitude estimate

120
100
80
60
40
20
0
0

35

depth estimate (km)

9

20

40

(c)

60

Time (s)

80

100

8
7
6
5

20

40

20

40

(d)

60

80

100

120

60

80

100

120

Time (s)

35

30

25

20

15
0

(b)

4
0

120

t0 estimate (sec)

Localization error (km)

140

20

40

60

Time (s)

80

100

120

30
25
20
15
10
0

Time (s)Figure 9: Results compiled from 15 independent runs for the period between 14:46:46 and14:48:46 on 11 March, 2011. The subfigures and included symbols are defined the same way
as in Figure 5.

408

AAppendix

Liu and Yamada

37

Algorithm 1: Outline of regularized Particle Filter for multiple seismic event detection.
The ”CONVERGED” criteria can be substituted with desired conditions, e.g. change
in estimates θ̂t−10 − θ̂t−1 < δ.
P F ← {}
Initialize thresholds τ , αInitialize bandwidth vector h0 ∈ Rnx
while not end do? Check for new eventZ ← list of station measurements that triggeredfor z ∈ Z dopr ← maxk L(z|θk )
if pr



 < τ then
N
,
z
←
RP
F
{θ
,
w
}
{θi , wi }N
i
i
i=1
i=1
for i = 1 → N do
Draw θi ∼ P (θ, z)Assign weights based on prior and z, wi ∼ P (θ, z),z
pf ← {θi , wi }N
i=1
P F ← P F ∪ pf
? Update weight, resample if needed
for pf ∈ P F do
{θi , wi }N
i=1 ← pf
for i = 1 → N do
wi ← wi L(z|θi )




N
{θi , wi }N
i=1 ← N ORM ALIZE {θi , wi }i=1
1[
Compute N
ef f ← PN w2
i=1

i

d
if N
ef f < α Nthen{θi , wi }i=1 ← RESAM P LE {θi , wi }N
i=1 , zfor i = 1 → N do
Draw  ∼ K from the Epanechnikov Kernel
Compute weighted empirical covariance matrix Sk of {θi , wi }N
i=1
Compute lower triangle Dk = chol(Sk ), Dk DkT = Skθi ← θi + h0 Dk 
? Check for termination
for pf ∈ P F do
{θi , wi }N
i=1 ← pf
if CON V ERGED {θi , wi }N
i=1 then
P F = P F − pf