NET-VISA: Network Processing Vertically Integrated Seismic Analysis

21 Pages·2013·1.78 MB·English

BulletinoftheSeismologicalSocietyofAmerica,Vol.103,No.2a,pp.–, ,doi:10.1785/0120120107
NET-VISA: Network Processing Vertically Integrated Seismic Analysis
by Nimar S. Arora, Stuart Russell,* and Erik Sudderth
Abstract Theautomatedprocessingofmultipleseismicsignalstodetectandlocal-
ize seismic events is a central tool in both geophysics and nuclear treaty verification.
Thispaperreportsonaprojectbegunin2009toreformulatethisprobleminaBayesian
framework. A Bayesian seismic monitoring system, NET-VISA, has been built com-
prisingaspatialeventpriorandgenerativemodelsofeventtransmissionanddetection,
aswellasaninferencealgorithm.Theprobabilisticmodelallowsforseamlessintegration
of various disparate sourcesofinformation. Appliedinthe context of the International
Monitoring System (IMS), a global sensor network developed for the Comprehensive
Nuclear Test Ban Treaty (CTBT), NET-VISA achieves a reduction of around 60% in
thenumberofmissedeventscomparedwiththecurrentlydeployedsystem.Italsofinds
events that are missed by the human analysts who postprocess the IMS output.
Introduction
The Comprehensive Nuclear Test Ban Treaty (CTBT), detections (as much as 90% of all detections are false); and
which bans all nuclear explosions on Earth, is gaining re- (3)manythousandsofdetectionsarerecordedperday,sothe
newed attention in light of growing worldwide interest in problem of proposing and comparing possible events (sub-
mitigating the risks of nuclear weapons proliferation and setsofdetections)isdaunting.Theseconsiderationssuggest
testing. To monitor compliancewith the treaty, the Prepara- that an approach based on probabilistic inference and com-
tory Commission for the Comprehensive Nuclear Test Ban binationofevidencemightbeeffective,andthispaperdem-
TreatyOrganization(CTBTO)hasinstalledasuiteofsensors onstrates that this is in fact the case. For example, such an
known as the International Monitoring System (IMS). The approach automatically takes into account nondetections as
IMS includes waveform sensor stations (seismic, hydroa- negative evidence for a hypothesized event, something that
classical methods typically do not do.
coustic,andinfrasound)connectedbyaworldwidecommu-
Theexistingnetworkprocessingalgorithminuseatthe
nications network to a centralized processing system in the
IDC,knownasglobalassociation(GA)(LeBrasetal.,1994),
InternationalDataCenter(IDC)inVienna.TheIDCoperates
uses various heuristics to cluster arrivals, and then deter-
continuouslyandinrealtime,performingstationprocessing
mines the location of the event in each cluster. The event
(analysis and reduction of raw seismic sensor data to detect
location algorithm is based on the original iterative linear
and classify signal arrivals at each station) and network
leastsquaresmethodofGeiger(1910,1912).Ofcourse,the
processing(associationofsignalsfromdifferentstationsthat
algorithmhasbeenenhancedconsiderablyovertheyears,for
havepresumablycomefromthesameevent).Perfectperfor-
example,bytheuseofsingularvaluedecompositiontosolve
manceremainswellbeyondthereachofcurrenttechnology:
the resulting matrix equations (Menke, 1989; Lay and Wal-
the IDC’s automated system, a highly complex and well-
lace, 1995) and by the use of azimuth and slowness to con-
tuned piece of software, misses nearly one-third of all seis-
strainthesolution(Magotraetal.,1987;Robertsetal.,1989).
miceventsinthemagnituderangeofinterest,andabouthalf However,inthewordsofMyersetal.(2007)(p.1049),“Seis-
of the reported events are spurious. A large team of expert miceventlocationis—atitscore—aminimizationofthedif-
analysts postprocesses the automatic bulletins to improve ference between observed and predicted arrival times.”
their accuracy to acceptable levels.
Further,alloftheseclassicalmethodsrelyonlyontheasso-
The IDC results indicate that the network processing ciated arrivals to locate an event. Even the multiple-event
problem is far from trivial. There are three primary sources
location algorithms, such as those due to Waldhauser and
of difficulty: (1) the travel time between any two points on Ellsworth (2000) and Myers et al. (2007), ignore data from
Earth and the attenuation of various frequencies and wave stations that fail to detect an event.
types are not known accurately; (2) each detector is subject Inthispaper,wepresentaprobabilisticmodelofseismic
to local noise that may mask true signals and cause false eventoccurrence,Pθ (seeTable1foralistofallmathemati-
cal notations), and another probabilistic model of seismic
*Also at Bayesian Logic, Inc., 32432 Lighthouse Way, Union City, detections triggered by events (or noise), Pϕ. We also de-
California94587. scribehowthemodelparametersθandϕareestimatedfrom
1
2 N. S. Arora, S. Russell, and E. Sudderth
(a)
Table 1
Mathematical Notations
Symbol Description
Pθ Priorprobabilitydistributionofevents
Pθ;n Priorprobabilityofthenumberofevents
Pθ;t,Pθ;d,Pθ;l,Pθ;m Priorprobabilityofeventtime,depth,
location,andmagnitude
e Setofevents;ei istheithevent (b)
ei,ei,ei,ei Eventtime,depth,location,and
t d l m
magnitude
Pϕ Probabilityofarrivalsgivenanevent
Pϕ;d Detectionprobabilityofanevent
Pϕ;z,Pϕ;s,Pϕ;a,Pϕ;h Probabilityofarrivalazimuth,slowness,
amplitude,andphasegivenevent
Λ Setofarrivalsassociatedtoevents
Λijk Thearrivalofeventi’sjthphaseatstation
k
Λitjk,Λizjk,Λisjk,Λiajk,Λihjk Arrivaltime,azimuth,slowness,
amplitude,andphase Figure1. Exampleofseismicwaveform(stationASAR,chan-
ζ Aspecialsymbolthatsignifiesthelackof nelSE),STA/LTA,andarrivals.Atthebottomofthelowerpanelthe
anarrival dottedlinesaretheautomaticallydetectedarrivals,whilethesolid
Pω Probabilityoffalsearrivals linesaretheanalystmarkedarrivals.Thecolorversionofthisfigure
Pω;n Probabilityofnumberoffalsearrivals is available only in the electronic edition.
Pω;t,Pω;z,Pω;s,Pω;a, Probabilityoffalsearrivaltime,azimuth,
Pω;h slowness,amplitude,andphase
ξ Allfalsearrivals.ξk arethefalsearrivals x(cid:5) (cid:4)argmaxPθ(cid:1)X (cid:4)x(cid:3)Pϕ(cid:1)Y (cid:4)yjX (cid:4)x(cid:3):
atstationk,andξklisoneofthesefalse EQ-TARGET;temp:intralink-;;313;468 x
arrivals
ξkl,ξkl,ξkl,ξkl,ξkl Falsearrivaltime,azimuth,slowness,
t z s a h Assuming that an algorithm can be devised to solve this
amplitude,andphase
Pγ Probabilitydistributionofcodaarrivals optimization problem, which turns out to be nontrivial, the
Pγ;d Probabilitythatacodaarrivalisdetected key determinant of the success of the Bayesian approach is
Pγ;z,Pγ;s,Pγ;a Probabilityofcodaazimuth,slowness, the accuracy of the probability models Pθ and Pϕ. These
andamplitudegivenpreviousarrival
models embody, in explicit form, elementary knowledge of
Pηaγ;h PCroodbaabarilriitvyalogfecnoedraatepdhabsyearrivalaorζif seismology and seismometry, as well as the residual uncer-
nocoda taintyinherentinthatknowledge.Addingtoandrefiningthe
ηa,ηa,ηa,ηa Codaarrival’stime,azimuth,slowness, knowledge embodied in the models reduces residual uncer-
t z s a
andamplitude tainty and improves system performance.
η Setofallcodaarrivals
Many other researchers have previously applied Baye-
A Setofallobservedarrivals
Δ Great-circledistancebetweeneventiand sian techniques to geophysical problems (see Duijndam,
ik
stationk 1988a,bforexamples).Mostoftheseapplicationsinvolvein-
ferenceoverafixedsetofcontinuous-valuedparameters;the
eventlocalizationprobleminparticularisaddressedbyMyers
historicaldataandpresentanalgorithmthatinfersthesetof etal.(2007).Thefullmonitoringprobleminvolvesinference
seismiceventsgiventheobserveddetections.Insimpleterms, overacombinatorialdiscretespace(thenumberofeventsand
let X be a random variable ranging over all possible collec- theassociationbetweeneventsandobservations),aswellasa
tions of seismic events, with each event defined by time, lo- continuousparameterspaceforeachevent.Inthisrespect,it
cation, depth, and magnitude. Let Y range over all possible resemblesthedataassociationproblemsarisinginmultitarget
signal observations at all seismic stations. Then Pθ(cid:1)X(cid:3) de- tracking (Bar-Shalom andFortmann,1988).
scribes a parameterized generative prior over events, and Our overall project, VISA (vertically integrated seismic
Pϕ(cid:1)YjX(cid:3)describeshowthesignalispropagatedandmeasured analysis), is divided into two stages. The first stage, NET-
(including travel time, selective absorption and scattering, VISA, is the subject of the current paper. As the name sug-
noise, artifacts, sensor bias, sensor failures, etc.). Given ob- gests, NET-VISA deals only with network processing and
servationsY (cid:4)y,weareinterestedintheposteriordistribu- reliesupontheIDC’spre-existingsignaldetectionalgorithm,
tionoverevents,P(cid:1)XjY (cid:4)y(cid:3),whichisgivenbyBayesrule: whichconvertstherawwaveformsintoasequenceofarriv-
P(cid:1)XjY (cid:4)y(cid:3)∝Pθ(cid:1)X(cid:3)Pϕ(cid:1)Y (cid:4)yjX(cid:3): als. Thus, the observations Y correspond to sets of arrivals
EQ-TARGET;temp:intralink-;;55;124 with measured attributes rather than raw waveform signals.
FortheCTBTmonitoringapplication,weareinterestedinthe An example of the signal processing is shown in Figure 1.
value x(cid:5) that maximizes this expression, that is, the most Figure1adisplaysafilteredseismicwaveform,andFigure1b
likely explanation for all the sensor readings: shows the typical short-term average to long-term average
NET-VISA: Network Processing Vertically Integrated Seismic Analysis 3
ratio(STA/LTA)usedtodetectthearrivals.Thedottedlinesin kilometers. Most detectors are sensitive to nanometer-scale
Figure1bshowtheautomaticallydetectedarrivals,whilethe displacements and are quite susceptible to noise.
solid lines show the analyst marked arrivals. The second A typical monitoring system runs station-processing
stage,SIG-VISA,willincorporateasignalwaveformmodel software that preprocesses the raw seismometer measure-
and thereby subsume the detection function. ments, filtering out information that is not of interest and
Asnotedpreviously,NET-VISAcomputesasinglemost computing short-term and long-term averages of the signal
likely explanation: a set of hypothesized events with their amplitude.Whentheratiooftheseaveragesexceedsafixed
associated arrivals, marking all other arrivals as noise. This threshold, an arrival (also called a trigger or onset) is an-
input–outputspecification,whilenotfullyBayesianinspirit, nounced.(PleaserefertoAllen,1978foranearlydescription
enablesdirectcomparisonwiththecurrentautomatedsystem ofthisapproach.)Variousparametersofthearrivalaremea-
bulletin, Standard Event List 3 (SEL3). Using the final sured: onset time, azimuth (direction from the station to the
expert-generated bulletin, LEB, as ground truth, we com- sourceof thewave), slowness (related to the angleof decli-
paredthetwosystemsonsevendaysofheld-outdata.NET- nation of the signal path), amplitude, signal-to-noise ratio
VISAreducesthenumberofmissedeventsby60%without
(SNR),etc.Basedontheseparameters,aphaselabelmaybe
introducingmorefalseevents.Furthermore,takingdatafrom
assigned to the arrival based on the standard International
the more comprehensive NEIC (National Earthquake Infor- Association of Seismology and Physics of the Earth’s
mation Center) bulletin as ground truth for the continental
Interior (IASPEI) phase catalog (Storchak et al., 2003). It is
UnitedStates,wefindthatNET-VISAisabletodetectevents
importanttonotethatnoneofthesemeasurementsareperfect
in the IMS data that are not in the LEB report produced by
and have uncertainty associated with them.
IDC’sexpertanalysts.Wefindsimilarresultsbycorroborat-
Theproblemthatweattempttosolveinthispaperisto
ingwithotherregionalnetworksaswell.Thus,NET-VISA’s
take a continuous stream of arrivals (with onset time, azi-
true performance may be higher than the LEB-based calcu-
muth, slowness, amplitude, SNR, and phase label) from the
lation would suggest. Independent IDC evaluation confirms
roughly120IMSseismicstationsasinputandproduceabul-
the gist of our findings, as well as demonstrates that NET-
letin,acontinuousstreamofeventsandassociationsbetween
VISA can handle large aftershock sequences such as the
events and arrivals. The parameters of an event are its lon-
one following the March 2011 Tohoku earthquake in near
gitude, latitude, depth, time, and magnitude (m or body-
real time. b
wave magnitude). Precise details of the metrics for evaluat-
Intherestofthispaperwedescribetheproblemandour
ingabulletinaregivenintheExperimental Resultssection.
approach in detail including some elementary seismology
(for readers who are not seismologists), the probability
model,theinferencealgorithm,theresults,andfinally,future
Generative Probabilistic Model
plans and conclusions.
A generative model is one that describes the relation-
The Seismic Association and Localization Problem ships among all variables of interest (here, events and ob-
servations) in such a way that complete samples can be
Seismiceventsaredisturbancesintheearth’scrust.Our
generatedfromthejointdistribution.Themodelisbuiltfrom
work is concerned primarily with earthquakes and explo-
severalcomponents,eachofwhichisaconditionalprobabil-
sions (nuclear and conventional), but other types of events,
itydistributiondescribinghowsomeaspectoftheunderlying
suchaswavesbreaking,treesfalling,oricefalling,maygen-
physical processoperates. Forexample, the conditional dis-
erateseismicwaves,too.Allsuchwavesoccurinavarietyof
tributionthat predicts the onset time of an arrival, given the
types:bodywavesthattravelthroughtheearth’sinteriorand
timeoftheeventthatcausedit,embodieswhatseismologists
surfacewavesthattravelonthesurface.Therearetwotypes
call the travel-time model as well as its associated residual
of body waves: compression or P waves and shear or S
uncertainty.
waves.Therearealsotwotypesofsurfacewaves:Loveand
The major components of our model are as follows:
Rayleigh.Further,bodywavesmaybereflectedoffdifferent
layers of the earth’s interior; these are labeled distinctly by • Events: This component describes the rate of occurrence,
seismologists. Each particular wave type generated by a geographical distribution, depth distribution, and magni-
giveneventiscalledaphase.Morethan100distinctphases tude distribution of seismic events (natural and man-
areidentifiedinstandardtables(Storchaketal.,2003);inour made).
work,weconsider14phasesmakingupthevastmajorityof • True arrivals: Each event is assumed to produce the 14
IMS seismic detections: P, Pn, PKP, S, PKPbc, PcP, pP, phasesnotedearlier;thedetectionprobabilitymodelspec-
Lg, PKPab, PKKPbc, Pg, Rg, Sn, and ScP. ifies the chance that a given phase from an event will be
Aftertravelingthroughtheearth,thesewavesarepicked detectedbyagivenstation;thevariouscomponentmodels
up in seismic stations as ground vibrations. Typically, seis- forarrivalattributes(onsettime,azimuth,slowness,ampli-
mic stations have either a single three-axis detector or an tude,andphaselabel)describethepredictedmeasurements
array of vertical-axis detectors spread over a scale of many resultingfromstationprocessingfordetectedphases.Note
4 N. S. Arora, S. Russell, and E. Sudderth
Event Rate and Time. Eventsare assumed to begenerated
bya time-homogeneous Poisson process with a rate param-
eter,λ .Ifeisasetofevents(ofsizejej),andT isthetime
e
periodunderconsideration(inseconds),thepriorprobability
of the number of events, Pθ;n(cid:1)·(cid:3), is given by
(cid:1)λ T(cid:3)jejexp(cid:1)−λ T(cid:3)
EQ-TARGET;temp:intralink-;df1;313;673Pθ;n(cid:1)jej(cid:3)(cid:4) e jej! e : (1)
Foreachevent,ei,theeventtime,ei,isuniformlydistributed
t
between0andT;inotherwords,thepriorprobabilitydensity
of the event time, Pθ;t(cid:1)·(cid:3), is given by
1
Figure2. Estimatingeventrate.Thecolorversionofthisfigure EQ-TARGET;temp:intralink-;df2;313;600Pθ;t(cid:1)eit(cid:3)(cid:4)T: (2)
is available only in the electronic edition.
The parameter, λ , is estimated from the average historical
e
eventfrequenciesasshowninFigure2.Theestimatedvalue
of λ is 0.001266 per second or 4.6 per hour.
e
EventLocation. Thelongitudeandlatitudeoftheithevent,
eil,aredrawnfromaneventlocationdensity,Pθ;l(cid:1)el(cid:3),onthe
surface of the Earth. This density is a mixture of a uniform
density(toallowforexplosionsanddenovoseismicactivity
anywhere on the surface) and a density estimate based on
historical event locations. The historical density is modeled
by a kernel density obtained by placing a kernel function
K (cid:1)·(cid:3) of some fixed width parameter b at every known
b;x
event location, so the overall mixture has the form
1 1 XH
Figure 3. Event location averageleave-one-out log likelihood EQ-TARGET;temp:intralink-;;313;409 Pθ;l(cid:1)el(cid:3)(cid:4):0014πR2(cid:6):999Hh(cid:4)1Kb;ghl(cid:1)el(cid:3);
versusbandwidth.Thecolorversionofthisfigureisavailableonly where R is the Earth’s radius, H is the number of historical
in the electronic edition. events, and g1;…;gH are the locations of those events. We
l l
use an exponential kernel of the form
that measurements may include errors; for example, the
phase label may not match the actual phase. K (cid:1)y(cid:3)(cid:4)1(cid:6)1=b2 exp(cid:1)−Δxy=b(cid:3) ;
• Falsearrivals:Themodeldescribestherateoffalse(noise- EQ-TARGET;temp:intralink-;;313;329 b;x 2πR2 1(cid:6)exp(cid:1)−π=b(cid:3)
generated)arrivalsateachstationandtheirattributedistri-
butions. where Δxy is the great-circle distance between locations x
• Coda arrivals: Any arrival may trigger further spurious and y, measured in radians. Notice that in the limit as
arrivals detected by the station-processing algorithms in b→∞, the kernel tends to the uniform density over the
the coda or tail of the waveform, due to additional peaks Earth’s surface.
in the arriving amplitude. Attributes of coda arrivals are Theoptimalvalueofbischosenbyleave-one-outcross
usually correlated with those of the triggering arrival. validation (LOOCV) over a random subset of 1000 event
locationsinthedataset.Weuseagridsearchtofindthevalue
Thefollowingsectionsdescribethesecomponentsinde-
of b within 0 and 2 that maximizes the LOOCV value. We
tail,givingtheirmathematicalformsandmethodsofparam-
initially use a grid size of 0.2 and refine the search with a
eter estimation.
grid size of 0.05. Figure 3 plots the results of this search.
Thebestvalueofbthatischosenis0.7.InFigure4weshow
Events the event location prior log(cid:7)Pθ;l(cid:1)·(cid:3)(cid:8), using all of the training
dataset and the optimal bandwidth.
In the following, we only consider events with body- Thedepthoftheevent,ei,isuniformlydistributedupto
d
wavemagnitudem 2orhigher.Allothereventsareconsid- a maximum depth D (700 km in our experiments); in other
b
erednoise,becausetheyaretoosmalltoberelevantfortreaty words, the probability density of the event depth, Pθ;d(cid:1)·(cid:3), is
monitoring purposes. given by
NET-VISA: Network Processing Vertically Integrated Seismic Analysis 5
(NotethatweareoverloadingPθ(cid:1)·(cid:3)torefertothedistribution
over the set of events as well as the distribution of a single
event.) Substituting from equations (1) and (5),
Yjej 1
EQ-TARGET;temp:intralink-;df6;313;697Pθ(cid:1)e(cid:3)(cid:4)exp(cid:1)−λeT(cid:3) Pθ;l(cid:1)eil(cid:3)Dλeλmexp(cid:7)−λm(cid:1)eim−2(cid:3)(cid:8):
i(cid:4)1
(6)
If we define the following single-event quantity, which is
independent of the arbitrary interval T,
1
EQ-TARGET;temp:intralink-;;313;606 P^θ(cid:1)ei(cid:3)(cid:4)Pθ;l(cid:1)eil(cid:3)Dλeλmexp(cid:7)−λm(cid:1)eim−2(cid:3)(cid:8):
Figure 4. Eventlocationlogdensity.Thecolorversionofthis
then we can simplify equation (6) to
figure is available only in the electronic edition.
Yjej
1 Pθ(cid:1)e(cid:3)(cid:4)exp(cid:1)−λeT(cid:3) P^θ(cid:1)ei(cid:3): (7)
EQ-TARGET;temp:intralink-;df3;55;524Pθ;d(cid:1)eid(cid:3)(cid:4)D: (3) EQ-TARGET;temp:intralink-;df7;313;551 i(cid:4)1
Asnoted in the Future Worksection, the time-homogeneity
and independence assumptions are violated by aftershock
EventMagnitude. TheGutenberg–Richterlaw(Gutenberg
phenomena.
and Richter, 1954) posits that the number of events with
magnitudemormoreis10timesthenumberofeventswith
magnitude m(cid:6)1 or more. In terms of the event magnitude True Arrivals
probability of an arbitrary event: Pθ;m(cid:1)eim ≥m(cid:3)(cid:4) An event can generate at most one true arrival of each
10·Pθ;m(cid:1)eim ≥m(cid:6)1(cid:3). We represent this prior knowledge phasetypeatastation.Whetherornotthearrivalisgenerated
in our model with an exponential distribution defined for depends on the detection probability.
ei ≥2 with rate λ (cid:4)log(cid:1)10(cid:3):
m m
DetectionProbability. Theprobabilitythataneventi’sjth
Pθ;m(cid:1)eim(cid:3)(cid:4)λmexp(cid:7)−λm(cid:1)eim−2(cid:3)(cid:8): (4) phase, 1≤j≤J, is detected by a station k, 1≤k≤K,
EQ-TARGET;temp:intralink-;df4;55;383
dependsonthephase,thestation,andtheevent’smagnitude,
Notethatbecauseweassumeaminimummagnitudeof2,the depth,anddistancetothestation(ei ,ei,andΔ ).Asnoted
m d ik
exponential distribution is shifted forward by this amount. in the Future Work section, such a model ignores source
Although man-made events may not follow the Guten- location and path effects as well as anisotropic radiation
berg–Richter law, their frequency is too low to affect the patterns.LetPjk (cid:1)e(cid:3)betheprobabilityofthisdetection.The
ϕ;d i
overall distribution significantly. phase- and station-specific detection distributions, Pjk (cid:1)·(cid:3),
ϕ;d
wereobtainedusinglogistic regression models, whichare a
Overall Event Prior. Under the assumption that the event standard tool for modeling true/false random variables. In
location,depth,time,andmagnitudeareindependentofeach suchmodels,aweightedlinearexpressionisformedfromthe
other, the probability of a single event, ei, is given by the inputs (and possibly additional features computed from the
product of the component terms: inputs), and the probability that the output is true; in this
case, the detection of the phase is given by applying a soft
EQ-TARGET;temp:intralink-;;55;229 Pθ(cid:1)ei(cid:3)(cid:4)Pθ;l(cid:1)eil(cid:3)Pθ;d(cid:1)eid(cid:3)Pθ;t(cid:1)eit(cid:3)Pθ;m(cid:1)eim(cid:3): thresholdlogisticfunctiontothevalueoftheweightedlinear
expression.Theneteffectisthatthelogoddsofdetectionare
Substituting from equations (2), (3), and (4), we get a linear function of the input features:
EQ-TARGET;temp:intralink-;df5;55;183Pθ(cid:1)ei(cid:3)(cid:4)Pθ;l(cid:1)eil(cid:3)D1T1λmexp(cid:7)−λm(cid:1)eim−2(cid:3)(cid:8): (5) EQ-TARGET;temp:intralink-;;313;180 log(cid:2)1−PjϕPk;dj(cid:1)ke(cid:1)i(cid:3)ei(cid:3)(cid:3)(cid:4) Xμwdjk·w(cid:1)eim;eid;Δik(cid:3);
ϕ;d w∈Fd
In our model, all the events are exchangeable and are gen-
erated independently, thus, where Fd is a set of feature functions such that for each
w∈F , w: R3 →R. Also, μwjk is theweight for the feature
d d
Yjej w.ThecompletesetoffeaturesisdefinedinTable2.Because
Pθ(cid:1)e(cid:3)(cid:4)Pθ;n(cid:1)jej(cid:3)·jej!· Pθ(cid:1)ei(cid:3): the event magnitude is one of the features, another way of
EQ-TARGET;temp:intralink-;;55;116
i(cid:4)1 thinking about the previous equation is that a unit increase
6 N. S. Arora, S. Russell, and E. Sudderth
Table 2
Detection Features*
Feature Value
(Intercept) 1
mag ei
m
depth ei
d
dist Δ
ik
dist0 N(cid:1)Δ ;0;5(cid:3)
ik
dist35 N(cid:1)Δ ;35;20(cid:3)
ik
dist40 N(cid:1)Δ ;40;20(cid:3)
ik
dist12520 N(cid:1)Δ ;125;20(cid:3)
ik
dist12540 N(cid:1)Δ ;125;40(cid:3)
ik
dist145 N(cid:1)Δ ;145;10(cid:3)
ik
dist170 N(cid:1)Δ ;170;20(cid:3)
ik
dist175 N(cid:1)Δ ;175;30(cid:3)
mag6 N(cid:1)eiik;6;5:5(cid:3) Figure5. ConditionaldetectionprobabilitiesforthePphaseof
m
mag68 N(cid:1)ei ;6;8(cid:3) surface events between 3 and 4 mb at station ASAR. The color
m
md (cid:1)7−ei (cid:3)·Δ version of this figure is available only in the electronic edition.
m ik
*List of features used for computing the probability of
detecting an arrival from an event i with magnitude ei , depth haveapositiveweightforevent magnitude forthePphase,
m
ei, and distance Δ from station k. Here N(cid:1)x;μ;σ(cid:3) is the then a station with very little data will also have a positive
d ik
standard Gaussian density with mean μ and standard deviation weight for this feature.
σmeasuredatx. In Figure 5 we show the model prediction that is, the
probability of detection, for one phase at a station.
intheeventmagnitudewouldresultintheoddsofdetection
increasingbyamultiplicativeconstant.Theexactconstantis Arrival Attributes. If event i’s jth phase is detected by a
the exponentiation of the corresponding weight of the mag- station k, we define Λijk as the corresponding arrival. For
nitude feature, and this is station dependent. In fact, if the notational convenience, we define Λijk (cid:4)ζ whenever this
weight was negative then the detections odds would de- phaseisnotdetected.Ourmodelspecifiesaprobabilitydis-
crease; this is, in fact, the case for the travel-time feature. tributionforthemeasuredattributesofthisarrival:theonset
alwaDysirpecotslsyibelestbimecaatiunsge mthaenyfeoaftutrhee swtaetiigohnt–sphμawdsjek ciosmnboi-t tΛimiajke,Λanitjdk,tthheeaaszsimiguntehdΛpizhjka,sethleabsleolwΛnihejks.sTΛhisjek,atrhreivaaml tpimliteudise
nations have very little data. To deal with this data sparsity
we used a hierarchical Bayesian procedure (Gelman et al.,
2004), which posits that the weight for a station–phase is EQ-TARGET;temp:intralink-;;313;367 Λitjk (cid:4)eit(cid:6)IjT(cid:1)eid;Δik(cid:3)(cid:6)ritjk;
drawn from a global prior for that phase. This global prior
is,inturn,drawnfromaweaklyinformativeprior,asfollows:
whereIj isthepredictionfromtheIASPEItravel-timemodel
T
(afunctionoftheeventdepthanddistancetothestation),and
μwjk∼N(cid:1)μwj;σwj(cid:3) μwj∼N(cid:1)0;100(cid:3) rijk is a random residual. The residual distribution accounts
EQ-TARGET;temp:intralink-;;55;306 d d d d fotr the inhomogeneities in the Earth’s crust, which allow
(cid:1)σwj(cid:3)−2∼Γ(cid:1)0:01;100(cid:3);
d seismicwavestotravelfasterorslowerthantheIASPEIpre-
diction. This distribution also accounts for any systematic
where N and Γ are the Gaussian and the gamma distribu- biases in picking seismic onsets from waveforms. Whereas
tions parameterized by their mean and standard deviation, most authors assume Gaussian residuals, an assumption
and shape and scale, respectively. Estimation of parameters implicitintheuseofquadraticcostfunctionsintheGAsys-
followsacoordinateascentprocedure.Foreachphasej,we tem,wefindexperimentallythattraveltimeandindeedmost
initializeμwj (cid:4)0andσwj (cid:4)1,andthenalternatelyoptimize other residuals are distributed according to a Laplacian
d d
μwjk:∀w,k,μwj:∀w,andσwj:∀wuntilconvergence.Ineach method. The parameters of a Laplacian distribution are
d d d
maximizationstep,theoptimalvalueofμwjkiscomputedby the mean μjk and scale bjk, and the distribution is given by
d t t
second-order optimization, whilethe remaining values have
a closed-form solution.
(cid:2) (cid:3)
For each phase, the previously described procedure 1 jΛijk−ei−I (cid:1)ei;Δ (cid:3)−μjkj
Pjk(cid:1)Λijkjei(cid:3)(cid:4) exp − t t T d ik t :
ensures that if a station has a lot of data (both detections EQ-TARGET;temp:intralink-;;313;176 ϕ;t t 2bjk bjk
t t
and nondetections), then its feature weights will be deter-
mined almost entirely by its own data, whereas feature
weights for data-poor stations tend toward a global average Similarly, the arrival azimuth and slowness follow a Lapla-
obtained from all stations. Thus, if most data-rich stations cian distribution:
NET-VISA: Network Processing Vertically Integrated Seismic Analysis 7
(cid:2) (cid:3)
1 jψ(cid:7)Λijk;G (cid:1)sk;ei(cid:3)(cid:8)−μjkj
Pjk (cid:1)Λijkjei(cid:3)(cid:4) exp − z z l l z ;
EQ-TARGET;temp:intralink-;;313;733 ϕ;z z 2bjk bjk
z z
(cid:2) (cid:3)
1 jΛijk−Ij(cid:1)ei;Δ (cid:3)−μjkj
Pjk (cid:1)Λijkjei(cid:3)(cid:4) exp − s s d ik s :
ϕ;s s 2bjk bjk
s s
Herethefunctionψcomputesthedifferenceintheobserved
angle Λijk and the angle computed from the geographical
z
function G , which depends on the station location, sk, and
z l
the event location, ei. Also, Ij is the slowness value com-
l s
puted from the IASPEI model for phase j. Technically, the
factthatthedomainofLaplacedistributionsistheentirereal
linemakestheminappropriateforazimuthandslownessvar-
iables; but the model is reasonable given the small scale of
Figure6. ArrivaltimedistributionaroundtheIASPEIprediction the residuals. It should be noted that the observedvalues of
forP-phasearrivalsatstationASAR.Thecolorversionofthisfig- azimuth and slowness referred to already include certain
ure is available only in the electronic edition. CTBTO-defined, station-specific corrections called slowness
azimuth site corrections (SASC).
The estimation of all the station-and-phase-specific
meansandscales,forexampleμjkandbjk,isbasedonahier-
t t
archical model.Inthismodel,eachmeanand scale isgener-
ated from a phase-global prior, which we describe for the
arrival time in the following equation:
rijk∼Laplace(cid:1)μjk;bjk(cid:3) μjk∼Laplace(cid:1)μj;bj(cid:3)
t t t t t t
EQ-TARGET;temp:intralink-;;313;465
(cid:1)bjk(cid:3)−1∼Γ(cid:1)1;βj(cid:3) μj∼Laplace(cid:1)0;100(cid:3)
t t t
(cid:1)bj(cid:3)−1∼Γ(cid:1):01;100(cid:3) (cid:1)βj(cid:3)−1∼Γ(cid:1):01;100(cid:3):
t t
Foreachphasej,westartbyinitializingμjk (cid:4)0,bjk (cid:4)1for
t t
all stations k, μj (cid:4)0, bj (cid:4)1, and βj (cid:4)1. Next, we itera-
t t t
tively optimize the values of μjk, bjk for each station, and
t t
theglobalvaluesμj,bj,andβj.Thisisrepeateduntilconver-
t t t
gence. Each of these optimization steps has a simple closed-
Figure7. Arrivalazimuthdistributionaroundthegeographical
form solution. A similar procedure is adopted for estimating
prediction (plus SASC correction) for P-phase arrivals at station theazimuthparameters,μjk andbjk,andtheslownessparam-
ASAR.Thecolorversionofthisfigureisavailableonlyintheelec- z z
tronic edition. eters,μjsk,andbjsk.Examplesforthesethreetypesofdistribu-
tions are shown in Figures 6, 7, and 8.
Theeffectofthehierarchicalmodelissimilarinspiritto
themodelfordetectionprobabilities.Foreachphase,stations
with a lot of data will compute the mean and scale almost
exclusively from their own data, whereas parameters for
data-poor stations will tend to reflect the global average.
The arrival amplitude Λijk is similar to the detection
a
probabilitybecauseitdependsonlyontheeventmagnitude,
depth, and distance to the station. We model the log of the
amplitudeviaalinearregressionmodelwithGaussiannoise,
X
rijk (cid:4)log(cid:1)Λijk(cid:3)− μwjk·w(cid:7)ei ;ei;Ij(cid:1)ei;Δ (cid:3)(cid:8);
EQ-TARGET;temp:intralink-;;313;188 a a a m d T d ik
w∈Fa
rijk∼N(cid:1)0;σjk(cid:3);
a a
where F is a set of feature functions (see Table 3), and
a
μwjk is the weight for feature w. This implies that
a
Figure 8. ArrivalslownessdistributionaroundtheIASPEIpre- (cid:5) (cid:6)
diction(plusSASCcorrection)forP-phasearrivalsatstationASAR. 1 (cid:1)rijk(cid:3)2 1
Tedhietiocno.lor version of this figure is available only in the electronic EQ-TARGET;temp:intralink-;;313;111 Pjϕk;a(cid:1)Λiajkjei(cid:3)(cid:4)p(cid:4)2(cid:4)(cid:4)π(cid:4)(cid:4)(cid:4)σjakexp − 2σajak2 Λiajk:
8 N. S. Arora, S. Russell, and E. Sudderth
Table 3
Amplitude Features*
Feature Value
(Intercept) 1
mag ei
m
depth ei
d
ttime Ij(cid:1)ei;Δ (cid:3)
T d ik
ttime0 N(cid:7)Ij(cid:1)ei;Δ (cid:3);0;50(cid:8)
T d ik
*Listoffeaturesusedforpredictingtheamplitudeofaphase,j,
arrivalfromanevent,i,atastation,k.
Inordertoestimatethefeatureweights,weuseahierarchical
model,asbefore,whichassumesthatforeachphasethefea-
ture weights at a station are drawn from a global prior: Figure9. ArrivallogamplitudedistributionforthePphaseofa
surfaceeventof3–4mbatstationASAR.Thecolorversionofthis
μwjk∼N(cid:1)μwj;σwj(cid:3) (cid:1)σjk(cid:3)−2∼Γ(cid:1)100;βj(cid:3) figure is available only in the electronic edition.
a a a a a
EQ-TARGET;temp:intralink-;;55;562
μwj∼N(cid:1)0;100(cid:3) (cid:1)σwj(cid:3)−2∼Γ(cid:1):01;100(cid:3)
a a
(cid:1)βj(cid:3)−1∼Γ(cid:1):01;100(cid:3):
a
Maximum a posteriori (MAP) inference of these weights is
along the lines of the method already described for other
parameters. Figure 9 has an example of a learned model.
Finally,thedistributionfortheassignedphaselabelΛijk
h
reflects the phase classification errors made by station
processing. We model this with a categorical distribution
whose parameters depends only on the true phase, j:
Pjk (cid:1)Λijkjei(cid:3)(cid:4)Pj(cid:1)aijk(cid:3):
EQ-TARGET;temp:intralink-;;55;404 ϕ;h h h h
Learning the categorical distributions Pj is a simple matter Figure 10. Arrival phase probability as function of the true
h
of counting with add-one smoothing. The learned distribu- phase.Thecolorversionofthisfigureisavailableonlyintheelec-
tionisplottedasaheatmapinFigure10.Theresultsindicate tronic edition.
that many phases are falsely classified as P by the station-
processing software.
The indicator variable on the condition Λijk (cid:4)ζ selects the
appropriate term.
Overall True Arrivals. We assume that all the measured
Theseconditionalindependenceassumptionsareimpor-
attributesofanarrivalareconditionallyindependentofeach
tantinreducingthenumberofmodelparameters,buttheyare
other, given the true event. Thus,
certainlyopentoquestion.Consideringarrivalsattwodiffer-
EQ-TARGET;temp:intralink-;;55;260 Pjϕk(cid:1)Λijkjei(cid:3)(cid:4)Pjϕk;t(cid:1)Λitjkjei(cid:3)Pjϕk;z(cid:1)Λizjkjei(cid:3)Pϕ;s(cid:1)Λisjkjei(cid:3) entstations,theassumptionimpliesthattheresidualsareun-
×Pjk (cid:1)Λijkjei(cid:3)Pjk (cid:1)Λijkjei(cid:3): correlated;thisisreasonablebecauseIMSstationsarewidely
ϕ;a a ϕ;h h separated. For arrivals of two different phases at the same
Further,assumingthatthearrivalsareallconditionallyinde- station,theassumptionofuncorrelatedresidualsismoreprob-
pendent, and keeping in mind that the observed arrivals Λ lematic. Our analysis of travel-time residuals for P and S
include by implication all the nondetections for every miss- phasesdoesshowsomecorrelation,butitisweakenoughthat
ing station–phase combination, we have ignoring itdoes not seriously compromise performance.
Each station k also generates a set ξk of false arrivals,
Yjej YJ YK
EQ-TARGET;temp:intralink-;dPf8;55;168 ϕ(cid:1)Λje(cid:3)(cid:4) i(cid:4)1j(cid:4)1k(cid:4)1(cid:7)1(cid:1)Λijk (cid:4)ζ(cid:3)(cid:7)1−Pjϕk;d(cid:1)ei(cid:3)(cid:8) tahreatdiess,carrirbievdalsbynoatnaosvsoercaiallteddiswtriibthutaionnyPevkωe.nTthpehansuem. bTehresoef
falsearrivalsisdistributedaccordingtoatime-homogeneous
(cid:6)1(cid:1)Λijk
Poisson process with rate λk:
f
≠ζ(cid:3)Pjk (cid:1)ei(cid:3)Pjk(cid:1)Λijkjei(cid:3)(cid:8): (8)
ϕ;d ϕ
(cid:1)λkT(cid:3)jξkjexp(cid:1)−λkT(cid:3)
Tbihlietyteorrmaidnestiedcetitohnepprroobdaubcitliitsyetiimtheesr tahenoanrrdievtaelcptiroonbapbriolbitay-. EQ-TARGET;temp:intralink-;;313;105 Pkω;n(cid:1)jξkj(cid:3)(cid:4) f jξkj! f :
NET-VISA: Network Processing Vertically Integrated Seismic Analysis 9
1 1 1
EQ-TARGET;temp:intralink-;;313;733 Pkω;t(cid:1)ξktl(cid:3)(cid:4)T; Pkω;z(cid:1)ξkzl(cid:3)(cid:4)M ; Pkω;s(cid:1)ξksl(cid:3)(cid:4)M ;
z s
where M and M are the range of values for azimuth and
z s
slowness, respectively. The log amplitude log(cid:1)ξkl(cid:3) of the
a
false arrival is generated from either a uniform distribution
withprobability0.1oramixtureoftwoGaussians,whichis
estimatedfromthedatawithastandardexpectationmaximi-
zation procedure. This distribution, fkω;a(cid:1)·(cid:3), at one station is
displayed in Figure 12.
False Arrivals
The resulting distributionon thefalsearrival amplitude
Figure 11. The average false arrival rate per hour at all the is given by
stations.Thecolorversionofthisfigureisavailableonlyintheelec-
1
tronic edition. EQ-TARGET;temp:intralink-;;313;561 Pkω;a(cid:1)ξkal(cid:3)(cid:4)fkω;a(cid:7)log(cid:1)ξkal(cid:3)(cid:8)ξkl:
a
Finally, the phase label ξkl assigned to the false arrival fol-
h
lows a categorical distribution, Pk (cid:1)ξkl(cid:3). This is estimated
ω;h h
via the posterior mean under a uniform Dirichlet prior
(see Fig. 13).
Overall, assuming the false arrival attributes are inde-
pendently generated, we have
EQ-TARGET;temp:intralink-;;313;455 Pkω(cid:1)ξkl(cid:3)(cid:4)Pkω;t(cid:1)ξktl(cid:3)Pkω;z(cid:1)ξkzl(cid:3)Pkω;s(cid:1)ξksl(cid:3)Pkω;a(cid:1)ξkal(cid:3)Pkω;h(cid:1)ξkhl(cid:3):
Because the false arrivals at a station are exchangeable, the
probability for the set ξk is
Yjξkj
Pkω(cid:1)ξk(cid:3)(cid:4)Pkω;n(cid:1)jξkj(cid:3)·jξkj! Pkω(cid:1)ξkl(cid:3)
EQ-TARGET;temp:intralink-;;313;402
l(cid:4)1
Figure 12. Amplitude distribution for false arrivals at station Yjξkj λk
ASAR.Thecolorversionofthisfigureisavailableonlyintheelec- (cid:4)exp(cid:1)−λkfT(cid:3) M Mf Pkω;a(cid:1)ξkal(cid:3)Pkω;h(cid:1)ξkhl(cid:3):
tronic edition. l(cid:4)1 z s
As before, we have overloaded Pkω to refer to a distribution
overasetofarrivalsaswellasasinglearrival.Asbefore,we
can simplify the equations by defining
λk
EQ-TARGET;temp:intralink-;;313;288 P^kω(cid:1)ξkl(cid:3)(cid:4)M Mf Pkω;a(cid:1)ξkal(cid:3)Pkω;h(cid:1)ξkhl(cid:3):
z s
Now,assumingthatthefalsearrivalsatdifferentstationsare
independent of each other, the probability for the complete
set ξ of all false arrivals at all stations is
(cid:2) (cid:3)
XK YK Yjξkj
EQ-TARGET;temp:intralink-;df9;313;215Pω(cid:1)ξ(cid:3)(cid:4)exp − λkfT P^kω(cid:1)ξkl(cid:3): (9)
k(cid:4)1 k(cid:4)1l(cid:4)1
Coda Arrivals
Figure 13. Phase distribution for all false arrivals. The color
Followingtheinitialonsetandpeak,thearrivingenergy
version of this figure is available only in the electronic edition.
from a seismic phase does not decay smoothly,butexhibits
many subsidiary peaks, some of which may even be larger
Thevalues of λk are displayed in Figure 11.If ξkl is one of than the initial peak. These subsidiary peaks may fool the
f
this set of false arrivals, its time ξkl, azimuth ξkl, and slow- STA/LTAdetectorintoannouncingadditionalarrivals,called
t z
nessξklaregenerateduniformlyovertheirrespectiveranges: codaarrivals;intheIMSdata,uptohalfofallarrivalsareof
s
10 N. S. Arora, S. Russell, and E. Sudderth
Figure 14. Codadetectionprobabilityasafunctionofthelog Figure 16. Coda azimuth difference versus triggering arrival.
The color version of this figure is available only in the electronic
amplitudeofthetriggeringarrival.Thecolorversionofthisfigureis
edition.
available only in the electronic edition.
Figure 15. Time delay for coda arrival after the triggering Figure 17. Codaslownessdifferenceversustriggeringarrival.
arrival.Thecolorversionofthisfigureisavailableonlyintheelec- The color version of this figure is available only in the electronic
tronic edition. edition.
thistype.AnexampleofcodaarrivalscanbeseeninFigure1
One problem that arises while training the model for
around 20 s after the main arrival, which is at offset 100 s.
coda arrivals is that the IDC analysts do not annotate coda
Onemightimaginethatcodaarrivalscanbelumpedinwith
arrivalsintheLEBbulletin,whichmakesithardtoestimate
falsearrivals,butitturnsoutthattheattributesofcodaarriv-
theparametersofPγ.Oursolutionistoheuristicallyannotate
alsarestronglycorrelatedwiththoseofthetriggeringarrival.
someoftheunassociatedarrivalsascodaarrivalsandusethis
If the coda arrivals are not modeled explicitly, then our in-
annotation to learn Pγ and also Pω.Our procedure searches
ferencewillenduphypothesizingadditionalspuriousevents
thetrainingdataforanyunassociatedarrivalswithin30sofa
asthemostlikelyexplanationformanyofthecodaarrivals.
prior arrival at the same station and with an azimuth and
Wemodeltheparametersofthecodaarrivalwithadis-
slowness within 50° and 10 s per degree, respectively, of
tributionPγ,whichisafunctionoftheparametersofthetrig- the prior arrival’s values, and to mark sucharrivals as coda.
gering arrival. Whether the triggering arrival was a false
Wemodelthedistributionofthetimedelaybetweenthe
arrival, or caused by an event, or itself triggered by another coda arrival and the triggering arrival as a gamma distribu-
arrival,isimmaterial.Wedefineηaasthearrivaltriggeredby tion (see Fig. 15):
arrival a (at the same station), or ζ if there is no such trig-
geredarrival.Now,theprobabilitythatanarrivalatriggersa ηa−a ∼Γ(cid:1)ρ;ν (cid:3); i:e:
t t t t
EQ-TARGET;temp:intralink-;;313;156 (cid:2) (cid:3)
coda arrival is given by Pγ;d(cid:1)a(cid:3), which is a function of the 1 ηa−a
amplitudeofthearrivala.WeestimatePγ;d withanonpara- Pγ;t(cid:1)ηatja(cid:3)(cid:4)Γ(cid:1)ρ(cid:3)νρt(cid:1)ηat −at(cid:3)ρt−1exp − t ν t :
metricmodelbydiscretizingthelogamplitudeofthetrigger- t t t
ingarrivalintobucketsofsize0.25between−4and10.This The differences in azimuth, slowness, and log amplitude of
distribution is displayed in Figure 14. Any points outside the coda versus the triggering arrival are all modeled as
these extreme values are mapped to the nearest bucket. Laplace distributions (Figs. 16, 17, and 18):