Elastic full-waveform inversion for earthquake source parameters

SUMMARY 2-D full-waveform inversion of double-couple earthquake sources is implemented. Temporally and spatially extended sources are represented by superposition of double-couples. The source parameters solved for are the spatial location, origin time, amplitude and orientation of each double-couple. The velocity and density distribution and source time function are assumed to be known a priori but may be arbitrarily complicated. The non-linear inverse problem is solved by iterative linear approximation. The Jacobian matrix elements for source depth and rupture angle are computed by wavefield extrapolation forward in time, while those for origin time and amplitude are computed analytically. A smoothing technique that results in faster convergence and avoids local minima associated with cycle skipping is applied at each iteration. A spatial sampling interval, between discrete sources, of one-quarter wavelength of the dominant shear wave is optimal for inversion if high uniqueness of the result is desired. The presence of a fault is inferred from the spatial continuity of the rupture solution, rather than being imposed a priori. The method is illustrated by successful application to three synthetic source models: a single double-couple, a single extended rupture and a double extended rupture. The resolutions of the source depth and origin time are higher, and their posterior covariances are lower than those of the amplitude and rupture angle at each source point. Source depth, origin time and amplitude are primarily determined by the data; the rupture angle is more strongly influenced by the a priori information.


INTRODUCTION
Estimations of earthquake source parameters from full-wavefield data, as linearized inverse problems (Stump & Johnson 1977;Spudich 1980), typically fall into one of two types. Estimations of rupture history (e.g. Wald, Helmberger & Hartzell 1990;Goldstein & Archuleta 1991;Hartzell, Stewart & Mendoza 1991;Mikumo & Miyatake 1978Wald & Somerville 1995) generally assume that the position and orientation of the fault plane are known a priori, and involve fitting of local or near-regional data. The solution is usually in terms of the superposition of a grid of virtual sources distributed in time and space over the fault plane. The second class of solution involves parametrization in terms of moment tensor components for an apparent source that is a point in space but complex in time, as recorded at teleseismic distances (e.g. Nabelek 1984;Campos et al. 1994). The latter has also been applied at shorter distances (Mao, Panze & Suhadolc 1994).
The algorithm presented below is unique in three aspects.
(1) The distribution of the propagation velocity used is specified on a grid and so is arbitrarily complicated. Most previous studies assume the earth to be very simple (typically 1-D) although a (very) few studies (e.g. Dreger & Helmberger 1993) have incorporated lateral velocity variations.
(2) The inversion is for the spatial and temporal distributions and orientations of the energy release. Thus, faults are inferred from the spatial continuity of the rupture solution, rather than being imposed a priori. The rupture velocity is spatially variable and is inferred from the temporal relation between initiation of adjacent energy releases.
(3) It uses a pseudo-spectral solution of the elastic wave equation for extrapolation, so automatically includes all body waves, surface waves and mode conversions.
The inversion uses standard linearized inversion technology and so provides resolution and covariance information for analysis of the quality of the solution. It incorporates wavenumber smoothing which produces solutions close to the global minimum by avoiding local minima associated with cycle skipping. It can be used for any data set that provides a reasonable sampling of the focal sphere of an earthquake.
Although this paper presents only elastic 2-D tests for local earthquake recordings (for data recorded within a few tens of kilometres from the epicentre), it is applicable to all scales, and can be directly extended to viscoelastic, anisotropic 3-D media by changing the wavefield extrapolation algorithm and the number of unknowns (but these are beyond the scope of this paper).

THEORY AND ALGORITHM
In 2-D, a double-couple source can be parametrized as 4x, 2, to, a, 0) 2 (1) where x and z are the horizontal and vertical positions respectively, to is the origin time, a is the amplitude of the double-couple force and 0 is the angle associated with the orientation of the double-couple (the local rupture direction). These are the main parameters in any source description; they are equivalent to a moment tensor solution for a pure doublecouple (for which the determinant of the tensor = 0), with the addition of source location and origin time. For a spatially and temporally extended earthquake, the composite source parameters correspond to a series of II double-couples, ... , s,)=(xI, zI, to,, a,, 4, ... , x,, zn, ton, an, on).
(2) As this is a discrete representation (necessitated by the computational grid), the total energy release per unit fault area is the integral over the energy associated with all the grid points in that area. The same energy release is provided by a small number of widely spaced virtual sources or by a larger number of very closely spaced ones, each of which is associated with a proportionately smaller energy. For inversion we define the source region to be larger than that which we anticipate exists, and solve for the source parameters of all points in this region. In the result, only the actual source region will correspond to significant energy. This inherently non-linear waveform inversion is solved as a sequence of linearized inversions. Define do as the observed wavefield and dk as the calculated wavefield at the kth iteration. Then is the residual wavefield at the kth iteration, ASk is the vector of model parameter perturbations, and J is the Jacobian matrix containing the derivatives of the calculated wavefield with respect to the source parameters.

14)
where Mi(Sk) denotes the displacement amplitude of the ith point in the wavefield calculated for the estimated source Sk, and t j is the jth element of the source parameter vector (2).
To stabilize the inversion process, the solution of (3) is obtained by the algorithm proposed by Tarantola & Valette (1982) and Jackson (1979Jackson ( , 1985 using where So is the starting guess for the source parameters, E is the data covariance matrix, superscript T denotes the matrix transpose and superscript -1 denotes the matrix inverse. Sk is the solution and D, is its a priori covariance (constraint) matrix at the kth iteration. Coefficent c ( is set to 1 for the first few iterations, and then to 0, as discussed below. The initial model (So) must contain the low-wavenumber features of each of the source parameters, as these lie in the null space of the frequency-and aperture-limited data. Intermediatewavenumber source information is extracted from the data by the inversion. Wavenumbers corresponding to frequencies higher than those in the data also correspond to the null space.
To reduce the possibility of convergence to a local minimum associated with cycle skipping, there must be some overlap between the highest wavenumber in So and the wavenumber corresponding to the lowest frequency in the data (see Xu, McMechan & Sun 1995). This is the basis of the wavenumber filter procedure described below.
The posterior covariance matrix of the parameters is Because the parameters are in different units and have different accuracies, normalization was applied to the covariance matrix for display; where oi is the standard deviation of the ith parameter. Eq. (5) was originally proposed to obtain a unique solution for an underdetermined problem, which is a joint solution of the data residuals (Ad), and the a priori constraints (So -Sk). The a priori constraints are used to damp the inversion in the first few iterations to avoid divergence from the neighbourhood of So, and then it is omitted (by setting CI to zero) for the following iterations, which is equivalent to setting the starting guess equal to the current estimated solution.
At the first iteration ( k = 0 ) S o -S k is undefined and so is set equal to the a priori standard deviation to provide appropriate damping. One assumption of eq. (5) is that the source parameters and data have Gaussian distributions with zero means; this may not be true (Rothman 1985). This inconsistency may be avoided by increasing the values of D k to relax the a priori constraints as iterations proceed; as this may lead to instability (Mora 1987), the a priori source covariances were increased only at the last few iterations where the residual curve becomes flat and model pertubations are small. This procedure also indicates which of the model parameters are primarily defined by the data, and which are primarily defined by the a priori coiisxaints. Iterations are performed until the rrns residual curve becomes substantially flat. The rms residual and the number of iterations needed to achieve it are data-dependent. The final level of the rms residual corresponds to the noise in the data. In the examples below, absolute convergence was never obtained, but acceptable solutions exist at the point of diminishing return (which occurs at 7-10 iterations). At the last iteration, singular value decomposition is performed to obtain the resolution and covariance matrices.

Calculation of E and Dk
The data covariance matrix E is fixed for all iterations, as it represents the a priori uncertainty in the data (i.e. the inherent ambient noise) plus the misfit associated with the composite of all effects in real sources and propagation that are not parts of the synthetic forward modelling. Assuming no correlation between amplitude uncertainties at each time sample, E is a diagonal matrix with each diagonal value equal to the variance of the corresponding data point. For synthetic data, E may be estimated as (Nabelek 1984), where N is the total number of data points, and I is an identity matrix. For field data, E can be replaced by real data variance values, or estimated by eq. (9) if no other information is available. The covariance (E) affects only the rate of convergence (through its damping effect) and minor details of the solution (assuming that the algorithm for forward modelling is able to reproduce the main source parameters), so a crude approximation is sufficient.
In waveform inversion, most local minima are associated with cycle skipping (Rothman 1985). This means that the correct correlation between the observed and the predicted waveforms will occur only when the starting model is sufficiently close to the correct solution. If the predicted waveform is more than a half-cycle out of phase with the observed data, parameter a:',",. updates will favour local correlation with the next adjacent peak, rather than the correct one. The a priori covariance matrix D, has diagonal elements that are the variances of each of the four unknowns at each source point. Parameters z and to primarily alter the phases of the wavefields (by time shifting) when perturbed, thus the standard deviations of z and to are chosen to be the halfwavelength and half-period, respectively, at the dominant frequency of the source time function (e.g. at the peak of the particle displacement or velocity spectrum, dependent on which type of data are being fitted). Parameters a and 0 primarily alter the amplitude of a wavefield when perturbed. The standard deviation of each source amplitude is set equal to the amplitude of the corresponding source in the starting guess (usually a constant for all), and the standard deviation of the rupture angle in the examples below is arbitrarily set to 10-15 degrees. The precise values used for these standard deviations are not crucial as their main function is to damp the solution as described above.

Calculation of the Jacobian
One obstacle for non-linear waveform inversion when a large number of parameters are being solved for is the computation of the Jacobian matrix. Using a finite-difference approach, the partial derivative (4) can be approximated by

M , ( S ) -M , ( S )
where model S' is the same as model S, with a small perturbation A t j of parameter t j , which is a representative element in vector (2). The size of parameter perturbations should be small to keep the approximation (10) accurate, but not so small that it produces numerical problems; we find that a 1.0   Horizontal position (km) MCx, t; m34l = BMCx, t; a(41. Eq. (14) is unstable when the amplitude a, approaches zero so numerically a lower amplitude limit was set to avoid floating point underflow.

Irnplementa tion
The pseudo-spectral method was used to extrapolate the wavefields. Compared with finite-difference solutions of the wave equation, the pseudo-spectral method requires, in each space dimension, as few as a quarter of the number of grid points of fourth-order finite-differencing for equivalent accuracy (Fornberg 1987). Free-surface (zero-stress) boundary conditions are applied on the top of the model (Kolsky 1963), and absorbing boundary conditions are applied to the sides and bottom of the grid (Cerjan et al. 1985).
For each example below, the computational grid is 256 points (horizontally) by 128 points (vertically); the spatial grid increment for the first two models is 0.2 km; for the third, 0.1 km. Synthetic two-component (vertical and horizontal) displacement seismograms are generated for every fourth grid point (exclusive of the edge grid points used for the absorbing boundary condition) on the free surface of each model. Seismograms are constructed by saving one time sample at each finite-difference time step at those grid points that correspond to receiver locations, and plotting these as functions of time.
The source definition is based on Aboudi (1971), with the first derivative of a Gaussian displacement (Kelly et al. 1976). The spatial distribution of a source is defined using a summation over a grid of point dislocations (e.g. Spudich & Archuleta 1987), which is the discrete form of the integral representation theorem (Burridge & Knopoff 1964;Backus & Gilbert 1968;Olson & Apse1 1982). In practice, a discrete distribution of sources will appear continuous if they are less than a quarter-wavelength apart (e.g. Sheriff & Geldart 1982). Conversely, the most parsimonious parametrization of a composite source consists of individual points that are a quarter-wavelength apart; this is an important consideration in maximizing resolution, minimizing inter-variable interactions, and decreasing the total number of unknowns. This is one of the points illustrated in the examples.
The examples are 2-D, in which a point source is equivalent to a line source in 3-D, and a finite-length fault represents a plane source in 3-D. The algorithm itself can be directly extended to 3-D sources by addition of the appropriate unknown parameters and by using data recorded on an areal array rather than along a line. Previous studies (e.g. Wald et al. Mikumo & Miyatake 1993) have shown that this areal array need not be dense or evenly sampled to provide a good solution; only a pertinent sampling of the focal sphere is needed.

SYNTHETIC EXAMPLES
In this section, we present a series of examples with increasing structural and source complexity. These illustrate the procedures and characteristics of the inversion algorithm for earthquake source parameters.

Single double-couple source in a homogeneous medium
The simplest model to begin with is a single double-couple source in a homogeneous half-space (Fig. 1, upper panel). The dominant frequency of the source is 1.0 Hz. Fig. 2  Iteration number Figure 9. Normalized root-mean-square wavefield displacement amplitude residuals as a function of iteration number. The dashed curve is the solution without, and the solid curve is the solution with, the wavenumber filter applied.
at each source point is fixed. The source amplitude distribution (in Fig. 2 ) after eight iterations shows that most of the source energy in the best-fit solution is correctly associated with the relative horizontal location of the actual source (at 25.5 km).
The other parameters at this location also converged to the correct values with high resolution and very low posterior covariance (Fig. 3). The amplitude resolution at all five points is between 0.7 and 0.9, and the posterior variances of the amplitudes are generally lower than those of the initial (a priori) constraints, indicating that the amplitudes are well resolved by the data. The variance of the amplitude of the source point (at 25.5 km horizontal position) in Fig. 2 is higher than that in So (Fig. 3); the information in the data overcomes the effect of the relatively poor amplitude estimate in So by finding a solution near the constraint provided by the a priori variance, after the allowable variances are increased in the last two iterations. Estimated values of z, t, and 0 after eight iterations at the locations that do not correspond spatially with the correct solution have low resolutions and high variances (Fig. 3) because the amplitudes associated with these locations are small (Fig. 2). The smaller the amplitude is at a point, the closer to the null space the other parameters are, and the more difficult they are to resolve from the data.
The resolution matrix (Fig. 4) represents the uniqueness of the solution. This matrix is asymmetric because of the nonzero a priori covariance Dk in eq. (6), which is contained within C, in the resolution computation in eq. (7). This asymmetry is a result of the non-reversibility of the directions of action of the a priori constraints. For example, box X indicates a strong dependance of dip angle on origin time, because the origin time is relatively tightly constrained by the a priori variance information (Fig. 3). On the other hand, box Y indicates a negligible dependance of origin time on dip angle, because the dip angle is only loosely constrained by the a priori information.
Corresponding relations apply to all other parameter pairs.
In the resolution matrix in Fig. 4, the parameters z and t (which primarily affect time shifts, so we refer to them as phase variables) of the points spatially corresponding to the correct solutions (at the centre of each submatrix) have smaller offdiagonal values than those of the amplitude variables a and 0. Phase variables carry traveltime information associated with low wavenumber variations in source parameters; these are well constrained by the data, so have a high resolution. The sensitivity of the amplitude of a source is proportional to the reciprocal of its amplitude (eq. 14). Thus, wavefield residuals are more sensitive to low-amplitude sources, which have higher resolution and lower variance. The two points with higher resolution in off-diagonal submatrix A (for the two locations Horizontal position (km) Figure 11. A two-segment temporally and spatially extended source (white circles) in the heterogeneous model (top) produces vertical and horizontal displacement seismograms (centre and bottom). The grey scale shows only the P-wave velocity; S velocites are obtained from these by multiplying by 3 ~ *''.
adjacent to the source) and the off-diagonal submatrices on Strong correlations exist between amplitude variables a and the bottom row (for the angle parameter) indicate that these $ of spatially adjacent points in the normalized posterior interactions have a uniqueness imposed upon them by the covariance matrix (Fig. 4); see submatrices B and C. Similarly, a priori constraints, rather than by the data. strong correlations exist between phase variables z and to (see submatrices D and E), because the spatial aperture of the data is limited. This is similar to the depth/velocity ambiguity that exists in velocity estimation in reflection seismology. Increasing the aperture will reduce correlations between the phase variables. The correlations between a and 6' are small (see submatrices F and G) because of the geometry of the double-couple; the wavefields have nodal planes that are controlled only by the rupture angle 0. The correlations between phase variables and amplitude variables are also small (see the unlabelled offdiagonal submatrices) because they deliver different types of information from the data. For the last two iterations of the eight carried out in this example, the values of the model a priori constraint covariances D, were increased by a factor of four. From iterations six to eight the solution does not change significantly, but the posterior covariances do. For the points that do not correspond spatially to the correct solution, the posterior covariances of z, t, and 0 are already higher than the initial variances (Fig. 3) used for iterations one to six. The posterior covariances of the amplitudes are still near, or less than, the initial variances. This indicates that after eight iterations the values of parameters z, to, a and 6' for the source point are resolved by the data. Parameters for the other points are strongly constrained by the a priori variances. For completeness, the wavefield residuals are shown in Fig. 5 and the normalized root-mean-square residuals are shown in Fig.6; both show that the source parameters obtained at iteration eight account for about 98 per cent of the data.

Reverse fault with finite rupture velocity in a homogeneous medium
Because of the heterogeneous mechanical properties (asperities or barriers) in a fault zone, rupture along a fault is not uniform (Aki 1979). Spatially and temporally extended earthquake sources may be simulated b p a series of double-couples placed locally parallel and perpendicular to the fault plane with an appropriate time and space distribution to represent (generally discontinuous) slip on a fault. The rupture starts from the lower end of the fault and travels upwards with a velocity of 85 per cent of that of the shearwave velocity of the medium. The dominant frequency is 1.0Hz. The synthetic seismograms are shown in Fig. 7. Compared with the seismograms for a single double-couple ( Fig. 1 ), these responses are more extended in time; both starting and stopping phases are visible in both P and S waves. 15 points with a separation of 0.28 km were used to represent the fault in both the modelling and the inversion (Fig. 8). After 10 iterations, the inversion converged (Fig. 9).
Cycle skipping is the main cause of local minima in waveform inversions (Rothman 1985). Cycle skipping occurs mainly in the higher wavenumber changes of the phase variables and is a consequence of spatial aliasing. Convergence to the correct (global) minimum can occur only for wavenumbers that are initially within a half-wavelength of that solution. Thus, application of a simple rectangular smoothing filter after each iteration helps to avoid local minima (Xu et al. 1995). The length of the smoothing filter is progressively reduced as iterations proceed to allow the higher wavenumbers to contribute to the solution only when they satisfy the half-wavelength criterion. Fig. 8 shows the result after 10 iterations with five-point filtering during the first two iterations, three-point during the third and fourth iterations, and no filter for the following iterations. The normalized residual curves (Fig. 9)  that convergence is faster by a factor of almost two when the wavenumber filter is applied. The resolution and covariance matrices of the solution after 10 iterations with the smoothing filter are shown in Fig. 10. The resolutions of the phase variables are higher and their covariances are lower than those of the amplitude variables. In this model, the distance between any two neighbouring source points is about 8 per cent of the dominant shear wavelength. Increasing the separation between sources will reduce the correlation between amplitudes and thus increase the uniqueness of the solution. This is demonstrated in the next example.

Two-segment reverse-fault model in a heterogeneous medium
Generally, the velocity model associated with a real earthquake is complicated, as a multifault system produces strong lateral velocity variations. Fig. 11 (upper panel) represents such a model, in which the main (normal) fault is reactivated as a thrust (after Letouzey, Werner & Marty 1990). The horizontal and vertical dimensions of the model are 25.2 and 12.0 km with a spatial grid increment of 0.1 km. The dominant source frequency is 2.0 Hz. The source contains two ruptured fault segments. The rupture velocity at each point is 85 per cent of the local shear-wave velocity of the medium. Seven doublecouples (one group of four and one group of three, separated by a three-point barrier) were used to generate the test data; 14 double-couples (Fig. 12) were used to parametrize the ruptures along the fault for the inversion. The distance between any two neighbouring points is about 25 per cent of the dominant wavelength for shear waves. Seismograms (Fig. 11) are simulated for recorders on the free surface. To test the inversion under more realistic conditions, 20 per cent of the traces were randomly removed, and 20 per cent of bandlimited (1-24 Hz) noise was added to the data.

Resolution Matrix Normalized Posterior Covariance Matrix
Depth Origin Time Amplitude Angle Depth OriginTime Amplitude Angle The inversion uses the damping and wavenumber filtering techniques described above to stabilize the inversion. The a priori variances of the model parameters were increased by a factor of four at the last two iterations (9 and 10) to test the dependence of the inverted model parameters on the data and on the a priori constraints. After 10 iterations (Fig. 12), all parameters are near the correct solution, but the estimated rupture angles are less accurate than the other parameters. The resolutions associated with the parameters of the two ruptured fault segments (see the left column in Fig. 13) have high values (mostly 20.8). The variances of the z, to and a values on the two ruptured segments are smaller than those specified a priori, but the variances of 0 for some points are beyond the a priori variances (see the right column in Fig. 13). The resolution matrix (Fig. 14, left) shows a clear diagonal trend, and the posterior covariance matrix (Fig. 14, right) shows little correlation between parameters. It is concluded that the phase terms z, to and amplitude a are well resolved from the data, but the dip angles 0 are mainly controlled by the a priori constraints.
The diagonal dominance of resolution in this example (Fig. 14) is significantly higher than that in the previous example (Fig. 10); the posterior covariance of all parameters is lower for the spatial locations with high source amplitudes. The uniqueness of this model is higher than that of the previous model because the separation between neighbouring points is a larger fraction of the dominant shear wavelength (about a factor of four). The correlation between the amplitudes of adjacent source points is high (as seen by the width of the covariance diagonal in Fig. 10). The width of the diagonal corresponds to the theoretical limit of resolution of a quarterwavelength (Sheriff & Geldart 1982). In Fig. 10, the source points are spaced more closely together than can be resolved at the frequencies generated by the source, so the corresponding off-diagonal elements in Fig. 14 are low. Sources more widely spaced would begin to appear as separate sources. A continuous fault rupture is optimally parametrized (with a minimum number of parameters) on a set of spatial points about one-quarter of a wavelength apart. For completeness, the wavefield residuals at iterations 0 and 10 are shown in Fig. 15, and the normalized root-mean-square residuals are in Fig. 16; the source parameters defined at iteration 10 account for about 70 per cent of the data, which is reasonable, considering that 20 per cent of the data is incoherent noise that ideally always remains in the residual. Often the amount of data available is considerably less than that used here. To evaluate the behaviour of inversion under these conditions, the 32 noisy traces in Fig. 11 were reduced to 10 traces by keeping only every third trace. The quality of the 10-trace solution (shown in Fig. 17), relative to that for 32 traces (in Fig. 12) is determined from the corresponding resolution and covariance matrices; compare Figs 18 and 14. For the 10-trace solution, both resolution and covariance are less diagonally dominant (but with the same overall pattern), indicating interactions between adjacent source values. The most prominent difference is in the resolution and variance of the estimates of the source amplitude; in the 10-trace solution, the resolution of the amplitude is lower (compare boxes labelled A) and the covariance is higher (compare boxes labelled B). It appears that although the amplitude estimates are close to the correct values (Fig. 17), they are less constrained by the reduced amount of data. On the other hand, the resolution and covariance of the phase parameters (depth and origin time) are almost unchanged (see boxes labelled C and D in Figs 14 and 18); these appear to be less sensitive to the amount of data. The dip angle is less well estimated than the other parameters for both 32 and 10 traces, and the resolution and variance are correspondingly low and high, respectively. Dip angle is probably most sensitive to the total aperture in which the data are recorded, and reducing the number of data traces with the same aperture has very little effect. The remaining boxes in Figs 14 and 18 display similar relations with more interdependence when the amount of data is reduced. Horizontal position (km) Figure 17. The solution of the same problem as in Fig. 12, but using 10 rather than 32 traces. See caption to Fig. 12 for details. The corresponding resolution and covariance matrices are in Fig. 18  of inversion of earthquake parameters from two-component data recorded on a free surface. The source time function could be solved as part of the inversion, so this assumption is not necessary, and was used here only for convenience. The inversion algorithm, combined with a progressive smoothing filter applied to the solution at each iteration, is stable and converges in a few iterations. The filter increases the rate of convergence by about a factor of two, and reduces the likelihood of the solution getting stuck in a local minimum due to cycle skipping.
A spatial sampling interval, between discrete sources, of onequarter of a wavelength of the dominant shear wave is optimal for inversion if high uniqueness of the result is desired; compare Figs 10 and 14. At each point, the resolutions of the source depth and origin time are higher and the posterior covariances are lower than those of the amplitude and the rupture angle. Source depth, origin time and amplitude are primarily determined by the data; the rupture angle is more strongly influenced by the a priori information. Increasing the a priori covariances of the model parameters at the last few iterations gives information on the dependence of the inverted parameters on the data.
The number of traces used per source (up to 40) is realistic. For example, see the data from the SMART-1 array in Goldstein & Archuleta (1991). Dreger & Helmberger (1991) suggest that full-wave inversion of data from one station can yield useful information. Also, through the IRIS/PASSCAL programme in the United States there are hundreds of portable three-component recorders available. The specific configuration of a line of recorders above an event is also straightforward to deploy in the case of aftershocks where good locations are available from previous events. In general, an areal rather than a linear deployment of recorders would be more characteristic; inversion of such data implies the extension of the current algorithm to 3-D. Fig. 17 shows that widely spaced (e.g. regional array) data should also be amenable to full-wave inversion. Borehole data may also be included as there is no restriction on receiver locations.
We have assumed in our inversion that the velocity structure is known a priori. Where such velocity information does not exist, or is poorly known, the uncertainties in the velocity will be translated into corresponding uncertainties or bias in the estimated source parameters. Fortunately, regions of high seismicity tend to be well instrumented and are the focus of a variety of refraction, tomography and reflection velocity studies. Examples in the United States include the Los Angeles Basin and the New Madrid seismic zone.
The inversion described above is independent of the scale of the problem. We used the local earthquake context and a 2-D recording configuration for convenience. Computational limitations may force lower frequencies to be used at larger distances.
The parametrization is not applicable for all events as only double-couples are included; it could be made more general by adding compensated linear vector dipole (CLVD) terms. Because of this limitation, the best-fit double-couple solution would be skewed and the final residual would be increased if the data included a significant CLVD contribution.
One advantage of using the pseudo-spectral solution of the wave equation is that it allows computation of the Jacobian matrix for arbitrarily complicated geological models including strong spatial velocity and density variations, anisotropy and viscoelasticity. All multiples, converted waves and surface waves are also automatically included. Other extrapolators such as finite-element or finite-difference methods may be substituted directly. The extension from 2-D to 3-D is also conceptually straightforward; the uncertainty associated with any parameter, and the interactions between parameters in 3-D, can be explicitly determined by analysis of the solution residuals, resolution and covariance, as illustrated here for the 2-D case. authors acknowledge the constructive comments of George Helffrich and two anonymous reviewers. This paper is Contribution No. 837 from the Program in Geosciences at the University of Texas in Dallas.