A Prospect of Earthquake Prediction Research

Statistic al Scienc e 2013, V ol. 28, No. 4, 521– 541 DOI: 10.1214 /13-STS439 c  Institute of Mathematical Statisti cs , 2013 A Prosp ect of Ea rthquak e Prediction Resea rch Y osihik o Ogata Abstr act. Earthquak es o ccur b ecause of abrupt slip s on faults d u e to accum u lated stress in the Earth’s cr u st. Because most of th ese faults and their mec han ism s are not readily apparen t, deterministic earth- quak e prediction is diﬃcult. F or eﬀect iv e prediction, complex co ndi- tions a nd uncertain e lement s m ust b e considered , whic h necessitates sto c hastic p rediction. In particular, a large amount of u ncertain t y lies in ident ifying wh ether abn ormal p henomena are precursors to large earthquak es, as well as in assigning urgency to the earthquak e. An y disco v ery of p oten tially useful information for earthquak e prediction is incomplete un less qu an titativ e mo d eling of risk is considered. Th ere- fore, this m anuscript describ es the prosp ect of earthquake predictabilit y researc h to realize pr actical op erational forecasting in the near future. Key wor ds and phr ases: Abnormal p henomena, aseismic slip, Ba yesian constrain ts, epidemic-t yp e aftersho c k sequence (ET AS) m o dels, h ier- arc h ical space–time ET AS mo dels, p robabilit y forecasts, p robabilit y gains, stress c hanges. 1. INTRODUCTION Through r emark able dev elopment s in solid Earth science since the late 1960s, our understandin g of earthquak es h as increased signiﬁcantly . Th e a v ail- abilit y of relev ant data has s teadily in creased as the study of ea rthqu ak es has pr ogressed remark- ably in g eophysic s. After ev ery ma jor earthquake, researc h ers h a v e elucidated imp ortan t seismic mec h- anisms asso ciated with it. Ho wev er, ev en though detailed an alysis and discussions ha v e b een con- ducted, large uncertain ties remain b ecause of dive r- sit y and complexit y of the earthquake phenomenon. Y osihiko O gata is Pr ofessor Emeritus, Institu t e of Statistic al Mathematics, Information and System R ese ar ch O r ganization, 10-3 Midori-cho, T achi kawa, T okyo 190-8562 and Institut e of In dustrial Scienc e, University of T okyo, Visiting Pr ofessor, 4-6-1 Komab a, Me gr o-Ku, T okyo 153-850 5 e-mail: o gata@ism.ac.jp . This is a n electronic r e print of the orig inal article published by the Institute of Mathematical Sta tistics in Statistic al Scienc e , 20 13, V ol. 28, No. 4 , 521– 5 41 . This reprint diﬀers from the original in pa gination and t yp ogr aphic detail. This leads to un ac h iev able c hallenges in determin- istic earthquak e prediction b ecause all diverse and complex scenarios must faithfu lly reﬂect the pro- cesses of earthquakes to b e considered for eﬀectiv e earthquak e pr ediction. On the other hand, sev eral tec hn iques for predict- ing earthquakes hav e b een p r op osed on the basis of anomalies of v arious t yp es; ho wev er, the eﬀec- tiv en ess of these tec hniques is contro v ersial (Jor- dan et al. ( 2011 )). Th erefore, ob jectivit y is required for suc h ev aluation; otherwise, argu m en ts p r esen ted ma y lac k merit. New prediction m o dels that claim to incorp orate p otentia lly useful in formation o ver those used in stand ard seismicit y mo dels s hould b e ev aluated to determine whether predictiv e p o we r is impro v ed. Earthquak e forecasting mod els should ev olv e in this manner . Recen tly , there has b een gro wing momentum for seismologist s to dev elop an organized research program on earthquak e pr edictabilit y . An interna- tional co op erativ e study kn o w n as Col lab or atory for the Study of Earthquake Pr e dictability (CSEP; http://w ww.csept esting.org/ ) is current ly und er w a y among countries prone to ma jor earthquakes 1 2 Y. OGA T A for exploring p ossib ilities in earthqu ak e prediction (e.g., Jordan ( 2006 )). An immediate ob j ectiv e of the pro ject is to encourage the d ev elopmen t of statisti- cal mo dels of seismicit y , su c h as those sub sequent ly discussed in Section 2 , and to ev aluate their predic- tiv e p erformances in terms of pr obabilit y . In addition, the C SEP stu d y aims to deve lop a scien tiﬁc inf rastructure to ev aluate statistical signif- icance and pr ob ability gain (Aki ( 1981 )) of v arious metho ds used to p r edict large earthquak es by using observ ed abnormalities su c h as seismicit y anomaly , transien t cru stal mov emen ts and electromagneti c anomaly . Here pr ob ability gain is deﬁned as the ra- tio of p robabilit y of a large earthquake estimated based on an anomaly to the underlying pr ob ab ility without anomaly . Section 3 describ es this imp ortan t concept, and then discus ses s tatistical p oin t-pro cess mo dels to examine th e signiﬁcance of causalit y of anomalies and also to ev aluate the pr obabilit y gains conditional on the anomalous eve nts. F or pred iction of large earthquak es with a h igher probabilit y gain, comprehensiv e studies of anoma- lous phenomena and observ ations of earthquak e mec h anisms are essential . Sev eral su c h studies are summarized in S ections 4 – 6 . P articularly , I ha v e b een in terested in elucidating abn ormal seismic ac- tivities and geod etic anomalies to apply th em for promoting forecasting abilities, as describ ed in these sections. 2. PROBABILITY F ORECASTING OF BASELINE S EISMICITY 2.1 Log-Lik eliho o d fo r the Evaluation S co re of Probabilit y F o recast Through rep eated revisions, CSEP atte mpts to es- tablish standard m o dels to predict probabilit y th at conform to v arious parts of the w orld. Here I mean the prediction/estimation of pr ob ab ility as p redict- ing/estimating conditional probabilities gi ve n the past history of earthquake s and other p ossible pr e- cursors. The lik eliho o d is u s ed as a reasonable measure of p r ediction p erf orm ance (cf. Boltzmann ( 1878 ); Ak aik e ( 1985 )). The ev aluation metho d for probabilistic forecasts of earthquake s by the log- lik eliho o d fun ction has b een pr op osed, discussed and imp lemen ted (e.g., Kagan and Jac kson ( 1995 ); Ogata ( 1995 ); Ogata, Utsu and Katsura, 1996 ; V er e- Jones ( 1999 ); Harte and V er e-Jones ( 2005 ); S c hor- lemmer et al. ( 2007 ); Zec har, Gerstenb erger and Rhoades, 2010 ; Nanjo et al. ( 2012 ); Ogata et al. ( 2013 )). In some suc h stud ies, the ev aluation score has b een r eferr ed to as (relativ e) entrop y , wh ic h is essen tially similar to the log-lik eliho o d. 2.2 Space–Time-Magnitude Fo recasting of Ea rthquak es Baseline mo d els should b e set to compare with and ev aluate all pr edictabilit y mo dels. Based on em- pirical la ws, we can pr edict s tand ard reference prob- abilit y of earthqu ak es in a space–time-mag nitude range on the b asis of the time series of present and past earthquak es. The framew ork of CSEP , w hic h has ev aluated p erformances of subm itted forecasts of resp ectiv e regions (Jordan ( 2006 ); Zechar, Ger- sten b erger and Rhoades, 2010 ; Nanjo et al. ( 201 1 )), is similar to that of the California Regional E arth - quak e Lik eliho o d Mo dels (REL M) pro j ect for spatial forecast (Field ( 2007 ); Sc horlemmer et al. ( 2010 )). Diﬀeren t s p ace–time mo dels were submitted to the CSEP Japan T esting Center at the Earthquake Re- searc h In stitute, Univ ersit y of T oky o, for the one- da y forecast applied to the testing region in Japan (Nanjo et al. ( 2012 )). This means that the mo d el forecasts the pr obabilit y of an earthqu ake at eac h space–time-mag nitude bin. Ho w ev er, the CSE P pro- to cols h av e to b e impr ov ed to those in terms of p oint- pro cesses on a con tin uous time axis for the ev aluation includin g a real-time forecast (Ogata et al. ( 2013 )). Almost all mo dels incorp orated the Gutenberg– Ric hter (G–R) law for forecasting the m agnitude factor, and tak e diﬀeren t v arian ts of the space–time epidemic-t yp e aftersh o c k sequence (ET AS) mo del (Nanjo et al. ( 2012 ), and Ogata et al. ( 2013 )). In the follo w ing sections, I will outline these mo dels. 2.2.1 Magnitude fr e quency distribution Guten- b erg and Rich ter ( 1944 ) determined th at the num- b er of earthqu ak es increased (decreased) exp onen- tially as their magnitude decreased (increased). De- scribing this theory in terms of p oint pro cesses, the in tensit y of magnitude M is λ 0 ( M ) = lim ∆ → 0 1 ∆ Prob( M < Magnitude ≤ M + ∆) (1) = 10 a − bM = Ae − β M for constants a and b . In other words, the magnitude of eac h earthquak e will ob ey an exp onen tial d is- tribution su c h that f ( M ) = β e − β ( M − M c ) , M ≥ M c , where β = b ln 10, and M c is a cutoﬀ magnitud e v alue ab o v e whic h all earthqu ak es are d etected. T ra- ditionally , the b -v alue had b een estimate d graph - EAR THQUAKE PREDICTION 3 Fig. 1 . T op p anel shows Delaunay tesse l lation up on which the pie c ewise line ar f unction is deﬁne d. The smo othness c on- str aint is p ose d i n the sum of squar es of inte gr ate d sl op es. Bottom p anel shows b-value estimates of the G–R f ormula in e quation ( 1 ) estimate d fr om the data f or e arthquakes of M ≥ 5 . 4 fr om the Harvar d University glob al CMT c atalo g ( http: // www. globalcmt. org/ CMTsearch. html ). One c onspicuous fe atur e is that b-values ar e lar ge in o c e anic ridge zones, but smal l in plate sub duction zones. ically , ho w ev er, more eﬃcien t estimation is p er- formed b y the maximum likelihoo d metho d . Utsu ( 1965 ) derived it by the moment metho d. Later, Aki ( 1965 ) demonstrated that this is a maximum- lik eliho o d estimate (MLE) and provided the error estimate. It sh ould b e noted that the magnitud es in most catalog s are giv en in the int erv al of 0.1 (dis- crete magnitude v alues), hen ce, care should b e tak en for av oiding the bias of the b estimate in lik eliho o d - based estimation pro cedur es (Utsu ( 1969 )). Although the co eﬃcien t b in a wide r egion is generally sligh tly smaller th an 1.0, Guten b erg and Ric hter ( 1944 ) further determined that the b -v alue v aries according to lo cation in sm aller seismic re- gions. The b -v alue v aries ev en within Japan and fur- ther v aries with time. T emp oral and s p atial b -v alue c h anges hav e attracted the atten tion of man y re- searc h ers ev er since Suyehiro ( 1966 ) rep orted a d if- ference b etw een b -v alues of foreshocks and after- sho c ks in a sequence. Here, we consider that β can v ary with time, space and space–time according to a function suc h as β ( t ), β ( x, y , z ) or β ( t, x, y ) . V arious non p aramet- ric smo othing algorithms such as ke rnel m etho ds ha v e b een prop osed (Wiemer and Wyss, 1997 ). Al- ternativ ely , the β v alue can b e parameterized by smo oth cubic sp lines (Ogata, Imoto and K atsur a, 1991 ; Ogata and K atsura ( 1993 )) or piece-wise lin- ear expansions on Delauna y tessellated space (Ogata ( 2011c ); see also Figure 1 , e.g.). In suc h a case, a 4 Y. OGA T A p enalized log-lik eliho o d (Go o d and Gaskins ( 1971 )) is used whereb y th e log-lik eliho o d function is asso- ciated with p enalt y fun ctions in whic h the co eﬃ- cien ts are constrained for smo othness of the β func- tion. F or the optimal estimation of the β fun ction, the wei ght s in the p enalt y function are ob j ectiv ely adjusted in a Ba ye sian f ramew ork, as su ggested by Ak aik e ( 1980a ). 2.2.2 Af tersho ck analysis and pr ob abilistic for e- c asting T ypical aftersho c k frequen cy deca ys accord- ing to the rev erse p ow er fu nction w ith time (Omori ( 1894 ); Utsu , 1961 , 1969 ). First, let N ( s , t ) b e the n umb er of aftersho c ks in an in terv al ( s, t ). Then, the o ccurrence r ate of aftersho cks at th e elapsed time t since the main sho c k is ν ( t ) = lim ∆ → 0 1 ∆ P { N ( t, t + d ∆) ≥ 1 | Mainsho c k at time 0 } (2) = K ( t + c ) p for constan ts K , c and p . This is kn own as the Omori–Utsu (O–U) la w . T raditionally , estimates of the p arameter p ha v e b een obtained since the study of Utsu ( 1961 ) in the follo wing m an n er. The n umb ers of aftersho cks in a unit time in terv al n ( t ) are ﬁ rst plotted against elapsed time on doubly logarithmic axes, and then are ﬁt to an asymp totic straigh t line. The slop e of this line is an estimate for p . Th e v alues of c can b e determin ed by measuring the b en d ing curv e im- mediately after the m ain sho c k. Such an analysis is based on th e time series of counte d num b ers of af- tersho c ks. By such a plot, we can ﬁnd aftersho ck sequences for w hic h the formula ( 2 ) lasts a long p e- rio d, more than 120 y ears, for example (Utsu ( 1969 ); Ogata ( 1989 ); Utsu, O gata and Matsu’ur a, 1995 ). T o eﬃcien tly estimate th e three parameters di- rectly on the basis of o ccurrence time r ecords of aftersho c ks, assuming nonstationary P oisson pro- cess with intensit y function ( 2 ), Ogata ( 1983 ) su g- gested the maximum-lik eliho o d metho d, wh ic h en- abled th e practical aftersho ck forecasting. Reasen- b erg and Jones ( 1989 ) pr op osed a pro cedure based on the join t in tensit y rate of time and magnitud e of aftersho c ks (Utsu ( 1970 )) according to the G–R la w ( 1 ), λ ( t, M ) = λ 0 ( M ) ν ( t ) (3) = 10 a + b ( M 0 − M ) ( t + c ) p ( a, b, c, p ; constant) , where M is the magnitude of an aftersho c k and t is the time follo wing a main sho c k of magnitude M 0 ; the parameters are indep enden tly estimated by the maxim um-lik eliho o d metho d for resp ectiv e empir i- cal la ws. After a large earthquake o ccurs, the Japan Me- teorolog ical Agency (JMA) and the United States Geolog ical S urve y (USGS) hav e under taken op era- tional probabilit y forecast of the aftersho cks. Ho w- ev er, the forecast is announ ced after the elapse of 24 h or more. T his is d ue to the d eﬁciency of after- sho c k data due to o verlapping of seismograms after the main sho ck. In particular, the p arameter a is crucial f or the early forecast, but diﬃcu lt to estimate in an early p erio d , whereas the other parameters can b e default v alues for the early forecast [Reasen b erg and Jones ( 1994 ); Earthqu ak e Researc h Committee (ER C), 1998 ]. T he diﬃ cu lt y is b ecause the parame- ter a can substan tially d iﬀer ev en if the magnitudes of th e main sho c ks are almost the same: for exam- ple, the num b ers of the aftersh o c ks of M ≥ 4 . 0 of t w o nearby main sho c ks of the same M6.8 diﬀer b y 6–7 times (JMA ( 2009 )). It is notable that th e strongest aftersho c ks o c- curred within 24 h in most sequences (JMA ( 2008 )). Therefore, d esp ite adv erse conditions du ring data collect ion, pr obabilistic aftersho ck f orecasts sh ould b e d eliv ered as so on as p ossible within 24 h after the main sho c k to mitigate secondary disasters in aﬀected areas. F or this p urp ose, it is necessary to estimate time- dep end en t m issing rates, or detection rates, of af- tersho c ks (Ogata and Katsura, 1993 , 2006 ; Ogata, 2005c ) b ecause they enable probabilistic forecasting immediately after the main sho c k (Oga ta, 2005c ; Ogata and Katsura ( 2006 )). The detection rate of earthquak es is describ ed by a probabilit y function q ( M ) of magnitude M such th at 0 ≤ q ( M ) ≤ 1 . The in tensit y λ ( M ) for actually observed magnitude fre- quency is describ ed b y λ ( M ) = λ 0 ( M ) q ( M ), corre- sp ond ing to th in ning or random deletion. An ex- ample of the detection rate fun ction is the cum u- lativ e of Gauss ian distribu tion or the so-called er- ror function q ( M ) = erf { M | µ, σ } . The parameter µ represents the magnitude at which earthquak es are d etected at a rate of 50%, and σ repr esen ts a range of magnitudes in whic h earthquakes are partially detected. Let a data set of magnitudes { ( t i , M i ); i = 1 , . . . , N } b e giv en at a p erio d imme- diately after the main sho ck. Assum e that the pa- rameters are time-dep en den t d u ring the p erio d suc h EAR THQUAKE PREDICTION 5 that λ ( t, M ) = 10 a + b ( M 0 − M ) ( t + c ) p q { M | µ ( t ) , σ } (4) with an impro ving detection r ate µ ( t ). An additional parametric approac h prop osed by Omi et al. ( 2013 ) uses the state–space repr esen tation metho d for r eal- time forecasting within the 24 h p erio d. 2.2.3 Epidemic-typ e aftersho ck se quenc e (ET AS) mo del The epidemic-t yp e aftersho ck sequence (ET AS) mo del describ es earthquake activit y as a p oint pro cess (Ogata ( 1986 ), 1988 ) and in clud es th e O–U la w for aftersh o c ks as a d escendan t p ro cess. This mo d el assumes that the bac kgrou n d seismicit y is a stationary Po isson pro cess with a constant oc- currence r ate or num b er of earthquake s p er da y , µ . The conditional in tensit y function of the pr o cess is describ ed by λ θ ( t | H t ) = µ + X { i : t i 0 , resp ectiv ely , w h ere σ 1 = 0 . 6709 an d σ 2 = 0 . 4456 (km). Supp ose that, at the curren t time, c | n s ho ws the stage where the n th earthqu ake ( n = 2 , 3 , 4 , . . . , # c ) has o ccurr ed in a cluster c , where # c is the num b er of all earthquak es in the cluster c . W e prop ose the forecasting probabilit y p c | n b y usin g the follo w ing logistic transformation: Set f = logit p = (1 − p ) /p , or p = 1 / (1 + e f ); th en, logit p c | n = logit µ ( x 1 , y 1 ) + 1 #( i < j ≤ n ) (23) × X i

Original Paper

Loading high-quality paper...

Comments & Academic Discussion

Loading comments...

Original Paper

Related Papers

Comments & Academic Discussion

Leave a Comment