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)cos(j)
(C A)gJ s
where 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 response, or some other consequence of the winds, may be responsible for the correlation. Of all the oceans, only pressure fluctuations over the southern hemisphere region is significantly correlated with polar motion. Although the three most significant regions, 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 atmospheric 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 Jpsinocos>e'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 nonlinear 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-floatingpoint 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.50 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 semiannual 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: 19801988, 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, 551554, 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.