A new technique in demodulation of normal modes

PHYSICS OFTHE EARTH AND P LAN ETA RY INTERIORS

_________

ELSEVIER

Physics of the Earth and Planetary Interiors 84 (1994) 139—160

A new technique in demodulation of normal modes Ramin Nawab

Philippe Lognonné

Laboratoire de Sismologie, Institut de Physique du Globe de Paris, 4 place Jussieu, 75252 Paris Cedex 05, France (Received 6 May 1993; revision accepted 8 December 1993)

Abstract Lateral heterogeneities of the Earth produce amplitude modulations of the so-called unresolved multiplets. A spectral fitting technique based on the polynomial interpolation of Hermite allows the retrieval of these slowly varying modulation functions. No restrictive hypothesis on coupling among multiplets is necessary. Modulation functions are constrained in time within a Q-cycle for data with typical signal-to-noise ratios. Sensitivity to noise is reduced by the introduction of cosine tapers and by a priori information on the modulation functions. The method then damps rapid oscillations caused by noise and becomes stable. These amplitude modulations can be used in a first stage for inversions where both phase and amplitude are necessary, and generalize basic mode observables such as the local frequency.

1. Introduction Lateral variations of the deep structure of the Earth can be retrieved by non-linear and simultaneous inversions of a large set of spectra of normal modes (see, e.g. Ritzwoller et al., 1986, 1988; Giardini et al., 1987, 1988). Such methods imply, however, theoretical approximations on the normal modes, for instance by neglecting coupling effects. This is also the case for stripping or stacking methods (see, e.g. Dziewonski and Gilbert, 1973a, b; Gilbert and Dziewonski, 1975; Buland et al., 1979), which extract from a global set of normal modes the eigenfrequencies under the isolated multiplet hypothesis. Other inversion methods, in contrast, extract in a first step some observables from the seismic _______

*

Corresponding author,

signals and then invert these observables to retrieve the aspherical structure. The simplest observable is the frequency of the maximum of an unresolvably split normal mode resonance peak and its deviation from the spherical frequency of the multiplet. As shown by Jordan (1978), this deviation, known as the location parameter (A), is, up to the first order of the perturbation theory of the Earth’s normal modes (see, e.g. Dahien, 1979) and for times shorter than the inverse of the frequency shift (short time approximation), related to the depth-averaged structure between the source and the receiver. Inversions of these observables have been performed by Silver and Jordan (1981), Masters et a!. (1982), Roult and Romanowicz (1984), and others. More exact expressions have been developed by Dahlen (1974), Woodhouse and Girnius (1982), Davis and Henson (1986) and Romanowicz and Roult (1986), including either coupling effects or higher-order

0031-9201/94/$07JJ0 © 1994 Elsevier Science B.V. All rights reserved SSDI 0031-9201(94)05030-2

140

R. Nawab, P. Lognonné/Physics of the Earth and Planetary Interiors 84 (1994) 139—160

asymptotic effects. However, the relationship between such observables and the Earth’s lateral variations remains based on a theory valid only for the first hours of the seismic signal, and for high angular order surface waves. Therefore, this approximation cannot be used for long-period, low angular order and high-Q normal modes. To overcome the short time approximation, Smith and Masters (1989) generalized Jordan’s location parameter, introducing a time-dependent location parameter for which Jordan’s A was only an initial value. Such an approximation is better than a simple time-independent location parameter, especially for the fundamental modes which have a low quality factor. However, it fails again around 10 h. Therefore, all these attempts have not been able to model the time domain envelope of the multiplets for large time intervals. The estimation of the multiplets’ envelope is necessary, however, to extract from the spectrum all the effects of the Earth’s asphericity, especially the coupling effects, either to first order, on the amplitude, or to second order, on the eigenfrequencies. They are also necessary in the retrieval of anelasticity lateral variations, where both phase and amplitude deviations are necessary. A first attempt towards a direct estimation of the modulating envelopes has been the introduction of the ‘complex demodulation’ technique (Tukey, 1961) in analysis of normal modes (Bolt and Brillinger, 1975, 1979; Hansen, 1980; Hansen and Bolt, 1980; Hansen, 1982a, b). This technique has been used to make a preliminary estimation of the mean complex frequency, amplitude and phase of multiplets before a non-linear regression scheme in the frequency domain for a more accurate estimation of these parameters. It has also been used as a tool in assessing the quality of the data. Indeed, the existence of instabilities in the demodulation is a clear sign of the presence of noise or of singlets beating, implying the violation of the adopted basic hypothesis of a simple resonant peak model for the multiplets. ‘Complex demodulation’ has also been used to estimate the bias in attenuation measurements caused by beating or by other phenomena such as the presence

of secondary sources of mode excitation (such as aftershocks; see e.g. Hansen and Bolt, 1980). Hansen (1982b) has reported a systematic use of this technique to calculate the mean complex frequency, amplitude and phase, and the corresponding statistical uncertainties for some fundamental spheroidal (53) and toroidal (13) multiplets. As, in the last decade, the quantity of digital very long period data has increased, and their quality has improved, it is now possible to use more sophisticated demodulation techniques to model the fine structure of the multiplets and the coupling effects, and hence to provide new constraints upon the Earth’s deep 3D structure. Park (1990) has used such a technique to perform a direct estimation of the envelope function A(t) for some coupled multiplets. This estimation is based on a multi-taper spectral analysis, using the n~rprolate tapers of Thomson (1982), adapted to studies of the Earth’s normal modes (see, e.g. Lindberg, 1986; Park et a!., 1987). The result of this ‘envelope inversion’ is a modulation estimate, A, that has minimum norm (i.e. II A II 2 minimum) and fits the K linear constraints arising from the application of the K first n~rprolate tapers to the time series. A further improvement to this ‘envelope inversion’ technique has been provided (see, e.g. Park, 1992; Park and Maasch, 1993) with, this time, an envelope estimate having a minimum second derivative (i.e. II D~~A 2 minimum, D being a first-difference operator). The retrieved modulation is the smoothest model for the oscillation envelope A(t). In the same perspective, we develop a technique to obtain direct estimated parameters of the modulation function. It consists in a damped linear least-squares cubic approximation of the envelope, made at regular time intervals, which leads to a set of modulation coefficients. The introduction of tapers and the use of a priori information on the envelope make this technique much more stable than the classical ‘complex demodulation’. This is also the main difference from the ‘envelope inversion’ described by Park (1992) and Park and Maasch (1993), which uses a simple roughness penalty. Another difference from the ‘envelope inversion’ technique is that,

R. Nawab, P. Lognonné/Physics of the Earth and Planetary Interiors 84 (1994) 139—160

Evolution of 0S24 ( M84C + L0256 ( real part: first order)

—

141

)

a

~t) 0~.

r.norinafl..d (Woodhouse. 1983) cubic approximation (Smith & Master.. 1989) .—. cubic short tim. approximation —

00

5

—

10 time

(

hours

)

Evolution of 0S24 ( M84C + L0256 (imaginary part: first order)

time

20

15

)

( hours )

Fig. 1. A comparison between three approximations to the exact real (a) and imaginary (b) parts of the modulation function and the true modulation for the 0S~multiplet computed for the Mexico earthquake of 19 September 1985 at Station ALE (IDA network), using Eqs. (2) and (3). The continuous curves represent the evolution of the exact function. The dashed curves represent the approximation by Woodhouse (1983), which consists of modelling the modulation by a sinusoidal having the same initial slope as the exact function. The dotted curves represent the approximation made by Smith and Masters (1989), which is a variant of the previous one with, this time, a sinusoidal having a time-dependent frequency. The remaining curves are the Taylor expansions of the modulation up to third order made at the origin. As can be seen, this last approximation fits well the exact modulation during the first 3 h of the signal. Our method is also based on such a cubic approximation renewed at regular time intervals. (Figure from Lognonné (1989).)

142

K Nawab, P. Lognonné/Physics of the Earth and Planetary Interiors 84 (1994) 139—160

instead of n’rr prolate tapers, classical tapers are used here. Even if the concentration of multiplet energy is less efficient, this does not affect the amplitude demodulation, at least for the frequency isolated normal modes. The retrieved modulation coefficients can be related to the eigenfrequencies and eigenfunctions of normal modes (the relationship being non-linear) and constitute for each mode and each source—receiver pair an observable of the seismogram. These observables can be used in a second stage to invert for the lateral variations of the Earth’s structure. We start by introducing the method, and show how modulation coefficients are related to the aspherical structure. Tests on synthetics are then presented, showing the stability of the method with respect to noise. We finally present results for some low angular order normal modes. We start with a description of the method.

2. The modulation function

.

Let us start with the expression of the surface acceleration after an earthquake with a source time function approximated by a Heaviside function, as given by Lognonné (1991): u(t)

=

Re~M: Vvk(rS) exr(w-kt)uk.

(1)

The summation is done over all normal modes k, 0’k are, respectively, the primal where Uk, Vkthe anddual eigenmode and their correeigenmode, sponding eigenfrequency. M is the seismic moment tensor. Eq. (1) can be rewritten as a double summation, the first one on all multiplets K of the Earth’s eigenspectrum and the second one on all singlets k of a given multiplet: u( t)

=

Re

~ AK( t)

exp(urK t)

(2)

K

where AK(t)

=

~

In Eq. (3), the summation is limited to all the singlets k of the multiplet K. The carrier exp(io~t) represents the fast variation, oscillating at the complex frequency o~, and AK(t) is the slowly varying modulation function of the multiplet, which includes all the effects related to lateral heterogeneities and depends also on the deviation of the mean multiplet frequency with respect to the complex frequency o-~used in the demodulation procedure. The spherical mean attenuation is then ineluded in the carrier wave which has a complex frequency. Therefore, in a spherically symmetric non-rotating anelastic isotropic (SNRAI) Earth model, the modulation function for a given multiplet will be a constant. Many alternatives have been put forward to describe the evolution of the modulation function AK(t). A first approach, based on the Born Series, was to develop the exponential terms in a power series, and evaluate the coefficients of this expansion in time (e.g. Tromp and Dahlen, 1990). The counterpart in the frequency domain is the decomposition of the spectrum in terms of higher-order moments, as defined by Jordan (1978). Unfortunately, this series converges too slowly. Other methods based on the renormalization technique have been developed, first of all by Woodhouse (1983). They have been extended to a time-dependent local frequency by Smith and Masters (1989). Fig. 1 shows a comparison between these three approximations and the true modulation function for the mode 0S24, for lateral heterogeneities given by the models M84A (Woodhouse and Dziewonski, 1984) and L02.56 (Dziewonski, 1984), obtained by using Eqs. (2) and (3), for the Mexico earthquake of 19 Septemher 1985. It is clear that these attempts do not fit the modulation function after 10 h, or if they do, they may not represent the true first moments of the multiplet. However, the modulation function is of great interest. Indeed, it separates the effect of lateral

M: Vvk(rS) exp[i(o’k—o~)t}uk.

kEK

(3)

heterogeneities from the spherical averaged structure which appears in the fast oscillating part exp(io~t). Hence it crystallizes all the information on the aspherical structure of the Earth and can be used in a second stage as a secondary

R. Nawab, P. Lognonné/Physics of the Earth and Planetary Interiors 84 (1994) 139—160

observable. The numerical method described in this paper allows us to extract this function from the data. It is based on the fact that a polynomial approximation of the modulation function is accurate for several hours (Fig. 1). The modulation function can then be approximated by a cubic interpolation at regular time intervals. With such an approach, we describe the modulation function by its amplitude A(1,) and derivative A(1,) at times 7, = (n 1) x AT, AT being the sampling interval of the modulation function. We have AT>> ~t, where 5t represents the time spacing between data points in the initial time series. In practice, the sampling rate per modulation (spm) varies from two to six. For a multiplet such as 0S6 which has1 a typical splitting width of ±5%c,using four samples per modulation leads to a AT=1/(4XAf~~11~) of the order of 13 h. Here Af5~1~~ represents half of the maximum splitting width of the multiplet (i.e. 5.2 p~Hz).We will start our discussion with modes which are isolated in the frequency domain and whose spectral peaks are not strongly affected by the flanks of the neighbouring resonant peaks. However, the proposed demodulation technique uses no restrictive theoretical hypothesis on the coupling between modes. Coupled modes can then be demodulated with the same technique, if the inversion of the modulations is carried out simultaneously for all relevant modes. Let us now consider one multiplet of the summation in Eq. (2) and its contribution to the signal, written as —

SK(t)

=

Re[AK(t)

exp(ioK0t)}

(4)

We omit hereafter the K subscript identifying the multiplet. If we cut S(t) into N1,,~segments, each of length AT, starting at T,, and ending at T,, + ~, its Fourier transform can then be written as

143

t

1,~ and t2~are, respectively, the times for the first and the last points in the nth segment. Let us now express the modulation A~(t)for the nth sequence as a cubic polynomial,

A~(t) where

=

a,,

+

f3nTn

+ y,,T~+ 6~T~

(7)

1,) (8) 1,) 13n’ Yn and c5,, are complex coefficients and an, can be expressed in terms of the amplitudes A(1~)and derivatives A(I~)as follows: a,, = A( 7~) (t

~

=

~

—

13,,

7,)AT y=3A(T+~) —3A(1~) AT[A(T) + 2A(T,,)] ~,,

=

(9)

AT[A(7~±1)+A(I~)] 2A(1~÷~) —

+2A(I~) It should be noted that the modulation function and its first derivative are here continuous. This interpolation method is known as the polynomial interpolation of Hermite. Other kinds of polynomial interpolation can be used. For example, we can chose to use cubic splines to interpolate the modulation function, putting implicit constraints on its second derivative. The analytical expression for the Fourier spectrum for the nth interval can then be easily expressed as S,,( w)

=

exp[ i(o’0 w) t1,,] X[A(~,)d~(w) +A(7,~1)d~(w) —

+A(I,)AI,d~(w)

Nh,,

S(w)

=

~ S,,(a)

(5)

+A(7,±i)AI,d~(w)]~t

(10)

n=1

where S,,(w) represents the discrete Fourier transform of the nth sequence; =

~A,,(t) t1~

exp[i(o-0 —w)t]6t

(6)

where the functions d~(w), d~(w),d~(w) and d~(w)are given in the Appendix A. These functions depend only on the spherical frequency of the multiplet and on the number of points per sequence. The discrete Fourier transform of the complete signal obtained by summation over all

R. Nawab, P. Lognonné /Physics of the Earth and Planetary Interiors 84 (1994) 139—160

144

intervals can be written in the following form: N

sensitive to noise. To gain stability, we introduce a priori information on the parameter vector, following Tarantola and Valette (1982a, b). Such a priori information will be detailed in the following section. The cost function becomes

1~,+1

S(w)

=

~ =

{A1,,(w) Re[A(7~)} 1

( ~)

[

+ A~

Im A( 1’)]

+A~,(w) Re[A(I~)J .

+A~(w).Im[A(T,,)]}

3. Inverting modulation parameters To invert the modulation parameters, we define the cost function in the least-squares sense, that is, d~(t) f ~-~lI

—

S~,(t)112 2

dt

(12)

where S~(t) is the filtered theoretical signal given by 1 +~ S~~(t)=~—~ [S(o)exp(iwt) +S(-w) exp(-i~t)1

f~IId~(t)—S~~(t)112 dt

1

(11)

where the partial derivatives of the signal with respect to the modulation coefficients are introduced. In these and in the following equations, Re and Im stand for the real and the imaginary parts of a given complex number. Eq. (11) represents the signal’s theoretical Fourier spectrum, and will now be used for inversion.

1

However, spectral fitting methods are very

,,iAt C~~A (14) Here, ~ is an ad-hoc constant whose value must be determined empirically, and +

.

A= {Re[A(Tm)1, Im[A(Tm)J, Re[A(Tm)],

Im[A(Tm)I}

for m = 1 to ~ + 1, is the parameter vector that we want to retrieve. The inverted modulation function that minimizes Eq. (14) is an envelope estimate that is as close as desired to some a priori envelope while yielding the best fit to the data. The problem being linear, the minimization of the cost function (14) leads easily, using expressions (11) and (13), to the following linear problem:

[ci’

•A =D

+

where S(w) is the theoretical spectrum and is given in Eq. (11). d~~(t) and S~~(t) are, respectively, the boxcar filtered data and the boxcar filtered theoretical signal. Aw will be generally chosen in such a way that most of the energy of the simultaneously demodulated modes is concentrated within the boxcar and other modes have a small energy in that frequency band. ~T and Tf represent, respectively, the starting and the ending times of the time window on which the fitting is performed. o~is the variance of the noise within the frequency window.

(16)

The elements of the covariance matrix Cd and the data vector D are, respectively, given by 1

(13)

(15)

[C~1ap=TTI

~10S~(t)

[

aA~

aS~(t) 0A~

dt (17)

and 1

[Dk]a

=

T

~[ ~f~ L

8S~,(t) ~AZ .d~~(t)Jdt

(18)

where A~denotes the components of vector A,,. As all numerical operations, such as filtering and Fourier transformations, are practically done by using a discrete Fourier transform, we will develop the results within that framework. Introducing the complex variables i~ = exp{i(w, +

K Nawab, P. Lognonné/Physics of the Earth and Planetary Interiors 84 (1994) 139—160

o.~)6tI and j3~= exp[i(w~ w~)~t I, we can simplify the expressions of the coefficients ~ and [Dk]a by performing the time integration analytically. We finally have —

1

=

N,,—l NF—l

-~---~Re E

145

the energy within the frequency window of width Au used for the inversion in (14), it is necessary to introduce tapers in our inversion. Let us consider the family of tapers which have T f(t)=a+bcos ~r T t T +ccos2~r —~j-— —

—

‘

(21)

+dcos[3~(~)J

+ {Aa(w~ . 1

for a, 13

[Dk]a

=

1,.

[A~(w.)]*}

. .

1 ~‘

,

1

—

—

I3jj

Nr

(19)

4, and

1 NF— 1 NF— 1 ~Re ~

XIAk(w.)~d~(w.)I

as their general expression, where a, b, c and d duration the signal. family of tapers inare real of constants and This T represents the total cludes all cosine tapers currently used (such as the Hanning, Hamming and Blackman—Harris tapers). With such tapers, the tapered signal can be written in the following form: S(t) =f(t)H(t t 1)A(t) exp(hr0t) (22) —

1

—

1

(m1)

where H(t) is the Heaviside function and t~ is the starting time of the time series. In a more expanded form, the expression for a tapered signal becomes

N T

—

1—

(13..)NT

S(t) =H(t

—

t1)A(t) b

1—13.. (20)

fora,13=1,...,4. In these equations, NF represents the number of points taken from the spectrum while filtering it and NT is the total number of points in the initial time series. The asterisk stands for the complex conjugate of a given quantity. The parameter’s covariance matrix, C~,will be described in the next section. The linear inverse problem could now be solved to retrieve the parameter’s vector A. However, the interference among neighbouring multiplets would still affect the results of the inversion. In fact, Dahlen (1982) showed how this interference leads to the distortion of the resonant peaks and to a significant shift of several microhertz in their central frequencies. As these effects are greatly reduced by the use of tapers, which concentrate

X

{a exp(i~0t)

—

~exp[i(ro

+

~)tJ

_~exp[i(uo_~)tI 2~r

c +

~XP

z

+

~

exp

2

ii u

—

—

exp i o-~+

2 —

d i -~-exp~

—~

t

T 3~r

°

d —

t

2

/

+

—i--

—

T / 3~—

t

(23)

The signal appears as a sum of seven sinusoidals with complex frequencies very close to that of the carrier wave. The expression of the partial derivatives of these contributions to S~,(t)

R. Nawab, P. Lognonné/Physics of the Earth and Planetary Interiors 84 (1994) 139—160

146

with respect to the real and imaginary parts of A(1,) and A(7,) can be easily obtained from Eqs. (10), (11) and (13) by replacing the complex frequency o~by the corresponding new frequencies o~+ (n~r/T). This would not be the case for more efficient tapers, such as the n~- prolate tapers (e.g. Park et a!., 1987). It should be noted, however, that as a + b + c + d equals zero for many usual tapers (e.g. the Hamming and Hanning tapers), the summation involves derivatives of the functions d~’(w),d~(w),d~(w)and d~(w), and that the corresponding numerical codes must be expressed carefully.

computing statistical terms such as
(

~~4*(Tk’+ 1) A* Tk~)I (A*(Tk,) I A(T,,)) which involve terms such as —

(~4*(Tk)A(Tk,)),

A

Tk))

(25)

~4*(Tk)~4(Tk,))1,

KA~(Tk~A(Tk))

(26)

These terms are computed using the classical expression of the mathematical expectation such as ( A* Tk) A( Tk~))

(

=1 f +~

+~

dP(Tk,pi,p2,...)

4. A priori parameter’s covariance matrix Defining a covariance matrix helps reduce the sensitivity to noise. Here we define an a priori assumption with respect to a slight deviation from

the spherical normal modes, for which the modulation function is constant (indeed, the mean spherical attenuation is also incorporated in the

carrier wave exp(io’0t) which possesses a complex frequency). Let Ad1 be a measure of the deviation of the temporal evolution of the modulation function with respect to that for a SNRAI Earth model. We define Ad1 as the following ~ + 2 vector: Ad1

‘4(Tm)1,

XdP(Tk,,p’l,p’2,...) XA*(Tk; p1, p2, )A(Tk~p’1, p’2,...) (27) . . .

Here d P represents the probability density of the random function A(t). p1, p2,. p,, are the ran. .

dom continuous variables constraining the tern-

poral evolution of A(t). We take here as random variables the singlets’ amplitudes at the station, ak, the shift from the central angular frequency of each singlet’s angular frequency, ~Wk, and

finally, ~ak, oftheeach shiftsinglet’s from theattenuation mean attenuation coefficient coeffi-

cient.

for m = 1, N111~. In the SNRAI case, where a~is the degenerate frequency of the multiplet, and neglecting the second-order non-linearity in our problem, Ad~ will equal zero. In the laterally heterogeneous

A first level in a priori assumptions (the level chosen in this paper) is to consider that the statistical distributions of the singlets’ amplitudes and complex frequencies are given by Gaussian laws, with variances o~ for the amplitude and o~and o-,, for the real and the imaginary parts of the complex frequency, respectively. The variances o~ and of,, will then be related to the maximum a priori splitting. This hypothesis seems to model fairly well the synthetic aspherical dispersion curves, which show a concentration of the aspherical frequencies near the central frequency, and a strong decrease in the number of

rotating models, this vector will have, in general,

dispersion branches when the value of the split-

non-zero values with a statistical distribution centred on zero (i.e. ~Ad1) = 0). The covariance matrix over the model space can then be expressed in terms of the components of Ad1 by

ting increases (see, e.g. Lognonné, 1989). The statistical terms constituting the elements of the a priori covariance matrix, C~,can then be computed under the Gaussian distribution as-

=

{Re[A(Tm+i)

{

Im A( Tm + Re

~)

—

A( Tm)],

[A( Tm)], Im [A( Tm)],..., (

Re [A TN,~1+ i)}~Im [A(TN

+1)1

}

(24)

K Nawab, P. Lognonné/Physics of the Earth and Planetary Interiors 84 (1994) 139—160

147

sumption for the complex frequencies and the

rameters recovered. We will now validate the

amplitudes of the singlets. Appendix C gives the expressions of these elements. A more sophisticated level in a priori assumptions for the distribution of the singlets’ complex frequencies and amplitudes may be reached by computing synthetic distributions for these quantities using realistic 3D Earth models. Other kinds of a priori assumptions may also be considered to construct the covariance matrix C~.For example, one can take (Ad2) = (A,,~

method with a battery of tests.

5. Numerical tests

We illustrate the method by some numerical tests which show the stability of the demodulation technique with respect to noise. The tests have been performed using a synthetic seismogram of the multiplet 1S4. The frequency and the quality factor of this mode are 1.17 mHz and 271, respectively, for the PREM model (Dziewonski and Anderson, 1981). Except for the first test (Fig. 2), the synthetics have been computed by assuming a splitting width of 15%o, representative of observed splitting widths for this multiplet. The distribution of singlets within this bandwidth was random. Thus the mean frequency of the synthetic does not satisfy the diagonal sum rule (Gilbert, 1971). In theory, the fastest modulation is produced by the most split singlet, and aliasing may occur when the deviation of a singlet’s frequency with respect to the frequency used in the demodulation is too high. The first tested parameter is the

A,,+~—A,,~1~A,,) as the new measure of the deviation of the temporal evolution of the modulation function with respect to its temporal evolution in the reference model. In this case, and for the laterally heterogeneous rotating models, Ad2 will have a zero value if the multiplet can be modelled by a simple resonant peak shifted with respect to its spherical frequency cr0 (i.e. A(t) = A0 e~~rt). This a priori assumption is, for the short time approximation, equivalent to using the squared amplitude of the second derivative of the modulation (i.e. II D~‘A II 2) as an a priori cost function (see, e.g. Park,1992; Park and Maasch, 1993). At this stage, the linear system described by Eq. (16) can be solved and the modulation pa-

1 S4

Amplitude Modulation Of 1S4

—a E

h

___i; 1.14’ 1.16’ 1.18

1.2

Frequency (mHz)

Phase Modulation Of 1S4 ,.._.~

f\

IC

Is

d

‘~‘/

/ft~c51X~L

0

12

24

36

Time (h)

48

80

0

12

24

36

48

60

Time (h)

Fig. 2. Influence of the choice of sampling rate for the modulation on the inversion results. The synthetic data used are depicted in (a), which represents the power spectrum of a single singlet shifted forward by an amount of 5%c with respect to the spherical frequency of the multiplet j54, indicated by a vertical arrow on the frequency axis. (b) represents the true (continuous thick curve) and the inverted (dashed curves) modulations of amplitude. (c) represents the phase differences (dashed curves) between the inverted phases and the true phase modulation. The horizontal continuous thick curve, in this figure, indicates the zero phase difference level (i.e. the level of a perfect retrieval of the true phase). The phase is accurately fitted in all cases (maximum relative error of 1.5 x 1O~),but the relative error in the amplitude recovery increases from 1% when using four samples per modulation to 10% when using only two samples per modulation. However, in the presence of noise, inversion using four samples per modulation can only be performed by using tapers and adding some a priori information to the procedure.

K Nawab, P. Lognonné/Physics of the Earth and Planetary Interiors 84 (1994) 139—160

148

54

Phase Modulation Of

Amplitude Modulation Of i

a

i___ 1.14

1,18IIlIJl 1.18

1S4

b

1S4

1.2

a6

18

30

42

54

6

18

30

42

54

Frequency (mHz) Time (h) Time (h) Fig. 3. Illustration of the sensitivity of the demodulation procedure to noise. (a) shows four synthetic amplitude spectra of the 1S4 multiplet. The synthetics were computed assigning to each singlet a random initial amplitude and a random shift in frequency and attenuation with respect to the spherical frequency of the multiplet. Here we imposed a splitting within the bandwidth known to contain all the singlets of this multiplet (i.e. ±7.5%e around the spherical frequency of the multiplet). The arrows on the horizontal axis represents the values of the random frequencies of the singlets. Real Gaussian noise from GEOSCOPE records has been added to the spectra. The noise content in each spectrum increases from the first synthetic (at the bottom), where SNR —~oc, to the last synthetic (at the top), where SNR = 3. The second and third spectra have a signal-to-noise ratio of eight and six, respectively. (b) and (c) show respectively the retrieved modulations of amplitude and phase. No taper is used here and the inversion was carried out without any constraint (ii = 0). As the signal becomes contaminated by an increasing amount of noise, the recovered modulation amplitudes (dashed curves in (b)) begin to oscillate around the true value (thick continuous curve in (b)). The initial amplitude is well determined in all cases.

number of sampling points per modulation (spm) necessary for the interpolation of the modulation function, that is, the number of points used to sample the shortest pseudo-period contained in

single singlet, with a frequency difference of 5%~ with respect to the frequency o~ used in the demodulation. We note that the error on the phase is very small (1.5 X iO~). This is not the case, however, for the amplitude. If two samples per modulation produce a relative error of 10%, four samples per modulation reduce the relative error on the amplitude estimation to less than

the modulation function (this being produced by the most split singlet). Figs. 2(a)—2(c) show the effect of the chosen spm on the recovery of the amplitude and phase, when the synthetic is a

Amplitude Modulation Of 1S4

a

b

3S2+1S4

/)~

(\

/

a

Phase Modulation Of 1S4

C

~

—

~LA’~ 1.1

1.2

Frequency (mHz)

1.8

6

18

30

42

Time (h)

54

6

18

30

42

54

Time (h)

Fig. 4. The effect of neighbouring multiplets on the recovery of the modulation for a given multipiet. (a) represents a series of four synthetic spectra, each of them computed for two neighbouring multiplets, namely 3S2 and j54. The gap between these two multiplets originally centred at their spherical frequencies is reduced from the bottom spectrum (sf0 = 60 iiHz) to the top spectrum, where the distance separating the two multiplets is 11.7 ~Hz. On the horizontal axis are also indicated the positions of the singlets’ frequencies used in computing the synthetic spectra of 1S4. No noise has been added to the spectra. As the gap separating the two multiplets becomes smaller, the recovered amplitudes (dashed curves in (b)) and phases (dashed curves in (c)) depart from the true values (thick continuous curves in (b) and (c)). No taper or a priori information has been used. The introduction of tapers can greatly improve the results by concentrating the multiplet’s energy in narrower frequency bands.

K Nawab, P. Lognonné/Physics of the Earth and Planetary Interiors 84 (1994) 139—160

1%. On the other hand, inversion with four sam-

band of the inversion, Aw, Max( I S(w) + n(w) I) SNR Max( I n w) I) Aw Aw for w0 ~ ~ +

ples per modulation can be done only by reconditioning the inverse problem with the a priori covariance matrix, that is, with an ~ value different from zero (see Eq. (14)), or by increasing the frequency band. For the frequency band used in

=

(

(28) 2 decreases from infinity, for the spectrum at the bottom of Fig. 3(a), to three, for the spectrum at the top of this same figure. It should be noted, nevertheless, that even for the worst signal-tonoise ratio (SNR 3), the phase is very well recovered for one half of a Q cycle. Results on the amplitude are, however, much more sensitive to the increasing noise. We have found that it is important to reduce the bandwidth used in the inversion: a large bandwidth creates strong oscillations in the modulation’s amplitude for noisy signals. In the same way, secondary peaks close to the inverted mode may produce the same ef1~t.This point is illustrated in Fig. 4. We see here that the demodulation is strongly affected by tIle neighbouring mode when the two peaks begin to overlap. It is then necessary to use tapers. Fig. 5 shows the efficiency of taper in concentration the energy in the —

the inversion, demodulation using six samples per modulation was highly unstable and is not represented here. Indeed, the sensitivity to noise increases with the number of modulation parameters inverted. Only demodulation with two or four samples per modulation can be achieved, For a low sampling rate of two samples per modulation without a priori constraints, ij 0, Figs. 3(a)—3(c) illustrate the stability of the recovery of the phase and the amplitude of A(t) when increasing noise is added to the synthetic signal. No taper is applied here. The frequency bandwidth used was of 23.4 ~.tHz, corresponding to ±10%o of the spherical frequency. Noise has been taken from a GEOSCOPE record at Station SSB, and is representative of typical very long period noise. The signal-to-noise ratio, SNR, defined as the ratio of the maximum of the signal plus noise amplitude spectrum to the maximum of the noise amplitude spectrum in the frequency

—

‘~

=

Amplitude Modulation Of 1S4

a E

1S4

/

.,‘-..‘.--—--—-

,~

I

I

I

A

I <

.~

11(11 ~

E

1.16

1.18

Frequency (mHz)

6

18

/

£~

_J

—,

[

1

__ 1.2

-~

liii ~I~I

i~

1$ 11)1 1.14

11—

liii ii1t ~i~ihi

~

1

I

~ I

Phase Modulation Of 1S4

I it

\_,—.\_

149

,/

a__ 30

42

Time (h)

&4

6

18

3~O 42 Time (h)

54

Fig. 5. The necessity of introducing tapers in the inversion procedure is here illustrated. (a) shows two synthetic amplitude spectra computed for the 1S4 multiplet alone. Amplified real Gaussian noise from GEOSCOPE records has also been added to the spectra. The amount of noise is even higher than in the examples of Fig. 3(a); SNR = 2.4. The spectrum at the top of this figure has been obtained after Hanning tapering the synthetic time series, and that at the bottom has been obtained without any use of taper. The vertical arrows on the horizontal axis indicate the values of the random frequencies of the singlets used in computing the synthetic spectra. (b) and (c) show the true and the recovered modulations of amplitude and phase, respectively. The results obtained with the taper (thick dashed curves in (b) and (c)) are much closer to the true modulations of amplitude and phase (thick continuous curves in (b) and (c)) than those obtained without any use of tapers (thin dashed curves in (b) and (c)). The frequency bandwidth used in the inversion was of ±8%o around the spherical frequency of the multiplet. For this frequency band the inversion of untapered data is highly unstable. In fact, tapering concentrates the energy of the multiplets in narrower frequency bands, enhancing the signal-to-noise ratio.

R. Nawab, P. Lognonné/Physics of the Earth and Planetary Interiors 84 (1994) 139—160

150

is4

Amplitude Modulation Of 1S4

Phase Modulation Of 1~4

Ea = 0.

1.14

1.16 1.18

1.2

Frequency (mHz)

6

Ti = 10—2

—~

18

30

42

54

66

Time (h)

6

18

30

42

54

66

Time (h)

Fig. 6. Illustration of the use of a priori information in the inversion. The test is done, as before, with a synthetic spectrum of the 1S4 multiplet plus amplified Gaussian noise (a). The resulting SNR = 2.4. The time series has been Hanning tapered before Fourier transformation. The vertical arrows on the horizontal axis indicate the positions of the random singlets used in computing the synthetic spectra. (b) and (c) show recovered and true modulations of amplitude and phase, respectively, and the true modulation is represented by thick continuous curves. The constraint parameter ~ varies from zero (no constraint put on results; oscillating thin dashed curves in (b) and (c)) to infinity (horizontal thin dashed curves in (b) and (c)), where some mean amplitude and phase is retrieved. The sampling rate used here was 2 spm. The best results (represented by thick dashed curves in (b) and (c)) are obtained for ~ = 10~2.As can be seen, constraining the inversion damps rapid oscillations caused by noise and by the use of tapers at the end-points of the time series.

vicinity of the multiplet’s central frequency. This fact stabilizes the demodulation when the inversion is performed within a smaller frequency band than in Fig. 3. The noise added to the synthetic is even greater than the noisiest example in Fig. 3 (here, SNR = 2.4). The taper used is a Hanning taper and the inversion is now performed only over a frequency bandwidth of 18.7 j.rHz, which

i S4

represents ±8%o of the mode frequency, and a sampling rate of four samples per modulation is used. The result obtained with the same frequency window but without taper is shown for comparison; it clearly appears to be too sensitive to noise. It should be noted that the inversion for untapered data becomes unstable for such frequency bandwidths. This is not the case for the

Amplitude Modulation Of 1S4

Phase Modulation Of 1S4

~

_______

I:i~~I1~2:5i~l Fig. 7. Same as in Fig.1.14 Frequency 6, but1.16 this time 1.16 (mHz) four 1.2 samples per 6 modulation 18 Time 30 were 42 (h) used 54 in66the inversion. 6 _______ 18The Time 30 real 42 (h) 54 66is represented modulation by thick continuous curves in (b) and (c). As in Figs. 6(b) and 6(c), oscillating thin dashed curves correspond to the results of the unconstrained inversion (o~= 0), and the horizontal thin dashed curves are the results of the infinitely constrained inversion ~ —~ The optimum ~ parameter was found to be 10—2 (thick dashed curves in (b) and (c)). The initial amplitude is even better recovered than in the previous test when only two samples per modulation were used.

oc).

K Nawab, P. Lognonnd/Physics of the Earth and Planetary Interiors 84 (1994) 139—160

151

inversion of tapered data, for which the fit of the inverted phase with respect to the synthetic one is fairly good during the first 36 h of the Q cycle, The use of tapers induces, however, poor reso-

increasing the weight ~ on the inversion and particularly the case of an infinite weight, ~ cc, for which the inversion is limited to a constant modulation function A(t) = A0, independent of

lution of the method for the amplitude during the first hours of the signal, owing to the zero amplitude of the taper at the beginning of the time series (Fig. 5(b)). This problem is resolved if some a priori information on the modulation function is added (see Eq. (14)), as shown in Fig. 6. This a priori information penalizes the rapid variations of the modulation function and the initial zero value of the modulation induced by the taper. It also stabilizes the amplitude of the modulation function at the end of the Q cycle, when the amplitude has decreased below the noise level, Also noteworthy in this figure is the effect of

time. We have found from numerical tests that the optimum value for ~ was between 10—2 and i0~ for sampling rate of 2 spm (thick dashed curves in Figs. 6(b) and 6(c). Fig. 7 shows the same results for a sampling rate of 4 spm. In this case, the optimum value for i~is between 10—2 and 10_i (thick dashed curves in Figs. 7(b) and 7(c)). The rapid variation of the amplitude at the beginning of the signal is even better recovered than in the previous case. All these numerical tests show that the method proposed in this paper for the demodulation of tapered normal mode spectra allows a very stable

OS9

—~

Os9

OS9

Os9

0 U

a, 0. (n

______

1~

a)

0 0~

______

‘O

a, N 0

E 0

z

____

I—

,~

I

1.555 1.575 1.595 Frequency (mHz)

.•

.

1.555 1.575 1.595 Frequency (mHz)

..

.

1.555 1.575 1.595 Frequency (mHz)

I

1.555 1.575 1.595 Frequency (mHz)

Fig. 8. A selection of 0S9 multiplet resonant peaks from our data base. Each peak corresponds to a different source—receiver couple. This multiplet is one of the best observed fundamental modes. The quality and number of records available allow us to extract a large amount of information about the eigenfrequencies and the attenuation rates of the 0S9 singlets from their modulation functions. Vertical ticks on the frequency axis indicate, for comparison, the positions of the singlets in the rotating hydrostatic (RH) model of Dahlen and Sailor (1979).

152

R. Nawab, P. Lognonné/Physics of the Earth and Planetary Interiors 84 (1994) 139—160

estimation of the phase modulation, and then of the effect of splitting induced by the asphericity and the Earth’s rotation. We will now apply our method to some low angular order normal modes.

6. Application to some normal modes We show here some results for two low angular order normal modes, namely 0S9 and 1S8. A more detailed study will be presented in a future paper, where we will use the extracted modulation functions to obtain information on the eigenfrequencies of the singlets and on the aspherical structure of the Earth. Figs. 8 and 9 show a set of normal modes’ power spectra for these two modes obtained on a selection of all earthquakes of

is8

magnitude greater than M~= 7 recorded at both GEOSCOPE (Romanowicz et al., 1991) and IDA (Agnew et al., 1976) stations. All data have been carefully deglitched and the tide has been removed with a time-fitting procedure before computing spectra. The Hanning taper has been used for the spectra, and 1Q cycle time series starting 6 h after the earthquake has been used (the first 6 h of the signal being saturated on the IDA instruments). Fig. 10 shows an example of demodulation for the mode 0S9 from a spectrum recorded on the vertical component of the IDA Station SUR after the Andreanof Island earthquake of 1986, Day 127. The modulation function is expressed during one Q cycle, starting 6 h after the source time. As shown in the power spectrum (Fig. 10(a)), the mode has a central frequency

~s8

0 0

a) ci CI) a) 0

0~

a) N 0

EI.0

z

ll]III~. 1.77 1.79 1.81 Frequency (mHz)

lilIl•

,

1.77 1.79 1.81 Frequency (mHz)

IIIIIHI

,

1.77 1.79 i.B1 Frequency (mHz)

llhIlI~ 1.77 1.79 1.81 Frequency (mHz)

Fig. 9. Same as in Fig. 8, but for the anomalously split overtone ~ This multiplet is also very well observed. Peaks may have their maximum at frequencies exceeding the RH model’s predictions, as a result of lateral heterogeneities in the lower mantle, where most of the energy of this mode is concentrated.

R. Nawab, P. Lognonné/Physics of the Earth and Planetary Interiors 84 (1994) 139—160

quake:aisl86l 27 station:SUR_Z mode:0S9

153

~ 1:~/~~’\ 1~’\\ ,

1.56 1.66 Frequency (mHz)

Amplitude Modulation

8

18

30

42

54

Phase Modulation

6

18

Time (h)

30

42

54

Time (h)

Fig. 10. An example of demodulation on real data; a constrained inversion with i~i= 0.1. (a), which represents the data used in the inversion, is the power spectrum of the 0S9 multiplet from the recorded signal at Station SUR (IDA network) after the Andreanof Island earthquake of 1986, Day 127. The arrows on the frequency axis indicate the positions of the singlets in the RH model of Dahlen and Sailor (1979), as a reference. The 0S9 multiplet appears as a single broadened resonant peak shifted backward with respect to its mean spherical frequency (f0 = 1.578281 mHz). This negative shift reflects itself in the retrieved phase (thick continuous line in (c)) which shows a negative trend. The recovered amplitude (thick continuous line in (b)) has, in contrast, a positive slope, implying a negative shift in the mean attenuation coefficient (ao) and a higher Q for this multiplet. We have obtained, for this record, the following least-squares mean shifts with respect to the PREM model predictions: 6f0 = — 1.55 j.tHz as the mean shift in the frequency; 6a0 = — 0.2536 ~Hz as the mean shift in the attenuation coefficient. This leads to an estimated mean frequency of 1.576731 mHz and an estimated mean quality factor of 338.2 for this multiplet. The results of this least-squares estimation are also represented (thin dashed curves in (b) and (c)).

negatively shifted with respect to its mean spherical frequency (f0 = 1.578281 mHz) used as carrier of the mode. This appears in Fig. 10(c),

quake:sumb77231 station:CMO_Z mode:0S9

where the phase clearly shows a negative trend. It should be noted that the amplitude for this example increases, which denotes an apparent nega-

Amplitude Modulation

Phase Modulation

iI/\~1,c~1’~/ 1.56 Frequency (mHz)

i.56

8

18

30

42

Time (h)

54

8

18

30

42

Time (h)

Fig. 11. Same as in Fig. 10, but for the Sumbawa earthquake of 1977, Day 231, recorded at Station CMO (IIDA network). In this example, s~= 0.1. The positive shift in the mean frequency corresponds to a positive trend in the recovered phase (as expected), and to a positive trend in the amplitude. The least-squares estimation of the mean shift in the complex frequency of O~9 with respect to the PREM model predictions is 8f0 = + 0.9068 ~Hz for the frequency and 6a0 = — 0.1045 ~Hz for the attenuation coefficient. This leads to an estimated mean frequency of 1.579188 mHz and an estimated mean quality factor of 335.3 for this multiplet. The results of this least-squares estimation are also represented (thin dashed curves in (b) and (c)). This estimation of the modulation function appears to be a satisfactory approximation to the retrieved modulation for the O~9 multiplet.

154

K Nawab, P. Lognonné/Physics of the Earth and Planetary Interiors 84 (1994) 139—160

quake:phi176229 station:HAL_2 mode:0S9

Amplitude Modulation

~± 1.56

1.56

6

16

Frequency (mHz)

30

42

54

Phase Modulation

~:T/TJTT.’ 6

18

Time (h)

30

42

54

Time (h)

Fig. 12. Same as in Figs. 10 and 11, but for the Philippine earthquake of 1976, Day 229, recorded at the IDA Station HAL. As in previous examples, the constraint parameter ~ = 0.1. As can be seen for this record, the shift in the mean frequency of the peak with respect to the spherical frequency of 0S9 is very small, and the resulting phase (thick continuous line in (c)) and amplitude (thick continuous line in (b)) modulations are complicated and cannot be considered as a result of a simple shift in the complex frequency of ~ A more sophisticated interpretation of the inversion results in terms of the singlets’ complex frequencies is needed. (c) has been plotted at the same scale as Fig. 10(c), for comparison.

tive shift in the attenuation factor a = w/2Q, and therefore a higher Q. If the observed peak is interpreted in terms of a simple shift in the central frequency, both real and imaginary parts of this frequency must be shifted positively. In this example, as well as in many other examples,

qua ke:sumb7723 1 station:RAR_Z mode:1S8

~aj~

1.78

1.6

1.82

some oscillations are found at the beginning and end of the time series. The thin dashed curves in Figs. 10(b) and 10(c) are the results of a mean square estimation of the modulation performed by approximating the multiplet by a single broadened resonant peak which is shifted with respect

~/“~\

Amplitude Modulation

6

18

30

42

54

Phase Modulation

6

18

30

42

54

Frequency (mHz) Time (Pt) Time (h) Fig. 13. Same as in Figs. 10—12, but for the anomalously split overtone ~ after the Sumbawa earthquake of 1977, Day 231, recorded at IDA Station RAR. The power spectrum in (a) shows a widely split multiplet. Two peaks are clearly visible. The arrows on the frequency axis indicate the positions of the singlets of this multiplet in the RH model of Dahlen and Sailor (1979), for comparison. The constrained inversion was carried out with ij = 1. The resulting modulations of phase (c) and of amplitude (b) are complicated. The oscillations in the recovered amplitude are due to the beating between singlets. A simple shift in the mean complex frequency of 1S8 cannot explain the recovered modulation.

K Nawab, P. Lognonné/Physics of the Earth and Planetary Interiors 84 (1994) 139—160

to its spherical mean position. We then solve the following least-squares problem:

j T~[A( t) .[A*(t)

i5o-

0A( t )]

—

+i5o~fA*(t)} dt Minimum

(29)

to retrieve the mean shift in the complex frequency from the estimated modulation function A(t) and its estimated first derivative A(t). This mean frequency shift is a quantity equivalent to Jordan’s location parameter (Jordan, 1978). The phase modulation is approximated by 4i(t) = 4~ + 6w0t and the amplitude modulation by I A(t) I = A0 exp( t). The results obtained for the example in Fig. 10 are f = 1.57673 1 mHz, as the mean square frequency, and Q = 338.2, as the mean square quality factor. Fig. 11 shows an example, for the vertical component of the IDA Station CMO after the Sumbawa earthquake of 1977, Day 231, which, in contrast, corresponds to a positive shift in the frequency and a negative shift in the attenuation coefficient for the 0S9 multiplet. Again, the least-squares approximation to the modulation function is indicated by thin dashed curves in Figs. 11(b) and 11(c). The recovered mean square frequency is, in this case, f = 1.579188 mHz, and the estimated mean square quality factor is Q = 335.2. Generally, this leastsquares approximation is valid when the modes are well shifted with respect to the spherical frequency. For modes which do not present a strong shift, the modulation functions appear to be much more complicated. This is the case in Fig. 12, which shows demodulation results for the 0S9 multiplet recorded at Station HAL for the Philippine earthquake of 1976, Day 229. In this figure, we observe, however, an apparent negative shift in the attenuation factor. The complexity of the amplitude modulation function can even be greater for more split modes such as 1S8. Fig. 13 shows an example for this multiplet. Here, the signal consists of two apparent resonant peaks. These two peaks induce beating of the amplitude, well observed in Fig. 13(b). Again, the leastsquares estimation of the mean frequency shift is meaningless in this case, and a more laborious explanation of the retrieved modulation in terms —

~a0

155

of the singlets’ complex eigenfrequencies is needed. The quality and number of data available will allow us to retrieve a large amount of information about multiplets’ splitting and coupling an therefore about the large-scale 3D structure of the Earth.

~•

Conclusion

We have described a new method which demodulates the effects of rotation and lateral heterogeneities on the normal modes. This method allows a good retrieval of the phase, even for data with a poor signal-to-noise ratio. For data with a high signal-to-noise ratio (SNR ~ 3), it allows a simultaneous retrieval of the amplitude modulations and then of the attenuation splitting and/or coupling effects. In contrast to previous methods, it does not assume that the mode is isolated. However, modes too close in frequency must be demodulated simultaneously. This can be done by minor changes in the procedure. The extracted modulation generalize the concept of local frequency and are well adapted to studies of normal modes. They constitute a secondary observable obtained from normal mode spectra. They can be used to invert the singlets for some low angular order normal modes, as will be described in a future paper. The modulation functions can also be related to the aspherical structure and then used to constrain the inner structure of the Earth. This will be the subject of a future contribution, in preparation.

Acknowledgements We thank B. Romanowicz, J.P. Montagner and G. Roult for critical remarks and discussions. We thank also the operators of the GEOSCOPE and IDA networks for providing high-quality data. We thank F. Pollitz and J. Park for constructive reviews. This research was carried out with support from the “Action Spécialisée programmée Tomographie”. This is IPGP Contribution 1281.

156

K Nawab, P. Lognonné/Physics of the Earth and Planetary Interiors 84 (1994) 139—160

References Agnew, D., Berger, J., Bulland, R., Farrell, W. and Gilbert, F. 1976. International deployment of accelerometers: a network for very long-period seismology. EOS, 57: 180—188. Bolt, BA. and Brillinger, D.R., 1975. Estimation of uncertainties in fundamental frequencies of decaying geophysical time series (abstract). EOS, 56: 403. Bolt, BA. and Brillinger, DR., 1979. Estimation of uncertainties in eigenspectral estimates from decaying geophysical time series. Geophys. J.R. Astron. Soc., 59: 593—603. Buland, R., Berger, J. and Gilbert, F., 1979. Observations from the IDA network of attenuation and splitting during a recent earthquake. Nature, 277: 358—362. Dahlen, F.A., 1974. Interference of the lateral heterogeneity of the earth from the eigenfrequency spectrum: a linear inverse problem. Geophys. J.R. Astron. Soc., 38: 143—167. Dahlen, F.A., 1979. The spectra of unresolved split normal mode multiplets. Geophys. J.R. Astron. Soc., 58: 1—33. Dahlen, F.A., 1982. The effect of data windows on the estimation of free-oscillation parameters. Geophys. JR. Astron. Soc., 69: 537—549. Dahlen, F.A. and Sailor, R.V., 1979. Rotational and elliptical splitting of the free-oscillation of the earth. Geophys. J.R. Astron. Soc., 58: 609—623. Davis, J.P. and Henson, I.H., 1986. Validity of the great circular average approximation for inversion of normal mode measurements. Geophys. J.R. Astron. Soc., 85; 69— 92. Dziewonski, A.M., 1984. Mapping the lower mantle: determination of lateral heterogeneity in P velocity up to degree and order 6. J. Geophys. Res., 89: 5929—5952. Dziewonski, AM. and Anderson, DL., 1981. Preliminary reference Earth model. Phys. Earth Planet. Inter., 25: 297—356. Dziewonski, AM. and Gilbert, F., 1973a. Identification of normal modes using spectral stacking and stripping. EOS, 54: 374. Dziewonski, A.M. and Gilbert, F., 1973b. Observations of normal modes from 84 recordings of the Alaskan earthquake of 1964 March 28—Il. Further remarks based on new spheroidal overtone data. Geophys. J.R. Astron. Soc., 35: 401—437. Giardini, D., Li, X.-D. and Woodhouse, J.H., 1987. Three-dimensional structure of the earth from splitting in freeoscillation spectra. Nature, 325: 405—411. Giardini, D., Li, X.-D. and Woodhouse, J.H., 1988. Splitting functions of long-period normal modes of the earth. J. Geophys. Res., 93: 13716—13742. Gilbert, F., 1971. The diagonal sum rule and averaged eigenfrequencies. Geophys. J.R. Astron. Soc., 23: 119—123. Gilbert, F. and Dziewonski, AM., 1975. An application of normal mode theory to the retrieval of structural parameters and source mechanisms from seismic spectra. Philos. Trans. R. Soc., Ser. A, 278: 187—269. Hansen, R.A., 1980. Simultaneous estimation of terrestrial eigenvibrations (abstract). EOS, 61: 298.

Hansen, R.A., 1982a. Simultaneous estimation of terrestrial eigenvibrations. Geophys. J.R. Astron. Soc., 70: 155—172. Hansen, R.A., 1982b. Observational study of terrestrial eigenvibrations. Phys. Earth Planet. Inter., 28: 27—69. Hansen, R.A. and Bolt, B.A., 1980. Variations between Q values estimated from damped terrestrial eigenvibrations. J. Geophys. Res., 85: 5237—5243. Jordan, T.H., 1978. A procedure for estimating lateral variations from low-frequency eigenspectra data. Geophys. J.R. Astron. Soc., 52: 441—455. Lindberg, CR., 1986. Multiple taper spectral analysis of terrestrial free oscillations. Ph.D. Thesis, University of California, San Diego. Lognonné, P., 1989. Mod~lisationdes modes propres de vibration dans une Terre anélastique et hétérogène: théorie et applications. These de Doctorat, Université de Paris VII. Lognonné, P., 1991. Normal modes and seismograms in an anelastic rotating Earth. J. Geophys. Res., 96: 20309— 20319. Masters, G., Jordan, T.H., Silver, PG. and Gilbert, F., 1982. Aspherical earth structure from fundamental spheroidal mode data. Nature, 298: 609—613. Park, J., 1990. Observed envelopes of coupled seismic free oscillations. Geophys. Res. Lett., 17: 1489—1492. Park, J., 1992. Envelope estimation for quasi-periodic geophysical signals in noise: a multitaper approach. In: AT. Walden and P. Guttorp (Editors), Statistics in the Environmental and Earth Sciences. Edward Arnold, London, pp. 189—219. Park, J. and Maasch, K.A., 1993. Plio-Pleistocene time evolution of the 100-kyr cycle in marine paleoclimate records. J. Geophys. Res., 98: 447—461. Park, J., Lindberg, C.R. and Thomson, D.J., 1987. Multipletaper spectral analysis of terrestrial free-oscillations: Part I. Geophys. J.R. Astron. Soc., 91: 755—794. Ritzwoller, M., Masters, G. and Gilbert, F., 1986. Observations of anomalous splitting and their interpretation in terms of aspherical structure. J. Geophys. Res., 91: 10203— 10228. Ritzwoller, M., Masters, G. and Gilbert, F., 1988. Constraining aspherical structure with low-degree interaction coefficients: application to uncoupled multiplets. J. Geophys. Res., 93: 6369—6396. Romanowicz, B. and Roult, G., 1986. First-order asymptotics for the eigenfrequencies of the earth and application to the retrieval of large-scale lateral variations of structure. Geophys. J.R. Astron. Soc., 87: 209—239. Romanowicz, B., Karczewski, J.F., Cara, M., Bernard, P., Borsenberger, J., Cantin, J.M., Dole, B., Fouassier, D., Koenig, J.C., Morand, M., Pillet, R., Pyrolley, A. and Rouland, D., 1991. The GEOSCOPE program: present status and perspectives. Bull. Seismol. Soc. Am., 81: 243— 264. Roult, G. and Romanowicz, B., 1984. Very long period data from the GEOSCOPE network: preliminary results on great circle averages of fundamental and higher Rayleigh and Love modes. Bull. Seismol. Soc. Am., 74: 2221—2243.

K Nawab, P. Lognonné/Physics of the Earth and Planetary Interiors 84 (1994) 139—160 Silver, P.G. and Jordan, T.H., 1981. Fundamental spheroidal mode observations of aspherical heterogeneity. Geophys. J.R. Astron. Soc., 64: 605—634. Smith, M.F. and Masters, G., 1989. Aspherical structure constraints from free-oscillation frequency and attenuation measurements. Geophys. J.R. Astron. Soc., 94: 1953—1976. Tarantola, A. and Valette, B., 1982a. Inverse problems = Quest for information. J. Geophys., 50: 159—170. Tarantola, A. and Valette, B., 1982b. Generalized inverse problems solved using the least squares criterion. Rev. Geophys. Space Phys., 20: 219—232. Thomson, D.J., 1982. Spectrum estimation and harmonic analysis. IEEE Proc., 70: 1055—1096. Tromp, J. and Dahlen, F.A. 1990. Free oscillations of a spherical anelastic earth. Geophys. J. Int., 103: 707—723. Tukey, J.W., 1961. Discussion emphasizing the connection between analysis of variance and spectrum analysis. Technometrics, 3: 1—29. Woodhouse, J.H., 1983. The joint inversion of seismic waveforms for lateral variations in Earth structure and earthquake source parameters. Proc. Enrico Fermi Int. Schl. Phys., 85: 366—397. Woodhouse, J.H. and Dziewonski, A.M., 1984. Mapping the upper mantle: three-dimensional modeling of each structure by inversion of seismic waveforms. J. Geophys. Res., 89: 5953—5986. Woodhouse, J.H. and Girnius, T.P., 1982. Surface waves and free-oscillations in a regionalized earth model. Geophys. JR. Astron. Soc., 68: 653—673.

with d~(w)= (1

One of the advantages of our method is that the Fourier transformation can be performed analytically. Introducing the complex variable z = exp[i(o~ w)~t],we can give the following expression for the spectrum:

—

3T~+ 2T~)F~(w)

) ( ~) + ( —3 + 6’r 1) F~ ( w) + 2F~( w) +

d~((0)

=

(—

6~

(3T~ — + (6T, +

d~(w)=

(3

(T1 +

—

(1

d~(w)

=

6 T ~ F~

—

6T~)F~(w)

6T1)F~’(0))— 2F~’(w)

—

2T~+ T~)F~’(W) —

(

4 T~+ 3 T ~) F~w)

F~”(w)

3T1)F~(w) +

T~+ T~)F~(w)

(—

+

+

2T~) F~’(0))

+( —2 +

(—

2’r~+ 3T~)F~( 0))

+( —1 + 3r1)F~(w)+ F~’(w)

(t1

(

—

T

n+i

)

(A4)

n)

The functions Ff’(w), F~(w),F~(w)and F(w) can be analytically evaluated using the following recurrence relationship: F~(0))= k

~ l~F~1(w) z~T w)] —

for k

=

2

4

‘

(AS)

—

We have

N~—1 k S~(w)=k~O ~ [z] XA(k~t)

(A3)

where T1

Appendix A An analytical expression for the Fourier spectrum

157

(Al).

Replacing in Eq. (7), a~,f3,,, y~and 3~by their values (given in Eq. (9)) and then replacing A(t) by its new expression in the preceding equation,

F’(o)=

F2 ~

=

S~(w)=exp[i(o-0—w)t1~}

N—i

~

l—[z]

[z]k=

k=O N —1 k~O ~ ~ 1k

1 —z =

(A2)

p1(z) 2

Na—i

~~ +A(T~~1)d~(w) +A(I~)AI,d~(w)

~)

(~) (~)

2(z)

F~(w)=

k~O~2[z]k Nn~i

F 4(w)

=

k0

=

p

8t

~

p3(z)

(A6)

R. Nawab, P. Lognonné/Physics of the Earth and Planetary Interiors 84 (1994) 139—160

158

[z] p1(z)

N,,±i

(N~ 1)

=

(1—z)

N

“N

[z]

—

2

—

(1—z)

the signal taking z’ = new complex variable.

2

exp[i(—o-~’

w)~t1as

—

the

z +

p2(z)=

(1

[z]

—

z)

N,,+3

(A7)

2

(IV~ 1)

2

Appendix B: Partial derivatives of the spectrum Using the analytical expression of the spec-

—

4

(1 + [z]

—

z)

trum at the nth segment, eq. (10), we obtain the expressions formodulation the derivatives of the signal with respect to the function parameters.

~i~c,2+ 41V;q) (1—z)4

N,+2

(—

N,,+i,, + [z]

Forn=lwetake

2

(3JV~ -

S~(o, T0) =A(T 1) ~d~(w) which gives

2N~-1)

3+z

+

(1—z)~ —“~—~1

(A8) (1—z)4

~T.d~(w) (B1)

—

d~(w)

8A(T1) 5kw) =d~(w) 8A(T 1) —

p3(z)

=

-[z]

+ [z]

3

N,, 4 (1 +—z)~

(IV~-1)

N,,+3

2 3 (4N~ 9N (1 —z)5

+

3I~+ 3)

(B2)

For 1
we take

—

S~(w,o~)=exp[i(u 0—w)t1~]

3+9N~2+3N~-3) N,, +2 + _____________________ [z] (-6N~ (1—z)~

[z]N,+1 (4N 3 +

XiA(T~)[d~(w)+

N,,-i

d~’(w)1 [z] j

+A(T~T~[d~() +

3N~23J~,l)

I.

(1—z)

(B3)

—[z]4—3z~+3[z]2+z

(A9)

~S~(w)

~A(T~) =exp[i(o0—o.)t1~]

(1—z)~ are functions of the complex variable exp[i(o0—w)5t1, and t

—

jj

which gives

(1—z)

+

di(co)1~ N,,-i [z]

[z]Ne(_1V,3) +

T

.

aS 1(w)

and

+A(T1)

z

=

<

d~.~(w) (0)) +

[z]

N,,—i

tin

Tn+i~Tn

(AlO)

In these relationships, Nn is the number of points in the nth sequence (i.e. the number of points p for which 1, ~ t~~ T,,÷1).Equivalent expressions to (A2)—(A9) can be obtained for the dual part of

a5~~(w) =exp[i(~0—w)t1~]

1

X \d3n(0)) +

d~’(w) [z]

N,,-i

~

)

(B4)

R. Nawab, P. Lognonné/Physics of the Earth and Planetary Interiors 84 (1994) 139—160

And finally, for n

+ 1 we take

~

=

21

21

~f f

=

SN

1,,,+i(0), °~)= exp[i(o-0

—

w)tiN1,,~l

159

m=O m’=O

dP(8am)dP(~5am,)

8am

Xexp( —~amTk)exp( —~am,Tg,)

x [A(TN~,,,+l) d~’(o)

X{ffIfdP[Rem)}

+A(TN1,,,+i) .~TN..d~ie’(w)}

dP[Im(am)]

(BS) which gives 8SN,+i(w,

~)

xdP[Re(am,)]

dP[Im(am,)]a~am,

XffdP(13Wm)

dP(i~(0m~)exp(—i5wmTk)

aA(T,,) exp[i(~o — w)ti N

=

xexP(i&omrTk~)}

1,,,] [d~int(w)]

(C4)

85N 1,,,+1(~, ~)

dP() represents the probability density of the enclosed variable. To simplify the above equation, we give the following identities:

~A(~) = exp[i(o0

—

C0)ti,N~,,,][d~’mt(w)j

(B6)

ff

dP(am) dP(amr)ama~~

a,,, a,,,~

0Re(a)

JJ

+ °~n(a)=

2o~, for m = m’

=

Appendix C: Evaluating elements of the covariance matrix

dP( am) dP( am’) ama~~

a,,, am

The elements of the covariance matrix are of the form

= ~Im(a) = 0, for m ~ m’ which, once introduced in (C4), lead to

(A*(Tk,±i)_A*(Tk,) I A(Tk+i) —A(Tk)>


KA*(Tk,+l) _A*(Tk,)

—

I A(T,,)>

21 =

<.A*(Tk,)

I A(T,

1))

(Cl)

~ m=O

f

(C5)

dP(~am)exp[—6am(Tk+Tk,)}

~

I

I

To compute these quantities, we have to evaluate statistical terms such as

x

(A*(Tk)A(Tk,)),
xl a ml 2

dP[Re(am)]

dP[Im(am)]

~Re(am)~Im(am)

(A*(Tk)A(Tk,)>, (C2)

xf

dP(6Wm) exp[i~wm(Tk~—Tk)J

We write down the modulation A(t) of a given multiplet as follows:

(C6)

2!

A( t)

am exp( h50mt)

=

(C3)

=

The mathematical expectation of the dot product A*(Tk) A(Tk’) can then be expressed in the following form: KA*(Tk) ~A(T~~)>

The evaluation of the integrals, under the hypothesis of Gaussian statistical distribution for the shifts in the angular frequency and attenuation, gives
(21+ 1)[Oe(a)+O~(a)I

160

R. Nawab, P. Lognonné/Physics of theEarth and Planetary Interiors 84 (1994) 139—160

1

2

And finally, for (A*(Tk)

X exp o~(Tk + Tk~) I 2 j

Tk)2

(C7)

I

2

the

=(21+

2 \ 1)(U~e(a)+UIm(a))

]

Setting o~e(Q)+ °~Im(am) = 2o~,we obtain the final expression

x (o~ +

+

°8a (Tk~ Tk ~2 —

—

—~r~(T~, — Tk)2

xexp

=2(2l+1)o~

2

1

J

O~a(Tk+Tk~)2l 1 __________ xexP[ I

cr~(Tk+Tk,)2l I 2

~

(Tk + Tk~)2o~)

KA*(Tk) ~A(Tk,))

xexP[

~A(T~,)>, we obtain


__________ 1 —

xexpI

following expression:

]

2]

xexpl —OL(Tk’ 2

1

Tk)2 I

—

(C8)

which, by setting again

O~~e(a) + UJm(a) =

(Cli)

2o~,gives

(A*(Tk) ~A(Tk’)>

]

For KA*(Tk).A(Tk,)> we obtain the following

expression: +(Tk+ KA*(Tk) =

Tk~)2(r:a)

~A(T~,)) xexP[

(21+ 1)(~e(a) +UIm(a))

—~(Tk,

2

x [o~(T~

—

xex~[

—o~~(Tk—1k,) I

I 2 u~(Tk

+

]

Tk~~2 1 I

which, by setting as before gives

K A* ( Tk) •A( Tk~)) = 2(21 + 1) cr,,~ x [~~( Tk Tk’) x exp

{

x ex~[

(C9)

2]

—

—

Tk)21

I

Tk’) + o~(Tk+ Tk~)] 21

xexP[

—

u~(Tk

URe(a)

=

2o~,

noted that to evaluate these elements we only have to know the variances of the shift in the angular frequency and in the attenuation and the variances of the amplitude ato~ theand station, i.e. o~,~ andsinglets’ o~,respectively. o~,, may be taken as the maximum a priori splitting in

+ —

0~a(

Tk +

frequency and attenuation, respectively. The sin-

Tk~)]

glets’ amplitude variance at the station, o~,can be evaluated by measuring the amplitude of the resonant peak. This is done by performing first a

Tk~)21

I j

2

o-~( Tk + Tk~)21 2

+ Uim(a)

Io~(Tk+Tk,)2l xexpl I (C12) I 2 I i The elements of the covariance matrix are linear combinations of the quantities evaluated above (eqs. (C8), (ClO) and (C12)). It should be

I

(ClO)

demodulation without constraint (,~= 0). From this first step, the mean amplitude is computed. The second step is then a constrained inversion using the covariance matrix Ci,.