Recursive stochastic deconvolution in the estimation of earthquake source parameters: synthetic waveforms

PHYSICS OFTHE EARTH AND PLAN ETA RY I NTERJORS

~ _________

ELSEVIER

Physics of the Earth and Planetary Interiors 86 (1994) 301—327

Recursive stochastic deconvolution in the estimation of earthquake source parameters: synthetic waveforms J. Qu

~‘,

T.L. Teng

Department of Geological Sciences, University of Southern California, Los Angeles, CA 90089-0740, USA Received 21 June 1993; revision accepted 16 March 1994

Abstract In this paper, a method of the linear minimum mean-squares error (LMMSE) solution for source inversion is presented in terms of a recursive algorithm. A covariance matrix of estimation error, as well as a resolution matrix are also computed through recursion. It is shown that this recursive solution corresponds to a stationary Kalman filtering estimation for a linear dynamic system, which makes it possible to perform satisfactorily in an environment where complete knowledge of the relevant signal characteristics is not available. In a stationary environment, our recursive solution converges to the optimum Wiener solution. In a rather straightforward manner, the multichannel deconvolution problem is translated into a set of recursive expressions. The procedures have been tested using a number of synthetic data sets, including a point and a complex source, with satisfactory results. It is found that the solution is improved recursively with each addition of new data. We have found further that it is the error-covariance matrix, not the resolution matrix, that gives a measurement of the recursive performance. Since the recursive scheme of LMMSE runs in a manner based on either block-by-block or sample-by-sample operation, the memory requirement can be quite small. For problems involving sparse matrices, the recursive algorithm leads to fast and efficient computation. This method is tested by examining the Sierra Madre earthquake (M~= 5.8) of 28 June 1991, California. This event is well-recorded by the broad-band TERRAscope array. The moment tensor inversion through the presented method indicates that the solution is improved recursively when new data become available. It was found that the later arrivals on the observed seismograms have very little influence on the solution while the inclusion of new data from different stations yields substantial improvement on the mechanism to a certain point where further addition of data will not make much difference to the resulting best double-couple decomposition. However, the content of double-couple components shows a striking increase from 50 to —~ 90% with the inclusion of more data from other stations. This result demonstrates clearly the robustness of our approach since the inadequacy of the source representation and the earth model in the moment—tensor inversion may be remedied by the inclusion of more data.

1. Introduction Geophysical inversion may be viewed as an attempt to fit the response of an idealized earth model to a finite set of observations. An earth model consists of a set of relations representing a particular mathematical abstraction of an observed process. These equations in turn depend on a certain number of *

Corresponding author.

0031-9201/94/$07.00 © 1994 Elsevier Science B.V. All rights reserved SSDI 0031-9201(94)02946-9

302

J. Qu, T.L. Teng/Physics of the Earth and Planetary Interiors 86 (1994) 301 —327

model parameters which we desire to estimate from the data. The model response consists of the synthetic data produced by a particular realization of the model. The purpose of inversion is to extract the model parameters by fitting the model response to observations. In general, model response is non-linearly related to model parameters. The corresponding inverse solution can be obtained by applying a certain optimization criterion to the linearized system, provided that an adequate initial model estimate is available. Lines and Treitel (1984) discussed the linearized least-squares inversion by both Marquardt—Levenberg and singular value decompositions. The least-squares approach, which involves the use of time averages, is essentially a deterministic method in which the statistical characteristics of observations may not easily be taken into account. However, inversion is usually applied to a set of noisy data from which information about a prescribed quantity of interest is extracted (Claerbout and Muir, 1973). Noise may arise from a variety of sources. For example, the data may have been derived through noisy sensors or may represent a useful signal component that has been corrupted on its way from the source to the receiver. Hartzell (1989) also discussed the broader definition of noise in order to explain the role of inadequate Green’s functions resulting from inaccurate model assumptions. As he suggested, any phase which is present in the data but not in the Green’s functions may be treated as noise. In this sense, seismic noise can include not just the electronic and cultural noise and high-frequency scattered energy, but also distinct arrivals from structures which are not included in the assumed velocity structure. Any a priori information about the statistics of the data can help to design an optimal filter to suppress the role of noise. Tarantola (1987) explained a useful way to describe estimated uncertainties in a data set in terms of a covariance operator, which gives not only the estimated variance for each particular datum, but also the estimated correlation between errors. Basically, Tarantola’s results correspond to a MAP estimation (maximizing a posterior probability function) by the assumption of a Gaussian distribution. It is shown (Mendel, 1987) that the MAP estimation equals the mean-squares estimation under the circumstance of a Gaussian distribution. In the mean-squares estimation, the assumption of the Gaussian distribution is unnecessary. One special but very popular case is the estimator linear minimum mean-squares error (LMMSE). In seismological inverse studies, it is referred to as the stochastic inverse solution (Franklin, 1970). The standard LMMSE estimator (referred to as batch LMMSE) requires a large amount of computational time and memory space when a large number of data and model parameters are involved. In order to save computer time and memory space, iterative solutions such as ART (Gordon et a!., 1970), which does not require matrix inversion, have been introduced to seismology (Mason, 1981; McMechan, 1983; Humphreys et a!., 1984). However, these iterative solutions do not give the error and resolution estimates for the model parameters, and sometimes suffer from convergence problems. The LMMSE estimator also requires a knowledge of statistical characteristics of the data. For a real-time operation, this information is often not known completely in advance. Geophysical inversion involves the estimation of parameters of a postulated earth model from a set of observations. The corresponding inverse solution may be derived by applying a certain optimization criterion. Without a priori information about the statistics of the data and model parameters, a least-squares solution can be obtained, which involves the use of time averages. Numerical instability often arises in this approach. With stochastic properties either available or estimated in certain ways, it is found that a mean-squares solution is preferable, which involves the use of ensemble averages. By the ergodic assumption, these two solutions are the same. Given more assumptions about the statistical behavior of observation and model parameters, a further appraisal of the solution can be achieved. For example, by assuming a Gaussian distribution, a mean-squares solution also yields the MAP solution which maximizes the a posterior probability function. However, in most general cases, the assumption of the Gaussian distribution is unnecessary in the design of estimators. In this paper, a recursive LMMSE is introduced, which bears the semblance of an adaptive filter. The self-designing adaptive filter relies for its operation on a recursive algorithm, which makes it possible for

J. Qu, T.L. Teng/Physics of the Earth and Planetary Interiors 86 (1994) 301 —327

303

the filter to perform satisfactorily in an environment where complete knowledge of the relevant signal characteristics is not available. The algorithm starts from some predetermined set of initial conditions, representing complete ignorance about the environment. Yet, in a stationary environment, after successive iterations of the algorithm, it converges to the optimum Wiener solution in some statistical sense. Recursive LMMSE estimator without matrix incursion gives the same solution as batch LMMSE. The error-covariance and resolution matrices are the by-products of the recursive algorithm. The computational advantages of the recursive LMMSE are obvious in the senses that: (1) numerical stability is an inherent characteristic of the recursive algorithm; (2) it is robust in that its rate of convergence is essentially insensitive to the eigenvalue-spread problem; (3) estimates of interest are updated on either a block-by-block or a sample-by-sample basis, in which additional data are easily processed without resorting to the previous data and the sparseness of the resulting matrix can be easily treated. The properties of the recursive LMMSE are particularly useful for geophysical inversion in which we need to estimate the model parameters from a large amount of observational data. In essence, the recursive LMMSE estimator derived here is a special case of Kalman filter for a stationary environment in which the state vector of the model assumes a constant value. Zeng (1991) discussed a special case of Kalman filter under the assumption of the Gaussian distribution and applied this method to invert for the fault slip distribution by a high-frequency isochron representation. Su (1992) used the same method to invert for the site amplification factors by a large amount of coda-wave data. Their results give a good perspective of how each additional data point improves the inverse solution in which the intermediate results can be stored incrementally without going back to the whole data set. By applying the recursive LMMSE estimator, this paper investigates the problem of deconvolving a multichannel signal contaminated with noise when the source function is known. A lot of deconvolution techniques already exist and contributions come from a variety of fields. In the geophysical literature the works of Rice (1962), Treitel and Robinson (1966), Helmberger and Wiggins (1971) and Deregowski (1971) have led to practical deconvolution algorithms. Basically, all of the methods work well when the data are accurate. For noisy data, different approaches have been developed. With these developments, it has been well appreciated that resolution must be sacrificed to some degree if the output of the convolution is to be statistically stable. Moreover, their algorithms, which trade off resolution for stability, have sometimes been derived via intuition and it is not clear whether significantly better results could be attained by optimizing the trade-off in a quantitative manner. Oldenburg (1981, 1982) derived a method to deconvolve a data set when the source wavelet is known approximately, using methods of Backus and Gilbert (1970) to generate localized averages of the model which are unique except for a statistical uncertainty caused by errors in the data. However, Oldenburg’s method required a lot of operator interaction to obtain the complete statistics for the solution. Sipkin (1982, 1986a,b) used Oldenburg’s method to recover a time-dependent moment-tensor source from waveform data. In Sipkin’s analysis, he used a multichannel signal enhancement method to estimate earthquake source parameters. This method is based on the recursive solution for the Toepliz matrix, in which formal estimates of the statistics for the final results are difficult to obtain. The advantages of the current method lie in the fact that once the data have been reviewed for unusable waveforms, the inversion is almost entirely automatic and involves no further operator interaction for decision making. This is very crucial in the real-time seismology. Furthermore, the error-covariance matrix and the resolution matrix for the solution are by-products of the recursive LMMSE. In this paper, deconvolution based on the recursive LMMSE has been tested using a number of synthetic full waveform data, including a point source and a complex source with finite-time duration. The method proposed was also examined by one real event: the Sierra Madre earthquake (M~ 5.8) of 28 June 1991, California. In our numerical experiment, the Green’s functions are considered to be the multi channel input, the moment rate tensor is viewed as the convolution filter operating on the input. Using the recursive LMMSE estimator, the elements of the moment rate tensor are obtained as the =

304

J. Qu, T.L. Teng/Physics of the Earth and Planetary Interiors 86 (1994) 301 —327

output from the optimum multichannel filter. Each element of the moment rate tensor can be an independent function of time; in addition, the duration of the source time function can be constrained by limiting the length of the filter. For example, the best-fitting point source in time may be found by limiting the filter to a length of one. Linear constraints can easily be incorporated.

2. Recursive LMMSE estimator In order to estimate unknown quantities (i.e. parameters or signals) from measurements and from other a priori information, we must begin with model parameters and express them in such a way that attention is focused on the explicit relationship between the unknown parameters and the measurements. The generic linear model can be written as Z(k)=H(k)x+V(k)

(1)

In this model, measurement Z(k) is a column vector. It can have column vectors as its components. k in the equation is an index to the measurement vector Z(k). We have Z(k)=col[z(k),z(k—1),...,z(k—N+i)},z(k—i),i=0,1,...,N—i is an rn-by-i vector, where col() is an operator which forms a column vector from its arguments. Therefore, Z(k) is by an Nm-by-i vector. The parameter vector x is an n-by-i vector. The observation matrix H(k) is an Nm-by-n matrix. Measurement noise V(k) is an Nm-by-i vector with zero-mean V(k)=col[v(k),v(z—i),...,v(k—N+i)],v(k—i),i=0,1,...,N—1 is an rn-by-i vector. To proceed, we introduce several concepts. First, we introduce the fundamental theorem of estimation (Mendel, 1987, p. 110). It is stated as follows: “The estimator that minimizes the mean-squares error is xMs(k) E{x Z(k)}, E{ } is the expectation operator”, where xMs(k) is the mean-squares estimator of parameter vector x from the available data up to index k as Z(k). Second, we introduce a matrix lemma. Suppose P,~,~m are two invertible matrices and a, b are two arbitrary constants. A is a matrix which is compatible in the following matrix multiplication. Then we have =

(2) where t denotes the transpose of a matrix; —1 denotes the matrix inverse. Matrix lemma Eq. (2) can be obtained by observing that (cL4p At + bP~). (aAtP,~’A+ bP;l)PmAt ~4tp~l ~

The orthogonality principle in estimation theory is the. direct extension of the fundamental theorem of estimation (Mendel, 1987, pp. iii—1i2). It is stated as: “Suppose f[Z(k)] is any function of the data z(k), then the error in the mean-squares estimator is orthogonal to f[Z(k)] in the sense that E{[x_~Ms(k)]ft[Z(k)}} 1(k)

=

x

—

=0

xMs(k) is defined as the estimation error”.

Now, suppose that we want to derive a LMMSE estimator G, i.e. XMS(k)GZ(k) +m~,m~=E{x} Combining the above expression with a special case of function f[Z(k)] nal expression, we have

E{[x—GZ(k) —m 1}[Z(k) —mz(k)]}

0

Z(k)

—

m~(k)in the orthogo-

J. Qu, T.L. Teng/Physics of the Earth and Planetary Interiors 86 (1994) 301 —327

305

which is equal to E{(x—m~)—G[Z(k) —mz(k)]}[Z(k) —m~(k)] 0 Therefore, we obtain the LMMSE estimator G

G =P1~(k)P~(k) The LMMSE estimator can then be written as follows =m~+P~~(k)P~(k)[Z(k) —m~(k)]

XMS(k)

mz(k) =E{Z(k)},

m~=E{x}

P~~(k) =E{[x—m~][Z(k) _mz(k)]t} P~~(k) =E{[Z(k) —m~(k)J[Z(k) _mz(k)]t}

(3)

Now, we are going to derive a recursive LMMSE estimator by means of the matrix lemma Eq. (2) and the generic linear model Eq. (i). To begin, we consider the case when one additional measurement z(k + i) becomes available z(k+i)=H(k+1)x±v(k+i) When this equation is combined with the earlier generic linear model Eq. (i), we obtain a new linear model, Z(k+i)=H(k+i)x+V(k+i) (4) where

Z(k + 1)

=

col[z(k + i)

I Z(k)]

H(k+i) H(k+i)=

H(k)

and V(k+i) =col{v(k+i)IV(k)] In the above notations, col(A B) means a larger column vector constructed from two smaller column

(~)

vectors A and B. means a new matrix which consists of two submatrices C and D. Putting the generic linear model Eq. (1) into the LMMSE estimator Eq. (3), we get

t(k)

and

P~~(k) =E{[x

—

m1][H(k)(x _mx)]t} =P~H

P~~(k) =E{[H(k)(x-m~)

+

V(k)} [H(k)(x-m

t(k)+R(k) 1)

where

P~=E{[x_mx][x_mx]t},

+

V(k)]t} =H(k)P~H

R(k) =E{V(k)V(k)t}

Therefore, the LMMSE estimator can be written as xMs(k) =mx+PxHt(k)[H(k)PxHt(k) When z(k

+

xMs(k

+R(k)J’[Z(k)

—m~(k)]

1) is available, the new LMMSE estimator becomes

+

i)

=

m~+P~Ht(k+ i)[H(k

+

i)P~Ht(k+ 1)

+

R(k + i)}

1

x [z(k

+

1)

—

m~(k+ i)]

306

J. Qu, TL. Teng/Physics of the Earth and Planetary Interiors 86 (1994) 301—327

To proceed further, we must assume that R(k) is a block-diagonal matrix, i.e. R(k

+

i)

diag[r(k

=

+

i) I R(k)}

(5)

By the matrix lemma Eq. (2) and subdivision structures of vector and matrix discussed above in Eqs. (4) and (5), we get ~(k + 1) =m~+[Ht(k + i)r’(k + i)H(k + i) +Ht(k)R1(k)H(k) +P’J~ x{Ht(k+ 1)r1(k+ i)[z(k+ i) —m

t(k)R1(k)[Z(k) 5(k+ 1)] +H where subscript MS is omitted for simplicity. Defining a new quantity P(k) as P’(k) =Ht(k)R1(k)H(k)

+m~(k)]}

+1’;’

We have P’(k

+

1) =P’(k) +Ht(k

Defining another matrix K(k

+

+

i)r’(k

+

i)H(k

+

i)

1) as

K(k + 1) =P(k + i)Ht(k + i)r’(k

+

i)

We obtain a recursive expression for the LMMSE estimator +

i) =~(k) +K(k + 1)[z(k + i) —H(k + i)~(k)]

K(k+i)=P(k+i)Ht(k+i)r’(k+i) P1(k + i) =P’(k) +Ht(k + i)r1(k

+

i)H(k + i)

(6)

Using the matrix lemma Eq. (2), above recursive expression can be rewritten as

K(k

+

i) =~(k) +K(k

+

+

1) =P(k)Ht(k

+

1){z(k

+

1) —H(k + i)2(k)J

1) {H(k + i)P(k)Ht(k + 1) + r(k +

P(k+i)=[I—K(k+i)H(k+1)JP(k)

1)1

1

(7)

The recursive equations, Eqs. (7), are equivalent to a stationary Kalman filtering estimate for a linear dynamic system. Mendel (1987, pp. 170—181) proves that the stationary Kalman filter is an infinite-length digital Wiener filter and the truncated stationary Kalman filter is a finite-impulse-response digital Wiener filter. Both filters minimize error variances. Therefore, our derived Eqs. (7) indeed give the optimum Wiener solution in a stationary environment. The recursive LMMSE estimator Eq. (7) is computationally very efficient since it only requires to invert a smaller matrix. For a special case where r(k + 1) is a scalar, no matrix inversion is needed. The sparseness of the matrix can also be taken into account in this recursiOn. The first form of the recursive LMMSE estimator Eq. (6) is also useful in deriving the initial conditions. Mendel (1987) discussed a possible choice of initial conditions for the recursive least-squares estimator. In an online application, we also apply some of his conclusions to the recursive LMMSE estimator. For example, the initial conditions can be given as 1 P1(0) ~(0)

=

=

—~I+Ht(0)r’(0)H(0)

P(0)

+

Ht(0)r1(0)z(0)]

(8)

where a is a very large number; r is a vector with very small values. Eqs. (8) are derived by Mendel (1973, pp. 101—106). When the initial values eqs. (8) are used in the recursive equations (7), then the

J. Qu, T.L. Teng/Physics of the Earth and Planetary Interiors 86 (1994) 301 —327

307

resulting estimates for model parameters are the very same ones that are obtained from the batch LMMSE estimator which makes use of the direct matrix inversion. It can be shown that P(k) is the covariance matrix of estimation error for the recursive LMMSE estimator (9) P(k) =E{[x—x~(k)}[x—~(k)]t} Using matrix lemma Eq. (2) and subdivision of vector and matrix in Eqs. (4) and (5), we can also derive an expression for the model resolution matrix M(k) M(k) =I—P(k)P;1

(10)

From this expression, we can easily obtain a recursive expression for M(k)

(11)

M(k+ 1) =M(k) +K(k+ i)H(k+ i)[I—M(k)]

The derivation of Eqs. (9), (10) and (11) is given in Appendix i. To implement the recursive LMMSE estimator, we need to know the measurement noise variance. Usually, it is assumed to be a constant matrix. The different choice of noise variance may affect the performance of recursion as shown later. We think it is a better idea not just to replace the noise variance with the data variance. It should also involve the uncertainty of the model parameters. Through matrix block manipulation, we can obtain a pseudo-recursive expression for r(k + 1) r(k+1)=r(k) + {[z(k + 1)— H(k

+l)2(k)]t[z(k +1)— H(k

+ 1)~(k)]—tr[H(k + 1)P(k)Ht(k + 1)]— r(k)} (k +1)dim[z(k+1)]

11

where dim() is the dimension of a vector; tr() is the trace of a matrix (see Appendix 2).

3. Multichannel deconvolution based on recursive LMMSE estimator The general case of multichannel input and multichannel output can be reduced to the following form ~gk(t)

u(t)

*

sk(t)

Now defining the column vectors g(t) =col[g 1(t), g2(t),...,g~(t)} and s(t) =col[s1(t),

s2(t),...,s~(t)]

we can rewrite the convolution model as u(t) =gt(t) ®s(t) where 0 denotes the summation over the convolution formed by each corresponding component of column vectors g(t) and s(t). In actual design, g(t) and s(t) have finite lengths. Now, supposing that the effective discrete form u(t), and g(t) have a length of N indexed from k 0, 1,. N 1, and the duration of s(t) lasts over M discrete time intervals, we get the discrete form of the convolution model N—i t(k i)s(i) (13) u(k) g . .,

—

=

~

—

308

J. Qu, T.L. Teng /Physics of the Earth and Planetary Interiors 86 (1994) 301 —327

This convolution equation may then be written in a matrix form as u(N—1) u(N— 2)

gt(N—1) gt(N— i)

gt(N—2) gt(N— 3)

u(M) u(M— 1) u(M—2)

gt(M) gt(M— 1) gt(M—2)

gt(M— i) gt(M—2) gt(M—3)

gt(i) gt(0)

gt(o) 0

=

u(i) u(0)

g~(N—M+1) gt(N—M)

gt(N—M) gt(N—M— 1)

.

gt(2)

gt(i)

.

.

gt(i)

gt(o)

.

.

g~(0)

0

.

: .

o .

0

0

s (0) s( i)

s(M-2) s(M-i) The convolution model has been translated into the generic linear model, i.e. Z(N— 1) =H(N— 1)x+ V(N— 1) where the actual measurement is contaminated by noise vector V(N 1). Observation matrix H(N 1) is defined by the above g matrix. Since we have obtained the generic linear model, the established recursive LMMSE estimator can be applied directly. This process is termed multichannel recursive stochastic deconvolution. —

—

4. Numerical experiment For the purpose of testing our recursive stochastic deconvolution, synthetic seismograms and Green’s functions were computed using the generalized reflection/ transmission method (Kennett and Kerry, 1979; Kennett, 1980) and flat-layered model for the San Bernardino mountains (Table 1). Both the seismograms and Green’s functions are of full responses in a frequency band from 0.03 to 2 Hz. The sampling rate is 0.2 s. Five hundred and twelve data points are computed for both the seismograms and the Green’s functions. We have selected seven stations, for each of which three-component synthetic seismograms have been computed and will be used as the input data. So, we have a total of 7 x 3 x 512 10752 data points. In our numerical experiment, the number of unknown model parameters ranges from six to 360. This is clearly an overdetermined problem. The experiment was conducted in the following three ways. First, we perform a single-channel deconvolution where impulse response is known for each station. This gives the scalar source time function with finite duration. Secondly, the six =

J. Qu, T.L. Teng/Physics of the Earth and Planetary Interiors 86 (1994) 301 —327

309

Table 1 Velocity model of the San Bernardino mountains (Magistrate, 1990) Thickness (km) 4.00 4.00 6.00 6.00 6.00 6.00 15.00 50.00 60.00

Density 3) (gcm 2.57 2.60 2.64 2.68 2.74 2.92 3.30 3.37 3.40 3.44

P velocity (kms1)

S velocity (kms1)

5.67

3.29 3.46 3.52 3.64 3.75 3.93 4.47 4.36 4.29 4.50

6.00 6.10 6.30 6.50 6.81 7.74

7.55 7.43 7.80

Synthetics by Flat-Layered Model with Finite Source Duration

5 a) E

.0

z

______

______

0

3—

_______

~

1

0.0

_

20.0

40.0 60.0 Time (sec)

80.0

100.0

Fig. 1. Synthetic seismograms computed using Kennett’s method for flat-layered models are obtained by convolving the impulse response with a finite source time function. Seven stations with three components are considered. Z, R, T in the figure denote the vertical, radial and transverse components at each station, respectively.

J. Qu, TL. Teng/Physics of the Earth and Planetary Interiors 86 (1994) 301 —327

310

moment rate tensors are regarded as the multichannel input. A multichannel recursive stochastic deconvolution is applied to noisy data to obtain the optimal filtered results for a best-fitting point source in time. Thirdly, each element of the moment rate tensor is treated as an independent time function. The duration of the source time function can be constrained by limiting the length of the filter. Then a multichannel recursive stochastic deconvolution is applied. In all the three tests, we demonstrate a very useful feature of the recursion by allowing a progressive addition of available data points. In most cases, the comparison is carried out on the inversion results based on data (1) from two stations and (2) from seven stations. From the inherent character of the recursion process, the computation can be terminated at any step and the results can be stored until additional data become available for further improvement of the estimation. These intermediate results give us a good perspective of how each additional observation improves our inverse solution. Because results of our inversion can be incrementally stored, we can always update our solution when new data points become available with no need of repeating the

Synthetics by Flat-LayeredModel with Finite Source Duration

~ 5 a) E z .0

0.

Cl)

2

T

1

R

z 0.0

20.0

40.0 60.0 Time (sec)

80.0

100.0

Fig. 2. Synthetic data are constructed by adding uniformly random noise to results in Fig. 1. The amplitudes for noise are allowed to reach 20% of the amplitudes of synthetic seismograms.

J. Qu, T.L. Teng/Physics of the Earth and Planetary Interiors 86 (1994) 301 —327

0.2

I

I

I

311

I

7 0.14 I 0.12 I 0.1 /•

0.08 ,0’

0.06 I

0.04 0.02~

.~

I

2

I

4

6 Time (sec)

8

10

12

Fig. 3. The deconvolution result for source time function by perfect data is shown in the figure as circles. The connecting dot-dashed line is the source time function. The correct solutions are obtained exactly. In our computation, the sampling rate is taken as 0.2 sec.

computation over the whole data set. The covariance matrix of estimation error and resolution matrix are the by-products of recursion. As shown later, the resolution matrix is not a good measurement of recursive performance. In all test cases, the resolution matrix is approaching the identity matrix, which makes it difficult to tell the improvement of estimation by the addition of more data points. Meanwhile, the covariance matrix of estimation error is indeed a very good measurement of recursive performance. The improvement in estimation is clearly seen for the addition of more data points. In Fig. 1, synthetic seismograms computed by using Kennett’s method for flat-layered models are obtained by convolving the impulse response with a~ finite source time function. Seven three-component stations are considered. Z, R and T in the figure denote the vertical, radial and transverse components at each station, respectively. In Fig. 2, synthetic data are contaminated by adding uniformly random noise. The amplitudes for noises are allowed to reach 20% of the maximum amplitudes of synthetic seismograms. In Fig. 3, the deconvolution result for the source time function by perfect data is shown in the figure as circles. The connecting dot-dashed line is the given known source time function. The solution is practically perfect. In Fig. 4, the deconvolution result for the source time function by the noise contaminated data is shown as a solid line. Connected circles give the exact solution. In this example, only data for the first two stations are utilized in the deconvolution. The influence of noisy data is reflected by large oscillation of estimated parameters around the exact solution. A lower order polynomial least-squares regression gives a good approximation to the exact solution. In Fig. 5, the resolution matrix of the recursive LMMSE estimator is shown, which corresponds to the results shown in Fig. 4. The diagonal shape of the resolution matrix in Fig. 5 shows that its values are very close to one, which is the resolution of standard least-squares estimate. So just for data from two stations, the resolution matrix has reached its saturation

J. Qu, T.L. Teng/Physics of the Earth and Planetary Interiors 86 (1994) 301 —327

312

0.4

I

I

I

I

2

4

6

8

0.3

0.2

0.1

0

-0.1

10

12

Time (sec) Fig. 4. Deconvolution result for source time function by contaminated data is shown as a solid line. Connected circles are the exact solution. In this example, only data from the first two stations are utilized in deconvolution.

0.8 0.6 0.4 0.2

6g

~ndex2020X~Index

Fig. 5. The resolution matrix of the recursive LMMSE estimator is shown, which corresponds to the results as shown in Fig. 4. The diagonal shape of the map demonstrates that the resolution matrix is very close to the identity matrix, which is the resolution of the standard least-squares estimate.

J. Qu, T.L. Teng/Physics of the Earth and Planetary Interiors 86 (1994) 301 —327

313

xlO4 42

0 -2.

~:x202Oxfldex~°

~:

Fig. 6. The three-dimensional mesh map for the covariance matrix of estimation error corresponds to the deconvolution result in Fig. 4. There is a scalar difference in the map.

value and cannot provide a further insight into the performance of recursion. In Fig. 6, a three-dimensional mesh map is shown for covariance matrix of estimation error which corresponds to the deconvolution result in Fig. 4. There is a scalar difference in the map. Fig. 7 gives a similar result to those in Fig. 4 0.251 I

I

I

I

I

2

4

6

8

10

0.2

0.15

0.1

0.05

0

~O.O5~

12

Time (sec) Fig. 7. Similar result to those in Fig. 4 except that data from all seven stations join the deconvolution. The smooth solid line is the exact solution.

314

.1. Qu, T.L. Teng /Physics of the Earth and Planetary Interiors 86 (1994) 301 —327

x 10

20

00

20

X-index

40

Fig. 8. Same configuration as in fig. 6 except that data from seven stations join the deconvolution. The improvement of the estimation error over Fig. 6 is clearly viewed.

1.25

~

2

4

6 Time (sec)

8

10

12

Fig. 9. This map shows the influence of a different choice of variance for observation noise. Here we do not use our derived recursive estimator for noise variance. A constant noise variance is used here. The deconvolution is for data from the first two stations. The smooth solid curve is the exact solution.

J. Qu, TL. Teng/Physics of the Earth and Planetary Interiors 86 (1994) 301 —327

0.4

I

I

I

I

315

I

0.3 0.2 0.1

0 -0.1

I

-0.2 ~0.30

12 Time (sec)

Fig. 10. Similar to Fig. 9 but for data from five stations. The improvement is not as good as that of Fig. 7 over Fig. 4. This map shows that our recursive estimator for noise variance indeed yields a better estimation of model parameters.

except that data from seven stations join the deconvolution. The smooth solid line is the exact solution. Fig. 8 is similar to Fig. 6 except that data from all seven stations join the deconvolution. The improvement of the estimation error over Fig. 6 is obvious. The covariance matrix of estimation error offers a good performance indicator of recursion. Fig. 9 shows the influence of different choices of variance for observation noise. Here we do not use our derived recursive estimator for the noise variance. Instead, a constant noise variance is given. The deconvolution is for data from the first two stations. The smooth solid curve is the exact solution. Fig. 10 is similar to Fig. 9 but for data from five stations. The improvement is not as good as that of Fig. 7 over Fig. 4. So we show here that our recursive estimator for noise variance indeed yields a better estimation of model parameters. Fig. 11 shows examples of Green’s functions at one station computed by using Kennett’s method for flat-layered models. g1 to g6 correspond to the six different combination of moment rate tensor. Fig. 12 shows synthetic seismograms of a given moment rate tensor assuming a step source time function. In Fig. i3, synthetic data are constructed in the say way as in Fig. 2. Fig. 14 shows the deconvolved best-fitting point source in time. Solid line denotes a double-couple focal mechanism for a given source. Dashed line denotes the deconvolved focal mechanism which is expressed in terms of both equal area projection and its best double-couple counterpart. In this figure, only data from the first two stations join the deconvolution. Fig. 15 has the same configuration as in Fig. 14 but for data from all seven stations. Improvement of the deconvolved focal mechanism is also obvious. However, it should be noted that while the systematic improvement of estimation can be achieved by including more available data points, whether it is worthwhile to do so in certain circumstances is another issue. For example, as shown in Fig. 14 and Fig. 15, although the improvement is clear, it is not substantial enough to be important. This difference may be important when we deal with the real data.

316

J. Qu, T.L. Teng/Physics of the Earth and Planetary Interiors 86 (1994) 301 —327

Examples of Green’s Functions

:4

__

0.0

20.0

40.0

60.0

80.0

100.0

Fig. ii. Examples of Green’s functions at one station computed by using Kennett’s method for flat-layered models. g

1 to g6

correspond to the six different combinations of the moment rate tensor.

Fig. i6 shows the deconvolved results when each element of the moment rate tensor is treated as an independent function of time. In this figure, data from the first two stations are included. Smoothed solid curves are the exact solutions of finite-duration moment rate tensor. Fig. 17 is the same as Figure i6 except that data from all seven stations are included. Again, we see how more data can quantitatively improve the estimation model parameters in a recursive manner. From the above numerical experiments, we conclude that our recursive stochastic deconvolution is very flexible and that it may be used either routinely or as a research tool for studying particular earthquakes in details. It requires very little operator interaction while the complete statistics for the solution can be obtained.

J. Qu, T.L, Teng/Physics of the Earth and Planetary Interiors 86 (1994) 301 —327

317

5. Results for the Sierra Madre earthquake of 28 June 1991, California Given the results of the above synthetic tests which illustrate the systematic contaminated circumstances, we can now apply the presented method to the inversion of moment tensors for a real event: the Sierra Madre earthquake of 28 June 1991, California. Fig. 18 displays the locations of the existing TERRAscope array stations. This figure was modified from Fig. 1 in Heimberger et a!. (1993). The coordinates of the TERRAscope stations and their geometric configuration with respect to the Sierra Madre earthquake are given in Table 2. The one-dimensional crustal structure over the study area is given in Table 3. Tables 2 and 3 are also given in the paper by Helmberger et al. (1993). Since only the earthquake source is to be studied, some constraints may be incorporated into the system of equations. The constraint that the source is due to purely deviatoric stress is imposed here before our recursive multichannel deconvolution is applied to obtain the moment-tensor solution. To determine the possible

Synthetics by Flat-LayeredModel ~

I

I’

I~

I

-___

0.0

20.0

40.0 60.0 Time (sec)

80.0

100.0

Fig. 12. Synthetic seismograms of given moment rate tensor assuming step source time function.

318

J. Qu, T.L. Teng/Physics of theEarth and Planetary Interiors 86 (1994) 301 —327 Synthetic Data with Uniform Random Noise Added I

7

—~

I

I

-~

5

__________

I_ 0.0

20.0

40.0 60.0 Time (sec)

80.0

100.0

Fig. 13. Synthetic data constructed in the same way as in Fig. 2.

influence of the inclusion of additional data, we gradually increase the number of the TERRAscope observations from one to four. At the fixed number of observations, the different lengths of data windows are also considered to assess the sensitivity of the inversion results to the complex later seismic arrivals. The modeling indicates that our data-adaptive optimal filters are consistent in tracking a reasonable source mechanism instead of being affected by the complex later arrivals. However, the addition of new data from other stations may, on the contrary, change the source geometry significantly. The moment-tensor solutions obtained from one, three and four stations are given in Table 4. Fig. i9(a), (b) and (c) show the corresponding source mechanisms with respect to the moment-tensor solutions involving the data from one, three and four stations. The source geometry of Fig. i9(a) differs from that of Fig. 19(b) and (c), while Fig. 19(b) and (c) have very similar patterns if the configuration of the source planes I and II can be interchanged. Meanwhile, the content of double-couple components increases steadily with the inclusion of new data from other stations from 51.74 to 66.80%, finally to 91.16% (see

J. Qu, T.L. Teng/Physics of the Earth and Planetary Interiors 86 (1994) 301 —327 *

Focal

Mechanism

319

*

for upper semisphere by equal area projection N

—

I

—

Ic

I

W~f,7

\

\‘

iI

V

S Fig. 14. Deconvolved best-fitting point source in time. Solid line gives the double-couple focal mechanism. The dashed line denotes the deconvolved focal mechanism which is expressed in terms of both equal area projection and its best double-couple counterpart. In this figure, only data from the first two stations join the deconvolution.

Table 4). We interpret these results as an indication of the robustness of our recursive stochastic deconvolution since the inclusion of more useful data will yield a more reliable earthquake source mechanism. This also illustrates that inadequate data can cause the distortion of source mechanism by introducing the non-double-couple components during.

6. Discussion and conclusion In this paper, a method is presented to express the LMMSE solution in terms of a recursive algorithm. A covariance matrix of estimation error as well as a resolution matrix are also computed through recursion. It is shown that our recursive solution corresponds to a stationary Kalman filtering estimation for a linear dynamic system, which makes it possible to perform satisfactorily in an environment where complete knowledge of the relevant signal characteristics is not available. The algorithm starts from some predetermined set of initial conditions, representing an ignorance about the environment. It is numerically stable owing to its inherent characteristics and robust in that its rate of convergence is essentially insensitive to the eigenvalue-spread problem. In a stationary environment, it converges to the optimum Wiener solution. In a rather straightforward manner, the multichannel deconvolution problem is translated into recursive expressions. This provides a direct way to study source parameters for particular earthquakes. The procedures have been tested using a number of synthetic data sets, including point and

320

.1. Qu, T.L. Teng/Physics of the Earth and Planetary Interiors 86 (1994) 301—327 *

Focal

Mechanism

*

for upper semisphere by equal area projection N

W

E

T

S Fig. 15. Same configuration as in Fig. 14 but for data from all seven stations. Improvement of the deconvolved focal mechanism is also obvious.

complex sources with finite durations, with satisfactory results. It is found that solution is improved systematically with the addition of more data. It is the error-covariance matrix not the resolution matrix that gives a measurement of the recursive performance. Since the recursive style of LMMSE runs in a manner based on either block-by-block or sample-by-sample operation, the memory requirement can be quite small. The sparseness of matrix can also be made full use of as in the iterative method. It leads to a fast and efficient inversion in certain situations. The method presented is also tested by examining a real event: the Sierra Madre earthquake (M~ 5.8) of 28 June 1991, California. This event is well-recorded by the broad-band TERRAscope array. The moment-tensor inversion through our method indicates that the solution is improved recursively when more data are available. The later arrivals on the observed seismograms have very little influence on the solution while the inclusion of new data from different stations yields substantial improvement on the mechanism to a certain point where further addition of data will not make much difference to the resulting source geometry. Meanwhile, the content of double-couple components shows a striking trend of increase with the inclusion of more data from different stations. This is an indication of the robustness of our recursive stochastic deconvolution since the inadequacy of the source representation and the earth model in the moment-tensor inversion may be remedied by the inclusion of more useful data to yield a more reliable earthquake source mechanism. The multichannel recursive stochastic deconvolution derived here inherits all advantages of the multichannel signal-enhancement (MSE) method discussed by Sipkin (1982, 1986a,b). Once the data =

J. Qu, TL. Teng/Physics of the Earth and Planetary Interiors 86 (1994) 301 —327

~ ~

15

~ 0

~

Time (see)

15

Time (sea)

15

15

Time (see)

:~+

321

Time (see)

:~.. I

15

Time (sec)

Time (see)

Fig. 16. Deconvolved results when each element of the moment rate tensor is treated as an independent function of time. In this figure, data from the first two stations are included. Smoothed solid curves are the exact solutions of the finite-duration moment rate tensor.

~

15

Time (see) 0.1

15

lime (see)

0.04 15

15

Time (see)

Time (see) 0.04

I :~ ~ ~

~

Time (see)

Time (see)

Fig. 17. Same as in Fig. 17 except that data from all seven stations are included.

j5

322

J. Qu, T.L. Teng/Physics of the Earth and Planetary Interiors 86 (1994) 301 —327

Sierra Madre Event: 28 June 1991

I

36

V ~&

-120

I

-

-118

-116

West Longitude Fig. 18. Map of southern California shows the locations of the TERRAscope stations and the Sierra Madre event indicated by the star. Note that four of the stations are essentially at the same range of 160 km: GSC, ISA, PFO, and SBC. This map is modified from Helmberger et al. (1993).

Table 2 Velocity model in the study area (Helmberger et al., 1993) Thickness (km)

Density 3) (gcm 2.40 2.60 2.80 3.00 3.40

5.50 9.50 19.50 5.00

P velocity (kmst)

S velocity (kms1)

5.50 6.30 6.70 7.80 7.85

3.00 3.60 3.80 4.30 4.40

Table 3 TERRAscope station locations (Heimberger et at., 1993) Station

Latitude (°N)

Longitude (°W)

Distance (km)

Azimuth (degree)

USC ISA PAS PFO SBC

35.300 35.643 34.148 33.609 34.442

116.810 118.480 118.172 116.455 119.713

158.2 159.6 20.6 159.6 159.4

44 344 232 117 277

J. Qu, TL. Teng/Physics of the Earth and Planetary Interiors 86 (1994) 301 —327

323

Table 4 Moment-tensor solutions from the different data sets TERRAscope stations Moment tensors M,.~ M, 1~ ~ ~ M5~ M~0 Principal axes T I P

(USC, ISA, PFO)

(GSC, ISA, PFO, SBC)

0.53

—1.24

—0.44

—0.20 1.73 —2.28 1.75 —2.26

—0.19 5.80 —1.02 1.16 —4.56

—0.01 4.32 —0.97 —0.43 —3.88

2.97° 0.94° —3.91°

5.94° —0.99° —4.95°

4.34° —0.19° —4.15°

27.57° 11.14° 59.90°

6.55° 14.57° 73.97°

3.10° 14.51° 75.15°

123.96° 219.86° 329.72°

92.43° 184.14° 338.90°

269.21° 178.40° 10.98°

3.44 51.74

5.45 66.80

4.25 91.16

a

Plunge T I P Azimuth T I P Best double-couple Scalar moment a Percentage double-couple Nodalplane 1/2 Strike Dip Rake a Scale = 1023 dyn-cm.

W

(GSC)

a

43 ‘/187° 73 ‘/20° —78°/—124°

15 ‘/166° 53 ‘/40° —72°/—113°

165 ‘/14° 50 ‘/44° —109°/—69°

for upper semisphere

for upper semisphere

for upper semisphere

by equal area projection

by equel area projecfion

by equal area projection

-~~‘~)E

(a) data from GSC

W

(b) datafrom GSC,ISA and FF0

E

W

E

Cc) data from GSC,ISA.FFOand SBC

Fig. 19. Deconvolved best-fitting point source in time. The deconvolved focal mechanism and its best double-couple counterpart are shown. (a) Only data from GSC station join the deconvolution. (b) Data from USC, ISA and PFO stations join the deconvolution. (c) Data from GSC, ISA, PFO and SBC stations join the deconvolution.

324

J. Qu, T.L. Teng/Physics of the Earth and Planetary Interiors 86 (1994) 301—327

have been reviewed for unusable waveforms, it is almost entirely automatic and involves no further decision making. Meanwhile, the major disadvantage of the MSE method (i.e. very difficult to obtain formal estimates of the statistics for the final estimation), is overcome in our method, since the error-covariance matrix as well as the resolution matrix are the by-products of the recursive LMMSE. The multichannel recursive stochastic deconvolution bears the semblance of an adaptive filter, which reduces the interaction of the operator to a very low level. This feature forms a sharp contrast with the multichannel vector-deconvolution method proposed by Oldenburg (1981,1982) which requires intensive machine-operator interaction. Therefore, we think that the recursive stochastic deconvolution offers a more straightforward, more flexible and more complete analysis when applied to the study of earthquake sources in detail.

Appendix 1: The derivation of Eqs. (9), (10) and (ii) In this appendix, we derive Eqs. (9), (10) and (11) in context. The LMMSE estimator can be written as (see Eq. 3) (A1.1) 1(k) =mx+PxHt(k){H(k)PxHt(k) +R(k)] 1[Z(k) -m~(k)] where subscript MS is omitted for simplicity, and Px=E[(x_mx)(x_mx)tI,

R(k) =E[V(k)V(k)t}

With Eq. (A1.1), we have ‘~k=E{[x—1(k)}

[x_1(k)]t}

=E({(x_m

1[Z(k) -mz(k)]} 1) _P~Ht(k)[H(k)P~Ht(k) +R(k)} x{(x-m~) _P~Ht(k)[H(k)P 1[Z(k) -m~(k)]}) 1Ht(k) +R(k)] =P~_P~Ht(k)[H(k)P~Ht(k)+R(k)} 1E{[Z(k) -mz(k)] [Z(k) mz(k)}t}

-

x[ii(k)P~Ht(k)+R(k)]lH(k)P~

(A1.2)

By the generic linear model Eq. (1), we obtain E{[Z(k) —mz(k)] [Z(k) _mz(k)]t} =H(k)P~Ht(k) +R(k)

(A1.3)

Putting Eq. (A1.3) into Eq. (A1.2), we have PkXxI’(k)[H(’c)”X”(1o)

+R(k)1H(k)P~.

(A1.4)

By the matrix lemma Eq. (2), Eq. (A1.3) can be written as ~k~’x =

[p_1+Ht(k)R_1(k)H(k)]l[J+Ht(k)R_1(k)H(k)p_Ht(k)R_1(k)jj(k)p]

=

This is the definition of P 1(k). So we get the equation P(k) =E{[x—1(k)][x—1(k)]t}

(9)

J. Qu, T.L. Teng /Physics of the Earth and Planetary Interiors 86 (1994) 301 —327

325

According to the definition of the model resolution matrix, we have M(k) =P~Ht(k)[H(k)P~Ht(k) +R(k)]

‘H(k)

(A1.5)

By the matrix lemma Eq. (2), Eq. (A1.5) is equal to M(k) =P~Ht(k)[H(k)P~Ht(k) +R(k)]1Ht(k)R1(k)H(k)

(A1.6)

By the subdivision structures of the vectors and matrices in Eqs. (4) and (5), we have r1(k + 1)

0

H(k + 1)

0

R1(k)

H(k)

M(k+ 1) =P(k+ 1)[Ht(k+ 1), Ht(k)]

=P(k+ 1)[Ht(k+ 1)r1(k+ 1)H(k+ 1) +Ht(k)R1(k)H(k)] By the definition of P1(k), we have M(k+ 1) =P(k+ 1)[P_1(k+ 1) —P’(k) +P1(k)

—~1I

Therefore, M(k+1)=I—P(k+1)P1

(A1.7)

i.e.

M(k)=I—P(k)P’

(10)

By Eq. (7), Eq. (A1.7) can be rewritten as

M(k+1)=I—[I—K(k+1)H(k+1)]P(k)P1 =

+K(k+1)H(k+1)[P(k)P;h]

[I_P(k)1ç1]

Considering Eq. (10), we obtain M(k

+

1) =M(k) +K(k

+

1)H(k

+

1)[I—M(k)]

(11)

Appendix 2: The derivation of Eq. (12) In this appendix, we derive Eqs. (12) in context. First, we define the measurement residual vector as

2(k) =Z(k) —2(k)

(A2.1)

where Z(k) is defined in Eq. (1), and

2(k) +H(k)1(k) Then E{2(k)Z(k)t} =E({H(k)[x_2(k)}

+ V(k)}{H(k)[x—1(k)]

+

V(k)}t)

=H(k)P(k)Ht(k) +R(k) We assume that the data noise correlation matrix R(k) has the following form

R(k) =r(k)I Then

E[2(k)t2(k)]

=

tr E[i(k)2(k)t]

=

tr[P(k)Ht(k)H(k)]

+r(k)k

(A2.2)

326

J. Qu, TL. Teng/Physics of the Earth and Planetary Interiors 86 (1994) 301 —327

Now, supposing we get the additional measurement z(k

E[i(k+1)t2(k+1)]

+

1), we have

+2(k)2(k)tI

=trE[2(k+1)2(k+1)

(A2.3)

Meanwhile, we have E[i(k + 1)ti(k

+

1)1

=

tr E{P(k + 1)[Ht(k

=

tr[P(k + 1)Ht(k

+

+

1), Ht(k)] [H(k±1)

1)H(k + 1) +P(k

+

J}

+

r(k + 1)(k +1)

1)Ht(k)H(k)]

+r(k+1)(k+1)

(A2.4)

To obtain a pseudo-recursive expression for r(k + 1), we approximate P(k + 1) by P(k) in Eq. (A.24). Then E[2(k + 1)t2(k

+

1)J =tr[P(k)Ht(k

+

1)H(k

+

1) +P(k)Ht(k)H(k)}

+r(k

+

1)(k + 1) (A2.4’)

From Eqs. (A2.2), (A2.3) and (A2.4’), we have

E[2(k

+

1)2(k + 1)]

+

r(k + 1)(k

+

1)

—

r(k)k

+

tr[H(k

+

1)P(k)Ht(k + 1)]

We know E[i(k+1)2(k+1)J=z(k+1)—2(k+1)=z(k+1)—H(k+1)1(k) Therefore {[z(k +1)— H(k + 1)~(k)]t[z(k + 1)— H(k + 1)~Q(k)]—tr[H(k + 1)P(k)Ht(k + 1)]— r(k)} r(k+1)°°r(k)+

(k+1)

(A2.5) Considering z(k + 1) as a vector, we can obtain similarly r(k+1)=r(k)+

{[z(k + 1)— H(k + 1)2(k)]t[z(k +1)— H(k + 1)~Q(k)]—tr[H(k + 1)P(k)Ht(k + 1)]— r(k)} .

(k+1)ditn[z(k+1)]

(12)

References Backus, G. and Gilbert, F., 1970. Uniqueness in the inversion of inaccurate gross earth data. Phil. Trans. R. Soc. London, 266A: 123—192. Claerbout, J.F. and Muir, F., 1973. Robust modeling with erratic data. Geophysics, 38: 826—844. Deregowski, S.M., 1971. Optimum digital filtering and inverse filtering in the frequency domain. Geophys. Prospect, 19: 729—768. Franklin, J.N., 1970. Well-posed stochastic extensions of ill-posed linear problems. J. Math. Anal. Appl., 31: 682—716. Gordon, R., Bender, R. and Herman, G.T., 1970. Algebraic reconstruction techniques (ART) for three-dimensional electron microscopy and X-ray photography. J. Theor. Biol., 29: 4741—4781. Hartzell, S., 1989. Comparison of seismic waveform inversion results for the rupture history of a finite fault: application to the 1986 North Palm Sprints, California, earthquake. J. Geophys. Res., 94: 7515—7534. Helmberger, D. and Wiggins, R.A., 1971. Upper mantle structure of midwestern United States. J. Geophys. Res., 76: 3229—3245. Helmberger, D. Dreger, D., Stead, R. and Kanamori, H., 1993. Impact of broadband seismology on the understanding of strong motions. Bull. Seismol. Soc. Am., 83: 830—850. Humphreys, E., Clayton, R.W. and Harger, B.H., 1984. A tomographic image of mantle structure beneath southern California. Geophys. Res. Lett., 11: 625—627. Kennett, B.L.N., 1980. Seismic waves in a stratified half space—Il: theoretical seismograms. Geophys. J. R. Astron. Soc., 61: 1—10.

J. Qu, TL. Teng/Physics of the Earth and Planetary Interiors 86 (1994) 301 —327

327

Kennett, B.L.N. and Kerry, N.J., 1979. Seismic waves in a stratified half space. Geophys. I. R. Astron. Soc., 57: 557—583. Lines, L.R. and Treitel, S., 1984. Tutorial: a review of least-squares inversion and its application to geophysical problems. Geophys. Prospect., 32: 159—186. Magistrate, H.W., 1990. Three-dimensional velocity structure of Southern California, Part II. Ph.D. Thesis, Caltech, Pasadena. Mason, I.M., 1981. Algebraic reconstruction of two-dimensional velocity inhomogeneity in the High Hazles seam of Thoresby colliery. Geophysics, 46: 298—308. McMechan, GA., 1983. Seismic tomography in boreholes. Geophys. J. R. Astron. Soc., 74: 601—612. Mendel, J.M., 1973. Discrete Techniques of Parameter Estimation: The Equation Error Formulation. Marcel Dekker, New York. Mendel, J.M., 1987. Lessons in Digital Estimation Theory. Prentice-Hall, Englewood Cliffs, NJ. Oldenburg, D.W., 1981. A comprehensive solution to the linear deconvolution problem. Geophys. J. R. Astron. Soc., 65: 331—357. Oldenburg, D.W., 1982. Multichannel appraisal deconvolution. Geophys. J. R. Astron. Soc., 69: 405—414. Rice, R.B., 1962. Inverse convolution filters. Geophysics, 27: 4—18. Sipkin, S.A., 1982. Estimation of earthquake source parameters by the inversion of waveform data: synthetic waveforms. Phys. Earth Planet. Inter., 30: 242—259. Sipkin, S.A., 1986a. Estimation of earthquake source parameters by the inversion of waveform data: global seismicity, 1981—1983. Bull. Seismol. Soc. Am., 76: 1515—1541. Sipkin, S.A., 1986b. Interpretation of non-double-couple earthquake source mechanisms derived from moment tensor inversion. J. Geophys. Res., 91: 531—547. Su, F., 1992. Study of earthquake source, propagation medium attenuation and recording site amplification using coda waves. Ph. D. Thesis, University of Southern California, Los Angeles, CA. Tarantola, A., 1987. Inverse Problem Theory. Elsevier, Amsterdam. Treitel, S. and Robinson, E., 1966. The design of high resolution digital filters. IEEE Trans. Geosci. Electron, GE-4: 25—38. Zeng, Y.H., 1991. Deterministic and stochastic modeling of the high frequency seismic wave generation and propagation in the lithosphere. Ph. D. Thesis, University of Southern California, Los Angeles, CA.