An effective chaos-geometric computational approach to analysis and prediction of evolutionary dynamics of the environmental systems: Atmospheric pollution dynamics

The present paper concerns the results of computational studying dynamics of the atmospheric pollutants (dioxide of nitrogen, sulphur etc) concentrations in an atmosphere of the industrial cities (Odessa) by using the dynamical systems and chaos theory methods. A chaotic behaviour in the nitrogen dioxide and sulphurous anhydride concentration time series at several sites of the Odessa city is numerically investigated. As usually, to reconstruct the corresponding attractor, the time delay and embedding dimension are needed. The former is determined by the methods of autocorrelation function and average mutual information, and the latter is calculated by means of a correlation dimension method and algorithm of false nearest neighbours. Further, the Lyapunov’s exponents spectrum, Kaplan-Yorke dimension and Kolmogorov entropy are computed. It has been found an existence of a low-D chaos in the time series of the atmospheric pollutants concentrations.


Introduction
The last decades have seen a great progress in the understanding, analysis, modelling and evet prediction of the evolutionary dynamics of nonlinear complex systems. Various methods and algorithms of the modern theory of dynamical systems and a chaos theory became a powerful tool in computational studying complex non-linear statistical systems [1][2][3][4][5][6][7][8][9][10][11][12][13][14]. Many studies in different fields of science and technique have appeared, where the chaos theory methods were applied to a great number of dynamical systems. The studies concerning non-linear behaviour in the time series of atmospheric constituent concentrations are sparse, and their outcomes are ambiguous. In ref. [5] there is an analysis of the NO 2 , CO, O 3 concentrations time series and is not received an evidence of chaos. Also, it was shown that O 3 concentrations in Cincinnati (Ohio) and Istanbul are evidently chaotic, and non-linear approach provides satisfactory results [6]. In Ref. [14] it has been fulfilled the detailed analysis of the NO 2 , CO, CO 2 concentration time series in the Gdansk region (Polland) and it has been definitely obtained the evidence of a chaos. Moreover it has been given a short-range forecast of atmospheric pollutants time evolution using non-linear prediction method. These studies show that chaos theory methodology can be applied and the short-range forecast by the non-linear prediction method can be satisfactory. Time series of concentrations are however not always chaotic, and chaotic behaviour must be examined for each time series.

The data for computational studying and method of testing chaos in time series
In our study, the nitrogen dioxide (NO 2 ) and sulphurous anhydride (SO 2 ) concentration data observed in the atmosphere of the Odessa city from 1976 till 2000 years [3]. In Refs. [2][3][4][11][12][13][14][15][16][17][18] it has been developed computational code for studying chaotic features of the complex non-linear systems and in details described a procedure of testing of the chaos elements in the corresponding time series. Here we are limited only by the key aspects. As usually, we consider scalar measurements s(n)=s(t 0 + nt) = s(n), where t 0 is a start time, t is time step, and n is number of the measurements. In a general case, s(n) is any time series, but here s(n) corresponds to an atmospheric pollutant concentration. The first fundamental step of modelling is in reconstruction of the corresponding phase space using as well as possible information contained in s(n). From the mathematical viewpoint, this procedure results in set of d-dimensional vectors y(n) replacing scalar measurements. One should further to operate with lagged variables s(n+), where  is some integer to be defined, results in a coordinate system where a structure of orbits in phase space can be captured. Using a set of the time lags to create a vector in d dimensions, y(n)=[s(n), s(n + ), s(n + 2),..,s(n +(d1))], the required coordinates are provided. The dimension d is defined as an embedding dimension, d E .
In Refs. [2][3][4] there are presented a few approaches to the choice of proper time lag. This point is important for the subsequent reconstruction of phase space. First approach [2] is to compute the linear autocorrelation function C L () and to look for that time lag where C L () first passes through 0. The alternative approach is based on using method of an average mutual information. Let us remind that the mutual information I of two measurements a i and b k is symmetric and non-negative, and equals to 0 if only the systems are independent. The average mutual information between any value a i from system A and b k from B is the average over all possible measurements of I AB (a i , b k ). In Ref. [4] it is suggested, as a prescription, that it is necessary to choose that  where the first minimum of I() occurs.
The fundamental goal of the d E calculation is in the further reconstruction of the Euclidean space R d large enough so that the set of points d A can be unfolded without ambiguity. The embedding dimension, d E , must be greater, or at least equal, than a dimension of the corresponding chaotic attractor, d A , i.e. d E > d A . The correlation integral analysis is one of the widely used techniques to investigate the signatures of chaos in a time series. This method is based on using the correlation integral, C(r). As usually, if the corresponding time series is characterized by an attractor, then the correlation integral C(r) is related to the radius r as where d is correlation exponent. In a case of the chaotic system the correlation exponent attains saturation with an increase in the embedding dimension. The saturation value of this exponent is defined as the correlation dimension (d 2 ) of the attractor (see details in Refs. [2][3][4][9][10][11][12][13][14][15][16][17][18][19][20][21][22]). The technique of application the correlation integral method (say, the Grassberger-Procaccia algorithm [9]) is presented in Refs. [2][3][4]8,9]. Another method for determining d E comes from asking the basic question addressed in the embedding theorem: when has one eliminated false crossing of the orbit with itself which arose by virtue of having projected the attractor into a too low dimensional space? This method is called as the method of false nearest neighbours [11]. As a rule, the simultaneous application of two methods provides more exact determination d E . It is noteworthy that the nearest integer above the saturation value provides the minimum or optimum embedding dimension for reconstructing the phase-space or the number of variables necessary to model the dynamics of the system. This concept can be applied, since the embedding dimension determined by both the correlation dimension method and the algorithm of false nearest neighbours are identical.
The further important step in studying the chaotic time series of the dynamical system is determination of predictability, which can be estimated by the Kolmogorov entropy. The Kolmogorov entropy is proportional to a sum of the positive Lyapunov's exponents. Let us remind that the Lyapunov's exponents spectrum is one of the fundamental dynamical invariants for non-linear system with chaotic behaviour.
According to definition, the Lyapunov's exponents are related to the eigenvalues of the linearized dynamics across the attractor. These parameters indicate the complexity of dynamics of the studied system. As usually, the positive values the Lyapunov's exponents show local unstable behaviour of the system, and respectively, their negative values show stable behaviour. The largest positive value of the Lyapunov's exponents determines some average prediction limit. Since the Lyapunov's exponents are defined as asymptotic average rates, they are independent of the initial conditions, and hence the choice of trajectory, and they do comprise an invariant measure of the attractor. An estimate of this measure is a sum of the positive Lyapunov's exponents. The estimate of the attractor dimension is provided by the conjecture d L and the Lyapunov's exponents are taken in descending order. The dimension d L gives values close to the dimension estimates discussed earlier and is preferable when estimating high dimensions. To compute Lyapunov's exponents, we use a method with linear fitted map [2,4,20], although the maps with higher order polynomials can be used too. time series remains positive. The values, where the autocorrelation function first crosses 0.1, can be chosen as , but in [1] it has been showed that an attractor cannot be adequately reconstructed for very large values of . So, before making up final decision we calculate the dimension of attractor for all values in Table 1.

Results
The outcome is explained not only inappropriate values of  but also shortcomings of correlation dimension method [2]. If algorithm [1] is used, then a percentages of false nearest neighbours are comparatively large in a case of large . If time lags determined by average mutual information are used, then algorithm of false nearest neighbours provides d E = 6 for all air pollutants.
10 20 Table 2 shows the calculated parameters: correlation dimension (d 2 ), embedding dimension (d E ), two Lyapunov exponents, E( 1 , 2 ), Kaplan-Yorke dimension (d L ), and average limit of predictability (Pr max , hours) for the NO 2 , SO 2 concentration time series in the Odessa region (for two measurement sites) during the period: Jan.-Dec., [1989][1990]. From the table 2 it can be noted that the Kaplan-Yorke dimensions, which are also the attractor dimensions, are smaller than the dimensions obtained by the algorithm of false nearest neighbours. It is very important to pay the attention on the presence of the two (from six) positive Lyapunov's exponents  i . This fact suggests that the system broadens in the line of two axes and converges along four axes that in the six-dimensional space. The time series of SO 2 at the site 2 have the highest predictability (more than 2 days), and other time series ( in particular, the NO 2 concentration) have the predictabilities slightly less than 2 days. by using the dynamical systems and chaos theory methods. A chaotic behaviour in the nitrogen dioxide and sulphurous anhydride concentration time series at several sites of the Odessa city is numerically investigated for the first time. As usually, to reconstruct the corresponding attractor, the time delay and embedding dimension were determined. The time delay was calculated on the basis of methods of autocorrelation function and average mutual information, and the embedding dimension was calculated by means of the correlation dimension method and algorithm of false nearest neighbours. Further, the Lyapunov's exponents spectrum, Kaplan-Yorke dimension and Kolmogorov entropy are calculated. Our computational study has shown an existence of a low-and high-D chaos in the atmospheric pollutants fluctuations dynamics in the Odessa. This conclusion allows further to develop the corresponding prediction models [13][14][15][16][17][18] for description of the temporal evolutionary dynamics of the air pollutants concentration in atmosphere of the industrial city.