www.advances.sciencemag.org/cgi/content/full/1/3/e1
500036/DC1
Supplementary Materials for
Crowdsourced earthquake early warning
Sarah E. Minson, Benjamin A. Brooks, Craig L. Glenn
ie, Jessica R. Murray, John O. Langbein,
Susan E. Owen, Thomas H. Heaton, Robert A. Iannucci
, Darren L. Hauser
Published 10 April 2015,
Sci. Adv.
1
, e1500036 (2015)
DOI: 10.1126/sciadv.1500036
This PDF file includes:
Text
Fig. S1. Background position noise for various GNSS
receivers found on
consumer devices.
Fig. S2. Spectra of drift of position time series f
or various GNSS receivers found
on consumer devices.
Fig. S3. Observed time series from consumer acceler
ometers and GNSS receivers.
Fig. S4. Hayward fault earthquake scenario.
Fig. S5. Epicenter location uncertainty for Hayward
fault scenario rupture.
Fig. S6. Tohoku-oki earthquake example.
Fig. S7. Epicenter location uncertainty for Tohoku-
oki earthquake.
Table S1. Description of observed GPS earthquake di
splacement time series
shown in Fig. 2B.
Table S2. Number of data used and detection respons
e times for Hayward fault
simulation.
Table S3. Locations of GEONET GPS stations used in
analysis of Tohoku-oki
earthquake.
References (
31
–
34
)
Crowdsourced Earthquake Early Warning
Supplemental Material
1. Operational Early Earthquake Warning
EEW is the subject of active research in many parts
of the world. There are few regions,
however, where EEW systems are currently operating
and sending information to users other than the
researchers developing the warning system (Fig. 1).
EEW operates throughout Japan and Taiwan. In
Mexico, the cities of Mexico City, Oaxaca City, Tol
uca, Chilpancingo, Acapulco, and Morelia currently
receive warnings (Juan Manual Espinosa Aranda, Shri
Vishna Singh, Carlos Valdes, personal
communication, 2014). Warnings are sent to the Mar
mara region of Turkey including Istanbul province
(Erdik Mustafa, personal communication, 2014). Fou
r nuclear power plants in Switzerland currently
receive warnings (John Clinton, personal communicat
ion, 2014) as do assorted nuclear facilities, dams,
and public and private emergency responders in Roma
nia (Constantin Ionescu, personal communication,
2014). An EEW system has recently begun operating
in the Chengdu region of China (Tun Wang,
personal communication, 2014.) Finally, while an E
EW system is currently in beta testing for most of
the west coast of the United States including Calif
ornia, Oregon, and Washington, at present only
California has external users of the U.S. EEW syste
m. All of these EEW systems utilize seismic data
only, although GNSS data are beginning to be incorp
orated into the U.S. and Japanese EEW systems
(Yusaku Ohta, personal communication, 2014) and the
Mexico system will soon use GNSS data as well
(Juan Manual Espinosa Aranda, Shri Vishna Singh, Ca
rlos Valdes, personal communication, 2014).
Background color in Fig. 1 is peak ground accelerat
ion with 10% probability of exceedance in 50 years
given a 475 year return period from the Global Seis
mic Hazard Assessment Program. Most regions of
the world with high seismic hazard have no EEW capa
bilities (Fig.1).
2. Empirical Noise Characterization
We assess a typical smartphone’s capability to dete
ct permanent surface displacements by first
empirically characterizing positioning noise from a
variety of different consumer GNSS receivers
(including those found in Google Nexus 5 smartphone
s, first-generation Apple iPad tablets, u-blox
receivers, and Garmin eTrex handheld receivers) as
well as Japan’s GEONET GNSS receivers. All
GNSS data analyzed are coarse acquisition (C/A) cod
e (
19
). For the time periods most pertinent to
EEW application (seconds to tens of seconds), we fi
nd that noise character depends both on the device
model version and whether it is stationary or movin
g (Fig. S1-S2). When stationary, epoch-to-epoch
repeatability (the ability to recover position at e
ach epoch) for modern devices (e.g., Nexus 5) is ~0.
01-
0.02 m and for older devices it is less than ~1 m (e
.g., iPad) (Fig. S1). When moving at a constant
velocity, the Nexus 5 repeatability is ~1.5 m (Fig.
S1). This strongly suggests that, when stationary,
the
Kalman filter onboard the smartphone assumes an ess
entially constant velocity model – an assumption
that would hold for most pedestrian navigation appl
ications where it would be desirable to minimize
accelerations due to bodily motion or device manipu
lation.
The repeatability for all of these devices depends
upon the time-span of the measurements
because they all have background noise that exhibit
s temporal correlations (i.e., the noise spectra ar
e
“red”). Consequently, the time series of position
estimates exhibit drift that we assess following (
31
)
(Fig. S2). At intervals of a few seconds, the raw
C/A code positions can resolve displacements of 10
cm; but for longer intervals (e.g., 60 seconds), th
eir sensitivity decreases to 2-3 m. (Note that the
re is
less drift in positions estimated from C/A code dat
a recorded at GEONET GPS stations than from C/A
code data recorded on a consumer receiver. This is
because scientific GPS stations have higher qualit
y
hardware. The most significant difference is proba
bly the reduction in multipathing made possible by
the choke ring antennas at scientific GPS sites.)
Satellite Based Augmentation System (SBAS)
corrections or phase smoothing significantly improv
es the precision of the consumer devices by a facto
r
or 10 or more.
3. Offset Tests
As discussed in the main text, we subjected the pla
tform to a series of displacements ranging
from ~0.1-2.0 m (Fig. 2, S3
)
and find that the raw C/A positions (with SBAS) fa
ithfully records the
entire displacement range. In contrast, the smartp
hone is essentially insensitive to the displacement
s
because its onboard data processing uses a Kalman f
ilter that removes most of the external motion to
which the device was subjected.
To assess the full potential of consumer grade GNSS
and INS devices for earthquake monitoring,
we constructed our own extended Kalman filter using
raw C/A code data from the u-blox receiver and
accelerometer data from the Nexus 5. The Kalman fi
lter uses a constant jerk state space model for
system dynamics (
21
). Given GNSS positions and three-component accele
rometer observations, the
filter simultaneously estimates displacement and ac
celerometer bias. Fig. 2A in the main text shows t
he
results of a 50 cm offset test and demonstrates tha
t the Kalman filtered solution could accurately rec
over
the displacement time history to which the receiver
was subjected even though the twice-integrated
(yielding displacement) accelerometer record demons
trates a large amount of drift over earthquake time
-
scales.
4. Earthquake Monitoring Capabilities of Consumer D
evices
Fig. 2B compares the empirical drifts we estimated
for a variety of instruments (Fig. S2) with
observed displacement time series for recent notabl
e damaging earthquakes (Table S1) (
22, 23, 25
). The
strong motion accelerometer drift was estimated by
differencing a twice-integrated observation of the
2011 Tohoku-oki earthquake recorded on the east com
ponent of K-net strong motion station MYG011
and a corresponding displacement time series from n
earby GEONET GPS station 0550. The consumer
accelerometer drifts were obtained by averaging at
62,225 observation locations the absolute value of
the difference between ground motion time series fr
om the
M
w
7 Hayward fault scenario and twice-
integrated corresponding simulated acceleration tim
e series using noise from a Nexus 5 smartphone’s
onboard accelerometer. Similarly, the Kalman filte
r drift time series were obtained by averaging at
those same observation locations the absolute value
of the difference between the ground motion time
series and corresponding simulated Kalman filter ti
me series computed by combining the same
acceleration records with simulated C/A code GNSS t
ime series in the same manner as Section 4. These
simulated time series were originally 90 sec durati
on and then extrapolated from 90 sec to 100 sec. B
oth
the simulated and empirical noise time series were
smoothed for plotting purposes. Note that noise
characteristics for data involving accelerometer ob
servations are very sensitive to the type of motion
to
which the instrument was subjected. This can be cl
early seen for the Kalman filter drift curve which
was calculated for an
M
w
7 earthquake. Once the ground shaking ceases, the
Kalman filter acts to damp
further drift of the estimated displacement.
Figure 2C shows the minimum magnitude observable (a
s a function of distance) with a signal-to-
noise ratio (SNR) of at least 10 for the devices in
Fig. 2B. The SNR for displacement observations wa
s
computed by assuming that the signal is the peak gr
ound displacement (PGD) as a function of
magnitude and distance (
27
), which we assume occurs at approximately the S-wa
ve arrival time. Thus
the noise is the amount of drift that accumulates b
etween the P-wave arrival and S-wave arrival
assuming wave speeds of 6 km/s and 3.5 km/s, respec
tively, where the drift accumulation follows the
curves shown in Fig. 2B. The SNR curve for the acc
eleration observed by a consumer accelerometer
was computed by assuming the amplitude of the signa
l was given by the peak ground acceleration
(PGA) (
27
), and that the noise over short time durations was
Gaussian with an amplitude of 0.031 m/s
2
based on the noise characteristics we observed over
short time periods for the accelerometer in a
Nexus 5 smartphone. As in Fig. 2B, the curves in F
ig. 2C were smoothed for plotting purposes.
5. Hayward Fault Rupture Simulation
We simulated a crowdsourced EEW response to a scena
rio magnitude 7 earthquake on
California’s Hayward fault using two different type
s of data. In the main text (Fig. 3), we presented
results from assuming that different fractions of t
he population contributed data from consumer device
s
containing both an accelerometer and a simple GNSS
C/A code receiver equipped with neither SBAS
differential corrections nor phase-smoothed L1 data
. (Almost all smartphones record acceleration and
C/A code data.) We then Kalman filtered the twice-
integrated acceleration time series and GNSS
position time series, and modeled the resulting dis
placement estimates. In Fig. S4, we show the resul
ts
from assuming different fractions of the population
contributed data from consumer devices with SBAS-
enabled GNSS receivers with phase smoothing but no
accelerometer. (This is the highest quality type o
f
GNSS receiver commonly found in consumer devices.)
The results in Fig. 3 and Fig. S4 are nearly
identical. Both types of positioning data yield ea
rthquake detection in 5 seconds (assuming 0.2% of t
he
population submits usable data), earthquake locatio
ns with errors <5 km and very accurate moment
magnitude estimates at all times for all levels of
participation.
Data Generation
The uncertainties of GNSS vertical positions are ty
pically much worse (at least a factor of two)
than horizontal positions (
19
). Accordingly, all of the modeling in this paper
was accomplished using
only the 2-D horizontal displacements calculated fr
om the east and north position time series.
To simulate the horizontal displacements we would o
btain from consumer devices located near
the rupture of a moderate magnitude earthquake, we
began with the predicted horizontal ground motion
time series for the
M
w
7 Hayward fault scenario rupture called hs+hn_n04_
hypoO_vr92_tr15 (
32
).
Smartphone GNSS (raw C/A code) and acceleration rec
ords were then simulated by adding to these time
series the empirical noise spectra of displacements
obtained from a stationary raw C/A code GNSS
receiver (Garmin eTrex) (
31
) and the accelerometer on a Nexus 5 smartphone. W
e then applied a
Kalman filter to combine the resulting C/A code and
twice-integrated acceleration time series. The
choice of the Garmin eTrex is a conservative one as
it is the noisiest of all the devices (Fig. S1-S2)
. We
also simulated consumer GNSS data with real-time co
rrections by adding to the original synthetic
displacement time series the noise characteristics
we had determined empirically for the stationary u-
blox receiver with SBAS and phase smoothing. These
two simulated datasets allow us to explore the
results obtainable from two hardware configurations
we are likely to encounter in a crowdsourced
earthquake monitoring system: one that provides GNS
S and accelerometer data without SBAS or phase
smoothing, and one that provides GNSS data with dif
ferential corrections from SBAS.
In the Kalman filter case, we computed east and nor
th component displacement time series with
a one second sampling rate from independent Kalman
filters using a constant jerk state space model for
system dynamics (
21
). Given a position time series and an acceleromet
er time series, the filter
simultaneously estimates displacement and the bias
for the accelerometer measurements. This is a
dynamic model that accounts for the changes in posi
tion, velocity, and acceleration that these devices
will be subjected to during an earthquake. Simulat
ed earthquake displacements only give three degrees
of freedom (x, y, and z displacement). For a full
inertial system simulation 6 degrees of freedom (3
displacements and 3 rotations) are required. We ma
de the simplifying assumption that the sensitive ax
es
of the accelerometer triad were always aligned with
the simulated motion directions of the earthquake.
With this assumption, because the accelerometers ar
e mutually orthogonal (as are the simulated
displacements), we can treat the x, y, and z displa
cements separately as three 1-D estimation problems
.
Observation locations were chosen to mimic the dist
ribution of data we might obtain if a certain
percentage of the population participated in crowds
ourced earthquake monitoring. Our approach was to
choose a random distribution of land-based observin
g locations within a 35 km-by-100 km rectangular
box centered on the Hayward fault with the number o
f observations equal to a percentage of the
population within that box (Fig. 3, S4). Since onl
y a small fraction of the population might choose t
o
participate in a crowdsourced EEW system and only a
small fraction of those participants may be
contributing usable data at any given time, we purp
osefully chose to consider cases where a very small
percentage of the population reported observations
of the earthquake. Specifically, we considered the
cases where we obtained data from 0.0125%, 0.025%,
0.05%, 0.1%, and 0.2% of the population in the
35 km-by-100 km region. These cases correspond to
294, 587, 1174, 2348, and 4696 observations,
respectively (Table S2).
Earthquake Detection
At each second, we counted the number of observing
devices that detected a trigger where we
determined that a device had been triggered if it a
nd the four devices nearest to it each have accumul
ated
a horizontal displacement in excess of 5 cm. This
criterion requires that the displacement signal be
spatially coherent as we expect earthquake displace
ments to be. We then declare that an earthquake ha
s
been detected when at least 100 observers simultane
ously report triggers. This criterion is reached i
n
less time with increasing numbers of observations s
imply because, with a denser set of observations, t
he
seismic waves do not have to travel as far to reach
100 observers (Table S2).
Locating the Epicenter
At each second after detection, we obtain an earthq
uake epicenter by fitting a power law to the
amplitudes of the horizontal offsets reported by al
l triggered observers and searching for the source
location that minimizes the L1 norm of the residual
s (Fig. S5). This estimation is not computationall
y
expensive and so the time required to obtain an epi
center depends on the time required for the
displacements to propagate from the earthquake rupt
ure to the receivers. Our errors in the epicenter
location are always less than 5 km (Figs. 3, S4).
Note that we use a power law regression of unknown
power: log
A
=
c
0
+
c
1
log
r
, where
A
is the
amplitude of the observed horizontal displacement,
r
is the distance between the observer and the
proposed epicenter, and
c
0
and
c
1
are the unknown linear regression parameters. Alt
ernatively, we
might have fixed
c
1
=−2
because surface offsets exhibit an
r
-square decay in the far-field (
5
). However,
we are using near-field observations for which the
r
-square decay does not necessarily hold. Therefore
,
we solve for the power of the decay in the linear r
egression.
Determining Earthquake Magnitude
We determine the magnitude at each second by first
calculating a finite fault slip model and then
computing the moment magnitude associated with that
slip model. The finite fault slip models are
computed using the Bayesian linear regression appro
ach of (
25
). In that paper, the authors solved for
the fault geometry and distributed slip model simul
taneously using three-component real-time GPS data.
However, here we only have the two horizontal compo
nents of motion. Our initial tests indicated that,
without the vertical component of deformation, we c
ould not constrain the fault orientation. Instead,
we
used the accurate epicenter locations obtained to j
ustify assuming that the rupture was located on the
Hayward fault, and then assumed a planar fault whos
e strike and dip represented the average strike and
dip of the Hayward fault in the region of the epice
nter (strike 325.0000°, dip 76.2309°). (For global
application, existing maps of the largest faults wo
uld need to be employed.) We sub-divided the fault
plane into one row of patches in the down-dip direc
tion and eight patches in the along-strike directio
n,
each patch having a length of 10 km on a side. Our
Green’s functions were computed for rectangular
dislocations in a homogeneous elastic half-space.
The moment magnitudes derived from the results of
our slip inversion at each second show high accurac
y and low latency (relative to the scenario
earthquake’s actual moment release as a function of
time) for all data sets tested (Figs. 3, S4).
6. Tohoku-oki Earthquake Model
Data
As part of routine collection of dual frequency, sc
ientific-grade GPS carrier phase data, the
GEONET GNSS stations also record the same pseudoran
ge data that are used for positioning in
consumer devices (
19
). To simulate real-time positioning using single
frequency data on a consumer
device, we estimated positions based only on C/A co
de pseudorange values on the L1 frequency (“C1”
observations) from a selection of GEONET GPS statio
ns (Table S3) using the pr2p (pseudorange to
position) algorithm developed by Dr. Mark Miller at
JPL as part of the GIPSY-OASIS software package
(
33
) (Fig. S6). The algorithm calculates the position
and receiver clock bias at each epoch by
minimizing a penalty function, the RMS of residuals
, for that epoch. An outlier detection algorithm i
s
triggered when the penalty function or individual r
esiduals do not satisfy solution tolerance criteria
. The
residual threshold for individual pseudorange obser
vations was 40 m, and the RMS pseudorange
threshold was 20 m. Once outliers are identified a
nd deleted, the penalty function is minimized again
.
Broadcast satellite orbits and satellite clock corr
ections were used. There are no troposphere or
ionosphere models applied. The estimation uses the
prior epoch’s position and receiver clock bias as
an
initial value for the next epoch’s position and clo
ck bias value, but there is no temporal smoothing
between position estimates.
Fig. S6 compares the displacement time series from
the pr2p processing of the GEONET data for
the Tohoku time series to displacement time series
obtained through scientific-quality processing that
takes advantage of the dual-frequency data these sc
ientific GPS stations collect. At displacements le
ss
than 0.5 m, the comparison is degraded by a subset
of far-field stations that experienced a spatially
confined ionospheric disturbance combined with loss
of satellite lock. However, at the vast majority
of
stations, the offsets estimated from the C/A code a
nd post-processed scientific-grade time series are
nearly equal (Fig. S6C).
Earthquake Detection
We use the same method for detecting the Tohoku-oki
earthquake as we did for the Hayward
fault simulation. However, because the estimated G
PS positions include ionospheric noise and spurious
jumps due to loss of satellite lock, we increase th
e minimum horizontal displacement required for
triggering to 75 cm. Further, rather than declarin
g that an earthquake has been detected when a certa
in
number of triggers have been received, we instead s
ay that we have detected an earthquake when the
number of triggers exceeds 5
σ
relative to the triggering rate during the pre-eve
nt period. This accounts
for the fact that, at any given time, some subset o
f GPS stations exhibit spurious displacements. Thi
s
threshold is reached at 77 sec after the origin tim
e.
Locating the Epicenter
We used the same approach to calculate the earthqua
ke epicenter as we used for the Hayward
fault simulation (Fig. S7).
Determining Earthquake Magnitude
At each second, we estimated the total magnitude re
lease for the Tohoku-oki earthquake by
estimating a finite fault slip model and calculatin
g the moment magnitude associated with that slip
model. We used the same inversion methodology to i
nvert for the slip distribution as we did for the
Hayward fault simulation combined with the 1-D laye
red elastic structure and 3-D curved fault
geometry of (
34
). The resulting magnitude release as a function o
f time closely matches the source-time
function obtained by a full kinematic inversion of
the Tohoku-oki earthquake rupture using scientific-
grade kinematic GPS, static GPS offsets, seafloor g
eodesy, and both near-field and far-field tsunami
data (
24
).
Tsunami Arrival Times
The tsunami arrival times in Fig. 4 were estimated
from the time of the first increase in sea level
recorded by the ocean-bottom pressure gauges, GPS w
ave gauges, and coastal tide and wave gauges
included in (
28
).
7. Determining Pre-Earthquake Position
Detecting, locating, and modeling the slip distribu
tion of an earthquake requires estimating the
observed offset at each second of time relative to
the sensor's position prior to the beginning of the
coseismic deformation. Note that the later the tim
e of the prior position, the less long-period drift
may
accumulate between the prior position and the senso
r's position during the earthquake. In fact, this
prior
position need not be from a time before the origin
time of the earthquake, only before the first
displacement at the sensor's location.
We leave for a future study determining the best wa
y to identify the pre-earthquake position. In this
paper, we presented results for earthquake detectio
n from two end-member cases. For the Tohoku-oki
earthquake, we used a pre-event position that was e
stimated from one minute of data observed 1,000
seconds before the current epoch. In contrast, for
the Hayward scenario, we used as our prior positio
ns
the instantaneous position of each observer at the
origin time.
Figure S1
.
Background position noise for various GNSS receiver
s found on consumer devices
.
Although each receiver is stationary, there is alwa
ys some uncertainty in its computed position. The
magnitude of the apparent horizontal displacement i
s shown for position time series computed using
C/A code pseudorange-based positioning for six scie
ntific-grade GEONET GPS stations in Japan, a
first-generation iPad, a Nexus 5 smartphone, and u-
blox and Garmin eTrex consumer GNSS receivers.
The position time series for the GEONET stations we
re computed using the pr2p program, part of the
GIPSY-OASIS software package. All other positions
were obtained directly from each receiver
according to each one’s onboard proprietary process
ing software. Note that the GNSS receiver in the
Nexus 5 smartphone applies a Kalman filter to the o
utput position time series and that the parameters
of
the filter differ depending on whether the receiver
appears to be stationary or in motion. The output
position time series for both the stationary and mo
ving cases are shown. (For the moving receiver,
positions are relative to actual position of the re
ceiver determined from scientific-grade GNSS.) Als
o,
note that the u-blox time series is phase-smoothed
and augmented with satellite-based differential
corrections (SBAS) that significantly reduce the no
ise level of its position estimates relative to the
other
receivers and processing schemes (with the exceptio
n of the heavily-filtered positions output by the
Nexus 5 in stationary mode).
Figure S2. Spectra of drift of position time series
for various GNSS receivers found on consumer
devices.
The drift can be characterized by computing the R
MS of (x(t+
τ
)-x(t)) as a function of various
intervals,
τ
, for position time series obtained as in Fig. S1.
Initially, the background noise is quantified
by computing a power spectral density (PSD) and fur
ther refined by characterizing the PSD as a
combination of a power-law and white noise processe
s. Using the parameters obtained for the power-
law and white noise, time series data are simulated
and the drift is evaluated to obtain a confidence
interval. For each set, the solid line is the nomi
nal RMS of the drift while the dashed lines represe
nt the
95% confidence interval determined from 200 simulat
ions.