Combination of Linear Inversion and Nonlinear Optimization for Hypocenter and Velocity Estimation

Abu M. Asad, John N. Louie and Satish K. Pullammanappallil
Seismological Laboratory, Mackay School of Mines, University of Nevada, Reno, NV 89557; URL: http://www.seismo.unr.edu

Presented as a poster to the American Geophysical Union Fall Meeting, San Francisco, Calif., December 9, 1994.

Introduction

Our objective is to locate an aftershock sequence in the Eureka Valley area of the Western Great Basin, and simultaneously estimate the local three-dimensional velocity structure. Occurrence of a magnitude 6.1 earthquake followed by several hundred aftershocks, recorded by permanent networks and portable stations, offered an opportunity for study of an area that has lacked detailed investigation. A better definition of the local structure would be a step forward in understanding small structures in the Great Basin, which in its turn would reduce the non-uniqueness of regional models.

(above) Location map of the California-Nevada border region showing major fault traces and epicenters of events located by UNR within 6 months of the May 19, 1993 Eureka Valley M6.1 earthquake. Red triangles show permanent seismic network stations maintained by UNR and Caltech.
(below) Detail map of Eureka Valley area showing locations of portable array stations set up by UNR (red triangles), and grid nodes used in 3-d velocity estimations (green grid). Sections A-A' and B-B' (below) show velocity results; colored dots are network epicenters.
(PDF file available)

We investigate structure and aftershock distribution with a parallel approach to the simultaneous hypocenter location and three-dimensional velocity estimation problem. We invert with both three-dimensional linear inversion (Thurber, 1983; Eberhart-Phillips, 1990) and three-dimensional nonlinear optimization by simulated annealing. (Pullammanappallil and Louie, 1994).

P-arrival (blue) and S-arrival (red) time picks from the UNR portable array, versus distance from epicenters given by the UNR permanent network.
(PDF file available)

Linear Inversion

We use Thurber's (1983) approximate ray tracing (ART) and pseudo-bending (PB) algorithm (Eberhart-Phillips, 1990) for forward modeling of P- and S-wave arrival times in an iterative damped least-squares inversion for hypocenters and three-dimensional velocity structure. The earth structure is represented in three dimensions by velocity defined at discrete grid points. The ART algorithm selects the path with the least travel time from a suite of circular arcs connecting the source and receiver. The iterative PB method fine tunes the ray path obtained by ART to approximate better the true ray path dictated by local velocity gradients (Um and Thurber, 1987; Eberhart-Phillips, 1990).

Non-Linear Inversion

Our approach for simultaneous determination of hypocenters and velocity structure without any a priori assumptions is an outgrowth of the simulated annealing optimization scheme of Pullammanappallil and Louie (1994). For the forward problem of travel time computation, we use Vidale's (1990) finite-difference scheme. An iterative Monte-Carlo optimization (simulated annealing) process solves the inverse problem. Each iteration comprises the steps:

Travel time computation through the current model and determination of the least-square error;
Simultaneous perturbation of origin time, hypocenter location, and velocity structure, and computation of new least square error;
If the error is less than that of the previous iteration, then the corresponding model is always accepted; if the error is larger, the model is conditionally accepted based on a probability criterion;
The previous steps are repeated until the annealing converges, where the difference in the least square error between successive models and probability of accepting new models become very small.

Diagram showing steps in four methods used to relocate Eureka Valley aftershocks and estimate 3-d velocities.
(PDF file available)

Comparative Relocation Results

Maps showing relocated epicenters from portable array picks, using respectively: a fixed 1-d velocity model; simultaneous nonlinear velocity and hypocenter optimization by simulated annealing; 3-d relocation through the 3-d optimization velocity model result; and linearized relocation and velocity adjustments starting with the nonlinear optimization results. Epicenters are plotted against topography (brown contours) and Bouguer gravity (blue contours) in the Eureka Valley area. Epicenters of the M6.1 main shock from 3 networks are shown as the open square, brown triangle, and brown diamond. A-B and C-D show the paths and the events plotted in the sections below.
(fixed 1-d & nonlinear, and fixed 3-d & linearized PDF files available)

Section C-D showing relocated epicenters from portable array picks, using respectively: a fixed 1-d velocity model; simultaneous nonlinear velocity and hypocenter optimization by simulated annealing; 3-d relocation through the 3-d optimization velocity model result; and linearized relocation and velocity adjustments starting with the nonlinear optimization results. Hypocenters (colored) are projected to the south-north line of section from within the pink areas on the epicenter maps above. Epicenters of the M6.1 main shock from 3 networks are shown as the open square, brown triangle, and brown diamond.
(fixed 1-d & nonlinear, and fixed 3-d & linearized PDF files available)

Section A-B showing relocated epicenters from portable array picks, using respectively: a fixed 1-d velocity model; simultaneous nonlinear velocity and hypocenter optimization by simulated annealing; 3-d relocation through the 3-d optimization velocity model result; and linearized relocation and velocity adjustments starting with the nonlinear optimization results. Hypocenters (colored) are projected to the west-east line of section from within the pink areas on the epicenter maps above. Epicenters of the M6.1 main shock from 3 networks are shown as the open square, brown triangle, and brown diamond.
(fixed 1-d & nonlinear, and fixed 3-d & linearized PDF files available)

Topographic (brown contours) and gravity (blue contours) map of Eureka Valley, with final epicenters (color; ASCII list available) and velocity grid locations (green) and section A-A' and B-B' locations (black). Velocities and hypocenters found by simulated-annealing optimization were used as a starting model for linearized velocity and hypocenter adjustments. Epicenters of the M6.1 main shock from 3 networks are shown as the open square, brown triangle, and brown diamond.
(PDF file available)

Sections A-A' and B-B', respectively, through the final 3-d velocity model. Velocities and hypocenters found by simulated-annealing optimization were used as a starting model for linearized velocity and hypocenter adjustments. Velocities below the red lines at about 9 km depth cannot be constrained by the simulated-annealing optimization, which produces the false velocity inversion. Note how the shallow low-velocity region follows the gravity low of the valley fill in the map above, and the deeper high velocities track the trend of the gravity high below the mountains to the southwest.
(PDF file available)

Conclusions

In an area like Eureka Valley, California, where no effective a priori information on the velocity structure is available, it is a difficult job to choose a starting velocity model that would satisfactorily locate the hypocenters and converge to the ``correct'' model in linear inversion schemes. Our tests show that nonlinear optimization by simulated annealing can render a worthwhile estimate of the initial velocity model, which can be used for joint hypocenter-velocity structure inversion by a linear scheme.
We relocated about 150 aftershocks using 673 P and 433 S arrival times. Though the hypocenters do not line up very well, we can mentally pass an almost vertical plane dipping eastwards through them. This plane is consistent with one of the nodal planes of the main shock focal mechanism. It is also consistent with the findings of radar interferometry work by Peltzer and Rosen (1995). The final velocity model is in tune with the gravity and topography.
An important observation of the study is that the random and independent perturbation of velocity and location parameters in the simulated annealing scheme can make a situation possible where we obtain good estimates of the velocity model, but estimates of hypocentral parameters that are not plausible. Therefore, we recommend a tandem combination of nonlinear and linear schemes for poorly-characterized areas.

References

Eberhart-Phillips, D., 1990, Three-dimensional P and S velocity structure in the Coalinga region, California: J. Geophys. Res., 95, 15342-15363.

Peltzer, G. , and Rosen, P., 1995 in press, Surface displacement of the May 17, 1993 Eureka Valley, California earthquake observed by SAR Interferometry: Science.

Pullammanappallil, S. K., and Louie, J. N., 1994, A generalized simulated-annealing optimization for inversion of first-arrival times: Bull. Seismol. Soc. Amer., 84, 1397-1409.

Thurber, C. H., 1983, Earthquake locations and three-dimensional crustal structure in the Coyote Lake area, central California, J. Geophys. Res., 99, 8226-8236.

Um, J., and Thurber, C., 1987, A fast algorithm for two-point seismic ray tracing, Bull. Seis. Soc. Am., 77, 972-986.

Vidale, J. E., 1990, Finite-difference calculation of travel times in three dimensions, Geophysics, 55, 521-526.