The correlation dimension D is commonly used to quantify the chaotic structure of atmospheric time series. The standard algorithm for estimating the value of D is based on finding the slope of the curve obtained by plotting ln C(r) versus ln r, where C(r) is the correlation integral and r is the distance between points on the attractor. An alternative, probabilistic method proposed by Takens is extended and tested here. This method is based on finding the sample means of the random variable (r/ρ)p[ln(r/ρ)]k, expressed as the conditional expected value E((r/ρ)p[ln(r/ρ)]k : r < ρ), for p and k nonnegative numbers.
The sensitivity of the slope method and of the extended estimators Dpk(ρ) for approximating D is studied in detail for three ad hoc correlation integrals and for integer values of p and k. The first two integrals represent the effects of noise or undersampling at small distances and the third captures periodic lacunarity, which occurs by definition when the ratio C(xρ)/C(ρ) fails to converge as ρ approaches zero. All the extended estimators give results that are superior to that produced by the most commonly used slope method. Moreover, the various estimators exhibit much different behavior in the two ad hoc cases: noise-contaminated signals are best diagnosed using D11(ρ), and lacunar signals are best studied using D0k(ρ), with k as large as possible in magnitude. Therefore, by using a wide range of values of p and k, one can infer whether degradation arising from noise or arising from lacunarity is more pronounced in the time series being studied, and hence, one can decide which of the estimates most efficiently approximates the correlation dimension for the series.
These ideas are applied to relatively coarse samplings of the Hénon, Lorenz convection, and Lorenz climate attractors that in each case are obtained by calculating the distances between pairs of points on two trajectories. As expected from previous studies, lacunarity apparently dominates the Hénon results, with the best estimate of D, D = 1.20 ± 0.01, given by the case D03(ρ). In contrast, undersampling or noise apparently affects the Lorenz convection and climate attractor results. The best estimates of D are given by the estimator D11(ρ) in both cases. The dimension of the convection attractor is D = 2.06 ± 0.005, and that of the climate attractor is D = 14.9 ± 0.1. Finally, lagged and embedded time series for the Lorenz convection attractor are studied to identify embedding dimension signatures when model reconstruction is employed.
In the last part of this study, the above results are used to help identify the best possible estimate of the correlation dimension for a low-frequency boundary layer time series of low-level horizontal winds. To obtain such an estimate, Lorenz notes that an optimally coupled time series must be extracted from the data and then lagged and embedded appropriately. The specific kinetic energy appears to be more closely coupled to the underlying low-frequency attractor, and so more nearly optimal, than is either individual wind component. When several estimates are considered, this kinetic energy series exhibits the same qualitative behavior as does the lagged and embedded Lorenz convective system time series. The series is either noise contaminated or undersampled, a result that is not surprising given the small number of time series points used, for which the best estimate is given by D11(ρ). The obtained boundary layer time series dimension estimate, 3.9 ± 0.1, is similar to the values obtained by some other investigators who have analyzed higher-frequency boundary layer time series. Although this time series does not contain as many points as might be required to accurately estimate the dimension of the underlying attractor, it does illustrate the requirement that in any estimate of the correlation dimension, a function of the measured variables must be chosen that is strongly coupled to the attractor.
Estimates of the correlation dimension D of a time series are commonly used for quantifying the chaotic complexity of a physical system such as the atmosphere (e.g., Grassberger and Procaccia 1983a,b; Eckmann and Ruelle 1985; Fraedrich 1986; Hense 1987; Henderson and Wells 1988; Tsonis and Elsner 1988; Fraedrich and Leslie 1989; Fraedrich et al. 1990; Baker and Gollub 1990; Sharifi et al. 1990; Yang 1991; Tsonis 1992; Zeng et al. 1993; Poveda-Jaramillo and Puente 1993; Orcutt and Arritt 1995; Weber et al. 1995). However, owing to the presence of both intrinsic and extrinsic sources of error, the merits of this procedure have been debated considerably in the literature (e.g., Nicolis and Nicolis 1984, 1987; Grassberger 1986, 1987; Smith et al. 1986; Theiler 1986, 1988; Procaccia 1988; Smith 1988; Tsonis and Elsner 1988, 1989, 1990; Pool 1989; Nerenberg and Essex 1990; Ruelle 1990; Essex and Nerenberg 1991; Zeng et al. 1992; Islam et al. 1993; Tsonis et al. 1993; Weber et al. 1995).
The method most commonly used to estimate the correlation dimension for atmospheric and other geophysical time series has been the mean slope method, which consists of approximating the dimension using mean slopes of the graph of ln C(r) versus ln r for small ranges of r, where C(r) is the correlation integral and r is the distance between time series points (e.g., Grassberger and Procaccia 1983a,b; Nicolis and Nicolis 1984; Eckmann and Ruelle 1985; Fraedrich 1986; Hense 1987; Grassberger 1986; Nese 1987; Simm et al. 1987; Henderson and Wells 1988; Fraedrich and Leslie 1989; Keppene and Nicolis 1989; Tsonis and Elsner 1989; Nerenberg and Essex 1990; Sharifi et al. 1990; Fraedrich et al. 1990; Albano et al. 1991; Cutler 1991; Lorenz 1991; Nerenberg et al. 1991; Zeng et al. 1992; Sundermeyer and Vallis 1993; Poveda-Jaramillo and Puente 1993; Weber et al. 1995). However, in actual applications, this method may be quite sensitive to undersampling or superimposed noise.
In this article we extend the novel alternative method for estimating the value of D that was introduced by Takens (1985) and that appears to be well suited for use with atmospheric and other natural time series. The extended method provides an infinite number of functionally independent estimates of the value of D, any one of which is superior to the ones given by the most commonly used mean slope method. Moreover, when the time series is contaminated by noise, these estimates display behavior that is much different from the behavior when the time series is characterized by lacunarity; lacunarity may be viewed as insufficiently damped, intrinsic irregularity in the chaotic structure as the scale decreases (Theiler 1988). This estimator behavior allows one to discern the best approximate value for the correlation dimension from a set of estimates provided by the extended method.
a. Sources of errors
There are numerous sources of error inherent in an approximation of the attractor dimension. These errors arise from such things as time series collection and preprocessing, embedding and model reconstruction, as well as the specific dimension algorithm used. Under time series collection, we include the frequency, number, continuity, and duration of the measurements (Nicolis and Nicolis 1984, 1987; Grassberger 1986, 1987; Fraedrich 1986; Smith et al. 1986; Essex et al. 1987; Smith 1988; Keppene and Nicolis 1989; Ramsey and Yuan 1989; Fraedrich et al. 1990; Nerenberg and Essex 1990; Ruelle 1990; Essex 1991; Zeng et al. 1992; Poveda-Jaramillo and Puente 1993; Tsonis et al. 1993; Orcutt and Arritt 1995; Weber et al. 1995) and contamination by noise (Franaszek 1984; Ben-Mizrachi et al. 1984; Theiler 1986; Simm et al. 1987; Osborne and Provenzale 1989; Theiler 1991); under preprocessing, we include low- or high-pass filtering (Hense 1987; Theiler 1991), interpolation and trend removal, and principal component analysis (Albano et al. 1988); under model reconstruction, we include the choice of variable or variables (Lorenz 1991; Islam et al. 1993; Poveda-Jaramillo and Puente 1993), embedding dimension (Mañé 1981; Takens 1981; Grassberger and Procaccia 1983a,b; Keppene and Nicolis 1989; Nerenberg and Essex 1990; Yang 1991; Islam et al. 1993; Poveda-Jaramillo and Puente 1993; Tsonis et al. 1993), and phase lag (Broomhead and King 1986; Fraser and Swinney 1986; Albano et al. 1991; Thomson and Henderson 1992; Tziperman et al. 1995); and under the dimension algorithm used, we include those proposed or used by Kaplan and Yorke (1979), Russell et al. (1980), Takens (1981, 1985), Greenside et al. (1982), Grassberger and Procaccia (1983a,b), Termonia and Alexandrowicz (1983), Holzfuss and Mayer-Kress (1986), Theiler (1987), Henderson and Wells (1988), Kember and Fowler (1992), Orcutt and Arritt (1995), Weber et al. (1995), and us in this paper. Yet another source of error is lacunarity of the correlation integral (Theiler 1988), which causes the Takens (1985) estimator to fail, though possibly only to a small degree. In this article, we consider noise, lacunarity, and embedding as sources of error in a time series.
b. Objectives of this study
To compare the sensitivities to superimposed noise of the mean slope method and of the Takens estimators and their extensions introduced here, we consider the case of a point moving randomly on the unit sphere in Euclidean three-space, perturbed by three-dimensional Gaussian noise. We find that, while the mean slope method is almost lethally sensitive to noise or undersampling, the other methods are much less sensitive. The least sensitive method is given by a certain extension of the Takens estimators; although the behavior of this extension is clearly indicated, we have no theoretical reason to explain this result. In all cases, our results illustrate the arguments of Ben-Mizrachi et al. (1984) and Simm et al. (1987) that the dimension to be expected from the noise-contaminated range is the embedding dimension M, which is greater than the actual correlation dimension; Essex (1991), however, points out that this contaminated estimate of the dimension could be smaller than the actual correlation dimension, as is the case in the presence of severe undersampling.
We also investigate the success of the extensions of the Takens estimators when the correlation integral is periodically lacunar (Theiler 1988). In this case a different extension of the Takens estimators seems to damp markedly the oscillations in the dimension estimate. As before, this is an empirical result for which we have no a priori explanation. We conclude from our analysis that, when both noise and periodic lacunarity contribute to a signal, it is wise to calculate many functionally independent estimates of the correlation dimension using several of the extensions of the Takens estimators. We demonstrate the efficacy of this procedure by finding several estimates of D for the Hénon (1976), Lorenz (1963) convection, and Lorenz (1991) climate attractors. Finally, we consider the issues of variable choice and embedding dimension when we apply the extended Takens estimators to the model-reconstructed Lorenz convection attractor and to observed low-frequency atmospheric boundary layer flow. These latter results were reported in a preliminary fashion in Fosmire et al. (1995).
2. Summary of the new algorithm
Any estimate of the correlation dimension D is based on the proportion C(r) of distances, less than a fixed distance r, between two points on an attractor in phase space. The value of C(r) is hypothesized to obey a power law of the form C(r) ∼ rD, and approximate slope methods for finding D may be regarded as implementations of the (exact) slope method that finds D as a limit given by
Of course, if C(r) is not differentiable, then one must seek an alternative procedure.
More general and robust formulas than (2.1) for the correlation dimension may be derived by applying that formula to a different integral, A(r) = ∫r0 C(s) ds. Because C(r) ∼ rD, it is easy to see that A(r) ∼ rD+1. Application of (2.1) yields
Using integration by parts, we obtain
and consequently, we have
At this point, we may invoke the very deep Karamata theory (Bingham et al. 1987) to conclude, for reasonable f (x), that
We use a probabilistic interpretation of (2.5) by recalling that C(r) is the probability (with respect to the ergodic measure) that a random pair of points on the attractor has distance less than or equal to r. Then the conditional expectation E(f(r/ρ) : 0 < r < ρ) is precisely the integral appearing on the left side of (2.5); by the Ergodic Theorem (Guckenheimer and Holmes 1983), this expectation may be found very rapidly by taking suitable (time) averages along a trajectory. As for the right side of (2.5), for many choices of f(x), the integral is easily found.
We may now make explicit our integral algorithm for estimating the correlation dimension D. Although we have used differentiation in deriving this algorithm, we do not use differentiation, either numerical or exact, when applying it; instead, we use the numerically more robust processes of integration or evaluation of averages (section 2b). Defining the particular conditional expected values
for nonnegative p and k, from (2.5) we obtain
in which Γ(k + 1) is the well-known gamma function (e.g., Abramowitz and Stegun 1964). These expected or average values are approximated by simply summing the quantities (r/ρ)p|ln(r/ρ)|k for distances less than ρ along the trajectory and dividing by the number N(ρ) of such distances. Then, for small values of ρ, the value of Epk(ρ) approximates the limit in (2.7). On setting this value equal to the right side of (2.7) and solving for D, we obtain the corresponding approximation Dpk(ρ) of D. When k = 0 or p = 0, we have the simpler formulas for the dimension estimates Dp0(ρ) and D0k(ρ) given by
Here, the two Takens (1985) estimators are D01 and D02. Finally, for positive, integer values of k, we may use the formula
As we demonstrate in section 3, the estimates Dp0(ρ), D0k(ρ), and Dpk(ρ) vary differently depending on whether noise, undersampling, or lacunarity plague a time series signal. The best estimates for noise-contaminated or undersampled signals are given by the case D11(ρ) for which p = k = l, while the best estimates for lacunarity-dominated signals are given by the case D0k(ρ) for which p = 0 and k is as large as possible; typically k = 3 is found to be adequate.
Given a finite sequence of n positive distances r1, r2, . . . , rn, the correlation integral C(r) in practice is estimated from the cumulative probability distribution function of these distances. The conditional expected values Ep0(ρ), E0k(ρ), and Epk(ρ) are approximated, respectively, by the corresponding sample means [Σri<ρ (ri/ρ)p]/N(ρ), [Σri<ρ |ln(ri/ρ)|k]/N(ρ), and [Σri<ρ (ri/ρ)p|ln(ri/ρ)|k]/N(ρ), where for a given ρ only the N(ρ) distances satisfying ri < ρ are included in the sums. Normally, these sums are approximated using 100 or 200 values of ρ that evenly partition the n original distances between the minimum rmin and maximum rmax values of interest. We use these means in the formulas (2.7)–(2.10) to obtain the dimension approximations Dp0(ρ), D0k(ρ), and Dpk(ρ) for rmin < ρ ≤ rmax. That is, we solve
for each value of ρ in the above range.
We may refine these formulas slightly to account for the fact that we are unable to sample small distances between 0 and the minimum distance rmin retained in the sample; this distance is not necessarily the smallest distance in the sample. We let a be a number such that aρ ≥ rmin for each ρ used in the above sums. Then we have, for example, that the approximation for Dp0(ρ) becomes
We use Newton’s method to find Dp0(ρ) with a first guess given by (2.11). We note that applying Newton’s method is not as laborious as it might appear because we may construct the recursion formula using a symbolic manipulator and then ask the manipulator to write in a BASIC or FORTRAN format the corresponding version of that (rather complicated) formula. Other formulas are given in the appendix.
In the following section, we use the above approach to investigate the dependence of the various estimates for D when the distance data are contaminated by errors.
3. Sensitivity of the estimators to noise and lacunarity
In this section we investigate the behavior of the estimators Dp0(ρ), D0k(ρ), and Dpk(ρ) as the values of p and k in (A.1)–(A.5) are varied for three idealized cases of C(r). In the first case, C(r) is the probability, with respect to a certain probability measure, that a point is within r units of the unit vector (1, 0, 0) in Euclidean space; in this case, strictly speaking, we investigate the approximations to the pointwise correlation dimension at the point (1, 0, 0). This probability measure is that given by the perturbation of the unit probability (surface) measure on the unit sphere by isotropic, three-dimensional Gaussian noise with standard deviation σ. As a second model of noise contamination at small r, we use a continuous ad hoc correlation integral C(r, δ) obeying power laws rc for r ≤ δ and rD for r ≥ δ, where c represents the embedding dimension with c > D (Ben-Mizrachi et al. 1984; Simm et al. 1987; and Theiler 1986, 1991). We use these two cases of C(r) to compare the sensitivities of the various estimators to noise.
In the third case, C(r) is an ad hoc, periodically lacunar correlation integral; a lacunar correlation integral is one for which the following limit does not exist for some 0 ≤ x ≤ 1: limρ→0 C(xρ)/C(ρ). Graphs of lacunar correlation integrals are consistently wavy, rather than linear, at any scale of r when ln C(r) is plotted as a function of ln r; such plots must become linear, rather than wavy, in the limit as r → 0 for the exact slope method (2.1) to work. One example of lacunar behavior can be found in Fig. 3a of Lorenz (1991), in which the correlation integrals based on a weakly coupled variable of the reconstructed climate attractor become increasingly more wavy as the embedding dimension increases; in contrast, the correlation integrals based on a strongly coupled variable do not exhibit this wavy behavior. For the lacunar case, it should not be surprising that the estimators introduced in Takens (1985), as well as those introduced here, do not apply (Theiler 1988). However, it is not hard to show that we may still use the extended method formulas to provide bounds on the dimension estimate. Thus, by applying various estimators to a lacunar case of C(r), we may form some idea of their sensitivities to lacunarity. Perhaps surprisingly, the sensitivities to lacunarity appear to be much the reverse of the sensitivities to noise. We have no analytic explanation for this empirical phenomenon, but it does suggest that these relative sensitivities may serve to indicate whether nonconvergence of a particular dimension calculation occurs because of noise contamination or because of lacunarity in the correlation integral. Once identified, the estimators that are least sensitive to the apparent obstacle, noise or lacunarity, should be used in a particular application.
a. Sensitivity to noise
We construct our correlation integral for the first case, the investigation of the sensitivity to noise, by letting p(x, y, z; x′, y′, z′; σ) be the isotropic Gaussian probability density centered at the point (x′, y′, z′) with standard deviation σ,
Then we average this density over all points (x′, y′, z′) on the unit sphere (i.e., the boundary of the unit ball). That is, we construct a new probability density P(x, y, z ; σ) by writing
where the integral is taken with respect to normalized surface area S on the unit sphere. Then we define our (pointwise) correlation integral as the probability that a given point lies within the distance r of the point (1, 0, 0),
The actual integrals are very laborious, and we use the symbolic manipulator Derive® to construct the function C(r, σ).
We take the scale of noise to be indicated by the value of σ. Then at large scales—large ρ in comparison with σ—the effects of noise on the dimension estimates are negligible and we expect the value of any dimension estimate to approach 2, which is the dimension of the surface of the ball, as the value of σ approaches zero; equivalently, holding the value of σ fixed, we expect the dimension estimate to approach 2 as the value of ρ becomes large. In contrast, as the value of ρ becomes small, the scale of the noise becomes large in comparison with ρ, the probability density approaches that of a Gaussian distribution in Euclidean three-space, and we expect the dimension estimate to approach the value 3. This property is consistent with the observations of Ben-Mizrachi et al. (1984), Simm et al. (1987), and Theiler (1986, 1991), who note that the effects of noise lead to larger-than-correct estimates for the correlation dimension. Because the actual dimension of our probability distribution is 2, cloaked by noise, we will regard as less sensitive to noise those estimates yielding values that lie closer to the value 2.
The correlation integral (3.3) is illustrated by the curve labeled σ = 0.05 in Fig. 1 for the case a = 0.01 and σ = 0.05. It can be seen that the curve has a single knee or point of maximum curvature at a value near r = 0.12. As the graph is plotted on a larger and larger scale, it approaches the graph of two straight lines, of slopes 3 and 2, that meet at the point of maximum curvature. This result can be seen by comparing the σ = 0.05 curve with the one labeled δ = 0.1, which is a curve having a slope of 3 for r < 0.1 and a slope of 2 for r > 0.1.
To compare the sensitivities of the various estimators to noise, we insert the function given by (3.3) into (2.6) with the lower limits replaced by aρ to determine analytically the expected values Ep0(ρ), E0k(ρ), and Epk(ρ). Then we use these formulas on the left sides of (A.1)–(A.5) and solve these equations using Newton’s method to find the approximate dimensions Dp0(ρ, σ), D0k(ρ, σ), and Dpk(ρ, σ). Of course, these approximations are also functions of a, but we hold a constant during each numerical run illustrated here. For the sake of comparison, we include a least squares estimate S(ρ) of the mean slope of the graph of ln C(r,σ) versus ln r on the interval [aρ, ρ]; that is, we find S(ρ) = m, where m and B are the values that minimize the integral
In Figs. 2–5 we display various graphs of Dpk(ρ, σ) and S(ρ) as functions of ρ with a = 0.01 and σ = 0.05.
In Fig. 2 we show the graph of the mean slope estimator S(ρ), the graphs of the original Takens estimators D0k(ρ, σ) (k = 1 or k = 2) as well as the next two Takens estimators D0k(ρ, σ) (k = 3 or k = 4), and the graphs of the first four moment estimators Dp0(ρ, σ) (p = 1, . . . , 4). All estimators approach the expected noise-dominated value of 3 as ρ approaches 0, but the approach to the noise-free value of 2 as the magnitude of ρ increases varies considerably with the choices for p and k in the dimension estimator. What we find in Fig. 2 appears to be typical of such noise-contaminated cases: the mean slope estimator (curve labeled Slope) is qualitatively much worse than any of the extended Takens estimators; then appear the pure Takens estimators D0k(ρ, σ) in reverse order by power (dashed curves), with the k = 1 estimator being the most accurate (curve labeled 1). Finally appear the moment estimators Dp0(ρ, σ) in order by power (solid curves labeled with their p values), with the highest moment estimator D40(ρ, σ) being the most accurate.
In Fig. 3 we show the graphs for the extended estimators Dlk(ρ, σ) corresponding to the first moment with p = 1 and k = 0, 1, ··· , 4. In this figure we see the first sign of a curious phenomenon: as the value of k decreases from 4 to 1, the accuracy of these estimators increases; however, the estimator for the value k = 0 (light solid curve) is not so accurate as the one for the value k = 1 (heavy solid curve). That is, the most accurate estimate is given by the case D11(ρ, σ).
In Fig. 4 we show the graphs of four mixed estimators D2k(ρ, σ) corresponding to the second moment, with p = 2 and k = 0, 2, 3, or 4; the case k = 1 is omitted as there is no solution to the Newton’s iteration for some values of ρ. Again the accuracy increases as k decreases from k = 4 to k = 2; the case k = 0 (light solid curve) is comparable to the case k = 3 (dashed curve labeled 3). Here, the most accurate estimate is given by the case p = k = 2 (heavy solid curve).
The mixed estimators Dpk(ρ, σ) for the third (p = 3) and fourth (p = 4) moments exhibit behavior very similar to that shown in Figs. 3 and 4, and so these results are not shown. The accuracy increases as k decreases, with the most accurate being given by k = 3 when p = 3 and k = 4 when p = 4. Also, as for the case p = 2 and k = 1, there are other combinations of p and k for which the Newton’s iteration does not converge for most values of ρ (e.g., k = 2 and p = 4). Generally, the best results are given by the mixed case p = k, which is summarized in Fig. 5; note the much finer dimension scale of the figure. As can be seen in this figure, the accuracy increases as k decreases, with the most accurate evidently given by D11(ρ, σ) for which p = k = 1.
As noted above, another method for modeling noise in a correlation integral C(r) has been suggested by Ben-Mizrachi et al. (1984), Simm et al. (1987), and Theiler (1986, 1991). Ben-Mizrachi et al. (1984) argue that during empirical observations the error in measuring small distances r overwhelms the C(r) distribution to produce one associated with a random distribution in the embedding space Rc at those small values. Consequently, C(r) scales as rc for small values of r, and as rD, with D the actual correlation dimension, at larger values of r. To represent this behavior using an ad hoc correlation integral C(r, δ), we write
It is clear from (3.5) that the graph of ln C(r, δ) versus ln r consists of two straight lines, one having slope c and the other D, that meet at r = δ. This graph is labeled δ = 0.1 in Fig. 1; clearly, this curve is qualitatively similar to that given by C(r, σ) in (3.3), and so we might expect to obtain similar approximations to the correlation dimension in these two cases. Our investigations of Dpk(ρ, δ) for the case c = 3, D = 2, a = 0.1, and δ = 0.2 (Fig. 6) do in fact reveal results very similar to those illustrated in Figs. 2–5: the best results for noise-contaminated time series are given again by the case D11(ρ, δ) (p = k = 1), with best pure-moment cases Dp0(ρ, δ) given by the largest considered value of p, the best original Takens estimator D0k(ρ, δ) given by the choice k = 1, and the worst overall estimate given by the mean slope S(ρ) (3.4).
b. Sensitivity to lacunarity
As mentioned above, a dimension calculation may fail to converge adequately because superimposed noise distorts the actual correlation integral, or, when noise is minimal, because the actual correlation integral itself may be lacunar. In this case, as shown in Theiler (1988), the original Takens estimators may fail to converge; his arguments apply equally well to the extended ones. From the results of Bingham et al. (1987), it is not hard to show that, in the perfectly noise-free case, convergence of these estimators is equivalent to the existence of limρ→0 C(xρ)/C(ρ) for all 0 ≤ x ≤ 1, which in turn is what we mean by nonlacunarity. However, the estimators may fail to converge by different amounts, and so we may seek those that exhibit the smallest error arising from lacunarity. To seek these estimators, we introduce an ad hoc, periodically lacunar correlation integral CL(r) given by
with aL = 0.5, bL = 1.0, and D = 2. This integral is plotted in Fig. 7; the curve is clearly wavy, making a value for the correlation dimension difficult to estimate from the local slope of the curve. Moreover, this curve resembles those found by Lorenz (1991, Fig. 3a) when he uses time series of a variable that is weakly coupled to the other variables in his dynamical system; as we note above, he does not find this wavy behavior when he uses a strongly coupled variable.
For the correlation integral CL(r), we find that, except with regard to the mean slope method, the numerical results for Dpk(ρ) (Figs. 8 and 9) are opposite of those obtained from the noise-contaminated cases of section 3a. As before, the mean slope method is the least accurate by a considerable margin (long-dashed curve in Fig. 8). But for the original Takens estimators D0k(ρ), the accuracy—as measured by the amplitude of the oscillation about the line D = 2—increases monotonically as k increases from k = 1 to k = 10 (solid curves). In contrast, as shown in Fig. 9, the pure moment cases Dp0(ρ) for p = 1 and the mixed-moment cases Dpk(ρ) for p = k = 1, 2, 3, or 4 (dashed curves) do not perform as well as even the worst of the original Takens estimators D01(ρ) (k = 1). Finally, the extended cases given by k ≠ 0 and p ≠ 0 (not shown) give better estimates as the value of k increases for a fixed value of p, but they are never better than the estimates given by the pure Takens estimators D0k(ρ).
As a tentative result from the numerical tests summarized in these two subsections, we suggest that a signature for lacunarity is that the variation (i.e., apparent error) decreases as k increases with p = 0 or as p decreases with k = 0, while the corresponding signature for noise contamination is that the apparent error decreases as k decreases with p = 0 or as p increases with k = 0, or this error is minimized when p = k. Thus, we might expect that the estimators that are least sensitive to lacunarity are those with p = 0 and k as large as possible, while the one least sensitive to noise contamination is the one with p = k = 1.
4. Hénon and Lorenz attractors
In this section we investigate whether the variations in the dimension estimates found in the previous section typify numerical nonlinear dynamical systems. Natural ones to consider are the Hénon (1976) attractor, which is already thought to be lacunar (Grassberger 1988); the Lorenz (1963) convection attractor, which may be structurally unstable (Williams 1979) and so, in numerical simulation, may amplify round-off error to the level of apparent noise, and the Lorenz (1991) climate attractor, which has a very large correlation dimension value of approximately 15. When estimating distances in this portion of the study, we completely specify the points in the appropriate two-, three-, or 21-dimensional phase space in order to avoid any possible problems associated with lagging and embedding single-variable time series; these other issues are discussed in section 5a. We also discard the initial portions of the trajectories in order to minimize the effects of transients on the dimension estimates. Finally, we use the simpler dimension formulas (2.11)–(2.14), thereby ignoring the small correction summarized in the appendix that arises from accounting for the minimum retained distance rmin.
Typically, and as we do in section 5b, correlation integrals are estimated by calculating distances between points on a single time series. For sufficiently long time series, this procedure requires an unreasonably large amount of computer memory, and so we use a different approach here and in section 5a: we follow a pair of points for a large number of iterations and let the correlation integral C(r) be the empirical probability distribution for the distances between the pair of points after the same number of iterations. It is not hard to show, with mild probabilistic conditions, that this procedure is equivalent to the usual one. The crucial practical consideration is that this procedure does not require a great deal of computer memory to implement, in contrast to the usual approach. One might even expect that the correlation integral given by our procedure would yield the dimension more accurately: if one were to generate, say, 5 × 107 distances with either method, then with the usual method these distances are drawn from a single trajectory only 104 iterations long, while with the present method, these distances are drawn from two trajectories 5 × 107 iterations long, which would seem to explore the attractor more thoroughly. Of course, the present method can be implemented only with computer-generated attractors, for which finding trajectories is relatively inexpensive.
a. Hénon attractor
We first consider the case of the standard Hénon system (Hénon 1976):
for which the parameter values are given by an aH = 1.40 and bH = 0.3. Here two trajectories are followed until 5 × 107 distances less than 1 are obtained; the first 2 × 105 pairs of points are not used, in order to eliminate any transient effects. Then these distances are tabulated using bins of width 0.005.
Applying the extended Takens estimators to this correlation integral C(r), we find in Fig. 10 that, when p = 0, the amplitude of the oscillation in the values of D0k(ρ) decreases as k increases from k = 1 to k = 3 (solid curves labeled with k values), while in the case of Dp0(ρ), the apparent accuracy increases as p decreases from p = 4 to p = 1 (dashed curves). The best estimate for D in this case, 1.20 ± 0.01, is given by the estimator D0k(ρ) for k = 3 in the scaling region 0.05 < ρ < 0.7. This approximation for D is compatible with the value 1.21 ± 0.01 reported by Grassberger and Procaccia (1983a,b) as well as with the value 1.199 ± 0.003 reported by Arneodo et al. (1987), who used 108 points to generate distances.
The changes in the curves of Fig. 10 as the values of p and k are varied match the lacunarity signature that we found in section 3b, and so we conclude that our approximation of the correlation dimension for the Hénon system suffers from apparent lacunarity rather than noise in our estimate of C(r). Indeed, this conclusion is plausible because the time series is computer generated so that the noise level is limited to that attributable to round-off error. This result also agrees with that of Grassberger (1988) when relatively coarse binning of the time series data is used.
b. Lorenz convection attractor
Our second example is the standard Lorenz (1963) convection attractor, given as a solution to the system
with parameters rLor = 28, bLor = 8/3, and σLor = 10. As our estimates for C(r) might be subject to noise contamination arising from undersampling, we consider the behavior of the dimension estimates for an increasingly large number of distances; if the dimension estimates stabilize in a scaling region of ρ as the number of distances increases, then we may conclude that the best possible estimate for D is found in that scaling region. In each case, the temporal differencing scheme implemented by Lorenz (1963) with a time step of 0.005 is used and the first 2 × 105 pairs of points are discarded. Then, only distances less than 5 are retained, as experience has shown that the scaling region tends to occur between ρ = 1 and ρ = 3; in all cases considered, bins of width 0.05 are used.
We first consider the case of 106 distances (Fig. 11a). It is difficult to conclude from this figure and our results in section 3 what the best estimate for D might be. In the region 2.5 < ρ < 3.5 the dimension estimates accumulate between 2.03 and 2.09, with the curve for D0k(ρ) when k = 1 showing the least variation. Mostly the behavior in this figure tells us that we do not have enough distances to obtain a usable estimate for the correlation dimension D.
In Fig. 11b we show the results for the case of 2 × 107 distances. Now we begin to see a recognizable pattern: the dimension estimates D0k(ρ) decrease in value monotonically as k decreases (upper solid curves), the estimates Dp0(ρ) further decrease monotonically as the value of p increases (dashed curves), and the estimates Dpk(ρ) attain their smallest values when p = k = 1 (lower solid curve); this curve is the only one that is level in a scaling region, given approximately by 1.9 < ρ < 2.7. Comparison of Fig. 11b with Figs. 2–5 clearly reveals the signature of dimension estimates contaminated by noise; it might therefore be no surprise that the curve for D11(ρ) gives the best apparent estimate, D(ρ) = 2.06, in this scaling region.
The case of 108 distances is shown in Fig. 11c. The same overall pattern as that found in Fig. 11b is obtained as the values of p and k are varied; here, however the dimension estimates vary less than they do in the previous case since there are five times as many distances binned. Consistent with the results for the smaller number of distances, and in spite of the fact that the time series are computer-generated, noise or undersampling appears to lead to uncertainty in the estimates for D, with the best estimates apparently given by the estimators D11(ρ) and D22(ρ); in the longer approximate scaling region of 0.9 < ρ < 2.4, these estimators again apparently yield the best value for D, D(ρ) = 2.06.
To demonstrate that the dimension estimates appear to converge in a scaling region of ρ as the number of distances increases, we show in Fig. 11d the results given by the best extended estimator D11(ρ) as the number of distances increases from 106 to l08. Clearly, the dimension estimates converge to D(ρ) = 2.06 ± 0.005 in the scaling region 1 < ρ < 2.5; moreover, 5 × 107 distances (equivalent to 1 × 104 points on a single time series) appear to be adequate for producing this estimate. Somewhat reassuringly, this estimate for D agrees with the usual one of 2.05 ± 0.01, which is given for example by Grassberger and Procaccia (1983a,b) and Nese et al. (1987).
To explain the apparent result in this section that stochastic noise might significantly contaminate a computer-generated numerical solution to a deterministic dynamical system, we offer two scenarios. First, it may be that for any of the estimators, the convergence as the trajectory lengthens may be exceptionally slow. This possibility is supported by the dramatic changes in the estimators that are revealed by comparing Figs. 11a–c as the pairs of trajectories are lengthened from those having 106 distances less than 5 to those having 2 × 107 distances to those having 108 distances. Here, even very long trajectories apparently suffer from undersampling at small distances less than 1 or so; this undersampling may then apparently produce a stochastic error in the estimates for these relatively small distances.
For the second scenario, we recall the geometric model of the Lorenz attractor constructed by Williams (1979). This dynamical system has the property that as the value of a control parameter is varied, the corresponding attractor varies through an uncountable family of dynamically inequivalent attractors. This phenomenon suggests that the Lorenz attractor is not structurally stable; in that case, round-off errors may in effect change the dynamical system being integrated into a nearby, but inequivalent, one. Because the trajectory is now approaching a slightly different attractor, the long-run effect is to magnify the original round-off error to produce a more noticeable effect; that is, the “attractor jumping” that results from small round-off errors may succeed in amplifying these very small stochastic errors to the level of detectable noise, which would be most apparent at small distances.
c. Lorenz climate attractor
As one final test of our correlation dimension algorithm, we consider an attractor of the 21-coefficient model studied by Lorenz (1991). For the case of the coupling parameter c = 0.1, the attractor has a Lyapunov dimension of 17. Although Lorenz does not find correlation dimension estimates using distances in the 21-dimensional phase space of the model, he does report a value of approximately 15 when he lags and embeds times series composed of 4 × 103 points of a function of the variables that is strongly coupled to all of the model variables. We thus take the value of 15 to be the one we seek from our algorithm.
Lorenz (1991) obtains his estimate of the correlation dimension by using a two-point slope method. In the results presented below we obtain a two-point slope estimate Ds(ρi) for each distance ρi given by
If our dimension estimators (2.8)–(2.10) are less sensitive to errors than is the two-point slope estimator Ds(ρi), then we should obtain estimates that are larger than those given by the two-point slope method but that are less than or equal to the Lyapunov dimension value of 17.
In our analysis we follow Lorenz (1991) by using a fourth-order Runge–Kutta scheme with a time step of 0.025 and by setting all other parameters to the values he uses for the case c = 0.1. The distance data we use are obtained by finding at each time step the distances in 21-dimensional phase space between two trajectories. For the first trajectory, all the variables except the first two are set to zero initially; these are set to X1(0) = 2 and Y1(0) = 1, respectively. The initial conditions for the second trajectory are found by adding 1 × 10−9 to the value of each variable in the first initial condition. The first 1 × 107 iterations of each trajectory are discarded as transients. We use every tenth point from the next 3 × 109 iterations of each trajectory and so generate 3 × 108 distances, which is the number of distances given by approximately 2.4 × 104 points from a single time series; this number of distances proves adequate for the Lorenz convection attractor discussed in section 4b. We tabulate the data using 200 bins of width 0.05. Finally, we use the formulas given by (2.11) and (2.14).
Our results are shown in Fig. 12. The two-point slope estimator Ds(ρi) (4.3), and the usual extended estimators D11(ρ), D40(ρ), D10(ρ), D01(ρ), and D02(ρ) are shown; D03(ρ) is not shown as its largest value is only 7.5. Here Ds(ρi), given by the double curve in the figure, is approximately 14.5 in a scaling region 2.4 < ρ < 2.7; this value of D is a bit smaller than that reported by Lorenz (1991). The remaining curves follow the order found above for undersampled attractors: the best estimator is D11(ρ), which has a scaling region similar to that for Ds(ρi) but which has a slightly larger value of 14.9 in this scaling region. This is the closest value to 17 in the set shown, although other estimators, particularly those in the family D1k(ρ), give values as large as 15.5 for large values of k. The next largest values are given by the estimators Dp0(ρ), with D40(ρ) ≥ D10(ρ) (dashed curves in Fig. 12). The estimators having the greatest errors are the pure Takens ones D0k(ρ), with the error increasing with k. We conclude that our method gives reliable results even for attractors having very large correlation dimension values; as with the Lorenz convection attractor, the best estimate is given by the estimator D11(ρ) when the attractor is undersampled at small distances.
5. PBL horizontal winds
In this section we apply the extended Takens estimators to horizontal wind observations taken in the boundary layer. Because we must now lag and embed the single observed time series, we must consider the issues of appropriate variable choice and embedding dimension. We first examine the dependence of the various estimators on the embedding dimension choice by using time series from the Lorenz convection attractor. Then we use these results to determine the best estimate for the dimension of the low-frequency PBL time series. We emphasize that the primary purpose of this part of the study is to demonstrate that the best estimate of the dimension for any given atmospheric time series is obtained by studying the results given by the family of extended Takens estimators. We are thus careful not to claim that the dimension value we obtain is an accurate estimate of the dimension of the underlying attractor, as we likely do not have enough time series points to justifiably make such a claim. Preliminary results are discussed in Fosmire et al.(1995).
a. Embedding dimension dependence
To obtain a guess for the required phase space embedding dimension M, investigators must use model reconstruction, in which they create several data vectors of greater and greater length m by lagging a single atmospheric series m times (e.g., Henderson and Wells 1988). Then they determine how the estimate for D varies with the vector length m. When the mean slope method for estimating the value of D is used, the magnitude of M is taken to be the smallest value of m for which the estimate of D stops changing.
In order to determine how the dimension estimates vary with the choice of embedding dimension guess m, we consider two time series of (4.2) for the standard Lorenz (1963) convection attractor, for which the correlation dimension D is found in the previous section to be 2.06 when the values of p and k are chosen correctly and an adequate number of time series values are used. Consistent with the analysis in the previous section, distances are obtained by simply calculating the Euclidean distance between the X components of these two trajectories at a fixed time. The parameter values are still rLor = 28, bLor = 8/3, and σLor = 10; the time step is the same, 0.005; and the same number of initial points, 2 × 105, are discarded from each trajectory.
For the case of embedding dimension guess m, we lag the X components of the time series by m time steps to produce a vector of length m. We normalize the distance values by m1/2 in order to keep distances between corresponding points approximately equal; we retain only the distance values that are less than 10, and we sample until we obtain 107 such distances for the case m = 10. For other embedding dimension guesses, we use this number of time steps, although in general the number of distance values used to calculate the expected values decreases slightly with m. We note that this procedure for obtaining distances requires large embedding dimension guesses before convergence of the results to a dimension value occurred.
As for the calculation using the full, unlagged time series (Figs. 11b,c), the best estimate is given by D11(ρ) for which p = k = 1 (not shown). In Fig. 13 we give a plot of the values of D11(ρ) as functions of the distance ρ for four embedding dimension guesses m. The dependence of D11(ρ) on the embedding dimension follows that typically found by others (Henderson and Wells 1988): when the magnitude of m is too small, the dimension estimate is also too small. In this case, the dimension curves also do not exhibit an unambiguous scaling region, with a general monotonic decrease in the dimension values as the value of ρ increases. In contrast, when the embedding dimension guesses m are increased to 30 and 40, a scaling region giving values for D11(ρ) near 2.06 can be clearly seen. The dimension estimates show some fluctuations, however, in contrast to those seen in Figs. 11b,c. Based on the results in Fig. 13, we may estimate the correlation dimension to be 2.06 ± 0.02, in agreement with the results of other studies (e.g., Grassberger and Procaccia 1983b; Nese et al. 1987) and the results in Fig. 11. In addition, a clear signature of an adequate embedding dimension can be found by examining the curves in Fig. 13: the monotonic decrease in dimension estimates becomes replaced with a wavy curve defining the scaling region and giving the correct estimate for D.
We next apply these results to boundary layer observations.
b. Wind observations
The wind measurements used here were obtained during the Appalachian Mountain Air Chemistry Experiment, AMACEing-88, held in July and August 1988 at the Black Moshanon field site at the Mid-State Airport adjacent to the Black Moshanon State Park about nine miles northwest of State College, Pennsylvania. The Mid-State Airport is located on a relatively flat Appalachian plateau. There is a significant drop in mean elevation about 2.5 km to the southeast. The elevation also drops, although not as markedly, about 2.5 km to the northwest. Black Moshanon Lake is just to the northeast of the airport. The Black Moshanon area may be generally categorized as having dendritic, complex terrain.
Data collection using an acoustic Doppler sodar began on 20 July at 2110 EDT and ended at 1810 EDT on 23 August. Owing to a software error, no data were collected between 1327 EDT on 8 August and 1137 EDT on 22 August. The data were archived every 15 min and represent 15-min averages of the winds from 50 m to 800 m in 50-m increments. In this study, the 100 m and 200 m U and V components of the wind were used, where U is the component along the strike of the nearby Appalachian ridge (from approximately 240°) and V is at right angles, from approximately 150°. The mean values of the U and V components at 100 m are 1.35 and 0.06 m s−1, with standard deviations of 1.93 and 2.69 m s−1, respectively, and at 200 m they are 2.43 and −0.14 m s−1, with standard deviations of 2.85 and 4.13 m s−1, respectively. Thus, the wind direction was typically along the strike of the Appalachian ridge, as expected.
Before the value of the correlation dimension D is estimated, the gaps in the time series are removed. In our case, the majority of the gaps involved one data point, with the largest gaps occurring near the beginning and end of the time series. These parts near the beginning and end of the series are not used. For the data in between, gaps of up to three consecutive points are augmented via linear interpolation; for gaps greater than three points, the gap is removed by concatenating the broken series together (Keppene and Nicolis 1989). Both of these changes to the dataset are made in only a few places, and as a result they do not appreciably affect the correlation dimension estimates. The resulting time series has 1795 points (approximately 6 × 106 distances) representing 15-min averages of the horizontal wind in the lower boundary layer.
Next, both the optimal lag for the model reconstruction step and the best time series variable must be determined. Following many investigators (e.g., Henderson and Wells 1988), we use the first zero crossing in the autocorrelation graph to determine the smallest lag yielding independent data. As advised by Lorenz (1991), we also seek the variable most closely coupled to the underlying dynamics on the attractor. In our analysis, we first estimate the correlation dimension values by calculating D10(ρ) separately for each wind component at each level, reasoning after Takens (1981) that any one time series should give us an adequate estimate of the correlation dimension for this low-frequency boundary layer flow. We find that different lags are needed for each time series (150 quarter hours for U component U100 at 100 m, 130 quarter hours for V100, 194 for U200, and 226 for V200), and that the values of D10(ρ) are markedly different, ranging from 2.6 for V200 to 3.2 for U200 and V100 to 3.9 for U100. Poveda-Jaramillo and Puente (1993) also obtain a range of values for the correlation dimension—between 4 and 7—when they use 24 30-min time series obtained from different wind components measured at different heights. Clearly, this variation calls into question whether any of these individual series are well coupled to the underlying dynamics (Lorenz 1991).
We now choose to combine all four variables into a single, more canonical variable, the specific horizontal kinetic energy, K = 2−1(U2100 + U2200 + V2100 + V2200). For the time series of K, the first zero crossing lag is 116. Values of Dp0(ρ), D0k(ρ), and Dpk(ρ) for p ≤ 4 and k ≤ 4 are calculated using embedding dimension guesses m up to 9, where the distances are normalized by m1/2, as for the Lorenz convection attractor discussed in section 5a. As before, we find that the best estimates for the Black Moshanon case are the same as those for undersampled signals: those given by the estimator D11(ρ) for which p = k = 1 or D22(ρ) for which p = k = 2. In Fig. 14a the usual family of extended estimators is shown for the Black Moshanon time series with m = 6; these curves clearly exhibit the same qualitative behavior as that found for the Lorenz convection attractor in Figs. 11b,c. As Fig. 14b reveals, the embedding dimension guess m = 6 seems to be the best choice: as with the model-reconstructed Lorenz convection attractor in Fig. 13, embedding dimension guesses that are too small produce a monotonically decreasing estimate for D, and embedding dimension guesses that are too large produce a wavy curve. From the results in Fig. 14b for m = 6, we estimate the correlation dimension for this low-frequency Black Moshanon data to be 3.9 ± 0.1; the value 3.9 is indicated by the dashed line. As we might expect, this estimate for D is close to that given by the 100 m U component U100.
These results are consistent with those of Lorenz (1991), who found that analysis of a single time series tended to underestimate the dimension value, but that the most strongly coupled variable tended to give the largest and closest-to-correct value. The largest value from the individual Black Moshanon series is given by the analysis of U100, indicating that it and the specific kinetic energy are about equally coupled to the dynamics of the attractor. Other evidence that boundary layer variables may be coupled differently to the underlying attractor is provided by Poveda-Jaramillo and Puente (1993), who examined 24 different 18 000-point, 30-min time series and found dimension values in the range 4–7. Perhaps surprisingly, the value 3.9 ± 0.1 we obtain here is similar to that found for much higher-frequency boundary layer flows by Henderson and Wells (1988) and Orcutt and Arritt (1995) who analyzed, respectively, 1500 points from 10 min and 108 000 points from 3 h of wind data that were measured at different times on the BAO tower. However, Tsonis and Elsner (1989) used 4000 points of different BAO data spanning 11 h and found a much larger estimate for D.
In contrast to the above results, Weber et al. (1995) used a million-point time series of 21-Hz sonic anemometer data spanning an entire day to conclude that there was no low-dimensional attractor characterizing the unstable boundary layer. However, their use of high-frequency data spanning a complete diurnal cycle may mean that they were sampling more than one boundary layer regime and hence that their dimension estimates could not converge to a well-defined value. This sort of behavior has been seen by other investigators; for example, Fraedrich (1986) found no apparent correlation dimension value when he used atmospheric data spanning entire years, but he did find such low values when he partitioned the data by season. Moreover, such rapid sampling may well record high-dimensional noise that masks the actual dimension of the attractor. It is clear that much further work is needed to determine the best means for sampling atmospheric data in order to ensure that the data adequately describe only one coherent structure type and hence only one underlying attractor.
6. Concluding remarks
In this article, we apply our new extended Takens (1985) algorithms to estimate the correlation dimension D of time series contaminated by errors such as those encountered when analyzing atmospheric data. These algorithms are very simple to implement for geophysical time series, as only tabulation of the weighted means in (2.11)–(2.14) is necessary, which is a relatively fast operation that requires little computer memory. Consequently, we eliminate the need to order by magnitude the set of distances, as is often done in the classic slope procedure (e.g., Albano et al. 1988). Moreover, we neither differentiate nor integrate the function C(r), but rather take appropriate means to find directly our approximations Dpk(ρ) of D. In addition, determining the best dimension estimate is aided significantly by the fact that each estimate exhibits behavior that may be identified with the source of contamination. Two major types of irregularities may occur: (i) at small distances r, stochastic noise arising from undersampling of the attractor may be present; or (ii) the correlation integral C(r) may be lacunar; that is, the ratio C(xρ)/C(ρ) may fail to converge as ρ approaches zero (Theiler 1988). The first type of irregularity is unavoidably present in any measured time series, while the second is a possible property of the attractor. Grassberger (1988), for example, shows numerically that the Hénon attractor may well be lacunar, although Theiler (1988) believes that it is not. The collective behavior of our estimators can be used to detect the occurrence of the two major types of irregularities and to produce the best estimate in each case. Finally, irregularities associated with under- or overembedding can also be identified and handled using the same estimators.
To determine how the estimates Dpk(ρ) given by the extended Takens estimators vary with k and p in the conditional expected value E((r/ρ)p|ln(r/ρ)|k : r < ρ) as well as how they compare with the estimates given by the mean slope method, we consider two ad hoc correlation integrals contaminated by noise and one correlation integral that is periodically lacunar. In all cases, the mean slope method exhibits the greatest sensitivity to problems in the time series; any value of k and p gives a better result than does the slope method. When we compare the results given by the extended estimators themselves, we find a significantly different dependence of the estimates on k and p for each type of time series problem: when the time series is dominated by stochastic noise at small distances, the estimates given by Dpp(ρ), in particular D11(ρ), provide the results least sensitive to noise; those given by Dp0(ρ) provide the next best, improving as p increases in value; and those given by D0k(ρ) provide the most sensitive results, worsening as k increases in value. Apparently, the logarithmic weighting determined by the value of k enhances the contribution of the noise when there is no weighting by powers of r; however, we are baffled as to why the best result is given by the case p = k. In contrast, when the correlation integral C(r) is periodically lacunar, the opposite sensitivity is seen: the best results are given by the case D0k(ρ), and they improve as the value of k increases; the next best are given by Dp0(ρ), with only small differences occurring as the value of p is changed; finally, the worst is given by Dpp(ρ). Clearly, when analyzing a time series of unknown characteristics, a large number of estimates Dpk(ρ) for different values of p and k may be required to determine the optimal estimate for D. As noted previously, we have only numerical evidence for our statements about sensitivity, and we have no analytic proof or explanation.
We next investigate whether the above results might apply to the Hénon (1976), Lorenz (1963) convection, and Lorenz (1991) climate attractors. For the Hénon attractor, the noise and lacunarity signatures found in the idealized cases indicate that the indeterminacies in the dimension arise more from lacunarity, or apparent lacunarity, than from noise in the system, just as we would expect from the fact that a computer model should have minimal noise contamination (Grassberger 1988). The optimal estimate for the correlation dimension is given by the estimator that shows the least variation in an interval of ρ known as the scaling region. As expected from the results of Grassberger (1988), lacunarity appears to dominate the Hénon results; the best estimate of D, D(ρ) = 1.20 ± 0.01, is given by the case D03(ρ) that typifies lacunar signals.
For the Lorenz (1963) convection system we also examine the behavior of the estimators as a function of the number of distances used to approximate the conditional expectations. Curiously, in spite of the fact that the time series is computer generated, in this case the uncertainty in the dimension estimates appears to be more related to noise than to apparent lacunarity. A signature for undersampling can be identified as well; when too few distances are used, none of the estimators Dpk(ρ) are approximately independent of ρ in any interval and so there is no scaling region. However, once enough distances are included in the approximation of the conditional expectations, a scaling region appears only for the case p = k. Thus, the signal of a noise-contaminated or undersampled time series is found, with the best estimate of D, D(ρ) = 2.06 ± 0.005, given by the case p = k = 1.
The last test attractor we study is the Lorenz (1991) climate attractor whose correlation dimension value is quite large, near 15. We again find that the best estimator is given by D11(ρ), and we indeed obtain a value of 14.9 ± 0.1, thereby showing that our method is not limited to low-dimensional cases.
To determine the proper embedding dimension, we study the model-reconstructed Lorenz convection attractor and find that underembedding is revealed by a monotonic decrease in the dimension estimate as ρ increases, while overembedding is revealed by a wavy behavior in the dimension estimate as ρ varies. The best embedding dimension is given by the occurrence of a scaling region wherein the dimension estimate varies only a little.
We next estimate the correlation dimension D for the low-frequency horizontal wind data measured at two levels (100 and 200 m) in the lower portion of the boundary layer at the Black Moshanon field site. As for the noise-dominated case using the entire series, the estimate D11(ρ) for the values p = k = 1 is best. Consistent with the results for the model-reconstructed Lorenz convection attractor, we find that the best embedding dimension guess, m = 6, is given by a scaling region for a range of distance values. As recommended by Lorenz (1991), we seek the variable that is the most strongly coupled to the attractor dynamics to obtain the best possible estimate of D from these boundary layer observations. We suggest that the specific kinetic energy K = 2−1(U2100 + U2200 + V2100 + V2200) might be such a variable, and we find that the largest estimate of the correlation dimension, 3.9 ± 0.1, is given by this variable. Not surprisingly, the individual wind components at the two levels give values that are less than or equal to this estimate. When investigating the chaotic properties of an atmospheric time series, we thus recommend that many variables be combined into a single one, such as the specific kinetic energy. We are currently applying these ideas to other measured and simulated boundary layer time series (e.g., Winstead et al. 1995).
We thank Dr. Dennis W. Thomson for making the PBL data available to us, Dr. Harry W. Henderson and two anonymous reviewers for their comments on an earlier draft of this manuscript, and Nathaniel S. Winstead for performing some of the numerical estimates of the attractor dimensions. This research was supported by the Office of Naval Research through Grants N00014-90-J-4012, N00014-93-1-0252, and N00014-96-1-0375.
Several Extended Takens Estimators
Here we list the formulas used to improve the estimates for the correlation dimension when we account for the unretained distances in the interval (0, aρ), where aρ ≥ rmin. By dropping the limit as ρ → 0, replacing D with Dpk(ρ), using aρ as the lower limit in the integral, and substituting (r/ρ)D for C(r) in (2.7), we obtain the following formulas with the aid of the symbolic manipulator Derive®. For brevity, we write Dpk instead of Dpk(ρ) in (A.2)–(A.6):
Current affiliation: Pacific NW Laboratories, Richland, Washington.
Corresponding author address: Dr. Hampton N. Shirer, Department of Meteorology, The Pennsylvania State University, 503 Walker Building, University Park, PA 16802-5013.
Current affiliation: Bell Geospace Inc., Lamont-Doherty Earth Observatory, Columbia University, Palisades, New York.