Chaos, Solitons & Fractals 45 (2012) 1086–1091
Contents lists available at SciVerse ScienceDirect
Chaos, Solitons & Fractals Nonlinear Science, and Nonequilibrium and Complex Phenomena journal homepage: www.elsevier.com/locate/chaos
Visibility graph approach to the analysis of ocean tidal records Luciano Telesca a,⇑, Michele Lovallo b, Jorge O. Pierini c a
National Research Council, Institute of Methodologies for Environmental Analysis – C.da S.Loja, 85050 Tito, Italy Agenzia Regionale per la Protezione dell’Ambiente di Basilicata, Via della Fisica 18 C/D – 85100, Potenza, Italy Universidad Nacional del Sur, Departamento de Física, CCT BB, Comisión de Investigaciones Científicas de la Provincia de Buenos Aires, C.C. 804, Florida 8000 Complejo CCT, B8000FWB Bahía Blanca, Argentina
b c
a r t i c l e
i n f o
Article history: Received 1 March 2012 Accepted 15 June 2012 Available online 20 July 2012
a b s t r a c t By using the recent method of the visibility graph, three time series of oceanic tide level in central Argentina were investigated. The degree distributions show a rich structure; in particular the maximum is due to the main periodic oscillations at 24 hours and 12 hours and higher harmonics. The degree distributions of the residuals (obtained removing from the original signals the cyclic components) suggest that the local effects, linked with the particular coastal conditions of the sites, are discernible for the degree k < 20, while the global effects, linked with linked with the more general and common atmospheric forcing and ocean current conditions, are visible for k > 100. Although a relationship between the spectral exponent a and the exponent of the degree distribution c of tidal signals can be recognized, this cannot be simply stated due to the very rich and complex structure of time dynamics of tides. The present study, even if still preliminary, show the importance of the visibility graph method in investigating the complex time dynamics of observational and experimental signals. Ó 2012 Elsevier Ltd. All rights reserved.
1. Introduction The seminal paper by Lacasa et al. [1] has stimulated an increasing interest in the application of the method of visibility graph (VG) to time series. Such method maps time series into networks, converting the properties of time series in properties of networks, and, viceversa, using features of networks to get information about the time series. In the VG approach a segment links any two values of the series that can be seen by each other, meaning that such segment is not broken by any other intermediate value of the series. In terms of graph theory, each value of the time series represents a node, and two nodes are connected if there exists visibility between them. The visibility criterion [1] can be mathematically defined as follows: two arbitrary data values (ta, ya) and (tb, yb) are visible to each
⇑ Corresponding author. E-mail address:
[email protected] (L. Telesca). 0960-0779/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.chaos.2012.06.003
other if any other data (tc, yc) placed between them fulfils the following constrain:
yc < yb þ ðya yb Þ
tb tc : tb ta
ð1Þ
The following properties always hold [1]: 1. Connection: each node is visible at least by its nearest neighbors (left and right); 2. Absence of directionality: no direction is defined in the links; 3. Invariance under affine transformations (rescaling of both axes and horizontal and vertical translations) of the time series. It was shown that the graph developed using the VG method transform periodic, random and fractal time series into regular, random and scale-free networks respectively [1–3]. Another property concerns more specifically the fractional Brownian motions (fBm), whose degree k distribution is a power-law, P(k) kc, where c is the powerlaw exponent that depends on the Hurst exponent H of the fBm through a linear relationship, c = 3 2H, leading to the relationship c = 4 b in terms of the spectral
L. Telesca et al. / Chaos, Solitons & Fractals 45 (2012) 1086–1091
exponent b of 1/fb noise [4]. The VG method has been applied in several scientific fields and recently to the analysis of time dynamics of earthquakes [5]. One of the most interesting advances in network theory is the mapping of a time series into a network. In the case examined in this paper such mapping is the VG that is based on the convexity of successive observations [2]. We can indicate as MVG the map from the space of time series T to the space of networks N, MVG: T N, where MVG associates to a time series x(t) a set of nodes and links, in which each node is an individual observation of the time series and the link is established according to a local convexity constraint between successive observations. Such proposed map is surjective, because at each time series x(t) always a network corresponds, but it is not injective, because two different time series, x(t) and ax(t) (the second one is just the scaled version of the first one), can be mapped into the same network. 2. Data analysis In this study we apply the VG method to the tidal time series recorded in three sites at Bahia Blanca (Argentina) estuary (Ingeniero White, Puerto Belgrano and Torre Mareografica) from October 8, 2002 to April 10, 2003 (Fig. 1). Fig. 2 shows the time series and Fig. 2 shows their power
1087
spectral densities. The power spectra reveal a rather rich structure of the tidal time variability. Main oscillations at 24 hours and 12 hours along with their harmonics (being the 6-h cycle is the higher one) are superimposed to approximately three frequency regions discernible at low, intermediate and high frequency bands (separated by the vertical arrowed bars in Fig. 3); the intermediate frequency region between about 24 hours and 1 hour can be approximated by a power law with scaling exponent a 5 (Ingeniero White), 4 (Puerto Belgrano) and 3.6 (Torre Mareografica). We applied the VG method to all the three tidal time series and Fig. 4 shows the distribution of the degree k. First of all, the complex and rich structure of the investigated time series is reflected in the non trivial degree distributions, which show several features: (i) for k < 20 the three distributions are vertically shifted from each other; in particular, Ingeniero White and Torre Mareografica are respectively the lower and the upper curve; this indicates that small k is less frequent at Ingeniero White than Torre Mareografica; (ii) for 20 < k < 60, the three curves appear horizontally shifted, with the maximum of the degree distribution at k = 36 for Ingeniero White, k = 41 for Puerto Belgrano and k = 46 for Torre Mareografica;
Fig. 1. Map of the tidal stations in central Argentina: Ingeniero White, Puerto Belgrano and Torre Mareografica in Bahia Blanca Estuary.
1088
L. Telesca et al. / Chaos, Solitons & Fractals 45 (2012) 1086–1091
Ingeniero White
Ingeniero White
12 h
6
10
6
5
10
4
10
5
6h
8h
24 h
3
10
2
4h
10
4
1
10
0
10
3
-1
S(f)
m
10
2
-2
10
-3
10
-4
10
1
α=5
-5
10
-6
10
0
-7
10
-8
10
(a)
-1 0.0
3
5.0x10
4
1.0x10
4
1.5x10
4
2.0x10
4
-10
10
4
2.5x10
(a)
-9
10
3.0x10
-6
10
t (10 min)
-5
10
-4
-3
10
-2
10
-1
10
10
-1
f(min ) Puerto Belgrano
5 4
S(f)
m
3 2 1 0
(b) -1 0.0
3
5.0x10
4
1.0x10
4
1.5x10
4
2.0x10
4
4
2.5x10
3.0x10
10 5 10 4 10 3 10 2 10 1 10 0 10 -1 10 -2 10 -3 10 -4 10 -5 10 -6 10 -7 10 -8 10 -9 10 -10 10 -11 10
6h
8h
24 h
4h
α=4.03
(b) -6
10
t (10 min)
Puerto Belgrano
12 h
6
6
-5
10
-4
-3
10
-2
10
-1
10
10
-1
f(min ) Torre Mareografica
12 h
5
Torre Mareografica
5
10
4
10
8h
3
4
10
6h
24 h
2
10
4h
1
10
3
0
m
10
-1
10
2
S(f)
-2
1
10
-3
10
-4
10
-5
10 0
α=3.6
-6
10
-7
(c)
10
-8
-1
10 0.0
3
5.0x10
4
1.0x10
4
1.5x10
4
2.0x10
4
2.5x10
4
3.0x10
t (10 min)
-10
10
0.000001
Fig. 2. Time series of oceanic tide in (a) Ingeniero White, (b) Puerto Belgrano and (c) Torre Mareografica.
(iii) for k > 60 the three distributions seem to collapse; (iv) a further relative maximum can be identified at k = 2, and is identical for all the sites. It seems reasonable that the most powerful periodicities, as identified in the power spectra (Fig. 3), should have an important role in depicting the particular form of the degree distributions. In order to better understand the particular shapes of the degree distribution of the three
(c)
-9
10
0.00001
0.0001
0.001
0.01
0.1
-1
f(min ) Fig. 3. Power spectral density of the data shown in Fig. 1: (a) Ingeniero White, (b) Puerto Belgrano and (c) Torre Mareografica.
tidal signals, we decomposed each tidal signal into two components: one periodical component (indicated as C hereafter), given by the sum of the main oscillations at 24 hours, 12 hours and their higher harmonics (8, 6 and 4 hours) (Fig. 5); and one residual component (indicated as R hereafter), obtained from the original tides after removing the cyclic component C (Fig. 6).
1089
L. Telesca et al. / Chaos, Solitons & Fractals 45 (2012) 1086–1091 3
Ingeniero White-original Puerto Belgrano-original Torre Mareografica-original
Ingeniero White
2
-1
10
1
0 -2
P(k)
C
10
-1 -3
10
-2
-3
-4
(a)
10
0.0 0
10
1
3
5.0x10
1.0x10
4
1.5x10
4
2.0x10
10
k
The degree distributions are not shaped as simple power law or simple decreasing exponential functions,
4
3.0x10
Puerto Belgrano
1
C
0
-1
-2
(b) -3 0.0
3
5.0x10
4
1.0x10
4
1.5x10
4
2.0x10
4
2.5x10
4
3.0x10
t (10 min) 2
Torre Mareografica
1
0
C
(i) the maximum of the degree distribution of the original signal is an effect of the periodical component C; in fact, the maximal value of the degree distribution of C is overlapped on that of the degree distribution of the original series; (ii) the further relative maximum, which in the distribution of the original series is located at k = 2, is an effect of the residual component R, whose degree distribution has the first maximum located at the same value k = 2; in particular the relative maximum at k = 2 for the R signals has the lower value for Ingeniero White and the higher one for Torre Mareografica; (iii) the periodical component C is responsible of the horizontal shift observed in the degree distributions of the original signals for 20 < k < 60; in fact, the comparison among the three degree distributions of C (Fig. 9) shows the same effect in the same k range; (iv) the residual component R is responsible of the vertical shift observed in the degree distributions of the original signals for k < 20; in fact, the comparison among the three degree distributions of R (Fig. 8) shows the same phenomenon in the same k range; (v) the range of k values is significantly larger for R, while is slightly smaller for C respect to the range of k values of the original series.
2.5x10
2
Fig. 4. Degree distribution of the tidal time series recorded at Ingeniero White (blue), Puerto Belgrano (red) and Torre Mareografica (green). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
We applied, then, the VG method to the components C and R of each tidal series. The analysis of periodic series has been performed by Luque et al. [11] and Nuñez et al. [12], who applied a modified VG method, called the horizontal VG. Fig. 7 shows for each station the comparison among the degree distributions obtained for the original signal, the periodical component C and the residual R. The following information can be deduced:
4
t (10 min)
2
10
4
-1
-2
(c) 0.0
3
5.0x10
4
1.0x10
4
1.5x10
4
2.0x10
4
2.5x10
4
3.0x10
t (10 min) Fig. 5. Periodical signals extracted from the original tidal signals shown in Fig. 2.
therefore their interpretation is not an easy task; however, regarding the degree distributions of the R signals, from a visual inspection it seem possible to recognize a k-range in which the power-law behavior holds, i.e. for k > 100. We calculated, then, for k > 100 the exponent c by using the maximum likelihood estimation (MLE). Recently it has been shown that the MLE is a very good estimator of the scaling parameters of a power-law distribution, and
1090
L. Telesca et al. / Chaos, Solitons & Fractals 45 (2012) 1086–1091 6
Ingeniero White
Ingeniero White original R C
-1
10
5
4
-2
P(k)
R
10 3
-3
10 2
-4
1
10
(a)
(a)
0 0.0
3
5.0x10
4
1.0x10
4
1.5x10
4
2.0x10
4
2.5x10
4
0
3.0x10
1
10
2
10
t (10 min)
10
3
10
k Puerto Belgrano
Puerto Belgrano original R C
5 -1
10 4
-2
P(k)
R
10 3
-3
10 2
-4
10
(b)
1
0.0
3
5.0x10
4
1.0x10
4
1.5x10
4
2.0x10
4
2.5x10
(b)
4
3.0x10
0
1
10
t (10 min)
2
10
10
3
10
k Torre Mareografica 4
Torre Mareografica original R C
-1
10 3
-2
2
P(k)
R
10
-3
10 1
(c)
0 0.0
3
5.0x10
4
1.0x10
4
1.5x10
4
2.0x10
4
2.5x10
-4
10
(c)
4
3.0x10
t (10 min)
0
10
Fig. 6. Residual signals obtained from the original tides after removing the cyclic shown in Fig. 5.
performs better than the least-square method (LSM) applied to the probability density function with constant or logarithmic bin width [10]. Therefore, using the MLE, the exponent c can be calculated by means of the following equation [6]
cMLE ¼ 1 þ n
" n X i
ki ln kmin 0:5
#1 ;
1
ð2Þ
where kmin is the smallest value of the node k for which the power law holds and n is the number of values k P kmin. The estimation error is given by [6]
2
10
10
3
10
k Fig. 7. Degree distributions for the three tidal signals along with those of the periodical C and residual R parts.
rMLE ¼
cMLE 1 pffiffiffi n
:
ð3Þ
In our case kmin = 100 and n = 1417, 1448 and 1453 for the residuals of Ingeniero White, Puerto Belgrano and Torre Mareografica respectively. The values for the exponent c is 3.48 ± 0.06 (Ingeniero White), 3.49 0.06 (Puerto Belgrano) and 3.48 0.06 (Torre Mareografica). The three values are almost identical due to the collapse of the degree distributions for k > kmin.
L. Telesca et al. / Chaos, Solitons & Fractals 45 (2012) 1086–1091
Residual Ingeniero White Puerto Belgrano Torre Mareografica
-1
10
-2
P(k)
10
-3
10
-4
10
0
10
1
2
10
3
10
10
k Fig. 8. Comparison among the three degree distributions of the residual parts.
Periodical Ingeniero White Puerto Belgrano Torre Mareografic -1
10
-2
1091
range temporal correlations between earthquake magnitudes [8,9]). The difference among the distributions could suggest that the range k < 20 may be more informative of the local tidal effects linked with the coastal characteristics of the site, while the range k > 100 may be more informative of the global tidal effects, linked with the more general and common atmospheric forcing and ocean current conditions. Furthermore, it should be recognized that in this study a quantitative relationship between the spectral exponent a and the degree distribution exponent c cannot be simply stated, neither empirically because only three observational signals were analysed, nor theoretically, because the signals present a very rich and complex structure, not simply classifiable as a fractional Brownian motion, for which such theoretical relationship was derived [4]. However, we observe that a increases with the decrease of c. 3. Conclusions A new approach based on the method of the VG was proposed for the investigation of time dynamics of tidal time series. In the future more investigations should be performed to link the several features in the VG degree distribution of tidal signals and to understand their links with patterns of their power spectra. Acknowledgements
P(k)
10
This study was supported by the Bilateral Project 20112012 CNR/CONICET ‘‘Development of models and time series analysis tools for oceanographic parameters analysis and forecasting’’.
-3
10
-4
10
References 0
10
1
2
10
10
k Fig. 9. Comparison among the three degree distributions of the periodical parts.
Note that the values of the scaling exponent deduced here from the three measuring sites are comparable to the value of the exponent obtained from the earthquake data [5]. In particular, the average exponent for the shuffled earthquake magnitude sequences is 3.18 which is almost identical to that of the original earthquake series (while in natural time analysis [7] the difference of the results obtained between the shuffled earthquake data and the original data clearly reveals the importance of long
[1] Lacasa L, Luque B, Ballesteros F, Luque J, Nuño JC. Proc. Natl. Acad. Sci. U.S.A. 2008;105:4972. [2] Donner RV, Small M, Donges JF, Marwan N, Zou Y, Xiang R, et al. Int. J. Bifurcation Chaos 2011;21:1019–48. [3] A.S.L.O. Campanharo, M.I. Sirer, R.D. Malgren, F.M. Ramos, L.A.N. Amaral, Plos One 6 e23378, 2011. [4] Lacasa L, Luque B, Luque J, Nuño JC. Europhys. Lett. 2009;86:30001. [5] Telesca L, Lovallo M. Europhys. Lett. 2012;97:50002. [6] Newmann MEJ. Contemp. Phys. 2005;46:323. [7] Varotsos PA, Sarlis NV, Skordas ES. Phys. Rev. E 2002;66:011902. [8] Varotsos PA, Sarlis NV, Skordas ES, Tanaka HK, Lazaridou MS. Phys. Rev. E 2006;74:021123. [9] Sarlis NV, Skordas ES, Varotsos PA. Phys. Rev. E 2009;80:022102. [10] Clauset A, Shalizi CR, Newmann NEJ. SIAM Rev. 2009;51:661–703. [11] Luque B, Lacasa L, Ballesteros FJ, Robledo A. PlosOne 2011;6:e22411. [12] A. Nuñez, L. Lacasa, E. Valero, J.P. Gómez, B. Luque, 2011, http:// arxiv.org/abs/1108.1693.