adg vol5 n02 Toz 279_287.pdf ANNALS OF GEOPHYSICS, VOL. 45, N. 2, April 2002 279 Detecting geomagnetic field nonlinearities by bispectral analysis and a phase coupling nonlinear technique Roberta Tozzi and Angelo De Santis Istituto Nazionale di Geofisica e Vulcanologia, Roma, Italy Abstract The Earth’s magnetic field varies over a wide range of characteristic times, say, from years to centuries, and more. In order to detect some nonlinear features of the geomagnetic field evolution we first apply a nonlinear spectral technique, i.e. bispectral analysis, to the secular variation of Hartland Geomagnetic Observatory (U.K.). Then, due to difficulties of bispectral analysis inherent in the shortness of data, we introduce a simpler, but more efficient, technique called Spectral Phase Analysis for Quadratic Coupling Estimation (SPAQCE). Both nonlinear spectral methods here applied are based on the presence of a phase (sum and difference) coupling in case of quadratic interactions between two constituent components of the physical system underlying the generation of the geomagnetic field. The application of SPAQCE to annual means of a U.K. combined geomagnetic time series allows us to discriminate nonlinear interactions between couples of characteristic times (specifically 5-6, 5-8 and 2-26 years) of the field generated in the outer core. 1. Introduction It is well known that to study a Gaussian process and to completely describe its statisti- cal properties, it is sufficient to use simple mathematical tools such as the autocorrelation function or the power spectrum. However, these methods are phase blind so, a general random additive process (for instance a coloured noise) cannot be completely discriminated from a random multiplicative process (such as non- linear cascade process) or a deterministic non- linear process (e.g., a chaotic phenomenon) throughout second order statistics or con- ventional Fourier power spectrum. In the past, this lack of appropriate analytical tools has led scientists to treat non-Gaussian pro-cesses (i.e. the greatest number of natural processes), as if they were Gaussian. Nevertheless, in recent decades, some powerful techniques have been developed to study nonlinear and non-Gaussian phenomena. Among them, an important role is played by the Higher Order Spectral Analysis (HOSA) or, equivalently, Higher Order Statistics (HOS) (e.g., Mendel, 1991). HOSA is an extension of second order spectral analysis (i.e. conventional Fourier analysis) to higher orders and is a technique to estimate higher order power spectra, the so-called polyspectra (Subba Rao and Gabr, 1984). Polyspectra discriminate nonlinear processes (i.e. processes characterised Mailing address: Dr. Roberta Tozzi, Istituto Nazionale di Geofisica e Vulcanologia, Via di Vigna Murata 605, 00143 Roma, Italy; e-mail: tozzir@ingv.it Key words bispectral analysis – nonlinearity – geomagnetic field secular variation 280 Roberta Tozzi and Angelo De Santis by a nonlinear relation between any perturbing input and the corresponding output of the system) because of their aptitude to reveal not only amplitude information, but also phase information. For instance, two random processes such as a Gaussian and a non-Gaussian process can be distinguished only by evaluating their third order spectrum (i.e. the bispectrum). In practice, considering that the bispectrum of a signal can be viewed in terms of a decomposition of the third statistical moment over frequency, it is easy to understand that the bispectrum of a Gaussian process is identically zero (e.g., fig. 1b). On the contrary, that related to a non-Gaussian process contains very interesting information such as the presence of some quadratic coupling. In this paper we introduce the basic concepts of bispectral analysis which, evidently, is the simpler case of HOSA together with its attractive application to a situation of geophysical interest, namely that of the geomagnetic field Secular Variation (SV). Due to the shortness of available geomagnetic data, the results from the ap- plication of the HOSA are not always reliable. In order to give more robustness to our results we will introduce a simpler, but more efficient technique that allows a sharper and even more consistent discrimination of any quadratic coupling. In the light of these considerations, the intention of this work is to open a new way to verify and even to give a quantitative description of the geomagnetic field SV nonlinear prop- erties. The main interest in this kind of analysis lies in the importance it would have for a better comprehension of the underlying dynamo processes that generate and sustain the field. Of course, nothing excludes the possibility to apply the new technique to other kinds of geophysical data, since typically most of them are very short time series. The next sections of this paper are arranged in the following way. The second section defines bispectral analysis and describes the main aspects of the theory. The third section illustrates the application of this methodology to synthetic as well as to experimental geomagnetic data, namely the SV data deduced as annual means first differences; then it shows a simpler alternative to bispectral analysis in the case of short time series, as is the case of SV data. Finally, in the last section we discuss the results provided by this technique. 2. Bispectral analysis Bispectral analysis is a powerful math- ematical tool to reveal the presence of quadratic couplings between two harmonics (different or not) each associated with a given frequency. From introductory considerations it follows that the bispectrum provides interesting information only if the third moment (or skewness) of a time series is non-zero. Recal-ling that skewness is connected to the asymmetry of the probability density function, we deduce that the bispectrum is a helpful tool to analyse only systems with asymmetric nonlinearities. If this condition is not satisfied, it is necessary to verify successive higher order spectra (i.e. the trispectrum); but higher order spectra are not here considered as the error of spectral estimation increases with the order of spectral analysis. Auto-bispectrum B(f i , f j ) of a stationary time series (e.g., Nayfeh and Balachandran, 1995) defined in time interval T corresponds to a sort of two dimensional Fourier spectrum and it is used to check out interactions between pairs of harmonics within a single time series; it is expressed as follows (2.1) where X(f) represents the Fourier amplitude spectrum of the time series x(t) at frequency f and X * (f) is the complex conjugate; the < > brackets stand for an ensemble averaging over many statistically similar realisations of the process under scrutiny; in general, the bispectrum (from this point onwards we will omit the preceding term «auto») defined in eq. (2.1) is a complex quantity. Evidently, B f fi j ,( ) = T X f f X f X f T i j i j * *= < +( ) ( ) ( ) >lim 1 281 Detecting geomagnetic field nonlinearities by bispectral analysis and a phase coupling nonlinear technique evaluation of X( f ) makes use of the Finite Fourier Transform algorithm, which takes into account the narrowness of the analysed series; in fact, when working with data coming from experimental measurements, we usually know x(t) only for a finite time T. In practical applications, it is customary to work only with the magnitude of the bispectrum that gives the strength of quadratic inter- action, so graphic representation of the bispectrum is essentially a two-dimensional function in terms of an f i -f j diagram. It is quite evident that the corresponding contours contain a high degree of symmetry owing to symmetrical properties of the bispectrum defined in (2.1). It is possible to demonstrate (Chandran and Elgar, 1994) that the knowledge of the bispectrum in the principal domain (or non-redundant area), cor- responding to the triangular region {f i , f j > 0}, is enough for its complete description. One of the methods to evaluate B(f i , f j ) is the so- called direct method based on a segment averaging approach. Some studies (Dalle Molle and Hinich, 1989) have shown that to calculate the nth order spectrum, the length of each segment should be about the (n-1)th root of the time series size. Collis et al. (1998) showed that the variance associated with the bispectrum is strictly related to the power spectrum properties of the signal; for this reason instead of the bispectrum of eq. (2.1), another spectral function is used to detect more punctually quadratic inter- actions, the so-called auto-bicoherence or simply bicoherence. This function, analogous to the bispectrum from which it derives, is rated (normalised) to the power at each contributing frequency; thus, the bicoherence varies between zero and one, just like the square of the correlation coefficient in linear analysis. The main difference is that bicoherence is useful to detect Quadratic Phase Coupling (QPC); the latter phenomenon occurs when two distinct harmonic components in a signal are nonlinearly related so that there is a phase correlation between components at frequencies f i and f j (MacDonald, 1989). Let us define the auto-bicoherence b2(f i , f j ) as the function that measures the fraction of power at frequency values of f = f i ± f j that is (2.2) From eq. (2.2), it is quite evident that, as already mentioned, bicoherence is a sort of normalised bispectrum and, in practice, it expresses fre- quency dependence of correlation. If b2( f i , f j ) = 1, there is a complete quadratic coupling between frequency components at fi and fj; if b 2( f i , f j ) = 0 there is no coupling at all while intermediate values of bicoherence correspond to partial quadratic coupling. 2.1. Nonlinear coupling in synthetic data In order to verify the potentiality of bispectral analysis, we estimated the bicoherence for some synthetic data, properly generated with known linear and/or nonlinear quadratic couplings. Consequently, we studied the time series representing white noise, nonlinear and chaotic sequences (i.e. Lorentz, Henon and logistic maps) and, of course, also periodic signals. These tests (here not completely presented) disclosed that this technique is really a good method not only to discriminate linear from nonlinear processes but also to estimate the periods of possible quadratically interacting modes. Bicoherences obtained from white noise (fig. 1a) are practically flat (fig. 1b). One of these experiments involved a very simple function of time t in terms of a linear combination of two sine functions I(t) = a1sin(2 f1t) + a2sin(2 f2t) (2.3) where a 1 and a 2 are two constant coefficients and f 1 and f 2 are two different given frequencies (we chose 0.64 and 1.11 Hz respectively). The function defined in (2.3) can be viewed as an input to a given system. If this system operates nonlinearly, its output contains nonlinear terms; b f fi j,( ) =2 B f f X f f X f X f i j i j i j , = ( ) < +( ) >< ( ) ( ) > 2 2 2 282 Roberta Tozzi and Angelo De Santis -0.4 -0.2 0.0 0.2 0.4 -0.4 -0.2 0.0 0.2 0.4 WHITE NOISE f 1 (Hz) f 2 ( H z ) 0.1750 -- 0.2000 0.1500 -- 0.1750 0.1250 -- 0.1500 0.1000 -- 0.1250 0.0750 -- 0.1000 0.0500 -- 0.0750 0.02500 -- 0.0500 0 -- 0.02500 0 25 50 75 100 125 150 175 200 -3 -2 -1 0 1 2 3 W hi te n oi se t (s) Fig. 1a,b. a) White noise (arbitrary amplitudes); b) the bicoherence of a white noise is a practically flat two- dimensional function. a b Fig. 2a,b. a) A period of the function I(t)2; b) bicoherence of the function I(t)2 obtained from I(t) defined in eq. (2.3) with a 1 = a 2 = 1, f 1 = 0.64 Hz and f 2 = 1.11 Hz. This is the case of a short data set (128 values); it shows that bicoherence gives a good estimation of the actual coupling frequencies although it is characterised by a large variance. -4 -2 0 2 4 -4 -2 0 2 4 (0.77,1.15)Hz Bicoherence (128 data) f 1 (Hz) f 2 ( H z ) 0.1225 -- 0.1400 0.1050 -- 0.1225 0.0875 -- 0.1050 0.0700 -- 0.0875 0.0525 -- 0.0700 0.0350 -- 0.0525 0.01750 -- 0.0350 0 -- 0.01750 0 10 20 30 40 50 60 0 1 2 3 4 I ( t )2 t (s)a b 283 Detecting geomagnetic field nonlinearities by bispectral analysis and a phase coupling nonlinear technique for example if the system squares the input, then the output, O(t) = I 2(t), includes three quadratic terms. A conventional (i.e. second order) power spectrum of the function I 2(t) shows four peaks (corresponding to the two interacting frequencies f1, f2 together with their sum and difference frequencies f 1 ± f 2 ), while the plot of bicoherence displays clearly the quadratic coupling between the two modes with frequencies f 1 and f 2 appearing in eq. (2.3); this result is illustrated in figs. 2b and 3 for the cases of two data sets of 128 and 1024 values, respectively. The com- parison of the results from the two analyses shows that bicoherence indeed reveals nonlinear couplings but that the length of the time series is a very critical point. In fact, in the former case the relative error is around 20%, while using more data it decreases to about 2%. This is the reason why in the next paragraph we will propose a method that is particularly suitable for scanty series. Of course, this is only a simple example to demonstrate that through the observation of bicoherence, quadratic phase couplings can be discriminated. The application of the method to synthetic data of different origin has shown that the bicoherence of signals associated with linear processes is always flat, whereas signals characterising nonlinear processes may produce both flat and non-flat bicoherence depending only on the presence of quadratic couplings. This means that a flat bicoherence may come from the signal analysis of both linear and nonlinear processes, but if bicoherence is not flat then this is necessarily due to the presence of quadratic phase couplings in the signal, typical of a non- linear process. -4 -2 0 2 4 -4 -2 0 2 4 (0.63, 1.12) Hz Bicoherence (1024 data) f 1 (Hz) f 2 ( H z ) 0.1225 -- 0.1400 0.1050 -- 0.1225 0.0875 -- 0.1050 0.0700 -- 0.0875 0.0525 -- 0.0700 0.0350 -- 0.0525 0.01750 -- 0.0350 0 -- 0.01750 Fig. 3. Bicoherence of the function I(t)2 obtained from I(t) defined in eq. (2.3) with a1 = a2 = 1, f1 = 0.64 Hz and f2 = 1.11 Hz. This is the case of a longer data set (1024 values); it shows that in this situation bicoherence gives a better estimation of the quadratic frequency coupling with a smaller variance than the case of fig. 2b. 284 Roberta Tozzi and Angelo De Santis 2.2. Spectral phase analysis for quadratic coupling estimation Synthetic tests showed that bispectral analysis might become inaccurate when analysing short time series, therefore to gain in faithfulness we developed a different spectral method. In the following, we will call it Spectral Phase Analysis for Quadratic Coupling Estimation or briefly SPAQCE, and use it essentially to verify information provided by HOSA. This is a nonlinear technique in the common sense that it is a technique suitable to detect nonlinear properties of time series. SPAQCE is a very direct and, in some ways, crude method that mainly applies the basic concept underlying QPC. We know that a QPC occurs when two characteristic processes are coupled by a nonlinear interaction producing a new signal characterised by a phase that is simply the sum and the difference of the phases of the original signals. In practice, we implemented a Fortran algorithm able to compare, for each pair of frequencies, the sum and the difference of their phases with the phases corresponding to the sum and difference of these two starting frequencies, this algorithm can also keep note of the corresponding spectral amplitudes. An example of a quadratic coupling in a short synthetic time series (composed of 128 data) produced by multiplying two sines with slightly different periods (5.5 and 8.5 years) is illustrated in fig. 4. Our software returns two series of data: one for the sum of the phases (white circles) and one for the difference (black squares). Only points common to both groups of data correspond to a quadratic coupling. An appropriate choice of a threshold in the spectral amplitudes discloses only the most significant phase coincidences. Practically this threshold is an appropriate percentage of all spectral amplitudes as appears in figs. 4 and 7. As fig. 4 shows, this method gives good results even for 0.0 0.1 0.2 0.3 0.4 0.5 0.0 0.1 0.2 0.3 0.4 0.5 SYNTHETIC SIGNAL Spectral Amplitude Threshold = 10% (9.2,5.5) year sum difference f 2 ( y e a r- 1 ) f 1 (year -1 ) Fig. 4. The method of SPAQCE applied to a short synthetic signal composed of 128 data and consisting of a product of two sines with periods of 5.5 and 8.5 years respectively. This figure evidences the quadratic phase coupling in terms of overlapping between black squares (representing couplings in the differences of couples of frequencies) and white circles (representing couplings in the sums of couples of frequencies). This diagram is based on the largest 10% spectral amplitudes. 285 Detecting geomagnetic field nonlinearities by bispectral analysis and a phase coupling nonlinear technique short time series and the relative error associated with this kind of procedure decreases to around 8% against 20% of HOSA. 3. The case of geomagnetic field From magnetohydrodynamic equations, it follows that the geomagnetic field enjoys deterministic and nonlinear properties. As a matter of fact, some recent studies (e.g., Barraclough and De Santis, 1997; Hongre et al., 1999) have shown some evidence that geomagnetic data contain clear nonlinearities. In the preceding section, we observed that HOSA is a very robust method to detect nonlinearities but it needs many data, which is not the case of annual mean datasets. Nevertheless, we can apply the technique of SPAQCE; with this kind of analysis we are, in principle, able to observe nonlinear couplings between the different harmonics components into which a signal could be split. Therefore, we have estimated the bicoherence and applied the SPAQCE technique as a confirmation to HOSA results to geomagnetic data. Analysed data consists of annual means of the Y component of the secular variation of the geomagnetic field measured in the U.K.; they were formed by a combination of Greenwich, Abinger and Hartland annual means obtained by comparing the annual means from three suc- cessive relocations of a British observatory (the most recent at around 51°N, 355°E), properly connected to coincide at their short overlapping running periods. Hereafter we will refer to this combination of time series simply as Hartland observatory data, i.e. with the name of the most recent observatory. Hartland observatory is characterised by a long first-rate tradition in monitoring the geomagnetic field. Therefore, the selection of this observatory was made on the basis of the quality of data; in particular we have taken into account of the apparent low degree of noise (statistically determined) and of the significant length of time series if compared with most observatories. Figure 5 shows the first- differences of Y component annual means, it is quite evident the potential presence of larger experimental errors for years preceding 1870; this is the reason why we have analysed only the last 128 years (1871.5 to 1999.5) of the combined time series. Considering that the aim of this work is to study eventual nonlinear properties of the Earth’s magnetic field of internal origin in terms of its secular variation, we have given particular relief to the Y component alone because, among the three geomagnetic components, it is the least dependent on the variations of the field of external origin. This fact can be explained expressing the Y component in terms of Spherical Harmonics (SH). It is a matter of fact that the magnetic effect of the solar activity perturbation can be represented as a combination of zonal harmonics and they are not contained in the SH expression of Y. Before starting with the bicoherence estimation, we checked the value of time series skewness verifying that is effectively non-zero. Equation (2.2) contains averaging over many statistical realisations of the same process. Since we have only one record of data, we replaced this average by the average on shorter segments of the time series (i.e. using the direct method). Therefore, after preliminary operations to obtain reliable data to be analysed with bispectral analysis we performed bicoherence for SV time series. Bicoherences of Y component shows very interesting results; a sharp peak corresponding to high values of bicoherence are clearly 1825 1850 1875 1900 1925 1950 1975 2000 -4 -3 -2 -1 0 1 2 Hartland SV Y [ n T /y e a r] year Fig. 5. Secular variation of the Y component of the geomagnetic field measured at Greenwich-Abinger- Hartland evaluated by first-differentiation. 286 Roberta Tozzi and Angelo De Santis -0.4 -0.2 0.0 0.2 0.4 -0.4 -0.2 0.0 0.2 0.4 (7,5) year 0.2625 -- 0.3000 0.2250 -- 0.2625 0.1875 -- 0.2250 0.1500 -- 0.1875 0.1125 -- 0.1500 0.0750 -- 0.1125 0.0375 -- 0.0750 0 -- 0.0375 f 1 (year -1 ) f 2 ( y e a r- 1 ) BICOHERENCE HARTLAND OBSERVATORY Fig. 6. Estimation of bicoherence for the Y component of geomagnetic field SV. It returns a quadratic coupling between 7 and 5 years. Fig. 7. The method of SPAQCE applied to the Y component of geomagnetic field SV from the combination of U.K. observatories reduced at Hartland location: three quadratic phase couplings are quite evident. QPC emerges from the superimposition of black squares and of white circles representing evidences of quadratic couplings in terms of differences and sums of couples of frequencies, respectively. The diagram is based on the 25% largest spectral amplitudes. 0.0 0.1 0.2 0.3 0.4 0.5 0.0 0.1 0.2 0.3 0.4 0.5 (2.4,25.6) year (5.5,8.5) years (5.1,5.8) years HARTLAND OBSERVATORY Spectral Amplitude Threshold = 25% sum difference f 2 ( y e a r- 1 ) f 1 (year -1 ) 287 Detecting geomagnetic field nonlinearities by bispectral analysis and a phase coupling nonlinear technique observed in fig. 6. According to the interpretation of bispectral analysis, this peak corresponds to a quadratic coupling between two periodic signals with periods of about 5 and 7 years. Results very similar to HOSA are obtained considering SPAQCE technique that returns some quadratic couplings (fig. 7), one of which is between two components with periods of 5.5 and 8.5 years. We consider the latter results more reliable on the basis of the estimated relative errors. 4. Conclusions The importance of this work lies in the consideration that bispectral analysis is a formal method that discloses out the presence of nonlinearities in a time series. It is frequently applied in physics with successful results (e.g., Hidalgo et al., 1995), but it has rarely been used to study geophysical phenomena (e.g., Clark and Bergin, 1997), especially those related to the geomagnetic field. The most important fact we deduce from our work is that bispectral analysis is a powerful method to detect nonlinearities in a sufficiently long time series. However, it may suffer, like other nonlinear techniques, from the shortness of time series such as the geomagnetic datasets. To overcome this problem we proposed a simple but more efficient nonlinear technique such as SPAQCE that gives results with a lower relative error. In fact, from the application of the method to synthetic data, quadratic couplings were always detected, thus confirming the nonlinearity of the investigated SV time series. In particular, the application of SPAQCE to annual means of a U.K. combined geomagnetic time series shows quadratic interactions between pairs of characteristic times, specifically 5-6, 5-8 and 2-26 years, providing important clues for the physical mechanisms generating the mag- netic field in the outer core. Certainly to reach more accurate results, more data and analyses are needed, and this will be the next step of this work, together with the search for the physical origin of the quadratic coupling in Y found here. At the moment, we can only make a suggestive hypothesis: for example it may be related to the variation of the length of the day, whose decadal variation is in turn probably correlated to electromagnetic and mechanic couplings between core and mantle. Another interesting hypothesis that would require more investigation is that of possible effects due to the differential rotation of the inner core. Acknowledgements We thank Drs. G.P. Gregori, V. Korepanov and Prof. A. Rivola for their valuable suggestions and comments, as well as Drs. G. Consolini and P. De Michelis for some discussion on the nonlinear characteristics of the geomagnetic field. REFERENCES BARRACLOUGH, D.R. and A. DE SANTIS (1997): Some possible evidence for a chaotic geomagnetic field from observational data, Phys. Earth Planet. Inter., 99, 207-220. CHANDRAN, V. and S. ELGAR (1994): A general procedure for the derivation of principal domains of Higher-Order Spectra, IEEE Trans. Sign. Proc., 42, 229-233. CLARK, R.R. and J.S. BERGIN (1997): Bispectral analysis of mesosphere winds, J. Atmos. Sol.-Terr. Phys., 59, 629-639. COLLIS, W.B., P.R. WHITE and J.K. HAMMOND (1998): Higher-Order Spectra: the bispectrum and trispectrum, Mech. Sys. Sign. Proc., 12, 375-394. DA L L E MO L L E, J.W. and M.J. HI N I C H (1989): The trispectrum, in Proceedings of the Workshop of Higher Order Spectral Analysis, 68-72. HIDALGO, C., R. BALBÍN, B. BRAÑAS, T. ESTRADA, I. GARCÍA-CORTÉS, M.A. PEDROSA, E. SÁNCHEZ and B. VAN MILLIGEN (1995): Nonlinear phenomena and plasma turbulence in fusion plasmas, Phys. Scr., 51, 624-626. HONGRE, L., P. SAILHAC, M. ALEXANDRESCU and J. DUBOIS (1999): Nonlinear and multifractal approaches of the geomagnetic field, Phys. Earth Planet. Inter., 110, 157-190. NAYFEH, A.H. and B. BALACHANDRAN (1995): Applied Nonlinear Dynamics (John Wiley & Sons Inc.), pp. 685. MACDONALD, G.J. (1989): Spectral analysis of time series generated by nonlinear processes, Rev. Geophys., 27, 449-469. MENDEL, J.M. (1991): Tutorial on higher-order statistics (spectra) in signal processing and system theory: theoretical results and some applications, Proc. IEEE, 79, 278-305. SUBBA RAO, T. and M.M. GABR (1984): An introduction to bispectral analysis and bilinear time series models, in Lecture Notes in Statistics, edited by D. BRILLINGER , S. FIENBERG, J. GANI, J. HARTIGAN and K. KRICKEBERG (Springer-Verlag), pp. 280.