A combined first-arrival travel time and reflection coherency
optimization approach to velocity estimation

John N. Louie and Sathish K. Pullammanappallil
Seismological Laboratory (174), Mackay School of Mines, University of Nevada, Reno

Geophysical Research Letters, vol. 24, no. 5 (March 1, 1997), p. 511-514.

Available electronically at: http://www.seismo.unr.edu/ftp/pub/louie/papers/grlcoh/grlcoh.html

Abstract

We present a nonlinear optimization method for two-dimensional velocity estimation in strongly heterogeneous media. We achieve this by designing a simulated-annealing algorithm to simultaneously maximize reflection coherency as well as minimize first-arrival travel time residuals. Relatively shallow velocity structure is constrained by first-arrival times, while coherency information from deeper reflections (if present) helps constrain deeper velocities. Pre-stack migration using the optimized velocities images reflections with greater continuity, compared with other velocity estimation methods. Using synthetic and real data we demonstrate that migrations through velocities obtained using our coherency criterion reconstruct basin-boundary structures more accurately than through those obtained using only first-arrival times. Images from shot gathers collected along COCORP Mojave line 5 across the Garlock fault in Cantil Valley, California show new evidence of links between Garlock fault branches and sub-horizontal structures in the upper crust. Imaging reflections and estimating velocities using coherency annealing is more effective and less expensive than picking reflection times from shot gathers recorded complex regions.

Introduction

Travel time inversion suffers from an inherent reflector depth-velocity ambiguity (Stork, 1992a, b). One way to resolve this ambiguity is to make use of additional information contained within the data. Including amplitude information is one such method. Pre-stack migration makes use of reflection amplitudes, but requires prior knowledge of velocities. Conventional velocity analysis methods such as semblance analysis (Toldi, 1989), focusing analysis (Yilmaz and Chambers, 1984), and migration moveout analysis (Al-Yahya, 1989) are done interactively. This could prove to be difficult in areas underlain by complex structures, particularly basin edges and narrow, fault-controlled basins.

Another approach to velocity estimation is by travel-time inversion. The depth of resolution is limited by the type of arrivals used- first arrivals constrain the relatively shallow velocities, to depths usually less than one-third the maximum receiver offset, while reflections provide deeper constraints. Usually travel time inversion with reflected arrivals is performed using times picked off one or a limited number of reflecting horizons, as in Pullammanappallil and Louie (1994). This leads to poor ray coverage in areas of the model not underlain by reflection picks. Moreover, it is often difficult to pick reflection times from noisy or complex land-recorded shot gathers.

To overcome some of these problems several workers have tried to combine pre-stack migration with reflection tomography. One such method casts the migration as a linearized inversion problem, but the need to solve a large linear system makes this method computationally expensive (Ehinger and Lailly, 1993). Jervis et al. (1993) use a genetic algorithm to estimate velocities. Instead of using travel time residual, their objective function is based on a differential semblance operator (Symes and Carazzone, 1991). Their approach requires iterative migration, and to keep the computation time tractable they have to define the velocities as spline functions.

We base our approach on the formalism of Landa et al. (1989). They perform velocity inversion in two steps. First, for a given velocity distribution, they map the depth to reflectors using zero-offset travel time information. Next, the velocity model is updated using the reflection coherency criterion. Their method requires the velocity model to be layered. The interface positions and velocity within each layer also need to be spline functions. Our method differs from this in that we use both travel time and amplitude information simultaneously, and we perform a pre-stack migration only once, using the final velocity model. We need not estimate any velocity distribution or reflector configuration prior to the optimization.

First-arrival times, easily picked, provide information about relatively shallow horizons, while coherency computed from reflection times and amplitudes should constrain velocities above deeper horizons. Performing a pre-stack migration through the velocity model obtained after the optimization images many reflections at the same time. Using simulated annealing for the optimization avoids the need for linearization and the consequent problem of the inversion being trapped in a local minimum. We demonstrate the utility of this technique using synthetic and real data.

Method

Our objective function is made up of two parts:
  1. The least square difference between the picked first-arrival times (tobs) and those calculated through the velocity model (tcal). We denote this as Ett and it is given by the expression,
    Ett = 1/n(sum l=1 to l=n [tobs(l) - tcal(l)]^2) (1)
    where n is the number of observations.

  2. Reflection coherency as defined by Landa et al. (1989). This is given by the equation:
    Scoh = sum i sum k=0 to k=K (sum j Uij[k deltat + tau])^2/(sum j (Uij[k deltat + tau])^2) (2)
    where Uij represents the seismic trace for the ith source and jth receiver, tau is the calculated travel time, deltat is the sample interval, and k deltat is the time window for the semblance calculation.

The objective of our optimization is to simultaneously minimize Ett (equation 1) and maximize Scoh (equation 2). We achieve this with a Monte-Carlo based technique called simulated annealing. The advantages that stem from the nonlinear nature of this optimization process have been documented by several workers (e.g., Sen and Stoffa, 1991; Pullammanappallil and Louie, 1993). By using it in our problem we demonstrate its ability to optimize multiple objective functions simultaneously. This allows us to utilize additional information from the data space, which could have a very different nature, and get more robust solutions in the model space.

Implementation

We extend the single-objective function method developed in Pullammanappallil and Louie (1993) to tackle our problem. We start with a constant velocity model and a flat reflector along the bottom of the model. Each iteration involves perturbing the velocity model and the reflector characteristics. We vary the length, shape, and position of the reflector. As shown by Pullammanappallil and Louie (1993; 1994) this Monte-Carlo process soon removes any dependence on or resemblance to the initial model, exploring the solution space very effectively to locate the global error minimum.

We use a fast finite-difference solution to the eikonal equation (Vidale, 1988) to compute the first-arrival times through the model, while a modification of this (computation of first arrival times from depth points to surface sources and receivers, as in Pullammanappallil and Louie, 1993) is used to get the reflection times. Once we have the reflection times, we compute Scoh using equation (2), which involves adding the amplitudes along the travel time trajectory.

A randomly altered model may be accepted based on four conditions:

  1. If Ett(i) < Ett(i-1) and Scoh(i) > Scoh(i-1), where (i) denotes the iteration number.

  2. If Ett(i) < Ett(i-1) and Scoh(i) < Scoh(i-1), with a probability Pcoh given as,
    Pcoh = exp((Scoh(i)-Smax)^qcoh DeltaScoh/Tcoh) (3)
    where DeltaScoh = Scoh(i-1) - Scoh(i), Tcoh and qcoh are empirical parameters, and Smax is the value of the objective function at the global maximum. We estimate Smax using synthetic tests, by constructing shot gathers through a known velocity model at receiver locations corresponding to the actual survey, and then running our optimization code with this model to estimate its maximum coherency, Smax. Although we can only roughly estimate Smax, its value affects only the speed of our optimization and not its outcome.

  3. If Ett(i) > Ett(i-1) and Scoh(i) > Scoh(i-1), with a probability Ptt given by,
    Ptt = exp((Ett(i)-Emin)^qtt DeltaEtt/Ttt) (4)
    where DeltaEtt = Ett(i-1) - Ett(i), Ttt and qtt are parameters determined empirically, and Emin is the value of the objective function at the global minimum, ideally equal to zero.

  4. If If Ett(i) > Ett(i-1) and Scoh(i) < Scoh(i-1), with a probability P, defined as:
    P = Pcoh x Ptt (5)
    where the above indicates the probability of accepting models with less reflection coherency provided they are first accepted with a probability of Ptt.

We determine the parameters Tcoh, Ttt, qcoh, and qtt with procedures described in detail by Pullammanappallil and Louie (1993; 1994). This involves doing several short runs of the optimization for fixed test values of these parameters, and computing the average least square error for each run. One set of short runs using a range of temperatures gives us the critical temperature Tc, from the run with the least error, which in turn determines how Ttt and Tcoh are varied during the optimization. Additional sets of short runs estimate the values of qtt and qcoh giving the least error.


Figure 1: Test of our optimization methods on recovering a known velocity model section from synthetic data. Louie and Qin (1991) developed the original velocity model, section (c), from their study of COCORP Mojave line 5 reflection data, gravity data, and Garlock fault (GF) tectonic models. For this test we computed acoustic finite-difference synthetic reflection records, maintaining the same source-receiver configuration as the COCORP survey. Section (a) is the result of simulated-annealing optimization of manual first-arrival time picks of the synthetic records (Pullammanappallil and Louie, 1994), which recovers the general features of the low-velocity basin between the southwest and east branches of the Garlock fault (SW GF and east GF, respectively) to about 3 km depth. Section (b) also includes automatic optimization of reflection coherency from the synthetic records. Note the improved match to the known basin shape, and that the optimization constrains velocities to 6 km depths.
(Click on image for printable Adobe Acrobat PDF file.)

Results

We first validate our technique on synthetic data. For compatibility with our data example below, we generate synthetic data from Louie and Qin's (1991) refraction velocity model across the Garlock fault at the location of COCORP Mojave line 5 in Cantil Valley, California (Figure 1c). Thus it models the low-velocity sediments associated with the narrow Cantil Basin, which is bounded by branches of the Garlock fault. We keep the source-receiver configuration identical to the actual survey, with 28 sources mostly towards the north end of the model, each recording into 96 receivers to the south. Maximum source-to-receiver offset is 10 km. The model is 29 km long and 10 km deep, and we use a grid spacing of 250 m in our optimizations. To include the reflection coherency criterion, we generate finite-difference synthetic seismograms through the velocity model. Louie and Qin (1991) show examples of these synthetic gathers. We also hand-pick the gathers for first-arrival times.

Figure 1b shows the synthetic velocity model reconstructed with our coherency criterion. As in Pullammanappallil and Louie (1994), we examine ray distributions traced through the velocity results, and the point-by-point deviations among suites of final models having similar errors, to outline what parts of the sections the optimizations have constrained. Comparing Figure 1b with the synthetic model obtained using only first-arrival time picks (Figure 1a) we see that the depth of the resolved velocities increases from 2 km to about 5 km. The shape of the basin also matches the original (Figure 1c) better when the coherency criterion is used. Next, we image reflections with a pre-stack Kirchhoff migration of the synthetic records through these velocities (Figure 2). The migration through the model shown in Figure 1b images both branches of the Garlock fault and connects them to the basin floor (Figure 2b), while the one through 1a misplaces the elements of the basin boundary, and fails to connect them.

Comparing the migrated image in Figure 2b with the migration through the original velocities (Figure 2c), we see that the main difference is in the sharpness of the basin bounding reflections. This can be attributed to smoothing of the optimized velocities. The lack of sources toward the south side of the model leads to more migration artifacts in that region. This synthetic test serves as a resolution study of the results we will obtain from our optimizations of the real COCORP Mojave line 5 data. The synthetic velocity sections (Figure 1) and the migrations of synthetic data derived from them (Figure 2) point out the strengths and shortcomings of both the velocity optimizations and the pre-stack migrations.


Figure 2: Pre-stack Kirchhoff imaging test of the optimized synthetic velocity models shown in Figure 1. Image (c), from Louie and Qin (1991), migrates finite-difference acoustic synthetic records through the same laterally variable velocities (Figure 1c) used to produce the synthetics. In this best case, the migration recovers the basin geometry almost perfectly, as well as deeper crustal layers on the left side where depth point coverage is good. It also shows that limitations in COCORP line 5 source and receiver geometries prevent proper placement of the crustal layers on the right side where coverage is poor, or good recovery of the steep southwest branch of the Garlock fault. Multiple reflections included in the synthetics also migrate incorrectly to below and outside the basin. Image (a), migrated using the velocities in Figure 1a and without knowledge of the true velocities, reconstructs individual elements of the basin boundaries without placing them correctly. Image (b), which uses the optimized coherency model of Figure 1b, is able to correctly place and link elements of the basin boundary and recovers both crustal layers. The misplaced multiple reflection and the artifacts of poor coverage remain.
(Click on image for printable Adobe Acrobat PDF file.)

To demonstrate the capabilities of our method with a real data set, we picked first-arrival times off 28 shot gathers of 96 traces each, recorded along COCORP Mojave line 5 across the Garlock fault in Cantil Valley, California. The raw shot-gathers are used in the reflection coherency calculation; Louie and Qin (1991) show examples of the raw line 5 data gathers. Figure 3a shows the velocity model using only first-arrival time picks. Only the well-resolved regions of the models are shown, and we get well-resolved velocities only down to 1.5 km. A pre-stack migration through this model is shown in Figure 4a. Lack of good velocity information below 1.5 km leads to poor focusing of coherent energy. Only the Garlock fault east branch is apparent (east GF in Figure 4), and its quality is poor compared to the results of Louie and Qin (1991).

Figure 3c shows the velocity model obtained when using the combined travel time-reflection coherency optimization. Pullammanappallil and Louie (1993) show an example of reflection picks from this data set. We see that the velocity field is smoother and the pre-stack migration through this (Figure 4c) images the east branch of the Garlock fault as well as Louie and Qin (1991), and in the same location as in Pullammanappallil and Louie (1993). The reflections become more coherent and continuous, and we even see hints of deeper reflections.


Figure 3: Tests of our optimizations using data recorded along the northern segment of COCORP Mojave line 5 crossing the Garlock fault in Cantil Valley, California, without knowledge of the true velocities. Only the well-resolved regions are shown (see text). (a) Velocity model optimized using only first-arrival time picks of the data records (Pullammanappallil and Louie, 1994). Velocity constraints become poor below 1.5 km. (b) Velocity model optimized using both hand-picked reflection travel times (Pullammanappallil and Louie, 1993) and picked first-arrival times. Velocities are resolved well only down to 4 km. (c) Velocity estimate optimized using picked first-arrival times together with the reflection coherency criterion. The coherency optimization did not require picking reflection times, unlike (b), and recovers velocity constraints to 5 km depth.
(Click on image for printable Adobe Acrobat PDF file.)

Finally we compare our coherency optimization results (Figures 3c and 4c) with those obtained when using only picked travel times, i.e., picked first arrivals and picked reflections (Figures 3a,b and 4a,b). We do this to see how well we are able to image the data reflections without picking the reflection times. Figure 3b and 4b show the velocities and migrated image obtained using only travel times. Comparing Figures 4b and 4c, we infer that the reflections are imaged better when we include the coherency criterion. Figure 4c is the first image to link the east branch of the Garlock fault (east GF) to the floor of Cantil Basin. A misplaced multiple reflection underlies the east branch. From the Basin floor, a coherent reflection continues 5 km south until it is lost in the artifacts of poor depth point coverage, and its interpretation can only be speculative. Given the data coverage, with receivers mostly to the south, it may not be surprising that we are not able to image the southwest branch of the Garlock fault (SW GF).


Figure 4: Reflection imaging by Kirchhoff migration of COCORP Mojave line 5 data records through the three velocity models shown in Figure 3. (a) Pre-stack migration through the optimized first-arrival model shown in Figure 3a. The data image recovers the reflectivity and dip of the east branch of the Garlock fault (east GF) established by Louie and Qin (1991). While there are hints of imaging of the bottom of Cantil Basin, the artifacts of poor depth point coverage and misplaced multiple reflections dominate the remainder of the image. (b) Migration through the optimized first-arrival and reflection time pick model shown in Figure 3b. The image recovers the Garlock fault east branch reflection (east GF) in essentially the configuration found by Louie and Qin (1991) and by Pullammanappallil and Louie (1993). The hints of additional reflections are still very poorly resolved. (c) Migration through the optimized first-arrival pick and reflection coherency criterion velocity model. Improved focusing and linking between basin-bounding reflections is apparent.
(Click on image for printable Adobe Acrobat PDF file.)

Inferences

Results from synthetic and real data show that including the reflection coherency criterion increases the depth of resolution of the obtained velocities. Consequently, prestack migration through velocities obtained by this method images both the relatively shallow and the deeper reflections. Our comparisons of this approach; i.e., reflection coherency maximization plus first-arrival time minimization, with reflection travel time plus first-arrival time minimization, show that the coherency criterion yields more robust results. With prestack migration, we are able to simultaneously image many reflections. In contrast, a travel-time-only inversion reconstructs just one or a few structures. Our method does not involve significant computations in addition to the travel times. Hence it is computationally comparable to the usual travel time inversion. Our velocity-model-based coherency criterion is able to image even data that contain what appears to be randomly-scattered ``noise.'' The method does not require picking of reflections, which makes it very attractive in areas such as fault zones and continental margins where the complex geology causes the recorded data to have complicated arrivals.

Standard processing of the COCORP Mojave reflection lines by Cheadle et al. (1986) revealed several strong, gently-dipping to sub-horizontal reflections distributed throughout the crust of the Mojave block, south of the Garlock fault. Examining shot records from line 3, which parallels the Garlock about 50 km south of the fault zone, Serpa and Dokka (1992) showed that one of the most prominent reflections was actually a side-swipe reflection from a south-dipping branch of the Garlock fault. On line 5, which crosses the Garlock at high angle as well as intersecting line 3 to the south, we can speculate from our result in Figure 4c that the observations of both Cheadle et al. (1986) and Serpa and Dokka (1992) could be correct. A south-dipping east branch of the Garlock (east GF) appears to flatten at 3.5 km depth and link to a sub-horizontal reflection that proceeds southward below the Mojave block for at least 5 km. Louie and Qin (1991) proposed that the Garlock has a listric geometry and flattens southward into a decollement below the Mojave block. Our results may well show direct evidence of such a structure.

Acknowledgements

This research was supported in part by the National Science Foundation under grants EAR-9296125, EAR-9416224, and EAR-9405534. COCORP Mojave line 5 data were kindly provided by Sid Kaufman of Cornell University. The manuscript was greatly improved with modifications suggested by Mrinal K. Sen and two anonymous reviewers.

References

Al-Yahya, K., Velocity analysis by iterative profile migration, Geophysics, 54, 718-729, 1989.

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., The deep crustal structure of the Mojave Desert, California, from COCORP seismic reflection data: Tectonics, 5, 293-320, 1986.

Ehinger, A., and P. Lailly, Prestack imaging by coupled linearized inversion, in Mathematical Methods in Geophysical Imaging, Proc. of the Society of Photo-Optical Instrumentation, San Diego, California, 2033, 66-87, 1993.

Jervis, M., P. L. Stoffa, and M. A. Sen, 2-D migration velocity estimation using a genetic algorithm, Geophys. Res. Lett., 20, 1495-1498, 1993.

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

Louie, J. N. and J. Qin, Structural imaging of the Garlock fault, Cantil Valley, California, J. Geophys. Res., 96, 14,461-14,479, 1991.

Pullammanappallil, S. K. and J. N. Louie, Inversion of seismic reflection traveltimes using a nonlinear optimization scheme, Geophysics, 58, 1607-1620, 1993.

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

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

Serpa, L., and Dokka, R. K., Geometry of the Garlock fault zone based on seismic reflection data: J. Geophys. Res., 97, 15,297-15,306, 1992.

Stork, C., Singular value decomposition of the velocity-reflector depth tradeoff, part 1: Introduction using a two-parameter model, Geophysics, 57, 927-932, 1992a.

Stork, C., Singular value decomposition of the velocity-reflector depth tradeoff, part 2: High-resolution analysis of a generic model, Geophysics, 57, 933-943, 1992b.

Symes, W. W., and J. J. Carazzone, Velocity inversion by differential semblance optimization, Geophysics, 56, 654-663, 1991.

Toldi, J. L., Velocity analysis without picking, Geophysics, 34, 191-199, 1989.

Vidale, J. E., Finite difference calculation of travel times, Bull. Seismol. Soc. Amer., 78, 2062-2076, 1988.

Yilmaz, O., and R. Chambers, Migration velocity analysis by wavefield extrapolation, Geophysics, 49, 1664-1674, 1984.