Abstract

Introduction

The objective of joint inversion is to minimize the least-square difference between the observed and calculated times through a model. This is cast as a  generalized nonlinear least-squares problem (Tarantola and Valette, 1982; Tarantola, 1987). Farra and Madariaga (1988) solve this problem iteratively using the Gauss-Newton method while Williamson (1990) uses an iterative subspace search method to do the same. van Trier (1990a, b) simultaneously inverts for the reflector geometry and interval velocities using a Gauss-Newton method.  Reflection times depend on both the velocity and the depth of the reflectors, making traveltime inversion a nonlinear process. All the above methods involve local linearization of the problem, and require the starting model to be close to the desired final solution. We demonstrate this requirement later using synthetic models.

In this paper we examine the use of a nonlinear optimization scheme, namely simulated annealing, to invert reflection traveltimes. We must rapidly compute traveltimes through models by a method avoiding raytracing. It employs a  fast finite-difference scheme based on a solution to the eikonal equation (Vidale, 1988). This accounts for curved rays and all types of primary arrivals, i.e., the fastest arrival, be it a direct arrival or a diffraction. We test the optimization process on synthetic models and compare its performance with linear inversion schemes that use curved rays. Synthetic examples and optimizations of data demonstrate that the method is not dependent on the choice of the initial model. Another advantage of the simulated-annealing algorithm is that it produces a suite of final models with comparable least-square error. This enables us to choose the model most likely to represent the geology of the region. We also use this property to determine the uncertainties associated with the model parameters obtained.

Lastly, we use the scheme to invert reflection times to image the Garlock  fault and the velocity within the Cantil Basin below Fremont Valley, California. The traveltime picks are made from the shot records of COCORP Mojave Line 5. Strong lateral velocity variations across the basin make it difficult to image the Garlock fault by standard seismic velocity analysis techniques. The use of prestack Kirchhoff migration (Louie and Qin, 1991), which accounts for curved rays and lateral velocity variations, requires prior knowledge of the velocities. This makes it important to use an inversion scheme that will image the velocities with minimum a priori constraints. In addition, the geometry of the Garlock fault at depth is not unequivocally known. The nonlinear optimization scheme we use can reconstruct velocities, as well as recover estimates of dip or length of the fault. We need not specify the length or dip prior to the inversion.

Methods

Traveltime calculation

Optimization by simulated annealing

We make use of a variation of this optimization process called generalized simulated annealing (Bohachevsky et al., 1986; Landa et al., 1989). Its basis is a Monte Carlo technique due to Metropolis et al. (1953). The algorithm essentially comprises the following steps:

1) Compute traveltimes through an initial model. Determine the least-square error, E₀. For any iteration i, we can define the least-square error (which is the objective function):
 

E_i = (1/n)(S(t_j^obs - t_j ^cal)²),                                       (1)

where n is the number of observations, j denotes each observation and t^obs and t^cal are the observed (data) and calculated traveltimes respectively. The summation goes from j=1 to j=n.

2) The next step is to perturb the model. To preserve the nonlinearity of the problem, we carry out simultaneous perturbations of the medium velocity and the reflector depth. We perturb the velocity by adding random sized boxes, followed by smoothing with averaging over four adjacent cells. The boxes can vary between one cell size and the entire model size. To perturb the reflector we add random length lines that are smoothed by averaging adjacent nodes. Again, the added lines can be as small as one grid spacing or as long as the length of the model. Following the perturbation, we compute the new least-square error E₁.

3) The following criteria are put to use to determine the probability P of accepting the new model:

 
P = 1 ; E₁ <= E₀                                                      (2)

P = P_c = exp{( E_min - E₁ ) ^q  DE / T}; E₁ > E₀          (3)

4) Repeat steps 3 and 4 till the annealing converges , where the difference in least-square error between successive models and the probability of accepting new models becomes very small. We discern the latter condition when there are a large number of iterations (50,000 or more) over which  no model is accepted.

Determination of the annealing parameters

The convergence of the inversion scheme is sensitive to the value of q and the rate of  cooling. The rate of cooling refers to the variation of parameter T with iteration to aid convergence. One iteration includes one perturbation of both velocity and reflector depth. Rothman (1986) found that simulated annealing finds the global minimum when one starts at a high temperature, rapidly cools to a lower temperature and then performs a number of iterations at that temperature. This low temperature, also referred to as the critical temperature (T_c), must be high enough to allow the inversion to escape local minima, but low enough to ensure that it settles in the global minimum.

We use a procedure based on the one developed by Basu and Frazer (1990) to determine the critical temperature. We determine the parameter before performing the inversion, and it may not be the same for a different inversion problem. The objective is to perform a series of short runs for a range of T values. In our case, 1000 iterations comprise a short run (a full run would be between 100,000 and 300,000 iterations).

1) First, for a fixed value of $T$, we compute the average least-square error, $E_{av}$, as shown below:

E_av(T) = (1/n) (S E_k )                                                     (4)

2) Next, repeat the above step for a range of T values. Typically, we vary T from 10^-6 to 10³ by factors of 10 (thus resulting in 10 values).

3) Then make a plot of E_av that we obtain for each T, versus log T.

4) The value of T which corresponds to the minimum E_av is the critical temperature, T_c (Figure 2a).

Having found T_c we can proceed to fix the cooling rate. The initial temperature is high (T = 10 in our inversions). This destroys any order in the model, so as not to bias the inversion. Within a few thousand iterations we rapidly decrease T to the critical temperature. On reaching T_c, we slow the rate of decrease of T, enabling the optimization process to settle in the deepest minimum. Figure 2b shows a typical cooling rate curve. Here, we keep the temperature constant at T_c from the 2000 th to the 50,000 th iteration, and then decrease the rate. After some trial runs we found this to work best, though it is still a guess.

To determine the other empirical parameter, q we do the following:
5) At the critical temperature, T_c, compute the average least-square error, E_av for a short run of the optimization process.

6) Repeat the above step for a suite of q values.

7) The value of q which corresponds to the minimum E_av will give the
optimum convergence (Figure 2c).

Results

We use two synthetic models in our study, namely the  three box model (Figure 3a) and the lateral velocity gradient model (Figure 5a). The models are 30 km long and 8 km deep. There are sharp velocity contrasts in the three box model, with slowness ranging from 0.75 s/km to 1.25 s/km. In contrast, the velocity variation is gradual in the lateral gradient model (patterned after Williamson, 1990). Here the slowness varies between 0.382 s/km and 0.418 s/km. In both cases the true reflector is fault-shaped, ranging between depths of 7 km and 5 km. There are 30 sources each with 30 receivers along the top of the model, thus resulting in 900 traveltime observations.

We first determine the critical temperature T_c. For the two synthetic cases T_c was equal to 0.1. Accordingly we fix the cooling rate. Finally, before performing the inversion, we determine the value of q. In our case q was equal to 2 for optimum convergence. Figure 2 illustrates the above steps.

Optimization results of the three box model

Optimization results of the lateral gradient model

The initial model is one of constant slowness, with a value of 0.3 s/km, and  has a flat reflector 6 km deep. During the inversion, the slowness can  assume values between 0.3 s/km and 0.5 s/km. Figures 5b and 5d  show the annealing results. The optimization scheme does a poor job of recovering the small lateral variations; it reconstructs only the  average velocity. On the other hand, it images the reflector depth and shape very well.  The recovered reflector is a smoothed version of the true reflector,  whereas the linearized scheme recovers a heavily averaged shape. The least-square error decreases from 2.97 s² for the initial model to 0.008543 s² for the final model.

Comparison with a linearized inversion scheme

We compare the performance of annealing with a linearized inversion scheme. In the linearized inversion, we linearize the nonlinear relationship between traveltime and slowness and the reflector depth. Generally, this linear scheme is cast in a matrix form and the solution obtained using a standard least-squares technique (Jackson, 1972; Wiggins, 1972). We use a singular value decomposition method to do this. Instead of solving for the perturbations to the model parameters, we compute a new model during each iteration (Shaw and Orcutt, 1985). We compute the traveltimes by the same method we used in the annealing. To speed up convergence and stabilize the inversion we smooth the velocity perturbations using a Laplacian operator (Lees and Crosson, 1989) and the reflector depth by a second-difference operator (Ammon et al., 1990).

Figures 3 and 5 compare the results we obtain from the linearized inversion with those from the nonlinear optimization scheme. Neither method recovers the velocity well when there is a gradual lateral variation in the true model (Figure 5).  The range of the reflection traveltimes through the lateral gradient model is only about 8 seconds while it is 23 seconds through the  three box model. Hence, the survey aperture might be too small to recover the slowness variations in the former case. Annealing does better in recovering the sharp features of the model - both in the case of velocity and of reflector shape. The smoothing applied during the linearized inversion scheme helps convergence, but smears the
velocities and averages out the reflector depth. Both methods reconstruct the reflector depth and shape fairly well close to the edges of the model. It is noteworthy that the velocity above these regions are not well recovered. This can be attributed to the inherent ambiguity in the simultaneous  reconstruction of velocity and reflector characteristics. It illustrates that velocity-depth trade-offs near the reflector are hard to resolve. In our tests the reflector depth is better resolved probably because it is described by fewer parameters (30 as compared to
240 velocity parameters). The second-difference smoothing applied during the linearized inversions and the averaging during annealing also helps in reconstructing the reflector shape better.

We also note that the convergence of the linearized inversions is initial model dependent (in the synthetic tests the starting model is close the true solution for the linearized inversions). To demonstrate this we do another synthetic test. We use the slowness model used by Louie and Qin (1991) for their prestack migrations (which will be put to further use later on) as the true model. The model has dimensions of 26 km by 7  km with a fault shaped reflector along the bottom as  shown in the top panel of Figure 6. The slowness varies between 0.2 s/km and 0.45 s/km.  We assume the reflector position and shape are known, and invert for only the velocities. Twenty-six receivers along the surface record traveltimes from twenty-six sources. Figure 6 shows inversion results of the linearized scheme using two different starting models. The solutions are different, though the least-square errors are comparable in both cases.  The inversion with an initial model of 0.2 s/km reconstructs the model better (a global minimum) than the one with a starting slowness of 0.4 s/km (a local minimum).

Figure 7 shows the optimization results by annealing. The final models we obtain by starting with two different models look similar. The random perturbations used during the annealing scheme make it independent of the initial model. It destroys any order in the starting model within the first  few iterations. The Laplacian smoothing applied during the inversion might explain the better reconstruction of the low velocity (high slowness) basin in the linearized scheme. Note that in both the methods, the region below the reflector is unconstrained due to absence of any ray coverage. Thus, even when we eliminate the nonlinearity due to simultaneous dependence of reflection traveltimes on the velocity field and reflector depth, by keeping the reflector position fixed, the linearized inversion scheme shows dependence on the initial model. Though the linearized inversion takes fewer iterations to converge, more time is spent in the matrix operations, making the total  computation time comparable in both cases.

For the nonlinear optimization scheme, the computation time is proportional to the size of the model and the number of reflector nodes. Models with 1500 cells and a reflector along the whole length took about 60 hours on a SPARCstation 2. With a grid spacing of 250 m and for a survey line 15 km long, this would mean we can invert for velocities and reflector depths down to 7 km in the above mentioned time.

Imaging the Garlock Fault, Cantil Basin

The left-lateral Garlock fault branches into two strands in the Fremont Valley, California, region. The eastern branch strikes almost east-west while the western branch strikes southwest-northeast. Cantil Basin itself is a pull-apart basin, and Aydin and Nur (1982) show that it is one of the largest of its kind in the San Andreas fault system.

The Consortium for Continental Reflection Profiling (COCORP) Mojave Line 5 collected seismic reflection data across the Garlock fault in Fremont Valley. Stacked sections of Line 5 displayed by Cheadle et al. (1986) show portions of the basin floor and the southern wall. But strong lateral velocity variations across the basin made the imaging inadequate. Louie and Qin (1991) used a prestack Kirchhoff migration, which accounts for curved rays and lateral velocity variations, on the Line 5 data. They infer that Cantil Basin was either a negative flower structure, or a detachment headwall with the Garlock fault flattening at depth. They prefer the detachment model since it allows a simple mechanism to explain the observed broadening of Cantil basin with fault movement over time.

Kirchhoff migration assumes prior knowledge about the velocity. Hence, an iterative application of this method would demand the initial model be close to the solution to keep it computationally inexpensive. Use of the nonlinear optimization will not require any constraints on the model. Since one of the goals is to determine the nature of the fault at depth, we do not put any constraints either on the location or the shape of the fault (the reflector). Resolving the structure of the Garlock fault and Cantil Basin would help constrain regional tectonics, mainly the configuration of the Mojave block.

Results of the traveltime inversion

The data we use are reflection times picked from shot gathers (Figure 8).  Louie and Qin (1991) compute synthetic seismograms to show that these are reflections off the east branch of the Garlock fault. The total number of observations is 621 reflection times, picked from 28 shot
gathers. In order to determine an optimum grid spacing to use in the inversions, we do a convergence test (not shown). We compute traveltimes through models with the same dimensions, but different grid spacings. For a range of grid spacings the optimum spacing is that for which smaller spacings did not appreciably alter the computed times (i.e.: less than the error in the data). The dimensions of the model were 14.4 km by 4 km with a grid spacing of 277.85 m. We incorporate elevation corrections (C. Ammon pers. comm.) by extrapolating the times to sources and receivers not on the grid, assuming plane wave incidence at the receivers.

Before embarking on the optimization process, it is essential to determine whether letting the reflector length vary hinders the inversion from
reaching the solution. For this study our true velocity model is the one used by Louie and Qin (1991) for the Cantil Valley region. The source-receiver configuration is same as the actual survey. The shape of the reflector corresponds to the dip of the east branch  of the Garlock fault shallowing at depth. It has a dip of about 50 degrees near the surface, changing to about 30 degrees at 2.5 km depth (Figure 9a). We do two optimization runs - one keeping the reflector length fixed and equal to the true length, and the other letting the length vary. From Figure 9b  we see that keeping the reflector length fixed hinders the optimization process. The velocity and reflector are reconstructed poorly. This is an  example of a bad choice of the objective function. The annealing process does a better job of reconstructing the velocities and reflector shape when the reflector length is allowed to vary (Figure 9c), suggesting that it is a more useful objective function. The optimization reconstructs
only the lower part of the reflector because the picked reflections see only this segment of the true reflector. We demonstrate this in a later section.

Hence during the data optimizations, we let the reflector length vary along with its position. The length can be as short as 277.85 m (the grid spacing) or as long as 14.4km (the length of the model). The reflector position can vary from near surface to 4 km depth. Hence, the reflector can assume any dip and the dip can change along its length. Each cell in the model can assume any velocity between 1.5 km/s and 6.3 km/s. This is the range of velocities expected in this region. We perform the optimization with three different constant velocity starting models - 2.5 km/s, 4.0 km/s and 8.3 km/s. We see that (Figure 10) the final model in all three cases are similar, even if we  start out with the unreasonably high velocity of 8.3 km/s. The features that are common to all three results would be the most reliable.

The inversion reconstructs a prominent low velocity basin (2.4 - 2.6 km/s) whose location corresponds to that of Cantil Basin. The basin seems to extend 1 - 2 km further south than suggested by Louie and Qin (1991)  and may be 2 to 2.5 km deep. The reflector is about a kilometer long and extends from a depth of 3 km to about 3.5 km. The dip of this segment of the reflector is between 20 to 30 degrees and it lies beneath the low velocity basin. If we interpret this as part of the Garlock fault, it indicates that the fault shallows at depth, since diffraction  studies (Louie and Qin, 1991) indicate that the fault dips about 45 degrees near the surface. The least-square error decreases from 3.477 s² for the constant velocity starting model to 0.0044465 s² for the final model. The traveltime residual between the final model and the data is less than 0.1 seconds for all the observations. The range of the observed traveltimes was between 2.25 and 3.65 seconds.

By using three different initial models, we allow the optimization process to sample different parts of the model space. It approaches the minimum least-square error from different directions. Hence, we can use the final models, with comparable least-square error, from each
of the three runs to estimate the uncertainty associated with the obtained model parameters. In our analysis, we use about 15 models from each run and thus a total of about 40 models. Figure 11a shows the reflector points sampled by the rays for these models along with the recovered reflector of the final model. Most of  the midpoints cluster around the recovered reflector. The standard deviation of the velocities (Figure 11b) range from 0.46 km/s to 0.85 km/s. We note here that these may not represent the actual deviations in the model parameters because of the nonlinear nature of the problem. It is difficult to separate out the contributions due to mislocation of the reflector and the errors in the velocities. But it does provide an idea about well- and poorly-resolved regions in the model. Our statistical analysis is over accepted models only. So there might be unconstrained regions in the model which are not perturbed, and hence have a small deviation. But, in general, the regions with higher standard deviation are parts of the model that have poor ray coverage. Hence, changing velocities in these regions do not affect the calculated traveltimes much and thus contribute minimally to the least-square error.

Our synthetic study using the velocity model constructed by Louie and Qin (1991)  for the region and the actual survey's source-receiver configuration (Figure 9c), indicate that the scheme recovers velocities in and around the low velocity basin to within 0.5 km/s to 1.0 km/s. The resolution decreases with distance from and depth below the basin. This approximately corresponds to the deviations we obtain from our statistical analysis on the models obtained using the real data. This study also shows that the optimization images only the lower part of the reflector, suggesting that the small segment of the reflector imaged during the data optimization is fairly robust.

We attribute the poor imaging of the reflector to the source-receiver configuration of the survey, that is, all the sources are towards the north side of the model and most of the receivers to the south. To confirm this, we determine the depth points along the reflector for every source-receiver pair in a manner similar to how we compute the traveltimes. We plot these along the true reflector in Figure 9a (delineated by a
square). This shows that the picked reflection rays see only the lower part of the reflector. Thus poor ray coverage causes the steeply-dipping
part not to be reconstructed and might also explain the poor reconstruction of velocities away from the basin.

Doing this study suggests what features of the model are robust and which are not. One has to note that we use only the primary reflected arrivals and a small subset of the data (only 28 shot gathers along Line 5) in our optimizations. In contrast, the prestack Kirchhoff migration (Louie and Qin, 1991) uses all types of arrivals and more data (228 shot gathers), to image the Garlock fault. So including more data in our  optimizations could enable us to image more of the fault.

Conclusions

We have presented a nonlinear optimization scheme to invert seismic reflection traveltimes, tested it on synthetic data and used it to image the Garlock fault and velocities within Cantil Basin.

Our synthetic study using the three box model and  lateral gradient model demonstrate the inherent ambiguity in recovering the velocity and reflector depth from reflection traveltimes. But the nonlinear optimization scheme avoids local linearization of the problem, allowing it to settle into a global (or deepest) minimum. The random perturbations that go into the optimization process make it independent of the initial model.
The method provides for easy implementation of constraints. It facilitates inversion for the reflector depth, its length and shape. The method does not involve matrix inversion; hence it is computationally tractable. The need to perform numerous iterations make the scheme computationally intensive when the depth of interest becomes greater than 10 km, for grid sizes we use in the above studies.

The velocities obtained by our optimization scheme within Cantil basin (2.4 - 2.6 km/s) are consistent with earlier studies using gravity data,
refraction data and Kirchhoff migration. The depth of the basin is about 2.5 km and seems to extend further south than suggested by previous studies. The synthetic study shows that the small segment of the reflector recovered by the optimization process is robust. This indicates that the dip of the Garlock fault shallows at 3.4 to 4 km depth. This is consistent with the detachment model suggested for the Mojave region by Louie and Qin (1991). Thus, using only the primary reflected arrivals and a small data set we were able to obtain velocities within the basin and image a segment of the fault.

Successful application of the nonlinear optimization scheme to image the  Garlock fault and Cantil basin suggests that the method is very useful in  regions where there are strong lateral velocity variations and where neither the fault shape nor depth are known. We demonstrate that inverting for the reflector length is crucial for the optimization to find the best solution. This feature makes the scheme more attractive, since generally we do not know the extent of subsurface structures. Though the recovered reflector is only a collection of depth points, it gives us an idea of the location and nature of the reflector. This information, along with the reconstructed velocity model, can be used in a prestack migration (which uses all arrivals) to recover finer details along the reflector.

Acknowledgements

The COCORP Mojave data were graciously provided by Sid Kauffman and Larry Brown of Cornell University. Part of S.K.P's research was supported by the Chevron Geophysics Fellowship while at Penn State University. The research was also supported by the National Science Foundation under grants EAR-9005302 at Penn State and EAR-9296125 at the University of Nevada. Chuck Ammon of U. C. Santa Cruz and Bob Clouser of Penn State  University provided valuable advice and discussions. The authors also  benefitted from comments by the Associate Editor Bill Harlan, by Matthew Brzostowski, and by two anonymous reviewers.

References

Ammon, C. J., and Vidale, J. E., 1990, Seismic travel-time tomography using combinatorial optimization techniques: Seismological  Research Letters, 61, 39.

Aydin, A., and Nur, A., 1982, Evolution of pull-apart basins and their scale independence: Tectonics, 1, 91-105.

Basu, A., and Frazer, L. N., 1990, Rapid determination of critical  temperature in simulated annealing inversion: Science, 249, 1409-1412.

Bishop, T. N., Bube, K. P., Cutler, R. T., Langan, R. T., Love, P. L., Resnick, J. R., Shuey, R. T., Spindler, D. A., and Wyld, H. W., 1985, Tomographic determination of velocity and depth in laterally varying media: Geophysics, 50, 903-923.

Bohachevsky, I. O., Johnson, M. E., and Stein, M. L., 1986, Generalized simulated annealing for function optimization: Technometrics, 28, 209-217

Bording, R. P., Gersztenkorn, A., Lines, L. R., Scales, J. A., and Treitel, S., 1987, Applications of seismic travel time tomography: Geophysical  Journal of the  Royal Astronomical Society, 90, 285-303.

Cheadle, M. J., Czuchra, B. L., Byrne, T., Ando, C. J., Oliver, J. E., Brown, L. D., Kaufman, S., Malin, P. E., and Phinney, R. A., 1986, The deep crustal structure of the Mojave Desert, California, from COCORP seismic reflection data: Tectonics, 5, 293-320.

Farra, V., and Madariaga, R., 1988, Non-linear reflection tomography: Geophysical Journal, 95, 135-147.

Jackson, D. D., 1972, Interpretation of inaccurate, insufficient, and inconsistent data: Geophysical  Journal of the  Royal Astronomical Society, 28, 97-109.

Kirkpatrick, S., Gelatt, C. D., Jr., and Vecchi, M. P., 1983, Optimization by simulated annealing: Science, 220, 1001-1005.

Kjartansson, E., 1980, Attenuation of seismic waves in rocks: Ph.D. thesis, Stanford University.

Landa, E., Beydoun, W., and Tarantola, A., 1989, Reference velocity model estimation from prestack waveforms: coherency optimization by simulated annealing: Geophysics, 54, 984-990.

Lees, J. M., and Crosson, R. S., 1989, Tomographic inversion for three-dimensional velocity structure at Mount St. Helens using earthquake data: Journal of  Geophysical Research, 94, 5716-5728.

Louie, J. N. and Qin, J., 1991, Structural imaging of the Garlock fault, Cantil Valley, California: Journal of  Geophysical Research, 96, 14,461-14,479.

Metropolis, N., Rosenbluth, A., Rosenbluth, M., Teller, A., and Teller, E., 1953, Equation of state calculations by fast computing machines: J. Chem. Phys., 21, 1087-1092.

Rothman, D. H., 1985, Non-linear inversion, statistical mechanics, and residual statics estimation: Geophysics, 50, 2784-2796.

Rothman, D. H., 1986, Automatic estimation of large residual statics corrections: Geophysics, 51, 332-346.

Sen, M. K., and Stoffa, P. L., 1991, Nonlinear one-dimensional seismic waveform inversion using simulated annealing: Geophysics, 56, 1624-1638.

Shaw, P. R., and Orcutt, J. A., 1985, Waveform inversion of seismic refraction data and applications to young Pacific crust: Geophys. J. Roy. Astr. Soc., 82, 375-414.

Stork, C., and Clayton, R. W., 1986a, Analysis of the resolution between ambiguous velocity and reflector position for travel time tomography: 56th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 545-550.

Stork, C., and Clayton, R. W., 1986b, Case study of travel time tomography applied to data sets containing lateral velocity variations: 56th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 550-553.

Stork, C., and Clayton, R. W., 1987, Applications of tomography to two data sets containing lateral velocity variations: 57th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 839-842.

Tarantola, A., 1987, Inverse Problem Theory, Methods for Data Fitting and Model Parameter Estimation: Elsevier Science Publ. Co., Inc.

Tarantola, A., and Valette, B., 1982, Generalized non-linear inverse problems solved using the least squares criterion: Rev. of Geophys. and Space Phys.,  20, 219-232.

van Trier, J., 1990a, Reflection Tomography after Depth Migration: Field Data Results: 60th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1279-1282.

van Trier, J. A., 1990b, Tomographic determination of structural velocities from depth migrated seismic data: Ph.D. thesis, Stanford Univ.

Vidale, J. E., 1988, Finite-difference calculation of travel times: Bull., Seis. Soc. Am., 78, 2062-2076.

Wiggins, R., 1972, The general inverse problem: implication of surface waves and free oscillations on Earth structure: Rev. of Geophys. and Space Phys., 10, 251-285.

Williamson, P. R., 1990, Tomographic inversion in reflection seismology: Geophys. J. Internat., 100, 255-274.