- Email: [email protected]
Transition matrix analysis of earthquake magnitude sequences
Chaos, Solitons and Fractals 24 (2005) 33–43 www.elsevier.com/locate/chaos
Transition matrix analysis of earthquake magnitude sequences Michele Lovallo, Vincenzo Lapenna, Luciano Telesca
*
Institute of Methodologies for Environmental Analysis (IMAA-CNR), National Research Council, C.da S.Loja, I-85050, Tito Scalo (PZ), Italy Accepted 28 July 2004
Abstract Estimation of complexity is a fascinating research topic in nonlinear signal and system analysis. Information theoretic functionals can be used to identify and quantify general relationships among variables; these relationships can be considered as the fingerprints of complexity. Up to now, the complexity of seismic sequences has been mostly related to the concept of self-similarity, suggesting that the earthquake dynamics can be interpreted as due to many components interacting over a wide range of time or space scales. This paper deals with a new idea of complexity of seismicity, focusing, in particular, on the transition probability between magnitudes. Using the Transition Matrix Method, a set of complexity parameters can be defined for earthquakes. Furthermore, the relationships among these parameters and those characterizing the earthquake magnitude dynamics have been analyzed in simulated and observational seismic sequences. Ó 2004 Elsevier Ltd. All rights reserved.
1. Introduction Qualifying and quantifying complexity of a time series has fascinated researchers for years. The interests in developing this idea have been fueled by addressing to issues of dynamical systems theory, statistics or formal mathematics and computer science [1], despite the fact that complexity itself is only value defined, and many alternatives have been proposed over the years [2–5]. Research in the directions suggested by the different scientific fields has led to several definitions of complexity, mainly in the realm of computational sciences. Among these, several can be cited: the Kolmogorov–Chaitin [6] algorithmic complexity, the Lempel–Ziv complexity [7], the logical depth Bennet [8], the effective measure complexity of Grassberger [9], the complexity of a system based on its diversity [10], the thermodynamic depth [11], etc. The common objective of all these definitions is to be operational in the context in which they are used, but to be general enough to be connected with our intuitive notion of complexity [1]. Information theoretic functionals such as the information entropy can be used to analyze the behaviour of nonlinear dynamical systems and time series [12–16]; furthermore they can be used to identify and quantify the complexity in time series. Information entropy is also used to reveal complexity in symbol sequences such as DNA sequences [17], and has represented a major advance in the description and comprehension of a wide range of phenomena, from geophysics to industrial processes [18,19].
*
Corresponding author. Tel.: +39 0971 427206; fax: +39 0971 427271. E-mail address: [email protected] (L. Telesca).
0960-0779/$ - see front matter Ó 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.chaos.2004.07.024
34
M. Lovallo et al. / Chaos, Solitons and Fractals 24 (2005) 33–43 6
b2>b1 5
b2 Log10N(M>ML)
4
b1
3
2
MC 1
0 -1
0
1
2
3
4
5
6
7
ML
Fig. 1. The Gutenberg–Richter law: the exponent b reveals the capability of an area to generate large earthquakes; lower b-value higher the probability of a large earthquake to occur. The completeness magnitude Mc can be estimated as the lowest magnitude above which the log-linear relationship is held.
In the study of complex systems, symbolic sequences play an important role, and most systems whose complexity we would like to quantify can be reduced to them. When analyzing a symbolic sequence the main issue is to extract the information it brings. In DNA sequences this would lead to the identification of sub-sequences codifying the genes and their specific functions; for a written text this would correspond to recognizing the language in which the text is written; for time series to extracting specific features and trends [5]. An earthquake magnitude sequence can be considered as a symbolic sequence. Each symbol of the sequence is a real value, the magnitude of the earthquake mi, extracted from the set [ML, . . ., MH] where ML is the lowest magnitude of the sequence, consistent with the completeness of the seismic catalogue, and MH is the highest magnitude of the same sequence. In seismology, it is well established the Gutenberg–Richter law [20] (GR), log10 ðN Þ ¼ a bM L , where N is the number of events whose magnitude M is equal to a threshold magnitude ML, and b is the exponent which informs about the productivity rate of the events. The constraint above ML is ML P Mc, where Mc is called the completeness magnitude, that it is the lowest magnitude, above which the seismic catalogue can be considered complete; the completeness of the catalogue is linked to the capability of the seismic network to recognize all the earthquakes. Mc can be estimated as the lowest magnitude under which the log-linear GR relation does not exist. Fig. 1 explains two cases of seismic areas, characterized by two different b-values: lower b-value, higher the probability to an earthquake with large magnitude to be generated. This law describes only the magnitude frequency of earthquakes occurring in a seismic area; but it reveals none about the properties of transition between one magnitude and another. The aim of this paper is the complexity of sequences of earthquakes, approached by means of the information entropy calculated by the transition matrix method applied to the sequence of magnitudes. This method has been used in the analysis of complexity of sequences of DNA [21], showing that some measures of complexity based on transition proportion matrix are of interest.
2. Transition matrix method The concept of transition matrix of a data sequence is well explained in the book of Davis [22]. Here we use this method to study sequences of earthquake magnitudes. For a given magnitude sequence m = m1m2 . . . mN we can construct the matrix A = (Tij), where Tij means the number of times a given magnitude being succeeded by another in the sequence. The dimension of the matrix is given by (MH ML) Æ 10 + 1, and it represents the number of different magnitudes ranging from ML to MH, if the step between one magnitude to another is 0.1. ML is the lowest magnitude of the sequence, consistent with the completeness of the seismic catalogue, while MH is the highest magnitude of the same sequence. A is called the transition frequency matrix of m, which is a concise way of expressing the incidence of one magnitude following another.
M. Lovallo et al. / Chaos, Solitons and Fractals 24 (2005) 33–43
35
For example, the transition frequency matrix for the magnitude sequence m ¼ 2:5; 2:5; 2:7; 2:5; 2:6; 2:8; 2:6; 2:7; 2:8; 2:8; 2:5; 2:5; 2:7; 2:6; 2:5; 2:8; 2:6; 2:5; 2:5 is the following: to
from
2:5 2:6
2:7
2:8
2:5
3
1
2
1
2:6
2
0
1
1
2:7
1
1
0
1
2:8
1
2
0
1
The tendency for one magnitude to succeed another can be emphasized by converting the frequency matrix to decimal fractions or percentages. Therefore, we can construct a matrix P = (Pij) by dividing each element by the sum of all entries in A. Such a matrix represents the relative frequency of all the possible types of transitions, and is called the transition proportion matrix of m. For the above example, the transition proportion matrix is To 2:5 from
2:6
2:7
2:8
2:5
0:17 0:06
0:11 0:06
2:6
0:11 0
0:06 0:06
2:7
0:06 0:06
0
0:06
2:8
0:06 0:11
0
0:06
P Since Ni;j¼1 P ij ¼ 1; 0 6 P ij 6 1, we can view Pij as the probability of one magnitude to succeed another. Let us now discuss a measure of complexity based on the statistical description of systems. In the case of transition of one magnitude to another, the system has different accessible states of transition. Each state is indicated by the corresponding probability Pij. It is possible to find a quantity measuring the amount of ‘‘information’’. If we denote #{Pij : Pij 5 0} = K as the number of probabilities which is not zero, and rewrite {Pij : Pij 5 0} as fP i gKi¼1 , then Shannons definition [23] of information entropy applies H ¼
K X
P i lnðP i Þ:
ð1Þ
i¼1
The Shannon entropy H measures the information needed to locate a system in a certain state, meaning that H is a measure of our ignorance about the system [24]. In our case the information entropy measures the average amount of information per magnitude transition, and, as nicely depicted in Puglisi et al. [25], it is an estimate of the ‘‘surprise’’ the source producing the sequence reserves to us. In which sense the entropy can be viewed as a measure of ‘‘surprise’’? Let us suppose that the ‘‘surprise’’ one feels in observing an event E occurring depends only on probability of E. If the event occurs with probability 1, our ‘‘surprise’’ is zero. If the event occurs with a very low probability, our ‘‘surprise’’ in observing it is proportionally large. Suppose the event E is a particular transition from one magnitude to another; for a single transition occurring with a probability Pi, the ‘‘surprise’’ is proportional to ln(Pi). For the set of K transitions, the expected amount of ‘‘surprise’’ is given by the formula (1). When Pi = 1/K, for i = 1, 2, . . ., K, i.e. the case of equilibrium state, the entropy H reaches its maximum value. When Pi = 1 for some i and Pj = 0 for j 5 i, we have H = 0. There is also a definition of disequilibrium [26] D, used as a measure of complexity in a K-system. The intuitive notion suggests that some kind of distance from an equiprobable distribution should be used. Two constraints have to be imposed on D : D > 0, in order that the complexity is positive, and D = 0 for equiprobable systems. The definition follows: 2 K X 1 D¼ Pi : ð2Þ K i¼1
36
M. Lovallo et al. / Chaos, Solitons and Fractals 24 (2005) 33–43
In the case of equilibrium D = 0. When Pi = 1 for some i and Pj = 0 for j 5 i, we have D assuming its maximum value. Lope-Ruiz et al. [27] proposed another statistical measure of complexity C, which is simply the interplay between the information stored in the system and its disequilibrium; it is defined as C ¼ H D:
ð3Þ
Now C = 0 for both the equilibrium state and the case of Pi = 1 for some i and Pj = 0 for j 5 i.
3. Results In order to investigate the relationship between seismic parameters (GR b-value and earthquake magnitude) and complexity parameters, we performed a numerical study on randomly generated seismic sequences; then we analyzed such relationship on observational earthquake series. We simulated earthquake magnitude sequences, in accordance with the GR law [20], N 10bM L , where N is the number of events with magnitude equal or larger than ML and b is the exponent of the power law, which, in observational cases, assumes values ranging from 0.5 to 2.0, depending on the seismic area. The GR law is an estimate of the probability density function of the earthquake magnitudes, furnishing their distribution frequency. The simulated earthquake magnitude sequences have been generated by means of the rejection method, which is a powerful general technique for generating random deviates. In our case they have distribution function p(ML) N/NTOT, where NTOT is the total number of events in the sequence, whose magnitude is larger or equal to Mc, being Mc the completeness magnitude. Fig. 2 presents the complexity parameters for earthquake magnitude sequences, generated with a completeness magnitude Mc = 2.0 and varying both the b-value of the GR law and the number of events N. Each plot shows the average between 10 runs. The error bars in each plot delimit the 1 r range around the average value. Fig. 2a shows the information entropy H as a function of the number of events N for different values of the GR b-value (from 0.5 to 2.0). The relationship between the entropy H and the number of events N is clearly nonlinear for all b-values. A small number of events N limits the possible transitions among the earthquake magnitudes; increasing the number of events N, the system can be found in any of its accessible transitions with the same probability. Therefore the information content is low for small N, while increases with N. The asymptotic values of the H reveal that the maximum of information decreases with the b-value; this indicates that earthquake sequences, with smaller b-value and typical of seismic areas characterized by high probability to generate larger shocks, convey larger amount of information. Fig. 2b shows the disequilibrium D. For all the considered b-values the disequilibrium decays with the number of events N, approaching to an asymptotic value for large N. It assumes relative high values for small N (there are few dominant states of transition among magnitudes), while D decreases with the increase of N (the number of equiprobable transitions increases). The variation with the b-value reveals that D increases with bvalue, indicating that the most equiprobable distribution of the transitions is assumed for low b-values. Fig. 2c shows the complexity C. This parameter behaves similarly to the disequilibrium D, increasing with the b-value and decreasing with the number of events N. We investigated the variation of the complexity parameters H, D and C with the threshold magnitude ML, that is the lowest magnitude of the sequence of earthquakes that follows a GR distribution. Fig. 3 shows the results for values of b = 0.5, 1.0, 1.5 and 2.0 respectively and N = 105 events, with varying the lowest magnitude ML P Mc, with Mc = 2.0. Since the increase of ML is equivalent to a decrease of the effective number of magnitudes, also the number of accessible transitions decreases, implying the decrease of H and the increase of D and C, consistently with the results shown in Fig. 2. In all examined cases we observe that the complexity parameters are almost constant for small threshold magnitude ML up to certain threshold; we can estimate that this threshold varies from approximately 4.0 (b = 0.5) to 3.2 (b = 2.0). Furthermore, lower the b-value slower the variation. We have investigated three different observational seismic catalogues: (i) shallow (depth 6 60 km) Italian seismicity from 1983 to 2002; (ii) shallow (depth 6 60 km) Italian aftershock-depleted seismicity from 1983 to 2002, obtained from the full catalogue by applying the method of Reasenberg [28]; (iii) shallow (depth 6 25 km) southern California seismicity. The Gutenberg–Richter analysis has yielded Mc = 2.5 and b 1.4 for the full Italian catalogue, Mc = 2.5 and b 1.5 for the aftershock-depleted Italian catalogue and Mc = 1.5 and b 0.94 for the southern California catalogue. Figs. 4–6 show the results for the three catalogues respectively. To verify the results we performed an analysis on surrogate sequences. Figs. 4–6 also show the averages among the values obtained from 10 randomly generated sequences, with the same number of events N and the same b-value as the observational ones. The good accordance between the complexity parameters calculated for the original sequences and those obtained from the surrogates is clearly visible. Fig. 7 shows the comparison between the complexity parameters calculated for full and depleted Italian catalogues, showing lesser entropy information and larger disequilibrium and complexity for the aftershock-depleted catalogue.
M. Lovallo et al. / Chaos, Solitons and Fractals 24 (2005) 33–43
37
b=0.5
6
Entropy
5
b=1.0
4
b=1.5
b=2.0 3 0.0
4
4
4
4.0x10
2.0x10
4
6.0x10
(a)
5
80 . x10
1.0x10
N
0.07
0.06
b=2.0 Disequilibrium
0.05
0.04
b=1.5
0.03
0.02
b=1.0 0.01
b=0.5 0.00 0.0
4
2.0x10
4
4
4.0x10
(b)
6.0x10
4
8 .0x10
5
1.0x10
N
0.24
0.20
b=2.0
Complexity
0.16
b=1.5 0.12
0.08
b=1.0
0.04
b=0.5 0.0
(c)
4
2.0x10
4
4
4.0x10
6.0x10
4
8 .0x10
5
1.0x10
N
Fig. 2. Variation of the complexity parameters with the number of events N and the b-value.
4. Discussion In many previous seismological studies, the complex behaviour of seismic sequences has been often linked to the concept of self-similarity, suggesting that their dynamics may be interpreted as due to many components
38
M. Lovallo et al. / Chaos, Solitons and Fractals 24 (2005) 33–43
N=10
b=0.5 b=1.0 b=1.5 b=2.0
5
6.5 6.0 5.5 5.0
Entropy
4.5 4.0 3.5 3.0 2.5 2.0 1.5 1.0 2
3
(a)
4
5
6
ML
N=10
b=0.5 b=1.0 b=1.5 b=2.0
5
0.30 0.28 0.26 0.24 0.22
Disequilibrium
0.20 0.18 0.16 0.14 0.12 0.10 0.08 0.06 0.04 0.02 0.00 2
3
(b) 0.60
4
5
6
ML
N=10
b=0.5 b=1.0 b=1.5 b=2.0
5
0.55 0.50 0.45
Complexity
0.40 0.35 0.30 0.25 0.20 0.15 0.10 0.05 0.00 2
(c)
3
4
5
6
ML
Fig. 3. Variation of the complexity parameters with the b-value and threshold magnitude ML (N = 105).
interacting over a wide range of time or space scales [30]. In this case, fractal tools are the most adequate methods to investigate and characterize these complex phenomena. In particular, tectonic processes have been generally considered to display power-law properties in the space-time domain [31–33], and the complexity of the earthquake
M. Lovallo et al. / Chaos, Solitons and Fractals 24 (2005) 33–43
39
original simulated
Italy Full catalogue 5
Entropy
4
3
2
1 2.5
3.0
3.5
(a)
0.30 0.28
4.0
4.5
5.0
ML
original simulated
Italy Full catalogue
0.26 0.24 0.22
Disequilibrium
0.20 0.18 0.16 0.14 0.12 0.10 0.08 0.06 0.04 0.02 0.00 2.5
3.0
3.5
(b)
0.7
4.0
4.5
5.0
ML original simulated
Italy Full catalogue
0.6 0.5
Complexity
0.4 0.3 0.2 0.1 0.0 -0.1 2.5
(c)
3.0
3.5
4.0
4.5
5.0
ML
Fig. 4. Complexity parameters for the full Italian catalogue.
dynamics has been revealed by the presence of space-time clustering behaviour at both short and long scales [34–36]. In this context, the magnitude distribution of an earthquake sequence has been only used to perform statistical analysis concerning the estimate of probability of occurrence of large events [29]. But, can an earthquake magnitude sequence convey additive information, concerning the complex behaviour of seismicity? An earthquake magnitude sequence is well described by the GR power-law [20], which is characterized by several parameters: (i) the b-value, which relates quantitatively to the probability that an earthquake with large magnitude
M. Lovallo et al. / Chaos, Solitons and Fractals 24 (2005) 33–43
5.4 5.2 5.0 4.8 4.6 4.4 4.2 4.0 3.8 3.6 3.4 3.2 3.0 2.8 2.6 2.4 2.2 2.0
Entropy
40
original simulated
Italy Depleted catalogue
2.5
3.0
(a)
3.5
4.0
4.5
ML original simulated
Italy Depleted catalogue
Disequilibrium
0.15
0.10
0.05
0.00 2.5
3.0
(b)
3.5
4.0
4.5
ML
0.50
original simulated
Italy Depleted catalogue
0.45 0.40
Complexity
0.35 0.30 0.25 0.20 0.15 0.10 0.05 2.5
(b)
3.0
3.5
4.0
4.5
ML
Fig. 5. Complexity parameters for the aftershock-depleted Italian catalogue.
can be generated; (ii) the completeness magnitude Mc, which indicates the lowest magnitude detectable by the seismic network; (iii) the threshold magnitude ML, which is directly involved in the GR power-law. Up to now these parameters have been exclusively used as statistical indicators of productivity of earthquakes of such magnitude and involved in seismic hazard evaluations [29]. In the analysis performed in the present study it is highlighted that the magnitude distribution of a seismic sequence can be informative of its complexity, which is defined in terms of the transition between one magnitude and another. Fundamentally, the complexity is related to the existence of equiprobable transitions from one magnitude to another.
M. Lovallo et al. / Chaos, Solitons and Fractals 24 (2005) 33–43 original simulated
Southern California
6.0
41
5.5
5.0
Entropy
4.5
4.0
3.5
3.0 1.0
1.5
2.0
2.5
3.0
(a)
3.5
4.0
4.5
5.0
5.5
ML original simulated
Southern California 0.050 0.045
Disequilibrium
0.040 0.035 0.030 0.025 0.020 0.015 0.010 0.005 1.0
1.5
2.0
2.5
3.0
(b)
3.5
4.0
4.5
5.0
5.5
ML original simulated
Southern California
0.16
0.14
Complexity
0.12
0.10
0.08
0.06
0.04 1.0
(c)
1.5
2.0
2.5
3.0
3.5
4.0
4.5
5.0
5.5
ML
Fig. 6. Complexity parameters for the Southern California catalogue.
The investigations performed on numerical and observational cases show that seismic areas ‘‘more capable’’ to generate large events (smaller b-value) are less complex than those characterized by higher frequency of small events (large b-value); larger b-value, less equiprobable the transitions between one magnitude and another. Furthermore, the complexity increases with the threshold magnitude, indicating a loss of equiprobability of the magnitude transitions with the increase of the threshold; this effect can be related to the number of events in the sequence, which becomes smaller augmenting the value of the threshold magnitude, thus limiting the equiprobability of the possible transitions among the earthquake magnitudes.
42
M. Lovallo et al. / Chaos, Solitons and Fractals 24 (2005) 33–43
5.0
depleted full
Italy
4.8 4.6 4.4 4.2 4.0
Entropy
3.8 3.6 3.4 3.2 3.0 2.8 2.6 2.4 2.2 2.0 2.5
3.0
3.5
(a) 0.12
4.0
4.5
5.0
ML depleted full
Italy
0.11 0.10 0.09
Disequilibrium
0.08 0.07 0.06 0.05 0.04 0.03 0.02 0.01 2.5
3.0
3.5
(b)
4.0
4.5
5.0
ML full depleted
Italy 0.26 0.24 0.22
Complexity
0.20 0.18 0.16 0.14 0.12 0.10 0.08 0.06 2.5
(c)
3.0
3.5
4.0
4.5
5.0
ML
Fig. 7. Comparison between the complexity parameters shown in Figs. 6 and 7.
5. Conclusions The application of information theoretical measures to the transition proportion matrix of the magnitudes of seismic sequences has revealed the following issues: (i) the ‘‘complexity’’ of a seismic sequence can be dealt focusing on the properties of transition from one magnitude to another; (ii) the ‘‘complexity’’ of a seismic sequence can be measured
M. Lovallo et al. / Chaos, Solitons and Fractals 24 (2005) 33–43
43
by means of three quantities, the information entropy H, the disequilibrium D and the complexity C; (iii) these three parameters show a nontrivial relationship with the length of the seismic sequence, with the b-value of the Gutenberg– Richter law, which describes the frequency of earthquake magnitude in a seismic area, and with the threshold magnitude ML; (iv) aftershocks cause increase of entropy information and decrease of disequilibrium and complexity in a seismic sequence. In the future, the time variation of complexity parameters will be performed in order to reveal the particular seismic patterns.
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36]
Bialek W, Nemenman I, Tishby N. Complexity through Nonextensivity.. Kolmogorov AN. Three approaches to the definition of the concept quantity of information. Probl Inform Trans 1965;1:1–7. Lloyd S, Pagels H. Complexity as thermodynamic depth. Ann Phys 1988;188:186–213. Crutchfield JP, Young K. Inferring statistical complexity. Phys Rev Lett 1989;188:105–8. Badii R, Politi A. Complexity-hierarchical structures and scaling in physics. Cambridge: Cambridge University Press; 1997. Kolmogorov AN. Three approaches to the quantitative definitions of information. Probl Inform Trans 1964;1(1):1–7; Chaitin GJ. Information, randomness and incompleteness. Singapore: World Scientific; 1990. Ziv J, Lempel A. On the complexity of finite sequences. IEEE Trans Inform Theory 1976;22:75–81. Bennet CH. Information, dissipation, and the definition of organization. In: Pines D, editor. Emerging syntheses in science. Reading, MA: Addison-Wesley; 1988. Grassberger P. Towards a quantitative theory of self-generated complexity. Int J Theor Phys 1986;25:907–38. Huberman BA, Hogg T. Complexity and adaptation. Physica D 1986;22:376–84. Lloyd S, Pagels H. Complexity as thermodynamic depth. Ann Phys (NY) 1988;188:186–213. Fraser AM. Information and entropy in strange attractors. IEEE Trans Inform Theory 1989;35:245–62. Palus M. Information theoretic test for nonlinearity in time-series. Phys Lett A 1993;175:203–9. Palus M. Testing for nonlinearity in weather records. Phys Lett A 1994;193:67–74. Palus M. Detecting nonlinearity in multivariate time-series. Phys Lett A 1996;213:138–47. Roulston MS. Estimating the errors on measured entropy and mutual information. Physica D 1999;125:285–94. Herzel H, Ebeling W, Schmitt A. Entropies of biosequences: the role of protein coding segments. Phys Rev E 1997;55:800–10. Kantz H, Schreiber T. In: Nonlinear time series analysis. Cambridge Nonlinear Science Series. Cambridge: Cambridge University Press; 1997. Abarbanel HDI. Analysis of observed chaotic data. New York: Springer; 1996. Gutenberg B, Richter CF. Frequency of earthquakes in California. Bull Seism Soc Am 1944;34:185–8. Yu Z-G, Chen G-Y. Rescaled range and transition matrix analysis of DNA sequences. Commun Theor Phys 2000;33(4):673–8. Davies JC. Statistics and data analysis in geology. New York, London, Sydney, Toronto: John & Sons, INC; 1973. Shannon CE, Weaver W. The mathematical theory of communication. Urbana, IL: University of Illinois Press; 1949. Torres ME, Gamero LG. Relative complexity changes in time series using information measures. Physica A 2000;286:457–73. Puglisi A, Bendetto D, Cagliotti E, Loreto V, Vulpiani A. Data compression and learning in time sequences analysis. Physica D 2003;180:92–107. Colis GNiE, Prigogine I. Self-organisation in nonequilibrium systems. New York: Wiley; 1977. Lopez-Ruiz R, Mancini HL, Calbet X. A statistical measure of complexity. Phys Lett A 1995;209:321–6. Reasenberg P. Second-order moment of central California seismicity, 1969–1982. J Geophys Res 1985;90:5479–95. Martinelli G, Albarello D. Main constraints for siting monitoring networks devoted to the study of earthquake related hydrogeochemical phenomena in Italy. Ann Geof 1997;XL:1505–25. Stanley HE, Amaral LAN, Goldberger AL, Havlin S, Ivanov PCh, Peng C-K. Statistical physics and physiology: monofractal and multifractal approaches. Physica A 1999;270:309–24. Godano C, Caruso V. Multifractal analysis of earthquake catalogues. Geophys J Int 1995;121:385–92. ¨ ncel AO, Main IG, Alptekin O ¨ , Cowie PA. Spatial variation in the fractal properties of seismicity in the north Anatolian fault O zone. Tectonophysics 1996;257:189–202. Wilson TH, Dominic J. Fractal interrelationships between topography and structure. Earth Surf Process Land 1998;23:509–25. Kagan YY, Jackson DD. Seismic gap hypothesis: ten years after. J Geophys Res 1991;96:21419–31. Telesca L, Cuomo V, Lapenna V, Macchiato M. Identifying space-time clustering properties of the 1983–1997 Irpinia-Basilicata (southern Italy) seismicity. Tectonophysics 2001;330:93–102. Telesca L, Cuomo V, Lapenna V, Macchiato M. Depth-dependent time-clustering behaviour in seismicity of southern California. Geophys Res Lett 2001;28:4323–6.
Transition matrix analysis of earthquake magnitude sequences Michele Lovallo, Vincenzo Lapenna, Luciano Telesca
*
Institute of Methodologies for Environmental Analysis (IMAA-CNR), National Research Council, C.da S.Loja, I-85050, Tito Scalo (PZ), Italy Accepted 28 July 2004
Abstract Estimation of complexity is a fascinating research topic in nonlinear signal and system analysis. Information theoretic functionals can be used to identify and quantify general relationships among variables; these relationships can be considered as the fingerprints of complexity. Up to now, the complexity of seismic sequences has been mostly related to the concept of self-similarity, suggesting that the earthquake dynamics can be interpreted as due to many components interacting over a wide range of time or space scales. This paper deals with a new idea of complexity of seismicity, focusing, in particular, on the transition probability between magnitudes. Using the Transition Matrix Method, a set of complexity parameters can be defined for earthquakes. Furthermore, the relationships among these parameters and those characterizing the earthquake magnitude dynamics have been analyzed in simulated and observational seismic sequences. Ó 2004 Elsevier Ltd. All rights reserved.
1. Introduction Qualifying and quantifying complexity of a time series has fascinated researchers for years. The interests in developing this idea have been fueled by addressing to issues of dynamical systems theory, statistics or formal mathematics and computer science [1], despite the fact that complexity itself is only value defined, and many alternatives have been proposed over the years [2–5]. Research in the directions suggested by the different scientific fields has led to several definitions of complexity, mainly in the realm of computational sciences. Among these, several can be cited: the Kolmogorov–Chaitin [6] algorithmic complexity, the Lempel–Ziv complexity [7], the logical depth Bennet [8], the effective measure complexity of Grassberger [9], the complexity of a system based on its diversity [10], the thermodynamic depth [11], etc. The common objective of all these definitions is to be operational in the context in which they are used, but to be general enough to be connected with our intuitive notion of complexity [1]. Information theoretic functionals such as the information entropy can be used to analyze the behaviour of nonlinear dynamical systems and time series [12–16]; furthermore they can be used to identify and quantify the complexity in time series. Information entropy is also used to reveal complexity in symbol sequences such as DNA sequences [17], and has represented a major advance in the description and comprehension of a wide range of phenomena, from geophysics to industrial processes [18,19].
*
Corresponding author. Tel.: +39 0971 427206; fax: +39 0971 427271. E-mail address: [email protected] (L. Telesca).
0960-0779/$ - see front matter Ó 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.chaos.2004.07.024
34
M. Lovallo et al. / Chaos, Solitons and Fractals 24 (2005) 33–43 6
b2>b1 5
b2 Log10N(M>ML)
4
b1
3
2
MC 1
0 -1
0
1
2
3
4
5
6
7
ML
Fig. 1. The Gutenberg–Richter law: the exponent b reveals the capability of an area to generate large earthquakes; lower b-value higher the probability of a large earthquake to occur. The completeness magnitude Mc can be estimated as the lowest magnitude above which the log-linear relationship is held.
In the study of complex systems, symbolic sequences play an important role, and most systems whose complexity we would like to quantify can be reduced to them. When analyzing a symbolic sequence the main issue is to extract the information it brings. In DNA sequences this would lead to the identification of sub-sequences codifying the genes and their specific functions; for a written text this would correspond to recognizing the language in which the text is written; for time series to extracting specific features and trends [5]. An earthquake magnitude sequence can be considered as a symbolic sequence. Each symbol of the sequence is a real value, the magnitude of the earthquake mi, extracted from the set [ML, . . ., MH] where ML is the lowest magnitude of the sequence, consistent with the completeness of the seismic catalogue, and MH is the highest magnitude of the same sequence. In seismology, it is well established the Gutenberg–Richter law [20] (GR), log10 ðN Þ ¼ a bM L , where N is the number of events whose magnitude M is equal to a threshold magnitude ML, and b is the exponent which informs about the productivity rate of the events. The constraint above ML is ML P Mc, where Mc is called the completeness magnitude, that it is the lowest magnitude, above which the seismic catalogue can be considered complete; the completeness of the catalogue is linked to the capability of the seismic network to recognize all the earthquakes. Mc can be estimated as the lowest magnitude under which the log-linear GR relation does not exist. Fig. 1 explains two cases of seismic areas, characterized by two different b-values: lower b-value, higher the probability to an earthquake with large magnitude to be generated. This law describes only the magnitude frequency of earthquakes occurring in a seismic area; but it reveals none about the properties of transition between one magnitude and another. The aim of this paper is the complexity of sequences of earthquakes, approached by means of the information entropy calculated by the transition matrix method applied to the sequence of magnitudes. This method has been used in the analysis of complexity of sequences of DNA [21], showing that some measures of complexity based on transition proportion matrix are of interest.
2. Transition matrix method The concept of transition matrix of a data sequence is well explained in the book of Davis [22]. Here we use this method to study sequences of earthquake magnitudes. For a given magnitude sequence m = m1m2 . . . mN we can construct the matrix A = (Tij), where Tij means the number of times a given magnitude being succeeded by another in the sequence. The dimension of the matrix is given by (MH ML) Æ 10 + 1, and it represents the number of different magnitudes ranging from ML to MH, if the step between one magnitude to another is 0.1. ML is the lowest magnitude of the sequence, consistent with the completeness of the seismic catalogue, while MH is the highest magnitude of the same sequence. A is called the transition frequency matrix of m, which is a concise way of expressing the incidence of one magnitude following another.
M. Lovallo et al. / Chaos, Solitons and Fractals 24 (2005) 33–43
35
For example, the transition frequency matrix for the magnitude sequence m ¼ 2:5; 2:5; 2:7; 2:5; 2:6; 2:8; 2:6; 2:7; 2:8; 2:8; 2:5; 2:5; 2:7; 2:6; 2:5; 2:8; 2:6; 2:5; 2:5 is the following: to
from
2:5 2:6
2:7
2:8
2:5
3
1
2
1
2:6
2
0
1
1
2:7
1
1
0
1
2:8
1
2
0
1
The tendency for one magnitude to succeed another can be emphasized by converting the frequency matrix to decimal fractions or percentages. Therefore, we can construct a matrix P = (Pij) by dividing each element by the sum of all entries in A. Such a matrix represents the relative frequency of all the possible types of transitions, and is called the transition proportion matrix of m. For the above example, the transition proportion matrix is To 2:5 from
2:6
2:7
2:8
2:5
0:17 0:06
0:11 0:06
2:6
0:11 0
0:06 0:06
2:7
0:06 0:06
0
0:06
2:8
0:06 0:11
0
0:06
P Since Ni;j¼1 P ij ¼ 1; 0 6 P ij 6 1, we can view Pij as the probability of one magnitude to succeed another. Let us now discuss a measure of complexity based on the statistical description of systems. In the case of transition of one magnitude to another, the system has different accessible states of transition. Each state is indicated by the corresponding probability Pij. It is possible to find a quantity measuring the amount of ‘‘information’’. If we denote #{Pij : Pij 5 0} = K as the number of probabilities which is not zero, and rewrite {Pij : Pij 5 0} as fP i gKi¼1 , then Shannons definition [23] of information entropy applies H ¼
K X
P i lnðP i Þ:
ð1Þ
i¼1
The Shannon entropy H measures the information needed to locate a system in a certain state, meaning that H is a measure of our ignorance about the system [24]. In our case the information entropy measures the average amount of information per magnitude transition, and, as nicely depicted in Puglisi et al. [25], it is an estimate of the ‘‘surprise’’ the source producing the sequence reserves to us. In which sense the entropy can be viewed as a measure of ‘‘surprise’’? Let us suppose that the ‘‘surprise’’ one feels in observing an event E occurring depends only on probability of E. If the event occurs with probability 1, our ‘‘surprise’’ is zero. If the event occurs with a very low probability, our ‘‘surprise’’ in observing it is proportionally large. Suppose the event E is a particular transition from one magnitude to another; for a single transition occurring with a probability Pi, the ‘‘surprise’’ is proportional to ln(Pi). For the set of K transitions, the expected amount of ‘‘surprise’’ is given by the formula (1). When Pi = 1/K, for i = 1, 2, . . ., K, i.e. the case of equilibrium state, the entropy H reaches its maximum value. When Pi = 1 for some i and Pj = 0 for j 5 i, we have H = 0. There is also a definition of disequilibrium [26] D, used as a measure of complexity in a K-system. The intuitive notion suggests that some kind of distance from an equiprobable distribution should be used. Two constraints have to be imposed on D : D > 0, in order that the complexity is positive, and D = 0 for equiprobable systems. The definition follows: 2 K X 1 D¼ Pi : ð2Þ K i¼1
36
M. Lovallo et al. / Chaos, Solitons and Fractals 24 (2005) 33–43
In the case of equilibrium D = 0. When Pi = 1 for some i and Pj = 0 for j 5 i, we have D assuming its maximum value. Lope-Ruiz et al. [27] proposed another statistical measure of complexity C, which is simply the interplay between the information stored in the system and its disequilibrium; it is defined as C ¼ H D:
ð3Þ
Now C = 0 for both the equilibrium state and the case of Pi = 1 for some i and Pj = 0 for j 5 i.
3. Results In order to investigate the relationship between seismic parameters (GR b-value and earthquake magnitude) and complexity parameters, we performed a numerical study on randomly generated seismic sequences; then we analyzed such relationship on observational earthquake series. We simulated earthquake magnitude sequences, in accordance with the GR law [20], N 10bM L , where N is the number of events with magnitude equal or larger than ML and b is the exponent of the power law, which, in observational cases, assumes values ranging from 0.5 to 2.0, depending on the seismic area. The GR law is an estimate of the probability density function of the earthquake magnitudes, furnishing their distribution frequency. The simulated earthquake magnitude sequences have been generated by means of the rejection method, which is a powerful general technique for generating random deviates. In our case they have distribution function p(ML) N/NTOT, where NTOT is the total number of events in the sequence, whose magnitude is larger or equal to Mc, being Mc the completeness magnitude. Fig. 2 presents the complexity parameters for earthquake magnitude sequences, generated with a completeness magnitude Mc = 2.0 and varying both the b-value of the GR law and the number of events N. Each plot shows the average between 10 runs. The error bars in each plot delimit the 1 r range around the average value. Fig. 2a shows the information entropy H as a function of the number of events N for different values of the GR b-value (from 0.5 to 2.0). The relationship between the entropy H and the number of events N is clearly nonlinear for all b-values. A small number of events N limits the possible transitions among the earthquake magnitudes; increasing the number of events N, the system can be found in any of its accessible transitions with the same probability. Therefore the information content is low for small N, while increases with N. The asymptotic values of the H reveal that the maximum of information decreases with the b-value; this indicates that earthquake sequences, with smaller b-value and typical of seismic areas characterized by high probability to generate larger shocks, convey larger amount of information. Fig. 2b shows the disequilibrium D. For all the considered b-values the disequilibrium decays with the number of events N, approaching to an asymptotic value for large N. It assumes relative high values for small N (there are few dominant states of transition among magnitudes), while D decreases with the increase of N (the number of equiprobable transitions increases). The variation with the b-value reveals that D increases with bvalue, indicating that the most equiprobable distribution of the transitions is assumed for low b-values. Fig. 2c shows the complexity C. This parameter behaves similarly to the disequilibrium D, increasing with the b-value and decreasing with the number of events N. We investigated the variation of the complexity parameters H, D and C with the threshold magnitude ML, that is the lowest magnitude of the sequence of earthquakes that follows a GR distribution. Fig. 3 shows the results for values of b = 0.5, 1.0, 1.5 and 2.0 respectively and N = 105 events, with varying the lowest magnitude ML P Mc, with Mc = 2.0. Since the increase of ML is equivalent to a decrease of the effective number of magnitudes, also the number of accessible transitions decreases, implying the decrease of H and the increase of D and C, consistently with the results shown in Fig. 2. In all examined cases we observe that the complexity parameters are almost constant for small threshold magnitude ML up to certain threshold; we can estimate that this threshold varies from approximately 4.0 (b = 0.5) to 3.2 (b = 2.0). Furthermore, lower the b-value slower the variation. We have investigated three different observational seismic catalogues: (i) shallow (depth 6 60 km) Italian seismicity from 1983 to 2002; (ii) shallow (depth 6 60 km) Italian aftershock-depleted seismicity from 1983 to 2002, obtained from the full catalogue by applying the method of Reasenberg [28]; (iii) shallow (depth 6 25 km) southern California seismicity. The Gutenberg–Richter analysis has yielded Mc = 2.5 and b 1.4 for the full Italian catalogue, Mc = 2.5 and b 1.5 for the aftershock-depleted Italian catalogue and Mc = 1.5 and b 0.94 for the southern California catalogue. Figs. 4–6 show the results for the three catalogues respectively. To verify the results we performed an analysis on surrogate sequences. Figs. 4–6 also show the averages among the values obtained from 10 randomly generated sequences, with the same number of events N and the same b-value as the observational ones. The good accordance between the complexity parameters calculated for the original sequences and those obtained from the surrogates is clearly visible. Fig. 7 shows the comparison between the complexity parameters calculated for full and depleted Italian catalogues, showing lesser entropy information and larger disequilibrium and complexity for the aftershock-depleted catalogue.
M. Lovallo et al. / Chaos, Solitons and Fractals 24 (2005) 33–43
37
b=0.5
6
Entropy
5
b=1.0
4
b=1.5
b=2.0 3 0.0
4
4
4
4.0x10
2.0x10
4
6.0x10
(a)
5
80 . x10
1.0x10
N
0.07
0.06
b=2.0 Disequilibrium
0.05
0.04
b=1.5
0.03
0.02
b=1.0 0.01
b=0.5 0.00 0.0
4
2.0x10
4
4
4.0x10
(b)
6.0x10
4
8 .0x10
5
1.0x10
N
0.24
0.20
b=2.0
Complexity
0.16
b=1.5 0.12
0.08
b=1.0
0.04
b=0.5 0.0
(c)
4
2.0x10
4
4
4.0x10
6.0x10
4
8 .0x10
5
1.0x10
N
Fig. 2. Variation of the complexity parameters with the number of events N and the b-value.
4. Discussion In many previous seismological studies, the complex behaviour of seismic sequences has been often linked to the concept of self-similarity, suggesting that their dynamics may be interpreted as due to many components
38
M. Lovallo et al. / Chaos, Solitons and Fractals 24 (2005) 33–43
N=10
b=0.5 b=1.0 b=1.5 b=2.0
5
6.5 6.0 5.5 5.0
Entropy
4.5 4.0 3.5 3.0 2.5 2.0 1.5 1.0 2
3
(a)
4
5
6
ML
N=10
b=0.5 b=1.0 b=1.5 b=2.0
5
0.30 0.28 0.26 0.24 0.22
Disequilibrium
0.20 0.18 0.16 0.14 0.12 0.10 0.08 0.06 0.04 0.02 0.00 2
3
(b) 0.60
4
5
6
ML
N=10
b=0.5 b=1.0 b=1.5 b=2.0
5
0.55 0.50 0.45
Complexity
0.40 0.35 0.30 0.25 0.20 0.15 0.10 0.05 0.00 2
(c)
3
4
5
6
ML
Fig. 3. Variation of the complexity parameters with the b-value and threshold magnitude ML (N = 105).
interacting over a wide range of time or space scales [30]. In this case, fractal tools are the most adequate methods to investigate and characterize these complex phenomena. In particular, tectonic processes have been generally considered to display power-law properties in the space-time domain [31–33], and the complexity of the earthquake
M. Lovallo et al. / Chaos, Solitons and Fractals 24 (2005) 33–43
39
original simulated
Italy Full catalogue 5
Entropy
4
3
2
1 2.5
3.0
3.5
(a)
0.30 0.28
4.0
4.5
5.0
ML
original simulated
Italy Full catalogue
0.26 0.24 0.22
Disequilibrium
0.20 0.18 0.16 0.14 0.12 0.10 0.08 0.06 0.04 0.02 0.00 2.5
3.0
3.5
(b)
0.7
4.0
4.5
5.0
ML original simulated
Italy Full catalogue
0.6 0.5
Complexity
0.4 0.3 0.2 0.1 0.0 -0.1 2.5
(c)
3.0
3.5
4.0
4.5
5.0
ML
Fig. 4. Complexity parameters for the full Italian catalogue.
dynamics has been revealed by the presence of space-time clustering behaviour at both short and long scales [34–36]. In this context, the magnitude distribution of an earthquake sequence has been only used to perform statistical analysis concerning the estimate of probability of occurrence of large events [29]. But, can an earthquake magnitude sequence convey additive information, concerning the complex behaviour of seismicity? An earthquake magnitude sequence is well described by the GR power-law [20], which is characterized by several parameters: (i) the b-value, which relates quantitatively to the probability that an earthquake with large magnitude
M. Lovallo et al. / Chaos, Solitons and Fractals 24 (2005) 33–43
5.4 5.2 5.0 4.8 4.6 4.4 4.2 4.0 3.8 3.6 3.4 3.2 3.0 2.8 2.6 2.4 2.2 2.0
Entropy
40
original simulated
Italy Depleted catalogue
2.5
3.0
(a)
3.5
4.0
4.5
ML original simulated
Italy Depleted catalogue
Disequilibrium
0.15
0.10
0.05
0.00 2.5
3.0
(b)
3.5
4.0
4.5
ML
0.50
original simulated
Italy Depleted catalogue
0.45 0.40
Complexity
0.35 0.30 0.25 0.20 0.15 0.10 0.05 2.5
(b)
3.0
3.5
4.0
4.5
ML
Fig. 5. Complexity parameters for the aftershock-depleted Italian catalogue.
can be generated; (ii) the completeness magnitude Mc, which indicates the lowest magnitude detectable by the seismic network; (iii) the threshold magnitude ML, which is directly involved in the GR power-law. Up to now these parameters have been exclusively used as statistical indicators of productivity of earthquakes of such magnitude and involved in seismic hazard evaluations [29]. In the analysis performed in the present study it is highlighted that the magnitude distribution of a seismic sequence can be informative of its complexity, which is defined in terms of the transition between one magnitude and another. Fundamentally, the complexity is related to the existence of equiprobable transitions from one magnitude to another.
M. Lovallo et al. / Chaos, Solitons and Fractals 24 (2005) 33–43 original simulated
Southern California
6.0
41
5.5
5.0
Entropy
4.5
4.0
3.5
3.0 1.0
1.5
2.0
2.5
3.0
(a)
3.5
4.0
4.5
5.0
5.5
ML original simulated
Southern California 0.050 0.045
Disequilibrium
0.040 0.035 0.030 0.025 0.020 0.015 0.010 0.005 1.0
1.5
2.0
2.5
3.0
(b)
3.5
4.0
4.5
5.0
5.5
ML original simulated
Southern California
0.16
0.14
Complexity
0.12
0.10
0.08
0.06
0.04 1.0
(c)
1.5
2.0
2.5
3.0
3.5
4.0
4.5
5.0
5.5
ML
Fig. 6. Complexity parameters for the Southern California catalogue.
The investigations performed on numerical and observational cases show that seismic areas ‘‘more capable’’ to generate large events (smaller b-value) are less complex than those characterized by higher frequency of small events (large b-value); larger b-value, less equiprobable the transitions between one magnitude and another. Furthermore, the complexity increases with the threshold magnitude, indicating a loss of equiprobability of the magnitude transitions with the increase of the threshold; this effect can be related to the number of events in the sequence, which becomes smaller augmenting the value of the threshold magnitude, thus limiting the equiprobability of the possible transitions among the earthquake magnitudes.
42
M. Lovallo et al. / Chaos, Solitons and Fractals 24 (2005) 33–43
5.0
depleted full
Italy
4.8 4.6 4.4 4.2 4.0
Entropy
3.8 3.6 3.4 3.2 3.0 2.8 2.6 2.4 2.2 2.0 2.5
3.0
3.5
(a) 0.12
4.0
4.5
5.0
ML depleted full
Italy
0.11 0.10 0.09
Disequilibrium
0.08 0.07 0.06 0.05 0.04 0.03 0.02 0.01 2.5
3.0
3.5
(b)
4.0
4.5
5.0
ML full depleted
Italy 0.26 0.24 0.22
Complexity
0.20 0.18 0.16 0.14 0.12 0.10 0.08 0.06 2.5
(c)
3.0
3.5
4.0
4.5
5.0
ML
Fig. 7. Comparison between the complexity parameters shown in Figs. 6 and 7.
5. Conclusions The application of information theoretical measures to the transition proportion matrix of the magnitudes of seismic sequences has revealed the following issues: (i) the ‘‘complexity’’ of a seismic sequence can be dealt focusing on the properties of transition from one magnitude to another; (ii) the ‘‘complexity’’ of a seismic sequence can be measured
M. Lovallo et al. / Chaos, Solitons and Fractals 24 (2005) 33–43
43
by means of three quantities, the information entropy H, the disequilibrium D and the complexity C; (iii) these three parameters show a nontrivial relationship with the length of the seismic sequence, with the b-value of the Gutenberg– Richter law, which describes the frequency of earthquake magnitude in a seismic area, and with the threshold magnitude ML; (iv) aftershocks cause increase of entropy information and decrease of disequilibrium and complexity in a seismic sequence. In the future, the time variation of complexity parameters will be performed in order to reveal the particular seismic patterns.
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36]
Bialek W, Nemenman I, Tishby N. Complexity through Nonextensivity.