Atmospheric Excitation of Polar Motion by John William Kuehne, B.M., B.S., M.A. Dissertation Presented to the Faculty of the Graduate School of the University of Texas at Austin in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy The University of Texas at Austin August 1996 Atmospheric Excitation of Polar Motion Dedicated to my Father, Robert A. Kuehne Preface Whatever Fortune brings, don't be afraid of doing things, (especially, of course, forKings) —A.A.Milne,NowWeAreSix This dissertation comes from many people in different ways. First is my supervisor, Clark Wilson, who made this research possible. Besides generous funding leading to a (probable) department record for student tenure and travel, Clark has allowed me to work unrestricted, and yet his knowledge, insight, and intuition into the polar motion problem has always guided this research. Stuart Johnson also figures prominently. Clark assigned him as an undergraduate assistant to this project, but unlike previous assistants he plowed through a pile of polar motion papers, and quickly left me feeling like the undergraduate. He lunged at the problem with verve and enthusiasm, making a number of important contributions to the 1993 The time that I spent working with Stuart was the most rewarding paper. and fun part of this study. My wife, Deborah Swart, while being a pillar in our personal relationship and growth, has also been important scientific co-supervisor. Much encouragement to an try something new or look at some data came from Debbie. father and I drove from This journey started during the Spring of 1984, when my Lexington, Kentucky to Austin to visit Hildegard and Rizer Everett (Hildegard and father not-so­ my are cousins). Their family has long been rooted in Texas, and he secretly wanted me to attend graduate school somewhere in Texas. My father did not live past 1984 to see me here at UT, but I think from our last conversations that he knew I was headed for Texas. Atmospheric Excitation of Polar Motion Publication No. John William Kuehne, Ph.D. TheUniversityofTexasatAustin, 1996 Supervisor: Clark R. Wilson Analysis of nonseasonal polar motion excitation and atmospheric mass equatorial angular momentum (EAM) over land for the period 1980-1993 reveals a clear pattern of high power and correlation during the northern hemisphere (NH) winter followed by low power and correlation during the NH summer. During the southern hemisphere winter there is significant correlation between the atmospheric EAM over midlatitude southern oceans and polar motion excitation indicating the existence of a dynamic atmosphere-ocean excitation. The atmospheric excitation power is too small to explain the large correlation during the NH winter. The momentum of wind probably accounts for the deficit in power. The implication of these results is that there are two main excitation sources each dominant at different seasons. Atmospheric mass redistribution over land forces polar motion during the NH winter, and a dynamic atmosphere-ocean response is important during the SH winter. The Chandler wobble frequency F and quality factor Q are estimated by least squares using polar motion data in combination with proxy excitations derived from Monte Carlo atmospheric data. tests show that this approach can provide useful traditional maximum likelihood which improvements over methods, use only observed polar motion. With less than a decade of good simultaneous polar motion and atmospheric data there is no new information about Q, but 90% confidence intervals for F are near ±0.004 cycles per year (cpy), making them comparable to maximum likelihood results based on nearly a century of historical polar motion data. Using a proxy excitation constructed from a linear combination of atmospheric EAM over land and SH ocean, the estimate of F is near 0.831 cpy (a period of 439.5 days), which is significantly different from the maximum likelihood value near 0.843 cpy (433 days). This difference can be understood in light of the assumptions underlying the two methods. Table of Contents Introduction 1 Atmospheric Excitation of Nonseasonal Polar Motion 6 Estimates ofthe Chandler Wobble Frequency and Q 25 46 Appendix A (Calculation of the Inverted Barometer Term) 48 Appendix B (Computation of Atmospheric Excitation Estimates) Appendix C (Neural Network System Identification) 60 88 Appendix D (Annual Component ofPolar Motion) Bibliography 91 Vita 94 Introduction One of the fundamental properties of the Earth which favors life is its rotation. distance Some other properties include atmosphere, oceans, tectonic activity, and from the Sun. The Earth must have attained its rapid rate of rotation early in its history. It is remarkable that all the planets except Venus and Mercury are rotating rapidly. The orientation of the Earth in a celestial reference frame is now routinely monitored with high accuracy using space geodetic techniques. These data are used to infer the orientation of the rotation axis in the celestial frame (precession and nutation), its spin rate (length-of-day), and its location in the terrestrial reference frame (polar motion). Gravitational torques of the moon, sun, and planets acting on the Earth's equatorial bulge slowly change the orientation of the angular momentum vector and the attending rotation axis. The resulting precessions and nutations, while ofconsiderable importance, are not the topic of this dissertation. The motions of the atmosphere, ocean, and liquid core change the rotation of the Earth. The conservation of angular momentum is the guiding principle for understanding how the rotation changes when one part of the Earth moves relative to solid Earth and another. For example, considering only the its atmosphere, the conservation principle allows one to determine the change in total atmospheric angular momentum by observing the change in the Earth's rotation (and vice versa) without knowledge ofhow the momentum is transferred. Polar motion was first observed nearly a century ago and meteorological as the excitation phenomena were immediately suggested source [Munk and MacDonald 1960]. Aside from the annual motion (see appendix D), the , meteorological hypothesis was unproved until Barnes, Hide, et al. [1983] showed that part of polar motion could be explained by fluctuations in atmospheric mass angular momentum. In conjunction with highly accurate space-geodetic polar motion data, their discovery was made possible by the advent of global circulation models, which produce estimates of pressure and winds at regular time intervals. Since their paper, a number of reports have appeared which suggest that between 9% and 25% of the polar motion variance is linearly related to atmospheric mass angular momentum, most notably Eubanks et al, [l9Bß]. To date there is no evidence that seismic events, or any solid-Earth mechanism, significantly excites polar motion. The explanation for the remaining polar motion variance is some combination of the following: Winds. The total atmospheric angular momentum is composed of two parts. The mass (or pressure) term is that part ofthe momentum which changes the inertia of the Earth, including its atmosphere. The mass term is estimated from surface pressure data, which is perhaps the best determined atmospheric parameter. The motion (or wind) term is the angular momentum of motion relative to the solid Earth. The wind term is difficult to estimate because it requires knowledge ofthe atmospheric velocity and density throughout the total volume of the atmosphere. Oceans. The atmospheric mass term is estimated by assuming that the oceans compensate local atmospheric pressure fluctuations, so that pressure fluctuations on the sea floor depend only on fluctuations of the mean atmospheric pressure over the This convenient entire ocean (the inverted barometer hypothesis). assumption probably is not true for periods less than about 5 days. The winds drive currents and change the pressure at the sea floor, but there are few current and pressure data available. Water storage variations. Variations is terrestrial water storage perturb the inertia of the Earth, but the details are unknown. Water storage variations cannot explain rapid polar motions, and may be important only at periods longer than a year (Kuehne and Wilson, 1991). Investigations in these categories have been hampered by lack of atmospheric, oceanic, and hydrologic data. Atmospheric modeling has been the most successful because the atmosphere is accessible and relevant to humans. Ocean models are which be driven by observations of surface winds, currents, and temperatures, may insufficient to model currents and variations at the sea floor. Surprisingly, pressure hydrologic variations are most difficult to model because there is no way to observe directly the amount of terrestrial water storage. The fundamental parameter - - precipitation is difficult to measure partly because it is not spatially correlated with itself. Hence good estimates of precipitation require a dense network of rain gauges. Evapotranspiration and runoff are also difficult to measure, and the required temporal integration of (precipitation-runoff-evapotranspiration) compounds the errors. The oceans are the most promising target to investigate for two reasons. First, they cover most of the planet, so that even slight deviations from the inverted barometer hypothesis might produce significant polar motion. Second, currents are variations at the sea partly driven by the wind and pressure surface, and global circulation models provide estimates of these data. Given that the oceans are partly driven by atmospheric processes which are observed and modeled, it is not entirely necessary to understand the physical processes occurring in the ocean to prove the hypothesis that the oceans cause polar motion. One can construct an empirical transfer function that transforms atmospheric data into polar motion for some time interval, and then test the success of the transfer function for a different time interval. This is system identification. The transfer function which I proposed to use is a kind of generalized non-linear function often called a neural network. Work on the neural network system identification showed that atmospheric angular momentum is repeatedly correlated with polar motion. However, this can be demonstrated directly by computing the correlation over short (150-day) overlapping windows of the time series. Such analysis of polar motion excitation and atmospheric mass equatorial angular momentum reveals a clear pattern of high power and correlation during the northern hemisphere (NH) winter followed by low power and correlation during the NH summer. The windowed correlation analysis also supports the hypothesis that the oceans are important in polar motion excitation. During the southern hemisphere winter there is significant correlation between the atmospheric mass angular momentum over midlatitude southern oceans and polar motion excitation indicating the existence of a dynamic atmosphere-ocean excitation, or some other effect of the winds. These fundamental and salient properties of the atmospheric and polar motion time series were not reported in the literature when this research began. With these discoveries, work on system identification ceased in favor of direct investigation into the excitation ofpolar motion. Appendix C contains a summary of the neural network system identification method used, and some examples. This dissertation is organized in two major sections and three appendices. The first section is a published paper in the Journal of Geophysical Research, and appears as Atmospheric excitation of non-seasonal polar motion by John Kuehne, Stuart Johnson, and Clark R. Wilson, v.98, pp. 19,973-19,978, 1993. This contains the fundamental results and interpretations of the research. The second section, which is published in the Journal of Geophysical Research as Estimation of the Chandler wobble frequency and Q by John Kuehne, Clark R. Wilson, and Stuart Johnson, v.lOl, pp. 13573-13,579, 1996, is an application of the results. The American Geophysical Union has granted permission for the use of this material for inclusion in this dissertation, including microfilm editions thereof. Appendix A, which appears in the first published paper, describes the calculation of the inverted barometer ocean Appendix B describes some new response. a atmospheric data available since the publication of the paper. Appendix C is complete derivation of the neural network filter, with some synthetic examples. Appendix D presents the results for the annual, i.e. seasonal, component of polar motion. Theoretical treatment of the physics of Earth rotation can be found in Munk and MacDonald [1960] and Eubanks [1993]. Atmospheric Excitation of Nonseasonal Polar Motion Abstract: Analysis of nonseasonal polar motion excitation and atmospheric mass 1980-1989 reveals a equatorial angular momentum (EAM) over land for the period and correlation clear pattern of high power during the northern hemisphere (NH) winter followed by low power and correlation during the NH summer. A special case ofthis pattern occurs for longer than 14 months (from January 1987 to March 1988) when the correlation throughout the NH summer remains statistically significant. During this epoch an average of 72% of the nonseasonal polar motion excitation power at frequencies between -30 and +l2 cycles/yr is linearly related to atmospheric EAM over land. During the southern hemisphere winter there is significant correlation between the atmospheric EAM over midlatitude southern oceans and polar motion excitation indicating the existence ofa dynamic atmosphere-ocean excitation. The atmospheric excitation power is too small to explain the large correlation during The the NH winter. The effects of winds probably account for the deficit in power. is that there are two main excitation sources each implication of these results dominant at different seasons. Atmospheric mass redistribution over land forces polar motion during the NH winter, and a dynamic atmosphere-ocean response is important during the SH winter. 1. Introduction A century after the discovery of the annual and 14-month components of polar motion by Chandler in 1892, the geophysical processes which excite polar motion at all frequencies remain enigmatic. Movement of air and water must be the excitation source of annual polar motion, although the details remain uncertain [Chao and Au, 1991]. At nonseasonal frequencies the variance of observed polar motion excitation exceeds the variance of atmospheric excitation, but significant correlation is observed [Eubanks et al., 1988]. We analyze polar motion data and accompanying estimates of polar motion excitation derived from meteorological data for the period 1980-1989. This is the decade in which two remarkable observational accomplishments became routine: of 5 space-geodetic determination of polar motion at periods days and less, with accuracy eventually exceeding 1 milliarc second (mas); and global determinations of on a atmospheric variations, providing values of various atmospheric parameters global grid every few hours. Our analysis of these data differs from that of previous investigations [e.g., Eubanks et al., 1988; Chao and Au, 1991; Preisig, 1992] mainly inthe useofsliding-windowestimatesofpower, correlation,andphase. Weillustrate anumberofimportantfeaturesaboutatmosphericexcitationofpolar motion, including the general nonstationary nature of the time series, a seasonal variation in power and correlation, evidence for a nonequilibrium ocean excitation, and a 14-month period during which 72% of the polar motion variance at frequencies from -30 to +l2 cycles per year (cpy) is related to atmospheric mass redistribution over land. These results demonstrate that during the study period, polar motion is principally excited by the motion ofair and water. 2. Excitation Functions Inferred Excitation from Polar Motion Data Conservation of angular momentum relates two quantities evaluated in the rotating geographic (or Earth-fixed) coordinate system: deviations in the angular momentum of the Earth (including atmosphere) and pole position, described by the complex functions m(t), respectively. The geographic coordinate system is defined by mutually orthogonal basis vectors ei, e2, and e3. ei and e 2 are in the equatorial plane with ei intersecting the Greenwich meridian and e 2 intersecting the meridian 90° east. e 3 intersects the geographic north pole. Quantities with subscripts 1,2, and 3 are components along these axes. The linearized Liouville equation is [Gross, 1992] i . dm(t) ,=+ 1 X(t) ;m(t) /?(A,o,r)e'sin is latitude, A is longitude (such that A increases eastward and is zero along the Greenwich meridian), g = 9.8 m/s 2 is the acceleration of gravity at the surface, 1.00 accounts for the effect of atmospheric loading, and the factor of 20.63 x 107 converts radians to mas. Equation (5) is typically evaluated as sums over the entire globe. However, in this study we obtain separate x’,’ f°r various geographic regions, especiallyfortheland oceanx™andsubsectionsofeach. The isostatic behavior of the combined ocean and atmosphere, sometimes called because we are interested the inverted barometer (IB) hypothesis, is not added to x^L in the separate contributions of the atmosphere over land and ocean. The IB term, computed from atmospheric data over land (see the appendix), is small compared to and its addition to X?L x?' does not significantly change our results. Significant , correlation between X? and x!'° during the southern hemisphere winter is evidence that a non-IB response is important, as explained in section 4. Because we do not have the gridded wind data to evaluate h(t), we use the and meridional winds ECMWF global x7'» computed from the zonal at pressure levels up to 10 mbar, again using the hydrostatic relationship to integrate over pressure.Dailyevaluationsofx,w/from 1986totheendof1989areusedinthisstudy. 3. Time Series Analysis estimates of the mean, trend, annual sinusoid, and The linear least squares semiannual sinusoid are subtracted from each time series to produce nonseasonal time series. The time series are resampled to a uniform 5-day increment after applying a zero-phase low-pass time domain filter with a cutoff at 0.8 of the Nyquist. The 5-day increment is chosen because there are no reliable variations present in the polar motion data at periods less than 3-5 days. Our method of time domain based analysis is upon sliding-window power, correlation, and phase estimates. This technique permits the analysis of time series whose statistical properties change through time. A window is a short segment of the time series which is prepared with zero mean and trend. Estimates ofthe variance a] , the squared magnitude of the complex correlation coefficient p], and phase 0 are T then computed for two time series %<5, as t» U(T) U(T) 'Lxx %x£, fUM n = e= r "'"j =1 arg %xXj 'Lx.x, t=L(t) J t=L(z) where the asterisk denotes complex conjugation, N is the number of 5-day data in the window,L(t)=T (Af/2) t/(r)=T+(iV/2)-1,andristhemiddletimeofthe , window,p]e[o,l] istheestimatedfractionoflinearlyrelatedvariancebetweentwo timeseries.If%t whererjisanoisetimeseriesuncorrelatedwith£then - , tt o are pliX’Z) is approximated by 2(£)Io](x) when p*(#,£), crj(n ) fall, or the 99% confidence level. This method is an improvement over the T standard approach which assumes that both time series are Gaussian, because the polar motion excitation time series is not statistically stationary. For 7V=3O (a 150-day window) the probability that any one p] will exceed 0.18 in the absence of true correlation is less than 1%. Frequency domain estimates of correlation and phase are based on the multiple taper method [ Thompson, 1982], which provides superior control of spectral bias and leakage. High-resolution spectral estimates X are computed by weighting the data fk w X, with orthogonal tapers, each with associated eigenvalue and computing the , k discrete Fourier transform of each resulting time series. The spectral estimates are then 2Wk XWk Xjk&jk =-n ——y~* "T Pf<-x,4)=-n ——\ w W W== t ‘fk ‘tk 'Yj '£i k \k=\ A*=i ) V*=i A*=i J r q o= x (X'o arg 'Lw s'fk ftfk J7 \k=l The number of tapers determines the trade-off between spectral resolution and q statistical confidence. For the spectral estimates with eight tapers presented in section 4theresolutionateachfrequency/incyclesperyear(cpy) is[/-3.36,/+3.36]cpy.The 99% confidence levels for p) are computed in the same manner as for the time domain estimates. In the absence of true correlation the probability of p) exceeding 0.5 at any one frequency is less than 1%. The interpretation of coherence is similar to the time domain correlation. The phase indicates temporal lag or lead between the time series as a function of frequency so that positive prograde phase means that x, lags • Spectral analysis of the complex time series produces a separation of prograde and ±,cot retrograde ( e ) components. This is especially appropriate for study of the Chandler wobble, which is a prograde normal mode. At other frequencies the separation may describe zonal prograde and retrograde motion that we observe in the atmospheric data. 4. Results WeevaluatetwodifferentversionsofxfL Thefirstisanintegrationoverallland. - The second version omits data from the Antarctic. Because this region is located x, at latitudes where the excitation is insensitive to mass redistribution (equation (5)) and the total area is small, the omission should have little effect. However, x? produces consistently higher correlation with polar motion excitation while retaining the variance of L suggesting that the omitted data are noisy. Trenberth and Olson , [l9BB] compared estimates of atmospheric quantities from the ECMWF and the National Meteorological Center (NMC) and reported discrepancies around Antarctica. This may indicate that the gridded-values over Antarctica are unreliable, presumably because there are insufficient observations from this region. Time Domain Estimates Examination of x? and ?L using the statistics a] and p] reveals a yearly Xt and correlation during the northern hemisphere winter recurrence of high power followed by low power and correlation during the summer (Figures 1 and 2). This result is insensitive to the window length N provided that it is shorter than 1 year. We use /V=3o (150 days). During NH winters, p](x"\xPL ) * s always statistically significant, greatly exceeding the 99% confidence level of 0.18, while during the summer it is insignificant except during 1987. Fig. 1. The variance of the polar motion excitation (solid) and mass atmospheric angular 2 PL momentum n >X ) - ForseveralNHwintersp]{Xm L)=0.7(oracorrelationcoefficientof0.84), implying that about 70% of a]{xPL ) is linearly related to o]{%m ) at these times. The phase is about zero when the correlation is significant except for isolated patches 2 2L early in the time series (Figure 3). However, o (x"') is larger than 12­ the total power of the time series remains unchanged, whereas decimation of the hour x^L reduces the total power by about 30%. However, aliasing of the observed the polar motion m and are different because of the resonant amplification near t Chandler frequency. Using x!’L as a model for x',\ we estimate that less than 5% of x the power of %"is due to aliasing. Analysis of high-frequency polar motion data may provide a better estimate of the aliased power. The wind contribution may also explain the power discrepancy. For some mPL 2PL wintertime pJ(^,^ ), g (X ) has ample power to explain polar motion t while retaining most of the original correlation (Figure 4). The importance of the wind contribution is also reported by Gross and Lindqwister [1992], Fig. 4. Complex correlation w squared g]{xpl + % ) (solid) and the variance ratio (dashed). Fig. 5. Complex correlation squared p2CT,*™') (solid) t 2 mPL and p (x,X ) (dashed) for r reference. Note that the ocean correlation peaks usually occur during the SH winter, 6 months out of phase with the land correlation peaks. Usually during the SH winter, p](,Xm P>X °) ls statistically significant. Virtually all of the correlation can be traced to the oceanic region roughly between 30° and 60° southern latitude (Figure 5). Between these latitudes, land covers less than 8% of the 2 P° total area. Denoting this subregion with a prime, o (x) also has a strong annual r When the correlation is pattern of high power during the SH winter (Figure 6). • significant, the phase averages so that x™ is about 20° east of X?We produce similar results with quantities computed from ocean surface wind data. Fig. 6. The variance (J2(XP 0 )» computed from atmospheric mass variations over the SFI oceans from 30° to 60° latitude. The variance is high during the SH winter. Foreach150-daywindowpast1984alinearcombinationofx^L Xt° and 1 >» of the form (Figure 7). The coefficients are computed to is significantly correlated with minimize the sum of the squared deviation from %? f°r each window centered at time T. The linear trend c t +d approximately accounts for variations longer than the T T s pL ?° xis the linear combination of % and X 150-day period of analysis. Thus % t t which best explains x? and maximizes the correlation during the 150-day period. Although the 99% confidence level increases from 0.18 to 0.23, > 0.23 for all times past 1984. Before 1984, p 2 is not significantly larger than T pl(Xm’X PL ) and we note generally that power, correlation, and high-frequency> spectral content for the polar motion and atmospheric excitation data all increase sometime after 1984. This indicates that the uncertainty of the data is relatively high before 1984. Fig. 7. Complex correlation squared p)(xm ,%s ). xl, is a t linear combination of atmospheric time series for land and ocean regions which minimizes the sum of the squared deviations from X? f°r each 150-day window centered on time 7. Becauseoftheyearlyvariationofcr2(x'n) ando]{xPL),thepresenceofnoisein 2 mPL either time series may produce a yearly variation in p (X )-We tested this >X r hypothesis by adding Gaussian white noise to %PL and computing p] between the original and corrupted time series. The addition of noise with an amplitude of several milliarc secondsreproducespartofthepatternofyearlyvariationinp2(x"',XPL)­ T However, such a noise amplitude in the atmospheric data implies large uncertainty in the planetary-scale atmospheric geopotential heights or, equivalently, spatially correlated noise. Because the noise characteristics of the atmospheric data are uncertain and the noise level of the polar motion data is well established from independent geodetic methods, we conclude that part of the yearly variation of p](Xm >X PL ) could be attributed to noise in the atmospheric data. However, some of theyearlyvariationisrealbecausep2(Xm>X?r) lsstatisticallysignificant(i.e., ator z above the 99% confidence level) past 1984. Frequency Domain Estimates Spectral estimates for various segments of the time series are all considerably different, supporting the time domain observations that the time series and their relationships are statistically nonstationary. The nonstationary behavior complicates time series analysis including the removal of trends and sinusoids and the estimation of power and coherence. The sliding-window procedure is a sensible way to handle the problem when the statistics of the time series change slowly with respect to the window length. The 1987 is a 14-month interval following January special case of the yearly patternofpowerandcorrelationwhen p\{%m^X?L)duringthesummeris higherthan average. The persistence of statistically significant windowed correlation suggests that the data are sufficiently stationary for spectral analysis of coherence and phase over the entire interval. This is the only period in the time series, greater than a few months duration, that is arguably suitable for spectral analysis. The graph of x? and PL’ X time series is unequivocal evidence that they are correlated (Figure 8). The t correlation coefficient between Re(^fL ) and is 0.67, and for the imaginary part it is 0.83. The magnitude of the complex coefficient is 0.78. The coherence spectrum suggests that about 60-80% of the polar motion and atmospheric excitation -30 to +l2 power is linearly related in the range cpy including the Chandler frequency of 0.843 cpy (Figures 9 and 10). At prograde frequencies between +l2 and +24 cpy the time series are essentially incoherent. The phase spectrum is close to zero and generally varies little at frequencies where the coherence is significant, which further supports the statistical significance ofthe coherence estimates. Fig. Ba. Rq(x!1L ) (dashed) and Re( x"’) (solid). The correlation / nr coefficient is 0.67, but X, has about 1/3 the variance of Spectral analysis of this portion of the time series is shown in Figures 9 and 10. Fig. Bb. \m{x',L ) (dashed) and Im(%™) (solid). The correlation coefficient is 0.84, but x^L has about 1/3 the variance of x?• Spectral analysis of this portion of the time series is shown in Figures 9 and 10. Fig 9. The squared coherence ofx!’L and x'," f°r January 1987 to March 19 1988. The 99% confidence level is 0.5. Fig. 10. The phase spectrum of X PL relative to x 7 • A positive t prograde (or negative retrograde) phase corresponds toa to lag of with respect XT'­ 5. Conclusions the and Our main result is discovery of a yearly recurrence of high power correlation during the NH winter followed by low power and correlation during the NH for the polar motion excitation and atmospheric excitation over land. The summer similarity between the interannual fluctuations is also striking. The variations power in the correlation and variance ratio (Figure 2) imply that there exist other excitations whose correlation with atmospheric excitation over land changes from year to year. The winds are a likely candidate for such additional excitation. Our second main result is the detection of yearly correlation between polar motion excitation and the atmospheric equatorial angular momentum over the midlatitude SH oceans during the SH winter. The correlation, although small, is statistically significant, hinting that a nonequilibrium response, winds, and wind-driven ocean currents mainly in the southern hemisphere contribute to the excitation during the NH summer. The average longitudinal difference of -20° may be related to the dynamic ocean response to pressure variations and winds. Determination of the oceanic contribution to polar motion will require a model of the dynamic ocean response more sophisticated than the inverted barometer model. The implication ofthese results is that there are two main excitation sources each dominant at different seasons. Atmospheric mass redistribution over land forces polar motion during the NH winter, and a dynamic atmosphere-ocean response is important during the SH winter. The persistent correlation between polar motion excitation and a linear combination of land and ocean excitations (Figure 7) further supports this hypothesis. Finally, we observe a 14-month period beginning in January 1987 when 72% of land at the polar motion variance is linearly related to atmospheric excitation over frequencies between -30 to +l2 cpy. The general shape of the coherence spectrum, remarkably flat with an abrupt drop at prograde frequencies between +l2 and +24 cpy, suggests that the spectral separation of prograde and retrograde frequencies is important and that the incoherent part of the spectrum is a genuine feature. At the Chandler frequency the fraction of linearly related power is about 0.7 (a correlation coefficient of 0.84). This is an important result because it unequivocally exceeds 99% confidence at the Chandler frequency and adjacent frequencies in a wide band. Thus it is good evidence that the Chandler wobble is meteorologically excited. Estimates of the Chandler Wobble Frequency and Q Abstract: The Chandler wobble frequency F and quality factor Q are estimated by least squares using polar motion data in combination with proxy excitations derived from atmospheric data. Monte Carlo tests show that this approach can provide useful improvements over traditional maximum likelihood methods, which use only observed polar motion. With less than a decade of good simultaneous polar motion and atmospheric data there is no new information about Q, but 90% confidence intervals for F are near ±0.004 cycles per year (cpy), making them comparable to maximum likelihood results based on nearly a century of historical polar motion data. Using our preferred proxy excitation we obtain an estimate of F near 0.831 cpy (a period of 439.5 days), which is significantly different from the maximum likelihood value near 0.843 cpy (433 days). This difference can be understood in light of differences in assumptions underlying the two methods. 1. Introduction This paper is concerned with the problem of estimating the physical parameters which describe the Earth's resonant Chandler wobble (CW). In a terrestrial reference frame, deviations in the equatorial angular momentum of the Earth %(t) X\(t) +i£2(o andpole position m(t) (t) are approximatelyrelated = = +im 2 by i dm(t) ,N = X(t) + m(t) G at G 2kF{\ + i / 2Q) = as a consequence of the conservation of angular momentum [Gross, 1992]. The 1 refers subscript to position along 0“ longitude, and subscript 2 refers to position along 90 East longitude, with the complex plane origin near the figure axis. One of the CW parameters is the resonant frequency F, which has approximate value of an or a 0.843 cycles per year (cpy) period near 433 days (d). F is proportional to the frequency of the Earth's rotation, and is a measure of the non-spherical nature of the Earth and its deformation in response to motion ofthe pole. The other CW parameter is the dimensionless quality factor Q whose reciprocal is proportional to the dissipation of energy near F. F and Q are whole-Earth properties which measure elastic and anelastic behavior at the lowest frequency normal mode observed for the Earth, and hence they are of general interest in geophysics. Eubanks [1993] provides a thorough reference for these and many other topics. Early values of F were obtained by harmonic analysis of m{t), but modern estimates of F and Q are based on a maximum likelihood (ML) method introduced by Jeffreys [1940, 1968] using the assumption that is a random Gaussian process for frequencies near F. After removal of seasonal and long period components from m(t), ML estimates of F and Q are those for which %(t) has minimum variance. Thus the ML method produces estimates of F and Q without explicit knowledge of X(t)­Barnes et al. [1983] and Eubanks et al. Beginning with papers by [l9Bß], evidence has accumulated that the excitation %(t) is related mainly to atmospheric variability. Polar motion is correlated with barometric pressure changes over a broad more than a range of periods extending from several weeks to year. The physical connection between polar motion and barometric pressure change over land is the 26 change in atmospheric angular momentum associated with atmospheric mass redistribution. Over the oceans the physical connection be different because the may tend as an oceans to respond to atmospheric surface loading inverted barometer, thereby annulling the effect at the sea floor. In a recent study, Kuehne et al. [1993] found that barometric over land are highly correlated with polar pressure changes motion during northern hemisphere (NH) winters. However, this correlation disappears during southern hemisphere (SH) winters, at which time pressure changes over SH oceans become correlated with polar motion. The physical connection between pressure changes over the oceans and polar motion remains uncertain given the inverted barometer assumption. Nevertheless, a linear combination of barometric pressure fluctuations from the two geographical regions land and SH oceans is correlated with observed polar motion throughout the year. Thus the linear combination may be useful as estimate for the full excitation. We examine a proxy whether the additional information in this and other proxy estimates ofthe excitation can improve estimates ofF and Q. Our approach is to solve a discrete version of (1) directly for F and Q via least squares (LS). We begin in Section 2 by describing the polar motion and atmospheric data, and develop proxy estimates for %(t) in section 3. In Section 4 we develop the LS estimator for the two CW parameters, with associated confidence intervals. With less than a decade of very high quality data available, our LS method has little new to sayabout Q, butoneestimateofFdifferssignificantlyfromMLvalues. 2. Data Polar Motion Data We use a subset of a polar motion time series called Space93, described in the annual report ofthe International Earth Rotation Service for 1993. The complex polar motion time series m is used to estimate x 7 with the two-point convolutional filter t [Wilson, 1985] • -inFT - = e 2 (f) dXd(f) , (5)(C-A)g JJ where 0 is latitude, X is longitude (increasing eastward from zero at the Greenwich meridian), and g = 9.8 m/s2 is the acceleration of gravity at the surface. Our method of estimating s reproduces almost identically the values provided separately by the ECMWF, but eliminates discontinuities associated with ECMWF model changes over 29 time. Equation 5 is evaluated by summations, usually over the entire globe. However, in this study we obtain separate integrals (summations) over land (x!’L ) and over SH (30° -60 S) oceans (X™)-This region of the oceans was found in Kuehne et al. to provide essentially all of the correlation between polar motion and atmospheric pressure variability over the oceans. 3. Calculation ofProxy Excitations We seek an estimate of the excitation x(t) containing all contributions (atmospheric mass and winds, oceanic mass and currents, ground water, etc.) to the terms in (3). In the absence of the necessary data, we construct proxy estimates of X(t) based on x?L and x!’° • The relationships among /L ;^o ,and %”1 are described in Kuehne et al. [1993] , (Figure 1). Briefly, is approximately equivalent to xf for the entire globe with an inverted barometer ocean response to pressure variations. The variance of x?L tends , to be larger during NH winters than during summers. x?L is well correlated with during NH winters, and essentially uncorrelated during summers. x!‘° has similar properties approximately 6 months out of phase, but correlation with x',” is generally weaker, and somewhat enigmatic because of the presumed inverted barometer behavior of the oceans. The x‘,’°-X? correlation arises mainly from three SH regions centered on longitudes 20°, 120°, 270°, all proximal to land areas. However, there is little correlation between X? and atmospheric pressure fluctuations over the adjacent land areas. The three regions do not contribute equally to the correlation at any given time. Rather, for any given SH winter correlation event, one region tends to dominate. 30 Figure 2 shows the locations of these three southern ocean regions which dominate the correlation with polar motion. Figure la. Correlation coefficients which form the .. - U.O basis for constructing the estimate The solid - . . | 0.6 \ L proxy 8 |\i JVi / IU 1 Ll line is the squared correlation I t between X? and °-4 " f]I\ j\ ) f\)[ a L 1 jl lI[V J 'll I|l V i computed for 150-day '\j 1 lit I !fill 1 Iill'll OP ft l'1' ' window centered on the date ry ft* tfi fWv % ffl\ \ffyfjVf " j , iJ i i P\„.f iAaA1f' r mwV i /V1 kv i 4 ’’ \}' y shown. The dashed line is the VV V, t[\j _ 1986 1987 1988 1991 1992 1993 correlation between an(i X? X™. The correlation for x?' is always high during the winter. Xt° follows the same rule 6 months out of phase, although the correlation is generally lower and less predictable. The 99% confidence level is 0.18. The 150-day variance for tends to be high during winter and low during the summer, and vice versa for xP,° • Figure lb. ls a Hnear combination and P° of %fL %t which is correlated with x? throughout the year. 31 —) Figure 2. Map of the number of times the squared correla­ -1 / : tion coefficients exceed the \ /­ 99% confidence level for the U period 1985-93. Correlation is y A between atmospheric angular J momentum (inertia term) over 140-47 9 41 82 "T j 30 6 111 23 70 / |/l67 0 - each region and polar motion excitation. The presumed iso­ 160= 7 static compensation of the 180 0 50 ioo 150 200 250 300 350 ocean and atmosphere should cancel the effect of atmo­ spheric pressure fluctuations of the oceans. Thus, a dynamic atmosphere-ocean re­sponse, or some other consequence of the winds, may be responsible for the correla­tion. Of all the oceans, only pressure fluctuations over the southern hemisphere region is significantly correlated with polar motion. Although the three most significant re­gions, centered on longitudes 20°, 120°, 270 are proximal to land, the atmospheric , angular momentum over the nearby land regions is less correlated with polar motion, which suggests that the effect responsible for the correlation occurs over the oceanic region. We construct three proxy estimates of x(0 from xfL and x™ The first is a • linear combination of xfL and X most like X? in a least squares sense, denoted by X*. The second is a portion of x?' which is well correlated with x? denoted by » where r is a subset of the times t in the full series. The third proxy excitation is the entire series x?L used in a weighted least squares sense, with weights equal to the correlation coefficient with X? determined over 150-day windows. Calculation of the Proxy Excitation X* The high NH wintertime variance and correlation of x?' tends to occur when the variance and correlation of x!’° ls low and vice versa. Hence a simple linear combination tends to be correlated with throughout the year (see Figure 1). Thus we construct a proxy excitation of the form - + ar+a x? a\x^L aixr°+ 3 A (6) where the coefficients are chosen to obtain a combination of the two atmospheric series which is most like the excitation determined from polar motion via (2). That is, we find the least solution to squares h 4°f 12=x: (7) X 1 x? if“'l h: i •• # • ' L*LJ for daily time series. The complex valued coefficients a, and a 2 serve to scale and a rotate (in longitude) the elements of their respective time series, while a 3 and A and account for a relative trend and offset, x? is initially estimated with F=0.843 cpy not critical because the series are all broadband. The <2=175, but this choice is solution of (7) is found from 8.6 years of data (3195 samples), and gives Because there a =0.96e89 i a 2 =0.36e~*ri , and both a 3 and aA nearly zero. are , relativelyfewparameters usedtofitthetwoatmosphericseriestox’t\theresulting proxy excitation should contain essentially the same number of degrees of freedom as either or X,°¦ 33 The proxy excitation obtained in this way has a broad band squared correlation and coefficient with x? of 0.51, which is significantly higher than the values 0.31 0.16 obtained separately for L and X™ ¦ln the frequency domain, coherence and phase spectra (Figure 3) are evidence that xt and X? are strongly correlated over virtually the entire bandwidth of the data. Retrograde frequencies are somewhat better correlated than prograde, although the high coherence at prograde frequencies between 0.02 and 0.025 cpd (40 to 50 d) is notable. Coherence remains significant across the band containing the Chandler frequency, and the stable phase spectrum further supports the coherence. The Figure 3. spectral coherence and phase, computed by multitaper a Fourier method, between x? and Xt The statistically • significant coherence is broadband, extending from -0.1 to 0.1 and at the cpy, Chandler frequency exceeds the 95% confidence level of 0.4. • PL Calculation ofProxy Excitations X T collection The least squares estimator developed in section 4 can be applied to a ofdisconnectedportionsof%?L Consequently,forthesecondproxyestimatex we . those times for which x*L is well correlated with at r select portions of the series m y We eliminate those data below a specified cutoff level of squared correlation . 34 coefficients, determined using a 150-day window centered on t as described in Kuehne et al. [1993]. We use squared correlation cutoff levels between 0 and 0.34 to form a set of 35 different estimates of x^• F°r example, when the cutoffis 0.018 (the 99% confidence level for the 150-day squared correlation coefficient in the null case) about 20% of PL discarded, mainly during NH summers. xt 4. Estimation of F and Q The Least Squares Estimator When independent estimates of x 7 an d are available, (2) can be written as a set of m t simultaneous linear equations mm C 0 t0 \Xt +TI2 -T/2 tO+TI20 mmC (^) tl+T/2 t-TI2 2_~ Xt x_ x inFr iaTie~ -ie~ e inFT where unknowns c. = and c. = are the complex coefficients in (2), 2 1 oT oT and xs the independent estimate of x? Our proposition is that if a good estimate i • t of Xt is useh on the right hand side of (8), then the solution of (8) ought to be superior to a ML estimate, which ignores the additional information contained in x, • We solve (8) by least squares using one of the proxy estimates of x, with additional parameters to account for a residual linear trend, a coordinate system offset, and sinusoidal deficiencies in the proxy at annual and semi-annual frequencies. As before, daily time series are used. The least squares solution is unweighted for the estimates x* and » and weighted by 150-day correlation coefficients when proxy 35 using the entire series x!’Las a estimate. LS estimates of the two CW proxy parameters are computed by dividing -c2 by ci to obtain e ,CTI , i.e., arg -Trtf 2kT log|-c2 /c, where we use the convention that lower case/and q denote estimates ofF and Q. Bias and Confidence Intervalsfor the LS Estimates Ensembles of synthetic data can be created using a random number generator in which the true values of the two CW parameters are known precisely, and the LS (or ML) estimator can be applied to the synthetic data to measure estimator performance in terms of bias and confidence interval width. This Monte Carlo approach is useful when knowledge of the excitation is approximate. The Monte Carlo results apply to noise-free data, and thus may underestimate the width of the confidence intervals for real data. However, they provide a valid way to compare the performance of LS and ML estimators. Before discussing the results ofthe Monte Carlo evaluation in detail, it is useful to examine the shape of the sum of squares (error-norm) surface as a function of trial estimates of Q and F. A point located at the minimum of this surface is a LS solution to (8). The surface is shown in Figure 4 for a noise-free synthetic data case in which zero no Q =175 and F=0.843 cpy. The error-norm is essentially with prominent minimum for a large range of trial values about the true Q. Thus accurate estimates of measure Q are difficult to obtain, even with perfect data. Because lIQ is a of the Chandler frequency spectral line width, a qualitative explanation for the problem in familiar one estimating Q is a in spectral analysis. A short time series (8.6 years is 36 only a few Chandler periods) does not contain the information needed to distinguish between a moderately narrow and an extremely narrow spectral peak at the Chandler frequency. In contrast, the shape of the error-norm surface in Figure 4 shows that F is easily determined. Figure 4. Typical error norm surface of (7) on a logio scale. Using (2), m is deconvolved t for f =0 843 and ‘ cpy .0001 t0X' m - •001 12=175. Then is repeatedly t deconvolved for a of " 01 range 450 variance of the difference with Q is computed. The surface " Frequency (cpy) °-9 bU thus graphically describes the behavior of the least squares solution of (7). For trial values of F near 0.843 the error is less than 0.001 for trial (2>loo, and on a linear scale the minimum of the surface would not be detected visually. This surface reveals one of the fundamental difficulties of estimating Q using the least squares procedure. If Q is large, then any larger value will do nearly as well in minimizing (7). In contrast F is well determined. These results are closely connected to the Monte Carlo confidence levels in Figure 4. To evaluate LS bias and confidence intervals for we first use (3) to transform the 8.6 year series m into % ™ with known F and Q. We then generate an ensemble of t to ensemble synthetic proxy excitations, x% by adding an of independent Gaussian noise time series, all having the variance needed to reduce the squared correlation coefficient between and % s to 0.51 (the same as between X? and the t actual proxy series X? )• From the ensemble of synthetic proxy series, we obtain an 37 are ensemble of estimates/and q via (8). In this process, both positive and negative q obtained because there is no positivity constraint in (8). Negative q, which are not physically meaningful, occur most often when Q is large and are simply discarded. A similar procedure is used to determine LS bias and confidence for each of the 35 versions ot x • r The LS section of table 1 summarizes Monte Carlo statistics applicable to estimates made with xtQ varies over the range 30-1000 with F held constant. The • precise value ofF used in these experiments is not critical. The LS estimate is labeled E[q].TheintervalfromL[q]toU[q]contains90%oftheensemble with5%ofthe q, valuesbelowL[q] and5%ofthevaluesaboveU[q].Ingeneral,thedistributionof q depends on the length of the time series and the squared correlation coefficient. When the squared correlation coefficient is zero, or the time series is short, q tends to be near zero. For the case considered (a squared correlation coefficient of 0.51), q approximately follows a beta distribution biased toward zero (Figure 5). For Q=3o the LS q is about 30, while for <2=loo the LS q is about 65. When Q>loo, LS q are hopelessly biased too low. Thus if Q is less than about 70, a reliable LS estimate is possibleusingx? dXf¦IfQIsmuchlarger,theLSestimatorcannotdistinguish an between the case of large and small Q. In contrast to these problems, the distribution 90% ofLSfisnormal,unbiased,andreasonablyindependentofbothFand Q,with a confidence interval of ±0.004 cpy for the case where we use %*. 38 Figure 5. An example of a 084 Monte Carlo experiment describedinthe textforQ-30. 835_ - & . The most likely q is about 28 § 083 ¦ with a sharp cutoff for smaller I" q. The probability of obtaining 0g25 a large q is small, but there is 082 no sharp cutoff. For Q> 100 „0.815 1 1 1 1 1 1 1 1 these effects are greatlyi ° J 10 20 30 40 50 60 70 80 90 100 110 . , _,, Q 1), and magnified (see Table we are left with double jeopardy: the extreme bias that we will most likely means underestimate Q, but neither can we assign a reasonable upper bound. Using Xt L 90% confidence for/is about ±0.006 cpy regardless of the correlation ’ cutoffchosen (Figure 6). Thus the confidence gained by omitting uncorrelated data is lost by using fewer data. Similarly, q derived from %PL are less reliable than for %,A . Although we do not compute them directly, confidence levels for the weighted LS estimate on must be similar to %pL x. For comparison with the LS estimator, we present Monte Carlo evaluations of Jeffreys' ML estimator described by Wilson and Vicente [1990]). The method developed by Jeffreys [1968] uses the sequence of Fourier coefficients at 6/7 cpy, computedfromnon-overlapping 14monthsegmentsofpolarmotiondata,asthebasis for the ML estimate. Precisely following the Monte Carlo procedure in Wilson and of Vicente, we examine the case of 9 years of daily data, and also the case of 86 years monthly data, corresponding to estimates made from the historical polar motion data available since the turn of the century. Again, the evaluations are made for the case of noise-free data. The results are summarized in the middle and right sections ofTable 1. ML q are biased toward zero for the 9-year case, and cannot be assigned useful 39 confidence intervals. We find from Table 1 that LS from 9 of data are similar q years in quality to 86 year ML estimates. ML/ are normally distributed and unbiased. The 90% confidence interval of ML/for the 9 year case (daily samples) is ±0.004 cpy , and ±0.005 cpy for the 86 year case which contains fewer total data. Wilson and Vicente obtained similar ML confidence intervals for/in their study, as shown in Table 2. Table 1 LS (8.6 years X?, daily) ML (9 years, da ily) ML (86 years, mon thly) 90%/: cpy±0.004 90%/ :±0.004 cpy 90%/:±0.005 cpy Q EM LM UM EM LM UM EM LM UM 30 25 25 40 19 16 22 35 20 60 50 45 30 120 22 18 24 50 30 100 70 50 40 230 23 19 25 60 35 140 150 75 55 925 23 20 26 70 40 215 500 85 70 1910 23 20 26 80 45 260 1000 105 75 2160 - - - 95 50 270 Table 1. Left section: approximate statistics (nearest multiple of 5) for the LS estimator using xt (8.6 years of synthetic data), based on 10,000 trials. Middle of section: approximate statistics (nearest digit) for the ML estimator for 9 years synthetic data and 1000 trials. Right section: approximate statistics (nearest multiple of 5) for the ML estimator for 86 years of synthetic data and 1000 trials. E[q] is the and the interval between L[q] and U[q\ contains 90% of the most likely value for q, trial values. The LS statistics are somewhat worse for which contains fewer data. 40 Table 2 Estimate q (90% range) /(cpy) LS(*,\ro,) 72 (30,500) 0.831 ±0.004 LS(^fL ,m,) (weighted LS) 77 (30,500) 0.836 ±0.006 ML(Wilson and Vicente [1990]) 179 (74,789) 0.843 ±0.004 ML(m,) - 0.840 ±0.004 Table 2. Confidence for x? is based on the assumption that contains the Chandler excitation. Confidence for is larger than ±0.004 cpy because it incorporates fewer data. Confidence for the ML estimates depends on the assumption that the excitation is a Gaussian random process at frequencies near F. Confidence levels for Q using x'rL and the weighted x 7 are not computed explicitly, and we state the confidence levels computed for X*•In both ML and LS cases, Qis poorly determined from 8.6 of data. years 5. Results excitations Table 2 gives the LS estimates for F and Q obtained using the proxy and % PL the ML results reported by Wilson and Vicente, and the ML results for , t m estimate x we find an LS value forFof 0.831 cpy with a 90% Using the proxy r a confidence interval of ±0.004 cpy, corresponding to Chandler period of about 439.5d. The raw estimate of Q is 52. From Table 1, the bias is about 20, and the corrected value for Q is about 72, with an approximate 90% confidence interval of 1 [30,500]. Confidence is determined from Table using the raw value of 52, which is above the 90% upper limit of <2=3o and below the 90% lower limit of £>=soo. The results for % PL using a range of squared correlation cutoffs between 0 and 0.35 are shown in Figure 6./reaches a maximum of .844 cpy near a cutoff of 0.18, 41 which is the 99% confidence level for squared correlation 150-day window. on a Nearly half of the data are omitted when the cutoff is 0.35, and as this value is approached the least squares solutions deteriorate, producing negatives values for q. Confidence levelsfor/areconstantacrosstherangeofcutoffs,implyingthereis no formal basis for choosing the cutoff. However, cutoffs for %PL in the 90%-99% range (0.10 to 0.18) provide a reasonable compromise between the desire to omit portions of the atmospheric series which are uncorrelated with polar motion, and the need to retain a sufficient quantity of data In this cutoff range /is between 0.837 and 0.844 . ±0.006 cpy cpy. Figure 6. / as a function of cutoff. L contains those Xt of data from portions xfL whose 150-day squared correlation with x? exceeds the cutoff value. At a cutoff of 0.35 approximately half the data are discarded. Confidence intervals for F increase in slightly with cutoff. Bias estimates of Q is severe for large cutoff. The weighted LS estimate of F is 0.836 cpy ±0.006 cpy using x with a bias-corrected estimate of 77 for Q. This estimate of F is not independent of estimates and are actually weighted LS estimates. The derived from because both x',’L weights for x?L are cither lor 0 (we either accept or reject a datum), whereas every 42 datum in ?L is used, but with correlation coefficient weights which range between 0 xt and 1. The ML estimate of F based on the same 8.6 years of polar motion data used in the LS estimate is 0.840 ±0.004 cpy. The ML Q estimate (19) is meaningless because, according to Table 2, a similar estimate would be obtained for virtually any physically reasonable Q. In the sense that their confidence intervals overlap, LS estimates of Q are consistent with the ML estimate made by Wilson and Vicente [1990] from the longer historical data (Table 2). The quality of the LS estimate, as measured by confidence are no new conclusions to be interval width, is not significantly better. Thus there made about Q from our study. The ML estimate of F derived from 8.6 years of data is consistent with the ML value from the much longer data set used by Wilson and Vicente [1990]. However, ourLS estimateofFderivedfromxtisnotconsistentwiththeMLvalueinthe sense that the confidence intervals do not intersect. To interpret this discordance, one must understand that estimates and confidence levels are conditioned on the validity of the underlying assumptions. The underlying ML assumption is that the excitation behaves as a Gaussian random process for frequencies near the Chandler frequency. The Gaussian assumption is a computational convenience to permit estimation of F and Q without explicit knowledge of the excitation, and it has no basis in observation. The excitation is linearly related to the true underlying LS assumption is that the proxy polar motion excitation, which is strongly supported by the significant broadband coherence and stable phase in Figure 3. Thus, discordance in the estimates may occur, given that the assumptions underlying the estimators are different. 43 A simple interpretation of the discordance evident in Table 2 is that the ML estimator of F is biased. To understand how this that the might occur, suppose excitation spectrum is not perfectly flat as the ML Gaussian assumption implies, but instead slopes upward between the Chandler frequency and the prograde annual frequency. This slope would cause the spectral peak of polar motion to be slightly asymmetric (skewed toward higher frequencies) with respect to its central value, the resonant Chandler frequency. The ML estimator, which uses a finite-width part of the spectrum to find the center of this Chandler frequency spectral peak, would then yield an estimate which is biased, that is, slightly too high. The LS estimator should have no comparable bias if the proxy for the excitation is a good one. We believe that xt I s a good proxy because of its broadband correlation with polar motion excitation (Figure 3). xt is preferred over either of and because it includes a contribution from oceanic regions, although this contribution is not fully understood. Thefore the Chandler frequency of 0.831 cpy (439.5 ±3 mean solar days) obtained It from x* is not only plausible but should be preferable to ML estimates. is a similar result of 439 interesting to note that Jeffreys [1940] obtained ±l2 mean solar days. We conclude that the LS approach can provide useful improvements over there traditional ML estimates. That is, is an important contribution from the additional information in the independent meteorological data used to determine the proxy excitation time series. The LS approach can be applied to recent space-geodetic polar motion data with about the same formal confidence levels for both F and Q as attained by the ML method with 86 years of historical polar motion data. Further with increased time series length, and with improvements in LS results will come 44 refined estimates of Xr Two such refinements may be possible with existing data. The first is to include atmospheric motion effects using wind data. The variable estimates among various meteorological analysis centers suggest that this will not be and easy, we have avoided using wind data for this reason. A second improvement would probably result from better understanding of the physical basis for correlation between and x?> which implies an improved understanding of the role of the oceans. Beyond these two effects, other data, especially pertaining to fluctuations in water storage among ice, ocean, and land reservoirs, will be required to improve estimates of • x, This time Figure 7. domain sample of X? (solid) and xt (dashed) is direct evidence that polar motion and atmo­spheric excitations are closely related. 45 Appendix A Calculation of the Inverted Barometer Term In his pioneering paper, Jeffreys [1916] developed a method for computing the polar motion excitation caused by atmospheric pressure fluctuations over the ocean. so that He assumed that the oceans are in isostatic equilibrium with the atmosphere each excess millibar of atmospheric pressure (relative to the mean pressure over the surface ~1 oceans) depresses the accompanying ocean cm. This idea is commonly called the inverted barometer (IB) hypothesis. In this appendix we essentially follow the prescription given by Jeffreys, which for the moment assumes that atmospheric and oceanic masses are conserved. An immediate consequence of the IB hypothesis is that fluctuations in ocean bottom pressure are everywhere the same and depend only on fluctuations in p' =(Mg-pLA )/A 0 1 L with which is the average atmospheric pressure over the oceans, pL the average the area of land and ocean, and M the total mass of the pressure over land, AL ,AQ atmosphere. Omitting the leading constants, the excitation integral is 46 sin0cos dS j*p(f)e,x Earth cos l*dS + dS =p' \sin $ (f)e Jpsinocose'A Ocean Land = dS -(Mg-p,A,)/A|sin0cos(f)e,x 0 Land + adS jpsin (f)cos (f)e Land of where the integral on the third line is a consequence a dS = 0 | sin (f)cos(f)e Earth Thus the IB effect can be computed from atmospheric data over land alone. One exception to conservation of atmospheric mass is the variable water vapor in the atmosphere. If the water vapor were concentrated entirely over land, the seasonal fluctuations would amount to little more than a millibar, and thus it is pressure probably a small correction. Also, the ocean mass is not constant. As Jeffreys a suggests, these problems belong to comprehensive study of precipitation, runoff, and evapotranspiration. 47 Appendix B Computation ofatmospheric excitation estimates Since this research began, many new meteorological data have become available. The old European Centre for Medium-Range Weather Forecasts (ECMWF) initialized geopotential data have been discontinued and supplanted by the new ECMWF are uninitialized data. The new ECMWF surface analysis provides direct estimates of surface pressure (and many other surface parameters) without requiring interpolation from geopotential heights, as described below. The NASA data provide surface pressure estimates 8 times per day on a finer resolution grid. Here is a comparison of angular momentum time series derived from these data which reveals important differences. Table 1 Agency Rate Epoch Resolution Data Symbol PL x 1980-89 2.5x2.5 ECMWF 2/day Initialized geopotential PLG x ECMWF 2/day 1985-93 2.5x2.5 Uninitialized geopotential PLS ECMWF 1985-93 2.5x2.5 Surface x 2/day pressure PLN NASA 8/day 1985-87 2x2.5 Surface pressure x are evaluated over land The atmospheric angular momentum time series in table 1 regions only. This is nearly the same as the integral over the entire globe using an inverted barometer ocean response, and corresponds to polar motion excitation better than the total integral. 48 and x!>LS Comparison ofx^LG The time series derived from these data do not have the same mean value (Figure 1). Certain systematic differences in the topography data can produce such an offset, such as a slight shift in the land-sea mask that may arise from the ambiguity of defining a coastline at the coarse resolutions used, without strongly affecting the zero-mean part of the time series. For this research the mean value is unimportant, and is subtracted before computing estimates of correlation and phase. Figure 1. Imaginary (90° east) components of the atmospheric mass momentum angular computed from the ECMWF surface data (xPLS), the ECMWF 1000 and 850 mbar geopotential heights (XPLG), and the NASA surface data (% PLN )i f°r 100 days. Far more serious are the large jumps and offsets in both real and imaginary °f which one is clearly visible in Figure 1. Several other large components in % offsets are evident when compared with x*LG (Figure 2), or %PLN These offsets are . clearly artifacts of the ECMWF analysis, and are not properties of the atmosphere. Also visible are smaller annual variations in geographic phase and magnitude of the difference. Aside from the offsets, x PL° and x!’LS are nearly identical. The coherence and phase estimates show that the time series are linear versions to within a scale factor 49 for periods between 4 and 2 days (Figure 3). At shorter periods the coherence years and phase remain highly significant indicators of linear relationship. The power ratio (Figure 4) is essentially unity except at periods less than 2 days. \% PLS PLG —— — . 3r—T T T T T T Figure2.-% The steps, which are caused by 2 5% sSo**** PLS tend wreck ' X to spectral ’ even estimates, though they . 2 might not be obvious in plots w c ' 1 1 -5 'ofthe whole time series. There cc is a small annual variation | implying imperfect agreement 0.5-in the annual sinusoid, which is a large component of both °° 1000 2000 3000 4o°o 5000 6000 time series. The agreement 1/2 daysfrom January 1, 1985 & also fluctuates annually, and there are several other features which reflect changes in the estimation procedures used at the ECMWF. The coherence and Figure 3. phase spectra for %PLS and PLG demonstrate that the time X series are essentially linear versions. These estimates are based on a 3000-day segment when PLS has no obvious % steps, either alone or in relation to the 3 other time series considered. 50 Figure 4. Power spectrum ?LS PL° ratio of / % The ratio • x to 0.5 is essentially unity up cpd. Because they are nearly identical except for the errors in the surface pressure data, XPLG (from the interpolated 1000 and 850 mbar geopotential heights) is more reliable than x‘,ILS (from the surface analysis). Even where an offset can be directly detected in XPLS , there is no way to remove it completely by detrending on either side of the step. Because of the noise-like statistics of the time series, a plot of x PLS may not reveal , that anything is askew, but the effect of the steps is fairly devastating to coherence and phase estimates at both low and high frequencies. Comparison of and x!‘LG Substitution of the new /c for does not significantly change the results or conclusions we make about polar motion, and comparison of these time series shows whythisisso.Themagnitudeof (FiEure5)issmall,asidefromanoffset and a small slope during the first 250 days. At periods from 4 days to 5 years the time are highly coherent. At periods shorter than 2 days the coherence drops below series never becomes the 99% confidence level of 0.5, especially the prograde part, but it 51 ’ totally insignificant (Figure 6). Xt L has about 6% more power than x‘,’L((Figure 7), with a fairly constant factor across the coherent band. Thus, for the polar power motion comparisons where all the time series are decimated to a 5-day sample interval, these two time series are practically identical except for a small constant factor. power ]y pl • Figure 5. -X™ Figure 6. Coherence and PL phase spectra for % and PLG The times series are es­ . X versions be­ sentially linear and 4 tween periods of 8 years days. The coherence is proba­ bly statistically significant at shorter periods. Note the sharp increase in coherence and stability in phase near 0.83 cpd, explained in the text and Figure 11. 52 7. Power Figure spectrum ?l LG ratio of X dx' • ' an Comparisonofx','LN dX? ( The most important difference between the NASA and ECMWF data is the the sampling rate of 8/day and 2/day respectively. In order to compute xl^~X',’LG > ECMWF data are interpolated to 8/day by zero-padding in the discrete Fourier transform. This forces the spectrum estimates of the interpolated time series to match that of the original data. The inverse transform recovers the original data, and the interpolated points are optimal in the spectral sense. As with the other estimates, mean values are different (Figure 1). The high frequency content of x!’lN contributes to the difference, and a periodic change in geographic phase is evident (Figure 8). The power spectrum of x!’LN (Figure 11) shows strong 1,2, and 3 cycle per day (cpd) lines. Visual inspection shows that %rLN has a strong 1 and 2 cpd variation. The t 3 cpd line is a harmonic of a non-sinusoidal 1 cpd variation. Between periods of 2 has and 4 days, x™ and X™ are essentially linear versions (Figure 9). %, LN years 53 about 15% more than x!‘LG hut in the coherent range of -.25 to .25 cpd the power > discrepancy is about 8% (Figure 10). Figures. I*™-*"*. Figure 9. Coherence and phase spectra for x?LN and PLG The time series are . X essentially linear versions between periods of 3 years and 4 days. Note the sharp increase in coherence and stability of phase at -0.83 cpd, as explained in the text and Figure 11. To bring this full circle, Figures 12-14 compare x™ with XfL (the old initialized ECMWF data). Except for the power ratio, the results are nearly identical to those for PLG PLN has about 8% more power than L but in the coherent of -.25 to X , range • %t ratio is about 1. 25 cpd, the power 54 10. Power Figure spectrum PLN PLG ratio of % I% . io'10 i r— i 1 1 , , , Figure 11. Power spectra for -, JO’2 I : xPLN (top) and a:™0 (bottom). - 14 both or duration of year |io‘/\ I I - ' io 16 the NASA data available to us. 10' 1"-The NASA data are sampled 8 -jj 3££ ; • ,, Cycles per day i .i times per day, and the ECMWF twice per day. The r -12_ I ECMWF 2/day 10 s ¦ A _ ECMWF data are interpolated £ ¦ jJ \r to 8 times per day to allow a ___—­ ¦ j coherence and phase estimate. 10 '' "*32 1 1234 Inthe -1to 1 the range cpd cydespe,day spectra are quite similar. The be aliased in the ECMWF data. The 1 and 2 NASA spectrum shows features that may cpd line spectra are plainly visible in the time series, and the 3 cpd line spectrum is a harmonic of a nonsinusoidal 1 cpd variation. Of special interest is the increase in continuous power near -0.83 cpd in both spectral estimates, with an attending increase in coherence and phase stability (Figure 9). This is one of the retrograde normal modes of the atmosphere which is present in the atmospheric angular momentum because of its harmonic purity. 55 W PU* Figure 12. -fL Coherence and Figure 13. phase spectra for %PLN and ?L for the NASA data X, available to us (1985-87). The time series are linear versions between periods of 3 and years 4 days. Note the increase in coherence of and stability phase near -0.83 cpd. 14. Ihe rigure power spectrumratioof%FLN!xPL . Summary ofComparisons The differences between the time series listed in table 1 fall into four classes: 1) Errors 2) Mean value 3) Broadband power 4) Spectral content at periods less than 4 days The most obvious errors are the steps in xfLS Some are obvious by direct . inspection, and others are revealed in comparison with the other time series. Aside from the steps, this time series is practically identical to %PLG Because the steps t. destroy spectral coherence at the lowest and highest frequencies, this time series and their gridded data should be avoided. All the time series have different mean values. Different treatment of topography data reduction could produce some of the difference, but I can only speculate about the cause of the discrepancies. The mean value itself is unimportant for comparison 57 with polar motion, as the first step in analyzing these data is removal of the mean and trend. PLN PLG Xt has 15% more broadband power than %( , and 8% more than x?L In the • coherent band of -.25 to .25 cpd, the discrepancy is reduced about 8%, so that PLN x matches in this range, while still retaining 8% more power than x' 'LG In this • t regard the old ECMWF data are closer to the NASA data the than the new ECMWF data. The spectral content ofall the time series differ at prograde and retrograde periods less than 4 days. The NASA data are sampled at 4 times the rate of the ECMWF data, and show spectral features that may be aliased, or incorrectly estimated in some other way, in the ECMWF time series. One feature that is common to all the power spectra is an increase in continuous power at -0.83 cpd, or a retrograde period of 1.2 days, caused by a barotropic normal mode of the atmosphere (Alquist, 1982). This is evident in all the coherence spectra, which increase to significant levels around -0.83 cpd, with adjacent frequencies at insignificant levels. This is of several low- one wavenumber modes that Alquist did not detect, but which exists theoretically. Eubanks [1993] notes that the 1.2 d and lOd modes are present in the atmospheric mass angular momentum because they are harmonically pure, and should excite polar motion at periods between and 20 days. In summary, at periods less than 4 days the 1 noise at certain and bands time series are practically except spectral lines corresponding to forced and free oscillations, where the power rises significantly above the continuous background power. There is good evidence that these time series are essentially linear versions of each other at periods longer than 4 days, up to periods of 2 years (the span of the 58 NASA data available to us). Because the polar motion data used in this study contain no reliable variations at periods less than 5 days, all the time series except %PLS produce essentially the same results when compared with the polar motion data. Thus the use of %PL in this research is justified, and at the same time the results can be confirmed using %PLN or . 59 Appendix C Neural network system identification This appendix describes a system identification filter which admits non-linear solutions in a general way, and which requires only minimal information about the functional form of the system. The filter is commonly described as an artificial neural network, but is more closely allied with statistics than with biological entities. The network filter successfully identifies a large class of non-linear systems. A system is something which takes an input signal and transforms it into an output signal. Audio amplifiers are systems which increase only the amplitude of a signal. Encryption machines are systems which change text into code in such a way that it is difficult to invert the code back into text simply by examining the code. A system is described by a filter, which is a special kind of function that operates on ordered data. A time series, which consists ofmeasurements of taken at some process a uniform time interval, is a natural example of such data. The task of system identification is to discover a filter that relates observations of the input and output of a system, and validate the results using input-output data that are not used to design the filter. One underlying assumption is that the system does not change. In certain applications, the system identification method should not require a statistically stationary input process. Indeed, we may wish to determine the of the system to an input process that does not much resemble the input response used to design the filter. On the other hand it is unreasonable to expect that process of a system to Gaussian white noise can be determined when the design the response data are sinusoidal. The amplifier is a perfectly linear system which is easy to 60 identify. A good encryption machine is a highly non-linear system which is difficult to identify. There are three main classes of applications of system identification: mapping, short-term prediction, and proof of relationship. Mapping is when the identified filter translates the input time series into the output time series with acceptable error at all times. Narendra and Parthasarathy [1990] discuss system identification and control in this context. Short-term prediction implies that the filter produces acceptable error for some finite time starting from correct initial conditions, after which time the prediction error increases. Lapedes and Farber [1987] and Tenorio and Wei-Tsih Lee a certain chaotic time series in this connection. Proof [1990] investigate prediction of error of relationship refers to a filter whose may be too large for mapping or prediction, but the filter significantly boosts spectral coherence for the validation data, thus showing that the input and output time series are related. This appendix is organized into four sections. The first is a brief introduction to linear system identification. The second is a description of the network filter and a simple method of using it for the system identification problem. The third is a set of demonstrations using specific filters and input data, and the last is the Fortran used. program Linear System Identification This appendix is not specifically about linear system identification, but the subject ground between linear and is a good starting point because there is much common network models in system identification. One approach to system identification is to that the filter is some linear combination of time-lagged inputs and outputs, assume i.e. 61 y(t) =c + a x(t) + djx(t 1)+...+b,y(t -1) + b y(t 2)+... la 2o where t is an integer time, x{t) is the input and y{t) is the output of the filter. This filter can describe derivatives and integrals, and more generally convolution and feedback, of a time series and hence is quite general. Observations of x(t) and y(t) leads to the matrix equation T x(l) x(0) y(0) y(-1) .¥ cl [y(l)~ 1 jc(2) x(l) y(l) y(0) a 0 y(2) . .. ... . <2, = 2 .. • ... . b, 1 x(n) x(n-1) y(n —1) y(n-2) b 2 y(n )_ . . or Gm =d. Practically a least squares solution is sought when the matrix is . overdetermined, and there is no linear combination of the columns which produces d The usual step taken is to find the linear combination of the columns which minimizes the 2-norm magnitude of the difference from d. The solution in=[c,aQia,...,bl,b]Tis usedtomodelthesystemas i2 - - +b 3 y(t)=c+dx(t)+ax{t I)+...bj(t-1) y(t 2)+... 2ox How well this works depends on many things, including: the linearity of the system; the number of columns forming G; the statistics of the input driving the 62 system; and noise in x(t) and y(t). For example, some non-linear systems are well described by a linear approximation when the input data are uniformly distributed noise in {o,l}, but are quite non-linear for input data that are zero-mean Gaussian white noise. Others behave quite linearly when the input is highly correlated, e.g. a sinusoid. If Gm = d is an appropriate model for the system, the number of time and lagged terms of x{t) y(t) can be determined experimentally. Noise, either measurement error or some other unwanted component, is always present in experimental data. As written above, the model treats noise and signal on an equal basis. Finally, the coefficients must be chosen so that the filter output does not grow uncontrollably. The coherence spectrum of the system input and output is a useful way of quickly assessing whether a linear model of (lb) will succeed in system identification. The coherence is a measure of the fraction of linearly related variance as a function of frequency. If the coherence spectrum is below the 99% confidence level across a broad range of frequencies, linear system identification will probably not succeed at reproducing the behavior of the system at those frequencies. Section 3 describes the method ofcomputing spectral estimatesusedinthis paper, andhowto interpret them. Another useful idea is the response of a system to an impulse. The impulse for a linear system completely describes the behavior of the system to any response input signal. This is not true for non-linear filters. A non-linear system that has an impulse response that rings forever in a complicated way is hard to identify. can be included in the linear In principle non-linear relationships, e.g. x(t)y(t -1), formulation. This works well for functions which primarily linear least squares are combinations of such terms. The linear method is fundamentally a parameter 63 estimation method for a known filter form, but for general system identification it is the filter function which requires estimation. Networkfilter This section describes a function which can approximate a large class of non­linear functions, and thus partly fulfills the requirement of being able to approximate the filter function rather than a set ofparameters to a known form. Define the function as F(Zj(t)) = f\YR b + F(z,-(0) = y(t) 3 kk VJ A: where f(w) =—-—4a l+ e~w =4b f[YQkia +AI K i vi J f\ fYpiJzjW+a> =a< 4c Vj J where P, Q, R, and a, ft, y are fixed parameters. The non-linear function f(w) 1 as w varies over the real numbers. As any has a sigmoidal shape between 0 and bounded function can be scaled and shifted into this interval, this is not a restriction. a, p, and determine the shift of the sigmoid along the w axis. The parameters y 64 variable the When the independent F(Zj(t)) takes form (x(t),x(t-1),...,y(r-1),y(7-2),...), then F(Zj(t)) is a filter. The original reason for using this form was to model a biological neural network. The idea is that each row of Ptj corresponds to a neuron whose input stimulus is z (t). f(w) is the activation function for the neuron, and the continues to the next process setofneuronscorrespondingtotherowsof Qki . The network model is a poor biological paradigm, but it has some interesting functional properties. F simply defines a surface in j dimensions. It is a non-linear mapping functionally composed of a network of identical elements, whose complexity is scalable by increasing the number of elements. Unlike the linear model, the parameters of each element have no easily identified meaning because they don't correspond in any simple way to zXt). The surface defined by F provides a way of explaining how the network model works in system identification. If the x(t) used to identify the system are representativeofthepossibleinputstothesystem,andFcanmap z}{t)accurately, then the surface so defined will act like a look-up table or interpolator for subsequent z (t). Even if these conditions are met, F may not represent the system at all. For - - be smooth continuous and differentiable but the example, the true surface may x{t) used to construct F may undersample it for subsequent x(t). For many non-linear systems, there exists an interpolation surface which is a good description of the system. If the system is linear, the best thing is to take something linear least methods to find the coefficients. For non-linear like (2), and use squares no reason to linearize, except that it gives the wrong answer a systems, there is million times faster than exploring non-linear solutions. 65 System Identification using F(Zj(t)) The most simple and least efficient method for determining the network parameters is based on a gradient search of the error = e5 (y(t)-y(t)f. Thus the gradient method described is a kind of least squares estimate. Although it is not a true steepest descent algorithm, it works well for the examples presented, allowing for the time required to perform hundreds to thousands of million-floating­point operations. It also can be written in a couple of pages of Fortran or C. To determine the the which network parameters, operational form of z (t), contains filter feedback terms y(t-l),y(t-2),..., is replaced with (x(t),x(t-1), ...y(t l),y(t -2)), which uses actual observations of the system output y{t). Thus Zj(t) corresponds to a row of G in (2). For a selected Zj(t) and y(t), the gradients de dedededede v——, t;——, v——, v—~, v— 6 dRdadpdy dP,i dQki k ik R and are subtracted from their corresponding estimates of P.., Qki , , fiy, k k, using a small step v. The Zj(t) are chosen randomly from the set used to design the is continued until filter and the process e is acceptably small for all the Zj{t) and y(t) pairs used. The gradient search does not guarantee a solution will be found, even if one exists, nor does it guarantee a solution will be optimal. And as already explained, 66 - a solution may not actually represent the system the crucial part of system identification is testing F with observations of x(t) and y(t) not used to construct F. In the following derivations, the gradients are expressed in terms of the observables Zj(t), y(t), and functions of them. Starting with (3), f\ w)=fXaa+r ¦andlettingr='LRA+r Vk) k de dedr d£ dedy . .~ , ... x_„ XI - === = ~~ 2iy =2(>> ? Yr^ Td YyTrK _ k which uses the chain rule, the fact that the network parameters are pairwise e~ independent,andthat/'=/(l-/)for/=l/(l+ w).Similarly, =2(y-y)y(\-y) 8 dy f\ = + , For (4b), b =f A and letting qt Y.Qkfl, +A. k Vi) ' d£ YL _ -­ Hi. f'<„ )a dQd y(t) cos(jc(0)+cos(y(t 1)) 2,20,10 x(t),y(t-l) =­ 2,30,15 x(t),y(t-1) The first filter in table 1 describes a system which takes the time derivative of a signal and applies a quadratic transfer function. Despite the linearity of the derivative, the filter is quite non-linear. For unit variance Gaussian noise, the coherence spectrum between x(t) and y(t) is below the 99% confidence level, meaning there is insignificant linear relationship (Figure 1.1). 70 riuer inpui ,„ i , 40, , mici uuipui , . -w—i* Figure 1 i1.1 a.First row: nniThe * nJljlklJilillijyb -2 -20 ,ilikll.l.iiLi HulUjJljLl.,.lil. iiiuß filter inPut is 500 independent Gaussian white noise samples 100 200 300 400 500 100 200 300 400 500 t. t with unit variance. The filter „Filter input jr~. J| 1 T Filter output T T y(t) = (x(t)-x(t-i)) 2 is the A# estimated squared time- Frequency Frequency ol a true power spectrum 1 Gaussian white noise I process ‘fmlm has equai power at ail 0 O204 0 °-204 The _ -frequencies. relatively Frequency Frequency J small variations in the estimated power spectrum are of the particular instantiation of white consequences noise, and the resolution of the spectral method. Third row: the spectral coherence is below the 99% confidence level of 0.5, and the phase apparently varies randomly in its domain of (-/r, n\. Hence there is no evidence of linearly related variance between x{t) and y(t) at anyfrequency. - The linear model y(t) = a x(t) + a x(t 1) fails utterly to describe the system, not ox , too surprisingly (Figure 1.2). A linear combination of the quadratic terms x(t) 2 x(t-\f, and x(t)x(t-1) correctly identifies the system, which only shows that the linear method can be used to identify the system when the system is known to within a linear combination. 71 IbUIIUi dflU llllcdl -« nuiudl IbUIICJ) allQ r (QdSMSQ) MClUdl icUj »—i• . r* J.—.—. ’ —; Figure 1.2 First row: after . ,— ,n 4 _ S-0-I _2.0.4 002 | 1 1 02 I designing the network filter F o:rJ jkLwdJ using the x(t) and y(t) data in \e 20406080100 20406080100 t t figure 1.1, the filter is tested 81| sl| with 500 independent A cc 0) 0) O5 O5 ia Gaussian white noise data, of ||(l | °° °o m which the first 100 are shQwn> Frequency Frequency Because the network. filter .has I2 ' |2 I All y \ \ a range of (0,1), the results are Q in this °-_2-| |\A I displayed range, ° °4° °4 although they Can be Scaled Frequency Frequency and shifted back to the proper domain. The network filter is identical to the actual filter. The linear model - = aa y(t) x{t)+ x(t 1)completelyfailsthetest.Explicitinclusionofthequadratic () x terms The is required for the linear method to identify the filter. Second row: coherence spectrum for the network filter F is unity, while that of the linear model remains statistically insignificant. Third Row: the filter F has zero phase, and the linear model phase is random, supporting the interpretation ofcoherence. Figure 1.3 The network is tested with a sinusoidal sweep, depicted in the top graph. The actual filter output on the lower graph is well estimated by the network filter, even of correlated though this sort signal is different from the Gaussian noise design data. Figure 1.4 The difference between the actual filter and the network filter, for both the design data (oS little improvement in Frequency Frequency the coherence and phase spectra from the design data. In contrast, the network filter correctly identifies the system when the input data are unit variance Gaussian noise (Figure 2.2). Not only is the rms error small, but the spectral coherence and phase are nearly perfect. Starting from an initial random set of most of the variance reduction network parameters, is accomplished in the first 100,000 iterations, or a few minutes on a machine that executes 1,000,000 floating point operations per second. After 1,000,000 iterations there is essentially no further variance reduction (Figure 2.3). 75 Figure 2.3 The total rms error between the network filter and the actual filter design data reaches nearly its final values during the first 100,000 iterations. As shown in the fourth example, there is no way to decide of the basis of such a graph as to when to cease iteration. A frequency-sweep sinusoid input is also correctly identified, even though this sort of correlated signal does not resemble the Gaussian noise process used to design the filter (Figure 2.4). Figure 2.4 The sinusoidal test is well identified by sweep the network filter F, even though this sort of correlated input does not resemble the Gaussian noise design data. Most of the rms error in this case is overshot along each sinusoid. Qualitatively this are error is consistent with the idea of an interpolation surface. The design data, which 76 uncorrelated and independent, are not likely to map in any detail the interpolation surface for a correlated signal. As the sinusoid reaches higher frequencies, the rms error decreases. Also note that the sinusoid input produces a small dynamic range from the filter relative to the range from the noise input design data. Hence the rms error for the sinusoidal part of the signal is a larger fraction of the dynamic range than it is for the design data, although the total rms error of each segment is similar. The format follows Figure 3.1 that of figure 1.1. Both x(t ) - and y(t) = x{t)x{t 1) +x(t l)cos(/ry(t 2)x(t)) The have fairly white spectra. time series are linearly related near the frequency of 0.2, based on the stable phase and increased coherence, but outside this band there is no evidence for relationship. The third filter in Table 1 an is an interesting non-linear example because it has with a 1. Because terms like apparently perfect impulse response time-lag of x{t)x{t -1) dominate the filter, their inclusion in the linear method identifies the system as well as the network method (Figure 3.2). Both linear and network methods similar to those discussed in the fourth fail the frequency-sweep test for reasons example below. The network method accomplishes most ofthe variance reduction in the first 100,000 steps. Both over-parameterization and a larger network produce identical results. Thus for both linear and network filters, identification only succeeds when the input time series is statistically similar to the Gaussian noise design data. 77 rvuiuai NiitJcH \—. \uuomcu/ aiiu _ rT^l r r 11 The format follows T -r—— Figure 3.2 „ r—¦ 1 |'—p —P °0 S o that of figure 1.2. The network 11 _n J . . ' 11 g I filter boost coherence above 20406080100 20406080100 1 t the 99% confidence level and ~ * 81| gl| produces nearly zero phase, as |o,VVWV*A does the linear model which o o o o — °o 0.2 °o includes x(t), x(t-1), 0.4 0A Frequency Frequency and all y(t —2), possible 2P ® The products of these terms. | o “¦-2 ¦ product terms dominate the -2 ° °4 ° ° 4 filter and henCe 3 linear Frequency Frequency ’ combination of these terms is a goodestimateofthefilter.Thenetworkfilterisbased on x(t),x(t-l), and y(t-2) only. Despite its simple form, the fourth example is the most difficult to identify of the filters listed in Table 1. The linear method fails, even when the cubic terms x(t)3 and 3 y(t -1) and fifth-order terms are supplied. The network filter produces significant variance reduction, boost in coherence, and stable phase enough to prove that the - two time series are related. The frequency sweep fails for low frequencies, but the agreement becomes better at high frequencies, where the rapid change more closely approximates the Gaussian design data. A plot of the frequency-sweep input and output demonstrates why this filter is difficult to identify (Figure 4.3). The input signal is smooth and correlated, but the output sporadically develops instabilities, characterized by large rapid fluctuations. The corresponding interpolation surface must also fluctuate rapidly over part of its domain. Hence this surface is difficult to especially with a Gaussian noise input, because the parts of the map map, corresponding to smooth changes of input signal are rarely encountered. Increasing 78 the number of network parameters does not decrease the rms error or increase coherence. The variance reduction is different from the previous examples, where most of the reduction occurs within the first 100,000 iterations. In this example, the reduction halts abruptly, but then continues suddenly after about 3 million iterations (Figure 4.4). This demonstrates that the gradient method used to estimate the network parameters is not optimal, and provides no indication as to when further iteration is fruitless. Figure 4.1 The format follows 1.1. Both that of figure x(t) and - = y(t) cos(x(t))+cos(y(t 1)) have fairly white spectra, and evidently there is no linear relationship between them. 79 r-\oiuai ~ lounui emu r MUIUdi lt>unui cliiu iintJdi ludsiieu; r r* 11 The format follows 04fr—Figure4.2 that °f fisure 1 2 The netw°rk 0.4f—i—r—i.—P—_ t°ititi*(I'iirifiy t°— . o -­ f'°!H *IT 'WpiyWl o filter p resembles the actual 20406080100 20406080100 t t filter TT TT = 81sly(0 cos(*(0) + cos(y(f -1)) || i‘pM ijUyJJ r.r;zi:r:: Frequency Frequency that F phase spectra prove partly ldentifies the system- I i 1 l/lilltl II Hh Ifl °-_2 °~-2 J 111/ The linear model, based on |jw| \ 0 °-2 04 0 °-2 04 x(t), y(t —1), and odd powers \ s\/ r Frequency Frequency ’ to 5, does not identify the filter, either qualitatively on the first row graph, or as measured by coherence and Figure 4.3 The total rms error between x(t) and y{t) during design of F. After several million iterations the error ceases to decrease, but shortly thereafter drops precipitously. The gradient method presented is not optimal, and provides no information when to cease as iteration. 80 21——-r T T t Figure 4.4 The first couple of 15 _ \ . | cycles of the frequency sweep I , test illustrate why this filter is 0.5­0 i i.' i lU'l' / I y^\I / difficult to identify. At first filter output is smooth and correlated like the input, but it -0.5­ l\ | f\ develops instabilities part-way i 1 I through the first sinusoid cycle <1 J U before smoothing out. These I I instabilities are set off by 0 20 40 60 80 100 minute changes in the slope of x(t). Phase space diagrams for a number of sinusoids show the instabilities are most complicated for low frequency sinusoids, and collapse into simpler states for higher frequencies. Conclusions The network filter is capable of identifying a variety of non-linear systems. For almost-linear systems the network filter performs at least as well as finding linear combinations of products and powers of the input data zXt). However, no significant change in coherence or phase is observed with the linear model for the filter - y(t)=cos(x(t))+cos(y(t 1)),whilethenetworkmodelcoherenceandphasespectra reveal that x(t ) and y(t) are related. This example suggests that the network filter is more robust for problems where the goal is to ascertain whether two time series are functionally related. The main disadvantage of the network filter is the slow algorithm for finding the network parameters. Faster algorithms are available, even for the steepest descent procedure. A sensible strategy in system identification is to try the linear method first before pursuing the network solution, especially to give some insight about the form of Zj(t). 81 program nnw c Language: completely standard Fortran 77. Channels 5 and 6 are c the standard input and standard output respectively. c Author: John Kuehne. Copyright Feb 1 992 John Kuehne c Neural network with two Description: hidden layers, quasi c i.e. steepest descent using word training parameters are updat­ c ed after the presentation of each pattern, not after the presen­ c tation of all patterns. c Directions: The initial matrices, patterns, and targets are c read and written as a single list of ascii numbers c without break. The dimensions are fixed at compile time via the c Parameter statement. The output matches the format of the in­ c put so the program is easy to restart. The idea is to use c Matlab as a front-end to setup the problem and examine the c results. Re: q 1 input lines, q 2 first layer nodes, q 3 second c layer nodes, q 4 output lines, q 5 patterns. For reading and writ-c ing the list of numbers, matrices are stored as the consecutive c column vectors, followed by the patterns which consist of con- c secutive lists of q 5 numbers for each of the q 1 input lines, c followed lastly by the targets which are written like the pat- c terns i.e. lists of q 5 numbers for each of the q 4 output lines. integer ql, q2, q3, q4, q 5 = parameter(q1=10,q2 =40,q3 =20,q4 =2,q5 99) real w 1 (q2,ql), w2(q3,q2), w3(q4,q3) real t2(q3), t3(q4) (qt12), real (qo12),02(q3), 03(q4) real pattern(ql ,q5), target(q4,qs) real d3(q4), d2(q3), dl (q2) real dwl(q2,ql), dw2(q3,q2), dw3(q4,q3) real dtl (q2), dt2(q3), dt3(q4) integer i, j, k, iter, pass!, pass2, s real n c Read initial matrices, pattern, and target from standard input, do 1j= 1,q1 do2i=l,q 2 read(s,*) w 1 (i,j) 2 continue 1 continue = do3j 1,q2 do4i=l,q 3 read(s,*) w2(i,j) 82 4 continue 3 continue do5j 1,q = do 6 i=l,q read(s,*) w3(i,j) 6 continue 5 continue do 7 i=l,q read(s *) tl(i) 7 continue do 8i=l,q read(s,*) t2(i) 8 continue do 9 i=l,q read(s,*) t3(i) 9 continue do 10j 1,q = do 11 i=l,q read(s,*) pattern(j,i) 11 continue 10 continue do 12j=1,q do 13i=l,q read(s,*) targetQ,i) 13 continue 12 continue iter = 50 n 0.1 = = do 14 passl 1 ,iter do 1 5 pass 2 = 1,100000 c select pattern 3 4 2 3 4 1 5 4 5 = s nint(rand(o)*float(qs-1 ))+1 First matrix-vector multiply (wl*pattern) and sigmoidal transfer, do 20 j = ol Q) = do 30 k ol Q) 1,q 2 wl GJ ) * pattern(l ,s) = 2,q1 = ol (j) + (wl (j,k) * pattern(k,s)) 30 continue ol G) = 1.0/ (1.0 + exp(-1.0 * (ol G) +tl G)))) 83 20 continue c Second matrix-vector multiply. Same as above, do 50 j=l,q 3 o2(j) = w2(j,l) *ol (1) = do60k 2,q2 = o2(j) o2(j) + (w2(j,k) * ol (k)) 60 continue = o2(j) 1.0 / (1.0 + exp(-1.0 * (o2(j) + t2(j)))) 50 continue c Third matrix-vector multiply. Same as above, do80j= 1,q4 = o3(j) w3Q,I) * o2(1) = do90 k 2,q3 = o3(j) o3(j) + (w3(j,k) * o2(k)) 90 continue = o3Q) 1.0 / (1.0 + exp(-1.0 * (o3Q) + t3Q)))) 80 continue c Backpropogate errors, = do 110j 1,q 4 = ­ d3Q) (target(j,s) o3Q)) * (o3Q)) * $ (1.0-o3(j)) 110 continue = do 130j 1,q3 d2(j) w3(l ,j) * d3(1) = = do 140 k 2,q4 = d2(j) d2(j) + (w3(k,j) * d3(k)) 140 continue - = d2(j) d2(j) * o2(j) * (1.0 o2(j)) 130 continue = do 155j 1,q 2 dl 0) = w2(l j) * d2(1) do 165 k 2,q3 = =+ dl (j) dl (j) (w2(k,j) * d2(k)) 165 continue = dl Q) d1(j)*o1(j)*(1-0-o1(j)) 155 continue c Estimate offset perturbations, do 170j 1,q 2 = = 0.0 dt 1(j) 170 continue do 190j 1,q 2 = = dt1(j) dt 1Q)+dlG) 190 continue 84 = do200j 1,q2 = t1(j) t1(j) + (dtl(j)*n) 200 continue do 210j= 1,q 3 dt2(j) = 0.0 210 continue = do230j 1,q3 dt2(j) = dt2(j) + d2Q) 230 continue do 240 j=l,q 3 = t2(j) t2(j) + (dt2(j) * n) 240 continue = do250j 1,q4 = 0.0 dt3(j) 250 continue = do270j 1,q4 dt3Q) = dt3G) + d3G) 270 continue = do280j 1,q4 = t3G) t3G) + (dt3G) * n) 280 continue Estimate weight perturbations for selected pattern do290j 1,q3 = = do300 i 1,q4 = dw3(i,j) d3(i) * o2G) 300 continue 290 continue = do292j 1,q 3 do293 i 1,q4 = = w3(i,j) w3(i,j) + (dw3(i,j) * n) 293 continue 292 continue do310j= 1,q2 do320 i = 1,q 3 dw2(i,j) d2(i) * ol G) = 320 continue 310 continue do312j 1,q2 = = do313 i 1,q3 = w2(i,j) w2(i,j) + (dw2(i,j) * n) 313 continue 312 continue = do330j 1,q1 do340 i 1,q2 = 85 = dwl (ij) d 1 (i) * pattern(j,s) 340 continue 330 continue do332j= 1,q1 do333 i = 1,q 2 = wl(ij) w1(ij)+(dwl(ij)*n) 333 continue 332 continue 15 continue 14 continue c c Write final matrices to standard output as a single list c = do400j 1,q 1 do 405 i=l,q 2 write(6,*) wl (i,j) 405 continue 400 continue do 410j= 1,q 2 do41 5i=l,q 3 write(6,*) w2(i,j) 415 continue 410 continue = do420j 1,q 3 do 425 i=l,q 4 write(6,*) w3(i,j) 425 continue 420 continue do 430 i=l,q 2 write(6,*) t 1(i) 430 continue do440i=l,q 3 write(6,*) t2(i) 440 continue do 450 i=l,q 4 write(6,*) t3(i) 450 continue write initial pattern and target as a single list c = do 460 j I,ql do470i=l,q 5 write(6 *) pattern(j,i) 86 470 continue 460 continue = do480j 1,q 4 do 490 i=l,q 5 write(6,*) target(j,i) 490 continue 480 continue end Appendix D Annual component ofpolar motion Throughout this dissertation emphasis has been placed on the continuous part of polar motion spectrum and special care has been taken to remove annual and semi­annual sinusoids. This section briefly describes the annual variation. The annual motion is of the pole a prominent part of polar motion. It was correctly understood by early workers as arising partly by accumulation of cold dense air over the northern hemisphere each winter [Munk and MacDonald 1960]. The , importance of the annual and Chandler frequencies in the literature is partly an artifact ofthe old optical data, which are reliable only at frequencies near the resonant reliable Chandler frequency of about +0.84 cpy (435d). Modern data are at periods between days and decades and the emphasis of modern investigations is no longer centeredon theannuallinespectrum,whichis describedbyjustfournumbers. 1 The annual component of polar motion is defined here as the cpy line spectrum of polar motion. It can be described by the sinusoidal variation along the real and imaginary axes of the geographic reference frame m{t) asin(cot+(]))+iftsin(cot+0(.) 1 = r and which describes an ellipse, where t is time, co is frequency, the remaining constants specify magnitude and phase. Using the identity = sin(x + y) cosO)sin(y) + sin(x)cos(y) (1) expands to the linear combination m(t) acos(cot)+bsin(cot)+i[ccos(cot)+dsin(cot)]. =2 The four parameters a,b,c,d can then be determined by the least squares projection of onto and the polar motion time series cos(mO sin(cut) for the real and imaginary 88 components. However, this representation is tied to the rather arbitrary coordinate system of Greenwich. More geophysical meaning can be had by projecting the polar motion time series onto prograde (west-to-east) and retrograde (east-to-west) circular components, i.e. im + Be~icot Ae3 . The complex constants A=A, +i\ and 5 = 5,+ iB2 are called the prograde and retrograde vectors. The magnitudes of these vectors are the same regardless of the choiceof0longitude,with |A|+ls|thelengthofthesemimajoraxisoftheellipse,and ||A| -15|| the length of the semiminor axis. The relationship between the Cartesian parameters a,b,c,d and the complex quantities A and 5 is A, ={a+d) /2, A,=(c-b)/ 2 ={a-d)l 2 5 =(b+c) /2. 4 B x2 The prograde and retrograde components are simply the coefficients of the discrete Fourier transform at ±1 cpy. The discussion above thus carries over to calculating the continuous part of the spectrum, producing the positive (prograde) and negative (retrograde) frequencies of spectral plots. The power spectrum, computed only from the magnitudes of the prograde and retrograde coefficients, is thus a coordinate-free representation. Using the old optical data and global barometer data, a number of studies have shown a significant discrepancy at 1 cpy between polar motion and atmospheric inverted barometer variations with an ocean. The prograde atmospheric vectors of Kuehne and Wilson [1990] and Wahr [1982], both based on historical data, are about 100°, which is considerably west of polar motion (65°), while the magnitudes are slightly larger than polar motion. The emphasis on prograde vectors in these studies using optical data is a result of the amplification near the prograde resonant Chandler 89 wobble. The retrograde component, being much farther from the Chandler frequency, is more susceptible to noise. The conclusion of Kuehne and Wilson [1990] is that there is another excitation source, and the results of the present study using the ECMWF data and space93 for the period 1985-1993 support that conclusion. The vectors for which is ¦> essentially an inverted barometer model, are comparable to those of Wahr [1982] and Kuehne and Wilson [1990]. Remarkably, the linear combination of land and SH ocean, described in Chapter 2, nearly matches the prograde and retrograde polar other motion vectors in both phase and magnitude, further suggesting that the excitation source is over the oceans, either as a non-inverted barometer response to pressure forcing, or some other effect of the winds. Table 1 Prograde Retrograde Phase Phase Study Magnitude Magnitude Wahr [1982] (1900-1958) -100 10.2E-8 -94 10.0E-8 KuehneandWilson [1990](1900-1985) -100 8.6E-8 -94 7.9E-8 X,FL (1985-1993) -101 6.0E-8 -106 6.7E-8 P°' (1985-1993) 70 3.3E-8 -149 3.9E-8 X -80 5.3E-8 -110 6.5E-8 X? (1985-1993) -70 6.9E-8 -119 3.3E8 -64 7.7E-8 -141 9.4E-8 X"' space93 (1985-1993) ILS (1900-1985) Note: phases are relative to January 1. 90 Bibliography Barnes, R.T.H., R. Hide, A.A. White, and C.A. Wilson, Atmospheric angular momentum fluctuations, length-of-day changes and polar motion, Proc. R. Soc. London Ser. A, 387, 31-73, 1983. Chao, 8.F., A.Y. Au, Atmospheric excitation of the Earth's annual wobble: 1980­1988, J. Geophys. Res., 96 (B4), 6577-6582, 1991. Eubanks, T.M., J.A. Steppe, J.O. Dickey, R.D. Rosen, Salstein, Causes of D.A. rapid motions of the Earth's pole, Nature, 334 (6178), 115-119, 1988. Eubanks, T.M., Variations in the orientation of the Earth, Contributions ofSpace Geodesy to Geodynamics: Earth Dynamics, Geodyn. Ser., vol. 24, edited by D.E. Smith and D.L. Turcotte, pp. 1-54, AGU, Washington, D.C., 1993. observations Gross, R.S., Correspondence between theory and of polar motion, Geophys. J. Int., 109, 162-170, 1992. Gross, R.S., U.J. Lindqwister, Atmospheric excitation of polar motion during the GIG '9l measurement campaign, Geophys. Res. Lett., 19 (9), 849-852, 1992. International Earth Rotation Service, Annual report for 1993, Observatoire de Paris, Paris 1993. Jeffreys, H., Causes contributory to the annual variation of latitude, Mon. Not. R. Astron. Soc, 76, 499-524, 1916. Jeffreys, H., The variation of latitude, Mon. Not. R. Astron. Soc., 100, 139-154, 1940. Jeffreys, H., The variation of latitude, Mon. Not. R. Astron. Soc., 141, 255-268, 1968. 91 Kuehne, J., C.R. Wilson, Terrestrial water storage and polar motion, J. Geophys. Res., 96,4337-4345,1991. Kuehne, J., C.R. Wilson, and S. Johnson, Atmospheric excitation of nonseasonal polar motion, J. Geophys. Res., 98, 19,973-19,978, 1993. Lapedes, A., R. Farber, How neural networks work, Proceedings of lEEE, Denver conference on neural nets 1987 (submitted). Los Alamos Nat. Lab. preprint LA­ UR-88-418. Lapedes, A., R. Farber, Nonlinear signal processing using neural networks: prediction and system modelling, Los Alamos National Laboratory, Lost Alamos, NM 87545. Narendra, K.S., K. Parthasarathy, Identification and control of dynamical systems using neural networks, lEEE Trans, on Neural Networks, vol. 1, no. 1, March 1990. Preisig, J.R., Polar motion, atmospheric angular momentum excitation and earthquakes - correlations and significance, Geophys. J. Int, 108, 161-178, 1992. Tenorio, M., W. Lee, Self-organizing netowork for optimum supervised learning, lEEE Trans, on Neural Networks, vol. 1, no.l, March 1990. Thompson, D.J., Spectrum estimation and harmonic analysis, Proc. of the lEEE, 70 (9), 1055-1096, 1982. Trenberth, K.E., J.G. Olson, An evaluation and intercomparison of global analyses from the National Meteorological Center and European Centre for Medium Range Weather Forecasts, Bull. Am. Meteorol. Soc., 69 (9), 1047-1057, 1988. Wilson, C.R., Discrete polar motion equations, Geophys. J. R. Astron. Soc., 80, 551­554, 1985. Wilson, C.R., R.O. Vicente, An analysis ofthe homogeneous ILS polar motion series, Geophys. J. R. Astron. Soc., 62, 605-6161, 1980. 92 Wilson,C.R.andR.O.Vicente,Maximumlikelihood estimatesofpolarmotion parameters, Variations in Earth Rotation, Geophys. Mono. Ser., vol. 59, edited by D.D. McCarthy and W.E. Carter, pp. 151-155, AGU, Washington, D.C., 1990. 93 The vita has been removed from the digitized version of this document.